## Abstract

An assessment is made here of the role played by the micropolar continuum theory on the cracked Brazilian disc test used for determining rock fracture toughness. By analytically solving the corresponding mixed boundary-value problems and employing singular-perturbation arguments, we provide closed-form expressions for the energy release rate and the corresponding stress-intensity factors for both mode I and mode II loading. These theoretical results are augmented by a set of fracture toughness experiments on both sandstone and marble rocks. It is further shown that the morphology of the fracturing process in our centrally pre-cracked circular samples correlates very well with discrete element simulations.

## 1. Introduction

The discrete element method (DEM) has become an increasingly popular choice for modelling discontinuous systems because of its success in describing complex deformation patterns without the need of sophisticated continuum constitutive laws. In the area of geomechanics, this strategy has proved to be particularly versatile in modelling particulate solids (soils, rocks, loose sands, etc.) and their associated failure mechanisms, and in the oil and gas industry, the DEM has been employed in studies of hydraulic fracturing, drilling tool/rock interaction and sand production in oil reservoirs (e.g. [1–3]).

Broadly speaking, the method is based on describing deformable bodies as large and dense collections of rigid bodies (sometimes referred to as *particles* or *grains*) that are in contact with each other. Typically, the particles have regular geometries (spheres in three dimension or discs in two dimension), their diameter being much smaller than the characteristic length of the continuum body that their collective behaviour aims to describe (figure 1). A fictitious linear-spring/dashpot/slider is employed to account for the interaction between any two particles in contact with each other. The stiffness of the elastic springs and the other ‘material’ constants for the sliders and the dashpots are not directly accessible from measurements, but are related to the bulk properties of the physical solid (Young's modulus, Poisson's ratio, etc.) and need to be calibrated through simulations of typical rock/soil mechanics tests such as, for example, bi-axial compression and indirect tension test [4–6]. Numerical simulations by DEM-type methods are dynamic in nature, the motion of the particles being governed by the usual linear and angular equilibrium equations of Newtonian mechanics, while the time integration is carried out with an explicit algorithm using a sufficiently small time step. Nonetheless, static behaviour can be easily described by introducing sufficiently large numerical damping (e.g. [7]).

The existence of vastly different multiple spatial scales is implicit in any DE model. At one end of the spectrum, there is the microscopic scale, i.e. the scale on which the particles are regarded as independent, rigid units interacting with each other; at the other end of the range, the continuum scale, mechanical response is governed by the behaviour of the entire assembly of particles in which the microscopic effects are smoothed out. An intermediate spatial scale that bridges the gap between the micro- and macroscopic scales already mentioned also exists, and it is associated with a suitably defined representative volume element (RVE) in a relevant homogenization process.

Owing to the inherent microstructure present in synthetic (i.e. virtual) rock materials modelled within the DE framework, starting fairly early on in the development of the subject, there have been numerous efforts (e.g. [8,9] and the references therein) attempting to establish the formal connection between DE solids and the general micropolar elastic theory [10–13]. This connection takes the form of continuum constitutive laws in which various material coefficients are derived explicitly in terms of grain and contact properties, thus allowing the formulation of boundary-value problems that can be solved more efficiently by standard numerical methods. Unfortunately, most of these studies ignore checking the validity of the results obtained by comparison with experiments or DEM simulations. Some notable exceptions are the recent works of Alassi & Holt [14] and Gerolymatou [15], who used Walton's [16] earlier homogenization scheme for random packings of elastic spheres. Homogenization in the context of Cosserat elasto-plasticity has also been considered by many authors (e.g. [17,18]).

One of the unique features of the DEM consists in the relative ease with which it can simulate crack initiation and propagation in the context of bonded assemblies of rigid spheres (three dimensions) or discs (two dimensions). This is made possible by identifying fractures with the progressive breaking of bonds between the rigid units. Although when taken at face value, the ensuing discontinuous crack propagation seems quite remote from the brittle fracture described by continuum theories, work done in the past decade has demonstrated that it is possible to reconcile the apparently two different points of view. Potyondy & Cundall [5] established a formal relationship between mode I fracture toughness and the microproperties of cubically packed bonded assemblies. Later, Moon *et al.* [19] developed a more general scheme to measure *K*_{IC} under random packing of non-uniform-size two-dimensional particles, by using both an energy balance approach and the collocation method based on the generalized Westergaard approach. Tavarez & Plesha [20] proposed an elasto-plastic modification of the usual contact bonds and performed simple fracture toughness numerical experiments demonstrating the possibility of achieving convergence with refinement of the DEM discretization. Applications of softening contact models and their implications to fracture and delamination in fibre-reinforced composites were studied by [21,22], and Hedjadzi *et al.* [23] made comparisons between finite element modelling and the DEM in regard to unstable crack propagation and branching in a vitreous dense biopolymer material. They found that the latter route led to closer agreement with both an existing analytical solution and the morphology of the cracks observed in the experiments.

