A new hybrid Lattice-Boltzmann method for thermal flow simulations in low-Mach number approximation

Guanxiong WANG (王冠雄), [START_REF] Draisey | Three-dimensional natural convection in a box: a numerical study[END_REF] Song ZHAO (赵崧), [START_REF] Draisey | Three-dimensional natural convection in a box: a numerical study[END_REF] Pierre BOIVIN, [START_REF] Draisey | Three-dimensional natural convection in a box: a numerical study[END_REF] Eric SERRE, [START_REF] Draisey | Three-dimensional natural convection in a box: a numerical study[END_REF] and Pierre SAGAUT [START_REF] Draisey | Three-dimensional natural convection in a box: a numerical study[END_REF] Aix-Marseille Univ., CNRS, Centrale Marseille, M2P2, Marseille, France (*Electronic mail: guanxiong.wang.isae@hotmail.com) (Dated: 11 April 2022)

A new low-Mach algorithm for the thermal Lattice Boltzmann method (LBM) is proposed aiming at reducing the computational cost of thermal flow simulations in the low Mach number limit. The well known Low Mach Number Approximation (LMNA) is adopted to accelerate the simulations by enlarging the time-step through re-scaling the pseudo acoustic speed to the same order of the fluid motion velocity. This specific process is inspired by the similarity between the artificial compressibility method and the isothermal LBM, and is further extended to its thermal counterpart. It must be emphasized that such low-Mach acceleration strategy is in a general form, thus can be easily applied to other compressible LB methods. The present method overcomes the drawback of classical PGS (Pressure Gradient Scaling) method due to the pressure gradient changing. The new algorithm is validated by various well-documented academic test cases in laminar (1D (one dimensional) gravity column, 2D (two dimensional) rising thermal bubble, 2D differentially heated square cavity) and turbulent (3D (three dimensional) Taylor-Green vortex and 3D heated cylinder) regimes. All the results show excellent agreement with the reference data and high computational efficiency.

I. INTRODUCTION

Low-speed thermal flows (low-Mach number limit) have been of scientific interest for a long time due to their fundamental importance in the study of buoyancy driven thermal flows and their relevance with many technological and geophysical processes. Cooling systems for automobiles and nuclear reactors, as well as heat exchangers for turbomachinery and electronic devices, are examples of typical engineering applications. Flow simulations based on the Boussinesq approximation (BO) have proliferated in the literature since the seminal work of Gordon Draisey et al. [START_REF] Draisey | Three-dimensional natural convection in a box: a numerical study[END_REF] .

However the BO approximation is only valid for small temperature differences [START_REF] Gray | The validity of the Boussinesq approximation for liquids and gases[END_REF] (i.e., ∆T /T 0 < 0.1) and all thermal properties such as the expansion coefficient, viscosity and heat conductivity are held constant. Nevertheless, in many applications, non-Boussinesq conditions with variable thermodynamic properties prevail due to large density or temperature differences [START_REF] Frohlich | Large departures from Boussinesq approximation in the Rayleigh-Bénard problem[END_REF] . In this circumstance, the fully compressible Navier-Stokes system must to be considered. However the system's inherent fast acoustic waves are expensive to model numerically and moreover irrelevant in most convection flows at low speed. When compared to implicit schemes, explicit time discretization requires a smaller time step to meet the Courant-Friedrich-Lewy (CFL) stability condition, while implicit systems may have convergence issues due to poor conditioning of the convective fluxes of Jacobian matrix. Furthermore, the solution of the discrete equations contains pressure fluctuations of the order of the Mach number Ma, whereas the continuous equations only have pressure fluctuations on the order of Ma 2 , resulting in a loss of accuracy [START_REF] Paolucci | On the filtering of sound from the Navier-Stokes equations[END_REF] Therefore, low Mach thermal flow simulations remain challenging despite ongoing development due to the issue of finding the appropriate trade-off between algorithm accuracy and numerical efficiency.

One common and simple strategy to improve the computational efficiency is the Pressure Gradient Scaling (PGS) method [START_REF] Ramshaw | Pressure gradient scaling method for fluid flow with nearly uniform pressure[END_REF] with an acceleration coefficient α PGS > 1. This method artificially reduces the effective speed of sound by a factor α PGS , which automatically increases the time step by nearly the same factor. However, the PGS method has an inherent disadvantage that stems from the modification of one dominant term in the governing equation, which can be a serious issue, particularly when the flow dynamics are controlled by external forces such as buoyancy [START_REF] Wang | Artificial acoustic stiffness reduction in fully compressible, direct numerical simulation of combustion[END_REF] . This is not the case of another well-suited approach for dealing with heat dominated flows at arbitrary large density differences -the so-called low Mach number approximation (LMNA). This approximation filters the acoustic waves from the Navier-Stokes equations using an asymptotic expansion in powers of the Mach number. In practice, at least two scaling methods are possible depending on the temperature variation [START_REF] Zank | Nearly incompressible hydrodynamics and heat conduction[END_REF][START_REF] Bayly | Density variations in weakly compressible flows[END_REF][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF] , while the truncating order O(Ma 2 ) is widely considered and adopted in this investigation [START_REF] Paolucci | On the filtering of sound from the Navier-Stokes equations[END_REF][START_REF] Chenoweth | Natural convection in an enclosed vertical air layer with large horizontal temperature differences[END_REF] . The pressure p(x,t) is then splitted into a thermodynamic pressure p(t), which varies with time but is spatially homogeneous, and a hydrodynamic pressure p h (x,t), i.e., p(x,t) = p(t) + O(Ma 2 )p h (x,t). The total pressure is substituted by the thermodynamic pressure p(t), except in the momentum equation, which decouples the pressure and density fluctuations and filters the sound waves from the system. Particularly, only thermodynamic pressure is taken into account in the law of state. The velocity field is no longer divergence free, and the spatial and temporal variations of density introduce additional nonlinearities into the equations. Nonetheless, similar numerical approaches can still be employed because these equations have the same mixed hyperbolic-parabolic character as the equations for incompressible flow. The low Mach equations have already been applied successfully to research on natural convection, combustion, astrophysical flows and others (see some relevant examples in Bouloumou et al. 11 ).

In this study, we try to include the LMNA strategy into the Lattice Boltzmann method (LBM) and use this approach to simulate low-Mach number thermal flows. The decision to choose LBM is straightforward because of its attractive inherent properties such as simplicity to implement and parallelize, which provides a high computational efficiency for unsteady complex flows [START_REF] Krüger | The Lattice Boltzmann method[END_REF][START_REF] Guo | Lattice Boltzmann method and its application in engineering[END_REF] .

