A reconstructed discontinuous Galerkin method for ...

Computers and Fluids 160 (2018) 26–41

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

A reconstructed discontinuous Galerkin method for compressible turbulent flows on 3D curved grids Xiaodong Liu a, Yidong Xia b, Hong Luo a,∗ a b

Department of Mechanical and Aerospace Engineering, North Carolina State University, Raleigh, NC 27695, USA Energy and Environment Science & Technology Directorate, Idaho National Laboratory, Idaho Falls, ID 83415, USA

a r t i c l e

i n f o

Article history: Received 8 March 2017 Revised 5 October 2017 Accepted 11 October 2017 Available online 20 October 2017 Keywords: Discontinuous Galerkin WENO RANS Curved elements

a b s t r a c t A third-order accurate reconstructed discontinuous Galerkin method, namely rDG(P1 P2 ), is presented to solve the Reynolds-Averaged Navier–Stokes (RANS) equations, along with the modified one-equation model of Spalart and Allmaras (SA) on 3D curved grids. In this method, a piecewise quadratic polynomial solution (P2 ) is obtained using a least-squares method from the underlying piecewise linear DG(P1 ) solution. The reconstructed quadratic polynomial solution is then used for computing the inviscid and the viscous fluxes. Furthermore, Hermite Weighted Essentially Non-Oscillatory (WENO) reconstruction is used to guarantee the stability of the developed rDG method. A number of benchmark test cases based on a set of uniformly refined quadratic curved meshes are presented to assess the performance of the resultant rDG(P1 P2 ) method for turbulent flow problems. The numerical results demonstrate that the rDG(P1 P2 ) method is able to obtain reliable and accurate solutions to 3D compressible turbulent flows at a cost slightly higher than its underlying second-order DG(P1 ) method. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction The discontinuous Galerkin (DG) methods, originally introduced for solving the neutron transport by Reed and Hill [1], have become popular for the solution of systems of conservation laws in recent decades. Nowadays, they are widely used in computational fluid dynamics (CFD), computational acoustics, and computational magneto-hydrodynamics (MHD) [2]. A lot of attractive features of DG methods have been listed in [3–7,11,22]. However, the DG methods also have a number of weaknesses that have yet to be addressed, e.g., how to reduce the high computational costs and how to develop more efficient time integration methods. In order to reduce the high costs associated to the DG methods, Dumbser et al. [8–10] introduced a new family of so-called reconstructed DG, termed Pn Pm schemes and referred to as rDG (Pn Pm ) in this paper. Pn indicates that a piecewise polynomial of degree of n is used to represent the underlying DG solution, and Pm represents a polynomial solution of degree of m (m ≥ n) that is reconstructed from the underlying Pn polynomial and used to compute the fluxes. The Pn Pm schemes can be constructed based on a few different algorithms, e.g., the recovery approach [12], the reconstruction approach [14,25], and the Gauss–Green approach [16,17],

∗

Corresponding author. E-mail address: [email protected] (H. Luo).

https://doi.org/10.1016/j.compfluid.2017.10.014 0045-7930/© 2017 Elsevier Ltd. All rights reserved.

all of which were proved to deliver the designed grid convergence of O (hm+1 ) [18]. Indeed, implicit methods can especially benefit from the use of rDG(Pn Pm ) methods as the costs can be substantially reduced in two aspects [19,20]. Firstly, fewer spatial integration points are required for evaluating the residual vector and Jacobian matrix. For instance, the third-order rDG (P1 P2 ) only needs 4 points for triangular boundary integral whereas the equivalent DG (P2 ) requires 7. Secondly, the Jacobian matrix of rDG (Pn Pm ) is based on the underlying DG (Pn ), and thus requires much less storage than the equivalent DG (Pm ). For example, for RANS-SA system, the memory needed for the diagonal part of the Jacobian matrix of rDG(P1 P2 ) is 576 word versus 3600 needed by DG (P2 ) for 3D cases. In our latest work, a rDG(P1 P2 ) method based on a Hierarchical WENO reconstruction has been successfully used to solve the compressible flows on unstructured grids, ranging from inviscid [19– 22,43] to viscous flows [13,20,28,43,45]. This rDG(P1 P2 ) method is designed not only to reduce the high computing costs of the DGM, but also to avoid spurious oscillations in the vicinity of strong discontinuities, thus effectively addressing the two shortcomings of the DGM. However, most practical applications belong to 3D turbulent flows. Therefore, the great success for rDG(P1 P2 ) has motivated us to extend this promising method to solution of the practical 3D turbulent flows. In recent years, development and application of high-order methods for computing high Reynolds number turbulent flows

X. Liu et al. / Computers and Fluids 160 (2018) 26–41

governed by the Reynolds-averaged Navier–Stokes (RANS) equations have been an active research topic all over the world, as demonstrated by the two European projects: ADIGMA and IDIHOM [49,50]. It is well known that the use of high-order methods for computing RANS problems is not commonplace, mainly due to the severe numerical instability by the highly non-smooth behavior between the turbulent and non-turbulent flow regions. Although it is still a challenging problem, there are several successful implementations for the one-equation Spalart and Allmaras (SA) model and the two-equation k − ω model. In [33] by Bassi and Rebay, they successfully solved the RANS equations with a modified k − ω model based on a high-order DG framework. In their work, the logarithm of ω rather than ω itself is used as the unknown, which has been found very useful to enhance stabilty of the method. Similarly, Hartmann et al. [34] developed a DG code for 3D turbulent flow computation with adaptive mesh refinement based on the k − ω model. In particular, some test cases with increasing complexity have been used to validate the solver. Burgess et al. [35] and Wang et al. [29] developed high-order DG methods for solving a fully coupled RANS-SA system respectively. The modified SA model is particularly designed to make the original SA model insensitive to negative values of turbulence working variables as numerical experiments show that turbulence working variable often drops several orders of magnitudes at the edge of the turbulent boundary layer. Ceze and Fidkowski [36] applied a high-order output-based adaptive solution technique to the 2D RANS equations closed with the modified negative SA model which originally proposed by Allmaras et al. [37]. Compared with uniform refinement at second order, high-order yields faster convergence. Nguyen et al. [39] implemented the SA model equation based on a DG framework with an artificial viscosity modification for SA equation. It is aimed to accomodate high-order RANS approximations on too coarse grids. Besides that, Crivellini et al. [38] analyzed and implemented a modified SA model based on high-order DG methods for incompressible flows. For modified SA model, they introduce an SA model implementation that dealt with negative ν˜ values by modifying the source and diffusion terms in the SA model equation only when the working variable or one of the model closure functions became negative. Zhou et al. [40] successfully solved RANS-SA system in a high-order correction procedure via reconstruction (CPR) framework. In their work, a high-order solver based on CPR has been developed for the Eikonal equation to compute the nearest distance to the wall. In addition, the use of high-order curved boundary elements is essential for high-order schemes to deliver an overall accurate solution [29]. Poor representation of the actual geometry could result in a significant amount of artificial entropy produced along the geometry surface, thus degrading the solution accuracy. Furthermore, since the elements are of high-aspect-ratio through the thin boundary layer for high Reynolds number flow, the interior elements are also required to be curvilinear to avoid the negative cells. The objective of this paper is to develop a reconstructed discontinuous Galerkin method for the solution of RANS-SA system on 3D curved grids. First, the quadratic curved elements generated by Gmsh [30] or simple agglomeration are shown satisfatory for the solution of high order reconstructed discontinuous Galerkin method. Second, grid convergence study on uniformly refined meshes with high-order reconstructed discontinuous Galerkin is reported about the 3D turbulent benchmark test cases with increasing complexity from NASA official website [47]. This provides direct and valuable comparison with second-order finite volume method. Since a WENO reconstruction has to be used for unstructured meshes, e.g.prism and tetrahedron, to guarantee the stability of rDG(P1 P2 ), it is necessary to demonstrate that WENO(P1 P2 ) could still guarantee the accuracy of the designed

27

