Free Access
Issue
A&A
Volume 578, June 2015
Article Number A106
Number of page(s) 15
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201425125
Published online 12 June 2015

© ESO, 2015

1. Introduction

Our lack of understanding of (magneto)hydrodynamic transport processes in stars has hampered progress in developing the stellar evolution theory since its earliest beginnings. One particular aspect of the problem is the mixing in the boundary layers between convection and radiative zones in stellar interiors, which is also known as the problem of convective overshooting. Despite the indisputable advance in numerical simulations, the problem remains extremely challenging owing to the extreme range of the length and time scales involved in it.

The set of physical mechanisms that provide mixing at a convective/stable interface very likely depends on the type of convection zone involved. Because it is exposed to outer space at the top, a convective envelope is driven by the cold plumes originating in the photosphere. It is quite possible that the plumes span the whole convection zone and even provide mixing at its bottom boundary (cf. Andrássy & Spruit 2013, and references therein). A convective core or shell is, on the other hand, fully embedded in the star, its stratification is much weaker, and the temperature fluctuations within it are much smaller. Consequently, a different set of physical mechanisms may dominate mixing at its boundary.

It has long been known that the kinetic energy of the low-Mach-number flow in a convective core (or shell) is so low that the convective motions are stopped within about one per cent of the pressure scale height once they enter the steep entropy gradient of the radiative zone (Roxburgh 1965; Saslaw & Schwarzschild 1965). The motions can reach much farther, though, if they are vigorous enough to flatten the radiative entropy gradient above the convective core. In this case, known as the process convective penetration, the motions gradually “erode” the radiative stratification on the thermal time scale until radiative diffusion stops any further advance of the erosion front (Shaviv & Salpeter 1973; van Ballegooijen 1982; Zahn 1991). Finally, the fluid parcels hitting the stable stratification always generates a spectrum of internal gravity waves, which may also provide a certain amount of mixing (Press 1981; Garcia Lopez & Spruit 1991; Schatzman 1996).

Several of the processes mentioned may operate at the convection zone’s boundary at the same time. Their effects on long time scales and at long distances from the boundary are very different. In full numerical hydrodynamic simulations, the restrictions on time scales that can be covered makes it difficult to disentangle these effects. Physical insight developed by different means is needed to extrapolate them to longer time scales and distances.

We take a closer look at one specific process operating at a convective/stable interface in the interior of a star. Thermal diffusion causes temperature fluctuations from the convection zone’s boundary to spread into the stable stratification. Temperature differences on surfaces of constant pressure set up a flow even in the absence of momentum transport by hydrodynamic stress. We call this process “differential heating”, explore the physics of it in an idealised set-up, and estimate what amount of mixing it could cause in the stellar interior.

2. The differential heating problem

2.1. Problem formulation and simplification

Consider a horizontal, solid surface with a stably-stratified fluid overlying it1. A temperature fluctuation imposed at the surface propagates into the fluid by a diffusive process and upsets the hydrostatic equilibrium. We investigate what the properties of the resulting flow are.

By replacing the convective/stable interface by a solid wall, we eliminate all the phenomena related to the inertia of the convective flows and the shear induced by them. This allows us to study the physics of differential heating in isolation. The upper boundary is taken far enough not to influence the flow. Next we introduce further assumptions to facilitate the mathematical description and the subsequent analysis of the problem:

  • (1)

    The flow is confined to a layer that is significantly thinner than the pressure scale height.

  • (2)

    The fluid is a chemically homogeneous, ideal gas.

  • (3)

    The Brunt-Väisälä (buoyancy) frequency of the stratification is constant.

  • (4)

    Thermal diffusivity is constant.

  • (5)

    The gravitational field is homogeneous.

  • (6)

    The differentially heated surface is flat and horizontal.

  • (7)

    The flow is constrained to two spatial dimensions.

Assumption (1) allows us to use the Boussinesq approximation and turns out to be justified. The chemical homogeneity that we assume in (2) is, at least for the nuclear-burning layers in a star, only realistic at the onset of the burning. The differential heating process above a convective core would weaken as the nuclear burning progresses owing to the increase in the mean molecular weight in the core. We focus on the chemically homogeneous case to keep the number of parameters tractable. We introduce (3) and (4) for the same reason. The Brunt-Väisälä frequency depends on the distance from the convective/stable interface in a real star. The constant frequency in our analysis can be thought of as a typical value for the layer influenced by differential heating. Finally, we add the last three assumptions to make our analysis more transparent and to reduce the computational costs of the numerical solutions. We have to keep in mind, however, that the constraint (7) might influence the stability properties of the flow, and thus some of our conclusions may not apply to the three-dimensional case.

The Boussinesq equations are (Spiegel & Veronis 1960)

where u is the fluid velocity, D/Dt = /∂t + u·∇ the Lagrangian time derivative, ρm and Tm are the mean density and temperature, respectively, p and T the pressure and temperature perturbations, respectively, g is the gravitational acceleration, k a unit vector pointed in the vertical direction, ν the kinematic viscosity, N the Brunt-Väisälä frequency, w the vertical velocity component, and ϰ the thermal diffusivity.

Equations (1)–(3) still contain several dimensional parameters. It is crucial to realise that there is a natural system of units for the differential heating problem that makes the equations dimensionless. The flow in this problem is set off by thermal diffusion in a stably stratified medium, therefore the inverse of the Brunt-Väisälä frequency, 1 /N (or a multiple of it), is a natural unit of time. Having made this choice, we can define a natural unit of distance as , which is a typical thermal-diffusion length scale on the time scale 1 /N. The dimensionless Boussinesq equations are then

where we omit any symbol to indicate the new units. We have also introduced a new pressure-like variable p = p′/ρm and the buoyancy acceleration ϑ = gT′/Tm, which we continue to call the “temperature fluctuation” in the rest of the paper, because that is the central concept in the differential heating process. The Prandtl number Pr = ν now becomes a measure of kinematic viscosity, because the new unit of diffusivity is ϰ. Equations (4)(6) are particularly well suited to theoretical studies since their solution is fully determined by the Prandtl number, the initial, and the boundary conditions.

The distance unit is rather short in stellar interiors, and it only weakly depends on the stratification. To see this, we express the Brunt-Väisälä frequency in terms of the more common stellar-structure parameters, (7)where Hp is the pressure scale height, ad the adiabatic temperature gradient, and the actual temperature gradient. Close to a convection zone’s boundary, we can write (8)where α ≈ 10-1 is a coefficient of proportionality and z the distance from the boundary (z> 0 in the stable stratification). When using Eqs. (7) and (8), the unit of distance can be expressed as (9)which is about 107 cm for values typical of a point close to the convective/stable interface (z ≈ 10-2Hp) in the core of a massive (10 M), main-sequence star (ϰ≈ 1010 cm2 s-1, α ≈ 10-1, g ≈ 105 cm s-2, Hp ≈ 1010 cm).

