Determination of the Stress State of a Piecewise Homogeneous Elastic Body with a Row of Cracks on an Interface Surface Subject to Antiplane Strains with Inclusions at the Tips

The stress state of a bimaterial elastic body that has a row of cracks on an interface surface is considered. It is subjected to antiplane deformations by uniformly distributed shear forces acting on the horizontal sides of the body. The governing equations of the problem, the stress intensity factors, the deformation of the crack edges, and the shear stresses are derived. The solution of the problem via the Fourier sine series is reduced to the determination of a singular integral equation (SIE) and consequently to a system of linear equations. In the end, the problem is solved in special cases with inclusions. The results of this paper and the previously published results show that the used approach based on the Gauss-Chebyshev quadrature method can be considered as a generalized procedure to solve the collinear crack problems in mode I, II, or III loadings.


Introduction
Stress analysis near a fracture in an elastic material is one of the most explored topics in solid mechanics.Calculation and stress analysis of engineering structures, particularly their connections and the determination of the stress and strain distribution fields of cracked bodies, have received attention from numerous investigators in recent years.The stress intensity factors, shear stresses, and crack opening displacements are major concepts that must be determined.The stress intensity factor is an important parameter that denotes the magnitude of the stress singularity.The singular order is a single real value, for example, 0.5 for a crack in a homogeneous material.The singular order of general interface corners may be real or complex.An asymptotic stress near the tip of a sharp interfacial corner is generally singular as a result of a mismatch between the materials' elastic constants.
In this paper, the problem of a piecewise homogeneous rectangular prismatic elastic body in an antiplane strain state due to antiplane forces is discussed.The plate is composed of two bonded dissimilar materials with a number of arbitrary collinear cracks on their interface surface.The aim of this work is the derivation of analytical expressions for the stress intensity factors of the cracks and the presentation of a new mathematical-numerical approach to solve singular integral equations related to the beginning of stresses at the tips of cracks in antiplane deformations; thus, a numerical method to calculate the SIFs of an interface crack between dissimilar materials is developed.Calculation of the SIFs for bimaterial plates in dimensionless form is conducted for several cases: a row of cracks, one and two cracks, and one crack with inclusions at the tips.The objectives of the present study are (i) to present a new method for calculating SIFs of interfacial cracks subject to antiplane loading and (ii) to investigate the influence of inclusion moduli on SIFs to reduce SIFs in cracked bodies and to direct us towards a method for repairs.
The existence of three-dimensional effects at cracks has been known for many years, but understanding has been limited, and for some situations, it still is.Understanding improved when the existence of corner point singularities and their implications became known for straight through-thethickness cracks [1,2].
It has been known for a long time that shear and antiplane fracture modes are coupled.This means that shear or antiplane loading of an elastic plate with a throughthe-thickness crack generates a coupled three-dimensional antiplane or shear singular stress state, respectively.These singular stress states (or coupled fracture modes) are currently largely ignored in theoretical and experimental investigations as well as in standards and failure assessment codes of structural components, in which it is implicitly assumed that the intensities of these modes as well as other three-dimensional effects are negligible in comparison with the stress fields generated by the primary modes (modes I, II, and III) [3,4].
The theoretical bases of fracture stresses are discussed in the literature [5][6][7], and the conclusions of numerous studies and investigations on the derivation of SIFs are categorized in [8][9][10].Most of these studies were performed on homogeneous plates.It is proven that, under certain circumstances, the three-dimensional governing equations of elasticity can be reduced to a system where a biharmonic equation and a harmonic equation have to be simultaneously satisfied.The former provides the solution of the corresponding plane problem, while the latter provides the solution of the corresponding out-of-plane shear problem [11].On the other hand, a mixed fracture mode under antiplane loading may also occur.This coupled fracture mode represents one of threedimensional phenomena that are currently largely ignored in numerical simulations and failure assessments of structural components weakened by cracks.It arises due to the boundary conditions on the plate-free surfaces, which negate the transverse shear stress components corresponding to classical mode III.Instead, a new singular stress state in addition to the well-known 3D corner singularity is generated.This singular stress state can affect or contribute significantly to the fracture initiation conditions [12,13].
Inclusions and cavities are also important in understanding the mechanical behavior of structures and are studied in several papers, for example, 2D linear elastic materials [14] and antiplane shear cracks [15].Photoelasticity and finite element methods have also been employed to study the interaction between collinear cracks, and good agreement was found [16].Photoelasticity is very helpful in investigating the stress state near inclusions.The results show that the singular stress field predicted by the linear elastic solution for an inclusion can be generated in reality with great accuracy [17].The inclusions form a thin material that constituted a rigid line inclusion, embedded in a linear elastic body to produce an inhomogeneous stress state.The experiments fully validate the stress state calculated for an elastic plate [18].
The mechanical behavior of thin inclusions is fundamental to the design of composite materials.It is realized that, for a given geometry and boundary condition,  III depends on the gradation of both the modulus of elasticity and Poisson's ratio [19].The cracked sandwich plate twist specimen is viable to characterize mode III fracture [20].
The stress intensity factors can be calculated by a pathindependent h-integral and through the virtual crack closureintegral method (VCCM) for numerical implementation [21].The present study is aimed at investigating the stress state of a piecewise homogeneous elastic body which has a row of collinear cracks in mode III.The numerical procedure based on the loading Gauss-Chebyshev quadrature method is applied.This approach can be used for other multicrack problems or more complicated types of loading [22,23].
At the end of this paper, we discuss the influence of the bimaterial nature of the body, the distance and geometry of the cracks, and the presence of inclusions at the tips on the characteristics of antiplane shear stresses and deformations.

