Stress Distribution and Mechanical Behaviour of Rock Mass Containing Two Openings Underground: Analytical and Numerical Studies

In this study, stress solution for rock mass containing two rectangular openings was calculated based on the Schwarz alternating method to investigate the stress distribution in rock mass around openings with different layouts. In addition, large-scale numerical models were further established for the two-opening system by means of the PFC-FLAC coupling method, in which the stress evolution, failure patterns, and acoustic emission (AE) events were presented. With the combination of analytical and numerical solutions, the interaction mechanism between two openings under different layouts was discussed from the perspective of stress and failure. The result shows that the confining stress within a certain range contributes to relieving tensile stress concentration around openings. The stress condition within the connecting area and coalescence pattern between two adjacent openings is dominated by their layout. Compared with small-size rock specimens in laboratory tests, the failure patterns around openings show a better agreement with the stress concentration characteristics determined by analytical stress solutions.


Introduction
There exist a large number of openings with different shapes and sizes underground in various mines, like roadways and tunnels in service and goafs waiting for filling or which have been abandoned, which may present potential environmental issues and security risks. Under high geostress, rock mass around openings or flaws may be subjected to high stress concentration, which may lead to initial failure triggering the instability of rock structures [1][2][3][4][5][6][7][8].
The analytical solution and failure behaviour for rock mass containing a circular opening have been comprehensively investigated by a wide range of literature. It has been reported that there are three kinds of failure, namely, tensile cracks, spalling cracks, and far-field cracks within the rock mass [9][10][11]. Under uniaxial compressive conditions without confining stress, the initial tensile cracks usually appear in the top and bottom of the openings, which are tensile stress concentration areas. Spalling cracks appear around the side walls of the openings because of high compressive stress concentration. Similar failure patterns also happened to rock mass containing a noncircular opening reported by other studies [12].
In addition, the analytical stress solution for plates containing a single opening has been also comprehensively investigated to explain the mechanical behaviour of rock mass around the opening by a lot of studies [13][14][15][16][17][18]. Based on the complex variable method, stress solution for elastic plates containing an opening with typical shapes, such as circular, elliptical, triangular, and rectangular shapes, has been calculated by Ukadgaonker and Awasare [19][20][21], Ukadgaonker and Rao [22], and Sharma [23]. Combining analytical solution and laboratory test, Wu and Ma [12] calculated the stress concentration factor on the periphery of a horseshoe-shaped opening and discussed its effect on the failure patterns of rock specimens. Tan et al. [24] considered six kinds of holes with typical shapes in practical rock engineering and discussed the shape effect on stress distribution and failure patterns around the openings based on the analytical solutions combined with numerical and experimental results.
When there exists more than one opening within an area underground, the adjacent openings may change the stress distribution around an opening significantly and thus lead to complicated failure patterns [25]. The interaction between two openings has also been discussed by many studies from the perspective of analytical stress solution. Zhang et al. [26,27] proposed a method to calculate stress solution for plates containing two or more openings with arbitrary shapes based on the Schwarz alternating principle. This method was further improved by Tan et al. [28] for a better solution accuracy, based on which the effects of opening layouts and confining stress on the stress distribution and failure patterns for two U-shaped openings and two arched openings were investigated, respectively. Besides analytical solutions, laboratory experiments for rock specimens containing two or more flaws or openings have also been widely carried out to understand the mechanical interaction mechanism between two defects. Zhou et al. [29] conducted uniaxial compression tests combined with acoustic emission (AE) and digital image correlation (DIC) techniques and then analyzed the fracture coalescence behaviour around marble specimens containing two rectangular openings. Zhao et al. [30] conducted a series of uniaxial compressive tests on rock specimens containing several circular openings and found that openings were strongly affected by each other with complicated failure patterns. Similar results were also presented by Lin et al. [31] when studying the crack initiation and coalescence patterns within rock specimens containing multiple openings under uniaxial compression. It can be concluded that the fracturing behaviour and specimen stability are significantly affected by opening layout including their distance and the connecting angle between them. However, in laboratory tests, the size of specimens is usually small relative to openings within them, which may lead to a strong boundary effect interfering with the mechanical behaviour of rock mass around openings. To address this problem, numerical modelling methods have been widely employed, which can avoid the boundary effect on interesting areas by establishing a large-scale numerical model. Among them, the finite element method (FEM) and finite difference method (FDM) are popular in terms of the solution for deformation problems in rock engineering but are not good at describing material fracturing behaviour. The discrete element method (DEM) is another leading numerical method dealing with rock failure problems, but DEM numerical modelling may give rise to huge calculation costs for the problem about a small interesting part in a large model with a great number of particles. In view of this, coupling the DEM method with the FEM method or FDM method is an interesting alternative promising the solution accuracy for mechanical behaviour of rock mass close to the opening within a large model with little effect from the model boundary. In particular, coupling simulation with particle flow code (PFC) and fast Lagrangian analysis of continua (FLAC) has been widely used in rock engineering [32][33][34].
In this study, analytical solutions as well as numerical simulations were conducted for a better understanding of the stress distribution and mechanical behaviour of rock mass around openings under a large-scale area. The influence of the connecting angle of two adjacent openings on stress distribution was discussed based on the analytical stress solutions calculated by the Schwarz alternating method. This influence on the failure patterns of surrounding rock mass was further investigated based on PFC-FLAC coupling numerical simulations.

