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

Question

使用FiPy的圆柱形几何中扩散方程的稳态解与从另一个软件(例如Mathematica)获得的解截然不同。

方程式为：

$0=\frac｛1｝｛r｝\frac{d｝{dr｝\left(\frac｝{T^｛1/2｝｝\frag｛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给出的解决方案是右边的：

绘制

然后，我尝试按照这样的非线性方程的建议进行扫描。所以我没有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｝{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的不同值，该解似乎是不稳定的。考虑到T形式的方程在T=0时具有无限的扩散率，这种不稳定性并不完全令人惊讶。

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

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

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