Turbulent spherical Rayleigh–Bénard convection: radius ratio dependence

Abstract

Direct numerical simulation (DNS) is performed to explore turbulent Rayleigh–Bénard convection in spherical shells. Our simulations cover six distinct values of radius ratio, $\eta = r_i/r_o = 0.2$, 0.3, 0.4, 0.5, 0.6 and 0.8, under the assumption of a centrally condensed mass with gravity profile $g \sim 1/r^{2}$; where $r_i$, $r_o$ and $r$ denote the inner shell radius, the outer shell radius and the local radial coordinate, respectively. The Prandtl number is kept constant at unity while the Rayleigh number ($Ra$) is varied from $3 \times 10^{3}$ to $5 \times 10^8$. Our primary aim is to analyze how the radius ratio influences the global transport properties and flow physics. To gain insights into the scaling behaviour of the Nusselt number ($Nu$) and the Reynolds number ($Re$) with respect to $Ra$ and $\eta$, we apply the Grossmann–Lohse (GL) theory (Grossmann & Lohse, J. Fluid Mech., vol. 407, 2000, pp. 27–56) to the system. It is observed that the scaling exponents for $Nu$ and $Re$ in relation to $Ra$ are more significant for smaller $\eta$ values, suggesting that the simulations with smaller $\eta$ reach the classical $Nu\sim Ra^{1/3}$ regime at a relatively lower $Ra$. This observation could also imply the systems with smaller $\eta$ might transition to the ultimate regime earlier at a smaller $Ra$. Based on our extensive DNS data, we establish that the thickness of the inner thermal boundary, $\lambda _{\vartheta }^{i}$, follows a scaling relationship of $\lambda _{\vartheta }^{i} \sim \eta ^{1/2}$. This relationship, in turn, leads to a scaling law for $Nu$ in the form of $Nu \sim f(\eta ) Ra^{\gamma }$, where the function $f(\eta )$ is defined as $f(\eta ) = {\eta ^{1/2}}/{(1+\eta ^{4/3})}$, and the exponent $\gamma$ depends on both $Ra$ and $\eta$. Additionally, we characterize and explain the asymmetry in the velocity field by introducing the separate Reynolds numbers for the inner and outer shells. The asymmetry of the kinetic and thermal energy dissipation rates in the inner and outer boundary layers (BLs) is also quantified.

1. Introduction

Thermal convection in spherical shells is ubiquitous in geophysical and astrophysical processes such as tropical convection at the Earth's atmospheric boundary (Hartmann, Moy & Fu 2001), mantle convection within the Earth (Davies & Richards 1992), thermal convection in Earth's outer core (Wicht & Sanchez 2019), deep convection in the molecular envelopes of gas giants (Aurnou et al. 2008) and thermal convection in the convective zone of the Sun (Hanasoge, Gizon & Sreenivasan 2016). In spherical shells, the curvature of the inner boundary differs from that of the outer boundary. Additionally, gravitational effects vary radially and are not constant (Phillips & Lambeck 1980). Consequently, the flow properties near the inner boundary are distinct from those at the outer boundary, and both are different to the bulk. Rayleigh–Bénard convection (RBC) within spherical shells (Gastine, Wicht & Aurnou 2015), where heat is supplied from the inner sphere and dissipated at the outer sphere, serves as a paradigmatic system for investigating thermal convection in spherical geometries. Similar to its planar counterpart, a fundamental aspect of studying spherical RBC is to ascertain how the system responds, such as determining dimensionless heat transfer (Nusselt number, $Nu$) and dimensionless momentum transfer (Reynolds number, $Re$), as functions of the input parameters (Rayleigh number, $Ra$, and Prandtl number, $Pr$). For a more comprehensive introduction to RBC in general, and specifically planar RBC, which has been extensively studied over recent decades, the interested reader is referred to the papers of Ahlers, Grossmann & Lohse (2009), Chillà & Schumacher (2012), Plumley & Julien (2019) and Xia et al. (2023).

To investigate RBC in spherical shell geometry, several laboratory experiments have been conducted. To eliminate vertical gravity and create a radially inward body force, fluids that can be influenced by electric fields are employed in micro-gravity environments. In these experiments, the radial body force is induced by the electric field to mimic buoyancy. A series of experiments called ‘GeoFlow’ have been performed at the International Space Station (ISS) (Egbers et al. 2003; Futterer et al. 2010, 2013; Zaussinger et al. 2018, 2020) to study the flow instability, pattern formation and transition to chaos of thermal convection in concentric, non-rotating and rotating spherical shells. Since these experiments focused on understanding the thermal convection in Earth's mantle, the Prandtl number considered was in the range of $40 < Pr < 200$. As a result of the limitation of the electric field strength and the fluid dielectric properties, the Rayleigh number was limited to a relatively small value, $Ra < 10^{7}$. In recent decades, several numerical studies (Zebib, Schubert & Straus 1980; Zebib et al. 1983; Bercovici et al. 1989; Bercovici, Schubert & Glatzmaier 1992; Tilgner 1996; Tilgner & Busse 1997; Choblet & Parmentier 2009; Choblet 2012; Gastine et al. 2015) have focused on non-rotating spherical shell thermal convection; some of which consider an infinite $Pr$ (Bercovici et al. 1992), simulating the Earth's mantle. Only a few of them investigated the scaling properties of the response parameters of the system as a function of the driving forces. For example, Tilgner (1996) studied $Nu(Ra, Pr)$ and $Re(Ra, Pr)$ scalings in the parameter space of $0.06 \leq Pr \leq 10$ and $4 \times 10^{3} \leq Ra \leq 8 \times 10^{5}$ at a fixed radius ratio $\eta =0.4$ and a gravity profile ($g(r) \sim r$). The scalings obtained were $Nu \sim Ra^{0.24}$ and $Re \sim Ra^{0.5}$. The scarcity of scaling studies on spherical RBC in the literature is a key motivation for this research. One of the primary goals of the present work is to explore a wide range of parameter space to examine the scaling behaviours of $Re$ and $Nu$ as functions of $Ra$ and $\eta$.

To quantify the thermal boundary layer (BL) asymmetry and scaling properties in spherical RBC, Gastine et al. (2015) numerically examined spherical thermal convection with different gravity profiles ($g(r) \in \{r/r_{o},1,(r_{o}/r)^{2}, (r_{o}/r)^{5} \}$) and a broad range of radius ratios ($0.2 \leq \eta \leq 0.95$), at $Pr=1$. In their approach, they assume the average plume density on the inner and the outer boundaries to be the same. Based on this assumption, Gastine et al. (2015) calculated the ratio of the inner and the outer BL thicknesses. The temperature drop within each of these BLs was also quantified. It is important to note that the scaling behaviours of $Nu(Ra)$ and $Re(Ra)$ were studied exclusively for the simulations with $Ra$ up to $10^9$ and at a specific radius ratio of $\eta = 0.6$. For all other radius ratios, the primary focus was on studying BL asymmetry, and as a consequence, only very low $Ra$ values were considered. Scaling relations in these cases were not established, leaving ample room for further research to explore the impact of radius ratio on the scaling of $Nu(Ra)$ and $Re(Ra)$.

Our primary objective in this study is to unravel the influence of the radius ratio on the scaling behaviour of dimensionless heat transport and flow velocity in spherical RBC. This investigation holds significance specifically in the context of the convective zones of different planets. These planetary convective zones exhibit varying sizes, leading to different radius ratios. For example, recent geodetic analyses of Mercury's interior have indicated the presence of a solid inner core and a liquid outer core, with the radius ratio falling within the range of 0.3–0.7 (Genova et al. 2019). Similarly, Earth's outer core and inner core possess a radius ratio of approximately 0.35 (Ahrens 1995). The situation becomes a bit more complex when considering gas giants like Jupiter and Saturn, where defining the convective zone depends on whether the core-metallic hydrogen interface or the core-molecular hydrogen interface is chosen as the inner boundary. Consequently, the radius ratio for Jupiter's convective zone varies between 0.55 and 0.75 (Guillot et al. 2004), and for Saturn, adopting the core-metallic hydrogen interface as the inner boundary results in a radius ratio spanning from 0.45 to 0.55 (Christensen & Wicht 2008; Vazan et al. 2016). Additionally, the study of gas giants is complicated by the uncertainty in their boundary conditions. Unlike the present study, gas giants may not have no-slip boundaries. It is important to acknowledge that the radius ratio values for other planets in our solar system remain uncertain. For instance, the inner core radius of Venus remains unclear (O'Neill 2021; Stähler et al. 2021). Mars, however, is highly suspected not to possess a solid inner core (Stähler et al. 2021). For the ice giants Uranus and Neptune, their interior structures remain a subject of controversy and uncertainty (Helled, Nettelmann & Guillot 2020). The radius ratios used in simulations and experiments pertaining to various planets in the Solar system are listed in table 1. In many geophysical flows, the effect of rotation must be taken into consideration. Recent studies (e.g. Gastine, Wicht & Aubert 2016; Long et al. 2020; Song et al. 2024a; Song, Shishkina & Zhu 2024b) indicate that the behaviour of heat transport in rotating convection approaches that of the non-rotating convection when the Rayleigh number is sufficiently high for a given Ekman number, such that the buoyancy effects dominate over the rotational effects.

Table 1. Radius ratios for the selected planets in the solar system. For Mercury and Earth, the radius ratio $\eta$ is defined as the ratio of the inner solid core radius and the outer liquid iron core radius (Ahrens 1995; Genova et al. 2019). The $\eta$ for Jupiter (Guillot et al. 2004) and Saturn (Christensen & Wicht 2008; Vazan et al. 2016) are chosen as the ratio between the radii of the inner metallic-core boundary and the outer upper atmosphere boundary.

In this study, we conduct 97 three-dimensional (3-D) DNS of spherical shell RBC with $0.2 \leq \eta \leq 0.8$ and $Ra$ up to $5 \times 10^{8}$, at $Pr=1$. We begin by looking at the flow structures in our simulations and report the existence of large-scale structures in spherical RBC. Using our DNS data, we validate the ratio of inner and outer boundary layer thicknesses as well as the mean temperature drop across thermal boundary layers given by Gastine et al. (2015). Building upon these relations, we derive an analytical relation for the radius ratio dependence of $Nu$. To further investigate the dependency of $Nu$ and $Re$ on $Ra$ for different radius ratios, we employ the Grossmann–Lohse (GL) theory (Grossmann & Lohse 2000, 2001; Gastine et al. 2015). GL theory is grounded on the assumption of laminar BLs and turbulent bulk flows. The thermal and kinetic energy dissipation rates are partitioned into boundary layer and bulk contributions. Scaling relations are established between dissipation rates and the global Reynolds number. Using the GL theory, we are able to elucidate the dependence of the local scaling exponents on $\eta$. In addition to the scaling behaviour of the dimensionless heat transport and flow velocity, we also aim to quantify the asymmetry of the convective flow by partitioning the spherical shell domain into inner and outer regions, with local Reynolds numbers defined for each segment. The ratio between the inner and outer shell Reynolds numbers is expressed as a function of $\eta$. Lastly, we explore the asymmetry of the BL dissipation rates – with the ratio between inner and outer BL dissipation rates also expressed as a function of $\eta$. All scaling relations are validated using our DNS data. The major contributions of the present study are as follows.

1. (i) The large-scale structures are observed in spherical RBC for the first time.

2. (ii) An analytical relation for the radius ratio dependence of $Nu$ is derived.

3. (iii) The relations for local scaling exponents for $Nu(Ra)$ and $Re(Ra)$ are investigated by using GL theory.

4. (iv) The asymmetry of the convective flow and the boundary layer dissipation rates are quantified as a function of radius ratio.

2. Model description

2.1. Governing equations

