Why is the Temperature curve looking like this at the end of the While loop?

Hello, This code should be giving out multiple monitoring parameters, such as the droplet's diameter, the time when the droplet undergoes full evaporation, the droplet's diameter and mass loss rates, and temperature of the particle. however my problem is with the latter, the curve showing the temperature of the particle is becoming noisy at the end of the loop when the evaporation is complete but I dont know why.

%%% DEFINITIONS %%%

clear;
Tb = 20+273;         % temperature of the bulk air surrounding the particle [C]
Tp(1) = 20+273;      % temperature of the particle's surface [C]
RHo = 0.5;           % Relative Humidity air     [%]
psat(1)= exp(16.7-(4060/(Tp(1) -37)));     % saturation pressure at room temp [kpa]
pinf(1) = RHo*2.31; % partial pressure in bulk [kpa]
M= 0.018;         % molar mass of water molecule [Kg.mol-1]
Ru = 8.314;      % Universal gas constant [Kg.mol.m2/s2.K]
cinf(1) = (pinf(1)*1000)*(M) / ((Ru)*(Tb));      % vapor concentration in air at 20 degrees
RHp = 1;         %Relative Humidity at surface particle [%]
p(1) = RHp*psat(1);    % partial pressure at the surface of particle [kpa]
st = 0.0727;     % surface tension at room temp [N/m]
den = 1000;      % density at room temp [kg/m3]
dp(1)=40000*10^-9; % initial diameter of the water particle [m]
Kr(1) = exp ((4*st*M)/(Ru*den*dp(1)*Tp(1))); % kelvin effect ratio
pd(1) =p(1);        % real partial pressure at the surface of the particle [kpa]
cp(1) =(pd(1)*1000)*(M)/((Ru)*(Tp(1))); % water vapor concentration at the particle's surface [Kg/m3]
t = 10^-6    % timestep [s]

t = 1.0000e-06

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]

kair = 0.0260

Kn(1) = (2*70*10^-9)/(dp(1));
Cm(1) = (1+Kn(1))/(1+(((4/3)+0.377)*Kn(1))+((4*Kn(1)^2)/3));

%%% INITIAL CONDITIONS %%%

Tp(1) = 20+273;
Kr(1) = exp ((4*st*M)/(Ru*den*dp(1)*Tp(1)));  
psat(1)= exp(16.7-(4060/(Tp(1) -37)));
pinf(1) = RHo*2.31;
p(1) = RHp*psat(1);
pd(1) =p(1);
cp(1) = (pd(1)*1000)*(M)/((Ru)*(Tp(1)));
cinf(1) = (pinf(1)*1000)*(M) / ((Ru)*(Tb));
mp(1) = (pi/6) *den*(dp(1))^3;                  %mass particle initially    
mpdot(1)=(2*pi)*D*dp(1)*Kr(1)*Cm(1)*(cp(1)-cinf(1));           %mass rate of the particle initially 
ddot(1) = -4*D*Kr(1)*Cm(1)*(cp(1) - cinf(1))/(den*dp(1));
t2(1) = 0;                                      %time interval initially
i = 0;

%%% ITERATIONS %%%

while ddot(i+1) < 0 
    i=i+1;
    % new diameter%
    
    dp(i+1) = dp(i) + t*ddot(i);   %% new diameter every iteration 
     mp(i+1) = (den*pi/6)*(dp(i+1))^3; %[Kg]
    
    
    % new temperature of the particle %
    A(i) = 2*pi*dp(i)*kair*(Tb-Tp(i));          % Qdot
    B(i) = 2*pi*dp(i)*L*D*Kr(1)*(cp(i)-cinf(i));         % mdot 
    C(i) = den*c*pi*1/6*(1/t)*(dp(i))^3;        % mass*cp/time
    Tp(i+1) = (A(i)-B(i))/(C(i)) + (Tp(i));     % solving the differential equation to get the new particle's surface temperature each iteration 
    
    %Iterated parameters%
    Kr(i+1)= exp((4*st*M)/(Ru*den*dp(i+1)*(Tb))); 
    Kn(i+1) = (2*70*10^-9)/(dp(i+1));
    Cm(i+1) = (1+Kn(i+1))/(1+(((4/3)+0.377)*Kn(i+1))+((4*Kn(i+1)^2)/3));
    psat(i+1)= exp(16.7-(4060/(Tp(i+1) -37)));
    pinf(i+1) = RHo*2.31;
    p(i+1) = RHp*psat(i+1);
    pd(i+1) =p(i+1);
    cinf(i+1) = (pinf(i+1)*1000)*(M) / ((Ru)*(Tb));
    cp(i+1) =(pd(i+1)*Kr(i+1)*1000)*(M)/((Ru)*(Tp(i+1)));
    t2(i+1) =t2(i) + t;
    mpdot(i+1) = (pi/6)*D*dp(i+1)*Kr(i+1)*Cm(i+1)*(cp(i+1)-cinf(i+1));  
    ddot(i+1)= -4*D*Kr(i+1)*Cm(i+1)*(cp(i+1) - cinf(i+1))/(den*dp(i+1));
   
