Assessing the mechanical responses for anisotropic multi-layered medium under harmonic moving load by Spectral Element Method (SEM)

Applied Mathematical Modelling 67 (2019) 22–37

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Assessing the mechanical responses for anisotropic multi-layered medium under harmonic moving load by Spectral Element Method (SEM) Lingyun You a,b, Kezhen Yan a,∗, Nengyuan Liu a, Tingwei Shi a, Songtao Lv c a

College of Civil Engineering, Hunan University, Changsha 410082, PR China Department of Civil and Environmental Engineering, Michigan Technological University, Houghton, MI 49931, USA c School of Traffic and Transportation Engineering, Changsha University of Science & Technology, Changsha 410004, PR China b

a r t i c l e

i n f o

Article history: Received 15 February 2018 Revised 30 September 2018 Accepted 9 October 2018 Available online 19 October 2018 Keywords: Multi-layered medium Transverse isotropy Spectral Element Method (SEM) Buchwald potential function Inverse Fast Fourier Transform (IFFT)

a b s t r a c t The article is concerned with the mechanical responses of anisotropic multi-layered medium under harmonic moving load. An analytical solution for two-dimensional anisotropic multi-layered medium subjected to harmonic moving load is devoted via Spectral Element Method (SEM), while the anisotropic property is approximated as transverse isotropy. Starting with the constitutive equations of transversely isotropic body and the governing equations of motion based on the loading properties. The analytical spectral elements in the wavenumber domain are obtained according to the principle of wave superposition and Fourier transformation. Then, the spectral global stiffness matrix of the multi-layered medium is derived by assembling the nodded stiffness matrices of all layers depended on the different interlayer conditions between the adjacent layers, i.e. sliding and bonded. The corresponding analytical solutions are achieved by taking the Fourier series and Inverse Fast Fourier Transform (IFFT) algorithm. Finally, some examples are given to validate the accuracy of the proposed analytical solution, and to demonstrate the impact of both anisotropy, top layer thickness, interlayer conditions, and loading properties (velocity and natural frequency) on the mechanical response of the multi-layered medium. © 2018 Elsevier Inc. All rights reserved.

1. Introduction Some modern composites, advanced materials, and traditional materials can be regarded as a multi-layered medium, such as flexible/rigid pavements, railway, geo-cell layers, microelectronic devices, and protective coating [1–6]. Accurate analysis of the mechanical response in these multi-layered medium are of great practical interest in many scientific and engineering applications, for example, in the study of mechanical responses of structures to moving loads. A large number of models and methods of analysis for the problem have been proposed. Semi-analytic methods based on integral transformation are widely used by Cole et al. [7] to analyze the ground vibration induced by impact and moving load. De Barros and Luco applied the numerical methods, such as the boundary element method or the finite element method, to study the dynamic responses due to impact and moving load [8]. Beskou and Theodorakopoulos gave a comprehensive review on the dynamic response of homogenous or layered half-space subjected to a moving load [9]. It should be noted that the above studies are

∗

Corresponding author. E-mail addresses: [email protected] (L. You), [email protected] (K. Yan).

https://doi.org/10.1016/j.apm.2018.10.010 0307-904X/© 2018 Elsevier Inc. All rights reserved.

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

23

confined to the homogeneous isotropic and linearly elastic continuous medium. However, almost all materials used in multilayered mediums exhibit apparent anisotropic properties, where anisotropy is a prerequisite of being direction-dependent. As opposed to with isotropy property, for example, asphalt mixture is composed of differently shaped granular materials, and its inner structure showed distinct anisotropic properties [10,11]. Wang et al. illustrated the stress fields of anisotropic and isotropic pavement under wheel loading through analytical solution and finite-element simulation in several cases, and they discovered that the anisotropic properties of asphalt mixture have significant effect on the pavement behavior [12]. The anisotropic properties of materials could be approximated as transverse isotropy, thus, transversely isotropic properties of materials can be viewed as a special case of anisotropy [13,14]. Also, the transversely isotropic model is more reasonable to describe the deformation of a layered medium as compared to an isotropic linear-elastic model. The interlayer condition of multi-layered medium is also important to its mechanical properties [15]. It is well known that soft polymers are widely used as adhesives and interlayers in multi-layered mediums for protection purposes. Stenzler and Goulbourne studied the behavior of several polyacrylate adhesives with varying microstructure, and the results illuminate that the interlayer microstructure of polyacrylate has positive effects on the impact behavior of multi-layered polymers [16]. Chupin et al. calculated the mechanical behavior of a viscoelastic layered medium to a moving load when interlayer slip was considered [17]. These studies shown that the interlayer conditions between the adjacent structure layers of multi-layered medium have significant influence on the mechanical behaviors. Therefore, the interlayer conditions should be considered in the process of multi-layered medium mechanical behavior analysis. In addition, taking into consideration the railway and road traffics background, the vibrations induced by moving loading is important because of the increasing speed of road and railway traffics. In such cases, the Rayleigh wave speed may be less than 100 m/s and the load speed may exceed the Rayleigh wave speed in the future [18,19]. For the load characteristics, there are different cases that can be considered in the mechanical response analysis of multi-layered mediums. As demonstrated in, constant amplitude, time-dependent amplitude, and harmonically time varying amplitude. Chahour et al. reported an alternative approach of dispersion curves to analyze resonant phenomena in the context of wave propagation induced by a harmonic load moving over a railway track coupled to a multi-layered poroviscoelastic medium [20]. Rajapakse and Wang derived the Green’s functions for transversely isotropic-elastic multi-layered medium subjected to the harmonic moving load [21]. The mechanical problems about multi-layered medium under the harmonic moving load are related to the wave propagation in the medium. Stoneley and Synge studied wave propagation in a transversely isotropic medium [22,23]. Ai et al. reported the analytical solutions for the dynamic response of a transversely isotropic multi-layered half-plane under axisymmetric load and moving load by using the analytical layer-element method [24,25]. Literature mentioned-above presented that many numerical models were used to investigate the mechanical response of transversely isotropic multilayered half-space subjected to a harmonic moving load. However, there are only a few research studies that exist about the transversely isotropic problem by semi-analytic method. Thus, it is still necessary to study these transversely isotropic problems further. Recently, the Spectral Element Method (SEM) developed by Rizzi and Doyle [26] has been a popular method to analyze the dynamic problem. SEM elegantly combines the excellence of the finite element method, dynamic stiffness method, and spectral analysis method, which includes the principle of wave superposition, the discretization of spatial domain, the assemblage theory of traditional finite element and the exactness of the stiffness matrix method [27]. In this approach, the system is solved by double summation over the involved frequencies and the wave numbers, thus alleviating the inconvenience of the numerical implementation of infinite integration [28]. An analytical layer element is used to describe a single layer in the spectral element technique, which not only reduces the computational requirement dramatically, but also demonstrates the numerical efficiency and stability due to the absence of positive exponential functions. The negative exponential function solution can effectively overcome the numerical overflow caused by the positive exponential function solution [24]. In addition, it eliminates the inconvenience of the numerical evaluation of contour integration between zero and infinity. Meanwhile, SEM is utilized for the analysis of the mechanical response of multi-layered mediums [29,30]. Based on the aforementioned literature, the objective of this study is to examine the mechanical responses of a twodimensional multi-layered transversely isotropic medium subjected to a harmonic moving load. In the mathematical derivation, Buchwald potential function will be applied to express the displacement and stress equations that could be figured out by solving the potential functions, which simplifies the calculation process. The negative exponential function in the proposed approach will be used to avoid numerical overflow problems. In the processes of solving the stress and displacement for the Fourier transform domain. The Fourier series expansion and the Inverse Fast Fourier Transform (IFFT) algorithm will be used, which could significantly reduce the operation time and improve efficiency. In addition, the effect analysis of the anisotropy, top layer thickness, interlayer conditions, and load properties (velocity and natural frequency) will be applied in the numerical discussion of this study. 2. Formulation of the problem and basic equations In this study, the mechanical response of a multi-layered half medium model under harmonic moving load will be analyzed, as shown in Fig. 1. The transversely isotropic can be observed as one special case of anisotropic [13,14], the material used in each layer was assumed as transverse isotropy, where Ev , Evh and Gv are the Young’s modulus in the vertical direction, Young’s modulus in the horizontal direction and shear modulus in the planes normal to the plane of transverse isotropy, respectively. μv and μvh are the Poisson’s ratios in the vertical and horizontal direction, respectively. In addition,

