# Numerical Solution for an Epicycloid Crack

Authors: Nik Mohd Asri Nik Long; Koo Lee Feng; Wong Tze Jin; Z. K. Eshkuvatov
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

A flat crack,Ω, is lying in a three-dimensional homogenous isotropic elastic solid subjected to shear loading. A mathematical formulation is developed based on the mixed boundary values for Ω such that the problem of finding the resulting force can be written in the form of hypersingular integral equation. Employing conformal mapping, the integral equation is transformed to a similar equation over a circular region, D. By making a suitable representation of hypersingular integral equation, the problem is reduced to solve a system of
linear equations. Numerical solution for the shear stress intensity factors, maximum stress intensity, and strain energy release rate is obtained. Our results give an excellent agreement to the existing asymptotic solutions.


---

## Body

## 1. Introduction

Crack problems play an important role in engineering application due to the fact that the presence of cracks may compromise the strength and toughness of structures. Hence, great efforts [1–10] have been made in solving the crack problems and it has been widely investigated since the pioneer work by Sneddon [11] for a penny-shaped crack. Integral transform method to the solution of a Fredholm integral equation of second kind and numerical approach was implemented by Kassir [12, 13] in solving the rectangular crack problem, while the classic collocation and Galerkin methods were applied by Ioakimidis [14] for solving the plane crack problem subjected to normal load, whereas a perturbation analysis and the complex potential method [15] were performed by Cotterell and Rice [16] to obtain the stress intensity factors for the curved and kinked crack subject to arbitrary tractions in an explicit and simple form. Recently, Wu [17] proposed the dual boundary element method to solve the antiplane crack problem, whilst Georgiadis and Gourgiotis [18] advocated distributed dislocation technique in solving crack problems within Cosserat elasticity with constrained rotations. Motivated by the work of Lazzarin and Zappalorto [19], Lazzarin et al. [20] investigated the stress fields close to a rectangular hole in a plate of finite thickness.Ioakimidis [21] introduced the concept of finite-part integrals and derived the hypersingular integral equation for a flat crack subjected to tensile pressure, where the unknown function is the crack displacement discontinuity while the right-hand terms are the applied tractions on the crack faces. This equation can be numerically solved effectively by using the Gaussian quadrature rules for finite-part integral. Hence, this concept had been advocated widely for the solution of crack problem and some of them can be found in [22–26].In this paper, the epicycloid crack problem is formulated into solving the hypersingular integral equation numerically for finding the stress intensity factors, maximum stress intensity and energy release rate for the crack subject to shear loading. Our computational results agree with the existing asymptotic solution.


## 2. Statement of Problem and Basic Equations

