# 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 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=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 = 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)