How to discretize a system of spatial and temporal differential equation for mass and heat transfer whit the Euler method?

5 views (last 30 days)
Jani
Jani on 11 Nov 2021
I want to make a model for the equations which are found in the attached article, and have some questions regarding them.
Here are the main second order, partial differential equations
And here them after I discretized them.
My 1) question: Is the discretization correct? And if so, than my code correctly describes it? What I mean is for delta t, and delta x, I should use dt and dx (the time and space increments)?
My code:
for j = 2:m-1 % Time Loop
for i= 2:n-1 % Space Loop
X(i,j+1) = D*((X(i+1,j)-2*X(i,j)+X(i-1,j))/(dx^2)+2/(r)*(X(i+1,j)-X(i-1,j)/(2*dx)))*dt + X(i,j)
T(i,j+1) = lambda_p*((T(i+1,j)-2*T(i,j)+T(i-1,j))/(dx^2)+2/(r)*(T(i+1,j)-T(i-1,j)/(2*dx)))*dt/(ro_p*c_p) + T(i,j)
end
end
My 2) question is, that some variables like p_sat, contains the result of the Temperature (T).
How to write it correctly? I tried this way, but I'm almost sure, that T(n,m) is not correct, but don't know how to write it.
p_sat = 1e3*10^(7.19621-1730.63/(233.426+(T(n,m)-273.15)));
My 3) question is that boundary conditions look like this
And I'm not sure if I wrote them correctly into my code, like this (the initials condition is that the Mositure contant (X) and the Temperature (T) equals to the initial X and T)?
% Initial and boundary conditions
X = zeros(n,m);
T = zeros(n,m);
X(1,1) = X_in;
T(1,1) = T_in;
% Boundary conditions
X(1,:) = 0;
X(n,:) = -k*(c_s-c_f)/(D*ro_p);
T(1,:) = 0;
T(n,:) = -(h*(T(n,m)-T_f) + k*(c_s-c_f)*dH_vap)/lambda_p;
Here is my full code which of course not runs :(
clear variables
close all
clear all
% 1. Parameters for the space
n = 20; % number of space steps
Lx = 3.3e-3; % max space (= radius)
dx = Lx/(n-1); % space step
x = 0:dx:Lx; % space increment
% 2. Parameters for the time
m = 100; % number of timesteps
tf = 100; % final time
dt = tf/(m-1); % time step
t = 0:dt:tf; % time increment
% Constants
r = 3.3e-3; % radius of particle (m)
lambda_p = 0.07; % Thermal conductivity of particle (W/mK)
T_in = 293; % Initial particle temperature (K)
X_in = 0.082; % Initial particla moisture content (kg H2O / kg dry matter)
dH_vap = 2.790e6; % Latent heat of vaporization (J/kg)
c_f = 0.01; % Drying gas water concentration (kg/m3)
T_f = 330; % Drying air temperature (K)
ro_p = 800; % Density of particle (kg/m3)
c_p = 1000; % Heat capacity of particle (J/kgK) %% ez egy változó amúgy
D = 3e-9; % Water diffusivity (m2/s)
k = 0.01; % Mass transfer coefficient (m/s)
h = 400; % Convective heat transfer coefficient (W/m2K)
M = 18; % Molar mass of water (kg/kmol)
R = 8314; % Universal gas constant (J/kmolK)
a_w = 0.7 % Water activity at the surface (-)
% Initial and boundary conditions
X = zeros(n,m);
T = zeros(n,m);
X(1,1) = X_in;
T(1,1) = T_in;
% Variables
p_sat = 1e3*10^(7.19621-1730.63/(233.426+(T(n,m)-273.15)));
c_sat = p_sat*M/(R*T(n,m));
c_s = a_w*c_sat;
% Boundary conditions
X(1,:) = 0;
X(n,:) = -k*(c_s-c_f)/(D*ro_p);
T(1,:) = 0;
T(n,:) = -(h*(T(n,m)-T_f) + k*(c_s-c_f)*dH_vap)/lambda_p;
% Implementation of the explicit method
for j = 2:m-1 % Time Loop
for i= 2:n-1 % Space Loop
X(i,j+1) = D*((X(i+1,j)-2*X(i,j)+X(i-1,j))/(dx^2)+2/(r)*(X(i+1,j)-X(i-1,j)/(2*dx)))*dt + X(i,j)
T(i,j+1) = lambda_p*((T(i+1,j)-2*T(i,j)+T(i-1,j))/(dx^2)+2/(r)*(T(i+1,j)-T(i-1,j)/(2*dx)))*dt/(ro_p*c_p) + T(i,j)
end
end
figure (1)
hold all
for i=1:4:numel(t)
plot(x,X(:,i),'linewidth',1.5,'DisplayName',sprintf('t = %1.2f',t(i)))
end
a = ylabel('Moisture');
a = xlabel('Radius (m)');
a=title('Moisture change in function of radius');
legend('-DynamicLegend','location','bestoutside');
grid;
figure (2)
[X, T] = meshgrid(x,t);
s2 = mesh(X',T',X);
set(s2,'FaceColor',[0 0 1],'edgecolor',[0 0 0],'LineStyle','--');
a = title('Moisture distribution in the function of Time and Radisus');
a = xlabel('Radius (m)');
a = ylabel('Time (s)');
Sorry for the long question, but I think they can be answered quickly if someone understands programming in MatLab better than I do.
(I tried to use pdepe, but whit it the results were crap)

Answers (0)

Community Treasure Hunt

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

Start Hunting!