high-order simulation of zero mass jets for ... - ECCM ECFD 2018

6th European Conference on Computational Mechanics (ECCM 6) 7th European Conference on Computational Fluid Dynamics (ECFD 7) 11-15 June 2018, Glasgow, UK

HIGH-ORDER SIMULATION OF ZERO MASS JETS FOR REDUCTION OF TURBULENT SKIN FRICTION Feng Xie1,2 , Jose Daniel P´ erez Mu˜ noz3 , Ning Qin4 1

2

Department of Mechanical Engineering, The University of Sheffield Mappin Street, Sheffield S1 3JD, UK [email protected]

Hypervelocity Aerodynamics Institute, China Aerodynamics Research and Development Center Mianyang, Sichuan, 621000, China 3


4


Key words: High order, Direct numerical simulation, Zero mass jets, Flow control, Turbulent drag reduction Abstract. A 5th order MUSCL scheme is implemented in an inhouse flow solver (SHEFFlow) to investigate a unique active flow control method, using spanwise zero mass jets to manipulate the near-wall turbulence structure. Friction drag reduction has been achieved for a fully developed turbulent channel flow at Reb = 2800, based on the bulk velocity and the half channel height. Direct numerical simulation (DNS) is employed to reveal the fundamental mechanisms of the drag reduction by accurately capturing the near-wall region. Statistical analysis has been used to process the data for the understanding of the flow physics. Implementing the 5th order MUSCL scheme in SHEFFlow has been numerically tested against the Taylor-Green vortex case and the classic turbulent channel flow case. Some examples of these tests include energy dissipation rate, turbulent profile comparisons, and probability density functions. After the verification of the high order method, SHEFFlow is used to simulate the flow control method by means of DNS. Drag reduction is achieved at a jet inclination angle between 60 ∼ 85◦ . It is observed that with a proper jet inclination angle, the friction drag reduction comes from the fact that the mean velocity profile is dramatically twisted by the spanwise zero-mass-jets. Applying the 5th order MUSCL scheme in SHEFFlow helps us to investigate the near-wall region by accurately capturing the turbulent structures.

Feng Xie, Jose Daniel Pérez Mu˜ noz, Ning Qin

1

INTRODUCTION

In the DNS simulations, it is common to see the usage of high order methods due to the advantages of high order accuracy and low dissipation/dispersion. With the same grid size, the typical second order schemes have less potential than the high order methods to fully resolve the critical flow structures. In the industry, most of the commercial CFD softwares are based on the finite volume method and unstructured grids, such as Fluent, CFD++, and CFX. However, since the requirement of high accuracy simulations is increasing, the development of high order methods is more and more fast, such as discontinuous Galerkin (DG), spectral volume (SV), k-exact/WENO/MUSCL finite volume and the correction procedure via reconstruction (CPR). Most of them are very hard to be implemented by simply changing the second order code except the MUSCL scheme. With DNS simulations, one of the most practical researches is turbulent friction drag reduction which can be applied in the design of modern vehicles. For example, the turbulent drag reduction is very important for commercial airplanes. If the friction drag can be reduced by 10% for commercial airplanes, an estimate of the general saving is 250 million of dollars per year [1]. Over decades, both passive and active methods have been developed to achieve the drag reduction, such as spanwise oscillating walls [2], rotating disks [3], plasma actuators [4], traveling waves [5], opposition control [6], the uniform wall blowing/suction [7], riblets [8] and wavy surfaces [9]. The last two are passive methods and the others are active methods. The active methods are preferred, which are more controllable and have higher drag reduction than the passive methods. In this paper, the ambition is to use the 5th order MUSCL scheme to investigate a unique flow control method, spanwise zero-mass-jets, which can reduce the turbulent friction drag in the channel. The spanwise zero-mass-jets are used for the first time to reduce turbulent friction drag, borrowing some ideas from previous active methods. The mechanisms of the drag reduction will be studied by analysing the flow fields. As the controlled region is in the viscous sublayer of the velocity profile, the flow structures cannot be well resolved and accurately captured by employing any turbulent models or subgrid models, so it is important to use direct numerical simulations in this study. 2

