FEM-DEMsimulation of Uniaxial Compressive Strength (UCS) laboratory tests

Finite/discrete element methods (FDEM) are hybrid numerical models that use algorithms to analyze the transition from continuous to discontinuous. This type of formulation allows modeling physical laboratory tests with greater proximity to real-ity. This article proposes to simulate the average behavior of a uniaxial compression test campaign. The tests were modeled and calibrated based on the strength and the fracture pattern using Geomechanica Inc. Irazu two-dimensional software. The simulated results were analyzed by the mean standard deviation of approximately 3000 elements in the middle third of the model, the same region where the clip gages are located in the physical test. The obtained results show that FDEM can replicate the laboratory test with great similarity.


Material and methods
The Geomechanica's Irazu software is a hierarchical modeling program combining the finite element and the discrete element methods. A detailed description of this methodology can be found in Munjiza (2004).
The steps of the simulations carried out by Irazu can be seen in Figure   1, including the pre and post-processing using open source tools such as Gmsh or Paraview.
After defining and discretizing the domain using 3-node triangular elements, and assigning the initial and boundary conditions, a second order finite-difference integration is applied to solve the equation of mo-tion at each time-step (Tatone and Grasselli, 2015): Where M is the mass diagonal matrix; x the nodal displacements vector; F the nodal forces vector and C is the damping diagonal matrix. The solution of the equation requires a viscous damping parameter for energy dissipation. The matrix C can be described as: (2) (1) This study proposes to simulate laboratory tests of a meta-andesite using a hybrid numerical model -Finite Element Method / Discrete Element Method (FEM-DEM or FDEM). Uniaxial and diametral compression tests were carried out following the suggestions of the International Society of Rock Mechanics (ISRM) using a specimen diameter of 36 mm (ISRM establishes 54 mm) due to the on-site drilling equipment.
Numerical methods can be classified into continuous, discontinuous and hybrid. The hybrid algorithms are generally used in rock mechanics, mainly for the analysis of rock mass stress-strain behavior (JING, 2003). The hybrid method FDEM was proposed by Munjiza et al. (1995) and modified for geomechanic applications by .
One of the currently available software for rock mechanical modeling is the Geomechanica's Irazu, used in this research.
The uniaxial compression strength (UCS) test is a physical test used to characterize the mechanical behavior of the intact rock. This test, described in Goodman (1989) and Brady and Brown (2004), is one of the most common procedures to directly define the intact rock strength, it's Young´s modulus and its Poisson ratio (Stefanizzi, et al., 2009).
The simulation of Laboratory tests are discussed by some authors. Stefanizzi, et al. (2009) use the ELFEN software to model uniaxial and diametral compression tests. Three-dimensional analyses of laboratory tests using FDEM code were carried out by Mahabadi, et al. (2014). Mardaliza, et al. (2018 used the LS-Dyna to investigate the mechanical behavior of rocks by comparing the uniaxial compression tests performed in the laboratory and numerical model outputs. The compression test simulation implemented in this study was based on the suggestions proposed by Tatone and Grasselli (2015). Initially, calibrations were carried out to define the mesh number of elements, the platten speed, the Viscous Damping factor, the contact penalties and finally the fracture energies. Subsequently, the meta-andesite UCS tests were simulated. The results were analyzed by the mean standard deviation of approximately 3000 elements in the middle third of the model, the same region where the clip gages are located in the physical test.

Results and discussion
The model calibration was based on the work of Tatone and Grasselli (2015), following the steps: (a) definition of the domain, (b) definition of the elements size and, discretization, (c) determination of the platen velocity, (d) selection of the viscous damping, (e) se-lection of contact penalties, and, (f) definition of fracture energies. Tatone and Grasselli (2015) used an algorithm that is based on the force and displacement of the platens' reaction. This study proposes a calibration methodology based on the Irazu's output file. For each time step, the software prints the σ1 mean and standard deviation of the selected elements. Before the specimen's fracture, the standard deviation mean from each time step was used as the calibration parameter to ensure the most uniform uniaxial loading condition.
Several properties obtained in the physical tests, such as density, Young's modulus, Poisson ratio, tensile strength, cohesion, and internal friction angle, were used in the UCS test modeling. The shape of the rupture also assists in the calibration of the model. Table 1 shows the results obtained in the physical compression tests of the meta-andesite.
Where μ is the viscous damping constant and I is the identity matrix (Tatone and Grasselli, 2015). A major advance in the FEM-DEM method, proposed by Munjiza (2004), is the algorithm for detection of contact between elements. The interaction of numerous discrete elements requires fast and accurate contact detections. The algorithm NO BINARY SEARCH (NBS) computes the contacts by investigating if there is an overlap between elements and, in this case, creating re-pulsive and frictional forces when there is a slip between them. The method assumes that these forces are finite. Mahabadi, et al. (2014) recommended using contact penalties: (i) normal compression penalty, (ii) tangential shear penalty and (iii) fracture penalty for the traction. The suggestion for contact penalties is, approximately, 1 to 100 higher than the modulus of elasticity of the material. The fractures are created according to the Morh-Coloumb failure criterion as a mechanical response to the induced stresses, thus allowing a transition from continuous to discontinuous . Irazu uses the concept of energy release rate (G) proposed by Irwin in 1956(Chau, 2013. The software defines 3 modes of fracture energy: the mode I for fracture opening (traction), the mode II for fracture shear and the mode III or mixed, when both phenomena occur. Figure 2 shows the modes I and II of fracture propagation. The results suggest the existence of two materials. A more resistant metaandesite (UCS > 110 MPa) and another one less resistant (UCS < 65 MPa). The simulations were carried out considering the high resistance meta-andesite due to the higher number of tests (CP-3, CP-4, CP-5, and CP-6).
Therefore, the specimens mean dimensions, deformability and strength parameters were redefined as shown in Table 2. The density was obtained by the ratio between the mass and the volume of each specimen. Thirty measurements were carried out. The cohesion and the friction angle were obtained from triaxial tests. The tensile strength was estimated by thirty diametral tensile tests. Table 2 summarizes the specimen and the platen dimensions and properties used in the construction of the model.  Figura 3 -Irazu uniaxial compressive strenght 2D model.
Since the same velocity is assigned to each nodal point of the platen, it can be considered infinitely rigid, non-destructible, regardless of its mechanical properties (Irazu, 2018).
The specimen was deemed to be homogenous and isotropic. The hy-pothesis of plane stress was considered since the friction between the specimen and the platens was considered null. Figure 3 shows the initial model.   Tatone and Grasselli (2015) discussed that the discretization sensitivity analysis of the UCS model must be carried out in conjunction with the sensitivity analysis of the size of the elements of the Brazilian test model, also simulated by the authors. In the case of the Brazilian tensile strength (BTS) model, the number of elements in each discretization is approximately 4 times lower than in the UCS. Therefore, the impact of discretization in the simulation of the indirect tensile strength is considerably greater as shown in Figure  5. The tensile strength doesn't show substantial variations for elements smaller than 0.75 mm. Hence, the chosen size of the elements was 0.75mm, which is the same as that proposed by Tatone and Grasselli (2015). The greater the number of elements of the model, the greater the freedom for the fracture trajectory, however, the greater the computational demand.
The integration time-increment was calculated according to the Irazu manual recommendations. All time-increments are smaller than the values recom-mended by Tatone and Grasselli (2015) to maintain the stability of the model. As a starting point and according to the suggestion of Tatone and Grasselli (2015), an automatic algorithm was used to define the size of the ele-ments, using the following values: 0.65, 0.75, 0.85, 0.90, 1, 1.5, 2 mm. Figure  4 shows that the simulated UCS does not present a significant variation for the different sizes used. The difference between the highest and lowest UCS (6.82 MPa) represents 4.4% of the average strength.

Platten displacement velocity analysis
The axial loading in the UCS simulation is achieved by the platen displacement at a predetermined constant velocity (m/s), as a boundary condition. The velocity was evenly divided between the upper and lower platens. For values smaller than 0.25 m/s the model shows a convergence to a value of UCS of 150.50 MPa, on average. It should be noted that 0.25 m/s is significantly higher than the displacement rate observed in the laboratory tests. The final velocity of 0.06 m/s was chosen in order to ensure the quasi-static condition, with a reasonable computational time. Figure  6 shows the results to the simulations to verify this convergence (Tatone and Grasselli, 2015).

Contact penalty analysis
Assuming a triangular element in a mass-spring-dashpot system, the co-efficient of critical viscous damping according to Irazu (2018) is estimated by: The analysis of the contact penalties was carried out by varying exponentially each penalty from 9.22x10 10 to 9.22x10 12 Pa as proposed by Tatone and Grasselli (2015). The magnitude of 10 13 was also simulated but did not generate a consistent result. The simulations show that the modulus of elasticity variation was 81 MPa. The contact penalties were chosen based on the lowest standard deviation of the principal stress (σ 1 ): Fracture Penalty 9.22x10 12 (Pa); Normal Penalty 9.22x 10 12 (Pa.mm) and; Tangential Penalty 9.22x 10 12 (Pa/mm) attaining a principal stress average standard deviation of 416748 Pa or 0.4 MPa. Table 4 shows the simulation results; the numbers in Fracture, Normal and Tangencial penalties are the expoents used in the simulations.
Where h is the size, ρ is the density and E is the modulus of elasticity of the elements. Tatone and Grasselli (2015) suggest that the viscous damping factor is related to the linearity of the stressstrain curve, mainly at its beginning. The lower the viscous damping factor, the less the linearity at the beginning of the stress-strain curve. This behavior was not observed in this study, probably due to the difference in the analysis used. Tatone and Grasselli (2015) analysis is based on platen reactions while in this study the analysis was based on the actual middle third of the specimen. However, the smaller the viscous damp-ing factor, the greater the principal stress (σ1) mean standard deviation. Table 3 presents the UCS simulations response to several viscous damping factors. Eq. 3 renders the value 1 as the critical viscous damping, which also presents the lowest σ1 standard deviation and, therefore, it was used in this study.

Uniaxial compression strength final simulation
The calibration of the fracture energies aimed to defined the two fracture modes: mode I and mode II (Tatone and Grasselli, 2015). The initial fracture ener-gies used are the ones suggested by Irazu. These lithology simulations showed that the predominant fracture mode for UCS tests is mode II, associated with the shear mechanism. The variation of the mode I fracture energies, corresponding to the traction, did not affect the simulated resistance as observed in Figure 7.
The physical parameters calibrated in the previous steps were used for the simulation of the uniaxial compression test of a specimen with the mean dimen-sions of the actual specimens used in laboratory tests. Figure 8 shows three A total of 82 scenarios were simulated to obtain the laboratory tests mean compressive strength of 134 MPa, finding a value of 45000 μN/mm for the mode II fracture energy. To refine this analysis and to find the conjugated mode I and mode II fracture energy, the results discussed by authors on the simulation of the Brazilian test were used. The combination of the results of the two calibrations yielded the fracture energies mode I and mode II of 16000 μN/mm and 40000 μN/mm, respectively. The uniaxial compressive strength calculated with these fracture energies was 134.17 MPa, which represents an error smaller than 0.15% when compared with the uniaxial compression strength obtained in the laboratory tests.
The Irazu manual suggests that the mode I should be 1 to 25 lower than mode II (Irazu, 2018). The ratio between the two energy fracture modes found in these calibrations is 2.5.

(a) (b) (c) (d)
moments of the simulation and a picture of a characteristic test result. Figure 8 a) shows the modeled specimen and platens used in the simulation. The softening shown in Figure 8 b) is predominantly related to the fracture energy mode II (Shear) in agreement with the graph in Figure 7. Both instances concur that fracture energy mode I (Traction) does not significantly influence the rupture and consequently, does not influence the modeled uniaxial compression resistance. The shearing fracture pattern of the calibrated model resembles the actual discontinuity observed in the specimens used in the physical test (CP_03), as can be seen in Figure 8 c) and 8 d). As the software used (Irazu) is a 2D modeling tool, some differences are expected when the model is compared against the laboratory tests.

