Effect of inherent microcrack populations on rock tensile fracture behaviour: numerical study based on embedded discontinuity finite elements

Inherent microcrack populations have a significant effect on the fracture behaviour of natural rocks. The present study addresses this topic in numerical simulations of uniaxial tension and three-point bending tests. For this end, a rock fracture model based on multiple intersecting embedded discontinuity finite elements is developed. The inherent (pre-existing) microcrack populations are represented by pre-embedded randomly oriented discontinuity populations. Crack shielding (through spurious locking) is prevented by allowing a new crack to be introduced, upon violation of the Rankine criterion, in an element with an initial crack unfavourably oriented to the loading direction. Rock heterogeneity is accounted for by random clusters of triangular finite elements representing different minerals of granitic numerical rock. Numerical simulations demonstrate the strength lowering effect of initial microcrack populations. This effect is substantially stronger under uniaxial tension, due to the uniform stress state, than in semicircular three-point bending having a non-uniform stress state with a clear local maximum of tensile stress.


Introduction
Rock behaviour under loading is to large extent influenced by heterogeneity. Rock heterogeneity manifests at many scales, ranging from grain level microdefects (micro-and mesoscale) to geologic macroscale of crustal faults. At the grain level, considered in this paper, it has two sources: the mineral mesostructure, i.e. the grain size, shape and material properties of the constituent minerals, and microdefects, i.e. microcracks and voids [1,14,21,24,29,52]. Microcracks originate naturally, during geologic rock forming processes, or are stress-induced via thermal or mechanical loading [16]. With respect to rock microstructure at grain level, microcracks occur at grain boundaries, inside the grains (intragranular cracks), or they can cut across several mineral grains and grain boundaries (transgranular cracks) [16]. In any case, microcracking has a detrimental effect on rock strength. In fact, the actual mechanism of brittle damage is microcracking, which realizes as stiffness and strength degradation.
Due to its importance in geology and rock engineering, the effect of microcracks on rock under various loading conditions has been widely investigated both experimentally [5,14,16,29,49,52] and numerically [5,11,13,16,18,20,22,25,26,53]. The focus in numerical modelling has been on the effect of microcracks under compression as it is the natural state of bedrock. However, as the tensile strength of rock is at least ten times smaller than the compressive strength, tensile fractures are frequently observed in geomechanical applications. Therefore, the microcrack influence on tensile behaviour is of primary interest. This is the topic of present numerical study.
Numerical modelling of rock fracture is a challenging task involving description of displacement discontinuities, i.e. cracks. There are basically two approaches to model rock fracture: the continuum-based methods, such as the finite element method (FEM), and the discontinuous or particle-based methods, such as the discrete/distinct element method (DEM). The particle-based methods are naturally superior to FEM in fracture modelling, due to the underlying discontinuum assumption, and has thus become a popular choice in numerical modelling of failure processes in brittle materials in general and rock fracture in particular [10,12,13,15,16,27,28,47,48,51]. These models are especially attractive, as, in many cases, their rock microstructure description includes explicit grain boundaries and grain boundary cracking. However, the critical drawback of the discrete methods is the computational effort required to keep track and update the particle contact configurations and neighbours. Classical FEM can model rock fracture only in the smeared sense, i.e. as localized deformation, by damage and/or plasticity models [9,34,39,42]. However, continuum models have the advantage of computational efficiency and relative simplicity of calibration and, in most cases, the physical meaning of material and model parameters. On the other hand, their most obvious shortcoming is the poor numerical description of fracture and fragmentation. For this reason, FEM has been enhanced or enriched during the last two decades to better describe discontinuities. This research has resulted in two classes of enrichment methods, which are the embedded discontinuity FEM [18,40,41], or EDFEM, and the extended FEM [2,3], or XFEM. EDFEM has the computational/implementational advantage over XFEM that extra degrees of freedom related to the enrichment of the displacement field can be eliminated by static condensation. EDFEM has been applied in modelling dynamic rock fracture analyses, e.g. by Saksala [35] and Saksala & Ibrahimbegovic [37].
EDFEM is adopted as the numerical method in this study, which investigates the weakening effect of inherent microcrack populations on rock tensile strength. This method is particularly suitable to account for inherent microcrack populations as pre-embedded randomly oriented discontinuities. However, a situation is very likely to occur where a finite element having an initial crack, unfavourably oriented with respect to the present loading axis, cannot open, resulting thus in spurious stress locking or artificial hardening of the element. In order to prevent this crack shielding effect, the concept of multiple intersecting discontinuities in a single element [30,33,37] is required. Saksala and Ibrahimbegovic [37] used this technique in the numerical analysis of thermally induced cracks but the second crack needed to prevent the spurious stress generation was assumed to be independent of the first crack. Here, this unrealistic assumption is mended by proper formulation of the intersecting discontinuities theory so that the cracks in the same finite element are coupled in a manner similar to failure surfaces in multisurface plasticity. The developed method is validated and applied in numerical simulations of heterogenous granitic rock under uniaxial tension and dynamic three-point bending of semicircular disc.
2 Theoretical formulation of the modelling approach 2.1 Finite element with two intersecting discontinuity lines Let a body X 2 R 2 be discretized with 3-node triangular constant strain (CST) elements. Assume now that the discretized body is split into parts by intersecting displacement discontinuities, i.e. cracks. Figure 1a illustrates an element with two displacement discontinuities C d1 and C d2 defined by the normals n 1 , n 2 and tangent vectors m 1 , m 2 .
As this element results in constant strain, the displacement jumps over the discontinuity lines are also assumed elementwise constant, which considerably simplifies the finite element implementation of the embedded discontinuity kinematics.
Assuming infinitesimal deformation kinematics justified by brittle nature of rock fracture, the displacement and the strain fields can be written as where a di denotes the displacement jump for crack i, and N i and u e i are the standard interpolation functions for the CST element and nodal displacements (i = 1,..,3 with summation on repeated indices), respectively. Moreover, H C di and d C di denote, respectively, the Heaviside function and its gradient, the Dirac delta function. The elementwise constant assumption of the displacement jump means that ra di 0, and thus (2) follows in a straightforward manner from taking gradient of (1). Moreover, the terms containing the Dirac's delta function, d C di n i a di ð Þ sym , in (2), are nonzero only when x 2 C di . Outside the discontinuity, this term is zero and can thus be neglected at the global level when solving the discretized equations of motion.
The task of function M C di in (1) is to restrict the effect of a di inside the corresponding finite element, i.e. a di 0 outside that element. This facilitates the finite element implementation of the kinematics as there is no need for special treatment of the essential boundary conditions. The ramp function u C di appearing in M C di is chosen, from among the combinations of the nodal interpolation functions, as illustrated in Fig. 1b, so that its gradient is as parallel as possible to the crack normal n i : The finite element formulation of the embedded discontinuity kinematics is based on the enhanced assumed strain concept (EAS) [3,30,43,44]. Thereby, the variation of the enhanced part of the strain in (2), i.e. the second and the third terms, is constructed in the strain space orthogonal to the stress field. Following, Mosler [30], the L 2 -orthogonality condition with a special Petrov-Galerkin formulation is applied, yielding the following expression for the weak form of the traction balance: Z where b di is an arbitrary variation of the displacement jump i, r is the stress tensor, t C di is the traction vector for crack i, A e is the area of the element, and l di is the length of the discontinuity C di . For CST element, the integrands in (4) are constants, which means that the weak traction continuity reduces to the strong (local) form of traction continuity. The final FE discretized form of the problem is then written as [30,36]: where € u j is the acceleration vector, N nodes is the number of nodes in the mesh, and N i is the interpolation function of node i. In addition,t is the traction defined on boundary C r . Equation (5) is the discretized form of the balance of linear momentum, while Eq. (6), with / di being the loading function defined in the next section, defines the elastic zone of stresses. In should be emphasized that this EAS-based formulation results in a simple implementation without the need to know explicitly neither the exact position of the discontinuity within the element nor its length.
Finally, the choice of the CST element coupled with the constant displacement jump assumption is justified. Namely, it leads to the simplest formulation, and thus implementation, of the embedded discontinuity theory as can be observed above. As for the other possibilities, such as using linear displacement jump with linear elements [23] or higher-order bulk elements with constant displacement jump [6,31], they lead to much more complicated finite element implementation: In the former option, the gradient of the displacement jump is nonzero, which leads to a more complicated expression of strain (Eq. (2)) and, consequently, the weak expression of traction balance (4) does not reduce to the strong form in Eq. (6). This would complicate the traction-separation model to be defined in the next section. In the latter option, the construction of ramp function (Eq. (3)) becomes more involved, and, as in the previous option, the weak expression of traction balance (4) does not reduce to the strong form either due to the nonlinear bulk part of the stress. However, these options may obviously perform better in many problems as well as alleviate the locking problems exhibited by the present combination. Notwithstanding, the locking problems can be alleviated by other means, as done in Sect. 2.4. Moreover, it will be shown in the numerical simulations that the linear bulk element-constant displacement jump combination works surprisingly well.

