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.
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 that n 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 , where ω is the angular frequency and can be expressed by the following equation:
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:where , , and are the normal stress in the r, θ, and z directions, respectively; is the shear stress in the plane ; ρ is the density of the medium; u and are the displacement components in r and z directions, respectively.
The displacements u and in the governing equations (2) can be expressed in the Helmholtz decomposition as follows:where and 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:where cij is the elastic constant of the crossanisotropic medium and , , , , , and χ are defined as follows:where and are Young’s modulus of the crossanisotropic medium in the horizontal and vertical directions, respectively; and represent Poisson’s ratios on the horizontal and vertical planes, respectively; 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):where .
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.with
Similarly, by substituting equation (3) into equation (7), the following equation can be obtained after algebraic simplification.with
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, interface zj of the layered medium displayed in Figure 1 is assumed to be imperfect, so the relationship between the stress on layer z = zj (, ) and the displacement of layer z = zj (, ) and the contact parameters (tr and tz) of the interface is defined as follows:where and 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:
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 interface z = 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:
The sum of the M 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.where is the flexibility coefficient of the medium in the translated domain; can be calculated by the following equations:
Expressing S (r) by means of Fourier–Bessel series, we obtain
5. Spectral Axisymmetric Element Formulation
The solution for 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 Figure 2, 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.
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 vector and stress vector in the single-layered medium are defined as follows:
Then, the state equation between and in the frequency-wave number domain can be expressed as follows:where is the characteristic matrix of the layered medium; stands for the first-order partial derivative of # with respect to the coordinate z; is the identity matrix; the matrices , , , and are defined as follows:
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:
5.2. Considering Incomplete Contact
The relationship between the surface force and the displacement of the nodes on the upper and lower sides of layer j (j = 1, 2, …, n) in a single-layered medium can be expressed according to the following equation:
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:where represents a 2 × 2 identity matrix and is defined as follows:
Considering the incomplete contact, T (zj − 0) and U (zj − 0) are eliminated based onwhere
It is noted that if , 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 Figure 3, 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:
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:with the eigenvalues and eigenvectors , and,where the real parts of all elements are positive.
Introducing a vector,
Then, equation (22) becomeswhere and are integration constants.
Substituting equation (34) into equation (33) yields
For a homogeneous half-space, the displacements and stresses must remain finite, which leads to . Applying Equation (35) results in the boundary condition at the surface (z = 0) of the half-space.
5.4. Assembling Structural Stiffness Matrix
The stiffness matrix of the entire system () 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:or simplified towhere is the external load vector of the node at the depth of zi, and,
is the stiffness matrix of the jth layer defined by Equation (40), and its value can be obtained from the following equation:
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. is a symmetric band matrix, and its inverse matrix is in equation (18). Therefore, the radial and vertical displacements in the spatial domain can be obtained by substituting displacements and 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.
The vertical stress can also be calculated by substituting equation (41) into equation (4) as follows:
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 () 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 , where and 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 and , respectively.
The calculated dimensionless vertical displacement and dimensionless vertical stress are delineated in Figures 4 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.
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 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 and , respectively. The displacement at the surface of the medium with a distance from 0 to 5 m was calculated at a dimensionless frequency .
The calculation results are plotted in Figure 7. 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.
(a)
(b)
7. Numerical Examples
7.1. Influence of tr 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.
(a)
(b)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Figure 8 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.
Figure 9 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.
Figure 10 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 Section 6.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.
(a)
(b)
Figure 12 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.
(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 Section 6.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.
Figure 13 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.
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 Section 6.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.
(a)
(b)
(a)
(b)
(a)
(b)
Figure 14 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.
Figures 15 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.
Appendix
A. Constructing Potential Function
The scalar potential function ξ1 is first solved. The Fourier analysis of the second-order differential Equation (8) is performed, and the following formula can be obtained by transforming the computational domain from the time domain into the frequency domain.
The variables with the caret symbol (ˆ) are the ones in the frequency domain. The solution of Equation (A.1) can be obtained by the method of separation of variables. Since is a function of two orthogonal axes r and z, the solution may consist of two independent functions of r and z:
Substituting Equation (A.2) into Equation (A.1), both sides of the equation are divided by and set equal to the constant −η2. Then, the two independent ordinary differential equations of and can be obtained.
If s = ηr, Equation (A.3) can be reduced to the Bessel equation by using the chain rule.
The solutions for Equation (A.5) are the Bessel functions of the first kind J0 and the second kind Y0; Y0 approaches infinity when r = 0. However, the solution should be discarded as the vibration is finite at the origin. is expressed by the following equation:where A1 is an unknown constant and η is the wave number in the radial direction. It is also assumed that the amplitude of the wave can be ignored when r = R (far away from the source), then the nonzero solution to equation (A.6) can be expressed in the following way:
This condition can be satisfied at the infinitely many positive root αm of function J0. If we let ηr = αm, that is, η = ηm = αm/R, the m solutions in the radial direction can be expressed as follows:
The solution for in equation (A.4) can be expressed through an exponential equation:where A2m is a constant; ; ηpzm represents the compression wave number as follows:
Then, the solution to equation (8) can be expressed as follows:
The potential function can also be solved in a similar way by transforming equation (9) from the time domain into the frequency domain with the Fourier transform.
Since is also a function of two orthogonal axes, r and z, the solution to is composed of two independent functions as expressed in the following equation:
The expression of the solution (equation (A.13)) is substituted into equation (A.12), and both sides of the equation are divided by , then both sides of the equation are set equal to the constant −η2, and the two independent ordinary differential equations can be calculated by the following equations:
If s = ηr, then by use of the chain rule, equation (A.14) reduces to Bessel’s equation as follows:
Then, the solution to equation (A.12) can be obtained as follows:where Bmn is a constant; J1 indicates the first-order Bessel function of the first kind; ηszm represents the number of the vertical shear waves and is calculated by the following equation:
Symbols
| Am and Bm: | Unknown constant coefficients |
| cij: | Elastic stiffness coefficients |
| Ehh and Evh: | Young’s moduli of the crossanisotropic medium in horizontal and vertical directions, respectively |
| and : | Correlation matrix between the displacement and stress vectors |
| : | Shear modulus on the vertical plane |
| : | Flexibility coefficients of the medium in the wave-number domain |
| : | The characteristic matrix of the layered medium |
| h: | Thickness of the layer medium |
| J0 and Y0: | The Bessel functions of the first kind and the second kind |
| J1: | The first-order Bessel function of the first kind |
| , , , and : | Properties of coefficient matrices of the layered medium |
| m: | Ordinal number of the wavenumber |
| M: | Number of the wavenumbers of summation |
| n: | Number of the layers of the layered medium |
| : | Axisymmetric harmonic load |
| (i = 0, 1, 2, …, n): | Load vector at the depth of zi |
| R: | Radius of the wave vanishing |
| (r, θ, z): | The cylindrical coordinate system |
| : | Stiffness matrix of a single homogeneous layer |
| tr and tz: | Horizontal and vertical contact parameters |
| t = diag {−1/tr − 1/tz}: | Interface contact parameters matrix |
| (j = 0, 1, 2, …, n): | Unknown tractions vector at the depth of zj |
| u and : | Radial and vertical displacements |
| (j = 0, 1, 2, …, n): | Displacement vector at the depth of zj |
| (i = 0, 1, 2, …, n): | Displacement vector of node i |
| : | The variables of the dual equation |
| zj (j = 1, 2, …, n): | The bottom depth of the jth layer |
| ω: | Angular frequency |
| ω0: | Dimensionless angular frequency |
| ρ: | Mass density of the layered medium |
| : | Normal stress in i-direction |
| τzr: | Shear stress on z-r plane |
| : | Fourier–Bessel transform coefficient |
| and : | Potential functions of the Helmholtz decomposition |
| and : | Poisson’s ratio on the horizontal and vertical plane, respectively |
| ηpzm (m = 1, 2, …, M): | Wave number of the compression wave |
| ηszm (m = 1, 2, …, M): | Wave number of the shear wave |
| χ = Ehh/Evh: | Ratio between Young’s modulus in horizontal and vertical directions |
| Δr: | Radius of uniform harmonic load |
| Ψ = [Ψ11 Ψ12; Ψ21 Ψ22]: | The eigenvectors matrix |
| Λ = diag {λi − λi}: | The eigenvalues matrix |
| αm: | Positive roots of function J0 |
| β: | Material damping ratio. |
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China under grant no. 51508203.