11.1 Synopsis

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:

(11.1a)

(11.1b)

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

11.1.1 Velocity Correction Scheme

The first approach uses a splitting/projection method where the velocity system and the pressure are typically decoupled. Splitting schemes are typically favoured for their numerical efficiency since the velocity and pressure are handled independently, requiring the solution of three (in two dimensions) elliptic systems of rank N (opposed to a single system of rank 3N solved in the Stokes problem). However, a drawback of this approach is the splitting scheme error which is introduced when decoupling the pressure and the velocity system, although this can be made consistent with the overall temporal accuracy of the scheme by appropriate discretisation of the pressure boundary conditions.

11.1.1.1 High order splitting scheme

In the original approach a stiffly-stable time integration was proposed in the work of Karniadakis, Israeli and Orszag [20]. This was then later fully analysed in the work of Guermond and Shen [17].

Briefly, high order splitting scheme was originally proposed in three steps involving explicit advection of the non-linear terms, followed by the solution of the pressure Poisson system and finally solving a Helmholtz problem to enforce the viscous terms and velocity boundary conditions. In the following however we briefly formulate this scheme as a two steps using a formulation outline by Guermond and Shen.

1. In the first step we formulate a weak pressure Poisson problem by taking the inner product over the solution domain Ω of equation (11.1a) with respect to the gradient of the test basis, ∇q, i.e.

   (11.2)

   where N(u) = u ⋅∇u. We recall that the term ∫ _Ω∇q ⋅∇p is the weak approximation to the Laplacian operator for pressure. To decouple this term from the velocity system a few steps are necessary. Using the identity

   

   we can enforce the divergence to be zero by setting the last term to zero. If we now integrate the 1st, 2nd and last term in equation (11.2) by parts we can obtain the weak pressure equation

   ∫ Ω∇q ⋅∇pn+1	= ∫ Ωq ∇⋅
   	-∫ ∂Ωq⋅n	(11.3)

   where ∂Ω is the boundary of the domain and we have used the factor that ∇⋅ (∇×∇×u) = 0. To get the final form of the weak pressure Poisson equation we can use a backward approximation of the time derivative to obtain

   (11.4)

   where ⁿ⁺¹ is an intermediate velocity upon which to decouple the system we impose that ∇⋅ⁿ⁺¹ = 0 and

   

   Finally we introduce a consistent extrapolation for the non-linear terms and the curl of vorticity terms of the form:

   

   A similar extrapolation can be used on the curl-curl term to end up with the final weak pressure approximation

   ∫ Ω∇q ⋅∇pn+1	= ∫ Ωq ∇⋅
   	-∫ ∂Ωq⋅n	(11.5)

   We note this can be recast into an equivalent strong form of the pressure Poisson equation of the form

   (11.6)

   with consistent Neumann boundary conditions prescribed as

   (11.7)

2. The second step is discretise equation (11.1a) at time level n + 1, use the pressure at n + 1 from the first step and solve for the velocity uⁿ⁺¹.

   In this step now approximate the time derivative using

   (11.8)

   which leads us to the Helmholtz problem

   (11.9)

This scheme is activated in the SolverInfo section with the SolverType specification:

11.1.1.2 Velocity Correction Scheme with a Weak Pressure formulation

As presented in the previous section in the work of Guermond and Shen [17] and subsequent work they formulate the pressure in a weak rather than strong form to obtain the pressure Poisson system. Therefore if we take the inner product of equation (11.1a) with respect to the gradient of the test space, ∇q, we obtain equation (11.2)

We again make the approximation

(11.10)

However this time we only integrate by parts the last term and do not integrate the non-linear term by parts. However we still need to enforce the condition that ∇⋅ⁿ⁺¹ = 0 and so we also integrate just this part of the time derivate by parts to arrive at a weak pressure system of the form:

∫ Ω∇q ⋅∇pn+1+	∫ ∂Ω0qn+1 ⋅n = ∫ Ω∇q ⋅ ( -N⋆,n+1)
	-∫ ∂Ω d⋃ ∂Ω0q ν(∇×∇×u)*,n+1 ⋅n + ∫ ∂Ωdqwn+1 ⋅n	(11.11)

where ∂Ω_d is the Dirichlet boundary conditions for the velocity and ∂Ω₀ is the outflow boundary.

This scheme is activated in the SolverInfo section with the SolverType specification:

11.1.1.3 Specifying pressure boundary conditions

