2.1. The dimensional equations

As a major assumption, the changes in the density of the fluid are considered to only affect the buoyancy term in the momentum conservation equation. Under the Boussinesq approximation, the conservation equations in terms of the dimensional variables, space coordinates $\hat {x}_{i}$, time $\hat {t}$, pressure $\hat {P}$, velocity $\hat {u}_{i}$ and temperature $\hat {T}$, are expressed as follows:

(2.1a)\begin{gather} \frac{\partial{\hat{u}_{i}}}{\partial{\hat{x}_{i}}}= 0, \end{gather}

(2.1b)\begin{gather}\hat{\rho}_{\infty}\left(\frac{\partial \hat{u}_{i}}{\partial \hat{t}}+\hat{u}_{j} \frac{\partial{\hat{u}_{i}}}{\partial {\hat x_{j}}}\right)={-}\frac{\partial({\hat{P}-\hat{P}_{\infty}})}{\partial{\hat{x}_{i}}}+\mu \frac{\partial^2 \hat{u}_{i}}{\partial {\hat{x}_{j}}^2} -\hat{\rho}_{\infty}\beta(\hat{T}-\hat{T}_{\infty})g_{i}, \end{gather}

(2.1c)\begin{gather}\frac{\partial \hat{T}}{\partial \hat{t}}+\hat{u}_{j} \frac{\partial{\hat T}}{\partial {\hat x_{j}}}= \alpha\frac{\partial^2 \hat{T}}{\partial {\hat{x}_{j}}^2}, \end{gather}

with $\hat \rho _\infty$, $\hat P_\infty$ and $\hat T_\infty$ the density, pressure and temperature in the stagnant flow region, sufficiently far away from the vertical wall. The parameters assumed constant in the equations above are the volumetric thermal expansion coefficient, $\beta$, the dynamic viscosity, $\mu =\hat \rho _\infty \nu$, with $\nu$ the kinematic viscosity, and the thermal diffusivity, $\alpha$. With the axes as in figure 2, the volume force per unit mass has components $g_i=-g \delta _{i1}$ with $g$ the gravitational acceleration and $\delta _{ij}$ the Kronecker index. The parameter controlling the thermal convection flow is the Rayleigh number $Ra$, defined as

(2.2)\begin{equation} Ra=\frac{g \beta (\hat T_w-\hat T_\infty)L^3}{\alpha \nu}, \end{equation}

where the temperature of the wall, $\hat T_w$, is maintained constant, and the plate height, $L$, is the macroscopic length scale of the problem. We also define the Grashof number, $Gr=Ra/Pr$, with $Pr=\nu /\alpha$ the Prandtl number, a property of the fluid. Given the presence of two characteristic length scales, a macroscopic and a microscopic one, the latter related to the periodicity $l$ of the microstructures present on the vertical surface, two problems will be set up. These two problems will be coupled at some distance from the wall, a distance that is asymptotically large when seen from the microscopic point of view and asymptotically small when seen from the macroscopic viewpoint.

Figure 2. Sketch of a general vertical rough surface, periodically micro-patterned, with notations and indication of microscopic and macroscopic domains.

2.2. The macroscale problem

To set the proper scales of the macroscopic problem we consider the fact that the motion of the fluid is generated by the buoyancy force; if $\mathcal {U}$ is the characteristic velocity of the fluid, we can write

(2.3)\begin{equation} \frac{\hat \rho_\infty \mathcal{U}^2}{L}\sim \hat \rho_\infty \beta (\hat T_w-\hat T_\infty) g. \end{equation}

We thus define the velocity scale $\mathcal {U}=\sqrt {\beta (\hat T_w-\hat T_\infty ) g L} = Gr^{1/2} \frac {\nu }{L}$ and normalize the velocity vector as

(2.4)\begin{equation} U_i=\frac{\hat u_i}{\mathcal{U}}. \end{equation}

The other dimensionless variables are defined as follows:

(2.5a–c)\begin{equation} X_i=\frac{\hat x_i}{L}, \quad t = \frac{\hat t \mathcal{U}}{L}, \quad P=\frac{\hat P-\hat P_\infty}{\hat \rho_\infty \mathcal{U}^2}, \quad \varTheta= \frac{\hat T-\hat T_\infty}{\hat T_w-\hat T_\infty}, \end{equation}

for the balance equations to become

(2.6a)\begin{gather} \frac{\partial{{U}_{i}}}{\partial{{X}_{i}}}= 0, \end{gather}

(2.6b)\begin{gather}\frac{\partial{U_i}}{\partial {t}}+ U_j \frac{\partial{U_i}}{\partial {X_j}}={-}\frac{\partial{P}}{\partial{X}_{i}}+\frac{1}{\sqrt{Gr}}\frac{\partial^2 {U}_{i}}{\partial {X}_{j}^2} +{\varTheta}\delta_{i1}, \end{gather}

(2.6c)\begin{gather}Pr \sqrt{Gr} \left(\frac{\partial{\varTheta}}{\partial {t}}+ U_j \frac{\partial{\varTheta}}{\partial {X_j}}\right)= \frac{\partial^2 {\varTheta}}{\partial {X}_{j}^2}. \end{gather}

These equations depend on only the macroscopic independent variables, $t$ and $X_i$, and must be solved subject to matching conditions at $X_2 \to 0$, together with $\varTheta =U_1=0$ and ${\partial U_2}/{\partial X_2}=0$ for $X_2 \to \infty$.

2.3. The microscale problem

The near-wall problem differs from the previous one in that the microscopic velocity scale is taken to be $\epsilon \, \mathcal {U}$, with $\epsilon =\ell /L\ll 1$. Also, the pressure scale for the near-wall flow is the viscous pressure, i.e. $\mu (\epsilon \,{\mathcal {U}})/l$. Dimensionless variables in the microscopic domain are introduced as follows:

(2.7a–e)\begin{equation} u_i=\frac{\hat u_i}{\epsilon \, \mathcal{U}}, \quad x_i=\frac{\hat x_i}{l}, \quad t = \frac{\hat t \mathcal{U}}{L}, \quad p=\frac{(\hat P-\hat P_\infty)\,L}{\mu \mathcal{U}}, \quad \theta= \frac{\hat T-\hat T_\infty}{\hat T_w-\hat T_\infty}. \end{equation}

The microscopic dimensionless equations are

(2.8a)\begin{gather} \frac{\partial u_i}{\partial x_i}=0, \end{gather}

(2.8b)\begin{gather}\epsilon \mathcal{R}_G \left(\frac{\partial u_i}{\partial t}+ u_j \frac{\partial u_i}{\partial x_j} \right)={-}\frac{\partial p}{\partial x_i}+ \frac{\partial^2 u_i}{\partial x_j^2}+\mathcal{R}_G \theta \delta_{i1}, \end{gather}

(2.8c)\begin{gather}\epsilon \mathcal{R}_G \,Pr \left(\frac{\partial \theta}{\partial t}+u_j \frac{\partial \theta}{\partial x_j}\right)=\frac{\partial^2 \theta}{\partial x_j^2}, \end{gather}

with the reduced Grashof number $\mathcal {R}_G$, defined by $\mathcal {R}_G= \epsilon \sqrt {Gr}$, assumed of order one. The microscale problem is bounded by the microstructured wall on one side; therefore, the following condition is imposed at this location:

(2.8d)\begin{equation} u_i=0,\quad \theta = 1\quad {\rm{for}} \ x_2=y_w, \end{equation}

with $y_{w}=y_{w}(x_1,x_3)$ the micro-patterned wall. A representative volume element must be chosen, of unit length along $x_1$ and $x_3$ (cf. figure 2), and periodic conditions are enforced for all dependent variables along these directions. On account of the scalings adopted for inner (i.e. near-wall) and outer problems, the conditions for $x_2\to \infty$ are

(2.8e)\begin{gather} -p\delta_{i2}+\left(\frac{\partial u_i}{\partial x_2}+\frac{\partial u_2}{ \partial x_i}\right)= S_{i2}, \end{gather}

(2.8f)\begin{gather}\frac{\partial \theta}{\partial x_2}=\epsilon \,\eta; \end{gather}

these amount to matching the components of the traction vector and of the heat flux between the two regions. For ease of notation, in the equations above we have introduced the following definitions for the macroscopic dimensionless stresses in the streamwise, normal and spanwise directions (respectively $S_{12}, S_{22}, S_{32}$) as well as the macroscopic dimensionless normal temperature gradient ($\eta$):

