Axisymmetric Torsion of an Elastic Layer Sandwiched between Two Elastic Half-Spaces with Two Interfaced Cracks

The class of problems related to the behavior of rigid disc inclusions embedded in bonded contact with an elastic medium, has been a subject of much interest in geomechanics, civil engineering, and applied mechanics. This work is motivated by both theoretical and practical interests in the problems of turbines disks, some pipes, and many industrial applications. It may give a better understanding of the behavior of foundations under external loads. In structure-medium interaction problems arising in foundation engineering, the foundation is usually modeled using a rigid or flexible inclusion having circular, strip, rectangular, or arbitrary shape. Nowadays, composites play a very important role in geomechanics engineering.

It is common knowledge that all existing structural materials contain different inter- and intra-component defects (cracks, delaminations, etc.) [1]. The problem of the torsion of an infinite elastic medium by a rigid inclusion (deeply embedded) was considered by Selvadurai [2, 3]. His results depend on the rotational and translational stiffnesses of the embedded rigid circular disc. The problem of the torsion of an elastic half- space was considered, at first, by Reissner and Sagoci [4]. They studied the static interaction of a rigid disc and an elastic isotropic half-space for which they obtained the solution by means of the spheroidal coordinates. The same problem was solved by Sneddon [5] using a different method. He used the Hankel transforms method for reducing the problem to a pair of dual integral equations. Collins [6] treated the torsional problem of an elastic half-space by assuming the displacement at any point in the half- space to be due to a distribution of wave sources over the part of the free surface in contact with the disc. The solution of the forced vibration problem of elastic layer of finite thickness when the lower face is either stress free or rigidly clamped was given by Gladwell [7]. Pak and Saphores [8] provided an analytical formulation for the general torsional problem of a rigid disc embedded in an isotropic half-space. Besides, Bacci and Bennati [9] considered the torsional of circular rigid disc adhered to the upper surface of an elastic layer fixed to an undefonnable support. More recently, Singh et al. [10] studied the torsional of a non-homogeneous, isotropic, half- space by rotating a circular part of its boundary surface. Cai and Zue [11] discussed the torsional vibration of a rigid disc bonded to a poroelastic multilayered medium.. Yu [12] studied the forced torsional oscillations inside the multilayered solid. The elastodynamic Green’s function of the center of rotation and a point load method were used to solve the problem. Pal and Mandal [13] considered the forced torsional oscillations of a transversely isotropic elastic half- space under the action of an inside rigid disc. A similar problem with the rocking rotation was solved later on by Ahmadi and Eskandari [14]. They used an appropriate Green’s function to write the mixed boundary -value problem posed as a dual integral equation. All these problems are based on the theory of Hankel integral transformation in order to bring the mixed boundary -value problem into a system of dual integral equations. Then, the corresponding solution is sought from an integral equation of Fredholm type.

The torsional of elastic layers with a penny -shaped crack was considered by some researchers. Sih and Chen [15] studied the problem of a penny-shaped crack in layered composite under a uniform torsional stress. The displacement and stress fields throughout the composite were obtained by solving a standard Fredholm integral equation of the second kind. Low [16] investigated a problem of the effects of embedded flaws in the form of an inclusion or a crack in an elastic half- space subjected to torsional deformations. The corresponding Fredholm integral equations were solved numerically by quadrature approach. The same method was used by Dhawan [17] for solving the problem of a rigid disc attached to an elastic half-space with an internal crack. By using Hankel and Laplace transforms and taking numerical inversion of Laplace transform, Basu and Mandal [18] treated the torsional load on a penny-shaped crack in an elastic layer sandwiched between two elastic half-spaces. The purpose of this article is to study an axisymmetric torsion of an embedded circular rigid disc in bonded contact with an isotropic elastic layer sandwiched between two elastic half-spaces with two penny-shaped cracks symmetrically located at each of the two interfaces between the layer and the half-spaces. A similar method was used in arecently published work by Madani and Kebli [20], dealing with the case of a penny-shaped crack problem in the interior of a homogeneous elastic material at the symmetry plane, under an axisymmetric torsion by two circular rigid discs symmetrically located in the elastic medium.