Two distinct regimes of differential heating can be expected, depending on the amplitude and the spatial scale of the temperature fluctuation imposed on the differentially heated surface. If the heating is strong enough, the heat transport is advection-dominated (i.e. the flow’s Péclet number is high), and the flow is generally unsteady. A similar phenomenon takes place right at the point where the convective flow leaves the unstable stratification, still retaining some positive temperature fluctuation. It quickly cools down as it rises in the stable medium, its temperature fluctuation turns negative, and the flow is brought to a halt. This is the place where we can impose a lower boundary condition for a much weaker kind of differential-heating-induced flow, in which diffusive heat transport plays a major role (i.e. the flow’s Péclet number is low). The latter case is the main focus of this paper. We show in Sect. 3 that such a flow is generally smooth and reaches a stationary state (to be specified in Sect. 3.1) even at rather high values of the Reynolds number, up to Re = 4 × 103. This allows us to gain some insight into the problem by exploring the scaling properties of the stationary differential-heating equations, which we do in the next section.

2.2. Analytical considerations

The stationary differential-heating problem is described in two dimensions by the set of equations (cf. Eqs. (4)–(6))

where x and z are the horizontal and vertical coordinates, respectively, with the z axis pointed against the gravitational acceleration vector, u is the horizontal velocity component, and w the vertical one. In what follows, we show how the characteristic properties of the stationary flow depend on the typical amplitude Θ and the typical horizontal length scale L of the heating applied.

Assume that there is a well-defined vertical length scale H in the differential heating flow pattern. Let us denote the typical horizontal and vertical velocities by U and W, respectively, and the typical pressure fluctuation by P. We then introduce a new set of variables , , û, ŵ, , and , which all reach values of the order of unity close to the differentially heated surface, and Upon making these substitutions in Eq. (10), we obtain (20)which implies the approximate relation (21)The horizontal momentum equation (Eq. (11)) attains the form (22)where Eq. (21) has been used, so the equality is only approximate. Nonetheless, we can see that the viscous terms are of the order of 1/Rex ≡ Pr/(UL) and 1/Rez ≡ Pr/(WH), where Rex and Rez are Reynolds-like numbers associated with horizontal and vertical motions, respectively. We introduce this unusual notation to characterise the relative contributions of the two viscous terms in the case of LH. We focus on this limit because it turns out to be the relevant one in stellar interiors (see Sect. 5). From now on, we assume Rex ≫ 1 and Rez ≫ 1. Equation (22) shows that pressure fluctuations are of the order of U2 in this high-Reynolds-number limit, so that we can estimate (23)The vertical momentum equation (Eq. (12)), with the substitutions defined above and Eqs. (21) and (23), becomes (24)and implies a close balance between the vertical component of the pressure gradient and the buoyancy-acceleration term provided that LH in addition to Rex ≫ 1 and Rez ≫ 1. This allows us to estimate (25)which is a plain, order-of-magnitude equality of the characteristic kinetic and potential energies. Finally, the energy equation (Eq. (13)) becomes (26)The diffusion terms in Eq. (26) are of the order of 1/Pex ≡ 1/(UL) and 1/Pez ≡ 1/(WH), where Pex and Pez are Péclet-like numbers associated with horizontal and vertical motions, respectively. We introduce them for the very same reason as we did in the case of Rex and Rez. Making use of Eqs. (21) and (25), we can put Eq. (26) into the form (27)which can be greatly simplified in the double limit of Pex ≫ Pez and Pez ≪ 1. In that case, the two terms in the parentheses on the right-hand side have to closely balance one another, so that we can estimate (28)and Eq. (27) becomes linear, (29)Equation (29) is a special case of the energy equation in the low-Péclet-number approximation of Lignières (1999).

Using Eq. (28), we eliminate H from Eq. (25) to get an estimate of U,L) and, with Eq. (21), also an estimate of W,L). The resulting relations also enable us to express Rex, Rez, Pex, and Pez as functions of Θ, L, and Pr. This way we obtain

One might be tempted to estimate the time scale τ of flow acceleration towards the stationary state directly from the buoyancy acceleration Θ provided by the temperature fluctuation imposed on the bottom boundary. It is crucial to realise that, as Eq. (24) shows, the buoyancy acceleration is almost completely compensated for by the vertical component of the pressure gradient in the case LH. It is only their difference that contributes to the vertical acceleration. We can, however, consider the horizontal acceleration provided by the horizontal component of the pressure gradient and write U/τP/LU2/L (see Eq. (23)). Using Eq. (30) we obtain (36)Finally, we would like to point out that the characteristic thermal-diffusion length scale corresponding to the time scale τ is τ1/2 ≈ Θ− 2/7L3/7, which scales with Θ and L in quite a different way than H does (see Eq. (28)). This comes about because our estimates take the back reaction of the flow on the temperature distribution into account.

2.3. Numerical solutions

The order-of-magnitude estimates derived in the preceding section assume that the flow is stationary and that there is a well-defined vertical length scale in the flow pattern. We performed a series of time-dependent, numerical simulations of the differential heating problem to confirm these assumptions and to determine how the solutions depend on the Reynolds number and how they decrease with height.

We have developed a specialised code dedicated to the study of the differential-heating problem, because the problem places rather high demands on the numerical scheme. For instance, it has to tackle the highly diffusive nature of the flow and its high aspect ratio and resolve a wide dynamic range within a single simulation box. The code is of the finite-difference type, and it solves the differential-heating equations on a collocated grid using a variant of the MacCormack integration scheme. The Poisson equation for pressure, which can be derived from Eqs. (4) and (5) (or (37), see below), is solved by a spectral method. Heat-diffusion terms are treated implicitly, again by a spectral method. In what follows, we discuss a few selected issues related to the numerical solution of the differential-heating equations that need to be borne in mind when interpreting our results. The reader interested in the details of the numerical scheme is referred to Appendix A.

We use periodic boundaries in the horizontal direction and force the shear stress and the vertical velocity component to vanish at the lower and upper boundaries of the computational domain. One could also use non-slip boundaries, but these are hardly more akin to the physical reality that motivated this study in the first place, so we omit this case. We impose a temperature fluctuation in the form ϑ(x,0) = Θsin(πx/L) at the bottom boundary and force the temperature fluctuation to vanish at the upper boundary. The parameters Θ and L can be identified with the same symbols as introduced in Sect. 2.2.

The high thermal diffusivity in the differential-heating problem forces us to use long implicit time steps for the heat-diffusion terms, which might have an adverse effect on the accuracy of the results. To show that this is not the case, we re-computed the simulations sr03, sr30, and Re1024 (see Tables 1 and 2 and Sect. 3), decreasing the time step by a factor of ten. This brings about a change in the velocity field, which is of the order of 0.1% in the cases sr03 and sr30 and of the order of 1% in the case of Re1024 (measured well away from the field’s zeroes). The reason for this insensitivity to the time step is the low Péclet number of the flow. Lignières (1999) shows that in the low-Pe regime, the energy equation can be approximated by a Poisson equation for the temperature fluctuation with w as a source term (see also our Eq. (29)). We do not use this approximation to make our code more versatile; instead, we naturally obtain a close equilibrium between the terms 2ϑ and w in Eq. (6) when the Péclet number is low. This equilibrium is reached so quickly that details of the evolution of ϑ towards the equilibrium become irrelevant.

Table 1

Parameters of the series of simulations sampling a patch of the parameter space { Θ,L } at the constant value of Re = 2.6 × 102.

Table 2

Parameters of the series of simulations with Re increasing at the fixed values of Θ = 10-3 and L = 100.