NUMERICAL METHODS

In this study, an inhouse flow solver (SHEFFlow) is used to do the investigation, which solves the 3D compressible Navier-Stokes equations with unstructured methods. The employed numerical methods are described below. The viscous flow of a Newtonian fluid is governed by the Navier-Stokes equations. Z Z Z Z Z ∂ W dV + [F − G] · dA = 0 (1) ∂t where, V is an arbitrary control volume, dA is the differential surface area, W is the variables vector, F and G are the inviscid and viscous fluxes, respectively. The variables

2


vector and the fluxes are:   ρ         ρvx   ρvy W = ,     ρvz       ρE

 ρv      ρvvx + pî F = ρvvy + pˆ j   ˆ  ρvvz + pk   ρvE + pv

          

,

G=

          

0 τ xi τ yi τ zi τ ij vj − q

          

where, ρ, v, E, and p are the density, velocity, total energy per unit mass, and pressure 2 of the fluid, respectively. The term τij = µ vi,j + vj,i − 3 vk,k δij is the viscous stress ˆ is the position vector. tensor, q = −κ ∂T is the heat flux vector, and x = xî + yˆ j + zk ∂x Total energy is related to the total enthalpy H by E = H − p/ρ, where H = h + |v|2 /2 and h = Cp T . The ideal gas equation, P = ρRT , plays the role of a constraint for the compressible Navier-Stokes equations. The Sutherland’s law is used to calculate the viscosity and thermal conductivity of the fluid. The viscous fluxes, G, are obtained by the second order central difference, whereas the inviscid fluxes, F , are using Roe scheme [10]. The left and right quantities (QL , QR ) for the flux F i are computed by the following 5th order MUSCL scheme without a limiter function [11]: QLi+1/2 =(2Qi−2 − 13Qi−1 + 47Qi + 27Qi+1 − 3Qi+2 )/60 QR i−1/2 =(−3Qi−2 + 27Qi−1 + 47Qi − 13Qi+1 + 2Qi+2 )/60

(2)

where Q = [p, vx , vy , vz , T ]T is the quantities vector. The limiter function is not used because it would introduce extra dissipation and damage the turbulent structure. The high order method is proved to be feasible to do DNS simulations [12]. Since SHEFFlow is a compressible flow solver, a preconditioning method is employed to compute the low-speed simulation. It is done by using a dual-time-step method proposed by [13], which has the following form: Z Z Z ∂ ∂ W dΩ + Γ QdΩ + (F − G)ndS = 0 (3) ∂t Ω ∂τ Ω S where Γ is the preconditioning matrix, t and τ are pseudo time and physical time, respectively. The time marching scheme is 3-stage second order Runge-Kutta. 3 3.1

VALIDATION Taylor-Green vortex

To test the accuracy and the performance of the 5th order MUSCL method in the DNS simulation, the DNS of Taylor-Green vertex at Re = 1600 is carried out. The geometry of this case is a periodic cubic with −πL x, y, z πL and it has 8 vortices inside. The mesh is uniform in every direction and periodic boundary conditions are applied in the

3


0.160

5123 Spectral 2563 3rd order Present 2563 5th order Present

0.140

5123 Spectral 2563 3rd order Present 2563 5th order Present

0.016

0.120 0.012

Ek

0.100 0.080

0.008

0.060 0.004 0.040 0.020

0.000 0

5

10

15

20

0

5

10

t

t

(a)

(b)

15

20

Figure 1: (a) and (b) are the comparisons of the history of non-dimensional kinetic energy Ek and energy dissipation rate between SHEFFlow and a spectral code on a mesh of 5123 grid points [14]. (t is non-dimensionalized by tc = L/V0 .) simulation. The initial flow field can be described by the following equations: x y z u =V0 sin cos cos , L L L y z x sin cos , v = − V0 cos L L L w =0, 2x 2y 2z ρ0 V02 cos + cos cos +2 . p =p0 + 16 L L L