rDG(P1 P2 ) method. Thus, third, the accuracy of WENO(P1 P2 ) is also validated by comparison with the output of rDG(P1 P2 ) without WENO reconstruction on subsonic turbulent flow cases, i.e. Zero Pressure Gradient Flat Plate. The outline of the rest of this paper is organized as follows. The governing equations are described in Section 2. The developed reconstructed discontinuous Galerkin method is presented in Section 3. The use of curved elements is introduced in Section 4. Extensive numerical experiments are reported in Section 5. Concluding remarks are given in Section 6.

2. Governing equations The conservation form of the compressible RANS equations with the modified one-equation Spalart–Allmaras (SA) turbulence model [44] is given as below,

∂ U ∂ F j (U ) ∂ G j (U ) + = +S ∂t ∂xj ∂xj

(1)

where the summation convention (j = 1, 2, 3) has been used. In Eq. (1), the conservative variables U are defined as

U = (ρ , ρ u j , ρ e, ρ ν˜ )T

(2)

where ρ , p and e denote the density, pressure and specific total energy of the fluid, respectively, uj is the velocity components of the flow in the coordinate direction xj , and ν˜ represents the turbulence working variable in the modified SA model. The invisid and viscous flux vector, i.e. F and G, are defined by

⎛

⎞

ρu j

⎛

⎜ρ ui u j + pδi j ⎟ , u j (ρ e + p ) ⎠ ρ u j ν˜

⎜

Fj = ⎝

Gk = ⎝

⎞

0

τi j

⎟

u i τi j + q j ⎠ 1 ∂ ν˜ μ σ (1 + ψ ) ∂ x

(3)

j

and the source term S is defined by

ST = 0

0

0

0

0

s

(4)

s = cb1 S˜μψ + σb2 ρ∇ ν˜ · ∇ ν˜ − cw1 ρ fw ( νψ )2 − σ1 ν (1 + d c

where

ψ )∇ ρ · ∇ ν˜ .

The Newtonian fluid with the Stokes hypothesis is valid under the current framework, since only air is considered. The viscous stress tensor τ is defined by

∂ ui ∂ u j 2 ∂ uk τi j = (μ + μt ) + − δ ∂ x j ∂ xi 3 ∂ xk i j

(5)

where δ ij is the Kronecker delta function, μ represents the molecular viscosity coefficient, which can be determined through Sutherland’s law, and μt denotes the turbulence eddy viscosity, which is given by:

μt =

ρ ν˜ fv1 i f ν˜ ≥ 0 0 i f ν˜ < 0

(6)

The heat flux vector qj , which is formulated according to Fourier’s law, is given by

q j = −c p

μ Pr

+

μt ∂ T P rt ∂ x j

(7)

where cp is the specific heat capacity at constant pressure , Pr is the nondimensional laminar Prandtl number, Prt is the turbulent Prandtl number, and the temperature of the fluid T is determined by T = ρpR .

28


Finally, the source term for the SA equation is introduced in detail. The parameters for the production and destruction terms of the modified SA turbulence model are given as,

S + Sˆ

S˜ =

S+

(cv3 −2cv2 )S−Sˆ

S=

where f v2 = 1 − f w = g(

Sˆ ≥ −cv2 S ˆ S < −cv2 S

if S(cv22 +cv3 Sˆ)

ψ

1+ψ fv1

6 1+cw 3

if

− → − ω ·→ ω,

− →

ω = ∇ × U,

(8)

SκT d

Sˆ = νψ f , κT2 d2 v2

f v1 =

ψ3

ψ 3 +cv31

,

)6 .

ψ=

0.05ln(1 + e20χ )

χ

χ ≤ 10 χ > 10

if if

(9)

where χ = νν˜ . The parameter ψ is designed to become zero as the turbulence working variable goes negative, thereby turning off the production, destruction, and dissipation terms to guarantee stabilities. The constants in the modified SA model to close the main flow equations are given by, c cb1 = 0.1355, σ = 23 , cb2 = 0.622, κT = 0.41, cw1 = κb12 + 1+cb2

σ

c w3 = 2,

c v1 = 7.1,

c v2 = 0.7,

c v3 = 0.9.

3.1. Reconstructed DG formulation for RANS-SA system The governing equation Eq. (1) is discretized using a discontinuous Galerkin finite element formulation. To formulate the discontinuous Galerkin method, we first introduce the following weak formulation, which is obtained by multiplying the above conservation law by a test function W, integrating over the domain , and then performing an integration by parts,

∂U W d + ∂t

=

G k nk d −

F k nk d −

Gk

Fk

∂W d + ∂ xk

∂W d ∂ xk

SW d

(10)

where (= ∂ ) denotes the boundary of , and n j the unit outward normal vector to a boundary. We assume that the domain is subdivided into a collection of non-overlapping elements e . We p introduce the following broken Sobolev space Vh

Vhp = {vh ∈ [L2 ( )]m : vh | e ∈ [Vpm ]

e

∀ e ∈ } ,

(11)

D = span xαi i : 0 ≤ α ≤ p, 0 ≤ i ≤ D ,

F k (U h )

Gk (U h , ∇ U h )nk W h d − S ( U h , ∇ U h )W h d

e

∂Wh d ∂ xk

Gk (U h , ∇ U h )

∂Wh d ∂ xk

, ∀W h ∈ Vhp

(13)

where U h and W h represent the finite element approximations to the analytical solution U and the test function W , respectively, and they are approximated by a piecewise polynomial function of degrees p, which are discontinuous at the cell interfaces. Assume that B is the basis of polynomial function of degrees p, this is then equivalent to the following system of N equations,

d U Bi d dt e h ∂ Bi v L R =− H in ( U , U , n ) B d

+ F k (U h ) d i k h h k ∂ xk

e e + H vkis (U Lh , ∇ U Lh + β rhL , U Rh , ∇ U Rh + β rhR , nk )Bi d −

e

e

+

e

Gk (U h , ∇ U h + R )

∂ Bi d ∂ xk

S ( U h , ∇ U h + R )Bi d

,1 ≤ i ≤ N

(14)

e

rh B i d = R=

LR

1 R (U − U Lh )nk Bi d , 2 h

1≤i≤N

rh

(15)

LR

where LR represents the common face shared by cell R and L. It is necessary to mention that, through our experiment, there exists lower bounds of the parameter β , e.g. 2 or 3, to ensure the stability of the method. This scheme is called discontinuous Galerkin method of degree p, or in short notation DG(p) method. In the DG framework, numerical polynomial solutions Uh in each element are expressed using either standard Lagrange finite element or hierarchical nodebased basis as below

Uh =

N

Ui (t )Bi (x ),

(16)

i

which consists of discontinuous vector-values polynomial functions of degree p, and where m is the dimension of the vector of conserved variables and p Vm

+

e

where N is the dimension of the polynomial space. The inviscid and viscous flux functions appearing in Eq. (13) are replaced by v numerical flux function Hin and Hvkis , respectively, where ULh and k R Uh are the solution polynomial at the left and right states of one cell interface. In the viscous flux, rh is the local lifting operator, β is called penalty parameter [46], and R is the global lifting operator [33]. The local and global lift operators are defined as following,

3. Discretizations

e

T

, c w2 = 0. 3,

e

g = r + cw2 (r 6 − r ), and

The parameter d denotes the distance from a specific position, i.e. a Gauss quadrature point, in the flow field to the nearest wall, which will be discussed in Section 4. The parameter ψ is designed for high-order discretization schemes to remove the effects of negative turbulence working variable to guarantee the robustness of the turbulence model. This parameter ψ is given by,

d U W d dt e h h =− F k (U h )nk W h d + +

1

6 g6 +cw 3

r = min(10.0, ˜ νψ 2 2 ),

,

p

Find U h ∈ Vh , such as

(12)

i=1

where α denotes a multi-index and D is the dimension of space. Then, we can obtain the following semi-discrete form by applying weak formulation on each element e .

