# Dynamic Response of Crossanisotropic Layered Media with Imperfect Interface under Axisymmetric Loads

Authors: Zejun Han; Chengchao Guo; Linqing Yang; Jin Zhang; Kai Sun
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

Considering the incomplete contact between the adjacent layers of crossanisotropic multilayered half-space media, an algorithm was derived based on the spectral element method (SEM) to determine the dynamic response at any point of the medium when an axisymmetric harmonic load is applied. The algorithm used two equivalent contact parameters to simulate the relationship between the discontinuity of the displacement and stress. The SEM was employed for the summation in the integral calculation to avoid the infinite integral in the integral transformation and greatly improve the efficiency of the solution to the dynamic response of the crossanisotropic layered half-space medium. The accuracy of the developed method is verified by the numerical examples as well as the efficiency of the developed method. Finally, an extensive parametric analysis of the contact parameters of the interface between the layers of the layered medium, namely, the point of the application of the external load, the excitation frequency, and the heterogeneous properties of the layered medium was conducted, providing a reliable numerical basis for engineering practice.


---

## Body

## 1. Introduction

Layered structures such as pavement, soil, and composite materials are common heterogeneous media, and the study of the dynamic response of layered mediums is widely applied to transportation, civil engineering, and material science. Considering the complexity of engineering calculations, in the early related literature, layered mediums were often assumed to be isotropic. However, many-layered media are anisotropic due to their construction technology or natural structure. With the development of computer technology, some scholars have proposed that crossanisotropic properties can be seen as anisotropic properties in a special case and that studying the dynamic response of layered mediums with crossanisotropic properties can reveal the dynamic response law of layered media more accurately, which conforms to the actual engineering situation [1–3].There are many literature works about the dynamic displacement and stress responses of isotropic media [4–9]. However, the assumption of isotropy cannot reflect the anisotropic properties of layered media in the horizontal and vertical directions. In an early work, Robert [10] proposed that the propagation of waves in a crossanisotropic medium is quite different from that in an anisotropic medium. Therefore, more scholars became interested in studying the dynamic response of crossanisotropic layered mediums.For instance, Buchwald [11] and Payton [1] examined the elastodynamics of crossanisotropic half-space bodies under surface loads. Rajapakse and Wang [12] made slight modifications to the method proposed by Buchwald [11] by using the three potential functions proposed by the latter to simplify the equation of motion into two coupled partial differential equations and a separate partial differential equation. They obtained Green’s function of the crossanisotropic medium under internal harmonic loads but failed to solve the problem completely. Later, Rahimian et al. [13] proposed a new, effective, elastodynamic, potential method and found the basic analytical solution to the displacement of a crossanisotropic half-space medium under surface loads by introducing two potential functions for the displacement. Khojasteh et al. [14] provided an analytical solution to the dynamic response of a crossanisotropic half-space body under internal harmonic loads by using this method.However, mediums have different degrees of heterogeneity along the vertical direction in many practical situations. Compared with homogeneous and half-space media, it is more difficult to determine the dynamic response of crossanisotropic layered mediums. Research on this issue has gradually increased in recent years, and methods such as numerical solutions, analytical solutions, and semianalytical and seminumerical solutions have emerged. Numerical solutions, mainly including the finite element method [15, 16] and the boundary element method [17, 18], usually require dividing a large number of meshes using incremental algorithms and calculating the mechanical responses of all the mesh elements at the same time, which demands powerful computation systems and huge computational cost; meanwhile, they have low efficiency. Rokhlin and Wang [19] presented an efficient recursive algorithm for wave propagation in multilayered, generally anisotropic media with the stiffness matrix method. Ba et al. [20] also solved the dynamic response of three-dimensional crossanisotropic saturated layered media under the action of concentrated force and pore pressure using the accurate stiffness matrix method. Ai et al. [21] and Ai and Li [22] proposed the analytical layer-element method to find the frequency domain solution for the crossanisotropic layered mediums and assembled the global stiffness matrix of the half-space body in the cylindrical coordinate system. Later, Ai et al. [23] further improved the method and obtained the dynamic response of an arbitrary 3-D crossanisotropic layered medium under moving loads. Ai et al. [24] analyzed the dynamic response of crossanisotropic saturated media under circular moving loads, in which the vertical and tangential loads can be simultaneously considered. In another work, Ba et al. [20] employed the analytical layer-element method to study the dynamic response of a three-dimensional crossanisotropic half-space body moving at a constant speed in the horizontal direction under a concentrated load. Ba et al. [25] presented an algorithm to obtain the dynamic responses of a crossanisotropic layer saturated under dynamic loads and pore pressure. In addition, Khojasteh et al. [26] used the dynamic potential method to determine the dynamic response of a double-layered crossanisotropic half-space body under the action of a concentrated load and then extended the algorithm to obtain 3-D Green’s function of a crossanisotropic multilayered half-space medium [27]. The thin-layer method, an effective semianalytical technique for calculating Green’s function of layered mediums, was used by Barbosa and Kausel [28] to establish the dynamic response of crossanisotropic layered media. Lin et al. [29] and Han et al. [30] obtained the dynamic response of a crossanisotropic layered foundation under arbitrary loads by using a precise integration algorithm and carried out a parameter analysis of the crossanisotropic properties of the foundation. Li et al. [31] obtained the dynamic response of a stratified crossanisotropic half-space with a poroelastic interlayer due to a buried moving source. Han et al. [32] presented a modified scaled boundary finite element method (SBFEM) to obtain the dynamic impedance matrix of the anisotropic layered media based on the continued-fraction method. The research mentioned previously on crossanisotropic layered mediums is a supplement to the research on crossanisotropic homogeneous bodies and truly reflects the mechanical behavior of layered media as the difference in the material properties of the layers has been considered. However, few works have taken the contact between the layers of layered media into account.At present, the contact between the layers of a layered medium is often assumed to be complete in engineering designs. However, in actual projects, it is almost incomplete due to the difference in the materials, the construction conditions, and the construction technology, so this assumption does not reflect the real conditions. The existing studies have shown that the contact between the layers seriously affects the dynamic response of layered structures [33]. Therefore, it is necessary to consider the influence of the incomplete contact between the layers of a layered medium on its dynamic response, so the examination of the response of crossanisotropic layered mediums can be more practical. The existing literature presents a variety of models such as the thin-layer model [34], the spring model (Goodman model) [35], and the density model [36] to deal with the problem of contact between the layers of layered mediums. Liu and Pan [37] proposed a transfer-matrix model (the dual variable and position method) based on a variety of traditional interlayer contact models, which is more suitable for thin-layered structures at a high excitation frequency.This paper takes account of the incomplete contact between the layers of a crossanisotropic layered medium and establishes a new SEM to calculate the displacement response and the vertical stress response of the layered medium under a vertical load. Compared with the previous numerical integration technique, the proposed algorithm has higher efficiency. Then, parameter analyses related to the point of the application of the external load, the excitation frequency, and the heterogeneous properties of the layered media have been carried out. Relevant results can provide a reliable theoretical and numerical basis for engineering practice.


## 2. Problem Description

We assumed thatn layers of the layered medium were located in the homogeneous half-space medium (see Figure 1), and each layer of the medium was presumed to be crossanisotropic. Thus, a cylindrical coordinate system (r, θ, and z) was established as shown in Figure 1, and the origin of the coordinate (O) was set on the surface of the layered medium with the z-axis vertically downward and the r-axis in the radial direction. The contact between the other adjacent layers was complete, except for the imperfect interface. The interface zj between layers j and j + 1 is imperfect with values on its upper and lower sides being indicated by −0 and +0. The displacements on the imperfect interface satisfied the discontinuous boundary conditions. An axisymmetric load was applied within the range of Δr. The expression of the load was defined as Pt=Sreiωt, where ω is the angular frequency and Sr can be expressed by the following equation:(1)Sr=q,for0⩽r⩽Δr,0,forr>Δr,Figure 1
The computational model of the layered medium.


## 3. Wave Equations for Axisymmetric Problem

### 3.1. Governing Equations