Plasticity inspired traction-separation model for tensile fracture
A glance at Eq. (5) and (6) reveals their formal similarity to the plasticity theory, which allows to recast the problem of solving the irreversible crack opening increment and the evolution equations in the computational plasticity format [30]. This means that the classical elastic predictor-plastic corrector split can be employed. The relevant model components, i.e. the loading function, softening rules and evolution laws are defined as where E is the elasticity tensor, j di ; _ j di are the internal variable and its rate related to the softening law q di for the discontinuity i, and r ti is the tensile strength while s d is the viscosity modulus. Parameter h di is the softening modulus of the exponential softening law, and parameter g d controls the initial slope of the softening curve and it is calibrated by the mode I fracture energy G Ic . Moreover, _ k di is crack opening increment. It should also be noticed that the loading function (7) has a shear term multiplied with the shear parameter b. Equation (11) shows the Kuhn-Tucker conditions imposing the consistency, meaning thus that this is the viscoplastic consistency formulation by Wang et al. [46].
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. Initial crack populations always present in natural rocks can be easily described in this approach as pre-embedded discontinuity populations with desired orientations and initial strengths.

A criterion to introduce the second crack
The crack shielding upon rotating stress states mentioned in Introduction occurs easily when modelling the inherent microcrack populations with random orientations in rock material. This shortcoming of the embedded discontinuity approach can be easily mended by introducing another crack in a finite element with a crack unfavourably oriented to the loading. A criterion modified from that by Saksala and Ibrahimbegovic [37] is employed here as: where n 1 is the normal of the initial crack, n * is the principal direction corresponding to the present major principal stress r 1 , exceeding the modified tensile strength r Ã t , which is a convex combination of the tensile strength of the initial crack, r ic t , and the intact tensile strength r t0 of the rock mineral. The strength of the initial crack is an adjustable parameter reflecting the integrity of initial cracks. Moreover, a new crack is introduced (only once) when the angle between the old crack normal and the new principal direction is greater than a ¼ acosðC a Þ, with C a being an adjustable parameter. A value 1= ffiffi ffi 2 p corresponding to a ¼ 45 , a compromise between 0 and 90°, is used in this study. It should be mentioned that slight alteration of this parameter did not result in significant differences in simulations.