where Bi are the finite element basis functions. As a result, the unknowns to be solved are the variables at the nodes. In the previous work [7,19–22,43], the piecewise polynomial solutions are represented using a linear Taylor series expansion at the cell centroid, which can be expressed as a combination of cell-averaged variables and their gradients at the cell centers regardless of the element shapes, shown in Eq. (17).

˜ + Ux,c B˜2 + Uy,c B˜3 + Uz,c B˜4 + Uxx,c B˜5 + Uyy,c B˜6 Uh = U +Uzz,c B˜7 + Uxy,c B˜8 + Uxz,c B˜9 + Uyz,c B˜10

(17)

The unknowns to be solved in this formulation are the cell˜ and their derivatives at the center of the cells, averaged variables U


B˜ j B˜ j B˜ j ∂U 1 = URz,i + URzz,i 4 + URxz,i 2 + URyz,i 3 ∂z j z i z i z i z i

shown in Eq. (18).

∂U x ∂x c ∂ 2U Uxx,c = x2 2 ∂x c ∂ 2U Uxy,c = xy ∂ x∂ y c Ux,c =

∂U y ∂y c ∂ 2U Uyy,c = y2 2 ∂y c ∂ 2U Uxz,c = xz ∂ x∂ z c Uy,c =

∂U z ∂z c ∂ 2U 2 Uzz,c = z 2 ∂z c 2 ∂ U Uyz,c = yz ∂ y∂ z c Uz,c =

where x = 0.5(xmax − xmin ), y = 0.5(ymax − ymin ), z = 0.5(zmax − zmin ) and xmax , ymax , zmax and xmin , ymin , zmin are the maximum and minimum coordinates in the cell e in x, y and z directions, respectively. The basis functions are defined in Eq. (19) as follows,

B˜2 1 B˜5 = 2 − 2 e B˜7 =

x − xc , B˜3 x ˜2 B2 d , e 2 ˜2 B4 d , e 2 B˜2 B˜4 d ,

B˜2 =

B˜24 1 − 2 e

B˜9 = B˜2 B˜4 −

1

e

e

y − yc z − zc , B˜4 = y z ˜2 B˜23 B3 1 B˜6 = − d 2 e e 2 1 B˜8 = B˜2 B˜3 − B˜2 B˜3 d

=

e

B˜10 = B˜3 B˜4 −

e

e

B˜3 B˜4 d (19)

=

˜R U i

+

URx,i B˜2

+URzz,i B˜7

+

+

URy,i B˜3

URxy,i B˜8

+

+

URz,i B˜4

URxz,i B˜9

+

+

URxx,i B˜5

+

URyy,i B˜6

URyz,i B˜10

(20)

In order to maintain the compactness of the DG methods, the reconstruction is required to involve only von Neumann neighborhood, i.e. the adjacent cells that share a face with the cell i under consideration. There are 10 degrees of freedom, and therefore 10 unknowns must be determined. The first 4 unknowns can be trivially obtained, by requiring the consistency of the rDG method with the underlying DG(P1 ) method: (1) The reconstruction scheme must be conservative, and (2) The values of the reconstructed first derivatives are equal to the ones of the first derivatives of the underlying DG solution at the centroid i. Due to the judicious choice of Taylor basis in our DG formulation, these 4 degrees of freedom simply coincide with the ones from the underlying DG solution, namely

˜R = U ˜ i, U i

URx,i = Ux,i ,

URy,i = Uy,i ,

URz,i = Uz,i

(21)

As a result, only 6 second derivatives need to be determined. This can be accomplished as follows. For any neighboring cell j, one requires

˜j = U

1

j

j

˜ R + UR B˜ j + UR B˜ j + UR B˜ j U i x,i 2 y,i 3 z,i 4

+URxx,i B˜5j + URyy,i B˜6j + URzz,i B˜7j + URxy,i B˜8j

j +URxz,i B˜9j + URyz,i B˜10 d j

B˜ j B˜ j B˜ j ∂U 1 = URx,i + URxx,i 2 + URxy,i 3 + URxz,i 4 ∂x j x i x i x i x i j j B˜ B˜ B˜ j ∂U 1 = URy,i + URyy,i 3 + URxy,i 2 + URyz,i 4 ∂y j y i y i y i y i

URxx,i

⎞

⎜ R ⎟ ⎜Uyy,i ⎟ ⎜ ⎟ ⎜ R ⎟ ⎜Uzz,i ⎟ ⎜ ⎟ A × ⎜ ⎟ ⎜URxy,i ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ UR ⎟ xz,i ⎝ ⎠ URyz,i

⎞ R R R ˜R ˜j ˜j ˜j ⎜U j − Ui d j + Ux,i j B2 d j + Uy,i j B3 d j + Uz,i j B4 d j ⎟ ⎜ ⎟ ⎜ ⎟ x i ˜ ˜ x,i ⎜ ⎟ Ux, j − U ⎜ ⎟ x j = ⎜ ⎟ (23) y ⎜ ⎟ i ˜ ˜ Uy, j − Uy,i ⎜ ⎟ y j ⎜ ⎟ ⎝ ⎠ z i

z j

As a result, the very same numerical polynomial solutions are used for arbitrary shapes of elements, which can be triangle, quadrilateral, and polygon in 2D, and tetrahedron, pyramid, prism, and hexahedron in 3D [23]. This makes the implementation of the so-called rDG (P1 P2 ) reconstruction straightforward [15,19– 22,24,25,43]. For this scheme, using the underlying linear polynomial DG(P1 ) solution in the neighboring cells, one can reconstruct a quadratic polynomial solution URi as follows:

URi

⎛

⎛

e

1

(22)

where the basis functions B are evaluated at the center of cell j, ˜ (x j , y j , z j ). The explicit expressions about B ˜ (x, y, z ) namely B˜ j = B could be found by Eq. (19). Finally, this group of equations can be written in a matrix form as follows:

(18)

B˜1 = 1,

29

˜ z,i ˜ z, j − U U

where

⎛ ⎜ ⎜ ⎜ A=⎜ ⎜ ⎝

j

B˜5j d j

j

B˜6j d j

j

B˜7j d j

j

B˜8j d j

j

B˜9j d j

j

j B˜10 d j

B˜2j

0

0

B˜3j

B˜4j

0

0

B˜3j

0

B˜2j

0

B˜4j

0

0

B˜4j

0

B˜2j

B˜3j

⎞ ⎟ ⎟ ⎟ ⎟. ⎟ ⎠

Similar equations can be written for all the cells connected to the cell i with a common face, which leads to a non-square matrix. The numbers of the face-neighboring cells for a tetrahedron, a pyramid, a prism and a hexahedron are four, five, five and six, respectively. Consequently, the size of the resulting non-square matrix is 16 × 6, 20 × 6, 20 × 6 and 24 × 6, respectively. In the present work, this over-determined linear system of 16, or 20, or 24 equations for 6 unknowns is solved in the least-squares sense using a normal equation approach to obtain the second derivatives of the reconstructed quadratic polynomial solution. The associated numerical details can be found in [19–22,28,43]. The final P2 solution is then obtained through a WENO reconstruction on the curvatures in each cell, and used for evaluating the inviscid flux. Note that as previous investigation indicates in [21], although the least-squares reconstructed rDG(P1 P2 ) method has been successfully used to solve the 2D compressible Euler equations on arbitrary grids [24], yet when used to solve the 3D compressible Euler equations on tetrahedral grids, it suffers from the so-called linear instability that is also observed in the second-order cell-centered finite volume methods [32]. Such linear instability can occur even when the 3D linear equations are being solved for smooth problems on tetrahedral grids. Indeed, this linear instability is attributed to the fact that the reconstruction stencils only involve von Neumann neighborhood [26]. The linear stability can be achieved using extended stencils, which will unfortunately sacrifice the compactness of the underlying DG methods. Alternatively, the non-linear ENO and WENO approaches can be used to achieve both the linear and non-linear stability. It has been demonstrated that rDG(P1 P2 ) method based on WENO reconstruction, termed as WENO(P1 P2 ) scheme, is able to deliver the designed third-order spatial accuracy for steady-state flow problems without much extra cost in computing time and storage than the

