PDE model not producing desired results

Currently trying to replicate the results from a research article on a model for supercritical fluid extraction using CO2, which takes the form of a system of Coupled PDEs. However i have not had any sucess thus far as the results that i have been obtaining are nowhere close to the values they should be. If anyone can provide any further assistance it would be much appreciated.

I have also attach the full article just in case i missout any details

The PDE model:

The Parameters (Currently investigating Rose for a flowrate of 1.08 Kg h-1)

dp = 2mm

delta/delta(o) = 0.01

How the result should look (trying to replicate the results for y vs z and x vs z at 50 min)

Function file

function dsdt = SC_CO2(t,s,n,dz,Ap,yo,ho,xt,u,rouf,rous,e,ec,ys1)

y = s(1:n);
x = s(n+1:2*n);

dydt = zeros(n,1); % setting the size of the array for dydt 
dxdt = zeros(n,1); % setting the size of the array for dxdt

 for i = 2:n % evaluation dxdt for each node from node 1 to n    
     
     h = ho/(((xt - x(i))/xt)+ 0.01);
     J = Ap*h*rouf*(yo - y(i));
     K1 = J/(rouf*e);
     
     dydt(i) = -u*((y(i) - y(i-1))/(dz)) + K1 ; % dydz  at nodes 2 to n-1 is solved using backwards difference approx
          
 end
    
 for i = 1:n % evaluation dxdt for each node from node 1 to n
      
     h = ho/(((xt - x(i))/xt)+0.01);
     J = Ap*h*rouf*(yo - y(i));
     K2 = -J/(rous*ec);
     
     dxdt(i) = K2;
        
 end
    
 %--------------------------------Boundary condition---------------------------------------------
 
 h = ho/(((xt - x(1))/xt)+ 0.01);
 J = Ap*h*rouf*(yo - y(1));
 K1 = J/(rouf*e);
 
 dydt(1) = -u*((y(1) - ys1)/(dz)) + K1 ; % dydz at node 1 is evaluated using backward difference approx where y(0) = ys1 = 0 is the boundary y value at z = 0
 
 %-------------------------------------------------------------------------------------------------------------------
 
 
 dsdt = [dydt ; dxdt];





end

Script File

clear
clc

u = 0.850E-3; % Velocity ms-1
rouf = 280;   % Fluid density kg m-3
rous = 612;   % Bulk dnesity of the not soluble part of concrete kg m-3
e = 0.3;      % Fluid volume fraction (dimensionless)
ec = 0.1;     % Concrete volume fraction (dimensionless)
xt = 0.33;    % Solute fraction in the concrete at transition between extraction regimes (dimensionless)
dp = 0.002;   % Glass particle diameter m
ho = 8.6E-7;  % Proportionality constant ms-1
yo = 0.0014;  % Initial solute fraction in the solvent (kg solute/kg solvent) 
xo = 0.96;    % Initial solute fraction in the concrete(kg solute/kg not soluable concrete)
ys1 = 0;      % Values of y at boundary node i.e y(0,t) = 0



%---------------------------Calculation of constant parameters for model-----------------------------

kp = yo/xt;         % Partition constant
Ap = (6*(1-e))/dp;  % Specific surface of mass transfer area m-1
ye = kp*xo;         % Solute fraction in the solvent in equilibrium with the concrete 

%----------------------------- z dimension nodal declaration --------------------------

Z = 0.02;     % space coordinate m
n = 101;      % No. of nodes
dz = Z/(n-1); % distance between nodes m
z = 0 :dz: Z  ;  
z = z.';

%------------------------------------Initial conditions-------------------------------------

yi = ones(n,1)*ye; % i.e y(z,0) = yo
xi = ones(n,1)*xo; % i.e x(z,0) = xo 

ystart = [yi ; xi];

%-----------------------------------Model Run Time-----------------------------------

t_final = 3000;
dt = 1;
tspan = 0 : dt : t_final;

%------------------------------Solving the Resulting ODE system-----------------------

[T,Y] = ode15s(@(t,s) SC_CO2(t,s,n,dz,Ap,yo,ho,xt,u,rouf,rous,e,ec,ys1), tspan , ystart);

plot(z,Y(3001,1:n)) % ploting  y vs z at 3000s = 50 min

Answer 1

[Torsten]

You need to impose a boundary condition for y at z=0.

Since the equations are hyperbolic in nature, central differencing is strictly forbidden. You'll have to approximate 1st order derivatives in space by an upwind scheme.

Why don't you stick to the method of characteristics that the author suggests to use ? It is the method-of-choice in the given context.

[Garion Charles]

Hey thanks for the advice I finally manage to figure out how to implement what you suggested and it help in obtaining the general shape but am now getting large relative errors between my graph and the top graph in my original question, could it have somthing to do with how i implemented my discretization scheme ?

I have tried increasing the number of nodes and also reducing the time step size and switching to other ode solvers but this is about the highest accuracy i am able to obtain

