Why is my ode15 differential equation solver acting like this?
2 ビュー (過去 30 日間)
古いコメントを表示
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);
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1197413/image.jpeg)
0 件のコメント
回答 (1 件)
Torsten
2022 年 11 月 17 日
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).
12 件のコメント
参考
カテゴリ
Help Center および File Exchange で Thermal Analysis についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!