Consider an arbitrary shaped crack,Ω, embedded in three-dimensional unbounded isotropic elastic body, Γ. Let the Cartesian coordinate (x,y,z) with origin O and Ω lie in the plane z=0. Assume that O is a point in Ω and the body force is absent. Let the crack edges be deformed by the application of equal and opposite constant shear stresses in the x and y directions, qx(x,y) and qy(x,y), and it is assumed that the z direction is traction-free; see Figure 1. Hence, in view of the shear load, the entire plane must be free from the normal stress; that is,
(1)τzz=0forz=0.
Thus, the stress field can be found by considering the half-space, z≥0, subject to the following mixed boundary condition on its surface z=0:
(2)τxz=μ1-νqx(x,y),(x,y)∈Ω,τyz=μ1-νqy(x,y),(x,y)∈Ω,ux(x,y,z)=uy(x,y,z)=0,(x,y)∈Γ∖Ω,
where τzz, τxz, and τyz denote the stress tensor, μ is shear modulus, and ν is Poisson’s ratio and the usual regularity requirements at the location away from the crack region,
(3)ui(x,y,z)=O(1R),τij(x,y,z)=O(1R),
where i,j=x,y,z, R→∞, R=(x-x0)2+(y-y0)2, ui are displacement vectors and τij is given by
(4)τij=cijkl∂ul∂xk;k,l=x,y,z,
where cijkl are the material moduli defined by
(5)cijkl=λδijδkl+μ(δikδjl+δilδjk),
where ν=λ/2(λ+μ). The δij is the Kronecker delta, defined as 1 if i=j and 0 if i≠j. And (5) must satisfy Hooke’s law linear elasticity symmetry conditions:
(6)cijkl=cjikl=cklij=cijlk.
The displacement vector, ui, is represented by Somigliana formula [27, 28]:
(7)um(x0,y0)=∫Ω‍[ui(x,y)]∑ijmF‍((x,y);(x0,y0))njdΩ,
where the component of Green function is
(8)∑ijmF‍((x,y);(x0,y0))=cijkl∂∂xkGlmF((x,y);(x0,y0)),(9)8πμGijF((x,y);(x0,y0))=18πμ(2Rδij-12(1-v)∂2R∂xi∂xj),
and [ui(x,y)] is the displacement discontinuity in ui across the crack,
(10)[ui(x,y)]=lim⁡(x0,y0)→(x,y)∈Ω+ui(x0,y0)-lim⁡(x0,y0)→(x,y)∈Ω-ui(x0,y0),
where nj is the unit normal vector, which is assumed to point into Ω. Equation (9) is known as Kelvin’s point-load solution. Substitute (5) and (9) into (7) and perform the integration by parts with respect to x and y, yielding a system of Cauchy principle-value integral equations [29]:
(11)qx(x0,y0)=-14π⨍Ω{α∂∂x(1R)+β∂∂y(1R)}dΩ,qy(x0,y0)=-14π⨍Ω{α∂∂y(1R)-β∂∂x(1R)}dΩ
for (x0,y0)∈Ω, where
(12)α=∂[ux]∂x+∂[uy]∂y,β=(1-ν)(∂[ux]∂y-∂[uy]∂x).
The resulting boundary terms which involved [ux(x,y)] and [uy(x,y)] are evaluated at the crack edge, assuming that there are no tractions applied onto the boundary. Consequently, [ux(x,y)] and [uy(x,y)] are zero; that is, these equations are to be solved subject to
(13)[ux(x,y)]=0,[uy(x,y)]=0for(x,y)∈∂Ω,
where ∂Ω is the boundary of Ω. Integrating (11) by parts and using condition (13) and making use of the relationship between Cauchy principle-value integral and hypersingular integral equations [30, 31],
(14)ddx⨍abf(t)(x-t)dt=-⨘abf(t)(x-t)2dt,t∈(a,b),
yield [29, 32]
(15)qx(x0,y0)=18π⨘Ω‍(2-ν+3νcos⁡2Θ)[ux]+3νsin2Θ[uy]8πR3dΩ,(16)qy(x0,y0)=18π⨘Ω‍3νsin2Θ[ux]+(2-ν-3νcos⁡2Θ)[uy]R3dΩ,
and the angle Θ is defined by x-x0=Rcos⁡Θ and y-y0=RsinΘ. The cross on the integral of (15) and (16) means the hypersingular, and it must be interpreted as a Hadamard finite part integral [31, 33, 34]. Multiplying (16) with complex j and adding to (15) lead to
(17)q(x0,y0)=18π⨘Ω‍(2-ν)w(x,y)+3νe2jΘw(x,y)¯R3dΩ,mmnmmmmmmmmmmmmmmmmmmm(x0,y0)∈Ω,
where q(x0,y0)=qx(x0,y0)+jqy(x0,y0), w(x,y)=[ux]+j[uy] is the unknown crack opening displacement, and the bar denotes the conjugation of  w(x,y)¯=[ux]-j[uy] and j2=-1. Equation (17) is to be solved subject to w=0 on ∂Ω and can be used for general crack problems under shear loading which is equivalent to those equations obtained in [27, 35]. Suppose the constant shear stress is applied on opposite crack surfaces at x direction, and then the general solution of (17) can be reduced into a single hypersingular integral equation:
(18)q(x0,y0)=18π⨘Ω‍2-ν+3νe2jΘR3w(x,y)dΩ,mmmmmmmimnmnm(x0,y0)∈Ω.Figure 1
Stresses acting on a plane.


## 3. Conformal Mapping and Epicycloid Cracks

