Time-dependent boundary conditions in adsorption process

[cge87]

Hello,
I want to use FiPy to simulate an adsorption process and have implemented the necessary equations. But I ran into 2 different problems when executing the calculations. I saw a similar discussion (#880) but I could not solve the issues based on the information. Further help to understand what is going on would be great!

Description of the system:
The process can be imagined using a pipe filled with adsorbent material and a fluid flowing through this pipe. 2 adsorbing substances are injected during a pre-specified time interval. These 2 substances compete for adsorption sites.
In the easiest case, the system can be modeled with 2 PDEs, one for each component, containing a transient term and a convection term:
TransientTerm + ConvectionTerm == 0
I wanted to model the injection time intervals of the 2 adsorbing substances with time-dependent boundary conditions "on the left side" of the pipe for variables var1 and var2. Below, you can find a code example:

#import
from fipy import *
from fipy.tools import numerix

# Define parameters
a1 = 0.8
a2 = 2.5
b = 0.00075
c1 = a1 / b
c2 = a2 / b
u = 0.025

#define a time variable
timeVar = Variable()

# Setting up the mesh
nx = 100
dx = 0.15
mesh = Grid1D(nx=nx, dx=dx)

# Define variables
var1 = CellVariable(mesh=mesh, name=r"$var1$", value=0.0)
var2 = CellVariable(mesh=mesh, name=r"$var2$", value=0.0)

# Define boundary conditions
#Left side
value_var1_left = (0.5 - 0.5*numerix.tanh(10*(timeVar-200)))*0.0015
var1.constrain(value_var1_left, mesh.facesLeft)
value_var2_left = (0.5 - 0.5*numerix.tanh(10*(timeVar-300)))*0.001
var2.constrain(value_var2_left, mesh.facesLeft)
#Right side
var1.faceGrad.constrain([0], mesh.facesRight)
var2.faceGrad.constrain([0], mesh.facesRight)

# Coefficients depending on the variables var1 and var2
Coeff_Transient_1 = 1 + (a1/(1+c1*var1+c2*var2)-var1*c1*a1/((1+c1*var1+c2*var2)**2))
Coeff_Transient_2 = 1 + (a2/(1+c1*var1+c2*var2)-var2*c2*a2/((1+c1*var1+c2*var2)**2))
Coeff_Convection_1 = u
Coeff_Convection_2 = u

# Equations
eq_1 = (TransientTerm(coeff=Coeff_Transient_1, var=var1) + VanLeerConvectionTerm(coeff=(Coeff_Convection_1,), var=var1) == 0)
eq_2 = (TransientTerm(coeff=Coeff_Transient_2, var=var2) + VanLeerConvectionTerm(coeff=(Coeff_Convection_1,), var=var2) == 0)

# Couple equations
eq = eq_1 & eq_2

# Set up time stepping
dt = 0.2
plot_interval = 10
timestep = 0
run_time = 600.0

# Viewer:
if __name__ == "__main__":
    viewer = Viewer(vars=(var1, var2), datamin=0., datamax=0.002)

while timeVar() < run_time:
    timeVar.setValue(timeVar() + dt)
    timestep += 1
    eq.solve(dt=dt)
    if (timestep % plot_interval ==0):
        if __name__ == '__main__':
            viewer.plot()

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

The 1st issue: I used the VanLeerConvectionTerm by chance and the time-dependent boundary conditions seem to work (in the sense that I at least see the "injection" stop). However, I would like to test other convection terms. But if I replace the VanLeerConvectionTerm e.g., by PowerLawConvectionTerm, the time-dependence of the left side boundary conditions is not captured any more. It just keeps the initial value. Any idea what is going on here?

The 2nd issue happens when I try to improve accuracy. I noticed that with the above code, the competition of the 2 substances for adsorption sites is not correctly accounted for. For example, the variable var1 needs to increase its value at the front (i.e. build-up of higher amount of the faster travelling substance 1 at the propagating front). This is accounted for if I change the above code to use "sweep": the front looks correct then. However, also here the time-dependence of the boundary conditions is not captured any more (even when I use the VanLeerConvectionTerm that seemed to work before). Any ideas what is going on and what to change? (Example code with sweep below)

Thanks in advance for any help and explanations!

#changed compared to first code version above: 
# "hasOld=True" for both variables
# "updateOld()" and sweep in the section where the model is solved

#import
from fipy import *
from fipy.tools import numerix

# Define parameters
a1 = 0.8
a2 = 2.5
b = 0.00075
c1 = a1 / b
c2 = a2 / b
u = 0.025

#define a time variable
timeVar = Variable()

# Setting up the mesh
nx = 100
dx = 0.15
mesh = Grid1D(nx=nx, dx=dx)

# Define variables
var1 = CellVariable(mesh=mesh, name=r"$var1$", value=0.0, hasOld=True)
var2 = CellVariable(mesh=mesh, name=r"$var2$", value=0.0, hasOld=True)

# Define boundary conditions
#Left side
value_var1_left = (0.5 - 0.5*numerix.tanh(10*(timeVar-300)))*0.0015
var1.constrain(value_var1_left, mesh.facesLeft)
value_var2_left = (0.5 - 0.5*numerix.tanh(10*(timeVar-300)))*0.001
var2.constrain(value_var2_left, mesh.facesLeft)
#Right side
var1.faceGrad.constrain([0], mesh.facesRight)
var2.faceGrad.constrain([0], mesh.facesRight)

# Coefficients depending on the variables var1 and var2
Coeff_Transient_1 = 1 + (a1/(1+c1*var1+c2*var2)-var1*c1*a1/((1+c1*var1+c2*var2)**2))
Coeff_Transient_2 = 1 + (a2/(1+c1*var1+c2*var2)-var2*c2*a2/((1+c1*var1+c2*var2)**2))
Coeff_Convection_1 = u
Coeff_Convection_2 = u

# Equations
eq_1 = (TransientTerm(coeff=Coeff_Transient_1, var=var1) + VanLeerConvectionTerm(coeff=(Coeff_Convection_1,), var=var1) == 0)
eq_2 = (TransientTerm(coeff=Coeff_Transient_2, var=var2) + VanLeerConvectionTerm(coeff=(Coeff_Convection_1,), var=var2) == 0)

# Couple equations
eq = eq_1 & eq_2

# Set up time stepping
dt = 0.2
plot_interval = 10
timestep = 0
run_time = 600.0

# Viewer:
if __name__ == "__main__":
    viewer = Viewer(vars=(var1, var2), datamin=0., datamax=0.002)

while timeVar() < run_time:
    timeVar.setValue(timeVar() + dt)
    timestep += 1
    var1.updateOld()
    var2.updateOld()
    for i in [1,2]:
        eq.sweep(dt=dt)

    if (timestep % plot_interval ==0):
        if __name__ == '__main__':
            viewer.plot()

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

[wd15]

The inlet boundary condition is a bit weird in FiPy. I'm not totally sure why it isn't working. However, the following code seems to work. It seems to capture the time dependence in the boundary conditions. Also, both the Van Leer and power law work. Obviously, the power law smooths things out a lot more. You don't seem to get the bump on x1 with the power law or it's a very small bump. Van Leer is probably giving the solution that you're looking for. BTW, FiPy isn't that great at hyperbolic equations in isolation.

These are the final curves using Van Leer.

In the following, x1 and x2 are constrained to 0 on the left and I've used a source of u * x1_left.divergence to account for the boundary. By setting x1 and x2 to 0 that creates a no inflow condition (you could also set the inlet velocity to 0) on the left for the convection terms and the source takes care of that. The x1_left face variable is updated at each time step.

#import
from fipy import *
from fipy.tools import numerix

# Define parameters
a1 = 0.8
a2 = 2.5
b = 0.00075
c1 = a1 / b
c2 = a2 / b
u = 0.025

#define a time variable
timeVar = Variable()

# Setting up the mesh
nx = 100
dx = 0.15
mesh = Grid1D(nx=nx, dx=dx)

# Define variables
x1 = CellVariable(mesh=mesh, name=r"$x1$", value=0.0, hasOld=True)
x2 = CellVariable(mesh=mesh, name=r"$x2$", value=0.0, hasOld=True)

# Define boundary conditions
#Left side
value_x1_left = (0.5 - 0.5*numerix.tanh(10*(timeVar-200)))*0.0015
x1.constrain(0, mesh.facesLeft)
value_x2_left = (0.5 - 0.5*numerix.tanh(10*(timeVar-300)))*0.001
x2.constrain(0, mesh.facesLeft)

x1_left = FaceVariable(mesh=mesh, value=0.0)
x1_left[0] = value_x1_left
x2_left = FaceVariable(mesh=mesh, value=0.0)
x2_left[0] = value_x2_left

#Right side
x1.faceGrad.constrain([0], mesh.facesRight)
x2.faceGrad.constrain([0], mesh.facesRight)

# Coefficients depending on the variables x1 and x2
Coeff_Transient_1 = 1 + (a1/(1+c1*x1+c2*x2)-x1*c1*a1/((1+c1*x1+c2*x2)**2))
Coeff_Transient_2 = 1 + (a2/(1+c1*x1+c2*x2)-x2*c2*a2/((1+c1*x1+c2*x2)**2))
Coeff_Convection_1 = u
Coeff_Convection_2 = u

# Equations
eq_1 = (TransientTerm(coeff=Coeff_Transient_1, var=x1) + VanLeerConvectionTerm(coeff=(Coeff_Convection_1,), var=x1) + u * x1_left.divergence == 0)
eq_2 = (TransientTerm(coeff=Coeff_Transient_2, var=x2) + VanLeerConvectionTerm(coeff=(Coeff_Convection_1,), var=x2) + u * x2_left.divergence == 0)

#eq_1 = (TransientTerm(coeff=Coeff_Transient_1, var=x1) + PowerLawConvectionTerm(coeff=(Coeff_Convection_1,), var=x1) + u * x1_left.divergence == 0)
#eq_2 = (TransientTerm(coeff=Coeff_Transient_2, var=x2) + PowerLawConvectionTerm(coeff=(Coeff_Convection_1,), var=x2) + u * x2_left.divergence == 0)

# Couple equations
eq = eq_1 & eq_2

# Set up time stepping
dt = 0.2
plot_interval = 10
timestep = 0
run_time = 600.0

# Viewer:
if __name__ == "__main__":
    viewer = Viewer(vars=(x1, x2), datamin=0., datamax=0.002)

while timeVar() < run_time:
    timeVar.setValue(timeVar() + dt)
    timestep += 1
    x1.updateOld()
    x2.updateOld()
    x1_left[0] = value_x1_left
    x2_left[0] = value_x2_left
    for i in [1,2]:
        eq.sweep(dt=dt)

    if (timestep % plot_interval ==0):
        if __name__ == '__main__':
            viewer.plot()

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