Nomenclature

θm, density inversion parameter;	a, aspect ratio;
Ra, Rayleigh number;	θ, dimensionless temperature;
Pr, Prandtl number;	t, dimensionless time;
T, temperature;	u, dimensionless velocity.
m, wave number;

1 Introduction

Penetrative convection refers to the phenomena whenever the convection in a thermally unstable fluid layer penetrates into the adjacent stable layers. It is commonly known to happen in many geophysical and astrophysical flows^[1]. Typical stellar convection zones are bounded by stably stratified regions, and the fluid motions in the convection zones can penetrate into these stably stratified regions, leading to momentum transport and internal gravity waves. It is also found that the zonal flow in Jupiter's atmosphere originates from the penetrative convection taking place in the deep hydrogen-helium interior^[2]. Although Earth's liquid core is unstable to convect, and generates a magnetic field^[3], the outermost part of Earth's core may be stably stratified^[4]. Thus, a detailed study of penetrative convection is necessary for the understanding of convection in many stars and planets.

The fact that water has a maximum density at T_m=4℃ makes it well suited to study penetrative convection. Consider a Rayleigh-Bénard convection (RBC)^[5-6] system using water as the working fluid, in which the temperature of the bottom plate is higher than T_m while the temperature of the top plate is lower than T_m. The fluid in the lower layer above 4℃ is convectively unstable, while the upper layer below 4℃ is gravitationally stable. When the Rayleigh number Ra of the unstable layer exceeds a critical value, the convection will happen and may penetrate into the upper stably stratified region.

The first important work on the subject of penetrative convection was performed in a layer of water, in which the bottom boundary was maintained at 0℃ and the top boundary was kept at a temperature higher than 4℃^[1]. Thus, a gravitationally unstable fluid layer lies below a layer which is stably stratified. For the convection onset, a stable finite amplitude solution and an unstable small amplitude solution are observed. The convective instability can take place at a finite amplitude, and the convection is stable at the values of Ra less than the critical value of the infinitesimal stability theory. Moore and Weiss^[7] used numerical simulations to extend the finite amplitude solution of penetrative convection to higher Ra, and gave a model explaining the finite amplitude instability with the consideration of the distortion of the mean temperature. Large and Andereck^[8] experimentally studied the onset and pattern formation in the penetrative RBC in water near 4℃, and confirmed the subcritical behavior. Due to the symmetric breaking, two different arrangements of hexagonal convection cells were observed. Hu et al.^[9] investigated the RBC of cold water near 4℃ in a cubical cavity with different thermal boundary conditions on the sidewalls. They found that there were multiple flow pattern coexistence and hysteresis phenomena during the flow pattern transition.

The above works are all performed in two-dimensional (2D) or three-dimensional (3D) rectangular containers. In recent years, the penetrative RBC in cylindrical containers has attracted much attention. Li et al.^[10-11] investigated the pattern formation of the RBC of cold water near 4℃ in a vertical cylindrical container. They showed that the critical Rayleigh number Ra_c increased with the increase in θ_m for the convection onset. Hu et al.^[12] investigated the aspect ratio dependence of the penetrative RBC in vertical cylindrical containers. Their results showed that the aspect ratio had a great effect on the flow patterns, and the number of rolls increased with the increase in the aspect ratio. More recently, Huang et al.^[13] studied the turbulent RBC of cold water near 4℃ in a vertical cylindrical container, and proposed the power law scalings between the average Nusselt number and Ra for different θ_m.

In the previous works, the subcritical behavior of finite amplitude solutions was observed in penetrative RBC. However, it is known that the axisymmetric convection occurs through supercritical pitchfork bifurcation within the Oberbeck-Boussinesq (OB) approximation^[14]. There are two kinds of axisymmetric convection with upward and downward flows at the center for a cylindrical cell. These two kinds of axisymmetric convection possess the same stability characteristics since they are conjugate to each other within the OB approximation^[15-16]. In the penetrative RBC, these two solutions are different, and also have distinctive bifurcation processes due to the symmetry breaking. However, few works are devoted to the bifurcations of finite amplitude solutions in the cylindrical penetrative RBC. The effect of θ_m on the stability of axisymmetric solutions is still unclear. In this paper, the onset of penetrative convection is studied, and the relationship of θ_m and Ra_c is investigated. More importantly, the bifurcation processes based on the axisymmetric solutions are investigated in detail. Besides, particular attention is paid to showing the differences between the stability of the two axisymmetric solutions by diagrams.

