基于FiPy反向求解回转窑颗粒流控制PDE的方法问询

我正在模拟回转窑内部颗粒流的动态行为，该过程的控制PDE是Descoins（2005）推导的动态Saemen方程。原稳态Saemen方程基于从z=L到z=0的倒置域建立，我想知道是否可以通过FiPy从回转窑出口到入口反向求解这个PDE，以及具体的实现方法是什么。

目前我已经参考NIST文档和示例编写了部分代码（尚未完成求解），目标是复现Descoins（2005）中展示的动态模拟效果，现有代码如下：

from fipy import CellVariable, Grid1D, TransientTerm, DiffusionTerm, Viewer, ConvectionTerm, FaceVariable,numerix
from builtins import range

Q = 3.055e-5
R = 0.061
L = 1
ang_theta = 0.645772
ang_beta  = 0.0523599
VeloRot   = 8.333e-2
np = 100

mesh = Grid1D(dx=L / np, nx=np)

H = CellVariable(name="H",
                   mesh=mesh,
                   value=0.001,
                   hasOld=True)

Ca=3*numerix.tan(ang_theta)*Q/(4*pi*VeloRot)
Cb=numerix.tan(ang_beta)/numerix.cos(ang_theta)
arg=R**2-(R-H)**2
dHdzL = Ca*(arg**(-3/2))-Cb

H.constrain(0.025,mesh.facesLeft)

dHdzL_facesRight = dHdzL.faceValue[mesh.facesRight.value]
H.faceGrad.constrain(dHdzL_facesRight, mesh.facesRight)

Fh = (2*H/R)-((H/R)**2)
Ut = 2*pi*VeloRot*R

coeffC = (Ut*numerix.tan(ang_beta)*(numerix.sqrt(Fh))*(numerix.sqrt(1-Fh))/numerix.sin(ang_theta))
coeffD = (Ut*R*numerix.arctan(ang_theta)*((Fh)**(1.5))/3)


eq1 = (TransientTerm(coeff=Fh,var=H) 
       == ConvectionTerm(coeff=coeffC,var=H)
       + DiffusionTerm(coeff=coeffD,var=H))
                                                                     
vi = Viewer(vars = H)

for t in range(500):
    H.updateOld()
    eq1.sweep(var=H,dt=1e-1)
    print(coeffD)
    vi.plot()

核心疑问 #

• 是否可以通过FiPy实现从回转窑出口到入口的反向PDE求解？
• 具体的实现方法是什么？

内容的提问来源于stack exchange，提问作者Barlonacio