Surface Roughness in RANS Applied to Aircraft Ice Accretion Simulation: A Review

: Experimental and numerical ﬂuid dynamics studies highlight a change of ﬂow structure in the presence of surface roughness. The changes involve both wall heat transfer and skin friction, and are mainly restricted to the inner region of the boundary layer. Aircraft in-ﬂight icing is a typical application where rough surfaces play an important role in the airﬂow structure and the subsequent ice growth. The objective of this work is to investigate how surface roughness is tackled in RANS with wall resolved boundary layers for aeronautics applications, with a focus on ice-induced roughness. The literature review shows that semi-empirical correlations were calibrated on experimental data to model ﬂow changes in the presence of roughness. The correlations for RANS do not explicitly resolve the individual roughness. They principally involve turbulence model modiﬁcations to account for changes in the velocity and temperature proﬁles in the near-wall region. The equivalent sand grain roughness (ESGR) approach emerges as a popular metric to characterize roughness and is employed as a length scale for the RANS model. For in-ﬂight icing, correlations were developed, accounting for both surface geometry and atmospheric conditions. Despite these research efforts, uncertainties are present in some speciﬁc conditions, where space and time roughness variations make the simulations difﬁcult to calibrate. Research that addresses this gap could help improve ice accretion predictions.


Introduction
Fluid dynamics studies, both experimental and numerical, intend to understand the fluid flow inside or above complex geometries.Among the applications often encountered are flows in pipes [1], flows around aerial or terrestrial vehicles [2,3] or flows and wind in urban zones [4].On the computational fluid dynamics (CFD) side, the simplest approach to model the walls and solid parts of the geometry is to assume smooth surfaces [5,6].The smooth approach has its limits.In practice, several geometries have irregular surfaces with roughness of various origins.Wear, deposit, for example in turbines [7], and ice accretion on an aircraft surface [8] cause roughness.Experimental studies showed different flow structures and behaviour between rough and smooth surfaces [9].The influence of roughness is mainly concentrated in the near-wall region and affects the turbulent boundary layer [10], particularly the inner layer.The boundary layer modification increases flow metrics such as the wall heat transfer and the skin friction [11].Experimental studies have tackled this for decades [12,13], highlighting heat transfers, and skin friction increases in the presence of roughness.Therefore, the mathematical models used in CFD should account for roughness to produce accurate predictions.
CFD simulations can either resolve the flow around individual roughness elements or model the effects of roughness.Direct Numerical Simulation (DNS) is one possible way to resolve a flow around surface roughness elements.DNS became popular over the last decade, with the increase in computational power [14,15].Large Eddy Simulation (LES) also resolves individual roughness elements, but with a lower computational cost [16,17].Nevertheless, the computational cost of DNS and LES limits their application to simple geometries and is impractical for industrial problems.The modeled approaches present better applicability and wider ranges of use compared to the resolved approaches.The Reynolds Averaged Navier-Stokes (RANS) is popular in an industrial context since it does not require as much computational power as DNS or LES [18].
The RANS models the effects of turbulent fluctuations with an additional stress term in the Navier-Stokes equations.This stress term is usually assumed proportional to the wall normal velocity gradient, and the proportionality factor is called the eddy viscosity.The eddy viscosity turbulence models, such as the Spalart-Allmaras (SA) model [19] or the k-ω model [20], compute the eddy viscosity by solving one (SA) or two (k-ω) transport equations.In a context of rough geometries with RANS, these turbulence models are commonly adapted to account for roughness [21,22].Experimental skin friction obtained using manufactured roughness patterns help calibrate the rough turbulence models [12].Sometimes, the rough turbulence models over predict the wall heat fluxes [23].This inaccuracy led to the development of so-called thermal correction models aiming to reduce the heat fluxes.Morency and Beaugendre [23], Chedevergne [24], Suga et al. [11], and Aupoix [25] derived thermal corrections to use with the rough turbulence models.Most of the thermal correction models increase the turbulent Prandtl number to reduce the wall heat fluxes.The correction depends on the roughness geometric characteristics.Apart from the rough turbulence models, the discrete element roughness method (DERM) is another modeled approach used in CFD [26,27].The DERM directly adds the roughness effects in the mass, momentum, and energy equations of the flow.The DERM adds source terms to the flow equations to account for blockage, drag, and heat transfer induced by roughness [28].Chedevergne recently used the DERM helped by experimental surface topography measurements to finely quantify the source terms in the model [29].The intensity of the additional source terms varies with the roughness geometric characteristics.
Experimental setups can characterize the roughness geometries.Modern measurement techniques allow a precise mapping of experimental surfaces, for example with laser measurements [30].Particle velocimetry methods are able to measure the velocity around rough surfaces [31], helping in the identification of typical flow structures.Roughness measurements show highly irregular roughness patterns in most engineering applications [32].These irregular patterns are usually non-trivial to input in the CFD models, since uncertainty persists.A roughness pattern requires quantifiable metrics to model it with CFD.Common metrics statistically characterize the height and/or spacing of the roughness elements on a surface [33], and are used to describe a given roughness pattern.For RANS turbulent flow description, one popular metric is the equivalent sand grain roughness (ESGR) height, k s , introduced almost one century ago by Nikuradse [34].The ESGR height allows the definition of a roughness pattern with a single metric.k s is the diameter of the packed spheres that would give the same skin friction as the actual real roughness pattern [12], obtained experimentally or geometrically resolved with DNS.k s is relevant for random and regular roughness patterns similar to sand grain paper.For regular roughness patterns, correlations related k s to some geometrical parameters of the roughness elements [35,36].However, this single metric can fail to predict the skin friction for highly irregular roughness patterns, for example, in a context of ice-induced roughness on an aircraft [37].
Aircraft icing is a specific application of flow over a rough surface.At the early stages, ice accretion creates surface roughness that increases the wall heat transfer and modifies the final ice shape [38].For CFD-based ice accretion simulations, the usual framework employs a RANS airflow solver, a droplet impingement solver, and an accretion solver [39].The RANS airflow solver must model surface roughness to increase the wall heat transfer and the wall skin friction [8].k s is used to characterize the roughness pattern in most ice accretion software, for example in FENSAP-ICE [40], IGLOO3D [41], or LEWICE [42].Ice-specific semi-empirical correlations are developed to relate the roughness pattern to the Fluids 2023, 8, 278 3 of 21 atmospheric conditions and the airfoil geometry [43][44][45][46].In most ice accretion software, the surface roughness imposed in the model is constant and uniform all over the surface [47], even if experimental observations show non-uniform roughness patterns over the airfoil surface [43,48,49].A recent study followed these observations and derive non-uniform roughness distributions [50] to obtain a heat transfer distribution closer to the experiment.Research efforts are still needed to improve the roughness models as some ice accretion simulations fail at predicting a satisfactory ice shape compared to the experiment [51].
This review paper investigates both RANS approaches with wall resolved boundary layers (denoted as wall resolved RANS for the rest of the paper) that account for roughness in aeronautics flows and the special case of ice-induced roughness in aircraft in-flight icing.The objective is to review (1) how roughness is modeled in RANS in the scope of the ESGR and ( 2) what the specific roughness treatments required for aircraft icing are.First, the roughness geometrical characterization is depicted, with a focus on the usual metrics and the ESGR approach.Second, the roughness impact on the turbulent flow structures is depicted, highlighting some rough turbulence model implementations for wall resolved RANS.Finally, the roughness impact on the in-flight aircraft ice accretion is detailed.This last section also highlights the semi-empirical correlations and models for ice-induced roughness description.

