path: root/typst/main.typ

#import "@preview/charged-ieee:0.1.3": ieee


#show: ieee.with(
  title: [Stokes flow - Particle interaction in low Reynolds Number environments],
  abstract: [
    This report describes the key overview of creeping flow, also called Stokes Flow, in the context of particle-particle and particle-fluid
    interactions relevant to porous media behaviour. Beginning from the Navier-Stokes equations, the transition to the Stokes Flow is assumed
    by a low Reynolds number $"Re"$. Solutions and more modern approaches are reviewed, including the lubrication effects and near-contact
    dynamics. The influence of these interactions on clogging and permeability reduction are discussed, followed by a brief overview of
    intended modelling techniques such as the Lattice-Boltzmann-Method(LBM), the Immersed Boundary Method(IBM) and the Homogenized Lattice-Boltzmann
    -Method(HLBM).
  ],
  authors: (
    (
      name: "Claudius Holeksa",
      department: [CSSR],
      organization: [NORCE Research AS],
      location: [Bergen, Norway],
      email: "clho@norceresearch.no"
    ),
  ),
  index-terms: ("Fluid Mechanics", "Review", "Writing Exercise for Fluid Dynamics", "Writing Assignment for Fluid Dynamics"),
  bibliography: bibliography("refs.bib"),
  figure-supplement: [Fig.],
)

= Introduction

Flow in porous subsurface structures often is dominated by extremely low velocity and high viscous fluids with low Reynolds numbers $"Re" << 1$.
In such a condition, the flow enters a state called the Stokes Flow or creeping flow, where inertial forces can be assumed to be zero and viscous
forces dominate the whole connected set.
This kind of flow is found in geological processes, in carbohydrate recovery, filtration systems and $"CO"_2$ storage where the pore-scale motion
of particles and their interaction with the surrounding porous medium strongly influence macroscopic properties such as the permeability of
rock layers.
Undestanding particle transport, bridging and in general clogging phenomena is critical for the effective use and management of subsurface
reservoirs.
The tiniest changes can alter the established path in a porous geometry and transport behaviour over time.
Thus, the behaviour of individual particles in Stokes Flow form the foundation upon which complex processes such as particle accumulation
in porous structures.

The study of the Stokes Flow has a long history, beginning with Stokes' study of flow past a rigid sphere@stokes1851.\
Subsequent work expanded his work to include a multitude of shapes and the interaction of multiple particles, different geometries and
boundary other conditions.
More recently, research was performed in the understanding of creeping flow, the particle transport and the clogging associated with it, and how those microscopic properties impact the macroscopic development. These carry implications which are essential in water treatment, medicine,
carbohydrate recovery and $"CO"_2$ storage. // CITE

This report delves into the central aspects governing the known behaviour of Stokes Flow and the key extension results towards past single and multiple spheres.
Explores the effect of it in porous media and examines the more recent particle bridging behaviour found in porous structures.


== Stokes Flow

// Write about Navier
The governing equations for incompressible flow arise from the Navier-Stokes formulation, which describes the conservation of momentum
and mass in a viscous fluid.
For a fluid with density $rho$, the dynamic viscosity $mu$ and the
velocity field $u(x,t)$ the Navier-Stokes formulation is primarly defined by the equations

#math.equation(
block: true,
$ rho D/(D t)u = ρ((∂)/(∂t)u + (u · ∇)u) = −∇p + μ∇²u + f,\
∇ · u = 0 $
)

The left side represents the material derivative and thus the inerial acceleration, while the right side contains the pressure gradient,
viscous diffusion and potential body forces such as gravity.

While we are often interested in Navier-Stokes flows, a generalized analytical solution does not exist,
most analytical approaches rely on a set of restrictions or assumptions.
One key parameter is the Reynolds number, previously mentioned before, defined as:

#math.equation(
block: true,
$ "Re" = (rho U L ) / mu
$
)

where $U$ is the characteristic velocity and $L$ a characteristic length.
When this number tends toward zero, the flow becomes dominated by viscous effects.
And since we are mostly interested in low velocity and high viscosity cases we are able to
disregard the inertial terms
Which is why we arrive at the Stokes equations

// Stokes
#math.equation(
block: true,
$ μ∇²u + f = ∇p,\
∇ · u = 0 $
)

The Stokes equations are linear, enabling the principle of superposition, which allows for the combination of multiple solutions to construct
complex flow fields.
The original study for simplicity also included the boundary conditions

#math.equation(
block: true,
$ v_n = 0 "at particles boundary",\
v arrow U "as" |x| arrow infinity
$
)

These conditions allowed for an principal analysis of the effects a particle had in a field.
In the next section we will specialize this for a sphere.

== Solid-Fluid single sphere interaction with Stokes flow

With the previous assumptions a solid sphere of radius $r$ moving with a relative velocity $v$ in an unbounded creeping flow,
we receive the drag force provided by Stokes@kundu2015fluid

#math.equation(
block: true,
$ F_S = 6 pi mu r v $
)

