Internally heated and fully compressible convection: Flow morphology and scaling laws

In stars and planets natural processes heat convective flows in the bulk of a convective region rather than at hard boundaries. By characterizing how convective dynamics are determined by the strength of an internal heating source we can gain insight into the processes driving astrophysical convection. Internally heated convection has been studied extensively in incompressible fluids, but the effects of stratification and compressibility have not been examined in detail. In this work, we study fully compressible convection driven by a spatially uniform heating source in two-dimensional (2D) and 3D Cartesian, hydrodynamic simulations. We use a fixed temperature upper boundary condition which results in a system that is internally heated in the bulk and cooled at the top. We find that the flow speed, as measured by the Mach number, and turbulence, as measured by the Reynolds number, can be independently controlled by separately varying the characteristic temperature gradient from internal heating and the diffusivities. Two-dimensional simulations at a fixed Mach number (flow speed) demonstrate consistent power at low wave number as diffusivities are decreased. We observe convection where the velocity distribution is skewed towards cold, fast downflows and that the flow speed is related to the length scale and entropy gradient of the upper boundary where the downflows are driven. We additionally find a heat transport scaling law which is consistent with prior incompressible work.

Figure 1

Left panel: Symbols indicate the simulations conducted in this work in $\mathcal{H}\text{−}\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}$ parameter space. Two-dimensional simulations are indicated with a solid circle. Three-dimensional simulations are indicated with the square outline, $\times$, and $+$ symbols. The marker color indicates the value of Ra and $\mathcal{H}$ for each point. The left side of each 2D point indicates the value of Ra (upper color bar) and the right side indicates the value of $\mathcal{H}$ (right color bar). Three-dimensional points are colored by the value of Ra. While these colors are redundant with the position on the plot, these colors are used consistently throughout the rest of the figures to refer to $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}$ and $\mathcal{H}$. Right panel: Initial conditions for $T$ and $\rho$ where dashed lines correspond to an adiabatic polytrope stratification and solid lines show the “Fast-IC” initial conditions, where we approximate the temperature structure of the final state and solve a boundary value problem for the density profile which satisfies hydrostatic equilibrium and mass conservation, see Appendix pp2. The “Fast-IC” profile shown is for $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}=2\times {10}^5$ and $\mathcal{H}=10.667$.

Figure 2

We show entropy fluctuations with the adiabatic entropy subtracted. For each panel we set the colorbar limits at the 3rd and 97th percentile values of entropy from that simulation. These minimum and maximum values are displayed above each panel. Note the substantial imbalance between the extent of the minimum and maximum values. We show 2D dynamics in the left column and 3D dynamics in the right column. In panels (a), (b), and (c) we show 2D dynamics at $\mathcal{H}=0.75$ and with $\text{Ra}/{\text{Ra}}_{\text{crit}}=2\times {10}^3,2\times {10}^5,$ and $2\times {10}^7$, respectively. As Rayleigh number increases, the flows become more turbulent. In panel (d) we show a simulation with the same Rayleigh number as panel (c) but with stronger internal heating at $\mathcal{H}=16.5$. In panel (e) we display the 3D simulation with the same parameters as panel (a) and in panel (f) we display the 3D simulation with the same parameters as panel (b). An animation of this figure can be found in the Supplemental Materials.

Figure 3

Time-averaged probability distribution functions (PDFs) of entropy fluctuations with the adiabatic entropy subtracted (left panel) and vertical velocity (right panel) are shown for 2D (green line) and 3D (orange line) simulations. Both simulations use $\text{Ra}/{\text{Ra}}_{\text{crit}}=2\times {10}^5$ and $\mathcal{H}=0.75$. The PDFs are averaged over 12 heating timescales. The entropy PDF is skewed towards negative values and is nearly identical between the 2D and 3D case. In contrast, the velocity PDFs exhibit different behavior in 2D and 3D. The 3D case has a narrower velocity distribution with more pronounced skew towards downward flows.

Figure 4

