Study on the Effect of Amygdales on Mechanical Characteristics of Basalt Using the Combined Finite-Discrete Element Method

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 Irazu^2D 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.

[figures omitted; refer to PDF]

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 $T_i$.

[figure omitted; refer to PDF]

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 $E_{k,y}$ is as follows:$$\begin{array}{}\text{(1)}& E_{k,y}={\displaystyle \sum _{i=1}^4{m}_i{v}_{i,y}^2},\end{array}$$where $m_i$ and $v_{i,y}$ are the nodal mass and nodal velocity at the time of yielding $T=T_y$, 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$$\begin{array}{}\text{(2)}& \Delta E_{k}t=E_{k}t-E_{k,y}.\end{array}$$

(3) The maximum increment of kinetic energy in the whole process of cohesive crack element failure from the yielding moment $T_y$ to the failing moment $T_f$ can be considered as the AE event energy $E_e$ . The detail calculation equation of AE event energy is as follows:$$\begin{array}{}\text{(3)}& E_e=\underset{T_y,T_f}{\max }\Delta E_{k}t.\end{array}$$

(4) Gutenberg [22] proposed the calculation method of the AE event energy magnitude $M_e$, the detail calculation equation of AE event energy magnitude is as follows:$$\begin{array}{}\text{(4)}& M_e=\frac{2}{3}\lg E_e-4.8.\end{array}$$

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.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

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.

[figures omitted; refer to PDF]

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.

[figure omitted; refer to PDF]

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.

[figures omitted; refer to PDF]

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 omitted; refer to PDF]

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.

[figures omitted; refer to PDF]

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.

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).

References

[1] T. Kato, T. Nishioka, "Analysis of micro--macro material properties and mechanical effects of damaged material containing periodically distributed elliptical microcracks," International Journal of Fracture, vol. 131 no. 3, pp. 247-266, DOI: 10.1007/s10704-004-4558-z, 2005.

[2] H. T. Yang, H. Lin, Y. X. Wang, R. H. Cao, J. T. Li, Y. L. Zhao, "Investigation of the correlation between crack propagation process and the peak strength for the specimen containing a single pre-existing flaw made of rock-like material," Archives of Civil and Mechanical Engineering, vol. 21 no. 2,DOI: 10.1007/s43452-021-00175-w, 2021.

[3] S. J. Xie, H. Lin, Y. X. Wang, R. H. Cao, R. Yong, S. G. Du, J. T. Li, "Nonlinear shear constitutive model for peak shear-type joints based on improved Harris damage function," Archives of Civil and Mechanical Engineering, vol. 20 no. 3, pp. 89-101, DOI: 10.1007/s43452-020-00097-z, 2020.

[4] Y.-X. Xiao, X.-T. Feng, G.-L. Feng, H.-J. Liu, Q. Jiang, S.-L. Qiu, "Mechanism of evolution of stress-structure controlled collapse of surrounding rock in caverns: a case study from the Baihetan hydropower station in China," Tunnelling and Underground Space Technology, vol. 51, pp. 56-67, DOI: 10.1016/j.tust.2015.10.020, 2016.

[5] B. Li, T. Li, N. Xu, F. Dai, W. Chen, Y. Tan, "Stability assessment of the left bank slope of the baihetan hydropower station, southwest China," International Journal of Rock Mechanics and Mining Sciences, vol. 104, pp. 34-44, DOI: 10.1016/j.ijrmms.2018.02.016, 2018.

[6] J. A. Hudson, E. T. Brown, F. Rummel, "Controlled failure of rock diss and rings loaded in diametral compression," International Journal of Rock Mechanics and Mining Sciences, vol. 9 no. 3, pp. 241-248, DOI: 10.1016/0148-9062(72)90025-3, 1972.

[7] J. A. Hudson, "Tensile strength and the Ring test," International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts, vol. 6 no. 1, pp. 91-97, DOI: 10.1016/0148-9062(69)90029-1, 1969.

[8] A. Lisjak, G. Grasselli, "A review of discrete modeling techniques for fracturing processes in discontinuous rock masses," Journal of Rock Mechanics and Geotechnical Engineering, vol. 6 no. 4, pp. 301-314, DOI: 10.1016/j.jrmge.2013.12.007, 2014.

[9] Y. Zhang, L. N. Y. Wong, "A review of numerical techniques approaching microstructures of crystalline rocks," Computers & Geosciences, vol. 115, pp. 167-187, DOI: 10.1016/j.cageo.2018.03.012, 2018.

[10] A. Lisjak, G. Grasselli, "A review of discrete modeling techniques for fracturing processes in discontinuous rock masses," Journal of Rock Mechanics and Geotechnical Engineering, vol. 6 no. 4, pp. 301-314, DOI: 10.1016/j.jrmge.2013.12.007, 2014.

[11] A. Munjiza, D. R. J. Owen, N. Bicanic, "A combined finite‐discrete element method in transient dynamics of fracturing solids," Engineering Computations, vol. 12 no. 2, pp. 145-174, DOI: 10.1108/02644409510799532, 1995.

