Неправильное решение уравнения диффузии в цилиндрических координатах с использованием 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 - справа:

участки

Затем я попытался развернуть, как рекомендуется для такого нелинейного уравнения, как это. Так что вместо 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, Эта нестабильность не совсем удивительна, учитывая, что уравнение в T форма будет иметь бесконечный коэффициент диффузии при T = 0.

Я подозреваю, что Mathematica делает примерно то же, что и FiPy в вашем первом наборе цифр, чтобы показать первый размах этой нелинейной задачи. Это не ответ; только первое предположение. Я ничего не знаю о решении PDE с помощью Mma, ни аналитически, ни численно.

Сюжет после одного цикла после вашего V Решение выглядит иначе, кстати, потому что вы не изменили начальное условие. Так должно быть:

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