- Email: [email protected]
A closed-form solution of beam on viscoelastic subgrade subjected to moving loads
Computers and Structures 80 (2002) 1–8 www.elsevier.com/locate/compstruc
A closed-form solution of beam on viscoelastic subgrade subjected to moving loads Lu Sun * ECJ Hall 6.10, Department of Civil Engineering, The University of Texas at Austin, Austin, TX 78712, USA Received 2 June 2000; accepted 30 September 2001
Abstract A closed-form solution of beam on viscoelastic subgrade subjected to moving loads is proposed in this paper. The Green’s function of the beam is obtained by means of Fourier transform. Theory of linear partial differential equation is used to express the deflection of the beam in terms of the Green’s function. Beam’s deflection is represented as a generalized integral using inverse Fourier transform. To evaluate this generalized integral analytically, poles of the integrand are identified using the theory of algebraic equation. Every pole and its order are isolated and given in a closed form. The theorem of residue is then applied to represent the generalized integral in the form of contour integral in the complex plane. Closed-form deflection and numerical computation of distinct cases are provided as an illustration of the proposed solution. Ó 2002 Published by Elsevier Science Ltd. Keywords: Beam; Moving loads; Viscoelastic subgrade; Green’s function; Fourier transform; Complex function; Residue theorem
1. Introduction The response of a beam to moving loads has been studied extensively over the past several decades [7]. Considerable studies have been conducted to investigate the deflection response of beams [1,2,6,8,11,17–19]. A growing interest on this topic is recently found in railway industry and highway industry due to the reason that this kind of physical model can be used as an idealized representation of dynamic interactions of railtrack/train and pavement/vehicle [9,10,14,16,21,22, 24,25]. An interesting manifestation of this subject is that when the subgrade is modeled as an elastic subgrade, a critical load velocity is found existing, which may cause a significant variation of a beam’s response [11]. Waves generated by a moving load with a supersonic velocity excite different patterns of response of the beam as they do when the load velocity is subsonic. This critical velocity is also predicted and observed by other researchers
by means of either electrical analogy using electrical device and experimental measurement [4] or theoretical analysis [16,20]. Since high-speed vehicles are getting more widely adopted as surface transportation carriers, considerable attention has been paid to the response of infrastructure under high-speed loading condition since 1990s [5,12,15,20,23]. In most previous studies the subgrade of a beam is usually modeled as an elastic Winkler foundation. Moving at the critical speed, the load will excite the beam to form an infinite deflection [11,20]. This idealization of subgrade is, however, not quite true because any physical system in reality possesses damping, which leads to energy dissipation. It can be expected that the unrealistic infinite deflection will disappear if the damping effect is included in the mathematical model of a beam subjected to moving loads. In this paper the effect of damping on the deflection response is taken into account. A closedform representation of beam deflection is presented. 2. The Green’s function
*
Tel.: +1-512-476-1906. E-mail address: [email protected] (L. Sun).
The governing equation of a flexible beam on a viscoelastic subgrade can be given by [22]
0045-7949/02/$ - see front matter Ó 2002 Published by Elsevier Science Ltd. PII: S 0 0 4 5 - 7 9 4 9 ( 0 1 ) 0 0 1 6 2 - 6
2
L. Sun / Computers and Structures 80 (2002) 1–8
EI
o4 y oy o2 y þ Ky þ C þ m 2 ¼ F ðx; tÞ 4 ox ot ot
ð1Þ
where yðx; tÞ represents the deflection of the beam, x represents the traveling direction of the moving load, and t represents time. Also, EI is the rigidity of the beam, E is Young’s modulus of elasticity, I is the moment of inertia of the beam, m is the unit mass of the beam, and F ðx; tÞ is external loads. A moving concentrated load can be given by F ðx; tÞ ¼ P dðx vtÞ
ð2Þ
where P is the amplitude of the applied load, and dðÞ is the Dirac-delta function, which is defined by Z 1 dðx x0 Þf ðxÞ dx ¼ f ðx0 Þ ð3Þ 1
Assume that the beam is at rest initially (i.e., initial deflection, velocity, and acceleration of the beam is zero). Eqs. (2) and (3) constitute the mathematical description of the problem of interest. According to the theory of linear partial differential equation [13,23], the Green’s function represents the fundamental solution of the equation in which the load is given by the Dirac-delta function. More specifically, the Green’s function of this problem is defined as the solution of Eq. (1) given that the applied load is characterized by Fd ðx; tÞ ¼ dðx x0 Þdðt t0 Þ
F½f ðnÞ ðtÞ ¼ ðixÞn F½f ðtÞ
Since the right-hand side of Eq. (7) Fed ðn; xÞ is the representation of Fd ðx; tÞ in the frequency domain, it is necessary to evaluate Fd ðx; tÞ in the frequency domain too. This can be achieved by taking two-dimensional Fourier transform on both sides of Eq. (4) Z 1 Z 1 Fed ðn; xÞ ¼ dðx x0 Þdðt t0 Þ 1
1
ð5Þ
f ðx; tÞ ¼ F1 ½f~ðn; xÞ Z 1 Z ¼ ð2pÞ2 1
1
e ðn; x; x0 ; t0 Þ ¼ exp½iðnx0 þ xt0 Þ ðEIn4 þ K G þ iCx mx2 Þ1
ð10Þ
The Green’s function given by Eq. (10) is represented in the frequency domain and must be converted to the time domain. To this end, taking the inverse Fourier transform on both sides of Eq. (10) gives Z 1 Z 1 Gðx; t; x0 ; t0 Þ ¼ ð2pÞ2 1
1
expfi½nðx x0 Þ þ xðt t0 Þ g
dn dx EIn4 þ K þ iCx mx2 ð11Þ Formula (11) is the Green’s function of the beam resting on a viscoelastic subgrade. The Green’s function serves as a fundamental solution of a partial differential equation. It is very useful when dealing with linear systems.
3. Integral representation of the solution f~ðn; xÞ exp½iðnx þ xtÞ dn dx
1
ð6Þ where F½ and F1 ½ are the Fourier transform and its inversion, respectively. To solve the Green’s function, taking two-dimensional Fourier transform on both sides of Eq. (1) gives 4e
ð9Þ
Here, the property of the Dirac-delta function, i.e., Eq. (3), is utilized to evaluate the integration in Eq. (9). Substituting the result of Eq. (9) into Eq. (7) and realizing that Eq. (7) is an algebraic equation, it is straightforward to see
ð4Þ
1
1
exp½iðnx þ xtÞ dx dt ¼ exp½iðnx0 þ xt0 Þ
Define two-dimensional Fourier transform and its inversion as f~ðn; xÞ ¼ F½f ðx; tÞ Z 1 Z 1 ¼ f ðx; tÞ exp½iðnx þ xtÞ dx dt
ð8Þ
e þ iCx G e mx2 G e ¼ Fed ðn; xÞ EIn G þ K G
ð7Þ
in which Fed ðn; xÞ is the Fourier transform of Fd ðx; tÞ. In Eq. (7) the deflection response yðx; tÞ has been replaced e ¼G e ðn; x; x0 ; t0 Þ to indicate the Fourier by the symbol G transform of the Green’s function. Also, the following property of Fourier transform is used in the derivation of Eq. (7)
According to Sun and Greenberg [23], the solution of Eq. (2) can be constructed by integrating the Green’s function in all dimensions Z t Z 1 F ðx0 ; t0 ÞGðx; t; x0 ; t0 Þ dx0 dt0 ð12Þ yðx; tÞ ¼ 1
1
Taking Eqs. (4) and (11) into Eq. (12) gives Z t Z 1 Z 1 Z 1 yðx; tÞ ¼ 1
1
1
1
P dðx0 vt0 Þ expfi½nðx x0 Þ þ xðt t0 Þ g
ð2pÞ2 EIðn4 þ K þ iCx mx2 Þ
dn dx dx0 dt0 where K ¼ K=EI and m ¼ m=EI. Note that
ð13Þ
L. Sun / Computers and Structures 80 (2002) 1–8
Z
1
dðx0 vt0 Þ expðinx0 Þ dx0 ¼ expðinvt0 Þ
ð14Þ
exp½iðx þ vnÞt0 dt0 ¼ 2pdðx þ vnÞ
ð15Þ
1
Z
t 1
Substituting Eqs. (14) and (15) into Eq. (13) and reapplying the property (3) one gets Z 1 P exp½inðx vtÞ dn ð16Þ yðx; tÞ ¼ 2pEI 1 n4 mv2 n2 þ K iCvn So far the integral representation of the dynamic deflection of the beam under a moving concentrated load has been obtained. Expression (16) can be further developed using complex analysis. In the following section, the theorem of residue is employed to evaluate the integration (16).
Before the integration (16) is further evaluated, it is necessary to determine the roots of the characteristic equation ð17Þ
Assume none of these parameters m, K, C and v being zero. Essentially, one seeks zeros of the polynomial at the left-hand side of the characteristic equation (17). Because several algebraic theorems [3] are used in the derivation of roots of Eq. (17), they are cited and stated below. Theorem 1 (fundamental theorem of algebra).P If pðnÞ is a polynomial of degree n P 1, that is, pðnÞ ¼ ni¼0 ai ni , with coefficients ai ði ¼ 1; 2; . . . ; nÞ real or complex numbers and an 6¼ 0, then pðnÞ has at least one zero; i.e., there exists a complex number n such that pðnÞ ¼ 0. Theorem 2. A polynomial of degree n P 1 has at most n zeros, counting multiplicity. Theorem 3 (Descartes’ rule of signs). The number np of positive real zeros of a polynomial pðnÞ is less than or equal to the number of variations nv in sign of the coefficients of pðnÞ. The number of negative real zeros of pðnÞ is at most equal to the number of sign changes in the coefficients of the polynomial pðnÞ. Moreover, the difference nv np is an even integer. Theorem 4. Every zero of a polynomial of degree n P 1 lies in the circular region defined by jnj 6 q, where the radius q is given by q ¼ 1 þ max jai =an j. 0 6 i 6 n1
Define abbreviations a0 ¼ K, a1 ¼ Cv, a2 ¼ mv2 . Due to Theorem 4 all zeros of the left-hand side of Eq. (17) (i.e., pðnÞ ¼ n4 a2 n2 ia1 n þ a0 ) lie in the open disk centered at the origin of the complex plan and whose radius q is given by q ¼ 1 þ max0 6 i 6 2 jai j. This result may be instrumental in locating the zeros when the numerical methods are used to compute the roots of Eq. (17). Since a coefficient of Eq. (17) is imaginary, according to Theorems 1 and 2, no real-valued roots exist for Eq. (17) and all four roots must be complex. In the following only imaginary roots and complex roots are considered analytically. Imaginary roots are discussed separately because in this case duplicated roots may exist, which is very crucial in determining the order of the poles in the generalized integral (16). 4.1. The condition of existing imaginary roots for Eq. (17) To derive the condition of exiting imaginary roots of Eq. (17), define the imaginary root n ¼ in0 ðn0 6¼ 0 and n0 2 RÞ. Substitution of the imaginary root into Eq. (17) gives
4. Nature of zeros of the characteristic equation
n4 mv2 n2 þ K iCvn ¼ 0
3
fðn0 Þ ¼ n40 þ a2 n20 þ a1 n0 þ a0 ¼ 0
ð18Þ
The objective here is to seek real roots of Eq. (18) and the condition under which these real roots exist. According to Theorem 3 one may conclude that two possibilities arise: (1) the polynomial f has two real negative zeros (could be duplicated) and one complex conjugate pair or (2) two complex conjugate pairs. Define polynomials f ðn0 Þ ¼ n40 and fj ðn0 Þ ¼ a2 n20 a1 n0 a0 ¼ 2 a2 ðn0 þ ðC=2mvÞÞ2 þ ðC =4mÞ K, ðj ¼ 1–3Þ with respect to n0 . These polynomials are plotted in Fig. 1. The first possibility corresponds to polynomials f1 and f2 , while the second possibility corresponds to polynomial f3 . Since n0 must be real, only the first possibility is considered. In order for (1) to be the case it is necessary and sufficient that the graph of f crosses the n0 axes twice and has a negative minimum at a point g < 0. Then, the polynomial gðn0 Þ ¼ fðn0 Þ0 ¼ 4n30 þ 2a2 n0 þ a1 must have a negative zero g, which is the case due to Theorem 3 again. Noting that the determinant of n30 þ a2 n0 =2 þ a1 =4 is ða1 =8Þ2 þ ða2 =6Þ3 > 0, only one real root exists for gðn0 Þ and this root is given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 a1 g ¼ þ ða1 =2Þ2 þ ð2a2 =3Þ3 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 a1 þ ða1 =2Þ2 þ ð2a2 =3Þ3 ð19Þ 2 Therefore the coefficients a0 , a1 , a2 must comply with the conditions g < 0 and
fðgÞ 6 0
ð20Þ
4
L. Sun / Computers and Structures 80 (2002) 1–8
theory of algebraic equations of the fourth order, the roots of Eq. (18) are identical to the roots of the quadratic equation system pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a1 n20 þ 2y a2 n0 þ y pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 0 ð24Þ 2 2y a2 n20
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a1 2y a2 n0 þ y þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 0 2 2y a2
ð25Þ
where y is any real root of the cubic equation y 3 a2 y 2 =2 a0 y þ a2 a0 =2 a21 =8 ¼ 0
ð26Þ
To determine roots for Eqs. (24) and (25), one needs to identify the real root y of Eq. (26). Using transform x ¼ y a2 =6, Eq. (26) is converted to zðxÞ ¼ x3 þ p0 x þ cq0 ¼ 0
ð27Þ
where p0 ¼ a0 a22 =12 and
Fig. 1. The root distribution of Eq. (18).
for f to have real (in fact, negative) zeros and p to have imaginary zeros. The first condition of Eq. (20) is equivalent to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 a1 þ ða1 =2Þ2 þ ð2a2 =3Þ3 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 a1 þ ða1 =2Þ2 þ ð2a2 =3Þ3 < ð21Þ 2 As can easily be seen, inequality (21) satisfies all the time, therefore confirming that this first condition always holds. Then the polynomial f has two negative zeros if and only if the second condition fðgÞ < 0 is true. The case fðgÞ ¼ 0 determines a negative double zero, which corresponds to a pole of the second order. To examine if poles of the third and fourth order exist for the integral (16) in the case of imaginary root, one needs to inspect the second and third order derivatives of the function fðn0 Þ. According to Eq. (18) it is straightforward to see f00 ðn0 Þ ¼ 12n20 þ 2a2 ¼ 0
ð22Þ
f000 ðn0 Þ ¼ 24n0 ¼ 0
ð23Þ
Clearly, none of the solution of Eqs. (22) and (23) is valid in this case (i.e., n0 2 R; n0 6¼ 0). Therefore, no higher order (greater than second order) poles exist. The following is contributed to actually determine the two negative roots of Eq. (18). According to the
q0 ¼
a0 a22 a3 a2 2 1: 3 108 8
Since zða2 =3Þ ¼ a21 =8 < 0 and zð1Þ > 0, clearly, at least one real root x > a2 =3 satisfies Eq. (27). In other words, at least one real root y > a2 =2 satisfies Eq. (26). The determinant D0 of Eq. (27) is given by D0 ¼ ðq0 =2Þ2 þ ðp0 =3Þ3 . According to the theorem of cubic equation, three roots of Eq. (27) are given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi q pffiffiffiffiffiffiffi 3 3 ð28Þ x1 ¼ q0 =2 þ D0 þ q0 =2 D0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 3 3 x2 ¼ w q0 =2 þ D0 þ w2 q0 =2 D0
ð29Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 3 3 q0 =2 þ D0 þ w q0 =2 D0 ð30Þ pffiffiffi pffiffiffi where w ¼ ð1 þ i 3Þ=2 and w2 ¼ ð1 i 3Þ=2. Denote the largest real root of xi ði ¼ 1–3Þ as xmax and let y ¼ xmax þ a2 =6. Substituting this y into Eqs. (24) and (25) gives all the four roots that also satisfy Eq. (18). The determinants a2 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiof Eqs. (24) and (25) are D1 ¼ 2y pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2a1 = 2y a2 Þ and D2 ¼ 2y a2 ð2a1 = 2y a2 Þ, respectively. As discussed earlier in this section, under the condition fðgÞ < 0, where g is given by Eq. (19), two negative roots exist for Eq. (18). Noticing D1 > D2 one concludes that only two real roots of Eq. (24) can be valid for Eq. (18), i.e. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi 2ðxmax a2 =3Þ D1 n0 ¼ ð31Þ 2 x3 ¼ w2
4.2. The condition of existing complex roots for Eq. (17) Define the complex root of Eq. (17) as n ¼ na þ inb ðna ; nb 2 R; and na nb ¼ 6 0Þ. Substituting this complex
L. Sun / Computers and Structures 80 (2002) 1–8
root into Eq. (17) and comparing the real and imaginary parts gives the following quadratic equation systems ( ðn2a n2b Þ2 4n2a n2b þ K þ Cvnb mv2 ðn2a n2b Þ ¼ 0 4na nb ðn2a
n2b Þ
Cvna 2mv2 na nb ¼ 0
ð32Þ Since na ¼ 0 corresponding to the imaginary root has been discussed in previous section, it is required here that na nb 6¼ 0. So the second equation of Eq. (32) is reduced to n2a ¼ n2b þ ðCv þ 2mv2 nb Þ=4nb Replacing gives
n2a
ð33Þ
in the first equation of Eq. (32) by Eq. (33)
64n6b þ 32a2 n4b þ 4ða22 4a0 Þn2b a21 ¼ 0
ð34Þ
Eq. (34) is an algebraic equation of sixth order. By putting u ¼ n2b into Eq. (34) one obtains a cubic equation in terms of u zðuÞ ¼ u3 þ
a2 2 a22 4a0 a2 u þ u 1 ¼0 2 16 64
ð35Þ
By Theorem 3 this equation has one and only one positive root. Using transformation b ¼ u a2 =6 one can convert Eq. (35) into zðbÞ ¼ b3 þ sb þ q ¼ 0
ð36Þ
where s ¼ ða22 þ 12a0 Þ=48 and q ¼ a0 a2 =24 a21 =64 a32 =864. Since Eq. (35) has only one positive root, it implies that only one real root of Eq. (36) satisfying b > a2 =6. Similar to foregoing analysis in Section 4.1, three roots of Eq. (36) can be given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi q pffiffiffiffi 3 3 b1 ¼ q=2 þ D þ q=2 D ð37Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi pffiffiffiffi 3 3 b2 ¼ w q=2 þ D þ w2 q=2 D
ð38Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi pffiffiffiffi 3 3 q=2 þ D þ w q=2 D
ð39Þ
b3 ¼ w2
where D ¼ ðq=2Þ2 þ ðs=3Þ3 , w and w2 are the same as defined for Eqs. (28)–(30). Denote bmax as the largest real root of bi ði ¼ 1–3Þ. nb is then given by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi nb ¼ u ¼ bmax þ a2 =6 ð40Þ
5
duplicated root, it must satisfy the derivative equation of Eq. (17), i.e., 4n3 2mv2 n iCv ¼ 0
ð42Þ
This equation is obtained by taking the first derivative of the left-hand side of Eq. (17). Substituting n ¼ na þ inb into Eq. (42) and rearranging the real and imaginary parts gives
2 na 3n2b mv2 =2 ¼ 0 ð43Þ 3n2a nb n3b mv2 nb =2 Cv=4 ¼ 0 The first equation in Eq. (43) gives n2a in terms of nb . Substitution of n2a into the second equation of Eq. (43) leads to n3b þ mv2 nb =8 Cv=32 ¼ 0
ð44Þ
Because nb must also satisfy Eq. (33), combination of Eq. (33) and the first equation of Eq. (43) gives p ffiffiffiffiffiffi 3 nb ¼ Cv=2 ð45Þ Substituting this result into Eq. (44) results in nb ¼
3C 4mv
ð46Þ
Noticing that nb respectively given by Eqs. (45) and (46) are inconsistent in their signs, it is concluded that the initial assumption of duplicated roots is not correct and no high order poles exist if roots of characteristic Eq. (17) are complex. To synthesize, the procedure of seeking the poles of the integrand of Eq. (16) involves the following steps: (1) Compute g using Eq. (19) and check fðgÞ using Eq. (18); (2) If fðgÞ > 0, no imaginary and high order (i.e., P 2) poles exist; (3) If fðgÞ ¼ 0, there is a second order pole existing, which is given by ig; (4) If fðgÞ < 0, there are two distinct imaginary poles and they are given by in0 where n0 is computed from Eq. (31); (5) Four distinct complex poles exist, which take the form of n ¼ na þ inb and can be obtained from Eqs. (40) and (41).
na can be obtained by solving Eq. (33) na ¼ ½n2b þ ða1 þ 2a2 nb Þ=4nb 1=2
ð41Þ
In Section 4.1 it is showed that a pole of second order may exist when zeros of the characteristic equation are imaginary. Here it is also needed to check if high order poles exist when zeros of the characteristic equation are complex. Assume that n ¼ na þ inb ðna 6¼ 0 and nb 6¼ 0Þ is a duplicated complex root of Eq. (17). Since it is a
5. A closed-form solution In general, four poles distributed in the complex n-plane exist for the generalized integral (16). Two cases need to be discussed. One case is the response of the beam in front of the moving load (i.e., x vt P 0). The other case is the response of the beam behind the moving load (i.e., x vt < 0).
6
L. Sun / Computers and Structures 80 (2002) 1–8
First the solution corresponding to x vt P 0 is considered. In this case, only those poles above the real axis ReðnÞ have contributions to the integral (16). Since no second order poles are present in the upper half n-plane, applying the theorem of residue gives Z 1 P exp½inðx vtÞ dn yðx; tÞ ¼ 2pEI 1 n4 mv2 n2 þ K iCvn "
# X p exp½inðx vtÞ 2pi res 4 ¼ 2pEI n mv2 n2 þ K iCvn Im n>0 ð47Þ Since all the poles are of the first order, it follows that
exp½inðx vtÞ res 4 2 2 n mv n þ K iCvn Im n>0 ¼
exp½inðx vtÞ 4n3 2mv2 n iCv
ð48Þ
Combining Eqs. (47) and (48) gives yðx; tÞ ¼
iP X exp½inðx vtÞ EI Im n>0 4n3 2mv2 n iCv
ð49Þ
In the case of x vt < 0, the same method applies. The result also remains similar form as given by Eq. (49) except that Im n < 0 and the residue of a second order pole should be included since this second order pole is present in the lower half n-plane. Let nl be the pole of the second order and denote n1 and n2 as the other two poles of the first order. The residue of this second order pole can be evaluated as
exp½inðx vtÞ res 4 2 2 n mv n þ K iCvn Im nl <0
d exp½inðx vtÞ ðn nl Þ2 4 ¼ lim n!nl dn n mv2 n2 þ K iCvn
d exp½inðx vtÞ ¼ lim n!nl dn ðn n1 Þðn n2 Þ " # iðx vtÞ ðnl n1 Þ þ ðnl n2 Þ ¼ ðnl n1 Þðnl n2 Þ ðnl n1 Þ2 ðnl n2 Þ2
exp½inðx vtÞ
ð50Þ
So far, all the solutions corresponding to different cases have been obtained. In the condition of beams resting on an elastic subgrade, it has been demonstrated in the previous work [1,20] that the wave excited by the moving load is propagating along the beam in two directions and it looks symmetric if the observer is traveling at the same speed with the moving load. These situations correspond to the solutions of x vt > 0 and x vt < 0, respectively. Mathematically, this requires the poles of the integrand are symmetrically distributed with respect to the real axis. For the elastic subgrade this symmetry always holds, but for the viscoelastic subgrade
it is not the case any more. The asymmetric poles result in distinct waves along the positive and negative longitudinal directions. 6. Maximal response approximation It is of interest to examine the extreme value of the dynamic response under a moving load. It has been demonstrated that the maximal response of beam on a viscoelastic subgrade occurs behind the moving load [11]. However, a good approximation of this maximal response will be the response of the position exactly beneath the moving load. Clearly, the approximate maximal response should be slightly less than the virtual maximum. To get the estimated maximal response, substituting x vt ¼ 0 into Eq. (49) gives iP X 1 yðx; tÞ ¼ ð51Þ EI Im n >0 4n3l 2mv2 nl iCv l Substituting nl ¼ na þ inb into Eq. (51) and canceling the imaginary terms yield P X Bl yðx; tÞ ¼ ð52Þ EI Im n >0 A2l þ B2l l
where Al ¼ 4n3a 12na n2b 2mv2 na and Bl ¼ 12n2a nb 4n3b 2mv2 nb Cv. The physical model described herein can be thought of as the simplified mathematical representation of a vehicle traveling along the railway track or highway pavement. Usually, the vehicle speed is far less than the critical speed at which elastic waves propagates in the beam. In practice, engineers are often interested in knowing the dynamic coefficient, a concept that reflects the dynamic effect with respect to the static response. Formula (52) can be useful to define an approximate dynamic coefficient. 7. Numerical computations Based on the analytical results obtained in previous sections, numerical computations are performed to illustrate the dynamic deflection response of beam to a moving load. Parameters used in calculation are EI ¼ 2:3 kN m2 , K ¼ 68:9 Mpa, m ¼ 48:2 kg/m, P ¼ 10:5 kN. Fig. 2 plots the deflection of beam at position x ¼ 0 versus time for damping coefficient C ¼ 230 kN s/m2 and different load velocities. It should be pointed out that in the following figures positive deflection means the bottom of the beam suffering with the tension stress. It is clear from Fig. 2 that load velocities have significant effect on the shape of the dynamic deflection and the maximum deflection. Higher speeds lead to sharper shapes and larger maximum deflection. Also, it can be seen that the dynamic response is asymmetric with
L. Sun / Computers and Structures 80 (2002) 1–8
7
8. Concluding remarks
Fig. 2. Dynamic deflection of the beam at x ¼ 0 versus time (C ¼ 230 kN s/m2 ).
respect to time t ¼ 0 and the maximum deflection does not occur beneath the moving load. Rather, the maximum deflection occurs at time a little bit later than time t ¼ 0, implying that maximum deflection response at x ¼ 0 appears after the load has passed through this point. This phenomenon is the consequence of the action of damping and the reaction of the beam to the moving load is delayed by damping effects. This hysteretic phenomenon can be further observed in Fig. 3 where the load velocity is fixed at v ¼ 30 m/s while the damping coefficient varies from C ¼ 0 kN s/m2 to C ¼ 2300 kN s/m2 . In the case of no damping, the deflection response is symmetric with respect to time t ¼ 0. However, this is no longer true if damping effect is taken into account. This conclusion is consistent with earlier theoretical analysis [11]. It also can be recognized that before the moving load reach the point x ¼ 0, a negative deflection is generated at this position of the beam. The magnitude of the negative deflection depends on damping coefficient, load velocity as well as other parameters. The maximum deflection of the beam decreases as the damping increase. For instance, the maximum deflection of beam with no damping is 0.725 mm, while the maximum deflection of beam with damping C ¼ 1150 kN s/m2 is 0.186 mm, only 25.7% of the former.
Fig. 3. Dynamic deflection of the beam at x ¼ 0 versus time (v ¼ 30 m/s).
In this paper, Fourier transform is used to solve the problem of steady-state response of a beam on a viscoelastic subgrade subjected to a moving constant point load. The solution is constructed in the form of the convolution of the Green’s function of the beam. The theorem of residue is employed to evaluate the generalized integral such that a closed-form solution can be achieved. All the different combinations of damping and speed are discussed and analytical solutions are presented. Numerical computations show that deflection response of beam with damping is asymmetric and the maximum deflection of a point of the beam occurs behind the moving load. Also, it is observed that damping and load velocity have similar effects on the deflection response of the beam.
Acknowledgements The author wants to thank the anonymous reviewers for their insights and valuable comments on the original manuscript, which enhance the contents and presentation of this paper. The author is also very grateful to one of the reviewers, whose suggestions lead to a concise description of roots of the algebraic equation. References [1] Achenbach JD, Sun C. Dynamic response of beam on viscoelastic subgrade. J Eng Mech, ASCE 1965;91:61. [2] Choros J, Adams GG. A steadily moving load on an elastic beam resting on a tensionless Winkler subgrade. J Appl Mech, ASME 1979;46(1):175–80. [3] Conte SD, de Boor C. Elementary numerical analysis: an algorithmic approach. 3rd ed. New York: McGraw-Hill; 1980. [4] Criner HE, McCann GD. Rails on elastic subgrades under the influence of high-speed traveling loads. J Appl Mech, ASME 1953;20(1):13. [5] Elattary MA. Moving loads on an infinite plate strip of constant thickness. J Phys D: Appl Phys 1991;24(4):541–6. [6] Florence AL. Traveling force on a Timoshenko beam. J Appl Mech, ASME 1965;32(2):351–8. [7] Fryba L. Vibration of solids and structures under moving loads. Groningen, Netherlands: Noordhoff International; 1977. [8] Hardy MSA, Cebon D. Importance of speed and frequency in flexible pavement response. J Eng Mech, ASCE 1994;120(3):463–82. [9] Huang CC. Traveling loads on a viscoelastic Timoshenko beam. J Appl Mech, ASME 1977;44(1):183–4. [10] Jezequel L. Response of periodic systems to a moving load. J Appl Mech, ASME 1981;48(3):613–8. [11] Kenney JT. Steady state vibrations of beam on elastic subgrade for moving loads. J Appl Mech, ASME 1954; 21(4):359–64.
8
L. Sun / Computers and Structures 80 (2002) 1–8
[12] Lee HP. Dynamic response of a beam with intermediate point constrains subject to a moving load. J Sound Vibr 1994;171(3):361–8. [13] Morse PM, Feshbach H. Methods of theoretical physics: Part I and II. New York: McGraw-Hill Book Company Inc; 1953. [14] Mulcahy TM. Steady-state response of a beam on idealized strain-hardening subgrades for a moving load. J Appl Mech, ASME 1973;42(4):1040–4. [15] Pan G, Atluri SN. Dynamic response of finite sized elastic runways subjected to moving loads: a coupled BEM/FEM approach. Int J Num Meth Eng 1995;38(18): 3143–66. [16] Patil S. Response of infinite railroad track to vibrating mass. J Eng Mech, ASCE 1988;114(4):688–703. [17] Saito H, Terasawa T. Steady-state vibrations of a beam on a Pasternak subgrade for moving loads. J Appl Mech, ASME 1980;47(4):879–83. [18] Steele CR. The finite beam with a moving load. J Appl Mech, ASME 1967;34(1):111.
[19] Steele CR. The Timoshenko beam with a moving load. J Appl Mech, ASME 1968;35(3):481–8. [20] Sun L. Theoretical investigations on vehicle-ground dynamic interaction. Final Report prepared for National Science Foundation of China, Southeast University, Nanjing, China. [21] Sun L, Deng X. Transient response for infinite plate on Winkler subgrade by a moving distributed load. Chin J Appl Mech 1997;14(2):72–8. [22] Sun L, Deng X. Dynamic analysis of infinite beam under the excitation of moving line loads. Appl Math Mech 1998;19(4):367–73. [23] Sun L, Greenberg B. Dynamic response of linear systems to moving stochastic sources. J Sound Vibr 2000;229(4): 957–72. [24] Yoder EJ, Witczak MW. Principles of pavement design. New York: Wiely; 1975. [25] Adams GG, Bogy DB. Steady solutions for moving loads on elastic beams with one-sided constraints. J Appl Mech 1975;42(4):800–5.
A closed-form solution of beam on viscoelastic subgrade subjected to moving loads Lu Sun * ECJ Hall 6.10, Department of Civil Engineering, The University of Texas at Austin, Austin, TX 78712, USA Received 2 June 2000; accepted 30 September 2001
Abstract A closed-form solution of beam on viscoelastic subgrade subjected to moving loads is proposed in this paper. The Green’s function of the beam is obtained by means of Fourier transform. Theory of linear partial differential equation is used to express the deflection of the beam in terms of the Green’s function. Beam’s deflection is represented as a generalized integral using inverse Fourier transform. To evaluate this generalized integral analytically, poles of the integrand are identified using the theory of algebraic equation. Every pole and its order are isolated and given in a closed form. The theorem of residue is then applied to represent the generalized integral in the form of contour integral in the complex plane. Closed-form deflection and numerical computation of distinct cases are provided as an illustration of the proposed solution. Ó 2002 Published by Elsevier Science Ltd. Keywords: Beam; Moving loads; Viscoelastic subgrade; Green’s function; Fourier transform; Complex function; Residue theorem
1. Introduction The response of a beam to moving loads has been studied extensively over the past several decades [7]. Considerable studies have been conducted to investigate the deflection response of beams [1,2,6,8,11,17–19]. A growing interest on this topic is recently found in railway industry and highway industry due to the reason that this kind of physical model can be used as an idealized representation of dynamic interactions of railtrack/train and pavement/vehicle [9,10,14,16,21,22, 24,25]. An interesting manifestation of this subject is that when the subgrade is modeled as an elastic subgrade, a critical load velocity is found existing, which may cause a significant variation of a beam’s response [11]. Waves generated by a moving load with a supersonic velocity excite different patterns of response of the beam as they do when the load velocity is subsonic. This critical velocity is also predicted and observed by other researchers
by means of either electrical analogy using electrical device and experimental measurement [4] or theoretical analysis [16,20]. Since high-speed vehicles are getting more widely adopted as surface transportation carriers, considerable attention has been paid to the response of infrastructure under high-speed loading condition since 1990s [5,12,15,20,23]. In most previous studies the subgrade of a beam is usually modeled as an elastic Winkler foundation. Moving at the critical speed, the load will excite the beam to form an infinite deflection [11,20]. This idealization of subgrade is, however, not quite true because any physical system in reality possesses damping, which leads to energy dissipation. It can be expected that the unrealistic infinite deflection will disappear if the damping effect is included in the mathematical model of a beam subjected to moving loads. In this paper the effect of damping on the deflection response is taken into account. A closedform representation of beam deflection is presented. 2. The Green’s function
*
Tel.: +1-512-476-1906. E-mail address: [email protected] (L. Sun).
The governing equation of a flexible beam on a viscoelastic subgrade can be given by [22]
0045-7949/02/$ - see front matter Ó 2002 Published by Elsevier Science Ltd. PII: S 0 0 4 5 - 7 9 4 9 ( 0 1 ) 0 0 1 6 2 - 6
2
L. Sun / Computers and Structures 80 (2002) 1–8
EI
o4 y oy o2 y þ Ky þ C þ m 2 ¼ F ðx; tÞ 4 ox ot ot
ð1Þ
where yðx; tÞ represents the deflection of the beam, x represents the traveling direction of the moving load, and t represents time. Also, EI is the rigidity of the beam, E is Young’s modulus of elasticity, I is the moment of inertia of the beam, m is the unit mass of the beam, and F ðx; tÞ is external loads. A moving concentrated load can be given by F ðx; tÞ ¼ P dðx vtÞ
ð2Þ
where P is the amplitude of the applied load, and dðÞ is the Dirac-delta function, which is defined by Z 1 dðx x0 Þf ðxÞ dx ¼ f ðx0 Þ ð3Þ 1
Assume that the beam is at rest initially (i.e., initial deflection, velocity, and acceleration of the beam is zero). Eqs. (2) and (3) constitute the mathematical description of the problem of interest. According to the theory of linear partial differential equation [13,23], the Green’s function represents the fundamental solution of the equation in which the load is given by the Dirac-delta function. More specifically, the Green’s function of this problem is defined as the solution of Eq. (1) given that the applied load is characterized by Fd ðx; tÞ ¼ dðx x0 Þdðt t0 Þ
F½f ðnÞ ðtÞ ¼ ðixÞn F½f ðtÞ
Since the right-hand side of Eq. (7) Fed ðn; xÞ is the representation of Fd ðx; tÞ in the frequency domain, it is necessary to evaluate Fd ðx; tÞ in the frequency domain too. This can be achieved by taking two-dimensional Fourier transform on both sides of Eq. (4) Z 1 Z 1 Fed ðn; xÞ ¼ dðx x0 Þdðt t0 Þ 1
1
ð5Þ
f ðx; tÞ ¼ F1 ½f~ðn; xÞ Z 1 Z ¼ ð2pÞ2 1
1
e ðn; x; x0 ; t0 Þ ¼ exp½iðnx0 þ xt0 Þ ðEIn4 þ K G þ iCx mx2 Þ1
ð10Þ
The Green’s function given by Eq. (10) is represented in the frequency domain and must be converted to the time domain. To this end, taking the inverse Fourier transform on both sides of Eq. (10) gives Z 1 Z 1 Gðx; t; x0 ; t0 Þ ¼ ð2pÞ2 1
1
expfi½nðx x0 Þ þ xðt t0 Þ g
dn dx EIn4 þ K þ iCx mx2 ð11Þ Formula (11) is the Green’s function of the beam resting on a viscoelastic subgrade. The Green’s function serves as a fundamental solution of a partial differential equation. It is very useful when dealing with linear systems.
3. Integral representation of the solution f~ðn; xÞ exp½iðnx þ xtÞ dn dx
1
ð6Þ where F½ and F1 ½ are the Fourier transform and its inversion, respectively. To solve the Green’s function, taking two-dimensional Fourier transform on both sides of Eq. (1) gives 4e
ð9Þ
Here, the property of the Dirac-delta function, i.e., Eq. (3), is utilized to evaluate the integration in Eq. (9). Substituting the result of Eq. (9) into Eq. (7) and realizing that Eq. (7) is an algebraic equation, it is straightforward to see
ð4Þ
1
1
exp½iðnx þ xtÞ dx dt ¼ exp½iðnx0 þ xt0 Þ
Define two-dimensional Fourier transform and its inversion as f~ðn; xÞ ¼ F½f ðx; tÞ Z 1 Z 1 ¼ f ðx; tÞ exp½iðnx þ xtÞ dx dt
ð8Þ
e þ iCx G e mx2 G e ¼ Fed ðn; xÞ EIn G þ K G
ð7Þ
in which Fed ðn; xÞ is the Fourier transform of Fd ðx; tÞ. In Eq. (7) the deflection response yðx; tÞ has been replaced e ¼G e ðn; x; x0 ; t0 Þ to indicate the Fourier by the symbol G transform of the Green’s function. Also, the following property of Fourier transform is used in the derivation of Eq. (7)
According to Sun and Greenberg [23], the solution of Eq. (2) can be constructed by integrating the Green’s function in all dimensions Z t Z 1 F ðx0 ; t0 ÞGðx; t; x0 ; t0 Þ dx0 dt0 ð12Þ yðx; tÞ ¼ 1
1
Taking Eqs. (4) and (11) into Eq. (12) gives Z t Z 1 Z 1 Z 1 yðx; tÞ ¼ 1
1
1
1
P dðx0 vt0 Þ expfi½nðx x0 Þ þ xðt t0 Þ g
ð2pÞ2 EIðn4 þ K þ iCx mx2 Þ
dn dx dx0 dt0 where K ¼ K=EI and m ¼ m=EI. Note that
ð13Þ
L. Sun / Computers and Structures 80 (2002) 1–8
Z
1
dðx0 vt0 Þ expðinx0 Þ dx0 ¼ expðinvt0 Þ
ð14Þ
exp½iðx þ vnÞt0 dt0 ¼ 2pdðx þ vnÞ
ð15Þ
1
Z
t 1
Substituting Eqs. (14) and (15) into Eq. (13) and reapplying the property (3) one gets Z 1 P exp½inðx vtÞ dn ð16Þ yðx; tÞ ¼ 2pEI 1 n4 mv2 n2 þ K iCvn So far the integral representation of the dynamic deflection of the beam under a moving concentrated load has been obtained. Expression (16) can be further developed using complex analysis. In the following section, the theorem of residue is employed to evaluate the integration (16).
Before the integration (16) is further evaluated, it is necessary to determine the roots of the characteristic equation ð17Þ
Assume none of these parameters m, K, C and v being zero. Essentially, one seeks zeros of the polynomial at the left-hand side of the characteristic equation (17). Because several algebraic theorems [3] are used in the derivation of roots of Eq. (17), they are cited and stated below. Theorem 1 (fundamental theorem of algebra).P If pðnÞ is a polynomial of degree n P 1, that is, pðnÞ ¼ ni¼0 ai ni , with coefficients ai ði ¼ 1; 2; . . . ; nÞ real or complex numbers and an 6¼ 0, then pðnÞ has at least one zero; i.e., there exists a complex number n such that pðnÞ ¼ 0. Theorem 2. A polynomial of degree n P 1 has at most n zeros, counting multiplicity. Theorem 3 (Descartes’ rule of signs). The number np of positive real zeros of a polynomial pðnÞ is less than or equal to the number of variations nv in sign of the coefficients of pðnÞ. The number of negative real zeros of pðnÞ is at most equal to the number of sign changes in the coefficients of the polynomial pðnÞ. Moreover, the difference nv np is an even integer. Theorem 4. Every zero of a polynomial of degree n P 1 lies in the circular region defined by jnj 6 q, where the radius q is given by q ¼ 1 þ max jai =an j. 0 6 i 6 n1
Define abbreviations a0 ¼ K, a1 ¼ Cv, a2 ¼ mv2 . Due to Theorem 4 all zeros of the left-hand side of Eq. (17) (i.e., pðnÞ ¼ n4 a2 n2 ia1 n þ a0 ) lie in the open disk centered at the origin of the complex plan and whose radius q is given by q ¼ 1 þ max0 6 i 6 2 jai j. This result may be instrumental in locating the zeros when the numerical methods are used to compute the roots of Eq. (17). Since a coefficient of Eq. (17) is imaginary, according to Theorems 1 and 2, no real-valued roots exist for Eq. (17) and all four roots must be complex. In the following only imaginary roots and complex roots are considered analytically. Imaginary roots are discussed separately because in this case duplicated roots may exist, which is very crucial in determining the order of the poles in the generalized integral (16). 4.1. The condition of existing imaginary roots for Eq. (17) To derive the condition of exiting imaginary roots of Eq. (17), define the imaginary root n ¼ in0 ðn0 6¼ 0 and n0 2 RÞ. Substitution of the imaginary root into Eq. (17) gives
4. Nature of zeros of the characteristic equation
n4 mv2 n2 þ K iCvn ¼ 0
3
fðn0 Þ ¼ n40 þ a2 n20 þ a1 n0 þ a0 ¼ 0
ð18Þ
The objective here is to seek real roots of Eq. (18) and the condition under which these real roots exist. According to Theorem 3 one may conclude that two possibilities arise: (1) the polynomial f has two real negative zeros (could be duplicated) and one complex conjugate pair or (2) two complex conjugate pairs. Define polynomials f ðn0 Þ ¼ n40 and fj ðn0 Þ ¼ a2 n20 a1 n0 a0 ¼ 2 a2 ðn0 þ ðC=2mvÞÞ2 þ ðC =4mÞ K, ðj ¼ 1–3Þ with respect to n0 . These polynomials are plotted in Fig. 1. The first possibility corresponds to polynomials f1 and f2 , while the second possibility corresponds to polynomial f3 . Since n0 must be real, only the first possibility is considered. In order for (1) to be the case it is necessary and sufficient that the graph of f crosses the n0 axes twice and has a negative minimum at a point g < 0. Then, the polynomial gðn0 Þ ¼ fðn0 Þ0 ¼ 4n30 þ 2a2 n0 þ a1 must have a negative zero g, which is the case due to Theorem 3 again. Noting that the determinant of n30 þ a2 n0 =2 þ a1 =4 is ða1 =8Þ2 þ ða2 =6Þ3 > 0, only one real root exists for gðn0 Þ and this root is given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 a1 g ¼ þ ða1 =2Þ2 þ ð2a2 =3Þ3 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 a1 þ ða1 =2Þ2 þ ð2a2 =3Þ3 ð19Þ 2 Therefore the coefficients a0 , a1 , a2 must comply with the conditions g < 0 and
fðgÞ 6 0
ð20Þ
4
L. Sun / Computers and Structures 80 (2002) 1–8
theory of algebraic equations of the fourth order, the roots of Eq. (18) are identical to the roots of the quadratic equation system pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a1 n20 þ 2y a2 n0 þ y pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 0 ð24Þ 2 2y a2 n20
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a1 2y a2 n0 þ y þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 0 2 2y a2
ð25Þ
where y is any real root of the cubic equation y 3 a2 y 2 =2 a0 y þ a2 a0 =2 a21 =8 ¼ 0
ð26Þ
To determine roots for Eqs. (24) and (25), one needs to identify the real root y of Eq. (26). Using transform x ¼ y a2 =6, Eq. (26) is converted to zðxÞ ¼ x3 þ p0 x þ cq0 ¼ 0
ð27Þ
where p0 ¼ a0 a22 =12 and
Fig. 1. The root distribution of Eq. (18).
for f to have real (in fact, negative) zeros and p to have imaginary zeros. The first condition of Eq. (20) is equivalent to qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 a1 þ ða1 =2Þ2 þ ð2a2 =3Þ3 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1=3 a1 þ ða1 =2Þ2 þ ð2a2 =3Þ3 < ð21Þ 2 As can easily be seen, inequality (21) satisfies all the time, therefore confirming that this first condition always holds. Then the polynomial f has two negative zeros if and only if the second condition fðgÞ < 0 is true. The case fðgÞ ¼ 0 determines a negative double zero, which corresponds to a pole of the second order. To examine if poles of the third and fourth order exist for the integral (16) in the case of imaginary root, one needs to inspect the second and third order derivatives of the function fðn0 Þ. According to Eq. (18) it is straightforward to see f00 ðn0 Þ ¼ 12n20 þ 2a2 ¼ 0
ð22Þ
f000 ðn0 Þ ¼ 24n0 ¼ 0
ð23Þ
Clearly, none of the solution of Eqs. (22) and (23) is valid in this case (i.e., n0 2 R; n0 6¼ 0). Therefore, no higher order (greater than second order) poles exist. The following is contributed to actually determine the two negative roots of Eq. (18). According to the
q0 ¼
a0 a22 a3 a2 2 1: 3 108 8
Since zða2 =3Þ ¼ a21 =8 < 0 and zð1Þ > 0, clearly, at least one real root x > a2 =3 satisfies Eq. (27). In other words, at least one real root y > a2 =2 satisfies Eq. (26). The determinant D0 of Eq. (27) is given by D0 ¼ ðq0 =2Þ2 þ ðp0 =3Þ3 . According to the theorem of cubic equation, three roots of Eq. (27) are given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi q pffiffiffiffiffiffiffi 3 3 ð28Þ x1 ¼ q0 =2 þ D0 þ q0 =2 D0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 3 3 x2 ¼ w q0 =2 þ D0 þ w2 q0 =2 D0
ð29Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi 3 3 q0 =2 þ D0 þ w q0 =2 D0 ð30Þ pffiffiffi pffiffiffi where w ¼ ð1 þ i 3Þ=2 and w2 ¼ ð1 i 3Þ=2. Denote the largest real root of xi ði ¼ 1–3Þ as xmax and let y ¼ xmax þ a2 =6. Substituting this y into Eqs. (24) and (25) gives all the four roots that also satisfy Eq. (18). The determinants a2 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiof Eqs. (24) and (25) are D1 ¼ 2y pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2a1 = 2y a2 Þ and D2 ¼ 2y a2 ð2a1 = 2y a2 Þ, respectively. As discussed earlier in this section, under the condition fðgÞ < 0, where g is given by Eq. (19), two negative roots exist for Eq. (18). Noticing D1 > D2 one concludes that only two real roots of Eq. (24) can be valid for Eq. (18), i.e. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi 2ðxmax a2 =3Þ D1 n0 ¼ ð31Þ 2 x3 ¼ w2
4.2. The condition of existing complex roots for Eq. (17) Define the complex root of Eq. (17) as n ¼ na þ inb ðna ; nb 2 R; and na nb ¼ 6 0Þ. Substituting this complex
L. Sun / Computers and Structures 80 (2002) 1–8
root into Eq. (17) and comparing the real and imaginary parts gives the following quadratic equation systems ( ðn2a n2b Þ2 4n2a n2b þ K þ Cvnb mv2 ðn2a n2b Þ ¼ 0 4na nb ðn2a
n2b Þ
Cvna 2mv2 na nb ¼ 0
ð32Þ Since na ¼ 0 corresponding to the imaginary root has been discussed in previous section, it is required here that na nb 6¼ 0. So the second equation of Eq. (32) is reduced to n2a ¼ n2b þ ðCv þ 2mv2 nb Þ=4nb Replacing gives
n2a
ð33Þ
in the first equation of Eq. (32) by Eq. (33)
64n6b þ 32a2 n4b þ 4ða22 4a0 Þn2b a21 ¼ 0
ð34Þ
Eq. (34) is an algebraic equation of sixth order. By putting u ¼ n2b into Eq. (34) one obtains a cubic equation in terms of u zðuÞ ¼ u3 þ
a2 2 a22 4a0 a2 u þ u 1 ¼0 2 16 64
ð35Þ
By Theorem 3 this equation has one and only one positive root. Using transformation b ¼ u a2 =6 one can convert Eq. (35) into zðbÞ ¼ b3 þ sb þ q ¼ 0
ð36Þ
where s ¼ ða22 þ 12a0 Þ=48 and q ¼ a0 a2 =24 a21 =64 a32 =864. Since Eq. (35) has only one positive root, it implies that only one real root of Eq. (36) satisfying b > a2 =6. Similar to foregoing analysis in Section 4.1, three roots of Eq. (36) can be given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi q pffiffiffiffi 3 3 b1 ¼ q=2 þ D þ q=2 D ð37Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi pffiffiffiffi 3 3 b2 ¼ w q=2 þ D þ w2 q=2 D
ð38Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffi pffiffiffiffi 3 3 q=2 þ D þ w q=2 D
ð39Þ
b3 ¼ w2
where D ¼ ðq=2Þ2 þ ðs=3Þ3 , w and w2 are the same as defined for Eqs. (28)–(30). Denote bmax as the largest real root of bi ði ¼ 1–3Þ. nb is then given by pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi nb ¼ u ¼ bmax þ a2 =6 ð40Þ
5
duplicated root, it must satisfy the derivative equation of Eq. (17), i.e., 4n3 2mv2 n iCv ¼ 0
ð42Þ
This equation is obtained by taking the first derivative of the left-hand side of Eq. (17). Substituting n ¼ na þ inb into Eq. (42) and rearranging the real and imaginary parts gives
2 na 3n2b mv2 =2 ¼ 0 ð43Þ 3n2a nb n3b mv2 nb =2 Cv=4 ¼ 0 The first equation in Eq. (43) gives n2a in terms of nb . Substitution of n2a into the second equation of Eq. (43) leads to n3b þ mv2 nb =8 Cv=32 ¼ 0
ð44Þ
Because nb must also satisfy Eq. (33), combination of Eq. (33) and the first equation of Eq. (43) gives p ffiffiffiffiffiffi 3 nb ¼ Cv=2 ð45Þ Substituting this result into Eq. (44) results in nb ¼
3C 4mv
ð46Þ
Noticing that nb respectively given by Eqs. (45) and (46) are inconsistent in their signs, it is concluded that the initial assumption of duplicated roots is not correct and no high order poles exist if roots of characteristic Eq. (17) are complex. To synthesize, the procedure of seeking the poles of the integrand of Eq. (16) involves the following steps: (1) Compute g using Eq. (19) and check fðgÞ using Eq. (18); (2) If fðgÞ > 0, no imaginary and high order (i.e., P 2) poles exist; (3) If fðgÞ ¼ 0, there is a second order pole existing, which is given by ig; (4) If fðgÞ < 0, there are two distinct imaginary poles and they are given by in0 where n0 is computed from Eq. (31); (5) Four distinct complex poles exist, which take the form of n ¼ na þ inb and can be obtained from Eqs. (40) and (41).
na can be obtained by solving Eq. (33) na ¼ ½n2b þ ða1 þ 2a2 nb Þ=4nb 1=2
ð41Þ
In Section 4.1 it is showed that a pole of second order may exist when zeros of the characteristic equation are imaginary. Here it is also needed to check if high order poles exist when zeros of the characteristic equation are complex. Assume that n ¼ na þ inb ðna 6¼ 0 and nb 6¼ 0Þ is a duplicated complex root of Eq. (17). Since it is a
5. A closed-form solution In general, four poles distributed in the complex n-plane exist for the generalized integral (16). Two cases need to be discussed. One case is the response of the beam in front of the moving load (i.e., x vt P 0). The other case is the response of the beam behind the moving load (i.e., x vt < 0).
6
L. Sun / Computers and Structures 80 (2002) 1–8
First the solution corresponding to x vt P 0 is considered. In this case, only those poles above the real axis ReðnÞ have contributions to the integral (16). Since no second order poles are present in the upper half n-plane, applying the theorem of residue gives Z 1 P exp½inðx vtÞ dn yðx; tÞ ¼ 2pEI 1 n4 mv2 n2 þ K iCvn "
# X p exp½inðx vtÞ 2pi res 4 ¼ 2pEI n mv2 n2 þ K iCvn Im n>0 ð47Þ Since all the poles are of the first order, it follows that
exp½inðx vtÞ res 4 2 2 n mv n þ K iCvn Im n>0 ¼
exp½inðx vtÞ 4n3 2mv2 n iCv
ð48Þ
Combining Eqs. (47) and (48) gives yðx; tÞ ¼
iP X exp½inðx vtÞ EI Im n>0 4n3 2mv2 n iCv
ð49Þ
In the case of x vt < 0, the same method applies. The result also remains similar form as given by Eq. (49) except that Im n < 0 and the residue of a second order pole should be included since this second order pole is present in the lower half n-plane. Let nl be the pole of the second order and denote n1 and n2 as the other two poles of the first order. The residue of this second order pole can be evaluated as
exp½inðx vtÞ res 4 2 2 n mv n þ K iCvn Im nl <0
d exp½inðx vtÞ ðn nl Þ2 4 ¼ lim n!nl dn n mv2 n2 þ K iCvn
d exp½inðx vtÞ ¼ lim n!nl dn ðn n1 Þðn n2 Þ " # iðx vtÞ ðnl n1 Þ þ ðnl n2 Þ ¼ ðnl n1 Þðnl n2 Þ ðnl n1 Þ2 ðnl n2 Þ2
exp½inðx vtÞ
ð50Þ
So far, all the solutions corresponding to different cases have been obtained. In the condition of beams resting on an elastic subgrade, it has been demonstrated in the previous work [1,20] that the wave excited by the moving load is propagating along the beam in two directions and it looks symmetric if the observer is traveling at the same speed with the moving load. These situations correspond to the solutions of x vt > 0 and x vt < 0, respectively. Mathematically, this requires the poles of the integrand are symmetrically distributed with respect to the real axis. For the elastic subgrade this symmetry always holds, but for the viscoelastic subgrade
it is not the case any more. The asymmetric poles result in distinct waves along the positive and negative longitudinal directions. 6. Maximal response approximation It is of interest to examine the extreme value of the dynamic response under a moving load. It has been demonstrated that the maximal response of beam on a viscoelastic subgrade occurs behind the moving load [11]. However, a good approximation of this maximal response will be the response of the position exactly beneath the moving load. Clearly, the approximate maximal response should be slightly less than the virtual maximum. To get the estimated maximal response, substituting x vt ¼ 0 into Eq. (49) gives iP X 1 yðx; tÞ ¼ ð51Þ EI Im n >0 4n3l 2mv2 nl iCv l Substituting nl ¼ na þ inb into Eq. (51) and canceling the imaginary terms yield P X Bl yðx; tÞ ¼ ð52Þ EI Im n >0 A2l þ B2l l
where Al ¼ 4n3a 12na n2b 2mv2 na and Bl ¼ 12n2a nb 4n3b 2mv2 nb Cv. The physical model described herein can be thought of as the simplified mathematical representation of a vehicle traveling along the railway track or highway pavement. Usually, the vehicle speed is far less than the critical speed at which elastic waves propagates in the beam. In practice, engineers are often interested in knowing the dynamic coefficient, a concept that reflects the dynamic effect with respect to the static response. Formula (52) can be useful to define an approximate dynamic coefficient. 7. Numerical computations Based on the analytical results obtained in previous sections, numerical computations are performed to illustrate the dynamic deflection response of beam to a moving load. Parameters used in calculation are EI ¼ 2:3 kN m2 , K ¼ 68:9 Mpa, m ¼ 48:2 kg/m, P ¼ 10:5 kN. Fig. 2 plots the deflection of beam at position x ¼ 0 versus time for damping coefficient C ¼ 230 kN s/m2 and different load velocities. It should be pointed out that in the following figures positive deflection means the bottom of the beam suffering with the tension stress. It is clear from Fig. 2 that load velocities have significant effect on the shape of the dynamic deflection and the maximum deflection. Higher speeds lead to sharper shapes and larger maximum deflection. Also, it can be seen that the dynamic response is asymmetric with
L. Sun / Computers and Structures 80 (2002) 1–8
7
8. Concluding remarks
Fig. 2. Dynamic deflection of the beam at x ¼ 0 versus time (C ¼ 230 kN s/m2 ).
respect to time t ¼ 0 and the maximum deflection does not occur beneath the moving load. Rather, the maximum deflection occurs at time a little bit later than time t ¼ 0, implying that maximum deflection response at x ¼ 0 appears after the load has passed through this point. This phenomenon is the consequence of the action of damping and the reaction of the beam to the moving load is delayed by damping effects. This hysteretic phenomenon can be further observed in Fig. 3 where the load velocity is fixed at v ¼ 30 m/s while the damping coefficient varies from C ¼ 0 kN s/m2 to C ¼ 2300 kN s/m2 . In the case of no damping, the deflection response is symmetric with respect to time t ¼ 0. However, this is no longer true if damping effect is taken into account. This conclusion is consistent with earlier theoretical analysis [11]. It also can be recognized that before the moving load reach the point x ¼ 0, a negative deflection is generated at this position of the beam. The magnitude of the negative deflection depends on damping coefficient, load velocity as well as other parameters. The maximum deflection of the beam decreases as the damping increase. For instance, the maximum deflection of beam with no damping is 0.725 mm, while the maximum deflection of beam with damping C ¼ 1150 kN s/m2 is 0.186 mm, only 25.7% of the former.
Fig. 3. Dynamic deflection of the beam at x ¼ 0 versus time (v ¼ 30 m/s).
In this paper, Fourier transform is used to solve the problem of steady-state response of a beam on a viscoelastic subgrade subjected to a moving constant point load. The solution is constructed in the form of the convolution of the Green’s function of the beam. The theorem of residue is employed to evaluate the generalized integral such that a closed-form solution can be achieved. All the different combinations of damping and speed are discussed and analytical solutions are presented. Numerical computations show that deflection response of beam with damping is asymmetric and the maximum deflection of a point of the beam occurs behind the moving load. Also, it is observed that damping and load velocity have similar effects on the deflection response of the beam.
Acknowledgements The author wants to thank the anonymous reviewers for their insights and valuable comments on the original manuscript, which enhance the contents and presentation of this paper. The author is also very grateful to one of the reviewers, whose suggestions lead to a concise description of roots of the algebraic equation. References [1] Achenbach JD, Sun C. Dynamic response of beam on viscoelastic subgrade. J Eng Mech, ASCE 1965;91:61. [2] Choros J, Adams GG. A steadily moving load on an elastic beam resting on a tensionless Winkler subgrade. J Appl Mech, ASME 1979;46(1):175–80. [3] Conte SD, de Boor C. Elementary numerical analysis: an algorithmic approach. 3rd ed. New York: McGraw-Hill; 1980. [4] Criner HE, McCann GD. Rails on elastic subgrades under the influence of high-speed traveling loads. J Appl Mech, ASME 1953;20(1):13. [5] Elattary MA. Moving loads on an infinite plate strip of constant thickness. J Phys D: Appl Phys 1991;24(4):541–6. [6] Florence AL. Traveling force on a Timoshenko beam. J Appl Mech, ASME 1965;32(2):351–8. [7] Fryba L. Vibration of solids and structures under moving loads. Groningen, Netherlands: Noordhoff International; 1977. [8] Hardy MSA, Cebon D. Importance of speed and frequency in flexible pavement response. J Eng Mech, ASCE 1994;120(3):463–82. [9] Huang CC. Traveling loads on a viscoelastic Timoshenko beam. J Appl Mech, ASME 1977;44(1):183–4. [10] Jezequel L. Response of periodic systems to a moving load. J Appl Mech, ASME 1981;48(3):613–8. [11] Kenney JT. Steady state vibrations of beam on elastic subgrade for moving loads. J Appl Mech, ASME 1954; 21(4):359–64.
8
L. Sun / Computers and Structures 80 (2002) 1–8
[12] Lee HP. Dynamic response of a beam with intermediate point constrains subject to a moving load. J Sound Vibr 1994;171(3):361–8. [13] Morse PM, Feshbach H. Methods of theoretical physics: Part I and II. New York: McGraw-Hill Book Company Inc; 1953. [14] Mulcahy TM. Steady-state response of a beam on idealized strain-hardening subgrades for a moving load. J Appl Mech, ASME 1973;42(4):1040–4. [15] Pan G, Atluri SN. Dynamic response of finite sized elastic runways subjected to moving loads: a coupled BEM/FEM approach. Int J Num Meth Eng 1995;38(18): 3143–66. [16] Patil S. Response of infinite railroad track to vibrating mass. J Eng Mech, ASCE 1988;114(4):688–703. [17] Saito H, Terasawa T. Steady-state vibrations of a beam on a Pasternak subgrade for moving loads. J Appl Mech, ASME 1980;47(4):879–83. [18] Steele CR. The finite beam with a moving load. J Appl Mech, ASME 1967;34(1):111.
[19] Steele CR. The Timoshenko beam with a moving load. J Appl Mech, ASME 1968;35(3):481–8. [20] Sun L. Theoretical investigations on vehicle-ground dynamic interaction. Final Report prepared for National Science Foundation of China, Southeast University, Nanjing, China. [21] Sun L, Deng X. Transient response for infinite plate on Winkler subgrade by a moving distributed load. Chin J Appl Mech 1997;14(2):72–8. [22] Sun L, Deng X. Dynamic analysis of infinite beam under the excitation of moving line loads. Appl Math Mech 1998;19(4):367–73. [23] Sun L, Greenberg B. Dynamic response of linear systems to moving stochastic sources. J Sound Vibr 2000;229(4): 957–72. [24] Yoder EJ, Witczak MW. Principles of pavement design. New York: Wiely; 1975. [25] Adams GG, Bogy DB. Steady solutions for moving loads on elastic beams with one-sided constraints. J Appl Mech 1975;42(4):800–5.