M Nguyen
email: nguyen@cerfacs.fr

J F Boussuge

P Sagaut

J C Larroya-Huguet

Large Eddy Simulation of a Thermal Impinging Jet using the Lattice Boltzmann Method Large Eddy Simulation of a Thermal Impinging Jet using the Lattice Boltzmann Method

A compressible Hybrid Lattice Boltzmann Method (LBM) solver is used to perform a wall-resolved Large Eddy Simulation (LES) of an isothermal axisymmetric jet issuing from a pipe and impinging on a heated flat plate at a Reynolds number of 23 000, a Mach number of 0.1, and an impingement distance of two jet diameters. The jet flow field statistics, Nusselt number profile (including the secondary peak), and shear stress profile were well reproduced. The azimuthal coherence of the primary vortical structures was relatively low, leading to no discernible temporal periodicity of the azimuthally-averaged Nusselt number at the location of the secondary peak. While local unsteady near-wall flow separation was observed in the wall jet, this flow-separation did not exhibit azimuthal coherence, and was not found to be the only cause of the thermal spots bluewhich lead to the secondary peak in the Nusselt number, as stream-wise oriented structures also played a significant role in increasing the local heat transfer.

I. INTRODUCTION

Impinging jet flows are widely used as a means of cooling or heating surfaces thanks to their high rates of convective heat transfer. Due to the thin boundary layers involved at the impinged surface, as well as the strong turbulent mixing, impinging jets can be several orders of magnitude more efficient as a method of heat transfer than parallel fluid flow [START_REF] Zuckerman | Jet impingement heat transfer: Physics, correlations, and numerical modeling[END_REF] . Heat transfer from an impinging jet can highly be localized, hence often necessitating multiple arrays of jets. Impinging jets can be found in diverse applications such as turbine blade cooling, turbine adaptive clearance control systems [START_REF] Ahmimache | Heat transfer and flow velocity study of a row of jets emerging from a perforated pipe at a low reynolds number[END_REF] , aircraft wing anti-icing configurations, and aircraft engine thrust reversers.

The turbulent axisymmetric impinging jet remains, despite its relative geometric simplicity, a widely studied academic test configuration that has been subject to many literature reviews [START_REF] Jambunathan | A review of heat transfer data for single circular jet impingement[END_REF][START_REF] Dewan | Recent trends in computation of turbulent jet impingement heat transfer[END_REF][START_REF] Carlomagno | Thermo-fluid-dynamics of submerged jets impinging at short nozzle-to-plate distance: A review[END_REF][START_REF] Shukla | Flow and thermal characteristics of jet impingement: comprehensive review[END_REF] , as it exhibits complex flow physics. It is characterized by three distinct regions [START_REF] Gauntner | Survey of literature on flow characteristics of a single turbulent jet impinging on a flat plate[END_REF] : a free jet zone, a stagnation zone, and a wall jet zone. In the free jet zone, the flow is sufficiently removed from the impacted surface such that the surface does not affect the flow dynamics. Further downstream, in the stagnation zone, the jet decelerates, leading to a region of high pressure, and the flow is deflected from the axial direction towards the radial direction. In the wall jet zone, the radial flow parallel to the wall predominates. Thanks to a favorable pressure gradient, the flow initially accelerates, until radial spreading and conservation of mass decelerate the flow.

One of the salient features of impinging jet flows with low nozzle-to-plate distances is the presence of two peaks of the Nusselt number with respect to the radius. A first peak is found directly in the vicinity of the impacted location, and a second one is found further downstream, in the wall jet. The use of traditional Reynolds-Averaged Navier-Stokes (RANS) approaches to properly simulate impinging jet flows has shown mixed results in reproducing this secondary peak, and often requires highly calibrated turbulence models and wall treatments [START_REF] Zuckerman | Jet impingement heat transfer: Physics, correlations, and numerical modeling[END_REF] . Thanks to increases in computational power, a number of large eddy simulations (LES) have been performed on this configuration. Hadžiabdić and Hanjalić 8 , Uddin et al. [START_REF] Uddin | LES simulations of an impinging jet: On the origin of the second peak in the Nusselt number distribution[END_REF] , Aillaud et al. [START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF] , Grenson and Deniau [START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] , and Colombié et al. [START_REF] Colombie | Large-eddy-simulation-based analysis of reynolds-stress budgets for a round impinging jet[END_REF] demonstrated that with a sufficiently resolved mesh, LES approaches could reliably reproduce the second peak, and could be used to extensively study the jet flow dynamics and their effects on heat transfer. The Lattice Boltzmann Method (LBM) has seen increasing scrutiny as a means of performing simulations with LES fidelity at a lower computational cost [START_REF] Löhner | Towards overcoming the LES crisis[END_REF] . Thanks to the collideand-stream algorithm, LBM calculations can often be significantly more rapid than with typical Navier-Stokes solvers for a given mesh size. The use of Cartesian octree meshes, although not without its disadvantages, can make pre-processing and mesh generation significantly easier to automate, greatly reducing the amount of working time needed by a highly trained specialist. LBM-LES approaches have been sucessfully applied to industrial applications [START_REF] Hou | Lattice-Boltzmann and Navier-Stokes simulations of the partially dressed, cavity-closed nose landing gear benchmark case[END_REF][START_REF] Daroukh | Low-speed turbofan aerodynamic and acoustic prediction with an isothermal lattice boltzmann method[END_REF] .

Classical LBM models are limited to weakly compressible isothermal flows [START_REF] Krüger | The Lattice Bolzmann Method[END_REF] . Recently, new compressible thermal LBM models [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Guo | An efficient lattice boltzmann method for compressible aerodynamics on d3q19 lattice[END_REF][START_REF] Guo | Improved standard thermal lattice boltzmann model with hybrid recursive regularization for commpressible laminar and turbulent 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 unified hybrid lattice-boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] have been developed that preserve the method's rapidity and accuracy while also being robust and practical to implement into an industrial solver. These models make use of the hybrid formulation, where the conservation of mass and momentum are computed using the collide-and-stream algorithm, and the conservation of energy is calculated using a finite volume or finite difference scheme. These two methods are coupled through the definition of the LBM equilibrium distribution function, which is made dependent on temperature. The hybrid method is combined with correction terms to eliminate the well known Mach error in the stress tensor, along with a Hybrid Recursive Regularized (HRR) collision operator [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann method-based large eddy simulation[END_REF] to improve stability.

The simulation of a 3D turbulent thermal impinging jet using an LBM-LES approach requires numerous features. The pipe inlet must mimic fully developed pipe flow, necessitating turbulence injection and non-reflecting boundary conditions. Due to the widely varying scales of the flow structures, local mesh refinement zones must be used to limit computational cost, requiring sophisticated grid-transition algorithms which can handle the sudden change in mesh size induced by the octree structure. The challenging nature of this test case makes it a useful one for the assessment of a numerical solver. The purpose of the present work is to demonstrate the improved maturity of LBM approaches for turbulent thermal test cases. Although this is an incompressible test case, it presents a useful step for the validation of the fully compressible solver for complex flows.

This paper is organized as follows. The hybrid LBM approach for LES, together with the treatment of boundary nodes, resolution domain transitions, and turbulent inflow conditions, are described in Section II. Section III describes the numerical setup of the test case. The flow field statistics, together with the relevant surface quantities, are then validated in Section IV. Lastly, an examination of the flow structures and their effects on the surface quantities is performed in Section V.

II. MODEL OVERVIEW

A. The Lattice Boltzmann Equation

