A Two-phase Model of Air Shock Wave Induced by Rockfall in Closed Goaf

In this paper, a two-phase model of air shock wave induced by rock-fall was described. The model was made up of the uniform motion phase (velocity was close to 0 m·s-1) and the acceleration movement phase. The uniform motion phase was determined by experience, meanwhile the acceleration movement phase was derived by the theoretical analysis.A series of experiments were performed to verify the twophase model and obtained the law of the uniform motion phase. The acceleration movement phase was taking a larger portion when height of rock-fall was higher with the observations. Experimental results of different falling heights showed good agreements with theoretical analysis values. Computational fluid dynamics (CFD) numerical simulation had been carried out to study the variation velocity with different falling height. The two-phase model could provide a reference and basis for estimating the air shock waves' velocity and designing the protective measures.


Introduction
In recent years, the failures of rock-roof and the collapse of surface happened more and more frequent which have severely damaged to people and mining [1][2][3].In these mine disasters, the air shock wave is one of typical forms with the characteristic of high velocity, high pressure, huge impact force and short effective duration [4].
With the purpose of accurately measuring the size and impact of this disaster, various models have been built to predict the velocity of the air shock wave [5][6].Chen et al. [5] modeled the shock process as compression process, which has a higher velocity than the reality.Zheng et al. [6] proposed that the shock process combined the compression and airflow process, but the model was not given the formulation.In order to estimate the velocity of air shock wave induced by rock-fall more accurately and conveniently, the two-phase model is raised.Then according to the prediction of air shock waves' velocity, the protective measures was designed to reduce the damage.
Experimental research-the most common method-can not only obtain the true value of the air shock waves' velocity, but also analyze the influencing factors of rock-fall height.To that end, Janovsky et al. [7] used methane and air as the tested material to analyze the velocity and their impact on the walls whose thickness is 14 and 29 cm thick walls in the process of exploration.Mixture air have also been adopted as pure air in many experiments [8][9][10].Xing et al. [11] revealed complex influences of goaf with different caved areas on the velocity and pressure of the air shock wave in different mining faces.Pennetier et al. [12] performed experiments of the air shock flushed the tunnel and presented evolution of the pressures and velocities of both the inlet and outlet tunnel in detail with the measurements on different blasted volume.Moreover, the test instruments used in experiments were not the same, such as thermal anemometer [13], different-pressure anemometer [14], ultrasonic anemometer [15].Ovink et al. [16] used the hot wire hot film Anemometer as a test instrument to study, the results acting conform to the actual values during travel over experimental material.
Euler-Euler CFD simulation as a common numerical method, has been widely applied to the simulation of gas flows [25][26][27].Zhu et al. [28] modeled 3D gas flows through Euler-Euler CFD simulation.Despite 3D models [29,30] were more precise than 2D ones [31,32], 2D models could still provide enough information and save time tremendously.Thus, in the case of the flow regime was not complex, a 2D model could be used to understand expecting results [33].Spicka et al. [34] presented a numerical simulation of gas flows by a two-phase model using a 2D model in the nuclear and car industries.The results show good agreement between experimental values and simulation ones.Furthermore, combination of CFD and user defined function (UDF) [35] has been carried out.Qin et al [36] used CFD and UDF to simulate small-scale experiments which better define movement forms and properties of mixture gas.
In this paper, we presented a two-phase model for air shock wave induced by rock-fall in closed goaf which could to predict the velocity more efficiently.Specifically, the model was composed of the uniform motion phase and the acceleration movement phase.The uniform motion phase was determined by experience, and the acceleration movement phase was derived by theoretical analysis.
Then a series of experiments were conducted to verify two-phase model and obtained the law of the uniform motion phase.CFD and UDF were used to stimulate the experiments by 2D models.Finally, we evaluated the effectiveness of our model and compared its performance with experimental and numerical simulation values.As a result , the two-phase model could estimate the velocity of air shock wave induced by rock-fall accurately and conveniently.Meanwhile, The two-phase model could provide a reference and basis for estimating the air shock waves' velocity and designing the protective measures.