It is a well-known fact that any numerical advection scheme either requires adding a so-called artificial-viscosity term to guarantee stability or it involves some viscous behaviour implicitly. In either case, the effective Reynolds number does not even come close to the astrophysical regime with current computing facilities. The artificial viscosity (be it explicit or implicit) thus exceeds the physical one by a wide margin, so it demands special attention.

Suppose we include an explicit viscous term as in Eq. (5) to model the artificial viscosity. Equations (32) and (33) show that for L ≈ 103 (equivalent to ~Hp in the astrophysical case mentioned in Sect. 2.1) we have Rez ≈ 10-4 Rex as a consequence of HL. Using equidistant grids with up to 103 grid points in each direction, we can achieve Rex ≈ 103. It follows that Rez ≲ 10-1 and the vertical momentum transport is dominated by the artificial-viscosity term. A value Rez ≫ 1 is, however, expected in stellar interiors. We use a simple workaround, replacing the viscous term Pr∇2u by the anisotropic form Prx2u/x2 + Prz2u/z2. The coefficients Prx and Prz are re-computed at each time step from the relations and , where hx and hz are the horizontal and vertical grid spacings, respectively, and Regrid is the Reynolds number on the grid scale. We performed a few numerical tests of the code on a convection problem to determine that the value Regrid = 4 is a conservative trade-off between the amount of viscous dissipation and the code’s stability, so we use this value in all the simulations presented here.

The anisotropic form of artificial viscosity enables us to reach Rex ≫ 1 and Rez ≫ 1 at the same time on a grid of reasonable size. We show in Sect. 3 that the solutions with constant values of Prx and Prz decay exponentially with height. The effective, local Reynolds number decreases in proportion to the flow speed, and the solutions quickly become dominated by the artificial viscosity. This would also happen in an (otherwise idealised) stellar interior at some point but that point would be much farther from the convection zone’s boundary. Therefore, we generalise the artificial-viscosity terms, and the momentum equation (Eq. (5)) in 2D becomes (37)where are Prx(z) and Prz(z) are proportional to eηz with η being an adjustable parameter. We set η = 0 when we are not interested in the precise vertical profiles and use η> 0 to suppress the viscous terms when examining how the solutions decrease with height. The latter case, η> 0, is a rather touchy problem because one has to increase η in a few steps, always using the (almost) stationary flow from a previous run as an initial condition for the next run. Overestimating the value of η can lead to a lack of viscous dissipation in some parts of the computational domain and a numerical instability ensues. Finally, the very goal that we want to achieve by this treatment, i.e. the flow dynamics’ being dominated by inertial terms at great heights, becomes an issue since such a flow evolves on the extremely long time scale corresponding to its low speed.

3. Results

3.1. The stationarity and structure of the flow

Our numerical investigation of the diffusion-dominated differential heating problem reveals that the flow reaches a stationary state at all values of the Reynolds number that we have been able to achieve (up to Re = 4 × 103). We use the rate of change of the quantity umax = max | u | (taken over the whole simulation box) as a convergence monitor and a basis of our criterion for deciding the flow’s stationarity. We show in Sect. 2.2 that the relevant dynamical time scale near the differentially heated boundary should be close to τ given by Eq. (36) (confirmed a posteriori, see below). We pronounce the flow stationary and stop the simulation once the condition (38)has been met at least for one time scale τ. A direct implementation of this condition would involve extrapolation from the time scale of one time step, Δt, to a much longer time scale τ, which would amplify the round-off noise by a factor of τ/ Δt ≫ 1. Instead, we approximate Eq. (38) by (39)where is the Euler-time-stepped solution of the equation . Thus, is a smoothed version of umax, lagging behind it approximately by τ in time.

The flow in all of the simulations presented in this paper is composed of several layers of overturning cells with flow speed rapidly decreasing from one layer to the next (see Figs. 4 or 6). We characterise the flow properties close to the differentially heated surface by a vertical length scale H, defined as the height above the hottest spot at which the flow first turns over (i.e. w(L/ 2,H) = 0), and the typical horizontal and vertical velocity components and , respectively, where the maxima are taken over the whole simulation box. The symbols H, U, and W can be identified with the same symbols as used in Sect. 2.2. The flow always reaches its maximal horizontal speed at the bottom boundary and the maximal vertical speed in the first overturning cell above the hot spot. The flow pattern is asymmetric, with the maximum downward flow speed (reached above the cold spot), max( − w), always lower than the maximal upward flow speed, max(w). We define the characteristic numbers Rex, Rez, Pex, and Pez in an analogous way to what is used in Sect. 2.2 with the difference that now we have two Prandtl-like numbers Prx and Prz instead of one Prandtl number Pr.

thumbnail Fig. 1

Dependence of the global flow characteristics on the heating amplitude Θ and length scale L. Circles show the values derived from numerical simulations (Table 1). Solid lines show the scaling relations (Eqs. (28), (30), (31), and (36)), normalised to fit all but the four simulations at L = 101, which are expected to deviate from the scaling relations.

Open with DEXTER

3.2. Scaling relations

We computed a grid of 16 simulations to verify the analytical relations derived in Sect. 2.2. All of these simulations, summarised in Table 1, have a resolution of 256 × 512, and the vertical grid spacing was adjusted so as to obtain Rex = Rez ≡ Re = 2.6 × 102. The decision to fix the value of Re is motivated by the fact that the flow pattern turns out to be scalable over a large part of the parameter space provided that Re = const. In other words, while changing the heating parameters Θ and L at Re = const. does change the amplitude and the vertical scale of the flow, the structure of the flow, as seen in a system of normalised coordinates x/L and z/H, remains unchanged (see Fig. 4). We show in Fig. 1 our numerical results as compared with the scaling relations fitted to all but the four data points at L = 101. The excluded data points do not comply well with the premise LH and are thus expected not to follow the scaling relations. Allowing only the constants of proportionality to change in the fitting process, we obtain The unexpectedly good fit is a result of the flow’s scalability.

thumbnail Fig. 2

Dependence of the maximal horizontal velocity on the Reynolds number. Simulation data (circles) are connected by the solid line to guide the eye. The scaling law umax ∝ Re0.054 is shown by the dashed line for comparison.

Open with DEXTER

The constants of proportionality in Eqs. (40)–(43), as well as the structure of the flow, depend on Re. We illustrate this dependence by computing a series of simulations with resolution ranging from 322 to 40962. This way, we cover about two orders of magnitude in the Reynolds number from 3 × 101 to 4 × 103 (again with Rex = Rez ≡ Re). We perform this experiment at L = const. and Θ = const., so any change in Re reflects a change in Prx and Prz. Nevertheless, we present the dependence on Re, because the scalability of the flow shows that the absolute values of Prx and Prz do not matter. Ideally, we should choose the heating parameters so as to have Θ ≪ 1 and L ≫ 1 as the scaling relations hold true in this limit (see Sect. 2.2).

Equation (43) shows, however, that the flow’s dynamical time scale becomes extremely long in the same limit, thus making any high-resolution computation unfeasible. Therefore we use Θ = 10-3 and L = 100, which still keeps the energy equation approximately linear (since Pex ≈ Pex ≈ 10-2), but we forgo having LH. Nevertheless, we expect the changes in the flow with increasing Re in this case to be similar to those that would be seen in a simulation with LH because of the energy equation’s being linear in both cases. All of this series of simulations, summarised in Table 2, reach the stationary state as defined by Eq. (39). Figure 2 shows that the maximum horizontal velocity in the computational domain slowly increases in proportion to Re0.054 in the high-Re regime. The flow also becomes increasingly asymmetric, as shown by the ratio of the maximum upward and downward flow speeds plotted as a function of Re in Fig. 3. The seemingly asymptotic trend changes at the highest Reynolds number considered, but we do not know the reason for this change.

