1

Download as pdf
Download as pdf
You are on page 1of 17
eran oul of Sls nd Suuctuses 160-18 (201) 45-61 Contents lists available at ScienceDirect * sie rn International Journal of Solids and Structures journal homepage: wiwwelseviensonvlocatesolst ELSEVIER Dynamic analysis of layered systems under a moving harmonic rectangular load based on the spectral element method Zhaojie Sun’, Cor Kasbergen®, Athanasios Skarpas"*, Kumar Anupam, Karel N. van Dalen*, Sandra MJ.G. Erkens* *Deproent of Eerie Sucre Fc of Cn gee and Geoscences Oe *Dapertmen of Cl nara and Evite Engng Cole of Engine Ki Ye Umveity fTecnsogy Stevie 1.2828 De ARTICLE INFO ABSTRACT Ire sone ecepteé 24 June 2019 frsle online 25 june 28 it order to design high-performance roadways robust Tool which can Compute the structural response ‘used by moving vehicles is necessary. Therefore, this paper proposes a special element method-based ‘model to accurately and etfectively predict the 3D dynatnie response of layered systems under a moving Toad. A layer special element and 2 semi-infinite special element ate developed to respectively mode 4 layer and a half-space, and the combinations ofthese two elements can simulate layered systems, The fizection, and this property is more dominant at higher velocities. In addition, the amplitudes of these vertical displacements are smaller i he loading frequency is higher or the loss factor is bigger. Finally the loading are and Poisson's ratio only have effect om the displacement amplitudes of point inthe ‘doze vicinity ofthe leading area. The proposed mode! is beneficial to the development of engineering ‘methods for pavement design and i a promising parameter back-aleulation engine for pavemest quai valuation. {© 2019 Elsevier Lid Allright reserved. 1. Introduction steady-state dynamic response of a railway track caused by mow Ing trains, through which an analytical expression of the steady= Roadways are important infrastructures and should be well designed, In order to ensure the performance, 2 clear understand ing ofthe response of roadways caused by moving vehicles is nec- essary. Theoretically, this problem can be regaréed as the dynamic analysis of semi-infinite or layered media caused by a moving load ‘which is generally solved by using either analytical or numerical methods Analytical methods generally give exact solutions to dynamic problems, and these methods are usually efficient. For exemple Eason (1965) investigated the stresses in a semi-infinite elas tic solid caused by moving surface forces with different load- ing conditions by using integral transforms. Vosirouknov and ‘Metrkine (2003) proposed a theoretical model to analyse the mat drs shape sunocucelt Sa c-73/8 208 hens ue A gis seed. state deflection of the rails was obtained, However, the analytical solutions are generally only valid for specific structural and load- ing configurations, and these solutions are usually difficult to culate because they often contain complicated integrals with sin- gular point ‘Numerical methods, such as the finite element method (FEM) and the boundary element method (BEM), ate powerful tools for the dynamic analysis of solid media with different structural, combinations and loading conditions. For instance, Zaghlosl ad White (1993) developed a three-dimensional dynamic finite el- fement program to analyse the behaviour of flexible pavements caused by loads moving at different velocities. “ndersen and Nielsen (2003) conducted boundary element analysis of the steady-state response of an elastic half-space caused by a surface ‘moving load. However, numerical methods are usually time and resource intensive, and numerical distortions may occur in some The limitations of analytical and numerical methods may hin- der their application in engineering, especially for the dynamic analysis of layered systems. Hence, a semi-analytical method called the spectral element method (SEM) (Doyle, 1997: Al-khowry etal 2002; Lee, 2009) is used in this paper to analyse the 3D dynamic response of layered systems caused by a moving load, The SEM is promising for efficient dynamic analysis because it has the ad- ‘vantages of both spectral analysis and finite element method. In the SEM, one element is sufficient to represent a whole layer be- cause of the exact description of mass distribution, which reduces the size of the system of dynamic equations and further increases the computational efficiency. Moreover, this method discretises the continuous integrals into series summations, which is more con- venient for numerical calculation. The SEM has been successfully used for analysing the 2D dynamic response of layered systems For example, You et al (2018) investigated the dynamic response of transversely isotropic pavement structures under an axisymmet~ ric impact loadin cylindrical coordinate system based on the SEM, Yan et al. 2018) applied the SEM to predict the dynamic response of 2D layered systems subject to a moving harmonic strip load However, the SEM has rarely been applied for the 3D dynamic analysis of layered systems under a moving hatmonie rectangular load, whichis the main focus of this study. “This paper includes the detailed mathematical formulation of a 3D dynamic model for layered systems under a moving harmonic rectangular load based on the SEM. The accuracy ofthis model has been verified both numerically and experimentally. The proposed ‘model can be used to analyse the 3D dynamic response of pave- iment structures caused by a moving harmonic rectangular load, ‘which contributes to the development of engineering methods for pavement design. Furthermore, this model could be combined with proper optimisation algorithms to back-calculae the parameters of Pavement structures by analysing the response, whichis useful for pavement quality evaluation 2. Model formulation In this section, the detailed formulation of a model which can predict the 3D dynamic response of elastic layered systems sub- jected to a uniformly moving, harmonically varying. evenly dis- tributed, rectangular surface load is presented. With considering the loading conditions caused by moving vehicles and structural parameters of pavement systems, this model can be used as a tool for structural design to ensure the durability. 211 Introduction of moving coordinate system As shown in Fig. in order to deal with the moving load prob- lem, itis convenient to introduce both a stationary Cartesian com ordinate system (OXYZ) and a moving Cartesian coordinate system. (xyz) Qones et a. 1988: Lefeuve-Mesgouez et al, 2000; Metrkine, 2004), The stationary coordinate system does not move and its ori fin is located atthe centre ofthe loading area when time is zero, ‘The moving coordinate system follows the load and its origin is located at the centre of the moving loading arca. The velocity of the moving load is assumed to he constant and is described by a vector €= [ce ¢ &F. The stationary coordinate vector is notated 3s X=[X YZ, and the moving coordinate vector is notated as x=[x y 2]! The relationship between these two coordinate vec tars can be expressed as follows; x=X-e a in which tis time. These two coordinate systems are coincident When t= ( Nocergen ad A Serpe Iterations ra of Seon Stare 140-89 (2019 45-67 z ‘Adéitionally, the partial derivatives in the two coordinate sys- tems have the following relationships for nonnegative integer 1 pe oe @ 4) | ® OR t+SH HOH where ¢ 2.2. Wave motion in a semi-infinite medium under a surface moving foad ‘As shown in Fig. 2, a homogeneous, isotropic, and linearelastic semi-infinite medium is subjected to a surface load which moves along the X-axis with a constant speed c. The corresponding wave ‘motion in this medium is considered fist. In the stationary coor dinate system (OX¥Z), the equations of motion for the medium can be expressed by Navier’s equation inthe absence of body forces: ay G+ Vo Us WVEU =p: @ inwhich Vo=[& Am det et is the Laplacian operator, UOK.0) = [Uy Uy Uz]" is the dis- placement vector, pis the mass density, 2 andy are Lamé con- ants defined by Young's modus F and Poissons rato». ‘an elegant approach to solve the Naviers equation is using the Helmhole decomposition, which exoreses a displacement field in the following form = Vo = Vox 6) Where @(X.1) is a scalar potential related to the P-wave, and HK.0)=[y Vy Yel" is a vector potential related to the S-wave. It can be seen that the three components of the dis placement vector are related to four other functions, the sealar Potential and the three components of the vector potential, Which indicates that an additional constraint condition is needed (Achenbach, 1999). The additional constraint condition can have different forms (Vostroukhov and Metrikine, 2003; Hung and Yang, 2001), but the solution is uniquely derermined by the govern ing equations and boundary conditions by virtue of the Unique ress theorem. In this paper, the Gauge condition Vo W(X. 1) = 0 is taken as the additional constraint condition, SI is the Del operator, 2 su, © Ksberzen and A. Sharpes ea inerationl Journ obs and Srctres 180-81 (2019) 45-61 * 2 (x= at.0.) 1 z The velocity vector of the load is ¢=[e 0 O}", which means movement with constant velocity ¢ along the X-axis. According to the relationship between the two coordinate systems, fi. (4) has the following form in the moving coordinate system - o(& S)u © 4:T is the Del operator in the moving co- Be the moving coordinate system, WOK.) = [tie uy uel! is the dise placement vector in the moving coordinate system. In the moving coordinate system, Eq. (5) has the following Ba VOrVKd a Ge wVT my uv in which Vala 3 ordinate system, VE is the Laplacian operator in where $6) and YOX.t) =[vx y. Yel! are the scalar potential and the vector potential in the moving coordinate system, respec- tively. The Gauge condition in the moving coordinate system reads vow By substituting fn, (7) into £. (6), the following uncoupled wave equations in the moving coordinate system ate obtained (For more details see Appendix A) ) in which cp = = Zpil/p isthe velocity of the P-wave, and ITIP isthe velocity of the S-wave. Im order to solve Eqs. (8) and (°), -m pair with respect to time is used! the following Fourier trans- 4 LP ae nent , Ge.0)= ge fae 0) 40 [ aee.0peaw mm) where iis the imaginary unit satisfying =—1, o is angular fre- quency, q(x.¢) is an arbitrary quantity in time domain, and 4tx.@) Js the corresponding quantity in frequency domain. After applying, the above forward Fourier transform, Eos. (7)-(9) become a) a) a) n which the “hat” means that these quantities are expressed in the frequency domain Im the Cartesian coordinate system, the solutions of, qs. (13) and (14) can be retrieved in exponential forms. Ac- cording to the phase matching principle (hao et al, 2015), different waves should have the same phase at the boundary (eg. the surface 2=0). Consequently, the P-wave and S-wave have the same wavenumbers not only in x-direction, but also in y-direction, Therefore the general expressions of G(x a) and P(x.) are: $4.0) = Aeeretire tbe as) daoy=[e ¢ ple Wwhete A.B, C D are unknown coefficients to be determined by the boundary conditions, ke isthe wavenumber in the x-direction, ky is the wavenumber in the y-rection, and kre and ks are respectively the wavenumbers in the z-direction for the Pawave and S-wave Note that the signs of kyr and ks, should be chosen carefully to ensure that the waves propagate andjor attenuate in the positive Zedirection, After substituting Fgs. (15) and (15) into Zs. (13) and (1), the expressions of kp, and ke, can be obtained: shxgetiy gritos as) = Crt gig a7 wg - ey gig (as) By substituting Fa. (12) into the expressions of the constitu- tive equations in frequency domain. the following relationships be- tween the stresses and the potentials in frequency domain ae ob- tained é » eR w#)\s #6 he Py saan -a( Ber Br Boone SB-8-28) (a9) 4“ 2S Rosberg a A, Sapa ea Internation! oul of So and Sractares 10-181 2009) 45-61 (20) We, (21) exon alo ad ay, (a) 5) svamoi[pbes 22-28 -(4-2)s] my we Waa auc 0) =u] ayes (24) in wiih Yas Vp and ae the tree components of $6 ta adn by combining the expressions of ey ad ey wth qs (17) and (28), the telationship between Lamé constants can be obtained +8 +26, - RB, Eee, 25) 23, Spectral element formulation In the SEM, the response of an element is determined by its total wave field, which is the superposition of wave fields origi- nating ftom different boundaries (AI-Khoury et al, 2001). In this method, the number of elements needed for a simulation is the same as the number of layers because one element is sufficient to simulate a whole layer, which makes it efficent for dynamic anal ysis of layered systems, In this section, a layer spectral element 1d a semi-infinite spectral element are formulated to simulate layer and a half-space, respectively. The combinations of these two spectral elements are capable of modelling different layered systems 231, Layer spectral element ‘As shown in Fig. 3(a). the layer spectral element consists of two parallel horizontal ectangular surfaces, which constrain the waves to propagate within the element, The element vertically covers the whole simulated layer, and it horizontally extends to a certain is tance after which the response caused by waves is negligible. In addition, the spectral element of a layer with thickness h is physi cally defined by two nodes located at (0,0, 0} and (0, 0,f), each of Which has three degrees of freedom. In the layer spectral element, the total potentials can be expressed as follows (A-Khoury et 2001; van Dalen et al, 2015) arf agertne « Apetnee] 6) Byentst 4 yee”) Feo) —eteeeo] Cheater y Creaeen on - Dyertoe Dye where A), Ar, By. By, Cy, Cy Dy, and Dp are the unknovin coef cients to be determined by boundary conditions, The first terms are the potentials of the wave fields originated from the top st face, while the second terms are the potentials of the wave fields reflected from the bottom surface, ‘The forward Fourier transform is applied to the Gauge condition in the moving coordinate system to obtain its spectral form. which is expressed as follows: v le (Ka) =0 (28) ‘Alter substituting Eg, (27) into Eq, (28), the following relationships ae obtained! kab 4G 29) xB +hyG Ds Ter (30) 2 su, © Ksberzen and A, Sarpos ea nection Journ os and rcs 180-81 (2019) 45-61 “” qs, (29) and (20) ate substituted into £9. (27) fist to decrease the number of unknown coefficients. Then Eqs. (26) and (27) ate substituted into 84. (12) to obtain the expressions of the displacements in frequency domain, which have the following forms: 0) tetera seins) Meet meta) BEB crete] ey 0%, ) = ~ieroreo [ss ene Ayelinecn) « a (eres — Beet) 4 BE Center Ge so)| (32) 0,0. 0) = ee [hg (Aye — Ae) — by (BeBe Bye) 4 (Gere 4. Celie“) (3) The displacements of the top node are notated as. f}. and d}, and the displacements of the bottom node are notated as @, 3, and @, Then, the coordinates of the nodes are substituted into Fqs.(21)-(33) fo obtain the nodal displacements, wiich can be expressed as follows g-fa ay Jn which the superscript “e" means the corresponding quantities are expressed for an element, 5 =[0} d} d) aa aBYt isthe nodal displacement vector of the element. &= [Ay Az By Bz, GF isthe unknown coefficient vector. and fis a frequency and wavenumber dependent matt which has the following expe ih ky ke e inh Ket Ke ohe Tes eye f, -ke he “ ied HEM ae BAR ak a“ Te Tse Tse set Hehe Kyeed ky myer ake By substituting Fqs. (26) and (27) into Eqs. (19)-(24) and considering Fa, (25), the expressions of stresses in frequency domain can be obtained! (Kg 1G ~ 2, +12.) (Aree + Agere») bat w) =e" vel 2c) HBR eg | i) (1 + 248, ~ (Aver + Ape) = prenhete ity 24k 2 ) =u ee BB) (get metoteh) - Be eee .»| 68) . yf (+R) (terns 4 Agee) an w(t Blk (Ave 4 Agen) (BeBe — Byes) F(t) = wee Oe a (a) AEA) tee Genin) . an ty| Bre ee! — Aye) + (HB +1) (Breet + Brea) (3s) FQ 0) = —we Me “fan 4 Cele“) | ” - angry] here(Ave“te® = Aneta) — Ika (BrerHet + Beaten) ; Fut 0) = pero (een Geto) | “ Based on the Cauchy stress principle, for a certain surface, the relationship between the surface traction vector t and the Cauchy stress matrix g can be expressed as follows: ay so 2 Sen Kosbergen a A, Sapa ea Internation! Jounal of So and Sractars 10-181 (2009) 45-61 where mis the outwatd unit normal vector ofthe suifce. The tractions of the top node are denoted as ff), and, and the tractons of the bottom node are denoted as . and. On the bass of Eo (41) the nodal trations have the following relationships with the nodal Cauchy stresses a), (42) (3) in which 4), 3, and 24 ae the Cauchy stresses of the top node, and 42, 83, and 43 are the Cauchy stresses ofthe bottom node “The nodal coordinates are substituted into Fos (25)-(40) to derive the nodal stresses, which are then Incorporated into Eqs (42) and (43) to obtain the expresions of nodal fractions, which are expressed as fa (aa) where =lé ff which fas the following form: matrix fF is the nodal traction vector of the element, fis a frequency and wavenumber Gepende [cher = 2kakseHO® = Dhay— —Dhakyeed B Re Dhykee — —DkyheseHnh Bett Dhaky keke sh -B ee ee MY akckese Boh Dhakre kako Dkk niger Dykre Kher ee ee R DiyksseHeh Dyke, —Dhaker et Dh with =H. 8, =, and =k 8 By combining Eqs" (34) and (44), the relationship between Which can be expressed as he nodal traction vector and the nodal displacement vector is obtained, Be (45) in which Appendix 8, ‘can be regarded as the element stifiness matrix, and the detailed expressions of its components are shown in 232, Semi-infinite spectral element {As shown in Fig. 3(B), the semiinfinite spectral element is composed of a horizontal rectangular surface, and physically defined by 2 node located at (0,0, 0) with three degrees of freedom. Inthe semi-infinite spectral element, the waves originated from the surface travel in the positive z-direction and no reflection occurs, which physically means that the energy is radiated away. Actually, the semi-infnive spectral element can be regarded as a special case of the layer spectral element that only contains the top sutface, which requites the coefficients of Ay, Bz, Cz, ad Dy in Ecs. (25) and (27) to be zero. Accordingly, the displacements of the semi-infinite spectral element in frequency domain can be expressed as fellows: a. RAK, Keck (0, 2) = ie MEM (eye? — hy Bye ype) (4s) At le Substituting the coordinates of the node, the nodal displacements can be expressed as: aky GS a Ts ay gs Wee (49) a “e & The stresses in fequency domain become: é itxe-ity| 2 ya eins — 2b g pings Z(H) ies 50) suc) =-prtel Ura, yee 2g ee BD «9 Gy (00) werteeo| (gap 28, 8) etn EI) gis 2 ge 6) Bex (X,c2) = pe-HotetHor| (KE + 1B — AG, Are“ + Dh seB eH — Dhaest] (52) (lH +R) auto) ~-neietfattyaerte METH Hs ne, HEB ene sn Tr 2 suc Ksberzen ad A. Sapos ea finernationl Journal os and rcs 180-81 (201) 45-61 ” Ge 0) = Were H9 [Dig keeAe RH + (KE — + Bee + Dh Certet] (4) Gn( 0) = — Wee Dhakyaye MY — Dhak Brer H+ (KE — KG REC] (55) ‘ter substituting the nodal coordinates and consider ba, (42), the nodal traction vector is expressed as # Dhak 2k) =I) ra, B)=n| akin Ig, ak [5] 6 B ARAL, 2h kee 2kaks, | LO By combining Es, (49) and (56), the relationship between the nodal traction vector and the nodal displacement vector is obtained! (E+E) Kehse kay re = ka) kal ‘a KES KG + hk: KEG + heekse FRE + Rooks | te ay ‘ (6+ Rhee + Base by a B We Beker + Bis oie a 57) 2 IE Ky + eek keke || ot & “ yk (ek ° BaR+ @ With 12 = 12 43 + kpkse = 24, Boundary conditions [As shown in Fig. 2, the external load is applied on the surface of the layered system in the positive Z-drection. The load is assumed to be a uniformly distributed traction over 2 rectangular area, the amplitude ofthe load varies with time. in the moving coordinate system, the loading area is fixed and it can be expressed as follows Pe(RY.E) = ho YPC (58) where p.(xy:t) is the traction applied in the positive z-direction, ho(x) isthe spatial distribution function of the traction without dimen- sion, p(t) is the loading history function of the traction with dimension of forcelarea. The spatial distribution function h(x) can be expressed as follows: hoGe.y) = He ~ DH ~ YD) (59) Jin which 1.) is the Heaviside function, 2x9 is the length of the loading area in the x-direction, and 2yo is the width of the loading area sm the y-direction According to £9s. (25)-(41), the distribution of tractions on the horizontal surfaces is in the form of e~¥e~%¥, o ha(xy) should be expressed in the same form to match the traction conditions, This can be achieved by using the Fourier sevies representation: 140.3) = OY lg Hote (60) where m isan integer that ranges from -M to M and n isan integer that ranges from —NV to N, where M and N should be large enough to ensure the accuracy of the representation. In addition, kim ~ my and kg ~ 7 /Yo, where 2X ts the length ofthe space window of interest inthe x-direction, and 2¥p is the corresponding wath in the y-direction, The dimensions of the space window should be large enough to cover the infiencing area of the applied Toad, The hana fy ate the Fourie coeicients defined as fllows hen - 2” ne (61) fea = 3 [Heo be (61) 1 pe Sy 62) in = 395 |, Ho betray (62) The moving load considered in this paper is harmonically varying. hence the loading history function pe) can be expressed as follows: BXt) = poe (63) in which py isthe amplitude of the traction, wg = 2xf with wp being the loading angular frequency and fq being the loading frequency. By applying the forward Fourier transform based on Fs, (10), the expression of pt) in the frequency domain is obtained: B(@) = pod (w ~ a) (64) where 8() Is the Dirac delta function. Hence, the expression of the applied surface traction in frequency domain is Psl8.¥.0) = BO) SD hale Horetoy (65) In addition, it can be concluded from Ec, (65) that the expressions of the potentials should be represented as summations over all Kae and kyy fo match the traction conditions 2 2 Sen Kosbergen a A, Sapa ea Internation! oul of So and Sractars 10-181 (2009) 45-61 25, Solution scheme ‘According to the SEM, the combination of several layer spec~ tral elements on top of a semi-infinite spectral element is capa ble of simulating 2 layered system. The numbering and assembling of these elements follow the same procedure as in the traditional FEM, However, because of the wavenumber dependence of the el- fement stfiness matrix in the SEM, the whole assembly process is done for each wavenumber combination. The total number of the odes in che spectral element model ofa layered system is notated as Land the global system of equations for a certain wavenumber Combination can be expressed as: TS) =k") Ge) (665) in which the superscript “mn” indicates that the quantities cor respond 0 a certain wavenumber combination of ken and Kye 72" (@) is the global nodal traction vector with dimensions 31 by 4.8") is the global stiffness matrix with éimensions 3! by 31 and 05" () isthe global nodal displacemes Sions 30 by ‘According to F. (65), the traction of the top node ean be ex- pressed as follows: 50) Dank vector with dimen- (0,0,0) (67) ‘Therefore, the global nodal traction vector for a certain wavenum- ber combination can be expressed as follows: BP) = BoMhmmh aes (68) here eis 231 by 1 unit vector with the third component being 1 ‘According (0 £4. (65), the global nodal displacement vector for a certain wavenumber combination is calculated by: Dye) = BeorhmmbinG™)-e5 (69) in which @™(o), the inverse of Rw), can be regarded as the transfer matik. The nodal displacement vectors for dillerent wavenumber combinations are summed to obtain the total nodal displacement vector caused by the applied oad, such that 0,60) = 10) 2 henlin@™ () -25 (70) where 00) isthe total nodal displacement vector. Thea, the in- verse Fourie transform is used to obtain the nodal displacements in time domain, hich can be expressed a: att) = Boe Fp ‘According to Eqs. (31)-(33), fora certain wavenumber combina tion, the displacement vector of the horizontal plane where a node is located equals the product ofthe nodal displacement vector and. the term e~Yex¥enbn, Therefore, the displacements of points on the nodal horizontal planes can be calculated as UEE™ 6.9.6) = Poti OY henfine Hote S™ (0) 85 (eo) -€5 my (72) in which UP*™ (xy. ¢) isa vector with dimensions 31 by 1, which Contains ihe displacement components of al the horizontal planes where nodes are located if the displacement field in a specific layer is desired, the following steps can be followed. First, one obtains the nodal displacement vector of this layer for certain wavenumber com= bination ffom £4. (68). Secondly, one calculates the correspond- ing coefcent vector via Es. (34). Then one substitutes these co- ficients into gs. (31)-(33) and sums over all the wavenum- ber combinations to compute the total displacement eld in fte- ‘quency domain. Finally, one applies the inverse Fourier transform Via Eo (12) to obtain the total displacement field within this layer in time domain. Corresponding stress and strain fields can also be calculated by using the constitutive equations. The procedure to determine the response fields in a half-space is the same 35 that for a layer. Ie should be highlighted that all the calculated response fields are steady-state Solutions, so they are varying over time with the Same frequency as the applied load. For a certain sesponse fleld (displacement field, stress field, or strain fleld),itcan be expressed, 2s follows: £00) = Foe on)e™ (2) ‘where f(x, ¢) isa certain response field, FAx wg) is the correspond ing time-independent quantity which is normally complex-valued For different loading history functions, a certain component of the response field vector has different forms: Re[ Fi. ea)e™'], ple) = pocos (wat) Iin[Fen)e™"], pO = pasin (oot) in which the subscript “k" represents the considered component of the corresponding vector, Re) denotes the real part of a complex term, and Im) denotes the imaginary part ofa complex term. ‘Assuming the loading history isin cosine form, Ea. (74) ean be revriten a follows according! Au&.1) = RefFie wo] e0s (ost) ~ | F ae)] sin oat) (75) 9, (75) indicates that a response field component equals to the real part (oF imaginary part) of corresponding time-independent ‘quantity at a specific time. In adlton, Eq (75) can also be written fe) { m4 JF. 0)] C05 aot + 8K, 20)] (78) where |F,(&..)| is the amplitude of vibration, and Q(X a) is the corresponding phase angle which satisfies tan[G, (x. 0»)] na It can be concluded from Fo, (75) that, in the moving coordinate system, any response quantity of a point is harmonically varying ‘with the same frequency as the applied load, but diferent points have different amplitudes and phase angles, which consequently forms a periodically varying profile over time. In this paper. all the results are presented in che moving coordinate system; cor responding results in the stationary coordinate system can be ob- tained based on the relationship between the coordinates. Working, in the moving or stationary co system should give equiva lent solutions, because the physical nature of the problem is coor dinate system independent (Louhshalam et al, 2013) ‘though the presented model is formulated for elastic layered systems, it can be combined with ¢ifferent damping models to simulate layered systems with damping. Note that the damping ‘models should be transformed to the moving coordinate system Additionally the presented model can handle different types of surface moving loads by changing the spatial distribution function and the loading history funetion of the applied load. In this paper, 2 hysteretic damping model defined in the frequency-wavenumbet domain related to the moving coordinate system is used to simu late the damping effect in the system by replacing Young's modu- lus E with F[1+insgn(o+ ck) im which n is the loss factor and san() is the signum function. In addition, in view of the practi- cal speeds of vehicles on roadways, all the considered velocities of the loaé are taken smaller than the Rayleigh wave speed in layered systems. 1 Sun. € Kesbergen a A, Sars ea ntenationel oul of So and Srastares 10-181 2018) 454 8 3, Model validation The accuracy ofthe presented model is validated in this section. [A first, this model is implemented in a computer program to com- pte the response of a layered system by executing the following steps: (1) For each wavenumber combination. it calculates the element sliffness matrices and assembles them to the global stiffness (2) Ie applies the boundary conditions and computes the global nodal displacement vector by solving the corresponding slobal system of equations (3) Ie calculates the response field within a certain layer on the basis of the nodal displacements and obtains the total re sponse field by summing all the contributions at different ‘wavenumber combinations. Then, two cases are used to compare the simulated results with corresponding boundary element solutions given by Andersen and Nielsen (2003), These two cases consider the surface deflections of 4 homogeneous half-space and a layered system caused by a mov Ing harmonic rectangular load. The points along the x-axis on the surface are considered in the result comparison, where specifically the corresponding amplitudes and phase angles of displacements Jn zdirection ur(¢0) are analysed. Note that the loading ampl- {ude used in the current paper is 10° times that in the reference literature to make the results comparable with realistic pavement response. Finally, the proposed model is validated by comparing simu lated results with field measurements, A pavement testing facility called LINTRACK (For more details see Appendix C) was used to ‘measure the strains of a pavement structure. The measured maxi- ‘mum longitudinal stains (in the moving éitection of loads) of the pavement structure are used for comparison with corresponding Simulated results. 31. Response of « homogeneous half-space under @ moving harmonic load This case considers the dynamic response of 2 homogeneous half-space caused by a harmonically varying load moving on its surface. The load is uniformly distributed over an area of 3 by 3 m?, and the amplitude is 1/9MPa (instead of 1/9Pa in the erature). The load varies at frequency of 40FHz and moves in the positive direction of the x-axis with velocities of 0, $0, 100, and 150m/s. The structural parameter values of the half-space are shown in Table 1, these parameter values are corresponding to some unsaturated sandy soll with moderate stfines. The amplitudes and phase angles of the displacements in = direction u(x,t) for points along the x-axis on the surface of the half-space are calculated by the presented SEM-based model. In of der to obtain converged solutions, 4096 x 4096 wavenumbers are tused, and this holes for all the results shown in this paper, The simulated results are compared with those given in the reference literature (Andersen and Nielsen, 2003) in Fig. 4. The comparison shows that the results calculated by these two methods ate almost Identical for different moving velocities, which proves the accuracy of the proposed semi-infinite spectral element. In addition, some observations can be made ‘ble (1) When the load does not move, the displacement amplitude curve along the x-axis is symmetric with respect to x=0 and the displacement amplitude is maximura at x=0. (2) When the load moves, the displacement amplitude curve slong the x-axis is asymmetric with respect to x0, The dis placement amplitudes at the points in front of the load de- ‘tease more rapidly than on the other side, and this trend is ‘more obvious if the moving velocity is higher. (3) When the moving velocity is increased, the position of the peak of the displacement amplitude curve along the x-axis shifts to the lef, and the maximum value is slightly higher (4) When the moving velocity is zero, the phase angle curve along the x-axis is symmetric with respect to x=0, How- ever, with inereasing moving velocity, the phase angles of u, at points in front of the loading area change more rapidly, fand consequently the phase angle curve is denser on this side, 32, Response of «layered system under a moving harmonic load ‘This case considers the dynamic response of a layered system. caused by a uniformly distributed harmonic load moving on its surface. The loading area and amplitude are the same as those in the case of the hal-space, while the loading frequency is 2012 and the moving velocities are 0, 25, 50, and 7Smjs in the posi- tive ditection of the x-axis. The layered system is composed of 2 horizontal layer with a certain thickness and a homogeneous hall space. The structural parameter values of the layered system are shown in Table 2. The parameter values of ths layered system cor- respond to two kinds of soll, and the soll in the layer is softer than, that in the half-space. “The displacements in the z-direction u(X,) at points along the -caxis on the surface of the layered system are computed by the presented SEM-based model, and the corresponding amplitudes and phase angles are compared with those given in the reference literature (Andersen and fielsen, 2002) in fig. 5, The comparison, indicates that the result calculated by the diferent methods have good agreement for different moving velocities, which confirms the accuracy of the proposed layer spectral element and its combina- tion with semi-infinite spectral element. Additionally, same obser vations can be made (2) The displacement amplitude curves along the x-axis on the layered system surface have similar changing trends 25 in the case of homogeneous half-space if the moving velocity is increased, However, the curves have some fluctuations for the layered system, which might be attributed to the com- plicated wave field in the layer spectral element, The half- space has higher stiffness than the layer above it, so the con- wibution of the reflected waves is pronounced, (2) The phase angle curves along the x-axis on the layered sys- tem surface are more complicated, but the changing trends ate similar to the case of homogeneous half-space when the ‘moving velocity is increased. Te results shown in this section indicate that the displace- ‘ment amplitude curve decreases moze rapidly in froat of the load ig afea, which is mote obvious at higher velocities. The reason, of this phenomenon is the uneven wave field distribution in the rl parameter values ofthe layered sistem 2 2000 30004594 tafe Py 2.5 € Kosergen and A Sharps eterna! Jura of Sods an Sacre 180-181 (2019) 45-61 as SEM 6 Refine Hn\\ xo) ei Phas angle of, (°) fae Ample of, (10* m) & bi SEM + Reteence i/ LA pees Amplitude of (10%) Phas anle f°) xo) / Ln sta) Amplitde of (10m) Phase angle of, eee aa 2 SEM © Reference z \ sa & Phase nae of) Amplitude of (10* m) 2 6 A xo) xm) Fie 4. comparsen comin 1 pins along the sans onthe hal space surfcecleated Wy dierent metieds at diferest mang velo (2) = Om, (b= 50m. 5 2 sun © Ksberzen and A. Shapes ea internation Journal of Sos end Srctres 180-81 (201) 45-61 ss ecm a B 7 Os bes? Sos lO TaD 5 2 &) 200 SEM 2 Retr rl i 6% / i my ff = / \ 7 sus 1» ' Phase angle of ,(°) Ample of (10° m) os. on f i : See TTT aS © son SES ae iad 190 \ v ie prs i \ \ E sw¥l se 7 is 0 ocean eee eenienetn stm) st Fig 5. Compatson of for pons along te sans onthe layered sistem surface caleuted by erent metheds at ferent moving velit: a} ¢= 01s. (b) = 255. (che nsomis na [€) e752 6 2 Se Kosberge a A, Sapa ea Internation! oul of So and Sractars 10-181 (2019) 45-61 Parareter ales ofthe ested pavement suctre Suberde 175350404 ninte Campion between the lite and measired maxim longi rine aes Tones ‘ian Togo sans (0-7) Vicinity of the loading area caused by the Doppler effect (Lefeuve- Mesgouee et al, 2002). The wavelengths of the waves in front of the loading atea are shorter while the wavelengths of the waves behind the loading area are longer. Hence, the moving load has a smaller influencing area in front ofthe load than behind it. 2.3. Comparison with fed measurements LINTRACK was used to measure the strains of an asphalt pave- iment structure which was designed fora heavily loaded motorway. ‘The first layer is porous asphalt concrete (PAC), the second layer is newly applied stone asphalt concrete (New STAC), the third layer is old STAC, the fourth layer is asphalt granulate cement (AGRAC), and the foundation is a thick and well-compacted sand subgrade ‘The parameter values ofthe tested pavement structure are shown in lable 3. Strain gauges were installed at the bottom of the first layer in the longitudinal direction (direction of movement). Dur ing the measurements the LINTRACK belt moved straight over the built-in strain gauges at a constant speed of 25ms. A constant force was applied on the tire, while the tire pressure was main- tained to be 900KPa ‘The maximum longitudinal strains of the pavement structure calculated by the presented model are compared with those mea- sured by the strain gauges, The results are shown in Table 4, which, indicates a good match between the simulated and measured data, and thus further proves the accuracy of the presented model 4. Response analysis of a pavement structure ‘This section focuses on a specific pavement structure subjected, to a surface moving load, and the parameter sensitivity analysis, and stress analysis are conducted, The reference loading conditions ate described as follows: + Auuniformiy distributed harmonically varying load moves in the positive direction of the x-axis on the surface of a pavement structure; + The moving velocity is ¢ + The loading frequency Is fy =20 Hz + The loading amplitude is pp =550 kPa + The dimensions of loading atea are 2x) =2yp =0.2683 m: + The dimensions of the space window are 2X =2Yp = 400m. 5 ms (80 kmh): ‘The amplitude of the total force applied on the surface is about 39.6 XN, which is comparable to the actual trafic load. The detailed Feference parameter values of a pavement structure are shown in ‘Table 5 Reference parameter values of «pavement suc lt ra = 41, Parameter sensitivity analysis By using single factor analysis, the sensitivity of the displace- ‘ment amplitude curve along the x-axis on the pavement surface to diferent parameters is investigated. The results are shown in Fig 6, In which the response of the reference structural configura tion to the reference loading conditions is shown in solid line. It is assumed that all the layers in the pavement structure have the same Poisson's ratio and loss Factor. AM. Sensitivity o moving velocity The displacement amplitude curves along the x-axis on the pavement surface caused by a load moving at different velocities (c=5, 25, and 45mjs) are shown in Fig. (a). The effect of the ‘moving velocity is similar to that observed in the previous section, However, fora realistic pavement structure the curves are smooth because of the relatively high structural si ness, and the Doppler effect is not as significant as that observed in the layered soll y= tems. In addition, within the range of analyses, the maximum of the displacement amplitude curve is slightly affected by the mov= ing velocity, 412, Sensitivity to loading frequency ‘The displacement amplitude curves along the x-axis on the pavement surface caused by a moving load with different load- ing frequencies (fy = 10, 20, and 30H2) ate shown in Fig. 6(b) The vertical displacement amplitudes of the surface points along the moving direction are smaller if the loading frequency is higher, Which might be the result of the damping mechanism playing a more pronounced roe, 413, Sensitivity to loading area The displacement amplitude curves along the x-axis on the pavement surface caused by a moving load with different loading teas (59 = 0.036, 0.072, and 0.108 m”) but the same amplitude of the total force (38.6kN) are shown in Fig. o(e) It can be seen that the maximum of the curve is higher ifthe loading area is smaller, whieh is caused by the increase of the traction amplitude. How ever, the differences appear only in the close vicinity of the load- ing area, the displacement amplitudes of points outside are almost identical, Therefore, if the applied force is the same, the effect of the loading area is localised in the close vicinity of the loa. 414, Sensitivity to loss factor ‘The surface displacement amplitude curves along the s-axis for pavement structures with different loss factors (n = 04, 02, and. (0.3) under reference loading conditions are shown in Fg. (2). It can be seen that the curve is slightly lower if the loss factor is higher. More energy is dissipated for a system with higher loss fac~ tor, which results in smaller displacements 415, Sensitivity to Poisson's ratio ‘The surface displacement amplitude curves along the x-axis for pavement structures with different Poisson's ratios (v = 025, 2 sun © Kasberzen and ‘Sarpes ea iverson Jura of ats and rctres 180-185 (2019) 45 2 (a) © = Sms ) © aes] i i g? 2? 2 2 © @ @ @ . & i : = = : : : © * = v= 025] : nu j i, 2 a 2 1 0 1 2 0.35, and 0.45) under reference loading conditions are shown in ig 6(e). The maximum of the curve is slightly smaller if the Pois- son's ratio is larger, while the displacement amplitudes of points cutside the loading area are almost unaffected. It should be highlighted that Fig. 6 only shows the amplitudes of u for points along the x-axis on the surface. In reality, all quan~ tities at all points are harmonically varying, 2s shown in Eq. (75) Furthermore, for a pavement structure with the reference struc tual configuration subjected to the reference loading conditions the profiles of us for points along different axes on the pavement surface for t=0 are shown in Fig. 7. The results show that the pro- fle of u, is asymmetric along the x-axis while symmetric along the a} moving ele 8) aang requeny aang axis, which means the Doppler effect appears only in the moving, direction. 42, Stress analysis For a pavement structure with the reference loading and struc- tural configuration, the stresses of points along the x-axis at depth 0.1 m are simulated by the presented model, The results for t= 0 are shown in Fig, 8, which indicates thatthe points under the load Ing area are most critical, For these points, the maximum stress component is oz, Which is followed by oan. zx and dy. In ad dition, the stress components of oy and oy, ate negligibly small ss 2.3 € Kosergen and A Stars entrain Jura of Sods on Sacre 180-181 (2019) 45-61 @* & + xo) vo) Fig 7 oles of uo pelts long diferent aes onthe pases sue fer £0: (2) als and (b) ye. @-* Ong 3° g° sm) xt © ( 5 ow a. sm) He) ig. Ssss of pints slong the sane 2 ep Om fr t= 0: (8) (8) () #8 A) ‘The stresses calculated by the presented model could be used for pavement structural design to ensure its durability. 5. Conclusions and recommendations ‘This paper proposes a SEM-based model which can be used to analyse the dynamic response of layered systems caused by a mov ing load, Based on the discussion shown in this paper, the follow ing conclusions can be drawn: (1) The proposed model is robust for the dynamic analysis of layered systems under a moving load, and this model is a potential tool for pavement structural design. (2) The displacement amplitude curves and phase angle curves along the moving direction are asymmetric when the load ‘moves, and this asymmetry is more dominant if the moving velocity is higher. The reason of this phenomenon is the i= homogeneous wave field distribution caused by the Doppler effect. However, the moving velocity only has slight effect on the maximum of the surface displacement amplitude curve within the scope of analysis, (3) The surface displacement amplitude curve will be lower if the loading frequency is higher or the loss factor is bigger, and the effect ofthe former is more dominant, (4) IC the amplitude of the applied total force is constant, the loading area only has influence on the displacement ampli- tudes of points in the close vicinity ofthe load, (5) The Poisson's ratio has slight effect on the maximum of the displacement amplitude curve, and it almost does not affect the displacement amplitudes of points outside the loading (6) The presented model is 2 promising parameter back- calculation engine for pavement quality evaluation, 2 su, € Ksbrzen ad A, Shapes ea fnernationl Journ Sos and ictus 180-81 (2019) 45-61 sa This paper proposed 2 3D dynamic model for elastic layered systems under a moving load, which is combined with a hysteretic damping model to analyse the response of a pavement structure caused by a moving harmonic load. In order to consider the frequency- dependent viscous effect in pavement structures tis recommended to use more suiteble damping models. Conflict of interest None. Acknowledgements This work is financially supported by the China Scholarship Council (No. 201608230114). The first author is grateful to the discussions with Mingiuan Zhao, who is 2 PhD candidate of the Department of Engineering Structures, Faculty of Civil Engineering and Geosciences, Delft University of Technology. ‘Appendix A In the moving coordinate system, the Navier’s equation has the fllowing form: G+ n)VV w+ nV ‘The Helmholtz decomposition of the displacement field is expressed as: BaVorVxy Substituting Ea. (0 ‘operators give into Eo. ((\1), considering the identities of V.Vp = V2 and V. Vx Y=, and interchanging the order of the [ 2 aay. 2 aay, . Vv) Lav? o-0( 3) $] [pe ¥-0(F- ) 4] (a3) This equation wil esate ifthe tems in the square brakes vanish, hence: ve- Jeno M) vy ae (as) with ep = VETER and eg = ATP W the velociy vector has the frm f= [e © Of shen the Fs (At) and (AS) become 2 1(@ ay CAB) vo 8(-e2)'o-0 as ny 2 (8 an ve-2(% ° an) ‘Appendix B The element stifness matrix ofthe layer spectral element can be expressed a follows: fy fs fy kB he i, i, Be Bee Be, i a i i a) in which mm guy te -apyentn gp asin wo 2 Kesbergn a A, Sapa ea Internation! Jounal of Sod and Sractares 10-18 2009) 45-61 eWay [2ke(1Q +19 + 392,)[e- Der As] — (2418 pak) eh )fe-2e Ha + 1] Fe eumene erlbr 89] — (42 1 — hs) n+ hs) (e-2HH 4 ens) | fog = ke f2[ (6 018) 08 1 20) lem eh] = ta) A hd HM] SB | satpaae[3(2 +18) 12ers — lth] — 2 4B — oaks) (H+ 3 — Zhe — 12.) (eB — es) ie Abre[kf +1BK3 + (KE + AKG) Jel 2B — hE IRE + BNR, + 20 ake + A) [e tee Ae Hat] a { bk HG + Hert eee AK) (+ RBAS + EHR, ~ 2 rae + [eR es 4} Dickey fit 4B 18, + Akko eh) — (AZ HB Let + tert } OR fk (l2 4 BR, 2h) [erHbo et eat] 4 Dh (2 6 BAe 2 fag = ~ Bua 2 ag pa yfebot Heh AOR A 5 MD ye MA] [ha KG + RIG — (243 — 18)I2,]fe rH 4 eB] (24 AB + aki) [1G +e + bk] feo 4-1] A | ett [ea = ete] — (09 +1} — hoc.) [ (Kj + Kr — Hi] (+ ec) } inky fle Ry -(@+8 = aie elem essay (09 +18 + Kosh.) (12 +13 + kash 1B) A | —abycke[3 (42 + AQ) = fete Hk — er] — (2 4 KB — hyd) (ME FIG = Akai =) (> fe, = 2m [Pre + BG + (48 + BYR Jetrnato — wuss ata be, abet ig, = 2 he (lS + Oyen AR) +g (A + RIG — 2s + BIB, + HB) [eH Satoh et 5 ek] is, = ue ey) oly Figg = Bie up feted — ett (Kea “1 yeke[e HIh — De Bah eA at (2 Bye Ha — tte Dah eit eka mhegksg[ere — eis Mah 4 ei Ah] is, — WH 118 + afl is, = Ba (a vf where A is defined a follows A ieee rks) [ebb + EE — hak) (ete = + B(KE + Ka) kesh ferHtest0® = af (4 gy? + big fern ‘Appendix ¢ LINTRACK is a pavement tester in the Faculty of Civil Engineering and Geosciences, Delft University of Technology. As shown in Fis. Cl it consists of a free-rolling wheel that moves forward and backward with a guidance system. The force applied on the wheel can be varied between 15 and 100KN and the moving speed can be changed between 0 and 20km/h. A fully automatic electronic control system makes it possible to run LINTRACK continuously with automatic data collection. Various measuring instruments (eg. strain gauges) can be built into test sections to collect necessary information about the response of a pavement structure 2 sun, € Ksberzen and A. Sharpes ea finerntionl our o os and Srctres 180-81 (201) 45-61 8 References ‘ebenbach JD. 1858. Wave Propagation a laste Sois. vier Scene Pblaber pewhoury Re Kebegen,C, Sepa, A Sasuwendraadl 2002 Foodie pe. tel eee fr sate ropaaion and parameter senlieaton im mul-aer Se nt J Sole Set 3, 4072-4091 seAhoury Betray, A Kasbergen, Cy Bastwerdeis. J, 2001, Special element fsiaton ne) Soles Set 38, 1605-1623. Aenean Nee SH 203 Bondar en says of he eye Deve Jf 1897 Wave Frapegaton in Scares: Spcal Anais Using Fs Di van Dalen KN Teowaas A Metoine, AV; Hevng J, 2015, Tanto radaion seed by fee lose thas moves aver te interac of we elas layers I Psat Seve 75-74 9812 san 1965, The sec peduced in & serif cl bya moving sae Tare nt) Eng So St eat Hong. Hi, Yang. YB. 200 laste waver in veel halespace generated by aoa vee lade Si Dy tah Eng 20 17, Janes DV, Le Houtee D.Fepow. AT Pei, M1998, Ground vibration inthe Vocnty aft meving Nemane rectangle on halspace far) Mech Sede 150156 ae 2008 Spectral Element Method in Stet Dynami: John Wey & Sons Singapore Letewe-wesgoue. Ge Houde, D. Feplow. AT. 2000, Ground vibration in the vcinty of 2 highspeed moving Ramone sp les. | Sound Vib 23, basis Lefeoe-eszouer. Fela. AT Le Noukée,D, 2002. Sac vation doe Sequence af hi reee moving hams recanglr loses. Sei Dyn Har, ing 22 158083 Louhghla A. Aibaran.M. Ulm. Fl. 2013, Bosse’ coiectre:dsiation- Mernr eflecon neve pavementvenice interstion Eng, Meet 10 Banos Merribine. AV, 2004 Steady state response of an infinite sting on @ neni tar volute funastion fo moving po las.) Seund Vib" 272, 1035- ‘one Vorrottow AY. Merikine, AV. 2003, Praia ppeted bea on 2 visc-e- Tat Tae a aed! fr yam analyss af higheped raloay Uk I Sete ee 4057235552 van! Ko Sb TY. Man) D018. Spectral element method for dynamic ce ‘sponse of mullayece hai ecu sbjectee to hasnene moving Teed. I P'ckomect 8 #518181 You Yan K- Hi ¥ Lt J Ge 20%, Spcral element method for dynamic response of tanverselyiseopicasphat pavement uncer impact Tose Rose ater Pavement es 1.203038, ‘xa van Dalen, KN, Biba, JM. Metine, AV. 2016 Semanal sl Vin forte natu fespese of «neil suc erabedeee ih boo neous halespace In: neato spose an Enionmenal Vibration ane ‘rasperaton Geogvoanuc, Sinzapee 9p. 338-388, ‘ah SN. White. 985 Use of tneedesional. mame fine element Dorian for ana of ele pavement. Transp. Res. Re 1388 50-58

You might also like