Pressure boundary condition in Stokes flow

[RadostW]

My end goal was to compute pressure drops in a geometry similar to example24, but in cylindrical symmetry rather than translational.

I think I understand how most of it works (figuring out that velocity and pressure dofs are just concatenated took me a bit), but I can't figure out how to put pressure boundary conditions rather than velocity boundary conditions.

When boundary conditions on velocity are prescribed everything works fine and the solution is parabolic velocity and linear pressure distribution. When pressure boundary conditions are prescribed, pressure on the boundary is as requested but flow field and pressure distribution is crazy - the flow stops being incompressible 😞 .

How to correctly give pressure boundary conditions in this problem?

Here is my code. I tried my best at making it a minimal example with mesh generation included.

#
# Generate mesh used
#

import numpy as np
import pygmsh
from skfem import MeshTri, Basis, ElementTriP1

floor_depth = 5.0
floor_width = 5.0
mesh_size = 1

box_points = [
    [floor_width, -floor_depth],
    [0, -floor_depth],
    [0, floor_depth],
    [floor_width, floor_depth],
]

# Create the geometry and mesh using pygmsh
with pygmsh.geo.Geometry() as geom:
    poly = geom.add_polygon(box_points, mesh_size=mesh_size)
    raw_mesh = geom.generate_mesh()

# Convert the mesh to a skfem MeshTri object and define boundaries
mesh = MeshTri(
    raw_mesh.points[:, :2].T, raw_mesh.cells_dict["triangle"].T
).with_boundaries(
    {
        "left": lambda x: np.isclose(x[0], 0),
        "right": lambda x: np.isclose(x[0], floor_width),
        "bottom": lambda x: np.isclose(x[1], -floor_depth),
        "top": lambda x: np.isclose(x[1], floor_depth),
    }
)

#
# Take mesh and solve Stokes equation on it
#

from skfem import ElementVector, ElementTriP2, FacetBasis, asm, bmat, solve, condense
from skfem.models.poisson import vector_laplace
from skfem.models.general import divergence

element = {"u": ElementVector(ElementTriP2()), "p": ElementTriP1()}
basis = {variable: Basis(mesh, e, intorder=3) for variable, e in element.items()}

velocity_bc = basis["u"].zeros()
pressure_bc = basis["p"].zeros()

#
# Toggle here to check velocity B.C. works
#
pressure_driven = True

if pressure_driven:
    pressure_bc[basis["p"].get_dofs("top")] = -100
    pressure_bc[basis["p"].get_dofs("bottom")] = 100.0

    velocity_boundary = basis["u"].get_dofs({"left", "right"})
    pressure_boundary = basis["p"].get_dofs({"top", "bottom"})
    combined_boundary = np.hstack([velocity_boundary, basis["u"].N + pressure_boundary])
else:
    inlet_basis = FacetBasis(
        mesh,
        element["u"],
        facets=np.hstack([mesh.boundaries["bottom"], mesh.boundaries["top"]]),
    )

    def parabolic(x):
        """return the plane Poiseuille parabolic inlet profile"""
        return np.stack([np.zeros_like(x[0]), 4 * x[0] * (floor_width - x[0])])

    velocity_bc = inlet_basis.project(parabolic)

    velocity_boundary = basis["u"].get_dofs({"left", "right", "top", "bottom"})
    combined_boundary = np.hstack([velocity_boundary])

boundary_conditions = np.hstack([velocity_bc, pressure_bc])

A = asm(vector_laplace, basis["u"])
B = -asm(divergence, basis["u"], basis["p"])
K = bmat([[A, B.T], [B, None]], "csr")

solution = solve(*condense(K, x=boundary_conditions, D=combined_boundary))
velocity, pressure = np.split(solution, K.blocks)

if __name__ == "__main__":

    mesh.draw(boundaries=True).show()
    basis["p"].plot(pressure, cmap="viridis").show()
    basis["u"].plot(velocity, cmap="viridis").show()

[kinnala]

Setting a pressure BC is not possible in the Stokes variational formulation you are using.

There exist some nonstandard variational formulations for the Stokes problem that allow such pressure conditions: https://www.ljll.fr/hecht/ftp/ff++days/2016/pdf/MS.pdf

I have never tried such formulations.

[RadostW]

Hmm... I just realized that velocity has a lot more degrees of freedom than pressure.

I know that inlet profile in my case is approximately parabolic, so I'll just use some test functions with small wave number and combine it with a least squares solver to guess a smooth velocity profile giving flat pressure distribution.

I'll post more code here if this approach is successful.