end 

%%% PLOTTING %%%

subplot (2,1,1)
plot (t2,Tp-273)
xlabel('time [s]'), ylabel('Particle Temperature [C]')
 subplot (2,1,2)
plot (t2 ,dp*10^9)
xlabel('time [s]'), ylabel('diameter [nm]')

so if you run this code, I will have two curves, one showing the temperature of the particle during evaporation, the other will be showing the diameter. in the temperature one at the end it shows like a straight line where the temperature starts increasing and decreasing abnormally, how to fix that?

[Torsten]

dp is almost zero - so you almost divide by 0 in the calculation of Tp:

C(i) = den*c*pi*1/6*(1/t)*(dp(i))^3;        % mass*cp/time
Tp(i+1) = (A(i)-B(i))/(C(i)) + (Tp(i));     % solving the differential equation to get the new particle's surface temperature each iteration 

[smith]

Alright, Is there anything can be done to avoid this noise?

[Torsten]

Sure. Don't plot the complete Tp and dp vectors.

Your update for Tp is strange. What is the underlying differential equation ?

[smith]

this.

Answer 1

[Torsten]

For clarity, I suggest the following code:

%%% DEFINITIONS %%%

clear;
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]

t = 10^-6    % timestep [s]
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]

%%% INITIAL CONDITIONS %%%

Tp(1) = 20+273;      % temperature of the particle's surface [C]
dp(1)=40000*10^-9; % initial diameter of the water particle [m]
mp(1) = (pi/6) *den*(dp(1))^3;                  %mass particle initially    
t2(1) = 0;                                      %time interval initially

i = 1;

%%% ITERATIONS %%%

while dp(i) > 4e-9        
    %Iterated parameters%
    Kr(i)= exp((4*st*M)/(Ru*den*dp(i)*(Tb))); 
    Kn(i) = (2*70*10^-9)/(dp(i));
    Cm(i) = (1+Kn(i))/(1+(((4/3)+0.377)*Kn(i))+((4*Kn(i)^2)/3));
    psat(i)= exp(16.7-(4060/(Tp(i) -37)));
    pinf(i) = RHo*2.31;
    p(i) = RHp*psat(i);
    pd(i) =p(i);
    cinf(i) = (pinf(i)*1000)*(M) / ((Ru)*(Tb));
    cp(i) = (pd(i)*Kr(i)*1000)*(M)/((Ru)*(Tp(i)));
    % new temperature of the particle %
    A(i) = 2*pi*dp(i)*kair*(Tb-Tp(i));          % Qdot
    B(i) = 2*pi*dp(i)*L*D*Kr(i)*(cp(i)-cinf(i));         % mdot 
    C(i) = den*c*pi*1/6*(1/t)*(dp(i))^3;        % mass*cp/time
    mpdot(i) = (pi/6)*D*dp(i)*Kr(i)*Cm(i)*(cp(i)-cinf(i));  
    ddot(i)= -4*D*Kr(i)*Cm(i)*(cp(i)-cinf(i))/(den*dp(i));
    %Updates
    t2(i+1) = t2(i) + t;
    Tp(i+1) = Tp(i) + (A(i)-B(i))/(C(i));     % solving the differential equation to get the new particle's surface temperature each iteration     
    dp(i+1) = dp(i) + t*ddot(i);   %% new diameter every iteration 
    mp(i+1) = (den*pi/6)*(dp(i+1))^3; %[Kg]
    i = i+1;
end 

%%% PLOTTING %%%

subplot (2,1,1)
plot (t2,Tp-273)
xlabel('time [s]'), ylabel('Particle Temperature [C]')
subplot (2,1,2)
plot (t2 ,dp*10^9)
xlabel('time [s]'), ylabel('diameter [nm]')

[smith]

But with this while condition, i will never get to see the lifetime of the particle, which is my objective, I still need the while ddot < 0. because once ddot breaks positive, this is where the particle has fully evaporated.

[Torsten]

The time when the particle has fully evaporated is when dp = 0. I just chose a slightly weaker condition to avoid division by zero in C.

