The Numerical Simulation Study of Pumping Airﬂow Driven by Wind Pressure for Single-and Multi-Room Buildings

: Pumping airﬂow occurs in single-sided natural ventilation buildings when the openings are located on the leeward side; the direction of airﬂow through the building then changes periodicity. In order to propose a calculation method of the ventilation rates for single-room buildings and analyze the pumping ventilation for multi-room buildings, the CFD method with an SST k-ω turbulence model was used to conduct numerical simulation in this study, which is veriﬁed by other experimental results. Firstly, the qualitative and quantitative characteristics of pumping ventilation were investigated for a single-room building with two openings. The results show that the steady-state method underestimates the ventilation ﬂow rates, and the unsteady-state method captures the microstructures of the ﬂow better. Secondly, the vortex-shedding and indoor airﬂow oscillation frequencies were analyzed based on transient simulation for a single-room building. It was found that both of them increase with air speed. Then, the factors affecting ventilation ﬂow rates were analyzed. A calculation method for dimensionless ventilation rates is proposed. Finally, the pumping ventilation for a multi-room building was numerically simulated. The ventilation rate for a middle room is greater than that for a corner room, and the ventilation rate for any room in a multi-room building is greater than the single-room building given the same room size and the same incoming wind speed. The ﬁndings of this paper are helpful for the design and evaluation of natural ventilation.


Introduction
Natural ventilation is an important passive technology for buildings, because its appropriate application can conserve energy and supply fresh outdoor air.Natural ventilation can be divided into single-sided ventilation and cross-ventilation.Many studies have shown that cross-ventilation is superior to single-side ventilation in the field of indoor ventilation [1,2], but in most cases, single-side ventilation is still very important in building design, as few buildings can be cross-ventilated due to internal zoning, fire codes, safety requirements and privacy issues.
Single-side ventilation systems can be categorized by their number of openings: a single opening (SS1) and multiple openings (SSn).In the former, air must enter and exit through the same openings, while in the latter, there is a clearer division of air entering and exiting the room [3].It has been found that multiple openings provide a better ventilation effect than a single opening [4].The unilateral natural ventilation model proposed by Warren [5] is suitable for the estimation of ventilation rates in most cases.In single-side ventilation, unsteady natural ventilation mechanisms are classified into continuous airflow, pulsating flow, eddy penetration and diffusion phenomena [6].
The ventilation is induced by the vortex shedding of the airflow around the building, which is referred to as the pumping airflow [7].It can be considered as a special type of eddy penetration.In 2016, Daish et al. [7] first analyzed the pumping airflow of the SS2 model, which has two openings located on the leeward side of the structure.The ventilation airflow in this instance oscillated with the Strohal frequency for time T, exhibiting a unique periodic fluctuation in flow.
In 2018, the pumping airflow was also observed by Zhong [8] in a 2D CFD simulation, and the impact of the building aspect ratio, incoming wind speed, and opening spacing on the frequency and ventilation flow rate of the pumping flow was examined.The two-dimensional simulation method, however, was unable to accurately capture threedimensional flow properties.In 2019, Zhong carried out quasi-3D numerical simulation, and they found that the increasing trend in the dimensionless ventilation flow rate disappears or even decreases as the wind speed increases to 5 m/s.The dimensionless ventilation flow rate increased as the aperture separation increased, but it had no effect on the vortexshedding frequency [9].In 2020, Daniel [10] pointed out that the studies proposed by Zhong et al. [8,9,11,12] did not consider the three-dimensional structure of the pumped ventilation flow in buildings.In 2021, Zhong [13] pointed out that the ventilation flow rate of the bottom and top floors of buildings is considerably different from that of the middle floors, because the ground and top floors of a building have a greater influence on vortex shedding.In 2021, Graça [14] found that there are particular wind angles that also cause pumping airflow, in addition to when the wind is perpendicular to the facade (single-sided rooms) or aligned with the building corner (corner rooms), when there were two openings in the same facade of a building.
In summary, pumping ventilation occurs when the two openings are located on the leeward side for single-room model buildings.The frequency change in the flow's direction is related to the shedding frequency of the vortex street.Moreover, it was found that the pumping ventilation volume was related to the aperture separation: for example, Daish [7] stated that the ventilation flow rate increased with s 1/2 (s is the aperture separation), while Kobayashi [15] and Daniel [10] concluded that the ventilation flow rate was linearly related to s .
There are different conclusions about the relationship.In addition, previous research did not consider a multi-room building, which is the more realistic case.Does pumping ventilation still exist, and what is the difference for rooms in different locations?Aiming to answer these questions and further understand pumping ventilation, this paper analyzed the qualitative and quantitative characteristics of the pumping airflow of a single-room building, and it also investigated the pumping airflow of a multi-room building.Figure 1 shows the flow chart of this paper.
Buildings 2023, 13, x FOR PEER REVIEW 2 of 21 The ventilation is induced by the vortex shedding of the airflow around the building, which is referred to as the pumping airflow [7].In 2016, Daish et al. [7] first analyzed the pumping airflow of the SS2 model, which has two openings located on the leeward side of the structure.The ventilation airflow in this instance oscillated with the Strohal frequency for time T, exhibiting a unique periodic fluctuation in flow.
In 2018, the pumping airflow was also observed by Zhong [8] in a 2D CFD simulation, and the impact of the building aspect ratio, incoming wind speed, and opening spacing on the frequency and ventilation flow rate of the pumping flow was examined.The twodimensional simulation method, however, was unable to accurately capture threedimensional flow properties.In 2019, Zhong carried out quasi-3D numerical simulation, and they found that the increasing trend in the dimensionless ventilation flow rate disappears or even decreases as the wind speed increases to 5 m/s.The dimensionless ventilation flow rate increased as the aperture separation increased, but it had no effect on the vortex-shedding frequency [9].In 2020, Daniel [10] pointed out that the studies proposed by Zhong et al. [8,9,11,12] did not consider the three-dimensional structure of the pumped ventilation flow in buildings.In 2021, Zhong [13] pointed out that the ventilation flow rate of the bottom and top floors of buildings is considerably different from that of the middle floors, because the ground and top floors of a building have a greater influence on vortex shedding.In 2021, Graça [14] found that there are particular wind angles that also cause pumping airflow, in addition to when the wind is perpendicular to the facade (single-sided rooms) or aligned with the building corner (corner rooms), when there were two openings in the same facade of a building.
In summary, pumping ventilation occurs when the two openings are located on the leeward side for single-room model buildings.The frequency change in the flow's direction is related to the shedding frequency of the vortex street.Moreover, it was found that the pumping ventilation volume was related to the aperture separation: for example, Daish [7] stated that the ventilation flow rate increased with  ′1 2  ⁄ ( ′ is the aperture separation), while Kobayashi [15] and Daniel [10] concluded that the ventilation flow rate was linearly related to  ′ .There are different conclusions about the relationship.In addition, previous research did not consider a multi-room building, which is the more realistic case.Does pumping ventilation still exist, and what is the difference for rooms in different locations?Aiming to answer these questions and further understand pumping ventilation, this paper analyzed the qualitative and quantitative characteristics of the pumping airflow of a single-room building, and it also investigated the pumping airflow of a multi-room building.Figure 1 shows the flow chart of this paper.Revealing the principles of the natural ventilation is a prerequisite for further designing the thermal environment and evaluating the diffusion of pollutants.This study can provide references for understanding the mechanisms and effects of the natural ventilation.Revealing the principles of the natural ventilation is a prerequisite for further designing the thermal environment and evaluating the diffusion of pollutants.This study can provide references for understanding the mechanisms and effects of the natural ventilation.The computational domain has a significant impact on the simulation results, and the blocking ratio is thought to be important [16], which meant that it should be less than 3% of the computational domain size.Additionally, for the building body, the computational domain's sides, the top, and the external exit boundary all need to be at least 5 H and 10 H away from the building, respectively.The computational domain is 186 m × 96 m × 9 m in size because only the quasi-3D building case was taken into account in this study, which meant that the height of the computational domain was equal to the height of the building, the side distance is 5 H, and the distance from the leeward side to the exit is 15 H.The blockage ratios of it are less than 3%.Due to the fact that the pumping flow involved two apertures and that the spacing between the two openings has an impact on the airflow, three cases with varying spacing are adopted in this paper (e.g., Figure 3).The quoted dimensionless aperture distance is

