Wave field of displacement in layered structures with internal cracks

1. Introduction

Layered structures are widely used across many sectors, including the oil, energy, chemical, aerospace, and shipbuilding industries, as well as civil engineering. External loads and prolonged operation can induce internal and interfacial cracks, delamination, and other defects in these materials, reducing their durability and leading to premature failure. The acoustic emission (AE) method detects such defects during their formation and growth by recording elastic waves that propagate through the material and cause surface oscillations [1]. In AE non-destructive testing (NDT), evaluating surface displacement parameters is important because AE transducers directly record these displacements during inspections.

Solving many scientific and technical problems requires studying the laws of elastic wave propagation in media with complex physical and mechanical properties and heterogeneous structures [2-5]. It is important to consider the influence of structural inhomogeneities (inclusions, cracks, cavities, etc.) when studying the material. Under various loads (impact, cyclic), defect initiation or development causes damage accumulation that threatens the object’s integrity [6, 7].

The determination of stress intensity factors under static loading serves as a basis for assessing a material’s resistance to fracture. The characteristics of the wave field in bodies with a layered structure are analyzed using different approaches, including the Fourier integral transform method and dual integral equations [8], a boundary-element-based method [4], the finite element method [9, 10], and the extended consecutive-interpolation finite element method [11].

At the same time, dynamic loads significantly impact the structural integrity of layered composite structures. The time-periodic nature of the applied load results in complex stress and displacement wave patterns throughout the body. Numerous studies have addressed stress analysis primarily focused on two-layer structures with various defects under dynamic loading [12]. Due to significant mathematical and computational challenges, limited research has addressed the propagation of elastic waves in multi-layered heterogeneous structures under dynamic loading. We are not aware of any analytical studies of the displacement field on the surface of layered structures with internal cracks.

The traditional analytical approach for calculating the stress-strain state of a cracked solid body consists of two stages: (1) determining the crack opening displacement (COD) functions under the applied load on the crack faces; and (2) calculating the dynamic stress intensity factors (DSIF) and displacements within the body. However, for complex structures, obtaining analytical solutions at the first stage is often not feasible, and the accuracy of numerical simulations depends heavily on the reliability of the model’s input parameters.

For the investigation of acoustic waves in infinite domains with defects, the boundary integral equation method has demonstrated good suitability [1]. Previously, we proposed a method for estimating surface displacements induced by crack opening, where the COD functions are prescribed in advance [13]. This simplified modeling approach, combined with the boundary integral equation (BIE) method, enables a qualitative assessment of displacements within the solid and on its surface.

This study aims to model the propagation of elastic waves induced by displacements of opposite faces of penny-shaped Mode III cracks in a three-layer structure using the boundary integral equation method, enabling a qualitative analysis of the resulting surface displacement field.

2. Mathematical problem

Let us consider the solid body from three isotropic layers $A$, $B$, $C$ of thickness $h_A$, $h_B$, $h_C$, respectively (Fig. 1). The materials of the layers have the mass densities ${\rho }_D$, Poisson’s ratios ${\mu }_D$, and shear moduli $G_D$, where $D=A,B,C$. Each layer contains a flat penny-shaped crack of radius $a_D$ with a fixed contour and area $S_D$. The cracks are parallel to the layer surface, and their centers lie on the same line. We chose the local coordinate system $O_D{x}_1{x}_2{x}_3^D$ in the region $S_D$ so that the origin of coordinates $O_D$ coincides with the center of the crack, and the axes $x_1$, $x_2$ are in the plane of the defect.

Fig. 1General geometry of the layered composite with the penny-shaped cracks

The opposite surfaces of the cracks $S_D^{\pm }$ are subjected to shear torsional displacements, which vary in time $t$ according to a harmonic law with frequency $\omega$ and are described by the COD functions:

where $r=\sqrt{x_1^2+x_2^2}$ is the radial coordinate; $i=\sqrt{-1}$ is the imaginary unit. It is assumed that the contour of the cracks is fixed, and there is no contact between the crack surfaces during displacement. It is necessary to determine the parameters of the wave field of displacements in the body based on the given COD functions.

Fig. 2Crack location and selection of the coordinate system in layer D

Fig. 3Diagram of the location of a disk-shaped crack in a half-space B

With the specified method of arrangement of cracks in the body and the type of opening of their surfaces, only horizontally polarized transverse waves (SH waves) will propagate in the composite [14]. The problem is reduced to solving differential equations with respect to the amplitude values of the displacement vector, which for layer $D$ (Fig. 2) have the form:

