Solving coupled Partial Differential Equations using method of Lines

64 ビュー (過去 30 日間)
Saad Abdul Samad Posharkar
Saad Abdul Samad Posharkar 2022 年 1 月 13 日
コメント済み: Torsten 2022 年 3 月 15 日
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

採用された回答

Torsten
Torsten 2022 年 1 月 13 日
編集済み: Torsten 2022 年 1 月 13 日
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
  16 件のコメント
Torsten
Torsten 2022 年 3 月 15 日
That was a close shave ! :-)

サインインしてコメントする。

その他の回答 (0 件)

Community Treasure Hunt

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

Start Hunting!

Translated by