Nonlinear Dynamic Analysis of SolutionMultiplicity of Buoyancy Ventilation in Two Vertically Connected Open Cavities with Unequal Heights

Solution multiplicity of natural ventilation in buildings is of much importance for personnel safety and ventilation design. In this paper, a new mathematical model of buoyancy pressure ventilation for two vertically connected open cavities is presented. Compared with the previous published papers studying two vertically connected open cavities with equal heights and hot source E2< 0 in the upper room, we study two vertically connected open cavities with unequal heights and hot source E2< 0 or E2> 0 in the upper room. By solving and analyzing the equilibrium points and characteristic roots of the differential equations, we analyze the stability of two systems with upward flow pattern and downward pattern and obtain the criteria to determine the stability and existence of solutions for two scenarios. According to these criteria, the multiple steady states of buoyancy ventilation in two vertically connected open cavities with unequal heights and variable strength of hot sources can be obtained.)ese criteria can be used to design buoyancy ventilation or natural exhaust ventilation systems in two vertically connected open cavities. Compared with two stable states of buoyancy ventilation existing in two vertically connected open cavities with equal heights in the previously published papers, we find that more stable states and unstable states of buoyancy ventilation exist in two vertically connected open cavities with unequal heights in our paper. Finally, bifurcation diagrams and the phase portraits for the two scenarios are given.


