Stokes Problems

Numerical methods for differential equations

Stokes problems

Stokes equation

We focus on 2D problems where u: Ω → ℝ^d, with d = 2, 3. The function y(x, y) = (ux, uy) and the gradient is represented as ∇ u, with the divergence div y = ∂ux/∂x + ∂uy/∂y = tr(∇ y).

The integral of the divergence of σ is given by:

∫ div σ · v dx = ∫∂Ω σ · v ds - ∫Ω ∇ · v dx , σ = Ω → ℝ^{d × d}

Stokes system

The conservation of momentum is expressed as div σ + ∇ g = 0.

We focus on fluids, specifically incompressible Newtonian fluids, where σ = 2μ sym∇ y - pI.

The equation is:

∫Ω (∇ u - pI) · v  dx + ∫Ω f · v dx - ∫∂Ω (∇ u - pI) · v ds = 0

∫Ω (∇ u : ∇ v - p div v + f · v dx = 0

∫Ω div y q dx = 0

Weak formulation

∫Ω y · u : ∇ v - p div u + f v dx = 0

∫Ω div u q dx = 0   ∀ q ∈ H⁰(Ω)

u ∈ H¹(Ω, ℝ^d) and u:Ω → ℝ^d

Introducing notation

U: H₀¹(Ω, ℝ^d)

Q: L²(Ω)

Numerical methods for differential equations

Stokes problems

• Stokes equation: We focus on 2D problems.
• Gradient: …
• Divergence: …
• Stokes system: Conservation of momentum.
• Incompressible Newtonian fluids: Symmetric gradient.

Strong form

This section covers vectorial and scalar equations, homogeneous Dirichlet boundary conditions, test functions, and the bilinear form with known terms.

The problem can be rewritten as:

Find (u,p) ∈ V x Q such that:

a(u,z) + b(u,p) = F(z) ∀z ∈ V + b(z,q) = 0 ∀q ∈ Q

Where:

F: x ∈ V -> F(x) = -∫_Ωf x dx

F(z) = -∫_Ωf x dx

|F(z)| = |∫_Ωf x z dx| ≤ ∫_Ω|f x| dx ≤ ||f||_L²Ω||x||_L²Ω ≤ ||f||_L²Ω||z||_V, ∀z ∈ V

|F(z)| is bounded, so it is a continuous functional.

Continuity

a(u,z) = ∫_Ω ∇u : ∇z dx

|a(u,z)| ≤ k ||u||_V||z||_V

Coercivity

a(u,z) ≥ α ||u||², α > 0

b(u,q) = -∫_Ωdiv u q dx

|b(u,q)| = |∫_Ωdiv q dx| = |∫_Ω|div u q| dx ≤ ||div u||_L²Ω||q||_L²Ω

∫_Ω|∂u_x/∂y|² dx = ∫_Ω|∂u_x/∂x| + |∂u_y/∂y|² + 2∫_Ω∂u_x∂u_y/∂x dy dx ≤ ||∇u||² + ||∂u_y/∂y||²

[Young inequality: 2ab ≤ a²+b²]

||div u||²_L²Ω ≤ 2||u||²_V

|b(u,q)| ≤ 2||div u||α ||q||α → bounded

b(u,q) is continuous.

Stability

z = u, p = q -> a(u,z) + b(u,p) = F(u)

b(u,p) = 0

a(u,u) = F(u)

α ||u||² ≤ a(u,u) = F(u) = ||f||_L²Ω||u|| ²

||u||² ≤ 1/α ||f||_L²Ω||u|| bounded

If f is small, u is of some order of magnitude small → stable

Stability of p?

INF-SUP condition:

^sup_z∈V b(z,q)

INF-SUP condition (Ladyzenskaja-Babuska-Brezzi) depends only on the shape of the domain. Property:

^sup_f∈sup|F|

Total variation diminishing method (TVD)

TV(^q∧_n+1) ≤ TV(^q∧_n)

Can not create new oscillations. For monotone initial data + TVD, we’ll have a monotone solution ∀ time.

(Actually, there are also non-TVD methods keeping monotonicity)

Introduce new ^¯q_iⁿ

Such that no over/under-shoots are introduced.

Min-mod slope: take slope computed with up and down wind, then apply minmod(a,b):

• a: if |a|<|b|, ab > 0
• b: if |b|<|a|, ab > 0
• 0: if ab ≤ 0

^¯q_iⁿ = minmod ( _{qiⁿ-qi-1ⁿ}/dx, _{qi+1ⁿ-qiⁿ}/dx )

Monotonized central differences limiter (2^nd order)

Comparison between up, down, centered.

Burgers equation (inviscid)

∂u/∂t + ∂F(u)/∂x = 0, F(x) = u²/2 ⇒ ∂u/∂t + u ∂u/∂x = 0

The first time at which 2 characteristic lines cross is:

Critical time tc = -Δx/(u₀(x₁+Δx) - u₀(x₁))

Computing lim: tc = min ( Δx/u₀(x) )

Equal area rule to define position of discontinuity:

Ś = _F(u)|_x2-_F(u)|_x1u2-u1 ⇒ for ^Burgers Ś = 1/2(u₂-u₁)

Entropy condition: F'(u₂) ≥ Ś ≥ F'(u₁)

Viscid burgers equation

∂u/∂t + u∂u/∂x = ε ∂²u/∂x² ⇒ u_ε(x,t), u(x,t) = lim u_ε(x,t)

Lax-Friedrichs, Upwind, Lax-Wendroff, Mac-Cormack:

v_j^n' = uⁿ – Δt / 2Δx [ F(u^j+1, n)+F(u^j, n) ]

uⁿ⁺¹_j = 1/2 [ (uⁿ_j+1+uⁿ_j) - Δt / 2Δx [ F(uⁿ_j+1)-F(u_j-1, n) ]

Uniqueness

∫_Ω-Δu + ∇p = f = 0

div u = 0

u = 0

If p ⇒ p + C, then ∇p ⇒ ∇ (p + c) = ∇p

Pressure is not defined in a unique way under Dirichlet boundary conditions.

p̄ ∈ Q: b(v, p̄) = 0 ∀ v ∈ V

∃_S div vₙ p̄ dx = 0 ⇒ ∫_S₂ vₙ p̄ ds - ∫_Ω v ∇p̄ dx = 0

Q₀ = {p̄ ∈ Q: b(v, p̄) = 0 ∀ v ∈ V} ⇒ Solution is unique up to an element of Q₀

For u, we have uniqueness.

Neumann boundary conditions

We provide the Neumann boundary condition: ∇n = 0 on Γ_N ⇒ (∇u - pIₙ) = 0.

It is sufficient to have a piece of Neumann (or Robin) boundary condition to make pressure unique.

Discretization of the problem

V_h ⊂ V, Q_h ⊂ Q, finite dimensional spaces

Find (u_h, p_h) ∈ V_h x Q_h such that:

a(u_h, v_h) + b(v_h, p_h) = F(v_h) ∀ v_h ∈ V_h

b(u_h, q_h) = 0 ∀ q_h ∈ Q_h

Properties

a(u_h, v_h) ≥ 1/_αh |u_h|²_V ∀ v_h ∈ V_h

sup_{vh ∈ Vh} b(v_h, q_h) / |v_h|_Vh ≥ β_h |q_h|_Qh ∀ q_h ∈ Q_h

Everything can be proved exactly as in continuous variables case.

Convergence

||u - u_h|| ≤ (1 + 1/_αh) inf_vh∈Vh ||u - U_h||_Vh + 1/β_h inf_qh∈Qh ||p - q_h||_Qh

Generally, taking u ∈ P^m approximation is such that:

||u - u_h|| ∼ h^m

||p - q_h|| ∼ h^m

||p - q_h|| ∼ h^m-1

||u - u_h|| ∼ h^m-1

∀_{z ∈ Vh}

∀_{qh ∈ Qh}

||p - q_h||₀ = 0

If V_h = V

If β_h = 0

If β_h ∞

V = span(φ_i)

Q_h = span(ψ_i)

ξ_e a (ψ_i, ψ_i) u_i

ξ_e b (ψ_i, ψ_j) u_i = 0

∀_{r ∈ Ωh} p(ϵ) + p₀ could be fixed, set a priori or a posteriori.

Iterative solvers

Initial guess |_ϵ=0 converges in any case g additional cost, fixed to unknown value or for MINRES, GMRES.

Direct solvers

Most of them are implemented with PA=LU, if properly done can solve even a singular matrix e.g. UMFPACK++ (FreeFem++).

Σ = ^T = Σ

If completely lacking a simp CGSpace for ψ, P must be compatible.

Algebraic counterpart of inf-sup

In principle, we could go through a "monolithic way":

A_u & B^TP = φ ⟹ u = -A_u^-1B^T

B(-A_{u^-1B^T ⟹} ⟹

Reduced the solution of the whole problem to only the solution of pressure S: knowing A_u positive definite dense non-singular matrix, must check _ker(B^T) = 0

For M ∈ R^n×m MV = 0, ker(M) is the set of all V satisfying the equation _ker(B^T) = 0

If ker(B^T) = 0, S is non-singular.

Applying B^T to a generic ξ e.g. from an algebraic point of view inf-sup can be checked with ker(B^T) = 0, or saying that _rank(B) = maximum.

Finite element spaces

When solving the Stokes problem, the main issue is using the correct finite element spaces. An alternative is possible: use a couple of fespaces not satisfying the inf-sup + STABILIZATION technique + stabilization parameter from a physical point of view. Introducing _ε s.t.

Not very good, _ε is a critical point!