ODE45 outputs only zeros
5 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Nader Mohamed
il 31 Ott 2021
Risposto: Bjorn Gustavsson
il 1 Nov 2021
I'm trying to use ode45 to model the movement of a piston. This works though a command (function command) that opens with respect of time. When I try to run the code the only output I get is zeros even if the output I expect for yy(:,1) should be something like the image attached.
I don't understand what's wrong with my code and why it outputs zeros instead of numbers.
t0 = 0;
t1 = 1;
t2 = 1.5;
tf = 3;
options = odeset('RelTol',1e-9,'AbsTol',1e-9);
[tt,yy] = ode45(@problem,[0,3],[0,0,1e-3,0]);
plot(tt,yy(:,1))
function z = command(t,t1,t2) %FUNCTION FOR COMMAND Z
z=zeros(size(t));
z(t>t1&t<t2)=(t(t>t1&t<t2)-t1)*t1/(t2-t1);
z(t>=t2)=t1;
end
function dydt = problem(t,y)
%------------------------------------------------------------------------------------%
data.fluid.rho = 890;
data.accumulator.V_N2 = 10e-3;
data.accumulator.P_N2 = 2.5e6;
data.accumulator.p0 = 21e6;
data.accumulator.gamma = 1.2;
data.accumulator.V0 = (data.accumulator.P_N2*data.accumulator.V_N2)/(data.accumulator.p0);
data.delivery.ka = 1.12;
data.delivery.kcv = 2;
data.delivery.D23 = 18e-3;
data.delivery.L23 = 2;
data.delivery.f23 = 0.032;
data.distributor.kd = 12;
data.distributor.d0 = 5e-3;
data.actuator.Dc = 50e-3;
data.actuator.Dr = 22e-3;
data.actuator.ma = 2;
data.actuator.xmax = 200e-3;
data.actuator.F0 = 1e3;
data.actuator.k = 120e3;
data.return.Dr = 18e-3;
data.return.L67 = 15;
data.return.f67 = 0.035;
data.tank.pT = 0.1e6;
data.tank.VT = 1e-3;
data.tank.kt = 1.12;
%------------------------------------------------------------------------------------%
t1 = 1;
t2 = 1.5;
z = command(t,t1,t2);
alpha = 2*acos(1 - 2*abs(z));
A23 = (data.delivery.D23^2/4)*pi;
Ar = (data.return.Dr^2/4)*pi;
Ad = ((data.distributor.d0^2/4)*(alpha - sin(alpha)))+sqrt(eps);
A4 = (data.actuator.F0 + data.actuator.k*y(1))/ (data.accumulator.p0 - (data.tank.pT*(1-0.3)));
A5 = (1-0.3)*A4;
Peq = (data.tank.pT*A5 + (data.actuator.F0 + data.actuator.k*y(1)))/A4;
PA = data.accumulator.p0 * (data.accumulator.V0/(data.accumulator.V_N2 - y(4)))^data.accumulator.gamma;
P1 = PA - 0.5*data.delivery.ka*data.fluid.rho*((y(2)*A4)/A23)^2;
P2 = P1 - 0.5*data.delivery.ka*data.fluid.rho*((y(2)*A4)/A23)^2;
P3 = P2 - data.delivery.f23*(data.delivery.L23/data.delivery.D23)*0.5*data.fluid.rho*((y(2)*A4)/A23)^2;
P7 = data.tank.pT + 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A5)/Ar)^2;
P6 = P7 + data.return.f67*(data.return.L67/data.return.Dr)*0.5*data.fluid.rho*((y(2)*A5)/Ar)^2;
% P4 = PA - (data.fluid.rho*(y(2)*A4)^2)*0.5*( (data.delivery.ka/A23^2) + (data.delivery.kcv/A23^2) + (data.distributor.kd/Ad^2) + (data.delivery.f23*(data.delivery.L23/(data.delivery.D23*A23^2))));
% P5 = data.tank.pT + (data.fluid.rho*(y(2)*A5)^2)*0.5*( (data.tank.kt/Ar^2) + (data.distributor.kd/Ad^2) + (data.return.f67*(data.return.L67/(data.return.Dr*Ar^2))));
if z == 0
P4 = Peq;
P5 = data.tank.pT;
elseif z > 0
if Ad < 1e-9
P4 = Peq;
P5 = data.tank.pT;
else
P4 = P3 - 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A4)/Ad)^2;
P5 = P6 + 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A5)/Ad)^2;
end
else
if Ad < 1e-9
P4 = Peq;
P5 = data.tank.pT;
else
P4 = P6 + 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A5)/Ad)^2;
P5 = P3 - 0.5*data.distributor.kd*data.fluid.rho*((y(2)*A4)/Ad)^2;
end
end
dydt = zeros(4,1);
%y(1) position; y(2) velocity
dydt(1) = y(2); % xa_dot = va;
dydt(2) = (1/data.actuator.ma)*((P4*A4) - (P5*A5) - (data.actuator.F0+(data.actuator.k*y(1)))); % va_dot = (1/ma)*((P4*A4) - (P5*A5) - (F_0+(k*xa)));
dydt(3) = y(2)*A5; % Vt_dot = Qout;
dydt(4) = y(2)*A4; % Vacc_dot = -Qin;
%BOUNDARY CONDITION FOR THE PISTON
% position cant go below zero nor xmax
%if velocity is positive when x=xmax put to zero
%if negative is positive when x=xmax put to zero
%acceleration must go to zero when piston arrives at x=0 and x=xmax
if y(1) <=0
y(1)=0;
end
if y(1)==0 && dydt(1)<0
dydt(1)=0;
end
if y(1) >= data.actuator.xmax
y(1) = data.actuator.xmax;
end
if y(1) == data.actuator.xmax && dydt(1)>0
dydt(1) = 0;
end
if (y(1) == data.actuator.xmax && dydt(2) > 0) || (y(1) == 0 && dydt(2) < 0)
dydt(1) = 0;
dydt(2) = 0;
end
% when the volume of liquid inside the accumulator is bigger than the
% initial one put to zero
if y(4) >=data.accumulator.V0
dydt(4) = 0;
end
end
0 Commenti
Risposta accettata
Bjorn Gustavsson
il 1 Nov 2021
When you have discontinuities in the ODEs (the derivative of command) then the safe approach is to separate the integration into sections where everything is continuous. You might get somewhere by reducing the max-stepsize of the ode-solever by setting MaxStep in the odeoptions-struct returned from odeset to something smaller, but really think about integrating in chunks where everyting behaves nicely.
HTH
0 Commenti
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!