Geometric Sensitivity and Chaos of Thermal Convection in a Horizontal Annular Cavity by the Lattice Boltzmann Method

1  Introduction

Thermal turbulence, also known as buoyancy-driven turbulence, is a fundamental type of turbulence that occurs widely in both natural and engineering systems. Compared with other typical turbulence types, such as wall-bounded turbulence and free shear turbulence, thermal turbulence exhibits richer and more universal flow instabilities due to the bidirectional coupling between velocity and temperature fields in the Navier–Stokes (NS) equations and the influence of the buoyancy force in the external forcing term. Among classical models of thermally driven flow instabilities, Rayleigh–Bénard (RB) convection has been the most extensively studied [1]. Research on RB systems began in the early 20th century with Bénard’s convection experiments [2], and was later advanced by Rayleigh in 1916 [3], who established a mathematical model for thermally driven convection based on the NS equations. Since then, investigations of the classical RB system have continued to the present day, with experiments and numerical simulations mutually complementing each other. Considerable progress has been made in understanding turbulent structures, heat transport, statistical properties of turbulent fluxes, and boundary-layer dynamics [4]. Furthermore, the RB framework has been extended to a variety of research directions, including inclined thermal turbulence [5], viscoelastic fluid thermal turbulence [6,7], multilayered convection [8], and lateral thermal convection [9].

In the classical RB convection model, the temperature gradient is aligned with the direction of gravity, and convection occurs when the Rayleigh number exceeds the critical value of 1708 [4]. In contrast, the horizontal annular configuration introduces curvature effects, causing the temperature gradient to distribute along the curved surface. Consequently, thermal convection in an annular cavity cannot be simply classified as classical RB convection or lateral thermal convection, as it exhibits distinct geometric characteristics and flow instability behaviors [10]. For both classical and extended RB convection models, geometric features play a crucial role in determining flow instabilities, and many studies have investigated flow and heat-transfer stability from this geometric perspective.

In classical RB systems, numerous studies have focused on the effects of geometric parameters such as the aspect ratio, inclination, and length-to-width ratio on thermal convection. For example, Huang and Xia [4] experimentally investigated the influence of cavity thickness on the reversal of large-scale circulation in a quasi-two-dimensional rectangular cavity. Their results showed that the confinement in the thickness direction enhances wall friction effects, which in turn weakens the intensity of the large-scale circulation and allows more plumes to pass through the cavity, thereby significantly increasing the flow reversal frequency. Xu et al. [7] conducted numerical simulations of convection in two-dimensional annular and rectangular cavities, revealing that due to differences in cavity geometry, the time required for plumes to self-organize into circulations in the 2D annular cavity is several to tens of times longer than that in rectangular cavities. Pan et al. [11] used numerical simulations to study the influence of different length-to-width ratios on heat transfer and flow characteristics in two-dimensional square cavities, with length-to-width ratios ranging from 0.3 to 8 and Rayleigh numbers from 5 × 103 to 108. Their study showed that, at small length-to-width ratios, the critical Rayleigh number for convection onset decreases with increasing temperature difference, whereas at large length-to-width ratios, it increases with increasing temperature difference. They also investigated the dependence of the Nusselt and Reynolds numbers on the cavity aspect ratio.

For Rayleigh–Bénard convection systems with lateral heating, researchers have primarily focused on the effects of geometric factors such as aspect ratio, heating method, and heating location on convection patterns. Batchelor [10] was among the first to consider such a system. Hu et al. [9] systematically investigated the effects of geometric confinement and fluid properties on the stability and heat transfer characteristics of thermal convection in laterally heated cavities through numerical simulations. By varying key parameters such as aspect ratio and fluid properties, they analyzed the evolution of flow structures and stability behavior, and demonstrated that both geometric constraints and fluid properties play crucial roles in determining the onset of instability and the overall heat-transfer performance of the system. Mac Huang and Zhang [12] introduced lateral wall heating into classical RB systems to investigate its influence on heat transfer and flow structures, and found that lateral heating significantly enhances the system’s heat-transfer efficiency.

Thermal convection in horizontal annular cavities exhibits both thermal and flow instabilities due to its unique geometric features under different computational conditions. Thermal instabilities are typically characterized by buoyancy effects, while flow instabilities are described by inertial forces, with their relative dominance and interaction varying across different parameter ranges. The radius ratio is the most critical geometric factor influencing convection patterns. Bishop et al. [13] conducted detailed smoke-visualization experiments to study the oscillatory characteristics of natural convection in horizontal annular spaces, with Grashof numbers ranging from 290 to 2.7 × 106 and four radius ratios of 1.23, 1.85, 2.46, and 3.69. Powe et al. [14] first mapped the flow patterns of thermal convection in horizontal annular cavities as functions of Rayleigh number and radius ratio R, using radius ratios of 1.2, 1.42, 1.6, 2.0, and 2.4, and identified three distinct unstable flow regimes. Usman et al. [15] numerically studied the chaotic transition of natural convection in a horizontal annular space at low Rayleigh numbers. They set parameters: radius ratio of 2, Prandtl number of 0.1, and Rayleigh number between 3500 and 5000. Their focus was on the evolution of the flow state from steady to transient, revealing that chaotic behavior can occur even at relatively low Rayleigh numbers under specific parameter conditions. Ahlers et al. [16] experimentally studied the effect of aspect ratio on heat transfer in a cylindrical Rayleigh-Bernard convection system. They systematically analyzed the change of Nusselt number with aspect ratio and demonstrated that geometric configuration plays a crucial role in determining overall heat transfer efficiency. Zhu and Zhou [17] numerically studied the flow structure of turbulent Rayleigh-Bernard convection in annular units with aspect ratios equal to or greater than 1. Their focus was on large-scale circulation and plume dynamics, showing that aspect ratio significantly affects the organization and evolution of turbulent structures. In a recent study, Wang et al. [18] numerically simulated centrifugal annular RB convection to examine the impact of radius ratio on system behavior, finding that under the combined influence of curvature and centrifugal forces, larger radius ratios lead to weaker convection intensity but enhanced heat-transfer efficiency, while smaller radius ratios result in greater deviation of the overall temperature from the arithmetic mean. Farkach et al. [19] applied the lattice Boltzmann method to study the velocity and temperature fields in annular cavities of different sizes, considering Rayleigh numbers from 104 to 5 × 105 and radius ratios of 5, 2.5, 1.67, and 1.25. Their results indicate that the Nusselt number increases with both Rayleigh number and radius ratio, leading to more vigorous convection.