Model Simplification and Simulation Settings
where s is the distance (in meters) between the centers of the two openings, and B is the building's width (in meters), using 0.25, 0.5, and 0.75 as examples.The computational domain has a significant impact on the simulation results, and the blocking ratio is thought to be important [16], which meant that it should be less than 3% of the computational domain size.Additionally, for the building body, the computational domain's sides, the top, and the external exit boundary all need to be at least 5 H and 10 H away from the building, respectively.The computational domain is 186 m × 96 m × 9 m in size because only the quasi-3D building case was taken into account in this study, which meant that the height of the computational domain was equal to the height of the building, the side distance is 5 H, and the distance from the leeward side to the exit is 15 H.The blockage ratios of it are less than 3%.Due to the fact that the pumping flow involved two apertures and that the spacing between the two openings has an impact on the airflow, three cases with varying spacing are adopted in this paper (e.g., Figure 3).The quoted dimensionless aperture distance is where s is the distance (in meters) between the centers of the two openings, and B is the building's width (in meters), using 0.25, 0.5, and 0.75 as examples.In this paper, the finite volume method was used to solve the governing equation.The governing equation was as follows: (2) The three-dimensional RANS equations were solved together with the k-ω shear stress transport (SST k-ω) model: This paper assumed that the external wind speed was stable and the supply air parameters were constant.In addition, the air in the fluid region was considered to be an incompressible ideal gas.The SIMPLE method was chosen to the pressure-velocity coupling, and the PRESTO!pressure interpolation format was chosen.The controlling equations were discretized as the second-order windward terms because second-order discretization has a higher precision than the first-order discretization.With the exception of the energy equation residuals, which were 1 × 10 −6 , the convergence criterion for the system of equations was 1 × 10 −9 for all parameter residuals.The steady-state simulation was carried out firstly, and then the results were adopted as the initial conditions in the transient calculation.
For transient calculations, a time step which is too large or too small will affect the results.The CFL number can be used to determine the appropriate time step. =  ×   (7) where l is the distance from wall to the first grid, and the value of that is 0.01 m in this study.The CFL value for the explicit solver is close to 1, and for the implicit solver, it will be increased to a higher number [17].The research [18] showed that the smaller time-step difference had a limited impact on the outcomes when the time-step length was

Numerical Simulation Settings Solution Settings
In this paper, the finite volume method was used to solve the governing equation.The governing equation was as follows: The three-dimensional RANS equations were solved together with the k-ω shear stress transport (SST k-ω) model: ∂ρω ∂t This paper assumed that the external wind speed was stable and the supply air parameters were constant.In addition, the air in the fluid region was considered to be an incompressible ideal gas.The SIMPLE method was chosen to the pressure-velocity coupling, and the PRESTO!pressure interpolation format was chosen.The controlling equations were discretized as the second-order windward terms because second-order discretization has a higher precision than the first-order discretization.With the exception of the energy equation residuals, which were 1 × 10 −6 , the convergence criterion for the system of equations was 1 × 10 −9 for all parameter residuals.The steady-state simulation was carried out firstly, and then the results were adopted as the initial conditions in the transient calculation.
For transient calculations, a time step which is too large or too small will affect the results.The CFL number can be used to determine the appropriate time step.
where l is the distance from wall to the first grid, and the value of that is 0.01 m in this study.The CFL value for the explicit solver is close to 1, and for the implicit solver, it will be increased to a higher number [17].The research [18] showed that the smaller time-step difference had a limited impact on the outcomes when the time-step length was substantially less than the vortex-shedding duration; hence, CFL may be considered as 10 [9].When the incoming wind speed is 1 m/s, ∆t is calculated as 0.1 s.Table 1 displays the time step changing with wind speed.A total of 10~20 iterations were run in each time step to ensure accurate calculation results.

Boundary Conditions
The velocity inlet boundary condition was adopted to the inlet, with U y values of 1 m/s, 2 m/s, 3 m/s, 4 m/s, 5 m/s, and the velocity in the other two directions was 0; the turbulence intensity at the inlet was set to 10%, and the turbulence integration length was set to 0.58 m [19].The outlet was set as outflow; the sides, the top and the bottom were set as symmetry; the opening was set as the interior, and all the other wall surfaces were set as no-slip.

