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

1 次查看（过去 30 天）
smith
smith 2022-10-21
评论： Torsten 2022-10-24
Ran in:
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?
  4 个评论
显示 2更早的评论
Torsten
Torsten 2022-10-21
编辑：Torsten 2022-10-21
Sure. Don't plot the complete Tp and dp vectors.
Your update for Tp is strange. What is the underlying differential equation ?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-10-21
编辑：Torsten 2022-10-21
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 个评论
显示 3更早的评论
Torsten
Torsten 2022-10-22
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
smith 2022-10-22
编辑：smith 2022-10-23
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.

请先登录，再进行评论。

更多回答（1 个）

Torsten
Torsten 2022-10-23
Ran in:
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 个评论
显示 2更早的评论
Torsten
Torsten 2022-10-24
编辑：Torsten 2022-10-24
@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
Torsten 2022-10-24
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];).

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by