Numerical Research on the Mixture Mechanism of Polluted and Fresh Air at the Staggered Tunnel Portals

: In longitudinal ventilation, circulating air is formed in portals for closely spaced twin tunnels, which causes mixing between the polluted air exhausted from one tunnel and the fresh air ﬂow of another tunnel, thus leading to the rising costs of ventilation system construction and operation. In this study, for the closely spaced tunnel with staggered inlet and outlet, the computational ﬂuid dynamics (CFD) numerical simulation method was adopted to reveal ﬂow characteristics of the circulating air as well as variation rules of the circulating air mixing ratio ϕ c with tunnel structure and operation parameters. Results show that both reducing inlet air velocity and increasing outlet air velocity and lateral distance can reduce the impact of the negative-pressure zone at the tunnel entrance on the jet ﬂow structure at the tunnel exit, thus weakening the circulating air. When the inlet is placed behind or aligned with the outlet (staggered distance ∆ l ≤ 0), ϕ c will increase linearly along with the increase of staggered distance; when the inlet is placed before the outlet ( ∆ l > 0), ϕ c will ﬁrst increase and then decrease with the increase of staggered distance. An expression to predict circulating air mixing ratio was created by sections. The predictions show a good correlation with the measurements and indicate that the front slope gradient of the tunnel portal is also one of the factors affecting the circulating air mixing ratio.


Introduction
With the current state of economic development and increasing transportation demands, more and more tunnels are being put into construction [1]. Tunnels can be referred to as a special type of semi-enclosed passageway. When vehicles pass through the tunnels, the exhaust gas will accumulate and therefore increase the concentration of harmful substances such as carbon monoxide, nitrogen oxides, and particulates, which will deteriorate the operating environment in the tunnel [2,3] and even cause a series of health problems such as dyspnea and cardiovascular diseases [4,5]. Usually, ventilation equipment is adopted to dilute the pollutant concentration in the tunnels [6][7][8]. In recent years, due to being restricted by landforms and architectural patterns, more and more twin tunnels with small spacing are being constructed [9,10]. During the operation of this kind of tunnel, not only do the pollutants of the outlet tunnel influence the local environment of the tunnel exit [11], but the pollutants also mix with new incoming air at the adjacent inlet to form circulating air, consequently intensifying the pollution level at the internal environment of the inlet tunnel and resulting in increasing costs of ventilation system construction and operation.
At present, many scholars have studied the circulating air at the portals of twin tunnels with small spacing. Through three-dimensional numerical simulations and field tests, Chen et al. [12] presented how the mixed volume of pollutants correlated with the tunnel spacing and airflow velocity at the inlet/outlet. They concluded that the mixed proportion of the polluted airflow could be reduced by increasing the outlet air velocity and lateral distance of the tunnel portals as well as decreasing the inlet air velocity. Through a model test, Tan et al. [13] established a testing ventilation system of a close tunnel model and conducted a series of experiments. The results showed that when inlet and outlet airflow velocity were equal, the mixed volume of pollutants was barely influenced by the airflow velocity but would decrease with the increase of the tunnel spacing. Based on the experimental data, they also obtained the quantitative relationship between the mixed volume of pollutants and inlet/outlet airflow velocity and tunnel spacing. Through numerical simulation, Yang et al. [14] further analyzed the influence of ambient wind on the pollutant diffusion at the tunnel portals. They concluded that under adverse ambient wind conditions, the degree of influence of the exhaust pollutant at outlet air tunnel portal on the inlet air tunnel lessened as the environmental air slowed down. According to previous findings [12,15], the circulating air can be effectively weakened by establishing midboards, green belts, or building horizontal ventilation tunnels.
It should be noted that many closely spaced twin tunnels with staggered inlet and outlet exist due to landform restrictions. For instance, for the Baziling Tunnel in the Hurong highway, the inlet is about 5 m behind the outlet. At the north opening of the #1 tunnel of the Zizhi Tunnel in Hangzhou, the inlet is about 19 m before the outlet; at the south opening of the #3 tunnel, the inlet is about 44 m before the outlet. The circulating air can be effectively influenced by staggering the inlet and outlet air tunnel portals. For example, by increasing the length of outlet air tunnel portal and repositioning it at a point 5 m away from the inlet, the proportion of polluted airflow mixed with fresh airflow will be lower than 5% [13]. However, previous research mostly focused on closely spaced twin tunnels with aligned inlet and outlet air tunnel portals, and studies relating to the circulating air of the staggered inlet and outlet are still inadequate.
The present research focuses on closely spaced twin tunnels with staggered inlet and outlet and uses the computational fluid dynamics (CFD) simulation technique to analyze the characteristics of circulating air, as well as to acquire the variation of the intake volume of pollutants caused by the circulating wind with the airflow velocity, tunnel spacing, and staggered distance. Based on these results, an expression for predicting the mixed proportion of the polluted airflow is suggested, which can provide an important reference for the economic and environmental protection evaluation of closely spaced twin tunnels design.

Governing Equations and Numerical Algorithm
The airflow and pollutant diffusion at the tunnel portals are influenced by various factors. In order to explore the inherent law of the flow field, the following assumptions are put forward in this paper: (1) vehicles are regarded as external and accidental factors, with no influence on the airflow; (2) natural wind is regarded as the external boundary at a velocity of 0 m·s −1 ; (3) the airflow temperature at all inlet and outlet boundaries is 300 K, without any volume source; (4) as the pollutant components are scalar with similar distributional patterns, CO concentration distribution is used to represent the pollutant diffusion.
The airflow velocity at the portals is much slower than the sound velocity, so it can be regarded as a steady turbulent flow of incompressible fluid [14]. This paper calculated the airflow at the inlet and outlet based on the Standard k-ε turbulence model. The pollutant diffusion was described through the species transport equation. The buoyancy-driven effects caused by density changes due to varied components of mixed gas were not considered. The transport equations of the turbulent kinetic energy k and turbulent dissipation rate ε in the Standard k-ε turbulence model and species transport equation are shown in Equations (1)- (3).
In the above formula, both the values of i and j range from 1 to 3; x i represents the coordinate components; ρ, µ, ν, and u are the mass density of air, turbulent viscosity, coefficient of kinematic viscosity, and air velocity; c s refers to the volume density of s; D s refers to the diffusion coefficient of s; S s refers to the quality of s generated from the chemical reaction in unit time and volume, namely, the productivity; C µ , σ k , σ , C 1 , C 2 , and C 3 are model coefficients that follow the principle of asymptotic proportionality [16], and their values are 0.09, 0.24, 0.15, 1.44, 1.92, and 0.2.
Equations (1)-(3) and the Reynolds-averaged Navier-Stokes equation [17,18] constituted the governing equations for the pollutant diffusion of closely spaced twin tunnels. The finite volume method [19] was used to discretize the governing equations. For the discretization of derivatives, the skewness-corrected linear upwind interpolation scheme (second-order accuracy) was applied to the convective terms, and second-order linear interpolation was used with a nonorthogonal correction scheme for the diffusive terms [20]. The solution of the velocity field and pressure field was conducted through a fully coupled block algorithm. Due to the simultaneous solution of momentum and continuity equations, implicit block coupling of pressure and velocity variables leads to faster convergence with respect to classical, loosely coupled, segregated algorithms of the SIMPLE family of algorithms [21].