(4)

where, p0 = 105527.625Pa, ρ0 = 1.225kg m−3 , V0 = 24.11m s−1 , L = 0.001m. The heat capacity ratio γ = 1.4 and Prandtl number Pr = 0.71. In addition, the value of p0 which is the Reynolds number is defined as Re = ρ0 Vµ0 L and the temperature T0 = Rρ 0 uniform for the initial flow field. The simulation duration is from t = 0 to 20tc , where tc = VL0 is the characteristic convective time. Both 3rd and 5th order of the MUSCL method are tested with 2563 grids and compared with a spectral code which is using 5123 grids, showing a tiny difference between the 5th order MUSCL method and the method [14, 15] on the evolution R Spectral 1 v·v of non-dimensional kinetic energy Ek = ρ0 Ω Ω ρ 2 dΩ and kinetic energy dissipation rate k = − dE in figure 1(a)(b). When the kinetic energy dissipation rate is high, the 3rd order dt MUSCL method has bigger than the other two before the maximum of the dissipation, meanwhile, the maximum dissipation happens earlier than the other two, which causes the kinetic energy of the 3rd order MUSCL method is lower than that of the other two. But after the maximum dissipation, the energy dissipation rate of the 3rd order MUSCL 4

Feng Xie, Jose Daniel Pérez Mu˜ noz, Ning Qin y x

ly

z lz Mean streamwise flow lx

Figure 2: Schematic of the turbulent channel flow. The dimensions of the computational domain are lx , ly and lz in the streamwise, wall-normal and spanwise directions, respectively. method is smaller than that of the other two. After t = 12.5, the 3rd order MUSCL method catches up with the other two, showing the three methods have similar kinetic energy. In the end, an homogeneous isotropic turbulence is obtained from the simulation of the Taylor-Green vortex which creates small scales and follows a decay phase. The validation of the Taylor-Green vortex proves that the 5th order MUSCL method is able to produce a reasonable evolution of the kinetic energy and the kinetic energy dissipation rate which are important to precisely simulate a turbulent flow. Moreover, the 5th order MUSCL method is better than the 3rd order, so it is necessary to use 5th order for the later simulations rather than using low order methods. 3.2

Fully developed turbulent channel flow

The fully developed turbulent channel flow shown in figure 2 is a classic case to validate the DNS ability of SHEFFlow. The results can be validated by the incompressible DNS data from reference [16] and their online database. For the boundary conditions, top and bottom are isothermal and no-slip wall conditions. Periodic conditions are applied on the streamwise and spanwise directions. The flow field is initialised by a parabolic velocity profile on which a random perturbation is applied. Due to the friction along with the walls and very small numerical dissipation, the mass flux rate would decrease without pressure difference in the streamwise direction. To keep the mass flow rate, a background force is introduced through the pressure derivative with respect to x, which is ∂p/∂x. This method was derived by Xu et al. [17]. According to the method, the pressure gradient can be split into two parts: ∂p/∂x = ∂pmean /∂x + ∂pf /∂x, where ∂pmean /∂x is the driving pressure gradient and ∂pf /∂x is the fluctuation of pressure gradient derived from the Navier-Stokes equations. In reference [17], a variable which is defined as = ∂pmean /∂x is employed to keep the mass flow rate in the simulation. is calculated as: 1 (m ˙ 0 − 2m ˙ n+m ˙ n−1 ) (5) n+1 = n − ∆τ where m ˙ 0 is the initial mass flux rate of the channel, m ˙ n is the mass flux rate at physical time step n, and ∆τ is the time of one physical time step. 5


20

4 3.5

DNS by et Kimal. & Moin Kim 1987aa DG-DES

Present

15

o