3.3. Flow at great heights

The flow speed in all our simulations quickly decreases with height. Figure 5 compares the vertical profiles of the root-mean-square (rms; computed in the x direction) vertical velocity component, wrms(z), in four simulations with widely disparate heating parameters (sr00, sr03, sr30, and sr33). We find that wrms decreases approximately as eβwz/H in a global sense with βw ≐ 1.5 almost independently of Θ and L. We use the values H(Θ,L) given by Eq. (40) instead of those measured in the simulations to normalise the z coordinate, because this brings the slopes much closer to one another. We have to keep in mind, though, that these flows are reasonably close to a stationary state only up to z/H ≐ 2.5, because our convergence criterion (Eq. (39)) is ignorant of the weak flow in the upper part of the simulation box, and consequently, that part of the flow is still slowly evolving when the computation is stopped.

thumbnail Fig. 3

Dependence on the Reynolds number of the flow’s asymmetry, characterised by the ratio of the maximum upward and downward flow speeds.

Open with DEXTER

thumbnail Fig. 4

Comparison of the flow structure in two simulations with very disparate heating parameters. The left panel shows simulation sr30 (Θ = 10-3, L = 101). The right panel shows simulation sr03 (Θ = 100, L = 104). In both cases, the vertical velocity component w, normalised to its maximal absolute value, is plotted on a split logarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a non-linear way to aid visualisation. The spatial coordinates are normalised by the characteristic length scales defined in Sects. 2.3 and 3.

Open with DEXTER

The simulations discussed so far use constant artificial-viscosity parameters Prx and Prz, which leads to a rapid decrease in the local Reynolds number with height (in proportion to the decreasing flow speed). We computed another two simulations, this time with Θ = 10-4 and L = 101. In the first one, we set Prx = const. and Prz = const. (the constant-Pr case hereinafter), just as we have done so far. In the other one, we set Prx ∝ eηz and Prz ∝ eηz as described in Sect. 2.3 to keep a local version of the Reynolds number approximately constant (the constant-Re case hereinafter). We increased the slope η from 0 in a few steps in order to make the ratio of the rms advection terms to the rms viscous terms, i.e. the local Reynolds number, as independent of height as possible; η = 2 turns out to be a good compromise in this case. There is a large-scale, residual variation by about a factor of four in the local Reynolds number, because the simple exponential profile of the artificial viscosity is not flexible enough to compensate for it. Using Fig. 2 we estimate that this variation can change the velocities by ~0.1 dex at most. Since our usual stopping condition, Eq. (39), cannot “sense” the weak flow at great heights, we judge the stationarity of the flow by comparing the rms values of the /∂t terms to the rms values of all other terms that appear in Eq. (37) and require the former to be significantly smaller than the latter. This way, we obtain the results summarised in Figs. 68. The constant-Pr flow can be considered stationary over the whole range shown, whereas the constant-Re flow is only stationary for z/H ≲ 3.8, because the topmost part of that flow evolves so slowly that a global oscillation develops before it has reached equilibrium (see Sect. 3.4 for details).

Figure 7 illustrates that the flow is somewhat faster at z/H> 1 in the constant-Re case, as could be expected from the massive increase in the local Reynolds number by as much as two orders of magnitude at z ≐ 3. Much more interesting is, however, that the overturning cells in the constant-Re case become apparently thinner with increasing height, hence with decreasing local temperature fluctuation. This observation suggests that the scaling relations derived in Sect. 2.2 could be used locally (see the dependence of H on Θ in Eq. (28)). Another piece of evidence for this hypothesis is shown in Fig. 8, in which we compare the relative rates of decrease in ϑrms(z), urms(z), and wrms(z). The envelope of ϑrms(z) can be approximated well by the function ϑe(z) ∝ eβϑz/H with βϑ = 1.7 for z/H ≲ 3. We then regard ϑe(z) as an estimate of the local temperature fluctuation and rewrite the scaling relations for the velocity components, Eqs. (30) and (31), to obtain their local versions, where ue(z) and we(z) are expected to be good envelope models of urms(z) and wrms(z).

thumbnail Fig. 5

Decline with height of the relative rms vertical velocity in four simulations with very different heating parameters. The flows are reasonably close to a stationary state only up to z/H ≐ 2.5.

Open with DEXTER

In other words, we expect ue(z) ∝ eβuz/H and we(z) ∝ eβwz/H with and . Indeed, these scalings turn out to be correct, as shown in Fig. 8. Similarly, we can produce a local version of Eq. (28), (46)where h(z) is a local, height-dependent estimate of a vertical length scale analogous to H. As a result, we expect h(z) ∝ eβhz/H with , i.e. a slow thinning of the overturning cells with increasing height, similar to what we observe in Figs. 68. This seemingly innocuous phenomenon has very grave consequences for the flow at great heights. Instead of fading out exponentially, it decreases even faster (see Fig. 8). We expand on this in Sect. 4.1 and derive a better model for the flow’s decline with height to show that the flow speed drops dramatically above a certain point.

thumbnail Fig. 6

Effect of two different artificial-viscosity prescriptions on the flow structure. The constant-Pr flow (Prx, Prz = const.; left panel) is compared with the constant-Re flow (Prx, Prz ∝ e− 2z; right panel). The vertical velocity component w is in both cases plotted on a split-logarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a non-linear way to aid visualisation.

Open with DEXTER

3.4. Late-time evolution of the flow

Having continued some of our simulations for as much as 104τ, we discover an intriguing phenomenon. At first, a horizontal mean shear flow develops on top of the differential-heating flow. Its amplitude grows, and the shear flow begins to oscillate at some point. Finally, the oscillation saturates at an amplitude ranging from ~10-3 to ~100 of the differential heating flow’s amplitude, depending on the parameters of the simulation. The oscillation’s period and development time strongly decrease with increasing Reynolds number. They do not seem to have an upper limit but approach 10τ at Re ≈ 103. This phenomenon most likely has a physical origin because decreasing the time step by a factor of ten does not affect the shear flow or its behaviour significantly. Any detailed study of this phenomenon is certainly beyond the scope of this paper, but our preliminary research suggests that it is unlikely to be a cumulative effect induced by internal gravity waves since it (1) also occurs in very small computational boxes, in which all internal-wave modes are over-damped by radiative diffusion; and (2) the temporal spectra of the average horizontal velocity are featureless at periods significantly shorter than that of the shear flow oscillation.

thumbnail Fig. 7

Effect of two different artificial-viscosity prescriptions on the rms vertical velocities. Note that the constant-Pr flow can be considered stationary over the whole range shown whereas the constant-Re flow is stationary only up to z/H ≐ 3.8.

Open with DEXTER

thumbnail Fig. 8

Comparison of the rms velocities and temperature fluctuations in the constant-Re case with two models approximating their global behaviour. Solid lines show urms (top), wrms (middle), and ϑrms (bottom). Dashed lines show the model, in which ue(z) ∝ eβuz/H, we(z) ∝ eβwz/H, and ϑe(z) ∝ eβϑz/H with , , and βϑ = 1.7. Dotted lines show the improved model given by Eqs. (52)(54) with γ = 1.3. The coefficients of proportionality have been adjusted for each variable independently.

