How to solve coupled partial differential equations with method of lines?

Question

Ari Dillep Sai 2022-5-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1730075-how-to-solve-coupled-partial-differential-equations-with-method-of-lines

编辑： Torsten 2024-11-4，19:17

IMG_20220530_194950.jpg

I am getting problem in solving the partial differential equations used for modelling of packed bed adsorption column. I have attached the equations which I need to solve.

I solved it using method of lines approach with the help of some code which I got on mathworks ask community. But the graph which I am getting is different from which I need to get.

Can anyone help me in solving these equations?

I am also providing the code which I used to solve.

% parameters
eps = 0.41;             % bed porosity
volflow = 5e-7;         % volumetric flow rate in m3/sec
dia = 0.022;            % bed internal dia in m
bedarea = pi*dia^2/4;
u = volflow/bedarea;    % superficial velocity in m/s
Dl = 0.08e-4;    % axial dipsersion coefficient in m2/sec
rhop = 1228.5;  % particle density 
Kl = 0.226;     % overall mass transfer coeff in sec-1
P = 102000;     % total pressure in pascals
R = 8.314;      % gas constant in jol/mol-k
T = 500;        % gas temp in K
yco2 = 0.2;     % co2 mole fraction
cfeed = yco2*P/(R*T);       % feed conc in mol/m3
K0 = 4.31e-9;           % in pascal-1
delh = -29380;          % heat of adsorption in j/mol
Keq = K0*exp(-delh/(R*T));      % equilibrium constant in pascal-1
% using method of lines
tf = 300;      % End time
nz = 400;       % Number of gateways
np = nz + 1;    % Number of gateposts
% Initial Concentrations
C = zeros(np,1);  
q = zeros(np,1);
c0 = [C; q];   
dt = tf/300;
tspan = 0:dt:tf;
[t, c] = ode15s(@rates, tspan, c0);
% Plot results
plot(t,c(:,np-4)/cfeed),grid
xlabel('time'), ylabel('C/Cfeed')
function dcdt = rates(~,c)
eps = 0.41;             % bed porosity
volflow = 5e-7;         % volumetric flow rate in m3/sec
dia = 0.022;            % bed internal dia in m
bedarea = pi*dia^2/4;
u = volflow/bedarea;    % superficial velocity in m/s
Dl = 0.08e-4;    % axial dipsersion coefficient in m2/sec
rhop = 1228.5;  % particle density 
Kl = 0.226;     % overall mass transfer coeff in sec-1
P = 102000;     % total pressure in pascals
R = 8.314;      % gas constant in jol/mol-k
T = 500;        % gas temp in K
yco2 = 0.2;     % co2 mole fraction
cfeed = yco2*P/(R*T);       % feed conc in mol/m3
K0 = 4.31e-9;           % in pascal-1
delh = -29380;          % heat of adsorption in j/mol
Keq = K0*exp(-delh/(R*T));      % equilibrium constant in pascal-1
qm = 5.09;          % max adsrobed conc in mol/kg
n = 0.429;          % toth isotherm parameter
L = 0.171;
nz = 400;       
np = nz + 1;   
dz = L/nz;
C = c(1:np);
q = c(np+1:2*np);
dCdt = zeros(np,1);
dqdt = zeros(np,1);
%  boundary condition at z=0
 C(1) = (1/(eps*Dl+u*dz))*((eps*Dl*C(2)+u*cfeed*dz));
for i = 2:np-1
       dCdt(i) = (Dl/dz^2)*(C(i+1)-2*C(i)+C(i-1)) - (u/(2*eps*dz))*(C(i+1)-C(i-1)) - (((1-eps)*rhop)/eps)*dqdt(i);
       dqdt(i) = Kl*((qm*Keq*C(i)*R*T/(1+(Keq*C(i)*R*T)^n)^(1/n)) - q(i)); 
end
% boundary conditions at z=L
dCdt(np) = dCdt(np-1); 
% Combine rates for the function to return to the calling routine
dcdt = [real(dCdt); real(dqdt)];
end