Derivation of the Singular Integral Equation for the General Form of the Problem
A piecewise homogeneous rectangular body in the Cartesian coordinate system Oxyz is considered as shown in Figure 1.
In the equations below, subscript 1 denotes rectangular plate  1 and subscript 2 denotes  2 .

󵄨 󵄨 󵄨 󵄨 󵄨𝑦=ℎ
=−ℎ 2 +0 =  2 () where the notation  ()  ( = 1, 2) for shearing stress components in plates   is used.Because the pair of vertical sides of the rectangle is rigidly fixed, the boundary value problem for the complex potential of state can be transformed to the problem in the half-plane, and its solution can be obtained in quadrature using conformal mapping [24,25].
Through the above assumptions, the body is under the antiplane strain condition, so the displacements of the crack edges are along the Oz-axis with the base on the Oxy surface.
It is required to calculate the dislocation density of the crack boundaries, the related stress intensity factors, and the distribution of shearing stresses in the regions that are outside of the crack system, with total length   = [0, ] \ .
For the development of the above-mentioned functions, the plate is divided into an upper rectangle ( 1 ) and a lower rectangle ( 2 ) on the Ox-axis with 0 ≤  ≤ , and we introduce a function   (, ) ( = 1, 2) as the only nonzero displacement component along the Oz-axis for both   ( = 1, 2) rectangular plates and separately investigate their individual elastic equilibria.In this manner, for rectangle ( 1 ), according to Hooke's law, the below boundary value problem is found: Suppose the notation in which the function  (0) + (), as mentioned before, is the shear stress loading on the cracks' top boundaries  + and () is the unknown tearing shear stress that acts on the noncracked interface surface system   .
According to [26], to solve (4) by means of the Fourier sine series, it is obtained that Hence from the inverse Fourier series we deduce Equation ( 7) is the Fourier sine series with coefficients given by (6).The reason for using the Fourier sine series is that boundary value problem (4) with the Laplacian and Dirichlet boundary conditions can be solved by separation of variables, and because the boundary conditions in the  variable are the Dirichlet boundary conditions, this generates a Fourier sine series with eigenfunctions sin(/) and eigenvalues (/) 2 ,  = 1, 2, . ... Considering ( 6) and (7), it follows that by taking the two parts of the differential equation ( 4), multiplying by sin(/), integrating it over the interval (0, ), and taking into consideration the boundary condition (4), we obtain the below boundary value problem: with Boundary value problem (8) can be solved by the following formula: For ( = 1, 2, . ..), we can conclude Using a similar approach, the boundary value problem related to the rectangular plate ( 2 ) is obtained: Then, considering the  1 (, 0) and  2 (, 0) displacements and the following functions: where Φ  , Ψ  , Ω  , and   are the Fourier sine series coefficients of the above functions that are similar to those introduced in ( 7) and ( 9), we obtained Substituting coefficients ( 14) into ( 11) and (12), a system of two linear equations with respect to Ω  and Ψ  is obtained.Solving the system, coefficients Ω  and Ψ  can be expressed by Φ  and   and, in turn, all of the functions in (13) can be expressed by the Φ  and   coefficients.In this case, as we expect, the required function (13) after some simple transforms is found in dimensionless form as follows (see Appendix): where the first part in the first integral when  =  is the main quantity in the Cauchy formulation.
The other parameters are and Δ  (, ℎ + , ℎ − ) =  sh(ℎ + )ch(ℎ − ) + ch(ℎ + )sh(ℎ − ) and the corresponding dimensionless parameters are In this manner the system of lengths () converts to system ( 0 ), and, through function ( 13), the dislocation density is imposed at the boundaries of the crack edges: Then, by applying key equation (15) to the system ( 0 ) and introducing the new variables regarding the unknown dislocation density  0 () from ( 18), the singular integral equation governing the problem is obtained as follows: Based on ( 16), the following influencing factors are derived: where   () and  −1 () are Chebyshev polynomials of the first kind and the second kind, respectively.The singular integral equations ( 20) and ( 21) in a condition that explains the continuity of displacements at the crack tips are expressed as which, based on the previously discussed dimensionless variables, converts to Hence after solving the singular integral equations ( 20) and ( 21) by considering condition (23) and verifying key equation (15) in the region outside of the crack system   0 , the antiplane shear stress () formula in the dimensionless form is obtained.

