Stokes Drag

Problem description

In this example, we solve the classical Stokes drag problem, which models the drag force experienced by a sphere moving through a viscous fluid. The governing equations are the Stokes equations:

\[\begin{align*} -\nabla p + \mu \Delta \mathbf{u} &= 0, \quad && \text{in } \Omega^c, \\ \nabla \cdot \mathbf{u} &= 0, \quad && \text{in } \Omega^c, \end{align*}\]

where:

• $\mathbf{u}$ is the velocity field,
• $p$ is the pressure,
• $\mu$ is the dynamic viscosity,
• $\Omega$ is the sphere, and $\Omega^c = \mathbb{R}^3 \setminus \overline{\Omega}$ is the fluid domain.

The boundary conditions are:

• $\mathbf{u} = \mathbf{U}$ on the sphere's surface, where $\mathbf{U}$ is the velocity of the sphere. This is a no-slip condition.
• $\mathbf{u} \to \mathbf{0}$ at infinity, which means that the fluid is at rest far away from the sphere.

The drag force experienced by the sphere is described by Stokes' law:

\[\mathbf{F}_d = -6\pi\mu R \mathbf{U},\]

where $R$ is the sphere's radius. This drag force, $\mathbf{F}_d$, is the primary quantity of interest in this example. We will compute it using Hebeker's formulation [9], which expresses the velocity field $\mathbf{u}$ as a combination of single- and double-layer potentials:

\[\mathbf{u}(\mathbf{x}) = \mathcal{D}[\boldsymbol{\sigma}](\mathbf{x}) + \eta \mathcal{S}[\boldsymbol{\sigma}](\mathbf{x}),\]

Here, $\boldsymbol{\sigma}$ is the unknown density, $\mathcal{S}$ and $\mathcal{D}$ denote the single- and double-layer potentials, respectively, and $\eta > 0$ is a coupling parameter, which we set to $\eta = \mu$ throughout this example.

Discretization

To discretize the boundary $\Gamma := \partial \Omega$, we employ a second-order triangular mesh created using Gmsh:

using Inti, Gmsh
meshsize = 0.4
R = 2.0
gmsh.initialize()
gmsh.option.setNumber("Mesh.MeshSizeMax", meshsize)
gmsh.model.occ.addSphere(0, 0, 0, R)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(2)
msh = Inti.import_mesh()
gmsh.finalize()

Info    : Meshing 1D...
Info    : [ 40%] Meshing curve 2 (Circle)
Info    : Done meshing 1D (Wall 0.000100689s, CPU 0.000101s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Sphere, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.022492s, CPU 0.02249s)
Info    : 412 nodes 838 elements
Info    : Meshing order 2 (curvilinear on)...
Info    : [  0%] Meshing curve 1 order 2
Info    : [ 30%] Meshing curve 2 order 2
Info    : [ 50%] Meshing curve 3 order 2
Info    : [ 70%] Meshing surface 1 order 2
Info    : [ 90%] Meshing volume 1 order 2
Info    : Surface mesh: worst distortion = 0.976791 (0 elements in ]0, 0.2]); worst gamma = 0.354711
Info    : Done meshing order 2 (Wall 0.00293652s, CPU 0.002936s)

Next we extract the Domain $\Gamma$ from the mesh, and create a Quadrature on it:

Ω = Inti.Domain(e -> Inti.geometric_dimension(e) == 3, Inti.entities(msh)) # the 3D volume
Γ = Inti.boundary(Ω) # its boundary
Γ_msh = view(msh, Γ)
Γ_quad = Inti.Quadrature(Γ_msh; qorder = 2) # quadrature on the boundary

With the quadrature prepared, we can now define the Stokes operator along with its associated integral operators. We use the FMM3D library to accelerate the evaluation of the integral operators:

using FMM3D
# pick a correction and compression method
correction = (method = :adaptive, )
compression = (method = :fmm, )

# define the Stokes operator
μ = η = 2.0
op = Inti.Stokes(; dim = 3, μ)

# assemble integral operators
S, D = Inti.single_double_layer(;
    op,
    target = Γ_quad,
    source = Γ_quad,
    compression,
    correction,
)

