Thermal convection over fractal surfaces

Abstract

We use well resolved numerical simulations with the lattice Boltzmann method to study Rayleigh–Bénard convection in cells with a fractal boundary in two dimensions for $Pr = 1$ and $Ra \in \left [10^7, 10^{10}\right ]$, where Pr and Ra are the Prandtl and Rayleigh numbers. The fractal boundaries are functions characterized by power spectral densities $S(k)$ that decay with wavenumber, $k$, as $S(k) \sim k^{p}$ ($p < 0$). The degree of roughness is quantified by the exponent $p$ with $p < -3$ for smooth (differentiable) surfaces and $-3 \le p < -1$ for rough surfaces with Hausdorff dimension $D_f=\frac {1}{2}(p+5)$. By computing the exponent $\beta$ using power law fits of $Nu \sim Ra^{\beta }$, where $Nu$ is the Nusselt number, we find that the heat transport scaling increases with roughness through the top two decades of $Ra \in \left [10^8, 10^{10}\right ]$. For $p$$= -3.0$, $-2.0$ and $-1.5$ we find $\beta = 0.288 \pm 0.005, 0.329 \pm 0.006$ and $0.352 \pm 0.011$, respectively. We also find that the Reynolds number, $Re$, scales as $Re \sim Ra^{\xi }$, where $\xi \approx 0.57$ over $Ra \in \left [10^7, 10^{10}\right ]$, for all $p$ used in the study. For a given value of $p$, the averaged $Nu$ and $Re$ are insensitive to the specific realization of the roughness.

1. Introduction

Thermal convection refers to fluid flows that are driven by buoyancy forces due to density variations, which in turn are caused by gradients in temperature (Chandrasekhar 2013). Such flows are ubiquitous in both the natural and engineering environments, and are key to understanding transport phenomena in the atmospheric boundary layer, in the outer core of Earth and in the outer layers of stars (Kadanoff 2001; Wettlaufer 2011) to name a few examples. The simplest setting in which thermal convection can be studied is classical Rayleigh–Bénard convection (RBC) in which a fluid is confined between two flat horizontal plates with the under side maintained at a higher temperature than the top (Rayleigh 1916). Applying the Boussinesq approximation to the Navier–Stokes equations, the dynamics of RBC is governed by three non-dimensional parameters: the Rayleigh number $Ra$, the ratio of buoyancy to viscous forces, the Prandtl number $Pr$, the ratio of the fluid's kinematic viscosity to its thermal diffusivity, and the aspect ratio $\varGamma$ of the flow domain.

Heat transport in a fluid at rest is due solely to thermal conduction and when convective motions ensue this transport is enhanced. The Nusselt number $Nu$, the ratio of total heat flux to conductive heat flux, is the quantitative measure of this enhancement. Determining the dependence of $Nu$ on $Ra$, $Pr$ and $\varGamma$ for asymptotically large values of $Ra$ has been a major goal of the studies of convection; see, e.g. Spiegel (1971), Kadanoff (2001), Ahlers, Grossmann & Lohse (2009), Chillà & Schumacher (2012) and references therein. Specifically, if $Nu$ is sought in terms of a power law $Nu = A(Pr, \varGamma )Ra^{\beta }$ then the goal is to determine the value of the exponent $\beta$ for $Ra \gg 1$.

