Hello,
i am working on advection dispersion in porous medium. I previously used
fipy for some complex boundary conditions and it worked nicely, but now
i was trying a 2d example with dispersion in two directions which gave
strange results. So i came back to the simplest problem in 1D (a fixed
velocity, a fixed dispersion coefficient) and it does not give the same
result as the analytical solution (also present in the file)
there is a problem of steepness at the front, i tried all convection
terms formulation : they give the same results, i tried also some
different formulation of the solver, always the same results
if somebody can help me
thanks in advance
olivier
# -*- coding: utf-8 -*-
"""
Test flow using fipy
Created on Fri Aug 01 09:38:14 2014
@author: olive
"""
from fipy import *
from fipy.terms.cellTerm import CellTerm
from scipy.special import erf,erfc
dx=1.;nx=50;Lx=dx*nx;x=linspace(dx/2.,Lx-dx/2,nx)
mesh = Grid1D(dx = dx, nx = nx)
K0 = 25 #m/h
H_init = 1. # m
e = 1. # thickness
poro = 0.25
H = CellVariable(name ="Head",mesh=mesh,value=H_init,hasOld=True)
H.constrain(1.5,where=mesh.facesLeft)
H.constrain(1.,where=mesh.facesRight)
K = CellVariable(mesh=mesh,value=K0)
(ImplicitDiffusionTerm(coeff = K.faceValue)).solve(var=H)
V=-K*H.grad/poro
C = CellVariable(mesh=mesh,value=0.,hasOld=True)
C.constrain(1.,where=mesh.facesLeft)
C.faceGrad.constrain(V[0].faceValue-(V[0]*C).faceValue,where=mesh.facesRight)
Dsp = .05 ; #V.mag*.05
eqTrans = TransientTerm() == ImplicitDiffusionTerm(coeff = Dsp) - \
PowerLawConvectionTerm(coeff=V)
tmax=30.
def calc(tmax):
dt=tmax/200
for i in range(200):
#print tstep,i
C.updateOld()
eqTrans.solve(var = C,dt = dt)
c=C.value
return c
c=calc(tmax)
#analytical solution (2nd order term omitted)
v=1.;t=tmax;Ys=4.
C_an = 0.5*erfc((x-v*t)/sqrt(4*Dsp*v*t))
plot(x,c)
plot(x,C_an)_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
