Boyouncy

1 Nonlinear dynamic analysis of solution multiplicity of
2
buoyancy ventilation in a typical underground
3 structure
4 Yanan Liua,b，Yimin Xiaoa,c*, Jianli Chend,
5 a School of Civil Engineering, Chongqing University, Chongqing 400045, China
6 b Key Laboratory of New Technology for Construction of Cities in Mountain Area of
7 Ministry of Education, Chongqing University, Chongqing 400045, China
8 c National Centre for International Research of Low-carbon and Green Buildings,
9 Chongqing University, Chongqing 400045, China
10 d National Renewable Energy Lab, USA

11
12 *Corresponding email: [email protected]
13 Postal address: School of Civil Engineering, Chongqing University, Chongqing, 400045, China
14 Abstract：
15 Buoyancy ventilation is widely used in underground buildings, such as underground
16 hydropower stations. Multiple solutions of buoyancy ventilation may exist in those
17 underground structures. In this study, we developed a transient model comprising an
18 ordinary differential equation system to describe buoyancy ventilation patterns in
19 typical two-zone underground structures. Additionally, the accuracy of the model was
20 validated. Nonlinear dynamical analysis was conducted to study multiple steady-state
21 airflow. According to mathematical derivation, the configuration of one local heat
22 source at the bottom corner introduces two stable solutions. The criterion to determine
23 the stability and existence of solutions for more general scenarios was developed. Using
24 this criterion, we obtained the multiple steady states of any two-zone underground
25 buildings for different stack height ratios and the strength ratios of the heat sources.
1
Pursuant to the DOE Public Access Plan, this document represents the authors' peer-reviewed, accepted manuscript.
The published version of the article is available from the relevant publisher.
1 This criterion can be adopted for the design of buoyancy ventilation or natural smoke
2 ventilation systems. Designers can change the height ratio of the stack or the heat ratio
3 of two zones to induce the desired ventilation patterns. Finally, a case study was
4 conducted with field measurements to demonstrate the use of the nonlinear dynamical
5 analysis method to investigate the multiple steady states of buoyancy ventilation.
6 Through the case study, we validated that the proposed criterion could produce the same
7 result as the nonlinear dynamical analysis.
8 Keywords：Nonlinear dynamics; buoyancy ventilation; multiple steady states；

9 Underground buildings
10 Nomenclature
𝐴1 Coefficient matrix for linearized differential equation system
𝑞1 Mass flow rate at zone 1, [kg/s]
𝑞2 Mass flow rate at zone 2, [kg/s]
𝑇1 Air temperature at zone 1, [K]
𝑇𝑎 Outdoor air temperature, [K]
𝑆1 Coefficient of Mass flow impedance at zone1
𝑆2 Coefficient of Mass flow impedance at zone2
𝑀1 Thermal mass at zone 1, [kg]
𝑀2 Thermal mass at zone 2, [kg]
𝐸1 Total heat gain at zone 1 , [kW]
𝐸2 Total heat gain at zone 2 , [kW]

𝐶𝑝 Specific heat of air,[kJ/(kg  K)]
t Time, [s]
𝑔 Gravitational acceleration, [m/s2]
𝐻1 Height of zone 1,[m]
𝐻2 Height of zone 2,[m]

𝐸2
𝜅 Heat ratio between two zones, 𝐸
1
∆𝑇1 Temperature difference between indoor and outdoor air at zone 1, [℃]
∆𝑇2 Temperature difference between indoor and outdoor air at zone 2, [℃]
𝐻2
𝛼 Height ratio between two zones, 𝐻
1
̅̅̅̅̅1
∆𝑇 Temperature difference between indoor and outdoor air at zone 1 in steady
state, [℃]
̅̅̅̅̅
∆𝑇2 Temperature difference between indoor and outdoor air at zone 1 in steady
state, [℃]
𝜌0 Ambient air density, [kg/m3]
𝜆 Eigenvalue of coefficient matrix

1
2 1. Introduction
3 Owing to the wide utilization of underground spaces, the evaluation of natural
4 ventilation in underground buildings has become significant [1-11]. Natural ventilation
5 as a passive strategy has been adopted in many underground structures, such as
6 underground shelters [1], warehouses [2],underground garages [3, 4], mines [5],
7 hydropower stations [8, 9], and underground roads and railway tunnels [10]. Aspects
8 such as air quality [3, 4], energy conservations [6, 7], computational models [5, 8, 9],
9 fire safety [5, 10], and ventilation performance [1, 2] have been investigated. Hence, as
10 a type of natural ventilation, buoyancy ventilation in underground structures is worth
11 investigating.
12 Buoyancy ventilation is the air flow driven by the air density difference caused by
13 air temperature difference. As early as 1954, Batchelor [12] began investigating the
14 buoyancy of airflow. Many factors, such as envelope heat transfer or internal heat
15 source, can result in thermal ventilation. Among them, various studies have been
16 performed regarding solar-chimney-induced thermal ventilation [13, 14], natural
17 ventilation reinforced by double-skin facade [15-20], indoor thermal plume, and single-
3
1 sided [21] or cross ventilation [22, 23] driven by the separate or combined effect of
2 buoyancy and wind pressure.
3 For the ventilation driven by the combined effect of buoyancy and wind pressure,
4 multiple solutions may exist. Hunt G.R & Linden.P.F [24] first raised the issue of
5 mutual reinforcement and confrontation of natural ventilation under wind pressure and
6 thermal pressure. However, they primarily studied the situation with the reinforced
7 effect of wind pressure and thermal pressure. A one-dimensional mathematical model
8 was proposed with a visual and quantitative comparative study performed using a
9 small-scaled brine model experiment. This study provides a theoretical basis and guide
10 for the calculation of natural ventilation of night cooling or air purging system of gas
11 leakage. As the first study that introduced the (concept or possibility) of multiple
12 solutions of building ventilation and smoke exhaust, Nitta [25] demonstrated that
13 solution multiplicity exists under specific room layouts and fan settings in the design
14 of smoke prevention and exhaust. This study presents the importance of solution
15 multiplicity to personnel safety and ventilation design, although its formation
16 mechanism is not elaborated.
17 Subsequently, the solution multiplicity of single-zone and double-opening buildings
18 under the combined effect of wind pressure and thermal pressure has attracted wide
19 attention. The existence of solution multiplicity was first investigated, where multiple
20 methods were employed to reproduce this phenomenon. Heiselberg et al. [26] analyzed
21 the multiple steady sates of the single-zone and double-opening buildings with wind
22 pressure and buoyancy confrontation by a salt water experiment and CFD simulation.
23 Based on this typical building configuration, Li & Delsante [27] established a complete
24 one-dimensional mathematical model while considering the effect of heat transfer in
25 the building envelope. Additionally, Li [28] reported that solution multiplicity could
26 exist in both inclined tunnels and two-story aboveground buildings, although the focus
27 of the study was wind pressure and buoyancy confrontation in a single-zone building.
28 Subsequently, the dynamical process and impact factors of multiple steady states were
1 investigated. Lishman & Woods [29] studied the effect of wind pressure changes on the
2 transition of multiple steady states in a single-zone building. Yuan & Glicksman [30-
3 32] studied the effects of different initial conditions on the formation of multiple steady
4 states in single-zone buildings under the combined effect of wind pressure and
5 buoyancy. The dynamic transition between different steady states was investigated
6 considering the effects of disturbance magnitude and action time. Gladstone et al. [33]
7 studied a single-zone building with distributed heat sources on the floor. The effects of
8 distributed roof cold source, outdoor temperature gradient, and outdoor wind pressure
9 on multiple steady states were discussed. In the study, they proposed a one-dimensional
10 model and compared it with experiments. Erhan & Hifzi [34] used IEA Annex 20 as an
11 example and discovered that different turbulence parameters may produce multiple
12 solutions in CFD simulations.
13 In addition, other studies focusing on the solution multiplicity of ventilation in single-
14 zone and multiopening buildings have been performed. Chenvidakarn & Woods [35]
15 and Durrani et al. [36] analyzed the solution multiplicity of a typical aboveground
16 building through a saltwater experiment employing one-dimensional model analyses
17 [35] and CFD simulations [36].The building contained one zone and three openings.
18 Two upper openings were the chimneys, while a lower large opening was the entrance
19 door. However, these studies primarily focused on the accuracy comparison of the LES
20 model with K-e model in CFD. The analytical model was established based on an
21 equilibrium state without considering the effect of thermal mass and time; therefore, an
22 in-depth analysis of the stability of multiple solutions was not performed. Chen & Li
23 [37] studied the buoyancy ventilation of a single-zone building with three horizontal
24 openings in different levels using theoretical analysis. For a specific geometric structure,
25 even with the same boundary conditions and geometric settings, the height of thermal
26 stratification might be higher or lower than the medium-level horizontal opening when
27 the initial conditions are different.
28 In a fire, the solution multiplicity of smoke flows is a significant topic that can guide
1 the design of smoke exhaust systems, thereby ensuring the safety of human evacuation.
2 Gong J. and Li Y. [38, 39] studied the solution multiplicity of smoke spread with the
3 effect of outdoor wind pressure in fires. The research involved both small-scaled
4 experiments and CFD simulations, including a single heat source in a typical single-
5 zone building and a single heat source in a two-zone building. A comparative study was
6 conducted using different forms of heat sources, such as point, line, and area heat
7 sources. In the study, they investigated the effects of different heat source locations on
8 solution multiplicity and presented a visualization experiment. Yang D. [40] analyzed
9 the confrontation between outdoor wind pressure and thermal pressure caused by fire
10 in an oblique straight tunnel. The one-dimensional model, which was based on the
11 transient energy balance and pressure balance of a single-zone building to establish a
12 nonlinear differential equation, was similar to that of Yuan & Glicksman [31], despite
13 the difference in geometry. Furthermore, the application scenario of the study was
14 different, in that it was suitable for the fire situation in the oblique straight tunnel instead
15 of a single-zone building. Furthermore, salt water experiments were conducted to
16 compare the entire process of formation and development of the two steady-state
17 solutions.
18 Additionally, studies regarding the solution multiplicity of ventilation in two-zone
19 spaces have been performed. Yang L. [41, 42] performed a detailed analysis of multiple
20 steady states and bifurcation of fluids in a two-zone building with four openings using
21 theoretical analyses and CFD simulations. Li [43] et al. investigated the buoyancy
22 ventilation in a two-story space with two heat sources and three openings. A
23 mathematical model was established using the nonlinear ordinary differential equation.
24 The effect of the heat source’s strength ratio on fluid bifurcation was analyzed. Yang D.
25 [44] analyzed the smoke exhaust spread in a tunnel with three entrances in fire scenarios
26 and concluded that six equilibrium states might exist. Subsequently, based on the energy
27 balance and pressure balance equations of the steady states, a mathematical model was
28 established to solve the smoke exhaust of each tunnel. CFD was also used to simulate
1 and reproduce some of the steady states. However, in the proposed one-dimensional
2 model that was based on equilibrium conditions, the transient development process of
3 smoke ventilation was not considered. Additionally, the stability of the solution was not
4 analyzed. Liu et al. [45] numerically studied the formation process of multiple steady
5 states in an underground building with two tunnel connecting to the outdoor
6 environment. The tunnels were set with equal heights with only one heat source at the
7 corner of the deep buried spaces. The authors used the two-stage CFD method to
8 reproduce the two steady states of the buoyancy ventilation in the underground building
9 by changing the initial conditions. This study primarily provided a method to
10 investigate the multiple steady states.