Suppose thatΩ is a penny-shaped crack, with radius a so that the crack occupies the region
(19)Ω={(r,θ):0≤r<a,-π≤θ<π},
where r and θ are polar coordinates, x=rcos⁡θ, and y=rsinθ.Now, letΩ be a simply connected domain in the z-plane defined as
(20)Ω={(r·θ):0≤r<ρ(θ),-π≤θ<π}
whose boundary has the polar equation r=1+cρ(θ), where ρ(θ) is bounded and piecewise continuous and c is a small positive parameter. Define ζ=seiϕ with |ζ|<1 such that the circular unit disc, D, is defined as
(21)D≡{(s,ϕ):0≤s<1,-π≤ϕ<π}.
Using the properties of Riemann Mapping theorem [36], a circular disc D is mapped conformally onto Ω by
(22)z=af(ζ)for|ζ|<1,
where ζ=ξ+iη=seiϕ, ζ0=ξ0+iη0=s0eiϕ0, x=au(ξ,η), and y=av(ξ,η). Let
(23)w(x(ζ),y(ζ))=a|f′(ζ)|-1/2ejδW(ξ,η),q(x(ζ0),y(ζ0))=a|f′(ζ0)|-3/2ejδ0Q(ξ0,η0),
and the analytic function f in (22) is known to exist for any simply connected domain Ω. Further, we assume that |f′(ζ)| is nonzero and bounded for all |ζ|<1. Define S, Φ, δ, and δ0 as
(24)SeiΦ=ζ-ζ0,f′(ζ)=|f′(ζ)|eiδ,f′(ζ0)=|f′(ζ0)|eiδ0.
Let z-z0=a(f(ζ)-f(ζ0))=ReiΘ such that, for small S, R≃aS|f′(ζ0)| and Θ≃Φ+δ0.A similar integral equation with (18) can be obtained by substituting (23) and (24) into (18); that is,
(25)Q(ξ0,η0)=2-ν+3νe2jΘ8π⨘DW(ξ,η)S3dξdη+2-ν8π⨍DW(ξ,η)K(1)(ζ,ζ0)dξdη+3ν8π∫D‍W(ξ,η)K(2)(ζ,ζ0)dξdη;(ξ0,η0)∈D,
where K(1)(ζ,ζ0) and K(2)(ζ,ζ0) are Cauchy type singular and weak singular kernel, respectively [29]:
(26)K(1)(ζ,ζ0)=|f′(ζ)|3/2|f′(ζ0)|3/2|f(ζ)-f(ζ0)|3ej(δ-δ0)-1|ζ-ζ0|3,K(2)(ζ,ζ0)=|f′(ζ)|3/2|f′(ζ0)|3/2|f(ζ)-f(ζ0)|3ej(2Θ-δ-δ0)-1|ζ-ζ0|3e2jΦ.
This transformed hypersingular integral equation (25) over a circular disc D is solved subject to W=0 on s=1.


## 4. Numerical Treatment

