Embedded LES of a turbulent thermal boundary layer over ice roughness

Contents lists available at ScienceDirect Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid Embedded LES of a turbulent thermal boundary layer over ice roughness Denis Sotomayor-Zakharov a,∗,1, Riccardo Gaudioso a,b,2, Mariachiara Gallia a,3 aInstitute of Fluid Mechanics, TU Braunschweig, Hermann-Blenk-Straße 37, Braunschweig, 38108, Niedersachsen, Germany bDepartment of Aerospace Science and Technology, Politecnico di Milano, Via La Masa 34, Milano, 20156, Italy A R T I C L E I N F O Keywords: Large eddy simulation Boundary layer Heat transfer Thermal Passive scalar Roughness Non-homogeneous Icing A B S T R A C T The numerical prediction of ice accretion via icing codes relies on the proper estimation of the heat transfer coefficient on rough ice geometries. To understand the heat transfer physics at play, direct numerical simulations (DNS) on rough surfaces can be performed, although this results in a very expensive option if geometries obtained from different icing conditions want to be analyzed. Large eddy simulation (LES) presents itself as a less expensive option to perform such studies, giving insight into the physics of turbulence, as well as opening the possibility for calibration of roughness models. The present study verifies and validates a setup to perform embedded LES (ELES) of a zero pressure gradient incompressible flow over a flat plate with ice roughness heated to a constant wall temperature. An experimental database is used, which provides the geometries of rough plates obtained from unwrapped scans of ice shapes generated on a NACA0012 airfoil. The low-speed flow over the flat plate presents a 𝑅𝑒𝐿= 3.85⋅105 and 𝑃 𝑟 = 0.729. The verification is carried out by analyzing the effects of the mesh resolution and the domain span on wall properties such as the skin friction coefficient and Stanton number. Additionally, an analysis of turbulence-related flow statistics is performed to guarantee the proper development of turbulence. The validation shows good agreement between ELES results and experimental data, especially for the Stanton number distributions, showcasing that this setup can be used for the study of heat transfer on ice roughness. 1. Introduction Icing codes are CFD tools that predict the location and shape of ice accretion on aircraft components exposed to icing conditions. Several icing codes exist such as DICEPS [1], FENSAP-ICE [2], LEWICE [3], PoliMIce [4] and IGLOO2D [5], which effectively predict ice shapes formed by the impingement of supercooled water droplets in rime ice conditions [6], with complete freezing of impinged water droplets. On the other hand, glaze ice conditions do not favor complete freezing right upon impact due to the lack of low enough freestream temperatures or due to large amounts of droplets impinging on the surface, resulting in the formation of thin water films. The calculation of the rate of ice accretion on these water films requires performing a local mass and thermal balance on the surface [7], from which accuracy issues arise due to the necessity of properly predicting the convective heat flux in the presence of rough ice surfaces. Studies such as the one of Hansman [8] note that the local rate of ice accretion on glaze ice conditions is almost proportional to the surface convective heat transfer ∗Corresponding author. E-mail address: [email protected] (D. Sotomayor-Zakharov). 1Ph.D. Research assistant, Institute of Fluid Mechanics, TU Braunschweig, Hermann-Blenk-Straße 37, 38108, Braunschweig, Germany. 2Ph.D. Student, Institute of Fluid Mechanics, TU Braunschweig, Hermann-Blenk-Straße 37, 38108, Braunschweig, Germany. Department of Aerospace Science and Technology, Politecnico di Milano, Via La Masa 34, 20156 Milan, Italy. 3Senior scientist, Institute of Fluid Mechanics, TU Braunschweig, Hermann-Blenk-Straße 37, 38108, Braunschweig, Germany. coefficient 𝐻𝑇 𝐶, meaning that icing codes need to be able to estimate this quantity with relatively high accuracy. This can be done through the use of RANS models that can take into account the ice roughness effects on the thermo-aerodynamic flow fields and surface properties such as the 𝐻𝑇 𝐶, although these ice roughness effects are yet not well understood. Research on the physics of flows over rough surfaces has led to experimental and numerical studies that look into correlating the statistical parameters of rough surfaces to a roughness-representative quantity, such as the equivalent sand-grain roughness height 𝑘𝑠, which is commonly used in roughness studies [9]. Said statistical parameters can be the root-mean-square of the surface fluctuations 𝑅𝑞, their skewness 𝑆𝑘, or an effective slope parameter, among other parameters that can be used to characterize roughness. Correlations such as the ones of Bons [10] and Flack and Schultz [11] have been elaborated and calibrated based on limited experimental databases. On the other hand, correlations such as the one of Forooghi et al. [12] are based on results of direct numerical simulations (DNS) of flows over rough https://doi.org/10.1016/j.compfluid.2025.106652 Received 20 December 2024; Received in revised form 30 March 2025; Accepted 23 April 2025 Computers and Fluids 297 (2025) 106652 Available online 5 May 2025 0045-7930/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). D. Sotomayor-Zakharov et al. surfaces. It is noted that in contrast to experiments, DNS provides a deeper insight into the physics of turbulence and heat transfer, which can be useful in the design of roughness correlations as evidenced in the studies of Thakkar et al. [13],Forooghi et al. [14],MacDonald et al. [15],Yang et al. [16,17], which have performed DNS on plane channel configurations featuring synthetic or experimental roughness. Plane channels constitute simplified domains with spanwise and streamwise periodicity, which results in the analysis of a surface with a unique homogeneous 𝑘𝑠 value. Ice roughness is inherently non-homogeneous as shown in several studies [18], i.e., 𝑘𝑠 and statistical parameters vary streamwise and the surface is exposed to growing boundary layers. This lack of streamwise periodicity can drastically increase the required computational resources for numerical studies. Cardillo et al. [19],Yuan and Piomelli [20] have performed DNS on boundary layers over homogeneous roughness, while no study has been carried out on non-homogeneous roughness surfaces. Since ice roughness presents different non-homogeneous distributions based on icing conditions and the geometry of the iced aircraft component, carrying a DNS for the study of several specific ice roughness cases could result in the usage of vast amounts of computational resources, making the elaboration and calibration of roughness correlations for icing applications very constrained. Large eddy simulation (LES) presents an option to drastically reduce the required computational resources relative to DNS, opening the possibility of studying different ice roughness cases. LES has been applied to flows over airfoils with accreted ice, being several of these works reviewed by Stebbins et al. [21], while more recent studies are found also in Wong et al. [22],Sheidani et al. [23], and Bornhoft et al. [24]. Nevertheless, these studies look towards the impact of ice accretion over aerodynamic forces and are not meant for calibration of roughness models or the calculation of 𝐻𝑇 𝐶 on ice roughness. Related to this last topic, the study of Sotomayor-Zakharov [25] presents a methodology to perform embedded LES (ELES) on scanned ice shapes obtained from icing experiments on an airfoil. In this numerical investigation, a LES region was created around the ice rough surface located at the leading edge of the airfoil in order to resolve the turbulent flow, while the rest of the flow domain was simulated via RANS. The expected enhancement of 𝐻𝑇 𝐶 due to roughness was obtained from the numerical simulations, although without any experimental 𝐻𝑇 𝐶 measurements to validate the approach. This prompted a high interest in validating this approach using an already available experimental database, which may include 𝐻𝑇 𝐶 measurements. The objective of the present study is to verify and validate an ELES setup for the simulation of boundary layers over non-homogeneous ice roughness. The experimental database of McCarrell et al. [26] is used as a reference for the numerical simulations since it presents unwrapped scans of ice shapes obtained on a NACA0012 airfoil, which result in flat plates with a non-homogeneous rough region. Additionally, the database presents measurements of skin friction 𝑐𝑓 and Stanton number 𝑆𝑡 (property analogous to 𝐻𝑇 𝐶) performed on these rough plates, which are used for validation purposes. The ELES are carried out using the software ANSYS FLUENT 21R2 and are based on the geometry and flow conditions of the experiments, resulting in the simulation of an incompressible zero pressure gradient (ZPG) flow over a heated rough plate. The temperature is treated as a passive scalar, therefore, it does not have any influence on the momentum transport, i.e., no effect on the velocity field or in fluid properties such as density or dynamic viscosity. The present manuscript is organized as follows: a detailed description of the numerical setup for the ELES is presented in Section 2, followed by a description of the tested cases and flow conditions for the study included in Section 3. Results related to the verification of the setup are presented in Section 4 which describe the effects of domain size and mesh resolution, while a comparison with the experimental measurements and validation of the setup is presented in Section 5. The conclusions of the study are presented in Section 6. Fig. 1. Experimentally obtained rough plates by McCarrell et al. [26]. Coloring displays local heights ℎ of roughness elements. Geometries: Left: 113 012.05 (low roughness). Right: 113012.04 (high roughness). Fig. 2. Schematic of the experimental setup used by McCarrell et al. [26] in the heated rough plate experiments at the BSWT. 2. Numerical setup 2.1. Experimental data The investigated surfaces were taken from the study of McClain et al. [27], which contains digital geometries of rough plates based on the unwrapped scans of ice shapes generated on a NACA0012 under icing conditions. For the experiment, these rough plates were scaled up 10 times and 3D printed in both plastic and aluminum. As an example, the 2 selected rough plates for the present numerical study are shown in Fig. 1, which present a size of 952.5 mm of streamwise length and 203.2 mm of span length. The streamwise non-homogeneous nature of ice roughness can be clearly observed, with geometry 113 012.05 presenting lower roughness levels than geometry 113012.04. In the experiment of McCarrell et al. [26], the rough plates were instrumented to perform measurements of the 𝐻𝑇 𝐶 via infrared cameras and surface heating, in addition to skin friction measurements, using a flow velocity 10 times slower to match the Reynolds number of the icing conditions after the scaling up procedure. The measurements were performed at the Baylor University Subsonic Wind Tunnel (BSWT) following the setup schematized in Fig. 2. The obtained 𝐻𝑇 𝐶 measurements are presented in the form of Stanton number 𝑆𝑡 vs. Reynolds number 𝑅𝑒𝑥 in Fig. 3 for the plastic and aluminum surfaces, being 𝑆𝑡 defined by Eq. (1), while 𝑅𝑒𝑥 is defined by Eq. (2). 𝑆𝑡 =𝐻𝑇 𝐶 𝜌 𝑐𝑃𝑈∞ (1) 𝑅𝑒𝑥=𝜌 𝑈∞𝑥 𝜇(2) Computers and Fluids 297 (2025) 106652 2 D. Sotomayor-Zakharov et al. Fig. 3. Measurements of 𝑆𝑡 vs. 𝑅𝑒𝑥 by McCarrell et al. [26] using geometries 113012.05 (low roughness) and 113012.04 (high roughness). Left: Plastic surface. Right: Aluminium surface. Fig. 4. Flow domain for the numerical simulations. Table 1 Conditions and flow properties for numerical simulations. 𝑈∞𝜌 𝜇 𝐿 𝑅𝑒𝐿 [m/s] [kg/m3] [Pa s] [mm] [–] 6.672 1.124 1.858 ⋅10−5 952.5 3.85 ⋅105 𝑇𝑏𝛥𝑇 𝑐𝑃𝜆 𝑃 𝑟 [K] [K] [J/kg K] [W/m K] [–] 300 10 1007 2.565 ⋅10−2 0.729 Here, 𝜌 is the fluid density, 𝑈∞ is the freestream velocity, 𝑐𝑃 is the fluid-specific heat capacity at constant pressure, 𝜇 is the fluid dynamic viscosity and 𝑥 is the streamwise position, with 𝑥= 0 mm being at the front of the rough plate. It is noted that the surface presents an unheated region extending up to 𝑥= 44 mm, equivalent to a 𝑅𝑒𝑥= 0.18 ⋅105, implying a delayed formation of the thermal boundary layer relative to the velocity boundary layer. The obtained data is compared with the theoretical values of 𝑆𝑡 for a laminar and turbulent boundary layer on a smooth surface with a constant heat flux and an unheated surface length 𝜉= 44 mm, being their definition presented in Eqs. (3) and (4), respectively [28], where 𝑃 𝑟 =𝑐𝑃𝜇∕𝜆 is the Prandtl number, with 𝜆 being the fluid thermal conductivity. The values of the fluid properties reproduce the air conditions of the experiment, listed in Table 1 in Section 3.1. 𝑆𝑡lam,HF = 0.453 𝑅𝑒−1∕2 𝑥𝑃 𝑟−0.66 [1 − (𝜉 𝑥)3∕4]−1∕3 (3) 𝑆𝑡turb,HF = 0.030 𝑅𝑒−1∕5 𝑥𝑃 𝑟−0.40 [1 − (𝜉 𝑥)9∕10]−1∕9 (4) Fig. 3 shows that the experimental data match the laminar correlation while presenting values above the smooth turbulent one after the roughness-induced transition takes place in the regions between 𝑅𝑒𝑥= 1⋅105 and 𝑅𝑒𝑥= 2 ⋅105. Nonetheless, in the cases with the aluminum surface, a sudden increase of 𝑆𝑡 can be observed at 𝑅𝑒𝑥= 1 ⋅105. It was discussed with the authors of the experimental work that the reason behind this was the high conductivity that aluminum presents, conveying heat horizontally towards the front regions of the rough plate. This is not an effect of an apparent boundary layer transition, but rather thermal conduction, disregarded in the present study. Therefore, several differences in the experimental measurements between both types of surfaces are attributed to this conductive effect being stronger in the aluminium surface in contrast to the plastic one. 2.2. Numerical domain The flow domain for the numerical simulations is depicted in Fig. 4, based on the reference experimental setup presented in the previous section. The domain height 𝛥𝐻 = 304.8 mm is preserved from the experiments, while the domain span 𝛥𝑍 is varied as a parameter in the present simulations, extending from 𝑧max =𝛥𝑍∕2 to 𝑧min = −𝛥𝑍∕2. The domain is divided into three subdomains: a Laminar subdomain, a LES subdomain, and a RANS subdomain, with interfaces between subdomains at locations 𝑥= 101.6 mm (𝑅𝑒𝑥= 0.41 ⋅105) and at 𝑥= 1117.6 mm (𝑅𝑒𝑥= 4.51 ⋅105). An analysis of the effect of the location of the interfaces on the simulation results is presented in Appendix. It can be seen that the 𝑥 coordinate is set to be zero at the beginning of the rough plate, which in the experiments extends up to 𝑥= 952.5 mm (𝑅𝑒𝑥= 3.85 ⋅105). This measure is taken as the reference length scale 𝐿= 952.5 mm. Nevertheless, the flat plate is extended with a smooth zone from 𝑥= 952.5 up to 𝑥= 4064.0 mm (𝑅𝑒𝑥= 16.41 ⋅105) for numerical reasons. The rough plate and its extension are no-slip walls set at a constant wall temperature 𝑇𝑤. This approach simplifies the analysis compared to applying a prescribed heat flux to the base of the rough plate, as in the experiment, which would involve modeling heat conductivity in the solid body. Although modeling solid body conduction in the simulation could capture the differences exhibited by the materials used in the experiments in the laminar region, this is out of the scope of the present study. Moreover, it is expected that considering an isothermal boundary condition instead of a constant heat flux one may have a strong influence on the results in regions where the boundary layer is laminar while presenting minimum influence on the turbulent boundary layer [28], which develops on the rough region. This is supported by comparing Eqs. (3) and (4) with their isothermal counterpart, defined in Eqs. (5) and (6). The 𝑆𝑡 in the laminar region differs by ∼36%, while in the turbulent region the difference falls down to ∼4%. The isothermal and constant heat flux trends are compared with experimental data and numerical results in Section 5. 𝑆𝑡lam,T = 0.332 𝑅𝑒−1∕2 𝑥𝑃 𝑟−0.66 [1 − (𝜉 𝑥)3∕4]−1∕3 (5) 𝑆𝑡turb,T = 0.0287 𝑅𝑒−1∕5 𝑥𝑃 𝑟−0.40 [1 − (𝜉 𝑥)9∕10]−1∕9 (6) A slip and adiabatic region before the rough plate is incorporated from 𝑥= −1016 mm up to 𝑥= 0 mm. The top wall of the entire flow domain is set as a symmetry condition (slip and adiabatic wall), as well as the side walls of the Laminar and the RANS subdomains, while the side walls of the LES subdomain are set as spanwise periodic interfaces.4 A homogeneous flow with velocity 𝑈∞, temperature 𝑇∞< 𝑇𝑤, and no turbulent content enters the flow domain via the Inlet. Results have 4Setting the side walls of the Laminar and RANS subdomains to spanwise periodic interfaces as with the LES does not have any impact on the numerical results in the LES subdomain, which covers the region of interest. Computers and Fluids 297 (2025) 106652 3 D. Sotomayor-Zakharov et al. Fig. 5. 𝛿lam and 𝛿turb vs. 𝑅𝑒𝑥. Vertical dashed lines ( ) correspond to the locations of the interfaces. shown that turbulence eventually develops due to roughness inside the LES subdomain due to a roughness-induced transition of the boundary layer. Since the flow presents a ZPG, the Outlet is set at zero relative static pressure as the Inlet. A boundary layer thickness 𝛿 is estimated following Eqs. (7) and Eqs. (8) for a laminar and turbulent boundary layer over a smooth plate [29], respectively. The thickness evolution along the plate length is presented in Fig. 5. The estimation delivers a 𝛿turb = 31 mm at the LES-RANS interface, implying a 𝛥𝐻∕𝛿turb ≈ 9.7 at this position, and indicating that the 𝛥𝐻 from the experimental setup is sufficient for the simulations to avoid boundary layer constraining. In addition, 𝛿turb needs to be considered when choosing 𝛥𝑍, and this aspect is analyzed in detail in Section 3.2. 𝛿lam = 4.91 𝑅𝑒−1∕2 𝑥𝑥(7) 𝛿turb = 0.38 𝑅𝑒−1∕5 𝑥𝑥(8) By subdividing the rough plate in streamwise sub-regions of extension 𝛥𝑥div = 25.4 mm (𝛥𝑅𝑒𝑥= 0.103 ⋅105), the r.m.s. of the spatial fluctuations of the rough geometry 𝑅𝑞 was obtained at each subregion as described in Sotomayor-Zakharov et al. [18]. Summarizing this method, 𝑅𝑞 is computed at the 𝑥 center position 𝑥𝑝 of a specific sub-region via Eq. (9), where ℎ is the location in the 𝑦 coordinate of the rough surface and ℎ𝑚 is the mean location in the 𝑦 coordinate of the rough surface in the current sub-region, computed via Eq. (10). 𝑅𝑞(𝑥𝑝) = ⎡⎢⎢⎢⎢⎣ 1 𝛥𝑥div 𝛥𝑍 𝑥𝑝+𝛥𝑥div 2 ∫ 𝑥𝑝−𝛥𝑥div 2 𝑧max ∫ 𝑧min (ℎ(𝑥, 𝑧) − ℎ𝑚(𝑥𝑝))2𝑑𝑥 𝑑𝑧⎤⎥⎥⎥⎥⎦ 1 2 (9) ℎ𝑚(𝑥𝑝) = 1 𝛥𝑥div 𝛥𝑍 𝑥𝑝+𝛥𝑥div 2 ∫ 𝑥𝑝−𝛥𝑥div 2 𝑧max ∫ 𝑧min ℎ(𝑥, 𝑧)𝑑𝑥 𝑑𝑧 (10) Fig. 6 presents the 𝑅𝑞 vs. 𝑅𝑒𝑥 distributions over the rough plate of the investigated surfaces. By using 𝑅𝑞 to quantify roughness levels, it can be observed that the roughness starts to slightly manifest at a location of 𝑅𝑒𝑥≈ 0.7⋅105, being this a distance of more than 20 𝛿lam from the Laminar-LES interface. On the other hand, roughness almost vanishes at 𝑅𝑒𝑥≈ 4 ⋅105, being this position located at slightly more than 4𝛿turb from the LES-RANS interface. Also, surface 113012.05 displays lower roughness levels in contrast to surface 113 012.04, explaining the differences that heat transfer measurements present in Fig. 3. Indeed, a more rough surface would imply a higher heat transfer enhancement in contrast to a smooth surface. 2.3. Mesh characteristics An adaptive cartesian mesh is generated in the LES subdomain, with cubic elements that progressively halve their edge size as they approach the boundary corresponding to the rough plate, resulting in different mesh level resolutions. The maximum mesh level starting at the top wall presents an edge size in 𝑥, 𝑦, and 𝑧 directions of 𝛥𝑙max = 6.4 mm Fig. 6. 𝑅𝑞 vs. 𝑅𝑒𝑥 of investigated surfaces. Vertical dashed lines ( ) correspond to the locations of the interfaces. Fig. 7. Mesh structure used for the LES subdomain. The mesh gets refined as it approaches the rough plate. (=𝛥𝑥max =𝛥𝑦max =𝛥𝑧max). At the rough wall, an edge size of 𝛥𝑙res = 0.4 mm (=𝛥𝑥res =𝛥𝑦res =𝛥𝑧res) was selected for the minimum mesh level, acting as the scale for the surface discretization. This results in 5 mesh levels, as shown in Fig. 7. Close to the rough plate, the cartesian mesh merges with a body-fitted boundary layer mesh, which presents a first layer thickness 𝛥𝑦min = 0.02 mm and a geometrical growth ratio of 1.15. The number of layers was selected by trying to match the last layer of the boundary layer mesh to half of 𝛥𝑙res. Additionally, while moving streamwise, the finer levels of the cartesian mesh tend to cover a bigger portion of the subdomain, accounting for the boundary layer growth, as shown in Fig. 8, which presents contours of time-averaged non-dimensional velocity |𝑢|∕𝑈∞ on an 𝑥-𝑦 plane cut of the LES subdomain of the case with the highest roughness levels (case 2, see Section 3.2). It can be seen that the boundary layer mainly resides in the 4th and 5th mesh levels, the finest of the LES subdomain. The Laminar and RANS subdomains are discretized with hexahedral elements, which grow geometrically in the 𝑦 direction from a first layer thickness 𝛥𝑦min = 0.02 mm at the slip region, rough plate and plate extension until a 𝛥𝑦 = 6.4 mm, keeping this size until reaching the top wall. Both domains present a spanwise discretization of 𝛥𝑧 = 6.25 mm. Also, in the Laminar subdomain, the mesh grows geometrically streamwise from 𝑥= 0 mm, from a value of 𝛥𝑥 = 0.8 mm until reaching 𝛥𝑥 = 6.4 mm at the Laminar-LES interface. On the other direction, from the interface towards the inlet, it grows from 𝛥𝑥 = 0.8 mm to values of 𝛥𝑥 ≈ 50 mm. In the RANS subdomain, the mesh grows geometrically streamwise from the LES-RANS interface starting with a value of 𝛥𝑥 = 6.4 mm towards a value of 𝛥𝑥 ≈ 150 mm at the outlet. It is noted that an analysis of the effect of the LES-RANS interface on the results is presented in Appendix. Both the Laminar and RANS subdomains present around 105 mesh elements, this being 1% of the amount of mesh elements present in the LES subdomain. The listed quantities are scaled in wall units relative to the case in analysis, being these presented in Table 2 in Section 3.2. Computers and Fluids 297 (2025) 106652 4 D. Sotomayor-Zakharov et al. Fig. 8. Contours of |𝑢|∕𝑈∞, along with mesh level distribution of the LES subdomain of case 2 (see Section 3.2). Black lines indicate the limits between mesh levels. The finer mesh levels normally grow from the rough plate while moving streamwise to account for the boundary layer growth. Table 2 Test cases for numerical simulations. ID Elements 𝛥𝑙res (Levels) 𝛥𝑙+ res 𝛥𝑦+ min 𝑅+ 𝑞,max 𝛥𝑍 (𝑆𝐹 ) 1 11.0 ⋅1060.4 mm (5) 4–20 0.1–0.5 70 50 mm (24%) 1-Ext. 22.4 ⋅1060.4 mm (5) 4–20 0.1–0.5 70 100 mm (12%) 1-Doub. 22.0 ⋅1060.4 mm (5) 4–20 0.1–0.5 70 100 mm (24%) 1-Coar. 3.3 ⋅1060.8 mm (4) 9–40 0.1–0.5 70 50 mm (24%) 1-Fine 36.6 ⋅1060.2 mm (6) 2–10 0.1–0.5 70 50 mm (24%) 2 11.8 ⋅1060.4 mm (5) 4–30 0.1–0.8 265 50 mm (24%) 1-𝜉11.0 ⋅1060.4 mm (5) 4–20 0.1–0.5 70 50 mm (24%) 2-𝜉11.8 ⋅1060.4 mm (5) 4–30 0.1–0.8 265 50 mm (24%) 2.4. Numerical schemes The ELES was performed considering an incompressible flow with constant properties, using the software ANSYS FLUENT 21R2. The intended physics to be solved are modeled via the Navier–Stokes equations presented in Eq. (11) for the momentum transport and in Eq. (12) for the thermal transport, where 𝑢𝑖 is the flow velocity in the 𝑖 direction, 𝑝 is the flow static pressure and 𝜃=𝑇𝑤−𝑇 is the flow relative temperature. This last property is treated as a passive scalar, therefore, having no influence on the momentum transport. 𝜕𝑢𝑖 𝜕𝑡 +𝜕𝑢𝑖𝑢𝑗 𝜕𝑥𝑗 = − 1 𝜌 𝜕𝑝 𝜕𝑥𝑖 +𝜇 𝜌 𝜕2𝑢𝑖 𝜕𝑥2 𝑗 (11) 𝜕𝜃 𝜕𝑡 +𝜕𝑢𝑗𝜃 𝜕𝑥𝑗 =𝜆 𝜌 𝑐𝑃 𝜕2𝜃 𝜕𝑥2 𝑗 (12) The Laminar subdomain does not present any turbulence model, and, given the absence of Inlet turbulence content, the boundary layer will remain laminar up to the Laminar-LES interface (𝑅𝑒𝑥= 0.41 ⋅105), solving Eqs. (11) and (12) in their current form. The RANS subdomain uses the SST k-𝜔 turbulence model [30] therefore modifying Eqs. (11) and (12) into RANS form. A constant turbulent Prandtl number of 𝑃 𝑟𝑡= 𝑐𝑃𝜇𝑡∕𝜆𝑡= 0.85 is used to calculate the turbulent thermal conductivity 𝜆𝑡 from the eddy viscosity 𝜇𝑡. The LES subdomain uses the WALE subgridscale (SGS) model of Nicoud and Ducros [31], therefore modifying Eqs. (11) and (12) into explicit LES form, allowing the solution of the flow field down to the surface without any extra roughness models, resulting in a wall-resolved LES. The study of Sotomayor-Zakharov [25] carries a numerical scheme analysis and validation on canonical cases, from which a subgrid-scale constant 𝐶𝑤= 0.325 is selected, as well as a SGS Prandtl number 𝑃 𝑟sgs =𝑐𝑃𝜇sgs∕𝜆sgs = 0.85 to calculate the SGS thermal conductivity 𝜆sgs from the SGS viscosity 𝜇sgs. This numerical schemes analysis establishes as well the usage of the non-iterative time advancement (NITA) solver with 2nd order implicit transient scheme [32], pressure-implicit with splitting operators (PISO) scheme for pressure-velocity coupling with neighbor correction [32], a 2nd order pressure scheme and a bounded central difference (BCD) scheme for convective terms, with a BCD boundedness parameter of Fig. 9. Variation of mean values of 𝑐𝑓 and 𝑆𝑡 vs. 𝑡′ inside the LES subdomain for case 1 (see Section 3.2). Vertical dashed lines ( ) correspond to the initial time 𝑡′ 𝑙 for sampling of flow statistics. 0.7 [33]. Lastly, the least-squares cell-based method is employed for the computation of spatial gradients. The ELES simulation is initialized from the results of a RANS simulation using the SST k-𝜔 𝛾-𝑅𝑒𝜃 transition model [34] in the whole flow domain.5 The converged RANS solution is perturbed to excite the initial development of turbulent content. This strategy is based on a synthetic turbulence generation technique provided by the solver [35], which generates a perturbed velocity field out of steady-state RANS results and it is used to reduce the time needed for the simulation to reach a statistically stable state. The ELES simulations are computed with a time-step 𝛥𝑡 that guarantees a maximum CFL <0.5, corresponding to a non-dimensional time-step 𝛥𝑡′=𝑈∞𝛥𝑡∕𝐿≈ 1.75 ⋅10−5. The computations were carried with a flow-time 𝑡𝐼′=𝑈∞𝑡𝐼∕𝐿≈ 2, which was long enough until averageable behavior could be observed on flow properties, such as on the mean values of 𝑐𝑓 and 𝑆𝑡 computed from the LES subdomain, as shown in Fig. 9 (obtained from case 1, see Section 3.2). Once this state was reached, flow statistics started to be computed during an extra amount of flow-time 𝑡′ 𝑆=𝑈∞𝑡𝑆∕𝐿≈ 2.5. 2.5. Analysis of wall properties The wall properties are computed from time-averaged flow properties at the rough surface. These quantities are obtained through the operation presented in Eq. (13), where 𝜓 is the time-average during a period 𝑡𝑆 of an arbitrary property 𝜓. 𝜓(𝑥, 𝑦, 𝑧) = 1 𝑡𝑠 𝑡𝐼+𝑡𝑆 ∫ 𝑡𝐼 𝜓(𝑥, 𝑦, 𝑧, 𝑡)𝑑𝑡 (13) Therefore, the time-averaged pressure 𝑝, shear stress 𝜏, and heat flux 𝑞 are extracted on the rough surface and then applied in Eq. (14) to 5The transition model is exploited to initialize the ELES with an already solved laminar boundary layer up to the position where roughness starts to appear. Computers and Fluids 297 (2025) 106652 5 D. Sotomayor-Zakharov et al. Fig. 10. 𝑅𝑞 vs. 𝑅𝑒𝑥 of cases based on 𝛥𝑍. Vertical dashed lines ( ) correspond to the locations of the interfaces. compute 𝜏𝑤 and in Eq. (15) to compute 𝑞𝑤. Here, 𝐴𝑤 corresponds to the wetted area of the rough surface in the current sub-region, 𝜏𝑥 is the 𝑥 direction component of 𝜏, and 𝜂𝑥 is the 𝑥 direction component of the normal of 𝑑𝐴𝑤. 𝜏𝑤(𝑥𝑝) = 1 𝛥𝑥div 𝛥𝑍 ∫ 𝐴𝑤(𝜏𝑥−𝑝 𝜂𝑥)𝑑𝐴𝑤(14) 𝑞𝑤(𝑥𝑝) = 1 𝛥𝑥div 𝛥𝑍 ∫ 𝐴𝑤 𝑞 𝑑𝐴𝑤(15) Finally, 𝑐𝑓 and 𝑆𝑡 are calculated for each sub-region via Eqs. (16) and (17), respectively, where 𝛥𝑇 =𝑇𝑤−𝑇∞. Following Eq. (1), the definition of 𝐻𝑇 𝐶 =𝑞𝑤∕𝛥𝑇 can be retrieved. 𝑐𝑓=2𝜏𝑤 𝜌∞𝑈2 ∞ (16) 𝑆𝑡 =𝑞𝑤 𝜌 𝑐𝑃𝑈∞𝛥𝑇 (17) 2.6. Analysis of flow statistics The analysis of the flow statistics is performed through the extraction of boundary layer profiles in the 𝑦 direction at locations 𝑥𝑝 with the intention of matching the profile with extracted wall properties for scaling at each sub-region. Therefore, the operation ⟨𝜓⟩ is presented for an arbitrary time-averaged property 𝜓, which stands for plane-averaging in the 𝑥 and 𝑧 directions. ⟨𝜓(𝑥𝑝, 𝑦)⟩=1 𝛥𝑥div 𝛥𝑍 𝑥𝑝+𝛥𝑥div 2 ∫ 𝑥𝑝−𝛥𝑥div 2 𝑧max ∫ 𝑧min 𝜓(𝑥, 𝑦, 𝑧)𝑑𝑥 𝑑𝑧 (18) Properties on the boundary layer profiles such as the velocity ⟨𝑢⟩+, temperature ⟨𝜃⟩+, and Reynolds and thermal stresses are scaled in wall units, indicated by the superscript (+), by means of the friction velocity 𝑢𝜏 and temperature 𝑇𝜏, which are computed at each sub-region via Eqs. (19) and (20), respectively. 𝑢𝜏=√𝜏𝑤 𝜌(19) 𝑇𝜏=𝑞𝑤 𝜌𝑐𝑃𝑢𝜏 (20) The wall-distance 𝑦+ is scaled via Eq. (21), noting that ℎ𝑚 is taken into account to represent that the profile begins at the mean roughness height of the sub-region. 𝑦+=𝜌 𝑢𝜏(𝑦−ℎ𝑚) 𝜇(21) 3. Test cases 3.1. Flow properties Based on the data of McCarrell et al. [26], temperatures of 𝑇∞= 295 K at the freestream and 𝑇𝑤= 305 K at the wall were selected for the current study. This leads to a bulk temperature 𝑇𝑏= 0.5 (𝑇𝑤+𝑇∞) = 300 K, from which the air properties are estimated and presented in Table 1. Such selection of parameters implies a temperature difference of 𝛥𝑇 = 10 K at the wall. Properties as 𝜌, 𝑐𝑃, 𝜇 and 𝜆 are kept constant after verifying that their change is around ±1.7%, ±0.0%, ±1.3% and ±1.5% for a ±5 K variation of 𝑇𝑏, respectively. Additionally, a Prandtl number 𝑃 𝑟 = 0.729 is obtained, as well as a Reynolds number based on the reference length of 𝑅𝑒𝐿= 3.85 ⋅105 and a Mach number of 𝑀𝑎∞= 0.02 justifying the incompressible flow assumption. 3.2. Selected cases The test cases considered for the present investigation are summarized in Table 2, where their designation (ID) is displayed. Case 1 is based on the experimental geometry 113 012.05, case 2 is based on 113012.04, and these are considered the reference cases. The selection of a small 𝛥𝑍 is encouraged to spare computational costs. Initially, the reference cases 1 and 2 are selected with a 𝛥𝑍 = 50 mm, covering a portion of the original geometry from 𝑧max = 25 mm to 𝑧min = −25 mm. However, two important aspects must be taken into consideration and are subjected to analysis. The first aspect relates to the effect of 𝛥𝑍 over 𝑅𝑞 distributions. Since the original experimental geometries present a total span of 𝛥𝑍 = 203.2 mm (see Section 2.1), the selection of a smaller 𝛥𝑍 has the potential of affecting flow properties if certain surface features are removed by this reduction or smoothed by any other additional procedure. This removal or smoothing can be observed by comparing the 𝑅𝑞 distributions of the total span geometry with the one with a smaller 𝛥𝑍. To analyze the effect that the selection of 𝛥𝑍 has on the flow properties, case 1-Ext. (Extended) with a 𝛥𝑍 = 100 mm is established (𝑧max = 50 mm to 𝑧min = −50 mm), effectively presenting twice the span than the reference case. Fig. 10 is presented, showing the effect of the choice of 𝛥𝑍 on the distributions of 𝑅𝑞 across 𝑅𝑒𝑥. It can be seen that the reference cases present lower values of 𝑅𝑞, while case 1-Ext. presents values closer to the trend of the original geometry. The main reason behind the differences between 𝑅𝑞 of the cases and the original geometry is an artificial procedure performed on the flow domain construction to ensure that the side walls of the LES subdomain would match for periodicity: the rough geometry close to the side walls is smoothed towards values of the local mean geometrical height in the 𝑦 direction using a sinusoidal function as a filter in the 𝑧 direction. The size of the filter at each side is 𝛥𝑍filter = 6 mm, which showed proper smoothing and good geometrical behavior. This implies that a fraction of the span of 𝑆𝐹 = 2 𝛥𝑍filter∕𝛥𝑍 is affected by the filter, being 𝑆𝐹 = 24% for 𝛥𝑍 = 50, resulting in a smoothing effect close to the side walls, and ultimately reducing 𝑅𝑞. The effect of the filter is mitigated as a larger 𝛥𝑍 = 100 mm is chosen, resulting in a 𝑆𝐹 = 12%. This leads to lesser 𝑅𝑞 decay, which is analyzed in case 1-Ext. On the other hand, further reduction of the span below 𝛥𝑍 = 50 would imply an increase on 𝑆𝐹 , which, for example, could reach values of 48% for a 𝛥𝑍 = 25 mm, having almost half of the domain filtered. It is noted that a second reason for the 𝑅𝑞 difference could be due to a strong reduction of 𝛥𝑍 below representative spatial wavelengths that form the roughness on the surface. The previously selected values of 𝛥𝑍 are not small enough for this effect to be prominent in comparison to the smoothing effect, although both affect the final 𝑅𝑞 distribution, potentially presenting an effect on flow properties such as the 𝐻𝑇 𝐶. The second aspect relates to the effect of 𝛥𝑍 over turbulence since small constraining span values may not allow the proper development of turbulent content on the boundary layer. The reference cases present Computers and Fluids 297 (2025) 106652 6 D. Sotomayor-Zakharov et al. Fig. 11. Cases with different mesh resolution 𝛥𝑙res. (a) Case 1-Coar. (b) Case 1 (reference). (c) Case 1-Fine. 𝛥𝑍 = 50 mm, which results in a 𝛥𝑍∕𝛿turb ≈ 1.6 at the LES-RANS interface, and a 𝛥𝑍∕𝛿turb ≈ 2.5 at the middle of the LES subdomain, where roughness is prominent. This 𝛥𝑍 is expected to be sufficient to capture the development of spanwise turbulence structures according to previous numerical and experimental studies. For example, Lund et al. [36] used a 𝛥𝑍∕𝛿turb ≈𝜋∕2 at the center of a LES ZPG turbulent boundary layer, while experiments from Tomkins and Adrian [37] found evidence of large structures in the logarithmic layer extending up to 0.6𝛿turb spanwise. Case 1-Ext. with a 𝛥𝑍 = 100 mm (𝛥𝑍∕𝛿turb ≈ 3.2) could be used to prove this aspect by comparing it to Case 1. However, the previously stated change on 𝑅𝑞 introduced by the spanwise extension could mislead the comparison of the resulting flow properties and statistics between case 1-Ext. and case 1. Therefore, to isolate the effect of 𝛥𝑍 on turbulent flow statistics, case 1-Doub. (Double) is introduced, resulting from the pairing of two identical case 1 geometries next to each other spanwise. This implies a 𝛥𝑍 = 100, which preserves the same 𝑅𝑞 vs. 𝑅𝑒𝑥 distribution as case 1, allowing the comparison of turbulent properties under domains with a different span but with the same roughness effects. In addition to the span analysis, a mesh study is carried out with a coarser and refined version of case 1, designated as case 1-Coar. and case 1-Fine, respectively. The coarsening/refinement is controlled by modifying the 𝛥𝑙res to twice/half its original value, resulting in the meshes shown in Fig. 11, which show the different cartesian mesh refinement levels. It was analyzed that at these resolutions, the choice of 𝛥𝑙res does not have any notable effect over 𝑅𝑞 distributions due to the effect of surface meshing. Finally, all previously stated cases consider that the entire rough plate and its extension are at 𝑇𝑤. Nevertheless, cases 1-𝜉 and 2-𝜉 are introduced with an unheated length of 𝜉= 44 mm as in the experiments to assess its effect on the results, treating the region between 𝑥= 0 mm and 𝑥= 44 mm as adiabatic. 4. Verification 4.1. Effect of mesh resolution 𝛥𝑙𝑟𝑒𝑠 The influence of 𝛥𝑙𝑟𝑒𝑠 on the results is analyzed by comparing case 1 with case 1-Coar. and case 1-Fine. Fig. 12 presents computed values of 𝑐𝑓 and 𝑆𝑡 vs. 𝑅𝑒𝑥 across cases. As a reference, the 𝑐𝑓 vs. 𝑅𝑒𝑥 distribution for a laminar (Eq. (22)) and turbulent (Eq. (23)) boundary layer over a smooth plate is presented, as well as the 𝑆𝑡 vs. 𝑅𝑒𝑥 distributions for a laminar (Eq. (5)) and turbulent (Eq. (6)) boundary layer over a smooth plate at a constant temperature. It can be seen that as 𝛥𝑙𝑟𝑒𝑠 is refined, the distributions of 𝑐𝑓 and 𝑆𝑡 tend to converge, being the greatest differences between mesh levels located at 𝑅𝑒𝑥≈ 1.8⋅105, corresponding to the location of the highest values of 𝑅𝑞 (see Fig. 6). All curves show that a laminar boundary layer is initially formed, which matches the theoretical smooth wall values and manages to cross the Laminar-LES interface (𝑅𝑒𝑥= 0.41 ⋅105). Transition to turbulence is induced by roughness inside the LES subdomain after 𝑅𝑒𝑥= 1 ⋅105, resulting in 𝑐𝑓 and 𝑆𝑡 distributions with values above theoretical smooth ones for turbulent flow. As 𝑅𝑞 decreases approaching the LES-RANS interface, Fig. 12. Influence of the mesh resolution 𝛥𝑙𝑟𝑒𝑠 on 𝑐𝑓 and 𝑆𝑡. Vertical dashed lines ( ) correspond to the locations of the interfaces. Table 3 Variation of 𝑐𝑓 and 𝑆𝑡 with mesh refinement. ID Mean(𝑐𝑓) Mean(𝑆𝑡) 1-Coar. 0.00717 0.00320 1 0.00737 0.00328 Var.(%) 2.76 2.45 1-Fine 0.00750 0.00334 Var.(%) 1.68 1.71 the turbulent flow recovers its smooth status, and wall properties match again the expected smooth distribution. 𝑐𝑓,lam = 0.664 𝑅𝑒−1∕2 𝑥(22) 𝑐𝑓,turb = 0.0574 𝑅𝑒−1∕5 𝑥(23) A quantification of the variation of the wall properties between a refined level with its respective coarser level is presented in Table 3, showing the mean values of 𝑐𝑓 and 𝑆𝑡 computed inside the LES subdomain. The variations of the mean values decrease as the mesh resolution is refined, demonstrating that the resolution of case 1 is appropriate for the study of the overall roughness effects on the flow. A qualitative analysis of the turbulent flow behavior is carried on between locations 𝑅𝑒𝑥= 1.8⋅105 and 𝑅𝑒𝑥= 2.0⋅105. Fig. 13 displays contours of non-dimensional instantaneous velocity |𝑢|∕𝑈∞ and time-averaged velocity |𝑢|∕𝑈∞ at the mid-plane of the flow domain (𝑧= 0 mm). It can be observed through |𝑢|∕𝑈∞ that as the mesh resolution is increased, more detailed fluctuating behavior is captured. Separation regions are observed after each roughness peak, and the highly packed roughness elements seem to not allow the flow to recover after separation, a behavior reported for d-type roughness [38–40]. From |𝑢|∕𝑈∞, it can be seen that the left separation bubble is affected Computers and Fluids 297 (2025) 106652 7 D. Sotomayor-Zakharov et al. Fig. 13. Contours of non-dimensional flow velocity for different mesh resolutions in the region of 𝑅𝑒𝑥= 1.8⋅105 to 𝑅𝑒𝑥= 2.0⋅105 and 𝑧= 0 mm. Left column: Instantaneous values |𝑢|∕𝑈∞. Right column: Time-averaged values |𝑢|∕𝑈∞. Fig. 14. Contours of non-dimensional flow temperature for different mesh resolutions in the region of 𝑅𝑒𝑥= 1.8⋅105 to 𝑅𝑒𝑥= 2.0⋅105 and 𝑧= 0 mm. Left column: Instantaneous values 𝜃∕𝜃∞. Right column: Time-averaged values 𝜃∕𝜃∞. Table 4 Locations for the extraction of boundary layer profiles. ID 𝑥 [mm] 𝑅𝑒𝑥⋅10−5 [–] 1 444.5 1.8 2 622.3 2.5 3 800.1 3.2 4 977.9 4.0 5 1079.5 4.4 by the mesh resolution between case 1-Coar. and case 1, justifying the difference of 𝑐𝑓 at 𝑅𝑒𝑥= 1.8⋅105 between these cases. In contrast, case 1 and case 1-Fine present similar flow characteristics, which results in the lower differences of 𝑐𝑓 between these cases observed in Fig. 12 and in Table 3. To analyze the temperature fields between 𝑅𝑒𝑥= 1.8⋅105 and 𝑅𝑒𝑥= 2.0⋅105, Fig. 14 is presented, displaying contours of non-dimensional instantaneous relative temperature 𝜃∕𝜃∞ and time-averaged relative temperature 𝜃∕𝜃∞ (𝜃∞=𝑇𝑤−𝑇∞) at the mid-plane of the flow domain (𝑧= 0 mm). As observed for the velocity fields, an increase in mesh resolution manages to capture more detailed behavior of spatial fluctuations of 𝜃∕𝜃∞. Still, contours of 𝜃∕𝜃∞ display a consistent behavior at different mesh resolutions, indicating small temperature gradients inside separation regions after roughness elements, while the highest temperature gradients are observed at the roughness peaks. By performing an analysis of the flow statistics, the boundary layer profiles in the 𝑦 direction across the flow domain are extracted at the 𝑥 locations listed in Table 4. Fig. 15. Mean turbulent profiles at different mesh resolutions. Dashed lines ( ) indicate the theoretical smooth plate velocity and temperature log-layers. Profiles of ⟨𝑢⟩+ and ⟨𝜃⟩+ are displayed in Fig. 15 for the three considered mesh resolutions, showing no major effects on almost all the profiles except on the 1st one, which presents higher values of ⟨𝑢⟩+ for case 1-Coar. This can be related to the lower 𝑐𝑓 prediction at 𝑅𝑒𝑥= 1.8⋅105 associated with such case. Still, the velocity and thermal shifts with respect to the log-law are clearly visible on the 2nd and 3rd profiles, as evidence of the roughness influence over the boundary layer, being this shift similar across mesh resolutions. The mean Reynolds and thermal stresses are analyzed in Fig. 16. Here, only the 2nd and 5th profiles are displayed, located in the rough zone and smooth zone, respectively, allowing the analysis of the impact of roughness on the flow statistics. The normal stresses ⟨𝑢′2⟩+, ⟨𝑣′2⟩+, ⟨𝑤′2⟩+ and ⟨𝜃′2⟩+ display a similar behavior at different resolutions for both profiles. Nevertheless, the 5th profile exhibits the expected anisotropic behavior close to the wall on smooth surfaces, where ⟨𝑢′2⟩+>⟨𝑤′2⟩+>⟨𝑣′2⟩+, reported by several studies as in Jiménez [39]. The 2nd profile exhibits the tendency to isotropy that roughness forces into the turbulent field, with ⟨𝑢′2⟩+≈⟨𝑤′2⟩+≈⟨𝑣′2⟩+. Specifically, the damping of the peak of ⟨𝑢′2⟩ is prominent due to the constraining nature of the roughness over the spanwise fluctuations, which is associated with the destruction of the buffer layer caused by the intrusion of large roughness elements, as also stated by Yuan and Piomelli [41]. Similar trends have been observed in the literature on sand-grain roughness by Cardillo et al. [19] and Yuan and Piomelli [20]. In the case of ⟨𝜃′2⟩+, a similar damping of the peak value due to roughness is observed, although its value increases around 𝑦∕𝛿𝜃= 0.5. The shear component ⟨𝑢′𝑣′⟩+ reveals sensitivity to the mesh on the 5th profile, possibly indicating a slightly different tendency to recover after the rough zone dictated by the different mesh resolutions. In particular, ⟨𝑢′𝑣′⟩+ is uniformly reduced across the boundary layer at the 2nd profile in contrast to the 5th one, manifesting once again the influence of increased local friction due to roughness. On the other hand, ⟨𝑣′𝜃′⟩+ shows similar peak values between the 2nd and 5th profile, although with a shifted peak location away from the wall due to roughness effects. Overall, the mesh resolution of case 1 manages to discretize the relevant spatial wavelengths of the rough surface and capture its interactions with the flow field, as well as not presenting strong variations on flow properties such as 𝑐𝑓 and 𝑆𝑡 distributions and turbulence statistics through further mesh refinement. This implies that this mesh level can be used for further analysis in the current study. 4.2. Effect of domain span 𝛥𝑍 The effect of 𝛥𝑍 is analyzed by comparing the results of case 1 with case 1-Ext. and case 1-Doub. It is worth reminding that 𝛥𝑍 can have mainly two effects on the results, namely the constriction of turbulent Computers and Fluids 297 (2025) 106652 8 D. Sotomayor-Zakharov et al. Fig. 16. Mean Reynolds and thermal stresses at different mesh resolutions. Only the 2nd ( ) and 5th profile ( ) are displayed. The wall distance is scaled with respect to the local boundary layer thickness, 𝛿99 or 𝛿, and thermal boundary layer thickness, 𝛿𝜃. Fig. 17. Influence of the domain span 𝛥𝑍 on 𝑐𝑓 and 𝑆𝑡. Vertical dashed lines ( ) correspond to the locations of the interfaces. content due to a small 𝛥𝑍 in relation to 𝛿, and variation on the values of 𝑐𝑓 or 𝑆𝑡 due to different 𝑅𝑞 distributions. Fig. 17 presents the results for 𝑐𝑓 and 𝑆𝑡 vs. 𝑅𝑒𝑥 distributions of case 1, case 1-Doub. and case 1Ext., together with theoretical values described by Eqs. (22), (23), (5) and (6). Case 1 and case 1-Doub. do not exhibit any notable differences, implying that 𝛥𝑍 = 50 mm allows the correct development of turbulent flow fields, although a detailed analysis of flow statistics is necessary to fully support this statement. On the other hand, case 1 and case 1Ext. display differences due to their different 𝑅𝑞 distributions, more prominent for 𝑐𝑓 in contrast to 𝑆𝑡, although the trends and orders of magnitude remain similar. Maximum peak values can be observed around 𝑅𝑒𝑥= 1.8⋅105 for all cases, with a recovery towards a smooth behavior around 𝑅𝑒𝑥= 4 ⋅105, where 𝛥𝑍∕𝛿≈ 2 and still inside the LES subdomain. Fig. 18. Mean turbulent profiles for different domain spans. Dashed lines ( ) indicate the theoretical smooth plate velocity and temperature log-layers. For a deeper analysis, the flow statistics of case 1, case 1-Doub. and case 1-Ext. are considered. Profiles of ⟨𝑢⟩+ and ⟨𝜃⟩+ are presented in Fig. 18. Case 1 and case 1-Doub. present similar profiles. The only sensible difference between these cases and case 1-Ext. is presented by the 1st profile of ⟨𝑢⟩+, corresponding to the zone before the first roughness peak. This is indeed reflected in the difference of the 𝑐𝑓 prediction at 𝑅𝑒𝑥= 1.8⋅105 showed in Fig. 17, where the higher friction captured by increasing the span is resulting in a lower ⟨𝑢⟩+ estimation. The better match in 𝑆𝑡 between the cases is found through the good agreement of the thermal profiles. Reynolds and thermal stresses are depicted in Fig. 19 for the 2nd and 5th profiles. In general, case 1 and case 1-Doub. display similar values across the stresses, supporting the claim that 𝛥𝑍 = 50 mm is large enough to sustain unconstrained turbulence development. The differences between case 1 and case 1-Ext. are more significant for the Reynolds stresses compared to the thermal ones. The largest discrepancies are observed in the quantities ⟨𝑢′2⟩+ and ⟨𝑢′𝑣′⟩+, where the predictions for case 1 and case 1-Doub. show notable variations. Thermal stresses match fairly well, although presenting the same behavior described in the mesh analysis, where roughness induces ⟨𝜃′2⟩+ to present reduced peak values next to the wall but larger values around 𝑦∕𝛿𝜃= 0.5, while ⟨𝑣′𝜃′⟩+ is shifted away from the wall. Once again, these insights justify the use of the numerical setup used for case 1 for further analysis and for validation. 5. Validation 5.1. Comparison with experimental data of 𝑐𝑓 Fig. 20 presents a comparison between the 𝑐𝑓 vs. 𝑅𝑒𝑥 distributions obtained from case 1 and case 1-Ext. with experimentally measured 𝑐𝑓 values,6 as well as the theoretical curves described by Eqs. (22) and (23). Qualitative resemblance can be observed between the ELES results, where the two peaks before and after 𝑅𝑒𝑥= 2⋅105 are captured by the simulations but overestimating their magnitudes. Nevertheless, case 1-Ext. seems to approach experimental results in a better way, being justified by the better description of 𝑅𝑞 distributions. Therefore, it can be stated that a main source of inaccuracy of numerical results with respect to the experimental data involves the usage of a 𝛥𝑍 that does not cover the entire experimental span of 𝛥𝑍 = 203.2 mm, although this result would remain valid respect to its own topography and 𝑅𝑞 distribution. 6In the results figure, the experimental data for 𝑐𝑓 was shifted in the 𝑥 direction an amount of −50.8 mm (2 in), based on discussions with the authors of the experimental study, who indicated that such a shift could exist even though it was not recorded in the paper. Indeed, this adjustment enhanced the qualitative alignment with the numerical results for both analyzed cases. Computers and Fluids 297 (2025) 106652 9