30


underlying second-order DG (P1 ) method, and also eliminate linear instability on unstructured tetrahedral grids [21]. In addition, in order to remove inherent oscillations in the underlying DG(P1 ) solutions, reconstructed discontinuous Galerkin method based on Hierarchical WENO reconstruction, termed as HWENO(P1 P2 ), is adopted [22]. In this paper, only turbulent benchmark test cases involving smooth flows are considered, therefore, we will only focus on rDG(P1 P2 ) and WENO(P1 P2 ). The verification of HWENO(P1 P2 ) for capturing the shock in the turbulent flows will be reported in a following paper. Regarding the discretization of the turbulence model equation, Burgess et al. [35] proposed 2 options for discretizing the turbulence model equation. In this paper, unlike the mean flow equation discretized with reconstructed discontinuous galerkin method, the modified SA model is just discretized with the underlying second-order DG(P1 ) method.

The spatial discretization of the governing equations leads to a system of ordinary differential equations (ODEs) in time and Eq. (14) can be written in an elemental semi-discrete form as

dU = R (U ) dt

(24)

For steady state problems, the spatially discretized governing equations Eq. (24) is integrated implicitly in time using BDF1, then we can arrive at

M

U n+1 − U n

t

= R U n+1

(25)

Since the governing equations are nonlinear, Eq. (25) is a system of nonlinear equations for the global solution vector U n+1 . In order to solve this type of equations, we can linearize the right-handside vector R with respect to the global solution vector U at the current time level n

R U n+1 ≈ R U n +

∂R ∂U

n

U n+1 − U n

(26)

n

where ∂∂ UR is the so-called Jacobian matrix of the system evaluated at time level n, and denoted symbolically as J(Un ), which involves the linearization of both the inviscid and viscous flux functions. If the right-hand-side term in Eq. (25) is replaced by Eq. (26), n n+1 and move ∂∂ UR U − U n to the left side, it then leads to a delta form of the equations

M − t

∂R ∂U

n

U=

∂ R ( U i , U j , ni j ) ∂U j

∂ R ( U i , U j , ni j ) ∂ Ui ∂ R(U i , U j , ni j ) M D= + t ∂ Ui j L=−

(28)

where, U, L and D represent the strict upper matrix, the strict lower matrix, and the diagonal matrix, respectively. Since reconstruction is not taken into consideration for Jacobian matrix, this makes implicit rDG(P1 P2 ) method suitable for parallization. Finally, BDF1 is used to drive the above equation to converge to steady state using a GMRES solver with LU-SGS preconditioner. The associated details about the LU-SGS precondioned GMRES for the solution of Eq. (27) could refer to [19,20,41]. 4. High-order curved meshes

3.2. Implicit temporal discretization

M

expressed as,

Un = R Un

(27)

where t is the time increment and U n = U n+1 − U n is the difference of global solution vector between time level n and n + 1. Note that if t tends to infinity, the scheme reduces to the standard Newton’s method with a property of quadratic convergence for solving a system of nonlinear equations. In this work, in order to guarantee the robustness of implicit solver, for the beginning time steps, small CFL No., e.g. 10, is adopted; then fixed bigger CFL, e.g. 10 0 0, is used to accelerate the convergence. In addition, in this work, acceleration techniques such as local time-stepping has been adopted. For rDG(P1 P2 ) method, the Jacobian matrix used in implicit temporal discretization is obtained by linearizing the RHS of underlying DG(P1 ) with respect to the solution vector U. In our work, The Jacobian is calculated through an auto differential tool [27]. And by using an edge-based data structure, the left hand side matrix is stored in upper, lower, and diagonal forms, which can be

4.1. Curved meshes For second-order accurate discretization used to solve the RANS-SA system, the geometry surface is usually represented by linear elements. However, it is known that the use of high-order curved boundary elements is essential for high-order schemes [29]. Therefore, in this paper, two methods were used to generate the quadratic curved meshes for high-order rDG method. The first is to use Gmsh [30] library to generate the quadratic curved meshes when possible. The second is based on the simple agglomeration of the linear meshes, which is very easy to implement, especially for structured hexahedral meshes. But for this method [31], it requires that the block size in each of the ijk-directions is a multiple of p then the mesh can be agglomerated by keeping every pth point in each direction, using these remaining points for defining polynomial functions of degrees p without any change of the position. The numerical algorithms previously described in Section 3.2 do not require any special modification for the use of curved meshes. Also the reconstruction procedure is not affected by the curved meshes. But the domain and boundary integrals are indeed affected by the use of the curved meshes. In the current implementation, the domain and boundary integrals are evaluated using Gauss quadrature formulas. The number of quadrature points used is chosen to integrate exactly polynomials of order of 2p on the reference element. Herein, let us take 2D as an example to explain the aforementioned. In the case of linear, quadratic, and cubic shape function, the domain integrals are evaluated using three, six, and twelve points for triangles and using four, nine, and sixteen points for quadrilaterals, respectively, and the boundary integrals are evaluated using two, three, and four points, respectively, for 2D. As described in Eq. (8), the parameter d denotes the distance from a Gaussian quadrature of the computational domain to the nearest wall surface [42]. Note that, due to the curved surface, special care has to be taken to calculate the distance. Assuming there exists one Gaussian quadrature point, P(x, y, z), and the nearest position, P0 (ξ 0 , η0 ), is located in the curved wall surface, then vector − → P P0 is vertical to the tangential plane through P0 (ξ 0 , η0 ) in the curved wall surface. It is convenient to give the tangential plane through parametric expression of the curved wall surface based on the geometric points. This leads to the non-linear sets of equations. Therefore, Newton method could be used to solve the non-linear equations, whereas, it requires a good initial estimate. In our work, some different points (ξ , η) in the reference space are used as the initial estimates to drive the Newton method rather than only one initial good estimate. Then with different initial estimates, maybe different roots could be obtained, which is determined by the nature of Newton method. As a result, the minimum distance among

X. Liu et al. / Computers and Fluids 160 (2018) 26–41 Table 1 Grid information for turbulent flow over a flat plate at Ma∞ = 0.2, and Re∞ = 5 × 106 . Grid no.

Dimensions

Off-wall spacing

Leading-edge spacing

y+

Level Level Level Level

2 × 35 × 25 2 × 69 × 49 2 × 137 × 97 2 × 273 × 193

8.32e−6 4.04e−6 2.00e−6 1.00e−6

1.62e−2 8.05e−3 4.01e−3 2.00e−3

1.70 0.83 0.41 0.21

1 2 3 4

all the distance from different initial guess positions could be selected as the distance from the integration point to the curved surface. 5. Numerical examples A number of numerical test cases with varying complexities have been presented in this section to validate and assess the performance of the modified SA turbulence model implemented in the rDG code. The test cases considered in this study are chosen from the NASA Langley Turbulence Modeling Resource Database [47], which provides multiple 2D and 3D test problems used for the validation and verification of turbulence models, as well as the reference results by several state-of-the-art CFD codes (e.g. FUN3D, and CFL3D). The grids provided by the database are suitable for use in many finite volume or finite difference solvers. Flow conditions and geometry are chosen to approximate those employed by Rumsey [47]. In this work, all the 2D test cases are run in 3D framework. A parallel implicit solver with the LU-SGS preconditioned GMRES algorithm is used for time marching, and the computation is continued until a steady-state solution is reached, where the residual of the mass equation drops 7 orders of magnitude if possible. The initial solution for third-order rDG(P1 P2 ) in each test case was taken from a fully converged solution of DG(P1 ). In addition the freestream turbulent viscosity has been set as ν˜ = 3ν∞ for all the computations in this section.

31

