# Computation of Energy Release Rates for a Nearly Circular Crack

Authors: Nik Mohd Asri Nik Long; Lee Feng Koo; Zainidin K. Eshkuvatov
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

This paper deals with a nearly circular crack,Ω in the plane elasticity. The
problem of finding the resulting shear stress can be formulated as a hypersingular
integral equation over a considered domain, Ω and it is then transformed into a
similar equation over a circular region, D, using conformal mapping.  Appropriate
collocation points are chosen on the region D to reduce the hypersingular integral
equation into a system of linear equations with (2N+1)(N+1) unknown coefficients,
which will later be used in the determination of energy release rate. Numerical
results for energy release rate are compared with the existing asymptotic solution
and are displayed graphically.


---

## Body

## 1. Introduction

The determination of energy release rate, a measurement of energy necessary for crack initiation in fracture mechanics, has stirred a huge interest among researchers, and different approaches have been applied. Williams and Isherwood [1] proposed an approximate method in terms of a mean stress to approximate the strain-energy release rates of finite plates. Sih [2] proposed the energy density theory as an alternative approach for fracture prediction. Hayashi and Nemat-Nasser [3] modelled the kink as a continuous distribution of infinitesimal edge dislocations to obtain the energy release rate at the onset of kinking of a straight crack in an infinite elastic medium subjected to a predominantly Mode I loading. Further, a similar method to [3] has also been adopted by Hayashi and Nemat-Nasser [4] to obtain the energy release rate for a kinked from a straight crack under combined loading based on the maximum energy release rate criterion. Gao and Rice [5] extended Rice's work [6] in finding the energy release rate for a plane crack with a slightly curved front subject to shear loading. While, Gao and Rice [7] and Gao [8] considered a penny-shaped crack as a reference crack in solving the energy release rate for a nearly circular crack subject to normal and shear loads. Jih and Sun [9] employed the finite element method based on crack-closure integral in calculating the strain energy release rate elastostatic and elastodynamic crack problems in finite bodies whereas Dattaguru  et al. [10] adopted the finite element analysis and modified crack closure integral technique in evaluating the strain energy release rate. Poon and Ruiz [11] applied the hybrid experimental-numerical method for determining the strain energy release rate. Wahab and de Roeck [12] evaluated the strain energy release rate from three-dimensional finite element analysis with square-root stress singularity using different displacement and stress fields based on the Irwin's crack closure integral method [13]. Guo  et al. [14] used the extrapolation approach in order to avoid the disadvantages of self-inconsistency in the point-by-point closed method to determine the energy release rate of complex cracks. Xie  et al. [15] applied the virtual crack closure technique in conjunction with finite element analysis for the computation of energy release rate subject to kinked crack, while interface element based on similar approach also adopted by Xie and Biggers [16] in calculating the strain energy release rate for stationary cracks subjected to the dynamic loading.In this paper, we focus our work on obtaining the numerical results for energy release rate for a nearly circular crack via the solution of hypersingular integral equation and compare our computational results with Gao's [8].


## 2. Formulation of the Problem