DNS by Kim & Moin DNSKim by Kim Moin1987 et& al. DNS by Kim & Moin DG-DES DG-DES Present DG-DES

3

o

+

UU

+

2.5 10

2 1.5

5

1 0.5

0

0 10

0

10 y + +

1

10

2

0

y

50

100 y+ +

150

y

(b) Root-mean-square velocity fluctuations normalized by the friction velocity: u0rms /uτ , 0 0 , ; vrms /uτ , , ; wrms /uτ , , .

(a) The mean streamwise velocity profiles.

Figure 3: Comparisons between SHEFFlow and Kim et al. [16] In this study, both dimensional and non-dimensional variables are used. The dimensional variables do not have any superscript. If the variables employ superscript +, it means the variables are non-dimensionalised in wall-units, which are defined as velocity V + = V /uτ , length L+ = ρuτ L/µ and time T + = ρu2τ T /µ. The validation is performed at Reτ = 180 based on the friction velocity uτ and the half channel height δ, using the physical time step ∆τ + = 0.25. The computational domain has lx /δ = 4π and lz /δ = 4π/3 for the streamwise length and spanwise width, respectively. The mean velocity is 42.0 m s−1 and the mean pressure is 16065.05 Pa. The temperature is 110.0 K and the half channel height is 0.001m. Moreover, the mesh resolution are ∆x+ = 8.84 and ∆z + = 4.49 in the streamwise and spanwise directions, respectively. ∆yw+ = 0.2 is the first wall-normal grid spacing from the walls and ∆yc+ = 4.13 is the wall-normal grid spacing at the centre of the channel. The mesh has 256 × 197 × 168 in the streamwise, wall-normal and spanwise directions respectively. + In the non-dimensional duration Ttotal = 2250, the time-averaged Cf is 8.20 × 10−3 compared to 8.18 × 10−3 from the DNS simulation of Kim et al. [16], which has 0.24% difference. Figure 3 shows a good agreement between the statistic results of SHEFFlow and the reference [16], including the mean streamwise velocity profiles and the root-meansquare velocity fluctuations. In this validation, by comparing the DNS results between SHEFFlow and Kim et al. [16] at Reb = 2800, it shows that SHEFFlow can accurately capture the turbulent structures and have the ability to investigate the drag reduction of skin friction by using zero-mass-jets in the channel flow.

6


4 4.1

SIMULATIONS OF THE CHANNEL WITH SPANWISE ZERO-MASSJETS Computational realization and flow conditions

In the previous sections, by using the high order method, SHEFFlow is able to do the DNS of turbulent channel flow which is one of the most classic cases to study the friction force and to investigate the unique flow control method of using spanwise zeromass-jets. The channel which has the spanwise zero-mass-jets will be called controlled channel to distinguish the smooth channel which is shown in subsection 3.2. The details of the computational model are shown in figure 4. To generate the spanwise zero-massjets, 12 steps are evenly put on the walls, which have 12 pairs of zero-mass-jets to obtain + oscillating layers in the near wall region. The height of the steps is δjet and the width + of the steps is Lstep . The zero-mass-jets come out from the sides of the steps. For a pair of zero-mass-jets, when one is blowing out, the other one is sucking in, showing on the left of figure 4. The jets have a distance L+ jets between a pair of zero-mass-jets and have an angle β which is against the incoming flow on the x-z plane. When β = 0◦ , the jets are fully along with the spanwise direction if there is no incoming flow. Because the Reynolds number of the zero-mass-jets is low, the zero-mass-jets can be assumed as laminar flow and the velocity profiles of the zero-mass-jets are parabolic. The velocity of the zero-mass-jets is in a sinusoidal function, 2π + + + t (6) Ujet = Wmax sin + Tosc + + are the velocity and the maximum velocity at the centre of the where, Ujet and Wmax + jet exits, respectively. Tosc and t+ are the period and the physical time of the zero-massjets, respectively. All of the variables are in wall units. According to the previous works of spanwise oscillating walls [2], a saving in power can reach 39% with the maximum + + wall velocity Wwall,max = 18 and the oscillating period Twall,osc = 125. In this study, the spanwise zero-mass-jets are used to emulate the similar mechanisms to the spanwise oscillating wall with the benefit of being more practical. Therefore, the period of the + zero-mass-jets is set to be Tosc = 125 and the maximum of the mean velocity of the + zero-mass-jets is set to be Wmean,max = 18. Then, the maximum velocity at the centre of + = 27 because the velocity distribution is parabolic at the jet exits. the jet exits is Wmax A pair of zero-mass-jets has two sides, called side.1 and side.2 in this study. At side.1, + + + + = Ujet cos β, the streamwise and spanwise velocities are Uside.1,x = −Ujet sin β and Uside.1,z + + respectively. At side.2, the streamwise and spanwise velocities are Uside.2,x = Ujet sin β + + and Uside.2,z = Ujet cos β, respectively. The simulation parameters of the controlled channel are displayed in Table 1. To make the controlled channel be able to compare with the smooth channel which does not have steps and jet exits, the computational domain of the controlled channel has the same streamwise and spanwise dimensions as the smooth channel. But, because of the steps on the walls, the controlled channel has two heights, either of which is not able to be the same with the height of the smooth channel. If one of the two heights is the same

