Heterogeneous Rock Simulation Using DIP-Micromechanics-Statistical Methods

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.*erefore, 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. *e 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. *e acoustic emission, strain energy release, and the effect of rock heterogeneities on the postpeak behavior process were investigated. *e numerical results are in good agreement with experimental data.


Introduction
Rock consists of crystals, grains, and cementitious material.Rock materials are usually made up of several different minerals.ese different individual minerals and components are usually distributed in the geomaterials.ey 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.e 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.ese 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].
e literature review indicates that the digital images have been used for the morphological features in many elds of sciences and engineering including biology, medical sciences, geography, civil engineering, and rock mechanics [8][9][10][11][12][13][14][15].In particular, Yue et al. [1,15] developed a DIP technique to establish the actual microstructures of geomaterials, which is called DIP-based nite 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 di erent individual components.
In the upscaling micromechanical model, a micro-to-macro transition called homogenization leads to evaluate the overall (e ective) elastic properties [20,21].
e advantage compared with macroscopic approaches is that the homogenized approach is able to systematically take into account the mineralogical composition in uences on the mechanical properties of rock material [22].e micromechanical approach tries to relate the physical mechanisms involved in the microstructure evolution and macroscopic behaviors observed in laboratory.
e heterogeneous material is considered as a matrix-inclusion composite.
e e ective properties of materials are obtained by an upscaling method based on the Eshelby inhomogeneous inclusion solution [20,22].
is 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.
e proposed statistical Mohr-Coulomb softening model was implemented into a DEM code.e rock behavior was simulated and the experimental stress-strain curve was reproduced numerically.Comparisons between numerical results and experimental data will be nally presented in order to show the capability of the proposed model to describe the main features of rock mechanical responses.
e acoustic emission, strain energy release, and the e ect of rock heterogeneities on postpeak behavior were investigated.

Rock Heterogeneity Investigation
e material studied is an extrusive porphyritic igneous rock called rhyodacite tu .
e mineralogical compositions, initial porosity, and natural water content of samples were rst investigated.In Figure 1, the petrographic microscopic thin section of a rhyodacite tu sample is shown.is rhyodacite tu 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 (μm − cm), the rhyodacite tu sample is composed of the ne grains of calcite, feldspar, and quartz embedded in the crystalline siliceous matrix.At the macroscopic scale (cm − dm), the rhyodacite tu constituted by the assembly of mineral grains and the crystalline siliceous matrix is considered as a homogeneous continuum.

Digital Image Processing.
e 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.eir 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.
e hue component (H) represents repression related to the dominant wavelength of the color stimulus.erefore, hue is the domain color perceived by human beings.
e saturation component (S) represents how strongly the color is polluted with white.e intensity component (I) stands for brightness or lightness and is irrelevant to colors.In general, hue, saturation, and intensity are obtained by di erent transformation formulate by converting numerical values of R, G, and B in the RGB color space to the HSI color space.e 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 di erent perceived colors in the rock sample are acquired according to the values of H, S, or I of individual pixels, and the di erent material properties (such as Young modulus) are speci ed for each pixel according to its catalog of minerals or colors.In theory, the material properties of di erent 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].
e HSI color 2 Advances in Civil Engineering space was used to investigate the rhyodacite tu rock sample microstructure characteristics.e 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 tu rock sample with the AA section at i 150 is shown.e 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 nd 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.
e 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.
e 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.
resholds 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.