24

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

Fig. 1. Schematic diagram of multi-layered medium.

considering the mechanical response of the multi-layered half medium model under harmonic moving load will be affected by material damping, the material damping ƞ was used in this model based on the complex damping principle [31], where ∗ = E (1 + 2η i ), E ∗ = E (1 + 2η i ), and G∗ = G (1 + 2η i ) are the Young’s modulus in the vertical direction, Young’s moduEhi vi vi hi vi vi lus in the horizontal direction, and the shear modulus in the planes normal to the plane of transverse isotropy, respectively, when the material damping was considered. The surface of the medium was subjected to a harmonic moving load F(x1 , t) in x1 direction, the load speed was c, and the length of loading area was 2l. It is noteworthy that the problem can be reduced to a plane problem when the medium was assumed as infinite in the plane normal to the plane of x1 -o-z1 . 2.1. Governing equations The relationships between strain and displacement for the transversely isotropic multi-layered medium could be presented in the transversely isotropic constitutive equations in the two-dimensional space [10], as shown in equation below.

⎧ ∂u ∂u ⎪ σx1 = C11 x1 + C13 z1 ⎪ ⎪ ∂ x ∂ z1 1 ⎪ ⎨ ∂ ux1 ∂ u z1 σz1 = C13 + C33 ∂x1 ∂ z1  ⎪ ⎪ ⎪ ux1 ∂ u z1 ∂ ⎪ ⎩τx1 z1 = C44 + ∂ z1 ∂ x1

(1)

where C11 = λn(1 − nμ2vh ), C13 = λnμvh (1 + μh ), C33 = λ(1 − μ2h ), C44 = G∗v , and λ is the lamé constant. λ can be presented as Eq. (2). μvh and μh are the Poisson’s ratios in the vertical and horizontal direction respectively, n is the transversely isotropic coefficient, and n = Evh /Ev .

λ = Ev∗ /[(1 + μh )(1 − μh − 2nμ2vh )]

(2)

In addition, the relationship of strain and displacement for the transversely isotropic multi-layered medium could also be presented in transversely geometric equations in the two-dimensional space, as shown in Eq. (3), while the equations of motion for the two-dimensional space of transversely isotropic media are shown in Eq. (4).

εx1 = ∂ ux1 /∂ x1 εz1 = ∂ uZ /∂ z1 γx1 Z = ∂ ux1 /∂ z1 + ∂ uz1 /∂ x1

(3)

⎧ 2 ⎪ ⎨ ∂ σ x 1 + ∂ τx 1 z 1 = ρ ∂ u x 1 ∂ x1 ∂ z1 ∂t2 2 ∂ σ ∂ τ u z1 ∂ z1 x1 z1 ⎪ ⎩ + =ρ ∂ z1 ∂ x1 ∂t2

(4)

2.2. Buchwald potential function

 The Buchwald potential functions and  were adopted as u = ∇ + ∇ ×  and ∇  = 0 in this study, where ∇ = ∂ 2 /∂ x21 + ∂ 2 /∂ z12 . Thus, the displacement equation demonstrated by Buchwald potential functions and  can be written as:

u x 1 = ∂ / ∂ x 1

u z1 = ∂  /∂ z 1

(5)

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

25

The stress equation depicted by Buchwald potential functions and  can be obtained by introducing Eqs. (3) and (5) into (1), as shown in equation below.

⎧ ∂ 2

∂ 2 ⎪ ⎨σz1 = C13 2 + C33 2 ∂x1 ∂ z1  2

∂ 2 ∂ ⎪ ⎩τx1 z1 = C44 + ∂ x1 ∂ z1 ∂ x1 ∂ z1

(6)

Substituting Eqs. (5) and (6) into Eq. (1), we can obtain:

⎧  ∂ 2

∂ 2

∂ 2 ∂ 2

⎪ ⎪ + C44 + (C13 + C44 ) =ρ ⎨ C11 2 2 ∂ x1 ∂ z1 ∂ z1  ∂t2  ∂ x21 2 2 ∂  ∂  ∂

∂ 2 ⎪ ⎪ ⎩ C33 2 + C44 2 + (C13 + C44 ) ∂ x ∂ z = ρ 2 ∂ z1 ∂ x1 ∂t 1 1

(7)

2.3. Coordinate transformation and Fourier transformation In order to simplify the calculation, the movement coordinate x-o-z will be introduced here, as shown in Fig. 1, while the coordinate x1 -o-z1 is the fix coordinate in the multi-layered medium. A harmonic moving load F(x1 , t) in x1 direction with the length of loading area was 2l, thus, the variables in the movement coordinate can be expressed as x = x1 -ct, z = z1 . The equations below can be obtained by applying the coordinate transformation for the vertical displacement (ux ), the horizontal displacement (uz ), and Buchwald potential functions and  .

⎧ iωt ⎪ ⎨ux1 (x1 − ct, z, t ) = ux (x, z, t )eiωt uz1 (x1 − ct, z, t ) = uz (x, z, t )e iωt ⎪ ⎩ (x1 − ct, z, t ) = φ (x, z, t )e iωt (x1 − ct, z, t ) = ψ (x, z, t )e

(8)

In addition, the first order and second order differential expressions of Buchwald potential functions and  with respect to time t and the differential expressions of stress and displacement on the coordinates x and z under the movement coordinate system are expressed as follows:

    ⎧ ∂

∂φ iwt ∂  ∂ψ iwt ⎪ ⎪ = i ωφ − c e = i ωψ − c e ⎪ ⎪ ∂t  ∂x ∂t  ∂x ⎪ ⎪ ⎪ ⎪ ∂ 2

∂φ ∂ 2φ ⎪ ⎪ = −ω2 φ − 2iωc + c2 2 eiwt ⎪ 2 ⎪ ∂t ∂x ∂x  ⎨  2 ∂ 2 ∂ψ 2 2∂ ψ = −ω ψ − 2iωc +c eiwt ⎪ 2 ⎪ ∂x ∂t ∂ x2 ⎪ ⎪ ⎪ ∂ ∂ ⎪ ⎪ (ux1 , σx1 , τx1 z1 , ,  ) = (ux , σx , τxz , ,  ) ⎪ ⎪ ∂ x ∂ x ⎪ 1 ⎪ ⎪ ∂ ∂ ⎩ (u , σ , τ , ,  ) = (uz , σz , τxz , ,  ) ∂ z 1 z1 z1 x1 z1 ∂z

(9)

where σ x , σ z , and τ xz are the stress components in the movement coordinate system. Substituting Eqs. (8) and (9) into Eq. (7), Eq. (10) can be obtained as:

