python - How to Couple Advection Diffusion Reaction PDEs with FiPy
问题描述
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.
解决方案
几个问题:
- 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
取一个向量系数。一种方法是
表示一维秩 1 变量convCoeff = g*(x-L/2) * [[1.]]
- 如果您声明适用于哪个
Variable
aTerm
,则必须对所有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)
似乎有效的完整代码在这里
推荐阅读
- javascript - 是否可以使用javascript获取mac地址?
- javascript - 角度文件上传 SyntaxError:位置 1 的 JSON 中的意外数字
- blogdown - blogdown 如何导入之前构建的 Rmd 文件
- c# - 如何在 C# 中使用 linq 动态添加列表中的多个项目
- php - 两张桌子并排
- angular - Angular 8:从本地存储读取,修改并保存回同一个文件
- c# - 保存的搜索以编程方式不返回任何结果
- powerbi - 与外部用户共享工作区、报告和应用程序
- node.js - 如何在nodejs中返回一个带有“for”循环和异步的数组?
- r - 具有不同的 .RMD 文件和输出文件名