Roughness Geometrical Characterization in Aerodynamics
This section gives an insight on how the presence of roughness impacts the flow structure in the boundary layer.The first subsection develops the roughness pattern characterizations classically employed in CFD, such as the statistical metrics.The second subsection focuses on the particular case of the ESGR height, a metric typically used in RANS applications.

Geometrical Parameters
In engineering, the spectrum of roughness patterns encountered varies from regular well-ordered configurations to highly irregular and uncertain geometries such as in ice accretion [52] or heterogeneous surface state [53].Roughness patterns are mathematically described by quantifiable parameters.Through the years, standards have normalized the description of surface states.The SAE J448 [54], the ISO 1302 norm [33,55,56], or the ISO 4287 norm [57] are among the international standards used.The parameters are defined either on a slice of the surface to represent the peaks and valleys amplitudes, or on an area to finely represent the spatial organization of the surface.The main metrics encountered for a surface description are:

•
Ra, the arithmetic mean height; Figure 1 represents a generic rough area and a slice of a rough surface and identifies the most common amplitude metrics.For a rough profile sampled with n points, each one with a height y i from the smooth mean-zero line, the expressions for the previous parameters are given in Equations ( 1)- (7).

𝑅𝑣 = |min (𝑦 )|
(3) =  +  In the cases of 3D roughness patterns covering an area, some authors introduced spatial parameters to account for the spacing and the wetted and/or frontal area of the roughness [58].The shape parameter Λ is a roughness parameter commonly used in the literature.The shape parameter suggested by Dirling [35] is a function of the average roughness element spacing d, the roughness height k, the frontal area of a roughness element Af, and the total windward wetted area of a roughness element As.Equation ( 8) presents the mathematical definition of the shape parameter defined by Dirling.