Before giving details about the organization of the paper, we sketch briefly some background that will help set the stage. The effect of couple stresses on the *stress concentration factor* (SCF) for a uniformly stretched rectangular plate containing a circular hole or an elastic space with a spherical cavity was studied by Mindlin [24], and by Muki & Sternberg [25] for a much wider class of scenarios. By way of example, in the case of a plate with a hole, Mindlin found that the classic elastic value SCF^{elas}=3 (independent of the size of the hole) was changed to an expression depending on Poisson's ratio and the radius of the hole, as well as a new material length-scale reflecting the existence of a microstructure. Furthermore, the upper bound for the new value, SCF^{mpol}, was found to be equal to 2.0, hence suggesting that couple stresses reduce the severity of stress concentrations around geometrical discontinuities. The more difficult task of finding the local, near-tip stress field (including the *stress-intensity factors*, *K*^{mpol}_{I} and *K*^{mpol}_{II}) for a crack in a material with couple stresses was also attacked with numerical methods by Sternberg & Muki [26]. Atkinson & Leppington [27] revisited that work and proposed a general analytic framework based on Wiener–Hopf techniques and matched asymptotics that succeeded in producing closed-form expressions for the stress-intensity factors in micropolar solids (including as a particular case the indeterminate couple-stress theory). In more recent times, their solution strategy has been adapted to the simpler case of mode III loading of a semi-infinite crack in a couple-stress material by Radi [28] and Zhang *et al.* [29].

With these comments in mind, we start our investigation in §2 with a quick review of some fracture toughness experiments that we carried out on two different rocks by using centrally cracked discs. The determination of the traditional *K*_{IC} and *K*_{IIC} is based on the formulae developed in [30], although our interest here is primarily on the fracture patterns developing in the tested rock samples. The equations of micropolar theory, which contain as a particular case those relevant for the couple-stress theory of Mindlin, are recapitulated in §3, in a form that allows for a natural interpretation of the numerical results reported later in the paper. To emphasize the close similarity between experiments and DE-type rock models, we include a representative sample of numerical simulations in §4; a more in-depth exploration of the link between these simulations and the experiments discussed in the next section will appear elsewhere [31]. Sections 5 and 6 are largely based on [27], and they include a number of new results for a semi-infinite crack in a micropolar material subjected to mode I and mode II deformations. The particular geometry of our fracture toughness experiments is accounted for in §7; based on the far-reaching analytical strategy of Atkinson & Leppington [27], we establish the leading-order effect that the couple stresses introduce in the finite geometry of the centrally cracked disc. The paper concludes with a discussion of our main findings, together with suggestions for further study.

## 2. Experimental work

The cracked Brazilian disc test (CBDT) is used frequently in rock mechanics to investigate mixed-mode fracturing (figure 2), and we choose to use this particular configuration as a starting point for our subsequent discussions. The rock specimens used in CBDT tests are circular, with a through-crack along one of their diameters, and a fracture is initiated by applying opposite compressive loads. The inclination angle *α* between the original position of the crack and the applied load determines the mode mixity of the resulting fracture process. Mode I requires the initial crack to be vertical, as seen in figure 2*a*, whereas mode II is obtained for an inclination angle that depends on the crack length and the radius of the specimen, as explained shortly.

Our experiments were carried out on circular samples of Lazonby sandstone and Carrara marble with nominal diameter and thickness of 52 mm and, respectively, 15 mm. Lazonby is a medium-grained, low-porosity sandstone, which has a consistent salmon-red colour enhanced by the sparkle of quartz grains, originating in Cumbria, UK. Both types of rocks have an unconfined compressive strength of approximately 110 MPa.

Rectangular notches of 2 mm width and 20 mm length were cut in the central part of each sample by machining them with a water jet. These notches were later sharpened down to 0.25-mm width, with the help of a 0.25 mm thick hard steel blade and an abrasive paste consisting of a mixture of water and silicon carbide (SiC) powder. Once the mixture was homogenized, it was applied to the notch tip using a syringe, and then the rock was sawed. The length of the narrow notch was roughly 2.5 mm, with the total length of the final central crack being about 25 mm. Figure 3 shows the notch produced by the water-jet cutting (left end) as well as the final sharp notch obtained by sawing (right end).

The rock samples were used in a compression setup in which the compressive load was transmitted through two parallel 10.0 mm thick tungsten platens. The machine employed in our experiments was an Instron 4505 with a maximum load of ±100.0 *kN*. The crosshead speed was set at 100.0 mm ⋅ min^{−1} in all of the tests, and a small piece of cardboard was placed between each circular sample and the tungsten platens to accommodate the rock discs. Data recording was made using a PC with an Advantech PCL-818h acquisition card that is capable of registering simultaneously 16 channels with a resolution of 12 bits and a speed of 10^{6} samples per second.

