Couple stresses and the fracture of rock.

Downloaded from http://rsta.royalsocietypublishing.org/ on April 27, 2015

Couple stresses and the fracture of rock rsta.royalsocietypublishing.org

Colin Atkinson1,2 , Ciprian D. Coman2 , Javier Aldazabal3

Research Cite this article: Atkinson C, Coman CD, Aldazabal J. 2015 Couple stresses and the fracture of rock. Phil. Trans. R. Soc. A 373: 20140120. http://dx.doi.org/10.1098/rsta.2014.0120

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

Subject Areas: materials science, mechanics Keywords: fracture mechanics, micropolar elasticity, singular perturbations, stress analysis Author for correspondence: Colin Atkinson e-mail: [email protected]

1 Department of Mathematics, Imperial College London,

South Kensington Campus, London SW7 2AZ, UK 2 Schlumberger Gould Research Center, High Cross, Madingley Road, Cambridge, UK 3 CEIT and Tecnun (University of Navarra) P. Manuel de Lardizabal 15, 20018 San Sebastian, Spain 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 closedform 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) 2015 The Author(s) Published by the Royal Society. All rights reserved.


(a)

(b)

.........................................................

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

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 373: 20140120

Figure 1. Examples of synthetic rock materials: spherical elements are used in the three-dimensional case (a), whereas discs/cylindrical elements correspond to the two-dimensional scenarios (b). This latter case is a version of the Taylor–Schneebeli material consisting of piled, small diameter, unit-length rods assumed to experience identical deformation in any transverse cross section. (Online version in colour.)

2


mpol

.........................................................

mpol

