Investigation of cracking behavior and mechanism of sandstone specimens with a hole under compression

International Journal of Mechanical Sciences 163 (2019) 105084

Contents lists available at ScienceDirect

International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci

Investigation of cracking behavior and mechanism of sandstone specimens with a hole under compression Hao Wu a,b, Guoyan Zhao a,∗, Weizhang Liang a a b

School of Resources and Safety Engineering, Central South University, Changsha 410083, PR China Department of Mining and Geological Engineering, The University of Arizona, Tucson 85721, USA

a r t i c l e

i n f o

Keywords: Crack Hole Digital image correlation (DIC) Acoustic emission (AE) Complex variable method

a b s t r a c t In this paper, the cracking behaviour of rocks with a hole loaded in compression was investigated both experimentally and theoretically. First, a series of uniaxial compression tests on prismatic red sandstone specimens containing a circular or a rectangular hole were performed using a rock mechanics testing system in conjunction with digital image correlation (DIC) and acoustic emission (AE) systems. Afterwards, the theoretical solutions of stress for a circular and a rectangular hole under different stress conditions were solved by combining the complex variable method and the Boxʹs complex algorithm. Experimental results show that the mechanical properties of rock specimens are significantly degraded by the hole and influenced by its shape. Crack evolution proceeds from primary tensile cracks through secondary tensile cracks and slabbing fractures to shear cracks, which can be visually displayed by the real-time deformation fields. The failure of specimens containing a hole under uniaxial loading is induced by the coalescence of the shear cracks and slabbing fracturing zones. Besides, several critical stresses associated with cracking behaviour were estimated based on a proposed method and the AE measurement. Theoretical investigation indicates that both hole shape and lateral stress ratio are important factors influencing the stress distribution, and the initiation mechanisms of different cracks can be well revealed according to the stress states around the holes under uniaxial compressive loads.

1. Introduction Cracks and pores (or holes) are two common types of defects in rock masses. Under external loads, stress concentrations are most likely to occur at tips and corners of these defects. This may result in the initiation and propagation of cracks. Currently, it is commonly recognized that the defects significantly weaken the mechanical properties of rocks and lead to more complex failure characteristics. To deeply understand the deformation and failure mechanisms of rock masses, it is essential to investigate the effect of defects on cracking behaviour of rock materials. In the past decades, extensive experimental and numerical investigations have been conducted on the mechanical responses and failure characteristics of rock or rock-like specimens with pre-fabricated cracks. Generally, for specimens containing a single crack under uniaxial compression or low confining stresses, wing cracks grow firstly from the crack tips along a curvilinear path with the increasing load, and eventually cease when they turn parallel to the maximum compression direction [1–4]. Next, secondary cracks, whose initiation directions are material dependent, appear later in most cases and are responsible for the macro failure [5,6]. Under the condition of high confining stresses, only secondary shear cracks occur. To reveal the propagation mechanisms

∗

of cracks initiated from the flaw tips, Wong and Einstein [7] conducted uniaxial compression tests on molded gypsum and Carrara Marble with a pre-existing crack, and classified the appeared cracks into seven types. Yang and Jing [8] also carried out similar tests on sandstone specimens, and proposed nine types of cracks. Ashby and Hallam [9] reported that wing cracks grow most easily in finite brittle specimen when the inclination of the embedded single crack is about 45°, and the failure is induced by axial splitting. Nemat-Nasser and Horii [10] claimed that the confining pressure is the most prominent factor in determining the failure modes. For specimens containing two cracks under uniaxial loading, research showed that crack coalescence depends not only on material properties but also on crack geometry and stress conditions [11]. For instance, Wong and Chau [12] observed three patterns of crack coalescence in sand-like specimens containing two parallel cracks with different crack inclinations, bridge angles and frictional coefficients, namely wing tensile mode, shear mode and mixed tensile-shear mode. Wong and Einstein [13] pointed out that the crack coalescence mode of gypsum specimens with two parallel cracks is dominated by the crack inclination, bridge angle and ligament length. Lee and Jeon [14] analyzed the effect of materials on crack coalescence in specimens containing a

Corresponding author. E-mail address: [email protected] (G. Zhao).

https://doi.org/10.1016/j.ijmecsci.2019.105084 Received 13 June 2019; Received in revised form 9 August 2019; Accepted 15 August 2019 Available online 17 August 2019 0020-7403/© 2019 Elsevier Ltd. All rights reserved.

H. Wu, G. Zhao and W. Liang

horizontal crack and an inclined crack. It was found that only tensile cracks are formed in PMMA specimens, while both tensile cracks and shear cracks are observed in gypsum specimens and granite specimens. Besides, Wang et al. [15] developed a novel conjugated bond-based peridynamic model to simulate the crack propagation and coalescence behaviour of rock specimens containing two prefabricated flaws. For multiple cracks, the cracking behaviour is relatively complicated, which is affected by many factors such as the inclination, length, number and distribution of cracks, bridge angle and ligament length, multi-crack intersection and material properties as well as loading conditions [16–20]. Overall, the formation of different failure modes is a result of competition of tensile cracks, shear cracks and mixed cracks. Furthermore, some scholars also explored the effects of 3-D cracks and loading methods on mechanical behavior of jointed specimens, which improves the understanding of the failure mechanisms of rock masses [21–27]. In addition to crack-like flaws, natural hole-like defects also have a pronounced influence on the mechanical behaviour of rock masses. In fact, underground excavations can be regarded as artificial hole flaws in rocks. In deep rock engineering, deformations and brittle failure of these openings increase because of the complex geological environment characteristics such as high stresses, high temperatures, high water pressures and strong mining disturbance [28–30]. As the mining depth increases, some special failure patterns occur frequently, e.g., rock bursts, spalling and zonal disintegration [31–34]. This poses huge challenges to the safety of excavations. Therefore, more and more attention has been paid to the induction mechanisms of these unconventional failure phenomena. Until now, some attempts have been made to study the failure behaviour of rock or rock-like materials embedded hole flaws subjected to compression. Hoek [35] performed a photoelastic experimental study on Columbia Resin 39 and glass specimens containing a circular hole in a biaxial stress field, and first identified the tensile cracks and remote cracks around the hole. These cracks have also been observed in rock materials under similar experimental conditions [36,37]. Based on physical and numerical models, it is concluded that three types of cracks are formed in the vicinity of the circular hole under uniaxial compression or low confining stresses; that is, tensile cracks on the roof and floor, slabbing fractures and remote cracks away from the hole [38]. The propagation directions of these cracks are parallel to the maximum compressive stress [39]. Dzik and Lajtal [40] found the tensile cracks propagate in a stable way for small cavities (radius < 20 mm) and turn to unstable propagation for large cavities (radius > 40 mm). By using a self-developed code, Tang et al. [41] modelled the effect of specimen width and hole diameter on splitting failure of brittle specimens containing a circular hole. Results indicated that cracks grow more easily in narrow specimens with a large hole than in wide ones with a small hole. To obtain the initiation stresses of the above three types of cracks, Carter et al. placed 19 strain gauges around a circular hole in granite specimens under uniaxial loading, and found that the critical stresses are highly dependent on the hole diameter [38,42]. However, the locations of the strain gauges are difficult to design due to the unknown propagation paths of the cracks in advance. Additionally, researchers,

International Journal of Mechanical Sciences 163 (2019) 105084

e.g., Mastin [43], Carter [44] and Martin [45,46], examined the effect of hole size on the required stress causing slabbing failure, and fitted a formula for the depth of the slabbing fracturing zone. Based on the acoustic emission (AE) technique, Fakhimi et al. [47] detected the locations of microcracks in Berea sandstone with a circular cavity under biaxial loading. Besides, the impacts of factors such as the heterogeneity, temperature, lateral stress ratio, dynamic load, tensile load, number and distribution of holes on the failure mechanisms were also investigated both experimentally and numerically [48–56]. Moreover, Li et al. [57], Yin et al. [58] and Fan [59] further explored the mechanical properties and cracking behaviour of specimens containing crack-hole combined flaws under loads. Compared to the number of investigations on circular cavities, there is relatively little information concerning crack evolution around noncircular holes in specimens [60]. In practical engineering, rectangular and other noncircular cross sections are widely used for excavations like roadways, tunnels and drifts. Therefore, it addresses a need to systematically investigate the cracking behaviour of non-circular hole-like flaws. Essentially, the initiation, propagation and coalescence of cracks are driven by stress. Thus, figuring out the stress distribution around the cavity is meaningful to reveal the crack initiation mechanism. But up to now, studies are still limited regarding the analytical stresses around complex shaped holes. Also, it should be emphasized that it is hard to accurately obtain the crack evolution process by an ordinary digital camera, because microcracks on digital photos are hardly recognized by naked eyes [61]. Hence, in this study, the digital image correlation (DIC) and AE techniques were jointly applied to investigate the cracking behaviour and damage characteristics of sandstone specimens containing a circular or a rectangular hole under uniaxial compression. After that, the complex variable theory coupled with an optimization algorithm was employed to find the asymptotic solutions for stresses around the holes under different lateral stresses, and then the cracking mechanisms were explained. 2. Experimental procedure 2.1. Rock material and specimen preparation In this study, the representative sandstone taken from Linyi city of China was selected as the experimental material. The sandstone is reddish brown in color. Its mineral composition and texture were examined by a petrographic microscope under the plane polarization light (PPL) and cross polarization light (XPL), which are shown in Fig. 1 [32]. Results indicated that the minerals in this sandstone mainly include K-feldspar (5%), zeolite (8%), calcite (9%), plagioclase (35%), quartz (42%), and other opaque minerals (1%). The rock can be classified as tuffaceous feldspar quartz sandstone with a fine-to-medium-grained sand-like texture and a massive structure. The mechanical parameter values of the dry density (𝛾), P-wave velocity (Vp ), uniaxial compressive strength (𝜎 d ), tensile strength (𝜎 t ), Young’s modulus (E), Poisson’s ratio (𝜇), cohesion (Cc ), friction angle (𝜙) and the fracture toughness (KIC ) are listed in Table 1. Fig. 1. Micrographs of red sandstone under: (a) PPL and (b) XPL. (Qtz-quartz; Pl-plagioclase; Cal-calcite).

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Table 1 Mechanical parameter values of red sandstone. 1