The calculation of the stress-intensity factors (and thus the mode mixity as a function of the angle required to obtain pure mode II conditions) is based on the formulae provided by Atkinson *et al.* [30]. Following their nomenclature, the usual stress-intensity factors *K*_{I} and *K*_{II} can be expressed in the form
2.1
where *F*_{I} and *F*_{II} are non-dimensional shape coefficients depending on the quantities indicated in equation (2.1). Atkinson *et al.* [30] provided expressions for these functions as complicated infinite series, although it was confirmed that in most situations of practical interest a five-term truncation of those series would suffice. It was also shown by the same authors that for 0<(*a*/*R*)<0.3 one can treat the BCDT problem as that of a small crack in an infinite medium, and this led to
2.2a
and
2.2b
As an example, we include the variation of the shape coefficients as a function of *α*, for *a*/*R*=0.5 (figure 4*a*) and *a*/*R*=0.2 (figure 4*b*), using both the five-term truncation from [30] (continuous line) and the simplified equations (2.2) (dashed lines).

By using equation (2.2a) it becomes immediately obvious why the CBDT can be used to determine fracture toughness in mode II. Indeed, setting *F*_{I}=0, one ends up with an equation for the inclination angle *α*≡*α*_{II}, which depends on the ratio *a*/*R*; the same argument can be applied for the more accurate expression of *F*_{I} given in [30], and we used those results to calculate the critical angle *α*_{II} at which the crack is under pure mode II loading. It was found that the corresponding values varied between 22.5° and 24° for the geometry used in our experiments.

Figure 5 shows the load-displacement dependence for all mode I (*a*) and mode II (*b*) tests carried out on marble samples. The nonlinear increase in load observed during the initial stage in all of the experiments is due to the accommodation of the samples on the aforementioned cardboard material ‘sandwiched’ between them and the tungsten platens. The behaviour displayed by the circular rock discs that failed under the same fracture mode is similar in both sets of results, and it is evident that the maximum load levels reached are comparable (largely due to the isotropy of the marble samples used). In the experiments, the difference observed between the rock response under mode I and mode II loading was almost imperceptible.

Two consecutive frames recorded in one of the mode I experiments for marble are included in figures 6 and 7. In figure 6*a*, the sample is still intact, although close visual inspection revealed the existence of two one-dimensional regions of high stress concentration emanating from both ends of the notch and extending vertically. The snapshot that appears in figure 6*b*, taken 1/10 of a second later, is radically different: now, one half of the sample has violently moved away, while the other remains jammed in the testing machine. The response associated with this behaviour is identified as the thick red curve in figure 5*a*. The sandstone counterpart of the information seen in figure 5 is reported in figure 8, similar observations being applicable here as well. However, in this case the failure load was clearly different for mode I and mode II tests, and we also note that once the load reached its peak value, it did not drop abruptly as in the tests on marble. This particular feature allowed us to take snapshots of the progressive failure of the rock discs. A representative sample of such results is included in figure 9*a*,*b*, it is easy to see how the two primary cracks have propagated from either notch tip, extending further to the edge of the circular configuration. For mode II (figure 9*b*), we have also consistently noted the development of secondary cracks that started at the outer rim of the sample and propagated inwards, roughly along the symmetry axis of the crack, before eventually reaching both tips of the original central crack.

By using equation (2.1) in conjunction with the numerical data from figures 5 and 8, we calculated the critical stress-intensity factors *K*_{IC} and *K*_{IIC} for both types of rocks. As each test was repeated under identical conditions, here we report mean values together with the corresponding standard deviation. For Lazonby sandstone it was found that *K*_{IC}=0.692 (s.d.=0.69%) and *K*_{IIC}=0.625 (s.d.=0.50%), whereas for Carrara marble, *K*_{IC}=1.257 (s.d.=0.75%) and *K*_{IIC}=1.529 (s.d.=0.55%); the unit for the stress-intensity factors is MPa ⋅ m^{1/2}.

## 3. The equations of micropolar theory

Next, we review the mathematical framework used for the study of the linear micropolar continuum discussed later in the paper. For the sake of definiteness, we shall tacitly consider a right-handed orthonormal system of coordinates whose axes are described by the unit vectors *e*_{i} (*i*=1, 2, 3). However, the outline of the theory included below is given using a coordinate-free approach.

The equilibrium equations for a linear micropolar continuum assume the form (e.g. [10] or [12])
3.1
where ** t** and

**are the asymmetric force-stress and couple-stress second-rank tensors characterized by the general constitutive equations 3.2a and 3.2b with**

*μ***and**

*ω***being the macro- and micro-axial rotation vectors, respectively, and the symbol ⊗ stands for the usual tensor product. In (3.1),**

*Φ*

**t**_{×}is the vector of the tensor

*(defined to be minus twice the axial vector of its skew-symmetric part), while the macro-rotation in equation (3.2a) is related to the infinitesimal displacement vector*

**t***by . Other notations used above are defined as follows:*

**u****represents the classic infinitesimal strain tensor , |**

*ε***| is its first principal invariant,**

*ε**I*represents the three-dimensional identity tensor, ∇⊗

**≡**

*Φ*

*Φ*_{j,i}

*e*

_{i}⊗

*e*

_{j}, and

**⊗∇≡(∇⊗**

*Φ***)**

*Φ*^{T}, with the superscript T indicating the transposition operation for second-rank tensors. The micro-rotation is an independent kinematical quantity that has no counterpart in the classic continuum theories. The constants λ and

*μ*are the usual Lame constants, and

*γ*,

*α*,

*η*and ℓ* are new parameters that characterize the constitutive behaviour of the micropolar solid.

An important particular case is that of (indeterminate) couple-stress theory that corresponds to the situation in which ** ω**≡

**, and which is obtained by formally taking the limit . Note that now ∇⋅**

*Φ***=∇⋅**

*Φ***≡**

*ω***0**, as the divergence of the vector of any tensor vanishes identically. Also, the last term in (3.2a) is no longer present, and that equation determines only the symmetric part of the force-stress tensor; furthermore,

*μ*in (3.2b) reduces to the deviator of the couple stress.

It is sometimes convenient to discuss the microstructural properties of polar media by introducing certain combinations of the material moduli that feature in (3.2). These are ℓ_{b}, ℓ_{t} and *N*, which, in the notation of this paper, are defined according to
3.3
where the last parameter (*N*) is non-dimensional and is known as the *coupling number*, whereas the other two have dimensions of length; ℓ_{b} and ℓ_{t} are characteristic lengths in bending and torsion, respectively. The coupling number reflects the degree of polarity present in the material [32]; the upper bound *N*=1 (obtained when *μ*/*γ*→0) corresponds to the indeterminate couple-stress theory mentioned earlier. Although ℓ_{b} and ℓ_{t} feature in many theoretical/numerical treatments available in the literature, experimental values of these parameters for solids are scarce. In Lakes [33], examples of such values are found for two foams with nearly spherical voids; ℓ_{b}=0.032 mm and ℓ_{t}=0.065 mm for a syntactic foam consisting of hollow glass microbubbles embedded in an epoxy matrix, with the mean diameter of the voids being 0.125 mm and the volume fraction equal to 0.468. The other material tested was a high-density polyurethane closed-cell foam, for which ℓ_{b}=0.327 mm and ℓ_{t}=0.62 mm; in this latter case, the mean diameter of the voids was 0.1 mm and the volume fraction 0.69.

The notations used in the literature on micropolar elasticity are rather confusing because often the same letter is used by different authors to denote different things. In table 1, we summarize the correspondence between the notations used in [12,27,34]. Above, we followed Bigoni & Drugan's [34] notation for the micropolar moduli, except for the parameter ℓ*, which in their notation is labelled ℓ. It is worth emphasizing that ℓ_{b} and ℓ_{t} will be the same whether one considers couple-stress or micropolar solids.

We are concerned exclusively with a plane-strain situation that involves deformations taking place in the 12-plane. In this particular case, the equations (3.1) and (3.2) admit a number of simplifications that will not be reiterated here ([12], pp. 173–178, for instance), and the solution of the linear micropolar elasticity is readily available with the help of Mindlin's stress potentials [24], *ϕ*≡*ϕ*(*x*_{1},*x*_{2}), *ψ*≡*ψ*(*x*_{1},*x*_{2}), defined according to
3.4
where **I**_{2} is the two-dimensional identity tensor in the 12-plane, * N*(

*ϕ*):=(

**∇**

^{2}

*ϕ*)

**I**_{2}−∇⊗∇

*ϕ*, and a similar relation holds for

*(*

**N***ψ*). Because of the plane-strain assumption the only non-zero components of the couple-stress tensor are

*μ*

_{ρ3}and

*μ*

_{3ρ}(

*ρ*=1,2), and in (3.4), we have tacitly assumed that

*μ*≡

*μ*

_{13}

*e*

_{1}+

*μ*

_{23}

*e*

_{2}as

*μ*

_{31}and

*μ*

_{32}do not feature in the equilibrium equations and can in fact be obtained directly from

*μ*

_{13}and

*μ*

_{23}, respectively (e.g. [12]).

The compatibility conditions for the strain measures (e.g. [10,12]) demand that the stress functions must satisfy the following set of partial differential equations
3.5a
and
3.5b
where
and (in our current notation)
Note that, in the couple-stress limit (i.e. ), we have ℓ→ℓ*. In [27], the key parameters for the micropolar theory were
3.6
and it is easy to see from table 1 that the limit of as is just (1−*ν*), which explains why in that work the micropolar solutions were obtained from the corresponding couple-stress results by performing the substitution
3.7
Relating the micropolar constitutive constants that appear in (3.2) to the micromechanical properties of discrete element solids can be traced to the contributions of Chang *et al.* [8] and Mühlhaus & Oka [35]. In general, the continuum theories that emerge from such analyses are of the higher order strain-gradient type, not necessarily directly related to the models discussed in this paper. However, for the sake of completeness, we illustrate below the correspondence between our micropolar parameters and the microproperties of a particular bonding strategy between the discrete element spheres in three dimension, in which each such sphere is elastically connected to its immediate neighbours. The links are idealized as two sets of three mutually perpendicular linear springs that act at the contact position between any two spheres that touch each other. The stiffnesses of these springs will be denoted by *K*, and their nature will be distinguished by the relevant subscripts. One set consists of a normal spring (*K*_{n}) and two identical shear (tangential) springs (*K*_{s}), and the other includes a torsional (*K*_{φs}) and two identical bending springs (*K*_{φn}). See [36,37] for further details.

The springs can be visualized by regarding the spheres as being bonded together through short linearly elastic isotropic cylinders of height *h* and radius *r*_{0}; such microcylinders can be seen in figure 1*a* (see also [38] for a similar concept implemented in the commercial software *PFC*3*D*). By following this idealization, the following expressions for the stiffnesses can be deduced
where *E*_{b} and *G*_{b} denote the Young modulus and the shear modulus, respectively.

By further assuming that the aforementioned collection of spheres is statistically homogeneous, and employing an averaging approach, it was claimed in [36,37] that the idealized particulate solid can be described by the equations of micropolar theory by explicitly identifying the micropolar moduli as follows:
3.8a
and
3.8b
where *f*:=*Zv*_{s}/(*πD*), with *Z* being the coordination number, *v*_{s} the volume fraction of the spheres and *D* their average diameter. As far as we are aware, no direct verification of these results exist with actual numerical simulations. Variants of () have been obtained by many others, but reviewing that body of work is beyond the scope of this paper.

An order-of-magnitude estimate of the micropolar moduli was presented in [36], in which the authors assumed that *E*_{b}∼2*G*_{b}, *b*∼*D*/2, *Zv*_{s}∼1.0, *G*_{b}∼10 GPa, *D*∼10^{−3} m and *h*/*D*∼0.1. With these assumptions, it was found that

## 4. DEM model of the CBDT

At this point, we shall digress briefly in order to illustrate the striking similarity between the behaviour of the centrally cracked discs considered in §2 and some DEM simulations. Although, from a strictly mathematical point of view, the cracks discussed in this paper are simply line singularities, the numerical experiments require avoiding the possible friction between the opposite crack faces. For this reason, the initial crack in the simulations looks very similar to those in the experiments (see also figure 2); the geometry and the loading conditions are actually identical.

Simulations were carried out with the two-dimensional version of the commercial code *PFC* [38], which provides several options in regard to the selection of the bonding between the particles. The results shown below correspond to the so-called *bonded-particle model* [5,7] that allows for the grains in the DE model to transmit both forces and moments to their immediate neighbours. A more detailed account of our modelling and simulation work will be reported elsewhere [31], so here we limit ourselves to the global ‘picture’ of the fracture process.

Figures 10 and 11 show two typical simulations of the CBDT using a hexagonal packing, each particle away from the boundaries of the sample being surrounded by six neighbours in the undeformed configuration. This particular choice was motivated by the recent results in [14]. In both figures, the left window shows the samples at some instant of time after the peak of the load versus cross-head displacement has been reached. Visual inspection of the mechanical failure experienced by the synthetic rock samples compares well with the experimental snapshots of the corresponding experiments reported earlier in this paper. A coarse-graining procedure (e.g. [39]) was used in the plots on the right to calculate the maximum principal stress in the discs, just before the bonds between the particles started to break; the dark shade of red in those plots represents sites of high tensile stress concentration, and it can be seen that their location correlates well with the subsequent failure processes of both rock specimens.

## 5. The tensile crack (mode I)

The analytical calculations presented in the remaining sections of the paper are largely based on the earlier work of the first author [27], so we shall avoid repeating similar theoretical arguments here; instead, only the final results will be stated. To facilitate the link with the work just cited, the same notation will be employed, although we do show in several places how a number of the formulae can be easily interpreted in light of the more natural notation of the current paper.

We start off by reviewing here the well-known classical elastic theory. In the polar coordinate system (*r*,*θ*) based at the tip of the crack, the local singular stresses have the form
5.1
and the corresponding displacement field is given by
5.2
where *E* is Young's modulus, *ν* denotes Poisson's ratio and *θ*=0 is directly ahead of the crack tip. The coefficient *K*_{I} (i.e. the stress-intensity factor in mode I) depends on the load and crack length, as well as the geometry of the specimen, but the local variation is universal.

For couple-stress theory, the near-tip singular crack field is given by (e.g. [27,40])
5.3a
5.3b
5.3c
5.3d
5.3e
5.3f
5.3g
where the constants *k*_{1} and *k*_{2} are proportional to the new stress-intensity factors that reflect the addition of couple stress (and are not to be confused with *K*_{I} or the bulk modulus *K*, etc.). In the above equations, one must employ the substitution (3.7) to obtain the corresponding micropolar results.

Local expansions for the singular fields at the crack tips have also been made by Diegele *et al.* [41]. It is a simple matter to get these expansions, but that is not our primary concern here. For completeness, we include below the well-known simplified formulae for the stresses directly ahead of the crack. To this end, reference will made to the usual right-handed local system of *x* and *y* coordinates attached to the tip of the crack (similar to the coordinates *X*_{1} and *X*_{2} in figure 2).

For the *elastic* case
5.4a
and
5.4b
where [[*φ*]] denotes the normal jump of the scalar field *φ* across the crack faces, formally defined as [[*φ*]]:=*φ*(*x*,*y*+0)−*φ*(*x*,*y*−0). Then the energy release rate as the crack advances is
5.5
In the *micropolar* case, by letting denote the mode I stress-intensity factor, we find
5.6a
and
5.6b
so the work done by *the field* *t*_{22} for virtual crack advance will be
5.7

Note, however, that this is not the only contribution to the energy release rate as the crack advances as the singular couple-stress field also does work. This is
5.8a
and
5.8b
where *k*_{3} is a factor that plays the role of a couple-stress intensity factor in the micropolar theory. These combine to give for the couple stress energy release rate , so that the total energy release rate is
5.9
The coefficients *K*_{I}, *K**_{I} and *k*_{3} are as yet unknown, but Atkinson & Leppington [27] have evaluated these for the problem of a semi-infinite crack with a special loading and also used singular perturbation theory for more difficult situations. The results are complicated, but one can write explicitly (in terms of complex integrals) those expressions for the micropolar case and compare it with the energy release rate for the corresponding elastic case, and hence find the effect of micropolar media on the flow of energy to the crack tip. Thus, if the resistance to fracture of the material remains the same as we increase the nature of the discrete medium, an estimate of this effect can, in principle, be made.

It is worth stressing that a consequence of these couple-stress (or micropolar) continuum theories is that there is an effect at the crack tip that changes the magnitude of the stress-intensity factor such that the limit as the couple-stress parameter (ℓ) tends to zero is *not* the same as that obtained from the purely elastic theory with ℓ=0 from the outset. However, the energy release rate given in equation (5.9) *does* tend to the elastic one recorded in (5.6). Thus,
5.10
In closing this section, we use equation (A.2) from appendix A to study the dependence of and on (ℓ*/*a*). To this end, it is important to observe that, in the notation of §3, the various parameters that enter into the expression on the right-hand side can be expressed as
5.11
Note also that Δ=1 corresponds to the couple-stress theory, so values of Δ close to unity will be indicative of the effects in that limiting case, whereas (e.g. Δ=5.0, etc.) are representative of the general micropolar theory. Weakly micropolar effects correspond to Δ≫1 and (ℓ*/*a*)≪1, which formed the object of the limited numerical comparisons reported in [27]. By fixing Δ, the only other free parameter left in the equation that we want to plot will be Poisson's ratio (*ν*), and we can then illustrate the aforementioned dependence for a range of values of *ν*. By letting be the micropolar quantity corresponding to *k*_{3} that appears in equation (5.8), we show its variation with (ℓ*/*a*) in figure 12, where the arrow indicates the direction in which *ν* increases. It is interesting to observe that decreases as (ℓ*/*a*) increases and, generally, has a small magnitude. Its largest values are associated with the couple-stress theory and, as expected, this effect drops down to zero as we approach the elastic limit (bottom part of the figure).

## 6. The shear crack (mode II)

We consider the semi-infinite crack loaded under shear by an internal shear stress , *x*_{1}<0 on *x*_{2}=0. From symmetry, we have the boundary conditions
and
This problem is simpler than the mode I case because there is now no work done by the couple stress as the crack advances incrementally. Using the Fourier transform, we can reduce the corresponding boundary-value problem to the functional equation ()
6.1
where
and
Equation (6.1) above can be routinely solved by the Wiener–Hopf technique. We omit the details (see Atkinson & Leppington [27] for the much more complicated mode I crack case).

The final results for the near crack-tip shear stress and displacement can be written
6.2
and
6.3
where
6.4
and we recall that *k*_{0}≡3−2*ν* and *β*_{1}≡1/(1−*ν*_{1}).

Note that in the above equation we replace (1−*ν*) as stated in equation (3.7). Also, from (3.6) only *c*_{1} and *b* involve the coefficient *γ*, which relates the couple stress to the gradient of the micropolar component *Φ*_{3}. Thus, (*c*_{1}/*a*) depends upon *γ*, but (*b*/*c*_{1}) does not.

It is interesting to compare the above results with an equivalent elastic result, i.e. one in which we just set *γ*=0, so *b*=0 and then the expressions in (6.2) and (6.3) become
6.5
and, respectively,
6.6
The above results also apply to the classical couple-stress problem in which, in place of *k*_{+}(i*c*_{1}/*a*), we have *k*_{+}(iℓ/*a*) and the underlined factor in (6.6) is replaced by 2(1−*ν*)/*μ*. In both cases, we see that the stress-intensity factor in (6.2) is not the same as that in (6.5) as (*c*_{1}/*a*) or (ℓ/*a*) tend to zero. However, the energy release rate for this crack (assuming it propagates straight ahead) is given by
6.7
where is the corresponding result for the elastic crack with local crack-tip stress field given by (6.5) and (6.6). Note that (also, when (ℓ/*a*)→0 in the couple-stress case). Let us mention briefly that due to the convention stated in equation (3.7), in the micropolar case
6.8
Also, note that, by using the correspondence explained in table 1, the limiting value of the *κ*_{0} defined in (6.8) as μ/γ→∞ turns out to be (3−2*ν*).

In figure 13, we provide a quantitative comparison between the standard elastic case and the micropolar theory. The strategy we follow is the same as that used in the plots shown in figure 12: a sample of three values for Δ, the quantity defined in equation (5.11), was chosen in conjunction with an additional four values for *ν*=0.1,0.2,0.3,0.4. Figure 13*a* shows the dependence of on (ℓ*/*a*). From the data included therein, it can be seen that (a) as expected, the couple-stress theory has a more pronounced effect on the flow of energy at the crack tip, and (b) in either case, increasing the Poisson's ratio reduces the effect of couple-stresses on the ratio . In figure 13*b*, the ratio between the *micropolar* mode II stress-intensity factor, *K*^{(mpol)}_{II}, and its elastic counterpart, *K*^{(elas)}_{II}, is plotted against the characteristic parameter (ℓ*/*a*). This ratio is obtained by formally using the near-field expressions stated in (6.2) and (6.5). Note that, rather unexpectedly, the micropolar stress-intensity factor is always bigger than the elastic version although, as the couple-stress effects get weaker (see bottom set of curves in figure 13*b*), the difference is less marked. Nevertheless, as figure 13*b* shows, the ratio of the energy release rates in the two cases tends to unity (irrespective of what the Poisson's ratio is).

