Evaluation of Riemann flux solvers for WENO reconstruction schemes: Kelvin–Helmholtz instability Omer San a,⇑ , Kursat Kara b a Interdisciplinary Center for Applied Mathematics, Virginia Tech, Blacksburg, VA, USA b Department of Aerospace Engineering, Khalifa University, Abu Dhabi, United Arab Emirates article info Article history: Received 22 April 2014 Received in revised form 17 April 2015 Accepted 30 April 2015 Available online 14 May 2015 Keywords: Euler equations WENO reconstruction schemes Approximate Riemann solvers Numerical dissipation Kelvin–Helmholtz instability Two-dimensional turbulence abstract Accurate and computationally efficient simulations of Euler equations are of paramount importance in both fundamental research and engineering applications. In this study, our main objective is to investi- gate the efficacy and accuracy of several Riemann solvers for high-order accurate weighted essentially non-oscillatory (WENO) reconstruction scheme as a state-of-the-art tool to study shear driven turbulence flows. The Kelvin–Helmholtz instability occurs when a perturbation is introduced to a continuous fluid system with a velocity shear, or where there is a velocity difference across the interface between two flu- ids. Here, we solve a stratified Kelvin–Helmholtz instability problem to demonstrate the performance of six different Riemann solvers’ ability to evolve a linear perturbation into a transition to nonlinear hydro- dynamic two-dimensional turbulence. A single mode perturbation is used for our evaluations. Time evo- lution process shows that the vortices formed from the turbulence slowly merge together since both energy and enstrophy are simultaneously conserved in two-dimensional turbulence. Third-, fifth- and seventh-order WENO reconstruction schemes are investigated along with the Roe, Rusanov, HLL, FORCE, AUSM, and Marquina Riemann flux solvers at the cell interfaces resulting in 18 joint flow solvers. Based on the numerical assessments of these solvers on various grid resolutions, it is found that the dis- sipative character of the Riemann solver has significant effect on eddy resolving properties and turbu- lence statistics. We further show that the order of the reconstruction scheme becomes increasingly important for coarsening the mesh. We illustrate that higher-order schemes become more effective in terms of the tradeoff between the accuracy and efficiency. We also demonstrate that AUSM solver pro- vides the least amount of numerical dissipation, yet resulting in a pile-up phenomenon in energy spectra for underresolved simulations. However, results obtained by the Roe solver agree well with the theoret- ical energy spectrum scaling providing a marginal dissipation without showing any pile-up at a cost of around 30% increase in computational time. Ó 2015 Elsevier Ltd. All rights reserved. 1. Introduction Accurate and efficient solutions of Euler equations have been subject to intensive research for at least five decades and many successful numerical methods have been proposed for solving them (e.g., see [1–12]). Euler equations are a system of non-linear hyperbolic conservation laws that govern the dynamics of compressible material such as gases or liquids at high pressures, for which the effects of body forces, viscous stresses and heat flux are neglected [13]. They are very important for many areas includ- ing astrophysics, weather and aerospace simulations. Astrophysical gas flows are often highly supersonic and require correct treatment of shocks and other discontinuities (e.g., contact discontinuities, across which the density and temperature jump but the pressure, and tangential discontinuities, across which the tangential velocity changes) [14]. For aeronautical applications, accurate computations of shock and vortex dominated flows are important for aerodynamics shape optimization and load calcula- tions during conceptual and preliminary design phases [15,16]. Traditionally, second-order accurate numerical methods are often preferred in solutions of Euler equations because of their simplicity and robustness. The main deficiency of these methods for accurate computation of vortex dominated flows is the numer- ical diffusion of vorticity to unacceptable levels [17]. Application of high-order methods can significantly alleviate this deficiency and be more efficient for the problems requiring high accuracy, such as wave propagation problems, vortex-dominated flows, large http://dx.doi.org/10.1016/j.compfluid.2015.04.026 0045-7930/Ó 2015 Elsevier Ltd. All rights reserved. ⇑ Corresponding author. E-mail address: omersan@vt.edu (O. San). Computers & Fluids 117 (2015) 24–41 Contents lists available at ScienceDirect Computers & Fluids journal homepage: www.elsevier.com/locate/compfluid