Originally designed for weakly compressible isothermal approaches, many attempts have been conducted to extend LBM capabilities to deal with more complex configurations such as high subsonic and supersonic flow dynamics [START_REF] Feng | Hybrid recursive regularized thermal Lattice Boltzmann model for high subsonic compressible flows[END_REF][START_REF] Feng | A compressible Lattice Boltzmann finite volume model for high subsonic and transonic flows on regular Lattices[END_REF][START_REF] Coratger | Large-eddy Lattice-Boltzmann modeling of transonic flows[END_REF][START_REF] Renard | Improved compressible hybrid Lattice Boltzmann method on standard Lattice for subsonic and supersonic flows[END_REF][START_REF] Farag | A pressure-based regularized Lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Zhao | Toward fully conservative hybrid Lattice Boltzmann methods for compressible flows[END_REF] , compressible thermal problems [START_REF] Feng | A three dimensional Lattice model for thermal compressible flow on standard Lattices[END_REF][START_REF] Feng | Regularized thermal Lattice Boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Wang | Large temperature difference heat dominated flow simulations using a pressure-based Lattice Boltzmann method with mass correction[END_REF] , turbulent flows [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large Eddy simulation[END_REF] , reactive flows [START_REF] Feng | A Lattice-Boltzmann model for low-Mach reactive flows[END_REF][START_REF] Muhammad | Development of combustion modelling within Lattice Boltzmann framework[END_REF] and others. For small temperature differences, thermal effects can be included using the BO approximation (see for example in litteratures [START_REF] He | A novel thermal model for the Lattice Boltzmann method in incompressible limit[END_REF][START_REF] Guo | A coupled Lattice BGK model for the Boussinesq equations[END_REF][START_REF] Guo | Thermal Lattice Boltzmann equation for low Mach number flows: Decoupling model[END_REF] ) which corresponds to a decoupling model and resolves an incompressible system with a linearized buoyancy force. In other configurations (corresponding to the fully compressible Navier-Stokes system), thermal effects are included either using another distribution function (Double Distribution Function -DDF approach) or either by solving an energy equation (hybrid approach) with a perfect gas law [START_REF] Feng | Regularized thermal Lattice Boltzmann method for natural convection with large temperature differences[END_REF][START_REF] Li | Coupling Lattice Boltzmann model for simulation of thermal flows on standard Lattices[END_REF] .

Most of these approaches investigate laminar flows and have been rarely extended to turbulent complex flows in low-Mach number limit.

Being intrinsically based on explicit time marching schemes, LBM's time-stepping is constrained by the Courant-Friedrich-Lewy (CFL) number (i.e., CFL ac ≡ max (|u|+c)∆t ∆x ≤ CFL max ). Thus, when using a compressible model to simulate a low-Mach number flow, the time-step is strongly constrained by the fast acoustic waves scale, which strongly penalizes the numerical efficiency of the solver [START_REF] Meister | Asymptotic single and multiple scale expansions in the low Mach number limit[END_REF] . To overcome this issue, several LB algorithms have been developed in the LMNA, by modifying classical isothermal LB algorithm, i.e., zeroth-order moment of the hydrodynamic pressure [START_REF] Filippova | Lattice-BGK model for low Mach number combustion[END_REF][START_REF] Lee | A Lattice Boltzmann algorithm for calculation of the laminar jet diffusion flame[END_REF][START_REF] Yamamoto | Simulation of combustion field with Lattice Boltzmann method[END_REF][START_REF] Seyed | Hybrid Lattice Boltzmann-finite difference model for low Mach number combustion simulation[END_REF][START_REF] Hosseini | Low-Mach hybrid Lattice Boltzmann-finite difference solver for combustion in complex flows[END_REF][START_REF] Cao | Variable property-based Lattice Boltzmann flux solver for thermal flows in the low Mach number limit[END_REF] . In such algorithms, additional terms need to be added to correctly recover Navier-Stokes equations in LMNA, thus making the resolved system complicated to understand and difficult to extend to other applications. This paper proposes a new algorithm for the hybrid LBM in the LMNA for simulating low-Mach number thermal flows. This low-Mach number algorithm only involves minor changes in the compressible LBM that makes it easy to implement and nearly independent of the thermal LBM used. The present method is evaluated on well-documented test cases of the literature for both simple laminar thermal flows and complex turbulent heat transfer problems. The 1D gravity column and the 2D rising bubble are considered to show the superiority of the new algorithm on the Pressure Gradient scaling commonly used to improve the computational efficiency of LBM. Then, a classical 2D non-Boussinesq natural convection flow is considered in a differentially heated square cavity at Rayleigh number 10 7 . The 3D Taylor-Green vortex flow at Reynolds number Re = 1.2 × 10 4 is chosen to assess the ability of the present method to evaluate viscous tensor. The final test case is the flow around a 3D heated cylinder at Re = 3900, to evaluate the performance of the present method to predict complex thermal turbulent flow.

The rest of the paper is organized as follows: §II introduces the general low-Mach number acceleration strategy adopted in this study and details the incorporation of low-Mach approach within the hybrid compressible LBM framework. Then, in §III the code is benchmarked with respect to numerical data from the literature for five test cases. At last, a conclusion with a summary of the speedup is drawn in §IV.

II. THE HYBRID LBM SOLVER IN THE LOW-MACH NUMBER APPROXIMATION

As mentioned in the introduction, the low-Mach number approximation (LMNA) of the Navier-Stokes equations is a very popular and efficient way to simulate low-Mach number flows at a macroscopic level (Appendix A). However, it involves a Poisson equation to resolve the hydrodynamic pressure p h which makes it not directly applicable in the LBM framework while keeping the locality feature.

We consider here the hybrid LBM for solving compressible flows and its extension to the low-Mach number approximation proposed in this paper.

A. The hybrid LBM solver for compressible flows

A complete description of the hybrid LBM is provided in literature. [START_REF] Farag | A pressure-based regularized Lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid Lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF][START_REF] Farag | Consistency study of Lattice-Boltzmann schemes macroscopic limit[END_REF] . In this method, the mass and momentum equations are solved using LBM, while the entropy equation is resolved using Finite Difference method (FDM). Three major steps detailed below are needed to obtain the variables at time step t + ∆t from the variables at t: a local collision step, a linear streaming step and the updating of macroscopic variables. [START_REF] Draisey | Three-dimensional natural convection in a box: a numerical study[END_REF] Thus, the three steps of the algorithm write as:

• Step I: Collision. The collision distribution function is calculated as

f col i (x,t) ≡ f eq i (x,t) + 1 - ∆t τ f neq i (x,t) + ∆t 2 F i (x,t) , (1)
where I.1 The equilibrium distribution f eq i is constructed in a regularized manner on the basis of the lattice, i.e.

f eq i = ω i    a (0),eq + H (1) i,α c 2 s a (1)
   , (2)
where ω i is the weight for direction i, c s the lattice sound speed (∆x/( √ 3∆t) in D3Q19 lattice), H the Hermite polynomial basis with α, β ,γ,δ being the coordinate indices.

They are all constants only depending on the lattice applied. a eq are the coefficients of the equilibrium distribution projected on the basis H . They can be different from one LBM kernel to another, but in any case, they only depend on local density ρ, velocity u, and a normalized temperature θ ≡ p ρc 2 s where p represent the pressure. Detailed expressions of the basis as well as the coefficients can be found in Farag et al. [START_REF] Farag | A pressure-based regularized Lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid Lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] .

