Dynamic Modulus of Elasticity as a Guideline for Deformation Evaluation by Using the Finite Element Method for Igneous Rocks in the Morelia Region

1. Introduction

Rocks play an increasingly critical role in construction, particularly in environments subjected to high loading rates and confining pressures. They serve as structural elements or aggregates in ceramic and cementitious matrices. Understanding their physical and mechanical properties is essential. Demand for stone aggregates near urban centers accelerates quarry exploitation near city peripheries. However, the lack of systematic experimentation and analysis of rock components often leads to suboptimal material use, depleting nearby quarry reserves prematurely.

Key indicators of stone quality include density (ρ), uniaxial compressive strength (UCS), and modulus of elasticity (Ed). Utilizing these parameters in deformation and displacement models helps identify optimal materials and quarry locations for efficient rock utilization based on strength characteristics. The mechanical and elastic properties of intact rock significantly influence engineering design across various sectors, from tunnels and dams to slope stability and historical building preservation [1]. Despite the crucial role of rock properties, direct testing remains cumbersome, expensive, and time-consuming due to the requirement for high-quality core samples [2,3].

Rocks are characterized through physical and mechanical tests which allow for the estimation of their behavior under engineering requirements. Consequently, specific research studies estimate various parameters in rocks such as uniaxial compressive strength, static and dynamic modulus of elasticity. Additionally, some studies employ methods like finite element analysis to estimate rock behavior. For example, disconnected vesicles in extensive basalt flows significantly affects properties such as compressive strength, modulus of elasticity, Poisson’s ratio, and sonic velocities [4].

Predictive models estimate rock parameters using data such as UCS, elasticity (E), and deformation modulus [5,6,7]. These models employ artificial neural networks for predicting deformation modulus [8]. Mechanical tests on various volcanic rocks have evaluated the relationships between chemical and mineralogical composition, microstructure, texture, and physico-mechanical properties [9]. Numerical analyses through the finite element method have been conducted to examine different adhesives, including cyanoacrylate, polyester, epoxy-based adhesive, and Portland cement paste [10]. Others have investigated the relationship between dynamic and static parameters using different methods for various rocks, finding that, the difference between dynamic and static elastic moduli can be attributed to microcracks and pores in the rocks [11,12,13]

The coupled discrete element-finite difference method (DEM-FDM) offers a numerical framework for analyzing the load-deformation behavior of soft soils stabilized by columns [14]. DEM estimated the elastic properties of damaged rocks [15], while the T-matrix model has shown potential for determining effective stone moduli using the full-wave sonic (FWS) tool [16]. Nondestructive methods have characterized stress-induced damage progression in three rock types by capturing their nonlinear elastodynamic responses [17] and material’s behavior [18]. Finite element methods (FEM) are used in conjunction with the software (Abaqus^®) to obtain rock parameters [19].

Alternatively, a laboratory compression test on granite was replicated using the grain-based discrete and finite element method (FDEM). The study focused on the solid-fragile fracture process, the damage behavior resulting from initiation and cracking, propagation before maximum failure, and the dependence on intrinsic factors such as rock type, grain size, distribution, mineralogy and rock texture [20]. Recently, the main features and capability of the field test-enriched finite element method to effectively simulate rock failure processes with different effects have been studied [21]. Similarly, a field test enriched finite element method has been proposed to simulate crack propagation and coalescence of brittle rocks under mixed-mode loading [22].

While existing studies provide detailed mechanical and physical property analyses, graphical representation of horizontal and vertical rock deformations remains underexplored. This study proposes a finite element method (FEM) model to visualize rock sample deformations under load, facilitating comparative analysis among quarries supplying construction materials in the city of Morelia, Michoacán, Mexico. Such detailed quarry characterization enables precise material specifications tailored for construction needs.

2. Methodology

The methodology proposed in this research, based FEM, requires critical parameters obtained from laboratory tests. These include ultrasonic pulse velocity (UPV), UCS, and density, which is an intrinsic property of the stone. In general, the proposed methodology allows for the determination of horizontal and vertical displacements of the stones, using essential laboratory test data for the proposed model.

2.1. Sampling

