Non-local criteria for the borehole problem: Gradient Elasticity versus Finite Fracture Mechanics

Two nonlocal approaches are applied to the borehole geometry, herein simply modelled as a circular hole in an infinite elastic medium, subjected to remote biaxial loading and/or internal pressure. The former approach lies within the framework of Gradient Elasticity (GE). Its characteristic is nonlocal in the elastic material behaviour and local in the failure criterion, hence simply related to the stress concentration factor. The latter approach is the Finite Fracture Mechanics (FFM), a well-consolidated model within the framework of brittle fracture. Its characteristic is local in the elastic material behaviour and non-local in the fracture criterion, since crack onset occurs when two (stress and energy) conditions in front of the stress concentration point are simultaneously met. Although the two approaches have a completely different origin, they present some similarities, both involving a characteristic length. Notably, they lead to almost identical critical load predictions as far as the two internal lengths are properly related. A comparison with experimental data available in the literature is also provided.


Introduction
Engineers and scientists have used classical continuum solid mechanics along with strength-based failure criteria since the middle of the nineteenth century, usually successfully. However, because of the lack of an internal length, this approach fails in predicting sizes effects, i.e. the change of the mechanical response when the spatial dimensions are scaled up or down, while the shape and all other characteristics are preserved. Understanding size effects is a problem of paramount importance when moving across the scales, e.g. going from laboratory specimen size to real size structures. Moreover, quoting Bažant [1], ''Scaling is a quintessential problem of every physical theory. If the scaling is not understood, a viable theory does not exist''.
In the twentieth century, the development of Linear Elastic Fracture Mechanics (LEFM), i.e. the theory based on Griffith infinitesimal energy balance during crack growth, allowed a first insight in the size effect explanation. Roughly speaking, LEFM is able to predict size effects because of the interplay between the strain energy stored in the structure, proportional to the volume, and the dissipated energy, proportional to the fracture surface. However, LEFM can be applied only to pre-cracked structures.
Among the different theories developed in the last decades to address the size effect issue (e.g. the fractal approach [2]), let us mention two main research directions, both falling within the framework of nonlocal models. The former one is non-local elasticity: the stress at a point depends on the strain in a region surrounding that point (strong non-locality or integral models) or on the strain at that point plus its derivatives (weak non-locality or gradient models). The latter one, retaining the classical linear elastic constitutive law, deals with non-local fracture criteria. The presence of an internal length characterizes both models: thus, they are able to detect size effects as well as to highlight phenomenologically the influence of the material microstructure at the macroscale.
Non-local and higher order elasticity theories date back to Eringen et al. [3] and Mindlin [4]. Accordingly, the stress state at a point depends also on the deformation state at neighbouring points [5,6]. Analogously, the first (or the second, etc.) gradient of the strain field may be considered as additional terms in the strain energy expression because they characterize the spatial differences in the elastic strain at neighbouring points. By including such terms, new modified constitutive equations appear and the relevant Gradient Elasticity (GE) theories are developed (see e.g. [7][8][9][10][11][12]). In other words, classical elasticity does not present any intrinsic length scale because it incorporates only the nearest neighbour interaction; an intrinsic length scale appears when the forces between particles are extended to include first, second, until n-th neighbour interactions. Thus, GE may be considered as a higher-order approximation of a fully non-local theory [9,13]. The addition of higher order terms in the governing equations takes into account (i) the nonlocality of the physical state of the system in an approximate way, and (ii) the physical nature of the interactions between the smallest structural particles of the material. The gradient coefficients, i.e. the factors multiplying the strain gradient terms, characterize the influence (or the importance) of these terms on the solution of the problem [10,11]. Noteworthy, in recent years, fractional non-local elastic models have been introduced which, by means of fractional calculus (i.e. changing the order of integro-differentiation), are able to span from integral to gradient non-local models [14][15][16].
Non-local fracture criteria date back to Novozhilov [17]. They are often referred to as theory of critical distances [18], and employed to assess the strength of notched components. The basic assumption of these approaches is that failure takes place when either a stress or an energy-related quantity at a finite distance from the notch tip reaches a critical value. Such a distance is the intrinsic length of the model, depending only on the material. However, these approaches fail in predicting the failure stress of a cracked or notched structure having size comparable to the critical distance. To overcome these incongruences, the coupled criterion of Finite Fracture Mechanics (FFM) was then introduced [19,20]. This model relies on the assumption of a discrete crack propagation under the simultaneous fulfilment of two conditions: a stress requirement and the discrete energy balance. Thus, the crack advance, i.e. the length characterizing the model, becomes a structural parameter, depending on both material and geometry. Being physically sound and efficient, FFM has been successfully applied to assess the strength of different materials, geometries and loading conditions [21]. Moreover, it has recently been proven that FFM provides close predictions to the well-consolidated Cohesive Zone Model, once the stress requirement and the cohesive law are properly matched [22][23][24].
In the present contribution, the attention is focused on the borehole problem (Fig. 1a), used as a test bench to compare GE and FFM approaches. The problem is chosen for its relevance in engineering applications, mainly geotechnical considering the inner pressure, but also mechanical if remote stresses are taken into account. Moreover, since the geometry is crack-less, LEFM is not efficient. On the contrary, GE and FFM are able to provide the critical value of the breakout pressure along with its dependence on the hole radius, i.e. the size effect. The predictions by the two approaches will be compared and checked against experimental data providing the breakout pressure for two different rock materials and different hole sizes [25].
The case of a circular hole under remote stresses has been already addressed both for what concerns GE [26] and FFM [27]. The originality of the present analysis lies in the solution to the inner pressure loading, a novelty for both approaches, as well as in the comparison between them, since, at the authors' best knowledge, it is the first time that Gradient Elasticity and Finite Fracture Mechanics are directly compared.