𝛾 (kg/m3 )

Vp (m/s)

𝜎 d (MPa)

𝜎 t (MPa)

E (GPa)

𝜇

Cc (MPa)

𝜙 (°)

KIC (MPa•m 2 )

2472.2

3174.5

99.3

5.3

24.4

0.26

19.0

40.4

0.6

Fig. 2. Schematic of fabricated sandstone specimens (all the dimensions are given in mm): (a) specimen IS; (b) specimen CS and (c) specimen RS.

Fig. 3. Schematic diagram of testing apparatus.

All specimens prepared for uniaxial compressive tests were extracted from one relatively homogeneous rock block, and were machined into prismatic shape with dimensions of 100 mm × 25 mm × 150 mm (length × width × height). Three groups of specimens were prepared, i.e., intact specimens (IS), specimens containing a circle hole (CS) and specimens containing a rectangular hole (RS), as shown in Fig. 2. Each group has three specimens. The holes in the center of the specimens were fabricated by high-pressure water jet cutting technology, and the cross-sectional area of the circular hole is basically equal to that of the rectangular hole. According to the standards proposed by the International Society for Rock Mechanics (ISRM), the non-flatness and non-perpendicularity of the grinding specimen surfaces were controlled within 0.02 mm and 0.001 radians respectively [62]. To minimize the friction effect, the loading ends of specimens were lubricated using vaseline before testing. 2.2. Experimental apparatus As shown in Fig. 3, the testing setup is mainly composed of three parts: the Instron 1346 rock mechanics testing system, the DIC system, and the AE system. All systems need to be started synchronously to obtain the relation between monitoring parameters. 2.2.1. Loading system Uniaxial compressive tests were carried out on an Instron 1346 rock mechanics testing system, which is able to conduct compression tests, tension tests, bending tests and shearing tests. The maximum load and

accuracy of this system are 2000 kN and ±0.5%, respectively. According to our previous research, we found that it is better to apply the load on specimens by the displacement-controlled loading mode with a constant displacement rate of 0.6 mm/min [63]. The axial displacements of specimens were monitored simultaneously by an extra linear variable differential transformer. 2.2.2. DIC system As a powerful optical technique, DIC is popularly used to measure the full-filed deformation in solid mechanics. Compared with the conventional deformation measurement techniques such as extensometer and strain gauges, this technique has significant advantages, e.g., high reliability, wide measuring range, non-contact, real-time, full-field, and noninterference [64,65]. Based on the principle that the DIC works when the relative value of the gray-scale stays the same, the movements of pixels can be tracked by matching the same pixels in two images before and after deformation. Thus, the strain and displacement fields can be obtained through image correlation and numerical computing. To capture the pixels reliably, it is necessary to make an artificial speckle field on specimen surfaces to be observed. First of all, these surfaces were cleaned by a degreaser and then were sprayed with white paint to form a thin uniform background. After the basecoat is dried, black paint was then strategically sprayed on it by a professional to generate randomly distributed speckles. Requirements of the speckle pattern can be found in reference [61]. Based on the sequential speckle images taken during tests, GOM correlate software was used to implement the deformation calculation, steps are as follows (Fig. 4): (1) import the continuously recorded speckle images into the GOM software;

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Fig. 4. Flow chart of deformation calculation by DIC technique. Table 2 Geometric sizes and mechanical properties of sandstone specimens under uniaxial compression. Specimen no.

Length (mm)

Width (mm)

Height (mm)

Density (kg m−3 )

𝜎 d (MPa)

𝜀p (‰)

E (GPa)

eb (kJ m−3 )

ea (kJ m−3 )

k

IS-1 IS-2 IS-3 CS-1 CS-2 CS-3 RS-1 RS-2 RS-3

100.1 99.6 100.7 99.9 100.0 100.0 100.2 100.2 100.2

25.0 24.9 24.8 25.1 25.0 25.9 24.8 25.9 24.9

150.1 150.5 150.7 150.4 149.3 150.2 149.8 150.2 150.6

2391.4 2427.5 2385.5 2417.9 2373.6 2362.4 2369.6 2415.7 2385.5

100.9 105.6 101.3 73.7 72.1 70.0 59.3 64.5 61.1

5.57 5.96 6.97 4.58 4.32 4.73 4.75 4.89 4.88

21.6 21.6 19.1 17.2 17.2 17.5 15.4 15.8 15.4

100.28 111.42 114.90 181.42 162.17 160.80 141.42 151.02 110.88

55.38 47.78 −88.33 68.98 66.64 50.48 19.23 40.27 27.09

1.51 1.62 / 0.63 0.66 0.66 0.71 0.69 0.88

(2) select a deformed image and the image before deformation as the target and reference respectively; (3) choose a region of interest (ROI) in the speckle images and create the surface component; (4) divide this region into many subsets with a certain size and record the coordinates of their centroids in the reference; (5) search for the matched subsets in the target image using correlation algorithm, and determine their centroidsʹ positions; and (6) calculate the displacement and strain through digital image processing. The DIC system consists of a white light source, a charged couple device (CCD) camera and an image collection system. The resolution and maximum sampling frequency of the CCD camera (Basler/piA2400-17 gm) are 2456 × 2058 pixels and 17 FPS, respectively. This camera under the illumination of a white fill light was placed perpendicular to the speckle surfaces of specimens at a distance of 1 m. To cover the whole speckle surface of the specimen, with the aid of an image acquisition software (Pylon Viewer), the width and height of the captured image were set to 1100 pixels and 1500 pixels, respectively. During loading, the images were continuously recorded at a rate of 15 frames/s.

2.2.3. AE system Acoustic emission refers to the phenomenon that acoustic wave radiation occurs when rock materials undergo irreversible changes in their internal structures, such as crack initiation and propagation. Therefore, AE is usually used to detect, locate and characterize damage during rock failure. In this study, a PCI-2 AE instrument produced by Physical Acoustics Corporation was applied to monitor the internal damage and cracking behaviour of specimens under uniaxial compression. Two pico-type AE sensors with a resonant frequency of 250 kHz were attached to the back surface of the specimen by adhesive tape, and vaseline was used as a coupling agent. It is worth noting that the one with satisfactory AE signals were selected for data analysis. To filter the noise generated by the surrounding environment, the detection threshold was set at 45 dB. The

continuous monitored AE signals amplified by two 2/4/6 preamplifiers (40 dB) were recorded and processed using AE-win software.

3. Experimental results 3.1. Strength and deformation characteristics Table 2 shows the geometrical sizes and test results of all specimens subjected to uniaxial loading. The average uniaxial compressive strength of group IS is 102.61 MPa, whilst that of groups CS and RS is 71.93 MPa and 61.63 MPa, respectively, which are 29.91% and 39.94% lower than that of the intact specimens. In respect of the Youngʹs modulus, group IS has the largest value (20.77 GPa), followed by groups CS (17.30 GPa) and RS (15.53 GPa) with reduction rates of 16.69% and 25.20%, respectively. Compared with the intact specimens, the peak strains of the specimens containing a hole drop by more than 20%. The values can be ordered from largest to smallest as: groups IS (6.17‰) ˃ RS (4.84‰) ˃ CS (4.54‰). It is concluded that the mechanical properties of specimens are significantly weakened by the prefabricated hole and the degradation amplitude is closely associated with its shape. Fig. 5 shows the axial stress-strain curves of three representative specimens under uniaxial compression. According to the characteristics of these curves, we divided the deformation process into five stages, i.e., the initial pores and fissures compaction stage (OA), elastic deformation stage (AB), stable deformation stage (BC), unstable deformation stage (CD) and post-peak stage (DE). In the initial stage (OA), the natural flaws in the specimens close rapidly under loads, which results in the downward concave curves. At stages AB and BC, the strain increases linearly with the axial stress. From point C, the curves turn into upward concaves because of the damage inside the specimens. At the last stage, the stresses fall to zero drastically, accompanied by violent rock ejections and crisp sounds. This suggests that the brittleness of this kind of rock is significant. Hereby an energy-based brittleness index proposed

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

the number of the yellow dots becomes larger and the average strain is 0.90‰. With the increase of the stress to 49%𝜎 d , the yellow dots gradually gather near the leading diagonal and right border of the specimen, but no macro cracks appear, and the average strain is 0.80‰. The reason is that the density and range of microcracks are not large enough. When the stress reaches 94%𝜎 d , a splitting tensile crack 1a parallel to the loading direction is observed from the lower-right corner of the specimen. At this time, the average strain increases to 1.52‰. At the peak stress, another splitting tensile crack 1b propagating along the loading direction appears. Meanwhile, a yellow shear band is formed along the diagonal and then grows into a shear crack 3. At this moment, the specimen remains intact with an average strain value of 1.47‰. At a stress of 79%𝜎 d after the peak, cracks 1a and 1b get connected, and the specimen has an average strain value of 2.01‰. Once they coalesce with the shear crack 3, the macro failure will occur.

Fig. 5. Stress-strain curves of representative specimens under uniaxial compression.

