Solving a System Consisting of PDEs & ODEs using Method of lines

25 views (last 30 days)
Mohammed Hamza on 28 Jun 2022
Commented: Mohammed Hamza on 29 Jun 2022
Hello,
I have recently tried to solve a mixed system of PDEs & ODEs, thanks to Torsten, he helped me somehow to get the code working. However, I came to a problem which started to give wierd solution using the same exact method and lines of code, so I have started an investegation, I have solved the same system using COMSOL Multiphysics just for the sake of experimenting, the results from COMSOL was spot on (a gif is attached). So I have Concluded that the code I haved used contained an issue that the I haven't been able to find through the past 2 weeks. I would be so thankful if anyone could help me with this problem.
This is the code that I have used:
L = 5; % length of reactor
v = linspace(0,L,100); % discretise the domain in 50 discretes
tspan = linspace(0,100,200); % given time span
n = numel(v);
% Dirichlet boundary condition at x=0
Ca0 = 9.3 *ones(n,1); % concentration of A
Cb0 = 0 *ones(n,1); % concentration of B
Ci0 = 1.033 *ones(n,1); % concentration of inert
T0 = 310 *ones(n,1); % Temperature of the reactor
X0 = 0 *ones(n,1); % Conversion rate
y0 = [Ca0;Cb0;Ci0;T0;X0];
% setting up the solver
[T,Y] = ode15s(@(t,y) solution(t,y,v,n),tspan,y0);
Y=Y.';
% saving results
Ca = Y( 1: n,:);
Cb = Y( n+1:2*n,:);
Ci = Y(2*n+1:3*n,:);
T = Y(3*n+1:4*n,:);
X = Y(4*n+1:5*n,:);
function DyDt = solution(t,y,v,n)
Cpa = 141 ; % Heat Capacity of Material A
Cpb = 141 ; % Heat Capacity of Material B
Cpi = 161 ; % Heat Capacity of inert Material
delHrxn = -6900; % Heat of Reaction
Ua = 5000 ; % Heat Transfare Coeffient
Ta = 310; % Ambiant Temperature
vo = 1.5774 ; % velocity of the flow
DCaDt = zeros(n,1); % Change in concentration of A with time
DCbDt = zeros(n,1); % Change in concentration of B with time
DCiDt = zeros(n,1); % Change in concentration of inert with time
DTDt = zeros(n,1); % Change in Temperature with time
%% Note that : DFiDv = DCiDv * vo
DFaDv = zeros(n,1); % Change in Flow of A with volume
DFbDv = zeros(n,1); % Change in Flow of A with volume
DFiDv = zeros(n,1); % Change in Flow of A with volume
DTDv = zeros(n,1); % Change in Temperature with volume
DXDv = zeros(n,1); % Change in Conversion rate with volume
Fa = zeros(n,1);
Fb = zeros(n,1);
Fi = zeros(n,1);
Fa(1) = 163*0.9*0.1;
Fb(1) = 0 ;
Fi(1) = 163*0.1*0.1;
k = zeros(n,1);
Kc = zeros(n,1);
rA = zeros(n,1);
Hc = zeros(n,1);
HcPv = zeros(n,1);
Ca = y(1:n); % concentration of A
Cb = y(n+1:2*n); % concentration of B
Ci = y(2*n+1:3*n); % concentration of inert
T = y(3*n+1:4*n); % Temperature
X = y(4*n+1:5*n); % Conversion rate
for i=2:n
Fa (i) = Ca(i)*vo ;
Fb (i) = Cb(i)*vo ;
Fi (i) = Ci(i)*vo ;
k (i) = 31.1*exp(7906 *(T(i) - 360)/(360*T(i)));
Kc (i) = 3.03*exp(-830.3*(T(i) - 333)/(333*T(i)));
rA (i) = -k(i)*9.3*(1 - (1 + 1/Kc(i))*X(i));
Hc (i) = Fa(i)*Cpa + Fb(i)*Cpb + Fi(i)*Cpi ;
HcPv (i) = Ca(i)*Cpa + Cb(i)*Cpb + Ci(i)*Cpi ;
DFbDv(i) = (Fb(i)-Fb(i-1))/(v(i)-v(i-1));
DFiDv(i) = (Fi(i)-Fi(i-1))/(v(i)-v(i-1));
DTDv(i) = (T(i)-T(i-1))/(v(i)-v(i-1));
DXDv(i) = (X(i)-X(i-1))/(v(i)-v(i-1));
end
% governing equations
for i=2:n
DCbDt(i) = -DFbDv(i) - rA(i) ;
DCiDt(i) = -DFiDv(i) ;
DTDt (i) = (Ua * (Ta-T(i)) - Hc(i) * DTDv(i) + (-rA(i))*(-delHrxn))/(HcPv(i));
DXDv (i) = -rA(i)/(163*0.9*0.1);
end
end