Gradient Elasticity
A special form of gradient elasticity is used in the present work, incorporating the Laplacian of the hydrostatic part of the strain tensor in Hooke's law. The appropriate constitutive equation reads [26]: where k; l are Lamé constants, d ij is the unit tensor (or identity tensor), e ij is the strain tensor, r ij is the stress tensor in indicial notation, and r 2 is the Laplacian operator. The gradient coefficient c multiplies the Laplacian of the hydrostatic part of the strain tensor; its value is related to a characteristic scale of the material, and its physical meaning can be attributed to the distance over which non-locality acts, smoothing high variations of the elastic stress field. Before applying the gradient model to the borehole problem, it is worth noting that GE yields either a set of fourth-order equations with the displacements as the sole unknowns or a set of coupled second-order equations with the displacements and the components of another variable (stresses, micro-deformations) as unknowns. In the former case (which is used in the present paper), the higher order of the governing equations means that higher-order boundary conditions must be formulated. In the latter case, additional boundary conditions accompany the auxiliary differential equations. Many studies on the issue of boundary conditions in gradient elasticity adopt second-order displacement derivatives. The reasons for this choice are based on energy arguments, removal of singularities, suppressing undesired boundary layer effects, or simply numerical convenience. For further information on this topic, see [13] and related references.
Let us now analyse the borehole problem, as depicted in Fig. 1a, i.e. a circular hole with radius R, under internal pressure p and remote biaxial (isotropic) tension r. Although in geotechnical applications the remote stress is negative, here we consider mainly traction since we focus on tensile failure. The consideration of a far-field tension rather than compression has been motivated by the case of drilling in an area of tensile remote tectonic stresses. Such areas are described in the work of Bott and Kusznir [28]. However, within the model assumptions, results are valid also for negative r; of course, higher inner pressures are required for cracking at the hole edge.
For the axisymmetric plane strain configuration at hand, shearing stresses and hoop displacement are null everywhere; the only displacement component different from zero is the radial one, u. Thus, in polar coordinates (r, h), the (classical) kinematic equations relating the radial and hoop strain to the radial displacement are: where we used the engineering notation, as in the following. Derivatives with respect to the angular coordinate h vanish, so that Eq. (1) yields: Fig. 1 The borehole problem: geometry, loadings and reference systems. Before (a) and after crack onset (b) Upon substitution of Eqs. (2) into (3) and (4), we get: The only meaningful equilibrium equation is the one in the radial direction, i.e. (neglecting the body forces): Upon substitution of Eqs. (5) and (6) into Eq. (7), a fourth-order ordinary differential equation for u is obtained [26]. Its general solution is derived as where A; B; C; D are integration constants, I 1 and K 1 are the first-order modified Bessel functions of first and second kind, respectively, and the characteristic (material) length l g has been introduced. It is a related to the gradient coefficient by: where m is the Poisson ratio. In order to determine the integration constants in Eq. (8), four boundary conditions are required: two classical boundary conditions, i.e. r r j r¼R ¼ Àp and r r j r!1 ¼ r, and two extra ones, commonly used with the gradient model, i.e. d 2 u=dr 2 r¼R ¼ 0 and d 2 u=dr 2 r!1 ¼ 0. Hence, the integration constants are as follows: where the dimensionless radius q g ¼ R=l g is introduced and K 0 is the zero-order modified Bessel function of the second kind. Hence, by Eqs. (8) and (10), the (radial) displacement field reads: with r ¼ r=R. The stress field components are then attained upon substitution of Eq. (11) into Eqs. (5) and (6): where the terms in square brackets at the right hand sides of Eqs. (12) and (13) represent the classical linear elastic solution, and the remaining parts the gradient ''corrections''. In Fig. 2 the hoop stress component is plotted vs. the radial coordinate for different hole radii in case of internal pressure alone (r = 0). For sufficiently large radii (i.e. for R approximately larger than l g ), the maximum stress over the whole domain is provided by the hoop stress at the hole edge (i.e. at r ¼ R): In the present gradient approach, the critical conditions are supposed to be achieved simply when the maximum normal stress r max attains a threshold value, that is, the material tensile strength r c (Rankine failure criterion). Equation (14) will be used in Sect. 4, where the results by GE and FFM are compared and checked against experimental data according to different loading conditions.
It is important to mention that for small radii (i.e. for R \ \ l g ), the actual hole size would be too small for the material to even ''feel'' it, and the related GE predictions could be not reliable.
In order to overcome this drawback, the present GE approach can be generalized by modifying the boundary conditions to Eq. (8) through an additional parameter, as done by Chen et al. [26] for circular holes under biaxial loading, or as discussed by Benvenuti and Simone [29] to catch size effects in micro/nano rods (see also [30]). Despite the process is arbitrary, a brief extension to present geometries will be given in Appendix 1 following the model presented in [26]: a correlative coefficient a will be introduced in the extra boundary conditions linking the radial strain gradient to that provided by the classical elastic solution.

