Large Eddy Simulation of a Muffler with the High-Order Spectral Difference Method

Large eddy simulation of a muffler with the high-order spectral difference method arXiv:1305.6681v1 [physics.flu-dyn] 29 May 2013 Matteo Parsani∗† Michael Bilka‡ Chris Lacor§ May 30, 2013 Abstract The combination of the high-order accurate spectral difference discretization on unstructured grids with subgrid-scale modelling is investigated for large eddy simulation of a muffler at Re = 4.64 · 104 and low Mach number. The subgrid-scale stress tensor is modelled by the wall-adapting local eddy-viscosity model with a cut-off length which is a decreasing function of the order of accuracy of the scheme. Numerical results indicate that although the highorder solver without subgrid-scale modelling is already able to capture well the features of the flow, the coupling with the wall-adapting local eddy-viscosity model improves the quality of the solution. 1 Introduction Throughout the past two decades, the development of high-order accurate spatial discretization has been one of the major fields of research in computational fluid dynamics (CFD), computational aeroacoustics (CAA), computational electromagnetism (CEM) and in general computational physics characterized by linear and nonlinear wave propagation phenomena. High-order accurate discretizations have the potential to improve the computational efficiency required to achieve a desired error level. In fact, compared with low order schemes, high order methods offer better wave propagation properties and increased accuracy for a comparable number of degrees of freedom (DOFs). Therefore, it may be advantageous to use such schemes for problems that require very low numerical dissipation and small error levels [1]. Moreover, since computational science is increasingly used as an industrial design and analysis tool, high accuracy must be achieved on unstructured grids which are required for efficient meshing. These needs have been the driving force for the development of a variety of higher order schemes for unstructured ∗ Current Institution: Computational Aerosciences Branch, NASA Langley Research Center, Hampton, VA 23681, USA (matteo.parsani@nasa.gov) † Division of Computer, Electrical and Mathematical Sciences & Engineering, King Abdullah University of Science and Technology, Thuwal, 23955-6900, KSA ‡ University of Notre Dame, Department of Aerospace and Mechanical Engineering, 365 Fitzpatrick Hall, Notre Dame, IN 46556-5637, USA (michael.bilka.1@nd.edu) § Department of Mechanical Engineering, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium (chris.lacor@vub.ac.be) 1 meshes such as the Discontinuous Galerkin (DG) method [2, 3], the Spectral Volume (SV) method [4], the Spectral Difference (SD) method [5, 6], the Energy Stable Flux Reconstruction [7] and many others. In this study we focus on a SD solver for unstructured hexahedral grids (tensorial cells). The SD method has been proposed as an alternative high order collocation-based method using local interpolation of the strong form of the equations. Therefore, the SD scheme has an important advantage over classical 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. Although the formulation of high-order accurate spatial discretization is now fairly mature, their application for the simulation of general turbulent flows implies that particular attention has still to be paid to subgrid-scale (SGS) models. So far, the combination of the SD method with SGS models for LES has not been widely investigated. In 2010, Parsani et al. [8] reported the first implementation in study of a two-dimensional (2D) third-order accurate SD solver coupled with the Wall-Adapting Local Eddy-viscosity (WALE) model [14] and a cut-off length which is a decreasing function of the order of accuracy. A successful extension of that approach to a three-dimensional (3D) second-order accurate SD solver has been reported in [12]. Very recently, Lodato and Jameson [13] have presented an alternative technique to model the unresolved scales in the flow field: A structural SGS approach with the WALE Similarity Mixed model (WSM), where constrained explicit filtering represents a key element to approximate subgrid-scale interactions. The performance of such an algorithm has been also satisfactory. In this study, we couple for the first time the approach proposed in [8] with a 3D fourthorder accurate SD solver, for the simulation of the turbulent flow in an industrial-type muffler at Re = 4.64 · 104 . The goal is to investigate if the coupling of a high-order SD scheme with a sub-grid closure model improves the quality of the results when the grid resolution is relatively low. The latter requirement is often desirable when a high-order accurate spatial discretization is used. 2 Physical model and numerical algorithm In this study the system of the Navier-Stokes equations for a compressible flow are discretized in space using the SD method and the subgrid-scale stress tensor is modelled by the WALE approach. 2.1 Filtered Navier-Stokes equations 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 |~ u| 2 1 P velocity vector field by E = γ−1 ρ + 2 , where γ is the constant ratio of specific heats and it is 1.4 for air in standard conditions. The system, written in divergence form and equipped with suitable initial-boundary condi- 2 tions, is   ∂w ∂w ~ ~ ~ ~ · ~f = 0, + ∇ · fC (w) − ~fD w, ∇w = +∇ (1) ∂t ∂t  T ˜, ρẼ where w = ρ, ρ~u is the vector of the filtered conservative variables and ~fC = ~fC (w)   ~ represent the convective and the diffusive fluxes, respectively. Here the and ~fD = ~fD w, ∇w symbols (·) and (˜·) represent the spatially filtered field and the Favre filtered field defined as ˜ = ρ~u/ρ. ~u In a general 3D (dim = 3) Cartesian space, ~x = [x1 , x2 , x3 ]T , the components of the flux ~ vector ~f w, ∇w = [f1 , f2 , f3 ]T are given by  ρũ1 sgs ρũ21 + P − σ̃11 + τ11 sgs ρũ1 ũ2 − σ̃21 + τ21 sgs ρũ1 ũ3 − σ̃31 + τ31    sgs sgs sgs ∂ T̃ ũ1 ρẼ + P − ũ1 (σ̃11 − τ11 + q1sgs ) − ũ2 (σ̃21 − τ21 ) − ũ3 (σ̃31 − τ31 ) − cP Pµr ∂x 1    ,     sgs sgs sgs ∂ T̃ ) − ũ2 (σ̃22 − τ22 ) − ũ3 (σ̃32 − τ32 ) − cP Pµr ∂x + q2sgs ũ2 ρẼ + P − ũ1 (σ̃12 − τ12 2    ,      f3 =      sgs sgs sgs ∂ T̃ + q3sgs ũ3 ρẼ + P − ũ1 (σ̃13 − τ13 ) − ũ2 (σ̃23 − τ23 ) − ũ3 (σ̃33 − τ33 ) − cP Pµr ∂x 3    ,   2.1.1 The wall-adapted local eddy-viscosity closure model    f1 =        f2 =     ρũ2 sgs ρũ1 ũ2 − σ̃12 + τ12 sgs 2 ρũ2 + P − σ̃22 + τ22 sgs ρũ2 ũ3 − σ̃32 + τ32 ρũ3 sgs ρũ1 ũ3 − σ̃13 + τ13 sgs ρũ2 ũ3 − σ̃23 + τ23 sgs ρũ23 + P − σ̃33 + τ33   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 [15]. Both momentum and 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 [15]. The interactions of τijsgs and qisgs with the resolved scales have to be modeled through a subgrid-scale closure model because they cannot be determined using only the resolved flow field w. The smallest scales present in the flow field and their interaction with the resolved scales have to be modeled through the subgrid-scale term τijsgs . The most common approach to model such 3 a tensor is based on the eddy-viscosity concept in which one assumes that the residual stress is proportional to a measure of the filtered local strain rate [15], which is defined as follows:   δij sgs sgs S̃kk . (2) τij − τkk δij = −2 ρ νt S̃ij − 3 In the WALE model, it is assumed that the eddy-viscosity νt is proportional to the square of the length scale of the cut-off length (or width of the grid filter) and the filtered local rate of strain. Although the model was originally developed for incompressible flows, it can also be used for variable density flows by giving the formulation as follows νt = (C∆)2 S̃ . (3) Here S̃ is defined as i3/2 d S̃ d S̃ij ij S̃ = h i5/2 h i5/4 , d d S̃ij S̃ij + S̃ij S̃ij h (4) d is the traceless symmetric part of the square of the resolved velocity gradient tensor where S̃ij ∂ ũi . Note that in Equation (3) ∆, i.e., the cut-off length, is an unknown function. g̃ij = ∂x j Often the cut-off length 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 2.2, where the SD method is discussed. 2.2 Spectral difference method Consider a problem governed by a general system of conservation laws given by Equation (1) and valid on a domain Ω ⊂ Rdim with boundary ∂Ω and completed with consistent initial and boundary conditions. The domain is divided into N non-overlapping cells, with cell index i. In order to achieve an efficient implementation of the SD method, all hexahedral cells in the physical domain are mapped into cubic elements using local coordinates ξ~ = [ξ1 , ξ2 , ξ3 ]T . ~ ~ Such a transformation is characterized by the Jacobian matrix ~J i with determinant det( ~J i ). Therefore, system (1) can be written in the mapped coordinate system as ~ ~ ~ ~ ξ ξ ξ ~ ∂f1,i ∂f2,i ∂f3,i ∂wiξ ~ ξ~ · ~fiξ , =− − − = −∇ ∂t ∂ξ1 ∂ξ2 ∂ξ3 (5) ~ ~ ~ ξ~ are the conserved variables and the generalized divergence where wiξ ≡ det( ~J i ) w and ∇ differential operator in the mapped coordinate system, respectively. For a (p + 1)-th-order accurate dim-dimensional scheme, N s solution collocation points with index j are introduced at positions ξ~js in each cell i, with N s given by N s = (p + 1)dim . Given the values at these points, a polynomial approximation of degree p of the solution in cell i can 4 be constructed. This polynomial is called  the solution polynomial and is usually composed of a set of Lagrangian basis polynomial Lsj ξ~ of degree p: Ns     X Wi,j Lsj ξ~ . Wi ξ~ = (6) j=1   Therefore, the unknowns of the SD method are the interpolation coefficients Wi,j = Wi ξ~js which are the approximated values of the conserved variables wi at the solution points. ~ ~ ξ~ ·~f ξ at the solution points is computed by introducing The divergence of the mapped fluxes ∇ a set of N f flux collocation points with index l and at positions ξ~f , supporting a polynomial ~ l ξ of degree p + 1. The evolution of the mapped flux vector ~f in cell i is then approximated by ~ ~ ξ , which is obtained by reconstructing the solution variables at the flux a flux polynomial F i ξ~ ~ at these points. The flux is also represented by a Lagrange points and evaluating the fluxes F i,l polynomial: Nf     X ξ~ ξ~ ~ ~ i,l ~ F Lfl ξ~ , (7) Fi ξ = l=1 where the coefficients of the flux interpolation are defined as    ~  F ~ iξ ξ~f , ξ~lf ∈ Ωi , ξ~ l ~ = F   i,l ξ~  ~ num F ξ~lf ∈ ∂Ωi . ξ~lf , (8) ξ~ ~ Here F num is the numerical flux vector at the cell interface. In fact, the solution at a face is in general not continuous and requires the solution of a Riemann problem to maintain conservation ~ ~ ξ · ~n ξ~ must be continuous between at a cell level (i.e., the flux component normal to a face F num two neighboring cells). Approximate Riemann solvers are typically used (e.g. Rusanov Riemann ~ ~ ξ is usually taken from the interior cell. solver). The tangential component of F num ~ ~ ξ~ · F ~ ξ in the solution points results in the Taking the divergence of the flux polynomial ∇ i following modified form of (5), describing the evolution of the conservative variables at the solution points: ~ dWi,j ~ ξ~ · F ~ ξ = Ri,j , ~ ·F ~i = − 1 ∇ =−∇ (9) i dt Ji,j j j ~ i is the flux polynomial vector in the physical space whereas Ri,j is the SD residual where F associated with Wi,j . This is a system of ODEs, in time, for the unknowns Wi,j . In this work, the optimized explicit eighteen-stages fourth-order Runge-Kutta schemes presented in [16] is used to solve such a system at each time step. 5 2.2.1 Solution and flux points distributions In 2007, Huynh [9] showed that for quadrilateral and hexahedral cells, tensor product flux point distributions based on a one-dimensional flux point distribution consisting of the end points and the Legendre-Gauss quadrature points lead to stable schemes for arbitrary order of accuracy. In 2008, Van den Abeele et al. [10] showed an interesting property of the SD method, namely that it is independent of the positions of its solution points in most general circumstances. This property implies an important improvement in efficiency, since the solution points can be placed at flux point positions and thus a significant number of solution reconstructions can be avoided. Recently, this property has been proved by Jameson [11]. 2.2.2 Cut-off length ∆ In Section 2.1.1 we have seen that in the WALE model the cut-off length ∆ is used to compute the turbulent eddy-viscosity νt , i.e., Equation (3). Following the approach presented in [8], for each cell with index i and each flux points with index l and positions ξ fl , we use the following definition of filter width !#1/dim  "  det(Ji,l ) 1/dim 1 ~~ = . (10) det J i ∆i,l = Ns Ns ξf l Notice that 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. Moreover, for a given mesh, the number of solution points depends on the order of the SD scheme, so that the grid filter width decreases by increasing the polynomial order of the approximation. 3 Numerical results The main purpose of this section is to evaluate the accuracy and the reliability of the fourthorder SD-LES solver for simulating a 3D turbulent flow in an industrial-type muffler. The results are compared with the particle image velocimetry (PIV) measurement performed at the Department of Environmental and Applied Fluid Dynamics of the von Karman Institute for Fluid Dynamics [17]. In Figure 1, the geometry of the muffler and its characteristic dimensions are illustrated, where the flow is from left to right. At the inlet, mass density and velocity profiles are imposed. The inlet velocity profile in the x3 direction is given by     r 1 1 d/2 u3 = umax − tanh 2.2 − . 2 2 d/2 r At the outlet only the pressure is prescribed. In accordance to the experiments, the inlet Mach number and the Reynolds number, based on maximum velocity at the inlet umax and the diameter of the inlet/outlet d (d = 4 cm), are set respectively to Minlet = 0.05 and Re = 4.64·104 . 6 0.625 d χ3 0.625 d d χ2 5d 7.5 d 5d 4.75 d Figure 1: Configuration of the 3D muffler test case. The flow is computed using fourth-order (p = 3) SD scheme on a grid with 36, 612 hexahedral elements which was generated with the open source software Gmsh [18]. Second-order boundary elements are used to approximate the curved geometry. The total number of DOFs is approximately 2.3 · 106 (i.e., 36, 612 · (p + 1)3 ). The maximum CFL number used for the computations started from 0.1 and increased up to 0.65. After the flow field was fully developed, time averaging is performed for a period corresponding to about 25 flow-through times. The computation is validated on the center plane of the expansion coinciding with the center planes of the inlet and outlet pipes using the PIV results from [17]. All of the measurements are taken on the symmetrical center plane of the muffler. The reference cross section corresponds to the entrance of the expansion chamber. It should be noted that the circular nature of the geometry acts as a lens causing a change in magnification in the radial direction (x2 ) which prevents from capturing images close to the wall. It is found that outside 1 cm from the wall the magnification effect is negligible and as the mean stream-wise direction is in the direction of constant magnification and has only little effect on the particle correlations no corrections are deemed necessary. In Figure 2, the non-dimensional mean velocity profile in the axial direction hũ3 i/umax is shown for four different cross sections in the expansion chamber, where the PIV measurements were done. In this figure, the PIV data are also plotted for comparison. Figure 3 shows the nondimensional Reynolds stress hu′2 u′3 i/u2max at the same cross sections. Although the high-order implicit LES is already able to capture well the features of the flow field, the use of the WALE model improves the results. In particular, when the SGS model is active, the local extrema of the time-averaged velocity profiles and the second-order statistical moment (which get fairly oscillatory by moving far away from the inlet pipe) are better captured. 4 Conclusions The fourth-order SD method in combination with the WALE model and the variable filter width performed well. The numerical results confirm that the model is correctly accounting for the unresolved shear stress computed from the resolved field, for the present internal flow. However, it should be noted that the SD discretization without subgrid-scale modelling also worked rather well, at least for the grid resolution used in this study. 7 (a) 1d downstream. (b) 4d downstream. (c) 6d downstream. (d) 7d downstream. Figure 2: Time-averaged velocity profile in the axial direction hũ3 i/umax at four cross sections in the expansion chamber, obtained with fourth-order (p = 3) SD-LES method. Comparison with experimental measurements (PIV) [17]. Work is currently under way to test both approaches for different orders of accuracy, grid resolutions and other realistic turbulent flows. We believe that the flexibility of the high-order SD scheme on unstructured grids together with the development of robust sub-grid closure models for highly separated flows and efficient grid generators for high-order accurate schemes will allow to perform LES of industrial-type flows in the near future. Acknowledgement The authors would like to thank Professor David I. Ketcheson and Professor Mark H. Carpenter for their support. This research used partially the resources of the KAUST Supercomputing Laboratory and was supported in part by an appointment to the NASA Postdoctoral Program at Langley Research Center, administered by Oak Ridge Associates Universities. These supports 8 (a) 1d downstream. (b) 4d downstream. (c) 6d downstream. (d) 7d downstream. Figure 3: Reynolds stress hu′2 u′3 i/u2max in the axial direction at four cross sections in the expansion chamber, obtained with fourth-order (p = 3) SD-LES method. Comparison with experimental measurements (PIV) [17]. are gratefully acknowledged. References [1] Wang, Z. J.: Adaptive High-order Methods in Computational Fluid Dynamics (Advances in Computational Fluid Dynamics), World Scientific Publishing Company (2011). [2] Busch, K., König, M., and Niegemann, J.: Discontinuous Galerkin methods in nanophotonics. Laser Photonics Rev. 5(6), 773–809 (2011). [3] Hesthaven, Jan S., and Warburton, Tim: Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Texts in Applied Mathematics, 54. Computational Science & Engineering, Springer Publishing Company, Incorporated (2007). 9 [4] Sun, Y., Wang, Z. J., and Liu Y.: Spectral (finite) volume method for conservation laws on unstructured grids VI: extension to viscous flow. J. of Comput. Phys., 215(1), 41–58 (2006). [5] May, G., and Jameson, A.: A spectral difference method for the Euler and Navier-Stokes equations on unstructured meshes. AIAA paper 2006-304. 44th AIAA Aerospace Sciences Meeting, Reno, Nevada, USA, January 9-12, 2006. [6] Sun, Y., Wang, Z. J., and Liu, Y.: High-order multidomain spectral difference method for the Navier-Stokes equations on unstructured hexahedral grids. Commun. Comput. Phys, 2(2), 310–333 (2007). [7] Castonguay, P., Vincent P., and Jameson, A.: Application of High-Order Energy Stable Flux Reconstruction Schemes to the Euler Equations. AIAA paper 2011-686. 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Orlando, Florida, USA, January 4-7, 2011. [8] Parsani, M., Ghorbaniasl, G., Lacor C., and Turkel, E.: An implicit high-order spectral difference approach for large eddy simulation. J. Comput. Phys., 229(14), 5373–5393 (2010). [9] Huynh, H. T.: A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods. AIAA paper 2007-4079. 18th AIAA Computational Fluid Dynamics Conference, Miami, Florida, USA, June 25-28, 2007. [10] Van den Abeele, K., Lacor, C., and Wang, Z. J.: On the stability and accuracy of the spectral difference method, J. Sci. Comput., 37(2), 162–188 (2008). [11] Jameson, A.: A proof of the stability of the spectral difference method for all orders of accuracy, J. Sci. Comput., 45(1-3), 348–358 (2010). [12] Parsani, M., Ghorbaniasl, G., and Lacor C.: Validation and application of an high-order spectral difference method for flow induced noise simulation. J. Comput. Acoust., 19(3), 241–268 (2011). [13] Lodato, G., and Jameson, A.: LES modeling with high-order flux reconstruction and spectral difference schemes. ICCFD paper 2201. 7th ICCFD Conference, Big Island, Hawaii, July 9-13, 2012. [14] Nicoud, F., and Ducros, F.: Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow Turbul. Combust., 62(3), 183–200 (1999). [15] Pope, Stephen B.: Turbulent flows. Cambridge University Press (2003). [16] Parsani, M., Ketcheson, David I., and Deconinck, W.: Optimized explicit Runge-Kutta schemes for the spectral difference method applied to wave propagation problems, SIAM J. Sci. Comput, 35(2):A957-A986 (2013). 10 [17] Bilka, M., and Anthoine, J.: Experimental investigation of flow noise in a circular expansion using PIV and acoustic measurements. AIAA paper 2008-2952. 14th AIAA/CEAS Aeroacoustics Conference, Vancouver, British Columbia, Canada, May 5-6, 2008. [18] Geuzaine, C., and Remacle J.-F.: Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Meth. Eng., 79(11), 1309– 1331 (2009). 11 View publication stats