In order to solve Richards equation:
$C(h)\frac{{\partial h}}{{\partial t}} = \frac{\partial }{{\partial z}}\left[
{K(h)\left( {\frac{{\partial h}}{{\partial z}} + 1} \right)} \right] + Q$
I try to define the convection term by K/C/h. Is that a correct expression?
But when I change the function to this form:
Req = (TransientTerm(coeff=1.0) ==
DiffusionTerm(coeff=K(h)/C(h)) +
VanLeerConvectionTerm (coeff=K(h)/C(h)/h))
The FiPy report me mistake.
How can I define the convection term of Richards equation, if not possible,
also let me know.
This is a similar example, which I want to make some change to solve Richards
equation:
The following script was improved by one of the FiPy authors on the
mailing list and now has the proper non-linear behavior, but crashes
when C(h) = 0, during infiltration. I am using the D = K/C
formulation, which blows up when C=0. I will keep working with this,
but I thought I would pass along the more correct script.
Kris
--
from fipy import *
L = 10.0
nx = 75
dx = L/nx
timeStep = dx/10.0
alpha = 0.5
hs = -1.0/alpha
n = 0.2
expt = -(2 + 3*n)
Ks = 10.0
qs = 0.33
qr = 0.13
steps = 150
mesh = Grid1D(dx=dx, nx=nx)
h = CellVariable(name="$\psi$", mesh=mesh, value=-100.0, hasOld=True)
h.constrain(hs, where=mesh.facesLeft)
h.faceGrad.constrain((5*Ks,), where=mesh.facesRight)
h_lo = (h < hs)
h_hi = (h >= hs)
C = (-(qs-qr)*n*(-alpha*h)**(-n)/h) * h_lo
C.name = "C"
K = (Ks*(-alpha*h)**expt) * h_lo + Ks * h_hi
K.name = "K"
Se = (qr + (qs-qr)*(-alpha*h)**(-n)) * h_lo + qs * h_hi
Se.name = "$S_e$"
# Richards' equation (positive up/right)
Req = (TransientTerm(coeff=1.0) ==
DiffusionTerm(coeff=K/C) +
VanLeerConvectionTerm(coeff=[[1.0]]))
viewer = Viewer(vars=Se, datamin=qr, datamax=qs+0.02)
viewer.plot()
for step in range(steps):
h.updateOld()
Req.solve(dt=timeStep,var=h)
res = 100.0
while res > 1.0E-5:
res = Req.sweep(dt=timeStep,var=h)
viewer.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