by Munoz et al. [66] was used to quantitatively characterize the brittle behaviour, which can be expressed as: 𝑒e 𝑘= (1) 𝑒𝑎 + 𝑒𝑏 Where ee =𝜎 d ˆ2/2E denotes the elastic energy density at the peak point, while eb and ea mean the energy densities before and after the peak stress, respectively. Note that the higher the value of k, the more significant the brittleness of the rock is. The values of k for different specimens are given in Table 2. For group IS, the average value is 1.57, while the values of groups CS and RS are 0.65 and 0.76, respectively. The results are agreeable with the slopes of the stress-strain curves at the post stage. Clearly, the brittleness of the specimens with a circular hole is the least significant among these specimens, leading to a relatively high deformation-bearing capacity. The invalid k value for the specimen IS-3 is attributed to the occurrence of the abnormal unloading. Besides, it is further observed in Fig. 5 that, for specimens containing a hole, stress drop occurs at stage CD, which is induced by the sudden appearance of cracks. 3.2. Crack evolution and failure patterns By using the DIC technique, the full-filed deformations of specimens can be derived based on the recorded speckle images during the tests. For each specimen, six deformed images corresponding to the five deformation stages and the peak stress were selected for fracture analysis. Literature indicates that cracks can be well reflected by strain localization zones, that is to say, strain concentration regimes represent potential macro-cracks [39]. Fig. 6 shows the major principal strain distributions of three representative specimens (IS-3, CS-3 and RS-1) at different deformation stages. Note that the numbers in the figures denote the patterns of cracks, and the superscript lowercase letters mean the emergence sequence of the same type of cracks. Besides, based on the crack orientations, macro cracks are divided into three categories: tensile cracks whose propagation directions are parallel to the maximum principal stress, shear cracks propagating obliquely with respect to the maximum principal stress, and tensile-shear mixed cracks that show both characteristics. 3.2.1. Specimen IS-3 As shown in Fig. 6a, the crack evolution of the specimen IS-3 at different stress levels is clearly presented. At a stress of 10%𝜎 d in the first stage, there are some randomly distributed yellow dots with an average strain value of 0.43‰. This results from the closure of the natural fissures and pores inside the specimen. When the stress approaches 30%𝜎 d ,

3.2.2. Specimen CS-3 Fig. 6b shows the cracking process of the specimen CS-3 under uniaxial compression. At low stresses, only the yellow dots as in specimen IS-3 are formed, and the number and average strain increase with increasing stress. For example, the strain increases from 0.28‰ at stress 7%𝜎 d to 0.30‰ at stress 26%𝜎 d . At stress 57%𝜎 d , two clear primary tensile cracks 1a and 1b on the roof and floor of the cavity are observed, with an average strain value of 0.47‰. Both of them propagate parallel to the loading direction in a stable way with the increase of loads. When the stress 89%𝜎 d is reached, it is observed that two vertical secondary tensile cracks 1c and 1d appear on the left side of the hole. They propagate rapidly towards the hole sidewalls along the loading direction. Besides, slabbing fractures emerge at the left sidewall at the same time, accompanied by many pieces of rocks being peeled off. At this time, the average strain is 0.69‰. With the increase of the axial stress to the peak stress, the secondary tensile cracks, such as 1c and 1e located at the upper left and lower right corners of the circular hole, have intersected with the slabbing fracturing zones, and the average strain increases to 0.97‰. Afterwards, two shear bands gradually emerge on the leading diagonal and eventually evolve into shear cracks 3a and 3b . Subsequently, the macro failure occurs abruptly when the shear cracks intersect with the hole. 3.2.3. Specimen RS-1 The major principal strain contours of the specimen RS-1 at different stress levels are presented in Fig. 6c. It can be calculated that the average values of the major principal strains at the five stress levels shown in Fig. 6c are 0.67‰, 1.12‰, 1.37‰, 1.19‰, 0.77‰ and 1.08‰. The progressive failure process is almost the same as in specimen CS-3, that is, cracks proceed from primary tensile cracks through slabbing fractures and secondary tensile cracks to shear cracks during loading. The essential difference between them is the cracking sequence of the secondary tensile cracks and the slabbing fractures. Since the brittle failure of specimens is sharp, it is difficult to capture enough images to characterize the complicated post-peak mechanical behavior. Furthermore, the horizontal displacement distributions of the three specimens at the same stress levels were also obtained by the DIC technique. As can be seen in Fig. 7, the variation of the displacement agrees well with the crack evolution. Taking specimen RS-1 as an example, at the beginning of loading, the deformation looks randomly distributed in the specimen. This is caused by the closure of the micro-defects. When the stress reaches 26%𝜎 d , a number of small deformation areas concentrate at two sides of the hole. At stress 57%𝜎 d , two local deformation zones with large displacements separated by two vertical lines (cracks 1a and 1b ) on the roof and floor of the cavity appear on the two sides. For the left side, the horizontal displacement is negative, while it is positive for the right one. With the stress increases to 89%𝜎 d , a vertical gradient line representing the secondary tensile crack 1c is formed at the upperright part of the hole. When the peak stress is approached, another clear

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Fig. 6. Major principal strain contours of specimens at different stress levels: (a) specimen IS-3; (b) specimen CS-3 and (c) specimen RS-1. (Notes: the strain is treated as positive when tensile and negative when compressive).

Fig. 7. Horizontal displacement distributions of specimens at different stress levels: (a) specimen IS-3; (b) specimen CS-3 and (c) specimen RS-1.

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Table 3 Several critical stresses of three representative specimens subjected to uniaxial compression.

Fig. 8. Failure patterns of specimens under uniaxial loading: (a) specimen IS-3; (b) specimen CS-3 and (c) specimen RS-1.

vertical gradient line at the lower-left part of the hole is observed, indicating the formation of the secondary tensile crack 1d . At a stress of 81%𝜎 d at the post-peak stage, two distinct diagonal lines representing the potential shear cracks 3a and 3b appear, and the displacement field is finally separated by these lines with a large displacement gradient until the macro failure occurs. From the above analysis, it can be concluded that the failure of specimens under uniaxial loading is a progressive process of crack nucleation, propagation and coalescence with each other. The failure patterns of the above three specimens are shown in Fig. 8. For the intact specimen, the ultimate failure is attributed to the coalescence of the splitting tensile cracks and shear cracks. Thus, it can be treated as tensile-shear mixing failure. In contrast, the failure of specimens containing a hole is mainly caused by the coalescence of the shear cracks and the slabbing fracturing zones. It is found that the primary and secondary tensile cracks are not responsible for the macro failure. Due to the unloading and end friction effects, the secondary tensile cracks may propagate further at the end of experiments, but their direction deflects when approaching the loading ends. Besides, few secondary shear cracks such as cracks 3c and 3d in Fig. 8c are also formed during failure. Some small rock debris drop from the speckle surfaces because of the high compressive stress. This leads to the formation of surface cracks such as the curving crack in the upper part of the intact specimen. The surface cracks have no effect on the failure mode since they do not run through the thickness of the specimens. 3.3. AE reaction and stress thresholds Fig. 9 shows the curves of the axial stress, AE counts and cumulative AE counts versus the time of the above three specimens. The AE characteristics at different deformation stages shown in Fig. 5 can be summarized as follows. At the first stage, few AE events are detected because the closure of the natural defects could not result in any or much energy release. In this stage, the AE acitivities of the intact specimen are more active than that of the specimens with a hole. When the elastic deformation stage (AB) is approached, the AE activities are quiet because there are only a few micro cracks. As a result, both the AE counts and accumulated AE counts increase linearly and slowly with the increase of the axial stress. Meanwhile, the strain energy is gradually stored in the specimen and can be calculated according to Ue = ee V (where Ue and V denote the strain energy and the volume of the specimen, respectively). At stage BC, compared with the previous two stages, the AE acitivies of specimens containing a hole are relatively active, and the accumulated AE counts increase exponentially. This suggests that the primary tensile cracks start to initiate and propage in a stable way. For the intact specimen, the AE counts nearly keep constant and the cumulative AE counts increase linearly, but no macro cracks are formed. The reason is that the density and range of the appeared microcracks are not large enough. From point C, it can be seen that evident AE events occur because of

Specimen no.

𝜎 a (MPa)

𝜎 b (MPa)

𝜎 c (MPa)

𝜎 d (MPa)

IS-3 CS-3 RS-1

24.2 17.3 15.8

44.8 19.4 22.7

92.2 59.6 52.6

101.31 70.01 59.30

the sudden appearance of cracks such as the splitting tensile cracks in the intact spcimens and the secondary tensile cracks in the pre-holed specimens. At the post-peak stage, the AE counts and cumulative AE counts increase abruptly in all specimens, indicating the coalescece of shear cracks with other cracks into macro failure. Moreover, it is found that the AE activities of the intact specimen are the most active, while the number of the AE events of the samples with a rectangular hole in each stage is the smallest. This is because the focusing of fracture development near the hole reduce the AE output by reducing the amount of diffuse fracture development throughout the intact rock. In a word, the variation of AE signals is in good agreement with the evolution of cracks. When the applied stresses on specimens reach certain levels, cracks will initiate, propagate and coalescence with neighbouring cracks. In the present study, several stress thresholds related to cracking behavior were defined as the crack closure stress (𝜎 a ), crack initiation stress (𝜎 b ) and crack damage stress (𝜎 c ), corresponding to the stresses at points A, B and C in Figs. 5 and 9, respectively. Based on the AE signals, the values of the latter two stresses can be easily determined, namely the stresses when the AE counts increase significantly. Considering the high human error in the traditional method for identifying the crack closure stress, we proposed a new method to determine its approximate value. Taking the specimen RS-1 as an example, the calculation procedures are as follows: (1) select two points with stresses of 20%𝜎 d and 30%𝜎 d (i.e., 11.9 and 17.8 MPa) on the stress-strain curve, marked as M and N respectively. (2) determine the slope of the line between the two points to be 15.43, and then plot a line 𝜎 = 15.43𝜀 from the origin point O (𝜎 and 𝜀 denote the axial stress and strain, respectively, see Fig. 10a). (3) calculate the axial strain difference between the plotted line and the stressstrain curve at the same axial stress. (4) plot the axial strain difference versus axial stress curve, and find the point where the strain difference value is the largest, denoted as A. The stress at point A (15.8 MPa) is the crack closure stress, as presented in Fig. 10b. Table 3 shows the results of the several critical stresses of the three specimens. The crack closure stresses of the specimens containing a hole are lower than that of the intact specimen. In respect of the crack initiation stress, the value of specimen IS-3 is 44% 𝜎 d . This is consistent with the results in literature [67,68]. In contrast, the values of specimens CS3 and RS-1 are 28% and 38% of their peak stresses, respectively. For the crack damage stress, the value of the specimen IS-3 is 91%𝜎 d , while that of the pre-holed specimens is 0.85 and 0.89 times their uniaxial compressive strength. Clearly, the stress thresholds of the specimens with a hole are much lower than that of the intact specimen. The reason why the crack initiation stress of the specimen RS-1 is a little larger than that of the specimen CS-3 is that the concentrated tensile stresses on the roof and floor of the holes are different (see Section 4). However, the crack damage stress of the specimen RS-1 is lower than that of the specimen CS-3, leading to an earlier failure. 4. Analysis of crack initiation mechanism 4.1. Stress distribution around a hole under compressive loads Studies on stresses around holes in infinite or finite plates have been long carried out. For circular holes and elliptic holes, the analytical solutions for surrounding stresses were found by Kirsch [69] and Inglis [70] using stress function method long ago, respectively. However, this

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Fig. 9. Axial stress, AE counts and cumulative AE counts versus time of specimens under uniaxial compression: (a) specimen IS-3, (b) specimen CS-3 and (c) specimen RS-1.

