Solving coupled Partial Differential Equations using method of Lines

I am working to solve mass balance equations for 4 connected chemical process control units for a continuous operation. The resulting one dimensional coupled and (non-coupled) Partial Differential Equations are aimed to be solved using Method of Lines with finite difference method. One process unit includes a coupled PDE. The general form of the equations follows convection – diffusion format (see attached pdf)

I came across the solution suggested by Torsten in the following question https://de.mathworks.com/matlabcentral/answers/371313-error-in-solving-system-of-two-reaction-diffusion-equations where a couple PDE-ODE is solved. I tried to go through the literature and extend the solution to my system of equations. To bench mark the solution approach as suggested, I tried to reproduce the result for the coupled system of PDEs from the book ‘A numerical primer for the Chemical Engineer’ by Edwin Zonderwan, CRC press {second edition, p 159}.

The code uses the pdepe solver and is attached alongwith, for the following equations (refer pdf attached).

The results from both the approaches differ significantly towards the end iterations. Could anyone suggest me what am I missing or any restriction on the solution approach? I am not sure how the second derivative is approximated at the first discrete as well. Thank you for all your help.

% Method of lines approach for four coupled pdes (refer pdf)
clear;clc
L = 10; % length of reactor
x = linspace(0,L,500); % discretise the domain in 500 discretes
tspan = linspace(0,20,2000); % given time span

n = numel(x);

% initial conditions
u0 = ones(n,1); % concentration of A
v0 = ones(n,1); % concentration of B
c0 = zeros(n,1); % concentration of C
d0 = zeros(n,1); % concentration of D

y0 = [u0;v0;c0;d0];

% setting up the solver 
[T,Y] = ode15s(@(t,y) numericalprimertubular(t,y,x,n),tspan,y0);

% saving results
myCA = Y(:,1:(n-10));
myCB = Y(:,n+1:(2*n)-10);
myCC = Y(:,2*n+1:(3*n)-10);
myCD = Y(:,3*n+1:(4*n)-10);

Using the pdepe solver for the same example

% Code from the book, 'A numerical primer for the Chemical
% Engineer' by Edwin Zonderman.
clc;clear;close all
z = linspace(0,10,30); % discretising the domain
t = linspace(0,20,30); % given time span
sol = pdepe(0,@pdex2pde,@pdex2ic,@pdex2bc,z,t); % calling the solver

% saving results
CA = sol(:,:,1); CB = sol(:,:,2); CC = sol(:,:,3); CD = sol(:,:,4);

function definitions:

% function definition for method of lines solution
function DyDt = numericalprimertubular(t,y,x,n)

vel = 1; % given velocity
Dax = 1e-4; % axial dispersion coefficient
k1 = 1; k2 = 1; % kinetic constants

u = zeros(n,1); % concentration of A
v = zeros(n,1); % concentration of B
c = zeros(n,1); % concentration of C
d = zeros(n,1); % concentration of D

DuDt = zeros(n,1);
DvDt = zeros(n,1);
DcDt = zeros(n,1);
DdDt = zeros(n,1);

DyDt = zeros(4*n,1);

xhalf = zeros(n-1,1);

DuDx = zeros(n-1,1);
DvDx = zeros(n-1,1);
DcDx = zeros(n-1,1);
DdDx = zeros(n-1,1);

D2uDx2 = zeros(n-1,1); 
D2vDx2 = zeros(n-1,1);
D2cDx2 = zeros(n-1,1);
D2dDx2 = zeros(n-1,1);


u = y(1:n);
v = y(n+1:2*n);
c = y(2*n+1:3*n);
d = y(3*n+1:4*n);


xhalf(1:n-1)=(x(1:n-1)+x(2:n))/2;

% danckwerts boundary condition
DuDx(1) = (u(1)-1)*(vel/Dax); 
DvDx(1) = (v(1)-1)*(vel/Dax);
DcDx(1) = (c(1)-0)*(vel/Dax);
DdDx(1) = (d(1)-0)*(vel/Dax);