2 Problem formulation

The problem is sketched in Fig. 1. The aspect ratio is defined as a=R/H. In the cylindrical coordinates, the domain is (r, φ, z)∈ [0, a]× [0, 2π]× [0, 1]. The cylinder is heated from the bottom with a temperature T_h, and cooled from the top with a temperature T_c. It should be noted that T_h>T_c. The side wall is adiabatic. All the boundaries of the container are rigid. All the material properties remain constant except for the density in the buoyancy term. The temperature-dependent density of cold water is given by Gebhart and Mollendorf^[17] as follows:

(1)

where

H, H²/ν, ν/H, and ρ (ν/H)² are used as the length, time, velocity, and pressure scales, respectively, where ν is the viscosity coefficient. The dimensionless temperature is defined as θ=(T-T_c)/(T_h-T_c). The dimensionless governing equations can be written as follows:

(2)

(3)

(4)

The dimensionless parameters controlling this problem are the Prandtl number Pr=ν/κ, the Rayleigh number Ra=(gγ(T_h-T_c)^q H³)/(νκ), and the density inversion parameter θ_m=(T_m-T_c)/(T_h-T_c), where κ is the thermal diffusivity, and g denotes the acceleration due to gravity. Here, Pr is fixed at 11.57.

The boundary conditions for the temperature and velocity can be expressed as follows:

(5)

(6)

(7)

The temperature distribution is θ(z)=1-z for the motionless conductive state.

3 Numerical method

An improved fractional-step method is used to solve the governing equations in 3D cylindrical coordinates with the second-order accuracy in time and space. This method is an improved version of the algorithm of Verzicco and Orlandi^[18], in which an additional pressure predictor^[19] is carried out at each time step. The results show that the improved method is more robust than the original version, especially for the Newton method.

In the linear stability analysis (LSA), the Jacobian-free Newton-Krylov method^[20] is used to obtain both the stable and unstable solutions. A first-order implicit/explicit Euler scheme is used to calculate the Stokes pre-conditioner, which is the matrix-free inversion of the preconditioned Jacobian for each Newton iteration^[21]. The BICGSTAB algorithm^[22] is used to solve the corresponding linear system, which is free from constructing the Jacobian matrix.

The linear stability equation for the infinitesimal perturbations u ' of a steady state u is

(8)

Thus, the eigenvalues λ_iof the matrix (N_u + L) control the linear stability of the steady state u:

(9)

For a given basic flow, we use the Arnoldi algorithm^[23] to calculate the leading eigenvalues. A small matrix is constructed by time stepping the linearized equations, which represents the action of the Jacobian matrix on the subspace of the leading eigenvectors. The leading eigenvalues and eigenvectors are obtained by the diagonalization of the small matrix.

These methods have been proven to be effective in our previous works^{[15, 24]}. Furthermore, the current numerical results with density inversion are also in good agreement with the results of Hu et al.^[12] and Li et al.^[10-11] (see Tables 1 and 2). Most of our simulations are performed with grids of 90×60×60 in the azimuthal, radial, and vertical directions, respectively.

4 Results and discussion

In the present system, the flow is quiescent at the conductive state at low Ra. When Ra increases above a critical value, the flow will bifurcate into the convective state. In a cylindrical container, the flow patterns at the convective state can be characterized by the azimuthal wave number m. The Ra_c and the corresponding m of the convection onset depend on the aspect ratio a and the density inversion parameter θ_m.

4.1 Effects of the density inversion parameter on convection onset

By the LSA, the convection onset is studied with 0.66≤ a≤ 2. Figure 2 shows the stability curves of the convection onset. At θ_m=0, the Ra_c and the critical wave numbers are similar to those under the OB approximation^[24]. The results also show that the critical wave numbers of the convection onset are affected by θ_m.

Figure 3 shows the normalized Ra_c for a=0.66, 1.1, and 2 and θ_m=0.4, 0.45, and 0.5. The normalizations are carried out by using the corresponding values at θ_m=0. It is seen that twelve different aspect ratios of each θ_m are uniformly distributed between 0.66 and 2. This indicates that, with different aspect ratios, the normalized Ra_c fluctuates within a limited range. Although the values of the normalized Ra_c fluctuate, no significant difference is found in the normalized Ra_c when 0.66≤ a≤ 2. The three curves in Fig. 3 are for the normalized Ra_c at a=0.66, 1.1, and 2. It is clear that all the three curves almost coincide with each other. This indicates that, the normalized Ra_c is strongly affected by θ_m, while the aspect ratio has a minor effect on the normalized Ra_c. In the parameter range of our study, the Ra_c of the convection onset increases when θ_m increases, and it increases more rapidly when θ_m>0.4.