Define(27)Akn(s,ϕ)=s|n|C2k+1|n|+(1/2)(1-r2)ejnϕ,Lhm(s,ϕ)=s|m|C2h+1|m|+(1/2)(1-r2)cos⁡mϕ
such that the orthogonal polynomials Akn(s,ϕ) and Lhm(s,ϕ) are satisfying the following relationship [37, Page 1054, 8.939.8]:
(28)∫Ω‍Akn(s,ϕ)Lhm(s,ϕ)sdsdϕ1-s2=Bknδkhδmn,
where the respective weight function is w(s)=(1-s)-1/2 and
(29)Bkn={π2Γ(2k+2)(2k+(3/2))(2k+1)![Γ(1/2)]2,n=0π2Γ(2k+2n+2)22n+1(2k+n+(3/2))(2k+1)![Γ(n+(1/2))]2,n≠0.
Write W(ξ,η) as a finite sum
(30)W(ξ,η)=∑n=-N1N1‍∑k=0N2‍s|n|WknC2k+1|n|+(1/2)(1-s2)ejnϕ.
Substituting (30) into (25) yields
(31)∑n,k‍Fkn(s0,ϕ0)Wkn=Q(ξ0(s0,ϕ0),η0(s0,ϕ0)),
where
(32)Fkn(s0,ϕ0)=-Ekn(2-ν+3νe2jΘ)Akn(s0,ϕ0)21-s02+2-ν8π∫D‍Akn(s,ϕ)K(1)(ζ,ζ0)dξdη+3ν8π∫D‍Akn(s,ϕ)K(2)(ζ,ζ0)dξdη;mmmmmnm0≤s≤1,0≤ϕ<2π.
The following formula [38] is useful in deriving (31):
(33)14π⨘ΩAkn(s,ϕ)R3dΩ=-EknAkn(s0,ϕ0)1-s02,
where
(34)Ekn=Γ(|n|+k+(3/2))Γ(k+(3/2))(|n|+k)!k!.
To determine the unknown coefficients, Wkn, multiply (31) by Lhm(s0,ϕ0) and integrate over D and using (28), leads to
(35)∑n,k‍Wkn(-2-ν+3νe2jΘ2δhkδ|m||n|mmmm+18π∫D‍Lhm(ζ0)∫D‍Akn(ζ)mmmmmmmmmmmminm×[(2-ν)K(1)(ζ,ζ0)mmmmmmmmmmmmminm+3νK(2)(ζ,ζ0)]dζdζ02-ν+3νe2jΘ2)=∫D‍Q(ζ0)Lhm(ζ0)dζ0;-N1≤m≤N1,0≤h≤N2
with the following notations:
(36)ζ0=ζ0(s0,ϕ0),dζ0=s0ds0dϕ0,Q(ζ0)=Q(ξ0,η0)=Q(s0cos⁡ϕ0,s0sinϕ0).
In evaluating the multiple integral in (35), we have used the Gaussian quadrature and trapezoidal formulas for the radial and angular directions with appropriate choice of collocation points (s,ϕ) and (s0,ϕ0). This effort leads to the (2N1+1)(N2+1)×(2N1+1)(N2+1) system of linear equations, Aw=b for the unknown coefficients of Wkn, where A=(aij) is a square matrix and w and b are vectors, and, solved numerically using LAPACK routine F07ASF (ZGETRS) in Numerical Algorithms Group (NAG).


## 5. Stress Intensity Factors, Maximum Stress Intensity, and Energy Release Rate

The sliding mode,K2(ϕ), and the tearing mode, K3(ϕ), stress intensity factors are defined as [39, 40]
(37)Kj(ϕ)=lim⁡r→a2πa-rVjw(x,y);j=2,3,
where Vj are constants. The maximum stress intensity, M(ϕ), is defined as
(38)M(ϕ)=[K2(ϕ)]2+[K3(ϕ)]2
while the energy release rate, G(ϕ), by Irwin’s relation subjected to shear loading is calculated from stress intensity factors and defined as
(39)G(ϕ)=(1-ν2)E[K2(ϕ)]2+(1+ν)E[K3(ϕ)]2,
where E, Young’s modulus, is a measurement of the stiffness of an isotropic elastic material and the relationship between E, ν, and μ is
(40)ν=E2μ-1.
Let a(ϕ)=|f(eiϕ)| and r=|f(seiϕ)|, followed by substituting (30) into (37), which leads to
(41)Kj(ϕ)=Vjlim⁡s→1-2π1-s|f′(ζ)|-1∑n,k‍WknAkn(s,ϕ);j=2,3,
where |f(eiϕ)-f(seiϕ)|=(1-s)|f′(eiϕ)| as s→1. Introduce
(42)Ykn(ϕ)=D2k+1|n|+1/2(0)cos⁡(nϕ),C2k+1|n|+1/2(1-s2)=1-s2D2k+1|n|+1/2(1-s2),
where Dmλ(x) is defined recursively by
(43)mDmλ(x)=2(m+λ-1)xDm-1λ(x)-(m+2λ-2)Dm-2λ(x);m=2,3,4,…,
with D0λ(x)=2λ and D1λ(x)=2λx. Substituting (42) into (41) yields
(44)Kj(ϕ)=2πVj|f′(eiϕ)|-1∑n,k‍WknYkn(ϕ),
where the unknown coefficients, Wkn, are obtained from (35).


## 6. Results and Discussion

Consider the conformal mapping [15]
(45)f(ζ)=ζ+cζm+1,
where m is an integer and c must satisfy -(1/m)≤c≤(1/m). The domain is circular if c=0 and has a smooth, regular boundary for 0≤(m+1)|c|<1. As (m+1)|c|→1, one or more cusps develop; see Figure 2 for various c and m, respectively.The domain off(ζ) with various m and c.
(a)
m
=
1
(b)
m
=
4Tables1, 2, and 3 show that our numerical scheme converges rapidly at a different point of the crack with only a small value of N=N1=N2 used. Tables 1 and 2 show numerical scheme for K2(ϕ) for m=1 with c=0.1 and c=0.45, respectively. Based on these two tables, it is evident that the convergence of stress intensity factors becomes slow as c increases, whilst Table 3 presents the numerical scheme for K3(ϕ), for m=-2 with c=0.1.Table 1
Numerical convergence for the sliding mode stress intensity factor,K2(ϕ), for m=1 when c=0.1.
N
K
2
(
0.00
)
K
2
(
π
/
4
)
K
2
(
π
/
2
)
K
2
(
3
π
/
4
)
K
2
(
π
)
0
1.0423
E
-
03
7.5379
E
-
04
6.9319
E
-
20
-
8.6663
E
-
04
-
1.2782
E
-
03
1
1.4538
0.9883
6.7158
E
-
17
-
1.6978
-
0.9411
2
1.3545
0.8915
6.7071
E
-
17
-
0.9971
-
0.9427
3
1.3333
0.8914
6.4321
E
-
15
-
0.9716
-
0.9426
4
1.3191
0.8637
-
0.1194
-
0.9709
-
1.2468
5
1.3191
0.8637
-
0.1194
-
0.9709
-
1.2468
6
1.3191
0.8637
-
0.1194
-
0.9709
-
1.2468Table 2
Numerical convergence for the sliding mode stress intensity factor,K2(ϕ), for m=1 when c=0.45.
N
K
2
(
0.00
)
K
2
(
π
/
4
)
K
2
(
π
/
2
)
K
2
(
3
π
/
4
)
K
2
(
π
)
0
5.9156
E
-
04
4.3514
E
-
04
4.35045
E
-
20
-
6.7348
E
-
04
-
1.2782
E
-
03
1
0.0000
1.0959
3.1941
E
-
19
-
8.7363
E
-
03
-
3.3449
E
-
02
2
1.7569
1.0959
6.3681
E
-
17
-
1.3444
-
0.3384
3
1.7624
1.0930
6.3681
E
-
17
-
1.3489
-
0.3143
4
1.5604
1.0945
6.3681
E
-
17
-
1.3511
-
1.3057
5
1.4613
0.9940
6.3681
E
-
17
-
1.3518
-
1.3020
6
1.4609
0.9812
6.3681
E
-
17
-
1.3518
-
1.3002
7
1.4011
0.8127
-
0.6744
-
1.3517
-
1.3993
9
1.3610
0.8115
-
0.6744
-
1.7515
-
1.2989
10
1.2030
0.8029
-
0.6317
-
1.9999
-
1.4986
11
1.1948
0.7846
-
0.6000
-
1.1489
-
1.2984
12
1.1849
0.7796
-
0.5913
-
1.1317
-
1.1983
13
1.1782
0.6912
-
0.5410
-
1.1245
-
1.1203
14
1.1782
0.5680
-
0.5302
-
1.0982
-
1.1782
15
1.1782
0.5680
-
0.5302
-
1.0982
-
1.1782
17
1.1782
0.5680
-
0.5302
-
1.0982
-
1.1782
16
1.1782
0.5680
-
0.5302
-
1.0982
-
1.1782Table 3
Numerical convergence for the tearing mode stress intensity factor,K3(ϕ), for m=-2 when c=0.1.
N
K
3
(
0.00
)
K
3
(
π
/
4
)
K
3
(
π
/
2
)
K
3
(
3
π
/
4
)
K
3
(
π
)
0
0.0000
-
6.330
E
-
04
-
9.5075
E
-
04
-
7.2780
E
-
04
-
1.3145
E
-
19
1
0.0000
-
0.7174
-
0.9211
-
0.5854
-
9.6785
E
-
17
2
0.0000
-
0.5440
-
0.9199
-
0.5854
-
9.6951
E
-
17
3
0.0000
-
0.5200
-
0.8775
-
0.200
-
1.04460
E
-
16
4
0.0000
-
0.5200
-
0.8775
-
0.5200
-
1.04460
E
-
16
5
0.0000
-
0.5200
-
0.8775
-
0.5200
-
1.04460
E
-
16
6
0.0000
-
0.5200
-
0.8775
-
0.5200
-
1.04460
E
-
16Figures3, 4, 5, and 6 display the comparison of asymptotic and numerical solutions for K2(ϕ) and K3(ϕ) stress intensity factors, maximum stress intensity, M(ϕ), and strain energy release rate, G(ϕ), respectively, for m=1 at c=0.1 and c=0.3. As demonstrated in these figures, our results seem to agree with those obtained by Gao [40] except at the cusps. As the cusps become sharper, the analytical result by Gao [40] does not work, and this gives rise to the difference between our and Gao’s [40] result. It can be seen that the stress intensity factors have local extremal values when the crack front is at cos⁡(ϕ)=±1 or sin(ϕ)=±1. Figure 7 presents the comparison of asymptotic and numerical solutions for K2(ϕ) and K3(ϕ) stress intensity factors, maximum stress intensity, M(ϕ), and strain energy release rate, G(ϕ), respectively, for m=2 at c=0.1. Figure 8 shows the variations of K2, K3, M, and G against ϕ for various of c≤0.3. Similar behavior can be observed for the solution of K2(ϕ), K3(ϕ), M(ϕ), and G(ϕ), for a different parameter of ν for f(ζ)+0.1ζ4, displayed in Figure 9.TheK2(ϕ) for f(ζ)=ζ+cζ2 at different c.
(a)
c
=
0.1
(b)
c
=
0.3TheK3(ϕ) for f(ζ)=ζ+cζ2 at different c.
(a)
c
=
0.1
(b)
c
=
0.3TheM(ϕ) for f(ζ)=ζ+cζ2 at different c.
(a)
c
=
0.1
(b)
c
=
0.3TheG(ϕ) for f(ζ)=ζ+cζ2 at different c.
(a)
c
=
0.1
(b)
c
=
0.3TheK2(ϕ), K3(ϕ), M(ϕ), and G(ϕ) for f(ζ)=ζ+cζ3 at c=0.1.
(a)
K
2
(
ϕ
)
(b)
K
3
(
ϕ
)
(c)
M
(
ϕ
)
(d)
G
(
ϕ
)TheK2(ϕ), K3(ϕ), M(ϕ), and G(ϕ) for f(ζ)=ζ+cζ3 at various c.
(a)
K
2
(
ϕ
)
(b)
K
3
(
ϕ
)
(c)
M
(
ϕ
)
(d)
G
(
ϕ
)TheK2(ϕ), K3(ϕ), M(ϕ), and G(ϕ) for f(ζ)=ζ+0.1ζ4 at various ν with μ=1.
(a)
K
2
(
ϕ
)
(b)
K
3
(
ϕ
)
(c)
M
(
ϕ
)
(d)
G
(
ϕ
)


## 7. Conclusion

The present work dealt with the epicycloid crack with the application of shear loading in fracture mechanics. To this end, the numerical solution for the stresses in such specimens is derived, based on the solution of the hypersingular integral equation, and the conformal mapping technique is adopted to transform the hypersingular integral equation over a circular region such that the equation is reduced into a system of linear equations and solved for the unknown coefficients. The stress intensity factors, maximum stress intensity, and strain energy release rate for the epicycloid crack subject to shear load are presented graphically. The proposed model and the obtained numerical results are in good agreement when compared to Gao [40].



---

*Source: 213478-2014-07-14.xml*