PDE solution not changing with time
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi,
As the title suggests. The code works without any errors or warnings. Furthermore, the values do go down as expected along the axial direction (horizontally in CFinal) but then they rise abnormally. Is there anyway I can closely monitor the abrupt rising of these values as well as the calculations done after the first time step to see why the values don't change?
I think that it's because of the equations. Specifically line 122 which is the boundary condition for diffusion in particle and the sole link to the fluid phase mass balance equation (line 127). I suspect the dCpdt(:,end) isn't exactly being shared with the two equations very well.
% This version uses the equations from 10.1186/s13068-016-0602-2
% This version shrinks the mass balance in the particle by using effective
% diffusion coefficient given by De = Dp + diff(Langmuir)*rho_b*Ds
% Data ====================================================================
eta = 0.78;
Ds = 0.001;
Dz = 9.89019E-05;
De = 0.001; %unused
u = 0.036;
C0 = 0.0437; %mg/L
rho = 997.3;
KL = 0.63; %1/mg
qmax = 154; %mg/g
kf = 0.027419644;
rho_b = 743;
% Axial and radial grid points ============================================
L = 0.04; %m
AxialGridPoints = 201;
IntAxGridPoints = AxialGridPoints - 1;
dz = L/AxialGridPoints;
rp = 1.1/1000; %m
PartGridPoints = 21;
IntPartGridPoints = PartGridPoints - 1;
dr = rp/PartGridPoints;
Radius = linspace(0,rp,PartGridPoints);
% Time ====================================================================
tspan = linspace(0,1440,100);
% Initial Conditions ======================================================
% reactor
C = zeros(AxialGridPoints,1);
C(1) = C0;
% Particle
Cp = zeros(AxialGridPoints, PartGridPoints);
% Packing data ============================================================
Initial = [C Cp];
Data = [eta,Ds,Dz,De,u,C0,rho,KL,qmax,kf,rho_b];
Spatial = [AxialGridPoints,IntAxGridPoints,PartGridPoints,...
IntPartGridPoints,Radius,dz,dr,rp];
% ODE solver ==============================================================
% options = odeset('OutputFcn', @odeOutputFcn);
[t, y] = ode15s(@(t, y) PoreDiffusion(t, y, Data, Spatial), tspan, Initial);
beep
% Result ==================================================================
CFinal = y(:,1:AxialGridPoints);
qFinal = y(:,AxialGridPoints+1:end);
plot(linspace(0,L,AxialGridPoints),[CFinal(2,:)])
% Loop ====================================================================
function dydt = PoreDiffusion(t,y,Data,Spatial)
%Unpacking data----------------------------------------------------
eta = Data(1);
Ds = Data(2);
Dz = Data(3);
De = Data(4);
u = Data(5);
C0 = Data(6);
rho = Data(7);
KL = Data(8);
qmax = Data(9);
kf = Data(10);
rho_b = Data(11);
AxialGridPoints = Spatial(1);
IntAxGridPoints = Spatial(2);
PartGridPoints = Spatial(3);
IntPartGridPoints = Spatial(4);
Radius = Spatial(5:length(Spatial)-3);
dz = Spatial(length(Spatial)-2);
dr = Spatial(length(Spatial)-1);
rp = Spatial(length(Spatial));
%------------------------------------------------------------------
% Sorting out input -----------------------------------------------
C = y(1:AxialGridPoints);
Cp = y(AxialGridPoints+1:end);
Cp = reshape(Cp,AxialGridPoints,PartGridPoints); %Rearranging into a grid
M = (3*(1-eta)*kf)/(rp*eta);
% coefficient = (2*kf*dr)/(Dp);
%Preallocating space
dCpdt = zeros(AxialGridPoints,PartGridPoints);
dCdt = zeros(AxialGridPoints,1);
% -----------------------------------------------------------------
for z = 2:IntAxGridPoints % z is the axial coordinate
% Particle
for r = 2:IntPartGridPoints
x = eta + (rho_b*qmax*KL)/(1 + KL*Cp(z,r))^2;
dCpdr = Cp(z,r)/2/dr - Cp(z,r)/2/dr;
d2Cpdr2 = (Cp(z,r+1) - 2*Cp(z,r) + Cp(z,r-1))/dr/dr;
L1 = (2*qmax*KL*KL*Ds*rho_b)/(1+KL*Cp(z,r))^3;
dCpdt(z,r) = De/x*d2Cpdr2 - L1/x*dCpdr^2 + (2*De)/(x*Radius(r))*dCpdr;
end
% Reactor
L2 = (3*De)/(2*dr) + kf;
dCpdt(:,end) = C(z)*kf/L2 + Cp(z,IntPartGridPoints)*(2*De/dr/L2) - Cp(z,IntPartGridPoints-1)*(De/2/dr/L2);
dCpdt(:,1) = 4*Cp(:,2)/3 - Cp(:,3)/3;
dCdz = C(z+1)/2/dz - C(z-1)/2/dz;
d2Cdz2 = (C(z+1) - 2*C(z) + C(z-1))/dz/dz;
dCdt(z) = -(u)*dCdz + (Dz)*d2Cdz2 - M.*(C(z) - dCpdt(z,end));
end
% Reassigning boundary conditions
dCdt(AxialGridPoints) = dCdt(IntAxGridPoints); %z = L
dCpdt = reshape(dCpdt,[],1);
dydt = [dCdt;dCpdt];
end
18 commentaires
Torsten
le 6 Sep 2023
Modifié(e) : Torsten
le 6 Sep 2023
Equation (3) with boundary condition (5) describe the diffusion from the surface to the interior of the particle.
To describe the rate of adsorption (adsorption isotherme), a simple model r_ads = constant * c*R*T is chosen where c is the concentration of the adsorptive in the bulk phase.
In your model, C_s (r_p) = c* would be an option where c* is the equilibrium concentration corresponding to the concentration of the adsorptive in the bulk phase.
Réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!