Isotropic damage model to alleviate locking problems
Low-order finite elements with embedded discontinuities suffer from crack locking and spreading problems under mixed mode I/II loading [33,38]. Fortunately, this problem can be alleviated, e.g., by the linear displacement discontinuity formulation [23], multiple intersecting discontinuities and rotating discontinuity approaches [30,33]. Here, a simple isotropic damage model, used by Jirasek and Zimmermann [19] with the rotating crack models to prevent spurious stress generation, is employed. In this method, a transition from embedded discontinuity model to isotropic damaging is made in an element with a discontinuity after the displacement jump has reached certain threshold opening. The damage variable is defined as where g d is defined in Sect. 2.2, a d is the history variable for the damage model being defined as the maximum value reached (during the loading process) by the norm of the sum of the displacement jumps a di . Moreover, a tr is the threshold value indicating the transition to the damage model. Finally, q tr =q 0 represents the ratio of the final (softened) and initial strength. By this model, spurious stresses are effectively eliminated, as the stress is multiplied with 1 À x (likewise to classical scalar damage models): It should also be noted that (16) is the final stress, not the trial stress.

Local and global solution methods
The flow of the whole solution process, including an explicit time marching scheme to solve the global problem, is illustrated in Fig. 2. The programming related parameter tag controls how many times an extra crack is introduced.
The update formulae required to perform the stress return mapping is based on the standard elastic predictorplastic corrector split. First, it is noted that by Eq. (10), Moreover,  the  trial  stress  violates  both  criteria, i.e.
Then, the update formulae are given in Box 1. In Box 1, Eq. (17) is The calculations related to this model are thus fully analogical to the plasticity models. It should be stressed that neither crack tracking algorithm is needed nor continuity over the element edges is imposed. Despite this, it will be shown in the numerical simulations that the model performs surprisingly well.

