Effect of Grain Level Anisotropy on Numerical Rock Fracture Behaviour under Dynamic Loading

The effect of rock crystal anisotropy is tested in numerical simulation of rock fracture at the laboratory sample level of scale. For this end, a numerical model based on the embedded discontinuity finite element method is applied here in simulations of uniaxial tension test on granitic rock. The rock material, consisting of Quartz (trigonal), Feldspar (triclinic) and Biotite (taken hexagonal here), is described as a linear elastic (up to fracture) heterogeneous and anisotropic material, which fails upon reaching the tensile strength of the individual constituent minerals. In the simulations of uniaxial tension test at two different levels of strain rate, the effect of averaged versus anisotropic elastic constants is tested. The results demonstrate that the effect of anisotropy is stronger at low strain rates (0.25 s here), resulting in macrocrack planes at different locations for the isotropic and anisotropic material descriptions, while at higher rates (25 s here) its notable effects disappear in the details of multiple macrocracks.


Introduction
Rock is a cemented, heterogenous aggregate of grains formed by different minerals with truly anisotropic material properties. At the mesoscopic level, these features are represented by the microdefects and mineral texture as well as by the inherent anisotropy of the constituent minerals [1]. While heterogeneity description is often included in the numerical studies of rock materials under dynamic loading, the anisotropy of the rock minerals is mostly ignored, and representative isotropic average values are used for elastic constants instead.
It is, however, interesting to test the effect of crystal anisotropy in numerical simulation of rock fracture at the laboratory sample level of scale. For this end, the fracture model based on the embedded discontinuity finite element method by Saksala [2] is applied here in simulations of uniaxial tension test on granitic rock. The rock material, consisting of Quartz (trigonal), Feldspar (triclinic) and Biotite (taken hexagonal here) minerals, is described as a linear elastic (up to fracture) heterogeneous and anisotropic material, which fails upon reaching the tensile strength of the individual constituent minerals. Upon fracture, a discontinuity (crack) is embedded in an element violating the Rankine criterion.
This model is applied in simulations of uniaxial tension test, at two different levels of strain rate, comparing the effect of averaged versus anisotropic elastic constants on the response of numerical rock. As grain anisotropy involves local directions, different cases concerning the grain orientation are tested.

Rock fracture model
A brief description of the rock fracture model based on the embedded discontinuity finite elements is given here. A more detailed treatment and comparison of different enrichment techniques of the finite element method to describe discontinuities (e.g. cracks) can be found in [3].

Finite element with an embedded displacement discontinuity
Consider a body in 3D space discretized with 4-node tetrahedral finite elements having a crack or displacement discontinuity, as illustrated in Figure 1b. The crack surface is represented by displacement discontinuity Γ , which is defined by the normal nd and tangent vectors m1, m2. As the 4-node element results in constant strain field, the displacement jump over the discontinuity plane is also assumed elementwise constant. Assuming infinitesimal deformation kinematics, justified by the brittle nature of rock fracture, the displacement and the strain fields can be written as where denotes the displacement jump, and and are the standard interpolation functions for the linear tetrahedron and nodal displacements (i = 1,..,4 with summation on repeated indices), respectively. Moreover, and denote, respectively, the Heaviside function and its gradient, the Dirac delta function. As the displacement jump is assumed elementwise constant, ∇ ≡ and thus (2) follows from taking gradient from (1). It should be noted that the term containing the Dirac's delta function, Γ ( ⊗ ) , in (2), is nonzero only when ∈ Γ . This term is zero outside the discontinuity, i.e. for all ∈ Ω ± , and can thus be neglected at the global level when solving the discretized equations of motion. In addition, the peculiar decomposition using functions Γ in (1) is used because it simplifies the treatment of the essential boundary conditions by limiting the influence of to a single element only. This is due to the fact that Γ ≡ 0 outside that element. The special function appearing in Γ is chosen as follows: In the present implementation of the embedded discontinuity kinematics, the discontinuity plane has no exact position inside the element where it is introduced. Moreover, continuity of the crack path over element boundaries is not imposed. The displacement jump is treated analogously to plastic strain. Actually, the problem of solving the crack opening increment and updating the stress can for formulated in an analogous manner to plasticity theory [4]. The finite element discretized form of the problem thus reads [4,5]: where ̈ is the acceleration vector, e is the elasticity tensor, Ni is the interpolation function of node i, is the stress tensor, is the traction vector, and is the loading function. In addition, ̂ is the traction defined on the part of boundary, Γ σ , i.e. a traction boundary condition.

Plasticity inspired traction-separation model
The problem of solving for the crack opening increment and the stress update is formulated as the corresponding problem in the plasticity theory. The relevant model components, i.e. the loading function, softening rules and evolution laws, are defined as: where ,̇ are the internal variable and its rate related to the softening law for the discontinuity, and t is the tensile strength while sd is the viscosity modulus. The exponential softening law in (8) is calibrated by the mode I fracture energy GIc, while parameter ℎ is the softening modulus and parameter controls the initial slope of the softening curve. Moreover, ̇ is crack opening increment. The evolution laws (10) have their equivalent counterparts in the plasticity theory. It should also be noticed that the loading function (7) has a shear term multiplied with the shear parameter . The equations (11) are the Kuhn-Tucker conditions imposing consistency. In this model, the strain rate dependence of the material is accommodated by the viscosity component.
A discontinuity (crack) is introduced in an element when the first principal stress exceeds the tensile strength of the material. The crack normal is parallel to the first principal direction, and once introduced, the crack orientation is kept fixed.