The governing equations of axisymmetric wave motion in the absence of body forces for a continuous medium in a cylindrical coordinate system are given by the following equation:(2)∂σr∂r+∂τzr∂z+σr−σθr=ρ∂2u∂t2,∂τzr∂r+∂σz∂z+τzrr=ρ∂2w∂t2,where σr, σθ, and σz are the normal stress in the r, θ, and z directions, respectively; τzr is the shear stress in the plane z−r; ρ is the density of the medium; u and w are the displacement components in r and z directions, respectively.The displacementsu and w in the governing equations (2) can be expressed in the Helmholtz decomposition as follows:(3)u=∂ξ1∂r−∂ξ2∂z,w=∂ξ1∂z+1r∂rξ2∂r,where ξ1 and ξ2 are the gradients of the scalar potential and curl of the vector potential, respectively.The constitutive equations for a crossanisotropic elastic medium is defined as follows:(4)σr=c11∂u∂r+c12ur+c13∂w∂z,σθ=c12∂u∂r+c11ur+c13∂w∂z,σz=c13∂u∂r+ur+c33∂w∂z,τzr=c44∂w∂r+∂u∂z,where cij is the elastic constant of the crossanisotropic medium and c11, c12, c13, c33, c44, and χ are defined as follows:(5)c11=Ehh1−χvvh21+vhh1−vhh−2χvvh2,c12=Ehhvhh+χvvh21+vhh1−vhh−2χvvh2,c13=Ehhvvh1−vhh−2χvvh2,c33=Evh1−vhh1−vhh−2χvvh2,c44=Gvh,χ=EhhEvh,where Ehh and Evh are Young’s modulus of the crossanisotropic medium in the horizontal and vertical directions, respectively; vhh and vvh represent Poisson’s ratios on the horizontal and vertical planes, respectively; Gvh is the shear modulus on the vertical plane; χ stands for the ratio of horizontal Young’s modulus of the crossanisotropic medium to the vertical one.The displacement wave equation can be obtained by substituting equation (4) into equation (2):(6)c11∇2−1r2u+c44∂2u∂z2+c13+c44∂2w∂r∂z=ρ∂2u∂t2,(7)c13+c44∂2∂r∂z+1r∂∂zu+c33∂2w∂z2+c44∇2w=ρ∂2w∂t2,where ∇2=1/r∂/∂rr∂/∂r.The displacement expressed by the scalar potential in equation (3) is substituted into equation (6), and the scalar potential needs to satisfy the following equation after algebraic simplification.(8)∂2ξ1∂r2+∂ξ1r∂r+∂2ξ1κ1∂z2=1cp12∂2ξ1∂t2,(9)∂2ξ2κ2∂z2+∂2ξ2∂r2+1r∂ξ2∂r−ξ2r2=1cz12∂2ξ2∂t2,with(10)κ1=c11c13+2c44,cp12=c11ρ,κ2=c11−c13−c44c44,cz12=κ2c44ρ.Similarly, by substituting equation (3) into equation (7), the following equation can be obtained after algebraic simplification.(11)∂2ξ1∂r2+∂ξ1r∂r+∂2ξ1κ3∂z2=1cp22∂2ξ1∂t2,(12)∂2ξ2κ4∂z2+∂2ξ2∂r2+1r∂ξ2∂r−ξ2r2=1cz22∂2ξ2∂t2,with(13)κ3=c13+2c44c33,cp22=κ3c33ρ,κ4=c44c33−c13−c44,cz22=c44ρ.It is necessary to solve potential functionsξ1 and ξ2 of equations (8) and (9) to avoid discussing the propagation paths of the waves in the r and z directions. Obviously, the solutions of equations (8) and (9) are similar to those of equations (11) and (12). Thus, for the sake of simplicity, only the solutions of equations (8) and (9) are presented. The detailed derivation process can be found in Appendix A.


### 3.2. Interlayer Contact

It is usually assumed that the contact between the layers of layered mediums is complete, that is, the displacement and stress are continuous at the interface. However, there is usually a certain relative slip between the adjacent layers of a layered medium due to the process conditions and natural structure of the layered mediums, that is, the contact between the layers is incomplete, and the interface is imperfect. In this work, interfacezj of the layered medium displayed in Figure 1 is assumed to be imperfect, so the relationship between the stress on layer z = zj (τ^zr, σ^z) and the displacement of layer z = zj (u^, w^) and the contact parameters (tr and tz) of the interface is defined as follows:(14)σ^zzj+0=σ^zzj−0,σ^zzj+0=tzw^zj+0−w^zj−0,τ^zrzj+0=τ^zrzj−0,τ^zrzj+0=tru^zj+0−u^zj−0,where tr and tz are the contact parameters of the interface and can be regarded as the stiffness modulus of the horizontal and vertical interfaces of the interface z = zj, respectively.Equation (14) can be expressed in a matrix form as follows:(15)u^zj−0w^zj−0τ^zrzj+0σ^zzj+0=10−1tr0010−1tz00100001u^zj+0w^zj+0τ^zrzj−0σ^zzj−0.Moreover, we assume that the modulus of the horizontal interface is the same in any direction of the horizontal plane, so it can be regarded as axisymmetric. We can deduce from the formula that the displacement of both sides of interfacez = zj is almost continuous when tr and tz approach infinity. On the contrary, the contact between both sides of interface z = zj is completely slidable when tr and tz approach zero. In addition, both sides of interface z = zj slide only in the horizontal direction when tr has a definite value and tz approaches infinity, while the vertical displacement is approximately continuous. The contact between the upper and lower surfaces in the horizontal direction is totally incomplete when tr is equal to zero.


## 3.1. Governing Equations

The governing equations of axisymmetric wave motion in the absence of body forces for a continuous medium in a cylindrical coordinate system are given by the following equation:(2)∂σr∂r+∂τzr∂z+σr−σθr=ρ∂2u∂t2,∂τzr∂r+∂σz∂z+τzrr=ρ∂2w∂t2,where σr, σθ, and σz are the normal stress in the r, θ, and z directions, respectively; τzr is the shear stress in the plane z−r; ρ is the density of the medium; u and w are the displacement components in r and z directions, respectively.The displacementsu and w in the governing equations (2) can be expressed in the Helmholtz decomposition as follows:(3)u=∂ξ1∂r−∂ξ2∂z,w=∂ξ1∂z+1r∂rξ2∂r,where ξ1 and ξ2 are the gradients of the scalar potential and curl of the vector potential, respectively.The constitutive equations for a crossanisotropic elastic medium is defined as follows:(4)σr=c11∂u∂r+c12ur+c13∂w∂z,σθ=c12∂u∂r+c11ur+c13∂w∂z,σz=c13∂u∂r+ur+c33∂w∂z,τzr=c44∂w∂r+∂u∂z,where cij is the elastic constant of the crossanisotropic medium and c11, c12, c13, c33, c44, and χ are defined as follows:(5)c11=Ehh1−χvvh21+vhh1−vhh−2χvvh2,c12=Ehhvhh+χvvh21+vhh1−vhh−2χvvh2,c13=Ehhvvh1−vhh−2χvvh2,c33=Evh1−vhh1−vhh−2χvvh2,c44=Gvh,χ=EhhEvh,where Ehh and Evh are Young’s modulus of the crossanisotropic medium in the horizontal and vertical directions, respectively; vhh and vvh represent Poisson’s ratios on the horizontal and vertical planes, respectively; Gvh is the shear modulus on the vertical plane; χ stands for the ratio of horizontal Young’s modulus of the crossanisotropic medium to the vertical one.The displacement wave equation can be obtained by substituting equation (4) into equation (2):(6)c11∇2−1r2u+c44∂2u∂z2+c13+c44∂2w∂r∂z=ρ∂2u∂t2,(7)c13+c44∂2∂r∂z+1r∂∂zu+c33∂2w∂z2+c44∇2w=ρ∂2w∂t2,where ∇2=1/r∂/∂rr∂/∂r.The displacement expressed by the scalar potential in equation (3) is substituted into equation (6), and the scalar potential needs to satisfy the following equation after algebraic simplification.(8)∂2ξ1∂r2+∂ξ1r∂r+∂2ξ1κ1∂z2=1cp12∂2ξ1∂t2,(9)∂2ξ2κ2∂z2+∂2ξ2∂r2+1r∂ξ2∂r−ξ2r2=1cz12∂2ξ2∂t2,with(10)κ1=c11c13+2c44,cp12=c11ρ,κ2=c11−c13−c44c44,cz12=κ2c44ρ.Similarly, by substituting equation (3) into equation (7), the following equation can be obtained after algebraic simplification.(11)∂2ξ1∂r2+∂ξ1r∂r+∂2ξ1κ3∂z2=1cp22∂2ξ1∂t2,(12)∂2ξ2κ4∂z2+∂2ξ2∂r2+1r∂ξ2∂r−ξ2r2=1cz22∂2ξ2∂t2,with(13)κ3=c13+2c44c33,cp22=κ3c33ρ,κ4=c44c33−c13−c44,cz22=c44ρ.It is necessary to solve potential functionsξ1 and ξ2 of equations (8) and (9) to avoid discussing the propagation paths of the waves in the r and z directions. Obviously, the solutions of equations (8) and (9) are similar to those of equations (11) and (12). Thus, for the sake of simplicity, only the solutions of equations (8) and (9) are presented. The detailed derivation process can be found in Appendix A.


## 3.2. Interlayer Contact

