Let $n$ be a power of 2 satisfying $n \asymp C_0  \frac{\log \frac{1}{\delta}}{\log \frac{1}{\eps}}$, then $||(S_0^+)^n|_{V_{\ell}}||_{\mathrm{op}} = ||S_0^+|_{V_{\ell}}||^n$. Therefore it follows for $\delta \ll 2^{-\ell}$ that
$$||S_0^+|_{V_{\ell}}||_{\mathrm{op}} \leq D^{\frac{1}{n}} \max\{ e^{\frac{\sigma_1\ell}{n}} \delta^{\frac{\sigma_2}{n}}, \delta^{-\frac{\gamma}{n}} e^{- \frac{\sigma_3 \ell}{n}}\},$$
for $D, \sigma_1, \sigma_2,\sigma_3 > 0$ absolute constants. We choose $\delta = e^{-\max\{ 1, \frac{2 \sigma_1}{\sigma_2}  \}\ell}$ so that $\delta \ll 2^{- \ell}$ and $e^{\frac{\sigma_1\ell}{n}} \delta^{\frac{\sigma_2}{n}} \leq e^{-\frac{\sigma_1 \ell}{n}}$.  We furthermore set $\gamma = \frac{\sigma_3}{2\max\{ 1, 2 \frac{\sigma_1}{\sigma_2} \}}$ and therefore $\delta^{-\frac{\gamma}{n}} e^{- \frac{\sigma_3 \ell}{n}} = e^{- \frac{\sigma_3 \ell}{2n}}$. With these choices, $||S_0^+|_{V_{\ell}}||_{\mathrm{op}} \leq D^{\frac{1}{n}}e^{-O(1)\frac{\ell}{n}}$. In addition we make $\ell$ large enough in terms of $c_1$ and $c_2$ such that $\delta$ becomes small enough for Proposition <ref> to hold. To conclude, it holds by construction that $\frac{\ell}{n} \asymp_{c_1,c_2} \log \frac{1}{\varepsilon}$ and therefore $e^{-O(1)\frac{\ell}{n}} = \varepsilon^{O_{c_1,c_2}(1)}$ and similarly $D^{\frac{1}{n}} = \varepsilon^{-\frac{O_{c_1,c_2}(1)}{\ell}}$, so choosing $\ell$ additionally larger than a further constant depending on $c_1$ and $c_2$, the claim follows.

§.§ Proof of Theorem <ref>