Computational Geometry, Domain, and Grid
The airflow within about a 450-m length, 200-m width, and 50-m height range outside the tunnel portals was adopted to simulate the pollutant diffusion. The tunnel section resembled the shape of a horseshoe. The two-lane and three-lane tunnels were considered separately, and the hydraulic diameters D were 6.1 and 8.6 m, respectively. In reference to a certain tunnel project, the front slope gradient of the tunnel portal was 15 • . The lateral distance of inlet and outlet d was set within the range of 5-20 m at 2.5-m intervals. The staggered distance ∆l was from −10 to 320 m at 5 m intervals. If ∆l < 0, the inlet would be positioned behind the outlet as shown in Figure 1a; if ∆l = 0, the outlet and inlet would be aligned as shown in Figure 1b; if ∆l > 0, the inlet would be positioned in front of the outlet as shown in Figure 1c.
In the above formula, both the values of i and j range from 1 to 3; xi represents the coordinate components; ρ, μ, ν, and u are the mass density of air, turbulent viscosity, coefficient of kinematic viscosity, and air velocity; cs refers to the volume density of s; Ds refers to the diffusion coefficient of s; Ss refers to the quality of s generated from the chemical reaction in unit time and volume, namely, the productivity; Cμ, σk, σ , C1 , C2 , and C3 are model coefficients that follow the principle of asymptotic proportionality [16], and their values are 0.09, 0.24, 0.15, 1.44, 1.92, and 0.2.
Equations (1)-(3) and the Reynolds-averaged Navier-Stokes equation [17,18] constituted the governing equations for the pollutant diffusion of closely spaced twin tunnels. The finite volume method [19] was used to discretize the governing equations. For the discretization of derivatives, the skewness-corrected linear upwind interpolation scheme (second-order accuracy) was applied to the convective terms, and second-order linear interpolation was used with a nonorthogonal correction scheme for the diffusive terms [20]. The solution of the velocity field and pressure field was conducted through a fully coupled block algorithm. Due to the simultaneous solution of momentum and continuity equations, implicit block coupling of pressure and velocity variables leads to faster convergence with respect to classical, loosely coupled, segregated algorithms of the SIMPLE family of algorithms [21].

Computational Geometry, Domain, and Grid
The airflow within about a 450-m length, 200-m width, and 50-m height range outside the tunnel portals was adopted to simulate the pollutant diffusion. The tunnel section resembled the shape of a horseshoe. The two-lane and three-lane tunnels were considered separately, and the hydraulic diameters D were 6.1 and 8.6 m, respectively. In reference to a certain tunnel project, the front slope gradient of the tunnel portal was 15°. The lateral distance of inlet and outlet d was set within the range of 5-20 m at 2.5-m intervals. The staggered distance Δl was from −10 to 320 m at 5 m intervals. If Δl < 0, the inlet would be positioned behind the outlet as shown in Figure 1a; if Δl = 0, the outlet and inlet would be aligned as shown in Figure 1b; if Δl > 0, the inlet would be positioned in front of the outlet as shown in Figure 1c. The whole domain was meshed by nonuniform hexahedral grids using the professional mesh generation software ICEM-CFD17.2 (ANSYS Inc. Pittsburgh, PA, USA, 2016), and the areas near the inlet and outlet were locally refined to ensure calculation efficiency, accuracy, and reliability. Figure  2 shows the grid map of a three-lane tunnel with a lateral distance d = 10 m. The whole domain was meshed by nonuniform hexahedral grids using the professional mesh generation software ICEM-CFD17.2 (ANSYS Inc. Pittsburgh, PA, USA, 2016), and the areas near the inlet and outlet were locally refined to ensure calculation efficiency, accuracy, and reliability. Figure 2 shows the grid map of a three-lane tunnel with a lateral distance d = 10 m.

Grid Dependence Study
Grid density exerts a great impact on the numerical simulation results. When the grid number is relatively small, the dispersion error will be very large; when the grid number is large, not only will the rounding error increase, but the calculation quantity will increase as well. Figure 3 shows that the velocity profile in the symmetry plane is 30 m away from the outlet under different grid orders of 0.08 × 10 6 , 0.47 × 10 6 , 1.05 × 10 6 , and 5.78 × 10 6 . It can be seen that the number of grids has a great effect on the velocity profile. When the number of grids reaches 1.05 × 10 6 , the velocity profile will no longer change with the increase of grid number.

Boundary Conditions
A Neumann boundary condition for pressure was used at the outer boundary of the computational domain. The tunnel walls and pavement were taken as an impenetrable, nonslip boundary. It should be noted that the Standard k-ε model as a turbulence model for high Re numbers is not applicable to the near-wall area where the flowing Re number is low. Therefore, in order to capture the near-wall region's flow structures, the kqRwallFunction and epsilonWallFunction, also called standard wall functions in OpenFOAM 2.1.1 (OpenCFD Ltd., Bracknell, UK, 2012) [20], were assigned to the turbulence kinetic energy and dissipation rate, respectively. Using the wall function allowed us to locate the first point adjacent to the wall in the logarithmic law instead of the sublayer. The calculated y+ value should be distributed between 30 and 300 in all simulations. According to y+, the distance of the near-wall node to the solid surface was calculated, which was also corrected continuously through computational simulation. Ultimately, it could be determined that y+ can be ensured between 32 and 132 under all computation conditions when the first layer grid thickness was 0.018 m.
The tunnel inlet and outlet adopted the Dirichlet boundary condition. According to the ventilation conditions of the actual tunnel and the possible working conditions during the project's operational phase, calculations were separately conducted with respect to the following working conditions: at outlet air velocity uout of 1.

Grid Dependence Study
Grid density exerts a great impact on the numerical simulation results. When the grid number is relatively small, the dispersion error will be very large; when the grid number is large, not only will the rounding error increase, but the calculation quantity will increase as well. Figure 3 shows that the velocity profile in the symmetry plane is 30 m away from the outlet under different grid orders of 0.08 × 10 6 , 0.47 × 10 6 , 1.05 × 10 6 , and 5.78 × 10 6 . It can be seen that the number of grids has a great effect on the velocity profile. When the number of grids reaches 1.05 × 10 6 , the velocity profile will no longer change with the increase of grid number.

