Fundamental Study on Fracture Process Analysis of Rocks under Quasi-Static Loading Based on Hybrid FEM-DEM Using Extrinsic Cohesive Zone Model

Yutaro Maeda; Sho Ogata; Daisuke Fukuda; Toru Inui; Hideaki Yasuhara; Kiyoshi Kishida

doi:10.2320/matertrans.MT-Z2024005

Abstract

In this study, FDEM (ECZM), which is the combined finite-discrete element method (FDEM) based on an extrinsic cohesive zone model (ECZM), was improved by incorporating a novel algorithm that makes it easier to insert the cohesive elements than the conventional adaptive remeshing and a scheme to avoid the generation of dormant cohesive element, which is generated as a result of inserting a single cohesive element within the FDEM mesh. The improved FDEM (ECZM) was implemented in a self-developed 2-dimensional code, and it was first validated against a numerical experiment assuming a virtual rock. Subsequently, the proposed FDEM was applied to model the uniaxial compression and the Brazilian tests of siliceous mudstones under quasi-static loading. The results of the proposed FDEM simulations reasonably captured the trends of fracturing and stress-strain curves observed in the experiments well. The proposed scheme can also be extended to parallel computation based on general-purpose-graphic-processing-unit (GPGPU) and 3-dimensional FDEM (ECZM) with minimal efforts.

This Paper was Originally Published in Japanese in J. Soc. Mater. Sci., Japan 71 (2022) 206–213. References were modified due to duplications.

Outline of the proposed Master–Slave algorithm which is the adaptive separation process to avoid the adaptive remeshing in the extrinsic cohesive zone model.

1. Introduction

When constructing deep underground structures, such as tunnels, underground power generation plants, and high-level radioactive waste disposal facilities, fracture initiation and propagation may occur around the cavities excavated in the rock mass depending on various factors, namely, the in-situ stress conditions and geological conditions. Fracture initiation and propagation cause a decrease in the stiffness/strength and an increase in the permeability of the rock mass around these cavities [1]. Therefore, for the long-term safety of deep underground structures, one of the most important challenges is to reasonably model and predict fracture initiation and propagation during the excavation stage. Many numerical methods for fracture propagation in rock masses have been developed [2], and hybrid analysis methods, that combine continuum and discontinuum methods, have recently been proposed. Among these hybrid analysis methods, the Combined Finite-Discrete Element Method (FDEM) [3], which is a hybrid of the Finite Element Method (FEM) and the Discrete Element Method (DEM), has attracted attention for its applicability to issues related to the rock fracturing process, from relatively small-scale issues, such as lab-scale rock fracturing experiments, to large-scale issues, such as the excavation of cavities and blasting [2, 4–11]. FDEM can accurately trace the following important processes: (i) continuous deformation of the rock mass, (ii) transition from continuum to discontinuum, i.e., fracture initiation, propagation, and consolidation, and (iii) complex contact processes on rock surfaces, including newly generated fractures. As for (ii), FDEM describes the material softening process observed in the fracture process zone (FPZ) near the crack tip using the Cohesive Zone Model (CZM) [12, 13], in which a cohesive element, with the initial thickness of zero, is inserted at the boundary of two continuum elements. Cohesive traction, that resists fracture opening and sliding, is applied according to the open and slide of the cohesive element. CZM represents FPZ, which is a plastic region near the crack tip, by replacing it with a virtual single crack, called a smeared crack [3], on which the cohesive traction acts as shown in Figs. 1 and 2.

Fig. 1

Outline of Cohesive Zone Model.

Fig. 2

Outline of Cohesive element.

