On the coupling between a high-order spectral difference method and large eddy simulation

V European Conference on Computational Fluid Dynamics ECCOMAS CFD 2010 J. C. F. Pereira and A. Sequeira (Eds) Lisbon, Portugal,14-17 June 2010 ON THE COUPLING BETWEEN A HIGH-ORDER SPECTRAL DIFFERENCE METHOD AND LARGE EDDY SIMULATION Matteo Parsani∗ , Ghader Ghorbaniasl∗ , Michael Bilka† and Chris Lacor∗ ∗ Vrije Universiteit Brussel, Department of Mechanical Engineering, Fluid Dynamics and Thermodynamics Research Group, Pleinlaan 2, 1050 Brussel, Belgium e-mails: mparsani@vub.ac.be, ghader.ghorbaniasl@vub.ac.be, chris.lacor@vub.ac.be † von Karman Institute for Fluid Dynamics, Department of Environmental and Applied Fluid Dynamics Chaussée de Waterloo 72, 1640 Rhode St.-Genèse, Belgium e-mail: bilka@vki.ac.be Key words: High-order Spectral Difference Method, Large Eddy Simulation, WallAdapting Local Eddy-Viscosity Model, Implicit LU-SGS algorithm. Abstract. The filtered fluid dynamic equations are discretized in space by a high-order spectral difference(SD) method coupled with large eddy simulation (LES) approach. The subgrid-scale stress tensor is modelled by the wall-adapting local eddy-viscosity model (WALE). We solve the unsteady equations by advancing in time using a second-order backward difference formula (BDF2). The nonlinear algebraic system arising from the time discretization is solved with the nonlinear lower-upper symmetric Gauss-Seidel (LU-SGS) algorithm. As test cases, the 2D turbulent flow past a square cylinder at Re = 2.2 × 104 and the 3D turbulent flow in a muffler at Re = 4.64 × 104 have been computed and the results compared with available solutions and obtained experimental data. A good agreement between the results is observed. 1 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor 1 INTRODUCTION Spatially high-order accurate compact schemes are now well established as accurate methods to solve a variety of flow problems. In Computational Fluid Dynamics (CFD), they are used for Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), Computational Aeroacoustics (CAA) etc., where the accurate resolution of small scales is required. Moreover, since CFD is increasingly used as an industrial design and analysis tool, high-order accuracy must be achieved on unstructured grids which are required for efficient meshing. These needs have been the driving force for the development of higher order schemes for unstructured meshes such as the Discontinuous Galerkin (DG) method,1, 2 the Spectral Volume (SV) method3, 4 and the Spectral Difference (SD) method.5, 6 All these methods use piecewise continuous functions as the solution approximation space. They are capable of achieving high-order accuracy on unstructured grids and they have a compact stencil, which makes them easily parallelizable. The SD method has an important advantage over the DG and SV methods, that no integrals have to be evaluated to compute the residuals, thus avoiding the need for costly high-order accurate quadrature formulas. Recently, there has been research on unifying several of the popular methods including the DG method, the SV method and the SD method with a technique that does not require the evaluation of the integral.7, 8 Although spatially high-order accurate numerical schemes guarantee the accurate resolution of small scales, their application for the simulation of general turbulent flows implies that particular attention has still to be paid to subgrid models. In this framework, we couple a high-order SD scheme on unstructured quadrilateral grids with the local eddy-viscosity (WALE) model9 to carry out large eddy simulations. The accuracy and the reliability of the of the new coupling between the SD method and the WALE model are tested by solving the 2D ”turbulent” flow around a square cylinder at Re = 2.2 × 104 with a third-order SD scheme and the 3D turbulent flow in a muffler at Re = 4.64 × 104 with a second-order scheme. The results are compared with LES and experiments reported in literature. Usually, explicit temporal discretizations such as multi-stage TVD (total variation diminishing) RungeKutta schemes1 are used to advance the solution in time. In general, explicit schemes and their boundary conditions are easy to implement, vectorize and parallelize, and require only limited memory storage. However, for large-scale simulations and especially for high-order solutions, the rate of convergence slows down dramatically, resulting in inefficient solution techniques. In addition, the solver must also be able to deal with the geometrical stiffness imposed by the Navier-Stokes grids where high-aspect ratios occur near walls. In the case of compressible solvers there is an additional stiffness when solving for low speed flows caused by the disparate eigenvalues of the system. Therefore, to exploit the potential of high-order methods efficient solvers are needed. In order to speed up convergence, a multigrid strategy and/or an implicit temporal discretization is required. Implicit schemes can advance the solution with significantly larger time 2 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor steps compared to explicit methods. However, they may be more expensive than explicit schemes if the algebraic solver employed is not efficient.12, 24 In the present study, the unsteady filtered fluid dynamic equations are solved by advancing in time using a secondorder backward difference formula (BDF2). The nonlinear algebraic systems arising from the time discretization is then solved with the nonlinear lower-upper symmetric GaussSeidel (LU-SGS) algorithm.10–13 The remainder of this article is organized as follows. A brief summary of the filtered fluid dynamic equation for a compressible flow and the description of the WALE model are given in Section 2. Section 3 is devoted to the description of the SD method and the treatment of the diffusive terms the second approach of Bassi and Rebay.14 In the same section the coupling of the SD method and the WALE model through a new definition of the grid filter width is given. The nonlinear LU-SGS algorithm combined with the BDF2 is described in Section 4. Section 5 deals with the numerical results, before finally drawing the conclusions in Section 6. 2 GOVERNING EQUATIONS FOR LARGE EDDY SIMULATION In this section the system of the LES fluid dynamic equations for a compressible flow are presented. The three physical conservation laws for a general Newtonian fluid, i.e. the continuity, the momentum and energy equations, are introduced using the following notation: ρ for the mass density, ~u ∈ Rdim for the velocity vector in a physical space with dim dimensions, P for the static pressure and E for the specific total energy which is related to the pressure and the velocity vector field by E= 1 P |~u|2 + , γ−1 ρ 2 where γ is the constant ratio of specific heats and it is 1.4 for air. Define a vector w of all the filtered conservative variables, i.e.   ρ w =  ρ~u˜  , ρẼ (1) (2) where the symbols (·) and (˜·) represent respectively the spatially filtered field and the Favre filtered field defined as ~u˜ = ρ~u/ρ. The system of the filtered fluid dynamic equations for a compressible flow, written in divergence form and equipped with suitable initialboundary conditions, is   ∂w ∂w ~ ~ ~ ~ · ~f = 0, ~ + ∇ · fC (w) − fD w, ∇w = +∇ (3) ∂t ∂t   ~fC (w) = [fC , g , hC ]T and ~fD w, ∇w ~ = [fD , gD , hD ]T represent the convective and the C diffusive fluxes, respectively. In a general 3D (dim = 3) Cartesian space, ~x = [x1 , x2 , x3 ]T , 3 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor the components of these fluxes are given by    ρũ1 ρũ2  ρũ2 + P   ρũ1 ũ2 1    2    ρũ1 ũ2 fC =   , gC =  ρũ2 + P      ρũ1 ũ3   ρũ2 ũ3  ũ1 ρẼ + P ũ2 ρẼ + P    fD =       gD =       hD =             , hC =      0 sgs σ̃11 + τ̃11 sgs σ̃21 + τ̃21 sgs σ̃31 + τ̃31 ∂ T̃ σ̃11 ũ1 + σ̃21 ũ2 + σ̃31 ũ3 + cP Pµr ∂x + q1sgs 1 0 sgs σ̃12 + τ̃12 sgs σ̃22 + τ̃22 sgs σ̃32 + τ̃32 ∂ T̃ σ̃12 ũ1 + σ̃22 ũ2 + σ̃32 ũ3 + cP Pµr ∂x + q2sgs 2 0 sgs σ̃13 + τ̃13 sgs σ̃23 + τ̃23 sgs σ̃33 + τ̃33 ∂ T̃ σ̃13 ũ1 + σ̃23 ũ2 + σ̃33 ũ3 + cP Pµr ∂x + q3sgs 3 ρũ3 ρũ1 ũ3 ρũ2 ũ3 2 ρũ  3+P  ũ3 ρẼ + P      ,     ,      ,      ,   where cP , µ, P r and T represent respectively the specific heat capacity at constant pressure, the dynamic viscosity, the Prandtl number and the temperature of the fluid. Moreover, σij represents the ij−component of the resolved viscous stress tensor defined as   δij σij = 2µ Sij − Smm i, j = 1, . . . , dim, (4) 3     1 ∂ ũi ∂ ũj ˜ Sij ~u = S̃ij = i, j = 1, . . . , dim, (5) + 2 ∂xj ∂xi where δij is the Kronecker Delta function. From the definitions of the fluxes components it is seen that both the momentum and the energy equations differ from the classical fluid dynamic equations only for two terms which take into account the contributions from the unresolved scales. These contributions, represented by the specific subgrid-scale stress tensor τijsgs and by the subgrid heat-flux vector defined qisgs , appear when the spatial filter is applied to the convective terms and they are defined as follows τijsgs = ρ (ug i uj − ũi ũj ) 4 i, j = 1, . . . , dim, (6) M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor   g qisgs = cP ρ T ui − T̃ ũi i = 1, . . . , dim. (7) The interactions of τijsgs and qisgs with the resolved scales have to be modeled through a subgrid-scale model because they cannot be determined using only the resolved flow field w. 2.1 The wall-adapted local eddy-viscosity model In the previous section, it was seen that the smaller scales and their interaction with the resolved scales have to be modeled through the subgrid-scale term τijsgs . The tensor τijsgs can be modeled at different levels of complexity. The most common approach is based on the eddy-viscosity concept in which one assumes that the residual stress is proportional to the filtered rate of strain, which is defined as follows:   δij sgs sgs τij − τkk δij = −2 ρ νt S̃ij − S̃kk = −2 ρ νt S̃ijD , (8) 3 In the wall-adapted local eddy-viscosity (WALE) proposed by Nicoud and Ducros,9 it is assumed that the eddy-viscosity νt is proportional to the square of the length scale of the filter and the filtered local rate of strain. Although the model was originally been developed for incompressible flows, it can also be used for variable density flows by giving the formulation as follows νt = (C∆)2 S̃ . (9) Here S̃ is defined as where S̃ijd is given by S̃ijd = h i3/2 d d S̃ij S̃ij S̃ = h i5/2 h i5/4 , d d S̃ij S̃ij + S̃ij S̃ij
δij 2 1 2 2 g̃ij + g̃ji − g̃kk , 2 3 with g̃ij2 = (10) ∂ ũi ∂ ũk . ∂xk ∂xj (11) Note that in Eq. (9) ∆, i.e. the width of the grid filter, is an unknown function. Often the grid filter width is taken proportional to the smallest resolvable length scale of the discretization. In the present work, the definition of the grid filter function is given in Section 3.1, where the spectral difference method is discussed. The WALE model is specifically designed to return the correct wall-asymptotic y +3 −variation of the subgrid-scale viscosity νt 9 and the constant model coefficient C can be adjusted so that the correct amount of subgrid dissipation is obtained. For the subgrid heat-flux vector qisgs , if an eddy diffusivity model is used, the following expression is obtained   νt Cp ∂ T̃ sgs g , (12) qi = Cp ρ T ui − T̃ ũi = −ρ P rt ∂xi 5 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor where the value of the turbulent Prandtl number P rt is set to 0.72 and the eddy-viscosity is computed by Eq. (9). 3 SPECTRAL DIFFERENCE METHOD Consider a problem governed by the general system of conservation laws given by Eq. (3) and valid on a domain V for the filtered conservative variable vector w defined in Eq. (2). The domain is divided into N non-overlapping cells, with cell index i. On each cell, a 3D mapped coordinate system ξ~ = [ξ1 , ξ2, ξ3 ]T is introduced, with the transformation to Cartesian coordinates defined by     x1,i x1,i (ξ1 , ξ2 , ξ3 )   ~xi =  x2,i  =  x2,i (ξ1 , ξ2 , ξ3 )  = ~xi ξ~ . (13) x3,i x3,i (ξ1 , ξ2 , ξ3 ) ~ The Jacobian matrix for this transformation is denoted as ~J i and the Jacobian determi~ nant as Ji . The fluxes projected in the mapped coordinate system ~f ξ are related to the Cartesian flux components ~f for the cell i by  ξ~    f f i i ~ ~ −1 ~ −1 ~  ~f ξ =  (14)  gξi  = Ji ~J i  gi  = Ji ~J i ~fi . i ~ hi hξ i Therefore, the governing equations (3) can be written in the mapped coordinate system as ~ ~ ~ ~ ~ ∂wξi ∂fiξ ∂gξi ∂hξi ∂ (Ji w) ~ ξ~ · ~fiξ , (15) ≡ =− − − = −∇ ∂t ∂t ∂ξ1 ∂ξ2 ∂ξ3 ~ with wξi ≡ Ji w the conservative variables in the mapped coordinate system. For a dimdimensional (p + 1)-th order scheme, N s (p, dim) solution points with index  j are intros s ~ ~ duced at positions ξj , supporting a set of Lagrangian basis polynomials Lj ξ of degree p. Based on these basis polynomials, one can approximate the solution in cell i with a p-th order polynomial as follows Ns   X   wi ≈ Wi ξ~ = Wi,j Lsj ξ~ , (16) j=1 where the conservative variables at the solution points Wi,j denote the solution variables of the SD method. The evolution of these variables is governed by Eq. (15) evaluated at the solution points. ~ ~ ξ~ ·~fiξ at the solution points, a set of To estimate the divergence of the mapped fluxes ∇ N f flux points with index l and at positions ξ~f , supporting a polynomial of degree p + 1, l 6 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor ~ ξ is introduced. The evolution of the mapped flux vector ~f in cell i is then approximated ~ ~ ξ , which is obtained by reconstructing the solution variables at the by a flux polynomial F i ~ ~ ξ at these points. To ensure a coupling between flux points and evaluating the fluxes F i,l the cells, a number of flux points need to lie at the faces or the corners of the cell. In order to maintain conservation at a cell level, the flux component normal to a face must be continuous between two neighbouring cells. Since the solution at a face is in general not continuous, this requires the introduction of Riemann solvers at those points. Two different approaches were discussed in Wang et al.5 The first approach involves the definition of multidimensional Riemann solvers, while the second one uses multiple 1D Riemann solvers. One can find more information in the references mentioned above. The flux polynomial is then defined by ~  ~ ξ ξ~ F i f = N X l=1   ~ ~ ξ Lf ξ~ , F i,l l (17)   where the Lfl ξ~ are the Lagrangian basis polynomials associated with the flux points. ~ ~ ξ~ · F ~ ξ in the solution points results in the Taking the divergence of the flux polynomial ∇ i following modified form of (15), describing the evolution of the conservative variables in the solution points: dWi,j ~ ·F ~i =−∇ dt j =− 1 ~ ξ~ ~ ξ~ ∇ · Fi = Ri,j , Ji,j j (18) where Ri,j is the SD residual associated to Wi,j . This is a system of ordinary differential equations in time for the solution unknowns Wi,j , which can be solved numerically using any classical method for such system. In the present paper, only meshes with quadrilateral (2D) and hexahedral (3D) are considered. For such cells, different sets of flux points are used for the different components ~ of the mapped flux vector. For the f ξ -component, a set of flux points that supports a polynomial of degree p + 1 in ξ1 and of degree p in ξ2 and ξ3 is defined. These flux ~ points are labeled the ‘ξ1 -flux points’ in the remainder of this paper. For the gξ - and ~ hξ -components, ‘ξ2 -’ and ‘ξ3 -flux points’ are introduced analogously. The solution and flux point distributions for third-order accurate quadrilateral SD cells are illustrated in Fig. 1. In Van den Abeele et al.15 it was shown that the distribution of the solution points has very little influence on the properties of the SD schemes, and in fact, for linear problems, the different distributions in the figures lead to identical results. It is easily understood that the rightmost distributions in the figures allow for a significant reduction in the solution reconstruction cost. Notice that for the third-order scheme, there is a free 7 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor 1 1 1 0 0 0 −1 −1 −1 0 0.58 α3 (a) General. 1 −1 −1 0 0.58 α3 (b) Locally 1D. 1 −1 0 0.58 α3 1 Most solution points at flux points, symmetrical. (c) Figure 1: Third-order quadrilateral SD cells. Solution points (◦) and ξ1 - (H) and ξ2 -flux points (N). parameter α3 for the flux point distributions. In general, the flux points are not uniquely defined for schemes with p > 1. The choice of the flux point distribution was shown to have an important influence on the stability and accuracy of the SD schemes in Van den Abeele et al.15 In Huynh et al.,7 it was proven that for quadrilateral and hexahedral cells, tensor product flux point distributions based on a 1D flux point distribution consisting of the end points and the Legendre-Gauss quadrature points, lead to stable schemes for arbitrary values of p. Consequently, in the present paper, symmetrical distributions with most solution points at flux points, shown in Fig. 1c, are used, where the flux point distribution is defined using Legendre-Gauss quadrature points. Finally, since only one component of the projected flux vectors is reconstructed in each set of flux points, a traditional 1D Riemann flux, like the Rusanov or the Roe flux, suffices for the treatment of the convective fluxes at face flux points in quadrilateral and hexahedral SD cells. In the present paper the second approach of Bassi and Rebay (BR2)14 is considered for   ~ the diffusive terms treatment. For the evaluation of the diffusive flux vector ~fD w, ∇w , the gradients of the conservative variables must be available at the flux points. Defining the vectors ~Jiξ1 , ~Jiξ2 and ~Jiξ3 as  ξ1  ~ ~J , ~J ξ2 , ~J ξ3 T , Ji ~J −1 (19) i = i i i ~ i (ξ) is obtained by computing the values at the a gradient approximation polynomial Φ solution points as follows: " # ~J ξ1 ∂ W̆i ~J ξ2 ∂ W̆i ~J ξ3 1 ∂ W̆ i i i i ~ ~ i,j , ≈ ∇w + + =Φ (20) i,j Ji,j ∂ξ1 ∂ξ2 ∂ξ3 j where W̆i is a polynomial of degree p + 1, defined by its values at the flux points. In an internal flux point, this is just the value of the polynomial Wi . In a face flux point, it is 8 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor equal to an average value Ŵ, of the two available values WL and WR , which is defined for the second approach of Bassi and Rebay as Ŵ = WL + WR . 2 (21)   ~ i ξ~ , given by The gradients in cell i are then approximated by Φ Ns   X   s ~ ~ ~ ~ ~ ∇w ≈ Φi ξ = Φi,j Lj ξ . i (22) j=1   ~ ~ The diffusive flux vector is approximated as fD Wi , Φi in an internal flux point. At a ~ L and Φ ~ R , of Φ ~ are available. An averaged value Φ ~ˆ is then used and face, two values, Φ it is defined as ~ ~ ~ ~ ~ˆ = ∇WL + ∇WR + ϕ ΛL + ΛR , Φ (23) 2 2 where ϕ defines the amount of damping added to the gradients. ϕ is always equal to ~ L and Λ ~ R associated to a face can be one in the present work. The lifting operators Λ interpreted as corrections to the gradients of the solution polynomials in the neighbouring cells. These lifting operators are polynomials defined in the neighbouring cells by their values at the solution points: h i ~ L(R),j = 1 ~ ξ~ δWL(R) . ∇ (24) Λ JL(R),j j In this expression, δW is a polynomial of degree p + 1, defined by its values at the flux points:   T  (W − W ) J ~~J −1 ~1 ξ~ ~1n l ∈ curr. face R,l L,l n (25) δWL(R),l = L(R),l  0 elsewise, ~ system. The normal where ~1nξ is the unit normal to the face in the mapped coordinate   ˆ ~ ~ component of the diffusive flux vector is thus evaluated as fD 0.5 (WL + WR ) , Φ · ~1n in a face flux point, where the unit normal vector in the physical space ~1n to a face points from the left (L) to the right (R) neighbouring cell. Unlike the Local Spectral Difference approach,13, 14 the BR2 approach is fully compact, as only the immediate neighbors are required for the computation of the residuals in a cell. 9 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor 3.1 Grid filter width for the subgrid-scale model In Section 2.1 it is seen that in the WALE model the grid filter width ∆ is used to compute the turbulent eddy-viscosity, i.e. Eq. (9). In general ∆ is an unknown function and it has to be defined to have a closed model. Often, the grid filter width is taken proportional to the smallest resolvable length scale of the discretization and for a general cell with index i is usually approximated by dim Y ∆i = hk k=1 !1/dim , (26) where hk is the size of the cell in the k−direction. However, at the face flux points, two values of ∆ are available, labeled ∆L and ∆R . Consequently, an averaged value for the filter width is generally used, ˆ i = ∆L + ∆R . ∆ (27) 2 For classical finite volume (FV) methods, Eq. (27) uniquely defines the grid filter width because for these schemes the flux points always lie at the cell face. However, the same reasoning applied to a general SD method implies a natural formulation of the grid filter width based on the Jacobian determinant of the transformation defined by Eq. (13). However, since in the SD scheme each cell has interior solution points and a high-order polynomial approximation occurs in the cell, it is natural to choose the filter width depending on the order of the polynomial. Therefore, in the definition of the grid filter width, the order of the polynomial is taken into account through the division of the Jacobian determinant by the number of solution points, i.e. N s (p, dim). Consequently, for each cell with index i and each flux points with index l and positions ξ fl , we propose the following new definition of filter width  1 ∆i,l = det Ns  ~~ Ji ξfl 1/dim =  Ji,l Ns 1/dim . (28) Eq. (28) uniquely defines the filter width for the internal flux points but for the face flux points two values of the Jacobians determinant are available, labeled JL,l and JR,l . Consequently, an averaged value is again used, i.e.   det ˆ i,l =  ∆   ~~ JL ξfl  + det 2 Ns  ~~ JR ξfl  1/dim    =  JL,l + JR,l 2 Ns 1/dim (29) Notice that with this approach, the cell filter width is not constant in one cell, but it varies because the Jacobian matrix is a function of the positions of the flux points. 10 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor Moreover, for a given mesh, the number of solution points depends on the order of the SD scheme, so that grid filter width vary by varying the order of the scheme. The proposed approach is valid for high-order approaches and it is consistent because the filter width is a function of the the polynomial order through the number of solution point. In fact, the grid filter width decreases by increasing the polynomial order of the approximation. 4 NONLINEAR LOWER-UPPER SYMMETRIC GAUSS-SEIDEL METHOD Applying the SD method to a set of convection-diffusion equations results in a semidiscretization in space of the form (18). This expression can be discretized in time using any classical method. If an explicit time stepping algorithm, like a Runge-Kutta method, is used, then, because of the Courant-Friedrichs-Lewy stability condition for the convective terms, there is a maximum value, proportional to the cell size and inversely proportional to the maximum wave propagation velocity, for the associated time step ∆t. This stability limit is often too restrictive, especially for high-order schemes and Navier-Stokes grids. Implicit time-integration schemes can be used to deal with these problems because they do not suffer from such a stability limit. In the present work, the second-order backward difference formula (BDF2) with variable time-step, τn2 1 + 2τn W|tn+1 − (1 + τn ) W|tn + W|tn−1 = ∆tn R|tn+1 , 1 + τn 1 + τn (30) is employed to integrate, in time, the system of ODEs which arises from the spatial discretization. In Eq. (30) n is the index of the time iteration, ∆tn = tn+1 − tn and n . The BDF2 is A−stable.16 This property is very important to solve systems of τn = ∆t∆tn−1 stiff ODEs which often arise from the discretization of the fluid dynamic equations with a spatially high-order numerical scheme. Expression (30) is a nonlinear algebraic system, which has to be solved each time iteration to find the solution at the next iteration tn+1 starting from the two previous time iterations tn and tn−1 . Writing (30) for one cell, with the current cell denoted by the subscript cc, the neighbouring cells that contribute to its residual denoted by the subscript nb, omitting the solution point index j, and linearizing the residual about time iteration tn , one obtains ! X ∂Rcc ∆Wcc ∂Rcc − β1 ∆Wcc + ∆Wnb = ∆tn ∂Wcc tn ∂W nb tn nb  (31) β2 Wcc |tn − Wcc |tn−1 + β1 Rcc |tn , with ∆W = W|tn+1 − W|tn , β1 = 1+τn 1+2τn and β2 = 11 τn2 . ∆tn (1+2τn ) Applying a (symmetric) M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor Gauss-Seidel (SGS) algorithm to solve this linear algebraic system results in  ∂Rcc −β1 ∂Wcc tn  I ∆Wm+1 = β1 + cc ∆tn X ∂Rcc Rcc |tn + ∆W∗nb ∂W nb tn nb   + β2 Wcc |tn − Wcc |tn−1 , ! (32) where the approximation of the linear system solution at the iteration m is written as ∆Wm , and the superscript ∗ signifies the latest available solution. ∂Rcc To avoid the storage of the off-diagonal block matrices ∂W , expression (32) is further nb manipulated to obtain  ∂Rcc −β1 ∂Wcc tn tn    ∂Rcc I m+1 ∗ ∗ ∆Wcc = β1 Rcc − + ∆Wcc ∆tn ∂Wcc tn   + β2 Wcc |tn − Wcc |tn−1 . (33) Defining δWm+1 = ∆Wm+1 −∆W∗cc = Wm+1 −W∗cc , the final expression of the nonlinear cc cc cc LU-SGS algorithm for the BDF2 with variable time step is then:    ∆W∗  I ∂Rcc ∗ cc δWm+1 = β R + β W | . (34) − W | + −β1 1 cc 2 cc tn cc tn−1 − cc ∂Wcc tn ∆tn ∆tn The inverse of the small Jacobian matrices in the left hand side of this expression can be computed using a LU decomposition at the beginning of each time iteration, which makes the solution of the small linear algebraic systems much more efficient during subsequent SGS sweeps. Note that this LU-SGS algorithm acts directly on the nonlinear algebraic system to be solved, which is the right hand side of expression (34). For this reason the algorithm is called nonlinear LU-SGS algorithm. Since only the diagonal block Jacobians should be stored, the total number of real variables N needed for these Jacobians, in 3D, is  2 N (p + 1)3 × # physical variables . (35) From expression (35), it is clear that the nonlinear LU-SGS method requires significantly less memory than the classical method that use the full Jacobian matrix (for instance the GMRES algorithm17 ), but the required amount still increases with p to the power six. 5 NUMERICAL RESULTS The new coupling between the SD method and LES is first tested by solving the 2D ”turbulent” flow past a square cylinder at Re = 2.2 × 104 . The obtained results are compared with the LES and experimental results reported in Bouris et al.22 Afterwards, 12 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor the proposed numerical approach is applied to the 3D turbulent flow in a muffler at Re = 4.64 × 104 and the statistics are compared with those obtained from experiments and other numerical simulations. The nonlinear system (34) is solved with multiple cell-wise symmetric forward and backward sweeps with a prescribed tolerance of 10−6 on the change of the L2 norm of the (m+1) solution variation δWcc and/or a maximum number of 100 symmetric Gauss-Seidel sweeps. The meshes were created using Gmsh23 which allows to have a second-order (p = 2) polynomial approximation of the boundary elements. 5.1 2D flow past a square cylinder The purpose of this test case is to evaluate the quality of the proposed numerical method by comparing the mean flow field and the turbulent kinetic energy results with some reference solution. In the literature, 2D and 3D large eddy simulations of the flow past a square cylinder have been performed, for instance see Murakami et al.18 and Breuer et al.21 Moreover, Bouris et al.22 have also performed 2D LES computations and compared the results with the experimental measurements of Lyn19 and Durao et al.20 It was demonstrated that the characteristic of this type of flow is its quasi-two-dimensional character and the presence of periodic vortex shedding from the front corners of the square rod which introduces a low-frequency variation of the velocity field behind the rod in addition to the high-frequency turbulence fluctuations. It has been stated in the past22 that 2D LES calculations are clearly inferior to 3D ones since certain important features of 3D turbulence are not resolved. The three dimensionality of turbulence cannot be questioned, however, in the present paper, we want to show that the importance of detailed simulation of the quasi-two-dimensional mechanisms can be achieved by performing a 2D LES with the new combination of SD and subgrid-scale model. In Fig. 2, the configuration of the test case is illustrated. The cylinder is placed on the length axis (y = 0) of the computational domain. At the left-hand-side boundary (the inflow), the flow is prescribed to be uniform with non-fluctuating velocity profiles. The same conditions are applied to the upper boundary and to the lower boundary. At the right-hand-side boundary (the outflow), only the pressure is prescribed. The Mach number is set to 0.05, so that the flow is almost incompressible. The Reynolds number based on the free-stream velocity | ~u∞ | and height of the square cylinder h is 2.2 × 104 . A mesh with 12622 quadrilateral elements is used. The flow is computed using third-order (p = 2) SD schemes with and without subgrid scale model. The time-step used for the computation starts from 0.00001 and increases linearly up to 0.0020. In the present test case, the 2D LES calculation with a second finite volume scheme and the Smagorisky-Lilly model of Bouris et al.22 will be used as a reference 2D LES solution. In the reference solution, a mesh with 350×300 points is employed. Therefore, the total number of degrees of freedom (DOFs) of the present SD calculation is practically the same as the number of DOFs used in Bouris et al.22 Table 1 shows the prediction of the dominant vortex shedding frequency in non13 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor x2 x2 = 6 h F ar f ield h Inlet x1 x1 = −5.5 h Outlet x1 = 15.5 h F ar f ield Figure 2: Configuration of the 2D square cylinder test case. dimensional form reported in literature22 as well as the results from the present simulations. In this table, the abbreviations (k − ε)∗ , RSE1 and RSE2 stand respectively for k −ε model with 2 layers, Reynolds stress equation with wall function and Reynolds stress equation with two layers (see Bouris et al.22 and the references therein.). The results with subgrid scale model are in good agreement with the experimental measurements while the solution without subgrid-scale model underestimates both the Strouhal number and the time averaged drag coefficient hcd i. Table 1: Numerical calculation of the Strouhal number and the mean drag coefficient from various turbulence simulations of the flow past a square cylinder. St (k − ε)∗22 RSE122 RSE222 2D LES22 0.124 0.136 0.159 0.134 Present SD no model 0.121 2.15 2.43 2.18 1.98 hcd i 1.179 Present Exps.19 ,20 SD-LES 0.133 0.132± 0.004 0.139 2.21 2.05-2.23 Time averaged results are obtained by integrating the data over approximately 20 shedding cycles and the mean center-line velocity is presented in Fig. 3. Although the SD scheme without subgrid-scale model overpredicts the value of the reverse velocity, it already gives a solution which has somewhat the same accuracy as does the 2D LES of Bouris et al.22 Moreover, it can be clearly seen that the modeling of the subgrid-scale stress tensor improves the accuracy of the results. In fact, a good agreement between the predicted results and experimental data of Durao et al.19 throughout the comparison domain shows that the quality of a high-order SD method increases when it is coupled with large eddy simulation. In the other parts of the comparison domain the method without subgrid-scale model is in a good agreement with LES of Bouris et al.22 Note that, the SD-LES method captures the peak of the mean stream-wise velocity considerably better than the others. Figure 4 shows the resolved total turbulent kinetic energy. One can compare the results with the experimental data and see that the SD method with subgrid scale model undershoot the peak experimental value somewhat less than the 2D LES of Bouris et al.22 14 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor 1 SD no model SD−LES Exp. (Durao et al.) Exp. (Lyn et al.) 2D LES (Bouris et al.) 3D LES (Murakami et al.) 3D LES (Breuer et al.) 0.8 0.6 hũ1 i |~ u∞ | 0.4 0.2 0 −0.2 −0.4 −2 0 Figure 3: Time averaged velocity 2 x1 h hũ1 i |~ u∞ | 4 6 8 along the center-line behind the square cylinder. 1 SD no model SD−LES Exp. (Durao et al.) Exp. (Lyn et al.) 2D LES (Bouris et al.) 3D LES (Murakami et al.) 0.8 ktot |~ u∞ |2 0.6 0.4 0.2 0 −0.2 −0.4 −2 0 Figure 4: Total fluctuation energy 2 x1 h ktot |~ u∞ | 2 4 6 8 along the center-line behind the square cylinder. and predicts the results well elsewhere. As it is observed, the present approach without subgrid scale model overshoots the peak experimental value less than does the 2D LES of Bouris et al.,22 whereas throughout most of the domain both results are on the top of each other. 5.2 3D flow in a muffler Mufflers are commonly used in a wide variety of applications. Industrial flow ducts as well as internal combustion engines frequently make use of silencing elements to attenuate the noise levels carried by the fluids and radiated to the outside atmosphere by the exhausts. Restrictive environmental legislation requires that silencer designers use high performance and reliable techniques. Various techniques are currently available for the modeling and testing of duct mufflers. Empirical, analytical and numerical techniques have been used and proved reliable under controlled conditions. Design of a complete muffler system is, usually, a very complex task. Each element is selected by considering its particular performance, cost and its interaction effects on the overall system performance 15 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor 0.625 d 0.625 d d and reliability. The main purpose of this work is to evaluate the accuracy and the reliability of the LES solver, performing the simulation of a confined flow in an industrial geometry. Therefore, the 3D turbulent flow in a muffler is considered as a test case. The results are compared with the experimental results obtained at the Department of Environmental and Applied Fluid Dynamics of the von Karman Institute for Fluid Dynamics, Belgium. 5d 7.5 d 5d 4.75 d Figure 5: Dimensions of the 3D muffler. In Fig. 5, the geometry of the muffler its characteristic dimensions are illustrated. This configuration was selected because is representative of the physical mechanisms involved in a generic muffler. At the inlet, according to a ”top-hat” velocity profile obtained from the experiments, total pressure, total temperature and direction of the flow are specified. At the the outlet only the pressure is prescribed. The Reynolds number based on maximum velocity at the inlet | ~uinlet |= umax and the diameter of the inlet/outlet d is 4.64 × 104 . In accordance to the experiments, the Mach number is set to 0.05. A very coarse grid with 125072 hexahedral elements and with a quadratic geometrical mapping is used (see Fig. 6). The flow is computed using second-order (p = 1) SD schemes with subgrid scale model. Therefore, the number of DOFs is approximately one million. The time-step used for the computation starts from 0.00001 and increases linearly up to 0.001. (a) x1-x2 Plane (b) x2-x3 Plane. Figure 6: Plane views of the grid used for the 3D LES of the muffler.23 The computation is validated on the center plane of the expansion coinciding with the center planes of the inlet and outlet pipes by means of particle image velocimetry (PIV). The laser sheet was generated by a 200 mJ Nd-YAG double-pulsed laser passing through 16 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor a spherical and cylindrical lens. A smoke generator was used to see the flow with fine oil droplets (1 µm). A digital CCD camera was used (PCO Sensicam, resolution - 1280x1024) to record the images with a region of interest limited to 1280x832. The pulse separation of the between the images of the PIV couple correspond to a maximum displacement of 6 pixels. In order to keep a high spatial resolution over the measurement domain, the measurements are divided into 3 zones and with a field of view of 100x55mm2 with a resolution of 12 pixels per millimeter. All of these measurements are taken on the symmetrical center plane of the muffler. It should be noted that the circular nature of the geometry acts as a lens causing a change in magnification which prevents from capturing images close to the wall. It is found that outside 1cm from the wall the magnification effect is negligible and as the mean flow direction is in the direction of constant magnification (i.e along the height of the cylinder) no corrections are deemed necessary. For each region, 1,000 image pairs have been acquired and pre-processed with the home-made cross-correlation algorithm WIDIM (Window Displacement Iterative Multigrid). This program is based on an iterative multigrid predictor-corrector method, handling the window distortion, for better resolution of shear flows, and the sub-pixel window displacement, to limit pixel-locking. The predictor-corrector method is then validated for each grid size if the signal-to-noise (SN) ratio is above 1.5.25 Mean quantities are then computed from the PIV data and profiles are extracted to compare with the LES. 3i is shown In Fig. 7, the non-dimensional mean velocity profile in the axial direction uhũmax for three different sections in the expansion chamber. In this figure, the PIV data and the LES solution obtained with FluentTM are also plotted for comparison. Moreover, in Fig. hu′ u′ i 8, the non-dimensional Reynolds stress u22 3 at the same sections is shown. It can be max seen that, although the mesh is very coarse, the mean velocity profile computed with the SD-LES, is in good agreement with the experiments throughout the comparison domain and is better than the FluentTM results. However, for the Reynolds stress, in some parts of the comparison domain, a deviation from the PIV data is observed. This is due to the fact that the mesh is too coarse (see Fig. 6) to resolve very well the fluctuations in the flow, even if a second-order SD scheme with subgrid scale model is used. A future work is needed to investigate the effect of the grid size and the turbulence modeling on the quality of the Reynolds stresses when the new numerical approach is applied to the chosen muffler test case. 6 CONCLUSIONS A new numerical approach has been proposed for space discretization of the filtered Navier-Stokes equations. The approach was based on a high-order Spectral Difference scheme coupled with large eddy simulation. The accuracy and the reliability of the solver are tested by solving the 2D turbulent flow around a square cylinder at Re = 2.2 × 104 and the 3D turbulent flow in a muffler at Re = 4.64 × 104. The new SD-LES approach offers a much more realistic flow picture than the standard LES. It has been shown that the SD-LES method improves the prediction of the unsteady flow features and reaches a 17 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor 1d downstream 0.75 0.5 0.5 0.5 hu˜3 i/umax 0.75 0.25 0.25 0 0 −0.25 −0.25 −0.25 0 x2 /d 1 2 −0.5 −2 3 (a) 1 d downstream −1 hũ3 i umax 1d downstream 0.02 x2 /d 1 2 −0.5 −2 3 3d downstream 0.015 ′ 0 −0.005 −0.005 −0.01 −0.01 −0.01 −0.015 −0.015 −0.015 x2 /d 1 2 3 (a) 1 d downstream Figure 8: Reynolds stress 3 ′ ′ 0 ′ ′ hu2 u3 i/u2max 0.005 hu2 u3 i/u2max 0.005 0 2 SD−LES LES−Fluent PIV 0.015 0.005 −1 1 6d downstream 0.01 −0.005 x2 /d 0.02 SD−LES LES−Fluent PIV 0.01 0 0 (c) 6 d downstream. 0.01 −0.02 −2 −1 at three sections in the expansion chamber. 0.02 SD−LES LES−Fluent PIV 0.015 0 (b) 3 d downstream. Figure 7: Time averaged velocity profile ′ 0.25 0 −1 SD−LES LES−Fluent PIV 1 0.75 −0.5 −2 hu2 u3 i/u2max 6d downstream SD−LES LES−Fluent PIV 1 hu˜3 i/umax hu˜3 i/umax 3d downstream SD−LES LES−Fluent PIV 1 −0.02 −2 −1 0 x2 /d 1 2 (b) 3 d downstream hu′2 u′3 i u2max 3 −0.02 −2 −1 0 x2 /d 1 2 (c) 6 d downstream. at three sections in the expansion chamber. quality that is comparable to experiments. Regarding the Reynolds stresses in the muffler, for better accuracy of the results, two issues need to be addressed - (a) grid fineness and (b) turbulence modeling. The overall performances of the SD-LES method have proven to be accurate and reliable in the simulation of fluid dynamic equations. Acknowledgement This research was funded by IWT under Project SBO 050163. This funding is gratefully acknowledged. 18 3 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor REFERENCES [1] B. Cockburn, C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework. Mathematics of Computation (1989), 52:411-435. [2] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin method for timedependent convection-diffusion systems. SIAM Journal of Numerical Analysis (1998), 35(6):2440-2463. [3] Y. Sun, Z.J. Wang, Y. Liu, Spectral (finite) volume method for conservation laws on unstructured grids VI: Extension to viscous flow. Journal of Computational Physics (2006), 215:41-58. [4] K. Van den Abeele, C. Lacor, An accuracy and stability study of the 2D spectral volume method. Journal of Computational Physics (2007), 226:1007-10026. [5] Z.J. Wang, Y. Liu, G. May, A. Jameson, Spectral difference method for unstructured grids II: extension to the Euler Equations. Journal of Scientific Computing (2007), 32(1):449-454. [6] Y. Sun, Z. J. Wang, Y. Liu, High-order multidomain spectral difference method for the Navier-Stokes equations on unstructured hexahedral grids. Communications in Computational Physics (2007), 2(2):310-333. [7] H. T. Huynh, A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods. 18th AIAA Computational Fluid Dynamics Conference, AIAA 2007-4079, Miami, Florida, 25-28 June 2007. [8] Z.J. Wang, Haiyang Gao, A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids. Journal of Computational Physics (2009); 228:81618186. [9] F. Nicoud, F. Ducros, Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbulence and Combustion (1999), 62(3):183-200. [10] A. Jameson, E. Turkel, Implicit schemes and LU decompositions. Mathematics of Computations (1981), 37:385-397. [11] Y. Sun, Z.J. Wang, Y. Liu, C.-L. Chen, Efficient implicit LU-SGS algorithm for high-order spectral difference method on unstructured hexahedral grids. 45th AIAA Aerospace Sciences Meeting, AIAA 2007-313, Reno, Nevada, 8-11 January 2007. [12] M. Parsani, K. Van den Abeele, C. Lacor, E. Turkel, Implicit LU-SGS algorithm for highorder methods on unstructured grid with p-multigrid strategy for solving the steady NavierStokes equations. Journal of Computational Physics (2010), 229:828-850. 19 M. Parsani, G. Ghorbaniasl, M. Bilka and C. Lacor [13] M. Parsani, G. Ghorbaniasl, C. Lacor, E. Turkel, An Implicit High-Order Spectral Difference Approach For Large Eddy Simulation. Journal of Computational Physics (2010), Accepted for Publication on March 25 2010. [14] F. Bassi, S. Rebay, GMRES discontinuous Galerkin solution of the compressible NavierStokes equations, in: K. Cockburn, Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, Springer, Berlin, 2000, 11:197208. [15] K. Van den Abeele, C. Lacor, Z.J. Wang, On the stability and accuracy of the spectral difference method. Journal on Scientific Computing (2008), 37(2):162-188. [16] A. Quarteroni, R. Sacco and F. Saleri Numerical mathematics, Springer, Berlin, 2002. [17] Y. Saad and M. H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing (1986), 7:856-869. [18] S. Murakami, A. Mochida, On turbulent vortex shedding flow past a square cylinder predicted by CFD. Journal of Wind Engineering and Industrial Aerodynamics (1995), 54:191-211. [19] D. Durão, M. Heitor, J. Pereira, Measurements of turbulent and periodic flows around a square cross-section cylinder. Experiments in Fluids (1988), 6:298-304. [20] D. Lyn, S. Einav, W. Rodi, J. Park, A laser Doppler velocimetry study of ensembleaveraged characteristics of the turbulent near wake of a square cylinder. Journal of Fluid Mechanics (1995), 304:285-319. [21] M. Breuer, M Poriquie, First experiences with LES of flows past bluff bodies. In 3rd International Symposium on Engineering Turbulence Modelling and Measurements, Elsevier, Amsterdam 1996. [22] D. Bouris, G. Bergeles, 2D LES of vortex shedding from a square cylinder. Journal of Wind Engineering and Industrial Aerodynamics (1999), 80:31-46. [23] C. Geuzaine, J.-F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering (2009), 79(11):1309-1331. [24] R.C. Swanson, E. Turkel, C.-C. Rossow, Convergence acceleration of Runge-Kutta schemes for solving the Navier-Stokes equations. Journal of Computational Physics (2007), 224:365388. [25] Fulvio Scarano, Particle Image Velocimetry: Development and Application. PhD Thesis, Università Degli Studi di Napoli ’Federico II’ (2000). 20