Grid Dependence Study
Grid density exerts a great impact on the numerical simulation results. When the grid number is relatively small, the dispersion error will be very large; when the grid number is large, not only will the rounding error increase, but the calculation quantity will increase as well. Figure 3 shows that the velocity profile in the symmetry plane is 30 m away from the outlet under different grid orders of 0.08 × 10 6 , 0.47 × 10 6 , 1.05 × 10 6 , and 5.78 × 10 6 . It can be seen that the number of grids has a great effect on the velocity profile. When the number of grids reaches 1.05 × 10 6 , the velocity profile will no longer change with the increase of grid number.

Boundary Conditions
A Neumann boundary condition for pressure was used at the outer boundary of the computational domain. The tunnel walls and pavement were taken as an impenetrable, nonslip boundary. It should be noted that the Standard k-ε model as a turbulence model for high Re numbers is not applicable to the near-wall area where the flowing Re number is low. Therefore, in order to capture the near-wall region's flow structures, the kqRwallFunction and epsilonWallFunction, also called standard wall functions in OpenFOAM 2.1.1 (OpenCFD Ltd., Bracknell, UK, 2012) [20], were assigned to the turbulence kinetic energy and dissipation rate, respectively. Using the wall function allowed us to locate the first point adjacent to the wall in the logarithmic law instead of the sublayer. The calculated y+ value should be distributed between 30 and 300 in all simulations. According to y+, the distance of the near-wall node to the solid surface was calculated, which was also corrected continuously through computational simulation. Ultimately, it could be determined that y+ can be ensured between 32 and 132 under all computation conditions when the first layer grid thickness was 0.018 m.
The tunnel inlet and outlet adopted the Dirichlet boundary condition. According to the ventilation conditions of the actual tunnel and the possible working conditions during the project's operational phase, calculations were separately conducted with respect to the following working conditions: at outlet air velocity uout of 1.

Boundary Conditions
A Neumann boundary condition for pressure was used at the outer boundary of the computational domain. The tunnel walls and pavement were taken as an impenetrable, nonslip boundary. It should be noted that the Standard k-ε model as a turbulence model for high Re numbers is not applicable to the near-wall area where the flowing Re number is low. Therefore, in order to capture the near-wall region's flow structures, the kqRwallFunction and epsilonWallFunction, also called standard wall functions in OpenFOAM 2.1.1 (OpenCFD Ltd., Bracknell, UK, 2012) [20], were assigned to the turbulence kinetic energy and dissipation rate, respectively. Using the wall function allowed us to locate the first point adjacent to the wall in the logarithmic law instead of the sublayer. The calculated y+ value should be distributed between 30 and 300 in all simulations. According to y+, the distance of the near-wall node to the solid surface was calculated, which was also corrected continuously through computational simulation. Ultimately, it could be determined that y+ can be ensured between 32 and 132 under all computation conditions when the first layer grid thickness was 0.018 m.
The tunnel inlet and outlet adopted the Dirichlet boundary condition. According to the ventilation conditions of the actual tunnel and the possible working conditions during the project's operational phase, calculations were separately conducted with respect to the following working conditions: at outlet air velocity u out of 1.5 m·s −1 , 3.0 m·s −1 , 4.5 m·s −1 , 6.0 m·s −1 , and 7.5 m·s −1 ; and inlet air velocity u in of 1.5 m·s −1 , 3.0 m·s −1 , 4.5 m·s −1 , 6.0 m·s −1 , and 7.5 m·s −1 . The numerical computations in this study were conducted by OpenFOAM, an opensource CFD code with flexible and extendable C++ libraries.

Validation
The discharge of the tunnel outlet pollutants was horizontal and was constrained by the ground boundaries, all of which constituted a typical three-dimensional wall jet problem. Rajaratnam [22] conducted a series of wall jet experiments in a semi-infinite three-dimensional space. This paper attempted to verify the numerical model of the airflow by utilizing the experimental findings of the jet flow over the circular nozzle. The diameter D of the circular nozzle was 0.00953 m and the jet flow outlet speed u 0 was 7.42 m·s −1 . The comparison between the numerical results and the measured values is shown in Figures 4 and 5, which respectively show the decay of nondimensional maximum velocity u max in the symmetry plane of the three-dimensional wall jet, and the dimensionless velocity profiles at the plane of symmetry at distances x = 0.165, 0.267, and 0.470 m from the jet flow outlet. In Figure 5, u y represents the speed at the height y above the ground, and b y is the length scale of the plane of symmetry [23]. It shows that CFD simulation results are consistent with the trend of the test data, with a maximum error of only 5.2%.

Validation
The discharge of the tunnel outlet pollutants was horizontal and was constrained by the ground boundaries, all of which constituted a typical three-dimensional wall jet problem. Rajaratnam [22] conducted a series of wall jet experiments in a semi-infinite three-dimensional space. This paper attempted to verify the numerical model of the airflow by utilizing the experimental findings of the jet flow over the circular nozzle. The diameter D of the circular nozzle was 0.00953 m and the jet flow outlet speed u0 was 7.42 m•s −1 . The comparison between the numerical results and the measured values is shown in Figures 4 and 5, which respectively show the decay of nondimensional maximum velocity umax in the symmetry plane of the three-dimensional wall jet, and the dimensionless velocity profiles at the plane of symmetry at distances x = 0.165, 0.267, and 0.470 m from the jet flow outlet. In Figure 5, uy represents the speed at the height y above the ground, and by is the length scale of the plane of symmetry [23]. It shows that CFD simulation results are consistent with the trend of the test data, with a maximum error of only 5.2%.  In addition to the above, this paper also used the scale model test results of the circulating air in closely spaced twin tunnels in the existing literature [13] to verify the species transport model. The scale model was 10-m long, and the section was a semicircle with a radius of 0.5 m. The lateral distance between the portals d = 4 m, and the front slope of the tunnel portal was perpendicular to the ground. The outlet air velocity uin was 3 m•s −1 and the air velocity uout was 1 m•s −1 , 3 m•s −1 , 5 m•s −1 , 7 m•s −1 , and 9 m•s −1 , respectively. The comparison between the numerical simulation and the scale test results is shown in Figure 6. Results show that change patterns are roughly consistent, with a maximum error of only 5.5%, which indicates that the mathematical model adopted by this paper

