ODE stop integration with imaginary numbers

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.

Answer 1

[Alan Stevens]

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.

[amy gravenstein]

Your observation is right... I added a check point to this to see which condition it runs and am seeing that is runs the "else" condition. So, I still encounter the imaginary number issue. I also adjusted the first equation to allow for units of meters only and the values are different (I am looking further into this to make sure it is accurate) but I still encounter the same problem.

[Alan Stevens]

Surely, both should have units of m/s, since drdt is a length per unit time.

[amy gravenstein]

I see! You are definitely right. Adding t to the denominator of the else condition slows the rate in which radius decreases as evaporation rate surpasses spill rate. This looks much better to me. However, I still cannot stop the ODE with the current format I am using...

[Alan Stevens]

Simply including t in the denominator sounds rather arbitrary! What is the derivation of the rate?

Wht do you mean by not being able to stop the ODE?

[amy gravenstein]

Change in radius can be directly related to the change in mass and density as the scenario plays out. I use those terms to calcuate the corresponding radius of the pool (as if it was a cylinder). This equation is still a work in progress, but my main issue is not being able to terminate the ODE after the mass term becomes negative. I gave the format I tried using in my original question. I need to stop it in order to prevent accumulating imaginary numbers.

[Alan Stevens]

Assuming rate of decrease of radius is proportional to rate of decrease of surface area when dMdt is negative you could reason that dA/dt = (dM/dt)/density*A/V where A is surface area and V is volume. Assuming that just the curved area of the cylinder is significant then 2pi*h*dr/dt = (dM/dt)/density*2pi*h*r/(pir^2h) or dr/dt = (dM/dt)/(density*pi*r*h).

Where you calculate dM/dt you could include something like

if h<hmin
    dMdt = 0;
else
    ...
end

This might help prevent the occurrence of complex numbers.

[amy gravenstein]

I did incorporate this but did not see the results I would need, probably because there are many ways in which dMdt can become imaginary. Is there a format I can use around the function to stop solving when a number becomes imaginary? I thought what I had in my original question would work but I think I have something wrong...

[Alan Stevens]

I've noticed another oddity in spillodel1. On entry you assign Y(4) to Evap, but on exit you set dY(4)/dt to Evap. This means Evap is interpreted differently between entry and exit.

[amy gravenstein]

is there a way to report Evap in spillmodel1 without putting it into the differential?

[Alan Stevens]

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]

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!