I.2 The off-equilibrium distribution f neq i
is evaluated on second and third order basis as

f neq i = ω i   H (2) i,αβ
where the coefficients of the non-equilibrium population on the basis a is evaluated through a hybrid recursive procedure introduced by Jacob et al. [START_REF] Jacob | A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large Eddy simulation[END_REF] .

I.3

The forcing term F i is added to recover correct effects of viscous and external force (e.g. gravity),

F i ≡ F E i + F g i , (4)
in which the force F E i correct the higher order defects of the lattice to maintain the required stress tensor as in Eq. (A3c). Its exact expression can be found in [START_REF] Farag | A pressure-based regularized Lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid Lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF][START_REF] Farag | Consistency study of Lattice-Boltzmann schemes macroscopic limit[END_REF] . The gravity is considered through

F g i = ω i   H (1) i,α c 2 s ρg α + H (2) i,αβ 2c 4 s (ρu α g β + ρu β g α )   . (5)
•

Step II: Streaming. The collision distribution is streamed to its corresponding neighbor, i.e.

f i (x,t + ∆t) = f col i (x -c i ∆t,t) (6)
where c i is the vector pointing from the current node to its i th neighbor.

•

Step III: Macroscopic variable update.

III.1

The entropy equation is discretized using a explicit Euler time stepping, together with the flux and viscous terms calculated with FDM.

III.2

The density and momentum fields are updated using the distributions

ρ(x,t + ∆t) = ρ(x,t) + ∑ i f i (x,t + ∆t) -f col i (x,t) (7a) ρu α (x,t + ∆t) = ∑ i c iα f i (x,t + ∆t) + ∆t 2 ρ(x,t + ∆t)g α (7b)
We propose now to adapt this hybrid LBM solver in the low-Mach number approximation.

B. The hybrid LBM solver in the low-Mach number approximation

The LMNA is adapted here to the compressible LBM by analyzing the classical iso-thermal LBM. In fact, in athermal LBM, the "pressure" is nothing but a scaled density such as

p LBM ≡ ρc 2 s ( 8
)
It is updated at each time step by

p LBM (x,t + ∆t) ≡ c 2 s ρ(x,t + ∆t) = c 2 s ρ(x,t) + Q ∑ i=1 f col i (x -c i ∆t,t) -f col i (x,t) ≈ p LBM (x,t) -∆tc 2 s ∂ ρu α ∂ x α (x,t) = p LBM (x,t) -∆t ∂ p LBM u α ∂ x α (x,t) = p LBM (x,t) -∆t u α ∂ p LBM ∂ x α (x,t) -ρc 2 s ∂ u α ∂ x α (9)
As pointed out by He et al. [START_REF] He | Comparison of the Lattice Boltzmann method and the artificial compressibility method for Navier-Stokes equations[END_REF] , under low-Mach assumption, this equation mimics the acoustic pressure transport except that the physical speed of sound c in the last term is substituted by the lattice sound speed c s . This makes the isothermal LBM an artificial compressibility method [START_REF] Joel | A numerical method for solving incompressible viscous flow problems[END_REF] where the transport of acoustic pressure fields are modified to adapt the permitted time step of the numerical schemes [START_REF] He | Comparison of the Lattice Boltzmann method and the artificial compressibility method for Navier-Stokes equations[END_REF] .

We now extend the same idea to thermal LBM. Unfortunately, the temperature/density gradients in low-Mach thermal flow fields can be relatively large, such that the flow dilatation ∂ (ρu α )/∂ x α might break the scale balance in equation (9). At this stage, we just remove its "thermal" part (-∂ ρ/∂t), and only apply its "hydrodynamic" part for the hydrodynamic pressure transport, i.e.

p h (x,t + ∆t) = p h (x,t) + c 2 s Q ∑ i=1 f col i (x -c i ∆t,t) -f col i (x,t) -[ρ(x,t + ∆t) -ρ(x,t)] ≈ p h (x,t) -∆tc 2 s ∂ (ρu α ) h ∂ x α ≈ p h (x,t) -∆t u α ∂ ρ h c 2 s ∂ x α + ρ h c 2 s ∂ u α ∂ x α + u h α ∂ ρc 2 s ∂ x α + ρc 2 s ∂ u h α ∂ x α (10)
Compared to the theoretical transport equation of the hydrodynamic pressure, one can find again that the speed of sound is substituted by c s . Similar to the athermal LBM, the CFL constrain of such system is now related to c s instead of c. In practice, we still want to conserve the low-Mach property, i.e.,

c s ≡ ∆x √ 3∆t ≥ 3|u| max . (11)
Now the time step is constrained by the flow speed instead of the speed of sound. In low-Mach number applications, this will give us a much larger ∆t that will accelerate the calculation. Similar to every other LMNA method, this approach only modifies small scale quantities under low-Mach conditions, so that it will not suffer from the problem of PGS method when external forces are under consideration.

Such strategy can be easily applied to any compressible LBM. The following algorithm shows how all variables are advanced from time step t to t + ∆t

•
Step I: Collision. It remains the same as for the compressible LBM, except that the reduced "temperature" for the low-Mach version is now defined as

θ ≡ p h ρc 2 s . (12)
•

Step II: Streaming. It remains the same as the compressible LBM.

•

Step III: Macroscopic variable updating.

III.1 The entropy equation ( A3d) is resolved in the same manner as in the compressible LBM.

III.2 The thermodynamic pressure p(t + ∆t).

In open systems, it can be set directly to the ambient pressure. For closed configurations, it can be evaluated through global mass conservation, see appendix B.

The density field is evaluated through the equation of state (EOS) (A3e), i.e.

ρ(x,t + ∆t) = p(t + ∆t) e s(x,t+∆t)/c v 1/γ (13)
The momentum fields are updated the same as in compressible solver.

The hydrodynamic pressure field is updated using equation ( 9), i.e.

p h (x,t +∆t) = p h (x,t)+c 2 s ∑ i f i (x,t + ∆t) -f col i (x,t) -[ρ(x,t + ∆t) -ρ(x,t)] (14)
This low-Mach algorithm involves three changes with respect on the compressible LBM of §II A:

• i) The hydrodynamic pressure is used instead of the pressure field equilibrium population evaluation.

• ii)The density is solved through EOS instead of using LBM populations.

• iii) The hydrodynamic pressure is transported in step III.

This means that it is almost independent of the original thermal LBM scheme. In practice, it will be successfully applied in this paper to the pressure-based, the modified density-based and the unified LBM framework [START_REF] Farag | A pressure-based regularized Lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid Lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] . This is a major benefit of the current implementation compared to other low-Mach LBM approaches based on iso-thermal LBM [START_REF] Seyed | Hybrid Lattice Boltzmann-finite difference model for low Mach number combustion simulation[END_REF][START_REF] Hosseini | Low-Mach hybrid Lattice Boltzmann-finite difference solver for combustion in complex flows[END_REF] . It also permits us to keep some ingredients of the original LBM without any modification. For instance, the turbulence and wall models [START_REF] Wang | Large temperature difference heat dominated flow simulations using a pressure-based Lattice Boltzmann method with mass correction[END_REF][START_REF] Coratger | Large-eddy Lattice-Boltzmann modeling of transonic flows[END_REF] keep working out in this low-Mach version.