Derivation of the Linear System of Equations
To solve the singular integral equations ( 20) and ( 21) under condition (23), the numerical solution for singular integral equations [7,27,28], which is based on the Gauss quadratic solutions for ordinary and singular integrals, is used.First, each section of [  ,   ] ( = 1, ) from system Λ 0 , by means of the technique of transitioning to new variables , (−1 ≤ ,  ≤ 1) can be converted to the section [−1, 1].The singular integral equations ( 20) and ( 21) can therefore be written as follows: The function () and kernel functions  0 (, ) and (, ) are expressed by relations (20) and (21).In this case, condition equation ( 23) may be expressed in the following form: Using the above approach, the singular integral equation ( 25) through condition equation ( 26) is reduced to a system of linear algebraic equations: In ( 25), (26), it is supposed that where  is an arbitrary natural number.Consider are the roots of the Chebyshev polynomials of the second kind  −1 () and of the first kind   ().
The opening displacement of each crack   = [  ,   ] is given by the equations: which can be expressed by means of the dimensionless variables as where It must be noted that the system of equations ( 27) can be solved by the Gaussian technique, and (31) also can be solved with high accuracy by means of the Gauss quadratic formula using (28).
At the end, the stress intensity factors, SIFs, for the antiplane mode of fracture at the crack tips [  ,   ] are obtained according to [24] ( = 1, ): which can be converted to the dimensionless form and be expressed through   () using formula (28).
Because of the large number of variable parameters involved in the problem, only three different cases are considered.
Case A: a plate with one symmetric crack: the crack length is assumed arbitrary, and antiplane shear loading is two point loads.The stress and dislocation states are studied for different crack sizes and shear modulus ratios.
Case B: a plate with two cracks: the cracks are located symmetrically and have arbitrary length and distance.The stress and dislocation states are studied for different crack sizes, distances between adjacent crack tips, and shear modulus ratios.
Case C: it is a plate with thin-walled inclusions at the crack tips, in which the sensitivity of SIFs to the size of the inclusion domain with respect to crack length and the effect of the shear modulus ratio of the inclusion material to the plate on the stress state are studied.

