使用 FiPy 围绕球体的斯托克斯流动的边界条件
Boundary conditions for stokes flow around a sphere using FiPy
我曾尝试使用 FiPy 解决围绕球体的 Stokes 流。为此，我选择了一个圆柱形二维网格（因为我的问题是轴对称的）。 z轴通过球心，网格尺寸为Lr x Lz。我用过的边界条件如下图所示：
我使用 Python 的 FiPy 库解决了上面的问题，请参阅下面的代码。
from fipy import *
from fipy.tools import numerix
from fipy.variables.faceGradVariable import _FaceGradVariable
viscosity = 5.55555555556e-06
pfi = 10000. #Penalization for being inside sphere
v0 = 1. #Speed far from sphere
tol = 1.e-6 #Tolerance
Lr=2. #Length of the grid
#No. of cells in the r and z directions
Nr=400
Nz=800
Lz=Lr*Nz/Nr #Height of the grid (=4)
dL=Lr/Nr
mesh = CylindricalGrid2D(nr=Nr, nz=Nz, dr=dL, dz=dL)
R, Z = mesh.faceCenters
r, z = mesh.cellCenters
#Under-relaxation factors
pressureRelaxation = 0.8
velocityRelaxation = 0.5
#Radius of the sphere
rad=0.1
#Distance to the center of the mesh (r=0, z=2)
var1 = DistanceVariable(name='distance to center', mesh=mesh, value=numerix.sqrt(r**2+(z-Lz/2.)**2))
#Pressure and pressure correction variables
pressure = CellVariable(mesh=mesh, value = 0., hasOld=True, name='press')
pressureCorrection = CellVariable(mesh=mesh, value = 0., hasOld=True)
#Cell velocities
zVelocity = CellVariable(mesh=mesh, hasOld=True, name='Z vel')
rVelocity = CellVariable(mesh=mesh,hasOld=True, name='R vel')
#face velocities
velocity = FaceVariable(mesh=mesh, rank=1)
velocityold = FaceVariable(mesh=mesh,rank=1)
#BOUNDARY CONDITIONS (no-flux by default)
zVelocity.constrain(v0, mesh.facesBottom)
zVelocity.constrain(v0, mesh.facesTop)
rVelocity.constrain(0., mesh.facesRight)
rVelocity.constrain(0., mesh.facesBottom)
rVelocity.constrain(0., mesh.facesTop)
pressureCorrection.constrain(0., mesh.facesBottom & (R < dL))
#Penalization term
pi_fi= CellVariable(mesh=mesh, value=0.,name='Penalization term')
pi_fi.setValue(pfi, where=(var1 <= rad) )
rFaces=numerix.array([]) #vertical faces
zFaces=numerix.array([]) #horizontal faces
#Number of cells in each processor
Nr_in_proc = mesh.nx
Nz_in_proc = mesh.ny
for zfcount in range(Nr_in_proc*(1+Nz_in_proc)) :
rFaces=numerix.append(rFaces,[False])
zFaces=numerix.append(zFaces,[True])
for rfcount in range(Nz_in_proc*(1+Nr_in_proc)) :
rFaces=numerix.append(rFaces,[True])
zFaces=numerix.append(zFaces,[False])
#Correct pressure gradient
pressure_correct_grad = CellVariable(mesh=mesh, rank=1)
pressure_correct_grad[0] = pressure.grad[0] - pressure / r
pressure_correct_grad[1] = pressure.grad[1]
#Correct pressure face gradient
pressure_correct_facegrad = FaceVariable(mesh=mesh,rank=1)
pressure_correct_facegrad0 = FaceVariable(mesh=mesh)
pressure_correct_facegrad0.setValue(pressure.faceGrad[0])
pressure_correct_facegrad0.setValue(pressure.faceGrad[0] - pressure.grad[0].arithmeticFaceValue + \
pressure_correct_grad[0].arithmeticFaceValue, where = zFaces)
pressure_correct_facegrad1 = FaceVariable(mesh=mesh)
pressure_correct_facegrad1.setValue(pressure.faceGrad[1])
pressure_correct_facegrad.setValue([pressure_correct_facegrad0.value, pressure_correct_facegrad1.value])
#Correct pressureCorrection gradient
pressureCorrection_correct_grad = CellVariable(mesh=mesh, rank=1)
pressureCorrection_correct_grad[0] = pressureCorrection.grad[0] - pressureCorrection / r
pressureCorrection_correct_grad[1] = pressureCorrection.grad[1]
#Correct pressureCorrection face gradient
pressureCorrection_correct_facegrad = FaceVariable(mesh=mesh,rank=1)
pressureCorrection_correct_facegrad0 = FaceVariable(mesh=mesh)
pressureCorrection_correct_facegrad0.setValue(pressureCorrection.faceGrad[0])
pressureCorrection_correct_facegrad0.setValue(pressureCorrection.faceGrad[0] - pressureCorrection.grad[0].arithmeticFaceValue + \
pressureCorrection_correct_grad[0].arithmeticFaceValue, where = zFaces)
pressureCorrection_correct_facegrad1 = FaceVariable(mesh=mesh)
pressureCorrection_correct_facegrad1.setValue(pressureCorrection.faceGrad[1])
pressureCorrection_correct_facegrad.setValue([pressureCorrection_correct_facegrad0.value, pressureCorrection_correct_facegrad1.value])
coeff = FaceVariable(mesh=mesh,value=1.)
#Navie Stokes equation (no inertia, cylindrical coordinates) + pressure correction equation
rVelocityEq = DiffusionTerm(coeff=viscosity) - pressure_correct_grad.dot([1.,0.]) - ImplicitSourceTerm(pi_fi + viscosity/r**2.)
zVelocityEq = DiffusionTerm(coeff=viscosity) - pressure_correct_grad.dot([0.,1.]) - ImplicitSourceTerm(pi_fi)
pressureCorrectionEq = DiffusionTerm(coeff=coeff) - velocity.divergence
#Matrix for Rhie-Chow interpolation
apr = CellVariable(mesh=mesh, value=1.)
apz = CellVariable(mesh=mesh, value=1.)
ap = FaceVariable(mesh=mesh, value=1.)
volume = CellVariable(mesh=mesh, value=mesh.cellVolumes, name='Volume')
contrvolume = R * dL * dL #Control volume for the faces
sweep=0.
#Residue from sweep methods
rres=1000.
zres=1000.
pres=1000.
cont=1000. #Checks if continuity equation is satisfied
pcorrmax=1000. #Max of pressure correction (from using SIMPLE algorithm)
pressure.updateOld()
pressureCorrection.updateOld()
rVelocity.updateOld()
zVelocity.updateOld()
while (rres > tol or zres > tol or pres > tol or cont > tol or pcorrmax > tol) :
sweep=sweep+1
#Solve the Navier Stokes equations to obtain starred values
rVelocityEq.cacheMatrix()
rres = rVelocityEq.sweep(var=rVelocity,underRelaxation=velocityRelaxation)
rmat = rVelocityEq.matrix
zVelocityEq.cacheMatrix()
zres = zVelocityEq.sweep(var=zVelocity,underRelaxation=velocityRelaxation)
zmat = zVelocityEq.matrix
#Update matrix with diagonal coefficients to be used in Rhie-Chow interpolation
apr[:] = -rmat.takeDiagonal()
apz[:] = -zmat.takeDiagonal()
ap.setValue(apr.arithmeticFaceValue,where=rFaces)
ap.setValue(apz.arithmeticFaceValue,where=zFaces)
#Update the face velocities based on starred values with the Rhie-Chow correction
#Final solution independent of the under-relaxation factor
velocity[0] = (rVelocity.arithmeticFaceValue + (volume / apr * pressure_correct_grad[0]).arithmeticFaceValue - \
contrvolume * (1. / apr).arithmeticFaceValue * pressure_correct_facegrad[0] + (1 - velocityRelaxation) * \
(velocityold[0] - rVelocity.old.arithmeticFaceValue))
velocity[1] = (zVelocity.arithmeticFaceValue + (volume / apz * pressure_correct_grad[1]).arithmeticFaceValue - \
contrvolume * (1. / apz).arithmeticFaceValue * pressure_correct_facegrad[1] + (1 - velocityRelaxation) * \
(velocityold[1] - zVelocity.old.arithmeticFaceValue))
#Boundary conditions (again)
velocity[0, mesh.facesRight.value] = 0.
velocity[0, mesh.facesBottom.value] = 0.
velocity[0, mesh.facesTop.value] = 0.
velocity[1, mesh.facesBottom.value] = v0
velocity[1, mesh.facesTop.value] = v0
#Solve the pressure correction equation
coeff.setValue(contrvolume * (1. / apr).arithmeticFaceValue, where=rFaces)
coeff.setValue(contrvolume * (1. / apz).arithmeticFaceValue, where=zFaces)
pressureCorrectionEq.cacheRHSvector()
pres = pressureCorrectionEq.sweep(var=pressureCorrection)
#Correct pressureCorrection gradient
pressureCorrection_correct_grad[0] = pressureCorrection.grad[0] - pressureCorrection / r
pressureCorrection_correct_grad[1] = pressureCorrection.grad[1]
#Correct pressureCorrection face gradient
pressureCorrection_correct_facegrad0.setValue(pressureCorrection.faceGrad[0])
pressureCorrection_correct_facegrad0.setValue(pressureCorrection.faceGrad[0] - pressureCorrection.grad[0].arithmeticFaceValue + \
pressureCorrection_correct_grad[0].arithmeticFaceValue, where = zFaces)
pressureCorrection_correct_facegrad1.setValue(pressureCorrection.faceGrad[1])
pressureCorrection_correct_facegrad.setValue([pressureCorrection_correct_facegrad0.value, pressureCorrection_correct_facegrad1.value])
#Update the pressure using the corrected value
pressure.setValue(pressure + pressureRelaxation * pressureCorrection )
#Correct pressure gradient
pressure_correct_grad[0] = pressure.grad[0] - pressure / r
pressure_correct_grad[1] = pressure.grad[1]
#Correct pressure face gradient
pressure_correct_facegrad0.setValue(pressure.faceGrad[0])
pressure_correct_facegrad0.setValue(pressure.faceGrad[0] - pressure.grad[0].arithmeticFaceValue + \
pressure_correct_grad[0].arithmeticFaceValue, where = zFaces)
pressure_correct_facegrad1.setValue(pressure.faceGrad[1])
pressure_correct_facegrad.setValue([pressure_correct_facegrad0.value, pressure_correct_facegrad1.value])
#Update the velocity using the corrected pressure
rVelocity.setValue(rVelocity - pressureCorrection_correct_grad[0] / apr * volume)
zVelocity.setValue(zVelocity - pressureCorrection_correct_grad[1] / apz * volume)
velocity[0] = velocity[0] - pressureCorrection_correct_facegrad[0] * contrvolume * (1. / apr).arithmeticFaceValue
velocity[1] = velocity[1] - pressureCorrection_correct_facegrad[1] * contrvolume * (1. / apz).arithmeticFaceValue
#Boundary conditions (again)
velocity[0, mesh.facesRight.value] = 0.
velocity[0, mesh.facesBottom.value] = 0.
velocity[0, mesh.facesTop.value] = 0.
velocity[1, mesh.facesTop.value] = v0
velocity[1, mesh.facesBottom.value] = v0
velocityold[0] = velocity[0]
velocityold[1] = velocity[1]
rVelocity.updateOld()
zVelocity.updateOld()
pcorrmax = max(abs(pressureCorrection.globalValue))
cont = max(abs(velocity.divergence.globalValue))
if sweep % 10 == 0 :
print ('sweep:', sweep,', r residual:',rres, ', z residual',zres, ', p residual:',pres, ', continuity:',cont, 'pcorrmax: ', pcorrmax)
代码在 140 次迭代后收敛。此代码中有很多行（对此感到抱歉），但其中很大一部分仅用于 correct the grad method for cylindrical coordinates in Fipy.
与我讨论过的大多数教授都建议我不要将 v=v0 设置为 z=Lz（不知道为什么）。相反，他们建议我在出口处使用 Neumann 边界条件（即 dvr/dz = 0 和 dvz/dz = 0）。我相信这是 boundary condition by default in FiPy，所以我所做的只是在我的代码中注释几行。
#zVelocity.constrain(v0, mesh.facesTop)
#rVelocity.constrain(0., mesh.facesTop)
#velocity[0, mesh.facesTop.value] = 0.
#velocity[1, mesh.facesTop.value] = v0
问题是我的代码在注释了这些行后不再收敛。 rVelocity 方程的残差 (rres) 去为 0，压力校正方程 (pres) 的残差也是如此。但是 while 循环中的剩余标准（zVelocity 方程的剩余误差、压力校正因子和速度散度）不会变为 0。
所以我的问题是：为什么将退出条件从 (vr=0,vz=v0) 更改为 (dvr/dz=0, dvz/dz=0) 导致收敛问题？
似乎设置velocity[1, mesh.facesTop.value] = v0可以确保流入和流出平衡，从而更容易实现连续性。现在，对于这个问题，
https://pages.nist.gov/pfhub/benchmarks/benchmark5-hackathon.ipynb/
建议零压修正值设置在出口附近。尝试使用您的代码似乎可以改善情况，
pressureCorrection.constrain(0., mesh.facesTop & (R < dL))
而评论 velocity[1, mesh.facesTop.value] = v0 得到的残差非常低。另外，设置
pressureCorrection.constrain(0., mesh.facesTop)
获得更低的残差，但这可能不是物理的。
这 fipy code (courtesy of @jeguyer) solves the problem above. It uses a source term 将单元格约束为零而不是使用边界约束。这也可能会给您带来额外的好处。
我曾尝试使用 FiPy 解决围绕球体的 Stokes 流。为此，我选择了一个圆柱形二维网格（因为我的问题是轴对称的）。 z轴通过球心，网格尺寸为Lr x Lz。我用过的边界条件如下图所示：
我使用 Python 的 FiPy 库解决了上面的问题，请参阅下面的代码。
from fipy import *
from fipy.tools import numerix
from fipy.variables.faceGradVariable import _FaceGradVariable
viscosity = 5.55555555556e-06
pfi = 10000. #Penalization for being inside sphere
v0 = 1. #Speed far from sphere
tol = 1.e-6 #Tolerance
Lr=2. #Length of the grid
#No. of cells in the r and z directions
Nr=400
Nz=800
Lz=Lr*Nz/Nr #Height of the grid (=4)
dL=Lr/Nr
mesh = CylindricalGrid2D(nr=Nr, nz=Nz, dr=dL, dz=dL)
R, Z = mesh.faceCenters
r, z = mesh.cellCenters
#Under-relaxation factors
pressureRelaxation = 0.8
velocityRelaxation = 0.5
#Radius of the sphere
rad=0.1
#Distance to the center of the mesh (r=0, z=2)
var1 = DistanceVariable(name='distance to center', mesh=mesh, value=numerix.sqrt(r**2+(z-Lz/2.)**2))
#Pressure and pressure correction variables
pressure = CellVariable(mesh=mesh, value = 0., hasOld=True, name='press')
pressureCorrection = CellVariable(mesh=mesh, value = 0., hasOld=True)
#Cell velocities
zVelocity = CellVariable(mesh=mesh, hasOld=True, name='Z vel')
rVelocity = CellVariable(mesh=mesh,hasOld=True, name='R vel')
#face velocities
velocity = FaceVariable(mesh=mesh, rank=1)
velocityold = FaceVariable(mesh=mesh,rank=1)
#BOUNDARY CONDITIONS (no-flux by default)
zVelocity.constrain(v0, mesh.facesBottom)
zVelocity.constrain(v0, mesh.facesTop)
rVelocity.constrain(0., mesh.facesRight)
rVelocity.constrain(0., mesh.facesBottom)
rVelocity.constrain(0., mesh.facesTop)
pressureCorrection.constrain(0., mesh.facesBottom & (R < dL))
#Penalization term
pi_fi= CellVariable(mesh=mesh, value=0.,name='Penalization term')
pi_fi.setValue(pfi, where=(var1 <= rad) )
rFaces=numerix.array([]) #vertical faces
zFaces=numerix.array([]) #horizontal faces
#Number of cells in each processor
Nr_in_proc = mesh.nx
Nz_in_proc = mesh.ny
for zfcount in range(Nr_in_proc*(1+Nz_in_proc)) :
rFaces=numerix.append(rFaces,[False])
zFaces=numerix.append(zFaces,[True])
for rfcount in range(Nz_in_proc*(1+Nr_in_proc)) :
rFaces=numerix.append(rFaces,[True])
zFaces=numerix.append(zFaces,[False])
#Correct pressure gradient
pressure_correct_grad = CellVariable(mesh=mesh, rank=1)
pressure_correct_grad[0] = pressure.grad[0] - pressure / r
pressure_correct_grad[1] = pressure.grad[1]
#Correct pressure face gradient
pressure_correct_facegrad = FaceVariable(mesh=mesh,rank=1)
pressure_correct_facegrad0 = FaceVariable(mesh=mesh)
pressure_correct_facegrad0.setValue(pressure.faceGrad[0])
pressure_correct_facegrad0.setValue(pressure.faceGrad[0] - pressure.grad[0].arithmeticFaceValue + \
pressure_correct_grad[0].arithmeticFaceValue, where = zFaces)
pressure_correct_facegrad1 = FaceVariable(mesh=mesh)
pressure_correct_facegrad1.setValue(pressure.faceGrad[1])
pressure_correct_facegrad.setValue([pressure_correct_facegrad0.value, pressure_correct_facegrad1.value])
#Correct pressureCorrection gradient
pressureCorrection_correct_grad = CellVariable(mesh=mesh, rank=1)
pressureCorrection_correct_grad[0] = pressureCorrection.grad[0] - pressureCorrection / r
pressureCorrection_correct_grad[1] = pressureCorrection.grad[1]
#Correct pressureCorrection face gradient
pressureCorrection_correct_facegrad = FaceVariable(mesh=mesh,rank=1)
pressureCorrection_correct_facegrad0 = FaceVariable(mesh=mesh)
pressureCorrection_correct_facegrad0.setValue(pressureCorrection.faceGrad[0])
pressureCorrection_correct_facegrad0.setValue(pressureCorrection.faceGrad[0] - pressureCorrection.grad[0].arithmeticFaceValue + \
pressureCorrection_correct_grad[0].arithmeticFaceValue, where = zFaces)
pressureCorrection_correct_facegrad1 = FaceVariable(mesh=mesh)
pressureCorrection_correct_facegrad1.setValue(pressureCorrection.faceGrad[1])
pressureCorrection_correct_facegrad.setValue([pressureCorrection_correct_facegrad0.value, pressureCorrection_correct_facegrad1.value])
coeff = FaceVariable(mesh=mesh,value=1.)
#Navie Stokes equation (no inertia, cylindrical coordinates) + pressure correction equation
rVelocityEq = DiffusionTerm(coeff=viscosity) - pressure_correct_grad.dot([1.,0.]) - ImplicitSourceTerm(pi_fi + viscosity/r**2.)
zVelocityEq = DiffusionTerm(coeff=viscosity) - pressure_correct_grad.dot([0.,1.]) - ImplicitSourceTerm(pi_fi)
pressureCorrectionEq = DiffusionTerm(coeff=coeff) - velocity.divergence
#Matrix for Rhie-Chow interpolation
apr = CellVariable(mesh=mesh, value=1.)
apz = CellVariable(mesh=mesh, value=1.)
ap = FaceVariable(mesh=mesh, value=1.)
volume = CellVariable(mesh=mesh, value=mesh.cellVolumes, name='Volume')
contrvolume = R * dL * dL #Control volume for the faces
sweep=0.
#Residue from sweep methods
rres=1000.
zres=1000.
pres=1000.
cont=1000. #Checks if continuity equation is satisfied
pcorrmax=1000. #Max of pressure correction (from using SIMPLE algorithm)
pressure.updateOld()
pressureCorrection.updateOld()
rVelocity.updateOld()
zVelocity.updateOld()
while (rres > tol or zres > tol or pres > tol or cont > tol or pcorrmax > tol) :
sweep=sweep+1
#Solve the Navier Stokes equations to obtain starred values
rVelocityEq.cacheMatrix()
rres = rVelocityEq.sweep(var=rVelocity,underRelaxation=velocityRelaxation)
rmat = rVelocityEq.matrix
zVelocityEq.cacheMatrix()
zres = zVelocityEq.sweep(var=zVelocity,underRelaxation=velocityRelaxation)
zmat = zVelocityEq.matrix
#Update matrix with diagonal coefficients to be used in Rhie-Chow interpolation
apr[:] = -rmat.takeDiagonal()
apz[:] = -zmat.takeDiagonal()
ap.setValue(apr.arithmeticFaceValue,where=rFaces)
ap.setValue(apz.arithmeticFaceValue,where=zFaces)
#Update the face velocities based on starred values with the Rhie-Chow correction
#Final solution independent of the under-relaxation factor
velocity[0] = (rVelocity.arithmeticFaceValue + (volume / apr * pressure_correct_grad[0]).arithmeticFaceValue - \
contrvolume * (1. / apr).arithmeticFaceValue * pressure_correct_facegrad[0] + (1 - velocityRelaxation) * \
(velocityold[0] - rVelocity.old.arithmeticFaceValue))
velocity[1] = (zVelocity.arithmeticFaceValue + (volume / apz * pressure_correct_grad[1]).arithmeticFaceValue - \
contrvolume * (1. / apz).arithmeticFaceValue * pressure_correct_facegrad[1] + (1 - velocityRelaxation) * \
(velocityold[1] - zVelocity.old.arithmeticFaceValue))
#Boundary conditions (again)
velocity[0, mesh.facesRight.value] = 0.
velocity[0, mesh.facesBottom.value] = 0.
velocity[0, mesh.facesTop.value] = 0.
velocity[1, mesh.facesBottom.value] = v0
velocity[1, mesh.facesTop.value] = v0
#Solve the pressure correction equation
coeff.setValue(contrvolume * (1. / apr).arithmeticFaceValue, where=rFaces)
coeff.setValue(contrvolume * (1. / apz).arithmeticFaceValue, where=zFaces)
pressureCorrectionEq.cacheRHSvector()
pres = pressureCorrectionEq.sweep(var=pressureCorrection)
#Correct pressureCorrection gradient
pressureCorrection_correct_grad[0] = pressureCorrection.grad[0] - pressureCorrection / r
pressureCorrection_correct_grad[1] = pressureCorrection.grad[1]
#Correct pressureCorrection face gradient
pressureCorrection_correct_facegrad0.setValue(pressureCorrection.faceGrad[0])
pressureCorrection_correct_facegrad0.setValue(pressureCorrection.faceGrad[0] - pressureCorrection.grad[0].arithmeticFaceValue + \
pressureCorrection_correct_grad[0].arithmeticFaceValue, where = zFaces)
pressureCorrection_correct_facegrad1.setValue(pressureCorrection.faceGrad[1])
pressureCorrection_correct_facegrad.setValue([pressureCorrection_correct_facegrad0.value, pressureCorrection_correct_facegrad1.value])
#Update the pressure using the corrected value
pressure.setValue(pressure + pressureRelaxation * pressureCorrection )
#Correct pressure gradient
pressure_correct_grad[0] = pressure.grad[0] - pressure / r
pressure_correct_grad[1] = pressure.grad[1]
#Correct pressure face gradient
pressure_correct_facegrad0.setValue(pressure.faceGrad[0])
pressure_correct_facegrad0.setValue(pressure.faceGrad[0] - pressure.grad[0].arithmeticFaceValue + \
pressure_correct_grad[0].arithmeticFaceValue, where = zFaces)
pressure_correct_facegrad1.setValue(pressure.faceGrad[1])
pressure_correct_facegrad.setValue([pressure_correct_facegrad0.value, pressure_correct_facegrad1.value])
#Update the velocity using the corrected pressure
rVelocity.setValue(rVelocity - pressureCorrection_correct_grad[0] / apr * volume)
zVelocity.setValue(zVelocity - pressureCorrection_correct_grad[1] / apz * volume)
velocity[0] = velocity[0] - pressureCorrection_correct_facegrad[0] * contrvolume * (1. / apr).arithmeticFaceValue
velocity[1] = velocity[1] - pressureCorrection_correct_facegrad[1] * contrvolume * (1. / apz).arithmeticFaceValue
#Boundary conditions (again)
velocity[0, mesh.facesRight.value] = 0.
velocity[0, mesh.facesBottom.value] = 0.
velocity[0, mesh.facesTop.value] = 0.
velocity[1, mesh.facesTop.value] = v0
velocity[1, mesh.facesBottom.value] = v0
velocityold[0] = velocity[0]
velocityold[1] = velocity[1]
rVelocity.updateOld()
zVelocity.updateOld()
pcorrmax = max(abs(pressureCorrection.globalValue))
cont = max(abs(velocity.divergence.globalValue))
if sweep % 10 == 0 :
print ('sweep:', sweep,', r residual:',rres, ', z residual',zres, ', p residual:',pres, ', continuity:',cont, 'pcorrmax: ', pcorrmax)
代码在 140 次迭代后收敛。此代码中有很多行（对此感到抱歉），但其中很大一部分仅用于 correct the grad method for cylindrical coordinates in Fipy.
与我讨论过的大多数教授都建议我不要将 v=v0 设置为 z=Lz（不知道为什么）。相反，他们建议我在出口处使用 Neumann 边界条件（即 dvr/dz = 0 和 dvz/dz = 0）。我相信这是 boundary condition by default in FiPy，所以我所做的只是在我的代码中注释几行。
#zVelocity.constrain(v0, mesh.facesTop)
#rVelocity.constrain(0., mesh.facesTop)
#velocity[0, mesh.facesTop.value] = 0.
#velocity[1, mesh.facesTop.value] = v0
问题是我的代码在注释了这些行后不再收敛。 rVelocity 方程的残差 (rres) 去为 0，压力校正方程 (pres) 的残差也是如此。但是 while 循环中的剩余标准（zVelocity 方程的剩余误差、压力校正因子和速度散度）不会变为 0。
所以我的问题是：为什么将退出条件从 (vr=0,vz=v0) 更改为 (dvr/dz=0, dvz/dz=0) 导致收敛问题？
似乎设置velocity[1, mesh.facesTop.value] = v0可以确保流入和流出平衡，从而更容易实现连续性。现在，对于这个问题，
https://pages.nist.gov/pfhub/benchmarks/benchmark5-hackathon.ipynb/
建议零压修正值设置在出口附近。尝试使用您的代码似乎可以改善情况，
pressureCorrection.constrain(0., mesh.facesTop & (R < dL))
而评论 velocity[1, mesh.facesTop.value] = v0 得到的残差非常低。另外，设置
pressureCorrection.constrain(0., mesh.facesTop)
获得更低的残差，但这可能不是物理的。
这 fipy code (courtesy of @jeguyer) solves the problem above. It uses a source term 将单元格约束为零而不是使用边界约束。这也可能会给您带来额外的好处。