It is usually assumed that the contact between the layers of layered mediums is complete, that is, the displacement and stress are continuous at the interface. However, there is usually a certain relative slip between the adjacent layers of a layered medium due to the process conditions and natural structure of the layered mediums, that is, the contact between the layers is incomplete, and the interface is imperfect. In this work, interfacezj of the layered medium displayed in Figure 1 is assumed to be imperfect, so the relationship between the stress on layer z = zj (τ^zr, σ^z) and the displacement of layer z = zj (u^, w^) and the contact parameters (tr and tz) of the interface is defined as follows:(14)σ^zzj+0=σ^zzj−0,σ^zzj+0=tzw^zj+0−w^zj−0,τ^zrzj+0=τ^zrzj−0,τ^zrzj+0=tru^zj+0−u^zj−0,where tr and tz are the contact parameters of the interface and can be regarded as the stiffness modulus of the horizontal and vertical interfaces of the interface z = zj, respectively.Equation (14) can be expressed in a matrix form as follows:(15)u^zj−0w^zj−0τ^zrzj+0σ^zzj+0=10−1tr0010−1tz00100001u^zj+0w^zj+0τ^zrzj−0σ^zzj−0.Moreover, we assume that the modulus of the horizontal interface is the same in any direction of the horizontal plane, so it can be regarded as axisymmetric. We can deduce from the formula that the displacement of both sides of interfacez = zj is almost continuous when tr and tz approach infinity. On the contrary, the contact between both sides of interface z = zj is completely slidable when tr and tz approach zero. In addition, both sides of interface z = zj slide only in the horizontal direction when tr has a definite value and tz approaches infinity, while the vertical displacement is approximately continuous. The contact between the upper and lower surfaces in the horizontal direction is totally incomplete when tr is equal to zero.


## 4. Summation of the Spectral Element Method

The superposition principle is the key theory of the SEM. Both equations (8) and (9) are linear homogeneous equations. Following the superposition principle, the solutions of equations (8) and (9) can be executed by summing an infinite number of wave numbers ηm (m = 1, 2, …, M) in Section 3.2:(16)ξ1r,z,ω=∑m=1MAme−iηpzmzJ0ηmre−iωt,(17)ξ2r,z,ω=∑m=1MBme−iηszmzJ1ηmre−iωt.The sum of theM wave numbers can be calculated by the Fourier–Bessel series. Equations (16) and (17) replace the infinite integral, which facilitates the summation of the wave numbers and improves the efficiency of the solution algorithm. For infinite integration, the integrand has singular points, especially if the damping ratio is very small, further complicating the calculations. However, the proposed method in this paper can effectively avoid this problem.The solution to the displacement of the axisymmetric problem can be obtained by substituting the general solutions forξ1 and ξ2 into the equation.(18)ur,z,twr,z,t=∑m=1MG^rzηmJ1ηmrG^zzηmJ0ηmrθ^meiωt,where G^ηm,z is the flexibility coefficient of the medium in the translated domain; θ^m can be calculated by the following equations:(19)θ^m=2R2J12αm∫0arJ0αmRr,dr=2ΔrαmRJ12αmJ1αmRΔr.ExpressingS (r) by means of Fourier–Bessel series, we obtain(20)Sr=∑m=1Mθ^mJ0ηmr.


## 5. Spectral Axisymmetric Element Formulation

The solution forG^ηm,z in Equation (18) in the translated domain can be found by solving a two-point boundary value problem. There are two situations that must be considered for the layered half-space medium: one is a homogeneous single-layered medium, and the other is a homogeneous half-space medium.


### 5.1. Homogeneous Single-Layered Medium

As shown in Figure2, the wave propagates on two parallel upper and lower boundary surfaces in a homogeneous single-layered medium. Both points have radial and vertical degrees of freedom. In the horizontal direction, it is deemed that the wave disappears when the propagation distance of the wave is equal to R.Figure 2
The homogeneous single-layered element.The wave response of any point in a single-layered medium is determined by the superposition of the incident wave and the reflected one. The external force produces the incident waves, which propagate and touch the vertical boundary to produce the reflected waves. The wave propagation can be decomposed into longitudinal and shear waves. Assuming a single-layered medium with the same material properties has parallel upper and lower boundaries. The displacement vectorU^ and stress vector T^ in the single-layered medium are defined as follows:(21)U^=u^w^T,T^=−τ^zr−σ^zT.Then, the state equation betweenU^ and T^ in the frequency-wave number domain can be expressed as follows:(22)V′^=HV^withV^=U^T^,H=−K22−1K21−K22−1−K11+K12K22−1K21+ρωn2IK12K22−1,where H is the characteristic matrix of the layered medium; #′ stands for the first-order partial derivative of # with respect to the coordinate z; I is the identity matrix; the matrices K11, K12, K21, and K22 are defined as follows:(23)K11=κm2d1100d44,K21=−K12T=iκm0d44d130,K22=d4400d33.The general solution to equation (22) in every single-layered medium can be obtained by the precise integration method (PIM) accurately and easily [29]. The dual vector form of the wave motion equation makes its solution easy to conduct and maintains numerical stability. Based on the solution of the differential equation (Equation (22)), for the conservative linear system, the following matrix form relationship between the displacements and stresses at the two ends of the layer can be obtained by the PIM:(24)U^zj−1T^zj=E^11mjE^12mjE^21mjE^22mjU^zjT^zj−1.


### 5.2. Considering Incomplete Contact

The relationship between the surface force and the displacement of the nodes on the upper and lower sides of layerj (j = 1, 2, …, n) in a single-layered medium can be expressed according to the following equation:(25)U^zj−1T^zj−0=E^11mjE^12mjE^21mjE^22mjU^zj−0T^zj−1.The relationship of the surface force with the displacements of the upper surface (z = zj − 0) and the lower surface (z = zj + 0) of interface z = zj can be given according to the following equation:(26)U^zj−0T^zj+0=It0IU^zj+0T^zj−0,where I represents a 2 × 2 identity matrix and t is defined as follows:(27)t=−1tr00−1tz.Considering the incomplete contact,T (zj − 0) and U (zj − 0) are eliminated based on(28)U^zj−1T^zj+0=E^11mj∗E^12mj∗E^21mj∗E^22mj∗U^zj+0T^zj−1,where(29)E^11mj∗=E^11mj+E^11mjtI−E^21mjt−1E^21mj,E^12mj∗=E^12mj+E^11mjtI−E^21mjt−1E^22mj,E^21mj∗=I−E^21mjt−1E^21mj,E^22mj∗=I−E^21mjt−1E^22mj.It is noted that ift=0, then equation (29) is reduced to the submatrices for the perfect interface case. If there is imperfect contact between the adjacent layers, we just need to use the matrices calculated in equation (29) to propagate the solution.


### 5.3. Homogeneous Half-Space Medium

As depicted in Figure3, the homogeneous half-space medium is similar to a carrier that transfers system energy, and the waves propagate in one direction without reflection. Two cases of the layered half-space are considered: the layered stratum on a rigid base and the layered stratum on a half-space. In the former case, the boundary condition is as follows:(30)U^zn=0.Figure 3
The homogeneous half-space medium.In case the multilayered soil overlies an elastic half-space, the radiative condition should be considered. To determine the boundary condition in this case, some treatment is conducted as described in the following paragraph.The eigenvalue problem of the homogeneous crossanisotropic half-space is solved as follows:(31)HΨ=ΨΛ,with the eigenvalues Λ and eigenvectors Ψ, and,(32)Λ=λi−λi,Ψ=Ψ11Ψ12Ψ21Ψ22,where the real parts of all elements λi are positive.Introducing a vector,(33)b=Ψ−1V.^Then, equation (22) becomes(34)b′=Λbwithbz=expλizexp−λizc1c2,where c1 and c2 are integration constants.Substituting equation (34) into equation (33) yields(35)U^T^=Ψ11Ψ12Ψ21Ψ22expλizc1exp−λizc2.For a homogeneous half-space, the displacementsU^ and stresses T^ must remain finite, which leads to c1=0. Applying Equation (35) results in the boundary condition at the surface (z = 0) of the half-space.(36)T^zn=S^∞U^znwithS^∞=Ψ22Ψ12−1.


### 5.4. Assembling Structural Stiffness Matrix

The stiffness matrix of the entire system (S^globalηm) can be assembled using the stiffness matrix of each layer and the half-space stiffness matrix with the same wave number after obtaining the stiffness matrices of the crossanisotropic single-layered medium and the homogeneous half-space medium. The assembly follows the classical assembly theory of the finite element stiffness matrix, and the equation for the case of the n-layered medium and the homogeneous half-space medium at the bottom can be defined as follows:(37)P^z0P^z1P^z2⋮P^zj⋮P^zn=S^11m1S^12m10⋯000S^21m1S^22m1+S^11m2S^12m2⋯0000S^21m2S^22m2⋯000⋮⋮⋮⋱⋮⋮⋮000⋯S^22mj−1+S^11mj⋯0000⋯⋮⋱S^12mn000⋯0S^21mnS^22mn+S^∞U^z0U^z1U^z2⋮U^zj⋮U^zn,or simplified to(38)P^=S^globalηmU,^where P^zii=0,1,2......,n is the external load vector of the node at the depth of zi, and,(39)P^zi=Tz0,i=0,−Tzi−1+Tzi,i≥1.S^mj is the stiffness matrix of the jth layer defined by Equation (40), and its value can be obtained from the following equation:(40)S^mj=E^12mj−1−E^12mj−1E^11mj−E^22mjE^12mj−1−E^21mj+E^22mjE^12mj−1E^11mj.If the layered stratum rests on a rigid base, the boundary conditions at the bottom interface of the layered stratum must satisfy Equation (30). Then, the dynamic (Equation (37)) for the entire n layers should be modified by eliminating the last (i.e., n + 1) row and the last (i.e., n + 1) column.The stiffness matrix in the SEM has the same characteristics as that in the general finite element method.S^globalηm is a symmetric band matrix, and its inverse matrix is G^ηm in equation (18). Therefore, the radial and vertical displacements in the spatial domain can be obtained by substituting displacements u^ and w^ in the wave number domain, which are determined by substituting the displacement calculated from equation (37) into equation (18) when a vertical axially symmetric load is applied.(41)ur,z,t=∑m=1Mu^mz,ηm,ωθ^mJ1ηmreiωt,wr,z,t=∑m=1Mw^mz,ηm,ωθ^mJ0ηmreiωt.The vertical stress can also be calculated by substituting equation (41) into equation (4) as follows:(42)σzr,z,t=c13∑m=1Mu^mz,ηm,ωθ^mηmJ0ηmreiωt+c33∑m=1M∂u^mz,ηm,ω∂zθ^mJ1ηmreiωt.