(2460×2460 LinearMaps.LinearCombination{StaticArraysCore.SMatrix{3, 3, Float64, 9}} with 2 maps:
  2460×2460 LinearMaps.FunctionMap{StaticArraysCore.SMatrix{3, 3, Float64, 9},true}(#11; issymmetric=false, ishermitian=false, isposdef=false)
  2460×2460 LinearMaps.WrappedMap{StaticArraysCore.SMatrix{3, 3, Float64, 9}} of
    2460×2460 SparseArrays.SparseMatrixCSC{StaticArraysCore.SMatrix{3, 3, Float64, 9}, Int64} with 133944 stored entries, 2460×2460 LinearMaps.LinearCombination{StaticArraysCore.SMatrix{3, 3, Float64, 9}} with 2 maps:
  2460×2460 LinearMaps.FunctionMap{StaticArraysCore.SMatrix{3, 3, Float64, 9},true}(#12; issymmetric=false, ishermitian=false, isposdef=false)
  2460×2460 LinearMaps.WrappedMap{StaticArraysCore.SMatrix{3, 3, Float64, 9}} of
    2460×2460 SparseArrays.SparseMatrixCSC{StaticArraysCore.SMatrix{3, 3, Float64, 9}, Int64} with 133944 stored entries)

Solution and drag force computation

We are now ready to set up and solve the problem. First, we define the boundary conditions (a constant velocity on the sphere):

using StaticArrays
v = 2.0
U = SVector(v,0,0)
f = fill(U, length(Γ_quad))

To solve the linear system, we will use the gmres function from IterativeSolvers. Since the function requires scalar types, we need to convert the vector-valued quantities into scalars and vice versa. We can achieve this by using reinterpret to convert between the vector of SVectors and a vector of Float64s types.

using IterativeSolvers, LinearAlgebra, LinearMaps
T = SVector{3, Float64} # vector type
L = I/2 + D + η * S
L_ = LinearMap{Float64}(3 * size(L, 1)) do y, x
    σ = reinterpret(T, x)
    μ = reinterpret(T, y)
    mul!(μ, L, σ)
    return y
end
σ  = zeros(T, length(Γ_quad))
σ_ = reinterpret(Float64, σ)
f_ = reinterpret(Float64, f)
_, hist = gmres!(σ_, L_, f_; reltol = 1e-8, maxiter = 200, restart = 200, log = true)
hist

Converged after 5 iterations.

Note that gmres converges in very few iterations, highlighting the favorable spectral properties of the Hebeker formulation for this problem.

Drag force computation

Now that we have the density σ, we can compute the drag force. As pointed out in [9, Theorem 2.4], the drag force of the body $\Omega$ is given by:

\[ \mathbf{F}_d = \eta \int_{\Gamma} \boldsymbol{\sigma} \, d\Gamma,\]

which can be approximated using our knowledge of σ and the quadrature Γ_quad:

drag = μ * sum(eachindex(Γ_quad)) do i
    return σ[i] * Γ_quad[i].weight
end

3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 150.7951842149006
  -2.860212291633564e-5
  -1.917912560031302e-5

A quick comparison with the analytical solution indicates a good agreement.

exact = 6π * μ * R * U
relative_error = norm(drag - exact) / norm(exact)
println("Relative error: ", relative_error)

Relative error: 8.379685140274065e-6

The relative error in this example is less than 1e-4, indicating that the numerical solution is very close to the analytical solution.

Visualization

Finally, to visualize the flow field, we need to evaluate our integral representation at points off the boundary. The easiest way to achieve this is to use IntegralPotentials, or the convenient SingleLayerPotential and DoubleLayerPotential wrappers:

𝒮 = Inti.SingleLayerPotential(op, Γ_quad)
𝒟 = Inti.DoubleLayerPotential(op, Γ_quad)
u(x) = 𝒟[σ](x) + η*𝒮[σ](x) - U # fluid velocity relative to the sphere

u (generic function with 1 method)

In the code above, we have created a function u that evaluates the velocity at any point x:

u(SVector(1,2,3))

3-element StaticArraysCore.SVector{3, Float64} with indices SOneTo(3):
 -1.0809566783839382
  0.08181469469358692
  0.12272296468843148

With u defined, we can visualize the flow field around the sphere. For this example we will simply sample points on a grid in the xz plane, and plot the velocity vectors at those points:

using Meshes
using GLMakie
L = 5
targets     = [SVector(x, 0, z) for x in -L:meshsize:L, z in -L:meshsize:L] |> vec
filter!(x -> norm(x) > 1.1 * R, targets) # remove points inside or close to the sphere
directions  = u.(targets)
strength    = norm.(directions)
fig = Figure(size = (1000, 800))
ax  = Axis3(fig[1, 1]; title = "Velocity field", aspect = :data, limits = ([-L, L], [-R, R], [-L, L]))
viz!(msh[Γ], showsegments=true)
arrows!(ax, Point3.(targets), Point3.(directions), arrowsize = 0.15, lengthscale = 0.4, arrowcolor = strength, linecolor = strength)
current_figure()
fig

Summary

This tutorial demonstrates how to solve the Stokes drag problem using the Inti library. The approach combines boundary integral equations with numerical quadrature and iterative solvers to compute the drag force on a sphere in a viscous fluid.