Why is my ode15 differential equation solver acting like this?

I want this code to solve me three differential equations, (dTp, ddot, dm), however for the first two I'm getting acceptable results, but for the third that should be giving me the mass of the particle each timestep, is increasing, which is contradicting because according to it's relative differential equation, the mass mp should be decreasing with time.

clear all
close all
tstart = 0;
tend =2.8;
Tp0 = 293;
dp0 = 40000*10^(-9);
mp0 = (pi/6) *1000*(dp0)^3;

options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[T,Y] = ode15s(@fun,[tstart,tend],[Tp0;dp0;mp0],options);

Warning: Failure at t=2.771985e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

mass_p = pi/6*1000*Y(:,2).^3;
figure(1)
yyaxis left
plot(T,Y(:,1))
yyaxis right
plot(T,Y(:,2))

figure(2)
%yyaxis left
plot (T,[mass_p,Y(:,3)])

function dy = fun(t,y)
Tp = y(1);
dp = y(2);
mp = y(3);
Tb = 20+273;         % temperature of the bulk air surrounding the particle [C]
RHo = 0.5;           % Relative Humidity air     [%]
M= 0.018;         % molar mass of water molecule [Kg.mol-1]
Ru = 8.314;      % Universal gas constant [Kg.mol.m2/s2.K]
RHp = 1;         %Relative Humidity at surface particle [%]
st = 0.0727;     % surface tension at room temp [N/m]
den = 1000;      % density at room temp [kg/m3]
D = 2.5*10^-5;   % diffusivity coefficient of water at room temp [m2/s]
L = 2.44*10^6;   % latent heat of vaporization for water at room temp
c = 4184;       %specific heat at room temp J
kair = 0.026;     % thermal conductivity for air at room temp [W/mK]

%Iterated parameters%
Kr= exp((4*st*M)/(Ru*den*dp*Tb));
Kn = (2*70*10^-9)/dp;
Cm = (1+Kn)/(1+((4/3+0.377)*Kn)+((4*Kn^2)/3)); % this is a correction factor Fuchs factor.
psat= exp(16.7-(4060/(Tp -37)));
pinf = RHo*2.31;
p = RHp*psat;
pd =Kr*p;
cinf = pinf*1000*M / (Ru*Tb);
cp = pd*1000*Kr*M/(Ru*Tp);


dTp =(((-L*Cm*D*(cp-cinf))-(kair*(Tp-Tb)))/(den*c*(dp^(2))*(1/12))); %differential for Temperature (Tp)
ddot= -4*D*Cm*(cp-cinf)/(den*dp);  %% differential for diameter (dp)
dmp = -(2*pi)*D*Cm*dp*(cp-cinf); %% differential for mass (mp)

dy = [dTp;ddot;dmp];
end

(cs is cp)

Answer 1

[Torsten]

Seems you have an error in the sign of the time derivative for mp:

dmp = -(2*pi)*D*Cm*dp*(cp-cinf); %% differential for mass (mp)

instead of

dmp = (2*pi)*D*Cm*dp*(cp-cinf); %% differential for mass (mp)

(see above).

[smith]

Hello Torsten,

upon plotting y(3), it will not be reaching zero, however if i do plot lets assume a predefined mass,

mass_p = pi/6*1000*Y(:,2).^3;

plotting mass_P, will give me a profile that actually reaches zero. Why is it not reaching zero in the former case?

I edited the code, check what i mean

[smith]

The graph in the edited code, shows the difference between both profiles, but why are the values from the differential equation solver arent reaching zero, but the one I'm generating explicitly are reaching zero?

[Torsten]

Did you look at the scale of the y-axis ? The error between the two curves (and the distace to 0) is approximately 1e-12 ... Not an error you are willing to accept ? Then lower RelTol and AbsTol in ODE15S, and the world is all right again (see above).

[smith]

I'm willing to accept that error, It just makes me curious to know the reasoning behind it no more.