Finite Fracture Mechanics
Let us now provide the solution by means of FFM, inspired by some previous works [31][32][33].
The coupled criterion of FFM [20,34] is based on the assumption of a finite crack advance. An abrupt crack appearance of length l is assumed to take place whenever two conditions are simultaneously fulfilled. The first is a stress condition, requiring the average stress over the crack advance to be higher than the tensile strength; the second is the discrete energy balance, requiring the available strain energy to be higher than the energy necessary to create the new crack surfaces.
For the geometry at hand (Fig. 1a), because of axial symmetry, crack onset may occur at any point along the hole edge. However, energy release per unit crack surface is highest when (only) two crack occur at the extremes of one diameter (Fig. 1b, a being the crack length), as it can be argued and checked by stress intensity factor (SIF) handbooks [35]. Thus, referring to Fig. 1b, the two FFM conditions may be written as: where G is the strain energy release rate and G c is the fracture energy. The crack onset load is the minimum one for which both inequalities (15) are fulfilled. It is worth mentioning that the first condition in Eq. (15) could be also expressed by imposing a condition straight to the stress field, and not to its average function. This leads to the FFM approach proposed by Leguillon [19]. The two FFM models yield different results, Eq. (15) always providing the most conservative predictions. By means of Irwin's relationship, the discrete energy balance can be expressed in terms of the SIF K I and its critical value, the fracture toughness K Ic . Hence, taking into account that cracks open under mode I loading, the latter inequality in Eq. (15) can be cast as: In the present FFM approach the material is assumed simply to be linear elastic. Thus, the stress field to be inserted in the stress condition is: coinciding with the first term of the gradient solution in Eq. (13). For the energy condition we need the SIF: where a ¼ a=R and The shape functions given by Eqs. (19) and (20) were evaluated through a boundary collocation method [35]: their accuracy is estimated to be better than 1%.
Let us now introduce Irwin's length l m : It is a material parameter, many studies proving its connection with the internal microstructure [18]. Thus, for FFM, l m plays a role analogous to the one played by l g in the GE approach. Similarly, the dimensionless radius q m = R / l m can be introduced.
Upon substitution of Eqs. (17) and (18) into Eqs. (15) and (16), we get: where l ¼ l=R is the crack advance normalized with respect to the radius and Analytical solutions for the integrals in Eq. (23) are provided in Appendix 2.
For positive geometries, i.e. for monotonically increasing SIF functions (18), the minimum load satisfying the system (22) is achieved when both inequalities are strictly fulfilled. Hence, Eq. (22) reverts to a system of two equations with two unknowns: the critical crack advancement l c (or, more precisely, its dimensionless counterpart l c ¼ l c =R) and the remote failure stress r f (once p is given) or the (breakdown) failure pressure p f (once r is given). On the other hand, for non-positive geometries, the search for the minimum load can be trickier (see, e.g., [27,36]). Further details will be given in Sect. 4.2 when considering the case of internal pressure alone. In any case, the FFM dimensionless failure load depends only on the parameter q m , exactly as in the GE model depends only on q g .