We consider RBC with Oberbeck–Boussinesq approximation in spherical shells with inner radius $r_{i}$ and outer radius $r_{o}$. The temperatures are fixed as $T_i$ and $T_o$ at the inner and the outer boundary, respectively. The mechanical boundary conditions are no-slip at both boundaries. The dimensionless equations were adopted by using the shell gap $d=r_{o}-r_{i}$, the viscous dissipation time scale $d^2/\nu$, the temperature difference between the inner and outer boundaries $\Delta T=T_i-T_o$, the momentum diffusive velocity $\nu /d$, and the characteristic pressure $\rho \nu ^2 /d^2$ as the reference scales. Gravity is non-dimensionalized using its reference value at the outer boundary $g_o=g(r_o)$. The dimensionless governing equations for the problem read

(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 + \frac{Ra}{Pr}gT\boldsymbol{e_r}+\nabla^2 \boldsymbol{u}, \end{gather}$$

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

Here, the symbols $\boldsymbol {u}$, $p$ and $T$ represent velocity, pressure and temperature, respectively. Additionally, $\boldsymbol {e_r}$ represents the unit vector in the radial direction. In our study, we focus on spherical shell convection under the influence of a centrally condensed mass with the gravity profile, $g(r) = (r_{o}/r)^2$ (Chandrasekhar 1981). By adopting this gravity profile, exact relations connecting dissipation rates with driving forces can be established (see § 2.3), enabling us to conduct scaling analysis effectively.

The dimensionless equations (2.1)–(2.3) are controlled by three input parameters – Rayleigh number, Prandtl number and the radius ratio, which are defined as

(2.4a–c)\begin{equation} Ra=\frac{\alpha g_o \Delta T d^3}{\nu \kappa}, \quad Pr=\frac{\nu}{\kappa}, \quad \eta=\frac{r_i}{r_o}, \end{equation}

respectively. Here, $\alpha$ is the thermal expansion coefficient, and $\nu$ and $\kappa$ are the viscous and thermal diffusivities, respectively.

2.2. Response parameters

In RB, there are two key response parameters in the system. One of them is the Nusselt number, $Nu$, which represents the dimensionless heat flux. Within the Oberbeck–Boussinesq approximation, one obtains

(2.5) \begin{equation} Nu = \frac{\overline{\langle u_rT \rangle_s}-\dfrac{1}{Pr}\dfrac{{\rm d} \vartheta}{{\rm d} r}}{-\dfrac{1}{Pr} \dfrac{{\rm d} T_c}{{\rm d} r}} ={-}\eta \frac{{\rm d} \vartheta}{{\rm d} r} \Bigg|_{r=r_i} ={-}\frac{1}{\eta} \frac{{\rm d} \vartheta}{{\rm d} r} \Bigg|_{r=r_o},\end{equation}

where $\overline {({\cdot })}$ represents the time average, and $\langle {\cdot } \rangle _s$ and $\langle {\cdot } \rangle$ represent the spatial average over the horizontal surface and the whole volume of the spherical shell, respectively. The time and horizontally averaged radial temperature profile is $\vartheta =\overline {\langle T \rangle _s}$, and $T_c$ is the conductive temperature profile which, for spherical shells with fixed thermal boundary condition, reads

(2.6)\begin{equation} T_c(r) = \frac{\eta}{(1-\eta)^2} \frac{1}{r} - \frac{\eta}{1-\eta}. \end{equation}

Another key response parameter is the Reynolds number,

(2.7)\begin{equation} Re = \overline{\sqrt{\langle u^2 \rangle}} = \overline{\sqrt{\langle u_r^2+u_{\theta}^2+u_{\phi}^2 \rangle}},\end{equation}

which represents the dimensionless velocity. The radial profile of the time and horizontally averaged horizontal velocity is

(2.8)\begin{equation} Re_h(r)=\overline{\sqrt{\langle u_{\theta}^2+u_{\phi}^2 \rangle}_s}. \end{equation}

2.3. Exact dissipation relations in spherical shells

There are two exact relations for the time- and volume-averaged kinetic energy dissipation rate $\epsilon _u$ and the thermal energy dissipation rate $\epsilon _{\vartheta }$. In spherical shells with a centrally condensed mass gravity profile $g(r) = (r_{o}/r)^2$, the exact relations can be written as (Gastine et al. 2015)

(2.9)$$\begin{gather} \epsilon_u = \frac{3}{1+\eta+\eta^2} \frac{Ra}{Pr^2} (Nu-1), \end{gather}$$

(2.10)$$\begin{gather}\epsilon_{\vartheta} = \frac{3\eta}{1+\eta+\eta^2} Nu. \end{gather}$$

In our dimensionless quantities, $\epsilon _u$ and $\epsilon _{\vartheta }$ are calculated as follows:

(2.11a,b) \begin{equation} \epsilon_u = \overline{\biggl\langle \left(\frac{\partial u_i}{\partial x_j} \right)^2 \biggr\rangle}, \quad \epsilon_{\vartheta} = \overline{\langle (\boldsymbol{\nabla} T )^2 \rangle}.\end{equation}

Using the relations (2.9), (2.10) and (2.11a,b), we can express the Nusselt numbers based on the kinetic and thermal energy dissipation rates as

(2.12) \begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle Nu_{\epsilon_u} = \dfrac{1+\eta+\eta^2}{3} \dfrac{Pr^2}{Ra} \epsilon_u +1,\\ \displaystyle Nu_{\epsilon_{\vartheta}} = \dfrac{1+\eta+\eta^2}{3\eta} \epsilon_{\vartheta}. \end{array}\right\}\end{equation}

The relative error between the estimates of the Nusselt number obtained from (2.5) and (2.12) is used as a measure of the resolution in our simulations. For more details, please refer to Appendix A.

2.4. Numerical settings and simulation parameters

In this work, we have used the magnetohydrodynamic code MagIC (Wicht 2002; Christensen & Wicht 2007; Lago et al. 2021) to solve (2.1)–(2.3). MagIC is a pseudo-spectral code in which all the unknown variables are expanded into complete sets of functions in radial and horizontal directions. Chebyshev polynomials are applied in the radial direction while spherical harmonic functions are used in the azimuthal and latitudinal directions. Since the velocity field is solenoidal under the Oberbeck–Boussinesq approximation, MagIC decomposes it into poloidal and toroidal components,

(2.13)\begin{equation} \boldsymbol{u} = \boldsymbol{\nabla} \times (\boldsymbol{\nabla} \times W \boldsymbol{e_r}) + \boldsymbol{\nabla} \times Z \boldsymbol{e_r}. \end{equation}

The velocity field $\boldsymbol {u}$, which has three components, can thus be replaced by two scalar fields which are the poloidal potential $W$ and the toroidal potential $Z$. The equations are time-stepped by advancing the nonlinear terms using an explicit second-order Adams–Bashforth scheme, while the linear terms are time-advanced using an implicit Crank–Nicolson algorithm. At each time step, the linear terms are calculated in the spectral space while the nonlinear terms are calculated in the physical space.

We conducted six sets of simulations with $\eta = 0.2$, $0.3$, $0.4$, $0.5$, $0.6$, $0.8$ and $Pr=1$. For $\eta =0.2, 0.3$ and $0.4$ simulations, $Ra$ is varied in the range $3 \times 10^3 \leq Ra \leq 5 \times 10^8$; for the $\eta =0.5$ and $0.8$ cases, we have $3 \times 10^{5} \leq Ra \leq 5 \times 10^8$. The simulations for the $\eta =0.6$ case are the same as those of Gastine et al. (2015). For more details about the simulation parameters, grid resolution and the balance of the turbulent kinetic energy budget, please refer to Appendix A.

3. Flow structures

A typical representation of the differences in the flow morphology with respect to $\eta$ is shown in figure 1. Three-dimensional contours of the temperature fluctuations $T'$ for a small radius ratio case ($\eta =0.2$) and a large radius ratio case ($\eta =0.8$) are plotted. The radial cuts are located inside the thermal BLs. There are many ways to define the thermal boundary layer (Long et al. 2020). In this work, we have adopted the slope method (Shishkina et al. 2010; Gastine et al. 2015) to define the thermal boundary layer thickness. It can be clearly seen that the flow structures mainly comprise sheet-like plumes originating from within the thermal BLs. The hot plumes (red) detach from the inner boundary and then expand in the middle of the domain, while the cold plumes (blue) detach from the outer boundary and expand, sinking downwards towards the inner shell. The rising hot plumes and the sinking cold plumes evolve into mushroom-type lobed structures in the bulk. The plume size near the inner boundary is smaller than that of the outer boundary, which is in line with the observation reported by Gastine et al. (2015). By comparing the plume morphology between different radius ratio cases, keeping $Ra$ constant, we find that the number of the plumes ejected from the boundaries is higher at a larger $\eta$. See figure 1. Moreover, the number of the sheet-like plumes is higher near the outer boundary than near the inner boundary. This phenomenon is more pronounced at a smaller $\eta$. For $\eta =0.2 \textrm { and } Ra=3 \times 10^{8}$, flow structures characterized by the normalized vertical velocity $u_{r}'=u_{r}/(\sqrt {Ra/Pr})$ as well as the temperature fluctuations $T^{'}=T-\overline { \langle T \rangle _{s}}$ at three different horizontal depths are shown in figure 2. These three layers are chosen to represent the mid-depth, a layer near the inner thermal boundary layer and a layer near the outer thermal boundary layer. We can see from both $u_{r}'$ and $T^{'}$ contour plots that the plume number is higher near the outer boundary than near the inner boundary.

Figure 1. Equatorial and meridional cuts of temperature fluctuations $T'$ for two selected cases. (a) $\eta =0.2, Ra=7 \times 10^{6}$, the inner radial cut is at $r=r_{i}+0.01d$ and the outer radial cut is at $r=r_{o}-0.02d$. (b) $\eta =0.8, Ra=7 \times 10^{6}$, the inner radial cut is at $r=r_{i}+0.02d$ and the outer radial cut is at $r=r_{o}-0.02d$. Colour levels range from $-0.3$ (blue) to $0.3$ (red).

Figure 2. Normalized vertical velocity $u_{r}'=u_{r}/(\sqrt {Ra/Pr})$ and temperature fluctuations $T^{'}=T-\overline { \langle T \rangle _{s}}$ at different horizontal layers for $\eta =0.2, Ra=3 \times 10^{8}$. (a) $u_{r}'$ at mid-depth, (b) $T^{'}$ at mid-depth, (c) $u_{r}'$ near the inner boundary layer at $r=r_{i}+0.0047d$, (d) $T^{'}$ near the inner boundary layer at $r=r_{i}+0.0047d$, (e) $u_{r}'$ near the outer boundary layer at $r=r_{i}+0.9770d$, (f) $T^{'}$ near the outer boundary layer at $r=r_{i}+0.9770d$.

The warm upflow (red) and cold downflow (blue) reveal the existence of the large-scale structures in the flow. Figure 3 shows normalized instantaneous vertical velocity $u_{r}'=u_{r}/(\sqrt {Ra/Pr})$ on the horizontal mid-plane at $Ra=7 \times 10^7$ and different $\eta$. As illustrated in figure 3, the higher the radius ratio, the greater the number of large-scale structures present in the flow. To reveal the dominant spherical harmonic degree (analogous to the dominant wavenumber in Cartesian coordinates) at each $\eta$ in figure 3, time-averaged kinetic energy spectra with respect to the spherical harmonic $l$ on the mid-plane are calculated as (Glatzmaier 2013)