Introduction
Nonlinear characteristics exist in systems with various research directions, such as chaotic circuit [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15], neural network [16][17][18][19][20][21], image encryption [22][23][24][25]. Natural ventilation also has nonlinear characteristics. Natural ventilation is increasingly recognized as an important energy saving technology that provides thermal comfort and a healthy indoor environment for buildings [26]. e main objectives of natural ventilation include the provision of adequate outdoor air, effective removal of indoor air, heat, smoke, and/or odour pollution, and comfortable range of air temperature and humidity. In many cases, buoyancy due to internal sources, such as people, heat generating equipment/ processes, or fire, may be properly utilized to facilitate ventilation flow and help achieve these objectives [27]. rough good architectural design and opening arrangement, the ideal airflow can be achieved to avoid the diffusion of pollutants or smoke from the original indoor source into the adjacent rooms/zones. Inappropriate airflow patterns in buildings can endanger human safety and comfort. In modern buildings with many rooms and complex configurations, smoke or pollutants in a portion of the building may pose a hazard to occupants and may migrate from the building's source to nearby or remote locations due to airflow [28]. Especially, in building fires, smoke is a serious hazard. Most deaths in building fires are due to smoke inhalation rather than burns. Fire smoke safety has always been a very important issue in architectural design, and smoke movement and control management have been paid more and more attention [29]. e airflow pattern that drives natural ventilation due to buoyancy is basically the same as that of smoke extraction [30], and similar principles and analytical methods used to predict natural ventilation can also be used to predict smoke or pollutant movement in building fires. For example, CONTAM, a multizone airflow and pollutant transport program, has been used in smoke movement and control design of buildings [29]. Awbi believes that ventilation designers are faced with two major problems: airflow control in space and indoor air movement [31,32]. When pollutants or smoke sources appear in a certain zone of the building, the ventilation direction or ventilation mode of the building plays an important role in controlling and avoiding bad airflow.
It is found that the air flow and temperature of indoor building thermal ventilation behave as nonlinear multiple solution states, which provide a reference for building designers to design building ventilation. Multiplicity of solutions means that under the same boundary conditions, there is more than one solution to the overall flow pattern in the room. Different initial conditions can lead to different solutions [33]. Zheng et al. investigated the nonlinear equations for ventilation flows in an enclosure with a single opening [34]. Yang et al. carried out theoretical and computational nonlinear dynamic analysis on natural ventilation of two-zone buildings [35,36]. Van Schijndel carried out a numerical study on the chaotic behavior of airflow in a normally ventilated room [37]. Xin et al. made a numerical study to reveal some bifurcated sequences that may lead to chaotic mixed convection in finite space [38]. Dillon et al. studied chaotic natural convection in a high rectangular cavity with nonisothermal walls [39]. Yee et al. paid particular attention to the separation of different nonlinear behaviors and spurious dynamics caused by some numerical uncertainties observed in actual CFD calculations [40]. Liu et al. studied the formation process of multiple steady states of the underground building of two tunnels connected with the outdoor environment [41]. e two tunnels are of the same height, and one of the heat sources is located in a corner of the underground space. e CFD method is used to reproduce two steady states of buoyancy ventilation of underground building. Liu et al. conducted a nonlinear dynamic analysis of the multiplicity of buoyancy ventilation solutions for typical underground structures [42]. e mathematical model of second-order nonlinear differential equation for underground structure with two zones and two tunnels connecting to outdoor environment is described.
As a special architectural structure, the cavity can exist in buildings, electronic equipment, etc [43]. Many papers about multiple solutions of air flow and temperature in cavity were published [39,44,45]. However, in the published papers, computational fluid mechanics method is adopted to solve the problem, which is relatively complex. Based on the equations of air flow equilibrium and heat equilibrium, Li established dynamic differential equations about temperature with two-layer open cavities and obtained two stable equilibrium points and two corresponding stable manifolds [46]. e hysteresis and Hopf bifurcation between stable manifolds are studied. However, the literature [46] assumes that the height of two cavities is the same. In practice, the height of the two cavities can be different; for example, the height of the first floor is often greater than that of the second floor in a building. erefore, it is of practical significance and novelty to study the air flow and temperature solutions of buoyancy ventilation in the cavities when the height of the cavity is not equal. In addition, the literature [46] only studies one situation corresponding to hot sources E 1 > 0 and E 2 < 0. In fact, there is another situation corresponding to hot sources E 1 > 0 and E 2 > 0. e situation occurs on fine day in winter, when sun shines into the second floor of the cavity which is equivalent of E 2 > 0. In this paper, we study the two-layer open cavities with unequal heights and hot sources E 1 > 0 and E 2 <0 or E 2 >0, the mathematical model is established, and multiple solutions are obtained. Compared with two stable solutions in literature [46], in our paper, more stable states and unstable states are obtained. e criterion for scenario 1(α is fixed and k is the control parameter) corresponds to solutions of three stable states and one unstable state; and criterion for scenario 2 ((k is fixed, α is the control parameter, and k > 0) is corresponding to solutions of ten stable states and three unstable states. e nomenclature of every variables and constants is shown in Figure 1.

Proposed Mathematical Model of Two Vertically Connected Open Cavities with Unequal Heights
Two vertically connected cavities with a top opening and a bottom opening are considered as shown in Figure 2. All walls and the partition between the two cavities are adiabatic. It is clear that two flow patterns, i. e., an upward flow shown in Figure 2(a) or a downward flow shown in Figure 2(b), are possible depending on the direction of the overall buoyancy force. We assume that the fluid temperature in each cavity is uniform. e magnitude of the heat source in the lower cavity is E 1 (>0), and the heat source in the upper cavity is E 2 (>0 or <0). Compared with [46], the height of two rooms is equal each other and heat source in the upper cavity E 2 is just less than zero, while in this paper, the heights of two rooms is unequal and heat source in the upper cavity E 2 is less than zero or larger than zero.
In Figure 2(a), the heat gain of the internal thermal mass should be equal to the heat released by the heat sources minus the heat loss through airflows. erefore, we can obtain the heat balance equations of the two cavities as shown in the following equations: According to the conservation of mass, we can obtain According to the flow loop method, the total pressure loss at ventilation vents on the flow loop is balanced by sum 2 Complexity of the buoyancy pressure [42,47].
e pressure balance equation can be obtained as follows: 1 where 1/A 12 2 � 1/A 1 2 + 1/A 2 2 and g i ′ � ((T i − T a )/T a )g. From equation (4), it is clear that In Figure 2(b), the heat gain of the internal thermal mass should be equal to the heat released by the heat sources minus the heat loss through airflows. erefore, we can obtain the heat balance equations of the two cavities shown in the following equations: According to the flow loop method, the total pressure loss at ventilation vents on the flow loop is balanced by sum of the buoyancy pressure [42,47]. e pressure balance equation can be obtained as follows:   Cavity 2 Complexity where 1/A 12 2 � 1/A 1 2 + 1/A 2 2 and g i � ((T i − T a )/T a )g. From equation (8), it is clear that Assuming that k � E 2 /E 1 , ������ � H 1 g/T a , and α � H 2 /H 1 , we can obtain simplified nonlinear differential equations as follows.
From Figure 2(a), according to equations (1), (2), and (5) and above assumptions about variables and parameters, the simplified differential equations can be obtained as follows: From Figure 2(b), according to equations (6), (7), and (9) and above assumptions about variables and parameters, the simplified differential equations can be obtained as follows:

Stability Analysis for Scenario 1(α Is Fixed and k Is Control
Parameter). e equilibrium point and the characteristic root should be obtained before the stability analysis. Figure 2(a). e steady state solutions (equilibrium points) of equations (10)∼(11) in Figure 2(a) can be denoted as (ΔT 01 , ΔT 02 ). In order to solve the equilibrium points, the right sides of equations (10)∼(11) can be equal to zero, and thus we can obtain the following equations:

Analysis of Equilibrium Points and Characteristic Roots in
After equation (14) plus equation (15), the following equation can be easily obtained: Substituting equation (16) into equation (14), we can obtain equilibrium points as follows: For convenience, we can make M � 1 + α + αk. To ensure that a real solution exists for equations (17)∼(18) and equilibrium points exist, we have When the characteristic determinant is equal to 0, the solution of the differential equation can be discussed by characteristic roots. e expression of characteristic determinant is shown as follows: Combining equations (10), (11), and (20), the following equation can be obtained: where J(Q) is the Jacobian matrix of equilibrium point Q, λ is the characteristic root, E is the identity matrix, and for convenience, we can make A � ����������� � ΔT 01 + αΔT 02 and B � ΔT 02 − ΔT 01 .
According to equations (17)∼(18), we can obtain the following equations: Combining equations (21)∼(25), we can obtain the characteristic equation of Figure 2(a) as follows: For convenience, we can make 4 Complexity Two characteristic roots of equation (26) can be obtained according to a quadratic formula: Since b > 0 and E 1 > 0, from equations (19) and (28), it is clear that c is always greater than 0.
In summary, according to equation (27), if the following expression is satisfied, λ 1 and λ 2 are all less than 0 and the equilibrium point is stable: If the following expression is satisfied, λ 1 and λ 2 are all greater than 0 and the equilibrium point is unstable:
According to equations (44)∼(45), we can obtain the following equations: Combining equations (47)∼(51), we can obtain characteristic equation of Figure 2(b) as follows: For convenience, we can make Two characteristic roots of equation (52) can be obtained according to a quadratic formula: Since b > 0 and E 1 > 0, from equations (43) and (54), it is clear that c is always greater than 0.
In summary, according to equation (53), if the following expression is satisfied, λ 1 and λ 2 are all less than 0 and the equilibrium point is stable: If the following expression is satisfied, λ 1 and λ 2 are all greater than 0 and the equilibrium point is unstable: Figure 2(b) (α Is Fixed and k Is Control Parameter)
In a word, from the above discussion, we can form a criterion shown in Table 1 to determine the stability of the system shown in Figure 2 when α is fixed and k is the control parameter. For any α > 0, according to the values of k in Table 1, we can evaluate whether there is an equilibrium point and the equilibrium point is stable or unstable in scenario 1(α is fixed and k is the control parameter).

Stability Analysis for Scenario 2 (k Is Fixed and α Is Control Parameter).
e characteristic equation of scenario 2 is the same as that of scenario 1. Figure 2 (1), we know that no equilibrium points exist in Figure 2(a), when the following inequalities (taking out the coefficient of α in expression (32)) are satisfied:
From the above discussion, when K is fixed and α is the control parameter, we can form a criterion shown in Table 2 to determine the stability of the system shown in Figure 2. According to the values of α in Table 2, we can know whether there is an equilibrium point and whether the equilibrium point is stable or not.
Compared with Tables 1-2 in this paper and theoretical analysis results in [44], it can be seen that more stable states and unstable states exist in two vertically connected open cavities with unequal heights in our paper.

Numerical Simulations
To further study the nonlinear dynamical behavior of the proposed model, we solve the differential equation numerically by the forth-order Runge-Kutta method. e area of each room is 3 m (wide) × 3 m (high). e height of lower room is fixed as 4 m. ere are three openings of 0.3 m 2 . One of the middle partitions is used to connect the two rooms, and the top and bottom partitions are used to two connect the two rooms to the outdoors, respectively. Section 4.1 describes the bifurcation diagrams and phase portraits of scenario 1, where the height ratio is fixed and the heat ratio is the control parameter; Section 4.2 describes the bifurcation diagrams and the phase portraits of scenario 2, where the heat strength ratio is fixed and the height ratio is the control parameter.

Fixed Height Ratio α and Different Heat Strength Ratio k.
For scenario 1, we use α � 1 as an example to study the effect of control parameter k on the stability of the system. e bifurcation diagram for ventilation flow rates is plotted in Figure 3(a), and the bifurcation diagram for air temperature is plotted in Figure 3(b). In Figure 3(a), there is only one stable equilibrium point for the flow rate in the zone between A d and B d , where the heat ratio is k < − 2. e air flow in this zone is downward and the temperature in the upper room is less than that in the lower room, as given in Figure 3(b). In the zone between E u and F u , there is also one stable equilibrium point for the flow rate. In the zone between E u and F u , when k is less than zero, there is a heat source in the lower room and a heat sink in the upper room and the air temperature in the upper room is lower than that in the lower room. e air temperature in the upper room equals to that in the lower room when k � 0. e heat sink in the upper room turns to a heat source when k further increases. e flow in the zone between C u and D u has two stable equilibrium points. When the flow rate is positive, the flow is upward and the air temperature in the lower room is higher than that in the upper room. When the flow rate is negative, the flow is downward and the same air temperature trend can be observed in Figure 3(b). e flow in the zones between B u and C u and between D d and E d has one stable equilibrium point and one unstable equilibrium point, respectively. At the point B u , the flow is upward and the temperatures in both rooms are infinite, which indicates that the upward flow is unstable. It can also be given that the upward flow between B u and C u is unstable. With the same principle, an unstable equilibrium point is generated from infinity at E d . e equilibrium point between D d and E d is also unstable. By solving the system equations numerically using the fourth-order Runge-Kutta method, we plot the phase portraits of the model. Based on the bifurcation diagram, six typical values of k are selected, which are − 2.5, − 1.85, − 1.75, − 0.56, − 0.52, and − 0.5. In Figure 4, we give the phase portraits for the six heat source ratios. When k � − 2.5 and k � − 0.5, the orbits from different initial conditions converge to one equilibrium point. When k � − 1.85 and k � − 0.52, a stable equilibrium point and an unstable equilibrium point can be seen. When k � − 1.75 and k � − 0.56, we see two stable equilibrium points and unstable periodic orbits may also be seen in both phase portraits. e numerical simulated phase portraits agree well with the theoretical analysis results as shown in Table 1.

Fixed Heat Strength and Different Height Ratio.
For scenario 2, we use four typical values of k, which are − 2, − 1.2, − 0.9, and − 0.5 to study the effect of control parameter α on the stability of the system. When k is fixed to − 2 and α varies from 0 to 1.5, the bifurcation diagram for the airflow rates and air temperatures are plotted in Figure 5(a) and Figure 5(b), respectively. In Figure 5(a), there are two stable equilibrium points between A u (or A d ) and B u (or B d ). e air temperature in the lower room is higher than that in the upper room when the flow is upward. When the flow is downward, the same trend can be observed. In Figure 5(a), it is shown that an unstable equilibrium point is generated from infinity at the point C u . is is shown in Figure 5(b) where the air temperature in both the lower room and the upper room is infinite. It can also be seen that the upward flow between B u and C u in Figure 5(a) is also unstable. ere is only one stable equilibrium point in the zone between C d and D d for the air flow rate. e air flow is downward. e air temperature in the lower room is higher than that in the upper room.
In Figure 6, we give the phase portraits for α selected as 0.5, 0.85, and 1.5, when k is fixed to − 2. When α � 0.5, we see two stable equilibrium points in the phase portrait. When α � 0.85, there is a stable equilibrium point and an unstable equilibrium point. When α � 1.5, there is only one stable equilibrium point. e numerical simulated phase portraits agree well with the theoretical analysis results as shown in Table 2. Table 2: Stability criterion for scenario 2(k is fixed and α is the control parameter). k α Figure 2(a) Figure 2 10 Complexity e bifurcation diagram for the airflow rates and air temperatures, when k is fixed to − 1.2 and α varies from 0 to 6, is plotted in Figures 7(a) and 7(b), respectively. In Figure 7(a), there are two stable equilibrium points between A u (or A d ) and B u (or B d ). In Figure 7(a), it is shown that an unstable equilibrium point is generated from infinity at the point C u .
is is shown in Figure 7(b) where the air temperature in both the lower room and the upper room is infinite. It can also be seen that the upward flow between B u and C u in Figure 7(a) is also unstable. ere is only one stable equilibrium point in the zone between C d and D d for the air flow rate. e air flow is downward. e air   temperature in the lower room is higher than that in the upper room.
In Figure 8, we give the phase portraits for α selected as 0.01, 1, 3, and 6, when k is fixed as − 1.2. When α � 0.01, we see a stable equilibrium point and an unstable equilibrium point in the phase portrait. When α � 1, there are two stable equilibrium points. When α � 3, there is one stable equilibrium point and one unstable equilibrium point. When α � 6, there is only one stable equilibrium point in the phase portrait. e numerical simulated phase portraits agree well with the theoretical analysis results as shown in Table 2. e bifurcation diagram for the airflow rates and air temperatures, when k is fixed to − 0.9 and α varies from 0 to 12, is plotted in Figures 9(a) and 9(b), respectively. In Figure 9(a), there is only one stable equilibrium point between A u and B u . In Figure 7(a), it is shown that an unstable equilibrium point is generated from infinity at the point B u .
is is shown in Figure 7(b) where the air temperature in both the lower room and the upper room is infinite. It can also be seen that the upward flow between B u and C u in Figure 7(a) is also unstable. In the zone between C u (or C d ) and D u (or D d ), there are two stable equilibrium points.
ere is one stable equilibrium point and one unstable equilibrium point in the zone between D d and E d for the air flow rate.
In Figure 10, we give the phase portraits for α selected as 0.05, 0.2, 1, and 12 and when k is fixed as − 0.9. When α � 0.05, we see only one stable equilibrium point in the phase portrait. When α � 0.2, there is one stable equilibrium point and one unstable equilibrium point. When α � 1, there are two stable equilibrium points. When α � 12, there is one stable equilibrium point and one unstable equilibrium point. e numerical simulated phase portraits agree well with the theoretical analysis as shown in Table 2. Figure 11 shows the bifurcation diagram for the airflow rates and air temperatures, when k is fixed to − 0.5 and α varies from 0 to 2, respectively. In Figure 11(a), there is only one stable equilibrium point between Au and Bu. In      Figure 11(a), it is shown that an unstable equilibrium point is generated from infinity at the point B d . is is shown in Figure 11(b) where the air temperature in both the lower room and the upper room is infinite. It can also be seen that the upward flow between B u and C u in Figure 11(a) is also unstable. In the zone between C u (or

Complexity 13
C d ) and D u (or D d ), there are two stable equilibrium points.
In Figure 12, we give the phase portraits for α selected as 0.5, 1.1, and 2 and when k is fixed to − 0.5. When α � 0.5, we see only one stable equilibrium point in the phase portrait. When α � 1.1 there is one stable equilibrium point and one unstable equilibrium point. When α � 2, there are two stable equilibrium points. e numerical simulated phase portraits agree well with the theoretical analysis results as shown in Table 2.

Conclusions
Nonlinear dynamical analysis was performed to study the solution multiplicity of buoyancy ventilation in two vertically connected open cavities with unequal heights and variable heat source strength. A new mathematical model of multiple solutions of buoyancy pressure ventilation for two vertically connected open cavities with unequal heights and hot source E2 < 0 or E2 > 0 in upper room was presented. By solving and analyzing the equilibrium points and characteristic roots of the differential equations, we analyze the stability of two systems with upward flow pattern and downward pattern and obtain the criteria to determine the stability and existence of solutions for two scenarios. e criteria for the stability and existence of equilibrium points were derived mathematically in detail about two scenarios. e criterion for scenario 1 (α � H 2 / H 1 is fixed and k � E 2 /E 1 is control parameter) is summarized in Table 1. e criterion for scenario 2 (k is fixed and α is control parameter) is summarized in Table 2. Finally, bifurcation diagrams and the phase portraits obtained by applying the fourth-order Runge-Kutta method are given.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest. Acknowledgments e work was supported by the National Natural Science Foundation of China (no. 51378184).