Consider the infinite isotropic elastic body containing a flat circular crack,Ω,  as in Figure 1, located on the Cartesian coordinate (x,y,x3) with origin O, and Ω lies in the plane x3=0. Let the radius of the crack, Ω be a and Ω={(r,θ):0≤r<a,-π≤θ<π}.Figure 1
Stresses acting on a circular crack.If the equal and opposite shear stresses in thex and y directions, q1(x,y) and q2(x,y), respectively, are applied to the crack plane, and it is assumed that the x3 direction is traction free, then in the view of shear load, the entire plane, must free from the normal stress, that is
(2.1)τ33(x,y,x3)=0forx3=0,
and the stress field can be found by considering the above problem subjected to the following mixed boundary condition on its surface, x3=0:
(2.2)τ13(x,y,x3)=μ1-νq1(x,y),(x,y)∈Ω,τ23(x,y,x3)=μ1-νq2(x,y),(x,y)∈Ω,u1(x,y,x3)=u2(x,y,x3)=0,(x,y)∈Γ∖Ω,
where τij is stress tensor, μ is shear modulus, ν is denoted as Poisson's ratio, and Γ is the entire x3=0. Also, the problem satisfies the regularity conditions at infinity
(2.3)ui(x,y,x3)=O(1R),τij(x,y,x3)=O(1R),i,j=1,2,3,R→∞,
where R is the distance
(2.4)R=(x-x0)2+(y-y0)2,(x0,y0)∈Ω.
Martin [17] showed that the problem of finding the resultant force with condition (2.2) can be formulated as a hypersingular integral equation
(2.5)18π⨘Ω(2-ν)w(x,y)+3νe2jΘw(x,y)¯R3dΩ=q(x0,y0),(x0,y0)∈Ω,
where w(x,y)=[u1(x,y)]+j[u2(x,y)] is the unknown crack opening displacement, q(x0,y0)=q1(x0,y0)+jq2(x0,y0), j2=-1, the w(x,y)¯=[u1(x,y)]-j[u2(x,y)], and the angle Θ is defined by
(2.6)x-x0=Rcos⁡Θ,y-y0=RsinΘ.
The cross on the integral means the hypersingular, and it must be interpreted as a Hadamard finite part integral [18, 19]. Equation (2.5) is to be solved subject to w=0 on ∂Ω where ∂Ω is boundary of Ω. For the constant shear stress in x direction, we have τ23=0 and [u2(x,y)]=0, hence, (2.5) becomes
(2.7)18π⨘Ω2-ν+3νe2jΘR3w(x,y)dΩ=q(x0,y0),(x0,y0)∈Ω.Polar coordinates(r,θ) and (r0,θ0) are chosen so that the loadings q(x,y) and q(x0,y0) can be written as a Fourier series
(2.8)q(x,y)=∑n=-∞∞qn(ra)ejnθ,q(x0,y0)=∑n=-∞∞qn(r0a0)ejnθ0,
where the Fourier components qn are j-complex. The j-complex crack opening displacement, w(x,y) and w(x0,y0), have similar expressions
(2.9)w(x,y)=∑n=-∞∞wn(ra)ejnθ,w(x0,y0)=∑n=-∞∞wn(r0a0)ejnθ0.
Without loss of generality, we consider a=1. Using Guidera and Lardner [20], the dimensionless function qn and wn can be expressed as
(2.10)qn(r)=r|n|∑k=0∞QknΓ(|n|+12)Γ(k+32)(|n|+k)!1-r2C2k+1|n|+12(1-r2),wn(r)=r|n|∑k=0∞WknΓ(|n|+1/2)k!Γ(|n|+k+3/2)C2k+1|n|+1/2(1-r2),
where the j-complex coefficients Qkn are known, Wkn are unknown, and Cmλ(x) is an orthogonal Gegenbauer polynomial of degree m and index λ, which is defined recursively by [21]
(2.11)(m+2)Cm+2λ(x)=2(m+λ+1)xCm+1λ(x)-(2λ+m)Cmλ(x),
with the initial values C0λ(x)=1 and C1λ(x)=2λx. For a constant shear loading, q(x,y)=-τ, the solution for a circular crack is obtainable.


## 3. Nearly Circular Crack

