1 Introduction

The buoyancy boundary layer above a heated inclined plate immersed in a stratified cold fluid was introduced by Prandtl^[1] as a model of mountain and valley winds (see Fig. 1). When the temperature difference between the hot plate and the cold fluid at the same height is a constant, the basic flow of the inclined buoyancy layer is a parallel flow. Such a boundary layer also appears above the inclined bottom of an ice-covered lake^[2].

When the inclination angle of the plate relative to the horizon is χ= 90°, the inclined buoyancy layer becomes a vertical buoyancy layer. Early linear stability analyses of the buoyancy layers were mainly concerned with two-dimensional normal modes^[3-5], although Squire's theorem does not apply here^[3]. Tao and Busse^[6] studied the three-dimensional stability of the inclined buoyancy layer, and found that the most unstable modes in some region of the parameter space are not the two-dimensional Tollmien-Schlichting (TS) waves but the oblique rolls.

The linear stability analysis only provides a sufficient condition of hydrodynamic instability, while the simplest necessary condition of instability is given by the energy stability analysis (the energy method). The energy stability of the vertical buoyancy layer was studied by Dudis and Davis^[7], but the energy stability of the inclined buoyancy layer is still unknown at present, and this is the first motivation of this paper.

When the Grashof number is less than the critical Grashof number of energy stability Gr_E, the energy of any disturbance decays monotonically, which precludes any transient disturbance energy growth, and the basic flow is therefore globally stable^[8]. In the Grashof number range between Gr_E and the critical Grashof number of linear stability Gr_L, the disturbance energy is possible to grow initially, and then the disturbance may decay eventually according to the linear stability or initiate a subcritical transition with the help of nonlinear effects. The transient disturbance energy growth in parallel shear flows can be as large as O(1 000) at subcritical Reynolds numbers^[9-11]. However, the maximum transient energy growth has never been evaluated in the inclined buoyancy layers, and this is the second motivation of the present study.

Both the critical Grashof number of the energy stability and the transient energy growth depend on the definition of the disturbance energy. The conventional energy method uses the kinetic energy of disturbance as the definition, while the generalized energy method defines a positive-definite Lyapunov functional as a generalized disturbance energy^[12-13], and may shrink the difference between Gr_E and Gr_L. In the case of the Rayleigh-Bénard convection, the energy stability and the linear stability theories share the same critical Rayleigh number^[14]. In this paper, we will use the generalized energy method to prove that the basic flow of the vertical buoyancy layer is globally stable to streamwise-independent disturbances.

In the following, we will introduce the mathematical model of the inclined buoyancy layer in Section 2, and extend the conventional energy stability analysis of the vertical buoyancy layer^[7] to the case of the inclined buoyancy layer in Section 3. The maximum transient growth at subcritical Grashof numbers in the inclined buoyancy layer will be calculated in Section 4. In Section 5, we will study the linear and nonlinear stabilities of the basic flow to streamwise vortices in the vertical buoyancy layer. The conclusions will be summarized in Section 6.

2 Mathematical model

A schematic of the inclined buoyancy boundary layer is shown in Fig. 1. The coordinates in the streamwise, wall-normal, and spanwise directions are x^*, y^*, and z^*, respectively. The asterisks denote dimensional variables. The coordinate in the −g direction is s^*, where g is the gravitational acceleration. The fluid in the far field is stably stratified with a linear temperature profile T _∞ ^* (s^*) = T _∞ ^* (0) + ζ^*s^*, where T _∞ ^* (0) is a reference temperature, and ζ^* > 0 is the magnitude of the temperature gradient in the stably stratified fluid. The temperature at the surface of the plate is assumed to be T_w ^*(s^*) = T _∞ ^* (s^*) + ∆T^*, where ∆T^* is the constant temperature difference between the hot plate and the cold fluid at the same height s^*.

The dimensional governing equations with the Oberbeck-Boussinesq approximation are^[6]

(1)

(2)

(3)

(4)

(5)

where U^* = (U^*, V^*, W^*) is the velocity of the fluid. T^* and P^* are the temperature and the pressure of the fluid, respectively. Parameters ρ₀, γ, υ, and κ are the reference density, the coefficient of thermal expansion, the kinematic viscosity, and the thermal diffusivity of the fluid, respectively.

