Fig 1. Hele-Shaw flow past a circle

Experiment

Simulation

Hele-Shaw flow #

Hele-Shaw flow is characterized as flow through the narrow gap between two parallel plates. The wall-contact induced viscous shear stresses lead to flow patterns that may be described as if the flow were ‘potential’. As the caption of this figure says, this may seem paradoxical: potential flows are those that are free of viscous stresses.

In the next two subsections we describe this effect mathematically. We start with the equations of viscosity dominated incompressible, creeping, flow (the Stokes equations), and show that they reduce to the potential flow equations for an inviscid fluid (the Darcy equations) in the case of fluid forced in between two parallel plates.

Stokes equations: from 3D to 2D #

When an incompressible fluid is characterized by a sufficiently low Reynolds number, the advective effects may be neglected and we retrieve the Stokes equations: $$ \begin{cases} \mu \Delta \mathbf{u} - \nabla p = \mathbf{0} \\ \nabla \cdot \mathbf{u} = 0 \end{cases} $$ where the time-derivative terms is dropped as we assumed that steady state is reached.

The “flattened” appearance of the flow in a narrow slit may lead one to naïvely conclude that Hele-Shaw flow should be modeled as a 2D Stokes goverened system. This is quite untrue: it is precisely in that third dimension that high velocity gradients occur, leading to high viscous drag forces. We thus ought to start our analysis from the three-dimensional Stokes equations: $$ \begin{cases} \mu ( \frac{\partial^2}{\partial x^2} u_x + \frac{\partial^2}{\partial y^2} u_x + \frac{\partial^2}{\partial z^2} u_x ) - \frac{\partial}{\partial x} p = 0 \\ \mu ( \frac{\partial^2}{\partial x^2} u_y + \frac{\partial^2}{\partial y^2} u_y + \frac{\partial^2}{\partial z^2} u_y ) - \frac{\partial}{\partial y} p = 0 \\ \mu ( \frac{\partial^2}{\partial x^2} u_z + \frac{\partial^2}{\partial y^2} u_z + \frac{\partial^2}{\partial z^2} u_z ) - \frac{\partial}{\partial z} p = 0 \\ \frac{\partial}{\partial x}u_x + \frac{\partial}{\partial y}u_y + \frac{\partial}{\partial z} u_z = 0 \end{cases} $$

Focusing first on the velocity dependence on the third dimension, we note that such wall-bounded channel flow has its own name: Poiseuille flow. The well-known corresponding velocity profile is a quadratic function that is zero at the two walls, and achieves its maximum at the center. Schematically, the velocity profile through the gap width looks like:

Making use of this knowlegde, we can write our full three-dimensional velocity field \(\mathbf{u} = \mathbf{u}(x,y,z)\) just in terms the two-dimensional center velocity \(\hat{\mathbf{u}} = (\hat{u}_x , \hat{u}_y , 0 )^\text{T} = \hat{\mathbf{u}}(x,y) \): $$ \mathbf{u}(x,y,z) = (1 - \frac{2}{\delta} z ) \, ( 1 + \frac{2}{\delta} z ) \, \hat{\mathbf{u}}(x,y) \,, $$ where \( \delta \) is the width of the gap. By construction, this assumed flow profile satisfies the no-slip boundary conditions at the two parallel plates: \( \mathbf{u}(x,y,\delta/2) = \mathbf{u}(x,y,-\delta/2) = \mathbf{0} \), and simply equals the center velocity at the midpoint: \( \mathbf{u}(x,y,0) = \hat{\mathbf{u}}(x,y) \).