In this section, the variables ρ, u, θ = T dim /T re f , p, and µ (density, velocity, temperature, pressure, and viscosity, respectively) are presented in their non-dimensional form unless otherwise noted by the subscript dim. They are nondimensionalized using a reference density ρ 0 , the local mesh size ∆x, the user-set reference temperature T re f , and the local time step ∆t = (∆xc s )/( RT re f ), where R is the specific gas constant and c s = 1/3 is the lattice constant such that nondimensional pressure is p = ρc 2 s θ . The Lattice Boltzmann collide-and-stream algorithm, when using the HRR collision model [START_REF] Jacob | A new hybrid recursive regularised bhatnagar-gross-krook collision model for lattice boltzmann method-based large eddy simulation[END_REF] , is as follows:

f i (x α + c iα ,t + 1) = f eq i (x α ,t)+ (1 - 1 τ )R( f neq i ) + 1 2 ψ i (x α ,t), (1)
where f i is the particle distribution function, with the subscript i representing the index of the discrete lattice velocities c iα . Here, the D3Q19 lattice is used. f eq i is the equilibrium distribution function. Its form plays a significant role in the physics reproduced by the LBM equation, and will be described in the following section. R( f neq i ) represents the off-equilibrium distribution reconstructed using the HRR collision operator. τ represents the nondimensional relaxation time τ = µ/(ρc 2 s ) + 1 2 . ψ i is a forcing term meant to correct errors that are introduced by the D3Q19 lattice (the O(Ma 3 ) error), as well as other errors introduced by f eq i .

B. Improved Density-based Isotropic Equilibrium

Eq. 1 is meant to reproduce the physics of the fully compressible Navier-Stokes equations obeying the ideal gas law p dim = ρ dim RT dim . It has been shown that there are many different choices of f eq i and ψ i that nominally produce the same macroscopic behavior, albeit with different stability properties. Initial models [START_REF] Feng | Hybrid recursive regularized thermal lattice boltzmann model for high subsonic compressible flows[END_REF][START_REF] Guo | An efficient lattice boltzmann method for compressible aerodynamics on d3q19 lattice[END_REF][START_REF] Renard | Improved compressible hybrid lattice boltzmann method on standard lattice for subsonic and supersonic flows[END_REF] used the standard density-based LBM equilibrium distribution functions, with temperature terms added to the second and third order moments, and with correction terms used to produce the correct viscous stress tensor. These models required significant amounts of additional dissipation in the form of a low value of the HRR parameter σ to stabilize the calculations. A pressure-based model [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] modified the zeroth order moment into a pressure-based equilibrium such that ∑ i f eq i = ρθ , greatly improving stability, but requiring a correction step to properly solve the mass conservation equation. Later, an improved model [START_REF] Guo | Improved standard thermal lattice boltzmann model with hybrid recursive regularization for commpressible laminar and turbulent flows[END_REF] was developed, which modified the calculation of the zeroth order moment for increased stability while retaining a density-based framework. Finally, Farag et al. [START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] performed an in-depth study showing that the pressure-based and density-based models can be unified into one single formalism. Here, an equilibrium distribution very similar to the improved-density based version of Farag et al. [START_REF] Farag | A unified hybrid lattice-boltzmann method for compressible flows: Bridging between pressure-based and density-based methods[END_REF] is used, although with certain modifications. An adjustment to the second order terms, developed by Bauer et al. [START_REF] Bauer | Letter to the editor: Truncation errors of the d3q19 lattice model for the lattice boltzmann method[END_REF] , is used to produce a greater level of isotropy.

As is typical for LBM approaches, Gauss-Hermite polynomials are used:

H (2) iαβ = c iα c iβ -c 2 s δ αβ , H (3)
iαβ γ = c iα c iβ c iγ -c 2 s c iα δ β γ + c iβ δ γα + c iγ δ αβ , (2)
The equilibrium distribution function is here split into four parts, corresponding to the zeroth to third order moments:

f eq i = ω i ρ( f eq,(0) i + f eq,(1) i + f eq,(2) i + f eq, (3) i
),

where ω i is the weight associated with each discrete velocity c iα . The first part, which contains the thermal information used to couple the LBM solver with the ideal gas law, is defined as:

f eq,0 i = 1+d i , where d i = ω 0 -1 ω 0 (θ -1), c i = (0, 0, 0) θ -1, else. (4)
The second part is unchanged from classic LBM approaches, and is equal to:

f eq,1 i = c iα c 2 s u α . (5)
The third part is based on the isotropic improvements of Bauer et al. [START_REF] Bauer | Letter to the editor: Truncation errors of the d3q19 lattice model for the lattice boltzmann method[END_REF] . It allows for the correct retrieval of certain fourth order moments. It is equal to:

f eq,2 i =          -u α u α , c i = (0, 0, 0) -3u α u α + 6(c iα u α ) 2 , c i ∈ {(±1, 0, 0), (0, ±1, 0), (0, 0, ±1)} -3 2 c 2 iα u 2 α + 9 2 (c iα u α ) 2 , else. (6
) Lastly, the fourth part is equal to:

f eq,3 i = 1 6c 2 s 3(H (3) ixxy + H (3) iyzz )(u x u x u y + u y u z u z )+ (H (3)
ixxy -H

(3) iyzz )(u x u x u y -u y u z u z )+ 3(H (3) ixzz + H (3) ixyy )(u x u z u z + u x u y u y )+ (H (3) ixzz -H (3) ixyy )(u x u z u z -u x u y u y )+ 3(H (3) iyyz + H (3) ixxz )(u y u y u z + u x u x u z )+ (H (3) iyyz -H (3) ixxz )(u y u y u z -u x u x u z ) . ( 7
)
The choice of the form of f eq i , together with the lack of discrete velocities in the D3Q19 lattice, leads to a number of errors in the viscous stress tensor, among them the well-known O(Ma 3 ) error. These are corrected via the forcing term ψ i in Eq. 1, equivalent to:

ψ i = -ω i H (2) iαβ 2c 4 s Ψ αβ , (8)
defined such that

Ψ α β = c 2 s u α ∂ (ρ(1 -θ )) ∂ x β + c 2 s u β ∂ (ρ(1 -θ )) ∂ x α + 2 3 δ αβ ρc 2 s ∂ u γ ∂ x γ -δ αβ c 2 s ∂ ρ(1 -θ ) ∂t + ∂ Err αβ γ ∂ x γ . (9)
The terms in Eq. 9 are discretized using a standard secondorder centered scheme. It can be shown through a Chapman-Enskog expansion (For further details, see Appendix A) that Eq. 1, together with the equilibrium distribution in Eq. 3, and the correction terms of Eq. 9 yield the compressible Navier-Stokes Equations:

∂ ρ ∂t + ∂ ρu α ∂ x α = 0, (10)
∂ ρu α ∂t + ∂ (ρu α u β + pδ αβ ) ∂ x β = ∂ ∂ x β τ αβ . ( 11
)
Where τ αβ is the viscous stress tensor:

τ αβ = µ( ∂ u α ∂ x β + ∂ u β ∂ x α - 2 3 ∂ u γ ∂ x γ ). (12)
The dynamic viscosity µ(θ ) varies with temperature according to Sutherland's law of viscosity.

C. Hybrid Entropy Equation

Because Eq. 10 and 11 do not calculate the evolution of energy, they must be coupled with an additional conservation equation. The approach used here is to use entropy s = c v ln( θ ρ γ-1 ), where γ is the ratio of specific heats. The equation for entropy is defined as:

∂ s ∂t + u α ∂ s ∂ x α = 1 ρθ ∂ ∂ x α (-q α ) + 1 ρθ τ αβ ∂ u α ∂ x β , (13)
Where the heat flux is defined as

q α = -λ f ∂ θ ∂ x α with the ther- mal conductivity λ f = c p µ
Pr . c p is the constant pressure specific heat capacity and Pr is the Prandtl number. The entropy equation is discretized using finite differences. The convective flux is computed using the MUSCL-Hancock method, identical to the one described by Farag et al. [START_REF] Farag | A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows[END_REF] , when a large stencil is available, while thermal conduction and viscous heat generation are computed using standard second order centered finite differences. The equation is advanced in time using a first order Euler approximation.

