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

Question

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

当然，边界条件固定在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，因为你正在解决 Eq. (43) 和 Eq. (52) 之间的问题。
No-flux 是 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)

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

Solving Fick's second law of diffusion in 1D sphere using FiPy

python

pde

fipy