14.1 Fundamental Theories of IncNavierStokesSolver

14.1.1 Governing Equations

A useful tool implemented in Nektar++ is the incompressible Navier Stokes solver that allows one to solve the governing equation for viscous Newtonians fluids governed by

+ u ⋅∇u	= -∇p + ν∇2u + f,	(14.1a)
∇⋅u	= 0,	(14.1b)

where u is the velocity, p is the specific pressure (including density) and ν the kinematic viscosity.

14.1.2 Time discretisation

Velocity-correction scheme

An efficient approach to solve the incompressible navier-stokes equations are splitting schemes. These schemes decouple the monolithic system of the original momentum and mass conversation into subsequent steps that separate solving into a Poisson problem and three velocity problems, assuming a three dimensional problem. One type of the splitting are the velocity-correction type schemes. They begin with an initial step that evaluates the advection term and subsequent steps to solve for the new pressure pⁿ⁺¹ and velocity vector uⁿ⁺¹, where we use the formulation uⁿ = u(tⁿ = nΔt) to denote the discrete time.

Within Nektar++ we consider a rotational velocity-correction scheme VelocityCorrectionScheme that enables higher-order time accuracy despite the splitting [46]. The first step evaluates the advection term N(u)ⁿ = [u ⋅∇u]ⁿ explicitly. Subsequently, the algorithms solves a pressure Poisson problem of the form

∇pn+1 = -- [u ⋅∇u]n
- ν∇×∇×un + fn+1,		(14.2)

where α_q are the coefficients of the backwards difference formula summation that approximates the timer derivative. The equation is supplemented with appropriate pressure boundary conditions to enable higher-order time accuracy [46]. These boundary conditions are

n+1 = n ⋅∑ q=0J-1α qun-q - [u ⋅∇u]n
- ν∇×∇×un + fn+1.		(14.3)

They are implemented within the Extrapolation class. The velocity equation uses the updated pressure pⁿ⁺¹ and solves a Helmholtz problem for each velocity component with

un+1 - ν∇2un+1 =	∑ q=0J-1α qun-q -∇pn+1
	- [u ⋅∇u]n + fn+1.	(14.4)

The right hand side for this equation is computed within ViscousForcing and the solution of the helmholtz problem happens within ViscousSolve via a call to MultiRegions::HelmSolve.

We also consider two different approaches to be unconditionally stable in the choice of time stepsize Δt. The sub-stepping scheme [64] is expressing the velocity-correction scheme in a semi-Lagrangian form. It expresses the material derivative in the Lagrangian frame of reference while the pressure and viscous terms are in the Eulerian frame as for the above scheme. Therefore, the scheme mainly changes the initial evaluation of the advection term, but does not change the subsequent pressure and velocity problems.

Starting with the material derivative in the Lagrangian form

(14.5)

Reference [46] proposes to discretise the unsteady term in time with a Jth order BDF as

(14.6)

where the u_d^n-q corresponds to the velocity at the departure point x_d = (x_d,y_d,z_d) and at time instance tⁿ.

To find the velocity field at the departure point, the unsteady advection equation is solved in an auxiliary pseudo-time t^n+1-q < τ < tⁿ⁺¹ where q stands for the total number of iterations executed in pseudo-time τ.

(14.7)

where is a complementary velocity field. The field represents the influence of convection on the examined flow, u is the divergence-free advection velocity. Having computed this updated velocity field, we evaluate the advection term for equation (14.2) and (14.3) to solve for the new pressure and velocity as in the semi-implicit form.

The sub-stepping does only change the evaluation of the advection term N and is implemented inside the class SubsteppingExtrapolate.

The implicit velocity-correction schemes VCSImplicit transforms the velocity equation to be unconditionally stable with respect to the CFL condition. The key difference to the semi-implicit scheme is a linearisation of the advection operator. The linearisation assumes that uⁿ⁺¹ ⋅∇uⁿ⁺¹ ≈ũ ⋅∇uⁿ⁺¹ where the advection velocity is approximated by either an extrapolation of the form ũ ≈∑ _q=0^J-1u^n-q [65] or an updated velocity that uses the new pressure pⁿ⁺¹ [28].

The velocity problem with linearisation becomes an advection-diffusion-reaction (ADR) problem. It is defined as

un+1 - ν∇2un+1 + ⋅∇un+1 = ∑ q=0J-1α qun-q
-∇pn+1 + fn+1.		(14.8)

This scheme is a child of the semi-implicit scheme and implemented in the class VCSImplicit with an appropriate extrapolation in ImplicitExtrapolate.

14.1.3 Spatial discretisation/Implementation

We implement the various time discretisations above with a spectral hp element method. Therefore, we transform the strong form equations for pressure and velocity to the weak form using the L2 inner product with basis functions ϕ. For the pressure equation (14.2), we do the inner product with the derivative of the basis ∇ϕ which leads to the equation

∫ Ω∇ϕ ⋅∇pn+1	= ∫ Ωϕ∇⋅
	-∫ Γ⋅n.	(14.9)

Then, taking the L2 inner product with the basis ψ for the velocity equation in equation 14.4 leads to the weak form

∫ Ω∇ψ ⋅∇un+1 + ψun+1 = -∫ Ωψ.

(14.10)

The final discrete data of the variables u,v,w,p are saved within an array of ExpListSharedPtr named m_fields. Each ExpListSharePtr in m_fields holds a reference to a ContField object that saves the data in coefficient Coeff and physical Phys space. For example, one can access the w-velocity data via m_fields[2]->GetPhys() which return an Array of the data in physical space (i.e. at the quadrature points) for each element. Note that the pressure is always the hint-most entry in m_fields. Also, it can separately be accessed via m_pressure which is simply reference to the last entry in m_fields.