In view of the axial symmetry of the problem, it is natural to consider (r, θ, z) the cylindrical polar co-ordinates. Here we consider an axisymmetric torsion of an embedded circular rigid disc with a radius b at the symmetry plane z = 0 in an isotropic and homogeneous layer of thickness 2h and with one material sandwiched between two half-spaces of the second material. Two penny -shaped cracks with a radius a were symmetrically located at the interface between the layer and the half- spaces z = ±h. The faces of the cracks are assumed to be stress free while the discs rotate with an equal angle ω about the z- axis passing through their center as shown in Figure .1.

Figure 1

Owing to symmetry about the z = 0 plane, it is sufficient to consider the problem in the upper half- space where z ≥ 0. In this case of an axisymmetric torsion problem, the displacement vector assumes the form (0, u_θ, 0) in the cylindrical polar coordinate system (r, θ, z). It is convenient to identify a layer region (superscript (1)) occupying the region r ∈ (0, ∞); z ∈ (0, h⁻) and a layer region (superscript (2)) occupying the region r ∈ (0, ∞); z ∈ (h⁻, ∞). The signs h⁺ and h⁻” denote the variables in the upper surface and the lower surface of plane z = h, respectively.

The only non-zero components of stress are given by

where u_θ = u_θ (r, z) and Gi is the shear modulus of the material. As the torsion of the homogeneous material is static, the displacement u_θ (r, z) must satisfy

By means of Hankel’s transformation, integral and its inverse given in [19],

And

The solution of equation (2) for the regions I (0 ≤ z ≤ h) and I I (z ≥ h) is expressed as

where λ is the transform variable, J1 is the Bessel function of the first kind of order one, and Ai and Bi are unknown functions

We consider the regularity conditions at infinity, the symmetry plane condition at z = 0, and the boundary and continuity conditions at the bonded interfaces z=h+$z={{h}^{+}}$and h−${{h}^{-}}$. Therefore, we find the following conditions

Applying the regularity conditions at infinity given in Eq. (4a), we obtain

where A₁(λ), B₁(λ), and A₂(λ) are arbitrary functions that need to be determined by satisfying the boundary and continuity conditions.

The boundary and the continuity conditions in Eqs. (4c) and (4e) lead to

The above equation implies

where γ=G1G2.$\gamma =\frac{{{G}_{1}}}{{{G}_{2}}}.$

The mixed boundary conditions Eqs. (4e), (4d), (4f), and (4b) are satisfied if A₁ and B₁ are solutions of the following dual integral equations:

Equations 8b and 8d are identically satisfied if we introduce the following representation:

where J12${{J}_{\frac{1}{2}}}$and J32${{J}_{\frac{3}{2}}}$are the Bessel functions of the first kind of order 12$\frac{1}{2}$and 32$\frac{3}{2}$, respectively. The unknown functions are given by the following equation:

where p(λ)=(1−γ1+γ)e−2λh+1andϕ(t)andψ(t)$p\left( \lambda \right)=\left( \frac{1-\gamma }{1+\gamma } \right){{e}^{-2\lambda h}}+1\,\,\,\,\,\text{and}\,\,\,\,\phi \left( t \right)\,\,\,\text{and}\,\,\,\psi \left( t \right)$are continuous unknown functions of t defined over two intervals 0≤t<a$0\le t<a$and 0≤t<b$0\le t<b$, respectively. Substituting A1(λ) and B1(λ) in Eqs. (8a) and (8c), we get

where

The expression for p(λ) approaches to 1 for large values of λ. Equation 14 can be converted to the Abel integral equation by means of the relation λJ1(λr)=1r2ddr[r2J2(λr)],$\lambda {{J}_{1}}\left( \lambda r \right)=\frac{1}{{{r}^{2}}}\frac{d}{dr}\left[ {{r}^{2}}{{J}_{2}}\left( \lambda r \right) \right],$and then, by taking into account the integral formula

we obtain the Abel equation corresponding to Eq. (14)

Next, we invert the last equation by applying the Abel transform formula

to obtain

For the left-hand side of the above equation, the integral is further simplified by using the following relationship:

we obtain the first Fredholm integral equation of the second kind

where

Following the similar procedure as before, Eq. (15) is reduced to the second Fredholm integral equation. Using the following formula

