PDE solution not changing with time

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

[Torsten]

Remember your axial length is 4 cm and your flow velocity is 3,6 cm/s. I wouldn't be surprised if your system gets stationary almost instantly.

[James99]

yes but that doesn't explain the abrupt rising of the values, in fact the values shouldn't rise above C0 at all. If I am acheiving breakthrough that fast then all the values in the last column of CFinal should be 0.0437. Also increasing the bed length should fix that.

By the way, my goal is to plot the last column of CFinal and it should be like a top left corner of a square but rounded.

[Torsten]

I don't have the time to look through your code thoroughly, but how should this line, e.g., be correct ?

dCdt(z) = -(u)*dCdz + (Dz)*d2Cdz2 - M.*(C(z) - dCpdt(z,end));

You subtract mol/(m^3 * s) from mol/m^3 here.

Same here:

dCpdt(:,1) = 4*Cp(:,2)/3 - Cp(:,3)/3;

Unit on the left-hand side: mol/(m^3 * s). Unit on the right-hand side: mol/m^3.

[James99]

The equation seems to be correct (see attached). Although dCp/dt does have the units of ug/L/s as opposed to C's ug/L, I wrote dCp/dt(z,end) to extract the values of concentration at the particle surface. Ideally I wanted it to be Cp(z,end) but since it needed to calculate the concentration in the particle from r = 0 until the r = rp-dr, I could only set the boundary condition after the calculations were done, that's what I thought. So setting the boundary condition that way would make sense. Anyway, I changed the code a bit. So you can check it out if you have time.

In any case, I just don't understand why the values repeat themselves and why there are abrupt rise in values about halfway through the axial direction.

% 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 ====================================================================
TimeStamps = 100;
t = linspace(0,1440,TimeStamps);


% Initial Conditions ======================================================
% reactor
C = zeros(1,AxialGridPoints);
C(1) = C0;

% Particle
Cp = zeros(1,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);
[t, y] = ode15s(@PoreDiffusion,t,Initial,[],Data,Spatial);
%beep

% Result ==================================================================

CFinal = y(:,1:AxialGridPoints);
CpFinal = y(:,AxialGridPoints+1:end); 



ggd = zeros(PartGridPoints,AxialGridPoints,TimeStamps);
%ggd contains 100 pages (1 for each time) with 21 rows for each particle
%coordinate and 201 columns for each axial coordinate

for f = 1:TimeStamps
    for j = 1:PartGridPoints
        ggd(j,:,f) = CpFinal(f, AxialGridPoints*(j-1) + 1:j*AxialGridPoints);
    end
end

plot(t,CFinal(:,end))

% 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;
                L1 = (2*qmax*KL*KL*Ds*rho_b)/(1+KL*Cp(z,r))^3;
                L2 = (3*De)/(2*dr) + kf;

                Cp(z,1) = 4*Cp(z,2)/3 - Cp(z,3)/3;
                Cp(z,end) = C(z)*kf/L2 + Cp(z,IntPartGridPoints)*(2*De/dr/L2) - Cp(z,IntPartGridPoints-1)*(De/2/dr/L2);
                
                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;
                dCpdt(z,r) = De/x*d2Cpdr2 - L1/x*dCpdr^2 + (2*De)/(x*Radius(r))*dCpdr;
            end

            % Reactor
            
            % 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) - Cp(z,end));
        end
        
      
        % Reassigning boundary conditions
        dCdt(AxialGridPoints) = dCdt(IntAxGridPoints); %z = L
                
        dCpdt = reshape(dCpdt,[],1);
        dydt = [dCdt;dCpdt];
end

[Torsten]

Where is the equation of the adsorption isotherm in your code which is needed for the driving force at the particle surface ?

[James99]

If you mean the equation for film diffusion, then that's integrated in the dCdt equation (M.*(C(z) - Cp(z,end)).

[Torsten]

I mean the equilibrium adsorption isotherm. It's independently given from measurements and describes the loading q_star when adsorptive and adsorbat are in equilibrium. As long as you don't have this isotherm for your system, you don't need to start simulating - it's the core of the simulation.

[James99]

I understand that the particle level adsorption is linked to the fluid phase equation using the adsorbate concentration at the outer surface of the adsorbent particle. So given that the two equations already have the concentration at the adsorbent surface term, I thought they could be directly linked. I would need the equilibrium equation if the pore diffusion equation was in terms of amount of adsorbate adsorbed. You can find this information in the DOI of the article on the first line of the code (it's commented).

I apologize for the late reply.

[Torsten]

In any case, you need the equilibrium concentration of the adsorbate as a function of the adsorptive concentration (and maybe temperature and pressure). How else should the maximum level of concentration of the adsorbate within the packed bed come out of your simulation if you don't feed the simulation with this input ?

Imagine one packed bed has a capacity of 0 kg adsorbate/kg solid and the other has a capacity of 10 kg adsorbate/kg solid for the same adsorptive concentration. How should this be a result of your simulation ? It's a characteristic of the adsorbens.

You cannot master to simulate something useful if you don't study the physical background first. That's what I already gave you as advice in a previous response.

https://books.google.de/books/about/Adsorptionstechnik.html?id=I2oaBQAAQBAJ&redir_esc=y

https://onlinelibrary.wiley.com/doi/10.1002/bbpc.19900940122

[James99]

Alright mate,

I think I do have pretty good notion of the physical processes but perhaps I missed something.

[James99]

Hello Torsten,

I have a couple of questions, I've seen multiple diffusion models for adsorption in spherical particles and they all seem to link with the fluid phase mass transfer with the concentration at the outer surface of the adsorbent particle. More specifically, the diffusion models' boundary condition at particle radius and liquid phase transport equation's uptake rate term. Both of the equations have this unknown and it can't be solved without it, so how would MATLAB solve this?

And for the second question, let me just rewrite the equations I've used in my code:

(Solved for dCp/dt)

Do you reckon these equations could be solved using pdepe?

And finally, thank you for all your help

[Torsten]

Both of the equations have this unknown and it can't be solved without it, so how would MATLAB solve this?

I don't understand what you mean.

To your second question:

You cannot use pdepe because your equations contain two spatial independent variables, namely z and r.

And the important thing to know would be the boundary condition for C_p at r = r_p.

Maybe the example "Transport and Adsorption" in the

COMSOL Multiphysics Model Library Manual

can help to understand things better.

[James99]

I don't know how to explain this, concentration on the surface of the particle (Cp|r = rp) is essentially an unknown which also connects the two equations. And it being part of the boundary condition of the second equation means that the second equation is not well defined. And since we don't have the values of Cp|r = rp, the first equation can't be solved. So how does MATLAB calculate Cp|r = rp? I hope you are understanding what I am trying to say.

Also, thank you for the link, appreciate it mate

[Torsten]

So how does MATLAB calculate Cp|r = rp?

You mean in your code ?

[James99]

Well yes, but also generally. How would MATLAB calculate that? I mean there are people who have written the code, pressumably in the correct way and I know there is no way to explicitly write an expression for Cp|r = rp.

[Torsten]

You must set how much adsorptive becomes adsorbate (by using the adsorption isotherm). Then the condition at the surface becomes defined.

Study my link. Everything is explained in the COMSOL example.

COMSOL is by the way the software best suited to handle your problem. Maybe ASPEN has similar features, I don't know.

[James99]

Very well, you were excellent, thank you!

[Torsten]

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.