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.

采纳的回答

Alan Stevens
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
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!
amy gravenstein
amy gravenstein 2020-11-23
This works well for calculating Evap. I would expect h (pool height) to decrease rapidly as the pool evaporates extremely quickly under most conditions. This matlab program is still a work in progress, but your help has provided much guidance!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by