I want to solve N-S equations with a model which water flows vertically from high position to low.
So, the inlet boundary is 'mesh.facesBack' with a assigned pressure or velocity, the outlet boundary is 'mesh.facesFront' with zero pressure.
I programme the code based on the example 'examples.flow.stokesCavity'.
The ideal result is the same 'zVelocity' along the vertical direction. However, I always get results which 'zVelocity' vary in the vertical direction and 'zVelocity' near outlet is bigger than it nears the inlet.
Are there something wrong with my fipy codes? Can anyone give me some advice?
from fipy import CellVariable, FaceVariable, Grid3D, DiffusionTerm, Viewer, TransientTerm, ConvectionTerm, TSVViewer
from fipy.tools import numerix
L = 20e-3
N = 5
dL = L / N
viscosity = 1e-6
#zVelocity boudary condition in mesh.facesBack
U = -0.01
#rho is density of water
rho = 1000.0
#0.8 for pressure and 0.5 for velocity are typical relaxation values for SIMPLE
pressureRelaxation = 0.8
velocityRelaxation = 0.5
mesh = Grid3D(nx=N, ny=N,nz=N ,dx=dL, dy=dL,dz=dL)
pressure = CellVariable(mesh=mesh, name='pressure')
pressureCorrection = CellVariable(mesh=mesh)
xVelocity = CellVariable(mesh=mesh, name='X velocity',hasOld=1)
yVelocity = CellVariable(mesh=mesh, name='Y velocity',hasOld=1)
zVelocity = CellVariable(mesh=mesh, name='Z velocity',hasOld=1)
#rank 1 FaceVariable for calculating the mass flux.
velocity = FaceVariable(mesh=mesh, rank=1)
#coeff_con is the coeff of convection term
coeff_con = FaceVariable(mesh=mesh,rank=1)
xVelocityEq = TransientTerm(var=xVelocity) + ConvectionTerm(coeff=coeff_con,var=xVelocity)== \
DiffusionTerm(coeff=viscosity,var=xVelocity) - 1./rho*pressure.grad.dot([1.,0.,0.])
yVelocityEq = TransientTerm(var=yVelocity) + ConvectionTerm(coeff=coeff_con,var=yVelocity)== \
DiffusionTerm(coeff=viscosity,var=yVelocity) - 1./rho*pressure.grad.dot([0.,1.,0.])
zVelocityEq = TransientTerm(var=zVelocity) + ConvectionTerm(coeff=coeff_con,var=zVelocity)== \
DiffusionTerm(coeff=viscosity,var=zVelocity) - 1./rho*pressure.grad.dot([0.,0.,1.])
ap = CellVariable(mesh=mesh, value=1.)
#apa = CellVariable(mesh=mesh,value=1.)
coeff = 1./ ap.arithmeticFaceValue*mesh._faceAreas * mesh._cellDistances
pressureCorrectionEq = DiffusionTerm(coeff=coeff,var=pressureCorrection) - velocity.divergence
from fipy.variables.faceGradVariable import _FaceGradVariable
volume = CellVariable(mesh=mesh, value=mesh.cellVolumes, name='Volume')
contrvolume=volume.arithmeticFaceValue
#slip boundary conditions. The inlet is 'facesBack' with a zVelocity, and the
#outlet is 'facesFront' with zero pressure.
xVelocity.constrain(0., mesh.facesRight | mesh.facesLeft)
yVelocity.constrain(0., mesh.facesTop | mesh.facesBottom)
zVelocity.constrain(U,mesh.facesBack)
pressure.constrain(0.0,mesh.facesFront)
pressureCorrection.constrain(0.0,mesh.facesFront)
##slip boundary conditions. The inlet is 'faceBack' with a assigned pressure,
##and the outlet is 'faceFront' with zero pressure
#xVelocity.constrain(0.,mesh.facesRight | mesh.facesLeft)
#yVelocity.constrain(0.,mesh.facesTop | mesh.facesBottom)
#zVelocity.grad.constrain(0.*mesh.faceNormals, mesh.facesFront | mesh.facesBack)
#pressure.constrain(20,mesh.facesBack)
#pressure.constrain(0.,mesh.facesFront)
#pressureCorrection.constrain(0.,mesh.facesBack | mesh.facesFront)
##no-slip boundary conditions. The inlet is 'facesBack' with a assigned pressure,
##and the outlet is 'faceFront' with zero pressure
#xVelocity.constrain(0.,mesh.facesLeft | mesh.facesRight | mesh.facesTop | mesh.facesBottom)
#yVelocity.constrain(0.,mesh.facesLeft | mesh.facesRight | mesh.facesTop | mesh.facesBottom)
#zVelocity.constrain(0.,mesh.facesLeft | mesh.facesRight | mesh.facesTop | mesh.facesBottom)
#zVelocity.grad.constrain(0.*mesh.faceNormals,mesh.facesFront | mesh.facesBack)
#pressure.constrain(20,mesh.facesBack)
#pressure.constrain(0.,mesh.facesFront)
#pressureCorrection.constrain(0.,mesh.facesBack | mesh.facesFront)
from builtins import range
dt = 1e-3
coupled_vel = xVelocityEq & yVelocityEq & zVelocityEq
#while t<1.0:
for t in range(100):
xVelocity.updateOld()
yVelocity.updateOld()
zVelocity.updateOld()
res = 1.0
rhs = 1.0
while (res>1e-7)|(rhs>1e-10):
#update the convection term coeff
coeff_con.setValue([xVelocity.arithmeticFaceValue.value,yVelocity.arithmeticFaceValue.value,zVelocity.arithmeticFaceValue.value])
coupled_vel.cacheMatrix()
res = coupled_vel.sweep(dt=dt,underRelaxation=velocityRelaxation)
mat = coupled_vel.matrix
## update the ap coefficient from the matrix diagonal
ap[:] = numerix.asarray(mat.takeDiagonal()[0:len(ap)])
#apa[:] = numerix.asarray(mat.numpyArray[0:len(ap)].sum(axis=1).flatten())
## cell pressure gradient
presgrad = pressure.grad
# face pressure gradient
facepresgrad = _FaceGradVariable(pressure)
#update the face velcoities with Rhie-Chow correction
velocity[0] = xVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
(presgrad[0].arithmeticFaceValue - facepresgrad[0])
velocity[1] = yVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
(presgrad[1].arithmeticFaceValue - facepresgrad[1])
velocity[2] = zVelocity.arithmeticFaceValue + contrvolume / ap.arithmeticFaceValue * \
(presgrad[2].arithmeticFaceValue - facepresgrad[2])
velocity[2, mesh.facesBack.value] = U
#velocity[...,mesh.facesLeft.value | mesh.facesRight.value | mesh.facesTop.value | mesh.facesBottom.value] =0.0
## solve the pressure correction equation
pressureCorrectionEq.cacheRHSvector()
pres = pressureCorrectionEq.sweep(var=pressureCorrection)
rhs = pressureCorrectionEq.RHSvector
rhs = max(abs(rhs))
res = max(res,pres)
## update the pressure using the corrected value
pressure.setValue(pressure + pressureRelaxation * pressureCorrection )
## update the velocity using the corrected pressure
xVelocity.setValue(xVelocity - pressureCorrection.grad[0] / ap * mesh.cellVolumes)
yVelocity.setValue(yVelocity - pressureCorrection.grad[1] / ap * mesh.cellVolumes)
zVelocity.setValue(zVelocity - pressureCorrection.grad[2] / ap * mesh.cellVolumes)
print ('t=',t,'res=',res,'rhs=',rhs)
TSVViewer(vars=zVelocity).plot()