Analytical Stress Solution of Plates
Containing a Single Opening 2.1. Mapping Function. Rock mass around openings in deep mines underground can usually be simplified to be an infinite elastic plate [35,36]. With conformal mapping, the region outside a unit circle in the ζ-plane can be mapped   3 Geofluids onto a simply connected region outside an opening in the z-plane by a mapping function in the following form: where C k = A k + iB k . Both A k and B k are real constants.
With C k truncated to m terms, the mapping function can be written as For m points on the unit circle in the ζ-plane with coordinates of ð1, θ j Þ, the coordinates of their mapping points are ðr j , α j Þ in the z-plane. According to the orthogonality of trigonometric functions, A k and B k can be calculated by [37,38] By sampling 2n points uniformly on the unit circle in the  5 Geofluids ζ-plane, the sampled points in two groups σ e,j ð1, θ e,j Þ and σ o,j ð1, θ o,j Þ can be expressed as By keeping submitting σ e,j and σ o,j into Equation (3) alternately and moving the calculated mapping points into the opening boundary during iterations, the optimal values of A k and B k can be quickly determined once the convergency is reached. For a rectangular opening whose size is  6 Geofluids 40 × 20 (width × height), the optimal mapping function calculated by this method can be expressed as 2.2. Principles of Schwarz Alternating Method. As shown in Figure 1, there is an infinite plate containing two openings, namely, Hole 1 and Hole 2. Mapping functions for them are Z 1 = ωðζ 1 Þ and Z 2 = ωðζ 2 Þ, respectively. Two coordinate systems, namely, the x 1 O 1 y 1 coordinate system and x 2 O 2 y 2 coordinate system, are defined in the plate. The origin of the former one is the centroid of Hole 1, and that of the later one is the centroid of Hole 2. For a point whose coordinate is z 1 in the former coordinate system and z 2 in the later one, the relation between them is z 1 = z 2 + c. The stress solution for the plate is calculated according to an improved Schwarz alternating method [28]. When there is only a hole (Hole 1) in the plate, the boundary stress condition can be expressed as [22] The two complex stress functions φ ð0Þ 1 ðζÞ and ψ ð0Þ 1 ðζÞ for the plate containing Hole 1 can be expressed as [39] B, B′, and C′ are constants relating to the far-field stress condition, which are expressed as where σ ∞ x , σ ∞ y , and τ ∞ xy are far-field vertical stress, horizontal stress, and shear stress, respectively.
The surface force on the suppositional boundary of unexcavated Hole 2 caused by Hole 1 can be calculated by where σ 2 is a point on the boundary of the mapping unit circle in ζ 2 plate for Hole 2 and γ 1 is its corresponding point in ζ 1 plate. The coordinate transformation from σ 2 to γ 1 is realized by s series of steps. Firstly, σ 2 is transformed into t 2 in the Z 2 plate by the mapping function z 2 = ω 2 ðσ 2 Þ. Next, t 2 in the x 2 O 2 y 2 coordinate system into η 1 in the x 1 O 1 y 1 coordinate system via η 1 = t 2 + c. The last step is the transformation of η 1 into γ 1 . Zhang and Lu [27] realized this transformation with the employment of inverse mapping function. Though it is feasible, this method may cause an evitable solution error. Alternatively, the optimization method suggested by Tan et al. [28], which promised little error of the coordinate transformation, is used in this study.