However, if i have a differential equation that includes another differential equation, just like the one attached, How do i include the dm/dt (mdot) inside the dTb/dt ?

clear all
close all
tstart = 0;
tend =3;
Tp0 = 293;
Tb0 = 293;
dp0 = 40000*10^(-9);
mp0 = (pi/6) *1000*(dp0)^3;

[T,Y] = ode15s(@fun,[tstart,tend],[Tp0;dp0;mp0;Tb0]);
mp = pi/6*1000*Y(:,2).^3;

figure(1)
grid on
yyaxis left
plot(T,Y(:,1))
xlabel('time [s]'), ylabel('Tp [C]')
yyaxis right
plot(T,Y(:,2))
xlabel('time [s]'), ylabel('dp [m]')

figure(2)
yyaxis left
plot (T,[mp,Y(:,3)])
xlabel('time [s]'), ylabel('mp [Kg]')
grid on



function dy = fun(t,y)
Tp = y(1);
dp = y(2);
mp = y(3);
Tb = y(4);         % temperature of the bulk air surrounding the particle [C]
RHo = 0.5;           % Relative Humidity air     [%]
M= 0.018;         % molar mass of water molecule [Kg.mol-1]
Ru = 8.314;      % Universal gas constant [Kg.mol.m2/s2.K]
RHp = 1;         %Relative Humidity at surface particle [%]
st = 0.0727;     % surface tension at room temp [N/m]
den = 1000;      % density at room temp [kg/m3]
D = 2.5*10^-5;   % diffusivity coefficient of water at room temp [m2/s]
L = 2.44*10^6;   % latent heat of vaporization for water at room temp
c = 4184;       %specific heat at room temp J
kair = 0.026;     % thermal conductivity for air at room temp [W/mK]
n =1; % number of particles
V = 1e-3; % volume of the air 
N = n/V; % unit per volume m-3 ( it is 1 cm3 in this case) 
denair =1.225 ; %density of air kg/m3
cpair = 1007.5; % J/kg.K







%Iterated parameters%
Kr= exp((4*st*M)/(Ru*den*dp*Tb));
Kn = (2*70*10^-9)/dp;
Cm = (1+Kn)/(1+((4/3+0.377)*Kn)+((4*Kn^2)/3));
psat= exp(16.7-(4060/(Tp -37)));
psatb = exp(16.7-(4060/(Tb -37)));
pinf = RHo*psatb;
p = RHp*psat;
pd =Kr*p;
cinf = pinf*1000*M / (Ru*Tb);
cp = pd*1000*Kr*M/(Ru*Tp);


dTp =(((-L*Cm*D*(cp-cinf))-(kair*(Tp-Tb)))/(den*c*(dp^(2))*(1/12)));
dTb = (((N*2*pi*dp*kair)*(Tp-Tb)) + (mdot*c*(Tp - Tb))) / ((denair*cpair)+(den*c))
ddot= -4*D*Cm*(cp-cinf)/(den*dp);
dmp = -(2*pi)*D*Cm*dp*(cp-cinf); 

dy = [dTp;ddot;dmp;dTb];
end

may you let me know what can be done in this case to replace this (mdot) by an actual dm/dt for the code to run? . the physical meaning for the dm/dt is the rate of mass lost every timestep.

Much appreciated torsten as always!

[smith]

the code won't run because of the identity "mdot" that I've used. I'm aware that it is not predefined, I'm just showing you it. theoritically "dmp" is what needs to be there.

[Torsten]

And why don't you insert it ?

dTp =(((-L*Cm*D*(cp-cinf))-(kair*(Tp-Tb)))/(den*c*(dp^(2))*(1/12)));
dmp = -(2*pi)*D*Cm*dp*(cp-cinf); 
dTb = (((N*2*pi*dp*kair)*(Tp-Tb)) + (dmp*c*(Tp - Tb))) / ((denair*cpair)+(den*c))
ddot= -4*D*Cm*(cp-cinf)/(den*dp);