Basic Assumptions
Based on the model of pump [5] and streaming [6], a model of two-phase was raised to estimate air shock wave induced by rock-fall.The model was composed of the uniform motion phase and the acceleration movement phase.The speed of the uniform motion phase was close to 0 m•s -1 on account of the pressure difference.The acceleration movement phase was derived by computational fluid mechanics and classical mechanics.Then, the assumptions of the acceleration movement phase were raised as follows: (a) Air was considered ideal gas and characterized by laminar flow conditions; (b) The process of air shock wave was identified as an isothermal process; (c) The goaf was considered as closed and the volume difference of it was neglected.

Modeling Equations
The governing equations were used in our 2D, isothermal model would be described.The law of conservation of momentum [37,38] is shown as follows: 2 3 Where ρ is the density of air, (Kg•m -3 );F is the body force, (N); μ is the viscosity of air, (Pa•s); v is the velocity of air, (m•s -1 ) .
As results of the above assumptions, values for μ is 0 and ρ is got from the ideal law, / p RT ρ = , which R is the universal gas constant and T is the absolute temperature respectively.The weight of falling roof was acted as the body force.Eq. ( 1) is simplified as: Where ρd is the density of the rock-fall, (kg•m -3 ); N is thickness of the rock-fall, (m); Pn is the pressure of the former in the goaf, (Pa); Pn-1 is the pressure of the latter in the goaf, (Pa); vz is the vertical velocity of air in the goaf, (m•s -1 ); g is the gravitational acceleration, (m•s -2 ).
In our model, the ideal gas state equation [38] is described: Where Vn is the falling volume of the latter, (m 3 ); Vn-1 is the falling volume of the former, (m 3 ).
Owing to the condition of the model, the initial falling height of the acceleration movement phase was not 0 m but the final height of the uniform motion phase.The volume can be calculated from the following equation: Where H is the height of goaf, (m); Hm is the final height of the uniform motion phase, (m); Hm is the falling height of the acceleration movement phase, (m); Hn is the latter falling height of the acceleration movement phase, (m); Hn-1 is the former falling height of the acceleration movement phase, (m).
On the basis of the Newton's Law, the velocity and the falling height can be obtained: ) Where vn is the latter velocity of air in the goaf, (m•s -1 ) ; vn-1 is the former velocity of air in the goaf , (m•s -1 ) ; ∆t is the time step, (s).The model is controlled by the three unknowns (P,v,H) in nature, and the three unknowns to be solved depend on one another, as shown in Eq. ( 6).Nowadays the implicit pressure explicit saturation(IMPES) scheme [39] was proposed to solve this problem.Meanwhile, that was an operator decomposition technique with which the unknowns ( P,H) is decoupled from the velocity v by lagging one time step behind calculations.This allowed for solving ( Pn,Hn) implicitly that gived vn-1 .Next ( Pn,Hn) was used to solve vn explicitly.The IMPES scheme was listed in Fig. 1[40].Because the outlet was linked with the goaf, the air shock wave in the outlet is applied to the Bernoulli equation [41], hence the air shock waves' velocity in the outlet can then be calculated from the following equation: Where vb is the velocity of the outlet, (m•s -1 ); pbn is the later pressure of the outlet, (Pa); pbo is the initial pressure of the outlet, (Pa); ρn is the density of air in the goaf, (kg•m -3 ); ρL is the density of air in the outlet, (kg•m -3 ).
According to the on-site experiment, the area of the falling roof was smaller than the area of goaf.Based on the air quantity balancing law, the pressure of the outlet can be received: Where S is the sectional area of goaf, (m 2 ); Sa is the sectional area of the falling roof, (m 2 ); Sb is the sectional area of the outlet, (m 2 ).

