PDE model not producing desired results

6 vues (au cours des 30 derniers jours)
Garion Charles
Garion Charles le 15 Avr 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

Réponse acceptée

Torsten
Torsten le 15 Avr 2022
Modifié(e) : Torsten le 15 Avr 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 commentaires
Garion Charles
Garion Charles le 24 Avr 2022
Modifié(e) : Garion Charles le 24 Avr 2022
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
Adam Ali le 20 Mai 2022
Hello Garion, were you able to rectify the issue?

Connectez-vous pour commenter.

Plus de réponses (1)

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

Produits


Version

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by