layer. Notice however that the eddy viscosity remains non-negative across the entire domain even if ν˜ becomes negative. It could be observed that, by comparison, the results using WENO(P1 P2 ) are almost identical to those using rDG(P1 P2 ). We need to mention that, the Hermite WENO reconstruction-based DG method WENO(P1 P2 ), is designed not only to enhance the accuracy of DG method but also to ensure linear stability of the rDG method on unstructured grids, e.g. prisms and tetrahedrons. Therefore, using this case with unstructured hexahedral mesh, we have furthure demontrated the WENO(P1 P2 ) could keep the accuracy of rDG method. Fig. 5 shows the grid convergence of the drag coefficient, i.e. CD . In the plot the x-axis is plotting the number of degrees of freedom, using logarithmic scale. From Fig. 5, we could observe some interesting phenomenon. First, DG(P1 ), rDG(P1 P2 ) and WENO(P1 P2 ) all could converge to the grid-independent solution. And it is obvious that rDG(P1 P2 ) is superior to DG(P1 ) in terms of number of degrees of freedom. This means that, third-order rDG(P1 P2 ) only requires less degrees of freedom to arrive at the grid-independent solution compared with DG(P1 ). Furthermore, the grid convergence of the drag coefficient is almost indistinguishable between rDG(P1 P2 ) and WENO(P1 P2 ), which demontrate the accuracy of WENO(P1 P2 ). Second, as a matter of fact, there exists a little difference between our converged solution and the reference solution of CFL3D and FUN3D. Obviously, even for well-validated CFL3D and FUN3D, there exists a little difference for the converged solution. This may be caused by different implementation of turbulence model equation. Compared with CFD3D, our converged solution agrees well with the reference solution up to 4 decimals, and the difference is less than 0.42%, which is acceptable. Finally, the typical convergence histories using different numerical methods on level 2 hexahedral grid are presented in Fig. 6. It could be observed that, the residual of DG(P1 ) and rDG(P1 P2 ) could drop 7 orders of magnitude, whereas the residual of WENO(P1 P2 ) keeps oscillating after it drops 4 orders of magnitude. This could be explained as follows. This phenomenon is sometimes called “convergence stall”, especially when nonlinear approaches like WENO schemes are used.

5.1. Zero pressure gradient flat plate

5.2. 2D turbulent flow over a 2D bump in a channel

The first test case to validate our implementation of the turbulence model, is a turbulent flow over a flat plate. For this turbulent flat plate case, Ma∞ = 0.2, and Re∞ = 5 × 106 based on unit length of the grid. For the flat plate, the linear mesh is accurate enough to represent the geometry. Therefore, the linear meshes from the NASA database are directly used for this case with highorder rDG method. The detailed information of such meshes, e.g., dimensions, is listed in Table 1. Fig. 1 shows the level 3 hexahedral grid used for computation. This grid consists of 136 cells in the x-direction, 1 cell in y-direction and 96 cells in z-direction, which has been specified in Table 1. Fig. 2 shows the contours of nondimensionalized turbulent eddy viscosity, i.e. μt /μ∞ , using rDG(P1 P2 ) on level 3 hexahedral grid, which is almost indistinguishable from the contours distribution from NASA database. In addition to DG(P1 ) and rDG(P1 P2 ), the numerical results using WENO(P1 P2 ) are also presented. It should be noted that, rDG(P1 P2 ) without WENO reconstruction on P2 works well for this case. The reason why the numerical results using WENO(P1 P2 ) is presented is that we just want to verify the accuracy of WENO(P1 P2 ) compared with rDG(P1 P2 ). Fig. 3 presents the logarithmic plot of friction coefficient distribution along the flat plate with DG(P1 ), rDG(P1 P2 ) and WENO(P1 P2 ) on level 3 hexahedral grid. Also, Fig. 4 presents the nondimensionalized turbulence eddy viscosity at x = 0.97 with DG(P1 ), rDG(P1 P2 ) and WENO(P1 P2 ) on level 3 hexahedral grid. As expected, the eddy viscosity drops to zero at the wall and approaches nearly zero rapidly at the edge of the boundary

This test case simulates a turbulent flow over a 2D bump in a channel, where Ma∞ = 0.2, and Re∞ = 3 × 106 . This test case is different from the flat plate test case above as it involves wall curvature, thus, pressure gradients is no longer equal to zero. Symmetric boundary conditions are used for the top and bottom of the channel, with exception of a smooth bump extending from x = 0.0 to x = 1.5, where adiabatic wall boundary condition is applied. The whole computational domain is [−25.0, 25.0] × [0.0, 5.0]. 3 hexahedral meshes used for this test case, are fully quadratic curved meshes generated by GMSH library. The details of the meshes are given in Table 2. Fig. 7 presents the nondimensionalized turbulent eddy viscosity contours, i.e. μt /μ∞ , using rDG(P1 P2 ) on level 3 hexahedral grid. It is essentially indistinguishable from the reference solution of CFL3D and FUN3D from NASA database. Fig. 8 presents the friction coefficient distribution along the 2D wall bump. By comparison with the reference solution from CFL3D and FUN3D on the finest mesh, i.e. 2 × 1409 × 641, it could be observed that, it is a little difficult to compute the skin friction accurately at the top of bump. Also, indeed there exists some oscillation about the skin friction coefficient distribution, which has been observed by other researchers [40,48]. From their work, DG(P1 ) combined with quadratic curved mesh leads to the oscillation, and higheroder scheme, e.g. DG(P2 ) combined with curved mesh could remove such oscillation. Figs. 9 and Fig. 10 show the grid convergence of lift coefficient CL , and drag coefficient CD for this case.

32


Fig. 1. The overview of the hexahedral grid (level 3)for the case, turbulence flow over a flat plate. Table 2 The quadratic meshes information for the turbulent flow over a 2D bump. Grid no.

Dimensions

Wall spacing (y plus)

No. of points defining bump

Level 1 Level 2 Level 3

2 × 56 × 31 2 × 111 × 61 2 × 221 × 121

4.00e − 5(4.8 ) 2.00e − 5(2.4 ) 1.00e − 5(1.2 )

75 149 297

Fig. 2. The eddy viscosity contours (nondimensionalized by freestream laminar viscosity μ∞ ) using rDG(P1 P2 ) method on level 3 hexahedral grid.

It could be observed that, both DG(P1 ) and rDG(P1 P2 ) could converge to the grid-independent solution , which agrees well with the reference solution from CFL3D and FUN3D. Also, it is quite obvious that, compared with the CFD3D and FUN3D, high-order rDG scheme just requires much less degrees of freedom to converge to the grid-independent solution. Finally, the logarithmic residual convergence histories using DG(P1 ) and rDG(P1 P2 ) on level 3 grid are given in Fig. 11, which shows that the residual of both schemes could almost reduce 7 orders of magnitude, demonstrating that the steady state solutions have been obtained.

Fig. 3. The logarithmic plot of surface friction coefficient distribution over the flat plate obtained by DG(P1 ), rDG(P1 P2 ) and WENO(P1 P2 ) methods on level 3 hexahedral grid. Table 3 The quadratic meshes information for the turbulent flow over a 2D NACA0012 airfoil. Grid no.

Dimensions



2 × 26 × 181 2 × 51 × 361 2 × 101 × 721

3.20e − 5(12.0 ) 1.60e − 5(6.0 ) 8.00e − 6(3.0 )

5.3. Turbulent flow over a NACA0012 airfoil This case considered is a turbulent flow over a NACA0012 airfoil, where Ma∞ = 0.15, and Re∞ = 6 × 106 based on the chod length of the airfoil. The computational meshes employed for this test case are given in Table 3, which are quadratic curved meshes generated by GMSH library. The far field boundary is located al-

most 480 chords away from the airfoil. The level 2 mesh is shown in Fig. 13. This test case could be run at different Angle of Attack (AoA). In our work, only the test cases for AoA = 00 , AoA = 100 , AoA = 150 will be calculated.


Fig. 4. Eddy viscosity ratio distribution at x=0.97, obtained by DG(P1 ), rDG(P1 P2 ), and WENO(P1 P2 ) methods.

33

Fig. 7. The eddy viscosity contours (nondimensionalized by freestream laminar viscosity μ∞ ) using the rDG(P1 P2 ) method along the 2D wall bump.