## 5.1. Homogeneous Single-Layered Medium

As shown in Figure2, the wave propagates on two parallel upper and lower boundary surfaces in a homogeneous single-layered medium. Both points have radial and vertical degrees of freedom. In the horizontal direction, it is deemed that the wave disappears when the propagation distance of the wave is equal to R.Figure 2
The homogeneous single-layered element.The wave response of any point in a single-layered medium is determined by the superposition of the incident wave and the reflected one. The external force produces the incident waves, which propagate and touch the vertical boundary to produce the reflected waves. The wave propagation can be decomposed into longitudinal and shear waves. Assuming a single-layered medium with the same material properties has parallel upper and lower boundaries. The displacement vectorU^ and stress vector T^ in the single-layered medium are defined as follows:(21)U^=u^w^T,T^=−τ^zr−σ^zT.Then, the state equation betweenU^ and T^ in the frequency-wave number domain can be expressed as follows:(22)V′^=HV^withV^=U^T^,H=−K22−1K21−K22−1−K11+K12K22−1K21+ρωn2IK12K22−1,where H is the characteristic matrix of the layered medium; #′ stands for the first-order partial derivative of # with respect to the coordinate z; I is the identity matrix; the matrices K11, K12, K21, and K22 are defined as follows:(23)K11=κm2d1100d44,K21=−K12T=iκm0d44d130,K22=d4400d33.The general solution to equation (22) in every single-layered medium can be obtained by the precise integration method (PIM) accurately and easily [29]. The dual vector form of the wave motion equation makes its solution easy to conduct and maintains numerical stability. Based on the solution of the differential equation (Equation (22)), for the conservative linear system, the following matrix form relationship between the displacements and stresses at the two ends of the layer can be obtained by the PIM:(24)U^zj−1T^zj=E^11mjE^12mjE^21mjE^22mjU^zjT^zj−1.


## 5.2. Considering Incomplete Contact

The relationship between the surface force and the displacement of the nodes on the upper and lower sides of layerj (j = 1, 2, …, n) in a single-layered medium can be expressed according to the following equation:(25)U^zj−1T^zj−0=E^11mjE^12mjE^21mjE^22mjU^zj−0T^zj−1.The relationship of the surface force with the displacements of the upper surface (z = zj − 0) and the lower surface (z = zj + 0) of interface z = zj can be given according to the following equation:(26)U^zj−0T^zj+0=It0IU^zj+0T^zj−0,where I represents a 2 × 2 identity matrix and t is defined as follows:(27)t=−1tr00−1tz.Considering the incomplete contact,T (zj − 0) and U (zj − 0) are eliminated based on(28)U^zj−1T^zj+0=E^11mj∗E^12mj∗E^21mj∗E^22mj∗U^zj+0T^zj−1,where(29)E^11mj∗=E^11mj+E^11mjtI−E^21mjt−1E^21mj,E^12mj∗=E^12mj+E^11mjtI−E^21mjt−1E^22mj,E^21mj∗=I−E^21mjt−1E^21mj,E^22mj∗=I−E^21mjt−1E^22mj.It is noted that ift=0, then equation (29) is reduced to the submatrices for the perfect interface case. If there is imperfect contact between the adjacent layers, we just need to use the matrices calculated in equation (29) to propagate the solution.


## 5.3. Homogeneous Half-Space Medium

As depicted in Figure3, the homogeneous half-space medium is similar to a carrier that transfers system energy, and the waves propagate in one direction without reflection. Two cases of the layered half-space are considered: the layered stratum on a rigid base and the layered stratum on a half-space. In the former case, the boundary condition is as follows:(30)U^zn=0.Figure 3
The homogeneous half-space medium.In case the multilayered soil overlies an elastic half-space, the radiative condition should be considered. To determine the boundary condition in this case, some treatment is conducted as described in the following paragraph.The eigenvalue problem of the homogeneous crossanisotropic half-space is solved as follows:(31)HΨ=ΨΛ,with the eigenvalues Λ and eigenvectors Ψ, and,(32)Λ=λi−λi,Ψ=Ψ11Ψ12Ψ21Ψ22,where the real parts of all elements λi are positive.Introducing a vector,(33)b=Ψ−1V.^Then, equation (22) becomes(34)b′=Λbwithbz=expλizexp−λizc1c2,where c1 and c2 are integration constants.Substituting equation (34) into equation (33) yields(35)U^T^=Ψ11Ψ12Ψ21Ψ22expλizc1exp−λizc2.For a homogeneous half-space, the displacementsU^ and stresses T^ must remain finite, which leads to c1=0. Applying Equation (35) results in the boundary condition at the surface (z = 0) of the half-space.(36)T^zn=S^∞U^znwithS^∞=Ψ22Ψ12−1.


## 5.4. Assembling Structural Stiffness Matrix

The stiffness matrix of the entire system (S^globalηm) can be assembled using the stiffness matrix of each layer and the half-space stiffness matrix with the same wave number after obtaining the stiffness matrices of the crossanisotropic single-layered medium and the homogeneous half-space medium. The assembly follows the classical assembly theory of the finite element stiffness matrix, and the equation for the case of the n-layered medium and the homogeneous half-space medium at the bottom can be defined as follows:(37)P^z0P^z1P^z2⋮P^zj⋮P^zn=S^11m1S^12m10⋯000S^21m1S^22m1+S^11m2S^12m2⋯0000S^21m2S^22m2⋯000⋮⋮⋮⋱⋮⋮⋮000⋯S^22mj−1+S^11mj⋯0000⋯⋮⋱S^12mn000⋯0S^21mnS^22mn+S^∞U^z0U^z1U^z2⋮U^zj⋮U^zn,or simplified to(38)P^=S^globalηmU,^where P^zii=0,1,2......,n is the external load vector of the node at the depth of zi, and,(39)P^zi=Tz0,i=0,−Tzi−1+Tzi,i≥1.S^mj is the stiffness matrix of the jth layer defined by Equation (40), and its value can be obtained from the following equation:(40)S^mj=E^12mj−1−E^12mj−1E^11mj−E^22mjE^12mj−1−E^21mj+E^22mjE^12mj−1E^11mj.If the layered stratum rests on a rigid base, the boundary conditions at the bottom interface of the layered stratum must satisfy Equation (30). Then, the dynamic (Equation (37)) for the entire n layers should be modified by eliminating the last (i.e., n + 1) row and the last (i.e., n + 1) column.The stiffness matrix in the SEM has the same characteristics as that in the general finite element method.S^globalηm is a symmetric band matrix, and its inverse matrix is G^ηm in equation (18). Therefore, the radial and vertical displacements in the spatial domain can be obtained by substituting displacements u^ and w^ in the wave number domain, which are determined by substituting the displacement calculated from equation (37) into equation (18) when a vertical axially symmetric load is applied.(41)ur,z,t=∑m=1Mu^mz,ηm,ωθ^mJ1ηmreiωt,wr,z,t=∑m=1Mw^mz,ηm,ωθ^mJ0ηmreiωt.The vertical stress can also be calculated by substituting equation (41) into equation (4) as follows:(42)σzr,z,t=c13∑m=1Mu^mz,ηm,ωθ^mηmJ0ηmreiωt+c33∑m=1M∂u^mz,ηm,ω∂zθ^mJ1ηmreiωt.


## 6. Numerical Example Verification

### 6.1. Crossanisotropic Layered Mediums