(3.1)\begin{equation} E_{kin}(l) =\frac{1}{2 {\rm \pi}r_{mid}^2 \overline{\langle \boldsymbol{| u |}^2 \rangle_{s}}} \sum_{m=0}^{l} l(l+1) \left[ \frac{l(l+1)}{r_{mid}^2} \overline{| W_{lm} |^2} + \overline{\left| \frac{{\rm d} W_{lm}}{{\rm d} r} \right|^2} + \overline{| Z_{lm} |^2}\right], \end{equation}

where $r_{mid}=(r_i+r_o)/2$ is the radius of the mid-plane, $m$ is the spherical harmonic order, and $W_{lm}$ and $Z_{lm}$ are the poloidal potential and the toroidal potential of the velocity in the spectral space, respectively (Christensen & Wicht 2007). In the above equation, the $m=0$ contribution entering in the summation is multiplied by one-half (Schwaiger, Gastine & Aubert 2021), see Appendix B for more details. Figure 4 shows kinetic energy spectra $E_{kin}(l)$ for the same cases as depicted in figure 3. The dominant harmonic degree $l_{max}$ corresponds to the peak of the spectra. A higher value of $l_{max}$ implies a greater number of large-scale structures. As seen from figure 4, at $Ra=7 \times 10^{7}$, the dominant harmonic degree varies monotonically from $l_{max}=1$ for $\eta =0.2$ to $l_{max}=10$ for $\eta =0.8$. This again indicates that the number of large-scale structures increases with $\eta$, consistent with what is observed in figure 3. Furthermore, although not plotted here, we find that this observation holds for any fixed $Ra$.

Figure 3. Hammer projection of the normalized instantaneous vertical velocity $u_{r}'=u_{r}/(\sqrt {Ra/Pr})$ at mid-depth $r_{mid}=(r_i+r_o)/2$. All the plots shown are at $Ra=7 \times 10^7$. The radius ratios vary as (a) $\eta =0.2$, (b) $\eta =0.3$, (c) $\eta =0.4$, (d) $\eta =0.5$, (e) $\eta =0.6$ and (f) $\eta =0.8$.

Figure 4. Kinetic energy spectra $E_{kin}(l)$ with respect to the harmonic degree $l$ at mid-depth $r_{mid}=(r_i+r_o)/2$. All the spectra shown are for $Ra=7 \times 10^7$ at different $\eta$. The radius ratios vary as (a) $\eta =0.2$, (b) $\eta =0.3$, (c) $\eta =0.4$, (d) $\eta =0.5$, (e) $\eta =0.6$ and (f) $\eta =0.8$. Here, $l_{max}$ is the harmonic degree corresponding to the peak of the spectrum.

Moreover, it is noted that $l_{max}$ decreases with increasing $Ra$ at a constant $\eta$, suggesting a reduction in the number of large-scale structures at higher $Ra$. This trend is illustrated in figure 5, which displays $u_{r}'$ on the horizontal mid-plane for $\eta =0.2$ at various $Ra$ values, and figure 6, showing the corresponding kinetic energy spectra. We can clearly see from figure 6 that the dominant harmonic degree varies from $l_{max}=2$ at $Ra=7 \times 10^{5}$ to $l_{max}=1$ at $Ra=3 \times 10^{8}$. Similar inferences can be made from figures 7 and 8 for $\eta =0.8$, where the dominant harmonic degree varies monotonically from $l_ {max}=12$ at $Ra=7 \times 10^{5}$ to $l_{max}=9$ at $Ra=3 \times 10^{8}$. The decrease in $l_{max}$ with rising $Ra$ is evident from these plots. Since $l_{max}$ corresponds to the dominant length scale as shown in (3.2), this phenomenon indicates that at a fixed $\eta$, the increase of $Ra$ relates to a decrease of $l_{max}$, which in turn translates to an increase of the dominant length scale $L_{dom}$.

Figure 5. Hammer projection of the normalized instantaneous vertical velocity $u_{r}'=u_{r}/(\sqrt {Ra/Pr})$ on the horizontal mid-plane $r_{mid}=(r_i+r_o)/2$ for $\eta =0.2$. $Ra$ varies as (a) $7 \times 10^5$, (b) $3 \times 10^6$, (c) $7 \times 10^6$, (d) $3\times 10^7$, (e) $7 \times 10^7$ and (f) $3 \times 10^8$.

Figure 6. Kinetic energy spectra $E_{kin}(l)$ with respect to harmonic degree $l$ on the horizontal mid-plane $r_{mid}=(r_i+r_o)/2$ for $\eta =0.2$. $Ra$ varies as (a) $7 \times 10^5$, (b) $3 \times 10^6$, (c) $7 \times 10^6$, (d) $3\times 10^7$, (e) $7 \times 10^7$ and (f) $3 \times 10^8$. Here, $l_{max}$ is the harmonic degree corresponding to the peak of the spectrum.

Figure 7. Hammer projection of the normalized instantaneous vertical velocity $u_{r}'=u_{r}/(\sqrt {Ra/Pr})$ on the horizontal mid-plane $r_{mid}=(r_i+r_o)/2$ for $\eta =0.8$. $Ra$ varies as (a) $7 \times 10^5$, (b) $3 \times 10^6$, (c) $7 \times 10^6$, (d) $3\times 10^7$, (e) $7 \times 10^7$ and (f) $3 \times 10^8$.

Figure 8. Kinetic energy spectra $E_{kin}(l)$ with respect to harmonic degree $l$ on the horizontal mid-plane $r_{mid}=(r_i+r_o)/2$ for $\eta =0.8$. The $Ra$ varies as (a) $7 \times 10^5$, (b) $3 \times 10^6$, (c) $7 \times 10^6$, (d) $3\times 10^7$, (e) $7 \times 10^7$ and (f) $3 \times 10^8$. Here, $l_{max}$ is the harmonic degree corresponding to the peak of the spectrum.

The degree $l_{max}$ is associated with the dominant length scale $L_{dom}$ by the definition of the characteristic wavelength (Backus, Parker & Constable 1996)

(3.2)\begin{equation} L_{dom} = \frac{2 {\rm \pi}r_{mid}}{\sqrt{l_{max}(l_{max}+1)}}, \end{equation}

which is associated with the length scale of the large-scale structures (Stevens et al. 2018). For comparison, we also look at the so-called integral length scale (Parodi et al. 2004)

(3.3)\begin{equation} L_{int} = \sum_{l} \frac{2 {\rm \pi}r_{mid}}{\sqrt{l(l+1)}} E_{kin}(l). \end{equation}

Figure 9 illustrates $L_{dom}$ and $L_{int}$ for $\eta =0.2$ and $0.8$ at various $Ra$. The integral length scale stays relatively constant at approximately 1.5 across different $Ra$ values. In planar RBC, Stevens et al. (2018) demonstrated that the integral length scale derived from kinetic energy spectra saturates at large aspect ratios ($\varGamma \geq 32$) to a value of $\sim$1.6.

Figure 9. (a) Integral ($L_{int}$) and dominant ($L_{dom}$) length scales with respect to $Ra$ for $\eta =0.2$ and $\eta =0.8$. (b) Integral length scale $L_{int}$ with respect to $\eta$ at different $Ra$.

