ODE stop integration with imaginary numbers
2 次查看(过去 30 天)
显示 更早的评论
Hello! I am modeling a cryogenic liquid spill in order to predict evaporation rate, pool mass accumulation, pool temperature, and pool radius. I incorporated a spill time (tspill) in addition to an overall run time to represent the time it takes until some can shut off the spill/leak. However, if this tspill is 10 seconds and I run the program anything more than maybe 40 seconds, it gets stuck calculating imaginary numbers in one of my for loops. I want to be able to stop the ODE integration before it begins calculating imaginary numbers. I have located this instance which is when either the mass term or radius term becomes less than zero.
I tried using these two formats to stop the integration:
function[value,isterminal,direction] = spillmodel1(t,Y)
value = [Y(1)<0; Y(2); Y(3)<0; Y(4)]; %stop when
isterminal = 1; %stop integration
direction = 0;
M=Y(1); %mass (kg)
Tpool=Y(2); %Tpool (K)
r=Y(3); %radius (m)
Evap = Y(4); %Evaporation (kg)
I think I struggle here to set the value terms to identify imaginary numbers.
And this when I call the function:
opts = odeset('Events',@spillmodel1);
[t,y]=ode45(@spillmodel1, tspan, yo, opts);
When I do this the code runs but it calculates incorrect data for all 4 y terms (they increase expontentially which is wrong). If i take this format out, it calculates correctly but the imaginary numbers are present in the data.
I will attach the code as it runs properly under certain conditions (run time less than 40 seconds, appropriate spill rate) without the ode stop integration formatting. If anyone can identify the error in my ode stop integration format please let me know.
Thank you.
0 个评论
采纳的回答
Alan Stevens
2020-11-18
Are you sure you have the units right here when dMdt <=0?
if dMdt > 0
drdt=sqrt(2*g*(h-hmin));
else
drdt = -M./(density.*pi*h*r);
end
When dMdt >0 the units of drdt are m/s. However, for the other condition it looks like drdt has units of m.
11 个评论
Alan Stevens
2020-11-20
You need to recalcuate it outside the function, from all the other parameters. I've attached a very coarse way of doing this. I've also added a plot of the height,h, which gets very small extremely rapidly. I've no feel for the validity or otherwise of any of this!
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!