Define the non-dimensional coordinates, time, velocity, temperature, and pressure as^[6]

where P _∞ ^*(s^*) = P _∞ ^* (0) + ρ₀g(0.5γζ^*s^*2 − s^*), and

The Grashof number and the Prandtl number are defined as

(6)

where χ is the inclination angle of the plate relative to the horizon.

Then, we have the non-dimensional governing equations and boundary conditions,

(7)

(8)

(9)

(10)

(11)

(12)

(13)

The basic flow is a parallel flow^[1],

(14)

The non-dimensional governing equations and boundary conditions of disturbance (u, v, w, T') = (U, V, W, T) − (U₀, 0, 0, T₀) are^[6]

(15)

(16)

(17)

(18)

(19)

(20)

(21)

where u, v, and w are the velocity components of disturbance corresponding to x, y, and z, respectively. T' and p are the temperature perturbation and the pressure perturbation, respectively.

3 Energy method

The disturbance energy is defined as

(22)

where λ is a positive parameter to be determined. As a result,

(23)

where

(u₁, u₂, u₃) ≡ (u, v, w), and (x₁, x₂, x₃) ≡ (x, y, z). The Einstein summation convention is used here. Only when λ=0.25PrGr² will the energy maximum problem be well posed, because the integral is not bounded by the dissipation terms ||∇u||²+(λ/Pr)||∇T'||^2[13]. Actually, if we set v = w = 0 and

where A is a positive parameter, then the disturbance is second-order continuously differentiable and divergence-free, and satisfies the boundary conditions. When A → + ∞ , we have

Therefore, the integral cannot be bounded by the dissipation terms.

Using λ=0.25PrGr² and introducing θ=0.5GrT', we obtain

(24)

and

(25)

As the coefficient of the third integral term in (23) vanishes, the inclination angle χ influences only through the definition of the Grashof number (6). If an inclined buoyancy layer and a vertical buoyancy layer have the same Grashof number (according to our definition), they will have the same for the same disturbance.

Define the critical Grashof number of energy stability as

(26)

The maximum is searched among all divergence-free disturbances satisfying the boundary conditions, and is evidently positive. The critical Grashof number of the energy stability is always finite, because is bounded by ^[15]. We also define the maximum growth rate of the disturbance energy as , and then we have

(27)

and

(28)

When Gr < Gr_E, we have , because is a strictly increasing function of Gr provided that Gr_E is finite. Therefore, according to (28), the disturbance energy will monotonically decay to zero, and the basic flow is asymptotically stable to any disturbance with a finite amplitude.

When calculating the maximum in (26), we obtain the same Euler-Lagrange equations as in the energy method of the vertical buoyancy layer^[7],

(29)

(30)

(31)

(32)

(33)

where μ and ψ(x, y, z) are the Lagrange multipliers. The minimum positive eigenvalue μ is just the critical Grashof number of the energy stability Gr_E.

In the following, we will introduce a discretization method different from the one used in Ref. [7] to solve the eigenvalue problem (29) -(33).

Expanding the disturbances in Fourier integrals

and introducing the wall-normal vorticity and , we have

(34)

(35)

(36)

(37)

(38)

(39)

where , and k = (α² +β²)^1/2. is introduced here to obtain two second-order equations instead of a fourth-order equation of . This D² method helps us to calculate eigenvalues more accurately^[16].

Assume

(40)

where the summation is over all integers m and n such that 0 ≤ n ≤ N, −n ≤ m ≤ n, and m+n is even. This discretization is an extension of the original method of descending exponentials, which was used in the linear stability analysis of the saline boundary layer near the surface of a porous medium^[17]. The original method of descending exponentials takes the form of , which is reminiscent of the general Dirichlet series in the number theory^[18] , where a_n and s are complex numbers, and λ_n > 0 are strictly increasing and tend to infinity.

Substituting (40) into (34) -(39) and neglecting terms that are O(exp(−(k +N + 1)y)) when y → + ∞ , the eigenvalue problem of ordinary differential equations (34) -(39) is reduced to a matrix eigenvalue problem, which is solved with the QZ function in MATLAB.

Given α and β, denoting the minimum positive eigenvalue of (34) -(39) by µ₁(α, β), we have

(41)

For Pr=0.72, the minimum in (41) is achieved at (α, β)=(0, 0.485 69) (see Fig. 2). Actually, the minimum in (41) is always achieved at α=0 within numerical accuracy for 0≤ Pr≤ 2. Therefore, the least stable modes of the energy stability are streamwise vortices. This is also the case for the energy stability of plane Poiseuille flow^[19-20] and plane Couette flow^[21]. When Pr > 20, μ₁ is almost independent of (β/α)^[7].

Define the spanwise wavenumber of the least stable mode as β_3D. Therefore, Gr_E = µ₁(0, β_3D). We also define the critical Grashof number of the energy stability for two-dimensional disturbances as^[7]

(42)

where α_2D is the streamwise wavenumber of the least stable two-dimensional disturbance in the energy method. In this section, we use N=15 in the discretization, which is accurate enough for the calculation of eigenvalues (see Table 1).

We calculate the critical wavenumbers of the least stable two-dimensional and three-dimensional modes in the energy method for Prandtl numbers between 0 and 100 (see Fig. 3(a)). The critical wavenumbers do not coincide with the result of Dudis and Davis^[7], because the critical wavenumbers listed in their paper were only accurate to 0.05, which was implied in their later work^[22]. The critical Grashof numbers they calculated nevertheless had only 5% relative errors compared with our results (see Fig. 3(b)). Not demanding better accuracy, they numerically integrated the eigenvalue equations only from y=0 to y = L_y = 8. To check accuracy of the present method, we use their method with L_y = 20 and identify the eigenvalues in Table 1 to be between 39.111 626 150 6 and 39.111 626 150 7, and between 54.383 008 825 and 54.383 008 826. Therefore, our algorithm is reliable, and our results are accurate enough. Note that the Grashof number in this paper is equivalent to 2R in Ref. [7].

We also calculate the critical Grashof numbers of the linear stability (Gr_L) for angles of inclination 30° and 90° with the Chebyshev-Tau method^[23] (see Fig. 3(b)). The concern of Dudis and Davis^[7] is confirmed, i.e., the critical Grashof numbers of the energy stability and the linear stability do not coincide (Gr_E < Gr_L). Actually, Dudis and Davis were not sure about it, because the only available linear stability result at that time concerned merely two-dimensional disturbances (β=0) ^[3]. On the contrary, the Gr_L curves plotted here concern all two-dimensional and three-dimensional disturbances (α²+β² > 0).

When Gr < Gr_E, the basic flow is stable to all disturbances with finite amplitudes; when Gr > Gr_L, there exists at least one infinitesimal disturbance increasing exponentially. When Gr_E < Gr < Gr_L, the disturbance energy is able to increase at the beginning for some initial disturbances, and then the linear stability theory implies that disturbances will eventually decay if nonlinear effects are negligible.

We use SIMSON^[24] to simulate the nonlinear evolution of disturbances in the inclined buoyancy layer for χ=30° and Pr=0.72. The discretization in the streamwise and spanwise directions uses a two-dimensional Fourier series expansion. Therefore, the solution is periodic in these directions. A Chebyshev series expansion is used in the wall-normal direction. The simulation is performed in a computational domain of size (L_x, L_y, L_z) = (4π, 30, 2π/β) using (64, 129, 64) spectral modes. The initial condition is the least stable normal mode in the energy method with (α, β)=(0, 0.485 69). Therefore, the disturbance energy will grow initially if Gr > Gr_E = 39.11. In this simulation, L_x can be arbitrary because the normal mode is independent of x. The evolution of the average disturbance energy E' = E/(L_xL_z) is shown in Fig. 4(a). When t=0, the growth rates of the average disturbance energy in the direct numerical simulation (DNS) are in accord with (see Fig. 4(b)), where the coefficients are calculated by numerically integrating (25), with the least stable normal mode substituted.

The contours of the wall-normal velocity v of the least stable normal mode in the energy method when Pr=0.72 are plotted in the zy-plane (see Fig. 5(a)). The profile is scaled so that the average disturbance energy E'=1. The contour levels are ± 0.1, ± 0.2, …, and the contours of negative values are plotted with dashed lines. The contours of v=0 are omitted for clarity.

We also plot the contours of v of the least stable two-dimensional normal mode in the energy method when Pr=0.72 in the xy-plane as shown in Fig. 5(b), and plot the contours of v of the least stable normal mode in the linear stability theory when χ=30° and Pr=0.72 in Fig. 5(c). Both profiles are scaled so that E'=1, and have maximal v at y ≈ 2.5. The profile of the two-dimensional least stable normal mode in the energy method, however, extends farther into the far field.

4 Transient energy growth

When Gr_E < Gr < Gr_L, the disturbance energy is able to grow initially, and then the disturbance may be sustained by the nonlinear effects or decay according to the linear stability theory. The transient energy growth exists even when the nonlinear effects are negligible, because the growth rate of the disturbance energy does not vary when the disturbance v is scaled to λv for any constant λ according to (24) and (25). In this section, we will calculate the maximum transient growth in the disturbance energy based on the linearized Navier-Stokes equation.

Formally, the solution to the linearized Navier-Stokes equation is v(t) = e^Ltv(0), where L is a linear operator depending on the streamwise and spanwise wavenumbers (α, β) of the disturbance. Define the energy growth function as^[10]

(43)

where E is the disturbance energy defined in (24). The maximum growth for all time is defined as^[10]

(44)

Define the maximum growth for all time and wavenumbers as^[11]

(45)

where (α_max, β_max) are the wavenumbers of the optimal disturbance.

To calculate the energy growth function, we first discretize the linearized Navier-Stokes equation with a Chebyshev polynomial expansion in the interval [0, 30] in the wall-normal direction and with a two-dimensional Fourier series expansion in the streamwise and spanwise directions. After that the linear operator L is approximated by a finite-dimensional matrix L', the eigenvalues and the eigenvectors of which are solved with the QZ function in MATLAB. In the end, the energy growth function is calculated with the singular value decomposition of the matrix e^L't. The detail of the procedure was described by Reddy and Henningson^[11].

We plot the energy growth function G(t) for (α, β)=(0.4, 0.4), χ=90°, Pr=0.72, and Gr=200 in Fig. 6(a). The Grashof number is between Gr_E = 39.1 and Gr_L = 202.35 for χ=90° and Pr=0.72. The maximum growth for all time is G_max = G(0.730) =11.19, i.e., the disturbance energy is able to grow to 11.19 times as large as the initial disturbance energy. Note that the points in the curve G (t)correspond to different initial disturbances. The contours of the wall-normal velocity v of the initial disturbance corresponding to G_max are plotted in Fig. 6(b). The profile is scaled so that the average disturbance energy E'=1. The contour levels are ± 0.1, ± 0.2, …, and the contours of negative values are plotted with dashed lines. The contours of v=0 are omitted for clarity.

The calculation of the maximum transient growth is based on the eigenvalue decomposition of the linearized Navier-Stokes equation. We examine the validity of this algorithm by performing DNS with SIMSON^[24] and comparing the evolution of the disturbance energy in the DNS with that predicted by E(t) = ||v(t)||²/2 = ||e^L'tv(0)||²/2. The simulation is performed in a computational domain of size (L_x, L_y, L_z) = (2π/α, 30, 2π/β) using (64, 513, 64) spectral modes, where (α, β)=(0.4, 0.4). The initial disturbances in the DNS are given by scaling the disturbance in Fig. 6(b) to E'=0.1 and E'=10.

When the initial average disturbance energy is E'=0.1, the evolution of the disturbance energy agrees well with that predicted by the eigenvalue decomposition of the linearized Navier-Stokes equation (see Fig. 7(a)). The xz-plane averaged disturbance velocity and temperature , which originate from the nonlinear effects, are negligible compared with the basic flow (14) (see Fig. 7(b)). The three-dimensional flow structure is still an oblique roll, which implies that the amplitudes of disturbances with (α, β)≠(0, 0), (0.4, 0.4) are relatively small (see Fig. 7(d)).

When the initial average disturbance energy is E'=10, although the three-dimensional flow structure is similar to that of E'(0) =0.1 (see Fig. 7(e)), the amplitudes of the xz-plane averaged disturbance velocity and temperature are comparable to those of the basic flow (see Fig. 7(c)). For example, the relative modification of the basic flow u_ave/U₀ is about 6% at y=π/4, where U₀ achieves its maximum, and is about 72% at y=2.6, where u_ave achieves its maximum. The xz-plane averaged flow differs from the basic flow, which is the base of the linear transient growth analysis. Therefore, the evolution of E' in the case E'(0) =10 is diverted from the linear theoretical prediction (see Fig. 7(a)).

We summarize the maximum growth for all time G_max and the corresponding t_max in the αβ-plane for χ=90°, Pr=0.72, and Gr=200 in Figs. 8(a) and 8(b). The maximum growth for all time and wavenumbers is S=11.25, which is achieved at (α_max, β_max)=(0.379, 0.449) and t_max = 0.708 2.

However, the optimal disturbances that have the maximum growth in the whole wavenumber space are not always three-dimensional oblique rolls. We calculate the maximum growth for all time and wavenumbers S(χ, Gr, Pr) in the χGr-plane when Pr=0.72 (see Fig. 9(a)), and find that the optimal disturbances are streamwise vortices when χ≤ 50° (see Figs. 9(b) and 9(c)). In Fig. 9(a), the upper and lower thick lines represent Gr_L and Gr_E, respectively. According to the definitions (43) and (44), G_max = + ∞ when Gr > Gr_L, and G_max = 1 when Gr ≤ Gr_E. Therefore, we only plot the contours of S when Gr_E ≤ Gr ≤ Gr_L. When Pr=0.72, the maximum growth in the disturbance energy over all angles of inclination and Gr_E ≤ Gr ≤ Gr_L is S=11.62, which is achieved at χ=90° and Gr = Gr_L = 202.35.

5 Weighted energy method

In the vertical buoyancy layer, the least stable modes of the energy stability are streamwise vortices, which was already known by Dudis and Davis^[7]. However, in the following, we will prove that the basic flow of the vertical buoyancy layer is linearly stable to all streamwise vortices no matter how large the Grashof number is.

Theorem 1 The basic flow of the vertical buoyancy layer is linearly stable to all streamwise vortices.

Substituting the temporal normal mode of linear stability into equations of streamwise-independent disturbance in the vertical buoyancy layer, neglecting nonlinear terms and eliminating pressure terms, we have

(46)

(47)

(48)

(49)

where . Multiplying (47) by the complex conjugate of and then integrating it over (0, + ∞ ), we have

Then, we have either λ < 0 or If λ < 0, the normal mode is damped. If , (46)–(49) become

Eliminating , we have

(50)

(51)

Multiplying (50) by the complex conjugate of and then integrating it over (0, + ∞ ), we have

is not always zero, otherwise, we have , and then (u, θ) = (0, 0) is not a nontrivial eigenmode. Because λ is one of the roots of a quadratic equation with positive coefficients, the real part of λ is negative, which means that the normal mode is damped. Thus, we have proven that the basic flow of the vertical buoyancy layer is linearly stable to all streamwise vortices.

Remark 1 This proof is inspired by the proof of "the solutions to the Squire equation are always damped"^[8].

Although streamwise vortices are always damped in the linear stability theory, the critical Grashof numbers of the energy stability (Gr_E) of streamwise vortices are finite for all Prandtl numbers in the vertical buoyancy layer (see Fig. 3(b)). When Gr > Gr_E, although streamwise vortices may grow transiently, they will decay eventually if the nonlinear effects are negligible. Because Gr_E depends on the definition of the disturbance energy (24), we can increase Gr_E by defining other positive-definite Lyapunov functionals as the generalized disturbance energy. With this generalized energy method, we will prove in the following that the basic flow of the vertical buoyancy layer is stable to finite-amplitude streamwise-independent disturbances (including the streamwise vortices in the linear stability analysis).

Theorem 2 The basic flow of the vertical buoyancy layer is stable to all streamwise-independent disturbances.

For a streamwise-independent disturbance (u, θ), we have

(52)

(53)

where , f=u, v, w, or θ. From (52), we know v, w → 0. When ||v|| is sufficiently small, the production terms on the right hand side of (53) are negligible, and then u and θ will finally tend to zero.

The preceding discussion inspires us to define the weighted disturbance energy as

(54)

where σ is a positive parameter that will be determined later. As a result,

(55)

Notice that

(56)

where the maxima are searched among all streamwise-independent divergence-free disturbances satisfying the boundary conditions. According to the definition (26), the right hand side of (56) is no larger than Gr_E⁻¹. Therefore, we have

(57)

Furthermore, from (55) and with the same variable substitution as in (56), we have , where the function is the maximum growth rate of the conventional disturbance energy in (27). It follows that

(58)

For any 0 < σ < (Gr_E/Gr)², we have , and the weighted disturbance energy will monotonically decay to zero. This conclusion applies to all finite-amplitude streamwise-independent disturbances. Thus, we have proven that the basic flow is always nonlinearly stable to streamwise-independent disturbances in the vertical buoyancy layer.

Remark 2 Joseph and Hung^[25] have proven that the parallel flows with prescribed boundary values (such as the plane Poiseuille flow and the plane Couette flow) are stable to all streamwise-independent disturbances. Their proof relied on the Poincaré inequality in a finite interval, , where C is a constant. However, there is no Poincaré inequality in the vertical buoyancy layer, which is infinite in the wall-normal direction. Therefore, we have to use a weighted energy method to prove Theorem 2.

To compare the different energy stabilities defined by (24) and (54), we simulate the evolution of the least stable mode in the energy method (α, β)=(0, 0.485 69) in the vertical buoyancy layer with SIMSON^[24], and plot the average disturbance energy E' and in Fig. 10 for Gr=200 and Pr=0.72, which are defined as

(59)

and

(60)

where Ω = [0, L_x] × [0, 30]× [0, L_z] with a resolution of (64, 513, 64) spectral modes. We choose L_z = 2π/0.485 69, whereas L_x can be arbitrary because the disturbance is independent of x.

Because we have known Gr_E = 39.11 when Pr=0.72, we choose σ=0.03 < (Gr_E/Gr)² in the definition (54), which leads to the above definition of . According to the proof of Theorem 2, although the conventional disturbance energy may grow transiently, the weighted disturbance energy with σ < (Gr_E/Gr)² will decay monotonically, which is confirmed in the simulation (see Fig. 10).

6 Conclusions

In this paper, we prove that the inclined buoyancy layer and the vertical buoyancy layer have the same critical Grashof number of energy stability. However, if the Grashof number is defined as Gr' = gγd^*3∆T^*/ν² = Gr/ sin χ^[4] instead of (6), the new critical Grashof number of the energy stability satisfies Gr_E'(χ) = Gr_E' (90°)/ sin χ. Furthermore, we take advantage of the exponential form of the basic flow, using the method of descending exponentials to discretize the eigenvalue equations in the energy stability analysis. The critical Grashof numbers and the critical wavenumbers are calculated with better accuracy than those listed in Ref. [7].

In the band between the energy stability curve and the linear stability curve in Fig. 3(b), the disturbance energy is able to increase transiently^[7]. We calculate the transient growth in the disturbance energy for Pr=0.72 based on the linearized Navier-Stokes equation. When χ ≤ 50°, the optimal disturbances that have the maximum energy growth are streamwise vortices. This conclusion also holds for 55° ≤ χ ≤ 90° and Gr≤ 110. However, when 55° ≤ χ ≤ 90° and Gr ≥ 120, oblique rolls with nonzero streamwise and spanwise wavenumbers are the optimal disturbances. The maximum growth in the disturbance energy over all angles of inclination at subcritical Grashof numbers is S=11.62 for χ=90°.

In the vertical buoyancy boundary layer, we prove that the basic flow is always linearly stable to streamwise vortices. We also define the weighted disturbance energy (54) instead of the original disturbance energy (24) to prove that the basic flow of the vertical buoyancy layer is stable to any finite-amplitude streamwise-independent disturbance. However, this definition of weighted energy cannot be extended either to streamwise-dependent disturbances or to the inclined buoyancy layer.

The energy stability of the inclined buoyancy layer is reminiscent of that of parallel shear flows, where the critical Reynolds numbers of energy stability are less than those of linear stability. Although continuous turbulence exists even when Re < Re_L in plane Poiseuille flows^[26], whether nonlinear effects are able to maintain a finite-amplitude disturbance in the inclined buoyancy boundary layer when Gr_E < Gr < Gr_L is worthy of further study.

Acknowledgements The work was carried out at National Supercomputer Center in Tianjin, and the calculations were performed on TianHe-1(A).