1. Introduction
Thermal convection can be generated by a thermal gradient in combination with a fluid property that varies with temperature, giving rise to a volumetric body force. A common example of such a volumetric force is buoyancy, driving natural convection through a fluid’s temperature-dependent density under gravity. This phenomenon was investigated in depth by Lord Rayleigh in the early 19th century, providing a non-dimensional parameter, the Rayleigh number,
$Ra$
. Early experiments on Rayleigh–Bénard convection by Davis (Reference Davis1922a
,Reference Davis
b
) demonstrated the role of natural convection in heat transport, characterised by an increase in Nusselt number,
${\textit{Nu}}$
, with
$Ra$
, following the relation
${\textit{Nu}} \sim {\textit{Ra}}^\gamma$
with
$\gamma = 1/4$
. (Throughout the text, the symbol
$\sim$
denotes ‘scales as’ and not `order of magnitude’.) Since then, numerous power-law relationships have been proposed, see Grossmann & Lohse (Reference Grossmann and Lohse2000).
According to the principle formulated by Archimedes, a cooler, higher dense fluid descends, while a warmer, less dense fluid ascends, which is true for fluids without density anomalies such as in water (Rayleigh Reference Rayleigh1916; Yoshikawa et al. Reference Yoshikawa, Tadie Fogaing, Crumeyrolle and Mutabazi2013). A similar mechanism is observed in magnetic and dielectric fluids, where temperature-dependent properties such as magnetisation or electric permittivity can induce convective motion. This phenomenon has been widely investigated within the fields of magnetohydrodynamics (Zhang & Schubert Reference Zhang and Schubert2000; Hollerbach Reference Hollerbach2009; Aubert, Gastine & Fournier Reference Aubert, Gastine and Fournier2017), ferrohydrodynamics (Shliomis Reference Shliomis1974; Bahiraei & Hangi Reference Bahiraei and Hangi2015) and electrohydrodynamics (Roberts Reference Roberts1969; Turnbull Reference Turnbull1969; Stiles, Lin & Blennerhassett Reference Stiles, Lin and Blennerhassett1993; Mutabazi et al. Reference Mutabazi, Yoshikawa, Fogaing, Travnikov, Crumeyrolle, Futterer and Egbers2016; Szabo et al. Reference Szabo, Meyer, Meier, Motuz, Sliavin and Egbers2022; Meyer et al. Reference Meyer, Yoshikawa, Szabo, Meier, Egbers and Mutabazi2023).
Convection in rotating spherical shells plays a central role in understanding the fundamentals of large-scale transport within planetary interiors (Pothérat & Horn Reference Pothérat and Horn2024) and atmospheres (Scolan & Read Reference Scolan and Read2017). Rotating convection in fluid interiors occurs in the liquid iron cores of terrestrial planets and gas giants, and rapidly rotating cool stars feature turbulent Rayleigh–Bénard convection. Numerical simulations, such as those by Gastine, Wicht & Aubert (Reference Gastine, Wicht and Aubert2016), often incorporate radial forcing terms following the inverse-square law to model gravity caused by a central concentration of mass. This is experimentally challenging, particularly in the presence of central force fields generated by magnetic or electric sources. Alternative approaches have also been discussed by Busse (Reference Busse1978) in the context of rapidly rotating convection in a spherical shell. However, experiments involving radial forcing, especially those using magnetic or electric fields in non-isothermal fluids, are influenced by natural convection when carried out under terrestrial conditions. Therefore, these experiments are typically performed in microgravity environments such as in the geophysical flow simulator GeoFlow’ (Futterer et al. Reference Futterer, Krebs, Plesa, Zaussinger, Hollerbach, Breuer and Egbers2013; Zaussinger et al. Reference Zaussinger, Haun, Szabo, Travnikov, Al Kawwas and Egbers2020).
In microgravity, GeoFlow can induce convection of a temperature sensitive dielectric fluid confined in a spherical shell by using the temperature dependence of the fluid’s electrical permittivity rather than its density. The GeoFlow experiment consists of a spherical shell where a high-frequency AC potential is applied to the inner sphere and the outer sphere is grounded, creating an electric field across the gap. The resulting dielectrophoretic (DEP) force varies with the fluid’s temperature-dependent permittivity, producing an artificial gravity-like central force field. This mechanism drives thermo-electrohydrodynamic (TEHD) convection, analogous to classical Rayleigh–Bénard convection (see Turnbull Reference Turnbull1969; Landau, Lifshits & Pitaevskiĭ Reference Landau, Lifshits and Pitaevskiĭ1984; Gastine et al. Reference Gastine, Wicht and Aubert2016). The GeoFlow experiment, together with the properties of the aforementioned dielectric fluid, is designed to mimic terrestrial buoyancy analogous to those of planetary systems. With a uniformly heated inner shell and a cooled outer shell, GeoFlow provided a testbed for core and mantle convection. To capture the influence of planetary rotation, the apparatus is mounted on a rotating table. The rotation rate and the imposed temperature difference thus served as key control parameters for simulating different planetary interiors. More details on the microgravity experiment on the international space station are found in Beltrame, Egbers & Hollerbach (Reference Beltrame, Egbers and Hollerbach2003), Egbers et al. (Reference Egbers, Beyer, Bonhage, Hollerbach and Beltrame2003), Travnikov, Egbers & Hollerbach (Reference Travnikov, Egbers and Hollerbach2003) and Zaussinger et al. (Reference Zaussinger, Haun, Szabo, Travnikov, Al Kawwas and Egbers2020).
A second microgravity experiment, AtmoFlow, is currently prepared for atmospheric studies. In contrast to GeoFlow, it employs independently rotating spherical shells, enabling the study of non-isothermal spherical Taylor–Couette flow under the influence of the DEP force. In addition, the polar-angle dependence of the thermal boundary conditions is implemented to better mimic planetary atmospheres. The present study combines elements of both experiments: we adopt the thermal boundary conditions of GeoFlow while including the differential rotation characteristic of AtmoFlow. A subsequent study will then address the polar-angle-dependent boundary conditions in detail in another future numerical study. More details about the AtmoFlow experiment are found in Zaussinger et al. (Reference Zaussinger, Canfield, Froitzheim, Travnikov, Haun, Meier, Meyer, Heintzmann, Driebe and Egbers2019), Travnikov & Egbers (Reference Travnikov and Egbers2021) and Travnikov et al. (Reference Travnikov, Szabo, Gaillard and Egbers2025).
Spherical Taylor–Couette flow, initially proposed for second-generation dynamo experiments, exhibits flow instabilities reminiscent of planetary dynamo regions, which resemble convection columns (Wicht Reference Wicht2014). Boundary curvature effects, determined by the aspect ratio,
$\varGamma$
, defined as the ratio of the spherical shell radii, may play a significant role in influencing flow instabilities, as investigated by Egbers & Rath (Reference Egbers and Rath1995), Hollerbach (Reference Hollerbach2003) and Hollerbach et al. (Reference Hollerbach, Futterer, More and Egbers2004, Reference Hollerbach, Junk and Egbers2006). In the presence of differential rotation, an angular momentum imbalance along the meridional direction can be observed. The fluid flows radially outwards due to a radially outward-oriented momentum transport in the equatorial region and was first investigated by Haberman (Reference Haberman1962) using an analytical approach to define the formation of recirculating cells in the spherical gap. This equatorial flow develops during inner sphere rotation and reverses direction at higher outer sphere rates
$\varOmega$
. The phenomenon is further explored by Ovseenko (Reference Ovseenko1963) and later refined by Munson & Joseph (Reference Munson and Joseph1971) using a pseudo-analytic method based on Legendre polynomial series. A Stewartson shear layer can develop at the tangent cylinder that contacts the inner spherical surface (Stewartson Reference Stewartson1966). This layer can act as a source of instabilities (Wicht Reference Wicht2014), and the outward-directed radial jet may become unstable, exhibiting non-axisymmetric flow patterns at larger values of
$\Delta {\varOmega }$
. For both inner and outer spherical Taylor–Couette (sTC) flows, Guervilly & Cardin (Reference Guervilly and Cardin2010) demonstrated that this behaviour persists at moderate values of
$\varOmega$
. At large
$\varOmega$
, the so-called basic flow is dominated by the Coriolis force and approaches a quasi-two-dimensional geostrophic solution (Wicht Reference Wicht2014). In the case of differential rotation where
${\varOmega }\gg \Delta {\varOmega }$
the equatorial jet is suppressed. The parameter space governing rotation is extensive, leading to a wide range of flow regimes and patterns, as demonstrated by Nakabayashi & Tsuchida (Reference Nakabayashi and Tsuchida1988, Reference Nakabayashi and Tsuchida2005).
The closest study exploring the interaction of differential rotation and convection is given by Feudel & Feudel (Reference Feudel and Feudel2021). Their focus, however, is on the investigation of the bifurcation structure of spherical shell convection under the influence of differential rotation. This model of spherical Rayleigh–Bénard convection is important for the interior dynamics and evolution of planets and stars. In contrast to their approach, we incorporate a central force field proportional to
$1/r^5$
, where
$r$
is the radius (see Travnikov et al. Reference Travnikov, Egbers and Hollerbach2003; Futterer et al. Reference Futterer, Egbers, Dahley, Koch and Jehring2010), consistent with the AtmoFlow and GeoFlow frameworks. Our focus with this study is not on planetary dynamo regions or further sTC regime classification, but rather on investigating the interaction between the outward-directed radial jet caused by angular momentum transport in the equatorial region and the influence of the DEP force by combining both basic flow regimes. The interplay of these forces leads to complex heat transfer phenomena, which we elucidate through analysis of the AtmoFlow experiment. By relating the experimental regime to the observed pattern formation and associated heat transfer, we advance the understanding of convective processes.
Section 2 presents the model formulation, including the problem geometry, governing equations, numerical method and diagnostic quantities of interest, followed by resolution checks for the chosen parameter set. Numerical results are provided in § 3, where we numerically investigate both basic flows of dielectric convection in a spherical shell. We also examine the quasi-stationary states of sTC flow in contrast to the numerical results of Wicht (Reference Wicht2014). After these basic flows, we investigate the interaction of rotation and DEP forcing. A discussion about the findings is given in the context of flow regime classification in § 4. Final remarks and conclusions are summarised in § 5.
2. Model formulation
A schematic of the spherical shell geometry is shown, with inner and outer radii denoted by
$R_i$
and
$R_o$
, respectively. The spherical shell is heated at
$R_i$
to
$T_i$
and cooled at
$R_o$
to
$T_o$
. An electric field,
$\boldsymbol{E}$
, is applied in the radial direction,
$\boldsymbol{r}$
. The meridional and azimuthal directions are denoted by
$\varphi$
and
$\theta$
. The inner and outer shells can rotate independently with angular velocities of
$\varOmega _i$
and
$\varOmega _o$
.
2.1. Problem geometry
We consider a spherical shell filled with a Boussinesq dielectric fluid, where both the inner and outer shells can rotate differentially about the
$z$
axis. Convective motion is driven by buoyancy, resulting from a fixed temperature difference,
$\Delta T=T_i-T_o$
and an electric tension applied between the inner and outer shells, which have radii
$R_i$
and
$R_o$
, respectively. The boundaries of the shell are impermeable, no slip and maintained at constant temperatures. To induce TEHD convection, an alternating electric tension,
$V(t)$
, is applied at the inner shell, while the outer shell is grounded. A schematic of the problem geometry is shown in figure 1. As part of the experimental frameworks GeoFlow and AtmoFlow, we adopt the aspect ratio of
$\varGamma = {R_i}/{R_o} = 0.7$
. Further details regarding these experiments and the fluid properties used can be found in Zaussinger et al. (Reference Zaussinger, Haun, Szabo, Travnikov, Al Kawwas and Egbers2020) and Gaillard et al., (Reference Gaillard, Szabo and Egbers2023), respectively.
2.2. Governing equations
The underlying mechanism of the volumetric force acting on a dielectric fluid under the influence of an electric field is described by the Korteweg–Helmholtz force density, given by
\begin{equation} \boldsymbol{f}_{\!E}= \rho _e \boldsymbol {E} -\frac {1}{2}\boldsymbol {E}^2\boldsymbol{\nabla } \epsilon + \boldsymbol{\nabla }\Biggl [ \frac {\rho }{2} \biggl (\frac {\partial \epsilon }{\partial \rho }\biggl )_T \boldsymbol{E}^2 \Biggr ] ,\end{equation}
where
$\rho _e$
is the electric charge density,
$\boldsymbol{E}$
the electric field,
$\epsilon$
the electrical permittivity,
$\rho$
the density and
$T$
the temperature. The Korteweg–Helmholtz force density consists of the three main components given on the right-hand side of (2.1): the electrophoretic (Coulomb) force, the DEP force and the electrostrictive force. A detailed derivation of the force density can be found in Landau et al. (Reference Landau, Lifshits and Pitaevskiĭ1984). When an alternating electric field is applied between the spherical shells under the frequency condition
\begin{equation} f \ll \tau _e^{-1}, \tau _\nu ^{-1}, \tau _\kappa ^{-1}, \tau _d^{-1}, \tau _m^{-1} ,\end{equation}
where
$\tau _e$
,
$\tau _\nu$
,
$\tau _\kappa$
,
$\tau _d$
,
$\tau _m$
are the characteristic time scales of charge relaxation, viscous dissipation, thermal dissipation, migration and diffusion, respectively, the fluid cannot respond to the rapid variations in
$\boldsymbol{E}$
and the Coulomb force does not influence the fluid motion (Turnbull Reference Turnbull1968; Jones Reference Jones1979; Yoshikawa et al. Reference Yoshikawa, Tadie Fogaing, Crumeyrolle and Mutabazi2013; Szabo & Egbers Reference Szabo and Egbers2025). For incompressible dielectric fluids with stationary boundaries, the electrostrictive force does not affect fluid motion and can be lumped into the pressure force. Hence, the electrostrictive force can be assimilated to the hydrostatic pressure in the presence of an electric field. Consequently, only the DEP force drives fluid motion. As we consider an electrohydrodynamic (EHD) Boussinesq dielectric fluid confined between two concentric spherical shells, an equation of state can be formulated to describe the temperature-dependent variation in electrical permittivity given by
\begin{equation} \epsilon =\epsilon _0 \epsilon _r \left [1-e(T-T_0)\right ]\! ,\end{equation}
where
$\epsilon _0$
is the free space electrical permittivity and equals to
$8.854\:\times 10^{-12}$
F m
$^{-1}$
,
$\epsilon _r$
the relative permittivity,
$e$
the coefficient of thermal permittivity variation and
$T_0$
the reference temperature where the electrical permittivity is defined. Under the frequency condition specified in (2.2), the fluid remains initially electroneutral, with no accumulation of space charge (
$\rho _e = 0$
) throughout the development of the convective flow under the applied electric field. This allows the fluid motion within a dielectric layer to be modelled using (2.1), provided that the condition
$L \gg \lambda _D$
is met. This ensures that charge transport by convection from the diffusion layer into the bulk is negligible, where
$\lambda _D$
is the Debye length and
$L$
the spherical gap width denoted by
$R_o-R_i$
which serves as the system’s characteristic length (Yoshikawa et al. Reference Yoshikawa, Tadie Fogaing, Crumeyrolle and Mutabazi2013). Under the above assumptions, the alternating electric tension
$V(t)=\sqrt {2}V_0 \sin (2\pi f t)$
can be expressed by its root-mean-squared value,
$V_0$
, representing the root mean square voltage, defined by
$V(t):V_0=\sqrt {\langle V^2(t)\rangle }$
, where the angle brackets denote time averaging over one period of the voltage signal. The expression of the effective electric tension from the AC source is justified for
$f\tau _\nu \gtrsim 100$
, as supported by prior studies (Turnbull Reference Turnbull1969; Chandra & Smylie Reference Chandra and Smylie1972; Hart et al. Reference Hart, Toomre, Deane, Hurlburt, Glatzmaier, Fichtl, Leslie, Fowlis and Gilman1986; Futterer et al. Reference Futterer, Krebs, Plesa, Zaussinger, Hollerbach, Breuer and Egbers2013; Yoshikawa et al. Reference Yoshikawa, Tadie Fogaing, Crumeyrolle and Mutabazi2013, Reference Yoshikawa, Kang, Mutabazi, Zaussinger, Haun and Egbers2020; Szabo et al. Reference Szabo, Meier, Meyer, Barry, Motuz, Mutabazi and Egbers2021). Using the approximation of temperature-dependent permittivity given in (2.3), the remaining DEP force term can be reformulated to yield a time-averaged expression that drives fluid convection as follows:
\begin{equation} \boldsymbol{f}_{\textit{DEP}}=-\frac {1}{2}\boldsymbol {E}^2\boldsymbol{\nabla } \epsilon (T)=\frac {\epsilon _0\epsilon _r e}{2} \boldsymbol {E}^2 \boldsymbol{\nabla } (T-T_0). \end{equation}
We introduce the non-dimensional formulation of the conservation equation for an incompressible, non-isothermal, rotating dielectric EHD Boussinesq fluid, formulated in a reference frame co-rotating with the outer spherical shell, following the approach of Wicht (Reference Wicht2014) and Gaillard et al. (Reference Gaillard, Szabo and Egbers2023). Using the spherical shell gap width,
$L$
, and the thermal diffusion time scale,
$\tau _\kappa =L^2/\kappa$
, the conservation equations for continuity, momentum and energy and Gauss’s law read
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} =0 , \end{align}
\begin{align}& \partial _t {\boldsymbol{u}} + ({\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla }) {\boldsymbol{u}} =- \boldsymbol{\nabla } P + {\textit{Pr}} \,{\nabla }^2 {\boldsymbol{u}} - Pr\, 2 {\varOmega } \: \boldsymbol{e}_z \times {\boldsymbol{u}} - {\textit{Pr}} \, {\textit{Ra}}_{{E}} \vert {\boldsymbol{E}}\vert ^2 \boldsymbol{\nabla } T ,\end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad \partial _t T + ({\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla })T = {\nabla} ^2 T , \end{align}
\begin{align}&\qquad\qquad\qquad \boldsymbol {\nabla }\boldsymbol{\cdot }\left [(1-\gamma _e T)\boldsymbol{E}\right ]=0 \quad \text{with} \quad \boldsymbol{E}=- \boldsymbol{\nabla } \phi ,\end{align}
using the boundary condition of
\begin{equation} \left \{\!\!\!\!\!\!\begin{array}{lllll} &\boldsymbol{u}=\Delta {\varOmega }\, Pr\sin (\theta )\, r \,\boldsymbol{e}_\varphi , \quad &T=1, \quad &\phi =1, \quad \text{at} \quad &r=\displaystyle \frac {\varGamma }{1-\varGamma }=2.3\overline {33},\\ &\boldsymbol{u}=\boldsymbol{0}, \quad &T=0, \quad &\phi =0, \quad \text{at} \quad &r=\displaystyle \frac {1}{1-\varGamma }=3.3\overline {33}, \end{array}\right . \end{equation}
where
$\kappa$
is the thermal diffusivity,
$\boldsymbol{u}$
the velocity,
$P$
the modified pressure containing both the lumped electrostrictive force and the centrifugal force due to rotation,
$\gamma _e$
the thermoelectric parameters defined as
$e\Delta T$
,
$r$
the radial distance and
$\boldsymbol{e}_z$
and
$\boldsymbol{e}_\varphi$
are the unit vectors in the axial and azimuthal directions, respectively. (One has to note, that the non-dimensional set of equations are based on the thermal time scale, which introduces the Prandtl number
${\textit{Pr}}$
as a prefactor in the inner boundary condition.) In (2.8),
$\phi$
is the electric potential of the effective electric field,
$\boldsymbol{E}$
, scaled by the normalised root-mean-squared value of the electric tension,
$V_0$
. The electric field is given by Gauss’s law via the thermoelectric parameter,
$\gamma _e$
, coupling the electric field with the temperature variation in the fluid.
The set of system equations is governed by six non-dimensional control parameters, the aspect ratio,
$\varGamma$
, the thermoelectric parameter,
$\gamma _e$
, the Prandtl number,
${\textit{Pr}}=\nu /\kappa$
, the quantity of the dimensionless rotation rate
\begin{equation} {\varOmega }=\frac {{\varOmega _o}\:L^2}{\nu }, \end{equation}
denoted as the outer boundary rotation rate corresponding to the rotating frame of reference in which the outer boundary is stationary and satisfies the no-slip condition, the differential rotation rate
\begin{equation} \Delta {\varOmega }=\frac {({\varOmega _i}-{\varOmega _o})\:L^2}{\nu } \end{equation}
and DEP forcing expressed by the electric Rayleigh number
\begin{equation} {\textit{Ra}}_{{E}}=\frac {\epsilon _0\epsilon _r \gamma _e V_{0}^2}{2\rho \nu \kappa } ,\end{equation}
where
$\rho$
is the fluid’s density. These parameters govern the convective and rotational motion. The non-isothermal contribution to centrifugal buoyancy is neglected, as both the rotational and DEP forces are several orders of magnitude greater. The parameter
$\varOmega$
corresponds to the inverse Ekman number, representing the ratio of Coriolis to viscous dissipation. The convective flow can be interpreted via the Péclet number,
$ \textit{Pe}=UL/\kappa$
, which is related to the Reynolds number by
$ \textit{Re}= \textit{Pe}\,\textit{Pr}^{-1}$
, where
$U$
is the dimensional velocity magnitude.
An additional non-dimensional parameter is the Rossby number, defined by Hollerbach (Reference Hollerbach2003) for spherical flow as
\begin{equation} Ro = \Delta {\varOmega } / {\varOmega } ,\end{equation}
which represents the ratio of differential rotation to the outer boundary rotation. This parameter is useful for classifying flow regimes observed in sTC systems (Wicht Reference Wicht2014).
2.3. Numerical technique
Computational simulations are performed using the open source finite volume method within the OpenFOAM ecosystem. A custom solver was developed and benchmarked, with further details found in Gaillard et al. (Reference Gaillard, Szabo, Travnikov and Egbers2025), while only a brief summary is provided here. Equation (2.5) to (2.8) are implicitly solved within an inner correction loop using the Semi-Implicit Method for Pressure-Linked Equations algorithm (Ferziger, Perić & Street Reference Ferziger, Perić and Street2020). Symmetric matrices are built of all discretised equations such as the electric field and the pressure matrix and then solved via a geometric algebraic multigrid solver using a V-cycle. The temperature matrix, being asymmetric, is solved using a preconditioned bi-conjugate gradient solver. Hexahedral cells are used to mesh the gap region of the spherical shell, with slightly skewed areas introduced to improve pressure matrix convergence. A time-step solution is considered acceptable when the overall residual of the implicit matrix solver of all solved equations falls below
$10^{-4}$
. Time discretisation is performed using the Crank–Nicolson scheme with a blend factor of 0.9, while spatial discretisation employs a central difference scheme. This combination introduces minimal numerical dissipation and maintains acceptable numerical dispersion, resulting in rapid convergence within a few iterations.
2.4. Diagnostics
To assess the influence of non-dimensional control parameters on flow and temperature distributions, we introduce several diagnostic quantities. Domain and boundary probes are used along with averaging over at least 300 time steps at fixed intervals. The time step varies with the flow conditions, as different parameter choices lead to different strength in velocities. To ensure consistency across all simulations, the time step is scaled according to the maximum Courant–Friedrichs–Lewy number. This approach ensures we cover slow-rotating simulations and enables consistent comparisons with higher-rotation cases. Averaging is denoted by
$\langle \boldsymbol{\cdot }\rangle$
, with subscripts
$V$
,
$A$
and
$t$
indicating volume, surface area and temporal averaging, respectively. The non-dimensional kinetic energy density is given by
\begin{equation} \tilde {E}_K = \frac {1}{2} \big\langle \boldsymbol{\vert u\vert }^2 \big\rangle _{V,t} ,\end{equation}
capturing the full non-isothermal convective and rotational flow, both of which contribute to heat transfer across the concentric spherical shells. Since the flow velocity is driven by both rotation and DEP forcing, a temporally and spatially averaged Reynolds number is defined by
\begin{equation} Re_c=\sqrt {2{\tilde {E}_K}}\,\textit{Pr}^{-1} \end{equation}
to facilitate comparison between Reynolds numbers across different studies, such as in Grossmann & Lohse (Reference Grossmann and Lohse2000), particularly in relation to basic flow regimes.
The heat transport is characterised by the Nusselt number and represents a ratio of the total heat,
$Q_{{t}}$
, to the conductive heat,
$Q_{{c}}$
, across the spherical shell surfaces
\begin{equation} Nu=\frac {Q_{{t}}}{Q_{{c}}} ,\end{equation}
where
$Q_{{c}}$
can be expressed in analytical form as
\begin{equation} Q_{{c}}=\frac {4\pi \varGamma }{(\varGamma -1)^2} ,\end{equation}
and
$Q_{{t}}$
is computed numerically as the surface integral of the temperature gradient in the direction normal to the shell surface
\begin{equation} Q_{{t}}=\int \boldsymbol{\nabla }T \boldsymbol{n}\,\mathrm{d}A ,\end{equation}
where
$\boldsymbol{n}$
is the unit normal vector to the spherical shell surface. Besides the well-known expression of
${\textit{Nu}}$
, an inflow Nusselt number,
${\textit{Nu}}^q$
, is derived using the convective radial velocity,
$u_r$
, and the conductive radial temperature,
$\partial _r T$
. A full derivation of
${\textit{Nu}}^q$
is found in the Appendix B and it is represented by
\begin{equation} Nu^q=\frac {\left (\varGamma -1\right )^2}{4\pi \varGamma }\big \langle r^2\left [ uT- \partial _r T\right ]\big \rangle _{A,t} \end{equation}
by considering that
${\textit{Nu}}$
at the shell boundaries must equal
${\textit{Nu}}^q$
, since the gap contains neither sources nor sinks of heat. The time-averaged heat flux between the shells must therefore be equal. Consequently, all thermal energy leaving the inner shell must pass through the gap and reach the outer shell, ensuring energy conservation.
2.5. Parameter choice and resolution checks
The parameter selection in this study is based on the conditions of the AtmoFlow experiment, a rotating spherical shell with a confined dielectric fluid in which the inner and outer boundaries can rotate independently and be differentially heated. Further details of the experiment are provided in Zaussinger et al. (Reference Zaussinger, Haun, Szabo, Travnikov, Al Kawwas and Egbers2020). A summary of the choice and range of the forcing intensity of the control parameter space is provided in table 1, consistent with the fluid properties used in the AtmoFlow experiment and in the study of Zaussinger et al. (Reference Zaussinger, Haun, Szabo, Travnikov, Al Kawwas and Egbers2020).
Choice and range of the forcing intensity given by the control parameters.
Ensuring sufficient numerical resolution is essential in global convection simulations (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010; King, Stellmach & Aurnou Reference King, Stellmach and Aurnou2012; Gastine et al. Reference Gastine, Wicht and Aubert2016), as inadequate resolution can affect the accuracy of the scaling exponents obtained in asymptotic studies (see Amati et al. Reference Amati, Koal, Massaioli, Sreenivasan and Verzicco2005). One effective approach for validating the chosen resolution is to compare kinetic and thermal energy rates against the time-averaged quantities. King et al. (Reference King, Stellmach and Aurnou2012) proposed using the Nusselt number as a consistency check for boundary resolution, while the bulk flow can be characterised by the volume-averaged kinetic energy density. This approach is equivalent in comparing thermal and kinetic energy dissipation rates, which are also commonly used for validation. Reported deviations are typically within
$1\,\%$
(see King et al. Reference King, Stellmach and Aurnou2012) and up to
$3\,\%$
(see Gastine et al. Reference Gastine, Wicht and Aubert2016).
To assess the numerical resolution, we conducted a systematic mesh refinement study with uniform refinement in all spatial directions for the largest parameter values considered, namely
${\textit{Ra}}_{{E}}={1.77\times 10^{6}}$
,
${\varOmega }=399.1$
,
$\Delta {\varOmega }=456.1$
and
$\gamma _e={6\times 10^{-2}}$
. The results are shown in figure 2, and the corresponding mesh parameters with results of
${\textit{Nu}}$
and
$\tilde {E}_K$
are summarised in table 3 in the appendix. This approach is equivalent to Stevens, Verzicco & Lohse (Reference Stevens, Verzicco and Lohse2010), King et al. (Reference King, Stellmach and Aurnou2012) and Gastine et al. (Reference Gastine, Wicht and Aubert2016). Figure 2 illustrates
${\textit{Nu}}$
at several grid resolutions. For grids exceeding
$1.5\times 10^6$
and
$5.2\times 10^6$
hexahedral cells, respectively,
${\textit{Nu}}$
and
$\tilde {E}_K$
show a relative change of the investigated quantities to the finest mesh below 1 %. Hence, this low errors results further confirm that the employed numerical resolutions are sufficient and fairly accurate.
Grötzbach (Reference Grötzbach1983) proposed that direct numerical simulation must resolve the Kolmogorov scale,
$\eta$
, and the Batchelor scale,
$\eta _T$
. The maximum grid spacing
$h=\max (\Delta x,\Delta y,\Delta z)$
should satisfy
\begin{align} h \leqslant \pi \eta &= \pi \left (\frac {{\textit{Pr}}^2}{{\textit{Ra}}_E\,\textit{Nu}}\right )^{1/4} \quad \text{for} \quad {\textit{Pr}} \leqslant 1, \end{align}
\begin{align} h \leqslant \pi \eta _T &= \pi \left (\frac {1}{{\textit{Ra}}_E\,\textit{Pr}\,\textit{Nu}}\right )^{1/4} \quad \text{for} \quad {\textit{Pr}} \geqslant 1, \end{align}
which are the commonly adopted resolution criteria (see Stevens et al. Reference Stevens, Verzicco and Lohse2010). Here, the standard Rayleigh number is replaced by the electric Rayleigh number,
${\textit{Ra}}_E$
. Since
${\textit{Pr}}\,\gg \,1$
in our case, the Batchelor scale is the resolution requirement, yielding
$h_{\eta _T}=0.024$
and a minimum of 7.4 million cells in the spherical shell volume. To adequately capture thermal and kinetic energy dissipation as well as boundary layers, we employ 12 million hexahedral cells, which provides a good compromise between computational cost and accuracy, with errors below 1 %.
Mesh independency for integral quantities: Nusselt number,
${\textit{Nu}}$
, (
) and kinetic energy density
$\tilde {E}_K$
(
).
3. Results
3.1. Conductive base state
When
${\textit{Ra}}_{{E}}=0$
, i.e. no electric potential is applied and both
${\varOmega }=0$
and
$\Delta {\varOmega }=0$
(no rotation), there is no flow (
$\boldsymbol{u}=\boldsymbol{0}$
) and hence the temperature in the spherical shell is at its conductive base state, which can be expressed by (2.17).
3.2. Basic flow
Before analysing the flow under combined rotational and DEP forcing, we first examine the basic flow regimes under each forcing separately. The two basic flows considered are DEP force-driven TEHD convection and non-isothermal sTC, as detailed in the following two sections.
3.2.1. Thermo-electrohydrodynamic convection in spherical shell
The basic flow of DEP force-driven convection in the spherical shell is considered under no rotation conditions, with both
${\varOmega } = 0$
and
$\Delta {\varOmega } = 0$
. The convective behaviour in pure TEHD convection within the spherical shell aligns with previous findings by Futterer et al. (Reference Futterer, Egbers, Dahley, Koch and Jehring2010) and Szabo et al. (Reference Szabo, Gaillard, Zaussinger and Egbers2023), where plume and sheet patterns are observed.
Isosurface plots of temperature,
$T$
, radial velocity,
$u_r$
, and meridional velocity,
$u_\theta$
, are shown in the top, middle and bottom rows, respectively, across equatorial, meridional and radial cross-sections. To illustrate heat transport across both boundaries, the heat flux,
$q$
, is plotted at the inner and outer shells alongside the temperature cross-sections in the top row. In addition to the cross-sectional views,
$u_r$
and
$u_\theta$
are also shown on a spherical shell located at a quarter radius, defined as
$R_{1/4} = (1 + 3\varGamma ) / (4 - 4\varGamma )$
. Colour scales are provided below each panel. All plots correspond to stationary shells, with increasing electric Rayleigh numbers,
${\textit{Ra}}_{{E}}$
, of
$1.59 \times 10^4$
,
$7.95 \times 10^4$
and
$1.77 \times 10^6$
for the (a,d,g), (b,e,h) and (c, f,i) columns, respectively.
The basic flow of DEP force-driven convection at the final simulation time is shown in figure 3. The top row displays the temperature field in the spherical shell along equatorial, meridional and radial planes. To highlight thermal transport, the heat flux,
$q$
, is plotted at the inner and outer boundaries, which in fact represents the temperature gradient. The middle row presents the radial velocity,
$u_r$
, on the same planes, together with a surface plot at the quarter gap,
$R_{1/4}=(1+3\varGamma )/(4-4\varGamma )$
. The bottom row shows the meridional velocity,
$u_\theta$
, on the corresponding planes used in the middle row.
At relative small electric Rayleigh numbers the flow can be characterised by large convective outward-oriented plumes, see figure 3(a) and a radial outward flow of hot fluid with an inward cold flow feeding the plumes, seen in figure 3(d), observed at the mid-latitudes and also observed in the meridional plane in figure 3(g). As the electric Rayleigh number increases, the large-scale plumes become increasingly irregular. Near the inner boundary, hot plumes rise from the inner boundary, while cold plumes sink from the outer boundary, as shown in figure 3(b). With stronger convection, the radial outflow of hot fluid (figure 3 c) develops irregularities, accompanied by recirculation of cold fluid from the outer boundary also evident in the meridional flow in figure 3(h). For the largest electric Rayleigh number, deep convection develops, as shown in figure 3(c). The heat is efficiently mixed, producing irregular, chaotic flows reminiscent of turbulent Rayleigh–Bénard convection. Within the gap, radial and meridional circulations contribute to strong mixing, with the maximum radial velocity nearly twice that of the meridional component (see figure 3 f,i).
Since convection is solely driven by the DEP force, quantified by the control parameter,
${\textit{Ra}}_{{E}}$
, (see figure 4
a), the Nusselt number,
${\textit{Nu}}$
, increases with
${\textit{Ra}}_{{E}}$
, while the influence of
$\gamma _e$
on heat transport remains marginal. A scaling law between the two parameters can be expressed as
\begin{equation} Nu \sim {\textit{Ra}}_{{E}}^{1/3} ,\end{equation}
which is consistent with the convective flow behaviour described in Grossmann & Lohse (Reference Grossmann and Lohse2000). A similar scaling relation can be applied to
$\tilde {E}_K$
as a function of
${\textit{Ra}}_{{E}}$
, as shown in figure 4(b), yielding
\begin{equation} \tilde {E}_K \sim {\textit{Ra}}_{{E}}^{1.16} ,\end{equation}
where the scaling of
$\tilde {E}_K$
can also be expressed in terms of the convective Reynolds number scaling by
$Re_c \sim {\textit{Ra}}_{{E}}^{0.58}$
, which is similar to the exponent of the `wind of turbulence’ introduced by Grossmann & Lohse (Reference Grossmann and Lohse2000) given by
$ \textit{Re}_c\sim {\textit{Ra}}^{0.44}$
.
Panels (a) and (b) respectively show the Nusselt number,
${\textit{Nu}}$
, and the kinetic energy density,
$\tilde {E}_K$
, as functions of electric Rayleigh number,
${\textit{Ra}}_{{E}}$
. The symbols (
$\circ$
), (
$\square$
), (
$\diamond$
) and (
$\times$
) represent
$\gamma _e$
values of 0.006, 0.015, 0.03 and 0.06. Temporal and spatial average associated error bars are smaller than the symbol size and thus omitted for clarity. The dashed lines in panels (a) and (b) represent the power-law scalings defined in (3.1) (
) and (3.2) (
), respectively.
3.2.2. Non-isothermal spherical Taylor–Couette flow
Isosurface plots of temperature,
$T$
, radial velocity,
$u_r$
, and meridional velocity,
$u_\theta$
, are shown in the top, middle and bottom rows, respectively, across equatorial, meridional and radial cross-sections. To illustrate heat transport across both boundaries, the heat flux,
$q$
, is plotted at the inner and outer shells alongside the temperature cross-sections in the top row. In addition to the cross-sectional views,
$u_r$
and
$u_\theta$
are also shown on a spherical shell located at a quarter radius, defined as
$R_{1/4} = (1 + 3\varGamma ) / (4 - 4\varGamma )$
. Colour scales are provided below each panel. All plots set the electric Rayleigh number to zero
${\textit{Ra}}_{{E}}=0$
with increasing rotation rate.
Panel (a) shows the relationship between
$\Delta {\varOmega }$
and
$\varOmega$
in a modified regime diagram based on the classification by Wicht (Reference Wicht2014). The dashed line (
) represents the Stewartson stability condition, defined by the Rossby number
${\textit{Ro}} = \Delta {\varOmega } / {\varOmega } = 1$
. The solid line (
) marks the non-axisymmetric stability threshold for
$\varGamma = 0.35$
, while line (
) provides an estimated boundary for non-axisymmetric stability at
$\varGamma = 0.7$
, inferred from the computed axisymmetric states for
${\textit{Ro}}\in [0.429,1.14]$
(
) and
${\textit{Ro}}\in ]0,0.429[\cup ]1.14,\infty [$
(
) and non-axisymmetric cases (
). Panel (b) illustrates the
$\varGamma$
-dependence of
${\varOmega }_i$
, adapted from Egbers & Rath (Reference Egbers and Rath1995), highlighting the transitions from laminar flow to secondary wave and turbulent regimes for
${\varOmega }=0$
.
In the basic flow of sTC systems at
${\textit{Ra}}_{{E}}=0$
, pattern formation can be classified following the approach of Wicht (Reference Wicht2014). Unlike Wicht’s configuration, the thermal boundary conditions in (2.9) are retained and the temperature field in (2.7) is explicitly solved. Consequently, the basic flow remains non-isothermal even in the absence of the DEP force.
The results of the basic flow of sTC is shown in figure 5 using the same cross-sections as introduced in § 3.2.1 and illustrates the thermal pattern formation associated with the shear flow when both spherical shells rotate in the radial and meridional directions at the same Rossby number. For relatively weak differential forcing (
${\varOmega } = 5.25$
,
$\Delta {\varOmega } = 3.72$
), the temperature isotherms bend more outward at the equator than at the poles, see figure 5(a). The shear flow advects hot fluid from the equatorial region toward the cold outer shell, where it splits symmetrically in the poleward directions, forming an equatorial plume. To satisfy continuity, the return flow emerges from both poles toward the inner shell, as seen in the radial and meridional velocity fields in figures 5(d) and 5(g). This shear-driven circulation enhances angular momentum transport and contributes to the heat flux across both shells. Increasing the rotation rate (
${\varOmega }=33.6$
,
$\Delta {\varOmega }=23.5$
) enhances the bending of temperature isotherms into plume-like structures, strengthens the radial outward-directed flow and intensifies the meridional circulation, with poleward motion at the outer shell and recirculation at the inner shell, see figures 5(b), 5(e) and 5(h). Since the flow is still predominantly antisymmetric, there are no plumes distributed in the meridional direction. These effects lead to a higher heat flux across both boundaries. For the strongest sTC forcing parameter (
${\varOmega } = 243.3$
,
$\Delta {\varOmega } = 172.1$
; figure 5
c), the flow symmetry breaks and the pattern shifts toward higher latitudes, forming columnar structures reminiscent of the tangent cylinder in sTC flow (Hoff & Harlander Reference Hoff and Harlander2019). At the equatorial plane, the radial velocity shows a reoccurring pattern with stable mode numbers, which also appears in the meridional velocity as a bent structure, see figures 5(f) and 5(i). Overall, the dynamics is now dominated by radial motion compared with the meridional component. A direct comparison with recent studies is not straightforward because of the different parameter regimes considered. Our simulations are performed at
$Ek=\varOmega ^{-1}=10^{-2}$
and
${\textit{Ro}} \approx 1$
, with an aspect ratio of
$0.7$
. In contrast, Wicht (Reference Wicht2014) investigated transitions over a much broader parameter range,
$Ek \in [10^{-8},4]$
, using an aspect ratio of
$0.3$
. Similarly, Barik et al. (Reference Barik, Triana, Hoff and Wicht2024) focused on more rapidly rotating regimes, with
$Ek \in [10^{-6},10^{-4}]$
and
${\textit{Ro}} \in [0.5,5]$
, also for an aspect ratio of
$0.3$
. These differences in forcing and geometry prevent a quantitative comparison as their results favour more the implications for angular momentum transport in astrophysical systems such as protoplanetary discs, planets and stars.
From a general point of view one can conclude, under weak rotational forcing, that viscous effects dominate, and the first instability appears as a radially outward-directed jet at the equator when
${\varOmega } = 0$
. Guervilly & Cardin (Reference Guervilly and Cardin2010) showed that this behaviour persists for moderate values of
$\varOmega$
. When
${\varOmega }\gt \Delta {\varOmega }$
, a tangent cylinder with a Stewartson shear layer develops. We adopt Wicht’s classification, illustrated in figure 6(a), where the solid line represents the non-axisymmetric stability line for
$\varGamma = 0.35$
. Above this line, non-axisymmetric instabilities arise. Egbers & Rath (Reference Egbers and Rath1995) demonstrated that boundary curvature affects transitions between laminar flow, secondary wave states and turbulence. Figure 6(e) shows the influence of
$\varGamma$
on regime classification. Consequently, for
$\varGamma = 0.7$
, non-axisymmetric instability is expected to occur at lower values of
$\Delta {\varOmega }$
compared with systems with
$\varGamma = 0.35$
. The results shown in figure 6(a) indicate that the onset of non-axisymmetric instability occurs earlier, where (
, 
) represents time-invariant axisymmetric cases and (
) denotes as non-axisymmetric transient cases. The thin dashed line marks the Stewartson stability condition, defined by
${\textit{Ro}}=\Delta {\varOmega } / {\varOmega }=1$
. Our results are consistent with those of Wicht (Reference Wicht2014), showing that, in the viscous regime, for large
${\textit{Ro}}$
(to the left of the Stewartson stability line), a radially outward-directed equatorial jet forms, while for small
${\textit{Ro}}$
(to the right of the Stewartson stability line), a tangent cylinder with a Stewartson shear layer dominates. For
${\textit{Ro}}=1$
a mix between equatorial jet and Stewartson shear layer is observed. In this study, we considered only the basic axisymmetric states (
) and investigated the influence of the DEP force.
In non-isothermal sTC flow, angular momentum transport can drive convection through shell rotation. Due to the meridional inclination angle,
$\theta$
, rotational forces are stronger near the equator than near the poles. Within the investigated parameter space, this imbalance drives the transport of warmer fluid, heated at the inner shell, toward the cooler outer shell. This occurs either via a radially outward jet at the equator or along the tangent cylinder through the Stewartson shear layer, depending on the value of
${\textit{Ro}}$
. In the regime near
${\textit{Ro}}\approxeq 1$
,
${\textit{Nu}}$
can be related to the differential rotation
$\Delta {\varOmega }$
, as illustrated in figure 7(a). The heat transport in this regime follows the scaling relation
\begin{equation} Nu \sim \Delta {\varOmega }^{1/2} ,\end{equation}
which corresponds to a larger scaling exponent,
$\gamma$
, than that in (3.1). It is also significantly larger than the value reported for natural convection for small
$Ra$
by Davis (Reference Davis1922a
), where
$\gamma = 1/4$
. The kinetic energy density is plotted in figure 7(b) and scales as
\begin{equation} \tilde {E}_K \sim \Delta {\varOmega }^{2} ,\end{equation}
indicating a significantly higher exponent
$\gamma =2$
compared with that observed in TEHD convection (see, (3.2)). One has to pay additional attention to the limitation of the scaling relations as the heat transport is expected to depend on the radial, and therefore on the poloidal, flow component. The relation between poloidal and toroidal flow scalings has been discussed in Wicht (Reference Wicht2014). Our results form (3.1) and (3.2) indicate a weaker scaling compared with (3.3) and (3.4), which reflects the role of heat dissipation as its origin. Moreover, one can expect that the scaling for sTC will vary across different parts of the instability; in particular, non-axisymmetric flows may contribute once they become significant.
) and (
).
3.3. Thermo-electrohydrodynamic convection in spherical Taylor–Couette flow
In this section, both basic flow regimes are combined to arrive at DEP force-driven TEHD convection in sTC. We first show the spatial and temporal evolution which is followed by analysis of the heat transfer. A discussion is provided in the following section.
3.3.1. Spatial and temporal evolution
Isosurface plots of temperature,
$T$
, radial velocity,
$u_r$
, and meridional velocity,
$u_\theta$
, are shown in the top, middle and bottom rows, respectively, across equatorial, meridional and radial cross-sections. To illustrate heat transport across both boundaries, the heat flux,
$q$
, is plotted at the inner and outer shells alongside the temperature cross-sections in the top row. In addition to the cross-sectional views,
$u_r$
and
$u_\theta$
are also shown on a spherical shell located at a quarter radius, defined as
$R_{1/4} = (1 + 3\varGamma ) / (4 - 4\varGamma )$
. Colour scales are provided below each panel. All plots correspond to a fixed rotation rate of
${\varOmega } = 33.6$
and differential rotation
$\Delta {\varOmega } = 23.5$
, with increasing electric Rayleigh numbers,
${\textit{Ra}}_{{E}}$
, of
$1.59 \times 10^4$
,
$7.95 \times 10^4$
and
$1.77 \times 10^6$
for the left, middle and right columns, respectively.
To illustrate the diversity of dynamical regimes resulting from various combinations of basic flows, we present a set of simulations with three different
${\textit{Ra}}_{{E}}$
values while keeping
${\varOmega }=33.6$
and
$\Delta \varOmega =23.5$
fixed. Figure 8, top row, displays the temperature field within the spherical shell using equatorial, meridional and radial cross-sections, as outlined in § 3.2.1.
In figure 8(a), for
${\textit{Ra}}_{{E}}=1.59\times 10^4$
, a radially outward-directed plume at the equator becomes evident that extends from the inner shell boundary to the outer shell. Upon reaching the outer boundary, the plume splits symmetrically toward both poles, attempting to form the tangent cylinder via the Stewartson shear layer. This is also indicated by the heat flux,
$q$
, as figure 8(d) shows a radially directed outward equatorial jet transporting hot fluid from the inner to the outer shell along the equatorial plane. Upon reaching the outer boundary, the flow bifurcates symmetrically toward both poles, forming the tangent cylinder through the Stewartson shear layer, as illustrated in figure 8(g). To satisfy continuity, a radially inward-directed flow develops at each pole (see figure 8
d), inducing a meridional return flow at the inner boundary directed towards the equator, as shown in figure 8(g). When
${\textit{Ra}}_{{E}}$
increases to
$7.95\times 10^4$
, see figure 8(b), the equatorial plume, initially directed radially outward, splits into multiple plumes and breaks the equatorial symmetry. Instead, the temperature field reveals symmetry breaking, indicating irregular convection. The radial velocity, depicted in figure 8(e), shows the equatorial plume dividing into distinct convective modes, forming inclined patterns along the spherical surface at
$R_{{1}/{4}}$
reminiscent of fish-bone structures. The meridional velocity field, see figure 8(h), exhibits a clear poleward flow, further emphasising the asymmetry. Both radial and azimuthal velocity profiles display strong poleward drafts at the outer boundary and irregular inward flows at the poles, accompanied by a meridional return flow at the inner boundary directed towards the equator. When
${\textit{Ra}}_{{E}}$
increases further to
$1.77\times 10^6$
, the temperature distribution in figure 8(c) indicates deep convection in the bulk where the fluid is almost isothermalised. Temperature fluctuations appear as small thermal plumes emerging from a thin boundary layer, forming long, thin sheets extending radially outward throughout the spherical domain. The boundary temperature gradient,
$q$
, also reflects deep convection, with plume and sheet structures indicating irregular convective motion within the spherical shell. The radial and meridional velocity fields shown in figures 8(f) and 8(i) confirm strong convection, suggesting irregular, chaotic radial flows reminiscent of turbulent Rayleigh–Bénard convection.
The basic flow regimes in figures 3 and 5 exhibit a radial outward-directed flow at the equator, however, both have a different driving mechanism the flow. In the non-rotating TEHD case, the flow is driven by the DEP force, whereas in the differentially rotating shell with
${\textit{Ra}}_E=0$
, shear flow induces angular momentum transport. Since both processes act in the same direction at the equator, the radial velocity is amplified. Moving poleward, the shear and associated angular momentum transport weaken. At mid-latitudes, shear-induced recirculation can suppress DEP-driven plume formation when the forcing is weak. Near the poles, where shear flow due to differential rotation is minimal, the strong DEP force dominates and produces polar plumes, as illustrated in the second column of figure 8.
While the first case with
${\textit{Ra}}_{{E}}=1.59\times 10^4$
yields a time-invariant solution, the other two cases displayed in figure 8 exhibit a time-dependent behaviour, which is illustrated by the space–time plots in figure 9, after the simulation has integrated long enough to reach a statistically equilibrated solution of the initial state given by the boundary conditions of (2.9). To assess temporal variations in temperature along the azimuthal direction,
$\varphi$
, (top row), and the meridional direction,
$\theta$
, (bottom row) the three-quarter gap defined by
$R_{3/4} = (3 + \varGamma ) / (4 - 4\varGamma )$
is used.
Figures 9(a) and 9(c) present the azimuthal and meridional temperature distributions, respectively, over time for
${\textit{Ra}}_{{E}}=7.95\times 10^4$
. The azimuthal space–time plot shows convective structures spanning from the inner to the outer shell, with thermal vacillations developing in the azimuthal direction. This behaviour is even clearer in the meridional space–time plot, where the vacillations appear irregular and are accompanied by fluctuating patterns near the poles. For the strong forcing case of
${\textit{Ra}}_{{E}}=1.77\times 10^6$
, the space–time plots in both the meridional and azimuthal directions (figures 9
b and 9
d, respectively) indicate deep convection throughout the bulk. The hot and cold patterns are randomly distributed and no coherent structures are clearly identifiable. Unlike Feudel & Feudel (Reference Feudel and Feudel2021), the present study does not perform a bifurcation analysis, nor does it resolve the complete bifurcation structure or unstable branches via path-following methods. Furthermore, whereas previous models assumed a linear gravity profile, we employ a radially scaled gravity proportional to
$1/r^{5}$
, as discussed above, which is limited to the electric field decay in TEHD convection. Despite similarities in set-up, the scope of the two studies differs substantially and a clear comparison is not applicable. From a qualitative perspective, the flow reported by Feudel & Feudel (Reference Feudel and Feudel2021) exhibits radially outward-oriented plumes that resemble the TEHD-induced plumes observed in the polar regions of the present study. Moreover, their equatorial cases display radially outward patterns that are deflected in a manner comparable to our case with
${\textit{Ra}}_{{E}}={7.95\times 10{^4}}$
,
$\varOmega ={33.6}$
and
$\Delta \varOmega ={23.5}$
. This deflection arises from the combined influence of the Coriolis force and angular momentum conservation.
Space–time plots of temperature,
$T$
, in the azimuthal direction,
$\varphi$
for
$\theta =0$
, at the equatorial cross-section (a,b), and in the meridional direction,
$\theta$
for
$\varphi =0$
, at the pole-to-pole cross-section (c,d), evaluated at the three-quarter gap defined by
$R_{3/4} = (3 + \varGamma ) / (4 - 4\varGamma )$
.
3.3.2. Heat transfer and scaling
After presenting three representative flow cases with progressively increasing DEP force in a steady-state non-isothermal sTC configuration, this study proceeds to a more quantitative analysis reporting on the diagnostics previously introduced, namely, the Nusselt number,
${\textit{Nu}}$
, and the kinetic energy density,
$\tilde {E}_k$
. This approach follows the methodology used in figures 4 and 7 to evaluate the basic flow states of TEHD convection and non-isothermal sTC, now extended by the two individual forcing parameters of
$\Delta {\varOmega }$
and
${\textit{Ra}}_{{E}}$
in figure 10.
Figure 10(a) shows
${\textit{Nu}}$
as a function of the
${\textit{Ra}}_{{E}}$
and rotation rate
$\Delta {\varOmega }$
. The results indicate that
${\textit{Nu}}$
increases independently with
${\textit{Ra}}_{{E}}$
and
$\Delta {\varOmega }$
, as well as when both parameters increase relatively simultaneously, consistent with the underlying basic flow states. However, for larger values of
${\textit{Ra}}_{{E}}$
, the influence of
$\Delta {\varOmega }$
on
${\textit{Nu}}$
is less pronounced and the increase in
${\textit{Nu}}$
is more dominated by larger
${\textit{Ra}}_{{E}}$
values regardless of any further increases in
$\Delta {\varOmega }$
. The kinetic energy density shown in figure 10(b) presents a different trend, where
$\tilde {E}_K$
increases with both
$\Delta {{\varOmega }}$
and
${\textit{Ra}}_{{E}}$
simultaneously. This suggests an alternative heat transport mechanism, which is further examined through the conductive and convective components of
${\textit{Nu}}^q$
discussed later.
The Nusselt number,
${\textit{Nu}}$
, and kinetic energy density,
$\tilde {E}_k$
as a functions of
$\Delta {\varOmega }$
and
${\textit{Ra}}_{{E}}$
in (a) and (b), respectively.
Panels (a) and (b) show the Nusselt number,
${\textit{Nu}}$
, and kinetic energy density,
$\tilde {E}_k$
, respectively, as a functions of
${\textit{Ra}}_{{E}}$
for three values of
$\Delta {\varOmega }={2.25}(\times )$
$\Delta {\varOmega }={3.90}\times 10{^1}(\square )$
and
$\Delta {\varOmega }={1.04}\times 10{^2}(\circ )$
. Results for the non-rotating case (
${\varOmega }=0$
and
$\Delta {\varOmega }=0$
) are included by (
) markers as a reference case where sTC flow is absent.
For clarity, figure 11(a) presents four representative cases with varying rotation rates under imposed DEP forcing. The plots show that, at low
${\textit{Ra}}_{{E}}$
, spherical shell rotation is capable of enhanced heat transport. As
${\textit{Ra}}_{{E}}$
increases, given the selected parameters, the flow becomes increasingly dominated by TEHD convection. In this regime, differential rotation cases converge toward the heat transfer characteristics by the basic flow of TEHD convection. Figure 11(b) also indicates the difference in the heat transport mechanism where
$\tilde {E}_K$
remains largely unaffected at moderate and high rotation rates. Only at the smallest rotation rate does TEHD convection contribute significantly to the kinetic energy when
${\textit{Ra}}_{{E}}$
is increasing. For moderate rotation, a noticeable increase in kinetic energy occurs only when
${\textit{Ra}}_{{E}}$
exceeds approximately
$3\times 10{^6}$
.
Panels (a) and (b) show the Nusselt number,
${\textit{Nu}}$
, and kinetic energy density,
$\tilde {E}_k$
, respectively, as a functions of
$\Delta {\varOmega }$
for three values of
${\textit{Ra}}_{{E}}={1.59\times 10^4}(\times )$
,
${\textit{Ra}}_{{E}}={3.18\times 10^5}$
axisymmetric (
$\square$
), non-axisymmetric (
) and
${\textit{Ra}}_{{E}}={8.84\times 10^5}$
axisymmetric (
$\circ$
), non-axisymmetric (
). Results for the case where TEHD convection is absent (
${\textit{Ra}}_{{E}}=0$
) are included by (
) markers as a reference case. The capital letters in (a) indicate the flow regimes: (A) represents the TEHD convection-dominated regime; (B) is the transitional regime, where both DEP force and angular momentum contribute comparably to heat transport; and (C) denotes the non-isothermal sTC regime, where strong shear flows due to differential rotation enhance heat transport across the spherical shell.
Figure 12 shows a similar approach by considering the case of an imposed
$\Delta {\varOmega }$
for four representative
${\textit{Ra}}_{{E}}$
. In figure 12(a), the Nusselt number,
${\textit{Nu}}$
, is shown as a function of
$\Delta \varOmega$
. For small
$\Delta \varOmega$
, heat transport is primarily driven by TEHD convection, referred to as regime (A). At moderate values of
$\Delta {\varOmega }$
, the heat transfer decreases, indicating a shift from radially driven DEP convection to shear-dominated flow induced by differential rotation. This transition produces an outward-directed equatorial jet, sustained by poleward flow along the inner sphere and returning toward the outer shell near the poles, influencing the DEP forced plume generation at mid-latitudes. This intermediate behaviour characterises the transition region, referred to as regime (B), and is identified by a slope fit when
${\textit{Nu}}$
decreases by more than 1 % relative to the preceding increment in
$\Delta {\varOmega }$
. At large
$\Delta {\varOmega }$
, the flow is dominated by strong shear flow and approaches the basic non-isothermal sTC state, corresponding to regime (C). The transition from regime (B) to (C) is marked by the recovery of
${\textit{Nu}}$
, indicating that DEP-driven TEHD convection vanishes and heat transport is governed primarily by the shear flow. The plot also includes non-axisymmetric rotation cases at large
$\Delta {\varOmega }$
(see figure 6
a). Hence, regime (C) is reached once
${\textit{Nu}}$
fully returns to its pre-reduction value or matches the level obtained at
${\textit{Ra}}_{{E}}=0$
. The kinetic energy density,
$\tilde {E}_k$
, shown in figure 12(e), indicates that convective heat transport is dominated by the DEP force at small
$\Delta {\varOmega }$
, converging to the non-isothermal sTC state with
${\textit{Ra}}_{{E}} = 0$
. A comparison of the observations reveals that the decrease in
${\textit{Nu}}$
is not reflected in
$\tilde {E}_k$
as no decrease in kinetic energy in trasition regime (B) is observable arising from a change in the underlying convective transport mechanism. This transition is examined in the following section, where the regime diagram in figure 12(a) is extended to include all investigated cases.
The current results show that increasing the relative rotation rate
$\Delta \varOmega$
leads to a decrease in the Nusselt number, while the kinetic energy still increases. Since the Rossby number remains of order unity throughout this investigation, the ratio between the solid-body rotation rate
$\varOmega$
and the relative rotation
$\Delta \varOmega$
scales proportionally. We aim to isolate the effect of solid-body rotation, which is known to influence the Nusselt number, as demonstrated by Gastine et al. (Reference Gastine, Wicht and Aubert2016), and to determine for the present configuration whether the observed behaviour arises primarily from solid-body rotation or from differential rotation. To this end, we investigate the behaviour of the Nusselt number and kinetic energy for the case
$\Delta \varOmega = 0$
with
$\varOmega \neq 0$
.
Although this analysis lies outside the primary scope of the present study, we performed six simulations covering the relevant parameter space for three different values of
${\textit{Ra}}_E$
. Because the convective flow is deflected orthogonally to the rotation axis, i.e. into the direction of the rotation axis, an influence on heat transport is observed. Figure 13(a) shows this effect for the three selected values of
${\textit{Ra}}_E$
: increasing rotation reduces the Nusselt number, most prominently for the smallest
${\textit{Ra}}_E$
, where
${\textit{Nu}}$
decreases from 2.9 at
$\varOmega = 12$
to 1.5 at
$\varOmega = 400$
. This reduction in
${\textit{Nu}}$
is accompanied by a decrease in kinetic energy, as shown in figure 13(b), with increasing
$\varOmega$
, noting that the system is formulated in the rotating frame of reference.
For larger
${\textit{Ra}}_E$
, the Nusselt number remains nearly constant, within the range
$[6.5,\,6.8]$
for
${\textit{Ra}}_E = {3.18\times 10^{5}}$
and
$[8.9,\,9.2]$
for
${\textit{Ra}}_E = {8.84\times 10^{5}}$
, indicating a weaker influence of the Coriolis force. Nevertheless, a small decrease in kinetic energy is still observed. The reduction in
${\textit{Nu}}$
highlighted in figure 12(a) cannot be attributed solely to the Coriolis force, instead, it reflects a change in the dominant transport mechanism, where angular momentum transport by shear flow increasingly opposes the DEP force at larger
$\Delta \varOmega$
. Further discussion is provided in § 4.
Panels (a) and (b) show the Nusselt number,
${\textit{Nu}}$
, and kinetic energy density,
$\tilde {E}_k$
, respectively, as a functions of
$\Delta {\varOmega }$
for three values of
${\textit{Ra}}_{{E}}={1.59\times 10^4}(\times )$
,
${\textit{Ra}}_{{E}}={3.18\times 10^5}$
(
$\square$
) and
${\textit{Ra}}_{{E}}={8.84\times 10^5}$
(
$\circ$
).
4. Discussion
Using the same classification introduced in the previous section, a comprehensive regime diagram is presented in the
$\Delta {\varOmega }$
–
${\textit{Ra}}_{{E}}$
space based on forcing intensity, as shown in figure 14. This diagram includes all the investigated cases. The onsets of the DEP force-dominated regime and non-isothermal sTC flow-dominated regime are confined between the transition lines where the minimum in
${\textit{Nu}}$
is also indicated. For small values of the
${\textit{Ra}}_{{E}}$
, the flow remains time invariant, and axis symmetric, even at large
$\Delta {\varOmega }$
, up to moderate
${\textit{Ra}}_{{E}}$
. At moderate values of both
${\textit{Ra}}_{{E}}$
and
$\Delta {\varOmega }$
, periodic flow states emerge, while irregular flow patterns develop at larger
${\textit{Ra}}_{{E}}$
reminiscent of spherical Rayleigh–Bénard convection.
Regime diagram in the
$\Delta {\varOmega }$
–
${\textit{Ra}}_{{E}}$
parameter space, where time-invariant, periodic and irregular flow states are denoted by (
), (
) and (
), respectively. The line (
) marks the onset of the transition from the DEP force-dominated regime (A), passing through the transitional regime (B), characterised by a drop in Nusselt number,
${\textit{Nu}}$
, to the non-isothermal sTC flow-dominated regime (C) found above line (
). The minimum in
${\textit{Nu}}$
defines the transition line (
).
To provide an indication of the decrease in
${\textit{Nu}}$
in the transitional regime, we recap the nature of both basic flows. The nature of sTC flow was analytically described by Haberman (Reference Haberman1962), Ovseenko (Reference Ovseenko1963) and Munson & Joseph (Reference Munson and Joseph1971). For a given co-rotating spherical shell system, their analyses show that as long as nonlinear effects remain limited, the flow consistently exhibits a radially outward motion at the equator, a meridional circulation cell and a polar radial inward-directed flow. Beyond the stability threshold, however, the flow can develop significantly more complex structures, as demonstrated by simulations in Wicht (Reference Wicht2014). The nature of TEHD convection in spherical systems is consistently dominated by radially oriented plumes or sheets, as demonstrated in previous studies by Futterer et al. (Reference Futterer, Hollerbach and Egbers2008, Reference Futterer, Egbers, Dahley, Koch and Jehring2010). Beyond the boundary layer, the Nusselt number in the fluid bulk can be characterised by the inflow Nusselt number,
${\textit{Nu}}^q$
, as defined in (2.19) where
${\textit{Nu}}=Nu^q$
at the boundary. As demonstrated in the appendix, the definition of the inflow Nusselt number remains constant across all radii. The decrease in the Nusselt number at the wall, as observed in the previous section, can be uniquely attributed to changes in the convective and conductive heat transfer in the bulk. This decrease in
${\textit{Nu}}$
is analysed using the time-averaged convective term,
$\langle r^2u_rT\rangle _t$
, and the conductive term,
$\langle -r^2\partial _rT\rangle _t$
, for a representative case with
${\textit{Ra}}_{{E}} = {3.18\times 10^{5}}$
, as shown in figure 12(a). The changes in both terms in the bulk illustrate the alteration in flow characteristics and explain the transition to the transitional regime (B), with parameters listed in table 2.
Representative selected forcing parameters to indicate the decrease of
${\textit{Nu}}$
in the transitional regime for constant
${\textit{Ra}}_{{E}}$
.
Figure 15 presents the time-averaged components in an azimuthal pole-to-pole cross-section. For small
$\Delta {\varOmega }$
, as shown in figure 15(a) (left), the convective component is distributed throughout the gap, with pronounced intensities near the poles and equator. The heat flux in figure 15(a) (right) indicates a relatively uniform conductive distribution along both boundaries. This suggests that DEP forcing dominates along the meridian under small rotation, corresponding to weak angular momentum transport. The conductive heat transported from the boundary to the bulk is evenly distributed over the meridian, as shown in the scope view. As
$\Delta {\varOmega }$
increases to the point where
${\textit{Nu}}$
reaches its minimum, see figure 15(b), the convective term (left) in the fluid bulk becomes more pronounced near the poles and equator. This pattern marks a transitional regime where angular momentum transport begins to dominate over TEHD convection. The conductive heat flux (right) at the boundary shows enhanced heat transport across the outer shell near the equator and poles, while transport at the inner shell near the equator decreases. In the scope view, at the outer boundary there is nearly no heat transport towards the boundary possible. As the flow travelling alongside the wall to the pole nearly reached the temperature of the boundary and the heat capacity is saturated, consequently,
${\textit{Nu}}$
declines as the dominant transport mechanism shifts. As
$\Delta {\varOmega }$
increases further,
${\textit{Nu}}$
recovers. The convective component, shown in figure 15(c) (left), nearly vanishes in the fluid bulk, but remains strongly concentrated in the equatorial region near the inner spherical shell. A less pronounced vertical structure suggests the presence of a tangent cylinder. This pattern indicates significant angular momentum transport from the inner towards the outer shell, particularly near the equator. The conductive component (right) mirrors the convective dynamics at the outer shell’s boundary, with intense heat flux observed at both the equator and the pole on the outer boundary, interrupted by two regions of weaker heat flux. The inner shell shows a more pronounced conductive heat flux in the pole region for
$\theta ={\pi }/{2}\pm {\pi }/{2}$
.
Azimuthal pole to pole slice of the isosurface of the time-averaged convective term,
$\langle r^2u_rT\rangle _t$
(left) and the conductive term,
$\langle -r^2\partial _rT\rangle _t$
(right) for constant
${\textit{Ra}}_{{E}}={3.18\times 10^{5}}$
for different rotations in: (a)
$\Delta {\varOmega }={1.43\times 10^{1}}$
, (b):
$\Delta {\varOmega }={3.90\times 10^{1}}$
and (c):
$\Delta {\varOmega }={1.04\times 10^{2}}$
.
5. Conclusion
This study examined the interaction of convection in a non-isothermal sTC flow under the influence of the DEP force. To investigate how the DEP force interacts with differential rotation, we first analysed the two basic flow states independently: TEHD convection in a spherical shell, and non-isothermal sTC flow. For the TEHD convection case, we derived scaling laws describing the overall heat transport, characterised by
${\textit{Nu}}$
as a function of
${\textit{Ra}}_{{E}}$
and
$\tilde {E}_K$
. These results were compared with established scaling frameworks for convection (Grossmann & Lohse Reference Grossmann and Lohse2000).
A similar methodology was applied to the non-isothermal sTC flow. Following the approach of Wicht (Reference Wicht2014), we classified observed regimes into axisymmetric time-invariant and non-axisymmetric time-dependent flow. For the axisymmetric regime with
${\textit{Ro}} \approxeq 1$
, scaling laws were developed for the
${\textit{Nu}}$
and
$\tilde {E}_K$
, and compared with convection results at low
$Ra$
(Davis Reference Davis1922a
).
When both driving forces were combined,
${\textit{Nu}}$
generally increased. However, the radial DEP force can counteract angular momentum transport, particularly near the equatorial region. Within the explored parameter space, this interaction led to the identification of a transitional regime (C) characterised by reduced heat transport. Two additional regimes were identified: regime (A), dominated by DEP force-driven TEHD convection, and regime (B), dominated by rotation-driven flow with suppressed TEHD convection. To support this classification, components of a derived inflow Nusselt number after the principle of Eckhardt, Grossmann & Lohse (Reference Eckhardt, Grossmann and Lohse2007) and Wang et al. (Reference Wang, Jiang, Liu, Zhu and Sun2022) are used to count for conductive and convective contributions, highlighting the distinct transport mechanisms in each regime.
These findings advance the understanding of coupled TEHD convection and rotational convection in spherical geometries, with implications for thermal management, geophysical fluid dynamics and the design of EHD systems in spherical or rotating configurations. This study has focused on the parameter space of the AtmoFlow spherical shell experiment. Future work could extend the parameter space by exploring a broader range of
${\textit{Pr}}$
,
$\varGamma$
and nonlinear, time-dependent regimes.
Acknowledgements
The authors want to thank U. Harlander, F. Feudel and F.-T. Schön for useful comments on the article.
Funding
The project AtmoFlow is supported by the BMWi via the German Space Administration (Deutsches Zentrum für Luft- und Raumfahrt) under Grants Nos. 50WP1709, 50WM1841, 50WM2141 and 50WM2441 and via the National High Performance Computing centre NHR with Grant No. bbi00021.
Declaration of interests
The authors report no known conflicts of interest that could have impacted the work given in this article.
Author contributions
Y.G.: conceptualisation, methodology, software implementation, formal analysis, validation, investigation, resources, data creation, writing original draft preparation, visualisation, funding acquisition.
P.S.B.S.: conceptualisation, methodology, validation, formal analysis, investigation, resources, data creation, review and editing, project administration, funding acquisition.
C.E.: project administration, funding acquisition.
All authors have read and agreed to the published version of the manuscript. All gave final approval for publication and agreed to be held accountable for the work performed therein.
Appendix A. Mesh independence study
Table 3 summarises the quantitative results of the mesh-independence study described in § 2.5, performed for the selected parameter largest forcing parameter expressed by the Rayleigh number,
${\textit{Ra}}_{{E}}={1.77\times 10^{6}}$
, with
${\varOmega }=399.1$
,
$\Delta {\varOmega }=456.1$
and
$\gamma _e={6\times 10^{-2}}$
.
Mesh-independence test showing the convergence of the integral values of the Nusselt number,
${\textit{Nu}}$
, and the kinetic energy density,
$\tilde {E}_K$
, as functions of the total number of grid cells.
Appendix B. Derivation of the in-flow Nusselt number
An in-flow Nusselt number,
${\textit{Nu}}^q$
, can be derived by the approach of Wang et al. (Reference Wang, Jiang, Liu, Zhu and Sun2022) using the analogy of for the angular momentum transport in cylindrical Taylor–Couette flow Eckhardt et al. (Reference Eckhardt, Grossmann and Lohse2007).
For clarity we introduced the spherical coordinate quantities for velocity given by
$u_r$
,
$u_\theta$
,
$u_\varphi$
\begin{align} u=u_r, v=u_\theta , w=u_\varphi ,\end{align}
and define a spherical surface which is time averaged and proves the commutativity for
$r$
\begin{equation} \langle \ldots \rangle _{A,t}=\int \frac {r^2 \sin (\theta )}{4\pi r^2}\text{d}\theta \text{d}\varphi \ldots =\frac {1}{4\pi }\int \sin (\theta )\text{d}\theta \int \text{d}\varphi \ldots .\end{equation}
This surface average is valid for flows provided that a sufficiently long timespan is resolved for the averaging. All flows are driven by differential rotation, which, in these regimes, is equatorially and rotationally symmetric. Consequently, the meridional and azimuthal components cancel out. The non-dimensional temperature of (2.7) can now be expressed using spherical coordinates
\begin{align}&\qquad\qquad\qquad\qquad \partial _tT+u\partial _r T+\frac {v}{r}\partial _\theta T+\frac {w}{r\sin \theta }\partial _\varphi T={\nabla} ^2T ,\end{align}
\begin{align}& \partial _tT+u\partial _r T+\frac {v}{r}\partial _\theta T+\frac {w}{r\sin \theta }\partial _\varphi T= \frac {1}{r^2}\partial _r\big (r^2\partial _r T\big )\nonumber \\&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad +\frac {1}{r^2\sin \theta }\partial _\theta \left (\sin (\theta )\partial _\theta T\right )+\frac {1}{r^2\sin (\theta )^2}\partial _\varphi ^2T. \end{align}
By adding the continuity equation in spherical coordinates multiplied by the temperature, one arrives at
\begin{align} T\boldsymbol{\cdot }\left (\frac {1}{r^2}\partial _r\big(r^2u\big)+\frac {1}{r\sin \theta }\partial _\theta (v\sin \theta )+\frac {1}{r\sin \theta }\partial _\varphi w\right )&=0\boldsymbol{\cdot }T, \end{align}
\begin{align} \frac {2uT}{r}+T\partial _ru+\frac {vT\cot \theta }{r}+\frac {T}{r}\partial _\theta v+\frac {T}{r\sin \theta }\partial _\varphi w&=0, \end{align}
and after some transformation at the following equations:
\begin{align}& \partial _tT+\frac {1}{r^2}\left (\partial _r\big(r^2uT\big)-\partial _r\big (r^2\partial _r T\big ) \right )+ \frac {1}{r\sin \theta }\left (\partial _\theta (\sin (\theta )vT) - \frac {1}{r}\partial _\theta \left (\sin (\theta )\partial _\theta T\right )\right ) \nonumber \\&\quad+ \frac {1}{r\sin \theta }\left (\partial _\varphi wT-\frac {1}{r\sin (\theta )}\partial _\varphi ^2T \right )=0, \end{align}
\begin{align}&\qquad \partial _tT+\frac {1}{r}\big (\partial _r\big (r^2\left [ uT- \partial _r T\right ]\big )\big )+ \frac {1}{\sin \theta }\left (\partial _\theta \left (\sin \theta \left [ vT - \frac {1}{r}\partial _\theta T\right ]\right )\right )\nonumber\\&\qquad\quad+ \frac {1}{\sin \theta }\left (\partial _\varphi \left ( wT-\frac {1}{r\sin (\theta )}\partial _\varphi T \right )\right )=0 .\end{align}
Applying the time and space averaging method as defined in Wang et al. (Reference Wang, Jiang, Liu, Zhu and Sun2022) and expressed in (B2), this simplifies to
\begin{equation} \big \langle \partial _r\big (r^2\left [ uT- \partial _r T\right ]\big )\big \rangle _{A,t}=0 .\end{equation}
Since
$r$
is commutative under integration, the integration with respect to
$r$
can be performed from either side. We define the resulting heat flux, normalised by the conductive heat flux, as the inflow Nusselt number, denoted by
${\textit{Nu}}^q$
, as follows:
\begin{equation} Nu^q=\frac {\left (\eta -1\right )^2}{4\pi \eta }\big \langle (r^2\left [ uT- \partial _r T\right ]\big \rangle _{A,t} .\end{equation}
Open access
