Next Article in Journal
Parameter Optimization and Experimental Study on Tool-Vibration-Assisted Pulsed Electrochemical Machining of γ-TiAl TNM Blades
Next Article in Special Issue
A Finite Element Model for Investigating Unsteady-State Temperature Distribution and Thermomechanical Behavior of Underground Energy Piles
Previous Article in Journal
Optimization of a Nature-Inspired Shape for a Vertical Axis Wind Turbine through a Numerical Model and an Artificial Neural Network
Previous Article in Special Issue
Evaluating the Influence of Fracture Roughness and Tortuosity on Fluid Seepage Based on Fluid Seepage Experiments
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Simulation of the Fracturing Process of Inclusions Embedded in Rock Matrix under Compression

1
School of Civil Engineering, Inner Mongolia University of Science and Technology, Baotou 014010, China
2
Inner Mongolia Autonomous Region Key Laboratory of Civil Engineering Safety and Durability, Baotou 014010, China
3
Department of Civil and Environmental Engineering, Brunel University London, London UB8 3PH, UK
4
State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China
5
Zhejiang Engineering Survey and Design Institute Group Co., Ltd., Ningbo 315012, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2022, 12(16), 8041; https://doi.org/10.3390/app12168041
Submission received: 23 June 2022 / Revised: 4 August 2022 / Accepted: 5 August 2022 / Published: 11 August 2022
(This article belongs to the Special Issue Fracture and Failure of Jointed Rock Mass)

Abstract

:
Typical parallel fractures are often observed in the outcrops of inclusions in the field. To reveal the failure mechanism of inclusions embedded in rock matrix, a series of heterogeneous models are established and tested based on the damage mechanics, statistical strength theory, and continuum mechanics. The results show that, with the spacing between two adjacent fractures decreasing, the stress is firstly transferred from negative to positive, then from positive to negative. Stress transition is profound for the fracture spacing. Meanwhile, three types of fractures, i.e., consecutive fracture, non-consecutive fracture, and debonding fracture, are found, which are consistent with the observed modes in the field. Multiple inclusions are often fractured easier than an isolated inclusion due to the stress disturbance between inclusions and newly generated fractures. Either in single or multiple inclusions, tensile stresses inside the inclusions are the main driving force for fracture initiation and propagation. Besides, although the material heterogeneity has a small effect on the stress variation, it has an evident impact on the fracturing mode of inclusions. The stiffness ratio is critical for the stress transition and failure pattern; the interface debonding occurs earlier than the fracture initiation inside the inclusion when the stiffness ratio is relatively high. Additionally, the inclusions content only affects the sequence of fracture initiation rather than the final fracture spacing pattern.

1. Introduction

Typical spacing fractures in inclusions, accompanied by the interface debonding, can be directly observed in the outcrops in the field. For instance, thousands of stiff inclusions with grotesque shapes and sizes protruding from the soft sandstone matrix are observed in Yehliu GeoPark located along the northern coastal line of Taiwan and is famous for a rare and stunning geological landscape. Meanwhile, several fractures may be sub-parallel to each other in a single inclusion, as shown in Figure 1a,b. Eidelman et al. (1992) [1] reported that the well-organized tensile fractures in pebbles of young conglomerates were found in the Salton trough, California, and in the Dead Sea rift, Israel (Figure 1c). McConaughy et al. (1999) [2] observed the Squaw Point outcrop in the Devonian black shale near Seneca Lake, NY, and pointed out that the in-situ stress configuration may contribute to the interaction between joints and embedded inclusions. Bessinger et al. (2003) [3] also observed the area along a beach on southeastern Vancouver Island, British Columbia, which is characterized by calcareous inclusions with a variety of shapes and sizes embedded in sandstone matrix. Both interface debonding between inclusion and matrix and multiple closely parallel fractures inside inclusions can be found in the field. Simultaneously, the following common characteristics are found: (a) Inclusion outcrops are often stiffer than the surrounding matrix; (b) the shape of most inclusions is ellipsoidal or spherical. The size varies significantly from centimeters to meters; (c) the fractures in a single inclusion are closely parallel to each other. The fracture density inside inclusions is higher than elsewhere in the matrix; (d) the fracturing pattern can be considered as the indicator of the regional tectonic stress.
Some studies have been conducted to reveal the fracturing mechanism of inclusions embedded in the matrix. Häfner et al. (2003) [4] created a geometrical inclusion-matrix model for the finite element analysis of concrete at multiple scales. Crack initiation and propagation depend on the maximum size, shape, and size distribution of natural aggregates within the mortar matrix. Janeiro et al. (2010) [5] conducted a series of experiments considering the effects of inclusion shape, size, stiffness, and strength on the cracking behavior of materials containing inclusions. Their experimental models include both single inclusion and double inclusions. However, most studies mainly focused on the fracturing behavior of the surrounding matrix rather than the embedded inclusion.
Some researchers have investigated the stress solution of an inclusion in solid. Eshelby (1957) [6] determined the stress and strain of an ellipsoidal inclusion in the solid with the help of a sequence of imaginary cutting, straining, and welding operations. He believed that the stresses inside a circular or an elliptical stiff inclusion were uniform. Jaeger and Cook (1979) [7] expressed an analytical solution to the stress in an inclusion. Quesada et al. (2009) [8] developed the brittle fracture model with a mixed criterion involving both energy and stress conditions for multiple failures in the inclusion to predict such fractures. He suggested that the material properties and the inclusion size play a significant role in the fracture events. Meng et al. (2012, 2014) [9,10] developed a MATLAB code to evaluate the Eshelby solution for the ellipsoidal inclusion problem and applied it to model compaction and shear-enhanced compaction bands.
In this study, to investigate the fracturing process in inclusions embedded in the matrix, the rock failure process analysis (RFPA) [11,12,13] is employed to reveal the failure mechanism associated with the inclusions in the matrix. The established RFPA models can simulate the progressive failure process of quasi-brittle materials without any assumptions on the locations where new cracks will initiate and the paths along which cracks will propagate [14,15,16]. The effectiveness of the RPPA code has been verified by the benchmarks [17,18,19]. To model the failure of rock material (or rock mass), the rock medium is assumed to be composed of many mesoscopic elements whose material properties are different from one to another and are specified according to a Weibull distribution [20,21]. A more detailed introduction of the code, coupled with its applications in simulating geological phenomena, has been illustrated in previous studies [22,23,24,25,26].
The main objectives of this study include (a) to give an intuitive interpretation of the fracturing process of single inclusion and cluster of inclusions under compression; (b) to provide supplementary information on the stress distribution and fracturing-induced stress re-distribution in the vicinity of inclusions that cannot be observed directly in-situ to reveal the spacing fracture mechanism of inclusions; and (c) to illustrate the influence of material properties (stiffness and heterogeneity) on the fracturing pattern of single inclusion and cluster of inclusions.

