A homogenization approach for buoyancy-induced flows over micro-textured vertical surfaces

Abstract

Asymptotic homogenization is employed to formulate upscaled effective boundary conditions at a smooth virtual surface for a natural-convection flow over a periodically roughened vertical wall, to bypass the expensive numerical resolution of flow and temperature fields near and within wall corrugations. Microscale problems are found by expanding near-wall variables in terms of a small parameter $\epsilon$, the ratio between the microscopic and the macroscopic length scales. The expressions of the upscaled velocity and temperature boundary conditions are provided up to second-order accuracy in $\epsilon$. The case of transverse square ribs is considered as a representative example. The classical Navier-slip condition for the streamwise and the spanwise velocity components is modified at second order by the gradient of the normal stress and the time derivative of the shear stress. The streamwise slip velocity is additionally corrected by a buoyancy term at first order and a temperature-gradient term at second order. The normal velocity at the virtual surface appears only as a second-order transpiration condition. A Robin-like condition for the temperature is found, where the wall temperature is corrected with a temperature-gradient term representing thermal slip. The accuracy levels and the applicability range of the effective conditions to mimic the macroscopic flow behaviour are investigated under laminar flow conditions, in comparison with results of full feature-resolving simulations. A formal validity limit for the approximation is sought in terms of a single accuracy criterion ($C$), which combines the effects of the Grashof number and the ribs’ density. The introduced model is further tested on different rib geometries.

1. Introduction and literature review

Natural convection over ribbed/finned surfaces is widely encountered in engineering applications, such as cooling of electronics and telecommunication devices, air solar collectors and gas-cooled nuclear reactors. Compared with forced convection, a system that depends on the natural-convection heat transfer regime has lower initial and running costs, less noise and vibrations, higher reliability, almost maintenance-free operations and better ability for use in hostile environments under dust, moist air, etc. On the other hand, the main problem facing designers is the low heat transfer coefficient of these systems. Due to the ever-growing trend of miniaturization of electronic components and the increase in power supply, higher heat generation rates per unit volume are encountered (Joshi, Willson & Hazard 1989). This trend has stimulated many investigations to enhance natural-convection cooling systems so that they can be effective at handling operation requirements. One intuitively appealing solution to enhance the heat transfer performance of these systems is to apply some sort of alteration or disturbance on the heated surface(s) in analogy to the well-established concept of heat transfer promotion by adding ribs/fins to surfaces exposed to forced convection (Bunker & Donnellan 2003; Chyu, Oluyede & Moon 2007; Han, Dutta & Ekkad 2012). However, studies on the effectiveness of adding surface alterations (ribs, interrupted fins, dimples, etc.) to vertical plates exposed to natural convection have not led yet to convincing guidelines, with some researchers reporting an improvement of up to 200$\%$ compared with the performance of plane vertical plates, and others who have found them useless or even of negative influence on the local and averaged heat transfer parameters (Bhavnani & Bergles 1990).

The need to better understand the interaction between the surface microstructure and the buoyancy-driven flow has motivated many experimental and numerical investigations to assess the usefulness and the feasibility of adding different types of protrusions to the heated surfaces in terms of their effects on the flow regime, the heat transfer characteristics and the mass of the cooling modules. Examples of some surface alterations/extensions, considered in previous investigations, are displayed in figure 1, including periodic (wavy, rounded, zigzag) corrugations (Kishinami et al. 1990; Bhavnani & Bergles 1991; Yao 2006; Hærvig & Sørensen 2020), steps (Bhavnani & Bergles 1990), two-dimensional ribs (Tanda 1997; Cavazzuti & Corticelli 2008) and different arrangements of fins (Guglielmini et al. 1987; Ahmadi et al. 2014; El Ghandouri et al. 2020). Various experimental techniques have been adopted for mapping the thermal field to assess the detailed heat transfer performance. Two-dimensional and three-dimensional feature-resolving numerical simulations have also demonstrated to be powerful tools for the acquisition of large amounts of data on thermal fields and flow regimes, particularly for complex configurations (Yao 2006; Cavazzuti & Corticelli 2008; Ahmadi et al. 2014; Hærvig & Sørensen 2020).

Figure 1. Schematic drawings of some surface alterations examined in the literature. The indicated geometries are (from top left to bottom right): sinusoidal waves (Bhavnani & Bergles 1991), complex waves (Yao 2006), convex–concave semi-circles (Kishinami, Saito & Tokura 1990), zigzag shaping (Hærvig & Sørensen 2020), steps (Bhavnani & Bergles 1990); transverse square ribs (Tanda 1997), transverse trapezoidal ribs (Cavazzuti & Corticelli 2008), rippled vertical fins (El Ghandouri et al. 2020), staggered arrangement of interrupted fins (Guglielmini, Nannei & Tanda 1987), in-line arrangement of interrupted fins (Ahmadi, Mostafavi & Bahrami 2014).