It is obvious that the stably stratified fluid layer will become thicker when θ_m increases, and the stably stratified fluid layer will suppress the convection onset of the lower unstable fluid layer. Therefore, the flow is more stable at higher θ_m. Supposing that there is an effective Rayleigh number in which the effect of the stably stratified fluid layer is considered. Based on the definition of Ra, the temperature difference and the effective height of the convective region should be reconsidered. For a penetrative RBC system, the effective temperature difference can be written as follows:

(10)

Now, consider the effective depth of the unstable fluid layer in the conductive solution. We note that, the convection in the unstable fluid layer will penetrate into the upper stably stratified fluid layer, and there is a penetrative depth in the penetrative RBC. We account for this effect by introducing a coefficient f in the following expression with the effective height of the convection layer:

(11)

The coefficient f of the penetrative depth is a function of θ_m and f ≥ 1. When θ_m=0, we have f=1, due to the restriction of the upper boundary on the penetrative depth. When θ_m increases, the stably stratified layer becomes thicker, and the influence of the upper boundary becomes weaker. We expect that f tends to be constant for large enough θ_m. is expressed as follows:

(12)

Without the consideration of the effect of side walls, the critical should be a constant value c, i.e.,

(13)

Thus, the relationship between θ_m and Ra_c can be written as follows:

(14)

Figure 4 shows the normalized critical Rayleigh numbers as a function of 1-θ_m for a=0.66, 1.1, and 2. The dashed lines in Fig. 4 denote the predicted normalized Ra_c based on Eq. (14) with f=1.5. The specific mathematical expression is

(15)

The results show that the normalized Ra_c distributes around the dashed lines when 1-θ_m < 0.6, while it tends to 1 when 1-θ_m>0.6. When 1-θ_m increases, the upper boundary will restrict the depth of the penetrative convection. Therefore, the normalized Ra_c tends to 1 when 1-θ_m increases to 1.

4.2 Bifurcations of the axisymmetric solutions

It is found that the conductive solution bifurcates to convective axisymmetric flow at a=1 when θ_m is from 0 to 0.5. Via the DNS with different initial conditions, we find that there are still two kinds of steady axisymmetric solutions for all considered θ_m with the central fluid moving upward (A) and downward (B), respectively, although the top-bottom symmetry is broken for the present system. The coexistence of the solutions A and B is observed in a certain range of Ra above the threshold of the first bifurcation. Figure 5 shows the typical temperature fields of the solutions A and B at Ra=3 800 and θ_m=0.2. Different from the prior results within the OB approximation, the stability properties of the two solutions are rather different due to the symmetry breaking. Moreover, when Ra increases, the bifurcation processes and flow patterns after the onset are also quite different.

Using the solutions A and B as the initial conditions, we first calculate Ra_c by gradually decreasing Ra in the DNS. Figure 6 shows the Ra_c of the solutions A and B for the first bifurcation. It is obvious that the Ra_c of the solution A is close to that of the solution B for small θ_m, but their difference becomes more evident when θ_m increases. For instance, at θ_m=0, the Ra_c of the solution A is 2 189, while it is 2 206 for the solution B. However, at θ_m=0.5, Ra_c is 15 878 for the solution A, while it is increased up to 18 829 for the solution B.

Figure 7 shows the bifurcation diagram near the threshold of the first bifurcation based on the vertical velocities of the monitoring point (φ=0, r=0.5, z=0.5) at θ_m=0.35. Different from the previous results in the penetrative convection, the coexistence of solutions with a small stable amplitude and a finite amplitude is observed, and the conductive solution bifurcates to the solution B with a small stable amplitude instead of the solution A with a finite amplitude when Ra increases. It is clearly shown that the values of | u_z | of the solutions A and B are not equal. Within the OB approximation, it has already been found that the two kinds of axisymmetric solutions appear through a supercritical pitchfork bifurcation due to the top-bottom symmetry, which are conjugate to each other. However, in the penetrative RBC, the onset of the axisymmetric convective solutions will occur through a trans-critical bifurcation, due to the fact that the top-bottom symmetry is broken because of the non-linear density and temperature relationship.

