An implicit dual-time stepping high-order nodal discontinuous Galerkin method for solving incompressible flows on triangle elements

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation 168 (2020) 173–214 www.elsevier.com/locate/matcom

Original articles

An implicit dual-time stepping high-order nodal discontinuous Galerkin method for solving incompressible flows on triangle elements M. Hajihassanpour 1 , K. Hejranfar ∗,2 Aerospace Engineering Department, Sharif University of Technology, Iran Received 27 December 2018; received in revised form 11 July 2019; accepted 27 August 2019 Available online 2 September 2019

Abstract In this work, a high-order nodal discontinuous Galerkin method (NDGM) is developed and assessed for the simulation of 2D incompressible flows on triangle elements. The governing equations are the 2D incompressible Navier–Stokes equations with the artificial compressibility method. The discretization of the spatial derivatives in the resulting system of equations is made by the NDGM and the time integration is performed by applying the implicit dual-time stepping method. Three numerical fluxes, namely, the local Lax–Friedrich, Roe and AUSM+-up are formulated and applied to assess and compare their accuracy and performance in the simulation of incompressible flows using the NDGM. Several steady and unsteady incompressible flow problems are simulated to examine the accuracy and robustness of the proposed solution methodology and they are the Kovasznay, backward facing step, NACA0012 airfoil, circular cylinder and two side-by-side circular cylinders. Indications are that the NDGM applied for solving the incompressible Navier–Stokes equations with the artificial compressibility approach and the implicit dual-time stepping method is accurate and robust for the simulation of steady and unsteady incompressible flow problems. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Nodal discontinuous Galerkin method; Incompressible flows; Artificial compressibility method; Implicit dual-time stepping method; Numerical flux; Triangle elements

1. Introduction Many practical applications in aerodynamic such as the airflow over cars, trains, motorcycles and etc. can be considered as incompressible flow problems. Moreover, in hydrodynamic, most of the flows are incompressible because the sound-speed in liquids is much greater than the liquid velocity. For instance, water flow over submarines, in pipes, in pumps and etc. often are incompressible. Therefore, the simulation of the incompressible flows is of prime importance. There is a difficulty for numerically simulating incompressible flows due to the decoupling of the continuity and momentum equations, known as the pressure–velocity coupling issue, and how to satisfy the continuity equation to ∗ Corresponding author.

E-mail address: [email protected] (K. Hejranfar). Ph.D. Candidate 2 Professor 1

https://doi.org/10.1016/j.matcom.2019.08.011 c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

174

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

be divergence free. Hence, each numerical method developed for computing incompressible flows should address this issue appropriately. Different solution methods have been applied in the literature to remedy this issue. Among these methods, the artificial compressibility approach introduced by Chorin [8] is an appropriate method which overcomes the pressure–velocity coupling problem for the solution of incompressible flows. In this approach, the incompressible governing equations are modified in a manner to be more like to the compressible formulation. For this aim, the time derivative of the density term is added to the continuity equation and by using an equation of the state and defining the artificial sound speed, the pressure and density variables are related and the unknown are the pressure and the velocity vector to be solved. Consequently, the eigenvalues of the resulting system of equations become real and therefore the type of the system of equations is changed to be hyperbolic. The previous studies have indicated that this approach is suitable for the simulation of the incompressible flow problems with relatively complex features and different applications (for example see Refs. [29,30,35]). Besides the difficulty in the modeling of incompressible flow, there are some other numerical difficulties that should be properly addressed. Most of the industrial applications dealing with incompressible flows have complex geometries, and therefore, this is valuable to develop a numerical method capable of handling complex geometries. Thus, working with unstructured grids is a very important property of a numerical method. Moreover, a numerical method which provides different orders of accuracy from low order to high-order accuracy, is worthwhile. Furthermore, in the past two decays, the computer architectures have been changed significantly and they are multicore today, and thus, the development of the numerical methods with good scalability is more important now. The discontinuous Galerkin (DG) method can be used on unstructured grids to handle complex geometries with different orders of accuracy, also it has a great scalability for parallel computing purposes. In the DG context, the variables can be estimated by either nodal [21,22] or modal [5,11] forms. The nodal form of the DG method has been used by Giraldo et al. [21] for the solution of the shallow water equations on the quadrilateral grids. The extension of the method to the triangular grids has been performed by Hesthaven and Warburton [31]and they have solved the Maxwell equations. Moreover, this form of the DG method has been used in combination with the different time integration schemes [33,57], parallel computing models [18,34,44], grid types [62] and applications [4,9,16,27,43,48,58,63]. Different numerical methods have been introduced in the literature for computing incompressible flows using the artificial compressibility method such as the finite-difference (FD) methods [51], the compact methods [28,50], the finite volume (FV) methods [36,42], the immersed boundary methods (IBM) [19], the smoothed particle hydrodynamics (SPH) methods [54] and the finite element (FE) methods [47]. In the context of the DG method, different methods such as the artificial compressibility flux [1,13,64,65], vorticity-stream function [40], projection methods [2,17], staggered semi-implicit method [59], density-based method [53] and hybrid methods [20,37,45] have been used for computing the incompressible flows. Note that in the works of Bassi et al. [1,13] they have used the artificial compressibility term only in the cell interface for the flux calculations. Then, they have used an exact solver for the Riemann problem (an iterative procedure) to obtain the fluxes between the elements. In the works of Zhang et al. [64,65] they have discretized the inviscid term based on the simplified artificial compressibility method and they have used the direct discontinuous Galerkin method for the calculation of the viscous term. In the present work, the artificial compressibility method is utilized to provide a truly incompressible flow solver and the resulting system of the equations is solved by the NDGM at all the solution points, not only at the flux points. In addition, three numerical fluxes, namely, the Lax–Friedrich, Roe and AUSM+-up are formulated and applied in the NDGM and the accuracy and performance of these numerical fluxes are also examined in comparison with each other. An important issue when applying the DG method is the discretization procedure of the time derivative terms. In the absence of the time derivative terms, the Newton–Raphson iterative method with quadratic convergence can be used for solving the nonlinear system of equations [49]. Note that in the NDGM with the artificial compressibility method used in the present study, the time derivative terms are retained to overcome the pressure–velocity coupling problem, and therefore, the Newton–Raphson iterative method cannot be applied here. Different time integration methods such as the implicit Euler [3], explicit Runge–Kutta [10] and exponential-time differencing [38] methods have been applied in the literature in the context of the DG method. It should be noted that the implicit Euler and exponential-time differencing methods provide a larger allowable CFL number in comparison with the explicit Runge–Kutta method. Although these two implicit time integration methods require lower number of iterations to attain the steady-state solution, they need a matrix inversion which makes the flow solver to be complicated. In [38], it was shown that the CPU time of the explicit third-order TVD Runge–Kutta method is more than the

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

175

implicit backward Euler and exponential-time differencing methods when a low-order polynomial is used, however, they are closed to each other for a high-order polynomial. Note also that the implementation of the explicit methods is much simpler than the implicit methods, and therefore, the explicit schemes are adapted to be applied for the time integration of the NDGM with the artificial compressibility method. It should be noted that the NDGM with the artificial compressibility method used in the present study is just appropriate for steady incompressible flows because the time derivative terms appeared in the formulation are the artificial ones. Here, the implicit dual-time stepping method which does not need a matrix inversion is incorporated in the formulation of the NDGM for the present solution method to be capable of accurately computing the unsteady incompressible flows. Both the steady and unsteady incompressible flows are simulated with the proposed solution methodology. The paper is prepared as follows: In Section 2, the incompressible Navier–Stokes equations based on the artificial compressibility method are given. The numerical method applied based on the nodal discontinuous Galerkin method, the numerical fluxes formulated, time integration and the boundary conditions used are presented in Section 3. The numerical results of several steady and unsteady incompressible flow problems together with the investigation of the effects of the polynomial degree and the numerical fluxes on the accuracy and performance of the solution method are given in Section 4. Finally, some conclusions are drawn in Section 5. 2. Governing equations For the two-dimensional incompressible flows, the continuity equation is written as follows: ∂u ∂v + =0 ∂x ∂y

(1)

and to obtain the artificial compressibility formulation introduced by Chorin [8], an artificial time derivative of the density is added to the continuity equation as ∂ρ ∂u ∂v + + =0 ∂τ ∂x ∂y

(2)

where τ is the artificial time, ρ is the density and u and v are the velocity components in the x- and y- direction, respectively. Then, using the chain rule we have ∂ρ ∂ρ ∂ p = ∂τ ∂ p ∂τ where p is the pressure. Now, by defining the artificial sound speed as β = artificial compressibility method is obtained as ∂v 1 ∂ p ∂u + + =0 2 β ∂τ ∂x ∂y

(3) √

∂p , ∂ρ

the continuity equation using the

(4)

The appearance of the pressure term in the continuity equation couples the velocity field to the pressure field and overcomes the pressure–velocity coupling problem. This is obvious from Eq. (4) that in the steady-state condition the divergence free velocity field is satisfied which is important in incompressible flow simulations. The parameter β determines the convergence rate of the solution and its value is selected between 1 to 10, or may be chosen variable [42]. In the present study, the parameter β is set to be constant and a value of 1.4 is considered for all the simulations performed. The resulting system of governing equations for the simulation of the incompressible flows by the artificial compressibility method can be written in the Cartesian coordinates and dimensionless form as ( ) ( ) ∂ β 2v ∂ p ∂ β 2u + + =0 (5) ∂τ ∂x ∂y ( 2 ) ( 2 ) ∂ u +p ∂u ∂ (uv) 1 ∂ u ∂ 2u + + = + (6) ∂τ ∂x ∂y Re ∂ x 2 ∂ y2 ( 2 ) ( ) ∂v ∂ (uv) ∂ v + p 1 ∂ 2v ∂ 2v + + = + (7) ∂τ ∂x ∂y Re ∂ x 2 ∂ y2

176

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