(2.9)\begin{gather} S_{i2}={-}Gr^{1/2} P \delta_{i2}+\left(\frac {\partial U_i}{\partial X_2}+\frac {\partial U_2}{\partial X_i}\right), \end{gather}

(2.10)\begin{gather}\eta=\frac{\partial \varTheta}{\partial X_2}. \end{gather}

Notice that both $S_{i2}$ and $\eta$ depend on only macroscopic variables; they represent the forcing of the outer flow on the near-wall state.

We still need to specify the asymptotic matching conditions which will eventually result in effective boundary conditions for the macroscopic problem, to be applied some distance from the microstructured wall. They are:

(2.11a)\begin{gather} \lim_{X_2\to 0} U_i = \lim_{x_2\to\infty} \epsilon u_i, \end{gather}

(2.11b)\begin{gather}\lim_{X_2\to 0} \varTheta = \lim_{x_2\to\infty} \theta. \end{gather}

3.1. Expansion of the inner variables

Asymptotic expansions in terms of the small parameter $\epsilon$ are introduced, and like-order terms are collected, leading to a hierarchy of problems. We impose

(3.1)\begin{equation} u_i=u_i^{(0)}+\epsilon\,u_i^{(1)}+\epsilon^2 \, u_i^{(2)}+ \cdots, \end{equation}

and likewise for $p$ and $\theta$. Furthermore, using the chain rule, we replace in the microscopic equations the term ${\partial }/{\partial x_i}$ by ${\partial }/{\partial x_i}+\epsilon \,({\partial }/{\partial X_i})$. The asymptotic expressions are plugged into (2.8a) to (2.8f) governing the microscale problem.

3.2. Reconstruction of the problem at different orders

The problems at the asymptotic orders of interest are given below.

3.2.1. The ${O}(\epsilon ^0)$ problem

We have

(3.2a)\begin{gather} \frac{\partial u_i^{(0)}}{\partial x_i}=0, \end{gather}

(3.2b)\begin{gather}-\frac{\partial p^{(0)}}{\partial x_i}+ \frac{\partial^2 u_i^{(0)}}{\partial x_j^2}+\mathcal{R}_G \theta^{(0)} \delta_{i1}=0, \end{gather}

(3.2c)\begin{gather}\frac{\partial^2 \theta^{(0)}}{\partial x_j^2}=0, \end{gather}

with boundary conditions

(3.2d)\begin{gather} u_i^{(0)}=0, \quad \theta^{(0)}=1 \quad {\rm at}\ x_2 = y_w, \end{gather}

(3.2e)\begin{gather}-p^{(0)}\delta_{i2}+\left(\frac{\partial u_i^{(0)}}{\partial x_2}+\frac{\partial u_2^{(0)}}{\partial x_i}\right)=S_{i2}, \quad \frac{\partial \theta^{(0)}}{\partial x_2}=0 \quad {\rm for} \ x_2\to\infty. \end{gather}

A solution of this problem can be sought by separation of variables, on account of the linearity of the system, for the solution to take the form

(3.3a,b)\begin{equation} u_i^{(0)}=\breve u_{ik}S_{k2} + u_i^{\dagger} \mathcal{R}_G, \quad p^{(0)}=\breve{p}_{k} S_{k2} + p^{\dagger} \mathcal{R}_G + P_0, \end{equation}

with $\breve u_{ik}$, $u_i^{\dagger}$, $\breve {p}_{k}$ and $p^{\dagger}$ tensors which depend on microscopic variables only, and $P_0$ an integration constant function only of $X_j$. After plugging the ansatz for the order-zero solution into the balance equations, it becomes clear that uniqueness conditions are needed for $\breve p_k$ and $p^{\dagger}$, which appear in the system only through their gradients. We enforce the vanishing of the integrals of $\breve p_k$ and $p^{\dagger}$ over a cubic cell of unit side length positioned sufficiently far from the wall (nominally for $x_2\to \infty$); this leads to the vanishing of $P_0$. It is also clear that we cannot stop the solution at this order, since the leading-order temperature solution is simply $\theta ^{(0)}=1$, i.e. the effect of the microstructure appears in the temperature at the next $\epsilon$ order.

The dynamic problem at ${O}(\epsilon ^0)$ yields the same equations for $\breve u_{ik}$ and $\breve p_k$ already given for the isothermal case by Bottaro & Naqvi (Reference Bottaro and Naqvi2020), so that we can anticipate that the first correction to the no-slip condition for the velocity will be a Navier-slip term. Such a leading-order problem reads

(3.4a)\begin{gather} \frac{\partial \breve u_{ik}}{\partial x_i}=0, \end{gather}

(3.4b)\begin{gather}-\frac{\partial \breve p_k}{\partial x_i}+ \frac{\partial^2 \breve u_{ik}}{\partial x_j^2}=0, \end{gather}

with

(3.4c)\begin{gather} \breve u_{ik}=0 \quad {\rm at} \ x_2 = y_w, \end{gather}

(3.4d)\begin{gather}-\breve p_k\,\delta_{i2}+\left(\frac{\partial \breve u_{ik}}{\partial x_2}+\frac{\partial \breve u_{2k}}{\partial x_i}\right)=\delta_{ik} \quad {\rm at} \ x_2\to\infty. \end{gather}

The ${\dagger}$ variables, which describe the effect of buoyancy on velocity and pressure fields, satisfy the steady system

(3.5a)\begin{gather} \frac{\partial u_{i}^{\dagger}}{\partial x_i}=0, \end{gather}

(3.5b)\begin{gather}-\frac{\partial p^{\dagger}}{\partial x_i}+\frac{\partial^2 u_{i}^{\dagger}}{\partial x_j^2}={-}\delta_{i1}, \end{gather}

with

(3.5c)\begin{gather} u_{i}^{\dagger}=0 \quad {\rm at} \ x_2 = y_w, \end{gather}

(3.5d)\begin{gather}- p^{\dagger} \delta_{i2}+\left(\frac{\partial u_{i}^{\dagger}}{\partial x_2}+\frac{\partial u_{2}^{\dagger}}{\partial x_i}\right)=0 \quad {\rm at} \ x_2 \to \infty. \end{gather}

As it will be shown later on, the problems can be further simplified when $x_3$-elongated wall ribs are examined, as in the case of riblets (Bechert & Bartenwerfer Reference Bechert and Bartenwerfer1989; Luchini, Manzo & Pozzi Reference Luchini, Manzo and Pozzi1991).

3.2.2. The ${O}(\epsilon ^1)$ problem

The equations at order $\epsilon$ are forced by the order-one state, i.e.

(3.6a)\begin{gather} \frac{\partial u_i^{(1)}}{\partial x_i}={-}\frac{\partial u_i^{(0)}}{\partial X_i}, \end{gather}

(3.6b)\begin{gather}-\frac{\partial p^{(1)}}{\partial x_i}+ \frac{\partial^2 u_i^{(1)}}{\partial x_j^2}+\mathcal{R}_G \theta^{(1)} \delta_{i1}= \frac{\partial p^{(0)}}{\partial X_i}-2 \frac{\partial^2 u_i^{(0)}}{\partial x_j \partial X_j}+\mathcal{R}_G \left(\frac{\partial u_i^{(0)}}{\partial t}+u_j^{(0)}\frac{\partial u_i^{(0)}}{\partial x_j}\right), \end{gather}

(3.6c)\begin{gather}\frac{\partial^2 \theta^{(1)}}{\partial x_j^2}={-}2\frac{\partial^2 \theta^{(0)}}{\partial x_j \partial X_j}+\mathcal{R}_G\, Pr\left(\frac{\partial \theta^{(0)}}{\partial t}+u_j^{(0)}\frac{\partial \theta^{(0)}}{\partial x_j}\right), \end{gather}

with boundary conditions

(3.6d)\begin{gather} u_i^{(1)}= \theta^{(1)}=0 \quad {\rm at}\ x_2 = y_w, \end{gather}

(3.6e)\begin{gather}-p^{(1)}\delta_{i2}+\left(\frac{\partial u_i^{(1)}}{\partial x_2}+\frac{\partial u_2^{(1)}}{\partial x_i}\right)={-}\left(\frac{\partial u_i^{(0)}}{\partial X_2}+\frac{\partial u_2^{(0)}}{\partial X_i}\right)\quad {\rm at} \ x_2\to\infty, \end{gather}

(3.6f)\begin{gather}\frac{\partial \theta^{(1)}}{\partial x_2}= \eta -\frac{\partial \theta^{(0)}}{\partial X_2} \quad {\rm at} \ x_2\to\infty. \end{gather}

