python - Solving a PDE backwards with FiPy - Stack Overflow

admin管理员组

文章数量:1349174

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

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

• python
• solver
• pde
• fipy

Share Improve this question Follow asked Apr 2 at 5:17 BarlonacioBarlonacio 1111 silver badge11 bronze badge New contributor Barlonacio is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct. 1

• The referenced paper is Descoins 2005 – jeguyer Commented 2 days ago

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 1

The code does not run because pi is not defined:

from fipy.tools.numerix import pi

and because coeffC needs to be a vector field, rather than a scalar field:

coeffC = (Ut*numerix.tan(ang_beta)*(numerix.sqrt(Fh))*(numerix.sqrt(1-Fh))/numerix.sin(ang_theta)) * [[1]]

Beyond that, I'm not following what you're asking about the exit to the entrance.

本文标签： pythonSolving a PDE backwards with FiPyStack Overflow