S Wilhelm

) J Jacob

P Sagaut

I. INTRODUCTION

Cartesian grids are attractive compared to body-fitted grids for the flow prediction around complex geometries. Automatic grid generation is indeed easier with Cartesian grids, and higher grid qualities are obtained around sharp edges and complex shapes. Moreover, numerical methods are more accurate and simpler to implement on Cartesian grids. [START_REF] Berger | Progress towards a Cartesian cut-cell method for viscous compressible flow[END_REF][START_REF] Tamaki | Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method[END_REF][START_REF] Hu | Towards efficient viscous modeling based on Cartesian methods for automated flow simulation[END_REF] The Lattice Boltzmann Method (LBM) is an alternative to the classical approach involving the resolution of the Navier-Stokes equations. It is typically used on Cartesian grids with cubic cells which contribute to its efficiency.

Accounting for the shape of an arbitrary solid body using Cartesian grids is mainly to be done, thanks to two methods: the cut-cell method and the Immersed Boundary Method (IBM). In the cut-cell grid, no cells are defined inside the solid body, and cells are cut at the intersection between the volumetric grid and the solid. [START_REF] Berger | Progress towards a Cartesian cut-cell method for viscous compressible flow[END_REF] In the immersed boundary method, the volumetric grid is defined in the whole computational domain, including inside the solid body. A force-field is then added to the governing flow equations to simulate the presence of the solid. [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Pinelli | Immersed-boundary methods for general finite-difference and finite-volume Navier-Stokes solvers[END_REF][START_REF] Favier | A lattice Boltzmann-immersed boundary method to simulate the fluid interaction with moving and slender flexible objects[END_REF] In both methods, the boundary nodes do not necessarily lie on the solid boundary. A reconstruction method for the velocity at the boundary nodes based on an interpolation step is then necessary. On cut-cell grids, extrapolation methods have been proposed to estimate the velocity at the boundary nodes from the velocity predicted at neighbouring fluid nodes. [START_REF] Berger | Progress towards a Cartesian cut-cell method for viscous compressible flow[END_REF][START_REF] Hu | Towards efficient viscous modeling based on Cartesian methods for automated flow simulation[END_REF] Similarly, interpolation methods are also used in IBM in which the velocity is reconstructed at solid nodes to recover the correct boundary condition at the wall. Examples can be found in the literature for Navier-Stokes computations [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF][START_REF] Jafari | Terrain effects on wind flow: Simulations with an immersed boundary method[END_REF] and for LBM. [START_REF] Pellerin | An implementation of the Spalart-Allmaras turbulence model in a multi-domain lattice Boltzmann method for solving turbulent airfoil flows[END_REF] For the prediction of turbulent flows at high Reynolds numbers relevant for real engineering applications, the a) Electronic mail: sylvia.wilhelm@univ-amu.fr Reynolds-Averaged Navier-Stokes (RANS) approach can be used in conjunction with Cartesian grids. Iaccarino and Verzicco [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF] studied the use of the Spalart-Allmaras (SA) turbulence model coupled with the IBM. A RANS-IBM methodology was also proposed by Capizzano [START_REF] Capizzano | Turbulent wall model for immersed boundary methods[END_REF] using the kg and the k -ω TNT turbulence models. Kalitzin and Iaccarino [START_REF] Kalitzin | Turbulence modeling in an immersedboundary RANS method[END_REF] used the k -ω, kg, and Spalart-Allmaras RANS turbulence models with the IBM approach for the prediction of the flow inside a turbine blade passage. The flow prediction around airfoils has been performed with the k -ω Shear Stress Transport (SST) turbulence model using the Navier-Stokes IBM approach [START_REF] Zhou | RANS simulation of high-Re turbulent flows using an immersed boundary method in conjunction with wall modeling[END_REF] and with the Spalart-Allmaras turbulence model using the cut-cell approach [START_REF] Berger | Progress towards a Cartesian cut-cell method for viscous compressible flow[END_REF] or using LBM with an immersed surface. [START_REF] Pellerin | An implementation of the Spalart-Allmaras turbulence model in a multi-domain lattice Boltzmann method for solving turbulent airfoil flows[END_REF][START_REF] Li | Non-body-fitted Cartesian-mesh simulation of highly turbulent flows using multi-relaxation-time lattice Boltzmann method[END_REF] When a RANS turbulence model is used on a Cartesian grid, the turbulent viscosity has to be determined at the boundary node. [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF] The turbulence quantities such as k and ω (or ) in two-equation models exhibit steep gradients and local extrema near the wall within the buffer layer. High accuracy of the numerical scheme and sufficient grid resolution are therefore required near the wall. According to Wilcox, [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF] the scaled distance to the wall y + should be lower than 2.5 for the integration of the k -ω model through the viscous sublayer and lower than 1 for the kmodel. As explained by Iaccarino and Verzicco, [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF] this is an issue for RANS simulations with the immersed solid boundary since the explicit resolution of the boundary layer with a Cartesian grid is too expensive at a high Reynolds number. Grid refinement in the wall normal direction only is not enabled using a uniform Cartesian grid, and isotropic grid refinement in the boundary layer may lead to an important number of cells, much higher than classical grids used for direct numerical simulation. Therefore, a wall model must be used inside the first off-wall node for both the velocity profile and the quantities related to the turbulence model.

Two methods are mainly used for wall modeling when addressing the issue of predicting the mean velocity profile Aix Marseille University, CNRS, Centrale Marseille, M2P2, Marseille, France

An explicit power-law-based wall model for lattice Boltzmann method-Reynolds-averaged numerical simulations of the flow around airfoils

In this paper, an explicit wall model based on a power-law velocity profile is proposed for the simulation of the incompressible flow around airfoils at high Reynolds numbers. This wall model is particularly suited for the wall treatment involved in Cartesian grids. Moreover, it does not require an iterative procedure for the friction velocity determination. The validation of this power-law wall model is assessed for Reynolds-averaged Navier-Stokes simulations of the flow around a two-dimensional airfoil using the lattice Boltzmann approach along with the Spalart-Allmaras turbulence model. Good results are obtained for the prediction of the aerodynamic coefficients and the pressure profiles at two Reynolds numbers and several angles of attack. The explicit power-law is thus well suited for a simplified near-wall treatment at high Reynolds numbers using Cartesian grids.

within the first cell. [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF] First, a simplified equation for the turbulent boundary layer (TBLE) equipped with an algebraic mixing length eddy viscosity model can be solved on an embedded grid between the wall and the first off-wall node; see, e.g., Ref. 10. Such an approach may be interpreted as a dramatically simplified version of the coupling between a boundary layer code (within the first cell) and a classical Navier-Stokes RANS solver. Second, the problem can be further simplified assuming that, when properly scaled, the tangential velocity profile follows a universal analytical algebraic law. Most existing wall models rely on a logarithmic-type profile. The tangential velocity is then given by a relation between the tangential velocity, the distance to the wall, and the friction velocity. Algebraic wall models are more often used than TBLE for RANS simulations because of their simplicity. For body-fitted grids, the no-slip boundary condition is thus replaced by the prescription of the wall shear stress. The latter is computed assuming that the first computed point lies within the inner layer and inverting the non-linear analytical velocity profile, thanks to an iterative procedure. It has been shown that the wall model should be consistent with the turbulence model used outside the boundary layer. Such wall models have been developed for body-fitted grids in RANS [START_REF] Kalitzin | Near-wall behavior of RANS turbulence models and implications for wall functions[END_REF][START_REF] Knopp | A grid and flow adaptive wallfunction method for RANS turbulence modelling[END_REF] and Large Eddy Simulation (LES). [START_REF] Kawai | Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy[END_REF][START_REF] Kawai | Wall-modeled large-eddy simulation of high Reynolds number flow around an airfoil near stall condition[END_REF] Along with the wall model for the velocity profile, wall models for the quantities of RANS turbulence models have been proposed. Kalitzin et al. [START_REF] Kalitzin | Near-wall behavior of RANS turbulence models and implications for wall functions[END_REF] derived the analytical solutions of turbulence equations of the Spalart-Allmaras, k -ω, kg, 2f RANS turbulence models in the viscous sublayer and logarithmic layer while tabulated functions are defined in the buffer layer. Knopp et al. [START_REF] Knopp | A grid and flow adaptive wallfunction method for RANS turbulence modelling[END_REF] proposed a near-wall solution for ω different from the one proposed in the work of Kalitzin et al. [START_REF] Kalitzin | Near-wall behavior of RANS turbulence models and implications for wall functions[END_REF] This solution for ω is used by Zhou, [START_REF] Zhou | RANS simulation of high-Re turbulent flows using an immersed boundary method in conjunction with wall modeling[END_REF] while k is determined to ensure consistency between the velocity profile and the turbulence model in the boundary layer. Similarly, Capizzano 10 calculated the boundary value for k from the mixing length eddy viscosity model used in the TBLE and the analytical expression of the second variable of the turbulence model. In the studies of Iaccarino and Verzicco 7 and Kalitzin and Iaccarino, [START_REF] Kalitzin | Turbulence modeling in an immersedboundary RANS method[END_REF] the authors took advantage of the linear behavior of the turbulent quantity in the Spalart-Allmaras model.