We must now substitute the results for $u_i^{(0)}$, $p^{(0)}$ and $\theta ^{(0)}$ into (3.6a) to (3.6f). As a first step, a solution for $\theta ^{(1)}$ is to be sought from the energy equation and the corresponding boundary conditions. Specifically, these equations read

(3.7a)\begin{gather} \frac{\partial^2 \theta^{(1)}}{\partial x_i^2}=0, \end{gather}

(3.7b)\begin{gather}\theta^{(1)}=0 \quad {\rm at}\ x_2 = y_w, \end{gather}

(3.7c)\begin{gather}\frac{\partial \theta^{(1)}}{\partial x_2}=\eta \quad {\rm at} \ x_2\to\infty. \end{gather}

Owing to linearity, the solution can be written as

(3.8)\begin{equation} \theta^{(1)}= \tilde{\theta} (x_i) \eta (X_i). \end{equation}

The new microscopic field $\tilde \theta$ solves the system

(3.9a)\begin{gather} \frac{\partial^2 \tilde \theta}{\partial x_i^2}=0, \end{gather}

(3.9b)\begin{gather}\tilde \theta=0 \quad {\rm at} \ x_2=y_w, \end{gather}

(3.9c)\begin{gather}\frac{\partial \tilde \theta}{\partial x_2}=1 \quad {\rm at} \ x_2\to\infty. \end{gather}

The equations governing the behaviour of $u_i^{(1)}$ and $p^{(1)}$ can be recast as follows:

(3.10a)\begin{gather} \frac{\partial u_i^{(1)}}{\partial x_i}={-}\breve u_{jk} \frac{\partial S_{k2}}{\partial X_j}\end{gather}

(3.10b)\begin{gather} -\frac{\partial p^{(1)}}{\partial x_i}+\frac{\partial^2 u_i^{(1)}}{\partial x_j^2}=\mathcal{R}_G^3 u_j^{\dagger} \frac{\partial u_i^{\dagger}}{\partial x_j}+\mathcal{R}_G^2 \left[\breve u_{jk} \frac{\partial u_i^{\dagger}}{\partial x_j}+ u_j^{\dagger} \frac{\partial \breve u_{ik}}{\partial x_j}\right]S_{k2}\nonumber\\ \hspace{6pc} +\mathcal{R}_G\left[\breve u_{jk} \frac{\partial u_{i \ell}}{\partial x_j}\right]S_{k2} S_{\ell 2}+ \mathcal{R}_G \breve u_{ik} \frac{\partial S_{k2}}{\partial t}\nonumber\\ \hspace{7pc} -\mathcal{R}_G \eta \tilde \theta \delta_{i1} + \breve p_k \frac{\partial S_{k2}}{\partial X_i} -2 \frac{\partial \breve u_{ik}}{\partial x_j} \frac{\partial S_{k2}}{\partial X_j}, \end{gather}

with boundary conditions

(3.10c)\begin{gather} u_i^{(1)}=0 \quad {\rm at} \ x_2=y_{w}, \end{gather}

(3.10d)\begin{gather}\left.\begin{array}{c@{}} \dfrac{\partial u_1^{(1)}}{\partial x_2}+\dfrac{\partial u_2^{(1)}}{\partial x_1}={-}\breve u_{1k} \dfrac{\partial S_{k2}}{\partial X_2} -\breve u_{2k} \dfrac{\partial S_{k2}}{\partial X_1}, \\ - p^{(1)}+2\dfrac{\partial u_2^{(1)}}{\partial x_2}={-}2 \breve u_{2k}\dfrac{\partial S_{k2}}{\partial X_2}, \\ \dfrac{\partial u_3^{(1)}}{\partial x_2}+\dfrac{\partial u_2^{(1)}}{\partial x_3}={-}\breve u_{3k} \dfrac{\partial S_{k2}}{\partial X_2} -\breve u_{2k} \dfrac{\partial S_{k2}}{\partial X_3} \quad {\rm at} \ x_2\to\infty. \end{array}\right\} \end{gather}

Again, a generic form of the solution can be sought, i.e.

(3.11a)\begin{gather} u_i^{(1)}=\dot u_{ijk} \left[\frac{\partial S_{k2}}{\partial X_j}\right] + \ddot u_{ik} [\mathcal{R}_G (S_{k2})^2]+ \ddot u_{i12} [\mathcal{R}_G S_{12}S_{22}]+\ddot u_{i13} [\mathcal{R}_G S_{12}S_{32}]\nonumber\\ \hspace{5pc} +\ddot u_{i23} [\mathcal{R}_G S_{22}S_{32}]+u_i' [\mathcal{R}_G \eta]+\bar{u}_{ik} [\mathcal{R}_G^2 S_{k2}] + u_i^ \ddagger [\mathcal{R}_G^3] +u^t_{ik} \left[\mathcal{R}_G \frac{\partial S_{k2}}{\partial t}\right], \end{gather}

(3.11b)\begin{gather} p^{(1)}=\dot p_{jk} \left[\frac{\partial S_{k2}}{\partial X_j}\right] + \ddot p_{k} [\mathcal{R}_G (S_{k2})^2]+ \ddot p_{12} [\mathcal{R}_G S_{12}S_{22}]+\ddot p_{13} [\mathcal{R}_G S_{12}S_{32}]\nonumber\\ \hspace{5pc} +\ddot p_{23} [\mathcal{R}_G S_{22}S_{32}]+p{'} [\mathcal{R}_G \eta]+\overline p_{k} [\mathcal{R}_G^2 S_{k2}] + p^ \ddagger [\mathcal{R}_G^3]+p^t_k\left[\mathcal{R}_G \frac{\partial S_{k2}}{\partial t}\right]. \end{gather}

Twenty-three decoupled systems of equations arise from substituting the preceding forms into (3.10a) to (3.10d). They are given in Appendix A.

3.2.3. Taking the temperature condition to higher order

Given that the macroscopic velocity at the matching surface is now available up to order $\epsilon ^2$ (cf. (2.11a)), it is advisable to do the same with the temperature. Employing the values of the dependent variables at the earlier orders, the microscopic energy equation at ${O}(\epsilon ^2)$ now reads

(3.12a)\begin{equation} \frac{\partial^2 \theta^{(2)}}{\partial x_i^2}= Pr \mathcal{R}_G \left[\tilde \theta \frac{\partial \eta}{\partial t}+ \breve u_{jk} \frac{\partial \tilde \theta}{\partial x_j} \eta S_{k2}+ u^{\dagger}_j \frac{\partial \tilde \theta}{\partial x_j} \eta \mathcal{R}_G\right] -2 \frac{\partial \tilde \theta}{\partial x_j}\, \frac{\partial \eta}{\partial X_j}. \end{equation}

The boundary conditions are:

(3.12b)\begin{equation} \theta^{(2)} = 0 \quad {\rm at} \ x_2=y_{w}, \quad \frac{\partial \theta^{(2)}}{\partial x_2}={-} \tilde \theta \frac{\partial \eta}{\partial X_2} \quad {\rm at} \ x_2 \to \infty. \end{equation}

The following general form for the solution of $\theta ^{(2)}$ may be assumed:

(3.13)\begin{equation} \theta^{(2)}= \theta^{!}_k \left[\frac{\partial \eta}{\partial X_k}\right] +\theta^*_k [Pr\, \mathcal{R}_G \eta S_{k2} ]+ \theta^{**} [Pr \,\mathcal{R}^2_G \eta ]+ \theta^t \left[Pr\,\mathcal{R}_G \frac{\partial \eta}{\partial t}\right]. \end{equation}

Eight decoupled systems of equations stem from substituting the latter form into (3.12a) and (3.12b); they are provided in Appendix B.

4.1. A synthesis of the microscopic results

The behaviours of the parameters of interest, those which contribute to the effective boundary conditions, are presented in figure 4, separating them into two groups according to their gradients in the $x_2$-direction (either positive or negative). At the matching interface ($x_2=y_\infty =5$), the variables contributing to the effective boundary conditions become independent of $x_1$ and take the following uniform values:

(4.1a,b)\begin{gather} \breve u_{11} = 5.0396, \quad \breve u_{33}= \tilde \theta = 5.0861, \end{gather}

(4.2a,b)\begin{gather}u_{1}^{\dagger}= \dot u_{112} ={-} \dot u_{211} = 12.7002, \quad \dot u_{332} ={-} \dot u_{233} = 12.9402, \end{gather}

(4.3a–c)\begin{gather}u_1^\prime= 43.0376, \quad u_{11}^t={-}42.4485, \quad u_{33}^t= \theta^t={-}43.8582. \end{gather}

Figure 4. Behaviour of the parameters of interest along a line in the $x_2$-direction which goes through the middle of the rib ($x_1=0$ with reference to figure 3). The displayed numerical results are for the case $e/\ell =0.25$ and $y_\infty =5$.

4.2. Effects of varying the matching interface location

The effect of changing the matching surface distance, $y_\infty$, on the values of the seven independent groups of effective parameters has been analysed with the aid of successive numerical simulations, varying $y_\infty$ from 2 to 6, as listed in table 1. An in-depth look into the table reveals that we have three categories of relations between the values of the microscopic parameters at the matching interface vs the location of the interface itself; specifically, linear, quadratic and cubic relations. Fitting the results, we propose the following expressions for the closure variables evaluated at $y_\infty$:

(4.3a,b)\begin{gather} \breve u_{11}= y_\infty +\lambda_x, \quad \breve u_{33}= \tilde \theta = y_\infty +\lambda_z, \end{gather}

(4.4)\begin{gather}\left.\begin{array}{c@{}} u_{1}^{\dagger}= \dot u_{112} ={-} \dot u_{211} = \dfrac{y_\infty^2}{2} + \lambda_x y_\infty + m_{12},\\[8pt] \dot u_{332} ={-} \dot u_{233} = \dfrac{y_\infty^2}{2} + \lambda_z y_\infty +m_{32}, \end{array}\right\} \end{gather}

