How to solve a differential equation system with time and space derivates? Where to write (i, j) (i-1, j) (i, j-1)?

Hi!

I am modelling the drying of granulates/particles based on the attached article and equation system which contains time and space derivates.

Based on this site I wrote a script which I think somewhat solves my system, but I don't quite understand where to write the (i, j) (i-1, j) (i, j-1)?

clear all, clc

% Timespan
tf = 100   % final time
n = 100    % number of timesteps
dt = tf/n  % timestep
tau = tf/n % timestep

%Spatial geometry
L = 2      % length
N = 5      % temperature nodes
h = L/(N-1)%length of section

% Set aside storage
H=zeros(N,n)       % Air Moisture content
Theta=zeros(N,n)   % Particle Temperature
T=zeros(N,n)       % Air Temperature
M=zeros(N,n)       % Particle Moisture content
t=zeros(1,n)

%Initial conditions
H(1,1) = 0.015
Theta(1,1) = 20
T(1,1) = 70
M(1,1) = 0.2

for i = 2:N  %space loop
    H(i,1) = 0.015
    Theta(i,1) = 20
    T(i,1) = 70
    M(i,1) = 0.2
    z(i) = (i-1)*h
    
for j = 2:n  %time loop
    
% consants and variables
A = 8.654e-5
B = 49.81
f = 1.8634
ro_p = 800
r = 0.008
a_s = 784
C_p = 1.465
C_w = 4.2e3
C_v(i-1,j) = 1859+0.236*T(i-1,j)
C_a(i-1,j) = 1003.4+0.178*T(i-1,j)
mu_a(i-1,j) = 1.691e-5+4.987e-8*T(i-1,j)-3.178e-11*(T(i-1,j))^(2)+1.319e-14*(T(i-1,j))^3
G_a = 40
alfa = 0.2755
beta = -0.34
h_a(i-1,j) = alfa*C_a(i-1,j)*G_a*a_s*(2*r*G_a/mu_a(i-1,j))^beta
h_fg(i,j-1) = (1094-0.57*Theta(i,j-1))*(1+4.35*exp(-28.25*M(i,j-1)))
K(i-1,j) = 9.0833e-5*(T(i-1,j))^2 + 2.9333e-3*T(i-1,j)
RH(i-1,j) = 1.626*H(i-1,j)/(1+1.608*H(i-1,j))*10^(5) / (6892.4*exp(54.633-12301.69/(1.8*T(i-1,j)+492)-5.169*log(1.8*T(i-1,j)+492)))
%M_e(i-1,j) = (-log(1-RH(i-1,j))/(A*(T(i-1,j)+B)))^(1/f)
M_e(i-1,j)=0.015

% differential equations 
dMdt = K(i-1,j)*T(i-1,j)*((M_e(i-1,j)-M(i,j-1))/M(i,j-1))
dHdt = (-ro_p/(G_a)*(M(i,j)-M(i,j-1)))/(2*h)
dThetadt = h_a(i-1,j)*(T(i-1,j)-Theta(i,j-1)) / (ro_p*C_p+ro_p*C_w*M(i,j-1)) - (h_fg(i,j-1)+(T(i,j)-Theta(i,j-1))*C_v(i-1,j)*G_a*(H(i,j)-H(i-1,j)))/((ro_p*C_p+ro_p*M(i,j-1)*C_w)*tau)
dTdt = -h_a(i-1,j)*(T(i-1,j)-Theta(i,j-1))/(G_a*C_a(i-1,j)+G_a*H(i-1,j)*C_v(i-1,j))

M(i,j) = M(i-1,j-1) + dt*dMdt
H(i,j) = H(i,j-1) + dt*dHdt
Theta(i,j) = Theta(i,j-1) + dt*dThetadt
T(i,j) = T(i,j-1) + dt*dTdt
end
end

figure(1)
subplot(2,2,1)
plot(1:n,Theta(2,:))
title('Air Temperature')
xlabel('t [s]');
ylabel('Theta [°C]');
subplot(2,2,2)
plot(1:n,T(2,:))
title('Particle Temperature')
xlabel('t [s]');
ylabel('T [°C]');
subplot(2,2,3)
plot(1:n,H(2,:))
title('Moisture of air')
xlabel('t [s]');
ylabel('H [kg/kg]');
subplot(2,2,4)
plot(1:n,M(2,:))
title('Moisture of particle')
xlabel('t [s]');
ylabel('M [kg/kg]');

Here are the equations:

And their solution (Btw the article solves the M by Runge-Kutta method, but I didn't understand because M doesn't contain (i,j) paramaters like the other properties, so I didn't solve whit it):

• If a differential equation contains a variable which is the function of that property (like the specific heat of the air (C_a) is the function of the Temperature (T)) than how to write C_a and T?

C_a(i,j) = 1003.4+0.178*T(i,j)       %or
C_a(i-1,j) = 1003.4+0.178*T(i-1,j)   %or
C_a(i,j-1) = 1003.4+0.178*T(i,j-1)   %???

Since in the article T is a derivate of space then write T(i-1,j) and if another properties like the Particle Temperature (theta) is a derivate of time then write Theta(i, j-1) after the constants, or just simply (i, j)?

• I use the Eulure integration for solving the system and in the case of Air temperature (T) I calculate dT. And after that what do I multiply the dT by? The change of time (dt) or the change of spatial value (h?). And for the ather properties like Particle temperature which is time derivate?

T(i,j) = T(i,j-1) + dt*dTdt    %or
T(i,j) = T(i,j-1) + h*dTdt     %??

• Lastly, I have to add the dt*dT (or h*dT) to the previeous value of the T. But how to write that previous value? And for the time derivative properties like Particle Temperature (theta)?

T(i,j) = T(i,j-1) + dt*dTdt    %or
T(i,j) = T(i-1,j) + dt*dTdt    %or
T(i,j) = T(i-1,j-1) + dt*dTdt  %??

Here is my script so far, ignore the constants and calculation the only thing I don't understand is that (i, j) (i, j-1) (i, j-1) thing. Thank you in advance.

Answer 1

[Sulaymon Eshkabilov]

First, that would be nice for simulation speed to suppress to display the results. To do that you'd need to put semicolon at the end of all calc lines and value assigned commands.

Second, put all constant values outside of the loop:

% consants and variables
A = 8.654e-5;
B = 49.81;
f = 1.8634
ro_p = 800;
r = 0.008;
a_s = 784;
C_p = 1.465;
C_w = 4.2e3;
G_a = 40;
alfa = 0.2755;
beta = -0.34;

That speeds up the simulation.

Here are the correct indexes (i, j):

C_a(i-1,j) = 1003.4+0.178*T(i-1,j);   % Correct one
T(i,j) = T(i,j-1) + dt*dTdt;          % Correct one

Note that all simulations start within loops starting at 2 and thus, i-1 for L space and j-1 for time space.

[Jani]

Thank you for you answer!

And if I understand correctly for a time derivate like Theta, the correct answer would be?:

C_a(i,j-1) = 1003.4+0.178*Theta(i,j-1);
Theta(i,j) = T(i-1,j) + h*dThetadt;

So in the second row, the 'dt' becomes 'h'? And the -1 swaps place?