Abstract
Rock as a natural material is heterogeneous. Rock material consists of minerals, crystals, cement, grains, and microcracks. Each component of rock has a different mechanical behavior under applied loading condition. Therefore, rock component distribution has an important effect on rock mechanical behavior, especially in the postpeak region. In this paper, the rock sample was studied by digital image processing (DIP), micromechanics, and statistical methods. Using image processing, volume fractions of the rock minerals composing the rock sample were evaluated precisely. The mechanical properties of the rock matrix were determined based on upscaling micromechanics. In order to consider the rock heterogeneities effect on mechanical behavior, the heterogeneity index was calculated in a framework of statistical method. A Weibull distribution function was fitted to the Young modulus distribution of minerals. Finally, statistical and Mohr–Coulomb strain-softening models were used simultaneously as a constitutive model in DEM code. The acoustic emission, strain energy release, and the effect of rock heterogeneities on the postpeak behavior process were investigated. The numerical results are in good agreement with experimental data.
1. Introduction
Rock consists of crystals, grains, and cementitious material. Rock materials are usually made up of several different minerals. These different individual minerals and components are usually distributed in the geomaterials. They usually have different physical and mechanical properties and responses under external loading. One of the most important factors affecting the mechanical behavior during the failure process is the inhomogeneities and internal microstructure of geomaterials. More realistic characterizations of the mechanical responses and failure of geomaterials under loading necessitate the consideration of the inhomogeneities and microstructures of the materials. In most of the mechanistic models, the composite geomaterials are always assumed to be homogeneous or piecewise homogeneous and their microstructure behavior is largely ignored [1].
In recent years, attempts have been made by many researchers to examine the behavior or response of geomaterials under loading by taking into account the effects of their material inhomogeneities and microstructures. The heterogeneity and microstructure of rock materials have been characterized by using statistical methods. In this method, the heterogeneity of rock is described by assigning different material properties to the simulated rock sample. These statistical tools can simulate numerically material inhomogeneities that are statistically equivalent to those of actual rock materials with known statistical parameters. Recently, Tang et al. [2, 3] carried out numerical investigations on micro-macro relationship of rock failure under uniaxial compression by taking into account the statistical material inhomogeneity. Using the method in Tang et al. [2, 3], Li et al. [4] further investigated the failure process of Hong Kong granite. Fang and Harrison [5, 6] developed a new technique to simulate the brittle fracture in heterogeneous rocks under compressive loading using commercial finite difference software.
It is usually difficult to adequately specify the statistical distribution parameters in order to reproduce real microstructures in rock. Some recent studies have shown that digital image processing (DIP) can be used to study and determine the rock heterogeneity [7].
The literature review indicates that the digital images have been used for the morphological features in many fields of sciences and engineering including biology, medical sciences, geography, civil engineering, and rock mechanics [8–15]. In particular, Yue et al. [1, 15] developed a DIP technique to establish the actual microstructures of geomaterials, which is called DIP-based finite element method.
In order to determine the matrix elastic properties of the studied rock sample, a micromechanical modeling of the mechanical behavior in the elastic regime was necessary. Although some researchers such as Zaitsev and Wittmann [16], Wittmann et al. [17], Bazant et al. [18], and Schlangen and van Mier [19] have proposed micromechanical models to take the inhomogeneity into account, these models are generally based on the fracture mechanics background without the sound and rigorous upscaling methods. In upscaling micromechanical theory, the rock material is considered as a composite material comprising different individual components. In the upscaling micromechanical model, a micro-to-macro transition called homogenization leads to evaluate the overall (effective) elastic properties [20, 21].
The advantage compared with macroscopic approaches is that the homogenized approach is able to systematically take into account the mineralogical composition influences on the mechanical properties of rock material [22]. The micromechanical approach tries to relate the physical mechanisms involved in the microstructure evolution and macroscopic behaviors observed in laboratory. The heterogeneous material is considered as a matrix-inclusion composite. The effective properties of materials are obtained by an upscaling method based on the Eshelby inhomogeneous inclusion solution [20, 22].
This paper is intended to present an incorporation of digital image processing, upscaling micromechanics, and statistical methods for the mechanical analysis of geomaterials by taking into account their actual inhomogeneities and microstructures.
The proposed statistical Mohr–Coulomb softening model was implemented into a DEM code. The rock behavior was simulated and the experimental stress-strain curve was reproduced numerically. Comparisons between numerical results and experimental data will be finally presented in order to show the capability of the proposed model to describe the main features of rock mechanical responses. The acoustic emission, strain energy release, and the effect of rock heterogeneities on postpeak behavior were investigated.
2. Rock Heterogeneity Investigation
The material studied is an extrusive porphyritic igneous rock called rhyodacite tuff. The mineralogical compositions, initial porosity, and natural water content of samples were first investigated. In Figure 1, the petrographic microscopic thin section of a rhyodacite tuff sample is shown. This rhyodacite tuff sample was cored at the depth of 113 m in order to site investigation of a civil underground project located in the north-west region of Tehran.
At the mesoscopic scale , the rhyodacite tuff sample is composed of the fine grains of calcite, feldspar, and quartz embedded in the crystalline siliceous matrix. At the macroscopic scale , the rhyodacite tuff constituted by the assembly of mineral grains and the crystalline siliceous matrix is considered as a homogeneous continuum.
2.1. Digital Image Processing
The digital image consists of a rectangular array of image elements or pixels. At each pixel, the image brightness is sensed and assigned with an integer value named as the gray level. Their gray levels have the integer interval from 0 to 255 and from 0 to 1, respectively. As a result, the digital image can be expressed as a discrete function f(i, j) in the i and j Cartesian coordinate system [1].
As an alternative to the RGB color space, the hue, saturation, intensity (HSI) color space may be used, as it is close to how humans perceive colors. The hue component (H) represents repression related to the dominant wavelength of the color stimulus. Therefore, hue is the domain color perceived by human beings. The saturation component (S) represents how strongly the color is polluted with white. The intensity component (I) stands for brightness or lightness and is irrelevant to colors. In general, hue, saturation, and intensity are obtained by different transformation formulate by converting numerical values of R, G, and B in the RGB color space to the HSI color space. The values of S and I vary from zero to one. But the value of H varies from 0 to 360, which can also be normalized to be from 0 to 1. Distinct microstructures (such as fractures and minerals) with different perceived colors in the rock sample are acquired according to the values of H, S, or I of individual pixels, and the different material properties (such as Young modulus) are specified for each pixel according to its catalog of minerals or colors. In theory, the material properties of different minerals or structures must be known based on mineralogical analysis of the rock sample, by this means, the relation between values of I (H or S) of the digital image pixels and their material properties can be uniquely established [7]. The HSI color space was used to investigate the rhyodacite tuff rock sample microstructure characteristics.
The commonly used image enhancement method called histogram equalization transformation and noise removal methods are adopted here. In Figure 2, an RGB image of the rhyodacite tuff rock sample with the AA section at i = 150 is shown. The variation of intensity component (I) with the j-coordinates at i = 150 (across the AA section) is illustrated in Figure 3.
For this original image in Figure 2, the numbers of the scanning lines along the i-axis and the j-axis are 555 and 416, respectively. Using the above image, we extract brightness levels along the line at i = 150 in Figure 2 and draw the variation of the brightness levels with the j-coordinate in Figure 3.
Figure 3 shows that there exist two major interface points along the i = 150 line. We find that the brightness level changes abruptly at the corresponding two positions in Figure 3.
Figure 3 shows that a majority of the feldspar minerals have higher brightness levels than those of the calcite minerals and matrices. Furthermore, the brightness levels of the calcite minerals are higher than those of the matrices.
The intensity component (I) variations diagram represents the change in mineral composition and heterogeneity in the rock microscopic section. A histogram of an image is used to display the distribution of brightness values in the image. It is a function to show, for each brightness level, the number of pixels in the image that have that brightness level. Figure 4 shows the histogram of the image brightness levels in Figure 1.
At each brightness level, the number in the vertical axis shows the number of the pixels that has the same brightness value in the image. We can divide the whole image pixels into four groups. Normally, the matrices in the image have low brightness levels and the feldspar minerals in the image have high brightness levels. The threshold value is a brightness level which is a boundary between two kinds of minerals. A trial-and-error method is used to adjust the threshold values so that the best results are obtained. Thresholds of the intensity component (I) values for these four components of the rock and their volume fractions based on the trial-error process are listed in Table 1.
2.2. Matrix Properties Determination Based on Micromechanics
Micromechanics investigates the behavior of the heterogeneities as well as their effects on the overall properties and performance of a material. An important task of micromechanics is to link mechanical relations on different length scales. The entire behavior of the microstructure is interpreted as the mechanical state of a material point on the macroscopic level which thereby is ascribed effective material properties. Such a micro-to-macro transition formally proceeds by appropriate averaging processes and is called homogenization as shown in Figure 5 [21].
Inhomogeneous material can be described by an equivalent homogeneous material. Based on Eshelby solution [23] of equivalent homogeneous material, the concentration tensor of each phase is constant. Therefore, the stiffness tensor of the heterogeneous media can be expressed as [20]where and are, respectively, the volume fraction and the concentration tensor of the rth inclusion family. The RVE is composed of an isotropic linear elastic matrix with stiffness tensor and of a random distribution of spherical inclusions with stiffness tensor . To evaluate the homogenized stiffness tensor , the inclusion concentration tensor must be determined. It is generally used from the Eshelby solution [23] and analytical schemes such as dilute, Mori–Tanaka [24], and self-consistent schemes to evaluate the inclusion concentration tensor in (1) for materials with random microstructure. It is shown that, unlike the dilute and the self-consistent schemes, the Mori–Tanaka scheme describes the in situ experimental data well [22]. Considering the previous results, the Mori–Tanaka scheme is the most suitable for a composite with matrix-inclusion morphology which is the case of the rhyodacite tuff.
In the Mori–Tanaka scheme (1973), therefore, the strain or stress field in the matrix is, in a sufficient distance from a defect, approximated by the constant field. The loading of each defect then depends on the existence of further defects via the average matrix strain or the average matrix stress. Fluctuations of the local fields, however, are neglected in this approximation of defect interaction. It follows that the localization tensor is then given by [22]
This leads to the following estimate of the effective stiffness tensor [22]:
Because of the isotropy of the constituents and the spherical shape assumption of inclusions, we have [22]where and are the shear and bulk moduli, respectively, of the phase r and and are the homogenized shear and bulk moduli. and are the volume fractions of the phases.
The mineralogical compositions of the rhyodacite tuff contain four main phases: calcite, feldspar, quartz, and matrix. It is organized in grains spread in a siliceous matrix. The first stage of the homogenization procedure is the definition of a representative elementary volume (r.e.v.). The observations led us to consider the rhyodacite tuff as a four-phase composite of the inclusion/matrix type in which we discern the calcite, feldspar, and quartz phases, assumed to be distributed individually in a siliceous matrix. The rhyodacite tuff sample can be represented by a four-phase composite with distinct mechanical properties. This material has a matrix/inclusion morphology with the phases randomly distributed, and the calcite, feldspar, and quartz minerals being embedded in the siliceous matrix. It is assumed a representative elementary volume containing four phases as shown in Figure 6.
The elastic properties of the crystalline siliceous matrix (; ) are not precisely known, and there is no direct measurement available. An “inverse method” is therefore used to determine these elastic properties from those which are known for the composite: the macroscopic elastic properties of the rhyodacite tuff were determined experimentally as and , and the elastic properties of calcite, feldspar, and quartz grains are determined from the existing data found in literature [22, 25]. For the calcite, and , for the feldspar and , and for the quartz and . The elastic properties of the matrix are calculated by solving the nonlinear (4) as and .
2.3. Statistical Distribution
Rock is a heterogeneous material. This heterogeneity causes rock in compression to fracture via the formation, extension, and coalescence of microcracks. Studies showed that the variation of mechanical properties can be explained statistically. In a general study on rock fracture, the Weibull distribution function was considered for heterogeneity description.
In this study, the rock is assumed to be composed of many elements of identical size, with the mechanical properties such as bulk and shear modulus of elements to conform to the Weibull distribution, so the mechanical parameters of every element are specified stochastically according to the given Weibull distribution defined in the following probability density function [21]:where u is the variable that follows the Weibull distribution and applied here to both bulk and shear modulus. m is the shape parameter describing the scatter of u and is a scale parameter expressing the average of the considered mechanical properties. In the Weibull distribution function, the shape parameter m plays a significant role. Generally, the higher the value for m, the smaller is the amount of heterogeneity in the model. Following the determination of the Young modulus of each mineral composing the rock sample and its volume fraction in the rock sample, the Weibull distribution function was fitted to the calculated probability density-Young modulus diagram as shown in Figure 7.
According to Figure 7, the homogeneity index (m) of the rhyodacite tuff sample was calculated to be equal to 6. Also, the mean value of minerals’ Young modulus was obtained 62 GPa.
3. The Numerical Simulation
3.1. The Proposed Model and Input Data
The plastic rock behavior is represented by the Mohr–Coulomb model with strain softening. The most well-known failure criterion for rock is the Mohr–Coulomb criterion. The criterion is a linear envelope touching the Mohr’s circles representing the magnitude of the maximum and minimum stresses at the moment of rock failure. The criterion states that the failure occurs if the magnitude of the shear stress on a specific plane reaches a critical threshold. The critical threshold is associated with both the cohesion of the rock grains at the plane of failure and friction resistance between them. The friction resistance of the failure surface is dependent on the normal stress imposed on the plane. The strain-softening behavior of rocks is governed by shrinking of the failure criterion with the advance of plastic deformation. The decline of rock strength with plastic strain is denoted as strain-softening behavior. The strain-softening model allows representation of material softening at postpeak behavior based on prescribed variations of the Mohr–Coulomb model properties (cohesion and friction) as functions of the deviatoric plastic strain after the onset of plastic yield [26].
In the plastic zone, it is supposed that strength parameters of rock mass decrease by bilinear function according to a softening parameter in comparison with the critical softening parameter () in the softening region and it reaches a minimum constant value in the residual region. The critical softening parameter controls the transition between the softening and residual stages.
In (6), represents one of the strength parameters , and and the subscripts and denote the peak and residual values, respectively [27]. The rock is assumed to be a heterogeneous material, and its mechanical properties are considered to conform to the Weibull distribution function. The mean values of bulk and shear modulus are specified according to real values obtained from laboratory tests. For simplification, it is assumed that the bulk and shear modulus have the same homogeneity index.
The proposed statistical Mohr–Coulomb strain-softening model was programmed within the C++ environment and was implemented into a commercial DEM code. Using the Weibull probability distribution function in a numerical simulation of a medium composed of many elements with different elastic properties, one can produce a heterogeneous material numerically. The proposed statistical Mohr–Coulomb strain-softening model used in the presented analysis was linked to a commercial DEM code as a separate constitutive model.
The studied rock is an extrusive porphyritic igneous rock called rhyodacite tuff. The mineralogical compositions and initial porosity of samples were first investigated. This rhyodacite tuff sample was cored at the depth of 113 m in order to site investigation of a civil underground project located in the north-west region of Tehran [28]. The complete stress-strain curve of the rhyodacite tuff under the UCS test condition is shown in Figure 8.
Hence, numerical simulation of the rhyodacite tuff uniaxial compression strength (UCS) test was performed with the proposed statistical Mohr–Coulomb strain-softening model. With regard to the experimental test, a summary of input data used in the numerical analysis is given in Table 2.
3.2. Geometry and Boundary Condition
Uniaxial compressive strength (UCS) test is the most widely conducted standard test on rock samples. The main objective of this test is to determine the peak strength , modulus of elasticity , and Poisson’s ratio of the rock material.
Moreover, employing the sophisticated servo-control testing machine, the complete stress-strain behavior of the rock can be determined in this test. Additionally, the shape of the stress-strain curve in the postpeak region is an indicative of rock breakage mechanism and its brittleness.
In order to verify the statistical Mohr–Coulomb strain-softening model, it was attempted to simulate the test condition as closely as possible. The sample shape, dimension, input material properties, and loading condition were selected similar to the test condition. The main objective was to reproduce the tested rock stress-strain curve numerically and delve into the sample failure mechanism in the postelastic range.
A plane stress condition was assumed for the analysis. It is understood that the actual problem has a 3D nature. But with regard to the 2D nature of the selected code, a two-dimensional plane slice was selected at the center of the sample and analyzed.
The complete fracture characteristic of a numerical specimen under uniaxial loading may be investigated only in a stable displacement-controlled test. The load is applied in a sequence of steps in the vertical direction through incremental axial displacement control at one end of the numerical sample in a quasi-static fashion, while the other end is prevented from vertical movement. The sample uniaxial loading was simulated imposing a velocity field in the range of static loading in accordance with the ISRM standard at the top of the model, while a zero vertical displacement was applied at the base. There are no constraints on the sides of the sample and the specimen sides are allowed to move in the horizontal and vertical directions. A view of the model geometry and employed boundary condition for the test condition is shown in Figure 9.
The numerical specimen was discretized with 4096 elements. The numerical specimen failure process takes place within it due to the heterogeneity of its properties.
As mentioned, the mean Young modulus for the entire numerical specimen is 62 GPa, but specified Young modulus of different elements is considerably different with this value and conform to the Weibull distribution function. Even with the same distribution parameters for the specimen, the spatial distribution of properties of elements may be stochastically different. Therefore, the spatial distribution of mechanical properties is shown in Figure 10.
However, the effect of the randomness due to the spatial distribution of mechanical properties of elements was studied thoroughly by Zhu and Tang [29].
3.3. The Simulation Results
In order to assess the local behavior of the sample, series of horizontal and vertical measuring points were placed throughout the model. Important variables such as stress, strain, and displacement components were monitored at these locations. The local stress-strain curves of some elements with different stiffness are shown in Figure 11.
In Figure 11, the local stress-strain behavior of elements composing the numerical sample with dissimilar stiffness is significantly different. The elements with higher stiffness bear more stress, so their stress states approach to the yield surface earlier. However, the elements with lower stiffness bear less stress, so their stress states approach to the yield surface later. The overall stress-strain curve of the entire numerical sample was shown and compared to the experimental result as shown in Figure 12. The comparison of numerical and experimental results shows that they are in agreement.
To study the influence of homogeneity indices (m) on macroscopic mechanical behavior, especially the postpeak region, three numerical specimens with homogeneity indices of 1, 3, and 6, respectively, are subjected to uniaxial compression loading while other input data remain as specified in Table 2.
As the homogeneity index increases, mechanical material properties become more homogeneous and approach that of the homogeneous body. The total envelope of the stress-strain curves for three numerical specimens with different homogeneity indices and experimental result can be seen in Figure 13.
The stress-strain curves of these numerical specimens are linear in the prepeak region and lose most of their load-carrying capacity in the postpeak region. As m increases, the strength and brittleness of the resulting stress-strain curve increases. Thus, the higher m is, the more brittle the stress-strain curve. The lower m is, the more ductile the model behavior. From the above results, we can conclude that the homogeneity index in this model controls the strength and ductility. This suggests that hard (brittle) rock has a higher homogeneity index than soft (ductile) rock.
3.4. AE and Released Strain Energy
Because of rock heterogeneity, some elements of the numerical specimen under loading reach to the failure criterion earlier than others. Their released strain energy is the origin of the acoustic emission in the rock-fracturing process. Acoustic emission (AE) can be used to detect the microscopic processes associated with heterogeneous rock fracture. Generally, AE events are not notable until the occurrence of nonlinearity in the stress-strain curve, and the rate at which the AE events appear changes with the development of fracture. The AE rate increases gradually with extension of the microcracks and increases rapidly as the microcracks link together. The rate maximizes when the final fracture planes form [9]. Monitoring acoustic emission (AE) event rates seems to be a good way to identify the initiation and propagation of microcracks in heterogeneous rock. In quasi-brittle materials such as heterogeneous rock, AE is predominantly related to the release of strain energy. Therefore, as an approximation, it is reasonable to assume that the AE counts are proportional to the number of failed elements and that the released strain energy by failed elements is all in the form of acoustic emissions [30]. By this means, in this model, each AE event corresponds to failure of an element. Therefore, the AE counts are accounted by the number of failed elements and the released strain energy is calculated based on energy computation of the entire numerical specimen. The released strain energy of the entire numerical specimen which is the difference between the work done at the boundary of the model and the total stored and plastic dissipated strain energies can be written as [31]where is the total stored strain energy, is the total change in potential energy, and is the plastic dissipated work. These energies described in (7) are calculated in an incremental fashion at each timestep from the stress, force, displacement, and strain changes. Based on the above assumptions, the cumulative AE counts and cumulative released strain energy can be realistically simulated with the abovementioned numerical model. The simulated stress-strain curve and the released strain energy for the numerical specimen with m = 6 under UCS test condition is shown in Figure 14.
The stress-strain curve of this numerical specimen is almost linear in the prepeak region, and then the stress-strain curve begins to deviate from linearity at stress about 50 MPa. Further increases in axial strain lead to a rapid released strain energy, and this process continues up to the peak strength. Finally, the strain energy release at the constant rate in the postpeak region and the stress-strain curve approaches a residual strength. The stress-strain curve as well as the AE counts during the fracture process of the numerical rock specimen under uniaxial compression test is shown in Figure 15.
It can be clearly seen that there are a few failed elements during the initial loading phase. But these failed elements release much less energy as shown in Figure 14. Therefore, the curve is nearly linear up to approximately 60% of the peak strength. Localization of deformations may appear due to cracking caused by these failed elements. After reaching the peak strength, the load-carrying capacity of rock drops considerably, followed by a long tail until fracture of the specimen.
The release of strain energy versus axial loading curves for the numerical specimens with different homogeneity indices is illustrated and compared in Figure 16.
The strain energy release of the heterogeneous specimen is less than the strain energy release of the homogeneous specimen. The heterogeneous numerical specimen releases its strain energy gradually and in a controlled manner. However, the homogeneous numerical specimen releases its strain energy abruptly at a stress level about the peak strength. The AE accumulation counts versus axial loading curves for the numerical specimens with different homogeneity indices are shown and compared in Figure 17.
Based on Figure 17, the more the heterogeneous numerical specimen, the sooner and more gradual and ductile the failure process was done. The more homogeneous rock fails abruptly because of its more brittleness.
4. Conclusion
The volume fractions of minerals in the microscopic thin section of the rhyodacite tuff sample were calculated precisely by digital image processing. The unknown matrix properties were determined based on Mori–Tanaka scheme in the framework of micromechanics. Then, the Weibull distribution function was fitted to the distribution of minerals’ Young modulus, and the homogeneity index was determined. Using the statistical Mohr–Coulomb strain-softening model, the rock behavior was simulated and the experimental stress-strain curve was reproduced numerically. From the numerical results, we can conclude that the homogeneity index in this model controls the strength and brittleness. The simulated AE and released strain energy during the loading process are dependent on the homogeneity index. The more the homogeneous numerical specimen, the more the AE and release of strain energy under UCS test condition.
Conflicts of Interest
The authors declare that they have no conflicts of interest.