## Abstract

In the last decade, progresses in additive manufacturing (AM) have paved the way for optimized heat exchangers whose disruptive design will heavily rely on predictive numerical simulations. However, due to typical roughness induced by AM, current wall models used in steady and unsteady 3D Navier–Stokes simulations do not take into account such characteristics. For the development and assessment of novel wall models for AM, a high-fidelity roughness-resolved large-eddy simulation (RRLES) database is built. This article describes the numerical setup and the methodology used for conducting RRLES, from surface generation to postprocessing. In addition, three different cases representing two printing directions plus a streamwise and spanwise isotropic case are investigated. While the roughness distributions are the same in the three cases, the effective slope (ES) is very different, and the impact of this parameter on turbulence and heat transfer is analyzed at different Reynolds numbers.

## 1 Introduction

In aeronautical engines, heat exchangers have a fundamental role especially in controlling the temperature of the liquid fuel and oil streams. The mass of heat exchangers is also critical to reduce the total weight of the engine. The design of compact and efficient heat exchangers (CHX) is therefore of paramount importance for the future aeronautical engines. Additive manufacturing (AM) offers the opportunity to explore new designs and boost the optimization of CHX [1]. For instance, as shown by Saltzman et al. [2], they obtained an improvement around 10% compared to heat exchangers built with traditional manufacturing.

Nonetheless, surface roughness generated with AM can be substantially larger than with traditional manufacturing techniques. Surface state of metal AM test samples analyzed in several studies have indeed underlined this key point [2–7]. Moreover, the wall roughness is not isotropic in space and varies according to the main direction of printing [4,5]. Some patterns called welding tracks can appear on the surface [3]. This typical roughness has a significant impact on the performances in terms of pressure drop and heat transfer capacity [2,6,7]. In addition, an important gap between performances obtained from simulations and/or empirical correlations and those measured a posteriori exists.

Many authors have been devoted to the study of the turbulent flow over different roughness types and surfaces experimentally [8] as well as numerically [9–11]. The presence of surface roughness causes an overall thickening of the boundary layer in addition to the modification of the streamwise structures. Secondary flow motions can also be observed due to spanwise inhomogeneities on the surface, and this induces a better flow mixing. In addition, the transition between laminar to turbulent regimes appears earlier than over smooth surfaces. Concerning the mean velocity profile, a downward shift compared to the log-law profile is observed [12,13]. The latter is commonly called the roughness function (Δ*U*^{+}). Roughness has also an impact on the Nusselt number since the increase of surface exchange enhances heat transfer [14].

*k*

_{s}[15,16,7]. This parameter

*k*

_{s}can be expressed as a function of the statistical roughness parameters such as the arithmetic average height (

*S*

_{a}), the root-mean-square height (

*S*

_{q}), the skewness (Sk), and the kurtosis (Ku). The two latter correspond, respectively, to the third and the fourth moments of the height distribution of the surface. Studies have underlined the predominance of the effective slope (ES), defined as the average slope of the surface height along the streamwise direction, in predicting roughness effects [10].

*h*(

*x*,

*y*) is the surface height, the mean plane is supposed to be at

*h*= 0, and

*L*

_{x}and

*L*

_{y}are, respectively, the streamwise and spanwise lengths.

A well-known threshold equal to *ES* ≈ 0.35 delimits two regimes: (i) the sparse or waviness regime below this limit and (ii) the roughness or dense one [17]. Nevertheless, such universal correlation that would cover the two regimes has not been found up to now. Furthermore, few articles propose numerical simulation of AM surfaces.

In the aim of building new rough-wall models dedicated to AM surfaces, the creation of a roughness-resolved large-eddy simulation (RRLES) database of representative channel flows emerges as a worthwhile process. The present article describes a methodology to perform parametric RRLES from a limited set of surface parameters to control the surface height distribution and the streamwise and spanwise spatial autocorrelation functions. Section 1 describes the numerical setup, configurations, and the used methodology for performing RRLES. The postprocessing and obtained results are presented and discussed in Sec. 2.

## 2 Numerical Setup

