ODE45 outputs only zeros

2 次查看(过去 30 天)
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

采纳的回答

Bjorn Gustavsson
Bjorn Gustavsson 2021-11-1
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 个)

Community Treasure Hunt

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

Start Hunting!

Translated by