On Wed, Jun 4, 2008 at 9:57 AM, asreeve <[EMAIL PROTECTED]> wrote:
> I'm still having problems with the convection term in my 3-D model. The flow
> conditions should be at steady state, but some cells divergence of the
> fluxes are far enough from zero to cause problems with the solution. If I've
> set the problem up correctly, the max and min concentration values should be
> 1700 and 180, respectively.
I played with the example a little bit. I am not sure why the values
cannot go below 180. After the very first sweep there are values below
180. This is consistent with a flow problem with non-uniform
velocities. In the attached file the diffusion is switched off by
setting the coefficient to 0. The velocity in the x and y directions
has also been set to 0. As the source terms are zero on the first
sweep, the problem reduces to a 1D flow problem. Look at the value
with the lowest value after one sweep, element 538. It has an entry
velocity in of 4.57618430e-05 and a velocity out of 9.72492493e-05.
Thus its value has to go below 180 to satisfy the equations.
Hopefully, the attached file will help explain what I'm trying to
explain. It will stop after one sweep and print a bunch of values.
> I'd note I get these problems what I add heterogeneity to the mesh by
> varying the permeability of cells.
In this case you have a fairly constant velocity in the z-direction
and only diffusion and thus no cells will be driven below 180.
> Things work fine when all the flow is in
> the Z direction (along the axis of the extruded tube), and not diverted in
> the x or y directions (where the grid is not orthogonal).
Not true. The example demonstrates that the values can get below 180
when the velocity in the x and y direction is switched off.
The values in any general flow problem are not necessarily bounded.
Hope this helps.
--
Daniel Wheeler
from fipy import *
def addLayer(verts):
addArray = verts.copy()
addArray[:] = 0
addArray[2] = length/slices
return verts + addArray
length=40. ##cm
slices=40
mesh=GmshImporter2D('./semicircle.msh').extrude(layers=slices, extrudeFunc=addLayer)
delZ=length/slices ##needed later in script
conc=CellVariable(name='solute conc',mesh=mesh,value=180.,hasOld=1)
###calculate head
head=CellVariable(name='head',mesh=mesh,value=110.,hasOld=1.)
Por=CellVariable(name='effective porosity',mesh=mesh,rank=0,value=.2)
SpStore=CellVariable(name='specific storage',mesh=mesh,rank=0,value=0.)
kond=CellVariable(name='conductivity',mesh=mesh,value=1.e-4)
###generate random numbers
#Dist=random.uniform(-2.,2.,mesh.getNumberOfCells())
###Set diffusion coeff
X,Y,Z=mesh.getCellCenters()
kond.setValue(1.e-6, where=((Z>20.)&(Z<30)&(X>2.54)))
#kond.setValue(1.e-4*(2**Dist))
Eq=TransientTerm(SpStore)==DiffusionTerm(kond.getHarmonicFaceValue())
###Set constant head. boundary conditions
X,Y,Z=mesh.getFaceCenters()
ZMax=max(Z)-0.005
bndry=(FixedValue(mesh.getFaces().where((Z<0.005)),120.),
FixedValue(mesh.getFaces().where((Z>ZMax)),100.))
###Calculate steady state distribution
delTime=1.e10
timeSteps=5
for step in range(timeSteps):
head.updateOld()
resid=1.
## print 'head',step,max(head),min(head)
while resid>1.e-14:
resid=Eq.sweep(
var=head,
solver=LinearCGSSolver(iterations=2000,tolerance=1.e-17),
boundaryConditions=bndry,
dt=delTime,
underRelaxation=None)
## print 'head',resid
make(head).plot()
##raw_input('wait')
##for ID in mesh._getCellFaceIDs()[...,846]:
## print ID, head.getFaceGrad()[...,ID]
hfg = head.getFaceGrad()
##print head.getFaceGrad()
##calculate diffusive flux
Flux=FaceVariable(name='specific discharge', mesh=mesh,
value=-hfg*kond.getHarmonicFaceValue())
Flux[0] = 0
Flux[1] = 0
Disp=FaceVariable(name='dispersion tensor',mesh=mesh,rank=1)
Disp.setValue([1.e-5,1.e-5,1.e-5])
Disp.setValue(0)
X,Y,Z=mesh.getFaceCenters()
bndryOn=(FixedValue(mesh.getFaces().where((Z<0.001)),1700.))
X,Y,Z=mesh.getCellCenters()
fluxInOn=CellVariable(name='inflow on',mesh=mesh)
ZValueIn=length/(slices*2.)+1.e-5
fluxInOn.setValue(1700.*Flux.getDivergence(),where=(Z<ZValueIn))
ZValueOut=length-length/(slices*2.)-1.e-5
fluxOut=CellVariable(name='outflow',mesh=mesh,value=0)
fluxOut.setValue(Flux.getDivergence(), where=(Z>ZValueOut))
fluxError=CellVariable(name='water balance error',mesh=mesh,value=0)
fluxError.setValue(Flux.getDivergence(), where=((Z>ZValueIn)&(Z<ZValueOut)))
conc=CellVariable(name='solute conc',mesh=mesh,value=180.,hasOld=1)
transportEqOn=(TransientTerm(coeff=Por)==
DiffusionTerm(coeff=Disp)
-UpwindConvectionTerm(coeff=Flux)
+fluxInOn
+ImplicitSourceTerm(coeff=fluxOut))
## +ImplicitSourceTerm(coeff=fluxError))###should not need this last term
delTime=600.
t=0
timeSteps=100
for step in range(timeSteps):
t+=delTime
print 'conc ',t,max(conc),min(conc)
conc.updateOld()
resid=1.
print step
while resid>1.e-14:
print 'conc before',conc[538]
print 'fluxInOn',fluxInOn[538]
print 'fluxOut',fluxOut[538]
print 'mesh.getCellCenters()[...,538]',mesh.getCellCenters()[...,538]
print 'mesh._getCellFaceIDs()[...,538]',mesh._getCellFaceIDs()[...,538]
print 'mesh.getCellVolumes()[538]',mesh.getCellVolumes()[538]
for ID in mesh._getCellFaceIDs()[...,538]:
print 'ID, Flux[...,ID]',ID, Flux[...,ID]
print 'mesh._getOrientedFaceNormals(...,ID)',mesh._getOrientedFaceNormals()[...,ID]
print 'mesh._getFaceAreas()[...,ID]',mesh._getFaceAreas()[ID]
print 'mesh.getFaceCellIDs()[...,ID]',mesh.getFaceCellIDs()[...,ID]
resid=transportEqOn.sweep(var=conc,
solver=LinearCGSSolver(iterations=500,
tolerance=1.e-16),
boundaryConditions=bndryOn,
dt=delTime)
print 'numerix.argmin(conc)',numerix.argmin(conc)
print 'conc after',max(conc),min(conc)
print 'conc[538]',conc[538]
print 'conc[510]',conc[510]
print 'conc[566]',conc[566]
raw_input("wait")
print 'conc', resid
#viewer=make(conc,limits={'datamin':180.,'datamax':1700.})
#viewer.plot()
