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 = Pr_x(z_l) and μ_l = Pr_z(z_l) 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⁽¹⁾ = u^(1d) − Δt∇(Δp)⁽¹⁾. The condition ∇·u⁽¹⁾ = 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⁽¹⁾ 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⁽¹⁾, 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⁽¹⁾, w⁽¹⁾, p⁽¹⁾, and ϑ⁽¹⁾. 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.

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 f_k,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 g_k,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 = ∇²ϑ by the equivalent equation ∂(ϑ − ζ) /∂t = ∇²(ϑ − ζ), where ζ is the static solution to the diffusion equation ∂ζ/∂t = ∇²ζ 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.

© ESO, 2015