Particle flow code simulation of intact and fissured granitic rock samples

Uxí Cstro-Filgueir, Lendro R. Alejno,b,*, Diego Ms Ivrs

a Department of Natural Resources and Environmental Engineering, University of Vigo, Vigo, Spain

b Department of Geology and Geological Engineering, Colorado School of Mines, Golden, CO, USA

c Swedish Nuclear Fuel and Waste Management Company (SKB), Stockholm, Sweden

d Division of Soil and Rock Mechanics, Royal Institute of Technology (KTH), Stockholm, Sweden

Keywords:Numerical methods Artificially fissured samples Rock mass behavior Particle flow code Parallel bond Flat-joint Smooth-joint

ABSTRACT This study presents a calibration process of three-dimensional particle flow code (PFC3D) simulation of intact and fissured granite samples.First,laboratory stress-strain response from triaxial testing of intact and fissured granite samples is recalled. Then, PFC3D is introduced, with focus on the bonded particle models (BPM). After that, we present previous studies where intact rock is simulated by means of flatjoint approaches, and how improved accuracy was gained with the help of parametric studies. Then,models of the pre-fissured rock specimens were generated, including modeled fissures in the form of“smooth joint” type contacts. Finally, triaxial testing simulations of 1 + 2 and 2 + 3 jointed rock specimens were performed. Results show that both elastic behavior and the peak strength levels are closely matched, without any additional fine tuning of micro-mechanical parameters. Concerning the postfailure behavior, models reproduce the trends of decreasing dilation with increasing confinement and plasticity. However, the dilation values simulated are larger than those observed in practice. This is attributed to the difficulty in modeling some phenomena of fissured rock behaviors, such as rock piece corner crushing with dust production and interactions between newly formed shear bands or axial splitting cracks with pre-existing joints.

1. Introduction

Reliably estimating the mechanical properties of rocks and rock masses under compression is a challenging issue in rock mechanics.This behavior can be studied by laboratory compression tests and numerical simulations. The authors reviewed a large number of standard triaxial compression tests on intact and fissured Blanco Mera granite samples in the laboratory of the University of Vigo as the database for forward modeling. According to the results of these tests,simulations with a three-dimensional(3D)particle flow code (PFC3D) were carried out. This code was chosen due to two features: (1) its capability to reliably model the stress-strain response of intact rock by means of the bonded particle model(BPM)as defined by Potyondy and Cundall(2004)and updated by Potyondy (2015), and (2) the possibility of modeling joints embedded in the rock following the strategy defined as synthetic rock mass (SRM) approach introduced by Mas Ivars et al. (2011).

Cundall and Strack (1979) showed that the distinct element method (DEM) could be used to model disc-shaped granular assemblies. Potyondy et al. (1996) demonstrated that a computational methodology of this type could be used for simulating deformation and fracture of rock. Based on these findings,Potyondy and Cundall(2004)proposed a particle flow code(PFC)in which the rock is modeled by a dense packing of different-sized ball-shaped particles bonded together at some contact points and their mechanical responses could be simulated by the DEM. The input micro-parameters include packing, particle size and size distribution, particle and bond stiffness and strength parameters.This produces emerging characteristics of the model being able to represent deformation mechanisms observed on rock samples as well as their resulting properties.

Hazzard et al. (2000) conducted two-dimensional (2D) simulation study on various types of rocks.They concluded that the failure responses observed, particularly in granite, were similar to those derived from laboratory tests. However, these 2D models were incapable of properly simulating some features of rock behaviors,such as formation and propagation of only one crack cluster to eventually make the sample fail. Park et al. (2004) compared the simulation results obtained from PFC in two dimensions (PFC2D)and universal distinct element code (UDEC) for a rock mass with discrete joints. In the PFC2D model, a sufficient number of both large- and small-scale discontinuities were set. The trends of the simulation results obtained with both codes were similar, and the decreases in both elastic modulus and peak strength with increasing joint density were observed.They also realized that the joint density influences the post-failure response, ranging from brittle to ductile as the number of discontinuities increases. These models tried to study the influence of joint occurrence on rock mass response, but they did not focus on simulating actual examples.

Holt et al. (2005) calibrated PFC results against controlled laboratory experimental results based on rock-like glass bead assemblies. The dependence of strength and elastic modulus on cement content observed in the laboratory was reliably reproduced by the model.Yoon(2007)proposed an advanced statistical approach able to calibrate contact-bonded particle assemblies based on “experimental design” and “optimization” for unconfined compression models. According to this complex approach, an optimum set of micro-parameters was computed and used in PFC2D simulations to fit elastic moduli and uniaxial peak strength for various rocks.Nevertheless, the post-failure behavior was not addressed.

Cho et al. (2007) performed another 2D calibration exercise based on the use of clumps (arbitrarily shaped rigid particles consisting of overlapping spheres/discs)to model a synthetic soft rock.They achieved a reasonable fitting of the unconfined strength,stiffness and strength envelope. Nevertheless, the non-linearity occurring in the initial load levels of the stress-strain curves in the laboratory was not modeled.Additionally,they did not analyze the post-failure behavior.