Open with DEXTER

4. Interpretation of the results

4.1. Improving the model at great heights

After picking up the threads of Sect. 3.3, we presently find that the diffusion-dominated, high-Re differential-heating flow actually decreases faster than exponentially with height. To see this, we make use of two results from Sect. 3.3. First, that the scaling relations derived in Sect. 2.2 have their local analogues, which hold within the flow (compare Eqs. (28), (30), and (31) with Eqs. (46), (44), and (45), respectively). Second, that the envelope of ϑrms(z) can be approximated well by the function ϑe(z) ∝ eβϑz/H at low heights, where βϑ is independent of Θ and L. This allows us to write (47)The characteristic vertical scale H is linked to the heating amplitude Θ by Eq. (28) and is thus relevant close to the differentially heated surface, where the typical temperature fluctuations are of the order of Θ. A straightforward generalisation of Eq. (47) is obtained by replacing H by the local, height-dependent estimate h(z) given by Eq. (46). Upon doing so, we have (48)where we have introduced the new variables , z′ = z/H, and h′(z) = h(z) /H. By Eqs. (28) and (46) we have (49)and Eq. (48) becomes (50)where β may differ slightly from βϑ, because we have used an order-of-magnitude relation in the last step. Equation (50) shows that decreases with a fairly constant slope over a few orders of magnitude, but the slope starts to change as soon as a wider dynamic range is considered. Since the slope is proportional to , Eq. (50) describes a runaway process. Indeed, the solution is (51)and vanishes at a finite height of . The constant must be very close to unity as and we can simplify Eq. (51) to obtain (52)where we have returned to the non-primed variables and introduced a new constant , which is a parameter to be adjusted to fit the numerical data. Using the local scaling relations, Eqs. (44)(46), we derive the functional dependencies

The functions ϑe(z), ue(z), and we(z) are shown in Fig. 8. The constants of proportionality in Eqs. (52)–(54) have been adjusted independently, but all three functions share the value γ = 1.3. The good fit indicates that our line of reasoning is probably correct.

Can we conclude that the flow stops at the finite height we have just derived? No, since the scaling relations only work in the high-Re regime. Provided that Re is high close to z = 0, the flow speed quickly decreases according to Eqs. (53) and (54) until Re ≈ 1 is achieved at some height z1 < 7H/γ. The weak flow above this point is supported by viscosity and gradually vanishes as z → ∞.

4.2. Allowing for a buoyancy-frequency gradient

So far, we have assumed that the flow occurs in a particularly simple type of thermal stratification – one characterised by a typical buoyancy frequency Ntyp = const. Nevertheless, we aim to apply our results to the immediate vicinity of a convection zone, i.e. to a medium, in that the buoyancy frequency rises continuously from zero to a finite value. In this section, we first show how to estimate the value of Ntyp in such a setting and then reapply the techniques developed in Sect. 4.1 to demonstrate how the varying buoyancy frequency affects the global flow field.

To do this, we have to recover the dependence of all the relevant flow properties on Ntyp by returning to a system of physical units. We recall that we use 1 /Ntyp as a unit of time and (ϰ/Ntyp)1/2 as a unit of distance, which implies that the unit of velocity is Ntyp)1/2 and the unit of acceleration (hence of ϑ) is . We use these conversion factors throughout this section without mentioning them further. The height of the bottommost overturning cell is by Eq. (40) (56)where we have introduced the index “ph” to indicate the use of physical units for quantities that are dimensionless in the rest of our analysis. The buoyancy frequency N is now an increasing function of zph and can be approximated by Eqs. (7) and (8), (57)The overall scale of the flow pattern is given by the bottommost overturning cell, which is thus the most important. Therefore we estimate Ntyp = N(Hph/ 2), i.e. (58)and combine Eqs. (56) and (58) to obtain (59)where we have also expanded Θph = gΔT/Tm to emphasise the dependence on the imposed temperature fluctuation ΔT/Tm. We use the sign in Eq. (59) and also in Eqs. (60)(62) below to indicate that we do not expect these estimates to be off by more than a few tens of percent. The dependence of Hph on the heating amplitude and length scale is somewhat weaker in Eq. (59) compared with Eq. (56), because Eq. (59) takes into account that any gain in the flow’s vertical extent brings about an increase in the typical buoyancy frequency, which in turn makes further penetration harder. This effect can also be seen when we express the characteristic velocity components and the flow’s dynamical time scale in physical units, where the exponents have slightly changed compared with Eqs. (41)–(43).

The spatial variation of N brings on a first-order effect, too; that is, the stratification offers less resistance to overturning in the bottom part of the flow field compared with the rest of it. We mimic this effect by using the flow’s excellent scaling properties under the assumption that the flow behaves locally as if N was constant. Our goal is to improve upon the envelope models of ϑrms(z), urms(z), and wrms(z) derived in Sect. 4.1 by taking the dependence of N on height into account.

Our starting point is Eq. (48) with the difference that now we define , z′ = zph/Hph and h′ = hph/Hph. We caution the reader that ϑe,ph refers to a model with N = N(z) and not to a direct translation of ϑe that appears in Eq. (48) to physical units. The local vertical length scale of the flow, h(z) given by Eq. (46), can be translated to physical units directly, (63)This equation, together with Eq. (56), implies (64)where N′ = N/Ntyp = (2z′)1/2 (see Eqs. (57) and (58)). It is evident that h diverges for N′ → 0+, i.e. z′ → 0+. This effect is purely artificial because the divergence occurs within the bottommost overturning cell of the flow, and the large-scale model we are developing here cannot capture such local phenomena. We ignore the divergence for now because only h′ − 1 appears in Eq. (48) and use the same procedure as in Sect. 4.1 to derive a generalised version of Eq. (50), (65)where the parameter β has absorbed all coefficients of the order of unity. Its value should still be of the order of unity, but it may be different in this model compared with the model developed in Sect. 4.1. By analogy to the derivation in Sect. 4.1, we can write the solution to Eq. (65) in the form (66)where we have also returned to the non-primed quantities, and is a parameter of the order of unity. The typical velocity components and the typical vertical vertical length scale can be estimated using the local scaling relations, Eqs. (44)(46), and (66). We obtain where an explicit dependence on N appears after the transition to physical units. These expressions diverge for z → 0+ where N → 0+ (see Eq. (57)), which is just another illustration of the envelope models’ inability to capture local phenomena (see also the discussion above). The bottommost part of the flow should in reality behave approximately as if it was in a medium with N = Ntyp = const., so we can cut off the problematic part of the N(z) profile and use, for example, the function (70)instead of N(z) in practical calculations. Doing so makes the right-hand sides of Eqs. (67)(69) converge to unity as zph → 0+.

Just as the results of Sect. 4.1 do not mean that the flow vanishes at a finite height, neither the results of this section mean that. Again, the sudden drop in the typical velocities predicted by Eqs. (67) and (68) only signifies that the flow undergoes a transition to the low-Re regime at a relatively low height. Equations (67) and (68) cease to be usable from that point on and the weak flow supported by viscosity gradually vanishes as zph → ∞.

5. Application to stellar conditions