[smith]

Exactly, when dp = 0 the evaporation has been attained, what I’m saying is that when ddot breaks positive, this indicates that the next diameter will be negative, logically this means that the second it breaks positive the diameter is going negative, in other words it is hitting zero. this is why I’m stopping it with that condition, I’m trying to be accurate with the lifetime simulation.

[Torsten]

Exactly, when dp = 0 the evaporation has been attained, what I’m saying is that when ddot breaks positive, this indicates that the next diameter will be negative, logically this means that the second it breaks positive the diameter is going negative, in other words it is hitting zero.

Before ddot breaks positive, dp will become zero. And this means that all your calculation will crash because of division by zero. You saw in your attempt above that this point leads to unphysical results in Tp. So just stop when dp is small enough to assume it has become 0.

this is why I’m stopping it with that condition, I’m trying to be accurate with the lifetime simulation.

You don't use any error estimates and an explicit Euler method for a highly stiff integration problem. Do you really think that making the difference between 1e-9 and 0 in particale diameter will be important compared to the errors you made during integration ?

Use ode15s to solve your 3 differential equations instead of your own solver.

[smith]

Noted, Thank you for your comments, much appreciated.

Hello Torsten, I've tried using ode15s to solve my differential equations but I'm not reaching to a solution, do you mind if you can show me how to do the mass one or temperature.

Answer 2

[Torsten]

tstart = 0;
tend = 2.735;
Tp0 = 293;
dp0 = 40000*10^(-9);
[T,Y] = ode15s(@fun,[tstart,tend],[Tp0;dp0]);
mp = pi/6*1000*Y(:,2).^3;
yyaxis left
plot(T,Y(:,1))
yyaxis right
plot(T,[Y(:,2),mp])

function dy = fun(t,y)
Tp = y(1);
dp = y(2);
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));
psat= exp(16.7-(4060/(Tp -37)));
pinf = RHo*2.31;
p = RHp*psat;
pd =p;
cinf = pinf*1000*M / (Ru*Tb);
cp = pd*Kr*1000*M/(Ru*Tp);
% new temperature of the particle %
A = 2*pi*dp*kair*(Tb-Tp);          % Qdot
B = 2*pi*dp*L*D*Kr*(cp-cinf);         % mdot 
C = den*c*pi*1/6*dp^3;        % mass*cp/time

dTp = (A-B)/C;  
ddot= -4*D*Kr*Cm*(cp-cinf)/(den*dp);

dy = [dTp;ddot];
end

[smith]

Much appreciated, I just have few questions regarding the understanding of the function.

1- at each iteration inside the function "fun", the values are changing for everything dependent on Tp and dp right?

2- what is the unit of time?

3- if i want to change the diameter of the particle dp, what other things i should be considerate to regulate along with that? because i tried changing the 4k to 1k and there was an error but i could not understand the issue behind it.

[Torsten]

1- at each iteration inside the function "fun", the values are changing for everything dependent on Tp and dp right?

The function "fun" implements the right-hand sides of your differential equations. Thus dTp and ddot are deduced from the actual Tp and dp (and maybe time itself).

2- what is the unit of time?

The unit of time is the unit you choose it to be according to your differential equations. I guess you chose "s" as unit for the physical parameters - so it's "s".

3- if i want to change the diameter of the particle dp, what other things i should be considerate to regulate along with that? because i tried changing the 4k to 1k and there was an error but i could not understand the issue behind it.

The initial diameter can be changed in the line

dp0 = 40000*10^(-9);

Since the particle diameter will become 0 earlier as for 4k, you will have to reduce "tend" in order to avoid division by zero. Make two or three attempts for "tend" and check the plot (as I did when I entered tend = 2.735).

If you want to invest more work, you can stop integration when the particle diameter becomes smaller than a given limit. For this, you will have to add an event function for ODE15S. Look up

https://de.mathworks.com/help/matlab/math/ode-event-location.html

for details on how to proceed.

[Torsten]

@smith added comment

Hello, Thank you for that explanation, I was rereading it again, I have two concerns left.

1- can you explain this line "mp = pi/6*1000*Y(:,2).^3;"

2- If i want to add another differential equation on there and new output, what is it to change in the code?

[Torsten]

1- can you explain this line "mp = pi/6*1000*Y(:,2).^3;"

It computes mp for the output times.

2- If i want to add another differential equation on there and new output, what is it to change in the code?

Add the initial value of the new component to the vector of initial values ([Tp0;dp0;newy0]), and add the derivative of the new component with respect to time to the vector of time derivatives (dy = [dTp;ddot;dnewy];).