The use of Cartesian grids further complexifies the problem since the velocity is not predicted at the boundary node and therefore has to be reconstructed. The wall model is thus used for high Reynolds number cases instead of the interpolation methods mentioned above for laminar flow simulations. [START_REF] Capizzano | Turbulent wall model for immersed boundary methods[END_REF] Moreover, the distance to the wall of the first off-wall node is highly variable on a Cartesian grid. This may induce numerical instabilities when using a wall model based on this distance. Many studies have been performed in the field of IBM to develop methodologies of application of a wall model on a Cartesian grid. They mainly consist of first applying the wall model to a point where the velocity is predicted in order to calculate the friction velocity from the wall model. The velocity at the boundary node is thus determined using the wall model with this calculated friction velocity. Research studies can be distinguished by the choice of the point to first apply the wall model which can be defined at a constant [START_REF] Capizzano | Turbulent wall model for immersed boundary methods[END_REF][START_REF] Kalitzin | Turbulence modeling in an immersedboundary RANS method[END_REF][START_REF] Zhou | RANS simulation of high-Re turbulent flows using an immersed boundary method in conjunction with wall modeling[END_REF] or variable [START_REF] Tessicini | Wall modeling for large-eddy simulation using an immersed boundary method[END_REF][START_REF] Roman | A simple wall-layer model for large eddy simulation with immersed boundary method[END_REF][START_REF] Kang | An improved near-wall modeling for large-eddy simulation using immersed boundary methods[END_REF] distance to the wall. Interpolation methods are again necessary to determine the velocity at this point.

As mentioned above, an iterative procedure is generally required to compute the friction velocity by inverting the algebraic wall function. Newton's method is usually used with several iterations for convergence. [START_REF] Tamaki | Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method[END_REF][START_REF] Zhou | RANS simulation of high-Re turbulent flows using an immersed boundary method in conjunction with wall modeling[END_REF] Alternatively, tabulated functions can be found [START_REF] Kalitzin | Near-wall behavior of RANS turbulence models and implications for wall functions[END_REF] by performing the inversion of the wall function once on a well-resolved numerical simulation and storing the inverse function in a look-up table. This iterative method is time-consuming and may introduce numerical errors in the determination of flow quantities at the boundary nodes.

In this work, an explicit wall model is proposed that do not rely on an iterative procedure for the determination of the friction velocity. Moreover, its use on Cartesian grids is very simple and enables the evaluation of the velocity at boundary nodes without the previous determination of the friction velocity. The purpose of this work is to assess the validity and the robustness of this wall model on the prediction of the aerodynamic flow around an airfoil at a high Reynolds number. The Lattice Boltzmann Method (LBM) is used for the prediction of the flow. The LBM is an appealing approach for Computational Fluid Dynamics (CFD) compared to the classical Navier-Stokes-based methods. The LBM is a local method particularly well suited for massively parallel simulations. [START_REF] Krüger | The Lattice Boltzmann Method[END_REF][START_REF] Pellerin | Equilibrium distributions for straight, curved, and immersed boundary conditions in the lattice Boltzmann method[END_REF] Contrary to the Navier-Stokes equations, it does not involve the computation of a non-local non-linear term. Moreover, the pressure is obtained by means of an equation of state easier to solve than the Poisson equation involved in classical methods. [START_REF] Chen | Lattice Boltzmann method for fluid flows[END_REF] LBM has become very popular for the prediction of subsonic aerodynamic flows. Low angles of attack on the airfoil are considered so that the flow remains attached and is resolved using a RANS approach based on the Spalart-Allmaras turbulence model.

The paper is organized as follows. The governing equations considered in this work are presented in Sec. II. The numerical method is then discussed in Sec. III. The proposed wall model and its implementation are described in Sec. IV. In Sec. V, the proposed model and its practical implementation are assessed considering the flow prediction around a two-dimensional (2D) clean airfoil at several angles of attack and two Reynolds numbers. Finally, conclusions of this study are given in Sec. VI.

II. MACROSCOPIC GOVERNING EQUATIONS

A. Navier-Stokes equations

We consider the three-dimensional Navier-Stokes equations for an incompressible flow,

∂u i ∂x i = 0, (1)
∂u i ∂t + ∂u i u j ∂x j = - 1 ρ ∂p ∂x i + 1 ρ ∂τ ij ∂x j , (2)
where u i is the velocity component in the x i direction, ρ is the density, p is the static pressure, and τ ij is the total stress tensor.

In this work, the RANS approach is used so that the macroscopic quantities in Eqs. ( 1) and ( 2) are mean quantities in the sense of Reynolds averaging. Note that the mean operator . is omitted for the sake of simplicity. Therefore, τ ij is the sum of the viscous and Reynolds stress tensors,

τ ij = 2(µ + µ t )S ij , ( 3
)
where µ is the dynamic viscosity of the fluid, µ t is the turbulent viscosity provided by the RANS turbulence model, and S ij is the strain rate tensor,

S ij = 1 2 ∂u i ∂x j + ∂u j ∂x i . (4)

B. Spalart-Allmaras turbulence model

The Spalart-Allmaras (SA) turbulence model is used for the RANS calculations in this work. This one-equation model is more versatile and powerful than algebraic models and easier to implement than two-equation models. Moreover, it is known for its robustness. [START_REF] Pope | Turbulent Flows[END_REF] This model was initially developed for the prediction of aerodynamic flows using dimensional analysis, Galilean invariance, and selected dependence on the molecular viscosity. [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] The SA model consists of a transport equation for the working viscosity ν defined from the eddy viscosity ν t . According to the Prandtl assumption, the eddy viscosity has a linear behavior in the logarithmic zone of the boundary layer: ν t = κu τ y, where y is the distance to the wall and u τ is the friction velocity. [START_REF] Kalitzin | Near-wall behavior of RANS turbulence models and implications for wall functions[END_REF] The working viscosity ν is defined in order to preserve this linear behavior ν = κu τ y in the logarithmic layer but also in the buffer layer and in the viscous sublayer. The eddy viscosity ν t is calculated from the working viscosity ν as follows:

ν t = νf v1 , (5)
with the damping function f 1 defined as

f v1 = χ 3 χ 3 + c 3 v1 with χ = ν ν , ( 6
)
where ν is the kinematic viscosity of the fluid and c 1 is a constant of the model. The Spalart-Allmaras model has been originally developed in 1992, [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] and several variants of the model have since been proposed. The baseline and negative SA model used in this work is described in this section. Two other variants of the original model have been tested which are described in the Appendix along with the results obtained with these variants.