Matrix Properties Determination Based on
Micromechanics.Micromechanics investigates the behavior of the heterogeneities as well as their e ects on the overall properties and performance of a material.An important task of micromechanics is to link mechanical relations on di erent length scales.e entire behavior of the microstructure is interpreted as the mechanical state of a material point on the macroscopic level which thereby is ascribed e ective 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 (A r ) is constant.erefore, the sti ness tensor of the heterogeneous media can be expressed as [20] where φ r and A r are, respectively, the volume fraction and the concentration tensor of the rth inclusion family.e RVE is composed of an isotropic linear elastic matrix with sti ness tensor C s and of a random distribution of spherical inclusions with sti ness tensor C r .To evaluate the homogenized sti ness tensor (C hom ), the inclusion concentration tensor (A r ) 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 A r 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 tu .
In the Mori-Tanaka scheme (1973), therefore, the strain or stress eld in the matrix is, in a su cient distance from a defect, approximated by the constant eld.e 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 elds, however, are neglected in this approximation of defect interaction.It follows that the localization tensor is then given by [22] A r I + P is leads to the following estimate of the e ective sti ness tensor [22]: Because of the isotropy of the constituents and the spherical shape assumption of inclusions, we have [22] where μ r and k r are the shear and bulk moduli, respectively, of the phase r and μ hom and k hom are the homogenized shear and bulk moduli.f r and f s are the volume fractions of the phases.e mineralogical compositions of the rhyodacite tu contain four main phases: calcite, feldspar, quartz, and matrix.It is organized in grains spread in a siliceous matrix.
e rst stage of the homogenization procedure is the de nition of a representative elementary volume (r.e.v.).e observations led us to consider the rhyodacite tu as a fourphase 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.e rhyodacite tu sample can be represented by a four-phase composite with distinct mechanical properties.
is 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.
e elastic properties of the crystalline siliceous matrix (K 0 ; μ 0 ) 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 tu were determined experimentally as K hom 16.6GPa and μ hom 15.7GPa, and the elastic properties of calcite, feldspar, and quartz grains are determined from the existing data found in literature [22,25].For the calcite, k 1 70.8GPa and μ 1 32.7GPa, for the feldspar k 2 54.2GPa and μ 2 25GPa, and for the quartz k 3 36.8GPaand μ 3 44.4GPa.e elastic 4 Advances in Civil Engineering properties of the matrix are calculated by solving the nonlinear ( 4) as E 0 25GPa and ϑ 0 0.3.

Statistical Distribution.
Rock is a heterogeneous material.is 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 speci ed stochastically according to the given Weibull distribution de ned 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 u 0 is a scale parameter expressing the average of the considered mechanical properties.In the Weibull distribution function, the shape parameter m plays a signi cant 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 tted to the calculated probability density-Young modulus diagram as shown in Figure 7.
According to Figure 7, the homogeneity index (m) of the rhyodacite tu sample was calculated to be equal to 6. Also, the mean value of minerals' Young modulus was obtained 62 GPa.

e Proposed Model and Input Data.
e plastic rock behavior is represented by the Mohr-Coulomb model with strain softening.e most well-known failure criterion for rock is the Mohr-Coulomb criterion.e 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.e criterion states that the failure occurs if the magnitude of the shear stress on a speci c plane reaches a critical threshold.e critical threshold is associated with both the cohesion of the rock grains at the plane of failure and friction resistance between them.e friction resistance of the failure surface is dependent on the normal stress imposed on the plane.e strain-softening behavior of rocks is governed by shrinking of the failure criterion with the advance of plastic deformation.
e decline of rock strength with plastic strain is denoted as strain-softening behavior.e 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 (c p ) in comparison with the critical softening parameter (c p * ) in the softening region and it reaches a minimum constant value in the residual region.e critical softening parameter controls the transition between the softening and residual stages.
In ( 6), μ represents one of the strength parameters φ, and C and the subscripts p and r denote the peak and residual values, respectively [27].e rock is assumed to be a heterogeneous material, and its mechanical properties are considered to conform to the Weibull distribution function.
e mean values of bulk and shear modulus are speci ed according to real values obtained from laboratory tests.For simpli cation, it is assumed that the bulk and shear modulus have the same homogeneity index.
e 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 di erent elastic properties, one can produce a heterogeneous material numerically.e proposed statistical Mohr-Coulomb strain-softening model used in the presented analysis was linked to a commercial DEM code as a separate constitutive model.e studied rock is an extrusive porphyritic igneous rock called rhyodacite tu .e mineralogical compositions and initial porosity of samples were rst investigated.is rhyodacite tu 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].e complete stress-strain curve of the rhyodacite tu under the UCS test condition is shown in Figure 8.
Hence, numerical simulation of the rhyodacite tu 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.

