Why is my ode15 differential equation solver acting like this?

Question

smith 2022-11-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1854293-why-is-my-ode15-differential-equation-solver-acting-like-this

编辑： Torsten 2022-11-21

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)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2022-11-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1854293-why-is-my-ode15-differential-equation-solver-acting-like-this#answer_1103083

在 MATLAB Online 中打开

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 个评论
显示 10更早的评论隐藏 10更早的评论

smith 2022-11-18

在 MATLAB Online 中打开

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!

Torsten 2022-11-19

编辑：Torsten 2022-11-20

在 MATLAB Online 中打开

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 2022-11-21

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.