Seismic Response Analysis of Deep Underground Roadways and Coal Pillars under the Influence of the Adjacent Goaf

In this work, a numerical study is conducted on the seismic response of deep-buried roadways in coal mines under the influence of goafs, and a 3D numerical model of the seismic response simulation of deep-buried roadways is established using the coupling model of the finite difference method and the distinct element method.*is model simulates the seismic response of different coal pillar widths and the seismic conditions of the deep-buried roadways under the influence of the adjacent goafs. *e deformation, stress distribution, and plastic area distribution of roadways and coal pillars are systematically studied, and the situations under the static load and the roadways, which are not affected by the goafs, are compared and analyzed. A reasonable width of the coal pillar is proposed on the basis of the stability of the roadway and the coal pillars. In the end, suggestions for the reasonable setting of coal pillars under seismic load are provided.


Introduction
e earthquake disaster causes serious damage to ground buildings and brings about serious effects to underground structures [1][2][3]. e deep underground mining often induces seismic activity, which is called mine earthquake [4][5][6][7]. Natural and mine earthquakes are not essentially different [8,9]. Natural earthquakes have caused different degrees of damage to coal mines [10,11]. Underground coal mines are gradually affected by mine earthquakes with increased coal mining depth [12,13]. In the limited deep space, mine earthquakes may cause a series of other coal mine disasters, such as rock explosion, coal and gas outburst [14], and water penetration [15][16][17][18][19][20]. erefore, the response of deep underground coal mine structure should be analyzed under the seismic load. e gob-side entry driving is a long-walled coal mining often used in the form of exploitation [21][22][23][24][25], and the roadway driving with coal pillars can effectively isolate the goaf to prevent water and harmful gases in the goafs into the roadway [26,27]. Feng et al. [28] established the relationship between microseismicity and the risk of rockburst and then proposed an early warning technique based on it to predict the rockburst in tunnels. By analyzing the monitored microseismic data, this method can provide a real-time early warning of the rockburst risk during excavation of the tunnel. e application of this method in some engineering cases shows that it can well predict the development of rock burst [29,30]. Li and Chai [31] studied the relationship between the width of coal pillars and the susceptible coal explosion in deep-buried roadways and proposed a reasonable coal pillar width setting to reduce the risk of the coal explosion. Considering that coal pillars play a supporting role in the active coal mining space and are subjected to different forms of dynamic action, the current mining stress of the coal mining face is the main dynamic form. Coal pillars are also affected by the process of isolating the goaf and the overlying rock strata sinking the goaf during the compaction process. Coal pillars are affected by dynamics, which often cause large deformations and bring about difficulties to the coal mining. Zhao et al. [32] and Wang et al. [33] analyzed the stress changes and the damage risks of coal pillars under dynamic and static loads with specific cases and proposed reasonable coal pillars under corresponding geological conditions. In previous studies, underground chambers or roadways often produced a certain amount of deformation under the action of earthquake loads. Wang et al. [34] analyzed the seismic response of underground gas storage salt caverns, which produce certain deformation and damage in the seismic response. In turn, the seismic response has a certain effect on the safety of gas storage salt caverns. is deformation does not pose a great threat to the roadway in the coal mining space of the coal mine, and the effect of the dynamic stress on the roadway far from the goaf is small [35][36][37]. However, the deformation of the roadway formed by the gob-side entry driving due to the isolation of coal pillars near the goaf is often severe. In the study of the dynamic response of isolated coal pillars under the action of the goaf, the goaf is often regarded as a continuous medium [38][39][40], which is often inconsistent with the actual situation of broken rocks in the adjacent goaf.
In this study, the finite difference method (FDM) and the distinct element method (DEM) are coupled, and a continuous medium model of isolated coal pillars and formations is established. A discrete medium model of the adjacent goaf is established using the DEM. e effect of the adjacent goaf under seismic load on the coal pillar and the response of coal pillars with different widths under the seismic load are analyzed. ese results should provide suggestions for the reasonable setting of coal pillars under the seismic load.