In order to specify the pressure boundary conditions given by equation (11.7) or for the equivalent conditions in the VCSWeakPressure scheme the USERDEFINEDTYPE condition “H” can be used. Therefore a zero velocity wall boundary condition on boundary region 0 in two-dimensions can be specified as

11.1.1.4 Outflow boundary conditions

The most straightforward outflow condition is to specify fully developed conditions of ∇uⁿ⁺¹ ⋅n = 0 and p = 0 which can be specified as

However when energetic vortices pass through an outflow region one can experience instabilities as identified by the work of Dong, Karnidakis and Chryssostomidis [9]. In this paper they suggest to impose a pressure Dirichlet outflow condition of the form

(11.12)

with a step function defined by S_o(n ⋅u) = (1 - tanh ), where u₀ is the characteristic velocity scale and δ is a non-dimensional positive constant chosen to be sufficiently small. f_b is the forcing term in this case the analytical conditions can be given but if these are not known explicitly, it is set to zero, i.e. f_b = 0. (see the test KovaFlow_m8_short_HOBC.xml for a non-zero example). Note that in the paper [9] they define this term as the negative of what is shown here so that it could be use used to impose a default pressure values. This does however mean that the forcing term is imposed through the velocity components u,v by specifying the entry VALUE (An example can be found in ChanFlow_m3_VCSWeakPress_ConOBC.xml). For the velocity component one can specify

(11.13)

This condition can be enforced using the USERDEFINEDTYPE “HOutflow”, i.e.

Note that in the moving body work of Bao et al. [4] some care must be made to identify when the flow over the boundary is incoming or outgoing and so a modification of the term

is replaced with

where the default values are given by θ = 1,α₁ = 0,α₂ = 0 and these values can be set through the parameters OutflowBC_theta, OutflowBC_alpha1 and OutflowBC_alpha2.

Dong has also suggested convective like outflow conditions in [8] which can be enforced through a Robin type specification of the form

(11.14)