LetΩ be an arbitrary shaped crack of smooth boundary with respect to origin O, such that Ω is defined as
(3.1)Ω={(r⋅θ):0≤r<ρ(θ),-π≤θ<π},
where the boundary of Ω, ∂Ω is given by r=ρ(θ). Let ζ=ξ+iη=seiφ with |ζ|<1 such that the unit disc is
(3.2)D≡{(s,φ):0≤s<1,-π≤φ<π}.
By the properties of Reimann mapping theorem [22], a circular disc D is mapped conformally onto Ω using z=af(ζ). This approach works for a general smooth star-shaped domain, Ω. For a particular application, let  f be an analytic function, simply connected in the domain Ω, |f′(ζ)| is nonzero and bounded for all |ζ|<1,
(3.3)f(ζ)=ζ+cg(ζ)withg(ζ)=ζm+1,
which maps a unit circle, D in the ζ-plane into a nearly circular domain Ω in the z-plane where c is a real parameter and r=ρ(θ) is the boundary of Ω. This domain has a smooth, regular boundary for 0≤(m+1)|c|<1.  As (m+1)|c|→1 one or more cusps develop; see Figure 2 with various choices of c.Figure 2
The domainΩ for f(ζ)=ζ+cζm+1 at different choices of c, m=2.Let(3.4)z-z0=a(f(ζ)-f(ζ0))=ReiΘ,
and define S and Φ as (3.5)ζ-ζ0=SeiΦ,dΩ=dxdy=a2|f′(ζ)|2dξdη=a2|f′(ζ)|2sdsdφ,
where x=au(ξ,η) and y=av(ξ,η) so that f=u+iv. Next, we define δ and δ0 as
(3.6)f′(ζ)=|f′(ζ)|eiδ,f′(ζ0)=|f′(ζ0)|eiδ0.
Set
(3.7)w(x(ζ),y(ζ))=a|f′(ζ)|-1/2ejδW(ξ,η),(3.8)q(x(ζ0),y(ζ0))=a|f′(ζ0)|-3/2ejδ0Q(ξ0,η0).
Substituting (3.5), (3.6), (3.7), and (3.8) into (2.7) gives(3.9)2-ν+3νe2jΘ8π⨘DW(ξ,η)S3dξdη+2-ν8π⨍DW(ξ,η)K(1)(ζ,ζ0)dξdη+3ν8π∫DW(ξ,η)K(2)(ζ,ζ0)dξdη=Q(ξ0,η0),(ξ0,η0)∈D,
where the kernel K(1)(ζ,ζ0) and K(2)(ζ,ζ0) are [17](3.10)K(1)(ζ,ζ0)=|f′(ζ)|3/2|f′(ζ0)|3/2|f(ζ)-f(ζ0)|3ej(δ-δ0)-1|ζ-ζ0|3,(3.11)K(2)(ζ,ζ0)=|f′(ζ)|3/2|f′(ζ0)|3/2|f(ζ)-f(ζ0)|3ej(2Θ-δ-δ0)-1|ζ-ζ0|3e2jΦ.
This hypersingular integral equation over a circular disc D is to be solved subject to W=0 on s=1, and the K(1)(ζ,ζ0) is a Cauchy-type singular kernel with order S-2, and the kernel K(2)(ζ,ζ0) is weakly singular with O(S-1), as ζ→ζ0 (see the appendix).We are going to solve (3.9) numerically. Write W(ξ,η) as a finite sum
(3.12)W(ξ,η)=∑n,kWknAkn(s,φ),
where Akn(s,φ) is defined by
(3.13)Akn(s,φ)=s|n|C2k+1|n|+1/2(1-s2)ejnφ,∑n,k=∑n=-N1N1∑k=0N2,N1,N2∈Z.
Introduce(3.14)Lhm(s,φ)=s|m|C2h+1|m|+1/2(1-s2)cos⁡⁡mφ,
where m,h∈ℤ. The relationship between these two functions, Akn(s,φ), and Lhm(s,φ) can be expressed as
(3.15)∫ΩAkn(s,φ)Lhm(s,φ)sdsdφ1-s2=Bknδkhδmn,
where δij is Kronecker delta and
(3.16)Bkn={2π4k+3,n=0,π2Γ(2k+2n+2)22n+1(2k+n+3/2)(2k+1)![Γ(n+1/2)]2,n≠0.
Both functions Akn(s,φ) and Lhm(s,φ) have square-root zeros at s=1.Krenk [23] showed that
(3.17)14π⨘ΩAkn(s,φ)R3dΩ=-EknAkn(s0,φ0)1-s02,
where
(3.18)Ekn=Γ(|n|+k+3/2)Γ(k+3/2)(|n|+k)!k!.Substituting (3.17) and (3.12) into (3.9) yields
(3.19)∑n,kFkn(s0,φ0)Wkn=Q(ξ0(s0,φ0),η0(s0,φ0)),
where
(3.20)Fkn(s0,φ0)=-Ekn(2-ν+3νe2jΘ)Akn(s0,φ0)21-s02+2-ν8π∫DAkn(s,φ)K(1)(ζ,ζ0)dξdη+3ν8π∫DAkn(s,φ)K(2)(ζ,ζ0)dξdη;0≤s≤1,0≤φ<2π.Next, define(3.21)Wkn=-W̃knG2k+1|n|+1/2EknBkn,
where G2k+1|n|+1/2=(2n+2k+1)!/(2k+1)!(2n)!. Multiply (3.19) by Lhm(s0,φ0), integrate over D and using (3.15), (3.19) becomes
(3.22)∑n,kW̃kn(-2-ν+3νe2jΘ2δhkδ|m||n|+Shkmn)=Qhm,-N1≤m≤N1,0≤h≤N2,
where
(3.23)Shkmn=18πEknBknEhmBhmThkmn,Thkmn=∫DLhm(ζ0)∫DAkn(ζ)H(ζ,ζ0)dζdζ0,Qhm=1EhmBhm∫DLhm(ζ0)Q(ζ0)dζ0,H(ζ,ζ0)=(2-ν)K(1)(ζ,ζ0)+3νK(2)(ζ,ζ0).In (3.22), we have used the following notation:  ζ0=ζ0(s0,φ0), dζ0=s0ds0dφ0, and Q(ζ0)=Q(ξ0,η0)=Q(s0cosφ0,s0sinφ0).In evaluating the multiple integrals in (3.22), we have used the Gaussian quadrature and trapezoidal formulas for the radial and angular directions, with the choice of collocation points (s,φ) and (s0,φ0) defined as follows:
(3.24)si=π4+π4∑i=1M1W(i),s0i=π4+π4∑i=1M1W0(i),φj=∑j=1M2jπM2,φ0j=∑j=1M2(j+0.5)πM2,
where W(i) and W0(i) are abscissas for si and s0i, respectively, M1 and M2 is the number of collocation points in radial and angular directions, respectively. This effort leads to the (2N1+1)(N2+1)×(2N1+1)(N2+1) system of linear equations
(3.25)AW̃=b̃,
where A is a square matrix, and W̃ and b̃ are vectors, W̃ to be determined.


## 4. Energy Release Rate

The energy release rate (measured inJM-2), G(φ) by Irwin's relation subject to shear load is defined as [7, 8]
(4.1)G(φ)=(1-ν2)E[KII(φ)]2+(1+ν)E[KIII(φ)]2,
where E, Young's modulus, a measurement of the stiffness of an isotropic elastic material and the relationship of E, ν and μ, is
(4.2)ν=E2μ-1,
and KII(φ) and KIII(φ), the sliding and tearing mode stress intensity factor, respectively, are defined as [5, 7, 8]
(4.3)Kj(φ)=lim⁡r→aVj2πa-rw(x,y),j=II,III,
where Vj are constants.Leta(φ)=|f(eiφ)|, r=|f(seiφ)|, and as s close to 1, (4.3) leads to
(4.4)Kj(φ)=lim⁡s→1-Vj2π(1-s)|f′(eiφ)|w(x,y),j=II,III.
Therefore, substituting (3.7) into (4.4) and simplifying gives
(4.5)Kj(φ)=Vj{|f′(eiφ)|-1∑n,kW̃knEknBknYkn(φ)},j=II,III,
where Ykn(φ)=D2k+1|n|+1/2(0)cos(nφ), and C2k+1|n|+1/2(1-s2)=1-s2D2k+1|n|+1/2(1-s2), where Dmλ(x) is defined recursively by
(4.6)mDmλ(x)=2(m+λ-1)xDm-1λ(x)-(m+2λ-2)Dm-2λ(x),m=2,3,…,
with D0λ(x)=2λ and D1λ(x)=2λx.Table1 shows that our numerical scheme converges rapidly at a different point of the crack with only a small value of N=N1=N2 are used.Table 1
Numerical convergence for the energy release rate,G(φ) for f(ζ)=ζ+cζ3 when c=0.1.
NG(0.00)G(π/4)G(π/2)G(3π/4)G(π)07.8676E-109.0123E-101.6067E-099.0123E-107.8676E-1017.2724E-108.9392E-101.3159E-098.9392E-107.2724E-1029.2668E-077.4652E-101.5649E-097.4652E-109.2668E-0730.0000E+000.0000E+006.3517E-100.0000E+000.0000E+0041.1859E-057.4041E-194.6709E-097.4041E-191.1859E-0553.1429E-038.8211E-049.2528E-068.8211E-043.1429E-0363.0421E-038.7908E-049.5791E-048.7908E-043.0421E-0371.5794E-038.4308E-049.2945E-048.4308E-041.5721E-0389.7557E-041.1903E-039.5001E-041.1903E-039.7557E-0499.7557E-041.1903E-039.5001E-041.1903E-039.7557E-04109.7557E-041.1903E-039.5001E-041.1903E-039.7557E-04Figures3, 4, 5, and 6 show the variations of G against φ for c=0.001, c=0.01, c=0.10, and c=0.30, respectively. It can be seen that the energy release rate has local extremal values when the crack front is at cos(φ)=±1 or sin(φ)=±1. Similar behavior can be observed for the solution of G(φ) for a different c and ν at c=0.1, displayed in Figures 7 and 8. Our results agree with those obtained asymptotically by Gao [8], with the maximum differences for m=2 are 3.6066×10-6, 4.7064×10-5, 5.3503×10-5, and 9.0000×10-5 for c=0.001, c=0.01, c=0.10, and c=0.30, respectively.Figure 3
The energy release rate,G(φ) for f(ζ)=ζ+0.001ζ3.Figure 4
The energy release rate,G(φ) for f(ζ)=ζ+0.01ζ3.Figure 5
The energy release rate,G(φ) for f(ζ)=ζ+0.1ζ3.Figure 6
The energy release rate,G(φ) for f(ζ)=ζ+0.3ζ3.Figure 7
The energy release rate,G(φ) for f(ζ)=ζ+cζ3 at various c.Figure 8
The energy release rate,G(φ) for f(ζ)=ζ+0.1ζ3 at various choices of ν with μ=1.


## 5. Conclusion

In this paper, the hypersingular integral equation over a nearly circular crack is formulated. Then, using the conformal mapping, the equation is transformed into hypersingular integral equation over a circular crack, which enable us to use the formula obtained by Krenk [23]. By choosing the appropriate collocation points, this equation is reduced into a system of linear equations and solved for the unknown coefficients. The energy release rate for the mentioned crack subject to shear load is presented graphically. Our computational results seem to agree with the asymptotic solution obtained by Gao [8].



---

*Source: 741075-2011-03-08.xml*