One of the first important attempts to model rock behavior by means of 3D particle codes was developed by Sch?pfer et al.(2009). They generated their models in PFC3D by means of random pores and unconnected contacts to simulate pre-existing cracks. They explored the dependence of the elasticity, strength and friction on the porosity,confining pressure and initial fracture density. The general trends observed in hard rocks in practice were indeed retrieved, but no comparison to actual test results was provided.

BPMs should fit rock behavior under different conditions;however, when calibrating the uniaxial strength, the tensile strength is typically overestimated. For instance, Potyondy and Cundall (2004) found a value of 7.2 for the ratio of uniaxial compressive strength (UCS) to uniaxial tensile strength (UTS) for the Lac du Bonnet granite models, while the ratio of the actual samples was 21.5. It was postulated that the overestimate for tensile strength was due to the use of round-shaped particles. It was also observed that the tensile strength increases with decreasing particle size (Cho et al., 2007). Indeed, the use of spherical rigid particles reduces the calculation time and the model is simpler.But this approach does not seem to completely represent the interlocked friction and dilation, which would require more complex geometrical distributions (Cho et al., 2007). For the case of PFC3D,the rotation of the particles is even less limited than that for PFC2D,which adds an additional moment to the bond and finally could contribute to contact failure. This freedom of rotation of the particles tends to decrease the UCS/UTS ratio, making the obtained yield surface unreliable. Potyondy and Cundall (2004) observed that the friction in the granite calibrated models in the 3D version was too low,and the strength envelope was very flat;both aspects were also confirmed by Castro-Filgueira et al. (2016). To solve this problem,Cundall and his co-workers(Potyondy and Cundall,2004;Potyondy, 2012, 2015, 2018a) developed more realistic contact models,in particular the flat-joint model,which will be used in this work to model intact and fissured granite sample behavior.

One of the significant improvements developed for the original BPM was the flat-joint contact model(Potyondy,2012,2018a).This contact model is based on the parallel-bonded model,where,once the contact breaks, the interface is reverted to a zero-strength interface, which does not withstand relative rotation. However, in a flat-joint model, the interface is segmented, where all the segments are jointed. As the junction of each segment breaks, the interface behavior evolves from a completely jointed state to a completely disjointed and frictional state, and since the flat-joint contact is not eliminated, even a completely failed interface is still able to withstand some rotation of the particles. This feature allows the model to represent the UCS/UTS ratios as observed for different ranges of both soft and hard rocks.The capabilities of the flat-joint model to simulate the behavior of intact joint have been demonstrated and highlighted by various authors including Wu and Xu (2016) for Jinping marble, Yang et al. (2019) for sandstone with PFC2D, or Peng et al. (2017) focusing on grain size.

There are also a number of particle code simulations focusing on the behavior of rough joints (e.g. Bahaaddini et al., 2013) and the actual behavior of in situ rock mass and its caving (e.g. Potyondy,2015; Rafiee et al., 2018). But these studies fall somewhat out the scope of this study, thus they are not reviewed here.

The main focus of this study is the calibration of the numerical models for intact and fissured rocks based on the flat-joint model.The study does not only account for elastic and strength parameters,but also analyze rock sample post-failure behavior,a topic that has not yet been dealt with sufficiently. The modeling exercise presented in this paper is founded in two previous studies carried out by the authors (Castro-Filgueira et al., 2016, 2017). It was first observed that it was possible to model intact samples following the BPM defined by Potyondy and Cundall (2004). As a first approach,the parallel-bond contact model was used;however,it was not able to reproduce the UCS/UTS ratio and the increase of peak strength with confining pressure(Castro-Filgueira et al.,2016).To solve this problem, a new study with the flat-joint contact model was also carried out, although the post-failure behavior was roughly reproduced.In order to improve the models and to understand how the micro-properties of the contact model influenced the macroscopic response,sensitivity analysis was carried out and presented in Castro-Filgueira et al. (2017). That study was the base that the present article builds upon.

2. Previous works in our laboratory

In the rock mechanics laboratory at the University of Vigo, a number of studies have been carried out to characterize the preand post-peak behaviors of intact and fissured rock samples of Blanco Mera granite- a hard granitic rock from Galicia(Spain). In the present study, tests on intact samples and two types of artificially fissured samples were considered to compare the modeling results. The fissured samples include one type with 1 sub-vertical joint and 2 sub-horizontal joints (1 + 2) and other type with 2 sub-vertical joints and 3 sub-horizontal joints (2 + 3), as depicted in Fig.1a. The jointed samples were considered to represent rock masses at laboratory scale. These samples were submitted to triaxial compression tests with unloading-loading cycles and up to a large strain level at various confining pressures. Relationships of axial stress versus axial and radial strains and volumetric strain versus axial strain for different types of samples under a confining pressure of 10 MPa are illustrated in Fig.1b and c, respectively.

Fig.1. (a) Sketches and pictures of granite samples tested in the laboratory;(b) Axial stress(σ1) versus axial (ε1) and radial(ε3) strains curves for three particular samples (intact,1 + 2 and 2 + 3) under a confining pressure of σ3 = 10 MPa; and (c) Corresponding volumetric strain (εv) versus axial strain (ε1) curves for the same samples, including the irrecoverable strain locus. Modified after Arzúa (2015).

