Finite element procedures for nonlinear structures in moving coordinates. Part II: Infinite beam under moving harmonic loads Vu-Hieu Nguyen, Denis Duhamel To cite this version: Vu-Hieu Nguyen, Denis Duhamel. Finite element procedures for nonlinear structures in moving coordinates. Part II: Infinite beam under moving harmonic loads. Computers and Structures, Elsevier, 2008, 86 (21-22), pp.2056-2063. ฀hal-00345915฀ HAL Id: hal-00345915 https://hal.archives-ouvertes.fr/hal-00345915 Submitted on 14 Dec 2008 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Finite element procedures for nonlinear structures in moving coordinates. Part II: Infinite beam under moving harmonic loads Vu-Hieu Nguyen, Denis Duhamel∗ Laboratoire Analyse des Matériaux et Identification, Institut Navier, Ecole Nationale des Ponts et Chaussées, 6 & 8 avenue Blaise Pascal, Cité Descartes, Champs sur Marne, 77455, Marne-La-Vallée Cedex 2, France Abstract This paper presents a numerical approach to the stationary solution of infinite EulerBernoulli beams posed on Winkler foundations under moving harmonic loads. The procedure proposed in Part 1 [1], which has been applied to consider the longitudinal vibration of rods under constant amplitude moving loads in moving coordinates, is enhanced herein for the case of moving loads with time-dependent amplitudes. Firstly, the separation of variables is used to distinguish the convection component from the amplitude component of the displacement function. Then, the stationary condition is applied to the convection component to obtain a dynamic formulation in the moving coordinates. Numerical examples are computed with a linear structure to validate the proposed method. Finally, nonlinear elastic foundation problems are presented. Key words: moving load, moving coordinates, Euler-Bernoulli beam, stationary, finite element method Preprint submitted to Computers and Structures 17 April 2008 1 Introduction The coupled model of a beam resting on a Winkler foundation under moving loads has been widely studied in various applications, especially in railway engineering. In the literature, various beam models such as Euler-Bernoulli beam [2] or Timoshenko beam [3][4][5] subjected to constant or varied amplitude loads have been considered. For analyzing the linear responses of this problem, most of works use the (semi-)analytical methods, which are based on transformation techniques (e.g. Fourier or Laplace transforms) to solve the problem in the frequency or wavenumber domain. The analytical analysis allows to determine the critical velocity and critical frequency of the moving loads which cause resonances of the system. It has also been shown that, the dynamic responses of structures under moving loads are very different when the speed of the load is smaller (subsonic case) or greater (supersonic case) than the critical velocity. When taking into account the nonlinear effects, numerical procedures such as the finite element analysis (FEA) should be used. For nonlinear dynamic structures under moving loads, a transient analysis should be performed where the load varies both with time and space. Consequently, using FEA for moving load problems requires larger meshes than for fixed load problems [6][7]. Especially for resolving high velocity problems (e.g., the supersonic cases), the number of elements required for computation becomes very important. In order to overcome this inconvenience, an alternative approach has been proposed Email addresses: vuhieu.nguyen@ponts.org (Vu-Hieu Nguyen), duhamel@lami.enpc.fr (Denis Duhamel). 2 in which the dynamic equations are reformulated in a moving reference frame which goes with the load’s position by using a variable transformation [8]. The advantage of this approach is that it requires smaller finite element meshes as the load is fixed in the moving coordinates. However, it has been noticed that when the speed of the load reaches the critical velocity, the stiffness matrix formulated in the moving coordinates becomes ill-conditioned because of an important convected term which occurs from the variable transformation (this convected term depends on the speed of the load and on the mass of the structure) [1][9]. In the previous paper [1], we have proposed a novel finite element formulation for the problem of an infinite bar under a constant moving load. The finite element formulation has been obtained through 2 steps: the dynamic equation is firstly discretized in time by using the generalized-α method and a variable transformation is then applied to this time-discretized equation for obtaining a static equation in the moving coordinates. It has been shown that the condition number of the stiffness matrix derived from the proposed method is significantly improved in the supersonic case. The purpose of this work is the simulation of nonlinear dynamic phenomena in railway track structures due to movement of a high speed train. It deals with analyzing the stationary response of an infinite Euler-Bernoulli beam posed on a nonlinear Winkler foundation under a harmonic load moving with a constant velocity. Both subsonic and supersonic cases are interested. Based on the similar technique as presented in [1], we will propose an enhanced procedure to consider the moving harmonic load problem. Since the load amplitude is no longer constant, the equation established in the moving coordinates becomes a time-dependent equation. Some comparisons of the analytical and finite element solutions will be shown to validate the proposed procedure. Then, several 3 examples of nonlinear foundations will be presented to illustrate the nonlinear effects on dynamic responses of structures under moving subsonic/supersonic and constant/harmonic loads. 2 Statement of the problem Consider an infinite Euler-Bernoulli beam with constant bending stiffness EI and mass per unit length m. This beam is posed on a Winkler foundation with spring stiffness k and viscous damping η per unit length (see Fig. A.1). The spring stiffness is assumed to be a linear or nonlinear function of the vertical displacement. A concentrated load (f0 (t)) moves with a constant speed v on the beam along x-axis direction. Taking into account the gravity force, the equation of motion of this system is written as m(ẅ(x, t) + g) + η ẇ(x, t) + kw(x, t) + EIw(4) (x, t) = −f0 (t)δ(x − vt), (1) where w(x, t) denotes the vertical displacement of a point x at instant t; ẇ and ẅ denote respectively the first- and the second- order derivatives of w with respect to time (i.e. velocity and acceleration); w(4) denotes the fourth-order derivative of w with respect to x; g is a constant denoting the gravitational acceleration; δ denotes the Dirac’s delta function. Since only the stationary solution of (Eq. 1) is interested in this paper, a transform of this problem into the moving coordinates is suggested. However, due to the numerical difficulties associated with the supersonic load (see [1]), the direct variable transformation (x∗ := x − vt) will not be used. Instead, 4 another technique will be proposed in the next section. 3 Dynamic equation in moving coordinates This section presents a procedure to establish and to resolve (Eq. 1) in the moving coordinates. The function w(x, t) is firstly separated into two independent components, of which one depends only on the load’s position (x − vt) and the other depends only on time (t). Next, a time discretized equation of the first component will be formulated by using the generalized-α method [10]. Then, a variable transformation will be applied to this time-discretized equation to obtain a dynamic equation in moving coordinates by assembling the components in time and in space. The finite element formulation can finally be established based on the last equation. 3.1 Variable separation Assume that the solution w(x, t) has the form: w(x, t) = wv (x − vt)wt (t), (2) where the function wv depends only on the load’s position (x − vt) and the function wt depends only on time. Hence, the first- and second derivatives of w(x, t) are given by: ẇ(x, t) = ẇv wt + wv ẇt , ẅ(x, t) = ẅv wt + 2ẇv ẇt + wv ẅt . (3) (4) 5 Substitution of (Eq. 3) and (Eq. 4) into (Eq. 1) leads to an equation of wv : m̃ẅv + η̃ ẇv + k̃wv + j̃wv(4) = −f0 δv − g̃, (5) where: m̃ = mwt , η̃ = 2mẇt + ηwt , (6) (7) k̃ = mẅt + η ẇt + kwt , j̃ = EIwt , g̃ = mg , δv = δ(x − vt) (8) (9) (10) (11) 3.2 Time discretized equation of wv Using the generalized-α method [11], (Eq. 5) may be expressed in the discretized form as = −f0 δv(n+α) − g̃, m̃ẅv(n+α) + η̃ ẇv(n+α) + k̃wv(n+α) + j̃wv(4) (n+α) (12) where wv(n+α) , ẇv(n+α) , ẅv(n+α) , δv(n+α) are respectively the approximation values of wv , ẇv , ẅv , δv at tn+α ∈ [tn , tn+1 ] obtained from linear interpolations of their values at tn and tn+1 : wv(n+α) = (1 − αf )wv(n+1) + αf wv(n) , (13) ẇv(n+α) = (1 − αf )ẇv(n+1) + αf ẇv(n) , (14) δ(n+α) = (1 − αf )δv(n+1) + αf δv(n) . (16) ẅv(n+α) = (1 − αm )ẅv(n+1) + αm ẅv(n) , (15) where αm , αf are two constants. Using the Newmark scheme, wv(n+1) and ẇv(n+1) can be expressed in the approximative form as 6 β2 ∆t2 (1 − β2 )∆t2 ẅv(n) + ẅv(n+1) , 2 2 + (1 − β1 )∆tẅv(n) + β1 ∆tẅv(n+1) , wv(n+1) = wv(n) + ∆tẇv(n) + (17) ẇv(n+1) = ẇv(n) (18) where β1 , β2 are two integration parameters and ∆t := tn+1 − tn denotes the time step size. Substitution of (Eqs. 13-18) into (Eq. 12) leads to the following discretized equation at tn+1 Ã1 wv(n+1) + Ã2 wv(4) − B̃1 wv(n) − B̃2 ẇv(n) − B̃3 ẅv(n) + B̃4 wv(4) (n+1) (n) = F1 δ(x − vtn+1 ) + F2 δ(x − vtn ) + G, (19) where Ã1 , Ã2 , B̃1 , B̃2 , B̃3 , B̃4 , F1 , F2 , G are the coefficients which are defined in function of αm , αf , β1 , β2 , ∆t and m̃, η̃, k̃, j̃ as Ã1 = (1 − αm )m̃ + (1 − αf )β1 ∆t η̃ + (1 − αf ) Ã2 = (1 − αf ) β2 ∆t2 k̃ 2 β2 ∆t2 j̃ 2 (20) (21) β2 ∆t2 k̃ 2 ¸ · ∆t η̃ ∆t B̃2 = (1 − αm )m̃ − [β2 − 2(1 − αf )β1 ] 2 B̃1 = (1 − αm )m̃ + (1 − αf )β1 ∆t η̃ − αf B̃3 = [(1 − β2 − αm )m̃ − (1 − αf )(β2 − β1 )∆t η̃] β2 ∆t2 j̃ 2 β2 ∆t2 f0 (t) F1 (t) = (1 − αf ) 2 β2 ∆t2 F2 (t) = αf f0 (t) 2 β2 ∆t2 G= g 2 B̃4 = αf (22) (23) ∆t2 2 (24) (25) (26) (27) (28) 7 3.3 Variables transformation At the time step tn+1 , we take a variable transformation x∗ = x − vtn+1 . Since in moving coordinates, the function wv∗ (wv∗ := wv (x − vtn+1 )) in the steady state no longer depends on time, we have the relations: wv (x − vtn+1 ) = wv∗ (x∗ ) , wv (x − vtn ) = wv∗ (x∗ + v∆t) := wv∗ (x∗+ ) , ẇv (x − vtn ) = −vwv∗ ′ (x∗ + v∆t) := −vwv∗ ′ (x∗+ ) , ẅv (x − vtn ) = v 2 wv∗ ′′ (x∗ + v∆t) := v 2 wv∗ ′′ (x∗+ ) . (29) (30) (31) (32) Applying these relations (Eqs. 29-32) into (Eq. 19) yields Ã1 wv (x) + Ã2 wv(4) (x) − B̃1 wv (x+ ) + v B̃2 wv′ (x+ ) −v 2 B̃3 wv′′ (x+ ) + B̃4 wv(4) (x+ ) = F1 (t)δ(x) + F2 (t)δ(x+ ) + G . (33) in which the indexes “n + 1” and “∗ ” have been omitted for simplifying purposes. 3.4 Variable assembly We note that the coefficients Ã1 , Ã2 , B̃1 , B̃2 , B̃3 , B̃4 can be also expressed as Ã1 = A1 wt + A3 ẇt + A4 ẅt , Ã2 = A2 wt , B̃1 = B1 wt + B5 ẇt − B8 ẅt , B̃2 = B2 wt − B6 ẇt , (34) (35) (36) (37) B̃3 = B3 wt − B7 ẇt , B̃4 = B4 wt , (38) (39) 8 where the coefficients Ai (i = 1 ÷ 4), Bj (j = 1 ÷ 8) depend on the structural parameters and are given in Appendix A. Now using the variable assembly wv (x)wt (t) = w(x, t), (40) one can eliminate the terms related to wt in the expressions of Ã1 , Ã2 , B̃1 , B̃2 , B̃3 , B̃4 and then, express the approximative equation in the moving coordinates as h A1 w(x, t) + A2 w(4) (x, t) − B1 w(x+ , t) i +v B2 w′ (x+ , t) − v 2 B3 w′′ (x+ , t) + B4 w(4) (x+ , t) h i + A3 ẇ(x, t) − B5 ẇ(x+ , t) − v B6 ẇ′ (x+ , t) + v 2 B7 ẇ′′ (x+ , t) + [A4 ẅ(x, t) + B8 ẅ(x+ , t)] = F1 (t) δ(x) + F2 (t) δ(x+ ) + G , (41) This represents a dynamic equation of w(x, t) and of w(x + v∆t) wherein the load’s position is now fixed. Next, we formulate the finite element discretization for this equation. 4 Numerical implementation 4.1 Dynamic finite element equation in moving coordinates Consider a domain Ω = [−L, L] with the boundary condition assumption: w(x = −L) = w(x = L) = w,x (x = −L) = w,x (x = L) = 0. The variational 9 formulation of (Eq. 41) may be expressed as ¯ ¯ ¯ ¯ ¯ Find w(x, t) : W → R × T such as ∀δw ∈ δW : ¯ ¯ ¯ ¯ ¯ ·D E D E ¯ ′′ ′′ ¯ A w(x, t), δw(x) + A w (x, t), δw (x) ¯ 1 2 Ω Ω ¯ ¯ ¯ ¯ D E D E ¯ ¯ − B1 w(x+ , t), δw(x) + vB2 w′ (x+ , t), δw(x) + ¯ Ω Ω ¯ ¯ ¯ ¯ ¯ D E D E ¸ ¯ 2 ′ ′ ′′ ′′ + v B w (x , t), δw (x) + B w (x , t), δw (x) ¯ 3 + 4 + ¯ Ω Ω ¯ ¯ ¯ ·D ¯ E D E ¯ ¯ + A3 ẇ(x, t), δw(x) − B ẇ(x , t), δw(x) 5 + ¯ Ω Ω ¯ ¯ ¯ ¯ D E ¸ D E ¯ 2 ′ ′ ′ ¯ − vB6 ẇ (x+ , t), δw(x) − v B7 ẇ (x+ , t), δw (x) ¯ Ω Ω ¯ ¯ ¯ ·D ¯ D E ¸ E ¯ ¯ + A4 ẅ(x, t), δw(x) + B ẅ(x , t), δw(x) 8 + ¯ Ω Ω ¯ ¯ ¯ ¯ D E D E ¯ ¯ , + G, δw(x) = F (t)δ(x) + F (t)δ(x ), δw(x) 1 2 + ¯ Ω Ω ¯ (42) where W and δW are the collections of the trial solutions and of the weighting solutions, respectively; the notation < ., . >Ω is defined by: < a, b >Ω = RL −L ab dx. The boundary conditions are already taken into account in this variational formulation. Performing the discretisation of (Eq. 42), one can obtain the dynamic finite element equation in the classical form Mẅ(t) + Cẇ(t) + Kw(t) = F(t), w(0) = 0 ; ẇ(0) = 0, (43) (44) where M, C, K denote respectively the global mass, damping and stiffness matrices; w(t), ẇ(t), ẅ(t) denote respectively the global nodal displacement, 10 velocity and acceleration vectors; F(t) denotes the exterior nodal force vector. This presents a system of differential equations with given initial conditions that may be solved by using direct time integration algorithms, e.g., classical implicit Newmark scheme [7,12] or the advanced scheme when including nonlinearities [11,13,14]. Using a vector of interpolation function N(x), one can formulate explicitly the elementary matrices Me , Ce , Ke and the force vector Fe of element e in [xi , xi+1 ] as follows: ³ ´ Me = Nt A4 N ³ ´Ω Ce = Nt A3 N e ³ t ³ t + vN ³ Ωe ´ K = N A1 N Ωe ³´ Ωe := ´ R xi+1 xi ³ − Nt B5 N+ ³ ′′ t ´ Ωe ³ ´ ´Ω + v N ³ ´ ′t 2 e , Ωe ′′ + N A2 N B2 N′+ e Ω Fe = Nt F1 δ(x) where . e ³ + Nt B8 N+ (45) ³ − vNt B6 N′+ Ωe ³ t − N B1 N+ ´ B3 N′+ e Ω ´ + Nt F2 δ(x+ ) Ωe ³ ´ Ωe ´ Ωe ′′ t ³ − v 2 N′ B7 N′+ (46) , e + N B4 N′′+ ³ + Nt G ´ Ωe , t ´ ´ Ωe , Ω (47) (48) (.)dx; x+ := x + v∆t and N+ := N(x + v∆t). Herein, the Hermite’s interpolation function [12] that is appropriate to the Euler-Bernoulli beam assumption can be used. Remark 1 In linear cases, analytical formulations of the elementary matrices can be expressed explicitly (see [15]). Otherwise, for nonlinear structures, numerical quadratures by Gauss points should be used to compute these matrices. Remark 2 If the load amplitude is constant, velocity and acceleration terms in moving coordinates should be zeros and one can check that (Eq. 43) is reduced to the static equation form Kw = F. 11 (49) 4.2 Parameter choices Assuming that the beam is discretized into elements with equal lengths, parameter studies show that the best approximation is obtained with the following values of the αm , αf , β1 , β2 , ∆t: αm = 0.25, αf = 0.5, β1 = 0.5, β2 = 0.5, ∆t = h/v where h denotes the size of an element. The error when using these parameters is estimated by second order function of element size h as follows: h2 v (3) 3m ′ e(h) = ηw(x, t) + 2mẇ(x, t) − D w (x, t) + O(h3 ), 6 2 µ ¶ (50) These values of parameters will be used for all numerical examples in the following section. 5 Absorbing boundary condition When simulating the wave propagation in infinite structures by using FEM with the usual finite boundary conditions, the modeled size must be sufficiently large in order to insure that the outgoing waves are totally attenuated before arriving at the boundaries. To reduce the mesh size without disturbing the solution due to the reflected wave from the boundaries, one may add artificial boundaries to absorb all outgoing energies. For this problem, we propose adding at two ends of the beam two layers in which viscous damping is introduced in a manner that all incident waves should be absorbed. Denoting the depth of each layer with δ, the viscous 12 damping introduced in this layer is assumed to be a function of x (Fig. A.2): η̃(x) =        η  if |x| < L0 ,  !4 Ã    |x| − L  0   η + γ δ (51) if |x| > L0 , where η, γ are two constant parameters chosen so that: (1) the damping in the absorbing layer should be sufficiently important so that the outgoing waves can be attenuated in this layer, and (2) the evolution of the damping along this layer absorbing layer should be sufficiently smooth so that there is no reflection waves caused by the suddenly damping changes. Figure (A.3) presents an example to demonstrate the feature of absorbing layers by considering an infinite Euler beam (EI = 108 N/m2 ) posing on a linear Winkler foundation without damping (k = 103 N m−2 , η = 0). This beam is subjected to an unit harmonic load applying at x = 0: f (t) = f0 eiωt , ω = 5.8 rad/s, f0 = 1 N . A domain [−250m, 250m] is used for finite element modeling and it is meshed by 100 uniform elements with the fixed boundary conditions at two ends. The absorbing layer depth is δ = 150m at each end. It may be checked that the absorbing layers allow to attenuate perfectly the waves propagating away from the load’s position and, hence, the solution in the interested zone is in good agreement with the analytical one [2]. When the system has no damping, the numerical solution obtained without using absorbing layers is completely incorrect (see the dashed curve in Fig. A.3) due to the reflected waves from the fixed boundaries. 13 6 Validation 6.1 Example 1: Linear cases The proposed formulations will be firstly validated with a linear problem in which analytical solutions are available for all velocities and frequencies of the load (see e.g [2] or [8]). As the structure is linear, using the transformation technique show that the critical values of the velocity (v0 ), of the frequency (Ω0 ) and of the damping (η0 ) depend on the structural characteristics [2] as follows v0 = s 4 4EIk ; (m)2 Ω0 = s k ; m √ η0 = 2 mk. (52) Let consider an European rail for high-speed trains with the following characteristics: section area S = 76.86 cm2 , section inertia I = 3060 cm4 , mass density ρ = 7850 kg/m3 , Young’s modulus E = 200 × 109 N/m2 , mass per length m = 60.34 kg/m, flexural rigidity EI = 6.12 × 106 N m2 . The stiffness of the Winkler foundation (k) is 1.6 × 107 N/m and its damping coefficient η is assumed to be equal to 0.1η0 . A concentrated harmonic force f0 (t) = P sin ωt (P = 10 kN ) is applied downward on the rail and moves rightward at constant speed v. The gravity is taken into account. Four load cases are taken into consideration: (a) constant subsonic load (v = 0.5v0 , ω = 0), (b) constant supersonic load (v = 1.5v0 , ω = 0), (c) harmonic subsonic load (v = 0.5v0 , ω = 0.5Ω0 ), (d) harmonic supersonic load (v = 1.5v0 , ω = 0.5Ω0 ). As mentioned earlier, if the load amplitude is constant, only resolution of a static system (49) is required to obtain the stationary response. Otherwise, 14 the direct time integration is required for moving harmonic load problems. In this example, the classical Newmark implicit scheme has been performed for resolving the system of transient equations (43). The integration parameters are: time step ∆t = 0.05 T (T = 2π/ω), number of time steps nt = 200, tmax = nt × ∆t; β1 = β2 = 0.5. Figures (A.4-A.9) present the numerical results obtained by using the proposed finite element formulation (Eqs. 43-48) in comparison with corresponding analytical solutions. One can state that the numerical results have good agreement with the analytical ones. Note that in cases (b) (Fig. A.6) and (d) (Fig. A.8), the valid solutions are found in zone [−20m, 20m]. The no-exact solutions outside this zone in two supersonic cases represent the effect of the absorbing layers that absorb the outgoing waves. The finite element meshes used for modeling the structure in all cases are presented in Table A.1. In the subsonic cases (Figs. A.4 and A.6), the waves are attenuated rapidly and the computations require neither a large number of elements nor large meshes. For the contrary, when v > v0 (supersonic), the waves are attenuated slowly when they propagate farther away from the load position (Figs. A.5 and A.8), therefore the absorbing layers would be used. The sizes of the absorbing layers used for each case are given in the Table A.1. Moreover, the displacement solutions in front of the load position (x > 0) has higher frequencies than those in the subsonic cases, so denser meshes are required for obtaining correct results. We depict also the displacements at point x = 0 versus time in the harmonic load cases to verify the fact that the stationary state may be obtained after several load time periods (Figs. A.7 and A.9). 15 Also, the well-known Doppler effect can be clearly observed in these figures, especially in the supersonic cases. 6.2 Example 2: nonlinear cases This example considers a weak nonlinear elastic foundation problem by assuming that the reaction force of the foundation R(w) is a cubic polynomial function of the vertical displacement of the beam (Fig. A.10): ³ ´ R(w) = k0 1 + γ0 w2 w, (53) where k0 and γ0 are two constant parameters. For the purpose of comparison, these parameters are chosen so that R(w)|w=W0 = kw|w=W0 (k is the spring stiffness in the linear case). Then we take k0 = k/2, γ0 = 1/W02 and W0 = 4 × 10−4 m. Other structural parameters used for this analysis are the same as in the linear case (Example 1). The gravity force is not considered in this case. The analyses are performed for 4 cases of loads similar to Example 1. The reference critical values v0 and Ω0 are the same as in the linear case. The finite element meshes are presented in Table A.1. In order to ensure the numerical stability, which is more difficult to obtain for nonlinear dynamic problems, we have used the α-generalized method [12] for the direct time integration (Eq. 43). The numerical parameters used in this example are: time step ∆t = 0.05 T (T = 2π/ω), number of time steps nt = 200, tmax = nt × ∆t; α = 0.01, β1 = β2 = 0.5. The displacements of the beam due to constant subsonic/supersonic moving loads are plotted in Figures (A.11-A.12). One can state that the displacement 16 amplitudes are only a little greater for the nonlinear structure in comparison with linear results. However, in the supersonic case, the wave lengths are very different from those in the subsonic cases. Figures (A.13-A.14) show important effects of the foundation nonlinearities on amplitude as well as on wavelength of displacements of the beam under harmonic moving loads. The displacements at point x = 0 versus time in both subsonic and supersonic cases presented in Figures (A.15) and (A.16) show that numerical solutions of nonlinear problem are stable and the stationary response of these nonlinear problems can be reached. 7 Conclusion In this paper, a numerical approach has been presented for studying dynamic responses of an infinite Euler-Bernoulli beam posed on a nonlinear Winkler foundation in a moving reference frame. In principle, one can apply a variable transformation to establish the problem in moving coordinates. In the moving coordinates, the stiffness matrix is modified to take into account the convective component that derived from the variable transformation. However, higher velocities of loads make the convective component more important, and cause the stiffness matrix to be ill-conditioned. When using the proposed approach to establish the finite element formulation, the convective component is moved partly out of the diagonal of the stiffness matrix. As a consequence, the system could be well-conditioned even when the load moves at supersonic velocities. The proposed formulation in moving co-ordinates greatly reduces the volume 17 of computation tasks, by taking advantage of the fact that the structure can be meshed similarly to the fixed load problems. Particularly, the problem becomes static when the load amplitude is constant. The numerical examples presented show that this procedure can be applied to resolve both linear and nonlinear problems with all values of velocity and frequency of loads. Significant dynamic effects has been stated when taking into account the nonlinearities of the foundation. For these simulations, the numerical parameters have chosen by assuming that the element sizes are identical. Further studies on the parameter choices should be carried out in case of a mesh with variable element sizes is needed. Similar formulations could be extended for analyzing other structures under moving loads, for example, when considering one-dimensional structures with numerous DOF [15] or 2D/3D solid structures. 18 A Analytical expressions of coefficients A1 , A2 , A3 , A4 , B1 , B2 , B3 , B4 , B5 , B6 , B7 , B8 . β2 ∆t2 k A1 = (1 − αm )m + (1 − αf )β1 ∆t η + (1 − αf ) 2 β2 ∆t2 A2 = (1 − αf ) EI 2 β2 ∆t2 A3 = 2(1 − αf )β1 ∆tm + (1 − αf ) η 2 β2 ∆t2 A4 = (1 − αf ) m 2 β2 ∆t2 B1 = (1 − αm )m + (1 − αf )β1 ∆t η − αf k 2 ¸ · ∆t η ∆t B2 = (1 − αm ) m − [β2 − 2(1 − αf )β1 ] 2 ∆t2 B3 = [(1 − β2 − αm ) m − (1 − αf )(β2 − β1 )∆t η] 2 β2 ∆t2 EI B4 = αf 2 β2 ∆t2 B5 = 2(1 − αf )β1 ∆tm − αf η 2 B6 = [β2 − 2(1 − αf )β1 ] ∆t2 m B7 = (1 − αf )(β2 − β1 ) ∆t3 m β2 ∆t2 m B8 = αf 2 (A.1) (A.2) (A.3) (A.4) (A.5) (A.6) (A.7) (A.8) (A.9) (A.10) (A.11) (A.12) References [1] V.-H. Nguyen, D. Duhamel, Finite element procedures for nonlinear structures in moving coordinates. Part I: infinite bar under moving axial loads, Computers and Structures 84 (21) (2006) 1368–1380. [2] L. Frýba, Vibration of Solids and Structures under Moving Loads, 3rd Edition, Thomas Telford, 1999. 19 [3] Y.-H. Chen, Y.-H. Huang, Dynamic stiffness of infinite Timoshenko beam on viscoelastic foundation in moving co-ordinate, International Journal for Numerical Methods in Engineering, 48 (2000) 1–18. [4] J.-F. Hamet, Railway noise : Use of Timoshenko model in rail vibration studies, Acta Acustica, 85 (1999) 54–62. [5] M. Hardy, The generation of waves in infinite structures by moving harmonic loads, Journal of Sound and Vibration 180 (4) (1995) 637–644. [6] J. Rieker, Y.-H. Lin, M. Trethewey, Discretisation consideration in moving load finite element beam models, Finite Element in Analysis and Design, 21 (1996) 129–144. [7] K. J. Bathe, Finite Element Procedures, Prentice Hall, 1996. [8] L. Andersen, S. Nielsen, P. Kirkegaard, Finite element modelling of infinite Euler beams on Kelvin foundations exposed to moving loads in convected coordinates, Journal of Sound and Vibration, 241 (4) (2001) 587–604. [9] V.-H. Nguyen, D. Duhamel, B. Nedjar, Calcul de structures soumises à une charge mobile par la méthode des éléments finis, XVème Congrès Français de Mécanique, Nancy, 3–7 Septembre, 2001 . [10] H. Hilber, T. Hughes, R. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Engineering & Structural Dynamics, 5 (1977) 283–292. [11] J. Chung, G. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation : the generalized-α method, ASME Journal of Applied Mechanics, 60 (1993) 371–375. [12] T. Hughes, The Finite Element Method - Linear Static and Dynamic Finite Element Analysis, Dover Publications, 2000. 20 [13] K. J. Bathe, M. M. I. Baig, On a composite implicit time integration procedure for nonlinear dynamics, Computers & Structures 83 (2005) 2513–2534. [14] K. J. Bathe, Conserving energy and momentum in nonlinear dynamics: A simple implicit time integration scheme, Computers & Structures 85 (2007) 437–445. [15] V.-H. Nguyen, Comportement dynamique de structures nonlinéaires soumises à des charges mobiles, Ph.D. thesis, Ecole Nationale des Ponts et Chaussées (Paris) (2002). 21 List of Figures A.1 Infinite Euler beam resting on Winkler foundation 25 A.2 Function η̃ in the absorbing layer 26 A.3 Absorbing layer effect: fixed harmonic load problem 26 A.4 Displacement of the beam posed on linear foundation under subsonic constant load 27 A.5 Displacement of beam posed on linear foundation under supersonic constant load 27 A.6 Displacement of beam (at t = tmax ) posed on linear foundation under subsonic harmonic load 28 A.7 Displacement in time at point x = 0 of beam posed on linear foundation under subsonic harmonic load 28 A.8 Displacement of beam (at t = tmax ) posed on linear foundation under a supersonic harmonic load 29 A.9 Displacement in time at point x = 0 of the beam posed on a linear foundation under supersonic harmonic load A.10 Reaction of the nonlinear elastic foundation 29 30 A.11 Displacement of beam posed on nonlinear elastic foundation under subsonic constant load 22 31 A.12 Displacement of beam posed on nonlinear elastic foundation under supersonic constant load 31 A.13 Displacement of beam (at t = tmax ) posed on nonlinear elastic foundation under a subsonic harmonic load 32 A.14 Displacement of beam (at t = tmax ) posed on nonlinear elastic foundation under supersonic harmonic load 32 A.15 Displacement in time at point x = 0 of beam posed on nonlinear elastic foundation under subsonic harmonic load 33 A.16 Displacement at point x = 0 of beam posed on nonlinear elastic foundation under supersonic harmonic load 23 33 List of Tables A.1 Computational parameters for linear cases 34 A.2 Computational parameters for nonlinear cases 34 24 z −∞ x f0 (t)δ(x − vt) +∞ k η Fig. A.1. Infinite Euler beam resting on Winkler foundation 25 η̃ η̃ = η + γ ¡ x−L0 ¢4 δ η L0 x δ absorbing layer Fig. A.2. Function η̃ in the absorbing layer −5 2 x 10 analytical FEM 1.5 solution obtained without using absorbing condition vertical displacement (m) 1 0.5 0 −0.5 −1 −1.5 absorbing layer absorbing layer −2 −250 −200 −150 −100 −50 0 X (m) 50 100 150 200 250 Fig. A.3. Absorbing layer effect: fixed harmonic load problem 26 v = 0.5 v0 ; ω = 0 −4 0.5 x 10 vertical displacement (m) 0 −0.5 −1 −1.5 −2 −2.5 −3 FEM analytical −10 −5 0 X (m) 5 10 Fig. A.4. Displacement of the beam posed on linear foundation under subsonic constant load v = 1.5 v0 ; ω = 0 −4 4 x 10 FEM analytical vertical displacement (m) 3 2 1 0 −1 −2 −3 −4 −40 −30 −20 −10 0 X (m) 10 20 30 40 Fig. A.5. Displacement of beam posed on linear foundation under supersonic constant load 27 v = 0.5 v0, ω = 0.5Ω0 −4 2 x 10 FEM analytical vertical displacement (m) 0 −2 −4 −6 −15 −10 −5 0 X (m) 5 10 15 Fig. A.6. Displacement of beam (at t = tmax ) posed on linear foundation under subsonic harmonic load −4 displacement at point x = 0 (m) 6 x 10 FEM analytical 4 2 0 −2 −4 −6 0 0.05 0.1 0.15 0.2 0.25 time (s) Fig. A.7. Displacement in time at point x = 0 of beam posed on linear foundation under subsonic harmonic load 28 v = 1.5 v0, ω = 0.5Ω0 −4 1 x 10 FEM analytical 0.5 vertical displacement (m) 0 −0.5 −1 −1.5 −2 −2.5 −3 −3.5 −40 −30 −20 −10 0 X (m) 10 20 30 40 Fig. A.8. Displacement of beam (at t = tmax ) posed on linear foundation under a supersonic harmonic load −4 x 10 FEM analytical displacement at point (x=0) (m) 1 0 −1 0 0.05 0.1 0.15 0.2 0.25 time (s) Fig. A.9. Displacement in time at point x = 0 of the beam posed on a linear foundation under supersonic harmonic load 29 R(w) k0 (1 + γw2 )w k k0 W0 w Fig. A.10. Reaction of the nonlinear elastic foundation 30 v = 0.5 v0 ; ω = 0 −4 2 x 10 nonlinear linear vertical displacement (m) 1 0 −1 −2 −3 −4 −10 −5 0 X (m) 5 10 Fig. A.11. Displacement of beam posed on nonlinear elastic foundation under subsonic constant load v = 1.5 v0 ; ω = 0 −4 4 x 10 nonlinear linear vertical displacement (m) 3 2 1 0 −1 −2 −3 −4 −40 −30 −20 −10 0 X (m) 10 20 30 40 Fig. A.12. Displacement of beam posed on nonlinear elastic foundation under supersonic constant load 31 v = 0.5 v0 ; ω = 0.5 Ω0 −4 vertical displacement (m) 2 x 10 nonlinear linear 1 0 −1 −2 −3 −4 −40 −20 0 X (m) 20 40 Fig. A.13. Displacement of beam (at t = tmax ) posed on nonlinear elastic foundation under a subsonic harmonic load v = 1.5 v0 ; ω = 0.5 Ω0 −4 vertical displacement (m) 2 x 10 nonlinear linear 1 0 −1 −2 −3 −4 −60 −40 −20 0 X (m) 20 40 60 Fig. A.14. Displacement of beam (at t = tmax ) posed on nonlinear elastic foundation under supersonic harmonic load 32 −3 x 10 v = 0.5 v0 ; ω = 0.5 Ω0 displacement at point (x = 0) (m) 1 nonlinear linear 0.5 0 −0.5 −1 0 0.05 0.1 0.15 0.2 0.25 time (s) Fig. A.15. Displacement in time at point x = 0 of beam posed on nonlinear elastic foundation under subsonic harmonic load −4 x 10 v = 1.5 v0 ; ω = 0.5 Ω0 displacement at point (x = 0) (m) 1 nonlinear linear 0.5 0 −0.5 −1 0 0.05 0.1 0.15 0.2 0.25 time (s) Fig. A.16. Displacement at point x = 0 of beam posed on nonlinear elastic foundation under supersonic harmonic load 33 Table A.1 Computational parameters for linear cases load type num. of elements domain length absorbing layers subsonic constant load 50 20 0 supersonic constant load 400 80 [−40m, −20m] & [20m, 40m] subsonic harmonic load 90 30 0 supersonic harmonic load 400 80 [−40m, −20m] & [20m, 40m] num. of elements domain length absorbing layers subsonic constant load 50 20 0 supersonic constant load 400 120 [−60m, −30m] & [30m, 60m] subsonic harmonic load 200 80 0 supersonic harmonic load 400 120 [−60m, −30m] & [30m, 60m] Table A.2 Computational parameters for nonlinear cases load type 34