we obtain the following Abel-type equation:

Now, we invert the above equation by applying the Abel transform formula to get

Using the following relationships

we get the second Fredholm integral equation of the second kind

with the kernel

To get a non-dimensionalized equation, from Eqs. (18) to (22), let us change the variables as follows

Next, we multiply the above two equations of the system by 2π4aωφ(au)and2π4bωψ(bu),$\frac{\sqrt{2\pi }}{4a\omega }\varphi \left( au \right)\,\text{and}\,\frac{\sqrt{2\pi }}{4b\omega }\psi \left( bu \right),$respectively, and using the following substitutions

we obtain the following equations:

where

As the kernels K, L, M, and N are continuous on the interval [1], the system of Fredholm integral equations can be solved by direct or iterative techniques [21]. The midpoint quadrature [22] is used to find the numerical solution for the system given by Eqs. (25a) and (25b). By dividing the interval [1] into N equal subintervals, so that the midpoints are u=um=2m−12s=un=2n−12m,n=1,2...N.,$u={{u}_{m}}=\frac{2m-1}{2}\,\,\,\,\,s={{u}_{n}}=\frac{2n-1}{2}\,\,\,\,m,n=1,\,\,\,2...N.,$and by introducing the following notations

We evaluate numerically the infinite integral K, L, M, and N using the Simpson rule. After solving the above system, the unknown coefficients can be obtained.

5.1

Stress intensity factors

The stress intensity factors at the edge of the crack and at the rim of the disc are defined, respectively, by the following equations:

(28a)KIIIa=limr→a+2π(r−a)τzθ1(r,z)|z=H$$K_{III}^{a}=\underset{r\to {{a}^{+}}}{\mathop{\lim }}\,\sqrt{2\pi \left( r-a \right)}\tau _{z\theta }^{1}\left( r,z \right){{|}_{z=H}}$$(28b)KIIIb=limr→b−2π(r−b)τzθ1(r,z)|z=0$$K_{III}^{b}=\underset{r\to {{b}^{-}}}{\mathop{\lim }}\,\sqrt{2\pi \left( r-b \right)}\tau _{z\theta }^{1}\left( r,z \right){{|}_{z=0}}$$

On the planes z = h for r ≥ a and z = 0, the expressions of stress are given by

(29)τzθ1(r,θ)=G1∫0∞λ32[1−e−2λhp(λ)(1+γ)∫0atϕ(t)J32(λt)dt+(eλh(1−1p(λ))−eλhp(λ))∫0btψ(t)J12(λt)dt]J1(λr)dλ$$\begin{align}& \tau _{z\theta }^{1}\left( r,\theta \right)={{G}_{1}}\int_{0}^{\infty }{{{\lambda }^{\frac{3}{2}}}}\left[ \frac{1-{{e}^{-2\lambda h}}}{p\left( \lambda \right)\left( 1+\gamma \right)} \right.\int_{0}^{a}{\sqrt{t}\phi \left( t \right){{J}_{\frac{3}{2}}}}\left( \lambda t \right)dt+ \\ & \left. \left( {{e}^{\lambda h}}\left( 1-\frac{1}{p\left( \lambda \right)} \right)-\frac{{{e}^{\lambda h}}}{p\left( \lambda \right)} \right)\int_{0}^{b}{\sqrt{t}\psi \left( t \right){{J}_{\frac{1}{2}}}\left( \lambda t \right)dt} \right]{{J}_{1}}\left( \lambda r \right)d\lambda \\ \end{align}$$(30)τzθ1(r,h)=G1∫0∞λ32∫0b[tψ(t)J12(λt)dt]J1(λr)dλ$$\tau _{z\theta }^{1}\left( r,h \right)={{G}_{1}}\int_{0}^{\infty }{{{\lambda }^{\frac{3}{2}}}}\int_{0}^{b}{\left[ \sqrt{t}\psi \left( t \right){{J}_{\frac{1}{2}}}\left( \lambda t \right)dt \right]}{{J}_{1}}\left( \lambda r \right)d\lambda $$

