使用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给出的解决方案是右边的解决方案:
然后我尝试按照这样的非线性方程推荐扫描。所以,我没有使用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解决方案不同!扫地给出了与以前相同的错误。
我不相信这个问题除了T = 0之外还有任何解决方案。此外,该解决方案似乎对初始条件和/或cte
的不同值不稳定。考虑到当t = 0时,T
形式的方程将具有无限的扩散性,这种不稳定性并不完全令人惊讶。
我怀疑Mathematica在你的第一组数字中大致正在做FiPy正在做的事情,这是为了显示这个非线性问题的第一次扫描。这不是答案;只是第一个猜测。但是,无论是分析还是数字,我都不知道用Mma解决偏微分方程。
在你的V
解决方案之后一次扫描后的情节看起来不同,顺便说一下,因为你没有调整初始条件。它应该是:
V.setValue( numerix.sqrt(1-r**2))