Validation
The discharge of the tunnel outlet pollutants was horizontal and was constrained by the ground boundaries, all of which constituted a typical three-dimensional wall jet problem. Rajaratnam [22] conducted a series of wall jet experiments in a semi-infinite three-dimensional space. This paper attempted to verify the numerical model of the airflow by utilizing the experimental findings of the jet flow over the circular nozzle. The diameter D of the circular nozzle was 0.00953 m and the jet flow outlet speed u0 was 7.42 m•s −1 . The comparison between the numerical results and the measured values is shown in Figures 4 and 5, which respectively show the decay of nondimensional maximum velocity umax in the symmetry plane of the three-dimensional wall jet, and the dimensionless velocity profiles at the plane of symmetry at distances x = 0.165, 0.267, and 0.470 m from the jet flow outlet. In Figure 5, uy represents the speed at the height y above the ground, and by is the length scale of the plane of symmetry [23]. It shows that CFD simulation results are consistent with the trend of the test data, with a maximum error of only 5.2%.  In addition to the above, this paper also used the scale model test results of the circulating air in closely spaced twin tunnels in the existing literature [13] to verify the species transport model. The scale model was 10-m long, and the section was a semicircle with a radius of 0.5 m. The lateral distance between the portals d = 4 m, and the front slope of the tunnel portal was perpendicular to the ground. The outlet air velocity uin was 3 m•s −1 and the air velocity uout was 1 m•s −1 , 3 m•s −1 , 5 m•s −1 , 7 m•s −1 , and 9 m•s −1 , respectively. The comparison between the numerical simulation and the scale test results is shown in Figure 6. Results show that change patterns are roughly consistent, with a maximum error of only 5.5%, which indicates that the mathematical model adopted by this paper In addition to the above, this paper also used the scale model test results of the circulating air in closely spaced twin tunnels in the existing literature [13] to verify the species transport model. The scale model was 10-m long, and the section was a semicircle with a radius of 0.5 m. The lateral distance between the portals d = 4 m, and the front slope of the tunnel portal was perpendicular to the ground.
The outlet air velocity u in was 3 m·s −1 and the air velocity u out was 1 m·s −1 , 3 m·s −1 , 5 m·s −1 , 7 m·s −1 , and 9 m·s −1 , respectively. The comparison between the numerical simulation and the scale test results is shown in Figure 6. Results show that change patterns are roughly consistent, with a maximum error of only 5.5%, which indicates that the mathematical model adopted by this paper demonstrates a high degree of accuracy and reliability and is suitable for calculating and simulating the circulating air at tunnel portals. demonstrates a high degree of accuracy and reliability and is suitable for calculating and simulating the circulating air at tunnel portals.  Figure 7 depicts the sketches of three-dimensional velocity profiles of jet flow at a three-lane tunnel exit. It can be observed that after the air is pushed out from the tunnel exit, due to the lateral fluctuation on the turbulence, the moving wind exchanges momentum with its surrounding stagnant air and drives the surrounding air molecules to move to produce an entrainment effect, thus continuously expanding the jet flow cross-section area. The three-dimensional wall jet can be categorized into the potential core region and decay region, depending on whether the maximum velocity umax decreases. The decay of nondimensional umax is shown by the solid line in Figure 8. It can be seen from Figures 7 and 8 that at the potential core region (where x ≤ 60 m), the jet flow fiercely entrains the surrounding air, thereby expanding the jet flow boundary and the flow rate altogether, while reducing the air velocity at the same time, but the umax remains the same. At the decay region (where x > 60 m), with further development of the jet-flow entrainment and mixing effects, the umax starts to reduce and the velocity profile also tends to flatten out.   Figure 7 depicts the sketches of three-dimensional velocity profiles of jet flow at a three-lane tunnel exit. It can be observed that after the air is pushed out from the tunnel exit, due to the lateral fluctuation on the turbulence, the moving wind exchanges momentum with its surrounding stagnant air and drives the surrounding air molecules to move to produce an entrainment effect, thus continuously expanding the jet flow cross-section area. The three-dimensional wall jet can be categorized into the potential core region and decay region, depending on whether the maximum velocity u max decreases. The decay of nondimensional u max is shown by the solid line in Figure 8. It can be seen from Figures 7 and 8 that at the potential core region (where x ≤ 60 m), the jet flow fiercely entrains the surrounding air, thereby expanding the jet flow boundary and the flow rate altogether, while reducing the air velocity at the same time, but the u max remains the same. At the decay region (where x > 60 m), with further development of the jet-flow entrainment and mixing effects, the u max starts to reduce and the velocity profile also tends to flatten out.

Flow Characteristic
demonstrates a high degree of accuracy and reliability and is suitable for calculating and simulating the circulating air at tunnel portals.  Figure 7 depicts the sketches of three-dimensional velocity profiles of jet flow at a three-lane tunnel exit. It can be observed that after the air is pushed out from the tunnel exit, due to the lateral fluctuation on the turbulence, the moving wind exchanges momentum with its surrounding stagnant air and drives the surrounding air molecules to move to produce an entrainment effect, thus continuously expanding the jet flow cross-section area. The three-dimensional wall jet can be categorized into the potential core region and decay region, depending on whether the maximum velocity umax decreases. The decay of nondimensional umax is shown by the solid line in Figure 8. It can be seen from Figures 7 and 8 that at the potential core region (where x ≤ 60 m), the jet flow fiercely entrains the surrounding air, thereby expanding the jet flow boundary and the flow rate altogether, while reducing the air velocity at the same time, but the umax remains the same. At the decay region (where x > 60 m), with further development of the jet-flow entrainment and mixing effects, the umax starts to reduce and the velocity profile also tends to flatten out.   When ∆l ≤ 0, the jet flow at the outlet deflects entirely under the influence of the negative-pressure zone at the inlet, as shown in Figure 9a-c. In the meantime, by referring to Figure 8a, it can be noticed that, different from the standard three-dimensional wall jet, the umax at the moment begins to attenuate when x > 0 (the bigger ∆l is, the quicker the attenuation will be); when ∆l > 0, we can learn from Figures 9d,f and 8b that the jet flow structure at the outlet when x ≤ ∆l and the law of umax attenuation are both similar to those of the standard three-dimensional wall jet. However, when x > ∆l, the jet flow will deflect under the influence of the negative-pressure zone at the inlet, for which umax will fall rapidly at a decay rate first increasing and then decreasing along with the increase of ∆l.