The volume change derived from shear distortion of a material is known as dilatancy. As such, dilatancy can be straight-forwardly implemented in the theory of plasticity, which was the only parameter in terms of plastic potential. The material obeys a normality rule if dilatancy coincides with friction angle.In this case,many complex mathematical problems on plasticity theory can be solved in an easier way,but this does not hold in practice(Alejano and Alonso, 2005). Vermeer and De Borst (1984) proposed an expression suitable for defining dilatancy for any stress-strain situation as follows:

Accordingly, dilatancy is the ratio between the plastic volumetric strain rate and a plastic parameter that controls the evolution of plasticity such as γp.To check how dilatancy evolves,one can use the graph representing volumetric plastic strain versus plastic axial strain. This curve is also known as the irrecoverable strain locus (Medhurst and Brown,1998), as illustrated in Fig.1c. When the trend of this curve goes towards verticality, this denotes a highly dilatant material. When the curve tends to be flat, this represents a material that hardly suffers volume increase after peak.

Fig.1c shows that in the initial yielding phase,the material tends to dilate more.This means that as far as the plasticity or damaging of the rock develops,the material tends to dilate less and less until a constant plastic volumetric strain is eventually observed. This is logical, since a material cannot infinitely dilate. To facilitate understanding,a black line corresponding to highly dilatant material and a gray line corresponding to a material suffering a low level of dilatation are shown in this figure.

One can see that lower dilatancy is observed for samples with increasing jointing (see Fig.1c) and that the same occurs for tests with increasing confining pressure. A more detailed and rigorous explanation of dilatancy can be found in Alejano and Alonso(2005)and Zhao and Cai (2010).

The results of the above-mentioned laboratory research are compiled in the work of Arzúa (2015) and others (Arzúa and Alejano, 2013; Arzúa et al., 2014; Alejano et al., 2017; Walton et al., 2018). The most complete set of strength and stiffness results of this rock mechanics program is presented in Alejano et al.(2017), where readers can consult the corresponding values. The most relevant observations regarding post-failure behavior on intact and fissured rock samples are reflected in Walton et al.(2018).

The obtained results in this context indicate that the stressstrain response of a rock mass is significantly controlled by its structure.Therefore,a continuous behavior can be obtained from a complete stress-strain curve, as shown in Fig. 2. This figure sketches the evolution of axial stress-axial strain behavior for a particular level of confining pressure for a typical test with unloading-loading cycles in black color, and over this line, the individual response of intact and fissured samples is represented(upper part of Fig.2).The peak strength decrease is associated with the occurrence of joints or fissures in intact rock samples. As the rock or rock mass is more fissured, the peak strength decreases;however,the residual strength seems to remain constant as pointed out by some authors (e.g. Walton et al., 2019), and the material is less stiff and its ability to dilate is more limited(lower part of Fig.2).

Fig. 2. Evolution of the behavior of intact and fissured samples with an increase in number of discontinuities,in parallel with the decrease in the quality of the rock mass in terms of axial stress-axial strain relation in the upper part and volumetric strainaxial strain relation in the lower part. Modified after Walton et al. (2018).

3. Simulation of intact samples

After carrying out the laboratory investigation,the next step was to simulate the behavior of this well-known rock. Previous attempts to model intact rock behavior by means of continuous methods were made showing moderately accurate results (Arzúa et al., 2013). However, fissured samples presented a number of relevant problems, thus the continuous modeling approach was discarded.First,intact samples have been calibrated in order to set them as basis of the fissured ones. After that, artificially fissured samples have been simulated creating both patterns of fractures:1 sub-vertical + 2 sub-horizontal and 2 sub-vertical + 3 subhorizontal. The modeling of fissured samples is presented in Section 4.

3.1. Initial simulations

In the study of Castro-Filgueira et al. (2016), the authors used the calibration approach suggested by Itasca(2014)adapted for the new version of the code PFC v.5, which incorporated the flat-joint model. Intact rock compression tests under various confining pressures were performed to calibrate the complete stress-strain curves obtained. A parallel bond material model was first tried,but it was unable to properly represent the frictional response observed. Then, the flat-joint model was adopted and the calibration approach suggested by Itasca (2014) was followed to obtain stress-strain curves, which were able to reproduce the main behavioral response of granite samples in relation to elasticity and peak strength. The fundamentals of the parallel bond, flat- and smooth-joint models are briefly reviewed in Castro-Filgueira et al.(2016) and detailed by Potyondy (2018b).

The final micro-mechanical inputs obtained in this first step of the numerical study were reported in Castro-Filgueira et al.(2016).Fig. 3 shows the numerical and laboratory complete stress-strain curves of compression tests on granite samples with confining pressures of 2 MPa and 10 MPa. This figure illustrates the elastic parameters (elastic modulus E and Poisson’s ratio ν), strengths(peak strength σpeak1 and residual strength) and post-failure parameters (dilatancy ψ and drop modulus M). The drop modulus M is defined as the mean negative slope of the σ1-ε1curve after peak strength in the first 50% of softening according to Arzúa and Alejano (2013).