𝛬
=     In the cases of 3D roughness patterns covering an area, some authors introduced spatial parameters to account for the spacing and the wetted and/or frontal area of the roughness [58].The shape parameter Λ is a roughness parameter commonly used in the literature.The shape parameter suggested by Dirling [35] is a function of the average roughness element spacing d, the roughness height k, the frontal area of a roughness element A f , and the total windward wetted area of a roughness element A s .Equation (8) presents the mathematical definition of the shape parameter defined by Dirling.
Sigal and Danberg [36] introduced a modified definition of the shape parameter by using a surface ratio to account for the roughness element proximity.The total smooth surface S, before roughness addition, and the total frontal area S f are accounted for in the expression (Equation ( 9)).The value of the exponent suggested by Sigal and Danberg was calibrated to match Schlichting's data [59].
The correlations of Dirling and Sigal and Danberg require the geometric characteristics of a single roughness element, such as A f or A s .These correlations are based on regular roughness elements such as spheres, cones, or rods.Therefore, the correlations may not be sufficient for irregular roughness with various geometries.To include such irregular cases, van Rij et al. [60] suggested another variant of the shape parameter definition.The total roughness windward wetted area S s is used in that case (Equation ( 10)).This definition requires detailed representations of the roughness pattern to extract S s .
These geometrical features are useful to calculate an important parameter for the rough turbulent flow: the equivalent sand grain roughness.

The Equivalent Sand Grain Roughness (ESGR)
To quantify the roughness effects, the ESGR constitutes a length scale often used in engineering fields.Nikuradse [34] carried out experiments with sand-grain-like roughness patterns.The roughness pattern was modeled as packed identical spheres whose diameter is called the sand grain roughness (SGR) height.This approach suits regular roughness patterns with a sand-grain texture well.For irregular roughness, Colebrook et al. [61] introduced the ESGR height, often denoted k s .The ESGR height relates non-uniform roughness patterns effects to the SGR height effects on the boundary layer [58].The mathematical definition of the ESGR is obtained knowing the flow velocity U at a normal distance y above the surface [59].u τ denotes the friction velocity and κ is the von Kármán constant.
Without the knowledge of U(y), the ESGR height is correlated to the geometrical configuration of the roughness elements.Dirling [35] established a correlation for cones or spheres elements.Bons [7,62] worked on correlations for turbine blades roughness.Sigal and Danberg [36] derived ESGR correlations for transversal roughness such as bars and rods, while van Rij et al. [60] developed correlations for 3D irregular roughness.Botros [13] recently extended the pioneer studies of Nikuradse and Colebrook for steel pipes.Table 1 presents some correlations used to compute the ESGR k s .

Roughness Regimes and Rough Turbulence Models
This section depicts roughness impacts on the flow structure and reviews how RANS models tackle roughness.The first subsection shows how the presence of roughness in a flow alters the boundary layer, defining the roughness regimes encountered.The second subsection presents some RANS approaches used to include the roughness in the computation, including the velocity and thermal corrections required to model the flow over a rough surface.

The Rough Flow Regimes
The roughness modifies the near-wall flow behavior compared to a smooth surface.Smooth surfaces were extensively studied for almost a century with the definition of the so-called law-of-the-wall [63][64][65].The smooth theory was extended to rough surfaces to face industrial applications.For the rough turbulent boundary layer, a roughness Reynolds k + s number is defined.The symbol ν appearing in Equation ( 12) stands for the kinematic viscosity of air.
Three main roughness regimes can be observed: Qualitatively, these regimes, illustrated in Figure 2, correspond to: • Roughness elements small compared to the viscous sublayer thickness (hydraulically smooth); • Roughness elements in the same order of thickness as the viscous sublayer thickness (transitionally rough); and • Roughness elements larger than the viscous sublayer thickness and ending inside the buffer or the logarithmic layer (fully rough) [10].The roughness effect of the flow field is no more limited to the viscous sublayer [66]. , and  , delimit the regimes and take different values depending on the type of roughness encountered.Nikuradse [34] was a pioneer in the field and suggested values for uniform sand grain roughness.Ligrani and Moffat [67]   k + s,smooth and k + s,rough delimit the regimes and take different values depending on the type of roughness encountered.Nikuradse [34] was a pioneer in the field and suggested values for uniform sand grain roughness.Ligrani and Moffat [67] suggested values for spherical roughness elements.Schultz and Flack: k + s,smooth = 2.5 and k + s,rough = 25 These theoretical and experimental developments are the basis for the CFD models aiming at accounting for surface roughness.The most common implementations used for RANS modeled approaches are reviewed in the following subsection.