[smith]

Right, overlooked something.

Do you have a clue how can i monitor the value of cp every iteration? is there a way to make an array of it's values with every output of dp tp etc

[Torsten]

tstart = 0;
tend =2.73;
Tp0 = 293;
Tb0 = 293;
dp0 = 40000*10^(-9);
mp0 = (pi/6) *1000*(dp0)^3;

[T,Y] = ode15s(@fun,[tstart,tend],[Tp0;dp0;mp0;Tb0]);
cp = zeros(size(T));
for i=1:numel(T)
  [~,cp(i)] = fun(T(i),Y(i,:));
end
plot(T,cp)

function [dy,cp] = fun(t,y)
Tp = y(1);
dp = y(2);
mp = y(3);
Tb = y(4);         % temperature of the bulk air surrounding the particle [C]
RHo = 0.5;           % Relative Humidity air     [%]
M= 0.018;         % molar mass of water molecule [Kg.mol-1]
Ru = 8.314;      % Universal gas constant [Kg.mol.m2/s2.K]
RHp = 1;         %Relative Humidity at surface particle [%]
st = 0.0727;     % surface tension at room temp [N/m]
den = 1000;      % density at room temp [kg/m3]
D = 2.5*10^-5;   % diffusivity coefficient of water at room temp [m2/s]
L = 2.44*10^6;   % latent heat of vaporization for water at room temp
c = 4184;       %specific heat at room temp J
kair = 0.026;     % thermal conductivity for air at room temp [W/mK]
n =1; % number of particles
V = 1e-3; % volume of the air 
N = n/V; % unit per volume m-3 ( it is 1 cm3 in this case) 
denair =1.225 ; %density of air kg/m3
cpair = 1007.5; % J/kg.K

%Iterated parameters%
Kr= exp((4*st*M)/(Ru*den*dp*Tb));
Kn = (2*70*10^-9)/dp;
Cm = (1+Kn)/(1+((4/3+0.377)*Kn)+((4*Kn^2)/3));
psat= exp(16.7-(4060/(Tp -37)));
psatb = exp(16.7-(4060/(Tb -37)));
pinf = RHo*psatb;
p = RHp*psat;
pd =Kr*p;
cinf = pinf*1000*M / (Ru*Tb);
cp = pd*1000*Kr*M/(Ru*Tp);


dTp =(((-L*Cm*D*(cp-cinf))-(kair*(Tp-Tb)))/(den*c*(dp^(2))*(1/12)));
dmp = -(2*pi)*D*Cm*dp*(cp-cinf); 
dTb = (((N*2*pi*dp*kair)*(Tp-Tb)) + (dmp*c*(Tp - Tb))) / ((denair*cpair)+(den*c));
ddot= -4*D*Cm*(cp-cinf)/(den*dp);

dy = [dTp;ddot;dmp;dTb];
end

[smith]

I tried adding this event function to stop the operation regardless of tend once the condition is met,

function [value, isterminal, direction] =eventFunction (t,y)
value = y(2);
isterminal = 1;
direction= 0;
end 

and then,

options =odeset('events',@eventFunction);

but, y(2) or dp, won't reach absolute zero, but something to the order of -23 which is zero, but the event function will not take that range, is there a way to set the event function's condition? (value of y(2) < 1% of dp1) for example or something around that

[Torsten]

value = y(2) - const

makes ODE15S stop as soon as dp = const.

You might want to choose

const = a*dp0

where a is a constant and dp0 is the initial diameter of the droplet or whatever.

[smith]

none of the parameters will be reaching a specific predicted constant that i can use, that's why i was hoping for a (+ve/-ve condition) or a specific range.

[Torsten]

The question is not if it reaches the constant exactly, but if it crosses the constant from above. The latter will be checked by ODE15S.

And this will be the case if the diameter of the droplet or whatever goes to zero.