Although the results were deemed to be reasonably representative, we decided to perform a sensitivity analysis with the objective of fine tuning these parameters and also focusing on the influence of input micro-parameters on resulting outputs.

3.2. Sensitivity analysis

Around 30 sets of parameters were used in this intact rock model and ten simulations (i.e. ten random particle packing realizations)were run for each sample to obtain the average results.The main variations of the parameters are briefly described in this section, which were reported in Castro-Filgueira et al. (2017).

Micro-parameter stiffness ratio (K, ratio between normal and shear stiffness in linear model and flat-joint model) was observed to have a large inverse relation with Young’s modulus,while microparameter contact cohesion showed a large direct influence on the peak and residual strengths (Castro-Filgueira et al., 2017). The contact tensile strength showed not only a direct influence on the tensile strength of the rock, but also an inverse influence on Poisson’s ratio.However, it does not affect the Young’s modulus or residual and peak strengths.The effective moduli of the particles and contacts showed to be directly proportional to the elastic Young’s modulus of the rock;however,this parameter does not significantly affect other macro-properties.As the contact friction angle and the friction coefficient increase, the peak and residual strengths also increase and Poisson’s ratio slightly increases,whereas the Young’s modulus and the tensile strength are unaffected. Interestingly, the stiffness ratio is the sole parameter clearly impacting all the macroproperties of the rock.As this parameter grows,the Poisson’s ratio also increases; however, all other macro-properties decrease(Castro-Filgueira et al., 2017).

3.3. Fine-tuned intact rock models

Based on the results of the two recent studies (Castro-Filgueira et al.,2016,2017),a new calibration of the flat-joint contact model was carried out in this study to fit the micro-properties listed in Table 1.They did not change significantly compared with previous proposal (Castro-Filgueira et al., 2017).

The sample size resolution (φv) - the ratio of sample width to average particle diameter - was set to 10. The sample width is 54 mm,and the maximum and minimum diameter ratio Dmax/Dminis equal to 1.66, where Dmaxis 6.74 mm and Dminis 4.06 mm. The 20%variability was considered for the tensile strength and cohesion of the flat-joint model (Potyondy, 2018a). The numerical simulations for both uniaxial and triaxial compression tests were carried out.The results of the elastic parameters(E and ν),strengthsand) and drop modulus M, obtained from simulations and laboratory, are presented in Table 2.

Fig. 3. Complete stress-strain curves of compression tests under confining pressures of (a) 2 MPa and (b) 10 MPa obtained from laboratory tests (solid line) and numerical simulations based on the flat-joint approach (dotted line). The graphs also illustrate the most relevant rock mechanics parameters to be analyzed.

Table 1Micro-properties of the flat-joint model (Itasca, 2014) to simulate the Blanco Mera granite.

Note:ρ-density;g0-installation gap;Nα-elements in circumferential direction;Nr-elements in radial direction;φB-bonded.φG-gapped fraction;Eb*-material effective modulus; Pm- material pressure; En* and E* - effective moduli for the linear and flat-joint models, respectively; Knand K - stiffness ratios for the linear and flat-joint models,respectively;μn-friction coefficient;σt_c-tensile strength;c- cohesion; φ - friction angle. For the micro-tensile strength and cohesion in the flat-joint model,average and standard deviation are provided instead of a constant value.

For a visual interpretation of the results, the complete stressstrain curves equivalent to those obtained in the previous study(Castro-Filgueira et al.,2017) are plotted in Fig.4, including all ten individual simulations, under confining pressures of 2 MPa and 10 MPa, obtained from laboratory tests and the new numerical simulations.In these graphs,although the pre-peak behavior fits in a rather accurate way, the post-peak behavior still shows little discrepancies.

From these results,it can be concluded that the flat-joint contact model is able to represent in a reasonable way the stress-strain behavior obtained in the laboratory for a wide range of confining pressures. However, the post-peak behavior at high confining pressure is not represented realistically. This discrepancy can be caused by the packing parameters;thus it could be interesting to make a sensitivity analysis of these parameters in relation to the material responses. A more detailed interpretation of these results is presented below together with that of the fissured samples.

4. Modeling of fissured samples

4.1. Fissured sample preparation

To create the fissured sample models, the SRM approach (Mas Ivars et al., 2011) was followed. This technique relies on the BPM(Potyondy and Cundall, 2004) to simulate the intact material,combined with the smooth-joint model, which represents the discontinuities. Therefore, the intact rock was used as the base material, and the joints were created using the smooth-joint model. The joints were incorporated in the models according to the design of the fractured samples created in the laboratory(see Fig. 5) and according to the geometric design features (see Table 3).

To simulate the behavior of the joints,the smooth-joint contact model was chosen. After consulting different bibliographic resources (Mas Ivars et al., 2011; Vallejos et al., 2013, 2014; Itasca,2014), the micro-properties chosen to simulate the artificially fissured samples are those listed in Table 4, with the calibrated values.

In the triaxial compression models,the sample is loaded axially in a single cylindrical wall,in which the upper and lower flat walls serve as loading plates. The velocity of the radial wall is servocontrolled to keep a predefined constant confining pressure.85°and 23.6°, respectively, in the 1 +2 fissured samples, and 78.3°and 22.7°, respectively, in the 2 + 3 fissured samples.

