# Analysis of separation in the roughness sublayer using DNS data and DANS/DEM modelling of roughness effects

François Chedevergne, Jiasheng Yang, Alexander Stroh, Pourya Forooghi

► To cite this version:

François Chedevergne, Jiasheng Yang, Alexander Stroh, Pourya Forooghi. Analysis of separation in the roughness sublayer using DNS data and DANS/DEM modelling of roughness effects. 2025. hal-04921678

HAL Id: hal-04921678

https://hal.science/hal-04921678v1

Preprint submitted on 30 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Analysis of separation in the roughness sublayer
using DNS data and DANS/DEM modelling of
roughness effects

François Chedevergne¹*, Jiasheng Yang²†, Alexander Stroh²†,
Pourya Forooghi³†

¹DMPE, ONERA, Université de Toulouse, Toulouse, France.

²Institute of Fluid Mechanics (ISTM), Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany.

³Department of Mechanical and Production Engineering, Aarhus University, Aarhus, Denmark.

*Corresponding author(s). E-mail(s): francois.chedevergne@onera.fr
ORCID:0000-0003-3634-8482;

Contributing authors: jiasheng.yang@kit.edu ORCID:0009-0007-1370-9085; alexander.stroh@kit.edu ORCID:0000-0003-0850-9883; forooghi@mpe.au.dk ORCID:0000-0001-9212-514X;

†These authors contributed equally to this work.

Abstract

From the recent DNS database (Yang et al., 2022) of channel flows with rough walls in the presence of heat transfer, the impact of the skewness of the roughness elevation map on the velocity and temperature profiles within the roughness sublayer is analysed. The separation zones observed near the wall in the sublayer are shown to play a significant role when the skewness is negative. The *k*-ω-based turbulence model (Chedevergne and Forooghi, 2020; Chedevergne, 2021, 2023), capable of capturing roughness effects and incorporating the Double Averaged Navier-Stokes (DANS) equations and the Discrete Element Method (DEM), is tested against this DNS database, showing some limitations in the description of the roughness sublayers, especially for configurations with negative skewness. To reproduce the observations made in the DNS database, the pressure gradient imposed in the simulated channel using the DANS/DEM model is adjusted

distributed under license CC-BY 4.0 https://creativecommons.org/licenses/by/4.0/
based on the distance to a reference wall in the roughness sublayers. Additionally, the increase in turbulent mixing observed in the DNS database for rough configurations with negative skewness is accounted for in the DANS/DEM model by modifying the source terms in the transport equations of the turbulent scalars with respect to the skewness, improving the prediction of the roughness effects.

**Keywords:** Double Averaged Navier-Stokes, Discrete Element Method, roughness effects

# 1 Introduction

The prediction of roughness effects in fluid mechanics simulations has long been identified as a challenge and has been a subject of study for many researchers since the founding work of Nikuradse (1937) and Schlichting (1937). There has been renewed interest in this issue recently, not least because of the growing importance of additive manufacturing technologies in industry. Although for some years now direct numerical simulations (DNS) databases (Busse et al., 2015; Forooghi et al., 2018; Thakkar et al., 2018; Jelly and Busse, 2019) give access to a fine description of the phenomena involved in wall-bounded flows in the presence of roughness, many questions remain unanswered. In particular, it is known (Flack and Schultz, 2010; Flack et al., 2019) from experimental work that the skewness $Sk$ of the roughness distribution plays an important role in the magnitude of the roughness effects produced on the flow. However, we do not know what mechanisms are at stake behind these findings. The DNS database established by Yang et al. (2022) very recently on a series of arbitrarily generated rough surface distributions makes it possible to study the influence of $Sk$ on roughness effects. This database is all the more valuable in that it also incorporates the effects of roughness on heat transfer, which is quite rare. The first objective of this study is therefore to analyse these DNS data in the region of the roughness sublayer in order to gain a better understanding of the effects induced by roughness, particularly when $Sk$ becomes negative.

In parallel with the recent advances made thanks to DNS, the modelling of roughness effects has greatly evolved over the last decade. The second objective of this study is to evaluate and improve if necessary a very recently proposed model (Chedevergne, 2021) that accurately describes the effects of roughness in a large number of situations, including a good description of flows in rough sub-layers. The comparisons of the model predictions with the DNS data of Yang et al. (2022) should highlight the limitations of the model and allow to propose appropriate corrections.

Historically, to take roughness effects into account in computational fluid dynamics, several approaches may be contemplated, the most widespread in an industrial context being that based on the resolution of the RANS (Reynolds Averaged Navier-Stokes) equations and the joint use of equivalent sand-grain type models (Wilcox, 1988; Hellsten and Laine, 1998; Knopp et al., 2009; Aupoix, 2015b) associated with eddy viscosity models. These models are based on Nikuradse's work, which highlighted the
link between the roughness function, which characterizes roughness effects on veloc-
ity profiles, and the sand grain height used in Nikuradse’s experiments (Nikuradse,
1937). Although undeniable advantages of these equivalent sand-grain models should
be acknowledged, such as their ease of implementation and robustness, there are still
major limitations in the use of these models due to their poor predictive capacity. On
the one hand, these models produce an artificial increase in turbulence in the wall
region and, on the other, the a priori determination of the main parameter of these
models, i.e. the equivalent sand-grain height, is extremely delicate. The sand-grain
height is given by empirical correlations (Musker, 1980-81; Sigal and Danberg, 1990;
van Rij et al., 2002; Schultz and Flack, 2009) that are still limited to a certain type of
application. Despite notable efforts (Forooghi et al., 2017) to come up with a universal
correlation, the quest for such an expression relating the rough surface topography to
the height of the equivalent sand grain remains unachieved and is, in some respects,
unachievable. More specifically, in the transitional rough regime, flow properties inter-
act in the relationship between the roughness function and the equivalent sand grain,
making the discovery of a universal correlation impossible in this regime.

In the RANS context, another approach, called the discrete element method (DEM), and this time based on the blockage factor concept put forward by Schlichting (1937), finally emerged from the late 1970s (Finson, 1975; Finson and Clarke, 1980) and into the 1980s (Finson, 1982; Lin and Bywater, 1982; Christoph and Pletcher, 1983; Taylor et al., 1985). The DEM enables to overcome the main limitations of equivalent sand-grain models and the most recent advances in this approach mean that it can now be used in a production environment. More specifically, the major development consisted into the coupling of the DEM with RANS equations with additional surface averaging. Therefore, the resulting equations (Aupoix, 2016), hereinafter referred to as DANS (Double Averaged Navier-Stokes), include a drag term that can be described using the discrete element approach. Accurate models for the drag term such as those of Chedevergne and Forooghi (2020) and of Kuwata and Kawaguchi (2019) compare favourably with DNS data. These models include the wall-normal distributions of both the viscous and pressure drag contributions throughout the boundary layer height covered by roughness. To complete the DANS/DEM models, dedicated closures for Reynolds stresses and dispersive stresses arising from the double averaging technique must be added. Comparisons with DNS data over rough academic configurations (Forooghi et al., 2018; Wu et al., 2019; Kuwata and Kawaguchi, 2019) at moderate Reynolds numbers proved the relevance of this type of approach to describe roughness effects with DANS/DEM models (Chedevergne, 2021; Kuwata et al., 2019). Lastly, there remains the delicate issue of determining the input parameters for these DANS/DEM approaches. As mentioned above, it is this link between rough surface topography and model input parameters that constitutes the main limitation of equivalent sand-grain approaches. Here, in the case of the discrete element method, input parameters may be purely of geometric origin. Introducing the RER concept (Representative Elementary Roughness), Chedevergne (2023) proposed an efficient methodology to determine the various input parameters from any rough surface topography. The methodology has thus made it possible to extend the validation of Chedevergne’s DANS/DEM model to high Reynolds number experimental configurations obtained on surfaces covered
with sandpaper.

To reach the two objectives described above, the DNS database of Yang et al. (2022) was thoroughly analysed to point out the influence of Sk on the wall transfer in presence of roughness. Section 2 describes the DNS database (Yang et al., 2022) and discusses identified specific points. Section 3 presents the DANS/DEM model (Chedevergne, 2021). To best reproduce the behaviors observed in the DNS, appropriate modifications in the DANS/DEM model (Chedevergne, 2021) were carried out. The obtained results are discussed in section 4 where the proposed modifications and outstanding questions are exposed.

## 2 DNS database description

The current roughness database is constructed through Direct Numerical Simulations (DNS) that are carried out in the preceding research of Yang et al. (2022). In this work, the concept of minimal-channel DNS (Chung et al., 2015; MacDonald et al., 2017) is employed and systematically investigated in terms of its applicability for the prediction of the dynamic and thermal roughness effects on irregular rough surfaces. The minimal channel approach enables to accurately predict the roughness functions $\Delta u^+$ and $\Delta \theta^+$ (Yang et al., 2023) for the velocity and temperature profiles respectively, and also the zero-plane displacement $\xi$. Nevertheless, it is noteworthy that the mean profiles depart from the ones obtained in full-sized DNS near the centerline of the channel, attributing to the nature of the minimal channels. Therefore, to facilitate comparisons with the current computations across the entire height of the channel, we only retain the simulations performed on a full-sized domain from all the computations done, which served as reference cases for the study by Yang et al. (2022).

