Multilayer medium behavior simulation in the presence of delamination defects

. The article is devoted to the study of the steady harmonic oscillation excitation processes by a surface source in a multilayer linearly elastic layered half-plane with an internal defect. The application of the work results is associated with the design and assessment of the multilayer structures state, such as geological structures, non-rigid pavements, layered composites. The study is based on the use of the boundary element method and the construction of boundary equations with respect to displacements at the defect boundary. Using the integral Fourier transform, fundamental solutions for a multilayer half-plane satisfying the stress-free conditions at the upper boundary of the domain are constructed. A feature of the construction is the good computational ability of the proposed algorithm with a large number of layers. The basic one is the construction of a solution for a homogeneous half-plane. This allows you to implement parallel computing. In the numerical implementation, a spline approximation of the defect boundary is used. Linear boundary elements are applied. As an example, a comparative analysis of the amplitude-frequency characteristics on the surface of the half-plane is given for cases of a limited defect, its absence and contact of layers without friction on the entire boundary. This allows us to draw conclusions about the possibility of identifying internal damage at the boundaries of layers by monitoring the wave field on the surface of a multilayer structure.


Introduction
The widespread use of multilayer structures under dynamic loads causes the interest of researchers in the study of wave excitation and propagation processes, the assessment of the multilayer media spectral properties, the interaction of surface and buried objects with layered structures [1][2][3].
It should be noted that numerical methods based on finite differences and finite elements in the problems of wave propagation in layered semi-bounded media have a number of limitations associated with the impossibility of exactly satisfying the radiation conditions at the outer boundary of the region.Numerical difficulties arise for multilayer media also when using infinite elements (IFE), perfectly matched layers (PLM) or absorbing boundary conditions (ABC).To solve this problem, a significant number of studies are devoted to a combination of finite (FEM) and boundary (BEM) element methods [4][5][6].
The FEM-BEM approach has also shown itself well when considered for the interaction problems of layered media with surface or partially buried objects [7,8], as well as in the presence of local cavities, inclusions or crack-like defects [9,10].
Of considerable interest for evaluating the dynamic properties of layered media is also the modeling of the conditions for the connection of layers at their interfaces.In [11][12][13][14] different variations of such conditions are considered.
In this paper, we consider the application of the boundary element method for modeling steady-state oscillations of a multilayer medium in the presence of a defect at one of the interfaces between the layers.The results of the numerical analysis of the medium amplitudefrequency characteristics on the top surface of the region are presented.

Materials and methods
The problem of steady-state harmonic oscillation excitation and propagation in a region  in the form of a multilayer half-plane is considered (see Fig. 1).The physical properties of the layered medium components   ,  = 1,2, . . .,  for an isotropic linear elastic material are determined by the elasticity modulus   , Poisson's ratio   and density   .
Wave fields for the layered medium points displacements are described by integral representations (1) obtained on the basis of the dynamic reciprocity theorem: Here, the components of the displacement vector () at the defect boundary  are unknown; ( 0 , ), ( 0 , ) -fundamental solution matrices for a layered foundation, implemented using the superposition principle according to the scheme below.In the local coordinate system for the  ℎ layer: (, ):  ∈ (−ℎ  , 0),  ∈ (−∞, ∞) the displacement amplitude functions under the action of a source concentrated in  0 have the form: Functions   () ( 0 , ),  = 1,2 satisfying the equations of motion, according to the proposed method, will be sought in the form   () =   (,1) +   (,2) +   (,3) .
Here, the terms   (,n) (, ),  = 1,2 of this representation are the solutions of the Lame equations for the half-plane with the parameters of the  ℎ medium, satisfying the stress boundary conditions: is the Lame coefficient.We represent the displacement vector  (,1) (, ) as a Fourier integral in terms of the stress vector  (,1) () transformants, denoted as  ̃(,1) (): Outline Γ is defined by the limiting absorption principle: it traverses the positive poles of the integrand from below, the negative ones -from above, and coincides with the real axis on the rest of the part.The matrix  (,1) elements look like: ,   ,   are the wave propagation velocities in the corresponding medium, expressed in terms of the elastic constants of the layers.
Determining the stress state of the layer as the sum of the corresponding solutions for two half-planes, we obtain: For the second group of terms we find: (,2) (, ) = (−1)   +1   (,1) (ℎ  − , ).
The introduced Fourier transforms for the stress functions  ̃(,) () of representations ( 2), ( 3) are unknown and must be determined from the matching conditions of the layered half-plane heterogeneous components to each other.Satisfying the equalities of the displacement and stress vector components when passing through the interfaces in the Fourier transforms, we obtain a linear algebraic equation system with 4 + 2 unknowns: where  ̃() is the common vector of unknown stresses for a multilayer structure.The found fundamental solutions have the property of the stress absence on the top surface  =   .This excludes unknown displacements on the given part of the region  boundary from representations (1).
Further directing the point  0 to the boundary of the area occupied by the foundation, we obtain a boundary integral equation system with respect to the vector :