FDM and DEM Constitutive Models
e computational resource required to model coal strata sequence using a distinct element code similar to PFC is enormous, especially in dynamic calculation. e goaf is filled with broken gangue, and the discrete medium in the goaf belongs to unbonded granular materials. erefore, using PFC to simulate the goaf is reasonable. e Mohr-Coulomb model is widely used in the seismic analysis in the nonlinear dynamic analysis in geotechnical engineering [41][42][43]. e continuum model is selected for the coal pillar and the adjacent complete stratum, and the Mohr-Coulomb model is selected for the constitutive relationship of the model. e contact connection model is used to simulate the discrete medium in the goaf without considering the internal fragmentation of the particles. e solid element (zone) in the FDM and the wall in the DEM are used for coupling. e coupling of FDM and DEM is required to reduce computation time while ensuring solution accuracy.
A model (Figure 1) is established to verify the feasibility of the coupling scheme. e length, width, and height dimensions of the overall model are 160 m × 5 m × 15 m in the x, y, and z directions in the software. e hollow area of 150 m × 5 m × 10 m is filled with the DEM and the FDM. e mechanical parameters related to the FDM and the DEM are shown in Table 1. FDM parameters are obtained from coal mine strata. e DEM refers to the data in literature [44]. In the coupling model of the FDM and the DEM, the DEM particle size is 0.2-0.5 m. A total of 24,816 particles (balls) are generated, and the number of physical units (zone) is 4,500. e number of solid elements (zone) is 12,000 in the model with only continuum.
After the model is established, the parameters are assigned, and the bottom and all sides of the model are fixed.
No stress or restriction is observed in the upper part. e self-weight stress balance is calculated, and the model displacement and the velocity are initialized to 0. e acceleration time range ( Figure 2) is entered from the bottom of the model and provides artificial waves for the software. e acceleration equation for the wave is where A is the amplitude and f is the frequency. e free-field boundary is applied around the model to reduce the reflection of the acceleration wave. e static boundary conditions are removed, and the dynamic response of the model is verified under the influence of no surrounding stress. e dimensional effects are considered. e distance at the y direction is short, and the free-field boundary is observed at both ends. e acceleration in the y direction is set to 0, and the acceleration in the x and z directions is entered in accordance with the actual size. e model dynamic calculation time is set to 10 s. e coupling of FDM and DEM under self-weight stress (Figures 3(a) and 3(b)) and dynamic load (Figures 3(c) and 3(d)) conditions and the vertical stress in the left pillar of the continuous medium models (Figures 3(b) and 3(d)) is extracted to verify that the response of the bulk and the continuum media to the cylinder differs. e coordinates of the pillars are x(0-5), y(0-5), and z(0-15), and the stress is extracted from x � 4.9 in the x plane. Considering that 5.0 is affected by displacement, some areas cannot be extracted due to stress. As shown in Figure 3, the vertical stress of the continuum (Figure 3(b)) shows a linear change along the zaxis as a whole under the condition of self-weight stress. However, the coupling model of the FDM and the DEM (Figure 3(a)) is affected by the contact stress of the loose body. Moreover, the local area is no longer linear, shows a fluctuation in the whole plane, and does not change uniformly in the direction of the y-axis. ese observations show that the DEM and the continuum have completely different forces on adjacent areas.
is difference is also evident under the dynamic load conditions. As shown in   Advances in Civil Engineering y direction is low at both ends, high in the middle, linear, and gradually stable from both ends to the middle. At the position of y � 2.5, the stress in the z direction increases first and then decreases. In the FDM and the DEM coupling model (Figure 3(c)), an irregularity between the self-weight stress and the continuous medium is observed. erefore,  Advances in Civil Engineering the dynamic load response of coal pillars under this discrete medium, which is different from the continuous medium, should be studied.

