使用 FiPy 求解圆柱坐标系中的扩散方程不正确

Question

使用 FiPy 对圆柱形几何扩散方程的稳态解与从其他软件（例如 Mathematica）获得的解有很大不同。

等式为：

$0 = \frac{1}{r}\frac{d}{dr}\left(\frac{r}{T^{1/2}}\frac{dT}{dr}\right ) + cte*T^{3/2}$

这意味着，通过使用 CylindricalGrid1D 网格，我们可以将等式写为：

mesh = CylindricalGrid1D(nr=100, dr=0.01, origin=0.0)
T = CellVariable(name='temperature', mesh=mesh, hasOld=True)
r = mesh.cellCenters()

#BC's
T.constrain(0., mesh.facesRight)
T.faceGrad.constrain(0.,mesh.facesLeft)

#initial temperature profile
T.setValue( 1-r**2)

eq = 0 == DiffusionTerm( coeff=T**(-1/2), var=T) + 20*ImplicitSourceTerm(coeff=T**(1/2), var=T)

viewer = Viewer(vars=T)
eq.solve()

viewer.plot()
raw_input(" Press <enter> to proceed...")

这里我设置了cte=20，但是不管这个值是多少，问题依旧。我得到的解决方案在左边，而 Mathematica 给出的解决方案在右边：

plots

然后我尝试按照推荐的方法扫描像这样的非线性方程。所以我没有 eq.solve()，而是：

current_residual = 1.0e100
desired_residual = 1e-5

while current_residual > desired_residual:

    current_residual = eq.sweep()
    T.updateOld()

但我得到错误：

/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py:66: RuntimeWarning: overflow encountered in square
  error0 = numerix.sqrt(numerix.sum((L * x - b)**2))
/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py:71: RuntimeWarning: overflow encountered in square
  if (numerix.sqrt(numerix.sum(errorVector**2)) / error0)  <= self.tolerance:
/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py:71: RuntimeWarning: invalid value encountered in double_scalars
  if (numerix.sqrt(numerix.sum(errorVector**2)) / error0)  <= self.tolerance:
/home/antonio/.local/lib/python2.7/site-packages/fipy/tools/numerix.py:966: RuntimeWarning: overflow encountered in square
  return sqrt(add.reduce(arr**2))
/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py:58: RuntimeWarning: overflow encountered in multiply
  b = b * (1 / maxdiag)
Traceback (most recent call last):
  File "stack.py", line 26, in <module>
    current_residual = eq.sweep()
  File "/home/antonio/.local/lib/python2.7/site-packages/fipy/terms/term.py", line 254, in sweep
    solver._solve()
  File "/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/scipySolver.py", line 61, in _solve
    self.var[:] = numerix.reshape(self._solve_(self.matrix, self.var.ravel(), numerix.array(self.RHSvector)), self.var.shape)   
  File "/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py", line 64, in _solve_
    permc_spec=3)
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 257, in splu
    ilu=False, options=_options)
RuntimeError: Factor is exactly singular

最后，我使用变量 V=T^{1/2} 以等价形式重写了初始方程。很容易看出，对于 V，等式变为

$0 = \frac{1}{r}\frac{d}{dr}\left(r\frac{dV}{dr}\right) + \frac{cte}{2}V ^3$

于是我使用了代码：

mesh = CylindricalGrid1D(nr=100, dr=0.01, origin=0.0)
V = CellVariable(name='V', mesh=mesh, hasOld = True)
r = mesh.cellCenters()

#BC's
V.constrain(0., mesh.facesRight)
V.faceGrad.constrain(0.,mesh.facesLeft)

#initial V profile
V.setValue( 1-r**2)

eqV = 0 == DiffusionTerm( coeff=1., var=V) + 20*0.5*ImplicitSourceTerm(coeff=V*V, var=V)

T = V*V
viewer = Viewer(vars=T)

eqV.solve()

viewer.plot()
raw_input(" Press <enter> to proceed...")

并且获得的配置文件非常相似，但是 y 轴上的值与第一个 FiPy 解决方案或 Mathematica 解决方案不同！扫一扫还是报同样的错误。

Answer 1

我不相信这个问题有除 T = 0 以外的任何解决方案。此外，该解决方案似乎对 cte 的初始条件 and/or 的不同值不稳定。这种不稳定性并不完全令人惊讶，因为 T 形式的方程在 T = 0 时具有无限扩散率。

我怀疑 Mathematica 所做的大致与 FiPy 在您的第一组图中所做的一样，即显示此非线性问题的第一次扫描。这不是答案；只是第一个猜测。不过，我对用 Mma 求解偏微分方程一无所知，无论是在分析上还是在数值上。

你的 V 解决方案后一次扫描后的情节看起来不同，顺便说一句，因为你没有调整初始条件。应该是：

V.setValue( numerix.sqrt(1-r**2))

使用 FiPy 求解圆柱坐标系中的扩散方程不正确

Solution to diffusion equation in cylindrical coordinates using FiPy not correct

python

differential-equations

cylindrical

fipy