(4.6) \begin{gather} \left.\begin{array}{c@{}} u_{1}^\prime\approx 2\left[\dfrac{y_\infty^3}{6} + m_{32} y_\infty^2 + \lambda_x y_\infty \right] + \mathcal{B}, \\ u_{11}^t\approx{-}2\left[\dfrac{y_\infty^3}{6} + m_{12} y_\infty^2 + \lambda_x y_\infty \right] +\mathcal{B}_{1t},\\ u_{33}^t\approx\theta^t={-}2\left[\frac{y_\infty^3}{6} + \lambda_z \frac{y_\infty^2}{2} + m_{32} y_\infty \right] +\mathcal{B}_{3t}.\end{array}\right\} \end{gather}

The dimensionless Navier-slip coefficients ($\lambda _x,\lambda _z$), surface permeability coefficients ($m_{12}, m_{32}$), velocity-flux sensitivity ($\mathcal {B}$) and time-fluctuation coefficients ($\mathcal {B}_{1t}, \mathcal {B}_{3t}$) are only dependent of the geometric parameters of the ribbed surface, $e/\ell$ in the case of square ribs. These coefficients can be calculated for any geometry of transverse ribs, once the microscopic numerical simulations are conducted with suitable values of $y_\infty$, and the results of the microscopic parameters at the matching interface are substituted in accurate fitting equations, or are extrapolated to $y_\infty = 0$.

Table 1. Microscopic results found at different values of $y_\infty$ for square ribs with $e/\ell =0.25$.

Simpler, accurate methods for the estimation of the coefficients of interest are proposed within the present framework. The Navier-slip coefficients can be calculated by running the simulations of the leading-order systems, forced by $S_{12}$ and $S_{32}$, with a suitable value of $y_\infty$ to get, respectively, the fields of $\breve u_{11}$ and $\breve u_{33}$; thereafter, the values of $\lambda _x$ and $\lambda _z$ can be found by averaging the corresponding field on the plane $x_2=0$. It is interesting that these same fields can then be employed to estimate the values of $m_{12}$ and $m_{32}$, making use of the numerical result pointed out by Bottaro & Naqvi (Reference Bottaro and Naqvi2020), i.e.

(4.7)\begin{gather} u_{1}^{\dagger}= \dot u_{112} ={-} \dot u_{211} = \int_{S_{cell}} \breve u_{11} \, {\rm d}\kern0.7pt x_1 \,{\rm d}\kern0.7pt x_2, \end{gather}

(4.8)\begin{gather}\dot u_{332} ={-} \dot u_{233} = \int_{S_{cell}} \breve u_{33} \,{\rm d}\kern0.7pt x_1 \,{\rm d}\kern0.7pt x_2, \end{gather}

with $S_{cell}$ the surface of the representative near-wall cell. The following values of the coefficients eventually arise when $e/\ell =0.25$:

(4.8a–b)\begin{gather} \lambda_x= 0.03975, \quad \lambda_z= 0.08609, \end{gather}

(4.9aa–b)\begin{gather}m_{12}= 0.002332 , \quad m_{32} = 0.009551, \end{gather}

(4.10a–c)\begin{gather}\mathcal{B}=0.0002399, \quad \mathcal{B}_{1t}={-}0.0000839, \quad \mathcal{B}_{3t}={-}0.0007794. \end{gather}

4.3. The formal expressions of the effective boundary conditions

The expressions of the microscopic dimensionless velocity components are now available up to ${O} (\epsilon ^1)$, while the microscopic dimensionless temperature ($\theta$) is known up to ${O} (\epsilon ^2)$. The values of the preceding quantities can be linked to the corresponding dimensionless macroscopic parameters at the matching interface, based on the concept of continuity of velocity (2.11a) and temperature (2.11b). In particular, it is convenient to enforce the conditions on the outer rim of the ribs, which amounts to specifying $x_2 = \epsilon X_2 = 0$ in the matching relations ((2.11a), (2.11b)), along with setting $y_\infty =0$ in the fits of the microscopic parameters (given in § 4.2) entering the effective boundary conditions. Eventually, we obtain

(4.11a)\begin{align} U_1|_{X_2=0} &= \epsilon [\lambda_x S_{12}+ m_{12} \mathcal{R}_G ]_{X_2=0} \nonumber\\ &\quad + \epsilon^2 \left[m_{12}\frac{\partial S_{22}}{\partial X_1}+ \mathcal{B} \mathcal{R}_G \frac{\partial \varTheta}{\partial X_2} + \mathcal{B}_{1t} \mathcal{R}_G \frac{\partial S_{12}}{\partial t} \right]_{X_2=0} + {O}(\epsilon^3), \end{align}

(4.11b)\begin{align} U_2|_{X_2=0} &={-}\epsilon^2 \left[ m_{12}\frac{\partial S_{12}}{\partial X_1}+ m_{32}\frac{\partial S_{32}}{\partial X_3}\right]_{X_2=0} + {O}(\epsilon^3), \end{align}

(4.11c)\begin{align} U_3|_{X_2=0} &= \epsilon \, \lambda_z S_{32}|_{X_2=0} + \epsilon^2 \left[ m_{32} \frac{\partial S_{22}}{\partial X_3}+ \mathcal{B}_{3t} \mathcal{R}_G \frac{\partial S_{32}}{\partial t}\right]_{X_2=0}+ {O}(\epsilon^3), \end{align}

(4.11d)\begin{align} \varTheta|_{X_2=0} &= 1+ \epsilon \, \lambda_z \left.\frac{\partial \varTheta}{\partial X_2}\right|_{X_2=0} + \epsilon^2 \mathcal{B}_{3t}\mathcal{R}_G\, Pr \left.\frac{\partial^2 \varTheta}{\partial X_2 \partial t}\right|_{X_2=0} + {O}(\epsilon^3). \end{align}

The no-slip conditions of the smooth surface are identically retrieved at ${O} (\epsilon ^0)$. The effective conditions for velocity are similar to those given by Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2020) and Bottaro & Naqvi (Reference Bottaro and Naqvi2020) for flow over rough surfaces without heat transfer. Nevertheless, the presence of the buoyancy terms, proportional to $\mathcal {R}_G$ and $\mathcal {R}_G ({\partial \varTheta }/{\partial X_2})$ in the equation of the velocity component $U_1$, and of the time-fluctuation terms in the equations of ($U_1, U_3, \varTheta$) should be highlighted. We emphasize that the presence of the buoyancy-related term is a first-order contribution to the effective condition for the streamwise velocity, $\hat U_1$, and is directly attributed to the assumption that the Grashof number is sufficiently large, i.e. $\epsilon \mathcal {R}_G = \epsilon ^2 \sqrt {Gr}$ is of ${O} (\epsilon ^1)$, and not ${O} (\epsilon ^2)$.

In dimensional terms, the conditions on the plane $\hat x_2 = 0$ read