RANS Implementations to Account for Roughness
The roughness effects mainly create a velocity deficit in the boundary layer compared to a smooth surface.This leads to a downward shift of the log-Law by a value ∆U + [58], as in Equation ( 13) where C is a constant.
This theoretical shift is confirmed by experimental studies for example from ONERA and Boeing [21], compared to the smooth profile [65], as illustrated in Figure 3.The velocity shift indicates an increase in the friction velocity, and by extension, in the wall shear stress and skin friction coefficient.The amplitude of the velocity shift is a function of the roughness parameters, mainly  .Therefore, RANS requires to accurately model the velocity shift ΔU + to accurately predict the rough skin friction.
To account for roughness in a RANS model, the turbulence models are adapted.A decade after its original publication, the Spalart-Allmaras (SA) turbulence model [19] received a rough extension [21].The original SA turbulence model solves one transport equation for the turbulence variable  .The computation of the turbulent variable enables the knowledge of the turbulent eddy viscosity µt, function of the wall distance, required for the RANS equations closure.The rough version of the SA turbulence model redefines the wall distance from its smooth value d to drough, using the ESGR.

𝑑
=  + 0.03 The velocity shift indicates an increase in the friction velocity, and by extension, in the wall shear stress and skin friction coefficient.The amplitude of the velocity shift is a function of the roughness parameters, mainly k + s .Therefore, RANS requires to accurately model the velocity shift ∆U + to accurately predict the rough skin friction.
To account for roughness in a RANS model, the turbulence models are adapted.A decade after its original publication, the Spalart-Allmaras (SA) turbulence model [19] received a rough extension [21].The original SA turbulence model solves one transport equation for the turbulence variable ν.The computation of the turbulent variable enables the knowledge of the turbulent eddy viscosity µ t , function of the wall distance, required for the RANS equations closure.The rough version of the SA turbulence model redefines the wall distance from its smooth value d to d rough , using the ESGR.
The main modification of the SA turbulence model is the wall boundary condition, which is ν wall = 0 for a smooth case.For the rough version, a gradient normal to the wall, with n being the normal direction, is imposed.
Similarly, the two-equation k-ω turbulence model [20], gained its roughness extensions through the years.Wilcox [70] was the first to derive a rough extension in 1998 before adding some corrections in 2006.The wall boundary condition for the parameter ω is altered.ρ is the fluid density and µ is the dynamic viscosity.
This model has the disadvantage of requiring a very fine mesh in the near wall [71] and is not fully compatible with the original SST formulation [72].Later, Knopp et al. [73] suggested a new formulation for the rough k-ω turbulence model to improve the weaknesses of the original rough extension.Apart from the SA and k-ω turbulence models, the k-ε model also received its rough extension [74].These rough versions predict accurate skin frictions compared to experimental data [23].This allows an evaluation of the friction velocity and of the roughness Reynolds number k + s (Equation ( 12)).Several studies showed that the rough velocity profile and the velocity shift ∆U + (see Figure 3) are functions of the roughness Reynolds number k + s [75][76][77].Equation (18) shows Nikuradse [34] estimation for the rough log-Law in a fully rough regime.
More recently, Radenac et al. [77] suggest using the relation developed by Grigson [76] based on the work of Colebrook [61] to directly estimate the velocity shift.
Another alternative, developed by Kays and Crawford [75], gives a variant to compute the velocity shift.
In all models, the velocity shift is zero in a hydraulically smooth regime.For very small roughness elements, the effects on the flow vanish due to the viscosity [10,78].Commercial CFD software with wall resolved RANS approaches, such as STAR-CCM+, implement the roughness effect such as Equation ( 21) [79].By default, the constants 5.23 and 0.253 in Equation ( 21) are set but tunable by the user, such as the values for k + s,smooth and k + s,rough .Even if the skin friction is accurate, the wall heat transfer is not correctly predicted for some rough cases.Past studies [11,23,25] reported an overestimation of the wall heat fluxes when the rough turbulence models are employed.Similar to the velocity profile, the temperature profile is shifted in the boundary layer [25] above roughness.Kays and Crawford [75] gave a formulation for the temperature shift ∆T + , where Pr t is the turbulent Prandtl number and C Pr is a function of the laminar Prandtl number Pr.
The prediction of the wall heat transfer in the RANS models requires to accurately predict the temperature shift ∆T + .
According to the investigations of Suga et al. [11], the heat flux overestimation comes from the trapped recirculating flow between the roughness elements.These recirculation zones create a thermal barrier reducing the heat flux, as schemed in Figure 4.
The prediction of the wall heat transfer in the RANS models requires to accurately predict the temperature shift ΔT + .
According to the investigations of Suga et al. [11], the heat flux overestimation comes from the trapped recirculating flow between the roughness elements.These recirculation zones create a thermal barrier reducing the heat flux, as schemed in Figure 4.The thermal correction models [11,23,25] intend to account for these recirculation zones.The common approach in the literature is to locally increase the turbulent Prandtl number Prt.This impacts both the thermal conductivity and the temperature shift (Equation ( 22)).The thermal correction models add ΔPrt to the turbulent Prandtl number, allowing the modification of the thermal conductivity.

