FEM-DEM simulation of Uniaxial Compressive Strength (UCS) laboratory tests

Abstract

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

Keywords:
Finite and Discrete element method; Hybrid numerical models; FDEM; FEM-DEM; UCS

1. Introduction

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 Mahabadi, et al. (2012). 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 and 2018). 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.

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

Figure 1
Simulation steps and processes used in Irazu (Irazu,2018).

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 motion at each time-step (Tatone and Grasselli, 2015):

(1) $\mathit{M}\frac{{\mathit{\partial }}^2\mathit{x}}{\mathit{\partial }{t}^2}+\mathit{C}\frac{{\mathit{\partial }}^2\mathit{x}}{\mathit{\partial }{t}^2}=\mathit{F}$

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) $\mathit{C}=\mathit{\mu }\mathit{l}$

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 repulsive 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 (Mahabadi, 2012). 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.

Figure 2
(a) Conceptual model for mode I and II of fractures, (b) the theoretical model of fracture and (c) Numerical model FDEM (Tatone and Grasselli, 2015).

2.1 Model calibration

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

3. Results and discussion

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.

The results suggest the existence of two materials. A more resistant meta-andesite (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).

Table 1
Samples uniaxial compressive strenght results.

Samples	Height (mm)	Diameter (mm)	Young's Modulus (GPa)	Poisson Ratio	UCS (MPa)
CP-1	106.75	36.00	67.55	0.246	37.53
CP-2	109.90	35.90	62.48	0.234	44.30
CP-3	112.30	35.90	85.16	0.213	139.77
CP-4	107.72	36.03	94.32	0.206	117.68
CP-5	99.05	35.80	97.20	0.171	135.85
CP-6	110.45	35.80	92.12	0.158	142.77
CP-7	108.15	35.85	55.41	0.241	63.35
Average	107.695	35.905	79.18	0.210	97.32

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.

Table 2
Parameters used in the simulation of meta-andesite uniaxial compressive strength.

Parameters	Rock	Platen
Height	107.38 mm	12 mm
Diameter	35.88 mm	37.5 mm
Density	2865 kg/m3	7000 kg/m3
Young's Modulus	92.2 GPa	200 GPa
Poisson Ratio	0.19	0.30
Cohesion	30.96 MPa	200 MPa
Friction Angle	38.23º	0º
Tensile Strenght	8.59 MPa	200 MPa

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 hypothesis of plane stress was considered since the friction between the specimen and the platens was considered null. Figure 3 shows the initial model.

Figura 3
Irazu uniaxial compressive strenght 2D model.

The integration time-increment was calculated according to the Irazu manual recommendations. All time-increments are smaller than the values recommended by Tatone and Grasselli (2015) to maintain the stability of the model.

3.1 Determination of the size of the elements

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 elements, 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.

Figure 4
Discretization size sensitivity analysis - UCS.

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.

Figure 5
Discretization size sensitivity analysis - BTS.

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

Figure 6
Platten final velocity analysis.

3.3 Viscous damping factor analysis

Assuming a triangular element in a mass-spring-dashpot system, the coefficient of critical viscous damping according to Irazu (2018) is estimated by:

(3) ${\mu }_{\mathit{c}}=2h\sqrt{\rho E}$

Where h is the size, r 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 stress-strain 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 damping factor, the greater the principal stress (s1) 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 s1 standard deviation and, therefore, it was used in this study.

Table 3
Viscous damping factor analysis.

Viscous Damping Factor	E (GPa)	Poisson Ratio	UCS (MPa)	Std. Dev. (MPa)
0.001	92159	0.192231	106.55	1.844544
0.01	92155	0.192541	111.97	1.816416
0.05	92164	0.192093	126.69	1.638905
0.1	92263	0.191469	134.78	1.361152
0.25	92173	0.191472	148.27	1.371358
0.5	92178	0.191137	152.78	0.7947753
0.75	92174	0.191369	150.69	0.4162142
1	92177	0.191302	150.65	0.4124939

3.4 Contact penalty analysis