Table 2Average results of the elastic parameters (Young’s modulus and Poisson’s ratio), strengths (peak and residual) and drop modulus from laboratory tests and numerical simulations.

Fig. 4. Complete stress-strain curves for confining pressures of (a) 2 MPa and (b) 10 MPa obtained from laboratory tests (blue) and numerical simulations using the flat-joint contact model (red) after fine tuning of parameters reported in Section 3.3.

Fig.5. Fractured samples created using PFC3D with the corresponding names of the fractures in both structural patterns.The dips of the sub-vertical and sub-horizontal joints are

Table 3Geometric properties of the smooth-joint contact model.

Table 4Micro-properties of the smooth-joint contact model (Itasca, 2014).

4.2. Simulation results

With the micro-properties for the intact rock presented in Table 1 and those of the joints shown in Table 4,simulations were run for 1+2 and 2+3 fissured samples under confining pressures of 2 MPa, 4 MPa, 6 MPa,10 MPa, and 12 MPa. For each of these 10 confinement levels,10 simulations were performed(i.e.10 random particle packing realizations).

Tables 5 and 6 present the results of the main geomechanical parameters obtained from these simulations and from the corresponding laboratory tests for the 1 + 2 and 2 + 3 fissured samples, respectively. Both tables include elastic parameters (E and ν) and strengthso have a general overview of the obtained results, Fig. 6 illustrates the complete stressstrain curves obtained for the intact, 1 + 2 and 2 + 3 fissured samples under three confinement levels (2 MPa, 6 MPa and 12 MPa). The general trends of increasing strength and stiffness with confining pressure and decreasing strength with fracturing can be readily observed.

Additionally, the decreasing dilation with confining pressure and fracturing is also deducible, based on the variations of the volumetric strain-axial strain curves in the right lower quarter of each graph. For each row in Fig. 6, the negative slope of the volumetric strain-axial strain curve diminishes with confining pressure, which indicates less dilatancy for increasing confining pressure. Also, the negative slope of the volumetric strain-axial strain curve, representative of dilatancy, diminishes for increasing jointing, indicating less dilatancy with increasing fracturing.

Table 5Average results of the elastic parameters(Young’s modulus and Poisson’s ratio)and strengths(peak and residual)derived from laboratory tests and numerical simulations for 1 +2 fissured samples.

Table 6Average results of the elastic parameters(Young’s modulus and Poisson’s ratio)and strengths(peak and residual)derived from laboratory tests and numerical simulations for 2 + 3 fissured samples.

Fig.6. Complete stress-strain curves for intact,1+2 and 2+3 fissured samples under confining pressures of 2 MPa,6 MPa and 12 MPa obtained from laboratory tests(blue)and numerical simulations using the flat-joint contact model (red).

4.3. Interpretation of results

In order to interpret the results in more detail, the outputs derived from numerical simulations and laboratory tests for intact and artificially fissured samples are interpreted altogether.

4.3.1. Elastic parameters

Young’s modulus was obtained as the slope of the axial strainaxial stress curve between 30% and 60% of the peak stress. The calculated values of Young’s modulus were plotted in relation to confining pressures (Fig. 7). As it can be observed, the results of intact samples from both laboratory tests and numerical simulations show a similar dependence with the confining pressure. On the other hand, the obtained results for 1 + 2 fissured samples show higher discrepancies, although the level of adjustment is quite reasonable with differences less than 4 GPa. Finally, the results for 2 + 3 fissured samples are quite similar in value and dependence with confining pressure.

Poisson’s ratio was calculated as the ratio of the radial to axial deformation, between 20% and 40% of the peak stress. It showed that the Poisson’s ratios of the intact samples tested in the laboratory kept almost constant. However, Poisson’s ratios of the artificially fissured samples decreased significantly when the confining pressure increased(Fig.8).This can be physically explained by the fact that, since the Poisson’s ratio is a representation of the radial deformation of the sample,the higher the confining pressure is,the more difficult it will be for the sample to deform.This behavior has been observed in the numerical simulations performed. One can observe that as the confining pressure increases, Poisson’s ratio tends to decrease slightly.

4.3.2. Strength variations

Fig.7. Graphical representation of Young’s modulus for each group of samples obtained from laboratory tests(gray)and numerical simulations(orange):(a)All values obtained,and(b) Average values for each confinement level. Intact rock,1 +2 and 2 + 3 fissured samples are represented by diamond, triangle and cross, respectively.

Fig.8. Graphical representation of Poisson’s ratio for each group of samples obtained from laboratory tests(gray)and numerical simulations(orange):(a)All values obtained,and(b) Average values for each confinement level. Intact rock,1 +2 and 2 + 3 fissured samples are represented by diamond, triangle and cross, respectively.

The peak strength (σpeak1 ) was obtained as the maximum value of the axial stress.The residual strength(σres1 )was estimated as the value of the axial stress for an axial strain of 0.03,typically more or less coinciding with a constant value of strength at the end of the laboratory tests and numerical simulations. It is noted that in the laboratory tests,the axial stress at the final testing stage is constant,whereas in the numerical simulations,further deformation lowers the residual strength level.This should be taken into account when interpreting the residual strength results. When it is impossible to continue the test until 0.03 of axial strain, the residual strength is estimated as the last recorded stress data.