Flow Characteristic
Moreover, it can be perceived from the comparison between Figure 9g-h that the bigger uout is, the bigger the inertia effect of the jet flow will be, and the less vulnerable it will be to the influence of the negative-pressure zone at the inlet and the fuller the jet flow will develop. In this way, the polluted air flowing into the adjacent air-inlet tunnel through the air-exhausting tunnel will be reduced. According to Figure 9g-i, when uin increases, the bigger the influencing area of the negative-pressure zone is, the bigger the damage caused to the flow structure at the outlet will be, thus resulting in more polluted air flowing into the adjacent air-inlet tunnel through the air-exhausting tunnel. Finally, through comparing Figure 9i with Figure 9j, the bigger d is, the smaller the damage caused by the negative-pressure zone at the inlet to the flow structure at the outlet will be and the fuller the jet flow will develop; when d = 20 m, there is almost no polluted air flowing into the adjacent air-inlet tunnel through the air-exhausting tunnel.
Therefore, it can be concluded that, for closely spaced twin tunnels with staggered inlet and outlet, the air flow discharged from their air-exhausting tunnels will diffuse transversely toward the air-inlet tunnels under the influence of the negative-pressure zone at the inlet, even causing inverse flow and forming the circulating air. Parameters including inlet air velocity uin, outlet air velocity uout, lateral distance of the tunnel portals d, and staggered distance ∆l are all important factors influencing the circulating air.  When ∆l ≤ 0, the jet flow at the outlet deflects entirely under the influence of the negative-pressure zone at the inlet, as shown in Figure 9a-c. In the meantime, by referring to Figure 8a, it can be noticed that, different from the standard three-dimensional wall jet, the u max at the moment begins to attenuate when x > 0 (the bigger ∆l is, the quicker the attenuation will be); when ∆l > 0, we can learn from Figure 9d,f and Figure 8b that the jet flow structure at the outlet when x ≤ ∆l and the law of u max attenuation are both similar to those of the standard three-dimensional wall jet. However, when x > ∆l, the jet flow will deflect under the influence of the negative-pressure zone at the inlet, for which u max will fall rapidly at a decay rate first increasing and then decreasing along with the increase of ∆l.
Moreover, it can be perceived from the comparison between Figure 9g-h that the bigger u out is, the bigger the inertia effect of the jet flow will be, and the less vulnerable it will be to the influence of the negative-pressure zone at the inlet and the fuller the jet flow will develop. In this way, the polluted air flowing into the adjacent air-inlet tunnel through the air-exhausting tunnel will be reduced. According to Figure 9g-i, when u in increases, the bigger the influencing area of the negative-pressure zone is, the bigger the damage caused to the flow structure at the outlet will be, thus resulting in more polluted air flowing into the adjacent air-inlet tunnel through the air-exhausting tunnel. Finally, through comparing Figure 9i with Figure 9j, the bigger d is, the smaller the damage caused by the negative-pressure zone at the inlet to the flow structure at the outlet will be and the fuller the jet flow will develop; when d = 20 m, there is almost no polluted air flowing into the adjacent air-inlet tunnel through the air-exhausting tunnel.
Therefore, it can be concluded that, for closely spaced twin tunnels with staggered inlet and outlet, the air flow discharged from their air-exhausting tunnels will diffuse transversely toward the air-inlet tunnels under the influence of the negative-pressure zone at the inlet, even causing inverse flow and forming the circulating air. Parameters including inlet air velocity u in , outlet air velocity u out , lateral distance of the tunnel portals d, and staggered distance ∆l are all important factors influencing the circulating air.

Variation Characteristics of Circulating Air Mixing Ratio
The circulating air mixing ratio is represented by φc = mc/mout, where mc is the amount of pollutants that is carried by the circulating air into the inlet tunnel from the outlet tunnel, and mout is the total amount of gaseous pollutants discharged from the outlet tunnel. Therefore, it can be perceived that φc indicates the degree of impacts of exhaust emission from the outlet tunnel on the inlet tunnel.

The Relationship between Circulating Air Mixing Ratio and the Air Velocity at the Inlet/Outlet
The flow field structure depicted in Figure 9g-i shows that both reducing uin and increasing uout can weaken the circulating air. Existing literature [13] has also shown that when uin = uout, for closely spaced twin tunnels with aligned inlet and outlet air tunnel portals, φc does not change with the change of air velocity value but is related to the ratio of uin and uout. Figure 10 shows the variation of φc with air velocity at the closely spaced twin tunnels with staggered inlet and outlet. It can be seen that regardless of the relative position of the inlet and outlet, when uin = uout, φc is not affected by the value of air velocity, and that the variation of φc does not exceed 0.3%.

Variation Characteristics of Circulating Air Mixing Ratio
The circulating air mixing ratio is represented by ϕ c = m c /m out , where m c is the amount of pollutants that is carried by the circulating air into the inlet tunnel from the outlet tunnel, and m out is the total amount of gaseous pollutants discharged from the outlet tunnel. Therefore, it can be perceived that ϕ c indicates the degree of impacts of exhaust emission from the outlet tunnel on the inlet tunnel.

The Relationship between Circulating Air Mixing Ratio and the Air Velocity at the Inlet/Outlet
The flow field structure depicted in Figure 9g-i shows that both reducing u in and increasing u out can weaken the circulating air. Existing literature [13] has also shown that when u in = u out , for closely spaced twin tunnels with aligned inlet and outlet air tunnel portals, ϕ c does not change with the change of air velocity value but is related to the ratio of u in and u out . Figure 10 shows the variation of ϕ c with air velocity at the closely spaced twin tunnels with staggered inlet and outlet. It can be seen that regardless of the relative position of the inlet and outlet, when u in = u out , ϕ c is not affected by the value of air velocity, and that the variation of ϕ c does not exceed 0.3%.

Variation Characteristics of Circulating Air Mixing Ratio
The circulating air mixing ratio is represented by φc = mc/mout, where mc is the amount of pollutants that is carried by the circulating air into the inlet tunnel from the outlet tunnel, and mout is the total amount of gaseous pollutants discharged from the outlet tunnel. Therefore, it can be perceived that φc indicates the degree of impacts of exhaust emission from the outlet tunnel on the inlet tunnel.