(4.12a)\begin{align} \hat u_1|_{\hat x_2 = 0} &\approx \underbrace{ \hat \lambda_{x} \left[\frac{\partial \hat u_1}{\partial \hat x_2}+\frac{\partial \hat u_2}{\partial \hat x_1}\right]_{\hat x_2 = 0} + \hat m_{12}\frac{\beta g (\hat T_w- \hat T_\infty)}{\nu}}_{\text{First order}} \nonumber\\ &\quad + \underbrace{\hat m_{12} \left[ \frac{\partial}{\partial \hat x_1}\left( - \frac{(\hat p - \hat p_\infty)}{\mu}+2\frac{\partial \hat u_2}{\partial \hat x_2}\right)\right]_{\hat x_2 = 0}+ \hat{\mathcal{B}} \frac{g \beta}{\nu} \left.\frac{\partial \hat T}{\partial \hat x_2}\right|_{\hat x_2 = 0}}_{\text{Second order}}\nonumber\\ &\quad + \underbrace {\hat{\mathcal{B}}_{1t} \frac{1}{\nu} \frac{\partial}{\partial \hat t}\left[\frac{\partial \hat u_1}{\partial \hat x_2}+\frac{\partial \hat u_2}{\partial \hat x_1}\right]_{\hat x_2 = 0}}_{\text{Second order}}, \end{align}

(4.12b)\begin{align} \hat u_2|_{\hat x_2 = 0} & \approx \underbrace{ -\hat m_{12} \frac{\partial}{\partial \hat x_1}\left[\frac{\partial \hat u_1}{\partial \hat x_2}+\frac{\partial \hat u_2}{\partial \hat x_1}\right]_{\hat x_2 = 0}-\hat m_{32} \frac{\partial}{\partial \hat x_3}\left[\frac{\partial \hat u_3}{\partial \hat x_2}+\frac{\partial \hat u_2}{\partial \hat x_3}\right]_{\hat x_2 = 0}}_{\text{Second order}}, \end{align}

(4.12c)\begin{align} \hat u_3|_{\hat x_2 = 0} &\approx \underbrace{\hat \lambda_{z} \left[\frac{\partial \hat u_3}{\partial \hat x_2}+\frac{\partial \hat u_2}{\partial \hat x_3}\right]_{\hat x_2 = 0}}_{\text{First order}}+ \underbrace{\hat m_{32} \,\frac{\partial}{\partial \hat x_3}\left[ - \frac{(\hat p - \hat p_\infty)}{\mu}+2\frac{\partial \hat u_2}{\partial \hat x_2}\right]_{\hat x_2 = 0}}_{\text{Second order}} \nonumber\\ &\quad +\underbrace{\hat{\mathcal{B}}_{3t} \frac{1}{\nu} \frac{\partial}{\partial \hat t}\left[\frac{\partial \hat u_3}{\partial \hat x_2}+\frac{\partial \hat u_2}{\partial \hat x_3}\right]_{\hat x_2 = 0}}_{\text{Second order}}, \end{align}

(4.12d)\begin{align} \hat T|_{\hat x_2 = 0} &\approx \underbrace{\hat T_w}_{\text{Zero order}} + \underbrace{\hat \lambda_z \left.\frac{\partial \hat T}{\partial \hat x_2}\right|_{\hat x_2 = 0}}_{\text{First order}} + \underbrace{\hat{\mathcal{B}}_{3t} \left.\frac{1}{\alpha} \frac{\partial^2 \hat T}{\partial \hat x_2 \partial \hat t}\right|_{\hat x_2 = 0}}_{\text{Second order}}. \end{align}

The dimensional groups of coefficients $(\hat \lambda _x, \hat \lambda _z)$, $(\hat m_{12}, \hat m_{32})$ and ($\hat {\mathcal {B}}, \hat {\mathcal {B}}_{1t}, \hat {\mathcal {B}}_{3t}$) are homogeneous to, respectively, a length, a surface area and a volume, and correspond to the product of their dimensionless counterparts times, respectively, $l$, $l^2$ and $l^3$. The conditions above represent the most important contribution of the present paper.

6.1. Isothermal vertical smooth surface case

A numerical calculation is first performed for a smooth isothermal surface at a plate Grashof number $Gr= 5.563 \times 10^8$ and a Prandtl number $Pr= 0.712$; this corresponds, for instance, to a buoyancy-driven air flow with $\hat T_\infty =18 \,^\circ \textrm {C}$, $\hat T_w=58 \,^\circ \textrm {C}$, $L= 0.5$ m and the fluid properties calculated at standard pressure and based on the film temperature $\hat T_f={(\hat T_w+\hat T_\infty )}/{2}$. Different purposes are targeted from this step: (i) estimation of the adequacy of the computational domain; (ii) validation of the CFD numerical scheme and of the inlet/outlet boundary conditions by comparing the results with available databases through the literature; (iii) the no-slip smooth surface case is equivalent to a homogenized simulation of the rough surface with zero-order effective conditions, so the results will help to monitor the accuracy enhancement when progressively higher-order approximations are used.

The computational domain and the boundary conditions are illustrated in figure 6. No-slip and constant temperature conditions are defined on the vertical wall; uniform pressure boundary conditions are imposed at the upper and the lower boundaries such that an equilibrium with the hydrostatic pressure head is satisfied. The width of the domain should be selected in such a way that the streamwise velocity smoothly vanishes at the far boundary at $X_2=S$, with the normal gradients of horizontal velocity and temperature smoothly decreasing to zero. This was checked by running the simulation with different values of the domain width, $S$, and monitoring a result of interest (the surface-averaged Nusselt number) until convergence was attained. The local Nusselt number ($Nu$) and its surface-averaged counterpart ($\overline {Nu}$) are defined for the smooth surface by

(6.1a)\begin{gather} Nu= \frac{-L}{\hat T_w-\hat T_\infty} \left.\frac{\partial \hat T}{\partial \hat x_2}\right|_{X_2=0}={-}\left.\frac{\partial \varTheta}{\partial X_2}\right|_{X_2=0}, \end{gather}

(6.1b)\begin{gather}\overline{Nu} =\int_0^1 -\left.\frac{\partial \varTheta}{\partial X_2}\right|_{X_2=0}\,{\rm d}X_1 . \end{gather}

As can be realized from figure 6(b), a domain width $S= 0.8$ appears to be sufficient; however, a value of $S= 2$ was used throughout the work to ensure the absence of spurious reflections from the outer boundary when testing microstructured walls and/or larger values of $Gr$. The two-dimensional grid is described in detail in Appendix C; eventually, the extrapolated value of the average Nusselt number is estimated to be 75.055 based on the conducted mesh-dependency study, also illustrated in the Appendix.

Figure 6. Computational domain with boundary conditions used for the numerical simulation of the natural convection over the isothermal vertical smooth plate ($X_2=0$). Panel (b) shows a graphical representation of the sensitivity of the solution to variations in the domain width. Here, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

The dimensionless temperature and streamwise velocity profiles are plotted across chosen normal sections distributed along the plate, as displayed in figure 7. The velocity and the temperature contours in the vicinity of the smooth wall are also shown, to highlight the development of the boundary layers. The peak of the velocity profile shifts away from the wall as $X_1$ increases, in qualitative agreement with the estimate of the classical Squire–Eckert theory (Lienhard & Lienhard Reference Lienhard IV and Lienhard V2019) according to which the velocity peaks at almost $\frac {1}{3}$ of the boundary layer thickness. At the same time, the temperature gradient at the wall is reduced with $X_1$. The latter effect is responsible for the decrease of the local Nusselt number ($Nu$) along the plate, plotted in figure 8. The distribution of the local Nusselt number is in perfect agreement with the corresponding reference results by Ostrach (Reference Ostrach1953). An analysis of Ostrach's results reveals that the Nusselt number ($Nu$) is related to the vertical position ($X_1$) via the expression

(6.2)\begin{equation} \frac{X_1 Nu}{\left(\dfrac{Gr}{4} X_1^3\right)^{0.25}} = f_n(Pr). \end{equation}

At a Prandtl number of 0.712, the function $f_n(Pr)$ was estimated to be almost 0.504. Therefore, (6.2) can be recast as an explicit relation between $Nu$ and $X_1$ at any fixed value of the Grashof number.

Figure 7. Smooth surface case: dimensionless velocity and temperature profiles across different normal sections distributed along the vertical plate. Contours representing the velocity and the thermal fields are also provided. Here, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

Figure 8. Smooth surface case: numerical prediction of the local Nusselt number distribution along the vertical plate, compared with the expected behaviour based on the similarity solution by Ostrach (Reference Ostrach1953).

