问题描述
我正在尝试解决传输反应问题,但根据方法我有不同的解决方案。我认为如果我试图求解耦合方程就会出现问题。
这些是 PDE:
我假设恒定温度(方程中的 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",value=10.,hasOld=True)
# Biochemical Oxygen Demand by nitrogenous Organic Matter
C_NBOD = CellVariable(name="concentration_NBOD",hasOld=True)
# Temperature (constant)
T = CellVariable(name="temperature",value=14.4,hasOld=True)
# Transport parameters
D_DO = FaceVariable(name='DO_diff',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)
取决于我是单独求解三个方程(eqDO.solve(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 的浓度。一起求解方程时,我可以这样定义前面的元素吗?或者,如果我有这些项,最好将方程一一求解?
解决方法
- 如果反应项在它们各自的控制方程之间完全相同,守恒性质可能会更好,但如果您扫描(您应该扫描),这可能无关紧要。
- 您需要sweep。方程之间存在非线性相关性。无论您是作为耦合系统求解还是连续求解方程都没有关系。任何出现在
Term
系数中的东西都必须被扫除。 - 当方程中每个
var=?
中的每个Term
指定相同的变量时,方程之间没有耦合,因此使用&
表示法没有意义。您可能只是使用更多内存来解决三个独立方程的问题。 - 即使您调整了一些
var=?
分配以引入耦合,您通常会因为一开始不耦合而更容易获得解决方案。耦合可以帮助收敛,但它经常对稳定性造成严重破坏。 - 类似地,明智地使用
ImplicitSourceTerm
可以 有助于收敛,但当您试图获得一组方程来求解时,通常只会使事情变得混乱。我会明确地编写这些源代码,例如- CBOD_reaction_coeff * C_CBOD * C_DO
,直到您知道一切正常。 -
描述了对 gradient 的约束,但
(mesh.facesTop * blahblah).divergence
强加了一个边界 flux。对于您的方程,通量为 。您应该使用C_DO.faceGrad.constrain(...)
。