基于FiPy的一维球对称菲克第二扩散定律数值求解方案

FiPy完全可以实现球对称一维径向扩散的求解，根本不需要手动推导坐标变换的离散格式，你现在的实现有几个关键的逻辑错位，改完就能跑出正确结果，核心说明如下：

原理校准 #

你要解的球坐标径向菲克第二定律：
$$\frac{\partial c}{\partial t} = \frac{D}{r^2}\frac{\partial}{\partial r}\left(r^2 \frac{\partial c}{\partial r}\right)$$
用FiPy自带的SphericalGrid1D网格时，完全不用手动给扩散项乘$r^2$的权重——网格离散的时候已经把不同r位置的球面面积、控制体体积的几何效应算进去了，直接写TransientTerm() == DiffusionTerm(coeff=D)的方程结构是对的，不存在形式不匹配的问题。之前参考的1D球网格实现踩坑，基本都是手动在笛卡尔网格上算几何权重才会出错，用官方的SphericalGrid1D根本没这个问题。

你当前代码的几个关键bug #

• 边界类型和需求错配：你一开始说球心、球面都是固定值边界，结果代码里给两个边界面都绑了梯度形式的通量约束，完全对不上。实际从物理上，球对称的球心r=0天然就是零通量（径向对称下不可能有物质穿过球心），FiPy的SphericalGrid1D在最左端的r=0面会自动应用这个零通量条件，根本不需要手动设置。
• 通量边界的物理关系搞反了：扩散通量的定义是$J = -D \cdot \frac{\partial c}{\partial r}$，你现在把算出来的表面通量直接塞给faceGrad，相当于把梯度值直接等于通量值，缺了除以扩散系数、调整符号的步骤，会让边界的物理效应直接放大十几个数量级，结果当然不对。
• 时间步长选得太离谱：你现在用1s的固定步长，虽然用的是隐式格式，但这么大的步长会让瞬态过程的精度完全丢失，甚至因为边界条件的错位导致结果发散。

正确的修正方式 #

对应不同的边界需求，按规则调整即可：

• 如果球面是固定值边界：直接用变量本身的约束语法，写X_ca.constrain(目标浓度, mesh.facesRight)，替换原来的faceGrad约束。
• 如果球面是固定通量边界：把梯度约束的数值改成-BoundaryR1_ca / D_ca，严格匹配通量和梯度的物理关系，也就是代码里写X_ca.faceGrad.constrain(-BoundaryR1_ca / D_ca, mesh.facesRight)。
• 时间步选合理值：你的网格尺寸大概是8.5e-8 m，按傅里叶数0.2的精度要求算，合适的时间步大概是100s量级，要捕捉早期瞬态的话可以再缩小。
  修正后的核心可运行代码：

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

nr = 100
r_ca = 8.5e-6
Iapp = -1.656
F = 96485.33212331001
S_ca = 1.1167
D_ca = 1e-14
Xinit_ca = 0.4952

# 计算球面的单位面积摩尔通量
BoundaryR1_ca = -Iapp/(S_ca*F)

mesh = SphericalGrid1D(nr=nr, Lr=r_ca) 
X_ca = CellVariable(mesh=mesh, name='Concentration cathode', value=Xinit_ca)
# 修正后的通量边界，匹配J=-D*grad的物理关系
X_ca.faceGrad.constrain(-BoundaryR1_ca / D_ca, mesh.facesRight)
# 球心r=0不需要额外设置，自动满足零通量
eq_ca = TransientTerm() == DiffusionTerm(coeff=D_ca)

tstep = 100
Nstep = 100

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

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

如果你确实需要在r=0位置设固定浓度（而不是自然零通量），只需要在左端面手动加一句X_ca.constrain(球心固定浓度值, mesh.facesLeft)就行，FiPy会自动处理对应的离散格式。

内容的提问来源于stack exchange，提问作者Oier Sasiain