By Evalf, and other Nutils contributors

officialStokesNavier-StokesTaylor-HoodRaviard-Thomas

Solves the lid driven cavity problem for stationary Stokes and Navier-Stokes flow. This benchmark problem consists of a square domain with fixed left, bottom and right boundaries and a top boundary that is moving at unit velocity in positive x-direction. Reference results can be found for instance at https://www.acenumerics.com/the-benchmarks.html.

The general conservation laws are:

Dρ/Dt = 0 (mass)
ρ Du_i/Dt = ∇_j(σ_ij) (momentum)

where we used the material derivative D·/Dt := ∂·/∂t + ∇_j(· u_j). The stress tensor is σ_ij := μ (∇_j(u_i) + ∇_i(u_j) - 2 ∇_k(u_k) δ_ij / δ_nn) - p δ_ij, with pressure p and dynamic viscosity μ, and ρ is the density.

Assuming steady, incompressible flow, we take density to be constant. Further introducing a reference length L and reference velocity U, we make the equations dimensionless by taking spatial coordinates relative to L, velocity relative to U, and pressure relative to ρ U^2. This reduces the conservation laws to:

∇_k(u_k) = 0 (mass)
Du_i/Dt = ∇_j(σ_ij) (momentum)

where the material derivative simplifies to D·/Dt := ∇_j(·) u_j, and the stress tensor becomes σ_ij := (∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij, with Reynolds number Re = ρ U L / μ.

The weak form is obtained by multiplication of a test function and partial integration of the right hand side of the momentum balance.

∀ q: ∫_Ω q ∇_k(u_k) = 0 (mass)
∀ v: ∫_Ω (v_i Du_i/Dt + ∇_j(v_i) σ_ij) = ∫_Γ v_i σ_ij n_j (momentum)

A remaining issue with this system is that its solution is not unique due to its strictly kinematic boundary conditions: since ∫_Ω ∇_k(u_k) = 0 for any u that satisfies the non-penetrating boundary conditions, any pressure space that contains unity results in linear dependence. Furthermore, the strong form shows that constant pressure shifts do not affect the momentum balance. Both issues are solved by arbitrarily removing one of the basis functions.

The normal velocity components at the wall are strongly constrained. The tangential components are either strongly or weakly constrained depending on the strongbc parameter. Since the tangential components at the top are incompatible with the normal components at the left and right boundary, the constraints are constructed in two steps to avoid Gibbs type oscillations, and to make sure that the non-penetrating condition takes precedence.

Depending on the compatible parameter, the script uses either a Taylor-Hood (False) or a Raviart-Thomas (True) discretization. In case of TH, the system is consistently modified by adding ∫_Ω .5 u_i v_i ∇_j(u_j) to yield the skew- symmetric advective term ∫_Ω .5 (v_i ∇_j(u_i) - u_i ∇_j(v_i)) u_j. In case of RT, the discretization guarantees a pointwise divergence-free velocity field, and skew symmetry is implied.