High-order discontinuous Galerkin solutions of three-dimensional incompressible RANS equations

Computers & Fluids 81 (2013) 122–133

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

High-order discontinuous Galerkin solutions of three-dimensional incompressible RANS equations Andrea Crivellini a,⇑, Valerio D’Alessandro a, Francesco Bassi b a b

Dipartimento di Ingegneria Industriale e Scienze Matematiche, Università Politecnica delle Marche, Via Brecce Bianche, 60100 Ancona (AN), Italy Dipartimento di Ingegneria Industriale, Università di Bergamo, Viale Marconi 5, 24044 Dalmine (BG), Italy

a r t i c l e

i n f o

Article history: Received 31 March 2012 Received in revised form 20 March 2013 Accepted 17 April 2013 Available online 2 May 2013 Keywords: Incompressible flow Discontinuous Galerkin method Artificial compressibility Riemann solver RANS equations

a b s t r a c t This paper presents the latest developments of the artificial compressibility flux Discontinuous Galerkin (DG) method introduced in [1], extended in [2] to natural convection flows, in [3] to unsteady flows and, more recently, in [4] to turbulent flows. Here we consider the three-dimensional incompressible Reynolds Averaged Navier–Stokes equations (RANS) coupled with the Spalart–Allmaras (SA) turbulence model. The development of efficient high-order RANS solvers is still a difficult task due to the extreme stiffness of the governing equations. For this reason the turbulence model here has been suitably modified, in the source terms and in the diffusion coefficient, in order to prevent unphysical conditions of the turbulent working variable and of one of the closure functions, which sometimes result in numerical instabilities. The reliability, accuracy and robustness of the method were assessed by computing several test cases in simple and real-life configurations: the flow over a sinusoidal bump, the flow field past a sphere in the supercritical regime, the flow field past a delta wing, and the flow around the DLR-F6 wing body transport configuration. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, DG methods have received increasing interest from the CFD community because they combine flexibility, stability, robustness, and accuracy. DG methods are based on polynomial approximations inside the computational elements with no global continuity requirement. As in continuous finite elements methods, DG methods can increase the order of the accuracy of the solution by increasing the degree of the polynomial approximation, and its accuracy is not affected by the use of the unstructured hybrid meshes usually adopted in complex and real-life problems. The discontinuous approximation between neighbouring elements is treated, as in finite volume methods, using numerical fluxes and the resulting compactness of the algorithm makes the method easily and efficiently implementable on parallel architectures. Despite these attractive advantages, the suitability of DG methods for analyzing complex engineering problems remains a challenging task. Complex fluid flow problems often involve turbulence. In the field of numerical computation of turbulent flows, Direct Numerical Simulation (DNS) and Large-Eddy Simulation (LES) are still prohibitive for many flows of engineering interest due to their large ⇑ Corresponding author. Tel.: +39 0712204436. E-mail addresses: [email protected] (A. Crivellini), [email protected] (V. D’Alessandro), [email protected] (F. Bassi). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.04.016

computational requirements. Thus the RANS approach, having a limited extra-computational load, remains the only way to analyze this kind of flow field. In this paper, several high-order solutions of three-dimensional problems, in simple and complex geometries, are presented for the turbulent and incompressible flow regime. All the computations were performed by a DG solver using a new stable and robust implementation of the RANS-SA equations. In the literature, only Lubon et al. [5], dealing with Detached Eddy Simulation (DES) based on the SA equation, have reported 3D DG solutions of the same model. Moreover the DG solutions of the RANS equations are limited to the compressible case; as far as the authors know, the present paper reports the first incompressible DG RANS solver in the literature. The SA model is here chosen since it allows good accuracy and robustness for aerodynamic flows, as demonstrated by several research papers, introducing only one additional equation. The model is also particularly attractive since the turbulent variable requires a resolution less than that needed to capture the velocity profiles at the wall. For a moderate high-order solution, thanks to the linear behaviour of the working variable close to the wall, the dimensionless first cell height, yþ c , can be of order 10, see [4]. Nevertheless, solving RANS equations using DG methods is still a difficult task due to the numerical stiffness of the problem coming from the highly non-linear nature of the source terms in the

123

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

turbulence model equation and from the grid stretching needed to solve the turbulent variable at the wall. Thus, in order to increase the stability and efficiency of the solver, the RANS and the turbulence model equations need to be handled with some form of limiting procedure or with the use of stabilizing terms. At the time of this writing, the authors are not aware of more than a few DG implementations of the RANS equations coupled with two turbulence models: the one-equation SA and the twoequation k–x models. Bassi et al. [6] solved the RANS and k–x equations with some modifications in the original model consisting of (i) rewriting the x equation in terms of log(x) instead of x, (ii) fulfilling the realizability constraints for the turbulent stresses (obtaining a lower bound for x), and (iii) enforcing a zero lower bound for the turbulent kinetic energy in the source terms and in the turbulent viscosity constitutive relation. Peraire et al. [7] solved the RANS and the SA turbulence model equations using an artificial viscosity term in order to stabilize the discretization of the turbulence model equation. Oliver and Darmofal [8,9] also developed a DG solver for the RANS equations and the SA model. They added an artificial dissipation, governed by an additional PDE, to avoid numerical oscillations in under-resolved regions. Landmann et al. [10] solved the RANS equations coupled with both the SA and k–x turbulence models. For the k–x model, they used the same approach adopted by Bassi et al. [6], while for the SA model, a strong instability for negative values of eddy viscosity was observed. Therefore they implemented a limiting technique that after each Newton iteration resets to zero negative values of the computed turbulence variable. Finally, Hartmann et al. [11] developed a DG code for turbulent flow computations using the same modifications proposed by Bassi et al. in [6] for the k–x equations. Here the governing equations are solved using the artificial compressibility flux DG method introduced in [1]. This scheme introduces a perturbation of the inviscid problem at the interface level in order to recover hyperbolicity of the system of equations and to obtain the numerical flux as the solution of the associated Riemann problem. As demonstrated in [1], the method guarantees, for polynomial approximations of degree k, a convergence rate of k + 1 for the velocity components and of k for the pressure. In [4] the method proved also to be extremely well suited for twodimensional RANS solutions using highly stretched and curved computational elements. This paper is organized as follows. Section 2 presents the governing equations and outlines the turbulence model modifications. Section 3 briefly describes the DG space discretization and the implicit solution algorithm. Finally Section 4 reports the numerical results of several test cases.

