Why is the Temperature curve looking like this at the end of the While loop?
1 view (last 30 days)
Show older comments
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]
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]
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?
4 Comments
Accepted Answer
Torsten
on 21 Oct 2022
Edited: Torsten
on 21 Oct 2022
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]')
5 Comments
Torsten
on 22 Oct 2022
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.
More Answers (1)
Torsten
on 23 Oct 2022
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
4 Comments
Torsten
on 24 Oct 2022
Edited: Torsten
on 24 Oct 2022
@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
on 24 Oct 2022
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];).
See Also
Categories
Find more on Thermal Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!