The Relationship between Circulating Air Mixing Ratio and the Air Velocity at the Inlet/Outlet
The flow field structure depicted in Figure 9g-i shows that both reducing uin and increasing uout can weaken the circulating air. Existing literature [13] has also shown that when uin = uout, for closely spaced twin tunnels with aligned inlet and outlet air tunnel portals, φc does not change with the change of air velocity value but is related to the ratio of uin and uout. Figure 10 shows the variation of φc with air velocity at the closely spaced twin tunnels with staggered inlet and outlet. It can be seen that regardless of the relative position of the inlet and outlet, when uin = uout, φc is not affected by the value of air velocity, and that the variation of φc does not exceed 0.3%.  Figure 10. Impacts of air velocity u on circulating air mixing ratio φc. Figure 10. Impacts of air velocity u on circulating air mixing ratio ϕ c . Figure 11 shows the variation of ϕ c with u in /u out when u in = u out . It can be seen that ϕ c increases gradually with the increase of u in /u out despite the relative position of the inlet and outlet. By further comparing the numerical calculation results of ϕ c in the tunnel with aligned inlet and outlet air tunnel portals with the scaled model test results in the literature [13], it can be found that the trend of the simulated data is congruent with the test results in the literature, but the values of the simulation are generally lower than the data in the literature. This is caused by the inconsistency between the front slope gradient (15 • ) of the tunnel portal for the numerical simulation in this paper and the front slope gradient (90 • ) of the scaled model in the previous study.
Appl. Sci. 2018, 8, x FOR PEER REVIEW 9 of 15 Figure 11 shows the variation of φc with uin/uout when uin ≠ uout. It can be seen that φc increases gradually with the increase of uin/uout despite the relative position of the inlet and outlet. By further comparing the numerical calculation results of φc in the tunnel with aligned inlet and outlet air tunnel portals with the scaled model test results in the literature [13], it can be found that the trend of the simulated data is congruent with the test results in the literature, but the values of the simulation are generally lower than the data in the literature. This is caused by the inconsistency between the front slope gradient (15°) of the tunnel portal for the numerical simulation in this paper and the front slope gradient (90°) of the scaled model in the previous study.  Figure 11. Impacts of velocity ratio uin/uout on circulating air mixing ratio φc.
In addition, according to Figures 10 and 11, ∆l and d have a great influence on the circulating air; with the decrease of ∆l and the increase of d, φc drops significantly.    In addition, according to Figures 10 and 11, ∆l and d have a great influence on the circulating air; with the decrease of ∆l and the increase of d, ϕ c drops significantly.  Figure 11 shows the variation of φc with uin/uout when uin ≠ uout. It can be seen that φc increases gradually with the increase of uin/uout despite the relative position of the inlet and outlet. By further comparing the numerical calculation results of φc in the tunnel with aligned inlet and outlet air tunnel portals with the scaled model test results in the literature [13], it can be found that the trend of the simulated data is congruent with the test results in the literature, but the values of the simulation are generally lower than the data in the literature. This is caused by the inconsistency between the front slope gradient (15°) of the tunnel portal for the numerical simulation in this paper and the front slope gradient (90°) of the scaled model in the previous study.  Figure 11. Impacts of velocity ratio uin/uout on circulating air mixing ratio φc.

The Relationship between Circulating Air Mixing Ratio and Tunnel Structure Parameters
In addition, according to Figures 10 and 11, ∆l and d have a great influence on the circulating air; with the decrease of ∆l and the increase of d, φc drops significantly.     When ∆l ≤ 0, the bigger the ∆l is, the more intense the impacts of the negative-pressure zone at the tunnel inlet on the jet flow at the tunnel outlet will be and the faster the u max will be attenuated ( Figure 8a). Therefore, ϕ c rises with the increase of ∆l/D. Take the curve d = 5 m in Figure 13b as an example. At ∆l/D = −1.2, the inlet is placed behind the outlet, and combined with Figure 9a, it can be found that the influence of the negative-pressure zone at the inlet on the jet structure at the outlet is very limited, and the jet flow fully develops. At this time, a negligible amount of polluted air permeates into the inlet and ϕ c is only 5.1%. With the increase of ∆l/D, the influence of the negative-pressure zone at the inlet on the jet structure at the outlet is gradually enhanced and ϕ c increases linearly. When ∆l/D = 0, the outlet and the inlet are aligned and the negative-pressure zone at the inlet causes obvious damage to the jet structure at the outlet, but the jet flow still develops to some extent. At this instant, the polluted air partially permeates into the inlet and the ϕ c is 19.3%.

The Relationship between Circulating Air Mixing Ratio and Tunnel Structure Parameters
When ∆l > 0, the decay rate of u max increases firstly with increasing ∆l before witnessing a decrease afterwards (Figure 8b). So ϕ c also rises firstly and then declines with increasing ∆l/D. When d = 5 m, ϕ c will reach the maximum at ∆l/D = 14.6, and the bigger u in /u out is, the greater the variation of ϕ c and ϕ c at unit ∆l/D will be, as shown in Figure 13a. With the gradual increase of d, ∆l/D also rises gradually when ϕ c reaches the maximum. For example, when d = 10 m, ϕ c will be the maximum at ∆l/D = 26.3, as shown in Figure 13b. When ∆l ≤ 0, the bigger the ∆l is, the more intense the impacts of the negative-pressure zone at the tunnel inlet on the jet flow at the tunnel outlet will be and the faster the umax will be attenuated ( Figure 8a). Therefore, φc rises with the increase of ∆l/D. Take the curve d = 5 m in Figure 13b as an example. At ∆l/D = −1.2, the inlet is placed behind the outlet, and combined with Figure 9a, it can be found that the influence of the negative-pressure zone at the inlet on the jet structure at the outlet is very limited, and the jet flow fully develops. At this time, a negligible amount of polluted air permeates into the inlet and φc is only 5.1%. With the increase of ∆l/D, the influence of the negativepressure zone at the inlet on the jet structure at the outlet is gradually enhanced and φc increases linearly. When ∆l/D = 0, the outlet and the inlet are aligned and the negative-pressure zone at the inlet causes obvious damage to the jet structure at the outlet, but the jet flow still develops to some extent. At this instant, the polluted air partially permeates into the inlet and the φc is 19.3%.
When ∆l > 0, the decay rate of umax increases firstly with increasing ∆l before witnessing a decrease afterwards (Figure 8b). So φc also rises firstly and then declines with increasing ∆l/D. When d = 5 m, φc will reach the maximum at ∆l/D = 14.6, and the bigger uin/uout is, the greater the variation of φc and φc at unit ∆l/D will be, as shown in Figure 13a. With the gradual increase of d, ∆l/D also rises gradually when φc reaches the maximum. For example, when d = 10 m, φc will be the maximum at ∆l/D = 26.3, as shown in Figure 13b.

Construction and Validation of Circulating Air Mixing Ratio Model
It can be learned from the analysis above that: (1) φc has no direct relation with the air velocity value but is related to the ratio of air velocity at the inlet and the outlet-uin/uout; (2) with the increase of d, φc decreases nearly exponentially; (3) when ∆l ≤ 0, φc increases linearly with the increase of ∆l/D; when ∆l＞0, however, φc will first increase and then decrease with the increase of ∆l/D; the bigger the d is, the bigger the ∆l/D will be when φc is the maximum. Therefore, the functional relationship between φc and uin/uout, d, ∆l, and D can be described in the following forms:   Figure 13. Impacts of ∆l/D on circulating air mixing ratio ϕ c.