6.2. The case of isothermal ribbed surface

A typical validation case is now considered. The developed asymptotic wall model is assumed to be reasonably accurate provided that $\epsilon$ is sufficiently small. In addition, limitations are imposed on the magnitude of the coefficient of the convective term in the normalized microscopic governing equations, $C= \epsilon ^2\, \sqrt {Gr}= \epsilon \, \mathcal {R}_G$, for convective effects to be absent in the leading-order problem but present at next order. For the basic validation case, we consider natural convection over an isothermal vertical plate with 168 transverse square ribs $(\epsilon = \frac {1}{168})$ with a pitch distance to rib height ratio ${l}/{e}=3.75$. The problem is characterized by a plate Grashof number $Gr= 5.563 \times 10^8$ and a Prandtl number $Pr= 0.712$. With these parameters, the value of the coefficient $C$ is 0.836. Results of the feature-resolving simulation and the homogenization-based calculations of the basic ribbed surface case are presented and compared.

6.2.1. Feature-resolving simulation of the ribbed surface case

The two-dimensional feature-resolving numerical simulation, where the details of the ribbed surface are captured by the grid, represents a necessary step for the validation of the homogenized model.

The computational domain is illustrated in figure 9, including the geometric details of the ribbed surface. The applied boundary conditions are the same as in the smooth surface case, taking into account that the no-slip velocity and temperature conditions are now imposed on a patterned surface, not on a plain one. The two-dimensional grid near the ribs is also shown, and the different grid refinement levels are stated. A near-wall region of thickness $5e$ is defined where a high mesh density is employed to capture the flow dynamics in the vicinity of the perturbed surface; however, the gradual growth of the mesh guarantees that the whole field is fairly well resolved. The number of two-dimensional cells given in the figure illustrates clearly the high computational cost of the fully featured simulation of the ribbed surface compared with requirements of the smooth surface case, described in Appendix C. The Nusselt number at any point on the ribbed surface is given by

(6.3a)\begin{equation} Nu= \frac{-L}{\hat T_w-\hat T_\infty}\times \left.\frac{\partial \hat T}{\partial \hat n}\right|_{wall}={-}\left.\frac{\partial \varTheta}{\partial n}\right|_{wall}, \end{equation}

where $\hat n$ denotes the dimensional distance in the surface-normal direction and $n={\hat n}/{L}$. A dimensional surface distance $\hat s$ is defined in such a way that it goes along the ribbed surface capturing its details, i.e. $\hat s$ goes from 0 to $L+ (2e \times N_{ribs})$ with $N_{ribs}={1}/{\epsilon }= {L}/{\ell }$ the number of ribs. Accordingly, the surface-averaged Nusselt number based on the projected area of the two-dimensional ribbed plate is defined as

(6.3b)\begin{equation} \overline{Nu} = \frac{1}{L}\int_0^{L+{2e}/{\epsilon}}-\left.\frac{\partial \varTheta}{\partial n}\right|_{wall}\,{\rm d} \hat s = \int_0^{1+{2e}/{\ell}}-\left.\frac{\partial \varTheta}{\partial n}\right|_{wall}\,{\rm d}s , \end{equation}

where $s={\hat s}/{L}$, and the value of ${e}/{\ell }$ represents the rib height to the pitch distance ratio. The given expression for $\overline {Nu}$ takes into account the surface area increase, with respect to the baseplate area, due to the presence of ribs. For the considered values of parameters ($Gr= 5.563 \times 10^8$, $Pr= 0.712$, $\epsilon = \frac {1}{168}$, ${l}/{e}=3.75$), the reported value of $\overline {Nu}$ was estimated based on Richardson's extrapolation of results for successively refined grids, and finally found to be 73.200 (compared with a value of $75.055$ for a corresponding smooth surface case). This finding suggests that adding ribs to the vertical surface deteriorates the total heat transfer rate, for the geometric parameters and flow conditions under study.

Figure 9. Description of the computational domain and two-dimensional grid structure used for the feature-resolving numerical simulation. The grid is shown for one unit in the near-wall region (of thickness $5e$), defined for the highest mesh density. Indicated mesh parameters for different refinement levels are: $N_{unit}$: number of cells for a single unit in the near-wall region; $N_{int.}$: number of cells in the whole near-wall region of thickness $5e$; $N_{total}$: total number of cells in the computational domain.

The fully featured simulation is described first to provide insight into the physics, before turning to the homogenized model. The patterns of the streamwise velocity, the normal velocity and the temperature are plotted over two distant regions along the plane surface tangent to the outer rims of the square ribs in order to show the behaviour of velocity and thermal fields near the leading edge and near the top of the plate, as displayed in figure 10. The fictitious surface at $X_2=0$ was specifically chosen for the plots as it represents the plane on which the effective conditions are imposed in the model simulations; therefore, monitoring the flow parameters along this surface is of interest. The contours of the velocity and the temperature near the wall are also shown so that details of the boundary layer can be captured. Velocity and temperature patterns are perturbed by the presence of the ribs and experience quasi-periodic behaviours along the vertical distance. By analysing one unit of the distributions shown in the plots, it is evident that the no-slip velocity and temperature conditions are typically satisfied at the physical surface of the rib whereas deviations occur in the inter-rib fluid region. Proceeding along the vertical direction, the average levels of both the streamwise velocity and the temperature increase, which is qualitatively similar to the smooth surface case. The deflections of the streamlines, due to the flow interaction with the surface protrusions, are directly reflected in perturbation of the normal velocity where the successive negative and positive fluctuations represent, respectively, the inward and outward normal flow through the inter-rib region. The characteristics of the flow structure and the way in which the heat transfer from the surface is accordingly affected are shown in figure 11. The flow behaviour close to the ribbed surface is visualized with the aid of streamlines in two distant regions along the vertical direction, so that the development of the flow can be monitored. Two distinct flow regimes are observed, a separation–reattachment–separation (SRS) regime and a full separation (FS) regime. For both patterns, the inter-rib region is characterized by the existence of two co-rotating vortices. At relatively low values of the local Grashof number $Gr_x=Gr\, X^3$, i.e. near the leading edge of the plate, the SRS flow regime is present where the low inertia of the mainstream allows the fluid to easily deflect in the normal direction and reattach to the surface of the baseplate, keeping the two eddies well isolated. In contrast, sufficiently away from the leading edge, the FS regime ensues as the increasing inertia of the accelerated stream hinders the normal deflection towards the baseplate, preventing the reattachment of the mainstream. As illustrated in the figure, the two vortices remain connected to each other via an outer belt-like stream that rotates in the same direction of both eddies, representing a separated entity that isolates the main flow from the baseplate in the inter-rib region.

Figure 10. Feature-resolving simulation of the ribbed surface: detailed behaviours of (a) streamwise velocity; (b) normal velocity; (c) temperature. Results are plotted along the vertical plane passing through the outer rims of the ribs, through two specific ranges of $X_1$. The contours of the fields are also provided. Here, ${Gr_L= 5.563 \times 10^8}$, $Pr= 0.712$.

Figure 11. Feature-resolving simulation of the ribbed surface: comparative description of the flow pattern and the behaviour of the Nusselt number through two distant regions: (a) near the leading edge; (c) near the end of the plate. The contours of $U_1$ are also shown; the colour map given in figure 10 is modified here so that the white portions within the grooves represent the regions with negative streamwise velocity, i.e. the backflow regions. Here, $\epsilon = \frac {1}{168}$, ${l}/{e}= 3.75$, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

The associated heat transfer behaviour is plotted in figure 11 in terms of detailed patterns of the local Nusselt number $Nu$. A quasi-periodic behaviour of the Nusselt number is observed while proceeding along the vertical plate, similarly to literature observations (Bhavnani & Bergles Reference Bhavnani and Bergles1990; Tanda Reference Tanda1997, Reference Tanda2008, Reference Tanda2017; Nishikawa et al. Reference Nishikawa, Otomo, Yoshida, Deguchi, Tsukamoto and Yamamoto2020). On a single-unit scale of analysis, it is evident that the heat transfer rate drastically drops just upstream and downstream of the square protrusion, a fact ascribed to the presence of the separation eddies that form a hot inactive zone in the vicinity of the rib where the thermal boundary layer thickening mitigates the heat transfer process. Conversely, the local Nusselt number peaks at some location within the inter-rib region as the mainstream reattaches to the surface of the baseplate. Even in the FS regime, the inter-rib peak is experienced since the mainstream still approaches the surface (without reattaching). The major peak of the local Nusselt number is present on the outer rim of the rib due to the considerable local thinning of the thermal boundary layer. From a macroscopic point of view, the average value of $Nu$ decreases away from the leading edge along with the development of the thermal boundary layer.