The analysis of the contact penalties was carried out by varying exponentially each penalty from 9.22x10¹⁰ to 9.22x10¹² Pa as proposed by Tatone and Grasselli (2015). The magnitude of 10¹³ 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 (σ₁): Fracture Penalty 9.22x10¹² (Pa); Normal Penalty 9.22x 10¹² (Pa.mm) and; Tangential Penalty 9.22x 10¹² (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.

Table 4
Contact penalty analysis.

Penalty			E (MPa)	Poisson Ratio	UCS (MPa)	Std. Dev. (MPa)
Fracture (Pa)	Normal (Pa.mm)	Tangencial (Pa/mm)
10	10	10	92161	0.192128	144.2	2.009624
10	10	11	92123	0.194039	145.01	3.855873
10	11	10	92154	0.192432	144.11	2.299027
10	11	11	92133	0.193674	147.73	4.402130
11	10	10	92176	0.191242	142.22	0.580680
11	11	10	92170	0.191518	141.56	1.287926
11	10	11	92180	0.191095	143.9	0.638572
11	11	11	92174	0.191412	148.67	0.612720
11	11	12	92166	0.191773	150.07	0.893334
11	12	11	92169	0.191633	150.18	0.699118
11	12	12	92179	0.191646	152.23	1.026024
12	12	12	92179	0.191119	150.66	0.416748
12	11	11	92179	0.191089	144.54	0.879189
12	11	12	92204	0.190837	112.23	1.675255
12	12	11	92181	0.191029	148.44	0.517088

3.5 Fracture energies analysis

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

Figure 7
Fracture energy analysis.

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.

3.6 Uniaxial compression strength final simulation

The physical parameters calibrated in the previous steps were used for the simulation of the uniaxial compression test of a specimen with the mean dimensions of the actual specimens used in laboratory tests. Figure 8 shows three moments of the simulation and a picture of a characteristic test result.

Figure 8
a) Model elastic behavior; b)Mode II predominant softening mode; c) Fractured UCS model d) Fractured sample in laboratory test - CP-03.

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.

4. Conclusions

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 calculated 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 (s1) 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.

References

• BRADY, B. H. G; BROWN, E. T. Rock mechanics for underground mining New York: Kluwer Academic Publish, 2004.
• CHAU, K. T. Analytic methods in geomechanics New York: CRC Press, 2013.
• GEOMECHANICA. Irazu 2d Tutorial Manual [S.l.: s. n.], 2018. p. 1-172.
• GOODMAN, R. E. Introduction to rock mechanics, New York: John Wiley, 1989.
• JING, L. A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. International Journal of Rock Mechanics and Mining, v. 40, n. 3, p. 283-353, 2003.
• LISJAK, A. et al. Numerical modelling of the anisotropic mechanical behaviour of Opalinus Clay at the laboratory-scale using FEM/DEM. Rock Mechanics and Rock Engineering, v.4, n.1, p. 187-206, 2014.
• MAHABADI, K. O. Investigating the influence of micro-scale heterogeneity and microstructure on the failure and mechanical behavior of geomaterials 2012. 224 f. Thesis (Doctor in Civil Engineering) - Univesity of Toronto, Toronto, 2012.
• MAHABADI, O. et al. Three-dimensional FDEM numerical simulation of failure processes observed in Opalinus Clay laboratory samples. Journal of Rock Mechanics and Geotechnical Engineering, v. 6, n. 6, p. 1-16, 2014.
• MAHABADI, O. K. et al. Y-Geo: new combined finite-discrete element numerical code for geomechanical applications. International Journal of Geomechanics, v. 12, n.6, Mar. 2012.
• MARDALIZAD, A. et al. Testing and numerical simulation of a medium strength rock material under unconfined compression loading. Journal of Rock Mechanics and Geotechnical Engineering, v. 10, n.2, p. 197-211, 2018.
• MUNIJZA, A.; OWEN, D. R. J.; BICANIC, N. A combined finite-discrete element method in transient dynamics of fracturing solids. Engineering Computations, v.12, n.2, p. 145-174, 1995.
• MUNJIZA, A. The combined finite-discrete element method London: Jon Wiley & Sons, 2004.
• STEFANIZZI, S. et al. Numerical modeling of standard rock mechanics laboratory tests using a finite/discrete element approach. In: CANADA-US ROCK MECHANICS SYMPOSIUM, 3.; CANADA ROCK MECHANICS SYMPOSIUM, 20., 2009, Toronto, Canada. Proceedings […] Toronto: RockEng09, 2009.
• TATONE, B. S. A.; GRASSELLI, G. A calibration procedure for two-dimensional laboratory-scale hybrid finite-discrete element simulations. International Journal of Rock Mechanics & Mining Sciences, v. 75, p. 56-72, 2015.

Publication Dates

• Publication in this collection
  30 Sept 2020
• Date of issue
  Oct-Dec 2020

History

• Received
  07 June 2019
• Accepted
  06 May 2020