Torsten on 28 Jun 2022
What's wrong with the results ? The fact that Ca becomes negative ?
L = 5; % length of reactor
v = linspace(0,L,100).'; % discretise the domain in 50 discretes
tspan = linspace(0,100,200); % given time span
n = numel(v);
% Dirichlet boundary condition at x=0
Ca0 = 9.3 *ones(n,1); % concentration of A
Cb0 = 0 *ones(n,1); % concentration of B
Ci0 = 1.033 *ones(n,1); % concentration of inert
T0 = 310 *ones(n,1); % Temperature of the reactor
X0 = 0 *ones(n,1); % Conversion rate
y0 = [Ca0;Cb0;Ci0;T0;X0];
% setting up the solver
[T,Y] = ode15s(@(t,y) solution(t,y,v,n),tspan,y0);
Y=Y.';
% saving results
Ca = Y( 1: n,:);
Cb = Y( n+1:2*n,:);
Ci = Y(2*n+1:3*n,:);
T = Y(3*n+1:4*n,:);
X = Y(4*n+1:5*n,:);
figure(1)
plot(v,Ca)
figure(2)
plot(v,Cb)
figure(3)
plot(v,Ci)
figure(4)
plot(v,T)
figure(5)
plot(v,X)
function DyDt = solution(t,y,v,n)
Cpa = 141 ; % Heat Capacity of Material A
Cpb = 141 ; % Heat Capacity of Material B
Cpi = 161 ; % Heat Capacity of inert Material
delHrxn = -6900; % Heat of Reaction
Ua = 5000 ; % Heat Transfare Coeffient
Ta = 310; % Ambiant Temperature
vo = 1.5774 ; % velocity of the flow
DCaDt = zeros(n,1); % Change in concentration of A with time
DCbDt = zeros(n,1); % Change in concentration of B with time
DCiDt = zeros(n,1); % Change in concentration of inert with time
DTDt = zeros(n,1); % Change in Temperature with time
%% Note that : DFiDv = DCiDv * vo
DFaDv = zeros(n,1); % Change in Flow of A with volume
DFbDv = zeros(n,1); % Change in Flow of A with volume
DFiDv = zeros(n,1); % Change in Flow of A with volume
DTDv = zeros(n,1); % Change in Temperature with volume
DXDv = zeros(n,1); % Change in Conversion rate with volume
Fa = zeros(n,1);
Fb = zeros(n,1);
Fi = zeros(n,1);
Fa(1) = 163*0.9*0.1;
Fb(1) = 0 ;
Fi(1) = 163*0.1*0.1;
k = zeros(n,1);
Kc = zeros(n,1);
rA = zeros(n,1);
Hc = zeros(n,1);
HcPv = zeros(n,1);
Ca = y(1:n); % concentration of A
Cb = y(n+1:2*n); % concentration of B
Ci = y(2*n+1:3*n); % concentration of inert
T = y(3*n+1:4*n); % Temperature
X = y(4*n+1:5*n); % Conversion rate
%for i=2:n
Fa (2:n) = Ca(2:n)*vo ;
Fb (2:n) = Cb(2:n)*vo ;
Fi (2:n) = Ci(2:n)*vo ;
k (2:n) = 31.1*exp(7906 *(T(2:n) - 360)./(360*T(2:n)));
Kc (2:n) = 3.03*exp(-830.3*(T(2:n) - 333)./(333*T(2:n)));
rA (2:n) = -k(2:n)*9.3.*(1 - (1 + 1./Kc(2:n)).*X(2:n));
Hc (2:n) = Fa(2:n)*Cpa + Fb(2:n)*Cpb + Fi(2:n)*Cpi ;
HcPv (2:n) = Ca(2:n)*Cpa + Cb(2:n)*Cpb + Ci(2:n)*Cpi ;
DFbDv(2:n) = (Fb(2:n)-Fb(1:n-1))./(v(2:n)-v(1:n-1));
DFiDv(2:n) = (Fi(2:n)-Fi(1:n-1))./(v(2:n)-v(1:n-1));
DTDv(2:n) = (T(2:n)-T(1:n-1))./(v(2:n)-v(1:n-1));
DXDv(2:n) = (X(2:n)-X(1:n-1))./(v(2:n)-v(1:n-1));
%end
% governing equations
%for i=2:n
DCbDt(2:n) = -DFbDv(2:n) - rA(2:n) ;
DCiDt(2:n) = -DFiDv(2:n) ;
DTDt (2:n) = (Ua * (Ta-T(2:n)) - Hc(2:n) .* DTDv(2:n) + (-rA(2:n)).*(-delHrxn))./(HcPv(2:n));
DXDv (2:n) = -rA(2:n)/(163*0.9*0.1);
%end
end
Mohammed Hamza on 29 Jun 2022
Thank you Torsten, you saved me weeks of work, I really appritiate it.

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by