Fig. 9 shows the peak strength results. The results of intact samples obtained from numerical simulations fit in a reasonable way with those obtained in the laboratory,with higher discrepancy for the UCS values.It is conclusive that the observed discrepancies are generally lower than the natural variability of the laboratory results. According to the assessment of obtained strengths of fissured samples, it is convenient to put forward that the differences between the laboratory and simulation results are less than 10 MPa. In general, the fittings tend to be better for 1 + 2 fissured samples than for 2 +3 ones, showing a higher variability.

The results of residual strength obtained from numerical simulations are similar to those obtained from laboratory tests, as shown in Fig.10. In addition, the results confirm that the residual strengths for both intact and fissured samples are almost the same.Therefore,the residual strength does not seem to be influenced by the effect of the rock structure.

4.3.3. Analysis of strength results with generalized Hoek-Brown criterion

Fig. 9. Graphical representation of the peak strength for each group of samples obtained from laboratory tests(gray) and numerical simulations (orange): (a) All values obtained,and (b) Average values for each confinement level. Intact rock,1 +2 and 2 + 3 fissured samples are represented by diamond, triangle and cross, respectively.

Fig.10. Graphical representation of the residual strength for each group of samples obtained from laboratory tests (gray) and numerical simulations (orange): (a) All values obtained, and (b) Average values for each confinement level. Intact rock,1 +2 and 2 + 3 fissured samples are represented by diamond, triangle and cross, respectively.

Although the first attempt is to fit the Hoek-Brown (H-B)failure criterion directly with the purpose of obtaining strength results, it is clear that the results of the fissured samples do not adapt to the H-B criterion for intact rock (Alejano et al., 2017).Therefore,the generalized H-B criterion(Hoek et al.,2002)is more suitable to assess the strength parameters because this criterion includes the geological strength index (GSI) (Hoek et al., 1995),which is an index of the rock mass quality related to the jointing or fracturing level.In this section,the peak and residual strengths will be evaluated. In both laboratory tests and numerical simulations,the intact and artificially fissured samples were considered together.

To apply the generalized H-B criterion,the H-B criterion(Hoek and Brown,1980) for intact rocks was fitted (Eq. (2)), in order to obtain the intact classical H-B rock parameters (the uniaxial compressive strength of the intact rock, σci, and a rock friction parameter, mi). To fit these parameters, the MATLAB software was employed for all the values obtained from numerical simulations and laboratory tests for intact samples. The obtained fitting parameters and their values are shown in Fig.11a.

Fig.11. Representation of average peak strength according to laboratory tests (black)and numerical simulations (red), and Hoek-Brown fitting parameters (Hoek et al.,2002) for intact and 1 +2 and 2 + 3 fissured samples.

with these values for each scenario (laboratory test or numerical simulation),the formula of the generalized H-B criterion was fitted for fissured samples.To do so,in the generalized formula(Eq.(3)),the rock mass parameters (Eqs. (4)-(6)) were shown as GSIdependent (Eq. (7)) by considering the H-B damage rock mass parameter D = 0, since the rock samples are not damaged.

Eq.(7)was applied to determination of the peak strengths of the artificially fissured samples(Fig.11),and of the residual strengths of the intact and artificially fissured samples (Fig.12), both from laboratory tests and numerical simulations. Once the strength envelopes are fitted, the obtained strength parameters and correlation coefficients(R2)are tabulated in Table 7.It should be noted that for the presented fitting parameters of the fissured samples, we have fixed the values of σciand mias derived from the fitting to intact rock results. This provides values of σci= 125 MPa and mi= 43.86 for the laboratory tests and σci=134.8 MPa and mi=40.44 for the numerical simulations.

As observed in Fig.11 and Table 7,relatively good fits with high correlation coefficients are obtained in all cases for peak strength(Fig.11). This allows obtaining equivalent GSI values from laboratory analogous to those of the numerical simulations, particularly for the case of the 1 +2 fissured samples. It is obvious that while the fitting for the 1+2 fissured samples is very good(Fig.11b),the strengths derived from numerical simulations are higher than those corresponding to laboratory tests for the 2 + 3 fissured samples (Fig. 11c). This divergence corresponds to an increase of almost 10 GSI units,as shown in Table 7.In any case,the comparison is considered reasonably good, taking into account the natural variability of the response of these fissured elements. This variability in strength results can be observed in Fig. 9a.

When comparing the residual strengths(Fig.12),the GSI values(around 50) for intact samples are quite similar. However, for fissured samples, there are discrepancies between the results. For laboratory results, the equivalent GSI is less than 50, while for simulation results, the equivalent GSI is larger than 50. This could be induced by the fact that in numerical simulations, the strength decrease rate is slower than that in laboratory tests.This is more for the higher confinement cases, as it can be checked in Fig. 6. Thus the residual strengths could be slightly overestimated by the numerical simulations.