Numerical work on flows over surfaces with complicated small-scale details including irregularities, roughness, porosity, etc. has been a challenge due to the high computational resources required to numerically discretize flow and temperature fields in the vicinity of the surface microstructures. In the present work, the multiscale homogenization approach is proposed to simplify the modelling of buoyancy-driven flows over periodically roughened vertical surfaces, while maintaining an acceptable level of accuracy. Asymptotic homogenization is an approach which targets the study of the macroscale behaviour of a medium which contains microscopic details, by replacing the rapidly varying properties related to the heterogeneity of the medium by equivalent homogeneous macroscopic properties (Babuška 1976). This technique can play a pivotal role when handling differential equations that govern physical problems with microscale fluctuations (Engquist & Souganidis 2008) which are characterized by some sort of periodicity or pseudo-periodicity. These problems can be computationally simplified by first solving ad hoc auxiliary systems of equations in a microscopic domain to evaluate the necessary upscaled conditions by means of averaging. The approach relies on the asymptotic expansion of the dependent variables in terms of a wisely chosen small parameter whose existence is related to the presence of well-separated scales, for instance a microscopic length scale ($\ell$) and a macroscopic length scale ($L \gg \ell$), so that the parameter $\epsilon = {\ell }/{L} \ll 1$ can be defined, and the solution of the problem can be sought up to different orders of accuracy in terms of $\epsilon$.

Flow over micro-textured surfaces represents a typical homogenization problem. Jiménez Bolaños & Vernescu (2017) have derived the Navier-slip effective condition for the Stokes flow over a rough surface via homogenization theory as a first-order corrector term to the no-slip condition of a smooth surface. Zampogna, Magnaudet & Bottaro (2019a) have pursued a generalization of the classical first-order Navier-slip condition (Navier 1823) over a rough surface by means of a third-order Navier-slip tensor. The homogenized model was pushed to second order by Lācis et al. (2020) with the introduction of a transpiration velocity, the normal velocity component at the fictitious interface, thus enhancing model predictions for a turbulent boundary layer over a rough surface. A further improvement has been added by Bottaro & Naqvi (2020), who sought a solution up to third-order accuracy. The range of applications subtended by homogenization theory is being continuously widened and enhancements to the basic formulation are ongoing. Zampogna et al. (2019b) have extended the theory to the study of the turbulent flow over compliant riblets, seeking reduction of the skin friction drag. Adjoint homogenization has been introduced by Bottaro (2019) as a method to take into account nonlinear effects within the microscopic region.

The work presented in this paper is a novel implementation of the multiscale homogenization technique to study natural-convection heat transfer over rough surfaces. The only previous contribution in this aspect was the work by Introïni, Quintard & Duval (2011), who applied the volume-averaging upscaling method to the study of the steady laminar buoyancy-driven flow over rough surfaces. However, their model suffered from some deficiencies that limit its applicability range. A critical assumption adopted by Introïni and collaborators was the neglect of buoyancy effects within the microscopic region, so that the momentum and energy conservation equations are decoupled. This assumption, despite being mathematically advantageous, limits the model applicability to cases in which the Rayleigh number characterizing the microscopic problem (based on the microscopic length scale and the temperature difference across the microscopic region) is sufficiently small. To satisfy this condition, the bulk Rayleigh number must be lower than some threshold value, and the roughness elements must be confined within the thermal boundary layer. In practical situations, high values of the Rayleigh number are often encountered. Moreover, the model developed by Introïni et al. (2011) is only accurate to first order in $\epsilon$.

In this paper, asymptotic homogenization is used to formulate expressions for the macroscopic velocity and temperature effective conditions at a virtual interface separating the microscopic and the macroscopic sub-domains. Unlike Introïni et al. (2011), the Boussinesq approximation is employed for the buoyancy term in the microscopic momentum equation to be linearly coupled with the energy equation. The dependent parameters are expanded asymptotically in powers of the small parameter $\epsilon ={\text {pattern periodicity }(\ell )}/{\text {plate length }(L)}$. The effective conditions for velocity and temperature are all sought up to second-order accuracy. In the next section, the governing equations and the boundary conditions of the problem are outlined, and domain decomposition is explained. In § 3, the microscopic region is considered where the asymptotic expansion of the dependent variables is defined, and the problem is reconstructed at different orders of $\epsilon$. For each order, generic forms of the solutions are assumed and auxiliary differential systems are formulated. Then, the case of transverse square ribs is discussed in § 4. The parameters of interest are determined via numerical solution of the auxiliary systems, and the effect of the matching surface location is considered. In § 5, a parametric study seeking the effect of varying the rib size to the pitch distance ratio on the different coefficients is presented. In § 6, the macroscale problem is considered by imposing the upscaled boundary conditions at a virtual vertical interface passing through the outer rims of the ribs; full feature-resolving simulations are also conducted to validate the predictions of the model. In § 7, the accuracy deterioration of the homogenized model is monitored with the increase of $\epsilon$ and/or the Grashof number, and the limit of validity of the approach is ascertained. Furthermore, the accuracy of the method is confirmed for different shapes of the roughness elements. In the concluding section, the main findings of the study are highlighted.

