# Why is my ode15 differential equation solver acting like this?

3 Ansichten (letzte 30 Tage)
smith am 17 Nov. 2022
Bearbeitet: Torsten am 21 Nov. 2022
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 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Antworten (1)

Torsten am 17 Nov. 2022
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)
dmp = (2*pi)*D*Cm*dp*(cp-cinf); %% differential for mass (mp)
(see above).
##### 12 Kommentare10 ältere Kommentare anzeigen10 ältere Kommentare ausblenden
smith am 21 Nov. 2022
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.
Torsten am 21 Nov. 2022
Bearbeitet: Torsten am 21 Nov. 2022
The question is not if it reaches the constant exactly, but if it crosses the constant from above. The latter will be checked by ODE15S.
And this will be the case if the diameter of the droplet or whatever goes to zero.

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Thermal Analysis finden Sie in Help Center und File Exchange

R2018a

### Community Treasure Hunt

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

Start Hunting!

Translated by