Geological Background and Site Details.
e Kongzhuang Coal Mine in Xuzhou, China, was selected as a numerical example. e longwall panels (7433 panels and 7435 panels) with width of 200 m and length of 100 m are selected. e coal seam thickness of the working face is 4.20 m∼5.10 m, and the average thickness is 4.60 m. e immediate roof of the coal seam is sandy mudstone, with a thickness of 0.60∼11.19 m, an average thickness of 3.90 m, gray-black, thin-layered, uneven sand content, and more plant fossils. e main roof is composed of siltstone and medium-grained sandstone. e thickness of siltstone is 1.87∼3.21 m with an average thickness of 1.50 m. e thickness of medium-grained sandstone is 5.23∼9.81 m with an average thickness of 6.66 m. e siltstone is gray-white; the main components are quartz, feldspar, calcareous cement, compact, fine-grained sand-like structure, and nearly horizontal bedding; for details, see Figure 4.

Modeling and Solutions.
e numerical model is established on the basis of the tail entry of the 7435 working face of the Kongzhuang Coal Mine. e 7435 working face is the first working face buried more than 1000 m deep in Kongzhuang Coal Mine. e plane relationship between 7433 and 7435 working faces is shown in Figure 5. During mining at the 7433 working face, if the 7435 working face is excavated, then an isolated coal pillar should be set between two working faces. e length along the dip of the two working faces is 150 m. us, the length along the dip of the 7433 goaf is 150 m, and the tail entry of the 7435 working face is shown in Figure 5. e width and the height of the roadway section are 5 and 4 m, respectively. e 3D sizes of the comprehensive strata histogram and model are shown in Figure 4. e size of the model in the x direction is 350 m, and the sizes in the z and y directions are 200 and 5 m, respectively. e layout of monitoring points around the roadway is shown in Figure 4 to monitor the displacement and the stress of the roadway and the coal pillars. e height dimension of the goaf, that is, the caving zone, uses statistical theoretical analysis to arrive at the regression formula of the fall zone on the basis of the highly measured data of the caving zone of a large number of goafs in China and the United States made by Bai et al. [45].
where H c is the height of the caving zone (m), h is the height mining (m), and c 1 and c 2 are the strata strength parameters. e values of c 1 and c 2 can be determined on the basis of literature [45]. e 7433 working face has a mining height of 4.5 m, and the calculated caving zone height is approximately 11.2 m and is regarded as 11 m. e working face has a length of 150 m. erefore, the size of the goaf is (i) X × Y × Z, 150 m × 5 m × 11 m by using the DEM for modeling. e DEM model ( Figure 6) is built. e randomly generated geometric model is imported from the external software to simulate the broken rocks in the goaf realistically, and the spherical particles are filled into the geometric model to make the surface of the model as close to the actual one as possible under stable model conditions. For the physical and mechanical parameters of the two models, the FDM parameters are determined in accordance with the geological conditions of field engineering and the mechanical properties of the rock formations. e DEM model parameters are determined in reference to literature [44], and the physical mechanical parameters are presented in Table 2. e total numbers of model units and DEM are 124 250 and 3866, respectively. e size of broken rock mass is 0.3-2 m, and the particles of the broken rock mass are randomly generated in this range. e porosity of the model is 0.3. In this discussion, the effect of particle breakage and size is ignored. e position of the 7435 tail entry in the model moves to the left in the dashed box of Figure 4 to change the widths of the coal pillar to 5,7,9,11,13,15,20,25,30,35, and 40 m, which indicates the construction of 11 models.
e grid size of the model is dense around the adjacent roadway and the goaf. e grid size is far from the roadway, and the goaf is appropriately enlarged to improve the accuracy of the calculation without affecting the calculation results.
e 7435 working face of the Kongzhuang Coal Mine has an elevation of −1017.50 m to −883.80 m. is study ignores the effects of the stratigraphic dip. e working face is arranged horizontally. e working face and the roadway are arranged at 1000 m. e rock mass density and the upper vertical stress are valued at 2500 kg/m 3 and 25 MPa, respectively, and the horizontal stress ratio is 1.5 during the static calculation of the model. e front, rear left, right, and the bottom are fixed at the model static calculation boundary conditions. e dynamic boundary condition of the model is the free-field boundary condition to reduce the reflection of waves.
e mechanical parameters of coal and rock can be found in Table 2, and the DEM mechanical parameters are the same as those in Table 1. e two sides of the 7435 return airway are supported by bolts, with length, diameter, and spacing of 2.4 m, 20 mm, and 1 m × 1 m, respectively. Four anchors are present on each section with the yield load and pretightening force of 150 and 80 kN, respectively. e roof is supported by bolts and cables with spacing of 1 m × 1 m. e anchor cables are arranged in the middle of the bolts near the middle point of the bolster and the roof of the roadway with anchor cable diameter, length, yield load, and pretightening force of 17.8 mm, 8.2, 600 kN, and 150 kN, respectively.