2. Governing equations The complete set comprising the incompressible RANS equations and the Spalart–Allamaras turbulence model can be written as follows:

@uj ¼ 0; @xj @ui @ @ þ ðui uj þ pdij Þ  ½2ðm þ mt ÞSij  ¼ 0; @xj @t @xj   ~ ~ @m @ 1 @ @m ~Þ  ¼ s; þ ðuj m n @t @xj r @xj @xj

ð1Þ

where p = P/q is the pressure divided by the density, m is the kinematic viscosity, dij is the Kronecker delta, Sij is the strain-rate tensor

  1 @ui @uj ; Sij ¼ þ 2 @xj @xi

and mt is the turbulent viscosity which is computed according to the m~ variable by

mt ¼ fv 1 maxð0; m~Þ:

ð2Þ

The diffusion coefficient n and the source terms s for the SA equation will be defined below. Even if the SA model was originally ~, its numerical solution can developed only for positive values of m admit unphysical negative values which often produce numerical instabilities and a blow-up of the computation. ~ values are Two different approaches for handling negative m implemented in our code. The first one relies on unpublished work of Allmaras himself and was implemented in a compressible DG solver by Oliver [9,12]. Recently it was successfully used by Crivellini et al. in [4] in the two-dimensional version of this incompressible DG solver. Following this approach, the diffusion coefficient and source terms are defined by