Fig. 10. Determination of crack closure stress for specimen RS-1: (a) schematic of axial displacement difference and (b) location of the point with the maximum axial displacement difference.

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Fig. 11. Conformal transformation: (a) external area of a complex shaped hole in z-plane; (b) external area of a unit circle in 𝜁 -plane.

method is not suitable for holes with complex geometries because of the difficulty in postulating the required stress functions. Literature indicated that the complex variable theory is particularly applicable to solve stress problems associated with complex shapes, multiply connected domains and high stress gradient [71–73]. In this section, the complex variable method combined with an optimization method was applied to find the asymptotic solutions for stresses around a circular and a rectangular hole. Similarly, the stress solutions of other shaped holes can also be derived by this method. 4.1.1. Complex variable theory To begin with, the rock materials are treated as continuous, homogeneous and isotropic media, and the underground opening is assumed to be deeply buried. Thus, the aforementioned stress problem can be simplified as a plane-strain problem. Muskhelishvili [74] systematically expounded the theory of complex function and developed a system of complex function equations of stress components for plane problems, namely [ ] { 𝜎𝑥 +𝜎𝑦 = 4Re 𝜑′ 1 (𝑧) [ ′′ ] (2) 𝜎𝑦 − 𝜎𝑥 + 2𝑖𝜏𝑥𝑦 = 2 𝑧𝜑 1 (𝑧) + 𝜓 ′ 1 (𝑧) where 𝜎 x , 𝜎 y and 𝜏 xy denote the stress components. z and i represent the complex variable and the imaginary unit, respectively. 𝜑1 (z), 𝜑1 ´(z), 𝜑1 ´´(z) and 𝜓 1 (z), 𝜓 1 ´(z), 𝜓 1 ´´(z) are the complex stress functions and their derivatives. Obviously, determining the two stress functions is the key to solving the stress components. According to the Riemann mapping theorem, the exterior of a hole in z-plane must be mapped to the exterior or interior of a unit circle in 𝜁 -plane by a mapping function z = w(𝜁 ) (z = x + iy, 𝜁 = 𝜉 + i𝜂), which is called conformal transformation (see Fig. 11) [75]. This makes it possible to solve the stress functions for complex shaped holes. As shown in Fig. 11, the hole suffers from a far-field stress: 𝜎 x ∞ =p, 𝜎 y ∞ =𝜆p and 𝜏 xy ∞ =0, 𝜆 denotes the lateral stress ratio. Defining the sign convention as positive for tensile stress and negative for compressive stress. Replacing Eq. (2) by the mapping function z = w(𝜁 ) and converting it into the form of polar coordinates yields { 𝜎𝜌 +𝜎𝜃 = 4Re[Φ(𝜁 )] [ ] 2 (3) 𝜎𝜃 −𝜎𝜌 + 2𝑖𝜏𝜌𝜃 = 𝜌22𝑤𝜁′ (𝜁) 𝑤(𝜁 )Φ′ (𝜁 ) + 𝑤′ (𝜁 )Ψ(𝜁 ) In Eq. (3), 𝑤(𝜁 ) and w’(𝜁 ) denote the conjugate and the first derivative of w(𝜁 ). Re(•) represents the real part of the complex variable. 𝜎 𝜌 , 𝜎 𝜃 and 𝜏 𝜌𝜃 are the three stress components (radial, tangential and shear stresses) in the orthogonal curvilinear system of a point in the z-plane, while (𝜌, 𝜃) denotes the polar coordinates of its mapping point in the 𝜁 plane. The two complex potential functions Φ(𝜁 ) and Ψ(𝜁 ) in Eq. (3) can be defined as below: 𝜑 ′ (𝜁 ) 𝜓 ′ (𝜁 ) Φ(𝜁 )=𝜑′1 (𝑧)= ′ , Ψ(𝜁 )=𝜓1′ (𝑧)= ′ (4) 𝑤 (𝜁 ) 𝑤 (𝜁 )

Here, 𝜑(𝜁 ) and 𝜓(𝜁 ) stands for two stress functions, which can be written as: { 𝜑(𝜁 )= − 2𝜋(11+𝜅) (𝑋 + 𝑖𝑌 ) ln 𝜁 + 𝐵𝑤(𝜁 ) + 𝜑0 (𝜁 ) ( ) (5) 𝜓 (𝜁 )= 2𝜋(1𝜅+𝜅) (𝑋 − 𝑖𝑌 ) ln 𝜁 + 𝐵 ′ + 𝑖𝐶 ′ 𝑤(𝜁 ) + 𝜓0 (𝜁 ) where X and Y are components of surface forces on the hole boundary, and their values are zero if the hole is unsupported. 𝜅 is a parameter that depends on material properties, 𝜅=(3−𝜇)/(1+𝜇); 𝜇 denotes the Poisson’s ratio of the surrounding rock mass. B=(𝜎 x ∞ +𝜎 y ∞ )/4=(𝜆+1)/4, Bʹ=(𝜎 y ∞ −𝜎 x ∞ )/2=(𝜆−1)/2, Cʹ=𝜏 xy ∞ =0; 𝜑0 (𝜁 ) and 𝜓 0 (𝜁 ) are two singlevalued analytic functions, which can be written in Laurent series form [72]: ∞ ⎧ 𝜑 (𝜁 ) = ∑ 𝑎 𝜁 − 𝑛 𝑛 ⎪ 0 𝑛=1 ∞ ⎨ ∑ ⎪𝜓0 (𝜁 )= 𝑏𝑛 𝜁 −𝑛 ⎩ 𝑛=1

(6)

The complex constant coefficients an and bn in Eq. (6) are real numbers if the geometry configuration of the hole exhibits symmetry about the x-axis or y-axis. In summary, the solution of the plane strain problem is equivalent to finding the above two single-valued analytic functions that satisfy the stress boundary condition shown below (𝜁 replaced by 𝜎 represents the point on the unit circle boundary). Generally, the power series method can be used to solve them. 𝜑0 (𝜎)+

𝑤(𝜎) 𝑤′ (𝜎)

( ) 𝜑′ 0 (𝜎) + 𝜓0 (𝜎) = −2𝐵 𝑤(𝜎) − 𝐵 ′ − 𝑖𝐶 ′ 𝑤(𝜎)

(7)

4.1.2. Optimization of mapping functions The conformal transformation function of any hole can be expressed as: 𝑧 = 𝜔(𝜁 ) = 𝑅(𝜁 +

∞ ∑ 𝑘=0

𝐶𝑘 𝜁 −𝑘 ),

|𝜁 | ≥ 1

(8)

where R denotes a positive real number that reflects the hole size. Ck (k = 1, 2, 3…∞) are complex constant coefficients related to the cavity shape. In most cases, the mapping function is accurate enough only to adopt a few terms of Ck [72]. C0 depends on the location of the hole in the coordinate system. From the above, the solution of the mapping function is reduced to finding the values of Ck and R. In this research, Ck are real numbers because both the circular and rectangular holes are symmetric about the x-axis (see Fig. 12), and the calculation steps are as follows: (i) First, divide the left half boundaries of the two holes into m (m = 30 in this study) equal parts, and thus the locations of (m + 1) sample points Aj (rj , 𝛼 j ) in an anti-clockwise direction are determined. Accordingly, (m + 1) mapping points Bj (1, 𝜃 j ) on the left half boundary of the unit circle can be formed, but their exact locations are unknown, that is, 𝜃 j needs to be solved as well.

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Fig. 12. Conformal mapping of the hole boundaries to the unit circle: (a) circular hole in z-plane; (b) rectangular hole in z-plane; (c) unit circle in 𝜁 -plane.

(ii) Next, based on Eq. (8), the relation of the point Aj and its mapping point Bj can be expressed as 𝑟𝑗 𝑒𝑖𝛼𝑗 = 𝑅(𝑒𝑖𝜃𝑗 +

∞ ∑ 𝑘=0

𝐶𝑘 𝑒−𝑖𝑘𝜃𝑗 )

(9)

(iii) Separate the real and imaginary parts of Eq. (9) with the aid of Euler’s formula, we have ∞ ⎧sin(𝛼 − 𝜃 ) + ∑ 𝐶 sin(𝛼 + 𝑘𝜃 ) = 0 𝑗 𝑗 𝑘 𝑗 𝑗 ⎪ 𝑘=0 [ ] ∞ ⎨ ∑ ⎪𝑟𝑗 = 𝑅 cos(𝛼𝑗 − 𝜃𝑗 ) + 𝐶𝑘 cos(𝛼𝑗 + 𝑘𝜃𝑗 ) ⎩ 𝑘=0

(10)

(iv) Substitute the starting point A1 (r1 , 0) and its mapping points B1 (1, 0) into Eq. (10), leads to /( ) 𝑛 −1 ∑ 𝑅 = 𝑟1 1+ 𝐶𝑘 (11) 𝑘=0

(v) If n terms of Ck are adopted, the question of satisfying Eq. (10) is equivalent to minimizing the following objective function f(Ck , 𝜃 j , R) with constraints: 𝑚 ( ) ∑ min 𝑓 𝐶𝑘 , 𝜃𝑗 , 𝑅 =

𝑠.𝑡. sin(𝛼𝑗 − 𝜃𝑗 ) + 𝑛∑ −1 𝑘=0

𝑛∑ −1 𝑘=0

𝑗=1

(

[ 𝑟𝑗 − 𝑅 cos(𝛼𝑗 − 𝜃𝑗 ) +

𝐶𝑘 sin(𝛼𝑗 + 𝑘𝜃𝑗 ) = 0

𝑘||𝐶𝑘 || < 1

𝑛∑ −1 𝑘=0

])2 𝐶𝑘 cos(𝛼𝑗 + 𝑘𝜃𝑗 )

⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭

(12) By substituting the coordinate values of (m + 1) sample points into Eq. (12), the values of Ck , 𝜃 j and R can be solved using an optimization method. Note that the closer the value of the objective function to zero, the more accurate the results are. In this paper, the powerful Boxʹs complex optimum method was first applied to solve the above problem [76]. The calculation procedures are described as follows: (i) Generate a polyhedron with N≥(n + 1) vertices, defined as an initial complex. Each vertex of the complex represents a set of possible solution. In this study, (0, 0, …, 0) was assigned to the starting vertex. For the remaining n vertices, their initial values can be randomly generated by ( ) 𝑥𝑖 𝑗 = 𝑎𝑖 + 𝜂𝑖 𝑗 𝑏𝑖 − 𝑎𝑖 (13) 𝑖 = 1, 2...𝑛, 𝑗 = 2, 3...𝑁 where (ai , bi )=(−1,1) is the interval of Ck . 𝜂 i j represents a pseudorandom number belonging to [0, 1]. (ii) Determine if these vertices meet the constraints. If there are q feasible vertices, then their centroid XC can be calculated by: 𝑋𝐶 =

𝑞 1∑ 𝑗 𝑋 𝑞 𝑗=1

(14)

In respect of the remaining (N-q) infeasible vertices, moving them to new locations corresponding to the midpoints of the lines between the centroid and the vertices, and then the original infea-

sible vertex Xq +1 is replaced by the new formed vertex X1 q +1 , namely ) 1( 𝐶 𝑋1 𝑞+1 = 𝑋 + 𝑋 𝑞+1 (15) 2 Repeat the above procedure until a satisfactory initial complex is formed. (iii) Evaluate the objective function values of the vertices and find the vertices Xj max and Xj min that have the maximum and minimum function values in the iteration. For each vertex, if 2 0.5 ∑ 𝑗 𝑗 { 𝑁1 𝑁 ≤ 𝑡𝑜𝑙𝑓 (the error tolerance tolf is 𝑗=1 [𝑓 (𝑋 ) − 𝑓 (𝑋min )] } 10−6 ), stop the iteration and proceed to the next step. If not, replace the worst vertex Xj max by its reflection vertex XR 1 max , defined as: ) ( (16) 𝑋𝑅1 max = 𝑋 (𝐶 ) + 𝜒 𝑋 (𝐶 ) − 𝑋 max where the reflection coefficient 𝜒 is 1.3 and X(C) denotes the centroid of the vertices excluding the worst vertex Xj max . Subsequently, checking whether the reflection vertex XR 1 max satisfies the variable interval and constraints. If yes and f(XR 1 max ) < f(Xj max ), replace Xj max by XR 1 max and maintain the remaining vertices of the complex. Otherwise, keep generating a new reflection vertex until it is eligible. The new reflection vertex can be obtained by the following equation: 𝑗 𝑋𝑅2 max = 𝜀𝑋𝑅1 max + (1 − 𝜀)𝑋min

(17)

where 𝜀 means the initial value of the reduction coefficient; it starts at 1 and is reduced by half in each reflection. Once a feasible reflection vertex of the worst point XR 1 max is found, go to the next iteration. (i) Return to the vertex Xj min as the optimization results. Thus, the values of the variables Ck , 𝜃 j and objective function can be derived. A MATLAB code was written to implement the above procedures, and the optimization results of the circular and the rectangular hole with different terms of Ck are listed in Table 4. The corresponding mapped shapes of holes are plotted in Fig. 13. Clearly, when their terms of Ck are 1 and 7, respectively, the predicted shapes are the most accurate. Therefore, the mapping functions are concluded as: 𝑧 = 𝑤(𝜁 ) = 10𝜁 𝑧 = 𝑤(𝜁 ) = 11.01𝜁 −

(18) 3.38 0.07 1.68 0.20 0.21 0.02 − − + + − 𝜁 𝜁2 𝜁3 𝜁4 𝜁5 𝜁6

(19)

4.1.3. Asymptotic solution for stress around a hole Circular hole: Substituting Eq. (18) into Eq. (7), the two single-valued analytic functions 𝜑0 (𝜁 ) and 𝜓 0 (𝜁 ) can be derived by comparing the coefficients of same powers of 𝜎 on the two sides of the equation. And then continue to substituting them into Eqs. (2)–(6), we have: ⎧ 2(𝜆−1)𝑝 cos2𝜃 ⎪𝜎𝜌 +𝜎𝜃 = (𝜆+1)𝑝 + [ 𝜌2 ] ⎪ 2 𝜎𝜌 + 2𝑖𝜏𝜌𝜃 = 1 − 𝜌2 + 𝜌34 (𝜆 − 1)𝑝cos2𝜃 + ⎨𝜎 𝜃 − ] ) ⎪ [( 2 3 ⎪+𝑖 𝜌2 + 1 − 𝜌4 (𝜆 − 1)𝑝 sin 2𝜃 ⎩

(𝜆+1)𝑝 𝜌2

(20)

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Table 4 Optimization results of Ck with different terms. Hole shape

Terms /n

C0

C1

C2

C3

C4

C5

C6

C7

Circle

1 2 3 4 5 6 3 4 5 6 7 8

0 0 0 0 0 0 0.0019 −0.0044 0.0012 −0.0020 0.0002 −0.0067

× 0 0 0 0 0 −0.3985 −0.3009 −0.3046 −0.3220 −0.3073 −0.3002

× × 0 0 0 0 −0.0316 0.0133 0.0034 0.0000 −0.0060 0.0002

× × × 0 0 0 × −0.1345 −0.1248 −0.1431 −0.1527 −0.1372

× × × × 0 0 × × 0.0017 −0.0080 0.0174 −0.0015

× × × × × 0 × × × 0.0433 0.0194 −0.0073

× × × × × × × × × × −0.0019 0.0027

× × × × × × × × × × × 0.0258

Rectangle

Fig. 13. Comparison between the actual and mapped shapes of holes with different terms of Ck : (a) circular hole; (b) rectangular hole.

By solving the above equation and performing coordinate system conversion, the stress components can be given as: ( ) ( ) 2 ⎧ (𝜆+1)𝑝 )𝑝 4𝑅2 3𝑅4 1 − 𝑅𝑟2 − (𝜆−1 1 cos2𝛼 − + ⎪𝜎 𝑟 = 2 2 4 2 ( ) ( 𝑟 4) 𝑟 ⎪ (𝜆+1)𝑝 (𝜆−1)𝑝 𝑅2 3𝑅 (21) = + + 𝛼 1 + 1 cos2 𝜎 ⎨ 𝛼 2 2) 𝑟2 𝑟4 ( ⎪ (𝜆−1)𝑝 2𝑅2 3𝑅4 1 + 𝑟2 − 𝑟4 sin2𝛼 ⎪𝜏𝑟𝛼 = 2 ⎩ It can be seen that Eq. (21) is completely agreeable with the famous Kirsch equation if their axes are in the same direction [69]. This indicates the complex variable method for stress solution is reliable. According to Eq. (21), the hoop stress distribution on the hole boundary is illustrated in Fig. 14. It is observed that the concentrated tensile stresses on the roof and floor gradually change to compressive stress as 𝜆 increases, whilst it is the opposite on the sidewalls. Besides, both the radial and hoop stresses around the hole can also be obtained. With the increase of the distance from the hole boundary, the radial stresses on the top and bottom locations of the hole gradually turn to p, while the tangential stresses progressively approach 𝜆p. In contrast, the situation is contrary on the two sides of the hole. In conclusion, when the distance exceeds five times the radius, the stress distribution is basically not affected by the hole. Rectangular hole: According to Eq. (19), the transformation forms of w(𝜎) can be written as: 𝑤(𝜎) = 11.01𝜎 − 3.38𝜎 −1 − 0.07𝜎 −2 − 1.68𝜎 −3 + 0.19𝜎 −4 + 0.21𝜎 −5 −0.02𝜎

−6

𝑤 (𝜎) = 11.01 + 3.38𝜎 ′

+0.13𝜎

−7

+0.13𝜎

−3

+ 5.04𝜎

−4

− 0.77𝜎

−5

− 1.07𝜎

𝑤(𝜎) = 11.01𝜎 −1 − 3.38𝜎 − 0.07𝜎 2 − 1.68𝜎 3 + 0.19𝜎 4 + 0.21𝜎 5 − 0.02𝜎 6 (22c) 𝑤′ (𝜎) = 11.01 + 3.38𝜎 2 +0.13𝜎 3 + 5.04𝜎 4 − 0.77𝜎 5 − 1.07𝜎 6 + 0.13𝜎 7 (22d)

(22a) −2

Fig. 14. Hoop stress distribution on the boundary of circular hole.

−6

(22b)

Substituting the above equations into Eq. (7), then the values of an and bn in Eq. (6) can be calculated using the powerseries expansion method [77], namely., a1 = (−3.25𝜆 + 6.27)p, a2 = (−0.03𝜆 + 0.15)p, a3 = (0.78𝜆 + 0.96)p, a4 = -(0.09𝜆 + 0.11)p,

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

of 1.25 m were placed on the roof of the hole, and the other ten monitoring points from L1 to L10 were arranged every 2.5 m on the left side of the hole. The horizontal stresses of points T and the vertical stresses of points L can be recorded by a self-written code. It is found that the numerical results are in good agreement with the theoretical results, which is shown in Fig. 17b. The slight difference results from the large size of the element in numerical model. This is because the numerical stress value is actually the value of the centroid of the element where the monitoring point is located. The error can be reduced using smaller element, but it will cause difficulties for computer operation [78].

4.3. Crack initiation mechanism

Fig. 15. Tangential stress distribution on the boundary of rectangular hole.

a6 = 0.01(𝜆 + 1)p, a7 = a8 =…= an = 0, b1 = a5 = −0.11(𝜆 + 1)p, (3.33𝜆 + 9.35)p, b2 = (0.08𝜆 − 0.09)p, b3 = (−3.46𝜆 + 5.73) p, …, b99 = (−1.99 × 10−6 𝜆 + 4.75 × 10−6 )p, b100 = (5.30 × 10−7 𝜆 −3.56 × 10−6 )p. Note that the number of items in an is six, but that in bn is infinite. It is found that, when n is large enough, the value of bn is approximately zero. Therefore, in this work, the terms of bn is set to 100, which has little effect on the solution of stresses. Introducing an and bn into source equations of 𝜑0 (𝜁 ) and 𝜓 0 (𝜁 ), and then substituting them into Eq. (5). The two complex stress functions are concluded as: [ ] .13)𝑝 2.75(𝜆 + 1)𝑝𝜁 + (5.42−4𝜁 .10𝜆)𝑝 − (0.05𝜆−0 + ⎧ 𝜁2 𝜑 𝜁 = 𝑝 ( ) ⎪ .04𝜆)𝑝 1)𝑝 (0.54+0.36𝜆)𝑝 )𝑝 − (0.06+0 − 0.05(𝜁𝜆+1 + 0.01(𝜁𝜆+ 5 6 ⎪ 𝜁3 𝜁4 (23) 7 . 66+5 . 02 𝜆) 𝑝 −0 . 06+0 . 04 𝜆) 𝑝 ( ( ⎨ ⎡5.51(𝜆 − 1)𝑝𝜁 − + +⎤ ⎪ 𝜁 𝜁2 ⎢ ⎥ 𝜓 𝜁 == 𝑝 ( ) −6 −7 +5.30×10 𝜆) ⎪ ⎢ (6.57−43.30𝜆)𝑝 + .... + 𝑝(−3.56×10 100 ⎥ ⎩ 𝜁 𝜁 ⎣ ⎦ Combining Eqs. (23), (4) and (3), the stress components 𝜎 𝜌 , 𝜎 𝜃 and 𝜏 𝜌𝜃 can be solved, and the stress distributions around the rectangular hole with different lateral stress ratios and distances are illustrated in Figs. 15 and 16. From the above figures, it is observed that the stress distribution around the rectangular hole is quite different from that of the circular hole. This indicates that the shape of the hole has a pronounced effect on the stress distribution. However, the variation laws of the surrounding stresses of the two holes with the lateral pressure coefficient and the distance are similar. 4.2. Validation of asymptotic stress solution For the circular hole, its stress solution obtained by the complex variable theory is strict and is exactly the same as the Kirsch equation. In contrast, the stress solution for the rectangular hole is asymptotic because the formulas of the single-valued analytic function 𝜓 0 (𝜁 ) and the mapping function were simplified in the calculation process. To further verify the reliability of the complex variable method, a numerical software Flac was used to simulate the stress distribution around the rectangular hole. According to the plane strain assumption, a prismatic model with dimensions of 120 m × 0.5 m × 120 m (L × W × H) was built, as shown in Fig. 17a. The mesh size is equal to the width of the model, leading to 206,440 nodes and 644,895 tetrahedral elements. The stress field was assumed as: 𝜎 x ∞ = 0, 𝜎 y ∞ = 20 MPa and 𝜏 xy ∞ =0, and the boundaries were fixed along x- or y- direction. To avoid plastic deformation of the model, the Young’s modulus was set large enough and the isotropic elastic constitutive model was adopted in this study. During the modelling, ten monitoring points (T1, T2, …, T10) with a spacing

As stated in the above experimental study, cracks of the pre-holed specimens under uniaxial loading proceed from primary tensile cracks through slabbing fractures and secondary tensile cracks to shear cracks. The initiation mechanisms of these cracks can be explained based on the stress distributions around the holes, and the details are described as below. Generally, under uniaxial loading, the intact specimens are stretched in both tangential and radial directions because of the lateral expansion. Thus, the splitting cracks such as cracks 1a and 1b in Fig. 6a can be formed gradually with the increasing load. Due to the end friction effect, the ends of the specimen are in a three-dimensional stress state. When the applied load increases to a certain value, a shear band will appear and eventually grow into a shear crack 3 with an inclination of (45°+ 𝜙/2), as shown in Fig. 8. Finally, the failure induced by the coalescence of the shear crack and the splitting tensile cracks occurs. For specimens containing a hole, the analytical stress distributions around the holes agree well with the strain distributions in the tests. As can be seen in Figs. 14 and 15, the maximum concentrated tensile stresses on the roofs (or floors) of the circular and the rectangular hole under uniaxial compressive loads (𝜆=0) are p and 0.92p, respectively. As the load increases, the tensile stresses increase accordingly. When they reach a critical value to overcome the cohesion of the specimens, the primary tensile cracks 1a and 1b shown in Fig. 6 start to initiate and propagate stably in a load-controlled manner. Obviously, the cracks around the circular hole appear earlier than that of the rectangular hole. This is why the crack initiation stress of the specimen CS-3 is smaller than that of the specimen RS-1. According to the asymptotic solutions for stresses, the concentrated tensile stresses decrease with the increase of the distance from the hole boundary, and gradually change to compressive stress. Once the tensile stress is reduced to zero, the primary tensile cracks will stop propagating. Based on this criterion, the theoretical lengths of the cracks for the circular and the rectangular holes can be calculated as 0.73r1 and 1.79 r1 respectively. By comparison, the average actual values of the specimens CS-3 and RS-1 measured by the DIC technique are 7.62 mm and 11.60 mm, respectively, which are basically agreeable with the theoretical results. In contrast, the two sides of the circular and the rectangular holes are concentrated by compressive stresses, with maximum values of 3p and 1.81p, respectively. When they reach the compressive strength of the rock, the slabbing fractures will emerge. As the primary tensile cracks propagate away from the cavity, the critical stress region shifts from the tips of the primary tensile cracks to a region on their either side. This leads to the nucleation of the secondary tensile cracks, and they develop rapidly along the loading direction and eventually intersect with the slabbing fracturing zones. Besides, compared with the specimen CS3, the secondary tensile cracks of the specimen RS-1 occur later than the slabbing fractures. The reason may be that, the long lengths of the primary tensile cracks take a long propagation time and the straight-walltype hole boundary is more prone to slabbing failure. For the initiation mechanism of the shear cracks in the pre-holed specimens, it is the same as that in the intact specimen.

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Fig. 16. Stress distributions at different distances from the boundary of rectangular hole: (a) one time the polar radius; (b) three times the polar radius and (c) five times the polar radius.

5. Discussion 5.1. Reliability of DIC technique As mentioned above, the cracking processes of specimens under uniaxial loading can be visually displayed using the DIC technique, and the strain localization zones on the deformed speckle images are good representations of the potential cracks. To verify the accuracy of the DIC

method for monitoring the strain, we strategically attached six strain gauges on the back surface of an intact specimen, three of which were used for horizontal strain and three for vertical strain. The locations of these strain gauges, denoted as A, B, …and F, are shown in Fig. 18a. A DH3818S static strain meter was used to record the real-time strains during the tests, and it worked simultaneously with other experimental apparatus. As a comparison, we placed six corresponding monitoring points (Aʹ, Bʹ, …, Fʹ) on the front surface of the sample, and their strain

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Fig. 17. Validation of the asymptotic solution of stress for rectangular hole: (a) locations of the monitoring points; (b) comparision between the theoretical and the numerical results.

5.2. Evolution characteristics of primary tensile cracks

Fig. 18. Validation of DIC technique: (a) locations of monitoring points on the back and ftont surfaces; (b) comparision of strain values between the strain gauge and DIC methods.

values can be derived using the GOM correlate software. The results of the two methods are compared in Fig. 18b. Clearly, there is a good correlation between them in the early stages of loading. At the end of the experiment, it is observed that the strain values derived by the DIC method are a little different from the actual values, which may be related to the anisotropy and heterogeneity of the rock. However, for the significant error of point Eʹ in the plastic deformation stage, it resulted from surface spalling. To sum up, the DIC technique is reliable and powerful for the analysis of cracking behavior.

During the tests on specimens with a hole, an interesting phenomenon that the primary tensile cracks experience a progressive process of ″open-close-reopen″ was observed; that is to say, the primary tensile cracks propagate steadily at stage BC, and then they close gradually after the appearance of the secondary tensile cracks at stage CD. However, they open again at the end of the experiments. This behavior can be explained as follows. First, the cracks initiate when the axial stress exceeds the crack initiation stress, and propagate along the loading direction with the increasing stress. Next, the secondary tensile cracks appear, leading to lateral compression on the primary tensile cracks. Therefore, the primary tensile cracks close gradually as the secondary tensile cracks propagate. When the failure occurs, the lateral compression diminishes. Thus, the tensile cracks reopen. To clearly reveal the propagation characteristic of the primary tensile cracks, four monitoring points with a spacing of 2 mm were evenly set on the propagation path of each crack to record the horizontal strain by the DIC technique. The strain-time curves of these points (i.e., R1, R2, R3 and R4 on the roof of the hole; F1, F2, F3 and F4 on the floor of the hole) in combination with the stress-time curves of the specimens CS-3 and RS-1 are shown in Fig. 19. Observations show that the propagation of the tensile cracks is progressive, namely the strain of the monitoring point near the hole increases first and its value is the largest. At stage OA, the strains are all approximately zero because no significant deformation occurs. When the stage AB is reached, the strains increase linearly with the increase of the axial stress. During the steady propagation of the primary tensile cracks at stage BC, the magnitude of the linear increase in strain becomes larger. When the secondary tensile cracks appear, the strains drop abruptly, indicating the closure of the primary tensile cracks. Once the failure occurs, the strain will increase dramatically. This suggests that the lateral constraint on the primary tensile cracks is removed and the cracks open again. To conclude, the laboratory tests and the theoretical analysis of the crack evolution increase the understanding of the failure mechanism of hard rock with a hole. Also, it can be seen that the mechanical properties of the rocks are degraded by the presence of the hole. Literature shows that the weaken effect depends on many factors such as the size, shape, number and distribution of holes, the material properties (e.g., lithology, porosity, and heterogeneity) of rocks and the loading conditions (uniaxial, biaxial and triaxial compression, shearing, tension and dynamic loads) [39,41,48–56,60,61,63,79,80]. However, for brittle rock specimen containing one hole under uniaxial loading, their crack evolution processes are the same, but the stress thresholds are different unless

H. Wu, G. Zhao and W. Liang

International Journal of Mechanical Sciences 163 (2019) 105084

Fig. 19. Stress and horizontal strain versus time for monitoring points: (a) specimen CS-3 and (b) specimen RS-1.

their shape are similar. In contrast, the stress distribution is independent of the hole size, but related to the hole shape. Note that the theoretical stress distribution around any shaped hole can also be obtained employing the above analytical method. In respect of the failure characteristics of rock tunnels or roadways in underground mines, the experimental results are applicable if their principal stresses are in the same direction. That is why spalling failure or rock burst often occur on the sidewalls [28–34]. Considering the various shapes and different stress states of actual tunnels in rock engineering, we will try to further investigate the corresponding mechanical behaviour of pre-holed specimens under triaxial loading in future research.

6. Conclusions In the present study, the mechanical properties and cracking behavior of sandstone specimens containing a circular hole or a rectangular hole loaded in compression were experimentally studied incorporating the DIC and AE techniques. Afterwards, the initiation mechanisms of the different types of cracks were analyzed based on the stress distributions obtained by the complex variable method. The experimental results indicate that the mechanical parameters including the uniaxial compressive strength, Young’s modulus, peak strain and brittleness index are remarkably weakened by the hole and influenced by its shape. Cracks in the pre-holed specimens under uniaxial compression evolve from primary tensile cracks through secondary tensile cracks and slabbing fractures to shear cracks, which can be clearly reproduced by the reliable DIC technique. The damage characteristics reflected by the AE signals agree well with the cracking behavior. Furthermore, several stress thresholds, i.e., the crack closure stress, crack initiation stress and crack damage stress, are estimated. The ultimate failure of the specimen with a hole is caused by the coalescence of the shear cracks and the slabbing fracturing zones. The theoretical results show that the stress distribution around the hole depends on the hole shape and lateral stress ratio, and the stress influence zone is about five times the polar radius. Under the condition of uniaxial compression, tensile stresses concentrate on the roofs and floors of the holes and gradually become compressive stresses as the distance from the hole boundary increases, while the compressive stress concentration occurs on the sides of the holes. This reveals the initiation mechanisms of cracks well. Besides, the evolution characteristic of ″open-close-reopen″of the primary tensile cracks is analyzed in detail based on the variation of the horizontal tensile strain, and the lengths

of the primary tensile cracks are also determined both experimentally and theoretically. Acknowledgements This work was financially supported by the National Natural Science Foundation of China (No. 51774321), the National Key Research and Development Program of China (No. 2018YFC0604606), and the Fundamental Research Funds for the Central Universities of Central South University (No. 2018zzts215). The authors are deeply indebted to Dr. Xiaoli Zhang and Dr. Xiangtai Zeng for their kind assistance with the complex variable theory. Also, we would like to express our sincere thanks to the reviewers for their constructive comments and the editors for their hard work. References [1] Bombolakis EG. Photoelastic stress analysis of crack propagation within a compressive stress field Ph.D. Thesis. Massachusetts Institute of Technology; 1963. [2] Lajtai EZ. Brittle fracture in compression. Int J Fract 1974;10:525–36. [3] Ingraffea AR, Heuze FE. Finite element models for rock fracture mechanics. Int J Num Anal Method Geomech 1980;4:25–43. [4] Huang JF, Chen GL, Zhao YH, Wang R. An experimental study of the strain field development prior to failure of a marble plate under compression. Tectonophysics 1990;175:283–90. [5] Hoek E, Bieniawski ZT. Brittle fracture propagation in rock under compression. Int J Fract 1965;1:137–55. [6] Boet A. The initiation of secondary cracks in compression. Eng Fract Mech 2000;66:187–219. [7] Wong LNY, Einstein HH. Systematic evaluation of cracking behavior in specimens containing single flaws under uniaxial compression. Int J Rock Mech Min Sci 2009;46:239–49. [8] Yang SQ, Jing HW. Strength failure and crack coalescence behavior of brittle sandstone samples containing a single fissure under uniaxial compression. Int J Fract 2011;168:227–50. [9] Ashby MF, Hallam SD. The failure of brittle solids containing small cracks under compressive stress states. Acta Metall 1986;34(3):497–510. [10] Horri H, Nemat-Nasser S. Brittle failure in compression: splitting, faulting and brittle-ductile transition. Philos Trans R Soc A 1986;319(1549):337–74. [11] Boet A, Einstein HH. Fracture coalescence in rock-type materials under uniaxial and biaxial compression. Int J Rock Mech Min Sci 1998;35:863–88. [12] Wong RHC, Chau KT. Crack coalescence in a rock-like material containing two cracks. Int J Rock Mech Min Sci 1998;35:147–64. [13] Wong LNY, Einstein HH. Crack coalescence in molded gypsum and carrara marble: part 1. Macroscopic observations and interpretation. Rock Mech Rock Eng 2009;42:475–511. [14] Lee H, Jeon S. An experimental and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression. Int J Solids Struct 2011;48:979–99. [15] Wang YT, Zhou XP, Shou YD. The modeling of crack propagation and coalescence in rocks under uniaxial compression using the novel conjugated bond-based peridynamics. Int J Mech Sci 2017;128–129:614–43.

H. Wu, G. Zhao and W. Liang [16] Zhou XP, Bi J, Qian QH. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws. Rock Mech Rock Eng 2015;48:1097–114. [17] Yang XX, Kulatilake PHSW, Jing HW, Yang SQ. Numerical simulation of a jointed rock block mechanical behavior adjacent to an underground excavation and comparison with physical model test results. Tunn Undergr Sp Technol 2015;50:129–42. [18] Cao RH, Cao P, Lin H, Pu CZ, Ou G. Mechanical behavior of brittle rock-like specimens with pre-existing fissures under uniaxial loading: experimental studies and particle mechanics approach. Rock Mech Rock Eng 2016;49:763–83. [19] Sagong M, Bobet A. Coalescence of multiple flaws in a rock-model material in uniaxial compression. Int J Rock Mech Min Sci 2002;39:229–41. [20] Bahaaddini M, Sharrock G, Hebblewhite BK. Numerical investigation of the effect of joint geometrical parameters on the mechanical properties of a non-persistent jointed rock mass under uniaxial compression. Comput Geotech 2013;49:206–25. [21] Dyskin AV, Sahouryeh E, Jewell RJ, Joer H, Ustinov KB. Influence of shape and locations of initial 3-D cracks on their growth in uniaxial compression. Eng Fract Mech 2003;70:2115–36. [22] Zhang ZN, Wang DY, Zheng H, Ge XR. Interactions of 3D embedded parallel vertically inclined cracks subjected to uniaxial compression. Theor Appl Fract Mech 2012;61:1–11. [23] Dang WG, Konietzky H, Chang LF, Frühwirt T. Velocity-frequency-amplitude-dependent frictional resistance of planar joints under dynamic normal load (DNL) conditions. Tunn Undergr Sp Technol 2018;79:27–34. [24] Dang WG, Konietzky H, Frühwirt T, Herbst M. Cyclic frictional responses of planar joints under cyclic normal load conditions: laboratory tests and numerical simulations. Rock Mech Rock Eng 2019. doi:10.1007/s00603-019-01910-9. [25] Li XB, Zhou T, Li DY. Dynamic strength and fracturing behavior of single-flawed prismatic marble specimens under impact loading with a split-hopkinson pressure bar. Rock Mech Rock Eng 2017;50(1):29–44. [26] Aliha MRM, Ayatollahi MR, Smith DJ, Pavier MJ. Geometry and size effects on fracture trajectory in a limestone rock under mixed mode loading. Eng Fract Mech 2010;77(11):2200–12. [27] Yang SQ, Dai YH, Han LJ, Jin ZQ. Experimental study on mechanical behavior of brittle marble samples containing different flaws under uniaxial compression. Eng Fract Mech 2009;76:1833–45. [28] Ranjith PG, Zhao J, Ju MH, De Silva RVS, Rathnaweera TD, Bandara AKMS. Opportunities and challenges in deep mining: a brief review. Eng PRC 2017;3:546–51. [29] He MC. Latest progress of soft rock mechanics and engineering in China. J Rock Mech Geotech Eng 2014;6:165–79. [30] Yin TB, Bai L, Li X, Li XB, Zhang SS. Effect of thermal treatment on the mode I fracture toughness of granite under dynamic and static coupling load. Eng Fract Mech 2018;199:143–58. [31] Liang WZ, Zhao GY, Wu H, Dai B. Risk assessment of rockburst via an extended Mabac method under fuzzy environment. Tunn Undergr Sp Techol 2018;83:533–44. [32] Gong FQ, Luo Y, Li XB, Si XF, Tao M. Experimental simulation investigation on rockburst induced by spalling failure in deep circular tunnels. Tunn Undergr Sp Technol 2018;81:413–27. [33] Tao M, Li XB, Wu CQ. 3D numerical model for dynamic loading-induced multiple fracture zones around underground cavity faces. Comput Geotech 2013;54:33–45. [34] Du K, Tao M, Li XB, Zhou J. Experimental study of slabbing and rockburst induced by true-triaxial unloading and local dynamic disturbance. Rock Mech Rock Eng 2016;49:3437–53. [35] Hoek E. Rock fracture under static stress conditions Ph.D. Thesis. University of Cape Town; 1965. [36] Lajtai EZ, Carter BJ, Duncan EJS. Mapping the state of fracture around cavities. Eng Geol 1991;31:277–89. [37] Martin CD. Strength of massive Lac du Bonnet granite around underground openings Ph.D. Thesis. University of Manitoba; 1993. [38] Carter BJ, Lajtai EZ, Petukhov A. Primary and remote fracture around underground cavities. Int J Numer Anal Methods Geomech 1991;15:21–40. [39] Zhou ZL, Tan LH, Cao WZ, Zhou ZY, Cai X. Fracture evolution and failure behaviour of marble specimens containing rectangular cavities under uniaxial loading. Eng Fract Mech 2017;184:183–201. [40] Dzik EJ, Lajtai EZ. Primary fracture propagation from circular cavities loaded in compression. Int J Fract 1996;79:49–64. [41] Tang CA, Wong RHC, Chau KT, Lin P. Modeling of compression-induced splitting failure in heterogeneous brittle porous solids. Eng Fract Mech 2005;72:597–615. [42] Carter BJ, Lajtai EZ, Yuan Y. Tensile fracture from circular cavities loaded in compression. Int J Fract 1992;57:221–36. [43] Mastin RJ. Development of borehole breakouts in sandstone Ph.D. Thesis. Stanford University; 1984. [44] Carter BJ. Size and stress gradients effects on fracture around cavities. Rock Mech Rock Eng 1992;25:167–86. [45] Martin DC. Seventeenth Canadian geotechnical colloquium: the effect of cohesion loss and stress path on brittle rock strength. Can Geotech J 1997;34:698–725. [46] Martin CD, Kaiser PK, Mccreath DR. Hoek–Brown parameters for predicting the depth of brittle failure around tunnels. Can Geotech J 1999;34:136–51. [47] Fakhimi A, Carvalho F, Ishida T, Labuz JF. Simulation of failure around a circular opening in rock. Int J Rock Mech Min Sci 2002;39:507–15.

International Journal of Mechanical Sciences 163 (2019) 105084 [48] Wang SY, Sloan SW, Sheng DC, Tang CA. Numerical analysis of the failure process around a circular opening in rock. Comput Geotech 2012;39:8–16. [49] Huang YH, Yang SQ, Tian WL, Zhao J, Ma D, Zhang CS. Physical and mechanical behavior of granite containing pre-existing holes after high temperature treatment. Arch Civ Mech Eng 2017;17:912–25. [50] Zhu WC, Liu J, Tang CA, Zhao XD, Brady BG. Simulation of progressive fracturing processes around underground excavations under biaxial compression. Tunn Undergr Sp Techol 2005;20(3):231–47. [51] Weng L, Li X, Taheri A, Wu QH, Xie XF. Fracture evolution around a cavity in brittle rock under uniaxial compression and coupled static–dynamic loads. Rock Mech Rock Eng 2018;51:531–45. [52] Haeri H, Khaloo A, Marji MF. Fracture analyses of different pre-holed concrete specimens under compression. Acta Mech Sin 2015;31:855–70. [53] Steen BVD, Vervoort A, Napier JAL. Observed and simulated fracture pattern in diametrically loaded discs of rock material. Int J Fract 2005;131:35–52. [54] Lajtai EZ, Lajtai VN. The collapse of cavities. Int J Rock Mech Min Sci 1975;12:81–6. [55] Lin P, Wong RHC, Tang CA. Experimental study of coalescence mechanisms and failure under uniaxial compression of granite containing multiple holes. Int J Rock Mech Min Sci 2015;77:313–27. [56] Jespersen C, Maclaughlin M, Hudyma N. Strength, deformation modulus and failure modes of cubic analog specimens representing macroporous rock. Int J Rock Mech Min Sci 2010;47:1349–56. [57] Li YP, Chen LZ, Wang YH. Experimental research on pre-cracked marble under compression. Int J Solids Struct 2005;42:2505–16. [58] Yin Q, Jing HW, Su HJ. Investigation on mechanical behavior and crack coalescence of sandstone specimens containing fissure-hole combined flaws under uniaxial compression. Geosci J 2018;22:825–42. [59] Fan X, Li KH, Hongpeng Lai HP, Xie YL, Cao RH, Zheng J. Internal stress distribution and cracking around flaws and openings of rock block under uniaxial compression: a particle mechanics approach. Comput Geotech 2018;102:28–38. [60] Li DY, Zhu QQ, Zhou ZL, Li XB, Ranjith PG. Fracture analysis of marble specimens with a hole under uniaxial compression by digital image correlation. Eng Fract Mech 2017;183:109–24. [61] Zhu QQ, Li DY, Han ZY, Li XB, Zhou ZL. Mechanical properties and fracture evolution of sandstone specimens containing different inclusions under uniaxial compression. Int J Rock Mech Min Sci 2019;115:33–47. [62] Ulusay R, Hudson JA. The complete ISRM suggested methods for rock characterization, testing and monitoring: 1974–2006. ISRM Turkish National Group; 2007. [63] Ma SW, Luo ZQ, Wu H, Qin YG. Mechanical properties of rock with intersection structures and its progressive failure mechanism. IEEE Access 2019;7(1):60920–30. [64] Song H, Zhang H, Fu D, Kang Y, Huang G, Ou C, Cai Z. Experimental study on damage evolution of rock under uniform and concentrated loading conditions using digital image correlation. Fatigue Fract Eng Mater Struct 2013;36:760–8. [65] Pan B, Qian K, Xie H, et al. Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review. Meas Sci Technol 2009;20:62001. [66] Munoz H, Taheri A, Chanda EK. Rock drilling performance evaluation by an energy dissipation based rock brittleness index. Rock Mech Rock Eng 2016;49:3343–55. [67] Cai M, Kaiser PK, Tasaka Y, Maejima T, Morioka H, Minami M. Generalized crack initiation and crack damage stress thresholds of brittle rock masses near underground excavations. Int J Rock Mech Min Sci 2004;41:833–47. [68] Hoek E, Martin CD. Fracture initiation and propagation in intact rock-A review. J Rock Mech Geotech Eng 2014;6:287–300. [69] Kirsch EG. Die theorie der elasitzität und die bedürfnisse der festigkeitslehre. Z Ver Dtsch Ing 1898;42:797–807. [70] Inglis CE. Stresses in a plate due to the presence of cracks and sharp corners. Proc Inst Nav Arch 1913;55:219–30. [71] Exadaktylosa GE, Stavropoulou MC. A closed-form elastic solution for stresses and displacements around tunnels. Int J Rock Mech Min Sci 2002;39:905–16. [72] Lu AZ, Zhang LQ. Complex function method on mechanical analysis of underground tunnel. Beijing: Science Press; 2007. [73] Sharma DS. Stress distribution around polygonal holes. Int J Mech Sci 2012;65:115–24. [74] Muskhelishvili NI. Some basic problems of the mathematical theory of elasticity. Amsterdam: Kluwer Academic; 1977. [75] Riemann Bernhard. Grundlagen für eine allgemeine Theorie der Functionen einer veränderlichen complexen Grösse. Göttingen: George-August-University of Göttingen; 1851. [76] Sgrott OL, Noriler D, Wiggers VR, Meier HF. Cyclone optimization by complex method and CFD simulation. Powder Technol 2015;277:11–21. [77] Chen ZY. Analytical method of rock mechanics analysis 1994. [78] Dang WG, Wu W, Konietzky H, Qian JY. Effect of shear-induced aperture evolution on fluid flow in rock fractures. Comput Geotech 2019;114:103151. doi:10.1016/j.compgeo.2019.103152. [79] Wu H, Zhao GY, Kulatilake PHSW, Liang WZ, Wang EJ. Fracturing behaviour of sandstone specimens with a cavity formed by intersecting excavations under compression: experimental study and numerical modelling. Strain 2019;103:242–53 e12316. [80] Zhong ZB, Deng RG, Lv L, Fu XM, Yu J. Fracture mechanism of naturally cracked rock around an inverted U-shaped opening in a biaxial compression test. Int J Rock Mech Min Sci 2018;103:242–53.