时间戳

首页 > 解决方案 > 具有非线性源项/势项的扩散反应 pde

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

问题描述

我第一次尝试用 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是常数。我已经用 de Diffusion 项和非线性源项编写了代码。但是对于 100 的范围，我对 c2 有一些奇怪的行为。也许是我写错的非线性源项？我使用 ImplicitSourceTerm 命令，我认为这将线性化该术语。我错过了什么吗？我必须通过 myslef 进行线性化吗？（就像泰勒一样？）

如何将平流扩散反应 PDE 与 FiPy和 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))

也尝试使用命令扫描做一些数值分析，如残差和范数。我在这里看到了关于扫描命令的一个很好的解释： Solver tolerance and residual error when using sweep function in FiPy

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

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 对 pdes 来解决吗？是否也可以对耦合的 pdes 执行命令扫描？

或者我错过了什么？

我非常感谢任何帮助或建议。我希望，我写的很清楚。十分感谢。

标签： pythonpdefipy

解决方案

这是您的代码的工作版本。我不确定它是否正确，但它似乎更明智。

使用显式来源

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

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为负数时才能正常工作[编辑]。有关更多说明，请参阅此示例。

[编辑]：实际上，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否则它运行得非常慢。

未来

从长远来看，线性化源和使用ImplicitSourceTerms 以及耦合方程以提高效率可能是一个好主意，但只有在完全理解和验证解决方案之后。请参阅此示例和此说明以获取帮助。

推荐阅读