When this assumed velocity profile is substituted into the three-dimensional Stokes equations, they produce: $$ \begin{cases} \mu ( \frac{\partial^2}{\partial x^2} \hat{u}_x + \frac{\partial^2}{\partial y^2} \hat{u}_x - \frac{8}{\delta^2} \hat{u}_x ) - \frac{\partial}{\partial x} p = 0 \\ \mu ( \frac{\partial^2}{\partial x^2} \hat{u}_y + \frac{\partial^2}{\partial y^2} \hat{u}_y - \frac{8}{\delta^2} \hat{u}_y ) - \frac{\partial}{\partial y} p = 0 \\ \frac{\partial}{\partial z} p = 0 \\ \frac{\partial}{\partial x}\hat{u}_x + \frac{\partial}{\partial y}\hat{u}_y = 0 \end{cases} $$ or, in vector notation: $$ \begin{cases} \mu \Delta \hat{\mathbf{u}} - \nabla \hat{p} = \frac{8 \mu}{\delta^2} \hat{\mathbf{u}} \\ \nabla \cdot \hat{\mathbf{u}} = 0 \end{cases} $$ where \( \hat{p} = \hat{p}(x,y) \), and the Laplace and gradient operators can now be understood as two-dimensional differential operators.

Darcy equations #

As may be observed from this dimensionally-reduced Stokes-like equation, the effect of the viscous wall-shear stress manifests itself as a drag force that scales linearly with the center velocity. In its current form, the equations are called the Darcy-Brinkman equations. However, if \( \mu \Delta \hat{\mathbf{u}} \ll \frac{8 \mu}{\delta^2} \hat{\mathbf{u}} \), as would be the case when \(\delta\) is much smaller than the other geometric features, then the remaining viscous stress may be neglected, and we find: $$ \begin{cases} \hat{\mathbf{u}} = - \frac{\delta^2}{8 \mu} \nabla \hat{p} \\ \nabla \cdot \hat{\mathbf{u}} = 0 \end{cases} $$ which are the Darcy equations.

The Darcy and Stokes equations leads to quite different solution fields. This can be observed in the figures below, for which the same streamline inflow locations have been chosen. In case of Stokes flow, the effect of the obstructing cylinder is noticible much further upstream.

Stokes

Darcy

Simulation #

Finite element formulation #

The focus of the current post is on the Hele-Shaw theory section. As such, the section of the mathematics behind the finite element approximation is kept relatively short, and may seem overly concise to the uninitiated. Other posts will focus more on this aspect.

To obtain a finite element approximation, we require a weak formulation of the partial differential equation at hand. The Sokes and Darcy equations are both systems of differential equations. In weak form, these lead to so-called mixed formulations. For the Darcy equations, we obtain: $$ \text{Find } (\hat{\mathbf{u}}, \hat{p} ) \in \mathbf{H}_{h}(\text{div}) \times L^2 \text{ s.t.}: \qquad\qquad\quad\qquad\quad\qquad\quad\\ \begin{cases} \int \big( \frac{8 \mu}{\delta^2} \hat{\mathbf{u}} \cdot \mathbf{v} - \hat{p} \, \nabla \cdot \mathbf{v} \big) \text{d}S = - \int g \, \mathbf{v}\cdot\mathbf{n} \, \text{d}L \quad \forall\, \mathbf{v} \in \mathbf{H}_0 (\text{div}) \\ \int q\,\nabla \cdot \hat{\mathbf{u}} \text{d}S = 0 \qquad\qquad\qquad\qquad\qquad\quad\,\,\, \forall\, q \in L^2 \end{cases} $$ where \( L^2 \) is the space of square integrable functions and \( \mathbf{H}(\text{div}) \) is the space of vector-valued functions whose divergence’s are functions in \( L^2 \). Members of \( \mathbf{H}(\text{div}) \) permit evaluation of their normal component on domain boundaries, whereby \( \hat{\mathbf{u}}\cdot\mathbf{n} = h \) represents an essential condition which is thus imposed on all functions in the space \( \mathbf{H}_{h} (\text{div}) \) (c.q. \( \mathbf{v}\cdot\mathbf{n} = 0 \) on \( \mathbf{H}_{0} (\text{div}) \)). The natural condition is the imposed pressure \( \hat{p}=g \), which is already substituted in the right-hand-side integral in the above equation. The chosen values for all boundary conditions are illustrated in the figure below.

The bottom of the computational domain is chosen slightly larger than the figure frame, as the streamlines in the original figure slightly curve out of the frame. This figure also shows the relative dimensions used for the computational domain. Units are purposely left out: the qualitative nature of the solutionfield (i.e., the streamslines) is not affected by the scale of the domain, nor by the magnitude of the inflow, or even by the value of the effective permeability \( \frac{8 \mu}{\delta^2} \).