Geofluids
In the first iteration, the boundary stress condition can be expressed as The stress function φ 1 2 ðζ 2 Þ and ψ 1 2 ðζ 2 Þ at the first iteration are the combination of φ 0 2 ðζ 2 Þ and φ ð1Þ 2 ðζ 2 Þ and that of ψ 0 2 ðζ 2 Þ and ψ ð1Þ 2 ðζ 2 Þ, respectively. In the next iteration, the redundant surface force around the boundary of Hole 1 is calculated by the current stress functions. Then, the iterations will be conducted alternately until the desired solution accuracy is achieved.

Calculation Models.
Based on the method introduced above, the stress solutions for plates containing two rectangular holes with different contacting angles are conducted and analysed. The schematic diagram of the two-hole system is shown in Figure 2. σ ∞ x was set as 1 MPa. The hole shape, size, and configuration are the same as those in the DEM simulations.
As shown in Figure 3, under the influence of the adjacent hole, shear stresses at the periphery of Hole 1 were no longer at symmetrical patterns. Such asymmetry is not obvious with the connecting angle of 0°. However, when the holes' connecting line is oblique, the difference between the stress on the side closer to the adjacent hole and on the other side is notable.
Combining Figure 3 with Figure 4 shows that, for monitoring points A and B which are closer to the adjacent hole, notable fluctuations of hoop stress are observed with the increase of the connecting angle. In contrast, hoop stress at points C and D which are further from the adjacent hole is relatively flat, indicating that they are less affected by the connecting angle.
This phenomenon demonstrates that interactive effects between adjacent holes mainly concentrate on the connecting area. With a lateral stress ratio of 0, roof and floor regions of Hole 1 are in intensified tensile stresses; when the lateral stress ratio rises up to the range of 0.5 to 1.0, the intensification level of tensile stress around holes drop dramatically. This illustrates that, under confining conditions, stability of the hole system will be considerably improved with the axial stress ratio being at a reasonable range. Nonetheless, when the axial stress ratio reaches 2.0, tensile stress intensification will reoccur at parts of the hole system, and continuous increasing of the axial stress ratio will impose negative effects on the stability of the hole system.
When λ is 0, tensile stress around Hole 1 is at its roof or floor, which reaches the maximum when the connecting angle is 45°. Accordingly, the stability of the hole system is minimized, and the initial failure is mainly characterized by tensile cracks at the roof of Hole 1 and the floor of Hole 2. Then, with connecting angle increasing, tensile levels gradually decline. When the connecting angle reaches 75°, stress at the roof of the hole converts to be compressive, while that at the floor remains tensile on relatively stable levels, which indicates that initial tensile failures might occur at the floor of Hole 1 or the roof of Hole 2. The discrete element method (DEM) has been used widely for studying the mechanism of rock's deformation and failure characteristics microscopically in the field of rock and soil mechanics, for it can reflect physical phenomena such as rock mass failure and fracture propagation [40,41]. The computing efficiency of DEM, however, compares unfavourably with the finite element method (FEM) or finite difference method (FDM) based on the continuous medium mechanics theory. When the model includes excessive particles, DEM is inappropriate for simulating large-scale rock engineering; the computation may be time-consuming and requires powerful computing capability that PC cannot afford. Whereas its limitation, a discrete-continuum method was adopted [42]. The surrounding rock mass around openings in the desired study area was stimulated by PFC, while other regions by FLAC, which combines virtues such as failure accuracy of DEM and time-saving of FDM.