11 In summary, despite the abovementioned studies, gaps still exist in predicting and
12 analyzing the multiplicity of buoyancy ventilation in underground buildings. First,
13 current studies are still limited to scenarios with single zones or two zones driven by
14 combined wind and buoyancy. More specifically, the building configurations and
15 driving forces are different. For the building configurations, most studies focus on
16 single-zone buildings with two openings, which differ from underground structures.
17 Typically, at least two tunnels are connected with the outdoor environment. Hence, at
18 least two zones exist for the underground structures provided that the height of the deep
19 buried rooms are neglected. For the driven forces, the combat between thermal
20 buoyancy and wind pressure is the main cause of solution multiplicity in previous
21 studies. However, underground buildings are not exposed to the outdoor environment,
22 and the wind pressure is not highly significant. By contrast, the heat transfer between
23 the indoor air and surrounding soil can affect the thermal pressure inside the tunnel.
24 Therefore, both zones may need to be considered as heat sources/sinks. The charge and
25 discharge process of heat from the soil to the tunnel can contribute to the solution
26 multiplicity of buoyancy ventilation in underground buildings. Additionally, heat is
27 released from the indoor environment, such as from equipment and human body;
28 therefore, the solution multiplicity of natural ventilation is driven by the combat of
1 thermal pressure between two tunnels. Next, the current studies are not generalizable.
2 To study multiple steady states, the CFD method or analytical method must be repeated
3 each time the geometry is changed. Especially for the CFD method, the initial
4 conditions must be changed and numerous CFD simulations must be conducted to study
5 the existence of multiple steady states for a single fixed geometry. This requires
6 significant computational resources and manpower. For the analytical method, the
7 strength of the buoyancy and the wind pressure are the main control parameters studied
8 to investigate solution multiplicity. In other words, buoyancy and wind pressure were
9 altered to investigate their effects on the performances of multiple steady states. This
10 may be useful for aboveground buildings if the geometry is fixed. However, for a more
11 general case, the height ratio of the tunnel may be an important control parameter.
12 In our study, we considered both the effects of the strength ratio of heat sources and
13 the different height ratios of tunnels simultaneously. To attain a deep understanding of
14 the solution multiplicity of underground buoyancy ventilation, we performed a
15 nonlinear dynamical analysis to study the formation mechanism of multiple steady
16 states. In addition to demonstrating the use of nonlinear dynamics to analyze the
17 buoyancy ventilation in underground buildings, our goal is to develop a criterion for
18 evaluating the multiple steady states for different buildings and heat source
19 configurations. This criterion is based on the strength ratio of the heat sources between
20 two zones and the tunnel height ratio of two zones. Once these two parameters are
21 determined, then whether a solution multiplicity exists for a typical underground
22 structure with two tunnels connected to the outdoor environment can be determined.
23 This is a straightforward method for the design and optimization of buoyancy
24 ventilation and smoke ventilation in underground buildings.
25 The organization of this paper is as follows: Section 1 presents the literature review
26 of previous studies regarding buoyancy ventilation and solution multiplicity with an
27 emphasis on the significance of our study; Section 2 provides the nonlinear dynamic
28 analysis of a typical underground structure, which includes the establishment of a one-
1 dimensional model, stability and existence analysis of underground buoyancy
2 ventilation, model validation through result comparison with previous literature, and
3 graphical presentation of multiple steady states through bifurcation diagram and phase
4 portrait; Section 3 presents a case study of an underground hydropower station to
5 demonstrate how nonlinear dynamical analysis is employed in real project applications.
6 Furthermore, the derived criterion is compared with the analysis results; Section 4
7 presents the conclusions, and current studies are summarized along with future studies.
8 2. Nonlinear dynamic analysis
9 To study the nonlinear dynamics of typical deep-buried underground buildings with two
10 openings, we first established the transient mathematical mode. Some assumptions
11 were made: (1) Each zone was well mixed;（2）Thermal mass was 1;（3）𝐸1 > 0; (4)
12 The mass flow impedance coefficient of the geometry was constant. The conservation
13 law was adopted to develop the transient model. As the driven force of the ventilation
14 is the thermal pressure in two tunnels, as illustrated in Fig. 1, we divided the building
15 into two zones. All heat transfer, including the envelope heat transfer and internal heat
16 source, was considered as one heat source in each zone. For analysis convenience, we
17 assumed the heat source in the left tunnel as positive, and the heat source of the right
18 tunnel could be either negative or positive such that we could discuss the scenario of
19 two opposing heat sources and the scenario of one heat sink and one heat source. As
20 shown in Fig. 1, there would be two realizations for the building configurations.
1 2.1 Description of mathematical model
2
3 Fig. 1. Schematics of a typical two-zone underground structure.
4 For realization 1, the conservation of mass is as follows:
5 𝑞1 = 𝑞2 （2-1）
6 In realization 1, the air flow enters from zone 1 to zone 2. Hence, the thermal pressure
7 in zone 1 will resist the airflow, while the thermal pressure in zone 2 will assist the air
8 flow. Based on the looped method, the sum of the buoyancy pressure balances the flow-
9 element pressure losses [8, 46].The pressure balance equation is as follows:

𝑇1 −𝑇𝑎 𝑇2 −𝑇𝑎
10 − 𝜌𝑎 𝑔𝐻1 + 𝜌𝑎 𝑔𝐻2 = 𝑆1 𝑞12 + 𝑆2 𝑞22 （2-2）
𝑇𝑎 𝑇𝑎
11 The heat gain of the internal thermal mass is equal to the heat released by the heat
12 sources minus the heat loss through airflows. The heat balance equation for zones 1 and
13 2 can be obtained:
𝑑𝑇1
14 𝑀1 𝐶𝑝 = −𝑞1 𝐶𝑝 (𝑇1 −𝑇𝑎 ) + 𝐸1 （2-3）
𝑑𝑡
𝑑𝑇2
15 𝑀2 𝐶𝑝 = −𝑞2 𝐶𝑝 (𝑇2 −𝑇1 ) + 𝐸2 （2-4）
𝑑𝑡
16 From Eqs. 2-1 and 2-2, we can obtain
10
𝑇 −𝑇 𝑇 −𝑇
− 1𝑇 𝑎 𝜌𝑎 𝑔𝐻1 + 2𝑇 𝑎 𝜌𝑎 𝑔𝐻2
1 𝑞1 = 𝑞2 = √ 𝑎 𝑎
（2-5）
𝑆1 +𝑆2
2 Substituting Eq. 2.5 into Eqs. 2-3 and 2-4 results in the following two-dimensional
3 nonlinear ordinary differential equation system:
𝑑𝑇1 − 1 𝑎 𝜌𝑎 𝑔𝐻1 + 2 𝑎 𝜌𝑎 𝑔𝐻2
= −√ 𝑇𝑎 𝑇𝑎
4 𝑀1 𝐶𝑝 𝐶𝑝 (𝑇1 −𝑇𝑎 ) + 𝐸1 （2-6）
𝑑𝑡 𝑆1+𝑆2
𝑑𝑇2 − 1 𝑎 𝜌𝑎 𝑔𝐻1 + 2 𝑎 𝜌𝑎 𝑔𝐻2
= −√ 𝑇𝑎 𝑇𝑎
5 𝑀2 𝐶𝑝 𝑑𝑡 𝑆1+𝑆2
𝐶𝑝 (𝑇2 −𝑇1) + 𝐸2 （2-7）
6 The prerequisite for this equation system is− 𝜌𝑎 𝑔𝐻1 + 𝜌𝑎 𝑔𝐻2 > 0, which
𝑇𝑎 𝑇𝑎
7 means that the thermal pressure in the right tunnel should be greater than that in the left
8 tunnel.
9 Similarly, we can obtain the mass balance equation for realization 2:
10 𝑞1 = 𝑞2 （2-8）
11 In realization 2, the air flow enters from zone 2 to zone 1. Hence, the thermal pressure
12 in zone 2 will resist the airflow, while the thermal pressure in zone 1 will assist the air
13 flow. We can obtain the pressure balance equation:

14 𝜌𝑎 𝑔𝐻1 − 𝜌𝑎 𝑔𝐻2 = 𝑆1 𝑞12 + 𝑆2 𝑞22 （2-9）
𝑇𝑎 𝑇𝑎
15 Furthermore, the heat balance equation for realization 2 is as follows:

𝑑𝑇1
16 𝑀1 𝐶𝑝 = −𝑞1 𝐶𝑝 (𝑇1 −𝑇2 ) + 𝐸1 （2-10）
𝑑𝑡
𝑑𝑇2
17 𝑀2 𝐶𝑝 = −𝑞2 𝐶𝑝 (𝑇2 −𝑇𝑎 ) + 𝐸2 （2-11）
𝑑𝑡
18 From Eqs. 2-8 and 2-9, we can obtain
𝑇1 −𝑇𝑎 𝑇 −𝑇
𝜌𝑎 𝑔𝐻1 − 2𝑇 𝑎 𝜌𝑎 𝑔𝐻2
𝑞1 = 𝑞2 = √ 𝑇𝑎
19 𝑎
（2-12）
𝑆1+𝑆2
20 Substituting Eq.2.12 into Eqs. 2-10 and 2-11 results in the following two-dimensional
21 nonlinear ordinary differential equation system:
𝑑𝑇1 𝜌𝑎 𝑔𝐻1 − 2𝑇 𝑎 𝜌𝑎 𝑔𝐻2
= −√ 𝑇𝑎
22 𝑀1 𝐶𝑝 𝑎
𝐶𝑝 (𝑇1 −𝑇2 ) + 𝐸1 （2-13）
𝑑𝑡 𝑆1 +𝑆2
11
𝑑𝑇2 𝜌𝑎 𝑔𝐻1 − 2𝑇 𝑎 𝜌𝑎 𝑔𝐻2
= −√ 𝑇𝑎
1 𝑀2 𝐶𝑝 𝑎
𝐶𝑝 (𝑇2 −𝑇𝑎 )+𝐸2 （2-14）
𝑑𝑡 𝑆1 +𝑆2
2 The prerequisite for this equation system is − 𝜌𝑎 𝑔𝐻1 + 𝜌𝑎 𝑔𝐻2 < 0, which
𝑇𝑎 𝑇𝑎
3 means that the thermal pressure in the left tunnel should be greater than that in the right
4 tunnel.
5 2.2 Stability and existence of the system
𝜌𝑎 𝑔𝐻1
𝐸2
6 Assuming 𝜅 = , ∆𝑇1 = 𝑇1 −𝑇𝑎 , ∆𝑇2 = 𝑇2 −𝑇𝑎 , 𝑛 = √ 𝑇𝑎
, 𝛼 = 𝐻2 /𝐻1 , 𝐸1 >
𝐸1 𝑆1+𝑆2
7 0, 𝐶𝑝 = 1, 𝑀1 = 𝑀2 = 1 , the nonlinear ordinary differential equation system can be
8 further simplified.
9 2.2.1 Stability analysis for scenario 1(𝜅 is fixed, 𝛼 is control parameter)
10 2.2.1.1 Stability analysis for 𝜅 > 0
11 First, we begin from the scenario where two heat sources are positive. For status 1,
𝑑∆𝑇1
12 𝑓1 (∆𝑇1 , ∆𝑇2 ) = = −𝑛√𝛼∆𝑇2 − ∆𝑇1 ∆𝑇1 + 𝐸1 （2-15）
𝑑𝑡
13
𝑑∆𝑇2
14 𝑓2 (∆𝑇1 , ∆𝑇2 ) = = −𝑛√𝛼∆𝑇2 − ∆𝑇1(∆𝑇2 − ∆𝑇1 ) + 𝐸2 （2-16）
𝑑𝑡
15 For realization 2, Eqs. 2-13 and 2-14 can be simplified as follows:

𝑑∆𝑇1
16 𝑓3 (∆𝑇1 , ∆𝑇2 ) = = −𝑛√∆𝑇1 − 𝛼∆𝑇2(∆𝑇1 − ∆𝑇2 ) + 𝐸1 （2-17）
𝑑𝑡
𝑑∆𝑇2
17 𝑓4 (∆𝑇1, ∆𝑇2 ) = = −𝑛√∆𝑇1 − 𝛼∆𝑇2 ∆𝑇2 + 𝐸2 （2-18）
𝑑𝑡
18
19 In summary，for the scenario with two positive heat sources (𝜅 > 0), when 0 <
1 1 5
20 𝛼< , the system has one stable fixed point in realization 2; when <𝛼< ,
1+𝜅 1+𝜅 4+5𝜅
21 the system has one unstable fixed point in realization 1 and a stable fixed point in
12
5 5+4𝜅
1 realization 2; when <𝛼< , the system has two stable fixed points; when
4+5𝜅 5𝜅
5+4𝜅 1+𝜅
2 <𝛼< , the system has one stable fixed point in realization 1 and an unstable
5𝜅 𝜅
1+𝜅
3 fixed point in realization 2; when < 𝛼, the system has one stable fixed point in
𝜅
4 realization 1 and no solution for realization 2. The detailed derivations are provided in
5 Appendices A1 & A2.
6 2.2.1.2 Stability analysis for 𝜅 < 0
7 When one heat source and one heat sink exist (𝜅 < 0), we still can utilize the
8 characteristic equation to evaluate the stability and existence of the nonlinear
9 differential equation system.
10 Assuming that no fixed point exists in realization 1, the following expression should be
11 complied:
−1 + 𝛼 + 𝛼𝜅 < 0
12 { 𝜅<0 （2-19）
𝛼>0
1
13 Therefore, 𝜅 ≤ −1 𝑎𝑛𝑑 𝛼 > 0 𝑜𝑟 − 1 < 𝜅 < 0 𝑎𝑛𝑑 0 < 𝛼 < 1+𝜅.
14 Assuming that a stable fixed point exists in realization 1, the following expression
15 should be complied:
−1 + 𝛼 + 𝛼𝜅 > 0
𝜅<0
16 { （2-20）
𝛼>0
−5 + 𝛼(4 + 5𝜅) > 0
4 5
17 Therefore, − 5 < 𝜅 < 0 and 𝛼 > 4+5𝜅.
18 Assuming that an unstable fixed point exists in realization 1, the following expression
−1 + 𝛼 + 𝛼𝜅 > 0
𝜅<0
20 { （2-21）
𝛼>0
−5 + 𝛼(4 + 5𝜅) < 0
4 1 4 1 5
21 Therefore, −1 < 𝜅 < − 5 and 𝛼 > 1+𝜅, or − 5 < 𝜅 < 0 and < 𝛼 < 4+5𝜅.
1+𝜅
13
2 complied:
1 + 𝜅 − 𝛼𝜅 < 0
3 { 𝜅<0 （2-22）
𝛼>0
1+𝜅
4 Therefore, 𝜅 < −1 and 0 < 𝛼 < .
𝜅
1 + 𝜅 − 𝛼𝜅 > 0
𝜅<0
7 { （2-23）
5 + (4 − 5𝛼)𝜅 > 0
𝛼>0
1+𝜅
8 Therefore, 𝜅 ≤ −1 and 𝛼 > , or −1 < 𝜅 < 0 and 𝛼 > 0.
𝜅
1 + 𝜅 − 𝛼𝜅 > 0
𝜅<0
11 { （2-24）
5 + (4 − 5𝛼)𝜅 < 0
𝛼>0
12 However, this expression system is not true.
13 In summary, for the scenario of one heat source and one heat sink (𝜅 < 0), different
4 4
14 situations should be considered: 𝜅 < −1, −1 < 𝜅 < − 5, and − 5 < 𝜅 < 0. Once the
15 interval of κ is fixed, we can evaluate the stability and existence of the fixed point
16 according to the value of α shown in Table 1.

17 Table 1
18 Criterion for scenario 1
𝜅 𝛼 Existence and stability Existence and stability for
for realization 1 realization 2

1 No Stable
(0,+∞) (0, )
1+𝜅
1 5 Unstable Stable
(1+𝜅 , 4+5𝜅)
14
5 5+4𝜅 Stable Stable
(4+5𝜅 , )
5𝜅
5+4𝜅 1+𝜅 Stable Unstable

( , )
5𝜅 𝜅
1+𝜅 Stable No
( , +∞)
𝜅
4 1 No Stable
(− 5 , 0) (0, 1+𝜅)
1 5 Unstable Stable
(1+𝜅 , 4+5𝜅)
5 Stable Stable
(4+5𝜅 , +∞)
4 1 No Stable
(−1, − ) (0, )
5 1+𝜅
1 Unstable Stable
(1+𝜅 , +∞)
1+𝜅 No No
(−∞, −1) (0, 𝜅
)
1+𝜅 No Stable
( , +∞)
𝜅
1 2.2.2 Stability analysis for scenario 2 (𝛼 is fixed, 𝜅 is control parameter)
2 In this scenario, the characteristic equation is the same as that of scenario 1.
3 Assuming no fixed point exists in realization 1, the following expression should be
4 complied:
−1 + 𝛼 + 𝛼𝜅 < 0
5 { （2-25）
𝛼>0
1−𝛼
6 Therefore, 𝛼 > 0 and 𝜅 < 𝛼
.
−1 + 𝛼 + 𝛼𝜅 > 0
9 { 𝛼>0 （2-26）
−5 + 𝛼(4 + 5𝜅) > 0
5−4𝛼
10 Therefore, 𝛼 > 0 and 𝜅 > .
5𝛼
15
−1 + 𝛼 + 𝛼𝜅 > 0
1 { 𝛼>0 （2-27）
−5 + 𝛼(4 + 5𝜅) < 0
1−𝛼 5−4𝛼
2 Therefore, 𝛼 > 0 and <𝜅< .
𝛼 5𝛼
4 complied:
1 + 𝜅 − 𝛼𝜅 < 0
5 { （2-28）
𝛼>0
1 1
6 Therefore, 0 < 𝛼 < 1 and 𝜅 < −1+𝛼 , or 𝛼 > 1 and 𝜅 > −1+𝛼.
1 + 𝜅 − 𝛼𝜅 > 0
9 {5 + (4 − 5𝛼)𝜅 > 0 （2-29）
𝛼>0
4 1 4 1 5
10 Therefore, 0 < 𝛼 ≤ 5 𝑎𝑛𝑑 𝜅 > −1+𝛼 or < 𝛼 < 1 𝑎𝑛𝑑 < 𝜅 < −4+5𝛼 ， or ≥
5 −1+𝛼
5
11 1 𝑎𝑛𝑑 𝜅 < −4+5𝛼 .
1 + 𝜅 − 𝛼𝜅 > 0
14 {5 + (4 − 5𝛼)𝜅 < 0 （2-30）
𝛼>0
4 5 5 1
15 Therefore, < 𝛼 ≤ 1 𝑎𝑛𝑑 𝜅 > −4+5𝛼, or 𝛼 > 1 𝑎𝑛𝑑 < 𝜅 < −1+𝛼.
5 −4+5𝛼
16 In summary, for scenario 2, different situations should be considered, i.e., 0 <

4 4
17 𝛼≤ , < 𝛼 < 1 , 𝛼 = 1 , and 𝛼 > 1 . Once the interval of 𝛼 is fixed, we can
5 5
18 evaluate the stability and existence of the fixed point according to the value of 𝜅
19 shown in Table 2.
20
21 Table 2
22 Criterion for scenario 2
𝛼 𝜅 Existence and stability Existence and stability for

16
for realization 1 realization 2
4 1 No No
(0, 5) (−∞, −1+𝛼)
1 1−𝛼 No Stable
(−1+𝛼 , )
𝛼
1−𝛼 5−4𝛼 Unstable Stable

( , )
𝛼 5𝛼
5−4𝛼 Stable Stable

( , +∞)
5𝛼
4 1 No No
(5 , 1) (−∞, −1+𝛼)
1 1−𝛼 No Stable
(−1+𝛼 , )
𝛼

( , )
𝛼 5𝛼
5−4𝛼 5 Stable Stable

( , −4+5𝛼)
5𝛼
5 Stable Unstable
(−4+5𝛼 , +∞)
1 (−∞, 0) No Stable
(0,0.2) Unstable Stable
(0.2,5) Stable Stable
(5, +∞) Stable Unstable
1−𝛼 No Stable
(1, +∞) (−∞, )
𝛼

( , )
𝛼 5𝛼
5−4𝛼 5 Stable Stable

( , )
5𝛼 −4+5𝛼
5 1 Stable Unstable
(−4+5𝛼 , −1+𝛼)
1 Stable No
(−1+𝛼 , +∞)
17
1 2.2.3 Stability analysis for scenario 3 (one heat source at the bottom of the building)
2
3 Fig. 2. Schematics of a typical two-zone underground structure with one local heat source.
5 As illustrated in Fig. 2，for status 1, the heat at the bottom releases to the right tunnel;
6 therefore,
𝑑∆𝑇1
7 𝑓1 (∆𝑇1 , ∆𝑇2 ) = = −𝑛√𝛼∆𝑇2 − ∆𝑇1 ∆𝑇1 （2-31）
𝑑𝑡
𝑑∆𝑇2
8 𝑓2 (∆𝑇1 , ∆𝑇2 ) = = −𝑛√𝛼∆𝑇2 − ∆𝑇1(∆𝑇2 − ∆𝑇1 ) + 𝐸1 （2-32）
𝑑𝑡
9 For status 2, the heat at the bottom releases to zone 1; therefore,

𝑑∆𝑇1
10 𝑓3 (∆𝑇1 , ∆𝑇2 ) = 𝑑𝑡
= −𝑛√∆𝑇1 − 𝛼∆𝑇2(∆𝑇1 − ∆𝑇2 ) + 𝐸1 （2-33）
𝑑∆𝑇2
11 𝑓4 (∆𝑇1, ∆𝑇2 ) = 𝑑𝑡
= −𝑛√∆𝑇1 − 𝛼∆𝑇2 ∆𝑇2 （2-34）
12 In summary, for the scenario of one heat source at the bottom of the building,
13 because the heat can enter either the left or the right tunnel, two steady states always
14 exist for this scenario. The parameter 𝛼 will not affect the stability and existence of
15 this building configuration; the detailed derivations are provided in Appendices B1 &
16 B2.
18
1 2.3 Model validation
2 To validate the model, we selected scenario 3 to compare the modeled results with
3 results from a previous study [44] because only this building configuration (one heat
4 source at the bottom with two adiabatic tunnels) has been reported in the literature. The
5 outdoor temperature was 288 K, the air density 1.225 𝑘𝑔/𝑚3 , C𝑝 1.0 𝑘𝐽/(𝑘𝑔 ∙ 𝐾),
6 heat source 1 kW, 𝐻1 5.5 m, 𝐻2 5.5 m, 𝑆1+2 37.2933 𝑘𝑔−1 ∙ 𝑚−1 , and gravity
7 acceleration g 9.81 𝑚/𝑠 2 . The comparison in temperature difference between the
8 proposed model and the validated results from [44] is illustrated in Fig. 3(a); the
9 maximum relative error is 13.2% when the strength of the local heat source is 100 W.
10 The maximum relative error for the flow rate is 15.9%, as indicated in Fig. 3(b), when
11 the strength of the local heat source is 100 W. It is clear that the relative error is small
12 for status 2 compared with status 1. For status 2, the airflow enters from the right stack
13 and is assisted by local buoyancy at the left corner, and displacement ventilation is
14 established. For status 1, the outdoor air enters from the left stack and combats with the
15 local buoyancy at the left corner of the room, and mixed ventilation is established.
16 Hence, the room height should be deducted from the stack height when the buoyancy
17 pressure is calculated in zone 2. After this ratification, the maximum relative error of
18 the flow rate is 12.31%, while the maximum relative error of the temperature difference
19 is 10.39%. In general, the modeled results agree well with the validated results from
20 the previous study.
21
19
6
o
5
4
S1 1&S2 2 Model
o
3 o S1 2&S2 1 Model
T C
o S1 2 Rectified Model
2 o o S2 1 Literature
o o S1 1&S2 2 Literature
1
S1 2 Literature
0 o o o o o o
1
0.0 0.2 0.4 0.6 0.8 1.0
E1 Kw
1 (a)
2
0.2
o
o o
0.1 o o
o o S1 airflow Literature
S2 airflow Literature
q kg s
0.0
S1 airflow Model
S2 airflow Model
0.1 S1 airflow Rectified Model
0.2
0.0 0.2 0.4 0.6 0.8 1.0
E1 Kw
3 (b)
4
5 Fig. 3. Modeled results validation: (a) Temperature comparison between the two-zone model and
6 previous CFD results; (b) Mass flow rates comparison between the two-zone model and previous
7 CFD results. (S1 indicates status1/realization 1, #-(1,2) indicates zone number)
8 2.4 Graphical representation
9 To further illustrate the nonlinear dynamical process, we used the fourth-order
10 Runge–Kutta method to numerically solve the differential equation system for different
20
1 scenarios. Different initial zone temperatures were used to produce the different
2 trajectories of the system in the phase portrait. We superimposed the vector field and
3 phase portrait together such that the formation of multiple steady states could be
4 observed more straightforwardly. The existence and stability of the nonlinear ordinary
5 equation system can be demonstrated by the phase portrait and vector field.
6 Furthermore, the results agree well with the analysis from the characteristic equation of
7 the differential equation system. We selected some typical cases to demonstrate this
8 dynamical process: Section 2.4.1 describes scenario 1, where the strength of the heat
9 ratio is fixed, and the height ratio is the control parameter; Section 2.4.3 presents
10 scenario 2, where the height ratio is fixed and the heat ratio is the control parameter;
11 Section 2.4.3 investigates scenario 3, where a single heat source exists at the bottom of
12 the building with two adiabatic tunnels.
13 2.4.1 Same heat source strength, different stack heights (scenario 1)
14 For scenario 1, we used 𝜅 = 1 as an example to graphically study the effect of
15 control parameter 𝛼 on the stability of the buoyancy ventilation system. All the
16 building configurations, including the building geometry, boundary conditions, and
17 thermal properties, are the same as those of the validated case except the height ratio
18 and number of heat sources. Subsequently, the bifurcation diagram, vector field, and
19 phase portrait were determined. As shown in Fig. 4.(a), when 0 < 𝛼 < 1/2 , the
20 ventilation system has only one stable solution for realization 2, in which the flow
1
21 decreases with the increase in 𝛼. When < 𝛼 < 5/9, the system has two solutions;
2
22 the solution for realization 1 is unstable, and that for realization 2 is stable. The flow
23 rates for realization 1 increases with 𝛼, and the flow rates for realization 2 decreases
24 with the increase in 𝛼. When 5/9 < 𝛼 < 9/5, the system has two stable solutions;
25 within this interval, both equilibrium solutions may occur. Under the same conditions,
26 the system can be transferred from one steady state to another steady state with the
21
9
1 appropriate disturbance. When < 𝛼 < 2 , realization 1 remains stable, while
5
2 realization 2 changes from stable to unstable. With the increase in 𝛼, the flow rates for
3 realization 1 increases, while that for realization 2 approaches zero gradually. When
4 𝛼 > 2, the system has only one stable equilibrium solution for realization 1. Therefore,
5 when 𝛼 is 1/2, 5/9, 9/5, or 2, bifurcation will occur when the stability and existence of
6 solutions of the underground buoyancy ventilation system change at these points.
7 Fig. 4(b) shows a bifurcation diagram of the two-zone air temperature as a function
8 of 𝛼. When 0 < 𝛼 < 1/2, only status 2 has a stable solution, and the fluid flows from
9 zone 2 to zone 1. With the increase in 𝛼, the temperatures of both zones increase, which
10 is consistent with the decrease in mass flow rate in Fig. 4(a). Because the mass flow
11 rate decreases and the heat strength of the heat sources remains constant, the indoor
1
12 temperature will increase such that the heat balance can be satisfied. When <𝛼<
2
13 5/9, the ventilation system has two solutions, realization 1 is an unstable solution; the
14 temperature difference between the two zones increases with 𝛼, and the temperature in
15 each zone decreases continuously. However, the state is not stable, while state 2 remains
16 stable. When 5/9 < 𝛼 < 9/5, both realizations are stable. In realization 1, with the
17 increase in 𝛼 , the temperature difference between the two zones does not change
18 significantly, but the height of tunnel 2 increases. Therefore, the thermal pressure of the
19 system increases, and the mass flow from zone 1 to zone 2 increases continuously,
20 which is consistent with the change trend of the mass flow in Fig. 4(a). In realization 2,
21 with the increase in 𝛼 , the temperature difference between the two zones does not
22 change significantly, while the indoor air temperature increases. Because the height of
23 tunnel 2 increases while that of tunnel 1 remains constant, the thermal pressure in tunnel
24 2 increases faster compared with that of tunnel 1. Therefore, the overall thermal
25 pressure in the system decreases, thereby resulting in a decrease in mass flow rate, as
26 shown in Fig. 4(a). When 9/5 < 𝛼 < 2, status 1 remains stable, while status 2 changes
27 from stable to unstable. With the increase in 𝛼, the temperature in the two zones for
28 realization 1 continues to decrease, while the temperature difference between the two
22
1 zones does not change significantly. Therefore, when 𝛼 increases, the height of tunnel
2 2 increases and the overall thermal pressure increases, thereby resulting in an increase
3 in the overall mass flow, as illustrated in Fig. 4(a). For realization 2, as the height of
4 tunnel 2 increases, the thermal pressure in tunnel 2 continues to increase, and the
5 direction of thermal pressure in tunnel 2 is opposite to the direction of flow for
6 realization 2. Therefore, the resistance of the thermal pressure ventilation system will
7 continue to increase. The airflow of the system will decrease continuously until 𝛼 is
8 approximately 2. At this time, the height of the tunnel 2 is twice that of tunnel 1, and
9 the ventilation of the system is approximately zero. At this time, the temperatures of
10 zones 1 and 2 are approximately infinite. It is impossible to achieve a near infinity
11 temperature, which represents an unstable solution. When 𝛼 > 2, no solution exists for
12 realization 2, while realization 1 is stable. For realization 1, the temperature difference
13 between the two zones is not obvious, while the air temperature in both zones decreases.
14 This is due to the increasing height of tunnel 2, which results in the increase in overall
15 thermal pressure in the system. The increase in the overall thermal pressure will
16 increase the overall ventilation, which is consistent with the change trend of airflow in
17 Fig. 4(a).
1.0
0.5
q kg s
0.0
0.5
S1 unstable S1 stable
S2 stable S2 unstable
1.0
0 1 2 3 4
18 (a)
23
60
S1 1 unstable
50 S1 2 unstable
S2 1 stable
40 S2 2 stable
S1 1 stable
T C
30 S1 2 stable
S2 1 unstable
20 S2 2 unstable
10
0
0 1 2 3 4
1 (b)
2 Fig. 4. Bifurcation diagram for scenario1 (𝜅 = 1): (a). Bifurcation diagram of mass flow rates;
3 (b). Bifurcation diagram of temperature difference.
4 The analysis and description process above were used to analyze how 𝛼 affects
5 the equilibrium results. In the vector field and phase portrait, we compared the
6 formation process of the system from the initial condition to the steady state under
7 different ranges of the control parameter 𝛼 . Based on the bifurcation diagram, five
8 typical values of 𝛼 were selected, which were 0.4, 0.54, 1.2, 1.9, and 3. In Fig. 5(a),
9 only a stable fixed point appeared; all different initial conditions converged to the same
10 steady state. As shown in Fig. 5 (b), a straight line ∆𝑇2 = ∆𝑇1 /𝛼 appeared, which
11 divided the vector field into two parts. The part above represents realization 1, where
12 the flow enters from the left tunnel and exits from the right tunnel. For the phase plane
13 below the straight line, airflow enters from the right tunnel and exits from the left tunnel,
14 which corresponds to realization 2. Unless the initial condition is exactly equal to the
15 fixed point in realization 1, the trajectories will converge to the stable fixed point in
16 realization 2. This can be realized by solving Eqs.2-15 and 2-16 or Eqs. 2-17 and 2-18
17 through the fourth-order Runge–Kutta method. For the initial condition above the
18 straight line (∆𝑇2 = ∆𝑇1 /𝛼), the flow direction changes when the trajectories cross the
19 straight line, as indicated in Fig. 5(b). However, for the initial condition below the
20 straight line, the flow will remain in the same direction until it converges to the same
24
1 fixed point. As shown in Fig. 5(c), for 𝛼 = 1.2, two stable fixed points appear in the
2 phase plane, where the initial conditions separated by the straight line will converge to
3 the corresponding fixed point. No flow direction changes will occur provided that no
4 large disturbance occurs. As shown in Fig. 5(d), two fixed points exist in the phase
5 plane. However, the fixed point in realization 2 is not stable. All the initial conditions
6 except the exact fixed point will follow the trajectories and converge to the fixed point
7 in realization 1. In fact, even though a fixed point exists in realization 2, the state cannot
8 exist in the real world because the initial condition cannot exactly the same as the fixed
9 point. Similar to that shown in Fig. 5(b), the initial conditions below the straight line
10 (∆𝑇2 = ∆𝑇1 /𝛼 ) will experience a direction change when the trajectories cross the
11 straight line, as illustrated in Fig. 5(d). The phase portrait indicated in Fig. 5(e) shows
12 that only a stable fixed point exists for the entire phase plane when 𝛼 = 3. All initial
13 conditions will eventually converge to the fixed point.

14
25
1
2 Fig. 5. Phase portrait and vector field for scenario 1: (a). Phase portrait and vector field for 𝜅 =
3 1, 𝛼 = 0.4; (b). Phase portrait and vector field for 𝜅 = 1, 𝛼 = 0.54;(c). Phase portrait and
4 vector field for 𝜅 = 1, 𝛼 = 1.2;(d).Phase portrait and vector field for 𝜅 = 1, 𝛼 = 1.9; (e).
5 Phase portrait and vector field for 𝜅 = 1, 𝛼 = 3
26
1 2.4.2 Same stack height (Scenario 2)
2 For scenario 2, we used 𝛼 = 1 as an example to graphically study the effect of
3 control parameter 𝜅 on the stability of the buoyancy ventilation system. As shown in
4 Fig. 6(a), when𝜅 < 0, only one stable solution exists for realization 2. When 0 < 𝜅 <
5 0.2, the ventilation system has a stable solution for realization 2 and an unstable solution
6 for realization 1. The flow rate for realization 2 is a fixed value, while that for realization
7 1 increases with 𝜅 . When 0.2 < 𝜅 < 5, the system has two stable solutions; both
8 equilibrium solutions may appear depending on the initial conditions. Under the same
9 conditions, the system can be transferred from one stable solution to another with the
10 appropriate disturbance. When 5 < 𝜅, state 1 remains stable, while state 2 changes from
11 stable to unstable. With the increase in 𝜅 , the flow rate for realization 1 increases.
12 Therefore, the flow bifurcation will occur when 𝜅 is 0, 0.2, or 5. The stability of the
13 solution changes at the corresponding points.
14 Fig. 6(b) shows a bifurcation diagram of the two-zone air temperature as a function
15 of 𝜅. When 𝜅 < 0, only realization 2 has a stable solution, where zone 2 contains a
16 heat sink and zone 1 acts as a heat source. Fluid flows from zone 2 to zone 1. With the
17 increase in the heat ratio, the air temperatures in both zones increase continuously, but
18 the temperature difference between the two zones remains constant, which is consistent
19 with the airflow results in Fig 6(a). When the temperature difference is constant,
20 because the height of the shaft does not change, the thermal pressure of the system
21 remains unchanged. Therefore, the mass flow rate remains unchanged. When 0 < 𝜅 <
22 0.2, an unstable solution exists for realization 1, while a stable solution for realization
23 2. For realization 1, the temperature difference between the two zones increases with
24 the increase in heat ratio. When 0.2 < 𝜅 < 5, two stable solutions exist. For realization
25 1, the temperature difference between the two zones continues increasing. However,
26 the air temperature of zone 1 decreases, while the temperature of zone 2 increases. The
27 heat released to zone 1 does not change with the increase in outdoor air flow rate;
27
1 therefore, the temperature decreases. For zone 2, although the mass flow rate increases,
2 the heat released to zone 2 increases with 𝜅 . Therefore, the temperature of zone 2
3 increases. When 5 < 𝜅, realization 1 remains stable, while realization 2 changes from
4 stable to unstable. For realization 1, the temperature difference between the two zones
5 continues to increase, which is due to the increase in mass flow rate in the interval of
6 0.2 < 𝜅 < 5. For realization 2, because the temperatures of both zones will increase to
7 infinite, which is not achievable, the equilibrium solution can be obtained.

q kg s
S1 unstable S1 stable
0.4 S2 stable S2 unstable
0.2
0.0 E2 E1
0.2
1 0 1 2 3 4 5 6
8 (a)
T C
40
S1 1 unstable S1 1 stable
S1 2 unstable S1 2 stable
S2 1 stable S2 1 unstable
30
S2 2 stable S2 2 unstable
20
10
0 E2 E1
1 0 1 2 3 4 5 6
9 (b)
10 Fig. 6. Bifurcation diagram for scenario 2(𝛼 = 1): (a). Bifurcation diagram of mass flow rates;
11 (b). Bifurcation diagram of temperature difference.
12 The analysis and description process above were used to study the effect of 𝜅 on
13 the equilibrium results. In the vector field and phase portrait, the typical value of 𝜅
28
1 was selected based on the bifurcation diagram to compare the dynamical process of the
2 system from the initial condition to the steady state. Fig. 7(a) describes that all the initial
3 values will follow the corresponding trajectory and converge to the same stable fixed
4 point when 𝜅 = −1. Fig. 7(b) reveals that two fixed points exist in the phase plane
5 when 𝜅 = 0.1 . The fixed point in realization 1 is unstable, while all the initial
6 conditions, except the fixed point for realization 1, will converge to the stable fixed
7 point in realization 2. For the initial value above the straight line (∆𝑇2 = ∆𝑇1 /𝛼), they
8 have to cross the straight line where the flow direction is changed. Fig.7(c) shows that
9 two stable fixed points exist at both sides of the straight line (∆𝑇2 = ∆𝑇1 /𝛼). All the
10 initial values will follow their trajectories and converge to the respective fixed points
11 unless the disturbances are imposed. Fig. 7(d) shows that two fixed points exist when
12 𝜅 = 6. The fixed point in realization 2 is unstable, while the fixed point in realization
13 1 is stable. Furthermore, the initial value below the straight line (∆𝑇2 = ∆𝑇1 /𝛼) must
14 experience a change in flow direction before it reaches the stable fixed point in
15 realization 1.
16
29
1
2 Fig. 7. Phase portrait for scenario 2: (a). Phase portrait and vector field for 𝜅 = −1, 𝛼 = 1; (b).
3 Phase portrait and vector field for 𝜅 = 0.1, 𝛼 = 1;(c). Phase portrait and vector field for 𝜅 =
4 1, 𝛼 = 1;(d). Phase portrait and vector field for 𝜅 = 6, 𝛼 = 1.
5 2.4.3 Single heat source at the bottom of the building (scenario 3)
6 For scenario 3, all of the parameters are the same as those of the validated case. As
7 proven in Section 2.2.3, no bifurcation exists in this building configuration. Hence, no
8 bifurcation diagram is available in this scenario. However, Fig. 3 shows the changes in
9 airflow rate and air temperature with the change in the strength of the local heat source.
10 Furthermore, the phase portrait and vector field were determined according to the fourth
11 Runge–Kutta method. As illustrated in Fig. 8, two stable fixed points exist in the phase
12 plane. The vector field is divided by the straight line (∆𝑇2 = ∆𝑇1 /𝛼), which represents
13 the balance in thermal pressure between zones 1 and 2.Above the straight line, all the
14 flow will reach the fixed point in realization 1, while the initial value below the straight
15 line will converge to the fixed point in realization 2. If the initial value is below the
30
1 straight line, but the fixed point in realization 1 is to be reached, then disturbance has
2 to be imposed such that the phase portrait can be changed for a period of time. When
3 the disturbance disappeared, if the air temperature is located at the top of the straight
4 line, then it can reach the fixed point in realization 1. The disturbance can be the wind
5 effect or mechanical fan power. In other words, this phase portrait is drafted according
6 to the fixed building configurations, which include the building geometry, boundary
7 conditions, and all physical properties. The patterns of the phase portrait and vector
8 field will change according to the building configuration. However, if the derived
9 criterion based on the heat ratios and height ratios is adhered, the general trend will
10 follow that of the criterion, e.g., the number of fixed points and their stability.
T2 C
20
10
0 T1 C
10
20
20 10 0 10 20
11
12 Fig. 8. Phase portrait and vector field for scenario 3 (one local heat source at the bottom).
13 3. Case study
14 In this section, the natural ventilation of a hydropower station in Xinjiang is used

15 as an example to analyze the polymorphism of buoyancy ventilation. As shown in Figs.
16 9 (a) and 9 (b), four generators were used in the hydropower station. The main heating
17 equipment in the plant includes generators, busbar cables, main transformers, lightings,
31
1 and control cabinets. As shown in Fig. 9(b), the transportation tunnel is on the left.
2 Generators and other main equipment are in the main workshop. The outlet tunnel of a
3 busbar cable and exhaust shaft is shared. The generator measures 3.6 m in diameter and
4 1.2 m in height; the main factory measures 36 m (L) × 9 m (W) × 8 m (H) ; the
5 transportation tunnel measures 3.1 m (L) × 4.5 m (W) in cross section, 170 m in length,
6 and 86 m in height; the exhaust tunnel measures 140 m in height, and the diameter of
7 the exhaust shaft is approximately 6 m.
10 Fig. 9. Basic information of the hydropower station.
11 Field measurements were conducted in the summer, and the censoring points and
12 setups are shown in Fig. 9(c). Testo 480 was used to capture the air temperature and air
13 velocity. The airflow rate was measured at the junction between the transportation
14 tunnel and the main factory. The mean mass flow rate is 3.22 m3 /s; the outdoor
15 temperature T0 is 23.2 °C; the temperature at the entrance of the factory T1 is
16 14.4 °C; T1 ′ 𝑎𝑛𝑑 T2 are 25.6 °C and 22.1 °C, respectively. The temperature changed
17 significantly in both the transportation tunnel and the exhaust tunnel. However, for the
18 nonlinear dynamical analysis, we still used the well-mixed assumption, because this
19 assumption will not significantly affect the solution multiplicity analysis. Based on
20 simple calculations, the heat loss through the transportation tunnel is 34.9 kW, the
32
1 overall heat gain from the factory is 44.4 kW, and the overall heat loss through the
2 exhaust tunnel is 13.9 kW; H1is 86 m and H2 is 140 m; hence, 𝛼 is 1.628, and S is
3 1.55 𝑘𝑔−1 ∙ 𝑚−1. The outdoor temperature is 23.2 °C, air density 1.22 kg ∙ 𝑚−3 , and
4 C𝑝 1.01 kJ/(𝑘𝑔 ∙ 𝐾). The thermal mass of the air inside the tunnel was considered
5 without considering other building elements. The thermal masses in the transportation
6 tunnel, factory, and exhaust tunnel are 2372, 2592, and 3956 kg, respectively.
7 For realization 1, the airflow entered from the transportation tunnel, and the thermal
8 mass and overall heat gain from the factory were transferred to zone 2; therefore, 𝑀1
9 is 2372 kg, 𝑀2 is 6548 kg, 𝐸1 is -34.9 kW, and 𝐸2 is 30.1 kW. Therefore, we have
𝑑∆𝑇1 𝜌𝑎 𝑔𝐻1
10 𝑀1 𝐶𝑝 = −√𝑇 √−∆𝑇1 + ∆𝑇2 𝛼 𝐶𝑝 ∆𝑇1 + 𝐸1 （3-1）
𝑑𝑡 𝑎 (𝑆1 +𝑆2)
𝑑𝑇2 𝜌𝑎 𝑔𝐻1
11 𝑀2 𝐶𝑝 = −√𝑇 √−∆𝑇1 + ∆𝑇2 𝛼 𝐶𝑝 (∆𝑇2 − ∆𝑇1 ) + 𝐸2 （3-2）
𝑑𝑡 𝑎 (𝑆1 +𝑆2)
12 The root of the characteristic equation of the differential equation is 𝜆1 =
13 −0.0034，𝜆2 = −0.0004321. Both roots are real distinct negative value; hence, a
14 stable fixed point exists for realization 1.
15 For realization 2, the airflow entered from the exhaust tunnel, and the thermal mass
16 and overall heat gain from the factory were transferred to zone 1; therefore, 𝑀3 is 4964
17 kg, 𝑀4 is 3956 kg, 𝐸3 is 9.5 kW, and 𝐸3 is -13.9 kW. We have

18 𝑀3 𝐶𝑝 = −√𝑇 √∆𝑇1 − 𝛼∆𝑇2(∆𝑇1 − ∆𝑇2 ) + 𝐸3 （3-3）
𝑑𝑡 𝑎 (𝑆1 +𝑆2)
19 𝑀4 𝐶𝑝 = −√𝑇 √∆𝑇1 − 𝛼∆𝑇2 ∆𝑇2 + 𝐸4 （3-4）
𝑑𝑡 𝑎 (𝑆1 +𝑆2)
20 The root of the characteristic equation of the differential equation is
21 𝜆1 = −0.001765，𝜆2 = −0.0005075 . Both roots are real distinct negative values;
22 hence, a stable fixed point exists for realization 2 as well.
23 Based on Eqs. 3-1, 3-2, 3-3, and 3-2, we numerically solved these differential
24 equation systems. The vector field and phase portrait are shown in Fig.10. The stability
25 and existence of the system obtained from the numerical computation is the same as
26 those from the nonlinear dynamical analysis. Two stable fixed points exist, and the
33
1 vector field is divided by the straight line (∆𝑇2 = ∆𝑇1/𝛼).
2
T2 C
20
10
3 0 T1 C
10
20
20 10 0 10 20
4 Fig. 10. Phase portrait and vector field for the dynamic system of buoyancy ventilation in the
5 hydropower station
6 In fact, we can use the derived criterion in scenario 1 to evaluate the stability and
7 existence of fixed points:

1+𝜅
8 When 𝜅 < −1: if 0 < 𝛼 < , no fixed point exists for realizations 1 and 2; if α
𝜅
1+𝜅
9 > , no fixed point exists for realization 1but a stable fixed point exists for realization
𝜅
10 2.
11 In scenario 1, we assume that 𝐸1 is constantly positive, while in this case study,
12 we have a negative value for 𝐸1 . Hence, both realizations 1 and 2 in the case study
13 corresponds to realization 2 in the criterion.
14 1 𝐸
For realization 1 in the case study, 𝐸1 is negative, 𝜅 = 𝐸 = −1.159468, and 𝛼 =
2
𝐻1
15 = 0.614286; 1+𝜅 = 0.1375358, that is α > 1+𝜅. Therefore, a stable fixed point
𝐻2 𝜅 𝜅
16 exists according to the criterion. At the same time, when 𝐻2 remains unchanged, if the
17 1+𝜅
height of the left traffic tunnel 𝐻1 is lower than 19.26 m (α < 𝜅
), this fixed point
18 will not exist.
34
1 𝐸
4
For realization 2 in the case study, 𝐸4 is negative, 𝜅 = 𝐸 = −1.46316， and 𝛼 =
3
2 𝐻2 1+𝜅 1+𝜅
= 1.627907; = 0.3165468，that is α > . Therefore, a stable fixed point
𝐻1 𝜅 𝜅
3 exists according to the criterion. At the same time, when 𝐻1 remains unchanged, if the
4 1+𝜅
height of the left traffic tunnel 𝐻2 is lower than 27.2 m (α < 𝜅
), this fixed point will
5 not exist.
6 From the analysis and verification above, it can be concluded that the criterion is
7 applicable to the evaluation of buoyancy ventilation polymorphism of two-zone
8 underground buildings. If the specific heat source ratio, height ratio, and flow direction
9 are known, then the existence and stability of the solution can be evaluated according
10 to the criterion. In the design stage, according to the criterion, the natural ventilation of
11 buildings can be optimized by selecting the appropriate height ratio of the shaft and the
12 distribution of heat sources to avoid unfavorable solutions and induce favorable
13 solutions.
14 4. Conclusions
15 In summary, nonlinear dynamical analysis was performed to study the buoyancy
16 ventilation of a typical two-zone underground building with different building
17 configurations. One dimensional model that described the buoyancy ventilation of the
18 two-zone underground buildings was proposed and validated by results from previous
19 studies. The model comprised two groups of mathematical nonlinear ordinary
20 differential equation systems. Three different scenarios were studied based on the
21 differential equation systems, and the corresponding criterion for the stability and
22 existence of fixed points for the underground buoyancy ventilation was derived
23 mathematically. The criterion was based on the heat ratio (𝜅) and height ratio (α). For
24 scenario 1 (𝜅 was fixed, control α ) and scenario 2 (α was fixed, control 𝜅 ), the
25 criterion is summarized in Tables 1 and 2. For scenario 3 (one heat source at the bottom
35
1 of the building, control 𝛼), two stable fixed points appeared.
2 Subsequently, the bifurcation diagram, phase portrait, and vector field for the three
3 scenarios were produced based on numerical computation by applying the fourth
4 Runge–Kutta method, which produced the similar result as the characteristic equations.
5 Finally, a case study was conducted based on a real project with field
6 measurements, which demonstrated the use of the nonlinear dynamical analysis method
7 to evaluate the stability and existence of fixed points for the buoyancy ventilation in a
8 hydropower station. This case study served as a validation case for the derived criterion,
9 which demonstrated its capability in predicting the existence and stability of fixed
10 points in the buoyancy ventilation of underground buildings.
11 Acknowledgement
12 The authors acknowledge the support from National Natural Science Foundation of China (NSFC)
13 (51678088, 51178482, 51578087)and National Key R&D Program of China（2017YFB0903700）.
14 Yanan Liu also especially thanks for the scholarship provided by China Scholarship Council (CSC
15 student ID: 201706050003).
16
17 Appendix A1 Dynamical analysis of Eqs.2-15 and 2-16
18 Denote the steady state solution(fixed point) of realization 1(Eqs.2-15 and 2-16)
19 ̅̅̅̅̅1，∆𝑇
is (∆𝑇 ̅̅̅̅̅2 )，we can obtain：
20 ̅̅̅̅̅2 − ̅̅̅̅̅
𝐸1 = 𝑛√𝛼∆𝑇 ∆𝑇1 ̅̅̅̅̅
∆𝑇1，（A1）
21 ̅̅̅̅̅2 − ̅̅̅̅̅
𝐸2 = −𝑛√𝛼∆𝑇 ̅̅̅̅̅1 − ̅̅̅̅̅
∆𝑇1 (∆𝑇 ∆𝑇2 ) （A2）
22 Eq.(A1)divided by Eq.(A2)：
𝐸1 ̅̅̅̅̅
∆𝑇1
23 = ̅̅̅̅̅ ̅̅̅̅̅ （A3）
𝐸2 ∆𝑇 −∆𝑇
2 1
𝐸2
24 Substituting = 𝜅 into Eq.(A3) and simplifying the equation, we obtain:
𝐸1
36
𝐸
1 ̅̅̅̅̅
∆𝑇2 = (1 + 𝐸2) ̅̅̅̅̅
∆𝑇1 = (1 + 𝜅) ̅̅̅̅̅
∆𝑇1 （A4）
1
𝑑∆𝑇1 𝑑∆𝑇2
2 Substituting Eq.(A4) into Eqs.2-15 and 2-16，Let and be zero，fixed point
𝑑𝑡 𝑑𝑡
3 can be obtained：
𝐸1 𝐸1
4 ̅̅̅̅̅1 = 2⁄3
∆𝑇 1⁄ 3 = 𝑛2⁄3(−𝐸 1⁄3 （A5）
𝑛 (−𝐸 1+𝐸1 𝛼+𝐸2 𝛼) 1 +𝐸1𝛼+𝐸1 𝛼𝜅)
(𝐸 +𝐸 ) 𝐸1 𝐸1 +𝐸1 𝜅
5 ̅̅̅̅̅
∆𝑇2 = 1𝐸 2 (𝑛2⁄3(−𝐸 ) = 𝑛2⁄3(−𝐸 （A6）
1 1 +𝐸1 𝛼+𝐸2 𝛼)1⁄3 1 +𝐸1𝛼+𝐸1 𝛼𝜅)
1⁄ 3
6 To analyze the stability of the nonlinear ordinary differential equations above, we
7 need to linearize the equations at the fixed point, and analyze their existence and
8 stability according to the characteristic equation. Eqs.2-15 and 2-16 can be written as
9 ̅̅̅̅̅1, ̅̅̅̅̅
Taylor expansions form at the fixed point (∆𝑇 ∆𝑇2 )：
10
𝜕𝑓 𝜕𝑓
11 ̅̅̅̅̅1 , ̅̅̅̅̅
𝑓1 (∆𝑇1 , ∆𝑇2 ) = 𝑓1 (∆𝑇 ∆𝑇2) + 𝜕∆𝑇1 |∆𝑇1=∆𝑇
̅̅̅̅̅ (∆𝑇1 − 𝑇̅1 ) + 𝜕∆𝑇1 |∆𝑇1=∆𝑇
̅̅̅̅̅ (∆𝑇2 − ̅̅̅̅̅
∆𝑇2 )
1 1 1 2
̅̅̅̅̅
∆𝑇2 =∆𝑇 2
̅̅̅̅̅
∆𝑇2 =∆𝑇 2
12 （A7）
𝜕𝑓 𝜕𝑓
13 ̅̅̅̅̅1 , ̅̅̅̅̅
𝑓2 (∆𝑇1 , ∆𝑇2 ) = 𝑓2 (∆𝑇 ∆𝑇2 ) + 𝜕∆𝑇2 |∆𝑇1 =∆𝑇
̅̅̅̅̅ (∆𝑇1 − 𝑇̅1 ) + 𝜕∆𝑇2 |∆𝑇1=∆𝑇
̅̅̅̅̅ (∆𝑇2 − ̅̅̅̅̅
∆𝑇2 )
1 1 1
2
̅̅̅̅̅
∆𝑇2 =∆𝑇 2 ̅̅̅̅̅
∆𝑇2 =∆𝑇 2
14 （A8）
15 For the steady state， and ̅̅̅̅̅1，∆𝑇
is equal to zero. Thus 𝑓1 (∆𝑇 ̅̅̅̅̅2 ) = 0，and
𝑑𝑡 𝑑𝑡
16 ̅̅̅̅̅1 ，∆𝑇
𝑓2 (∆𝑇 ̅̅̅̅̅2) = 0.
𝜕𝑓1 𝜕𝑓1 𝜕𝑓2
17 Assuming𝛼11 = | ̅̅̅̅̅ ，𝛼12 = | ̅̅̅̅̅ ，𝛼21 = | ̅̅̅̅̅ ，and 𝛼22 =
𝜕∆𝑇1 ∆𝑇1 =∆𝑇1 𝜕∆𝑇2 ∆𝑇1 =∆𝑇1 𝜕∆𝑇1 ∆𝑇1 =∆𝑇1
̅̅̅̅̅
∆𝑇2 =∆𝑇 2 ̅̅̅̅̅
∆𝑇2 =∆𝑇 2 ̅̅̅̅̅
∆𝑇2 =∆𝑇 2
𝜕𝑓2
18 | ̅̅̅̅̅ , Eqs. (A7) and (A8) can be written as follows:
𝜕∆𝑇2 ∆𝑇1 =∆𝑇1
̅̅̅̅̅
∆𝑇2 =∆𝑇 2
19 𝑓1 (𝑇1, 𝑇2 ) = 𝛼11 (𝑇1 − 𝑇̅1 ) + 𝛼12 (𝑇2 − 𝑇̅2 ) （A9）
20 𝑓2 (𝑇1, 𝑇2 ) = 𝛼21 (𝑇1 − 𝑇̅1 ) + 𝛼22 (𝑇2 − 𝑇̅2 ) （A10）
21 Written as matrix form：

𝑑𝑇1
𝑑𝑡 𝑇
22 [𝑑𝑇 ] = 𝐴1 [ 1 ] + 𝐵 (A11)
2 𝑇2
𝑑𝑡
𝛼11 𝛼12 −𝛼11𝑇̅1 − 𝛼12 𝑇̅2
23 Where 𝐴1 = [𝛼 𝛼22 ], and 𝐵 = [ ]. To evaluate the existence and
21 −𝛼21𝑇̅1 − 𝛼22 𝑇̅2
24 stability of the solution, we have to obtain the eigenvalue of matrix 𝐴1 .
37
1 The characteristic equation of realization 1 can be easily obtained as follows:
(𝐸 (−1+𝛼+𝛼𝜅))2⁄3
𝑛√ 1 2⁄ 3 (−5+𝛼(4+5𝜅))𝜆
𝑛 3
2 𝜆2 + + 2 𝑛 4⁄3(𝐸1 (−1 + 𝛼 + 𝛼𝜅))2⁄3 = 0 (A12)
2(−1+𝛼+𝛼𝜅)
3
𝐸
4 Assuming 𝜅 = 𝐸2 > 0, to ensure that a real solution exists for Eqs. (A5) and (A6),
1
5 we have −𝐸1 + 𝐸1 𝛼 + 𝐸1 𝛼𝜅 > 0. Because the heat source in zone 2 is positive as well,
1
6 we have 𝜅 > 0. Therefore, 𝛼 > 1+𝜅.
(𝐸 (−1+𝛼+𝛼𝜅)) 2⁄3
𝑛√ 1 2⁄ 3 (−5+𝛼(4+5𝜅))
𝑛 3
7 Assuming𝛽 = , 𝛾 = 2 𝑛 4⁄3(𝐸1 (−1 + 𝛼 + 𝛼𝜅))2⁄3 .The
2(−1+𝛼+𝛼𝜅)
8 root of the characteristic equation (Eq.A12) is the eigenvalue of matrix 𝐴1 . The
9 eigenvalue is given as 𝜆1,2 = −𝛽 ± √𝛽2 − 4𝛾. Since𝐸1 (−1 + 𝛼 + 𝛼𝜅) > 0, and 𝑛 >
3
10 0, thus 𝛾 = 2 𝑛 4⁄3 (𝐸1 (−1 + 𝛼 + 𝛼𝜅))2⁄3 > 0. Therefore, 𝜆1,2 are determined by the
11 real part of 𝛽. When 𝛽 > 0，the real part of the eigenvalue is negative; otherwise，
12 the real part of the eigenvalue is positive. To ensure the solution of this differential
13 equation is stable, the root of the characteristic equation should be less than zero. That
(𝐸 (−1+𝛼+𝛼𝜅)) 2⁄3
𝑛√ 1 2⁄ 3 (−5+𝛼(4+5𝜅))
𝑛
14 is −𝛽 = − < 0. Therefore, −5 + 𝛼(4 + 5𝜅) > 0，that
2(−1+𝛼+𝛼𝜅)
5
15 is 𝛼 > 4+5𝜅.
1 5
16 In summary，when < 𝛼 < 4+5𝜅, the real parts of two eigenvalues are positive，
1+𝜅
5
17 the fixed point is unstable；when 𝛼 > ，the real parts of two eigenvalues are
4+5𝜅
1
18 negative，the fixed point is stable；when 𝛼 < 1+𝜅， no fixed point exists.
19 Appendix A2 Dynamical analysis of Eqs.2-17 and 2-18
20 Denote the steady state solution(fixed point) of realization 2(Eqs.2-17 and 2-18) is
38
1 ̅̅̅̅̅1 ，∆𝑇
(∆𝑇 ̅̅̅̅̅2 )，we can obtain：
2 ̅̅̅̅̅1 − ̅̅̅̅̅̅̅
𝐸1 = 𝑛√∆𝑇 𝛼∆𝑇2 (∆𝑇̅̅̅̅̅1 − ̅̅̅̅̅
∆𝑇2 )，（A13）
3 ̅̅̅̅̅1 − ̅̅̅̅̅̅̅
𝐸2 = 𝑛√∆𝑇 𝛼∆𝑇2 ̅̅̅̅̅
∆𝑇2 （A14）
4 Eq.(A13) divided by Eq.(A14)：

𝐸1 ̅̅̅̅̅
∆𝑇1 −∆𝑇 ̅̅̅̅̅
2
5 = ̅̅̅̅̅
（A15）
𝐸2 ∆𝑇2
𝐸 1
6 Substituting 𝐸21 = 𝜅 into Eq.(A15)，we can obtain:
7 ̅̅̅̅̅1 = (1 + 𝐸1) ̅̅̅̅̅

∆𝑇
1
∆𝑇2 = (1 + 𝜅) ̅̅̅̅̅
∆𝑇2 （A16）
𝐸 2
8 Substituting Eq.(A16) into Eqs.2-17 and 2-18，Let and be zero，fixed point
𝑑𝑡 𝑑𝑡
9 can be obtained：
1
𝐸2 𝐸1 (1+𝜅)𝜅
10 ̅̅̅̅̅1 = (1 + 𝐸1 ⁄𝐸2 ) 2⁄3
∆𝑇 = 𝑛2⁄3(𝐸 （A17）
𝑛 (𝐸 1 +𝐸2 −𝐸2 𝛼)1⁄3 1 +𝐸1 𝜅−𝐸1 𝛼𝜅)
1⁄ 3
𝐸2 𝐸1 𝜅
11 ̅̅̅̅̅
∆𝑇2 = 𝑛2⁄3(𝐸 1⁄ 3 = 𝑛2⁄3(𝐸 1⁄ 3 （A18）
1 +𝐸2−𝐸2 𝛼 ) 1 +𝐸1 𝜅−𝐸1 𝛼𝜅)
12 The characteristic equation of realization 2 can be easily obtained as follows:
(𝐸 (1+𝜅−𝛼𝜅))2⁄3
𝑛(5+(4−5𝛼)𝜅)√ 1 2⁄ 3
𝜆
𝑛 3
13 𝜆2 + + 2 𝑛 4⁄3 (𝐸1 (1 + 𝜅 − 𝛼𝜅))2⁄3 = 0 （A19）
2(1+𝜅−𝛼𝜅)
14
𝐸
15 Given 𝜅 = 𝐸2 > 0, to ensure that a real solution exists for Eqs. (A17) and (A18), we
1
16 have 𝐸1 + 𝐸1 𝜅 − 𝐸1 𝛼𝜅 > 0, and we consider heat source in zone 2 is positive as well,

1+𝜅
17 we have 𝜅 > 0. Therefore, 𝛼 < 𝜅
.
(𝐸 (1+𝜅−𝛼𝜅)) 2⁄3
𝑛(5+(4−5𝛼)𝜅)√ 1 2⁄3 𝑛 3
18 Assuming𝛽 = , and 𝛾 = 2 𝑛 4⁄3 (𝐸1 (1 + 𝜅 − 𝛼𝜅))2⁄3 .
2(1+𝜅−𝛼𝜅)
19 The root of Eq.(A19) is the eigenvalue of matrix𝐴1 . The eigenvalue is given as 𝜆1,2 =
3
20 −𝛽 ± √𝛽2 − 4𝛾 . Since 𝐸1 (1 + 𝜅 − 𝛼𝜅) > 0, and 𝑛 > 0 , thus 𝛾 = 2 𝑛 4⁄3 (𝐸1 (1 +
21 𝜅 − 𝛼𝜅))2⁄3 > 0. Hence, 𝜆1,2 are determined by the real part of 𝛽. When 𝛽 > 0，
22 the real part of the eigenvalue is negative; otherwise，the real part of the eigenvalue is
39
1 positive. To ensure the solution of this differential equation is stable, the root of the
2 characteristic equation should be less than zero. That is −𝛽 =

2⁄3
(𝐸 (1+𝜅−𝛼𝜅))
𝑛(5+(4−5𝛼)𝜅)√ 1 2⁄ 3 𝜆
𝑛
3 − < 0. Therefore, 5 + (4 − 5𝛼)𝜅 > 0 ， that is 𝛼 <
2(1+𝜅−𝛼𝜅)
5+4𝜅
4 .
5𝜅
5+4𝜅 1+𝜅
5 In summary，when <𝛼< ，the real parts of two eigenvalues are positive，
5𝜅 𝜅
5+4𝜅
6 the fixed point is unstable；when 𝛼 < ，the real parts of two eigenvalues are
5𝜅
1+𝜅
7 negative，the fixed point is stable；when 𝜅
< 𝛼，no fixed point exists.
8 Appendix B1 Dynamical analysis of Eqs.2-31 and 2-32
9 Denote the steady state solution(fixed point) of realization 1(Eqs.2-31 and 2-32) is
10 ̅̅̅̅̅1 ，∆𝑇
(∆𝑇 ̅̅̅̅̅2 )，we can obtain：
11 ̅̅̅̅̅1 = 0
∆𝑇 （B1）
2⁄ 3
𝐸
12 ̅̅̅̅̅
∆𝑇2 = 𝑛2⁄13𝛼 1⁄3 （B2）
13 After linearization of Eqs. 2-31 and 2-32 at the fixed point，the characteristic equation
14 of realization 1 can be easily obtained as follows：

̅̅̅̅̅
5𝑛√𝛼∆𝑇 2𝜆 3
15 𝜆2 + ̅̅̅̅̅2 = 0
+ 2 𝑛 2 𝛼∆𝑇 （B3）
2
̅̅̅̅̅
5𝑛√𝛼∆𝑇 3
16 Assuming𝛽 = 2 ̅̅̅̅̅2 , by solving Eq.(B3), two eigenvalue of
, and 𝛾 = 2 𝑛 2 𝛼∆𝑇
2
17 ̅̅̅̅̅2 > 0 𝑎𝑛𝑑 𝑛 > 0, hence, 𝛾 =

matrix𝐴1 are given as𝜆1,2 = −𝛽 ± √𝛽2 − 4𝛾.Since 𝛼∆𝑇
3
18 ̅̅̅̅̅2 > 0. Therefore, 𝜆1,2 are determined by the real part of 𝛽 . Since −𝛽 =
𝑛 2 𝛼∆𝑇
2
̅̅̅̅̅
5𝑛√𝛼∆𝑇 2
19 − < 0 always holds，two eigenvalues are always negative，and a stable fixed
2
20 point always exists for realization 1.
40
1 Appendix B2 Dynamical analysis of Eqs.2-33 and 2-34
2 ̅̅̅̅̅1 ，∆𝑇
Assuming the fixed point of realization 2(Eqs.2-33 and 2-34) is (∆𝑇 ̅̅̅̅̅2), we
3 obtain:
2⁄ 3
4 ̅̅̅̅̅1 = 𝐸12⁄3
∆𝑇 （B4）
𝑛
5 ̅̅̅̅̅
∆𝑇2 = 0 （B5）
6 After linearization of Eqs.2-33 and 2-34 at fixed point，the characteristic equation of
7 realization 2 can be easily obtained as follows：

̅̅̅̅̅
5𝑛√∆𝑇 1𝜆 3
8 𝜆2 + + 2 𝑛 2 ̅̅̅̅̅
∆𝑇1 = 0 （B6）
2
̅̅̅̅̅
5𝑛√∆𝑇 3
9 Assuming 𝛽 = 1
, and 𝛾 = 2 𝑛 2 ̅̅̅̅̅
∆𝑇1 , by solving Eq.(B6), two eigenvalue of
2
2⁄ 3
𝐸
10 matrix 𝐴1 are given as 𝜆1,2 = −𝛽 ± √𝛽2 − 4𝛾 . Since ̅̅̅̅̅
∆𝑇1 = 12⁄3 > 0 and 𝑛 > 0 ,
𝑛
3
11 thus 𝛾 = 2 𝑛 2 ̅̅̅̅̅
∆𝑇1 > 0. Therefore, 𝜆1,2 are determined by the real part of 𝛽 . Since
̅̅̅̅̅
5𝑛√∆𝑇 1
12 −𝛽 = − < 0 always holds，two eigenvalues are always negative，and a stable
2
13 fixed point always exists for realization 2.

14
15 Figure captions
16 Fig. 1. Schematics of a typical two-zone underground structure.
17 Fig. 2. Schematics of a typical two-zone underground structure with one local heat source.
18 Fig. 3. Modeled results validation: (a) Temperature comparison between the two zone model and
19 previous CFD results; (b) Mass flow rates comparison between the two zone model and previous
20 CFD results.(S1 indicates status1/realization 1, #-(1,2) indicates zone number )
21 Fig. 4. Bifurcation diagram for scenario1 (𝜅 = 1): (a). Bifurcation diagram of mass flow rates; (b).
22 Bifurcation diagram of temperature difference.
23 Fig. 5. Phase portrait and vector field for scenario 1: (a).Phase portrait and vector field for 𝜅 =
24 1, 𝛼 = 0.4; (b).Phase portrait and vector field for 𝜅 = 1, 𝛼 = 0.54;(c).Phase portrait and vector field
41
1 for 𝜅 = 1, 𝛼 = 1.2 ;(d).Phase portrait and vector field for 𝜅 = 1, 𝛼 = 1.9 ; (e).Phase portrait and
2 vector field for 𝜅 = 1, 𝛼 = 3
3 Fig. 6. Bifurcation diagram for scenario 2(𝛼 = 1): (a). Bifurcation diagram of mass flow rates; (b).
4 Bifurcation diagram of temperature difference.
5 Fig. 7. Phase portrait for scenario 2: (a).Phase portrait and vector field for 𝜅 = −1, 𝛼 = 1; (b).Phase
6 portrait and vector field for 𝜅 = 0.1, 𝛼 = 1 ;(c).Phase portrait and vector field for 𝜅 = 1, 𝛼 =
7 1;(d).Phase portrait and vector field for 𝜅 = 6, 𝛼 = 1.
8 Fig. 8. Phase portrait and vector field for scenario 3(one local heat source at the bottom).
9 Fig. 9. Basic information of the hydro power station.
10 Fig. 10.Phase portrait and vector field for the dynamic system of buoyancy ventilation in the
11 hydropower station.
12
13 Table captions
14 Table 1. Criterion for scenario 1
15 Table 2. Criterion for scenario 2
16 References
17 1. Mukhtar, A., K.C. Ng, and M.Z. Yusoff, Design optimization for ventilation shafts of
18 naturally-ventilated underground shelters for improvement of ventilation rate and thermal
19 comfort. Renewable Energy, 2018. 115: p. 183-198. DOI:
20 https://doi.org/10.1016/j.renene.2017.08.051.
21 2. Porras-Amores, C., et al., Natural ventilation analysis in an underground construction:
22 CFD simulation and experimental validation. Tunnelling and Underground Space
23 Technology, 2019. 90: p. 162-173. DOI: https://doi.org/10.1016/j.tust.2019.04.023.
24 3. Liu, Z., et al., On-site assessments on variations of PM2.5, PM10, CO2 and TVOC
25 concentrations in naturally ventilated underground parking garages with traffic volume.
26 Environmental Pollution, 2019. 247: p. 626-637. DOI:
27 https://doi.org/10.1016/j.envpol.2019.01.095.
28 4. Zhao, Y., et al., Seasonal patterns of PM10, PM2.5, and PM1.0 concentrations in a
29 naturally ventilated residential underground garage. Building and Environment, 2017.
30 124: p. 294-314. DOI: https://doi.org/10.1016/j.buildenv.2017.08.014.
42
1 5. Cheng, L.H., T.H. Ueng, and C.W. Liu, Simulation of ventilation and fire in the
2 underground facilities. Fire Safety Journal, 2001. 36(6): p. 597-619. DOI:
3 https://doi.org/10.1016/S0379-7112(01)00013-3.
4 6. Alkaff, S.A., S.C. Sim, and M.N. Ervina Efzan, A review of underground building
5 towards thermal energy efficiency and sustainable development. Renewable and
6 Sustainable Energy Reviews, 2016. 60: p. 692-713. DOI:
7 https://doi.org/10.1016/j.rser.2015.12.085.
8 7. Mukhtar, A., M.Z. Yusoff, and K.C. Ng, The potential influence of building optimization
9 and passive design strategies on natural ventilation systems in underground buildings:
10 The state of the art. Tunnelling and Underground Space Technology, 2019. 92: p. 103065.
11 DOI: https://doi.org/10.1016/j.tust.2019.103065.
12 8. Liu, Y., et al., A network model for natural ventilation simulation in deep buried
13 underground structures. Building and Environment, 2019. 153: p. 288-301. DOI:
14 https://doi.org/10.1016/j.buildenv.2019.01.045
15 9. Li, A., X. Gao, and T. Ren, Study on thermal pressure in a sloping underground tunnel
16 under natural ventilation. Energy and Buildings, 2017. 147: p. 200-209. DOI:
17 https://doi.org/10.1016/j.enbuild.2017.04.060.
18 10. Li, Y.Z. and H. Ingason, Overview of research on fire safety in underground road and
19 railway tunnels. Tunnelling and Underground Space Technology, 2018. 81: p. 568-589.
20 DOI: https://doi.org/10.1016/j.tust.2018.08.013.
21 11. Li, M., S.M. Aminossadati, and C. Wu, Numerical simulation of air ventilation in super-
22 large underground developments. Tunnelling and Underground Space Technology, 2016.
23 52: p. 38-43. DOI: https://doi.org/10.1016/j.tust.2015.11.009
24 12. Batchelor, G., Heat convection and buoyancy effects in fluids. Quarterly journal of the
25 royal meteorological society, 1954. 80(345): p. 339-358. DOI:
26 https://doi.org/10.1002/qj.49708034504
27 13. Khanal, R. and C. Lei, Solar chimney—A passive strategy for natural ventilation. Energy
28 and Buildings, 2011. 43(8): p. 1811-1819. DOI:
29 https://doi.org/10.1016/j.enbuild.2011.03.035
30 14. Zhai, X.Q., Z.P. Song, and R.Z. Wang, A review for the applications of solar chimneys in
31 buildings. Renewable and Sustainable Energy Reviews, 2011. 15(8): p. 3757-3767. DOI:
32 https://doi.org/10.1016/j.rser.2011.07.013
33 15. Zhou, J. and Y. Chen, A review on applying ventilated double-skin facade to buildings in
34 hot-summer and cold-winter zone in China. Renewable and Sustainable Energy Reviews,
35 2010. 14(4): p. 1321-1328. DOI: https://doi.org/10.1016/j.rser.2009.11.017
36 16. Shameri, M.A., et al., Perspectives of double skin façade systems in buildings and energy
37 saving. Renewable and Sustainable Energy Reviews, 2011. 15(3): p. 1468-1475. DOI:
39 17. Saroglou, T., et al., A study of different envelope scenarios towards low carbon high-rise
40 buildings in the Mediterranean climate - can DSF be part of the solution? Renewable and
41 Sustainable Energy Reviews, 2019. 113: p. 109237. DOI:
43
1 18. Ghaffarianhoseini, A., et al., Exploring the advantages and challenges of double-skin
2 façades (DSFs). Renewable and Sustainable Energy Reviews, 2016. 60: p. 1052-1065.
3 DOI: https://doi.org/10.1016/j.rser.2016.01.130
4 19. De Gracia, A., et al., Numerical modelling of ventilated facades: A review. Renewable and
5 Sustainable Energy Reviews, 2013. 22: p. 539-549. DOI:
7 20. Barbosa, S. and K. Ip, Perspectives of double skin façades for naturally ventilated
8 buildings: A review. Renewable and Sustainable Energy Reviews, 2014. 40: p. 1019-1029.
9 DOI: https://doi.org/10.1016/j.rser.2014.07.192
10 21. Stabat, P., M. Caciolo, and D. Marchio, Progress on single-sided ventilation techniques
11 for buildings. Advances in Building Energy Research, 2012. 6(2): p. 212-241. DOI:
12 https://dx.doi.org/10.1080/17512549.2012.740903
13 22. Lo, L.J. and A. Novoselac. Effect of indoor buoyancy flow on wind-driven cross
14 ventilation. in Building Simulation. 2013. Springer. DOI: https://doi.org/10.1007/s12273-
15 012-0094-3
16 23. Stavridou, A.D. and P.E. Prinos, Natural ventilation of buildings due to buoyancy assisted
17 by wind: Investigating cross ventilation with computational and laboratory simulation.
18 Building and Environment, 2013. 66: p. 104-119. DOI:
20 24. Hunt, G. and P. Linden, The fluid mechanics of natural ventilation—displacement
21 ventilation by buoyancy-driven flows assisted by wind. Building and Environment, 1999.
22 34(6): p. 707-720. DOI: https://doi.org/10.1016/s0360-1323(98)00053-5
23 25. Nitta, K. Variety modes and chaos in smoke ventilation by ceiling chamber system. in
24 Proceedings of the Sixth International IBPSA Conference, Kyoto, Japan. 1999.
25 26. Heiselberg, P., et al., Experimental and CFD evidence of multiple solutions in a naturally
26 ventilated building. Indoor Air, 2004. 14(1): p. 43-54. DOI:
27 https://doi.org/10.1046/j.1600-0668.2003.00209.x
28 27. Li, Y. and A. Delsante, Natural ventilation induced by combined wind and thermal forces.
29 Building and Environment, 2001. 36(1): p. 59-71. DOI: https://doi.org/10.1016/s0360-
30 1323(99)00070-0
31 28. Lishman, B. and A.W. Woods, On transitions in natural ventilation flow driven by
32 changes in the wind. Building and Environment, 2009. 44(4): p. 666-673. DOI:
34 29. Li, Y., et al., Some examples of solution multiplicity in natural ventilation. Building and
35 Environment, 2001. 36(7): p. 851-858. DOI: https://doi.org/10.1016/s0360-
36 1323(01)00011-7
37 30. Yuan, J. and L. Glicksman, Multiple steady states in a combined buoyancy and wind
38 driven natural ventilation system: necessary conditions and initial values. Proceedings of
39 indoor air, 2005: p. 1207-1212.
40 31. Yuan, J. and L.R. Glicksman, Transitions between the multiple steady states in a natural
41 ventilation system with combined buoyancy and wind driven flows. Building and
44
1 Environment, 2007. 42(10): p. 3500-3516. DOI:
3 32. Yuan, J. and L.R. Glicksman, Multiple steady states in combined buoyancy and wind
4 driven natural ventilation: The conditions for multiple solutions and the critical point for
5 initial conditions. Building and Environment, 2008. 43(1): p. 62-69. DOI:
7 33. Gladstone, C. and A.W. Woods, On buoyancy-driven natural ventilation of a room with a
8 heated floor. Journal of fluid mechanics, 2001. 441: p. 293-314. DOI:
9 https://doi.org/10.1017/s0022112001004876
10 34. Pulat, E. and H.A. Ersan, Numerical simulation of turbulent airflow in a ventilated room:
11 Inlet turbulence parameters and solution multiplicity. Energy and Buildings, 2015. 93: p.
12 227-235. DOI: https://doi.org/10.1016/j.enbuild.2015.01.067
13 35. Chenvidyakarn, T. and A. Woods, Multiple steady states in stack ventilation. Building and
14 Environment, 2005. 40(3): p. 399-410. DOI:
16 36. Durrani, F., M.J. Cook, and J.J. McGuirk, Evaluation of LES and RANS CFD modelling
17 of multiple steady states in natural ventilation. Building and Environment, 2015. 92: p.
18 167-181. DOI: https://doi.org/10.1016/j.buildenv.2015.04.027
19 37. Chen, Z.D. and Y. Li, Buoyancy-driven displacement natural ventilation in a single-zone
20 building with three-level openings. Building and Environment, 2002. 37(3): p. 295-303.
21 DOI: https://doi.org/10.1016/s0360-1323(01)00021-x
22 38. Gong, J. and Y. Li, Solution multiplicity of smoke flows in a simple building. Fire Safety
23 Science, 2008. 9: p. 895-906. DOI: https://doi.org/10.3801/iafss.fss.9-895
24 39. Gong, J. and Y. Li, Smoke flow bifurcation due to opposing buoyancy in two horizontally
25 connected compartments. Fire Safety Journal, 2013. 59: p. 62-75. DOI:
26 https://doi.org/10.1016/j.firesaf.2013.03.014
27 40. Yang, D., et al., Multiple patterns of heat and mass flow induced by the competition of
28 forced longitudinal ventilation and stack effect in sloping tunnels. International Journal of
29 Thermal Sciences, 2019. 138: p. 35-46. https://doi.org/10.1016/j.ijthermalsci.2018.12.018
30 41. Yang, L., et al., Nonlinear Dynamic Aalysis of Natural Ventilation in a Two-Zone
31 Building: Part B—CFD Simulations. HVAC&R Research, 2006. 12(2): p. 257-278. DOI:
32 https://doi.org/10.1080/10789669.2006.10391178
33 42. Yang, L., P. Xu, and Y. Li, Nonlinear dynamic analysis of natural ventilation in a two-
34 zone building: Part A—Theoretical analysis. HVAC&R Research, 2006. 12(2): p. 231-
35 255. DOI: https://doi.org/10.1080/10789669.2006.10391177
36 43. Li, Y., et al., Flow bifurcation due to opposing buoyancy in two vertically connected open
37 cavities. International journal of heat and mass transfer, 2006. 49(19-20): p. 3298-3312.
38 DOI: https://doi.org/10.1016/j.ijheatmasstransfer.2006.03.016
39 44. Yang, D., et al., Multiple steady states of fire smoke transport in a multi-branch tunnel:
40 Theoretical and numerical studies. Tunnelling and Underground Space Technology, 2017.
41 61: p. 189-197. DOI: https://doi.org/10.1016/j.tust.2016.10.009
45
1 45. Liu, Y., et al., The formation of multi-steady-states of buoyancy ventilation in
2 underground building. Tunnelling and Underground Space Technology, 2018. 82: p. 613-
3 626. DOI: https://doi.org/10.1016/j.tust.2018.09.008
4 46. Axley, J., Multizone airflow modeling in buildings: History and theory. HVAC&R
5 Research, 2007. 13(6): p. 907-928. DOI:
6 https://doi.org/10.1080/10789669.2007.10391462
46

Boyouncy

Uploaded by

Copyright:

Available Formats

Boyouncy

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Boyouncy

Uploaded by

Copyright:

Available Formats

1 Nonlinear dynamic analysis of solution multiplicity of

5 a School of Civil Engineering, Chongqing University, Chongqing 400045, China

6 b Key Laboratory of New Technology for Construction of Cities in Mountain Area of

7 Ministry of Education, Chongqing University, Chongqing 400045, China

8 c National Centre for International Research of Low-carbon and Green Buildings,

9 Chongqing University, Chongqing 400045, China

10 d National Renewable Energy Lab, USA

8 Keywords：Nonlinear dynamics; buoyancy ventilation; multiple steady states；

𝑞1 Mass flow rate at zone 1, [kg/s]

𝑞2 Mass flow rate at zone 2, [kg/s]

𝑇1 Air temperature at zone 1, [K]

𝑇𝑎 Outdoor air temperature, [K]

𝑆1 Coefficient of Mass flow impedance at zone1

𝑆2 Coefficient of Mass flow impedance at zone2

𝑀1 Thermal mass at zone 1, [kg]

𝑀2 Thermal mass at zone 2, [kg]

𝐸1 Total heat gain at zone 1 , [kW]

𝐸2 Total heat gain at zone 2 , [kW]

𝑔 Gravitational acceleration, [m/s2]

𝐻1 Height of zone 1,[m]

𝐻2 Height of zone 2,[m]

𝜌0 Ambient air density, [kg/m3]

𝜆 Eigenvalue of coefficient matrix

3 Owing to the wide utilization of underground spaces, the evaluation of natural

4 ventilation in underground buildings has become significant [1-11]. Natural ventilation

5 as a passive strategy has been adopted in many underground structures, such as

10 a type of natural ventilation, buoyancy ventilation in underground structures is worth

16 performed regarding solar-chimney-induced thermal ventilation [13, 14], natural

2 buoyancy and wind pressure.

7 effect of wind pressure and thermal pressure. A one-dimensional mathematical model

15 multiplicity to personnel safety and ventilation design, although its formation

16 mechanism is not elaborated.

17 Subsequently, the solution multiplicity of single-zone and double-opening buildings

24 one-dimensional mathematical model while considering the effect of heat transfer in

12 solutions in CFD simulations.

13 In addition, other studies focusing on the solution multiplicity of ventilation in single-

16 building through a saltwater experiment employing one-dimensional model analyses

27 the initial conditions are different.

8 solution multiplicity and presented a visualization experiment. Yang D. [40] analyzed

11 transient energy balance and pressure balance of a single-zone building to establish a

15 of a single-zone building. Furthermore, salt water experiments were conducted to

18 Additionally, studies regarding the solution multiplicity of ventilation in two-zone

5 states in an underground building with two tunnel connecting to the outdoor

9 by changing the initial conditions. This study primarily provided a method to

10 investigate the multiple steady states.

8 2. Nonlinear dynamic analysis

4 For realization 1, the conservation of mass is as follows:

9 element pressure losses [8, 46].The pressure balance equation is as follows:

16 From Eqs. 2-1 and 2-2, we can obtain

3 nonlinear ordinary differential equation system:

9 Similarly, we can obtain the mass balance equation for realization 2:

13 flow. We can obtain the pressure balance equation:

15 Furthermore, the heat balance equation for realization 2 is as follows:

18 From Eqs. 2-8 and 2-9, we can obtain

21 nonlinear ordinary differential equation system:

5 2.2 Stability and existence of the system

7 0, 𝐶𝑝 = 1, 𝑀1 = 𝑀2 = 1 , the nonlinear ordinary differential equation system can be

9 2.2.1 Stability analysis for scenario 1(𝜅 is fixed, 𝛼 is control parameter)

10 2.2.1.1 Stability analysis for 𝜅 > 0