具有非线性源项/势项的扩散反应 pde

Question

我第一次尝试用 FiPy 写几个扩散反应方程。

我的方程式适用于三种不同的浓度 c_i，仅包含扩散部分和源（LateX 输入）：

\partial_t c_1 = \nabla \cdot (D_1 \nabla c_1) + A_1 c_{1} + B_1 c_2^4 + W_1 \
\partial_t c_2 = \nabla \cdot (D_2 \nabla c_2) + c_{1}(A_2 + 2c_1) + c_2^4 + W_2\
\partial_t c_3 = \nabla \cdot (D_3 \nabla c_3) + A_3 c_{1} - W_3 \

系数D_i、A_i、B_i、P_i和W_1是常量。我已经用去扩散项和非线性源项编写了代码。但是对于 100 的范围，我得到 c2 的一些奇怪行为。也许是非线性源项，我写错了吗？我使用 ImplicitSourceTerm 命令，我认为这将线性化该术语。我错过了什么吗？我必须通过 myslef 进行线性化吗？（就像泰勒？）

和 https://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.coupled.html

from fipy import *
from fipy import CellVariable, Variable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, LinearLUSolver, Viewer
import math
from fipy.tools import numerix

# Konstanten
F       = 96485.3399  #C /mol
R       = 8.314472    # J/(mol*K) 
T       = 273.18      # K
alpha   = 0.5
c_std   = 1.           #  mol/s
c1_sat  = 1.           #  mol/s
eta_Zn  = 1.           # V
kIc = 6.3      # mol/s Von Reaktion I an der Kathode
kIa = 5.       # mol/s Von Reaktion I an der Anode
kII = 0.25        # mol/s von Reaktion II
mu_1I   = 1.    # Von Zinkat in Reaktion I
mu_1II  = -2.   # Von Zinkat in Reaktion II
mu_2I   = -4.   # Von OH in Realtion I
mu_2II  = 2.    # Von OH in Reaktion II
mu_3II  = 1.   # Von H2O in Reaktion II
ep1 = 1.    # Von Zinkat
ep2 = 1.   # Von OH
ep3 = 1.   # Von H20
D1 = .5     # Von Zinkat
D2 = 1.     # Von OH
D3 = 0.1    # Von H20

## Mesh
L = 10.
nx = 1000
mesh = Grid1D(dx=L/1000, nx=nx)
x = mesh.cellCenters[0]

## Initial Conditions
c1 = CellVariable(name="c1", mesh=mesh, value=1., hasOld=True)
c2 = CellVariable(name="c2", mesh=mesh, value=1., hasOld=True)
c3 = CellVariable(name="c3", mesh=mesh, value=1., hasOld=True)

## Boundary Conditions
c1.constrain(2., mesh.facesLeft)
c1.constrain(0., mesh.facesRight)
c2.constrain(0., mesh.facesLeft)
c2.constrain(2., mesh.facesRight)
c3.constrain(0., mesh.facesLeft)
c3.constrain(2., mesh.facesRight)
c2.faceGrad.constrain(1., where=mesh.facesRight)

# Definition Konstanten fuer SourceTerm
Zn1 = (kIc/c_std)*math.exp(-(1-alpha)*(F/(R*T))*eta_Zn)+kII
Zn2 = (kIa/(c_std)**4)*math.exp(alpha*(F/(R*T))*eta_Zn)
Zn3 = kII*c1_sat
OH1 = 4*(kIc/c_std)*math.exp(-(1-alpha)*(F/(R*T))*eta_Zn)+2*kII
OH2 = 4*(kIa/(c_std)**4)*math.exp(alpha*(F/(R*T))*eta_Zn)
OH3 = 2*kII*c1_sat
H1 = kII
H2 = 0.
H3 = kII*c1_sat

sourceCoeff1 = (Zn1*c1) + (Zn1*c2**4) + Zn3
sourceCoeff2 = c1*(OH1 + 2*c1) + c2**4 + OH3
sourceCoeff3 = c1*H1 - H3

eq1 = (TransientTerm(var=c1) == DiffusionTerm(coeff=D1, var=c1) + ImplicitSourceTerm(coeff=sourceCoeff1, var=c1))
eq2 = (TransientTerm(var=c2) == DiffusionTerm(coeff=D2, var=c2) + ImplicitSourceTerm(coeff=sourceCoeff2, var=c2))
eq3 = (TransientTerm(var=c3) == DiffusionTerm(coeff=D3, var=c3) + ImplicitSourceTerm(coeff=sourceCoeff3, var=c3))

eqn = eq1 & eq2 & eq3


vi = Viewer((c1, c2, c3))

for t in range(50):
    c1.updateOld()
    c2.updateOld()
    c3.updateOld()
    eqn.solve(dt=1.e-3)
    vi.plot()

我也在尝试添加潜在术语，但我不明白这是否可能。类似的东西（与上面相同，但有潜在的部分）