7


y x Side.1 + δjet

z ly

L+ step lx

β + Ujet

Side.2

L+ jet

lz Mean streamwise flow

+ + Ujet = Wmax sin

2π + + t Tosc

+ Figure 4: The schematic of the controlled channel with spanwise zero-mass-jets. L+ x , Ly and L+ z are the streamwise, wall-normal and spanwise dimensions of the computational domain, respectively. The steps length and the height of the jet exits are in wall units based on the friction velocity in the smooth channel at Reb = 2800.

Mean Mach number M ab Mean velocity Ub Mean pressure Pb Mean temperature Tb Mean density ρb Reference length δ Nx × Ny × Nz Sx × Sy × Sz lx /δ × ly /δ × lz /δ + + δstep × L+ step × Ljet ∆x+ × ∆(yw+ ∼ yc+ ) × ∆z + + + ∆x+ s × ∆ys × ∆zs Number of devices

0.211 42.0 m s−1 16065.1 Pa 110.0 K 0.509 kg m−3 0.001m 256 × 198 × 168 256 × 10 × 20 4π × 1.98 × 34 π 2 × 35.9 × 89.8 8.84 × (0.2 ∼ 4.09) × 4.49 8.84 × 0.2 × 4.49 12

Table 1: Simulation parameters of the fully developed turbulent channel flow with steps. ∆x and ∆z are the collocation resolutions parallel to the wall, using Nx and Nz cells. ∆yw is the first wall-normal cell spacing from the wall and ∆yc is the wall-normal cell spacing at the centre of the channel. Ny is the number of cells in the wall-normal direction. In between a pair of two jets, Sx , Sy and Sz are the number of cells in streamwise, normal and spanwise direction, respectively. The friction velocity in the smooth channel at Reb = 2800 is used to convert the variables in wall units.

8


with the height of the smooth channel, the other height will be different and the area of the streamwise cross-section will be not the same with the area of the streamwise cross-section of the smooth channel. Then, because both the controlled channel and the smooth channel are set to have the same streamwise mass-flux, the mean velocities will be different between the controlled channel and the smooth channel in the streamwise direction. Consequently, Reb will be different between the controlled channel and the smooth channel. Therefore, the two heights of the controlled channel are adapted to ensure that both the controlled channel and the smooth channel have the same area of the streamwise cross-section. Since the areas of the streamwise cross-sections are the same, the mean velocity can be the same in the streamwise direction because the mass fluxes of the controlled channel and the smooth channel are fixed to be a same constant in the streamwise direction, employing the method introduced in Subsection 3.2. The mesh resolution of the controlled channel is similar to that of the smooth channel to make sure that the investigation of the unique flow control method is under DNS. 4.2