In this section, chosen configurations and methodology are addressed. In order to ensure periodic channels, recycling conditions are applied on velocity and temperature. This process is particularly exposed here.

### 2.1 Configurations and Meshes.

First of all, rough surface (RSG) and body-fitted unstructured mesh (RRMG) generators have been developed [18]. The former creates triangulated surfaces stored in the STL format that are used by the latter to compute level set functions. The level set functions are then employed to refine an existing 3D unstructured mesh and to create a conforming body-fitted mesh with controlled cell size, quality, and size gradation. The resulting meshes are suitable for performing RRLES with higher-order finite-volume schemes. The RRMG has been developed within the YALES2 code. YALES2 is an unstructured low-Mach number code developed at CORIA for the direct numerical simulations and large-eddy simulations (LES) in complex geometries [19]. It heavily relies on the parallel volume and surface mesh adaptation.

The RSG/RRMG are utilized to generate rough planar channels but may also be applied to more complex shapes such as cylindrical tubes or plate/fin configurations. RRLES are conducted for periodic channels of size 8*H* x 3*H* x 2*H* in the streamwise (*x*), spanwise (*y*), and crosswise (*h*) directions with *H* being the half height of the channel, which is equal to 1.0 mm. This domain size was confirmed to be sufficient to have a negligible impact of the periodic recycling on the wall friction [20]. In addition, the half height was selected in order to be close to hydraulic diameters and channel heights encountered for some AM experiments as the L-2x-In case in the study by Stimpson et al. [14].

Concerning the chosen configurations, three different cases representing two printing directions plus a streamwise and spanwise isotropic case have been selected for the analysis (Fig. 1). A smooth channel is also considered as a reference. Roughness parameters that are targeted in this article are the root-mean-square height (*S*_{q}), the height distribution skewness (Sk) and kurtosis (Ku), and the ES. While the surface height distributions are the same in the three cases, the effective slope is different due to the different space auto-correlation functions.

In addition of printing direction, the idea is to generate surfaces that mimic roughness height distribution encountered in AM. An overall range of roughness parameters values can be found in the literature even though there are differences between upward- and downward-facing surfaces [2–7]. For this purpose, values of the latter for *S*_{q}, Sk, and Ku are set, respectively, to 20 *μ*m, around 0.2, and around 4.0. Final values are exposed in Table 1, and probability density functions of height are plotted in Fig. 2 for the three geometries.

With the RRMG, conforming body-fitted meshes are obtained from these geometries. A balance on the element count has to be reached in order to resolve sufficiently the flow while minimizing the cost. In addition, the meshing process has to keep the topology intact. The mesh generation process starts from a Cartesian grid, which is tessellated and adapted numerous times.

The initial numbers of elements in each direction (*N*_{x}, *N*_{y}, *N*_{h}) for the Cartesian grid are fixed in order to initially get isotropic cells of 20 *μ*m. The cell size gradient is equal to 0.1, and the cell size on the rough surface is set to Δ*y*_{w} = 7 *μ*m. The final number of mesh elements for RRLES is about 130 million cells.

All characteristics and performance of meshing process are summarized in Tables 2 and 3. The mesh for the smooth case is a stretched cartesian grid with the dimensionless wall grid spacing of $\Delta x+=13.7$, $\Delta y+=5.5$, $\Delta h+=[0.5,7.3]$ in accordance with operating conditions.

For each configuration, the mesh size is fixed for all the tested Reynolds numbers. The mesh resolution for the present study is discussed in the Appendix. It must be noted that the dimensionless roughness *S*_{a}/Δ*y*_{w} is only larger than two, which may seem insufficient for the wall resolution. However, with body-fitted grids, this latter has to be sufficient for capturing the wall roughness (curvature and position as illustrated in Fig. 3) and the wall-normal velocity gradient.

### 2.2 Statistics and Performance.