6.3. The macroscopic homogenization-based simulations

The effect of the surface microstructure on the behaviour of the buoyancy-driven stream is replaced here by the implementation of the homogenized effective boundary conditions on the plane at $X_2=0$ (refer to figure 12). As the present work targets the validation of the model on the steady-state solution of a two-dimensional laminar flow, the effective conditions can be simplified by neglecting the time-derivative terms and the gradients in the spanwise direction. The dimensionless conditions up to the second order in all variables thus read

(6.4a)\begin{gather} U_1|_{X_2=0} = \epsilon \, [\lambda_x S_{12}+ m_{12} \mathcal{R}_G ]_{X_2=0} + \epsilon^2 \left[m_{12}\frac{\partial S_{22}}{\partial X_1}+ \mathcal{B} \mathcal{R}_G \frac{\partial \varTheta}{\partial X_2} \right]_{X_2=0}+{O}(\epsilon^3), \end{gather}

(6.4b)\begin{gather}U_2|_{X_2=0} ={-}\epsilon^2 \left[ m_{12}\frac{\partial S_{12}}{\partial X_1}\right]_{X_2=0}+{O}(\epsilon^3), \end{gather}

(6.4c)\begin{gather}\varTheta|_{X_2=0} = 1+ \epsilon \, \lambda_z \left.\frac{\partial \varTheta}{\partial X_2}\right|_{X_2=0} + {O}(\epsilon^3), \end{gather}

Based on the parametric study presented in § 5, at ${l}/{e}=3.75$, the following values of the model coefficients are found:

(6.5a–d)\begin{equation} \lambda_x= 0.03791, \quad \lambda_z= 0.08404, \quad m_{12}= 0.002125, \quad \mathcal{B}= 0.0002247. \end{equation}

Since the ribbed surface is impermeable, the transpiration velocity is zero on average and its inclusion is not significant under laminar flow conditions; this was tested and confirmed in the present work.

Figure 12. Effective boundary conditions for streamwise velocity and temperature in comparison with the running-average results of the feature-resolving simulation. Here, $\epsilon = \frac {1}{168}, {l}/{e}= 3.75$, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

The set-up of the homogenization-based macroscopic simulations is similar to the set-up of the smooth surface case with regard to the computational domain, the grid structure, the refinement levels and the boundary conditions, except for replacing the no-slip velocity and temperature conditions by the effective conditions ((6.4a) to (6.4c)) on a virtual wall in $X_2=0$. It is comforting that the macroscopic simulations reach mesh independence for grids which are more than 30 times coarser as compared with the fully featured case, while providing accurate predictions of the surface-averaged Nusselt number (the metric being evaluated in the grid-dependence study). For the considered flow and geometric conditions ($Gr= 5.563 \times 10^8$, $Pr= 0.712$, $\epsilon = \frac {1}{168}$, ${l}/{e}=3.75$), the converged values of $\overline {Nu}$ with first-order and second-order conditions are, respectively, $73.1667$ and $73.1618$. In comparison with the fully featured result, the errors of the homogenized models are, respectively, $-$0.045% and $-$0.052%. It is worth restating that $\overline {Nu}$ of the smooth surface case is 2.54$\%$ larger than in the fully featured ribbed case.

The results which can be achieved from the homogenized simulations illustrate the macroscopic behaviour of velocity and temperature fields; clearly, these results should be interpreted as spatially averaged values, whereas the detailed patterns near the wall are unavailable from the model simulations. For this purpose, the validation of the present approach is done by comparing the results of the macroscopic simulations with the running-average values of the fully featured fields over streamwise distances equal to the periodicity of the pattern of the surface structure. For instance, the running-average value of the dimensionless velocity $U_1$ at an arbitrary point ($X_1=a, X_2=b$) is computed as

(6.6)\begin{equation} \langle U_1\rangle |_{X_1=a,X_2=b} = \frac{1}{\epsilon} \int_{a-{\epsilon}/{2}}^{a+{\epsilon}/{2}}U_1(X_1,b)\,{\rm d}{X_1}. \end{equation}

The numerical predictions of $U_1$ and $\varTheta$ resulting from the macroscopic simulations with the first-order accurate and the second-order accurate boundary conditions are extracted at the fictitious boundary in $X_2=0$ to explicitly assess the accuracy of the expressions given in (6.4a) to (6.4c). The homogenized results are plotted in figure 12 in comparison with the corresponding running-average values of the feature-resolving simulation. It is clear that the present model can qualitatively predict the difference of the results from the no-slip values. The results show perfect agreement of the effective temperature estimates, apparently insensitive to the mild deviations observed for the predictions of the slip velocity. This fact may be attributed to the absence of strong nonlinearities, i.e. the coupling between the velocity and the thermal fields is weak.

In order to show how the effect of the homogenized conditions propagates from the virtual wall to the flow domain, the profiles of streamwise velocity and temperature are plotted across two normal sections and compared with the corresponding running-average profiles (figures 13 and 14). It is noticeable that, in the present case, the effect of the surface inhomogeneities on the flow field is moderate. Another point is that the predictions based on first- and second-order conditions are almost indistinguishable from one another to graphical accuracy, due to the very small value of $\epsilon$. The normal gradients of $\varTheta$ along the fictitious boundary, represented by the slopes at $X_2=0$ of the $\varTheta$ profiles, were used to obtain the macroscopic behaviour of the Nusselt number along the plate (cf. (6.1a)). The results are presented in comparison with the corresponding running-average values from the fully featured simulation in figure 15. It can be realized that, under the present conditions, the ribs on the surface have a very mildly unfavourable effect on the heat transfer rate. It is very important to highlight that the present approach is only able to model the temperature-gradient-based heat transfer from the matching interface, while the convective contribution, resulting from the product of normal velocity and temperature, is not accounted for, since the fluctuations of the normal velocity cannot be resolved by the homogenized model under laminar flow conditions. The applicability of the model is, therefore, limited here to cases in which convective effects through the fictitious plane are negligible. This is assumed to be valid in the absence of strong nonlinearities that may occur for large values of $\epsilon$ or in the presence of turbulence.

Figure 13. Homogenized-model predictions of the streamwise velocity profiles across two normal sections in comparison with the running-average results of the fully featured simulation. Here, $\epsilon = \frac {1}{168}$, ${l}/{e}= 3.75$, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

Figure 14. Validation of the homogenized-model results of the temperature profiles across two normal sections in comparison with the running-average results extracted from the fully featured simulation. Here, $\epsilon = \frac {1}{168}$, ${l}/{e}= 3.75$, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

Figure 15. Homogenized-model predictions of the Nusselt number in comparison with the running-average results of the feature-resolving simulation, based on the normal temperature gradient along the plane $X_2=0$. Here, $\epsilon = \frac {1}{168}$, ${l}/{e}= 3.75$, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

7.1. Effects of the increase in $\epsilon$ at a given Grashof number

The simulations of the macroscopic problem are now conducted for increasing values of the parameter $\epsilon$ in the effective boundary conditions ((6.4a) to (6.4c)), starting from $\epsilon =\frac {1}{84}$ up to $\epsilon =\frac {1}{10}$, at a constant value of the Grashof number ($Gr= 5.563 \times 10^8$) and for the values of the model coefficients at ${\ell }/{e}= 3.75$ (cf. (6.5a–d)), in order to monitor the deterioration of the model with the increase of the controlling parameters $\epsilon$ and $C$. First, a validation database has been built by running the fully featured simulations with the corresponding numbers of ribs (from 84 to 10). The running-average fields obtained from these simulations are presented in a comparative manner in Appendix E. The results of the macroscopic simulations with first-order and second-order accurate homogenized effective conditions are validated by comparing the streamwise velocity profiles and the temperature profiles across a normal section taken at $X_1=0.5$ with the corresponding running-average patterns from the fully resolved numerical simulations, cf. figures 16 and 17. The purpose is to ascertain the validity range of the asymptotic model away from the conditions ($\epsilon =\frac {1}{168}, C= 0.836$) discussed in § 6. In general, the predictions of the present approach concerning velocity and temperature fields are reliable below $\epsilon =\frac {1}{21}$ at the given Grashof number. It will be argued later that the reliability range becomes wider at lower values of the Grashof number. The accuracy of the temperature predictions is better than the velocity predictions, especially above the mentioned limit where the boundary conditions at second order are able to produce better results in comparison with the first-order conditions.

