1 Introduction

The reality of problems on stress concentration in elastic bodies weakened by cracks causes the necessity of mathematical modeling of the material’s heterogeneity in the elements of building construction and in machine components. The necessity of the creation of mathematical methods pushes the development of new mathematical apparatus from the other side^[1-9]. The problems of torsion on elastic cylinders weakened by the cracks were investigated sufficiently both in the case of the circular cracks and for boundary cracks^[10-15].

The problems on the stress state of multilayered cylinders^[16-18] were investigated to a lesser degree. In Ref. [16], the problem on stress concentration near the interface cylindrical crack during torsion of a two-layered cylinder is reduced to the integral Cauchy’s equation. In Ref. [17], the torsion problem for the two layered cylinder with a circular crack is reduced to the integral Fredholm equation of the second kind. In the mentioned works, the problems were solved for the case of a two-layered cylinder. The solutions in both papers were done numerically. These approaches, leading to the numerical solution of the integral equations, are most common in the papers where the authors consider the problems for a multi-layered cylinder with a crack.

The advantage of the proposed method consists of the reduction of the number of constants regardless of the number of layers and the presentation of the mechanical characteristics in the form of the uniformly convergent series.

The main difficulty in solving stress problems for a multilayered media is the necessity to determine the unknown constants of the layers. Their quantity is proportional to the number of the layers and to the order of the corresponding differential equation^[19-20]. The effective method of the layered medium elasticity problem solving was proposed by Popov and Grivnyak^[21]. The novelty of the proposed paper consists of the application of this method leading to the reduction of the number of constants regardless of the number of layers. The solution to the problem is constructed for an arbitrary number of layers. With the help of this method, the problem of the torsion of a multilayered spherical medium^[22] and the problem of the torsion of elastic conical multilayered bodies^[23] were solved. The same approach was used for the construction of the exact solution for the case of a multilayered cylinder under the influence of torsion loading^[24].

In the proposed paper, the construction of the solution to the torsion problem for a circular finite multilayered elastic cylinder weakened by a circular crack is done with the approach given in Ref. [21]. The layers of the cylinder are the coaxial cylinders differing from each other by elastic constants. The conditions of an ideal contact are fulfilled between them. The circular crack is located in the internal layer. The problem’s solution is constructed for an arbitrary number of layers. The obtained singular integro-differential equation is solved using the orthogonal polynomials method. This method enables a solution to be found with regard to the real singularity of an unknown function and to calculate the stress intensity factor (SIF).

2 Problem’s statement

The elastic cylinder 0≤r≤a, -π≤φ≤π, 0≤z≤h (here (r, φ, z) is the cylindrical coordinate system) consists of n coaxial elastic cylindrical layers r_j-1≤r≤r_j, 0≤z≤h, j=1, 2, …, n, r₀=0, r_n=a. Every layer has its own shear module G_j and the Poisson’s coefficient µ_j, j=1, 2, …, n.

The axisymmetric tangent stress is applied to the lateral surface of the cylinder

(1)

Here and hereinafter the upper indexes near the stresses or the displacement symbols denote the number of a layer, p(z) is the known function. The bottom face of the cylinder is fixed

(2)

The upper face is free from stresses,

(3)

The layers of the cylinder are in the conditions of ideal mechanical contact, so the stresses and the displacements are continuous functions when passing through the layers

(4)

The circular crack 0 < r < c, z=d, c < r₁ parallel to the cylinder’s faces is located in the internal layer. The branches of the crack are free from stresses, and the displacement is discontinuous,

(5)

(6)

One needs to find the formulas for the cylinder’s stress state and SIF, when the displacements of each layer should satisfy the torsion equation

(7)

3 Reduction of problem to one-dimensional boundary problem and its solution in transformation’s domain

Before reduction of the problem to a one-dimensional boundary valued problem with the help of the integral transformation method, one needs to reformulate the boundary conditions in terms of the displacement, taking into consideration the known relations between stresses and the tangent displacement^[25], . The obtained problem is reformulated with dimensionless variables. The new variables and the designations are taken as with this aim. The following discontinuous boundary value problem is obtained:

(8)

(9)

(10)

(11)

The finite Fourier’s transformation is applied to the boundary problem (8)-(11) with regard to the variable ζ

(12)

This integral transformation is applied to the equation (8) when j=1 by the generalized scheme. The advantage of this method is in the possibility to use the scheme when the function has discontinuity with regard to the integration’s variable^[26] (see Appendix A). The usual classical scheme of integral transformation method is applied for other equations when j≥2. First two conditions from (9) are satisfied automatically during this procedure. One obtains instead the equation (8) in the transformations domain

(13)

(14)

After the transformation (12), the boundary condition (10) and the conditions of the conjugation (11) take the form