The DNS of incompressible turbulent flow over irregular roughness are performed in fully developed turbulent plane channels, where the flow is driven by a constant pressure gradient (CPG). The roughness are installed on both the upper and lower walls. The channel half-height, defined as the distance between the deepest trough in the roughness and the center plane of the channel, is denoted as $H$ and will be used to normalize the geometric scales. The geometry of the simulation domain is $L_x \times L_y \times L_z = 8H \times 2H \times 4H$. Periodic boundary conditions are applied in the streamwise and spanwise directions. The Immersed Boundary Method (IBM), based on the approach described in Goldstein et al. (1993), is employed to impose the no-slip and constant temperature boundary condition on the rough surfaces by introducing an external volume force field $f_{IBM}$ and $f_{\theta}$ directly into the momentum and energy equations, respectively. The Navier-Stokes equation and the energy conservation equation write:

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot (\mathbf{u}\mathbf{u}) = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} - \frac{1}{\rho}\Pi\mathbf{e}_x + \mathbf{f}_{IBM}, \qquad (2) $$

$$ \frac{\partial \theta}{\partial t} + \nabla \cdot (\mathbf{u}\theta) = \alpha \nabla^2 \theta + Q + \mathbf{f}_{\theta}. \qquad (3) $$
Where $\mathbf{u} = (u, v, w)^{\top}$ is the velocity vector and $\Pi$ is the CPG source term to drive the flow. Here $p, \mathbf{e}_x, \rho$ and $\nu$ are pressure fluctuation, streamwise basis vector, density and kinematic viscosity, respectively. The friction Reynolds number is defined as $\text{Re}_{\tau} = u_{\tau}(H - k_{\text{md}})/\nu$, where $u_{\tau} = \sqrt{\tau_w}/\rho$ and $\tau_w = -\Pi \cdot (H - k_{\text{md}})$ are the friction velocity and wall shear stress, respectively. $\theta$ represents the relative wall temperature, namely $\theta = T - T_w$. $Q = -ud\theta/dx$ is the temperature source term to achieve developed temperature field. The constant Prandtl number is set to $Pr = 0.7$. In the present DNS database, $\Pi$ is selected to achieve $\text{Re}_{\tau} = 500$ across all the considered cases. The roughness melt-down height denoted by $k_{\text{md}}$ is the mean roughness height measured from the deepest trough. The above equations are solved using the pseudo-spectral solver SIMSON (Chevalier et al., 2007), where wall-parallel directions are discretized in Fourier space, while in wall-normal direction Chebyshev discretization is employed. The grid resolution of the simulation is designed to fully resolve the smallest roughness wavelength present. Typically, eight grid points in each wall-parallel direction are used to capture the smallest roughness wavelength, i.e., $\Delta_{x,z} \le \lambda_{\text{min}}/8$, where $\lambda_{\text{min}}$ denotes the smallest roughness wavelength in the artificially generated roughness. Additionally, the simulation resolution for wall-bounded turbulent flow is compared to the viscous length scale, ensuring $\Delta_{x,z}^{+} < 5$ and $\Delta_{y,k}^{+} < 2$, where $\Delta_{y,k}^{+}$ represents the inner-scaled grid size at the roughness peak. Consequently, the simulation resolution for the present DNS cases is $N_x \times N_y \times N_z = 900 \times 401 \times 480$.

Yang et al. generated artificial roughness leveraging the methodology proposed by Prez-R fols and Almqvist (2019), where the wall-parallel properties of the topography is controlled through the calibration of the associated in-plane roughness height power spectrum (PS), and the wall-normal properties of the topography using roughness height probability density function (p.d.f.)

Specifically, three types of p.d.f. – distinguished by the values of the skewness $Sk$ – were utilized. On the other hand, the roughness PS is characterized by two power-law slope $p = -1$ or $-2$, therefore the PS is formulated as a function of the wavenumber $q$, i.e. $q^p$. For a given p.d.f., varying the slope of the PS effectively means giving more or less importance to longer wavelengths, and ultimately altering the effective slope $\phi$ of the rough surface generated. According to Schultz and Flack (2009) the effective slope has a significant impact on the roughness effects. The final parameter varied by Yang et al. is the position $\lambda_c$ of the long-wavelength cut-off in the PS. This parameter, as expected, can directly influence the effective slope of the generated surface as well. Through grid-searching the combinations of these parameters, twelve different types of roughness with varying configurations were generated and simulated.

In the following content, the naming convention $Xnm$ is used to designate each of these twelve configurations. Herein, $X$ is a character indicating the type of p.d.f. (G for Gaussian, i.e. $Sk = 0$, $P$ for positive skewness, i.e. $Sk = 0.48$, and $N$ for negative skewness, i.e. $Sk = -0.48$). The digits $n$ and $m$ represent the opposite of the PS slope, i.e. $n = -p$ and the large cutoff wavelength ($m = 5\lambda_c/h$, where $h$ is the channel height), respectively. Consequently, the characteristics of the roughness configurations associated with each topography are summarized in table 1. Flow characteristics are given through $\Delta u^+$ the roughness function, $Re_b$ the bulk Reynolds number, $\xi$ the distance
<table><caption>**Table 1:** Characteristics of the roughness configurations</caption><thead><tr><th>Topography</th><th>Δu<sup>+</sup></th><th>Re<sub>b</sub></th><th>ξ<sup>+</sup></th><th>h<sup>+</sup></th><th>k<sub>r</sub><sup>+</sup></th><th>φ</th><th>y<sub>th</sub>/k<sub>r</sub></th><th>β<sub>m</sub></th><th>d<sub>M</sub><sup>+</sup></th></tr></thead><tbody><tr><td>P14</td><td>7.33</td><td>5073</td><td>40.74</td><td>524.80</td><td>60.17</td><td>0.54</td><td>0.36</td><td>0.48</td><td>61.87</td></tr><tr><td>P18</td><td>7.23</td><td>5102</td><td>40.52</td><td>532.48</td><td>61.91</td><td>0.52</td><td>0.36</td><td>0.48</td><td>65.52</td></tr><tr><td>P24</td><td>6.99</td><td>5234</td><td>39.03</td><td>533.80</td><td>62.06</td><td>0.42</td><td>0.36</td><td>0.48</td><td>85.91</td></tr><tr><td>P28</td><td>6.57</td><td>5433</td><td>37.80</td><td>532.85</td><td>61.95</td><td>0.37</td><td>0.36</td><td>0.48</td><td>95.86</td></tr><tr><td>G14</td><td>6.67</td><td>5317</td><td>47.34</td><td>534.39</td><td>61.95</td><td>0.54</td><td>0.50</td><td>0.51</td><td>58.54</td></tr><tr><td>G18</td><td>6.56</td><td>5345</td><td>47.25</td><td>532.48</td><td>61.91</td><td>0.52</td><td>0.50</td><td>0.51</td><td>60.46</td></tr><tr><td>G24</td><td>6.30</td><td>5489</td><td>45.80</td><td>533.80</td><td>62.06</td><td>0.41</td><td>0.50</td><td>0.51</td><td>78.24</td></tr><tr><td>G28</td><td>5.94</td><td>5661</td><td>44.90</td><td>532.85</td><td>61.95</td><td>0.36</td><td>0.50</td><td>0.51</td><td>89.09</td></tr><tr><td>N14</td><td>6.14</td><td>5495</td><td>53.10</td><td>539.83</td><td>58.63</td><td>0.56</td><td>0.66</td><td>0.54</td><td>55.58</td></tr><tr><td>N18</td><td>5.84</td><td>5510</td><td>52.83</td><td>540.24</td><td>58.68</td><td>0.52</td><td>0.66</td><td>0.54</td><td>57.38</td></tr><tr><td>N24</td><td>6.09</td><td>5661</td><td>51.48</td><td>540.94</td><td>58.76</td><td>0.42</td><td>0.66</td><td>0.54</td><td>74.54</td></tr><tr><td>N28</td><td>5.51</td><td>5801</td><td>50.63</td><td>539.83</td><td>58.59</td><td>0.37</td><td>0.66</td><td>0.54</td><td>84.16</td></tr></tbody></table>