For planar geometries, if one assumes that the dimensional heat flux becomes independent of the depth of the cell as $Ra \rightarrow \infty$, then one obtains $Nu \sim Ra^{1/3}$. This is the so-called classical theory of Priestley (1954), Malkus (1954) and Howard (1966). However, if one assumes that the dimensional heat flux becomes independent of the molecular properties of the fluid when $Ra \rightarrow \infty$, then one obtains $Nu \sim (Pr \, Ra)^{1/2}$. This ‘mixing length’ theory is originally due to Spiegel (1963) and such scaling behaviour – with possible logarithmic corrections (Kraichnan 1962; Chavanne et al. 1997) – is now often referred to as the ultimate regime of thermal convection. The scaling $Nu \sim Ra^{1/2}$ is also an upper limit (uniformly in $Pr$) to the asymptotic heat transport scaling as $Ra \rightarrow \infty$ for no-slip fixed-temperature boundaries whether they are flat (Howard 1963; Doering & Constantin 1996) or corrugated, i.e. textured but sufficiently smooth (Goluskin & Doering 2016). (For flat no-slip boundaries at infinite Prandtl number the best known upper bound corresponds to the classical scaling $Nu \lesssim Ra^{1/3}$ within logarithmic corrections (Constantin & Doering 1999; Doering, Otto & Reznikoff 2006; Otto & Seis 2011). In a wide range of studies at ${O}(1)$ Prandtl number, the exponent $\beta$ is found to vary between $2/7$ (Verzicco & Camussi 2003; Johnston & Doering 2009; Stevens, Verzicco & Lohse 2010; Urban, Musilová & Skrbek 2011; Urban et al. 2012; Doering, Toppaladoddi & Wettlaufer 2019; Iyer et al. 2020) and $1/3$ (Niemela et al. 2000; Niemela & Sreenivasan 2003; Verzicco & Camussi 2003; Niemela & Sreenivasan 2006; Stevens et al. 2010; Urban et al. 2011, 2012; Doering et al. 2019; Iyer et al. 2020). Several experiments have reported $\beta > 1/3$ (Chavanne et al. 1997; He et al. 2012); however, because of the diversity of scalings reported for overlapping ranges of $Ra$, those findings await independent confirmation (Urban et al. 2012; He et al. 2013; Urban et al. 2013; Skrbek & Urban 2015; He, Bodenschatz & Ahlers 2016).

The key difference between the classical ($\beta = 1/3$) and the ultimate ($\beta = 1/2$) theories principally lies in the role played by the thermal boundary layers. In the former regime, thermal boundary layers presumably limit the rate of transport and hence control it (Howard 1966). In the latter regime, the transport of heat is predominantly due to convective motions (Kraichnan 1962; Spiegel 1963). Indeed, these regimes have been observed in recent experiments on radiatively driven convection (Lepot, Aumaître & Gallet 2018; Bouillaut et al. 2019). Hence, it is necessary to investigate the role of thermal boundary layers in turbulent convection to determine the asymptotic high Rayleigh number heat transport.

Motivated by the studies that used surface roughness to probe the boundary layers in turbulent shear flows, Shen, Tong & Xia (1996) studied turbulent thermal convection experimentally in a cell whose top and bottom surfaces were covered with pyramidal roughness elements of aspect ratio $2$, where the aspect ratio is the element width to height. They observed that roughness led to the emission of a larger number of plumes compared to that in convection over smooth surfaces, and that when $Ra$ was above a certain threshold, $Nu$ increased by $20\,\%$ compared to its value for smooth surfaces. However, the value of $\beta \approx 2/7$ was found to be the same as that for planar surfaces for the range of $Ra$ considered. In later experiments, Du & Tong (1998, 2000) concluded similarly. Several subsequent studies, however, report that roughness does lead to an increase in $\beta$ from its planar value (Roche et al. 2001; Qiu, Xia & Tong 2005; Stringano, Pascazio & Verzicco 2006; Tisserand et al. 2011; Salort et al. 2014; Wei et al. 2014; Toppaladoddi, Succi & Wettlaufer 2015a; Wagner & Shishkina 2015; Toppaladoddi, Succi & Wettlaufer 2017; Zhu et al. 2017).

The first study to use roughness to manipulate the interaction between the boundary layers and the outer region to attain the ultimate regime was that of Roche et al. (2001) who studied convection experimentally in a cylindrical cell covered by V-shaped grooves on all sides. They observed that when the thickness of the thermal boundary layers becomes smaller than the amplitude of roughness, $\beta$ attains a value of $0.51$ for $Ra = [2 \times 10^{12}, 5 \times 10^{13}]$. Later, Toppaladoddi et al. (2015a, 2017) used DNS in two dimensions to systematically manipulate this interaction by varying the wavelength of sinusoidal upper and/or lower surfaces at a fixed amplitude. They discovered the existence of an optimal wavelength at which $\beta$ is maximized, and that for wavelengths much smaller and much greater than the optimal wavelength $\beta$ attains its planar value. They also found that $\beta = 0.483$ for the optimal wavelength when both top and bottom surfaces are corrugated (Toppaladoddi et al. 2017). Their findings that roughness wavelength modulation leads to optimal heat transport and results in $\beta \approx 0.5$ for a certain wavelength were subsequently confirmed by experimental (Xie & Xia 2017) and numerical (Zhu et al. 2017) studies, although it has been suggested that $\beta$ can decrease again at even higher $Ra$ (Zhu et al. 2017). More recently, Zhu et al. (2019) reported $Nu \sim Ra^{1/2}$ for $Ra = [10^8, 10^{11}]$ over corrugated surfaces with just three characteristic length scales.

The central physical issue we are addressing here is as follows. As emphasized above, the regimes of determining the exponent $\beta$ centre around the interaction of the thermal boundary layers and the core flow. As the Rayleigh number increases the thermal boundary layers thin. Indeed, as first noted by Niemela & Sreenivasan (2006), one can understand the results of Roche et al. (2001) as a transition between a regime where the groove depth is less than the thermal boundary-layer thickness to a regime where the groove depth is larger than the boundary thickness. Thus, as emphasized by Toppaladoddi et al. (2017), when a given experiment or simulation has a fixed roughness geometry, the boundary-layer core flow interaction may evolve as the Rayleigh number increases. It is for this reason that surfaces with a spectrum of roughness length scales are of interest.

Although we have considerable understanding of the effects of periodic corrugation on plume production and heat transport, it is still not clear a priori if these results can be used to describe the effects of fractal roughness. Indeed, there have been far fewer studies on turbulent convection over multi-scale surfaces, the earliest being that of Villermaux (1998) who theoretically considered the effects of fractal surfaces with power-law distributed amplitudes. Villermaux (1998) argued that given a regime where $Nu \sim Ra^{2/7}$ for smooth boundaries, the effective exponent increased from $2/7$ to $1/3$ with increasing degree of roughness. Ciliberto & Laroche (1999) studied the effects of power-law distributed fractal surfaces on the heat transport experimentally and found larger $\beta$ values of $0.35$ and $0.45$ depending on the distribution of roughness amplitudes. Those studies motivate our own.

In this work we consider the effects of one fractal boundary on the dynamics and bulk transport properties of turbulent Rayleigh–Bénard convection and address the following questions: (i) What are the effects of fractal surface roughness on the heat transport? (ii) How sensitive is the heat transport to the details of the roughness realization? (iii) Can one infer the characteristic length scale(s) of roughness from a study of its effects on the flow? We do this using well resolved two-dimensional (2-D) numerical simulations using the lattice Boltzmann method. The choice of the domain and roughness properties is motivated by our aim to understand the interactions between rough Arctic sea ice and the underlying ocean.

2. Governing equations

The spatial domain in our study (in dimensional units) is $(x,z) \in [0,L]\times [0,h(x)]$ where $0 < h(x) \le H$ is the vertical height of the layer at horizontal position $x$. We model thermal convection via the Oberbeck–Boussinesq equations (Rayleigh 1916; Chandrasekhar 2013), non-dimensionalizing the system using the length scale $H$, the free-fall velocity scale $u_0 = \sqrt {g \alpha {\rm \Delta} T H}$, where $g$ is acceleration of gravity, $\alpha$ is the coefficient of thermal expansion, ${\rm \Delta} T$ is the temperature difference between the bottom and top boundaries and the free-fall time scale $t_0 = H/u_0$.

The equations and boundary conditions for the dimensionless velocity, temperature and pressure fields, $\boldsymbol {u}(\boldsymbol {x},t) = [u(\boldsymbol {x},t), w(\boldsymbol {x},t)]$, $T(\boldsymbol {x},t)$ and $p(\boldsymbol {x},t)$, are

(2.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

(2.2)\begin{gather}\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} = -\boldsymbol{\nabla} p + T \boldsymbol{k} + \sqrt{\frac{Pr}{Ra}} \nabla^2 \boldsymbol{u}, \end{gather}

(2.3)\begin{gather}\frac{\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \sqrt{\frac{1}{Ra \, Pr}} \nabla^2 T, \end{gather}

(2.4a,b)\begin{gather}\boldsymbol{u} = 0 \quad \text{and} \quad T = 1 \quad \text{at}\ z = 0, \end{gather}

(2.5a,b)\begin{gather}\boldsymbol{u} = 0 \quad \text{and} \quad T = 0 \quad \text{at}\ z = h(x). \end{gather}

The Rayleigh number is $Ra = \alpha g {\rm \Delta} T H^3/\kappa \nu$, where $\nu$ is the kinematic viscosity and $\kappa$ is the thermal diffusivity of the fluid, and the Prandtl number is $Pr = \nu /\kappa$. The fractal boundaries are such that $0.9 \le h(x) \le 1$. All variables are periodic in the horizontal $x$ direction, and the aspect ratio of the domain is $\varGamma = L/H$.

The bulk heat transport is measured by the Nusselt number,

(2.6)\begin{equation} Nu = \frac{\left\langle\overline{w^* T^*}\right\rangle - \kappa \left\langle\dfrac{\partial \bar{T}^*}{\partial z^*}\right\rangle}{\kappa {\rm \Delta} T/H}, \end{equation}

evaluated across horizontal layers in the cell. Here, the superscript $*$ indicates the variable is dimensional, and $\overline {(\cdot )}$ and $\left \langle \cdot \right \rangle$ indicate horizontal and time averages, respectively. We compute $Nu$ at eight different heights in the cell and report the average value over these locations.

We use the lattice Boltzmann method (Benzi, Succi & Vergassola 1992; Chen & Doolen 1998; Succi 2001) to solve the governing equations numerically. The principal reason for this choice of numerical method is the ease with which one can impose the boundary conditions for the velocity and temperature fields on complicated domains (Succi 2001). The code used here has been tested extensively in Toppaladoddi, Succi & Wettlaufer (2015b) for different fluid flow problems (Clever & Busse 1974; Lipps 1976; Rozhdestvensky & Simakin 1984) and was previously used to study turbulent convection over planar and corrugated (i.e. smooth but non-flat) upper and lower boundaries (Toppaladoddi et al. 2015a, 2017).

We performed extensive checks on spatio-temporal convergence with the fractal boundaries used in the study. The spatial resolutions used in our simulations were such that the Kolmogorov length scale was resolved with at least two grid points everywhere in the domain (see table 1) and the boundary layer was resolved with at least eight grid points (see table 2). The simulations were run for sufficiently long times to attain a stationary state and the statistics were collected over the last $200$ time units, except for the $Ra = 10^{10}$ cases where the statistics were collected over the last $100$ time units. Details of these tests are provided in appendix A.

We should note here that the depth of the cell $H$ is one of the many choices for the characteristic length scale for this geometry. However, this choice only affects the pre-factor in power-law scalings for Nusselt and Reynolds numbers with Rayleigh number, and not the exponent. (See appendix B.)

3. Roughness profiles

Following Rothrock & Thorndike (1980), we consider upper boundary functions $h(x)$ to be ‘rough’ when they are continuous but not differentiable. The increments in $h(x)$ are given by the Hölder condition,

(3.1)\begin{equation} \lim_{{\rm \Delta} x \to 0} \frac{\left|h(x+ {\rm \Delta} x) - h(x)\right|}{({\rm \Delta} x)^{\gamma}} = C, \end{equation}

where $C$ is an ${O}(1)$ constant and $0 < \gamma \le 1$ is the Hölder exponent. Functions are Lipschitz continuous with a bounded derivative only when $\gamma = 1$. The power spectral density (PSD) of $h(x)$ for all non-zero wavenumbers $k$ decays as $\sim k^{p}$, where $p = -2 \gamma - 1$ (Rothrock & Thorndike 1980). This characteristic decay of the PSD is a common feature shared by many natural and artificial surfaces, and thus can be used to classify different classes of rough surfaces (Sayles & Thomas 1978; Rothrock & Thorndike 1980).

To generate roughness profiles for the upper surface with the desired spectral properties for our simulations, we use the truncated Steinhaus series (Rothrock & Thorndike 1980):

(3.2)\begin{equation} h(x) = h_0 + A \sum_{k=1}^{\mathcal{K}} \left(-p-1\right)^{1/2} k^{p/2} \cos(k x + \phi_k), \end{equation}

where $\mathcal {K}$ is the maximum wavenumber and $\phi _k$ are independent random variables uniformly distributed in $[0, 2{\rm \pi} ]$. It is clear that the PSD of $h(x)$ in (3.2) scales as $\sim k^p$ up to the cutoff wavenumber $\mathcal {K}$.

Figure 1 shows roughness functions for different values of $p$ generated using (3.2). As $p$ is increased from $-3.0$ to $-1.5$, $h(x)$ becomes rough on smaller scales. It is also intuitively clear from figure 1 that with increasing roughness, as $\mathcal {K} \rightarrow \infty$, $h(x)$ tends to be more space filling than a 1-D curve but less space filling than a 2-D surface. Hence these curves are fractals in the limit $\mathcal {K} \rightarrow \infty$, with fractal or Hausdorff dimension $D_f = 2 - \gamma$ (Rothrock & Thorndike 1980).We use (3.2) to generate the rough upper surfaces $h(x)$ for the simulations. All the rough surfaces used in this study have $\mathcal {K} = 100$ and the values of $h_0$ and $A$ are chosen such that their maximum and minimum amplitudes, measured from the top of the cell, are $0\,\%$ and $10\,\%$ of the depth of the cell, respectively. This implies that the upper portion of the fractal boundary coincides with the top flat surface. (See figure 2.)

Figure 1. Functions used for the upper surface of the convecting domain generated using (3.2) for different values of $p$ and for $\mathcal {K} = 100$. The degree of roughness increases as the value of $p$ increases. The curves are vertically displaced by $2$ units to improve their visibility.

Figure 2. Temperature fields at $t = 100$ for $Ra = 2.15 \times 10^9$ and (a) $p = -3.0$, (b) $p = -2.0$ and (c) $p = -1.5$. (See also the movies in the supplementary material available at https://doi.org/10.1017/jfm.2020.826.)

From figure 1 it is clear that there is a distribution of amplitudes associated with the fractal curves; hence, it is not a priori clear what choice of the characteristic length scale would be appropriate. The length scale chosen in our study is the depth of the cell $H$. We present a detailed discussion of this point in appendix B, where we also show that a different choice of length scale simply leads to a uniform rescaling of $Nu$ and $Ra$ values for any given topography.

4. Results

The simulation results are for $Pr = 1$ and $\varGamma = 2$. The simulations, except for $Ra = 10^{10}$, ran to at least $t = 330$ to allow adequate spin up, and in all cases the statistics were obtained for the last $200$ time units. The simulations for $Ra = 10^{10}$ were run for at least $t = 180$, and the statistics were collected over the last $100$ time units. See appendix A for details.

4.1. Temperature fields

Figure 2(a–c) shows the snapshots of the temperature fields for $Ra = 2.15 \times 10^9$ and $p = -3, -2$ and $-1.5$, respectively. Focusing on the region close to the rough upper surfaces, for $p = -3$ (which is comparatively smooth) plumes are emitted only from a fraction of the surface and the temperature field is qualitatively similar to that for convection over flat walls (e.g. Johnston & Doering 2009). As seen in figures 2(b) and 2(c), however, as $p$ increases so too do the number of roughness elements triggering more plume generation: roughness enhances the coupling between the boundary layer and the core flow. Moreover, as seen in the case of a periodically corrugated upper surface (Toppaladoddi et al. 2015a), the enhanced emission of cold plumes decreases the mean interior temperature relative to the planar surface case.

4.2. Variation of heat flux with roughness properties

The $Nu(Ra)$ data are shown in figure 3(a–c) for $Ra \in [10^7, 10^{10}]$ and (a) $p = -3.0$, (b) $p = -2.0$ and (c) $p = -1.5$, respectively. In these figures, for a given value of $p$, the simulations for the whole $Ra$ range were performed using the same realization of the fractal boundary. The scaling fits (i.e. linear least squares of the logarithms) are for $Ra \in [10^8, 10^{10}]$ and are shown as dashed lines in these figures. When $p$ increases from $-3.0$ to $-1.5$, the scaling exponent increases from $\beta = 0.288$ to $\beta = 0.352$. The power-law fit to the $Nu(Ra)$ data for $p=-3.0$ gives $Nu = 0.125 \times Ra^{0.288 \pm 0.005}$, which is remarkably close to $Nu = 0.138 \times Ra^{0.285}$ obtained for a similar $Ra$ range for flat boundaries in two dimensions (Johnston & Doering 2009). This suggests that, for this range of $Ra$, the fractal boundary corresponding to $p=-3.0$ is hydrodynamically smooth with respect to heat transport, even though it lies at the border between smooth and rough boundaries (Rothrock & Thorndike 1980). For $p = -2.0$, the power-law fit gives $Nu = 0.055 \times Ra^{0.329 \pm 0.006}$, which is surprisingly close to $Nu = (0.0525 \pm 0.006) \times Ra^{0.331 \pm 0.002}$ obtained for $Ra \in [10^{10}, 10^{15}]$ for flat boundaries in a slender cylinder with $\varGamma = 0.1$ (Iyer et al. 2020). And, lastly, for $p = -1.5$ the power law is $Nu = 0.037 \times Ra^{0.352 \pm 0.011}$, which is remarkably close to $Nu = 0.034 \times Ra^{0.359}$ obtained for (scaled) wavelength, $\lambda = 0.154$ for a sinusoidally corrugated boundary in two dimensions over $Ra \in [4 \times 10^6, 2.5 \times 10^9]$ (Toppaladoddi et al. 2015a). In fact, Toppaladoddi et al. (2015a) found that $\lambda = 0.154$ was the optimal wavelength that maximized heat transport for their geometry. Evidently, heat transport increases with increasing degree of roughness, i.e. with $p$. This is in qualitative agreement with the results of Villermaux (1998), who also found that the scaling exponent increases with increasing roughness. However, the increase in $\beta$ is due to a change in the dynamics brought about by the introduction of surface roughness, principally the increase in plume production (Stringano et al. 2006; Toppaladoddi et al. 2015a, 2017), and not due to increased surface area. See appendix C for details. We also note that the power-law fits to the whole range of $Ra$ give: (a) $p = -3.0$: $Nu = 0.197 \times Ra^{0.267 \pm 0.015}$; (b) $p = -2.0$: $Nu = 0.111 \times Ra^{0.30 \pm 0.02}$; and (c) $p = -1.5$: $Nu = 0.069 \times Ra^{0.321 \pm 0.023}$.

Figure 3. Plots of $Nu(Ra)$ versus $Ra \in [10^7, 10^{10}]$ for (a) $p = -3.0$, (b) $p = -2.0$ and (c) $p = -1.5$. Circles denote data from simulations and the dashed lines are the linear least-squares fits of $\log Nu$ to $\log Ra$ over the range $Ra \in [10^8, 10^{10}]$. The power laws are fit for the range $Ra \in [10^8, 10^{10}]$. For (a) $p=-3.0$, $Nu = 0.125 \times Ra^{0.288 \pm 0.005}$; (b) $p = -2.0$, $Nu = 0.055 \times Ra^{0.329 \pm 0.006}$; and (c) $p = -1.5$, $Nu = 0.037 \times Ra^{0.352 \pm 0.011}$. The error bar on each $Nu$ data point represents the standard deviation of the averaged $Nu$ calculated from eight different horizontal sections, and the uncertainties in $\beta$ are the $95\,\%$ confidence intervals.

The goodness of power-law fits can be tested by computing the residuals of the actual data from the fit data and examining them for any systematic curvature. These are shown in figure 4. It is apparent that the residuals do not exhibit any curvature for $Ra \in [10^8, 10^{10}]$; hence, it can be concluded that the power laws describe the data well and more general fits are not required.

Figure 4. Residuals of the power-law fits shown in figure 3 for $p = -3.0, -2.0$ and $-1.5$. Here, $Nu_{fit}$ are the values of the Nusselt number obtained from the power-law fits for the range $Ra \in [10^8, 10^{10}]$.

We can estimate an effective hydrodynamic length scale for the roughness amplitude by examining where the $Nu(Ra)$ curves for different values of $p$ intersect. When boundary variations are present, the flow is not influenced by the roughness until the boundary layers become smaller than the amplitude of roughness (Roche et al. 2001). Once this is achieved, as $Ra$ increases further the direct effects of roughness are associated with the increased number of plumes produced and the concomitant augmentation in $Nu$ (Roche et al. 2001). We use similar ideas to estimate the effective roughness of the fractal boundaries. As seen in figures 1–3, the additional roughness structure introduced as $p$ is increased is associated with the increase in $\beta$. If one takes the $Nu-Ra$ curve for $p = -3$ as the benchmark case, then the intersection of this curve with that for a larger value of $p$ gives the value of the effective amplitude at which the transition to enhanced heat transport occurs. The choice of this reference is because $p=-3$ corresponds to $\gamma = 1$, representing the border between ‘smooth’ and ‘rough’ surfaces (Rothrock & Thorndike 1980). Thus the effects of any additional roughness (see figure 1) can be conveniently studied with respect to the surface for $p=-3$. This transition happens at $Ra \approx 2.15 \times 10^8$, and the value of $Nu$ at this point is $\approx 31$ (see the $Nu$ values for the fourth realizations ($r = 4$) in tables 3–5). Hence, the transition occurs when the effective amplitude of roughness $h_f$ over the surfaces with $p = -2$ and $p = -1.5$ first exceeds the boundary layer thickness $\delta _T$ for the curve with $p = -3$, so that the roughness elements protrude outside of the boundary layer and interact with the interior of the flow. Using the planar-wall estimate of $Nu$, we estimate the cross-over scaling as

(4.1)\begin{equation} \delta_T = h_f \approx \frac{1}{2 Nu} = 0.016. \end{equation}

Thus, the effective amplitude of the roughness for surfaces with $p=-2$ and $p=-1.5$ is approximately $2\,\%$ of the depth of the cell.

4.3. Sensitivity of $Nu$ to details of roughness realization

To investigate the effects of a given roughness realization on heat transport, we computed $Nu(Ra)$ for four different realizations for each value of $p$. To generate each realization for a fixed $p$, we have used different values of $\phi _k$. However, the first realization for all $p$ values have the same set of $\phi _k$. Similarly for second, third, and fourth realizations. The $Nu(Ra)$ curves from these simulations are shown in figure 5.

Figure 5. Plots of $Nu(Ra)$ data for four different realizations for each p: (a) $p = -3.0$, (b) $p = -2.0$ and (c) $p = -1.5$. The first realizations for the different values of p are generated using the same set of $\phi_k$ (see (3.2)). The error bar on each $Nu$ data point represents the standard deviation of the averaged $Nu$ calculated from eight different horizontal sections.

It is seen from figures 5(a)–5(c) that for a fixed $p$, the $Nu$ depends primarily on the $Ra$ with very little dependence on the realization. Hence, to a good approximation, the heat transport for the fractal surfaces used here depends only on the latter's statistical properties, i.e. $p$, and in turn on $D_f$. Hence, this suggests that the scaling exponents in § 4.2 depend uniquely on $p$.

Furthermore, to compute higher-order moments, we have run simulations for $Ra = 10^8$ and $t = 875$ for all the roughness realizations. The maximum variations in the means of $Nu(t)$ measured at $z=0$ between ensemble members for $p = -3.0, -2.0$ and $-1.5$ are $3.3\,\%$, $1\,\%$ and $0.2\,\%$, respectively. Similarly, the maximum variations in the standard deviations for $p = -3, -2$ and $-1.5$ are $5.4\,\%$, $16\,\%$ and $9.1\,\%$, respectively. The variations in the higher-order moments (skewness and kurtosis) are relatively larger. This shows that the mean of $Nu(t)$ is less sensitive to the details of the roughness than its higher-order moments.

4.4. Reynolds number and its sensitivity to the details of the roughness realization

In addition to considering the bulk heat transport, we also studied the behaviour of the bulk Reynolds number ($Re$) with $Ra$ and $p$ to further characterize the response of the flow. The Reynolds number is

(4.2)\begin{equation} Re = \frac{U_0 H}{\nu}, \end{equation}

where $U_0$ is a velocity scale, the choice of which is not unique. Previous studies over smooth (Niemela et al. 2001; Qiu & Tong 2001; Sun & Xia 2005; Niemela & Sreenivasan 2006) and regular rough surfaces (Wei et al. 2014) have either constructed $U_0$ based on the depth of the cell and the dominant frequency of oscillations of the large-scale circulation, or used a root-mean-squared (r.m.s.) velocity deduced from single-point measurements. We take $U_0 = U_{rms}$, where $U_{rms}$ is the bulk averaged r.m.s. velocity computed over all the nodes in the domain.

Figure 6 shows $Re(Ra)$ data along with power-law fits $Re \sim Ra^{\xi }$ for the three different $p$. Unlike $\beta$, the exponent $\xi$ characterizes scaling behaviour of $Re$ over three full decades of $Ra$. Moreover, $Re(Ra)$ is substantially less sensitive to details of the roughness: $\xi \approx 0.57$ for all three values of $p$ and the prefactor variation among the three values of $p$ is less than 8 %. This suggests that the strength of the velocity variations in the cell is set by the large scale properties of the boundary profile that are present for the smooth surface with $p=-3$, and that smaller scale roughness does not appreciably affect $\xi$. Recent observation of $Re \sim Ra^{0.617}$ scaling for turbulent convection over flat boundaries in two dimensions (Wan et al. 2020) is consistent with this suggestion.

Figure 6. Plot of $Re(Ra)$ versus $Ra \in [10^7, 10^{10}]$ for $p = -3.0$, $p = -2.0$ and $p = -1.5$. Symbols denote data from simulations and the dashed lines are the linear least-squares fits of $\log Re$ to $\log Ra$ for the whole $Ra$ range. The power-law fits for the different values of p are: $p=-3.0$, $Re = 0.094 \times Ra^{0.571 \pm 0.018}$; $p = -2.0$, $Re = 0.087 \times Ra^{0.576 \pm 0.022}$; and $p = -1.5$, $Nu = 0.091 \times Ra^{0.571 \pm 0.017}$. The uncertainties in the values of $\xi$ are the $95\,\%$ confidence intervals.

Note that $\xi = 0.5$ corresponds to the (dimensional) r.m.s. fluid speed being proportional to the free-fall velocity across the cell, $u_0 = \sqrt {g \alpha {\rm \Delta} T H}$. Because the boundary temperatures are fixed, $g \alpha {\rm \Delta} T$ is the maximal buoyancy acceleration of any fluid element, so suitably conspiratorial flow configurations would be required to sustain $\xi > 0.5$ as $Ra \rightarrow \infty$. Such is the case in coherent steady (albeit unstable) convection between stress-free boundaries where $\xi \rightarrow 2/3$ as $Ra \rightarrow \infty$ (Chini & Cox 2009).

In order to understand the effects of details of roughness realization on the variation of $Re(Ra)$, we performed an analysis similar to that reported for $Nu(Ra)$ in § 4.3. Figure 7 shows $Re(Ra)$ data for all $p$ and the different realizations. The realizations used here are the same as those used in § 4.3. Figure 7 clearly suggests that $Re(Ra)$ is independent of the details of the roughness realizations and the value of $p$ itself. This is further supported by the fact that the power-law fits to $Re(Ra)$ data for all four different realizations for each $p$ and $Ra \in [10^7,2.15\times 10^9]$ give: (a) $p=-3.0$: $Re = 0.073 \times Ra^{0.584}$; (b) $p = -2.0$: $Re = 0.069 \times Ra^{0.588}$; and (c) $p = -1.5$: $Re = 0.068 \times Ra^{0.589}$. Hence, unlike $\beta$, $\xi$ is independent of the roughness geometries used in this study.

Figure 7. Plot of $Re(Ra)$ data for all $p$ and the different realizations denoted r1–r4. The power-law fits shown here are those reported in figure 6 for $Ra \in [10^7, 10^{10}]$.

5. Conclusions

We have systematically studied turbulent thermal convection in two-dimensional domains with a fractal upper boundary for $Ra \in [10^7, 10^{10}]$ using the lattice Boltzmann method. The fractal nature of the boundaries is characterized by their spectral exponent $p=2 D_f -5$ representing the degree of roughness, where $D_f$ is the Hausdorff dimension of the boundary function. Simulations with roughness exponents $p = -3.0, -2.0$ and $-1.5$ revealed the following.

1. (i) With increasing roughness, the fractal boundaries provide an increasing number of sites for the generation of plumes. Hence, at fixed $Ra$ the plume production increases with increasing $p$.

2. (ii) The $Nu \sim Ra^{\beta }$ power-law fit exponent $\beta$ for the range $Ra \in [10^8, 10^{10}]$ increased from $0.288$ to $0.352$ as $p$ increased from $-3.0$ to $-1.5$. Heat transport increased with roughness for larger $Ra$, in qualitative agreement with Villermaux (1998). This increase in $\beta$ is due to a change in the dynamics that results from the enhanced interactions between the rough boundary and the inner flow, through the increase in plume production (Stringano et al. 2006; Toppaladoddi et al. 2015a, 2017). The increased surface area over the boundaries also increases heat transport, but only by increasing the pre-factor in the power law and not $\beta$. (See appendix C.)

3. (iii) The fractal surfaces used in the experiments of Ciliberto & Laroche (1999) were built with glass spheres such that the amplitudes were power-law distributed, i.e. $P(h) \sim h^{\varLambda }$. They reported $\beta = 0.35$ for $\varLambda = -2$ and $\beta = 0.45$ for $\varLambda = -1$. Hence, $\beta$ increased with increasing $\varLambda$, which represents the degree of roughness. Hence, our findings are qualitatively consistent with the results of Ciliberto & Laroche (1999).

4. (iv) The following observations can be made based on the analysis of $Nu-Ra$ data: (a) the fractal boundary for $p=-3.0$ is hydrodynamically smooth for heat transfer, as both the exponent and prefactor in the $Nu-Ra$ power law are approximately those observed in convection over flat boundaries (Johnston & Doering 2009); (b) both the prefactor and exponent of the $Nu-Ra$ power law for $p = -2.0$ correspond surprisingly well, albeit perhaps only fortuitously, to those reported for convection over smooth surfaces in a small aspect ratio 3-D cylindrical geometry (Iyer et al. 2020); and (c) the $Nu-Ra$ power law for $p=-1.5$ is remarkably close to the one obtained for the optimal wavelength of a corrugated sinusoidal boundary that maximizes heat transport (Toppaladoddi et al. 2015a).

5. (v) Using the roughness curve for $p = -3$ as the reference profile, we estimated the effective amplitudes of the roughness curves for $p=-2$ and $-1.5$ that lead to increased heat transport. This is about $2\,\%$ of the depth of the cell for $p = -2$ and $p = -1.5$.

6. (vi) The averaged $Nu$ values for a fixed $p$ depend primarily on $Ra$ and have a very weak dependence on the roughness realization. However, the higher-order moments are more sensitive to the details of roughness realizations.

7. (vii) The Reynolds numbers based on the r.m.s. velocity computed over all fluid nodes scaled as $Re \sim Ra^{\xi }$, with $\xi \approx 0.57$, for all three values of $p$ studied here. Perhaps surprisingly, the bulk intensity of the flow was substantially less sensitive to small-scale details in the roughness profiles than the heat transfer.

8. (viii) Like the averaged $Nu$ values, the averaged $Re$ values for a fixed $p$ depend primarily on $Ra$, and have a weak dependence on the roughness realization.

9. (ix) To a good approximation, the exponent $\beta$ is solely a function of $p$ (and in turn $D_f$), whereas $\xi$ is independent of $p$.

These simulations demonstrate the feasibility of studying turbulent flows over fractal walls using numerical simulations. Importantly, they provide a framework to study heat transport in high $Ra$ convection that can reveal the influence of interactions between the boundary layers and core flow. Namely, we know that such interactions are important for the $Nu(Ra)$ behaviour and that as $Ra$ increases, boundary layers thin and so too will the size of roughness elements that trigger plume production. For a given fractal surface, only a fraction of the roughness elements are driving boundary-layer instability and that fraction changes with $Ra$. Therefore fractal surfaces that enhance plume production and heat transport must also optimize the fraction of the ‘active’ surface roughness elements. However, although a fractal surface reveals finer details with increasing resolution, all numerical simulations are ultimately limited by finite resolution. Hence there will always be details of the surface that the simulated flow would not be able to sense. This leads naturally to the question of how one can represent the effects these unresolved details of roughness on the turbulent flows, a perennial conundrum in all manner of flows adjacent to surfaces.