For collecting statistics, each RRLES is split into two steps. The initialization step is performed during a given transient time, and then statistics are accumulated. For each step, the durations are summarized in Table 4. For the second step, the dimensionless time *t*^{+} defined by *t*^{+} = *t u*_{τ}/*H* is around 20 in average. Interestingly, due to inhomogeneities of the surface, the time to establish a fully developed turbulent channel flow is reduced compared to smooth cases.

RRLES CPU costs are presented in Table 5. For clarity, the CPU hours are averaged over all Reynolds numbers cases for each configuration.

### 2.3 Numerics and Models

#### 2.3.1 Working Hypotheses.

In each case, incompressible flow simulations are performed. The chosen target bulk Reynolds number range is the fully developed turbulent flow. This is why RRLES are performed at Re = 5000, Re = 8000, Re = 17,000, and Re = 25,000. The fluid kinematic viscosity is set to $\nu =1.51710\u22125m2/s,$ and the maximal Courant–Friedrichs–Lewy (CFL) number used is equal to 0.8. The wall-adapting local Eddy-viscosity (WALE) subgrid scale model is retained as it is widely used for LES of boundary layer flows [21]. A fourth-order central finite-volume scheme is used, and the four-step fourth-order scheme TFV4A is applied for velocity and scalar transport prediction [22].

#### 2.3.2 Methodology, Boundary Conditions.

The automatic generation of fully periodic channels is very challenging, and ensuring periodic meshes is complex for unstructured grid. The periodicity is though necessary to reach statistical convergence of the flow. Instead of imposing strict periodic boundary conditions in the streamwise and spanwise directions, a Lagrangian recycling method has been developed. This method does not require to add a body force as source term in Navier–Stokes momentum equation to compensate the wall friction. This approach can also be used as a precursor in spatially developing boundary layers by increasing *L*_{out}. The idea behind this recycling method is simply to use a time-varying inlet boundary condition with a velocity coming from a distant plane in the channel. This three-step recycling process is based on passive Lagrangian particles as illustrated in Fig. 4. A flowrate is imposed at the inlet and for optimizing performances, and the recycling is done every *N* time steps, *N* > 1. A linear interpolation in time is performed between two recycled planes. No rescaling is applied on velocity profile as such profile is unknown a priori. Thus, contrary to what is proposed in the study by Xiao et al. [23], for instance, no target or corrections on velocity field are applied as input of this method. In this article, we select *N* = 20 as it gives the better compromise between performances and velocity interpolation errors.

The complete process is described here:

Lagrangian particles are created at the inlet of the grid (spheres at the inlet on the left in Fig. 4) and are translated to the recycling plane.

The target field at the recycling plane, for instance, velocity, is interpolated for translated particles (spheres at the recycling plane in Fig. 4)

This is the relocation step: particles are transferred at the inlet, and the field at this location is updated.

The computational domain is then defined with a recycling zone and a buffer area to avoid any influence of the outlet boundary treatment onto the recycled velocity (Fig. 5). Thus, the location of the recycling plane is primordial and should be set at a given distance from the outlet within the domain. For two parallel planes, this distance was found to be at least equal to the length between both planes, and in our cases, the latter corresponds to 2*H*. Indeed, perturbations of the velocity field due to the outlet boundary condition may be recycled and injected at the inlet if *L*_{out} is below this threshold.

In order to decrease interpolation margin errors between the inlet and the recycling plane, both locations should be the equivalent at the wall. Both the generated rough surface and conformal mesh respect this periodicity condition only in terms of wall topology. A no-slip wall boundary condition is applied on rough planes, and the other sides in the spanwise direction are considered as slip walls. This methodology is also applied in the smooth configuration. Validation of the numerical methods may be found in the Appendix.

Finally, from STL generation to postprocessing, all calculation steps are integrated into a workflow tool. This allows to manage automatically series of runs on a distant super-computer.

## 3 Simulation Analysis

This section details the different tools used for the subsequent analysis of the RRLES.

### 3.1 Nondimensional Velocity and Temperature.

*h*

_{w}being the minimal height of the surface. The variable

*φ*corresponds to the

*x*−

*y*plane porosity, which is the ratio between the

*x*−

*y*plane surface occupied by the fluid and the total

*x*−

*y*plane area.

*u*

_{τ}(Eq. (7)) is based on the difference between average pressure at the inlet and at the recycling plane as exposed in Fig. 6. The shear Reynolds number is then calculated as Re

_{τ}=

*u*

_{τ}

*H*/

