# PDE model not producing desired results

1 view (last 30 days)
Garion Charles on 15 Apr 2022
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

Torsten on 15 Apr 2022
Edited: Torsten on 15 Apr 2022
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.
##### 2 CommentsShowHide 1 older comment
Adam Ali on 20 May 2022
Hello Garion, were you able to rectify the issue?

Adam Ali on 20 May 2022
Hello Garion, were you able to rectify the issue with the graphs?

R2018b

### Community Treasure Hunt

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

Start Hunting!

Translated by