Case A: Plate with One Crack
A rectangular plate in the Cartesian coordinate system Oxyz is considered, which is composed of an upper rectangle  1 with modulus of rigidity  1 , length , and height ℎ and a lower rectangle  2 with modulus of rigidity  2 , length the same as , and height ℎ.On the interface surface line of the two rectangular plates in the interval 0 ≤  ≤ , there is a central crack with length equal to 2 (Figure 2): The boundaries of the crack indeed have not any loadings.Furthermore, the vertical edges  = 0 and  =  of rectangles   ( = 1, 2) are restrained, and the upper and lower horizontal edges  = ±ℎ at the midpoints (/2, ±ℎ) are stressed by antiplane shear forces equal to   (1)        =ℎ−0 =  (2)        =−ℎ+0 =  ( −  2 ) (0 <  < ) (35) in which  ()  ( = 1, 2) are antiplane shear stresses on the top and bottom boundaries of plates   and () is the Dirac delta function.
It is necessary to determine the dislocation field around the crack boundaries, the stress intensity factors, SIFs, and the shearing stresses on the interface surface   = [0, ] \  outside of the crack region.
The singular integral equation (SIE) governing the current problem can be found in the previous section, particularly by observing the dislocation field around the crack boundaries as follows: in which   (, 0) are the displacement components of region { = ±0, 0 ≤  ≤ } of rectangles   ( = 1, 2) in the direction of Oz, so the singular integral equation of this boundary value problem from [29] is written as follows: where   () and  −1 () are Chebyshev polynomials of the first kind and the second kind, respectively.
To solve the singular integral equation (37), the Gaussian quadratic solution for ordinary and singular Cauchy integrals can be used.By substitution, Using the approach shown here and also in [7,27,28], solving the singular integral equation (37) leads to a system of linear algebraic equations as follows: where  is an arbitrary natural number and   and   are the roots of Chebyshev polynomials of the second kind  −1 () and the first kind   (), respectively.We can substitute the below parameters to express the system of equations (39) in a simpler form:  (41) By assumption of the below expression for the right hand side of ( 41) (0)  ( = 1, ) is found.The solution of system ( 26) is thereby obtained as follows: For the derivation of SIF, the general formula from [7, 30] may be used; that is, The above formula can be converted to a dimensionless form as follows: (45)

Case B: Plate with Two Cracks
In this section, a rectangular plate in the Cartesian coordinate system Oxyz with two cracks is considered.On the interface surface of the two segments in the interval 0 ≤  ≤  there are two central cracks located symmetrically at  = ∪ 2 ≃1 [  ,   ] that have equal lengths (Figure 3): The boundaries of the cracks have no traction; furthermore, the vertical edges of plates   ( = 1, 2) at  = 0 and  =  are clamped, and the upper and lower horizontal edges  = ±ℎ are loaded by antiplane distributed shear loading (), so that  (1)        =ℎ−0 =  (2)        =−ℎ+0 =  () (0 <  < ) , in which  ()  ( = 1, 2) are the antiplane shear stresses on the top and bottom boundaries of segment   .
The displacement of the crack boundaries , the stress intensity factors SIFs, and the shear stresses at the interface surface   = [0, ] \  outside of the cracks are determined.Moreover, it is necessary to investigate the influence of adjacent crack tips  2 and  1 .
Suppose that the upper and lower edges of the plate are stressed by antiplane shear forces  so that () = ( − /2) in which () is the known Dirac delta function.By this assumption, the function () is calculated as follows: To calculate the shear stress, making use of the variables shown before and variable  and taking into consideration (38), the shear stress is concluded: From the symmetry of the problem due to axis  = /2 and taking into consideration only the right hand crack with end points at  =  2 and  =  2 , the stress intensity factors are defined as follows [3,4,7,26]: where () are the shear stresses according to (31).