⎧ 2 ∂2 ∂ ∂ 2ψ ⎪ 2 ∂ 2 ⎪ φ + (C13 + C44 ) 2 = 0 ( C − ρ c ) + ρω + C + 2 i ρω c 44 ⎨ 11 2 2 ∂x ∂x ∂z ∂z

2 2 ⎪ ∂ ∂ ∂ ∂ 2φ ⎪ ⎩ (C44 − ρ c2 ) 2 + ρω2 + C33 2 + 2iρωc ψ + (C13 + C44 ) 2 = 0 ∂x ∂x ∂z ∂x

(10)

2.3.1. Fourier transformation for Eq. (10) In order to realize the transformation of second order partial differential equation, i.e. Eq. (10), from the space domain to wave number domain, the Fourier transform was adopted, the transformation result is shown in equation below.

⎧ 2 ¯ 2 ¯ ⎪ ⎨θ 2 − k3 ξ 2 φ¯ + d φ + k1 d ψ = 0 d z2 d z2 2 ¯  ⎪ ⎩ θ 2 − ξ 2 ψ¯ + k2 d ψ − k1 ξ 2 φ¯ = 0 2

(11)

dz

where k1 = (C13 + C44 )/C44 , k2 = C33 /C44 , k3 = C11 /C44 , and ϑ2 = [ρ (ω − ξ c)2 ]/C44 . ξ is the Fourier parameter, φ¯ and ψ¯ are the Fourier transforms with respect to the and  . Set the exponential solution of Eq. (11) as φ¯ = Ae−ϕ z and ψ = Be−ϕ z , and

26

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

Fig. 2. Two-node spectral element.

substituting φ¯ and ψ into Eq. (11), the general solution of Eq. (11) based on the exponential formula can be presented as:

φ = Aβ1 e−ϕ1 z + Bβ2 e−ϕ2 z ψ = A e −ϕ 1 z + B e −ϕ 2 z

(12)

where χ (ξ ) = [ν (ξ /θ )2 − (1 + k2 )]2 − 4k2 [k3 (ξ /θ )4 − (1 + k3 )(ξ /θ )2 ], ν = 1 + k2 k3 − k3 2 , and ϕ1,2 =  √ θ [ν (ξ /θ )2 − (1 + k2 ) ± χ ]/2k2 . 2.3.2. Fourier transformation for Eqs. (5) and (6) Similarly, the Fourier transformation was adopted in Eqs. (5) and (6) based on the movement coordinate x-o-z, we obtained following equations:

u¯ x = iξ φ¯

(13-a)

dψ¯ dz

(13-b)

u¯ z =

∂ 2 ψ¯ ∂ z2  dψ¯

σ¯ z = −C13 ξ 2 φ¯ + C33  τ¯xz = C44 iξ

dφ¯ + dz dz

(13-c)

(13-d)

3. Analytical stiffness matrices for the multi-layered medium 3.1. Construction of the spectral element matrices 3.1.1. Two-node spectral element In the study of multi-layered mediums, we assume every structure layer (except the bottom layer) is a two-node spectral element, according to the concepts of finite element method. The two-node spectral element was schematically shown in Fig. 2. As shown in Fig. 2 and implied by its name, the two nodes of the spectral element are constructed at the center points of the upper and lower boundary surfaces of the structure layer. The length of the spectral element is considered infinite compared to its elemental thickness. The wave diffracted in the z-direction can be regarded as the superposition of a direct wave and a reflected wave, according to the principle of SEM. Therefore, the potential function in the spectral element can also be regarded as the superposition of the direct wave potential function and the reflected wave potential function:

φ¯ = Aβ1 e−ϕ1 z + Bβ1 e−ϕ1 (h−z) + C β2 e−ϕ2 z + Dβ2 e−ϕ2 (h−z) ψ¯ = Ae−ϕ1 z + Be−ϕ1 (h−z) + C e−ϕ2 z + De−ϕ2 (h−z)

(14)

where h is the thickness of the spectral element, while A, B, C, and D are four different coefficients. In Eq. (14), for both φ¯ and ψ , the first exponential function of the potential function represents the direct wave emitted from z = 0 and the second exponential function represents the reflected wave emitted from z = h. Substituting Eq. (14) into Eqs. (13-a) and (13-b), the horizontal and vertical displacements of the two-node spectral element can be expressed as Eq. (15), where the superscript j ( j) ( j) in u¯ x and u¯ z represents the sequence number of the node in the spectral element. The construction node in the two-node spectral element is located at the top and bottom. The condition z = 0 for the upper node j = 1 and the condition z = h for the bottom node j = 2 were substituted into Eq. (15). The displacement of two-node spectral element can be expressed as the dot production of a 4 × 4 stiffness matrix K¯1 and a coefficient vector (A, B, C, D)T , as shown in Eq. (16).

u¯ x( ) = Aiξ β1 e−ϕ1 z + Biξ β1 e−ϕ1 (h−z ) + Ciξ β2 e−ϕ2 z − Diξ β2 e−ϕ2 (h−z ) j u¯ z( ) = −Aϕ1 e−ϕ1 z + Bϕ1 e−ϕ1 (h−z ) − C ϕ2 e−ϕ2 z + Dϕ2 e−ϕ2 (h−z ) j

(15)

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

27

Fig. 3. One-node spectral element.

⎛

⎞

⎡

u¯ x( ) iξ β1 ⎜u¯ z(1) ⎟ ⎢ −ϕ1 ⎜ (2 ) ⎟ = ⎣ ⎝u¯ x ⎠ iξ β1 e−ϕ1 h (2 ) ϕ 1 e −ϕ 1 h u¯ z 1

iξ β1 e−ϕ1 h ϕ 1 e −ϕ 1 h iξ β1 −ϕ1

iξ β2 −ϕ2 iξ β2 e−ϕ2 h ϕ 2 e −ϕ 2 h

⎤⎛ ⎞

⎛ ⎞

iξ β2 e−ϕ2 h A A ϕ2 e−ϕ2 h ⎥⎜B ⎟ ¯ ⎜B ⎟ ⎦⎝C ⎠ = K1 ⎝C ⎠ iξ β2 D D −ϕ2

(16)

Similarly, the vertical stress and the shear stress can also be obtained by substituting Eq. (14) into Eqs. (13-c) and (13-d), as shown in equation below.

σ¯ x( j ) = Aς1 e−ϕ1 z + Bς1 e−ϕ1 (h−z) + C ς2 e−ϕ2 z + Dς2 e−ϕ2 (h−z) τ¯xz( j ) = −Aiξ δ1 e−ϕ1 z + Biξ δ1 e−ϕ1 (h−z) − Ciξ δ2 e−ϕ2 z + Diξ δ2 e−ϕ2 (h−z)

(17)

where ςi = (C33 ϕi2 − C13 β j ξ 2 ), δi = C44 ϕi (1 + βi ), i = 1, 2. Substituting the condition z = 0 for the upper node j = 1 and the condition z = h for the bottom node j = 2 into Eq. (17), the stress of the two-node spectral element can be written as the dot product of a 4 × 4 stiffness matrix K¯2 and a coefficient vector (A, B, C, D)T , as shown in equation below.

⎛ (1 ) ⎞ ⎡ σ¯ z ς1 ⎜τ¯xz(1) ⎟ ⎢ −iξ δ1 ⎜ (2 ) ⎟ = ⎣ ⎝σ¯ z ⎠ ς 1 e −ϕ 1 h (2 ) −iξ δ1 e−ϕ1 h τ¯xz