2. Fracturing Behaviors of Single Inclusion Embedded in Rock Matrix

2.1. Stress Distribution Inside an Inclusion

Before investigating the fracture spacing behavior in an isolated inclusion, the stress distribution is firstly illustrated by using the RFPA code. Figure 2 is the employed numerical model in this study, in which an isolated inclusion is embedded in the rock matrix. The geometry of the whole model is 1200 mm × 1200 mm and the diameter of the inclusion is 200 mm. The specimen is meshed into 300 × 300 = 90,000 mesoscopic elements. The vertical compressive far-field stresses are applied to the top boundary. The model bottom is fixed in the Y-direction. The mechanical parameters for the rock sample are listed in Table 1. Eidelman and Reches (1992) [1] have mentioned that the Poisson ratio υi of the inclusion has a small effect on the stress solution, whereas the Poisson ratio υr of the matrix has a profound effect. Thus, the Poisson ratio υi of the inclusion is constant, i.e., 0.25 in the study. In this section, the homogeneity index, m, is chosen to be 100.0, representing the assumption that the rock sample is a homogeneous material. The influence of material heterogeneity on fracturing behavior will be discussed in detail in the following sections.
Jaeger and Cook (1979) [7] proposed that the analytical solution of stress inside an inclusion under the constraints can be expressed as follows:
σ ( max ) i = σ ( max ) r [ K ( X r + 2 ) + X i ] K [ X r + 1 ] 2 [ 2 K + X i 1 ] [ K X r + 1 ] + σ ( min ) r [ K ( X r 2 ) X i + 2 ] K [ X r + 1 ] 2 [ 2 K + X i 1 ] [ K X r + 1 ]
σ ( min ) i = σ ( max ) r [ K ( X r 2 ) X i + 2 ] K [ X r + 1 ] 2 [ 2 K + X i 1 ] [ K X r + 1 ] + σ ( min ) r [ K ( X r + 2 ) + X i ] K [ X r + 1 ] 2 [ 2 K + X i 1 ] [ K X r + 1 ]
K = μ i / μ r ,   X r = 3 4 υ r ,   X i = 3 4 υ i ,  
where σ(max)r and σ(min)r are the maximum and minimum compressive remote stress, respectively. σ(max)i and σ(min)i are the maximum and minimum compressive stress inside the inclusion, as shown in Figure 2. μi and μr are the shear moduli, and υi and υr are the Poisson ratios of the inclusion and the matrix, respectively.
Figure 3 shows the analytical result and numerical result of stress varying with the Poisson ratio of the matrix. Figure 3a is the stress variation for the case under uniaxial compressive loading, i.e., σ(min)r = 0. Figure 3b is the stress variation for the case under biaxial compressive loading, i.e., σ(min)r = 0.002 MPa/step, σ(max)r = 0.02 MPa/step. It is shown that the numerical solutions agree well with the analytical solutions in both cases. The value of σ(min)i/σ(max)r inside the inclusion is changed from positive to negative with the increase in the Poisson ratio υr of the matrix. It means the minimum stress in the inclusion transforms from compression into tension.

2.2. Stress Field between Two Adjacent Fractures

Like Figure 2, a similar model under uniaxial compression conditions is employed. The inclusion contains two pre-existing fractures. The mechanical parameters of the model are the same as those in Table 1. The Poisson ratio of the inclusion is 0.25 and the Poisson ratio of the matrix is 0.45 in bold. With the defined model, the minimum principal stress is examined as a function of the spacing between the two adjacent fractures.
As shown in Figure 4, with the decrease in the spacing between the adjacent fractures, the stress ratio, σ(min)r/σ(max)r, is firstly transferred from negative to positive, then transferred from positive to negative. In other words, the state of minimum principal stress, σ(min)i, changes twice, i.e., from tensile stress to compressive stress, and then from compressive stress to tensile stress. The results obtained here are slightly different from those about the spacing fractures in multiple-layered rocks [27,28,29,30,31]. In multiple-layer rocks, when the fracture spacing to layer thickness ratio reaches the critical value, the stress state transition from tensile to compressive, which defines the condition of fracture saturation. However, in the model employed in this study, when the spacing is relatively large (for instance, 200 mm in Figure 4a) or even if there are no fractures in the inclusion, the minimum principal stress is tensile. When the spacing is small enough (for instance, 40 mm in Figure 4e), the stress in the inclusion is still tensile, but the magnitude of the stress is so small that there are no more fractures that can be infilled between the adjacent two fractures.