As stated above, in spherical RBC, it is observed that the integral length scale saturates around the value of $1.5$ when $\eta \geq 0.2$, as shown in figure 9(b). This convergence value of $L_{int}$ in spherical RBC is very close to that observed in the planar case. We provide a detailed, qualitative explanation on this observation in Appendix C. Moreover, the dominant length scale increases with increasing $Ra$ for both $\eta =0.2$ and $\eta =0.8$, indicating a growth in the scale of large-scale structures with increasing $Ra$. Notably, data for $\eta =0.2$ and $\eta =0.8$ exhibit a noticeable convergence of $L_{dom}$ around $L_{dom} \approx 3.3$. At $Ra=3 \times 10^{8}$, the dominant length scales for both $\eta =0.2$ and $0.8$ closely resemble those observed by Stevens et al. (2018) at $Ra=10^{8}$ and aspect ratios larger than 32. Additionally, we found that the large-scale structures in our simulations do not vary with time. As examples, figures 10 and 11 show the time evolution of the temperature fluctuation $T^{'}=T-\overline { \langle T \rangle _{s}}$ at $\eta =0.2$ and $\eta =0.8$, respectively. The shape and position of the large-scale structures remain unchanged for time intervals exceeding $150$ convective time units. In fact, we examined our simulations for more than $500$ convective time units and observed no variation in the large-scale structures. Furthermore, a small set of additional simulations revealed that the number of large-scale structures is independent of the initial conditions. Multiple states, as have been observed in two-dimensional simulations (Wang et al. 2020), were not observed in our 3-D cases.

Figure 10. Time series of the temperature fluctuations $T^{'}=T-\overline { \langle T \rangle _{s}}$ at the mid-depth for $\eta =0.2$ at different time. Here, $t_{conv} = \sqrt {d/(g_{0} \alpha \Delta T)}$ is the convective time unit. (a) $t_0$, (b) $t_0+45.8t_{conv}$, (c) $t_0+95.5t_{conv}$ and (d) $t_0+165.1t_{conv}$.

Figure 11. Time series of the temperature fluctuations $T^{'}=T-\overline { \langle T \rangle _{s}}$ at the mid-depth for $\eta =0.8$ at different time. Here, $t_{conv} = \sqrt {d/(g_{0} \alpha \Delta T)}$ represents for the convective time unit. (a) $t^{'}_0$, (b) $t^{'}_0+54.7t_{conv}$, (c) $t^{'}_0+109.3t_{conv}$ and (d) $t^{'}_0+173.1t_{conv}$.

4. Global scaling laws of response parameters

4.1. Asymmetry of temperature and velocity profiles

RBC in spherical shell geometry is significantly different from its planar counterpart, chiefly due to the curvature of the bounding surfaces. The presence of curvature in spherical RBC leads to asymmetry of the radial profiles of temperature and velocity fluctuations. Owing to the thermal shortcut (Ahlers et al. 2009), the temperature drop mostly happens inside the thermal BLs. As a result, the properties of BLs are essential, especially in the quantification of the scaling relations of the global response parameters. In this work, we use the slope method (Shishkina et al. 2010) to define the kinetic as well as the thermal BL thicknesses for consistency with the BL analyses of Prandtl (1905) and Blasius (1907). The definitions of the boundary layer thicknesses are the same as those of Gastine et al. (2015).

Time and horizontally averaged radial temperature profiles $\vartheta (r)$ exhibit a strong asymmetry in spherical shell convection. Figure 12 shows $\vartheta (r)$ and $Re(r)$ at different radius ratios. At small radius ratios, the differences in curvature and gravity between the inner and outer boundaries are significantly larger. This results in a stronger asymmetry in both the temperature as well as the velocity profiles. For example, as illustrated in figure 12(a), the inner BL temperature drop is approximately $57\,\%$ of the total temperature drop at $\eta =0.8, Ra=7 \times 10^{7}$. However, at $\eta =0.2, Ra=7 \times 10^{7}$, the inner BL temperature drop constitutes approximately $90\,\%$ of the total temperature drop, implying a stronger asymmetry in the temperature field. For the velocity field, represented by $Re(r)$, the profile is almost symmetrical when $\eta =0.8$. As $\eta$ becomes smaller, the difference between $Re(r)$ in the vicinity of the inner and the outer boundaries becomes larger. At $\eta =0.2$, the peak value of $Re(r)$ near the inner boundary is approximately twice that of the value near the outer boundary, as shown in figure 12(b). The asymmetry of the velocity field also illustrates that the flow is more turbulent in the inner half-shell than in the outer half-shell.

Figure 12. Radial profiles of time and horizontally averaged (a) temperature $\vartheta (r)$ and (b) Reynolds number $Re(r)$, at different $\eta$. All the simulations presented here are at $Ra=7 \times 10^7$ and $Pr=1$.

For the temperature field, the asymmetry is quantified by the temperature drop $\Delta \vartheta _{i}$ ($\Delta \vartheta _{o}$) and the thermal BL thickness $\lambda _{\vartheta }^{i}$ ($\lambda _{\vartheta }^{o}$) at the inner (outer) BL. We assume that the temperature drop occurs only inside the BLs, which leads to

(4.1)\begin{equation} \Delta \vartheta_{i} + \Delta \vartheta_{o} = 1,\end{equation}

so that we have

(4.2a,b)\begin{equation} \Delta \vartheta_{i} = 1 - \vartheta_{m}, \quad \Delta \vartheta_{o} = \vartheta_{m}, \end{equation}

where $\vartheta _{m}=\vartheta ((r_i+r_o)/2)$ is the time and horizontally averaged temperature at the mid-depth. Gastine et al. (2015) gives an analytical expression for $\Delta \vartheta _{i}, \Delta \vartheta _{o}$ and $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }^{o}$. For the present paper to be self-contained, an outline of the scaling arguments used by Gastine et al. (2015) is presented below. From (2.5), we obtain

(4.3)\begin{equation} \frac{\Delta \vartheta_{i}}{\lambda_{\vartheta}^{i}} = \frac{1}{\eta^2} \frac{\Delta \vartheta_{o}}{\lambda_{\vartheta}^{o}}.\end{equation}

Gastine et al. (2015) proposed that the average plume density in the inner and outer BLs should be the same, which yields

(4.4)\begin{equation} \rho_{p}^{i} = \rho_{p}^{o} \rightarrow \frac{\alpha g_{i} \Delta \vartheta_{i} {\lambda_{\vartheta}^{i}}^{5}}{\nu \kappa} = \frac{\alpha g_{o} \Delta \vartheta_{o} {\lambda_{\vartheta}^{o}}^{5}}{\nu \kappa},\end{equation}

where $\rho _{p}^{i}$ ($\rho _{p}^{o}$) is the average plume density in the inner (outer) boundary, $g_{i}(g_{o})$ being the corresponding values of gravitational acceleration. Combining (4.1)–(4.3) and (4.4), Gastine et al. (2015) finally derive

(4.5a–c)\begin{equation} \Delta \vartheta_{i} = \frac{1}{1+\eta^{5/3}\chi_{g}^{1/6}}, \quad \Delta \vartheta_{o} = \vartheta_{m} = \frac{\eta^{5/3} \chi_{g}^{1/6}}{1+\eta^{5/3} \chi_{g}^{1/6}}, \quad \frac{\lambda_{\vartheta}^{i}}{\lambda_{\vartheta}^{o}} = \frac{\eta^{1/3}}{\chi_{g}^{1/6}}, \end{equation}

where $\chi _{g}=g(r_i)/g(r_o)$ denotes the ratio of the gravitational acceleration between the inner and the outer boundary. For the $g(r)$ profile considered in the present work, $\chi _{g}=1/\eta ^{2}$. Hence, (4.5a–c) can be written as

(4.6a–c)\begin{equation} \Delta \vartheta_{i} = \frac{1}{1+\eta^{4/3}}, \quad \Delta \vartheta_{o} = \vartheta_{m} = \frac{\eta^{4/3}}{1+\eta^{4/3}}, \quad \frac{\lambda_{\vartheta}^{i}}{\lambda_{\vartheta}^{o}} = \eta^{2/3}. \end{equation}

However, there are some other models for $\Delta \vartheta _{i}, \Delta \vartheta _{o}$ and $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }$ (e.g. Wu & Libchaber 1991; Zhang, Childress & Libchaber 1997; Wang et al. 2022). Gastine et al. (2015) shows that the model specified in (4.5a–c) is better than the others when $Pr=1$. Thus, we employ (4.5a–c) and (4.6a–c) in our work. Defining the asymmetry factor $\chi$ by employing the temperature drop at the inner BL,

(4.7)\begin{equation} \chi = \frac{2 \Delta \vartheta_{i}}{\Delta \vartheta} = \frac{2}{1+\eta^{4/3}},\end{equation}

where $\Delta \vartheta$ is the total temperature drop across the spherical shell. With $Pr=1$, the ratio of the inner and the outer kinetic BL thickness should be the same as that for thermal BL thicknesses (Shishkina et al. 2010), such that

(4.8)\begin{equation} \frac{\lambda_{u}^{i}}{\lambda_{u}^{o}} = \frac{\lambda_{\vartheta}^{i}}{\lambda_{\vartheta}^{o}} = \eta^{2/3}. \end{equation}

In figure 13, (4.7) and (4.8) are successfully validated against our DNS data.

Figure 13. (a) Asymmetry factor $\chi$ and (b) ratio of the inner and outer BL thicknesses, with respect to $\eta$. Symbols represent DNS data. Dashed curves represent analytical predictions given by (4.7) and (4.8).

From (4.6a–c) and the $Nu$ definition in (2.5), we can directly obtain the following relations:

(4.9a,b)\begin{equation} \lambda_{\vartheta}^{i} = \frac{\eta}{1+\eta^{4/3}} \frac{1}{Nu}, \quad \lambda_{\vartheta}^{o} = \frac{\eta^{1/3}}{1+\eta^{4/3}} \frac{1}{Nu}. \end{equation}

Thus, we can define a normalized thermal boundary layer thickness

(4.10)\begin{equation} \tilde{\lambda}_{\vartheta} = \frac{1+\eta^{4/3}}{\eta} \lambda_{\vartheta}^{i} = \frac{1+\eta^{4/3}}{\eta^{1/3}} \lambda_{\vartheta}^{o} = \frac{1}{Nu}.\end{equation}

Figure 14 shows a validation of (4.10) and it shows that the relation perfectly aligns with our DNS data.

Figure 14. Normalized thermal boundary layer thickness $\tilde {\lambda }_{\vartheta }$ as a function of $Nu$ at different $\eta$. Symbols represent DNS data, while the dashed black line denotes $\tilde {\lambda }_{\vartheta } = 1/Nu$.

4.2. Radius ratio dependence of Nusselt number

In spherical RBC, the geometrical asymmetry present in the system is expected to have an effect on the scaling law of the dimensionless heat transport. Figure 15(a) shows $Nu$ as a function of $Ra$ at different $\eta$. In the considered $Ra$ range, the smaller $\eta$ cases have a smaller value of $Nu$ at a fixed $Ra$. As $\eta \rightarrow 1$, $Nu$ scaling will asymptotically approach the planar case with an infinite aspect ratio and the differences in the exponent become smaller at larger $\eta$; for $\eta \geq 0.4$, the differences in $Nu$ at different $\eta$ are vanishingly small and almost indistinguishable in the plot. Moreover, $Ra$ dependence of $Nu$ also varies with $\eta$. For example, the data fit for $Nu$ at $\eta =0.2$ is $Nu |_{\eta =0.2}=0.101 Ra^{0.308}$, while as for the $\eta =0.8$ case, it reads $Nu |_{\eta =0.8}=0.156 Ra^{0.291}$. Taking into consideration the effect of $\eta$ on $Nu$ as well as the $\eta$-dependence of the scaling exponent with respect to $Ra$, a power law of the form $Nu=f(\eta ) Ra^{\gamma (\eta,Ra)}$ is deemed appropriate (Ahlers et al. 2009).

Figure 15. (a) Nusselt number, $Nu$, and (b) the compensated Nusselt number, $Nu (1+\eta ^{4/3})/(\eta ^{1/2})$, as functions of $Ra$ at different $\eta$. Symbols represent DNS data and the dashed black line corresponds to the equation given by (4.14).

Now we come to the scaling for $Nu$ with radius ratio dependence. It is well known that the width of the BL is related to $Nu$. Equation (2.5) can be written as

(4.11)\begin{equation} Nu ={-}\eta \frac{{\rm d} \vartheta}{{\rm d} r} |_{r=r_i} \approx \eta \frac{\Delta \vartheta_{i}}{\lambda_{\vartheta}^{i}}.\end{equation}

Equation (4.6a–c) gives us the $\eta$-dependence of the temperature drop in the inner BL (Gastine et al. 2015), while the radius ratio dependence of $\lambda _{\vartheta }^{i}$ remains to be found.

Within the BL, the viscous term balances the advection term in (2.2), which leads to

(4.12)\begin{equation} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} \sim \nu \nabla^2 \boldsymbol{u} \Rightarrow \frac{U}{L} \sim \nu \frac{1}{\delta^2},\end{equation}

where $U$ is the velocity of the large-scale circulation, $L$ is the associated length scale and $\delta$ represents the length scale relevant inside the BL. Inside the inner BL, if we assume $L=L_{i} \sim (r_{i}/r_{o}) d$ and $\delta = \tilde {\lambda }_{u}^{i} \approx \tilde {\lambda }_{\vartheta }^{i}$, we have

(4.13)\begin{equation} \tilde{\lambda}_{\vartheta}^{i} \sim \left( \frac{\nu d}{U} \right)^{1/2} \eta^{1/2} \Rightarrow \lambda_{\vartheta}^{i} = \frac{\tilde{\lambda}_{\vartheta}^{i}}{{\rm d}} \sim Re^{1/2} \eta^{1/2},\end{equation}

where $\lambda _{u}^{i}, \lambda _{\vartheta }^{i}$ are the normalized inner viscous and thermal BL thicknesses, respectively; whereas $\tilde {\lambda }_{u}^{i}, \tilde {\lambda }_{\vartheta }^{i}$ represent their dimensional counterparts. The shell thickness $d=r_o-r_i$ is used to non-dimensionalize $\tilde {\lambda }_{u}^{i}$ and $\tilde {\lambda }_{\vartheta }^{i}$. Here, $Re=Ud/\nu$ is the global Reynolds number as defined in (2.7). In our simulations, it is observed that $Re$ has only a weak dependence on $\eta$, as shown in figure 16. Hence, it is reasonable to regard $Re$ as being independent of the radius ratio. We will investigate the radius ratio dependence of $Re$ more precisely by applying the GL theory (Grossmann & Lohse 2000) in subsequent sections.

Figure 16. Reynolds number $Re$ as a function of the Rayleigh number $Ra$ at different $\eta$. Symbols represent DNS data, while the dashed black line denotes $Re=0.34Ra^{0.48}$.

Figure 17 shows that the DNS data collapse well with $\lambda _{\vartheta }^{i} \eta ^{-1/2} = 3.3 Ra^{-0.295}$, which validates the result that $\lambda _{\vartheta }^{i}$ is proportional to $\eta ^{1/2}$. Finally, using the $\eta$-dependence of $\Delta \vartheta _{i}$ and $\lambda _{\vartheta }^{i}$ from (4.6a–c) and (4.13), respectively, in (4.11), we get the following scaling relation for the Nusselt number:

(4.14)\begin{equation} Nu \approx \eta \frac{\Delta \vartheta_{i}}{\lambda_{\vartheta}^{i}} \sim \frac{\eta^{1/2}}{1+\eta^{4/3}} Ra^{\gamma}.\end{equation}

Here, $\gamma \approx 0.295$ if we fit a single scaling exponent in the whole data range. A similar exponent has also been reported by Brown et al. (2005) and Funfschilling et al. (2005) in the cylindrical geometries. Figure 15 also shows that (4.14) exhibits a good agreement with our DNS data.

Figure 17. (a) Inner thermal BL thickness $\lambda _{\vartheta }^{i}$ and (b) compensated inner BL thickness $\lambda _{\vartheta }^{i} \eta ^{-1/2}$, as functions of $Ra$ at different $\eta$. Symbols represent DNS data, while the dashed black line represents the fit $\lambda _{\vartheta }^{i} \eta ^{-1/2} = 3.3 Ra^{-0.295}$, which follows from the scaling argument (4.13).

4.3. Dissipation analysis

It must be stressed that in (4.14), although $\gamma \approx 0.295$ agrees well with the DNS data if fitting with one single scaling exponent in the whole $Ra$ range, the local effective heat transfer scaling exponent varies with $Ra$ (Grossmann & Lohse 2000, 2001). To investigate $\eta$ and $Ra$ dependence of $Nu$ and $Re$ in more detail, we apply the GL theory to our dataset. There are basically two key assumptions in GL theory (Grossmann & Lohse 2000, 2001; Ahlers et al. 2009). One is that there exists a large-scale circulation (LSC) with only one typical velocity scale $U$. Another is that the kinetic BLs are scaling-wise characterized by a single effective thickness $\lambda _{u}$ regardless of the positions along the boundaries. For the first assumption, we observe that there exist large-scale structures in our simulations, as demonstrated in § 3. These large-scale structures can be regarded as indicative of the LSC. Moreover, we use the time- and volume-averaged velocity (2.7) as the typical velocity scale $U$, and use the $Re$ based on this velocity scale to derive the scaling relations for $\epsilon _{u}^{bulk}(Re), \epsilon _{u}^{BL}(Re), \epsilon _{\vartheta }^{bulk}(Re)$ and $\epsilon _{\vartheta }^{BL}(Re)$. These scaling relations are validated by our DNS data, as will be shown in this section. For the assumption of kinetic BLs, Gastine et al. (2015) have already shown that the radial profile of the horizontal velocity looks very similar to the Prandtl–Blasius solution. We also validated that in our simulations. Hence, it would be reasonable to state that in the parameter space explored in our simulations, these two key assumptions of the GL theory are also valid in spherical shell RBC.

The basic idea of GL theory is to separate the whole flow into two parts, i.e. the laminar part and the turbulent part. The laminar part resides in the vicinity of the boundaries while the turbulent part concentrates in the bulk region. Following this idea, the kinetic and thermal energy dissipation rates can be separated into two contributions, which are the BL contribution and the bulk contribution. Therefore, the two dissipation rates can be expressed as

(4.15a,b)\begin{equation} \epsilon_u = \epsilon_{u}^{BL} + \epsilon_{u}^{bulk}, \quad \epsilon_{\vartheta} = \epsilon_{\vartheta}^{BL} + \epsilon_{\vartheta}^{bulk}.\end{equation}

Here, we use the subscripts $BL$ and $bulk$ for the BL and bulk contributions, respectively. The bulk and BL contributions are defined as follows:

(4.16a,b) $$\begin{gather} \epsilon_{u}^{bulk} = \frac{4 {\rm \pi}}{V} \int_{r_i+\lambda_{u}^{i}}^{r_o-\lambda_{u}^{o}} \overline{\biggl\langle \left(\frac{\partial u_i}{\partial x_j} \right)^2 \biggr\rangle}_{s} r^{2}\,\mathrm{d}r, \quad \epsilon_{\vartheta}^{bulk} = \frac{4 {\rm \pi}}{V} \int_{r_i+\lambda_{\vartheta}^{i}}^{r_o-\lambda_{\vartheta}^{o}} \overline{\langle (\Delta T )^2 \rangle}_{s} r^{2}\,\mathrm{d}r, \end{gather}$$

(4.17) $$\begin{gather} \left. \begin{array}{@{}l@{}l@{}} \displaystyle \epsilon_{u}^{BL} = \dfrac{4 {\rm \pi}}{V} \int_{r_i}^{r_i+\lambda_{u}^{i}} \overline{\biggl\langle \left(\dfrac{\partial u_i}{\partial x_j} \right)^2 \biggr\rangle}_{s} r^{2}\,\mathrm{d}r + \dfrac{4 {\rm \pi}}{V} \int_{r_o-\lambda_{u}^{o}}^{r_o} \overline{\biggl\langle \left(\dfrac{\partial u_i}{\partial x_j} \right)^2 \biggr\rangle}_{s} r^{2}\,\mathrm{d}r, \\ \displaystyle \epsilon_{\vartheta}^{BL} = \dfrac{4 {\rm \pi}}{V} \int_{r_i}^{r_i+\lambda_{\vartheta}^{i}} \overline{\langle (\Delta T )^2 \rangle}_{s} r^{2}\,\mathrm{d}r + \dfrac{4 {\rm \pi}}{V} \int_{r_o-\lambda_{\vartheta}^{o}}^{r_o} \overline{\langle (\Delta T )^2 \rangle}_{s} r^{2}\,\mathrm{d}r, \end{array}\right\} \end{gather}$$

where $V=(4/3) {\rm \pi}(r_{o}^3-r_{i}^3)$ is the volume of the spherical shell. Assuming the flow in the bulk region is highly turbulent, Grossmann & Lohse (2000) proposed the following scalings for $\epsilon _{u}^{bulk}$ and $\epsilon _{\vartheta }^{bulk}$ at a fixed $Pr$:

(4.18a,b) \begin{equation} \epsilon_{u}^{bulk} \sim Re^{3}, \quad \epsilon_{\vartheta}^{bulk} \sim Re. \end{equation}

Figures 18 and 19 show the bulk dissipation rates as a function of $Re$ at different $\eta$. The data fits for $Ra \geq 10^{5}$ yield $\epsilon _{u}^{bulk}$ scales between $\epsilon _{u}^{bulk} \sim Re^{2.72}$ (at $\eta =0.2$) and $\epsilon _{u}^{bulk} \sim Re^{2.80}$ (at $\eta =0.8$), and $\epsilon _{\vartheta }^{bulk}$ scales between $\epsilon _{\vartheta }^{bulk} \sim Re^{0.66}$ (at $\eta =0.2$) and $\epsilon _{\vartheta }^{bulk} \sim Re^{0.71}$ (at $\eta =0.6$). For both dissipation rates, our data show considerable deviations from the predictions of GL theory. These deviations can mainly occur because of two potential reasons: (1) the bulk regions are still not turbulent enough for the considered Reynolds number range, $80 < Re < 7000$; and (2) the bulk region defined in this paper is the region outside the BLs, which, as a result of plume ejection from the inner and outer boundaries, may contain plumes that are mostly laminar and can be considered as detached BLs (Grossmann & Lohse 2004). The presence of these plumes in the bulk region will decrease the scaling exponents of the bulk contributions in the two dissipation rates. Similar deviations have also been reported, for example, by Calzavarini et al. (2005).

Figure 18. (a) Bulk contribution of kinetic energy dissipation rate $\epsilon _{u}^{bulk}$ and (b) the compensated $\epsilon _{u}^{bulk} Re^{-3}$, as functions of $Re$. Symbols represent the DNS data and dashed lines correspond to the fitting relations.

Figure 19. (a) Bulk contribution of the thermal energy dissipation rate $\epsilon _{\vartheta }^{bulk}$ and (b) the compensated $\epsilon _{\vartheta }^{bulk} Re^{-1}$, as functions of $Re$. Symbols represent the DNS data and the dashed lines correspond to the fitting relations.

Based on the Prandtl–Blasius BL theory and the fact that the BL thickness is significantly smaller than the gap length, Grossmann & Lohse (2000) derived the following scaling relations for $\epsilon _{u}^{BL}$ and $\epsilon _{\vartheta }^{BL}$ at $Pr=1$,

(4.19a,b) \begin{equation} \epsilon_{u}^{BL} \sim Re^{5/2}, \quad \epsilon_{\vartheta}^{BL} \sim Re^{1/2}. \end{equation}

BL contributions for kinetic and thermal energy dissipation rates are shown in figures 20 and 21, respectively. As shown in figure 20, the best fit to DNS data with $Ra \geq 10^{5}$ yields $\epsilon _{u}^{BL} \sim Re^{2.54}$ (at $\eta =0.8$) and $\epsilon _{u}^{BL} \sim Re^{2.62}$ (at $\eta =0.2$). These fits are very close to the theoretical scaling of $\epsilon _{u}^{BL} \sim Re^{5/2}$. However, our data suggest the local scaling of $\epsilon _{u}^{BL}$ becomes larger with increasing $Ra$. It is even more clear in the compensated plot, as shown in figure 20(b). The local exponents of $\epsilon _{u}^{BL}(Re)$ initially decrease with $Re$ and then increase with $Re$, with the transitional $Re$ being approximately $500$. This observation reveals that a single power law is not sufficient for describing the $\epsilon _{u}^{BL}(Re)$ relation and there exists a secondary $Re$ dependence on $\epsilon _{u}^{BL}$. For $\epsilon _{\vartheta }^{BL}$ scaling, the best fit to DNS data with $Ra \geq 10^{5}$ shows $\epsilon _{\vartheta }^{BL}$ scales between $\epsilon _{\vartheta }^{BL} \sim Re^{0.57}$ (at $\eta =0.4$) and $\epsilon _{\vartheta }^{BL} \sim Re^{0.61}$ (at $\eta =0.2$), which are close to the expected exponents. It is observed that $\epsilon _{\vartheta }^{BL}$ scaling exhibits a similar behaviour to $\epsilon _{u}^{BL}$; that the local exponent becomes larger at a higher $Re$. Although this is less obvious for $\epsilon _{\vartheta }^{BL}$, we can still observe a gradual steepening of the slope when $Re$ increases, as shown in figure 21(b). This increase in the local exponent with respect to $Re$ points to the deviations from the simple power law predicted by GL theory for the BL dissipation rates. Similar deviations have also been reported by Gastine et al. (2015), as a result of the difficulty to perfectly separate the BL from the bulk.

Figure 20. (a) BL contribution of the kinetic energy dissipation rate $\epsilon _{u}^{BL}$ as a function of $Re$. Since the data fits at different $\eta$ are indistinguishable from each other, only the fit for the $\eta =0.8$ case is displayed on the plot. (b) Compensated $\epsilon _{u}^{BL} Re^{-5/2}$ as a function of $Re$. Symbols represent the DNS data and the dashed lines correspond to the fitting relations.

Figure 21. (a) BL contribution of the thermal energy dissipation rate $\epsilon _{\vartheta }^{BL}$ and (b) the compensated $\epsilon _{\vartheta }^{BL} Re^{-1/2}$, as functions of $Re$. Symbols represent the DNS data and the dashed lines correspond to the fitting relations.

Despite potential deviations of the local assumptions from the current DNS data, GL theory remains remarkably effective for our system across various radius ratios, as explained in detail below. Following GL theory, the two dissipation rates can be represented by the sum of bulk and BL contributions,

(4.20)\begin{equation} \left. \begin{array}{@{}c@{}} \hat{\epsilon}_{u} = \alpha_{1}(\eta) Re^{5/2} + \alpha_{2}(\eta) Re^{3}, \\ \hat{\epsilon}_{\vartheta} = \beta_{1}(\eta) Re^{1/2} + \beta_{2}(\eta) Re, \end{array}\right\}\end{equation}

where the bulk and BL contributions for viscous and thermal energy dissipation rates from (4.18a,b) and (4.19a,b) are used. Here, we denote the curve fits by the tilded variables and the simulation data by the plain variables. The coefficients $\alpha _{1}, \alpha _{2}, \beta _{1}$ and $\beta _{2}$ depend only on geometry (Grossmann & Lohse 2000), which means they only depend on $\eta$ in a spherical shell configuration. To fit (4.20) to our data, we proceed as

(4.21) \begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle \log_{10} \epsilon_{u} = \log_{10}\left[\alpha_{1}(\eta) \exp \left(\dfrac{5}{2} \dfrac{\log_{10}Re}{\log_{10}e}\right) + \alpha_{2}(\eta) \exp\left(3 \dfrac{\log_{10}Re}{\log_{10}e} \right)\right], \\ \displaystyle \log_{10} \epsilon_{\vartheta} = \log_{10} \left[\beta_{1}(\eta) \exp\left(\dfrac{1}{2} \dfrac{\log_{10}Re}{\log_{10}e} \right) + \beta_{2}(\eta) \exp\left(\dfrac{\log_{10}Re}{\log_{10}e} \right) \right]. \end{array}\right\}\end{equation}

For each radius ratio, the coefficients $\alpha _1, \alpha _2$, $\beta _1$ and $\beta _2$ are determined by applying the least squares method to $\log _{10}(\epsilon _{u})$–$\log _{10}(Re)$ and $\log _{10}(\epsilon _{\vartheta })$–$\log _{10}(Re)$ data obtained from our DNS database. The least squares estimate of these coefficients at different $\eta$ are listed in table 2.

Table 2. Least squares estimate of coefficients $\alpha _1, \alpha _2$, $\beta _1$ and $\beta _2$ in (4.20).

In figures 22 and 23, the curve fits represented by (4.21) are shown. Figure 22(a) shows $\epsilon _{u}$ as a function of $Re$. Figure 22(b) shows $\epsilon _{u}$ normalized by the curve fit $\hat {\epsilon _{u}}$. The curve fits show good agreement with the DNS data in both the log-log and the semi-log plots. Figure 23(a) shows the thermal energy dissipation rate $\epsilon _{\vartheta }$ and the corresponding fits $\hat {\epsilon _{\vartheta }}$ at all $\eta$, while figure 23(b) shows the $\epsilon _{\vartheta }$ data normalized by $\hat {\epsilon _{\vartheta }}$. The curve fits in this case also agree well with the DNS data.

Figure 22. (a) Volume-averaged kinetic energy dissipation rate $\epsilon _{u}$ as a function of $Re$. Symbols represent DNS data at different $\eta$. Dashed light and dark red lines represent the predictions for $\eta = 0.2$ and $\eta = 0.8$, respectively, given by (4.20). (b) Ratio between volume-averaged kinetic energy dissipation rates obtained from the DNS data ($\epsilon _{u}$) and the predictions ($\hat {\epsilon _{u}}$) given by (4.20), as a function of $Re$.

Figure 23. (a) Volume-averaged thermal energy dissipation rate $\epsilon _{\vartheta }$ as a function of $Re$. Symbols represent DNS data, while the dashed lines with the corresponding colours denote predictions given by (4.20). (b) Ratio between volume-averaged thermal energy dissipation rates obtained from the DNS data ($\epsilon _{\vartheta }$) and the predictions ($\hat {\epsilon _{\vartheta }}$) given by (4.20), as a function of $Re$.

Using (2.3) in (4.20), we can write

(4.22) \begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle \dfrac{3}{1+\eta+\eta^2} \dfrac{Ra}{Pr^2} (Nu-1) = \alpha_{1}(\eta) Re^{5/2} + \alpha_{2}(\eta) Re^{3}, \\ \displaystyle \dfrac{3 \eta}{1+\eta+\eta^2} Nu = \beta_{1}(\eta) Re^{1/2} + \beta_{2}(\eta) Re, \end{array}\right\}\end{equation}

with the coefficients $\alpha _1, \alpha _2, \beta _1$ and $\beta _2$ at each $\eta$ estimated as introduced above. At a fixed $Pr$ and $\eta$, from (4.22), we can numerically estimate $Nu$ and $Re$ as functions of $Ra$. Furthermore, to quantify the trend of the exponents, we can also define the local effective exponents as

(4.23a,b)\begin{equation} \alpha_{eff} = \frac{\partial \ln Nu }{\partial \ln Ra}; \quad \beta_{eff} = \frac{\partial \ln Re }{\partial \ln Ra}. \end{equation}

The radius ratio dependence of $Nu \sim f(Ra)$ is shown in figure 24. In figure 24(a), a good agreement between the DNS data and the GL theory predictions is observed. It can be seen from figure 24(a) that in the considered $Ra$ range, $Nu$ is small for the smaller radius ratios. The data corresponding to $\eta =0.2$ and $0.3$ clearly show this trend. For the rest of $\eta$, the values of $Nu$ are very close to each other. A look at the GL theory predictions reveals that $Nu$ for $\eta =0.3$ exceeds those for other $\eta$ in the range $7 \times 10^{8} < Ra < 3 \times 10^{9}$, and it reaches the largest value at $Ra=10^{10}$. We can also see that the GL theory predictions for $\eta =0.2$ will exceed all the other radius ratio cases (except $\eta =0.3$) in the range $4 \times 10^9 < Ra < 9 \times 10^{9}$. This trend can be explained clearly from the behaviour of the local effective exponents of $Nu(Ra)$. As shown in figure 24(b), at a fixed $Ra$, the smaller radius ratio ($\eta = 0.2, 0.3$) cases have a larger local effective exponent. A larger local effective exponent indicates a faster increase in $Nu$ with respect to $Ra$. Based on the magnitude of the local effective exponents, $Nu$ at $\eta =0.2$ is expected to finally exceed the $\eta =0.3$ curve beyond a certain $Ra$. However, the local effective exponents for the $\eta = 0.5, 0.6$ and $0.8$ cases are indistinguishable, which is consistent with the observation that the response of the system asymptotes to the planar case as $\eta \rightarrow 1$.

Figure 24. (a) $NuRa^{-1/3}$ as a function of $Ra$. Symbols represent DNS data; solid lines with corresponding colours are numerical results calculated from (4.22). (b) Local exponents for $Nu(Ra)$ using GL theory predictions given by (4.23a,b). The data are plotted at all radius ratios.

In figure 24(a), the data show that for all $\eta$, $Nu$ approaches $Nu \sim Ra^{1/3}$ at $Ra \approx 10^{9}$. By using the GL theory, we can predict the transitional $Ra$, where $Nu(Ra)$ attains the local scaling $Nu \sim Ra^{1/3}$. In figure 24(b), by looking at the local effective exponents predicted by the GL theory, we find that the $\eta =0.2$ curve reaches the $1/3$ local exponent at approximately $Ra = 1.4 \times 10^{8}$, while the $\eta =0.3$ and $0.4$ curves reach this local exponent at approximately $Ra = 4.1 \times 10^{8}$ and $8.6 \times 10^{8}$, respectively. For $\eta =0.5, 0.6$ and $0.8$, the transitional $Ra$ are so close that they are indistinguishable in the plot. The transitional $Ra$ for $\eta =0.5, 0.6$ and $0.8$ systems is roughly at $Ra = 10^{9}$. The behaviour of the transitional $Ra$ indicates that the transition to an enhanced heat transport regime would occur at a much lower $Ra$ in spherical shells than in Cartesian or cylindrical cells. For example, experiments by Roche et al. (2010) in cylindrical cells show an enhanced scaling of $Nu \sim Ra^{0.33}$ for $Ra = 5 \times 10^{11}$, which is much higher than the transitional $Ra$ observed in our system, although the thermal boundary conditions in their experiments are different from our work. In essence, GL theory predicts that the smaller radius ratio system may reach the $Nu \sim Ra^{1/3}$ scaling at a relatively lower $Ra$. This prediction may also imply that the smaller radius ratio system may reach the ultimate regime at a lower $Ra$. Furthermore, figure 24(b) shows that the differences in the local effective exponent at different $\eta$ are smaller at larger $\eta$. This is again in line with our expectation that the system will asymptotically approach the planar case with an infinite aspect ratio as $\eta \rightarrow 1$.

The radius ratio dependence of the compensated $Re \sim f(Ra)$ is shown in figure 25. The predictions of GL theory are in a good agreement with the DNS data. It is interesting to note that the behaviour of the $Re(Ra)$ relation is quite different for the smaller radius ratio ($\eta =0.2, 0.3$ and $\eta =0.4$) cases compared with the larger radius ratio ($\eta =0.5, 0.6$ and $\eta =0.8$) cases. The $Re(Ra)$ curves for the $\eta =0.5, 0.6$ and $0.8$ cases follow the same trend and are almost parallel to each other. However, the $Re(Ra)$ curves for $\eta =0.2, 0.3$ and $0.4$ show different trends and are observed to have some cross-overs with the $Re(Ra)$ curves at larger $\eta$. For example, at $Ra = 3 \times 10^{5}$, the $\eta =0.2$ case has the lowest $Re$, which increases with increasing $Ra$, faster than the rest of the cases. As mentioned above, some cross-overs with the rest of the $Re$ curves are also observed. The predictions from the GL theory in figure 25(a) indicate that at $Ra \approx 10^{10}$, $Re$ for the $\eta =0.2$ case will increase to be the largest. Different trends for different radius ratio cases can be explained by looking at the local effective exponents for $Re$ (see (4.23a,b)) in figure 25(b). At a fixed $Ra$, the local effective exponent is the largest at $\eta =0.2$, followed by the $\eta = 0.3$ and $0.4$ cases. For the $\eta \geq 0.5$ cases, the local effective exponents are indistinguishable from each other in figure 25(b). Similar to the larger $\alpha _{eff}$ for $Nu$, a larger local effective exponent $\beta _{eff}$ for $Re$ indicates a faster increase in $Re$ with respect to $Ra$, which is consistent with the steep rise of $Re$ for the $\eta =0.2$ cases as shown in figure 25(a). Furthermore, GL theory predictions of $\beta _{eff}$ indicate that the smaller radius ratio system may reach the ultimate scaling for $Re(Ra)$ at a relatively lower $Ra$, similar to what we observe for $Nu(Ra)$. Moreover, figure 25(b) shows that the differences in $\beta _{eff}$ at different $\eta$ are smaller for larger $\eta$, in line with the fact that the system approaches the planar case as $\eta \rightarrow 1$.

Figure 25. (a) $ReRa^{-1/2}$ as a function of $Ra$. Symbols represent DNS data; solid lines with corresponding colours are the GL theory predictions given by (4.22). (b) Local exponents for $Re(Ra)$ using GL theory predictions given by (4.23a,b). The data are plotted at all radius ratios.

In the preceding analysis, it is natural to question how the scaling behaviour of $Nu$ and $Re$ relates to the degree of supercriticality $Ra/Ra_c$ instead of $Ra$, where $Ra_{c}$ is the critical $Ra$ for the onset of convection. Plotting the scaling relations $Nu\sim Nu(Ra/Ra_c)$ and $Re\sim Re(Ra/Ra_c)$, the corresponding local exponents for different $\eta$ revealed that the use of $Ra/Ra_c$ as an independent variable does not change any of our conclusions and adds nothing new or interesting to the discussion.

5. Asymmetry of the velocity field

In contrast to the planar RBC, time and horizontally averaged radial profiles of temperature and velocity in spherical shells are asymmetric. For example, the temperature drop inside the inner shell thermal BL and the Reynolds number near the inner boundary are larger than those near the outer boundary. These radial asymmetries are directly inherited from the geometrical asymmetry present in the system. Therefore, quantifying the asymmetry of the radial profiles as a function of the radius ratio is important. For the mean temperature profile $\vartheta (r)$, the properties in which we are interested are the thermal BL thicknesses $\lambda _{\vartheta _i}$ and $\lambda _{\vartheta _o}$, as well as the mean temperature drop within the BLs $\Delta \vartheta _i$ and $\Delta \vartheta _o$. The subscripts ‘$i$’ and ‘$o$’ refer to the inner and the outer BLs, respectively. The ratios $\lambda _{\vartheta _i}/\lambda _{\vartheta _o}$ and $\Delta \vartheta _i / \Delta \vartheta _o$ provide good measures to quantify the asymmetry of the mean temperature profile. The theoretical estimates from Gastine et al. (2015), stated in (4.6a–c), agree well with our DNS data; see figure 13.

To quantify the asymmetry of the Reynolds number $Re(r)$, we divide the domain into an inner shell and an outer shell at the middle plane. The Reynolds numbers corresponding to the inner shell ($Re_i$) and the outer shell ($Re_o$) are defined as