The characterized samples were obtained from 10 quarries of stone material located in the surroundings of Morelia (Figure 1). The names of the material quarries stones, sample codes and stone types are given in Table 1.

An image of the El Colegio quarry showing the stones used for crushing is shown in Figure 2. Samples were obtained and 11 cubes (5 x 5 x 5 cm) were cut from each quarry (Figure 3). The cubes were used for characterization (UPV and UCS); the data obtained were used for modelling.

2.2. Characterization

Characterizations were conducted determine the structural properties of the samples. X-ray fluorescence was performed using ANDREAH equipment. SEM was performed using JEOL equipment (model JSM-6400), which has a Bruker X-ray microanalyzer, energy dispersive spectrum (EDS) and the S150A SPUTTER COATER metallizer. The X-ray diffraction was performed by means of a SIEMENS diffractometer, model D5000 X-ray diffractometers, 30 kV, 20 mA at 1.5406 Angstrom wavelength tube.

To determine the physical-mechanical properties, density (ρ) was obtained using a pycnometer, for which the sample was tagged, saturated in water for 24 hours and dried superficially before measuremen[23]. UCS was assessed using the Forney universal testing machine on cubic samples (5 x 5 x 5 cm) evaluated using Equation 1.

UCS=P/A

(1)

where UCS is the uniaxial compressive strength (kgf/cm2), P is the maximum applied load (kgf) and A is the cross-sectional area (cm2). The Forney testing machine is shown in Figure 4a and the cubic specimen is shown in Figure 4b [24,25].

2.3. Determination of the Ultrasonic Pulse Velocity

In rock mechanics, ultrasonic pulse wave velocity tests are becoming increasingly common due to their non-destructive nature, speed, high precision, and low cost [26].

The dynamic modulus of elasticity was determined by means of ultrasonic pulse velocity using the equipment (model 58E0049/B brand CONTROLS); with a frequency scale from 24 to 150kHz.

Based on ASTM D2845-08, the dynamic modulus of elasticity [27] can be determined by means of equation 2.

$$Ed=\frac{\rho V_s^2(3V_p^2-4V_s^2)}{V_p^2-V_s^2}$$

(2)

Where Ed is dynamic modulus of elasticity (kg/cm2), ρ is density (kg/m3), Vs is secondary wave velocity (m/s) and Vp is primary wave velocity (m/s), based on the work of González De Vallejo and Ferrer [28] for this type of rocks (basalt and andesite), Vs≈2/3Vp and ν≈0.25.

2.4. Modeling

For the development of this model, an almost linear behavior with local isotropy was taken, for which the governing constitutive equation was proposed (Equation 3). Transverse isotropy indicates that at each point in the rock, there is an axis of rotational symmetry. The rock exhibits isotropic properties in the plane perpendicular to this axis, known as the plane of transverse isotropy [29]. Some researchers have studied the influence of igneous rock on seismic images, based on the wave equation that originates from Hooke’s law and Newton’s second law for isotropic, ideally elastic, heterogeneous [30]. Additionally, other researchers have studied rock masses according to the characteristic of elastic wave propagation in the medium and the application of the elastic wave method in rock mass engineering considering a quasi-isotropic crack mass [31]. In this research the study is performed on the rock matrix therefore it is justified to consider it isotropic as in the aforementioned researches. To estimate the horizontal and vertical displacements based on a four-node element, equations (4), (5) and (6) were used based on the references Kattan (2008), as a generalization of the rectangular four-node element [32]. Thus, the governing constitutive equation will be given by (Equation 3).

$$\left\{\sigma \right\}=\left[D\right]\left\{ϵ\right\},$$

(3)

Where σ is the stress matrix, D is the material property matrix and ϵ is the deformation matrix. For plane stresses, D takes the simple form:

$$D=\frac{E}{1-\nu }\left[\begin{array}{ccc}1& \nu & 0\\ \upsilon & 1& 0\\ 0& 0& \frac{1-\upsilon }{2}\end{array}\right],$$

(4)

Where E is the modulus of elasticity and ν is Poisson’s ratio. For the basalt and andesite samples, for ν an average value of 0.25 was taken based on reference [28]. The isoparametric element stiffness matrix for a bilinear quadrilateral element with sides can then be written as