A multilayered medium model constructed by Khojasteh et al. [27] was used as a numerical example to verify the applicability of the algorithm to crossanisotropic mediums. The parameters of the layered medium are listed in Table 1. The medium includes three homogeneous layers, and the bottom is a homogeneous half-space medium. Layer 3 is an isotropic medium, and the other layers and the half-space medium are crossanisotropic. A vertical, uniformly distributed, harmonic load with an amplitude (p) was applied to a disk with a radius (Δr) of 0.0705h centered at the point (0, 0, 4h), and the displacement and stress responses of the crossanisotropic layered medium at r = 0 and at different depths were calculated at a dimensionless frequency (ω0) of 2 and 4. The frequency was made dimensionless using equation ω0=ωh/Gvh1/ρ1, where Gvh1 and ρ1 are the elastic parameter and the mass density of the first layer, respectively. The dimensionless depth was also defined as z/Δr. The calculated dimensionless vertical displacement and dimensionless stress were expressed in Gvh1wh/p and σzh2/p, respectively.Table 1
The material parameters of the crossanisotropic layered medium.
Layerh (m)Ehh (MPa)Evh (MPa)μhh = μvhGvh (MPa)ρ (kg/m3)1h2.53.00.251.001.022h3.04.00.251.401.134h5.05.00.252.001.3Half-space—7.56.00.252.501.5The calculated dimensionless vertical displacement and dimensionless vertical stress are delineated in Figures4 and 5. Moreover, the comparison of the calculation data with the results of Khojasteh et al. [27] verifies the accuracy of the algorithm for solving the dynamic response of isotropic layered mediums.Figure 4
The real and imaginary parts of the dimensionless vertical displacement at origin versus the dimensionless depth.Figure 5
The real and imaginary parts of the dimensionless vertical stress at origin versus the dimensionless depth.


### 6.2. Verifying Models of Contact between Different Layers

To verify the correctness of the calculation method proposed herein, the layered medium model proposed by Liu and Pan [37] was used as an example to calculate the displacement response of the crossanisotropic layered medium when incomplete contact between the layers was considered, and the results were analyzed and compared with those obtained from the dual variable and position (DVP) method. A double-layered half-space medium model was also constructed as shown in Figure 6, and we assumed that the upper crossanisotropic layered medium was over the homogeneous elastic half-space body. The material properties of each layer are also presented in Figure 6. The mass density of all the materials is equal to 5000 kg/m3, and the damping ratio equals 5%; the thickness of the upper medium is h = 1.0 m. Three models with different interface contact parameters were used to calculate the displacement response; Model I: tr = 1 × 1010 N/m3 and tz = 3 × 1010 N/m3; Model II: tr = 1 × 1010 N/m3 and tz = ∞; Model III: tr = ∞ and tz = 3 × 1010 N/m3. A vertical, uniformly distributed, harmonic load with an amplitude of p=1.0 Pa was applied on the surface of the medium. The radius of the uniformly distributed load Δr = 1.0 m. The dimensionless frequency and the calculated dimensionless vertical displacement were expressed through ω0=ωh/c44/ρ and c44uzh/ρ, respectively. The displacement w at the surface of the medium with a distance from 0 to 5 m was calculated at a dimensionless frequency ω0=1.Figure 6
A crossanisotropic double-layered half-space medium with an imperfect interface.The calculation results are plotted in Figure7. Then, the displacement response at the surface of the medium, i.e., (x, y, z) = (0.5, 0, 0), at different load frequencies was calculated. In addition, our calculation results (see Figure 7) are in good agreement with the ones determined by the DVP algorithm in the existing literature [37], verifying that the algorithm developed in this work is suitable for crossanisotropic layered mediums when incomplete contact between the layers is considered. It is worth noting that the proposed method in this paper is much more efficient than the algorithm using an infinite integral. It takes only 0.35 seconds for each frequency point by adopting the proposed method, while the algorithm using infinite integrals does the same calculation in about 60 seconds.Figure 7
The vertical displacement of the medium calculated by imperfect interface models I, II, and III at a fixed surface point, i.e., (x, y, z) = (0.5, 0, 0): (a) as a function of the dimensionless distance at a fixed dimensionless frequency of one (ω0 = 1); (b) as a function of the dimensionless frequency.
(a)(b)


## 6.1. Crossanisotropic Layered Mediums

A multilayered medium model constructed by Khojasteh et al. [27] was used as a numerical example to verify the applicability of the algorithm to crossanisotropic mediums. The parameters of the layered medium are listed in Table 1. The medium includes three homogeneous layers, and the bottom is a homogeneous half-space medium. Layer 3 is an isotropic medium, and the other layers and the half-space medium are crossanisotropic. A vertical, uniformly distributed, harmonic load with an amplitude (p) was applied to a disk with a radius (Δr) of 0.0705h centered at the point (0, 0, 4h), and the displacement and stress responses of the crossanisotropic layered medium at r = 0 and at different depths were calculated at a dimensionless frequency (ω0) of 2 and 4. The frequency was made dimensionless using equation ω0=ωh/Gvh1/ρ1, where Gvh1 and ρ1 are the elastic parameter and the mass density of the first layer, respectively. The dimensionless depth was also defined as z/Δr. The calculated dimensionless vertical displacement and dimensionless stress were expressed in Gvh1wh/p and σzh2/p, respectively.Table 1
The material parameters of the crossanisotropic layered medium.
Layerh (m)Ehh (MPa)Evh (MPa)μhh = μvhGvh (MPa)ρ (kg/m3)1h2.53.00.251.001.022h3.04.00.251.401.134h5.05.00.252.001.3Half-space—7.56.00.252.501.5The calculated dimensionless vertical displacement and dimensionless vertical stress are delineated in Figures4 and 5. Moreover, the comparison of the calculation data with the results of Khojasteh et al. [27] verifies the accuracy of the algorithm for solving the dynamic response of isotropic layered mediums.Figure 4
The real and imaginary parts of the dimensionless vertical displacement at origin versus the dimensionless depth.Figure 5
The real and imaginary parts of the dimensionless vertical stress at origin versus the dimensionless depth.


## 6.2. Verifying Models of Contact between Different Layers

To verify the correctness of the calculation method proposed herein, the layered medium model proposed by Liu and Pan [37] was used as an example to calculate the displacement response of the crossanisotropic layered medium when incomplete contact between the layers was considered, and the results were analyzed and compared with those obtained from the dual variable and position (DVP) method. A double-layered half-space medium model was also constructed as shown in Figure 6, and we assumed that the upper crossanisotropic layered medium was over the homogeneous elastic half-space body. The material properties of each layer are also presented in Figure 6. The mass density of all the materials is equal to 5000 kg/m3, and the damping ratio equals 5%; the thickness of the upper medium is h = 1.0 m. Three models with different interface contact parameters were used to calculate the displacement response; Model I: tr = 1 × 1010 N/m3 and tz = 3 × 1010 N/m3; Model II: tr = 1 × 1010 N/m3 and tz = ∞; Model III: tr = ∞ and tz = 3 × 1010 N/m3. A vertical, uniformly distributed, harmonic load with an amplitude of p=1.0 Pa was applied on the surface of the medium. The radius of the uniformly distributed load Δr = 1.0 m. The dimensionless frequency and the calculated dimensionless vertical displacement were expressed through ω0=ωh/c44/ρ and c44uzh/ρ, respectively. The displacement w at the surface of the medium with a distance from 0 to 5 m was calculated at a dimensionless frequency ω0=1.Figure 6
A crossanisotropic double-layered half-space medium with an imperfect interface.The calculation results are plotted in Figure7. Then, the displacement response at the surface of the medium, i.e., (x, y, z) = (0.5, 0, 0), at different load frequencies was calculated. In addition, our calculation results (see Figure 7) are in good agreement with the ones determined by the DVP algorithm in the existing literature [37], verifying that the algorithm developed in this work is suitable for crossanisotropic layered mediums when incomplete contact between the layers is considered. It is worth noting that the proposed method in this paper is much more efficient than the algorithm using an infinite integral. It takes only 0.35 seconds for each frequency point by adopting the proposed method, while the algorithm using infinite integrals does the same calculation in about 60 seconds.Figure 7
The vertical displacement of the medium calculated by imperfect interface models I, II, and III at a fixed surface point, i.e., (x, y, z) = (0.5, 0, 0): (a) as a function of the dimensionless distance at a fixed dimensionless frequency of one (ω0 = 1); (b) as a function of the dimensionless frequency.
(a)(b)


## 7. Numerical Examples

### 7.1. Influence oftr and tz on Displacement and Stress Responses