While devising a finite element approximation of the above weak formulation, one must take care to choose an appropriate pair of discrete function spaces. Not every choice of vector-valued velocity and scalar-valued pressure finite element spaces is “stable” and produces accurate results, or even produces a solveable system of equations. Typical compatible spaces for the Darcy equations, are a \( P \)-th order Brezzi-Douglas-Marini (BDM) space for the velocity and a \( ( P-1 ) \)-th discontinuous Galerkin (DG) space for pressure, or \( P \)-th order Raviart-Thomas (RT) space paired with a \( ( P-1 ) \)-th order DG space. In the following implementation, we make use of BDM elements. The mesh size used is illustrated in the figure above.

Code implementation #

This simulation is performed with the finite element method, with the library FEniCS. The following code outputs the required data:

from dolfin import *
from mshr import *

# Physical parameters
D = 1
L = 6
Ht = 1.62 # height till top
Hb = 1.9 # height till bottom

# Numerical parameters
resolution = 120

# Create mesh
domain = Rectangle(Point(-L/2,-Hb),Point(L/2,Ht)) - Circle(Point(0,0), D/2)
mesh = generate_mesh(domain, resolution)

# Create finite element spaces
BDM = FiniteElement("BDM", mesh.ufl_cell(), 1)
DG  = FiniteElement("DG", mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, BDM * DG)

# Define test, trial and solution function
u, p = TrialFunctions(V)
v, q = TestFunctions(V)
u_sol = Function(V)

# Define boundary conditions
def boundary_noflow(x):
    return x[1] < -Hb+DOLFIN_EPS or x[1] > Ht - DOLFIN_EPS \
        or x[0]**2+x[1]**2 < (D/2)**2+DOLFIN_EPS
def boundary_inflow(x):
    return x[0] < -L/2+DOLFIN_EPS
bcs = [DirichletBC(V.sub(0), Constant((0,0)), boundary_noflow),\
      DirichletBC(V.sub(0), Constant((1,0)), boundary_inflow),]

# Weak formulation
a = (dot(u,v) + div(v)*p + div(u)*q)*dx
Lrhs = - Constant(0)*q*dx

# Solve and export
solve(a == Lrhs, u_sol, bcs)
(u,p) = u_sol.split()
u.rename('u','u')
File("sol.pvd") << u

"""
Created on Oct 16 2024

@author: vishal, dokken
 
https://fenicsproject.discourse.group/t/hele-shaw-flow-past-a-circle-in-dolfinx/14014
"""

import numpy as np
import ufl
import gmsh
import dolfinx
from petsc4py import PETSc
from dolfinx import fem, mesh, io
from dolfinx.fem.petsc import LinearProblem
from ufl import Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner
from basix.ufl import element, mixed_element
from mpi4py import MPI
from dolfinx import default_scalar_type, plot
from dolfinx.fem import (functionspace, dirichletbc)
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio, XDMFFile)


# Physical parameters
D = 1
Ht = 1.62 # height till top
Hb = 1.9 # height till bottom

L = 6
H = Ht+Hb
c_x = 0
c_y = 0
r = D/2

# Numerical parameters
resolution = 0.1  # Target characteristic length factor #0.05


model_rank = 0
mesh_comm = MPI.COMM_WORLD

gdim=2

# Create gmsh geometry
gmsh.initialize()

# Create channel
rectangle = gmsh.model.occ.addRectangle(-L/2, -Hb, 0, L, H)

# Create circle
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)


gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()



# # Set mesh resolution
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", resolution)
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")


# Add physical tag for bulk
fluid_marker = 1
volumes = gmsh.model.getEntities(dim=gdim)
assert (len(volumes) == 1)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")

# Save generated Mesh
gmsh.write("mesh2.msh")

# Convert gmsh mesh to dolfinx mesh
mesh_obj, _, _ = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()


# Function spaces for the velocity and for the pressure
k = 1
Q_el = element("BDM", mesh_obj.basix_cell(), k)
P_el = element("DG", mesh_obj.basix_cell(), k - 1)
V_el = mixed_element([Q_el, P_el])
V = functionspace(mesh_obj, V_el)

