Efficient SVV stabilized triangular spectral element methods for incompressible flows of high Reynolds numbers

In this paper, we propose a spectral vanishing viscosity method for the triangular spectral element computation of high Reynolds number incompressible flows. This can be regarded as an extension of a similar stabilization technique for the standard spectral element method. The difficulty of this extension lies in the fact that a suitable definition of spectral vanishing viscosity operator in non-structured elements does not exist, and it is not clear that if a suitably defined spectral vanishing viscosity provides desirable dissipation for the artificially accumulated energy. The main contribution of the paper includes: 1) a well-defined spectral vanishing viscosity operator is proposed for non-standard spectral element methods for the Navier-Stokes equations based on triangular or tetrahedron partitions; 2) an evaluation technique is introduced to efficiently implement the stabilization term without extra computational cost; 3) the accuracy and efficiency of the proposed method is carefully examined through several numerical examples. Our numerical results show that the proposed method not only preserves the exponential convergence, but also produces improved accuracy when applied to the unsteady Navier-Stokes equations having smooth solutions. Especially, the stabilized triangular spectral element method efficiently stabilizes the simulation of high Reynolds incompressible flows.


Introduction
The spectral-element method is a high-order variational method which combines the geometric flexibility of finite-elements with the high accuracy of spectral methods. It possesses several good computational properties, such as fast evaluation due to its tensor product structure, diagonal mass matrices, and suitability for parallel computation. However, in order to use fast tensor product summation, the standard spectral-element method is usually restricted to quadrilateral (2D) or hexahedron (3D) partitions. This partition requirement is difficult to meet in certain situations like complex geometries and adaptive meshing. One way to overcome this drawback is to allow triangular (2D) or tetrahedron (3D) partitions, which are more flexible in handling complex geometries. Progress has been made in this direction in recent years, and some triangular spectral methods of different types have been developed based on (i) approximations by polynomials through mapping (see, e.g., [1][2][3][4][5]); (ii) approximations by polynomials using special nodal points such as Fekete points (see, e.g., [6][7][8][9]); and (iii) approximations by non-polynomial functions (see, e.g., [10][11][12]). In this paper we will focus on the triangular spectral method (TSM) based on approximations by rational polynomials. This type of TSM enjoys full tensor product structure, thus allowing fast summation in matrix-vector multiplications. It was first proposed and analyzed for elliptic problems in [12], then extended to the Stokes problem in [13]. Subsequently, an unstructured nodal spectral-element method for the Navier-Stokes equations was developed in [14]. This method is based on elementwise rational approximations in triangular or tetrahedral, and can be easily implemented using nodal basis. The aim of the paper is to investigate the applicability of the TSM to simulate incompressible flows of high Reynolds number. This is not a trivial task from our past experience in using standard rectangular spectral-element methods, and suitable dissipation mechanics need to be introduced to stabilize the calculation.
Indeed, high Reynolds number flows are difficult to compute, especially when using spectrally accurate numerical schemes as spectral methods [15]. This difficulty is linked to the fact that spectral methods are much less numerically diffusive than low-order methods. It has been found; see, e.g. [15,16], that the non-artificially dissipated energy accumulates at the high spatial frequencies may cause unstable calculation. One way to overcome this difficulty is to use stabilization techniques, as proposed in [15,16]. However in stabilizing the spectral element method, it is important to maintain the spectral accuracy of the algorithm. In standard rectangular spectral element methods, the spectral vanishing viscosity (SVV) approach has been proved to be an efficient stabilization method, which possesses the properties of the inter-element continuity and the spectral accuracy for smooth solutions. The SVV method was initially developed for the resolution of nonlinear conservation laws using standard Fourier spectral methods [17]. This method has been extensively investigated for periodic problems [18][19][20][21][22][23] and non-periodic problems in the framework of Legendre pseudo-spectral method [23][24][25][26][27][28][29]. Attempts to apply the SVV approach in different (but all standard) spectral element methods for simulation of incompressible flows have been made in [15,26,30].
The current work aims at making an attempt to apply the SVV approach in triangular spectral element methods (TSEM). The key is to define suitable SVV operators allowing fast evaluation on one side, and stabilization effect while maintaining the spectral accuracy on the other side. This is not at all obvious since, in standard spectral element methods, there exist several possibilities to define SVV operator in rectangular elements. For example, as pointed out in [15,30], some definitions of tensor product form may induce non-desired dissipative terms, and contrastly some definitions may have no dissipation effect in high frequencies at some direction. It is notable that a similar SVV approach has been used in [31] which considered the quadrilateral and triangular grids based on approximations by polynomials through mapping.
The purpose of this paper is first to well define the SVV operator in triangular domain. Then we show how to implement the SVV term in the framework of TSEM. The implementation has to take into account complex multidimensional geometries and vector valued functions, which make this point non-trivial. Finally, we will check the capabilities of the SVV-stabilized TSEM by simulating high Reynolds number incompressible flows in terms of accuracy and stability.
The paper is organized as follows. In the next section, we define the SVV operator in triangular domain and propose a new stabilized TSEM based on the SVV approach for the unsteady Navier-Stokes equations. The implementation details of the SVV-stabilized TSEM are given. Section 3 presents some numerical results which validate the accuracy and efficiency of the proposed method. Finally, the paper ends with some concluding remarks in Section 4. Extension to 3D case is presented in the Appendix.