The flow in a layer of thickness hph(zph) and vertical velocity we,ph(zph) at distance zph from the boundary overturns a passive tracer in it on a time scale τm = hph/we,ph. This suggests an effective diffusion coefficient Deffhphwe,ph. For the first layer above the boundary, this is (71)At distance zph, Eqs. (68) and (69) give (72)where Ntyp is given by Eq. (58) and we have replaced N(zph) in Eqs. (68) and (69) by given by Eq. (70) as discussed in Sect. 4.2. The constant γ is of the order of unity but cannot be constrained further by our present analysis. It determines the maximum height zmax,ph that the mixing process can reach, zmax,ph = (9 /γ)7/9Hph.

For a specific example, consider the boundary of the core convection zone in a 10 M zero age main sequence star. This environment is characterised by α = d(∇ad − ∇)/d(zph/Hp) = 0.14, a thermal diffusivity ϰ= 5.9 × 1010 cm2 s-1, a gravitational acceleration g = 1.1 × 105 cm s-2 and a pressure scale height Hp = 2.9 × 1010 cm. A mixing-length estimate for convection in the core produces temperature fluctuations ΔT/Tm ≈ 10-6 on a horizontal length scale LphHp. Equation (59) then predicts that the typical vertical length scale is Hph ≈ 2 × 108 cm = 7 × 10-3Hp. The typical vertical velocity (Eq. (61)) is Wph ≈ 5 × 101 cm s-1. These numbers imply Pez = (WphHph)/ϰ ≈ 2; i.e., the bottom part of the flow is located right at the transition between the regions of advection-dominated and diffusion-dominated heat transport. This is not a coincidence, because we are modelling the region where heat leaks from the convective eddies, allowing them to turn over and sink back to the convection zone. Such a flow has to have Pez ≈ 1. Therefore, the effective diffusivity close to the convection zone, Deff(0) in Eq. (72), is of the same order as the diffusivity of heat ϰ. Diffusivities that are several orders of magnitude smaller than ϰ can be important on the long nuclear time scale. The maximum height reached by the differential heating process on this time scale can thus be approximated by zmax,ph. Assuming γ = 1 we obtain zmax,ph ≈ 4 × 10-2Hp.

Equation (72) is likely to be somewhat of an overestimate of the actual mixing rate of the differential-heating process. The layers mix on the hydrodynamic time scale in their interiors, but as long as they are stationary, transport of the tracer between layers takes place by diffusion. As in the case of semiconvective layering (cf. Spruit 2013), this reduces the effective mixing rate to the geometric mean of the microscopic diffusion coefficient κt of the tracer and the estimate (72).

More significantly, the picture is complicated by the time dependence of the convective heat source. For the 10 M example, only the bottommost part of the flow can approach the stationary flow speed before the heating pattern changes because the dynamical time scale τph ≈ 5 × 106 s (Eq. (62) with the parameter values stated above) is of the same order as the convective overturning time scale in the core. This is likely to lead to some form of averaging, reducing the effective amplitude of the source. The level of this effect can probably be investigated with a time-dependent simulation.

6. Summary

Various observations show that there is a need for some additional mixing at the interfaces between the convective and radiative layers of stars. Even processes that are too weak to be detectable in numerical hydrodynamic simulations need to be considered as candidate sources of this mixing, because the nuclear time scale on the main sequence is so much longer than the dynamical time scale of convection, and cumulative effects are likely to play an important role.

In this work, we have investigated one such weak process, which we call “differential heating”. The differential heating process occurs when radiative diffusion transports a temperature fluctuation from the boundary of a convection zone into the neighbouring stable stratification. The resulting perturbation of hydrostatic equilibrium triggers a weak flow, which may provide mixing up to some distance from the convection zone. We investigated the flow that is driven by a static temperature fluctuation varying sinusoidally along the solid horizontal boundary of a stably stratified, thin layer of gas. This low-Péclet number problem (i.e. a slow flow dominated by thermal diffusion) turns out to be intrinsically non-linear, in the sense that the horizontal structure of the flow is asymmetric. Even for symmetric boundary conditions, the upflow is narrower than the downflowing part for the flow, and the shape of the flow pattern is nearly independent of the amplitude of the driving temperature perturbation.

A few additional assumptions (Sect. 2) allow us to describe the problem by a set of dimensionless equations, the solution to which depends (apart from the boundary and initial conditions) only on the Prandtl number. We analysed these differential-heating equations for their scaling properties under the assumption that the flow is stationary (Sect. 2.2). An astrophysically interesting corner of the parameter space is characterised by Rex ≫ 1, Rez ≫ 1, Pex ≫ Pez, and Pez ≪ 1. (The x and z directions have to be distinguished because such flow has a high aspect ratio.) In this limit we derive a set of simple relations (Eqs. (28) and (30)–(36)) to describe how the global flow properties depend on the heating amplitude Θ and length scale L. We find, in particular, that the characteristic vertical length scale H depends only weakly on the heating parameters (Eq. (28)).

We developed a dedicated numerical code to solve the equations. The main difficulties are related to the highly diffusive nature of the flow, its high aspect ratio, and the need to resolve a wide dynamic range in the flow amplitude within the computational box (as much as five orders of magnitude). The flow in our two-dimensional, time-dependent simulations reaches a stationary state at all values of the Reynolds number that we have been able to achieve (up to Re ≡ Rex = Rez = 4 × 103). The flow is always composed of several layers of overturning cells, the shape of which depends only on the Reynolds number and not on the heating length scale L and amplitude Θ. This property makes the flow scaleable in the sense that the flow field corresponding to some heating parameters L1, Θ1 can be stretched in space and scaled in amplitude to get a good approximation of the flow field corresponding to a different set of heating parameters L2, Θ2 provided that Re is in both cases the same. This is also the reason the scaling relations derived in Eq. (2.2) fit the simulation data remarkably well at Re = const. (see Fig. 1). Increasing the Reynolds number has little influence on the flow speed, but it makes the flow pattern increasingly asymmetric.

We decrease the artificial-viscosity coefficients in the code with height in order to keep the Reynolds number approximately the same in every layer of flow cells. The numerical data show that the global scaling relations derived in Sect. 2.2 have their local analogues, which can be used within the flow. The flow speed’s decrease with height, being locally exponential, steepens with the decreasing flow amplitude according to the local scaling relations. Based on this we derive a model of the flow’s dependence on height, which closely fits the numerical data over the whole dynamic range that we have been able to cover (as much as five orders of magnitude, see Fig. 8). The model shows that the flow speed drops abruptly to a negligible value at a finite height. The local scaling relations also allow us to generalise our results to the more realistic case, in which the buoyancy frequency N increases with height (see Sect. 4.1).

We illustrated the typical scales associated with the stellar differential-heating process with the example of the convective core of a 10 M zero-age main sequence star (see Sect. 5). We approximate the mixing due to the differential-heating flow by an “effective” diffusion coefficient Deff, which is of the order of the diffusivity of heat near the convection zone and decreases with height according to Eq. (72). The mixing relevant for stellar evolution extends about 4% of the pressure scale height above the convection zone.

The main findings of the paper are:

  • (1)

    The flow has a cellular structure and reachesa stationary state at all values of the Reynoldsnumber that we have been able to achieve (up toRe = 4 × 103).

  • (2)

    Both global and local properties of the flow can be described by a set of simple analytical relations.

  • (3)

    The flow speed drops abruptly to a negligible value at a finite height above the source of heating.

  • (4)

    The mixing relevant for stellar evolution extends about 4% of the pressure scale height above the convection zone of a 10 M zero-age main sequence star.