III. VALIDATION TESTS

The proposed low-Mach thermal Lattice Boltzmann solver is now assessed by a variety of test cases related to low-Mach thermal flows. In order to show the efficiency of the present solver (LMNA), numerical results are compared to simulations obtained using a fully compressible solver, i.e., the classical Hybrid Recursive Regularization pressure-based HRR-p solver [START_REF] Farag | A pressure-based regularized Lattice-Boltzmann method for the simulation of compressible flows[END_REF][START_REF] Farag | A unified hybrid Lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] , and using the PGS method [START_REF] Ramshaw | Pressure gradient scaling method for fluid flow with nearly uniform pressure[END_REF] (The PGS method is realized by multiplying pressure by a factor 1/α 2

PGS

). The following validation tests are investigated:

• One-dimensional gravity column and 2D rising thermal bubble for which the fluid dynamics is dominated by pressure gradient and external force, i.e., buoyancy.

• 2D natural convection in a closed square cavity at Rayleigh number Ra = 10 7 .

• 3D Taylor Green Vortex flow at Re = 1.2 × 10 4 .

• 3D turbulent flow around a heated cylinder at Re = 3900.

A. One-dimensional gravity column

The one-dimensional gravity column is firstly considered to evaluate the ability of LMNA on capturing the balance between hydrodynamic pressure gradient with external force. The configuration consists of a 1D vertical column filled by air with pressure uniformly distributed. The gravity drives the fluid motion in the beginning till a converged static state is attained with a constant pressure gradient ∇p = ρ 0 g.

In order to emphasize the computational efficiency of LMNA, simulations are also performed using PGS method and the fully compressible (HRR-p) solver. For the latter, the time step is fixed at ∆t = 1.65 × 10 -4 s, whereas LMNA and PGS allow to consider a time step ten times larger equal to ∆t = 1.65 × 10 -3 that speeds up the computations of the same factor. Fig. 1 shows the pressure profiles obtained from LMNA, PGS and the fully compressible (HRR-p) solver simulations, the latter being considered as the reference solution. The LMNA result matches perfectly the reference profile, while the PGS result amplifies by a factor of about 100 compared to the reference profile due to the re-scaled parameter α PGS = 10.

B. 2D rising thermal bubble

The rising thermal bubble configuration is dominated by gravity and pressure with thermal effect. The numerical configuration consists of a hot thermal bubble with a core temperature equal to 1200K that emerges in an atmosphere temperature of 300K. The width and height of the computational domain is respectively 20m and 15m. The bottom surface is configured to be a non-slip adiabatic wall, while the side and top boundaries are set to outlets with given pressure.

The maximum time step adopted by the fully compressible solver is equal to ∆t = 5.77×10 -5 s, whereas LMNA and PGS allow to consider a time step around 25 times larger equal to ∆t = 1.44 × 10 -3 that speeds up the computations of the same factor. Buoyancy with a gravity g = 9.81m s -2

drives the fluid motion, the bubble rises and deforms during the rising process. The temperature contours and streamlines are shown on Fig. 2. The LMNA simulation provides a result very similar to the reference solution with only negligible differences. However, the PGS simulation crashes using the same time step as LMNA. These results are consistent with the drawback of PGS as pressure plays a dominant role balancing the external force and viscous effects during the thermal bubble rising process. Consequently, the change of pressure gradient in PGS method leads to either incorrect results, i.e., one-dimensional gravity column in §III A or simulation crash as in this test case. By contrast, LMNA is able to provide the same results as with the fully compressible (HRR-p) solver, but with a time step between 10 and 25 times larger that allows much more efficient computations in these two low-Mach number configurations.

C. 2D natural convection inside a square cavity

This validation test has been intensively investigated in the literature both experimentally and numerically [START_REF] Salat | Experimental and numerical investigation of turbulent natural convection in a large air-filled cavity[END_REF][START_REF] Mergui | Convection naturelle en cavité carrée différentiellement chauffée: investigation expérimentale à Ra= 1, 69× 109[END_REF][START_REF] Wan | A new benchmark quality solution for the buoyancydriven cavity by discrete singular convolution[END_REF][START_REF] Massarotti | Characteristic-based-split (CBS) algorithm for incompressible flow problems with heat transfer[END_REF][START_REF] Vierendeels | Numerical study of natural convective heat transfer with large temperature differences[END_REF][START_REF] Vierendeels | Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences[END_REF] . Many benchmarks exist and solutions are available in the literature [START_REF] Vierendeels | Numerical study of natural convective heat transfer with large temperature differences[END_REF][START_REF] Vierendeels | Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences[END_REF] difference is considered to be beyond the Boussinesq validity limit, with the two vertical walls heated respectively at T 1 = 240K and T 2 = 960 K. This corresponds to a normalized temperature difference equal to ε = T 1 -T 2 T 1 +T 2 = 0.6, which is indeed far beyond the Boussinesq validity limit ε < 0.1 [START_REF] Wan | A new benchmark quality solution for the buoyancydriven cavity by discrete singular convolution[END_REF][START_REF] Massarotti | Characteristic-based-split (CBS) algorithm for incompressible flow problems with heat transfer[END_REF] . The fluid considered is air, and the dynamic viscosity is computed using the Sutherland's The thermal physics is characterized by a dimensionless number:

law µ(T ) = µ 0 ( T T * ) 3/2 T * +S T +S with T * = 273K, S = 110.5K, µ 0 = 1.68 × 10 -5 kg m -1 s -2 .
Ra = gρ 2 0 ∆T L 3 Pr T 0 µ 2 0 ( 15
)
where ρ 0 , µ 0 is the density and dynamic viscosity at mean temperature

T 0 = T 1 +T 2 2 = 600K. The
Rayleigh number Ra = 10 7 is chosen as it's related to the strongest heat convection effects over the range 10 3 -10 7 which is widely investigated [START_REF] Vierendeels | Numerical study of natural convective heat transfer with large temperature differences[END_REF][START_REF] Vierendeels | Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences[END_REF] . The grid resolution is N = L ∆x = 400 in both directions, as proposed by Wang et al. [START_REF] Wang | Large temperature difference heat dominated flow simulations using a pressure-based Lattice Boltzmann method with mass correction[END_REF] . The time steps adopted in the present LMNA and the fully compressible solver (HRR-p) are respectively equal to ∆t = 2.32 × 10 -5 s and ∆t = 4.17 × 10 -7 s. There is thus a factor of about 55.6 in the time steps of the two simulations showing again the superior efficiency of the LMNA in this type of low-Mach number flow configuration.

The heat flux at the heated wall is estimated by the Nusselt number and by its average over the heated wall, defined as follows:

Nu(y) = L λ 0 ∆T λ ∂ T ∂ x | wall Nu = 1 L y=L y=0 Nu(y)dy (16)
Time evolution of the average Nusselt number is plotted in Fig. 4, over a total dimensionless time U re f t/L = 250, long enough to ensure the convergence as proposed by Wang et al. [START_REF] Wang | Large temperature difference heat dominated flow simulations using a pressure-based Lattice Boltzmann method with mass correction[END_REF] , where

U re f = Ra 0.5 µ(T 0 ) ρ 0 L ≈ 1.1ms -1
is the reference velocity at Ra = 10 7 . As shown in Fig. 4, the average Nusselt number provided by the LMNA simulation (Nu = 16.19)converges to the value of the benchmark solution (Nu = 16.24) [START_REF] Vierendeels | Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences[END_REF] that corresponds to a difference less than 0.3%.

The flow features including the temperature distribution and velocity profiles are displayed in Fig. 5. The boundary layer close to the heated wall is well captured by LMNA as shown in Fig.

5(c).

The high heat flux due to high temperature variations in the thermal boundary layer (see Fig.

5(b)

) is also correctly predicted compared to the benchmark solution [START_REF] Vierendeels | Benchmark solutions for the natural convective heat transfer problem in a square cavity with large horizontal temperature differences[END_REF] and the one obtained with the fully compressible (HRR-p) solver [START_REF] Wang | Large temperature difference heat dominated flow simulations using a pressure-based Lattice Boltzmann method with mass correction[END_REF] .

D. Taylor Green Vortex

The investigation of Taylor Green Vortex (TGV) is aimed at assessing the ability of the present LMNA method of predicting the viscous tensor. The results and the computational efficiency of the classical pressure-based method reported by Tayyab et al. [START_REF] Muhammad | Development of combustion modelling within Lattice Boltzmann framework[END_REF] are included as a reference.

Pure fluid component (air) is considered and the Reynolds number investigated in this study is

Re = 1.2 × 10 4 .
The flow field and hydrodynamic pressure field are initialized using an analytic N-S and Poisson solution respectively given by:

u x (x, y, z) = u 0 sin ( 2πx L ) cos ( 2πy
L ) cos ( 2πz L ) (17a) u y (x, y, z) = u 0 cos ( 2πx L ) sin ( 2πy
L ) cos ( 2πz L ) (17b) u z (x, y, z) = 0 (17c) p h (x, y, z) = 1 16 ρ 0 u 2 0 cos ( 4πx L ) + sin ( 4πy
L ) 2 + cos ( 4πz L ) (17d)
The initial velocity is u 0 = 1m s -1 , and the viscosity is set to µ = 0.625×10 -3 kg m -1 s -1 . The thermodynamic pressure is the atmospheric pressure, i.e., P 0 = 101325.0 Pa, and the temperature is at T = 300K. Thus, density is taken as ρ = P 0 /(rT ) = 1.172 kg m -3 where r = 287.15m 2 s -2 The time evolution of the z-component of the vorticity is illustrated in Fig. 6. As can be seen, the vortex roll-up happens in the early stage followed by stretch and breakdown of these vortex while preserving laminar property. Laminar turbulence transition is then observed, and the coherent structures are developed to fully turbulence structures. At the end of the flow evolution, the turbulence is dissipated. These results compare very favourably to the results of Tayyab et al. [START_REF] Muhammad | Development of combustion modelling within Lattice Boltzmann framework[END_REF] which are performed using a classical pressure-based method.

K
In order to quantitatively analyze the viscous effect, the overall kinetic energy and its dissipation rate are evaluated in the next. The dissipation rate is of great importance because it's associated with the stress tensor. The overall kinetic energy is defined as the integral of the of kinetic energy over the whole domain, and the dissipation rate of kinetic energy is the temporal derivative of the overall kinetic energy, i.e.,

κ(t) = 1 2 V (u 2 x + u 2 y + u 2 z )dV ε(t) = ∂ κ(t) ∂t (18)
where the u i are the velocity components with i ∈ (x, y, z) and V is the integration over the whole domain. Fig. 7 shows the time evolution of these two normalized quantities. Excellent agreements are obtained compared to the results of reference of Tayyab et al. [START_REF] Muhammad | Development of combustion modelling within Lattice Boltzmann framework[END_REF] . Furthermore, the x and ycomponents of the velocity profile at the central line are plotted at time 12.11T c in Fig. 8. The LMNA results match very well the solution of reference [START_REF] Muhammad | Development of combustion modelling within Lattice Boltzmann framework[END_REF] . In addition, even the small turbulent structures are well captured by the LMNA, as shown on the x-component of the velocity in Fig.

8(a)

. All theses results show that the present LMNA method is able to simulate turbulent flow and laminar turbulence transition processes while keeping high computational efficiency.

E. Flow around a circular heated cylinder at Reynolds number 3900

This section is dedicated to study the flow passing through a cylinder at Reynolds number Re = 3900 with isothermal solid wall in order to assess the capability of the present method in the LES (Large Eddy Simulation) framework [START_REF] Sagaut | Large Eddy simulation for incompressible flows: an introduction[END_REF] . A mesh convergence study is conducted firstly on an isothermal case where the temperature of solid wall is the same as the inflow. The results are benchmarked with experiments [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] and DNS (Direct Numerical Simulation) [START_REF] Dong | A combined Direct Numerical Simulation-Particle Image Velocimetry study of the turbulent near wake[END_REF] . In addition, the LES results reported by [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF][START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF] are also included to give a reference. Then, the thermal case is investigated with temperature difference between the solid wall and the inflow of ∆T = 300K. The thermal effects on the flow features and turbulent heat flux are studied comparing with the LES results of Sircar et al. [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF] and Jogee et al. [START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF] . In this study, we focus on the near wake region, for which a variety of experimental and numerical results are available in the literature for validation.

FIG. 9. Sketch of the grid refinement in the xy plane with seven levels of grid refinement. Black and blue lines show respectively the refinement boxes and the edges of each grid level.

Numerical setup

The sketch of the computational domain is shown in Fig. 9 together with the mesh. The cylinder has a diameter of D = 0.01m, and its center is chosen as the origin. The x, y, z directions are respectively the streamwise, transverse and spanwise direction regarding the inflow. The mesh in the transverse and spanwise directions as well as boundary conditions are symmetric with a spanwise extension of L z = 3D which is closed to the conventional L z = πD adopted by most LES studies [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF][START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF][START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF][START_REF] Dmitry A Lysenko | Large-eddy simulation of the flow over a circular cylinder at Reynolds number 3900 using the OpenFOAM toolbox[END_REF][START_REF] Mani | Computational study of optical distortions by separated shear layers and turbulent wakes[END_REF] . The transverse width is L y = 20D for each side, and the distances to the inlet and outlet are respectively equal to 15 and 50 diameters, which is sufficiently large to reduce the effects of spurious pressure wave reflecting from the boundary. The inlet velocity and temperature boundary conditions are constant equal to U ∞ = 5.968 m s -1 and T ∞ = 300 K, respectively. Zero pressure gradient boundary is used for the outlet, and symmetric transverse top and bottom wall.