Having established that $||S_0^{+}|_{V_{\ell}}||_{\mathrm{op}}$ is small for $L \geq L(c_1,c_2)$, we aim to convert this to an estimate that $||S_0^{+}|_{\bigoplus_{\ell \geq L}V_{\ell}}||_{\mathrm{op}}$ is also small. We use that the spaces $S_0^{+}V_{\ell}$ and $V_{\ell'}$ are almost orthogonal for $\ell \neq \ell'$ as shown in Lemma <ref>.

The Lie algebra of $K$ is denoted $\mathfrak{k}$ and we also write $\lambda_K$ for the Lie algebra representation induced by the regular representation $\lambda_K$ on $K$. Indeed, for a smooth function $\varphi$ on $K$ the function $(\lambda_K(X)\varphi)(k) = \lim_{t \to 0} \frac{1}{t}(\varphi(e^{-tX}k) - \varphi(k))$ with $X \in \mathfrak{k}$ and $k \in K$ is the directional derivative of $\varphi$ in the direction $-X$.

As in <cit.>, we use an argument based on partial integration to show that $S_0^{+}V_{\ell}$ and $V_{\ell'}$ are almost orthogonal. For a general manifold there is no suitable partial integration formula. However, for compact Lie groups we overcome this issue by exploiting that the Laplacian acts as a scalar on functions on $L^2(K)$ induced by the representation $\pi_{\gamma}$. Indeed, for a fixed orthonormal basis $X_1,\ldots , X_{\dim K}$ of $\mathfrak{k}$ recall that the Casimir element is defined as $\triangle = - \sum_{i} X_i \circ X_i$. We then use as replacement to partial integration that

\begin{equation}\label{PartialIntegration}
\langle \varphi_1, \lambda_K(\triangle)\varphi_2 \rangle = \sum_i \langle \lambda_K(-X_i) \varphi_1, \lambda_K(X_i) \varphi_2 \rangle.
\end{equation}

In order to give a suitable estimate for (<ref>), we first analyse $||\lambda_K(X) \varphi||_2$ for $X \in \mathfrak{k}$.

Let $\ell \geq 0$ and $\varepsilon > 0$. Then for $\varphi \in V_{\ell},g \in B_{\varepsilon}$ and $X \in \mathfrak{k}$ of unit norm,
$$||\lambda_K(X) \varphi||_2 \ll 2^{\ell} ||\varphi|| \quad\quad \text{and} \quad\quad ||\lambda_K(X) (\rho_0^+(g) \varphi)||_2 \ll (1 + O(\varepsilon^{O(1)}))2^{\ell}||\varphi||_2.$$

Without loss of generality we assume that $X \in \mathfrak{t}$. Fix $\gamma \in \overline{C} \cap I^{*}$. The eigenvalues of the operator $\pi_{\gamma}(e^{tX}) - \mathrm{Id}$ can be calculated as $e^{t\gamma'(X)} - 1$ for $\gamma'$ the various weights of the representation $\pi_{\gamma}$. Therefore the operator $\pi_{\gamma}(X) = \lim_{t \to 0} \frac{1}{t}(\pi_{\gamma}(e^{tX}) - \mathrm{Id})$ has eigenvalues $\gamma'(X)$. Let $v_1, \ldots , v_n$ be an orthonormal basis of eigenvectors of $\pi_{\gamma}(X)$. Then the functions $\psi(k) = \sqrt{d_{\gamma}} \langle \pi_{\gamma}(k)v_i, v_j \rangle$ for $k \in K$ satisfy $(\lambda_K(X)\psi)(k) = \sqrt{d_{\gamma}}\langle \pi_{\gamma}(k) v_i , \pi_{\gamma}(X)v_j \rangle = (\gamma'(X) \psi)(k)$. The first claim follows as $||\gamma'(X)|| \ll ||\gamma|| \leq 2^{\ell}$ and by decomposing the function $\varphi$ as a sum of functions of the form $\psi$.

For the second claim recall that $\rho_0^+(g)\varphi = \sqrt{\alpha_g'} \cdot (\varphi \circ \alpha_g)$ and therefore

\begin{equation}\label{DirectionalDerivativeMain}
\lambda_K(X)(\rho_0^+(g)\varphi) =  \left(\lambda_K(X) \sqrt{\alpha_g'} \right)\cdot (\varphi \circ \alpha_g) + \sqrt{\alpha_g'}\cdot \lambda_K(X)(\varphi \circ \alpha_g).
\end{equation}

To deal with the first term of (<ref>), since $\alpha_g'$ is a smooth polynomial perturbation of the identity, it follows that $||\lambda_K(X) \sqrt{\alpha_g'}||_{\infty} \leq (1 + O(\varepsilon^{O(1)}))$ and furthermore using integration by substitution, $||\varphi \circ \alpha_g||_2 \ll (1 + O(\varepsilon^{O(1)}))||\varphi||_2$. For the second term of (<ref>), we use the chain rule and the first step to conclude that $||\lambda_K(X)(\varphi \circ \alpha_g)||_2 \ll (1 + O(\varepsilon^{O(1)}))2^{\ell}||\varphi||_2$, concluding the lemma.

We now apply (<ref>) to prove the following lemma.

For $\varphi_{\ell_1} \in V_{\ell_1}$ and $\varphi_{\ell_2} \in V_{\ell_2}$ with $\ell_1 \neq \ell_2$ and $g \in B_{\varepsilon}$,
$$|\langle \rho_0^+(g)\varphi_{\ell_1}, \varphi_{\ell_2} \rangle| \ll (1 + O(\varepsilon^{O(1)})) 2^{-|\ell_1 - \ell_2|}||\varphi_{\ell_1}||_2 ||\varphi_{\ell_2}||_2. $$

Without loss of generality we assume that $\ell_2 > \ell_1$. Denote by $\psi\in V_{\ell_2}$ the function such that $\lambda_K(\triangle) \psi = \varphi_{\ell_2}$. Then by Lemma <ref>, $||\psi||_2 \ll 2^{-2\ell_2}||\varphi_{\ell_2}||_2$. Using then (<ref>) and Lemma <ref>,

\begin{align*}
|\langle \rho_0^+(g)\varphi_{\ell_1}, \varphi_{\ell_2} \rangle| &= |\langle \rho_0^+(g) \varphi_{\ell_1}, \lambda_K(\triangle) \psi \rangle | \\ &= \bigg|  \sum_{i} \langle \lambda_K(-X_i) \rho_0^+(g)\varphi_{\ell_1}, \lambda_K(X_i)\psi \rangle \bigg| \\ &\leq \sum_{i} ||\lambda_K(-X_i)(\rho_0^+(g)\varphi_{\ell_1})||\, ||\lambda_K(X_i)\psi|| \\
&\ll (1 + O(\varepsilon^{O(1)}))2^{\ell_1 + \ell_2} ||\varphi_{\ell_1}||_2 ||\psi||_2 \\
&\ll (1 + O(\varepsilon^{O(1)}))2^{\ell_1 - \ell_2} ||\varphi_{\ell_1}||_2 ||\varphi_{\ell_2}||_2.
\end{align*}

We conclude this section by proving Theorem <ref> by combining Proposition <ref> and Lemma <ref>.

(of Theorem <ref>)
By Proposition <ref>, there is $\varepsilon_0 = \varepsilon_0(c_1,c_2) > 0$ and $L = L(c_1,c_2) \in \Z_{\geq 1}$ such that $||S_0^+|_{V_{\ell}}||_{\mathrm{op}} \leq \varepsilon^{O_{c_1,c_2}(1)}$ for $\ell \geq L$. Let $\varphi \in \bigoplus_{\ell \geq L} V_{\ell}$ and let $N\geq 1$ to be determined later. Then

\begin{align*}
||S_0^+ \varphi||_2^2 &\leq \sum_{\ell, \ell' \geq L} |\langle S_0^+ \pi_{\ell} \varphi, S_0^+ \pi_{\ell'}\varphi \rangle| \\
&= \sum_{|\ell - \ell'|\leq N} |\langle S_0^+ \pi_{\ell} \varphi, S_0^+ \pi_{\ell'}\varphi \rangle| + \sum_{|\ell - \ell'| > N} |\langle S_0^+ \pi_{\ell} \varphi, S_0^+ \pi_{\ell'}\varphi \rangle|,
\end{align*}
where both of the sums are with $\ell, \ell' \geq L$. For the first of these two terms one uses the conclusion of Proposition <ref>,

$$\sum_{|\ell - \ell'|\leq N} ||S_0^+ \pi_{\ell} \varphi ||\,  ||S_0^+ \pi_{\ell'}\varphi ||\leq N \sum_{\ell \geq L}  ||S_0^+ \pi_{\ell}\varphi||_2^2 \leq N \varepsilon^{O_{c_1,c_2}(1)}||\varphi||_2^2.$$

Lemma <ref> is used to bound the second term:

\begin{align*}
\sum_{|\ell - \ell'| >  N} |\langle S_0^+ \pi_{\ell} \varphi, S_0^+ \pi_{\ell'}\varphi \rangle|  &\ll \sum_{|\ell - \ell'| >  N}  2^{-|\ell - \ell'|} ||\pi_{\ell} \varphi|| \, ||\pi_{\ell'}\varphi ||  \\ &\ll \sum_{|\ell - \ell'| >  N} 2^{-|\ell - \ell'|} ||\pi_{\ell}\varphi||_2^2 \\ &\ll 2^{-N} \sum_{\ell \geq L} ||\pi_{\ell} \varphi||_2^2  = 2^{-N}||\varphi||_2^2.
\end{align*}
Therefore it follows that $||S_0^+ \varphi||_2 \leq \sqrt{N \varepsilon^{O_{c_1,c_2}(1)} + 2^{-N}}||\varphi||_2$. Setting $N =  \log \frac{1}{\varepsilon}$ implies the claim of the theorem.

§.§ Smoothness of the Furstenberg Measure

In this section we prove Theorem <ref>, which we restate here for convenience of the reader.

(Theorem <ref>)
Let $G$ be a non-compact connected simple Lie group with finite center. Let $c_1, c_2 > 0$ and $m \in \Z_{\geq 1}$. Then there is $\eps_m = \eps_m(G,c_1,c_2) > 0$ depending on $G,c_1,c_2$ and $m$ such that every symmetric and $(c_1,c_2, \varepsilon)$-Diophantine probability measure $\mu$ with $\varepsilon \leq \varepsilon_m$ has absolutely continuous Furstenberg measure with density in $C^m(\Omega)$.

By Corollary <ref>, we know that the Furtstenberg measure is absolutely continuous if we choose $\eps_m$ small enough, i.e. there is $\FurstenbergDensity \in L^2(\Omega)$ such that $d\FurstenbergMeasure = \FurstenbergDensity d\Haarof{\Omega}$. In order to prove Theorem <ref>, we use the smoothness condition from Lemma <ref> for $\psi_F$. Indeed, for $P_{\ell}$ the projection from $L^2(K)$ to $V_{\ell}$, it suffices to show
$$||P_{\ell}\psi_F||_2 \leq 2^{-(s + 1)\ell}$$
for $s > m + \frac{1}{2} \dim K$ and $\ell$ large enough.

By the characterization of the Furstenberg measure, for any $n \geq 1$ it holds that $\FurstenbergMeasure = \mu^{*n} * \FurstenbergMeasure$ and therefore for $\varphi \in L^2(K)$,

\begin{align}
|\langle \FurstenbergDensity, \varphi \rangle| &= \bigg|\int \varphi \, d\FurstenbergMeasure\bigg| \nonumber \\ &= \bigg|\int\int \varphi(g.k) \, d\mu^{*n}(g)d\FurstenbergMeasure(k)\bigg| \nonumber \\ &\leq \bigg|\bigg|  \int \varphi \circ \alpha_g \, d\mu^{*n}(g)  \bigg|\bigg|_{\infty}. \label{ReductionPhi_nEstimate}
\end{align}

We thus study the $L^{\infty}$-norm of the function
$$\Phi_n = T_0^n \varphi =\int \varphi \circ \alpha_g \, d\mu^{*n}(g).$$
We will use  Corollary <ref> to give $L^2$-estimates of $\Phi_n$. In order to convert these estimates to an $L^{\infty}$-bound, we use Agmon's inequality (cf. <cit.> chapter 13), which we introduce for compact Lie groups.

(Agmon's Inequality for Compact Lie Groups).
Let $K$ be a compact Lie group. Then there is $t \in \Z_{\geq 2}$ depending on $K$ such that for any $\varphi \in C^{\infty}(K)$,
$$||\varphi||_{\infty} \ll ||\varphi||_2^{1/2} ||\varphi||_{H^t}^{1/2}.$$

For $M \in \R_{> 0}$ to be determined, we group together the contribution of the representations with $||\gamma|| \leq M$ and $||\gamma|| > 0$. Indeed, by (<ref>), for $k \in K$,

\begin{align*}
\varphi(k) &=  \sum_{\gamma \in \overline{C}\cap I^{*}} \sum_{i,j = 1}^{d_{\gamma}} d_{\gamma}^{1/2}a_{ij}^{\gamma} \chi_{ij}^{\gamma}(k) \\
&= \sum_{||\gamma|| \leq M} \sum_{i,j = 1}^{d_{\gamma}} d_{\gamma}^{1/2}a_{ij}^{\gamma} \chi_{ij}^{\gamma}(k)  + \sum_{||\gamma|| > M} \sum_{i,j = 1}^{d_{\gamma}} d_{\gamma}^{1/2}a_{ij}^{\gamma} \chi_{ij}^{\gamma}(k)  \\
&= \sum_{||\gamma|| \leq M} \sum_{i,j = 1}^{d_{\gamma}} d_{\gamma}^{1/2}a_{ij}^{\gamma} \chi_{ij}^{\gamma}(k)  + \sum_{||\gamma|| > M} \sum_{i,j = 1}^{d_{\gamma}} \lambda_{\gamma}^t\lambda_{\gamma}^{-t} d_{\gamma}^{1/2} a_{ij}^{\gamma} \chi_{ij}^{\gamma}(k),  \\
\end{align*}
where in the last line we multiplied the second term by $1 = \lambda_{\gamma}^t\lambda_{\gamma}^{-t}$ for some $t\in \Z_{\geq 0}$.  By Cauchy-Schwarz and using Lemma <ref>, the first term can be bounded by $||\varphi||_2 \sqrt{\sum_{||\gamma||\leq M, i,j} d_{\gamma}} \ll M^{C}||\varphi||_2$, where $C$ is a constant depending on $K$. For the second term, we choose $t$ large enough such that $\sqrt{\sum_{||\gamma|| > M,i,j} \lambda_{\gamma}^{-2t}d_{\gamma} } \ll M^{-C}.$ Again using Cauchy-Schwarz, the second term is bounded by $M^{-C}||\varphi||_{H^t}$. The claim is implied by setting $M = (\frac{||\varphi||_{H^t}}{||\varphi||_2})^{1/2C}$.

For $\varphi \in V_{\ell}$ set $\Phi_n = T_0^n \varphi$. Let $\gamma \in \overline{C}\cap I^{*}$ and $r \in \Z_{\geq 1}$. Then it holds for $\widehat{\Phi_n}(\gamma) = \pi_{\gamma}(\Phi_n)$,
$$||\widehat{\Phi_n}(\gamma)||_{\mathrm{op}} \ll_r 2^{O(1)\ell - r\ell}(1 + \varepsilon)^{O(1)nr} ||\gamma||^{O(1)r} ||\varphi||_2.$$

Let $v_1, \ldots , v_{d_{\gamma}}$ be an orthonormal basis of $\pi_{\gamma}$. Then

\begin{align}
||\widehat{\Phi_n}(\gamma)||_{\mathrm{op}} &\leq d_{\gamma} \sup_{1 \leq i \leq d_{\gamma}} ||\widehat{\Phi}(\gamma)v_i|| \nonumber \\ &\leq  d_{\gamma} \sup_{1 \leq i,j \leq d_{\gamma}} |\langle \widehat{\Phi}(\gamma) v_i, v_j \rangle| \nonumber \\ &= d_{\gamma} \sup_{1 \leq i,j \leq d_{\gamma}} |\langle \Phi_n, \chi^\gamma_{ij} \rangle| \nonumber \\
&\leq d_{\gamma} \sup_{\substack{g \in \supp{\mu^{*n}} \\ 1 \leq i,j \leq d_{\gamma}}} |\langle \varphi \circ \alpha_g, \chi^{\gamma}_{ij}\rangle|. \label{Phi_nBound}
\end{align}

Notice further that for $g \in B_{\varepsilon}$ and a further $\gamma' \in \overline{C} \cap I^{*}$ and $1 \leq i',j' \leq d_{\gamma}$,

\begin{align*}
|\langle \chi^{\gamma'}_{i'j'} \circ \alpha_g, \chi^{\gamma}_{ij} \rangle| &= \frac{\lambda_{\gamma}^r}{\lambda_{\gamma}^r}	|\langle \chi^{\gamma'}_{i'j'} \circ \alpha_g, \chi^{\gamma}_{ij} \rangle| \\
&= \frac{1}{\lambda_{\gamma}^r} |\langle \chi^{\gamma'}_{i'j'} \circ \alpha_g, \lambda_K(\triangle)^r\chi^{\gamma}_{ij}\rangle| \\
&\leq \frac{1}{\lambda_{\gamma}^r} \sum_{i_1, \ldots , i_r} |\langle \lambda(-X_{i_1})\cdots \lambda(-X_{i_r}) \chi_{i'j'}^{\gamma'} \circ \alpha_g , \lambda(X_{i_1})\cdots \lambda(X_{i_r})\chi^{\gamma}_{ij} \rangle| \\
&\ll_r \lambda_{\gamma}^{-r}(1 + \varepsilon)^{O(1)nr} ||\gamma'||^r\, ||\gamma||^r \\
&\ll_r (1 + \varepsilon)^{O(1)nr} ||\gamma'||^r\, ||\gamma||^{-r}
\end{align*}
where for the penultimate line one argues as in Lemma <ref> and in the last line we use Lemma <ref>. Similarly, it holds that $|\langle \chi^{\gamma'}_{i'j'} \circ \alpha_g, \chi^{\gamma}_{ij} \rangle| \ll_r (1 + \varepsilon)^{O(1)nr} ||\gamma'||^{-r}\, ||\gamma||^{r}$.
Then using the decomposition
$$\varphi = \sum_{2^{\ell-1}\leq ||\gamma'|| < 2^{\ell}}\sum_{i',j' = 1}^{d_{\gamma}} d_{\gamma'}^{1/2}a^{\gamma'}_{i'j'}\chi_{i'j'}^{\gamma'}$$
we conclude

\begin{align*}
|\langle \varphi \circ \alpha_g, \chi^{\gamma}_{ij}\rangle| &\leq \sum_{\gamma',i',j'} d_{\gamma'}^{1/2} |a_{i'j'}^{\gamma'}|\, |\langle \chi^{\gamma'}_{i'j'} \circ \alpha_g, \chi^{\gamma}_{ij} \rangle| \\
&\leq 2^{O(1)\ell}||\varphi||_2 \sup_{\gamma',i',j'} |\langle \chi^{\gamma'}_{i'j'} \circ \alpha_g, \chi^{\gamma}_{ij} \rangle| \\
&\ll_r 2^{O(1)\ell - r\ell}(1 + \varepsilon)^{O(1)nr} ||\gamma||^r ||\varphi||_2.
\end{align*}
This implies the claim by (<ref>) and using Lemma <ref>.

(of Theorem <ref>)
Let $\varphi \in V_{\ell}$ be of unit norm and write $\Phi_n = T_0^n \varphi$. It suffices to prove for $\varepsilon < \varepsilon_m$ and some some $n \geq 1$ that

\begin{equation}\label{AimDecayBound}
||\Phi_n||_{\infty} \leq 2^{-(s + 1)\ell},
\end{equation}

where $s$ is a constant depending on $G$ and $m$. Indeed, if (<ref>) holds, then by (<ref>),
$$||P_{\ell} \psi_F||_2 \ll 2^{O(1)\ell}2^{-(s + 1)\ell},$$
which satisfies the smoothness condition from Lemma <ref> for $s$ large enough depending on $G$ and $m$.

We will use Agmon's inequality to prove (<ref>). Notice first that for the fixed $t \in \Z_{\geq 2}$ from Lemma <ref>,

\begin{align*}
||\Phi_n||_{H^t}  &= ||\lambda_K(\triangle)^{t/2} \Phi_n||_2 \\ &\leq \sup_{g \in \mathrm{supp}\mu^{*n}} ||\lambda_K(\triangle)^{t/2}(\varphi \circ \alpha_g)||_{\infty} \\ &\leq ||\lambda_K(\triangle)^{t/2} \varphi||_{\infty}(1 + \varepsilon)^{O(1) n} \leq 2^{A\ell},
\end{align*}
for an absolte constant $A$ and where we choose $n = \frac{1}{10E_2\varepsilon}\ell$ for $E_2$ a fixed constant to be determined later.

We next bound $||\Phi_n||_2$. In order to do so, we decompose $\Phi_n$ into a low and high frequency part:

\begin{align*}
\Phi_n = \Phi_n^{(1)} + \Phi_n^{(2)} \quad\text{where}\quad \Phi_n^{(1)} = \sum_{||\gamma|| \leq L(c_1,c_1)} \sum_{i,j}^{d_{\gamma}} d_{\gamma}^{1/2}\widehat{(\Phi_n)}^{\gamma}_{ij} \chi_{ij}^{\gamma}.
\end{align*}

Then for $n \geq 1$, exploiting Corollary <ref>

\begin{align}
||\Phi_{n}||_2 &\leq \bigg|\bigg| \int \Phi_{n-1}^{(1)} \circ \alpha_g \, d\mu(g) \bigg|\bigg|_{\infty} + \bigg|\bigg| \int \Phi_{n-1}^{(2)} \circ \alpha_g \, d\mu(g) \bigg|\bigg|_{2} \nonumber  \\
&\leq ||\Phi_{n-1}^{(1)}||_{\infty} + \frac{1}{2} ||\Phi_{n-1}^{(2)}||_2 \label{Phi_nIteration}
\end{align}
Using Lemma <ref>, it follows for all $m \leq n$ and $r \geq 1$,
$$||\Phi_{m}^{(1)}||_{\infty} \ll_r 2^{O(1)\ell - r\ell}(1 + \varepsilon)^{O(1)nr}L(c_1,c_2)^{O(1)r}||\varphi||_2.$$
Iterating (<ref>), there are absolute constants $E_1,E_2,E_3 \geq 1$ such that
$$||\Phi_n||_2 \ll_r (n2^{E_1\ell - r\ell}(1 + \varepsilon)^{E_2nr}L(c_1,c_2)^{E_3r} + 2^{-n})||\varphi||_2.$$
By Lemma <ref>, it therefore follows that
$$||\Phi_n||_{\infty} \ll_r (n2^{(E_1 + A)\ell  - r\ell}(1 + \varepsilon)^{E_2nr}L(c_1,c_2)^{E_3r} + 2^{-n})||\varphi||_2.$$

Setting the parameters suitably, the proof is concluded. Indeed, choose for instance
$$r = 2(s + 1) + E_1 + A  +100$$
and $n = \frac{1}{10E_2\varepsilon}\ell$. For $s$ large enough and choosing $\varepsilon$ small enough in terms of $r$ and $s$ the claim (<ref>) holds for large $\ell$ (depending on $s$ and $\varepsilon$).


††thanks: Currently employed by Cambridge Flow Solutions Ltd., United Kingdom.

# Large Eddy Simulations of Flow over Additively Manufactured Surfaces: Impact
of Roughness and Skewness on Turbulent Heat Transfer

Himani Garg  Lund University, Department of Energy Sciences, SE-22100 Lund,
Sweden  Guillaume Sahut  Lund University, Department of Energy Sciences,
SE-22100 Lund, Sweden  Erika Tuneskog  Chalmers University of Technology ,
Department of Industrial and Materials Science , SE-41296 Gothenburg, Sweden
Karl-Johan Nogenmyr  Siemens Energy AB, SE-61231 Finspång, Sweden  Christer
Fureby<EMAIL_ADDRESS><EMAIL_ADDRESS>Lund University,
Department of Energy Sciences, SE-22100 Lund, Sweden

###### Abstract

## Abstract

Additive manufacturing creates surfaces with random roughness, impacting heat
transfer and pressure loss differently than traditional sand-grain roughness.
Further research is needed to understand these effects. We conducted high-
fidelity heat transfer simulations over three-dimensional additive
manufactured surfaces with varying roughness heights and skewness. Based on an
additive manufactured Inconel 939 sample from Siemens Energy AB, we created
six surfaces with different normalized roughness heights,
$R_{a}/D=0.001,0.006,0.012,0.015,0.020,$ and $0.028$, and a fixed skewness,
$s_{k}=0.424$. Each surface was also flipped to obtain negatively skewed
counterparts ($s_{k}=-0.424$). Simulations were conducted at a constant
Reynolds number of 8000 and with temperature treated as a passive scalar
(Prandtl number of 0.71). We analyzed temperature, velocity profiles and heat
fluxes to understand the impact of roughness height and skewness on heat and
momentum transfer. The inner-scaled mean temperature profiles are of larger
magnitude than the mean velocity profiles both inside and outside the
roughness layer. This means the temperature wall roughness function,
$\Delta\Theta^{+},$ differs from the momentum wall roughness function, $\Delta
U^{+}$. Surfaces with positive and negative skewness yielded different
estimates of equivalent sand-grain roughness for the same $R_{a}/D$ values,
suggesting a strong influence of slope and skewness on the relationship
between roughness function and equivalent sand-grain roughness. Analysis of
the heat and momentum transfer mechanisms indicated an increased effective
Prandtl number within the rough surface in which the momentum diffusivity is
larger than the corresponding thermal diffusivity due to the combined effects
of turbulence and dispersion. Results consistently indicated improved heat
transfer with increasing roughness height and positively skewed surfaces
performing better beyond a certain roughness threshold than negatively skewed
ones.

Large eddy simulations, Wall-bounded turbulence, Heat transfer, Additive
Manufacturing

##  I Introduction

Additive Manufacturing (AM) is being pursued for the demanding task of
manufacturing gas turbine hot gas path components. The challenge comes from
precisely maintaining the component at the highest temperature the material
can tolerate given the load they are exposed to, such that the benefits in
increased thermal efficiency are harvested. To develop such components,
precise control over heat transfer and pressure losses in the cooling system
is crucial. To this end, the studies detailed in this and previous papers have
been conducted [1; 2].

AM surfaces, particularly Powder Bed Fusion Laser Beam (PBF-LB), exhibit high
average surface roughness, typically falling within the range of
$R_{a}=5-25\mu$m. The high surface roughness of PBF-LB surfaces stems from the
powder-based manufacturing process, where a high-power laser selectively melts
and solidifies layers of powder material to fabricate components [3]. The
intricate interplay of melt pool dynamics, laser angle of incident, print
orientation, and other print process parameters contributes significantly to
the formation of surface roughness during the AM process [4]. The topography
of PBF-LB surfaces is characterized by the presence of partly sintered powder
particles or larger spatter particles from the melt pool on the surface,
accompanied by an underlying waviness originating from the melt pools. These
particles manifest as sharp peaks in surface roughness measurements, often
resulting in positively skewed surfaces. In addition, high peak density or
variations in melt pool characteristics can result in valley-dominated
surfaces.

Flow and heat transfer in channels have been studied for over a century. For
rough and smooth walls, reliable correlations are available for pressure
losses [5; 6; 7; 8]. For heat transfer, the situation is more complex. Heat
transfer in internal air-cooled systems is governed by Newton’s law of
cooling, $Q=A\alpha(T_{w}-T_{cool})$. Typically, in a gas turbine context, the
temperature difference between the cooling flow ($T_{cool}$) and the wall
temperature ($T_{w}$) is given. However, the cooling engineer can optimize the
design by adjusting the exposed surface area ($A$) and the heat transfer
coefficient ($\alpha$). The heat transfer coefficient is favored by "thin flow
passages," characterized by the hydraulic diameter, $D_{h}$, which can be seen
directly in the definition of the Nusselt number, $Nu={\alpha D_{h}/\lambda}$,
where $\lambda$ is the thermal conductivity of the fluid. Maximizing the
surface area while minimizing the hydraulic diameter can be achieved with a
cooling design that relies on in-wall cooling (or mini-) channels.

Smooth channels are thoroughly explored, and well-established correlations are
found (predominately Dittus-Bölter and Gnielinski [9; 10]). Regarding the heat
transfer dynamics over rough surfaces, various empirical studies by Nunner
[11], Dipprey & Sabersky [12], and Kays & Crawford [13] have investigated
predictive correlations for the Stanton number ($St$). Nunner [11] conducted
pioneering experiments on a surface roughened with two-dimensional transverse
ribs, providing an empirical expression for $St$ in terms of the Reynolds
number, Prandtl number, and the ratio of rough to smooth skin friction
coefficients. Additionally, Dipprey & Sabersky [12] proposed a semi-analytical
expression for $St$ based on the law of the wall similarity, incorporating the
Prandtl number, skin friction coefficient, and inner-scaled equivalent sand-
grain roughness, $k_{s}$. They addressed the challenge of incorporating the
effects of $k_{s}$ on $St$ by leveraging experimental data from a pipe with
sand-grain roughness. Subsequent studies have refined and modified these
expressions based on experimental and numerical investigations [14; 15].
Furthermore, due to challenges in obtaining accurate temperature fields within
the roughness sublayer, the underlying physics of how wall roughness
influences heat transfer remains ambiguous.

The effect of surface roughness has also been investigated numerically [16;
17; 18; 19; 20; 21; 22; 23; 24], with various methods. Several studies have
proposed roughness models to account for roughness in laminar flows.
Advancements in computer technology have facilitated Direct Numerical
Simulations (DNS) for studying turbulent heat transfer over resolved rough
surfaces [25; 26; 27; 28; 29]. MacDonald et al. [27] conducted DNS
specifically on turbulent heat transfer over sinusoidally rough surfaces,
indicating that the correlation function proposed by Dipprey & Sabersky [12]
accurately captured the influence of equivalent sand-grain roughness on $St$.
Analysis of instantaneous temperature fields revealed that dissimilarities
between heat and momentum transfer were primarily due to pressure drag effects
on the rough wall, enhancing momentum transfer but not heat transfer. Peeters
& Sandham [26] explored DNS for grit-blasted surfaces, affirming the efficacy
of Dipprey and Sabersky’s correlation [12] in predicting $St$ for such
surfaces. Notably, dissimilarities between heat and momentum transfers were
evident in the recirculation zone behind roughness elements, where effective
Prandtl numbers increased rapidly within the rough surface. Recent DNS studies
[27; 28] have further investigated the scaling behavior of the temperature
roughness function, $\Delta\Theta^{+}$, analogous to wall roughness function,
$\Delta U^{+}$, and observed notable deviations from the established $\Delta
U^{+}$ correlation, particularly at higher $k_{s}^{+}$ values.

In recent times, substantial efforts have been dedicated to investigating the
relationship between topological roughness parameters and the consequent
increase in friction factor ($c_{f}$). However, the surface roughness
generated by AM techniques exhibits spatial non-uniformity, influenced by
various factors such as profile curvature, layer thickness, laser power,
sample orientation, metallic composition, and particle size. As a result, AM
roughness differs considerably from regular, random, and artificial roughness.
To date, convective heat transfer in AM-made mini-channels remains
understudied, with limited experimental and numerical investigations [30; 31;
32; 33; 29]. Snyder et al. [31] demonstrated successful tailoring of surface
roughness through AM process parameters, enhancing the performance of a
generic micro-channel cooling design. More recent studies by Garg et al. [34;
29] have largely focused on systematically exploring the relationship between
AM roughness characteristics and $k_{s}$ to reveal their impact on turbulent
heat transfer using the similarity function, all-normal Reynolds stress, heat
flux, and probability density functions by using wall resolved Large Eddy
Simulations (LES). It has been explored how turbulence statistics is affected
by the presence of the roughness elements, particularly "peaks" dominant
surfaces, and the highly non-uniform heat transfer was shown to predominately
take place at the "peaks" of the wall roughness. In this paper, we continue
these studies and aim to show higher-order statistics of the roughness
characteristics that affect heat transfer by directly comparing "peak" and
"valleys" dominant surfaces. To the best of the authors’ knowledge, this study
represents the first instance of using LES to simulate turbulent heat transfer
in "peak" and "valley" dominated AM rough wall pipe flows with grid-conforming
three-dimensional roughness elements.

Figure 1: Illustration of the test specimen. (a) Naming of different levels.
(b) Test specimen in relation to gas flow, reacoater, and laser.

##  II Geometry creation

###  II.1 Experimental

The test specimen used for surface roughness measurements was manufactured
using PBF-LB on an EOS M400-4 system at Siemens Energy AB, Finspång. The
powder feedstock was Inconel 939, with a powder particle size distribution
ranging from 15 to 45 $\mu$m. Printing parameters were set according to
Siemens Energy standard and the inert gas was argon. The build direction was
set from bottom to top, as shown in Fig. 1. The specimen had a height of 30mm
and a maximum diameter of 15mm, featuring a total of 41 surfaces with
orientations at -60∘ (downskin), 0∘, 60∘ (upskin), and 90∘. Each orientation
(except 0∘) comprised a minimum of 8 surfaces. Surface measurements were
performed on each surface using both an optical high-resolution light
microscope, the Leica DM6 M, with focus variation, and an SEM, the Jeol JSM-
IT500. Post-processing, encompassing profile, and areal evaluation were
carried out using the Leica Map DCM software (Release: 7.4.8964).

Table 1: Summary of roughness parameters based on roughness height, $y$ (assuming mean line is at $y=0$), with continuous formulations, where $L_{x}$ is the surface length, and $p$, the probability density function. Parameter |  Description | Continuous formulation
---|---|---
$R_{a}$ |  Arithmetic average height [m] | $\frac{1}{L_{x}L_{z}}\int_{0}^{L_{x}}\int_{0}^{L_{z}}\left|y(x,z)\right|dzdx$
$R_{q}$ |  Root-mean-square roughness height [m] | $\sqrt{\frac{1}{L_{x}L_{z}}\int_{0}^{L_{x}}\int_{0}^{L_{z}}\left\\{y(x,z)\right\\}^{2}dzdx}$
$R_{z}$ |  Maximum height between the highest peak and the deepest valley of the profile [m] | $\left|\min y_{i}\right|+\left|\max y_{i}\right|$
$s_{k}$ |  Skewness [-] | $\frac{1}{R_{q}^{3}}\int_{-\infty}^{+\infty}y^{3}p(y)dy$
$k_{u}$ |  Kurtosis [-] | $\frac{1}{R_{q}^{4}}\int_{-\infty}^{+\infty}y^{4}p(y)dy$
Figure 2: Visualization of additive manufacturing rough surface height map,
$y(x,z)$, extracted from the measurement of 3D printed microchannels at
Siemens Energy AB. The surfaces in the top row have $s_{k}=0.424$; those in
the bottom row have $s_{k}=-0.424$ value. From left to right, the sampled
surface is the same, but the mean roughness height increases.

###  II.2 Rough surface characterization

Our goal in this study is to directly isolate the effects of roughness height
and skewness on turbulent heat transfer. To achieve this, we needed to tightly
control the roughness characteristics of the surfaces we investigated by using
a specific rectangular area from the downskin surface at position 5 (see Fig.
1) as a base. We then resized this base area to create six samples with
varying roughness heights, while keeping the underlying surface pattern
constant. The original surface was positively skewed, i.e., dominated by
peaks. To create samples with negative skewness, we simply flipped these
original surfaces upside down. In total, twelve rough surfaces were generated
based on the data provided by Siemens Energy AB, shown in Fig. 2. Flipping the
surfaces essentially transformed the peaks into valleys, all while maintaining
consistent roughness parameters. Table 1 summarizes the key characteristics
derived from the probability density function (PDF) of the surface height.
These characteristics include the average roughness height, $R_{a}$, the
maximum peak-to-valley height, $R_{z}$, the skewness factor, $s_{k}$, and the
kurtosis factor, $k_{u}$, all calculated using the height function, $y(x,z)$,
over sampling lengths $L_{x}$ and $L_{z}$ in the streamwise and spanwise
directions, respectively. We denote surfaces with positive skewness (peak-
dominated) as PS and those with negative skewness (valley-dominated) as NS for
brevity.

Table 2: Statistical quantities for the twelve rough surfaces considered. The descriptions and analytical expressions are given in Table 1. Here $D$ is the characteristic diameter of the rough pipe (defined later). Case |  $R_{a}/D$ [-] |  $R_{z}/D$ [-] |  $s_{k}$ [-] |  $k_{u}$ [-]
---|---|---|---|---
PS1 |  0.001 |  0.010 |  0.424 | 3.054
PS2 |  0.006 |  0.052 |  0.424 | 3.054
PS3 |  0.012 |  0.103 |  0.424 | 3.054
PS4 |  0.015 |  0.127 |  0.424 | 3.054
PS5 |  0.020 |  0.164 |  0.424 | 3.054
PS6 |  0.028 |  0.234 |  0.424 | 3.054
NS1 |  0.001 |  0.010 |  -0.424 | 3.054
NS2 |  0.006 |  0.053 |  -0.424 | 3.054
NS3 |  0.012 |  0.103 |  -0.424 | 3.054
NS4 |  0.015 |  0.127 |  -0.424 | 3.054
NS5 |  0.020 |  0.164 |  -0.424 | 3.054
NS6 |  0.028 |  0.234 |  -0.424 | 3.054

Figure 2(a) displays surfaces with a positive value of $s_{k}$, while Fig.
2(b) depicts surfaces with a negative value of $s_{k}$. Going from left to
right, the normalized roughness height ($R_{a}/D$) ranges from 0.001 to 0.028.
The roughness parameters, determined using the code developed by Garg et al.
[34; 29], are detailed in Table 2. Notably, the values of $k_{u}$ remain
constant across all rough surfaces examined in this study.

###  II.3 Geometries and meshes from experiment

A rectangular portion of the downskin surface at position 5, as shown in Fig.
1, was extracted and written in STL format. Figure 3 shows the transformation
of this rough plane to a rough pipe. The planar surface is mirrored along the
$x$-axis. This operation ensures that both sides of the surface parallel to
this axis are identical. Then a rotation around the $x$-axis is applied to all
points of the STL file. Along the closing line, all pairs of points merge
perfectly thanks to the previous mirroring step, as shown in Fig. 3.

Figure 3: A planar rough surface mirrored with respect to the black line and
then wrapped around the same direction to produce a rough pipe.

Practitioners can easily bend a planar surface using commercial and open-
source CAD software. However, merging the points along the closing line is
quite time-consuming and error-prone. Indeed, this operation relies on a
threshold distance to decide whether two points on the closing line should be
merged. The determination of this threshold value for a given STL file is
cumbersome.

To ease the process of creating rough cylinders, the authors have developed an
open-source Fortran code available on GitHub
111https://github.com/CoffeeDynamics/STLRoughPipes. The code runs on Unix
systems and takes as input parameters the STL file containing the planar rough
surface and a tolerance factor used to merge the points on the closing line.
Optionally, the user can specify a roughness factor to rescale the surface
roughness. Several values of tolerance factors may need to be tested in a
trial-and-error approach in order to find the most appropriate one. A
tolerance factor of $0.002$ was found accurate to properly merge the points on
the closing line for all the pipes considered in this work. The flexibility of
the code is a major asset for the creation of rough pipes from AM rough
surfaces. Moreover, the ability to rescale the surface roughness leads to
straightforward parametric studies of the impact of the roughness height and
skewness on turbulence and heat transfer. The code produces an STL file
containing the rough pipe, as shown in Fig. 3, and related statistical
quantities, as detailed in Table 1 (only the ones of interest for the present
work). This code has been used in our previous studies [34; 29].

Figure 4: Illustration of the mesh for PS6 rough surface.

We used snappyHexMesh, a meshing tool included in the OpenFOAM software (OF7),
to create meshes for the volumes within the cylinders. As an example, Fig. 4
shows the mesh generated for the PS6 rough surface. To precisely capture the
surface roughness, we refined the mesh near the walls. However, due to the
intricate surface details and the complexity of cells close to the walls,
obtaining exact values for non-dimensional wall distances like $x^{+}$,
$r^{+}$, and $z^{+}$ proved challenging. Therefore, we present estimated
average values for these quantities in our current meshes. Across all surfaces
at a Reynolds number ($Re_{b}$) of 8000, the average values of $x^{+}$,
$r^{+}$, and $z^{+}$ near the wall are all found to be below 2.35. In the
channel center, $r^{+}$ values range from 3.36 to 10.90 for PS1 to PS6 and NS1
to NS6 roughness configurations. The maximum values of skewness, aspect ratio,
and non-orthogonality for cells across PS1 to PS6 and NS1 to NS6 roughness are
observed to be between 0.72 and 1.59, 10.37 and 13.41, and $51.56^{\circ}$ and
$65.33^{\circ}$, respectively.

##  III Numerical methodology

###  III.1 Governing equations

This study involves conducting wall-resolved LES of thermal flow in pipes
roughened with AM structures. The flow is assumed to be fully incompressible,
Newtonian, and heated with a uniform wall heat flux $q_{w}$. The governing
equations are the spatially filtered mass, momentum, and energy conservation
equations, where the filtered quantities are denoted by overbars, with a
subgrid-scale (SGS) model introduced to represent small-scale unresolved
turbulence effects. The filtered LES equations are expressed as follows:

$\frac{\partial\overline{u}_{i}}{\partial x_{i}}=0,$ (1)
$\frac{\partial\overline{u}_{i}}{\partial t}+\frac{\partial}{\partial
x_{j}}(\overline{u}_{i}\overline{u}_{j})=-\frac{1}{\rho}\frac{\partial\overline{p}}{\partial
x_{i}}-\frac{\partial\tau_{ij}}{\partial x_{j}}+\frac{\partial}{\partial
x_{j}}\left[\nu\left(\frac{\partial\overline{u}_{i}}{\partial
x_{j}}\right)\right],$ (2) $\frac{\partial\overline{T}}{\partial
t}+\frac{\partial\left(\overline{T}\bar{u}_{i}\right)}{\partial
x_{i}}=-\frac{\partial q_{i}}{\partial x_{i}}+\frac{\partial}{\partial
x_{i}}\left[\frac{\nu}{Pr}\frac{\partial\overline{T}}{\partial x_{i}}\right],$
(3)

where $\overline{{u}}_{i}$ is the filtered velocity, $\overline{p}$ the
filtered pressure, $\nu$ the kinematic viscosity, and
$\tau_{ij}=\overline{u_{i}u_{j}}-\overline{u_{i}}\,\overline{u_{j}}$ the
subgrid-scale stress tensor, modeled using the Boussinesq hypothesis,
$\tau_{ij}=-2\rho\nu_{sgs}\bar{S}_{ij}$, where $\bar{S}_{ij}$ is the resolved
rate-of-strain tensor and $\nu_{sgs}$ is the subgrid viscosity. Additionally,
$T$ denotes the temperature, $\alpha=\nu/Pr$, the thermal diffusivity, $Pr$,
the Prandtl number, and
$q_{i}=\overline{u_{i}T}-\overline{u_{i}}\overline{T}$, the subgrid heat flux
vector. The subgrid viscosity is modeled using the Wall-Adapting Local Eddy-
viscosity (WALE) model [36], given by:

$\nu_{t}=(c_{w}\Delta)^{2}\frac{(S^{d}_{ij}S^{d}_{ij})^{3/2}}{(\overline{S}_{ij}\overline{S}_{ij})^{5/2}+(S^{d}_{ij}S^{d}_{ij})^{5/4}},$
(4)

where $c_{w}\simeq 0.325$ is a model parameter, $\Delta$ is the cell filter
size,
$S^{d}_{ij}=\frac{1}{2}(\overline{g}^{2}_{ij}+\overline{g}^{2}_{ji})-\frac{1}{3}\delta_{ij}\overline{g}^{2}_{ij}$
is the traceless symmetric part of the square of the velocity gradient tensor,
with $\overline{g}_{ij}=\partial\overline{u}_{i}/\partial x_{j}$ being the
filtered velocity gradient. The subgrid heat flux is modeled with a gradient
transport model [37], with $\alpha_{sgs}=\nu_{sgs}/Pr_{t}$, where $Pr_{t}$ is
the turbulent $Pr$ number.

With a constant heat flux boundary condition, the fluid’s bulk mean
temperature linearly increases with axial distance, complicating the
implementation of periodic boundary conditions for the temperature. To address
this, we adopt the approach suggested by Kasagi et al. [38], introducing a
transformed temperature $\Theta=T-T_{w}$ in the energy equation [29], where
$T_{w}$ is the temperature at the wall. The advantage of this method lies in
the periodicity of $\Theta$ in the streamwise direction [39; 40]. Under fully
developed flow conditions, the energy equation is modified to incorporate an
internal heat source as follows,

$\frac{\partial\overline{\Theta}}{\partial
t}+\frac{\partial\left(\overline{\Theta}\overline{u}_{i}\right)}{\partial
x_{i}}=-\frac{\partial q_{i}^{*}}{\partial x_{i}}+\frac{\partial}{\partial
x_{i}}\left[\frac{\nu}{Pr}\frac{\partial\overline{\Theta}}{\partial
x_{i}}\right]+\frac{q_{w}u_{\tau}}{\rho
c_{p}\nu_{w}}\frac{\overline{u}_{x}}{\left<\overline{u}_{x}\right>},$ (5)

where $u_{\tau}$ is the friction velocity and
$q_{i}^{*}=\overline{u_{i}\Theta}-\overline{u_{i}}\overline{\Theta}$, the
transformed subgrid heat flux vector.

###  III.2 Numerical methods

In this study, OF7 was used to simulate turbulent pipe flow with AM rough
walls. The LES equations were solved numerically using a second-order cell-
centered discretization scheme for convective and diffusive fluxes within the
finite volume method. Time stepping involved the implicit Adams-Bashforth
method [41], with a maximum Courant-Friedrichs-Lewy (CFL) number of 0.5 to
ensure numerical stability. The Pressure Implicit Splitting of Operators
(PISO) algorithm [42] was used to achieve pressure-velocity coupling, with
three corrector steps to minimize discretization errors. The resulting
pressure equation was solved using the Generalized Geometric-Algebraic
Multigrid (GAMG) method [43] with a Diagonal Incomplete Cholesky (DIC)
smoother [44]. For velocity solutions, the Preconditioned BiConjugate Gradient
(PBiCG) solver [45] was employed along with a Diagonal Incomplete Lower-Upper
(DILU) preconditioning method [44]. Additionally, to address non-orthogonality
in our roughness-conforming meshes, two inner loop correctors were
implemented, considering that the maximum non-orthogonality in all meshes
remained below $70^{\circ}$.

###  III.3 Initial and boundary conditions, dimensionless numbers

The simulations assume fully developed, turbulent, incompressible, and
Newtonian flow with constant properties ($\rho=1$ kg m-3 and $c_{p}=1005$ J
kg-1K-1). The effects of gravitational acceleration forces are neglected. The
value of $Pr$ was set to 0.71 and $Pr_{t}$ was set to 0.85. The simulations
explored the influence of two key factors: roughness height and surface
skewness. We simulated six different normalized roughness heights ($R_{a}/D$ =
0.001 to 0.028) on surfaces with two distinct skewness values ($s_{k}=-0.424$
and $0.424$), where $D$ is the characteristic diameter of the rough pipe
defined as $D=S/\left(\pi L_{x}\right)$ with $S$ being the total rough surface
area of the pipe and $L_{x}$ the pipe length. The $Re_{b}$ was kept constant
at $8000$ throughout all simulations. Here, $Re_{b}$ is defined as
$Re_{b}=U_{b}D/\nu$, where $U_{b}$ is the bulk velocity of the fluid. For easy
reference, we labeled the surfaces with positive skewness ($s_{k}=0.424$) as
PS followed by a number indicating the roughness height (e.g., PS1 for
$R_{a}/D=0.001$). Similarly, surfaces with negative skewness ($s_{k}=-0.424$)
were labeled NS followed by a number (e.g., NS6 for $R_{a}/D=0.028$).

In analyzing wall-bounded flows, the friction velocity scale, $u_{\tau}$, is
determined a posteriori using the pressure gradient, $\Delta p/L_{x}$,
calculated as $u_{\tau}=\sqrt{0.25D\Delta p/(\rho L_{x})}$. The friction
temperature, $\Theta_{\tau}$, is defined as $\Theta_{\tau}={q_{w}}/{\rho
c_{p}u_{\tau}}$. A key parameter, the turbulent friction $Re$ number,
$Re_{\tau}$, is defined as $Re_{\tau}=u_{\tau}D/(2\nu)$. In this study,
$Re_{\tau}$ ranges from 260 to 568 for NS1 to NS6 and 260 to 602 for PS1 to
PS6. Additionally, the roughness Reynolds number,
$k_{s}^{+}={k_{s}u_{\tau}}/{\nu}$, is considered, where $k_{s}$ is estimated
based on the roughness function, $\Delta U^{+}$ [34; 29]. This yields a range
of $k_{s}^{+}$ between 1 and 518 for negatively skewed surfaces and between 1
and 536 for positively skewed surfaces. Notably, in flows over roughened
walls, $u_{\tau}$ and the skin friction coefficient, $c_{f}$, reflect the
total wall drag, encompassing pressure and viscous drag instead of solely the
skin-friction drag.

The computational domain length, $L_{x}$, is set to eight times the pipe
diameter, i.e., $L_{x}=8D$, to mitigate periodicity effects in all considered
statistics [46; 47; 34; 29; 48]. The axial, radial, and azimuthal directions
are denoted by $x$, $r$, and $\theta$, respectively, with corresponding
velocity components $u_{x}$, $u_{r}$, and $u_{\theta}$. The effective wall-
normal distance is denoted by $r=r^{0}-d$, where $r^{0}$ is the average pipe
radius and $d$ is the zero-plane displacement [34; 29]. Periodic boundary
conditions are applied in the streamwise direction, while no-slip boundary
conditions are imposed at the wall for all velocity components. The zero-
Neumann boundary condition is used for pressure, with an additional forcing
term incorporated to address the pressure gradient effect within the context
of periodic boundary conditions. The boundary condition for the transformed
temperature is simply $\Theta=0$ at $r=0$ and $r=D$.

Statistical representation involves angular brackets, $\left<...\right>$, for
time-averaged quantities and primes, $\left(...\right)^{\prime}$, for
fluctuating quantities. The reported data correspond to a dataset extending to
$100t_{\text{ftt}}$, where $t_{\text{ftt}}=L_{x}/U_{b}$ is the flow-through
time.

##  IV Results and discussions

###  IV.1 Instantaneous flow

(a)

(b)

Figure 5: (a) Skin friction factor ($c_{f}$, colormap varies from pink
(minimum) to yellow (maximum)) and (b) heat transfer enhancement factor
($NuPr^{-1/3}$, colormap varies from blue (minimum) to red (maximum)), as a
function of mean roughness height, $R_{a}$, for (i) positively and (ii)
negatively skewed rough surfaces. (iii) Zooms are also included for clarity.

In the context of heat transfer over rough surfaces, the combined effect of
increasing surface roughness and skewness can result in a synergistic
enhancement of heat transfer. From a statistical sense, the surfaces with
positive skewness can be considered peak-dominated surfaces, and the ones with
negative skewness as valley-dominated surfaces. To explore the influence of
peaks and valleys alongside their height impact on the friction factor
$c_{f}=\Delta pD/(2L_{x}\rho U_{b}^{2})$ and Nusselt-Prandtl ratios, firstly,
the instantaneous flow visualizations of $c_{f}$ (absolute values) and
$NuPr^{-1/3}$ are shown in Fig. 5.

In Fig. 5(a) and 5(b), moving from left to right, the $s_{k}$ values are
fixed, but roughness height increases. Panels (i) and (ii) show the results
for positively and negatively skewed roughness, respectively. For small values
of $R_{a}/D$, the distribution of $c_{f}$ resembles the one obtained for
smooth pipe flows, and the impact of $s_{k}$ is negligible. With increasing
$R_{a}/D$ values, all surfaces show elevated $c_{f}$. Higher $c_{f}$ values at
peaks indicate that the roughness elements are inducing more substantial
energy losses in the flow. These energy losses are associated with the work
done by the fluid to overcome the increased resistance caused by the rough
surface. In geometries transitioning from peaks to cavities at a fixed
$R_{a}/D$ ratio, changes in the distribution of $c_{f}$ along the surface are
observed (see zooms included in panel (iii)). In peak geometries, such as
protrusions, higher $c_{f}$ values are typically found at the front head of
the peak due to flow separation and the formation of recirculation zones.
However, when the geometry is rotated to form a cavity, the flow behavior
alters, leading to flow reattachment downstream of the cavity. Consequently,
the downstream side of the cavity experiences increased turbulence levels and
higher momentum transfer near the surface, resulting in higher $c_{f}$ values
on the downstream side of the cavity. These variations in flow behavior and
boundary layer characteristics contribute to the observed changes in the
distribution of $c_{f}$ as well as in heat transfer along the surface of the
geometry, see Fig. 5(b). The heat transfer enhancement factor exhibited a
pronounced upward trend as the surface roughness increased. The presence of
roughness elements altered the near-wall flow patterns, leading to increased
turbulence and improved heat transfer rates. Furthermore, in peak geometries
where $c_{f}$ values are higher at the front head of the peak, the enhanced
turbulence near the surface promotes better mixing of the fluid, resulting in
higher convective heat transfer coefficients. Consequently, $NuPr^{-1/3}$
tends to be relatively higher on the upstream side of the peak compared to the
downstream side.

###  IV.2 Mean velocity and temperature

Figure 6 presents inner-scaled profiles of the streamwise mean velocity,
$\left<u_{x}^{+}\right>$, for positively and negatively skewed surfaces. The
non-dimensional wall distance, $r^{+}$, is measured from the plane at which
the total drag acts. This involves shifting the velocity profiles for rough
surfaces by removing the negative velocities, commonly referred to as zero-
plane displacement [49; 34; 29]. Additionally, the plot includes numerical
results of mean streamwise velocity profiles for LES of smooth pipe flow at
the same bulk $Re$ number ($Re_{b}=8000$) for comparison and standard log-
layer profiles for smooth walls (solid green line). In Figs. 6(a) and 6(b), it
is observed that increasing $k_{s}^{+}$ results in a downward shift of the
logarithmic profile of $\left<u_{x}^{+}\right>$, indicating enhanced momentum
transfer by the wall roughness. For small values of $R_{a}/D$, i.e., in the
hydraulically smooth regime as $k_{s}^{+}\approx 1$, the PS1 and NS1 profiles
are very close to the one for smooth surfaces, S0, as expected. The impact of
$s_{k}$ is little for PS2, PS3, NS2, and NS3 surfaces as the downward shift is
found to be roughly the same. For the larger value of $R_{a}/D>0.015$, the
downward shift in the logarithmic profile is found to be larger for positively
skewed surfaces compared to negatively skewed ones, signifying larger momentum
transfer by the wall roughness in the former ones.

In Figs. 7(a) and 7(b) mean temperature profiles, $\Theta^{+}$, are shown for
all surfaces with $s_{k}=0.424$ and $s_{k}=-0.424$, respectively. Similar to
$\left<u_{x}^{+}\right>$, an increasing downward shift is observed for
$\Theta^{+}$ with increasing wall roughness height in the region $r^{+}\geq
20$, while a slight upward shift is observed for $r^{+}\leq 15$, signifying
increased heat transfer close to the surface. For small roughness
($R_{a}/D=0.001$, $k_{s}^{+}=1$), the $\Theta^{+}$ profiles coincide with the
ones for S0 irrespective of the skewness, as expected. Furthermore, the
downward shift in the $\left<u_{x}^{+}\right>$ profile is consistently larger
than that of the $\Theta^{+}$ profile at the same $k_{s}^{+}$ value,
consistently with previous DNS studies [26; 50; 28; 51; 29]. The direct
comparison of $\Theta^{+}$ profiles in Figs. 7(a) and 7(b) for comparable
$R_{a}/D$ value shows that peak-dominated surfaces ($s_{k}=0.424$) have larger
impact on heat transfer enhancement than the valley-dominated ones
($s_{k}=-0.424$). In other terms, surfaces with the same value of $R_{a}/D$
but different $s_{k}$ lead to slightly different predictions of $k_{s}^{+}$,
and hence, larger roughness height and peak-dominated surfaces result in
higher turbulence levels and eventually higher heat transfer rates.

(a)

(b)

Figure 6: Mean velocity profile for (a) positively skewed and (b) negatively
skewed surfaces in the inner wall units.

(a)

(b)

Figure 7: Mean temperature profile for (a) positively skewed and (b)
negatively skewed surfaces in the inner wall units.

The magnitudes of the downward shift can be measured by using the roughness
function ($\Delta U^{+}$) and temperature difference function
($\Delta\Theta^{+}$). Figure 8(a) shows $\Delta U^{+}$ and Fig. 8(b) shows
$\Delta\Theta^{+}$ as a function of $k_{s}^{+}$, indicating the changes in the
mean profiles due to wall roughness. The results for $s_{k}=0.424$ are shown
in blue circles and the ones with $s_{k}=-0.424$, in red triangles. The chosen
scaling ensures that the roughness function collapses with Colebrook and
Nikuradse’s sand-grain data in the fully rough regime, as detailed by Garg et
al. [34; 48]. It is important to note that $k_{s}$ must be dynamically
determined for each specific rough surface and does not represent a simple
geometric length scale of the roughness. Consequently, in this study, the
value of $k_{s}^{+}$ is not known in advance and has been computed by
extracting $\Delta U^{+}$ from Fig. 6(a) for positively skewed surfaces and
Fig. 6(b) for negatively skewed surfaces and then using a Colebrook-type
roughness function, $\Delta U^{+}=2.44\log{\left(1+0.26k_{s}^{+}\right)}$ [5].

Figure 8(a) illustrates the progressive increase of $\Delta U^{+}$ with
$k_{s}^{+}$, eventually approaching the asymptotic limit for the fully rough
regime. The fully rough regime for the present rough surface is observed in
the range of $80<k_{s}^{+}<537$, which is somewhat comparable to sand-grain
roughness with $k_{s}^{+}\simeq 70$ [7]. The literature reveals similar
roughness functions obtained for sand-grain roughness at approximately the
same $Re$ number, even though the roughness heights significantly differ [52;
53; 54; 55; 34; 29]. Interestingly, for $k_{s}^{+}=1$, peak-dominated (PS1)
and valley-dominated (NS1) surfaces result in the same $\Delta U^{+}$ values,
as seen from $\left<u^{+}\right>$ profiles. With increasing $k_{s}^{+}$
values, we found that surfaces with $s_{k}>0$ result in a larger value of
$\Delta U^{+}$ compared to the one with $s_{k}<0$ even if the surfaces share
the same $R_{a}/D$ values. This suggests that the roughness height alone is
inadequate for scaling the momentum and heat transfer deficit resulting from
surface roughness.

Similarly, the temperature difference function ($\Delta\Theta^{+}$) is
computed and shown in Fig. 8(b) for $s_{k}>0$ (blue circles) and $s_{k}<0$
(red triangles). Globally, $\Delta\Theta^{+}$ initially increases with
$k_{s}^{+}$, but the rate of increase diminishes as $k_{s}^{+}$ becomes
larger. This behavior aligns reasonably well with Kays and Crawford’s
correlation [13] (dashed pink line):
$\Delta\Theta^{+}=\frac{Pr_{t}}{\kappa}\log(k_{s}^{+})-3.48\frac{Pr_{t}}{\kappa}-1.25{k_{s}^{+}}^{0.22}Pr^{0.44}+\beta(Pr),$
where $\kappa=0.4$ is the von Kármán constant and $\beta(Pr)=5.6$ represent
the log-law intercept. Additionally, we compare the results of
$\Delta\Theta^{+}$ for different rough surfaces, including three-dimensional
irregular roughness [28; 56] with positive and negative $s_{k}$ values and
three-dimensional AM roughness [29] with positive $s_{k}$ values. In the limit
of small $k_{s}^{+}\leq 10$, the value of $\Theta^{+}\approx 0$ indicates
negligible impact of roughness height and $s_{k}$ on heat transfer
enhancement. Interestingly, for $k_{s}^{+}>10$, we observe an increasing trend
of $\Delta\Theta^{+}$ against $k_{s}^{+}$, which is found to be consistent
with the reference data [56; 28; 29], despite variations in roughness type,
$R_{a}$, and $s_{k}$ values from previous studies. This suggests that
topological parameters, such as roughness type (AM or artificial), $s_{k}$,
and $ES$, have little impact on the $\Delta\Theta^{+}$ trend against
$k_{s}^{+}$, especially in the fully rough regime. However, differences in the
absolute value of $\Theta^{+}$ are noticeable for AM roughness of Garg et al.
[29], potentially due to variations in $Re_{b}$ and $s_{k}$ values used in
their simulations. Given that $\Delta U^{+}$ and $\Delta\Theta^{+}$ reflect
enhancements of the momentum and heat transfer, respectively, the observation
of $\Delta U^{+}>\Delta\Theta^{+}$ indicates that the wall roughness increases
the momentum transfer more than the heat transfer.

(a)

(b)

Figure 8: (a) Roughness function, $\Delta U^{+}$ and (b) temperature
difference, $\Delta\Theta^{+}$, as a function of equivalent sand-grain
roughness Reynolds number, $k_{s}^{+}$, for positively skewed (blue circles)
and negatively skewed (red triangles) surfaces.

###  IV.3 Effective turbulent Prandtl number

According to Kader [57] the log-region of the mean temperature profile
$\Theta$ in the absence of wall roughness can be described as

$\Theta^{+}(r^{+})=\frac{Pr_{t}}{\kappa}\log(r^{+})+\beta(Pr).$ (6)

In this expression,
$\beta(Pr)=\left(3.85Pr^{1/3}-1.3\right)^{2}+\left(Pr_{t}/\kappa\right)\log\left(Pr\right)$.
The literature typically suggests $Pr_{t}$ values in the range of $0.85$ to
$1$. Spalart and Strelets [58] proposed a method for determining an effective
Prandtl number within a recirculation bubble as
$Pr_{eff}=\nu_{eff}/\alpha_{eff}$, where
$\nu_{eff}=-\left<{u_{i}^{\prime}u_{j}^{\prime}}\right>\left<S_{ij}\right>$
and
$\alpha_{eff}=-\left<{u_{i}^{\prime}\Theta^{\prime}}\right>\left<\partial{\Theta}/\partial
x_{i}\right>/\left(\left<\partial{\Theta}/\partial x_{i}\right>\right)^{2}$.
Far from the wall, $Pr_{eff}\approx Pr_{t}$. Profiles of $Pr_{eff}$ are
displayed in Fig. 9 for positively and negatively skewed surfaces. $Pr_{eff}$
shows similarity in cases such as S0, PS1, PS2, NS1, and NS2, where
$k_{s}^{+}<10$. However, for larger $k_{s}^{+}$ values, $Pr_{eff}$
demonstrates a peak just below the plane $y/\delta=0$, which is equally
pronounced for both positively and negatively skewed surfaces. This peak
signifies a reduced effective thermal diffusivity, suggesting that regions
$r/R<0$ act as a thermal resistance. Moving away from the wall, $Pr_{eff}$
decreases with wall distance, reflecting the unmixedness of the scalar, as
discussed in more detail in [59; 60]. Additionally, far from the wall,
simulations with roughness show a slightly smaller $Pr_{eff}$ compared to
$Pr_{t}$ in smooth channel cases, especially for high $k_{s}^{+}$ with
$s_{k}=-0.424$. These results indicate that roughness indeed influences
$Pr_{t}$, with the influence of $s_{k}$ on these outcomes deemed insignificant
for the surfaces examined in this study. Moreover, the reason for the rise in
$Pr_{eff}$ within the rough wall can also be found in the contour maps of
Reynolds shear stress, $\left<u^{\prime}_{x}u_{r}^{\prime}\right>^{+}$ and
wall-normal component of heat flux,
$\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$ in Fig. 10. From Fig. 10, it
is evident that $\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$ and
$\left<u^{\prime}_{x}u_{r}^{\prime}\right>^{+}$ exhibit a similar trend but
there is a noticeable difference just behind the roughness crests, where
$\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$ is somewhat smaller than
$\left<u^{\prime}_{x}u_{r}^{\prime}\right>^{+}$. This is also reported by
Kuwata [28], attributing it to recirculation zones behind the roughness crest,
which had a detrimental effect on heat transfer. Moreover, the abrupt
fluctuations in $Pr_{eff}$ in the vicinity of the bottom ($r/R<0.01$) are
associated with the effective eddy diffusivity $\nu_{eff}$. In this region,
the local equilibrium state of turbulence breaks down; enhanced pressure
diffusion acts as an energy source near the bottom of the rough surface,
leading to an enhancement of the wall-normal Reynolds stress and the Reynolds
shear stress [61].

(a)

(b)

Figure 9: Effective Prandtl number plotted in outer units for (a) positively
skewed surfaces and (b) negatively skewed surfaces. Results for the benchmark
case, i.e., smooth surface, are also included for reference. The dashed red
line shows $Pr_{t}=0.85$ fixed in the simulations. Figure 10: Contour maps of
wall-normal heat flux component,
$\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$, and Reynolds shear stress,
$\left<u^{\prime}_{x}u^{\prime}_{r}\right>^{+}$, in a $x-y$ plane for
$R_{a}/D=0.028$ and $s_{k}=0.424$. Red dashed lines indicate the region where
the difference between $\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$ and
$\left<u^{\prime}_{x}u^{\prime}_{r}\right>^{+}$ is apparent.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11: The wall-normal Reynolds shear stress,
$\left<u^{\prime}_{x}u^{\prime}_{r}\right>$, streamwise velocity fluctuations,
$\left<u^{\prime}_{x}u^{\prime}_{x}\right>$, normalized by $u_{\tau}^{2}$ and
streamwise heat flux, $\left<u^{\prime}_{x}\Theta^{\prime}\right>$ and wall-
normal heat flux, $\left<u^{\prime}_{r}\Theta^{\prime}\right>$, normalized by
$u_{\tau}T_{\tau}$. Solid and dashed lines represent negatively skewed and
positively skewed surface results, respectively.

###  IV.4 Impact of skewness and roughness height on Reynolds stresses and
heat fluxes

The temperature fluctuations, $\Theta^{\prime}$, and velocity fluctuations,
$u_{x}^{\prime}$ and $u_{r}^{\prime}$, give rise to significant mean Reynolds
stresses and heat fluxes. Recent studies by Garg et al. [29] have shown the
decreasing correlation between streamwise velocity fluctuations,
$\left<u^{\prime}_{x}u^{\prime}_{x}\right>$ and streamwise component of heat
flux, $\left<u^{\prime}_{x}\Theta^{\prime}\right>$, with increasing roughness
height. However, the influence of roughness topologies, such as peak-dominated
roughness with $s_{k}>0$ and valley-dominated ones with $s_{k}<0$, is still
unclear. Thus, the effect of $s_{k}$ on the correlation between
$\left<u^{\prime}_{x}\Theta^{\prime}\right>$ and
$\left<u^{\prime}_{x}u^{\prime}_{x}\right>$ and between
$\left<u^{\prime}_{r}\Theta^{\prime}\right>$ and
$\left<u^{\prime}_{r}u^{\prime}_{x}\right>$ with increasing roughness height
is investigated in this section.

In Fig. 11, we compare $\left<u^{\prime}_{x}u^{\prime}_{x}\right>^{+}$ with
$\left<u^{\prime}_{x}\Theta^{\prime}\right>^{+}$ and
$\left<u^{\prime}_{x}u^{\prime}_{r}\right>^{+}$ with
$\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$, for the range of
$R_{a}/D=0.001-0.028$, yielding a range of $k_{s}^{+}=1-536$. All solid and
dashed lines represent negatively skewed and positively skewed surfaces,
respectively. All cases show that irrespective of the roughness height and
skewness values, the general trend between
$\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$ and
$\left<u^{\prime}_{x}u^{\prime}_{r}\right>^{+}$ is practically identical, with
a slight difference in the magnitude,
$\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$ being slightly smaller than
$\left<u^{\prime}_{x}u^{\prime}_{r}\right>^{+}$. This is consistent with
$Pr_{eff}$ being greater than one in the vicinity of the wall and the results
shown in Fig. 10, indicating the recirculation zone behind the roughness wall
acting as a thermal resistance. With an increase in roughness height, we
observe that the peak value of these quantities shifts outwards after a
critical value of $R_{a}/D>0.006$, i.e., outside the hydraulically smooth
regime. For larger $R_{a}/D$ values, the impact of $s_{k}$ starts to become
visible for $R_{a}/D>0.012$ as the peak location starts moving outward for the
positively skewed surfaces compared to the negatively skewed ones, which means
positively skewed surfaces generate larger friction, consistently with our
earlier results. This outward shift can also be explained by calculating the
modified diameter, $D^{mod}_{0}=D_{0}-d$, where $D_{0}$ is the smooth pipe
diameter. The estimation of normalized $D^{mod}_{0}$ is shown in Fig. 12. Here
we can see that the “effective” diameter is about 10% larger for the negative
skewness than the positive one (at $Ra/D_{0}=0.032$), having a profound
influence on the overall flow. Effectively, the positive one squeezes the flow
towards the center, which causes a larger outward shift of the peak for
positively skewed surfaces compared to negative ones. A slight difference
between $\left<u^{\prime}_{x}u^{\prime}_{x}\right>$ and
$\left<u^{\prime}_{x}\Theta^{\prime}\right>$ is seen in the magnitude. The
general trend is very similar for these two quantities as well. Interestingly,
the difference between $\left<u^{\prime}_{x}\Theta^{\prime}\right>$ and
$\left<u^{\prime}_{x}u^{\prime}_{x}\right>$ is a non-monotonic function of
$R_{a}/D$.

Figure 12: Estimation of modified rough surface diameter due to negative
velocities for positively and negatively skewed surfaces.

###  IV.5 Impact of skewness and roughness height on effective heat transfer

Figure 13(a) presents the friction factor as a function of $R_{a}/D$ for
various rough surfaces with $s_{k}=0.424$ and $s_{k}=-0.424$ values,
normalized with respect to the numerical results obtained for smooth pipe flow
using the Moody chart. The plots cover a $R_{a}/D$ number range of 0.001 to
0.028. Due to the lack of experimental validation for rough surfaces, we
employed LES with the WALE model for smooth pipe flow at the same mesh
resolution as rough pipes and matched Reynolds numbers. Previous studies by
Garg et al. [48; 29] utilized LES where the WALE model results for smooth pipe
flow were used as a validation tool, showing excellent agreement with
theoretical predictions from the Moody chart, with a deviation of less than 3%
(refer to Fig. 10(a) dashed black and solid green lines of Garg et al. [29]).
This validation process instills confidence in the accuracy and robustness of
our numerical approach for predicting the friction factor behavior in smooth
pipe flow. By extension, it reinforces the reliability of our numerical
results for turbulent flow over rough surfaces, given the utilization of the
same numerical approach.

In Fig. 13(a), the normalized $c_{f}/c^{0}_{f}$ values indicate a consistent
upward trend with increasing roughness height ($R_{a}/D$ increases from 0.001
to 0.028), regardless of $s_{k}$ values. For PS1 and NS1 with $R_{a}/D=0.001$
and $k_{s}^{+}\approx 1$, representing hydraulically smooth roughness, $c_{f}$
is negligibly impacted, leading to $c_{f}/c_{f}^{0}$ approaching 1. In the
transition regime of PS2 to PS4 and NS2 to NS4, where $k_{s}^{+}$ ranges from
9 to 50, roughness elements disrupt the laminar flow, inducing additional drag
and increasing resistance compared to smooth pipe flow, thus explaining the
observed increasing behavior of $c_{f}/c_{f}^{0}$. The influence of $s_{k}$
remains insignificant up to $R_{a}/D=0.03$. However, as $R_{a}/D$ surpasses
this threshold, rough surfaces PS5, PS6, NS5, and NS6 yield $k_{s}^{+}$ values
between $100$ and $536$, indicating a fully rough regime where intensified
turbulence significantly increases $c_{f}$. The impact of $s_{k}$ becomes
pronounced, with positively skewed surfaces displaying higher $c_{f}$ than
negatively skewed ones, indicating a greater imbalance between turbulence
effects and roughness-induced drag in peak-dominated surfaces. Surprisingly,
up to considerable $R_{a}/D$ values, the impact of $s_{k}$ on $c_{f}$ is
negligible.

In Fig. 13(b), the normalized $Nu$ number is shown as a function of $R_{a}/D$
for various rough surfaces, normalized against the Nusselt number obtained for
smooth pipe flow, $Nu^{0}$, using the Dittus-Boelter correlation
$Nu^{0}=0.023Re_{b}^{0.8}Pr^{0.4}$. LES with the WALE model results for smooth
pipe flow were utilized as a validation tool, exhibiting excellent agreement
with theoretical predictions, with a deviation of less than 4% (refer to Fig.
10(b) dashed black and solid green lines in Garg et al. [29]). Similar to the
prediction of $c_{f}$, the observed behavior indicates potent heat transfer
enhancement within the considered range of $R_{a}/D$ due to roughness-induced
disturbances. In the hydraulically smooth regime, the impact of protrusions
and valleys is negligible on the $Nu$ number, leading to $Nu/Nu^{0}$
approaching 1, regardless of $s_{k}$. As the transition roughness regime is
entered, $Nu$ values begin to increase, resulting in potent heat transfer
enhancement, with the impact of $s_{k}$ in this regime found to be negligible.
However, in the fully rough regime, $Nu$ values exhibit a consistent upward
trend, with the influence of $s_{k}$ becoming apparent. The results
demonstrate that surfaces dominated by peaks ($s_{k}>0$) result in higher heat
transfer enhancement compared to those dominated by valleys ($s_{k}<0$).

To quantify the heat transfer enhancement effectiveness for the current
roughness and assess the influence of surface topology, we calculated the
thermal performance factor ($TPF$), representing the ratio of the relative
change in heat transfer rate to the change in friction factor:
$TPF=\left(Nu/Nu^{0}\right)/\left(c_{f}/c_{f}^{0}\right)^{1/3}$, illustrated
in Fig. 13(c). Three distinct regimes emerge. Firstly, the hydraulically
smooth regime, where all rough surfaces perform equivalently to smooth ones.
Secondly, the transitionally rough regime, where $TPF$ notably increases,
remains between 1 and 1.5, signifying superior performance of all roughness
configurations compared to smooth surfaces in transitional flow. The influence
of $s_{k}$ becomes noticeable for $R_{a}/D=0.015$. Lastly, in the fully rough
regime, the impact of roughness height and skewness becomes markedly
pronounced. Surfaces with positive $s_{k}$ values reveal a larger and more
monotonic increase in $TPF$ compared to those with negative $s_{k}$ values.
This highlights the critical importance of surface topology, where positively
skewed surfaces outperform negatively skewed ones for similar roughness
heights.

(a)

(b)

(c)

Figure 13: Estimation of (a) skin friction factor, $c_{f}$, (b) Nusselt
number, $Nu$, and (c) thermal performance factor, TPF, as a function of
$R_{a}/D$. All values are normalized by the benchmark case, i.e., the smooth
pipe flow. The dashed black line represents the expected smooth pipe flow
results.

###  IV.6 Effect of surface anisotropy on the anisotropy of the Reynolds
stress tensor

The net anisotropy of the Reynolds stresses is commonly quantified using the
second, $II_{b}$, and third, $III_{b}$, invariants of the normalized
anisotropy tensor, $b_{ij}$, given by [62]

$b_{ij}=\frac{\left<\overline{u^{\prime}_{i}u^{\prime}_{j}}\right>}{\left<\overline{u^{\prime}_{k}u^{\prime}_{k}}\right>}-\frac{1}{3}\delta_{ij}.$
(7)

The state of anisotropy can then be characterized with the two variables
$\eta$ and $\xi$ defined as

$\eta^{2}=-\frac{1}{3}II_{b}$ (8)

and

$\xi^{3}=-\frac{1}{2}III_{b}.$ (9)

The Reynolds stress tensor’s feasible states are confined within a triangular
region in the ($\xi$, $\eta$) plane, known as the Lumley triangle. Distinct
turbulence scenarios can be distinguished by examining the two invariants of
$b_{ij}$ at the Lumley triangle’s theoretical extremes. Garg et al. [48]
provided a comprehensive analysis of these states. Figures 14 and 15 depict
Lumley triangle plots for rough surfaces with specific values of $s_{k}$
(0.424 and -0.424, respectively), along with a reference plot for smooth pipe
flow (Figs. 14(a) and 15(a)). Data are extracted from a $y-z$ plane slice at a
consistent location for all surfaces, as illustrated in the inset of Fig.
14(a). Utilizing PDF post-processing, cell-centered values of $\xi$ and $\eta$
are visualized with a colormap to correspond with the wall-normal distance. In
the smooth-wall turbulent pipe flow reference case, near-wall flow closely
resembles a two-component state along the upper boundary of the Lumley
triangle, progressing towards the one-component state at the triangle’s upper-
right apex. Maximum anisotropy occurs at approximately $r^{+}\approx 8$,
beyond which anisotropy diminishes. For $r^{+}>8$, the ($\xi$, $\eta$) curve
aligns closely with the right boundary of the Lumley triangle, indicating
proximity to an axisymmetric, rod-like state, but not fully reaching the
maximum isotropy state at the triangle’s base summit as $r^{+}$ increases. For
smaller $R_{a}/D$ values, the behavior of surfaces PS1 and NS1 closely
resembles that of smooth surface flows, irrespective of $s_{k}$ values.

For rough surfaces with $R_{a}/D>0.01$, regardless of $s_{k}$ values,
turbulent states occupy various positions within the Lumley triangle, except
for the plain strain condition marked by the dashed line (see Figs. 14(b-f)
and 15(b-f)). In the deepest valleys of the surfaces, the Reynolds stress
anisotropy tensor tends toward a strongly anisotropic, one-component state. As
$r^{+}$ increases, the results diverge from the trajectory observed for
smooth-wall pipe flow on the $(\xi,\eta)$-map. Across all rough surfaces (PS1
to PS6 and NS1 to NS6), the flow converges towards the left side of the Lumley
triangle, attaining an axisymmetric disk-like state at the roughness mean
plane. Here, the streamwise and azimuthal Reynolds stresses exhibit comparable
magnitudes. Anisotropy tends to center around the axisymmetric expansion and
the two-component limit, indicating that one stress component predominates or
that two components are similar in magnitude. Such axisymmetric, disk-like
states of the Reynolds stress anisotropy tensor are characteristic of mixing
layers. Similar behavior has been noted in turbulent flows over transverse bar
roughness [63], k-type roughness [64], irregular roughness [55], and recently
in AM roughness [34]. Beyond the roughness mean plane, the trajectory shifts
back towards the right side of the triangle, resembling an axisymmetric, rod-
like state, mirroring the behavior of the smooth-wall case once the wall-
normal coordinate surpasses the maximum roughness height. The most prevalent
anisotropic states include axisymmetric expansion, one-component, two-
component, and two-component axisymmetric states, with the likelihood of a
specific turbulent state increasing with higher $R_{a}/D$.

The influence of $s_{k}$ values on Reynolds stress anisotropy is evident. For
instance, comparing the PS4 surface in Fig. 14(d) with the NS4 surface in Fig.
15(d), both sharing the same $R_{a}/D$ but differing $s_{k}$ values,
highlights this impact. The disparity in turbulent states between these
surfaces stems from distinct flow interactions with surface irregularities.
Surfaces characterized by positive and negative $s_{k}$ values represent two
rough surface types: one dominated by peaks and the other by valleys,
respectively. Peak-dominated surfaces exhibit protruding irregularities,
promoting turbulence generation through flow acceleration, vorticity, and
separation, resulting in a more intense and anisotropic turbulent flow.
Conversely, valley-dominated surfaces have recessed features that dampen
turbulence intensity by inducing flow deceleration and recirculation,
resulting in a less intense and more isotropic turbulent flow. Thus, the
physical mechanisms driving turbulence generation and anisotropy differ
between peak-dominated and valley-dominated surfaces, explaining the observed
differences in turbulent states.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 14: Anisotropy-invariant mapping of turbulence in (a) $R_{a}/D=0$
turbulent smooth pipe flow and rough pipe flow compiled from the present LES
data at: (b) $R_{a}/D=0.001$, (c) $R_{a}/D=0.012$, (d) $R_{a}/D=0.015$, (e)
$R_{a}/D=0.020$, (f) $R_{a}/D=0.028$, for positively skewed surfaces. The data
points for each case are based on all cells in the domain at $x/D=4$ and
colored with normalized wall distance values, $r^{+}$. Colormap varies from
purple (minimum) to red (maximum).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 15: Anisotropy-invariant mapping of turbulence in (a) $R_{a}/D=0$
turbulent smooth pipe flow and rough pipe flow compiled from the present LES
data at: (b) $R_{a}/D=0.001$, (c) $R_{a}/D=0.012$, (d) $R_{a}/D=0.015$, (e)
$R_{a}/D=0.020$, (f) $R_{a}/D=0.028$, for negatively skewed surfaces. The data
points for each case are based on all cells in the domain at $x/D=4$ and
colored with normalized wall distance values, $r^{+}$. Colormap varies from
purple (minimum) to red (maximum).

##  V Conclusions

We conducted a detailed examination of how wall roughness influences turbulent
heat transfer in Additively Manufactured (AM) rough surfaces using roughness-
resolved high-fidelity Large Eddy Simulations (LES) in OpenFOAM 7. We created
six configurations of rough pipes from a single actual AM surface, maintaining
fixed skewness and kurtosis while varying roughness height distributions.
Additionally, we flipped these six surfaces to produce surfaces with fixed
roughness height and kurtosis but with negative skewness. Precise spatial
measurements at a constant bulk Reynolds number $Re_{b}=8000$ enabled us to
numerically estimate the roughness function for all cases, which was then
utilized to approximate the equivalent sand-grain roughness height, $k_{s}$.
The temperature was treated as a passive scalar with a Prandtl number of 0.71,
neglecting buoyancy effects. Consistent with our previous research [34; 29],
we found that wall roughness affects heat transfer and momentum differently.
The temperature and momentum wall rough functions ($\Delta\Theta^{+}$, $\Delta
U^{+}$) differed significantly, with the former being notably smaller than the
latter. This discrepancy arises from high-temperature fluid from the bulk
region penetrating the roughness layer, resulting in a larger wall-scaled mean
temperature profile compared to the mean temperature profile, which is
predominantly negative due to pressure effects within the roughness sublayer.
Consequently, normalized temperature values in the bulk region are larger than
the normalized mean velocity values. The discrepancy between
$\Delta\Theta^{+}$ and $\Delta U^{+}$ directly challenges the validity of the
Reynolds analogy under fully rough conditions, consistent with existing
literature regardless of roughness nature, as the fact that $\Delta
U^{+}>\Delta\Theta^{+}$ can be attributed to the influence of pressure on the
velocity field without a corresponding mechanism for the thermal field. This
was further corroborated by the rapid increase in the effective Prandtl number
($Pr_{eff}$) within the rough wall, where the effective thermal diffusivity,
due to combined turbulence and dispersion effects, is significantly smaller
than the effective diffusivity within the rough wall. With increasing surface
roughness height, $Pr_{eff}$ values increase near the roughness sublayer,
indicating a reduction in the Reynolds analogy. The influence of skewness was
found to be insignificant on the overall behavior. Furthermore, while the
wall-normal Reynolds shear stress
($\left<u^{\prime}_{v}u^{\prime}_{r}\right>^{+}$) and heat flux
($\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$) decreased with larger wall
roughness height, their magnitudes remained similar for different values of
surface roughness height ($R_{a}/D$). However, it is worth noting that the
magnitude of $\left<u^{\prime}_{r}\Theta^{\prime}\right>^{+}$ was slightly
smaller than $\left<u^{\prime}_{x}u^{\prime}_{r}\right>^{+}$, especially in
the region just behind the roughness crest, highlighting the influence of
recirculation bubbles on reducing heat transfer. The impact of skewness became
more apparent for larger values of $R_{a}/D>0.006$. For a fixed $R_{a}/D$
value, negatively skewed surfaces exhibited smaller $k_{s}^{+}$ values
compared to positively skewed ones, indicating less turbulence in surfaces
dominated by cavities and consequently less heat transfer. This was further
quantified by visualizing turbulence states of
$\left<u^{\prime}_{x}u^{\prime}_{r}\right>^{+}$, where peak-dominated surfaces
demonstrated a higher probability of flow acceleration compared to valley-
dominated ones. Evaluations of global features such as friction factor
($c_{f}/c_{f}^{0}$) and Nusselt number ($Nu/Nu^{0}$) demonstrated a clear
dependence on surface roughness height and skewness. The influence of surface
skewness was only evident on $c_{f}/c_{f}^{0}$ in the fully rough regime,
while its effect on $Nu/Nu^{0}$ was noticeable in the transitional rough
regime. Thermal performance evaluation further indicated that for efficient
heat transfer enhancement, positively skewed surfaces outperformed negatively
skewed ones in both transitional and fully rough flow regimes.

## Acknowledgement

The authors greatly appreciate the financial support provided by Vinnova under
project number "2020-04529". Computer time was provided by the Swedish
National Academic Infrastructure for Supercomputing (NAISS), partially funded
by the Swedish Research Council through grant agreement no. 2018-05973. The
authors also appreciate funding acquisition partners Setrab and Siemens Energy
AB, and Sebastian Ritcher for conducting the surface roughness measurements.

## Author declarations

### Conflict of interest

The authors have no conflicts of interest to disclose.

### Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * Bunker [2017] R. S. Bunker, “Evolution of turbine cooling,” in  _Turbo Expo: Power for Land, Sea, and Air_ , Vol. 50770 (American Society of Mechanical Engineers, 2017) p. V001T51A001.
  * Adair, Kirka, and Ryan [2019] D. Adair, M. Kirka, and D. Ryan, “Additive manufacture of prototype turbine blades for hot-fired engine performance validation trials,” in  _Turbo Expo: Power for Land, Sea, and Air_ , Vol. 58677 (American Society of Mechanical Engineers, 2019) p. V006T24A012.
  * Chowdhury  _et al._ [2022] S. Chowdhury, N. Yadaiah, C. Prakash, S. Ramakrishna, S. Dixit, L. R. Gulta, and D. Buddhi, “Laser powder bed fusion: A state-of-the-art review of the technology, materials, properties & defects, and numerical modelling,” Journal of Materials Research and Technology (2022).
  * Ngo  _et al._ [2018] T. D. Ngo, A. R. Kashani, G. Imbalzano, K. T. Nguyen, and D. Hui, “Additive manufacturing (3d printing): A review of materials, methods, applications and challenges,” Composites Part B: Engineering (2018).
  * Colebrook  _et al._ [1939] C. F. Colebrook, T. Blench, H. Chatley, E. H. Essex, J. R. Finniecome, G. Lacey, J. Williamson, and G. G. Macdonald, “Turbulent flow in pipes, with particular reference to the transition region between the smooth and rough pipe laws,” J. Inst. Civ. Eng. 12, 393–422 (1939).
  * Schlichting [1937] H. Schlichting,  _Experimental investigation of the problem of surface roughness_ , 823 (National Advisory Commitee for Aeronautics, 1937).
  * Nikuradse [1933] J. Nikuradse, “Laws of flow in rough pipes,” (1933).
  * Clauser [1954] F. H. Clauser, “Turbulent boundary layers in adverse pressure gradients,” J. Aeronaut. Sci. 21, 91–108 (1954).
  * Dittus [1930] F. W. Dittus, “Heat transfer in automobile radiators of the tubular type,” Univ. of California Pub., Eng. 2, 443–461 (1930).
  * Gnielinski [1975] V. Gnielinski, “New equations for heat and material transfer in turbulent flow pipes and channels,” Research in Engineering A 41, 8–16 (1975).
  * Nunner [1958] W. Nunner,  _Heat transfer and pressure drop in rough tubes_ (Atomic Energy Research Establishment, 1958).
  * Dipprey and Sabersky [1963] D. Dipprey and R. Sabersky, “Heat and momentum transfer in smooth and rough tubes at various prandtl numbers,” Int. J. Heat Mass Transf. 6, 329–353 (1963).
  * Kays, Crawford, and Weigand [1980] W. M. Kays, M. E. Crawford, and B. Weigand,  _Convective heat and mass transfer_ , Vol. 4 (McGraw-Hill New York, 1980).
  * Kays, Crawford, and Weigand [2005] W. M. Kays, M. E. Crawford, and B. Weigand,  _Convective heat and mass transfer_ , Vol. 4 (McGraw-Hill New York, 2005).
  * Bons [2002] J. P. Bons, “$st$ and $c_{f}$ augmentation for real turbine roughness with elevated freestream turbulence,” in  _ASME Turbo Expo._ , Vol. 36088 (2002) pp. 349–363.
  * Li  _et al._ [2007] Z. Li, Y.-L. He, G.-H. Tang, and W.-Q. Tao, “Experimental and numerical studies of liquid flow and heat transfer in microtubes,” Int. J Heat Mass Transf. 50, 3447–3460 (2007).
  * Pelević and van der Meer [2016] N. Pelević and T. H. van der Meer, “Heat transfer and pressure drop in microchannels with random roughness,” Int. J Thermal Sci. 99, 125–135 (2016).
  * Lu  _et al._ [2020] H. Lu, M. Xu, L. Gong, X. Duan, and J. C. Chai, “Effects of surface roughness in microchannel with passive heat transfer enhancement structures,” International Journal of Heat and Mass Transfer 148, 119070 (2020).
  * Ansari and Zhou [2020] M. Q. Ansari and G. Zhou, “Influence of structured surface roughness peaks on flow and heat transfer performances of micro-and mini-channels,” Int. Commun. Heat Mass Transf. 110, 104428 (2020).
  * Croce, D’agaro, and Nonino [2007] G. Croce, P. D’agaro, and C. Nonino, “Three-dimensional roughness effect on microchannel heat transfer and pressure drop,” Int. J Heat Mass Transf. 50, 5249–5259 (2007).
  * Xiong and Chung [2010] R. Xiong and J. N. Chung, “Investigation of laminar flow in microtubes with random rough surfaces,” Microfluidics and Nanofluidics 8, 11–20 (2010).
  * Kadivar, Tormey, and McGranaghan [2022] M. Kadivar, D. Tormey, and G. McGranaghan, “Cfd of roughness effects on laminar heat transfer applied to additive manufactured minichannels,” Heat Mass Transf. , 1–15 (2022).
  * Kadivar, Tormey, and McGranaghan [2023] M. Kadivar, D. Tormey, and G. McGranaghan, “A comparison of rans models used for cfd prediction of turbulent flow and heat transfer in rough and smooth channels,” International Journal of Thermofluids 20, 100399 (2023).
  * Garg, Wang, and Fureby [2023] H. Garg, L. Wang, and C. Fureby, “Large-eddy simulations of separated flow and heat transfer in a ribbed channel,” in  _ICHMT Digital Library Online_ (Begel House Inc., 2023).
  * Forooghi  _et al._ [2018] P. Forooghi, A. Stroh, P. Schlatter, and B. Frohnapfel, “Direct numerical simulation of flow over dissimilar, randomly distributed roughness elements: A systematic study on the effect of surface morphology on turbulence,” Phys. Rev. Fluids 3, 044605 (2018).
  * Peeters and Sandham [2019] J. W. R. Peeters and N. D. Sandham, “Turbulent heat transfer in channels with irregular roughness,” Int. J. Heat Mass Transf. 138, 454–467 (2019).
  * MacDonald  _et al._ [2019] M. MacDonald, N. Hutchins, D. Lohse, and D. Chung, “Heat transfer in rough-wall turbulent thermal convection in the ultimate regime,” Phys. Rev. Fluids 4, 071501 (2019).
  * Kuwata [2021] Y. Kuwata, “Direct numerical simulation of turbulent heat transfer on the reynolds analogy over irregular rough surfaces,” Int. J. Heat Fluid Flow 92, 108859 (2021).
  * Garg, Wang, and Fureby [2024] H. Garg, L. Wang, and C. Fureby, “Heat transfer enhancement with additively manufactured rough surfaces: Insights from large-eddy simulations,” Phys. Fluids 36 (2024).
  * Stimpson  _et al._ [2017] C. K. Stimpson, J. C. Snyder, K. A. Thole, and D. Mongillo, “Scaling roughness effects on pressure loss and heat transfer of additively manufactured channels,” J. Turbomach. 139 (2017).
  * Snyder and Thole [2020] J. C. Snyder and K. A. Thole, “Tailoring surface roughness using additive manufacturing to improve internal cooling,” J. Turbomach. 142 (2020).
  * McClain  _et al._ [2021] S. T. McClain, D. R. Hanson, E. Cinnamon, J. C. Snyder, R. F. Kunz, and K. A. Thole, “Flow in a simulated turbine blade cooling channel with spatially varying roughness caused by additive manufacturing orientation,” J. Turbomach. 143 (2021).
  * Favero  _et al._ [2022] G. Favero, G. Berti, M. Bonesso, D. Morrone, S. Oriolo, P. Rebesan, R. Dima, P. Gregori, A. Pepato, A. Scanavini, and S. Mancin, “Experimental and numerical analyses of fluid flow inside additively manufactured and smoothed cooling channels,” Int. Commun. Heat Mass Trans. 135, 106128 (2022).
  * Garg  _et al._ [2023] H. Garg, L. Wang, G. Sahut, and C. Fureby, “Large eddy simulations of fully developed turbulent flows over additively manufactured rough surfaces,” Physics of Fluids 35 (2023).
  * Note [1] https://github.com/CoffeeDynamics/STLRoughPipes.
  * Nicoud and Ducros [1999] F. Nicoud and F. Ducros, “Subgrid-scale stress modelling based on the square of the velocity gradient tensor,” Flow Turbul. Combust. 62, 183–200 (1999).
  * Moin  _et al._ [1991] P. Moin, K. Squires, W. Cabot, and S. Lee, “A dynamic subgrid-scale model for compressible turbulence and scalar transport,” Phys. Fluids A 3, 2746–2757 (1991).
  * Kasagi, Tomita, and Kuroda [1992] N. Kasagi, Y. Tomita, and A. Kuroda, “Direct numerical simulation of passive scalar field in a turbulent channel flow,” ASME J. Heat Trans. 114, 598–606 (1992).
  * Kozuka, Seki, and Kawamura [2009] M. Kozuka, Y. Seki, and H. Kawamura, “DNS of turbulent heat transfer in a channel flow with a high spatial resolution,” Int. J. Heat Fluid Flow 30, 514–524 (2009).
  * Lluesma-Rodríguez, Hoyas, and Perez-Quiles [2018] F. Lluesma-Rodríguez, S. Hoyas, and M. J. Perez-Quiles, “Influence of the computational domain on DNS of turbulent heat transfer up to Reτ= 2000 for Pr= 0.71,” Int. J. Heat Mass Transf. 122, 983–992 (2018).
  * Butcher [2016] J. C. Butcher,  _Numerical methods for ordinary differential equations_ (John Wiley & Sons, 2016).
  * Weller  _et al._ [1998] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby, “A tensorial approach to computational continuum mechanics using object-oriented techniques,” Comp. Phys. 12, 620–631 (1998).
  * Ruge and Stüben [1987] J. W. Ruge and K. Stüben, “4. algebraic multigrid,” in  _Multigrid Methods_ (1987) pp. 73–130.
  * Golub and Van Loan [1996] G. H. Golub and C. F. Van Loan, “Matrix computations (Johns Hopkins studies in mathematical sciences),” Matrix Comput. (1996).
  * Press  _et al._ [2007] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical recipes 3rd edition: The art of scientific computing_ (Cambridge university press, 2007).
  * Lozano and Jiménez [2014] D. A. Lozano and J. Jiménez, “Effect of the computational domain on direct simulations of turbulent channels up to re $\tau$= 4200,” Phys. Fluids 26, 011702 (2014).
  * Lluesma, Hoyas, and Perez [2018] R. F. Lluesma, S. Hoyas, and M. J. Q. Perez, “Influence of the computational domain on dns of turbulent heat transfer up to re$\tau$= 2000 for pr= 0.71,” Int. J. Heat Mass Transf. 122, 983–992 (2018).
  * Garg  _et al._ [2024] H. Garg, L. Wang, M. Andersson, and C. Fureby, “Large eddy simulations of turbulent pipe flows at moderate reynolds numbers,” Physics of Fluids 36 (2024).
  * Raupach [1994] M. R. Raupach, “Simplified expressions for vegetation roughness length and zero-plane displacement as functions of canopy height and area index,” Boundary Layer Meteorol. 71, 211–216 (1994).
  * MacDonald, Hutchins, and Chung [2019] M. MacDonald, N. Hutchins, and D. Chung, “Roughness effects in turbulent forced convection,” J. Fluid Mech. 861, 138–162 (2019).
  * Forooghi, Stripf, and Frohnapfel [2018] P. Forooghi, M. Stripf, and B. Frohnapfel, “A systematic study of turbulent heat transfer over rough walls,” Int. J. Heat Mass Transf. 127, 1157–1168 (2018).
  * Busse, Thakkar, and Sandham [2017] A. Busse, M. Thakkar, and N. D. Sandham, “Reynolds-number dependence of the near-wall flow over irregular rough surfaces,” J. Fluid Mech. 810, 196–224 (2017).
  * Thakkar, Busse, and Sandham [2017] M. Thakkar, A. Busse, and N. Sandham, “Surface correlations of hydrodynamic drag for transitionally rough engineering surfaces,” J. Turb. 18, 138–169 (2017).
  * Barros and Christensen [2019] J. M. Barros and K. T. Christensen, “Characteristics of large-scale and superstructure motions in a turbulent boundary layer overlying complex roughness,” J. Turb. 20, 147–173 (2019).
  * Busse and Jelly [2020] A. Busse and T. O. Jelly, “Influence of surface anisotropy on turbulent flow over irregular roughness,” Flow Turb. Combust. 104, 331–354 (2020).
  * Kuwata [2022] Y. Kuwata, “Reynolds number dependence of turbulent heat transfer over irregular rough surfaces,” Phys. Fluids 34 (2022).
  * Kader [1981] B. A. Kader, “Temperature and concentration profiles in fully turbulent boundary layers,” Int. J. Heat Mass Transf. 24, 1541–1544 (1981).
  * Spalart and Strelets [2000] P. R. Spalart and M. K. Strelets, “Mechanisms of transition and heat transfer in a separation bubble,” J. Fluid Mech. 403, 329–349 (2000).
  * Guezennec, Stretch, and Kim [1990] Y. Guezennec, D. Stretch, and J. Kim, “The structure of turbulent channel flow with passive scalar transport,” Studying Turbulence Using Numerical Simulation Databases. 3: Proceedings of the 1990 Summer Program , 127–138 (1990).
  * Abe and Antonia [2017] H. Abe and R. A. Antonia, “Relationship between the heat transfer law and the scalar dissipation function in a turbulent channel flow,” J. Fluid Mech. 830, 300–325 (2017).
  * Yuan and Piomelli [2014] J. Yuan and U. Piomelli, “Roughness effects on the reynolds stress budgets in near-wall turbulence,” J. Fluid Mech. 760, R1 (2014).
  * Raupach and Shaw [1982] M. R. Raupach and R. H. Shaw, “Averaging procedures for flow within vegetation canopies,” Bound. Layer Meteorol. 22, 79–90 (1982).
  * Ashrafian and Andersson [2006] A. Ashrafian and H. I. Andersson, “The structure of turbulence in a rod-roughened channel,” Int. J. Heat Fluid Flow 27, 65–79 (2006).
  * Smalley  _et al._ [2002] R. Smalley, S. Leonardi, R. Antonia, L. Djenidi, and P. Orlandi, “Reynolds stress anisotropy of turbulent rough wall layers,” Exp. Fluids 33, 31–37 (2002).