49 views (last 30 days)

Show older comments

Follow up from my previous post about my ODE system using the runge-kutta 4th order method.

The main body of the code looks like this:

(it could be far cleaner I'm sure, but I am a novice at MATLAB so I didn't want to change much in case It results in errors)

fD_p=@(t,D_p,Tp,u_p) (-2.*kc(t,D_p,u_p).*(C_s(t,Tp)-C_inf))./rhoP;%differential eqn for diameter change of droplet

%dT_p/dt

fTp=@(t,D_p,Tp,u_p) (6.*h_m(t,D_p,u_p)*(Ta-Tp)-(6.*L(t,Tp).*kc(t,D_p,u_p).*(C_s(t,Tp)-C_inf)))./(rhoP*cpP.*D_p);% differential eqn for temperature change of droplet

%du_p/dt, velocity differential "u_p"

fu_p=@(t,D_p,Tp,u_p) -18*muG.*u_p./(rhoP.*D_p.^2)+g.*(1-(rhoG./rhoP)); % differential eqn for velocity of droplet

fh_p=@(t,u_p) u_p; % differential for vertical height travelled as a function of time and velocity

%step size

h=0.0001; % step size

tfinal=18.4; % final time

N=tfinal/h; % number of iterations

%update loop

for i=1:N

%update time

t(i+1)=t(i)+h;

%update D_p,T_p and u_p

%Runge-Kutta 4th order code

k1D_p= fD_p(t(i), D_p(i), Tp(i),u_p(i));

k1Tp= fTp(t(i), D_p(i),Tp(i),u_p(i));

k1u_p= fu_p(t(i), D_p(i), Tp(i),u_p(i));

k1h_p=fh_p(t(i),u_p(i));

k2D_p= fD_p(t(i)+h/2, D_p(i)+h/2*k1D_p, Tp(i)+h/2*k1Tp, u_p(i)+h/2*k1u_p);

k2Tp= fTp( t(i)+h/2, D_p(i)+h/2*k1D_p, Tp(i)+h/2*k1Tp, u_p(i)+h/2*k1u_p);

k2u_p= fu_p(t(i)+h/2, D_p(i)+h/2*k1D_p, Tp(i)+h/2*k1Tp, u_p(i)+h/2*k1u_p);

k2h_p= fh_p(t(i)+h/2, u_p(i)+h/2*k1u_p);

k3D_p= fD_p(t(i)+h/2, D_p(i)+h/2*k2D_p, Tp(i)+h/2*k2Tp, u_p(i)+h/2*k2u_p);

k3Tp= fTp( t(i)+h/2, D_p(i)+h/2*k2D_p, Tp(i)+h/2*k2Tp, u_p(i)+h/2*k2u_p);

k3u_p= fu_p(t(i)+h/2, D_p(i)+h/2*k2D_p, Tp(i)+h/2*k2Tp, u_p(i)+h/2*k2u_p);

k3h_p= fh_p(t(i)+h/2, u_p(i)+h/2*k2u_p);

k4D_p= fD_p(t(i)+h,D_p(i)+h*k3D_p,Tp(i)+h*k3Tp, u_p(i)+h*k3u_p);

k4Tp= fTp(t(i)+h, D_p(i)+h*k3D_p,Tp(i)+h*k3Tp, u_p(i)+h*k3u_p);

k4u_p= fu_p(t(i)+h, D_p(i)+h*k3D_p,Tp(i)+h*k3Tp, u_p(i)+h*k3u_p);

k4h_p= fh_p(t(i)+h, u_p(i)+h*k3u_p);

%update D_p value

D_p(i+1)=D_p(i)+h/6*(k1D_p+ 2*k2D_p + 2*k3D_p +k4D_p);

%update Tp value

Tp(i+1)=Tp(i)+h/6*(k1Tp+ 2*k2Tp + 2*k3Tp +k4Tp);

%update u_p value

u_p(i+1)=u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p);

%update h_p value

h_p(i+1)=h_p(i)+h/6*(k1h_p+ 2*k2h_p + 2*k3h_p +k4h_p);

end

The other vairables in the functions which call at time t,such as

h_m(t,D_p,u_p) % heat transfer coefficient, depends on the diameter and velocity of particle at time t.

..are calculated outside the R4K loop as they arent differentials, they depend on the parameters at each time step t.

Now on to the issue at hand. If you look at the velocity differential, eventually it starts to decrease and get close to 0 which is expected. The issue is, the iteration keeps going past 0, into the negatives, resulting in the other 3 differentials to break. Is there a way to stop the iteration for it,such that the differential becomes a constant at or as close to 0 as possible?

I cant assign a loop for each differential since they all depend on each other,so I'm a bit confused on what to do..

Many thanks

Abdu

Alan Stevens
on 27 Feb 2021

Have you tried replacing

u_p(i+1)=u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p);

with

u_p(i+1)= max(u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p),0);

Alan Stevens
on 28 Feb 2021

Try setting

flag = 1;

before your for loop, then have

if (u_p(i)>0) && (flag == 1)

u_p(i+1)=u_p(i)+h/6*(k1u_p+ 2*k2u_p + 2*k3u_p +k4u_p);

else

u_p(i+1) = 0;

flag = 0;

end

where you calculate u_p(i+1).

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

Start Hunting!