Numerical examples
This section presents the numerical simulations that demonstrate and validate the model. First, the model behaviour, especially the second crack scheme, is demonstrated with a two-element mesh. Then, the model prediction is validated against an experimental result in a mixed mode I/mode II crack propagation. Finally, the model applied to investigate the effect of initial crack populations in uniaxial tension and three-point bending of half disc. All the simulations are carried out with a self-written MATLAB code.

Material properties and model parameters
The material properties and model parameters used in the initial crack population simulations are given in Table 1. The numerical rock consists of Quartz (33%), Feldspars (59%) and Biotite (8%) minerals, and most of their properties are taken from Zhao et al. [53], Mahabadi [10], and Vasquez et al. [45]. The rock heterogeneity is described by random clusters of finite elements assigned with the material properties of the constituent minerals.
The viscosity values in Table 1 are small enough so as not to cause notable strain rate hardening effect in low rate simulations, and high enough to secure robust solution of Eq. (17). More precisely, the diagonal entries of matrix G should be positive to have a solution that fulfils the consistency conditions. However, if the element size is relatively large and, consequently, softening modulus h di very large, the diagonal entry may be negative in the inviscid case (s d = 0). Nevertheless, in the viscid case with explicit time marching, the term s d /Dt is large enough to guarantee positivity of the diagonal entry and hence the robust solution of (17).
Finally, it is noted that the number of material and model parameters in this model is exceptionally small.

Two-element model simulation:
demonstration of the second crack scheme The first simulation demonstrates the model behaviour with a two-element model shown Fig. 3a. The purpose is to show how the second crack scheme works. For this end, element 1 in the mesh has an initial crack, as illustrated in Fig. 3c, with a strength half of the intact material. The feldspar material properties and model parameters in Table 1 are used. The simulation results with the constant velocity of 1 mm/s are presented in Fig. 3. Without an initial crack, the model response with this 2element model is linear elastic until the tensile strength (11 MPa) is reached upon which cracks (discontinuities) with a normal into y-direction are introduced into both elements. It should be noted that (red) cracks are plotted at the centroids of the elements for demonstration purposes only; the cracks have no specific location in this model, as discussed in Sect. 2.1. The crack initiation triggers the exponential softening process (described in Sect. 2.1) shown as the blue curve in Fig. 3e. When an initial crack is pre-embedded with the orientation almost parallel to the loading direction, the strength of the numerical specimen, being 31 MPa, is clearly overestimated. This is the crack shielding effect. Had the crack been perfectly parallel to the loading direction, the strength of the specimen would be infinite, as the single crack in element 2 is not enough to cause this mesh to lose its global load bearing capacity. Enabling the second crack to initiate in element 1, according to criterion (12), clearly mends this shortcoming, as observed in Fig. 3e, where the strength of the 2-element model is virtually equal to the case of no initial cracks. However, the softening response is somewhat stiffer than that without initial cracks, which is due to the softening modulus calibration scheme: g d ¼ r Ã t =G Ic through Eq. (13).

Mixed mode I/mode II crack propagation in a pre-notched PMMA specimen
This experiment is used, e.g. by Pham et al. [32] to validate a phase-field model. The PMMA specimen with dimensions in mm is shown in Fig. 4. A very sharp crack, of length 13.47 mm in the present case, is created by impacting a razor blade against the tip of 4 mm wide and 22.86 mm long notch. The specimen is loaded by imposing constant velocity at points A and B at rate 4E-4 mm/s. The thickness of the sample is 3 mm. This problem is used here to validate the fracture model. The computational mesh is shown in Fig. 4. The mesh is refined in the area of interest so that the average mesh size is 0.15 mm therein. The sharp crack is modelled by preembedded displacement discontinuities with (practically) zero strength along the line where the sharp crack is located. Furthermore, it is practically impossible to carry out the numerical experiment at the experimental loading rate due to the explicit time marching of the present approach. Therefore, a loading rate 12 mm/s is applied in The results in Fig. 5 show that the present embedded discontinuity model predicts the experimental crack path with a good accuracy, albeit with slight deviation upon closing to the hole. It is reminded that neither crack path tracking algorithm is used nor crack path continuity over adjacent elements is imposed. Moreover, the force (per thickness of the plate) resisting the crack opening displacement (COD) is 8.8% higher than the experimental (Fig. 5d). It should also be noted that representing the sharp crack in the specimen by pre-embedded discontinuities does not really model a crack tip, as the embedded discontinuity approach has no such enrichment (in contrast to some formulations of XFEM [17]). However, the crack tip enrichment may not always be necessary, as shown by Borja [3] in the context of a sufficiently long edge crack. Furthermore, the present formulation does not explicitly locate the discontinuities inside the elements. There is thus some ambiguity as to interpret the exact initial crack length in the present approach. In view of these aspects of the model, its prediction is surprisingly good. It can thus be concluded that the model is validated with respect to both the ability of the pre-embedded discontinuities to represent real cracks and its performance in predicting the crack path under mixed mode I/mode II loading.