Periodic boundary condition is provided in the spanwise direction. Isothermal no-slip boundary condition is imposed at the cylinder with two values of temperature equal to T w = 300 K and T w = 600 K to study the impact of the heat transfer.

Details on the grid refinement and related numerical parameters adopted in this study are given in Table . I. As Cartesian mesh is adopted by LBM, the mesh size is divided by a factor 2 for each grid refinement with the coarsest mesh size ∆x max /D = 0.5. The block refinement is the combination of cylinders and cuboids, as can be seen in Fig. 9 (cycle and rectangular for the 2D illustration). Three resolutions are chosen for this mesh convergence study with six, seven and eight grid levels. The number of mesh points are respectively equal to 1.88 millions, 6. For compressible flows, the Favre-filtering is more convenient to avoid the SGS (Sub-Grid Scale) terms introduced in continuity equation [START_REF] Garnier | Large Eddy simulation for compressible flows[END_REF][START_REF] Pino Martin | Subgrid-scale models for compressible large-eddy simulations[END_REF] . A Vreman sub-grid scale model is adopted [START_REF] Aw Vreman | An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications[END_REF] with the model constant C s = 0.18. A constant turbulent Prandtl number is chosen equal to Pr t = 0.85 as suggested [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF][START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF][START_REF] Ba Kader | Heat and mass transfer laws for fully turbulent wall flows[END_REF] . [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] . The results obtained with the middle mesh agree well with experimental data, except visible discrepancy of the maximum value of mean velocity in the re-circulation region. Finally, an excellent agreement is achieved with the fine mesh. In particular, the present LMNA result gives a re-circulation length L rc /D equal to 1.46, compared to 1.51 given by Parnaudeau et al. [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] and 1.47 given by Dong et al. [START_REF] Dong | A combined Direct Numerical Simulation-Particle Image Velocimetry study of the turbulent near wake[END_REF] that corresponds to an error of around 3% and less than 1%, respectively (see in Table II).

In the next, details of the flow features obtained for the finest mesh for the isothermal case are discussed including the mean flow, the turbulent characteristics and some aerodynamic coefficients. (a the laminar-turbulent transition. Thus, the large scale vortices break down to small scale vortices in the wake region. et al. [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] , the re-circulation length is correctly estimated by Sircar et al. [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF] with L rc,sircar /D = 1.48

) Q(D 2 /U 2 ∞ ) = 1 (b) Q(D 2 /U 2 ∞ ) = 100 FIG. 11. Instantaneous iso-surfaces of the Q criterion Q * = Q(D 2 /U 2 ∞ ) = 1 (
and by the present LMNA with L rc,present /D = 1.46. However, the overestimation of Jogee et al. [START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF] results in overestimating the peak value of the mean streamwise Reynolds stress ⟨u⟩/U ∞ , and in shifting the profiles of these quantities to the downstream direction. Moreover, the peak values of the mean streamwise velocity and corresponding Reynolds stress given by the present LMNA are closest to the experimental data. Furthermore, even in the relatively far wake region, i.e.,

x/D ∈ [2.5, 4.5], present method gives best agreement with the experimental data. The far wake region, i.e., x/D > 4.5, won't be intensively studied in this investigation because of the lack of experimental results. The deviation of Jogee's results are very significant compared to the experimental [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] and DNS [START_REF] Ma | Dynamics and low-dimensionality of a turbulent near wake[END_REF] data. However, the discrepancy between the present LMNA results and the LES results of Sircar et al. [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF] , or the experimental [START_REF] Parnaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF] and DNS [START_REF] Ma | Dynamics and low-dimensionality of a turbulent near wake[END_REF] results is rather small.

The streamwise Reynolds stress profiles show two peaks (Fig. 14

Non-isothermal results at ∆T = 300K

The cylinder is now heated to T w = 600 K. In order to investigate the influence of heat transfer on the flow features, the turbulent heat flux and temperature fluctuations are analyzed quantitatively at the same probe locations as in previous section. The non-isothermal LES results of Sircar et al. [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF] and Jogee et al. [START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF] at the same temperature difference ∆T = 300 K are included as references. The fine mesh defined in former section is used for the non-isothermal simulations here.

The time step is fixed at ∆t * = 4.49 × 10 -4 . Let's notice that similar computations with the fully compressible solver requires a time step 16.4 times smaller equal to ∆t * = 2.74 × 10 -5 . Table III shows that the re-circulation bubble length becomes larger when the cylinder is heated. This is probably caused by the temperature increase in the detached shear layer. As can be seen Compared to the isothermal simulations, the bubble length is increased of 21.4 % and 8.5 % for Jogee et al. [START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF] and Sircar et al. [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF] , respectively, indicating that the Jogee's LES results are more sensitive to the thermal effects than the Sircar ones. In the present simulation the increase is of 15.8%. However, as the isothermal bubble length predicted by the Jogee's LES results was already The sensitivity with thermal effect of the turbulent statistics is shown in Fig. 17(c) and Fig.

17 In the next, the thermal characteristics including the temperature distribution and turbulent heat flux will be discussed. previous description. It is probably due to the weaker influence of the thermal effects on the flow features.

These results show that the new LMNA provides promising results for complex turbulent isothermal configurations and consistent explanation for heat transfer problems with large temperature differences. Moroever, with a time step more than 15 times larger than for a fully compressible LBM solver, the computational efficiency is strongly improved.

IV. CONCLUSION

In this paper, we have proposed a new hybrid Lattice Boltzmann method for dealing efficiently with Low Mach thermal flows. To overcome the time step constraint of the fully compressible LBM solvers due to acoustic waves, the low-Mach number approximation has been considered and adapted to the compressible hybrid LBM solver developed by the team [START_REF] Farag | A unified hybrid Lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] . The present algorithm overcomes the drawback of the PGS method related to the modification of the pressure gradient. In addition, the present algorithm can be easily implemented to improve any density based or pressure based compressible solvers, which makes it more flexible than other existing LB algorithms based on LMNA.

The performance of present method are assessed by two classes of cases. The first category of test scenarios includes flows in laminar regime, i.e., one dimensional gravity column, 2D rising thermal bubble and 2D natural convection for which all the results agree well with the reference.

The second category consists of 3D Taylor Green Vortex and 3D flow passing through a heated cylinder in turbulent regime. Excellent agreements are obtained especially on the average flow features with respect to reference results of the literature. Some discrepancies are observed for the turbulent quantities of 3D heated cylinder compared to other LES results. The results of present method are between the LES results of Sircar et al. [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF] and Jogee et al. [START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF] and closer to the one of Sircar. This indicates that current solution is more sensible to the turbulent thermal effect than Sircar et al. [START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF] and less sensible than Jogee et al. [START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF] . This disparity is most likely due to the differences in turbulent models and mesh configurations adopted in these LES investigations.

Furthermore, table IV summarizes the computing speed-ups of all the test cases employed in this investigation. As can be seen, the computational efficiency is increased by at least a factor of ten.