*ν*. and the quantity $he+$ is defined as $he+=heu\tau /\nu $:

Concerning the friction factor, the Fanning definition *f* = 2(*u*_{τ}/*U*_{b})^{2} is used with *U*_{b} the bulk velocity. All these quantities are monitored at each iteration in the LES simulation.

*T*

_{p}is the temperature imposed at a wall and

*T*

_{∞}is the bulk temperature. This scalar is considered passive, and this hypothesis is valid if the temperature difference has no significant impact on the density, which is assumed here. This is why, the temperature can be replaced by this dimensionless scalar. The equation for this scalar is as follows:

*D*

_{z}includes the molecular and turbulent diffusivities.

### 3.2 Budget Equations.

In the present RRLES, different turbulence budget equations have been computed, dumped, and stored. This is the case for the mean kinetic equation (MKE) (Eq. (9)). Quantities denoted with the bracket 〈〉 symbol are time averaged, and *u*′ corresponds to the fluctuating component of velocity within the Reynolds decomposition. The quantity *ν*_{t} refers to the turbulent viscosity, which is modeled through the WALE model.

*MKE*:

*S*

_{ij}= (1/2)((∂

*u*

_{i}/∂

*x*

_{j}) + (∂

*u*

_{j}/∂

*x*

_{i})).

## 4 Results and Discussion

In this section, the wall-unit velocity and temperature profiles as well as the momentum and kinetic energy balance equations described in the previous section are analyzed. This analysis should provide a better understanding of the impact of the wall roughness on the flow and especially the effective slope parameter.

### 4.1 Impact on Turbulence.

A first qualitative analysis can be done by visualizing the vortices that are generated over the rough surface. The same Q-criterion iso-contours are plotted for the three different cases in Fig. 7 at a Reynolds number of 8000. At a first glance, it is clearly noticed that the parallel case promotes elongated vortices in the flow direction, while the two other cases also feature small vortices trapped between valleys. These trapped vortices are more coherent in the front case due to the orthogonal wavy pattern of the wall. Even if the roughness height distribution is the same, important differences are visually perceptible.

Mean streamwise velocity profiles are plotted in Fig. 8 for the three cases at Reynolds numbers of 8000 and 17,000 as well as the smooth wall configuration. The impact of ES is clearly visible as the velocity profile is shifted downward with the increasing ES. Differences are also noticeable between isotropic and anisotropic surfaces. Indeed, if the ES in isotropic case is lower than the front configuration, the mean velocity near the wall has a steeper increase for the isotropic case. These results indicate that the orientation of the wall roughness has a nonlinear impact on the mean velocity in the channel. Moreover, increasing the Reynolds number tends to shift downward the profile in the log-law region as expected.

Concerning the fluctuating velocity, streamwise and spanwise components are plotted for Reynolds numbers Re = 8000 and Re = 17,000 in Figs. 9 and 10, respectively. For the streamwise fluctuating velocity *u*_{rms}, and considering the two anisotropic cases, the ES has an impact on the maximum value, and this is latter shifted to the flow stream when changing the roughness direction. This peak is also less pronounced for higher ES value. As for the mean velocity, the isotropic case has a slightly different behavior, and while its ES value is in between the two anisotropic cases, the *u*_{rms} profile is not in between the two of the anisotropic cases confirming the nonlinear behavior of the rough wall with ES.

For spanwise velocity fluctuations, the maximum of the profile is less affected by the geometry, but its location follows the same trend as the streamwise fluctuations. Increasing the Reynolds number to 17, 000 tends to shift the peaks to the channel center, reduce the sharpness of the peaks for both velocity fluctuations, and slightly increase the maximum value of peaks for spanwise fluctuations.

The reduction of the anisotropy of the velocity fluctuations observed in each case in comparison with the smooth configuration is consistent with previous studies [9,24]. It has been shown that this is caused by the redistribution of the mean roughness wake kinetic energy into turbulence, generating additional normal and spanwise stresses [25]. The streamwise turbulence is also mainly converted into wake energy as the streaks are destroyed by the roughness elements. Increasing the Reynolds number contributes to intensify this phenomenon.

The presence of the wake is slightly visible on the mean velocity profiles when they become negative in the near-wall region for the isotropic and front cases. Mean kinetic energy balance plotted in Fig. 11(compared to Eq. (9) for each budget term) shows that this region is dominated by the work of the pressure drag against the flow direction. Those observations are no more valid for the parallel case since the surface topology does not create significant re-circulation zones, letting the viscous drag to be dominant. Thus, the decreasing of the turbulence anisotropy and the pressure loss are less pronounced for this case.

For each Reynolds number and all roughness types, the roughness velocity function Δ*U*^{+} is studied. The latter is evaluated at $he+=100$ and compared to the boundary layer log law. Results are summarized in Table 6. The first remark is that Δ*U*^{+} is amplified with higher Reynolds number. Values are mainly higher in the front case except at Re = 25,000 where the isotropic case has a larger value. These results point out that the additive manufacturing direction has a strong impact on the mean velocity profile although roughness height distributions are the same.

*f*

_{Darcy}= 4

*f*. Therefore, in Fig. 12, friction factors from the Colebrook equation are divided by a factor of 4 for getting comparative results. Two relative roughness

*S*

_{a}/

*D*

_{h}values (

*S*

_{a}/

*D*

_{h}= 10

^{−2}and

*S*

_{a}/

*D*

_{h}= 5.10

^{−3}) for Colebrook formula are plotted as the relative roughness of the present cases is intermediate:

*S*

_{a}/

*D*

_{h}= 7.8 10

^{−3}.

As expected, for each Reynolds number, the friction factor is higher with the increasing ES. For the front and isotropic cases, values are above the Colebrook correlation ones and tend to slightly increase with the Reynolds number. However, the trend for parallel case is reversed: the curve is below Colebrook expectations and decreases. There is a clear distinction in the friction factor behavior when varying the printing direction. Therefore, this correlation is inadequate for describing the tendency and evaluating a priori the friction factor. It also highlights again the need for better correlations for this type of roughness.

### 4.2 Impact on Temperature.

The roughness influence on heat transfer is also studied with this RRLES database. Indeed, modifications to the turbulent boundary layer induced by the roughness inhomogeneities can modify the temperature boundary layer and potentially the mean temperature field. As explained earlier, a passive scalar $Z\xaf$ is transported and averaged (Eq. (8)). Instantaneous $Z\xaf$ fields for the three rough cases are illustrated in Fig. 13. This figure shows that the temperature fluctuations observed close to the wall have wavelengths similar to those of the wall roughness. This is particularly visible in the front and isotropic cases. In the parallel case, the chosen cut cannot represent the wall roughness as the roughness features are aligned with the cut plane. As a result, less short wavelengths are visible in this cut close to the wall.

In order to quantify if temperature transport by turbulence is enhanced by wall roughness, mean temperature profiles are plotted for two different Reynolds numbers in Fig. 14. At Re = 8000, the temperature transport in the boundary layer is enhanced by wall roughness as the mean temperature gradient increases toward a constant value in the whole channel height. The three rough cases have a slightly different behavior close to the wall. It seems that the isotropic case is the configuration that is the closest to the smooth case. For the two anisotropic cases, the presence of coherent structures in the roughness valleys seems to slightly increase the temperature wall gradient and thus the wall heat transfer.

For Re = 17,000, the conclusions are different and more intuitive. Indeed, the parallel roughness and smooth cases have almost the same temperature profiles as it would be expected due to the alignment of the roughness elements with the flow. The front and isotropic cases have also the same temperature profile with a clear shift toward the constant gradient imaginary profile. In these cases, the wall roughness promotes temperature transport by turbulence. For this higher Reynolds case, the dependence of the wall heat transfer with the effective slope is better recovered. Again, the different nonlinear behaviors of the temperature with both the effective slope and the Reynolds number illustrate the challenge of finding a Nusselt number correlation for rough walls.

## 5 Conclusions

This article presents roughness-resolved large-eddy simulations that are representative of the flow obtained in additive-manufactured heat exchangers. The aim of these simulations is to provide a rich database that will ease the development and assessment of rough-wall models. The most challenging task in the building of this database is the generation of representative rough surfaces and conformal unstructured meshes. Three different configurations of parallel rough plane channels with the same roughness height distribution but different effective slopes have been chosen and modeled. The impact of the effective slope parameter, which is directly linked to the alignment of the wall roughness with the flow, has a strong impact on the flow topology, velocity, and temperature profiles, as expected. In these cases, the existing empirical correlations find their limits, and new correlations are needed.

## Acknowledgment

The STREAM project has received funding from the Clean Sky 2 Joint Undertaking (JU) (Grant No. 865378). The JU receives support from the European Union’s Horizon 2020 research and innovation program and the Clean Sky 2 JU members other than the Union. This article reflects only the authors’ view, the Joint Undertaking is not responsible for any use that may be made of the information it contains. This work was granted access to the HPC resources of CINES/TGCC/IDRIS under the allocation 2020-A0092B06880 and 2021-A0112B06880 made by GENCI and of CRIANN under the allocation 2012006.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

*f*=fanning friction factor

*H*=half height channel

- $Z\xaf$ =
passive scalar

- Re
_{τ}= shear Reynolds number

*S*_{a}=arithmetic average height

*S*_{q}=root-mean-square height

*S*_{h}=maximum distance between peaks and valleys

- AM =
additive manufacturing

- (C)HX =
(compact) heat exchanger

- ES =
effective slope

- Ku =
Kurtosis

- Pr =
Prandtl number

- Re =
Reynolds number

- Sk =
skewness

## Appendix: Numerical Methods Assessment

We validate our numerical methods by reproducing the classical Re_{τ} = 180 pressure-gradient driven periodic channel flow DNS test case from Ref. [27] (refered to as KMM hereafter). The geometry is *L*_{x} = 4*πH*, *L*_{y} = 2*πH*, *L*_{h} = 2*H* with the same RRLES direction denomination, and the dimensionless mesh resolution is $\Delta x+=8.6$, $\Delta y+=8$, $\Delta h+=[0.38,5.3]$. A first test (T1) is done with the periodic condition methodology, and a second one (T2) uses the recycling boundary condition. For the latter, we set the recycling plane at a distance 2*H* above the outlet in the longitudinal direction and the input flowrate is chosen to impose the same bulk velocity measured in Ref. [27]. Results are summarized in Fig. 15. Good agreement between tests and references shows that our numerical schemes and recycling boundary condition are appropriate to perform infinite channel flow simulations.

We discuss now the RRLES mesh size used in this article. Figure 16 compares the results obtained from the Re = 17,000 isotropic configuration with the mesh described in Table 2 (denoted M2) and those obtained with a coarser mesh (M1) whose cell size is 10 *μ*m. It is seen that no significant differences may be observed between the two simulations, on both velocity first- and second-order moments. This suggests that results will be probably weakly improved with a mesh finer than M2, for the analyzed moments.

Thus, we consider that the dimensionless wall normal resolution *y*^{+} obtained with the mesh M2 for this case is a satisfactory reference. Due to the surface alternating of peaks and valleys, making the wall normal mean velocity gradient to vary in space, it is more relevant to observe the statistical distribution of the wall resolution than its mean value alone to evaluate the mesh quality. Figure 17 shows the probability density function of *y*^{+} obtained for the considered Reynolds numbers. Obviously, the distributions for the two lowest Reynolds numbers are satisfactory. On the other hand, the distribution for Re = 25,000 is larger than our limit. Nevertheless, it is similar to the one obtained for Re = 17,000 with M1. According to the previous results, we can consider that M2 is adequate for the simulations performed in this article. The analysis for the front and parallel surfaces (not shown here) leads to the same conclusions.