Initial crack population simulations
Now that the model is validated, the main topic of the study, i.e. the effect of the inherent microcrack populations on the tensile strength, can be addressed. The finite element mesh and the boundary conditions are shown in Fig. 6. The average element size is 1 mm. The purpose of the horizontal line the middle height is to attract nodes so that the elements therein have some structuration to be used later. Figure 6 also shows three numerical rock samples with the heterogeneous rock mineral texture described by random clusters of finite elements representing the rock  forming minerals. These are obtained by random permutation (by randperm command of Matlab) of an array of dimensions 1 9 10,980 (number of elements in the mesh) having integers 1, 2, and 3 (coding the minerals), each with the amount of occurrence (in the array) corresponding to the percentage of the minerals in the rock. More precisely, the present rock has 33% of Quartz, which means that this array has 0.33 9 10,980 & 3623 entries of integer 1. As each entry location in this array implicitly codes the global element number in the mesh (e.g., 35 th entry is the global element number 35, which is dictated by the triangulation algorithm of the mesh generator), a spatially heterogeneous representation of rock follows. This simple description, originally introduced by Tang [42], of rock heterogeneity performs well enough [24], when compared to more sophisticated methods based on petrographic image analysis of the rock sample and the mesostructure generation thereof, in the present kind of application, where the research focus is on the sample level behaviour, i.e. the failure mode and strength, instead of grain level behaviour. When the grains and grain boundary effects are of concern, a grain-based or particle modelling approach is required, see, e.g. Wang et al. [48]. Finally, it is reminded that the material properties of the rock forming minerals are given in Table 1.

Uniaxial tension test on intact rock: Effect of mesostructure
The first simulation concerns the uniaxial tension test on intact rock, i.e. in the absence of initial crack populations. Moreover, the effect of mineral mesostructure is tested with the numerical rocks shown in Fig. 6. Figure 7 shows the simulation results for v 0 = 8 mm/s (= 0.1 s -1 ).
The results show that cracks initiate all over the sample but only some of them open enough so as to coalesce and form the final experimental transversal splitting mode [50]. The three numerical rock samples show practically identical stress-strain response, the tensile strength being 9.8 MPa, but the failure modes differ in details and location. The double crack failure mode, where two cracks initiate at the edges of the sample, at different levels of height, and propagate inwards, is realized with NumRock3. The reason for the strikingly similar stress-strain response of the three samples is the relatively small element size which smoothens out the effect of heterogeneity description. Figure 7b shows the cumulative number of cracks initiated during the test on NumRock1. The crack initiation starts when the strain is 1E-4 corresponding to the stress of 7.8 MPa, which is 80% of the tensile strength. The final number of cracks is 1342, which is 12% of the elements.