2. Governing equations and domain decomposition

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. Asymptotic analysis of the microscale problem

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 (2020), 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 1989; Luchini, Manzo & Pozzi 1991).

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. The case of transverse square ribs

As an example of the implementation of the theory, the case of transverse square ribs is considered so that the auxiliary systems can be significantly simplified. In particular, because of invariance along $x_3$, all auxiliary problems simplify considerably (with derivatives $\partial /\partial x_3$ set to zero), and only two-dimensional Stokes-like (or Laplace-like or Poisson-like) problems remain to be solved in the $(x_1, x_2)$ plane, subject to periodic conditions along $x_1$. A sketch of the microscopic representative volume element is provided in figure 3. Some of the microscopic problems admit trivial solutions. For instance, it is easy to find that in the elementary cell it is $\breve u_{12}=\breve u_{22}=\breve u_{13}=\breve u_{23}=\breve u_{31}=\breve u_{32}=u_3^{\dagger} =0$, plus $\breve p_2 = -1$ and $\breve p_3=0$. The systems which do not have a simple solution have been solved numerically by using the STAR-CCM+ multi-physics software (version 15.06.007-R8), by successfully refining the grid until fully grid-converged states are found, for varying dimensions of the cell along $x_2$. Detailed numerical results of the reduced auxiliary systems relative to the ${O} (\epsilon ^0)$, ${O} (\epsilon ^1)$ and ${O} (\epsilon ^2)$ problems are presented as supplementary material available at https://doi.org/10.1017/jfm.2022.320, for a rib size to periodicity ratio, $e/l$, equal to $0.25$, and matching interface location positioned at $x_2=y_\infty = 5$.

Figure 3. Sketch of a unit cell in the microscopic domain, indicating coordinates and geometric parameters.

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 (2020), 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. (2020) and Bottaro & Naqvi (2020) 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.

5. The role of rib height to pitch distance ratio: parametric study

From a practical point of view, it is advantageous to generate a database of the values of the seven dimensionless, geometry-dependent coefficients of interest, to cover a wide range of rib height to pitch distance ratios, $e/l$, in order to enable the direct use of the effective boundary conditions for the macroscopic problems. In this study, the ratio was varied within the range $0.025 \leq e/l \leq 0.8$. For each value of $e/l$, the procedure described in § 4.2 for the accurate estimation of the coefficients was followed. The resulting database is presented in tabular form (table 2) and graphically in figure 5. It is clear that all model coefficients peak, in magnitude, within the range $e/\ell =0.1$ to $0.3$, which implies significant velocity and thermal slip. All coefficients tend to zero as $e$ tends to zero or approaches $\ell$, for the effective boundary conditions at $x_2=0$ to become no slip and isothermal wall.

Figure 5. The behaviour of the upscaled coefficients of interest with the change of rib height to pitch distance ratio. The curves are fitted on the basis of kriging interpolation.

Table 2. The upscaled coefficients of interest for different rib height to pitch distance ratios.

6. Macroscale behaviour of the flow

In this section, attention is given to validation of the effective conditions obtained in § 4, with the upscaled coefficients calculated for the case of square ribs. The macroscale problem is considered, with the governing equations given in § 2. Since the ribs are elongated in the transverse direction, and since only the case of laminar flow is considered, there is no need to resolve the spanwise direction; the problem can be simplified to its two-dimensional form in the ($X_1,X_2$) plane. In addition, steady-state solutions are targeted for validation purposes. Three types of simulations have been carried out: (i) natural convection over a vertical smooth surface; (ii) full feature-resolving natural convection over a vertical ribbed surface; (iii) homogenized problem with effective boundary conditions at a virtual wall. For each simulation, the computational domain, the boundary conditions and the grid structure are explained in detail later in this section. As for the case of the microscopic problems, we have found it convenient to carry out the simulations with STAR-CCM+. The second-order upwind formulation has been adopted for the spatial discretization of all fields, with the calculation of the gradients based on a hybrid Gauss–least squares method. The SIMPLE scheme has been employed for the pressure–velocity coupling.

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 2019) 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 (1953). 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 (1953).

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 1990; Tanda 1997, 2008, 2017; Nishikawa et al. 2020). 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. Applicability range and limit of validity of the model