with our boundary conditions specified as

#math.equation(
block: true,
$ v_n = 0 "at" a = r,\
v arrow U "as" a arrow infinity
$
)

where according to Proudman and Pearson@proudman_pearson_1957 at a Reynolds number of 0.05 the predicted drag is two percent lower
than the possibly more correct value by Proudman and Pearson@proudman_pearson_1957.

#math.equation(
block: true,
$ F = F_S ( 1 + (3/8)"Re" + (9/40)("Re")^2 *\
(log "Re" + gamma + (5/3) log 2 -\
(323/360)) + (27/80)*"Re"^3 log "Re" ) $
)

An simpler version of this is Osen correction@proudman_pearson_1957,@kundu2015fluid which is described by

#math.equation(
block: true,
$ F = F_S ( 1 + (3/8)"Re") $
)


//Using spherical coordinates $(r,theta,phi.alt)$ and no external force we receive the velocity components

//#math.equation(
//block: true,
//$ v_r = U cos theta (1- (3r)/(2a) + r^3 / (2a^3)) $
//$ v_theta = - U cos theta (1- (3r)/(4a) + r^3 / (4a^3)) $
//)

== Lubrication Forces for Near Contact

When spheres or a sphere near a straight wall approach each other with a small separation $h<<r$ the flow in the gap is dominated by
lubrication forces. For a sphere approaching a straight wall with velocity $U$ normal to the wall we
have 

#math.equation(
block: true,
$ F_(L) tilde.eq (6 pi mu r^2 U) / h (1-(9h)/(16r) + ...) $
)


With two spheres of equal radius r and otherwise the same conditions approaching

#math.equation(
block: true,
$ F_(L) tilde.eq (6 pi mu r^2 U) / h (1-(h)/(5r) + ...) $
)

These forces diverge the smaller $h$ becomes.
This is due to the viscous stress increasing on small films of
fluid which then in turn dominate the interaction.

== Volume Averaged Navier-Stokes (VANS)

For the suspension of multiple particles Laurez et al@MAYA2024 suggests the use of a volume averaged Navier-Stokes equation where
the a porostiy field $epsilon$ is introduced and defined as

#math.equation(
  block:true,
$ epsilon = 
{
  cases(
    1",      if the cell is fully fluid",
    (0,1)", if the cell fluid-solid mix",
    0",      if the cell is fully solid."
  )
 $
)

The mass balance equation for the fluid phase then reads

#math.equation(
block:true,
$
partial (epsilon rho_f) / (partial t) + nabla ( epsilon rho_f v^f) = 0
$
)
with a fluid momentum balance equation of the form

#math.equation(
block:true,
$
partial(epsilon rho_f v^f) / (partial t) + nabla (epsilon rho_f v^f v^f ) = - epsilon rho_f g + epsilon nabla\
(mu_f (nabla v^f + (nabla v^f)^T) - epsilon^2 mu_f ( v^f - #overline[v]^p) / K). $
)

where $p$ is the pressure, $mu_f$ is the fluid dynamic viscosity, $g$ is the gravitational acceleration, and $#overline[v]^p$ is the averaged
particle velocity on the grid with K being the cell-permeability with distinctions in unresolved and resolved particles@MAYA2024.


== Clogging of porous structures

Modifying a poiseulle channel flow by narrowing the channel we can observe mainly three different ways how the porosity is reduced.
By sieving, bridging and aggregation of surface deposition.
The least interesting phenomenon compared to the other two is sieving which results purely from a single particle having
too large of a volume to fit into the
pore. While bridiging occurs when multiple particles arrive at a passage at the exact time, forming a stable arch which blocks the flow.
Aggregation involves the deposition of particles gradually narrowing the channel.
And while sieving and bridging results from purely mechanical forces aggregation results from electrostatic interactions between particles in the fluid and surface boundaries at micrometer scales. In Laurez study@laurez2025bridging this
effect is disregarded due to the scale of the problem and accumulation of particles occuring gradually.

#figure(
image("narrow_channel.png"),
caption: [One example geometry from Laurez@laurez2025bridging describing the general narrow setup]
)

With the narrow throat diameter $D$ and the pore radius $d_p$ the clogging probability starts to occur spuriously around
$2.43 < D/d_p < 5.26$
while higher ratios allow for the unobstructed passage of particles, lower ratio begin to form stable arches
and ratios below $D/d_p < 1$ cause immediate sieving.

#figure(
  image("bridging.png"),
  caption: [Particles forming an arch visualized by Laurez@laurez2025bridging]
)
#figure(
  image("prob_bridge.png"),
  caption: [Probability 𝑃 (𝑠) of clogging after s escape events for (a) 5, (b) 50, (c) 500, (d) 5000 particles in round constriction pore geometry. @laurez2025bridging]
)

For the aggregational aspects other forces are significant. Derjaguin-Landau-Verwey-Overbeek (DLVO) theory is one which unifies
Van-Der-Waal forces and electromagnetic repulsion where the final particle-particle force reads as