Construction and Validation of Circulating Air Mixing Ratio Model
It can be learned from the analysis above that: (1) ϕ c has no direct relation with the air velocity value but is related to the ratio of air velocity at the inlet and the outlet-u in /u out ; (2) with the increase of d, ϕ c decreases nearly exponentially; (3) when ∆l ≤ 0, ϕ c increases linearly with the increase of ∆l/D; when ∆l>0, however, ϕ c will first increase and then decrease with the increase of ∆l/D; the bigger the d is, the bigger the ∆l/D will be when ϕ c is the maximum. Therefore, the functional relationship between ϕ c and u in /u out , d, ∆l, and D can be described in the following forms: Based on the CFD simulation results, setting ∆l = 0 as the demarcation point, the fitted piecewise expressions of ϕ c can be obtained as follows via the software Matlab with the gradient descent method: A comparison between calculated results of Equation (5) and simulated values of ϕ c is given in Figure 14. The goodness of fit R new is 0.99. A comparison between calculated results of Equation (5) and simulated values of φc is given in Figure 14. The goodness of fit Rnew is 0.99. The circulating air mixing ratio model (Equation (5)) was validated by the measured data of two actual tunnels with small spacing. Both of these tunnels are one-way and double-lane tunnels with a horseshoe-shaped cross section. The structural parameters of the two tunnels are shown in Table 1. The field measurement was conducted on a windless day. Because of the toxicity of CO, this test replaced CO with CO2, which is less harmful, cheaper, and easier to obtain and measure. The pollutant-releasing device of the model test mainly included a CO2 cylinder, a CO2 two-stage pressure-reducing valve (better pressure stabilizing effect), and a vortex flow meter. In the experiment, CO2 gas was released at about 100 m away from the outlet into the outflow tunnel so as to ensure that CO2 was well mixed at the outlet section. The locations of measuring points in #1 tunnel are shown in Figure 15. The jet velocity of the jet fan was regulated in the tunnel in order to change inlet and outlet air velocity. The air velocity of the inlet and outlet tunnel was obtained by a Testo425 thermal anemometer (Testo AG, Lenzkirch, Germany). In fact, the air velocity of a single point on the cross section of the tunnel cannot represent the average air velocity of the entire cross section. Thus, several The circulating air mixing ratio model (Equation (5)) was validated by the measured data of two actual tunnels with small spacing. Both of these tunnels are one-way and double-lane tunnels with a horseshoe-shaped cross section. The structural parameters of the two tunnels are shown in Table 1. The field measurement was conducted on a windless day. Because of the toxicity of CO, this test replaced CO with CO 2 , which is less harmful, cheaper, and easier to obtain and measure. The pollutant-releasing device of the model test mainly included a CO 2 cylinder, a CO 2 two-stage pressure-reducing valve (better pressure stabilizing effect), and a vortex flow meter. In the experiment, CO 2 gas was released at about 100 m away from the outlet into the outflow tunnel so as to ensure that CO 2 was well mixed at the outlet section. The locations of measuring points in #1 tunnel are shown in Figure 15. A comparison between calculated results of Equation (5) and simulated values of φc is given in Figure 14. The goodness of fit Rnew is 0.99. The circulating air mixing ratio model (Equation (5)) was validated by the measured data of two actual tunnels with small spacing. Both of these tunnels are one-way and double-lane tunnels with a horseshoe-shaped cross section. The structural parameters of the two tunnels are shown in Table 1. Table 1. Structural parameters of the tunnels. The field measurement was conducted on a windless day. Because of the toxicity of CO, this test replaced CO with CO2, which is less harmful, cheaper, and easier to obtain and measure. The pollutant-releasing device of the model test mainly included a CO2 cylinder, a CO2 two-stage pressure-reducing valve (better pressure stabilizing effect), and a vortex flow meter. In the experiment, CO2 gas was released at about 100 m away from the outlet into the outflow tunnel so as to ensure that CO2 was well mixed at the outlet section. The locations of measuring points in #1 tunnel are shown in Figure 15.  The jet velocity of the jet fan was regulated in the tunnel in order to change inlet and outlet air velocity. The air velocity of the inlet and outlet tunnel was obtained by a Testo425 thermal anemometer (Testo AG, Lenzkirch, Germany). In fact, the air velocity of a single point on the cross section of the tunnel cannot represent the average air velocity of the entire cross section. Thus, several air velocity measuring points were set on the cross section according to fluid's characteristics of velocity distribution on the cross section of the tunnel, with the distribution of the serial numbers and the effective areas being shown in Figure 16a. Therefore, the weighted average value of the effective areas covered by the air velocities at the measuring points was the average air velocity u of the cross section of the tunnel. In other words, u = ΣV j ∆F j ΣF , where ∆F j and V j are the area of each part and the air velocity of the measuring point within each area. Obtaining the CO 2 concentration of the inlet and outlet tunnel by using CO 2 concentration transmitters, the arrangements of concentration measurement points are shown in the Figure 16b. As for #1 tunnel, H = 1.5 m, X = 7.0 m; as for #2 tunnel, H =1.8 m, X = 5.5 m. The arithmetic mean value of the CO 2 concentrations at these measuring points is the average concentration of the cross section of the tunnel. air velocity measuring points were set on the cross section according to fluid's characteristics of velocity distribution on the cross section of the tunnel, with the distribution of the serial numbers and the effective areas being shown in Figure 16a. Therefore, the weighted average value of the effective areas covered by the air velocities at the measuring points was the average air velocity u of the cross section of the tunnel. In other words,