(

mð1 þ vÞ v P 0;   m 1 þ v þ 12 v2 v < 0;

n¼

s¼

 8  2  < mdv cb1 þ crb2 2  c w1 fw k r

cb2 : c Sm m~2 b1 ~g n þ c w1 2 þ r d

ð3Þ

~ @m ~ @m @xj @xj

~ @m ~ @m @xj @xj

v P 0;

ð4Þ

v < 0;

2

v where g n ¼ 1  1000 ; d is the minimum distance from the wall, 1þv2 and the v variable is the dimensionless eddy viscosity:

m~ v¼ : m

ð5Þ

~ modifications were built taking in consideration These negative m the continuity of both the function and the first derivative with re~, i.e. the resulting Jacobian matrix is continuous even at spect to m m~ ¼ 0. The second approach tried here is inspired by the DG implementation of the k–x turbulence model in Bassi et al. [6]. The ~ variable is manipulated with the same approach equation for the m used there for the turbulent kinetic energy. Thus, whenever a neg~ value is reached, the source term is set to zero while the difative m fusion coefficient is fixed at the molecular viscosity value:

 n¼

mð1 þ vÞ v P 0; m v < 0;

(  2  cb1 mv s¼

ð6Þ

  cw1 fw þ crb2 0 v < 0: k2 r

d

~ @m ~ @m @xj @xj

v P 0;

ð7Þ

Unlike the previous case, using the Bassi approach the Jacobian of ~ ¼ 0. the source and diffusion terms are no longer continuous at m To completely define the PDE system, the following closure functions are now introduced. 3

fv 1 ¼ v3 vþc3 ; ð v1 Þ

f v 2 ¼ 1  ð1þvvfv 1 Þ ;

h 6 i16 1þc g ¼ r þ cw2 ðr 6  rÞ; f w ¼ g g6 þcw3 ; 6 w3   r max r < 0; r ¼ m~2 2 ; r¼ eSk d minðr  ; r max Þ r  P 0;

ð8Þ

where rmax is a defined positive constant and e S is a function of both ~: the vorticity magnitude, S, and of the turbulent working variable m

m~ e S ¼ S þ 2 2 fv 2 ; k d

Xij ¼

  1 @ui @uj ;  2 @xj @xi

S¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Xij Xij :

ð9Þ

Finally, in order to completely define the SA model, the standard closure constants are introduced:

124

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

cb1 ¼ 0:1355; cb2 ¼ 0:622; cv 1 ¼ 7:1; r ¼ 2=3; cw1 ¼ ckb12 þ ð1þcr b2 Þ ; cw2 ¼ 0:3;

cw3 ¼ 2;

ð10Þ

k ¼ 0:41:

m~ v ¼ ; 2 e SðkdÞ2 SðkdÞ m þ vfv 2

and, since fv2 is negative for 1.003 6 v 6 18.4, when S(kd)2/m 6 6.008 the r function can have a negative value. Moreover, the r function can be safely limited because, quoting the original Spalart–Allmaras paper [13], ‘‘the region r > 1 is exercised only in adverse pressure gradients, and then rarely beyond r = 1.1’’. All things considered, the new definition consists in limiting r to a positive value, rmax, and it handles negative and large values in the same way because the original r, for a positive v, is negative only if it goes through a singular point with r ? ±1. Furthermore, the new definition of the r function avoids the change of sign of the production and destruction source terms that can sometimes produce numerical instabilities. A detailed analysis of the behaviour of the source terms, as functions of r and v, is given in [4]. The new r formulation, as will be more clear in the numerical results section, greatly improves the stability of our DG solver without entailing significant differences in the computed results. At last note that with meaningful, thus positive, values of ~ and r, the resulting turbulence model is the standard SA both m model without trip term, thus implemented in a fully turbulent mode.

def

sv t ¼ v þ  nþ þ v   n ;

ð12Þ

Here, u 2 RM denotes the state vector of M primitive variables, ^ 2 RM is ~gT ; s 2 RM is the vector of source terms, and u u ¼ fp; ui ; m ^ ¼ fp=a2 ; ui ; m ~gT where a2 is a standard artificial comdefined by u pressibility coefficient. The vectors Fc ; Fv 2 RM  RN denote the inviscid and viscous flux functions, respectively, and N is the spatial dimension. The entries in s, Fc, Fv can be found by comparison with Eq. (1). Even if the above equation appears to be a standard artificial compressibility method note that, as will be clear later, the artificial compressibility flux DG method is suitable also for pure incompressibity, thus with a2 ? 1. In order to construct a DG discretization of the PDE equations here considered, we introduce a triangulation T h ¼ fKg of an approximation Xh of X, that is we partition Xh into a set of ne non-overlapping elements K (not necessarily simplices). We denote by E 0h the set of internal element faces, by E @h the set of boundary element faces, and by E h ¼ E 0h [ E @h their union. We moreover set

[ e2E 0h

e;

C@h ¼

[

e;

Ch ¼ C0h [ C@h :

def

sqt ¼ qþ nþ þ q n ;

ð15Þ

where the ± superscripts indicate the traces of the solution over an edge shared by two adjacent elements K+ and K, see Fig. 1. Similarly, for all scalar, vector or tensor quantities such that a (possibly two-valued) trace is available on e, we introduce the average operator

ðÞþ þ ðÞ : 2

def

fg ¼

ð16Þ

The DG formulation of Eq. (12) then requires finding uh1 ; . . . ; uhM 2 Uh such that

Z

/h

^h @u dx  @t

Z

$h /h  Fðuh ; $h uh þ rðsuh tÞÞdx þ

Z

s/h t Ch

Xh

  b F uh ; ð$h uh þ ge re ðsuh tÞÞ dr Z ¼ /h sðuh ; $h uh þ rðsuh tÞÞdx

3. Discontinuous Galerkin discretization The complete set of equations for the incompressible RANS with the SA turbulence model are implemented in our code as follows, where we use a compact form.

ð14Þ

for some polynomial of degree k P 0, where denotes the space of polynomials of global degree at most k on the element K. Note that taking advantage of the flexibility ensured by DG methods, unstructured hybrid grids (consisting of hexahedral, pyramidal, tetrahedral and prismatic elements) can be easily handled adopting the same polynomial Pk ðKÞ space for all the elements. In our code we have implemented a hierarchical and orthogonal set of shape functions, defined in the physical space, that should overcome the ill-conditioning of the elemental mass matrices for higher-order polynomials of high aspect ratio and curved elements. This set is obtained using a modified Gram–Schmidt procedure assuming as a starting point a set of monomial functions. When dealing with DG discretizations, it is also convenient to introduce the jump, Eq. (15), trace operators. In particular for all vector quantities v and scalar quantities q on a generic internal face e 2 E 0h we define

Xh

C0h ¼

component

def

ð11Þ

^ @u þ $  Fc ðuÞ þ $  Fv ðu; $uÞ ¼ sðu; $uÞ @t

each

uhi 2 Uh ¼ f/h 2 L2 ðXÞ : /h jK 2 Pk ðKÞ8K 2 T h g

Looking at the closure functions of Eq. (8), it is quite easy to see that in our high-order SA model implementation, the function r is modified by introducing the realizability condition proposed in [4]. In the original SA model, this quantity was defined to be r = l2/(k d)2, where l is a mixing length inspired by the algebraic model, thus it should be always positive. However, the original formulation of r can be written as

r¼

sume the following space setting for uhi ¼ uh1 ; . . . ; uhM of the numerical solution uh:

ð17Þ

Xh

for each function /h 2 Uh. F defines the sum of the inviscid and viscous numerical flux functions. The lifting operator re, which is assumed to act componentwise on the jumps of uh, is defined as the solution of the problem

Z

/h  re ðv Þdx ¼ 

Z

f/h g  v dr;

8/h 2 ½Uh N ;

v 2 ½L1 ðeÞN ;

e

Xh

ð18Þ and the function r is related to re by the equation def

rðv Þ ¼

X

re ðv Þ:

e2E h

ð13Þ

e2E @h

The solution is approximated on T h by a piecewise polynomial function, possibly discontinuous on element interfaces, i.e. we as-

Fig. 1. Normals and local frame at quadrature point P on edge e.

ð19Þ

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

In our DG method, the inviscid and viscous parts of the numerical F are treated independently. Taking in consideration the pure flux b incompressible case, a2 ? 1, the key idea in handling the inviscid numerical flux is to reduce the problem of the flux computation to the solution of a Riemann problem as in the compressible case. In order to recover the hyperbolicity of the equations, the incompressibility constraint is relaxed by adding an artificial compressibility term to the continuity equation. At each quadrature point P on Ch we therefore solve the Riemann problem for the equations

1 @p @u þ ¼ 0; c2 @t @x @u @ðu2 þ pÞ þ ¼ 0; @t @x @ v @uv þ ¼ 0; @t @x ~ @um ~ @m þ ¼ 0; @t @x

ð20Þ

with initial datum

( u¼

uh

if x < 0;

uþh

if x < 0;

where x denotes a locally defined axis oriented as the normal vector n+ pointing out of K+ and located in such a way that x = 0 at P (see Fig. 1), u is the velocity component along the x axis, and v is the vector gathering the tangential components of the velocity field. Denoting by u⁄ the solution of the Riemann problem on the space–time line x/t = 0, we finally set

  b F c uh ¼ Fc ðu Þ: The details of the procedure adopted to determine the state u⁄ for the general case of the coupled Navier–Stokes and a general scalar equation are thoroughly discussed in Appendix A of [1]. The numerical viscous flux is given by

  def b F v uh ; ð$h uh þ ge re ðsuh tÞÞ ¼ fFv ðuh ; $h uh þ ge re ðsuh tÞÞg;

ð21Þ

where ge is a stability parameter which must be greater than the number of faces of the elements, see e.g. [14]. All integrals appearing in the spatially discretized problem Eq. (17) are computed by means of appropriate Gauss integration rules. The theoretical degree of integration must take into account the degrees of (i) the polynomial approximation, (ii) the mapping, and (iii) the Jacobian of the mapping. With the purpose of reducing the computational effort of this operation, cheaper non-product formulae, if the required degree is available in [15], are preferred to the tensor product ones. 3.1. Boundary conditions The DG discretization is best suited for a weak enforcement of boundary conditions. This can be easily achieved by properly defining a boundary state which, together with the internal state, allows computing the numerical fluxes and the lifting operator on the portion C@h of the boundary Ch. For Dirichlet-type boundary data this can easily be achieved by properly defining a boundary state ubh which, together with the internal state uþ h , is suitable for the evaluation of the the numerical F c and b F v , and the function rþ fluxes, b e . Moreover the averages f$h uh g; fre ðsuh tÞg are set in b F v equal to the internal values. The wall-type boundary conditions have been implemented by defining the on the exterior of the boundary faces as

boundaryþstate  ~ ~bh ¼ $m ~þ ubh ¼ pþ ; uþ and $m i ; m h . In this case the external boundary state exactly replaces u h in the jump operators, in the numerical fluxes and in the lifting operator.

125

3.2. Time discretization and linear system solution The discrete problem corresponding to Eq. (17) can be written as

M

dU þ RðUÞ ¼ 0; dt

ð22Þ

where U is the global vector of unknown degrees of freedom and M is the global block diagonal mass matrix This is a system of (nonlinear) ODEs that, in many applications, requires implicit schemes in order to be solved efficiently. If a backward Euler scheme is adopted, the following equation is obtained:



M @RðUn Þ þ DU ¼ RðUn Þ: @U Dt

ð23Þ

@RðUn Þ The Jacobian matrix @U can be regarded as an NK  NK block sparse matrix where NK is the number of elements in T h and the rank of each block is M  N KDOF , where N KDOF is the number of degrees of freedom for each variable in the generic element K. Thanks to the DG discretization here adopted, the degrees of freedom of a generic element K are only coupled with those of the neighbouring elements and the number of nonzero blocks for each (block) row K of the Jacobian is therefore one more than the number of elements surrounding the element K. The Jacobian matrix of the DG discretization is computed analytically and therefore, even with very large time steps, the method can achieve quadratic convergence in the computation of steady state solutions because Eq. (23) is identical, for Dt ? 1, to one iteration of the Newton method applied to the steady discrete problem. For steady state computation, the growth of the local (to the grid element) CFL number is related to the residual vector L2 and L1 norms by

CFL ¼ min f ¼ max

  CFLmin ; CFL max fb

! kRðUn ÞkL2 kRðUn ÞkL1 ; ; kRðU0 ÞkL2 kRðU0 ÞkL1

ð24Þ

where CFLmin, CFLmax and b are user-defined input parameters. To solve Eq. (23), we resort to both the matrix-based and matrix-free preconditioned GMRES (Generalized Minimal RESidual) linear solvers available in the PETSc library [16], the software upon which our DG code relies for the purpose of parallelization (see [17]). The SPMD (single process multiple data) strategy is based on a grid partitioning accomplished by means of the METIS package, [18]. Each processor owns the data related to its local portion of the grid and the data on remote processors are accessed through the MPI paradigm. Thanks to the compactness of our DG method, the only the data owned by the near neighbour elements at the partition boundaries needs to be shared by different processes. The final consideration of this section regards the a2 value in Eq. (12). When a2 – 1, a time derivative for the pressure is introduced in the continuity equation, but here only with the aim of reducing the matrix condition number. The use of a finite a2 value is particularly relevant during the first pseudo time steps when, since the residuals are large, Dt is small and the number of Krylov iterations required to reach the prescribed convergence level for the liner system can decrease significantly. However, once again, we want to point out that this algorithm simply defines, together with Eq. (24), the globalization strategy of the Newton–Krylov algorithm and it is not a standard artificial compressibility method. The a coefficient can safely differ from the artificial compressibility coefficient c, moreover the code can work with a2 ? 1 and this means with the block corresponding to the pressure degree of freedoms null inside the M mass matrix.

126

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

4. Numerical results This section presents several high order DG computations of steady-state problems in simple and in realistic configurations. All the numerical tests were performed in order to asses the accuracy, the reliability, and the robustness of the proposed RANS-SA implementation. All the solutions were obtained on distributed memory parallel machines: the computations with lower memory consumption were run on a Linux Cluster, with eight AMD Opteron based nodes for a total of 64 CPU cores operating at 2.3 GHz, while the more expensive ones were performed on an IBM SP6 supercomputing facility at CINECA. All the meshes here used have quadratic mappings of the element faces, thus second order geometric computational elements are employed. The linear system derived from the DG discretization was solved using the restarted GMRES algorithm preconditioned with an additive Schwarz method (ASM), with one strip of overlapping elements, adopting an incomplete lower–upper ILU(0) factorization of the partitioned matrix owned by each process. For the iterative solver, the following properties were set: a Krylov subspace of 120 vectors, 240 maximum inner iterations, and a relative residual convergence equal to 105. The initial flow field of all the P1 solutions was set to uniform ^ , Eq. (12), and CFLmin = 101, freestream conditions with a = 10 in u 20 CFLmax = 10 , and b = 1 in Eq. (24). For computational efficiency, each converged lower-order solution was taken as the initial flow field of the next higher-order solution. Thanks to the good initial field approximation, the CFLmin parameter could be gradually increased to 10, raising the degree of the polynomial approximation. This is a very important behavioural feature, since the greater number of Newton steps is requested at the lower polynomial orders when a lower computational load is involved. The artificial compressibility parameter of the Riemann solver is here always set to 1, however, as can be seen in [4], the results are not significantly affected by the c2 value. Even the choice of the ge value in Eq. (21) has, in our computational experience, little effect on the results and a value smaller of theoretical stability limit can be safely adopted. Here, the results were obtained with ge ranging

from 4, for almost all the computations, to 9, for the DLR-F6 test case. Finally, if not stated otherwise, the following results were obtained using the modified r closure function, with rmax = 10, and ~ case. the Bassi approach for the negative m 4.1. Flow over a 3D sinusoidal bump in a channel (Re = 3  106) The first case considered is the 3D sinusoidal bump in a channel defined analytically in Eq. (25) and represented in Fig. 2a; this is the simplest test case since, the wall bump being smooth, it has no flow separation. The wall bump, extending from x = 0 to x = 1.5 at the two sides of the computational domain y = 0 and y = 1, is defined by

8 0 6 x0 < 0:3; > < 0;   x0  p 4 ; 0:3 6 x 6 1:2; z ¼ 0:05  sin p0:9  3:0 0 > : 0; 1:2 < x0 6 1:5;

ð25Þ

x ¼ x0 þ 0:3  ½sinðpyÞ4 ; where x0 is any given x location on the plane located at y = 0. Symmetry conditions are imposed, as well as in all other boundaries, between the inflow/outflow conditions and the wall condition at the bump. The Reynolds number per unit of length is Re = 3  106. There is no experiment associated with this flow field, i.e. it is a purely numerical verification problem. It was developed by Rumsey in [19] as part of a NASA effort aimed at improving the consistency, verification, and validation of turbulence models widely used in the aerospace community. Four grids of hexahedral elements were used for this test case: the first two have a different near wall spacing, with ne = 4550, while the other two are finer, ne = 9790 and ne = 18,200, respectively. In order to solve accurately the boundary layer profile, the grids were built stretching the elements near the wall bump: one  coarse grid is built to obtain O yþ 10, while all the other grids c have a yþ c value of an order of magnitude of about 1. The dimenus þ þ sionless firstpcell ffiffiffiffiffiffiffiffiffiffiffiheight yc value is here evaluated as yc ¼ m yc , where us ¼ sw =q is the friction velocity, sw is the viscous stress tangential to the wall and yc is the centroid distance of the element next to the wall.

Fig. 2. Three-dimensional sinusoidal bump, Re = 3  106.

127

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

Fig. 3. Three-dimensional sinusoidal bump, Re = 3  106, ne = 18,200. P3 solution.

This case was computed up to P3 for the finest grid, see Fig. 3 for a representation of the result obtained, while for all the other grids the solution was computed up to P4 . The DG results were compared with those reported in [20], of the well established finite volume code CFL3D on a very fine 65  705  321 structured grid. Fig. 2b, showing the pressure distribution on the bump, proves that the reference solution and our results are almost indistinguishable. This is a crucial issue for the proposed high-order implementation of the RANS-SA equations since the negative modifications almost do not effect the final solution and they only improve the robustness of the code. Tables 1 and 2 report the force coefficients for all the computations performed. The DG and the CFL3D codes produce very close results, confirming the irrelevance of the turbulence model manipulations to the final solution. Note that only in this test case, involving no flow separation, can the results with the original r closure function be easily computed. For this reason, they are shown here. However the Tables show that the results obtained with the new and the original r definition are really close to each other. Finally, Tables 3 and 4 report the force coefficients obtained with the two ne = 4550 grids. These results confirm what we verified in [4] as regards the turbulent flat plate problem: for high-order computations, the dimensionless height yþ c of the elements next to the wall can be taken to be of order 10.

Table 1 Three-dimensional sinusoidal bump, drag coefficient – CFL3D: Cd = 0.003589. Cd

P1

P2

P3

P4

~ ne = 4550 Standard r, Allmaras m ~ ne = 9790 Standard r, Allmaras m ~ ne = 9790 Modified r, Allmaras m ~ ne = 9790 Modified r, Bassi m ~ ne = 18,200 Standard r, Allmaras m

0.004427 0.004266 0.004233 0.004232 0.003493

0.003188 0.003173 0.003157 0.003165 0.003464

0.003568 0.003566 0.003567 0.003568 0.003619

0.003580 0.003580 0.003581 0.003578 –

Table 2 Three-dimensional sinusoidal bump, lift coefficient – CFL3D: Cl = 0.025005. Cl

P1

P2

P3

P4

~ ne = 4550 Standard r, Allmaras m ~ ne = 9790 Standard r, Allmaras m ~ ne = 9790 Modified r, Allmaras m ~ ne = 9790 Modified r, Bassi m ~ ne = 18,200 Standard r, Allmaras m

0.026218 0.026560 0.026410 0.026474 0.025277

0.024956 0.024961 0.024939 0.024946 0.024594

0.024633 0.024634 0.024630 0.024634 0.024564

0.024559 0.024560 0.024560 0.024558 

Table 3 ~, ne = 4550. Three-dimensional sinusoidal bump, drag coefficient – modified r, Bassi m Cd   ’1 O yþ  c ’ 10 O yþ c

P1

P2

P3

P4

0.004401

0.003180

0.003568

0.003578

0.002911

0.003130

0.003396

0.003546

Table 4 ~, ne = 4550. Three-dimensional sinusoidal bump, lift coefficient – modified r, Bassi m Cl   ’1 O yþ  c ’ 10 O yþ c

P1

P2

P3

P4

0.026148

0.024940

0.024633

0.024558

0.025748

0.024980

0.025071

0.024764

4.2. Flow past a sphere (Re = 1.14  106) The flow field around a sphere, at the supercritical regime Re = 1.14  106, was investigated, see Fig. 4 for a representation of our solution, and the results compared with experimental measurements performed by Achenbach in [21] and DES computations of Costantinescu and Squires [22]. The solutions were computed on two grids, consisting of 2994 and 11,456 hexahedral elements, built, since two symmetry conditions were adopted, only on a quarter of a spherical domain. These grids were generated by a process of radial extrusion from the anisotropic unstructured quadrilateral elements on the surface of the sphere and they differ also in the radius of the far-field condition: 100 and 30 times the sphere diameter, respectively. The first cell height at the wall defines a maximum yþ c value of about 10 in both cases. The solutions were computed, using 128 CPU cores, up to the P6 polynomial approximation on the ne = 2994 grid and up to the P5 on the ne = 11,456 grid. This test case was selected to validate the 3D DG solver even if, thanks to its geometric simplicity, the solution is expected to be axisymmetric. However, the flow field that develops around the sphere possesses two aspects that are difficult to accurately compute: the massive separations and the transition from laminar to turbulent inside the boundary layer. As regards the first topic, numerical experiments have revealed ~ modifications results in a degraded that the use of only the negative m convergence history which can be significantly improved by adopting the new definition of r in Eq. (8). As a matter of fact, we have verified that, in the latter case, the Newton steps required to achieve full convergence up to P4 are more than halved compared to what is needed with the standard formulation of r, see Fig. 5.

128

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

Fig. 4. Sphere, Re = 1.14  106, ne = 2994. P6 solution.

Fig. 5. Sphere, Re = 1.14  106, ne = 11,456. Residuals convergence history.

As already mentioned, the flow field past a sphere in the supercritical regime is characterized by a boundary layer which becomes turbulent before its detachment from the wall surface. Thus it is quite easy to understand as a simple treatment of the boundary layer, fully turbulent from the stagnation point, can lead to relevant differences between the measured and computed pressure and skin friction coefficients. Thus the effect of the laminar to turbulent transition was investigated and compared to the fully turbulent solution. Here the transition was introduced in the simplest way, i.e. specifying the position at which the source terms of the SA equation become active. Let (r, u, h) be the radial, polar, and azimuthal direction, respectively. Then the transition point was forced in the same location as was observed in the experimental measurements, i.e. u = 92.5°. Fully turbulent computations accurately predict the pressure distribution up to the minimum, but downstream, the separation results are not in very good agreement with the experimental measurements, Figs. 6a and 7a . The effect of forcing the transition at a specified location does not produce important modification of the pressure distribution, Figs. 6b and 7b. Only a quite small reduction around the suction peak and a small growth in the separated region are in fact noticeable for cp. In contrast, the skin friction coef-

ficient distribution is really sensitive to the transition: Fig. 8a and b display how the fully turbulent treatment of the boundary layer leads to a large discrepancy between the numerical results and the experimental measurements and the agreement is clearly improved by forcing the transition. Tables 5 and 6 report, with or without the imposed turbulent transition, the drag coefficients obtained, for both the grids, at different polynomial orders. It is quite easy to observe that, for each case, these values tend to converge. The fully turbulent treatment in the computations yields, for the finest solutions, a drag coefficient Cd = 0.10804 on the grid with ne = 11,456 and a Cd = 0.10806 value on the grid with ne = 2994. These values are lower than the drag measurements, Cd = 0.12. Nevertheless these results are in good agreement with DES computations, using the SA model close to the wall, performed by Costantinescu and Squires, [22], since they have reported Cd = 0.102 ± 0.011. As observed in the same paper, fixing the transition produces a slightly delayed boundary layer separation. This behaviour contributes to a lower drag force, hence the agreement between the experimental and numerical data is less satisfactory, but this result is still in good agreement with the time averaged DES results obtained when forcing the same transition point.

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

129

Fig. 6. Sphere, Re = 1.14  106, ne = 2994. Pressure coefficient.

Fig. 7. Sphere, Re = 1.14  106, ne = 11,456. Pressure coefficient.

Table 5 Sphere, Re = 1.14  106. Fully turbulent computation. Exp. data Cd = 0.12, [21]; fully turbulent DES: Cd = 0.102 ± 0.011, [22]. Cd

P1

P2

P3

P4

P5

P6

ne = 2944 ne = 11,456

0.29723 0.21246

0.13362 0.11086

0.11245 0.10914

0.10885 0.10798

0.10737 0.10804

0.10806 

Table 6 Sphere, Re = 1.14  106. Laminar-to-turbulent transition at u = 92.5°. Exp. data Cd = 0.12, [21]; DES-imposed transition point: Cd = 0.080 ± 0.0075, [22]. Cd

P1

P2

P3

P4

P5

P6

ne = 2944 ne = 11,456

0.21479 0.14947

0.10590 0.07169

0.07382 0.06914

0.07591 0.06770

0.06654 0.06772

0.06637 –

4.3. Flow field past VFE-2 configuration (Re = 3  106, a = 13°) In this subsection, we discuss the results concerning the flow past a delta wing, defined by Chu and Luckring in [23], in a configuration called VFE-2 that corresponds to the geometry of a 65° sweep angle and to a blunt leading edge. Here we consider the wing operating at Re = 3  106 with an angle of attack of a = 13°; the solutions will be compared with experimental measurements performed by Konrath et al. [24] at M1 = 0.4. The computations were performed, up to P5 , on a grid with 13,816 hexahedral elements obtained by agglomerating 2197 elements of a multi-block structured grid supplied by Crippa [25]. The resulting first cell height at the wall guarantees a maximum yþ c of the order of 10. The P1 and P2 computations were run using 16 CPU cores, the P3 solution was obtained using 28 CPU cores, and for the P4 and P5 polynomial approximations we used 128 CPU

130

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

Fig. 8. Sphere, Re = 1.14  106. Skin friction coefficient.

Fig. 9. VFE-2, Re = 3  106, a = 13°. Convergence histories.

cores. In this case study, convergence up to the machine precision ~ can be reached only using the Bassi negative modifications for m coupled with the modified r function, Eq. (8). The use of the Allmaras corrections, for instance, does not permit the full convergence of the P1 approximation, as shown in Fig. 9. Fig. 10 represents our P5 solution while Fig. 11 shows the comparison between the numerical and experimental data for several of the spanwise stations defined in Fig. 10b. These results show a satisfactory prediction of the main flow features that develop around the VFE-2 configuration. However, in our computations, the leading-edge separation, forming downstream the outer primary vortex, is strongly delayed compared to the experimental data, and there is no evidence of an inner vortex structure, [26]. At least partially, this solution behaviour is related to the incompressible regime of our simulations. It is in fact known that the compressibility effect promotes leadingedge separation, see [27]. Moreover, in a further investigation of

the delta wing flow, which will be the subject of a forthcoming paper, we have clearly verified that the fully turbulent approach has a strong effect on the results, i.e. the solution quality significantly improves when taking into account the turbulent transition in the region close to the wing apex. Tables 7 and 8 compare our results to those obtained by Bassi et al. [28] using a compressible k–x DG code sharing almost all the numerical aspects with this incompressible solver. The results are in good agreement, especially those concerning the lift coefficient. 4.4. Flow around the DLR-F6 wing body transport configuration (Re = 106, a = 1°) The last case presented in this paper is the DLR-F6 wing body configuration without fairing, see Fig. 12. DLR-F6 is a simplified wing and fuselage geometry typical of current subsonic commer-

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

Fig. 10. VFE-2, Re = 3  106, a = 13°. P5 solution.

Fig. 11. VFE-2, Re = 3  106, a = 13°. Numerical and experimental pressure distributions.

131

132

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133

Table 7 VFE-2, Re = 3  106, a = 13°. Drag coefficient. Cd

P1

P2

P3

P4

P5

SA model k–x model

0.10230 0.10363

0.10009 0.09077

0.08487 0.09026

0.07803 0.08970

0.07467 –

Table 8 VFE-2, Re = 3  106, a = 13°. Lift coefficient. Cl

P1

P2

P3

P4

P5

SA model k–x model

0.46873 0.48091

0.49091 0.47301

0.48689 0.48026

0.47598 0.47962

0.46851 –

cial aircraft; it was developed as part of an AIAA push to evaluate CFD codes for aircraft drag prediction. The DLR-F6 case was introduced and extensively studied with numerical and experimental techniques in [29,30] and even high-order DG simulations were already performed in [17,31]. The original case is characterised by a transonic flow with shock waves on the wings. Nevertheless, here we have used an incom-

pressible solver in order to demonstrate the suitability of the method in complex geometry configurations. More in depth, we computed a case with Re = 106 and an angle of attack a = 1°. For this test case, in order to stress our numerical method, we set a ? 1 in Eq. (12), thus during the pseudo-time steps solution strategy the equations are handled as purely incompressible. The solutions were computed, up to P4 polynomial approximation, on a grid consisting of 50,618 hexahedral elements characterised, at this flow regime, by a maximum yþ c first cell dimensionless height value on the order of 1. The P1 and P2 computations were run using 28 CPU cores, the P3 solution was computed using 128 CPU cores, and the P4 solution, represented in Fig. 13, was computed using 256 CPU cores. It is worth noting that also in this case the convergence up to the machine precision can be reached only using the Bassi negative m~ modification coupled with the new r closure function. Moreover for this very large problem it is interesting to observe in Fig. 12b the low number of Newton steps required to reach full convergence. This behaviour puts in evidence the efficiency of the solver developed in this paper. In Table 9, the computed force coefficients are reported: the lift coefficient, as the polynomial order increases,

Fig. 12. DLR-F6, Re = 106, a = 1°.

Fig. 13. DLR-F6, Re = 106, a = 1°. P4 solution.

A. Crivellini et al. / Computers & Fluids 81 (2013) 122–133 Table 9 DLR-F6, Re = 106, a = 1°. Lift and drag coefficients.

Cd Cl

P1

P2

P3

P4

0.03005 0.36129

0.02423 0.44927

0.02889 0.44371

0.02641 0.44294

clearly tend to converge, while the drag coefficient is still slightly oscillating. However, its value is one order of magnitude smaller than that of the lift. Both values are compatible with reference data, even if here the flow regime is clearly completely different. 5. Conclusions The RANS-SA equations implementation proposed in [4] has been successfully validated for three-dimensional problems characterised by large separations or complex geometries. In this paper we have verified that the use of the modified r closure function introduced in [4] greatly improves the reliability and the stability of our high-order solver. A second key ingredient of the proposed implementation refers to the SA equation modifica~ case. Two different aptions developed to handle the negative m proaches have been investigated, and we saw that the simpler Bassi ones provide better performance even if in this case the ~ ¼ 0. This behaviour Jacobian matrix is no longer continuous at m ~ < 0, the Bassi modificais probably related to the fact that when m tions reduce the SA model to a simple advection–diffusion equation, while the Allmaras approach still involves stiffer non-linear source terms. Future research will be devoted to applying the solver to Detached Eddy Simulations. Acknowledgements We acknowledge the CINECA Award N. HP10C4OSL5 YEAR 2011 under the ISCRA initiative for the availability of high performance computing resources and support. References [1] Bassi F, Crivellini A, Di Pietro DA, Rebay S. An artificial compressibility flux for the discontinuous Galerkin solution of the incompressible Navier–Stokes equations. J Comput Phys 2006;218(2):794–815. [2] Bassi F, Crivellini A. A high-order discontinuous Galerkin method for natural convection problems. In: Wesseling, P, Oñate E, Périaux J, editors. Electronic proceedings of the ECCOMAS CFD 2006 conference, Egmond aan Zee, The Netherlands; September 5–8 2006. TU Delft. [3] Bassi F, Crivellini A, Di Pietro DA, Rebay S. An implicit high-order discontinuous Galerkin method for steady and unsteady incompressible flows. Comput Fluids 2007;36(10):1529–46. Special Issue Dedicated to Professor Michele Napolitano on the Occasion of his 60th Birthday. [4] Crivellini A, D’Alessandro V, Bassi F. A Spalart–Allmaras turbulence model implementation in a discontinuous Galerkin solver for incompressible flows. J Comput Phys 2013;241:388–415. [5] Wagner S, Lubon C, Kessler M. Turbulence modeling and Detached Eddy simulation with a high-order unstructured discontinuous Galerkin code. Notes Numer Fluid Mech Multidiscip Des 2010;112:143–50. [6] Bassi F, Crivellini A, Rebay S, Savini M. Discontinuous Galerkin solution of the Reynolds averaged Navier–Stokes and k–x turbulence model equations. Comput Fluids 2005(34):507–40. [7] Persson PO, Nguyen N, Peraire J. RANS solutions using high order discontinuous Galerkin methods. In: 45th AIAA aerospace sciences meeting and exhibit, Reno, Nevada; January 8–11 2007 [AIAA-2007-914].

133

[8] Oliver TA, Darmofal DL. An unsteady adaptation algorithm for discontinuous Galerkin discretizations of the RANS equations. In: 18th AIAA computational fluid dynamics conference, Miami, Florida; June 25–28 2007 [AIAA-20073940]. [9] Oliver TA, Darmofal DL. Impact of turbulence model irregularity on high-order discretization. In: 47th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition, Orlando, Florida; January 5–8 2009 [AIAA-2009-953]. [10] Kessler M, Wagner S, Landmann B, Krämer E. A parallel, high-order discontinuous Galerkin code for laminar and turbulent flows. Comput Fluids 2008;37(4):427–38. Turbulent Flow and Noise Generation. [11] Hartmann R, Held J, Leicth T, Prill F. Discontinuous Galerkin methods for computational aerodynamics – 3D adaptive flow simulation with the DLR PADGE code. Aerosp Sci Technol 2010;14(7):512–9. [12] Oliver TA. A high-order, adaptive, discontinuous Galerkin finite element method for the Reynolds-averaged Navier–Stokes equations, PhD thesis, M.I.T.; 2008. [13] Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic flows. La Rech Aérospatiale 1994;1:5–21. [14] Arnold DN, Brezzi F, Cockburn B, Marini D. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J Numer Anal 2002;39(5):1749–79. [15] Cools R. An encyclopædia of cubature formulas. J Complex 2003;19:445–53. [16] Balay S, Buschelman K, Gropp WD, Kaushik D, Knepley MG, McInnes LC, et al. PETSc Web page; 2001. . [17] Crivellini A, Bassi F. An implicit matrix-free discontinuous Galerkin solver for viscous and turbulent aerodynamic simulations. Comput Fluids 2011;50(1):81–93. [18] Karypis G, Kumar V. METIS, a software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices. Technical Report Version 4.0, University of Minnesota, Department of Computer Science/Army HPC Research Center; 1998. [19] Rumsey CL, Smith BR, Huang GP. Description of a website resource for turbulence modeling verification and validation. In: 40th AIAA computational fluid dynamics conference and exhibit, Chicago, IL; June 28–July 1 2010 [AIAA2010-4742]. [20] Rumsey CL. Consistency, verification, and validation of turbulence models for Reynolds-averaged Navier–Stokes applications. In: 3rd European conference for aerospace sciences; 2009. [21] Achenbach E. Experiments on the flow past spheres at very high Reynolds numbers. J Fluid Mech 1972;54(565). [22] Costantinescu G, Squires K. Numerical investigation of flow over a sphere in the subcritical and supercritical regimes. Phys Fluids 2004;16(5):1449–65. [23] Chu J, Luckring JM. Experimental surface pressure data obtained on 65° delta wing across Reynolds and Mach number ranges, Technical Report 4645, NASA; February 1996. [24] Konrath R, Klein C, Engler RH, Otter D. Analysis of PSP results obtained for the VFE-2 65° delta wing configuration at sub- and transonic speed. AIAA-Paper, (2006-0060); 2006. [25] Crippa S. Advances in vortical flow prediction methods for design of deltawinged aircraft. PhD thesis, KTH, Royal Institute of Technology, School of Engineering Sciences (SCI), Aeronautical and Vehicle Engineering, Stocholm, Sweden; 2008. [26] Luckring JM, Hummel D. What was learned from the new VFE-2 experiments? In: 46th AIAA aerospace sciences meeting and exhibit, Reno, NV, United States; January 7–10 2008 [AIAA]. [27] Luckring JM. ICAS 2004, Reynolds number, compressibilty, and leading-edge blutness effect on delta-wing aerodynamics. In: 24th International congress of the aeronautical sciences, Yokohama, Japan; August 29–September 3 2004 [Japan Society for the Aeronautical and Space Sciences]. [28] Bassi F, Botti L, Colombo A. VFE-2 Delta wing, medium-radius LE, subsonic case. CFD test report U1-b/UNIBG, project IDIHOM: industrialisation of highorder method – a top – down approach; February 2012. [29] Laflin K, Klausmeyer SM, Zickuhr T, Vassberg JC, Wahls RA, Morrison JH, et al. Data summary from second AIAA Computational Fluid Dynamics drag prediction workshop. J Aircraft 2005;42(5):1165–78. [30] Vassberg JC, Tinoco EN, Mani M, Brodersen OP, Eisfeld B, Morrison JH, et al. Summary of the third AIAA CFD drag prediction workshop. AIAA Paper 20070260; January 2007. [31] Bassi F, Botti L, Colombo A, Crivellini A, Franchina N, Ghidoni A, et al. Very high-order accurate discontinuous Galerkin computation of transonic turbulent flows on aeronautical configurations. In: Kroll N, Bieler H, Deconinck H, Couaillier V, Ven H, Sorensen K, editors. ADIGMA – a European initiative on the development of adaptive higher-order variational methods for aerospace applications. Notes on numerical fluid mechanics and multidisciplinary design, vol. 113. Springer-Verlag; 2010. p. 25–38.