FiPy Transport-reaction 一项取决于两个不同方程上的两个变量

Question

我正在尝试解决传输反应问题，但根据不同的方法，我有不同的解决方案。我认为如果我试图求解耦合方程就会出现问题。

这些是偏微分方程：

我假设温度恒定（方程式中的 T）以及速度场恒定（vx、vy）。

如你所见，反应项中有一个元素取决于两个变量，并且存在于两个不同的变量中（CBOD的降解取决于氧气CDO的浓度，以及氧气的浓度取决于 CBOD 降解的量）。

这是我的代码：

# Geometry
Lx = 2  # meters
Ly = 2  # meters
nx = 41 # nodes
ny = 41 # nodes

# Build the mesh:
mesh = Grid2D(Lx=Lx, Ly = Ly, nx=nx, ny=ny)
X,Y = mesh.cellCenters


# Main variable and initial conditions

#Velocity field (constant):
Vf = CellVariable(name='velocity_field',
                 mesh = mesh,
                 value = [vx, vy])

# Dissolved oxygen concentration:
C_DO = CellVariable(name="concentration_DO", 
                 mesh=mesh, 
                 value=0.,
                 hasOld=True)
C_DO.setValue(9.5, where=(Y >= Ly - 0.5))
C_DO.setValue(9.25, where=(Y < Ly - 0.5) & (Y >= Ly - 1.0))
C_DO.setValue(8.9, where=(Y < Ly - 1.0) & (Y >= Ly - 1.5) )
C_DO.setValue(8.8, where=(Y < Ly - 1.5) & (Y >= Ly - 2.0))

# Biochemical Oxygen Demand by Carbonaceous Organic Matter
C_CBOD = CellVariable(name="concentration_CBOD", 
                 mesh=mesh, 
                 value=10.,
                 hasOld=True)

# Biochemical Oxygen Demand by Nitrogenous Organic Matter
C_NBOD = CellVariable(name="concentration_NBOD", 
                 mesh=mesh, 
                 value=10.,
                 hasOld=True)

# Temperature (constant)
T = CellVariable(name="temperature", 
                 mesh=mesh,
                 value=14.4,
                 hasOld=True)


# Transport parameters
D_DO = FaceVariable(name='DO_diff', mesh=mesh, value=1.)
D_DO.constrain(0., mesh.exteriorFaces)

D_CBOD = 1.
D_NBOD = 1.


## Reaction & source terms

# DO
O_r = 1.025
K_r = 1.

# CBOD:
O_CBOD = 1.047
K_CBOD_0 = 0.2
K_CBOD = K_CBOD_0 / DOsat
CBOD_reaction_coeff = K_CBOD * (O_CBOD ** (T - 20))

# NBOD:
O_NBOD = 1.047
K_NBOD_0 = 0.2
K_NBOD = K_NBOD_0 / DOsat
NBOD_reaction_coeff = K_NBOD * (O_NBOD ** (T - 20))

# Boundary conditions
### fixed flux, atmospheric exchange, included in the main equation.

# Equations definition:

# DO transport-reaction
eqDO = (TransientTerm(var = C_DO) == 
        DiffusionTerm(coeff=D_DO, var = C_DO) 
        - ConvectionTerm(coeff=Vf, var=C_DO)
        + ImplicitSourceTerm(coeff= -1 * CBOD_reaction_coeff * C_CBOD, var=C_DO) 
        + ImplicitSourceTerm(coeff= -1 * NBOD_reaction_coeff * C_NBOD, var=C_DO) 
        + (mesh.facesTop * (K_r * (O_r ** (T.faceValue - 20)) * ((14.652 - 0.41022 * T.faceValue + 0.007991 * T.faceValue ** 2 - 0.000077774 * T.faceValue ** 3) - C_DO.faceValue))).divergence))
    
# CBOD transport-reaction
eqCBOD = (TransientTerm(var = C_CBOD) == 
          DiffusionTerm(coeff=D_CBOD, var = C_CBOD) 
          - ConvectionTerm(coeff=Vf, var=C_CBOD) 
          + ImplicitSourceTerm(coeff= -1 * CBOD_reaction_coeff * C_DO, var=C_CBOD))

# NBOD transport-reaction
eqNBOD = (TransientTerm(var = C_NBOD) == 
          DiffusionTerm(coeff=D_NBOD, var = C_NBOD) 
          - ConvectionTerm(coeff=Vf, var=C_NBOD) 
          + ImplicitSourceTerm(coeff= -1 * NBOD_reaction_coeff * C_DO, var=C_NBOD))

eqQ = eqDO & eqCBOD & eqNBOD

# PDESolver hyperparameters
steps = 230 
dt = 1e-2

for step in range(steps):
  C_DO.updateOld()
  C_CBOD.updateOld()
  C_NBOD.updateOld()

  eqQ.solve(dt=dt)

取决于我是否分别求解三个方程([=30=](dt=dt), eqCBOD.solve(dt=dt), eqNBOD.solve(dt=dt)), 或者在系统中耦合 (eqQ.solve(dt=dt))，我得到不同的结果（网格中的分布相同，但值不同）。我不知道我是否可以在两个不同的方程式中使用不同变量的这个术语；例如：

eqDO = ... + ImplicitSourceTerm(coeff= -1 * CBOD_reaction_coeff * C_CBOD, var=C_DO) <--- Is this OK?
eqCBOD = ... + ImplicitSourceTerm(coeff= -1 * CBOD_reaction_coeff * C_DO, var=C_CBOD) <--- Is this OK?

我想一起求解 CBOD、NBOD 和 DO 的浓度。一起解方程时可以这样定义前面的元素吗？或者如果我有这些项，一个一个地解方程会更好吗？

Answer 1

如果反应项在各自的控制方程之间完全相同，则守恒特性可能更好，但如果您扫过（您应该扫过），则可能无关紧要。
你需要sweep。方程之间有 non-linear 依赖关系。无论您是作为耦合系统求解还是连续求解方程都无关紧要。任何显示为 Term 的系数的东西都必须被清除。
当方程中每个 Term 中的每个 var=? 都指定相同的变量时，方程之间没有耦合，因此使用 & 符号没有意义。您只是使用更多的内存可能会更糟糕地求解三个独立方程。
即使您调整了一些 var=? 分配以引入耦合，您通常会通过一开始不耦合来更容易地获得解决方案。耦合可以帮助收敛，但它通常会对稳定性造成严重破坏。
类似地，明智地使用 ImplicitSourceTerm 可以帮助收敛，但当您试图获得一组方程式来求解时，往往只会混淆事情。我会明确地写下这些来源，例如 - CBOD_reaction_coeff * C_CBOD * C_DO 直到你知道一切正常。
$B(y=L_y) = \frac{\partial C_{DO}}{\partial y}$ 描述了对 gradient 的约束，但是 (mesh.facesTop * blahblah).divergence 强加了边界 flux。对于您的方程式，通量为 $D_{DO} \nabla C_{DO} - C_{DO}\vec{v}$ 。你应该使用 C_DO.faceGrad.constrain(...).

FiPy Transport-reaction 一项取决于两个不同方程上的两个变量

FiPy Transport-reaction with one term depending on two variables on two different equations

fipy