Error in step size and i dont know why. I need help..... please.

4 views (last 30 days)
Marcel
Marcel on 3 Jul 2025
Answered: Star Strider on 3 Jul 2025
Ran in:
clc;
clear;
close all;
% Startwerte
y0 = [0; 0; 0; 0; 0; 0;30e5; 30e5; 30e5; 288.15]; % Anfangswerte für x1,x2,z,v1,v2,vz,p2,p1,p3,T
% Zeitintervall
tspan = linspace(0,10,500);
% ODE lösen
% ODE-Optionen setzen
options = odeset('RelTol', 1e-5, 'AbsTol', 1e-8)
options = struct with fields:
AbsTol: 1.0000e-08 BDF: [] Events: [] InitialStep: [] Jacobian: [] JConstant: [] JPattern: [] Mass: [] MassSingular: [] MaxOrder: [] MaxStep: [] MinStep: [] NonNegative: [] NormControl: [] OutputFcn: [] OutputSel: [] Refine: [] RelTol: 1.0000e-05 Stats: [] Vectorized: [] MStateDependence: [] MvPattern: [] InitialSlope: []
[t, y] = ode113(@damper, tspan, y0,options);
Warning: Failure at t=6.505343e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t.
function dydt = damper(t,y)
% Variablen
x1 = y(1);
x2 = y(2);
z = y(3);
v1 = y(4);
v2 = y(5);
vz = y(6);
p2 = y(7);
p1 = y(8);
p3 = y(9);
T = y(10);
% Anregung des Hauptkolbens
f = 0.1; % Frequenz [Hz]
A = 0.01; % Amplitude [m]
omega = 2*pi*f;
%x= A * cos(omega*t);
%v = -A * omega * sin(omega*t);
%a = -A * omega^2 * cos(omega*t);
x = A * sin(omega*t);
v = A * omega * cos(omega*t);
a = -A * omega^2 * sin(omega*t);
%A soll bei 0 Starten und entsprechend der Zeitschritte maximal in der Frequenz die Wete annhemen
% Konstanten
FR = 10;
L1 = 30e-3;
L2 = 50e-3;
L3 = 54e-3;
kappa = 1.4;
M1= 0.256;
M2= 0.06;
c1= 60;
c2= 60;
k1= 1600;
k2= 200;
Fi1= 20;
Fi2= 20;
m1 = 0.008;
m2 = 0.006;
% Geometrie und Parameter berechnen
A1 = KolbenStange();
A2 = Kolben();
Atk = FlaechTrennkolben();
Ap1 = Durchlasskompression1();
Ap2 = Durchlasskompression2();
% Flächen Ventile (jetzt mit aktuellem x1 und x2 als Werte)
AV1 = FlaecheVentil1(x1);
AV2 = FlaecheVentil2(x2);
% Druckdifferenzen
dp_12 = dp12(p1,p2);
dp_23 = dp23(p2,p3);
% Temperaturabhängige Größen
mu = Viskositaet(T);
RHO = Dichte(T);
betta = Kompressionsmodul();
% Strömungs- und Widerstandskoeffizienten
REcv_val = REcv(T,p1,p2);
CDV_val = CDV(p1,p2,T);
% Volumina mit aktuellen x,z (für V2 auch z)
V1 = Volumen1(x);
V2 = Volumen2(x,z);
V3 = Volumen3(z);
if V1 <= 0 || V2 <= 0 || V3 <= 0
error('Negatives oder Nullvolumen bei t=%f: V1=%.3e, V2=%.3e, V3=%.3e', t, V1, V2, V3);
end
% Volumenströme (hier brauchst du auch aktuelle x1, x2)
QV1_val = QV1(p1,p2,T,x2); % z.B. x2 als Ventilposition
QV2_val = QV2(p1,p2,T,x1);
QL_val = Ql(p1,p2,T,v1);
% Ventildrücke
pv1 = Ventildruck1(p1,p2,T,x1);
pv2 = Ventildruck2(p1,p2,T,x2);
% Dämpferkraft
FD = damperkraft(p1,p2,a);
% Kraft auf Shim
Fp1_val= Fp1(p1,p2,T,x2);
Fp2_val= Fp2(p1,p2,T,x1);
% Moment auf Fm
Fm1_val= Fm1(p1,p2,T,x2);
Fm2_val= Fm2(p1,p2,T,x1);
% Differentialgleichungen
dydt = zeros(10,1);
dydt(1) = y(4); % dx1/dt = v1
dydt(2) = y(5); % dx2/dt = v2
dydt(3) = y(6); % dz/dt = vz
% Beispielhafte Kraftgleichungen (platzhalter - anpassen!)
dydt(4) = (-c1*y(5)-k1*y(1)+Fp1_val+Fm1_val+Fi1+m1*a)/m1;
dydt(5) = (-c2*y(6)-k2*y(2)+Fp2_val+Fm2_val+Fi2+m2*a)/m2;
dydt(6) = (Atk * (y(7) - y(9)) - FR); % vz Ableitung
% Druckänderungen
dydt(7) = ( - QV1_val - QL_val + A2*v - Atk*y(6)/(betta*(A2*(L2 - x) + Atk*y(3))));
dydt(8) = (QV2_val + QL_val - A1*v/(betta*(A1*(L1 + x))));
dydt(9) = -kappa * (y(9) * (-Atk*y(6)) / (Atk*(L3 - y(3))));
dydt(10) = 0; % Temperaturänderung (Platzhalter)
end
function Atk = FlaechTrennkolben()
D = 35.7e-3; % Durchmesser Trennkolben
Atk = (pi/4)*(D^2);
end
function A1 = KolbenStange()
D1 = 35.9e-3;
d1 = 11e-3;
A1 = (pi/4)*(D1^2) - (pi/4)*(d1^2);
end
function A2 = Kolben()
D1 = 35.9e-3;
A2 = (pi/4)*(D1^2) ;
end
function Ap1 = Durchlasskompression1()
%ist aktiv bei pvq
dp1= 4.95e-3 ;
Ap1 = (pi/4)*(dp1^2)*4;
end
function Ap2 = Durchlasskompression2 ()
% ist Aktiv bei pv2
dp2 = 3.3e-3 ;
Ap2 = (pi/4)*(dp2^2)*8;
end
function AV1 = FlaecheVentil1(x1)
alpha = 0.4; %?3
Dmaxs = 40e-3;
AV1 = alpha * pi * Dmaxs * x1;
end
function AV2 = FlaecheVentil2(x2)
alpha = 0.4;
Dmaxs = 40e-3;
AV2 = alpha * pi * Dmaxs * x2;
end
%% Temperaturabhängige Funktionen
function RHO = Dichte(T)
alphat = 1e-3;
rho0 = 840;
T0 = 288.15;
RHO = rho0 / (1 + alphat*(T - T0));
end
function mu = Viskositaet(T)
C = 3717;
mu0 = 0.01;
T0 = 288.15;
mu = mu0 * exp(C*(1/T - 1/T0));
end
function betta = Kompressionsmodul()
betta0 = 6.6e-10;
r = 17.985e-3;
e = 200e9;
t = 1.6e-3;
betta = betta0 + 2*r/(e*t);
end
function QV1_val = QV1(p1,p2,T,x2)
AV2 = FlaecheVentil2(x2);
CDV1 = 0.5;
pv1 = Ventildruck1(p1,p2,T,x2);
RHO = Dichte(T);
QV1_val = AV2 * CDV1 * sign(p2 - pv1) * sqrt(2 * abs(p2 - pv1) / RHO);
end
function QV2_val = QV2(p1,p2,T,x1)
AV1 = FlaecheVentil1(x1);
CDV2 = 0.5;
pv2 = Ventildruck2(p1,p2,T,x1);
RHO = Dichte(T);
QV2_val = AV1 * CDV2 * sign(pv2 - p1) * sqrt(2 * abs(pv2 - p1) / RHO);
end
function q = Ql(p1,p2,T,v)
b = 0.00035;
s = 10e-2;
D1 = 35.9e-3;
eta = 0.01 * exp(3717 * (1./T - 1/288.15));
q = ((abs(p2 - p1) * b^3) / (12 * s * eta) + v * b / 2) * pi * D1;
end
%Kräfte
function Fp1_val = Fp1(p1,p2,T,x2)
Ap1 = Durchlasskompression1();
pv2=Ventildruck2(p1,p2,T,x2);
Fp1_val=Ap1*(pv2-p1);
end
function Fp2_val = Fp2(p1,p2,T,x1)
Ap2= Durchlasskompression2();
pv1 = Ventildruck1(p1,p2,T,x1);
Fp2_val= Ap2*(pv1-p2);
end
function Fm1_val = Fm1(p1,p2,T,x2)
cf= 0.3;
QV2_val = QV2(p1,p2,T,x2);
Ap1= Durchlasskompression1();
RHO= Dichte(T);
pv2 = Ventildruck2(p1,p2,T,x2);
Fm1_val= ((cf*RHO*QV2_val^2)*sign(pv2-p1))/Ap1;
end
function Fm2_val = Fm2(p1,p2,T,x1)
cf= 0.3;
QV1_val = QV1(p1,p2,T,x1);
Ap2= Durchlasskompression2();
RHO= Dichte(T);
pv1= Ventildruck1(p1,p2,T,x1);
Fm2_val= ((cf*RHO*QV1_val^2)*sign(pv1-p2))/Ap2;
end
%% Reynoldszahl und Cd-Funktionen
function REcv_val = REcv(T,p1,p2)
mu = Viskositaet(T);
Lv = 2e-2;
RHO = Dichte(T);
dp = abs(p2-p1);
vstr = sqrt(2 * dp / RHO);
REcv_val = RHO * vstr * Lv / mu;
end
function CDV_val = CDV(p1,p2,T)
dv = 10e-3;
mu = Viskositaet(T);
RHO = Dichte(T);
nu = mu / RHO;
REcv_val = REcv(T,p1,p2);
dp_12 = abs(p2 - p1);
CDVmax = 0.7;
CDV_val = CDVmax * tanh((2*dv/nu*REcv_val)*sqrt(2*(dp_12/RHO)));
end
%% Volumenfunktionen
function V1 = Volumen1(x)
L1 = 30e-3;
A1 = Kolben();
V1 = A1 * (L1 + x);
end
function V2 = Volumen2(x,z)
A2 = KolbenStange();
Atk = FlaechTrennkolben();
L2 = 50e-3;
V2 = A2 * (L2 - x) + Atk * z;
end
function V3 = Volumen3(z)
Atk = FlaechTrennkolben();
L3 = 54e-3;
V3 = Atk * (L3 - (z));
end
%% Dämpferkraft
function FD = damperkraft(p1,p2,a)
M1 = 0.256;
FR = 10;
A1 = KolbenStange() ;
A2 = Kolben();
FD = sign(p1-p2)*FR - p1 * A1 + p2 * A2 + M1 * a ;
end
%% Ventildruckfunktionen
function pv1 = Ventildruck1(p1,p2,T,x1)
AV1 = FlaecheVentil1(x1);
Ap1 = Durchlasskompression1();
CDC1 = 0.6;
CDV1 = CDV(p1,p2,T);
pv1 = (CDC1^2 * Ap1^2 * p1 + CDV1^2 * AV1^2 * p2) / (CDC1^2 * Ap1^2 + CDV1^2 * AV1^2);
end
function pv2 = Ventildruck2(p1,p2,T,x2)
AV2 = FlaecheVentil2(x2);
Ap2 = Durchlasskompression2;
CDC2 = 0.4;
CDV2 = CDV(p1,p2,T);
pv2 = (CDC2^2 * Ap2^2 * p1 + CDV2^2 * AV2^2 * p2) / (CDC2^2 * Ap2^2 + CDV2^2 * AV2^2);
end
%% Druckdifferenzen
function dp = dp12(p1,p2)
dp = abs(p2 - p1);
end
function dp = dp23(p2,p3)
dp = abs(p3 - p2);
end
% Plot p1, p2 und p3 über die Zeit
figure;
% p1 aus y(:,8)
subplot(3,1,1);
plot(t, y(:,8), '-o', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('p1 [Pa]');
title('Druck p1 über die Zeit');
grid on;
% p2 aus y(:,7)
subplot(3,1,2);
plot(t, y(:,7), '-o', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('p2 [Pa]');
title('Druck p2 über die Zeit');
grid on;
% p3 aus y(:,9)
subplot(3,1,3);
plot(t, y(:,9), '-o', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('p3 [Pa]');
title('Druck p3 über die Zeit');
grid on;
% Volumen V1, V2, V3 aus den Ergebnissen berechnen
V1 = zeros(length(t),1);
V2 = zeros(length(t),1);
V3 = zeros(length(t),1);
for i = 1:length(t)
x = y(i,1); % x1 in deinem System ist y(:,1)
z = y(i,3);
V1(i) = Volumen1(x);
V2(i) = Volumen2(x,z);
V3(i) = Volumen3(z);
end
% Plot V1, V2 und V3 in einem Diagramm
figure;
plot(t, V1, '-o', 'LineWidth', 1.5);
hold on;
plot(t, V2, '-x', 'LineWidth', 1.5);
plot(t, V3, '-s', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('Volumen [m^3]');
title('Volumen V1, V2 und V3 über die Zeit');
legend('V1','V2','V3');
grid on;
% Volumenströme QV1, QV2 und QL aus den Ergebnissen berechnen
QV1_val = zeros(length(t),1);
QV2_val = zeros(length(t),1);
QL_val = zeros(length(t),1);
for i = 1:length(t)
p1 = y(i,8);
p2 = y(i,7);
T = y(i,10);
x1 = y(i,1);
x2 = y(i,2);
v1 = y(i,4);
QV1_val(i) = QV1(p1,p2,T,x2);
QV2_val(i) = QV2(p1,p2,T,x1);
QL_val(i) = Ql(p1,p2,T,v1);
end
% Plot der Volumenströme in einem Diagramm
figure;
plot(t, QV1_val, '-o', 'LineWidth', 1.5);
hold on;
plot(t, QV2_val, '-x', 'LineWidth', 1.5);
plot(t, QL_val, '-s', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('Volumenstrom [m^3/s]');
title('Volumenströme QV1, QV2 und QL über die Zeit');
legend('QV1','QV2','QL');
grid on

Sign in to answer this question.

Answers (1)

Star Strider
Star Strider on 3 Jul 2025
Ran in:
The warning simply means that ode113 has encountered a singluarity in one of your differential equations. It is not obvious to me what is causing that.
I added a save call for 'dydt' and printed out the derivatives just before the failure. I suggest that you add any other variables to the save call and then add appropriate variables after the load call to display them.
clc;
clear;
close all;
% Startwerte
y0 = [0; 0; 0; 0; 0; 0;30e5; 30e5; 30e5; 288.15]; % Anfangswerte für x1,x2,z,v1,v2,vz,p2,p1,p3,T
% Zeitintervall
tspan = linspace(0,10,500);
% ODE lösen
% ODE-Optionen setzen
options = odeset('RelTol', 1e-5, 'AbsTol', 1e-8)
options = struct with fields:
AbsTol: 1.0000e-08 BDF: [] Events: [] InitialStep: [] Jacobian: [] JConstant: [] JPattern: [] Mass: [] MassSingular: [] MaxOrder: [] MaxStep: [] MinStep: [] NonNegative: [] NormControl: [] OutputFcn: [] OutputSel: [] Refine: [] RelTol: 1.0000e-05 Stats: [] Vectorized: [] MStateDependence: [] MvPattern: [] InitialSlope: []
[t, y] = ode113(@damper, tspan, y0,options);
Warning: Failure at t=6.505343e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t.
load('Check_DEs.mat')
format longE
disp('Derivatives: ')
Derivatives:
dydt
dydt = 10×1
1.0e+00 * -1.491541833379143e+10 1.564429053897165e-01 5.148700178849773e-04 -1.442389690172004e+13 -2.789991148555638e+05 -9.851629603492505e+00 -1.402757225602631e+07 1.210199568451514e+09 4.004525350597411e+04 0
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
function dydt = damper(t,y)
% Variablen
x1 = y(1);
x2 = y(2);
z = y(3);
v1 = y(4);
v2 = y(5);
vz = y(6);
p2 = y(7);
p1 = y(8);
p3 = y(9);
T = y(10);
% Anregung des Hauptkolbens
f = 0.1; % Frequenz [Hz]
A = 0.01; % Amplitude [m]
omega = 2*pi*f;
%x= A * cos(omega*t);
%v = -A * omega * sin(omega*t);
%a = -A * omega^2 * cos(omega*t);
x = A * sin(omega*t);
v = A * omega * cos(omega*t);
a = -A * omega^2 * sin(omega*t);
%A soll bei 0 Starten und entsprechend der Zeitschritte maximal in der Frequenz die Wete annhemen
% Konstanten
FR = 10;
L1 = 30e-3;
L2 = 50e-3;
L3 = 54e-3;
kappa = 1.4;
M1= 0.256;
M2= 0.06;
c1= 60;
c2= 60;
k1= 1600;
k2= 200;
Fi1= 20;
Fi2= 20;
m1 = 0.008;
m2 = 0.006;
% Geometrie und Parameter berechnen
A1 = KolbenStange();
A2 = Kolben();
Atk = FlaechTrennkolben();
Ap1 = Durchlasskompression1();
Ap2 = Durchlasskompression2();
% Flächen Ventile (jetzt mit aktuellem x1 und x2 als Werte)
AV1 = FlaecheVentil1(x1);
AV2 = FlaecheVentil2(x2);
% Druckdifferenzen
dp_12 = dp12(p1,p2);
dp_23 = dp23(p2,p3);
% Temperaturabhängige Größen
mu = Viskositaet(T);
RHO = Dichte(T);
betta = Kompressionsmodul();
% Strömungs- und Widerstandskoeffizienten
REcv_val = REcv(T,p1,p2);
CDV_val = CDV(p1,p2,T);
% Volumina mit aktuellen x,z (für V2 auch z)
V1 = Volumen1(x);
V2 = Volumen2(x,z);
V3 = Volumen3(z);
if V1 <= 0 || V2 <= 0 || V3 <= 0
error('Negatives oder Nullvolumen bei t=%f: V1=%.3e, V2=%.3e, V3=%.3e', t, V1, V2, V3);
end
% Volumenströme (hier brauchst du auch aktuelle x1, x2)
QV1_val = QV1(p1,p2,T,x2); % z.B. x2 als Ventilposition
QV2_val = QV2(p1,p2,T,x1);
QL_val = Ql(p1,p2,T,v1);
% Ventildrücke
pv1 = Ventildruck1(p1,p2,T,x1);
pv2 = Ventildruck2(p1,p2,T,x2);
% Dämpferkraft
FD = damperkraft(p1,p2,a);
% Kraft auf Shim
Fp1_val= Fp1(p1,p2,T,x2);
Fp2_val= Fp2(p1,p2,T,x1);
% Moment auf Fm
Fm1_val= Fm1(p1,p2,T,x2);
Fm2_val= Fm2(p1,p2,T,x1);
% Differentialgleichungen
dydt = zeros(10,1);
dydt(1) = y(4); % dx1/dt = v1
dydt(2) = y(5); % dx2/dt = v2
dydt(3) = y(6); % dz/dt = vz
% Beispielhafte Kraftgleichungen (platzhalter - anpassen!)
dydt(4) = (-c1*y(5)-k1*y(1)+Fp1_val+Fm1_val+Fi1+m1*a)/m1;
dydt(5) = (-c2*y(6)-k2*y(2)+Fp2_val+Fm2_val+Fi2+m2*a)/m2;
dydt(6) = (Atk * (y(7) - y(9)) - FR); % vz Ableitung
% Druckänderungen
dydt(7) = ( - QV1_val - QL_val + A2*v - Atk*y(6)/(betta*(A2*(L2 - x) + Atk*y(3))));
dydt(8) = (QV2_val + QL_val - A1*v/(betta*(A1*(L1 + x))));
dydt(9) = -kappa * (y(9) * (-Atk*y(6)) / (Atk*(L3 - y(3))));
dydt(10) = 0; % Temperaturänderung (Platzhalter)
if t >= 6.505342e-01
save('Check_DEs.mat','dydt')
end
end
function Atk = FlaechTrennkolben()
D = 35.7e-3; % Durchmesser Trennkolben
Atk = (pi/4)*(D^2);
end
function A1 = KolbenStange()
D1 = 35.9e-3;
d1 = 11e-3;
A1 = (pi/4)*(D1^2) - (pi/4)*(d1^2);
end
function A2 = Kolben()
D1 = 35.9e-3;
A2 = (pi/4)*(D1^2) ;
end
function Ap1 = Durchlasskompression1()
%ist aktiv bei pvq
dp1= 4.95e-3 ;
Ap1 = (pi/4)*(dp1^2)*4;
end
function Ap2 = Durchlasskompression2 ()
% ist Aktiv bei pv2
dp2 = 3.3e-3 ;
Ap2 = (pi/4)*(dp2^2)*8;
end
function AV1 = FlaecheVentil1(x1)
alpha = 0.4; %?3
Dmaxs = 40e-3;
AV1 = alpha * pi * Dmaxs * x1;
end
function AV2 = FlaecheVentil2(x2)
alpha = 0.4;
Dmaxs = 40e-3;
AV2 = alpha * pi * Dmaxs * x2;
end
%% Temperaturabhängige Funktionen
function RHO = Dichte(T)
alphat = 1e-3;
rho0 = 840;
T0 = 288.15;
RHO = rho0 / (1 + alphat*(T - T0));
end
function mu = Viskositaet(T)
C = 3717;
mu0 = 0.01;
T0 = 288.15;
mu = mu0 * exp(C*(1/T - 1/T0));
end
function betta = Kompressionsmodul()
betta0 = 6.6e-10;
r = 17.985e-3;
e = 200e9;
t = 1.6e-3;
betta = betta0 + 2*r/(e*t);
end
function QV1_val = QV1(p1,p2,T,x2)
AV2 = FlaecheVentil2(x2);
CDV1 = 0.5;
pv1 = Ventildruck1(p1,p2,T,x2);
RHO = Dichte(T);
QV1_val = AV2 * CDV1 * sign(p2 - pv1) * sqrt(2 * abs(p2 - pv1) / RHO);
end
function QV2_val = QV2(p1,p2,T,x1)
AV1 = FlaecheVentil1(x1);
CDV2 = 0.5;
pv2 = Ventildruck2(p1,p2,T,x1);
RHO = Dichte(T);
QV2_val = AV1 * CDV2 * sign(pv2 - p1) * sqrt(2 * abs(pv2 - p1) / RHO);
end
function q = Ql(p1,p2,T,v)
b = 0.00035;
s = 10e-2;
D1 = 35.9e-3;
eta = 0.01 * exp(3717 * (1./T - 1/288.15));
q = ((abs(p2 - p1) * b^3) / (12 * s * eta) + v * b / 2) * pi * D1;
end
%Kräfte
function Fp1_val = Fp1(p1,p2,T,x2)
Ap1 = Durchlasskompression1();
pv2=Ventildruck2(p1,p2,T,x2);
Fp1_val=Ap1*(pv2-p1);
end
function Fp2_val = Fp2(p1,p2,T,x1)
Ap2= Durchlasskompression2();
pv1 = Ventildruck1(p1,p2,T,x1);
Fp2_val= Ap2*(pv1-p2);
end
function Fm1_val = Fm1(p1,p2,T,x2)
cf= 0.3;
QV2_val = QV2(p1,p2,T,x2);
Ap1= Durchlasskompression1();
RHO= Dichte(T);
pv2 = Ventildruck2(p1,p2,T,x2);
Fm1_val= ((cf*RHO*QV2_val^2)*sign(pv2-p1))/Ap1;
end
function Fm2_val = Fm2(p1,p2,T,x1)
cf= 0.3;
QV1_val = QV1(p1,p2,T,x1);
Ap2= Durchlasskompression2();
RHO= Dichte(T);
pv1= Ventildruck1(p1,p2,T,x1);
Fm2_val= ((cf*RHO*QV1_val^2)*sign(pv1-p2))/Ap2;
end
%% Reynoldszahl und Cd-Funktionen
function REcv_val = REcv(T,p1,p2)
mu = Viskositaet(T);
Lv = 2e-2;
RHO = Dichte(T);
dp = abs(p2-p1);
vstr = sqrt(2 * dp / RHO);
REcv_val = RHO * vstr * Lv / mu;
end
function CDV_val = CDV(p1,p2,T)
dv = 10e-3;
mu = Viskositaet(T);
RHO = Dichte(T);
nu = mu / RHO;
REcv_val = REcv(T,p1,p2);
dp_12 = abs(p2 - p1);
CDVmax = 0.7;
CDV_val = CDVmax * tanh((2*dv/nu*REcv_val)*sqrt(2*(dp_12/RHO)));
end
%% Volumenfunktionen
function V1 = Volumen1(x)
L1 = 30e-3;
A1 = Kolben();
V1 = A1 * (L1 + x);
end
function V2 = Volumen2(x,z)
A2 = KolbenStange();
Atk = FlaechTrennkolben();
L2 = 50e-3;
V2 = A2 * (L2 - x) + Atk * z;
end
function V3 = Volumen3(z)
Atk = FlaechTrennkolben();
L3 = 54e-3;
V3 = Atk * (L3 - (z));
end
%% Dämpferkraft
function FD = damperkraft(p1,p2,a)
M1 = 0.256;
FR = 10;
A1 = KolbenStange() ;
A2 = Kolben();
FD = sign(p1-p2)*FR - p1 * A1 + p2 * A2 + M1 * a ;
end
%% Ventildruckfunktionen
function pv1 = Ventildruck1(p1,p2,T,x1)
AV1 = FlaecheVentil1(x1);
Ap1 = Durchlasskompression1();
CDC1 = 0.6;
CDV1 = CDV(p1,p2,T);
pv1 = (CDC1^2 * Ap1^2 * p1 + CDV1^2 * AV1^2 * p2) / (CDC1^2 * Ap1^2 + CDV1^2 * AV1^2);
end
function pv2 = Ventildruck2(p1,p2,T,x2)
AV2 = FlaecheVentil2(x2);
Ap2 = Durchlasskompression2;
CDC2 = 0.4;
CDV2 = CDV(p1,p2,T);
pv2 = (CDC2^2 * Ap2^2 * p1 + CDV2^2 * AV2^2 * p2) / (CDC2^2 * Ap2^2 + CDV2^2 * AV2^2);
end
%% Druckdifferenzen
function dp = dp12(p1,p2)
dp = abs(p2 - p1);
end
function dp = dp23(p2,p3)
dp = abs(p3 - p2);
end
% Plot p1, p2 und p3 über die Zeit
figure;
% p1 aus y(:,8)
subplot(3,1,1);
plot(t, y(:,8), '-o', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('p1 [Pa]');
title('Druck p1 über die Zeit');
grid on;
% p2 aus y(:,7)
subplot(3,1,2);
plot(t, y(:,7), '-o', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('p2 [Pa]');
title('Druck p2 über die Zeit');
grid on;
% p3 aus y(:,9)
subplot(3,1,3);
plot(t, y(:,9), '-o', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('p3 [Pa]');
title('Druck p3 über die Zeit');
grid on;
●
% Volumen V1, V2, V3 aus den Ergebnissen berechnen
V1 = zeros(length(t),1);
V2 = zeros(length(t),1);
V3 = zeros(length(t),1);
for i = 1:length(t)
x = y(i,1); % x1 in deinem System ist y(:,1)
z = y(i,3);
V1(i) = Volumen1(x);
V2(i) = Volumen2(x,z);
V3(i) = Volumen3(z);
end
% Plot V1, V2 und V3 in einem Diagramm
figure;
plot(t, V1, '-o', 'LineWidth', 1.5);
hold on;
plot(t, V2, '-x', 'LineWidth', 1.5);
plot(t, V3, '-s', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('Volumen [m^3]');
title('Volumen V1, V2 und V3 über die Zeit');
legend('V1','V2','V3');
grid on;
% Volumenströme QV1, QV2 und QL aus den Ergebnissen berechnen
QV1_val = zeros(length(t),1);
QV2_val = zeros(length(t),1);
QL_val = zeros(length(t),1);
for i = 1:length(t)
p1 = y(i,8);
p2 = y(i,7);
T = y(i,10);
x1 = y(i,1);
x2 = y(i,2);
v1 = y(i,4);
QV1_val(i) = QV1(p1,p2,T,x2);
QV2_val(i) = QV2(p1,p2,T,x1);
QL_val(i) = Ql(p1,p2,T,v1);
end
% Plot der Volumenströme in einem Diagramm
figure;
plot(t, QV1_val, '-o', 'LineWidth', 1.5);
hold on;
plot(t, QV2_val, '-x', 'LineWidth', 1.5);
plot(t, QL_val, '-s', 'LineWidth', 1.5);
xlabel('Zeit [s]');
ylabel('Volumenstrom [m^3/s]');
title('Volumenströme QV1, QV2 und QL über die Zeit');
legend('QV1','QV2','QL');
grid on
.

Categories

Find more on Pulse and Transition Metrics in Help Center and File Exchange

Tags

Products


Release

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!