Simulation results

In the turbulent channel flow, for an instantaneous result, the drag coefficient never is constant at a single point on the wall, even the drag coefficient is space averaged in the limited area of the walls. The history of the space-averaged drag coefficients is shown in Figure 5b. Therefore, the drag coefficients are also time averaged in a sufficient time to consider the overall effects, which are statistically convergent. To make sure that the added steps in the channel only have little influence on the drag, a baseline model is simulated with the zero-mass-jets off. The drag coefficient of the baseline model is 8.22 × 10−3 from t+ = 3750 to 7343 which is long enough to do 32.7 sweeps over the whole channel. The result of the baseline model is 0.49% different from that of the smooth channel, which indicates that the influence of the steps is small and the change of the drag coefficient is mainly because of the zero-mass-jets. The drag reduction is shown in Figure 5a, indicating that the zero-mass-jets can obtain positive drag reduction between β = 60◦ and β = 85◦ . The value of gross drag reduction (DR) compares the drag coefficients between the controlled channel and the smooth channel, which is defined as DR(%) = 100

Cf,smooth − Cf,controlled Cf,smooth

(7)

When the value of the drag reduction is negative, it means that the drag is increased. Moreover, from β = 70◦ to β = 80◦ , the values of the drag reduction are similar with a maximum value of 19.7% at β = 75◦ . However, the drag is increased 104.0% at β = 0◦ , which gives the maximum drag coefficient in the results. When the jet angle is between 0◦ and 60◦ , the drag is almost linearly reduced, meanwhile, when β = 90◦ , the zero-mass-jets have little influence as they do not inject momentum in the flow. Figure 5 suggests that a remarkable drag reduction can be obtained with a proper jet angle and the amount of the drag reduction can be controlled by changing the jet angle. Moreover, for the baseline model with the jets off and the controlled channel with β = 75◦ , 9


24

Baseline model (Jets off) Jets on with β = 0◦ Jets on with β = 75◦

20

Cf × 103

DR (%)

16

12

8

4

0 0

β(◦ )

200

400

600

800

1000

1200

t+ = ρu2τ t/µ

(a)

(b)

Figure 5: (a) Correlation between β and the level of gross drag reduction. (b) The history of drag coefficients in the different conditions of the zero-mass-jets.

u0 /uτ

(a) Baseline model (jets off)

(b) Jets on with β = 75

0 Figure 6: Iso-surface of λ+ 2 = −2 coloured by u /uτ .

10

◦


the λ+ 2 iso-surface of the instantaneous flow fields is shown in Figure 6, which can visually reveal the cores of the vortices. Obviously, the baseline model has fewer vortices than the controlled channel with β = 75◦ , which means the zero-mass-jets introduces more eddies into the channel flow. 5

CONCLUSIONS