where ${\mathrm{\Delta }}_3={\partial }^2/\partial x_1^2+{\partial }^2/\partial x_2^2+{\partial }^2/\partial x_3^{D2}$ is the three-dimensional Laplace operator; ${\mathit{u}}^{\left(D\right)}(u_1^{\left(D\right)},u_2^{\left(D\right)},0)$ is the displacement vector amplitude; $k_2^{\left(D\right)}=\omega /c_2^{\left(D\right)}$ is the wave number; $c_2^{\left(D\right)}=\sqrt{G_D/{\rho }_D}$ is the velocity of transverse elastic wave in the layer $D$.

Eq. (2) must be solved for three groups of boundary conditions, which we formulate only for the amplitudes of shear displacements and tangential stresses.

Group І: boundary conditions on crack surfaces:

Group ІІ: the conditions of ideal mechanical contact are met on the interface surfaces of the layers:

Group ІІІ: the lower surface of layer $A$ and the upper surface of layer $C$ (Fig. 1) are free from loads:

The total displacement field in layer $D$ is given as the sum:

where $u_{10}^{\left(D\right)}$ are the displacements from oscillations of points on the lower surface $S_{10}^D$ of the layer; $u_{20}^{\left(D\right)}$ are the displacements from oscillations of points on the upper surface $S_{20}^D$ of the layer; $u_D$ are the displacements from the opening of the crack surfaces $S_D^{\pm }$.

To solve the problem, we use the BIE method. We write the displacement in the form:

where:

are the Helmholtz potentials; unknown densities ${\beta }_{jk0}^{\left(D\right)}$ characterize the displacement of surface points $S_{k0}^{\left(D\right)}$ of the layer D; $\left|x-\xi \right|=\sqrt{(x_1-{\xi }_1{)}^2+(x_2-{\xi }_2{)}^2+(x_3^D{)}^2}$ is the distance between the fixed point $x(x_1,x_2,x_3^D)$ and the point $\xi ({\xi }_1,{\xi }_2,0)$ of integration areas $S_{k0}^D$ or $S_D$. The representation of the displacements in the form Eqs. (7), (8) automatically satisfies the Helmholtz Eq. (2) and the Sommerfeld radiation conditions at infinity [15].

3. Problem solution

The solution algorithm of the problem in general formulation is as follows. In the first step, using Eqs. (6)-(8), the displacements in layer $D$ relative to the coordinate system $O{x}_1{x}_2{x}_{310}^D$, associated with the lower surface of the layer, and relative to the coordinate system $O{x}_1{x}_2{x}_{320}^D$, associated with the upper surface of the layer, were written. Next, the boundary conditions Eq. (4) and (5) of the problem, considering the obtained representations, were written. Boundary conditions form a system of 12 two-dimensional BIEs of the Helmholtz potential type for the unknown densities ${\beta }_{jk0}^{\left(D\right)}$, $j,k=\mathrm{1,2}$, $D=A,B,C$. The right-hand sides of the equations are determined through known given COD functions. To solve the BIEs, we apply to them the two-dimensional integral Fourier transform with respect to the variables $x_1$, $x_2$, the convolution theorem [16], and some values of the Fourier transforms. As a result, we obtained a system of linear algebraic equations; the solutions of the system are the Fourier transform expressions of the desired functions. After applying the inverse Fourier transform to these solutions, we obtained integral representations of the desired functions. To determine the displacements in layer $C$, a series of algebraic transformations was performed using the values of some integrals [16, 17].

The final expression for angular displacements in layer $C$ initiated by the CODs in the layers $A$, $B$, $C$ of the layered composite is written as:

4. Results and discussions

To verify the proposed numerical method and analyze the basic patterns, a numerical solution was obtained for the case of a layer perfectly bonded to a half-space (Fig. 3). A half-space is formed here, provided that two layers ($B$ and $A$ in Fig. 1) have identical properties. An analogue of this structure may be a layered composite formed by applying a coating (for example, galvanic) to a structurally homogeneous element of significantly larger geometric dimensions than the coating’s thickness. To proceed to such a problem from the general representation Eq. (9) for displacements, we assume that $Z_A=$ 0, $h_B=d_1^B=\infty$, $E_B=$ 0, $e^{-d_1^B{R}_2^{\left(B\right)}}=$ 0, $Y_{1A}=Y_{2A}=$ 2, $Y_{1B}=Y_{2B}=$ 1, ${\mu }_B={\mu }_C$, ${\rho }_B={\rho }_C$.

The integral representation for displacements on the surface of layer $C$ takes the form:

where $\mathrm{\Gamma }\left(\tau \right)=-Z_B{E}_C\exp (-d_2^B{R}_2^{\left(B\right)})$, ${\mathrm{\Delta }}_1\left(\tau \right)=Z_C{Y}_{2C}+Z_B{Y}_{1C}$, $Z_D=G_D{R}_2^{\left(D\right)}$, $D=B,C$.

Note that for identical materials of the half-space and the layer, relation (10) reduces to displacement relations for a half-space with a crack at a depth $d=h_C+d_2^B$ [18].

For $G_C>G_B$ the function ${\mathrm{\Delta }}_1$ in the denominator of Eq. (10) has no real zeros, and the semi-infinite integral was split into three integrals with the following limits during the calculation $[0;k_2^{\left(C\right)}]$, $[k_2^{\left(C\right)};k_2^{\left(B\right)}]$ and $[k_2^{\left(B\right)};+\infty )$. The choice of radical branches $R_2^{\left(D\right)}$ is determined by the need to satisfy the Sommerfeld radiation conditions:

The integrals were solved using Gaussian elimination on the first two intervals with 40 nodes and Laguerre quadrature on the third interval with 60 nodes [19].

For $G_C<G_B$, the semi-infinite integral Eq. (10) was divided into three intervals with the following limits: $[0;k_2^{\left(B\right)}]$, $[k_2^{\left(B\right)};k_2^{\left(C\right)}]$ and $[k_2^{\left(C\right)};+\infty )$. For the first two intervals, the integral was calculated using the trapezoidal method along the Sommerfeld contour with 40 equally spaced nodes [20]. On the third interval, the integral was calculated using Laguerre quadrature with 60 nodes.

According to Eq. (10), the normalized displacement amplitudes relative to the layer thickness were determined for different wave numbers and the parameter $G=G_C/G_B$, which characterizes the stiffness ratio of the layers’ materials.

Fig. 4 shows the dependence of normalized displacements on the surface of layer $C$ (Fig. 3) on the wave number at the observation point of $r=0.5h_C$ for different $G$ and various crack depths at the half-space. The crack has a radius of $a=0.5h_C$, curve 1 – $d_2^B=2h_C$; curve 2 – $d_2^B=5h_C$; curve 3 – $d_2^B=10h_C$. The graphs show that, regardless of the ratio of the materials’ stiffnesses, the maximum displacement amplitudes on the layer surface decrease with increasing crack depth in the half-space. If the layer material has greater stiffness, then with its increase, the local maximum of the amplitudes of displacements on its surface increases, shifting towards lower frequencies (Fig. 4(a)). If the material of the half-space has greater stiffness, then with its increase, the local maximum of the amplitudes of displacements on its surface decreases, shifting towards lower frequencies (Fig. 4b).

Fig. 5 shows the dependence of the normalized displacements on the distance $r$ from the observation point for different $G$, with a fixed oscillation frequency of $k_2{h}_C=8$. The crack depth is equal to $d_2^B=h_C,$ the crack radius of $a=0.5h_C$. The distance from the epicenter of the fracture $r$ varied in the range 0,…, 10 with a step of 0.2.

Fig. 4Dependence of normalized displacements on the layer C surface at the point above the crack contour on the wave number for various crack depths d2B in the half-space for different ratios of the stiffness of the layer and half-space materials: а) G= 2; b) G= 1/2

a)

b)

Fig. 5Dependence of normalized displacements on the distance r to the observation point on the layer C surface for different ratios of stiffness of composite components

a) Curve 1 – $G=$2; curve 2 – $G=$10; curve 3 – $G=$100

b) Curve 1 – $G=$1/2; curve 2 – $G=$1/10; curve 3 – $G=$1/100

Regardless of the ratio of the body components’ stiffnesses, the amplitudes of displacements on its surface at the selected frequency are greatest if the distance from the epicenter of the fracture is equal to the layer thickness; with further distance from the epicenter, they gradually decrease and disappear. The increase in the stiffness of the layer and half-space materials reduces the displacement amplitudes. For example, with a 5-fold increase in the stiffness of the half-space, the maximum amplitude of the displacements decreases by 25 % (curves 1 and 2 in Fig. 5(a)), and with a 5-fold increase in the stiffness of the layer, the maximum amplitude of the displacements decreases by 38 % (curves 1 and 2 in Fig. 5(b)).

5. Conclusions

The solution to the 3D dynamic problem of determining the displacement field on the surface of an elastic three-layered composite for steady oscillations of the surfaces of internal penny-shaped Mode III cracks is obtained.

Numerically, the dependence of displacement on the two-layer composite’s surface on the elastic characteristics of the components, the frequency of crack surface displacement, the distance between the observation point and the fracture epicenter, and the crack depth has been established.