DEM-FDM Coupling Numerical Modelling
The basic principle of PFC/FLAC coupling is the exchange of data including forces and velocities between elements from FLAC and particles from PFC along the interfaces. The communication between the codes of PFC and FLAC is achieved by FISH socket connection that resembles TCP/IP transmission over the Internet. It allows the exchange of data in FISH arrays in binary on the same machine or two different machines on the same network without the defect of data loss. A dominant approach of PFC/FLAC coupling described by the Itasca consulting group, shown in Figure 5, has been adopted by the majority [43]. In this approach, the data structure that supports  Step ( 10 Geofluids communication is defined as a segment list. Each segment corresponds to the edge of a FLAC zone along the inclusion boundary. A set of controlled particles is identified in PFC2D by finding the particles that intersect the segments, and each controlled particle is associated with only one segment. The time steps in both codes are identical, which is achieved by running FLAC in static mode and PFC2D with differential density scaling. When a computing cycle completes, the forces at the segment end points are applied directly to the corresponding FLAC grid points and the velocities of the FLAC grid points are passed to PFC2D via the corresponding segment end points; the linear interpolation procedure is used to map these velocities to the PFC2D controlled particles. The approach above requires neither special modelling skills nor small-sized mesh around the boundary, making PFC/FLAC coupling convenient as well as reducing the number of meshes in the FLAC model. The basic preconditions for this approach, however, are that data along the boundary should be linear distribution rather than erratic fluctuation. It is achievable only if the boundary is far enough from the source of stress disturbance (such as goaf underground), which means that the size of the PFC model should be large enough to realize PFC/FLAC coupling with acceptable coupling precision. In view of this, as shown in Figure 6, a refined approach that transfers the data exchange mode from "segment-particle" to "grid-particle" was proposed. Regular particles were created one by one along the boundary, which accurately matched to certain grids on the other side. Instead of being calculated by the linear interpolation algorithm, the velocities set in each particle were derived from unique corresponding grids during the data exchange process. Though this approach requires higher mesh density in FLAC as well as complex PFC modelling approaches, the increased meshes in FLAC have limited effects on computing speed since FDM, showing an effi-ciency advantage over DEM, was applied by FLAC. In the actual test, the particle number in PFC is responsible for much of the long computing time of FDM-DEM coupling simulation while the mesh number in FLAC for little. This approach reaches the maximum coupling precision under the condition that the particle number is determined, with computing time considered at the same time. Figure 7. The distance between two-hole centroids d = 20 m. The connecting angle α is defined as the angle of the line segment between centroids of two holes with respect to the horizontal direction. Figure 7 illustrates the layout of two rectangular openings in numerical models. Combined goafs in model 1 were on the same horizontal plane. The connecting angle α is the angle of the line segment between centroids of two holes with respect to the horizontal direction. In this study, connecting angles from 0°to 90°with an interval of 15°are considered. The size of the rectangular opening and of the DEM models is 40 m × 20 m and 120 m × 120 ðlength × widthÞ, respectively. The shortest distance d between two openings is 1 m. The mesoparameters of numerical specimens are listed in Table 1. 3.3. Stress and AE Characteristics. Acoustic emission (AE) activity is closely related to evolution and propagation of rock defects. Recording AE numbers of goaf models during the failure process under compression can provide significant information about the goaf's instability behaviour, which is conducive to revealing the failure mechanism of rock mass around openings. In numerical simulation, AE events originating from material failure can be accessed by tracking bond breakage between particles. Particles are associated with each other by a defined bond between them. When the bond strength of contact is overcome, strain 11 Geofluids energy stored in the particle contacts will be released as kinetic energy in the form of a seismic wave that could be an indication of an AE activity, based on which, the simulation of AE can be reached [44]. This method has been widely adopted in the field of rock mass engineering. Although the particle size and number are limited by the computing capability of computers, AE characteristic in rock revealed by PFC numerical simulation is consistent with that revealed by laboratory tests, which is very significant to the study of AE behaviours in rock under different conditions [45].

Numerical Models. The layout of two rectangular openings is illustrated in
The stress condition at the DEM model boundaries is relatively stable, which cannot reflect the stress change in the rock mass around the openings. The stress of a monitoring point, which is 60 m right above and away from the middle point of the line segment between the centroids of the two openings, was recorded. The stress curves as well as acoustic emission numbers obtained in the simulations on different models are shown in Figure 8. At the stages of initial loading, each model is mostly in linear elastic deformation, with little AE events observed. Then, AE events become more and more frequent with slight stress curve fluctuation in the prepeak stage. Some typical points with significant AE burst and stress change are marked in each curve. It can be seen that there is an obvious correlation between AE event and stress evolution. Almost every AE burst happens with corresponding stress fluctuation, which may indicate sudden appearance of failure within rock mass. Large numbers of AE events emerged in the postpeak stage with intensive failure, but because of the confinement of the far-field rock mass, postpeak stress remains stable to the end of calculation.
3.4. Failure Patterns. The ultimate failure patterns around openings with three representative connecting angles are plotted in Figure 9. When two openings exist at the same horizontal level (α = 0°), cracks initiate at the top of the connecting area between them. Then, intensive failure appears in the connecting area, indicating that the connecting area lost its bearing capacity. Tensile cracks initiate at the top and the bottom of the damaged connecting rock mass, respectively, and propagate along the loading direction. At the same time, failure is observed in the corners of the openings. Later, tensile cracks from the damaged pillar continue to develop; both side walls of the goaf are squeezed to deform and spalled subsequently, resulting in sloughing at last. Lots of failure appears around openings and causes the instability of the surrounding rock mass. When α = 45°, failure only exists in a narrow area between the floor of the upper goaf and the roof of the lower goaf; adjacent side walls between two openings (the left side wall of the upper goaf and the right of the lower) both remain almost intact. When θ = 90°, failure forms at the side walls of the lower opening and spreads around them. Inclined cracks from the failure zone are observed. All the failure zones, however, mainly concentrate around the lower opening, contrary to what happened to the upper goaf, the region around which remains almost intact with only little failure near the opening's corners.
The analytical stress solutions around Hole 1 for the three representative kinds of two-opening systems under uniaxial compressive stress conditions, whose failure patterns are given in Figure 9, are plotted in Figure 10. It can be seen that compressive stress concentration appears at opening corners and tensile stress is observed in the roof and floor of the openings in all cases, but stress distribution and concentration level vary from case to case. When the connecting angle is 0°or 90°, two openings are symmetrical to each other about X or Y axis. Accordingly, the stress in the two-opening system also shows symmetrical distribution. When the connecting angle is 45°, the stress concentration at the upper left corner of Hole 1 close to the adjacent hole is much stronger than that in the two-opening system with other connecting angles. Combining Figures 9 and 10 shows that when α = 0°, the highest compressive stress concentration appears in the connecting area, which leads to the coalescence between two openings. In contrast, when α = 90°, the stress around the connecting area is relatively lower than that in other areas and no failure is found in this area in numerical simulation.

Discussion
This study emphasizes the stress distribution and failure patterns of rock mass surrounding openings within a large-scale area without model boundary effect. As the numerical model size is much greater than the opening size, the failure patterns of rock mass around openings in numerical models are in line with the stress concentration areas determined by analytical stress solutions. Though rock mass is a kind of elastic-plastic material, for hard rock with strong brittleness, it is mainly subjected to elastic deformation in the prepeak stage. Therefore, the analytical elastic solutions for plates containing openings can be an important reference of failure predication within the surrounding rock mass.
There have been a great number of studies on the mechanical behaviour of rock mass containing openings or flaws based on laboratory experiments. However, mechanical tests for rock specimens may lead to an evitable boundary effect, which may significantly interfere with the stress distribution and fracturing evolution of rock mass around openings. Many laboratory tests show that the rock specimens are finally broken by tensile or shear cracks connecting opening corners and specimen boundaries. Figure 11 presents the failure patterns of rock specimens containing two rectangular openings under uniaxial compression tests conducted by Zhou et al. [29]. Comparing Figure 11 with Figure 9 shows that the coalescence patterns between openings in rock specimens agree well with the corresponding large-scale numerical models. However, for all rock specimens, the instability is dominated by macro shear or tensile cracks from openings and through the whole specimen, which are suppressed in the large-scale numerical models.

Geofluids
For the one hand, strong interaction between the rock specimen boundaries and openings may promote crack propagation from openings to boundaries. For another, rock mass in the large-scale area may also contribute to confining the development of cracks. Therefore, the laboratory tests may exaggerate the influence of deep-buried openings on rock mass stability. The large-scale numerical modelling via the PFC-FLAC coupling method is expected to be an important addition to the study of the mechanical behaviour of deepburied surrounding rock mass.

Conclusions
In this study, the Schwarz alternating method based on complex variable theory was used to calculate the stress solution for the two-opening system. In addition, PFC-FLAC coupling numerical modelling was employed to study the mechanical behaviour of rock mass containing two rectangular openings with different connecting angles underground. With the combination of numerical and analytical solutions, the effect of connecting angles on the stress distribution and failure characteristics of rock mass around openings was discussed. The main conclusions of this paper include the following: (1) According to the analytical stress solutions for elastic plates containing two rectangular openings, the existence of confining stress within a certain range contributes to relief of the tensile stress concentration around the openings, which is supposed to help improve the stability of rock mass around openings. With the increase of lateral pressure coefficient, the stress field changes gradually and tensile stress concentration may appear again in new areas once it reaches a certain level (2) The opening layout mainly affects the stress distribution in the connecting area between two adjacent openings, which further dominates the coalescence pattern and stability of the two-opening system. Under uniaxial compression, the connecting area is a high compressive stress concentration area in which appears coalescence failure during numerical tests when the two openings are arranged horizontally. When the connecting angle is 90°, the connecting area is a lowstress area and remains intact during the numerical test (3) The AE characteristics in the numerical models show that massive failure mainly happens in rock mass in the postpeak stage, which follows the coalescence between the two openings. In the prepeak stage, rock 14 Geofluids mass around openings is mainly subjected to elastic deformation. The positions of initial failure always agree well with the stress concentration areas determined by analytical solutions

Data Availability
The figures and tables used to support the findings of this study are included in the article.