Search code examples
pythonfipyfluid-dynamics

Boundary conditions for stokes flow around a sphere using FiPy

I have tried to solve the Stokes flow around a sphere using FiPy. To do that, I chose a cylindrical 2-D mesh (since my problem is axisymmetric). The z-axis passes through the center of the sphere, and the size of the mesh is Lr x Lz. The boundary conditions I have used are shown in the figure below:
enter image description here

I solved the problem above using FiPy library for Python, see the code below. 


    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)

The code converges after 140 iterations. There are many lines in this code (sorry about that), but a great part of them are only to correct the grad method for cylindrical coordinates in Fipy.
Most of the professors with whom I discussed advised me not to set v=v0 at z=Lz (not sure why). Instead, they have suggested me to use Neumann boundary conditions at the exit (i.e., dvr/dz = 0 and dvz/dz = 0). I believe this is the boundary condition by default in FiPy, so all I did was commenting a few lines in my code. 


    #zVelocity.constrain(v0, mesh.facesTop)
    #rVelocity.constrain(0., mesh.facesTop)
    #velocity[0, mesh.facesTop.value] = 0.
    #velocity[1, mesh.facesTop.value] = v0

The problem is that my code no longer converges after commenting these lines. The residual error of the rVelocity equation (rres) goes to 0, and so does the residual error of the pressure correction equation (pres). But the remaining criteria in the while loop (residual error for the zVelocity equation, pressure correction factor and velocity divergence) do not go to 0.
So my question is: why changing the exit condition from (vr=0,vz=v0) to (dvr/dz=0, dvz/dz=0) is causing a convergence issue?

Solution

  • It seems that setting velocity[1, mesh.facesTop.value] = v0 ensures that the inflow and outflows are balanced making it easier to achieve continuity. Now, for this problem,

    https://pages.nist.gov/pfhub/benchmarks/benchmark5-hackathon.ipynb/

    it's suggested that the zero pressure correction value is set near the outlet. Trying that with your code seems to improve things,

    pressureCorrection.constrain(0., mesh.facesTop & (R < dL))

    whilst commenting velocity[1, mesh.facesTop.value] = v0 gets quite low residuals. Also, setting

    pressureCorrection.constrain(0., mesh.facesTop)

    gets even lower residuals, but that might not be physical.

    This fipy code (courtesy of @jeguyer) solves the problem above. It uses a source term to constrain a cell to be zero rather than using a boundary constraint. That might also give you an additional benefit.