Неправильное решение уравнения диффузии в цилиндрических координатах с использованием FiPy - PullRequest
Новая
       72

Неправильное решение уравнения диффузии в цилиндрических координатах с использованием FiPy

0 голосов
/ 13 февраля 2019

Стационарное решение уравнения диффузии в цилиндрической геометрии с использованием 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!Подметание дает ту же ошибку, что и раньше.

1 Ответ

0 голосов
/ 28 февраля 2019

Я не уверен, что у этой проблемы есть какое-либо решение, кроме T = 0. Кроме того, это решение кажется нестабильным при различных значениях начального условия и / или cte.Эта нестабильность вовсе не удивительна, учитывая, что уравнение в форме T будет иметь бесконечный коэффициент диффузии при T = 0.

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

График после одного цикла после вашего V решения выглядит иначе, потому что вы не настроили начальное условие,Должно быть: 

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