Furthermore, CZM is classified into the Intrinsic CZM (ICZM) [10, 14] and the Extrinsic CZM (ECZM) [10, 14] depending on the timing of the insertion of the cohesive elements. Most of the FDEM codes currently being developed are based on ICZM. In ICZM, the cohesive elements are inserted at the start of the simulation, making the insertion easier to implement and the calculation easier to parallelize. However, it has been pointed out that FDEM based on ICZM (hereinafter referred to as FDEM (ICZM)) causes a degradation in the accuracy of the calculation of the stress in the intact rock. This is because CZM was originally developed to model the FPZ near the crack tip, but FDEM (ICZM) must also artificially represent the elastic behavior before material softening [10]. In ECZM, on the other hand, the cohesive elements are not inserted at the start of the simulation; they are only inserted at the boundaries of the continuum elements that satisfy the prescribed criteria according to the states of the stresses and deformations in the continuum elements. Therefore, it is seen that the above-mentioned drawbacks of ICZM can be avoided, although its implementation is more complicated than that of ICZM, and several analyses using FDEM based on ECZM (hereinafter referred to as FDEM (ECZM)) have been reported [10, 15]. In addition, in the case of its extension to coupled problems related to heat and fluid, FDEM (ECZM) can calculate the heat transfer and fluid flow in intact rock with the same accuracy as the general FEM and can robustly evaluate the fracture initiation and propagation associated with heat expansion and seepage force. In contrast, in order to accurately calculate the heat conduction in intact rock within the framework of FDEM (ICZM), it is necessary to set the heat exchange coefficient which is experimentally difficult to evaluate for all cohesive elements. Previous studies have shown that the calculation accuracy varies largely depending on the value of the heat exchange coefficient [16]. Consequently, FDEM (ECZM) is expected to be a more useful analysis tool than FDEM (ICZM) for examining fracture initiation/propagation and associated coupled problems in rock masses.

However, the applicability of FDEM (ECZM) to fracture initiation and propagation issues in rock masses is not yet sufficiently understood, including the fundamental issues to be discussed below. The following is a list of items that require further study.

(1) Method of inserting cohesive elements
- - An algorithm without complicated adaptive remeshing should be developed for the insertion process.
- - When a single cohesive element is inserted in isolation, the geometrical constraints of the surrounding elements may restrict the displacement of the cohesive element, resulting in the non-physical suppression of fracture propagation, as shown in Fig. 3.
(2) Verification of the reproducibility of rock fracturing tests under quasi-static loading
- - Although the reproducibility of dynamic fracturing tests on rock has been verified in a previous study [10], verification cases of rock fracturing tests under quasi-static loading are deficient.

Fig. 3

Problem description when CE4 is isolated.

Based on the above background, a novel FDEM (ECZM) with a cohesive element insertion method, which does not require remeshing and does not suppress non-physical fracture propagation, is developed in this paper. Moreover, through simulations of rock fracturing experiments under quasi-static loading conditions, the reproducibility of the simulator is verified.

2. Treatment of CZM in FDEM

In FDEM, the above processes (i)–(iii) are described based on a continuous deformation model of a continuum element, CZM, that deals with the relative opening/closing of cohesive elements, the sliding of cohesive elements, and contact mechanics, respectively, and the equations of motion are solved for the nodes by the explicit time integration scheme. The only difference between FDEM (ICZM) and FDEM (ECZM) is the treatment of CZM. The purpose of this paper is to provide a fundamental study on the application of FDEM (ECZM) for process (ii). To this end, an explanation of FDEM (ICZM) is firstly necessary. In the following, the difference between FDEM (ICZM) and FDEM (ECZM) is explained. It is noted that this paper focuses on a two-dimensional (2D) problem for simplicity, but the findings from this study can easily be extended to three-dimensional (3D) problems. In addition, a plane strain condition is assumed in the present study.

2.1 Overview of FDEM (ICZM) [11]

In FDEM (ICZM) for 2D problems, a four-node cohesive element (hereinafter referred to as CE4) with zero initial thickness is inserted at the boundary of two three-node triangular elements (hereinafter referred to as TRI3s) that describe process (i). Normal cohesive traction σ^coh [N/m²] and tangential cohesive traction τ^coh [N/m²] (In this paper, tensile stress is positive.) are applied to the two relative surfaces composing CE4, according to the opening amount, o (The opening is positive.), and the sliding amount, |s|, of CE4, respectively, as shown in Fig. 2. σ^coh and τ^coh are modelled by the tensile softening law and shear softening law based on the Mohr-Coulomb shear strength model [11], respectively, as shown in Fig. 4. In ICZM, CE4s are inserted at the boundaries of all TRI3s at the start of the simulation, and the cohesive tractions, expressed by the following equations, are applied:

\begin{equation} \sigma^{\text{coh}} = \left\{ \begin{array}{cl} \dfrac{2o}{o_{\text{p}}}f_{\text{t}} & \text{if $o < 0$}\\ \left[\dfrac{2o}{o_{\text{p}}} - \left(\dfrac{o}{o_{\text{p}}} \right)^{2} \right]f(D)f_{\text{t}} & \text{if $0 \leq o \leq o_{\text{p}}$}\\ f(D)f_{\text{t}} & \text{if $o_{\text{p}} \leq o \leq o_{\text{t}}$} \end{array} \right., \end{equation}

(1)

\begin{equation} \tau^{\text{coh}} = \left\{ \begin{array}{cl} \left[\dfrac{2|s|}{s_{\text{p}}} - \left(\dfrac{|s|}{s_{\text{p}}} \right)^{2} \right]f_{\text{s}} & \text{if $0 \leq |s| \leq s_{\text{p}}$}\\ f_{\text{s}} & \text{if $s_{\text{p}} \leq |s| \leq s_{\text{t}}$} \end{array} \right., \end{equation}

(2)

\begin{equation} f_{\text{s}} = - \sigma^{\text{coh}}\tan \phi + f(D)c, \end{equation}

(3)

where o [m] and |s| [m] are the opening and sliding of the two surfaces of CE4, respectively, o_p [m] and s_p [m] are the artificial elastic limits of o and |s|, respectively, and f_t [N/m²], f_s [N/m²], c [N/m²], and ϕ [°] are the microscopic tensile strength, shear strength, cohesion, and internal friction angle, respectively. It is noted that the word “microscopic” means the scale of a grain in contrast to the scale of a specimen. D is the damage variable indicating the damage state of CE4 and f(D) is the softening function that determines the shape of the softening curves. o_p and s_p are defined using element size h [m] and penalty coefficient P_coh [N/m²], as given by o_p = 2hf_t/P_coh and s_p = 2hc/P_coh. Shape function f(D) and damage variable D (0 ≤ D ≤ 1) are defined by the following equations using tensile damage variable D_o and shear damage variable D_s:

\begin{equation} D_{\text{o}} = \min\left(1, \frac{o - o_{\text{p}}}{o_{\text{t}} - o_{\text{p}}} \right)\quad \text{if $o_{\text{p}} \leq o \leq o_{\text{t}}$, otherwise 0}, \end{equation}

(4)

\begin{equation} D_{\text{s}} = \min\left(1, \frac{|s| - s_{\text{p}}}{s_{\text{t}} - s_{\text{p}}} \right)\quad \text{if $s_{\text{p}} \leq |s| \leq s_{\text{t}}$, otherwise 0}, \end{equation}

(5)

\begin{align} f(D) & = \left[1 - \frac{\alpha + \beta - 1}{\alpha + \beta}\exp \left(D\frac{\alpha + \gamma \beta}{(\alpha + \beta)(1 - \alpha - \beta)} \right) \right]\\ &\quad \times [\alpha (1 - D) + \beta (1 - D)^{\gamma}], \end{align}

(6)

\begin{equation} D = \min\big(1,\sqrt{D_{\text{o}}{}^{2} + D_{\text{s}}{}^{2}}\big), \end{equation}

(7)

where α [-], β [-], and γ [-] are parameters that determine the curve shape; their values are set as in a previous study [11]. o_t [m] and s_t [m] are the critical opening and sliding, respectively, which separate microscopic and macroscopic fractures. They are uniquely determined from Mode I and Mode II fracture energy and the shape of f(D) [11].

Fig. 4

Tensile (left) and shear (right) softening curves. Proposed algorithm for the computation of cohesive traction when CE4 inserted in FDEM (ECZM).

2.2 Overview of FDEM (ECZM)

In FDEM (ECZM), CE4s are not inserted into the boundaries of TRI3s at the start of the simulation; therefore, the calculation when CE4s do not exist is equivalent to that of an explicit FEM analysis. Then, at each time t, CE4 is adaptively inserted only at the boundary of two TRI3s that satisfy the failure criteria. Two problems are discussed in terms of the previous study on FDEM (ECZM) [10]. As to the failure criteria for inserting CE4 into ECZM, the average value σ_ave of the Cauchy stress tensor of two TRI3s adjacent to the element boundary of TRI3s is calculated. Letting n be the outward unit normal vector of the surface on which σ_ave acts, the normal stress (σ_n = σ_aven) and the shear stress (τ_n = σ_ave − σ_aven⊗n) acting on the boundary are calculated (⊗ is the tensor product). When the normal stress reaches the microscopic tensile strength assigned to the boundaries or the shear stress reaches the microscopic shear strength based on the Mohr-Coulomb model, the failure criteria are assumed to have been satisfied, and CE4 is inserted. Here are the two problems with this approach. The first problem is that when a homogeneous microscopic strength distribution is given for all the element boundaries, CE4 is not inserted unless σ_n and τ_n have reached their respective designated levels of strength. In recent FDEM analyses, unstructured meshes are often used. This leads to heterogeneity being unintentionally included in the model even if no heterogeneity is introduced because σ_n and τ_n are dependent on n. In such cases, even if adopting a heterogeneous model, e.g., one in which the heterogeneous distribution of strength is given at the boundaries based on the Weibull distribution [17], the effect of introducing a heterogeneous model is lost. This is because the heterogeneous effect due to the above-mentioned unstructured meshes has already been included. In addition, the timing of the fracture initiation, i.e., the CE4 insertion, is strongly dependent on the mesh. The second problem is the complicated remeshing algorithm [18] that is applied for the CE4 insertion process. It leads to the inapplicability of the sequential calculation of this algorithm to the large-scale parallel computation proposed by the authors [11] due to the sequential characteristics of the remeshing algorithm, which is a bottleneck in the analyses. Therefore, this paper proposes a novel FDEM (ECZM) to solve the above two problems.

2.2.1 Failure criteria at element boundary

The failure criteria at the element boundary are based on the principal stresses, as expressed in the following equation. They are used instead of the criteria based on σ_n and τ_n described above.

\begin{equation} \left\{ \begin{array}{c} F_{1} \equiv \sigma_{1} - f_{\text{t}}\\ F_{2} \equiv - \sigma_{3} + \dfrac{1 + \sin \phi}{1 - \sin \phi}\sigma_{1} - 2c\cos \phi \end{array} \right., \end{equation}

(8)

where σ₁ [N/m²] and σ₃ [N/m²] are the maximum and minimum principal stresses, respectively, and CE4 is inserted when either of the two TRI3s across the element boundary satisfies F₁ > 0 or F₂ > 0. In this way, the mesh dependency of the timing of the CE4 insertion can be reduced. However, FDEM (ECZM) requires a seamless transition with the CE4 insertion from the boundary stresses (σ_n and τ_n) acting just before the CE4 insertion to the cohesive tractions (σ_coh and τ_coh). For this purpose, the following methods are adopted in this study. At first, it is important to note that, before the CE4 insertion, the relevant element boundary is not separated and the relative displacement of the element boundary surfaces is zero, and that, after the CE4 insertion, the relative displacement starts to appear at the boundary. Therefore, as shown in Fig. 4, the cohesive tractions are calculated using the sum of the actual relative displacements after the CE4 insertion (o^E and |s^E|) and the nominal displacements (o₀ [m] and |s₀| [m]) as the relative displacements (o and |s|) in eqs. (1) and (2). The nominal displacements are defined from boundary stresses σ_n (= σ₀) and τ_n (= τ₀) when eq. (8) is satisfied, i.e., when CE4 is inserted, according to the following equations:

\begin{equation} o_{0} = \min \left(o_{\text{p}},\frac{2h\sigma_{0}}{P_{\text{coh}}} \right), \end{equation}

(9)

\begin{equation} |s_{0}| = \min \left(s_{\text{p}},\frac{2h|\tau_{0}|}{P_{\text{coh}}}\right). \end{equation}

(10)

In this way, discontinuity does not occur between the boundary stresses and the cohesive tractions upon the CE4 insertion.

Finally, in previous FDEMs, both ICZM and ECZM, the critical opening and sliding (o_t and s_t) in the softening zone of CE4 were determined from Mode I and Mode II fracture energy and the shape of f(D), as described above [11]. In the present approach, the critical displacements of CE4 (o_t and s_t) are constant, i.e., independent of the element size, but the relative displacements (o, |s|) themselves are dependent on the element size, resulting in the fracture propagation behavior being dependent on the element size. In order to solve this dependency on the element size, apparent opening strain ε_open and shear strain ε_slip are calculated in this study from the relative displacements of CE4 (o, |s|), respectively, as shown in the following equations:

\begin{equation} \varepsilon_{\text{open}} = \frac{o - o_{\text{p}}}{h}, \end{equation}

(11)

\begin{equation} \varepsilon_{\text{slip}} = \frac{|s| - s_{\text{p}}}{h}. \end{equation}

(12)

Furthermore, using unknown parameter λ [-] and the Young’s modulus of the rock around CE4, E [N/m²], it is assumed that CE4 completely fails and a macroscopic fracture is formed when ε_open reaches λ(f_t/E) or ε_slip reaches λ(c/E). The value of λ is set to 100 based on a preliminary study. The establishment of a method for determining λ according to the rock type and scale of the target object is a subject for future work.

In order to take account of the heterogeneity of rocks, which are brittle materials, the Weibull distribution is used to statistically distribute the mechanical properties in the analysis domain [18]. Specifically, among the mechanical properties, only Young’s modulus E and each microscopic strength (tensile strength f_t and cohesion c) are distributed based on the Weibull distribution. E is given to TRI3, and f_t and c are given to the TRI3 boundary or CE4. Assuming that the mechanical properties of each calculation element are u (= E, f_t, c), the index values of the distribution of the mechanical properties are u^ref (= E^ref, f_t^ref, c^ref) and the shape parameter indicating the heterogeneity of the Weibull distribution is m. The distribution of the mechanical properties is given according to the probability density function shown in the following equation:

\begin{equation} \varphi (u) = \frac{m}{u^{\text{ref}}}\left(\frac{u}{m}\right)^{m - 1} \exp \left[- \left(\frac{u}{u^{\text{ref}}} \right)^{m} \right]. \end{equation}

(13)

2.2.2 Algorithm of CE4 insertion

Next, a new method is proposed to overcome the highly sequential process of the CE4 insertion, namely, the Master–Slave algorithm (hereinafter referred to as the M–S method) [19]. An overview of the process is shown in Fig. 5. In the M–S method, all the TRI3 element boundaries are virtually separated at the start of the simulation, similar to that for FDEM (ICZM), as shown in Figs. 5(a) and (b). Then, for each group of nodes separated from the common nodes before the virtual separation, one node is designated as the Master node (M-node) and the remaining nodes are designated as Slave nodes (S-nodes), as shown in Fig. 5(b). Focusing on a nodal group consisting of an M-node and its corresponding S-nodes (hereinafter referred to as a common nodal group), the nodal masses and all the nodal forces calculated for each TRI3 with S-nodes are gathered from each S-node to the M-node. In this way, all the nodal information is retained by the M-node before the insertion of CE4 for the group nodes. After that, the equations of motion are solved for the M-node only, and its nodal coordinates are updated. Then, the updated information on the M-node (i.e., coordinates and nodal velocity) is copied to the S-nodes belonging to it. As shown in Fig. 5(d), to determine the M-node, a vector connecting the current position of a common nodal group to the geometric center of each TRI3, which has the separated nodes that make up the common nodal group, is considered. The angle θ between this vector and the reference axis is calculated for all the separated nodes. The node with the smallest θ is designated as the M-node and the other nodes are designated as the S-nodes. In this approach, CE4 can be newly inserted very easily by updating the relation of the M-node and the S-nodes (hereinafter referred to as the M–S nodal relation) for each common nodal group, as described below. The following explanation uses the example of newly inserting two CE4s into the element boundaries containing the nodes located in the center of Fig. 5. In this case, each common nodal group is traced counterclockwise in the order of its associated θ value, and the separated node, immediately after crossing the CE4 boundary to be inserted, is registered as a new M-node. Then, all the separated nodes, before crossing the other CE4, are registered as S-nodes corresponding to the M-node, and the same process is repeated until one round is completed. As a result, the number of M nodes is equal to the number of inserted CE4s. In this approach, even when three or more CE4s exist in the area surrounding the nodes, no special conditional branching is required, and processing can be performed for each common nodal group.

Fig. 5

Outline of Master-Slave algorithm.

In addition, while it is necessary in the conventional adaptive remeshing algorithm [18] to perform special and extremely complicated processing, such as the connection of CE6 (a six-node cohesive element in the 3D case) by modifying the connectivity of the mesh, the proposed algorithm only requires the modification of the M-S nodal relation for each common nodal group, making the implementation extremely simple. It should be noted that the large-scale parallel computation itself is not the target of the discussion in this paper, and the above algorithm is implemented as a trial within the framework of the sequential computation. A drawback to this approach, compared to the CE4 insertion process with remeshing, is that the random-access memory (RAM) used here is expected to increase because the nodes are separated from the onset of the simulation, although virtually, as in ICZM. However, the amount of RAM used by FDEM, which is based on an explicit solver, is very small compared to that of implicit solver-based methods. Therefore, even considering the extra amount of RAM needed to adopt this approach, the total amount of required RAM is about the same as that of FDEM (ICZM). The approach is also innovative in that only a small amount of programming is required to add it to the FDEM (ICZM) code. It is important to note that this small amount of effort leads to the avoidance of the complicated remeshing process that has been the cause of the extremely limited development and reporting of FDEM (ECZM) in the past.

2.2.3 Problem with CE4 insertion

As mentioned above, in the FDEM (ECZM) based on the TRI3 mesh, a CE4 is inserted when a TRI3 boundary has satisfied the fracture criteria. However, as shown in Fig. 3, when a single CE4 is inserted at one of the element boundaries in isolation, this CE4 is geometrically unable to open/close or slide, and the effect of the failure of this cohesive element is not reflected even though the failure criteria have been locally satisfied. This problem is inevitable in FDEM (ECZM) using TRI3 in 2D and in 4-node tetrahedral elements in 3D, but was not addressed in the previous study [10]. A possible solution to this problem is to use higher-order elements with intermediate nodes on the element boundaries. However, this approach is undesirable in FDEM, where the element size must be small enough regardless of the element type, because it requires a large amount of computation and significant modifications to the existing FDEM codes. In this paper, therefore, as shown in Fig. 6 by the thick solid line, when a CE4 is inserted and no CE4 exists near it, additional CE4s are inserted at the element boundaries indicated by the dashed lines in the figure, even if they do not satisfy the failure criteria. For these additional surrounding CE4s, as discussed in Section 2.2.1, a seamless transition from the element boundary stress acting immediately before the CE4 insertion to the cohesive traction is applied. It is noted that since these CE4s do not satisfy the failure criteria, they are in the artificial elastic region (o < o_p, 0 < |s| < s_p) shown in Fig. 4, and no softening of the cohesive traction occurs. The reason for selecting all the surrounding element boundaries for the insertion of CE4s is not only because the insertion is easy to implement, but also because it is difficult to reasonably select the boundaries at which additional elements are to be inserted when applied to a 3D problem or when the grain distribution is considered in detail even for 2D problems.

Fig. 6

Schematic of approach inserting additional CE4s adjacent to element boundary where meets failure criteria.

3. Reproduction Analyses

In this chapter, the FDEM (ECZM) implemented in this study is used to perform analyses to reproduce the Brazilian tensile strength (BTS) test and uniaxial compression strength (UCS) test on siliceous mudstone to validate the simulator.

3.1 Overview of reproduction analyses

Reproduction analyses of the BTS and UCS tests were conducted using siliceous mudstone sampled from the depth of 350 m at the Horonobe Underground Research Laboratory of the Japan Atomic Energy Agency (JAEA). The rock specimens were 30 mm in diameter and 29.4 mm in depth for the BTS test, and 30 mm in diameter and 56.3 mm in height for the UCS test. A 2D model, representing the specimens used in each experiment, was created in these analyses, as shown in Fig. 7. The model was divided by an unstructured mesh consisting of a TRI3 with an average element size of 0.5 mm.

Fig. 7

Numerical models of a Brazilian test simulation (left) and Uniaxial compression test simulation (right).

For the analytical parameters, common values were set for both tests, as shown in Table 1. In a previous study [20], it was reported that the value of the shape parameter m of the Weibull distribution, that can reasonably reproduce fracture initiation and propagation, is between 1.2 and 5.0 when the Weibull distribution is used for fracture analyses considering the heterogeneity of the mechanical properties of rocks. Based on this, m was set to 3.0 in these analyses. The properties obtained from the experiments were a density of 1835 kg/m³, Young’s modulus of 1.601 GPa, Poisson’s ratio of 0.17, Brazilian tensile strength of 2.393 MPa, and uniaxial compression strength of 16.69 MPa. Young’s modulus was calculated as the slope of the stress value at two points, 45% and 55%, of the UCS test. The index values of the mechanical properties (E^ref, f_t^ref, and c^ref) used for the Weibull distribution were set to be similar between the experimental and simulation results. The density and Poisson’s ratio were set to be equal to the aforementioned experimental values. The microscopic internal friction angle was set to be 26° based on the results of triaxial compression tests using rock masses sampled at various depths similar to that (350 m) at the Horonobe Underground Research Laboratory [21]. As for the loading rate for the BTS and UCS tests, apparent strain rates of 1.2 × 10⁻⁶%/step and 6.2 × 10⁻⁷%/step, which are equal to the displacement velocity of 0.05 m/s (3.5 × 10⁻¹⁰ m/step), were applied to the upper and lower loading plates, respectively. In the explicit method, where time increment Δt is limited by the Courant condition, it is difficult to use the experimental strain rate of 0.1%/min, which is equal to the displacement velocity of 0.01 mm/s (7.0 × 10⁻¹⁴ m/step) in the case of the UCS test, from the viewpoint of the total analysis time. Therefore, in order to approximately model the quasi-static loading conditions in an explicit method, a critical artificial viscous damping scheme [11] is introduced in these analyses to reduce the undesirable wave effects, such as stress reflections, and the loading rate is sufficiently lowered. As an indication that a quasi-static state has been approximately achieved, the standard at which the response obtained from the analysis is considered to be quasi-static when the ratio of the total kinetic energy to the total strain energy in the rock specimen is less than 0.05 for the elastic behavior of the entire rock, is widely adopted [22, 23]. In the present simulations, the loading rate is considered to be sufficiently low because the aforementioned ratio of the energies is about 10⁻⁵, much smaller than 0.05 which is considered to be the standard for achieving the quasi-static condition. It should be noted that a quantitative evaluation of the fracture propagation speed is not possible due to the introduction of artificial viscous damping, which artificially reduces the stress wave behavior.

Table 1 Parameters used in FDEM simulations.

3.2 Analysis results and discussion

Firstly, the results of the reproduction analysis of the BTS test are explained. Figure 8 shows the indirect tensile stress–axial strain curve obtained from the numerical simulation with the experimental results for the tensile strength, using a dashed line. The figure shows that the peak strengths obtained from the simulation and the experiment are close to each other. Next, the distribution of damage variable D, which indicates the fracture propagation process, is shown in Fig. 9(a). This figure shows the fracture initiation and propagation process from Point A to Point B for the indirect tensile stress–axial strain curve in Fig. 8. Figure 9(a) shows that a microscopic fracture occurs near the upper and lower loading plates, and then a macroscopic fracture occurs in the center of the specimen in the axial direction, leading to failure of the specimen. Figure 9(b) shows the fracture modes for the inserted CE4 at Point B. The blue and red lines in the figure correspond to the CE4 inserted by satisfying the tensile and shear failure criteria, respectively. The figure shows that shear failure is dominant near the loading plates, while tensile failure is dominant in the center of the specimen, and the failure of the entire specimen is attributed to the propagation of a macroscopic fracture caused by the tensile failure. In addition, compared to the fracture patterns in the experiment shown in Fig. 9(c), although the number of generated fractures near the loading plates is slightly overestimated, the fracture patterns obtained from the simulation reproduce the experimental results well.

Fig. 8

Indirect tensile stress versus axial strain curve.

Fig. 9

Results of FDEM simulation for Brazilian test. (a) Fracture process, (b) Failure modes, (c) Fracture pattern in experiment.

The following section describes the results of the reproduction analysis of the UCS test. Firstly, a comparison of the results obtained from the simulation and the experiment for the axial stress–axial strain curve is given in Fig. 10. The figure shows that the experimental results deviate from the simulation results because the stiffness at the beginning of the experiment is small due to the nonlinear behavior of the axial stress–axial strain curve. In the experiment, plastic deformation, which is induced by the reloading of the specimens that were unloaded from the in-situ stress state, and the pre-existing fractures, tend to cause strain at the beginning of the experiment. However, they are not taken into account in the analyses, which may be the reason for the deviation. The simulation results, including the pre- and post-peak loading processes, show a relatively good agreement with the experimental results, especially in terms of the strength and strain at the peak and Young’s modulus. The simulation results show nonlinear behavior with a decrease in stiffness near the peak strength, due to the effect of the transition of a large number of CE4 from the pre-softening (non-softening) process to the softening process as a result of the microscopic fracture initiation and propagation.

Fig. 10

Axial stress versus axial strain curve.

Next, the distribution of damage variable D, which indicates the fracture propagation process, is shown in Fig. 11(a). The figure shows the fracture initiation and propagation process from Point A to Point B in the axial stress–axial strain curve in Fig. 10. Comparing the simulation results with the experimental results in Fig. 11(c), the processes observed in the experiments have been captured well, such as the formation of a shear band from the upper right to the lower left portions of the specimen and the fracture branching behavior toward the lower right portion of the specimen. In terms of the fracture mode in Fig. 11(b), shear failure, which is shown by red lines in the figure, is dominant overall.

Fig. 11

Results of FDEM simulation for uniaxial compression test. (a) Fracture process, (b) Fracture mode, (c) Fracture pattern in experiment.

4. Conclusions

The FDEM based on ECZM, i.e., FDEM (ECZM), has been recognized for its advantages as a fracture initiation and propagation analysis method for rocks and rock masses, but no fundamental study on its actual applicability has been sufficiently conducted. In this study, a method of improvement for the problems encountered with the conventional method has been proposed as part of the fundamental study lacking in the conventional method. Specifically, a new FDEM (ECZM) was proposed that does not require remeshing; it introduces a cohesive element insertion method that does not cause any non-physical fracture propagation restrictions. In addition, the analyses of rock fracturing tests, i.e., Brazilian tensile strength and uniaxial compression strength tests, on siliceous mudstone using the proposed method, showed that the mechanical responses, including the fracture patterns and axial stress–axial strain curves observed in the experiments, were captured well, confirming the reproducibility of the proposed method for rock fractures under quasi-static loading conditions.

Acknowledgments

This work was supported by JSPS KAKENHI (Grant Number 20K14826), and the research grants funded by the Takahashi Industrial and Economic Research Foundation, MAEDA ENGINEERING FOUNDATION, FUKADA GEOLOGICAL INSTITUTE and ENEOS TONENGENERAL RESEARCH/DEVELOPMENT ENCOURAGEMENT & SCHOLARSHIP FOUNDATION, Japan. Their support is gratefully acknowledged.

REFERENCES

Corresponding author

Register with J-STAGE for free!