$$\left[K\right]\left\{U\right\}=\left\{F\right\},$$

(5)

Where U is the vector of global nodal displacements, K is the matrix of global stiffnesses and F is the vector of global nodal forces.

To obtain the matrix K, the stiffness matrices of each element (k1, k2, k3,...kn) are assembled using Equation 6.

$$\left[k\right]=t{\int }_{-1}^1{\int }_{-1}^1[B{]}^T\left[D\right]\left[B\right]\left[J\right]d\xi d\eta ,$$

(6)

Where t is the thickness of the material, D is the matrix of the material properties, J is the Jacobian matrix corresponding to the transformation between the physical and computational spaces. Bi (B1, B2, B3, B4) are calculated using the bilinear quadrilateral element with natural coordinates of the following Equation [32].

$$\left[B\right]=\frac{1}{\left|J\right|}\left[B_1{B}_2{B}_3{B}_4\right]$$

(7)

The following equation is used to obtain Bi (B1, B2, B3 and B4):

$$\left[B_i\right]=\left[\begin{array}{cc}a\frac{\partial Ni}{\partial \xi }-b\frac{\partial Ni}{\partial \eta }& 0\\ 0& c\frac{\partial Ni}{\partial \eta }-d\frac{\partial Ni}{\partial \xi }\\ c\frac{\partial Ni}{\partial \eta }-d\frac{\partial Ni}{\partial \xi }& a\frac{\partial Ni}{\partial \xi }-b\frac{\partial Ni}{\partial \eta }\end{array}\right]$$

(8)

Where a, b, c, and d are calculated by the following equations

$$a=\frac{1}{4}\left[y_1\left(\xi -1\right)+y_2\left(-1-\xi \right)+y_3\left(1+\xi \right)+y_4\left(1+\xi \right)\right]$$

$$b=\frac{1}{4}\left[y_1\left(\eta -1\right)+y_2\left(1-\eta \right)+y_3\left(1+\eta \right)+y_4(-1-\eta )\right]$$

$$c=\frac{1}{4}\left[y_1\left(\eta -1\right)+y_2\left(1-\eta \right)+y_3\left(1+\eta \right)+y_4(-1-\eta )\right]$$

(9)

$$d=\frac{1}{4}\left[y_1\left(\xi -1\right)+y_2\left(-1-\xi \right)+y_3\left(1+\xi \right)+y_4(1-\xi )\right]$$

The following equation is used to calculate Ni (N1, N2, N3 and N4):

$$N1=\frac{1}{4}(1-\xi )(1-\eta )$$

$$N2=\frac{1}{4}(1+\xi )(1-\eta )$$

$$N3=\frac{1}{4}(1+\xi )(1+\eta )$$

(10)

$$N4=\frac{1}{4}(1-\xi )(1+\eta )$$

To obtain the Jacobian matrix J the following equation is used.

$$\left|J\right|=\frac{1}{8}\left[\begin{array}{cccc}{x}_1& x_2& x_3& x_4\end{array}\right]\left[\begin{array}{cccc}0& 1-\eta & \eta -\xi & \xi -1\\ \eta -1& 0& \xi +1& -\eta -\xi \\ \eta -\xi & -\xi -1& 0& \eta +1\\ 1-\xi & \xi +\eta & -\eta -1& 0\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]$$

(11)

The coordinates were obtained from the bilinear element in the physical coordinate (Figure 5).

3. Results and Discussion

The results from the X-ray fluorescence characterization performed on all of the samples from the proposed quarries can be seen in Figure 6. The predominant compound in all samples was SiO₂ (between 50 and 60 %) followed by Al2O₃ (20 %). All the samples present a similarity in the percentages of the compounds, no differences of compounds were observed among the materials. The importance of the above geochemical data allows the chemical classification of rocks and the understanding of the geological process using different types of diagrams. Conventionally, geochemical data used for these purposes are divided into major and trace elements. Major elements are those elements that predominate in any rock such as Si, Ti, Al, Fe, Mn, Mg, Ca, Na and K [33,34]. With the development of new techniques that allowed the determination of trace and ultra-trace elements, it became essential to determine the chemical composition of rocks over a wide range of concentrations (wt. % to μg g-1) [33]. In this research work, the Total Alkali-Silica (TAS) diagram (Figure 7) was used which consists in classifying according to their ratio between SiO₂ wt. % content and the sum of NaO₂ wt. % and KO₂ wt. %.