Left panel: Time- and horizontally averaged vertical energy flux profiles for the 3D simulation with $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}=2\times {10}^4$ and $\mathcal{H}=10.667$ plotted against height. The imposed flux from the source term is shown with the solid black line and the adiabatic flux [defined in Eq. (25)] with the dashed black line. The time-averaged total flux $\overline{F_{\mathrm{total}}}$ is shown with the orange line which converges to imposed heat flux $\overline{F_{\mathrm{imposed}}}$. The green line shows the conductive flux $\overline{F_{\mathrm{cond}}}$ which is adiabatic at the lower boundary and interior. The conductive flux becomes large in the superadiabatic upper boundary layer. We show the convective flux $\overline{F_{\mathrm{conv}}}$ with the purple line. We also display the kinetic energy flux, a component of the convective flux, with a dashed pink line. Due to the asymmetry introduced by stratification the KE flux is negative. In the right panel we show the entropy gradient. The entropy gradient is zero everywhere except the superadiabatic upper boundary. The integral of the entropy gradient, denoted by the shaded gray region, is the change in entropy, $\mathrm{\Delta }s$, which is important in predicting the Mach number of the flows as we discuss in Sec. 3d.

Figure 5

Time evolution of Reynolds number (left panel) and Nusselt number (right panel). We display time traces for three 2D simulations at $\mathcal{H}=0.75$ at $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}=1.9\times {10}^3,1.9\times {10}^5$, and $1.9\times {10}^7$. The $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}=1.9\times {10}^7$ simulation uses Fast-IC, and the $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}=1.9\times {10}^3$ simulation uses an adiabatic polytrope for its initial condition. We show both adiabatic (orange line) and fast initial conditions (brown line) for the simulation at $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}=1.9\times {10}^5$ in the Reynolds-number plot. The simulations using Fast-IC reach a thermally converged state quickly. At high Ra, the instantaneous value of Nu becomes chaotic, so we also plot the rolling time average over 40 heating timescales which is denoted by the darker colored trace.

Figure 6

Upper left panel: The scaling of average Mach number against $\mathcal{H}$, where color represents $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}$. Dots represent 2D simulations and the $\times$, box, and $+$ symbols represent 3D simulations. We display a ${\mathcal{H}}^{1/2}$ power law with a dashed black line which approximates the scaling observed in the data. Upper right panel: The scaling of average Reynolds number against $\mathcal{H}$, where color represents $\mathrm{Ra}/{\mathrm{Ra}}_{\mathrm{crit}}$. Lower left panel: The scaling of Mach number with Rayleigh number where color denotes $\mathcal{H}$. Note that in 3D there is a downward trend in Mach number with increasing Rayleigh number. The power laws found by Ref. [42] (AB17) are shown by gray lines, with the solid line showing the 2D power law from AB17 and the dashed line showing the 3D power law. Lower right panel: The scaling of average Reynolds number against supercriticality where color denotes $\mathcal{H}$. We show a ${\mathrm{Ra}}^{1/2}$ power law with a dashed black line.

Figure 7

We show Nusselt number against Rayleigh number for our suite of simulations. We find a $\text{Nu}\propto {(\text{Ra}/{\text{Ra}}_{\text{crit}})}^{1/5}$ scaling, which is consistent with the incompressible internally heated simulations from Goluskin and Spiegel [34].

Figure 8

In the left panel we display a scatter plot of $\mathrm{\Delta }s$ against $\text{Ra}/{\text{Ra}}_{\text{crit}}$. Colors correspond to values of $\mathcal{H}$, and squares, $\times$, and $+$ symbols correspond to 3D runs. Circles denote 2D runs, with large circles indicating a 2D run with a corresponding 3D run at the same control parameters. In the middle panel we show a scatter plot of $\mathrm{\Delta }s$ against $\mathcal{H}$, with the same marker styling as in the left panel but with colors now corresponding to Rayleigh number. In these two panels we show that the scaling of $\mathrm{\Delta }s$ with our control parameters is consistent between 2D and 3D simulations and the small differences that are seen are not from a change in the power-law exponent. In the right panel we show how Mach number scales with $\mathrm{\Delta }s$. We find a scaling consistent with $\text{Ma}\propto {\left(\mathrm{\Delta }s\right)}^{1/2}$ for 3D simulations and for 2D simulations at a fixed value of Ra. However, when we increase Ra in 2D simulations $\mathrm{\Delta }s$ decreases. We show an inlay plot with just the simulations with $\mathcal{H}=0.75$. We see that for 3D, all simulations fall on the $\text{Ma}\propto {\left(\mathrm{\Delta }s\right)}^{1/2}$ trend line, but the 2D simulations do not, with increased Ra leading to smaller $\mathrm{\Delta }S$ and slightly decreased Ma.