The contact conditions of the interface between the adjacent layers affect the dynamic response of the layered medium. Thus, this section discusses the impact of the interface contact parameters, namely,tr and tz, on the dynamic responses of the crossanisotropic multilayered medium. The material properties of the layers of the crossanisotropic layered medium used for the solution analysis are listed in Table 2. The layered medium includes two homogeneous layers resting on a homogeneous elastic half-space. We assumed that the contact between layer 1 and layer 2 is incomplete, but the contact between the second layer and the homogeneous half-space medium is complete. The calculation examples were divided into six cases (Case 1 to Case 6) according to the values of tr and tz, and the interface contact parameters of the six cases are tabulated in Table 3. A vertical, uniformly distributed, harmonic load with an amplitude of 1.0 Pa was applied to the disk with a radius of 0.15 m on the surface of the medium. The displacement and stress responses of the medium at a distance (i.e., r) in the range of 0 to 5h from the center of the disc on the planes with a depth (i.e., z) of 0, 0.15, and 0.45 m were calculated at a dimensionless frequency of 2, i.e., ω0 = 2. The real and imaginary parts of the calculated dimensionless displacement and stress are delineated in Figures 8–10.Table 2
The material properties of the layers of the crossanisotropic layered medium.
Layerh (m)Ehh (MPa)Evh (MPa)μhhμvhGvh (MPa)ρ (kg/m3)β (%)10.15300035000.350.3012702400520.306007000.330.3025423002Half-space—40600.450.4018.519002Table 3
The interface contact parameters (tz and tr) of the six cases.
CaseCase 1Case 2Case 3Case 4Case 5Case 6tr (N/m3)3.2 × 1073.2 × 1083.2 × 1093.2 × 10103.2 × 10113.2 × 1012tz (N/m3)9.6 × 1079.6 × 1089.6 × 1099.6 × 10109.6 × 10119.6 × 1012Figure 8
The vertical displacement of the medium versus the dimensionless distance for the different cases at a fixed dimensionless frequency of 2 (ω0 = 2) and a depth of zero (z = 0 m): (a) the variation of the real part; (b) the variation of the imaginary part.
(a)(b)Figure 9
The vertical displacement and the vertical stress of the medium versus the dimensionless distance for the different cases at a fixed dimensionless frequency of two (ω0 = 2) and a depth of 0.15 m (z = 0.15 m): (a) and (c) are the variation of the real parts; (b) and (d) are the variation of the imaginary part.
(a)(b)(c)(d)Figure 10
The vertical displacement and the vertical stress of the medium versus the dimensionless distance for the different cases at a fixed dimensionless frequency of two (ω0 = 2) and a depth of 0.45 m (z = 0.45 m): (a) and (c) are the variation of the real parts; (b) and (d) are the variation of the imaginary parts.
(a)(b)(c)(d)Figure8 depicts the displacement response of the medium at a distance within the range of 0 to 5h from the center of the disc when the load is applied on the surface of the layered medium. It can be seen that the peak of the displacement curve is generally higher when the contact parameters are small, and at r = 0, only the real part of Case 2 is larger than that of Case 1. The peak values of the real and imaginary parts of the displacement curves gradually decrease with an increase in the contact parameters. The curves of Cases 4–6 basically overlap, that is, the cases in which the contact parameters are greater than or equal to those in Case 4 can simulate complete contact. Indeed, the contact between the layers of the medium improves further, so the first layer of the medium is more constrained at the bottom, resulting in a gradual decrease in the peak of the displacement curve.Figure9 illustrates the dynamic displacement and stress responses of the surface of the medium when the contact between the layers of the medium is incomplete. It is clear that the displacement and stress amplitude of the layer tends to enlarge as the contact parameters increase; that is, the influence of the load of the layered medium on the bottom of the incomplete contact layer gradually intensifies as the contact between the layers of the medium further improves.Figure10 shows the displacement and stress responses of the medium at a depth of 0.45 m. It is obvious that the load on the surface of the layered medium has less influence on the inner layers of the medium when the contact parameters are small, i.e., in Case 1. Moreover, the displacement and the stress generated by the load applied on the plane increase greatly when the contact parameters enlarge, i.e., in Case 2. However, the displacement and stress responses of this layer tend to decline gradually as the contact parameters continue to increase because it is difficult for the load on the surface of the medium to generate a dynamic response to the inner layers of the medium through incomplete contact when the contact parameters are small. On the one hand, as the contact parameters increase, the impact of the load on the inner layers of the medium intensifies; on the other hand, as the contact parameters continue to increase, the integrity of the layered medium strengthens, and the dynamic response of the entire layered medium gradually declines under the same load.


### 7.2. Effect of the Application Point of Load on Displacement and Stress Responses

The dynamic response generated in the layered medium is different when the load of the same magnitude was applied to the different layers of the layered medium. The model of the three-layer medium with a homogeneous half-space medium at the bottom described in Section6.1 was used to analyze the influence of the load applied to the different depths of the layered medium on the distribution of the displacement and stress responses. To this end, a vertical, uniformly distributed, harmonic load was applied to the different depths of the layered medium at r = 0. Specifically, the load was applied within the range of the disc with a radius of Δr = 0.0705h and at a depth of 0, 1h, 2h, 3h, 4h, and 5h, respectively. The displacement and stress responses of the layered medium were calculated at a dimensionless frequency of two (ω0 = 2.0). Figure 11 delineates the distributions of the displacement and stress responses of the medium in a depth range of 0 to 7h at r = 0. It can be seen that the real part of the displacement of the medium at the position where the load is applied is singular due to the small radius of the uniform load, while its imaginary part is finite and continuous when the load is applied. As the depth to which the load is applied increases, the displacement response of the medium surface gradually decreases; that is, the effect of the load on the medium surface gradually lessens. The distribution of the stress response of the medium at r = 0 is similar to that of the displacement of the medium. The real part of the stress response is also singular when the load is applied, but the real parts of the stress on the upper and lower surfaces of the plane to which the load is applied are in opposite numerical directions. The imaginary part is also limited and continuous when the load is applied to the plane and equals zero when the load is applied to the ground surface. As the depth to which the load is applied increases, the first peak of the distribution curve along the vertical direction gradually declines.Figure 11
The real and imaginary parts of the vertical displacement at (a)r = 0 and (b) r = h versus the dimensionless depth at a fixed dimensionless frequency of two (ω0 = 2).
(a)(b)Figure12 depicts the distributions of the vertical displacement and the vertical stress of the medium at r = h. According to the displacement distribution, the real part of the displacement is no longer singular on the plane to which the load is applied when there is a certain horizontal distance from the position where the load is applied; however, there is a kink, and the imaginary part is still finite and continuous. As the depth to which the load is applied increases, the impact of the load on the surface of the layered medium gradually softens. The distribution of the stress response along the vertical direction at r = h is similar to that at r = 0.Figure 12
The real and imaginary parts of the vertical stress at (a)r = 0 and (b) r = h versus the dimensionless depth at a fixed dimensionless frequency of two (ω0 = 2).
(a)(b)


### 7.3. Influence of Frequency on Displacement and Stress Responses

This section considers the influence of the interlayer contact on the dynamic response of the medium using the model of the three-layered medium with a homogeneous half-space body at the bottom described in Section6.1 to explore the impact of the frequency on the dynamic response of the layered media. The parameters of the layered medium are listed in Table 1. Moreover, we assume that the contact between layer 1 and layer 2 is incomplete; that is, there is an imperfect interface between the two layers with tr = 1 × 1010 N/m3 and tz = 3 × 1010 N/m3. The contact between the remaining layers and the adjacent layered medium or the uniform half-space body is complete. A vertical, uniformly distributed, harmonic load with an amplitude of 1.0 Pa was applied to the surface of the medium within the range of the disc with a radius of 0.0705h. The displacement and stress responses at different depths and r = 0 were calculated at a dimensionless frequency of 0, 4, 6, and 10.Figure13 displays the distributions of the vertical displacement and the vertical stress in a depth range of 0 to 2h at r = 0. It can be seen that the displacement response of the surface of the medium gradually decreases as the frequency rises, and the positions where the first peaks of the real part and the imaginary part of the displacement emerge gradually move up. The real and imaginary parts of the displacement at a depth of h (z = h) are shifted due to the incomplete contact between layer 1 and layer 2, while the real and imaginary parts of the stress are continuous. The displacement response at different load frequencies approaches zero as the depth to which the load is applied increases. In a depth range of 0–h, the rate of change in the displacement-depth curve declines with an increase in the depth to which the load is applied, and the curve gradually flattens. In a depth range of h to 2h, the displacement of the medium gradually approaches zero. The rate of change in the displacement-depth curve is low, and the curve approximates to a straight line; this indicates that when there is an incomplete contact between the layers of a layered medium, the load on its surface has a small effect on the stress and displacement responses of the medium below the depth of the imperfect interface; this also implies that the frequency has a negligible impact on the stress and displacement of the medium.Figure 13
The real and imaginary parts of the vertical displacement and the vertical stress versus the dimensionless depth at different dimensionless load frequencies.The distribution of the vertical stress response of the layered medium in a depth range of 0–2h at r = 0 is similar to that of the vertical displacement of the medium, and the medium surface has a similar stress response at different load frequencies. The curve of the real part of the stress still oscillates with a small amplitude in a depth range of h to 2h when the dimensionless frequency equals 10, which is because the displacement of the adjacent layers of the medium at the imperfect interface is discontinuous, while the stress on the medium is still continuous. The imaginary part of the stress on the medium is equal to zero when the dimensionless frequency is 10. As the frequency increases, the peak of the imaginary part of the stress rises, and the depth of the point where the peak appears is different when the dimensionless frequency equals 4, 6, and 10.


### 7.4. Influence of Heterogeneous Properties of Mediums on Displacement and Stress Responses