Stabilized triangular spectral element method
This section is devoted to designing our SVV-stabilized triangular spectral element method for the Navier-Stokes equations. Hereafter we use boldface letters to denote vectors, vector functions, or vector spaces. Let be an open, connected and bounded domain in R 2 with boundary ∂ assumed to be Lipschitz continuous. We use L 2 ( ) to denote the space of square integrable functions in . The inner product of The norm and semi-norm of H 1 ( ) are denoted by u 1, and |u| 1, respectively. Let H 1 0 ( ) be the space of all functions in H 1 ( ) having vanishing trace on ∂ . Let us denote the velocity vector by u, the ratio between the pressure and the (constant) density by p, and let f be a forcing known term. The Navier-Stokes equation reads: subject to appropriate initial and boundary conditions. In the above equations, D t u denotes the material (Lagrangian) derivative of u with respect to time t, which can be expressed by ∂ t u + u · ∇u. ν is the dimensionless viscosity (the inverse of the Reynolds number).

Triangular spectral method
To clearly explain the idea, we start with a description of the spectral method in a single triangular domain : The weak formulation of the Navier-Stokes Eq. 1 in the triangular domain reads: find In order to well define the triangular spectral element approximation in space to the above weak problem, we will need some notations. The one-to-one transformation between and the square := (−1, 1) 2 is given by the Duffy mapping x = F(ξ ): with its inverse ξ = F −1 (x) from to by (ξ , η) is often referred to as collapsed coordinate system or the Duffy coordinates. It is an easy matter to compute the Jacobian determinant, denoted by J, of the mapping F: We associate a function u in with a function u in through The following formulas for the gradient operators will be useful: The approximation space to be used consists of the rational functions generated by polynomials in the reference square through the Duffy transform. Define the rational function R(x, y) in : whereR mn (ξ , η) be the polynomial in defined by: with J α,β k (ζ ), ζ ∈ being the Jacobi polynomial of degree k. Define the approximation spaces and their transformations as follows: Let ξ p , p = 0, 1, · · · , N, be the Legendre-Gauss-Lobatto points associated to L N , i.e., zeros of (1 − z 2 )L N (z); ω p , p = 0, 1, · · · , N, be the corresponding weights. We then define the discrete inner product (·, ·) N on : (9) where J is defined in (5). Let X N and M N be the approximation spaces: We now consider the rational spectral approximation to (2): Find u N ∈ X N and p N ∈ M N , such that It is notable that, similar to the standard spectral method for the Stokes problem, the pressure approximation space used in (11) is two degrees less than the velocity approximation space X N . This is to satisfy the well-known discrete inf-sup condition, which is necessary to avoid spurious pressure modes.

SVV stabilization
We first define the spectral vanishing operator in P N ( ), denoted by S, using the Legendre basis by [23], m N = N/2 [26], or N − 2 [32]. It is desirable to use a smooth variation forŜ n as: Then we define the SVV term − N ∂ x (S(∂ x u N )), which is written in weak form as follows: Note that the SVV term may be made symmetric: with the following definition of S 1/2 : The SVV operator in the 2D reference domain is defined in the following way. For 2 , which is defined in (24), the SVV term reads . Now we turn to define the SVV operator in the triangular domain . For u N , v N ∈ Q N ( ) 2 , we use the Duffy mapping (3) to associate the functions u N and v N through (6). Doing so allows to define the SVV operator by where G is the Jacobian of the mapping (4): We are now in a position to propose our SVV-stabilized TSM in single domain as follows: find u N ∈ X N , p N ∈ M N , such that where the stabilization term V N is defined in (12). In practice, it is highly beneficial to have the original diffusion term and the stabilization term combined together. Thus we propose to introduce the term Finally, the SVV-stabilized TSM for the Navier-Stokes equations reads: find u N ∈ X N , p N ∈ M N , such that

