Problems of interaction longitudinal shear waves with V-shape tunnels defect

The problem of determining the two-dimensional dynamic stress state near a tunnel defect of V-shaped cross-section is solved. The defect is located in an infinite elastic medium, where harmonic longitudinal shear waves are propagating. The initial problem is reduced to a system of two singular integral or integro-differential equations with fixed singularities. A numerical method for solving these systems with regard to the true asymptotics of the unknown functions is developed.


Introduction
At present, there exist numerous solutions of two-dimensional dynamical problems of the theory of elasticity for bodies with defects shaped as a segment of a straight line or an arc of a smooth curve. Examples of solving such problems can be found in [1][2][3][4][5][6][7][8]. However, the real defects can have corner points, be piecewise smooth, intersect one another, and bifurcate. Practically, the problems of determining the dynamical stress state in a neighborhood of such defects have not yet been solved. This is related to the difficulties arising when they are solved by using the method of boundary integral equations widely applied nowadays and the problems are reduced to singular integro-differential or hypersingular equations with fixed singularities. Most of the studies deal with static problems for bodies whose defects have corner points. First, we mention [9][10][11][12], where the exact solution was obtained by the Wiener-Hopf method, and an accurate value of the stress intensity factor was determined. These solutions and the results of [12] indicate that the presence of kernels with fixed singularities influences the singularity of the solutions in neighborhoods of the ends of integration intervals. The stress state near branched, broken, and edge cracks was also determined in [14][15][16][17][18][19], where the integral equations were numerically solved by the method of mechanical quadratures. This method is based on the application of the Gauss-Chebyshev quadrature formulas that give a root singularity of the solutions. In this case, the true asymptotics of the solution is not taken into account, or an additional condition leading to a singularity weaker than the root one is imposed on the solution. Another drawback of the numerical methods used in those works is the formal application of the Gauss-Chebyshev quadrature formulas to integrals with fixed singularities. As a result, the convergence is rather slow (in order to obtain results with an error less than 0.1%, several dozens of collocation points should be used). In the present paper, we solve the problem of determining the dynamical stress intensity factors for tunnel defects (cracks or thin rigid inclusions) which are V-shaped in a cross-section under the interaction with a harmonic longitudinal shear wave. The problem is reduced to a system of three singular integro-differential or integral equations solved by the method of collocations. Within this method, the true singularity of the solution is considered, and special quadrature formulas are used to calculate the integrals with fixed singularities.

Statement of the problem and its reduction to a system of singular integral equations
We consider the unbounded isotropic elastic medium under the conditions of antiplane deformation with a through V-shaped defect in a cross-section in the plane Oxy (figure 1). It may be a crack or a thin rigid inclusion. The cross-section of the defect occupies two segments what are of length 2d k and make angles α k , k = 1, 2, with the axis Ox. The defect interacts with a plane longitudinal shear wave. Its front makes an angle θ 0 with the axis Ox and causes the displacements in the medium along the axis Oz. In (1), the symbols G and ρ stand for the shear modulus and the medium density, and ω is the frequency of oscillations. The time dependence is given by the factor exp(−iωt) that is omitted here and below. Let the vector of displacements caused by the wave reflected from the defect have a single component W (x, y) which is nonzero under an antiplane deformation. In the coordinate system Oxy, it satisfies the Helmholtz equation To specify the boundary conditions on the defects, we associate the coordinate system with each of their segments. Let these systems O k x k y k , k = 1, 2 be centered at the middles of segments, and let its centers have coordinates O k (d k cos α k , d k sin α k ). Let be the displacement in the coordinate systems related to the kth segment of the defect. In the case of a thin rigid inclusion, in view of its small thickness, we formulate the corresponding boundary conditions on its median surface. Assume that the condition of perfect contact is realized between the body and the inclusion where On the surface of the inclusion, the tangential stresses are discontinues and their jumps are denoted by In formula (4), c is the amplitude of longitudinal oscillations of the inclusion in the equation of its motion as a solid body. In the case of harmonic oscillations, this equation becomes where h is the inclusion thickness and ρ 0 is the inclusion density.
In the case of a crack, its surface is assumed to be free form stresses, In addition, the displacement on the crack surface has a jump denoted by Also, the continuity of displacements along the edge yields We start solving the posed problem by constructing a discontinuous solution of the Helmholtz equation in the coordinate system O k x k y k , k = 1, 2, for each link of the defect [6]. In the case of an inclusion, this is a solution with jumps (5): In the case of a crack, we construct the discontinuous solution with jumps (8): In formulas (9), (10), we denote where H (1) 0 (z) is the Hankel function. Then the displacements of the reflected wave field can be given in the form where In order to finally determine the displacements of the diffraction field, we should find the unknown jumps (5) or (8). To this end, we use conditions (4) or (7). First, we consider the case of inclusion. To determinate jumps (5) from the boundary conditions (4), we obtain a system of singular integral equations. For the singular component of their kernels contain a singularity in the form of the Cauchy kernel, instead of (4), we should use [7] the conditions The result of substitution of (3), (11), and (9) into (12) is the following singular integral equation with an auxiliary condition: In system (13), we introduce the notation In the case of a crack, the following system of singular integro-differential equations for jumps (8) and their derivatives is obtained after substitution of (3), (10) and (11) into (7): We used the following notation in system (15): As is seen in (14) and (16), the functions q kl (τ, ξ) have singularities at τ = ±1, ξ = ±1. The functions R lk , U lk , R lk (τ, ξ), D l (τ ), k, l = 1, 2, are bounded for −1 ≤ τ , ξ ≤ 1.

