1

I tried to solve 1D coupled PDEs for an advection-diffusion-reaction problem with the Matlab function Pdepe (https://www.mathworks.com/help/matlab/ref/pdepe.html). This function is not working properly in my case of a high advection term as compared to the diffusion term. Therefore, I searched and found this option of using the Python library FiPy to solve my PDEs system. My initial conditions are u1=1 for 4*L/10

My coupled equations are of the following form:

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

I tried to write it by combining the FiPy examples (examples.convection.exponential1DSource.mesh1D, examples.levelSet.advection.mesh1D, examples.cahnHilliard.mesh2DCoupled).

The following lines are not a working example but my first attempt to write the code. This is my first use of FiPy (out of the tests and examples of the documentation), so this might seem to miss the point completely for the regular users.

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()

Thank you for any suggestion or correction. If you happen to know a good tutorial for this specific kind of problem, I would be happy to read it, since I did not find anything better than the examples in the FiPy documentation.

4

1 回答 1

2

几个问题：

  • 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 变量
  • 如果您声明适用于哪个Variablea Term，则必须对所有Terms 执行此操作，因此请编写，例如， 
    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
  • ACellVariable必须用声明hasOld=True才能调用updateOld()，所以 
    u1 = CellVariable(name="u1", mesh=mesh, value=u10, hasOld=True)
u2 = CellVariable(name="u2", mesh=mesh, value=u20, hasOld=True)

似乎有效的完整代码在这里

于 2019-02-14T17:49:21.967 回答