Uniaxial tension test on rock with a specified initial crack
Before testing the effect of random initial crack populations on the uniaxial tension response of numerical rock, it is instructive to investigate the effect of a single crack of variable length. This is tested here with the model shown in Fig. 6 and the NumRock1. The horizontal crack is introduced as pre-embedded discontinuities with (almost) zero strength located at the middle of the left edge. Simulation results for 1, 5, and 20 mm initial crack lengths are shown in Fig. 8.
All the predicted failure modes attest the general transversal splitting of the numerical samples albeit with differing details, to some extent, despite the fact that the same numerical rock sample was used in each case. This is because in each case of initial crack length (blue in Fig. 8c), the new cracks (red in Fig. 8c) initiate at different location of the numerical rock mesostructure (Fig. 6) and thus the stress state causing the initiation of the new cracks is different in each case. However, they do exhibit some general, mesostructure dictated features, such as the deviation of the crack path from the straight line indicated by the blue arrow in Fig. 8c. The weakening effect of these cracks, when referred to the tensile strength with no initial crack, 9.8 MPa, are 4.1%, 31.6%, and 63.3% for crack lengths 1, 5, and 20 mm, respectively. These percentages exceed the pure geometric cross-sectional area decrease     effect on the nominal stress, which are 2.5%, 12.5% and 50% for these crack lengths.

Uniaxial tension test on rock with random initial crack populations
Next, the effect of randomly oriented initial crack populations always present in natural rocks is simulated. In order to have to worst case scenario, i.e. maximum weakening effect, the pre-embedded cracks are given 1% of the strength of the corresponding finite element in NumRock1. This zero-strength assumption of initial cracks was adopted also by Hamdi et al. [16]. Moreover, the initial crack orientations are assumed to conform to uniform distribution of crack angles (a in Fig. 7c) between 0 and p. Figures 9, 10, 11 and 12 show the simulation results when, respectively, 1, 10, 40 and 80% of the elements in the mesh contain an initial crack. The loading rate is v 0 = 8 mm/s (= 0.1 s -1 ) as above.
The weakening effect with 1% of the elements having a randomly oriented initial crack is 4.1%, which is, coincidentally, the same as in the case of having a 1 mm horizontal initial crack at the left edge of the sample (Sect. 3.3.2). The number of initial cracks per unit area (or crack number density), N A = 0.034 mm -2 (number of cracks/specimen area), is so small that the characteristics of the stress-strain response and the cumulative crack number (which has no contribution from the initial cracks) are similar to the case without initial cracks. Moreover, some extra cracks have been introduced according to the scheme in Sect. 2.3. In most cases, their orientation is almost horizontal but some deviations with almost vertical alignment have occurred (Fig. 9d). The reason for these is the local disturbances in the stress state due to heterogeneity and initial cracks.
When the initial crack number density is 0.34 mm -2 (10% of the elements have a crack), weakening effect is 23.5% (tensile strength = 7.5 MPa). Moreover, the softening part of the stress-strain response (Fig. 10b) is more ductile. The corresponding cumulative number of cracks starts to grow much earlier, in loading, than in the cases of no initial cracks (Fig. 7) and 1% of initial cracks (Fig. 9). This is a realistic feature, as the initial cracks cause stress concentrations thus serving as initiation sites to new cracks. A notable feature is also that the number of new cracks decrease as the initial crack number density increases.
Upon increasing the initial crack number density to 1.37 mm -2 (40% of the elements have a crack), the weakening effect increases to 50% (the tensile strength = 4.9 MPa). Otherwise, the features of the simulation results in Fig. 11 have similar characteristics to those in Fig. 10. A notable difference is that the stress-strain response display quite slow softening as more energy is dissipated in opening cracks.
When the crack number density is 2.75 mm -2 (80% of elements have a crack), the tensile strength of the sample is 2.5 MPa, i.e. weakening effect is 75.5%. The average stress-strain response is highly ductile, which is reflected in the blurred failure mode (Fig. 12a). Moreover, in this extremely cracked case, the number of second cracks, introduced in elements which have unfavourably oriented initial cracks, is larger than the number of new cracks (see Fig. 12d). Without the second crack scheme introduced in this paper, severe crack shielding effect would have taken place in this case.
The effect of initial crack number density on the stiffness modulus and the tensile strength is plotted in a graph in Fig. 13. The stiffness exhibits a linear degradation effect while the tensile strength shows a nonlinear dependency on the crack number density.

Uniaxial tension test on rock with random initial crack populations: effect of crack population
In the simulations above, a single initial crack population was used at each level of crack number density. As their orientations and locations in the numerical rock are random, it is necessary to check what is the effect of the randomness. For this end, the uniaxial tension test is carried out on two additional numerical rock samples having the NumRock1 as their mineral mesostructure and 10% of elements having a randomly oriented initial crack. The simulation results are shown in Fig. 14.
The effect of the random initial crack population is in the details while the general features are similar, as can be observed in Figs. 10 and 14. The tensile strength of the sample at the initial crack number density level of 0.34 mm -2 varies from 7.4 to 7.8 MPa (Fig. 14c) while the stiffness modulus is virtually identical. Therefore, the results in the previous Section are representative.