## 7. Central crack in the disc

In this section, we take advantage of the results derived by Atkinson *et al.* [30] for the CBDT to account for the presence of the length-scale parameter ℓ, which here will be assumed to be very small (i.e. 0<ℓ≪1). With this assumption in place, one can then use singular perturbation theory to deal with the particular geometry of the CBDT. Only a leading-order analysis will be pursued here.

### (a) Mode I

For this problem the near-crack-tip field has the form given by equations (5.1) and (5.2) of §5, where *K*_{I} corresponds to the expression recorded in (2.1). According to our brief discussion of the micropolar and couple-stress equations in §3, it is clear that in both theories the relationship between the two-dimensional couple-stress vector ** μ** and the micro-rotation

**depends on a numerical factor proportional to ℓ**

*Φ*^{2}. It also transpires from a quick glance at table 1 that both

*b*and

*c*

_{1}defined in (3.6) are , whereas the ratio

*c*

_{1}/

*b*is independent of ℓ.

For most of the region of the disc, the limiting ‘outer’ solution is obtained by setting ℓ=0, and is simply the elastic solution (without couple stresses) obtained by Atkinson *et al.* [30]. However, this outer solution fails near the crack tips, where the effect of the couple stresses must be properly accounted for.

For the tensile crack, this has been treated by Atkinson & Leppington [27], and we briefly describe it here. First, an inner coordinate system is used near each edge, with origin based at the crack tip (*x*_{1}=*a*, *x*_{2}=0, say); see figure 2 for further details. One writes
7.1
so that with this scaling, no explicit dependence on ℓ appears in the equilibrium equations of the classic couple-stress theory or in the micropolar equations.

