Voronoi-Based Discrete Element Analyses to Assess the Influence of the Grain Size and Its Uniformity on the Apparent Fracture Toughness of Notched Rock Specimens

The fracture toughness reflects the rock resistance to crack propagation, and therefore represents an important parameter for rock fracture assessments. From a strict point of view, the real fracture toughness (KIC\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\it {K}_{\mathrm{\it IC}}$$\end{document}) corresponds to a cracked situation in which the notch radius is theoretically equal to zero. However, most of the defects in rocks have a finite radius and, therefore, should be studied as notch-type defects. Here, the notch effect is numerically studied together with the influence of the grain size and the sorting coefficient (grain size uniformity) on the apparent fracture toughness (KIN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\it {K}_{\mathrm{\it IN}}$$\end{document}). To this end, several four-point bending tests with different U-shaped notch radii, mean grain sizes and degrees of uniformity in grain size and shape have been simulated using the Discrete Element Method. In order to represent the grains of the rocks, the Voronoi tessellation is used to create randomly sized and distributed polygonal blocks. These Voronoi polygons have been defined, on the one hand, by an average edge length of 1, 2 and 3 mm, and, on the other hand, by a different number of iterations (n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n$$\end{document}) in the relaxation process during the generation of the polygons, which defines the grain size uniformity. The numerical analyses performed and the interpretation of the results show a clear notch effect in all the studied cases, as the apparent fracture toughness (KIN\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\it {K}_{\mathrm{\it IN}}$$\end{document}) increases with notch radius. Finally, the obtained stress fields at the notch tip have been compared to those obtained from the traditional finite element method. Four-point bending tests with U-shaped notches simulated using the Discrete Element Method. Different notch radii, mean grain sizes and degrees of uniformity are studied. Results show a clear notch effect, increase of apparent fracture toughness with notch radius. Interpretation of the results using the Theory of Critical Distances. A linear relation between the critical distance of the rock and the grain size is observed. The critical distance slightly increases when less uniform grains are studied. Four-point bending tests with U-shaped notches simulated using the Discrete Element Method. Different notch radii, mean grain sizes and degrees of uniformity are studied. Results show a clear notch effect, increase of apparent fracture toughness with notch radius. Interpretation of the results using the Theory of Critical Distances. A linear relation between the critical distance of the rock and the grain size is observed. The critical distance slightly increases when less uniform grains are studied.