4.3 Stability of the axisymmetric solutions

We now turn to study the stability properties of the axisymmetric solutions. Figure 8 shows the pattern evolutions of the solutions A and B after the convection onset. The stability curves are calculated by the LSA with either decreasing or increasing Ra. The base flows for Figs. 8(a) and 8(b) are the solutions A and B, respectively. The results of the stability analysis show that the bifurcations of both the axisymmetric solutions to the 3D convection are stationary. In Fig. 8(a), it is shown that a stable region of the solution A (AS) is bounded by two stability curves. Taking the solution A in the AS region as the base flow, the decrease in Ra will result in a conductive state, while, when Ra increases, the solution A will bifurcate into the steady 3D convection. In Fig. 8(b), a region of the stable solution B is non-existent (NE), and two stable regions (AS) and two unstable regions (AU) of the solution B are separated by four stability curves. Resembling the prior work within the OB approximation for low Pr^[16], we here find a stable region of the solution B beyond the secondary bifurcation at Pr = 11.57.

From Fig. 8, we can also find that, when θ_m increases, the values of Ra_c of both the solutions increase rapidly, and θ_m will affect the azimuthal wave number m. Moreover, the effect of θ_m on the pattern evolutions of the solution A is considerably different from that of the solution B. Specifically, the Ra_c of the solution B is more sensitive to the increase in θ_m than that of the solution A. For instance, the Ra_c of the solution A is 21 226 at θ_m=0 and 38 328 at θ_m=0.5, but the highest Ra_c of the solution B is 13 858 at θ_m=0 and 142 395 at θ_m=0.5.

In addition, the critical m of the solution A is more sensitive to the change in θ_m. When θ_m changes from 0 to 0.5, the critical m of the solution A changes from 1 to 4, to 3, and then to 1 sequentially. While for the solution B, the critical m of the secondary bifurcation remains to be 2 until θ_m reaches 0.5, and then becomes 3. In the stable region of the solution B beyond the secondary bifurcation, when θ_m changes from 0 to 0.5, the critical m is 1 when Ra increases, while is 2 when Ra decreases.

The existence of multiple active patterns can lead to complicated nonlinear flow evolutions. At Ra=37 000 and θ_m=0.4, the flow transition process from the solution A obtained by the DNS is shown in Fig. 9(a), where the time series of the azimuthal velocity of a monitoring point (φ=0, r=0.5, z=0.5) is given. For this parameter combination, the solution A is unstable, and the m=3 mode grows in the amplitude and saturates due to the nonlinear effect. This is consistent with the results of the LSA, which shows that the most unstable mode of the solution A is of the azimuthal wave number m=3 for θ_m=0.4. The saturated m=3 pattern is actually unstable, and will further transit to a m=2 pattern as time goes on. The typical flow structures appearing in this transition process are shown in Figs. 9(b), 9(c), and 9(d).

5 Conclusions

In summary, the bifurcation processes are studied in detail by the LSA and DNS in the penetrative RBC for cold water near its maximum density in the cylindrical containers at 0.66≤ a≤ 2.0 and 0≤ θ_m≤ 0.8. For the convection onset, Ra_c increases with the increase in θ_m, and it increases more rapidly when θ_m>0.4. The effect of the aspect ratio on the normalized Ra_c is limited in the parameter range of our study. Furthermore, the relationship between θ_m and the normalized Ra_c is formulated as follows:

which is in good agreement with the results of the LSA when θ_m>0.4.

In such a top-bottom symmetry breaking system, we still identify two kinds of qualitatively different steady axisymmetric solutions beyond the first bifurcation. The onset of the axisymmetric convection occurs through a trans-critical bifurcation rather than the supercritical pitchfork bifurcation within the OB approximation. Moreover, the subcritical behavior is observed for the solution A, while is absent for the solution B. Interestingly, we find that the pattern evolution based on the solution A is considerably different from that of the solution B, and a stable region of the solution B beyond the secondary bifurcation is identified at Pr=11.57. Furthermore, the Ra_c of the solution B is more sensitive to the increase in θ_m than that of the solution A, while the critical m of the solution A is more sensitive to the variation of θ_m than that of the solution B.