Grid Division and Independence Verification
The domain was constructed using the hexahedral gird (Figure 4).The grid was refined in the opening and the near-wall.For the case of s = 0.5, three grid schemes with the grid numbers of 2,079,003, 3,167,446, and 4,851,384 were selected for independence verification.The corresponding y+ values for the three-grid method were 45, 30 and 15, and the minimum distances of the grids from the wall were 0.015 m, 0.01 m, and 0.005 m, respectively.The monitoring line Y = 3 m, and Y = 9 m in the plane with a Z value of 4.5 m is chosen for the sample set (Figure 5).The results with different numbers of grids are shown in Figure 6.As can be seen, the changing rules of the wind speed outside the building showed the same trend.Based on the results, the grid number of 3,167,446 (Figure 4) is finally selected for the case, as a further increase in the mesh density does not produce obvious changes in the numerical solutions.
Buildings 2023, 13, x FOR PEER REVIEW 5 of 21 substantially less than the vortex-shedding duration; hence, CFL may be considered as 10 [9].When the incoming wind speed is 1 m/s, ∆ is calculated as 0.1 s.Table 1 displays the time step changing with wind speed.A total of 10~20 iterations were run in each time step to ensure accurate calculation results.

Boundary Conditions
The velocity inlet boundary condition was adopted to the inlet, with Uy values of 1 m/s, 2 m/s, 3 m/s, 4 m/s, 5 m/s, and the velocity in the other two directions was 0; the turbulence intensity at the inlet was set to 10%, and the turbulence integration length was set to 0.58 m [19].The outlet was set as outflow; the sides, the top and the bottom were set as symmetry; the opening was set as the interior, and all the other wall surfaces were set as no-slip.

Grid Division and Independence Verification
The domain is constructed using the hexahedral gird (Figure 4).The grid is refined in the opening and the near-wall.For the case of s′ = 0.5, three grid schemes with the grid numbers of 2,079,003, 3,167,446, and 4,851,384 were selected for independence verification.The corresponding y+ values for the three-grid method were 45, 30 and 15, and the minimum distances of the grids from the wall were 0.015 m, 0.01 m, and 0.005 m, respectively.The monitoring line Y = 3 m, and Y = 9 m in the plane with a Z value of 4.5 m is chosen for the sample set (Figure 5).The results with different numbers of grids are shown in Figure 6.As can be seen, the changing rules of the wind speed outside the building showed the same trend.Based on the results, the grid number of 3,167,446 (Figure 4) is finally selected for the case, as a further increase in the mesh density does not produce obvious changes in the numerical solutions.

Validation of CFD Methodology
In 2003, Jiang et al. [20] conducted a wind tunnel experiment for natural ventilation.The study focused on a small (0.25 cubic) building model with a window opening size 0.084 m × 0.125 m (length × height), as shown in Figure 7.There is a workspace with the area of 2.0 m × 2.0 m, and the height 1.0 m, which is in the wind tunnel.The velocity measurements were conducted on a grid plane at the center of the building.In order to verify the CFD method used in this paper, Jiang's case was chosen for the numerical simulation.The downstream length was set as 8 H; the upstream length was 4 H; and the lateral length was 4 H on both sides of the building in the computational zone.The total number of the grids was 969,570.The boundary conditions and the wind speed were consistent with the literature, and the turbulence model was the SST k-ω.The numerical results are shown to have good agreement with the experimental data, and the relative error is less than 20% in the most regions (Figure 8).Therefore, the CFD method adopted in the paper was validated for the numerical simulation of the natural ventilation.

Validation of CFD Methodology
In 2003, Jiang et al. [20] conducted a wind tunnel experiment for natural ventilation.The study focused on a small (0.25 cubic) building model with a window opening size 0.084 m × 0.125 m (length × height), as shown in Figure 7.There is a workspace with the area of 2.0 m × 2.0 m, and the height 1.0 m, which is in the wind tunnel.The velocity measurements were conducted on a grid plane at the center of the building.In order to verify the CFD method used in this paper, Jiang's case was chosen for the numerical simulation.The downstream length was set as 8 H; the upstream length was 4 H; and the lateral length was 4 H on both sides of the building in the computational zone.The total number of the grids was 969,570.The boundary conditions and the wind speed were consistent with the literature, and the turbulence model was the SST k-ω.The numerical results are shown to have good agreement with the experimental data, and the relative error is less than 20% in the most regions (Figure 8).Therefore, the CFD method adopted in the paper was validated for the numerical simulation of the natural ventilation.The study focused on a small (0.25 cubic) building model with a window opening size 0.084 m × 0.125 m (length × height), as shown in Figure 7.There is a workspace with the area of 2.0 m × 2.0 m, and the height 1.0 m, which is in the wind tunnel.The velocity measurements were conducted on a grid plane at the center of the building.In order to verify the CFD method used in this paper, Jiang's case was chosen for the numerical simulation.The downstream length was set as 8 H; the upstream length was 4 H; and the lateral length was 4 H on both sides of the building in the computational zone.The total number of the grids was 969,570.The boundary conditions and the wind speed were consistent with the literature, and the turbulence model was the SST k-ω.The numerical results are shown to have good agreement with the experimental data, and the relative error is less than 20% in the most regions (Figure 8).Therefore, the CFD method adopted in the paper was validated for the numerical simulation of the natural ventilation.