1

Equivalently, the stably stratified fluid could be placed under the differentially heated surface. The role of hot and cold spots on the surface would be reversed in this case. We discuss only one case for the sake of concreteness.

Acknowledgments

We would like to thank Achim Weiss for the 10 M model, Ewald Müller and Maxime Viallet for enlightening discussions on numerical hydrodynamics, and the anonymous referee for critical comments that improved the overall presentation of the text.

References

Online material

Appendix A: Numerical methods

Appendix A.1: Integration scheme

We have adapted the standard MacCormack method to suit our specific problem. In the simplest case of a one-dimensional vector q of conserved quantities being advected on and equidistant grid with a spacing of Δx, MacCormack’s method can be written as where is the value of q at the kth grid point and the nth time step, Δt the time step, f(q) the flux function, and we use the convention that any parenthesised upper index refers to a sub-step of the method instead of a time-step index. The method is linearly stable provided that the CFL condition Δt ≤ Δx/ρ(A) is met, where A is the Jacobian matrix of the flux vector and ρ(A) is the largest characteristic value of A. Non-linear stability typically requires the addition of some form of artificial viscosity. MacCormack’s method is second-order accurate both in space and time.

We discretise Eqs. (4), (6), and (37) on a collocated, two-dimensional grid of M × N cells with constant cell spacing x, Δz). The two spatial dimensions and the presence of source terms in the equations forces us to significantly extend the basic MacCormack scheme. We begin by advecting the vector of variables q = (u,w,ϑ) in both spatial directions using Strang splitting, where we have written out the explicit form of the flux terms. The indices k and l refer to the position along the x and z axes, respectively. We proceed by adding the source terms to the momentum equations, where we use second-order-accurate central differences to keep up with the order of accuracy of the advection scheme, νl = Prx(zl) and μl = Prz(zl) are the coefficients of our anisotropic artificial-viscosity prescription (see Sect. 2.3), and μl + 1/2 = (μl + μl + 1)/2. The new velocity field is, in general, slightly divergent. We correct for this divergence by subtracting the gradient of a pressure-correction field, u(1) = u(1d) − Δtp)(1). The condition ·u(1) = 0 leads to a Poisson equation for the pressure correction, (A.9)Since we use central differences to compute partial derivatives, the discrete form of the Laplace operator in Eq. (A.9) should be derived by applying the central differences twice. That would, however, lead to a sparse operator and cause odd-even-decoupling problems on our collocated grid. Therefore we use the standard compact Laplacian and solve the approximate pressure-correction equation (A.10)Equation (A.10) is solved by a spectral solver, see Sect. A.3. Having computed the pressure correction, we apply it to the velocity field, The approximate nature of the pressure-correction equation (Eq. (A.10)) causes ·u(1) to be small, but non-zero. Practical experience has shown that the residual divergence is negligibly small in the flows analysed in this paper provided that the boundary conditions are treated properly, see Sect. A.2. We should also write at this point, but our numerical tests have shown that the residual divergence in the velocity field becomes much smaller if we set , so we use the latter form. The next step is to integrate the remaining two terms in the energy equation. We begin by adding the w term, (A.13)where its latest available value, w(1), has been used. The diffusion sub-step is given by the implicit equation (A.14)which is also solved by a spectral solver, see Sect. A.3. We have thus completed the first step of the MacCormack scheme, analogous to Eq. (A.1), and obtained the new variables u(1), w(1), p(1), and ϑ(1). The second step, which we do not do not go into detail on, differs from the first one at two points. First, advection is done using backward-space flux differencing, as in Eq. (A.2) (compare with Eq. (A.1)). Second, the pressure field is updated in this step, i.e. . The final step of the MacCormack’s scheme, Eq. (A.3), is used in the same form, with q = (u,w,ϑ). We also update the pressure field in the same way, , so that we obtain an estimate of the pressure field for the next time step.

thumbnail Fig. A.1

Effect of three different methods of treating the pressure at the solid top and bottom boundaries. In all three panels, the vertical velocity component, w, is plotted on a split logarithmic colour scale. We use Θ = 10-3, L = 101, constant kinematic viscosity and set the resolution to only 16 × 64 to make the spurious oscillations visible. We obtain the result plotted in the left panel using the simple symmetry conditions for pressure (Eqs. (A.19) and (A.20)). Preceding the pressure-gradient computation by third-order pressure extrapolation to the ghost cells reduces the oscillations’ amplitude by a factor of ~100 (middle panel). Increasing the extrapolation order to six brings about another decrease by a factor of ~30 in the oscillations’ amplitude (right panel). The pressure gradient is in all three cases computed by the second-order central differences in the whole computational domain.

Open with DEXTER

Finally, there is a simple way of increasing the accuracy of the scheme at a given grid resolution, which we use. The MacCormack method contains a built-in asymmetry: Eqs. (A.1) and (A.2) show that it always starts with forward-space flux differencing and continues with backward-space flux differencing. The two flux-differencing methods can be reversed, obtaining a “reverse” MacCormack method, without decreasing the order of accuracy of the overall scheme. We compute every time step using both the “direct” and the “reverse” methods and use the arithmetic average of the estimates given by the two methods.

Appendix A.2: Boundary conditions

The treatment of boundaries is restricted by our decision to use spectral solvers, which do not allow changing the differentiation operators anywhere in the computation domain. We use the ghost-cell technique for this reason. The boundary conditions we impose on the differential-heating flow are summarised in Sect. 2.3. The periodic boundaries in the horizontal direction are trivial to implement. The solid boundaries on the top and bottom of the computational domain, however, require much more care. We implement them using reflective boundary conditions for the velocity vector, so that the imaginary walls are located at l = −1/2 and at l = N − 1/2. The conditions imposed on u also eliminate any shear on the boundary. The pressure field is required to be symmetric with respect to the solid boundaries, The conditions imposed by Eqs. (A.15)–(A.20) can easily be shown to be consistent with the pressure-correction equation (Eq. (A.10); sum both sides over k = 0, 1,...,M and l = 0, 1,...,N). They typically do, however, bring about a cusp in the pressure field along the normal to the walls. The resulting discontinuity in the vertical pressure gradient then propagates to the rest of the domain and can be seen as a low-amplitude oscillatory field superimposed on the true pressure field (see the left panel of Fig. A.1). We tried to cure this problem by changing the discretisation of the vertical-gradient operator at the walls, so that the ghost cells would not be used when computing the pressure gradient. This solution has met with very little success, most likely because the abrupt change in the operator brings about an abrupt change in the discretisation error so the problem remains. Quite surprisingly, preceding the pressure-gradient computation by high-order pressure extrapolation to the ghost cells has turned out to be an effective solution, able to eliminate nearly all of the spurious oscillations (see the middle and right panels of Fig. A.1). We therefore use sixth-order extrapolation in the simulations with constant artificial viscosity and increase the extrapolation order to ten when we let the artificial viscosity decrease with height. This technique cannot be viewed, however, as an all-purpose solution, because it is likely to be too unstable to be useful when computing highly turbulent flows.