(15)

The general solutions of (13) and (14) are

(16)

(17)

where I₁(x) and K₁(x) are the modified Bessel’s functions, and Φ_k(ρ, ξ) is the fundamental function^[27] of (13)

(18)

This function satisfies all conditions for a fundamental function: it is limited when ρ=0, satisfies the homogenous equation on each interval [0; t) and (t; c₁], is continuous on each interval; its derivative has the jump equal to , by passing through the line ρ=t.

One needs to find the constants A_1k, A_jk, B_jk in the solutions (16)-(17). Let us use the method allowing the use of recurrent formulas for their estimation. This method was proposed first in Ref.[28]. The vectors and matrices are input with this aim,

where t_jk=αβ_kc_j. With the help of these matrices and vectors, the conjugation conditions (16) can be written in the forms

(19)

(20)

where . The formulas

were used during the obtaining of these correspondences^[29]. One gets after the solving of the matrix equations (19) and (20)

(21)

It is used during the calculation of the inverse matrix R_jk^-1 that I₁(x) K₂(x)+I₂(x) K₁(x)=x^-1. Finally, the inverse matrix is found in the form

For the reduction of the calculations, the following matrices are taken

With the help of these matrices, the formulas (21) can be rewritten as

(22)

Now one can find consecutively

After the construction of the matrix

the coefficients’ vector for the nth layer will be expressed through the coefficient’s vector of the first layer

The found coefficients A_nk and B_nk are substituted at the boundary conditions (15), which with regard to the formula (16), take the form

So one can find that

With the help of the recurrent formulas (22), one can consecutively find the coefficients X_2k, X_3k and so on. The expressions for the displacements’ transformations (10) are found. The inverse integral transformation (12) is used to obtain the displacement u_j(ρ, ζ) formulas.

(23)

The displacements must be constructed completely if the function χ₁(ζ) is to be found.

4 Obtaining and solving singular integral equation with regard to function χ₁(ζ)

One must use the condition (11) on the crack to find the function χ₁(ζ), with this aim the expression (23) should be transformed previously. The formula 1.445(1) from Ref. [30] is applied to calculate the sum of the series in the expression

The displacement’s expression (23) can be written in the form

The obtained expression has the unknown function χ₁(ρ). One needs to differentiate last equality with regard of variable ζ, to take ζ=δ -0 and to equal this expression to zero to satisfy the condition (11). The known identity is used during the simplification^[29]. As a result one gets

The integral equation of the first kind is obtained after passing to ζ=δ -0 and equating the expression to zero

The extraction of the weakly convergence part leads to the equation

Changing the variables in the obtained formula is done

(24)

Here, is the discontinuous Weber-Schafheitlin’s integral, R_i(τ, η) (i=1, 2) are the regular functions

This equation is solved approximately with the help of the orthogonal polynomial method^[26]. This is possible because of the spectral correspondence existence (see A11.5)^[26]

where are the Jacobi polynomials. The solution of (24) is constructed as the series expansion for the Jacobi polynomials

(25)

The expansion (25) is substituted at the equation (24), in accordance with the classical scheme of the orthogonal polynomials method, the equation (24) is multiplied by the expression and is integrated with regard of the variable τ. The orthogonal property of the Jacobi’s polynomials is used

As a result, the infinite system of the linear integral equations is obtained

(26)

Here,

These integrals were calculated in the obvious forms (see Appendix B):

The regularity of the obtained system (26) can be proved by the method proposed in Ref. [26]. The coefficients x_m were taken in the form . The system (26) was transformed to the system with the symmetrical matrix of the coefficients

(27)

where

5 SIF calculation

It is needed to estimate the value of such important mechanical characteristic as the SIF near the crack. With this aim one uses the formula

SIF is written in the form

(28)

with regard to the formula obtained earlier for the displacement u₁(ρ, ζ) and the changing of variables. The summands having the finite limit when τ → 1 + 0 are omitted in (28). After substitution of the expression (25) into the formula for SIF one can rewrite SIF in the form

(29)

In the last integral τ > 1, the direct use of the spectral correspondence for its calculation is impossible. Its calculation and the estimation of its behavior when τ → 1 + 0 is shown in Appendix C, where the next formula is proved

With the help of this formula, SIF is finally found in the form

(30)

6 Numerical results

The concretization of the proposed method was done for the case of a two-layered medium n=2. The next correspondences are written,

on the base of the formulas constructed for the general case. The calculations were done for the cylinder with the ratio of the cylinder’s radius to the height α=a: h=1: 4. The external loading is distributed along the upper part of the lateral cylinder’s surface by the law

where P_k=β_k^-1 cosβ_kb.