Different layered media have various heterogeneous properties which affect their dynamic response. This section uses the model of the three-layered medium with the homogeneous half-space body at the bottom described in Section6.1 to analyze the displacement and stress responses of layered mediums with different heterogeneous properties under the same load. The parameters of the medium are presented in Table 1. The medium was divided into four cases (Cases 1–4) according to the heterogeneous properties of the medium; the properties of the medium in Case 1, that is, a three-layer medium (layers 1–3) with a homogeneous half-space body (the half-space) at the bottom, are listed in Table 1. Case 2 is a double-layer medium (layers 1 and 2) with a homogeneous half-space body at the bottom and the properties of which are the same as those of layer 3 in Case 1 (see Table 1). Case 3 is a single-layer medium (layer 1) with a homogeneous half-space body at the bottom; the properties of the half-space medium are the same as those of layer 2 in Case 1. Case 4 is composed of a crossanisotropic homogeneous half-space medium, and its properties are the same as those of layer 1. A vertical, uniformly distributed, harmonic load with an amplitude of 1.0 Pa was applied to the surface of the medium within the range of a disc and Δr = 1.0 m; the displacement and stress responses at the surface were obtained at r values of 0, 0.5h, and h in a dimensionless frequency range of 0–10. The calculation results are depicted in Figures 14–16.Figure 14
The real part and imaginary parts of (a) the vertical displacement and (b) the vertical stress at a fixed surface point (r, θ, z) = (0, 0, 0) for the different cases versus the dimensionless frequency.
(a)(b)Figure 15
The real and imaginary parts of (a) the vertical displacement and (b) the vertical stress at a fixed surface point (r, θ, z) = (0.5h, 0, 0) for the different cases versus the dimensionless frequency.
(a)(b)Figure 16
The real and imaginary parts of (a) the vertical displacement and (b) the vertical stress at a fixed surface point (r, θ, z) = (h, 0, 0) for the different cases versus the dimensionless frequency.
(a)(b)Figure14 delineates the distributions of the vertical displacement of and the vertical stress on the surface of the layered medium at r = 0 in a dimensionless frequency range of 0–10. It can be seen that the real and imaginary parts of the displacement tend to decline as the frequency rises in the range of 0–10. The curves of the displacement response in Cases 1–4 gradually coincide when the frequency is in the high range of 6–10, indicating that the impact of the heterogeneous properties on the dynamic response of the layered medium gradually softens as the frequency enlarges and is small within a high-frequency band.In addition, among Cases 1–4, the curves of the displacement response in Cases 2 and 3 coincide, indicating that the difference in the heterogeneous properties of the layered medium in Cases 2 and 3 has a small influence on the dynamic response of the medium. Case 4 has a quite different curve of the displacement response compared to the other cases, implying that the heterogeneous properties of the medium have a significant effect on its dynamic response. The dimensionless frequency at which the peaks of the real and imaginary parts of the displacement of the medium appear is smaller in Case 4 than in Cases 1–3, and the peak of the imaginary part of the displacement of the medium is obviously larger in Case 4 than in Cases 1–3.Figures15 and 16 illustrate the distributions of the vertical displacement of and the vertical stress on the surface of the layered medium in a dimensionless frequency range of 0–10 at r values of 0.5h and h, respectively. The trend of the variation in the displacement and stress responses of the surface of the medium at r values of 0.5h and h is similar to the one at r = 0.


## 7.1. Influence oftr and tz on Displacement and Stress Responses

The contact conditions of the interface between the adjacent layers affect the dynamic response of the layered medium. Thus, this section discusses the impact of the interface contact parameters, namely,tr and tz, on the dynamic responses of the crossanisotropic multilayered medium. The material properties of the layers of the crossanisotropic layered medium used for the solution analysis are listed in Table 2. The layered medium includes two homogeneous layers resting on a homogeneous elastic half-space. We assumed that the contact between layer 1 and layer 2 is incomplete, but the contact between the second layer and the homogeneous half-space medium is complete. The calculation examples were divided into six cases (Case 1 to Case 6) according to the values of tr and tz, and the interface contact parameters of the six cases are tabulated in Table 3. A vertical, uniformly distributed, harmonic load with an amplitude of 1.0 Pa was applied to the disk with a radius of 0.15 m on the surface of the medium. The displacement and stress responses of the medium at a distance (i.e., r) in the range of 0 to 5h from the center of the disc on the planes with a depth (i.e., z) of 0, 0.15, and 0.45 m were calculated at a dimensionless frequency of 2, i.e., ω0 = 2. The real and imaginary parts of the calculated dimensionless displacement and stress are delineated in Figures 8–10.Table 2
The material properties of the layers of the crossanisotropic layered medium.
Layerh (m)Ehh (MPa)Evh (MPa)μhhμvhGvh (MPa)ρ (kg/m3)β (%)10.15300035000.350.3012702400520.306007000.330.3025423002Half-space—40600.450.4018.519002Table 3
The interface contact parameters (tz and tr) of the six cases.
CaseCase 1Case 2Case 3Case 4Case 5Case 6tr (N/m3)3.2 × 1073.2 × 1083.2 × 1093.2 × 10103.2 × 10113.2 × 1012tz (N/m3)9.6 × 1079.6 × 1089.6 × 1099.6 × 10109.6 × 10119.6 × 1012Figure 8
The vertical displacement of the medium versus the dimensionless distance for the different cases at a fixed dimensionless frequency of 2 (ω0 = 2) and a depth of zero (z = 0 m): (a) the variation of the real part; (b) the variation of the imaginary part.
(a)(b)Figure 9
The vertical displacement and the vertical stress of the medium versus the dimensionless distance for the different cases at a fixed dimensionless frequency of two (ω0 = 2) and a depth of 0.15 m (z = 0.15 m): (a) and (c) are the variation of the real parts; (b) and (d) are the variation of the imaginary part.
(a)(b)(c)(d)Figure 10
The vertical displacement and the vertical stress of the medium versus the dimensionless distance for the different cases at a fixed dimensionless frequency of two (ω0 = 2) and a depth of 0.45 m (z = 0.45 m): (a) and (c) are the variation of the real parts; (b) and (d) are the variation of the imaginary parts.
(a)(b)(c)(d)Figure8 depicts the displacement response of the medium at a distance within the range of 0 to 5h from the center of the disc when the load is applied on the surface of the layered medium. It can be seen that the peak of the displacement curve is generally higher when the contact parameters are small, and at r = 0, only the real part of Case 2 is larger than that of Case 1. The peak values of the real and imaginary parts of the displacement curves gradually decrease with an increase in the contact parameters. The curves of Cases 4–6 basically overlap, that is, the cases in which the contact parameters are greater than or equal to those in Case 4 can simulate complete contact. Indeed, the contact between the layers of the medium improves further, so the first layer of the medium is more constrained at the bottom, resulting in a gradual decrease in the peak of the displacement curve.Figure9 illustrates the dynamic displacement and stress responses of the surface of the medium when the contact between the layers of the medium is incomplete. It is clear that the displacement and stress amplitude of the layer tends to enlarge as the contact parameters increase; that is, the influence of the load of the layered medium on the bottom of the incomplete contact layer gradually intensifies as the contact between the layers of the medium further improves.Figure10 shows the displacement and stress responses of the medium at a depth of 0.45 m. It is obvious that the load on the surface of the layered medium has less influence on the inner layers of the medium when the contact parameters are small, i.e., in Case 1. Moreover, the displacement and the stress generated by the load applied on the plane increase greatly when the contact parameters enlarge, i.e., in Case 2. However, the displacement and stress responses of this layer tend to decline gradually as the contact parameters continue to increase because it is difficult for the load on the surface of the medium to generate a dynamic response to the inner layers of the medium through incomplete contact when the contact parameters are small. On the one hand, as the contact parameters increase, the impact of the load on the inner layers of the medium intensifies; on the other hand, as the contact parameters continue to increase, the integrity of the layered medium strengthens, and the dynamic response of the entire layered medium gradually declines under the same load.


## 7.2. Effect of the Application Point of Load on Displacement and Stress Responses