Function file

function dsdt = SC_CO2(t,s,n,dz,Ap,kp,yo,ho,xt,u,rouf,rous,e,ec,ys1)

y = s(1:n);
x = s(n+1:2*n);

dydt = zeros(n,1); % setting the size of the array for dydt 
dxdt = zeros(n,1); % setting the size of the array for dxdt

 for i = 2:n % evaluation dydt for each node from node 2 to n  
             
     
     if x(i) < xt
         ye = kp*x(i);  %kp*x(i)
         h = ho/(((xt - x(i))/xt)+ 0.01);
         J = Ap*h*rouf*(ye - y(i));
         K1 = J/(rouf*e);
         
         dydt(i) = -u*((y(i) - y(i-1))/(dz)) + K1 ; % dydz  at nodes 2 to n-1 is solved using backwards difference approx 
         
     else
         
         ye = yo;
         h = 2*ho/u;
         J = Ap*h*rouf*(ye - y(i));
         K1 = J/(rouf*e);
         
         dydt(i) = -u*((y(i) - y(i-1))/(dz)) + K1 ;  
     end
     
         
          
 end
    
 for i = 1:n % evaluation dxdt for each node from node 1 to n
     
     if x(i) < xt
         ye = kp*x(i);
         h = ho/(((xt - x(i))/xt)+ 0.01);
         J = Ap*h*rouf*(ye - y(i));
         K2 = -J/(rous*ec);
         
         dxdt(i) = K2;
         
     else
         
         ye = yo;
         h =  2*ho/u;
         J = Ap*h*rouf*(ye - y(i));
         K2 = -J/(rous*ec);
         
         dxdt(i) = K2;
         
     end
     
              
 end
    
 %--------------------------------Boundary condition---------------------------------------------
 
 if x(1) < xt
     
     ye = kp*x(1);
     h = ho/(((xt - x(1))/xt)+ 0.01);
     J = Ap*h*rouf*(ye - y(1));
     K1 = J/(rouf*e);
     dydt(1) = -u*((y(1) - ys1)/(dz)) + K1; % dydz at node 1 is evaluated using backward difference approx where y(0) = ys1 = 0 is the boundary y value at z = 0
 else
     
    ye = yo;
    h =  2*ho/u;
    J = Ap*h*rouf*(ye - y(1));
    K1 = J/(rouf*e); 
    dydt(1) = -u*((y(1) - ys1)/(dz)) + K1; % dydz at node 1 is evaluated using backward difference approx where y(0) = ys1 = 0 is the boundary y value at z = 0
 end
 
 
 
 
 %-------------------------------------------------------------------------------------------------------------------
 
 
 dsdt = [dydt ; dxdt];





end

Script file

% Main script file for the determination of y and x along z axis at time t

clear
clc

u = 0.850E-3; % Velocity ms-1
rouf = 280;   % Fluid density kg m-3
rous = 612;   % Bulk dnesity of the not soluble part of concrete kg m-3
e = 0.3;      % Fluid volume fraction (dimensionless)
ec = 0.1;     % Concrete volume fraction (dimensionless)
xt = 0.33;    % Solute fraction in the concrete at transition between extraction regimes (dimensionless)
dp = 0.002;   % Glass particle diameter m
ho = 8.6E-7;  % Proportionality constant ms-1
yo = 0.0014;  % Initial solute fraction in the solvent (kg solute/kg solvent) 
xo = 0.96;    % Initial solute fraction in the concrete(kg solute/kg not soluable concrete)
ys1 = 0;      % Values of y at boundary node i.e y(0,t) = 0


%---------------------------Calculation of constant parameters for model-----------------------------

kp = yo/xt;         % Partition constant
Ap = (6*(1-e))/dp;  % Specific surface of mass transfer area m-1
        

%----------------------------- z dimension nodal declaration --------------------------

Z = 0.16;     % space coordinate m
n = 501;      % No. of nodes
dz = Z/(n-1); % distance between nodes m
z = 0 :dz: Z  ;  


%------------------------------------Initial conditions-------------------------------------

yi = ones(n,1)*yo; % i.e y(z,0) = yo
xi = ones(n,1)*xo; % i.e x(z,0) = xo 

ystart = [yi ; xi];

%-----------------------------------Model Run Time-----------------------------------

t_final = 3000;
dt = 1;
tspan = 0 : dt : t_final;

%-----------------------------Solving the Resulting ODE system-----------------------

[T,Y] = ode45(@(t,s) SC_CO2(t,s,n,dz,Ap,kp,yo,ho,xt,u,rouf,rous,e,ec,ys1), tspan , ystart);

plot(z,Y(3001,1:n)) % ploting  y vs z at 3000s = 50 min

[Adam Ali]

Hello Garion, were you able to rectify the issue?

Answer 2

[Adam Ali]

Hello Garion, were you able to rectify the issue with the graphs?