% not sure about the applicability of these conditons !
D2uDx2(1) = 1.0/(xhalf(1)-x(1))*(u(2)-u(1))/(x(2)-x(1)); 
D2vDx2(1) = 1.0/(xhalf(1)-x(1))*(v(2)-v(1))/(x(2)-x(1)); 
D2cDx2(1) = 1.0/(xhalf(1)-x(1))*(c(2)-c(1))/(x(2)-x(1));
D2dDx2(1) = 1.0/(xhalf(1)-x(1))*(d(2)-d(1))/(x(2)-x(1));

for i=2:n-1
  DuDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(u(i+1)-u(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(u(i)-u(i-1)))/(x(i+1)-x(i-1));
  D2uDx2(i) = (xhalf(i)*(u(i+1)-u(i))/(x(i+1)-x(i))-xhalf(i-1)*(u(i)-u(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
  
  DvDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(v(i+1)-v(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(v(i)-v(i-1)))/(x(i+1)-x(i-1));
  D2vDx2(i) = (xhalf(i)*(v(i+1)-v(i))/(x(i+1)-x(i))-xhalf(i-1)*(v(i)-v(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
  
  DcDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(c(i+1)-c(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(c(i)-c(i-1)))/(x(i+1)-x(i-1));
  D2cDx2(i) = (xhalf(i)*(c(i+1)-c(i))/(x(i+1)-x(i))-xhalf(i-1)*(c(i)-c(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
  
  DdDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(d(i+1)-d(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(d(i)-d(i-1)))/(x(i+1)-x(i-1));
  D2dDx2(i) = (xhalf(i)*(d(i+1)-d(i))/(x(i+1)-x(i))-xhalf(i-1)*(d(i)-d(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
end

% governing equations 
for i=1:n-1
  DuDt(i) = Dax*D2uDx2(i) - vel*DuDx(i) - k1*u(i)*v(i);
  DvDt(i) = Dax*D2vDx2(i) - vel*DvDx(i) - k1*u(i)*v(i)-k2*v(i)*c(i);
  DcDt(i) = Dax*D2cDx2(i) - vel*DcDx(i) + k1*u(i)*v(i) - k2*v(i)*c(i);
  DdDt(i) = Dax*D2dDx2(i) - vel*DdDx(i) + k2*v(i)*c(i);
end

% boundary condition at z = L
DuDt(n) = 0.0;
DvDt(n) = 0.0;
DcDt(n) = 0.0;
DdDt(n) = 0.0;

DyDt = [DuDt;DvDt;DcDt;DdDt];

end

% pdepe function definition
function [c,f,s] = pdex2pde(z,t,C,dCdz)

k1 = 1.00; % kinetic rate constant
k2 = 1.00; % kinetic rate constant
vz = 1.00; % superficial velocity
Dz = 1e-4; % dispersion coefficient
c = [1;1;1;1];
f = [Dz*dCdz(1);
Dz*dCdz(2);
Dz*dCdz(3);
Dz*dCdz(4)];
s = [-vz*dCdz(1)-k1*C(1)*C(2);
-vz*dCdz(2)-k1*C(1)*C(2)-k2*C(2)*C(3);
-vz*dCdz(3)+k1*C(1)*C(2)-k2*C(2)*C(3);
-vz*dCdz(4)+ k2*C(2)*C(3)];
end

function [pl,ql,pr,qr] = pdex2bc(zl,Cl,zr,Cr,t)
Cin =[1 1 0 0]; % inital concentrations
Dz =1e-4; % dispersion coefficient
pl = [Cl(1)-Cin(1);Cl(2)-Cin(2);Cl(3)-Cin(3);Cl(4)-Cin(4)];
ql = [-1;-1;-1;-1];
pr = [0;0;0;0];
qr = [1/Dz;1/Dz;1/Dz;1/Dz];
end

function C0 = pdex2ic(z);
C0=[1 ;1 ;0 ;0];
end

Answer 1

[Torsten]

For the pdepe version, better use

function [pl,ql,pr,qr] = pdex2bc(zl,Cl,zr,Cr,t)
  Cin =[1 1 0 0]; % inital concentrations
  vz = 1.00; % superficial velocity
  pl = [Cl(1)-Cin(1);Cl(2)-Cin(2);Cl(3)-Cin(3);Cl(4)-Cin(4)];
  ql = [-vz;-vz;-vz;-vz];
  pr = [0;0;0;0];
  qr = [1;1;1;1];
end

I adopted the method-of-lines version:

function main
L = 10; % length of reactor
x = linspace(0,L,100); % discretise the domain in 500 discretes
tspan = linspace(0,20,2000); % given time span
n = numel(x);
% initial conditions
u0 = ones(n,1); % concentration of A
v0 = ones(n,1); % concentration of B
c0 = zeros(n,1); % concentration of C
d0 = zeros(n,1); % concentration of D
y0 = [u0;v0;c0;d0];
% setting up the solver 
[T,Y] = ode15s(@(t,y) numericalprimertubular(t,y,x,n),tspan,y0);
Y=Y.';
% saving results
myCA = Y(1:n,:);
myCB = Y(n+1:2*n,:);
myCC = Y(2*n+1:3*n,:);
myCD = Y(3*n+1:4*n,:);
figure(1)
plot(x.',[myCA(:,100),myCA(:,500),myCA(:,1000),myCA(:,1500),myCA(:,2000)])
figure(2)
plot(x.',[myCB(:,100),myCB(:,500),myCB(:,1000),myCB(:,1500),myCB(:,2000)])
figure(3)
plot(x.',[myCC(:,100),myCC(:,500),myCC(:,1000),myCC(:,1500),myCC(:,2000)])
figure(4)
plot(x.',[myCD(:,100),myCD(:,500),myCD(:,1000),myCD(:,1500),myCD(:,2000)])
end
% function definition for method of lines solution
function DyDt = numericalprimertubular(t,y,x,n)

persistent iter
if isempty(iter)
  iter = 0;
end
iter = iter + 1;

vel = 1; % given velocity
Dax = 1e-4; % axial dispersion coefficient
k1 = 1; k2 = 1; % kinetic constants
u = zeros(n,1); % concentration of A
v = zeros(n,1); % concentration of B
c = zeros(n,1); % concentration of C
d = zeros(n,1); % concentration of D
DuDt = zeros(n,1);
DvDt = zeros(n,1);
DcDt = zeros(n,1);
DdDt = zeros(n,1);
DyDt = zeros(4*n,1);
xhalf = zeros(n-1,1);
DuDx = zeros(n-1,1);
DvDx = zeros(n-1,1);
DcDx = zeros(n-1,1);
DdDx = zeros(n-1,1);
D2uDx2 = zeros(n-1,1); 
D2vDx2 = zeros(n-1,1);
D2cDx2 = zeros(n-1,1);
D2dDx2 = zeros(n-1,1);
u = y(1:n);
v = y(n+1:2*n);
c = y(2*n+1:3*n);
d = y(3*n+1:4*n);
xhalf(1:n-1)=(x(1:n-1)+x(2:n))/2;
% danckwerts boundary condition at x=0
DuDx(1) = (u(1)-1)*(vel/Dax); 
DvDx(1) = (v(1)-1)*(vel/Dax);
DcDx(1) = (c(1)-0)*(vel/Dax);
DdDx(1) = (d(1)-0)*(vel/Dax);
% not sure about the applicability of these conditons !
D2uDx2(1) = ((u(2)-u(1))/(x(2)-x(1))-DuDx(1))/((x(2)-x(1))/2);
D2vDx2(1) = ((v(2)-v(1))/(x(2)-x(1))-DvDx(1))/((x(2)-x(1))/2);
D2cDx2(1) = ((c(2)-c(1))/(x(2)-x(1))-DcDx(1))/((x(2)-x(1))/2);
D2dDx2(1) = ((d(2)-d(1))/(x(2)-x(1))-DdDx(1))/((x(2)-x(1))/2);
%D2uDx2(1) = 1.0/(xhalf(1)-x(1))*(u(2)-u(1))/(x(2)-x(1)); 
%D2vDx2(1) = 1.0/(xhalf(1)-x(1))*(v(2)-v(1))/(x(2)-x(1)); 
%D2cDx2(1) = 1.0/(xhalf(1)-x(1))*(c(2)-c(1))/(x(2)-x(1));
%D2dDx2(1) = 1.0/(xhalf(1)-x(1))*(d(2)-d(1))/(x(2)-x(1));
for i=2:n-1
  DuDx(i) = (u(i)-u(i-1))/(x(i)-x(i-1));
  %DuDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(u(i+1)-u(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(u(i)-u(i-1)))/(x(i+1)-x(i-1));
  D2uDx2(i) = (xhalf(i)*(u(i+1)-u(i))/(x(i+1)-x(i))-xhalf(i-1)*(u(i)-u(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
  
  DvDx(i) = (v(i)-v(i-1))/(x(i)-x(i-1));
  %DvDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(v(i+1)-v(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(v(i)-v(i-1)))/(x(i+1)-x(i-1));
  D2vDx2(i) = (xhalf(i)*(v(i+1)-v(i))/(x(i+1)-x(i))-xhalf(i-1)*(v(i)-v(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
  
  DcDx(i) = (c(i)-c(i-1))/(x(i)-x(i-1));
  %DcDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(c(i+1)-c(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(c(i)-c(i-1)))/(x(i+1)-x(i-1));
  D2cDx2(i) = (xhalf(i)*(c(i+1)-c(i))/(x(i+1)-x(i))-xhalf(i-1)*(c(i)-c(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));

  DdDx(i) = (d(i)-d(i-1))/(x(i)-x(i-1));  
  %DdDx(i) = ((x(i)-x(i-1))/(x(i+1)-x(i))*(d(i+1)-d(i))+(x(i+1)-x(i))/(x(i)-x(i-1))*(d(i)-d(i-1)))/(x(i+1)-x(i-1));
  D2dDx2(i) = (xhalf(i)*(d(i+1)-d(i))/(x(i+1)-x(i))-xhalf(i-1)*(d(i)-d(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));
end
% boundary condition at z = L
DuDx(n) = 0.0;
DvDx(n) = 0.0;
DcDx(n) = 0.0;
DdDx(n) = 0.0;
D2uDx2(n) = -(u(n)-u(n-1))/(x(n)-x(n-1))/((x(n)-x(n-1))/2);
D2vDx2(n) = -(v(n)-v(n-1))/(x(n)-x(n-1))/((x(n)-x(n-1))/2);
D2cDx2(n) = -(c(n)-c(n-1))/(x(n)-x(n-1))/((x(n)-x(n-1))/2);
D2dDx2(n) = -(d(n)-d(n-1))/(x(n)-x(n-1))/((x(n)-x(n-1))/2);
% governing equations 
for i=1:n
%for i=1:n-1
  DuDt(i) = Dax*D2uDx2(i) - vel*DuDx(i) - k1*u(i)*v(i);
  DvDt(i) = Dax*D2vDx2(i) - vel*DvDx(i) - k1*u(i)*v(i) - k2*v(i)*c(i);
  DcDt(i) = Dax*D2cDx2(i) - vel*DcDx(i) + k1*u(i)*v(i) - k2*v(i)*c(i);
  DdDt(i) = Dax*D2dDx2(i) - vel*DdDx(i) + k2*v(i)*c(i);
end
% boundary condition at z = L
%DuDt(n) = 0.0;
%DvDt(n) = 0.0;
%DcDt(n) = 0.0;
%DdDt(n) = 0.0;
DyDt = [DuDt;DvDt;DcDt;DdDt];
if mod(iter,100) == 0
  iter = 0;  
  disp(t);
end  
end

[Saad Abdul Samad Posharkar]

Thank you so much Torsten for the quick response.

May I know how you adopted the boundary conditions, particularly the second derivative at the nodes, as this impacted the solution too, apart from the discresation scheme. I would appreciate if you could share any literature which states/outlines the approach for a particular problem.

[Torsten]

Discretization schemes are taken from

https://opus4.kobv.de/opus4-zib/frontdoor/index/index/year/1995/docId/191

If you use the activated discretization for DuDx

DuDx(i) = (u(i)-u(i-1))/(x(i)-x(i-1)) 

you should use more grid points since this method increases numerical diffusion. The advantage is that it is stable and avoids oscillations of the solution. This is relevant for flows with high Peclet numbers:

https://www.researchgate.net/publication/325531692_Mathematical_relationship_between_grid_and_low_Peclet_numbers_for_the_solution_of_convection-diffusion_equation/link/5b12b442a6fdcc4611dea495/download

[Saad Abdul Samad Posharkar]

Thank you once again, for your detailed and patient reply. This helps a lot!

[Torsten]

Since I have not the chance to run pdepe, maybe you could report after testing whether the solutions from both solvers now match better.

[Saad Abdul Samad Posharkar]

The solution from both the approaches match much better now, close or similar values.

Also, adapting the boundary conditions for the pdepe solver did not change the solution. The values are exactly the same.

[Torsten]

Also, adapting the boundary conditions for the pdepe solver did not change the solution. The values are exactly the same.

By chance, vz = 1 in your case. But if "vz" changes, this should make a difference.

[Saad Abdul Samad Posharkar]

Yes, I did notice this. However, I didn't understand the omission of the dispersion coefficient from the boundary condition.

I was playing around with the pdepe boundary function and I noticed no change in the solution for the following (p+qf = 0 format), as well. Is this because of the geometry of the problem? (I am not very much familiar with the pdepe solver and will read on its documentation)

function [pl,ql,pr,qr] = pdex2bc(zl,Cl,zr,Cr,t)
Cin =[1 1 0 0];
Dz =1e-4;
pl = [Cl(1)-Cin(1);Cl(2)-Cin(2);Cl(3)-Cin(3);Cl(4)-Cin(4)];
ql = [-1;-1;-1;-1];
pr = [0;0;0;0];
qr = [Dz;Dz;Dz;Dz];
end

[Torsten]

The equation to set the boundary conditions is

p + q*f = 0

You set

f = Dz*du/dx

So in order to implement

du/dx = vz/Dz * (u - u_in)

or

- vz*(u-u_in) + Dz * du/dx = 0

at the left boundary point, you have to determine p and q such that

p + q * Dz*du/dx = -vz * (u - u_in) + Dz * du/dx

thus

ql = 1 and  pl = -vz * (u-u_in)

Thus a "vz" is missing in your settings for the boundary condition to the left.

In order to implement

Dz * du/dx = 0

at the right boundary point, you can write

p + q*f = 0 + 1* Dz*du/dz,

thus

pr=0, qr=1

.

[Saad Abdul Samad Posharkar]

Thanks a ton! After reading the pdepe documentation, I got some more clarity as well. These conversations have been very helpful and informative!

[Torsten]

For consistency, you should change

DuDx = zeros(n-1,1);
DvDx = zeros(n-1,1);
DcDx = zeros(n-1,1);
DdDx = zeros(n-1,1);
D2uDx2 = zeros(n-1,1); 
D2vDx2 = zeros(n-1,1);
D2cDx2 = zeros(n-1,1);
D2dDx2 = zeros(n-1,1);

to

DuDx = zeros(n,1);
DvDx = zeros(n,1);
DcDx = zeros(n,1);
DdDx = zeros(n,1);
D2uDx2 = zeros(n,1); 
D2vDx2 = zeros(n,1);
D2cDx2 = zeros(n,1);
D2dDx2 = zeros(n,1);

[Saad Abdul Samad Posharkar]

Surely, I'll incorporate this. Thank you so so much.

[Torsten]

Please note that the approximation of second derivatives in the code above is wrong.

You will have to replace

D2uDx2(i) = (xhalf(i)*(u(i+1)-u(i))/(x(i+1)-x(i))-xhalf(i-1)*(u(i)-u(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));

by

D2uDx2(i) = ((u(i+1)-u(i))/(x(i+1)-x(i))-(u(i)-u(i-1))/(x(i)-x(i-1)))/(xhalf(i)-xhalf(i-1));

(Same for D2vDx2(i),D2cDx2(i) and D2dDx2(i))

[Saad Abdul Samad Posharkar]

Thank you Torsten for this update. As mentioned in the question, I was confused about this bit but as I got a relatively good convergence, I did not dive much deeper into the implementation of xhalf at the individual nodes.

Making the changes, I did not observe any change in the values/speed of computation.

[Torsten]

That's surprising - the result of the computation should change since the xhalf - terms are deleted.

But maybe your diffusion coefficients are that small that you can see no effect in the concentrations.

[Saad Abdul Samad Posharkar]

Quiet possible as there is orders of magnitude difference b/w the convection and diffusion term. The convection dominating, thus a plug type behaviour along the length.

[Torsten]

That was a close shave ! :-)