The dynamic response generated in the layered medium is different when the load of the same magnitude was applied to the different layers of the layered medium. The model of the three-layer medium with a homogeneous half-space medium at the bottom described in Section6.1 was used to analyze the influence of the load applied to the different depths of the layered medium on the distribution of the displacement and stress responses. To this end, a vertical, uniformly distributed, harmonic load was applied to the different depths of the layered medium at r = 0. Specifically, the load was applied within the range of the disc with a radius of Δr = 0.0705h and at a depth of 0, 1h, 2h, 3h, 4h, and 5h, respectively. The displacement and stress responses of the layered medium were calculated at a dimensionless frequency of two (ω0 = 2.0). Figure 11 delineates the distributions of the displacement and stress responses of the medium in a depth range of 0 to 7h at r = 0. It can be seen that the real part of the displacement of the medium at the position where the load is applied is singular due to the small radius of the uniform load, while its imaginary part is finite and continuous when the load is applied. As the depth to which the load is applied increases, the displacement response of the medium surface gradually decreases; that is, the effect of the load on the medium surface gradually lessens. The distribution of the stress response of the medium at r = 0 is similar to that of the displacement of the medium. The real part of the stress response is also singular when the load is applied, but the real parts of the stress on the upper and lower surfaces of the plane to which the load is applied are in opposite numerical directions. The imaginary part is also limited and continuous when the load is applied to the plane and equals zero when the load is applied to the ground surface. As the depth to which the load is applied increases, the first peak of the distribution curve along the vertical direction gradually declines.Figure 11
The real and imaginary parts of the vertical displacement at (a)r = 0 and (b) r = h versus the dimensionless depth at a fixed dimensionless frequency of two (ω0 = 2).
(a)(b)Figure12 depicts the distributions of the vertical displacement and the vertical stress of the medium at r = h. According to the displacement distribution, the real part of the displacement is no longer singular on the plane to which the load is applied when there is a certain horizontal distance from the position where the load is applied; however, there is a kink, and the imaginary part is still finite and continuous. As the depth to which the load is applied increases, the impact of the load on the surface of the layered medium gradually softens. The distribution of the stress response along the vertical direction at r = h is similar to that at r = 0.Figure 12
The real and imaginary parts of the vertical stress at (a)r = 0 and (b) r = h versus the dimensionless depth at a fixed dimensionless frequency of two (ω0 = 2).
(a)(b)


## 7.3. Influence of Frequency on Displacement and Stress Responses

This section considers the influence of the interlayer contact on the dynamic response of the medium using the model of the three-layered medium with a homogeneous half-space body at the bottom described in Section6.1 to explore the impact of the frequency on the dynamic response of the layered media. The parameters of the layered medium are listed in Table 1. Moreover, we assume that the contact between layer 1 and layer 2 is incomplete; that is, there is an imperfect interface between the two layers with tr = 1 × 1010 N/m3 and tz = 3 × 1010 N/m3. The contact between the remaining layers and the adjacent layered medium or the uniform half-space body is complete. A vertical, uniformly distributed, harmonic load with an amplitude of 1.0 Pa was applied to the surface of the medium within the range of the disc with a radius of 0.0705h. The displacement and stress responses at different depths and r = 0 were calculated at a dimensionless frequency of 0, 4, 6, and 10.Figure13 displays the distributions of the vertical displacement and the vertical stress in a depth range of 0 to 2h at r = 0. It can be seen that the displacement response of the surface of the medium gradually decreases as the frequency rises, and the positions where the first peaks of the real part and the imaginary part of the displacement emerge gradually move up. The real and imaginary parts of the displacement at a depth of h (z = h) are shifted due to the incomplete contact between layer 1 and layer 2, while the real and imaginary parts of the stress are continuous. The displacement response at different load frequencies approaches zero as the depth to which the load is applied increases. In a depth range of 0–h, the rate of change in the displacement-depth curve declines with an increase in the depth to which the load is applied, and the curve gradually flattens. In a depth range of h to 2h, the displacement of the medium gradually approaches zero. The rate of change in the displacement-depth curve is low, and the curve approximates to a straight line; this indicates that when there is an incomplete contact between the layers of a layered medium, the load on its surface has a small effect on the stress and displacement responses of the medium below the depth of the imperfect interface; this also implies that the frequency has a negligible impact on the stress and displacement of the medium.Figure 13
The real and imaginary parts of the vertical displacement and the vertical stress versus the dimensionless depth at different dimensionless load frequencies.The distribution of the vertical stress response of the layered medium in a depth range of 0–2h at r = 0 is similar to that of the vertical displacement of the medium, and the medium surface has a similar stress response at different load frequencies. The curve of the real part of the stress still oscillates with a small amplitude in a depth range of h to 2h when the dimensionless frequency equals 10, which is because the displacement of the adjacent layers of the medium at the imperfect interface is discontinuous, while the stress on the medium is still continuous. The imaginary part of the stress on the medium is equal to zero when the dimensionless frequency is 10. As the frequency increases, the peak of the imaginary part of the stress rises, and the depth of the point where the peak appears is different when the dimensionless frequency equals 4, 6, and 10.


## 7.4. Influence of Heterogeneous Properties of Mediums on Displacement and Stress Responses

Different layered media have various heterogeneous properties which affect their dynamic response. This section uses the model of the three-layered medium with the homogeneous half-space body at the bottom described in Section6.1 to analyze the displacement and stress responses of layered mediums with different heterogeneous properties under the same load. The parameters of the medium are presented in Table 1. The medium was divided into four cases (Cases 1–4) according to the heterogeneous properties of the medium; the properties of the medium in Case 1, that is, a three-layer medium (layers 1–3) with a homogeneous half-space body (the half-space) at the bottom, are listed in Table 1. Case 2 is a double-layer medium (layers 1 and 2) with a homogeneous half-space body at the bottom and the properties of which are the same as those of layer 3 in Case 1 (see Table 1). Case 3 is a single-layer medium (layer 1) with a homogeneous half-space body at the bottom; the properties of the half-space medium are the same as those of layer 2 in Case 1. Case 4 is composed of a crossanisotropic homogeneous half-space medium, and its properties are the same as those of layer 1. A vertical, uniformly distributed, harmonic load with an amplitude of 1.0 Pa was applied to the surface of the medium within the range of a disc and Δr = 1.0 m; the displacement and stress responses at the surface were obtained at r values of 0, 0.5h, and h in a dimensionless frequency range of 0–10. The calculation results are depicted in Figures 14–16.Figure 14
The real part and imaginary parts of (a) the vertical displacement and (b) the vertical stress at a fixed surface point (r, θ, z) = (0, 0, 0) for the different cases versus the dimensionless frequency.
(a)(b)Figure 15
The real and imaginary parts of (a) the vertical displacement and (b) the vertical stress at a fixed surface point (r, θ, z) = (0.5h, 0, 0) for the different cases versus the dimensionless frequency.
(a)(b)Figure 16
The real and imaginary parts of (a) the vertical displacement and (b) the vertical stress at a fixed surface point (r, θ, z) = (h, 0, 0) for the different cases versus the dimensionless frequency.
(a)(b)Figure14 delineates the distributions of the vertical displacement of and the vertical stress on the surface of the layered medium at r = 0 in a dimensionless frequency range of 0–10. It can be seen that the real and imaginary parts of the displacement tend to decline as the frequency rises in the range of 0–10. The curves of the displacement response in Cases 1–4 gradually coincide when the frequency is in the high range of 6–10, indicating that the impact of the heterogeneous properties on the dynamic response of the layered medium gradually softens as the frequency enlarges and is small within a high-frequency band.In addition, among Cases 1–4, the curves of the displacement response in Cases 2 and 3 coincide, indicating that the difference in the heterogeneous properties of the layered medium in Cases 2 and 3 has a small influence on the dynamic response of the medium. Case 4 has a quite different curve of the displacement response compared to the other cases, implying that the heterogeneous properties of the medium have a significant effect on its dynamic response. The dimensionless frequency at which the peaks of the real and imaginary parts of the displacement of the medium appear is smaller in Case 4 than in Cases 1–3, and the peak of the imaginary part of the displacement of the medium is obviously larger in Case 4 than in Cases 1–3.Figures15 and 16 illustrate the distributions of the vertical displacement of and the vertical stress on the surface of the layered medium in a dimensionless frequency range of 0–10 at r values of 0.5h and h, respectively. The trend of the variation in the displacement and stress responses of the surface of the medium at r values of 0.5h and h is similar to the one at r = 0.


## 8. Conclusions

In this paper, a new SEM was proposed to determine the dynamic response of the crossanisotropic layered mediums under a dynamic load. The incomplete contact between the adjacent layers of the multilayered medium is considered in the proposed algorithm. A double summation algorithm was utilized in this model instead of the infinite integral in the traditional algorithm. The efficiency of the inverse transformation has greatly improved. The numerical examples verified the accuracy and efficiency of the algorithm in this paper. Furthermore, an extensive parametric analysis of the contact parameters of the interface between the layers of the layered medium, namely, the point of application of the external load, the excitation frequency, and the heterogeneous properties of the layered medium, was conducted. The results demonstrate that(1)
The contact parameters of the interface between the layers of the layered medium have a significant influence on its dynamic response. The effect of the load on the dynamic response of the layered medium is significantly lessened when the contact parameters are small. However, when they gradually increase, the contact between the layers of the medium becomes closer, so the displacement of the upper and lower surfaces of the contact layer tends to be continuous.(2)
The dynamic response of the point at which the load is applied to the layered medium is also significant; indeed, the dynamic response near the load is correspondingly larger. As the vertical distance from the load rises, the impact of the load on the medium is gradually reduced.(3)
Although both the excitation load and the heterogeneous properties of the layered medium have a significant influence on the dynamic response of the layered medium, they act through different mechanisms.



---

*Source: 7697131-2022-12-10.xml*