The local inner problem rescaled with respect to ℓ is that of a semi-infinite crack. Precise boundary conditions at infinity for the inner problems are determined by matching with the outer solution. Since both inner and outer approximations are to hold in common regions ℓ≪*r*≪*a*<1 near each edge, the two approximations must be asymptotically equivalent there.

After some analysis (see [27] for full details), it can be shown that the near-crack-tip stress and displacement has the form (5.6a) and (5.8a), with
7.2
and
7.3
where *d* and *c*_{0} are as defined in appendix A, and *K*_{I} is the elastic stress-intensity factor given by (2.1). These results are strictly true only in the limit 0<ℓ≪1.

### (b) Mode II

The solution strategy for this situation follows closely what was said in the last section. We begin with the elastic solution in the disc (i.e. the outer problem) from which we calculate an elastic *K*_{II}; again, this corresponds to our equation (2.1) that was originally obtained in [30].

As in the previous section, we also rescale the coordinates about the crack tip *x*_{1}=*a*, *x*_{2}=0 (figure 2), and hope to match the limiting form of the outer solution, which near the point *x*_{1}=*a* is
7.4a
and
7.4b
(Here we use the limiting elastic form of the micropolar equations.)

In the limit ℓ→0+, the dependent variables are rescaled as
where *X*_{α} (*α*=1,2) are as in (7.1), and *U*_{α}, *Σ*_{α}, *T*_{ij}, *Ω* are the new dependent variables. The upshot of this rescaling is that the dependence on ℓ in the governing equations is completely eliminated.