Implementation based on nodal basis
In this subsection, we give the details of the implementation of the SVV stabilization term T N (·, ·). The approach described here follows what is usually done when a nodal basis is chosen.
For notation convenience, we denote by u 1 ∈ X N and v 1 ∈ X N the first component of u N and v N respectively, respectively. The first component of T N (u N , v N ), denoted by T 1 N , can be written as where ξ pq = (ξ p , η q ), ω pq = ω p ω q . G 1 , G 2 , and G 3 are three geometric factors, defined as
Expressing u 1 on this basis, i.e., 0N h N (η), and choosing the test function v 1 ∈ X N to be each of the above basis functions, we arrive at the matrix statement of T 1 N , denoted still by T 1 N : Here D s = T 1/2 D with D being the Legendre differentiation matrix. The matrix form of the operator T 1/2 is defined by where M is the passage matrix from physical space to Legendre spectral space.

SVV stabilization in TSEM
We now briefly describe how to set up SVV-stabilized TSEM with triangle and rectangle mixed partition. Let be an open bounded polygonal domain, which is decomposed as: Let F k denote the mapping from the reference domain to k . In this case, the velocity and pressure approximation spaces are: where The SVV-stabilized TSEM in this spectral element case can be written in the same way as (13), with the SVV term taking now the element-wise sum as where with G k being the Jacobian of the mapping F −1 k , and J k the Jacobian determinant of the mapping F k .

Numerical results and discussions
In this section, we will investigate the convergence and stability property of the SVVstabilized TSEM. We first consider the Navier-Stokes equations with a stiff analytical solution to verify the accuracy. For this fabricated analytical solution, we observe that the SVV-stabilized TSEM not only keeps exponentially accurate, but also improves the accuracy as compared to the standard SEM. Then the flow through a backward facing step is simulated with the SVV-stabilized TSEM to show the stabilization capability of the method.
We use the classical semi-implicit splitting scheme for the temporal discretization, in which the non-linear convection terms are treated explicitly using third-order Adams-Bashforth; see, e.g., [33]. The spatial discretization makes use of the C 0 -TSEM with SVV stabilization, described in Section 2.4. The C 0 -continuity of the velocity across the spectral elements is naturally enforced by using a nodal basis. Similar to C 0 -finite element methods for the Navier-Stokes equations, the global matrix system can be built by assembling the local matrix system, together with the continuity conditions for the velocity, which is accomplished by requiring that the neighboring solutions share the same nodal values at the element interfaces [14].

Unsteady Navier-Stokes equations
First we investigate the accuracy of our stabilized triangular spectral element method with ν being fixed to be 0.01. In order to demonstrate the accuracy of the SVV-stabilization method for the Navier-Stokes Eq. (1), we fabricate the exact solution as follows: u(x, y, t) = sin(t) sin(x) cos(y), v(x, y, t) = − sin(t) cos(x) sin(y), p(x, y, t) = sin(t) sin(x) sin(y).
The triangular spectral element mesh with N = 8 in the computational domain := 2 is shown in Fig. 1. Figure 2 shows the velocity and pressure errors in H 1 , L ∞ and L 2 norms versus the polynomial degree N at t = 1 obtained with the SVV-stabilized TSEM and the standard SEM respectively. In this test we set the parameters m N = N − 2 and N = 1/N. We observe that the SVV stabilized-SEM not only preserves the exponential convergence rate, but also is more accurate than the standard SEM. This phenomenon is in contrast with what we have obtained in the previous linear elliptic equation. The reason of this desirable result is that finite mode approximation to the solution causes accumulation of the undissipated energy at the high-frequency modes of the numerical solution through the nonlinear convection terms. These spurious high-frequency modes resulting from aliasing effects of the convection terms may be damped when the SVV is activated in the SVV-stabilized TSEM, leading to more accurate results. It is worth mentioning that a similar result has been obtained for the SVV-stabilization in the standard SEM [15].
The error history during the time integration is displayed in Fig. 3, in which the errors   are compared between the SVV-stabilized TSEM and the standard SEM. It is observed that the accuracy of the SVV-stabilized TSEM is higher than the standard SEM during all time. The similar result is observed in Fig. 4 with ν = 0.001.