The SIF calculation was done for different ratios of the external and internal layers shear modulus G₁: G₂=2: 1, G₁: G₂=1: 1, G₁: G₂=1: 2 (the curvature 1, 2, and 3 correspondently) and for the cylinder with the dimensionless geometrical cylinder’s parameters c₁=r₁: a, c′=c: a, δ=d: h and b=0.8.

In Fig. 2 the dependence of the value K=K_Ⅲ: P^* from the ratio of the crack’s radius to the radius of the first layer X=c′: c₁ is demonstrated. The crack is located at the middle of the cylinder’s height δ=0.5, both layers have equal depth c₁=0.5. In Fig. 3, the dependence of the value K from the height of the crack’s location δ with crack’s radius c′=0.5c₁ and c₁=0.5 is shown. The influence of the first layer’s radius c₁ when δ=0.5 and c′=0.5c₁ on the SIF value is investigated and represented in Fig. 4. For numerical results with three significant digits after the comma, it was sufficient to take six equations in the system (27). In the case when the radius of the crack was more than 90% of the first layer radius, the calculations lost stability.

7 Conclusions

The analyses of the numerical results allow us to conclude

(i) The SIF values increase, when the ratio of the first layer shear’s module to the second one increases.

(ii) The SIF values increase with the increase of the crack’s radius with respect to the first layer’s radius. The increase is greater when the ratio of the first layer’s shear modulus to the second is greater.

(iii) The maximal values of the SIF are fixed when the crack is located in the middle part of the cylinder. In the case where the crack is nearer to the upper or to the bottom bases of the cylinder, the SIF value is much less.

(iv) Increasing the first layer’s radius and decreasing the second layer’s width leads to an increase of the SIF value. This increase is greater when the ratio of G₁ to G₂ is increased.

Appendix A

Applying integral transformation to torsion equation by generalized scheme

Let us consider the case when j=1. The crack is located at ζ=δ. The unknown function u₁(ρ, ζ) is a discontinuous function when ζ=δ, and it satisfies the equation

The interval of the integration (0; 1) will be formed as the sum of the two intervals (0; δ₁ -0) and (δ₁ + 0; 1). The first two summands in the equation are transformed in the standard way

In the last summand the integration by the parts will be realized with the intervals (0; δ -0) and (δ + 0; 1) As a result the jumps of the function and its derivative when ζ=δ will appear

Here, the conditions (4) were used. The application of the integral transformation leads to the equation in the transformation’s domain

Appendix B

Special method of integral’s calculation

This method is applicable to the integrals of the Mellin conjugation’s type^[31]

After application of Mellin’s integral transformation one gets

where K_j^*(s) are the Mellin’s transformations of the functions K_j(x), j=1, 2. The integrals are found with the help of the inverse integral transformation and the residual theorem.

This approach can be used, when the interval of the integration is (0; 1). One needs to continue the function K₁(x) by zero on the interval (1; ∞) in this case. Let us consider the integral

It has the form

after the changing the variables

Let us take , where and K₂(x)=.

The Mellin’s transformations of these functions are equal to and correspondently.

Then,

(31)

According to the inverse transformation’s formula

This integral is calculated with the help of the residual theorem, after closing of the contour to the right half-plane and with regard to the simple poles

One obtains after returning to y variable,

It can be written finally using the well known expansion of Bessel’s function into the series ^[29],

For the second integral’s calculation

(32)

one should use another well known equality connecting the modified Bessel’s function with Bessel’s function I₁(x)=-iJ₁(ix)^[28]:

Appendix C

Calculation of asymptotical formula

The approach used in Appendix B, is applied to the integral of the more general form

where

One gets the integral after changing the variables ,

where

The Mellin’s transformation of the function K₁(ξ) is^[31]

The Mellin’s transformation of the function K₂(ξ) is found with the help of formula 6.561(14)^[30]

The Mellin’s transformation K^*(s)=K₁^*(s) K₂^*(s) is found also. With the help of the inverse formula it can be written

This integral is calculated with the help of the residual theorem after closing the integration contour in the right half plane. One obtains with regard to the simple poles ,

where

The final formula is constructed after using the hypergeometrical Gauss’s function

It is needed to find the integral for the SIF’s calculation with the formula (29). Let us use the identity , where operator . Then, we get

Let us take γ=λ=1, ν=0, β= in the constructed formula

The differential operator D is applied to this result with the help of the formula I.1.18^[30],

when . As a result, the correspondence is found

Further use of the hypergeometrical Gauss’s function analytical continuing in the neighborhood of τ=1 + 0 with regard to 9.131(2)^[29],

leads to the asymptotical formula

Compliance with ethical standards Conflict of interest: the authors declare that they have no conflict of interest.

Acknowledgements The authors would like to express their sincere appreciation to S.P. DYKE for his understanding and suggestions in editing of this article.