Measurement of mixed mode interfacial strengths with cementitious materials

Engineering Fracture Mechanics xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Measurement of mixed mode interfacial strengths with cementitious materials ⁎

Maryam Tabatabaeia,b, Arash Dahi Taleghania, , Nasim Alemb a b

Department of Energy and Mineral Engineering, Pennsylvania State University, University Park, PA 16802, USA Department of Materials Science and Engineering, Pennsylvania State University, University Park, PA 16802, USA

A R T IC LE I N F O

ABS TRA CT

Keywords: Wellbore cementing Zonal isolation Rotational disk Theory of elasticity Theory of dislocation Frictional slip

In the oil and gas industry, a successful zonal isolation, which is the primary goal of cementing, is achieved by a strong bond at the cement-rock and the cement-casing interfaces. We introduce a universal test in conjunction with an analytical solution to measure the mixed mode interfacial strength of cementitious materials at the casing-cement or rock-cement interfaces. To this end, the most practical method is employed, which is the direct measurement of adhesion at the interface. This experimental setup is consisted of a rotational disk composed of two identical halves. Each half includes a cylindrical hole at its center, which can be filled by cement slurry, a rock cylinder or a steel bar. Therefore, different interfaces i.e. cement-cement, cement-steel, and cement-rock can be examined. After curing, the disk is subjected to two diametral point loads using a compressive frame. By adjusting the angle between the load direction and the overlapping interface, we can examine different combinations of modes I and II . On the other hand, an analytical elasticity solution is developed to calculate the critical shear strength that would initiate sliding between the interfaces. Analytical results show a non-uniform distribution of shear tractions along the sliding interface. Having the peak load obtained from experiment, derived formulations are used to evaluate the bonding strength of cementitious materials corresponding to different loading modes of pure shear, the combination of tension/compression and shear. The value of the current analytical derivation can be appreciated by comparison with the typical approach that assumes a uniform traction distribution along the sliding interface. We find that assuming uniform distribution will overestimate the shear strength at the cement-casing interface. The analytical solution in conjunction with the experimental results can also be used to calibrate cement interfacial friction properties.

1. Introduction In primary cementing of oil and gas wells, cement is placed in the annulus between the casing (steel pipe located inside the borehole) and the formation rock to provide an isolation between the hydrocarbon-bearing zones and the other layers penetrated by the well [1]. The safe production of the well can be achieved if the cement integrity is preserved, by preventing formation of any crack inside the cement matrix or debonding at the cement interfaces. One of the serious challenges jeopardizing the cement integrity is the failure of the cement sheaths’ bonding at either the cement-casing interface or the cement-rock interface. This debonding may result in the sustained casing pressure (SCP). United States Mineral Management Service (MMS) defines SCP as a pressure buildup at

⁎

Corresponding author. E-mail address: [email protected] (A. Dahi Taleghani).

https://doi.org/10.1016/j.engfracmech.2019.106739 Received 24 July 2019; Received in revised form 16 September 2019; Accepted 16 October 2019 0013-7944/ © 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: Maryam Tabatabaei, Arash Dahi Taleghani and Nasim Alem, Engineering Fracture Mechanics, https://doi.org/10.1016/j.engfracmech.2019.106739

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

Nomenclature d c a e P x

y α

σx

σy

tα τα

σyy τxy fΩ

bx ξ

diameter of disk radius of holes at the center of semi-cylinders depth of holes at the center of semi-cylinders length of dislocation distribution from the origin two diametral compressive point loads axis along the overlapping interface of semi-cylinders axis perpendicular to the overlapping interface angle between the load direction and the interface of two semi-cylinders normal stress component due to the external loading along x -axis normal stress component due to the external loading along y -axis normal traction due to the external loading along the interface of two semi-cylinders shear traction due to the external loading along the interface of two semi-cylinders interface normal traction due to glide dislocation shear stress component in the x − y Cartesian coordinates friction coefficient of the region Ω, Ω ≡ C for the

Bx (ξ ) δ (x ) μΩ ν EΩ

t N (x )

S (x )

φ σr σθ τrθ

cement and Ω ≡ D for the disk discrete glide dislocation Burgers vector, (bx , 0) location of the glide dislocation along the sliding interface density of the glide dislocation distribution at point ξ Dirac's delta function shear modulus of the region Ω, Ω ≡ C for the cement and Ω ≡ D for the disk Poisson’s ratio Young’s modulus of elasticity of the region Ω, Ω ≡ C for the cement and Ω ≡ D for the disk time variable of the Fourier transform total normal traction (due to the external loading and dislocation distribution) at point x along the overlapping interface total shearing traction (due to the external loading and dislocation distribution) at point x along the overlapping interface Airy stress function radial stress component in the polar coordinate tangential stress component in the polar coordinate shear stress component in the polar coordinate

the top of the well’s annulus which is potentially due to the migration of gas through the cracks and channels in the cement but behind the casing [2]. According to MMS (now known as BOEM), SCP can endanger wells life by losing hydrocarbon reserves, polluting the water column with leaking hydrocarbons, and even resulting in an underground blowout in the most serious failure of the production casing [2]. A study performed by MMS in 2003 concluded that nearly 52% of all wells in the Outer Continental Shelf (OCS) of the Gulf of Mexico exhibited SCP [3]. On the other hand, about 25% of the wells in Marcellus shale in 2011 were reported as

Sand/gravel/clay

Sand/clay

Cement sheath

Sand Casing

Fig. 1. A schematic of a wellbore with cement and casing is shown in the left picture. The right picture presents the formation of radial cracks and partial/annular micro-cracks around the casing-cement/cement-rock interfaces. 2

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

having SCP [4]. In addition, frequent pressure and temperature fluctuations in the wellbore during regular production operations may severely compromise the wellbore integrity. This issue intensifies in unconventional and injector wells, especially for wells undergoing thermal shocks and hydraulic fracturing treatments. For instance, hydraulic fracturing operations result in large pressure fluctuations which may occur hundreds of times during a single treatment. Ingraffea et al. [5] examined 41, 381 conventional and unconventional oil and gas wells drilled in Pennsylvania and reported the loss of zonal isolation at the casing-cement interface as one possible mechanism for methane migration into the atmosphere and/or into groundwater aquifers. High permeable cement sheath, formation of radial cracks and poor interfacial bonding were mentioned of some reasons failing the cement barriers [5]. Carter and Evans [6] proposed an experimental setup to measure shear and hydraulic strength bonds between cement and pipe to characterize limits of cement isolation. They defined shear bond strength as the required force to initiate sliding of the pipe with respect to the cement sheath divided by the outer surface area of the pipe. The cement hydraulic bonding was determined as the required pressure applied at the cement-pipe interface to initiate leakage [6]. Therefore, the quality of cement bonding with the casing and the surrounding rock becomes imperative to prevent gas leakage [7]. Achieving appropriate cement barriers particularly needs strong shear bonds and tensile strengths at the cement-casing and cement-rock interfaces to prevent formation of micro-annuli (interfacial cracks outside the casing) [1,8]. Fig. 1 schematically shows scenarios deteriorating cement integrity by debonding at either the cement-casing interface or the cement-rock interface, formation of radial cracks, and channelization through cement matrix. Therefore, quantification of interfacial fracture resistance is a key to assess the cement integrity. There are some studies in the literature devoted to analytically derive the stress disturbance due to the presence of partial/annular micro-cracks [9,10]. However, in this study, we present a method for measuring the mixed mode interfacial strength of the cement-casing and the cement-rock interfaces. First, an experimental setup is designed and described. Then, an analytical solution is investigated to calculate the critical shear and tensile strengths of the cementitious materials bonding. The experimental results combined with the analytical analyses provide a comprehensive investigation of the interfacial strengths and frictional properties of the cement bonding. Ayatollahi and Sistaninia [11] used a centrally cracked Brazillian disk to predict the mode II fracture toughness of rocks in terms of their mode I fracture toughness. They proposed a new mode II fracture criterion by incorporating not only the singular stress term but also the second and the third non-singular stress terms, to improve predictions of mode II fracture toughness values. Wang and Suo [12] developed an experimental method to measure the interfacial toughness of substrate/epoxy systems as a function of loading phase. They proposed Brazil-nut-sandwiches and used epoxy to glue two identical halves. Substrates were made of steel, aluminum, brass, and plexiglass. Different phases of loading were controlled by the angle θ between the interface and the diametral compression. For a homogeneous disk (in the absence of the epoxy interlayer), from fracture mechanics, we know that the mode I of fracture (pure tension) happens when θ = 0° . For the above-mentioned homogeneous disk, Wang and Suo [12] considered that the mode II of fracture (pure shear) occurs at θ ≈ 25° . Dorogoy and Banks-Sills [13] used finite difference to look at the effect of loading angle on the stress intensity factors of a centrally cracked Brazilian disk subjected to a pair of concentrated compressive loads. They [13] found that at a loading angle 23° < θ < 24° , the crack faces start coming into contact beginning at both crack tips and gradually close by approaching θ to 30°. Atkinson et al. [14] employed the dislocation theory to solve centrally cracked circular disk considering different crack lengths. They [14] obtained a range of loading angles between 21.3° and 27.2° for the pure mode II behavior of disks with the ratio of crack length to the diameter of disk varying from 0.6 to 0.3. Dong et al. [15] studied the stress distribution within both uncracked and centrally cracked circular disks under two diametral compressive forces. For the case of an uncracked circular disk, they [15] developed a series expansion for the stress distribution and listed θ = 30° for the pure mode II failure. For the central cracked circular disks with the ratio of crack length to the diameter of disk equal to 0.1, 0.2 , and 0.3, they obtained θ as 29. 67°, 28. 72°, and 27.23°, respectively. However, none of these tests were designed to measure interfacial properties as it specifically studied here. Hence, we have to also examine the effect of the frictional slip along the interface. Recently, Zhang et al. [16] extended the idea of sandwiched Brazilian disk to examine the shear bond strength between cement slurry and steel pipe. Their proposed Brazilian disk is composed of two identical halves in which the cement cylinder is fixed in one half and a steel coin is placed on the other half of the Brazilian disk. Then, they adjusted the angle between the interface and the diametral compression to be 22° to include both mode I and II . They [16] tested the effect of the residual drilling fluids as well as the surface conditions of steel pipe on the cement bond strength. In this study we first use theory of elasticity to derive the loading angle at which the shear mode failure happens for an uncracked homogeneous circular disk, obtained as θ = 30° which is in agreement with Ref. [15]. We include the effect of sliding interfaces utilizing the theory of dislocation on the derivation of the interfacial fracture mechanics of cementitious materials. Formulations are derived for cases when sliding parts are made of identical materials as well as bi-material interfaces. For the validation of this part of derivation, we consider two case studies: (1) first, identical and bi-material interfaces are separately solved. Then, formulations derived for the bi-material interfaces are solved for a special case of identical materials and compared with those directly obtained from the formulation corresponding to the case of identical materials. (2) Using the derived formulation, the problem of the frictional slip between a layer and a semi-infinite substrate subjected to a normal load is solved, and results are compared with those presented by Comninou and Schmueser [17]. Finally, experiments in conjunction with the derived formulations are used to obtain mixed mode bonding strengths of cementitious materials. Further validation of the proposed methodology is sought by comparing the shear bonding strength obtained from the current study with those available in the literature [16]. 2. Problem definition: rotational disk In this study, we develop a general-purpose technique for measuring binding at cement interfaces with casing and formation rock. 3

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