Input and Processing of Seismic Waves.
In the calculation process, the Northridge, Loma, Kobe, and artificial seismic waves are used [46]. e first 20 of the four seismic waves are extracted as input dynamic loads, and the four seismic waves      Figure 9 shows the vertical displacement change curve of monitoring points 3 and 7, which is the largest displacement change of all monitoring points:  Figures 10 and 11 show the maximum main stress wave cloud diagrams of the roadway and coal pillars under the static load and Northridge seismic waves. More models of four seismic waves exist, and only the Northridge seismic wave is selected for comparative analysis because of space limitation. As shown in Figures 10 and 11, a low stress zone is formed around the roadway under static and dynamic conditions, and this result is due to the deformation of the excavated chamber to unload the pressure. ese findings are similar to those of previous studies. Under the static load, the stress concentration in coal pillars with width of 5 and 7 m is low. is case is conducive to the stability of the coal pillars. When the width is increased from 9 m to 20 m, the stress concentration in the coal pillar gradually increases. e stress value gradually increases from −10 MPa at 9 m coal column to −35 MPa when the width of the coal pillar is increased. e symbol here is opposite to the axis direction when the model is built; thus, this stress distribution is negative. In the range of 9-20 m, the high-stress areas formed on the sides of the goafs and the side of the roadway gradually converge, and this condition forms a stress core with a high value, which is extremely detrimental to the stability of the coal pillars, especially in the process of mining. No convergence occurs between the stress on the side of the roadway under the coal pillars with width of 20-40 m and the stress on the sides of the goafs due to the distance, which is beneficial to the stability of the coal pillars. However, the wide coal pillars reduce the rate of resource extraction especially for the scarce types of coal to improve the recovery rate. e use of narrow coal pillars entry protection is relatively economical.

Stress Analysis
Under the action of the Northridge seismic wave, the stress value of the whole model is lower than that of the static load. is condition indicates that the stress concentration in the coal pillars is not evident under the dynamic action. However, the stress concentration area of the whole coal pillars is still distributed on the solid coal side of the roadway and inside the coal pillars but may gradually spread to the top and the bottom. e high value of local stress is at most around 20 MPa, which is relatively safe for the stability of coal pillars. However, as mentioned earlier, the displacement of the coal pillar side of the roadway increases as the earthquake continues, and this condition may cause the risk of instability for coal pillars that are subject to continuous or frequent movement. erefore, the narrow coal pillars cannot ensure the safety of the coal pillars and roadway under the dynamic action and should be increased to 20 m or more to ensure safety. Figure 12 shows Advances in Civil Engineering 11 respectively. As shown in Figure 12, the vertical stress is higher than the vertical stress under 5 m wide coal pillars and is roughly distributed in the range of −25 MPa to −40 MPa. Under the 5 m wide coal pillar, the vertical stress in the coal pillars ranges from 10 MPa to −15 MPa. Given the plastic damage to the 5 m wide coal pillar, the overlying strata and roof can no longer be supported, and this condition shifts the stress. Figure 12 shows that the vertical stresses in the 5 and 20 m wide coal pillars are unevenly distributed in areas closer to the goaf. Evident unevenness is observed when the coal pillar is close to the goaf. With the increased distance from the goaf, the stress gradually becomes uniform, which indicates that the loose medium of the goaf causes uneven stress of the coal pillar in the close range from the goaf, and the influence of the position of the distance is gradually reduced. Figures 13 and 14 show the distribution of plastic regions under the static load and Kobe seismic waves, respectively. e blue area represents the area where no plastic damage has occurred, the red area represents the area where shear damage has occurred, and the purple area represents the area where stretching and shearing damages have occurred. As shown in Figures 13 and 14, the plastic area under the static load occurs around the roadway, and the plastic area under dynamic action expands along the coal seam in the horizontal direction with a striped distribution.

