Как связать PDE адвекционной диффузионной реакции с FiPy - PullRequest
Новая
       91

Как связать PDE адвекционной диффузионной реакции с FiPy

0 голосов
/ 30 декабря 2018

Я пытался решить 1D-связанные ФДЭ для задачи адвекции-диффузии-реакции с помощью функции Matlab Pdepe (https://www.mathworks.com/help/matlab/ref/pdepe.html). Эта функция не работает должным образом в моем случае высокого адвекционного члена по сравнению с диффузионным членомПоэтому я искал и нашел эту возможность использования библиотеки Python FiPy для решения моей системы PDE. Мои начальные условия u1 = 1 для 4 * L / 10

Мои связанные уравнения имеют следующую форму:

du1 / dt = d / dx (D1 * du1 / dx) + g * x * du1 / dx - mu1 * u1 / (K + u1) * u2

du2 / dt = d / dx (D2 * du2 / dx) + g * x * du2 / dx + mu2 * u1 / (K + u1) * u2

Я попытался написать это, комбинируя примеры FiPy (examples.convection.exponential1DSource.mesh1D), examples.levelSet.advection.mesh1D, examples.cahnHilliard.mesh2DCoupled).

Следующие строки не рабочий пример, но моя первая попытка написания кода. Это мое первое использование FiPy (изтесты и примеры документации), так что это может показатьсяМазь полностью для постоянных пользователей. 

from fipy import *

g = 0.66
L = 10.
nx = 1000
mu1 = 1.
mu2 = 1.
K = 1.
D1 = 1.
D2 = 1.

mesh = Grid1D(dx=L / 1000, nx=nx)

x = mesh.cellCenters[0]
convCoeff = g*(x-L/2)

u10 = 4*L/10 < x < 6*L/10
u20 = 1.

u1 = CellVariable(name="u1", mesh=mesh, value=u10)
u2 = CellVariable(name="u2", mesh=mesh, value=u20)

## Neumann boundary conditions
u1.faceGrad.constrain(0., where=mesh.facesLeft)
u1.faceGrad.constrain(0., where=mesh.facesRight)
u2.faceGrad.constrain(0., where=mesh.facesLeft)
u2.faceGrad.constrain(0., where=mesh.facesRight)

sourceCoeff1 = -1*mu1*u1/(K+u1)*u2
sourceCoeff2 = 1*mu2*u1/(K+u1)*u2

eq11 = (TransientTerm(var=u1) == DiffusionTerm(coeff=D1, var=u1) + ConvectionTerm(coeff=convCoeff))
eq21 = (TransientTerm(var=u2) == DiffusionTerm(coeff=D2, var=u2) + ConvectionTerm(coeff=convCoeff))

eq12 = ImplicitSourceTerm(coeff=sourceCoeff1, var=u1)
eq22 = ImplicitSourceTerm(coeff=sourceCoeff2, var=u2)

eq1 = eq11 & eq12
eq2 = eq21 & eq22

eqn = eq1 & eq2
vi = Viewer((u1, u2))

for t in range(100):
    u1.updateOld()
    u2.updateOld()
    eqn.solve(dt=1.e-3)
    vi.plot()

Спасибо за любые предложения или исправления.Если вам случится знать хорошее руководство для решения этой конкретной проблемы, я был бы рад прочитать его, поскольку я не нашел ничего лучше, чем примеры в документации FiPy.

1 Ответ

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

Несколько проблем:

  • python цепные сравнения не работают в numpy и, следовательно, не работают в FiPy.Итак, напишите 
    u10 = (4*L/10 < x) & (x < 6*L/10)
    Далее, это делает u10 полем логических значений, что приводит к путанице в FiPy, поэтому пишите 
    u10 = ((4*L/10 < x) & (x < 6*L/10)) * 1.
    или, что еще лучше, пишите 
    u1 = CellVariable(name="u1", mesh=mesh, value=0., hasOld=True)
u2 = CellVariable(name="u2", mesh=mesh, value=1., hasOld=True)
u1.setValue(1., where=(4*L/10 < x) & (x < 6*L/10))
  • ConvectionTerm с векторным коэффициентом.Один из способов получить это - 
    convCoeff = g*(x-L/2) * [[1.]]
    , представляющий переменную 1-го ранга
  • . Если вы объявляете, к какой Variable a Term относится, вы должны сделать это для всех Term s, поэтому напишитеНапример, 
    ConvectionTerm(coeff=convCoeff, var=u1)
  • ConvectionTerm(coeff=g*x, var=u1) не представляет g * x * du1 / dx.Он представляет собой d (g * x * u1) / dx.Итак, я полагаю, вы захотите, чтобы 
    ConvectionTerm(coeff=convCoeff, var=u1) - ImplicitSourceTerm(coeff=g, var=u1)

  • ImplicitSourceTerm(coeff=sourceCoeff1, var=u1 не представляло -1*mu1*u1/(K+u1)*u2, скорее оно представляет -1*mu1*u1/(K+u1)*u2*u1.Итак, для лучшей связи между уравнениями напишите 

    sourceCoeff1 = -mu1*u1/(K+u1)
sourceCoeff2 = mu2*u2/(K+u1)

... ImplicitSourceTerm(coeff=sourceCoeff1, var=u2) ...
... ImplicitSourceTerm(coeff=sourceCoeff2, var=u1) ...

  • Как указано в комментариях @ wd15, вы объявляете четыре уравнения для двух неизвестных.& не означает «сложить два уравнения вместе» (что можно сделать с помощью +), скорее это означает «решить эти два уравнения одновременно».Итак, напишите 

    sourceCoeff1 = mu1*u1/(K+u1)
sourceCoeff2 = mu2*u2/(K+u1)

eq1 = (TransientTerm(var=u1) 
       == DiffusionTerm(coeff=D1, var=u1) 
       + ConvectionTerm(coeff=convCoeff, var=u1) 
       - ImplicitSourceTerm(coeff=g, var=u1) 
       - ImplicitSourceTerm(coeff=sourceCoeff1, var=u2))
eq2 = (TransientTerm(var=u2) 
       == DiffusionTerm(coeff=D2, var=u2) 
       + ConvectionTerm(coeff=convCoeff, var=u2) 
       - ImplicitSourceTerm(coeff=g, var=u2) 
       + ImplicitSourceTerm(coeff=sourceCoeff2, var=u1))

eqn = eq1 & eq2
  • A CellVariable должно быть объявлено с hasOld=True, чтобы вызвать updateOld(), поэтому 
    u1 = CellVariable(name="u1", mesh=mesh, value=u10, hasOld=True)
u2 = CellVariable(name="u2", mesh=mesh, value=u20, hasOld=True)

Полный код, которыйкажется, работает здесь

...