To measure mixed mode strengths of cementitious materials at the above-mentioned interfaces, the rotational disk shown in Fig. 2(a) is used. As it can be seen from Fig. 2(a), the rotational disk is composed of two semi-cylinders with a circular hole at its center. One hole is filled with the cement, and the other one can be occupied by cement, a steel piece, or rock depending on what interfacial properties need to be studied. Then, the disk is loaded using the compressive frame shown in Fig. 2(b). For more illustration, the fixture used to apply two compressive point loads P to the disk is also given in Fig. 2(c). Using theory of elasticity, it can be shown that tensile (mode I ) and pure shear (mode II ) failures happen at α = 0° and 30°, respectively. Details of the analysis are provided in the following section. Therefore, by changing the angle α within the interval of 0° and 30°, different modes of loading at the interface can be measured: pure tension, pure shear, and a combination of tension and shear. On the other hand, failure modes due to the combination of shear and compression can be examined by changing the angle α within the interval of 30–90°. This set up also provides the opportunity to examine the effect of cement-casing and cement-rock sliding interfaces. The effect of sliding of interfaces on the distribution of tractions along the interface is also analytically examined in the next section. The analytical solution developed here enables a quick assessment of cement interface properties obtained from laboratory experiments. 3. Theoretical framework In this section, an analytical solution in the elasticity framework is developed based on the testing apparatus and testing procedure. First, the stress distribution as well as the normal and shear tractions along the sliding interface are analyzed for a rotational disk subjected to two diametral compressive point loads. To this end, the principle of superposition is employed to calculate stress field and interfacial tractions. Then, the theory of dislocations is utilized to include the effect of the frictional slip of interfaces on the distribution of normal and shear tractions. The combination of test results and analytical solution is then used for determining the strength characteristics of different interfaces under mixed mode loadings. 3.1. Superposition: Theory of elasticity In this section, the stress distribution within the rotational disk under compressive point loads is analyzed. Let’s consider a circular disk with the diameter d under two equal and opposite forces P acting along the diameter, as shown in Fig. 3(a) which is the problem illustrated in Fig. 2. Using the principle of superposition, the problem (I) explained in Fig. 3(a) can be solved by the superposition of problems (II) , (III) , and (IV) illustrated in Fig. 3(b–d), respectively. Problems (II) and (III) are related to the concentrated force P acting at a point on the planar boundary of half-space. Problem (IV) is a disk under an external load X(ρ , ϕ) . This loading, X(ρ , ϕ) , is calculated in such a way that its superposition with the stress distribution at the circumference of imaginary (dashed) disks in problems (II) and (III) results in a stress distribution similar to problem (I) . The radius ρ and the angle ϕ are shown in Fig. 3(d). Solving problem (I) using the above-mentioned superposition, we can obtain the stress distribution within the circular disk under two diametral compressive point load, P in terms of α and ρ as following,

σx =

