FIPY 求解 PDE 和 ODE

问题描述

我有一组耦合方程，其中主偏微分方程（时间和位置 z 的函数）给出为：

第二个方程是这样的：

其中 k_m = f(q) 和 q^* = f(c)。如您所见，第二个方程是 ODE（q 不直接依赖于空间）。我发现很难编写代码来耦合这两个方程。到目前为止，对于我忽略第二个方程并采用：q = A*c，其中 A 是某个常数的简单情况，我能够简化并求解以下对流扩散方程：

使用以下代码：

from fipy import Variable,FaceVariable,CellVariable,Grid1D,ExplicitDiffusionTerm,TransientTerm,DiffusionTerm,Viewer,AdvectionTerm,PowerLawConvectionTerm,VanLeerConvectionTerm
from fipy.tools import numerix

#define the grid
L = 3.
nx = L * 512
dx = L/nx
mesh = Grid1D(dx=dx,nx=nx)

# create the variable and initiate it's value on the mesh
conc = CellVariable(name="Conc",mesh=mesh,value=0.)

# physical parameters
Dapp = 1e-7 
u = 0.1
A = 0.85
e = 0.4
F = (1-e)/e

# provide the simplified coefficients
DiffCoeff = Dapp/(1+A*F)
ConvCoeff = ((u/(1+A*F)),)

#Boundary conditions
valueLeft = 1
valueRight = 0.
conc.constrain(valueLeft,mesh.facesLeft)
conc.faceGrad.constrain(valueRight,where=mesh.facesRight)

# define the equation
eqX = TransientTerm() == (DiffusionTerm(coeff=DiffCoeff) - VanLeerConvectionTerm(coeff=ConvCoeff)) 

# time stepping parameters
timestepDuration = 0.001
steps = 50000

from tqdm import tqdm
for step in tqdm(range(steps),desc="Iterating..."):
    eqX.solve(var=conc,dt=timestepDuration)
    # plot every 5000 iterations
    if step%5000 == 0: 
        viewer.plot()

有人可以帮助将对流扩散方程与 fipy 框架中的 ODE 耦合。我对如何取右手边有点困惑，在有限体积的意义上应该只是一个源术语。

_{（https://www.codecogs.com/latex/eqneditor.php 用于生成 Latex 方程）}

解决方法

我不清楚你的两个方程中 q 和 q* 之间的区别是什么。

对于未知数 q、q*、c 和 u，您有两个方程。我不知道 F 和 D 采取什么形式。

在我看来，你的公式只有一两个等式。

您能说出变量的含义吗（例如 u == 速度、c == 浓度、q == 热通量）？如果你说出你想解决的物理问题，也许答案会更容易。

,

这是你的问题，用一些可能会给你线索的东西重新处理。

from fipy import Variable,FaceVariable,CellVariable,Grid1D,ExplicitDiffusionTerm,TransientTerm,DiffusionTerm,Viewer,AdvectionTerm,PowerLawConvectionTerm,VanLeerConvectionTerm
from fipy.tools import numerix

#define the grid
L = 3.
nx = L * 512
dx = L/nx
mesh = Grid1D(dx=dx,nx=nx)


# create the variable and initiate it's value on the mesh
conc = CellVariable(name="Conc",mesh=mesh,value=0.)

# physical parameters
Dapp = 1e-7
u = 0.1
A = 0.85
e = 0.4
F = (1-e)/e

# provide the simplified coefficients
DiffCoeff = Dapp
ConvCoeff = ((u),)

#Boundary conditions
valueLeft = 1
valueRight = 0.
conc.constrain(valueLeft,mesh.facesLeft)
conc.faceGrad.constrain(valueRight,where=mesh.facesRight)

def q_star_func(conc):
    return some_conc_expression

q_star = Variable(q_star_func(conc))
q_star_old = Variable(q_star_func(conc))

# define the equation
eqX = TransientTerm() + F * (q_star - q_star_old) / dt == (DiffusionTerm(coeff=DiffCoeff) - VanLeerConvectionTerm(coeff=ConvCoeff))

# time stepping parameters
timeStepDuration = 0.001
steps = 50000

q_value = 1.

def k_m(q):
    return some_q_expression

viewer = Viewer(conc)
from tqdm import tqdm
for step in tqdm(range(steps),desc="Iterating..."):
    q_star.set_value(q_star_func(conc))
    eqX.solve(var=conc,dt=timeStepDuration)

    q_star_old.set_value(q_star.value)
    q_value = (q_value + timeStepDuration * k_m(q) * q_star.value) / (1 + timeStepDuration * k_m(q))

    # plot every 5000 iterations
    if step%5000 == 0:
        viewer.plot()

您可以为 q_star 的瞬态项添加源项，并将 ODE 添加到循环中。 q_star 必须是一个变量，因为它需要在方程内部更新。 q_value 似乎对 conc 的主要方程没有任何反馈，这有点奇怪。需要定义函数 q_star_func 和 k_m。如果您想要更复杂的方法来求解 q，您也可以使用 Scipy 的 ODE 求解器。您需要通过求解 q_star 的主要方程随时间收集 conc 的值。由于主方程与 q 无关。