Variable and Term definition for structurally complicated PDEs

[zwang1995]

Hello all,

I am currently working with a complicated dynamic adsorption process using FiPy. I have trouble in representing equations using FiPy terms, due to the intricate equation structures. To help solve the problem, I only give a small part of the entire PDE system as a simple example, as shown below.

Overall mass balance (Eq1):
$$\frac{\partial P}{\partial t} = - P \frac{\partial v}{\partial z} - v \frac{\partial P}{\partial z}$$
Ergun equation (Eq2):
$$-\frac{\partial P}{\partial z} = k_a v + k_b v^2$$
$$k_a = \frac{150 \mu (1 - \epsilon)^2}{4 r_p^2 \epsilon^3}$$
$$k_b = \frac{1.75 (1 - \epsilon)}{2 r_p \epsilon^3} \frac{P}{R T_0}MW$$

All variables except for $P$ and $v$ are constant. The original PDEs don't have an equation containing a time derivative of $v$ (i.e., $\frac{\partial v}{\partial t}$), thus as far as I know $v$ cannot be considered as a variable in the FiPy framework.

It is well known that, the solution to a quadratic equation can be obtained by
$$x = \frac{-b \pm \sqrt{b^2 - 4 a c}}{2a}$$ (see Wikipedia quadratic equation).

After reformulating Eq2, a quadratic equation can be obtained ($k_a$ and $k_b$ are always positive, and basiclly $\frac{\partial P}{\partial z}$ should be negative due to the pressure drop):
$$k_b v ^ 2 + k_a v + \frac{\partial P}{\partial z} = 0$$

Therefore, $v$ can be represented as a function of $P$:
$$v = \frac{-k_a + \sqrt{k_a^2 - 4 k_b \frac{\partial P}{\partial z}}}{2 k_b}$$
Note: for ease of implementation only $+$ is considered here instead of $\pm$.

In this way, Eq2 can be eliminated. The PDE system has only one equation and one variable $P$.
$$\frac{\partial P}{\partial t} = - P \frac{\partial v}{\partial z} - v \frac{\partial P}{\partial z} = -\frac{\partial (vP)}{\partial z}$$
$$v = \frac{-k_a + \sqrt{k_a^2 - 4 k_b \frac{\partial P}{\partial z}}}{2 k_b}$$
$$k_a = \frac{150 \mu (1 - \epsilon)^2}{4 r_p^2 \epsilon^3}$$
$$k_b = \frac{1.75 (1 - \epsilon)}{2 r_p \epsilon^3} \frac{P}{R T_0}MW$$

In this way, I attempted to solve this problem using FiPy. My current concern is to figure out if this structurally complicated PDE system is correctly represented. This would be of great importance in helping formulate the original PDE system and find a potential solution to it.

Here I attach my code for debugging.

from fipy import CellVariable, Grid1D, TransientTerm, ConvectionTerm, Viewer

nx = 50
dx = 0.02
mesh = Grid1D(nx=nx, dx=dx)
P = CellVariable(name="Pressure", mesh=mesh, value=10000., hasOld=True)

P.constrain(400000., where=mesh.facesLeft)  # left boundary condition
P.faceGrad.constrain(0, where=mesh.facesRight)  # right boundary condition

# Ergun Equation
R = 8.314  # Universal gas constant [J/mol/K : Pa*m^3/mol/K]
T_0 = 313.15  # Temperature [K]
mu = 1.72e-5  # Viscosity of gas [Pa*s]
epsilon = 0.37  # Void fraction
r_p = 1e-3  # Radius of the pellets [m]
MW = 0.03042  # Molecular weight of gas [kg/mol]
k_a = 150. * mu * (1. - epsilon) ** 2. / 4. / r_p ** 2 / epsilon ** 3
k_b = 1.75 * (1. - epsilon) / (2. * r_p * epsilon ** 3) * (P / (R * T_0) * MW)
v = (-k_a + (k_a ** 2 - 4 * k_b * P.grad) ** 0.5) / (2 * k_b)

# Overall Mass Balance
eqOMB = TransientTerm(var=P) == (- ConvectionTerm(coeff=[[1.]] * v, var=P))

timeStepDuration = 0.01
steps = 20

viewer = Viewer(vars=P, datamin=0., datamax=1000000.)
viewer.plot()

desired_residual = 1e-2
for step in range(steps):
    P.updateOld()
    eqOMB.sweep(dt=timeStepDuration)
    viewer.plot()

Unfortunately, there are problems in solving the equation caused by the calculation of $v$. Nevertheless, I would like to know if such a strategy of representing $v$ as a function of $P$ is acceptable to FiPy. If yes, then I can proceed to solve the mathematical problems in the calculation of $v$.

Alternatively, I was wondering if there is any solution that can consider $v$ as a variable in FiPy. Therefore, the issues in the calculation of $v$ can be eliminated. This would definitely be better than my current strategy.

Let me know if any aspect was not clearly presented. Thank you :-)

[guyer]

A time derivative is not necessary for a FiPy equation, but I agree it's probably better to eliminate the Ergun equation as you have done.

I think what you have is OK and the reason that the solutions look wrong is because the time step is too large. Advective equations must obey the CFL condition, such that nothing is advected across more than one grid space in a time step. At $t=0$, $\Delta x / |v| \approx 0.00016$. If I set timeStepDuration = 0.0001, I get smoothly evolving solutions.

Additional things that might improve the solution:

• A VanLeerConvectionTerm tends to preserve shapes better for purely convective equations.
• Since v is a function of P.grad, sweeping may produce more accurate solutions.

[zwang1995]

Thank you so much for these comments and suggestions! I will implement these improvements and move on to the original PDE system.