Abstract
Amygdaloidal basalt, as a heterogeneous rock, is widely exposed at Baihetan hydropower station, China. The geometric effect of amygdales needs further studies and quantifying the shape, orientation, and statistical distribution of amygdales plays an important role in the laboratory and numerical experiments. Therefore, digital image processing (DIP) was first utilized to build a heterogeneous model (HM) to calibrate against the laboratory test results. Then, the heterogeneous models (HMs) with prescribed geometric features were generated by the inverse Monte-Carlo (IMC) algorithm. The uniaxial compression experiments based on HMs were conducted to study the mechanism of the crack initiation and propagation in the amygdaloidal basalt. The tensile fractures were mainly occurred in the matrix, and the shear fractures were mainly occurred in the amygdales. With the increase in the elliptic coefficient of amygdales, the uniaxial compressive strength (UCS) showed a linear growth trend. With the increase in the orientation of amygdales, the UCS exhibited a “V-shaped” distribution characteristic. This paper provides a numerical method for studying the mechanical properties of rocks with flaws.
1. Introduction
Natural rocks and rock-like materials with preexisting flaws and cracks may exhibit unique damage evolution characteristics, and they may reduce the mechanical properties and strength of rocks [1, 2]. The crack initiation and growth of rocks and rock-like materials during uniaxial compression are significantly influenced by preexisting flaws, cracks, and soft inclusions. Investigating the damage evolution process of various rock samples is beneficial to implement better designs of underground geomechanical engineering [3–5]. More importantly, quantifying the statistical characteristics of preexisting flaws and cracks plays important role in the laboratory and numerical experiments, and they have been considered to investigate the influence on damage evolution vital discrete fracture network (DFN) and grain-based model (GBM). The related studies exhibited that tensile cracks caused by local stress concentration around preexisting flaws and cracks affected significantly the deformation characteristics and led to the brittle failure of rock material [6, 7].
The numerical methods based on continuum mechanics cannot truly characterize the preexisting cracks and flaws in heterogeneous rocks and simulate the crack behavior during the fracture process [8, 9]. In addition, the traditional numerical methods based on noncontinuum mechanics cannot effectively capture the continuous-discontinuous trans-scale fracture process in rock masses [10]. The combined finite-discrete element method (FDEM) [11], as an advanced numerical method, can simulate the physical phenomena during the fracture process. Mahabadi [12] investigated the effect of mesostructure on the fracturing and mechanical behavior in rock materials based on the FDEM code, Y2D. Cai [13] conducted the Brazilian splitting experiments by ELFEN software and investigated the influence of the friction coefficient of structures on the crack initiation and growth of rock specimens. Zhang et al. [14, 15] deeply studied the fracture characteristics and failure evolution process of amygdaloidal basalt under different confining stresses based on the continuous-discontinuous model. However, the effect of the microgeometric characteristics of amygdales on the crack behaviors during the failure process of amygdaloidal basalt needs to be explored deeply.
This framework investigated the damage evolution in amygdaloidal basalt based on laboratory test results and FDEM. First, digital image processing (DIP) [16, 17] was utilized to build a heterogeneous model (HM) to calibrate against the laboratory test results. Then, the heterogeneous models (HMs) with prescribed geometric features were generated by the inverse Monte-Carlo (IMC) algorithm. The uniaxial compression experiments based on HMs were conducted to study the mechanism of the crack initiation and propagation in the amygdaloidal basalt.
2. HM in FDEM
2.1. Basic Theory of FDEM
This study utilizes the FDEM-based Irazu2D software. In the FDEM, the modeling zone is composed of plane strain triangular elements and cohesive crack elements (CCEs). Each cohesive crack element is located between a pair of adjacent triangular elements. The elastic deformation is simulated by plane strain triangular elements, and the calculation method is based on the assumption that the triangular elements are in state of linear elastic and small deformation. The contact forces are calculated between all contacting couples that overlap in space. A distributed penalty function method and a Coulomb-type friction law of contact forces are, respectively, utilized to calculate the repulsive forces and the frictional forces of triangular elements. Three types of failure, including mode I (tensile), mode II (shear), and mixed mode III (tensile-shear) failure, can occur in the CCEs when CCEs reach the failure threshold states. After triangular elements are separated, they will change into discrete media undergoing finite displacement and rotation. The fracture criterion and constitutive behavior of cohesive crack elements in three failure modes are exhibited in Figure 1.
(a)
(b)
(c)
(d)
(e)
(f)
Because many relevant constitutive models and mechanical principles based on laboratory test results about crack initiation and growth are considered in the framework of FDEM, cracks grow freely under the constraint of the topological mesh based on the stress and strain condition. Therefore, the method can explicitly capture the continuous-discontinuous trans-scale fracture process and has been used in the analysis and simulation of the progressive fracture process of brittle rock. More detailed FDEM theory can be found in [19].
2.2. FDEM Simulation of Acoustic Emission
In FDEM, the elastic deformation is simulated by plane strain triangular elements, and the strain energy generated in the process of elastic deformation is stored in plane strain triangular elements when the deformation of brittle rock materials is simulated. The new cracks can occur when CCEs reach the failure threshold states, and the strain energy stored in plane strain triangular elements will release gradually. The released strain energy is dissipated in various forms of energies. The dissipated energies are mainly composed of fracture energy, frictional energy, and kinetic energy. The brittle rock can generate acoustic emission signals in the rupture process and the kinetic energy can be considered as acoustic emission event energy.
A cohesive crack element rupture can be considered as an acoustic emission event. Therefore, the mass center coordinates of failure CCEs can be regarded as the AE source location. The acoustic emission event energy is measured by kinetic energy increment of the ruptured cohesive crack element from entering the yield state to complete failure.
Because the time interval of the failure process in rock materials is very short, the CCEs also have this characteristic. The initiation time is determined via analyzing the kinetic energy change of the CCEs. The stress state and kinetic energy evolution of a cohesive crack element are typically exhibited in Figure 2. The time when the kinetic energy of the CCE reaches a peak is defined as AE initiation time .
The calculation algorithm of AE energy was developed for the particle flow code (PFC), a particle-based DEM software [20]. The detail calculation principles are as follows [21]:(1)When a cohesive crack element reaches the failure threshold state, the detail calculation equation of kinetic energy is as follows:where and are the nodal mass and nodal velocity at the time of yielding , respectively.(2)The kinetic energy is monitored by the change of CCE displacement (the element fails). The increment of kinetic energy is calculated at each time step as(3)The maximum increment of kinetic energy in the whole process of cohesive crack element failure from the yielding moment to the failing moment can be considered as the AE event energy . The detail calculation equation of AE event energy is as follows:(4)Gutenberg [22] proposed the calculation method of the AE event energy magnitude , the detail calculation equation of AE event energy magnitude is as follows:
2.3. Heterogeneous Model (HM) Generation
The Baihetan hydropower station [15] is located on the lower reaches of the main stream of the Jinsha River. The large-scale underground civil engineering, including the powerhouse cavern and diversion tunnel, is partially located in the rock mass of amygdaloidal basalt (as shown in Figure 3). Amygdales are common in this amygdaloidal basalt (as shown in Figure 4), and this heterogeneity caused by flaws significantly effects the strength characteristic and damage mechanism. Therefore, investigating the effect of amygdales on the fracture mechanism and damage evolution plays a great significant role in deeply revealing the damage evolution characteristic and predicting the risk level of rockburst in the amygdaloidal basalt surrounding the underground engineering.
(a)
(b)
(a)
(b)
In this work, digital image processing (DIP) was utilized to obtain the geometric outlines of amygdales (as shown in Figure 5(b)). In image processing, the image is processed by boundary enhancement filtering in Python. The number of ellipses was 51 and the total area of amygdales in the basalt was about 15.4%. The outlines of amygdales were approximately drawn using elliptic curves and import into the Gmsh environment via a C++ subroutine, and the heterogeneous model (HM) was wholly completed in Gmsh software. Finally, the HM was meshed into triangular elements via unstructured Delaunay algorithm embedded in Gmsh, to create the model compatible with Irazu based on FDEM.
(a)
(b)
(c)
To study the effect of ellipticity and orientations of amygdales on the mechanical characteristics, the inverse Monte-Carlo algorithm is used to generate elliptic curves with an artificial setting geometric feature in this study. Centroids of ellipses are generated by a random Poisson point process. More detailed inverse Monte-Carlo algorithm can be found in [23].
In FDEM, the sample changes into a two-dimensional rectangle, and the size of the rectangle is consistent with the laboratory test sample [14]. The size of the simplified two-dimensional model is 50 mm × 100 mm, and the height of the rectangle is 100 mm and the width of the rectangle is 50 mm. The two loading plates located on the upper and lower ends are rigid and are loaded opposite each other at a constant velocity of 0.1 m/s (as shown in Figure 6). Although the simulated loading rate is obviously greater than the loading rate commonly utilized in laboratory tests, the loading rate has been proven that, at the above loading rate, quasi-static conditions can be ensured [12]. The simulation adopts the plane strain model. The nodal displacement of loading plates and nodal force of the upper and lower ends of the rock sample need to be measured to calculate the uniaxial compression stress and axial strain, and the nodal displacement of the middle on both sides of the model needs to be measured to calculate the lateral strain.
2.4. Microparameters’ Verification
In the FDEM-HM, the boundaries between amygdales and matrix are simulated using the CCEs. This setting is able to reflect the difference in the microstrength characteristics of the boundary between two rock materials. Therefore, this method can reflect the heterogeneity caused by amygdales and capture spatial distribution characteristics of fracture evolution under the different elliptic coefficient and orientation conditions.
There are many mesoparameters in Irazu software when the simulation needs to be carried out, mainly including the basic physical and mechanical parameters of the rock materials (density, Young’s modulus, and Poisson’s ratio), the mesoparameters of the CCEs, and contact penalty values between elements. The laboratory test results can partially provide reference values of the mesoparameters. Many scholars have carried out lots of three-point bending tests to study the fracture toughness. The fracture toughness of brittle rocks is measured to calculate the mode I fracture energy based the theory method. The model II fracture energy is about ten times the mode I fracture energy in this study. The contact penalty values are generally set at 1–10 times larger than Young’s modulus. The simplified numerical model of amygdaloidal basalt is composed of matrix and amygdales regarded as homogenous and isotropic materials. The laboratory test results of cryptocrystalline basalt [15] provide reference values for the matrix. The strength of the amygdale flaws is lower than the matrix, and the mesoparameters of amygdales are appropriately weakened. This study utilizes a trial-and-error approach [24] based on the laboratory test results and weakened mesoparameter values. The input parameter values are regarded as reasonable until the FDEM numerical results are basically consistent with the laboratory test results. The mesoparameter values are exhibited in Table 1. The comparison of macroscopic mechanical parameters between the experiment [14] and simulation is shown in Table 2, including Young’s modulus, uniaxial compressive strength (UCS), Poisson’s ratio, and uniaxial tensile strength (UTS). The failure phenomena in the laboratory [25] and the simulation are shown in Figure 7.
(a)
(b)
(c)
(d)
3. Effect of Mesoheterogeneity on the Strength of Amygdaloidal Basalt
Bubeck et al. [26] thought that quantifying the shape, orientation, and statistical distribution of pores plays an important role in the laboratory experiments.
According to the characteristics of Baihetan amygdaloidal basalt, the following two numerical testing schemes of uniaxial compression are proposed. As shown in Table 3, the effect of mesoheterogeneity on strength and deformation are analyzed from the two aspects of ellipticity and orientation.
3.1. Effect of the Elliptic Coefficient of Amygdale
The amygdale content of each sample is about 15.4%, the number of amygdales in each sample is 60, and the orientation of amygdales in each sample obeys a normal distribution. The area of each amygdale is the same, and the elliptic coefficients of amygdales in the five samples are π, 3.2, 3.3, 3.4, and 3.5.
The relationships between the UCS, Young’s modulus, and elliptic coefficient are exhibited in Figure 8. With the increase of the elliptic coefficient, the UCS shows a linear growth trend and Young’s modulus shows a “V-shaped” characteristic distribution. The reason for this trend and distribution is that the stress field surrounding a circular hole is more uniform than that surrounding an elliptical hole in the process of compression, and oblique elliptical holes may also improve the mechanical properties of the sample [26].
Figure 9(a) shows the growth paths of cracks under different elliptic coefficient of amygdale conditions. Some researchers have investigated the rock materials containing filled holes. This study will deeply study the mechanical behavior incorporating previous research results. Li et al. [27] conducted many laboratory experiments and RFPA simulations about hydraulic fracturing of sandy conglomerate to investigate failure evolution and found that the new cracks usually propagated across the weakened gravels. The weakened gravels mainly underwent tensile failure. Heap et al. [25] found that cracks generally propagated between pores in andesite samples. The holes in rocks are mainly the best crack propagation paths, and the numerical results of this study verify this fact.
(a)
(b)
(c)
The AE evolution and AE events energy under different elliptic coefficients of amygdale conditions are exhibited in Figures 9(b) and 9(c). The zones where the AE events initiate are significantly different for the five samples. High-energy AE events are located near the amygdales. The shear cracks and energy release are closely related.
3.2. Effect of Amygdale Orientation
The amygdale content of each sample is about 15.4%, the orientation of amygdales in each sample remains the same, the number of amygdales in each sample is 60, and the elliptic coefficient of amygdales is 3.3. The orientations of amygdale in the five samples are 0°, 30°, 45°, 60°, and 90°.
The relationships between the UCS, Young’s modulus, and orientation are exhibited in Figure 10. As the orientation of amygdales increases, the UCS exhibits a “V-shaped” distribution characteristic, and Young’s modulus increases. Bubeck et al. [26] found that specimens that contain flat pores were weaker if compression was applied parallel to the short axis than if compression was applied parallel to the long axis. Yan et al. [28] found that the UCS of samples exhibited a “V-shaped” trend with increasing bedding angle in chlorite schist. The results of this paper have a certain degree of the engineering reference value.
Figure 11(a) shows the growth paths of cracks. The cracks generally propagated across the amygdales. The inclination angles of cracks change significantly as the amygdale orientation changes. When the orientations of the amygdales are 0° and 90°, the inclination angles of most fractures are parallel to the loading direction; when the orientation of the amygdales is 60°, the inclination angles of most cracks are along the orientation direction of the amygdales.
(a)
(b)
(c)
The AE evolution and AE events’ energy under different orientation of amygdale conditions are exhibited in Figures 11(b) and 11(c). The zones where the AE events initiate are significantly different for the five samples. AE events initiate in the amygdales, and the growth paths tend to be oriented toward the amygdales. High-energy AE events are located near the amygdales. The shear cracks and energy release are closely related.
4. Conclusions
This paper studies the effect of ellipticity and orientations of amygdales on the strength and deformation characteristics of amygdaloidal basalt. The main conclusions are summarized as follows:(1)Macroscopic fractures are caused by extension and coalescence of local cracks. The sample are mainly penetrated by some longitudinal macroscopic cracks, resulting in the flaking away of debris. The fracture mechanism of amygdaloidal basalt is splitting mechanism caused by tensile cracks, and the phenomena in the laboratory tests verify this fact. Simultaneously, the tensile cracks mainly occur in the matrix, and the shear cracks almost occur in the amygdales.(2)As the elliptic coefficient of amygdale increased, the UCS showed a linear growth trend, and Young’s modulus showed a “V-shaped” characteristic distribution. Weak holes provided the best crack propagation path.(3)As the amygdale orientation increased, the UCS exhibited a “V-shaped” characteristic distribution and Young’s modulus increased. The cracks generally propagated across the amygdales and the inclination angles of the cracks changed significantly with the orientation of the amygdales.
Data Availability
The data used to support the findings of this study were included in the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the National Natural Science Foundation of China (Grant no. 41877256), the Key Research Program of the Chinese Academy of Sciences (KFZD-SW-423), and the Natural Science Foundation of Hubei Province, China (Grant no. 2019CFB268).