to the zero-displacement plane and $h^+$ the dimensionless channel height. The roughness distribution is characterized by $k_r$ the maximum elevation, $\phi$ the effective slope while $y_{th}$ the threshold distance, $\beta_m$ the mean blockage and $d_M$ the maximum diameter are parameters characterising the roughness in the DANS/DEM model described in section 3. The topography of the G24 case is shown in figure 1 as an instance. The upper section of the figure provides a zoomed-in view that includes the averaged longitudinal velocity distribution $\langle u \rangle^{f1}$. This underscores the presence of near-wall zones where, on average, flow recirculation can be observed. Additionally, the isocontour at $y = y_0$ is plotted in blue, with $y_0$ such that $\langle u \rangle (y_0) = 0$. The three types of p.d.f. used in the present work and the roughness profiles for the surfaces P24, G24 and N24 along a y-z cutting plane are displayed in figure 2. Furthermore, the impact of the variation in skewness on the effective slope $\phi$ as well as on the position of the zero-displacement plane $\xi$ is evident. However, it is noteworthy that $\xi$ is evaluated a posteriori via Jackson's method (Jackson, 1981) and therefore depends on the simulated flow and not merely on the geometry of the rough wall.

Yang et al. (2022) investigated the distribution of surface force with the help of the sheltering effect. By revisiting the original sheltering model developed by Yang et al. (2016), initially adapted to rectangular prismatic rough elements, Yang et al. (2022) noted that a significant portion of the rough surface does not substantially contribute to drag force due to its location in the "shadowed areas" of the up-stream extreme peaks. This central observation is to be considered in the light of the concepts already introduced by McClain et al. (2006) and taken up again by Chedevergne (2023), where identifying the so-called "dead water zones" is crucial to avoid their contribution to the drag penalty prediction. The sheltering effect is illustrated in the left part of the figure 3 for the cases P24, G24 and N24. A key point to be drawn from these observations is that the longitudinal force distribution $\langle F_{d_x} \rangle$ is not necessarily correlated with the mean velocity profile $\langle u \rangle$. We can see from the central part of the figure 3 that the drag force increases in absolute value while the mean velocity is

<footnote n="1">¹Notations are defined in section 3</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1:** Surface topography for configuration G24 ($Sk = 0$). The time-averaged velocity profile and the isocontour indicating $\langle u \rangle (y_0) = 0$ are displayed in the zoomed-in subfigure.</figcaption></figure>

negative. We recall that $y_0$ is the position where $\langle u \rangle (y_0) = 0$. Close to the wall, for values of $y$ less than $y_0$, a very large part of the roughness distribution is sheltered. The flow is mainly separated, which explains why $\langle u \rangle$ is negative, but in the downstream part of these sheltered zones the flow can locally be oriented downstream, generating force contributions $F_{d_x}$ that are opposite to the mean flow direction. Interestingly, we observe a different behavior for the vertical component $\langle F_{d_y} \rangle$ of the roughness force which correlates well with the mean longitudinal velocity $\langle u \rangle$, with the force $\langle F_{d_y} \rangle$ only starting to increase for $y > y_0$. In the right part of the figure 3, the impinging point locations of the sheltering effect model, corresponding to the initial locations of the green parts in the roughness distributions, were collected to plot the associated probability functions. It clearly shows that some important contributions to $\langle F_{d_x} \rangle$ can be found below $y_0$.
<figure><img src="image_3.png" /><figcaption>**Fig. 2:** Probability density function for the three configurations G24, P24 and N24 and one-dimensional roughness profiles along a constant *z* with dashed lines indicating the threshold values *y*<sub>th</sub>.</figcaption></figure>

Note that the data are extracted from computations obtained on minimal channels <fnref n="2" />
and for which we observed slight differences in the separation zones compared to
computations on full channels, particularly for the y₀ location. However, conclusions
drawn from these data are equally valid for full-channel computations discussed in the
rest of the paper.

3 System of equations and the Representative
Elementary Roughness

From the RANS equations, by application of the surface averaging, also sometimes
referred as the volume averaging (Whitaker, 1986), the DANS equations are obtained.
In the following, $\bar{q}$ and $q'$ designate the mean and fluctuating part of any physical
quantity $q$ obtained with the Reynolds averaging, while $\langle q \rangle^f$ and $\tilde{q}$ correspond to the
intrinsic surface averaging. The surface averaging $\langle q \rangle^f$ is given by $\frac{1}{S_f} \int_{S_f} q ds$, where
$S_f$ is the part of the surface $S$ open to the flow. $\langle \Box \rangle^f$ is called the intrinsic average
while the standard average $\langle \Box \rangle$ can easily be found using the relationship $\langle \Box \rangle^f = \beta \langle \Box \rangle$.

<footnote n="2"><sup>2</sup>Case M2 detailed in Yang et al. (2022). Extraction of the force components was not realized for computations on full channels.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3:** Normalized force $f_x$ and roughness distribution profile for N24 (upper), G24 (middle) and P24 (lower) at $z = 0.2H$. Sheltering effect modelled by Yang et al. (2016) is illustrated by gray dashed lines, the unsheltered surfaces are marked by green lines. Velocity profiles are represented by blue lines on the middle column, the fragment of positive velocity being illustrated by dashed line. The wall normal position $y^+$ where the mean velocity profile intersect $\langle u \rangle^+ = 0$ is marked with a black dashed dotted line, angle bracket $\langle \cdot \rangle$ indicates wall-parallel averaging (see section 3). The p.d.f. distributions of height of the impinging point of the sheltering are shown on the right column.</figcaption></figure>

The surface ratio $\beta = S_f/S$ is the blockage factor due to the roughness and is of critical importance in the DANS equations. For an incompressible fluid, the resulting
continuity, momentum and energy equations read:

$$ \frac{\partial \langle \beta \overline{u_i} \rangle^f}{\partial x_i} = 0 \\ \frac{\partial \langle \overline{u_i} \rangle^f}{\partial t} + \langle \overline{u_j} \rangle^f \frac{\partial \langle \overline{u_i} \rangle^f}{\partial x_j} = -\frac{1}{\rho} \frac{\partial \langle \overline{p} \rangle^f}{\partial x_i} + \frac{1}{\beta} \frac{\partial}{\partial x_j} \left( \nu \frac{\partial \beta \langle \overline{u_i} \rangle^f}{\partial x_j} \right) \\ -\frac{1}{\beta} \frac{\partial}{\partial x_j} \left( \beta \langle \overline{u_i' u_j'} \rangle^f + \beta \langle \overline{u_i \overline{u_j}} \rangle^f \right) - \frac{\nu}{\beta} \frac{\partial \beta}{\partial x_j} \frac{\partial \langle \overline{u_i} \rangle^f}{\partial x_j} - \overline{F_{d_i}} \quad (4) $$

$$ \frac{\partial \langle \overline{h_t} \rangle^f}{\partial t} + \langle \overline{u_j} \rangle^f \frac{\partial \langle \overline{h_t} \rangle^f}{\partial x_j} = \frac{1}{\rho} \frac{\partial \langle \overline{p} \rangle^f}{\partial t} + \frac{1}{\beta} \frac{\partial}{\partial x_j} \left( \nu \langle \overline{u_i} \rangle^f \frac{\partial \beta \langle \overline{u_i} \rangle^f}{\partial x_j} + \kappa \frac{\partial \beta \langle \overline{h_t} \rangle^f}{\partial x_j} \right) \\ -\frac{1}{\beta} \frac{\partial}{\partial x_j} \left( \beta \langle \overline{u_i} \rangle^f \langle \overline{u_i' u_j'} \rangle^f + \beta \langle \overline{h_t' u_j'} \rangle^f - \beta \langle \overline{h_t \tilde{u}_j} \rangle^f \right) - \frac{\kappa}{\beta} \frac{\partial \beta}{\partial x_j} \frac{\partial \langle \overline{h_t} \rangle^f}{\partial x_j} - \overline{F_t} $$

where $u_i$ is the velocity component in the direction $x_i$, $p$ the pressure, $h_t$ the total enthalpy, $\rho$ the density and $\nu$ and $\kappa$ the kinematic viscosity and the thermal diffusivity, respectively. The force $F_{d_i}$ due to roughness and the corresponding energy source term $F_t$ are given by:

$$ F_{d_i} = \frac{1}{S_f} \left( \int_{I_{fs}} \frac{\tilde{p}}{\rho} n_i ds - \int_{I_{fs}} \nu \frac{\partial \tilde{u}_i}{\partial x_k} n_k ds \right); \quad F_t = - \frac{1}{S_f} \int_{I_{fs}} \kappa \frac{\partial \tilde{h}_t}{\partial x_k} n_k ds \quad (5) $$

$I_{fs}$ designates the fluid/solid interface. In the following, the $\bar{\square}$ symbol for Reynolds averaged quantities is dropped to lighten the notation. More details on the derivation of the double averaged equations can be found in Aupoix (2016).

The streamwise component of the force $F_{d_i}$ is the drag force and it is modelled through a drag coefficient $C_d$. For channel flows, such as those considered in the DNS of Yang et al. (2022), the longitudinal drag component $F_{d_x}$ is related to the drag coefficient $C_d$ according to the relation:

$$ F_{d_x} = (1 - \beta) \frac{2\nu}{\pi d^2} C_d \|\langle u \rangle\|^f \|\langle u \rangle\|^f = f_d \langle u \rangle^f \quad (6) $$

We introduce $f_d$, a frequency related to drag, which will be reused later in turbulent closures.

In order to develop the different closures induced by the introduction of the Reynolds and volume averagings, the roughness distribution under consideration is reduced to limited number of parameters to control the closure relations, following the principles of the discrete element method alluded in the introduction. Therefore, Chedergne (2023) introduced the notion of Representative Elementary Roughness (RER), enabling a roughness distribution to reduce to only two independent parameters. From the digitalized topography of a rough surface and by defining a reference plane, simple image processings enable to extract the distribution in the normal direction of the
blockage factor $\beta$ and of the wetted surface $P_w$, i.e. the length of the fluid/solid interface in each plane parallel to the reference plane. With the assumption that the RER is a rough element with a circular cross-section, the wetted parameter is then used to compute the diameter $d$ of the RER at each altitude $y$. The spacing between RER elements $L_x \times L_z$ is supposed to be squared, i.e. $L_x = L_z$. Given the relationship existing between the spacing $L_x \times L_z$, the diameter $d$ and the blockage factor $\beta$, which reads $\beta = 1 - \frac{\pi d^2}{4L_x L_z}$, the two last parameters of the RER, $L_x$ and $L_z$, can be computed.

At last, the RER definition requires the position of a threshold, denoted $y_{th}$, below which the RER properties are considered to be fixed. Because there may be near-wall areas which do not participate in transfer, we have introduced the threshold position to account for this characteristics of the roughness distributions. As an example, in figure 1 we can expect that below the iso-contour $y = y_0$ there is negligible contributions of the roughness distributions. In that sense, $y_{th}$ can be viewed as an estimate of $y_0$. Practically, the threshold value $y_{th}$ chosen to be equal to the position of the wetted perimeter maxima. The positions of the thresholds are plotted on rough profiles in figure 1 for three different configurations. Obviously, the type of p.d.f. used for a given rough configuration, i.e. the associated skewness value, has a direct influence on the threshold value $y_{th}$ (see table 1).

The closure of equations set (4) starts with the modelling of the drag coefficient $C_d$. From a series of experiments by Žukauskas (1972), Chedevergne and Forooghi (2020) derived a local expression for the drag coefficient $C_d$. Žukauskas data were obtained by measuring the pressure drop in channels equipped with identical vertical tubes opposing to the upstream flow. The data retained to establish the expression of $C_d$ correspond to a staggered arrangement of circular cylinders. The main conclusion from these experiments is that the pressure drop is firstly controlled by the Reynolds number $Re_d$ defined with the diameter $d$ of the cylinders and the spatially averaged velocity noted $u$. The pressure drop is then directly related to the drag coefficient $C_d$ generated by the RER. Indeed, the RER is assimilated to a superposition of independent infinitesimal cylinders with a local diameter $d$. In other words, the model proposed by Chedevergne and Forooghi (2020) neglects the curvature of the RER in the wall-normal direction and only retain the evolution of the diameter $d$ in the wall-normal direction. The drag coefficient $C_d$ is given by:

$$C_d = 1.5 (\alpha\beta)^2 \xi \qquad (7)$$

In the expression (7), $\alpha$ is a geometrical parameter of the rough surface, given by $\alpha = \frac{L_z}{L_z - d}$, where $d$ is the diameter of the RER, and $L_z$ the transverse spacing between roughness elements. Function $\xi$, characterizing the pressure drop in Žukauskas experiments, is shown (Chedevergne and Forooghi, 2020) to only depends on the Reynolds
number $Re_d = \frac{\langle u \rangle^f d}{\nu}$:

$$ \log \xi = (0.58f - 0.86) \log (\beta Re_d) + 1.82 - 1.1f, \forall \beta Re_d \le 116883 \quad (8) $$

$$ \xi = 0.2, \forall \beta Re_d > 116883 $$

with f a blending function:

$$ \text{if } \beta Re_d \le 60, f = 0 $$

$$ \text{else if } 60 < \beta Re_d \le 200, f = 1 - \frac{60}{\beta Re_d} \quad (9) $$

$$ \text{else when } 200 < \beta Re_d \le 116883, f = 1 - \frac{1}{\left(\frac{\beta Re_d}{10}\right)^{0.4} + \left(\frac{\beta Re_d}{10000}\right)^{2.78}} $$

The description and the validation of the drag model are detailed in Chedevergne and Forooghi (2020). To fully appreciate the limitations of the model presented here, it is essential to note that the drag term (5), which is by definition a surface integral, is taken into account here via a volume source term. In other words, the "local" contributions of the drag term are distributed homogeneously and continuously over the entire volume occupied by the roughness. Furthermore, the use of Zukauskas' experimental data leads to a model for the drag force, which is necessarily directed in the direction of the flow. As discussed in section 2 and illustrated with figure 3, this is not necessarily the case with DNS data. Zukauskas' data (Žukauskas, 1972) correspond to pressure drop measurements in a channel in which circular cylinders are arranged. In these experiments, the force exerted on the cylinders is then aligned with the upstream velocity. Consequently, the resulting model of Chedevergne and Forooghi (2020) can only reflect this situation. This is reflected in the relationship (6) which shows that the sign of the force $F_{d_x}$ is directly given by that of the velocity of the flow $\langle u \rangle^f$.

An eddy viscosity approach involving the Boussinesq hypothesis is used to model the Reynolds stresses that entering eq. (4). The $k-\omega$ SST model of Menter (1994) was selected as a base model and enhanced with additional source terms $F_k$ and $F_\omega$ in the respective transport equations of the turbulent scalars $k$ and $\omega$, to capture the effects of roughness on turbulence.