The second part of the integrals in Eq. (29) converge quickly as their limits r → a and r → b automatically vanishes; however, the limits of the other two integrals analyzed asymptotically as follows. Using the relation J1(λR)=−1RddRJ0(λR),${{J}_{1}}\left( \lambda \,R \right)=-\frac{1}{R}\frac{d}{dR}{{J}_{0}}\left( \lambda \,R \right),$we obtain

(31)τzθ1(r,θ)=−G1∫0atϕ(t)dt∫0∞λ12p(λ)(1+γ)J32(λt)J0(λr)dλ−G1∫0∞λ32[e−2λhp(λ)(1+γ)∫0atϕ(t)J32(λt)dt+(eλh(1−1p(λ))−eλhp(λ))∫0btψ(t)J12(λt)dt]J1(λr)dλ$$\begin{align}& \tau _{z\theta }^{1}\left( r,\theta \right)=-{{G}_{1}}\int_{0}^{a}{\sqrt{t}\phi \left( t \right)dt\int_{0}^{\infty }{\frac{{{\lambda }^{\frac{1}{2}}}}{p\left( \lambda \right)\left( 1+\gamma \right)}{{J}_{\frac{3}{2}}}\left( \lambda t \right){{J}_{0}}\left( \lambda r \right)d\lambda -{{G}_{1}}\int_{0}^{\infty }{{{\lambda }^{\frac{3}{2}}}\left[ \frac{{{e}^{-2\lambda h}}}{p\left( \lambda \right)\left( 1+\gamma \right)} \right.}}} \\ & \left. \int_{0}^{a}{\sqrt{t}\phi }\left( t \right){{J}_{\frac{3}{2}}}\left( \lambda t \right)dt+\left( {{e}^{\lambda h}}\left( 1-\frac{1}{p\left( \lambda \right)} \right)-\frac{{{e}^{\lambda h}}}{p\left( \lambda \right)} \right)\int_{0}^{b}{\sqrt{t}\psi \left( t \right){{J}_{\frac{1}{2}}}\left( \lambda t \right)dt} \right]{{J}_{1}}\left( \lambda r \right)d\lambda \\ \end{align}$$(32)τzθ1(r,h)=−G1∫0btψ(t)dt∫0∞λ12J12(λt)J0(λr)dλ$$\tau _{z\theta }^{1}\left( r,h \right)=-{{G}_{1}}\int_{0}^{b}{\sqrt{t}}\psi \left( t \right)dt\int_{0}^{\infty }{{{\lambda }^{\frac{1}{2}}}}{{J}_{\frac{1}{2}}}\left( \lambda t \right){{J}_{0}}\left( \lambda r \right)d\lambda $$

For large values of λ, we use the following asymptotic behavior of the Bessel function of the first Jv(λ)≅2λπcos(λ−vπ2−π4)${{J}_{v}}\left( \lambda \right)\cong \sqrt{\frac{2}{\lambda \pi }}\cos \left( \lambda -v\frac{\pi }{2}-\frac{\pi }{4} \right)$and using the following integral formulas for the first infinite integral in the right part of the Eqs. (31) and (32), respectively,