[12] O. Mahabadi, "Investigating the influence of micro-scale heterogeneity and microstructure on the failure and mechanical behaviour of geomaterials," ,DOI: 10.13140/2.1.2218.9762, 2012.

[13] M. Cai, "Fracture initiation and propagation in a Brazilian disc with a plane interface: a numerical study," Rock Mechanics and Rock Engineering, vol. 46 no. 2, pp. 289-302, DOI: 10.1007/s00603-012-0331-1, 2013.

[14] C. Q. Zhang, Z. J. Liu, Y. B. Pan, Y. Gao, H. Zhou, G. J Cui, "Influence of amygdale on crack evolution and failure behavior of basalt," Engineering Fracture Mechanics, vol. 226,DOI: 10.1016/j.engfracmech.2019.106843, 2020.

[15] Z. J Liu, C. Q. Zhang, C. S. Zhang, Y. Gao, H. Zhou, Z. R. Chang, "Deformation and failure characteristics and fracture evolution of cryptocrystalline basalt," Journal of Rock Mechanics and Geotechnical Engineering, vol. 11 no. 5, pp. 990-1003, DOI: 10.1016/j.jrmge.2019.04.005, 2019.

[16] X. Tan, H. Konietzky, W. Chen, "Numerical simulation of heterogeneous rock using discrete element model based on digital image processing," Rock Mechanics and Rock Engineering, vol. 49 no. 12, pp. 4957-4964, DOI: 10.1007/s00603-016-1030-0, 2016.

[17] X. Hu, N. Xie, Q. Z. Zhu, "Modeling damage evolution in heterogeneous granite using digital image-based grain-based model," Rock Mechanics and Rock Engineering, vol. 53, pp. 4925-4945, DOI: 10.1007/s00603-020-02191-3, 2020.

[18] S. R. Zhang, S. L. Qiu, P. F. Kou, S. J. Li, P. Li, S. Q. Yan, "Investigation of damage evolution in heterogeneous rock based on the grain-based finite-discrete element model," Materials, vol. 14 no. 14,DOI: 10.3390/ma14143969, 2021.

[19] A. Lisjak, O. Mahabadi, L. He, B. S. A. Tatone, P. Kaifosh, S. A. Haque, G. Grasselli, "Acceleration of a 2d/3d finite-discrete element code for geomechanical simulations using general purpose GPU computing," Computers and Geotechnics, vol. 100, pp. 84-96, DOI: 10.1016/j.compgeo.2018.04.011, 2018.

[20] J. F. Hazzard, R. P. Young, "Simulating acoustic emissions in bonded-particle models of rock," International Journal of Rock Mechanics and Mining Sciences, vol. 37, pp. 867-872, DOI: 10.1016/s1365-1609(00)00017-4, 2000.

[21] A. Lisjak, Q. Liu, Q. Zhao, O. K. Mahabadi, G. Grasselli, "Numerical simulation of acoustic emission in brittle rocks by two-dimensional finite-discrete element analysis," Geophysical Journal International, vol. 95, pp. 423-443, DOI: 10.1093/gji/ggt221, 2013.

[22] B. Gutenberg, "The energy of earthquakes," Journal of the Geological Society, vol. 112,DOI: 10.1144/gsl.jgs.1956.112.01-04.02, 1956.

[23] G. Mollon, J. Zhao, "“Fourier–Voronoi-based generation of realistic samples for discrete modelling of granular materials," Granular Matter, vol. 14, pp. 621-638, DOI: 10.1007/s10035-012-0356-x, 2012.

[24] B. S. A. Tatone, G. Grasselli, "A calibration procedure for two-dimensional laboratory-scale hybrid finite-discrete element simulations," International Journal of Rock Mechanics and Mining Sciences, vol. 75, pp. 56-72, DOI: 10.1016/j.ijrmms.2015.01.011, 2015.

[25] M. J. Heap, Y. Lavallee, L. Petrakova, P. Baud, T. Reuschle, N. R. Varley, D. B. Dingwell, "Microstructural controls on the physical and mechanical properties of edifice‐forming andesites at Volcan de Colima, Mexico," Journal of Geophysical Research: Solid Earth, vol. 119 no. 4, pp. 2925-2963, DOI: 10.1002/2013jb010521, 2014.

[26] A. Bubeck, R. J. Walker, D. Healy, M. Dobbs, D. A. Holwell, "Pore geometry as a control on rock strength," Earth and Planetary Science Letters, vol. 457, pp. 38-48, DOI: 10.1016/j.epsl.2016.09.050, 2017.

[27] G. Li, C. A. Tang, L. C. Li, "Three-dimensional micro flow-stress-damage (FSD) model and application in hydraulic fracturing in brittle and heterogeneous rocks," Key Engineering Materials, vol. 985, pp. 581-584, DOI: 10.4028/www.scientific.net/kem.452-453.581, 2011.

[28] S. Q. Yan, Q. K. Wang, H. F. Wang, S. L. Qiu, Z. Q. Zeng, Y. Fang, "Strength control factors of chlorite schist under schistose structure," International Journal of Design & Nature and Ecodynamics, vol. 15 no. 5, pp. 631-637, DOI: 10.18280/ijdne.150503, 2020.