Figure 9

Left panel: Convective power spectra of 2D simulations accumulated over 40 heating timescales. We plot the power spectra of the mid-$z$ plane velocity vector for a set of simulations with constant $\mathcal{H}$. We find consistent power at low wave number for these simulations. The less turbulent simulations enter the viscous regime at lower $k_x$ than the more turbulent simulations. We plot canonical $k^{-5/3}$ and $k^{-3}$ power laws and note that the $k^{-3}$ power law better matches the energy cascade in our simulations. The sharp vertical cutoff of the spectra is found at the highest possible wave number given the simulation's resolution. Right panel: Convective power spectra of 3D simulations accumulated over two heating timescales. We plot power spectra of the mid-$z$ plane velocity vector for a suite of simulations at constant $\mathcal{H}$. The low-wave-number power decreases with increasing Ra.

Figure 10

Left panel: Example of a thermal equilibrium solution with $\mathcal{H}=10.667$ to be used as input to the linear stability analysis. We plot temperature, $T$ (green), and density, $\rho$ (orange), against height, $z$. Middle panel: Example of growth rates with $\mathcal{H}=10.667$ plotted against horizontal wave number on the $x$ axis and Rayleigh number on the $y$ axis. Black dots indicate the interpolated value of Ra where the growth rate is zero for the given wave number. Right panel: Values of ${\text{Ra}}_{\text{crit}}$ plotted against $\mathcal{H}$. We find that ${\text{Ra}}_{\text{crit}}$ is roughly constant below $\mathcal{H}\sim 1$, and increases for $\mathcal{H}\gtrsim 1$.

Figure 11

An example of time traces of Reynolds number (left panel) and $\mathrm{\Delta }s$ (right panel) for simulations with $\mathrm{Ra}=2\times {10}^5{\mathrm{Ra}}_{\mathrm{crit}}$ and $\mathcal{H}=0.75$ with “Fast-IC” (orange) and adiabatic initial conditions(green). Gray horizontal lines show the time- and volume-averaged values reported in the Appendix pp3 [66]. Arrows show the time where a rolling average over 500 data points has converged to within $3\%$ of the the accepted value.

Figure 12

[(a) and (c)] An example of entropy fluctuations at $t=t_{{\mathrm{d}}_{\mathrm{fast}}}$ for simulations with $\mathrm{Ra}=2\times {10}^5{\mathrm{Ra}}_{\mathrm{crit}}$ and $\mathcal{H}=0.75$ with adiabatic initial conditions (a) and “Fast-IC” (c). The run with adiabatic initial conditions has reached convective onset; however, cold downflow plumes stall out part way through the atmosphere, becoming more buoyant (redder) than the nearby fluid. As time progresses convective downflows reach further into the atmosphere. In the run with “Fast-IC” the downflow plumes cross the domain immediatly after convective onset and the dynamics quickly reach a converged state. The colorbar is calculated from the adiabatic initial conditions case in the same manner as in Fig. 2 and is used in all dynamics panels. In panel (e) we show probability density functions (PDFs) for entropy fluctuations averaged over 12 heating timescales for both adiabatic-IC and Fast-IC simulations at $t=t_{{\mathrm{d}}_{\mathrm{fast}}}$. The adiabatic-IC run has not converged at this time as evidenced by the discrepancy betweeen the PDFs. In the right column [(b), (d), and (f)] we show dynamics and PDFs at $t=6000$, which is after dynamics and structure have converged for both choices of initial conditions. The dynamics of the adiabatic-IC simulation (b) are comparable to the dynamics of the Fast-IC run (d). In panel (f) we show the PDFs for both runs where the distributions have converged, indicating that the dynamics have reached a statistically similar final state.