D. Subgrid Scale Model

The unresolved subgrid scales (SGS) are modeled using the Boussinesq approximation, thus modifying the form of the relaxation parameter to include SGS viscosity such that τ = µ+µ sgs ρc 2 s + 1 2 . The SGS viscosity is calculated using the Shear-Improved Smagorinsky Model (SISM), developed by Lévêque et al. [START_REF] Lévêque | Shear-improved smagorinsky model for large-eddy simulation of wall-bounded turbulent flows[END_REF] . This model is similar to the standard Smagorinsky Model, but seeks to improve its overlydissipative nature near walls by subtracting the magnitude of the mean shear from the instantaneous strain-rate tensor when calculating the eddy-viscosity. The mean shear is calculated using weighted exponential smoothing. The entropy equation is closed using a subgrid heat flux q α,SGS = -λ SGS

∂ θ ∂ x α , with λ SGS = c p µ SGS
Pr SGS and a constant SGS Prandtl number set to Pr SGS = 0.90.

E. Boundary Treatment in LBM

In the Lattice Boltzmann method, fluid nodes whose velocity distributions intersect a simulation domain boundary (henceforth called boundary nodes) require special treatment, as they lack at least one incoming distribution function. This becomes particularly important when complex geometries are involved, as the boundary nodes may be located at an arbitrary distance from the boundaries. Many different methods can be found in the literature, including partial reconstruction methods [START_REF] Dorschner | Grad's approximatino for moving and stationary walls in entropic lattice boltzmann simulations[END_REF] and interpolated bounce-back schemes [START_REF] Bouzidi | Momentum transfer of a boltzmann-lattice fluid with boundaries[END_REF] . Here, the full-reconstruction method is used.

First, macroscopic variables ρ, u, θ are set at boundary nodes using an interpolation procedure, which makes use of fictitious nodes located in the fluid domain whose macroscopic information is calculated using an inverse distanceweighing (IDW) algorithm. This allows for the equilibrium distribution function to be fully defined. Subsequently, finite difference approximations are used to evaluate the shear-stress tensor, which yields the off-equilibrium distribution functions, and the forcing terms, which are used to correct the Mach error. For inlets and outlets, the Local One Dimensional Inviscid (LODI) relations are used to apply characteristic boundary relations. These are calculated using finite differences. For more details, the reader is invited to consult the paper of Feng et al. [START_REF] Feng | Solid wall and open bound-ary conditions in hybrid recursive regularized lattice boltzmann method for compressible flows[END_REF] .

F. Transition Algorithms

Since complex turbulent flows involve a wide variety of length-scales, mesh elements of non-uniform size are necessary so that finer mesh elements can be used in zones of interest. For most LBM approaches, octree mesh structures are used, as is the case here. This has significant advantages: nodes in a coarse zone are calculated half as often as nodes in a fine zone, leading to important increases in computational speed. However, the interfaces between different refinement zones present certain difficulties. There are numerous methods of treating the communication between two zones of different refinement. Here, the Direct Coupling (DC) algorithm of Astoul et al. is used [START_REF] Astoul | Lattice boltzmann method for computational aeroacoustics on non-uniform meshes: A direct grid coupling approach[END_REF] .

This method calculates the values at co-located coarse and fine nodes by taking information from both the coarse and fine part of the grid and makes use of LBM collision invariance relations, together with a Newton-Raphson algorithm, to produce distribution functions that greatly reduce interface discontinuity and spurious noise. Since in the current study, f eq i is also a function of θ , this method must be adapted to include temperature. Entropy s is calculated prior to the application of the DC algorithm, using a coarse stencil. For asynchronous time steps, a temporal interpolation is used. This entropy value is then used in the Newton agorithm of the DC method used to calculate the distribution functions. For fine nodes that are not co-located with coarse nodes, spatial interpolations are used.

G. Anisotropic Turbulent Inflow

Proper turbulence injection conditions are crucial to the simulation of impinging jet flows. For large-eddy simulations, both the correct mean velocity profile together with the instantaneous injection of realistic turbulent fluctuations are required. Turbulence injection remains a highly active area of research, involving numerous different methods that have been covered by numerous literature reviews [START_REF] Wu | Inflow turbulence generation methods[END_REF][START_REF] Tabor | Inlet conditions for large eddy simulation: A review[END_REF][START_REF] Dhamankar | Overview of turbulent inflow boundary conditions for large-eddy simulations[END_REF] . Methods include precursor simulations, internal mapping, random Fourier sums, digital filtering, synthetic eddy methods, and volume forcing methods. In the literature of turbulence resolving impinging jet simulations, Hadžiabdić and Hanjalić 8 used a precursor simulation, while Uddin et al. [START_REF] Uddin | LES simulations of an impinging jet: On the origin of the second peak in the Nusselt number distribution[END_REF] used a dig-ital filter. Both Dairay et al. (whose simulation was meant to simulate a converging nozzle) [START_REF] Dairay | Direct numerical simulation of a turbulent jet impinging on a heated wall[END_REF] and Aillaud et al. (whose simulation represented a long pipe injection) 10 used a random Fourier sum that generated homogeneous isotropic turbulence. Grenson and Deniau used a combination of a random Fourier producing homogeneous isotropic turbulence, as well as a vortex ring generator at the walls meant to produce boundary layer turbulence [START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] . For the present simulation, an approach that uses a random Fourier sum, along with a modification to produce non-homogeneous, anisotropic flow, will be used.

The mean velocity profile, along with the Reynolds stress tensor components, are taken from the pipe flow DNS at Re b = 24 580 of Wu et al. [START_REF] Wu | Direct numerical simulation of a 30r long turbulent pipe flow at r+ = 685: large and very large-scale motions[END_REF] . The simulation data, which can be found online in wall units, is re-scaled with friction velocity u τ = 2.01 m s -1 ). The turbulent fluctuations are generated with the following steps similar to the one described by Shur et al. [START_REF] Shur | Synthetic turbulence generators for rans-les interfaces in zonal simulations of aerodynamic and aeroacoustic problems[END_REF] :

• Homogeneous isotropic turbulent fluctuations based on the Von-Karman Pao spectrum are generated. This approach is the similar to the one described by Bailly and Juvé [START_REF] Bailly | A stochastic approach to compute subsonic noise using linearized euler's equations[END_REF] . However, the amplitude of each mode ûn is normalized such that the total turbulent kinetic energy is equal to unity. This yields the homogeneous isotropic fluctuations defined as:

v HIT (x,t) = 2 3 2 N ∑ n=1 ûn cos[k n •(x-tu c )+φ HIT,n +ω HIT,n t]σ n . (14)
Where k n , is the wavevector whose magnitude, u c is the bulk convective velocity (in this case (0, -36.0, 0) m s -1 , φ HIT,n is the phase, and ω HIT,n is a frequency dependent on k n , and σ n is a unit vector defining the direction of the fluctuation.

• The fluctuations are then scaled using a i j (x), which is the Cholesky decomposition of the Reynolds stress tensor R i j (x) such that:

    √ R xx 0 0 R xy /a xx R yy -a 2 yx 0 R xz /a xx (R yz -a yx a zx )/a yy R zz -a 2 zx -a 2 zy     .
The result is that the fluctuations, which are now nonhomogeneous and anisotropic, are equivalent to:

u i (x,t) = a i j (x)v HIT j (x,t). (15)
It is important to note that although Eq. 14 theoretically produces fluctuations with zero divergence

∂ u α ∂ x α = 0 if k n • σ n = 0,
this is no longer the case when the anisotropic scaling is employed. This produced significant amounts of non-physical pressure fluctuations in the pipe. These are attenuated using a sponge zone near the pipe inlet, relaxing the values of density to the mean using a moving average, while the velocity vector is allowed to evolve freely. The implementation of the sponge zone is described by Feng, et al. [START_REF] Feng | Solid wall and open bound-ary conditions in hybrid recursive regularized lattice boltzmann method for compressible flows[END_REF] .

The integral length scale used to generate the Von Karman-Pao turbulence spectrum is L i = 0.4D. In the current approach, although the intensity of the fluctuations are nonhomogeneous and anisotropic, the length scales remain homogeneous and isotropic. Because the inflow data of Wu et al. is in cylindrical coordinates, and because the meshes used in the pipe-flow DNS and the current jet-flow simulation are not the same, both a transformation to Cartesian coordinates, together with an interpolation to match the grid nodes at the inlet, are required in order to assign each boundary node the correct data. These are performed at the beginning of the simulation using basic coordinate transformations, and a 1-dimensional inverse distance weighing (IDW) scheme based on the radial position, respectively.

III. SETUP A. Flow Configuration

The flow configuration, as shown in Fig. 1, is composed of a cylindrical pipe of inner diameter D = 0.01 m and of height 5D embedded into a larger cylindrical domain of diameter 24D and total height 7D. It is an unconfined, axisymmetric, isothermal, subsonic jet, with an impingement height H = 2D. The origin of the coordinate system is placed at the very center of the impacted plate, at the stagnation point. The axis of the flow is in the negative y-direction. The dynamic viscosity

µ at T = 300K is set to ≈ 1.84 • 10 -5 kg • m -1 • s -1 .

B. Boundary Conditions

At the pipe inflow, a mean velocity profile in the negative y-direction is imposed such that it recovers the shape of the DNS of Wu et al [START_REF] Wu | Direct numerical simulation of a 30r long turbulent pipe flow at r+ = 685: large and very large-scale motions[END_REF] with a bulk speed U b = 36 m s -1 . Resolved non-homogeneous anisotropic turbulent fluctuations are superposed onto the mean velocity profile. The temperature of the jet is T j = 300 K, yielding a bulk Mach number Ma b equal to ≈ 0.1 and a bulk Reynolds number Re b = 23 000. The walls of the pipe are isothermal, with a temperature T pipe = 300K. As the jet is unconfined, a choice must be made regarding the boundary condition used for the top boundary. Here, a coflow velocity with an intensity of 2% of the bulk velocity is used, which is similar to the axial velocity far away from the jet recorded by Tummers, et al. [START_REF] Tummers | Turbulent flow in the near field of a round impinging jet[END_REF] . Both inflow velocities are at a uniform temperature T ∞ = 300K. The impingement plate is a no-slip isothermal wall condition at T w = 330 K. Finally, the sides of the domain are given a pressure outlet condition at p ∞ = 101325.0 Pa.

C. LBM Model Parameters

The free HRR parameter σ is set to 0.99. T re f , which determines the timestep through the relation ∆t

= ∆xc s √ RT re f
, is set to

D. Meshing Strategy

LBM approaches typically use Cartesian mesh elements as an inherent feature of the method, which is here coupled with the above-mentioned octree refinement algorithms allowing for different levels of fidelity depending on the region of interest. The domain is split into separate resolution domains (RD), with RD1 using the finest grid spacing, RD2 using a grid spacing that is 2x greater, and so on.

Although the Cartesian octree approach allows for significantly greater rapidity in meshing for the user, it also presents an important constraint. The use of cubic mesh elements prohibits the use of grid stretching, which is often used in zones where the user as a priori knowledge of significant anisotropic non-homogeneity, such as in boundary and shear layers. In the literature for impinging jet LES, grid stretching ratios on the order of 10 at the impacted plate are common [START_REF] Hadžiabdić | Vortical structures and heat transfer in a round impinging jet[END_REF][START_REF] Uddin | LES simulations of an impinging jet: On the origin of the second peak in the Nusselt number distribution[END_REF][START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF][START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] . Thus, reaching a given wall-normal mesh spacing using a Cartesian mesh can require far more cells to cover a given volume. The result is that the present simulation uses meshes that are significantly larger than the ones used in comparable impinging jet LES using Navier-Stokes approaches.

Three meshes are treated in this study. Mesh 1 sets the grid spacing for RD1 to ∆x/D = 0.005, leading to a pipe wall ∆r + = 6.4, and an impact plate ∆y + = 8.45. Mesh 2 adds an additional zone of refinement, with the grid spacing for RD1 set to ∆x/D = 0.0025. RD1 is placed at the pipe walls, in the initial free jet shear layer, and at the impact plate. The remaining refinement zones are identical to those of Mesh 1. This additional zone leads to a division of the pipe wall and impact plate normal wall distance by two. Finally, Mesh 3 is identical to Mesh 2, but with a wholesale refinement of the entire domain, with the grid spacing for RD1 set to ∆x/D = 0.002. Since this is an octree mesh, the grid spacing in every single direction of every resolution domain is reduced by 20%, nearly doubling the amount of points.

The properties of the three meshes used are summarized in Table I fine" nodes: the points in RD1 are counted as 1 equivalent fine point, those of RD2 are counted as 1/2 of an equivalent fine point, and those of RD3 are counted as 1/4 of an equivalent fine point, etc. The mesh is Cartesian, so that locally, ∆x = ∆y = ∆z. The result is that at the impact plate, the wallnormal distance is invariant relative to φ , but the streamwise and spanwise spacings are not. In the pipe, the streamwise spacing is invariant relative to φ , while the wall-normal and spanwise spacings are not.

Figure 2 shows a view Mesh 2, together with instantaneous vorticity magnitude, at the principal zone of interest in the main section of the domain, as well as the pipe exit. The zones where RD1 is employed are shown with the zoomed in boxes. The radial locations sampled to validate against the experimental data are shown with dashed lines. In the pipe, three levels of mesh resolution are used in order to maximize the mesh resolution at the pipe walls. The pipe flow core is meshed using RD3, whereas the pipe boundary layer is meshed using RD2 and RD1 As the pipe boundary layer transitions to become the free-jet shear layer, RD1 is extended 0.5D downstream of the pipe exit, while RD2 is extended from the shear layer all the way to the impact plate. RD2 and RD1 cover the entire surface of the impact plate for R/D ≤ 3.

Figure 3 shows the azimuthally averaged value of the local wall-normal mesh size ∆y + along the surface of the impinged plate. The maximum value of is 8.45, 4.68, and 3.82 for Meshes 1, 2, and 3 respectively. These values are higher than the ones recommended by Choi and Moin [START_REF] Choi | Grid-point requirements for large-eddy simulation: Chapman's estimates revisited[END_REF] for a typical wallresolved LES boundary layer, although ∆y + for Mesh 3 is similar to the value found in the simulation of Aillaud et al. [START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF] . The streamwise and spanwise mesh sizes are, on the other hand, significantly smaller than the standard boundary-layer requirements. While the streamwise, spanwise, and wallnormal mesh requirements for an attached, standard boundary are well attested in the literature, these values are not well attested for impinging jet LES, and high degrees of grid stretching can lead to poor results [START_REF] Hadžiabdić | Vortical structures and heat transfer in a round impinging jet[END_REF][START_REF] Uddin | LES simulations of an impinging jet: On the origin of the second peak in the Nusselt number distribution[END_REF] .

E. Simulation Strategy

The total simulation time corresponds to roughly 40 = 0.65. f is the frequency of the passage of the primary structures, which will be discussed in Section V. Statistics were taken over a period of 25 τ conv or approximately 16 τ St . This is similar to the time used to construct statistics in other impinging jet LES [START_REF] Hadžiabdić | Vortical structures and heat transfer in a round impinging jet[END_REF][START_REF] Uddin | LES simulations of an impinging jet: On the origin of the second peak in the Nusselt number distribution[END_REF][START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF] . For the velocity and surface profiles, the statistics are azimuthally averaged.

IV. SIMULATION STATISTICS A. Pipe Velocity Profiles

Figure 4 shows the mean and RMS velocities, along with the Reynolds shear stress, at the pipe exit located at y = 2D, which is 5D downstream of the inlet. They are compared against the periodic DNS of Wu et al. [START_REF] Wu | Direct numerical simulation of a 30r long turbulent pipe flow at r+ = 685: large and very large-scale motions[END_REF] at Re = 24 580. The mean velocity profile shows excellent agreement with the reference DNS. The axial fluctuations also show good agreement, with only minor differences shown. More noticeable discrepancies, on the order of 20%, can be seen for the radial and azimuthal fluctuations, along with the Reynolds shear stress.

Although there is a significant increase in the total number of points between the first, second, and third meshes, the three meshes employed show largely similar results. The level of discretization in the jet core, where discrepancies with the DNS are more visible, remain similar for all three meshes (∆(x, y, z) + = 12.8 for Meshes 1 and 2, ∆(x, y, z) + = 10.24 for Mesh 3). Mesh 2 adds an an additional zone of refinement in the pipe boundary layer while conserving the same refinement level in the pipe core (since octree meshes are used, an additional refinement zone leads to a considerable increase in

the total number of mesh nodes), leading to a noticeable improvement in the Reynolds shear stress. This improvement is confirmed by Mesh 3, which has the same refinement zones as Mesh 2, but has a different baseline mesh size to Mesh 2 (all of the mesh nodes of Mesh 3 have length 4/5 the length of the equivalent nodes of Mesh 2, leading to another multiplication by 2 of the mesh elements).

Although the amplitude of the fluctuations are nonhomogeneous and anisotropic, and are imposed by reference data, the length scales are still generated using homogeneous isotropic assumptions. The DNS of Wu et al. [START_REF] Wu | Direct numerical simulation of a 30r long turbulent pipe flow at r+ = 685: large and very large-scale motions[END_REF] performed a detailed analysis of the streamwise and azimuthal spectra and found large degrees of variation at different flow locations. Any single choice of length-scales will therefore produce significant errors. Several tests were made by varying the lengthscales used to generate the Von Karman Spectrum. In this study, larger length-scales were ultimately chosen, as it was found that smaller length scales led to an unacceptable decay of the turbulence intensity further downstream. Furthermore, the relatively short entrance length of 5D used to maintain a reasonable simulation time also contributed to the discrepancies.

Overall, the agreement with the reference pipe flow can be considered fair. However, more advanced turbulence injection methods capable of further reducing the entrance length would be highly useful for future LBM simulations, as resolving attached, equilibrium turbulent boundary layers is particularly expensive for LBM approaches relative to body-fitted Navier-Stoke solvers, as this type of flow topology can be resolved in many Navier-Stokes approaches via highly extensive grid stretching.

B. Jet Velocity Profiles

The flow dynamics are validated on experimental data using laser Doppler anemometry (LDA) provided by Tummers et al. [START_REF] Tummers | Turbulent flow in the near field of a round impinging jet[END_REF] . Their experiment took detailed measurements at eight radial stations r/D = 0, 0.25, 0.5, 0.75, 1, 1.5, 2, and 2.5.

Figure 5 shows the flow dynamics at three stations r/D = 0, 0.25, and 0.5, corresponding to each row. The five columns represent mean radial velocity u r , mean axial velocity u y , Reynolds radial stress u 2 r , Reynolds axial stress u 2 y , and Reynolds shear stress u y u r , respectively.

At r/D = 0, directly in the core of the jet, the quantities are quite close to the experimental data. The turbulence intensities remain mostly constant from the nozzle exit until around y/D ≈ 0.4, indicating that this region is in the potential core of the free jet, as remarked upon by Hadžiabdić and Hanjalić [START_REF] Hadžiabdić | Vortical structures and heat transfer in a round impinging jet[END_REF] . Both the Reynolds axial and radial stresses are underestimated as the flow approaches the impingement wall for all meshes. This is likely due to the choice of length scales in the synthetic turbulence injection, combined with the relatively short pipe length. It is likely that turbulence in the center line was not able to adapt sufficiently quickly to fully developed pipe flow conditions. It is possible that the finer mesh produced smaller turbulent structures, leading to a slight underestimation of the axial Reynolds stress at r/D = 0. The fact that sometimes grid refinement does not lead automatically to improvement of the results in LES has already been reported in several configurations. This is mainly due to partial error cancellation between the numerical scheme and the subgrid model.

At r/D = 0.25, halfway between the jet core and the shear layer, the jet has deflected, transferring momentum along the wall in the radial direction. The agreement with the exper-iment is shown to be quite good, although there are slight discrepancies in the Reynolds stresses, notably u y u r for high values of y/D. Meshes 2 and 3, which both use RD3 for r/D = 0.25, y/D > 0.2, show similar degrees of discrepancy despite refinement from ∆(x, y, z)/D jet = 0.01 to ∆(x, y, z)/D jet = 0.008. A further verification was performed with a fourth mesh where RD2 (∆(x, y, z)/D jet = 0.005) was used in this zone (not shown here), with no noticeable change to the discrepancies for the Reynolds stress at y/D > 1. As shown in Fig. 4, all of the meshes slightly overestimate the Reynolds shear stress at r/D = 0.25 of the nozzle exit. This overestimation in the pipe may contribute to the continued discrepancies the free jet seen in Fig. 5.

The third row of Figure 5 shows the flow dynamics at r/D = 0.5, the position directly in the shear layer between the jet and the surrounding flow. The axial velocity downstream of the lip increases, as the shear layer thickens. The mean radial velocity along the wall has accelerated significantly from the previous station, doubling in magnitude. The strong gradients between the jet and the surrounding flow produces high intensities of turbulence in the free jet relative to the values seen in the jet core at r/D = 0 and r/D = 0.25. Although the overall results are good, they present some discrepancies, with a slight overestimation in the magnitude of the mean axial velocity, as well as overestimation of the all fluctuations in the free-jet. Mesh 1, with its coarse resolution of the pipe boundary layer and free jet shear layer, exhibits the largest discrepancies, showing a significant overestimation in the Reynolds stresses, while Mesh 2 and 3 display less marked discrepancies. The overestimation of the fluctuations have been observed in other LES, [START_REF] Hadžiabdić | Vortical structures and heat transfer in a round impinging jet[END_REF][START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] overestimation in the width of the shear layer. Figure 6 shows the development of the jet profiles further downstream, in the wall jet. At r/D = 0.75, the mean radial velocity continues to accelerate, now exceeding U b at its peak. All of the results for all of the meshes are shown to be in good agreement.

At the station r/D = 1, the wall jet has become noticeably thinner compared to r/D = 0.75, and has continued to accelerate, reaching a maximum value of 1.10U b for the experiment and for all of the meshes. The agreement with the experiments remains very good, with only a moderate overestimation for all three meshes seen for the Reynolds shear stress.

As the wall jet develops, by the time it reaches r/D = 1.5 it has begun to slightly decelerate, with a maximum radial velocity declining from about to about 1.02U b . At the same time, the wall jet has also become more turbulent. The magnitudes of the maximum axial, radial fluctuations, along with the Reynolds shear stress, see significant increases. Compared to the experimental results, the discrepancies remain relatively small, with only the Reynolds shear stress once again slightly overestimated.

At r/D = 2, the wall jet has decelerated in earnest, reaching a maximum radial velocity just below 0.8U b . At the same time, the turbulent intensity of the jet has continued to increase, reaching a maximum (relative to the other stations sampled) for the axial and radial fluctuations, as well as in the Reynolds shear stress. Notably, this radial station is also approximately the location of the secondary peak of the Nusselt number, which will be discussed in Section IV C.

Lastly, at r/D = 2.5, the wall jet has both thickened and decelerated. This time, the maximal magnitude of the turbulent fluctuations has decreased relative to the previous station.

When taking into account the level of uncertainty in the experimental Reynold stress statistics is between 5 -8% [START_REF] Tummers | Turbulent flow in the near field of a round impinging jet[END_REF] , the overall agreement of the flow statistics with the experiments is very good, particularly in the wall jet. Even the coarsest mesh, which has rather poor radial resolution at the pipe walls (∆r + as high as 9.05) and jet shear layer, was able to produce largely satisfactory results for almost every quantity at almost every station. It can be hypothesized that the excellent azimuthal mesh resolution, as well as the excellent axial (in the free-jet) and radial (in the wall-jet) resolution plays a role here.

C. Heat Transfer and Shear Stress at the Impinged Plate

The temporal and azimuthally averaged Nusselt number, which signifies the ratio of convective to conductive heat transfer, is defined as:

Nu t,az (r) = h t,az (r)D λ f , (16)
where h(r, φ ,t) is the heat transfer coefficient defined as:

h(r, θ ,t) = q w (r, θ ,t) T re f -T w , (17)
and where q w is the local wall heat transfer that takes into account the fluid thermal conductivity λ f , which varies according to Sutherland's law, as well as SGS thermal conductivity λ SGS :

q w (r, φ ,t) = (λ f (r, φ ,t) + λ sgs (r, φ ,t)) ∂ T ∂ n (r, φ ,t) (18)
The choice of λ f for the denominator in Eq. 16 can vary when using compressible solvers with variable viscosity. Otero-Pérez and Sandberg [START_REF] Otero-Pérez | Compressibility and variable inertia effects on heat transfer in turbulent impinging jets[END_REF] , in their study of Mach number and temperature gradient effects on impinging jet heat transfer, noted that the variation of λ f has rarely been evoked in the literature (as this test case has mostly been done for incompressible flows). They state that either λ f (T ) or λ f constant can be used in the calculation of the Nusselt number. In this particular case, the low Mach number and temperature gradient lead to only minor differences. Here, the denominator in Eq. 16 will be taken as λ f ,∞ , or the thermal conductivity of air at T ∞ , as is done in the study of Aillaud et al. [START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF] . Figure 7a shows the temporal and azimuthally averaged Nusselt number for all three meshes along the impingement surface up to r/D = 3, along with the experimental results found by Baughn and Shimizu [START_REF] Baughn | Heat transfer measurements from a surface with uniform heat flux and an impinging jet[END_REF] at Re = 23 750, as well as those found by Fénot et al. [START_REF] Fenot | Local heat transfer due to several configurations of circular air jets impinging on a flat plate with and without semi-confinement[END_REF] at Re = 23 000. Impinging jet Nusselt number results exhibit a wide degree of variation across different experiments in the literature, and here the LES results are found to be largely in between the two results, with a slight overestimation in the stagnation zone. The Nusselt number in the vicinity of the jet impact for all three meshes is close to those of Baughn and Shimizu, but the decline along the wall jet region is slightly more marked, reaching levels slightly below those of Fénot et al. at around r/D = 1.25. Mesh 1 exhibits a noticeably weaker second peak, but it still attains a level of heat transfer similar to the one observed by Fénot et al. Mesh 2 and Mesh 3 show a strong second peak reaching the level of Baughn and Shimizu. Overall, taking into account the large variations in the Nusselt number seen in experiments, the agreement between the simulations and the experiments can be considered good. It is interesting to note that Mesh 1, despite the wall-normal mesh resolution that falls far below the levels typically recommended (∆y + = 8.45), is still able to recover a rather accurate secondary peak. It is likely that the high radial and azimuthal resolution plays a significant role here, as it has been hypothesized that the phenomena that produce the secondary peak of the Nusselt number show significant stream and spanwise non-homogeneity [START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] .

Figure 7b shows the normalized wall shear stress compared with the results of Tummers et al. [START_REF] Tummers | Turbulent flow in the near field of a round impinging jet[END_REF] . The results for Meshes 2 and 3 show good agreement, although there are important discrepancies between r/D = 0.5 and r/D = 1, though the small number of experimental points make it difficult to evaluate the degree of discrepancy in this zone. In contrast, Mesh 1 largely underestimates the near-wall shear stress. The LES of Aillaud et al. [START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF] similarly found that the wall shear stress was more sensitive to mesh resolution than the Nusselt number.

The flow-field statistics and surface quantities are overall in good agreement with the experimental literature for for Meshes 2 and 3. With the exception of the mean wall shear

, (• Tummers, et al. [START_REF] Tummers | Turbulent flow in the near field of a round impinging jet[END_REF] ), ( Fénot, et al. [START_REF] Fenot | Local heat transfer due to several configurations of circular air jets impinging on a flat plate with and without semi-confinement[END_REF] ), ( Baugh and Shimizu 40 ) stress, even Mesh 1 is capable of producing very satisfactory statistical results at a fraction of the cost. We will next perform an examination of the unsteady flow structures, in order to determine whether the LBM approach used here is capable of producing the same qualitative flow features seen in other LES in the literature. This investigation is performed with Mesh 2.

V. UNSTEADY FLOW CHARACTERISTICS

A. Primary Structures

The presence of azimuthally coherent vortical structures in the free-jet shear layer is well attested for jets issuing from converging nozzles. However, these structures are not always seen for long pipe jets, both in experiments and in unsteady simulations. Mi et al. [START_REF] Mi | Influence of jet exit conditions on the passive scalar field of an axisymmetric free jet[END_REF] compared the flow-field downstream of the nozzle for both a converging nozzle injection and a long pipe injection at Re = 16 000, and found that there were no periodic azimuthally coherent structures in the early stages of the free jet for the long pipe injection jet. By contrast, the experiment of Grenson et al. [START_REF] Grenson | Investigation of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number: An extensive experimental approach[END_REF] at Re = 60 000 revealed the formation of large scale structures immediately downstream of the nozzle. The LES of Hadžiabdić and Hanjalić [START_REF] Hadžiabdić | Vortical structures and heat transfer in a round impinging jet[END_REF] showed structures with weak azimuthal coherence immediately downstream of the nozzle, while the simulation of Uddin et al. [START_REF] Uddin | LES simulations of an impinging jet: On the origin of the second peak in the Nusselt number distribution[END_REF] did not show any azimuthally coherent structures in the free jet at all. By contrast, the simulations of Aillaud et al. [START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF] and of Grenson and Deniau 11 produced highly coherent structures in the free jet.

Figure 8 shows an instantaneous snapshot of the pressure isosurfaces, colored by distance from the impingement plate. The structures in the free-jet, together with the vortex ring at the wall jet, display a moderate degree of azimuthal coherence. In order to identify the frequency of azimuthally coherent structures, probes are placed in the domain, with an azimuthal resolution of ∆φ = π/4, and an axial resolution of ∆y/D = 0.1. The outputs are taken at a frequency equivalent to St = 867.5. To capture the azimuthally coherent structures, the probes are azimuthally averaged.

Figure 9 shows a map of the dominant frequencies of pressure along the y-axis at two radial locations r/D = 0.5 and r/D = 1.25. It can be seen that the frequency of the primary structures in the jet shear layer and wall jet is equal to approximately St = 0.65.

Figure 10 shows a view "beneath" the impinged plate. The primary structures, represented by the isosurfaces of static pressure, are colored in blue and are located at approximately r/D = 1.5. Iso-surfaces of the Q-criterion are colored in green. The zones of flow separation at the impact plate, defined as regions where the radial wall shear stress τ w,r < 0, are colored in black. The effect of the impact of the primary structures near the wall has been subject to some variation within the literature. The DNS of Dairay, et al. [START_REF] Dairay | Direct numerical simulation of a turbulent jet impinging on a heated wall[END_REF] (using a converging nozzle injection) and the LES of Aillaud, et al. [START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF] showed a secondary ring vortex, produced by the passage of the primary vortex, that showed a high degree of azimuthal coherence. By contrast, the LES of Uddin et al. [START_REF] Uddin | LES simulations of an impinging jet: On the origin of the second peak in the Nusselt number distribution[END_REF] did not show any near-wall separation, whereas the LES of Grenson and Deniau [START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] showed unsteady separation, but in the form of isolated spots rather than as a coherent ring. The experimental study of Yadav and Agrawal [START_REF] Yadav | Effect of vortical structures on velocity and turbulent fields in the near region of an impinging turbulent jet[END_REF] observed secondary vortices for Re = 3 050 and 5 700, but not for Re = 10 000. In the present LES, the nearwall flow separation behavior is similar to the one shown in the study of Grenson,et al. 11 . While a roughly coherent ring of structures of azimuthal vorticity can be seen near the primary structures via the Q-criterion, only small patches of the boundary layer show flow separation.

F sep (r,t) = 1 2π 2π 0 min[τ w,r (r, φ ,t), 0]rdφ (20)
Figure 11a shows the evolution of the azimuthally averaged Nusselt number over time and space. At r/D = 1.25, the passage of the primary structure, separated by the characteristic time τ St , is associated with a dip in the Nusselt number, forming a front of weak thermal exchanges similar to the one described by Dairay et al. [START_REF] Dairay | Direct numerical simulation of a turbulent jet impinging on a heated wall[END_REF] and Aillaud et al. [START_REF] Aillaud | Secondary peak in the Nusselt number distribution of impinging jet flows: A phenomenological analysis[END_REF] . However, in contrast to these simulations and the one of Grenson and Deniau [START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] , no distinguishable periodic mode of strong thermal exchanges at the radial location of the secondary peak can be seen.

Examining Figure 11b, the radial shear stress at r/D < 0.1 is nearly zero, corresponding to the stagnation zone of the jet, before increasing precipitously. For r/D < 1, the azimuthal average is largely steady. Beyond this point, periodic declines in the shear stress can be seen with the the same characteristic time τ St , linking the passage of the primary structures with lower radial shear stress. This phenomenon can be linked to the unsteady separation. It can be seen from Figure 11c that the flow separation occurs regularly with roughly the frequency of the passage of the primary structures. Since the zone of flow separation does not form a coherent ring, as shown in Figure 10, this does not cause the azimuthally averaged instantaneous radial shear stress to go below zero. The flow separation begins at roughly r/D = 1.2, in agreement with the experimental results observed by Tummers et al. [START_REF] Tummers | Turbulent flow in the near field of a round impinging jet[END_REF] . Similarly, the location where the separation criteria is at most intense, at r/D ≈ 1.5 also corresponds to the local minimum of the Nusselt number, as was also observed by Grenson and Deniau [START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] .

B. Near-Wall Aerothermal Dynamics

In order to evaluate the ability of the LBM solver to reproduce the fine near-wall physics, an examination of the instantaneous near-wall flow phenomena is performed. The twodimensional shear-stress vector τ w can be considered as projections of the near-wall velocity onto the surface of the plate. Thus, a criterion of near-wall 2-dimensional divergence can be defined as: Since this flow is considered to be incompressible, positive 2D divergence in the near-wall r -φ plane indicates flow from above "downwashing" onto the plane.

div 2D (τ w ) = ∂ τ w ∂ r + ∂ τ w ∂ rφ (22)
Figure 12 shows the instantaneous distribution of the Nusselt number for one quarter of the impacted plate 0 < φ < π/2. The top image highlights zones of flow separation through contour lines of τ w,r < 0, while the bottom image highlights zones of downwash through contour lines of div 2D (τ w ). It can be seen that there are far more thermal spots than there are spots of local flow separation, indicating that the principal increases in local heat transfer do not directly originate from unsteady separation. By contrast, almost all thermal spots correspond to the contour lines of 2D divergence. Since the divergence is necessarily compensated by fluid coming from above, and since this fluid is generally colder than the fluid near the impinged plate, greater 2D divergence leads to a greater value of the Nusselt number. This is strongly in line with the results reported by Grenson and Deniau [START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] , who denoted that all thermal spots were associated with the divergence of the friction lines.

To examine the near-wall flow conditions leading to the divergence of the 2D shear-stress vectors at the wall, and the subsequent increase in thermal exchanges, it is instructive to inspect the 3D near-wall flow field.

Figure 13a shows the topology of an "impingement" type thermal spot. On the left, the colormap shows that Nusselt number at this spot high, and the friction lines in the vicinity of the spot, together with the colormap of the radial shear stress as seen in Figure 13a(b), indicate a local flow reversal near the wall. As previously shown in Figure 10, the flow separation is only a "patch", rather than forming a ring. In fact, the thermal spot to the right and above the "impingement" spot shows streamlines diverging in the azimuthal, not radial direction. Finally, the 3-D view on the right shows the ry plane intersecting the thermal spot. Flow separation of the boundary layer ejects hot fluid upward while bringing cold fluid towards the plate, thinning the temperature boundary layer and enhancing thermal exchanges.

Figure 13b shows the formation of a near-wall streamwise flow structure leading to a thermal spot, with iso-surfaces of the Q-criterion colored in green. On the left, the vorticity is first oriented in the azimuthal direction, with boundary layer separation. The center image shows the structure breaking down and stretching in the radial direction, forming a streamwise structure. The friction lines underneath the streamwise structure diverge in the azimuthal direction. On the right, the intersection with the φy plane shows the up and down-wash induced by the streamwise structure. This phenomenon is similar to the one seen for the "impingement" structure, but with the mixing occurring in a plane perpendicular to the flow direction. Not all streamwise structures are formed by the breakdown of an impingement structure, and further analysis is required to better characterize their formation.

The hybrid compressible LBM approach is able to produce fine near-wall flow dynamics that are similar to the ones observed by Grenson and Deniau 11 , despite lacking the same de- An LES of an axisymmetric thermal impinging jet issuing from a pipe was performed using a compressible Lattice Boltzmann Method. Using an improved-density based thermal equilibrium, an entropy equation resolved by finite differences, cut-cell full reconstruction boundary treatment, direct-coupling transition algorithms, and non-homogeneous anisotropic turbulence injection, very good jet flow-field statistics and surface quantities were obtained for the two finer grids, while good results, with the exception of the wall shear stress, were obtained for the coarser one.

The simulations performed in this article also made use of anisotropic, non-homogeneous turbulence injection, producing satisfactory results for a reasonably short entrance length. However, the non-homogeneous anisotropic fluctuations produced non-physical pressure fluctuations that had to be damped via a sponge zone, and better results might be obtained with different approaches to turbulence injection that do not rely on homogeneous isotropic length-scales.

The simulation produced a free-jet whose vortex structures displayed a mild degree of azimuthal coherence. The vortex structures in the wall-jet led to local unsteady flow separation at the wall, although generally these zones of flow separation were relatively isolated. The thermal spots at the wall involving flow separation were significantly less common than thermal spots that did not involve flow separation. The near-wall flow structures leading to the thermal spots were similar to the ones observed by Grenson and Deniau [START_REF] Grenson | Large-eddy simualtion of an impinging heated jet for a small nozzle-to-plate distance and high reynolds number[END_REF] .

Although it can be seen here that an LBM approach can be a valid choice for this test case, it is by no means established here that it is a conclusively superior choice. The use of Cartesian meshes means that in regions of high flow anisotropy (in this case: the pipe boundary layer, the free jet shear layer, the impact zone, and the wall-jet), wall-resolved LBM approaches require significantly greater numbers of mesh elements compared to traditional Navier-Stokes approaches that can make use of grid stretching. However, this significantly larger mesh is not without its advantages. Attaining reasonably high wallnormal or shear-layer resolution using this approach also results in excellent mesh resolution in the streamwise and spanwise direction, which may be a factor in the high quality of the flow-field results seen in this paper.

The hybrid compressible LBM approach has thus shown that it is sufficiently mature to perform nearly incompressible test cases of moderate complexity. An interesting subject of further investigation would be the simulation of compressible impinging jets involving higher Mach numbers and/or temperature differences.

Appendix A: The Chapman-Enskog development Let us start from the discrete velocity Boltzmann BGK equation. As was the case in Section II, the variables are here presented in non-dimensional form.

∂ f i ∂t + c iα ∂ f i ∂ x α = - 1 τ ( f i -f eq i ) + ψ i (A1)
The moments of f eq i for the D3Q19r Lattice (which fails to properly recover the third order equilibrium moment), are defined as: Equation A1 is expanded to zeroth and first order according to a parameter ε related to the Knudsen number, such that the temporal operator is defined as ∂ ∂t ≈ ∂ ∂t (0) + ε ∂ ∂t (1) . The equations at order zero and one are, respectively:

∂ f (0) i ∂t (0) + c iα ∂ f (0) i ∂ x α = - 1 τ ( f (1)
i ) + ψ i (A6) ∂ f (1) i ∂t (0) + ∂ f (0) i ∂t (1) + c iα ∂ f (1) i ∂ x α = - 1 τ ( f (2)
i ) (A7)

The zeroth, first, and second order moments of Eq. A6 are: The zeroth and first order moments of Eq. A7 are:

∂ ρ ∂t [START_REF] Zuckerman | Jet impingement heat transfer: Physics, correlations, and numerical modeling[END_REF] = 0 (A11)

∂ ρu α ∂t [START_REF] Zuckerman | Jet impingement heat transfer: Physics, correlations, and numerical modeling[END_REF] + ∂ Π

(1) αβ ∂ x β = 0 (A12) Thus, the second order moment of f

(1) i is:

Π The spatial derivative of the third order equilibrium moment with error is defined as:

Π eq * αβ γ ∂ x γ = ∂ ρu α u β u γ ∂ x γ + c 2 s ( ∂ ρu α ∂ x β + ∂ ρu β ∂ x α + ∂ ρu γ ∂ x γ δ αβ ) + ∂ Err αβ γ ∂ x γ (A15)
Inserting this term into Eq. A14, and performing some simplifications, we obtain The second moment Ψ αβ of the forcing term ψ i can thus be set to:

Π (1) αβ = -τ c 2 s ( ∂ ρu α ∂ x β + ∂ ρu β ∂ x α ) -u α ∂ ρc 2 s θ ∂ x β + u β ∂ ρc 2 s θ ∂ x α + ∂ ρc
Ψ αβ =

FIG. 1 :

 1 FIG. 1: Schematic view of the computational domain

  convective times τ conv = D U b or 26 characteristic times τ St = D StU b based on the Strouhal number of St = f D U b

FIG. 2 :FIG. 3 :

 23 FIG. 2: Zoomed X -Y cross-section of instantaneous non-dimensional vorticity magnitude, as well as nodes of Mesh 2, showing the additional zones of refinement that are not present in Mesh 1 in the boxes. The location of the refinement zones are identical for Mesh 2 and Mesh 3. The dashed lines represent the radial positions r/D = 0, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, and 2.5, which are used to validate the flow statistics in Section IV.

2 bFIG. 4 :

 24 FIG. 4: Flow statistics at the outlet of the nozzle y/D = 2, showing mean flow, velocity fluctuations, and Reynolds shear stress. ( Mesh 1), ( Mesh 2), ( Mesh 3), ( DNS of Wu, et al. 34 )

6 -u y u r /U 2 b- 2 FIG. 5 :

 6225 FIG. 5: Azimuthally averaged flow statistics along y in the jet for r/D = 0 (first row), r/D = 0.25 (second row), and r/D = 0.5 (third row). ( Mesh 1), ( Mesh 2), ( Mesh 3), (• Tummers, et al. 37 )

FIG. 7 :

 7 FIG. 7: Azimuthally averaged Nusselt number (a) and azimuthally averaged normalized wall-shear stress (b) as a function of radial position. ( Mesh 1), ( Mesh 2), (Mesh 3), (• Tummers, et al.[START_REF] Tummers | Turbulent flow in the near field of a round impinging jet[END_REF] ), ( Fénot, et al.[START_REF] Fenot | Local heat transfer due to several configurations of circular air jets impinging on a flat plate with and without semi-confinement[END_REF] ), ( Baugh and Shimizu 40 )

FIG. 8 :

 8 FIG. 8: Iso-contours of instantaneous p = p ∞ -100 Pa showing the primary structures in the free and wall jet, along with the instantaneous Nusselt number distribution. Solid black lines represent r/D = 1.0, r/D = 2.0, and r/D = 3.0

FIG. 9 :

 9 FIG. 9: Iso-contours of the PSD of the azimuthally averaged pressure signal in the (St, y/D) plane at r/D = 0.5 (top) and r/D = 1.25 (bottom)

FIG. 12 :

 12 FIG. 11: Spatiotemporal map of key azimuthally averaged surface quantities

FIG. 13 :

 13 FIG. 13: Near-wall flow topology around thermal spots

  iβ f eq i = Π eq αβ = ρu α u β + ρc 2 s θ δ αβ (A4) ∑ i c iα c iβ c iγ f eq i =Π eq * αβ γ = ρu α u β u γ + ρc 2 s (u α δ β γ + u β δ αγ + u γ δ αβ )+ Err αβ γ (A5)

∂

  (ρu α u β + ρc 2 s θ δ αβ ) ∂ x β = 0 (A9) ∂ (ρu α u β + ρc 2 s θ δ αβ ) ∂t (0) + ∂ Π eq * αβ ) + Ψ αβ (A10)