In this paper, it proves that the 5th order MUSCL scheme is successfully implemented in SHEFFlow to do the direct numerical simulations and the unique flow control method is investigated to reduce the turbulent skin friction by the spanwise zero-mass-jets. By using direct numerical simulations, the viscous layer is accurately simulated without any turbulent models. A summary of the preceding investigation is listed below. • The 5th order MUSCL scheme has good performance of the kinetic energy evolution and the kinetic energy dissipation rate by simulating the Taylor-Green vortex. Based on the same mesh, the 5th order MUSCL scheme is more accurate than the 3rd order MUSCL scheme, showing a good ability on the decaying simulation of the homogeneous isotropic turbulence. • By using the 5th order MUSCL scheme, SHEFFlow simulated the fully developed turbulent channel flow at Reb = 2800. The result of the turbulent skin friction is only 0.24% different from the standard. Moreover, the turbulent structures can be clearly and accurately captured, showing good agreement by comparing the mean streamwise velocity profiles, the root-mean-square velocity fluctuations, the shear stresses and the normalised PDFs of three directions with the standard. • A unique flow control method is introduced, employing the zero-mass-jets in the spanwise direction. Comparing to the smooth channel, the turbulent skin friction can be reduced when the jets angle is between 60◦ and 85◦ , however, it is increased when the jets angle is between 0◦ and 60◦ . Moreover, the biggest drag reduction happens at 75◦ of the jets angle. References [1] Michael J Walsh, William L Sellers III, and Catherine B Mcginley. Riblet drag at flight conditions. Journal of Aircraft, 26(6):570–575, 1989. [2] Maurizio Quadrio and Pierre Ricco. Critical assessment of turbulent drag reduction through spanwise wall oscillations. Journal of Fluid Mechanics, 521:251–271, 2004. [3] Pierre Ricco and Stanislav Hahn. Turbulent drag reduction through rotating discs. Journal of Fluid Mechanics, 722:267–290, 2013. [4] O Mahfoze and S Laizet. Skin-friction drag reduction in a channel flow with streamwise-aligned plasma actuators. International Journal of Heat and Fluid Flow, 66:83–94, 2017. 11


[5] Davide Gatti and Maurizio Quadrio. Reynolds-number dependence of turbulent skinfriction drag reduction induced by spanwise forcing. Journal of Fluid Mechanics, 802:553–582, 2016. [6] Bing-Qing Deng, Chun-Xiao Xu, Wei-Xi Huang, and Gui-Xiang Cui. Strengthened opposition control for skin-friction reduction in wall-bounded turbulent flows. Journal of Turbulence, 15(2):122–143, 2014. ¨ u, and Philipp Schlatter. Effect of [7] Yukinori Kametani, Koji Fukagata, Ramis Orl¨ uniform blowing/suction in a turbulent boundary layer at moderate reynolds number. International Journal of Heat and Fluid Flow, 55:132–142, 2015. [8] M Walsh. Turbulent boundary layer drag reduction using riblets. In 20th aerospace sciences meeting, page 169, 1982. [9] Sacha Ghebali, Sergei I Chernyshenko, and Michael A Leschziner. Turbulent skinfriction reduction by wavy surfaces. arXiv preprint arXiv:1705.01989, 2017. [10] Philip L Roe. Approximate riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 43(2):357–372, 1981. [11] Kyu Hong Kim and Chongam Kim. Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows: Part i: Spatial discretization. Journal of Computational Physics, 208(2):527–569, 2005. [12] Chung-Gang Li, Makoto Tsubokura, and Keiji Onishi. Feasibility investigation of compressible direct numerical simulation with a preconditioning method at extremely low mach numbers. International Journal of Computational Fluid Dynamics, 28(610):411–419, 2014. [13] Jonathan M Weiss and Wayne A Smith. Preconditioning applied to variable and constant density flows. AIAA journal, 33(11):2050–2057, 1995. [14] Zhijian Wang, Krzysztof Fidkowski, Rémi Abgrall, Francesco Bassi, Doru Caraeni, Andrew Cary, Herman Deconinck, Ralf Hartmann, Koen Hillewaert, Hung T Huynh, et al. High-order cfd methods: current status and perspective. International Journal for Numerical Methods in Fluids, 72(8):811–845, 2013. [15] Xiaodong Liu, Lijun Xuan, Yidong Xia, and Hong Luo. A reconstructed discontinuous galerkin method for the compressible navier-stokes equations on threedimensional hybrid grids. Computers & Fluids, 152:217–230, 2017. [16] John Kim, Parviz Moin, and Robert Moser. Turbulence statistics in fully developed channel flow at low reynolds number. Journal of Fluid Mechanics, 177:133–166, 1987. [17] Xiaofeng Xu, Joon Sang Lee, and Richard H Pletcher. A compressible finite volume formulation for large eddy simulation of turbulent pipe flows at low mach number in cartesian coordinates. Journal of Computational Physics, 203(1):22–48, 2005. 12