Fig. 5. The grid convergence of drag coefficient obtained by DG(P1 ), rDG(P1 P2 ), and WENO(P1 P2 ) methods.

Fig. 8. The friction coefficient distribution along the 2D wall bump.

Fig. 6. Convergence histories using DG(P1 ), rDG(P1 P2 ), and WENO(P1 P2 ) on level 2 hexahedral grid.

Fig. 9. Grid convergence for lift coefficient CL .

34


Fig. 12. The nondimensionalized turbulent eddy viscosity using rDG(P1 P2 ) on level 2 hexahedral grid for AoA = 100 .

Fig. 10. Grid convergence for drag coefficient CD .

Fig. 11. Convergence histories using DG(P1 ) and rDG(P1 P2 ) on level 3 hexahedral grid.

First, the numerical results for AoA = 00 are shown in this section. For AoA = 00 , the computed lift coefficient CL is approximately equal to zero. The grid convergence of drag coefficient CD is shown in Fig. 14. The converged values for the lift and drag coefficient are comparable to the reference value of CFL3D from NASA database. It is obvious that, with higher discretization order, the rDG(P1 P2 ) requires much less degrees of freedom to get the grid-independent solution compared with DG(P1 ). Also, the logarithmic residual convergence histories using DG(P1 ) and rDG(P1 P2 ) on level 2 hexahedral grid are presented in Fig. 15, which shows that the residual of both schemes could almost reduce 7 orders of magnitude. Second, the numerical results for AoA = 100 are shown herein. For AoA = 100 , the grid convergence of lift coefficient CL and drag coefficient CD is shown in Fig. 16. The converged values for the lift and drag coefficient are comparable to the reference value of CFL3D from NASA database. It is obvious that, with higher discretization order, the rDG(P1 P2 ) requires much less degrees of freedom to get grid-independent solution compared with DG(P1 ). In addition, for AoA = 100 , the nondimensionalized turbulent eddy viscosity, i.e. μt /μ∞ using rDG(P1 P2 ) on level 3 hexahedral grid, is presented in Fig. 12. The wake is resolved very well using rDG(P1 P2 ).

Fig. 13. The quadratic level 2 hexahedral mesh for the turbulent flow over 2D NACA0012 airfoil.


35

Table 4 Computed lift and drag coefficients for the NACA0012 with Ma∞ = 0.15, AoA = 150 , and Re∞ = 6 × 106 using DG(P1 )and rDG(P1 P2 ) on level 3 hexahedral grid.

Fig. 14. The grid convergence of drag coefficient CD for turbulent flow over 2D NACA0012 airfoil for AoA = 00 .

Schemes

CL

CD

DG(P1 ) rDG(P1 P2 ) Reference value

1.5484 1.5496 1.5446 − 1.5642

0.02130 0.02110 0.02073 − 0.02159

Third, the numerical results for AoA = 150 are shown in this section. For AoA = 150 , the lift coefficient CL and drag coefficient CD on level 3 hexahedral grid are shown in Table 4. The converged values for the lift and drag coefficient are comparable to the reference value range from NASA database. From NASA database, for this test case, the numerical results are from several independent CFD codes: CFL3D, FUN3D, NTS, JOE, SUMB and so on. Thus, there exists one reference solution range, which has been given in Table 4. In addition to the lift and drag coefficient, the friction coefficient Cf distribution at different AoA on the third mesh is also presented in Fig. 17. Our computed solution agrees well with the reference one from CFL3D. Finally, the drag polar for this case is given in Fig. 18. It can be observed that, the numerical results of our rDG schemes match well with the experimental data, demonstrating the accuracy of the developed rDG method for turbulent flows.

5.4. Turbulent flow over a 3D sinusoidal bump in a channel

Fig. 15. Convergence histories using DG(P1 ), rDG(P1 P2 ) on level 2 hexahedral grid for AoA = 00 .

This case is the 3D sinusoidal bump in a channel defined analytically in the NASA database. This is the simplest 3D test case, where the wall bump is smooth and it has no flow separation. The Reynolds number per unit of length is Re∞ = 3 × 106 and the Mach number is equal to 0.2. For this 3D case, the mesh is from official website. The first 3 quadratic curved meshes for DG(P1 ) method are obtained just by simple agglomeration of the linear hexahedral grids of the NASA TMR website [47], without any change of the high-order nodes. The detailed information about the quadratic curved meshes and the original linear ones are given in Table 5. Due to the spanwise variation added for this 3D sinusoidal bump test, in order to guarantee the stabiltiy of rDG(P1 P2 ), it requires a little finer resolution along the spanwise direction. Thus, the fourth quadratic mesh is obtained by agglomeration based on the linear

Fig. 16. The grid convergence of lift coefficient CL and drag coefficient CD for turbulent flow over 2D NACA0012 airfoil for AoA = 100 .

36


Fig. 18. Drag polar on level 3 mesh for the turbulent flow over 2D NACA0012 airfoil.

Fig. 17. The friction coefficient Cf distribution for the turbulent flow over 2D NACA0012 airfoil for different AoA.

Table 5 The quadratic meshes information for 3D turbulent flow over a 3D bump. Grid no.

Dimensions

Original linear meshes


Level Level Level Level

5 × 45 × 21 9 × 89 × 41 17 × 177 × 81 31 × 181 × 41

9 × 89 × 41 17 × 177 × 81 33 × 353 × 161 61 × 361 × 81

1.60e − 5(2.0 ) 8.00e − 6(1.0 ) 4.00e − 6(0.5 ) 2.15e − 6(0.26 )

1 2 3 4

mesh generated by the author. The level 1 quadratic hexahedral mesh is presented in Fig. 19. First, the numerical results based on the first 3 quadratic meshes are presented here using DG(P1 ). The grid convergence of the lift and drag coefficient, CL and CD , is shown in Fig. 20. From Fig. 20a, it could be observed that, DG(P1 ) scheme just requires less degrees of freedom to arrive at the grid-independent CL . And from Fig. 20b, DG(P1 ) scheme is a little superior to CFL3D and FUN3D in terms of number of degrees of freedom. In addition, the gridindependent lift and drag coefficients of our rDG schemes agree well with the reference solution from CFL3D and FUN3D. The convergence of the pressure and viscous drag coefficients, CDp and CDv , using DG(P1 ) is shown in Fig. 21. Second, the numerical results based on the fourth quadratic mesh are presented here using both DG(P1 )and rDG(P1 P2 ). The calculation with this mesh is intended to make a direct comparison between DG(P1 )and rDG(P1 P2 ). The lift and drag coefficients, CL and CD , are shown in Fig. 22. Fig. 22a presents the lift coefficient, CL , which is almost identical for DG(P1 )and rDG(P1 P2 ). Whereas, the drag coefficient, CD , using rDG(P1 P2 ) which is shown in Fig. 22b, is approaching closer to the reference solution from CFL3D and FUN3D. This indicates that, compared with the reference solution, on the same mesh, rDG(P1 P2 ) could provide more accurate solution. In addition, the pressure and viscous drag coefficients of our rDG schemes are also presented in Fig. 23. The convergence histories using different numerical methods on level 4 hexahedral grid are presented in Fig. 24. It can be observed that, the residual of DG(P1 ) and rDG(P1 P2 ) could drop 7 orders of magnitude, demonstrating that the steady state solutions are reached. The overall Cp distribution (axially, on the body along y = 0) on the fourth mesh is given in Fig. 25. It can be observed that our solution agrees well with the refrence solution from FUN3D. 5.5. 3D turbulent flow over a 3D hemisphere cylinder NASA TMR also provides a test case for a turbulent flow over a smooth body of revolution in 3D. This case is designed primarily for numerical analysis of turbulence model simulations; e.g. convergence properties, effect of order of accuracy, etc. The Reynolds number per unit of length is Re∞ = 1.3779 × 107 and the Mach number is equal to 0.6. For this 3D case, the quadratic curved meshes are obtained just by simple agglomeration of the linear


