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

4 ビュー (過去 30 日間)
Marcel
Marcel 2025 年 7 月 3 日
回答済み: Star Strider 2025 年 7 月 3 日
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

サインインしてこの質問に回答する。

回答 (1 件)

Star Strider
Star Strider 2025 年 7 月 3 日
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
.

カテゴリ

Help Center および File Exchange で Pulse and Transition Metrics についてさらに検索

タグ

製品


リリース

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by