Scaling in internally heated convection: a unifying theory

We offer a unifying theory for turbulent purely internally heated convection, generalizing the unifying theories of Grossmann and Lohse (2000, 2001) for Rayleigh--B\'enard turbulence and of Shishkina, Grossmann and Lohse (2016) for turbulent horizontal convection, which are both based on the splitting of the kinetic and thermal dissipation rates in respective boundary and bulk contributions. We obtain the mean temperature of the system and the Reynolds number (which are the response parameters) as function of the control parameters, namely the internal thermal driving strength (called, when nondimensionalized, the Rayleigh--Roberts number) and the Prandtl number. The results of the theory are consistent with our direct numerical simulations.

Thermally driven turbulence is omnipresent in nature and technology. The thermal driving can be thanks to the temperature boundary conditions such as in Rayleigh-Bénard convection -a flow in a container heated from below and cooled from above [1][2][3] -or in horizontal convection [4][5][6] or vertical convection [7][8][9][10], where parts of the top, bottom, or sidewalls of the container are set at different temperatures. However, the thermal driving can also be thanks to internal heating, where the temperature field is driven by some forcing in the bulk. In many cases in nature, both ways of driving play a role at the same time.
E.g., this holds for the Earth's mantle due to the driving through the hot inner core of the Earth and an additional driving due to the decay of radioactive material inside the core, producing heat [11][12][13][14][15][16][17][18]. Thus, in the Earth's mantle, about 10-20% of the heat is transferred from the core, while the rest occurs due to the internal heating [18]. The internal heating dominates also in the atmosphere of Venus [19,20], which is heated up due to the absorption of sun light. One more example is the formation of Pluto's polygonal terrain, which is caused not only by convection of Rayleigh-Bénard type [21,22], but also by internally heated convection [23]. And, of course, internally heated convection is relevant in many engineering applications, e.g., liquid-metal batteries [24,25].
To obtain a theoretical understanding of thermally driven turbulence including the cases of mixed thermal driving, it is mandatory to first understand the pure and well defined cases, namely on the one hand turbulent Rayleigh-Bénard convection (RBC, exclusively driven by the heated and cooled plates), and on the other hand turbulent purely internally heated convection (IHC, exclusively driven by thermal forcing of the temperature field in the bulk), see [26][27][28][29]. For the former case, Grossmann and Lohse (GL) have developed a unifying theory [30][31][32][33][34], with which the heat transfer and the degree of turbulence can quantitatively be described as function of the control parameters, in excellent agreement with the experimental and numerical data over a range of more than 7 orders of magnitude in the control parameters Ra and P r. Later this theory was also extended to horizontal convection [5] (HC) and double diffusive convection [35]. GL arguments were also applied to IHC, to estimate the bulk temperature for small and moderate P r [27]. A complete theory, however, does not yet exist for purely internally heated convection.
The objective of the present work is to apply the reasoning of GL's theory to the case of purely internally heated convection and to develop a unifying theory for this case. In addition, we perform direct numerical simulations (DNS) of turbulent purely internally heated convection over a large range of control parameters and compare the DNS results with the theoretical predictions. The DNS are conducted in two dimensions (2D), as (i) the theory is based on Prandtl's equations, which are also 2D in spirit, as (ii) 2D and 3D thermally driven turbulence show very close analogies with respect to the integral quantities, in particular for large Prandtl numbers P r ≥ 1 [36], and as (iii) otherwise, due to unavoidable limitations in available CPU time, we could explore only a much smaller portion of the parameter space.
In RBC, next to the geometric aspect ratio Γ of the sample (the ratio between lateral and vertical extensions), the control parameters of the system are the temperature difference between top and bottom wall (in dimensionless form, the Rayleigh number) and the ratio between kinematic viscosity ν and thermal diffusivity κ, namely the Prandtl number P r = ν/κ. The response of the system consists of the heat flux from bottom to top (in dimensionless form, the Nusselt number N u) and the degree of turbulence (in dimensionless form, the Reynolds number Re). In IHC, instead of the Rayleigh number, the dimensionless driving strength Rr of the temperature field takes the role of the second control parameter, next to P r. It is often called Rayleigh-Roberts number (and that is why we use the abbreviation Rr) and will be defined below. The main response parameter, next to Re, is the mean temperature which the bulk achieves thanks to the internal driving. This is related to the heat fluxes into the top and bottom plates; note that they are different from each other. So the objective of this paper is to explain how the mean temperature and the Reynolds number in turbulent IHC depend on Rr and P r, for large enough aspect ratio Γ of the sample.
The flow in IHC is confined between two parallel plates with distance L, with the gravitational acceleration g ≡ −ge z acting orthogonally to these plates. The underlying dynamical equations within the Boussinesq approximation are the compressibility condition ∂ i u i = 0, for the velocity field u(x, t), the kinematic pressure field p(x, t), and the reduced temperature Here T plate is the temperature of both top and bottom plates, β is the thermal expansion coefficient, δ ij the Kronecker delta and Ω the constant bulk driving of the temperature field, which in non-dimensional form is called Rayleigh-Roberts number The equations (1) and (2)  The main responses of the system can be expressed in terms of the mean temperature ∆ ≡ θ(x, t) V achieved in the system, where the average · V is over volume and time. The nondimensional form of this response parameter is The other main nondimensional response parameter is the Reynolds number Re = U L/ν, Obviously, due to the internal heating, the heat flux (or in dimensionless formQ(z) ≡ Q(z)/(ΩL)) in the system is not constant as in RBC, but depends on the height z. Here, · means average in time and in a plane of constant z.
However, a simple integration of eq. (2) over the full volume yields that the quantitỹ is constant for all z and equals Q 0 is thus a further dimensionless response parameter of the system. Eq. (7) implies that the dimensionless heat fluxQ is non-positive at z = 0. Applying eq. (6) at z = L gives the dimensionless flux at z = L, Relations (7) and (8) immediately show that 0 ≤Q 0 ≤ 1.
As it is well known [1,37], in RBC exact relations between the time and volume averaged thermal and kinetic dissipation rates, . This phase space for IHC is the analogous one to the one of RBC of [31] and the one of HC in [5].
the dimensionless control and response parameters N u, Ra, and P r can be obtained from multiplying the thermal advection equation with θ(x, t) and the Navier-Stokes equation with u i (x, t) and subsequent Gauss integration and time and space averaging. Here we apply the same procedure to Eqs. (2) and (1) and obtain As u ≥ 0 is non-negative by definition, we can now further restrain the magnitude of the dimensionless heat flux through the bottom plate: 0 ≤Q 0 ≤ 1/2. Just as the corresponding relations in RBC, also here, Eqs. (9) and (10) relate the averaged thermal and kinetic dissipation rates with the dimensionless control (Rr, P r) and response (∆,Q 0 ) parameters.
The key idea of the GL theory [30,31] is to split the kinetic and thermal dissipation rates  into contributions from the corresponding boundary layers (BLs) and bulks, and to apply the respective scaling relations for those (i.e., for u,BL , u,bulk , θ,BL and θ,bulk ), based on boundary layer theory and Kolmogorov's theory for fully developed turbulence in the bulk. The introduced scaling regimes I, II, III, and IV correspond to BL-BL, bulk-BL, BL-bulk, and bulk-bulk dominance in u and θ , respectively. Here one should also take into account mean thicknesses of the thermal BLs (λ θ ) and viscous BLs (λ u ). The cases λ θ < λ u (large P r) and λ θ > λ u (small P r) correspond to different scaling regimes, and therefore we assign the subscripts u and to regimes I, II, III and IV, which indicate the upper-P r and ower-P r cases, respectively. Equating u and θ to their estimated either bulk or BL contributions and employing the classical Prandtl scaling relations for the BL thicknesses λ θ and λ u [38], one in principle obtains eight theoretically possible scaling regimes. Regimes II u and III are rather small, because, e.g., in II u , it is expected that λ θ ≥ λ u due to the BL-dominance in θ , but on the other hand, λ θ ≤ λ u should hold due to the large P r. By similar arguments, regime III is also small.
The value of u,bulk is estimated as which is relevant in the u -bulk dominating regimes II , IV and IV u , while the value of θ,bulk is estimated as which is relevant in the θ -bulk dominating regime IV . For large P r (regimes III u and IV u ), the thermal BL is embedded into the kinetic one and therefore in eq. (11), the magnitude of the velocity of the flow, which carries the temperature in the bulk, is reduced from U to The kinetic dissipation rate in the BL is ∼ ν(U/λ u ) 2 . Hence, which is relevant in the u -BL dominating regimes I , I u and III u . As in [30,31], the factor λ u /L accounts for the volume fraction of the kinetic BL. With increasing P r, λ u saturates to ∼ L, so this factor becomes one (just as argued in [31]), which yields For small Ra or very large P r, this leads to special regimes I < ∞ and III ∞ on top of, respectively, I u and III u . In III ∞ , also θ,bulk scales differently to (11), namely as The thermal dissipation rate in the BL scales as ∼ κ(∆/λ θ ) 2 , which is relevant in the θ -BL dominating regimes I , I u , I < ∞ and II . This (again with the volume fraction factor) leads to In the limiting regimes I , II and I < ∞ , it holds λ u /λ θ ∼ P r 1/2 [41,42], while in regime I u it holds λ u /λ θ ∼ P r 1/3 , all just as in the classical Prandtl-Blasius-Pohlhausen theory [38].
Equating u and θ to their estimated bulk or BL contributions, we obtain the scalings of∆ and Re in IHC, which are summarized in Tab. I and sketched in Fig. 1.
As already mentioned above, the very same idea was already applied to horizontal convection [5]. Interestingly enough, even a formal analogy between IHC and HC exists, out of which we could have already derived the scaling relations of Table I and Fig. 1. The reason for this formal analogy is that the relations obtained for θ and u (see Eqs. (5) and (6) of [5]) formally resemble the corresponding relations (9) and (10) here. For the first equation this becomes particular obvious when writing θ = L 2 κ Ω 2∆ = κ L 2 ∆ 2∆−1 and for the second when realizing that ( 1 2 −Q 0 ) is only a dimensionless factor between 0 and 1/2. Then one sees immediately that the role of the control parameter Ra in HC is taken by that of the control parameter R in IHC and the role of the response parameter N u in HC is taken by that of the (inverse) response parameter∆ −1 in IHC. All derived scaling relations in the different limiting regimes of HC can directly be taken over. The corresponding values for IHC give the same results as obtained above and have already been shown in Table I and To check these predictions of the GL theory generalized to IHC, we have performed 2D DNS according to Eqs. (1) and (2) with the corresponding BCs. We chose an aspect ratio of Γ = 2 for the laterally periodic box. The numerics have been validated by making sure that the exact relations (10) are fulfilled. Simulations were performed using the second-order staggered finite difference code AFiD [43,44]. This code has already been extensively used to study RBC, see, e.g. [45,46].
The parameter combinations (Rr, P r) for which we performed simulations are shown in the parameter space of Fig. 2a. A typical snapshot of the temperature field together with the mean temperature profile for one parameter combination are displayed in Fig. 2b.
One can see the stably-stratified layer near the bottom plate. The interaction of the upper convection zone and the lower stably-stratified region leads to the so-called penetrative convection [47,48]. The mean temperature profile, which, as expected and typical for IHC, displays top-down asymmetry.
The results for the response parameters∆ and Re as functions of the control parameters Rr and P r are shown in Figs. 3 and 4. As can be seen, in general, there are no pure scaling laws over the simulated range, but smooth crossovers from one regime to the other, very similarly as in RBC [34], reflecting the key idea of the unifying theory by [30,31]. We first discuss the dependences for the dimensionless mean temperature∆, see Figs. 3a, b. As a function of P r (Fig. 3b), for all Rr the transition from∆ ∼ P r −1/10 of regime I l to the P r-independence of regime I < ∞ can clearly be seen. The more turbulent regimes IV u, are not yet realized, as the driving is not strong enough. This is also reflected in the Rr dependencẽ ∆ ∼ Rr −1/5 reflecting that of regimes I u, . No indication to a stronger dependence as typical for the more turbulent regimes IV u, can yet be seen. This is also seen in the dependences of the Reynolds number (Fig. 3 c, d): For small P r ≤ 1, it goes like Re ∼ Rr 2/5 as in regimes I u, . For large P r = 10 the results are consistent with Re ∼ Rr 1/2 as in regime I < ∞ . This scaling should go hand in hand with the scaling∆ ∼ Rr −1/4 for the dimensionless mean temperature, but as seen from Fig. 3a, those data are presumably better described bỹ ∆ ∼ Rr −1/5 . Finally, on the P r-dependence of Re: As seen from Fig. 3d, for all Rr the data show a transition from the Re ∼ P r −4/5 scaling of regimes I u, to the Re ∼ P r −1 scaling of regime I < ∞ , consistent with the corresponding transition for∆ in Fig. 3b.
All these results are consistent with our unifying theory, which however goes much beyond  One sees from Fig. 4a thatQ 0 only weakly depends on Rr in the present parameter range; this behaviour has also been found before in [28]. Fig. 4b illustrates that much less heat is transported outwards from the bottom plate with increasing P r. The smallQ 0 for large P r is due to the less efficient shear-driven mixing of the fluid near the bottom plate.
In conclusion, in the spirit of the prior unifying theories for RBC [30,31] and for HC [5], in this paper we have developed a unifying theory of IHC for the scaling of the mean temperature and the Reynolds number as functions of the control parameters Rr and P r.
The main result is visualized in Fig. 1. We have shown that the 2D DNS results are consistent with this theory, though the numerically accessible regimes are still dominated by the boundary layers, and not all predictions of the theory can already be tested at this point.
Also 3D DNS are still outstanding. We have furthermore pointed towards the formal analogy between IHC and HC and it will be illuminating to explore this analogy also numerically.