Rock anisotropy and heterogeneity description
A general granitic rock is considered in this study. The numerical rock, consisting of Quartz (33%), Feldspar (59%) and Biotite (8%), is described as heterogeneous, anisotropic linear (up to fracture) elastic material. The crystal systems for these minerals are trigonal (-Quartz), triclinic (Plagioclase Feldspar) and monoclinic (Biotite) [6][7][8][9]. However, Biotite is considered here as pseudo-hexagonal and the hexagonal values measured by Alexandrov and Ryzhov [6] are used. The corresponding elasticity matrices are [9]: The elastic constants are given in the mineral crystal coordinates. In order to transfer them to the rock sample coordinate system, i.e. the global coordinates (XYZ) of the model, the following formula is used where is the 66 coordinate transformation matrix (for stress vector in Voigt notation) consisting of the products of the entries of A, which is the 33 direction cosine matrix of the angles between the original and new coordinate axes.
The rock heterogeneity is described by random clusters of finite elements assigned with the material properties of the constituent minerals. More precisely, the element numbers are first assigned with integers from 1 to 3 corresponding to the percentage of the constituent minerals in the numerical rock. Then, the sets representing different minerals are randomly mapped into the mesh resulting in a spatially mesoscopic description of heterogeneity.

Numerical examples
Numerical simulations are presented here. First, the model and material parameters are given. The simulations are carried out with a self-written Matlab code.

Model and material parameters
The elasticity constants of the rock constituent minerals are given in Table 1. The rest of the model and material parameters are given in Table 2  With the homogenized values of Young's modulus and Poisson's ratio, isotropic elasticity tensor is used. The viscosity value given in Table 2 is so low it does not cause notable loading rate hardening at low rates. Finally, the shear control parameter in the loading surface (7) is set to 0.5.

Uniaxial tension tests on numerical rock
Uniaxial tension test simulations at two different levels of strain rate are carried out. The finite element mesh consisting of 206617 linear tetrahedrons is shown in Fig. 2 along with the rock mineral mesostructure. Three cases concerning the grain orientation are considered. In the first (Case1), all the grains (the crystal axes) conform to the global coordinate system. The second (Case2) is the randomly rotated grains scheme while in the third case (Case3) the crystal c-axis (polar axis) is aligned at 45  to the global z-axis. The simulation results for the isotropic case at strain rate 0.25 s -1 are shown in Fig. 2. The predicted failure mode in Fig. 2c of the isotropic specimen is the experimental transverse splitting [10] with a single primary failure plane. Moreover, some secondary crack opening localization can be observed. There are cracks all over the specimen (Fig. 2d) but only part of them open so as to form the failure plane. The resulting tensile strength is 8.2 MPa. The simulation results for the anisotropic sample are shown in Fig. 3.   Fig. 3. Simulation results with anisotropic elasticity at strain rate 0.25 s -1 .
In the anisotropic cases, the failure mode is, again, the transverse splitting but at a location different from that in the isotropic case. In addition, the type of anisotropy, i.e. the different orientations of local crystal coordinate systems, has some influence on the details of the crack plane. This effect is naturally reflected in the average stress-strain responses. The case with the crystal z-axis aligned 45  to the loading axis (Case3) deviates the most from the other cases in that it shows the stiffest loading and the most ductile softening part. However, the tensile strengths are virtually the same for all cases. Finally, the simulation is repeated at higher loading rate with the results being shown in Fig. 4.   Fig. 4. Simulation results with isotropic and anisotropic elasticity at strain rate 25 s -1 .
Experimentally, the 100 times higher loading rate, 25 s -1 , results in multiple cracking of the specimen and to a DIC (dynamic increase factor of tensile strength) of about 2 [10]. However, the simulated results in Fig. 4 attest only the multiple cracking mode while the tensile strength is almost the same as that at the lower strain rate. The predicted post peak softening is, however, much more ductile at the higher strain rate reflecting the fact that more energy is dissipated in the specimen due to multiple macrocracking. The strain rate hardening effect could be improved to match the experimental DIC by adjusting the viscosity modulus. However, within the present model this would make the details of the failure mode blurred as the energy dissipated in the numerical rock would increase. Moreover, this issue, being out of the scope, is not considered here any further.

Conclusion
A numerical study was presented on the effect of mineral crystal anisotropy on the dynamic fracture behaviour of granitic rock. The simulations of the uniaxial tension test allow to draw the conclusion that, compared to results with averaged isotropic elasticity, accounting for anisotropy has relatively more pronounced effect on the failure mode and stress-strain response at low strain rate (0.25 s -1 ) than at high strain rate (25 s -1 ). At the higher strain rate, relative differences in the failure modes disappear due to the abundancy of details. Therefore, mineral crystal anisotropy of rock forming minerals can be neglected at higher loading rates.