(33)∫0∞cos(λt)J0(λr)dλ={1(r2−t2)t<r0t>rand∫0∞sin(λt)J0(λr)dλ={0t<r1(t2−r2)t>r$$\begin{align}& \int_{0}^{\infty }{\cos }\left( \lambda t \right){{J}_{0}}\left( \lambda r \right)d\lambda =\left\{ \begin{matrix}\frac{1}{\sqrt{\left( {{r}^{2}}-{{t}^{2}} \right)}} & t<r \\0 & t>r \\\end{matrix} \right.\text{and} \\ & \,\,\,\,\int_{0}^{\infty }{sin}\left( \lambda t \right){{J}_{0}}\left( \lambda r \right)d\lambda =\left\{ \begin{matrix}0 & t<r \\\frac{1}{\sqrt{\left( {{t}^{2}}-{{r}^{2}} \right)}} & t>r \\\end{matrix} \right. \\ \end{align}$$

As p(λ) → 1 as λ → ∞, Eqs. (31) and (32) become

(34a)τzθ1(r,θ)=G1γ+12πddr∫0aϕ(t)(r2−t2)dt+R1(r)$$\tau _{z\theta }^{1}\left( r,\theta \right)=\frac{{{G}_{1}}}{\gamma +1}\sqrt{\frac{2}{\pi }\frac{d}{dr}}\int_{0}^{a}{\frac{\phi \left( t \right)}{\sqrt{\left( {{r}^{2}}-{{t}^{2}} \right)}}}dt+{{R}_{1}}\left( r \right)$$(34b)τzθ1(r,h)=−G12πddr∫0bψ(t)(t2−r2)dt$$\tau _{z\theta }^{1}\left( r,h \right)=-{{G}_{1}}\sqrt{\frac{2}{\pi }\frac{d}{dr}}\int_{0}^{b}{\frac{\psi \left( t \right)}{\sqrt{\left( {{t}^{2}}-{{r}^{2}} \right)}}}dt$$

where

(35)R1(r)=−G1∫0∞λ32[e−2λhp(λ)(1+γ)∫0atϕ(t)J32(λt)dt+(eλh(1−1p(λ))−eλhp(λ))∫0btψ(t)J12dt]J1(λr)dλ$$\begin{align}& {{R}_{1}}\left( r \right)=-{{G}_{1}}\int_{0}^{\infty }{{{\lambda }^{\frac{3}{2}}}}\left[ \frac{{{e}^{-2\lambda h}}}{p\left( \lambda \right)\left( 1+\gamma \right)}\int_{0}^{a}{\sqrt{t}\phi \left( t \right){{J}_{\frac{3}{2}}}\left( \lambda t \right)dt}+ \right. \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left. \left( {{e}^{\lambda h}}\left( 1-\frac{1}{p\left( \lambda \right)} \right)-\frac{{{e}^{\lambda h}}}{p\left( \lambda \right)} \right)\int_{0}^{b}{\sqrt{t}\psi \left( t \right){{J}_{\frac{1}{2}}}dt} \right]{{J}_{1}}\left( \lambda r \right)d\lambda \\ \end{align}$$

Now integrating by parts, we get

(36a)τzθ1(r,0)=G1γ+12π[aϕ(a)r(r2−a2)∫0atϕ′(t)(r2−a2)dt]+R1(r)$$\tau _{z\theta }^{1}\left( r,0 \right)=\frac{{{G}_{1}}}{\gamma +1}\sqrt{\frac{2}{\pi }}\left[ \frac{a\phi \left( a \right)}{r\sqrt{\left( {{r}^{2}}-{{a}^{2}} \right)}}\int_{0}^{a}{\frac{t{\phi }'\left( t \right)}{\sqrt{\left( {{r}^{2}}-{{a}^{2}} \right)}}dt} \right]+{{R}_{1}}\left( r \right)$$

We note that the infinite integrals in the preceding expressions are convergent throughout the medium except at the singular points r → a⁺, which occupy the crack boundary.

(36b)τzθ1(r,h)=G12π[bψ(b)r(b2−r2)−∫rbtψ′(t)(r2−t2)dt]$$\tau _{z\theta }^{1}\left( r,h \right)={{G}_{1}}\sqrt{\frac{2}{\pi }}\left[ \frac{b\psi \left( b \right)}{r\sqrt{\left( {{b}^{2}}-{{r}^{2}} \right)}}-\int_{r}^{b}{\frac{t{\psi }'\left( t \right)}{\sqrt{\left( {{r}^{2}}-{{t}^{2}} \right)}}dt} \right]$$

In this case, the integral in the above relation converges quickly and the integral is bonded as r → b⁻. As a result, we obtain a square root singularity at r = b and the constant ψ(b) is the measure of the strength of singularity at the vicinity of the rigid inclusion. By using the following transformations ϕ(a)=4aω2πΦN,$\phi \left( a \right)=\frac{4a\omega }{\sqrt{2\pi }}{{\Phi }_{N}},$ψ(b)=4bω2πΨN.$\psi \left( b \right)=\frac{4b\omega }{\sqrt{2\pi }}{{\Psi }_{N}}.$.

we obtain the stress intensity factor at the edge of the crack and at the rim of the disc

(37a)KIIIa=−4G1ωa(1+γ)πΦN$$K_{III}^{a}=\frac{-4{{G}_{1}}\omega \sqrt{a}}{\left( 1+\gamma \right)\sqrt{\pi }}{{\Phi }_{N}}$$(37b)KIIIb=4G1ωb(1+γ)πΨN$$K_{III}^{b}=\frac{4{{G}_{1}}\omega \sqrt{b}}{\left( 1+\gamma \right)\sqrt{\pi }}{{\Psi }_{N}}$$

Figure 2 shows the variation of the normalized stress intensity factor KIIIa$K_{III}^{a}$at the edge of the crack defined by Eq. (37a) against c for various values of the layer thickness H = 1, 0.75, 0.5, 0.25 calculated using the shear modulus ratio γ = 1. It is observed that the values of stress intensity factor increase and attain its maximum values at c = 1, and with the increase in the value of c, the stress intensity factor decreases. In addition, the effect of the axial distance between the crack and the disc H on the stress intensity factor is also shown in this figure. The increase in the distance H induces the decrease in stress intensity factor for all the values of parameter c. Figure 3 illustrates the variation of the stress intensity factor at the crack because of variations in shear modulus ratio γ for different values: the normalized crack size c = 0.25, 0.5, 0.75, 1 and the layer thickness H = 1. We observe from the figure that as the shear parameter γ increases, the stress intensity factor decreases for all values of c.

Figure 2

Figure 3

Figure 4 illustrates the variation of the normalized stress intensity factor KIIIb$K_{III}^{b}$at the edge of the rigid inclusion defined by Eq. (37b) versus c for H = 1, 0.75, 0.5, and 0.25. Relatively, small variation for smaller values of c and considerable variation for larger values of c are observed. Also, the interaction between the cracks and the rigid disc is greater when the cracks are closer to the disc. In addition to the interaction, the stress intensity factor KIIIb$K_{III}^{b}$KIIIb$K_{III}^{b}$increases as the crack radius increases. From the

Figure 4

formulation and the presented figures, conclusions may be deduced:

Singularity at the edge of the internal crack and the internal rigid inclusion is observed (The results seem to agree with the previous works (Low [16], Dhawan [17], and Madani and Kebli [20]).

The Mode III stress intensity factor KIIIb$K_{III}^{b}$at the edge of the crack is negative, and it decreases with increasing γ.

There is considerable interaction between the cracks and the rigid inclusion when a/b is large and the cracks are close to the rigid inclusion.

5.3

Displacement and stress fields

The results for the variation of the normalized displacement u⁽ⁱ⁾(ρ, ξ)/ωa and stress τ⁽ⁱ⁾(ρ, ξ)/Gⁱωa with ρ = r/b are shown graphically in Figures 5–8 for the different values of the dimensionless axial distances ξ = z/b. For each region, five different axial distances are selected as I (ξ =0; H/5; 2H/5; 3H/5; H) and I I (ξ = H ; 6H/5; 7H/5; 8H/5; 2H), with the particular values of the height H = 1, the dimensionless

Figure 5

crack size c = 1 and the shear parameter γ = 1 and γ = 2. The variation of the normalized displacements is shown in Figures 5 and 6. We notice that the displacements in the two regions increase at first, reach maximum values at ρ = c, and then decrease out of the disc band with increasing ρ. The distribution of the shear stresses in the elastic medium are also discussed and shown in Figures 7 and 8. The stresses initially increase, attain its maximum values, and then with the increase in the value of ρ, the stresses go on decreasing.

Figure 6

Figure 7

Figure 8

In this article, an axisymmetric torsion of a rigid disc embedded in the interior of a homogeneous elastic layer sandwiched between two half-spaces containing two interface cracks is analytically addressed. By using the Hankel integral transformation and its inverse, the mixed boundary value problem is reduced to a system of dual integral equations, which are reduced to a Fredholm integral equation system of the second kind. The numerically computed results of the displacements, the stresses, and the stress intensity factors are presented graphically for some dimensionless parameters. The numerical results show that the discontinuities around the crack and the inclusion cause a large increase in the stresses that decay with distance from the disc loaded. Furthermore, the dependence of the stress intensity factor on the disc size, the distance between the crack and the disc and the nonhomogeneity parameter is observed.