Dynamic three-point bending of a semicircular disc
Dynamic bending test of a semi-circular disc using the Split Hopkinson Pressure Bar is a relatively recent method to measure the dynamic flexural tensile strength of rocks and other brittle solids proposed by Dai et al. [7]. This test is different from the low-velocity uniaxial tension test addressed above in two respects. First, it is dynamic and, second, the stress state in the specimen is not uniform but has a maximum at the middle of the straight edge, which is the spot for failure initiation. Therefore, the effect of initial crack populations is expectedly different from that in uniaxial tension. The computational model is shown in  [8].
The principle of the computational model, illustrated in Fig. 15a, is as follows. A compressive stress wave, induced by the impact of the striker bar, is simulated as an external stress pulse, r i (t). The incident and transmitted bars are modelled with two-node standard bar elements. The contacts between the bars and the half disc are modelled by contact mechanics. This means that kinematic impenetrability constraints are imposed between the bar end nodes and the half disc nodes at the location of force P 1 and at the supporting pin nodes P 2 /2 in form u bar,xu n,x = b n . Here, u bar,x and u n,x denote the axial degrees of freedom of the bar node and a rock contact node n, respectively, and b n is the distance between the bar end and rock boundary node.
The thickness and diameter of the half disc are L SC-= 16 mm and D SC = 40 mm, respectively. The distance between the supporting pins is 21.8 mm. The length of the incident and the transmitted bar is 1500 mm. The load pulse of form r i (t) = A p sin(xt), where t is time, and x = 2p/T with T = 160 ls, is applied with amplitude A p-= 25 MPa. The loading rate realized in the simulation is measured as the secant of the linear part of the flexural tensile stress vs. time curve. This curve is are calculated from the support reaction forces, P 1 , P 2 , by relation [8] where P is the measured or simulated force (here taken as (P 1 ? P 2 )/2), B and R are the thickness and radius of the half disc, and Y is a function of dimensionless distance S/ 2R with S being the distance between supporting pins. Dai et al. [8] calibrated Y using finite element analysis and found that Y = 5.14 for the present sample dimensions. The value obtained by Eq. (18) is called flexural tensile strength by Dai et al. [8]. Simulation results with the three numerical samples without initial cracks are shown in Fig. 16.
The failure mode in this experiment is the axial splitting of the half disc and the consequent rigid body rotation of the disc quarters about the contact point [4,7,8]. This is predicted by the present method, as attested in Fig. 16b. In a successful design and execution of the experiment, the contact forces should be in balance, i.e. P 1 = P 2 , which is reasonably well achieved in the present simulation (see Fig. 16d). Results in Fig. 16c demonstrate the rock mesostructure has an effect on the details of the macrocrack. Moreover, some deviation occurs in the predicted flexular tensile strength, which is about 25 MPa. This deviation is clearly due to occurrence of different minerals, in each mesostructure, at the location of the maximum tensile stress. As the loading rate in these simulations was 600 GPa/s, the predicted flexular tensile strength is too low, the experimental one being 35 MPa at this loading rate [8]. Actually, the quasi-static flexural tensile strength for Laurentian granite is 25 MPa, while the direct tensile strength for that rock is 12.8 MPa [8]. It should be reminded that the viscosity modulus, which dictates the strain rate sensitivity of the present approach, was set so low that it does not cause any strain rate hardening effects. Moreover, the stress rate of 600 GPa/s translates, using E = 78 GPa, to strain rate of 7.7 s -1 , which is low enough to preclude the influence of inertia effects. In this sense the simulation results represent, to some extent, the quasi-static case and the model prediction is thus acceptable. The effect of initial crack populations is tested next. Figure 17 shows the simulation results for NumRock1 ( Fig. 15c) with initial crack number densities 1.37 mm -2 (10% cracks) and 5.47 mm -2 (40% cracks).
The effect of initial crack number density is substantially weaker here than it was in the uniaxial tension test above. However, the mesh size is different here as well.
Therefore, in order to make the comparison more appropriate, the concept of crack density should be used instead of crack number density as the latter does not take the crack size into account. The crack density is calculated by [14]  where a is the crack half length. However, the crack length does not appear explicitly in the present implementation of the displacement discontinuity model, which is the reason why the notion ''crack number density'' has been used until now in this paper. Nonetheless, as the discontinuity is inside a finite element in the mesh, the finite element size can be used here. Thereby, the crack number density N A-= 0.34 mm -2 in the uniaxial tension simulations corresponds to c = 0.085 (a = 0.5 mm), for which the weakening effect was 23.5%. Now here, the crack number density N A = 1.37 mm -2 gives, using a = 0.25 mm, c = 0.0856 but the weakening effect is only 5.3%. Even at the level N A = 5.47 mm -2 , corresponding to c = 0.34, the weakening effect of initial cracks on the flexural tensile strength is only 12.3%. Repeating the simulation (with NumRock1) for two additional random initial crack populations did not change this result much, as can be observed in Fig. 18. However, with the higher initial crack number density, the failure modes get quite blurred with crack branching.
One reason for the weaker initial crack effect in this problem is the nonuniform stress state in the specimen. Consequently, a favourably oriented crack needs to be at the maximum spot of the tensile stress to have the maximal weakening effect. However, this depends on the random initial crack population. Another reason is that that the element size here is half of the size used in the uniaxial tension simulations, which means that the initial crack length is also half of the that used therein. Therefore, the effect of using the same element size, i.e. 1 mm, is tested as the final simulations. Figure 19 shows the simulation results.
Without initial cracks, the results with the coarser mesh are similar, the flexular tensile strength being about 25 MPa, to those with the finer mesh (Fig. 16). In this respect, the mesh effect is negligible. However, the weakening effect is clearly stronger here as it was with the finer mesh. Indeed, here the weakening effect, at the crack number density level 1.37 mm -2 meaning now c = 0.342 (a = 0.5 mm), ranges from 17.8 to 22.7% (Fig. 19h), when it was only 12.3% with the finer mesh. In uniaxial tension, weakening effect at this level initial of crack density was 50%, reflecting the different nature of the stress state in these two tests.
The crack density concept, while accounting for the crack length, does not account for the stress concentration at the crack tip. This is clearly demonstrated here: at c = 0.34, the finer mesh resulted in weakening effect of about 12%, while the coarser mesh gave a weakening effect of about 20%. This result is in accordance with the fracture mechanics theory saying that longer cracks have stronger influence, as also demonstrated in Sect. 3.3.2.