References
The discretization element sizes defined for the model are in accordance with the research of Tatone and Grasselli (2015). The size sensitivity analysis must be carried out in conjunction with the Brazilian test which presented a number of elements, approximately, 4 times smaller than in the uniaxial compression strength simulation. The smoothing factor of 0.75 was used in these simulations.
Platte velocities lower than 0.25 m/s returned a constant uniaxial compression resistance of approximately 150.5 MPa, which is also in agreement with the values proposed in the article of Tatone and Grasselli (2015). The chosen speed for each platten was 0.03 m/s resulting in a final velocity of 0.06 m/s, assuring a quasi-static condition while achieving a reasonable simulation time.
The viscous damping factor was simulated using the multiples of the cal-culated critical value. The nonlinearity at the beginning of the stress-strain curve as discussed by Tatone and Grasselli (2015) was not observed in the analyses carried out. Therefore, the mean standard deviation of the principal stress (σ1) was used as a target. The lowest mean standard deviation (0.4 MPa) was attained for critical viscous damping factor. This result corroborates with the viscous damping used by Tatone and Grasselli (2015).
The analysis of the contact penalties did not show significant variations in the modulus of elasticity. The difference between the minimum and maximum values found was 81 MPa. Therefore, the choice of contact penalties was also based on the principal stress mean standard deviation (0.4 MPa) in this case.
The calibration of the fracture energies required a greater number of simulations. The values found for the calibration of fracture energies are 16000 μN/mm for mode I (Traction) and 4000 μN/mm for mode II (Shear). Mode II is 2.5 higher than Mode I, which corroborates with the suggestions of the Irazu manual.
The predominance of softening mode 2 (shear) at the final stages of the UCS modeling substantiate the results found in the calibration of the fracture energies, where the mode I (Traction) fracture energy has no significant influence on the UCS strength. The fracture pattern identified in the model (shearing) is consistent with the fracture pattern showed by the physical test specimens. Since the software used in the simulations is a two-dimensional modeling tool, the pattern of the rupture can show some differences as compared to the physical tests, due to the simplification to a 2D condition.