ς 1 e −ϕ 1 h iξ δ1 e−ϕ1 h ς1 iξ δ1

ς2 −iξ δ2 ς 2 e −ϕ 2 h −iξ δ2 e−ϕ2 h

⎤⎛ ⎞ ⎛ ⎞ ς 2 e −ϕ 2 h A A iξ δ2 e−ϕ2 h ⎥⎜B ⎟ ⎜B ⎟ ¯ = K2 ⎝ ⎠ C ς2 ⎦⎝C ⎠ D D iξ δ2

(18)

Owing to the finite element method and the Cauchy stress principle [32], the boundary tractions, stress, and displacement of the two-node spectral element can be presented as:

⎛

⎞

⎛

⎞

⎛

⎞

1 1 1 −σ¯ z( ) u¯ x( ) T¯z( ) ⎜T¯xz(1) ⎟ ⎜ −τ¯xz(1) ⎟ ⎜u¯ z(1) ⎟ ⎜ (2) ⎟ = ⎜ (2) ⎟ = K¯two−node ⎜ (2) ⎟ ⎝T¯z ⎠ ⎝ σ¯ z ⎠ ⎝u¯ x ⎠ 2 (2 ) (2 ) ¯ τ¯xz u¯ z( ) Txz

(19)

where K¯two−node = N K¯2 K¯1−1 , and N can be written as:

⎡

⎤

−1

⎢

⎥ ⎦

−1

N=⎣

1

(20)

1 3.1.2. One-node spectral element The one-node spectral element, which is the bottom layer, was schematically shown in Fig. 1. As shown in Fig. 3 and implied by its name, the one-node spectral element just has one single boundary at the top of the layer and extends infinitely in all other directions [33]. One-node spectral element is essentially an infinite element [34]. If the depth of the bottom layer in the multi-layered medium is sufficiently large, the damping of the layer-material itself will cause the wave to disappear before reaching the boundary and there is no reflected wave in the depth direction. Therefore, the one-node spectral element can be used to represent an infinite half-space media, i.e. the bottom layer of the multi-layered medium. The solution of the one-node potential function was represented as Eq. (12), which degenerated from Eq. (14) when the coefficient B and D are zero, since there is no reflected wave in the depth direction. Meanwhile, the coefficients B and D in Eqs. (16) and (18) are also zero. For the one-node spectral element, z = 0, so its stiffness matrix of displacement K¯3 and the stiffness matrix of stress K¯4 can be obtained as:



u¯ x( ) 1 u¯ z( ) 1





=

iξ β1 −ϕ1

 (1 )  ς1 σ¯ z = −iξ δ1 τ¯xz(1)

iξ β2 −ϕ2

 

ς2 −iξ δ2

A B

 

= K¯3

  A B

A B

(21)

 

A = K¯4 B

(22)

28

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

where ς j = (C33 ϕ 2j − C13 β j ξ 2 ), δ j = C44 ϕ j (1 + β j ), j = 1, 2. Owing to the Cauchy stress principle [32], the boundary tractions, stress, and displacement of the one-node spectral element can be presented as:



1 T¯z( ) 1 ( T¯xz )





=

−σ¯ z( ) 1 −τ¯xz( ) 1





= K¯one−node

u¯ x( ) 1 u¯ z( ) 1



(23)

where K¯one−node = −K¯4 K¯3−1 . 3.2. Construction of the global stiffness matrix The global stiffness matrix can be obtained by the matrix stacking. There are two different methods of assembly for the different interlayer condition, i.e. sliding condition and bonded condition. 3.2.1. Bonded condition Suppose the two-node spectral element stiffness matrix K¯ i and K¯ i+1 are:

⎡

ki11 i ⎢ki K = ⎣ 21 ki31 ki41

ki12 ki22 ki32 ki42

ki13 ki23 ki33 ki43

⎤

ki14 ki24 ⎥ ki34 ⎦ ki44

⎡

K

i+1

+1 ki11 i+1 ⎢k = ⎣ 21 +1 ki31 +1 ki41

+1 ki12 i+1 k22 +1 ki32 +1 ki42

+1 ki13 i+1 k23 +1 ki33 +1 ki43

⎤

+1 ki14 i+1 k24 ⎥ +1 ⎦ ki34 i+1 k44

(24)

The contact surface equilibrium equation when the interlayer condition was completely continuous, i.e. bonded condition, can be presented as Eq. (25). Meanwhile, the continuity deformation equation can be written as Eq. (26). The boundary condition in the contact surface when the contact surface is without boundary tractions was shown in Eq. (27).

⎧ i i i k ux1 + ki32 uiz1 + ki33 uix2 + ki34 uiz2 = Tz2 ⎪ ⎨ 31 i i i i i i i i i

k41 ux1 + k42 uz1 + k43 ux2 + k44 uz2 = Txz2 +1 i+1 +1 i+1 +1 i+1 i+1 + ki12 uz1 + ki13 ux2 + ki14 uz2 = Tz1 ⎪ki11+1 uix+1 1 ⎩ i+1 i+1 i+1 i+1 i+1 i+1 i+1 i+1 i+1 k21 ux1 + k22 uz1 + k23 ux2 + k24 uz2 = Txz1

(25)

uix2 = uix+1 1

+1 uiz2 = uiz1

(26)

i i+1 Tz2 = Tz2

i i+1 Txz2 = Txz2

(27)

Substituting Eqs. (25)–(27) into Eq. (24), the superposition stiffness matrix of two adjacent layers when the interlayer condition was completely continuous presented in Eq. (28), while Eq. (28) also can be simplified as Eq. (29) when the i + 1 layer as the one-node spectral element.

⎡

ki11 ⎢ki21 ⎢ i ⎢k31 ⎢ ki ⎢ 41 ⎣0 0

⎡

ki11 ⎢ki21 ⎣k i 31 ki41

⎤⎧ i ⎫ ⎧ i ⎫ ux1 ⎪ Tz1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ i i ⎪ ⎥⎪ ⎪ ⎪ uz1 ⎪ ⎪ ⎪ ⎪ ⎨ Txz1 ⎪ ⎬ i+1 ⎥⎨ i+1 ⎬ k14 ⎥ ux1 0 = i+1 ⎥ i+1 k24 ⎥⎪uz1 ⎪ ⎪ 0 ⎪ ⎪ i+1 ⎪ ⎪ i+1 ⎪ +1 ⎦⎪ ⎪ ⎪ ⎪ ⎪T ⎪ ⎪ ki34 u ⎪ ⎪ ⎩ x2 ⎪ ⎭ ⎩ z2i+1 ⎭

ki12 ki22 ki32 ki42 0 0

ki13 ki23 +1 ki33 + ki11 i+1 i k43 + k21 +1 ki31 +1 ki41

ki14 ki24 +1 ki34 + ki12 i+1 i k44 + k22 +1 ki32 +1 ki42

ki12 ki22 ki32 ki42

ki13 ki23 +1 ki33 + ki11 i+1 i k43 + k21

i ki14 ui ⎪ ⎨ x1 ⎪ ⎬ ⎪ ⎨Tz1i ⎪ ⎬ ki24 ⎥ uiz1 Txz1 = i+1 ⎦ i+1 i k34 + k12 ⎪ux1 ⎪ ⎪ 0 ⎪ ⎩ i+1 ⎭ ⎩ ⎭ +1 0 ki44 + ki22 uz1

0 0 +1 ki13 i+1 k23 +1 ki33 +1 ki43

⎤⎧

0 0

+1 ki44

⎫

⎧

+1 uiz2

(28)

Txz2

⎫