∂ p ∂ x α = 0, (A3b)
∂ ρu α ∂t + ∂ ρu α u β ∂ x β = - ∂ p h ∂ x α + 1 Re ∞ ∂ Π αβ ∂ x β + 1 Fr 2 ∞ ρg α , (A3c)
ρT ∂ s ∂t + u α ∂ s ∂ x α = - 1 Re ∞ Pr ∞ ∂ q α ∂ x α , (A3d)
p = ρrT , s = c v ln p ρ γ . ( A3e
)
Together with the reference Reynolds, Froude and Prandtl numbers

Re ∞ ≡ ρ ∞ u ∞ x ∞ µ ∞ , Fr ∞ ≡ u ∞ √ g ∞ x ∞ , Pr ∞ ≡ c p ∞ µ ∞ λ ∞ , (A4)
In the above system, the thermodynamic pressure p is spatial independent according to equation (A3b). Also, one additional equation for the hydrodynamic pressure p h is necessary to resolve this set of equations. Typically, this is achieved through a Poisson equation constructed from the mass and momentum equations [START_REF] Joel H Ferziger | Computational methods for fluid dynamics[END_REF] .

FIG. 1 .

 1 FIG. 1. Comparison of the pressure profiles obtained from LMNA (a) and PGS (b) simulations and with the fully compressible (HRR-p) solver which is considered as the reference.

FIG. 2 .

 2 FIG. 2. Comparisons of 2D temperature contours together with streamlines showing the rising bubble. For each plot, the left hand side corresponds to the fully compressible (HRR-p) solver simulation considered as the reference. (a) on the right hand side, LMNA simulation. (b) on the right hand side, PGS simulation which is crashed.

  FIG. 3. Sketch of the 2D differentially heated square cavity together with boundary conditions and a Cartesian mesh.

  Then the thermal conductivity is calculated by λ (T ) = µ(T )C p Pr , where C p being the specific heat capacity of air, and Pr being the Prandtl number equal to 0.71.

  FIG. 5. Profiles of temperature ((a)-(b)) and velocity (c) at Rayleigh number Ra = 10 7 . Comparison between the present LMNA result (solid line), the fully compressible (HRR-p) result (dashed lines), and the benchmark solution 47 (markers).

FIG. 6 .

 6 FIG. 6. Time slots of the iso-surfaces of the z-component of the vorticity at 5.0s (a), 8.0s (b), 12.1s (c), 19.0s (d).

FIG. 7 .

 7 FIG. 7. Kinetic energy (a) and kinetic dissipation rate (b) evolution. T c = 1s, the kinetic energy and its dissipation rate are normalized respectively by the initial value and maximum value of corresponding quantity.

  FIG. 8. Profiles of x (a) and y (b) components of velocity at time t = 12.11T c .

52 millions and 25 .

 25 16 millions. The finest grid sizes normalized by the cylinder diameter distributed along the isothermal cylinder wall are correspondingly ∆x min /D = 1.56 × 10 -2 , ∆x min /D = 7.81 × 10 -3 and ∆x min /D = 3.91 × 10 -3 . The numbers of iterations required for a dimensionless time U ∞ t * /D = 1 are 1115, 2227 and 4464, respectively, that corresponds to CPU times (by hour) of 1.44, 30.7 and 194.68, respectively. In order to obtain statistically converged results, 500 FOT (Flow Over Time), .i.e, t = 500t * (around 100 vortex shedding cycles) is chosen as initial transient time. Then the turbulent statistics and mean flow validation are realized from data collected over 150 FOT, i.e., t = 150t * (around 30 vortex shedding cycles). In this section, brackets ⟨•⟩ and (• ′ ) denote time-averaged quantities and the corresponding fluctuations, respectively.

2 .

 2 Isothermal results at ∆T = 0 K It's observed by Parnaudeau et al. 49 that the main flow features and turbulent characteristics in the wake region is determined by the re-circulation length L rc , especially the transverse velocity and Reynolds stress. The length of this re-circulation bubble is defined as the distance from the rear stagnation point of the cylinder to the position where the sign change of the mean streamwise velocity along the centerline in the wake region. So the grid convergence study is focused on the prediction of this re-circulation zone. The mean streamwise velocity ⟨u⟩ normalized by the inlet velocity U ∞ and the streamwise Reynolds stress ⟨u ′ u ′ ⟩ normalized by U 2 ∞ are plotted in Fig. 10. The negative mean velocity corresponds to the re-circulation zone as shown in Fig. 10(a). The overestimated re-circulation length obtained by the coarse mesh simulation pushes the streamwise Reynolds stress to the downstream direction compared to the experimental data by Parnaudeau et al.

Fig. 11 shows

 11 Fig. 11 shows the laminar and turbulent structures of the flow. At this Reynolds number, the flow remains laminar along the cylinder wall and close to the cylinder just after the boundary layer separation. Large scale vortices develop at around two diameters from the rear stagnation point and are convected by the mean flow thus shedding alternatively from the upper and lower shear layers forming the well-known Von Karman street. The instability of the shear layer contributes to

  FIG. 11. Instantaneous iso-surfaces of the Q criterion Q * = Q(D 2 /U 2 ∞ ) = 1 (a), and Q(D 2 /U 2 ∞ ) = 100 (b), colored by the velocity magnitude.

Fig. 12

 12 Fig. 12 shows the time-averaged velocity field and the associated streamlines. As expected, the mean flow is globally symmetric to the axial line in the xy-plane, indicating that the solution is convergence over the computational time. The re-circulation zone is confined within the top and bottom shear layers. Two transverse lines in the re-circulation region at positions X/D = 1.06 and X/D = 1.54 are probed where the shear layer dynamics and turbulent statistics are quantitatively analyzed.As shown in Fig.13(a) and Table.II, compared to the experimental results of Parnaudeau

FIG. 12 .

 12 FIG. 12. Streamlines and velocity field of the mean flow at Re = 3900. Isothermal case.

  FIG. 14. (a, b) Mean streamwise and transverse velocity components. (c, d) Streamwise and transverse Reynolds stress components at position x/D = 1.06 and 1.54. The reference data from experiments PIV case II of Parnaudeau et al. 49 and DNS case II of Ma et al. 61 are included for comparison.

  (c)) which correspond to the detached shear layers. The instability induced by shear phenomena produces highest perturbation near the shear layers. The streamwise Reynolds stress are three times larger at the shear layer than the values in the center of the re-circulation bubble. The two peaks show a diffusive trend as the probe moves downstream, which is probably influenced by the oscillations of the alternative vortex shedding phenomena. The values of the transverse Reynolds stress components are around three times larger than the streamwise ones in the center of the re-circulation bubble because of the oscillations in the transverse direction of the vortex shedding phenomena.

  FIG. 15. Color maps of the mean velocity magnitude (a) and temperature (b).

  FIG. 16. Profiles along the centerline in the wake region of the mean streamwise velocity (a) and streamwise Reynolds stress (b) components. Non-isothermal LES results of Sircar et al. 51 and Jogee et al. 52 at the same temperature difference ∆T = 300 K are included as references.

  FIG. 17. Mean streamwise and transverse velocity components (a, b), streamwise and transverse Reynolds stress components (c, d) at x/D = 1.06 and 1.54. Non-isothermal LES results of Sircar et al. 51 and Jogee et al. 52 at the same temperature difference ∆T = 300 K are included as references.

  FIG. 18. Mean temperature profiles (a) and RMS of temperature fluctuations (b) at x/D = 1.06 and x/D = 1.54. Non-isothermal LES results of Sircar et al. 51 and Jogee et al. 52 at the same temperature difference ∆T = 300 K are included as references.

  Fig. 15(b). The temperature fluctuations peaks correspond to the peaks on the mean temperature profiles due to the temperature amplitude and the strong instability at the shear layers. In the present LMNA results, the mean temperature at the re-circulation bubble is slightly higher than the LES results of Sircar et al. and Jogee et al., whereas the temperature fluctuations are very close. Fig. 19 shows the color map of the streamwise ⟨u ′ T ′ ⟩/U ∞ T ∞ and the transverse turbulent heat flux ⟨v ′ T ′ ⟩/U ∞ T ∞ where the positive (red) and negative (blue) zones are highlighted.