Fig.12. Representation of residual strength according to laboratory tests (black) and numerical simulations (red), and Hoek-Brown fitting parameters (Hoek et al., 2002)for intact and 1 +2 and 2 + 3 fissured samples.

Table 7Results of generalized Hoek-Brown criterion fitting for the peak and residual strengths of intact and artificially fissured samples from laboratory tests and numerical simulations.

In Fig.12, the correlation coefficients for residual strength tend to be lower than those for peak strength. This can be interpreted not only by considering that the residual criterion does not fit the proposed H-B equation well, but also by the influence of some erratic results (outliers) that can bias the mean values, producing less regular results in this way. It is noted that for the residual strength fitting(Fig.12),the correlation coefficients(around 0.9)of the fissured samples are higher than those(around 0.8)of the intact samples.This is attributed to two facts. First, the residual strength may not fit the H-B failure criterion as well as intact strength does.Second,the residual response of fissured samples seems to be more regular than that of the intact ones.

The results ofresidual strengthsshowequivalent GSIvalues of37-57,and apparently the residual strengths are similar for all laboratory tests and numerical simulations.This is possible evidence that once broken after deformation,all the samples behave in a similar way at a certain strength level. Similar results were obtained by Cai et al.(2007) based on rock mass observations at field scale and Gao and Kang (2016) based on distinct element models. These considerations led Alejano et al.(2017)and Walton et al.(2019)to conclude that the residual strength of the intact rock can reach similar values of the residual strength of a rock mass under particular circumstances.

Another curious observation is that in all laboratory and numerical cases,the residual strengths under low confining pressures(2-4 MPa) are below the fitting line, whereas they are above this line for larger confining pressures (10 MPa or more). This effect seems to be clearer for the numerical cases(Fig.12).This may be an indication that the H-B strength criterion may not well fit the residual strengths.Interestingly,this trend coincides with the results of Bahrani and Kaiser(2013).They put forward that if one assumes a constant value for GSI, the generalized H-B failure criterion systematically tends to significantly underestimate the strength of fissured rock masses for high confinement levels at field scale.

The simulated residual strengths are considered reasonably representative of those observed in the laboratory, although some aspects of the simulation could be fine-tuned in order to improve the models. It is noted, however, that the observed residual strengths are the results of a number of complex failure mechanisms interacting among them (shear-band and axial splitting occurrence interacting with pre-existing fractures, rock piece corner crushing and rock piece movement,etc.),thus the variability observed can be rated as acceptable.

4.3.4. Post-failure behavior

As shown in Fig.13,the variations in volumetric deformation,i.e.less dilatation under greater confining pressure, are recovered naturally in the simulation as the damaging process progresses.Note also that at the strains up to 0.005 or 0.01,the comparison of results is reasonably approximated; subsequently, the volumetric deformations tend to increase well over the values observed for the laboratory samples.

Numerical models were able to show a reasonably realistic representation of some of the micro-mechanisms that lead to ultimate failure of samples.In the case of intact rock samples,failure is produced in the laboratory tests,especially for the most confined cases,associated with the formation of a shear band.In less confined cases,failure is originated by axial splitting phenomena.In the numerical models performed, it has been observed that the numerical simulations represent only the breakage of the contacts by shear stress.

In Fig.14, the numerical and actual final states of samples are presented for the case of intact rock under confining pressures of 2 MPa and 10 MPa. In this figure, it can be seen that failure occurs following a shear plane,both for low and high confining pressures.This is due to the fact that the tensile strength of the contacts is,on average, 4 MPa and their cohesion is 132 MPa. This is the reason why traction fractures are produced earlier and in greater quantity(see Fig. 14c). The final macro-fracture of the samples can be described as a clear shear band for the sample with 10 MPa of confining pressure both for the numerical model and laboratory tests(Fig.14b). However, for the 2 MPa confinement case, an axial splitting mechanism could better reflect the laboratory observations on the broken sample, which is also consistent with the simulation response(Fig.14a).

On the other hand, for the case of pre-fissured samples, failure tends to partially follow pre-existing joints.Unfortunately,it is not always possible to clearly identify these mechanisms due to material deterioration in the laboratory. However, it is generally observed that failure tends to mainly follow pre-existing fissure planes, typically in the less confined cases. On the other hand, the failure surface typically combines newly formed shear bands with pre-existing fissures in the more confined cases.Similar trends can be derived from the PFC3D models. In general, for 1 + 2 fissured samples, failure tends to follow the pre-existing joints for the less confined cases and the occurrence of shear bands is observed in the more confined cases.In Fig.15a and b,it is observed how the failure follows more the pre-existing planes for the less confined cases than for the more confined ones.

For the 2 + 3 fissured samples, in Fig. 15d, e and f, one can observe the clear occurrence of a newly formed shear band that passes through intact pieces of the sample. Interestingly, this phenomenon is observed for the physical laboratory sample(Fig.15f) as well as for the numerical model vector representation(Fig. 15e). Anyhow, this shear band combines with pre-existing joints to produce the ultimate failure surface with residual strength in the sample. Fig. 15b and e indicates that for the high confined case,there is a tendency for the occurrence of new shear bands across which the ultimately failure would take place. This could also help to explain the unexpected high values of residual strength observed in the laboratory and simulation results in line with practical observations (Bahrani and Kaiser, 2013).