(5.1a,b)\begin{equation} Re_{i} = \frac{1}{V_{i}} \int_{r_{i}}^{r_{mid}}Re(r)4{\rm \pi} r^{2} \mathrm{d}r, \quad Re_{o} = \frac{1}{V_{o}} \int_{r_{mid}}^{r_{o}}Re(r)4{\rm \pi} r^{2} \mathrm{d}r. \end{equation}

Here, $r_{mid}=(r_{i}+r_{o})/2$ is the radius at mid-depth, while $V_{i}, V_{o}$ are the inner and the outer shell volumes, respectively. From figure 12(b), it can be seen that the time and horizontally averaged $Re$ is higher near the inner boundary in comparison with that of the outer boundary. A smaller radius ratio results in a larger difference, pointing to a clear $\eta$-dependence. The ratio $Re_{i}/Re_{o} \sim f(\eta )$ is used to quantify the asymmetry of the velocity profile.

To define the local Rayleigh numbers, $Ra_i$ and $Ra_o$ for the inner and outer sub-shells, respectively, we assume symmetric temperature drops: $2\Delta \vartheta _{i}$ for the inner sub-shell and $2\Delta \vartheta _{o}$ for the outer sub-shell. With this definition, we could write (Tisserand et al. 2011; Wei et al. 2014; Zhu, Verzicco & Lohse 2017)

(5.2) \begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle Ra_{i} = \dfrac{\alpha g_{o} (2\Delta \vartheta_{i}) d^{3}}{\nu \kappa}=\dfrac{2\Delta \vartheta_{i}}{\Delta \vartheta} \dfrac{\alpha g_{o} \Delta \vartheta d^{3}}{\nu \kappa} = \chi Ra, \\ \displaystyle Ra_{o} = \dfrac{\alpha g_{o} (2\Delta \vartheta_{o}) d^{3}}{\nu \kappa}=2 \left(1-\dfrac{\Delta \vartheta_{i}}{\Delta \vartheta} \right) \dfrac{\alpha g_{o} \Delta \vartheta d^{3}}{\nu \kappa} = (2-\chi) Ra, \end{array}\right\} \end{equation}

as the Rayleigh numbers for the inner and outer shells, respectively. Here, $\chi = 2\Delta \vartheta _{i}/\Delta \vartheta$ is the asymmetry factor defined in (4.7).

The general form of $Re(Ra)$ and $Re_{i}(Ra_{i})$ relations can be written as

(5.3a,b)\begin{equation} Re = h(\eta) Ra^{\beta(\eta,Ra)}, \quad Re_{i} = h_{i}(\eta) Ra_{i}^{\beta_{i}(\eta,Ra_i)}. \end{equation}

From figure 12(b), it is clear that the flow in the inner shell is more turbulent than in the outer shell, and thus the global Reynolds number $Re$ is dominated by the inner shell contributions. Figure 26 shows $Re$, $Re_i$ and $Re_o$ as functions of $Ra$, $Ra_i$ and $Ra_o$, respectively. Since the global Reynolds number $Re$ is dominated by the contributions from the inner shell, the $Re(Ra)$ and $Re_i(Ra_i)$ curves in figure 26 are indistinguishable for all $\eta$. Hence, we can conveniently write $h(\eta )=h_{i}(\eta )$ and $\beta (\eta,Ra)=\beta _{i}(\eta,Ra_i)$. Here, we use $\beta (\eta,Ra) \approx 1/2$ as the estimation of the $Re(Ra)$ scaling. This scaling works well in our simulations, as shown in figure 25, and has also been observed in many other experiments and simulations (Tilgner 1996; Brown, Funfschilling & Ahlers 2007; Chillà & Schumacher 2012). From (5.2) and (5.3a,b), we can write

(5.4)\begin{equation} \frac{Re_{i}}{Re} = \left(\frac{Ra_{i}}{Ra} \right)^{\beta} \approx \left(\frac{Ra_{i}}{Ra} \right)^{{1}/{2}} = \chi^{{1}/{2}}.\end{equation}

From the definitions of $Re_i$ and $Re_o$, it follows that

(5.5)\begin{equation} Re {\cdot} V = Re_{i} {\cdot} V_{i} + Re_{o} {\cdot} V_{o},\end{equation}

where $V=(4/3) {\rm \pi}(r_{o}^3-r_{i}^3)$ is the total shell volume, and $V_i=(4/3) {\rm \pi}(r_{m}^3-r_{i}^3)$ and $V_{o}=(4/3) {\rm \pi}(r_{o}^3-r_{m}^3)$ are the inner and outer shell volumes, respectively. Combining (4.7), (5.4) and (5.5), we arrive at the following relation:

(5.6)\begin{equation} \frac{Re_i}{Re_o} = \frac{\left( \dfrac{2}{1+\eta^{{4}/{3}}}\right)^{{1}/{2}}(\eta^2 +4\eta+7)}{8(\eta^2+\eta+1)-\left(\dfrac{2}{1+\eta^{{4}/{3}}}\right)^{{1}/{2}} (7\eta^2+4\eta+1)}. \end{equation}

Equation (5.6) represents the ratio of the inner and outer shell Reynolds numbers as a function of the radius ratio, which can be regarded as a measure to quantify the asymmetry of the velocity field in our system.

Figure 26. Global, inner shell and outer shell Reynolds numbers $(Re, Re_i, Re_o)$ as functions of the global, inner shell and outer shell Rayleigh numbers $(Ra, Ra_i, Ra_o)$, respectively. (a) $\eta =0.2$, (b) $\eta =0.3$, (c) $\eta =0.4$, (d) $\eta =0.5$, (e) $\eta =0.6$ and (f) $\eta =0.8$.

Equation (5.6) is validated against our DNS data in figure 27. It can be observed that the predictions by (5.6) agree reasonably well with the data. We can also see from the data that $Re_{i}/Re_{o}$ exhibits a dependence on $Ra$, albeit weak, especially for the small radius ratio cases. As $\eta \rightarrow 1$, $Re_i/Re_o \rightarrow 1$, indicating the asymmetry of the velocity field will become smaller at a higher $\eta$, as expected.

Figure 27. Ratio of the inner and outer shell Reynolds numbers, $Re_{i}/Re_{o}$, as a function of $\eta$. Symbols represent DNS data. The dashed horizontal line corresponds to $Re_{i}/Re_{o} = 1$.

6. Asymmetry of energy dissipation rates

As a result of the asymmetry in the thermal and viscous BLs, the energy dissipation rates can be divided into the partial contributions from the two asymmetric BLs and an interior bulk region,

(6.1) \begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle \epsilon_{u} = \epsilon_{u}^{bulk}+\epsilon_{u}^{BL,inner}+\epsilon_{u}^{BL,outer}, \\ \displaystyle \epsilon_{\vartheta} = \epsilon_{\vartheta}^{bulk}+\epsilon_{\vartheta}^{BL,inner}+\epsilon_{\vartheta}^{BL,outer}, \end{array}\right\} \end{equation}

where we use subscripts bulk, BL to represent the bulk and BL contributions, respectively. The subscripts inner and outer denote contributions from the inner BL and the outer BL, respectively. The inner and the outer BL kinetic energy and thermal energy dissipation rates respectively are

(6.2) \begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle \epsilon_{u}^{BL,inner} = \dfrac{4 {\rm \pi}}{V} \int_{r_i}^{r_i+\lambda_{u}^{i}} \overline{\biggl\langle \left(\dfrac{\partial u_i}{\partial x_j} \right)^2 \biggr\rangle}_{s} r^{2}\,\mathrm{d}r, \\ \displaystyle \epsilon_{u}^{BL,outer} = \dfrac{4 {\rm \pi}}{V} \int_{r_o-\lambda_{u}^{o}}^{r_o} \overline{\biggl\langle \left(\dfrac{\partial u_i}{\partial x_j} \right)^2 \biggr\rangle}_{s} r^{2}\,\mathrm{d}r, \end{array}\right\} \end{equation}

and

(6.3) \begin{equation} \left. \begin{array}{@{}c@{}} \displaystyle \epsilon_{\vartheta}^{BL,inner} = \dfrac{4 {\rm \pi}}{V} \int_{r_i}^{r_i+\lambda_{\vartheta}^{i}} \overline{\left\langle \left(\Delta T \right)^2 \right\rangle}_{s} r^{2}\,\mathrm{d}r, \\ \displaystyle \epsilon_{\vartheta}^{BL,outer} = \dfrac{4 {\rm \pi}}{V} \int_{r_o-\lambda_{\vartheta}^{o}}^{r_o} \overline{\left\langle \left(\Delta T \right)^2 \right\rangle}_{s} r^{2}\,\mathrm{d}r. \end{array}\right\} \end{equation}

The relative contributions of the thermal and kinetic energy dissipation rates from the bulk and BLs vary with $Ra$ and $\eta$. We explore this in the following subsections.

6.1. Asymmetry of thermal energy dissipation rate

Figure 28 shows the relative contributions of the bulk and BL regions for the thermal energy dissipation rate as a function of $Ra$, at different $\eta$. The total BL contribution is always larger than the bulk contribution for the considered $Ra$ range. With increasing $Ra$, the bulk contribution will also increase. In figure 28, it is seen that for the thermal energy dissipation rate, the inner BL contribution is always larger than the outer BL contribution. The bulk and the total BL contributions as functions of both $Ra$ and $\eta$ are shown in figure 29. The relative contributions of the bulk and the sum of the two BLs are almost independent of $\eta$.

Figure 28. Contributions of the bulk and BL regions to the thermal energy dissipation rate as functions of $Ra$. Curves are plotted for all $\eta$. Solid circles represent the bulk contribution, while the solid squares represent the total BL contributions. Open downward triangles and open upward triangles denote the inner and the outer BL contributions, respectively. (a) $\eta =0.2$, (b) $\eta =0.3$, (c) $\eta =0.4$, (d) $\eta =0.5$, (e) $\eta =0.6$ and (f) $\eta =0.8$.

Figure 29. Relative contributions of the bulk and the BLs (inner $+$ outer) to the thermal energy dissipation rate as a function of $Ra$ at different $\eta$.

In figure 23, for $\eta \geq 0.3$, we observe cross-over points between the bulk and the outer BL contribution curves. From the GL theory, it is known that the energy dissipation rates will gradually transition from BL-dominant to bulk-dominant, as $Ra$ increases. At small $Ra$, the flow is quasi-laminar with a thicker BL and a less turbulent bulk, implying a larger fraction of the BL dissipation and a relatively smaller fraction of the bulk dissipation. As $Ra$ increases, the BL shrinks in size and the bulk becomes more turbulent, leading to an increase in the fraction of the bulk dissipation. Since the outer BL contribution is smaller than the inner BL contribution for the thermal energy dissipation rate, the bulk contribution will first exceed the outer BL contribution, followed by the inner BL contribution, and will finally exceed the total BL contribution at an even larger $Ra$. This explains the cross-overs observed in the thermal dissipation rate data. As a result of the limited $Ra$ range considered in this simulation campaign, the cross-over point between the inner BL dissipation and the bulk dissipation has not been observed yet. However, in figure 28, we can clearly observe a trend that the bulk contribution would exceed the inner BL contribution at an even higher $Ra$.

In figure 28, it is also observed that $Ra$ at the cross-over point is larger for higher $\eta$ cases. For example, for $\eta =0.3$, the bulk contribution exceeds the outer BL contribution at approximately $Ra=3 \times 10^{3}$, while for $\eta =0.8$, the cross-over point is at $Ra > 5 \times 10^{8}$. Since both $\epsilon _{\vartheta }^{BL}$ and $\epsilon _{\vartheta }^{bulk}$ stay almost the same with increasing $\eta$, and $\epsilon _{\vartheta }^{BL,outer}$ increases, the cross-over point shifts towards the higher $Ra$ end of the plot.

