fipy - 当我求解 NS 方程时，奇怪的速度结果似乎不能满足连续性

Question

我想用一个水从高处垂直流向低处的模型来求解 NS 方程。

因此，入口边界是具有指定压力或速度的“mesh.facesBack”，出口边界是具有零压力的“mesh.facesFront”。

我根据示例“examples.flow.stokesCavity”编写代码。

理想的结果是沿垂直方向的“zVelocity”相同。但是，我总是得到“zVelocity”在垂直方向上变化的结果，并且出口附近的“zVelocity”大于入口附近的结果。

我的fipy代码有问题吗？谁能给我一些建议？

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

标签： fipy