2.3. Failure Analysis of Single Inclusion under Uniaxial Compression

In this section, a numerical model similar to Figure 2 is employed to investigate the fracturing behavior in a single inclusion subjected to uniaxial compression. The homogeneity index m of inclusion is 3.0, representing the assumption that the inclusion is a heterogeneous material. The Poisson ratio of the inclusion is 0.25 and the Poisson ratio of the matrix is 0.45. All other mechanical parameters are the same as those listed in Table 1. The fracture-induced AE energy accumulation under compression is shown in Figure 5.
Figure 6 shows the sequence of fracture infilling during the failure of an inclusion modeled in Figure 2. Figure 6a shows the acoustic emission (AE) distribution with red and white circles indicating tensile and compressive failure of the elements at the current loading step, while the black circles are used to indicate the damaged elements at previous loading steps. The AE events correlate well with the number of failed elements. The failure process is accompanied by a rapid increase in the AE event count [32,33,34]. Figure 6b shows the evolution of the fracture and maximum principal stress field. The grey scale indicates the relative magnitude of the stress within the elements. The elements with a lighter shade of grey have relatively higher stresses. The dark elements in Figure 6b represent the nucleated flaw. Fractures form through the connection of flaws.
According to the sequence of the fractures in Figure 6, the failure process of inclusion can be divided into three stages. In the first stage, the initial flaws are obviously initiated at the interface between the inclusion and matrix. In step 49, the macroscopic fracture starts at the interface, and several isolated fractures initiate randomly at the location inside the inclusion where the local tensile stress reaches its local tensile strength depending on the high heterogeneity of the inclusion. The AE events increase slightly because only a small number of elements with lower strength meet the failure criterion. In the second stage, from step 50 to step 80, the isolated fractures at the interface soon coalesce and propagate to form multiple sub-parallel fractures. There are sudden and rapid increases in the number of failure events. Each sharp increase in the number of AE events corresponds to the formation of a new macroscopic fracture. Almost all the fractures show tortuous paths which depend on the distribution of the heterogeneity in the material. As new fractures are continuously infilled between the earlier formed fractures, the spacing between adjacent fractures decreases gradually. At the last stage, as shown in step 100, the spacing between adjacent fractures is so closely that no more new fractures are infilled. In this stage, the AE events generally remain small and stable until new fractures nucleate and propagate.
In Figure 7, three types of fracture, i.e., consecutive fracture, non-consecutive fracture, and debonding fracture, can be observed. The first type of fracture is that two cracks initiate from the interface, then propagate and meet each other to form a macroscopic consecutive fracture. The macroscopic consecutive fracture completely cut through the inclusion. The non-consecutive fractures include the fractures initiated from the interface and the fractures initiated inside the inclusion. These non-consecutive fractures did not completely pass through the inclusion body. Parts of the non-consecutive fractures can develop into consecutive fractures. The debonding fracture is the cracks formed around the interface. Although the interface between the inclusion and matrix is assumed to be perfectly bonded, the delamination effect may occur for the obvious change of the local stress near the inclusion boundary [23]. Although the path of macroscopic fractures is almost flexuous with small flaws such as branches because of the heterogeneity of the material, the numerical results have a good conformity with the field observations shown in Figure 7.
The variation of minimum principal stress and the fractures location along the line A–A’ is illustrated graphically in Figure 8. Those fractures formed in the previous loading step are shown with solid black lines and each newly infilled fracture at the current step is shown with a solid red line.
When the far-field compression load is small at step 49, the stress of each element along the line A–A’ is tensile as shown in Figure 8a. At step 50, an interface debonding fracture (fracture No.1) first occurs in the inclusion. Once a fracture occurs, the stress is transferred and stress concentration is induced near the adjacent region of the previously formed fracture, as shown in Figure 8b. Consequently, a consecutive fracture (fracture No.2) and a non-consecutive fracture (fracture No.3) appear near fracture No.1 where the local tensile stress reaches its peak strength. At step 60, a non-consecutive fracture No.6 occurs between No.4 and No.5 at the point of tensile stress concentration, see Figure 8c.
The stress transition from tensile stress to compressive stress can be seen in the zone between fractures No.3 and No.5. Moreover, in Figure 8d, the zone transforms into a compression stress zone between fractures No.3 and No.5, and the tensile stress between fractures No.6 and No.10 changes to be smaller enough and a new fracture does not always initiate at the middle point between two previously generated fractures. The high-stress concentration zone heralds new fractures there between fractures No.7 and No.8 and between No.9 and No.4, which can be seen in Figure 8e such as fractures No.11, No.12, and No.13. However, with the increase in compression from step 100 to step 115, there are no new fractures forming in the inclusion and the spacing is so small enough that the stress between adjacent fractures is almost tensile such as between No.12 and No.8, No.8 and No.9. Even though the compression stress in the zone between fractures No.13 and No.11 and the zone between fractures No.11 and No.4, it is too small to induce new fractures further.

2.4. Effect of Material Properties on the Fracturing Behaviors of Inclusion

In this section, the effect of material properties of inclusion and matrix and the heterogeneity on the fracturing behavior of the inclusion will be discussed. The employed numerical model is still the same as Figure 2. The mechanical properties are shown in Table 2.

2.4.1. Effect of the Material Heterogeneity on Fracturing Behavior