Distribution Analysis of the Plastic Zone
Under the static load, the area of the plastic zone in the range of 5-15 m is insignificantly reduced, but the plastic area within the coal pillars always runs through the entire coal pillars with the widening of the coal pillars. erefore,  Advances in Civil Engineering 13 increasing the width of the coal column from 5 m to 15 m does not reduce the plastic damage, which is affected by the goaf. When the width of the coal pillar increases to 20 m, the plastic area within the coal pillar is gradually reduced. For the coal mining, the coal pillar without plastic damage has higher carrying capacity than that with plastic damage and is relatively safe for coal mining. Wide coal pillars are often used for coal mining in the past. e distribution of plastic areas under the influence of seismic waves in a single roadway not affected by the goaf is shown in the no gob of Figure 14. A certain range of shear plastic failure areas occur around the roadway, the range of which increases with the continuous occurrence of the earthquake. e distribution of plastic regions under the action of Kobe seismic waves is different from that of static load. With increased width of the coal pillar, the length of the plastic area does not decrease significantly. is finding implies that the plastic area expands with the low strength of the coal seam with the continuous seismic movement, whereas the coal column, which plays a supporting role, completely causes plastic damage. As mentioned earlier, the amount of deformation of the coal pillar where plastic damage occurs gradually increases as the earthquake continues. After the plastic damage occurs, the carrying capacity of the coal pillar is greatly reduced. At this point, the risk can only be reduced by increasing the width of the coal pillar or by strengthening the support of the coal pillar. e coal pillars under the frequent or the continuous

Conclusions
e seismic response analysis of coal pillars and roadways of adjacent goafs was conducted using the numerical method by utilizing the return airway of the 7435 working face in the Kongzhuang Coal Mine as a case study. e goaf was established using a different bulk model from the past, which was more in line with the actual situation of coal mining. Under the FEM and the DEM, the seismic response analysis was performed on coal pillars and roadways under different coal pillar widths and four seismic waves. is numerical method can be used to determine the width of coal pillars in gob-side entry driving. e results show that the effect of the bulk medium on the surrounding continuous medium was simply verified because of the effect of the dispersal contact stress, which led to the nonuniform distribution of stress with the continuous medium part of the bulk body and was confirmed in subsequent calculations. For coal pillars and roadways in adjacent goafs that were subject to frequent or sustained dynamic load, a great difference between the static load and a single roadway not affected by the goaf was observed. e deformation of coal pillars and roadways increased with the continuous dynamic load. Narrow coal pillars showed this characteristic, and this trend was reflected in the coal pillar deformation. e stress concentration in the coal pillar tended to weaken compared with the static load with the continuous earthquake motion, and the high-stress area gradually spread to the top and the floors of the coal pillar. e range of plastic areas in roadways and coal pillars gradually expanded along the direction of weak coal seams as the earthquake motion continued. For the coal pillars and roadways in the adjacent goafs that were subject to frequent or sustained dynamic loads, the width of the coal pillar was increased to ensure the stability of roadways and coal pillars, and the width of the coal pillar should be calculated and determined according to the geological background and site.

Data Availability
All data, models, or code generated or used during the study are available in a repository or online in accordance with funder data retention policies.

Conflicts of Interest
e authors declare that they have no conflicts of interest.