[Python] Solving a PDE backwards with FiPy

I'm trying to simulate the dynamic behaviour of the granular flow inside a rotary kiln. The PDE describing it was derived by Descoins (2005) and it is expressed as follows:



The original equation (stationary Saemen's equation) was formulated under an inverted domain (from z=L to z=0). Is it possible to solve this PDE from exit to entrance of the rotary kiln using FiPy? And how...?

For now, I wrote the problem based on the documentation and examples from NIST(not solving):

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()


The desired response was demonstrated by Descoins(2005):


Continue reading...