Fig. 19. Level 1 quadratic hexahedral mesh for turbulent flow over a 3D sinusoidal bump in a Channel.

Fig. 20. The grid convergence of lift and drag coefficient using DG(P1 ) on the first 3 quadratic meshes for turbulent flow over a 3D sinusoidal bump in a Channel.

37

38


Fig. 21. The grid convergence of pressure and viscous drag coefficients using DG(P1 ) on the first 3 quadratic meshes for turbulent flow over a 3D sinusoidal bump in a Channel.

Fig. 22. The lift and drag coefficient using DG(P1 )and rDG(P1 P2 ) on the fourth quadratic mesh for turbulent flow over a 3D sinusoidal bump in a Channel.

Fig. 23. The pressure and viscous drag coefficients using DG(P1 ) and rDG(P1 P2 ) on the fourth quadratic mesh for turbulent flow over a 3D sinusoidal bump in a Channel.


Fig. 24. The logarithmic residual using DG(P1 ) and rDG(P1 P2 ) on level 4 quadratic mesh for turbulent flow over a 3D sinusoidal bump in a Channel.

39

achieve the same accuracy compared with DG(P1 ) scheme. However, with mesh refined, it could be observed that, DG(P1 ) scheme outperforms the second-order CFL3D. Also, on level 3 quadratic curved mesh, the solution for WENO(P1 P2 ) is approaching closer to the converged reference value from CFL3D. Fig. 28 shows the pressure coefficient distribution axially. It is obvious that, our numerical solutions agree quite well with the experimental data and the reference solution from CFL3D. Fig. 29 shows the friction coefficient distribution axially. Indeed, there exist a little difference between our numerical solution and the reference solution from CFL3D. However, it is necessary to point out that, the reference solution from CFL3D is based on the finest mesh, i.e. 161 × 289 × 129. Fig. 30 shows the cell-averaged nondimensionalized turbulent eddy vicosity contours using WENO(P1 P2 ) on level 3 quadratic mesh. It could be observed that, the turbulent eddy vicosity grows gradually along the flow direction near the hemisphere wall. Finaly, the logarithmic residual convergence histories using DG(P1 ) and WENO(P1 P2 ) on level 3 grid are given in Fig. 31, which shows that the residual of DG(P1 ) schemes could reduce 7 orders of magnitude. We also note that, similar to that in the first test case, the phenomenon ’convergence stall’ also occurs for this case with WENO(P1 P2 ). And we need to mention that, the phenomenon ’convergence stall’ is about the residual convergence, rather than about the results.

6. Conclusion and outlook

Fig. 25. The pressure coefficient Cp distribution on level 4 quadratic mesh for turbulent flow over a 3D sinusoidal bump in a Channel. Table 6 The quadratic meshes information for the 3D turbulent flow over a 3D hemisphere cylinder. Grid no.

Dimensions

Original linear meshes

Wall spacing


11 × 19 × 21 21 × 37 × 17 41 × 73 × 33

21 × 37 × 17 41 × 73 × 33 81 × 145 × 65

2.76e − 5 1.24e − 5 5.90e − 6

hexahedral grids from the TMR website, without any change of the high-order nodes. The detailed information about the quadratic curved meshes and the original linear ones is given in Table 6. Then level 3 quadratic curved mesh is presented in Fig. 26. Fig. 27 presents the grid convergence of drag coefficient, CD using DG(P1 ) and WENO(P1 P2 ) respectively. It is WENO(P1 P2 ) rather than rDG(P1 P2 ) is used for this calculation. WENO reconstruction on P2 is required to ensure the stability in this case, as the hybrid grids are used, where tetrahedral cells exist at the leading edge. Only the numerical results obtained by WENO(P1 P2 ) on level 3 quadratic mesh are presented. First, it seems that for the coarse mesh, the second-order CFL3D is superior to high-order DG(P1 ) scheme, since it only requires a little less degrees of freedom to

A reconstructed discontinuous Galerkin (rDG(P1 P2 )) method has been developed and extended for the solution of the ReynoldsAveraged Navier–Stokes equations along with the modified oneequation Spalart and Allmaras turbulence model on 3D uniformly refined curved grids. In this paper, two efficient methods, i.e. GMSH library and simple agglomeration, have been used to generate the quadratic curved meshes for high-order rDG(P1 P2 ) method. The numerical experiments demonstrate that both methods could provide satifactory quadratic curved meshes for highorder methods. Especially, simple agglomeration strategy based on linear meshes provide an easy, attractive yet reliable direction to generate the curved meshes. The developed rDG(P1 P2 ) method has been used to compute a variety of turbulent flow problems to assess its accuracy and efficiency with increasing complexity. The numerical experiments indicate that the lift and drag coefficients obtained from rDG(P1 P2 ) solutions agree well with the reference values from CFL3D and FUN3D, demonstrating that the rDG(P1 P2 ) method provides an attractive alternative for computing turbulent flows around complex geometries at a cost slightly higher than its underlying second-order DG(P1 ) method. Second, WENO reconstruction has been adopted to guarantee the linear stability of rDG(P1 P2 ) method, termed as WENO(P1 P2 ). Based on the grid convergence study for the first test case, i.e. Zero Pressure Gradient Flat Plate, it is obvious that, compared with rDG(P1 P2 ) method, WENO(P1 P2 ) can still keep the accuracy of the solution. For the last test case, i.e. 3D Turbulent FLow over a 3D Hemisphere Cylinder, WENO reconstruction has to be used to guarantee the linear stability. The numerical results demonstrate that WENO(P1 P2 ) can provide reliable and accurate solution. Our future work is to conduct a systematic study for a number of test cases to establish whether the rDG methods can actually outperform or have any advantage over the second-order finite volume methods for the compressible turbulent flows.

40


Fig. 26. Level 3 quadratic hexahedral mesh for turbulent flow over a 3D hemisphere cylinder.

Fig. 27. The grid convergence of drag coefficient CD using DG(P1 ) and WENO(P1 P2 ).

Fig. 29. The friction coefficent distribution using DG(P1 ) and WENO(P1 P2 ) on level 3 quadratic mesh.

Fig. 30. The cell-averaged nondimensionalized turbulent eddy vicosity contours using WENO(P1 P2 ) on level 3 quadratic mesh.

Fig. 28. The pressure coefficent distribution using DG(P1 ) and WENO(P1 P2 ) on level 3 quadratic mesh.


Fig. 31. The logarithmic residual convergence histories using DG(P1 ) and WENO(P1 P2 ) on level 3 quadratic mesh.