Experimental Results and Discussion
The straight line and a rapid growth curve.It could also be said that the model of twophase was possible to estimate the air shock waves' velocity.
At the phase of uniform motion, as the tank was closed, a huge pressure difference would be engendered in the tank with the iron-on salver's way down and prevented the air shock wave from the tank.At the phase of straight line, reducing the bound of pressure difference was a major factor in the rapid growth of velocity.
Therefore, the realistic air shock wave was made up of the phase of acceleration movement and uniform motion.Meanwhile, the peak values of three velocity curves were not in direct proportion to its release heights, it is known that the closed goaf could reduce the air shock waves' hazards but not remove it.According to the experiment, huge pressure difference was generated with the iron-on salver falling.The air was regarded as incompressible due to the process of the experiment was transient.The upper part of the iron-on salver raised the zone of negative pressure and the bottom presented positive pressure, the outlet kept normal pressure [42].The air of bottom would burst in the upside and outlet.If the pressure difference between upper apart and bottom could absorb enough air, the siphon phenomenon would form on the downside.As a result of this, the stage of uniform motion would appear.Otherwise, the stage of acceleration movement would turn up.
From Fig 5, it can be known that the velocity of 0.9 m and 1.2 m show the same trend as 0.6 m, only that the phase of uniform motion is more obvious and the peak value is much higher.The velocity of 1.2 m release height was largest, followed by 0.9 m and 0.6 m, were 21.3488 m•s -1 , 26.1003 m•s -1 and 28.4508 m•s -1 respectively, and the time of phase of uniform motion were 0.23 s, 0.27 s and 0.31 s respectively.
The turning points of uniform motion phase and acceleration movement phase respectively were located 0.2645 m, 0.3645 m, 0.4805 m in the above release heights on account of Newton's Laws.From Eq.( 6), the pressure difference was inverse proportional to the falling height.As a result of this, the phase of the acceleration movement gets enlarge.This indicates that the phase of acceleration movement is taking a larger portion and uniform motion is smaller when the height of rock-fall is higher.However, increasing the height of rock-fall was strengthening the role function time of the pressure difference.Hence, the higher falling roof was, the low the growth rate of the acceleration movements' phase.Because of occupied a larger portion, the higher falling roof was, the higher the air shock waves' velocity, but not in direct proportion to the falling roof.Hence, the higher falling roof was, the slower the growth rate of the acceleration movements' phase.Because of occupied a larger portion, the higher falling roof was, the higher the air shock waves' velocity was, but not in direct proportion to the falling roof.

The Uniform Motion Phase
From