The other end of the crack at *x*_{1}=−*a* becomes *X*_{1}=−2*a*/ℓ, which disappears to (−∞) in the formal limit ℓ→0+. Hence, in the inner approximation, we have the following boundary conditions to leading order
Matching with the outer solution requires
and
The next step is to apply the complex Fourier transform to the governing equations. To this end, we define
, and note that is analytic in some upper half-plane with Im(*s*)>0, while enjoys the same property in some lower half-plane with Im(*s*)<0. After some algebra, one can show that the boundary conditions imply
this equation holding on the real *s*-axis, and
This functional equation is rearranged as
7.6
where we have written . Here has a branch cut from 0 to i∞, while has a branch cut from 0 to −i∞, both square-root functions being positive when . The function *Q*(*s*) above, defined by both sides of the equation, is analytic, except possibly at *s*=0. For |*s*|≫1 each side is bounded because of the edge conditions for the stress *T*_{21} and displacement *U*_{1}, as *X*_{1} tends to zero from the positive and negative sides, respectively. In the Fourier transformed variables, the matching conditions as *X*_{1}→±∞ become,
and
Hence we find
Note that *k*_{−}(0)*k*_{+}(0)=*k*(0)=1.

Taking the limit as in the solution of equation (7.6) in the respective half-planes of regularity, and then inverting the resulting transforms, gives
and
Hence, reverting to the original local system of coordinates at the crack tip, we get as *x*_{1}→*a*+
7.7a
and
7.7b
where *k*_{0} is given by equation (6.8) and *k*_{+}(0) is defined according to
7.8
To summarize, in the last two sections, we have succeeded in finding expressions similar to equation (2.1) of §2 for the general micropolar solid, namely,
7.9
where
and all quantities that appear above have already been defined (see also appendix A). A similar expression to (7.9) can be given for *k*_{3}, the stress-intensity factor for the *Φ*_{3}-component of the rotation vector in the micropolar theory, but we omit it in the interest of brevity.