References [1] Reed WH, Hill T. Triangular mesh methods for the neutron transport equation; 1973. Los Alamos scientific laboratory report. LA-UR73-479. [2] Cockburn B, Karniadakis GE, Shu CW. The development of discontinuous galerkin method. Discontinuous Galerkin Methods, Theory, Computation, and Applications Lecture Notes in Computational Science and Engineering, vol 11; 20 0 0. [3] Bassi F, Rebay S. High-order accurate discontinuous finite element solution of the 2d Euler equations. J Comput Phys 1997;138:251–85. [4] Bassi F, Rebay S. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations. J Comput Phys 1997;131(2):267–79. [5] Warburton TC, Karniadakis GE. A discontinuous Galerkin method for the viscous MHD equations. J Comput Phys 1999;152(2):608–41. [6] Rasetarinera P, Hussaini M. An efficient implicit discontinuous spectral Galerkin method. J Comput Phys 2001;172(2):718–38. [7] Luo H, Baum JD, Löhner R. A discontinuous Galerkin method based on a Taylor basis for the compressible flows on arbitrary grids. J Comput Phys 2008;227(20):8875–93. [8] Dumbser M, Balsara DS, Toro EF, Munz CD. A unified framework for the construction of one-step finite volume and discontinuous galerkin schemes on unstructured meshes. J Comput Phys 2008;227(18):8209–53. [9] Dumbser M, Zanotti O. Very high order PNPM schemes on unstructured meshes for the resistive relativistic MHD equations. J Comput Phys 2009;228(18):6991–7006. [10] Dumbser M. Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier–Stokes equations. Comput Fluids 2010;39(1):60–76. [11] Sudirham JJ, JJW VDV, Damme V, RMJ. Space-time discontinuous Galerkin method for advection–diffusion problems on time-dependent domains. Appl Numer Math 2006;56(12):1491–518. [12] van Leer B., Lo M., van Raalte M. A discontinuous Galerkin method for diffusion based on recoveryAIAA Paper 20 07-408320 07. [13] Luo H, Luo L, Nourgaliev R, Mousseau VA, Dinh N. A reconstructed discontinuous Galerkin method for the compressible Navier–Stokes equations on arbitrary grids. J Comput Phys 2010;229(19):6961–78. [14] Luo H, Xia Y. A class of reconstructed discontinuous Galerkin methods in computational fluid dynamics. In: Proceedings of the international conference on mathematics and computational methods applied to nuclear science and engineering, Brazil; May 2011. [15] Luo H., Xia Y., Nourgaliev R., Cai C. A class of reconstructed discontinuous Galerkin methods for the compressible flows on arbitrary gridsAIAA Paper 2011-01992011. [16] Zhang L, Wei L, He L, Deng X, Zhang H. A class of hybrid DG/FV methods for conservation laws i: basic formulation and one-dimensional systems. J Comput Phys 2012;231(4):1081–103. [17] Zhang L, Wei L, He L, Deng X, Zhang H. A class of hybrid DG/FV methods for conservation laws II: two-dimensional cases. J Comput Phys 2012;231(4):1104–20. [18] Luo H., Xiao H., Nourgaliev R., Cai C. A comparative study of different reconstruction schemes for a reconstructed discontinuous Galerkin method on arbitrary gridsAIAA paper 2011-38392011. [19] Xia Y, Luo H, Nourgaliev R. An implicit hermite WENO reconstruction-based discontinuous Galerkin method on tetrahedral grids. Comput Fluids 2014;96:406–21.

41

[20] Xia Y, Luo H, Frisbey M, Nourgaliev R. A set of parallel, implicit methods for a reconstructed discontinuous Galerkin method for compressible flows on 3d hybrid grids. Comput Fluids 2014;98:134–51. [21] Luo H, Xia Y, Li S, Nourgaliev R. A hermite WENO reconstruction-based discontinuous Galerkin method for the Euler equations on tetrahedral grids. J Comput Phys 2012;231(16):5489–503. [22] Luo H, Xia Y, Spiegel S, Nourgaliev R, Jiang Z. A reconstructed discontinuous Galerkin method based on a hierarchical WENO reconstruction for compressible flows on tetrahedral grids. J Comput Phys 2013;236:477–92. [23] Xia Y., Frisbey M., Luo H., Nourgaliev R. A WENO reconstruction-based discontinuous Galerkin method for compressible flows on hybrid gridsAIAA paper 2013-05162013. [24] Luo H, Luo L, Nourgaliev R. A reconstructed discontinuous Galerkin method for the Euler equations on arbitrary grids. Commun Comput Phys 2012;12(5):1495–519. [25] Luo H, Luo L, Nourgaliev R, Mousseau VA, Dinh N. A reconstructed discontinuous Galerkin method for the compressible Navier–Stokes equations on arbitrary grids. J Comput Phys 2010;229(19):6961–78. [26] Haider F, Croisille JP, Courbet B. Stability analysis of the cell centered finite-volume MUSCL method on unstructured grids. Numir Math 2009;113(4):555–600. [27] Xia Y., Luo H., Nourgaliev R. An implicit reconstructed discontinuous Galerkin method based on automatic differentiation for the compressible flows on tetrahedral grids, AIAA paper 2013-06872013. [28] Liu X, Xuan L, Xia Y, Luo H. A reconstructed discontinuous Galerkin method for the compressible Navier–Stokes equations on hybrid grids. Comput Fluids 2017;152:217–30. [29] Wang L, Kyle Anderson W, Erwin T, Kapadia S. High-order discontinuous Galerkin method for computation of turbulent flows. AIAA J 2015;53(5):1159–71. [30] Geuzaine C, Remacle J. Gmsh: a 3-d finite element mesh generator with built-in pre-and post-processing facilities. Int J Numer Meth Eng 2009;79(11):1309–31. [31] Hartmann R., Leicht T. High-order unstructured grid generation and discontinuous Galerkin discretization applied to a 3d high-lift configuration, AIAA Paper 2015-08192015. [32] Xia Y., Liu X., Luo H. A finite volume method based on a WENO reconstruction for compressible flows on hybrid gridsAIAA Paper 2014-09392014. [33] Bassi F, Rebay S. Discontinuous Galerkin solution of the reynolds-averaged Navier–Stokes and κ − ω turbulence model equations. Comput Fluids 2005;34:507–40. [34] Hartmann R, Held J, Leicht T. Adjoint-based error estimation and adaptive mesh refinement for the RANS and κ − ω turbulence model equations. J Comput Phys 2011;230(11):4268–84. [35] Burgess N., Mavriplis D. Robust computation of turbulent flows using a discontinuous Galerkin methodAIAA Paper 2012-04572012. [36] Ceze M.A., Fidkowski K. High-order output-based adaptive simulations of turbulent flow in two dimensionsAIAA Paper 2015-15322015. [37] Allmaras SR, Johnson FT, Spalart P. Modifications and clarifications for the implementation of the Spalart–Allmaras turbulence model. In: Proceedings of international conference on fluid dynamics; 2012. [38] Crivellini A, DAlessandro V, Bassi F. A Spalart–Allmaras turbulence model implementation in a discontinuous Galerkin solver for incompressible flows. J Comput Phys 2013;241(11):388–415. [39] Nguyen N.C., Persson P., Peraire J. RANS solutions using high order discontinuous Galerkin methodsAIAA Paper 20 07–291420 07. [40] Zhou C., Wang Z. CPR high-order discretization of the RANS equations with the SA modelAIAA Paper 2015-12862015. [41] Luo H, Baum JD, Löhner R. A fast, matrix-free implicit method for compressible flows on unstructured grid. J Comput Phys 1998;146(2):664–90. [42] Tucker PG, Rumsey CL, Spalart PR, Bartels RE, Biedron R. Computations of wall distances based on differential equations. AIAA J 2005;43:539–49. [43] Cheng J, Liu T, Luo H. A hybrid reconstructed discontinuous Galerkin method for compressible flows on arbitrary grids. Comput Fluids 2016;139:68–79. [44] Moro D., Nguyen N.C., Peraire J. Navie–Stokes solution using hybridizable discontinuous Galerkin methodsAIAA Paper 2011-34072011. [45] Xia Y, Liu X, Luo H, Nourgaliev R. A third-order implicit discontinuous Galerkin method based on a hermite WENO reconstruction for time-accurate solution of the compressible navierstokes equations. Int J Numer Methods Fluids 2015;79(8):416–35. [46] Arnold DN, Brezzi F, Cockburn B, Marini L. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J Numer Anal 2002;39(5):1749–79. [47] Rumsey C. Turbulence modeling resource, NASA Langley Research Center, 2011, http://turbmodels.larc.nasa.gov/. [48] Schrock C.R., Benek J.A., Galbraith M.C., Knapke R.D., Turner M. Evaluation of a discontinuous Galerkin implementation of RANS and Spalart Allmaras turbulence modelAIAA Paper 2015-00562015. [49] Kroll N, Bieler H, Deconinck H, Couaillier V, van der Ven H. ADIGMA - A European initiative on the development of adaptive higher-order variational methods for aerospace applications. Notes on numerical fluid mechanics and multidisciplinary design, 113. Springer; 2010. [50] Kroll N, Hirsch C, Bassi F, Johnston C, Hillewaert K. IDIHOM: industrialization of high-order methods – a top-down approach. Notes on numerical fluid mechanics and multidisciplinary design, 128. Springer; 2015.