No. Length L/m Cross-Sectional Area A/m 2 D/m d/m Δl/m Front Slope of Tunnel
, where ΔFj and Vj are the area of each part and the air velocity of the measuring point within each area. Obtaining the CO2 concentration of the inlet and outlet tunnel by using CO2 concentration transmitters, the arrangements of concentration measurement points are shown in the Figure 16b. As for #1 tunnel, H = 1.5 m, X = 7.0 m; as for #2 tunnel, H =1.8 m, X = 5.5 m. The arithmetic mean value of the CO2 concentrations at these measuring points is the average concentration of the cross section of the tunnel. Based on the measured air velocity and CO2 concentration, φc can be calculated according to the following formula: In the formula, Cin and Cout is the CO2 concentration of the cross section of the inlet and outlet tunnels, respectively; C0 refers to the local concentration of CO2.
Meanwhile, φc under different uin/uout were predicted according to the relation model established in this paper as well as the model proposed by Tan [13]. The prediction and test results are shown in Figure 17. It shows that, calculated by measured data, φc gradually increases with the increasing of uin/uout. The variation tendency of the φc value predicted by Tan's formula agrees with the measurement; however, the predicted φc value of #1 twin tunnels with staggered inlet and outlet is generally smaller than the measured value, and the predicted φc value of #2 twin tunnels with aligned inlet and outlet air tunnel portals is generally bigger than the measured value. Hence, it can be concluded that Tan's formula is only applicable to tunnels the entrance slope gradient of which is 90° and have aligned inlet and outlet air tunnel portals, not to tunnels with staggered inlet and outlet or with relatively small entrance slope gradients. Unlike Tan's formula, the φc value of #1 twin tunnels with staggered inlet and outlet predicted by the model built in this paper agrees with the variance tendency of the measured value, showing a good numerical match, while the maximum relative error is only 7.2%. Compared with Tan's formula, the φc value of #2 twin tunnels with staggered inlet and outlet predicted by this model is closer to the actual measurement. Based on the measured air velocity and CO 2 concentration, ϕ c can be calculated according to the following formula: In the formula, C in and C out is the CO 2 concentration of the cross section of the inlet and outlet tunnels, respectively; C 0 refers to the local concentration of CO 2 .
Meanwhile, ϕ c under different u in /u out were predicted according to the relation model established in this paper as well as the model proposed by Tan [13]. The prediction and test results are shown in Figure 17. It shows that, calculated by measured data, ϕ c gradually increases with the increasing of u in /u out . The variation tendency of the ϕ c value predicted by Tan's formula agrees with the measurement; however, the predicted ϕ c value of #1 twin tunnels with staggered inlet and outlet is generally smaller than the measured value, and the predicted ϕ c value of #2 twin tunnels with aligned inlet and outlet air tunnel portals is generally bigger than the measured value. Hence, it can be concluded that Tan's formula is only applicable to tunnels the entrance slope gradient of which is 90 • and have aligned inlet and outlet air tunnel portals, not to tunnels with staggered inlet and outlet or with relatively small entrance slope gradients. Unlike Tan's formula, the ϕ c value of #1 twin tunnels with staggered inlet and outlet predicted by the model built in this paper agrees with the variance tendency of the measured value, showing a good numerical match, while the maximum relative error is only 7.2%. Compared with Tan's formula, the ϕ c value of #2 twin tunnels with staggered inlet and outlet predicted by this model is closer to the actual measurement.

Conclusions
(1) For closely spaced twin tunnels with staggered inlet and outlet, the air flow discharged from their air-exhausting tunnels would diffuse transversely toward the air-inlet tunnels under the influence of the negative-pressure zone at the inlet, even causing inverse flow and forming circulating air. Parameters including inlet air velocity uin, outlet air velocity uout, lateral distance of the tunnel portals d, and staggered distance ∆l were all important factors influencing the circulating air.
(2) Both reducing inlet air velocity uin and increasing outlet air velocity uout and lateral distance d can reduce the impact of the negative-pressure zone at the tunnel entrance on the jet flow structure at the tunnel exit, thus weakening the circulating air.
(3) When the inlet is placed behind or aligned with the outlet (staggered distance ∆l ≤ 0), under the influence of the negative-pressure zone at the inlet, the jet flow at the outlet deflects entirely. The shorter the longitudinal distance between the inlet and the outlet (the larger the ∆l), the more the jet flow deflects; the quicker the jet flow maximum velocity umax in the symmetry plane decreases, the higher the circulating air mixing ratio grows.
(4) When the inlet is before the outlet (∆l > 0), the inlet's negative-pressure zone hardly influences the decreasing law of maximum velocity umax in the symmetry plane of jet flow between the inlet and the outlet but changes the decreasing law of the jet flow's umax in other zones. With the increase of longitudinal distance between the inlet and the outlet (increase of ∆l), the decreasing velocity of umax manifests a law that is at first increasing and then decreasing. Thus, the circulating air mixing ratio also increases first and then decreases.
(5) The quantitative relation formula of the circulating air mixing ratio for small-spaced twin tunnels with staggered inlet and outlet with air velocity ratio uin/uout, lateral distance of the tunnel portals d, staggered distance ∆l, and tunnel portal hydraulic diameter D is created by sections. The predictions show a good correlation with measurements. This formula supplements the limitation of the circulating air mixing ratio prediction formula which can be applied only in tunnels with aligned inlet and outlet. (6) According to the relation model, predictions in the present research have a good correlation with measured results for tunnels with entrance slope gradients of 15° but do not relate well with those of 20°, which indicates that the front slope gradient of the tunnel portal is also one of factors affecting the circulating air mixing ratio. More in-depth research and further discussions of this will be conducted in the future.

Conclusions
(1) For closely spaced twin tunnels with staggered inlet and outlet, the air flow discharged from their air-exhausting tunnels would diffuse transversely toward the air-inlet tunnels under the influence of the negative-pressure zone at the inlet, even causing inverse flow and forming circulating air. Parameters including inlet air velocity u in , outlet air velocity u out , lateral distance of the tunnel portals d, and staggered distance ∆l were all important factors influencing the circulating air.
(2) Both reducing inlet air velocity u in and increasing outlet air velocity u out and lateral distance d can reduce the impact of the negative-pressure zone at the tunnel entrance on the jet flow structure at the tunnel exit, thus weakening the circulating air.
(3) When the inlet is placed behind or aligned with the outlet (staggered distance ∆l ≤ 0), under the influence of the negative-pressure zone at the inlet, the jet flow at the outlet deflects entirely. The shorter the longitudinal distance between the inlet and the outlet (the larger the ∆l), the more the jet flow deflects; the quicker the jet flow maximum velocity u max in the symmetry plane decreases, the higher the circulating air mixing ratio grows.
(4) When the inlet is before the outlet (∆l > 0), the inlet's negative-pressure zone hardly influences the decreasing law of maximum velocity u max in the symmetry plane of jet flow between the inlet and the outlet but changes the decreasing law of the jet flow's u max in other zones. With the increase of longitudinal distance between the inlet and the outlet (increase of ∆l), the decreasing velocity of u max manifests a law that is at first increasing and then decreasing. Thus, the circulating air mixing ratio also increases first and then decreases.
(5) The quantitative relation formula of the circulating air mixing ratio for small-spaced twin tunnels with staggered inlet and outlet with air velocity ratio u in /u out , lateral distance of the tunnel portals d, staggered distance ∆l, and tunnel portal hydraulic diameter D is created by sections. The predictions show a good correlation with measurements. This formula supplements the limitation of the circulating air mixing ratio prediction formula which can be applied only in tunnels with aligned inlet and outlet. (6) According to the relation model, predictions in the present research have a good correlation with measured results for tunnels with entrance slope gradients of 15 • but do not relate well with those of 20 • , which indicates that the front slope gradient of the tunnel portal is also one of factors affecting the circulating air mixing ratio. More in-depth research and further discussions of this will be conducted in the future.

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