Fig. 20 FIG. 21 .

 2021 Fig. 20 shows the temperature fluctuations, and the streamwise heat flux distribution along the center line. As can be seen, the RMS of temperature fluctuations and turbulent heat flux attain a peak value near the rear point of the cylinder, showing that the thermal fluctuations are most significant at this position. The turbulent heat flux is positive in the re-circulation region, and switches to a negative value outside the re-circulation region. This means that the heat transfer is from the heated cylinder to the cold fluid in the re-circulation zone. The position where the turbulent heat flux becomes zero is very close to the margin of re-circulation zone, and this phenomenon is consistent for all the three LES results: L rc,Sircar < L rc,current < L rc,Jogee . The streamwise and transverse turbulent heat flux profiles along the two probe lines are plotted in Fig. 21. The heat transfer direction is represented by the sign of the turbulent heat flux. The steep gradient of temperature at the shear layer produces the four peaks of each profile with opposite sign near each other. This is consistent with the color map of Fig. 19, where the highlighted blue zone and red zone are distributed alternatively near the shear layer. Moreover, the transverse turbulent heat flux at center line (y/D = 0) is zero because the transverse velocity along the center line is zero inside the re-circulation region. Lastly, the discrepancy between the results of Sircar et al. and the present ones is not significant for both the peak values and the profile shapes. However, the deviation with the LES results of Jogee et al. is considerably larger, which is consistent with

TABLE I .

 I Details of meshes and computational setup, t * = D/U ∞ = 1.68 × 10 -3 s is the FOT (Flow Over Time) and it's used as a dimensionless quantity.

		Mesh Total cell points (M) Grid level	∆x min /D	∆t min /t *	CFL u+c t w (h)
		Coarse	1.88	6	1.56 × 10 -2 8.97 × 10 -4	0.684	1.44
	Isotherm	Middle	6.52	7	7.81 × 10 -3 4.49 × 10 -4	0.684	30.7
		Fine	25.16	8	3.91 × 10 -3 2.24 × 10 -4	0.684 194.68
	Non-isotherm	Fine	25.16	8	3.91 × 10 -3 4.49 × 10 -4	0.84	44.26

TABLE II .

 II Experimental and numerical results of the literature together with present LMNA results. EXP stands for experiment.

	Case	Model	-Cpb	⟨C d ⟩	St	L rc /D
	Parnaudeau et al. 49	EXP (PIV)	-	-	0.208	1.51
	Dong et al. 50	EXP (PIV)	-	-	0.212	1.47
	Norberg et al. 59,60	EXP (Multiple) 0.85-0.95	0.93-1.05	0.21	1.4-1.6
	Ma et al. 61	DNS (case II)	0.84	-	0.219	1.59
	Dong et al. 50	DNS	0.93-1.04	-	0.206-0.210	1-1.18
	Lysenko et al. 53	LES	0.91	0.97	0.209	1.67
	Mani et al. 54	LES	0.86	0.99	0.206	-
	Franke and Frank et al. 62	LES	0.85-0.94 0.978-1.005	0.209	1.34-1.64
	Kravchnko and Moin et al. 63	LES	0.94	1.04	0.210	1.35
	Parnaudeau et al. 49	LES	-	-	0.208	1.56
	Jogee et al. 49	LES	-	-	-	1.68
	Sircar et al. 51	LES	0.91	1.048	0.221	1.48
	Coarse mesh	LES	0.81	1.06	0.212	1.67
	Middle mesh	LES	0.86	1.04	0.223	1.41
	Fine mesh	LES	0.84	1.04	0.223	1.46
	Fine mesh (non-isothermal)	LES	0.79	1.00	0.212	1.69

TABLE III .

 III Recirculation lengths for the isothermal and the non-isothermal cases. Present LMNA results are given together with LES predictions of the literature by Sircar et al.[START_REF] Sircar | Turbulent flow and heat flux analysis from validated Large Eddy simulations of flow past a heated cylinder in the near wake region[END_REF] and Jogee et al.[START_REF] Jogee | Large-Eddy simulation of nonisothermal flow over a circular cylinder[END_REF] .

		L rc /D
	Case			Bubble length change (%)
		Isothermal	∆T = 300K
	Jogee et al. 52	1.68	2.04	21.4
	Sircar et al. 51	1.480	1.606	8.5
	Current	1.46	1.69	15.8

TABLE IV .

 IV Computations speedup for the different test cases of the paper. The speedup is evaluated with respect to fully compressible LBM computations.

	Then, the leading order Navier-Stokes equations read 2
	∂ ρ ∂t	+	∂ ρu α ∂ x α	= 0 ,		(A3a)
	Test cases					Speedup
	1D Gravity column at Ma =	√ gL c ≈ 0.029	10.0
	2D Rising bubble at Ma = 0.018	25.0
	2D Natural convection at Ma = 0.0035	55.6
	3D Taylor Green Vortex at Ma = 0.0029	70.3
	3D Heated cylinder at Ma = 0.017	16.4

The equations in the algorithm are normalized by LBM units, i.e. x ∞ ≡ ∆x, t ∞ ≡ ∆t.

The superscript 0 for the variables with 0 t h order is neglected here.

ACKNOWLEDGEMENTS

This work was performed within ProLB solver. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. We also acknowledge the supports of ANR, Renault, Airbus and SafranTech by the Industrial Chair Program ALBUMS (ANR-CHIND-18-ALBUMS).

Appendix A: The low Mach number approximation of the Navier-Stokes equations

The low-Mach number acceleration strategy based on the scaling of fully compressible Navier-Stokes equations in low-Mach number limit is presented in this section. One of the main challenges when simulating low-Mach number thermal flows using a compressible solver is to keep high computational efficiency despite the non-dominant acoustic scale fluctuations. A method to cure this issue is based on a scale analysis, which divides the fully compressible Navier-Stokes equations by a leading order fluctuations thermal dynamics and a non-dominant acoustic scale.

Based on the asymptotic analysis, any variable φ appears in the fully compressible Navier-Stokes system will be expanded as