Geometry and Boundary Condition.
Uniaxial compressive strength (UCS) test is the most widely conducted standard test on rock samples.e main objective of this test is to determine the peak strength (σ c ), modulus of elasticity (E), 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 strainsoftening model, it was attempted to simulate the test condition as closely as possible.e sample shape, dimension, input material properties, and loading condition were selected similar to the test condition.e 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 twodimensional plane slice was selected at the center of the sample and analyzed.
e complete fracture characteristic of a numerical specimen under uniaxial loading may be investigated only in a stable displacement-controlled test.e 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.e sample uniaxial loading was simulated imposing a velocity eld 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.
ere 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.
e numerical specimen was discretized with 4096 elements.e 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 speci ed Young modulus of di erent elements is considerably di erent with this value 6 Advances in Civil Engineering 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 di erent.erefore, the spatial distribution of mechanical properties is shown in Figure 10.
However, the e ect of the randomness due to the spatial distribution of mechanical properties of elements was studied thoroughly by Zhu and Tang [29].

e 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.e local stress-strain curves of some elements with di erent sti ness are shown in Figure 11.
In Figure 11, the local stress-strain behavior of elements composing the numerical sample with dissimilar sti ness is signi cantly di erent.e elements with higher sti ness bear more stress, so their stress states approach to the yield surface earlier.However, the elements with lower sti ness bear less stress, so their stress states approach to the yield surface later.
e overall stress-strain curve of the entire numerical sample was shown and compared to the experimental result as shown in Figure 12. e comparison of numerical and experimental results shows that they are in agreement.
To study the in uence 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 speci ed in Table 2.
As the homogeneity index increases, mechanical material properties become more homogeneous and approach that of the homogeneous body.e total envelope of the stress-strain curves for three numerical specimens with di erent homogeneity indices and experimental result can be seen in Figure 13.e 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.us, the higher m is, the more brittle the stress-strain curve.e 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.is suggests that hard (brittle) rock has a higher homogeneity index than soft (ductile) rock.

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.eir 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.e AE rate increases gradually with extension of the microcracks and increases rapidly as the microcracks link together.e rate maximizes when the nal 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.erefore, 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.erefore, 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.e released strain energy of the entire numerical specimen which is the di erence between the work done at the boundary of the model and the total stored and plastic dissipated strain energies can be written as [31] W where U c is the total stored strain energy, U b is the total change in potential energy, and W p is the plastic dissipated work.ese 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.e simulated stressstrain curve and the released strain energy for the numerical specimen with m 6 under UCS test condition is shown in Figure 14.
e stress-strain curve of this numerical specimen is almost linear in the prepeak region, and then the stressstrain 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.e 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.Advances in Civil Engineering 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.erefore, 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.
e release of strain energy versus axial loading curves for the numerical specimens with different homogeneity indices is illustrated and compared in Figure 16.
e strain energy release of the heterogeneous specimen is less than the strain energy release of the homogeneous specimen.e 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.e 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.e more homogeneous rock fails abruptly because of its more brittleness.

Conclusion
e volume fractions of minerals in the microscopic thin section of the rhyodacite tuff sample were calculated precisely by digital image processing.
e unknown matrix properties were determined based on Mori-Tanaka scheme in the framework of micromechanics.
en, 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.
e simulated AE and released strain energy during the loading process are dependent on the homogeneity index.e more the homogeneous numerical specimen, the more the AE and release of strain energy under UCS test condition.

Figure 1 :
Figure 1: Petrographic microscopic thin section of the rhyodacite tu sample.

Figure 2 :Figure 3 :Figure 4 :
Figure 2:e original image with the j-coordinates at i 150 (section AA).

Figure 7 :Figure 6 :
Figure 7: Fitting the Weibull distribution function to the probability density-Young modulus diagram.

Figure 13 :Figure 14 :Figure 15 :
Figure13: e total envelope of the stress-strain curves for three numerical specimens with di erent homogeneity indices and experimental result.

6 Figure 16 :Figure 17
Figure 16: e release of strain energy versus axial loading for the numerical specimens with different homogeneity indices.

Table 1 :
reshold values and volume fractions of the rock.