(29)

3.2.2. Sliding condition The continuity deformation equation was Eq. (30) when the interlayer condition was sliding, while boundary condition in the contact surface when the contact surface without boundary tractions can be obtained as below. +1 uiz2 = uiz1

i i+1 Tz2 = Tz1

(30) i i+1 Txz2 = 0 Txz1 =0

(31)

Substituting Eqs. (30) and (31) into Eqs. (24) and (25), the superposition stiffness matrix of two adjacent layers when the interlayer condition was sliding presented in Eq. (32), while Eq. (32) also can be simplified as Eq. (33) when the i + 1

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

layer as the one-node spectral element.

⎡

ki11 ⎢ki21 ⎢ki ⎢ 31 ⎢ki41 ⎢ ⎢0 ⎣ 0 0

⎡

ki11 ⎢ki21 ⎢ i ⎢k31 ⎣ ki 41 0

ki12 ki22 ki32 ki42 0 0 0

ki13 ki23 ki33 ki43 0 0 0

ki12 ki22 ki32 ki42 0

ki13 ki23 ki33 ki43 0

ki14 ki24 +1 ki34 + ki12 ki44 +1 ki22 i+1 k32 +1 ki42

0 0 +1 ki11 0 +1 si21 i+1 s31 +1 si41

⎤⎧ i ⎫ ⎧ i ⎫ u Tz1 ⎪ ⎪ ⎪ ⎪ uxi 1 ⎪ ⎪ ⎪ i ⎪ ⎪ ⎪ ⎪ T ⎥⎪ ⎪ ⎪ ⎪ z1 xz1 ⎪ ⎪ ⎪ +1 ⎥⎪ i ⎪ ⎪ ⎪ ⎪ ⎪ u ⎨ ⎬ ⎨ ⎬ ki14 0 2 ⎥ ix+1 0 0 ⎥ u = 1 ⎥ xi+1 ⎪ 0 ⎪ +1 ⎥⎪ ⎪ ⎪ ⎪ ki24 uz1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎦ i+1 ⎪ i+1 ⎪ i+1 ⎪ ⎪ ⎪ ⎪ Tz2 k34 ⎪ ux2 ⎪ ⎪ ⎪ ⎩ ⎭ ⎩ ⎭ i+1 i+1 i+1

0 0

0 0

+1 ki13 0 +1 ki23 i+1 k33 +1 ki43

k44

uz2

⎤⎧ i ⎫ ⎧ i ⎫ ki14 Tz1 ⎪ ⎪ ⎪ uxi 1 ⎪ ⎪ ⎪ i ⎪ ⎥⎪ ⎨ uz1 ⎪ ⎬ ⎪ ⎨Txz1 ⎬ ki24 i+1 ⎥ i i k34 + k12 ⎥ ux2 = 0 ⎪ ⎦⎪ ⎪ ⎪ ⎪ ⎪ 0 ⎪ ⎪ ki44 ui+1 ⎪ 1 ⎪ ⎩ xi+1 ⎭ ⎩ 0 ⎭ i+1

0 0 +1 ki11 0 +1 ki21

k21

29

(32)

Txz2

(33)

uz1

3.2.3. Boundary conditions and the global stiffness matrix format As shown in Fig. 1, the surface of the medium was subjected to a harmonic moving load F(x1 , t) in x1 direction, and the load speed was c, as well as the length of loading area was 2l. The response of each point in the multi-layered medium must be completely attenuated at the beginning of the next load cycle. Therefore, a sufficient time-interval T = 2L should be taken between two adjacent harmonic moving loads. The harmonic moving load acting on the medium upper surface, z = 0, at any time t can be expressed as Eq. (34), and the boundary condition of the upper surface can be written as Eq. (35).

P (x ) =

qeiωt /2l |x| ≤ l 0l ≤ |x| ≤ L

(34)

σz |z=0 = −P (x ) τxz |z=0 = 0

(35)

The multi-layered medium structure, as shown in Fig. 1, can be seen as the combination between a number of two-node spectral elements and a one-node spectral element. The global stiffness matrix of the multi-layered medium is composed of the spectral element stiffness matrix of each layer. According to the principle of combination of two adjacent stiffness matrices, i.e. bonded condition and sliding condition, the global stiffness matrix of the multi-layered medium in this study was represented by the combination of stiffness matrix of spectral elements in each layer as shown in equation below.

⎛

⎜ ⎜ ⎜ ⎜ K¯Global (z, ξ ) = ⎜ ⎜ ⎜ ⎝

1 K¯two−node





1 K¯two−node

⎞









j−1 K¯two−node j ¯ Ktwo−node

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(36)

Therefore, we can obtain the expression of boundary tractions, displacements, and global stiffness matrix of each node of multi-layered medium as follows:



K¯Global (z, ξ )

⎧ (1 ) ⎫ u¯ x ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ u¯ z(1) ⎪ ⎪ ⎪ ⎪ ⎪. ⎪ ⎪ ⎨ ⎬  .

⎪ ⎪ ⎪ . ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪u¯ x(n) ⎪ ⎪ ⎪ ⎩ (n ) ⎪ ⎭ u¯ z

=

⎧ (1 ) ⎫ P¯ ⎪ ⎪ ⎪ ⎪ ¯z(1) ⎪ ⎪ ⎪ P ⎪ ⎪ ⎪ ⎪ ⎪. xz ⎪ ⎪ ⎨ ⎬ .

⎪ ⎪ ⎪ . ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪P¯z(n) ⎪ ⎪ ⎪ ⎩ ¯ (n ) ⎪ ⎭

(37)

Pxz

3.2.4. Solutions in the space-domain and time-domain SEM uses the finite element space discretization theory to realize the transformation from the wavenumber domain to the space domain [35]. The computational efficiency can be improved using the discretization theory, meanwhile, the extremely high computational precision can also be maintained. In this study, Fourier series expansion was mainly used to transform the wavenumber domain to the space domain. Fourier series expansion can disperse the potential function, the global stiffness matrix, and the boundary conditions at 2L wavelength into 2m calculation points with the same spacing, where each calculation point was presented as wavenumber ξ m , m is discrete calculation order of points [36]. The solution of displacement and stress, i.e. Eqs. (15) and (17), were expressed in the form of series superposition, as shown in Eqs. (38) and (39). In addition, the global stiffness matrix can also be expressed in the form of superposition

30

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

series since both the displacement and stress can be expressed in the form of superposition series. Thus, the global stiffness matrix can be rewritten as Eq. (40).

⎧ ⎪ ⎨u¯ x = ⎪ ⎩u¯ z =

2M



m=1 2M



m=1

⎧ ⎪ ⎨σ¯ z =



2M m=1

−Am ϕ1m e−ϕ1m z + Bm ϕ1m e−ϕ1m (h−z ) − Cm ϕ2m e−ϕ2 z + Dm ϕ2m e−ϕ2m (h−z )



2M m=1 2M

⎪ ⎩τ¯xz =

Am iξm β1m e−ϕ1m z + Bm iξm β1m e−ϕ1m (h−z ) + Cm iξm β2m e−ϕ2m z − Dm iξm β2m e−ϕ2m (h−z )

Am ς1m e−ϕ1m z + Bm ς1 e−ϕ1m (h−z ) + Cm ς2m e−ϕ2m z + Dm ς2 e−ϕ2m (h−z )



m=1

K¯Global (z, ξm )

.

⎪ ⎪ ⎪ . ⎪ ⎪ ⎪ ⎪ (n ) ⎪ ⎪ ⎪ ⎪ ¯ u x ⎪ ⎪ ⎩ (n ) ⎪ ⎭