In this section, the results of several numerical simulations are presented to assess the deterioration of the accuracy of the proposed technique with the increase of the small parameter $\epsilon = {\ell }/{L}= {1}/{N_{ribs}}$ and the coefficient of the microscopic momentum-convective term $C= \epsilon ^2\, \sqrt {Gr}= \epsilon \, \mathcal {R}_G$.

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 2019; Bottaro & Naqvi 2020; Lācis et al. 2020), 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$.

8. Conclusions

A homogenization-based model is proposed for the study of the heat transfer by free convection over regularly microstructured vertical surfaces. The approach provides a computationally cheap alternative to the standard feature-resolving simulations in the cases where the macroscopic behaviour of the flow is of interest, and it has been adopted in the past for the case of rough, micro-textured surfaces, in the absence of thermal effects. The procedure, eventually, yields parameters needed to enforce equivalent velocity and temperature boundary conditions at a plane virtual surface, up to second order in terms of a small parameter $\epsilon$, the ratio of the pattern periodicity, $\ell$, to the total plate length, $L$. Importantly, the effective boundary conditions derived here do not contain any empirical parameters.

As a typical implementation of the theory, the model is applied to the case of two-dimensional square ribs. The auxiliary systems are then reduced to either two-dimensional Stokes-like problems or Laplace-like or Poisson-like problems, which either admit trivial solutions or require a numerical solution in a periodic representative cell of the microscopic domain. The parameters contributing to the effective conditions belong to seven independent groups, i.e. the numerical solution of only seven auxiliary problems is sufficient to completely retrieve the effective conditions. The results are then extrapolated from distant matching surfaces to the plane passing through the outer edges of the ribs, beyond which the macroscopic simulation is intended to be performed. The most significant finding of the procedure is the proposed form of the effective boundary conditions. For the streamwise slip velocity, a buoyancy term acts as a corrector to the classical Navier-slip condition at first order, while pressure-gradient, temperature-gradient and time-derivative terms appear at second order. A Robin boundary condition appears for the temperature effective condition, where a normal temperature-gradient term, with a coefficient identical to Navier's spanwise slip coefficient, corrects the uniform wall temperature. The spanwise slip velocity and the transpiration velocity are also considered, to allow for example the usage of the model in turbulent flow cases where the spanwise and the normal velocity fluctuations are to be resolved in direct or large-eddy numerical simulations (Bottaro 2019; Lācis et al. 2020). A parametric study is conducted to investigate the effect of varying the rib size to pitch distance ratio on the values of the coefficients.

The efficiency of the proposed first- and second-order accurate conditions in modelling the effect of the surface microstructure on the macroscopic behaviour of the flow has been tested by comparing the obtained thermal and velocity fields with the corresponding results of full feature-resolving simulations at different values of $\epsilon$ and the Grashof number; the case of tiny square ribs is first considered for validation purposes, while other geometries are studied at a later stage for accuracy confirmation. All simulations have been conducted for laminar flow conditions at a constant Prandtl number equal to $0.712$ (air). It is shown that the expensive mesh requirements for resolving complex inter-rib flow structures, associated with the SRS regime at low values of $Gr_x$ and the FS regime at high values of $Gr_x$, can be significantly alleviated when the model is employed. A significant result is that the accuracy of the model can be linked to the single parameter $C= \epsilon ^2\, \sqrt {Gr}$ which measures the significance of the energy flux within the microscopic domain. A value of $C\approx 40$ is the critical limit below which the model is believed to yield acceptable predictions.

The dependence of the accuracy of the proposed model on a single parameter combining the effects of $\epsilon$ and $Gr$ renders the approach applicable to large values of the Grashof number, provided that $\epsilon$ is sufficiently small, i.e. the number of ribs is adequately large. The upscaling model described in this work represents a more versatile version of the effective conditions for natural convection over ribbed surfaces in comparison with the earlier model by Introïni et al. (2011) which neglected the buoyancy effect in the microscopic region and reported a single validity-limiting value of $Gr=10^7$. In addition, the asymptotic homogenization method employed here represents a rigorous tool to formally advance in the order of accuracy. Second-order accurate boundary conditions are attained, which provides an enhancement to the validity range of Introïni's first-order approach.

This work opens up several perspectives, related for example to the accuracy and applicability limit of the model in the case of turbulent natural convection over ribbed surfaces. It would also be interesting to develop an optimization strategy to find optimal wall micro-patterns, able to maximize heat transfer from the surface. The procedure described can be easily extended to the case of weakly conducting or adiabatic corrugation elements. This will constitute the object of future investigations.