The Flow Field
The transient simulation can describe the small structures in turbulent better than steady-state calculation.For example, as shown in Figure 9 (incoming wind speed was 1

Comparison of Steady-State and Transient Simulation 2.2.1. The Flow Field
The transient simulation can describe the small structures in turbulent better than steady-state calculation.For example, as shown in Figure 9 (incoming wind speed was 1 m/s), there was an obvious vortex street in the building wake in the transient simulation compared to the steady-state results.Figure 10 depicts the y-direction velocity cloud charts at the openings on the lee side of the room, where the steady-state result did not have a symmetric image an velocity distribution was rather uniform.The unsteady-state calculation revealed tha velocity distribution was symmetric about z = 4.5 m, where the minimum velocity located at one of the openings and the maximum velocity was located near the wa Figure 10 depicts the y-direction velocity cloud charts at the openings on the leeward side of the room, where the steady-state result did not have a symmetric image and the velocity distribution was rather uniform.The unsteady-state calculation revealed that the velocity distribution was symmetric about z = 4.5 m, where the minimum velocity was located at one of the openings and the maximum velocity was located near the wall.In general, there is a significant difference in the calculation results between the transient calculation and the steady-state calculation because the latter result is the average value.Figure 10 depicts the y-direction velocity cloud charts at the openings on the leew side of the room, where the steady-state result did not have a symmetric image and velocity distribution was rather uniform.The unsteady-state calculation revealed tha velocity distribution was symmetric about z = 4.5 m, where the minimum velocity located at one of the openings and the maximum velocity was located near the wal general, there is a significant difference in the calculation results between the trans calculation and the steady-state calculation because the latter result is the average va

The Ventilation Flow Rates
The ventilation flow rates are also compared here (Table 2).The steady-state results could be obtained directly; however, the ventilation flow rate is the average o total ventilation volume for each time step with the entire period in the non-statio calculation.It can be seen that the steady-state calculation underestimates the ventilation rates of the building.The transient simulation method is more accurate because it t into account the fluctuation of the wind around the buildings.Therefore, the metho the unsteady-state simulation was selected in the following calculation.

The Ventilation Flow Rates
The ventilation flow rates are also compared here (Table 2).The steady-state case results could be obtained directly; however, the ventilation flow rate is the average of the total ventilation volume for each time step with the entire period in the non-stationary calculation.It can be seen that the steady-state calculation underestimates the ventilation flow rates of the building.The transient simulation method is more accurate because it takes into account the fluctuation of the wind around the buildings.Therefore, the method of the unsteady-state simulation was selected in the following calculation.

Periodic Characteristics of Pumping Airflow
The detailed flow field and the periodic characteristics of the typical pumping airflow are analyzed in this section.The case is the same as that in Section 2.2.1.Figure 11

Periodic Characteristics of Pumping Airflow
The detailed flow field and the periodic characteristics of the typical pumping airflow are analyzed in this section.The case is the same as that in Section 2.2.1.Figure 11  The flow around buildings in a natural wind field is similar to the flow around a blunt body.As shown in Figure 11, the incoming air flow separated on the building's windward surface and flowed to its two sides.As a result of the pressure difference, a clockwise vortex formed on the left side of the surface.Next, a separation was created at the backward corner The flow around buildings in a natural wind field is similar to the flow around a blunt body.As shown in Figure 11, the incoming air flow separated on the building's windward surface and flowed to its two sides.As a result of the pressure difference, a clockwise vortex formed on the left side of the surface.Next, a separation was created at the backward corner of the left surface, and the vortex gradually grew to the entire leeward surface.
The airflow at the back of the building was forced to enter the building through the opening on the right side, create a circulating airflow that followed the vortex's direction, and exit through the opening on the left.The airflow on the opposite side, however, was exactly the opposite.It formed an anti-clockwise vortex on the right side surface.This caused the airflow to enter the room through the left side opening after separating and attaching itself to the leeward surface at the rear corner point of the lower right side.It then formed a circulating airflow that flowed out of the right side opening in the same direction as the vortex in an anti-clockwise direction.
Vortices were produced and shed alternatively on the either side of the square building, finally forming a regular vortex street in the wake flow.As a result, the fluctuating ventilation named as the pumping airflow was generated.Figure 12 describes the ventilation flow rates of the building changing with time.During the first 200 s, the changing amplitude of the ventilation flow rates was increasing gradually because the vortex shedding of the airflow was still in the developmental stage.After that, the change in amplitude decreased to form the regular fluctuations similar to the sinusoidal form.From Figure 13, it can be discovered that not all of the ventilation flow rates change with the sine wave, especially when the opening spacing was adjusted.However, all the cases clearly displayed the features of a regular fluctuation.The flow around buildings in a natural wind field is similar to the flow around a blunt body.As shown in Figure 11, the incoming air flow separated on the building's windward surface and flowed to its two sides.As a result of the pressure difference, a clockwise vortex formed on the left side of the surface.Next, a separation was created at the backward corner of the left surface, and the vortex gradually grew to the entire leeward surface.
The airflow at the back of the building was forced to enter the building through the opening on the right side, create a circulating airflow that followed the vortex's direction, and exit through the opening on the left.The airflow on the opposite side, however, was exactly the opposite.It formed an anti-clockwise vortex on the right side surface.This caused the airflow to enter the room through the left side opening after separating and attaching itself to the leeward surface at the rear corner point of the lower right side.It then formed a circulating airflow that flowed out of the right side opening in the same direction as the vortex in an anti-clockwise direction.
Vortices were produced and shed alternatively on the either side of the square building, finally forming a regular vortex street in the wake flow.As a result, the fluctuating ventilation named as the pumping airflow was generated.Figure 12 describes the ventilation flow rates of the building changing with time.During the first 200 s, the changing amplitude of the ventilation flow rates was increasing gradually because the vortex shedding of the airflow was still in the developmental stage.After that, the change in amplitude decreased to form the regular fluctuations similar to the sinusoidal form.From Figure 13, it can be discovered that not all of the ventilation flow rates change with the sine wave, especially when the opening spacing was adjusted.However, all the cases clearly displayed the features of a regular fluctuation.

Pumping Airflow in a Single-Room Building
In this part, the numerical simulation was used to examine the features of the pumping airflow with different incoming wind speed in a single-room building.The vortex-shedding frequency and the ventilation flow rates were analyzed.

Fluctuating Frequency
Vortex-Shedding Frequency

Pumping Airflow in a Single-Room Building
In this part, the numerical simulation was used to examine the features of the pumping airflow with different incoming wind speed in a single-room building.The vortex-shedding frequency and the ventilation flow rates were analyzed.

Fluctuating Frequency Vortex-Shedding Frequency
Four monitoring points were selected in order to analyze the vortex-shedding frequency, as shown in Figure 14.Point b was located at the center of the opening 1 with a height of 4.5 m.The time-domain and the frequency spectrum plots of the static pressure at point a with the incoming wind speed as 4 m/s are presented in Figure 15a,b, respectively.As shown in the time-domain plot, the static pressure changed periodically.From the frequency spectrum plot, it was possible to determine that the shedding frequency was 0.0825 Hz, at which the power spectral density achieves its maximum value.

Pumping Airflow in a Single-Room Building
In this part, the numerical simulation was used to examine the features of the pumping airflow with different incoming wind speed in a single-room building.The vortex-shedding frequency and the ventilation flow rates were analyzed.

Fluctuating Frequency
Vortex-Shedding Frequency Four monitoring points were selected in order to analyze the vortex-shedding frequency, as shown in Figure 14.Point b was located at the center of the opening 1 with a height of 4.5 m.The time-domain and the frequency spectrum plots of the static pressure at point a with the incoming wind speed as 4 m/s are presented in Figure 15a,b, respectively.As shown in the time-domain plot, the static pressure changed periodically.From the frequency spectrum plot, it was possible to determine that the shedding frequency was 0.0825 Hz, at which the power spectral density achieves its maximum value.

Pumping Airflow in a Single-Room Building
In this part, the numerical simulation was used to examine the features of the pumping airflow with different incoming wind speed in a single-room building.The vortex-shedding frequency and the ventilation flow rates were analyzed.

Fluctuating Frequency
Vortex-Shedding Frequency Four monitoring points were selected in order to analyze the vortex-shedding frequency, as shown in Figure 14.Point b was located at the center of the opening 1 with a height of 4.5 m.The time-domain and the frequency spectrum plots of the static pressure at point a with the incoming wind speed as 4 m/s are presented in Figure 15a,b, respectively.As shown in the time-domain plot, the static pressure changed periodically.From the frequency spectrum plot, it was possible to determine that the shedding frequency was 0.0825 Hz, at which the power spectral density achieves its maximum value.The vortex-shedding frequencies of the four points with different wind speed are presented in Table 3.It was evident from the simulated data that with the same wind speed, the vortex-shedding frequencies essentially stayed constant at different places.The frequency at point b increased by 12.5% more than the other three places only at the case wind speed of 5 m/s.For various wind speeds, the frequency of vortex shedding at the same place increased as the wind speed rose.As a dimensionless number used to describe fluid flow, the Reynolds number (Re) represents the ratio of inertial forces to viscous or frictional forces acting on a moving object.
In the research of the wind field around buildings, the dimensionless number Strauhl number (St) was also utilized to investigate the flow frequencies.
where d was the characteristic length of the object (m); µ was the dynamic viscosity (N•m/s 2 ); U was the velocity (m/s); and f was the frequency (Hz).
Due to the fact that the relationship between the Reynolds number and Strauhl number in the vortex-shedding frequency research was not linear, Roshko proposed another dimensionless number, Roshko number (Ro) [21]: where ν was the kinematic viscosity (m 2 /s).Previous studies had given the formulae for the corresponding Ro numbers as follows [10,22]: where Equation (11) was obtained in the flow field across a cylinder with a low Re number, while Equation ( 12) was obtained for the 3D square cylinder model.Based on the numerical simulation results and combined with the Fourier transform, f and St were obtained.Using Equations ( 8) and ( 10), the Ro number and Re number were determined, and the results are provided in Table 4. Daniel [10] discovered that the St values measured during vortex shedding in the 2D model case ranged from 0.07 to 0.21.In this study, the values of St in the quasi-3D model fell within the range of 0.12 to 0.14 and in a bandwidth of roughly 0.015.The result is relatively reasonable compared to the research of Daniel [10].The following equation was developed by fitting the Re and Ro data using the least squares method, as shown in Figure 16. 4 × 10 5 < Re < 2.1 × 10 6 Ro = 0.098Re + 4207.2 (13)

Indoor Airflow Oscillation Frequency
During the simulations, the instantaneous ventilation flow rates of two openings were monitored.The results with wind speed as 1 m/s are shown in Figure 17, that positive ventilation values indicated the air flow into the room and negative values showed the ai flow out.It could be seen that the ventilation flow rates of two openings were the same but the flow direction was opposite, and all of their sinusoidal properties changed with time.

Indoor Airflow Oscillation Frequency
During the simulations, the instantaneous ventilation flow rates of two openings were monitored.The results with wind speed as 1 m/s are shown in Figure 17, that positive ventilation values indicated the air flow into the room and negative values showed the air flow out.It could be seen that the ventilation flow rates of two openings were the same but the flow direction was opposite, and all of their sinusoidal properties changed with time.Indoor Airflow Oscillation Frequency During the simulations, the instantaneous ventilation flow rates of two openings were monitored.The results with wind speed as 1 m/s are shown in Figure 17, that positive ventilation values indicated the air flow into the room and negative values showed the air flow out.It could be seen that the ventilation flow rates of two openings were the same but the flow direction was opposite, and all of their sinusoidal properties changed with time.Figure 18 depicts the oscillation frequency of the airflow with various wind speeds, showing there was a significant association between the oscillation frequency of the airflow and the wind speed.As the wind speed increased, the frequency of the oscillations inside the building rose.The airflow oscillation frequency was basically the same as the vortexshedding frequency around the building; all of them varied monotonically with the wind speed, as shown in Table 3 and Figure 18.Moreover, if the dimensionless opening distances s changed from 0.25 to 0.5, the simulation results showed that the frequency was only different by about 5%.It could be concluded that the opening distance had almost no impact on the oscillation frequency of the indoor airflow.
airflow and the wind speed.As the wind speed increased, the frequency of the oscillations inside the building rose.The airflow oscillation frequency was basically the same as the vortex-shedding frequency around the building; all of them varied monotonically with the wind speed, as shown in Table 3 and Figure 18.Moreover, if the dimensionless opening distances s′ changed from 0.25 to 0.5, the simulation results showed that the frequency was only different by about 5%.It could be concluded that the opening distance had almost no impact on the oscillation frequency of the indoor airflow.

Ventilation Flow Rates
In the analysis of natural ventilation, the most important parameter is the ventilation flow rates.In this study, the dimensionless ventilation flow rates could be obtained by the following equation.

𝑄𝑄′ =
(14) where UR is the reference wind speed (m/s), which is considered equal to the incoming wind speed in this study; Aref is the reference opening area (m 2 ).
The reference opening area is the effective opening area, which can be calculated by the formula: where A1 and A2 are the areas (m 2 ) of the two openings.In this study, A1 = A2 = A, so   = /√2.According to the literature [23], the mean and fluctuating ventilation made up the instantaneous ventilation.The fluctuating ventilation rates can be predicted by the standard deviation as follows: where σQ is the fluctuating ventilation (m 3 /s), Nm is the time step, Q is the monitored ventilation (m 3 /s), and   is the mean ventilation (m 3 /s), which is the average of the sequence of instantaneous ventilation changes over time.In addition, ′ and σQ′ represent the dimensionless numbers of the mean and fluctuating ventilation, respectively.

The Influence of Different Wind Speeds
The ventilation flow rates with different wind speeds are shown in Table 5.From the results, it can be seen that the dimensionless mean ventilation ′ is larger than the dimensionless fluctuating ventilation σQ′ for all the five cases.It may indicate that

Ventilation Flow Rates
In the analysis of natural ventilation, the most important parameter is the ventilation flow rates.In this study, the dimensionless ventilation flow rates could be obtained by the following equation.
where U R is the reference wind speed (m/s), which is considered equal to the incoming wind speed in this study; A ref is the reference opening area (m 2 ).The reference opening area is the effective opening area, which can be calculated by the formula: where A 1 and A 2 are the areas (m 2 ) of the two openings.In this study, According to the literature [23], the mean and fluctuating ventilation made up the instantaneous ventilation.The fluctuating ventilation rates can be predicted by the standard deviation as follows: where σ Q is the fluctuating ventilation (m 3 /s), N m is the time step, Q is the monitored ventilation (m 3 /s), and Q m is the mean ventilation (m 3 /s), which is the average of the sequence of instantaneous ventilation changes over time.In addition, Q and σ Q represent the dimensionless numbers of the mean and fluctuating ventilation, respectively.

The Influence of Different Wind Speeds
The ventilation flow rates with different wind speeds are shown in Table 5.From the results, it can be seen that the dimensionless mean ventilation Q is larger than the dimensionless fluctuating ventilation σ Q for all the five cases.It may indicate that although there are periodic changes in the direction of ventilation airflow, the average ventilation intensity is greater than the fluctuating ventilation in terms of the overall pumping ventilation effect on the room.Moreover, the simulated results show that almost all the ventilation flow rates increased with the wind speed except the dimensionless fluctuating air volume under the condition of wind speed as 1 m/s.The probable reason is that there was more airflow (which on the vortex scale was smaller than the opening size) entering or exiting the room from the openings at lower wind speeds, which increased the dimensionless fluctuating air volume.Based on the ventilation flow rates, incoming wind speed, and the effective opening area, the following formula can be obtained: The plot of the relationship after the data fitting and the Pearson's correlation coefficient are shown in Figure 19, and its Pearson's correlation coefficient (R 2 ) is 0.993.It indicates that there is a strong linear correlation between the ventilation flow rates and the wind speeds.
although there are periodic changes in the direction of ventilation airflow, the average ventilation intensity is greater than the fluctuating ventilation in terms of the overall pumping ventilation effect on the room.Moreover, the simulated results show that almost all the ventilation flow rates increased with the wind speed except the dimensionless fluctuating air volume under the condition of wind speed as 1 m/s.The probable reason is that there was more airflow (which on the vortex scale was smaller than the opening size) entering or exiting the room from the openings at lower wind speeds, which increased the dimensionless fluctuating air volume.Based on the ventilation flow rates, incoming wind speed, and the effective opening area, the following formula can be obtained: The plot of the relationship after the data fitting and the Pearson's correlation coefficient are shown in Figure 19, and its Pearson's correlation coefficient (R 2 ) is 0.993.It indicates that there is a strong linear correlation between the ventilation flow rates and the wind speeds Warren [5] previously proposed a unilateral natural ventilation model, which was considered the most accurate model suitable for single-sided ventilation [5,24].In addition, Daish [7] carried out special research about the pumping ventilation.Therefore, the two models were applied here to compare the simulated results.Using the simulation results as the base value, the relative error was calculated as: where ε is the relative error (%); Qm is the ventilation flow rates volume obtained by existing model (m 3 /s); and Qn is the numerically simulated ventilation flow rate (m 3 /s).The comparison results are shown in Table 6.The model of the Warren underestimated the ventilation flow rates to a greater extent with the average error about 70%.The results of QDaish were closer to the numerical data than the Warren model, with the relative error smaller than 15%, except for the wind speed of 1 m/s.It indicated that the Warren model Warren [5] previously proposed a unilateral natural ventilation model, which was considered the most accurate model suitable for single-sided ventilation [5,24].In addition, Daish [7] carried out special research about the pumping ventilation.Therefore, the two models were applied here to compare the simulated results.Using the simulation results as the base value, the relative error was calculated as: where ε is the relative error (%); Q m is the ventilation flow rates volume obtained by existing model (m 3 /s); and Q n is the numerically simulated ventilation flow rate (m 3 /s).
The comparison results are shown in Table 6.The model of the Warren underestimated the ventilation flow rates to a greater extent with the average error about 70%.The results of Q Daish were closer to the numerical data than the Warren model, with the relative error smaller than 15%, except for the wind speed of 1 m/s.It indicated that the Warren model was not suitable for the case of the pumping airflow.Comparing the results of Equation (17) and the model of Daish, as shown in Table 7, it could be found that the fitted function drawn in this study had a better ability to predict the pumping ventilation flow rates since there was only about a 10% relative error between the calculated results of the function obtained after fitting and the results measured by Daish [7].The aforementioned data showed that the model that only considered the windward or the parallel wind was no longer applicable for the leeward pumping flow, and it was possible to assume that the fitted function had a better ability to predict the pumping ventilation.

The Influence of Different Opening Spacings
Simulations were carried out for three cases with different opening spacings under the condition of a wind speed of 1 m/s.The dimensionless ventilation for the simulation results is given in Table 8 along with data from related studies in the literature [10,14].The correlation between the two variables of the three sets of data was calculated to obtain the corresponding Pearson's correlation coefficient (R 2 ) to determine the pattern of change between the dimensionless ventilation and the dimensionless opening distance.The results demonstrated that the data in this study and those in Danie [10] followed the same general pattern of law: namely, that the dimensionless ventilation Q was linear with the dimensionless aperture distance s , and the Pearson's correlation coefficients (R 2 ) were both 0.99, which indicated a strong correlation.Additionally, Daish had demonstrated through an experiment that in contrast to the other three listed above, the dimensionless ventilation Q was linearly proportional to the dimensionless aperture distance s .
It should be mentioned that the simulation model or experimental subjects used in the table were different, which may account for the conclusions' discrepancy.Zhong [8] chose a two-dimensional model for the opening distance problem, whereas this paper used a quasi-3D model.The hanging window's opening angle was 14 • in the study of Danie [10], whereas it was fully open in this work.However, all studies explained that the ventilation flow rates and the opening distances have a linear relationship.
Based on the simulation data, the following equation can be received.
where U R is the reference wind speed (m/s); A eff is the effective opening area (m 2 ); S is the opening spacing (m); and B is the building width (m).
The numerical simulation results were compared with the predicted values of the fitted function shown in Equation (19).The results showed that the formula can estimate the ventilation flow rates suitably under different opening sizes with the average relative error of 6.27%.Given a constant building size, the ventilation flow rates will increase with the opening sizes.

Pumping Airflow in Multi-Room Buildings
The pumping airflow of a multi-room building is analyzed in this section.The building model 2 is shown in Figure 20 with three rooms on each floor.There are a total of nine separate rooms, with a size of 18 m × 6 m × 9 m, and two openings on the leeward.The size of opening is 1 m × 1 m, and the distance between their center point is 3 m.The numerical simulation results were compared with the predicted values of the fitted function shown in Equation (19).The results showed that the formula can estimate the ventilation flow rates suitably under different opening sizes with the average relative error of 6.27%.Given a constant building size, the ventilation flow rates will increase with the opening sizes.

Pumping Airflow in Multi-Room Buildings
The pumping airflow of a multi-room building is analyzed in this section.The building model 2 is shown in Figure 20 with three rooms on each floor.There are a total of nine separate rooms, with a size of 18 m × 6 m × 9 m, and two openings on the leeward.The size of opening is 1 m × 1 m, and the distance between their center point is 3 m.
The size of the external computational domain is 186 m × 108 m × 9 m.The gridding, boundary conditions, turbulence model, discrete formatting, and parameter settings for simulation were identical to those in the aforementioned sections.The wind speed is assumed to be 1 m/s.
Considering the symmetry of the building model, four rooms were analyzed.The openings of the rooms were colored yellow, green, red, and blue, respectively, and these are called room 1, room 2, room 3, and room 4.  The size of the external computational domain is 186 m × 108 m × 9 m.The gridding, boundary conditions, turbulence model, discrete formatting, and parameter settings for simulation were identical to those in the aforementioned sections.The wind speed is assumed to be 1 m/s.
Considering the symmetry of the building model, four rooms were analyzed.The openings of the rooms were colored yellow, green, red, and blue, respectively, and these are called room 1, room 2, room 3, and room 4.

The Flow Characteristics
The simulation calculation settings in this section are the same as those in Section 2. The pumping flow characteristics are the same as well.However, the simulation results showed that the indoor airflow oscillation frequency was smaller than that of building 1 (0.0224 Hz); their values were 0.0073 Hz, 0.0075 Hz, 0.0075 Hz, and 0.0079 Hz, respectively.
The exterior flow fields of the single-room and multi-room building are depicted in Figure 21.It was possible to see that with the same incoming wind speed, the building body's expansion caused a gradient increase in the velocity.In a result, a larger vortex formed and the velocity at the openings increased.

The Flow Characteristics
The simulation calculation settings in this section are the same as those in Section 2. The pumping flow characteristics are the same as well.However, the simulation results showed that the indoor airflow oscillation frequency was smaller than that of building 1 (0.0224 Hz); their values were 0.0073 Hz, 0.0075 Hz, 0.0075 Hz, and 0.0079 Hz, respectively.
The exterior flow fields of the single-room and multi-room building are depicted in Figure 21.It was possible to see that with the same incoming wind speed, the building body's expansion caused a gradient increase in the velocity.In a result, a larger vortex formed and the velocity at the openings increased.The indoor flow fields of the four rooms at t = 900 s are examined in Figure 22.It was discovered that for the side rooms (rooms 1 and 2), there was a greater difference in velocity magnitude at the openings, and the airflow entered the room at an angle to the windows; for the middle rooms 3 and 4, there was a smaller difference in velocity The indoor flow fields of the four rooms at t = 900 s are examined in Figure 22.It was discovered that for the side rooms (rooms 1 and 2), there was a greater difference in velocity magnitude at the openings, and the airflow entered the room at an angle to the windows; for the middle rooms 3 and 4, there was a smaller difference in velocity magnitude at the openings and the airflow flowed in and out of the windows almost parallel to the openings.

Ventilation Flow Rates
Table 9 lists the ventilation flow rates in the various rooms.According to the table, room 4 was the sole instance in which the dimensionless fluctuation ventilation flow rates are larger than the dimensionless mean value with a 20% difference.While the dimensionless fluctuating ventilation flow rates in rooms 1 and 2 were much less than the average value with differences of 36% and 27%, respectively, and the dimensionless average ventilation flow rates in room 3 were comparable to the fluctuating value with a difference of roughly 7%.The greatest value of the dimensionless fluctuating ventilation rate (room 4) was 1.3 times higher than the minimum value (room 1) in these four rooms, and the maximum value of the dimensionless mean ventilation rate (room 4) was 39% different from the lowest value (room 1).

Ventilation Flow Rates
Table 9 lists the ventilation flow rates in the various rooms.According to the table, room 4 was the sole instance in which the dimensionless fluctuation ventilation flow rates are larger than the dimensionless mean value with a 20% difference.While the dimensionless fluctuating ventilation flow rates in rooms 1 and 2 were much less than the average value with differences of 36% and 27%, respectively, and the dimensionless average ventilation flow rates in room 3 were comparable to the fluctuating value with a difference of roughly 7%.The greatest value of the dimensionless fluctuating ventilation rate (room 4) was 1.3 times higher than the minimum value (room 1) in these four rooms, and the maximum value of the dimensionless mean ventilation rate (room 4) was 39% different from the lowest value (room 1).
In summary, room 1 had the lowest ventilation flow rate, while room 4 had the greatest.The reason is the vortex produced by the side walls of the building.Moreover, it can be seen from Table 9 that the ventilation flow rates of single-room buildings are smaller than those of multi-room buildings under the condition of same incoming wind speed and same building height.The reason may still be the flow split caused by the sidewall.

Figure 1 .
Figure 1.The flow chart of the study.Figure 1.The flow chart of the study.

Figure 1 .
Figure 1.The flow chart of the study.Figure 1.The flow chart of the study.

2 . 1 .
Model Simplification and Simulation Settings 2.1.1.Computational Geometry and Domain A three-story building was taken as the model 1 (e.g., Figure 2).The total height of the building H is 9 m, and there is one room on each floor.The dimension of each room is 6 m × 6 m × 3 m (B × B × 1/3 H).The natural ventilation room is located on the second floor with two leeward openings of the same size: 1 m × 1 m.The centers of the two openings are 3 m apart and 1.5 m from the respective closer sides of the building.The three rooms in the building are not connected to each other and are separated by floors or ceilings.Buildings 2023, 13, x FOR PEER REVIEW 3 of 21

2. 1 . 1 .
Computational Geometry and DomainA three-story building was taken as the model 1 (e.g., Figure2).The total height of the building H is 9 m, and there is one room on each floor.The dimension of each room is 6 m × 6 m × 3 m (B × B × 1/3 H).The natural ventilation room is located on the second floor with two leeward openings of the same size: 1 m × 1 m.The centers of the two openings are 3 m apart and 1.5 m from the respective closer sides of the building.The three rooms in the building are not connected to each other and are separated by floors or ceilings.

Figure 3 .
Figure 3. Schematic diagram of different opening distances.

Figure 4 .
Figure 4. Schematic diagram of grid division.Figure 4. Schematic diagram of grid division.

Figure 4 .
Figure 4. Schematic diagram of grid division.Figure 4. Schematic diagram of grid division.

Figure 6 .
Figure 6.Airflow velocity distribution in the plane with Z = 4.5 m.(a) Y = 3 m.(b) Y = 9 m.2.1.3.Validation of CFD Methodology In 2003, Jiang et al. [20] conducted a wind tunnel experiment for natural ventilation.The study focused on a small (0.25 cubic) building model with a window opening size 0.084 m × 0.125 m (length × height), as shown in Figure7.There is a workspace with the area of 2.0 m × 2.0 m, and the height 1.0 m, which is in the wind tunnel.The velocity measurements were conducted on a grid plane at the center of the building.In order to verify the CFD method used in this paper, Jiang's case was chosen for the numerical simulation.The downstream length was set as 8 H; the upstream length was 4 H; and the lateral length was 4 H on both sides of the building in the computational zone.The total number of the grids was 969,570.The boundary conditions and the wind speed were consistent with the literature, and the turbulence model was the SST k-ω.The numerical results are shown to have good agreement with the experimental data, and the relative error is less than 20% in the most regions (Figure8).Therefore, the CFD method adopted in the paper was validated for the numerical simulation of the natural ventilation.

21 Figure 7 .
Figure 7. Schematic diagram of wind tunnel experiment room.Figure 7. Schematic diagram of wind tunnel experiment room.

Figure 7 .
Figure 7. Schematic diagram of wind tunnel experiment room.Figure 7. Schematic diagram of wind tunnel experiment room.

Figure 8 .
Figure 8. Validation of CFD methodology.(a) Schematic diagram of the location of wind tunnel experimental data measurements.(b) Numerical simulation of the velocity profile along the measurement line (red line) compared to the experimental values (black solid dots).

Figure 10 .
Figure 10.Velocity distribution on the vertical surface of the opening.(a) Transient calculation (b) Steady-state calculations.

Figure 10 .
Figure 10.Velocity distribution on the vertical surface of the opening.(a) Transient calculations.(b) Steady-state calculations.
displays the streamline distribution of the flow field with the cross-section situated at Z = 4.5 m.Buildings 2023, 13, x FOR PEER REVIEW 9 of 21

Figure 14 .Figure 14 .
Figure 14.Schematic of the location of the different monitoring points.

Figure 14 .Figure 15 .
Figure 14.Schematic of the location of the different monitoring points.

Figure 17 .
Figure 17.Time-domain plot of ventilation for two openings.

Figure 18
Figure 18 depicts the oscillation frequency of the airflow with various wind speeds showing there was a significant association between the oscillation frequency of the

Figure 16 .
Figure 16.Linear relationship between Ro number and Re number.

Figure 17 .
Figure 17.Time-domain plot of ventilation for two openings.

Figure 18
Figure 18 depicts the oscillation frequency of the airflow with various wind speeds, showing there was a significant association between the oscillation frequency of the

Figure 17 .
Figure 17.Time-domain plot of ventilation for two openings.

Figure 18 .
Figure 18.Frequency of airflow oscillations in different locations with different wind speeds.

Figure 18 .
Figure 18.Frequency of airflow oscillations in different locations with different wind speeds.

Figure 19 .
Figure 19.The relationship of ventilation flow rates and wind speed.

Figure 19 .
Figure 19.The relationship of ventilation flow rates and wind speed.

Figure 20 .
Figure 20.Schematic diagram and computational domain of the multi-room building.3.2.1.The Flow CharacteristicsThe simulation calculation settings in this section are the same as those in Section 2. The pumping flow characteristics are the same as well.However, the simulation results

Figure 20 .
Figure 20.Schematic diagram and computational domain of the multi-room building.

Figure 20 .
Figure 20.Schematic diagram and computational domain of the multi-room building.

Table 1 .
Time steps at different wind speeds.

Table 1 .
Time steps at different wind speeds.

Table 2 .
Results for different cases.

Table 2 .
Results for different cases.

Table 3 .
Frequency of vortex shedding at 4 points under different wind speeds.

Table 4 .
Parameter values with different wind speeds.Therefore, it could be judged that there was also a good linear relationship between the Ro number and Re number in this study case, which was similar to that in other research.
Figure 16.Linear relationship between Ro number and Re number.

Table 5 .
Dimensionless mean and fluctuating ventilation flow rates with different wind speeds.

Table 5 .
Dimensionless mean and fluctuating ventilation flow rates with different wind speeds.

Table 6 .
Ventilation flow rates obtained from different calculation models.

Table 7 .
Ventilation flow rates obtained by fitted equations.

Table 8 .
Dimensionless ventilation flow rates with different opening distances.