Despite the extensive studies on Rayleigh–Bénard convection and annular convection systems, most existing work has primarily focused on heat transfer characteristics or low-Rayleigh-number flow stability, and systematic investigations of the geometric sensitivity and chaotic characteristics of annular convection from the perspective of nonlinear dynamics remain limited. Therefore, the present study aims to investigate the geometric sensitivity and nonlinear dynamical evolution of natural convection in horizontal annular cavities using the lattice Boltzmann method (LBM) [20], building on the thermal LBM with a double-distribution-function model previously applied by Xu et al. to study the stability of natural convection in eccentric annular cavities [21]. In particular, the transition from steady to chaotic regimes is systematically analyzed by combining the maximum Lyapunov exponent, runs test, and phase-space trajectories, with the results further compared with those obtained in eccentric annular systems to analyze the stability of convective solutions. The main novelties of this study can be summarized as follows: (1) A systematic identification of critical Rayleigh numbers corresponding to different dynamical regimes (steady, periodic, quasi-periodic, and chaotic) under varying radius ratios; (2) A combined application of nonlinear dynamics tools to characterize chaotic convection in annular systems; and (3) A comparative analysis between concentric and eccentric annular configurations, highlighting the stabilizing effect of geometric symmetry. The results show that decreasing the radius ratio significantly lowers the critical Rayleigh numbers for flow transitions, indicating enhanced geometric sensitivity. In addition, concentric annular configurations require higher Rayleigh numbers to reach the same dynamical states compared with eccentric annuli, demonstrating the stabilizing effect of geometric symmetry. These contributions provide new insights into the role of geometry in controlling flow instability and chaotic behavior in thermally driven convection systems.

2  Problem Description

2.1 Physical Model

Natural convection heat transfer in a horizontal concentric annular cavity is numerically investigated in the present study. A two-dimensional physical model is adopted, as shown in Fig. 1. The outer and inner radii of the annulus are denoted by Ro and Ri, respectively. Four radius ratios, Ro/Ri = 2.6, 1.6, 1.4, and 1.2, are considered and are referred to as Model 1, Model 2, Model 3, and Model 4, respectively. Specifically, the selected radius ratios are derived from the classic study by Powe et al. [14], in this study, we selected four representative radius ratios that can cover a wide range of geometric constraints, from narrow to wide. The inner and outer circular walls are maintained at constant temperatures Th and Tc, respectively, with Th > Tc, generating buoyancy-driven flow within the enclosure. The dimensionless temperature is defined as θ=(T−Tc)/(Th−Tc), and the gravitational acceleration is denoted by G. The center of the annular cavity is located at (65, 65) in the computational domain. To analyze the temporal evolution of flow and thermal characteristics, four monitoring points are selected at (66.9, 90.4), (66.9, 107), (65, 121), and (66.9, 120.2).

Figure 1: Schematic diagram of the physical model.

The working fluid inside the enclosed cavity is assumed to have a Prandtl number of Pr = 0.7, corresponding to air under standard conditions. The flow is considered incompressible, but satisfies the Boussinesq approximation, in which density variations are neglected except in the buoyancy term.

2.2 Double-Distribution-Function Model

The physical problem under consideration is a two-dimensional natural convection heat transfer system. Under the condition of relatively small temperature variations, the fluid is treated as incompressible, and the Boussinesq approximation is adopted. All the governing equations and related formulations employed in this study follow the standard forms of natural convection theory and the lattice Boltzmann method, as widely documented in the literature [20]. Accordingly, the governing macroscopic equations can be expressed as follows:

∇⋅u=0,(1)

∂u∂t+∇⋅(uu)=−∇p+υ∇2u+gβ(T−T0),(2)

∂T∂t+∇⋅(uT)=χ∇2T.(3)

In the above equations, υ denotes the kinematic viscosity and χ represents the thermal diffusivity. Both parameters are related to the dimensionless relaxation times τf and τT as follows:

υ=13c2(τf−12)δt,χ=13c2(τT−12)δt.(4)

The numerical simulations are performed using the coupled double-distribution-function lattice Boltzmann model (CLBGK) proposed by Guo et al. [22], and the discrete velocity space is constructed based on the D2Q9 lattice. The evolution equation for the velocity field in the D2Q9 model, including the external force term, is given as follows:

fα(x+cαδα,t+δt)−fα(x,t)=−1τf[fα(x,t)−fαeq(x,t)]+δtFα.(5)

Here, cα=cei, c=Δx/δt, Δx and δt represent the lattice spacing and the time increment, respectively, while τf denotes the dimensionless relaxation time describing the relaxation of the distribution function toward the local equilibrium state.

A modified equilibrium distribution function is then constructed as follows:

fαeq={ρ−5p3c2+ρs0(u),  α=0p3c2+ρsα(u),  α=1,2,3,4p12c2+ρsα(u),  α=5,6,7,8,(6)

sα(u)=ωα[cα⋅ucs2+9(cα⋅u)22cs4−3u22cs2].(7)

In the above expression, ωα denotes the weighting factor, with ω0=4/9, ω1−4=1/9, ω5−8=1/36. The macroscopic density ρ and velocity u are obtained from the velocity distribution functions, and the corresponding macroscopic variables satisfy:

ρ=∑α=08fαeq,(8)

ρu=∑α=08cαfα+δt2F,(9)

p=ρ3c25[∑α=08fα+s0(u)].(10)

The velocity field and the temperature field are coupled through the external force term Fα. Accordingly, the external force term associated with the equilibrium temperature distribution function is defined as follows:

Fα=(1−12τf)ϖα[3eα−uc2+9(eα⋅u)c4cα]⋅F.(11)

In this expression, F represents the gravitational body force, which, under the Boussinesq approximation, is expressed as follows:

F=−gβ(T−T0).(12)

Here, g denotes the gravitational acceleration, and β is the thermal expansion coefficient.

The evolution of the temperature field is calculated using the following approach:

Tα(x+cαδt,t+δt)−T(x,t)=−1τT[Tα(x,t)−Tαeq(x,t)].(13)

In this case, the temperature distribution function is discretized using the D2Q9 lattice, which provides higher numerical accuracy compared with the four-velocity model. The corresponding equilibrium distribution function is given as follows:

Tαeq=ωαT[1+3cα⋅uc2+9(cα⋅u)22c4−3u22c4].(14)

The macroscopic temperature is obtained from the above expression as follows:

T=∑α=08Tα.(15)

For the natural convection heat transfer problem, two dimensionless parameters characterizing the system are defined as follows:

Pr=υa,Ra=gβΔTH3Prυ2.(16)

Here, ΔT denotes the temperature difference between the hot and cold walls, and H is the characteristic length. In general, when the thermal expansion coefficient β is unknown, the dimensionless relaxation time τf and τT cannot be fully determined. In this case, a third dimensionless parameter, the Mach number (Ma), is introduced:

Ma=uccs.(17)

Here, uc is the characteristic velocity, defined as uc=gβΔTH. For the incompressible flow considered in this study using the D2Q9 model, the Mach number is set to Ma = 0.1. This allows the determination of:

τf=12+MaH3PrcδtRa,(18)

τT=12+3MaHcδtRaPr.(19)

The above numerical model, through the Chapman–Enskog expansion, recovers the macroscopic governing equations given in Eqs. (1)–(3).

2.3 Boundary Conditions

The physical boundaries considered in this study are irregular curved surfaces, while the computational grid is uniformly and symmetrically distributed. As a result, the boundary points rarely coincide exactly with the lattice nodes, and appropriate interpolation or approximation is required to ensure that the boundary conditions are satisfied. In this work, a curved boundary treatment combining spatial interpolation with non-equilibrium extrapolation, as proposed by Guo et al. [22], is employed. The basic idea is to decompose the distribution function at the boundary node into equilibrium and non-equilibrium components. The equilibrium part is replaced by a virtual equilibrium distribution function, while the non-equilibrium component is obtained through interpolation from neighboring fluid nodes. Fig. 2 presents a local magnification of the curved boundary. Let eα denote the lattice velocity pointing from the fluid node toward the solid node, and eα− denote the lattice velocity in the opposite direction. During the collision-streaming cycle of the particles, collisions occur only at the nodes themselves, whereas streaming transfers the corresponding distribution functions along the lattice velocity directions to neighboring nodes. Considering the streaming process between a fluid node xf and a boundary node xb, the unknown distribution function entering the fluid node from the boundary needs to be determined. The unknown distribution function fα+(rb,t) at the boundary node can be expressed as:

fα+(rb,t)=fα+.eq(rb,t)+fα+,neq(rb,t).(20)

Figure 2: Schematic diagram of a curved boundary.

In the above expression, fα+(rb,t) represents the equilibrium distribution function at the boundary node, which is explicitly given as follows:

fα+,eq(rb,t)=Feqα(ρb,ub,t).(21)

Here, ρb and ub are unknown and need to be obtained by interpolation from the three nodes xw, xf, and xff.

When q≥qc:

ρ(rb)=ρ(rf),ub=uw+(q−1)ufq.(22)

When q<qc:

ρ(rb)=ρ(rf),ub=uw+(q−1)ufq+(1−q)2uw+(q−1)uff1+q.(23)

In the above expression, the variable q represents the relative distance from the fluid node to the actual physical boundary, defined as q=|rf−rw|/|rf−rb|,0≤q≤1. The parameter qc serves as the criterion for the interpolation scheme, with a value in the range (0, 1); here, qc = 0.75.

In Eq. (20), fα+,neq(rb,t) represents the non-equilibrium component at the boundary node. Since solid nodes do not participate in the collision process, the non-equilibrium part is obtained by interpolation from the distribution functions at the fluid nodes xf and xff.

When q≥qc:

fα+,neq(rb,t)=fα(rf,t)−fαeq(rf,t).(24)

When q<qc:

fα+,neq(rb,t)=q[fα(rf,t)−fαeq(rf,t)]+(1−q)[fα(rff,t)−fαeq(rff,t)].(25)

Once fα+(rb,t) and fα+,neq(rb,t) are obtained according to the interpolation scheme, the unknown distribution function fα+(rb,t) at the boundary node rb can be determined. This approach achieves second-order accuracy in both space and time and exhibits good numerical stability.

2.4 Nusselt Number

The local Nusselt number is defined as [23]:

NuL=hHχ=−∂Θ∂N|wall(26)

Because the model boundary is a curved surface, the wall normal temperature gradient cannot be directly calculated through mesh nodes. Therefore, this paper uses multiple sampling points for interpolation and then fits the gradient. The specific method is as follows:

1 Five sampling points are uniformly selected along the in-wall normal direction. The temperature value of each sampling point can be obtained from its four adjacent grid nodes through bilinear interpolation. The obtained radial temperature distribution T(r) is fitted to T(r)=ar2+br+c using a least-squares quadratic polynomial. The wall normal temperature gradient is obtained by analytical differentiation:

∂T∂n|wall=±dTdr|r=rw=±(2arw+b)(27)

The sign in the formula is defined according to the normal direction of the wall (negative for the outer wall and positive for the inner wall).

2 Select multiple angles uniformly along the entire circumference to calculate NuL, obtain the local Nusselt number distribution NuL(ϕ), and then obtain Nu¯ by integral averaging:

Nu¯=12π∫02πNuL(ϕ)dϕ(28)

The inner wall surface Nu¯ is essentially the arithmetic sum of the contributions of each local heat source, while the outer wall surface Nu¯ reflects the overall heat transfer efficiency of the system. The Nusselt numbers given below in this paper refer to the outer wall surface Nu¯.

2.5 Validation

This paper uses the Coupled Double Distribution Function Model (CLBGK model) proposed by Guo et al. [22] to numerically simulate the natural convection heat transfer process of concentric rings. To verify the reliability and applicability of the established program and mathematical model, the calculation conditions were set to Ra=4.8×104 and Pr=0.7. The dimensionless temperature distribution calculated by the program was compared and analyzed with the experimental results of Kuehn and Goldstein [24]. As shown in Fig. 3, the dense region inside the concentric rings is concentrated in the outer boundary layer, with a large vertical temperature gradient, exhibiting an upward thermal plume pattern. The isothermal line densities at the bottom of the inner and outer rings are high, with good symmetry and a basically consistent distribution trend.

Figure 3: Distribution of isotherms. Left: reference; Right: this paper.

Subsequently, the natural convection of the concentric rings was simulated. The radial dimensionless temperature distribution of the concentric rings under calculation conditions Ra=5×104 and Pr=0.7 is compared with the experimental results in reference [24], as shown in Fig. 4. This figure shows that the simulation results and experimental results have very small errors and are basically consistent, confirming the correctness of the program.

Figure 4: Compares the radial dimensionless temperature distribution with that in reference [24].

3  Numerical Determination of Chaotic States

3.1 Maximum Lyapunov Exponent Map

Chaotic systems exhibit extreme sensitivity to initial conditions, which is manifested by the exponential divergence of nearby trajectories in phase space. Given an initial condition x0 and a neighboring point x0+δ0 with a very small initial separation δ0, let δn denote the separation after n iterations. If |δn|≈|δ0|enλ, then λ is defined as the Lyapunov exponent. Taking the logarithm and noting δn=fn(x0+δ0)−fn(x0), one obtains:

λ=1nln|(fn)(1)(x0)|.(29)

The Lyapunov exponent is defined based on the exponential growth rate of an initially small perturbation. Let the initial perturbation be |x0−y0|=δ0. After one iteration, the perturbation becomes |x1−y1|=|f(x0)−f(y0)|=A1|x0−y0|, where A1=|f(x0)−f(y0)||x0−y0|. Since the initial perturbation is very small, x0→y0 is negligible, and thus A1=(f1)(1)(x0). After two iterations, the perturbation is |x2−y2|=|f(x1)−f(y1)|=A2|x0−y0|, with A2=|f(x1)−f(y1)||x1−y1|∗|f(x0)−f(y0)||x0−y0|, yielding A2=f′(x0)f′(x1)=(f2)(1)(x0). By continuing this process, one obtains:

An=|∏i=0n−1f(1)(xi)|=(fn)(1)(x0).(30)

As n→∞, the Lyapunov exponent can be obtained as:

λ=limn→∞{1n∑i=0n−1ln|f(1)(xi)|}.(31)

It should be noted that the above formulation of the Lyapunov exponent is derived for discrete dynamical systems. In the present study, the Lyapunov exponent is evaluated from time series obtained from the LBM simulations. The dimensionless temperature at selected monitoring points is used to construct the time series. A phase-space reconstruction is carried out using the time-delay embedding method, where the delay time is determined from the autocorrelation function and the embedding dimension is selected based on the false nearest neighbors criterion. The maximum Lyapunov exponent is then computed using a standard algorithm that evaluates the average exponential divergence rate of neighboring trajectories in the reconstructed phase space. This procedure enables a consistent application of Lyapunov analysis to CFD data and ensures that the computed exponent reflects the intrinsic dynamical characteristics of the system.

Any attractor has at least one negative Lyapunov exponent. Stable fixed points and limit cycles cannot possess positive Lyapunov exponents, whereas a chaotic attractor has at least one positive Lyapunov exponent. Therefore, the Lyapunov exponent can be used to determine whether a system is in a chaotic state. Fig. 5 shows the spectra of the maximum Lyapunov exponents for Models 1 and 2 under different Rayleigh numbers (Ra). It can be seen that when Ra = 2.5 × 106 and Ra = 5.5 × 105, the maximum Lyapunov exponents of both systems are positive. Based on these Lyapunov exponents, it can be concluded that the two seemingly disordered systems have indeed developed into chaotic states at the corresponding Rayleigh numbers.

Figure 5: Lyapunov exponent maps of Models 1 and 2 at different Rayleigh numbers.

3.2 Run Test for Randomness

The runs test is a method for testing randomness, also referred to as a test of continuity. The basic idea is to select a threshold value within a dataset to divide the data into two categories. The threshold can be the mean, median, mode, or other system-related variables. If n consecutive data points belong to the same category, this segment of data is called a run (n ≥ 1). The runs test evaluates whether preceding observations influence subsequent ones. The null hypothesis H0 assumes that the values of the test variable occur randomly. By calculating the test statistic Z and comparing it with the corresponding values from the standard normal distribution at a given significance level, a significance test is performed to either accept or reject the null hypothesis.

Let U denote the total number of runs in a dataset, and let m and n be the numbers of observations in the first and second categories, respectively. When m + n = N > 20, the normal approximation with continuity correction can be used to calculate the Z value as follows:

Z=U+0.5−1−2mn/N2mn(2mn−N)/N2/(N−1).(32)

The corresponding asymptotic significance p value can be obtained from the standard normal distribution table using the calculated Z value. By comparing p with the chosen significance level σ, the null hypothesis is rejected if P<σ, indicating that the test variable is not random.

As shown in Tables 1 and 2, the asymptotic significance p values of the monitored temperatures were calculated for the cases of radius ratio R = 2.6 with Rayleigh number Ra = 2.5 × 106 and radius ratio R = 1.6 with Ra = 5.5 × 105. In the runs test, a threshold must be selected to divide the time series into two categories: values above and below the threshold. Considering that the time series of a natural convection system in nonlinear oscillatory and chaotic states may be asymmetrically distributed, using a single statistic as the threshold could introduce bias. Therefore, the mean, median, and mode of the series were respectively chosen as thresholds for the runs test. The mean reflects the overall statistical level of the series; the median is robust against skewed distributions and outliers; and the mode corresponds to the most frequently occurring state, representing the dominant dynamical feature of the system. In all cases, the resulting p values were 0 (far below the commonly used significance level of 0.05), indicating that the data series significantly deviate from randomness. In other words, the system solutions are not random fluctuations but exhibit clear intrinsic regularity or correlation, which is the characteristic of chaotic attractors. This further confirms that the two systems have developed into chaotic states at the corresponding Rayleigh numbers. It should be noted that the reported p values equal to zero do not indicate an exact probability of zero, but rather result from numerical truncation in the statistical calculation. Owing to the extremely large absolute values of the Z statistics, the corresponding tail probabilities are far below the minimum resolution of double-precision floating-point arithmetic and are therefore reported as zero. This implies a very strong rejection of the null hypothesis and confirms that the time series exhibits significant non-randomness and intrinsic dynamical correlation.

The significance level is commonly set at 0.05, which implies that even if the null hypothesis is true, there is a 5% probability of rejecting it due to random error. In other words, if the p value is less than 0.05, there is statistically sufficient evidence to conclude that the observed results are not caused by random fluctuations but instead reflect some actual effect or trend. That is, the data have a 95% probability of conforming to the expected random binomial distribution, and the null hypothesis is rejected when p < 0.05.

4  Stability of Natural Convection Solutions in Concentric Annular Cavities at Different Radius Ratios

4.1 Stable Steady-State Solution

Fig. 6 shows the temperature and flow fields of the four models under their respective Rayleigh number conditions. Panels (a) and (b) correspond to Model 1 at Ra = 5 × 105, showing the temperature and flow fields; panels (c) and (d) correspond to Model 2 at Ra = 1.5 × 105; panels (e) and (f) correspond to Model 3 at Ra = 1 × 104; and panels (g) and (h) correspond to Model 4 at Ra = 2 × 103. From the overall distribution patterns, it can be observed that all four models form well-defined left–right symmetric structures at the given Rayleigh numbers. Neither the temperature fields nor the flow fields exhibit any symmetry-breaking phenomena, indicating that natural convection can still sustain a stable and symmetric circulatory pattern under the current Rayleigh number conditions.

Figure 6: Temperature and flow fields of each model at the corresponding Rayleigh numbers: (a,c,e,g) temperature fields; (b,d,f,h) flow fields.

From the dimensionless temperature–time curves at the monitoring points shown in Fig. 7, it can be observed that the temperature gradually increases during the initial stage and eventually reaches a stable value, indicating that the system has attained a steady state. This is further corroborated by the phase-space trajectories at the monitoring points in Fig. 8. For Models 1 and 2, the trajectories in phase space eventually collapse to a single point in the phase plane, demonstrating that the solutions converge to a stable steady state. Similarly, the phase-space trajectories of Models 3 and 4 also converge to a fixed point, known as a steady attractor, with the orbits being drawn toward this point. This indicates that the system state no longer changes and has reached a stable condition. Therefore, under the corresponding Rayleigh numbers, all models achieve state convergence after a brief transient period. The system solutions are unique and stable, exhibiting no temporal oscillations or complex dynamical behavior.

Figure 7: Dimensionless temperature at the monitoring points of each model at the corresponding Rayleigh numbers.

Figure 8: Phase portraits in the u–v plane at the monitoring points of each model at the corresponding Rayleigh numbers: (a) Model 1; (b) Model 2; (c) Model 3; (d) Model 4.

4.2 Periodic Oscillatory Solution

As the Rayleigh number increases, the flow intensity and heat transfer capability within the annular space are significantly enhanced, gradually destabilizing the original steady-state equilibrium and driving the system toward periodic oscillatory behavior. When the Rayleigh numbers reach Ra = 1.5 × 106 for Model 1, Ra = 3 × 105 for Model 2, Ra = 1.25 × 104 for Model 3, and Ra = 2.5 × 103 for Model 4, the temperature and flow fields of the models are shown in Fig. 9. Panels (a), (c), (e), and (g) illustrate the temperature fields, while panels (b), (d), (f), and (h) show the corresponding flow field distributions. Compared with the symmetric structures observed at lower Rayleigh numbers, all four models now exhibit pronounced asymmetry in both temperature and flow fields, indicating that the system has transitioned from a stable convective mode to an oscillation-dominated unsteady state.

Figure 9: Temperature and flow fields of each model at the corresponding Rayleigh numbers: (a,c,e,g) temperature fields; (b,d,f,h) flow fields.

To further reveal the dynamical characteristics of the oscillations, Fig. 10 presents the dimensionless temperature–time curves at spatial monitoring points. The results indicate that, following a brief initial transient, the temperatures evolve into regular oscillations with fixed periods and stable amplitudes, confirming that the system has entered a periodic natural convection state. The corresponding phase-space trajectories are shown in Fig. 11, where the trajectories of all models form closed limit-cycle structures. This further demonstrates that the system’s dynamical attractor has transformed from a steady-state attractor to a limit-cycle attractor, representing a typical periodic oscillatory solution.

Figure 10: Dimensionless temperature at the monitoring points of each model at the corresponding Rayleigh numbers.

Figure 11: Phase portraits in the u–v plane at the monitoring points of each model at the corresponding Rayleigh numbers: (a) Model 1; (b) Model 2; (c) Model 3; (d) Model 4.

4.3 Quasi-Periodic Oscillatory Solution

By further increasing the Rayleigh number beyond the periodic oscillatory state, the flow structures and heat transfer characteristics of the system continue to undergo significant changes. Fig. 12 shows the temperature and flow fields for Model 1 at Ra = 2.1 × 106, Model 2 at Ra = 4.5 × 105, Model 3 at Ra = 1.8 × 104, and Model 4 at Ra = 5 × 103. Panels (a), (c), (e), and (g) correspond to the temperature fields of the four models, while panels (b), (d), (f), and (h) display the respective flow fields. The results indicate that, as the driving force is further enhanced, the internal flow intensity and structural complexity of all models increase, with both temperature and flow fields exhibiting more complex unsteady features.

Figure 12: Temperature and flow fields of each model at the corresponding Rayleigh numbers: (a,c,e,g) temperature fields; (b,d,f,h) flow fields.

The dimensionless temperature–time histories at the monitoring points, shown in Fig. 13, reveal that the temperatures of all models exhibit irregular oscillations, indicating that the system has deviated from single-frequency periodic dynamics. The heat transfer process is further intensified and evolves toward a multi-frequency coupled complex pattern.

Figure 13: Dimensionless temperature at the monitoring points of each model at the corresponding Rayleigh numbers.

Further analysis of the phase-space trajectories in Fig. 14 can more clearly reveal the dynamic evolution process of the system. For Model 1, the temperature curve superficially appears to retain periodic-like oscillations, but its phase trajectory forms a continuously contracting two-dimensional toroidal structure rather than collapsing to a single point, indicating that Model 1 has entered a quasi-periodic oscillatory state. For Models 2, 3, and 4, the phase-space trajectories are concentrated within closed band-like regions and gradually converge to a set of trajectories with a two-dimensional toroidal shape, which is a typical manifestation of quasi-periodic attractor in phase space. Therefore, under the corresponding Rayleigh number conditions, all four models have transitioned from periodic oscillations to quasi-periodic oscillations, with the system’s dynamical characteristics evolving from single-frequency to higher-order multi-frequency coupled behavior as the Rayleigh number increases.

Figure 14: Phase portraits in the u–v plane at the monitoring points of each model at the corresponding Rayleigh numbers: (a) Model 1; (b) Model 2; (c) Model 3; (d) Model 4.

4.4 Chaotic Solution

With a further increase in the Rayleigh number, the flow structures and heat transfer characteristics of the system become highly complex. Fig. 15 presents the temperature and flow fields for Model 1 at Ra = 2.5 × 106 and Model 2 at Ra = 5.5 × 105, where panels (a) and (c) show the temperature fields, and panels (b) and (d) depict the corresponding flow fields. It can be observed that the convective structures within the annular space have completely lost symmetry, with the flow patterns becoming strongly distorted and highly unsteady, indicating that the system has departed far from the regular dynamics of periodic or quasi-periodic states.

Figure 15: Temperature and flow fields of each model at the corresponding Rayleigh numbers: (a,c) temperature fields; (b,d) flow fields.

The dimensionless temperature–time curves at the monitoring points, shown in Fig. 16, reveal pronounced unpredictability and irregular fluctuations, indicating a progression toward fully unpredictable behavior.

Figure 16: Dimensionless temperature at the monitoring points of each model at the corresponding Rayleigh numbers.

To further elucidate the dynamical mechanisms, the phase-space trajectories in Fig. 17 show extreme instability, with neighboring trajectories repelling and diverging from each other. Nevertheless, overall trajectories outside the attractor still converge toward it, indicating that the system’s state remains bounded. This behavior exhibits the typical characteristics of chaos, with the solution structure becoming increasingly complex and Hopf bifurcations occurring again. The attractor at this stage is referred to as a chaotic attractor. The calculation of the Lyapunov exponent at this Rayleigh number confirms that the system has entered a chaotic state, and the results of the runs test provide statistically robust support for this conclusion. In Fig. 18, under chaotic conditions, the average Nusselt number on the inner wall surface exhibits irregular and aperiodic violent fluctuations, further confirming that the system has entered a chaotic thermal convection state.

Figure 17: Phase portraits in the u–v plane at the monitoring points of each model at the corresponding Rayleigh numbers: (a) Model 1; (b) Model 2.

Figure 18: Temporal evolution of Nu¯ in different models.

4.5 Stability Analysis of Solutions in Comparison with an Eccentric Annulus

Fig. 19 illustrates the distributions of critical Rayleigh numbers for the four models under steady, periodic, quasi-periodic, and chaotic states. In the figure, the color shading from light to dark corresponds to steady, periodic, quasi-periodic, and chaotic solutions, respectively. It can be observed that, for a given dynamical state (e.g., steady or periodic), the critical Rayleigh number decreases significantly as the radius ratio decreases. That is, models with smaller radius ratios are more prone to transition to unstable states or more complex oscillatory behaviors under lower driving forces.

Figure 19: Critical Rayleigh number diagrams for different radius ratio models in corresponding states.

This trend can be explained based on the system’s geometric characteristics: when the radius ratio decreases, the annular space becomes narrower, concentrating the buoyancy-driven forces. The convection paths are shortened, facilitating the accumulation of temperature and velocity gradients, which more rapidly disrupts the symmetric structure of the system and induces instabilities. Therefore, for the same solution type, models with smaller radius ratios correspond to lower critical Rayleigh numbers, whereas models with larger radius ratios exhibit stronger geometric stability and require higher Rayleigh numbers to trigger the same flow states.

This trend is consistently observed across the four typical dynamical states—steady, periodic, quasi-periodic, and chaotic—highlighting the significant influence of the radius ratio on the bifurcation pathways and instability evolution of natural convection in annular spaces.

Table 3 compares the critical Rayleigh numbers of Model 1 in this study with those of an eccentric annulus model with the same radius ratio reported in [21] under different solution types. The eccentric annulus in [21] is defined by the geometric parameters ε=(r,φ), where the eccentricity is given by r = L/(Ro − Ri) = 0.625 and the deflection angle is φ = 180°, with L representing the distance between the centers of the inner and outer circles.

From the data in Table 3, it is evident that for all typical dynamical states—steady, periodic, quasi-periodic, and chaotic—the critical Rayleigh numbers of the concentric annulus are significantly higher than those of the eccentric annulus. This difference primarily arises from the effect of geometric symmetry. The concentric annulus possesses strict axial symmetry, resulting in relatively uniform temperature and flow fields. Thermal disturbances in natural convection are thus less likely to break this inherent symmetry, delaying the onset of instability and requiring higher Rayleigh numbers for bifurcations or transition to more complex oscillatory states.

In contrast, the eccentric annulus exhibits spatially non-uniform temperature and velocity gradients due to the offset of the inner circle, naturally breaking symmetry and making the flow more susceptible to instabilities. As a result, under the same dynamical state, the critical Rayleigh numbers of the concentric annulus are generally higher than those of the eccentric annulus.

5  Conclusions

This study systematically investigated the geometric sensitivity and nonlinear dynamic transition of natural convection within a toroidal cavity, providing new insights into the generation of chaos and the role of geometric parameters in flow stability. The Lattice Boltzmann method (LBM) was used to numerically simulate natural convection and heat transfer in concentric toroidal models with four different geometric dimensions. The instability of the system was assessed using the maximum Lyapunov exponent, run test, and phase space analysis. The main conclusions are summarized below:

(1)   With increasing Rayleigh number, the system successively undergoes transitions from a steady state to a periodic state, a quasi-periodic state, and a chaotic state, exhibiting complex nonlinear dynamic behavior.

(2)   The radius ratio has a significant impact on flow stability. A decrease in the radius ratio leads to a lower critical Rayleigh number, indicating increased sensitivity of the system to geometric constraints. A larger radius ratio enhances the geometric stability of the system, requiring a stronger driving force to trigger instability.

(3)   Compared to concentric rings, the temperature gradient in the eccentric annular region is non-uniform due to the offset of the inner circle. This makes it easier to induce flow instability and allows the system to achieve the same type of dynamic structure at a lower Rayleigh number.

Acknowledgement: Not applicable.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Haojie Ju and Ming Zhao; methodology, Haojie Ju and Ming Zhao; software, Haojie Ju; validation, Haojie Ju; formal analysis, Haojie Ju; investigation, Haojie Ju; resources, Ming Zhao; data curation, Haojie Ju; writing—original draft preparation, Haojie Ju; writing—review and editing, Haojie Ju and Ming Zhao; visualization, Haojie Ju; supervision, Ming Zhao; project administration, Ming Zhao; funding acquisition, Ming Zhao. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data available on request from the authors. The data that support the findings of this study are available from the corresponding author, [Ming Zhao], upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Chillà F, Schumacher J. New perspectives in turbulent Rayleigh-Bénard convection. Eur Phys J E. 2012;35(7):58. doi:10.1140/epje/i2012-12058-1. [Google Scholar] [PubMed] [CrossRef]

2. Bénard H. Les tourbillons cellulaires dans une nappe liquide. J Phys Theor Appl. 1901;10(1):254–66. [Google Scholar]

3. Rayleigh L. LIX. On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. Lond Edinb Dublin Philos Mag J Sci. 1916;32(192):529–46. doi:10.1080/14786441608635602. [Google Scholar] [CrossRef]

4. Huang SD, Xia KQ. Effects of geometric confinement in quasi-2-D turbulent Rayleigh-Bénard convection. J Fluid Mech. 2016;794:639–54. doi:10.1017/jfm.2016.181. [Google Scholar] [CrossRef]

5. Weiss S, Ahlers G. Effect of tilting on turbulent convection: cylindrical samples with aspect ratio. J Fluid Mech. 2013;715:314–34. doi:10.1017/jfm.2012.520. [Google Scholar] [CrossRef]

6. Yin C, Wang C, Wang S. Thermal instability of a viscoelastic fluid in a fluid-porous system with a plane Poiseuille flow. Appl Math Mech Engl Ed. 2020;41(11):1631–50. doi:10.1007/s10483-020-2663-7. [Google Scholar] [CrossRef]

7. Xu A, Chen X, Xi HD. Tristable flow states and reversal of the large-scale circulation in two-dimensional circular convection cells. J Fluid Mech. 2021;910:A33. doi:10.1017/jfm.2020.964. [Google Scholar] [CrossRef]

8. Schindler F, Eckert S, Zürner T, Schumacher J, Vogt T. Collapse of coherent large scale flow in strongly turbulent liquid metal convection. Phys Rev Lett. 2022;128(16):164501. doi:10.1103/physrevlett.128.164501. [Google Scholar] [PubMed] [CrossRef]

9. Hu H, Qin JJ, Zheng LY, Li XL, Zhao BX. Geometric confinement and fluid properties effects on stability and heat transfer of thermal convection in lateral heated cavities. Int J Heat Fluid Flow. 2026;119:110312. doi:10.1016/j.ijheatfluidflow.2026.110312. [Google Scholar] [CrossRef]

10. Batchelor GK. Heat transfer by free convection across a closed cavity between vertical boundaries at different temperatures. Q Appl Math. 1954;12(3):209–33. doi:10.1090/qam/64563. [Google Scholar] [CrossRef]

11. Pan X, Yu W, Choi JI. Effects of aspect ratio on Rayleigh-Bénard convection under non-Oberbeck-Boussinesq effects in glycerol. Eur Phys J Plus. 2023;138(12):1096. doi:10.1140/epjp/s13360-023-04672-0. [Google Scholar] [CrossRef]

12. Mac Huang J, Zhang J. Side-heated Rayleigh-Bénard convection. J Fluid Mech. 2024;999:A35. doi:10.1017/jfm.2024.968. [Google Scholar] [CrossRef]

13. Bishop EH, Carley CT, Powe RE. Natural convective oscillatory flow in cylindrical annuli. Int J Heat Mass Transf. 1968;11(12):1741–52. doi:10.1016/0017-9310(68)90017-3. [Google Scholar] [CrossRef]

14. Powe RE, Carley CT, Bishop EH. Free convective flow patterns in cylindrical annuli. J Heat Transf. 1969;91(3):310–4. doi:10.1115/1.3580158. [Google Scholar] [CrossRef]

15. Usman M, Son JH, Park IS. A low-Rayleigh transition into chaos for natural convection inside a horizontal annulus at Prandtl number 0.1. Int J Heat Mass Transf. 2021;179(2):121658. doi:10.1016/j.ijheatmasstransfer.2021.121658. [Google Scholar] [CrossRef]

16. Ahlers G, Bodenschatz E, Hartmann R, He X, Lohse D, Reiter P, et al. Aspect ratio dependence of heat transfer in a cylindrical Rayleigh-Bénard cell. Phys Rev Lett. 2022;128(8):084501. doi:10.1103/PhysRevLett.128.084501. [Google Scholar] [PubMed] [CrossRef]

17. Zhu X, Zhou Q. Flow structures of turbulent Rayleigh-Bénard convection in annular cells with aspect ratio one and larger. Acta Mech Sin. 2021;37(8):1291–8. doi:10.1007/s10409-021-01104-z. [Google Scholar] [CrossRef]

18. Wang D, Jiang H, Liu S, Zhu X, Sun C. Effects of radius ratio on annular centrifugal Rayleigh–Bénard convection. J Fluid Mech. 2022;930:A19. doi:10.1017/jfm.2021.889. [Google Scholar] [CrossRef]

19. Farkach Y, Derfoufi S, Ahachad M, Bahraoui F, Mahdaoui M. Numerical investigation of natural convection in concentric cylinder partially heated based on MRT-lattice Boltzmann method. Int Commun Heat Mass Transf. 2022;132:105856. doi:10.1016/j.icheatmasstransfer.2021.105856. [Google Scholar] [CrossRef]

20. Li Q, He YL, Wang Y, Tao WQ. Coupled double-distribution-function lattice Boltzmann method for the compressible Navier-Stokes equations. Phys Rev E Stat Nonlin Soft Matter Phys. 2007;76(5 Pt 2):056705. doi:10.1103/PhysRevE.76.056705. [Google Scholar] [PubMed] [CrossRef]

21. Xu YF, Zhang HM, Zhao M. Determinism and chaos of natural convection in eccentric annulus. Chin J Comput Phys. 2025;42(1):28–37. (In Chinese). doi:10.37649/aengs.2008.14206. [Google Scholar] [CrossRef]

22. Guo Z, Shi B, Zheng C. A coupled lattice BGK model for the Boussinesq equations. Numer Methods Fluids. 2002;39(4):325–42. doi:10.1002/fld.337. [Google Scholar] [CrossRef]

23. Li Z, Yang M, Zhang Y. Heat transfer characteristics of natural convection in an annulus: correlation of Nusselt number with Rayleigh number and geometric parameters. Int J Heat Mass Transf. 2016;94:222–31. [Google Scholar]

24. Kuehn TH, Goldstein RJ. An experimental study of natural convection heat transfer in concentric and eccentric horizontal cylindrical annuli. J Heat Transf. 1978;100(4):635–40. doi:10.1115/1.3450869. [Google Scholar] [CrossRef]