where Re = L ref Vref /ν denotes the Reynolds number where L ref , Vref and ν are the characteristic length, the characteristic velocity and the kinematic viscosity, respectively. Note that the above formulation with the artificial compressibility method is just suitable for steady incompressible flows because of appearing the artificial time τ in the formulation. For unsteady computations, one can implement the implicit dual-time stepping method in the problem formulation with the artificial compressibility method. In the present study, both the steady and unsteady flows are considered. 3. Numerical method 3.1. Nodal discontinuous Galerkin method In this section, the nodal discontinuous Galerkin method (NDGM) is applied for solving the steady incompressible flows with the artificial compressibility method on triangle elements and its extension for the unsteady incompressible flows will be given in Section 3.4. In the NDGM, the computational domain Ω , which has the boundary ∂Ω , is divided into K non-overlapping cells. The cells are the straight-sided triangles with arbitrary sizes and they are denoted by D i and therefore, Ω can be approximated by the union of D i Ω=

K ⋃

Di

(8)

i=1

The local approximate solution for each element can be expressed as the nodal or modal description. The nodal description as a polynomial of degree N is expressed as follows: u˜ i (x, y, t) =

Np ∑

( ) u˜ i x ij , y ij , t l ij (x, y)

(9)

j=1

where N p = (N + 1) (N + 2) /2 is the number of nodes in an element. However, the modal description is u˜ i (x, y, t) =

Np ∑

ˆ u ij (t) ψ ij (x, y)

(10)

j=1

where ψ j , ˆ u j and l j represent the basis function, the expansion coefficient and the shape function, respectively, and the superscript ˜ denotes the approximated solution. This is more convenient to express the local approximate solution in the computational space than the physical space and then to map the results to the physical space. The local approximate solutions in the nodal and modal forms are u˜ i (r, s, t) =