+ pn+1	= --ν(∇×∇×u)*,n+1 + N*,n+1⋅n
	-fn+1 + E(n,u*,n+1 + pn+1n - ν(∇⋅u*,n+1)n	(11.15)

11.1.1.5 Substepping/subcycling the Velocity Correction Scheme

It is possible to use different time steps in the velocity correction scheme using a substepping (also known as subcycling) [41] or auxiliary semi-Lagrangian approach [49]. Originally the scheme was proposed by Maday, Patera and Ronquist who referred to as an operator-integration-factor splitting method [28]

A schematic of the approach can be understood from figure 11.1.1.5 where we observe that smaller time steps can be used for the explicit advection steps whilst a larger overall time step is adopted for the more expensive implicit solve for the diffusion operator. More details of the implementation can be found in [49] and [41]. In the following sections we outline the parameters that can be used to set up this scheme. Since the explicit part is advanced using a DG scheme it is necessary to use a Mixed_CG_Discontinuous expansion with this option.

11.1.1.6 Approximation spaces for the velocity correction scheme

For well resolved simulations it appears that often using the same polynomial space for the pressure and velocity does give suitable answer but this does not satisfy the so-called LBB or inf-sup condition. Therefore, it is potentially better to specify an equivalent of the Taylor Hood approximation and use one higher polynomial order for velocity than the pressure with a continuous expansion. To specify this type of expansion you can use an expansion section of the form:

In the above example the “u,v” fields are specified to have a polynomial order of 7 using a modified expansion. Implicitly this form of the expansion definition uses a quadrature order of 9. The above definition then also uses a modified expansion for pressure but of polynomial order 6. Since currently for this solver to run we need to use a consistent quadrature order for both the velocity and pressure fields we specify the MODIFIEDQUADPLUS1 to tell the solver to use an additional quadrature point and therefore also use 9 quadrature points in each 1D direction for the pressure.

In other cases it is sometimes useful to run with an even higher quadrature order, for example to handle highly deformed elements where the Jacobian is represented by a polynomial expansion. This can be done by using a more detailed definition of the expansion of the form:

In this example we have specified an 8th order expansion for “u,v” and a 7th order expansion for “p”. The BasisType is given as “Modified_A, Modified_B” which is for a triangular expansion (note that for a quadrilateral expansion it would have been “Modified_A,Modified_A”) and so the number of quadrature points in this case is 9 in the first direction which uses Gauss-Lobatto-Legendre points but only 8 in the second direction since this uses a Gauss-Radau formula with α = 1,β = 0 weights (see [21] for details on why).

11.1.2 Immersed Boundary Methods: Smoothed Profile Method

The usual way to solve any PDE requires the definition of a well defined domain where the solution is to be determined. Thus, for complex geometries, the meshing process may get cumbersome and, in any case, the solver will have to the meshing process may get cumbersome, and likely struggle to avoid the presence of highly deformed and skewed elements. In addition, when solving cases with moving boundaries, the mesh has to be updated every time step to follow the shape of the boundaries, leading to a very resource and time-consuming simulation that limits the capabilities of the solver.

Immersed Boundary Methods may be very useful in these situations, where the definition of the boundaries requires a very complex mesh to reach convergence. The main idea behind them is the use of a forcing term in the incompressible Navier-Stokes equations in such a way that the mesh does not necessariliy follow the boundaries. The solution in the regions falling outside the boundaries is simply that of the boundaries, forcing the flow to behave as if there were a real object even if the mesh does not represent it. The method presented here is an adaptation of the Smoothed Profile Method [35] extended to a high-order semi-implicit splitting scheme [26, 47]. This method ensures that the no-slip, no-penetration and incompressibility constraints are mathematically enforced. Starting from the incompressible Navier-Stokes equations, the term f_s is added to the right hand side:

(11.16a)

(11.16b)

The definition of this term depends on the method but, for the Smoothed Profile Method (SPM), it is related to a shape function Φ(x,t) valued 0 in the fluid domain and 1 outside. It is usually defined as:

(11.17)

being ξ a scaling factor [47] and d(x,t) a function representing the distance to the boundary (positive inside the body, negative inside the fluid). If the case to be simulated includes more than one immersed boundaries, the final shape function is calculated by adding the individual ones as long as they do not overlap:

(11.18)

The approach followed during the implementation in Nektar++ is an extension of the Velocity Correction Scheme, using the final velocity obtained with this method as an intermediate velocity to determine the value of the forcing term. The initial equation (11.16a) is slightly modified and integrated in time by means of a stiffly-stable scheme and, then, split into different smaller parts that are solved separately:

(11.19)

(11.20a) (11.20b) (11.20c) (11.20d)

The new term f_s is defined as follows:

(11.21)

where α_q, β_q and γ₀ are coefficients of the stiffly-stable time integration method and u_p is the velocity of the points that lay outside the boundaries. Thus, the new term is just an acceleration proportional to the difference between the expected and the intermediate velocity, forcing the flow to follow the shapes defined by Φ and u_p. Some transformations in these expressions lead to the final SPM equations:

1. Advection and external forces:

   (11.22)

2. Pressure:

   (11.23a) (11.23b)

3. Viscous effects:

   (11.24)

4. SPM pressure:

   (11.25a) (11.25b)

5. SPM force:

   (11.26)

Since the term ∇p_p in the last equation may induce a velocity slightly different to u_p inside the bodies, it may be changed to (1 - Φ)∇p_p, cancelling this effect but adding some compressibitly to the flow [26]:

(11.27)

11.1.2.1 Input files

The session file follows the same rules as any other incompressible Navies-Stokes solver. However, there are some additional parameters that must be supplied when using this approach. First, the property SolverType must be set to SmoothedProfileMethod, while the immersed boundaries are defined in a function called ShapeFunction. Besides, the property ForceBoundary can be set to True if equation (11.27) is preferred over (11.26). Thus, the TIMEINTEGRATIONSCHEME and SOLVERINFO section could be similar to:

In addition, somewhere in the CONDITIONS section this function must appear:

As a brief guideline, to define a cylinder of radius 1 and center at the point (0,0) according to expression (11.17), the "Phi" field in the ShapeFunction function should be:

where the scaling coefficient has been set to 0.04. The variable names are compulsory, being Phi the shape of the bodies, and Up, Vp and Wp functions representing the velocity field inside them. The attribute USERDEFINEDTYPE is compulsory only if the functions depend on time, when it must be set to "TimeDependent".

For immersed boundaries with geometries that cannot be represented by means of analytical functions, an .stl binary file can be supplied as well. However, the geometry file must be first converted to .fld format with the phifile module of FieldConvert. The simplest way to proceed is by issuing the command:

where the value of scale corresponds to the coefficient ξ in equation (11.17). More details can be found in section 5.6.38. In any case, it is important to remark that this functionality is still under development.

In the ShapeFunction block of the session file, the line <E VAR="Phi" ... /> indicates that the immersed bodies are defined by the function introduced in VALUE, while a line like the following:

must be used when the Φ field is defined in an .stl file previously converted to .fld format. In this case, the solver only supports non-moving geometries and the attribute USERDEFINEDTYPE="TimeDependent", if specified, will not be used.

11.1.3 Direct solver (coupled approach)

The second approach consists of directly solving the matrix problem arising from the discretization of the Stokes problem. The direct solution of the Stokes system introduces the problem of appropriate spaces for the velocity and the pressure systems to satisfy the inf-sup condition and it requires the solution of the full velocity-pressure system. However, if a discontinuous pressure space is used then all but the constant mode of the pressure system can be decoupled from the velocity. When implementing this approach with a spectral/hp element discretization, the remaining velocity system may then also be statically condensed to decouple the so called interior elemental degrees of freedom, reducing the Stokes problem to a smaller system expressed on the elemental boundaries. The direct solution of the Stokes problem provides a very natural setting for the solution of the pressure system which is not easily dealt with in a splitting scheme. Further, the solution of the full coupled velocity system allows the introduction of a spatially varying viscosity, which arise for non-Newtonian flows, with only minor modifications.

We consider the weak form of the Stokes problem for the velocity field u = [u,v]^T and the pressure field p:

(11.28a)

(11.28b)

where the components of A,B and C are ∇ϕ_b,ν∇ub, ∇ϕ_b,ν∇ui and ∇ϕ_i,ν∇ui and the components D_b and D_i are q,∇ub and q,∇ui. The indices b and i refer to the degrees of freedom on the elemental boundary and interior respectively. In constructing the system we have lumped the contributions form each component of the velocity field into matrices A,B and C. However, we note that for a Newtonian fluid the contribution from each field is decoupled. Since the interior degrees of freedom of the velocity field do not overlap, the matrix C is block diagonal and to take advantage of this structure we can statically condense out the C matrix to obtain the system:

(11.29)

To extend the above Stokes solver to an unsteady Navier-Stokes solver we first introduce the unsteady term, ∂u∕∂t, into the Stokes problem. This has the principal effect of modifying the weak Laplacian operator ∇ϕ,ν∇u] into a weak Helmholtz operator ∇ϕ,ν∇u) -λ(ϕ,u where λ depends on the time integration scheme. The second modification requires the explicit discretisation of the non-linear terms in a similar manner to the splitting scheme and this term is then introduced as the forcing term f. For more details see [1, 42].

11.1.4 Linear Stability Analysis

Hydrodynamic stability is an important part of fluid-mechanics that has a relevant role in understanding how an unstable flow can evolve into a turbulent state of motion with chaotic three-dimensional vorticity fields and a broad spectrum of small temporal and spatial scales. The essential problems of hydrodynamic stability were recognised and formulated in 19th century, notably by Helmholtz, Kelvin, Rayleigh and Reynolds.

Conventional linear stability assumes a normal representation of the perturbation fields that can be represented as independent wave packets, meaning that the system is self-adjoint. The main aim of the global stability analysis is to evaluate the amplitude of the eigenmodes as time grows and tends to infinity. However, in most industrial applications, it is also interesting to study the behaviour at intermediate states that might affects significantly the functionality and performance of a device. The study of the transient evolution of the perturbations is seen to be strictly related to the non-normality of the linearised Navier-Stokes equations, therefore the normality assumption of the system is no longer assumed. The eigenmodes of a non-normal system do not evolve independently and their interaction is responsible for a non-negligible transient growth of the energy. Conventional stability analysis generally does not capture this behaviour, therefore other techniques should be used.

A popular approach to study the hydrodynamic stability of flows consists in performing a direct numerical simulation of the linearised Navier-Stokes equations using iterative methods for computing the solution of the associated eigenproblem. However, since linearly stable flows could show a transient increment of energy, it is necessary to extend this analysis considering the combined effect of the direct and adjoint evolution operators. This phenomenon has noteworthy importance in several engineering applications and it is known as transient growth.

In Nektar++ it is then possible to use the following tools to perform stability analysis:

• direct stability analysis;
• adjoint stability analysis;
• transient growth analysis;

11.1.4.1 Direct stability analysis

The equations that describe the evolution of an infinitesimal disturbance in the flow can be derived decomposing the solution into a basic state (U,p) and a perturbed state U + εu′ with ε « 1 that both satisfy the Navier-Stokes equations. Substituting into the Navier-Stokes equations and considering that the quadratic terms u′⋅∇u′ can be neglected, we obtain the linearised Navier-Stokes equations:

(11.30a)

(11.30b)

The linearised Navier-Stokes equations are identical in form to the non-linear equation, except for the non-linear advection term. Therefore the numerical techniques used for solving Navier-Stokes equations can still be applied as long as the non-linear term is substituted with the linearised one. It is possible to define the linear operator that evolved the perturbation forward in time:

(11.31)

Let us assume that the base flow U is steady, then the perturbations are autonomous and we can assume that:

(11.32)

Then we obtain the associated eigenproblem:

(11.33)

The dominant eigenvalue determines the behaviour of the flow. If the real part is positive then there exist exponentially growing solutions. Conversely, if all the eigenvalues have negative real part then the flow is linearly stable. If the real part of the eigenvalue is zero, it is a bifurcation point.

11.1.4.2 Adjoint Stability Analysis

The adjoint of a linear operator is one of the most important concepts in functional analysis and it plays an important role in understanding transition to turbulence. Let us write the linearised Navier-Stokes equation in a compact form:

(11.34)

The adjoint operator H^* is defined as:

(11.35)

Integrating by parts and employing the divergence theorem, it is possible to express the adjoint equations:

(11.36a)

(11.36b)

The adjoint fields are in fact related to the concept of receptivity. The value of the adjoint velocity at a point in the flow indicates the response that arises from an unsteady momentum source at that point. The adjoint pressure and the adjoint stream function play instead the same role for mass and vorticity sources respectively. Therefore, the adjoint modes can be seen as a powerful tool to understand where to act in order to ease/inhibit the transition.

11.1.4.3 Transient Growth Analysis

Transient growth is a phenomenon that occurs when a flow that is linearly stable, but whose perturbations exhibit a non-negligible transient response due to regions of localised convective instabilities. This situation is common in many engineering applications, for example in open flows where the geometry is complex, producing a steep variation of the base flow. Therefore, the main question to answer is if it exists a bounded solution that exhibit large growth before inevitably decaying. Let us introduce a norm to quantify the size of a perturbation. It is physically meaningful to use the total kinetic energy of a perturbation on the domain Ω. This is convenient because it is directly associated with the standard-L2 inner product:

(11.37)

where σ = . This is no other than the singular value decomposition of A(τ). The phenomenology of the transient growth can be explained considering the non-normality of the linearised Navier-Stokes evolution operator. This can be simply understood using the simple geometric example showed in figure 11.1.4.3. Let us assume a unit-length vector f represented in a non-orthogonal basis .This vector is defined as the difference of the nearly collinear vectors Φ₁ and Φ₂. With the time progression, the component of these two vectors decrease respectively by 20% and 50%. The vector f increases substantially in length before decaying to zero. Thus, the superposition of decaying non-orthogonal eigenmode can produce in short term a growth in the norm of the perturbations.

11.1.5 Steady-state solver using Selective Frequency Damping

To compute linear stability analysis, the choice of the base flow, around which the system will be linearised, is crucial. If one wants to use the steady-state solution of the Navier-Stokes equations as base flow, a steady-state solver is implemented in Nektar++. The method used is the encapsulated formulation of the Selective Frequency Damping method [19]. Unstable steady base flows can be obtained with this method. The SFD method is based on the filtering and control of unstable temporal frequencies within the flow. The time continuous formulation of the SFD method is

(11.38)

where q represents the problem unknown(s), the dot represents the time derivative, NS represents the Navier-Stokes equations, χ ∈ ℝ₊ is the control coefficient, q is a filtered version of q, and Δ ∈ ℝ₊^* is the filter width of a first-order low-pass time filter. The steady-state solution is reached when q = q.

The convergence of the method towards a steady-state solution depends on the choice of the parameters χ and Δ. They have to be carefully chosen: if they are too small, the instabilities within the flow can not be damped; but if they are too large, the method may converge extremely slowly. If the dominant eigenvalue of the flow studied is known (and given as input), the algorithm implemented can automatically select parameters that ensure a fast convergence of the SFD method. Most of the time, the dominant eigenvalue is not know, that is why an adaptive algorithm that adapts χ and Δ all along the solver execution is also implemented.

Note that this method can not be applied for flows with a pure exponential growth of the instabilities (e.g. jet flow within a pipe). In other words, if the frequency of the dominant eigenvalue is zero, then the SFD method is not a suitable tool to obtain a steady-state solution.