Results and discussions
To solve the equation system ( 5) by the boundary element method on the curve  with the length , parameterization of the arc length  is introduced: The curve  is divided into  boundary elements by a network of nodal points  ℎ = {  = ℎ, ℎ =   ,  = 0,1, . . ., }.The following notations are introduced: The boundary is approximated by cubic splines.On each element   , a local coordinate system is introduced by the parameter : At   the  ℎ component of the unknown vector of unknowns can be represented as: (), where the basis function system   () determines the approximation nature of the unknowns on the element.In this paper, power functions   () =  −1 were used, as well as polynomials:  = 2,  1 = 0.5(1 − ),  2 = 0.5(1 + ) -linear boundary elements.
The midpoints of boundary elements and docking nodes were chosen as collocation points.
The resulting linear algebraic equation system was numerically inverted by the Gaussian factorization method.
To carry out a comparative analysis, we considered problems with rigid contact on the entire interface (the absence of a defect -task 1), as well as with contact of layers without friction on the entire interface (task 2).In these cases, in the displacement representation (1), the integral over the border of the defect  disappears.Problems (1) and (2) differ from each other only in the form of the function ( 0 , ), where the term   (,3) is absent.
Consider the nature of a three-layer medium amplitude-frequency characteristics, determined by a set of parameters (see

Conclusions
Summarizing the results of the study, we can note the following: 1.When a slip zone appears at the layer interface, an increase in the first maximum of the observation points vertical displacement amplitude-frequency characteristic on the surface of a multilayer medium is observed, in our case, in the vicinity of 10 Hz frequency.
2. An increase in the amplitude of oscillations in a higher frequency region is noted.In this case, the presence of maximum values in this region is determined with greater certainty by the behavior of the selected observation point phase-frequency characteristic.
3. The different sensitivities of the displacement vector components to changes in the layer adhesion conditions leads to the need to consider the amplitude-frequency characteristics   ,   not of each of them separately, but in the form of the horizontal to vertical component moduli ratio |  | |  | ⁄ .

Fig. 1 .
Fig. 1.Study area D Oscillations are generated by an intensity source () distributed in a limited area Ω on the top surface of the medium  = −  with the frequency .The multiplier (−) is omitted below.There is a defect  in the form of delamination at the boundary between the layers.The conditions of rigid contact are set on the remaining media conjugation boundaries.The physical properties of the layered medium components   ,  = 1,2, . . .,  for an isotropic linear elastic material are determined by the elasticity modulus   , Poisson's ratio   and density   .Wave fields for the layered medium points displacements are described by integral representations (1) obtained on the basis of the dynamic reciprocity theorem:

Figure 2
Figure2shows the amplitude-frequency characteristics of the vertical displacement of the medium point with coordinates  = −6,  = 10 (on the top surface) with rigid adhesion of all structure (1) layers, for the case of the upper layer contact without friction (2) and for the case of a defect at the interface of layers 2 and 3 with length of 4 m (3).An increase in the quality factor of the resonance is observed in the vicinity of 10 Hz frequency.

Fig. 2 .
Fig. 2. Amplitude-frequency characteristics of the vertical medium point displacement

Table 1 .
Parameters of layers