\partial_t c_1 = \nabla \cdot (D_1 \nabla c_1) + \nabla \cdot \biggl(P_1c_1\nabla\Phi\biggr) + A_1 c_{1} + B_1 c_2^4 + W_1 \
\partial_t c_2 =  \nabla \cdot (D_2 \nabla c_2) + \nabla \cdot \biggl(P_2c_2\nabla\Phi\biggr) + c_{1}(A_2 + 2c_1) + c_2^4 + W_2\
\partial_t c_3 = \nabla \cdot (D_3 \nabla c_3) + \nabla \cdot \biggl(P_3c_3\nabla\Phi\biggr) + A_3 c_{1} - W_3 \

我分析了来自 https://www.ctcms.nist.gov/fipy/examples/phase/generated/examples.phase.binaryCoupled.html 尝试添加像 Diffusionterm 这样的潜在术语，但我总是得到错误：

SolutionVariableNumberError: 解变量和方程的数量不同。

我所期望的，但也许我遗漏了什么。

请参阅该部分代码：

eq1 = (TransientTerm(var=c1) == DiffusionTerm(coeff=D1, var=c1) + ImplicitSourceTerm(coeff=sourceCoeff1, var=c1)) + DiffusionTerm(coeff=D1phi, var=phi
eq2 = (TransientTerm(var=c2) == DiffusionTerm(coeff=D2, var=c2) + ImplicitSourceTerm(coeff=sourceCoeff2, var=c2)) + DiffusionTerm(coeff=D2phi, var=phi)
eq3 = (TransientTerm(var=c3) == DiffusionTerm(coeff=D3, var=c3) + DiffusionTerm(coeff=D3phi, var=phi) + ImplicitSourceTerm(coeff=sourceCoeff3, var=c3))

也尝试使用命令sweep做一些数值分析，比如残差和范数。我在这里看到了关于扫描命令的很好的解释：

那部分代码（仅用于带有源的扩散项）

dt = 1.e5

solver = LinearLUSolver(tolerance=1e-10)

c1.updateOld()
c2.updateOld()
c3.updateOld()
res = 1.
initialRes = None

while res > 1e-4:
     res = eq.sweep(dt=dt, solver=solver)
     if initialRes is None:
         initialRes = res
     res = res / initialRes

然而我没有得到任何结果或错误，必须手动停止该过程。

总而言之，我的问题是：有可能用带有扩散项和潜在项的 FiPy couple pdes 来解决吗？耦合 pdes 也可以执行命令扫描吗？

或者我遗漏了什么？

非常感谢任何帮助或建议。我希望，我写清楚了。非常感谢。

Answer 1

这是一个working version of your code。不知道对不对，但是好像比较有道理

使用明确的来源

所做的主要更改是删除隐式来源术语并用显式来源替换它们。例如第一个等式

eq1 = (TransientTerm(var=c1) == DiffusionTerm(coeff=D1, var=c1) + ImplicitSourceTerm(coeff=sourceCoeff1, var=c1))

现在

eq1 = (TransientTerm() == DiffusionTerm(coeff=D1) + sourceCoeff1)

在 OP 代码的所有三个等式中使用 ImplicitSourceTerm 都是不正确的。要了解其用法，请使用 value * phi 形式的来源，其中 phi 是要求解的变量。在这种情况下，正确的来源是 ImplicitSourceTerm(value)。 phi 的乘法是隐式假定的。此外，这仅在 value 为负数时才能正常工作 [edit]。有关更多说明，请参阅 this example。

[编辑]：实际上，ImplicitSourceTerm 处理正值和负值，因此如果值为正，则明确添加源。

解耦方程

最好先解开方程式，以便更好地理解发生了什么。当在单个方程中耦合的空间导数中有多个变量时，耦合方程最有效。

在这种情况下，根据源项耦合方程式可能很有用，但仅出于效率目的且仅在验证问题解决方案之后。新代码现在以非耦合形式求解方程。

for t in range(100):
    c1.updateOld()
    c2.updateOld()
    c3.updateOld()
    for sweep in range(5):
        res1 = eq1.sweep(c1, dt=dt, solver=solver)
        res2 = eq2.sweep(c2, dt=dt, solver=solver)
        res3 = eq3.sweep(c3, dt=dt, solver=solver)
        print('res1', res1)
        print('res2', res2)
        print('res3', res3)

    vi.plot()

使用扫描循环

由于方程在源代码中是非线性的，因此最好使用内部循环并打印残差以检查上面代码中每个时间步的收敛性。

使用不同的求解器

为了使事情正常进行，似乎需要与 FiPy 选择的默认求解器不同的求解器。这通常取决于所使用的求解器套件。在这种情况下，我使用的是 Python 3，所以我可能使用的是 Scipy 求解器。为了让它工作，我使用了 LinearLUSolver 否则它运行非常慢。

未来

I long 运行将源线性化并使用 ImplicitSourceTerm 以及耦合方程以提高效率可能是个好主意，但只有在解决方案完全理解和验证。请参阅 this example and this explanation 寻求帮助。

具有非线性源项/势项的扩散反应 pde

Diffusion Reaction pde with nonlinear source term / Potential Term

python

pde

fipy

使用明确的来源

解耦方程

使用扫描循环

使用不同的求解器

未来