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

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

Answer 1

[Torsten]

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

[Mohammed Hamza]

Yeah, Ca becomes negative is one thing. The other thing is the temperature profile is not as the solution that I should get ( Have a look on the gif that I have attached).

[Torsten]

If the attached gif are the COMSOL results, you must have used different equations therein.

The reaction rate rA should approach 0 as Ca approaches 0. This does not seem to be the case in the above code. And the reaction rate determines everything - also the temperature profile.

[Mohammed Hamza]

I don't think so, in COMSOL I have opend a general mathmatical model and applied the same equations. To further check & to make sure about my equations implementation, I'll attach all equations used in order to solve the COMSOL model.

[Torsten]

Which boundary conditions did you impose in COMSOL ?

I've worked with this software quite a long time and I know that one has to specify boundary conditions at both ends (l=0 and l=L) which is not appropriate for your kind of hyperbolic equations from above.

One error in the above code is that you try to solve

dX/dt = -rA/(163*0.9*0.1)

instead of

dX/dv = -rA/(163*0.9*0.1)

[Mohammed Hamza]

As mentioned in my code, I have used Dirichlet boundary condition which required me to insert just the prescribed value for each state.

[Torsten]

Can't be true at l=L.

Or - if you prescribed Dirichlet conditions at l=L - you'll get a nonphysical solution from COMSOL.

[Mohammed Hamza]

I couldn't see the error in my code for this:

%%%%%%%%%%%%%%%%%%%%%%%%%

One error in the above code is that you try to solve

dX/dt = -rA/(163*0.9*0.1)

instead of

dX/dv = -rA/(163*0.9*0.1)

%%%%%%%%%%%%%%%%%%%%%%%%%

[Torsten]

The vector DyDt you return to ODE15S is the vector of time derivatives of the variables.

So although you name the last part DXDv, integration will be with respect to t, not with respect to v.

[Mohammed Hamza]

Oh, sorry I have just understood what you have meant.

[Mohammed Hamza]

Sorry Torsten for taking your time, I have tried to fix it in this way but it seems I made it worse could you guide me into the right way of doing it?

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


y0 = [Ca0;Cb0;Ci0;T0;];

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

     0

     0

     0

     0

    2.2847

    7.6275

    7.6275

    7.6275

    7.6275

   13.3706

   13.3706

   13.3706

   13.3706

   19.4882

   19.4882

   19.4882

   19.4882

   83.9923

   83.9923

   83.9923

   83.9923

   89.7179

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,:);
 

% filename = 'EX8-4_Temp.gif'; % file to store image in
% for i=1:2:length(t) % only look at every 10th time step
%     plot(v,T(:,i));
% %     ylim([310 330])
%     xlim([0   5])
%     xlabel('Rx volume')
%     ylabel('Temperature')
%     title('time = ',t(i))
%     text(10,10,sprintf('t =',t(i))); % add text annotation of time step
%     frame = getframe(1);
%     im = frame2im(frame); %convert frame to an image
%     [imind,cm] = rgb2ind(im,256);
% 
%     % now we write out the image to the animated gif.
%     if i == 1
%         imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
%     else
%         imwrite(imind,cm,filename,'gif','WriteMode','append');
%     end
% end
% close all; % makes sure the plot from the loop above doesn't get published


function DyDt = solution(t,y,v,n)

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

    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
    Fa0    = 163*0.9*0.1;
    Fb0    = 0          ;
    Fi0    = 163*0.1*0.1;


    DCaDt = zeros(n,1); % Change in concentration of A
    DCbDt = zeros(n,1); % Change in concentration of B
    DCiDt = zeros(n,1); % Change in concentration of inert
    DTDt  = zeros(n,1); % Change in Temperature

    DFaDv = zeros(n,1);
    DFbDv = zeros(n,1);
    DFiDv = zeros(n,1);
    DTDv  = zeros(n,1);
    X     = 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

    IC = zeros(n,1); % Initial values for the dependent variables T, and x
    Vspan = [0 5]; % Range for the independent variable i.e.volume of the reactor
    [V,sv] = ode15s(@(V,sv) ConvRate(V,sv,n,Fa0,T), Vspan, IC);
    X = sv(end,:);


    Fa    = zeros(n,1);
    Fb    = zeros(n,1);
    Fi    = zeros(n,1);

    Fa(1)    = Fa0;
    Fb(1)    = Fb0;
    Fi(1)    = Fi0;

    k     = zeros(n,1);
    Kc    = zeros(n,1);
    rA    = zeros(n,1);
    
    Hc    = zeros(n,1);
    HcPv  = zeros(n,1);

 
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));  
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));
end


 
DyDt = [DCaDt;DCbDt;DCiDt;DTDt;];
if mod(iter,100) == 0
  iter = 0;  
  disp(t);
end  
end





function output = ConvRate(V,s,n,Fa0,T)
    X = s;
    k     = zeros(n,1);
    Kc    = zeros(n,1);
    rA    = zeros(n,1);
    dXdV  = zeros(n,1);

    
    for i=2:n

    % Explicit equations
    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));
    
    
    % Differential equations
    dXdV (i) = 0 - (rA(i) / Fa0);
    end


    output = dXdV;  

end

[Torsten]

IC = 0; % Initial values for the dependent variables T, and x
Vspan = v; % Range for the independent variable i.e.volume of the reactor
[~,sv] = ode15s(@(V,sv) ConvRate(V,sv,Fa0,v,T), Vspan, IC);
X = sv(:,1);

function dXdV = ConvRate(V,X,Fa0,v,T)
  T_actual = interp1(v,T,V);
  k  = 31.1*exp(7906  *(T_actual - 360)/(360*T_actual));
  Kc  = 3.03*exp(-830.3*(T_actual - 333)/(333*T_actual));
  rA = -k*9.3*(1 - (1 + 1/Kc)*X);
  dXdV = - rA / Fa0; 
end

[Mohammed Hamza]

Thank you Torsten, you saved me weeks of work, I really appritiate it.