$$ F_k = c_k f_d k \quad ; \quad F_\omega = c_{\omega_s} \left( c_{\omega_1} + c_{\omega_2} \frac{1-\varphi}{\varphi^{c'_{\omega_2}}} \right) \left( 1 + c_{\omega_3} e^{\frac{-1}{1-\beta_m}} \right) f_d \omega \quad (10) $$

Constants are $c_k = -2$, $c_{\omega_1} = 2.7$, $c_{\omega_2} = 0.1$, $c'_{\omega_2} = 1$ and $c_{\omega_3} = 30$. Using $f_d$ to model the source terms $F_k$ and $F_\omega$ allows us to propose simple closures, which are patterned on the behavior of the drag term. In particular, this implies that the source terms
no longer act above the roughness crest. The expression for $F_\omega$ (10) is depending on the mean blockage factor $\beta_m$, computed from the blockage evolution in the RER, and the effective slope $\phi$, reflecting the key role played by these two parameters on the turbulence induced by roughness. During the model elaboration (Chedevergne, 2021), the rough surfaces considered have led to fix $c'_{\omega_2}$ to 1. Since these considered topography were composed of academic rough surfaces, it is anticipated that when considering more realistic rough distributions, slight modifications of the constant $c'_{\omega_2}$ may be beneficial. Similarly, an additional coefficient $c_{\omega_s}$ is incorporated in this study to address the impact of negative skewness values on turbulence induced by roughness. This coefficient is always set to 1 for configurations Pnm and Gnm where the skewness is positive or null.

The eddy viscosity is deduced from $k$ and $\omega$ as done in the SST model Menter (1994) and the Reynolds stresses are deduced from the Boussinesq relation

$$ -\langle u'_i u'_j \rangle^f = \nu_t \left( \frac{\partial \langle u_i \rangle^f}{\partial x_j} + \frac{\partial \langle u_j \rangle^f}{\partial x_i} \right) - \frac{2}{3} \langle k \rangle^f \delta_{ij} $$

It is worth mentioning that $k$ and $\omega$ are volume averaged quantities and thus we have $k = \langle k \rangle$ and $\omega = \langle \omega \rangle$, which is not true in the DNS. Additionally, in the DANS/DEM model $\langle k \rangle$ and $\langle \omega \rangle$ were specially designed to assess the shear stress component $\langle u'v' \rangle$ (see eq. 15) through the Boussinesq relation since it is the main contribution entering the streamwise momentum equation for wall-parallel flows, $y$ being the wall-normal direction. Therefore, in the roughness sublayer, it is not expected that $\langle k \rangle$ be representative of the mean turbulent kinetic energy observed in the DNS which contributions come from the wake zones behind roughness elements and that involves high levels of diagonal stresses.

At last, as it is usually done in RANS turbulence models, a constant turbulent Prandtl number $Pr_t$ is used to relate the turbulent thermal conductivity $\kappa_t$ to the eddy viscosity $\nu_t$ and to deduce the turbulent heat flux from a Fourier relation $-\beta \langle u'_i \theta' \rangle = \kappa_t \frac{\partial \theta}{\partial u_i}$ with $\kappa_t = \frac{\nu_t C_p}{Pr_t}$ and where $\theta = T_w - T$ is the temperature difference with the wall. In this relation, the perfect gaz and incompressibility hypothesis were used to simplify the expression.

To compute the source term in the energy equation, Chedevergne (2021) used the relation:

$$ F_t = 4\kappa N u \frac{(1-\beta)}{d^2} \langle \theta \rangle^f = f_t \langle \theta \rangle^f \quad (11) $$

The drag coefficient and the local Stanton number $St = \frac{Nu}{RePr}$ are found to be related by the following relation $St/C_d \approx 0.09$ when analyzing the DNS database of Forooghi et al. (2018). Ultimately, this gives:

$$ f_t = \frac{2\pi}{11} f_d \quad (12) $$

In the case of 1D sheared flow, the dispersive stress and dispersive heat flux entering the momentum and energy equations (4) respectively, were closed with relations (Chedevergne, 2021) obtained from dimensional considerations on mixing lengths. The relation
for the dispersive shear stress reads:

$$ -\overline{\overline{\tilde{u}\tilde{v}}} = F_{disp} \nu_t \frac{\partial \sqrt{k}}{\partial y} \qquad (13) $$

with $F_{disp}$, a damping function, which restricts the extent of action of the dispersive stress above the roughness top $k_r$:

$$ F_{disp} = e^{-\left(\frac{y}{k_r} - 1\right)} \qquad (14) $$

The contributions of the dispersive shear stress on the momentum balance is generally low (Toussaint et al., 2020) and it is ignored in most DEM models Taylor et al. (1985); Kuwata et al. (2019); Stripf et al. (2008); Hanson et al. (2019). Nevertheless, the dispersive stresses can have a significant impact in some configurations, in particular at low Reynolds number, and it is thus interesting to be able to include them in a comprehensive modelling strategy for roughness effects.

The DANS/DEM model described above was extensively validated on both academic rough configurations (Chedevergne, 2021) for which DNS data were available. Therefore, only limited Reynolds numbers were tested, up to $Re_τ = 500$. Additional validation at high Reynolds numbers ($\delta^+ \approx 30000$) on two boundary layer configurations were also performed (Chedevergne, 2023) but only comparisons of velocity profiles in the inner region were performed. To further illustrate the performances of the DANS/DEM model on both the momentum and heat transfer upon rough walls, the experiments from Hosni et al. (1991) were simulated. It corresponds to boundary layers at moderate Reynolds numbers flowing over a flat plate covered with hemispheres (base diameter $d = 1.27$ mm) without pressure gradient. Velocity profiles, friction coefficient $C_f$ and Stanton number $St$ were carefully determined from the measurements. Two freestream velocity were tested, i.e. 11.8 m/s and 58.1 m/s. The DANS/DEM model was implemented in a boundary layer code from ONERA, named CLCET, already used by Aupoix (2015a,b), and which uses self-adaptative grids and several levels of grids to ensure convergence. Velocity profiles extracted at fixed Reynolds number values $R_{\delta_2}$, based on the momentum thickness $\delta_2$, are plotted against the experimental data for the two freestream velocity in the left part of figure 4. The corresponding temperature profiles are also depicted in the figure. The friction coefficient $C_f$ and the Stanton number $St$ are drawn in the right part of figure 4 with respect to $R_{\delta_2}$. Although measurements are not available in the roughness sublayers, computed profiles are in very good agreements with the experiments. Some discrepancies can be seen on the evolution of $C_f$ and $St$ but given the experimental uncertainty inherent to these measurements, the overall agreement is very satisfactory. These results shows the efficiency of the DANS/DEM model in a boundary layer context, at moderate to high Reynolds numbers (up to $\delta^+ \approx 10000$ and $R_{\delta_2} \approx 30000$). Beyond the good recovery of the roughness function $\Delta u^+$, the model is shown to behave nicely in the wake region, in a similar fashion to the original $k-\omega$ SST since the addition of the source terms in the DANS/DEM model only impact the roughness sublayer.

In all the following only fully developed channel flows will be considered to conform to
<figure><img src="image_5.png" /><figcaption>Fig. 4: (left) Mean velocity (blue) and temperature (orange) profiles computed with the DANS/DEM model (lines) and compared to the experimental results of Hosni et al. (1991) (symbols, △: 11.8 m/s and ◦: 58.1 m/s). Black dashed lines are the linear profile (⟨u⟩<sup>f+</sup> = y<sup>+</sup> and the logarithmic profile ⟨u⟩<sup>f+</sup> = 1/0.39 ln y<sup>+</sup> + 4.3 (Nagib and Chauhan, 2008) corresponding to a smooth wall configuration. (right) Friction coefficient C<sub>f</sub> (blue) and Stanton number St (orange) from the DANS/DEM model (lines) compared to measurements (symbols).</figcaption></figure>

the DNS configurations of Yang et al. (2022). A dedicated 1D code is used to solve the set of equations presented above but turned dimensionless. The code shares the same numerical ingredient as the boundary layer code CLICET. Using wall units, defined with the kinematic viscosity ν and friction velocity u<sub>τ</sub> and marked with a sign +, the momentum and energy equations in a channel configuration thus reduces to:

$$
\begin{align}
0 &= \frac{1}{Re_{\tau}} + \frac{1}{\beta} \frac{\partial}{\partial y^{+}} \left( \beta(1 + \nu_{t}^{+}) \frac{\partial \langle u \rangle^{f^{+}}}{\partial y^{+}} \right) + \frac{\langle u \rangle^{f^{+}}}{\beta} \frac{\partial^{2} \beta}{\partial y^{+}} - \frac{1}{\beta} \frac{\partial \overline{\tilde{u}\tilde{v}}^{+}}{\partial y^{+}} - F_{d_x}^{+} \nonumber \\
0 &= \frac{\langle u \rangle^{f^{+}}}{u_{b}^{+} Re_{\tau}} + \frac{1}{\beta} \frac{\partial}{\partial y^{+}} \left( \beta(1 + \kappa_{t}^{+}) \frac{\partial \langle \theta \rangle^{f^{+}}}{\partial y^{+}} \right) + \frac{\langle \theta \rangle^{f^{+}}}{\beta} \frac{\partial^{2} \beta}{\partial y^{+}} - \frac{1}{\beta} \frac{\partial \overline{\tilde{\theta}\tilde{v}}^{+}}{\partial y^{+}} - F_{t}^{+} \tag{15}
\end{align}
$$

The computations are driven by the pressure gradient which reads once turned dimen-
sionless $\frac{1}{Re_τ}$. The temperature difference in wall units is denoted $\theta^+$ and is given by
$\theta^+ = \frac{T-T_w}{T_τ}$, with $T_τ$ the friction temperature defined from the wall heat flux $\phi_w$,
i.e. $T_τ = \frac{-\phi_w}{\rho C_p u_τ}$. $C_p$ is the thermal capacity at constant pressure and in eq. (15) $u_b$
designates the bulk velocity in the channel.

**4 Results and discussions**

As a starting point, the twelve rough configurations were simulated by implementing
the DANS/DEM approach with the methodology outlined above. A good agreement
with the DNS reference was obtained for the mean velocity profiles. The results were
<figure><img src="image_6.png" /><figcaption>**Fig. 5:** Mean velocity (blue) and temperature (orange) profiles computed with the initial version of the DANS/DEM model (solid lines) and compared to DNS results (symbols).</figcaption></figure>

slightly improved by adjusting the impact of the effective slope $\varphi$ on the source term of the $\omega$-equation (10). The currently utilized DNS database (Yang et al., 2022) is particularly well designed to investigate the impact of $\varphi$ on turbulence over roughness. Therefore, the coefficient $c'_{\omega_2}$ was adjusted from 1 to 2. It is worth mentioning that this modification maintains consistency with previous findings (Chedevergne, 2021) while improving the agreement with the dataset analyzed in this study. These initial calculations also showed that it was necessary to adjust the ratio $St/C_d$, which drives the thermal closure in eq. (12). The ratio was changed to $St/C_d = 0.105$ yielding $f_t^+ = 0.65 f_d^+$. The turbulent Prandtl number value was fixed at $Pr_t = 0.85$ for all computations.

The results obtained for three representative configurations are shown in figure 5. The results show excellent agreement for the roughness configurations with positive or zero skewness values, i.e., the cases initiating with $G$ and $P$, in terms of both velocity profiles $\langle u \rangle^+$ and temperature profiles $\langle \theta \rangle^+$. However, for cases with negative $Sk$, significant underestimations of the effects of rough surfaces on the velocity and temperature profiles are observed. Lower values of the roughness functions $\Delta u^+$ and $\Delta \theta^+$ are predicted by the calculations compared to those in the DNS data. Potential additional effects present in these topographical configurations may not be adequately captured by the DANS/DEM model in its current form. Furthermore, as stated in section 3, the drag force (6) being aligned with $\langle u \rangle$ prevent the DANS/DEM model from reproducing separated zones in the near-wall region. This can be seen in the enlarged views of figure 5.

To elucidate the origin of the discrepancies observed in the configurations with negative skewness, we begin by examining the separation zones close to the reference wall, where the mean $\langle u \rangle$ profile shows negative values. Close-ups of the velocity profiles taken from the DNS data for all cases are depicted in figure 6. As can be observed, a reduction in $\varphi$ or an increase in skewness results in a corresponding decrease in both the size and intensity of the separation zones. The positions of the values of $y_{th}$, whose values in wall units vary slightly depending on the case considered because of $u_\tau$, whereas $y_{th}$ depends only on $Sk$ (see table 1), were plotted in figure 6 with thick grey lines covering the various values found. Although $y_{th}$ does not follow the
<figure><img src="image_7.png" /><figcaption>**Fig. 6:** Mean velocity (blue) and temperature (orange) profiles in the near wall region taken from the DNS of Yang et al. (2022). The vertical grey areas give the location of $y_{th}$ in the DANS/DEM model for the corresponding cases.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7:** Mean turbulent kinetic energy profiles in linear scale (linear) and in semi-logarithmic scale (right) taken from the DNS of Yang et al. (2022).</figcaption></figure>

$\varphi$ dependence observed on the near-wall profiles, the definition of $y_{th}$ seems coherent with respect to the length of the separated zones. In the DANS/DEM model, the threshold value $y_{th}$ delimits the zones that contribute to the momentum and energy balances from those that are not supposed to participate in transfer, such as the near-wall recirculation zones. This assumption holds well for configurations with positive or zero skewness values. However, in cases with negative skewness, these separation zones are suspected of enhancing turbulent mixing. To validate this assumption, we present the turbulent kinetic energy profiles as functions of both $y/H$ and $(y-\xi)^+$ in figure 7. For simplicity, we only consider cases X14 and X28, where X can be G, P, or N. Cases X18 and X24 are always bounded by the other two cases X14 and X28 with extreme $\varphi$ values, since varying the $p$ or $\lambda_c$ effectively modifies the effective slope $\varphi$. The left panel of the figure displays the shift and alteration of the $\langle k \rangle$ peak across the different cases, while the right panel highlights the loss of similarity of the profiles within the roughness sublayer.
The left plot of figure 7 clearly demonstrates that as the skewness decreases, the wall-normal position of the turbulent kinetic energy peak approaches the roughness height $k_r$. Additionally, it can be observed that increasing $\phi$ diminishes the intensity of the $\langle k \rangle$ peak and enhances turbulent mixing within the rough sublayer. The right panel of figure 7 presents the profiles plotted against $(y - \xi)^+$, emphasizing potential similarities across different cases, particularly within the roughness sublayer. It is evident that the overall shape of the $\langle k \rangle$ profiles is preserved for the positive and zero Sk cases. Nevertheless, this trend does not persist when Sk is in the negative region, as the $\langle k \rangle$ peak shifts to lower values of $(y - \xi)^+$. A more intense $\langle k \rangle$ zone also appears above this peak, suggesting a rise in turbulent mixing.

Figure 6 and 7 pointed out two phenomenons observed in the DNS and that are not accounted for in the DANS/DEM model. For separations zones in figure 6, negative Skewness values magnify the intensity and extent of the recirculation whereas in figure 7 the Skewness is shown to trigger the turbulence mixing enhancement observed in the rough sublayer when Sk becomes negative. In both cases, Sk has a determining role in the occurrence and intensity of these two phenomenons. However, the mechanisms at play in both situations being different, the modelling of both in the DANS/DEM model is addressed independently.

The existence of near-wall separation zones is directly related to the force exerted on the rough surface in the wall-normal direction. This non-zero force, $F_{dy}$, can be taken into account in our parallel approach ($\langle v \rangle = 0$) by applying a non-constant pressure gradient $\partial \langle p \rangle^f / \partial x$ under the roughness peak. In particular, negative values of $\partial \langle p \rangle^f / \partial x$ are expected to trigger flow separation. While the overall velocity profile across the channel height is little affected by the application of such a near-wall pressure gradient, the temperature profile may be more significantly perturbed. Heat transfer in both types of simulations (DNS and DANS/DEM) is indeed governed by an energy source $-\langle u \rangle^f / u_b \partial \langle p \rangle^f / \partial x$. This relation indicates that when a negative pressure gradient is applied where negative velocity values are obtained, the source term is positive, promoting the rise of the temperature profile. The influence of Sk on the pressure gradient within the roughness height is considered through the use of $y_{th}$, which was shown to be directly connected to Sk. A correction function $\alpha_p$ is introduced so that the pressure gradient is now $\alpha_p \partial \langle p \rangle^f / \partial x$. DNS results of MacDonald et al. (2018) containing extractions of pressure profiles were analyzed to define the shape of the function $\alpha_p$, which reads:

$$ \alpha_p = 1 - 4 \cos \left( 0.8 \pi \frac{y}{k_r} \right) \exp \left[ -2 \left( \frac{y - y_{th}}{y_{th}} \right)^2 \right] \quad (16) $$

In figure 8, the function $\alpha_p$ obtained for case G24 is plotted. Practically, in eq. 15, the term $1/Re_\tau$ is replaced by $\alpha_p/Re_\tau$. Since $1/h \int_0^h \alpha_p dy \approx 1$, the mean pressure gradient is preserved, as is the friction Reynolds number characterizing the simulation.

Accounting for the induced turbulent mixing due to separations is more challenging. A simple solution involves modifying the source term $F_\omega$ of the $\omega$-equation since it
<figure><img src="image_9.png" /><figcaption>Fig. 8: Function $\alpha_p$ used in the DANS/DEM model for case G24.</figcaption></figure>

controls the effects of roughness on turbulence. Although already substantial, the DNS database Yang et al. (2022) does not allow for accurately defining a correction to $F_\omega$ to include the impact of the separation zones in the DANS/DEM model. A careful review of the literature did not reveal any other data capable of enriching this database. Additionally, the enhancement of turbulent mixing due to separation zones is also expected to depend on the Reynolds number. In this respect, the database of MacDonald et al. (2018) provides some interesting information but does not enable the desired correction to be defined. To properly define a correction dedicated to this issue, specially designed additional DNS data are probably required. A better understanding of these phenomena and the possibility of identifying the most influential associated parameters are of major importance in modelling roughness effects, given that no model currently exists to tackle this specific issue. In the present study, a simple constant coefficient, noted $c_{\omega_s}$, is applied in front of the term $F_\omega$ for the $Nnm$ cases with negative $Sk$ to demonstrate the potential of a correction directly applied to the source term of the $\omega$ equation. The best results are obtained with $c_{\omega_s} = 0.5$ for configurations N14 and N18, $c_{\omega_s} = 0.65$ for case N24, and $c_{\omega_s} = 0.6$ for case N28.

To complete the analysis of the data and to estimate the relevance of the DANS/-DEM model, a focus is made on the evolution of the turbulent Prandtl number across the channel height. From their DNS database, MacDonald et al. (2018) already found a unique evolution of the turbulent Prandtl number $Pr_t$, independently of the roughness topography. As for smooth channel flow, $Pr_t$ is not constant in the channel height but satisfactory agreement are obtained using constant $Pr_t$ in eddy viscosity models. The turbulent Prandtl number for a given configuration is computed with the following relation:

$$ Pr_t = \frac{\overline{\langle u'v' \rangle} \frac{\partial \langle \theta \rangle^f}{\partial y}}{\overline{\langle v'\theta' \rangle} \frac{\partial \langle u \rangle^f}{\partial y}} \quad (17) $$

Figure 9 presents the evolution of $Pr_t$ for all the considered rough configurations. The figure highlights a distinct clustering of turbulent heat flux profiles with respect
<figure><img src="image_10.png" /><figcaption>**Fig. 9:** Turbulent Prandtl number $Pr_t$ (left), Reynolds shear stress – $\langle u't' \rangle$ (middle) and turbulent heat flux – $\langle v'\theta' \rangle$ (right) profiles extracted from the DNS of Yang et al. (2023). Blues lines are for Gnm configurations with $Sk = 0$, green lines for Pnm configurations with $Sk > 0$ and orange lines for Nnm cases with $Sk < 0$. Line styles correspond to *n* and *m* values and are indicated in each figure.</figcaption></figure>

to the rough configuration compared to that observed for the shear stress profiles. Two distinct clusterings are clearly visible in figure 9 for profiles – $\langle u't' \rangle$ and – $\langle v'\theta' \rangle$. Specifically, the three configurations P18, P28, and N28 exhibit turbulent heat flux – $\langle v'\theta' \rangle$ separated from the others. Since $Pr_t$ does not change between the different configurations, the consequence of these different clusterings between – $\langle u't' \rangle$ and – $\langle v'\theta' \rangle$ is a modification of the relative slopes of the velocity and temperature profiles in the logarithmic region from one case to another. From a modelling viewpoint, this suggests that $F_t$ should not be directly related to $F_{d_x}$. An additional contribution coming from $F_{d_y}$ may help recover the observed trends. A possible relationship could also exist between these observed behaviors for temperature and velocity profiles in the logarithmic zone and the existence of separation zones near the wall, which reinforces the need for more in-depth analysis of dedicated DNS data.

Figure 10 displays the final results obtained for the velocity and temperature profiles $\langle u \rangle^+$ and $\langle \theta \rangle^+$ for all twelve configurations tested. The special trend of the temperature profile of configuration N28 in the logarithmic zone cannot be reproduced with the present formulation of the DANS/DEM model. Similar but less pronounced effects are also visible in cases P18 and P28. The discrepancies between the reference DNS data and the DANS/DEM computations observed near the centerline of the channel are due to an overestimation of the eddy viscosity, as it was already discussed by Chedevergne (2021). The $k-\omega$ SST model, on which the DANS/DEM model relies, is known to reach too high values of eddy viscosity near the channel centerline, and this overestimation is even more pronounced when rough surfaces are considered. The DANS/DEM profiles tend to curve in this region, while DNS profiles exhibit a more straight-line behavior in a semi-logarithmic diagram. In a boundary layer configuration, such as that considered in figure 4, the wake region is well captured by the DANS/DEM model since the underlying $k-\omega$ SST model performs well in this region. This confirms that the observed discrepancies in channel configurations can be attributed to the performance of the $k-\omega$ SST model, and that the additional closures introduced in the DANS/DEM model to represent roughness effects are not responsible for the observed differences. With the exception of this particular point, the agreement between the DNS data and the DANS/DEM results is good for most
<figure><img src="image_11.png" /><figcaption>**Fig. 10:** Velocity (blue) and temperature (orange) profiles for the twelve configurations of Yang *et al* (2022). Solid lines are the DANS/DEM results and symbols are the DNS data.</figcaption></figure>
configurations, especially for cases with positive or zero skewness values, for which no
correction is applied to $F_\omega$ (i.e., $c_{\omega_s} = 1$), proving the predictive capabilities of the
DANS/DEM approach.

**5 Conclusion**

The DNS database established by Yang et al. (2022) shed particular light on the
influence of skewness on the distribution of roughness heights. These simulations high-
light the presence of recirculation zones in the very near-wall region, whose intensity
and extent increase with the skewness *Sk*. The surface-averaged streamwise velocity
profile therefore takes negative values in these zones, while it is also observed that
the amplitude of the streamwise force due to roughness, *F*<sub>d,x</sub>, starts growing. Recent
advances (Chedevergne and Forooghi, 2020; Chedevergne, 2023) in the discrete ele-
ment approach to tackle roughness effects in flows have led to the definition of the
DANS/DEM model (Chedevergne, 2021). Initial validation of this model has been
conclusive, but some aspects still need to be perfected. By construction, the DANS/-
DEM model is unable to reproduce the behaviors exhibited in the DNS of Yang et al.
(2022) in the near-wall region because of some strong hypotheses related to the drag
model (Chedevergne and Forooghi, 2020). However, these limitations prove to be of
little detriment to the prediction of roughness functions Δ*u*<sup>+</sup> and Δ*θ*<sup>+</sup>. For example,
the DANS/DEM model described in section 3 gives very satisfactory results for all the
cases in the DNS database of Yang et al. (2022) for which *Sk* is positive or zero. We
note that the model, as it stands, is unable to reproduce the near-wall recirculation
zones. These zones, which are assumed to have little influence on the overall balance
of roughness effects, are simply taken into account by means of a threshold *y*<sub>th</sub> in
the DANS/DEM model, below which the roughness parameters are considered con-
stant. It is shown that considering a non-constant pressure gradient profile ∂⟨*p*⟩<sup>f</sup>⁄∂*x
allows reproducing the near-wall separated regions and the corresponding increase in
the temperature profile. In terms of roughness functions Δ*u*<sup>+</sup> and Δ*θ*<sup>+</sup>, improvements
are small and confirm the weak impact of these recirculation zones on the velocity and
temperature profiles. On the other hand, for configurations with negative *Sk*, there
are substantial differences between the DANS and DNS calculations. When we look
at the turbulent kinetic energy profiles, we see that these are modified in the rough-
ness sublayer when *Sk* becomes negative. This indicates that the near-wall zones,
which were thought not to participate in the transfers, do in fact have an effect on
the turbulent mixing as *Sk* becomes negative. The ⟨*k*⟩ peak is shifted towards the
wall, and we observe the amplification of a second peak located above the roughness.
From a modelling point of view, one possible option to take this effect into account
is to modify the source term in the ω equation. Ad-hoc modifications to the constant
*c*<sub>ωs</sub> in equation (10) are proving effective. However, the limited extent of the current
DNS database makes it impossible to develop a rigorous correction. Furthermore, the
current DNS data are not sufficient to obtain a comprehensive understanding of the
mechanisms involved in configurations with negative skewness. Therefore, to further
enrich the model, particularly regarding the near-wall behaviors and the associated
heat transfer, multiple DNS computations on various rough surfaces, judiciously cho-
sen, are required. Another route that could be followed in parallel is to derive the
DANS/DEM approach on the basis of a Reynolds stress model, as initiated by Kuwata
et al. (2019). The DNS data in the roughness sublayers should then be carefully ana-
lyzed to develop closure relations for each of the Reynolds stress transport equations
and the dissipation equation, as has been done with the present model based on the
$k-\omega$ turbulence model.

**Acknowledgements.** The authors are grateful to M. MacDonald and N. Hutchins for sharing their DNS data (MacDonald et al., 2018). This work has been funded within the frame of the Clean Aviation Joint Undertaking, Multi-MW Hybrid-Electric Propulsion System for Regional Aircraft, being part of the Horizon Europe research and innovation funding programme of the European Commission. A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grants open access conditions.

Declarations

* Funding: This work has been funded within the frame of the Clean Aviation Joint Undertaking, Multi-MW Hybrid-Electric Propulsion System for Regional Aircraft, being part of the Horizon Europe research and innovation funding programme of the European Commission.

* Conflict of interest/Competing interests: The authors do not have any competing interests to declare that are relevant to the content of this article.

*   Ethics approval: This study does not involve any research with human participants and/or animals, so no ethical approval was required.

* Informed consent: This study does not involve any research with human participants animals, so no informed consent was required.

* Author contribution: F.C. performed the DANS/DEM developments and computations. J.Y.,A.S. and P.F. performed the DNS. F.C. and Y.J. performed detailed analysis of the DNS data in the rough sublayers. All authors reviewed the manuscript

*   Data availability Statement: DNS data that support the find-
    ings of this study have been deposited in the following address:
    https://publikationen.bibliothek.kit.edu/1000150842

References

Aupoix, B.: Improved heat transfer predictions on rough surfaces. International Journal of Heat and Fluid Flows **56**, 160–171 (2015) https://doi.org/10.1016/j.ijheatfluidflow.2015.07.007

Aupoix, B.: Roughness Corrections for the $k-\omega$ Shear Stress Transport Model: Status and Proposals. Journal of Fluids Engineering **137**, 021202 (2015) https://doi.org/10.1115/1.4028122
Aupoix, B.: Revisiting the discrete element method for predictions of flows over rough surfaces. Journal of Fluids Engineering, 31205 (2016) https://doi.org/10.1115/1.4031558

Busse, A., Ltzner, M., Sandham, N.D.: Direct numerical simulation of a turbulent flow over a rough surface based on a surface scan. Computers & Fluids 116, 129–147 (2015) https://doi.org/10.1016/compfluid.2015.04.008

Chung, D., Chan, L., MacDonald, M., Hutchins, N., Ooi, A.: A fast direct numerical simulation method for characterising hydraulic roughness. Journal of Fluid Mechanics 773, 418–431 (2015) https://doi.org/10.1017/jfm.2015.230

Chedevergne, F., Forooghi, P.: On the importance of the drag coefficient modelling in the double averaged Navier-Stokes equations for prediction of the roughness effects. Journal of Turbulence 21(8), 463–482 (2020) https://doi.org/10.1080/14685248.2020.1817465

Chedevergne, F.: A double-averaged Navier-Stokes $k-\omega$ turbulence model for wall flows over rough surfaces with heat transfer. Journal of Turbulence 22(11), 713–734 (2021) https://doi.org/10.1080/14685248.2021.1973014

Chedevergne, F.: Modeling rough walls from surface topography to double averaged Navier-Stokes computation. Journal of Turbulence 24(1-2), 36–56 (2023) https://doi.org/10.1080/14685248.2023.2173760

Christoph, G.H., Pletcher, R.H.: Predictions of rough-wall skin friction and heat transfer. AIAA Journal 21(4), 509–515 (1983)

Chevalier, M., Schlatter, P., Lundbladh, A., Henningson, D.: SIMSON–A pseudo-spectral solver for incompressible boundary layer flow. Tech. Rep. TRITA-MEK 2007:07, Royal Institute of Technology, Stockholm, Sweden, 1–100 (2007)

Finson, M.L., Clarke, A.S.: The Effect of Surface Roughness Character on Turbulent Reentry Heating. AIAA Paper 80-1459 15th Thermophysics Conference, Snowmass, Colorado (1980)

Finson, M.: A Reynolds stress model for boundary layer transition with application to rough surfaces. Interim scientific report, Physical Sciences Inc. Wakefield, Massachusetts (August 1975)

Finson, M.L.: A Model for Rough Wall Turbulent Heating and Skin Friction. AIAA Paper 82-0199 20th Aerospace Science Meeting, Orlando, Florida (1982)

Flack, K.A., Schultz, M.P.: Review of hydraulic roughness scales in the fully rough regime. Journal of Fluid Engineering 132, 041203–104120310 (2010)
Flack, K.A., Schultz, M.P., Barros, J.M.: Skin friction measurements of systematically-varied roughness: Probing the role of roughness amplitude and skewness. Flow, Turbulence and Combustion **104**(2-3), 317-329 (2019) https://doi.org/10.1007/s10494-019-00077-1

Forooghi, P., Stripf, M., Frohnapfel, B.: A systematic study of turbulent heat transfer over rough walls. International Journal of Heat and Mass Transfer **127**, 1157–1168 (2018) https://doi.org/10.1016/j.ijheatmasstransfer.2018.08.013

Forooghi, P., Stroh, A., Magagnato, F., Jakirlic, S., Frohnapfel, B.: Toward a universal roughness correlation. Journal of Fluids Engineering **139**(12), 121201 (2017) https://doi.org/10.1115/1.4037280

Goldstein, D., Handler, R., Sirovich, L.: Modeling a no-slip flow boundary with an external force field. Journal of Computational Physics **105**(2), 354–366 (1993)

Hosni, M.H., Coleman, H.W., Taylor, R.P.: Measurements and calculations of rough-wall heat transfer in the turbulent boundary layer. International Journal of Heat and Mass Transfer **34**(4/5), 1067–1082 (1991)

Hanson, D.R., Kinzel, M.P., McClain, S.T.: Validation of the discrete element roughness method for predicting heat transfer on rough surfaces. International Journal of Heat and Mass Transfer **136**, 1217–1232 (2019) https://doi.org/10.1016/j.ijheatmasstransfer.2019.03.062

Hellsten, A., Laine, S.: Extension of the $k-\omega$ shear-stress transport turbulence model for rough-wall flows. AIAA Journal **36**(9), 1728–1729 (1998)

Jackson, P.: On the displacement height in the logarithmic velocity profile. Journal of Fluid Mechanics **111**, 15–25 (1981)

Jelly, T.O., Busse, A.: Reynolds number dependence of reynolds and dispersive stresses in turbulent channel flow past irregular near-gaussian roughness. International Journal of Heat and Fluid Flow **80**, 108485 (2019) https://doi.org/10.1016/j.ijheatfluidflow.2019.108485

Knopp, T., Eisfeld, B., Calvo, J.B.: A new extension for $k-\omega$ turbulence models to account for wall roughness. International Journal of Heat and Fluid Flow **30**, 54–65 (2009)

Kuwata, Y., Kawaguchi, Y.: Direct numerical simulation of turbulence over resolved and modeled rough walls with irregularly distributed roughness. International Journal of Heat and Fluid Flow **77**, 1–18 (2019) https://doi.org/10.1016/j.ijheatfluidflow.2019.02.009

Kuwata, Y., Suga, K., Kawaguchi, Y.: An extension of the second moment closure model for turbulent flows over macro rough walls. International Journal of Heat
and Fluid Flow **77**, 186–201 (2019) https://doi.org/10.1016/j.ijheatfluidflow.2019.04.003

Lin, T.C., Bywater, R.J.: Turbulence models for high-speed, rough-wall boundary layers. AIAA Journal **20**(3), 325–333 (1982)

MacDonald, M., Chung, D., Hutchins, N., Chan, L., Ooi, A., Garcia-Mayoral, R.: The minimal-span, channel for rough-wall turbulent flows. journal of fluid mechanics **816**, 5–42 (2017) https://doi.org/10.1017/jfm.2017.69

McClain, S.T., Collins, S.P., Hodge, B.K., Bons, J.P.: The importance of the mean elevation in predicting skin friction for flow over closely packed surface roughness. Journal of Fluids Engineering **128**, 579–586 (2006)

Menter, F.: Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal **32**(8), 1598–1605 (1994)

MacDonald, M., Hutchins, N., Chung, D.: Roughness effects in turbulent forced convection. Journal of Fluid Mechanics **861**, 138–162 (2018) https://doi.org/10.1017/jfm.2018.900

Musker, A.J.: Universal roughness functions for naturally-occurring surfaces. Transactions of the Canadian Society of Mechanical Engineering **6**(1), 1–6 (1980-81)

Nagib, H., Chauhan, K.: Variations of von krmn coefficient in canonical flows. Physics of Fluids **20**(10) (2008) https://doi.org/10.1063/1.3006423

Nikuradse, J.: Laws of flows in rough pipes. Technical Memorandum 1292, NACA, Washington (April 1937)

Prez-R fols, F., Almqvist, A.: Generating randomly rough surfaces with given height probability distribution and power spectrum. Tribology International **131**, 591–604 (2019) https://doi.org/10.1016/j.triboint.2018.11.020

Schlichting, H.: Experimental investigation of the problem of surface roughness. Technical Memorandum 823, NACA, Washington (April 1937)

Sigal, A., Danberg, J.E.: New correlation of roughness density effect on the turbulent boundary layer. AIAA Journal **28**(3), 554–556 (1990)

Schultz, M., Flack, K.: Turbulent boundary layers on a systematically varied rough wall. Physics of Fluids **21**(015104) (2009) https://doi.org/10.1063/1.3059630

Stripf, M., Schulz, A., Bauer, H.-J.: Modeling of rough-wall boundary layer transition and heat transfer on turbine airfoils. Journal of Turbomachinery **130**(2) (2008) https://doi.org/10.1115/1.2750675
Thakkar, M., Busse, A., Sandham, N.D.: Direct numerical simulation of turbulent chaneel flow over a surrogate for nikuradse-type roughness. Journal of Fluid Mechanics **837**(R1) (2018)

Taylor, R., Coleman, H., Hodge, B.: Prediction of turbulent rough-wall skin friction using a discrete element approach. Journal of Fluids Engineering **107**, 251–257 (1985)

Toussaint, D., Chedevergne, F., Léon, O.: Analysis of the different sources of stress acting in fully rough turbulent flows over geometrical roughness elements. Physics of Fluids **32**(7), 075107 (2020) https://doi.org/10.1063/5.0010771

Žukauskas, A.: Heat transfer from tubes in crossflow. Advances in Heat Transfer, vol. 8, pp. 93–160. Elsevier (1972). https://doi.org/10.1016/S0065-2717(08)70038-8

van Rij, J.A., Belnap, B.J., Ligrani, P.M.: Analysis and experiments on three-dimensional, irregular surface roughness. Journal of Fluid Engineering **124**, 671–677 (2002)

Wu, S., Christensen, K.T., Pantano, C.: Modelling smooth- and transitionally rough-wall turbulent channel flow by leveraging innerouter interactions and principal component analysis. Journal of Fluid Mechanics **863**, 407453 (2019) https://doi.org/10.1017/jfm.2018.899

Whitaker, S.: Flows in porous media I: A theoretical derivation of Darcy’s law. Transport in Porous Media **1**, 3–25 (1986)

Wilcox, D.C.: Reassessment of the scale-determining equation for advanced turbulence models. AIAA Journal **26**(11), 1299–1310 (1988)

Yang, J., Stroh, A., Chung, D., Forooghi, P.: Direct numerical simulation-based characterization of pseudo-random roughness in minimal channels. Journal of Fluid Mechanics **941** (2022) https://doi.org/10.1017/jfm.2022.331

Yang, X.I., Sadique, J., Mittal, R., Meneveau, C.: Exponential roughness layer and analytical model for turbulent boundary layer flow over rectangular-prism roughness elements. Journal of fluid mechanics **789**, 127–165 (2016) https://doi.org/10.1017/jfm.2015.687

Yang, J., Velandia, J., Bansmer, S., Stroh, A., Forooghi, P.: A comparison of hydrodynamic and thermal properties of artificially generated against realistic rough surfaces. International Journal of Heat and Fluid Flow **99**, 109093 (2023) https://doi.org/10.1016/j.ijheatfluidflow.2022.109093