Solution of the problem in the case of an inclusion
The presence of the fixed singularity τ = ±1, ξ = ±1 in the singular component of the system integral equations (13) affects the asymptotics of its solutions in neighborhoods of the points τ = ±1. The singularity of the solutions in neighborhoods of these points can be determined by analyzing the asymptotic properties of singular integrals or by investigating the symbol of the singular kernel [13]. As a result, we established that the unknown function is to be sought in the form where ψ k (τ ), k = 1, 2, satisfy the Hölder condition, and the exponents are then we obtain ψ 0k (±1) = 0. This is why we can assume that where g k (τ ) are new unknown functions. Substituting (19) and (18) into (17), we obtain the following representation of the solutions of integrals equation (13): Here σ 1 = λ 2 = 1 − δ, σ 2 = λ 1 = 1/2. Further, the approximate method of solution is based on approximation of the functions g k (τ ) by the interpolation polynomials of degree (n − 1): where Q n1 (τ ) = P 1−δ,1/2 n (τ ), Q n2 (τ ) = P 1/2,1−δ n (τ ) are the Jacobi polynomials and τ m1 , τ m2 are roots of these polynomials. Then, for the singular integral with Cauchy kernel, the following quadrature formula holds [20]: where ξ j1 , ξ j2 (j = 1, . . . , n + 1) are roots of the Jacobi functions of the second kind J 1/2,1−δ n (ξ), J 1−δ,1/2 n (ξ), and A mk are the coefficients of the corresponding Gauss-Jacobi quadrature formula. The use of representation (20) requires to calculate the following integrals with Cauchy kernel: We obtain their values by the method described in [21] and based on their transformation to the Mellin convolution. Further, applying the theorem on convolution, we can represent these integrals as the sum of residuals at the poles of the integrands. Then we calculate the integrals with a fixed singularity Let 0 < r 1 < 1 be a positive number. In the case where 1 ± ξ k > r 1 , the integral (24) is not singular, and it can be calculated by the Gauss-Jacobi quadrature formulas with the corresponding weight function [22]. If 1±ξ k → 0, then we should replace ϕ k (τ ) by expression (20) and use a method based on using the Mellin integral transformation [21]. We represent these integrals in terms of the following power series which converge for 0 ≤ 1 ± ξ k < r 1 < 1: The integrals are computed by a method similar to that used to calculate integrals (23).
To calculate the integrals with a logarithmic singularity, we use representation (20) and approximate the function g k (τ ) by the interpolation polynomial (21), which we preliminarily transform according to the Darboux-Christoffel identity [22]. Thereafter we obtain the formula where Here Ψ(x) is the logarithmic derivative for the Γ function.  The obtained formulas for the singular integrals (22), (25), (26) and Gauss-Jacobi quadrature formulas applied to the integrals with regular kernels allow us to replace (13) and (6) by a system of linear algebraic equations. As a result of solving this system, we obtain g m1 = g(τ m1 ), g m2 = g(τ m2 ), ψ 1 (±1), ψ 2 (±1), c 0 . Then the approximate solution of the system is determined by formulas (20) and (21). One of the important characteristics of the stress state near the inclusion is the stress intensity factor (SIF). It is determined from the known asymptotic representation [23] of stresses in neighborhoods of the end of the inclusion. We obtain the following approximate formulas for the dimensionless SIF: On the basis of the obtained formulas, the dependence of SIF on the dimensionless wave number κ 0 = κ 2 d is investigated. For this purpose, the inclusion with links of equal length d symmetrical with respect to the axis Oy was considered (figure 2). The links of inclusion come from the origin of coordinates and make an angle 2β between them. Realizing the proposed method of solution, we first perform numerical investigation of its practical convergence. The calculation was carried out for the following data: Under these conditions, the equality K 1 = K 2 = K is executed. As results of computation by formula (27), we construct graphs of the dependence of the absolute value of the dimensionless SIF k = K/(G √ d) on the dimensionless wave number, which is shown in figure 3. The curves correspond to the number of interpolation nodes n = 5, 10, 15, 20 in approximation (21). We can see that, in the region of low frequencies κ 0 < 2, to attain a fairly high accuracy (the relative error does not exceed 1%), it suffices to use n = 5 interpolation nodes. In the whole frequency range 0 ≤ κ 0 ≤ 5, to obtain values of SIF with an error below 0.1%, it suffices to have n = 20 interpolation nodes in formula (21). The Let us use quadrature formulas (29), (31), (34), (35), the Gauss-Jacobi quadrature formulas, and the roots of Jacobi polynomials as points of collocation. Then system (15) can be reduced to a system of linear algebraic equations for g mk = g k (τ mk ). By solving this system, we obtain the following formulas for the approximate values of SIF: .
The numerical realization of the proposed method is carried out for the crack shown in figure 3. The calculation was carried out for the same data as in the case of inclusion. The results of studying the practical convergence of the proposed method are given in figure 6 (21). These results show that, in the region of low frequencies κ 0 < 1, to attain a fairly high accuracy (the relative error does not exceed 1%), it suffices to use n = 5 interpolation nodes. In the whole considered frequency range 0 ≤ κ 0 ≤ 6, to obtain values of SIF with an error below 0.1%, it suffices to have n = 20 interpolation nodes in formula (21). We also numerically studied the influence of the angle β between the crack links on the dependence of SIF on the frequency. The results of calculations are given in figure 7. The curves correspond to values: 1 -β = 15 • , 2 -30 • , 3 -45 • , 4 -60 • , 5 -87 • and are constructed at the angle of wave propagation θ 0 = 90 • . We can see that, at low frequencies 0 ≤ κ 0 ≤ 2, as the angle increases, the value of SIF also increases. The largest values of SIF are observed as the crack is transformed in a segment of the straight line β ≈ 90 • . For the further increase in the frequency 2 ≤ κ 0 ≤ 6, the dependence becomes more complicated and has many maxima and minima. This complication is related to the multiple reflections of waves from the sides of the angle made by the cracks links.

Conclusions
In the presence of angular points in cracks and thin inclusions the method of integral equations reduce boundary value problems to singular integral or integro-differential equations with fixed singularities. For the numerical solution of such equations one must take into account the real property of unknown functions and apply special quadrature formulas for singular integrals. It ensures rapid convergence of numerical method of solution. On the example of V-shaped defects it can be seen the shape and type of defect substantially influence on stress state in their region. In particular, the angle formed by the sides of the defect substantially influence on the SIF dependence on the frequency. So, in the low frequency region. In the case of crack, the largest values of the SIF are observed for strait lines cracks. If the defect is inclusion for this angle values of the SIF will be the smallest. The presence of maximus resonance type in the region of high frequencies caused by the multiple reflection of the waves from the defects sides is shown.