Flow through a backward facing step
In this example, we apply the SVV-TSEM to simulate the flow through a backward facing step by solving the Navier-Stokes Eq. 1 in the domain shown in Fig. 5. The Dirichlet boundary condition g and initial condition u 0 in (1) for the backward facing step flow are specified as follows: The Reynolds number is defined as Re =ū(H − h)/ν, whereū is the average velocity at the entrance, H and h are respectively the height of the outlet and the entrance. Thus H−h is the height of the step. The expansion ratio of the step is h : H = 2 : 3. It is seen from (17) that the profile of the inflow boundary condition is taken parabolic. The outlet boundary is taken far away from the step (22 step heights) to avoid possible artificial reflection. The mixed triangular and quadrangular spectral-element mesh is shown in Fig. 6. An enlarged view of some elements near the step is displayed in Fig. 7.
In our calculations, we fix the element number to K = 40 and the polynomial degree to N = 8. The simulation is performed for two Reynolds numbers: Re = 10000 and Re = 50000 (corresponding to ν = 2.666666 −05 and ν = 5.333332 −06 respectively). We would like to emphasize that the calculation by the standard SEM for the flow at Re = 50000 is unstable. Therefore we only present the numerical results obtained with the SVV-stabilized TSEM. In Figs. 8 and 9 we present the vector of the velocity, streamline contours, and pressure contours captured at t = 1.42 for Re = 10000 and Re = 50000 respectively. Enlarged views of the computed solutions near the step are given in Figs. 10, 11, and 12 respectively. The well resolved flow structures can be readily observed in these figures, and compare well with the existing results in the literature. This test confirms the accuracy and the stabilization effect of the proposed SVV-TSEM for simulating incompressible flows at high Reynolds numbers.

Conclusions
We have proposed a SVV stabilization for non-standard spectral element methods based on triangular or triangular/quadrangular mixed meshes. The main goal is to design a stable high order method for the spatial discretization of the Navier-Stokes equations. The contribution of this paper is threefold: Firstly, the SVV operator is suitably defined for a triangular domain; Secondly, we have shown how to efficiently implement the SVV operator in the framework of triangular spectral element methods. The implementation had to take into account complex multidimensional geometries and vector valued functions, which made this point non-trivial. Thirdly, the efficiency of the proposed stabilized TSEM was verified through several numerical examples, including linear elliptic equation, Navier-Stokes equations in simple domains, and incompressible flows at high Reynolds numbers. The accuracy, convergence order, and stability property of the method were carefully examined. Our numerical finding is: 1) the proposed SVV-stabilized TSEM keeps the exponential convergence as standard SEM for smooth solutions; 2) the SVVterm may cause a loss of accuracy for linear elliptic equations, thus is not suggested to use. In fact, SVV-stabilization is usually unnecessary for linear equations; 3) SVV-stabilized TSEM significantly improved both the accuracy and stability of the incompressible flow simulation at high Reynolds number. In particular, SVV-stabilized TSEM allows to stably simulate flows of high Reynolds numbers for which the standard SEM fails. We would like to emphasize that the method developed in the paper is directly applicable to practical flow simulations in aerodynamics, such as airfoil flows and turbulence, as one has done with the traditional spectral element method. One more point worth mentioning is that it was reported in [34] that the use of some SVV kernels may lead to undesirable property. We will investigate this phenomenon in future work.

Appendix: Extension to the three-dimensional case
We describe in this appendix some basic ingredients for developing non classical 3-D spectral-element methods for the Stokes equation.
The one-to-one transformation from to is given by the 3-D Duffy mapping x = F(ξ ): with its inverse ξ = F −1 (x) from to by We collect below some properties of the 3-D Duffy mapping: , Inversely, From the above, one easily finds the determinant of the Jacobian for (18): Therefore, we have and inversely, It follows from (21) that if u ∈ H 1 ( ), then necessarily A direct calculation from (21) gives (∇u, ∇v) We shall use rational functions generated by polynomials in the reference cube through the 3-D Duffy transform. More precisely, letR lmn (ξ , η, ζ ) be the polynomial in defined by:R we define the rational function R(x, y, z) in by the Duffy transformation of R lmn (ξ , η, ζ ), i.e., The approximation spaces and their transformations will be as follows: By the properties of the Jacobi polynomials, we can easily verify the following orthogonality relation: R lmn (x, y, z)R l m n (x, y, z)dxdydz Any function u ∈ L 2 ( ) can be expressed as with the coefficient u lmn given by On the other hand, we haveũ ∈ L 2 ( ) if and only if u ∈ L 2 ( ), where the weight function (ξ, η, ζ ) := (1 − η)(1 − ζ ) 2 64 is the Jacobian defined in (20). Similarly, we have the following expression forũ: with u lmn given in (25) or equivalently in the alternative form: u lmn = 1 γ lmn ũ(ξ , η, ζ )R lmn (ξ , η, ζ ) (ξ, η, ζ )dξ dηdζ .

A.2 Spectral-element method in 3D domains with tetrahedron partition
We now turn to describe the case of general 3D domains with partition. Let be a polyhedron which is partitioned into a number of tetrahedrons: We assume that the partition is conforming in the usual sense. In this case, the velocity and pressure approximation spaces are defined by: where F k is the mapping from to k .
Then the SVV-stabilized tetrahedron spectral-element method can be constructed and implemented in a similar way as 2D case.