对流扩散反应PDEs与FiPy的耦合

Question

我试图用Matlab函数Pdepe（https://www.mathworks.com/help/matlab/ref/pdepe.html）求解对流扩散反应问题的一维耦合偏微分方程。这个函数在我的高平流项和扩散项的情况下不能正常工作。 因此，我搜索并找到了这个使用Python库FiPy来解决PDEs系统的选项。 我的初始条件是4*L/10的u1=1

我的耦合方程如下：

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

我试着结合一些小例子来写这篇文章(examples.conversion.exponential1DSource.mesh1D，examples.levelSet.advision.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文档中的示例更好的东西。在

Answer 1

几个问题：

• pythonchained comparisons不能在numpy中工作，因此不能在FiPy中工作。所以，写吧

  u10 = (4*L/10 < x) & (x < 6*L/10)
  

  此外，这使得u10成为一个布尔字段，这会混淆FiPy，因此 写 ^{pr2}$ 或者，更好的是，写信

  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.]]
  

  它表示1D rank-1变量
• 如果您声明VariableaTerm适用于哪一个，则必须对所有Term进行声明，因此请写下，例如。，

  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
  

• 为了调用CellVariable，必须用hasOld=True声明updateOld()，因此

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

似乎可以工作的完整代码是here