Figure 8 (a-d) shows SEM images of the AGC, El Colegio, Comanja and Huiramba quarries. The samples showed a surface composed of different shapes such as grains or slabs and the morphology observed in the four images is similar, the morphology between crushed and volcanic samples is very similar. In the four images a different change in relief is observed, presenting hollows in most cases, this could be interpreted as the porosity of the material.

The physical and mechanical properties were estimated as proposed in the experimental section, the values of density, uniaxial compressive strength and modulus of elasticity are showed in Table 2. The samples with higher densities present higher uniaxial compressive strength (G and S), coming from volcanic material quarries. Samples G and S present the highest values of UCS (110.10 MPa) and density (2567.00 kg/m3) respectively. Volcanic quarries H, C, K, M, CO, T, CU and J showed the lowest values in UCS, density and primary waves (Vp) compared to samples from crushed stone quarries, due to the fact that the uncrushed samples are weaker and more abrasive.

Figure 9 shows the linear correlation between UCS and Ed for the specimens analyzed in this study; the correlation coefficient was R2 = 0.95, which implies an excellent correlation [37,38]. Some previous research works have validated this correlation between uniaxial compressive strength and dynamic modulus of elasticity [39,40], the above data confirms the use of Equations 4 to 6 for an isotropic material.

With the values of maximum applied load and dynamic modulus of elasticity, equations were discretized and coded in MATLAB using 5041 nodes per mesh to model the deformations and displacements. The application of the load at the top was considered, producing a value of uniaxial compressive strength, horizontal and vertical displacements. The values of the displacements obtained were amplified 15 times to appreciate the deformation in the images, which are shown in Figure 10. The value of the dynamic modulus of elasticity obtained allowed classifying the rock and estimating its mechanical behavior.

Figure 10 a-j shows the 5cm x 5cm meshes that represent the profile of each cubic sample and the deformations are observed once the maximum load is applied, the above is obtained by the proposed FEM methodology, the effect of the maximum load on the sample is observed by superimposing the mesh without deformation (gray) and the deformed mesh of each sample (blue). Samples e, i and j in Figure 10 presents the smallest deformations, while sample c presents the largest deformation when the load is applied. The deformations match with the values obtained from the physical and mechanical properties obtained with the characterizations (Table 2).

A comparison of the values of the mean displacements (x ̅ and y ̅) of each stone quarry is observed in Figure 11. The values in shear stress displacement obtained are equal to less than 0.02 cm similar to those obtained by Vasconcelos and Lourenço 2009 [41] and Pinho et al. 2021 [42].

The results obtained are especially important because these quarries provide all the raw materials for the construction of civil infrastructure. The raw materials (volcanic and crushed) were analyzed, the results showed that the samples of crushed materials present better mechanical and physical properties, such as density and uniaxial compressive strength. In addition, the modulus of elasticity is higher compared to the volcanic samples.

Tests performed on different quarries, such as X-ray fluorescence (Figure 6) indicate that the samples come from similar flows, the elements and chemical composition are similar, the SiO₂ values are between 50 and 60% for all quarries.

5. Conclusions

Different samples of stone quarry used for construction near the capital of Morelia, Michoacán, were analyzed. Physical and mechanical tests were performed and the values of density, modulus of elasticity and uniaxial compressive strength were compared. Results from model showed a clear difference between the crushed and volcanic quarry samples (crushed samples Ed 53035.1 MPa on average compared to volcanic samples Ed 9284.31 MPa on average). The samples from the crushed quarries showed lower average deformations and displacements (0.004 cm for the crushed samples and 0.018 cm for the volcanic samples). The sample with the best properties for construction activity corresponds to the El Colegio quarry.

The mathematical model proposed in this work can be applied to comprehend the deformations of the stones used in other regions and to determine the stone quarries with the best properties for construction use. The application of this methodology could aid decision makers by identifying the unexploited quarries that comply the construction requirements and so stopping the overexploitation of actual working quarries.