Introduction
A comprehensive understanding of rock fracture processes is a major issue of interest in many engineering fields such as civil engineering (e.g., slopes, foundations), underground engineering (e.g., tunneling, mining) or energy engineering (e.g., gas-oil extractions, coal gasification, geothermal energy). The fracture initiation is affected to a great extent by the boundary conditions that are not linked to the rock 1 3 mass or rock matrix itself (e.g., external loads, presence of water). However, other microstructural and macrostructural aspects such as rock composition and mineralogy (e.g., grain size, grain bonding) or the presence of different scale defects also play a key role in the fracture processes (e.g., Hoek 1968;Palmström 1995;Hudson and Harrison 1997;Jaeger et al. 2007).
Defects are generally classified as crack-type defects, those with a theoretical vanishing (or at least negligible) root radius, or as notch-type defects, those with a finite and non-negligible root radius. Rock fracture mechanics (e.g., Whittaker et al. 1992;Aliabadi 1999;Jaeger et al. 2007) traditionally addresses different practical applications such as rock cutting, hydraulic fracturing or underground excavations from a conservative perspective, assuming that the analysed stress risers behave as crack-type defects, based on the use of the conventional Stress Intensity Factor (SIF). However, that approach might be too conservative in many practical situations, since notch-type defects generate less demanding stress fields than crack-like defects and, therefore, develop a higher load-bearing capacity (e.g. , Neuber 1958;Peterson 1959;Pluvinage 1998;Taylor 2007). This is what is generally called the notch effect, which is numerically evaluated in this work. This paper is based on two previous works of the authors (Justo et al. 2017(Justo et al. , 2020b. In the first one, the notch effect of four different rocks was experimentally studied based on the application of the Theory of Critical Distances (TCD) through a parameter called the critical distance ( L ), which is assumed by the TCD as an intrinsic material parameter with length units. The physical meaning of L is still a fundamental challenge among researchers (Taylor 2017), but it is independent of the geometrical features of the stress concentrator and is related to the size of the dominant source of microstructural heterogeneity in the material (e.g., Askes and Susmel 2015). A common source of microstructural heterogeneity in the case of rocks is the grain size (e.g., Taylor 2017). In fact, Justo et al. (2017) confirmed that there is a linear relation between the critical distance of rocks and their mean grain size. Aiming to go deeper into those results, the second work was focused on numerically simulating the influence of grain size on the fracture behaviour of rocks, namely on the main mechanical properties (tensile strength, fracture toughness, Young's modulus and Poisson's ratio) and on the notch effect (through the analysis of the apparent fracture toughness). To do so, several discrete numerical analyses were performed to simulate unconfined compression tests, Brazilian tests, and four-point bending tests (as those carried out in the laboratory by Justo et al. 2017) with variable mean grain sizes. Particles were defined by means of Voronoi tessellations with relatively uniform shape in all the cases (with a sorting coefficient very close to unity), representing polygonal grains and ideal non-porous, crystalline and isotropic rocks. The modelled materials were ideal zero porosity crystalline rocks, whose parameters have the same order of magnitude as those of the Macael marble tested in Justo et al. (2017). To make the computational cost feasible, the modelled grains (1-3 mm) were much larger than those of the Macael marble (average grain size of 335 µm), and consequently, the intention was not to reproduce the same behaviour observed in the laboratory but to represent models of rock-like materials within a realistic order of magnitude.
Here, it is intended to extend those works and consider not only the influence of the mean grain size but also the grain size distribution (i.e., sorting coefficient) as a variable. Following the methodology used by Justo et al. (2020b), additional discrete numerical models have been constructed in this work to simulate the same type of Single Edge Notched Bend (SENB) specimens with variable notch radii and subjected to four-point bending conditions (i.e., mode I loading conditions), keeping the same mean grain sizes constant but varying the degree of uniformity of the grains (i.e., increasing the sorting coefficient). These models represent ideal non-porous, crystalline and isotropic rocks with different grain sizes and sorting coefficients, and allow to obtain, for each case, the stress state in the surrounding of the notches at the onset of crack initiation and propagation. The stresses derived from the Discrete Element Method (DEM) are interpreted according to the TCD, which is used to analyse the variation of the apparent fracture toughness with the notch radius. Finally, those DEM stress fields used for the analyses of the notch effect are compared to those obtained by the finite element method (FEM), to compare discrete and continuum approaches and to validate the considered methodology.
A limitation of this study is that only intergranular fracture is considered. This may be valid or representative for some specific rocks or materials, where intergranular fracture is the only fracture micromechanism (e.g., Ortiz and Suresh 1993). However, in many other rocks, such as in the marble used here as a reference, intergranular fracture coexists with transgranular fracture. If transgranular fracture was simulated (e.g., Hofmann et al. 2015;Peng and Wong 2017), the fracture toughness would be reduced, because in some cases, the fracture would occur earlier through the grains. Besides, the notch effect would also be reduced, because the presence of a grain at the notch tip would not restrict fracture as much as if the grains were unbreakable, as occurs in this paper. Finally, Li et al. (2019Li et al. ( , 2021) present a detailed discussion of the methods to consider transgranular fractures and their influence. For example, finite-discrete element methods (FDEM) allow for transgranular fracturing and also for advanced contact algorithms (e.g., Zhao et al. 2018;Wang et al. 2021a). Likewise, in DEM, multi-scale tessellations can be applied to simulate transgranular fracturing (e.g., Wang and Cai 2018). A novel technique, such as the grain-based finite-discrete element methods (GB-FDEM) used by Abdelaziz et al. (2018) or Li et al. (2020), would be suitable to consider transgranular fractures in this problem and extend the present work.

Background
The notch effect in rocks, especially in crystalline rocks, depends on the grain size as a consequence of its relationship with L . For example, to evaluate when the notch root radius is negligible, this is compared with L , since it has been observed for different materials that the notch effect is generally negligible when the notch radius is smaller than L . Using laboratory tests to study the influence of the grain size on the notch effect is extremely difficult, because rocks have many different features, such as inhomogeneities, complex grain size distributions, grain aspect ratio, porosity, etc. Numerical analyses, and in particular the DEM, are available suitable tools at present for an in-depth analysis of this problem.
Block-based DEMs allow the rocks to be modelled as an assemblage of blocks (grains) with different boundary conditions, which makes them an appropriate tool for the problem addressed in this work. In the literature, four methods are typically used to simulate grain structure, namely disk-shaped grains (e.g., Potyondy and Cundall 2004), square-shaped grains (e.g., Li and Konietzky 2015), triangular grains (e.g., Kazerani 2013) and polygonal grains (e.g., Kazerani and Zhao 2010). The polygonal grain structure appears to be a more realistic representation of the microstructure of rock-type materials (especially crystalline rocks). The conventional polygonal structure is usually generated using the Voronoi tessellation technique (e.g., Gao et al. 2016). Here, the grains are generated using the Voronoi tessellation, which allows a random distribution of the grains to be defined with a controlled mean size. This technique provides blocks with similar shapes to the grains of zero porosity crystalline rocks, similar to those observed by the authors for different marbles (e.g., Justo et al. 2017). Voronoi-based discrete models have been successfully used to simulate rock-like materials under different considerations, for instance, for compression, Brazilian and fracture toughness tests (Chen et al. 2015). The influence of grain size, pore size and mineral composition have also been studied using this technique (e.g., Li et al. 2017a;Chung et al. 2019;Liang et al. 2021). Voronoi tessellations were first extended to 3D models by Ghazvinian et al. (2014), who simulated crack damage development in brittle rocks. Thus, this technique has proven to be appropriate for the problem under study. Finally, Lisjak and Grasselli (2014) and Zhang and Wong (2018) present reviews of these numerical techniques and the introductions of Li et al. (2019) and Wang et al. (2021b) summarise some of the latest numerical advances.
Different criteria can be found in the literature to predict fracture loads of notched components subjected to mode I loading (e.g., Kipp and Sih 1975;Carpinteri 1987;Seweryn 1994;Gómez et al. 2000;Lazzarin and Zambardi 2001;Yosibash et al. 2004;Taylor 2004) or even mixed mode loading conditions (e.g., Papadopoulos and Paniridis 1988;Seweryn and Mróz 1995;Yosibash et al. 2006;Berto et al. 2007). Among them, the most widely used criteria at the moment are probably the following ones: • The Cohesive Zone Model (CZM), first proposed by Barenblatt (1959) and Dugdale (1960) to describe the stress fields and fracture processes near the defect tip. • Finite Fracture Mechanics (FFM), based on the assumption that the crack grows by finite steps determined by a condition of consistency of both energy and stress requirements (Carpinteri et al. 2008). • The Strain Energy Density (SED) criterion, an energybased approach that combines the elementary volume proposed by Neuber (1958) and the local Mode I concept of Erdogan and Sih (1963). • The Theory of Critical Distances (TCD), first proposed by Neuber (1958) and Peterson (1959), although it has not been until the last decades with the development of finite element analyses that this methodology has been scientifically applied and further developed (e.g., Taylor and Wang 2000;Taylor 2001Taylor , 2007Susmel and Taylor 2003;Cicero et al. 2012).
One of the greatest advantages of the TCD consists of the possibility of obtaining (semi-) analytical results to correctly perform static assessments without any loss of accuracy (Cornetti et al. 2016), and without the need for significant computational efforts as commonly required by the SED criterion or by the CZM, for example. However, despite its simplicity and potential for the analysis of fracture processes as demonstrated for a wide range of materials (e.g., metals, ceramics, polymers and composites), scarce work can be found on the application of the TCD in the particular field of rock mechanics. Some examples of the application of the TCD on rock mechanics are, for instance, the works of Lajtai (1972), Ito and Hayasi (1991), Ito (2008) or Schwartzkopff et al. (2017). The authors have recently applied the TCD for the fracture assessment of several rocks under different temperature and loading conditions (Justo et al. 2017(Justo et al. , 2020a(Justo et al. , 2021, and given the successful results, the evaluation of the notch effect in this work is also performed using this local failure criterion. All the failure criteria indicated above, including the TCD, use the fracture toughness ( K IC ) for the fracture assessment of structural components. Indeed, this parameter represents the residual strength of a cracked component to crack propagation. By definition, the fracture toughness is an intrinsic property of the material. It addresses crack problems where the notch radius can be assumed to be equal to zero and, therefore, where no notch effect is deployed. However, when notch-type defects are analysed, an apparent fracture toughness ( K IN ) should be considered instead of the strictly real value of K IC , as demonstrated, for example, by the authors in previous works (Cicero et al. 2014;Justo et al. 2017). With this, the interpretation of the notch effect is performed in this work by analysing the variation of the apparent fracture toughness with the notch radius.
For the correct application of the TCD, the stresses around the notch tip must be assessed. In the particular cases of mode I loading conditions, as those considered in this work, the stresses are evaluated along the bisector plane of the notch tip and, to do so, numerical methods provide a suitable tool. The different geological processes that rocks have undergone over millennia makes them highly heterogeneous materials. However, despite this generally accepted condition, it is a common practice to consider the rock as a continuum in order to simplify the rock mechanics analyses. In fact, applications of the TCD found in the literature are generally based on continuum stress assessments regardless of the analysed material, even for rock-type materials. This continuum approach might be a suitable option when global responses are of interest (e.g., Eberhardt et al. 2004), but it is less appropriate for applications where detailed information is required for more accurate assessments. For example, Gui et al. (2016), investigated the grain size effect simulating Brazilian disks and unconfined compression tests using distinct element analyses, and they reported that larger particle sizes produce a higher stiffness and strength of the intact rock. Li et al. (2017b) also used distinct element analyses to simulate crack initiation and propagation of a granular rock. Liu et al. (2018) numerically studied the so-called Fracture Process Zone (FPZ) by means of the DEM, addressing size effect and particle size. Likewise, Wang et al. (2019) studied the influence of mineral heterogeneity on the thermomechanical behaviour of rocks also using the DEM. Considering all of the above, it is clear that detailed microstructural analyses require discrete numerical approaches.

Numerical Analyses
Two types of numerical analyses have been considered in this work. The first one uses the DEM to simulate the rock specimens as an assemblage of grains. Idealized crystalline and non-porous rocks are modelled, in which grain size and uniformity are the analysed variables. This microstructure could correspond to a marble, such as the Macael marble analysed by Justo et al. (2017). In fact, this work takes as a reference the experimental results obtained in Justo et al. (2017), for the definition of the parameters of the models. However, it is not the purpose of this research to simulate the exact behaviour of the previously analysed Macael marble (which has a relatively smaller grain size than the studied models), but to model ideal rock-like materials with comparable and realistic properties.
The proposed ideal rock simulations allow the effect of the grain size and the degree of uniformity to be studied in isolation, without the influence of other possible variables. Thus, the interpretation of the results becomes straightforward and generalizable conclusions can be reported. The second approach is based on the use of the FEM and aims to compare the stress fields around the notches and, in particular, at the bisector plane, where stresses are assessed according to the TCD.

Discontinuum Approach
In this work, the Universal Distinct Element Code, UDEC v6.00 (Itasca 2010), is used to numerically study the influence of the grain size and its uniformity on the apparent fracture toughness of notched rock specimens. This code refers to a particular DEM scheme that uses deformable contacts and an explicit, time-domain solution. It consists of a blockbased method that models the rock masses as an assembly of blocks with interfaces (or contacts), allowing the simulations of rock fracturing. The rock behaviour is described by the interaction between these blocks. Thus, the contact laws govern the macroscopic response of the material according to normal and tangential behaviour at the interfaces. The DEM allows finite displacements and rotations of discrete bodies (including complete detachment), and considers the interaction of the fractured rock fragments.
To represent the grains of the rock, the Voronoi tessellation is used to generate randomly sized polygonal blocks. These polygons are defined by an average edge length ( l ), which has been varied in this work to consider the effect of the grain size. Likewise, the blocks (grains) are discretised into constant-strain finite-difference triangular zones defined by a maximum edge length ( e ). These zones make the blocks deformable. Here, e = l has been considered to keep the proportion of zones within the blocks roughly constant and, on average, each grain is approximately subdivided into seven zones in all the cases. The authors have not performed mesh sensitivity analyses of the influence of the parameter e , information about the influence of the l∕e ratio may be found in Fabjan et al. (2015). At the same time, the size and shape uniformity of the Voronoi polygons is defined by the number of iterations ( n ) in the relaxation process during the generation of the mesh (Itasca 2010), in such a way that the increasing number of iterations leads to Voronoi polygons of more uniform size and shape. n has been set to 30, to provide a relatively uniform distribution of the grain size, and also to 1 to provide less uniform grains, as observed in Fig. 1. This figure represents the considered Voronoi tessellations with average edge lengths ( l ) of 1, 2 and 3 mm, as well as with different degrees of uniformity by considering n = 30 and n = 1.
For a detailed description of the modelled grain conditions, Fig. 2 provides the size distribution curves for all the considered combinations of l and n , in which the grain diameter corresponds to the diameter of ideal circular grains with the same area as the actual modelled grains. Both frequency and cumulative frequency are represented in percentage.
It is observed that those grains modelled with n = 30 present a very narrow variation of the grain size, i.e. nearly uniform distribution. In contrast, those curves in Fig. 2 corresponding to n = 1 keep the same mean grain size as in the models with n = 30 but present a larger variation of the grain size. With this, Table 1 summarises some statistical values of the modelled grains, namely the mean grain diameters (corresponding to the diameters of ideal circular grains with an equivalent area), their standard deviations, the first quartile (Q1), the third quartile (Q3) and the sorting coefficients defined as the ratio between the first and third quartiles (Q1/Q3). This latter parameter indicates that the size distribution is relatively uniform for n = 30 as the coefficient is close to 1. In contrast, as the sorting coefficient increases for n = 1, the size distribution of the grains is less uniform.
Regarding the definition of the constitutive models, the grain behaviour has been defined in this work by a linearly elastic isotropic model, using the parameters indicated in Table 2.
On the other hand, the grain boundaries or contacts have been defined by the Coulomb slip model with residual strength. This model provides a linear representation of joint stiffness and yield limit, and is based upon elastic stiffness, frictional, cohesive and tensile strength properties, and dilation characteristics common to rock joints (Itasca 2010). This residual strength simulates the displacement-weakening of the joint by loss of frictional, cohesive and/or tensile strength at the onset of shear or tensile failure. That is, when a joint is fractured, the joint tensile strength, the joint friction angle and the joint cohesion are set to residual values. Table 3 gathers the used parameters (Justo et al. 2020b) for the Coulomb slip constitutive model and the residual values, as well as the joint normal and shear stiffnesses, which are zone size dependent (Itasca 2010). It should be highlighted that the ratio between the cohesion and the tensile strength values indicated in Table 3 looks quite unusual and would be wrong considering classical constitutive laws in continuum mechanics, where cohesion comprises both adhesion between particles and interlocking. However, since polyhedral Voronoi elements and only intergranular fracture are being considered, the interlocking is very strong. Thus, the contact cohesion has been reduced to get realistic strength values. The tensile strength test models reported in Justo et al. (2020b) show only tensile failure along the grain boundaries, leading to a single vertical macro-fracture. In contrast, the uniaxial compression test models showed tensile fracturing, but also some shear fracturing in the Fig. 1 Representation of the Voronoi tessellations for an average edge length ( l ) equal to 1, 2 and 3 mm and for n equal to 30 and 1 boundaries, which indicates that the parameters used lead to realistic failure modes.
The definition of these parameters is clarified by the authors in a previous work (Justo et al. 2020b). They correspond to an ideal zero porosity crystalline rock that loosely resembles a Macael marble tested by the authors (Justo et al. 2017). Since the modelled grains (1-3 mm) are much larger than those of the Macael marble (average grain size of 335 µm), a detailed calibration was not considered necessary and only a reasonable macroscopic response of the same order of magnitude as that of the marble was searched.   Table 4 shows, for comparison purposes, the uniaxial compression strength ( c ) and tensile strength ( u ) values of the Macael marble together with those obtained from numerical models. In particular, the results of the uniaxial compression test, Brazilian test and direct tensile test models with an average edge length ( l ) of 1, 2 and 3 mm (and n = 30) are included in this table. The emergent values of c are slightly lower than those of the Macael marble, while the emergent values u of the Brazilian test models are slightly higher than those of the laboratory tests. No direct tensile tests were performed in the laboratory for the Macael marble, so the emergent values of u obtained from the direct tensile test models cannot be directly compared. However, the tensile strength of the Macael marble was also characterised by the authors (Justo 2020) by means of four-point bending and three point bending tests, obtaining u values of 17.53 MPa and 13.94 MPa, respectively. Thus, although the tensile strength values derived from the direct tensile test models are relatively high (possibly because of the selected contact tensile strength and cohesion values indicated in Table 3), they are still within a realistic order of magnitude.
Based on the aforementioned meshing criterion and constitutive models, several four-point bending tests have been numerically simulated. These models consist of 180 × 30 mm size specimens under plane strain conditions, with variable U-shaped notch radii ( ) of 1, 2, 3, 4, 7, 10 and 15 mm to assess the notch effect, all of them having a notch length equal to half of the height (i.e., 15 mm). The geometry of the models is depicted in Fig. 3. As observed in the figure, the Voronoi tessellation has only been generated within a specified range around the notch, where the fracture is assumed to start. The rest of the model is analysed by means of the Finite Difference Method (FDM) using the same material constitutive model for the whole geometry. For the sake of simplicity, the same values of the bulk ( K ) and shear ( G ) moduli have been used in the continuum domain as for the grain-based domain. From a strict point of view, the elastic parameters ( K and G ) of the continuum domain should be the emergent properties of the Grain-Based Models (GBM). Since these are not initially known, the same values as those of the grain-based domain were used. The models were not recalculated with the emergent values of K and G, because the differences were not large and it was assumed that, for isostatic problems as those studied in this work, their influence is small. Besides, different sizes have been checked for the DEM region, providing similar stress fields around the notch tip. For this reason, the Voronoi tessellation has only been applied within a 30 × 30 mm size region as shown in Fig. 3. The objective of simplifying the construction of the models by only introducing the Voronoi tessellation at the center of the specimen models was to reduce the computational cost. In this manner, the models are computationally efficient without a significant influence on the results.
The random arrangement of the Voronoi polygons (grains) produces a scatter of the results, as in reality. For this reason, a repetitiveness of six models has been considered for each mean grain size, degree of uniformity and notch radius combination, only varying the randomly generated Voronoi blocks.   The simulations have been executed under displacement control, specifying a constant vertical velocity at the upper loading points (Fig. 3) sufficiently small to minimise shocks to the system, and zero vertical velocity conditions at the supports. Thus, quasi-static calculations have been performed under mode I loading. Considering these boundary conditions, the failure load F (with force/depth length units) has been calculated for each of the individual numerical models (see the scheme in Fig. 4). Considering the analysed specimens as beams, the bending moment ( M ) between the loading points is constant in a four-point bending configuration as depicted in Fig. 4: where s is the span between the outer supports. In parallel, the horizontal stress law along the bisector of the notch can be equated to a pair of forces ( f ) with a lever arm ( z ), or in other words, to a bending moment ( M ). This bending moment is related to that in Eq. (1) as follows: Thus, the failure load ( F ) is derived from the bending moment generated at the bisector of the notch tip just in the calculation step prior to the appearance of the first crack (Justo et al. 2020b): The obtained F values are collected in Table 5 for the different notch radii ( ), grain sizes ( l ) and degrees of uniformity ( n).
These models assume that the cracks will propagate only along the grain boundaries. Thus, only intergranular failure is being considered. As a consequence, the minimum radius of the aforementioned notches has been limited to the grain size in each case ( ≥ l ). This restriction is intended to avoid the possibility of a notch not having at least one grain boundary at the tip, which could cause the breakdown of the calculation model. Figure 5 shows, as an example, some representative numerical models for each grain size after the appearance of the cracks at the notch tip, once the failure load has been exceeded.

Continuum Approach
The influence of the grain size cannot be directly assessed using FEM analyses, since the materials are modelled as a continuum instead of as an assemblage of blocks. As a consequence, the failure criterion established for the discrete numerical models (crack initiation when shear or tensile strength of the grain contacts is reached) cannot be imposed in a continuum approach. However, in this work, the stress fields of both methods are compared, aiming to analyse the extent of the differences between the two types of analyses.
In this case, a finite element code called PLAXIS 2D (2017), has been used to model the test specimens. The simulated rock samples correspond to the same geometry as that shown in Fig. 3, including the same notch radii. The finite element mesh has been refined in the region surrounding the notches as depicted in Fig. 6, with a view to avoiding any possible influence of the mesh. For this reason, no repetitiveness of the models has been considered in this case, as the influence of the mesh is considered to be negligible.
On the other hand, a linear elastic constitutive model has been used to simulate the rock and, therefore, only two parameters are required: the Young's modulus ( E ) and the Poisson's ratio ( v ). These deformational parameters correspond to the macro-scale behaviour of the rocks. Thus, they cannot be derived from the parameters in Table 2 that define the behaviour of the grains. The parameters used in the FEM models have been obtained from the unconfined compression test simulated by the authors in a previous work (Justo et al. 2020b) by means of DEM, considering the same idealised non-porous crystalline rock material with different grain sizes. These parameters are collected in Table 6 and correspond to the emergent elastic properties for each of the analysed grain sizes.
With regards to the boundary conditions of the FEM models: all the contours have been set free except for the supporting points and loading points. In this case, the simulations are not performed under displacement control. Instead, to obtain comparable stress fields to those of the DEM approach, the failure loads derived from the DEM models (Table 5) have been directly introduced in the

Analytical Interpretation
This research studies the fracture of U-notched rock specimens. To this end, the fracture analysis is equated to a situation in a cracked component where the apparent fracture   Srawley and Gross (1976), for Single Edge Notched Bend (SENB) specimens, as those modelled in this work and shown in Fig. 3: where F is the failure load of each of the simulated models (Table 5), h is the specimen height ( h = 30 mm) and Y is a compliance non-dimensional factor that only depends on the geometry of the specimen and is given by the following expression: Fig. 5 Examples of formed cracks in DEM models for an average edge length ( l ) equal to 1, 2 and 3 mm and for n equal to 30 and 1 Fig. 6 Representation of the simulated four-point bending test finite numerical models s o and s i represent the spans between the outer supporting rollers (150 mm) and the inner loading points (50 mm), respectively, and 0 is the relative crack length defined as the ratio between the initial notch length (15 mm) and the total height (30 mm) of the specimen ( 0 = 0.5 ). With all this, given the four-point bending test configuration and geometry depicted in Fig. 3, the Y factor is equal to 10.16. The analytical interpretation of these numerically obtained results is performed using the TCD, and more specifically the Line Method (LM). A detailed description of this methodology is presented by Taylor (2007). The LM is a local failure criterion based on the stress law over a certain distance ( d ) from the notch tip. It states that failure occurs when the average stress over d is equal to the inherent strength ( 0 ) of the material, as represented in Fig. 7.
The distance d is related to a material characteristic parameter called the critical distance ( L ) which, in the case of rocks, is in the order of a few millimeters (e.g., Cicero et al. 2014). It can easily be demonstrated that d = 2L (e.g., Taylor 2007). Thus, the failure criterion can be written as follows when the LM is considered: Based on this expression and considering the stress distribution function along the bisector of the notch tip proposed in the past by Creager and Paris (1967), the LM of the TCD provides the following analytical solution for the assessment of the apparent fracture toughness: This expression is used to assess the notch effect (i.e. the variation of K IN with the notch radius) and the influence of (6) the grain size and uniformity on it, using as a reference the numerical results derived from Eq. (4). Variation of the apparent fracture toughness with the notch radius for different grain sizes ( l ) of 1 mm, 2 mm and 3 mm, for the particular case of n = 1 1 3

Notch Effect
The notch effect is graphically shown in Figs. 8 and 9 documenting the variation of the apparent fracture toughness with the notch radius for n = 30 (Fig. 8), representing a relatively uniform grain structure, and n = 1 (Fig. 9), for less uniform grain structures. Each of the dots in the graphs corresponds to the individual results of the apparent fracture toughness ( K IN ) calculated with Eq. (4), while the dashed (Fig. 8) and solid ( Fig. 9) lines stand for the best-fit curves according to Eq. (8), which corresponds to the LM of the TCD. Three different curves are depicted in both graphs, for grain sizes of 1 mm (A), 2 mm (B) and 3 mm (C). The statistical method considered to fit the curves to the data points is the least squares method in all the cases, leaving fracture toughness ( K IC ) and the critical distance ( L ) as free variables.
In all the cases, the notch effect looks clear, since the apparent fracture toughness shows a continuous increment with the notch radius. This increment seems to lessen as the notch radius increases. However, the notch size after which the notch effect is negligible has not been captured with the analysed range of notches. In parallel, the notch effect is also supposed to be negligible below a certain notch radius as stated by Taylor (2017). However, this is neither observed in Figs. 8 and 9. The physical reason explaining this phenomenon is probably related to the stress concentration region around the notch tip. The smallest notches develop the highest stress concentration at the notch tip, probably affecting only a small number of grains, while the stresses in the case of the largest notches are not so concentrated just at the tip and are somehow distributed among more grains in the surrounding of the notch tip. Consequently, if the notch radius is sufficiently small, differences in stress concentrations would probably be insignificant compared to the scatter of the numerical or experimental results and, in the opposite situation of a theoretically infinite notch radius, the notch concentration would be null and would behave as a section reduction rather than as a stress riser.
Likewise, a clear displacement of the three curves is directly observed when analysing both Figs. 8 and 9. Similar trends are observed in both cases. The curves move upwards as the grain size increases, which implies an increment of the fracture toughness. On the other hand, the curves slightly flatten when the grain size increases, which seems to indicate that the notch effect is less significant when the grain size is larger.
With this, Table 7 gathers the corresponding results of both K IC and L variables derived from the adjustment of Eq. (8) for each of the analysed combinations of l and n . The coefficient of determination ( R 2 ) is also given in Table 7 to specify the goodness of fit of each of the curves with respect to the mean numerical values. It is observed that, in general terms, the value of the fracture toughness ( K IC ) is roughly constant for both values of the parameter n . Thus, K IC seems to be influenced only by the average grain size in statistical terms. On the other hand, the value of the critical distance ( L ) presents a certain decrease when moving from n = 30 to n = 1, or in other words, it reduces when grain size is less uniform. These trends are consistent for the three analysed grain sizes with l equal to 1, 2 and 3 mm. This is easily observed in Fig. 10, which gathers in a single plot the best-fit curves of the LM of the TCD represented in Figs. 8 and 9. It is observed that the three curves A, B and C tend to move downwards and to slightly flatten for more uniform grains  . 10 Comparison between the best-fit curves of the variation of the apparent fracture toughness with the notch radius represented in Figs. 8 and 9 for different grain sizes ( l ) of 1 mm, 2 mm and 3 mm, for both n = 30 and n = 1 (i.e., n = 30), which gives as a result the values of K IC and L collected in Table 7. Similar responses have been observed in previous experimental studies as those performed by Justo et al. (2017), in four different rocks. Figure 11 represents the case of the Macael marble as an example. The dots correspond to the individual test results while the dashed line indicates the best-fit curve according to Eq. (8), the same as in Figs. 8, 9 and 10. Although it is not the purpose of this work to reproduce the same behaviour, both the experimental and the numerical results offer similar trends within a relatively similar order of magnitude.
Finally, comparing in Fig. 12, the values of the critical distances and the mean grain size, there seems to be a linear relation among them. This relation should be considered in qualitative terms, since the performed numerical simulations are an idealisation of the real problem. For example, 3D effects are neglected here, as plane strain conditions are being considered. Besides, only intergranular fractures are considered. Taylor (2017) related the critical distance of different materials with clearly distinguishable microstructural distances such as the grain size in the case of rocks. He reported that in most of the cases, L lies between 1 and 10 times this distance. Although an accurate definition of the relation between L and the grain size still requires further research, this research shows that the correlation exists for zero porosity crystalline rocks. This relation could allow simplified finite element analyses to be performed for the fracture assessment of rocks considering the influence the grain size through the critical distance.

Comparison of the Stress Fields
According to the TCD, the stresses are evaluated along the bisector plane of the notch tip in the case of mode I loading conditions, as those analysed in this work. For this reason, the stresses normal to this plane corresponding to both the DEM and FEM analyses are compared here. Figure 13 gathers some representative curves for both approaches, including the stress laws for = 3 mm (Fig. 13a), = 7 mm (Fig. 13b) and = 15 mm (Fig. 13c), all of them for a mean grain size of 1 mm and 3 mm and for the particular case of n = 30. In general terms, good agreement is observed between the solid curves corresponding to the continuum approach and the dotted curves corresponding to the discrete approach. The illustrated DEM curves correspond to a particular random case and can vary in each model depending on the actual position of the grains. Besides, 20 history points were considered along the bisector plane to represent the latter DEM curves in all the cases.
If the stresses of FEM and DEM results are compared for r = 0 mm, in general, the maximum stress at the notch tip seems to be slightly higher when FEM is used, although the difference is practically negligible when the grain size is sufficiently small. Besides, the stepped appearance of the dotted curves softens as the grain size decreases. In fact, in those curves corresponding to the models with a grain size of 1 mm (Fig. 13a1, b1, c1), the differences between the DEM and FEM curves are insignificant.
The observed influence of the grain size close to the notch is somehow absorbed by the LM of the TCD, which evaluates the stresses along a certain distance ( 2L ) from the notch tip instead of considering the maximum stress at the tip. This distance is assumed to be sufficiently small not to be affected by the boundary of the model. According to the critical distance values obtained from the numerical models and gathered in Table 7, this hypothesis is valid.
The consequences of the notch effect can be directly observed in the plots of Fig. 13. The stress concentrations around the notch tip are relatively higher for the smallest notch radii, as expected. When the notch radius is sufficiently large (e.g., 15 mm), the notch effect tends to vanish and the  Fig. 14 Comparison of the numerically obtained DEM and FEM stress contours ( xx ) around the notch tips for ρ = 3 mm, ρ = 7 mm and ρ = 15 mm, both for grain sizes of 1 mm and 3 mm and for the particular case of n = 30 stresses approximate to those corresponding to a simple section reduction with no appreciable stress intensification.
Finally, aiming to support the results observed in Fig. 13, Fig. 14 represents the horizontal stress contours near the notch tip considering both the FEM and the DEM models. These stresses correspond to the moment just prior to failure and, as in Fig. 13, they stand for the same representative cases, i.e., for = 3 mm, = 7 mm and = 15 mm. All of them are for a mean grain size of 1 mm and 3 mm and for the particular case of n = 30. In consistence with Fig. 13, good agreement is observed between the stresses of the continuum and discontinuum approaches, not only along the bisector planes but also throughout the surroundings of the notches. The influence of the presence of grains on the stress contours is more pronounced for the largest grains (i.e., 3 mm), which confirms the staggered shape of the curves (see Fig. 13a2, b2 and c2). In contrast, for the smallest analysed grains (i.e., 1 mm), the stress variation is appreciably smoother and very close to that obtained with the continuum approach.

Conclusions
This work studies the influence of the grain size and its uniformity on the apparent fracture toughness ( K IN ) of an idealised non-porous crystalline rock. To this end, different four-point bending tests with variable notch radii, grain sizes and degrees of uniformity have been numerically studied, using block-based numerical models based on DEM. The main limitation of these models is the fact that only intergranular fracture is considered, not accounting for the possible coexistence between intergranular and transgranular fracturing. The interpretation of the results is based on the use of the LM of the TCD, which analyses the stress fields along the bisector of the notch tip. Finally, these stresses have been compared with those obtained from a traditional continuum approach.
Based on the obtained results, the following conclusions should be highlighted: • The notch effect is clear for the range of analysed notch radii regardless of the grain size or uniformity, as the apparent fracture toughness increases with the notch radius in all the cases. • The best-fit curves of K IN corresponding to the LM of the TCD show that the fracture toughness increases with the grain size. Additionally, it seems that the notch effect is relatively softened with the grain size as well, as the curves slightly flatten when the grain size increases. • The fracture toughness does not vary with the degree of uniformity, if the mean grain size is the same. On the other hand, the critical distance slightly decreases when less uniform grains are analysed. • There seems to be a linear relation between the critical distance of the rock and the grain size. However, the obtained relation should be understood in qualitative terms, as the problem studied corresponds to an idealised situation in which the 3D effect is neglected and only intergranular fracture is considered. • DEM and FEM provide similar stress fields at the notch tip when the grain size is sufficiently small. The staggered stress curves associated to the discrete approach soften as the grain size decreases. In general terms, good agreement between the two types of models have been obtained for the studied range of notch radii and grain sizes.
With all this, the TCD provides a suitable tool for the fracture assessment of rocks with notch-type defects. Likewise, discrete numerical approaches allow detailed analyses to be performed when precise information is required, as the influence of the grain size and its uniformity for example. However, in general terms, finite element analyses provide accurate approximations when the global behaviour is studied and the influence of microstructural aspects (e.g., grain size) can be neglected.