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

25 views (last 30 days)
Mohammed Hamza
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 ;
DFaDv(i) = (Fa(i)-Fa(i-1))/(v(i)-v(i-1));
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
DCaDt(i) = -DFaDv(i) + rA(i) ;
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
DyDt = [DCaDt;DCbDt;DCiDt;DTDt;DXDv];
end

Accepted Answer

Torsten
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 ;
DFaDv(2:n) = (Fa(2:n)-Fa(1:n-1))./(v(2:n)-v(1:n-1));
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
DCaDt(2:n) = -DFaDv(2:n) + rA(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
DyDt = [DCaDt;DCbDt;DCiDt;DTDt;DXDv];
end
  12 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by