To quantify the asymmetry of the thermal energy dissipation rate in the inner and the outer BLs, we estimate

(6.4)\begin{equation} \epsilon_{\vartheta}^{BL,inner} = \frac{4 {\rm \pi}}{V} \int_{r_i}^{r_i+\lambda_{\vartheta}^{i}} \overline{\langle (\Delta T )^2 \rangle}_{s} r^{2}\,\mathrm{d}r \approx \biggl(\frac{\Delta \vartheta_{i}}{\lambda_{\vartheta}^{i}}\biggr)^{2} \frac{4{\rm \pi} r_{i}^{2}}{V} \lambda_{\vartheta}^{i}\end{equation}

and

(6.5)\begin{equation} \epsilon_{\vartheta}^{BL,outer} = \frac{4 {\rm \pi}}{V} \int_{r_o-\lambda_{\vartheta}^{o}}^{r_o} \overline{\langle (\Delta T )^2 \rangle}_{s} r^{2}\,\mathrm{d}r \approx \left(\frac{\Delta \vartheta_{o}}{\lambda_{\vartheta}^{o}} \right)^{2} \frac{4{\rm \pi} r_{o}^{2}}{V} \lambda_{\vartheta}^{o}.\end{equation}

Combining the above equations with (4.6a–c) leads to

(6.6)\begin{equation} \frac{\epsilon_{\vartheta}^{BL,inner}}{\epsilon_{\vartheta}^{BL,outer}} \sim \left(\frac{\Delta \vartheta_{i}}{\Delta \vartheta_{o}} \right)^{2} \frac{\lambda_{\vartheta}^{o}}{\lambda_{\vartheta}^{i}} \eta^{2} = \eta^{{-}4/3}.\end{equation}

In figure 30, it can be seen that (6.6) agrees well with the DNS data.

Figure 30. Ratio of the inner and outer BL contributions of thermal energy dissipation rate as a function of $\eta$. The dashed line denotes $\epsilon _{\vartheta }^{BL,inner}/\epsilon _{\vartheta }^{BL,outer} = 1$. Open symbols represent the DNS data.

6.2. Asymmetry of kinetic energy dissipation rate

Relative contributions of the kinetic energy dissipation rate $\epsilon _{u}^{bulk}$, $\epsilon _{u}^{BL}$, $\epsilon _{u}^{BL,inner}$ and $\epsilon _{u}^{BL,outer}$ are shown in figure 31. The kinetic energy dissipation rate is bulk dominant since the bulk contribution is always larger than the total BL contribution. For the $\eta \geq 0.3$ cases, $\epsilon _{u}^{bulk}$ keeps increasing with increasing $Ra$. However, at $\eta = 0.2$, $\epsilon _{u}^{bulk}$ will first increase and then decrease at a higher $Ra$. Furthermore, figure 32 shows that there is no $\eta$ dependence on the fractions of the bulk and the total BL contributions for kinetic energy dissipation rate.

Figure 31. Bulk and BL (inner $+$ outer) contributions of the kinetic energy dissipation rate as functions of $Ra$ for different radius ratios. Solid circles represent the bulk contribution and the solid squares represent the total BL contributions. Open downward and upward triangles denote inner and outer BL contributions, respectively. Dashed lines with corresponding colours are added to improve the readability of the plot. (a) $\eta =0.2$, (b) $\eta =0.3$, (c) $\eta =0.4$, (d) $\eta =0.2$, (e) $\eta =0.6$ and (f) $\eta =0.8$.

Figure 32. Relative contributions of the bulk and the BLs (inner $+$ outer) to the kinetic energy dissipation rate as the function of $Ra$ and different $\eta$.

In contrast to the thermal energy dissipation rate, the outer BL contribution is always larger than the inner BL contribution for the kinetic energy dissipation rate. Following a similar procedure as in § 6.1, the inner and the outer BL kinetic energy dissipation rates are estimated as

(6.7) \begin{equation} \epsilon_{u}^{BL,inner} = \frac{4 {\rm \pi}}{V} \int_{r_i}^{r_i+\lambda_{u}^{i}} \overline{\biggl\langle \left(\frac{\partial u_i}{\partial x_j} \right)^2 \biggr\rangle}_{s} r^{2}\,\mathrm{d}r \approx \left(\frac{Re_{i}}{\lambda_{u}^{i}} \right)^{2} \frac{4{\rm \pi} r_{i}^{2}}{V} \lambda_{u}^{i}\end{equation}

and

(6.8) \begin{equation} \epsilon_{u}^{BL,outer} = \frac{4 {\rm \pi}}{V} \int_{r_o-\lambda_{u}^{o}}^{r_o} \overline{\biggl\langle \left(\frac{\partial u_i}{\partial x_j} \right)^2 \biggr\rangle}_{s} r^{2}\,\mathrm{d}r \approx \left(\frac{Re_{o}}{\lambda_{u}^{o}} \right)^{2} \frac{4{\rm \pi} r_{o}^{2}}{V} \lambda_{u}^{o},\end{equation}

respectively. Using (5.6) in combination with (6.7) and (6.8), we get

(6.9) \begin{align} \frac{\epsilon_{u}^{BL,inner}}{\epsilon_{u}^{BL,outer}} \sim \left(\frac{Re_{i}}{Re_{o}} \right)^{2} \frac{\lambda_{u}^{o}}{\lambda_{u}^{i}} \eta^{2} =\frac{\eta^{4/3}\left(\dfrac{2}{1+\eta^{{4}/{3}}}\right)^{{1}/{2}}(\eta^2 +4\eta+7)}{8(\eta^2+\eta+1)-\left(\dfrac{2}{1+\eta^{{4}/{3}}}\right)^{{1}/{2}} (7\eta^2+4\eta+1)} < 1.\end{align}

Equation (6.9) ensures $\epsilon _{u}^{BL,outer} > \epsilon _{u}^{BL,inner}$.

Figure 33 shows $\epsilon _{u}^{BL,inner} / \epsilon _{u}^{BL,outer}$ as a function of $\eta$. The theoretical estimation of (6.9) predicts the trend reasonably well. However, there are significant deviations from the data, particularly at smaller $\eta$. The smaller the radius ratio, the larger the scatter in the $\epsilon _{u}^{BL,inner} / \epsilon _{u}^{BL,outer}$ data with respect to $Ra$. The deviations come from the $Ra$ dependence of $\epsilon _{u}^{BL,inner} / \epsilon _{u}^{BL,outer}$ especially in the small $\eta$ cases, since both $Re_i/Re_o$ and $\lambda _{u}^{i}/\lambda _{u}^{o}$ have $Ra$ dependence as illustrated in figures 13 and 27. Further investigations are still needed for demonstrating the $Ra$-dependence of $\epsilon _{u}^{BL,inner} / \epsilon _{u}^{BL,outer}$.

Figure 33. Ratio of the inner and outer BL contributions of kinetic energy dissipation rate as a function of $\eta$. Open symbols represent the DNS data.

It must be stressed here that owing to the dominance of $\epsilon _{\vartheta }^{BL} \textrm { and } \epsilon _{u}^{bulk}$, our simulations lie in regime II of the GL theory, in which the boundary layer and the bulk contributions dominate the thermal dissipation rate and the kinetic energy dissipation rate, respectively (Grossmann & Lohse 2000).

7. Conclusion

Three-dimensional DNS for spherical shell thermal convection with $0.2 \leq \eta \leq 0.8$, $Ra \leq 5 \times 10^{8}$ and $Pr=1$ were conducted in this study. We used the centrally condensed gravity profile ($g \sim (r_0/r)^{2}$) in our simulations so that the exact relations for the dissipation rates exist. A qualitative look at the flow structures reveals the existence of large-scale structures. At a fixed radius ratio, the number of large-scale structures decreases with increasing $Ra$, while at a fixed $Ra$, the number of large-scale structures increases with $\eta$. The integral length scale at the horizontal mid-depth was found to be saturated at approximately 1.5 for $\eta \geq 0.2$.

From the BL scaling analysis, it was found that the inner thermal BL thickness scales as $\lambda _{\vartheta }^{i} \sim \eta ^{1/2} Ra^{-0.295}$ in our simulations. Using the definition of $Nu$ as well as the asymmetry of the thermal BLs of Gastine et al. (2015), we derived the radius ratio dependence of the Nusselt number as $Nu \sim ({\eta ^{1/2}}/({1+\eta ^{4/3}})) Ra^{\gamma }$, where $\gamma \approx 0.295$ for the fit with one single exponent in the whole data range considered. However, it is well known that the local effective scaling exponent varies with respect to $Ra$ (Grossmann & Lohse 2000, 2001; Xu, Bajaj & Ahlers 2000; Ahlers et al. 2009; Zhu et al. 2018). To quantify how the local exponent changes with $Ra$, we employ the GL theory to investigate the non-power law effects of $Nu(Ra), Re(Ra)$ scalings. Larger local effective exponents $\alpha _{eff}$ and $\beta _{eff}$ for $Nu \sim Ra^{\alpha _{eff}}$ and $Re \sim Ra^{\beta _{eff}}$ are observed in lower $\eta$. From the prediction of the GL theory, the transition point to an enhanced heat transport regime ($\alpha _{eff} \geq 1/3$) is found to occur at a lower $Ra$ for a smaller $\eta$.

To assess the asymmetry in the velocity field, we introduce $Re_i$ and $Re_o$ for the inner and outer shells, respectively. The dominance of $Re_i$ in the global Reynolds number $Re$ is evident from the similarity in the scaling exponents and prefactors of $Re(Ra)$ and $Re_i(Ra_i)$ relations. Leveraging this similarity, we establish an empirical relation $Re_{i}/Re_{o}(\eta )$ as a quantification of the velocity field's asymmetry.

The investigation extends to the asymmetry of thermal and kinetic energy dissipation rates. Regarding thermal energy dissipation, it is observed that the inner BL contribution surpasses the outer BL contribution. For higher radius ratio cases, $Ra$ at the cross-over point between $\epsilon _{\vartheta }^{bulk}$ and $\epsilon _{\vartheta }^{BL,outer}$ increases, which is explained by the system becoming less asymmetric at higher $\eta$. An analytical relation $\epsilon _{\vartheta }^{BL,inner}/ \epsilon _{\vartheta }^{BL,outer}(\eta )$ is derived and validated against DNS data. For kinetic energy dissipation, the inner BL contribution is smaller than the outer BL contribution, with a provided relation for $\epsilon _{u}^{BL,inner}/\epsilon _{u}^{BL,outer}(\eta )$. Although the provided relation reasonably predicts the observed trend, its accuracy diminishes at smaller $\eta$. The deviations are attributed to the strong $Ra$-dependence. Additionally, we find that, for both thermal and kinetic energy dissipation rates, the fractions of bulk and total BL contributions remain independent of $\eta$.

Further investigations are needed to answer some open questions. For example, further simulations at $Ra > 10^{9}$ should be performed to really show the transition to the $Nu \sim Ra^{1/3}$ regime and beyond. Larger scaling exponents for $Nu(Ra), Re(Ra)$ at smaller $\eta$ were obtained from this study. A physical explanation for this behaviour will be sought in the future research. It still needs to be verified whether or not the smaller $\eta$ cases will transition to the ultimate regime at smaller $Ra$, and is expected to be the subject of the future editions of this work. Moreover, whether the scaling behaviour for $Nu(Ra), Re(Ra)$, as well as the the asymmetry of the velocity field $Re_{i}/Re_{o}(\eta )$ proposed in this work, will hold in other gravity profiles and different $Pr$ remains to be investigated in the future.