Thermal Convection over Fractal Surfaces

We use well resolved numerical simulations to study Rayleigh-B\'enard convection in cells with a fractal boundary in two dimensions for $Pr = 1$ and $Ra \in \left[10^7, 2.15 \times 10^9\right]$. The fractal boundaries are functions characterized by power spectral densities $S(k)$ that decay with wavenumber 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$ in power law fits $Nu \sim Ra^{\beta}$, where $Nu$ and $Ra$ are the Nusselt and the Rayleigh numbers, we observe that heat transport increases with roughness. For $p$ $= -3.0$, $-2.0$ and $-1.5$ we find, respectively, $\beta = 0.256, 0.281$ and $0.306$. For a given value of $p$ we observe that the mean heat flux is insensitive to the specific realization of the roughness


INTRODUCTION
Thermal convection refers to fluid flows that are driven by buoyancy forces due to density variations, which in turn are effected by gradients in temperature [1]. 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 [2,3] 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 [4]. Applying the Boussinesq approximation to the Navier-Stokes equations, the dynamics of RBC are governed by three non-dimensional parameters: the Rayleigh number Ra, the ratio of buoyancy to viscous forces; the Prandtl number P r, the ratio of the fluid's kinematic viscosity to its thermal diffusivity, and the aspect ratio Γ 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 N u, the ratio of total heat flux to conductive heat flux, is the quantitative measure of this enhancement. Determining the dependence of N u on Ra, P r, and Γ for asymptotically large values of Ra has been a major goal of the studies of convection; see, e.g., [2,5,6] and refs therein. Specifically, if N u is sought in terms of a power-law N u = A(P r, Γ) Ra β then the goal is to determine the value of the exponent β for Ra ≫ 1.
The key difference between the classical (β = 1/3) and the ultimate (β = 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 [9]. In the latter regime, the transport of heat is predominantly due to convective motions [10,11]. Indeed, these regimes have been observed in recent experiments on radiatively driven convection [32,33]. 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 et al. [34] studied turbulent thermal convection experimentally in a cell whose top and bottom surfaces were covered with pyramidal roughness elements of aspect ratio (defined as the ratio of their width to height) 2. 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, N u increased by 20% compared to its value for smooth surfaces. However, the value of β ≈ 2/7 was found to be the same as that for planar surfaces for the range of Ra considered. The later experiments of Du and Tong [35,36] concluded similarly. Several subsequent studies, however, have shown that surface roughness does lead to an increase in β from its planar value [37][38][39][40][41][42][43][44][45].
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. [37] 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, β attains a value of 0.51 for Ra = 2 × 10 12 , 5 × 10 13 . Later, Toppaladoddi et al. [44,45] used DNS in two-dimensions (2D) 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 β is maximized, and that for wavelengths much smaller and much greater than the optimal wavelength β attains its planar value. They also found that β = 0.483 for the optimal wavelength when both top and bottom surfaces are corrugated [45]. Their findings were subsequently confirmed by experimental [46] and numerical [47] studies, although there is suggestion that the exponent can decrease again at even higher Ra [47]. More recently, Zhu et al. [48] reported observing N u ∼ Ra 1/2 for convection over rough surfaces consisting of three characteristic amplitudes for Ra = 10 8 , 10 11 .
The central physical issue we are addressing here is as follows. As emphasized above, the regimes of determining the exponent β center 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 and Sreenivasan [26], one can understand the results of Roche et al. [37] 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. [45], 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 could 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 [49] who theoretically considered the effects of fractal surfaces with power-law distributed amplitudes. That study suggested that the increase in the effective area of the surfaces leads to enhancement in the values of N u compared to its planar values, and that the effective exponent for N u(Ra) increased from 2/7 to 1/3 with increasing degree of roughness. Ciliberto and Laroche [50] studied the effects of power-law distributed fractal surfaces on the heat transport experimentally and found larger β 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: (1) What are the effects of fractal surface roughness on the heat transport? (2) Is the heat transport sensitive to the details of the roughness distribution? (3) 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 2D numerical simulations using the Lattice Boltzmann Method.

GOVERNING EQUATIONS AND NUMERICAL METHOD
The spatial domain in our study is ( ≤ H is the vertical height of the layer at horizontal position x. We model thermal convection via the Oberbeck-Boussinesq equations [1,4], non-dimensionalizing the system using the length scale H, the free-fall velocity scale u 0 = √ g α ∆T H where g is acceleration of gravity, α is the coefficient of thermal expansion, ∆T is the temperature difference between the bottom and top boundaries, and the free-fall time scale The equations and boundary conditions for the dimensionless velocity, temperature, and pressure fields u(x, t), T (x, t), and p(x, t) are The Rayleigh number Ra = α g ∆T H 3 /κ ν, where ν is the kinematic viscosity and κ is the thermal diffusivity of the fluid, the Prandtl number P r = ν/κ. All variables are periodic in the horizontal x direction, and the aspect ratio of the domain is Γ = L/H.
The bulk heat transport is measured by the Nusselt number, evaluated at horizontal layers in the cell. Here the superscript * indicates the variable is dimensional, and (·) and · indicate, respectively, horizontal and time averages. We compute N u at eight different heights in the cell and report the average value over these locations. The standard deviation from the mean over these locations was less than 2% in the simulations. We use the Lattice Boltzmann Method [51][52][53] 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 [53]. The code used here has been tested extensively in [54] for different fluid flow problems [55][56][57] and was previously used to study turbulent convection over planar and corrugated (i.e., smooth but non-flat) upper and lower boundaries [44,45].
We performed further checks on spatio-temporal convergence with fractal rough boundaries used in this study. Convergence tests were performed for the case Ra = 2.15 × 10 9 and p = −1.5 (See Eq. 8) using three resolutions: (N x , N z ) = (2800, 1600), (2400, 1200), and (2000, 1000). The lowest resolution simulation was run until t = 583 and the highest resolution simulation until t = 325; the statistics were collected for the last 200 time units. The maximum variation in N u obtained from these runs was 1%, quantifying the degree of convergence of N u.

ROUGHNESS PROFILES
Following Rothrock and Thorndike [58], 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 where C is an O(1) constant and 0 < γ ≤ 1 is the Hölder exponent. Functions are Lipschitz continuous with bounded derivative only when γ = 1. The power spectral density (PSD) of h(x) for all non-zero wavenumbers k decays as ∼ k p , where p = −2γ −1 [58]. This characteristic decay of the PSD is a common feature shared by many natural and artificial surfaces [58,59], and thus can be used to classify different classes of rough surfaces [58].
To generate roughness profiles for the upper surface with the desired spectral properties for our simulations, we use the so-called truncated Steinhaus series [58]; where K is the maximum wavenumber and φ k are independent random variables uniformly distributed in [0, 2π]. The factor (−p − 1) 1/2 is included to make the variance independent of p [58]. It is clear that the PSD of h(x) in (8) scales as ∼ k p up to the cutoff wavenumber K. Figure 1 shows roughness functions for different values of p generated using equation (8). 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 K → ∞, h(x) tends to be more space filling than a 1D curve but less space filling than a 2D surface. Hence these curves are fractals in the limit K → ∞, with fractal or Hausdorff dimension D f = 2 − γ [58]. We use equation (8) to generate the rough upper surfaces h(x) for the simulations. All the rough surfaces used in this study have K = 100 and are such that their maximum amplitude, measured from the top of the cell, is 10% of the depth of the cell.

RESULTS
The simulation results are for P r = 1 and Γ = 2. Simulations described below 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. cusing 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., 22]. As seen in figures 2(b) and (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 [44], the enhanced emission of cold plumes decreases the mean interior temperature relative to the planar surface case.

Variation of heat flux with roughness properties
The N u(Ra) data are shown in figures 3(a)-(c), along with their scaling fits (i.e., linear least squares of the logarithms) to N u = A Ra β , for Ra = 10 7 , 2.15 × 10 9 and p = −3, −2 and −1.5. When p increases from −3 to −1.5 the value of β increases from 0.256 to 0.306, thereby demonstrating the enhanced heat flux at increasingly large Ra. However, in this range of Ra, relative to the case of periodic roughness with varying wavelength and fixed amplitude [44], the increase in β for these fractal boundaries is modest. Whereas plume production is imposed at a single wavelength in the periodic case, which leads to an optimization of heat transport, the distribution of asperates in the fractal case is associated with enhanced plume production on only a fraction of the surface.
Another apparent feature in figures 3(b) and (c) is the change in the local slope of the N u(Ra) data. Although the range of Ra studied here does not permit fitting local scaling exponents to the data with a break in the slope, we can estimate an effective hydrodynamic length scale for the roughness amplitude. When boundary variations are present, the flow is not influenced by the roughness until the boundary layers become smaller than the amplitude of roughness. 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 N u [37]. We use similar ideas to estimate the effective roughness of the fractal walls. As seen in figures 1 and 2, the additional roughness structure introduced as p is increased is associated with the increase in β. If one takes the N u − 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 γ = 1, representing the border between "smooth" and "rough" surfaces [58]. Thus the effects of any additional roughness (see figure 1) can be conveniently studied with respect to the surface for p = −3. From figures 3(b) and (c) this transition happens at Ra ≈ 2.15 × 10 8 , and the value of N u at this point is ≈ 31. 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 δ 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 N u, we estimate the crossover scaling Thus the effective amplitude of the roughness for surfaces with p = −2 and p = −1.5 is about 2% of the depth of the cell.
We also note that the increase of the scaling exponent β from 0.256 to 0.306 as the degree of roughness is increased is in qualitative agreement with the theoretical study of Villermaux [49], who suggested that β should increase from 2/7 to 1/3 as the degree of roughness increases.

Sensitivity of heat flux to details of the roughness distribution
To investigate the effects of the details of the roughness on the statistics of the heat flux we examine the time series of the instantaneous Nusselt number, N u(t), at the flat bottom wall for four different realizations of top boundary roughness for each p. Figures 4(a)-(c) show the probability density functions for the N u(t) data for different p. For each case presented we consider Ra = 10 8 and the length of the N u(t) time series is 200.
Two features are apparent: (i) the distributions are monomodal and reasonably symmetric about their peaks, and (ii) there are slight, possibly significant, variations between different realizations for each p. The maximum variations in the means of N u(t) between ensemble members for p = −3, −2, and −1.5 are, respectively, 1%, 3% and 4.2%. Similarly, the maximum variations in the variances for p = −3, −2, and −1.5 are, respectively, 17%, 6%, and 6.5%. The variations in the higher-order moments (skewness and kurtosis) are relatively larger, which might partly be due to insufficient convergence. This suggests that the mean of N u(t) is less sensitive to the details of the roughness than its higher-order moments, i.e., at a given value of p the time averaged bulk heat transport is less sensitive to details of the roughness realization than fluctuations in the bulk heat transport.

Reynolds number
In addition to considering the bulk heat transport, we also studied the behavior of the bulk Reynolds number (Re) with Ra and p to further characterize the response of the flow. The Reynolds number is where U 0 is a velocity scale, but the choice of U 0 is not unique. Previous studies over smooth [26,[60][61][62] and regular rough surfaces [41] 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 rootmean-squared (RMS) velocity deduced from single-point measurements. We take U 0 = U rms , where U rms is the bulk averaged RMS velocity computed over all the nodes in the domain. Figures 5(a) -5(c) show Re(Ra) data along with power-law fits Re ∼ Ra ξ for the three different p. Unlike β, the exponent ξ characterizes convincing scaling behavior of Re as a function of Ra over more than two decades in the Rayleigh number. Moreover, Re(Ra) is substantially less sensitive to details of the roughness: ξ ≈ 0.59 for all three values of p and the prefactor variation among the three values of p is of the order of the variation of the individual data from their scaling fits. This suggests that Note that ξ = 0.5 would imply that the characteristic (dimensional) fluid speed corresponds to the free-fall velocity across the cell, u 0 = √ g α ∆T H. The robustly larger value of ξ observed in our simulations likely indicates a transient Re(Ra) scaling behavior: we expect ξ to be bounded from above by 0.5 as Ra → ∞. Indeed, because the boundary temperatures are fixed, gα∆T is the maximal buoyancy acceleration of any fluid element so an alternative mechanism, i.e., suitably conspiratorial flow configurations, would be required to sustain substantially higher bulk averaged speeds.

CONCLUSIONS
We systematically studied turbulent thermal convection in domains with a fractal upper boundary for Ra ∈ 10 7 , 2.15 × 10 9 in two-dimensions using the Lattice Boltzmann Method. The fractal nature of the boundaries are characterized by their spectral exponent p = 2D f − 5 representing the degree of roughness where D f is the Hausdorff dimension of the boundary function. Simulations with roughness exponents p = −3, −2 and −1.5 revealed the following: 1. 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. The N u ∼ Ra β power-law fit exponent β increased from 0.256 to 0.306 as p varied from −3 to −1.5.
Heat transport increased with roughness for larger Ra in qualitative agreement with [49].
3. Based on the increase in β with p, we conclude that the enhancement in the interaction between flow in the boundary and bulk regions is not as strong as that for convection over periodic corrugation with an optimally tuned wavelength [44]. This might be attributed to the absence of a dominant wavelength in the roughness geometries.
4. 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 apparently about 2% of the depth of the cell for p = −2 and p = −1.5.

5.
Instantaneous heat flux measurements revealed mono-modal and approximately symmetric distributions about their means. They also showed mean values that are less sensitive than the higher-order fluctuations to the details of individual statistical realizations of the roughness.
6. The Reynolds numbers based on the RMS velocity computed over all fluid nodes scaled as Re ∼ Ra ξ , with ξ ≈ 0.59, 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.
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 at the core of the N u(Ra) behavior 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 have finite resolution so there will always be details of the surface that the flow would not be able to sense. This leads to question of how one can represent the effects these unresolved details of roughness on the tur- Close examination of figures 3(a) -3(c) reveals the convexity of the N u(Ra) data. To test the goodness of the power-law fits and quantify the the convexity we fit higher-order polynomials to the log(N u)-log(Ra) data and examined the residuals and the implied local slopes β = d log(N u)/d log(Ra).
Figures 6(a) -6(c) show the residuals for linear, quadratic, and cubic fits to the log(N u)-log(Ra) data for the different p's. It is apparent from the figures that: (a) the residuals of the linear fit have systematic curvature indicating convexity in the data, and (b) the residuals of the quadratic and cubic fits are, in comparison, structureless, indicating that they fit the data better than pure scaling and are, essentially, statistically equivalent.
Local exponents obtained by differentiating the log(N u)-log(Ra) fits to the data are shown in figures 7(a) -7(c). While there is clear systematic increase in the local exponents for each of the three values of p, the trends of increasing exponent with increasing p indicated by the single scaling-fit exponent β are still apparent at higher values of Ra.