The original version of the SA model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF][START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] without the laminar suppression term and the trip term is first considered. Indeed, as explained in the work of Allmaras et al., [START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[END_REF] the laminar term is negligible in fully turbulent flows with large free-stream values for ν. In Ref. 29, a version of the SA model is also proposed for the case χ < 0. On coarse grids and during unphysical transient states, the discrete solution of the SA model can lead to negative values of ν. [START_REF] Oliver | A high-order, adaptive, discontinuous Galerkin finite element method for the Reynolds-averaged Navier-Stokes equations[END_REF] A modification of the production, destruction, and diffusion terms of the original SA model is proposed to avoid negative values of the working viscosity. The model is thus as follows:

• For χ > 0, the equation for the working viscosity is

∂ ν ∂t + ∂ νu i ∂x i = c b1 S ν Production + 1 σ       ∂ ∂x i (ν + ν) ∂ ν ∂x i + c b2 ∂ ν ∂x i 2      Diffusion -c w1 f w ν d w 2 Destruction , (7)
where d w is the distance from the nearest wall. Similar to the eddy viscosity, the vorticity magnitude Ω is modified into S through a damping function f v 2 in order to preserve its log-layer behavior S = u τ /(κy) in the whole boundary layer. However, similarly for ν, the discrete solution can lead to negative S values. A correction is proposed in the work of Allmaras et al. [START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[END_REF]

such that S =        Ω + S if S ≥ -c v2 Ω Ω + Ω(c 2 v2 Ω+c v3 S) (c v3 -2c v2 )Ω- S if S < -c v2 Ω , ( 8
) with S = ν κ 2 d 2 w f v2 (9)
and

f v2 = 1 - χ 1 + χf v1 . ( 10
)
The function f w is defined to accelerate the decaying behavior of the destruction term in the outer region of the boundary layer,

f w = g 1 + c 6 w3 g 6 + c 6 w3 1 6
,

g = r + c w2 (r 6 -r), r = min ν Sκ 2 d 2 w , 10 . (11)
• For χ < 0, the following equation is solved instead:

∂ ν ∂t + ∂ νu i ∂x i = c b1 (1 -c t3 )Ω ν Production + 1 σ       ∂ ∂x i (ν + νf n ) ∂ ν ∂x i + c b2 ∂ ν ∂x i 2      Diffusion + c w1 ν d w 2 Destruction , (12)
with

f n = c n1 + χ 3 c n1 -χ 3 . ( 13
)
Moreover, the turbulent viscosity ν t is set to 0 if χ is negative.

The model constants are σ = 2/3, κ = 0.41, c b1 = 0.1355, c b2 = 0.622,

c w1 = c b1 κ 2 + 1 + c b2 σ , c w2 = 0.3, c w3 = 2, c v1 = 7.1, c v2 = 0.7, c v3 = 0.9, c n1 = 16. (14)

III. NUMERICAL METHOD

A. The lattice Boltzmann method

In this work, the Lattice Boltzmann Method (LBM) [START_REF] Krüger | The Lattice Boltzmann Method[END_REF][START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF][START_REF] Guo | The Lattice Boltzmann Method and its Applications in Engineering[END_REF] is used for the numerical resolution of the macroscopic governing equations ( 1) and ( 2). This approach describes the dynamics of a group of particles that collide and propagate over a discrete lattice based on the Boltzmann equation at the mesoscopic level. This equation describes the evolution of a particle distribution function f = f ( x, c, t) that represents the probability density of particles with velocity c at time t and position x. In LBM, the Boltzmann equation is discretized in the velocity space on a DdQq (d dimensions and q velocities) lattice.

In this work, the D3Q19 lattice given by

c α =          (0, 0, 0), α = 0 (±1, 0, 0), (0, ±1, 0), (0, 0, ±1), α = 1 -6 (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1), α = 7 -18 (15)
is used. The associated lattice Boltzmann equation with the time step ∆t and space step -→ ∆x = c α ∆t is given by 33

f α ( x + c α ∆t, t + ∆t) -f α ( x, t) = ∆t Ω α ( x, t), ( 16
)
where Ω α is the collision operator.

The macroscopic quantities (density ρ, momentum ρu i , and momentum flux tensor Π ij ) are then computed from the velocity moments of the distribution functions as follows:

ρ = q-1 α=0 f α , ( 17a
)
ρu i = q-1 α=0 f α c αi , (17b)
Π ij = q-1 α=0 f α c αi c αj . ( 17c
)
The collision operator Ω α is usually interpreted as a relaxation toward an equilibrium state. Using the Bhatnagar-Gross-Krook (BGK) model, [START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one component systems[END_REF] one has

Ω α = - 1 τ (f α -f eq α ), (18)
with τ being the relaxation time and f eq α being the local equilibrium function given by

f eq α = ρω α 1 + c αi u i c 2 s + 1 2c 4 s Q αij u i u j , (19)
where c s is the speed of sound, ω α are the weighting coefficients of the D3Q19 scheme, and

Q αij = c αi c αj -c 2 s δ ij .
In order to ensure good numerical stability and convergence properties in high Reynolds number flows, the hybrid recursive regularized BGK collision operator proposed in the work of Jacob et al. [START_REF] Jacob | A new hybrid recursive regularized Bhatnagar-Gross-Krook collision model for lattice-Boltzmann-method based large-eddy simulation[END_REF] is used. The distribution function is split into the sum of the equilibrium f eq α and the nonequilibrium f neq α distribution functions, which leads to the following regularized collision operator:

Ω α = - 1 τ f neq α , (20)
where

f neq α = σ HRR f neq,LBM α + (1 -σ HRR ) f neq,FD α , (21)
with

f neq,FD α = - τ ρω α 2c 2 s Q αij ∂u j ∂x i + ∂u i ∂x j , f neq,LBM α = Q αij 2c 4 s α c αi c αj (f α -f eq α ), (22)
with σ HRR = 0.98 in the present study.

B. Discretization of the SA model

The equation of the Spalart-Allmaras model is solved for χ = ν/ν. A second-order accurate finite difference scheme is used for the resolution of the equations of the Spalart-Allmaras model. The temporal derivative is discretized using the first-order Euler explicit scheme. A hybrid upwindcentered scheme is used for the advection term. Gradients and Laplacian operators are calculated using a centered scheme.

The eddy viscosity ν t calculated by the SA model is introduced in the relaxation time of the regularized BGK model

such as τ = 1 2 + ν + ν t c 2 s
in order to take into account turbulence effects in the collision operator.

C. Near-wall treatment

In this work, a Cartesian cut-cell grid is used to solve the LBM scheme in the fluid domain. The lattice Boltzmann method needs particular treatment at boundary nodes which are the nodes that do not have all their neighboring nodes in the fluid domain. The distribution functions in directions of these missing neighboring nodes are unknown. The distribution functions at the boundary nodes need therefore to be reconstructed to complete the LBM scheme. There are two main approaches for this reconstruction: the bounce-back approach or the wet boundary conditions. [START_REF] Krüger | The Lattice Boltzmann Method[END_REF] In this paper, an approach belonging to the second category is used in which the distribution functions at boundary nodes are reconstructed based on macroscopic variables and their derivatives. This approach is similar to the one proposed by Verschaeve and Müller. [START_REF] Verschaeve | A curved no-slip boundary condition for the lattice Boltzmann method[END_REF] The distribution function at a boundary node is calculated as the sum of the equilibrium f eq α and non-equilibrium f neq α distribution functions given by Eqs. (19) and (22). The non-equilibrium distribution function is calculated according to f neq,FD α . The missing distribution function can then be determined from the macroscopic variables u and ρ at the boundary node. The reconstruction method of the macroscopic values at a boundary node is illustrated in Fig. 1, where is the boundary node of interest and and nodes are the neighboring fluid nodes of at which all required information is available. We address here high Reynolds number flows. As explained in the Introduction, it would be too expensive to explicitly resolve the boundary layer using a Cartesian grid. A wall model is thus used to evaluate the tangential velocity at the boundary node. The wall model is a function g defining the universal rescaled tangential velocity profile u + in the turbulent boundary layer as a function of the normalized distance to the wall y + , u + = g(y + ),

with

u + = u u τ and y + = yu τ ν , ( 24
)
where u is the tangential velocity at the distance y to the wall, ν is the kinematic viscosity of the fluid, and u τ is the friction velocity related to the wall shear stress such as τ w = ρu 2 τ . Malaspinas and Sagaut [START_REF] Malaspinas | Wall model for large-eddy simulation based on the lattice Boltzmann method[END_REF] proposed a method to use a wall model with the LBM approach for a body-fitted grid. As explained in the Introduction, methods of application of a wall model with a Cartesian grid have been developed principally for the IBM. Methods for the velocity reconstruction at the boundary node that use only fluid points can thus be adapted in this work. A common practice consists of defining a fictitious node on the line normal to the wall passing through the boundary node and located at an arbitrarily fixed distance to the wall, thereby imitating a body-fitted grid. [START_REF] Capizzano | Turbulent wall model for immersed boundary methods[END_REF][START_REF] Kalitzin | Turbulence modeling in an immersedboundary RANS method[END_REF][START_REF] Bernardini | On the suitability of the immersed boundary method for the simulation of high-Reynolds-number separated turbulent flows[END_REF] This fictitious node, called Ref point in this work, is represented with a • in Fig. 1. This point is located at a distance of 2.5dx from the wall in the present work, where dx is the local grid size. The tangential velocity at the Ref point is calculated using the modified Shepard's Inverse Distance Weighting (IDW) interpolation of tangential velocities resolved at the and nodes. It is reminded that IDW originates in the minimization of an error cost function related to a measure of deviations between tuples of interpolating points and tuples of interpolated points, leading to an optimized interpolation procedure that can be interpreted as a generalization of Lagrange interpolation.

The wall model ( 23) is then applied at the Ref point in order to determine the friction velocity u τ . An iterative procedure is often needed to determine u τ from Eq. ( 23). Then, using this obtained friction velocity u τ , the wall model ( 23) is applied at the boundary node to determine the tangential velocity u.

The density at the boundary node is calculated by using IDW interpolation of densities at the and nodes.

Finally, the normal velocity at the boundary node is equated to the velocity of the solid boundary in the normal direction, which is equal to zero for a stationary wall.

Macroscopic values are thus known at the boundary node. The velocity gradients at the boundary node are calculated using one-sided finite difference discretization in the fluid. The distribution function at the boundary node can then be calculated using Eqs. ( 19) and (22).

IV. A NEW EXPLICIT POWER-LAW-BASED WALL MODEL

A. Velocity profile in a turbulent boundary layer

As explained in Sec. III C, a wall model describes the profile of the scaled tangential velocity in the turbulent boundary layer. In the viscous sublayer, the velocity profile is linear. The logarithmic law of the wall is classically used to describe the velocity profile in the boundary layer outside the viscous sublayer,

u + = 1 κ ln(y + ) + C, ( 25
)
where κ ∼ 0.41 is the von Karman constant and C is a constant. This is the most popular model for engineering applications. An iterative procedure is necessary to determine the friction velocity from the non-linear relation (25) with known tangential velocity. However, several authors have questioned the validity of the logarithmic law of the wall and in particular the universal character of its parameters κ and C, which are assumed to be independent of the Reynolds number. [START_REF] Buschmann | Recent developments in scaling of wall-bounded flows[END_REF][START_REF] George | Is there a universal log law for turbulent wall-bounded flows?[END_REF] Marusic et al. [START_REF] Marusic | On the logarithmic region in wall turbulence[END_REF] have shown, based on experiments in the laboratory boundary layer, pipe flow, and atmospheric surface layer, that the mean velocity and streamwise turbulent intensities have a universal logarithmic behavior at Reynolds numbers ranging from 2 × 10 4 to 6 × 10 5 . The associated coefficients in Eq. ( 25) are κ = 0.39 and C = 4.3.

A review of the existing laws of the wall and scaling proposed by different authors for a turbulent boundary layer is done in the work of Buschmann and Gad-el Hak. [START_REF] Buschmann | Recent developments in scaling of wall-bounded flows[END_REF] In particular, it is shown that the boundary layer can also be described by a power-law model. Prandtl proposed a power law for the whole boundary layer at a low Reynolds number,

u + = A(y + ) 1 7 , ( 26
)
where A is a coefficient which may depend on the Reynolds number. A velocity profile following a 1/n-th power law was also found in Nikuradse's experiments for a pipe flow. The exponent 1/n of the law depends on the Reynolds number and decreases with increasing Reynolds number. [START_REF] Pope | Turbulent Flows[END_REF] According to Prandtl, the logarithmic law of the wall is the limit of the power law at high Reynolds numbers. [START_REF] Schlichting | Boundary-Layer Theory[END_REF] The powerlaw model with n = 7 and constant A = 8.3 was used, for example, by Werner and Wengle [START_REF] Werner | Large-eddy simulation of turbulent flow over and around a cube in a plate channel[END_REF] in LES. The power-law model of Werner and Wengle [START_REF] Werner | Large-eddy simulation of turbulent flow over and around a cube in a plate channel[END_REF] was used in the work of Chang et al. [START_REF] Chang | Simulations of laminar and turbulent flows over periodic hills with immersed boundary method[END_REF] for the prediction of turbulent flows over periodic hills in LES using an IBM on a Cartesian grid. Results similar to those obtained with a TBLE were observed. Murakami et al. [START_REF] Murakami | Three-dimensional numerical simulation of air flow around a cubic model by means of large eddy simulation[END_REF] also used a power-law model with n = 4 on a flat plate and n = 2 on a cube. More recently, Afzal [START_REF]Power law and log law velocity profiles in fully developed turbulent pipe flow: Equivalent relations at large Reynolds numbers[END_REF] showed that both the power-law and the logarithmic law velocity profiles are solutions of the matching between the asymptotic forms of the law of the wall and the defect law in the overlap region of the boundary layer. On the other hand, some authors of experimental studies found that both the power law and logarithmic law exist and describe different parts of the boundary layer. [START_REF] Zagarola | Log laws or power laws: The scaling in the overlap region[END_REF] In this work, the Spalart-Allmaras turbulence model is used for the RANS simulations. The following consistent wall model for velocity was derived in the study of Allmaras et al.:

29 u + (y + ) = B + c 1 ln((y + + a 1 ) 2 + b 2 1 ) -c 2 ln((y + + a 2 ) 2 + b 2 2 ) -c 3 Arctan(y + + a 1 , b 1 ) -c 4 Arctan(y + + a 2 , b 2 ), (27)
where the constants are defined in Ref. 29. In particular, the SA model asymptotically predicts the log-law at large Reynolds numbers. However, this model is complex and requires the use of an iterative procedure to compute the friction velocity. Figure 2 displays a comparison of the power-law model of Werner and Wengle [START_REF] Werner | Large-eddy simulation of turbulent flow over and around a cube in a plate channel[END_REF] with the wall model of the SA model and the logarithmic law of Marusic et al. [START_REF] Marusic | On the logarithmic region in wall turbulence[END_REF] The power law is in agreement with the logarithmic laws for 50 < y + < 1000 with less than 2% of difference in the u + velocity predicted by the power law and the logarithmic law. FIG. 2. Comparison of the logarithmic and power law models.

On the other hand, Abe and Antonia [START_REF] Abe | Relationship between the energy dissipation function and the skin friction law in a turbulent channel flow[END_REF] have recently obtained the relation between the skin friction C f and the Reynolds number Re b in a turbulent channel flow with the bulk velocity U b . This relation, reported in Eq. ( 28), is obtained from a logarithmic model for the velocity in the channel flow,

1 C f = 1.80 ln(Re b C f ) -0.163. ( 28
)
Using the same method as in the study of Abe and Antonia, [START_REF] Abe | Relationship between the energy dissipation function and the skin friction law in a turbulent channel flow[END_REF] the corresponding relation is obtained for the power law, u

+ = A(y + ) B , C f = 2 * 8 B * (B + 1) 2 A 2 1 B+1 Re -2B B+1 b . ( 29
)
An iterative procedure is required in order to determine C f from Eq. ( 28), whereas ( 29) is explicit. Figure 3 presents the dependence of the skin friction C f to the Reynolds number Re b in a turbulent channel flow using a logarithmic law or a power law for the velocity in the channel flow. Both models have a similar behavior for large Reynolds numbers Re b > 10 3 .

The logarithmic and power law wall models are in agreement with the range of Reynolds numbers and y + values considered in this work.

Several formulas for the parameter A and the exponent B = 1/n have been proposed, based on experimental data [START_REF] Zagarola | Log laws or power laws: The scaling in the overlap region[END_REF] or by identification of the power and logarithmic laws. [START_REF]Power law and log law velocity profiles in fully developed turbulent pipe flow: Equivalent relations at large Reynolds numbers[END_REF][START_REF] Castro-Orgaz | Power-law velocity profile in turbulent boundary layers: An integral Reynolds-number dependent solution[END_REF] Castro-Orgaz and Dey 49 used a power law to model the entire boundary layer with parameters calculated by analogy with a logarithmic law and a law of the wake. These parameters are functions of the Reynolds number. Moreover, formulas have been proposed to take into account the pressure gradient and roughness effects in the power-law model. [START_REF] Castro-Orgaz | Hydraulics of developing chute flow[END_REF][START_REF] Zoumakis | The dependence of the power-law exponent on surface roughness and stability in a neutrally and stably stratified surface boundary layer[END_REF] Breuer et al. [START_REF] Breuer | Development of wall models for les of separated flows using statistical evaluations[END_REF] proposed an extended version of the wall model of Werner and Wengle, [START_REF] Werner | Large-eddy simulation of turbulent flow over and around a cube in a plate channel[END_REF] taking into account the effect of the streamwise pressure gradient for LES calculations. Additional improvements were also proposed using the artificial eddy viscosity concept and a new definition of the relative thickness of the FIG. 3. Comparison of the skin friction for the logarithmic and power law models in a turbulent channel flow.

viscous sublayer. This extended wall model was applied to a periodic arrangement of hills involving separated flows. The prediction of the separated zone was slightly improved compared to the original model, but the improvement linked to the inclusion of the pressure gradient effect was small.

B. Power-law model for the turbulent boundary layer

In this work, a power-law model is used to describe the velocity profile in the turbulent boundary layer because of its simplicity of implementation and its ease of use. The 1/n power-law model with constant parameters (independent of the Reynolds number) is evaluated in this work,

u + =      y + if y + ≤ y + c A(y + ) B if y + ≥ y + c , (30)
with B = 1/7. y c is the height of the viscous sublayer and y + c = y c u τ /ν. According to Werner and Wengle, [START_REF] Werner | Large-eddy simulation of turbulent flow over and around a cube in a plate channel[END_REF] y + c = 11.81. The parameter A can be calculated by continuity of the velocity profile at y + c ,

y + c = A(y + c ) B ⇐⇒ A = (y + c ) 1-B . ( 31
)
For y + c = 11.81 and B = 1/7, one has A ≈ 8.3 which is in agreement with the value used by Werner and Wengle. [START_REF] Werner | Large-eddy simulation of turbulent flow over and around a cube in a plate channel[END_REF] The expression for the friction velocity u τ is explicitly obtained from (30), written in non-scaled units, at a distance y > y c to the wall, u(y)

u τ = A yu τ ν B , (32)
leading to u τ = u(y)

1 1+B A -1 1+B y -B 1+B ν B 1+B , ( 33
) which is an explicit relation that do not require a non-linear iterative method to compute the skin friction.

It is also convenient to check whether the point lies within the viscous sublayer or not in a fast and non-expensive way. To this end, one introduces the local Reynolds number Re l (y) = yu(y)/ν = y + u + (y + ). The critical Reynolds number is defined at the top of the viscous sublayer by Re c = Re l (y c ) = (y + c ) [START_REF] Tamaki | Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method[END_REF] and is known since y + c is fixed. If the local Reynolds number at the first off-wall point is lower than Re c , then the linear law is used. Otherwise, the exponential law is selected.

Observing that for two points in the fluid such as y 1 ≥ y c and y 2 ≥ y c ,

Re l (y 1 ) Re l (y 2 ) = y 1 u(y 1 ) y 2 u(y 2 ) = y 1 y 2 1+B , ( 34
)
the height of the viscous sublayer y c can be estimated as follows:

y c = y 2 Re c Re l (y 2 ) 1 1+B . ( 35
)
C. An explicit near-wall treatment for LBM

For the application of the power-law in LBM simulations, the reference point Ref defined in Sec. III C is used. The power law can be applied at the reference point,

u + Ref =      y + Ref if y + Ref ≤ y + c A(y + Ref ) B if y + Ref ≥ y + c . ( 36
)
If y + Ref ≤ y + c , then both the reference point and the boundary node lie within the viscous sublayer and

u(y) u Ref = y y Ref ⇐⇒ u(y) = u Ref y y Ref . ( 37
)
If y + Ref ≥ y + c , we suppose that the boundary node is also in the power-law region (i.e., y + ≥ y + c ). By dividing the expression for u + (y + ) at the boundary node [Eq. (30)] by the one for u + Ref [Eq. ( 36)], the expression for the non-scaled tangential velocity at the boundary node is directly obtained,

u(y) u Ref = y y Ref B ⇐⇒ u(y) = u Ref y y Ref B . ( 38
)
Contrary to the logarithmic law of the wall (25), the friction velocity u τ is not needed for the calculation of the velocity at the boundary node. Note that the knowledge of the coefficient A is also not necessary. The use of a wall model on a Cartesian grid is thus simplified by using a power-law model.

D. Boundary condition for the Spalart-Allmaras turbulence model

The boundary conditions for the turbulent quantities should be consistent with the RANS turbulence model used to resolve the flow in order to avoid grid dependency of the solution. [START_REF] Kalitzin | Near-wall behavior of RANS turbulence models and implications for wall functions[END_REF][START_REF] Knopp | A grid and flow adaptive wallfunction method for RANS turbulence modelling[END_REF] Moreover, according to Iaccarino and Verzicco, [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF] the resolution of the turbulent quantities of a RANS model is an important issue when using an immersed solid boundary method. The accuracy of the flow prediction is influenced by the behavior of these quantities in the vicinity of the wall. The linear behavior of ν in the whole boundary layer (see Sec. II B) is thus very attractive from a numerical point of view when used in conjunction with a Cartesian grid. The near-wall resolution of ν is thus less sensitive to the near-wall the angle of attack grid resolution than for other turbulence models and is consistent with the interpolations used in the wall treatment of Cartesian grids. [START_REF] Kalitzin | Turbulence modeling in an immersedboundary RANS method[END_REF] When the Spalart-Allmaras model is used along with the power-law model defined in Sec. IV B, u τ is calculated using Eq. ( 33) after the velocity is calculated at the boundary node. The value of χ is then obtained with the relation χ = κy + at the boundary node according to the definition of ν in the Spalart-Allmaras model presented in Sec. II B. This boundary condition is thus consistent with the SA model.

V. APPLICATION TO THE TURBULENT FLOW AROUND A 2D AIRFOIL

To validate the proposed power-law wall model and its numerical implementation, the flow around a NACA 23012 airfoil at a chord-based Reynolds number Re = 16 × 10 6 and a Mach number Ma = 0.18 is computed. Three angles of attack, namely, α = 0 • , 5 • , and 10 • , are considered. Calculations are performed using the baseline and negative SA model presented in Sec. II B. The dependency versus the SA variants is investigated in the Appendix. For the LBM simulation, the speed of sound is c s = 340 m/s and the kinematic viscosity is ν = 1.5 × 10 -5 m 2 /s.

A. Calculation setup

The present simulations are carried out using a timemarching method that converges toward a steady state for the studied range of angles of attack.

The computational domain is illustrated in Fig. 4. Its size in the spanwise direction (z) corresponds to one cell size at the coarsest level of refinement which is equal to 8.7%c, where c is the airfoil chord. Periodic boundary conditions are used in that direction. The simulation can thus be considered as a pseudo-two-dimensional case.

The free-stream velocity is imposed at the inlet of the computational domain with an angle corresponding to the angle of attack. A constant density condition is prescribed at the outlet.

Sponge layers are defined according to Xu and Sagaut 53 at the inlet and outlet with, respectively, a fixed velocity and fixed density to prevent spurious reflection of disturbances and ensure convergence toward a steady state. It consists of adding a source term to the LBM equation which is defined as follows:

F i = σ SL (x)(Φ(t + ∆t) -Φ target ), ( 39
)
where Φ is the considered physical variable that should reach the value Φ target in the sponge layer. In this work, target values of density and velocity are defined as their corresponding farfield values. σ SL is the absorbing strength defined as follows:

σ SL (x) = 3125(L sponge -x)(x -x 0 ) 4 256(L sponge -x 0 ) 5 , ( 40
)
where L sponge is the width of the sponge layer and x 0 is its starting position. The solution field in the sponge layer can thus be expressed as

Φ(t + ∆t) = Φ(t + ∆t) -σ SL (x)(Φ(t + ∆t) -Φ target ). ( 41
)
L sponge corresponds to four chord lengths from the inlet and the outlet in this study. Moreover, for the case with an angle of attack equal to 0 • , the upper and lower boundaries are defined as frictionless walls with sponge layers on the density of two chord widths.

For the other angles of attack, the lower surface is defined as an inlet with inclined velocity and the upper surface is an outlet condition with imposed density. For the SA model, χ = 0 is imposed on the airfoil surface. The free-stream value χ = 2 is imposed at the inlet to ensure a fully turbulent calculation.

B. Verification and validation

A grid convergence study is performed first, and simulations are validated against Xfoil results. Xfoil is a solver written by Drela [START_REF] Drela | XFOIL 6.94 User Guide[END_REF] dedicated to the design and analysis of 2D subsonic airfoils. The viscous formulation for the calculation of the boundary layer and wake available in Xfoil is used in this work. It is based on two equations from the integral boundary layer formulation. Results in Xfoil are generated for the same Reynolds and Mach numbers as for the LBM simulation. Furthermore, Xfoil uses the envelope e n transition criterion to determine where laminar-turbulent transition occurs on the airfoil depending on the amplification factor of the frequency that causes transition. For this purpose, the userspecified parameter n crit , linked to the log of this amplification factor, should be adapted in Xfoil depending on the disturbance level of the flow around the airfoil. Since the Xfoil results are used to validate the RANS calculations in fully turbulent mode, a small value of n crit , n crit = 0.001, is used to bypass the transition mechanism in Xfoil, as recommended in the documentation. [START_REF] Drela | XFOIL 6.94 User Guide[END_REF]

Computational grids

A multi-domain uniform grid is used to refine the grid around the airfoil. Two types of refinement are used: rectangular boxes and offsets from the airfoil surface to efficiently refine the grid in the near-wall region, as shown in Fig. 4. At each refinement level, the grid size is halved in each direction.

Three different computational grids have been designed to evaluate the grid influence on the results (see Table I). Grid A is the one presented in Fig. 4. Grid B is based on grid A with additional refinement at both the leading and trailing edges by addition of an additional offset on 10% of the chord at both edges (see Fig. 5). Grid C consists of grid A with grid refinement in the near-wake region over two chord lengths, as shown in Fig. 5 for α = 0 • . This refinement region is adapted for each angle of attack, following the wake angle behind the airfoil.

The simulations are performed over about 17 flow through times to ensure the convergence toward a steady state. This corresponds to 240 000 iterations on grids A and C and 480 000 iterations on grid B. The time step is indeed twice smaller for the calculation on grid B due to the additional refinement level. The convergence of the aerodynamic lift C L and drag C D coefficients is shown in Figs. 678 As can be seen in Fig. 10, regions where the boundary nodes lie in the viscous sublayer (y + < 11.81) are very sparse. Moreover, y + Ref is generally higher than 100 except at the leading and trailing edges, near-stagnation points, where the use of a wall model is not justified. Note that the same observations have been made for α = 5 • and 10 • . A model for the viscous sublayer is thus not necessary, and only the power-law model of Eq. ( 30) is used.

Aerodynamic coefficients

The aerodynamic performance of an airfoil is classically evaluated using the lift (C L ) and drag (C D ) coefficients calculated as follows:

C L = L 1 2 ρ ∞ V 2 ∞ ce , C D = D 1 2 ρ ∞ V 2 ∞ ce , ( 42
)
where c is the airfoil chord, e is the span of the airfoil considered in the simulation, and V ∞ and ρ ∞ are, respectively, the free-stream velocity and density. L and D are the lift and drag forces applied on the airfoil, respectively, with

L = F a • e Y , ( 43
)
D = F a • e X , ( 44
) where e X and e Y are the unit vectors, respectively, in the direction parallel and normal to the free-stream velocity vector. F a is the sum of the pressure and viscous forces applied on the airfoil,

F a = S a -p ndS + S a ( τ : n)dS, ( 45
)
where τ is the viscous stress tensor and n is the outward pointing unit vector normal to the airfoil surface S a . Different methods for the force calculation applied on a body in a fluid are presented in the study of Van Dam. [START_REF] Van Dam | Recent experience with different methods of drag prediction[END_REF] The surface integration consists in calculating the force directly from Eq. ( 45). However, when using Cartesian grids, the pressure and stress tensor are calculated at the boundary nodes which do not systematically lie on the airfoil surface. These quantities have to be extrapolated on the vertices of the surface grid of the airfoil to perform the integration in Eq. ( 45). The force acting on a body can also be obtained with far-field integration where this force is derived from the global conservation law of momentum applied to a control volume encompassing the body. According to Chao and Van Dam, [START_REF] Chao | Airfoil drag prediction and decomposition[END_REF] the accurate prediction of the drag using surface integration necessitates an accurate prediction of the pressure distribution on the airfoil surface. They compared both methodologies of calculation, surface and far-field integration, in a RANS simulation on a body-fitted grid. They obtained fairly good agreement between surface and far-field integration with very well predicted pressure profile. However, the extrapolation performed for the surface integration when using Cartesian grids may lead to oscillations of the stress tensor and pressure distribution on the surface and thus to errors in the surface integration. [START_REF] Mei | Force evaluation in the lattice Boltzmann method involving curved geometry[END_REF] In this work, the far-field integration is used to preclude the use of extrapolated values on the airfoil surface. This estimation method of the drag is thus not based on interpolated values on the airfoil surface.

A control volume V c is defined around the airfoil. Contrary to Chao and Van Dam, [START_REF] Chao | Airfoil drag prediction and decomposition[END_REF] the inlet and side surfaces were not moved to infinity to restrict the far-field integration to the outlet surface only. All surfaces of the control volume are considered in the present work. The inlet surface, denoted by index in, is located approximately one chord upstream of the airfoil when the outlet surface, denoted by index out, is put 3.5 chords downstream the airfoil. The upper (index 1) and lower (index 2) surfaces are defined at about one chord from the airfoil. The conservation law of momentum applied to the control volume V c delimited by surface S = S in + S 1 + S 2 + S out for the RANS steady-state solution is (46) where n is the outward pointing normal to the surface S, u is the mean velocity vector, p is the mean pressure, and τ is the viscous stress tensor. The force acting on the airfoil is thus

S ρ u( u• n)dS = Σ F = -F a + S -p ndS + S ( τ : n)dS,
F a = - S ρ u( u• n)dS + S -p ndS + S ( τ : n)dS. (47)
It is worth noting that in Xfoil, the lift coefficient is calculated by surface pressure integration, whereas the drag coefficient is calculated using the Squire-Young formula at the last point in the wake. A comparison of drag coefficients calculated using the Squire-Young formula or surface integration in RANS body-fitted simulations of flows around airfoils is presented in the study of Coder and Maughmer. [START_REF] Coder | Numerical validation of the Squire-Young formula for profile-drag prediction[END_REF] Very good agreement was found between both methods, and the error was within 2%-3% in the low drag range.

Figure 11 compares the lift and drag coefficients obtained by far-field integration for the three angles of attack on the three grids compared with Xfoil results. Very similar results are obtained on grids A and C. This means that the grid refinement in the wake has no effect on the aerodynamic coefficients' prediction for the considered cases. The prediction of C L is good for α = 0 • and 5 • with the three grids. For α = 10 • , C L is underestimated with grids A and C and better predicted with grid B. Similarly, grid refinement at the leading and trailing edges improves the drag coefficients' prediction for α = 0 • and 5 • which are well predicted with grid B. C D is yet slightly overestimated at α = 10 • with grid B. However, the validity of the Xfoil results, and of the Squire-Young formula, at this higher angle of attack may be called into question. Note that we have checked that the estimation of C L and C D is not dependent on the size of control volume V c .

Pressure coefficient

For a finer validation of the numerical results, the pressure coefficient distributions obtained with the present method are compared with Xfoil results for the three angles of attack under consideration. The pressure coefficient C p is calculated as follows:

C p = p -p ∞ 1 2 ρ ∞ V 2 ∞ , (48)
where p ∞ is the free-stream pressure. As explained in Sec. V B 2, the values of p calculated at the boundary nodes have to be extrapolated on the vertices of the surface grid of the airfoil in order to calculate the C p coefficient. The pressure coefficient distributions obtained with grids A, B, and C are compared with Xfoil results in Fig. 12. They are almost identical which is in agreement with the observation made in Sec. V B 2. Large variations of C p are observed at the trailing edge on grids A and C which are underlined in Fig. 12 for α = 0 • . This is reduced by refining the grid in this region as shown by the results on grid B. The airfoil has indeed a blunt trailing edge which may lead to significant variations of C p if the grid is not fine enough in this region. For α = 0 • , the pressure profile C p on grid B is in agreement with the one predicted by Xfoil. For α = 5 • , the pressure coefficient is well predicted with the three grids despite the slight underestimation of the maximum of -C p . The prediction of C p for α = 10 • is improved near the trailing edge on grid B, as for the other angles of attack, but also near the leading edge.

Small wiggles are observed on the C p distributions, in particular, at α = 0 • . This is linked to the boundary treatment when using Cartesian grids and has been observed in several studies using the cut-cell method [START_REF] Berger | Progress towards a Cartesian cut-cell method for viscous compressible flow[END_REF][START_REF] Hu | Towards efficient viscous modeling based on Cartesian methods for automated flow simulation[END_REF] and immersed boundary techniques [START_REF] Zhou | RANS simulation of high-Re turbulent flows using an immersed boundary method in conjunction with wall modeling[END_REF] in the Navier-Stokes framework or LBM. [START_REF] Li | Non-body-fitted Cartesian-mesh simulation of highly turbulent flows using multi-relaxation-time lattice Boltzmann method[END_REF][START_REF] Mei | Force evaluation in the lattice Boltzmann method involving curved geometry[END_REF] In particular, Hu et al. [START_REF] Hu | Towards efficient viscous modeling based on Cartesian methods for automated flow simulation[END_REF] have shown that these oscillations are linked with the interpolation method used for the solid boundary treatment. In order to avoid such numerical instabilities, (38)] and its version modified as in the work of Tamaki et al. [START_REF] Tamaki | Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method[END_REF] [Eq. ( 50)].

Capizzano [START_REF] Capizzano | Turbulent wall model for immersed boundary methods[END_REF] proposed to assume a linear behavior of the velocity near the wall when using a Cartesian grid. Based on this idea, Tamaki et al. [START_REF] Tamaki | Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method[END_REF] modified the use of the wall model to obtain a linear tangential velocity profile between the boundary node and the reference (Ref ) point. The velocity at the boundary node is thus calculated according to

u + = u + Ref - ∂f ∂y + (y + Ref ) (y + Ref -y + ), ( 49
)
where f is the function of the wall model. This modification has been applied to the power-law wall model used in this work with f (y + ) = A(y + ) B . Using Eq. ( 49), the calculation of the velocity at the boundary node using Eq. ( 38) is replaced by

u(y) = u Ref + B u Ref y Ref (y -y Ref ). (50)
The original wall model [Eq. (38)] and the modified wall model according to Tamaki et al. [START_REF] Tamaki | Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method[END_REF] [Eq. (50)] are compared in Fig. 13. With the modified wall model, the velocity profile is indeed linear between the wall (y = 0) and the Ref point (y = 2.5dx). For a given distance to the wall y o of the boundary node, the velocity obtained with the modified wall model u m is higher than that with the original wall model u o .

The pressure coefficient obtained with grid B using this modified wall model is presented for α = 0 • in Fig. 14. A smoother profile is indeed obtained with this method. The same observation has been made for α = 5 • and 10 • . Note that the same results can also be obtained by smoothing the wall distance of the boundary nodes in the wall model equation. However, as shown in Fig. 15, the aerodynamic coefficients obtained are not correct, and the drag coefficient is in particular overestimated. Smooth pressure profiles can thus be obtained on Cartesian grids but at the expense of the force prediction.

In order to evaluate the difference between the original and the modified model, the distance y m is defined in Fig. 13. y m corresponds to the distance to the wall of a fictitious node where the velocity would be u m by application of the original wall model. Applying the modified wall model to boundary nodes at distance y o to the wall is thus equivalent to applying the original wall model to boundary nodes at distance y m to the wall. Figure 16 flow physics is thus modified by supposing a linear velocity profile in the boundary layer or smoothing the wall distance. In the latter case, the resulting solution can be interpreted as the solution around a slightly corrugated profile which is obtained by moving each point of the profile by a distance equal to the difference between the modified y m and the real y o distance to the wall. A portion of the resulting modified profile is shown in Fig. 17. This results in a rough-like profile which may explain the higher drag obtained with the modified wall model.

As a conclusion, good results are obtained on the prediction of integral values of C L and C D and of the pressure profile C p with proper grid definition at the leading and trailing edges. Grid refinement in the wake does not appear to have a significant influence in the present cases. The far-field integration on grid B yields an estimation of the aerodynamic coefficients that is in very satisfactory agreement with Xfoil results. Moreover, as shown in the Appendix, the power-law model used in this work is robust in regard to different versions of the Spalart-Allmaras model.

C. Application to a lower Re number case

The method presented in Secs. V A and V B is now applied to a lower Reynolds number case in order to further validate the power-law model. The Reynolds number is Re = 1.8 × 10 6 which corresponds to the experiment of Broeren et al. [START_REF] Broeren | Validation of 3-D ice accretion measurement methodology for experimental aerodynamic simulation[END_REF] In order not to change the Mach number, Ma = 0.3, the Reynolds number is modified by scaling the airfoil chord. At this Reynolds number, grid A defined in Sec. V B leads to y + values between 0 and 80 with values principally lower than 50 for x/c > 0.5. These values are low and may not be adapted for the use of the power-law wall model defined in Eq. (38). A new grid, grid A0, is defined by removing one refinement offset around the airfoil, thus leading to y + values mainly between 20 and 120.

The aerodynamic coefficients obtained by far-field integration for several angles of attack are compared in Fig. 18 with the Xfoil results and with the experimental measurements published in the study of Broeren et al. [START_REF] Broeren | Validation of 3-D ice accretion measurement methodology for experimental aerodynamic simulation[END_REF] The prediction of C L with grids A and A0 is very similar and in agreement with Xfoil results and is slightly overestimated compared with the experimental measurements for high angles of attack. The drag coefficient C D obtained with Xfoil is overestimated compared to the experimental measurements. At this lower Reynolds number, Xfoil seems to be less reliable. For angles of attack higher than 5 • , C D is well predicted with grid A0, whereas it is overestimated with grid A. For α ≥ 5 • , the y + value with grid A is mostly lower than 50 over the airfoil, whereas it has been shown in Sec. IV A that the power-law model is valid for y + ≥ 50. This may explain why better results are obtained with a coarser mesh at high angles of attack. The experimental measurements present an almost constant value of C D for α ≤ 5 • . This plateau is not reproduced with the calculations on grids A and A0 in which C D is decreasing with the angle of attack. Drag values closer to the experiment are thus obtained with grid A as a consequence of the overestimation of C D with this grid at high angles of attack. Note however that the C D values for α ≤ 5 • are very low so that a small prediction error may seem important since it is of the same order of magnitude as the C D value.

The experimental measurements of the pressure coefficient C p are available for α = 6.2 • and α = 9.3 • . Comparison of the C p profiles obtained on grids A and A0 with the Xfoil results and experimental measurements is shown in Fig. 19. First, the Xfoil results and experimental measurements are in agreement except at the trailing edge where -C p obtained with Xfoil is slightly lower than that experimentally measured. Calculations on grids A and A0 lead to very similar results which are close to the experiment expect for the slight underestimation of the maximum of -C p .

In conclusion, good results are once more obtained at this lower Reynolds number, thus further validating the use of the power-law model for the flow prediction around airfoils.

VI. CONCLUDING REMARKS

An explicit wall model based on a power-law velocity profile for the description of the turbulent boundary layer is proposed for high Reynolds number flow RANS simulations. This model is particularly well suited for the near-wall treatment using Cartesian grids since the evaluation of the friction velocity is not required for the reconstruction of the velocity at the boundary nodes. Moreover, contrary to most of the commonly used wall models that rely on a logarithmic velocity profile, no iterative procedure is needed to calculate the friction velocity if necessary for the computation of the boundary condition of a RANS model equation or for post-processing purpose.

LBM-based RANS simulations of the flow around an airfoil at high Reynolds numbers have been performed to assess the proposed model and its numerical implementation on uniform grids. The Spalart-Allmaras turbulence model has been used. Verification and validation steps have been conducted in order to assess the reliability of the numerical simulations. A grid convergence study based on three grids with refinement in regions of interest, specifically the leading and trailing edges and the wake of the airfoil, has been carried out at Re = 16 × 10 6 . Grid refinement near the leading and trailing edges, where stagnation points occur, has a significant effect on the accuracy of the results. Numerical simulations at Re = 16 × 10 6 and Re = 1.8 × 10 6 are validated against results from Xfoil and experimental measurements looking at integral quantities, i.e., C L and C D aerodynamic coefficients, and a local quantity, the C p distribution.

Very satisfactory results are obtained at both Reynolds numbers for angles of attack ranging from 0 • to 10 • . The methodology presented in this work, namely, the LBM-based RANS simulation with the power-law wall model, is thus validated for the prediction of attached aerodynamic flows at a high Reynolds number.

The power-law model is thus promising for simplified near-wall treatment in LBM with Cartesian grids. Spalart proposed another way, compared to the negative SA model (see Sec. II B), to handle negative values of the modified vorticity S which was not published but is referred to as the "f v 3 " formulation in the study of Allmaras et al. [START_REF] Allmaras | Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[END_REF] and is presented in the study of Deck et al. [START_REF] Deck | Development and application of Spalart-Allmaras one equation turbulence model to threedimensional supersonic complex configurations[END_REF] The model is the same as the equation for χ > 0 given in Sec. II B but with a modification of the definition of the modified vorticity in order to enforce S ≥ 0,

S = f v3 Ω + ν κ 2 d 2 w f v2 , with f v2 = 1 + χ c v2 -3 , f v3 = (1 + χf v1 )(1 -f v2 ) χ , (A1)
with c v 2 = 5. All other constants remain the same as in the original model. Moreover, to prevent numerical problems, Deck et al. [START_REF] Deck | Development and application of Spalart-Allmaras one equation turbulence model to threedimensional supersonic complex configurations[END_REF] proposed to clip χ by taking max( χ, 10 -4 ). This replaces the use of Eq. ( 12) since χ is then never negative.

According to Rumsey et al., [START_REF] Rumsey | CFD sensitivity analysis of a modern civil transport near buffet-onset conditions[END_REF] a delay of the laminarturbulent transition is observed with this f 3 version compared to the original SA model even without trip terms (fully turbulent simulation) for Reynolds numbers from one million to 10 × 10 6 .

b. Wall-distance-free version of the SA model

The original SA model relies on the distance to the nearest wall in the destruction term. The evaluation of this distance may be difficult to compute and not univocally defined for realistic geometries, as illustrated in the study of Weiss and Deck. [START_REF] Weiss | On the coupling of a zonal body-fitted/immersed boundary method with ZDES: Application to the interactions on a realistic space launcher afterbody flow[END_REF] Menter [START_REF] Menter | Eddy viscosity transport equations and their relation to the k-ε model[END_REF] transformed the two-equation kmodel into a one-equation model in order to establish a connection between one-and two-equation models. He thus obtained a local equation without dependence on the distance to the wall. Rahman et al. [START_REF] Rahman | Wall-distancefree version of Spalart-Allmaras turbulence model[END_REF] only the wall-distance-free feature is retained, and the transport equation is

∂ ν ∂t + ∂ νu i ∂x i = C b S ν Production + ∂ ∂x i ν + ν t σ ∂ ν ∂x i Diffusion -C w ∂ ν ∂x i 2 Destruction . (A2)
The wall-distance dependency is removed, and the first derivative of ν is introduced in the destruction term. The diffusion term is also modified since the term involving the c b2 constant is removed. Moreover, the diffusion term involves the eddy viscosity ν t as in the equation obtained by Menter. [START_REF] Menter | Eddy viscosity transport equations and their relation to the k-ε model[END_REF] In this work, the damping functions f v 1 , f v 2 , and f v 3 , as well as the definition of S, are unchanged according to the SA-f v 3 model. For the numerical resolution of Eq. (A2), ν t is replaced by its expression ν t = f v1 ν in the diffusion term. The constant C b in the production term is equated to the constant c b1 of the original model, C b = c b1 = 0.1355. The formulation proposed by Rahman et al. [START_REF] Rahman | Wall-distancefree version of Spalart-Allmaras turbulence model[END_REF] is used for calculating C w ,

C w = C b C µ (1 + C b ), (A3)
where C µ = 0.09. This leads to C = 1.4. The other constants are the same as in the original model.

Results

Results obtained with these two variants of the SA model are now compared with those obtained with the baseline and negative SA model to check the robustness of the proposed wall modeling approach. In this section, "SA baseline" refers to the variant used in this work and presented in Sec. II B, while "SA fv3" and "SA fv3 wall distance free" refer to variants presented in Subsections 1 a and 1 b of the Appendix, respectively. Calculations are conducted for the high Reynolds number case Re = 16 × 10 6 with grid B presented in Sec. V B.

a. Convergence toward a steady-state solution

The convergence of the aerodynamic coefficients with the three variants of the SA turbulence model is compared in Figs. 20-22. The speed of convergence and the converged values are observed to be identical with the three variants.

b. Aerodynamic coefficients and pressure distribution

The aerodynamic coefficients C L and C D are presented in Fig. 23, along with the C p distributions in Fig. 24. All variants lead to very similar results irrespective of the angle of attack.

The three studied variants of the Spalart-Allmaras turbulence model lead to almost identical results, proving the robustness of the methodology proposed in this paper.

FIG. 1 .

 1 FIG. 1. Wall treatment strategy used for the LBM calculation; : boundary node of interest, and : neighboring fluid nodes of , ◽ and : other boundary nodes, and •: reference point for the application of the wall model; the dashed line surrounds the fluid nodes used to calculate the velocity at the reference point Ref.

FIG. 4 .

 4 FIG. 4. Computational grid. (a) Computational domain; c is the airfoil chord. (b) Refinement offsets around the airfoil.

FIG. 5 .

 5 FIG. 5. Grid refinement (a) at the leading edge, (b) at the trailing edge, and (c) in the wake for the angle of attack α = 0 • .

FIG. 7 .

 7 FIG. 7. Convergence history of lift coefficient C L (a) and drag coefficient C D (b) for the angle of attack α = 5 • .

  FIG. Comparison of the y + and y + Ref values for (a) grid A and (b) grid B on the airfoil for the angle of attack α = 0 • .

Figure 9

 9 Figure 9 displays the y + values at boundary nodes and y + Ref at the reference points in the case α = 0 • . The wall model is indeed applied at both types of points. y + at boundary nodes is highly fluctuating due to the staircase-like form of the grid around the airfoil. This value varies around 200 with maximum values locally reaching 500 at the leading edge for grid A. On grid B, the y + value at boundary nodes is mostly lower than 200 at the trailing and leading edges and is the same as for grid A elsewhere. At the reference point, y + Ref depends on the local grid size ∆x and the friction velocity since y + Ref = 2.5∆xu τ /ν. y + Ref reaches high values especially at the leading edge on grid A. Refining the grid at the leading and trailing edges in grid B divided y + Ref by two in these refinement regions.

FIG. 11 .

 11 FIG. 11. (a) Lift coefficient C L and (b) drag coefficient C D calculated with far-field integration for grids A, B, and C.

FIG. 13 .

 13 FIG. 13. Comparison of the velocity profiles in the boundary layer obtained with the original model and with the model modified according to Tamaki et al.; 2 u o and u m are, respectively, the velocities obtained with the original and the modified wall models for a boundary node at a distance y o from the wall; y m is the distance to the wall required to obtain the velocity u m with the original wall model.

FIG. 14 .

 14 FIG.14. Pressure coefficient C p at the angle of attack α = 0 • with grid B obtained with the power-law wall model modified as in the work of Tamaki et al.[START_REF] Tamaki | Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method[END_REF]

FIG. 15 .

 15 FIG. 15. (a) Lift coefficient C L and (b) drag coefficient C D calculated with far-field integration for grid B with the original power-law wall model [Eq.(38)] and its version modified as in the work of Tamaki et al.[START_REF] Tamaki | Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method[END_REF] [Eq. (50)].

  FIG. 16. Comparison of the y + o and y + m values for grid B at the angle of attack α = 0 • .

FIG. 18 .

 18 FIG.18. (a) Lift coefficient C L and (b) drag coefficient C D calculated with far-field integration for grids A and A0 compared with the experimental measurements (Exp.) of Broeren et al.[START_REF] Broeren | Validation of 3-D ice accretion measurement methodology for experimental aerodynamic simulation[END_REF] and Xfoil results.

FIG. 19 .

 19 FIG. 19. Pressure coefficients C p for grids A and A0 at the angles of attack α = 6.2 • (a) and 9.3 • (b).

1 .

 1 Variants of the SA model a. SA model with the f v3 term (SA-f v3 )

  FIG. 20. Convergence history of lift coefficient C L (a) and drag coefficient C D (b) for the angle of attack α = 0 • .

FIG. 22 .

 22 FIG. 22. Convergence history of lift coefficient C L (a) and drag coefficient C D (b) for the angle of attack α = 10 • .

FIG. 24 .

 24 FIG. 24. Pressure coefficients C p obtained with the three variants of the SA model at the angles of attack α = 0 • (a), 5 • (b), and 10 • (c).

TABLE I .

 I Characteristics of grids A, B, and C.

			Grid B (leading and	Grid C
		Grid A trailing edges refined)	(wake refined)
	Smallest cell size (%c) 0.07	0.035	0.07
	Number of nodes (in million)	10	15.8	23-26 depending on

ACKNOWLEDGMENTS

This work was performed using HPC resources from GENCI-TGCC/CINES (Grant No. 2017-A0012A07679). This work was supported by the French project CLIMB, with the financial support of BPIFrance (Project No. P3543-24000), in the framework of the program "Investissement d'Avenir: Calcul Intensif et Simulation Numérique." This work was carried out using the ProLB solver.