=

u¯ z

(38)





−Am iξm δ1m e−ϕ1m z + Bm iξm δ1 e−ϕ1m (h−z ) − Cm iξm δ2m e−ϕ2m z + Dm iξm δ2 e−ϕ2m (h−z )

⎧ ⎫ ⎪u¯ x(1) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪u¯ z(1) ⎪ ⎪ ⎪ ⎪

⎪ ⎨. ⎪ ⎬





(39)

⎧ (1 ) ⎫ P¯z ⎪ ⎪ ⎪ ⎪ ⎪ (1 ) ⎪ ⎪ ⎪ P¯xz ⎪ ⎪ ⎪ ⎪ ⎪ ⎨. ⎪ ⎬ .

(40)

⎪ ⎪ ⎪ . ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ¯z(n ) ⎪ ⎪ ⎪ P ⎪ ⎩ ¯ (n ) ⎪ ⎭ Pxz

The surface stress function of multi-layered mediums, i.e. boundary tractions, can be transformed as Eq. (41) by the Fourier series expansion, where ξ m = π m/L, and Fm is the space dispersion coefficient. The space dispersion coefficient can be presented as Eq. (42).

Pz (x ) =

Fm =

∞ !

Fm eiξm x

(41)

m=−∞

qb/Bm = 0 q/π m sin (π mb/B )m = 0

(42)

In addition, the basis of Inverse Fast Fourier Transform (IFFT) was used to achieve the transformation of the frequency domain to time domain, where the time dispersion coefficient F(t) as shown in equation below.

F (t ) =

N !

Fn eiωn x

(43)

n=1

Therefore, the displacement solutions in the physical domain of the global stiffness matrix and the boundary conditions can be obtained by the double superposition of wavenumbers and frequencies, as shown in equation below.

"

u(x, z, t ) =

N 2M ! ! n=1

#

−1 K¯Global

(ξm , z )Fm e

i ξm x

Fn eiωn t

(44)

m=1

where M and N represent the total number of superposition required for wavenumber and frequency. 4. Numerical results and discussion Based on the analytical solution for the two-dimensional anisotropic, i.e. transversely isotropic, multi-layered medium subjected to harmonic moving load obtained above. Several examples were put forward to verify the accuracy of the proposed analytical solution, and to analyze the effect of interlayer conditions and transversely isotropic properties on the multi-layered medium at different moving speed of load and natural frequency of load. Three subsections were merged into this section. The first subsection focuses on the implementation and validation of proposed solution, while the impact of the transversely isotropic properties, interlayer conditions, and top layer thickness were studied in detail in Sections 4.2–4.4. It is noteworthy that the civil engineering background is used in Sections 4.2–4.4, i.e. simulating typical three-layered flexible pavement, and the mechanical behaviors of multi-layered medium with transversely isotropic-elastic layer-1, isotropic-elastic layer-2 and isotropic-elastic layer-3 were calculated with harmonic moving load. Some studies reported that the influence of transversely isotropic properties of the asphalt mixture layer, i.e. layer-1, is greater than that of the subbase layer and base layer, i.e. layer-2 and layer-3 [3,30]. In this study, the transversely isotropic behavior of the asphalt mixture, i.e. layer-1, and the adopted suitable parameters of harmonic moving load and structure model were used to study the influence of the transversely isotropic properties, interlayer conditions, and top layer thickness on the two-dimensional multi-layered medium. The schematic diagram of a two-dimensional three-layered medium subjected to harmonic moving load is shown in Fig. 4, where the load density q is 400 kN/m2 , the length of load 2l is 0.2 m,

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

31

Fig. 4. Schematic diagram of three-layered medium. Table 1 Material properties of three-layered medium. Layers

Ev (MPa)

Evh (MPa)

μv

μvh

ƞ (%)

ρ (kg/m3 )

h (m)

Layer-1 Layer-2 Layer-3

1400 400 50

n · Ev 1 400 50

0.25 0.25 0.35

0.25 0.25 0.35

5.0 5.0 5.0

2300 2100 1800

0.2 0.3 Infinite

Fig. 5. Comparison with the example in Ref. [24].

and the load cycle 2L is 200 m. In addition, the transversely isotropic property of the layer-1, i.e. asphalt mixture layer, was considered in Sections 4.2–4.4, where the detailed parameters of multi-layered medium were presented in Table 1. After a number of calculations, we found that N = 2048 and M = 2400 can meet the calculation accuracy requirements in the two-dimensional multi-layered medium subjected harmonic moving load. 4.1. Implementation and validation of proposed solution The validation of the proposed analytical solution was performed by comparison with Refs. [24] and [37]. The following two examples in Refs. [24] and [37] were calculated to verify the accuracy of the proposed analytical solution. Ai and Ren computed the standard deflection, i.e. standard vertical displacement, on the surface of layer-1 subjected to the linear moving load [24]. The standard deflection can be represented as u+ z = Gv uz /qa. Fig. 5 presents the comparison of the calculated standard deflection horizontal distribution curve obtained from the proposed analytical solution and all four cases standard vertical displacement from Ref. [24]. It can be seen that the curves obtained from the proposed analytical solution are a good match with the ones stemming from Ref. [24]. In addition, Zuo et al. reported their finding on the vertical displacement change of Gibson foundation under linear moving load [37]. The two-dimensional transversely isotropic multi-layered model of this study was degenerated into isotropic model, i.e. transversely isotropic coefficient n = 1.0. The proposed analytical solution was simulated using the same parameters by Zuo et al., the model parameters were set as follow: E = 290 MPa, μ = 0.45, ρ = 1650 kg/m3 , model depth h = 5.0 m, load length 2b = 4.0 m, and the load density q = 400 kPa. The comparison of the standard vertical displacement varies with the depth obtained from proposed analytical solution and

32

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

Fig. 6. Comparison with the example in Ref. [37].