In this study, homogeneity indexes of 1.5, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, and 20.0 were selected. The materials with higher m values represent more homogeneous materials, whereas those with lower m values are more heterogeneous. It can be seen in Figure 9 that the effect of heterogeneity on the stress inside inclusion is negligible. The absolute value of the stress ratio of σ(min)i/σ(max)r is close to 0.30, which agrees well with the analytical solution as shown in Figure 3.
However, the impact of heterogeneity on fracturing mode is not negligible. Figure 10 shows the failure process of inclusion with a homogeneity index of 1.5, 4.0, and 10.0, respectively. For the heterogeneous inclusion (e.g., m = 1.5), the location of initial failure is disordered and scattered. The final fracturing path is very flexuous and too many flaws are formed inside the inclusion. However, for the relatively homogeneous rock (e.g., m = 10.0), the fractures initiated from one side of the circular interface are kept propagating along the direction of maximum principal stress. The final fracturing path is relatively smooth. The number of flaws and consecutive fractures in the inclusion is relatively small.

2.4.2. Effect of Stiffness on Fracturing Behavior

Figure 11 shows the fractures development process of inclusion with different Ei. The stiffness of matrix (Er = 5 GPa) and other parameters remain the same, for a relatively soft inclusion, Ei = 5 GPa, the fractures occur firstly inside the inclusion, then the fractures are kept expanding to reach the inclusion interface along the direction of maximum principal stress. For a relatively hard inclusion embedded in a soft matrix, Ei = 50 GPa, most of the flaws are initiated in the interface at an earlier stage. There are several sub-parallel fractures in the inclusion extending from the interface. Finally, too many non-consecutive fractures are formed in the inclusion, because there are too many flaws initiated from the interface.
Figure 12 shows the different fracturing modes of inclusion with different matrix stiffness Er. In the case with a relatively soft surrounding matrix, Ei = 100 GPa, Er = 2 GPa, the debonding fractures firstly formed at an early stage. Then some fractures initiated from the interface are kept propagating to connect each other to be non-consecutive and/or consecutive fractures. While for the inclusion embedded in a relatively hard matrix, the debonding fractures along the interface are not easier to form. Most of the fractures are initiated from the central vicinity of the inclusion, then are kept extending along the direction of maximum principal stress. The number of final macroscopic fractures is relatively lower.

3. Fracturing Behavior of Multiple Inclusions

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

3.1. Failure Process of Multiple Inclusions

In rock strata in the field, there are always a large number of embedded inclusions with different sizes, shapes, and orientations. In this section, a 2D numerical model shown in Figure 13 is employed to investigate the fracturing behavior of multiple inclusions. In this model, a number of circular inclusions with different sizes are randomly distributed in the matrix. The model geometry and the boundary conditions are similar to that in Figure 2. The axial load is increased by 0.01MPa/step and the other mechanical parameters are the same as those listed in Table 3. The inclusion is assumed to be heterogeneous with the homogeneity index of m = 3.0, while the matrix is relatively homogenous, m = 20.0.
Figure 14a illustrates the fracturing process in detail, with the acoustic emission (AE) distribution. Moreover, Figure 14b is the fracturing process presented with maximum principal stress. Due to stress concentration, the fractures often initiate from the interface or near the cluster area with several inclusions adjacent to each other. With the increase in load, the fractures are kept expanding in inclusions, and more and more inclusions are fractured. Finally, almost all of the inclusions are fractured with different extents of damage.
The distribution of inclusions naturally formed in the matrix is random. It causes dense inclusion areas and sparse inclusion areas. Figure 15a,b shows the field observations and numerical simulation of inclusion fractures, respectively. It can be seen that the consecutive sub-parallel fractures mainly are in dense inclusion areas, while in sparse inclusion areas, non-consecutive fractures and debonding fractures can also be found. The numerical results agree well with the field observations as shown in Figure 15.

3.2. Effect of Stiffness on the Fracturing Behavior

Figure 16 is the numerical results with different inclusion stiffness. It is found that, when the Ei/Er ratio is lower than 1.0, the internal fractures and even the matrix failure are always the common modes, and few debonding fractures are found. Only a number of consecutive and non-consecutive fractures are formed in some inclusions, even though the compressive load is relatively high (the compressive load imposed on the top boundary of the model is 4.2 MPa at step 420). When the Ei/Er ratio is higher than 1.0, mixed fracturing modes are found, and consecutive sub-parallel fractures are easily formed with the increasing Ei/Er ratio. For example, in the case of Ei/Er ratio = 10.0, a large number of consecutive and non-consecutive fractures are formed when the compressive loading is 2.26 MPa at step 226.

3.3. Effect of Inclusion Content on the Fracturing Behavior

To analyze the fracturing behavior of the model with different content of inclusions, three values are selected, C = 15%, C = 30%, C = 45%. All the mechanical parameters and boundary conditions are kept constant as those in Figure 17.
Figure 17 is the corresponding numerical results. At the initial loading stage (for instance, at step 35), for the case containing more inclusions, the high inclusion content enhances the disturbing and interaction among inclusions. More scattered flaws are found in the case with inclusions content C = 45%. Nevertheless, the inclusion content has no evident impact on the final fracturing mode. Regardless of whether the inclusion content is low or high, the overall fracturing mode for the three cases is similar. In the final loading state, the three types of fractures, interface debonding fractures, consecutive fractures, and non-consecutive fractures, can be observed in the three cases.

4. Conclusions