#math.equation(
block:true,
$ F_"ij"^"DLVO"_"pp" = ( - A_"iLj" r_"ij" / (6h^2_"ij") + 64 pi epsilon_0 epsilon_r r_"ij" kappa (k_B T / (Z e)) ^2\
tanh ( Z e psi / (4k_B T) ) ^2 e^(kappa h_"ij") ) n_"ij"
$
)

with $lambda$ being the characteristic wavelength of the interaction, $k_B$ the Boltzmann constant, $T$ the absolute temperature, "Z" is the valence
of the electrolyte, "e" is the electron charge, $A_"iLj"$ the Hamaker constant of the particle, i, the particle j and the liquid medium L,
$R_i$,$R_j$ are the particle radii with $h_"ij" << R_i,R_j$, the equivalent radius $r_"ij" = R_i R_j / (R_i + R_j)$, the dielectric constant
$epsilon_0 epsilon_r$, the surface potential $psi_i$ of the particle, $kappa^(-1) = sqrt((epsilon_0 epsilon_r k_B T ) / (2 dot 10^3 N_A e^2 I_S))$
is the Debye screening length with $I_S$ being the ionic electrolyte strength.
The surface potential is assumed equal $psi_i = psi_j =: psi$.

#math.equation(
block:true,
$ F_"ij"^"DLVO"_"pw" = ( - (A_"iLj" R_"i" ( 1+28 (h_"ij" / lambda))) / (6h^2_"ij" ( 1 + 14 h_"ij" / lambda )^2)\
+2 pi epsilon_0 epsilon_r R_i kappa ((2 psi_i psi_j e^(-kappa h_"ij") - (psi_i^2 + psi_j^2) e^(-2 kappa h_"ij") )/(1-e^(-2 kappa h_"ij")))
) n_"ij"
$
)
where both form a basis to approach aggregation simulation. @MAYA2024
While the potential again contains a singularity with zero distance these can be avoided by shoft shell spring collision models on deeper intersection.

== Discussion

The particles forming arches are often inherently unstable. Slight fluctuations in flow rate, local velocity gradients or the particle shape
can impact the particle arch formation and restore partial flow.
As aggregation progresses, the particle bridge arrive at less unstable states.
Which in turn makes a stochastic analysis of particles in different conditions@laurez2025bridging based on different width to gap ratios interesting.

Laurez@laurez2025bridging explores in detail how probable clogging mechanisms are. While these are quite in depth, opportunities
in relation to multiphase fluid interaction with the solid structures arise as well since his work focuses on a single phase fluid.

Another point of interest is a variation in the inbound boundary pressure condition as Laurez looks into a constant inbound pressure.
Laurez mentions that pressure changes are naturally introduced into pores within complex geometries which makes the observation of a single pore throat with only a constant pressure unsatisfying.
Similar gaps can be seen in the other clogging mechanisms which might provide additional research targets analogous to the gaps in regard to the bridging phenomena.
A comparison with a Lattice-Boltzmann-Method modell is promising.
The required forces within the particles are accounted for now.
Those include Van-Der-Wall forces, electromagnetic repulsion, lubrication, collision, drag, gravity, buoyancy.
Of which gravity, drag, collision and buoyancy where already accounted for.

// Insert png

== Lattice-Boltzmann-Method in an overview

Direct numerical simulation of creeping flow around many interacting particles in realistic porous structures is computationally challenging using
conventional CFD methods.
The Lattice-Boltzmann-Method offers an efficient mesoscopic alternative.
Instead of solving the Navier-Stokes equations directly, LBM is a celullar automaton on an equidistant grid based on a stochastic particle
distribution among discretized velocities, so called descriptors.

#figure(
image("d3q27.jpg"),
caption: [An example of a discrete 27 velocity descriptor.]
)

Through these stochastic behaviours, the macroscopic quantities such as velocity and pressure emerge naturally, recovering the Stokes flow
or Navier-Stokes flow with one caveat, the incompressibility is not fully recovered and is thus weakly compressible.

LBM is particularly suitable for simulating flow in complex porous geometries, because boundary conditions can be implemented directly at the
cell level, making it ideal for digital rock simulations and is often employed in various fields.
The equisdistant grid and the properties of a cellular automaton make it a prime target for acceleration due to the local
dependencies @krueger2017lbm.

To couple a particle dynamics system two major approaches exist within Lattice Boltzmann.
The Immersed Boundary Method(IBM)@Peskin_2002 and the Homogenized Lattice-Boltzmann-Method(HLBM)@KRAUSE2017HLBM.
The IBM approach is managed by calculating the coupling based on Lagrangian nodes from the point of view of the particles,
requiring the recovery of the macroscopic variables and also distributing the force back.
While the approach by HLBM incorporates most the of the coupling into the collision modell.
The key difference between those two being that the coupling in IBM happens from the point of view of the cell while
the exchange in HLBM happens from the cell grids. There are key similarities to VANS due to the identical model of fluid averaging.