TABLE I :

 I Summary of the properties of the three meshes used 290K. This results in a bulk Courant-Friedrichs-Lewy (CFL) number of 0.77 based on U b and T ∞ .

  ρc 2 represents the term meant to eliminate the errors introduced by the D3Q19 lattice and is defined as:

		s	2 3		∂ u γ ∂ x γ	δ αβ + c 2 s (u α	∂ ρ(θ -1) ∂ x β	+ u β	∂ ρ(θ -1) ∂ x α	)+
	∂ ρc 2 s (θ -1) ∂t (0)	δ αβ +	∂ Err αβ γ ∂ x γ	,
										(A18)
	where	∂ Err αβ γ ∂ x γ
								∂ ρu x u x u x	∂ ρu x u y u z	∂ ρu x u y u z	
	∂ Err αβ γ ∂ x γ	= -	  	∂ x ∂ ρu x u y u z ∂ z ∂ ρu x u y u z	∂ z ∂ ρu y u y u y ∂ y ∂ ρu x u y u z	∂ y ∂ ρu x u y u z ∂ x ∂ ρu z u z u z	  	(A19)
									∂ y	∂ z	∂ z
	Thus, yielding the traceless viscous stress tensor found in the
	Navier-Stokes equations:
		Π (1) αβ = -τρc 2 s	∂ u α ∂ x β	+	∂ u β ∂ x α	-	2 3	∂ u γ ∂ x γ	δ αβ	(A20)

ACKNOWLEDGEMENTS

The authors would like to thank Florian Renard and Thomas Astoul for their insightful advice. Acknowledgments are also expressed to Safran Aircraft Engines and ANRT/CIFRE convention N°2019/1220 for their financial support, as well as for access to high performance computing resources. This work has been carried out using ProLB, a lattice-Boltzmann solver currently being developed by a consortium of French universities, research institutions, and industrial partners.

The authors have no conflicts of interest to disclose.