We require the temperature fluctuation ϑ to have a fixed sinusoidal profile at the bottom boundary and to vanish at the upper boundary, which translates into

Appendix A.3: Spectral solvers

We use spectral methods to solve the two equations involving the Laplace operator, the Poisson equation for the pressure-correction equation (Eq. (A.10)) and the implicit heat-diffusion equation (Eq. (A.14)). We express both the knowns and unknowns as linear combinations of the Laplacian’s eigenfunctions that comply with the desired boundary conditions. The solution procedure is then much simplified and effective, provided that the transform to the eigenfunction basis can be computed efficiently.

In case of the pressure-correction equation (Eq. (A.10)), we use the linear transform (A.23)and its inverse (A.24)to transform any field fk,l to an array of complex amplitudes and back. We can see that the basis functions in Eq. (A.24) are periodic in k and even around l = −1/2 and l = N − 1/2; i.e., they comply with our boundary conditions on the pressure field (see Sect. A.2). Upon using the spectral decomposition defined by Eq. (A.24) on both sides of the pressure-correction equation (Eq. (A.10)), we readily obtain its solution in the wavenumber space, (A.25)where we have omitted the upper indices because the expression applies to both steps of the MacCormack scheme, Ŝk,l is the transformed right-hand side of Eq. (A.10). The eigenvalues λm,n of the Laplacian are (A.26)and can be pre-computed. We set λ0,0 to a large number to prevent division by zero and make the undetermined component vanish.

In case of the heat-diffusion equation (Eq. (A.14)), we use the linear transform (A.27)\newpage

and its inverse (A.28)to transform any field gk,l to an array of complex amplitudes ĝm,n and back. We can see that the basis functions in Eq. (A.28) are periodic in k and odd around l = −1/2 and l = N − 1/2; i.e., they comply with our boundary conditions on the temperature field in case of a vanishing heating amplitude (see Sect. A.2). To allow for an arbitrary heating profile at the bottom boundary, we take out the known boundary term from the Laplacian on the right-hand side of Eq. (A.14) and treat it as a source term. One can show that it is the same as replacing the diffusion equation ∂ϑ/∂t = 2ϑ by the equivalent equation (ϑζ) /∂t = 2(ϑζ), where ζ is the static solution to the diffusion equation ∂ζ/∂t = 2ζ with the desired boundary conditions (ζ can be pre-computed for a fixed heating profile). The boundary conditions on the difference ϑζ are then identically zero, and the spectral decomposition defined by Eq. (A.28) can be used. This way we obtain an explicit expression for the solution of the implicit Eq. (A.14) in the wavenumber space, (A.29)where the eigenvalues Λm,n of the Laplacian are (A.30)and can be pre-computed. An equation analogous to Eq. (A.29) relates to .

In the practical implementation, we use the FFTW library (Frigo & Johnson 2005) to compute the transforms in Eqs. (A.23), (A.24), (A.27), and (A.28). We combine standard, one-dimensional transforms of different kinds to obtain the non-standard, two-dimensional transforms that we need. Namely, Eq. (A.23) is implemented as a series of DCT-II transforms over the rows of the input array, after which the columns of the resulting array are transformed by a series of DTF transforms. The backward transform (Eq. (A.24)) is then computed by a series of DFTs followed by a series of DCT-IIIs. The transforms for the diffusion equation (Eqs. (A.27) and (A.28)) are implemented in the same way, but simply replacing the DCT-IIs by DST-IIs and DCT-IIIs by DST-IIIs. The transforms from the FFTW library do not include the normalisation factor (2MN)-1.

All Tables

Table 1

Parameters of the series of simulations sampling a patch of the parameter space { Θ,L } at the constant value of Re = 2.6 × 102.

Table 2

Parameters of the series of simulations with Re increasing at the fixed values of Θ = 10-3 and L = 100.

All Figures

thumbnail Fig. 1

Dependence of the global flow characteristics on the heating amplitude Θ and length scale L. Circles show the values derived from numerical simulations (Table 1). Solid lines show the scaling relations (Eqs. (28), (30), (31), and (36)), normalised to fit all but the four simulations at L = 101, which are expected to deviate from the scaling relations.

Open with DEXTER
In the text
thumbnail Fig. 2

Dependence of the maximal horizontal velocity on the Reynolds number. Simulation data (circles) are connected by the solid line to guide the eye. The scaling law umax ∝ Re0.054 is shown by the dashed line for comparison.

Open with DEXTER
In the text
thumbnail Fig. 3

Dependence on the Reynolds number of the flow’s asymmetry, characterised by the ratio of the maximum upward and downward flow speeds.

Open with DEXTER
In the text
thumbnail Fig. 4

Comparison of the flow structure in two simulations with very disparate heating parameters. The left panel shows simulation sr30 (Θ = 10-3, L = 101). The right panel shows simulation sr03 (Θ = 100, L = 104). In both cases, the vertical velocity component w, normalised to its maximal absolute value, is plotted on a split logarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a non-linear way to aid visualisation. The spatial coordinates are normalised by the characteristic length scales defined in Sects. 2.3 and 3.

Open with DEXTER
In the text
thumbnail Fig. 5

Decline with height of the relative rms vertical velocity in four simulations with very different heating parameters. The flows are reasonably close to a stationary state only up to z/H ≐ 2.5.

Open with DEXTER
In the text
thumbnail Fig. 6

Effect of two different artificial-viscosity prescriptions on the flow structure. The constant-Pr flow (Prx, Prz = const.; left panel) is compared with the constant-Re flow (Prx, Prz ∝ e− 2z; right panel). The vertical velocity component w is in both cases plotted on a split-logarithmic colour scale. The length of the velocity vectors (arrows) is scaled in a non-linear way to aid visualisation.

Open with DEXTER
In the text
thumbnail Fig. 7

Effect of two different artificial-viscosity prescriptions on the rms vertical velocities. Note that the constant-Pr flow can be considered stationary over the whole range shown whereas the constant-Re flow is stationary only up to z/H ≐ 3.8.

Open with DEXTER
In the text
thumbnail Fig. 8

Comparison of the rms velocities and temperature fluctuations in the constant-Re case with two models approximating their global behaviour. Solid lines show urms (top), wrms (middle), and ϑrms (bottom). Dashed lines show the model, in which ue(z) ∝ eβuz/H, we(z) ∝ eβwz/H, and ϑe(z) ∝ eβϑz/H with , , and βϑ = 1.7. Dotted lines show the improved model given by Eqs. (52)(54) with γ = 1.3. The coefficients of proportionality have been adjusted for each variable independently.

Open with DEXTER
In the text
thumbnail Fig. A.1

Effect of three different methods of treating the pressure at the solid top and bottom boundaries. In all three panels, the vertical velocity component, w, is plotted on a split logarithmic colour scale. We use Θ = 10-3, L = 101, constant kinematic viscosity and set the resolution to only 16 × 64 to make the spurious oscillations visible. We obtain the result plotted in the left panel using the simple symmetry conditions for pressure (Eqs. (A.19) and (A.20)). Preceding the pressure-gradient computation by third-order pressure extrapolation to the ghost cells reduces the oscillations’ amplitude by a factor of ~100 (middle panel). Increasing the extrapolation order to six brings about another decrease by a factor of ~30 in the oscillations’ amplitude (right panel). The pressure gradient is in all three cases computed by the second-order central differences in the whole computational domain.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.