# Get subspace of V
V0 = V.sub(0)

Q, _ = V0.collapse()



# Define test, trial and solution functshoes that fits feet designion

(u, p) = ufl.TrialFunctions(V)
(v, q) = ufl.TestFunctions(V)


 
# Define boundary conditions

def boundary_noflow(x):    
    walls = np.logical_or(np.isclose(x[1], -Hb), np.isclose(x[1], Ht))
    obstacle = np.isclose((x[0]-c_x)**2 + (x[1]-c_y)**2, r**2)
    return np.logical_or(walls, obstacle)
 
def boundary_inflow(x):
    return np.isclose(x[0], -L/2)

fdim = mesh_obj.topology.dim - 1
facets_noflow = mesh.locate_entities_boundary(mesh_obj, fdim, boundary_noflow)
dofs_noflow = fem.locate_dofs_topological((V0, Q), fdim, facets_noflow)


# Walls

# Define a function for the inflow velocity on the full vector space V
u_noflow_func = fem.Function(Q)

# Define a function to set the inflow profile
def noflow_profile(x):
    values = np.zeros((mesh_obj.geometry.dim, x.shape[1]))
    values[0, :] = 0.0  # Set x-component of the velocity
    return values

# Interpolate the inflow profile into the function
u_noflow_func.interpolate(noflow_profile)

# Apply the Dirichlet BC using the function
bc_noflow = fem.dirichletbc(u_noflow_func, dofs_noflow, V0)

# InFlow
facets_inflow = mesh.locate_entities_boundary(mesh_obj, fdim, boundary_inflow)
dofs_inflow = fem.locate_dofs_topological((V0, Q), fdim, facets_inflow)

# Define a function for the inflow velocity on the full vector space V
u_inflow_func = fem.Function(Q)

# Define a function to set the inflow profile
def inflow_profile(x):
    values = np.zeros((mesh_obj.geometry.dim, x.shape[1]))
    values[0, :] = 1.0  # Set x-component of the velocity
    return values

# Interpolate the inflow profile into the function
u_inflow_func.interpolate(inflow_profile)

# Apply the Dirichlet BC using the function
bc_inflow = fem.dirichletbc(u_inflow_func, dofs_inflow, V0)


bcs = [bc_noflow, bc_inflow]

#Defining the variational problem

f = - fem.Constant(mesh_obj, PETSc.ScalarType(0))


# Define the bilinear form correctly
a = (ufl.inner(u, v) + ufl.div(v)*p + ufl.div(u)*q) * ufl.dx
L = -f*q*ufl.dx

#solving the equation

problem = LinearProblem(
    a,
    L,
    bcs=bcs,
    petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"},
)
try:
    w_h = problem.solve()
except PETSc.Error as e:  # type: ignore
    if e.ierr == 92:
        print("The required PETSc solver/preconditioner is not available. Exiting.")
        print(e)
        exit(0)
    else:
        raise e

sigma_h, u_h = w_h.split()

# Saving the Output

with io.XDMFFile(mesh_obj.comm, "Hele-Shaw_flow/mesh.xdmf", "w") as file:
    file.write_mesh(mesh_obj)
    file.write_function(u_h)


gdim = mesh_obj.geometry.dim
V1 = fem.functionspace(mesh_obj, ("Discontinuous Lagrange", 2, (gdim,)))
u1 = fem.Function(V1, dtype=default_scalar_type)
u1.interpolate(sigma_h)
    
with dolfinx.io.VTXWriter(mesh_obj.comm, "album1.bp", u1, engine="BP4") as vtx:
    vtx.write(0.0)

Visualization #

Visualization is performed in Paraview, with the streamline and tube filters. The seeds of the streamlines are manually placed at the correct inflow points. The tube filter permits a variable tube width, as specified by a minimum value, a maximum scaling w.r.t that minimum value, and a scalar-field between whose minimal and maximal values is interpolated. For the latter, we specify a simple \(\frac{1}{\sqrt{|\mathbf{u}|}}\) function, representing the streching and broadening of the paint streak with increasing resp. decreasing flow velocity. Below follow screenshots of the visualization settings, and some further visualizations of the qualitative behavior of the solution.