Np ∑

) ( u˜ i r ij , s ij , t l j (r, s)

(11)

ˆ u ij (t) ψ j (r, s),

(12)

j=1

u˜ i (r, s, t) =

Np ∑ j=1

respectively. On the above two equations, r = [−1, 1] and s = [−1, 1] are the spatial coordinates in the computational space. The computational space is a right angle triangle. The nodes in the computational space are defined somehow to minimize the Lebesgue constant defined as Λ = max

Np ∑ ⏐ ⏐ ⏐l j (x, y)⏐

(13)

j=1

These nodes along the edge of the triangle correspond to the Legendre–Gauss–Lobatto nodes in one-dimension [27,32]. The basis function is obtained in some way that the canonical basis function ψm (r, s) = r i s j ,

(i, j) ≥ 0; i + j ≤ N i m = j + (N + 1) i + 1 − (i − 1) , (i, j) ≥ 0; 2

(14) i+j≤N

(15)

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

177

to be orthonormalized by using a Gram–Schmidt process, which yields √ ψm (r, s) = 2Pi(0,0) (a) P j(2i+1,0) (b) (a − b)i (16) 1+r − 1, b = s (17) a=2 1−s and Pn(Π ,Υ ) is the nth order Jacobi polynomial and for Π = Υ = 0 it is Legendre polynomial. A mapping is used to connect the nodes in the physical and computational spaces to each other, which is r +s 1 r +1 2 s+1 3 x (r, s) = − v + v + v (18) 2 x 2 x 2 x r +s 1 r +1 2 s+1 3 v + v + v (19) y (r, s) = − 2 y 2 y 2 y where (vx1 , v 1y ), (vx2 , v 2y ) and (vx3 , v 3y ) are the Cartesian coordinates of the vertices of the triangle in the physical space. The Jacobian is defined as the ratio of the triangle area in the physical space to one in the computational space, and it is ∂x ∂y ∂x ∂y − (20) J= ∂r ∂s ∂s ∂r The Vandermonde matrix υ connects the basis function to the shape function through this relation υ T l (r, s) = ψ (r, s)

(21)

which the Vandermonde matrix is defined as υi, j = ψ j (ri , si ). Now, the shape function can be obtained from the above equation. An appropriate application of the NDGM to the governing equations requires changing the second order partial differential terms to the first order ones using the change of the variables, so that ( ) ( ) ∂ β 2v ∂ p ∂ β 2u + + =0 (22) ∂τ ∂x ∂y ( 2 ) ( ) ∂ u +p ∂u ∂ (uv) 1 ∂q1 ∂q2 + + = + (23) ∂τ ∂x ∂y Re ∂ x ∂y ( 2 ) ( ) ∂v ∂ (uv) ∂ v + p 1 ∂q3 ∂q4 + + = + (24) ∂τ ∂x ∂y Re ∂ x ∂y ∂u (25) q1 = ∂x ∂u q2 = (26) ∂y ∂v q3 = (27) ∂x ∂v q4 = (28) ∂y Here, the details of the implementation of the NDGM to the governing equations are described for Eqs. (23) and (25) and the extensions to the remaining equations are straightforward. By defining E = u 2 + p and G = uv and applying the weighted residual method to these two equations give ( )] ∫∫ [ i ∂ u˜ ∂ E˜ i ∂ G˜ i 1 ∂ q˜1 i ∂ q˜2 i + + − + φ i d xdy = 0 (29) ∂τ ∂x ∂y Re ∂ x ∂y ] ∫∫ [ ∂ u˜ i i q˜1 − φ i d xdy = 0 (30) ∂x where φ i (x, y) denotes the test function. The nodal version of the mass and stiffness matrices can be obtained if the test and shape functions are set equal to each other, and therefore ( )] ∫∫ [ i ∂ u˜ ∂ E˜ i ∂ G˜ i 1 ∂ q˜1 i ∂ q˜2 i + + − + lmi d xdy = 0 (31) ∂τ ∂x ∂y Re ∂ x ∂y

178

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

] ∫∫ [ ∂ u˜ i i q˜1 i − l d xdy = 0 ∂x m Using the chain rule gives ( ) ( ( ) ( ∫∫ [ i i i i i ) i i ˜ ˜ ∂ u ˜ 1 ∂ q˜1 i lmi 1 ∂l i ∂ G l ∂ q˜2 i lmi ∂ E l m m i i ∂lm i ∂l m ˜ ˜ lm − + q˜1 i m + + − E +G + ∂τ ∂x ∂y ∂x ∂y Re ∂x ∂y Re ∂x )] ∂l i d xdy = 0 +˜ q2 i m ∂y ] ∫∫ [ ∂ u˜ i lmi ∂l i lmi q˜1 i − + u˜ i m d xdy = 0 ∂x ∂x and the surface term is reduced to the flux term by using the divergence or Gauss’ theorem as follows: ( ( ∫∫ [ ∮ [( i i ) i )] i ˜i 1 ∗ i ∂l m i ∂u i ∂lm i ∂l m i ∂lm ˜ ˜ lm − E +G + q˜2 nˆ ix E˜ i lmi + q˜1 d xdy = − ∂τ ∂ x ∂ y Re ∂ x ∂ y D ∂D ) )] 1 ( i i∗ i i ˜ i∗ i i i∗ i +nˆ y G lm − nˆ q˜1 lm + nˆ y q˜2 lm ds Re x ] ∮ [ ∫∫ [ ] ∂l i ∗ nˆ ix u˜ i lmi ds lmi q˜1 i + u˜ i m d xdy = ∂x ∂D D

(32)

(33)

(34)

(35)

(36)

where ˆ n = (nˆ x , nˆ y ) is the unit outward pointing normal and ds is an infinitesimal segment along the edge of the triangle. The substitution of Eq. (9) into Eqs. (35) and (36) yields ( ) ⎡ Np ∫∫ ∂ u˜ i x ij , y ij , t l ij (x, y) ( ( ∑ ) ∂l i (x, y) ⎣lmi (x, y) − E˜ i x ij , y ij , t l ij (x, y) m ∂τ ∂x D j=1 ) ( ( ) ( ) ∂l i (x, y) 1 ∂l i (x, y) +G˜ i x ij , y ij , t l ij (x, y) m + q˜1 i x ij , y ij , t l ij (x, y) m ∂y ∂x ⎤ Re 3×N p f ∮ ) (37) [( i ∑ ( ) ) ∂l (x, y) ⎦ ∗ ( +˜ q2 i x ij , y ij , t l ij (x, y) m d xdy = − nˆ ix E˜ i x ij , y ij , t l ij (x, y) lmi (x, y) ∂y ∂D j=1 ) ( ( ) ( ) 1 ∗ ∗ +nˆ iy G˜ i x ij , y ij , t l ij (x, y) lmi (x, y) − nˆ ix q˜1 i x ij , y ij , t l ij (x, y) lmi (x, y) )] Re ) ∗ ( +nˆ iy q˜2 i x ij , y ij , t l ij (x, y) lmi (x, y) ds, 1 ≤ m ≤ N p ] Np ∫∫ [ ∑ ) ( ) ( ∂l i (x, y) d xdy lmi (x, y) q˜1 i x ij , y ij , t l ij (x, y) + u˜ i x ij , y ij , t l ij (x, y) m ∂x D j=1 (38) 3×N p f ∮ ] [ ∑ ( i i ) i i i i∗ = nˆ x u˜ x j , y j , t l j (x, y) lm (x, y) ds, 1 ≤ m ≤ N p j=1

∂D

where N p f = N + 1 is the number of points on an edge of the triangle. In these equations, the numerical fluxes are denoted by the superscript ∗ and this form of the equations is known as the weak form. Now, the following integrals have to be solved ∫∫ ∫∫ [i ] [ ] i i lm (x, y) l j (x, y) d xdy = J lm (r, s) l j (r, s) dr ds = J i M (39) ] [ ( )] ∫ ∫ D[ ∫ ∫ ∂l i (x, y) ∂lm (r, s) ∂r ∂lm (r, s) ∂s l ij (x, y) m d xdy = J i l j (r, s) + dr ds ∂ x ∂r ∂ x ∂s ∂x D ( ) (40) ∂r ∂s = J i Sr + Ss ∂x ∂x

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

( ] )] ∫∫ [ ∫∫ [ ∂l i (x, y) ∂lm (r, s) ∂r ∂lm (r, s) ∂s l j (r, s) l ij (x, y) m d xdy = J i dr ds + ∂ y) ∂r ∂y ∂s ∂y D ( ∂r ∂s = J i Sr + Ss ∂y ∂y

179

(41)

The mass M and two stiffness matrices Sr and Ss are computed in the computational space, so the superscript i which stands for the ith triangle in the physical space is omitted in these matrices. Note that the dimension of the mass matrix in the left-hand side of the discretized equations is N p × N p while in the right-hand side (from the boundary terms) is N p × 3N p f . The mass and stiffness matrices can be evaluated either analytically or numerically by utilizing the Gaussian quadrature rule and here, they are evaluated analytically. By substitution Eqs. (39) to (41) into Eq. (37) the semidiscrete form of Eq. (23) will be obtained. Moreover, the substitution of Eqs. (39) to (41) into Eq. (38) results the discrete form of Eq. (25). Therefore, applying such a procedure for Eqs. (22) to (28) results in the semidiscrete form of the governing equations used for computing the incompressible flows as follows: i

∂˜ Q i =˜ R ∂τ i i where ˜ Q = ( p˜ i , u˜ i , v˜ i ) is the solution vector and ˜ R is the right-hand side vector.

(42)

3.2. Numerical flux To have a stable and accurate method, the numerical flux terms in the discretized form of the governing equations ∗ have to be calculated correctly. Therefore, for instance, the numerical fluxes in Eqs. (37) and (38) namely E˜ i and ∗ ∗ ∗ ∗ G˜ i as the inviscid fluxes and q˜1 i , q˜2 i and u˜ i as the viscous fluxes should be computed suitably. For the inviscid fluxes, different methods have been introduced in the literature such as the local LF (Lax–Friedrich), Roe, AUSM, HLL, HLLC and etc. [27,52]. Here, the local LF method, Roe and AUSM+-up are formulated to be used for the inviscid fluxes. The local LF or Rusanov numerical flux is a simple one ∗and it is widely used in the literature. If one considers ∗ i∗ i ˜i = ˆ the inviscid flux as H ni .˜ F , where the flux vector ˜ F , then the local LF numerical flux is defined as r iz 1 ( ( i ) ) ∗ i ˜ li + H ˜ ir + Λ ˜ ˜ i = {{ H ˜ i }} + Λ ˜ Q = H Ql − ˜ Qr (43) H 2 2 2 r iz ˜ i }} and ˜ Q are the average and jump, respectively. The variable Λ determines the maximum value of where {{ H the eigenvalues as Λ = max(Vli + cli , Vri + cri ) (44) √ √ where V = u 2 + v 2 is the velocity magnitude, ci = V i 2 + β 2 is the effective sound speed, and the subscripts l and r stand for the left and right boundaries∗ (edge of the triangle), respectively. ˜ i , the Roe average values for u, To calculate the Roe numerical flux H ¯ v¯ and p¯ should be calculated. These values based on the values in the interface are ˜ ui + ˜ u ri u¯ i = l (45) 2 v˜ i + v˜ri v¯ i = l (46) 2 pi + p˜ri p¯ i = l (47) 2 Moreover, the averaged eigenvalues λ and the eigenvectors σ are needed. The eigenvalues of the governing equations is ⎡ ⎤ i √Vn ⎢ i ⎥ i V + V i 2 + β 2⎥ λ =⎢ (48) ⎣ n √ n ⎦ 2 i Vn − Vni + β 2

180

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 1. Sequence of the grids used to solve the Kovasznay problem with the characteristic lengths, (a) h = 0.5, (b) h = 0.25 and (c) h = 0.125.

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

181

Fig. 2. Spectral convergence of the solution for the Kovasznay problem based on the (a) L ∞ - and (b) L 2 -norm errors.

where Vni = u¯ i nˆ ix + v¯ i nˆ ix is the velocity normal to the interface. The matrix of the eigenvectors (columns are eigenvectors) is ⎡ ⎤ 1 1 √ 0 − √ ⎢ ⎥ ⎢ ⎥ 2 Vni 2 + β 2 2 Vni 2 + β 2 ⎢ ⎥ ( ) ( ) ⎢ ⎥ √ √ ⎢ ⎥ 2 2 2 ⌢i 2 ⌢i ⎥ i i i 2 i 2 ⎢ v¯ Vn + Vn + β + β n y u¯ Vn − Vn + β + β n x ⎥ ⌢i ⎢ n i y σ = ⎢− (49) ⎥ ( ) ( ) ⎢ ⎥ 2 2 V i2 + β2 2 V i2 + β2 2β 2β ⎢ Vni + β 2 ⎥ ) ) ( ( √ n √ n ⎢ ⎥ ⎢ ⎥ 2 2 2 ⌢i i i i 2 i 2 v ¯ V + V + β + β n v ¯ V − V + β + β 2⌢ n ix ⎥ ⎢ ⌢i n n y n n nx ⎣ ⎦ − 2 ( 2 ) ( 2 ) i 2 2 i 2 2 i 2 Vn + β 2β Vn + β 2β Vn + β

182

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 3. Effects of the value of β on the L 2 -norm error of the solution for the Kovasznay problem for the flow variables, (a) u-velocity, (b) v-velocity, (c) pressure and (d) ∇.V .

and its inversion is ⎡

⎤ u¯ i ⌢ n iy

σ

i −1

v¯ i ⌢ n ix

− ⎢ ( ) ⎢ √ ⎢ i i2 + β2 ⎢ − V − V =⎢ n n ⎢ ( ) √ ⎣ 2 i i 2 − Vn + Vn + β

−(β 2 ⌢ n iy

+ v¯

i

−

v˜li

β 2⌢ n ix

+ u¯

β 2⌢ n ix

β 2⌢ n iy

β 2⌢ n ix

β 2⌢ n iy

The jump vector related to the primitive variables is ⎤ ⎡ p˜ri − p˜li ⎢ ⎥ Θ i = ⎣u˜ ri − u˜ li ⎦ v˜ri

Vni )

i

Vni ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

(50)

(51)

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

183

Fig. 4. Effect of the value of β on the L ∞ -norm error of the solution for the Kovasznay problem for the flow variables, (a) u-velocity, (b) v-velocity, (c) pressure and (d) ∇.V .

Then, the wave strength vector can be expressed as ai = σ

i −1

Θi

(52)

and finally, the Roe numerical flux can be calculated through the following relation: 3

∑ ⏐⏐ i ⏐⏐ ∗ ˜ i = {{ H ˜ i }} − 1 H a i ⏐λ ⏐ σ⃗¯ i 2 k=1 k k k

(53)

Different modifications and extensions for the classical AUSM have been made [27]. In this study, AUSM+-up is considered: ∗

˜i = H ˜i H

∗c

˜i +H

∗p

(54)

184

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 5. Calculated L 2 -norm error of the solution for the Kovasznay problem for the polynomial degrees, (a) N = 1, (b) N = 2, (c) N = 3 and (d) N = 4.

˜i where H given by ˜i H ˜i H

∗c

∗p

∗c

˜i and H

∗p

are the numerical convective flux and the numerical pressure flux, respectively, and they are

⎡ ⎤ 0 i ⎢ i⎥ = p¯ ⎣ˆ nx ⎦

(55)

ˆ n iy i

i = c¯i M ψ l/r

(56)

where p, ¯ c, ¯ M, the subscript (.)l/r and ψ for the ith triangle are defined as follows: − p¯ = pl P+ 5 (Ml ) + pr P5 (Mr ) + Pu

(57)

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

185

Fig. 6. Calculated L ∞ -norm error of the solution for the Kovasznay problem for the polynomial degrees, (a) N = 1, (b) N = 2, (c) N = 3 and (d) N = 4.

1 (cl + cr ) 2 − M = M+ 4 (Ml ) + M4 (Mr ) + M p { (.)l i f M ≥ 0 (.)l/r = (.)r i f M < 0 ⎡ ⎤ 1 ψ = ⎣u ⎦ v

c¯ =

± and P± 5 , Pu , M4 , M p , Ml and Mr are expressed as ⎧ 1 ± ⎨ M1 i f |M| ≥ 1 M P± = 5 ) ⎩ ±( M2 ±2 − M ∓ 16α MM∓ i f |M| < 1 2 − Pu = −2K u c¯2 P+ 5 (Ml ) P5 (Mr ) (Mr − Ml )

(58) (59) (60) (61)

(62) (63)

186

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 7. Grid used to solve the Kovasznay problem with h = 0.125 for increasing the resolution of the results, (a) illustrating the main (thick lines) and fine (thin lines) grids together and (b) the fine grid.

⎧ ⎨

M± i f |M| ≥ 1 1 ( ) ⎩ ± ∓ M2 1 ∓ 16δM2 i f |M| < 1 ( ) p −p r l M p = −K p max 1 − γ M˜ 2 , 0 c¯2 Vn Ml = l c¯ Vnr Mr = c¯

M± 4 =

(64) (65) (66) (67)

and ( ) ˜ 2 = 1 Ml 2 + Mr 2 M 2

(68)

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

187

Fig. 8. Computed u-velocity contours for the Kovasznay problem with h = 0.125 and N = 4 presented by using the (a) first and (b) second approaches.

1 (M ± |M|) 2 1 2 M± 2 = ± (M ± 1) 4

M± 1 =

(69) (70)

3 The constants are δ = 18 , K p = 14 , γ = 1, K u = 43 and α = 16 . For the viscous fluxes, different methods have been introduced in the literature such as the central, local DG (LDG) and internal penalty (IP) fluxes [32]. Here, the stabilized central flux is used and it can be defined as

q y ∗ q˜1 i = {{˜ q1 i }} − Γ u˜ i

(71)

188

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

q y ∗ q˜2 i = {{˜ q2 i }} − Γ u˜ i u˜

i∗

(72)

i

(73)

= {{u˜ }}

and Γ is a constant parameter. This parameter is added to the flux to improve both the accuracy and stability of the numerical method. However, in [32] it is shown that Γ = 0 results in accurate solutions for the compressible flows without shocks (such as the steady shear flow), so here this value is considered and the numerical results obtained show that this value is reliable for having a stable and accurate solution of the incompressible flows with the artificial compressibility method. 3.3. Boundary conditions The effects of the boundary conditions (BCs) into the simulation are made through the numerical flux terms related to the boundary elements. Therefore, the right value of the boundary edges should be modified appropriately to implement the BCs. Here, three types of the boundary conditions (BCs) are considered in the simulations, namely, the no-slip, inflow and outflow BCs. These BCs are comprised of the Dirichlet and Neumann BCs. The relations for the Dirichlet BCs are Q˜ ri = 2Q i − Q˜ li ⏐i ⏐i ∂ Q˜ ⏐⏐ ∂ Q˜ ⏐⏐ ⏐ = ⏐ ∂n ⏐ ∂n ⏐ r

(74) (75)

l

and for the Neumann BCs (∂ Q i /∂n = 0) are Q˜ ri = Q˜ li ⏐i ⏐i ∂ Q˜ ⏐⏐ ∂ Q˜ ⏐⏐ ⏐ =− ⏐ ∂n ⏐ ∂n ⏐ r

(76) (77)

l

where Q can be the pressure or u- or v-velocity variable and Q i is the exact boundary value. For the no-slip BCs at the wall boundary, the Dirichlet BCs are applied for the velocities (u i = v i = 0) and the Neumann BCs are used for the pressure (∂ pi /∂n = 0). For the inflow boundary, the Dirichlet BCs are used for the velocities (u i and v i are assumed to be equal to their freestream values) and the Neumann BCs are applied for the pressure. For the outflow boundary, the Neumann BCs are utilized for the velocities and the Dirichlet BCs are used for the pressure ( pi is assumed to be constant). 3.4. Time integration For the simulation of the unsteady incompressible flows, the implicit dual-time stepping method is applied. For this aim, a term related to the real time has to be added to Eq. (42) as follows: i ˜i ∂W ∂˜ Q i + =˜ R ∂t ∂τ

(78)

˜ i = (0, u˜ i , v˜ i ). By applying the second-order backward differencing method for the real time derivative where W terms the following equation can be obtained: m

i ˜i ˜ i − 4W ˜i + W ˜i ∂W 3W ∂˜ Q i i ˜ ˜ =R − =R − ∂τ ∂t 2∆t ′i

m−1

′ =˜ R

i

(79)

where ˜ R is the right-hand side vector which includes the real time derivative terms. Note that the accuracy of the solution in the artificial-time does not depend on the choice of the time integration method for discretizing the artificial-time derivative terms and because of simplicity of the implementation, the explicit methods are used here for this aim. In the present study, two explicit methods, namely, the third-order TVD Runge–Kutta and Euler

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

189

methods are formulated and applied for the discretization of the artificial-time derivative terms. The explicit Euler method applied to Eq. (79) is written as follows: ( in n ) in i n+1 ′i ˜ R ˜ Q , τi (80) =˜ Q + ∆τ i ˜ Q and the explicit third-order TVD Runge–Kutta scheme applied to Eq. (79) is given by ( in n ) in i (1) ′i ˜ R ˜ Q , τi (81) Q + ∆τ i ˜ Q =˜ ] [ ) ( (1) n 1 ˜in (1) i (2) ′i ˜ Q = (82) 3 Q + Q i + ∆τ i ˜ R Q i , τ i + ∆τ i 4 ( [ )] 1 ˜in 1 n (2) (2) i n+1 ′i ˜ Q = Q i , τ i + ∆τ i Q + 2 Q i + 2∆τ i ˜ R (83) 3 2 where the superscripts n and m in the above relations indicate the iteration stage connected to the discretization of the artificial and real time derivative terms, respectively. The local time-stepping procedure is used for calculating ∆τ i as ( )( i ) 2 rD i ∆τ = C F L min ∆r (84) 3 λmax i where ∆r denotes the spacing between the grid points in the computational space and in the one-dimensional space, r D is the characteristic length ⏐√ of the triangle√(the area⏐ of the triangle divided by the radius of the inscribed ⏐ ⏐ circle is used here) and λmax = ⏐ β 2 + u 2 + v 2 + u 2 + v 2 ⏐ is the maximum eigenvalue in an element. Here, max

the Courant–Friedrichs–Lewy (C F L) number is selected to be C F L = 0.1 for all the simulations. In the local time-stepping procedure, each element has a distinct time and time-step and such a procedure helps to converge the solution with fewer time steps. The study has indicated that both the explicit time integration methods used for the discretization of the artificialtime derivative terms in the formulation of the NDGM with the artificial compressibility method are stable and they provide identical results. Although the third-order TVD Runge–Kutta method gives a larger CFL number, the computation time of this time integration method is higher than the explicit Euler method. Therefore, all the computations performed in the present study are based on the explicit Euler method for the time integration in the artificial time. 4. Numerical results Herein, several test cases are simulated to assess the accuracy and robustness of the developed incompressible flow solver with the artificial compressibility method based on the NDGM. At first, the Kovasznay flow problem is solved to investigate the order of accuracy and performance of the solution method. For this aim, the polynomial degrees up to four are considered. The effects of three numerical fluxes used on the accuracy and performance of the NDGM is also investigated for this test case. The accuracy of the solution method based on the NDGM involving different boundary conditions are examined for the flow over the backward facing step problem. Then, to show the accuracy and capability of the solution methodology for computing the steady and unsteady flows, the NACA 0012 airfoil and the circular cylinder are simulated. Finally, as a complicated problem the flow over two side-by-side circular cylinders is computed to demonstrate the capability of the solution method based on the NDGM in handling complex geometries. The effects of the polynomial degrees on the solution are also examined for the problems simulated. The programming language of the developed code is based on C++, the unstructured grids are based on the Delaunay triangulation and the post-processing of the results is performed by the Tecplot software. 4.1. Kovasznay flow problem The Kovasznay flow problem has the analytical solution and thus it is appropriate for assessing the accuracy of the solution of the NDGM on unstructured grids. In this problem, all the boundaries are Dirichlet and therefore only the effects of the spatial discretization will be examined. The steady-state solution for this problem is as follows [1]: u (x, y) = 1 − exp(Φx) cos(2π y)

(85)

190

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214 Table 1 L 2 and L ∞ -norm errors and convergence rate of the u-velocity for the Kovasznay problem solved by using the local LF flux. N

h = 0.5

h = 0.25

h = 0.125

Convergence rate

6.66E−02 2.14E−02 1.74E−03 2.58E−04

1.81E−02 1.43E−03 8.55E−05 4.30E−06

1.9 3.1 4.6 5.1

2.20E−01 7.38E−02 1.08E−02 2.21E−03

6.59E−02 9.43E−03 7.87E−04 4.60E−05

1.5 2.6 3.8 4.6

L 2 -norm 1 2 3 4

2.52E−01 1.07E−01 5.01E−02 5.37E−03

1 2 3 4

5.39E−01 3.57E−01 1.45E−01 2.64E−02

L ∞ -norm

Table 2 L 2 and L ∞ -norm errors and convergence rate of the u-velocity for the Kovasznay problem solved by using the Roe flux. N

h = 0.5

1 2 3 4

2.73E−01 1.26E−01 5.50E−02 6.30E−03

h = 0.25

h = 0.125

Convergence rate

8.41E−02 2.52E−02 2.01E−03 3.34E−04

2.36E−02 1.75E−03 1.02E−04 5.18E−06

1.8 3.1 4.5 5.1

3.52E−01 8.14E−02 1.23E−02 2.70E−03

9.55E−02 1.18E−02 8.01E−04 5.44E−05

1.5 2.5 3.9 4.6

L 2 -norm

L ∞ -norm 1 2 3 4

7.37E−01 4.05E−01 1.82E−01 3.26E−02

λ exp(Φx) sin(2π y) (86) 2π 1 (87) p (x) = − exp (2Φx) + C 2 where C is an arbitrary parameter and here it sets to be two and the parameter Φ depends on Re as √ Re Re2 Φ= − + 4π 2 (88) 2 4 and Re = 40 is used here. The purpose of computing the current test case is to study the performance of the NDGM. For this aim, the grid used for the analysis of the convergence rate is shown in Fig. 1 where different grids with three characteristic lengths, i.e., h = 0.5, h = 0.25 and h = 0.125 are used. The errors based on the local LF, Roe and AUSM+-up numerical fluxes related to these three grids and different polynomial degrees are presented in Tables 1–3, respectively. Note that the L 2 - and L ∞ -norm errors are calculated based on the u-velocity component computed by the NDGM in comparison with the analytical solution. It is observed that the expected convergence rate is obtained for the NDGM with all the three numerical fluxes. In [32], this is mentioned that the optimal convergence rate is typically close to O(h N +1/2 ). In Fig. 2, the spectral convergence of the NDGM is illustrated on the triangle elements. It is also indicated that increasing the polynomial degree results in an exponential decay of the error. The effects of the value of artificial sound speed, 1 ≤ β ≤ 27, on the accuracy of the solution are examined and the corresponding results are illustrated in Figs. 3 and 4. Figure 3 gives the L 2 -norm errors obtained by considering the calculated velocity components, pressure and divergence of the velocity compared to their exact ones and it is observed that the errors are reduced by employing higher polynomial degrees. The study shows that the parameter β v (x, y) =

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

191

Fig. 9. Comparison of the u-velocity contours for the Kovasznay flow problem with different polynomial degrees. Table 3 L 2 and L ∞ -norm errors and convergence rate of the u-velocity for the Kovasznay problem solved by using the AUSM+-up flux. N

h = 0.5

1 2 3 4

2.48E−01 1.30E−01 5.92E−02 6.37E−03

h = 0.25

h = 0.125

Convergence rate

8.41E−02 2.65E−02 2.09E−03 4.06E−04

2.34E−02 1.77E−03 1.07E−04 5.61E−06

1.7 3.1 4.6 5.1

3.52E−01 1.30E−01 1.22E−02 3.86E−03

1.05E−01 1.10E−02 1.05E−03 6.81E−05

1.4 2.8 3.9 4.4

L 2 -norm

L ∞ -norm 1 2 3 4

7.39E−01 5.50E−01 2.20E−01 3.23E−02

does not significantly affect the accuracy of the solution. The errors calculated based on the velocities are relatively independent with respect to the value of β and a higher β results in a lower error based on the divergence of the velocity and a higher error based on the pressure. The effects of the value of β on the L 2 - and L ∞ -norm errors are similar (see Fig. 4), unless the errors calculated based on the L ∞ -norm are higher than those based on the L 2 -norm. Now, the spatial convergence rate of the solution of the nodal discontinuous Galerkin method (NDGM) is examined by comparison of the numerical results for the u-velocity, v-velocity, pressure and divergence of the velocity with the analytical ones and the L 2 - and L ∞ -norm errors for the different polynomial degrees are given in Figs. 5 and 6. Here, the local Lax–Friedrich is used as the numerical flux. It can be seen from these results that the expected spatial convergence rate is obtained for the NDGM. Again, the L ∞ -norm errors are higher than the L 2 -norm ones. Here, two approaches for illustrating the numerical results are investigated. The first approach is to use the values in the center of the main grid (Fig. 1(c)) and the second approach is to use FEM’s property of the NDGM in conjunction with the generation of the sub-cells to improve the resolution of the results (Fig. 7). In the second approach, each cell is divided into 16 sub-cells (Fig. 7(a)) that results in a finer grid (Fig. 7(b)), then the values are interpolated into this finer grid. The differences between these two approaches are illustrated in Fig. 8. It is shown that the resolution of Fig. 8(b) related to the second approach is better than Fig. 8(a) related to the first approach. Note that in the high-order FD or FV methods values of the variables between the solution points are undetermined

192

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 10. Contours of the (a) u-velocity, (b) pressure, (c) v-velocity and (d) the streamlines for the Kovasznay problem.

while in the FE or NDG methods it is available analytically and this is an important feature of the NDGM. In Fig. 9, the contours of the u-velocity are shown for different polynomial degrees. In this figure, the second approach is used for illustrating the results and the main grid characteristic length is h = 0.125. It is shown that N = 1 results in a linear polynomial and this is definitely obvious in this figure. Moreover, FV’s property of the NDGM can be shown because there are jumps at the cell interface (the solution between the elements can be discontinuous) and it can be observed clearly for N = 1. The study shows that the results obtained with N = 4 are more accurate in comparison with other polynomial degrees and they are close to the exact solution. In fact, increasing the polynomial degrees results in a more accurate solution. Note that all the subsequent results are based on the second approach unless otherwise mentioned. The computed flow field shown by the contours and streamlines are presented in Fig. 10 that are based on the grid with h = 0.125 and N = 4. In Fig. 11(a) convergence history of the NDGM based on the local LF, Roe and AUSM+-up numerical fluxes is presented. In this figure, the grid with h = 0.125 and N = 4 are used. It is shown that the type of the numerical flux used does not significantly affect the convergence history, so the effect of the numerical fluxes used on the

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

193

Fig. 11. Convergence history of the solution for the Kovasznay problem for (a) different numerical fluxes (h = 0.125 and N = 4), (b) different grids (the Roe numerical flux and N = 4) and (c) different polynomial degrees (h = 0.125 and the Roe numerical flux).

194

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 12. (a) Geometry including defined important flow parameters and (b) the grid used to solve the backward facing step problem.

Table 4 Computational time (s) of the solution for the Kovasznay problem. h = 0.5 h = 0.25 h = 0.125

N =1

N =2

N =3

N =4

714.8 1117.8 2689.4

848.8 1350.4 3264.4

989.2 1680.4 3887.4

1182 2074.2 5038

convergence history is negligible. The effects of the grid size on the convergence history is illustrated in Fig. 11(b). The study shows that a finer grid requires more iterations to reach the desired error criteria. Moreover, the effect of the polynomial degree on the convergence history is given in Fig. 11(c) and it is shown that the convergence rate is decreased by increasing the polynomial degree. Since the local LF numerical flux provides accurate results with a lower computational cost compared to the Roe and AUSM+-up numerical fluxes, it is applied in the subsequent simulations. As shown in Table 4, increasing the polynomial degree results in an increase of the computational cost. Since in relatively fine grids fourth-order solutions (N = 3) give relatively accurate results with lower computational cost compared to the fifth-order solutions (N = 4), thus in the next test cases, the simulations are performed with N = 3, unless otherwise indicated.

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

195

Fig. 13. Computed flow field shown by the streamlines for the backward facing step problem at (a) Re = 200, (a) Re = 400, (a) Re = 600 and (a) Re = 800.

4.2. Flow over backward facing step This problem has a relatively complex flow field and thus it is suitable for the assessment of high-order accurate methods which should deal with different types of the boundary conditions and also have to resolve the flow field features accurately. In addition, the computation of the incompressible flow over the backward facing step problem examines the accuracy and capability of the solution method in computing the internal flows. The geometry of this problem and important flow field characteristics are shown in Fig. 12(a). in which the step height, the total height of the channel, the primary vortex reattachment point, the secondary vortex separation point, the secondary vortex reattachment point and the length of the channel are denoted by h, H , L 1 , L 2 , L 3 and L, respectively. At the inflow,

196

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 14. Comparison of the wall shear distribution for the (a) lower and (b) upper walls of the backward facing step problem at different Reynolds numbers.

a parabolic velocity distribution is given as follows: (H − y) (y − h) u = 6u 0 (89) (H − h)2 where the constant parameters are u 0 = 1, h = 0.5, H = 1 and L = 25. The present computations based on the NDGM are performed for the backward facing step for different Reynolds numbers and the comparisons of the results obtained with the available numerical results are made. A grid with 4038 triangle elements, as shown in Fig. 12(b), is used for the flow simulations. The computed results shown by the streamlines are illustrated in Fig. 13. It is shown that increasing the value of the Reynolds number results in the more complex flow field in such a way that the recirculation zones become larger. The distribution of the vorticity along the lower and upper walls of the channel are shown in Fig. 14(a) and (b), respectively. These numerical results

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

197

Fig. 15. Effect of the polynomial degree on the wall shear distribution for the lower and upper walls of the backward facing step problem at Re = 600.

Fig. 16. Convergence history of the solution based on the L 2 -norm error of the u-velocity for the backward facing step problem at different Reynolds numbers.

are compared with those of the compact finite-difference method [28] which show excellent agreement. One of the advantages of the NDGM is that it allows calculating the flow variables in an element everywhere because it has the properties of both the finite element method (FEM) and the finite volume method (FVM) together. Therefore, here the vorticity at the center of the boundary face of the boundary elements is calculated although there is not any point on the center of the face of each triangle when N = 3 is used for the simulation. The common way for determining the position of the separation and reattached points is to check where the gradient of the tangential velocity (u) with respect to the normal direction of the wall (y) becomes zero. These values can be obtained from Fig. 14 and they are presented in Table 5, where excellent agreement is observed between the numerical results obtained by the NDGM and those of Refs. [25,28]. The effect of the polynomial degree on the accuracy of the solution for the flow at Re = 600 is illustrated in Fig. 15. It is shown that increasing the polynomial degree from

198

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 17. Grid with 5600 elements and the minimum edge length 0.01 used to solve the flow over the NACA 0012 airfoil. Table 5 Comparison of the lengths of the primary and secondary vortices for the backward facing step problem. Re

L1 NDGM

Ref. [28]

Ref. [32]

NDGM

Ref. [28]

Ref. [32]

NDGM

Ref. [28]

Ref. [32]

200 400 600 800

2.6694 4.3234 5.3731 6.0906

2.6660 4.3221 5.3725 6.0953

2.6634 4.3312 5.3963 6.1048

– 4.0022 4.3771 4.8593

– 4.0084 4.3748 4.8530

– 4.0512 4.4037 4.8623

– 5.2027 8.1156 10.4833

– 5.1865 8.1070 10.4733

– 5.1461 8.0949 10.4739

L2

L3

N = 1 to 2 considerably improves the accuracy of the results whereas it is slighter from N = 2 to 3. The results presented for the backward facing step flow show that the NDGM with N = 3 gives more accurate results. The convergence history of the solution for different Reynolds numbers is illustrated in Fig. 16 and it is shown that more iterations are required for the convergence of the solution at higher Reynolds numbers. As shown in Fig. 12, the grid used is relatively coarse, especially in the y-direction. Note that calculating the flow features near and on the walls such as the shear stress is very important for the determination of the separation and reattachment points, especially for the secondary vortex. The reattachment point of the secondary vortex at Re = 600 for N = 1, 2 and 3 are 7.3747, 8.0868 and 8.1156, respectively, which shows increasing the polynomial degree results in a more accurate solution and it facilitates the use of coarser grids. This investigation demonstrates that the NDGM can be used effectively to simulate the incompressible flow problems based on the artificial compressibility formulation considering different types of boundary conditions and it can provide reasonable results even on coarse grids. 4.3. Flow over the NACA0012 airfoil In the previous problem, the solution domain was straight and simple, however, in the current problem the airfoil shape causes the solution domain to be more complex than the previous test case. The separated and non-symmetric flow over the airfoil make this problem appropriate for a better assessment of the accuracy and robustness of the NDGM. The grid used for this problem is shown in Fig. 17. Here, a grid with 5600 triangle elements and minimum length edge 0.01 is used for the numerical simulation and the outer domain diameter is 20c where c is the airfoil chord. The initial condition for the flow variables is considered to be the freestream condition. The computed flow

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

199

Fig. 18. Computed flow field shown by (a) the u-velocity contours and (b) the streamlines for the flow over the NACA 0012 airfoil at Re = 500 and AO A = 10◦ .

200

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 19. Comparison of the surface pressure coefficient distribution for the flow over the NACA 0012 airfoil at Re = 500 and AO A = 10◦ .

field over the NACA 0012 airfoil at the Re = 500 and the angle of attack AO A = 10◦ shown by the contour lines of the u-velocity and the streamlines are illustrated at Fig. 18(a) and (b), respectively. As shown in this figure, the flow separates at the upper surface (the suction side) of the airfoil and reattaches near the trailing edge. In Fig. 19, the distribution of the pressure coefficient over the airfoil surface obtained by different polynomial orders is illustrated and the numerical results are compared with those obtained by Hafez [24] which show good agreement. It is shown that increasing the polynomial degree results in a more accurate solution. In Fig. 20, the flow field features shown by the contours of the v-velocity for some elements near the leading edge of the airfoil are illustrated. It can be seen in this figure that the solution points will be 3, 6 and 10 for N = 1, 2 and 3, respectively. As indicated in Fig. 20, N = 1 results in a linear distribution for the v-velocity in each element whereas for N = 2 and 3 it will not be linear. It can also be seen that the v-velocity goes to zero near the airfoil surface to satisfy the no-slip boundary conditions. The flow over the NACA0012 airfoil at Re = 800 and AO A = 20◦ is now simulated and this condition has a complex unsteady flow feature. This simulation is performed by using the real time step size ∆t = 0.1. The flow field indicated by the streamlines are illustrated in Fig. 21 at t = 5. As obvious in this figure, there are two circulating flow zones on the suction side of the airfoil, each one with the distinctive separation and reattachment points. For the first circulating flow zone, the flow separates at x = 0.063 and reattachs at x = 0.419. The separation of the flow close to the leading edge is an attribute of the stall which occurs at the high angle of attack considered here. For the second circulating flow zone, the separation and reattachment points take place at x = 0.544 and 0.829. In addition, there is a circulating flow zone separated from the airfoil body which is an indication of the vortex shedding. In Fig. 22, the time-variation of the v-velocity at (x, y) = (1.1, 0) computed by the NDGM is compared with the other available numerical results [24,26]. As shown in this figure, the agreement between the numerical results is satisfactory, so the NDGM with the implicit dual-time stepping method can be used to accurately simulate the unsteady flows over the curved boundaries with relatively complex flow features. . 4.4. Flow over the circular cylinder Now, both the steady and unsteady flows over the circular cylinder are simulated. For the steady flow, the Reynolds numbers are 20 and 40 and for the unsteady one Re = 100 is considered. A grid with 2934 elements and the minimum edge length 0.05 is generated, as shown in Fig. 23. The far-field boundary is supposed as a circle with a diameter 20D where D = 1 is the diameter of the circular cylinder. The streamlines and pressure contours of the flow over the circular cylinder computed by the NDGM at Re = 20 and 40 are shown in Fig. 24. This is

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

201

Fig. 20. v-velocity contours and the solution points for different polynomial degrees for some elements located on the leading edge of the NACA 0012 airfoil at Re = 500 and AO A = 10◦ .

clear in this figure that the higher Reynolds number causes a bigger circulating flow zone behind the cylinder at the steady-state condition. The comparison of the surface pressure coefficient distribution obtained by the NDGM at Re = 40 with the numerical [29] and experimental [23] results is made in Fig. 25 which shows a good agreement. To further evaluate the numerical results obtained by the NDGM, the drag coefficient Cd , the pressure coefficient at the aft (θ = 0) and rear (θ = π ) stagnation points, the length of the circulating flow zone normalized by the radius of the cylinder 2L/D and the separation angle θs are compared with the available numerical and experimental results [12,14,26,29,46,60], as listed in Table 6, which exhibits good agreement. The convergence history of the solution for this problem at Re = 20 and 40 is illustrated in Fig. 26 and it is shown that the solution method at a higher Reynolds number requires more iterations to reach the steady-state condition with a prescribed convergence criteria. The unsteady flow over the circular cylinder is now simulated at Re = 100 and the simulation is performed by utilizing the real time step size ∆t = 0.125. The streamlines of the flow field obtained by the NDGM for a time period T are illustrated in Fig. 27. The time-variation of the lift and drag coefficients is shown in Fig. 28 which depicts the mean drag and lift coefficients of 1.345 and 0.325, respectively, and the Strouhal number

202

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 21. Computed flow field shown by the streamlines for the flow over the NACA 0012 airfoil at Re = 800, AO A = 20◦ and t = 5.

Fig. 22. Comparison of the time history of the v-velocity at x = 1.1 for the flow over the NACA 0012 airfoil at Re = 800 and AO A = 20◦ .

St =

fD u0

= 0.165 where f is the frequency of the vortex shedding. The comparison of the present results with the

available numerical results [6,7,26,29,41,55] is made in Table 7 and the agreement is reasonable. The convergence history of the solution method by applying the implicit dual-time stepping method for this unsteady problem is shown in Fig. 29. It can be observed that for the first time steps more iterations are needed for the convergence rate of the inner loops with the convergence criteria considered (here, 10−6 ) and the convergence rate is faster for the last time steps. This flow problem with a relatively complex flow feature shows that the solution method proposed based on the NDGM with the artificial compressibility approach and the implicit dual-time stepping method can be used for accurately simulating unsteady incompressible flows.

203

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 23. Grid with 2834 elements and the minimum edge length 0.05 used to solve the flow over the circular cylinder. Table 6 Comparison of the flow field parameters for the circular cylinder problem. Re

Author

Method

Cd

−C p (0)

C p (π)

2L/D

θs

20

Present solution Hejranfar and Parseh Hejranfar and Ezzatneshan Dennis and Chang Nieuwstadt and Keller Coutanceau and Bouard Tritton

N-S N-S CFDLBM N-S N-S Experiment Experiment

2.071 2.025 2.021 2.045 2.053 – 2.2

0.5419 0.5449 0.5465 0.589 0.582 – –

1.3442 1.2938 1.2659 1.269 1.274 – –

1.843 1.843 1.848 1.88 1.786 1.86 –

43.86 43.62 43.58 43.7 43.37 44.8 –

40

Present solution Hejranfar and Parseh Hejranfar and Ezzatneshan Dennis and Chang Nieuwstadt and Keller Coutanceau and Bouard Tritton

N-S N-S CFDLBM N-S N-S Experiment Experiment

1.545 1.514 1.515 1.522 1.550 – 1.65

0.486 0.477 0.481 0.509 0.554 – –

1.203 1.166 1.154 1.144 1.117 – –

4.443 4.537 4.510 4.690 4.357 4.26 –

53.89 53.64 51.86 53.80 53.34 53.50 –

Table 7 Comparison of the drag and lift coefficients and the Strouhal number for the circular cylinder problem at Re = 100. Author

Cd

Cl

St

Present solution Hejranfar and Parseh Hejranfar and Ezzatneshan Liu et al. Calhoun Russell and Wang Chiu et al.

1.345 ± 0.0096 1.33 ± 0.0097 1.336 ± 0.0114 1.35 ± 0.012 1.35 ± 0.014 1.38 ± 0.007 1.35 ± 0.012

±0.325 ±0.3255 ±0.3354 ±0.339 ±0.3 ±0.322 ±0.303

0.165 0.164 0.163 0.164 0.175 0.169 0.167

4.5. Flow over two side-by-side circular cylinders The steady and unsteady flows over two side-by-side circular cylinders at several gap sizes gs are now simulated in this section. The gap size gs measures the distance from the center to center of the cylinders and this parameter is

204

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 24. Computed flow field shown by the pressure contours and the streamlines for the flow over the circular cylinder at (a) Re = 20 and (b) Re = 40.

non-dimensionalized based on the cylinder diameters (see Fig. 30(a)). The flow over the two side-by-side cylinders, especially at relatively small gap sizes results in a large circulating flow zone which highlights the ability of the high order accurate numerical methods such as the NDGM in handling the complex geometries and resolving the complex flow field features. Here, the diameter of the cylinders is D = 1 and different gap sizes are considered which are gs = 1.05, 1.1, 1.2, 1.4, 1.5, 1.75, 2, 2.5 and 3. The generated grids for this problem with different gap sizes are shown in Fig. 30(b–j). The generated grids have the minimum edge length size 0.025 and the number of elements increases from 5713 for gs = 1.05 to 7042 for gs = 3 and the far-field boundary is defined by a circle with the radius 50. As mentioned in the Kovasznay problem, one of the most important properties of the NDGM is that it has both the features of the FVM and FEM. The FEM’s feature of the NDGM let us to calculate the variables in an element wherever it is required. For example, if the results obtained by N = 3 are considered then for showing the streamlines, as mentioned in the Kovasznay problem, two approaches are used. These two approaches are shown in Fig. 31 for the two side-by-side circular cylinder problem with Re = 20 at the gap size 1.05. In the first approach

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

205

Fig. 25. Comparison of the surface pressure coefficient distribution for the flow over the circular cylinder at Re = 40.

Fig. 26. Convergence history of the solution for the flow over the circular cylinder at Re = 20 and 40.

shown in Fig. 31(a), the value of the variables at the center of the elements (the main grid) are used for showing the streamlines, however, in the second approach shown in Fig. 31(b), the sixteen triangular elements are generated in each element (the finer grid) and the interpolated variables in the new finer grid are used to obtain the streamlines. It can be seen in this figure that by the first approach the position of the rear stagnation point of the vortex take places at x = 5.75, however for the second approach this is at x = 6.95 which is in a good agreement with x = 6.85 obtained by Vakil and Green [61]. Note that with using the first approach (this approach is very common in CFD works) the position of the rear stagnation point of the vortex for both N = 1 and 3 will be similar and it can be misleading. Therefore, in the NDGM special attention should be paid for the post-processing because by choosing an inconvenient approach for presenting the numerical results, the valuable information obtained by applying the high-order polynomials may be spoiled. All the steady-flow simulations are performed at Re = 20 and the numerical results based on the NDGM for the pressure field and the streamlines are illustrated in Fig. 32. It is observed by increasing the gap size from 1.05

206

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 27. Computed flow field shown by the streamlines for the flow over the circular cylinder at Re = 100 and different times.

to 1.5, the distance along the x-axis between the center of the vortices and cylinders increases and at gs = 1.75 there is no vortex behind the cylinders. For the gap sizes from 2 to 3, there are the vortices attached on the cylinder bodies and the vortices for gs = 2.5 and 3 are larger than the case gs = 2. The numerical results for the drag, lift, viscous drag and pressure drag coefficients obtained by the NDGM are compared with those reported by Vakil and Green [61], as shown in Fig. 33 and given in Table 8, which show good agreement. To further examine the accuracy the solution method, the unsteady flow over the two side-by-side circular cylinders at Re = 100 and gs = 3.0 is simulated. The unstructured grid for this problem is the same as the steady case described before and the unsteady computations are performed by using the real time step size ∆t = 0.125. The time-variation of the drag and lift coefficients for the upper cylinder are shown in Fig. 34 and the mean and variation of the drag and lift coefficients and the Strouhal number are listed and compared with the available numerical results [15,39,56], as given in Table 9, which indicates the agreement is satisfactory. The computed flow field shown by the streamlines for a time period is given in Fig. 35 and it is shown that the behavior of the vortices on the surface of the upper and lower circular cylinders is anti-phase [15]. Indications are that the complex unsteady flow field structures are reasonably resolved by applying the NDGM with the artificial compressibility approach and the implicit dual-time stepping method and thus the solution methodology proposed is capable of accurately simulating the incompressible flows over the complicated problems.

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

207

Fig. 28. Time history of the drag and lift coefficients for the flow over the circular cylinder at Re = 100.

Fig. 29. Convergence history of the solution by applying the dual-time stepping method for the flow over the circular cylinder at Re = 100.

5. Conclusion A high-order nodal discontinuous Galerkin method (NDGM) is developed and applied for computing the steady and unsteady incompressible flows. For this aim, the two-dimensional incompressible Navier–Stokes equations with the artificial compressibility approach and the implicit dual-time stepping method are considered and the resulting system of equations is solved using the NDGM with different polynomial degrees on triangle elements. Three numerical fluxes are formulated and considered to be applied in the NDGM and their effects on the accuracy and convergence rate of the solution are also investigated. Different incompressible flow problems are computed to assess the accuracy and robustness of the NDGM. The conclusions and remarks of the present study are as follows: 1. To properly account coupling between the continuity and momentum equations and so that the continuity equation to be divergence free, the artificial compressibility method is considered in the formulation for the incompressible flow simulations by applying by the NDGM. The solution methodology proposed by applying

208

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 30. (a) Geometry, and the grid used to solve the flow over the two side-by-side circular cylinders with the gap sizes, (b) gs = 1.05, (c) gs = 1.1, (d) gs = 1.2, (e) gs = 1.4, (f) gs = 1.5, (g) gs = 1.75, (h) gs = 2.0, (i) gs = 2.5 and (j) gs = 3.0.

the implicit dual-time stepping method is capable of accurately computing the unsteady flows. Note that the solution procedure proposed here can be considered as a high-order accurate truly incompressible flow solver.

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

209

Fig. 31. Computed flow field shown by the streamlines for the flow over the two side-by-side circular cylinders at Re = 20 with the gap size gs = 1.05 for (a) the main grid (the first approach) and (b) interpolated into the fine grid (the second approach). Table 8 Computed drag, lift, viscous drag and pressure drag coefficients for the two side-by-side cylinders problem at Re = 100. Gap size

Cd

Cl

Viscous drag coefficient

Pressure drag coefficient

1.05 1.10 1.20 1.40 1.50 1.75 2.00 2.50 3.00

1.614 1.644 1.707 1.851 1.928 2.093 2.176 2.216 2.209

1.378 1.383 1.379 1.271 1.154 0.792 0.562 0.379 0.289

0.393 0.436 0.513 0.641 0.695 0.794 0.837 0.860 0.860

1.221 1.207 1.194 1.211 1.233 1.299 1.338 1.355 1.347

210

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 32. Computed flow field shown by the streamlines and 11 equally spaced pressure contours from 0.5 to 1.5 for the two side-by-side circular cylinders at Re = 20 with the gap sizes, (a) gs = 1.05, (b) gs = 1.1, (c) gs = 1.2, (d) gs = 1.4, (e) gs = 1.5, (f) gs = 1.75, (g) gs = 2.0, (h) gs = 2.5 and (i) gs = 3.0.

Fig. 33. Comparison of the drag, lift, viscous drag and pressure drag coefficients for the two side-by-side circular cylinders at Re = 20 with different gap sizes.

2. The study indicates that the results obtained by applying the NDGM for the different incompressible flow problems exhibit good agreement with the analytical and other high-order accurate numerical solutions. The solution method applied can properly model the physics of the incompressible flows, as indicated from the numerical results obtained and it can be used for an accurate and efficient computation of the incompressible flows with complex flow features.

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

211

Fig. 34. Time history of the drag and lift coefficients of the upper circular cylinder for the flow over the two side-by-side circular cylinders at Re = 100 with the gap size gs = 3.0. Table 9 Comparison of the drag and lift coefficients for the upper cylinder in the two side-by-side circular cylinders problem at Re = 100. Author

Cd

Cl

St

Present solution Liu et al. Sarvghad–Moghadam and Nooredin Ding et al.

1.48 ± 0.035 1.53 ± 0.038 1.49 ± 0.020

0.12 ± 0.45 0.12 ± 0.49 0.13 ± 0.05

0.179 0.179 0.183

1.56 ± 0.038

0.13 ± 0.25

0.182

3. In this study, the LF, Roe and AUSM+-up numerical fluxes are used in the solution algorithm of the NDGM and the accuracy and performance and these three numerical fluxes are compared with each other. It is indicated that these three numerical fluxes provide similar results and their convergence history are also nearly the same. The computational cost of these three numerical fluxes in conjunction with the NDGM is also similar because the most contribution of the computational cost comes from the NDGM discretization of the system of equations, not from the numerical flux calculations. 4. One of the main advantages of the solution algorithm based on the NDGM is that it can be implemented on unstructured grids with different orders of accuracy for reasonably computing the flow over complex geometries. The flow over the two side-by-side cylinders, as a complicated problem, is successfully simulated by the NDGM on the triangle elements and it highlights the capability of the NDGM in handling complex geometries. 5. The FEM’s property of the NDGM allows improving the resolution of the results. By the generating sub-cells in the main grid a finer grid can be obtained and then the pre-described polynomial will give the values in the finer grid. This method called the second approach and the numerical results obtained show that this approach considerably improves the resolution of the results. Therefore, in the NDGM a great attention should be paid for the post-processing because by choosing an improper approach for presenting the results, the valuable information obtained by implementing the high-order polynomials may be spoil. 6. It is indicated that the solution methodology based on the NDGM to solve the incompressible Navier–Stokes equations with the artificial compressibility method is accurate and robust in computing the steady and unsteady incompressible flows. The high accuracy solutions obtained by the application of the NDGM on unstructured grids can be considered as benchmark solutions to evaluate the other unstructured incompressible flow solvers.

212

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

Fig. 35. Computed flow field shown by the streamlines for the flow over the two side-by-side circular cylinders at Re = 100 with the gap size gs = 3.0.

Acknowledgment The authors would like to thank Sharif University of Technology for supporting this research. References [1] F. Bassi, A. Crivellini, D.A. Di Pietro, S. Rebay, An artificial compressibility flux for the discontinuous Galerkin solution of the incompressible Navier–Stokes equations, J. Comput. Phys. 218 (2) (2006) 794–815. [2] F. Bassi, A. Crivellini, D.A. Di Pietro, S. Rebay, An implicit high-order discontinuous Galerkin method for steady and unsteady incompressible flows, Comput. & Fluids 36 (10) (2007) 1529–1546. [3] F. Bassi, S. Rebay, GMRES Discontinuous Galerkin solution of the compressible Navier–Stokes equations, Discontin. Galerkin Methods (2000) 197–208. [4] O. Bou Matar, P.-Y. Guerder, Y. Li, B. Vandewoestyne, K. Van Den Abeele, A nodal discontinuous Galerkin finite element method for nonlinear elastic wave propagation, J. Acoust. Soc. Am. 131 (5) (2012) 3650–3663. [5] A. Burbeau, P. Sagaut, C.-H. Bruneau, A problem-independent limiter for high-order Runge–Kutta discontinuous Galerkin methods, J. Comput. Phys. 169 (1) (2001) 111–150. [6] D. Calhoun, A cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions, J. Comput. Phys. 176 (2) (2002) 231–275.

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

213

[7] P.H. Chiu, R.K. Lin, T.W.H. Sheu, A differentially interpolated direct forcing immersed boundary method for predicting incompressible Navier–Stokes equations in time-varying complex geometries, J. Comput. Phys. 229 (12) (2010) 4476–4500. [8] A.J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2 (1) (1967) 12–26. [9] K.A. Cliffe, E.J.C. Hall, P. Houston, Hp-adaptive discontinuous Galerkin methods for bifurcation phenomena in open flows, Comput. Math. Appl. 67 (4) (2014) 796–806. [10] B. Cockburn, C.-W. Shu, TVB RUnge–Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework, Math. Comp. 52 (186) (1989) 411. [11] B. Cockburn, C.-W. Shu, The Runge–Kutta discontinuous Galerkin method for conservation laws V, J. Comput. Phys. 141 (2) (1998) 199–224. [12] M. Coutanceau, R. Bouard, Experimental determination of the main features of the viscous flow in the wake of a circular cylinder in uniform translation. Part 1. Steady flow, J. Fluid Mech. 79 (02) (1977) 231. [13] A. Crivellini, V. D’Alessandro, F. Bassi, Assessment of a high-order discontinuous Galerkin method for incompressible three-dimensional Navier–Stokes equations: Benchmark results for the flow past a sphere up to Re=500, Comput. & Fluids 86 (2013) 442–458. [14] S.C.R. Dennis, G.-Z. Chang, Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100, J. Fluid Mech. 42 (03) (1970) 471. [15] H. Ding, C. Shu, K.S. Yeo, D. Xu, Numerical simulation of flows around two circular cylinders by mesh-free least square-based finite difference methods, Internat. J. Numer. Methods Fluids 53 (2) (2006) 305–332. [16] F. Fambri, M. Dumbser, Spectral semi-implicit and space–time discontinuous Galerkin methods for the incompressible Navier–Stokes equations on staggered Cartesian grids, Appl. Numer. Math. 110 (2016) 41–74. [17] E. Ferrer, D. Moxey, R.H.J. Willden, S.J. Sherwin, Stability of projection methods for incompressible flows using high order pressure-velocity pairs of same degree: Continuous and discontinuous Galerkin formulations, Commun. Comput. Phys. 16 (03) (2014) 817–840. [18] R. Gandham, D. Medina, T. Warburton, GPU Accelerated discontinuous Galerkin methods for shallow water equations, Commun. Comput. Phys. 18 (01) (2015) 37–64. [19] A. Gilmanov, F. Sotiropoulos, A hybrid Cartesian/immersed boundary method for simulating flows with 3D, geometrically complex, moving bodies, J. Comput. Phys. 207 (2) (2005) 457–492. [20] G. Giorgiani, S. Fernández-Méndez, A. Huerta, Hybridizable discontinuous Galerkin with degree adaptivity for the incompressible Navier–Stokes equations, Comput. & Fluids 98 (2014) 196–208. [21] F.X. Giraldo, J.S. Hesthaven, T. Warburton, Nodal high-order discontinuous Galerkin methods for the spherical shallow water equations, J. Comput. Phys. 181 (2) (2002) 499–525. [22] F.X. Giraldo, M. Restelli, A study of spectral element and discontinuous Galerkin methods for the Navier–Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases, J. Comput. Phys. 227 (8) (2008) 3849–3877. [23] A.S. Grove, F.H. Shair, E.E. Petersen, An experimental investigation of the steady separated flow past a circular cylinder, J. Fluid Mech. 19 (01) (1964) 60. [24] M. Hafez, A. Shatalov, M. Nakajima, Improved numerical simulations of incompressible flows based on viscous/inviscid interaction procedures, Comput. Fluid Dyn. 2006 (2009) 309–314. [25] K. Hejranfar, E. Ezzatneshan, A high-order compact finite-difference lattice Boltzmann method for simulation of steady and unsteady incompressible flows, Internat. J. Numer. Methods Fluids 75 (10) (2014) 713–746. [26] K. Hejranfar, E. Ezzatneshan, Implementation of a high-order compact finite-difference lattice Boltzmann method in generalized curvilinear coordinates, J. Comput. Phys. 267 (2014) 28–49. [27] K. Hejranfar, M. Hajihassanpour, A high-order nodal discontinuous Galerkin method for solution of compressible non-cavitating and cavitating flows, Comput. & Fluids 156 (2017) 175–199. [28] K. Hejranfar, A. Khajeh Saeed, Implementing a high-order accurate implicit operator scheme for solving steady incompressible viscous flows using artificial compressibility method, Internat. J. Numer. Methods Fluids (2010). [29] K. Hejranfar, K. Parseh, Application of a preconditioned high-order accurate artificial compressibility-based incompressible flow solver in wide range of Reynolds numbers, Internat. J. Numer. Methods Fluids (2017). [30] K. Hejranfar, K. Parseh, Preconditioned characteristic boundary conditions based on artificial compressibility method for solution of incompressible flows, J. Comput. Phys. 345 (2017) 543–564. [31] J.S. Hesthaven, T. Warburton, Nodal high-order methods on unstructured grids, J. Comput. Phys. 181 (1) (2002) 186–221. [32] J.S. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods, Texts Appl. Math. (2008). [33] A. Kanevsky, M.H. Carpenter, D. Gottlieb, J.S. Hesthaven, Application of implicit–explicit high order Runge–Kutta methods to discontinuous-Galerkin schemes, J. Comput. Phys. 225 (2) (2007) 1753–1781. [34] A. Klöckner, T. Warburton, J. Bridge, J.S. Hesthaven, Nodal discontinuous Galerkin methods on graphics processors, J. Comput. Phys. 228 (21) (2009) 7863–7882. [35] D. Kwak, C. Kiris, J. Dacles-Mariani, An assessment of artificial compressibility and pressure projection methods for incompressible flow simulations, Lecture Notes in Physics, 177–182. [36] E. Lee, H.T. Ahn, A reconstruction-based cell-centered high-order finite volume method for incompressible viscous flow simulation on unstructured meshes, Comput. & Fluids 170 (2018) 187–196. [37] C. Lehrenfeld, J. Schöberl, High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows, Comput. Methods Appl. Mech. Engrg. 307 (2016) 339–361. [38] S.-J. Li, L.-S. Luo, Z.J. Wang, L. Ju, An exponential time-integrator scheme for steady and unsteady inviscid flows, J. Comput. Phys. 365 (2018) 206–225.

214

M. Hajihassanpour and K. Hejranfar / Mathematics and Computers in Simulation 168 (2020) 173–214

[39] K. Liu, D. Ma, D. Sun, X. Yin, Wake patterns of flow past a pair of circular cylinders in side-by-side arrangements at low Reynolds numbers, J. Hydrodyn. Ser. B 19 (6) (2007) 690–697. [40] J.-G. Liu, C.-W. Shu, A high-order discontinuous Galerkin method for 2D incompressible flows, J. Comput. Phys. 160 (2) (2000) 577–596. [41] C. Liu, X. Zheng, C. Liao, C. Sung, T. Huang, C. Liu, X. Zheng, C. Liao, C. Sung, T. Huang, Preconditioned multigrid methods for unsteady incompressible flows, in: 35th Aerospace Sciences Meeting and Exhibit, 1997. [42] A.G. Malan, R.W. Lewis, P. Nithiarasu, An improved unsteady, unstructured, artificial compressibility, finite volume scheme for viscous incompressible flows: Part I. Theory and implementation, Internat. J. Numer. Methods Engrg. 54 (5) (2002) 695–714. [43] E.D. Mercerat, N. Glinsky, A nodal high-order discontinuous Galerkin method for elastic wave propagation in arbitrary heterogeneous media, Geophys. J. Int. 201 (2) (2015) 1101–1118. [44] A. Modave, A. St-Cyr, W.A. Mulder, T. Warburton, A nodal discontinuous Galerkin method for reverse-time migration on GPU clusters, Geophys. J. Int. 203 (2) (2015) 1419–1435. [45] N.C. Nguyen, J. Peraire, B. Cockburn, An implicit high-order hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations, J. Comput. Phys. 230 (4) (2011) 1147–1170. [46] F. Nieuwstadt, H.B. Keller, Viscous flow past circular cylinders, Comput. & Fluids 1 (1) (1973) 59–71. [47] P. Nithiarasu, An efficient artificial compressibility (AC) scheme based on the characteristic based split (CBS) method for incompressible flows, Internat. J. Numer. Methods Engrg. 56 (13) (2003) 1815–1845. [48] M. Paipuri, S. Fernández-Méndez, C. Tiago, Comparison of high-order continuous and hybridizable discontinuous Galerkin methods for incompressible fluid flow problems, Math. Comput. Simulation 153 (2018) 35–58. [49] M. Paipuri, C. Tiago, S. Fernández-Méndez, Coupling of continuous and hybridizable discontinuous Galerkin methods: Application to conjugate heat transfer problem, J. Sci. Comput. 78 (1) (2018) 321–350. [50] K. Parseh, K. Hejranfar, Assessment of characteristic boundary conditions based on the artificial compressibility method in generalized curvilinear coordinates for solution of the Euler equations, Comput. Methods Appl. Math. 18 (4) (2018) 717–740. [51] A. Pentaris, K. Nikolados, S. Tsangaris, Development of projection and artificial compressibility methodologies using the approximate factorization technique, Internat. J. Numer. Methods Fluids 19 (11) (1994) 1013–1038. [52] J. Qiu, B.C. Khoo, C.-W. Shu, A numerical study for the performance of the Runge–Kutta discontinuous Galerkin method based on different numerical fluxes, J. Comput. Phys. 212 (2) (2006) 540–565. [53] S.M. Renda, R. Hartmann, C. De Bartolo, M. Wallraff, A high-order discontinuous Galerkin method for all-speed flows, Internat. J. Numer. Methods Fluids 77 (4) (2014) 224–247. [54] F. Rouzbahani, K. Hejranfar, A truly incompressible smoothed particle hydrodynamics based on artificial compressibility method, Comput. Phys. Comm. 210 (2017) 10–28. [55] D. Russell, Z. Jane Wang, A cartesian grid method for modeling multiple moving objects in 2D incompressible viscous flow, J. Comput. Phys. 191 (1) (2003) 177–205. [56] H. Sarvghad-Moghaddam, N. Nooredin, B. Ghadiri-Dehkordi, Numerical simulation of flow over two side-by-side circular cylinders, J. Hydrodyn. Ser. B 23 (6) (2011) 792–805. [57] W. Shao, J. Li, Three time integration methods for incompressible flows with discontinuous Galerkin Boltzmann method, Comput. Math. Appl. 75 (11) (2018) 4091–4106. [58] W. Shaoa, J. Li, Three time integration methods for incompressible flows with discontinuous Galerkin Boltzmann method, Comput. Math. Appl. 75 (2018) 4091–4106. [59] M. Tavelli, M. Dumbser, A staggered semi-implicit discontinuous Galerkin method for the two dimensional incompressible Navier–Stokes equations, Appl. Math. Comput. 248 (2014) 70–92. [60] D.J. Tritton, Experiments on the flow past a circular cylinder at low Reynolds numbers, J. Fluid Mech. 6 (04) (1959) 547. [61] A. Vakil, S.I. Green, Two-dimensional side-by-side circular cylinders at moderate Reynolds numbers, Comput. & Fluids 51 (1) (2011) 136–144. [62] D. Wirasaet, S. Tanaka, E.J. Kubatko, J.J. Westerink, C. Dawson, A performance comparison of nodal discontinuous Galerkin methods on triangles and quadrilaterals, Internat. J. Numer. Methods Fluids 64 (1336) (2010) 10–12–1362. [63] A. Zadehgol, M. Ashrafizaadeh, S.H. Musavi, A nodal discontinuous Galerkin lattice Boltzmann method for fluid flow problems, Comput. & Fluids 105 (2014) 58–65. [64] F. Zhang, J. Cheng, T. Liu, A direct discontinuous Galerkin method for the incompressible Navier–Stokes equations on arbitrary grids, J. Comput. Phys. 380 (2019) 269–294. [65] F. Zhang, J. Cheng, T. Liu, A high-order discontinuous Galerkin method for the incompressible Navier–Stokes equations on arbitrary grids, Internat. J. Numer. Methods Fluids 90 (5) (2019) 217–246.