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.
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-the-thickness 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 through-the-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–7], and the conclusions of numerous studies and investigations on the derivation of SIFs are categorized in [8–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 . On the other hand, a mixed fracture mode under antiplane loading may also occur. This coupled fracture mode represents one of three-dimensional 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  and antiplane shear cracks . Photoelasticity and finite element methods have also been employed to study the interaction between collinear cracks, and good agreement was found . 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 . 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 .
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, depends on the gradation of both the modulus of elasticity and Poisson’s ratio . The cracked sandwich plate twist specimen is viable to characterize mode III fracture .
The stress intensity factors can be calculated by a path-independent h-integral and through the virtual crack closure-integral method (VCCM) for numerical implementation . 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.
2. 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 and subscript 2 denotes .
The top rectangular plate has a rigidity modulus , length , and height , and the bottom rectangular plate has a rigidity modulus , length , and height . On the bonding surface of the two materials in the interval , there are arbitrary collinear cracks with total length and individual lengths ; that is,The upper and lower boundaries of these cracks are subjected to antiplane shear stresses , which act as varying internal loadings of openings:
The edges and of rectangles () are restricted, and the edges of the rectangles and are subjected to antiplane shear stresses: (), where the notation 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 .
For the development of the above-mentioned functions, the plate is divided into an upper rectangle () and a lower rectangle () on the Ox-axis with , and we introduce a function () as the only nonzero displacement component along the Oz-axis for both () rectangular plates and separately investigate their individual elastic equilibria. In this manner, for rectangle (), according to Hooke’s law, the below boundary value problem is found:Suppose the notationin which the function , 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 .
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 and eigenvalues , .
Considering (6) and (7), it follows that by taking the two parts of the differential equation (4), multiplying by , integrating it over the interval , 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 , we can conclude
Using a similar approach, the boundary value problem related to the rectangular plate is obtained:
Then, considering the and 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 areand and the corresponding dimensionless parameters areIn this manner the system of lengths converts to system , 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 and introducing the new variablesregarding the unknown dislocation density from (18), the singular integral equation governing the problem is obtained as follows:
Based on (16), the following influencing factors are derived:where and 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 , the antiplane shear stress formula in the dimensionless form is obtained.
3. 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 from system , by means of the technique of transitioning to new variables can be converted to the section . The singular integral equations (20) and (21) can therefore be written as follows:
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 thatwhere is an arbitrary natural number. Consider are the roots of the Chebyshev polynomials of the second kind 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
At the end, the stress intensity factors, SIFs, for the antiplane mode of fracture at the crack tips are obtained according to  : 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 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 : 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 : 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.
4. Case A: Plate with One Crack
A rectangular plate in the Cartesian coordinate system Oxyz is considered, which is composed of an upper rectangle with modulus of rigidity , length , and height and a lower rectangle with modulus of rigidity , length the same as , and height . On the interface surface line of the two rectangular plates in the interval , there is a central crack with length equal to (Figure 2):
The boundaries of the crack indeed have not any loadings. Furthermore, the vertical edges and of rectangles () are restrained, and the upper and lower horizontal edges at the midpoints are stressed by antiplane shear forces equal to in which 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 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 are the displacement components of region of rectangles () in the direction of Oz, so the singular integral equation of this boundary value problem from  is written as follows:where and 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 and the first kind , respectively. We can substitute the below parameters to express the system of equations (39) in a simpler form: Hence, the system of equations (39) is reduced to
By assumption of the below expression for the right hand side of (41)
is found. The solution of system (26) is thereby obtained as follows:
The above formula can be converted to a dimensionless form as follows:
5. 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 there are two central cracks located symmetrically at that have equal lengths (Figure 3):The boundaries of the cracks have no traction; furthermore, the vertical edges of plates () at and are clamped, and the upper and lower horizontal edges are loaded by antiplane distributed shear loading , so thatin which 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 outside of the cracks are determined. Moreover, it is necessary to investigate the influence of adjacent crack tips and .
Suppose that the upper and lower edges of the plate are stressed by antiplane shear forces so that 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 and taking into consideration only the right hand crack with end points at and , the stress intensity factors are defined as follows [3, 4, 7, 26]:where are the shear stresses according to (31).
6. 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 and possessing a shear modulus is considered. The elastic body is rigidly clamped at 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 , it has a through-the-thickness crack in the shape of the strip ; , . Shear forces of equal intensities are acting in opposite directions along the Oz-axis on the upper (+) and lower (−) areas of the crack edgesof the crack. Additionally, at the tips the edges of the crack are joined by thin-walled inclusions with shear modulus deforming by the Winkler model (1867) that act as linear elastic springs: (Figure 4) .
The prismatic body subjected to the above-mentioned shear forces is in a state of antiplane deformation in the direction of the Oz-axis on the basic plane Oxy. The main rectangle with the central crack is a cross section of the body at the plane Oxy.
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 of rectangle then according to (4), the following equations are obtained:After some simple transformations and calculations, the following equations can be derived:withTransforming 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.
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 can be expressed in the following way:and the function from (57) is obtained:It is obvious that the function with respect to the symmetry of line in this special case and, consequently, the function 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.
The numerical analysis of the main characteristics of the stated problem can be carried out for the considered special case.
7. Numerical Results and Discussions
To solve the system of equations (40) that is summarized in system (41) with the parameters shown in (42), , factor (which represents the cracking length) equal to , , and can be used. By means of Chebyshev polynomials, the deformation parameter is obtained from (43). For example, leads to solving matrix equation (26).
The dimensionless SIF according to formula (50) using the vector , and considering the assumptions , , , , and for different crack lengths is calculated. The results are shown in Table 1, indicating that longer cracks lead to larger SIFs.
The shear stress from (49) in dependence of and taking the rigidity ratios of the two materials 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.
To solve the system of linear equations (27) for the problem with two cracks we use iteration in order that the results converge to the order . Calculation is based on function (27) , , , and , and (25), (37), and (39).
Parameters and are supposed to be equal to and in order to determine the influence of the distance between the two adjacent cracks on their fracture characteristics. The shear stress is calculated according to (49), in which function is obtained from (37), supposing the quantity . Making use of the above calculated stresses, the variation curves of can be drawn for , 0.3, 0.5, and 1.0, as shown in Figure 7.
The crack dislocation density as a function of the edges’ opening displacements may be found through formula (31). The 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