Concluding remarks
A method to study the effect of inherent microcrack populations on the rock tensile strength was developed in this paper. Present rock fracture model based on the embedded discontinuity finite elements has the computational efficiency of plasticity models. Moreover, the plasticity inspired formulation of the softening model with two discontinuities per finite element effectively prevents the spurious stress locking (crack shielding) of the element upon loading in the direction parallel to the initial crack. Despite the facts that neither crack tracking algorithm was employed, nor crack path continuity was imposed, the method predicts the crack path and the force-COD response in mixed mode I/II propagation with an impressive accuracy.
The simulations of a rock under quasi-static uniaxial tension and dynamic three-point bending of a semicircular disc yield the following conclusions: 1. In uniaxial tension, the weakening effect of initial cracks on the stiffness depends linearly on the crack density while nonlinear dependence on the tensile strength was attested. 2. In uniaxial tension, the higher the initial crack density, the more nonlinear is the pre-peak stress-strain response and the more ductile is the post-peak response. 3. The weakening effect of randomly oriented initial cracks is much weaker in three-point bending than in uniaxial tension due to the non-uniform stress state in the former and the uniform stress state in the latter loading. More precisely, if the average crack length is 1 mm, the weakening effect of initial crack density 0.34 was about 50% in uniaxial tension when it was only about 20% in three-point bending. 4. In accordance with fracture mechanics theory, the model predicted stronger weakening effect with longer initial cracks in three-point bending: at the crack density level 0.34, the weakening effect was about 12% when the crack length was 0.5 mm, while the crack length of 1 mm resulted in 20% weakening effect. A drawback of the present method is that the length of initial cracks depends on the element size. This shortcoming could, however, be mended by developing a technique to extend the initial cracks over the neighbouring elements. Another further urgent development of the method is an extension to 3D case, which should give more realistic results.