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 Reference Hartmann, Moy and Fu2001), mantle convection within the Earth (Davies & Richards Reference Davies and Richards1992), thermal convection in Earth's outer core (Wicht & Sanchez Reference Wicht and Sanchez2019), deep convection in the molecular envelopes of gas giants (Aurnou et al. Reference Aurnou, Heimpel, Allen, King and Wicht2008) and thermal convection in the convective zone of the Sun (Hanasoge, Gizon & Sreenivasan Reference Hanasoge, Gizon and Sreenivasan2016). 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 Reference Phillips and Lambeck1980). 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 Reference Gastine, Wicht and Aurnou2015), 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 (Reference Ahlers, Grossmann and Lohse2009), Chillà & Schumacher (Reference Chillà and Schumacher2012), Plumley & Julien (Reference Plumley and Julien2019) and Xia et al. (Reference Xia, Huang, Xie and Zhang2023).
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. Reference Egbers, Beyer, Bonhage, Hollerbach and Beltrame2003; Futterer et al. Reference Futterer, Egbers, Dahley, Koch and Jehring2010, Reference Futterer, Krebs, Plesa, Zaussinger, Hollerbach, Breuer and Egbers2013; Zaussinger et al. Reference Zaussinger, Haun, Neben, Seelig, Travnikov, Egbers, Yoshikawa and Mutabazi2018, Reference Zaussinger, Haun, Szabo, Peter, Travnikov, Al Kawwas and Egbers2020) 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 Reference Zebib, Schubert and Straus1980; Zebib et al. Reference Zebib, Schubert, Dein and Paliwal1983; Bercovici et al. Reference Bercovici, Schubert, Glatzmaier and Zebib1989; Bercovici, Schubert & Glatzmaier Reference Bercovici, Schubert and Glatzmaier1992; Tilgner Reference Tilgner1996; Tilgner & Busse Reference Tilgner and Busse1997; Choblet & Parmentier Reference Choblet and Parmentier2009; Choblet Reference Choblet2012; Gastine et al. Reference Gastine, Wicht and Aurnou2015) have focused on non-rotating spherical shell thermal convection; some of which consider an infinite $Pr$ (Bercovici et al. Reference Bercovici, Schubert and Glatzmaier1992), 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 (Reference Tilgner1996) 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. (Reference Gastine, Wicht and Aurnou2015) 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. (Reference Gastine, Wicht and Aurnou2015) 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. Reference Genova2019). Similarly, Earth's outer core and inner core possess a radius ratio of approximately 0.35 (Ahrens Reference Ahrens1995). 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. Reference Guillot, Stevenson, Hubbard and Saumon2004), 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 Reference Christensen and Wicht2008; Vazan et al. Reference Vazan, Helled, Podolak and Kovetz2016). 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 Reference O'Neill2021; Stähler et al. Reference Stähler2021). Mars, however, is highly suspected not to possess a solid inner core (Stähler et al. Reference Stähler2021). For the ice giants Uranus and Neptune, their interior structures remain a subject of controversy and uncertainty (Helled, Nettelmann & Guillot Reference Helled, Nettelmann and Guillot2020). 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 Reference Gastine, Wicht and Aubert2016; Long et al. Reference Long, Mound, Davies and Tobias2020; Song et al. Reference Song, Kannan, Shishkina and Zhu2024a; Song, Shishkina & Zhu Reference Song, Shishkina and Zhu2024b) 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.
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. (Reference Gastine, Wicht and Aurnou2015). 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 Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Gastine et al. Reference Gastine, Wicht and Aurnou2015). 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.
(i) The large-scale structures are observed in spherical RBC for the first time.
(ii) An analytical relation for the radius ratio dependence of $Nu$ is derived.
(iii) The relations for local scaling exponents for $Nu(Ra)$ and $Re(Ra)$ are investigated by using GL theory.
(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
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 Reference Chandrasekhar1981). 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
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
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
Another key response parameter is the Reynolds number,
which represents the dimensionless velocity. The radial profile of the time and horizontally averaged horizontal velocity is
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. Reference Gastine, Wicht and Aurnou2015)
In our dimensionless quantities, $\epsilon _u$ and $\epsilon _{\vartheta }$ are calculated as follows:
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
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 Reference Wicht2002; Christensen & Wicht Reference Christensen and Wicht2007; Lago et al. Reference Lago, Gastine, Dannert, Rampp and Wicht2021) 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,
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. (Reference Gastine, Wicht and Aurnou2015). 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. Reference Long, Mound, Davies and Tobias2020). In this work, we have adopted the slope method (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010; Gastine et al. Reference Gastine, Wicht and Aurnou2015) 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. (Reference Gastine, Wicht and Aurnou2015). 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.
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 Reference Glatzmaier2013)
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 Reference Christensen and Wicht2007). In the above equation, the $m=0$ contribution entering in the summation is multiplied by one-half (Schwaiger, Gastine & Aubert Reference Schwaiger, Gastine and Aubert2021), 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$.
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}$.
The degree $l_{max}$ is associated with the dominant length scale $L_{dom}$ by the definition of the characteristic wavelength (Backus, Parker & Constable Reference Backus, Parker and Constable1996)
which is associated with the length scale of the large-scale structures (Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018). For comparison, we also look at the so-called integral length scale (Parodi et al. Reference Parodi, von Hardenberg, Passoni, Provenzale and Spiegel2004)
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. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) 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.
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. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) 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. Reference Wang, Verzicco, Lohse and Shishkina2020), were not observed in our 3-D cases.
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. Reference Ahlers, Grossmann and Lohse2009), 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. Reference Shishkina, Stevens, Grossmann and Lohse2010) to define the kinetic as well as the thermal BL thicknesses for consistency with the BL analyses of Prandtl (Reference Prandtl1905) and Blasius (Reference Blasius1907). The definitions of the boundary layer thicknesses are the same as those of Gastine et al. (Reference Gastine, Wicht and Aurnou2015).
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.
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
so that we have
where $\vartheta _{m}=\vartheta ((r_i+r_o)/2)$ is the time and horizontally averaged temperature at the mid-depth. Gastine et al. (Reference Gastine, Wicht and Aurnou2015) 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. (Reference Gastine, Wicht and Aurnou2015) is presented below. From (2.5), we obtain
Gastine et al. (Reference Gastine, Wicht and Aurnou2015) proposed that the average plume density in the inner and outer BLs should be the same, which yields
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. (Reference Gastine, Wicht and Aurnou2015) finally derive
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
However, there are some other models for $\Delta \vartheta _{i}, \Delta \vartheta _{o}$ and $\lambda _{\vartheta }^{i}/\lambda _{\vartheta }$ (e.g. Wu & Libchaber Reference Wu and Libchaber1991; Zhang, Childress & Libchaber Reference Zhang, Childress and Libchaber1997; Wang et al. Reference Wang, Jiang, Liu, Zhu and Sun2022). Gastine et al. (Reference Gastine, Wicht and Aurnou2015) 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,
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. Reference Shishkina, Stevens, Grossmann and Lohse2010), such that
In figure 13, (4.7) and (4.8) are successfully validated against our DNS data.
From (4.6a–c) and the $Nu$ definition in (2.5), we can directly obtain the following relations:
Thus, we can define a normalized thermal boundary layer thickness
Figure 14 shows a validation of (4.10) and it shows that the relation perfectly aligns with our DNS data.
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. Reference Ahlers, Grossmann and Lohse2009).
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
Equation (4.6a–c) gives us the $\eta$-dependence of the temperature drop in the inner BL (Gastine et al. Reference Gastine, Wicht and Aurnou2015), 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
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
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 Reference Grossmann and Lohse2000) in subsequent sections.
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:
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. (Reference Brown, Nikolaenko, Funfschilling and Ahlers2005) and Funfschilling et al. (Reference Funfschilling, Brown, Nikolaenko and Ahlers2005) in the cylindrical geometries. Figure 15 also shows that (4.14) exhibits a good agreement with our DNS data.
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 Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001). 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 Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Ahlers et al. Reference Ahlers, Grossmann and Lohse2009). 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. (Reference Gastine, Wicht and Aurnou2015) 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
Here, we use the subscripts $BL$ and $bulk$ for the BL and bulk contributions, respectively. The bulk and BL contributions are defined as follows:
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 (Reference Grossmann and Lohse2000) proposed the following scalings for $\epsilon _{u}^{bulk}$ and $\epsilon _{\vartheta }^{bulk}$ at a fixed $Pr$:
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 Reference Grossmann and Lohse2004). 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. (Reference Calzavarini, Lohse, Toschi and Tripiccione2005).
Based on the Prandtl–Blasius BL theory and the fact that the BL thickness is significantly smaller than the gap length, Grossmann & Lohse (Reference Grossmann and Lohse2000) derived the following scaling relations for $\epsilon _{u}^{BL}$ and $\epsilon _{\vartheta }^{BL}$ at $Pr=1$,
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. (Reference Gastine, Wicht and Aurnou2015), as a result of the difficulty to perfectly separate the BL from the bulk.
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,
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 Reference Grossmann and Lohse2000), which means they only depend on $\eta$ in a spherical shell configuration. To fit (4.20) to our data, we proceed as
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.
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.
Using (2.3) in (4.20), we can write
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
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$.
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. (Reference Roche, Gauthier, Kaiser and Salort2010) 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$.
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. (Reference Gastine, Wicht and Aurnou2015), 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
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. Reference Tisserand, Creyssels, Gasteuil, Pabiou, Gibert, Castaing and Chillà2011; Wei et al. Reference Wei, Chan, Ni, Zhao and Xia2014; Zhu, Verzicco & Lohse Reference Zhu, Verzicco and Lohse2017)
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
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 Reference Tilgner1996; Brown, Funfschilling & Ahlers Reference Brown, Funfschilling and Ahlers2007; Chillà & Schumacher Reference Chillà and Schumacher2012). From (5.2) and (5.3a,b), we can write
From the definitions of $Re_i$ and $Re_o$, it follows that
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:
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.
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.
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,
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
and
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$.
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
and
Combining the above equations with (4.6a–c) leads to
In figure 30, it can be seen that (6.6) agrees well with 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.
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
and
respectively. Using (5.6) in combination with (6.7) and (6.8), we get
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}$.
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 Reference Grossmann and Lohse2000).
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. (Reference Gastine, Wicht and Aurnou2015), 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 Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2001; Xu, Bajaj & Ahlers Reference Xu, Bajaj and Ahlers2000; Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018). 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.
Acknowledgements
We thank T. Gastine for providing the flow fields for $\eta =0.6$ simulations. We are also grateful to D. Lohse, O. Shishkina, J. Wicht and A. Tilgner for their illuminating discussions. All the simulations have been conducted on the HPC systems of Max Planck Computing and Data Facility (MPCDF) as well as the National High Performance Computing (NHR@ZIB and NHR-Nord@Göttingen).
Funding
We gratefully acknowledge the financial support from the Max Planck Society, the German Research Foundation through grants 521319293, 540422505 and 550262949, the Alexander von Humboldt Foundation, the International Max Planck Research School for Solar System Science at the University of Göttingen (IMPRS, Solar System School) and the Daimler Benz Foundation.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Simulation data and grid resolution.
An appropriate resolution should achieve a balance of the turbulent kinetic energy budget, which is the balance between the volume-averaged buoyancy power $\mathcal {P}_{buoy}$ and volume-averaged kinetic energy dissipation rate $\epsilon _{u}$,
We show two examples at $\eta =0.2$ and $\eta =0.8$ with $Ra=3 \times 10^{7}$ in figure 34. The time derivative of the total kinetic energy $\textrm {d}E_{kin}/\textrm {d}t$ is also shown in the figure. As it is clearly seen, the volume-averaged generation of kinetic energy from the radial buoyancy flux is balanced by the volume-averaged viscous dissipation rate for both the cases. The simulation details are shown in table 3.
Appendix B. Kinetic energy spectra calculation
In MagIC code, the kinetic energy spectra $E_{kin}(l)$ are calculated as
where the $m=0$ contribution is pre-multiplied by $1/2$. In the above equation, $W_{lm}(r)$ and $Z_{lm}(r)$ are the poloidal potential and toroidal potential of velocity in the spectral space, respectively. They are defined by
and
where $W(r, \theta, \phi )$ and $Z(r, \theta, \phi )$ are the poloidal potential and toroidal potential of velocity $\boldsymbol {\boldsymbol {u}}$, respectively. In physical space, the velocity $\boldsymbol {\boldsymbol {u}}$ can be written as
Note that both $W(r, \theta, \phi )$ and $Z(r, \theta, \phi )$ are real functions, so that we have
where the asterisk denotes the complex conjugate. From (B5a,b), we can get
and
Combining (B1), (B6a,b) and (B7), we could get
where the domain of harmonic order $m$ entering the summation is $[ -l,l ]$. Here, the $m=0$ contribution in (B8) is equally weighted with others. Both (B1) and (B8) are correct and they are equivalent to each other. In this work, we have adopted (B1) to compute the kinetic energy spectra.
Appendix C. Effective horizontal wavenumber and the integral length scale
The integral length scale converges at approximately the value of 1.5 in the spherical RBC, which is very close to the convergence value $(\sim 1.6)$ in the planar cases (Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018). This could be qualitatively explained from the kinetic energy spectra with respect to the so-called effective horizontal wavenumber. The effective horizontal wavenumber can be defined as follows.
In spherical coordinates, the spherical harmonics $Y_{l}^{m}(\theta, \phi )$ are the orthogonal eigenfunctions of the horizontal Laplacian operator in $\theta$ and $\phi$,
where $-l (l+1)/r^2$ is the eigenvalue of $Y_{l}^{m}(\theta,\phi )$. In analogy with the eigenvalues in the Cartesian coordinates,
where $k_{h}^{2}=k_{x}^{2}+k_{y}^{2}$ is the horizontal wavenumber. Comparing (C1) and (C2), we could then define an effective horizontal wavenumber in spherical coordinates as
Correspondingly, a horizontal wavelength can be invoked as
To obtain $E_{kin}^{'}(k_{H})$ spectra from $E_{kin}(l)$ spectra, we proceed as
Figures 35 and 36 show $E_{kin}(l)$ and $E_{kin}^{'}(k_{H})$ spectra for different $\eta$ at $Ra=7 \times 10^{7}$ and $Ra=3 \times 10^{8}$, respectively. Although $E_{kin}(l)$ spectra do not, all the $E_{kin}^{'}(k_{H})$ spectra collapse regardless of $\eta$ considered. The reason for the collapse is that it gets rid of the effect of the mid-depth radial coordinate, which is different for different $\eta$.
In planar RBC, Stevens et al. (Reference Stevens, Blass, Zhu, Verzicco and Lohse2018) found that $E_{kin}(k_{h})$ spectra collapse when aspect ratio $\varGamma \geq 32$. It can be reasoned that in spherical RBC, there also exists a threshold $\eta _{c}$ such that $E_{kin}^{'}(k_{H})$ spectra collapse when $\eta \geq \eta _{c}$. Since $E_{kin}^{'}(k_{H})$ spectra collapse for all values of $\eta$ considered in this work, this $\eta _{c}$ must be smaller than $0.2$. In the limit $\eta \to 1$, the behaviour of the system, including $E_{kin}^{'}(k_{H})$ spectra, will asymptotically approach the planar RBC with an infinite aspect ratio. The collapse in the $E_{kin}^{'}(k_{H})$ spectra for all $\eta \geq 0.2$, as shown in figures 35 and 36, is analogous to the collapse of the spectra in the planar case. Therefore, the integral length scale in simulations is more or less the same as the asymptotic value in planar RBC. This could qualitatively explain the reason why our $L_{int}$ converges to the similar value as observed in planar RBC.