In this study, the RFPA code was applied to conduct a series of numerical simulations to investigate the stress state and fractures evolution in single and multiple inclusions embedded in the rock matrix. The following conclusions can be drawn:
(1)
The spacing between two existing fractures directly governs the stress redistribution during fracture filling. With the decreasing of the spacing between the adjacent fractures, the stress is firstly transferred from negative to positive, then transferred from positive to negative. When the fracture spacing reaches the critical value, the tensile stress becomes so small that the magnitude of the stress is lower than the tensile strength of the inclusion. Consequently, there are no more fractures that can be infilled between the adjacent two fractures. Numerical results provide an intuitive way to see the fracturing process and fracture-induced stress re-distribution of inclusions embedded in the rock matrix, which cannot be observed directly in the field.
(2)
Three types of fractures may occur in either single or multiple inclusions, i.e., consecutive, non-consecutive, and debonding fractures. However, the path of macroscopic fractures may be flexuous with small flaws such as branches because of the heterogeneity of material. Meanwhile, although the effect of heterogeneity on the stress inside inclusion is negligible, the impact of heterogeneity on fracturing mode is not negligible. Fractures typically initiate at local concentrations of tensile stress around flaws. Since flaws produce greater stress concentrations, the location of fracture initiation depends on the distribution of the flaws as well as the magnitude of maximum principal tension stress.
(3)
The ratio Ei/Er affects the fracturing process of inclusions significantly. It is found that, when Ei/Er ratio is lower than 1.0, the internal fractures even the matrix failure are always the common mode, and few debonding fractures are found. Only a few numbers of consecutive and non-consecutive fractures are formed in some inclusions, even though the compressive load is relatively high. When the Ei/Er ratio is higher than 1.0, mixed fracturing modes are found. With the increasing Ei/Er ratio, a larger number of consecutive and non-consecutive fractures sub-parallel fractures are easily formed with a relatively lower compressive load.
(4)
The high inclusion content will enhance the disturbing and interaction between inclusions. Therefore, more scattered flaws can be found at the initial loading stage of the model containing multiple inclusions. Nevertheless, the inclusion content has no evident impact on the final fractured mode.

Author Contributions

Data curation, C.Y.; formal analysis, N.W.; funding acquisition, C.Y.; investigation, C.Y., P.X. and X.B.; methodology, B.G.; supervision, B.G.; validation, N.W.; writing—original draft, C.Y.; writing—review and editing, B.G. All authors have read and agreed to the published version of the manuscript.

Funding

