使用 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 将单元格约束为零而不是使用边界约束。这也可能会给您带来额外的好处。