Figure 16. Predictions of the homogenized model in comparison with the running averages of the fully featured results of the dimensionless streamwise velocity profiles at $X_1=0.5$ for different values of $\epsilon$. Here, ${l}/{e}= 3.75$, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

Figure 17. Predictions of the homogenized model in comparison with the running averages of the fully featured results of the dimensionless temperature profiles at $X_1=0.5$ for different values of $\epsilon$. Here, ${l}/{e}= 3.75$, ${Gr= 5.563 \times 10^8}$, $Pr= 0.712$.

From the practical point of view, the most important factor is the surface-averaged Nusselt number. The behaviour of $\overline {Nu}$ with the increase of $N_{ribs}={1}/{\epsilon }$ is shown in figure 18. Both fully featured and homogenized-model results show a reduced heat transfer for the ribbed surface despite the increase in surface area (compared with a smooth flat plate), for the values of parameters considered here. It is obvious that the level of accuracy of the model predictions is even better than that relative to velocity and temperature profiles. It is also noteworthy that improved predictions of $\overline {Nu}$ by shifting up to the second-order conditions are not systematically guaranteed.

Figure 18. Effect of the number of the ribs ($N_{ribs}={1}/{\epsilon }$) on the surface-averaged Nusselt number (corrected to account for the increase in surface area in the case of feature-resolving simulations of ribbed surfaces, cf. (6.3b)). The parameters are ${l}/{e}= 3.75$, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

The accuracy of the homogenization-based models is reported in a more quantitative manner in table 3. For the velocity and temperature profiles shown in figures 16 and 17, root-mean-square (r.m.s.) deviations between the results of the macroscopic simulations and the results of the reference fully featured simulations are defined. The r.m.s. deviations of the profiles are calculated over a normal distance between $X_2=0$ and $X_2=0.02$. For instance, the r.m.s. deviation of a modelled velocity profile ($U_{mod}$ vs $X_2$) relative to the corresponding fully featured one ($U_{FF}$ vs $X_2$) is defined as

(7.1)\begin{equation} \text{r.m.s. deviation}=\sqrt{\frac{1}{0.02} \int_0^{0.02} \left(\frac{U_{mod}-U_{FF}}{U_{FF}}\right)^2 \,{\rm d}X_2}. \end{equation}

The errors on the predictions of the surface-averaged Nusselt number relative to the fully featured estimations are also shown in the table.

Table 3. Error estimations of the homogenized-model predictions for the velocity and temperature profiles at $X_1=0.5$ and surface-averaged Nusselt number. The fully featured case is used as a reference. The deviations of the results of the smooth-wall case compared with the rough case are also provided. In all cases: ${l}/{e}= 3.75$, $Gr= 5.563 \times 10^8$, $Pr= 0.712$.

7.2. Effect of the Grashof number at a given $\epsilon$

The observed deterioration of the predictions at relatively large values of $\epsilon$ is not explicitly related to the increase in $\epsilon$; rather, it is due to the associated increase of the convective coefficient $C= \epsilon ^2\, \sqrt {Gr}$ beyond a critical limit. In many circumstances (Bottaro Reference Bottaro2019; Bottaro & Naqvi Reference Bottaro and Naqvi2020; Lācis et al. Reference Lācis, Sudhakar, Pasche and Bagheri2020), the theory has been validated for $\epsilon$ up to 0.2. Here, we set $\epsilon =0.1$ and show that by reducing the Grashof number (and thus, $C$), the accuracy of the model improves. The macroscopic simulations are now set at a Grashof number of $7.509 \times 10^6$ (instead of $5.563 \times 10^8$), which results in a decrease of the convective coefficient $C$ from $235.860$ to $27.402$. Figure 19 demonstrates that even at first order, the effective conditions now provide a very good match with fully featured simulation results. The same occurs for the temperature distribution along the virtual interface ($X_2=0$) and the behaviour of the local Nusselt number (figure 20).

Figure 19. Predictions of the homogenized model in comparison with the fully featured running-average results of the dimensionless streamwise velocity and temperature across a normal section at $X_1= 0.5$, for two values of $Gr$: (a,b) $Gr=5.563 \times 10^8$; (c,d) $Gr=7.509 \times 10^6$. For both cases, $\epsilon = \frac {1}{10}$, ${l}/{e}= 3.75$, $Pr= 0.712$.

Figure 20. Predictions for: (a) the temperature distribution along the vertical fictitious surface ($X_2=0$); (b) the performance of the local Nusselt number. Patterns are given at two values of $Gr$. In both cases: $\epsilon = \frac {1}{10}$, ${l}/{e}= 3.75$, $Pr= 0.712$.

7.3. Limit of validity of the approximation

It has been argued in § 7.2 that the accuracy of the proposed homogenization-based model may be linked to a single controlling parameter ($C$) that combines the effects of $\epsilon$ and $Gr$. Therefore, it is advantageous to define a limiting value of $C$ below which the predictions of the presented model are assumed to be reliable. Based on analysis of the accuracy levels shown in table 3, the critical value of $C$ is expected to be around 40; below this value, r.m.s. deviations of the predicted velocity and temperature profiles are, respectively, below $20\,\%$ and $10\,\%$, and the absolute error on the predicted $\overline {Nu}$ is less than $4\,\%$. To validate this estimate, the simulation of the macroscopic problem has been carried out for the case of a vertical surface roughened with only five square ribs, i.e. $\epsilon =0.2$ (relatively large), at $Gr = 9.386\times 10^5$ so that the accuracy of the model at a value of $C=38.752$ can be checked. The geometry of the ribs is characterized by a value of ${\ell }/{e}=3.75$; the model coefficients given in (6.5a–d) are used. The accuracy of the model is assessed through comparative analysis of velocity and temperature predictions across a normal section at $X_1=0.5$ (figure 21). Although the velocity predictions in the near-wall region are not perfect, especially with the first-order conditions, the temperature results are almost identical to the fully featured running-average behaviour. From a practical point of view, the reliability of the thermal field predictions is sufficient to consider the model acceptable under the given condition, i.e. $C \lesssim 40$.

Figure 21. Predictions of the homogenized model in comparison with the fully featured running-average results of the dimensionless temperature profiles at $X_1=0.5$. Here, $\epsilon =0.2$, ${l}/{e}= 3.75$, $Gr= 9.386\times 10^5$, $Pr= 0.712$.

7.4. Model validation on different rib geometries

The robustness of the introduced model is further checked by employing the effective boundary conditions to study natural-convection heat transfer for different configurations of the roughness pattern, particularly those shown in figure 22. Feature-resolving simulations and model calculations are conducted for the natural-convection flow over an isothermal vertical surface roughened with $40$ ribs, at $Gr= 10^8$ and $Pr= 0.712$. Some complex flow structures within the inter-rib regions, captured by the full simulations, are displayed in the figure: these structures highlight the need of very well resolved (and expensive) simulations, highlighting the advantage of implementing equivalent, upscaled boundary conditions. In all cases examined the virtual wall is positioned on a plane passing through the tips/crests/outer rims of the ribs. It is noticeable that the FS regime captured for blunt shapes (geometries B, C, D) differs qualitatively from that displayed in figure 11, with a single separation eddy between neighbouring ribs, isolating the mainstream from the baseplate.

Figure 22. The different roughness patterns considered for model validation: (a) half-circles with a diameter of 0.4$\ell$; (b) rectangular ribs with a width of 0.2$\ell$ and a height of 0.3$\ell$; (c) inward-curved ribs with a width of 0.5$\ell$ and a height of 0.45$\ell$ – circular curvature is tangential to baseplate; (d) right-angle triangles with a width of 0.1$\ell$ and a height of 0.7$\ell$. The microscopic length scale, $\ell$, indicates the pattern periodicity. Streamlines are displayed near and within the wall corrugations, for a specified range of $X_1$.

The model predictions of the surface-averaged Nusselt number ($\overline {Nu}$) are presented in table 4; they are in good agreement with the corresponding results obtained from the full feature-resolving simulations, with a maximum deviation of less than $2.5\,\%$ (detected for configuration A). For the flow conditions under investigation, all rib geometries considered reduce the heat transfer performance of the natural-convection system, with respect to smooth-wall case.

Table 4. Macroscopic coefficients and predictions of $\overline {Nu}$ for the studied rib geometries. the Nusselt number is calculated considering $\epsilon =\frac {1}{40}$, $Gr= 10^8$ and $Pr= 0.712$.