Due to difficulties in removing the broken samples from the plastic sleeve, sometimes the failure mechanism cannot be clearly identified. It has also been observed that broken samples of the same rock under the same confining pressure present different failure patterns. These make it difficult to rigorously interpret the results.However,general observations of failure mechanisms from laboratory tests and numerical simulations are consistent. Moreover, results of the presented numerical models are more realistic than those of other continuous simulation methodologies.

5. Discussion

PFC is capable of simulating the dependency of the macroproperties of intact and fissured rock samples on confining pressure.From the presented results,it is evident that the behavior of a rock mass is controlled to a large extent by its structure.The fitting of the elastic and peak parameters is quite good. The residual strength is not so accurate in line with other approaches.Therefore,although the development of these simulations is time-consuming in preparation and calculation, these results show that PFC3D is a simulation tool with capacities greater than other existing tools to simulate this rock behavior.

Fig.13. Response of volumetric deformation with confining pressure in three cases tested: (a) Intact rock samples; (b) 1 +2 fissured samples; and (c) 2 + 3 fissured samples.

Fig.14. Representation of the intact samples by displacement vectors,where the shear band(left)and the contacts that break due to shear failure occurring in the sample(center;in red color) can be observed, compared to the fracture of a real sample (right), under confining pressures of (a) 2 MPa and (b) 10 MPa (Castro-Filgueira et al., 2016); and (c)Representation of a sample with traction (green) and shear (maroon) fractures.

The post-failure behavior of fractured rock samples is complex and influenced by various aspects,such as formation of shear bands and tensile cracks, the interactions of these with pre-existing cracks, as well as the crushing of corners of rock blocks and the formation of dust material that in turn interacts with the previous phenomena (Hoek, 1983). That is the reason why it is difficult to obtain good results in the simulation. In any case, the two most marked trends of dilatancy, which indicate their decrease with confining pressure and level of damage or plastification suffered by the samples, are reflected by the models, albeit in a non-accurate way.

The decreasing trend to dilate with the level of jointing of the rock sample is also reflected in different models performed, in the same way as it was observed in the laboratory, and in line with Hoek and Brown (1997) proposal on rock mass behavior with increasing number of joints (decreasing GSI). However,it does not seem that other discrete element based approaches (Lisjak and Grasselli, 2014) or SRM models based on lattice-springs (Bastola and Cai, 2020) have obtained much more realistic results. Looking to the future improvement of modeling in this respect,the packing parameters should be revised. With a closer packing, new microproperties may be found capable of reproducing intact rock results and they can be later used in fissured sample models.

6. Conclusions

The main object of this study was to model the behavior of artificially fissured samples tested in the laboratory. To carry out these simulations,intact samples were first simulated with the flatjoint contact model. After that, the artificially fissured samples were modeled, based on the flat-joint contact model as the base material, and the smooth-joint contact model to simulate the joints. Once the fissured models were calibrated, the results of intact and artificially fissured samples obtained both in laboratory tests and numerical simulations were interpreted and compared.

The results show that the obtained results for the pre-peak behavior of the laboratory tests and simulations are quite similar.The observed numerical parameters represent rock mass trends well. This suggests that PFC3D models are able to represent the dependence of the elastic and strength parameters on the confining pressure.As for intact rocks,adjustment of the peak strength levels is rather good and those of residual strength are reasonably adjusted. The modeling results match well with the averages obtained in the laboratory and are completely within the range of the observed values.

However, there are certain discrepancies in the behavior of the post-failure volumetric deformation. Indeed, the volumetric strain obtained from numerical simulations is typically larger than that observed in laboratory tests. For this, the study has proposed an improvement in the capabilities available for more accurately modeling the laboratory stress-strain response of intact and fissured rock samples. The peak strength and elastic parameters obtained in laboratory tests and numerical simulations match well with each other. It is likely that these parameters tend to be the most relevant ones in the design of many actual rock engineering projects. Moreover, this is the first particle-based numerical study where post-failure evolution branches have been carefully presented for intact and fissured rock samples and for a wide range of confining pressures.These branches have shown to follow the same trends as those observed in practice, although there is room for further improvement.

Fig.15. Representation of the 1+2 fissured samples(a,b and c)and the 2+3 fissured samples(d,e and f)by displacement vectors under confining pressures of 2 MPa(left)and 10 MPa (center) obtained from simulations with the PFC3D and derived in laboratory tests (right).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

Itasca Consulting Group is acknowledged by the authors for inviting the first author in their Itasca Education Partnership Program,who additionally provided her with the license of PFC3D.The University of Vigo is acknowledged for financing part of the first author’s PhD studies. The second author thanks the Spanish Ministry of Economy and Competitiveness for funding of the project‘Deepening on the behaviour of rock masses: Scale effects on the stress-strain response of fissured rock samples with particular emphasis on post-failure’, awarded under Contract Reference No.RTI2018-093563-B-I00, partially financed by means of European Regional Development Funds from the European Union (EU). Piedad García and Melissa Formisano are kindly acknowledged for language revision of a version of this manuscript.