Here I used time interms of seconds. So when I plot the graph the x-axis is also in seconds and my breakthrough time which I am getting is around 35 seconds. But in the experimental values for showed that the breakthrough time for the graph is around 300 minutes.

so please someone help whether there is any mistake which I am doing in my code.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Davide Masiello 2022-5-31

What are z+ and z- in the flux boundary condition?

Ari Dillep Sai 2022-5-31

C at z+ corresponds to the concentration after entering into the column and C at z- corresponds to the concentration before entering into the column which is inlet concentration

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Davide Masiello 2022-5-31

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1730075-how-to-solve-coupled-partial-differential-equations-with-method-of-lines#answer_974770

编辑：Davide Masiello 2022-5-31

在 MATLAB Online 中打开

Hi @Ari Dillep Sai, take a look at the code below.

The breakthrough now is found to be sometimes after 4000 s (or ~67 mins), so not the 300 mins you have suggested, but maybe a bit better than the 35 s you were getting before.

Also note the following improvements to the code:

1 - I have passed all the constants as external parameters, so you don't have to define them again in the function.

2 - I have written the method of lines without use of for loops, but with wise use of logical indexing.

3 - ode15s solves a system of DAEs now, which is necessary given the algebraic nature of the boundary conditions for the diffusion equation.

I am not sure this will work for you but it might be a good starting point to develop your final solution.

I suggest you go throught the equations and the values of the parameters again, and ensure that everything is correct.

clear,clc

% Parameters

p.eps = 0.41; % bed porosity

p.volflow = 5e-7; % volumetric flow rate in m3/sec

p.dia = 0.022; % bed internal dia in m

p.bedarea = pi*p.dia^2/4;

p.u = p.volflow/p.bedarea; % superficial velocity in m/s

p.Dl = 0.08e-4; % axial dipsersion coefficient in m2/sec

p.rhop = 1228.5; % particle density

p.Kl = 0.226; % overall mass transfer coeff in sec-1

p.P = 102000; % total pressure in pascals

p.R = 8.314; % gas constant in jol/mol-k

p.T = 500; % gas temp in K

p.yco2 = 0.2; % co2 mole fraction

p.cfeed = p.yco2*p.P/(p.R*p.T); % feed conc in mol/m3

p.K0 = 4.31e-9; % in pascal-1

p.delh = -29380; % heat of adsorption in j/mol

p.Keq = p.K0*exp(-p.delh/(p.R*p.T)); % equilibrium constant in pascal-1

p.qm = 5.09; % max adsrobed conc in mol/kg

p.n = 0.429; % toth isotherm parameter

p.L = 0.171; % bed length

% using method of lines

tf = 300*60; % End time

p.nz = 400; % Number of nodes

z = linspace(0,p.L,p.nz); % Space domain

dz = diff(z); p.dz = dz(1); % Space increment

% Initial conditions

c = zeros(p.nz,1);

q = zeros(p.nz,1);

y0 = [c;q];

tspan = [0 tf];

% Mass matrix (for DAE system)

M = eye(2*p.nz);

M(1,1) = 0;

M(p.nz,p.nz) = 0;

options = odeset('Mass',M,'MassSingular','yes');

[t,y] = ode15s(@(t,y)rates(t,y,p),tspan,y0,options);

% Plot results

plot(t/60,y(:,p.nz)/p.cfeed)

grid

xlabel('time (mins)')

ylabel('C/Cfeed')

title('Dimensionless concentration at ''z = L''')

hold on

plot(t/60,0.05*ones(size(t)),'--r')

legend('C/Cfeed','breakthrough')

function out = rates(~,y,p)

% Variables

c = y(1:p.nz,1);

q = y(p.nz+1:2*p.nz,1);

% Space derivatives

dcdz = zeros(p.nz,1);

d2cdz2 = zeros(p.nz,1);

dcdz(2:end-1) = (c(3:end)-c(1:end-2))/(2*p.dz);

d2cdz2(2:end-1) = (c(3:end)-2*c(2:end-1)+c(1:end-2))/(p.dz^2);

% Saturation

qstar = p.qm*p.Keq*c*p.R*p.T./((1+(p.Keq*c*p.R*p.T).^(p.n)).^(1/p.n));

% Governing equations

dqdt = p.Kl*(qstar-q);

dcdt = p.Dl*d2cdz2-p.u*dcdz/p.eps-p.rhop*(1-p.eps)*dqdt/p.eps;

% Boundary conditions

dcdt(1) = p.eps*p.Dl*(c(2)-c(1))/p.dz+p.u*(p.cfeed-c(1));

dcdt(end) = c(end)-c(end-1);

% Output

out = [dcdt;dqdt];

end

24 个评论
显示 22更早的评论隐藏 22更早的评论

Una 2023-1-1

编辑：Torsten 2023-1-1

在 MATLAB Online 中打开

Dear Davide,

I have the same mass transfer equation as Ari but with different isotherm. I have tried to follow your code but i get an error massage. I dont know what i am doing wrong. I hope you can help me.

clear all
clc
% Design parameters of the adsorbent bed
p.eps         = 0.87;                                    % bed porosity
%p.volflow     = 2.86;                                    % volumetric flow rate in m3/sec
%p.dia         = 5.05;                                    % bed internal dia in m
%p.bedarea     = pi*p.dia^2/4;                            % bed area in m2
p.u           = 3;                                       % superficial velocity in m/s
p.rhop        = 500;                                     % particle density in kg/m3
p.r_1         = 0.000525;                                % particle radius in m
p.r_2         = 0.000585;                                % particle radius in m                          
p.L           = 0.3;                                     % bed length in m
% General parameters
p.Sc          = 2.03e-2;                                 % schmidt number                   
p.De          = 1e-11;                                   % effective diffusivity 
p.Dm          = 1.6e-6;                                  % molecular diffusivity of co2 in air in m2/s
p.Dl          = p.Dm + (p.u*p.r_1*2)^2/(192*p.Dm);       % axial dipsersion coefficient in m2/sec
%p.rohg        = 1.2;                                     % density of air in kg/m3
%p.viskog      = 18.37e-6;                                % dynamic viscosity of air in Ns/m2
%p.Re          = 190;
%p.sh          = 2+1.1*p.Re^0.6+p.Sc^0.33;
%p.kf          = (p.sh*p.Dm)/p.r_1;
p.Kl          = 0.005;
%p.P           = 101325;                                   % total pressure in pascals
p.R           = 8.314;                                    % gas constant in J/mol-k
p.T           = 298;                                      % gas temp in K
%p.yco2        = 0.0004;                                   % co2 mole fraction
p.cfeed       = 16;                     % feed conc in mol/m3
% Parameters of the isotherm (depend on the adsorbent used)
p.qsat1       = 0.16;
p.b1          = 4.89;
p.qsat2       = 3.50;
p.b2          = 22;
p.b3          = 0.00062;
% using method of lines
tf      = 300*60;                                        % End time
p.nz    = 300;                                           % Number of nodes
z       = linspace(0,p.L,p.nz);                          % Space domain
dz      = diff(z); 
p.dz    = dz(1);                                         % Space increment
% Initial conditions
c       = zeros(p.nz,1);
q       = zeros(p.nz,1);
y0      = [c;q];   
tspan   = [0 tf];
% Mass matrix (for DAE system)
M               = eye(2*p.nz);
M(1,1)          = 0;
M(p.nz,p.nz)    = 0;
options         = odeset('Mass',M,'MassSingular','yes');
[t,y]           = ode15s(@(t,y)rates(t,y,p),tspan,y0,options);
ans = 1×2
   300   300
Error using odearguments
@(T,Y)RATES(T,Y,P) must return a column vector.

Error in ode15s (line 153)
    odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
% Plot results
plot(t,y(:,p.nz)/p.cfeed)
grid
xlabel('time (mins)')
ylabel('C/Cfeed')
title('Dimensionless concentration at ''z = L''')
hold on
function out = rates(~,y,p)
% Variables
c       = y(1:p.nz,1);
q       = y(p.nz+1:2*p.nz,1);
% Space derivatives
dcdz            = zeros(p.nz,1);
d2cdz2          = zeros(p.nz,1);
dcdz(2:end-1)   = (c(3:end)-c(1:end-2))/(2*p.dz);
d2cdz2(2:end-1) = (c(3:end)-2*c(2:end-1)+c(1:end-2))/(p.dz^2);
% Isotherm
if p.R*p.T*c < 0.16
    qstar=(p.qsat1*p.b1*p.R*p.T*c)/(1+p.b1*p.R*p.T*c);
else 
    qstar=(p.qsat2*p.b2*(p.R*p.T*c-0.16))./(1+p.b2*(p.R*p.T*c-0.16)) + p.b3*p.R*p.T*c;
end
size(qstar)
% Governing equations
dqdt = p.Kl*(qstar-q); 
dcdt = p.Dl*d2cdz2-p.u*dcdz/p.eps-p.rhop*(1-p.eps)*dqdt/p.eps;
       
% Boundary conditions
dcdt(1)     = p.eps*p.Dl*(c(2)-c(1))/p.dz+p.u*(p.cfeed-c(1));
dcdt(end)   = c(end)-c(end-1);
% Output
out = [dcdt;dqdt];
end

Una 2023-1-2

编辑：Una 2023-1-4

在 MATLAB Online 中打开

Thank you veery much Torsten!

I am not sure whether there is a continuous transation. I have taken the isotherm from a peer reviewed paper that works with the same adsorbent as I have to . But now the results seems to be reasonable. Since there is no experimental data available for that isotherm I am trying to validate the code with another isotherm. But now i am am getting another error.

clear all
clc
% Design parameters of the adsorbent bed
p.eps         = 0.4;                                     % bed porosity
%p.volflow     = 2.86;                                    % volumetric flow rate in m3/sec
%p.dia         = 5.05;                                    % bed internal dia in m
%p.bedarea     = 20; %pi*p.dia^2/4;                       % bed area in m2
p.u           = 7.06e-2;                                  % superficial velocity in m/s
p.rhop        = 880;                                     % particle density in kg/m3
p.dp          = 0.00026*2;                              % particle diameter in m
p.Rp          = 0.00026;                                  % particle radius in m
p.L           = 0.01;                                      % bed length in m
% General parameters
%p.Sc          = 2.03e-2;                                 % schmidt number                   
%p.De          = 1e-11;                                   % effective diffusivity 
p.Dm          = 1.3e-5;                                   % molecular diffusivity of co2 in air in m2/s
p.Dl          =(0.45+0.55*p.eps)*1.3e-5+0.35*p.Rp*p.u;    % axial dipsersion coefficient in m2/sec
%p.rohg        = 1.2;                                     % density of air in kg/m3
%p.viskog      = 18.37e-6;                                % dynamic viscosity of air in Ns/m2
%p.Re          = p.rohg*p.u*p.dp/p.viskog;
%p.sh          = 2+1.1*p.Re^0.6+p.Sc^0.33;
%p.kf          = (p.sh*p.Dm)/p.dp;
p.Kl          = 0.003; %1/(((p.Rp/3*p.kf)+(p.Rp^2/15*p.De))*(p.eps/1-p.eps));
%p.P           = 101325;                                  % total pressure in pascals
p.R           = 8.314;                                    % gas constant in J/mol-k
p.T           = 298.15;                                      % gas temp in K
%p.yco2        = 0.0004;                                  % co2 mole fraction
p.cfeed       = 0.016;                                    % feed conc in mol/m3
% Parameters of the isotherm (depend on the adsorbent used)
p.b0          = 93;                                   % Toth parameter in 1/bar
p.delh        = 95.3e3;                               % heat of adsorption in J/mol
p.T0          = 353.15;                               % refernece temperature in K
p.T           = 298.15;                                  % temperature in K
p.b           = p.b0*exp((p.delh/(p.R*p.T0))*((p.T0/p.T)-1));
p.ns0         = 3.40;                                 % Toth perameter in mol/kg
p.exe         = 0;                                    % dimensionless constant 
p.ns          = p.ns0*exp(p.exe*(1-p.T/p.T0));        % max adsrobed conc in mol/kg
p.t0          = 0.37;                                 % Toth parameter
p.alp         = 0.33;                                 % dimensionless constant
p.t           = p.t0+p.alp*(1-p.T0/p.T);
% using method of lines
tf      = 3000;                                        % End time
p.nz    = 400;                                        % Number of nodes
z       = linspace(0,p.L,p.nz);                       % Space domain
dz      = diff(z); 
p.dz = dz(1);                                         % Space increment
% Initial conditions
c       = zeros(p.nz,1);
q       = zeros(p.nz,1);
y0      = [c;q];   
tspan   = [0 tf];
% Mass matrix (for DAE system)
M               = eye(2*p.nz);
M(1,1)          = 0;
M(p.nz,p.nz)    = 0;
options         = odeset('Mass',M,'MassSingular','yes');
[t,y]           = ode15s(@(t,y)rates(t,y,p),tspan,y0,options);
% Plot results
plot(t,y(:,p.nz)/p.cfeed)
grid
xlabel('time (s)')
ylabel('C/Cfeed')
title('concentration at ''z = L''')
hold on
function out = rates(~,y,p)
% Variables
c       = y(1:p.nz,1);
q       = y(p.nz+1:2*p.nz,1);
% Space derivatives
dcdz            = zeros(p.nz,1);
d2cdz2          = zeros(p.nz,1);
dcdz(2:end-1)   = (c(3:end)-c(1:end-2))/(2*p.dz);
d2cdz2(2:end-1) = (c(3:end)-2*c(2:end-1)+c(1:end-2))/(p.dz^2);
% Isotherm
qstar = p.ns*p.b*c*p.R*p.T./((1+(p.b*c*p.R*p.T).^(p.t)).^(1/p.t));
% Governing equations
dqdt = p.Kl*(qstar-q); 
dcdt = p.Dl*d2cdz2-p.u*dcdz-p.rhop*(1-p.eps)*dqdt/p.eps;
       
% Boundary conditions
dcdt(1)     = p.Dl*(c(2)-c(1))/p.dz+p.u*(p.cfeed-c(1));
dcdt(end)   = c(end)-c(end-1);
% Output
out = [dcdt;dqdt];
end

Sona 2023-1-9

在 MATLAB Online 中打开

Thank you Torsten for your fast reply.

function main
cFeed= 0.017;            % Feed concentration in mol/m3
L = 3.2;              % Column length in m
t0 = 0;                % Initial Time in s            
tf= 5000;             % Final time in s
dt= 0.5;              % Time step
z = [0:0.01:L];        % Mesh generation
t = [t0:dt:tf];       % Time vector
n = numel(z);         % Size of mesh grid
%%%%%%Initial Conditions and boundary condition %%%%%%%$%
c0 = zeros(n,1);        % c = 0 at t =0 for all z
c0(1)= cFeed;             % c = cFeed at z =0 and t>0
q0 = zeros(n,1);        % t = 0, q = 0 for all z
y0 = [c0 ; q0];         
%%%%%ODE15S Solver%%%%%%
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);
%%%%%% Plots %%%%%
plot(T,Y(:,n)/cFeed)
legend('c/cFeed')
figure
plot(T,Y(:,2*n))
legend('q')
end
function DyDt=MyFun(~, y, z, n)
epsilon  = 0.5;                         % bed porosity
volflow = 2.86;                        % volumetric flow rate in m3/s
bedarea = 20;                         % bed area in m2
v= volflow/bedarea;                         % superficial velocity in m/s
rho= 55.4;                                    % particle density in kg/m3
Rp = 2.25e-3;                                 % particle radius in m
Dm = 1.3e-5;                                  % molecular diffusivity in m2/s
Dl = (0.45+0.55*epsilon)*Dm+0.35*Rp*v;        % axial dipsersion coefficient in m2/sec
k  = 3.8e-3;
R  = 8.314;                                   % gas constant in J/mol-k
T  = 298;                                     % gas temp in K
% isotherm 
b0  = 2.25e4;                                   % Toth parameter in 1/bar
delh = 6e4;                                      % heat of adsorption in J/mol
T0 = 296;                                      % refernece temperature in K
T  = 298;                                      % temperature in K
b  = b0*exp((delh/(R*T0))*(T0/T-1));
ns0 = 1.97;                                     % Toth perameter in mol/kg
exe = 2.37;                                     % dimensionless constant 
ns  = ns0*exp(exe*(1-T/T0));                    % max adsrobed conc in mol/kg
t0 = 4.22e-1;                                  % Toth parameter
alp= 9.49e-1;                                  % dimensionless constant
t = t0+alp*(1-T0/T);
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
DyDt = zeros(2*n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
for i=2:n-1
 DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1));
  D2cDz2(i) = ((c(i+1)-c(i))/(z(i+1)-z(i))-(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = 0; 
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));
% Set time derivatives at z=0
qstar = ns*b*c(1)*R*T./((1+(b*c(1)*R*T).^(t)).^(1/t));             %%%% I have substituted p with R*T*c
DqDt(1) = k*(qstar - q(1) );
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
    % q* 
 qstar = ns*b*c(i)*R*T./((1+((b*c(i)*R*T).^(t))).^(1/t));
    %Equation: dq/dt = k(q*-q)
    DqDt(i) = k*(qstar - q(i) );
      %Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*roh*dq/dt
      DcDt(i) = Dl*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i);
  end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

Torsten 2023-1-9

编辑：Torsten 2023-1-9

在 MATLAB Online 中打开

Strange curve. Seems it only knows two values: 0 and 1.9...

Maybe a mistake of units in the setting of b0 ? Shouldn't it be in unit 1/Pa ?

epsilon = 0.5; % bed porosity

volflow = 2.86; % volumetric flow rate in m3/s

bedarea = 20; % bed area in m2

v= volflow/bedarea; % superficial velocity in m/s

rho= 55.4; % particle density in kg/m3

Rp = 2.25e-3; % particle radius in m

Dm = 1.3e-5; % molecular diffusivity in m2/s

Dl = (0.45+0.55*epsilon)*Dm+0.35*Rp*v; % axial dipsersion coefficient in m2/sec

k = 3.8e-3;

R = 8.314; % gas constant in J/mol-k

T = 298; % gas temp in K

% isotherm

b0 = 2.25e4; % Toth parameter in 1/bar

delh = 6e4; % heat of adsorption in J/mol

T0 = 296; % refernece temperature in K

T = 298; % temperature in K

b = b0*exp((delh/(R*T0))*(T0/T-1));

ns0 = 1.97; % Toth perameter in mol/kg

exe = 2.37; % dimensionless constant

ns = ns0*exp(exe*(1-T/T0)); % max adsrobed conc in mol/kg

t0 = 4.22e-1; % Toth parameter

alp= 9.49e-1; % dimensionless constant

t = t0+alp*(1-T0/T);

c = 0:0.01:1;

qstar = ns*b*c*R*T./(1+(b*c*R*T).^t).^(1/t);

plot(c,qstar)

Torsten 2023-1-10

编辑：Torsten 2023-1-10

在 MATLAB Online 中打开

You should replace your function by a vectorized version.

Note that I changed

Dcdz(n) = 0

to

Dcdz(n) = (c(n)-c(n-1))/(z(n)-z(n-1));

and

b0 = 2.25e4; % Toth parameter in 1/bar

to

b0 = 2.25e4*1e-5; % Toth parameter in 1/Pa

function DyDt=MyFun(t, y, z, n) 
epsilon  = 0.5;                               % bed porosity
volflow = 2.86;                               % volumetric flow rate in m3/s
bedarea = 20;                                 % bed area in m2
v= volflow/bedarea;                           % superficial velocity in m/s
rho= 55.4;                                    % particle density in kg/m3
Rp = 2.25e-3;                                 % particle radius in m
Dm = 1.3e-5;                                  % molecular diffusivity in m2/s
Dl = (0.45+0.55*epsilon)*Dm+0.35*Rp*v;        % axial dipsersion coefficient in m2/sec
k  = 3.8e-3;
R  = 8.314;                                   % gas constant in J/mol-k
T  = 298;                                     % gas temp in K
% isotherm
b0  = 2.25e4*1e-5;                            % Toth parameter in 1/Pa
delh = 6e4;                                   % heat of adsorption in J/mol
T0 = 296;                                     % refernece temperature in K
T  = 298;                                     % temperature in K
b  = b0*exp((delh/(R*T0))*(T0/T-1));
ns0 = 1.97;                                   % Toth perameter in mol/kg
exe = 2.37;                                   % dimensionless constant
ns  = ns0*exp(exe*(1-T/T0));                  % max adsrobed conc in mol/kg
t0 = 4.22e-1;                                 % Toth parameter
alp= 9.49e-1;                                 % dimensionless constant
t = t0+alp*(1-T0/T);
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
DyDt = zeros(2*n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
% Calculate spatial derivatives
Dcdz(2:n) = (c(2:n)-c(1:n-1))./(z(2:n)-z(1:n-1));
D2cDz2(2:n-1) = ((c(3:n)-c(2:n-1))./(z(3:n)-z(2:n-1))-(c(2:n-1)-c(1:n-2))./(z(2:n-1)-z(1:n-2)))./(zhalf(2:n-1)-zhalf(1:n-2));
% Calculate D2cDz2 at z=L for boundary condition dc/dz = 0
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*Dcdz(n);
% Set time derivatives 
qstar = ns*b*c(1:n)*R*T./(1+(b*c(1:n)*R*T).^(t)).^(1/t);             %%%% I have substituted p with R*T*c
DqDt = k*(qstar - q );
DcDt(1) = 0.0;
DcDt(2:n) = Dl*D2cDz2(2:n) - v*DcDz(2:n) - ((1-epsilon)/(epsilon))*rho*DqDt(2:n);
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

Sona 2023-1-10

编辑：Torsten 2023-1-11

在 MATLAB Online 中打开

I get an error:

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in main>MyFun (line 71)

Dcdz(2:n) = (c(2:n)-c(1:n-1))./(z(2:n)-z(1:n-1));

function main
cFeed= 0.017;         % Feed concentration in mol/m3
L = 3.2;              % Column length in m
t0 = 0;               % Initial Time in s            
tf= 5000;             % Final time in s
dt= 0.5;              % Time step
z = [0:0.001:L].';    % Mesh generation
t = [t0:dt:tf];       % Time vector
n = numel(z);         % Size of mesh grid
%%%%%%Initial Conditions and boundary condition %%%%%%%$%
c0 = zeros(n,1);        % c = 0 at t =0 for all z
c0(1)= cFeed;           % c = cFeed at z =0 and t>0
q0 = zeros(n,1);        % t = 0, q = 0 for all z
y0 = [c0 ; q0];         
%%%%%ODE15S Solver%%%%%%
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);
%%%%%% Plots %%%%%
plot(T,Y(:,n)/cFeed)
legend('c/cFeed')
figure
plot(T,Y(:,2*n))
legend('q')
end
function DyDt=MyFun(t, y, z, n) 
epsilon  = 0.5;                               % bed porosity
volflow = 2.86;                               % volumetric flow rate in m3/s
bedarea = 20;                                 % bed area in m2
v= volflow/bedarea;                           % superficial velocity in m/s
rho= 55.4;                                    % particle density in kg/m3
Rp = 2.25e-3;                                 % particle radius in m
Dm = 1.3e-5;                                  % molecular diffusivity in m2/s
Dl = (0.45+0.55*epsilon)*Dm+0.35*Rp*v;        % axial dipsersion coefficient in m2/sec
k  = 3.8e-3;
R  = 8.314;                                   % gas constant in J/mol-k
T  = 298;                                     % gas temp in K
% isotherm
b0  = 2.25e4*1e-5;                            % Toth parameter in 1/Pa
delh = 6e4;                                   % heat of adsorption in J/mol
T0 = 296;                                     % refernece temperature in K
T  = 298;                                     % temperature in K
b  = b0*exp((delh/(R*T0))*(T0/T-1));
ns0 = 1.97;                                   % Toth perameter in mol/kg
exe = 2.37;                                   % dimensionless constant
ns  = ns0*exp(exe*(1-T/T0));                  % max adsrobed conc in mol/kg
t0 = 4.22e-1;                                 % Toth parameter
alp= 9.49e-1;                                 % dimensionless constant
t = t0+alp*(1-T0/T);
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
DyDt = zeros(2*n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
% Calculate spatial derivatives
DcDz(2:n) = (c(2:n)-c(1:n-1))./(z(2:n)-z(1:n-1));
D2cDz2(2:n-1) = ((c(3:n)-c(2:n-1))./(z(3:n)-z(2:n-1))-(c(2:n-1)-c(1:n-2))./(z(2:n-1)-z(1:n-2)))./(zhalf(2:n-1)-zhalf(1:n-2));
% Calculate D2cDz2 at z=L for boundary condition dc/dz = 0
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*DcDz(n);
% Set time derivatives 
qstar = ns*b*c(1:n)*R*T./(1+(b*c(1:n)*R*T).^(t)).^(1/t);            
DqDt = k*(qstar - q );
DcDt(1) = 0.0;
DcDt(2:n) = Dl*D2cDz2(2:n) - v*DcDz(2:n) - ((1-epsilon)/(epsilon))*rho*DqDt(2:n);
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

Zahra 2024-11-4，13:58

the boundary conditions should be written for dcdz not dcdt! To me, it seems that these equations need to be deiscretized at z direction and then be solved by using ode solvers. This is what I am trying to do for my model.

Torsten 2024-11-4，19:16

编辑：Torsten 2024-11-4，19:17

在 MATLAB Online 中打开

the boundary conditions should be written for dcdz not dcdt!

The boundary condition is inserted into the differential equation. Thus you get an expression for dcdt that uses the set boundary condition.

Example:

Equation

dT/dt = a*d^2T/dz^2

with

dT/dz = 0

at the right endpoint z_n.

This gives a differential equation in z_n as

dT/dt @z_n ~ 
a * (dT/dz @z_n - dT/dz @z_(n-1/2))/(dz/2) =   % Use boundary condition dT/dz @z_n = 0 here
-2*a * dT/dz @z_(n-1/2) / dz ~ 
-2*a * (T @z_n - T @z_(n-1) )/dz^2

请先登录，再进行评论。

Answer 2

Torsten 2022-5-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1730075-how-to-solve-coupled-partial-differential-equations-with-method-of-lines#answer_974420

编辑：Torsten 2022-5-30

在 MATLAB Online 中打开

The sink term

- (((1-eps)*rhop)/eps)*dqdt(i);

for the concentration of your gaseous species is always 0 because you specify dqdt(i) after dCdt(i) in the for-loop.

Thus the breakthrough is at the time instant

t_break = L/(u/eps)

4 个评论
显示 2更早的评论隐藏 2更早的评论

Ari Dillep Sai 2022-5-30

Can you please let me know the mistakes in the code?

Ari Dillep Sai 2022-5-31

编辑：Ari Dillep Sai 2022-5-31

Hi Torsten,

Can you help me in finding mistakes in my code with the respect to the boundary conditions and modelling equations?

请先登录，再进行评论。

How to solve coupled partial differential equations with method of lines?

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

24 个评论
显示 22更早的评论隐藏 22更早的评论

更多回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

How to solve coupled partial differential equations with method of lines?

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

24 个评论 显示 22更早的评论隐藏 22更早的评论

更多回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

24 个评论
显示 22更早的评论隐藏 22更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论