𝑃𝑟 ,
=  + ∆ Suga et al. [11] derived the following expression for the turbulent Prandtl number correction, where εs and εy are the local turbulent kinetic energy at the ESGR height and The thermal correction models [11,23,25] intend to account for these recirculation zones.The common approach in the literature is to locally increase the turbulent Prandtl number Pr t .This impacts both the thermal conductivity and the temperature shift (Equation ( 22)).The thermal correction models add ∆Pr t to the turbulent Prandtl number, allowing the modification of the thermal conductivity.
Pr t, rough = Pr t + ∆Pr t Suga et al. [11] derived the following expression for the turbulent Prandtl number correction, where ε s and ε y are the local turbulent kinetic energy at the ESGR height and at the height y from the wall, respectively.
C 0 = 5.5 The thermal correction model from Suga et al.only requires one roughness parameter, k s , to compute the thermal correction.
Morency and Beaugendre [23] derived a 2-parameter thermal correction model by blending the works from Ligrani et al. [80], Dipprey and Sabersky [81], and Owen and Thomson [82], and the integral method from Kays and Crawford [75].The general form of ∆Pr t is inspired from the thermal correction model from Aupoix [25].
The expression of F is given by Equation (27).
The factor g damps the correction in a transitionally rough regime.The limits used to separate the hydraulically smooth (g = 0), transitional, and fully rough (g = 1) regimes corresponds to the k + s,smooth = 5 and k + s,rough = 70 values suggested by Nikuradse [34].A third thermal correction model is the 3-parameter model of Aupoix [25].The turbulent Prandtl number correction ∆Pr t is the same as in Equation ( 26), but with an alternate definition for the factor F, function of the velocity shift ∆U + , and a geometrical correction parameter S corr .S corr corresponds to the ratio between the rough wetted area (i.e., accounting for the wetted area of the roughness elements) and the smooth wetted area, with the through below the recirculation zone (see Figure 4) neglected.The factor F is expanded in Equation (28).
The factors A and B are defined in the following Equations ( 29) and (30).
∆U + being a function of k + s , the correction suggested by Aupoix requires three roughness parameters: k s ; k; and S corr .
The RANS models depicted in this subsection predict both wall shear stress and wall heat flux over roughness.The next section gives an insight into a specific application of rough surfaces: aircraft icing.

The Specific Case of Aircraft Icing
In-flight ice accretion simulation commonly relies on RANS and rough turbulence to model airflow.This section first depicts the models used in aircraft ice accretion simulations and how surface roughness plays a key role in the process.The second subsection reviews the roughness semi-empirical correlations developed to determine the roughness pattern above ice surfaces.