(1/2 − (ρ / d) cosα ) (1/2 + (ρ / d) cosα ) ρ 2 2P ⎧ ⎤⎫ 1 − ⎛ ⎞ sin2 α ⎡ + 2 + 1/4 − (ρ / d )cosα ) 2 2 + 1/4 + (ρ / d )cosα ) 2 ⎥ ⎬ ⎢ (( / ) (( / ) ρ d ρ d πd ⎨ d ⎝ ⎠ ⎣ ⎦⎭ ⎩

(1)

σy =

(1/2 − (ρ / d) cosα )3 (1/2 + (ρ / d) cosα )3 2P ⎧ ⎤⎫ 1−⎡ + 2 + 1/4 − (ρ / d )cosα ) 2 ⎢ (( / ) (( / )2 + 1/4 + (ρ / d )cosα )2 ⎥ πd ⎨ ρ d ρ d ⎣ ⎦⎬ ⎭ ⎩

(2)

τxy = −

(1/2 − (ρ / d) cosα )2 (1/2 + (ρ / d) cosα )2 2P ρ ⎤ ⎛ ⎞ sinα ⎡ − 2 + 1/4 − (ρ / d )cosα ) 2 ⎢ (( / ) (( / )2 + 1/4 + (ρ / d )cosα )2 ⎥ ρ d ρ d πd ⎝ d ⎠ ⎣ ⎦

(3)

in which, σx and σy are normal stress components along x - and y -directions, respectively, and τxy is the shear stress component in the π Cartesian coordinates. Considering θ = α − 2 and using Eqs. (1)–(3), the normal traction, tα and the shear traction, τα along the

Fig. 2. Geometry of the rotational disk (a) which is loaded using a compressive frame (b). The fixture used to apply two compressive point loads P to the disk is shown in (c). α is the angle between the load direction and the interface of two semi-cylinders. 4

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

Fig. 3. Stress analysis of the problem shown in (a) by the superposition of the following problems: (b) and (c) which are Boussinesq problems, halfspaces subjected to the point load, P acting on their planar boundaries, and (d) the disk under external traction, X(ρ, ϕ) . X(ρ, ϕ) is determined so that its superposition with the stresses obtained from (b) and (c) on the periphery of the imaginary (dashed) disks leads to the boundary tractions of (a).

interface with the angle α are calculated as

σx − σy σx + σy σx − σy cos2θ − τxy sin2θ = cos2α + τxy sin2α − + 2 2 2 2 2P (2cos(2α ) − 1 − 16(ρ/ d)2 + 32(ρ/ d) 4 [2 + cos(2α )] − 256(ρ / d)6cos(2α ) + 256(ρ /d )8) = πd (16(ρ / d ) 4 − 8(ρ / d )2cos(2α ) + 1)2

tα =

σx + σy

τα = −

σx − σy 2

sin2θ + τxy cos2θ =

σx − σy 2

sin2α − τxy cos2α =

4P sin (2α )(4(ρ / d )2 − 1)2 (4(ρ / d )2 + 1) πd (16(ρ / d ) 4 − 8(ρ / d )2cos (2α ) + 1)2

(4)

(5)

where the angle α is shown in Fig. 3(a). Details of solving problems (II)-(IV), respectively, as shown in Fig. 3(b–d) and the derivation of Eqs. (1)–(3) are given in Appendix A. From Eqs. (1)–(3), it is found that pure shear, mode II of failure, happens at α = 30° with the 2P corresponding shear stress of τxy = 3 πd , which is in agreement with Dong et al. [15] findings. For further illustration, the stress distribution as well as the normal and shear tractions along the interfaces with α = 0°, 30°,and 90° calculated from Eqs. (1)–(5) are presented in Fig. 4(a–c), respectively. 3.2. Sliding interfaces: Theory of dislocation In this section, the effect of sliding along the interface on the stress distribution within the rotational disk is studied. To examine the problem of the frictional slip, the slip zones are represented by continuous distributions of edge dislocations with Burgers vectors acting on the slip plane, as shown in Fig. 5. For this purpose, we employ Burgers vectors corresponding to a layer with the thickness a sliding on an elastic half-plane [17–21]. We analyze cases when the layer and substrate are identical as well as when they have different material properties. Therefore, we can study the effect of sliding at the cement-steel, cement-rock, and cement-cement interfaces on the stress distribution. In the following, the interfacial traction components, τxy and σyy , corresponding to a discrete edge dislocation with the Burgers vector bx located on the layer/substrate interface ( y = 0 plane) at x = ξ are given. For the case when layer/substrate are made of the same materials [17–20],

τxy (x , 0) =

2μbx ⎧ 1 12a2 (x − ξ ) 64a4 (x − ξ ) ⎫ x−ξ − + − 2 2 2 + (x − ξ ) 2]2 2 + (x − ξ ) 2]3 ⎬ 4 ( ) [4 [4 π (κ + 1) ⎨ x ξ a x ξ a a + − − ⎭ ⎩ 5

(6)

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

(a)

(b)

(c) Fig. 4. The stress distribution and normal and shear tractions along the interfaces with α = 0°, 30°, 90° .

Fig. 5. Cartesian coordinates: x -axis along the gliding interface with y -axis perpendicular to the interface. Antisymmetric distribution of dislocations along the gliding interface is shown with positive edge dislocations (⊥) and negative edge dislocations (⊤ ). Two different Burgers vectors Bx (ξ ) dξ and Bx∗ (ξ ) dξ are introduced for the cement region (0 ≤ |x| ≤ c ) and for the disk region (c < |x| ≤ e ), respectively.

6

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

σyy (x , 0) =

16μbx ⎧ 3a3 16a5 ⎫ − 2 + 2]2 2 + (x − ξ ) 2]3 ⎬ π (κ + 1) ⎨ [4 a + ( x − ξ ) [4 a ⎩ ⎭

(7)

and for the case when layer/substrate are different materials [21]

τxy (x , 0) =

μ1 (1 + η) bx ⎧ 2 − 1−η π (κ1 + 1)(1 − β 2) ⎨ a ⎩x − ξ

σyy (x , 0) =

μ1 (1 + η) bx ⎧−2πβδ (x − ξ ) + 1 − η a π (κ1 + 1)(1 − β 2) ⎨ ⎩

∫0

∞

(ξ − x ) t ⎫ N22 dt exp (−t )sin ⎬ a D ⎭

∫0

∞

(8)

(ξ − x ) t ⎫ N12 dt exp (−t )cos ⎬ a D ⎭

in which, t is the time variable of the Fourier transform, μ is the shear modulus, κ = 3 − 4ν for the plane strain problem andκ = for the plane stress problem, ν is the Poisson’s ratio and

(9) 3−ν 1+ν

η=

μ 2 (κ1 + 1) − μ1 (κ2 + 1) μ 2 (κ1 + 1) + μ1 (κ2 + 1)

(10)

β=

μ 2 (κ1 − 1) − μ1 (κ2 − 1) μ 2 (κ1 + 1) + μ1 (κ2 + 1)

(11)

D = (1 − β 2) exp (2t ) + (η2 − β 2) exp (−2t ) − 4(1 + β )(η − β ) t 2 − 2(η − β 2)

(12)

N12 = 4[(1 + β 2) t 2 + β ] exp (t ) − 2β (1 + η) exp (−t )

(13)

N22 = 2{2(1 + β )[(1 + β ) t − (1 − β )] t + 1 + β 2} exp (t ) − 2(η + β 2) exp (−t )

(14)

To consider an array of continuously distributed edge dislocations, bx is replaced by Bx (ξ ) dξ , in which Bx (ξ ) is the density of the dislocation distribution. We first derive tractions along the sliding interface for the case of identical materials in Section 3.2.1. Then, the derivation for bi-material interfaces is given in Section 3.2.2. To verify the presented derivation, the solution given in Section 3.2.2 is obtained for a special case of identical materials which should coincide with the results obtained from Section 3.2.1. For further validation of our mathematical approach, we also solve the frictional slip between a layer and a semi-infinite substrate subjected to a normal load. This problem has been already solved by Comninou and Schmueser [17]. Therefore, for the sake of verification, we first calculate the shear traction along the sliding interface of the above-mentioned problem and, then, compare our results with those presented in Ref. [17]. Details of the calculation are presented in Appendix B. We now turn to the derivation of tractions along the sliding interface of the rotational disk for both identical and bi-material cases. 3.2.1. Identical materials Using Eqs. (4)–(7), the total normal traction, N (x ) and the total shearing traction, S (x ) due to the external loading and dislocation distribution are obtained as

N (x ) = tα + σyy =

16μbx 2P (2cos (2α ) − 1 − 16(x / d )2 + 32(x / d ) 4 [2 + cos (2α )] − 256(x / d )6cos (2α ) + 256(x / d )8) + πd (16(x / d ) 4 − 8(x / d )2cos (2α ) + 1)2 π (κ + 1) 3a3 16a5 ⎧− ⎫ + 2 + (x − ξ ) 2]2 2 + (x − ξ ) 2]3 ⎬ ⎨ [4 a [4 a ⎩ ⎭

(15)

S (x ) = −τα + τxy =−

2μbx 4P sin (2α )(4(x / d )2 − 1)2 (4(x / d )2 + 1) + πd (16(x / d ) 4 − 8(x / d )2cos (2α ) + 1)2 π (κ + 1)

x−ξ 12a2 (x − ξ ) 64a4 (x − ξ ) ⎫ ⎧ 1 − − + 2 + (x − ξ ) 2 2 + (x − ξ ) 2]2 2 + (x − ξ ) 2]3 ⎬ ⎨ x ξ 4 a [4 a [4 a − ⎭ ⎩

(16)

Traction components should always satisfy the relation |S (x )| = f |N (x )|, which results in the following Cauchy singular integral equation,

μC ⎧ (κ C + 1) ⎨ ⎩

∫0

c

Bx (ξ ) dξ + ξ−x

∫0

c

μD ⎧ Bx (ξ )  (x , ξ ) dξ ⎫ + D ⎬ κ ( + 1) ⎨ ⎭ ⎩

∫c

where 7

e

Bx∗ (ξ ) dξ + ξ−x

∫c

e

P Bx∗ (ξ )  (x , ξ ) dξ ⎫ = −  (x , α ) ⎬ d ⎭

(17)

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

 (x , ξ ) =

12a2 (x − ξ + 2af ) 12a2 (x + ξ + 2af ) 64a4 (x − ξ + 2af ) x−ξ x+ξ 1 + + − − 2 + 2 4a + (x − ξ )2 4a + (x + ξ )2 [4a2 + (x − ξ )2]2 [4a2 + (x + ξ )2]2 [4a2 + (x − ξ )2]3 x+ξ 64a4 (x + ξ + 2af ) [4a2 + (x + ξ )2]3

−

(18)

 (x , α ) =

(16(x / d ) 4

(4(x / d )2 − 1) − 8cos (2α )(x / d )2 + 1)2

(f − 2sin (2α ) − 2fcos (2α ) + 4f (x / d )2 (5 − 2cos(2α )) + 16(x / d ) 4 (f + 2sin(2α ) − 4fcos (2α )) + 64f (x / d )6)|f = f D orf C

μD ,

fD

κ D,

μC ,

κC,

(19)

fC

f is the friction coefficient. and correspond to the material properties of disk and and correspond to the material properties of cement. f = f C when 0 ≤ |x| ≤ c , and f = f D when c < |x| ≤ e . Note that two different Burgers vectors Bx (ξ ) dξ and Bx∗ (ξ ) dξ are introduced for the cement region (0 ≤ |x| ≤ c ) and for the disk region (c < |x| ≤ e ), respectively. Considering plane stress conditions for the disk, κ =

3−ν . 1+ν

Using the relation μ =

E , 2(1 + ν )

we have

(κ C + 1) μC

=

8 EC

and

(κ D + 1) μD

=

8 , ED

where E C and ED are the

Young’s modulus of elasticity of cement and disk, respectively. By utilizing the Gauss-Chebyshev quadrature developed by Erdogan and Gupta [22] and Erdogan, Gupta and Cook [23], the following Burgers vectors corresponding to the slip at the cement interface and at the disk interface are introduced,

Bx (ξ ) = −

4P w1 (ξ ) φ1 (ξ ), 0 ≤ ξ ≤ c dE C

(20)

Bx∗ (ξ ) = −

4P w2 (ξ ) φ2 (ξ ), c < ξ ≤ e dED

(21)

where w1 (ξ ) = (1 − ξ 2)1/2 and w2 (ξ ) = (1 − ξ )−1/2 (1 + ξ )1/2 . The behavior of Burgers vector within the cement region is considered to be bounded, while that at the end points of the disk region ( x = ± e points shown in Fig. 5) is considered to be singular. Introducing Eqs. (20) and (21) into Eq. (17) and solving it, we can calculate Bx (ξ ) and Bx∗ (ξ ) . The details about solving the singular integral Eq. (17) is given in Appendix C. Subsequently, shear tractions on both sides of the sliding interfaces are calculated as

S (sm) =−

P⎧2 1 (sm, α )| f C = 0 + d ⎨π ⎩

n

n

2

∑ ⎡⎢ (1 − ri ) φ1 (ri ) ⎛

1 − 1 (ri, sm )| f C = 0 ⎞ ⎤ − ⎠⎥ ⎝ sm − ri ⎦

⎜

i=1

⎣ n+1

⎟

∑ ⎡⎢ 2(1 + i) φ2 (i) ⎛

1 ⎫ + 2 (i, m )| f C = 0 ⎞ ⎤ , 0 ≤ x ≤ c ⎠⎥ ⎝ i − (δ1/ δ2 ) m + σ / δ2 − δ1/ δ2 ⎦⎬ ⎭

⎜

i=1

⎣ 2n + 1

⎟

(22)

S (m ) =−

P⎧2  2 (m , α )| f D = 0 + d ⎨π ⎩

n

2

n

∑ ⎡⎢ 2(1 + i) φ2 (i) ⎛

1 − 3 (i, m )| f D = 0 ⎞ ⎤ − ⎠⎥ ⎝ m − i ⎦ ⎜

i=1

⎣ 2n + 1

⎟

∑ ⎡⎢ (1 − ri ) φ1 (ri ) ⎛

1 ⎫ +  4 (ri, sm )| f D = 0 ⎞ ⎤ , c < x ≤ e ⎠⎥ ⎝ ri + 1 − (δ2/ δ1 ) sm − σ / δ1 ⎦⎬ ⎭ ⎜

i=1

⎣ n+1

⎟

(23)

3.2.2. Bi-material configuration Employing Eqs. (4), (5), (8) and (9), the sliding tractions along the interface of bi-materials are derived as

N (x ) = tα +  (x ),

(24)

where

{

1−η

}

⎧ ⎪C −βδ (x − ξ ) Bx (ξ ) dξ + 2πa K1 (ξ , x ) Bx (ξ ) dξ , 0 ≤ |x| ≤ c  (x ) = 3a3 16a5 ⎨ 16μ D − B ∗ (ξ ) dξ + 2 B ∗ (ξ ) dξ , c < |x| ≤ e ⎪ π (κ D + 1) [4a2 + (x − ξ )2]2 x [4a + (x − ξ )2]3 x ⎩

(25)

S (x ) = −τα +  (x )

(26)

{

}

and

in which, 8

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

{

}

⎧C 1 Bx (ξ ) dξ − 1 − η K2 (ξ , x ) Bx (ξ ) dξ , 0 ≤ |x| ≤ c 2πa ⎪ π x−ξ  (x ) = Bx∗ (ξ ) dξ x−ξ 12a2 (x − ξ ) ⎨ 2μ D ∗ ∗ ⎪ π (κ D + 1) x − ξ − 4a2 + (x − ξ )2 Bx (ξ ) dξ + [4a2 + (x − ξ )2]2 Bx (ξ ) dξ − ⎩

{

C=

64a4 (x − ξ ) [4a2 + (x − ξ )2]3

}

Bx∗ (ξ ) dξ , c < |x| ≤ e ,

2μ1 (1 + η) (κ1 + 1)(1 − β 2) ∞

K1 (ξ , x ) =

∫0

K2 (ξ , x ) =

∫0

∞

(27)

(28)

(ξ − x ) t N12 dt exp (−t )cos a D

(29)

(ξ − x ) t N22 dt . exp (−t )sin a D

(30) c

Satisfying the relation |S (x )| = f |N (x )| and using ∫−c Bx (ξ ) dξ = 0 , we obtain

1 C⎧ ⎨ π ⎩

∫−c Bxx (−ξ ) ξdξ c

⎧1 ⎨ ⎩π

∫c

e

−

1−η 2πa

Bx∗ (ξ ) dξ 1 − x−ξ π

D

∫−c K2 (ξ , x ) Bx (ξ ) dξ − 12−πaη f ∫−c K1 (ξ , x ) Bx (ξ ) dξ ⎫⎬ + (κ D2μ+ 1)

∫c

c

c

⎭

e

 (x ,

2P = ⎛ ⎞  (x , α ) ⎬ πd ⎝ ⎠ ⎭

ξ ) Bx∗ (ξ ) dξ ⎫

(31)

where  (x , ξ ) and  (x , α ) have been already given in Eqs. (18) and (19), respectively. Burgers vectors are considered as following

Bx (ξ ) =

P w1 (ξ ) φ1 (ξ ), 0 ≤ |ξ| ≤ c d

(32)

Bx∗ (ξ ) =

P w2 (ξ ) φ2 (ξ ), c < |ξ| ≤ e d

(33)

in which, w1 (ξ ) and w2 (ξ ) are equal to those used in Section 3.2.1. By substituting Eqs. (32) and (33) into the Eq. (31) and solving the resultant equation, the shear traction is calculated along the sliding interface for the case when different materials are placed in halves of the rotational disk,

S (sm) =−

2P P 1 (sm, α )|f = 0 + C πd d

n

n

∑ (1 − ri ) φ1 (ri ) ⎧

∑ 2(1 + i) φ2

i=1

i=1

2

n+1

2μD P c (1 − η) 1 − K21 (ri, sm) ⎫ + D ⎬ ⎨ (κ + 1) d 2a ⎭ ⎩ sm − ri

2n + 1

1 (i) ⎧ − 2 (i, m )|f = 0 ⎫ , 0 ≤ |x| ≤ c ⎬ ⎨ ⎭ ⎩ (c /(δ2 d )) m − i − σ / δ2

(a) Calculated from the derivation for indentical m

(34)

lated from the derivation for bi-materials.

Fig. 6. Shear traction along the sliding cement-cement interface. 9

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

S (m ) =− P d

2P P  2 (m , α )|f = 0 + C πd d n

∑ i=1

n

2

2μD c (1 − η) 1 K22 (ri, sm) ⎫ + D − ⎬ ⎨ (κ + 1) 2a ⎭ ⎩ ((δ2 d )/ c ) sm + σd/ c − ri

∑ (1 − ri ) φ1 (ri ) ⎧ i=1

n+1

2(1 + i) 1 φ (i) ⎧ − 3 (i, m )|f = 0 ⎫ , c < |x| ≤ e ⎬ 2n + 1 2 ⎨ − i  m ⎭ ⎩

(35)

Our approach for solving Eq. (31) is presented in Appendix D. As a validation, the formulation presented here for bi-materials is solved for identical materials and the obtained results are compared with those calculated from the formulation derived for identical materials in Section 3.2.1. To this end, it is assumed that both sides of the rotational disk are filled with the cement, and the mechanical properties of system are determined by μ1 = μ 2 = 2.17 , κ1 = κ2 = 2.4 , μD = 78.74 , and κ D = 1.92. The geometry of the rotational disk, as shown in Fig. 5, is determined by d = 5.08cm , a = 1.00cm , and c = 1.27cm . The distribution of shear traction is calculated along the slip plane with α = 30° . Considering f C = 0.25 and f D = 0.1, the shear traction is obtained from Eqs. (22) and (23) and plotted in Fig. 6(a), and is calculated from Eqs. (34) and (35) and shown in Fig. 6(b). As it is seen, a reasonable agreement is obtained between Fig. 6(a) and (b). For the sake of comparison, the shear traction corresponding to the mode II failure, calculated from Eq. (5), is also included in Fig. 6(a) and (b). For the purpose of showing a numerical example, the shear traction along the sliding cement-steel and the cement-sandstone interfaces is also calculated from Eqs. (34) and (35). The mechanical properties of cement and steel are considered to be as same as those used to obtain results shown in Fig. 6. The mechanical properties of sandstone are determined by μ1 = 4.17 and κ1 = 2.33. The distribution of shear traction along the slip plane for the cases of the cement-steel interface with α = 20°, f C = 0.25 and f D = 0 and the cement-rock interface with α = 25°, f C = 0.3 and f D = 0.1 are calculated and plotted in Fig. 7(a) and (b), respectively. 4. Experimental framework 4.1. Experimental set-up and measurement Cement slurries are prepared using API class-G cement with the water to cement weight ratio of 0.44 . To mimic the bottom-hole conditions, cement samples are cured under 190° F and 3000 psi for 24 hr . Then, the rotational disk is subjected to a diametral compression with the displacement control loading at the rate of 0.1 mm/min . To this end, the unconfined compression frame shown in Fig. 2(b) is employed. 4.2. Results and discussion Using the fixture shown in Fig. 2(c), two diametral compressive point loads, P are applied to the rotational disk at different angles with respect to the direction of the overlapping interface to adjust different modes of loading. Modes I , II , and the mixed mode of I and II happen when α = 0° , α = 30° , and 0° < α < 30° , respectively. After adjusting α , the system is loaded to measure the maximum value of P . In the following, the interfacial shear strengths of the cement-cement and the cement-casing interfaces are examined experimentally. Therefore, α is adjusted to be 30°. Fig. 8(a) and (b) respectively illustrate the behavior of cement-cement and cement-steel casing bonding under mode II of loading. From Fig. 8(a) and (b), the maximum point load P is obtained as 33.85 and 6.84 kN , respectively. Then, the ultimate shear bonding strength can be calculated from Eq. (22) for the cement-cement interface and from Eq.

(a)

,

(b)

,

Fig. 7. Shear traction along the slip plane for the (a) cement-steel interface with α = 20° and (b) the cement-rock interface with α = 25° . 10

Engineering Fracture Mechanics xxx (xxxx) xxxx

35

7

30

6

25

5

20

4

P (kN)

P (kN)

M. Tabatabaei, et al.

15

3

10

2

5

1

0

0 0

0.5 1 Displacement (mm)

1.5

0

0.05

(a) cement-cement bonding,

0.1 0.15 0.2 Displacement (mm)

0.25

(b) cement-casing bonding,

Fig. 8. Interfacial shear strength. Table 1 Shear bond strength at different interfaces: cement-cement, cement-casing, and cement-rock. Interface

α = 30° S (x = 0),

Cement-Cement Cement-Casing

MPa

1.470 0.271

(34) for the casing-cement interface. Mechanical properties of the above-mentioned set-ups are considered the same as those solved and presented in Figs. 6 and 7. Shear bond strengths calculated from Eqs. (22) and (34) are given in Table 1. As mentioned earlier in the Introduction section, Zhang et al. [16] measured the shear bond strength at the cement-steel interface considering different surface roughness for the steel casing including polished, sandblasted, and rusted. They prepared cement slurries using class H cement with the water to cement weight ratio of 0.37 . Their prepared cement slurry was composed of antifoam, dispersant, fluid loss liquid additive, retarder, and potassium chloride. Zhang et al. [16] conducted experiments with the load angle α = 22° (α has shown in Fig. 3) including both failure modes of I and II . In the current study, we derived that the pure shear mode happens at α = 30° and considered the shear strength at this load angle. On the other hand, as mentioned in Section 4.1, we cure cement at the bottom hole conditions, 190°F and 3000 psi for 24 hr , while the cement is not cured in Ref. [16]. We employ Eqs. (22) and (34) to obtain the shear bond strength of cement for the identical and bi-material configurations, respectively. As shown in Figs. 6 and 7, the distribution of the shear traction along the interface is not uniform. Zhang et al. [16] formulated the shear bond strength considering a uniform distribution of shear traction along the interface as follows,

S=

P cosα πc 2

(36)

where c and P have been illustrated in Fig. 5. In the following, the shear strength obtained from current study considering a nonuniform distribution and Ref. [16] considering a uniform distribution are compared. As given in Table 1, we obtain S = 0.271 MPa for the shear bond strength at the cement-casing interface. Zhang et al. [16] reported S = 0.45 ± 0.2 MPa for the corresponding shear bond strength when the surface condition of steel is polished, equivalent to our experiments. It is found that considering a uniform distribution of shear traction along the interface will overestimate the magnitude of the shear bond strength in comparison with the consideration of a non-uniform distribution. 5. Concluding remarks An experimental set-up in conjunction with an analytical solution is introduced to measure the mixed mode interfacial strength of cementitious materials at both casing-cement and rock-cement interfaces. The idea is benefitting from the stress distribution within a circular disk subjected to two diametral compressive point loads. The disk is made of two identical halves in such a way that by rotating the overlapping interface with respect to the loading direction, different modes of loading, particularly, shear and combination of shear and tension can be studied at the interface. Each half is composed of an empty cylinder at its center which can be filled by cement, steel bar, or rock. Therefore, different interfaces of cement-cement, cement-steel casing, and cement-rock, can be examined. A compressive frame is used to apply diametral point loads to the disk and to measure the peak load, P where the identical halves start sliding. An analytical solution is developed using dislocation theory to consider the effect of the frictional slip on the mixed mode interfacial strength of cementitious materials. Finally, by having P and using analytical formulation, the shear bonding strength is obtained. The validation of the presented formulations is verified by comparing the results obtained from current study with those presented in the literature. Results given in the literature are based on considering a uniform distribution of shear traction 11

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

along the sliding interface, while the current developed formulation results in a non-uniform shear traction along the overlapping interface. We found that consideration of a uniform interfacial shear traction will end up with a larger shear strength than the consideration of a non-uniform shear traction. Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgement The authors would like to thank the college of Energy and Mineral Sciences at the Pennsylvania State University to support them during this research. Appendix A. Solution of problem (I) illustrated in Fig. 3(a) Applying Boussinesq solution [24] to the problems (II) and (III) , it is found that any point, M inside the semi-infinite region (and, consequently, at points on the periphery of the imaginary disks) is subjected to a simple compression, σr in the radial direction. Therefore, the stress components corresponding to the problem (II) are

σr = −

2P cosθ , σθ = τrθ = 0 π r

(A-1)

where σr , σθ , and τrθ are radial, tangential, and shear stresses, respectively. Similarly, the stress components for the problem (III) are

σr = −

2P cosθ1 , σθ = τrθ = 0 π r1

(A-2)

Noting that r , θ , r1, and θ1 are illustrated in Fig. 3(b) and (c). When point M is located on the periphery of imaginary disks, r and r1 becomes normal to each other and, subsequently, the following relationship can be used

cosθ cosθ1 1 = = r r1 d

(A-3)

Therefore, two equal compressive stresses of magnitude 2P/(πd ) are obtained at any point M on the circumference of imaginary disks. Since point M is an arbitrary point on the periphery of imaginary disks, we can conclude that a normal compressive force of constant magnitude 2P/(πd ) should be applied to the periphery of the disk to reproduce the same stress distribution induced by the pair of the point loads. Assuming zero external forces at the circumference of the real disk shown in Fig. 3(a), the normal compression, obtained from problems (II) and (III) , should be superposed with a uniform tension of the same magnitude 2P/(πd ) . Consequently, X(ρ , ϕ) corresponding to problem (IV) is obtained as the normal tensile traction of the uniform intensity 2P/(πd ) . Employing the theory of the Airy stress function [25] to the problem (IV) , the general Airy stress function, φ corresponding to Axisymmetric problems is written as follows

φ = Alogr + Br 2logr + Cr 2 + D

(A-4)

in which, four constants of integration are determined from the boundary conditions. The stress components corresponding to this Airy stress function are

σr = σθ =

1 ∂φ A = 2 + B (1 + 2logr ) + 2C r ∂r r ∂ 2φ ∂r 2

=−

(A-5)

A + B (3 + logr ) + 2C r2

(A-6) (A-7)

τrθ = 0

r 2logr

result in infinite stress components when r = 0 . Therefore, As it is seen from Eqs. (A-5) and (A-6), two terms of logr and constants A and B vanish, leading to σr = σθ = 2C which is a constant. Finally, the components of stress can be calculated as

σr =

2P 2P , σθ = , τrθ = 0 πd πd

(A-8)

The stress components in the Cartesian coordinates can be written by transforming stress components from polar coordinate system using the following relations,

σx = σrcos 2 θ + σθsin2 θ − 2τrθsinθcosθ σy = σr

sin2

θ+

σθcos 2

(A-9)

θ + 2τrθsinθcosθ

(A-10)

τxy = (σr − σθ )sinθcosθ + τrθ (cos 2 θ − sin2 θ)

(A-11) 12

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

where σx and σy are the normal stress components along x - and y -directions, respectively, and τxy is the shearing stress component in the Cartesian coordinates [25]. Applying Eqs. (A-9) through (A-11), the stress components in the Cartesian coordinates corresponding to the problems (II − IV ) can be derived in the following forms, Problem (II) :

σx = −

2P cosθ 2 sin θ π r

(A-12)

σy = −

2P cosθ 2 cos θ π r

(A-13)

2P cosθ cosθsinθ π r

(A-14)

τxy =

Problem (III) :

σx = −

2P cosθ1 2 sin θ1 π r1

(A-15)

σy = −

2P cosθ1 2 cos θ1 π r1

(A-16)

2P cosθ1 cosθ1sinθ1 π r1

(A-17)

τxy = −

Problem (IV) :

σx =

2P 2P 2 2P cos 2 ϕ + sin ϕ = πd πd πd

(A-18)

σy =

2P 2 2P 2P sin ϕ + cos 2 ϕ = πd πd πd

(A-19)

τxy = 0

(A-20)

Finally, the superposition of problems (II) , (III) , and (IV) leads to the following stress components in the Cartesian coordinates for the problem (I) , Problem (I)≡ Problem (II)+ Problem (III)+ Problem (IV) :

2P cosθ 2 cosθ1 2 ⎫ 2P sin θ + sin θ1 + σx = −⎛ ⎞ ⎧ ⎬ r1 πd ⎝ π ⎠⎨ ⎩ r ⎭

(A-21)

2P cos3 θ cos3 θ1 ⎫ 2P + + σy = −⎛ ⎞ ⎧ ⎬ r r πd ⎝ π ⎠⎨ 1 ⎩ ⎭

(A-22)

2P cos2 θ cos2 θ1 τxy = ⎛ ⎞ ⎧ sinθ − sinθ1⎫ ⎬ r1 ⎝ π ⎠⎨ ⎩ r ⎭

(A-23)

Using Eqs. (A-21)–(A-23), the stress state at the center of the disk can be obtained as

Fig. A1. (a) Mohr circle corresponding to the stress state at the center of disk, (b) plane of shear mode II. 13

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

σx =

2P 6P , σy = − , τxy = 0 πd πd

(A-24)

For more clarification, the Mohr circle corresponding to the stress state at the center of the disk is given in Fig. A1(a). As it can be 2P seen from Fig. A1(a), the shear mode occurs at β = 60° which corresponds to α = 30° in Fig. A1(b). τps = 3 πd presents the pure shear ° on the plane α = 30 in Fig. A1(b) where pure mode II loading for the fracture ensues. From Fig. A1, it can also be found that at any angle α between 0 and 30° (0 < β < 60° ), the load P induces both tension and shear components. Normal and shear tractions along the interface with the Angle α The normal and shear tractions along the interface with the angle α (Fig. A1(b)) is presented herein. To this end, the geometry given in Fig. A1(b) is utilized to derive the following identities,

()

⎧r = d ρ d ⎪ ⎪ ⎪ ⎪ sinθ = − ⎨ ⎪ ⎪ ⎪ cosθ = ⎪ ⎩ ⎧ r1 = d ⎪ ⎪ ⎪ ⎪ sinθ1 =

2

+

1 4

−

( ) cosα ρ d

ρ sinα d 2 ρ ρ 1 + − cosα 4 d d

()

()

()

ρ 1 − cosα 2 d 2 ρ ρ 1 + − cosα 4 d d

()

() ρ 2 d

()

⎨ ⎪ ⎪ ⎪ cosθ1 = ⎪ ⎩

()

1 4

+

+

(A-25)

( ) cosα ρ d

ρ d

( ) sinα ( ) + + ( ) cosα + ( ) cosα ( ) + + ( ) cosα ρ 2 d

1 2 ρ 2 d

ρ d

1 4

ρ d

1 4

ρ d

(A-26)

Replacement of Eqs. (A-25) and (A-26) into the Eqs. (A-21)–(A-23), relations of (4) and (5) are derived. Appendix B. Proof – frictional slip at the layer/substrate interface due to a normal load A layer of thickness a located on a semi-infinite substrate is under a uniform compressive traction with magnitude p0 on the free surface ( y = a ), as shown in Fig. B1. A tensile concentrated force P applied on the free surface is increased to emerge the slip of layer on the substrate. Due to the symmetry of the problem, slip zones are symmetrically located on b ≤ |x| ≤ c . The mechanical properties of layer and substrate are described with the same shear modulus μ and Poisson’s ratio ν . Employment of the Flamant’s solution [25], the tractions along the layer/substrate interface ( y = 0 ) are obtained as 0 σxy (x , 0) = −

2P a3x 2 πa (a + x 2)2

0 σyy (x , 0) = −p0 +

(B-1)

a4

2P πa (a2 + x 2)2

(B-2)

Introducing a continuous distribution of edge dislocations with the density of Bx (ξ ) at y = 0 and utilizing Eqs. (6) and (7) result in the total normal traction, N (x ) and the total shear traction, S (x ) as follows

Fig. B1. A layer of thickness a located on a substrate, and the expected slip zones. 14

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

0 S (x ) = σxy (x , 0) + τxy (x , 0) = −

2μBx (ξ ) dξ ⎧ 1 x−ξ 12a2 (x − ξ ) 64a4 (x − ξ ) ⎫ 2P a3x − + + − 2 2 2 + (x − ξ ) 2]2 2 + (x − ξ ) 2]3 ⎬ πa (a2 + x 2)2 π (κ + 1) ⎨ x ξ 4 a ( x ξ ) [4 a [4 a + − − ⎭ ⎩ (B-3)

0 N (x ) = σyy (x , 0) + σyy (x , 0) = −p0 +

a4

16μBx (ξ ) dξ 2P + πa (a2 + x 2)2 π (κ + 1)

3a3

16a5

⎧− ⎫ + 2 2 2 ⎨ [4a2 + (x − ξ )2]3 ⎬ ⎭ ⎩ [4a + (x − ξ ) ]

(B-4)

Satisfying |S (x )| = f |N (x )| ends up with the following relation, c

∫−c

Bx (ξ ) dξ + x−ξ

∫−c ⎡⎢− 4a2 +x −(x ξ− ξ )2 c

+

12a2 (x − ξ + 2af ) 64a4 (x − ξ + 2af ) ⎤ − Bx (ξ ) dξ 2 2 2 [4a + (x − ξ ) ] [4a2 + (x − ξ )2]3 ⎥ ⎦

⎣ 2Pa2 (af + x ) ⎤ π (κ + 1) ⎡ = , 0 ≤ |x| ≤ c −fp0 + ⎢ ⎥ 2μ π (a2 + x 2)2 ⎦ ⎣

(B-5)

where f is the firction coefficient at the layer/substrate interface. Eq. (B-5) can be rewritten in the following form

∫b

c

Bx (ξ ) dξ + ξ−x

∫b

c

K (x , ξ ) Bx (ξ ) dξ =

2Pa2 (af + x ) ⎤ π (κ + 1) ⎡ ,b≤x≤c fp0 − ⎢ ⎥ 2μ π (a2 + x 2)2 ⎦ ⎣

(B-6)

where

K (x , ξ ) =− +

12a2 (x − ξ + 2af ) 12a2 (x + ξ + 2af ) 64a4 (x − ξ + 2af ) x−ξ x+ξ 1 + − − + 2 + 2 2 2 2 2 2 2 2 2 x+ξ 4a + (x − ξ ) 4a + (x + ξ ) [4a + (x − ξ ) ] [4a + (x + ξ ) ] [4a2 + (x − ξ )2]3 64a4 (x + ξ + 2af ) [4a2 + (x + ξ )2]3

(B-7)

Considering Bx as a bounded function

Bx (r ) =

P (κ + 1) w (r ) φ (r ) 2μa

(B-8)

with

w (r ) = (1 − r 2)1/2

(B-9)

the discretized form of Eq. (B-6) becomes n

∑ i=1

fp a 1 − ri2 f + δsm + σ 1 2 + K (ri, sm) ⎤ = 0 − , m = 1, ⋯, n + 1 φ (ri ) ⎡ ⎥ ⎢ P π [1 + (δsm + σ )2]2 n+1 ⎦ ⎣ ri − sm

(B-10)

where

ri = cos(iπ /(n + 1)), i = 1, 2, ⋯, n

(B-11)

sm = cos(π (2m − 1)/(2(n + 1))), m = 1, 2, ⋯, n + 1

(B-12)

Fig. B2. Comparison of the shearing tractions for the compressive force with δ = 0.2 and f = 0.25. 15

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

K (ri, sm) 12(sm − ri + 2f / δ ) 12(ri + sm + 2σ / δ + 2f / δ ) sm − ri ri + sm + 2σ / δ 1 − 2 + + − 2 ri + sm + 2σ / δ δ [4/ δ 2 + (sm − ri )2]2 δ [4/ δ 2 + (ri + sm + 2σ / δ )2]2 4/ δ 2 + (sm − ri )2 4/ δ 2 + [ri + sm + 2σ / δ ]2 64(sm − ri + 2f / δ ) 64(ri + sm + 2σ / δ + 2f / δ ) + 4 + 4 δ [4/ δ 2 + (sm − ri )2]3 δ [4/ δ 2 + (ri + sm + 2σ / δ )2]3 (B-13)

=−

δ = (c − b)/(2a)

(B-14)

σ = (c + b)/(2a)

(B-15)

Then, φ (ri ), i = 1, ⋯, n , are obtained by solving Eq. (B-10). Having φ (ri ) , the shear traction, S (x ) given in Eq. (B-3) can be calculated from the following relation

S (sm) = −

P⎧2 (δsm + σ ) + a ⎨ π [1 + (δsm + σ )2]2 ⎩

n

∑ i=1

1 − ri2 1 ⎫ + K (ri, sm )|f = 0 ⎤ , m = 1, ⋯, n + 1. φ (ri ) ⎡ ⎥ ⎢ − n+1 r s i m ⎦⎬ ⎣ ⎭

(B-16)

Shearing traction for the case of compressive force with δ = 0.2 and f = 0.25 is calculated from Eq. (B-16) and compared with the corresponding result presented by Comninou and Schmueser [17] in Fig. B2. Appendix C. Solution of Eq. (17) and the derivation of Eqs. (22) and (23) Eq. (17) can be rewritten as

∫0

Bx (ξ ) ED dξ + C ξ−x E

c

∫c

e

Bx∗ (ξ ) dξ + ξ−x

∫0

c

Bx (ξ )  (x , ξ ) dξ +

ED EC

∫c

e

Bx∗ (ξ )  (x , ξ ) dξ = −

8 P  (x , α ), EC d

(C-1)

To solve the singular integral Eq. (C-1), we first rewrite the equation in a dimensionless form and then employ the method developed by Erdogan and Gupta [22]. Therefore, the following changes of variables are introduced: x

⎧ d = δ1 (s + 1), 0 ≤ x ≤ c ⎨ ξ = δ1 (r + 1), 0 ≤ ξ ≤ c ⎩d

(C-2)

where

δ1 = x ⎧d ⎨ξ ⎩d

c 2d

(C-3)

= δ2 s + σ , c < x ≤ e = δ2 r + σ , c < ξ ≤ e

(C-4)

where

δ2 =

e−c e+c , and σ = . 2d 2d

(C-5)

Introducing Eqs. (C-2)–(C-5) into Eq. (C-1) leads to the following equations,

1 π

D

∗

D

Bx (r ) E E dr + ∫−1 Br x−(rs) dr + π1 ∫−1 Bx (r ) 1 (r , s) dr + πE ∫ Bx∗ (r ) 2 (r , s) dr C ∫ −1 r − (δ1/ δ2 ) s + σ / δ2 − δ1/ δ2 πE C −1 1

=−

1

1

1

8P 1 (s, α ), 0 ≤ x ≤ c πE Cd

(C-6)

and

1 π

∗

C

C

E Bx (r ) E dr + ∫−1 Br x−(rs) dr + π1 ∫−1 Bx∗ (r ) 3 (r , s) dr + πE ∫ Bx (r )  4 (r , s) dr D ∫ −1 r + 1 − (δ2/ δ1 ) s − σ / δ1 πED −1 1

=−

1

1

1

8P  2 (s, α ), c < x ≤ e πEDd

(C-7)

where

1 (r , s ) =

1 s−r r+s+2 12a2 (s − r + 2af C /(δ1 d )) 12a2 + 2 2 2 − 2 2 2 − 2 2 + 2 2 2 2 2 2 2 2 2 r+s+2 4a /(δ1 d ) + (s − r ) 4a /(δ1 d ) + (r + s + 2) δ1 d (4a /(δ1 d ) + (s − r ) ) δ1 d (r + s + 2 + 2af C /(δ1 d )) 64a4 (s − r + 2af C /(δ1 d )) 64a4 (r + s + 2 + 2af C /(δ1 d )) + 4 4 − 4 4 2 2 2 2 2 2 2 2 2 3 (4a /(δ1 d ) + (r + s + 2) ) δ1 d (4a /(δ1 d ) + (s − r ) ) δ1 d (4a2 /(δ12 d 2) + (r + s + 2)2)3

16

(C-8)

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al.

2 (r , s ) =

1 (δ / δ ) s + δ1/ δ2 − r − σ / δ2 (δ / δ ) s + δ1/ δ2 + r + σ / δ2 + 2 2 21 2 − 2 2 21 2 (δ1/ δ2 ) s + δ1/ δ2 + r + σ / δ2 4a /(δ2 d ) + ((δ1/ δ2 ) s + δ1/ δ2 − r − σ / δ2 )2 4a /(δ2 d ) + ((δ1/ δ2 ) s + δ1/ δ2 + r + σ / δ2 )2 −

12a2 ((δ1/ δ2 ) s + δ1/ δ2 − r − σ / δ2 + 2af C /(δ2 d )) 12a2 ((δ1/ δ2 ) s + δ1/ δ2 + r + σ / δ2 + 2af C /(δ2 d )) 64a4 + 2 2 + 4 4 δ22 d 2 [4a2 /(δ22 d 2) + ((δ1/ δ2 ) s + δ1/ δ2 − r − σ / δ2 )2]2 δ2 d [4a2 /(δ22 d 2) + ((δ1/ δ2 ) s + δ1/ δ2 + r + σ / δ2 )2]2 δ2 d

((δ1/ δ2 ) s + δ1/ δ2 − r − σ / δ2 + 2af C /(δ2 d)) 64a4 ((δ1/ δ2 ) s + δ1/ δ2 + r + σ / δ2 + 2af C /(δ2 d )) − 4 4 2 2 2 2 3 [4a /(δ2 d ) + ((δ1/ δ2 ) s + δ1/ δ2 − r − σ / δ2 ) ] δ2 d [4a2 /(δ22 d 2) + ((δ1/ δ2 ) s + δ1/ δ2 + r + σ / δ2 )2]3

(C-9)

3 (r , s ) =

s−r r + s + 2σ /δ2 1 12a2 (s − r + 2af D /(δ2 d )) 12a2 + 2 2 2 − 2 2 2 − 2 2 + 2 2 2 2 2 2 2 2 2 r + s + 2σ /δ2 δ2 d [4a /(δ2 d ) + (s − r ) ] δ2 d 4a /(δ2 d ) + (s − r ) 4a /(δ2 d ) + (r + s + 2σ /δ2 ) (s + r + 2σ /δ2 + 2af D /(δ2 d )) 64a4 (s − r + 2af D /(δ2 d )) 64a4 (s + r + 2σ / δ2 + 2af D /(δ2 d )) + 4 4 − 4 4 , 2 2 2 2 2 2 2 2 2 3 [4a /(δ2 d ) + (s + r + 2σ /δ2 ) ] δ2 d [4a /(δ2 d ) + (s − r ) ] δ2 d [4a2 /(δ22 d 2) + (s + r + 2σ /δ2 )2]3

(C-10)

 4 (r , s ) =

r + (δ2/ δ1 ) s + σ / δ1 + 1 1 (δ / δ ) s + σ / δ1 − r − 1 12a2 + 2 2 22 1 − 2 2 2 − 2 2 r + (δ2/ δ1 ) s + σ / δ1 + 1 δ1 d 4a /(δ1 d ) + ((δ2/ δ1 ) s + σ / δ1 − r − 1)2 4a /(δ1 d ) + (r + (δ2/ δ1 ) s + σ / δ1 + 1)2 ((δ2/ δ1 ) s + σ / δ1 − r − 1 + 2af D /(δ1 d )) 12a2 (r + (δ2/ δ1 ) s + σ / δ1 + 1 + 2af D /(δ1 d )) 64a4 + 2 2 + 4 4 2 2 2 2 2 [4a /(δ1 d ) + ((δ2/ δ1 ) s + σ / δ1 − r − 1) ] δ1 d [4a2 /(δ12 d 2) + (r + (δ2/ δ1 ) s + σ / δ1 + 1)2]2 δ1 d ((δ2/ δ1 ) s + σ / δ1 − r − 1 + 2af D /(δ1 d )) 64a4 (r + (δ2/ δ1 ) s + σ / δ1 + 1 + 2af D /(δ1 d )) − 4 4 2 2 2 2 3 [4a /(δ1 d ) + ((δ2/ δ1 ) s + σ / δ1 − r − 1) ] δ1 d [4a2 /(δ12 d 2) + (r + (δ2/ δ1 ) s + σ / δ1 + 1)2]3

(C-11)

1 (s, α ) =

(16δ14 (s

(4δ12 (s + 1)2 − 1) (f C − 2sin (2α ) − 2f C cos (2α ) − 4δ12 f C (s + 1)2 [2cos (2α ) − 5] + 16δ14 + 1) 4 − 8δ12 cos (2α )(s + 1)2 + 1)2

(s + 1) 4 [f C + 2sin (2α ) − 4f C cos (2α )] + 64δ16 f C (s + 1)6)

(C-12)

 2 (s, α ) =

(4(δ2 s + σ )2 − 1) [f D − 2sin (2α ) − 2f D cos (2α ) + 4f D (δ2 s + σ )2 (5 − 2cos (2α )) + 16(δ2 s + σ ) 4 (16(δ2 s + σ ) 4 − 8cos (2α )(δ2 s + σ )2 + 1)2 (f D + 2sin (2α ) − 4f D cos (2α )) + 64f D (δ2 s + σ )6]

(C-13)

Considering the Burgers vectors introduced in Eqs. (20) and (21) and using the related Gauss-Chebyshev integration formula [26], the singular integrals are calculated using the following discrete forms,

1 π

n

1

∫−1 f (r )(1 − r 2)1/2dr ≈ ∑ i=1

(1 − ri2 ) f (ri ) n+1

(C-14)

with

ri = cos(iπ /(n + 1)), i = 1, 2, ⋯, n ⎧ ⎨ ⎩ sm = cos(π (2m − 1)/(2(n + 1))), m = 1, 2, ⋯, n + 1

(C-15)

and

1 π

n

1

∫−1 f ()(1 − )−1/2 (1 + )1/2d  ≈ ∑ i=1

2(1 + i) f (i) 2n + 1

(C-16)

with

⎧ i = cos(π (2i − 1)/(2n + 1)), i = 1, 2, ⋯, n ⎨ ⎩ m = cos(π (2m)/(2n + 1)), m = 1, 2, ⋯, n + 1

(C-17)

Therefore, the discritized form of the integral Eqs. (C-6) and (C-7) are obtained as below, n

2

∑ ⎧ (1 − ri ) φ1 (ri ) ⎛

1 + 1 (ri, sm) ⎞ ⎫ + ⎠⎬ ⎝ ri − sm ⎭

⎜

⎟

n

∑ ⎧ 2(1 + i) φ2 (i) ⎛

1 + 2 (i, m ) ⎞ ⎫ ⎠⎬ ⎝ i − (δ1/ δ2 ) m + σ / δ2 − δ1/ δ2 ⎭ ⎜

⎨ ⎨ i = 1 ⎩ 2n + 1 ⎩ n+1 2 = 1 (sm, α ), m = 1, 2, ⋯, n + 1, 0 ≤ x ≤ c π

i=1

⎟

(C-18)

17

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al. n

∑ ⎧ 2(1 + i) φ2 (i) ⎛

1 + 3 (i, m ) ⎞ ⎫ + ⎠⎬ ⎝ i − m ⎭ ⎜

⎟

n

2

∑ ⎧ (1 − ri ) φ1 (ri ) ⎛

1 +  4 (ri, sm) ⎞ ⎫ ⎠⎬ ⎝ ri + 1 − (δ2/ δ1 ) sm − σ / δ1 ⎭

⎜

⎨ ⎨ i = 1 ⎩ 2n + 1 i=1 ⎩ n + 1 2 =  2 (m , α ), m = 1, 2, ⋯, n + 1, c < x ≤ e π

⎟

(C-19)

Eqs. (C-18) and (C-19) are solved within MATLAB to obtain φ1 and φ2 , resulting in the determination of Bx (ξ ) and

Bx∗ (ξ ) .

Appendix D. Solution of Eq. (31) We normalize Eq. (31) using the following changes of variables

⎧ x = cs, 0 ≤ |x| ≤ c ⎨ ⎩ ξ = cr , 0 ≤ |ξ| ≤ c

(D-1)

x

⎧ d = δ2 s + σ , c < |x| ≤ e ⎨ ξ = δ2 r + σ , c < |ξ| ≤ e ⎩d

(D-2)

δ2 and σ have been previously determined by Eq. (C-5). Applying Eqs. (D-1) and (D-2) into the Eq. (31), the following relations are obtained, 1 C⎧ ⎨ π ⎩

1

⎧1 ⎨ ⎩π 1 C⎧ ⎨ ⎩π

D

2μ − η) − η) C f ∫ K11 (r , s ) Bx (r ) dr ⎫ + D ∫−1 Br x−(rs) dr + c (12πa ∫−1 K21 (r , s) Bx (r ) dr + c (12πa −1 ⎬ (κ + 1) 1

1

∫−1

Bx∗ (r ) 1 dr + r + σ / δ2 − (c /(δ2 d )) s π

1

1

∫−1 5 (r ,

⎭ 2P = −⎛ ⎞  3 (s, α ), 0 ≤ |x| ≤ c, ⎬ ⎝ πd ⎠ ⎭

s ) Bx∗ (r ) dr ⎫

(D-3)

2μD

− η) − η) D f ∫ K12 (r , s ) Bx (r ) dr ⎫ + D ∫−1 r − ((δ2 dB)/x (cr))s − σd/c dr + c (12πa ∫−1 K22 (r , s) Bx (r ) dr + c (12πa −1 ⎬ (κ + 1) 1

1

1

⎭

⎧1 ⎨ ⎩π

1

∫−1

Bx∗ (r ) 1 dr + r−s π

2P Bx∗ (r ) 3 (r , s ) dr ⎫ = −⎛ ⎞  2 (s, α ), c < |x| ≤ e −1 ⎬ ⎝ πd ⎠ ⎭

∫

1

(D-4)

in which, ∞

K11 (r , s ) =

∫0

K12 (r , s ) =

∫0

K21 (r , s ) =

∫0

K22 (r , s ) =

∫0

∞

∞

∞

N12 c (r − s ) t exp (−t )cos dt D a

(D-5)

N12 (cr − δ2 ds − σd) t exp (−t )cos dt D a

(D-6)

N22 c (r − s ) t exp (−t )sin dt D a

(D-7)

N22 (cr − δ2 ds − σd) t exp (−t )sin dt D a

(D-8)

5 (r , s ) =

1 (c /(δ d )) s − r − σ / δ2 r + σ / δ2 + (c /(δ2 d )) s 12a2 + 2 2 2 2 − 2 2 2 − 2 2 r + σ / δ2 + (c /(δ2 d )) s 4a /(δ2 d ) + ((c /(δ2 d )) s − r − σ / δ2 )2 4a /(δ2 d ) + (r + σ / δ2 + (c /(δ2 d )) s )2 δ2 d ((c /(δ2 d )) s − r − σ / δ2 + 2af C /(δ2 d )) 12a2 (r + σ / δ2 + (c /(δ2 d )) s + 2af C /(δ2 d )) 64a4 + 2 2 + 4 4 [4a2 /(δ22 d 2) + ((c /(δ2 d )) s − r − σ / δ2 )2]2 δ2 d [4a2 /(δ22 d 2) + (r + σ / δ2 + (c /(δ2 d )) s )2]2 δ2 d ((c /(δ2 d )) s − r − σ / δ2 + 2af C /(δ2 d )) 64a4 (r + σ / δ2 + (c /(δ2 d )) s + 2af C /(δ2 d )) − 4 4 [4a2 /(δ22 d 2) + ((c /(δ2 d )) s − r − σ / δ2 )2]3 δ2 d [4a2 /(δ22 d 2) + (r + σ / δ2 + (c /(δ2 d )) s )2]3

(D-9)

 3 (s, α ) =

(4(c 2/ d 2) s 2 − 1) (f C − 2sin (2α ) − 2f C cos (2α ) + 4f C (c 2/d 2) s 2 (5 − 2cos(2α )) − 8cos (2α )(c 2/ d 2) s 2 + 1)2

(16(c 4/ d4 ) s 4

+ 16(c 4/ d4 ) s 4 (f C + 2sin(2α ) − 4f C cos (2α )) + 64f C (c 6/d6) s 6)

(D-10)

We employ the Legendre-Gauss quadrature algorithm to calculate the infinite integrals given in Eqs. (D-5)–(D-8). 3 (r , s ) and  2 (s, α ) have been previously presented in Appendix C. Considering the Burgers vectors of Bx (ξ ) and Bx∗ (ξ ) given, respectively, in Eqs. (32) and (33), Eqs. (D-3) and (D-4) can be rewritten as

18

Engineering Fracture Mechanics xxx (xxxx) xxxx

M. Tabatabaei, et al. n

C∑ i=1

(1 − ri2 ) 2μD c (1 − η) C c (1 − η) 1 + f K11 (ri, sm) ⎫ + D φ (ri ) ⎧ K21 (ri, sm) + ⎬ ( 2 2 κ − + 1) a n+1 1 ⎨ r s a m ⎭ ⎩ i

n

∑ 2(1 + i) φ2 i=1

2n + 1

1 2 (i) ⎧ + 2 (i, m ) ⎫ = −⎛ ⎞ 1 (sm, α ), 0 ≤ |x| ≤ c, ⎬ ⎨ ⎝π ⎠ ⎭ ⎩ i + σ / δ2 − (c /(δ2 d )) m n

C∑ i=1

(D-11)

ri2 )

(1 − 2μD c (1 − η) D c (1 − η) 1 f K12 (ri, sm) ⎫ + D + φ (ri ) ⎧ K22 (ri, sm) + ⎬ ( 2 (( )/ ) / 2 κ a − − + 1) n+1 1 ⎨ r δ d c s σd c a i m 2 ⎭ ⎩

1 2 (i) ⎧ + 3 (i, m ) ⎫ = −⎛ ⎞  2 (m , α ), c < |x| ≤ e ⎬ ⎨ ⎝π ⎠ ⎭ ⎩ i − m Eqs. (D-11) and (D-12) are solved within MATLAB to determine Burgers vectors of Bx (ξ ) and

n

∑ 2(1 + i) φ2 i=1

2n + 1

(D-12)

Bx∗ (ξ ) .

References [1] Nelson EB, Guillot D. Well cementing. 2nd ed. Texas: Schlumberger; 2006. [2] Wojtanowicz AK, Nishikawa S, Xu R. Diagnosis and remediation of sustained casing pressure in wells. Final report submitted to US Department of Interior Minerals Management Service; 2001. [3] Xu R, Wojtanowicz AK. Pressure buildup test analysis in wells with sustained casing pressure. J Nat Gas Sci Engng 2017;38:608–20. https://doi.org/10.1016/j. jngse.2016.12.033. [4] Williams RH, Khatri DK, Keese RF, Le Roy-Delage S, Roye JM, Leach DLR, et al., Expanding cement system (FECS) successfully provides zonal isolation across marcellus shale gas trends. Can Unconv Resour Conf 2011; SPE-149440-MS. doi: 10.2118/149440-MS. [5] Ingraffea AR, Wells MT, Santoro RL, Shonkoff SBC. Assessment and risk analysis of casing and cement impairment in oil and gas wells in Pennsylvania, 2000–2012. PNAS 2014;111:10955–60. https://doi.org/10.1073/pnas.1323422111. [6] Carter LG, Evans GW. A study of cement-pipe bonding. J Petrol Technol 1964;16:157–61. https://doi.org/10.2118/764-PA. [7] Kiran R, Teodoriu C, Dadmohammadi Y, Nygaard R, Wood D, Mokhtari M, et al. Identification and evaluation of well integrity and causes of failure of well integrity barriers (A review). J Nat Gas Sci Engng 2017;45:511–26. [8] Wang W, Dahi Taleghani A. Three-dimensional analysis of cement sheath integrity around Wellbores. J Petrol Sci Engng 2014;121:38–51. https://doi.org/10. 1016/j.petrol.2014.05.024. [9] Tabatabaei M, Dahi Taleghani A. Randomly distributed interfacial arc cracks within the inclusion-inhomogeneity-matrix system. Meccanica 2017;52:1123–42. https://doi.org/10.1007/s11012-016-0442-y. [10] Dahi Taleghani A, Klimenko D. An analytical solution for microannulus cracks developed around a wellbore. J Energy Resour Technol 2015;137:062901https:// doi.org/10.1115/1.4030627. [11] Ayatollahi MR, Sistaninia M. Mode II fracture study of rocks using Brazilian disk specimens. Int J Rock Mech Min Sci 2011;48:819–26. [12] Wang J-S, Suo Z. Experimental determination of interfacial toughness curves using Brazil-Nut-Sandwiches. Acta Metall Mater 1990;38:1279–90. https://doi.org/ 10.1016/0956-7151(90)90200-Z. [13] Dorogoy A, Banks-Sills L. Effect of crack face contact and friction on Brazilian disk specimens—a finite difference solution. Engng Fract Mech 2005;72:2758–73. https://doi.org/10.1016/j.engfracmech.2005.05.005. [14] Atkinson C, Smelser RE, Sanchez J. Combined mode fracture via the cracked Brazilian disk test. Int J Fract 1982;18:279–91. https://doi.org/10.1007/ BF00015688. [15] Dong S, Wang Y, Xia Y. Stress intensity factors for central cracked circular disk subjected to compression. Engng Fract Mech 2004;7–8:1135–48. https://doi.org/ 10.1016/S0013-7944(03)00120-6. [16] Zhang Z, Scherer GW, Prud’homme RK. Adhesion and bonding between steel pipe and cement/spacer/mud system. FraMCoS-9 2016. https://doi.org/10.21012/ FC9.313. [17] Comninou M, Schmueser D. Frictional slip between a layer and a substrate caused by a normal load. Int J Engng Sci 1980;18:131–7. https://doi.org/10.1016/ 0020-7225(80)90012-9. [18] Dundurs J, Mura T. Interaction between an edge dislocation and circular inclusion. J Mech Phys Solids 1964;12:177–89. https://doi.org/10.1016/00225096(64)90017-1. [19] Dundurs J, Sendeckyj GP. Behavior of an edge dislocation near a bimetallic interface. J Appl Phys 1965;36:3353. https://doi.org/10.1063/1.1702981. [20] Dundurs J. Mathematical theory of dislocations. In: Mura T, editor. Elastic interaction of dislocations with inhomogeneities. New York: American Society of Mechanical Engineers; 1969. p. 70–113. [21] Comninou M, Dundurs J. Partial closure of cracks at the interface between a layer and a half space. Eng Frac Mech 1983;18:315–23. https://doi.org/10.1016/ 0013-7944(83)90142-X. [22] Erdogan F, Gupta GD. On the numerical solution of singular integral equations. Quart Appl Math 1972;29:525–34. https://doi.org/10.1090/qam/408277. [23] Erdogan F, Gupta GD, Cook TS. Numerical solution of singular integral equations. In: Sih GC, editor. Mechanics of fracture: methods of analysis and solutions of crack problems. Netherland: Noordhoff Int. Publ; 1973. p. 368–425. [24] Boussinesq J. Application des Potentiels á Létude de léquilibre et due Mouvement des Solides Elastique. Paris: Gauthier Villars; 1885. [25] Timoshenko SP, Goodier JN. Theory of elasticity. New York: McGraw-Hill; 1970. [26] Abramowitz M, Stegum IA. Handbook of mathematical functions: with formulas, graphs, and mathematical tables, Washington D.C.: National Bureau of Standards Applied Mathematics Series 55; 1964.

19