我尝试使用 FiPy 解决球体周围的斯托克斯流。为此，我选择了一个圆柱形二维网格（因为我的问题是轴对称的）。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 次迭代后收敛。这段代码中有很多行（对此感到抱歉），但其中很大一部分只是为了纠正 Fipy 中柱坐标的grad方法。
与我讨论过的大多数教授都建议我不要将 v=v0 设置为 z=Lz（不知道为什么）。相反，他们建议我在出口处使用 Neumann 边界条件（即 dvr/dz = 0 和 dvz/dz = 0）。我相信这是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 ) 的残差也变为 0。但是while循环中的其余标准（zVelocity 方程的残余误差、压力校正因子和速度散度）不会变为 0。
所以我的问题是：为什么将退出条件从 ( vr=0 , vz=v0 ) 更改为( dvr/dz =0 , dvz/dz=0 ) 导致收敛问题？