Ref. [37]. Fig. 6 shows that the curves obtained from proposed analytical solution are also a good match with the ones stemming from Ref. [37]. Therefore, the analytical solution described in Sections 2 and 3 is validated and will be applied in Sections 4.2–4.4. 4.2. Impact of transversely isotropic property on the multi-layered medium Standard vertical displacement of the surface of layer-1 is a representative of the mechanical response index of the multilayered medium. The standard vertical displacement of the surface of layer-1 was selected as a major mechanical index to analyze the dynamic response for two-dimensional anisotropic multi-layered mediums subjected to harmonic moving load in this study. In order to evaluate the effect of transversely isotropic property or transversely isotropic coefficient (n = Eh /Ev 1 ) on the standard vertical displacement of the surface of layer-1 for multi-layered medium at different moving speed and natural frequency of load. The transversely isotropic coefficient of layer-1 selected were 1/4, 1/2, 3/4, and 1. The simulation results were separated in two subsections: (i) at different moving speeds of the load, and (ii) at different natural frequencies of the load. It is noteworthy that the moving speeds were selected as 0 m/s, 20 m/s, and 40 m/s, and the natural frequencies selected were 8 Hz, 32 Hz, and 64 Hz in this study. The simulation results are presented in Figs. 7 and 8. 4.2.1. At different moving speed of loads Fig. 7 shows the standard vertical displacement horizontal distribution curves of the surface of layer-1 at different moving speed and 8 Hz natural frequency of load. It can be seen that the standard vertical displacement increases with the decrease of the transversely isotropic coefficient n, and the influence degree of transverse isotropy on the standard vertical displacement of the surface of layer-1 decreases from the load center x = 0 to side. There is a critical point x = 3.4l when the moving speed is increased to 20 m/s, as shown in Fig. 7(a), where the standard vertical displacement on the left side of the critical point increases with the decrease of the transversely isotropic coefficient, while the standard vertical displacement on the right side of the critical point decreases with the decrease of transversely isotropic coefficient. Furthermore, when the moving speed increases to 40 m/s, the critical point moving forward to x = 2.9l. The above-mentioned phenomenon proves that the critical point will appear when the load changes from static to moving, and the bigger the velocity, the closer the critical point is to the load center. In addition, taking comparison with the isotropic condition (n = 1), the peak standard vertical displacement of anisotropic condition (e.g. n = 1/4) was increased by 12.9%, 11.6%, and 10.2% when the load velocities were 0 m/s, 20 m/s, and 40 m/s, respectively. It can be concluded that the transversely isotropic properties of the layer-1 have a significant impact on the peak value of standard vertical displacement of the surface of layer-1, and these effects increase with the decrease of load velocity. 4.2.2. At different natural frequency of loads Fig. 8 presents the standard vertical displacement horizontal distribution curves of the surface of layer-1 at different natural frequencies and moving speed 20 m/s of loads. There is also a critical point x = 3.4l when the natural frequencies of the load was 8 Hz, as shown in Fig. 8(a), the standard vertical displacement on the left side of the critical increases with the decrease of the transversely isotropic coefficient n, while the trend of standard vertical displacement on the right side of the critical point is opposite to the left side. With the increase of the load frequency to 32 Hz and 64 Hz, the decrease of the transversely isotropic coefficient causes the inclination of the vertical displacement horizontal distribution curve. There is no critical point, which shows a no trend reversal point. By comparing with the isotropic condition (n = 1), the peak value of standard vertical displacement of anisotropic condition (e.g. n = 1/4) was increased by 11.6%, 8.1%, and 6.9% when the natural frequencies of load were 8 Hz, 32 Hz, and 64 Hz, respectively. It can be seen that the lower the natural frequency of

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

33

Fig. 7. Effect of transversely isotropic property on the standard vertical displacement of the surface of la yer-1 for multi-layered medium at different moving speed of load.

load, the greater the influence of the transversely isotropic coefficient n on the peak vertical displacement of the surface of layer-1. 4.3. Impact of interlayer condition on the multi-layered medium Fig. 9(a) shows the standard vertical displacement horizontal distribution curves of the surface of layer-1 for the bonded condition between the layer-1 and layer-2. Fig. 9(b) depicts the standard vertical displacement horizontal distribution curves of the surface of layer-1 for the sliding condition between the layer-1 and layer-2. It is easy to see that the curves for the bonded condition are almost the same as the curves for the sliding condition, and the critical point also occurs at about x = 3.4l for both sliding and bonded conditions. Taking comparison with the isotropic condition (n = 1), the peak value of standard vertical displacement of anisotropic condition (e.g. n = 1/4) were increased by 11.6% and 15.6% for the bonded condition and sliding condition, respectively. It also can be seen that the effect of the transversely isotropic property of layer-1 on the standard vertical displacement of the surface of layer-1 was greater when the interlayer condition between layer-1 and layer-2 was sliding, compared with the bonded condition. 4.4. Impact of the top layer thickness on the multi-layered medium The effect of the top layer thickness, i.e. layer-1, on the maximum standard vertical displacement of the surface of layer-1 for multi-layered medium under the moving load (c = 20 m/s and ω = 8 Hz) is presented in Fig. 10. Five different thicknesses of the top layer were employed in this numerical discussion, i.e. h = 0.16 m, 0.18 m, 0.20 m, 0.22 m, and 0.24 m. Fig. 10 shows that the smaller the thickness of the top layer, the larger the maximum standard vertical displacement at the same transversely isotropic coefficient (n) condition. Taking comparison with h = 0.24 m, the maximum standard vertical displacement

34

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

Fig. 8. Effect of transversely isotropic property on the standard vertical displacement of the surface of layer-1 for multi-layered medium at different natural frequencies of load.

of the thinner condition (e.g. h = 0.16 m) was increased by 7.45%, 7.23%, 7.08%, and 6.05% for n = 1/4, 1/2, 3/4, and 1, respectively. The result of this comparison shows that the influence of thickness on the maximum standard vertical displacement was decreased with the increase of the transversely isotropic coefficient of layer-1. In addition, by comparing it with the isotropic condition (n = 1), the maximum standard vertical displacement of the anisotropic condition (e.g. n = 1/4) were increased by12.34%, 11.80%, 11.60%, 11.41%, and 10.88% for h = 0.16 m, 0.18 m, 0.20 m, 0.22 m, and 0.24 m, respectively. It is easy to see that the maximum standard vertical displacement was more sensitive to changes in the transversely isotropic coefficient when the top layer thickness is smaller. 5. Summary and conclusions Based on the three basic equations of anisotropy, i.e. transverse isotropy, this study combined the Buchwald potential function, moving coordinate transformation, and Fourier transform theory to develop the displacement and stress solution in the transversely isotropic spectral elements. Meanwhile, according to the superposition theory of waves, the element stiffness matrix of transversely isotropic two-node spectral element and one-node spectral element was deduced. By combining the principle of two kinds of spectral elements (that is, the layers are sliding and bonded), the matrices were assembled into a global stiffness matrix of the multi-layered medium. The application of Fourier series expansion and IFFT algorithm to achieve the dynamic response from the wavenumber-frequency domain to real domain transformation. In addition, the consistency and accuracy of the transversely isotropic algorithm in this study were verified by comparison with the literature. At the same time, the impact of transversely isotropic properties, different natural frequencies and moving speed of the load, interlayer conditions, and top layer thickness on the standard vertical displacement of the surface of the layer-1 were analyzed. The conclusion is as follows:

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

Fig. 9. Effect of interlayer condition on the standard vertical displacement of the surface of layer-1 for multi-layered medium.

Fig. 10. . Effect of the top layer thickness on the standard vertical displacement of the surface of layer-1 for multi-layered medium.

35

36

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