The Acceleration Movement Phase
According to experimental set-up and Eq.( 10), the coefficients for both conditions, the different falling heights were summarized in Table 1.Used to MATALAB to solve Eq.(1) to Eq.( 9 existing, led to a large impact on absolute error of initial points.The higher the air shock waves velocity, the less the impact of external wind.As results of this, the prediction precision of the two-phase model was increased.This results show good agreement between experimental values and theory results.This demonstrates despite the small differences due to experimental error, that the model of two-phase provides an accurate way to evaluate it.
The calculation formula of the two-phase model follows: ( )

Numerical Simulation
Computational fluid dynamics (CFD) [43] was an important component of fluid mechanics that applied numerical modeling to analyze the flow problems.User defined function(UDF) [44] approved the user to write C programs that could be linked with CFD solver to promote the performance of the code.Meanwhile, the simulation process was governed by the fundamental laws of fluid mechanics.

Simulation Models
In CFD simulation, air movement in the goaf and diffusive motion process followed Darcy's law, Hooke's law and Langmuir's equation respectively.Goaf and the outlet face were regarded as rigid wall , air was regarded as ideal gas, and roof falling process was regarded as an isothermal process.In this experiment, a laminar model was used to solve the process of air movement.Besides, the rock-fall was defined as moving body, and the UDF was used to compile the rock-fall.In order to more accurately simulating rock falling, the option of Six DOF was selected to define the rock-fall.The parameters of rock-fall came from the experimental apparatus.Then air outlet boundary of numerical model setting: PRESSURE-OUTLET; boundary conditions obtained from field.And three different groups including 0.6 m, 0.9 m, 1.2

Simulation Results and Discussion
The simulation results of velocity using three groups with different falling height were shown in Considering the objective impact force and energy were ignored, the influence of the pressure difference was enlarged.Consequently, the time of uniform motion phase got extended.Therefore, the location of the turning points were much lower than the experimental results'.The air was flushed to outlet immediately on the acceleration movement phase, considering the hypothesis of the air was identified as an incompressible gas.Hence, the ratio of numerical simulations' curves were much larger than the experimental results' on the acceleration movement phase.This analysis shows that, the air shock wave induced by rock-fall was composed of the uniform motion phase and the acceleration movement phase.It verified the effectiveness of the two-phase model.

Conclusions
In this paper, we had proposed and verified the semi-empirical model of the air shock waves' velocity induced by rock-fall using theoretical analysis and experiments.
We had evolved a two-phase model which was composed of the uniform motion phase(speed was close to 0 m•s -1 ) and the acceleration movement phase through theoretical analysis.Satisfactory agreement had surveyed in the field of the velocity and the variation trend of the air shock wave between the model and experimental results, which predicted the correctness and effectiveness of the model.By

Figure 1 .
Figure 1.A solution procedure for the two-phase model based on operator decomposition An experimental set-up was built to investigate the air shock wave induced by rock-fall in closed goaf, as shown in Fig 2. The experimental set-up consist of the air shock wave generator (Fig 3), the hot wire hot film Anemometer (HWFA) of IFA300 (Fig 4), data collection system and other components.

Figure 3 .
Figure 3.The air shock wave generator.Figure 4. The hot wire hot film

Figure 4 .
Figure 3.The air shock wave generator.Figure 4. The hot wire hot film Anemometer of IFA300.
velocity of the air shock wave after releasing the iron-on salver in the different falling heights were shown in Fig 5.As shown in Fig 5, the three curves followed the same trend, it was known that the velocity curves were composed of a Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 23 January 2017 doi:10.20944/preprints201611.0083.v2

Fig 5 ,Where
it could be seen that there was a linear relation between the height and time of the uniform motion phase, and the falling height.The relation between both variables could be expressed as: Hm corresponds to the height of the uniform motion phase, in meters (m); a, b and c are the parameters of the curve; and H is the falling height, in meters (m).The values for a, b and c for the uniform motion phase are 0.11025, 0.19551 and 0.0866761 m respectively, with a coefficient of determination R 2 equal to 0.99.
) by the above coefficients, the experimental values and theoretical values of the air shock waves' velocity in the outlet could be appearing on Fig 6.As shown in Fig 6, the experimental values and theoretical values followed the same trend.The absolute errors between experimental results and theory analysis were shown in Fig 7, for different falling heights.The absolute errors were almost less than 10%, and the absolute errors of cut-off velocity were 1.87%, 0.50% and 3.26% respectively.Because the initial velocity was too little and the external wind was Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 23 January 2017 doi:10.20944/preprints201611.0083.v2

Fig 8 .
As shown in Fig 8, the simulation results of air shock wave were composed of a straight line and a rapid growth curve.It can be known that the air shock waves' velocity of 0.9 m and 1.2 m show the same trend as 0.6 m from Fig 6.Meanwhile, the time of the uniform motion phase were 0.32 s, 0.39 s and 0.46 s respectively, and the maximum velocity of three groups were 13.2089 m•s -1 , 19.6315 m•s -1 and 25.0314 m•s -1 respectively.Compared with the simulation values of the air shock wave with the same falling height in experimental results, the values in experimental results were larger than that in numerical simulations.And the duration time of the uniform motion phase were less than the numerical simulation results.

Figure 8 .Figure 9 .
Figure 8. Simulation results of velocity using three groups with different falling roof.(A) 0.6 m. (B) 0.9 m. (C) 1.2 m.The velocity distribution at the turning point by three groups were shown in Fig 9.The maximum velocity were mostly distributed at the sides of the rock-fall, such as falling height with 0.6 m in Fig 9(A), 0.9 m height in Fig 9(B) and falling height with 1.2 m in Fig 9(C).Meanwhile, the values at upper part of rock-fall were larger than the values at bottom of rock-fall.The appearance of phenomenon might interpretation the role of pressure difference.

Table 1 .
The coefficients for acceleration movement phase

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 23 January 2017 doi:10.20944/preprints201611.0083.v2 m
were selected for the model verification.Group number represented the height of rock-fall.For example, 0.6 m represented that height of rock-fall was 0.6 m.