使用 FiPy 求解一维球体中的菲克第二扩散定律

Question

我正在尝试使用 fipy 解决球体 PDE 的第二扩散定律。我已经检查了文档，但是没有这种情况的示例，所以我想知道是否真的可以这样做，因为我没有成功地获得足够的方程定义结构。我考虑了方位角和天顶角对称性，因此我需要求解的方程如下所示。

当然，边界条件固定在 r=0 和 r=R 固定值，并且初始浓度也是已知的。我也尝试遵循here 中提出的想法，但没有得到任何明确的结果。欢迎任何想法。

我目前使用的代码如下：

from fipy import Viewer, TransientTerm, DiffusionTerm, SphericalGrid1D, CellVariable
from builtins import range

nr = 100 #number of cells on the mesh
r_ca = 8.5e-6 #sphere radius

Iapp = -1*1.656 #Current [A] discharge
F = 96485.33212331001 #Faraday constant
S_ca = 1.1167 #Surface of cathode
D_ca = 1e-14 #Diffusion coef.

Xinit_ca = 0.4952 #[p.u.]
Xinit_an = 0.7522

BoundaryR0 = 0 #Fixed flux at r=0
BoundaryR1_ca = -Iapp/(S_ca*F) #Fixed flux at r=r_ca

mesh = SphericalGrid1D(nr=nr, Lr=r_ca) 
X_ca = CellVariable(mesh=mesh, name='Concentration cathode', value=Xinit_ca) 
X_ca.faceGrad.constrain(BoundaryR1_ca, mesh.facesRight) #Fixed flux boundary condition
X_ca.faceGrad.constrain(BoundaryR0, mesh.facesLeft)

eq_ca = TransientTerm() == DiffusionTerm(coeff=D_ca) # dX/dt = D/r^2 * d/dr(r^2*dX/dr)

tstep = 1 #s
Nstep = 1000

for step in range(Nstep):
    eq_ca.solve(var=X_ca, dt=tstep)

viewer = Viewer(vars=X_ca, datamin=0.45, datamax=0.8)

Answer 1

发生了一些事情：

• 一些求解器不喜欢球形网格，可能是因为单元体积范围很大。 SciPy LinearLUSolver 似乎有效。明智的预处理可能会对其他求解器有所帮助。
• 等式。 (45) 在the paper you linked below 中定义了 flux，但您正在约束梯度。您需要除以扩散率。
• X_ca 以[stoichiometry] 为单位，但BoundaryR1_ca 以mol/(m**2 * s) 为单位（或mol/m**4 除以D_ca 后。您需要除以C_ca_max 以得到[stoichiometry]/m，因为你'正在解决公式 (43) 和公式 (52) 之间的问题。
• 无通量是 FiPy 的自然边界条件，因此无需限制在 mesh.facesLeft。
• mesh.facesRight 处的梯度应限制为向量（即使是一维向量）。

通过以下更改，我得到的行为看起来与Mayur et al. 中的图 7 一致。

diff --git a/spherical.py b/spherical.py
index e6c2969..1eee346 100644
--- a/spherical.py
+++ b/spherical.py
@@ -1,4 +1,4 @@
-from fipy import Viewer, TransientTerm, DiffusionTerm, SphericalGrid1D, CellVariable
+from fipy import Viewer, TransientTerm, DiffusionTerm, SphericalGrid1D, CellVariable, LinearLUSolver
 from builtins import range
 
 nr = 100 #number of cells on the mesh
@@ -13,19 +13,22 @@ Xinit_ca = 0.4952 #[p.u.]
 Xinit_an = 0.7522
 
 BoundaryR0 = 0 #Fixed flux at r=0
-BoundaryR1_ca = -Iapp/(S_ca*F) #Fixed flux at r=r_ca
+BoundaryR1_ca = -Iapp/(51410*D_ca*S_ca*F) #Fixed flux at r=r_ca
 
 mesh = SphericalGrid1D(nr=nr, Lr=r_ca) 
 X_ca = CellVariable(mesh=mesh, name='Concentration cathode', value=Xinit_ca) 
-X_ca.faceGrad.constrain(BoundaryR1_ca, mesh.facesRight) #Fixed flux boundary condition
-X_ca.faceGrad.constrain(BoundaryR0, mesh.facesLeft)
+X_ca.faceGrad.constrain([BoundaryR1_ca], mesh.facesRight) #Fixed flux boundary condition
 
 eq_ca = TransientTerm() == DiffusionTerm(coeff=D_ca) # dX/dt = D/r^2 * d/dr(r^2*dX/dr)
 
 tstep = 1 #s
 Nstep = 1000
+solver = LinearLUSolver()
+
+viewer = Viewer(vars=X_ca, datamin=0.45, datamax=0.8)
 
 for step in range(Nstep):
-    eq_ca.solve(var=X_ca, dt=tstep)
+    eq_ca.solve(var=X_ca, dt=tstep, solver=solver)
+    if step % 100 == 0:
+        viewer.plot()
 
-viewer = Viewer(vars=X_ca, datamin=0.45, datamax=0.8)