问题描述
我想将网格(偏微分方程)与复杂的混合反应器(常微分方程)连接起来。我想取一点值,求解 ODE,然后重新输入网格中的值(将在其中求解 PDE)。
Clang
但是,我看到 .setValue 之后的变量没有用于新的迭代。我知道由于扩散而存在一些差异,但我认为这是不正确的:
# Geometry
Lx = 2 # meters
Ly = 2 # meters
nx = 101 #int((Lx / 0.05) + 1) # nodes
ny = 101 #int((Ly / 0.05) + 1)
# Build the mesh:
mesh = Grid2D(Lx = Lx,Ly = Ly,nx = nx,ny = ny)
X,Y = mesh.cellCenters
xx,yy = mesh.faceCenters
C0_DO = 10. #mg/L # dissolved Oxygen initial concentration
C0_OM = 5.
T0 = 28. #mg/L # # Temperature
print(f'DO saturation concentration: {calc_DOsat(T0)} mg/L')
# dissolved oxygen:
C_DO = CellVariable(name="concentration_DO",mesh=mesh,value=C0_DO,#hasOld=True
)
# Temperature
T = CellVariable(name="temperature",value=T0,#hasOld=True
)
C_OM = CellVariable(name="concentration_OM",value=C0_OM,#hasOld=True
)
# Transport parameters (diffusion constants)
# dissolved oxygen diffusion constant
D_DO = 2. #m2/s
D_DO = FaceVariable(name='DO_diff',value=D_DO) # This is necessary to set a fixed-flux boundary condition
D_DO.constrain(0.,mesh.exteriorFaces) # This is necessary to set a fixed-flux boundary condition
D_OM = 2.
# Reaction and source terms
O_R = 1. # Re-aeration temperature correction factor (-)
Y_R = 1. #reaeration_rate() # Re-aeration rate (1/day)
DOsat = (14.652 - 0.41022 * T + 0.007991 * T ** 2 - 0.000077774 * T ** 3) # Oxygen saturation concentration (mg/L)
xi_SOD = 0.1 # Sediment Oxygen Demand (BC) (1/day)
# Oxygen transport-reaction
eqDO = (TransientTerm(var = C_DO) == # Transient term
DiffusionTerm(coeff = D_DO,var = C_DO) # Diffusion term
+ (mesh.facesTop * (Y_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 #Top Boundary Condition (fixed-flux)
+ (mesh.facesBottom * (xi_SOD * C_DO.faceValue)).divergence # Bottom Boundary Condition (fixed-flux)
)
eqOM = (TransientTerm(var = C_OM) == # Transient term
DiffusionTerm(coeff = D_OM,var = C_OM) # Diffusion term
)
eqTOT = eqDO & eqOM
def reactor_model(C,t):
k1 = 0.5
k2 = 0.1
C1 = C[0]
C2 = C[1]
dC1dt = -k1 * C1
dC2dt = -k1 * C1
dCdt = [dC1dt,dC2dt]
return dCdt
# PDESolver hyperparameters
steps = 10
sweeps = 2
dt = 1e-1
reactor_act = True
t = time.time()
for step in range(steps):
C_DO.updateOld()
C_OM.updateOld()
for sweep in range(sweeps):
eqTOT.sweep(dt=dt)
if reactor_act:
# Reactor
r_in = [C_DO([(Lx/2),(0)]).mean(),# dissolved Oxygen reactor input
C_OM([(Lx/2),(0)]).mean() # Organic Matter reactor input
]
print('Reactor input:',r_in)
# solve ODE
C = odeint(reactor_model,r_in,np.linspace(0,dt,10))
r_out = C[-1]
print('Reactor output:',r_out)
# Reinput ODE
C_DO.setValue(r_out[0],where=(Y >= Ly - mesh.dy / 2 ) & (X <= (Lx/2 + mesh.dx / 2)) & (X > (Lx/2 - mesh.dx / 2)))
C_OM.setValue(r_out[1],where=(Y >= Ly - mesh.dy / 2 ) & (X <= (Lx/2 + mesh.dx / 2)) & (X > (Lx/2 - mesh.dx / 2)))
C_DO.setValue(r_out[0],where=(Y <= mesh.dy / 2 ) & (X <= (Lx/2 + mesh.dx / 2)) & (X > (Lx/2 - mesh.dx / 2)))
C_OM.setValue(r_out[1],where=(Y <= mesh.dy / 2 ) & (X <= (Lx/2 + mesh.dx / 2)) & (X > (Lx/2 - mesh.dx / 2)))
elapsed = time.time() - t
print('Calculation completed')
print(f'Runtime: {elapsed/60} minutes')
有什么方法可以在迭代之间更改 CellValues 吗?
谢谢。
解决方法
暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!
如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@)