如何将平流扩散反应 PDE 与 FiPy 结合

Question

我试图用 Matlab 函数 Pdepe ( https://www.mathworks.com/help/matlab/ref/pdepe.html )解决一维耦合偏微分方程对流扩散反应问题。与扩散项相比，在我的高平流项的情况下，此功能无法正常工作。因此，我搜索并找到了使用Python库FiPy来解决我的PDEs系统的这个选项。我的初始条件是 u1=1 for 4*L/10我的耦合方程具有以下形式：du1/dt = d/dx(D1 * du1/dx) + g * x * du1/dx - mu1 * u1 / (K + u1) * u2du2/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.66L = 10.nx = 1000mu1 = 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/10u20 = 1.u1 = CellVariable(name="u1", mesh=mesh, value=u10)u2 = CellVariable(name="u2", mesh=mesh, value=u20)## Neumann boundary conditionsu1.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)*u2sourceCoeff2 = 1*mu2*u1/(K+u1)*u2eq11 = (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 & eq12eq2 = eq21 & eq22eqn = eq1 & eq2vi = Viewer((u1, u2))for t in range(100):    u1.updateOld()    u2.updateOld()    eqn.solve(dt=1.e-3)    vi.plot()感谢您的任何建议或更正。如果您碰巧知道针对此类特定问题的好教程，我很乐意阅读它，因为我没有找到比 FiPy 文档中的示例更好的内容。

Answer 1

几个问题：

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

代表一维 rank-1 变量

如果你声明哪个VariableaTerm适用，你必须为所有的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必须用 with 声明hasOld=True才能调用updateOld()，所以

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

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