Я попытался решить поток Стокса вокруг сферы, используя FiPy. Для этого я выбрал цилиндрическую двумерную фигуру me sh (поскольку моя проблема - axymymmetri c). Ось z проходит через центр сферы, а размер me sh равен Lr x Lz. Граничные условия, которые я использовал, показаны на рисунке ниже:
Я решил вышеупомянутую проблему, используя библиотеку FiPy для Python, см. Код ниже.
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 итераций. В этом коде много строк (извините за это), но большая часть из них только для исправляет метод grad для цилиндрических координат в Fipy.
Большинство из профессора, с которыми я обсуждал, советовали мне не устанавливать v = v0 при z = Lz (не знаю почему). Вместо этого они предложили мне использовать граничные условия Неймана на выходе (т.е. 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 ). Но остальные критерии в , в то время как l oop (остаточная ошибка для уравнения zVelocity, поправочный коэффициент давления и расхождение скорости), не go до 0.
Так что мой вопрос : почему изменение условия выхода с ( vr = 0 , vz = v0 ) на ( dvr / dz = 0 , dvz / dz = 0 ) вызывает проблему сходимости?