The quantitative assessment of the formulae (7.9) is given in figure 14. As those expressions are leading-order results, the dependence on the small parameter (ℓ*/*a*) is unfortunately not captured. However, we can fix *ν* (equal to 0.1, 0.2, 0.3, 0.4, as before), and then vary . For mode I, the corresponding results appear in figure 14*a*, where the arrow indicates the direction of increasing of the Poisson's ratio. It is apparent that, as Δ≫1, both stress-intensity factor ratios asymptote towards different constant values, although neither of these is equal to unity, hence reinforcing the remarks made earlier in relation to equation (5.10). The situation is somewhat different in figure 14*b* as there the ratio of the *K*_{II} factors for the elastic and micropolar solid, respectively, tends to a limiting value *not equal to one*, as Δ→∞ (compare with the results illustrated in figure 13).

## 8. Concluding remarks

The brittle fracture of rock has been investigated by three contrasting approaches. Our attention has been directed mostly to the CBDT, which provides a robust and relatively simple means for calculating mode I and mode II fracture toughness of rock, this numerical factor being determined via an elastic theory of load transmission through the rock and the interaction of the pre-cut crack with the sample boundary. A complementary approach involved describing the fracture of those samples using a discrete model of the rock as distributions of rigid discs with specified contact rules between them. The qualitative comparison of the fracture in these discrete models with the observed behaviour in the experiments is evocative, even when the properties of the discs are not necessarily tuned to the actual rock properties. A facile explanation might be ‘well, the complete fracture is a macroscopic phenomenon, and so overrides the detailed picture of the rock, even though there are temporary local variations’.

The third approach we used was to consider the discrete model as behaving like a linear micropolar continuum, and then solve the fracture problem for the micropolar continuum. This last step is considerably more difficult since in the discrete model, one can follow a fracture path fairly easily, whereas the analogy of the couple stress (micropolar continuum) turns out to be much more complicated. Nevertheless, it was shown that by adopting the micropolar continuum framework and by using singular perturbation arguments, one can provide a quantitative comparison between the new (i.e. micropolar) stress-intensity factors and their purely elastic counterparts.

## A. Appendix A. Additional formulae

The formulae below are deduced from the paper by Atkinson & Leppington [27] and are obtained by performing the substitution (3.7) in their corresponding results; interested readers may want to refer to the aforementioned work for full details.

The values of the parameters *k*_{1} and *k*_{3} in (5.3) and (5.8), corresponding to the micropolar continuum, will be denoted by and , respectively. Lengthy arguments give the expressions of these coefficients as
A 1
and
A 2
where
The integrand in these last two integrals corresponds to
with
Here *k*_{−}(*s*) satisfies
and, finally, *L* in (A.1) corresponds to
The constant *c*_{0} that enters in the equations (7.2) and (7.3) is defined by

## Acknowledgements

C.A. and C.D.C. thank Schlumberger Gould Research Center (Cambridge, UK) for giving them permission to publish this work. The authors would also like to express their gratitude to Profs A. Martin-Meizoso and J.M. Martinez-Esnaola (San Sebastian, Spain) for their advice regarding the mechanical tests, and to Dr F. Guillard (Sydney, Australia) for several discussions regarding [39].

## Footnotes

One contribution of 19 to a theme issue ‘Fracturing across the multi-scales of diverse materials’.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.