This study was financially supported by the Inner Mongolia Autonomous Region Natural Science Foundation of China, China (Grant No. 2020BS05017) and Building Science Institute Open fund of Inner Mongolia University of Science and Technology, China (JYSJJ-2021Q06).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Eidelman, A.; Reches, Z.E. Fractured pebbles—A new stress indicator. Geology 1992, 20, 307–310. [Google Scholar] [CrossRef]
  2. McConaughy, D.T.; Engelder, T. Joint interaction with embedded concretions: Joint loading configurations inferred from propagation paths. J. Struct. Geol. 1999, 21, 1637–1652. [Google Scholar] [CrossRef]
  3. Bessinger, B.; Cook, N.G.W.; Myer, L.; Nakagawa, S.; Nihei, K.; Benito, P.; Suarez-Rivera, R. The role of compressive stresses in jointing on Vancouver Island, British Columbia. J. Struct. Geol. 2003, 25, 983–1000. [Google Scholar] [CrossRef]
  4. Häfne, S.; Eckardt, S.; Könke, C. A geometrical inclusion-matrix model for the finite element analysis of concrete at multiple scales. In Proceedings of the 16th IKM, Weimar, Germany, 10–12 June 2003. [Google Scholar]
  5. Janeiro, R.P.; Einstein, H.H. Experimental study of the cracking behavior of specimens containing inclusions (under uniaxial compression). Int. J. Fract. 2010, 164, 83–102. [Google Scholar] [CrossRef]
  6. Eshelby, J.D. The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems. Proc. R. Soc. Lond. 1957, 241, 376–396. [Google Scholar]
  7. Jaeger, J.C.; Cook, N.G.W. Fundamentals of rock mechanics. Third edition. Sci. Paperb. Vol. 1979, 9, 251–252. [Google Scholar]
  8. Quesada, D.; Leguillon, D.; Putot, C. Multiple failures in or around a stiff inclusion embedded in a soft matrix under a compressive loading. Eur. J. Mech. A Solids 2009, 28, 668–679. [Google Scholar] [CrossRef]
  9. Meng, C.; Pollard, D.D. Eshelby’s solution for ellipsoidal inhomogeneous inclusions with applications to compaction bands. J. Struct. Geol. 2014, 67, 1–19. [Google Scholar] [CrossRef]
  10. Meng, C.; Heltsley, W.; Pollard, D.D. Evaluation of the Eshelby solution for the ellipsoidal inclusion and heterogeneity. Comput. Geosci. 2012, 40, 40–48. [Google Scholar] [CrossRef]
  11. Tang, C. Numerical simulation of progressive rock failure and associated seismicity. Int. J. Rock Mech. Min. Sci. 1997, 34, 249–261. [Google Scholar] [CrossRef]
  12. Zhu, W.C.; Tang, C.A. Micromechanical Model for Simulating the Fracture Process of Rock. Rock Mech. Rock Eng. 2004, 37, 25–56. [Google Scholar] [CrossRef]
  13. Li, G.; Tang, C.A. A statistical meso-damage mechanical method for modeling trans-scale progressive failure process of rock. Int. J. Rock Mech. Min. Sci. 2015, 74, 133–150. [Google Scholar] [CrossRef]
  14. Wang, Y.Y.; Gong, B.; Tang, C.A.; Zhao, T. Numerical study on size effect and anisotropy of columnar jointed basalts under uniaxial compression. Bull. Eng. Geol. Environ. 2022, 81, 41. [Google Scholar] [CrossRef]
  15. Wang, Y.Y.; Gong, B.; Tang, C.A. Numerical investigation on fracture mechanisms and energy evolution characteristics of columnar jointed basalts with different model boundaries and confining pressures. Front. Earth Sci. 2021, 9, 763801. [Google Scholar] [CrossRef]
  16. Tang, S.B.; Tang, C.A. Crack propagation and coalescence in quasi-brittle materials at high temperatures. Eng. Fract. Mech. 2015, 134, 404–432. [Google Scholar] [CrossRef]
  17. Feng, X.H.; Gong, B.; Tang, C.A.; Zhao, T. Study on the non-linear deformation and failure characteristics of EPS concrete based on CT-scanned structure modelling and cloud computing. Eng. Fract. Mech. 2022, 261, 108214. [Google Scholar] [CrossRef]
  18. Chen, B.P.; Gong, B.; Wang, S.Y.; Tang, C.A. Research on zonal disintegration characteristics and failure mechanisms of deep tunnel in jointed rock mass with strength reduction method. Mathematics 2022, 10, 922. [Google Scholar] [CrossRef]
  19. Liang, Z.Z.; Gong, B.; Li, W. Instability analysis of a deep tunnel under triaxial loads using a three-dimensional numerical method with strength reduction method. Tunn. Undergr. Space Technol. 2019, 86, 51–62. [Google Scholar] [CrossRef]
  20. Weibull, W. The Phenomenon of Rupture in Solids; Generalstabens Litografiska Anstalts Förlag: Stockholm, Sweden, 1939. [Google Scholar]
  21. Basu, B.; Tiwari, D.; Kundu, D.; Prasad, R. Is Weibull distribution the most appropriate statistical strength distribution for brittle materials? Ceram. Int. 2009, 35, 237–246. [Google Scholar] [CrossRef]
  22. Li, L.C.; Tang, C.A.; Fu, Y.F. Influence of heterogeneity on fracture behavior in multi-layered materials subjected to thermo-mechanical loading. Comput. Mater. Sci. 2009, 46, 667–671. [Google Scholar] [CrossRef]
  23. Li, G.; Wang, K.; Gong, B.; Tao, Z.G.; Du, K. A multi-temporal series high-accuracy numerical manifold method for transient thermoelastic fracture problems. Int. J. Solids Struct. 2021, 230–231, 111151. [Google Scholar] [CrossRef]
  24. Gong, B.; Wang, S.Y.; Sloan, S.W.; Shen, D.C.; Tang, C.A. Modelling rock failure with a novel continuous to discontinuous method. Rock Mech. Rock Eng. 2019, 52, 3183–3195. [Google Scholar] [CrossRef]
  25. Gong, B.; Tang, C.A.; Wang, S.Y.; Bai, H.M.; Li, Y.C. Simulation of the nonlinear mechanical behaviors of jointed rock masses based on the improved discontinuous deformation and displacement method. Int. J. Rock Mech. Min. Sci. 2019, 122, 104076. [Google Scholar] [CrossRef]
  26. Liang, Z.Z.; Gong, B.; Wu, X.K.; Zhang, Y.B.; Tang, C.A. Influence of principal stresses on failure behavior of underground openings. Chin. J. Rock Mech. Eng. 2015, 34, 3176–3187. [Google Scholar]
  27. Bai, T.; Pollard, D.D. Fracture spacing in layered rocks: A new explanation based on the stress transition. J. Struct. Geol. 2000, 22, 43–57. [Google Scholar] [CrossRef]
  28. Bai, T.; Pollard, D.D.; Gao, H. Explanation for fracture spacing in layered materials. Nature 2000, 403, 753. [Google Scholar] [CrossRef] [PubMed]
  29. Bao, C.; Tang, C.A.; Cai, M.; Tang, S. Spacing and failure mechanism of edge fracture in two-layered materials. Int. J. Fract. 2012, 181, 241–255. [Google Scholar] [CrossRef]
  30. Li, L.C.; Tang, C.A.; Wang, S.Y. A Numerical Investigation of Fracture Infilling and Spacing in Layered Rocks Subjected to Hydro-Mechanical Loading. Rock Mech. Rock Eng. 2012, 45, 753–765. [Google Scholar] [CrossRef]
  31. Tang, C.A.; Liang, Z.Z.; Zhang, Y.B.; Chang, X.; Tao, X.; Wang, D.G.; Zhang, J.X.; Liu, J.S.; Zhu, W.C.; Elsworth, D. Fracture spacing in layered materials: A new explanation based on two-dimensional failure process modeling. Am. J. Sci. 2008, 308, 49–72. [Google Scholar] [CrossRef]
  32. Gong, B.; Wang, Y.Y.; Zhao, T.; Tang, C.A.; Yang, X.Y.; Chen, T.T. AE energy evolution during CJB fracture affected by rock heterogeneity and column irregularity under lateral pressure. Geomat. Nat. Hazards Risk 2022, 13, 877–907. [Google Scholar] [CrossRef]
  33. Wu, N.; Liang, Z.Z.; Zhang, Z.H.; Li, S.H.; Lang, Y.X. Development and verification of three-dimensional equivalent discrete fracture network modelling based on the finite element method. Eng. Geol. 2022, 306, 106759. [Google Scholar] [CrossRef]
  34. Tang, C.A.; Xu, X.H.; Kou, S.Q.; Lindqvist, P.A.; Liu, H.Y. Numerical investigation of particle breakage as applied to mechanical crushing—Part I: Single-particle breakage. Int. J. Rock Mech. Min. Sci. 2001, 38, 1147–1162. [Google Scholar] [CrossRef]