and KII ) for a finding the local, near-tip stress field (including the stress-intensity factors, KI 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 KIC and KIIC 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

3


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 KIC 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 SCFelas = 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, SCFmpol , 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


The cracked Brazilian disc test (CBDT) is used frequently in rock mechanics to investigate mixedmode 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 2a, 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 106 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 KI and KII can be expressed in the form a P a , j ∈ {I, II}, (2.1) Kj = Fj α, R Rt π where FI and FII 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 FI 1 − 4 sin2 α + 4 sin2 α(1 − 4 cos2 α) and

a 2 R

a 2 sin(2α). FII 2 + (8 cos2 α − 5) R

(2.2a)

(2.2b)

.........................................................

2. Experimental work

4


on the far-reaching analytical strategy of Atkinson & Leppington [27], we establish the leadingorder 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.


(a)

(b)

5

P

P

x1

X2

X2 2a

x2

2R

x2 a

P

P

Figure 2. (a,b) Geometry of the CBDT. The orientation of the crack (controlled by the angle α) is variable: mode I corresponds to α ≡ αI = 0◦ , whereas mode II loading is achieved for α ≡ αII , where αII is obtained by solving a trigonometric equation that depends on the ratio (a/R). Further details are included in the text. (Online version in colour.)

Figure 3. Crack-tip detail in CBDT, as used in our experiments. (Online version in colour.)

As an example, we include the variation of the shape coefficients as a function of α, for a/R = 0.5 (figure 4a) and a/R = 0.2 (figure 4b), 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 FI = 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 FI 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

.........................................................

X1


x1 X1


(a)

(b)

4

4

6

FI

FI , FII

FI , FII

2

0

−2

0

−2

0°

25°

50° a

75°

100°

0°

25°

50° a

75°

100°

Figure 4. The variation with α of the geometrical coefficients in the stress-intensity factors formula (2.1). The plots on (a) correspond to a/R = 0.5, whereas on (b) a/R = 0.2. Continuous lines are used for showing the five-term truncated result from [30], with the dashed lines representing the small-crack approximation given by equation (2.2). (Online version in colour.)

5

5

4

4 load (kN)

(b) 6

load (kN)

(a) 6

3 2

3 2

1

1

0

0 0

3 6 9 crosshead displacement (mm)

12

0


12

Figure 5. Load-displacement plots using the data recorded during the mode I (a) and mode II (b) tests on the marble samples. Each curve corresponds to a different test, and the colour coding in the two sets of plots is independent of each other. (Online version in colour.) 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 6a, 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 6b, 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

.........................................................

FII

FII 2


FI


7

(b)

1 cm

(b)

1 cm

(a)

Figure 7. Same as per figure 6, except that in this case, both halves of the rock sample moved apart suddenly within milliseconds of each other. (Online version in colour.) (b) 3

3

2

2

load (kN)

load (kN)

(a)

1

1

0

0 0


12

0


12

Figure 8. Same as per figure 5 for the sandstone tests: mode I (a) and mode II (b). (Online version in colour.) identified as the thick red curve in figure 5a. 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 9a,b, it is easy to see how the two primary cracks have

.........................................................

Figure 6. Two snapshots, just before (a) and after (b) the fracture of one of the marble samples subjected to pure mode I loading. In this case, after the fracture failure had occurred, one half of the rock disc was violently pushed away sideways, while the other continued to exert a non-zero load-bearing capacity. The load-displacement curve associated with this scenario corresponds to the thick red curve in figure 5a. (Online version in colour.)


1 cm

(a)


(a)

8

1 cm

propagated from either notch tip, extending further to the edge of the circular configuration. For mode II (figure 9b), 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 KIC and KIIC 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 KIC = 0.692 (s.d. = 0.69%) and KIIC = 0.625 (s.d. = 0.50%), whereas for Carrara marble, KIC = 1.257 (s.d. = 0.75%) and KIIC = 1.529 (s.d. = 0.55%); the unit for the stress-intensity factors is MPa · m1/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 ei (i = 1, 2, 3). However, the outline of the theory included below is given using a coordinatefree approach. The equilibrium equations for a linear micropolar continuum assume the form (e.g. [10] or [12])

V · t = 0,

V · μ + t× = 0,

(3.1)

where t and μ are the asymmetric force-stress and couple-stress second-rank tensors characterized by the general constitutive equations t = λ|ε|I + 2με − γ [I ∧ (ω − Φ)]

(3.2a)

μ = α(V · Φ)I + 4μ2∗ [V ⊗ Φ + η(Φ ⊗ V)],

(3.2b)

and 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 t (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 u by ω = − 12 V ∧ u. Other notations used above are defined as follows: ε represents the classic infinitesimal strain tensor ε ≡ 12 (V ⊗ u + u ⊗ V), |ε| is its first principal invariant, I represents the three-dimensional identity tensor, V ⊗ Φ ≡ Φ j,i ei ⊗ ej , and Φ ⊗ V ≡ (V ⊗ Φ)T , with the superscript T indicating the transposition operation for secondrank 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.

.........................................................

Figure 9. Crack paths in two representative experiments on Lazonby sandstone: mode I fracture (a) and mode II (b). (Online version in colour.)


1 cm

(b)


Table 1. Correspondence of the micropolar constants in various referenced works.

9

Atkinson & Leppington [27]

2α

γ

K

β

α

—

μ η

γ +ε

4μ2

γ

(γ − ε)/(γ + ε)

η

—

λ

λ

λ

μ

μ

μ + K /2

ν

ν

[2(μ/γ ) + 1]

c1

—

b in equation (2.12)

..........................................................................................................................................................................................................

—

..........................................................................................................................................................................................................

η

2

.......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

ν1 in equation (2.13)

..........................................................................................................................................................................................................

1/2

.......................................................................................................................................................................................................... .....................................................................

in equation (2.3)

..........................................................................................................................................................................................................

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 V · Φ = V · ω ≡ 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 b := ∗ ,

t := ∗ 2(1 + η)

and

N :=

γ , 2μ + γ

(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 couplestress 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

.........................................................

Bigoni & Drugan [34]


Nowacki [12]

..........................................................................................................................................................................................................


of the linear micropolar elasticity is readily available with the help of Mindlin’s stress potentials [24], φ ≡ φ(x1 , x2 ), ψ ≡ ψ(x1 , x2 ), defined according to (3.4)

where I2 is the two-dimensional identity tensor in the 12-plane, N(φ) := (∇ 2 φ)I2 − V ⊗ Vφ, 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 e1 + μ23 e2 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 ∇ 4 φ = 0,

2∗ ∇ 4 φ − ∇ 2 φ = 0

(3.5a)

and

V[ΓL2 (∇ 2 ψ) − ψ] = ΓR2 V(∇ 2 φ) ∧ e3 , where

ΓL2

≡

2 ,

for micropolar

2∗ ,

for couple stress

and (in our current notation)

,

(3.5b)

ΓR2 ≡ 2(1 − ν)2∗ ,

μ . 2 := 2∗ 1 + 2 γ

Note that, in the couple-stress limit (i.e. γ → ∞), we have → ∗ . In [27], the key parameters for the micropolar theory were c21 =

γ (μ + K ) , K (2μ + K )

b2 =

γ 2(2μ + K )

and ν1 =

λ , 2λ + 2μ + K

(3.6)

and it is easy to see from table 1 that the limit of (1 − ν1 )b2 /c21 as γ → ∞ is just (1 − ν), which explains why in that work the micropolar solutions were obtained from the corresponding couplestress results by performing the substitution (1 − ν) →

(1 − ν1 )b2 c21

.

(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 (Kn ) and two identical shear (tangential) springs (Ks ), 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 r0 ; such microcylinders can be seen in

.........................................................

μ = Vψ,


t = N(φ) + N(ψ) ∧ e3 ,

10


figure 1a (see also [38] for a similar concept implemented in the commercial software PFC3D). By following this idealization, the following expressions for the stiffnesses can be deduced Ks =

π r20 Gb , h

Kϕn =

π r40 Gb 2h

and Kϕs =

π r40 Eb , 4h

where Eb and Gb 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: 1 λ = f (Kn − Ks ), 5 and

3 1 μ = f Kn + Ks , 5 2

1 4μ2∗ = f (Kϕn + 4Kϕs ), 5 η=

Kϕn − Kϕs , Kϕn + 4Kϕs

γ = fKs

1 α = f (Kϕn − Kϕs ), 5

(3.8a)

(3.8b)

where f := Zvs /(π D), with Z being the coordination number, vs 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 (3.8) 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 Eb ∼ 2Gb , b ∼ D/2, Zvs ∼ 1.0, Gb ∼ 10 GPa, D ∼ 10−3 m and h/D ∼ 0.1. With these assumptions, it was found that Kn ∼ 150 MPa · m, 4μ2∗

2

∼ 3.0 kPa · m ,

Ks ∼ 75 MPa · m, μ ∼ 15 GPa,

Kϕn ∼ Kϕs ∼ 10 N · m,

λ ∼ 5 GPa,

γ ∼ 20 GPa

and α ∼ 0.

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.

.........................................................

π r20 Eb , h


Kn =

11


(a)

(b)

12

2.65

−2.65 −2.65

0 x1

2.65

−2.65

0 x1

2.65

Figure 10. The CBDT for mode II loading of the crack. (a) The failure pattern of the circular synthetic rock material and (b) the contour plot displays the distribution of the maximum principal stress, just before the fracture started to propagate. The darker regions indicate the maximum stress concentration experienced by the circular sample. (Online version in colour.) (a)

(b) 2.65

x2

0

−2.65 −2.65

0 x1

2.65

−2.65

0 x1

2.65

Figure 11. Similar to figure 10, except that here the central crack is under pure mode I loading. (Online version in colour.)

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

⎤ ⎡ 3θ θ ⎢− cos 2 + 5 cos 2 ⎥ ⎡ ⎤ ⎢ ⎥

trr θ ⎥ 1⎢ KI 3θ ⎢ ⎥ ⎢ ⎥ (5.1) + 3 cos · ⎢ cos ⎥, ⎣tθθ ⎦ = √ ⎢ ⎥ 4 2 2 2π r ⎢ ⎥

trθ ⎣ ⎦ 3θ θ + sin sin 2 2

.........................................................

0


x2


13

where E is Young’s modulus, ν denotes Poisson’s ratio and θ = 0 is directly ahead of the crack tip. The coefficient KI (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])

1 3θ θ − sin(θ) sin , (5.3a) t11 = −(1 − 2ν)k1 (2r)−1/2 cos 2 2 2

1 − 2ν 3θ θ − sin(θ) sin , (5.3b) t22 = (3 − 2ν)k1 (2r)−1/2 cos 2 2(3 − 2ν) 2

1 − 2ν 3θ θ + sin(θ) cos , (5.3c) t12 = −4(1 − ν)k1 (2r)−1/2 sin 2 8(1 − ν) 2

1 3θ , (5.3d) t21 = − (1 − 2ν)k1 (2r)−1/2 sin(θ) cos 2 2

1 θ , (5.3e) μ13 = − k2 (2r)−1/2 sin 2 2

θ 1 μ23 = k2 (2r)−1/2 cos (5.3f ) 2 2

1 θ , (5.3g) k2 (2r)−1/2 sin and Φ3 = 8γ 2 where the constants k1 and k2 are proportional to the new stress-intensity factors that reflect the addition of couple stress (and are not to be confused with KI 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 X1 and X2 in figure 2). For the elastic case KI (5.4a) t22 = √ x−1/2 , for x > 0, 2π and

8(1 − ν 2 ) KI (−x)1/2 , · [[u2 ]] = √ E 2π

for x < 0,

(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 G1 =

KI2 1−ν (1 − ν 2 ) ≡ KI2 · . E 2μ

(5.5)

In the micropolar case, by letting KI∗ denote the mode I stress-intensity factor, we find K∗ t22 = √ I x−1/2 , 2π and

for x > 0

K∗ 8(1 − ν1 )(−x)1/2 , [[u2 ]] = √ I · 2π (2μ + K )[1 + 2(1 − ν1 )b2 /c21 ]

(5.6a)

for x < 0,

(5.6b)

.........................................................

(5.2)


and the corresponding displacement field is given by

⎤ ⎡ θ 1 3θ 5 − 4ν cos − cos (1 + ν) ⎢ 2 2 2 2 ⎥ ur ⎢

⎥ , = KI 7 θ 3θ ⎦ 1 uθ E ⎣ − − 4ν sin + sin 2 2 2 2


so the work done by the field t22 for virtual crack advance will be 1 − ν1

.

(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 k3 μ23 = √ x−1/2 , 2π

for x > 0

(5.8a)

and [[Φ3 ]] =

k3 4 · √ (−x)1/2 , γ 2π

for x < 0,

(5.8b)

where k3 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 Gc = k32 /(2γ ), so that the total energy release rate is G =G∗ +

k32 . 2γ

(5.9)

The coefficients KI , KI∗ and k3 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 G given in equation (5.9) does tend to the elastic one recorded in (5.6). Thus, KI = lim KI∗ , →0

but GI = lim G . →0

(5.10)

In closing this section, we use equation (A 2) from appendix A to study the dependence of k¯ 1 and k¯ 3 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

∗ μ c1 2(1 − ν) = 1/2 , β1 = , κ0 = 1 + and ≡ 1 + 2 . (5.11) a a 1−ν γ 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 | − 1| = O(1) (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 k¯ 3 be the micropolar quantity corresponding to k3 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 k¯ 3 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).

.........................................................

(2μ + K )[1 + 2(1 − ν1 )b2 /c21 ]


G1∗ = K1∗2 ·

14


(a)

(b)

0.20

15

2.5 k1

k3 0.10 2.0 0.05

1.5 1.0

0

0.1

0.2 * /a

0.3

0

0.4

0.1

0.2 * /a

0.3

0.4

Figure 12. The variation of the micropolar stress-intensity factors k¯1 and k¯3 defined in (A 2) from appendix A, with the nondimensional parameter (∗ /a). In both windows, there are three families of curves that correspond to = 20.0 (top), = 5.0 (middle) and = 1.2 (bottom). Each one of these three sets consists of four curves obtained for ν = 0.1, 0.2, 0.3, 0.4; the arrow shows the direction in which this parameter increases. (Online version in colour.)

6. The shear crack (mode II) We consider the semi-infinite crack loaded under shear by an internal shear stress t21 = −t0 exp(x1 /a), x1 < 0 on x2 = 0. From symmetry, we have the boundary conditions t22 = σ23 = 0 on x2 = 0, (∀)x1 , and u1 = 0 on x2 = 0, x1 > 0 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 (s ∈ C) t0 a 2μ + K ¯t21+ (s) − F(s)u¯ 1− (s) on the axis Im(s) = 0, (6.1) = 1 + isa 2(1 − ν1 ) where

|s| , F(s) := −|s| 1 + 4(1 − ν1 )b2 s2 1 − β0

and ¯t21+ (s) :=

∞ 0

t21 (ζ , 0) exp(isζ ) dζ

and u¯ 1− (s) :=

β0 ≡ s 2 +

1

1/2

c21

0 −∞

u1 (ζ , 0) exp(isζ ) dζ .

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 t21 = (π x1 )−1/2

t0 a1/2 , k+ (ic1 /a)

x1 > 0, x2 = 0

(6.2)

and u1 =

4(1 − ν1 ) t0 a1/2 (−x1 )1/2 · , · √ k0 k+ (ic1 /a) 2μ + K π

x1 < 0, x2 = 0,

(6.3)

.........................................................

0.15


3.0


where

1

(ζ − is1 )

tan

−1

4ζ 3 2 (β1 − 4ζ )(1 − ζ 2 )1/2

16 dζ ,

(6.4)

and we recall that k0 ≡ 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 c1 and b involve the coefficient γ , which relates the couple stress to the gradient of the micropolar component Φ3 . Thus, (c1 /a) depends upon γ , but (b/c1 ) 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 t21 = (π x1 )−1/2 t0 a1/2 ,

x1 > 0, x2 = 0,

(6.5)

and, respectively, u1 =

(−x1 )1/2 4(1 − ν1 ) 1/2 · t0 a , √ 2μ + K π

x1 < 0, x2 = 0,

(6.6)

The above results also apply to the classical couple-stress problem in which, in place of k+ (ic1 /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 (c1 /a) or (/a) tend to zero. However, the energy release rate for this crack (assuming it propagates straight ahead) is given by 1 , (6.7) G = GE · k0 [k+ (ic1 /a)]2 where GE is the corresponding result for the elastic crack with local crack-tip stress field given by (6.5) and (6.6). Note that lim(c1 /a)→0 (k0 )1/2 k+ (ic1 /a) = 1 (also, when (/a) → 0 in the couplestress case). Let us mention briefly that due to the convention stated in equation (3.7), in the micropolar case 2(1 − ν1 )b2 . (6.8) k0 ≡ 1 + c21 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 13a shows the dependence of G /GE 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 G /GE . In (mpol)

figure 13b, the ratio between the micropolar mode II stress-intensity factor, KII

, and its elastic

(elas) KII ,

counterpart, 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 13b), the difference is less marked. Nevertheless, as figure 13b 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.

.........................................................

0

−1


1 ln{k+ (s1 )} ≡ − π


(a)

(b) 1.6

KII(mpol) / KII(elas)

E

/

0.85

0.75

1.4 1.3 1.2 1.1

0.65

0 0.05

0.15

0.25

0.35

0.45

1.0

* /a

0 0.05

0.15

0.25

* /a

0.35

0.45

Figure 13. Dependence on (∗ /a) of G /GE (a), and of the ratio of mode II stress-intensity factors for the elastic and micropolar case (b). The values of the parameters used to obtain these curves are the same as in figure 12, with the arrows pointing in the direction of increase of the Poisson’s ratio ν. (Online version in colour.)

(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 KI 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 c1 defined in (3.6) are O(), whereas the ratio c1 /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 (x1 = a, x2 = 0, say); see figure 2 for further details. One writes x1 − a = X1 ,

x2 = X2 ,

as → 0+,

(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 c0 L ∗ (7.2) KI = KI (1 − d) + 1−d and k3 = −KI

ic0 L , 1−d

(7.3)

where d and c0 are as defined in appendix A, and KI is the elastic stress-intensity factor given by (2.1). These results are strictly true only in the limit 0 < 1.

.........................................................

1.5


0.95

17


(b) Mode II

18

and u1 (x1 , 0) ∼

KII (x1 − a)−1/2 , (2π )1/2

as x1 − a → 0+

KII 4(1 − ν1 ) · (a − x1 )1/2 , 2μ + K (2π )1/2

as a − x1 → 0 + .

(7.4a)

(7.4b)

(Here we use the limiting elastic form of the micropolar equations.) In the limit → 0+, the dependent variables are rescaled as a 1/2 uα (x1 , x2 ) = (a)1/2 Uα (X1 , X2 ), τij (x1 , x2 ) = Tij (X1 , X2 ), a 1/2 Ω(X1 , X2 ), σα (x1 , x2 ) = (a)1/2 Σα (X1 , X2 ), ω(x1 , x2 ) = where Xα (α = 1, 2) are as in (7.1), and Uα , Σα , Tij , Ω 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 x1 = −a becomes X1 = −2a/, which disappears to (−∞) in the formal limit → 0+. Hence, in the inner approximation, we have the following boundary conditions to leading order Σ2 = 0

on X2 = 0, (∀)X1 ;

U1 = 0 on X1 > 0, X2 = 0;

T22 = 0

on X2 = 0, (∀)X1 ;

T21 = 0 on X1 < 0, X2 = 0.

Matching with the outer solution requires T21 (X1 , 0) ∼ and U1 (X1 , 0) ∼

KII −1/2 X , (2aπ )1/2 1

as X1 → +∞,

KII 4(1 − ν1 ) · (−X1 )1/2 , 2μ + K (2aπ )1/2

as X1 → −∞.

The next step is to apply the complex Fourier transform to the governing equations. To this end, we define ∞ 0 T21+ (s) := T21 (X1 , 0) exp(isX1 ) dX1 and U1− (s) := U1 (X1 , 0) exp(isX1 ) dX1 −∞

0

(s ∈ C), and note that T21+ is analytic in some upper half-plane with Im(s) > 0, while U1− enjoys the same property in some lower half-plane with Im(s) < 0. After some algebra, one can show that the boundary conditions imply T21+ (s) = −

2μ + K |s|k(s)U1− (s), 2(1 − ν1 )

this equation holding on the real s-axis, and k(s) := 1 + 4(1 − ν1 )

b c1

2

s2 1 −

|s| . (1 + s2 )1/2

This functional equation is rearranged as T21+ 1/2 s+ k+ (s)

=−

(2μ + K ) 1/2 s k− (s)U1− (s) =: Q(s), 2(1 − ν1 ) −

(7.6)

.........................................................

τ21 (x1 , 0) ∼


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 KII ; 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 x1 = a, x2 = 0 (figure 2), and hope to match the limiting form of the outer solution, which near the point x1 = a is


1/2 1/2

1/2

1/2

and U1− ∼ −

as s → 0,

KII 2(1 − ν1 ) −3/2 exp(iπ/4)s− , · (2a)1/2 2μ + K

Hence we find Q(s) =

A , s

with A ≡

as s → 0.

KII exp(iπ/4) . · 1/2 k+ (0) (2a)

Note that k− (0)k+ (0) = k(0) = 1. Taking the limit as |s| → ∞ in the solution of equation (7.6) in the respective half-planes of regularity, and then inverting the resulting transforms, gives T21 ∼ and U1 ∼

KII −1/2 X , (2aπ )1/2 k+ (0) 1

as X1 → 0+

4(1 − ν1 ) KII (−X1 )1/2 , · (2aπ )1/2 k0 k+ (0) 2μ + K

as X1 → 0−.

Hence, reverting to the original local system of coordinates at the crack tip, we get as x1 → a+ KII (x1 − a)−1/2 , (2π )1/2 k+ (0)

(7.7a)

4(1 − ν1 ) KII (a − x1 )1/2 , · (2π )1/2 k0 k+ (0) 2μ + K

(7.7b)

τ21 (x1 , 0) ∼ and u1 (x1 , 0) ∼

where k0 is given by equation (6.8) and k+ (0) is defined according to 4ζ 3 1 11 −1 k+ (0) = exp − dζ . tan π 0 ζ (β1 − 4ζ 2 )(1 − ζ 2 )1/2

(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, a P a mpol , j ∈ {I, II}, (7.9) = Hj (, ν)Fj α, Kj R Rt π where

⎧ ⎪ ⎨(1 − d) + c0 L , if j ≡ I, 1−d Hj (, ν) := 1 ⎪ , if j ≡ II, ⎩ k+ (0)

and all quantities that appear above have already been defined (see also appendix A). A similar expression to (7.9) can be given for k3 , 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 ∈ (1, +∞). For mode I, the corresponding results appear in figure 14a, where the arrow indicates the direction of increasing of the Poisson’s ratio. It is apparent that, as 1, both stress-intensity factor ratios

.........................................................

KII −1/2 exp(iπ/4)s+ , (2a)1/2

T21+ ∼

19


where we have written |s| = s+ s− . Here s− has a branch cut from 0 to i∞, while s+ has a branch cut from 0 to −i∞, both square-root functions being positive when s ∈ R+ . 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 T21 and displacement U1 , as X1 tends to zero from the positive and negative sides, respectively. In the Fourier transformed variables, the matching conditions as X1 → ±∞ become,


6

1.5 KI*/KI

KII(mpol) / KII(elas)

k3/ KI (K*I /KI)

5 4

k3 /KI

3 2

1.4 1.3 1.2 1.1

1 0

20

5

10

15 D

20

25

30

1.0

5

10

15 D

20

25

30

Figure 14. (a) We use the equations (7.2) and (7.3) to plot the dependence of the corresponding stress-intensity factor ratios, as increases (recall that in the limit → ∞ the micropolar effects diminish). (b) We show a similar result for the shear crack (see equations (7.4a) and (7.7a)). In both windows, the arrows show the direction in which ν increases; each set of curves corresponds to ν = 0.1, 0.2, 0.3, 0.4. (Online version in colour.)

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 14b as there the ratio of the KII 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. 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].

.........................................................

(b) 1.6


(a) 7


Appendix A. Additional formulae

21

k¯ 3 = where J

c 1

a

:=

1 2π i

1 0

c i 1 , J 1−d a

m(−iζ ) dζ , ζ + (c1 /a)

d :=

(A 2) 1 2π

1 m(−iζ ) dζ . 0

The integrand in these last two integrals corresponds to m(−iζ ) = − with β1 ≡ Here k− (s) satisfies

c21

8ζ 3/2 (1 − ζ 2 )1/2 (1 + ζ )1/2 · k− (−iζ ), β 1 − β2 ζ 2 + β3 ζ 4

b2 (1 − ν1 )

,

β2 ≡ β1 + 8,

2 . β3 ≡ 8 1 + β1

4ζ 3 1 1 k− (s) −1 −1 ≡− log (ζ + is) tan dζ k0 π 0 (1 − ζ 2 )1/2 (β1 − 4ζ 2 )

and, finally, L in (A 1) corresponds to 4ζ 3 i 1 i i 1 −1 tan ζ m(−iζ ) dζ . + dζ − L≡ π 0 2 2π 0 (1 − ζ 2 )1/2 (β1 − 4ζ 2 ) The constant c0 that enters in the equations (7.2) and (7.3) is defined by 1 1 m(−iζ ) dζ . c0 ≡ 2π i 0 ζ

References 1. Block G, Jin H. 2009 Role of failure mode on rock cutting dynamics. SPE Annu. Tech. Conf. Exhibition, New Orleans, Louisiana, 4–7 October, 2009. SPE-124870-MS. Dallas, TX: SPE. 2. Deng S. 2014 Simulation of shale-proppant interaction in hydraulic fracturing by the discrete element method. Int. J. Rock Mech. Min. Sci. 70, 219–228. (doi:10.1016/j.ijrmms.2014.04.011) 3. O’Connor RM, Torczynski JR, Preece DS, Klosek JT, Williams JR. 1997 Discrete element modeling of sand production. Int. J. Rock Mech. Min. Sci. 34, 231.e1–231.e15. 4. Potyondy D. 2012 The bonded-particle model as a tool for rock mechanics research and application: current trends and future directions. In The 7th Asian Rock Mech. Symp., Seoul, Korea, 15–19 October 2012, pp. 1–33. Seoul: The Korean Science and Technology Center. 5. Potyondy D, Cundall PA. 2004 A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 41, 1329–1364. (doi:10.1016/j.ijrmms.2004.09.011) 6. Yang B, Jiao Y, Lei S. 2006 A study on the effects of microparameters on macroproperties for specimens created by bonded particles. Eng. Comp. 23, 607–631. (doi:10.1108/026444006 10680333) 7. Sun Q, Wang G. 2013 Mechanics of granular matter. Southampton, UK: WIT Press. 8. Chang CS, Ma L. 1991 A micromechanical-based micropolar theory for deformation of granular solids. Int. J. Solids Struct. 28, 67–86. (doi:10.1016/0020-7683(91)90048-K)

.........................................................

and


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 k1 and k3 in (5.3) and (5.8), corresponding to the micropolar continuum, will be denoted by k¯ 1 and k¯ 3 , respectively. Lengthy arguments give the expressions of these coefficients as c ic1 c1 L 1 J + (1 − d) + J (A 1) k¯ 1 = − 1−d a a a


22 .........................................................


9. Suiker ASJ, deBorst R, Chang CS. 2001 Micro-mechanical modelling of granular materials. I. Derivation of a second-gradient micro-polar constitutive theory. Acta Mech. 149, 161–180. (doi:10.1007/BF01261670) 10. Dyszlewicz J. 2004 Micropolar theory of elasticity. Berlin, Germany: Springer. 11. Mindlin RD, Tiersten HF. 1962 Effects of couple-stresses in linear elasticity. Arch. Rational Mech. Anal. 11, 415–448. (doi:10.1007/BF00253946) 12. Nowacki W. 1970 Theory of micropolar elasticity. Vienna, Austria: Springer. 13. Toupin RA. 1964 Theories of elasticity with couple-stress. Arch. Rational Mech. Anal. 17, 85–112. (doi:10.1007/BF00253050) 14. Alassi HT, Holt R. 2011 Relating discrete element method parameters to rock properties using classical and micropolar elasticity theories. Int. J. Numer. Anal. Methods Geomech. 36, 1350–1367. (doi:10.1002/nag.1056) 15. Gerolymatou E. 2014 A micromechanically derived anisotropic micropolar constitutive law for granular media: elasticity. Int. J. Numer. Anal. Methods Geomech. 38, 1761–1775. (doi:10.1002/nag.2275) 16. Walton W. 1987 The effective elastic moduli of a random packing of spheres. J. Mech. Phys. Solids 35, 213–226. (doi:10.1016/0022-5096(87)90036-6) 17. Cambou B, Dubujet P, Emeriault F, Sidoroff F. 1995 Homogenization for granular materials. Eur. J. Mech. A Solids 14, 255–276. 18. D’Addetta GA, Ramm E, Diebels S, Ehlers W. 2004 A particle center based homogenization strategy for granular assemblies. Eng. Comp. 21, 360–383. (doi:10.1108/02644400410519839) 19. Moon T, Nakagawa M, Berger J. 2007 Measurement of fracture toughness using the Distinct Element Method. Int. J. Rock Mech. Min. Sci. 44, 449–456. (doi:10.1016/j.ijrmms.2006.07.015) 20. Tavarez FA, Plesha ME. 2006 Discrete element method for modelling solid and particulate materials. Int. J. Numer. Method Eng. 70, 379–404. (doi:10.1002/nme.1881) 21. Sheng Y, Yang D, Tan Y, Ye J. 2010 Microstructure effects on transverse cracking in composites laminae by DEM. Compos. Sci. Technol. 70, 2093–2101. (doi:10.1016/j.compscitech.2010.08.006) 22. Yang D, Sheng Y, Ye J, Tan Y. 2010 Discrete element modelling of the microbond test of fiber reinforced composites. Comput. Mater. Sci. 49, 253–259. (doi:10.1016/j.commatsci.2010.05.003) 23. Hedjazi L, Martin CL, Guessasma S, Della Valle G, Dendievel R. 2012 Application of the discrete element method to crack propagation and crack branching in a vitreous dense biopolymer material. Int. J. Solids Struct. 49, 1893–1899. (doi:10.1016/j.ijsolstr.2012.03.030) 24. Mindlin RD. 1963 Influence of couple-stresses on stress concentrations. Exp. Mech. 1, 1–7. (doi:10.1007/BF02327219) 25. Muki R, Sternberg E. 1965 The influence of couple-stresses on singular stress concentrations in elastic solids. J. Appl. Mech. 16, 611–648. 26. Sternberg E, Muki R. 1967 The effect of couple-stresses on the stress concentration around a crack. Int. J. Solids Struct. 3, 69–95. (doi:10.1016/0020-7683(67)90045-5) 27. Atkinson C, Leppington FG. 1977 The effect of couple stresses on the tip of a crack. Int. J. Solids Struct. 13, 1103–1122. (doi:10.1016/0020-7683(77)90080-4) 28. Radi E. 2008 On the effects of characteristic lengths in bending and torsion on mode III crack in couple-stress elasticity. Int. J. Solids Struct. 45, 3033–3058. (doi:10.1016/j.ijsolstr.2008.01.010) 29. Zhang L, Huang Y, Chen JY, Hwang KC. 1998 The mode III full-field solution in elastic materials with strain gradient effects. Int. J. Fract. 92, 325–348. (doi:10.1023/A:1007552621307) 30. Atkinson C, Smelser RE, Sanchez J. 1982 Combined mode fracture via the cracked Brazilian disc test. Int. J. Fract. 18, 279–291. 31. Atkinson C, Coman CD, Aldazabal J. In preparation. DEM modelling of mode I and mode II fracture toughness using the cracked Brazilian disc test. 32. Lakes RS. 1995 Experimental methods for study of Cosserat elastic solids and other generalized elastic continua. In Continuum models for materials with microstructure (ed. H Muhlhaus), pp. 1–22. New York, NY: J. Wiley & Sons. 33. Lakes RS. 1986 Experimental microelasticity of two porous solids. Int. J. Solids Struct. 22, 55–63. (doi:10.1016/0020-7683(86)90103-4) 34. Bigoni D, Drugan WJ. 2007 Analytical derivation of Cosserat moduli via homogenization of heterogeneous elastic materials. ASME J. Appl. Mech. 74, 741–753. (doi:10.1115/1.2711225) 35. Mühlhaus HB, Oka F. 1995 A continuum theory for random packings of elastic spheres. In Fracture of brittle, disordered materials (eds G Baker, BL Karihaloo), pp. 285–299. London, UK: E & FN Spon.


23 .........................................................


36. Pasternak E, Dyskin AV. 2014 On a possibility of reconstruction of Cosserat moduli in particulate materials using long waves. Acta Mech. 225, 2409–2422. (doi:10.1007/s00707014-1132-2) 37. Pasternak E, Mühlhaus HB, Dyskin AV. 2004 On the possibility of elastic strain localisation in a fault. Pure Appl. Geophys. 161, 2309–2326. (doi:10.1007/s00024-004-2565-7) 38. Itasca Inc. 2008 PFC3D, particle flow code. Minneapolis, MN: Itasca Consulting Group Inc. 39. Guillard F, Forterre Y, Pouliquen O. 2014 Lift forces in granular media. Phys. Fluids 26, 043301. 40. Sih GC, Liebowitz H. 1968 Mathematical theories of brittle fracture. In Fracture: an advanced treatise (ed H Liebowitz), pp. 68–188, ch. 2. New York, NY: Academic Press. 41. Diegele E, Elsaber R, Tsakmakis C. 2004 Linear micropolar elastic crack-tip fields under mixed mode loading conditions. Int. J. Fract. 129, 309–339. (doi:10.1023/B:FRAC. 0000049492.13523.5a)

Ras and the awd couple.

Porosity, permeability and 3D fracture network characterisation of dolomite reservoir rock samples.

Two carboxypeptidase counterparts from rock bream (Oplegnathus fasciatus): molecular characterization, genomic arrangement and immune responses upon pathogenic stresses.

Behaviorally-based couple therapies reduce emotional arousal during couple conflict.

Letter: The importance of the couple.

Fracture of the lateral tubercle of the posterior talar process caused by a rock-climbing fall: a case report.

The stresses and strains of tensegrity.

Aggression and love in the relationship of the couple.

Heightening in couple therapy.

Life stresses and depression in the menopause.

The right and left ventricles: the odd couple.

Protest rock and drugs.

Stresses and rewards within the psychiatrist's family.

Chip-off-the-old-rock: the study of reservoir-relevant geological processes with real-rock micromodels.

The twin condition and the couple effects on personality development.

The Chinese-Western Intercultural Couple Standards Scale.

The Norway Couple Project: Lessons Learned.

The stresses of the surgical residency.

Men together: understanding the gay couple.

End-of-life care: lifting the rock.

The parental couple relationship in child and adolescent mental health.

Couple: David Park.

Optomechanics: a strained couple.

Occupational stresses and job satisfaction.