Adsorption Coupled System

[agoncalvescoel]

Hi all, I'm learning how to use FiPY and wanted to solve the following adsorption column problem:

(transient term dq is missing rho in the formula)

The code for now looks like:

from fipy import Variable, FaceVariable, CellVariable, Grid1D, ExponentialConvectionTerm, ExplicitDiffusionTerm, ImplicitSourceTerm, TransientTerm, DiffusionTerm,  PowerLawConvectionTerm, Viewer
from fipy.tools import numerix
import matplotlib.pyplot as plt
import numpy as np

Dz = 5.4E-4 # m^2/s
vDot = 1.83E-5 # m^3/s
dBed = 0.03 # m
epsBed = 0.808
u = vDot/(epsBed*3.14*(dBed/2)**2) # m/s
print("velocity "+ str(u))
rhoP = 2505.4 # kg/m^3
kG = 2.87E-3 # 1/s
qS = 0.35 # mol/kg
kL = 31.72 # m^3/mol
x0 = 162E-6 # molar fraction of adsorbate on inlet
rhoGas = 1.1012 # kg/m^3
MWH2S = 34.082*1E-3 # kg/mol, molecular weight of H2S
c0 = x0*rhoGas/MWH2S # mol/m^3

convCoeff = (u,) # velocity vector, m/s
L = 0.1
nx = 200
dx = L/nx
m = Grid1D(nx=nx, dx=dx)
timeStepDuration = 0.9 * dx**2 / (2*Dz)
steps = 1000

C = CellVariable(name = "C", mesh=m, hasOld=True, value = 0.)
q = CellVariable(name ="q", mesh=m, hasOld=True, value = 0.)

C.constrain(c0, where=m.facesLeft)

eq0 = TransientTerm(var=C) == DiffusionTerm(coeff=Dz, var=C) -  ExponentialConvectionTerm(coeff=convCoeff, var=C) - ((1-epsBed)/epsBed)*(q-q.old)/timeStepDuration
eq1 = TransientTerm(var=q) == kG*qS*kL*C/(1+kL*C) - ImplicitSourceTerm(coeff=kG, var=q)

eqn = eq0 + eq1

vi = Viewer((C, q), datamax = c0)

from builtins import range
for step in range(steps):
    C.updateOld()
    q.updateOld()
    eqn.solve(dt=timeStepDuration)
    vi.plot()

from fipy import input
if __name__ == '__main__':
    input("Press <return> to proceed...")

The output for C has the format of the curve I would expect. But the values increase quicker than I do when I solve the same model in COMSOL. And q does not correspond to COMSOL at all, and doesn't vary either with space or time. Can someone help me to what I'm possibly missing here?

COMSOL:

[guyer]

• eqn = eq0 + eq1 adds two equations together, producing the single equation

  TransientTerm(var=C) + TransientTerm(var=q) == DiffusionTerm(coeff=Dz, var=C) -  ExponentialConvectionTerm(coeff=convCoeff, var=C) - ((1-epsBed)/epsBed)*(q-q.old)/timeStepDuration + kG*qS*kL*C/(1+kL*C) - ImplicitSourceTerm(coeff=kG, var=q)

  rather than the coupled pair of equations you wanted. I'm surprised this doesn't generate an error, because you're not specify the single variable you want this single equation to solve for.

  Use eqn = eq0 & eq1.

• Potentially better coupling can be obtained with

  eq0 = (TransientTerm(var=C)
         == DiffusionTerm(coeff=Dz, var=C) 
         - ExponentialConvectionTerm(coeff=convCoeff, var=C)
         - TransientTerm(coeff=((1-epsBed)/epsBed), var=q))
  eq1 = (TransientTerm(var=q)
         == ImplicitSourceTerm(coeff=kG*qS*kL/(1+kL*C), var=C)
         - ImplicitSourceTerm(coeff=kG, var=q))

• Because eq1 is nonlinear in C, you should sweep. Replace

  --- a/langmuir.py
  +++ b/langmuir.py
  @@ -42,7 +42,8 @@ from builtins import range
   for step in range(steps):
       C.updateOld()
       q.updateOld()
  -    eqn.solve(dt=timeStepDuration)
  +    for sweep in range(4):
  +        res = eqn.sweep(dt=timeStepDuration)
       vi.plot()
   
   from fipy import input

  I picked four sweeps because the residual doesn't really get any smaller after that.

• FiPy's default boundary condition is zero flux ($D_z\partial C/\partial z - \vec{u} C = 0$), not zero gradient ($\partial C/\partial z = 0$). Add
  C.faceGrad.constrain([[0.]], where=m.facesRight)

• You can take much bigger time steps

Beyond that, I don't know why this doesn't produce the same results as in COMSOL. I would take a close look to ensure that all of the parameters are the same. I definitely don't see anything here that would give rise to the structure that COMSOL shows in q.

[agoncalvescoel]

Thanks for the help how how to use the tool better.
I have tried with your suggestions, and now q does change within space and time. But I would say still "quicker" in time, I'm checking the units I used in both. (I'm following the model from the paper http://dx.doi.org/10.1016/j.cej.2015.12.075)

The breakthrough curve given by the last run is, for C0 = 162 ppm

From the paper, I would expect the following:

[agoncalvescoel]

I do agree that the result for q in COMSOL is a bit off, I'm re-running it

[guyer]

Thanks for the reference. I'll take a look later today.

[agoncalvescoel]

The q looks maybe a problem of boundary conditions in COMSOL, because the bumps from the beginning are reduced by increasing the mesh points.

[guyer]

I'm still working through coefficients, but I think one of the coefficients below Eq. (12) isn't right in the paper you linked. I think it should be

$$a = \begin{bmatrix} 0 & 0 \\ 0 & S \end{bmatrix}$$

[agoncalvescoel]

Yep, I have corrected that when I ran too. I didn't ran it dimensionless like the model, actually, I adapted everything

[guyer]

Eq. (1) is not dimensionally consistent. Every term except the last is in units of $\mathrm{\frac{mol / m^3}{s}}$. The last term, $\frac{1 - \varepsilon}{\varepsilon}\frac{\partial q}{\partial t}$, is in units of $\mathrm{\frac{mol / kg}{s}}$.

[agoncalvescoel]

This is the outputs I get from running with the density :/

Not really sure on what I'm missing here.

[guyer]

If I change the final term in eq0 to TransientTerm(coeff=((1-epsBed)/epsBed) * rhoP, var=q), I get stable solutions with much slower evolution in $C$. Based on what you're seeing, I suspect you're taking too large a time step.

When I said "You can take much bigger time steps", I meant maybe a factor of 100. The criterion you use, timeStepDuration = 0.9 * dx**2 / (2*Dz), is appropriate for an explicit diffusion problem. FiPy solves diffusion implicitly, so this stability criterion does not apply. FiPy is unconditionally stable (although not unconditionally accurate) for diffusion. You have a convection-diffusion problem, so your time step is limited by the CFL limit that says your convection velocity field cannot transport anything more than one grid point in a time step, so $\Delta t \approx \Delta x / u$. We typically limit the time step to 10% of that value, although I find (and the Wikipedia article supports) that larger values are stable. For your problem, timeStepDuration = 0.05 seems OK. Accuracy is left as an exercise for the reader.

[guyer]

In addition to the "bumps", your COMSOL solution appears to constrain $q$ to zero at the left-hand side. I don't see this constraint in the paper you linked, your problem statement, or your code.

[guyer]

What is "breakthrough" exactly? The term is used throughout the paper, but never defined AFAICT. Is it $C(x=L)$, $q(x=L)$, or something else?

[agoncalvescoel]

Thanks for the inputs. I'm running it again now, let's see how it will behave.
I'm sorry, breakthrough curves is the curve of the concentration (C) at x=L in respect to time (t)

[agoncalvescoel]

Hi hi, I got to run the script, but it's taking really long.... more than 24 hrs with the sweep, I deleted the sweep without damaging results, and now it takes 10 hrs to compute. But it still take too long to compute. I have found the codes used for iterating the time steps, but on the documentation is missing the information on how to use those. Which is the default? How can I change it?

[guyer]

Sorry for the delay in replying.

I have found the codes used for iterating the time steps, but on the documentation is missing the information on how to use those. Which is the default? How can I change it?

I don't know what you're referring to