Figure 1. Examples of morphology of inclusions and the parallel fractures. (a) An outcrop of isolate inclusion with fractures (left) and its geologic sketch map (right), at Yehliu Geopark, Taiwan. (b) Several fractured inclusions embedded in islet at Yehliu Geopark, Taiwan. (c) Detailed map of fractured inclusions, Arava, Israel, photo by Eidelman and Reches (1992) [1].
Figure 1. Examples of morphology of inclusions and the parallel fractures. (a) An outcrop of isolate inclusion with fractures (left) and its geologic sketch map (right), at Yehliu Geopark, Taiwan. (b) Several fractured inclusions embedded in islet at Yehliu Geopark, Taiwan. (c) Detailed map of fractured inclusions, Arava, Israel, photo by Eidelman and Reches (1992) [1].
Applsci 12 08041 g001
Figure 2. Numerical model with a single circular inclusion embedded in matrix under compression conditions. σ(max)r is the load of 0.02 MPa/step applied in Y-direction of specimen. The yellow dotted rectangle (1/3 area of entire model) refers to the zone for making analysis of fracturing mode and stress distribution. The dotted red line A–A’ is the cross-section of the inclusion for the stress analysis in the following sections.
Figure 2. Numerical model with a single circular inclusion embedded in matrix under compression conditions. σ(max)r is the load of 0.02 MPa/step applied in Y-direction of specimen. The yellow dotted rectangle (1/3 area of entire model) refers to the zone for making analysis of fracturing mode and stress distribution. The dotted red line A–A’ is the cross-section of the inclusion for the stress analysis in the following sections.
Applsci 12 08041 g002
Figure 3. Analytical result and numerical result of normalized stress versus different Poisson ratio of matrix. (a) σ(min)r/σ(max)r = 0, υi = 0.25; (b) σ(min)r/σ(max)r = 0.1, υi = 0.25.
Figure 3. Analytical result and numerical result of normalized stress versus different Poisson ratio of matrix. (a) σ(min)r/σ(max)r = 0, υi = 0.25; (b) σ(min)r/σ(max)r = 0.1, υi = 0.25.
Applsci 12 08041 g003
Figure 4. The stress transition between two adjacent fractures with different spacing. (a) 200 mm. (b) 160 mm. (c) 120 mm. (d) 80 mm. (e) 40 mm. (Note: the left is the location of fractures along the diameter A–A’ crossing the inclusion and the right is the corresponding stress ratio along the cross-section A–A’ between two fractures.).
Figure 4. The stress transition between two adjacent fractures with different spacing. (a) 200 mm. (b) 160 mm. (c) 120 mm. (d) 80 mm. (e) 40 mm. (Note: the left is the location of fractures along the diameter A–A’ crossing the inclusion and the right is the corresponding stress ratio along the cross-section A–A’ between two fractures.).
Applsci 12 08041 g004
Figure 5. Fracture-event counts and AE-released energy accumulation versus remote compressive loading.
Figure 5. Fracture-event counts and AE-released energy accumulation versus remote compressive loading.
Applsci 12 08041 g005
Figure 6. Numerically simulated failure mode associated with both maximum stress field and AE distribution in the inclusion (a) AE events; (b) maximum principal stress.
Figure 6. Numerically simulated failure mode associated with both maximum stress field and AE distribution in the inclusion (a) AE events; (b) maximum principal stress.
Applsci 12 08041 g006aApplsci 12 08041 g006b
Figure 7. Field observation and numerical results of the fractured single inclusion. (a) An inclusion with closely spaced fractures, Vancouver Island, British Columbia (Bessinger et al., 2003 [3]). (b) One inclusion in the top left corner of Figure 1c [1]. (c) Numerical results of fractured single inclusion. The red wave lines indicate the debonding fractures, the green lines indicate the non-consecutive fractures and the blue lines indicate the consecutive fractures.
Figure 7. Field observation and numerical results of the fractured single inclusion. (a) An inclusion with closely spaced fractures, Vancouver Island, British Columbia (Bessinger et al., 2003 [3]). (b) One inclusion in the top left corner of Figure 1c [1]. (c) Numerical results of fractured single inclusion. The red wave lines indicate the debonding fractures, the green lines indicate the non-consecutive fractures and the blue lines indicate the consecutive fractures.
Applsci 12 08041 g007
Figure 8. Fracturing sequence and stress distribution inside inclusion. The right is the local map of minimum principal stress inside each element along the diameter A–A’. The left is the corresponding complete graph of the stress field. Each pre-existing fracture is shown by a solid black line and each infilling fracture at the current step is shown with a solid red line. (a) 0.98 MPa. (b) 1.0 MPa. (c) 1.2 MPa. (d) 1.6 MPa. (e) 2.3 MPa.
Figure 8. Fracturing sequence and stress distribution inside inclusion. The right is the local map of minimum principal stress inside each element along the diameter A–A’. The left is the corresponding complete graph of the stress field. Each pre-existing fracture is shown by a solid black line and each infilling fracture at the current step is shown with a solid red line. (a) 0.98 MPa. (b) 1.0 MPa. (c) 1.2 MPa. (d) 1.6 MPa. (e) 2.3 MPa.
Applsci 12 08041 g008
Figure 9. Effect of heterogeneity on the stress inside inclusion.
Figure 9. Effect of heterogeneity on the stress inside inclusion.
Applsci 12 08041 g009
Figure 10. Fracturing behavior of inclusion with different homogeneity index of inclusion.
Figure 10. Fracturing behavior of inclusion with different homogeneity index of inclusion.
Applsci 12 08041 g010aApplsci 12 08041 g010b
Figure 11. The effect of inclusion stiffness on the fracturing behavior inside inclusion. It is assumed the value of homogeneity index, m = 3.0 for the inclusion representing relatively heterogeneous materials and m = 20 for the matrix representing relatively homogeneous materials. Ei and Er are the elastic modulus of the inclusion and the matrix, respectively. In the case, Er = 5 GPa.
Figure 11. The effect of inclusion stiffness on the fracturing behavior inside inclusion. It is assumed the value of homogeneity index, m = 3.0 for the inclusion representing relatively heterogeneous materials and m = 20 for the matrix representing relatively homogeneous materials. Ei and Er are the elastic modulus of the inclusion and the matrix, respectively. In the case, Er = 5 GPa.
Applsci 12 08041 g011aApplsci 12 08041 g011b
Figure 12. The effect of matrix stiffness on the fracturing behavior inside inclusion. In the case, Ei = 100 GPa.
Figure 12. The effect of matrix stiffness on the fracturing behavior inside inclusion. In the case, Ei = 100 GPa.
Applsci 12 08041 g012
Figure 13. Sketch map of numerical model with multiple circular inclusions embedded in matrix under compressive condition established with RFPA. σ(max)r is the load of 0.01 MPa/step applied in Y-direction of specimens. The zone with a yellow dotted rectangle (3/4 area of entire model) is split out to make analysis of fracturing mode and stress distribution.
Figure 13. Sketch map of numerical model with multiple circular inclusions embedded in matrix under compressive condition established with RFPA. σ(max)r is the load of 0.01 MPa/step applied in Y-direction of specimens. The zone with a yellow dotted rectangle (3/4 area of entire model) is split out to make analysis of fracturing mode and stress distribution.
Applsci 12 08041 g013
Figure 14. Numerically simulated failure mode associated with both maximum principal stress field and AE distribution in the inclusions. (a) AE events; (b) maximum principal stress.
Figure 14. Numerically simulated failure mode associated with both maximum principal stress field and AE distribution in the inclusions. (a) AE events; (b) maximum principal stress.
Applsci 12 08041 g014
Figure 15. Comparison between field observations and numerical results. (a) The mapped section of the field site, Vancouver Island, British Columbia, and fracturing predominantly constrained to the inclusions. Debonding fractures and multiple, internal fractures within the inclusions are observed [3]. (b) Simulation results of multiple fractured inclusions at step 200, and the red wave line indicates the debonding fractures, the green wave line indicates the non-consecutive fractures and the blue line indicates the consecutive fractures.
Figure 15. Comparison between field observations and numerical results. (a) The mapped section of the field site, Vancouver Island, British Columbia, and fracturing predominantly constrained to the inclusions. Debonding fractures and multiple, internal fractures within the inclusions are observed [3]. (b) Simulation results of multiple fractured inclusions at step 200, and the red wave line indicates the debonding fractures, the green wave line indicates the non-consecutive fractures and the blue line indicates the consecutive fractures.
Applsci 12 08041 g015
Figure 16. Numerical results of fracturing behavior with different inclusion stiffness.
Figure 16. Numerical results of fracturing behavior with different inclusion stiffness.
Applsci 12 08041 g016aApplsci 12 08041 g016b
Figure 17. Numerical results of fracturing behavior with different inclusion content.
Figure 17. Numerical results of fracturing behavior with different inclusion content.
Applsci 12 08041 g017aApplsci 12 08041 g017b
Table 1. Physico-mechanical parameters employed in the model with single inclusion.
Table 1. Physico-mechanical parameters employed in the model with single inclusion.
ParametersInclusionMatrix
Homogeneity index (m)100100
Young’s modulus (E)I/GPa1005
Uniaxial compressive strength(fc0)/MPa50100
Tensile strength (ft0)/MPa2.510
Poisson ratio (ν)0.250.25, 0.3, 0.35, 0.4, 0.45, 0.5
Friction angle (φ)/°3030
Table 2. Physico-mechanical parameters employed in simulation for analysis of the effect of properties on the failure behavior of single inclusion.
Table 2. Physico-mechanical parameters employed in simulation for analysis of the effect of properties on the failure behavior of single inclusion.
ParametersInclusionMatrix
Homogeneity index (m)1.5, 2, 3, 4, 6, 8, 10, 15, 201.5, 2, 3, 4, 6, 8, 10, 15, 20
Young’s modulus (E)/GPa5, 50, 100, 1502, 5, 50, 150
Uniaxial compressive strength (fc0)/MPa50100
Tensile strength (ft0)/MPa2.510
Poisson’s ratio (ν)0.250.45
Friction angle (φ)/°3030
Table 3. Physico-mechanical parameters employed in the models with multiple inclusions.
Table 3. Physico-mechanical parameters employed in the models with multiple inclusions.
ParametersInclusionMatrix
Homogeneity index (m)320
Elastic modulus (E)/GPa3, 5, 20, 50, 1002, 5, 20, 100, 150
Uniaxial compressive strength (fc0)/MPa50100
Tensile strength (ft0)/MPa2.510
Poisson’s ratio (ν)0.250.45
Friction angle (φ)/°3030
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Yu, C.; Gong, B.; Wu, N.; Xu, P.; Bao, X. Simulation of the Fracturing Process of Inclusions Embedded in Rock Matrix under Compression. Appl. Sci. 2022, 12, 8041. https://doi.org/10.3390/app12168041

AMA Style

Yu C, Gong B, Wu N, Xu P, Bao X. Simulation of the Fracturing Process of Inclusions Embedded in Rock Matrix under Compression. Applied Sciences. 2022; 12(16):8041. https://doi.org/10.3390/app12168041

Chicago/Turabian Style

Yu, Chaoyun, Bin Gong, Na Wu, Penglei Xu, and Xiankai Bao. 2022. "Simulation of the Fracturing Process of Inclusions Embedded in Rock Matrix under Compression" Applied Sciences 12, no. 16: 8041. https://doi.org/10.3390/app12168041

APA Style

Yu, C., Gong, B., Wu, N., Xu, P., & Bao, X. (2022). Simulation of the Fracturing Process of Inclusions Embedded in Rock Matrix under Compression. Applied Sciences, 12(16), 8041. https://doi.org/10.3390/app12168041

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop