How to add 2 other equations which depends on Tp change in the ODE.
1 view (last 30 days)
Show older comments
These two equations (psatv =exp (A-B/(Tp+C))*133.322 and Cvs = (psatv*Nl)/(1000*R*Tp) have parameter Tp that solves by ode3 equation. My questions is how to write the coding for these two equations that depends on the solution of the ode3.
My full coding so far:
clc;
clear;
close all;
%parameters
R = 8.314; %gas constant J/mol K - the temperature must be absolute!
Tp = 298.15; %Initial temperature of the droplet/particle K
Nl = 18.015; %Molecular weight of the water, g/mol
A = 18.3036; %Antoine constant
B = 3816.44; %Antoine constant
C = -46.13; %Antoine constantTps = 304.4; %Temperature of particle at surface, Tps=Twb Are you sure?
Dp = 0.00152; %Diameter of droplet/particle, m
rp = Dp/2; %Radius of the droplet/particle, m
vp = 0.75; %Velocity of droplet/particle, m/s
rhog = 1.028; %Density of drying gas, kg/m3
ug = 2.052*10^-5; %Viscosity of drying gas, kg/m s
Dv = 2.78*10^-5; %Vapour air diffusion coefficient, m2/s
Cvi = 0.01024; %Vapour concentration in the drying gas, kg/m3
kg = 0.03; %Thermal conductivity of drying gas, W/m K
cpp = 4009; %Specific heat capacity of the droplet/particle, J/kg K
Tg = 343.15; %Temperature of drying gas, K
hfg = 2.33*10^6; %Latent heat of vapourisation, J/kg
rhod = 1052; %Density of the droplet, kg/m3
rhol = 997; %Density of the water, kg/m3
Dls = 1.5*10^-10; %Diffusivity of water into solid,m2/s
%Eqns
psatv =exp (A-B/(Tp+C))*133.322;
Cvs = (psatv*Nl)/(1000*R*Tp);
Re = (Dp*vp*rhog)/ug;
Sc = ug/(rhog*Dv);
Pr = cpp*ug/kg;
kc = (2+(0.6*Re^0.5*Sc^1/3))*Dv/Dp;
alpha = (2+(0.6*Re^0.5*Pr^1/3))*kg/Dp;
mp = rhod*(4*pi*(rp)^3)/3;
%Define the equations
%ode1 = @(t,y) -4*pi*rp^2*kc*(Cvs-Cvi);
%ode2 = @(t,y) ode1/(4*pi*rhol*rp^2);
%ode3 = @(t,y)(-4*pi*rp^2*alpha*t*(Tp-Tg)-hfg*(ode1))/(mp*cpp);
%f=[ode1;ode2;ode3]
format long
%f= @(t,y) [-4*pi*rp^2*kc*(Cvs-Cvi);
% (-4*pi*rp^2*kc*(Cvs-Cvi))/(4*pi*rhol*rp^2);
% (-4*pi*rp^2*alpha*(t)*(Tp-Tg)-hfg*(-4*pi*rp^2*kc*(Cvs-Cvi)))/(mp*cpp)]
% The functions
%y(1)=ml
%y(2)=rp
%y(3)=Tp
f= @(t,y) [-4*pi*(rp)^2*kc*(Cvs-Cvi);
(-4*pi*(rp)^2*kc*(Cvs-Cvi))./(4*pi*rhol*(y(2))^2);
(-4*pi*(rp)^2*alpha*(t)*(y(3)-Tg)-hfg*(-4*pi*(rp)^2*kc*(Cvs-Cvi)))/(mp*cpp)]
% time discretization
dt = 0.1;
t0 = 0;
tf = 100;
h=dt;
t = t0:dt:tf;
tspan=[t0 tf];
%Define Initial conditions
cond1 = 1.93e-6;
cond2 = 0.00076;
cond3 = 298.15;
y0=[cond1;cond2;cond3];% the initial value of the vector y
%Call the Euler function
[ T, Y ] = EULER( f, t0, tf, y0, h);
disp([' ' 'time' ' ' 'Ml (Y(1))' ' ' 'rp (Y(2))' ' ' 'Tp (Y(3))'])
disp([T(:) Y(:,1) Y(:,2) Y(:,3)])
figure(1)
plot(t,Y(:,1), 'r', 'linewidth', 2)
title('Droplet Mass versus Time')
xlabel('Time,t')
ylabel('Mass,kg')
grid on
figure(2)
plot(t,Y(:,2), 'g', 'linewidth', 2)
title('Droplet Radius versus time')
xlabel('Time,t')
ylabel('Droplet size,m')
grid on
figure (3)
plot(t,Y(:,3), 'b', 'linewidth', 2)
title('Temperature versus time')
xlabel('time,t')
ylabel('Temperature')
grid on
EULER
function [ T, Y ] = EULER( f, t0, tf, y0, h)
% [ T, Y ] = EULER( f, t0, tf, y0, h)
% uses Euler's method to solve a system of differential equations
% It's evaluated by f Starting at time instant t0 with initial condition,
% x0, a solution is approximated at time tf moving in time step h.
% T contains the evaluation times and
% Rows of Y contain the corresponding estimates for y(t).
T = t0 : h : tf;
Y = zeros(length(T),length(y0));
Y(1,:) = y0';
for i = 2 : length(T)
Y(i,:) = Y(i-1,:)+f(T(i-1),Y(i-1,:))'*h;
end
end
0 Comments
Answers (1)
Nipun
on 23 May 2024
Hi Faridatul,
I understand that you are trying to solve a system of ordinary differential equations where the parameters psatv and Cvs depend on the solution of the ODE, specifically the temperature Tp.
I recommend using a MATLAB function such as "ode45" instead of the custom "Euler" method for better accuracy.
For more information on "ode45" MATLAB function, refer to the following MathWorks documentation: https://in.mathworks.com/help/matlab/ref/ode45.html#d126e1117222
Hope this helps.
Regards,
Nipun
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!