The Ice Accretion Process: A Roughness-Dependent Phenomenon
Ice accretion studies address the safety concerns induced by ice build-up on aircraft surfaces [83,84].The studies apply to both ice-induced aerodynamic degradation [85,86] and shedding debris [87,88].Initially carried out experimentally [38,45,89], ice accretion studies have been supplemented by numerical simulations with the increase of computational power [52,90,91].
For numerical simulations, Messinger [92] pioneer works develop a well-known model to compute the ice growth on an aircraft surface.The Messinger model has been improved through the years to include features such as unsteadiness, shear-stress driven liquid film, or partial differential equation formulation [93][94][95][96].The original Messinger model computes the mass and the energy balance in a control volume composed of a moving liquid film and an ice layer, as shown in Figure 5.The white arrows in Figure 5 symbolize the main contributions to the mass (Equation ( 31)) and the energy balance (Equation ( 32)) where  (kg/s) and  (J/s mass flux rate and energy flux rate, respectively.The subscripts correspond to t illustrated in Figure 5.

𝑚 + 𝑚
The accretion can be seen as a phase change problem where the liquid film releasing latent heat during the process.In a control volume, the water phase i mined by the internal energy.If the energy in the control volume is lower than th fication energy, the water is solid.If the energy in the control volume is higher t fusion energy (the solidification energy plus the latent heat of phase change), the liquid.When the total energy of the control volume is between the energy of solidi and the energy of fusion, ice and liquid water coexist, and a glaze ice state is ob The contributions remain almost identical among the ice accretion models, with so ditions such as the radiation or the conductivity through ice [47,97].The variables i .
The accretion can be seen as a phase change problem where the liquid film freezes, releasing latent heat during the process.In a control volume, the water phase is determined by the internal energy.If the energy in the control volume is lower than the solidification energy, the water is solid.If the energy in the control volume is higher than the fusion energy (the solidification energy plus the latent heat of phase change), the water is liquid.When the total energy of the control volume is between the energy of solidification and the energy of fusion, ice and liquid water coexist, and a glaze ice state is observed.The contributions remain almost identical among the ice accretion models, with some additions such as the radiation or the conductivity through ice [47,97].The variables in Equations ( 31) and (32) were described by Özgen and Canibek [98].Among the variables, several depend on quantities obtained at the airflow computation step, such as the heat transfer coefficient htc obtained from the wall heat flux φ, the wall temperature T, and the recovery temperature T rec .
More precisely:

•
The evaporation/sublimation mass and energy variables depend on the heat transfer coefficient [99]; The convective heat loss is a function of the heat transfer coefficient;

•
The runback velocity is a function of the skin friction coefficient in the shear-stress driven liquid film models [94,97].
Both heat transfer and wall shear stress are functions of the surface roughness pattern and impact the computed ice shape.The development of surface roughness in the early stages then drives the ice accretion process [100].This surface roughness is usually uncertain, both in space and time, raising a difficulty to correctly model it.Experimental measurements using photogrammetry and 3D scanning highlight a spatial variation of the roughness, with a smoother zone around the stagnation point [49,101].Additional studies allow the monitoring of the temporal evolution of the roughness characteristics [102].To highlight the spatial roughness variation, the experimental measurements from Han and Palacios [101] were conducted on a NACA0012 airfoil and Baghel et al. [49] used the High Altitude Pseudo Satellites (HAPS) airfoil.Variation of measured experimental Ra (Equation ( 1)) and Rq (Equation ( 2)) against the chord fraction on the upper part of the airfoil are shown in Figure 6.The x-axis origin corresponds to the stagnation point location.Figure 7 highlights the temporal roughness variation obtained on a 0.33 m chord, 0.3 m span NACA0012 wing at −5 • angle of attack [102].Figure 7 plots Ra (Equation ( 1)) against the exposure time for various freestream velocities (V) and total temperatures (T total ).Figures 6 and 7 also give indications on the magnitude of roughness height usually encountered in aircraft ice accretion.This surface roughness only characterizes the very first moments of the ice accretion.The roughness models the surface irregularities of size similar to the boundary layer thickness.The ice then grows on the rough airfoil, driven by the heat transfer and shear stress, and can reach large size and thickness.In the case of ice accretions larger than the boundary layer thickness, strategies such as geometry update [103] or the immersed boundary method [104] are required to account for the wall macroscopic deformation.The multi-layer approach allows updating the airflow during the simulation, as the accretion is computed layer by layer.In most cases, the single layer approach is preferred.The single layer approach only relies on the initial flow field obtained on the rough airfoil without any large ice accretion.In some cases, the single and multi-layer approach can lead to similar results [47].
Figure 8 illustrates the influence of the thermal correction model on the ice shape in numerical simulations.Figure 8 highlights how the changes in the heat transfer coefficient (Equation ( 33)), resulting from the thermal correction chosen, impact the final ice shape on a NACA0012 airfoil in glaze conditions [105].The results are presented for the smooth SA turbulence model [19], the rough SA without thermal correction [21], the rough SA with the Aupoix thermal correction [25], and for the rough SA with Morency and Beaugendre thermal correction [23].
as the accretion is computed layer by layer.In most cases, the single laye preferred.The single layer approach only relies on the initial flow field ob rough airfoil without any large ice accretion.In some cases, the single an approach can lead to similar results [47].    Figure 8 illustrates the influence of the thermal correction model on the ice sh numerical simulations.Figure 8 highlights how the changes in the heat transfer coe (Equation ( 33)), resulting from the thermal correction chosen, impact the final ice on a NACA0012 airfoil in glaze conditions [105].The results are presented for the s SA turbulence model [19], the rough SA without thermal correction [21], the rou with the Aupoix thermal correction [25], and for the rough SA with Morency and gendre thermal correction [23].Given the uncertainty and complex roughness patterns in the experimental setups, empirical correlations were developed to take surface roughness into account in the numerical simulations.These correlations are reviewed in the next section.

Empirical Roughness Correlations in Ice Accretion Simulations
Usual roughness correlations developed for sand grain paper, manufactured roughness elements, and more regular roughness patterns must be adjusted for irregular iceinduced roughness.Ice-specific roughness height and ESGR correlations were developed in an attempt to mimic the experimental observations.The early correlations assumed a uniform roughness pattern for the sake of simplification and kept the same roughness for all atmospheric conditions [106].Gent et al. [107] suggested a roughness height function Given the uncertainty and complex roughness patterns in the experimental setups, empirical correlations were developed to take surface roughness into account in the numerical simulations.These correlations are reviewed in the next section.

Empirical Roughness Correlations in Ice Accretion Simulations
Usual roughness correlations developed for sand grain paper, manufactured roughness elements, and more regular roughness patterns must be adjusted for irregular iceinduced roughness.Ice-specific roughness height and ESGR correlations were developed in an attempt to mimic the experimental observations.The early correlations assumed a uniform roughness pattern for the sake of simplification and kept the same roughness for all atmospheric conditions [106].Gent et al. [107] suggested a roughness height function of the liquid water content (LWC), measuring the mass of liquid water droplet in a given volume of dry air (usually expressed in g/m 3 ), the freestream velocity, the ambient air temperature, and the airfoil chord length c.Ruff [44] suggested a correlation for a uniform ESGR estimation.Shin and Bond [45] extended this correlation to include the dependency to the median volume diameter (MVD) of the water droplets, noticing that this dependency vanishes for MVD below 20 µm.The correlations from Ruff [44] and Shin and Bond [45] take the general form given by Equation (34).The coefficients A, B, C, and D are functions of the LWC, the freestream temperature, the freestream velocity, and the MVD (for [45] only), respectively.
Recent works by CIRA still use the Shin and Bond uniform relation to determine k s [108].The commercial software FENSAP-ICE initially used Shin and Bond uniform correlation [46,109], or a constant user-defined k s .These uniform roughness correlations are still available in the FENSAP-ICE package and are regularly used [110].FENSAP-ICE later received a non-uniform roughness beading model accounting for the roughness variability in space and time [111][112][113].Li et al. [114] also developed a non-uniform roughness distribution in FENSAP-ICE, accounting for the state of the liquid water.In case of a continuous liquid film, the roughness height depends on the film height h w as in Equation ( 35), where g is the gravitational constant, τ w is the wall shear stress, and µ w is the water viscosity.In the case of scattered beads, the gravitational forces and water surface tension drive the roughness height.
Fortin et al. [115] also derived a specific roughness height correlation function of the local wall shear stress and local water film thickness.Anderson and Shin [116] also developed a non-uniform roughness height correlation, based on the local freezing fraction f.The freezing fraction is the proportion of incoming water (runback and impingement) that actually freezes in a control volume.This roughness height correlation, used in the ice accretion software LEWICE [117], is given in Equation (36).
Fortin [43] recently developed a logarithmic correlation for the ESGR height that includes the local collection efficiency β, the water latent heat of fusion L f , the water fusion temperature T 0 , and the recovery temperature T rec alongside the liquid water content LWC (Equation (37)).
To account for the spatial variability of the roughness height, and following the LEWICE correlation (Equation ( 36)), Han and Palacios [50] derived a parabolic roughness distribution with respect to the wrapped distance s from the stagnation point.This correlation, with the smooth width w and the icing limit l as parameters, takes the form of Equation (38).
The parameters w, l, and k max in Equation ( 38) are mainly dependent on the local freezing fraction, the accumulation parameter [50], and the airfoil leading edge radius.The correlation somehow fits the experimental roughness height, as plotted in Figure 9.The recent advances made in ONERA with the IGLOO3D suite combine both experimental measurements and empirical correlation for the roughness estimation [41].IG-LOO3D relies on the thermal correction model from Aupoix [25] (see subsection III.B, Equations ( 28)-( 30)).Scorr and the roughness height k are measured from experimental short-time-exposure accretions and ks is then computed using the semi-empirical correlation from Flack and Schultz [58].The research work at NASA with the GlennICE code adopts a slightly different approach, where the roughness metric Rq is directly computed using the freezing fraction, the pressure coefficient, and the total temperature.The heat transfer coefficient is then augmented by a correction function of Rq [118].
Despite all the research efforts to develop roughness correlations for numerical simulations, existing roughness models fail to predict some experimental ice shapes.The roughness models, as well as the accretion models, are probably not universal and their results depend on both test case and code used [47,51].Data-driven calibration methods and uncertainty quantification techniques are emerging solutions that could establish roughness patterns specifically tailored for a given case/code configuration.Such methods are already used for roughness estimation in other engineering fields [119].Up to recently, only a few published works utilize the methods to calibrate roughness patterns for ice accretion prediction.Although the calibrated roughness patterns lack universality, the methods open new perspectives for future investigations of roughness models and their impact on ice shapes.

Conclusions
This paper reviewed how roughness is tackled in wall resolved RANS and the representation of surface roughness in ice accretion simulations.Focusing on the RANS turbulence models with equivalent sand grain roughness (ESGR) approach, the survey established that specific adjustments are required for realistic predictions of both heat transfer and skin friction over surface roughness.These requirements are justified by the boundary layer flow modification induced by roughness, especially the viscous sub-layer.Semi-empirical correlations were developed through decades to model the roughness geometry and ease its addition in the CFD models.The ice accretion phenomenon is a particular case of a rough surface, where the roughness pattern is highly uncertain, irregular, and The recent advances made in ONERA with the IGLOO3D suite combine both experimental measurements and empirical correlation for the roughness estimation [41].IGLOO3D relies on the thermal correction model from Aupoix [25] (see subsection III.B, Equations ( 28)-( 30)).S corr and the roughness height k are measured from experimental short-time-exposure accretions and k s is then computed using the semi-empirical correlation from Flack and Schultz [58].The research work at NASA with the GlennICE code adopts a slightly different approach, where the roughness metric Rq is directly computed using the freezing fraction, the pressure coefficient, and the total temperature.The heat transfer coefficient is then augmented by a correction function of Rq [118].
Despite all the research efforts to develop roughness correlations for numerical simulations, existing roughness models fail to predict some experimental ice shapes.The roughness models, as well as the accretion models, are probably not universal and their results depend on both test case and code used [47,51].Data-driven calibration methods and uncertainty quantification techniques are emerging solutions that could establish roughness patterns specifically tailored for a given case/code configuration.Such methods are already used for roughness estimation in other engineering fields [119].Up to recently, only a few published works utilize the methods to calibrate roughness patterns for ice accretion prediction.Although the calibrated roughness patterns lack universality, the methods open new perspectives for future investigations of roughness models and their impact on ice shapes.

Conclusions
This paper reviewed how roughness is tackled in wall resolved RANS and the representation of surface roughness in ice accretion simulations.Focusing on the RANS turbulence models with equivalent sand grain roughness (ESGR) approach, the survey established that specific adjustments are required for realistic predictions of both heat transfer and skin friction over surface roughness.These requirements are justified by the boundary layer flow modification induced by roughness, especially the viscous sub-layer.Semi-empirical correlations were developed through decades to model the roughness geometry and ease its addition in the CFD models.The ice accretion phenomenon is a particular case of a rough surface, where the roughness pattern is highly uncertain, irregular, and difficult to measure.Ice accretion simulation codes need to account for ice roughness.Several icing-specific correlations relate roughness parameters to various atmospheric conditions and local ice properties on surfaces.The most complex correlations attempt to predict the roughness variability in space and time.A research gap was highlighted, showing that the roughness correlations may fail to predict correct ice shapes in some cases.This opens the perspective of using uncertainty quantification tools to improve the roughness pattern prediction for ice accretion simulations and detach the roughness estimation from the historical semi-empirical correlations.As a recommendation, other gaps still need attention such as a better understanding of the roughness time evolution.

Figure 1 .
Figure 1.(A) Example of rough area, with fluid flowing above, and (B) a rough profile with the usual metrics.

Figure 1 .
Figure 1.(A) Example of rough area, with fluid flowing above, and (B) a rough profile with the usual metrics.

Figure 2 .
Figure 2. (A) Hydraulically smooth; (B) Transitionally rough; and (C) Fully rough.The dashed line corresponds to the viscous sublayer thickness and the solid line represents the rough surface.

Figure 2 .
Figure 2. (A) Hydraulically smooth; (B) Transitionally rough; and (C) Fully rough.The dashed line corresponds to the viscous sublayer thickness and the solid line represents the rough surface.

Figure 3 .
Figure 3. Velocity shift in the case of rough walls [21] and comparison with the smooth case [65].

Figure 3 .
Figure 3. Velocity shift in the case of rough walls [21] and comparison with the smooth case [65].

Figure 4 .
Figure 4. Illustration of the trapped fluid between roughness elements.

Figure 4 .
Figure 4. Illustration of the trapped fluid between roughness elements.

Figure 5 .
Figure 5.Typical control volume for ice accretion simulations.The white arrows in Figure 5 symbolize the main contributions to the mass balance (Equation (31)) and the energy balance (Equation (32)) where .m (kg/s) and .Q (J/s) are the mass flux rate and energy flux rate, respectively.The subscripts correspond to the ones illustrated in Figure 5.

Fluids 2023, 8 ,
x FOR PEER REVIEW