Case C: Plate with Inclusions at the Crack Tips
In this section, attention is specifically paid to the effect of inclusions at the crack tips on the stress intensity factors (SIFs).
A prismatic elastic body Ω with a rectangular cross section in Cartesian coordinates Oxyz occupying an area Ω = {0 ≤  ≤ ℓ; −ℎ ≤  ≤ ℎ; −∞ <  < ∞} and possessing a shear modulus  is considered.The elastic body is rigidly clamped at  = 0 and  =  and loaded by shear forces equal to () acting both in the positive and in negative directions of the Oz-axis by the horizontal  = ±ℎ.Furthermore, on the symmetry plane  = 0, it has a through-the-thickness crack in the shape of the strip  = { = 0; ℓ/2 −  <  < ℓ/2 + ; −∞ <  < ∞}, ( < ℓ/2).Shear forces of equal intensities  0 () are acting in opposite directions along the Oz-axis on the upper (+) and lower (−) areas of the crack edges of the crack.Additionally, at the tips the edges of the crack are joined by thin-walled inclusions with shear modulus  0 deforming by the Winkler model (1867) that act as linear elastic springs:  =  (Figure 4) [31].
It is necessary to determine the dislocation density on the crack edges, the SIFs, the shear contact stresses on the edges of the inclusion, and the shear stresses outside the crack along the surface of its location.This is formally similar to the Dugdale-Barenblatt model for a central crack containing yielding as confined and localized narrow plastic zones, which shows the effect of yielding on the crack length [32,33], but, in this problem, which is a linear elastic fracture mechanics problem, we investigate the effect of inclusions at the end areas of crack on the mechanical behavior of crack tips, decreasing the dislocation density and antiplane SIFs.
The component   = (, ) in the direction of the Ox-axis is the only nonzero component of the displacement, and   and   are the only nonzero stress components.Therefore, the problem can be mathematically stated as a boundary value problem in the following way: ( Let us again consider the following functions based on [18,34]: where Φ  , Ψ  , Ω  , and   are sine Fourier series coefficients as before.Inserting these functions into the displacements and stresses on the edges of the crack at the interval [0, ] of rectangle  2 then according to (4), the following equations are obtained: After some simple transformations and calculations, the following equations can be derived: Transforming the above equation, the following equation is derived: As mentioned above, the determinative singular integral equation (SIE) can be reduced to a system of linear equations.For this purpose, applying the Heaviside function this equation can be expressed in the following way: with the notations explained above.Supposing the determinative SIE (39) with conditions of ( 25) and ( 40) is reduced to the system of linear equations For this case, the crack edges are free of shearing forces, and concentrated shearing forces are acting on the horizontal sides of the rectangular plate; that is, where () is a certain Dirac delta function.Additionally, the following equation can be obtained: Taking into consideration the above-mentioned definition, the function f() can be expressed in the following way: and the function () from (57) is obtained: It is obvious that the function  0 () with respect to the symmetry of line  = /2 in this special case and, consequently, the function  0 () according to (61) are odd functions.Therefore, the components of the second integrals in (56) and (58) containing polynomial  or  in arguments tend to zero, so that the above-mentioned equations and the kernel-matrix   of the system of equations (62) are simplified.
It must be emphasized that the expressions of functions f() and () from ( 65), (66) on various intervals are used in the equations.
The numerical analysis of the main characteristics of the stated problem can be carried out for the considered special case.

Numerical Results and Discussions
To solve the system of equations (40) that is summarized in system (41) with the parameters shown in (42), The dimensionless SIF  (0) III according to formula (50) using the vector   ,  = 10 and considering the assumptions  = /3, /4, /8, /16, and /32 for different crack lengths is calculated.The results are shown in Table 1, indicating that longer cracks lead to larger SIFs.
The shear stress  0 () from (49) in dependence of  and  taking the rigidity ratios of the two materials  =  1 / 2 = 0.1, 0.3, 0.5, 0.8, 1.0, 2, 5, 10 is obtained and the graphs are shown in Figure 5.The dimensionless crack dislocations that depend on parameters  and  are calculated, and the curves are represented in Figure 6.This figure shows that the maximum dislocation of a crack is at midpoint of its length, and it has a zero value at the crack tips, the points where maximum shear occurs.
Parameters  and  are supposed to be equal to  = /3 and  = /4; /8; /16; /32; /64 in order to determine the influence of the distance between the two adjacent cracks on their fracture characteristics.The shear stress  0 () is calculated according to (49), in which function () is obtained  from (37), supposing the quantity  0 = 0.01.Making use of the above calculated stresses, the variation curves of  0 () can be drawn for  = 0.1, 0.3, 0.5, and 1.0, as shown in Figure 7.
The crack dislocation density as a function of the edges' opening displacements Ψ (0)  () may be found through formula (31).The Ψ (0)  () curves are shown in Figure 8 again for varying  and .At the end, using (50), the stress intensity factors in the dimensionless form can easily be calculated, as shown in Figure 9 and Tables 2 and 3, for the above parameters.The variation of SIF  III based on parameter  and  = 0.3, 0.5, and 1.0 is shown in Figure 10 to recognize the state of the SIFs under the change of the lengths of cracks.
For the case of  = 1.0, the cracks approaching each other (that means a smaller value of ) lead us to construct curves as shown in Figure 11 for  0 () and Figure 12 for Ψ (0)  ().The main characteristics of the crack problem with inclusions in the mentioned special case can be calculated through numerical analysis.For calculation it is considered as some    engineering practical ratios for the shear modulus quantities  =  0 / that are  = 0.0 (crack without repair),  = 0.5, 1.0 (repair with the same material), and  = 2.0, 5.0, and 10.0 (for very high rigid repair materials).By solving system of linear equations (62) which leads to solving a matrix (10 × 10) and obtaining vectors   (10 × 1), the shear stress  0 () can be calculated from (58) and the results are shown in Figure 13.
Dislocation functions ψ0 () = 2 0 () were investigated using relation (31), as shown in Figure 14 in which the maximum displacement of crack boundaries occurred at midpoint and obviously with zero value at crack tips.
The stress intensity factors SIFs  (0) III calculated through (50) are shown in Figure 17, which presents the decreasing trend of the  III curve when the ratio of  =  0 / increases.It also shows that the crack tip repairing by this method can reduce the SIF by approximately 50 percent, which means that this approach is very effective to control the crack propagations in cracked plates and it is also a treatment for singularities at the crack tips, defects, and holes.Figures 15  and 16 show that the reduction of the magnitude  III is not  high when we use a very rigid material for tip repair.For example, for  = 2.0 and  = 10.0 (ratio = 5), the reduction of  III is only 18 percent (1.98/2.42).
The above results coincide with experimental observations of high concentrations and singularities in the stress fields within elastic materials [17,18,35,36].

Conclusion
In the present paper, the singular integral equations governing the piecewise homogeneous elastic plate problem subject to uniform remote antiplane shear loading have been considered.To determine the antiplane shear stresses, the crack edge dislocation densities, and the mode III stress intensity factors, a new mathematical-numerical calculation  is developed.To perform elastic analyses and investigations on the state of stresses at the crack tips, the dimensions of the cracks, the distances between adjacent tips, and the influence of shear modulus ratios between the two materials have been varied and studied.The governing singular integral equation of a problem of stress-strain state of a piecewise homogeneous elastic prismatic body of a rectangular cross section, when there is a system of an arbitrary finite number of collinear cracks on the interface line of dissimilar materials, is obtained.For problem solving, the method of mechanical quadrature is used.It allowed conducting a detailed analysis of the main mechanical characteristics of the stated problem and revealed their dependence on the mechanical and geometrical parameters of the problem.For the case of cracks coming together it has been shown that the stress intensity factors  III grow based on two parameters, the total distance between their far tips and the closeness of their near tips.The technique presented in this paper can be used to solve a class of problems associated with cracking in bimaterial interfaces.For a crack with inclusions at the tips, it is observed that the insertion of a material at the crack tips avoids singularities, reduces the antiplane SIF  III , and strengthens and controls the crack propagation in the region at the tips.This is an effective method and does not require using a material with a very high shear rigidity value  0 .Meanwhile, the crack opening displacement (COD) and the shear stresses at the crack tips decrease.This suggests a method to repair cracked plates and members.
Practice shows that the method of mechanical quadrature is a very effective method and results in a very good convergence.The data of Table 2 confirms this.The accuracy of the method in predicting the intensity factor may be verified by a comparison with experimental measurements, carried out by a photoelasticity method, and by commercial finite element software.The drawn conclusions provide meaningful reference for the analysis of SIFs in mode III.Numerical results show the influence of ratios of shear moduli on the stress state.
The results of this paper and the previously published results show that the used approach based on the Gauss-Chebyshev quadrature method can be considered as a generalized procedure to solve the collinear crack problems in mode I, II, or III loadings.

Figure 1 :
Figure 1: Rectangular piecewise homogeneous elastic body with several collinear cracks.

Figure 4 :
Figure 4: Plate with inclusions at crack tips.

Table 1 :
Dimensionless SIF  (0) III as a function of  and .