FIPY 求解与 ODE 耦合的 PDE

Question

我有一组耦合方程，其中主 PDE（时间和位置 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 方程）}

Answer 1

这是你的问题，可能会给你一些线索。

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 求解器。您需要通过求解conc 的主要方程来收集q_star 的值。方程不必一起求解，因为主方程独立于q。

Answer 2

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

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

在我看来，你的方程式少了一到两个。

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