Discussion of results
GE and FFM models developed in the previous section will be now implemented and compared by considering separately the following loading cases: (i) remote tensile stress r (p = 0); (ii) internal pressure p (r = 0); (iii) remote tensile stress and internal pressure r ? p. Then, for case (ii), we will provide a comparison with experimental data from the literature.
As stated before, both approaches are based on a length, l g for GE and l c for FFM. Whereas the former results a material property, being involved in the nonlocal constitutive relationship (through Eqs. (1) and (9)), the latter is a function of the hole radius R, and the loading conditions, r and/or p, too. In other words, the FFM crack advance can be expressed by the relationship l c ¼ f ðR; r=pÞl m where l m is the material length provided by Eq. (21). Note that similar arguments hold also for the length scale parameter related to the Phase Field approach, whose connection with l c was investigated in [37].
This being clarified, in order to compare GE and FFM results, a proportionality coefficient between the internal lengths l g and l m is thus introduced: 4.1 Uniform biaxial tension (p = 0) Applying Rankine criterion and setting p = 0 into Eq. (14), the dimensionless failure remote stress according to GE is: For what concerns FFM, the failure remote stress is provided by Eq. (22) setting p = 0. Since the geometry is positive, the inequalities in Eq. (22) become equations. Simple analytical manipulations yield an equation in a unique unknown, the dimensionless crack advance: whose solution l ¼ l c must be replaced into, e.g., the first equation of the system (22) to get the failure remote stress: The dimensionless remote stresses at failure according to GE and FFM are plotted in Fig. 3: a good agreement is achieved by setting b = 0.24. For very large sizes, the strength is reduced by a factor 2, which coincides with the classical concentration factor. Decreasing the hole size, the strength increases in an almost coincident way according to the two models. Only for very low hole sizes the GE and FFM predictions start to differ from each other. Indeed, FFM provides coherently that the failure stress r f tends to tensile strength r c as the radius R vanishes. On the contrary, according to GE, the stress concentration factor becomes unrealistically lower than unity for very small radii, tending finally to 0 as R ! 0, i.e. an infinite strength is predicted. Therefore, the GE model is no longer valid below a cutoff size. This threshold size is approximately equal to GE internal length. In other words, for GE to work, it is required that, approximately, R [ l g (q g [ 1). However, this is not a big concern, since most of the experimental data are expected to fulfill this condition.
Finally, although the parameter b has been fitted to match the solutions by the two models, it is worth emphasizing the almost perfect agreement of GE and FFM predictions for q m [ 0.5 (i.e. q g [ 2), as clearly shown in Fig. 3. This nice feature is rather surprising given the completely different origin of the two nonlocal models, which, therefore, corroborate each other.

Pressurized hole (r = 0)
In this case, according to Eq. (14) and setting r = 0, the dimensionless failure pressure for GE is: for any q g . For what concerns FFM, the geometry is locally positive and globally negative, since the SIF (18) is first increasing (starting from zero) and then tending to zero as R tends to infinity. Thus, the failure pressure is again given by Eq. (22) (setting r = 0), but now two scenarios are possible. For radii R larger than 0.40 lm (q m [ 0.40), the minimum load fulfilling the two inequalities of system (22) is still achieved at the interception point between the curves obtained replacing the inequalities with the corresponding equations, as shown in Fig. 4a (where R = 2 lm). Analytical manipulations yield the following equation for the dimensionless crack advance: whose solution l ¼ l c must be replaced into the first equation of Eq. (22) to get the failure stress: On the other hand, for radii R below 0.40 lm (q m \ 0.40) , the minimum load fulfilling both the inequalities of the system (22) is not given by the interception point: it is given by the minimum of the discrete energy balance, as depicted in Fig. 4b (where R = 0.1 lm). This minimum is achieved when the crack advance l c equals 1.29 R. Replacing this value in the second equation of system (22) (along with r = 0) and evaluating I pp ðl ¼ 1:29Þ by the corresponding expression in Appendix 2 (Eq. (40)), we get: which can be rewritten equivalently as FFM and GE predictions (by assuming b = 0.24 in Eq. (24), as before) are reported in Fig. 5: the matching is again excellent except for vanishing radii. It should be noted that the ranges where GE and FFM nearly coincide are generally those of practical interest, as it will be shown in Sect. 4.4.
The critical FFM advance for the two cases analyzed before (p = 0 and r = 0) is plotted in Fig. 6. As can be seen, l c results a structural parameter depending on the hole radius R, besides on the material length l m . For both loading cases, l c /l m tends to 2/ (1.12 2 )p as the size increases. In this case, fracture is stress-governed and the crack advance is provided by the energy condition. The same happens at small sizes for p = 0, l c /l m tending to 2/p. On the other hand, for

Remote stress and internal pressure
In this case, according to Eq. (14), GE states that critical conditions are attained whenever the following equation is satisfied: for any q g . Once the dimensionless ratio q g is fixed, Eq. (32) defines the safety domain in the (p, r) plane, see Fig. 7. Looking at Eq. (32), it is clear that the domain has a triangular shape, the upper border simply being the straight line connecting the two points defining the critical conditions analysed in the previous subsections.
For what concerns FFM, the general analysis can be quite complex since, for small radii and large pressures, the SIF is first increasing, then decreasing and finally increasing again (vs. a). However, restricting the analysis to sufficiently large holes (q m [ 0.40) and positive remote stresses, the geometry remains positive and, hence, the failure load is simply achieved at the interception point, as in Fig. 4a. Let us assume that the remote tensile stress is c times the inner pressure, i.e. r = c p. Then Eq. (22) can be managed to yield the following equation in the finite crack increment: whose solution l ¼ l c must be replaced into the first equation of the system (22) to get the failure stress: Fixing the hole size (i.e. fixing q m [ 0.40) and solving Eq. (33) for different c values, we get the safety domain according to FFM.
Taking as before b = 0.24, in Fig. 7 we plotted the safety domain according to GE and FFM for q m equal to 0.5, 1, 2, ?, the last value coinciding with classical linear elasticity along with Rankine criterion. Differently from GE, the border of the safety domain is (weakly) non-linear according to FFM, the nonlinearity increasing for decreasing q m . This is due to the discrete energy balance, which is quadratic in the load. Nevertheless, the two approaches provide close predictions also in this more general case. In order to verify the soundness of the GE and FFM approaches, we provide a comparison with experimental data. The comparison is about the case of inner pressure alone, i.e. the case we investigated in subsection 4.2.
Cuisiat and Haimson [25] investigated the hole size effect on hydraulic fracturing breakdown pressure under zero far-field stresses by testing two rock materials: Lac du Bonnet granite and Indiana limestone. The material strength r c was measured experimentally, resulting in 8.1 MPa and 9.6 MPa, respectively. The authors applied the point stress criterion to the data, evaluating the stress field at a finite critical distance l c from the hole edge: l c was fitted to 0.30 mm for granite, and to 0.24 mm for limestone. Taking into account that, for the point method, l c ¼ l m =2p [18], we can get l m % 1:87 mm and l m % 1:48 mm for granite and limestone, respectively. FFM results are plotted in Fig. 8, together with GE results, implementing l g % 0:449 mm and l g % 0:355 mm via Eq. (24), b = 0.24.
As it can be seen in Fig. 8, the agreement between the two approaches and the experimental data is more than satisfactory. It should be added that the fitted results for the internal length l m would lead to the following estimations on the fracture toughness: K Ic ¼ 0:35 MPa ffiffiffiffi m p for granite, and K Ic ¼ 0:37 MPa ffiffiffiffi m p for limestone. Indeed, they turn out sensibly lower (2-3 times) than those generally evaluated experimentally, as it often happens when dealing with rock materials [39]. This aspect would require further investigations.

Conclusions
Gradient Elasticity and Finite Fracture Mechanics were applied to the borehole problem. Two different loading conditions were considered and investigated: uniform biaxial loading and internal pressure. This latter case represents a novelty for both approaches. The following conclusions can be drawn: • the two models have a completely different origin.
GE nature is non-local in the constitutive law and local in the failure criterion; it considers the (sizedependent) stress concentration factor as the governing failure parameter. FFM nature is local in the constitutive law, but non-local in the failure criterion: it considers the fulfilment of two (stress and energy) conditions for fracture onset. It is more difficult to obtain the stress field according to GE, but, once the stress field is known, GE is more straightforward, since it requires the implementation of just one single equation. • both models are based on a characteristic (internal) length, l g for GE and l m for FFM. They take into account the effect of the material microstructure at the macroscale. While l g is generally fitted on the experimental data, l m is achieved a priori once K Ic and r c are known. • by considering the following proportionality law l g & 0.24 l m , it has been shown that GE and FFM predictions are in excellent agreement over the range of practical engineering interest. This a key point, the independence of the two approaches strengthening each other. Moreover, it follows that l g can be linked straightforwardly to the brittleness of the material in future GE applications: and E is the Young's modulus.
• for vanishing hole sizes, FFM catches coherently the transition to a plain geometry. On the other hand, GE presents some drawbacks as the internal length becomes comparable to the hole radius so that a threshold size has to be considered below which the model can not be applied. • the excellent agreement with experimental data related to two different rock materials (Lac du Bonnet granite and Indiana limestone) proves the soundness of the present formulation.
Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement.

Declarations
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix 1
Let us modify the two extra-boundary conditions of Eq. where is the radial displacement component according to (classical) linear elasticity, and the correlative coefficient a ranges between 0 and 1, i.e. 0 B a B 1 [26]. Through some analytical manipulations, similar to those presented in Sect. 2, the maximum hoop stress at the hole edge expressed by Eq. (14) reverts to: r max ¼ r h j r¼R ¼ r þ ðr þ pÞ 1 À 2 ð1 À aÞK 0 ðq g Þ q g K 1 ðq g Þ þ K 0 ðq g Þ " # : The parameter a has to be chosen to fulfil the asymptotic condition for vanishing hole sizes, R ! 0. It is immediate to check that setting a = 0.5 satisfies this requirement for both geometries under investigation: accordingly, the stress concentration factor Biaxial tension (a) and pressurized hole (b), failure stress / critical pressure vs. hole radius: comparison between FFM and GE results by setting a = 0.5 in Eq. (36) reduces to 1 (r max ! r) in case of biaxial loading (p = 0), whereas it diverges (r max ! 0) in case of a pressurized hole (r = 0). The comparison with FFM results is shown in Fig. 9: both approaches catch theoretically the size effects even for small holes, although the differences are not negligible (nearly 15% for R/l m = 0.1, if p = 0). Finally, according to the approach described above, Eq. (32) reduces to: r r c 2 1 À ð1 À aÞK 0 ðq g Þ q g K 1 ðq g Þ þ K 0 ðq g Þ