(i) The influence of transversely isotropic properties on the standard vertical displacement of the upper surface of the multi-layered medium change with the verity of the moving speed of load or load velocity. At the moving speed of 0 m/s for the load, the standard vertical displacement increases with the decrease of the transversely isotropic coefficient, and the critical point will appear when the static load changes to moving load at c = 0 m/s. The standard vertical displacement on the left of the critical point increases with the decreases of transversely isotropic coefficient, while the value of the standard vertical displacement on the right side of the critical point decreases with the decrease of the transversely isotropic coefficient. The greater the velocity of load, the closer the critical point to the load center. In addition, the transversely isotropic property of the layer-1 has a greater impact on the standard vertical displacement peak value, and that impact increased when the load velocity decrease. (ii) The influence of transversely isotropic properties on the standard vertical displacement of the upper surface of the multi-layered medium also changes with the natural frequency of the load. There is no critical point, i.e. no trend reversal point, with the increase of the load frequency from 8 Hz to 32 Hz and 64 Hz. Meanwhile, the lower the natural frequency of load, the greater the influence of the transversely isotropic coefficient n on the peak vertical displacement of the surface of layer-1. (iii) The effect of the transversely isotropic property of layer-1 on the standard vertical displacement of the surface of layer-1 was greater when the interlayer condition was sliding or/and the top layer thickness was smaller. Acknowledgments This work is financially supported by the National Natural Science Foundation of China (NSFC Nos. 51278188, 50808077, and 51778224) and Young Teacher Growth of Hunan University. The first author also acknowledges the financial support from the China Scholarship Council (CSC) under No. 201606130 0 03. The authors are sincerely grateful for their financial support. In addition, the manuscript has received the written quality improvement assistance from Michigan Tech Multiliteracies Center during the revisions. The views and findings of this study represent those of the authors and may not reflect those of NSFC, Hunan University, and CSC. References [1] S.M. Tafreshi, O. Khalaj, A. Dawson, Repeated loading of soil containing granulated rubber and multiple geocell layers, Geotext. Geomembr. 42 (1) (2014) 25–38. [2] Y.-X. Zhan, et al., Dynamic analysis of pavement and multi-layered transversely isotropic saturated ground system subjected to a rectangular moving load, in: Proceedings of the 2015 International Conference on Mechanics and Mechanical Engineering (MME2015), 2016. [3] L. You, et al., Spectral element solution for transversely isotropic elastic multi-layered structures subjected to axisymmetric loading, Comput. Geotech. 72 (2016) 67–73. [4] H.D.J.N. Aldridge, R.M. Brigham, Load carrying and maneuverability in an insectivorous bat: a test of the 5% “rule” of radio-telemetry, J. Mammal. 69 (2) (1988) 379–382. [5] A. Khanna, A. Kotousov, Stress analysis of a crack in a fiber-reinforced layered composite, Compos. Struct. 118 (2014) 139–148. [6] Z.Y. Ai, W.J. Liu, B.K. Shi, Quasi-static interaction between non-uniform beams and anisotropic permeable saturated multilayered soils with elastic superstrata, Appl. Math. Model. 53 (2018) 400–418. [7] J. Cole, Stresses produced in a half-plane by moving loads, J. Appl. Mech. 25 (1958) 433–436. [8] F. De Barros, J. Luco, Stresses and displacements in a layered half-space for a moving line load, Appl. Math. Comput. 67 (1–3) (1995) 103–134. [9] N.D. Beskou, D.D. Theodorakopoulos, Dynamic effects of moving loads on road pavements: a review, Soil Dyn. Earthq. Eng. 31 (4) (2011) 547–567. [10] N. Liu, et al., Influence of interface conditions on the response of transversely isotropic multi-layered medium by impact load, J. Mech. Behav. Biomed. Mater. 77 (2018) 485–493. [11] L. You, Z. You, K. Yan, Effect of anisotropic characteristics on the mechanical behavior of asphalt concrete overlay, Front. Struct. Civil Eng. (2018) 1–13. [12] L. Wang, et al., Anisotropic properties of asphalt concrete: characterization and implications for pavement design and analysis, J. Mater. Civ. Eng. 17 (5) (2005) 535–543. [13] Y. Feng, et al., Measurements of mechanical anisotropy in brain tissue and implications for transversely isotropic material models of white matter, J. Mech. Behav. Biomed. Mater. 23 (2013) 117–132. [14] F.J. Rodriguez-Fortuno, A.V. Zayats, Repulsion of polarised particles from anisotropic materials with a near-zero permittivity component, Light: Sci. Appl. 5 (1) (2016) e16022. [15] L. You, et al., Impact of interlayer on the anisotropic multi-layered medium overlaying viscoelastic layer under axisymmetric loading, Appl. Math. Model. 61 (2018) 726–743. [16] J.S. Stenzler, N. Goulbourne, Effect of polyacrylate interlayer microstructure on the impact response of multi-layered polymers, Time Dependent Constitutive Behavior and Fracture/Failure Processes, 3, Springer, 2011, pp. 241–258. [17] O. Chupin, et al., Influence of sliding interfaces on the response of a layered viscoelastic medium under a moving load, Int. J. Solids Struct. 47 (25–26) (2010) 3435–3446. [18] C.c. Madshus, A. Kaynia, High-speed railway lines on soft ground: dynamic behaviour at critical train speed, J. Sound Vib. 231 (3) (20 0 0) 689–701. [19] K. Yan, et al., Spectral element method for dynamic response of multilayered half medium subjected to harmonic moving load, Int. J. Geomech. 18 (12) (2018) 04018161. [20] K. Chahour, G. Lefeuve-Mesgouez, A. Mesgouez, Spectral analysis of a railway track in contact with a multilayered poroviscoelastic soil subjected to a harmonic moving load, Soil Dyn. Earthq. Eng. 64 (2014) 24–37. [21] R. Rajapakse, Y. Wang, Elastodynamic Green’s functions of orthotropic half plane, J. Eng. Mech. 117 (3) (1991) 588–604. [22] R. Stoneley, The seismological implications of aeolotropy in continental structure, Geophys. Suppl. Mon. Not. R. Astron. Soc. 5 (8) (1949) 343–353. [23] J. Synge, Elastic waves in anisotropic media, Stud. Appl. Math. 35 (1–4) (1956) 323–334. [24] Z.Y. Ai, G.P. Ren, Dynamic analysis of a transversely isotropic multilayered half-plane subjected to a moving load, Soil Dyn. Earthq. Eng. 83 (2016) 162–166. [25] Z.Y. Ai, Z.X. Li, N.R. Cang, Analytical layer-element solution to axisymmetric dynamic response of transversely isotropic multilayered half-space, Soil Dyn. Earthq. Eng. 60 (2014) 22–30. [26] S. Rizzi, J. Doyle, A spectral element approach to wave motion in layered solids, J. Vib. Acoust. 114 (4) (1992) 569–577. [27] M. Dehghan, M. Abbaszadeh, A. Mohebbi, Legendre spectral element method for solving time fractional modified anomalous sub-diffusion equation, Appl. Math. Model. 40 (5–6) (2016) 3635–3654.

L. You et al. / Applied Mathematical Modelling 67 (2019) 22–37

37

[28] D. Komatitsch, J.-P. Vilotte, The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. Seismol. Soc. Am. 88 (2) (1998) 368–392. [29] K. Yan, H. Xu, L. You, Analytical layer-element approach for wave propagation of transversely isotropic pavement, Int. J. Pavement Eng. 17 (3) (2016) 275–282. [30] L. You, et al., Spectral element method for dynamic response of transversely isotropic asphalt pavement under impact load, Road Mater. Pavement Des. 19 (1) (2018) 223–238. [31] Y. He, et al., An efficient finite element method for computing modal damping of laminated composites: theory and experiment, Compos. Struct. 184 (2018) 728–741. [32] L. Sultanov, R. Davydov, Mathematical modeling of large elastic-plastic deformations, Appl. Math. Sci. 8 (57–60) (2014) 2991–2996. [33] H.S. Lee, Viscowave – a new solution for viscoelastic wave propagation of layered structures subjected to an impact load, Int. J. Pavement Eng. 15 (6) (2014) 542–557. [34] Gharti, H.N. and J. Tromp, A Spectral-Infinite-Element Solution of Poisson’s Equation: An Application to Self-gravity. arXiv preprint arXiv:1706.00855, 2017. [35] Chan, J. and J.A. Evans, Multi-patch Discontinuous Galerkin Spline Finite Element Methods for Time-Domain Wave Propagation. arXiv preprint arXiv:1708. 02972, 2017. [36] R. Al-Khoury, et al., Spectral element technique for efficient parameter identification of layered media. I. Forward calculation, Int. J. Solids Struct. 38 (9) (2001) 1605–1623. [37] Y.-h. Zuo, C.-j Xu, Y.-q. Cai, Dynamic analysis of Gibson soil medium on bedrock under moving load, J. Vib. Eng. 3 (2005) 017.