[go: up one dir, main page]

Academia.eduAcademia.edu
RaoCh12ff.qxd 10.06.08 13:48 Page 856 Aurel Boreslav Stodola (1859–1942) was a Swiss engineer who joined the Swiss Federal Institute of Technology in Zurich in 1892 to occupy the chair of thermal machinery. He worked in several areas including machine design, automatic controls, thermodynamics, rotor dynamic, and steam turbines. He published one of the most outstanding books, namely, Die Dampfturbin, at the turn of the century. This book discussed not only the thermodynamic issues involved in turbine design but also the aspects of fluid flow, vibration, stress analysis of plates, shells and rotating discs, thermal stresses and stress concentrations at holes and fillets, and was translated into many languages. The approximate method he presented for the computation of natural frequencies of beams has become known as the Stodola method. (Photo courtesy of Applied Mechanics Reviews.) C H A P T E R 1 2 Finite Element Method 12.1 Introduction The finite element method is a numerical method that can be used for the accurate solution of complex mechanical and structural vibration problems [12.1, 12.2]. In this method, the actual structure is replaced by several pieces or elements, each of which is assumed to behave as a continuous structural member called a finite element. The elements are assumed to be interconnected at certain points known as joints or nodes. Since it is very difficult to find the exact solution (such as the displacements) of the original structure under the specified loads, a convenient approximate solution is assumed in each finite element. The idea is that if the solutions of the various elements are selected properly, they can be made to converge to the exact solution of the total structure as the element size is reduced. During the solution process, the equilibrium of forces at the joints and the compatibility of displacements between the elements are satisfied so that the entire structure (assemblage of elements) is made to behave as a single entity. The basic procedure of the finite element method, with application to simple vibration problems, is presented in this chapter. The element stiffness and mass matrices, and force vectors are derived for a bar element, a torsion element, and a beam element. The transformation of element matrices and vectors from the local to the global coordinate system is presented. The equations of motion of the complete system of finite elements and the 856 RaoCh12ff.qxd 10.06.08 13:48 Page 857 12.2 EQUATIONS OF MOTION OF AN ELEMENT 857 incorporation of the boundary conditions are discussed. The concepts of consistent and lumped mass matrices are presented along with a numerical example. Finally, a computer program for the eigenvalue analysis of stepped beams is presented. Although the techniques presented in this chapter can be applied to more complex problems involving two- and three-dimensional finite elements, only the use of one-dimensional elements is considered in the numerical treatment. 12.2 Equations of Motion of an Element For illustration, the finite element model of a plano-milling machine structure (Fig. 12.1a) is shown in Fig. 12.1(b). In this model, the columns and the overarm are represented by triangular plate elements and the cross slide and the tool holder are represented by beam elements [12.3]. The elements are assumed to be connected to each other only at the joints. The displacement within an element is expressed in terms of the displacements at the corners or joints of the element. In Fig. 12.1(b), the transverse displacement within a typical element e is assumed to be w(x, y, t). The values of w, (0w)/(0x), and (0w)/(0y) at joints 1, 2, and 3—namely w(x1 y1, t), (0w)/(0x)(x1, y1, t), (0w)/(0y)(x1, y1, t), Á , (0w)/(0y)(x3, y3, t)—are treated as unknowns and are denoted as w1(t), w2(t), w3(t), Á , w9(t). The displacement w(x, y, t) can be expressed in terms of the unknown joint displacements wi(t) in the form w(x, y, t) = a Ni(x, y) wi(t) n (12.1) i=1 where Ni(x, y) is called the shape function corresponding to the joint displacement wi(t) and n is the number of unknown joint displacements (n = 9 in Fig. 12.1b). If a distributed w7(t) w1(t) w(x, y, t) w9(t) Tool holder Overarm w3(t) Column Column 1 Element e Cross-slide Cutter w2(t) w4(t) w6(t) Plate elements 2 w5(t) 3 f(x, y, t) e Beam elements Bed Fy y w8(t) Fz Fx Cutting forces z x (a) Plano-milling machine structure FIGURE 12.1 Finite element modeling. (b) Finite element model RaoCh12ff.qxd 858 10.06.08 13:48 Page 858 CHAPTER 12 FINITE ELEMENT METHOD load f (x, y, t) acts on the element, it can be converted into equivalent joint forces fi(t) (i = 1, 2, Á , 9). If concentrated forces act at the joints, they can also be added to the appropriate joint force fi(t). We shall now derive the equations of motion for determining the joint displacements wi(t) under the prescribed joint forces fi(t). By using Eq. (12.1), the kinetic energy T and the strain energy V of the element can be expressed as # 1 :# : T = WT [m] W (12.2) 2 1: : V = WT [k] W (12.3) 2 where # w1(t) w1(t) dw1/dt # w2(t) dw2/dt w2(t) # ! : . . . W= e . u = e . u W = e . u, . . . # dwn/dt wn(t) wn(t) and [m] and [k] are the mass and stiffness matrices of the element. By substituting Eqs. (12.2) and (12.3) into Lagrange’s equations, Eq. (6.44), the equations of motion of the finite element can be obtained as $ : : : [m]W + [k] W = f (12.4) $ : : where f is the vector of joint forces and W is the vector of joint accelerations given by $ d2w1/dt2 w1 $ w2 d2w2/dt2 $ : . . W = e . u = e u . . . $ d2wn/dt2 wn Note that the shape of the finite elements and the number of unknown joint displacements may differ for different applications. Although the equations of motion of a single element, Eq. (12.4), are not useful directly (as our interest lies in the dynamic response of the assemblage of elements), the mass matrix [m], the stiffness matrix [k], and the joint force : vector f of individual elements are necessary for the final solution. We shall derive the element mass and stiffness matrices and the joint force vectors for some simple onedimensional elements in the next section. 12.3 Mass Matrix, Stiffness Matrix, and Force Vector 12.3.1 Bar Element Consider the uniform bar element shown in Fig. 12.2. For this one-dimensional element, the two end points form the joints (nodes). When the element is subjected to axial loads f1(t) and f2(t), the axial displacement within the element is assumed to be linear in x as u (x, t) = a (t) + b (t)x (12.5) RaoCh12ff.qxd 10.06.08 13:48 Page 859 12.3 MASS MATRIX, STIFFNESS MATRIX, AND FORCE VECTOR 859  , E, A Joint 1 u1(t) u(x, t) f1(t) f (x, t) u2(t) x Joint 2 x f2(t) l FIGURE 12.2 Uniform bar element. When the joint displacements u1(t) and u2(t) are treated as unknowns, Eq. (12.5) should satisfy the conditions u (l, t) = u2 (t) u (0, t) = u1(t), (12.6) Equations (12.5) and (12.6) lead to a (t) = u1(t) and a (t) + b (t)l = u2 (t) or b(t) = u2(t) - u1(t) l (12.7) Substitution for a(t) and b(t) from Eq. (12.7) into Eq. (12.5) gives x x b u1(t) + u2(t) l l (12.8) u (x, t) = N1(x) u1(t) + N2(x) u2(t) (12.9) u (x, t) = a1 - or where N1(x) = a1 - x b, l N2(x) = x l (12.10) are the shape functions. The kinetic energy of the bar element can be expressed as 1 0 u(x, t) 2 f dx rA e 2 L0 0t l T(t) = = = l 1 x du2(t) 2 x du1(t) + a b f dx rA e a 1 - b 2 L0 l dt l dt 1 rAl # 2 # # # (u1 + u1u2 + u22) 2 3 (12.11) RaoCh12ff.qxd 860 10.06.08 13:49 Page 860 CHAPTER 12 FINITE ELEMENT METHOD where du1(t) # u1 = , dt du2(t) # u2 = , dt r is the density of the material and A is the cross-sectional area of the element. By expressing Eq. (12.11) in matrix form, T(t) = where # 1 :# T u (t) [m] : u (t) 2 (12.12) # # u1(t) f u (t) = e # u2(t) : and the superscript T indicates the transpose, the mass matrix [m] can be identified as [m] = rAl 2 c 6 1 The strain energy of the element can be written as 1 d 2 (12.13) 1 0 u(x, t) 2 f dx EA e 2 L0 0x l V(t) = 2 1 1 1 EA e - u1(t) + u2(t) f dx 2 L0 l l l = = 1 EA 2 (u1 - 2u1u2 + u22) 2 l (12.14) where u1 = u1(t), u2 = u2(t), and E is Young’s modulus. By expressing Eq. (12.14) in matrix form as V(t) = where u (t) ! u(t) = e 1 f u2(t) 1 ! T ! u (t) [k] u (t) 2 and the stiffness matrix [k] can be identified as [k] = (12.15) ! u(t)T = 5u1(t) u2(t)6 EA 1 c l -1 -1 d 1 (12.16) RaoCh12ff.qxd 10.06.08 13:49 Page 861 12.3 MASS MATRIX, STIFFNESS MATRIX, AND FORCE VECTOR The force vector f = e : 861 f1(t) f f2(t) can be derived from the virtual work expression. If the bar is subjected to the distributed force f (x, t), the virtual work dW can be expressed as dW(t) = = L0 l L0 l = a f(x, t) du(x, t) dx f(x, t) e a 1 - L0 + a l f(x, t) a1 - x x b du1(t) + a b du2(t) f dx l l x b dx b du1(t) l x f(x, t) a b dx b du2(t) l L0 l (12.17) By expressing Eq. (12.17) in matrix form as : ! dW(t) = du(t)T f (t) K f1(t) du1(t) + f2(t) du2(t) the equivalent joint forces can be identified as f1(t) = L0 l f(x, t) a1 - x b dx l x f2(t) = f(x, t) a b dx l L0 l t (12.18) (12.19) Consider a uniform torsion element with the x axis taken along the centroidal axis, as 12.3.2 shown in Fig. 12.3. Let Ip denote the polar moment of inertia about the centroidal axis and Torsion Element , Ip, G, A Joint 1 Joint 2 x x 1(t) f1(t)  (x, t) f (x, t) 2(t) f2(t) l FIGURE 12.3 Uniform torsion element. RaoCh12ff.qxd 862 10.06.08 13:49 Page 862 CHAPTER 12 FINITE ELEMENT METHOD GJ represent the torsional stiffness (J = Ip for a circular cross section). When the torsional displacement (rotation) within the element is assumed to be linear in x as u(x, t) = a(t) + b(t)x (12.20) and the joint rotations u1(t) and u2(t) are treated as unknowns, Eq. (12.20) can be expressed, by proceeding as in the case of a bar element, as u (x, t) = N1(x) u1(t) + N2(x) u2(t) (12.21) where N1(x) and N2(x) are the same as in Eq. (12.10). The kinetic energy, the strain energy, and the virtual work for pure torsion are given by 1 0u(x, t) 2 f dx rIp e 2 L0 0t l T(t) = (12.22) 1 0u(x, t) 2 f dx GJ e 2 L0 0x (12.23) f(x, t) du (x, t) dx (12.24) l V(t) = dW(t) = L0 l where r is the mass density and f(x, t) is the distributed torque per unit length. Using the procedures employed in Section 12.3.1, we can derive the element mass and stiffness matrices and the force vector: [m] = [k] = rIp l 2 c 6 1 GJ 1 c l -1 f (t) L0 f = e 1 f = d f2(t) : 12.3.3 Beam Element 1 d 2 -1 d 1 f(x, t) a1 - l (12.25) (12.26) x b dx l x f(x, t) a b dx l L0 l t (12.27) We now consider a beam element according to the Euler-Bernoulli theory.1 Figure 12.4 shows a uniform beam element subjected to the transverse force distribution f (x, t). In this case, the joints undergo both translational and rotational displacements, so the unknown joint displacements are labeled as w1(t), w2(t), w3(t), and w4(t). Thus there will be linear joint forces f1(t) and f3(t) corresponding to the linear joint displacements w1(t) and w3(t) and rotational joint forces (bending moments) f2(t) and f4(t) corresponding to the rotational 1 The beam element, according to the Timoshenko theory, was considered in Refs. [12.4–12.7]. RaoCh12ff.qxd 10.06.08 13:49 Page 863 12.3 MASS MATRIX, STIFFNESS MATRIX, AND FORCE VECTOR f1(t) f3(t) w(x, t) f2(t) 863 w1(t) f4(t) f(x, t) w2(t) w3(t) w4(t) x x , A, I, E l Joint 1 Joint 2 FIGURE 12.4 Uniform beam element. joint displacements w2(t) and w4(t), respectively. The transverse displacement within the element is assumed to be a cubic equation in x (as in the case of static deflection of a beam): w(x, t) = a(t) + b(t)x + c(t)x2 + d(t)x3 (12.28) The unknown joint displacements must satisfy the conditions 0w (0, t) = w2(t) 0x t 0w (l, t) = w4(t) 0x w(0, t) = w1(t), w(l, t) = w3(t), (12.29) Equations (12.28) and (12.29) yield a(t) = w1(t) b(t) = w2(t) 1 c(t) = 2 [ - 3w1(t) - 2w2(t)l + 3w3(t) - w4(t)l] l 1 d(t) = 3 [2w1(t) + w2(t)l - 2w3(t) + w4(t)l] l (12.30) By substituting Eqs. (12.30) into Eq. (12.28), we can express w(x, t) as w(x, t) = a1 - 3 + a3 x2 x2 x3 x x3 2 + 2 b w (t) + a + b lw2(t) 1 l l2 l3 l2 l3 x3 x2 x3 x2 2 b w (t) + a + b lw4(t) 3 l2 l3 l2 l3 (12.31) This equation can be rewritten as w(x, t) = a Ni (x)wi (t) 4 i=1 (12.32) RaoCh12ff.qxd 864 10.06.08 13:49 Page 864 CHAPTER 12 FINITE ELEMENT METHOD where Ni(x) are the shape functions given by x 2 x 3 N1(x) = 1 - 3a b + 2a b l l x 2 x 3 N2(x) = x - 2l a b + la b l l x 2 x 3 N3(x) = 3a b - 2a b l l (12.33) (12.34) (12.35) x 2 x 3 N4(x) = - l a b + la b l l (12.36) The kinetic energy, bending strain energy, and virtual work of the element can be expressed as T(t) = l # 1 :# T 1 0w(x, t) 2 : f dx K w (t) [m]w (t) rAe 2 L0 2 0t 1 1 ! ! 0 2w(x, t) 2 V(t) = f dx K w (t)T [k]w (t) EI e 2 2 L0 2 0x (12.37) l dW(t) = L0 (12.38) l : ! f(x, t) dw (x, t) dx K dw (t)T f (t) (12.39) where r is the density of the beam, E is Young’s modulus, I is the moment of inertia of the cross section, A is the area of cross section, and w1(t) w (t) ! w(t) = µ 2 ∂ , w3(t) w4(t) dw1/dt #: dw /dt w (t) = µ 2 ∂ dw3/dt dw4/dt dw1(t) dw2(t) ! ∂, dw(t) = µ dw3(t) dw4(t) f1(t) f (t) f (t) = µ 2 ∂ f3(t) f4(t) : By substituting Eq. (12.31) into Eqs. (12.37) to (12.39) and carrying out the necessary integrations, we obtain 156 rAl 22l [m] = ≥ 420 54 - 13l 22l 4l2 13l - 312 54 13l 156 - 22l -13l -3l2 ¥ -22l 4l2 (12.40) RaoCh12ff.qxd 10.06.08 13:49 Page 865 12.4 TRANSFORMATION OF ELEMENT MATRICES AND VECTORS 12 EI 6l [k] = 3 ≥ -12 l 6l L0 fi (t) = 12.4 - 12 - 6l 12 -6l 6l 4l2 -6l 2l2 865 6l 2l2 ¥ - 6l 4l2 (12.41) i = 1, 2, 3, 4 (12.42) l f(x, t) Ni(x) dx, Transformation of Element Matrices and Vectors As stated earlier, the finite element method considers the given dynamical system as an assemblage of elements. The joint displacements of an individual element are selected in a convenient direction, depending on the nature of the element. For example, for the bar element shown in Fig. 12.2, the joint displacements u1(t) and u2(t) are chosen along the axial direction of the element. However, other bar elements can have different orientations in an assemblage, as shown in Fig. 12.5. Here x denotes the axial direction of an individual element and is called a local coordinate axis. If we use u1(t) and u2(t) to denote the joint displacements of different bar elements, there will be one joint displacement at joint 1, three at joint 2, two at joint 3, and 2 at joint 4. However, the displacements of joints can be specified more conveniently using reference or global coordinate axes X and Y. Then the displacement components of joints parallel to the X and Y axes can be used as the joint displacements in the global coordinate system. These are shown as Ui(t), i = 1, 2, Á , 8 in Fig. 12.5. The joint displacements in the local and the global coordinate system for a Y U8(t) u2(t) U4(t) 3 u2(t) 1 u2(t) x U3(t) x u1(t) 2 U2(t) U7(t) Load u1(t) 2 4 4 U6(t) x x u1(t) 1 u2(t) U1(t) u1(t) 3 X U5(t) FIGURE 12.5 A dynamical system (truss) idealized as an assemblage of four bar elements. RaoCh12ff.qxd 866 10.06.08 13:49 Page 866 CHAPTER 12 FINITE ELEMENT METHOD U2j (t) u2(t) j U2i(t) e u1(t) Y x  i X U2j1(t) U2i1(t) x  local coordinate axis X,Y  global coordinate axes u1(t), u2(t)  local joint displacements U2i1(t), ... , U2j (t)  global joint displacements FIGURE 12.6 Local and global joint displacements of element e. typical bar element e are shown in Fig. 12.6. The two sets of joint displacements are related as follows: u1(t) = U2i - 1(t) cos u + U2i(t) sin u u2(t) = U2j - 1(t) cos u + U2j(t) sin u (12.43) ! ! u(t) = [l] U(t) (12.44) These can be rewritten as where [l] is the coordinate transformation matrix given by [l] = c cos u 0 sin u 0 0 cos u 0 d sin u (12.45) ! ! and u(t) and U(t) are the vectors of joint displacements in the local and the global coordinate system, respectively, and are given by u (t) ! u(t) = e 1 f, u2(t) U2i - 1(t) ! U (t) U(t) = µ 2i ∂ U2j - 1(t) U2j (t) It is useful to express the mass matrix, stiffness matrix, and joint force vector of an element in terms of the global coordinate system while finding the dynamical response of the RaoCh12ff.qxd 10.06.08 13:49 Page 867 12.4 TRANSFORMATION OF ELEMENT MATRICES AND VECTORS 867 complete system. Since the kinetic and strain energies of the element must be independent of the coordinate system, we have T(t) = # # 1 :# T 1 :# : u (t) [m] : u (t) = U (t)T [ m ] U(t) 2 2 (12.46) V(t) = ! 1 ! T 1 ! ! u(t) [k] u(t) = U(t)T [k]U(t) 2 2 (12.47) where [m] and [k] denote the# element mass and stiffness matrices, respectively, in the : global coordinate system and U (t) is the vector of joint velocities in the global coordinate # : system, related to u (t) as in Eq. (12.44): # # : : (12.48) u (t) = [l]U (t) By inserting Eqs. (12.44) and (12.48) into (12.46) and (12.47), we obtain T(t) = # # 1 :# T T 1 :# : : U(t) [l] [m][l] U(t) K U(t)T[ m ] U(t) 2 2 (12.49) V(t) = : 1: : 1: T T U(t) [l] [k][l] U(t) K U(t)T[k] U(t) 2 2 (12.50) Equations (12.49) and (12.50) yield [ m ] = [l]T[m][l] (12.51) [k] = [l]T[k][l] (12.52) Similarly, by equating the virtual work in the two coordinate systems, ! : : ! dW(t) = d u(t)T f (t) = d U(t)T f (t) (12.53) : we find the vector of element joint forces in the global coordinate system f (t): : : f (t) = [l]T f (t) (12.54) Equations (12.51), (12.52), and (12.54) can be used to obtain the equations of motion of a single finite element in the global coordinate system: $ ! : : [ m ] U(t) + [k] U(t) = f (t) (12.55) Although this equation is not of much use, since our interest lies in the equations of motion : of an assemblage of elements, the matrices [ m ] and [k] and the vector f are useful in deriving the equations of motion of the complete system, as indicated in the following section. RaoCh12ff.qxd 868 12.5 10.06.08 13:49 Page 868 CHAPTER 12 FINITE ELEMENT METHOD Equations of Motion of the Complete System of Finite Elements Since the complete structure is considered an assemblage of several finite elements, we shall now extend the equations of motion obtained for single finite elements in the global system to the complete structure. We shall denote the joint displacements of the complete structure in the global coordinate system as U1(t), U2(t), Á , UM(t) or, equivalently, as a column vector: U1(t) U2(t) . ! U ' (t) = f . v . UM(t) For convenience, we shall denote the quantities pertaining to an element e in the assemblage by the superscript e. Since the joint displacements of any element e can ! be identified ! in the vector of joint displacements of the complete structure, the vectors U(e)(t) and U ' (t) are related ! ! U(e)(t) = [A(e)] U (12.56) ' (t) where [A(e)] is a rectangular matrix composed of zeros and ones. For example, for element 1 in Fig. 12.5, Eq. (12.56) becomes U1(t) 1 ! U (t) 0 U(1)(t) K µ 2 ∂ = ≥ U3(t) 0 U4(t) 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 U1(t) 0 U2(t) . 0 ¥f . v 0 . 0 U8(t) (12.57) The kinetic energy of the complete structure can be obtained by adding the kinetic energies of individual elements # E # 1: : T = a U (e)T [ m ] U (e) e=1 2 (12.58) where E denotes the number of finite elements in the assemblage. By differentiating Eq. (12.56), the relation between the velocity vectors can be derived: # #! : (t) U(e)(t) = [A(e)] U (12.59) ' Substitution of Eq. (12.59) into (12.58) leads to T = #! 1 E # ! T (e) T (e) (e) U ' [A ] [ m ] [A ] U ' a 2 e=1 (12.60) RaoCh12ff.qxd 10.06.08 13:49 Page 869 12.6 INCORPORATION OF BOUNDARY CONDITIONS 869 The kinetic energy of the complete # ! structure can also be expressed in terms of joint velocities of the complete structure U ' T = #! 1 # ! (e)T [M] U U ' ' 2 (12.61) where [M ' ] is called the mass matrix of the complete structure. A comparison of Eqs. (12.60) and (12.61) gives the relation2 (e) T (e) (e) [M ' ] = a [A ] [ m ][A ] E (12.62) e=1 Similarly, by considering strain energy, the stiffness matrix of the complete structure, [K ' ], can be expressed as (e) (e) T (e) [K ' ] = a [A ] [k ] [A ] E (12.63) e=1 Finally the !consideration of virtual work yields the vector of joint forces of the complete structure, F ': E ! : (e) T (e) F ' = a [A ] f (12.64) e=1 Once the mass and stiffness matrices and the force vector are known, Lagrange’s equations of motion for the complete structure can be expressed as $! ! : [M] (12.65) ' U ' + [K] U = F ' ! Note that the joint force vector F ' in Eq. (12.65) was generated by considering only the distributed loads acting on the various elements. If there is any concentrated load act! ing along the joint displacement Ui(t), it must be added to the ith component of F . ' 12.6 Incorporation of Boundary Conditions In the preceding derivation, no joint was assumed to be fixed. Thus the complete structure is capable of undergoing rigid body motion under the joint forces. This means that [K '] is a singular matrix (see Section 6.12). Usually the structure is supported such that the 2 An alternative procedure can be used for the assembly of element matrices. In this procedure, each of the rows and columns of the element (mass or stiffness) matrix is identified by the corresponding degree of freedom in the assembled structure. Then the various entries of the element matrix can be placed at their proper locations in the overall (mass or stiffness) matrix of the assembled system. For example, the entry belonging to the ith row (identified by the degree of freedom p) and the jth column (identified by the degree of freedom q) of the element matrix is to be placed in the pth row and qth column of the overall matrix. This procedure is illustrated in Example 12.3. RaoCh12ff.qxd 870 10.06.08 13:49 Page 870 CHAPTER 12 FINITE ELEMENT METHOD displacements are zero at a number of joints, to avoid rigid body motion of the structure. A simple method of incorporating the zero displacement conditions is to eliminate ! the cor] [K ] F responding rows and columns from the matrices [M and and the vector ' . The final ' ' equations of motion of the restrained structure can be expressed as $ : : : [M] U + [K] U = F N*NN*1 N*N N*1 N*1 (12.66) where N denotes the number of free joint displacements of the structure. Note the following points concerning finite element analysis: 1. The approach used in the above presentation is called the displacement method of finite element analysis because it is the displacements of elements that are directly approximated. Other methods, such as the force method, the mixed method, and hybrid methods, are also available [12.8, 12.9]. 2. The stiffness matrix, mass matrix, and force vector for other finite elements, including two-dimensional and three-dimensional elements, can be derived in a similar manner, provided the shape functions are known [12.1, 12.2]. 3. In the Rayleigh-Ritz method discussed in Section 8.8, the displacement of the continuous system is approximated by a sum of assumed functions, where each function denotes a deflection shape of the entire structure. In the finite element method, an approximation using shape functions (similar to the assumed functions) is also used for a finite element instead of the entire structure. Thus the finite element procedure can also be considered a Rayleigh-Ritz method. 4. Error analysis of the finite element method can also be conducted [12.10]. Analysis of a Bar EXAMPLE 12.1 Consider a uniform bar, of length 0.5 m, area of cross section 5 * 10-4 m2, Young’s modulus 200 GPa, and density 7850 kg/m3, which is fixed at the left end, as shown in Fig. 12.7. a. b. Find the stress induced in the bar under an axial static load of 1000 N applied at joint 2 along u2. Find the natural frequency of vibration of the bar. Use a one-element idealization. 1 2 u1 0.5 m FIGURE 12.7 Uniform bar with two degrees of freedom. u2 RaoCh12ff.qxd 10.06.08 13:49 Page 871 12.6 INCORPORATION OF BOUNDARY CONDITIONS 871 Solution a. Using the stiffness matrix of a bar element, Eq. (12.16), the equilibrium equations can be written as AE 1 c l -1 - 1 u1 f d e f = e 1f 1 u2 f2 (E.1) - 1 u1 f de f = e 1 f 1 u2 1000 (E.2) With A = 5 * 10-4, E = 2 * 1011, l = 0.5, f2 = 1000, Eq. (E.1) becomes 2 * 108 c 1 -1 where u1 is the displacement and f1 is the unknown reaction at joint 1. To incorporate the boundary condition u1 = 0, we delete the first scalar equation (first row) and substitute u1 = 0 in the resulting Eq. (E.2). This gives 2 * 108 u2 = 1000 or u2 = 500 * 10-8 m (E.3) The stress (s) - strain (e) relation gives s = Ee = E u2 - u1 ¢l = Ea b l l where ¢l = u2 - u1 denotes the change in length of the element and (E.4) ¢l indicates the strain. l Equation (E.4) yields s = 2 * 1011 a b. 500 * 10-8 - 0 b = 2 * 106 Pa 0.5 (E.5) Using the stiffness and mass matrices of the bar element, Eqs. (12.16) and (12.13), the eigenvalue problem can be expressed as AE 1 c l -1 r Al 2 -1 U1 d e f = v2 c 1 U2 6 1 1 U1 de f 2 U2 (E.6) where v is the natural frequency and U1 and U2 are the amplitudes of vibration of the bar at joints 1 and 2, respectively. To incorporate the boundary condition U1 = 0, we delete the first row and first column in each of the matrices and vectors and write the resulting equation as r Al AE U2 = v2 (2) U2 l 6 or v = 3(2 * 1011) 3E = = 17, 485.2076 rad/s 2 Br l B 7850 (0.5)2 (E.7) ■ RaoCh12ff.qxd 872 10.06.08 13:49 Page 872 CHAPTER 12 FINITE ELEMENT METHOD Natural Frequencies of a Simply Supported Beam EXAMPLE 12.2 Find the natural frequencies of the simply supported beam shown in Fig. 12.8(a) using one finite element. Solution: Since the beam is idealized using only one element, the element joint displacements are the same in both local and global systems, as indicated in Fig. 12.8(b). The stiffness and mass matrices of the beam are given by 12 EI 6l (1) ≥ [K ' ] = [K ] = l3 -12 6l 156 rAl 22l (1) ≥ [M ' ] = [M ] = 54 420 -13l 6l 4l2 - 6l 2l2 -12 - 6l 12 - 6l 22l 4l2 13l - 3l2 54 13l 156 -22l 6l 2l2 ¥ - 6l 2 4l -13l -3l2 ¥ - 22l 2 4l (E.1) (E.2) and the vector of joint displacements by W1 w (1) 1 ! W2 w (1) W ∂ K µ 2(1) ∂ ' = µ W3 w3 W4 w (1) 4 (E.3) The boundary conditions corresponding to the simply supported ends (W1 = 0 and W3 = 0) can be incorporated3 by deleting the rows and columns corresponding to W1 and W3 in Eqs. (E.1) and (E.2). l (a) w1(1)  W1 w2(1)  W2 w(x, t) w3(1)  W3 w4(1)  W4 x 1 1 2 2 l (1)  l (b) FIGURE 12.8 3 Simply supported beam. The bending moment cannot be set equal to zero at the simply supported ends explicitly, since there is no degree of freedom (joint displacement) involving the second derivative of the displacement w. RaoCh12ff.qxd 10.06.08 13:49 Page 873 12.6 INCORPORATION OF BOUNDARY CONDITIONS 873 This leads to the overall matrices [K] = [M] = 2EI 2 c l 1 1 d 2 rAl3 4 c 420 - 3 -3 d 4 (E.4) (E.5) and the eigenvalue problem can be written as c 2EI 2 c l 1 rAl3v2 1 4 d c 2 420 -3 -3 W 0 d d e 2f = e f 4 W4 0 (E.6) By multiplying throughout by l/(2EI), Eq. (E.6) can be expressed as c 2 - 4l 1 + 3l 1 + 3l W 0 d e 2f = e f 2 - 4l W4 0 (E.7) where l = rAl4v2 840EI (E.8) By setting the determinant of the coefficient matrix in Eq. (E.7) equal to zero, we obtain the frequency equation ` 2 - 4l 1 + 3l 1 + 3l ` = (2 - 4l)2 - (1 + 3l)2 = 0 2 - 4l (E.9) The roots of Eq. (E.9) give the natural frequencies of the beam as 1 7 or l2 = 3 or l1 = v1 = a v2 = a 120EI 1/2 b rAl4 2520EI 1/2 b rAl4 (E.10) (E.11) These results can be compared with the exact values (see Fig. 8.15): v1 = a 97.41EI 1/2 b , rAl4 v2 = a 1558.56EI 1/2 b rAl4 (E.12) ■ RaoCh12ff.qxd 874 10.06.08 13:49 Page 874 CHAPTER 12 FINITE ELEMENT METHOD Stresses in a Two-Bar Truss EXAMPLE 12.3 Find the stresses developed in the two members of the truss shown in Fig. 12.9(a), under a vertical load of 200 lb at joint 3. The areas of cross section are 1 in.2 for member 1 and 2 in.2 for member 2, and the Young’s modulus is 30 * 106 psi. Y 1 Element 1 3 10 in. Element 2 200 lb 5 in. X 2 10 in. (a) U2 U1 1 1 x X U6 3 U4 x U3 2 2 X (b) FIGURE 12.9 Two bar truss. F 3  200 lb U5 RaoCh12ff.qxd 10.06.08 13:49 Page 875 12.6 INCORPORATION OF BOUNDARY CONDITIONS 875 Solution Approach: Derive the static equilibrium equations and solve them to find the joint displacements. Use the elasticity relations to find the element stresses. Each member is to be treated as a bar element. From Fig. 12.9(a), the coordinates of the joints can be found as (X1, Y1) = (0, 10) in.; (X2, Y2) = (0, 0) in.; (X3, Y3) = (10, 5) in. The modeling of the truss as an assemblage of two bar elements and the displacement degrees of freedom of the joints are shown in Fig. 12.9(b). The lengths of the elements can be computed from the coordinates of the ends (joints) as l(1) = 5(X3 - X1)2 + (Y3 - Y1)261/2 = 5(10 - 0)2 + (5 - 10)261/2 = 11.1803 in. l(2) = 5(X3 - X2)2 + (Y3 - Y2)261/2 = 5(10 - 0)2 + (5 - 0)261/2 = 11.1803 in. (E.1) The element stiffness matrices in the local coordinate system can be obtained as [k(1)] = A(1)E (1) 1 c -1 l(1) = 2.6833 * 106 c [k(2)] = A(2)E (2) 1 c -1 l(2) = 5.3666 * 106 c (1)(30 * 106) 1 -1 c d = 11.1803 -1 1 -1 d 1 (2)(30 * 106) 1 -1 c d = 11.1803 -1 1 -1 d 1 1 -1 1 -1 -1 d 1 -1 d 1 (E.2) The angle between the local x-coordinate and the global X-coordinate is given by cos u1 = sin u1 = X3 - X1 l(1) Y3 - Y1 cos u2 = sin u2 = l(1) 10 - 0 = 0.8944 11.1803 t for element 1 5 - 10 = - 0.4472 = 11.1803 X3 - X2 l(2) Y3 - Y2 l(2) = 10 - 0 = 0.8944 11.1803 t for element 2 5 - 0 = 0.4472 = 11.1803 (E.3) = (E.4) RaoCh12ff.qxd 876 10.06.08 13:49 Page 876 CHAPTER 12 FINITE ELEMENT METHOD The stiffness matrices of the elements in the global (X, Y) coordinate system can be derived as [k (1)] = [l(1)]T[k(1)][l(1)] 1 0.8 6 = 2.6833 * 10 E - 0.4 - 0.8 0.4 2 - 0.4 0.2 0.4 - 0.2 5 - 0.8 0.4 0.8 - 0.4 6 0.4 1 - 0.2 U 2 - 0.4 5 0.2 6 (E.5) 4 0.4 0.2 - 0.4 - 0.2 5 - 0.8 - 0.4 0.8 0.4 6 - 0.4 3 - 0.2 U 4 0.4 5 0.2 6 (E.6) [k (2)] = [l(2)]T[k(2)][l(2)] 3 0.8 = 5.3666 * 106 E 0.4 - 0.8 - 0.4 where [l(1)] = c = c [l(2)] = c = c cos u1 0 0.8944 0 sin u1 0 0 cos u1 - 0.4472 0 cos u2 0 sin u2 0 0.8944 0 0.4472 0 0 d sin u1 0 d - 0.4472 0 0.8944 0 cos u2 0 0.8944 (E.7) 0 d sin u2 0 d 0.4472 (E.8) Note that the top and right-hand sides of Eqs. (E.5) and (E.6) denote the global degrees of freedom corresponding to the rows and columns of the respective stiffness matrices. The assembled stiffness (1) (2) matrix of the system, [K ' ] can be obtained, by placing the elements of [k ] and [k ] at their proper places in [K ' ], as 1 0.8 - 0.4 [K ] = 2.6833 * 106 H ' 2 - 0.4 0.2 3 4 0.8 0.4 - 0.8 - 0.4 - 0.8 0.4 1.6 0.8 - 1.6 0.4 - 0.2 - 0.8 5 - 0.8 0.4 - 1.6 - 0.8 (0.8 +1.6) (- 0.4 +0.8) 6 0.4 1 - 0.2 2 - 0.8 3 - 0.4 4 X (- 0.4 +0.8) 5 (0.2 +0.4) 6 (E.9) RaoCh12ff.qxd 10.06.08 13:49 Page 877 12.6 INCORPORATION OF BOUNDARY CONDITIONS 877 The assembled force vector can be written as FX1 FY1 ! FX2 F v ' = f FY2 FX3 FY3 (E.10) where, in general, (FXi, FYi) denote the forces applied at joint i along (X, Y) directions. Specifically, (FX1, FY1) and (FX2, FY2) represent the reactions at joints 1 and 2, while (FX3, FY3) = (0, - 200) lb shows the external forces applied at joint 3. By applying the boundary conditions U1 = U2 = U3 = U4 = 0 (i.e., by deleting the rows and columns 1, 2, 3, and 4 in Eqs. E.9 and E.10), we get the final assembled stiffness matrix and the force vector as 5 2.4 [K] = 2.6833 * 106 c 0.4 ! F = e 0 5 f -200 6 6 0.4 5 d 0.6 6 (E.11) (E.12) The equilibrium equations of the system can be written as ! ! [K]U = F (E.13) ! U5 where U = e f. The solution of Eq. (E.13) can be found as U6 U5 = 23.2922 * 10-6 in., U6 = - 139.7532 * 10-6 in. (E.14) The axial displacements of elements 1 and 2 can be found as U1 U2 u1 (1) t e f = [l(1)] d u2 U5 U6 0.8944 = c 0 = e - 0.4472 0 0 f in. 83.3301 * 10-6 0 0.8944 0 dµ - 0.4472 0 0 ∂ 23.2922 * 10-6 -6 - 139.7532 * 10 (E.15) RaoCh12ff.qxd 878 10.06.08 13:49 Page 878 CHAPTER 12 FINITE ELEMENT METHOD U3 u1 (2) U e f = [l(2)]µ 4 ∂ u2 U5 U6 0.8944 = c 0 = e 0.4472 0 0 0.8944 0 f in. - 41.6651 * 10-6 0 dµ 0.4472 0 0 ∂ 23.2922 * 10-6 - 139.7532 * 10-6 (E.16) The stresses in elements 1 and 2 can be determined as s(1) = E(1)P(1) = E(1) = (30 * 106)(83.3301 * 10-6) = 223.5989 psi 11.1803 s(2) = E(2)P(2) = = E(1)(u2 - u1)(1) ¢l(1) = l(1) l(1) (E.17) E(2)(u2 - u1)(2) E(2) ¢l(2) = (2) l l(2) (30 * 106)( - 41.6651 * 10-6) = - 111.7996 psi 11.1803 (E.18) where s(i) denotes the stress, e(i) represents the strain, and ¢l(i) indicates the change in length of element i (i = 1, 2). ■ 12.7 Consistent and Lumped Mass Matrices The mass matrices derived in Section 12.3 are called consistent mass matrices. They are consistent because the same displacement model that is used for deriving the element stiffness matrix is used for the derivation of mass matrix. It is of interest to note that several dynamic problems have been solved with simpler forms of mass matrices. The simplest form of the mass matrix, known as the lumped mass matrix, can be obtained by placing point (concentrated) masses mi at node points i in the directions of the assumed displacement degrees of freedom. The concentrated masses refer to translational and rotational inertia of the element and are calculated by assuming that the material within the mean locations on either side of the particular displacement behaves like a rigid body while the remainder of the element does not participate in the motion. Thus this assumption excludes the dynamic coupling that exists between the element displacements and hence the resulting element mass matrix is purely diagonal [12.11]. RaoCh12ff.qxd 10.06.08 13:49 Page 879 12.7 CONSISTENT AND LUMPED MASS MATRICES 879 12.7.1 Lumped Mass Matrix for a Bar Element By dividing the total mass of the element equally between the two nodes, the lumped mass matrix of a uniform bar element can be obtained as 12.7.2 Lumped Mass Matrix for a Beam Element In Fig. 12.4, by lumping one half of the total beam mass at each of the two nodes, along the translational degrees of freedom, we obtain the lumped mass matrix of the beam element as [m] = rAl 1 c 2 0 1 rAl 0 ≥ [m] = 2 0 0 0 0 0 0 0 d 1 0 0 1 0 (12.67) 0 0 ¥ 0 0 (12.68) Note that the inertia effect associated with the rotational degrees of freedom has been assumed to be zero in Eq. (12.68). If the inertia effect is to be included, we compute the mass moment of inertia of half of the beam segment about each end and include it at the diagonal locations corresponding to the rotational degrees of freedom. Thus, for a uniform beam, we have I = rAl3 1 rAl l 2 a ba b = 3 2 2 24 (12.69) and hence the lumped mass matrix of the beam element becomes 1 0 rAl [m] = F 2 0 0 12.7.3 Lumped Mass Versus Consistent Mass Matrices 0 l2 a b 12 0 0 0 0 0 1 0 0 0 l2 a b 12 V (12.70) It is not obvious whether the lumped mass matrices or consistent mass matrices yield more accurate results for a general dynamic response problem. The lumped mass matrices are approximate in the sense that they do not consider the dynamic coupling present between the various displacement degrees of freedom of the element. However, since the lumped mass matrices are diagonal, they require less storage space during computation. On the other hand, the consistent mass matrices are not diagonal and hence require more storage space. They too are approximate in the sense that the shape functions, which are derived using static displacement patterns, are used even for the solution of dynamics problems. RaoCh12ff.qxd 880 10.06.08 13:49 Page 880 CHAPTER 12 FINITE ELEMENT METHOD The following example illustrates the application of lumped and consistent mass matrices in a simple vibration problem. Consistent and Lumped Mass Matrices of a Bar EXAMPLE 12.4 Find the natural frequencies of the fixed-fixed uniform bar shown in Fig. 12.10 using consistent and lumped mass matrices. Use two bar elements for modeling. Solution: The stiffness and mass matrices of a bar element are AE 1 c l -1 [k] = [m]c = [m]l = rAl 2 c 6 1 rAl 1 c 2 0 -1 d 1 (E.1) 1 d 2 (E.2) 0 d 1 (E.3) where the subscripts c and l to the mass matrices denote the consistent and lumped matrices, respectively. Since the bar is modeled by two elements, the assembled stiffness and mass matrices are given by 1 1 AE [K ] = C 1 ' l 0 2 -1 1 +1 -1 1 2 rAl [M ] = C1 ' c 6 0 1 2 1 1 rAl [M ] = C0 ' l 2 0 0 1 U1 U2 Element 1 2 U3 l L FIGURE 12.10 Fixed-fixed uniform bar. -1 2 -1 +2 1 3 0 1 2 rAl 1S 2 = C1 6 2 3 0 1 4 1 0 1S 2 +1 0 3 0 1 1 rAl 0S 2 = C0 2 1 3 0 0 2 0 0 0S 1 2 Element 2 l 3 0 1 1 AE -1 S 2 = C -1 l 1 3 0 0 -1 S 1 (E.4) (E.5) (E.6) RaoCh12ff.qxd 10.06.08 13:49 Page 881 12.8 EXAMPLES USING MATLAB 881 The dashed boxes in Eqs. (E.4) through (E.6) enclose the contributions of elements 1 and 2. The degrees of freedom corresponding to the columns and rows of the matrices are indicated at the top and the right-hand side of the matrices. The eigenvalue problem, after applying the boundary conditions U1 = U3 = 0, becomes [[K] - v2[M]] 5U26 = 506 (E.7) ƒ [K] - v2[M] ƒ = 0 (E.8) The eigenvalue v2 can be determined by solving the equation which, for the present case, becomes ` rAl AE [2] - v2 [4] ` = 0 l 6 with consistent mass matrices (E.9) with lumped mass matrices (E.10) and ` rAl AE [2] - v2 [2] ` = 0 l 2 Equations (E.9) and (E.10) can be solved to obtain vc = 3E E = 3.4641 A rl2 A rL2 (E.11) vl = 2E E = 2.8284 A rl2 A rL2 (E.12) These values can be compared with the exact value (see Fig. 8.7) E A rL2 v1 = p (E.13) ■ 12.8 Examples Using MATLAB Finite Element Analysis of a Stepped Bar EXAMPLE 12.5 Consider the stepped bar shown in Fig. 12.11 with the following data: A1 = 16 * 10-4 m2, A 2 = 9 * 10-4 m2, A 3 = 4 * 10-4 m2, Ei = 20 * 1010 Pa, i = 1, 2, 3, ri = 7.8 * 103 kg/m3, i = 1, 2, 3, l 1 = 1 m, l 2 = 0.5 m, l 3 = 0.25 m. Write a MATLAB program to determine the following: a. b. Displacements u1, u2, and u3 under load p3 = 1000 N Natural frequencies and mode shapes of bar RaoCh12ff.qxd 882 10.06.08 13:49 Page 882 CHAPTER 12 FINITE ELEMENT METHOD A1, E1,  1 A2, E2,  2 u0 A3, E3,  3 u1 u3 u2 l1 l2 p3 l3 FIGURE 12.11 Stepped bar. Solution: The assembled stiffness and mass matrices of the stepped bar are given by [K ] = H ' A 1E1 l1 - A 1E1 l1 0 0 -A 1E1 l1 A 1E1 A 2E2 + l1 l2 - A 2E2 l2 0 2r1A 1l1 1 r1A 1l1 D [M ' ] = 6 0 0 r1A 1l1 2r1A 1l1 + 2r2A 2l2 r2A 2l2 0 0 0 -A 2E2 l2 A 3E3 A 2E2 + l2 l3 - A 3E3 l3 0 - A 3E3 l3 A 3E3 l3 0 r2A 2l2 2r2A 2l2 + 2r3A 3l3 r3A 3l3 X 0 0 T r3A 3l3 2r3A 3l3 (E.1) (E.2) The system matrices [K] and [M] can be obtained by incorporating the boundary condition u0 = 0— that is, by deleting the first row and first column in Eqs. (E.1) and (E.2). a. The equilibrium equations under the load p3 = 1000 N are given by ! ! [K]U = P (E.3) where A 2E2 A 1E1 + l1 l2 -A 2E2 [K] = G l2 0 - A 2E2 l2 A 3E3 A 2E2 + l2 l3 -A 3E3 l3 0 - A 3E3 W l3 A 3E3 l3 (E.4) RaoCh12ff.qxd 10.06.08 13:49 Page 883 12.8 EXAMPLES USING MATLAB u1 ! U = c u2 s , u3 b. 883 0 ! P = c0 s 1000 The eigenvalue problem can be expressed as ! ! c[K] - v2[M] dU = 0 (E.5) where [K] is given by Eq. (E.4) and [M] by 2r1A 1l1 + 2r2A 2l2 1 [M] = C r2A 2l2 6 0 r2A 2l2 2r2A 2l2 + 2r3A 3l3 r3A 3l3 0 r3A 3l3 S 2r3A 3l3 The MATLAB solution of Eqs. (E.3) and (E.5) is given below. %------ Program Ex12_5.m %------Initialization of values------------------------A1 = 16e4 ; A2 = 9e4 ; A3 = 4e4 ; E1 = 20e10 ; E2 = E1 ; E3 = E1 ; R1 = 7.8e3 ; R2 = R1 ; R3 = R1 ; L1 = 1 ; L2 = 0.5 ; L3 = 0.25 ; %------Definition of [K]--------------------------------K11 = A1*E1/L1+A2*E2/L2 ; K12 = A2*E2/L2 ; K13 = 0 ; K21 = K12 ; K22 = A2*E2/L2+A3*E3/L3 ; K23 = A3*E3/L3 ; K31 = K13 ; K32 = K23 ; K33 = A3*E3/L3 ; K = [ K11 K12 K13; K21 K22 K23; K31 K32 K33 ] %-------- Calculation of matrix P = [ 0 0 1000]’ U = inv(K)*P (E.6) RaoCh12ff.qxd 884 10.06.08 13:49 Page 884 CHAPTER 12 FINITE ELEMENT METHOD %------- Definition of [M] ------------------------M11 = (2*R1*A1*L1+2*R2*A2*L2) / 6; M12 = (R2*A2*L2) / 6; M13 = 0; M21 = M12; M22 = (2*R2*A2*L2+2*R3*A3*L3) / 6; M23 = R3*A3*L3; M31 = M13; M32 = M23; M33 = 2*M23; M= [M11 M12 M13; M21 M22 M23; M31 M32 M33 ] MI = inv (M) KM = MI*K %-------------Calculation of eigenvector and eigenvalue-------------[L, V] = eig (KM) >> Ex12_5 K = 360000000 680000000 320000000 680000000 360000000 0 0 320000000 320000000 P = 0 0 1000 U = 1.0e005 * 0.3125 0.5903 0.9028 M = 5.3300 0.5850 0 0.5850 1.4300 0.7800 0 0.7800 1.5600 0.2000 0.1125 0.0562 0.1125 1.0248 0.5124 0.0562 0.5124 0.8972 MI = KM = 1.0e+008 * 1.7647 4.4542 2.2271 1.6647 9.0133 6.5579 0.5399 4.9191 4.5108 0.1384 0.7858 0.6028 0.6016 0.1561 0.7834 0.3946 0.5929 0.7020 L = RaoCh12ff.qxd 10.06.08 13:49 Page 885 12.8 EXAMPLES USING MATLAB 885 V = 1.0e+009 * 1.3571 0 0 0 0.1494 0 0 0 0.0224 >> ■ Program for Eigenvalue Analysis of a Stepped Beam EXAMPLE 12.6 Develop a MATLAB program called Program17.m for the eigenvalue analysis of a fixed-fixed stepped beam of the type shown in Fig. 12.12. Solution: Program17.m is developed to accept the following input data: xl(i) = length of element (step) i xi(i) = moment of inertia of element i a(i) = area of cross section of element i bj(i, j) = global degree of freedom number corresponding to the local jth degree of freedom of element i e = Young’s modulus rho = mass density The program gives the natural frequencies and mode shapes of the beam as output. Natural frequencies of the stepped beams 1.6008e+002 6.1746e+002 2.2520e+003 7.1266e+003 Mode shapes 1 1.0333e002 1.8915e004 1.4163e002 4.4518e005 2 3.7660e003 2.0297e004 4.7109e003 2.5950e004 3 1.6816e004 1.8168e004 1.3570e003 2.0758e004 4 1.8324e004 6.0740e005 3.7453e004 1.6386e004 W1 W2 W3 W4 W5 W7 W6 W8 2  2 3  3 1  1 l1  40 l2  32 l3  24 x, X E  30  106 psi,   0.283 lb/in3 FIGURE 12.12 Stepped beam. ■ RaoCh12ff.qxd 886 10.06.08 13:49 Page 886 CHAPTER 12 FINITE ELEMENT METHOD 12.9 C++ Program An interactive C++ program, called Program17.cpp, is given for the eigenvalue solution of a stepped beam. The input and output of the program are similar to those of Program17.m. Eigenvalue Solution of a Stepped Beam EXAMPLE 12.7 Find the natural frequencies and mode shapes of the stepped beam shown in Fig. 12.12 using Program17.cpp. Solution: The input data are to be entered interactively. The output of the program is shown below. NATURAL FREQUENCIES OF THE STEPPED BEAM 160.083080 617.459700 2251.975785 7126.595358 MODE SHAPES 1 2 3 4 0.0103332469 0.0037659939 0.0001681571 0.0001832390 0.0001891485 0.0002029733 0.0001816828 0.0000607403 0.0141625494 0.0047108916 0.0013569926 0.0003745272 0.0000445137 0.0002594971 0.0002075837 0.0001638648 ■ 12.10 Fortran Program A Fortran program called PROGRAM17.F is given for the eigenvalue analysis of a stepped beam. The input and output of the program are similar to those of Program17.m. Eigenvalue Analysis of a Stepped Beam EXAMPLE 12.8 Find the natural frequencies and mode shapes of the stepped beam shown in Fig. 12.12 using PROGRAM17.F. Solution: The output of the program is given below. NATURAL FREQUENCIES OF THE STEPPED BEAM 0.160083E+03 0.617460E+03 MODE SHAPES 1 0.103333E01 2 0.376593E02 3 0.167902E03 4 0.182648E03 0.225198E+04 0.189147E03 0.202976E03 0.181687E03 0.607247E04 0.712653E+04 0.141626E01 0.471097E02 0.135672E02 0.373775E03 0.445258E04 0.259495E03 0.207580E03 0.163869E03 ■ RaoCh12ff.qxd 10.06.08 13:49 Page 887 REVIEW QUESTIONS 887 REFERENCES 12.1 O. C. Zienkiewicz, The Finite Element Method (4th ed.), McGraw-Hill, London, 1987. 12.2 S. S. Rao, The Finite Element Method in Engineering (3rd ed.), Butterworth-Heinemann, Boston, 1999. 12.3 G. V. Ramana and S. S. Rao, “Optimum design of plano-milling machine structure using finite element analysis,” Computers and Structures, Vol. 18, 1984, pp. 247–253. 12.4 R. Davis, R. D. Henshell, and G. B. Warburton, “A Timoshenko beam element,” Journal of Sound and Vibration, Vol. 22, 1972, pp. 475–487. 12.5 D. L. Thomas, J. M. Wilson, and R. R. Wilson, “Timoshenko beam finite elements,” Journal of Sound and Vibration, Vol. 31, 1973, pp. 315–330. 12.6 J. Thomas and B. A. H. Abbas, “Finite element model for dynamic analysis of Timoshenko beams,” Journal of Sound and Vibration, Vol. 41, 1975, pp. 291–299. 12.7 R. S. Gupta and S. S. Rao, “Finite element eigenvalue analysis of tapered and twisted Timoshenko beams,” Journal of Sound and Vibration, Vol. 56, 1978, pp. 187–200. 12.8 T. H. H. Pian, “Derivation of element stiffness matrices by assumed stress distribution,” AIAA Journal, Vol. 2, 1964, pp. 1333–1336. 12.9 H. Alaylioglu and R. Ali, “Analysis of an automotive structure using hybrid stress finite elements,” Computers and Structures, Vol. 8, 1978, pp. 237–242. 12.10 I. Fried, “Accuracy of finite element eigenproblems,” Journal of Sound and Vibration, Vol. 18, 1971, pp. 289–295. 12.11 P. Tong, T. H. H. Pian, and L. L. Bucciarelli, “Mode shapes and frequencies by the finite element method using consistent and lumped matrices,” Computers and Structures, Vol. 1, 1971, pp. 623–638. REVIEW QUESTIONS 12.1 Give brief answers to the following: 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. What is the basic idea behind the finite element method? What is a shape function? What is the role of transformation matrices in the finite element method? What is the basis for the derivation of transformation matrices? How are fixed boundary conditions incorporated in the finite element equations? How do you solve a finite element problem having symmetry in geometry and loading by modeling only half of the problem? Why is the finite element approach presented in this chapter called the displacement method? What is a consistent mass matrix? What is a lumped mass matrix? What is the difference between the finite element method and the Rayleigh-Ritz method? How is the distributed load converted into equivalent joint force vector in the finite element method? 12.2 Indicate whether each of the following statements is true or false: 1. For a bar element of length l with two nodes, the shape function corresponding to node 2 is given by x/l. RaoCh12ff.qxd 888 10.06.08 13:49 Page 888 CHAPTER 12 FINITE ELEMENT METHOD 2. The element stiffness matrices are always singular. 3. The element mass matrices are always singular. 4. The system stiffness matrix is always singular unless the boundary conditions are incorporated. 5. The system mass matrix is always singular unless the boundary conditions are incorporated. 6. The lumped mass matrices are always diagonal. 7. The coordinate transformation of element matrices is required for all systems. 8. The element stiffness matrix in the global coordinate system, [k], can be expressed in terms of the local matrix [k] and the coordinate transformation matrix [l] as [l]T[k][l]. 9. The derivation of system matrices involves the assembly of element matrices. 10. Boundary conditions are to be imposed to avoid rigid body motion of the system. 12.3 Fill in each of the following blanks with appropriate word: 1. In the finite element method, the solution domain is replaced by several _____. 2. In the finite element method, the elements are assumed to be interconnected at certain points known as _____. 3. In the finite element method, an _____ solution is assumed within each element. 4. The displacement within a finite element is expressed in terms of _____ functions. 5. For a thin beam element, _____ degrees of freedom are considered at each node. 6. For a thin beam element, the shape functions are assumed to be polynomials of degree _____. 7. In the displacement method, the _____ of elements is directly approximated. 8. If the displacement model used in the derivation of the element stiffness matrices is also used to derive the element mass matrices, the resulting mass matrix is called _____ mass matrix. 9. If the mass matrix is derived by assuming point masses at node points, the resulting mass matrix is called _____ mass. 10. The lumped mass matrices do not consider the _____ coupling between the various displacement degrees of freedom of the element. 11. Different orientations of finite elements require _____ of element matrices. 12.4 Select the most appropriate answer out of the choices given: 1. For a bar element of length l with two nodes, the shape function corresponding to node 1 is given by x x x (b) (c) a 1 + b b l l l 2. The simplest form of mass matrix is known as (a) lumped mass matrix (b) consistent mass matrix (c) global mass matrix 3. The finite element method is (a) an approximate analytical method (b) a numerical method (c) an exact analytical method (a) a 1 - RaoCh12ff.qxd 10.06.08 13:49 Page 889 PROBLEMS 889 4. The stiffness matrix of a bar element is given by EA 1 1 EA 1 -1 EA 1 0 c d c d c d (b) (c) l 1 1 l -1 1 l 0 1 5. The consistent mass matrix of a bar element is given by rAl 1 0 rAl 2 1 rAl 2 -1 c d c d c d (a) (b) (c) 6 1 2 6 -1 2 6 0 1 6. The finite element method is similar to (a) Rayleigh’s method (b) the Rayleigh-Ritz method (c) the Lagrange method 7. The lumped mass matrix of a bar element is given by rAl 1 rAl 2 1 1 0 d c d c (a) rAl c (b) (c) 0 1 6 1 2 2 0 (a) 0 d 1 8. The element mass matrix in the global coordinate system, 3m4, can be expressed in terms of the element mass matrix in local coordinate system [m] and the coordinate transformation matrix [l], as (a) 3m4 = [l]T[m] (b) 3m4 = [m][l] (c) 3m4 = [l]T[m][l] 12.5 Match the items in the two columns below. Assume a fixed-fixed bar with one middle node: rAl 2 1 rAl 1 0 AE 1 -1 c d, [m]c = c d, [m]l = c d Element matrices: [k] = l -1 1 6 1 2 2 0 1 Steel bar: E = 30 * 106 lb/in.2, r = 0.0007298 lb - sec2/in.4, L = 12 in. Aluminum bar: E = 10.3 * 106 lb/in.2, r = 0.0002536 lb - sec2/in.4, L = 12 in. 1. Natural frequency of steel bar given by lumped mass matrices 2. Natural frequency of aluminum bar given by consistent mass matrices 3. Natural frequency of steel bar given by consistent mass matrices 4. Natural frequency of aluminum bar given by lumped mass matrices (a) 58,528.5606 rad/sec (b) 47,501.0898 rad/sec (c) 58,177.2469 rad/sec (d) 47,787.9336 rad/sec PROBLEMS The problem assignments are organized as follows: Problems Section Covered 12.1, 12.2, 12.4 12.3 12.5, 12.7 12.6, 12.9 12.4, 12.8, 12.10–12.30 12.4 12.5 12.6 Topic Covered Derivation of element matrices and vectors Transformation matrix Assembly of matrices and vectors Application of boundary conditions and solution of problem RaoCh12ff.qxd 890 10.06.08 13:49 Page 890 CHAPTER 12 FINITE ELEMENT METHOD Problems Section Covered 12.31, 12.32 12.7 12.33–12.35 12.36 12.3, 12.37–12.40 12.41–12.42 12.8 12.9 12.10 — Topic Covered Consistent and lumped mass matrices MATLAB programs C++ program Fortran programs Design projects 12.1 Derive the stiffness matrix and the consistent and lumped mass matrices of the tapered bar element shown in Fig. 12.13. The diameter of the bar decreases from D to d over its length. , E D d x l FIGURE 12.13 12.2 Derive the stiffness matrix of the bar element in longitudinal vibration whose cross-sectional area varies as A(x) = A 0e - (x/l), where A 0 is the area at the root (see Fig. 12.14). A0 e(x/l) O x x l FIGURE 12.14 12.3 Write a computer program for finding the stresses in a planar truss. 12.4 The tapered cantilever beam shown in Fig. 12.15 is used as a spring to carry a load P. (a) Derive the stiffness matrix of the beam element. (b) Use the result of (a) to find the stress induced in the beam when B = 25 cm, b = 10 cm, t = 2.5 cm, l = 2 m, E = 2.07 * 1011 N/m2, and P = 1000 N. Use one beam element for idealization. 12.5 Find the global stiffness matrix of each of the four bar elements of the truss shown in Fig. 12.5 using the following data: Nodal coordinates: (X1, Y1) = (0, 0), (X2, Y2) = (50, 100) in., (X3, Y3) = (100, 0) in., (X4, Y4) = (200, 150) in. Cross-sectional areas: A 1 = A 2 = A 3 = A 4 = 2 in.2. Young’s modulus of all members: 30 * 106 lb/in.2. RaoCh12ff.qxd 10.06.08 13:49 Page 891 PROBLEMS 891 z B/2 P B/2 O t b/2 l x b/2 FIGURE 12.15 12.6 Using the result of Problem 12.5, find the assembled stiffness matrix of the truss and formulate the equilibrium equations if the vertical downward load applied at node 4 is 1000 lb. 12.7 Derive the stiffness and mass matrices of the planar frame element (general beam element) shown in Fig. 12.16 in the global XY-coordinate system. Y U5 U6 Joint 2 U4 U2 U3 U1 Joint 1 O X FIGURE 12.16 A frame element in global system. 12.8 A multiple-leaf spring used in automobiles is shown in Fig. 12.17. It consists of five leaves, each of thickness t = 0.25 in. and width w = 1.5 in. Find the deflection of the leaves under a load of P = 2000 lb. Model only a half of the spring for the finite element analysis. The Young’s modulus of the material is 30 * 106 psi. RaoCh12ff.qxd 892 10.06.08 13:49 Page 892 CHAPTER 12 FINITE ELEMENT METHOD P P Eye t Center bolt w 5 5 5 5 5 FIGURE 12.17 A multiple-leaf spring. 12.9 Derive the assembled stiffness and mass matrices of the multiple-leaf spring of Problem 12.8 assuming a specific weight of 0.283 lb/in.3 for the material. 12.10 Find the nodal displacements of the crane shown in Fig. 12.18 when a vertically downward load of 1000 lb is applied at node 4. The Young’s modulus is 30 * 106 psi and the crosssectional area is 2 in.2 for elements 1 and 2 and 1 in.2 for elements 3 and 4. Y 2 4 2 1000 lb 3 25 1 1 3 1 1 2 1 2 X 50 100 4 50 2 1 2 100 FIGURE 12.18 12.11 Find the tip deflection of the cantilever beam shown in Fig. 12.19 when a vertical load of P = 500 N is applied at point Q using (a) a one element approximation and (b) a two-element approximation. Assume l = 0.25 m, h = 25 mm, b = 50 mm, E = 2.07 * 1011 Pa, and k = 105 N/m. 12.12 Find the stresses in the stepped beam shown in Fig. 12.20 when a moment of 1000 N-m is applied at node 2 using a two-element idealization. The beam has a square cross section 50 * 50 mm between nodes 1 and 2 and 25 * 25 mm between nodes 2 and 3. Assume the Young’s modulus as 2.1 * 1011 Pa. RaoCh12ff.qxd 10.06.08 13:49 Page 893 PROBLEMS , A, I, E X 893 h Q X b k Section X – X l FIGURE 12.19 1000 N-m 1 2 3 0.40 m 0.25 m FIGURE 12.20 12.13 Find the transverse deflection and slope of node 2 of the beam shown in Fig. 12.21 using a two element idealization. Compare the solution with that of simple beam theory. P  , A, I, E 1 2 3 3l 4 l 4 FIGURE 12.21 12.14 Find the natural frequencies of a cantilever beam of length l, cross-sectional area A, moment of inertia I, Young’s modulus E, and density r, using one finite element. 12.15 Using one beam element, find the natural frequencies of the uniform pinned-free beam shown in Fig. 12.22. , A, I, E l FIGURE 12.22 RaoCh12ff.qxd 894 10.06.08 13:49 Page 894 CHAPTER 12 FINITE ELEMENT METHOD 12.16 Using one beam element and one spring element, find the natural frequencies of the uniform, spring-supported cantilever beam shown in Fig. 12.19. 12.17 Using one beam element and one spring element, find the natural frequencies of the system shown in Fig. 12.23. , A, I, E k m   Al El l3 m l FIGURE 12.23 12.18 Using two beam elements, find the natural frequencies and mode shapes of the uniform fixedfixed beam shown in Fig. 12.24. , A, I, E L 2 L 2 FIGURE 12.24 12.19* An electric motor, of mass m = 100 kg and operating speed = 1800 rpm, is fixed at the middle of a clamped-clamped steel beam of rectangular cross section, as shown in Fig. 12.25. Design the beam such that the natural frequency of the system exceeds the operating speed of the motor. 2m FIGURE 12.25 *An asterisk denotes a problem with no unique answer. RaoCh12ff.qxd 10.06.08 13:49 Page 895 PROBLEMS 895 12.20 Find the natural frequencies of the beam shown in Fig. 12.26, using three finite elements of length l each. , A, I, E l l l FIGURE 12.26 12.21 Find the natural frequencies of the cantilever beam carrying an end mass M shown in Fig. 12.27, using a one beam element idealization. , A, I, E M M  10 Al l FIGURE 12.27 12.22 Find the natural frequencies of vibration of the beam shown in Fig. 12.28, using two beam elements. Also find the load vector if a uniformly distributed transverse load p is applied to element 1. y, Y , E, A1, I1 , E, A2, I2 x, X l1 l2 FIGURE 12.28 12.23 Find the natural frequencies of a beam of length l, which is pin connected at x = 0 and fixed at x = l, using one beam element. RaoCh12ff.qxd 896 10.06.08 13:49 Page 896 CHAPTER 12 FINITE ELEMENT METHOD 12.24 Find the natural frequencies of torsional vibration of the stepped shaft shown in Fig. 12.29. Assume that r1 = r2 = r, G1 = G2 = G, Ip1 = 2Ip2 = 2Ip, J1 = 2J2 = 2J, and l1 = l2 = l. 1, G1, Ip1, J1 2, G2, Ip2, J2 l1 l2 FIGURE 12.29 12.25 Find the dynamic response of the stepped bar shown in Fig. 12.30(a) when its free end is subjected to the load given in Fig. 12.30(b). P(t) Area  4A Area  A P0 P(t) l l O (a) t0 t (b) FIGURE 12.30 12.26 Find the displacement of node 3 and the stresses in the two members of the truss shown in Fig. 12.31. Assume that the Young’s modulus and the cross-sectional areas of the two members are the same with E = 30 * 106 psi and A = 1 in.2. 2 10 in. 3 1 1 lb 25 in. FIGURE 12.31 RaoCh12ff.qxd 10.06.08 13:49 Page 897 PROBLEMS 897 12.27 The simplified model of a radial drilling machine structure is shown in Fig. 12.32. Using two beam elements for the column and one beam element for the arm, find the natural frequencies and mode shapes of the machine. Assume the material of the structure as steel. z 400 mm 2.4 m 0.4 m 415 mm Cross section of column Arm 2m A Column 550 mm 15 mm x 350 mm Cross section of arm FIGURE 12.32 A radial drilling machine structure. 12.28 If a vertical force of 5000 N along the z-direction and a bending moment of 500 N-m in the xz-plane are developed at point A during a metal cutting operation, find the stresses developed in the machine tool structure shown in Fig. 12.32. 12.29 The crank in the slider-crank mechanism shown in Fig. 12.33 rotates at a constant clockwise angular speed of 1000 rpm. Find the stresses in the connecting rod and the crank when the X X X X p 0.5 in. 3 in. 2 in. 0.5 in. Section X – X FIGURE 12.33 A slider-crank mechanism.  RaoCh12ff.qxd 898 10.06.08 13:49 Page 898 CHAPTER 12 FINITE ELEMENT METHOD pressure acting on the piston is 200 psi and u = 30°. The diameter of the piston is 12 in. and the material of the mechanism is steel. Model the connecting rod and the crank by one beam element each. The lengths of the crank and connecting rod are 12 in. and 48 in., respectively. 12.30 A water tank of weight W is supported by a hollow circular steel column of inner diameter d, wall thickness t, and height l. The wind pressure acting on the column can be assumed to vary linearly from 0 to pmax, as shown in Fig. 12.34. Find (a) the bending stress induced in the column under the loads, and (b) the natural frequencies of the water tank using a one beam element idealization. Data: W = 10,000 lb, l = 40 ft, d = 2 ft, t = 1 in., and pmax = 100 psi. W Water tank pmax Column l FIGURE 12.34 12.31 Find the natural frequencies of the stepped bar shown in Fig. 12.35 with the following data using consistent and lumped mass matrices: A 1 = 2 in.2, A 2 = 1 in.2, E = 30 * 106 psi, r = 0.283 lb/in.3, and l1 = l2 = 50 in. A1 U1 A2 U2 l1 U3 l2 FIGURE 12.35 12.32 Find the undamped natural frequencies of longitudinal vibration of the stepped bar shown in Fig. 12.36 with the following data using consistent and lumped mass matrices: l1 = l2 = l3 = 0.2 m, A 1 = 2A 2 = 4A 3 = 0.4 * 10-3 m2, E = 2.1 * 1011 N/m2, and r = 7.8 * 103 kg/m3. RaoCh12ff.qxd 10.06.08 13:49 Page 899 DESIGN PROJECTS A1 A2 899 A3 x, X O l1 l2 l3 FIGURE 12.36 12.33 Consider the stepped bar shown in Fig. 12.11 with the following data: A1 = 25 * 10-4 m2, A2 = 16 * 10-4 m2, A3 = 9 * 10-4 m2, Ei = 2 * 1011 Pa, i = 1, 2, 3, ri = 7.8 * 103 kg/m3, i = 1, 2, 3, l1 = 3 m, l2 = 2 m, l3 = 1 m. Using MATLAB, find the axial displacements u1, u2, and u3 under the load p3 = 500 N. 12.34 Using MATLAB, find the natural frequencies and mode shapes of the stepped bar described in Problem 12.33. 12.35 Use Program17.m to find the natural frequencies of a fixed-fixed stepped beam, similar to the one shown in Fig. 12.12, with the following data: Cross sections of elements: 1, 2, 3: 4– * 4–, 3– * 3–, 2– * 2– Lengths of elements: 1, 2, 3: 30–, 20–, 10– Young’s modulus of all elements: 107 lb/in.2 Weight density of all elements: 0.1lb/in.3 12.36 Use Program17.cpp to solve Problem 12.35. 12.37 Use PROGRAM17.F to solve Problem 12.35. 12.38 Write a computer program for finding the assembled stiffness matrix of a general planar truss. 12.39 Generalize the computer program of Section 12.10 to make it applicable to the solution of any stepped beam having a specified number of steps. 12.40 Find the natural frequencies and mode shapes of the beam shown in Fig. 12.12 with l1 = l2 = l3 = 10 in. and a uniform cross section of 1 in. * 1 in. throughout the length, using the computer program of Section 12.10. Compare your results with those given in Chapter 8. (Hint: Only the data XL, XI, and A need to be changed.) DESIGN PROJECTS 12.41 Derive the stiffness and mass matrices of a uniform beam element in transverse vibration rotating at an angular velocity of Æ rad/sec about a vertical axis as shown in Fig. 12.37(a). Using these matrices, find the natural frequencies of transverse vibration of the rotor blade of a helicopter (see Fig. 12.37b) rotating at a speed of 300 rpm. Assume a uniform rectangular cross section 1– * 12– and a length 48– for the blade. The material of the blade is aluminum. 12.42 An electric motor weighing 1000 lb operates on the first floor of a building frame that can be modeled by a steel girder supported by two reinforced concrete columns, as shown in Fig. 12.38. If the operating speed of the motor is 1500 rpm, design the girder and the columns such that the fundamental frequency of vibration of the building frame is greater than the operating speed RaoCh12ff.qxd 900 10.06.08 13:49 Page 900 CHAPTER 12 FINITE ELEMENT METHOD R l x O Beam element. , A, I, E (a) Rotor blade (b) FIGURE 12.37 Motor h Girder b Cross section of girder 9 ft Columns d Cross section of columns 18 ft 9 ft FIGURE 12.38 of the motor. Use two beam and two bar elements for the idealization. Assume the following data: Girder: E = 30 * 106 psi, Columns: E = 4 * 106 psi, r = 8.8 * 10-3 lbm/in.3, r = 2.7 * 10-3 lbm/in.3 h/b = 2 RaoCh13ff.qxd 10.06.08 13:57 Page 901 Jules Henri Poincaré (1854–1912) was a French mathematician and professor of celestial mechanics at the University of Paris and of mechanics at the Ecole Polytechnique. His contributions to pure and applied mathematics, particularly to celestial mechanics and electrodynamics, are outstanding. His classification of singular points of nonlinear autonomous systems is important in the study of nonlinear vibrations. (Photo courtesy of Dirk J. Struik, A Concise History of Mathematics, 2nd ed. Dover Publications, New York, 1948) C H A P T E R 1 3 Nonlinear Vibration 13.1 Introduction In the preceding chapters, the equation of motion contained displacement or its derivatives only to the first degree, and no square or higher powers of displacement or velocity were involved. For this reason, the governing differential equations of motion and the corresponding systems were called linear. For convenience of analysis, most systems are modeled as linear systems, but real systems are actually more often nonlinear than linear [13.1–13.6]. Whenever finite amplitudes of motion are encountered, nonlinear analysis becomes necessary. The superposition principle, which is very useful in linear analysis, does not hold true in the case of nonlinear analysis. Since mass, damper, and spring are the basic components of a vibratory system, nonlinearity into the governing differential equation may be introduced through any of these components. In many cases, linear analysis is insufficient to describe the behavior of the physical system adequately. One of the main reasons for modeling a physical system as a nonlinear one is that totally unexpected phenomena sometimes occur in nonlinear systems—phenomena that are not predicted or even hinted at by linear theory. Several methods are available for the solution of nonlinear vibration problems. Some of the exact methods, approximate analytical techniques, graphical procedures, and numerical methods are presented in this chapter. 901 RaoCh13ff.qxd 902 13.2 10.06.08 13:57 Page 902 CHAPTER 13 NONLINEAR VIBRATION Examples of Nonlinear Vibration Problems The following examples are given to illustrate the nature of nonlinearity in some physical systems. 13.2.1 Simple Pendulum Consider a simple pendulum of length l, having a bob of mass m, as shown in Fig. 13.1(a). The differential equation governing the free vibration of the pendulum can be derived from Fig. 13.1(b): $ ml2 u + mgl sin u = 0 (13.1) For small angles, sin u may be approximated by u and Eq. (13.1) reduces to a linear equation: $ (13.2) u + v20 u = 0 v0 = 1g/l21/2 (13.3) u1t2 = A 0sin1v0t + f2 (13.4) where The solution of Eq. (13.2) can be expressed as where A 0 is the amplitude of oscillation, f is the phase angle, and v0 is the angular frequency. The values of A 0 and f are determined by the initial conditions and the angular frequency v0 is independent of the amplitude A 0. Equation (13.4) denotes an approximate solution of the simple pendulum. A better approximate solution can be obtained by using a two-term approximation for sin u near u = 0 as u - u3/6 in Eq. (13.1) $ u3 ml2 u + mgl au - b = 0 6 $ u + v201u - 61 u32 = 0 or O  l ..  ,  T Inertia moment .. ml 2 m mg mg cos   mg sin  (a) FIGURE 13.1 (b) Simple pendulum. (13.5) RaoCh13ff.qxd 10.06.08 13:57 Page 903 13.2 EXAMPLES OF NONLINEAR VIBRATION PROBLEMS f (x) f (x) x O (a) Soft spring FIGURE 13.2 903 O x (b) Hard spring Nonlinear spring characteristics. It can be seen that Eq. (13.5) is nonlinear because of the term involving u3 (due to geometric nonlinearity). Equation (13.5) is similar to the equation of motion of a spring-mass system with a nonlinear spring. If the spring is nonlinear (due to material nonlinearity), the restoring force can be expressed as f(x), where x is the deformation of the spring and the equation of motion of the spring-mass system becomes $ (13.6) mx + f1x2 = 0 If df/dx1x2 = k = constant, the spring is linear. If df/dx is a strictly increasing function of x, the spring is called a hard spring, and if df/dx is a strictly decreasing function of x, the spring is called a soft spring as shown in Fig. 13.2. Due to the similarity of Eqs. (13.5) and (13.6), a pendulum with large amplitudes is considered, in a loose sense, as a system with a nonlinear elastic (spring) component. 13.2.2 Mechanical Chatter, Belt Friction System Nonlinearity may be reflected in the damping term as in the case of Fig. 13.3(a). The system behaves nonlinearly because of the dry friction between the mass m and the moving belt. For this system, there are two friction coefficients: the static coefficient of friction 1ms2, corresponding to the force required to initiate the motion of the body held by dry friction; and the kinetic coefficient of friction 1mk2, corresponding to the force required to maintain the body in motion. In either case, the component of the applied force tangent to the friction surface (F) is the product of the appropriate friction coefficient and the force normal to the surface. The sequence of motion of the system shown in Fig. 13.3(a) is as follows [13.7]. The mass is initially at rest on the belt. Due to the displacement of the mass m along with the belt, the spring elongates. As the spring extends, the spring force on the mass increases until the static friction force is overcome and the mass begins to slide. It slides rapidly towards the right, thereby relieving the spring force until the kinetic friction force halts it. The spring then begins to build up the spring force again. The variation of the damping RaoCh13ff.qxd 904 10.06.08 13:57 Page 904 CHAPTER 13 NONLINEAR VIBRATION x Dry friction  Friction force, F (x.) k m v Belt Roller  Roller v O (a) FIGURE 13.3 Velocity of mass, m(x.) (b) Dry friction damping. force with the velocity of the mass is shown in Fig. 13.3(b). The equation of motion can be expressed as $ # mx + F1x2 + kx = 0 (13.7) # where the friction force F is a nonlinear function of x, as shown in Fig. 13.3(b). # For large values of x, the damping force is positive (the curve has a positive slope) and # energy is removed from the system. On the other hand, for small values of x, the damping force is negative (the curve has a negative slope) and energy is put into the system. Although there is no external stimulus, the system can have an oscillatory motion; it corresponds to a nonlinear self-excited system. This phenomenon of self-excited vibration is called mechanical chatter. 3.2.3 Variable Mass System Nonlinearity may appear in the mass term as in the case of Fig. 13.4 [13.8]. For large deflections, the mass of the system depends on the displacement x, and so the equation of motion becomes d # 1mx2 + kx = 0 dt (13.8) Note that this is a nonlinear differential equation due to the nonlinearity of the first term. FIGURE 13.4 Variable mass system. RaoCh13ff.qxd 10.06.08 13:57 Page 905 905 13.3 EXACT METHODS 13.3 Exact Methods An exact solution is possible only for a relatively few nonlinear systems whose motion is governed by specific types of second-order nonlinear differential equations. The solutions are exact in the sense that they are given either in closed form or in the form of an expression that can be numerically evaluated to any degree of accuracy. In this section, we shall consider a simple nonlinear system for which the exact solution is available. For a single degree of freedom system with a general restoring (spring) force F(x), the free vibration equation can be expressed as $ x + a2F1x2 = 0 (13.9) where a2 is a constant. Equation (13.9) can be rewritten as d #2 1x 2 + 2a2F1x2 = 0 dx (13.10) Assuming the initial displacement as x0 and the velocity as zero at t = t0, Eq. (13.10) can be integrated to obtain # x2 = 2a2 Lx # ƒ x ƒ = 12a e x0 F1h2dh or Lx x0 F1h2 d h f 1/2 (13.11) where h is the integration variable. Equation (13.11), when integrated again, gives 1 t - t0 = 12a L0 x e F1h2 dh f dj Lj x0 1/2 (13.12) where j is the new integration variable and t0 corresponds to the time when x = 0. Equation (13.12) thus gives the exact solution of Eq. (13.9) in all those situations where the integrals of Eq. (13.12) can be evaluated in closed form. After evaluating the integrals of Eq. (13.12), one can invert the result and obtain the displacement-time relation. If F(x) is an odd function, F1-x2 = - F 1x2 (13.13) By considering Eq. (13.12) from zero displacement to maximum displacement, the period of vibration t can be obtained: t = 4 12a L0 x0 e F1h2 dh f dj Lj x0 1/2 (13.14) For illustration, let F1x2 = xn. In this case Eqs. (13.12) and (13.14) become 0 dj 1 n + 1 a A 2 L0 1xn0 + 1 - jn + 121/2 x t - t0 = (13.15) RaoCh13ff.qxd 906 10.06.08 13:57 Page 906 CHAPTER 13 NONLINEAR VIBRATION and 0 dj 4 n + 1 a A 2 L0 1xn0 + 1 - jn + 121/2 By setting y = j/x0, Eq. (13.16) can be written as x t = dy 4 n + 1 1 t = a 1xn0 - 121/2 A 2 L0 11 - yn + 121/2 (13.16) 1 (13.17) This expression can be evaluated numerically to any desired level of accuracy. 13.4 Approximate Analytical Methods In the absence of an exact analytical solution to a nonlinear vibration problem, we wish to find at least an approximate solution. Although both analytical and numerical methods are available for approximate solution of nonlinear vibration problems, the analytical methods are more desirable [13.6, 13.9]. The reason is that once the analytical solution is obtained, any desired numerical values can be substituted and the entire possible range of solutions can be investigated. We shall now consider four analytical techniques in the following subsections. 13.4.1 Basic Philosophy Let the equations governing the vibration of a nonlinear system be represented by a system of n first order differential equations1 # : ! ! ! : (13.18) x 1t2 = f 1x, t2 + ag1x, t2 ! ! where the nonlinear terms are assumed to appear only in g1x, t2 and a is a small parameter. In Eq. (13.18) x1 x2 . ! x = f . v, . xn and 1 dx1/dt dx2/dt . # : x = f . v, . f11x1, x2, Á , xn, t2 f21x1, x2, Á , xn, t2 . : ! f 1x, t2 = f v . . dxn/dt fn1x1, x2, Á , xn, t2 g11x1, x2, Á , xn, t2 g21x1, x2, Á , xn, t2 . ! ! g1x, t2 = f v . . gn1x1, x2, Á , xn, t2 Systems governed by Eq. (13.18), in which the time appears explicitly, are known as nonautonomous systems. On the other hand, systems for which the governing equations are of the type #! : ! x 1t2 = f1x2 + ag1x2 where time does not appear explicitly are called autonomous systems. RaoCh13ff.qxd 10.06.08 13:57 Page 907 907 13.4 APPROXIMATE ANALYTICAL METHODS The solution of differential equations having nonlinear terms associated with a small parameter was studied by Poincaré [13.6]. Basically, he assumed the solution of Eq. (13.18) in series form as ! ! ! ! ! x1t2 = x01t2 + ax11t2 + a2 x21t2 + a3 x31t2 + Á (13.19) The series solution of Eq. (13.19) has two basic characteristics: 1. As a : 0, Eq. (13.19) reduces to the exact solution of the linear equations #! : ! x = f 1x, t2. 2. For small values of a, the series converges fast so that even the first two or three terms in the series of Eq. (13.19) yields a reasonably accurate solution. The various approximate analytical methods presented in this section can be considered to be modifications of the basic idea contained in Eq. (13.19). Although Poincaré’s solution, Eq. (13.19), is valid for only small values of a, the method can still be applied to systems with large values of a. The solution of the pendulum equation, Eq. (13.5), is presented to illustrate the Poincaré’s method. Solution of Pendulum Equations. Equation (13.5) can be rewritten as $ x + v20x + ax3 = 0 where x = u, v0 = 1g/l21/2, and a = - v20/6. Equation (13.20) is known as the free Duffing’s equation. Assuming weak nonlinearity (i.e., a is small), the solution of Eq. (13.20) is expressed as x1t2 = x01t2 + ax11t2 + a2x21t2 + Á + anxn1t2 + Á (13.20) where xi1t2, i = 0, 1, 2, Á , n, are functions to be determined. By using a two-term approximation in Eq. (13.21), Eq. (13.20) can be written as that is, (13.21) $ $ 1x 0 + ax12 + v201x0 + ax12 + a1x0 + ax123 = 0 $ $ 1x0 + v20 x02 + a1x1 + v20 x1 + x302 + a213x20 x12 + a313x0 x212 + a4x31 = 0 (13.22) If terms involving a2, a3, and a4 are neglected (since a is assumed to be small), Eq. (13.22) will be satisfied if the following equations are satisfied: $ (13.23) x0 + v20 x0 = 0 $ x 1 + v20 x1 = - x30 (13.24) x01t2 = A 0 sin 1v0t + f2 (13.25) The solution of Eq. (13.23) can be expressed as RaoCh13ff.qxd 908 10.06.08 13:57 Page 908 CHAPTER 13 NONLINEAR VIBRATION In view of Eq. (13.25), Eq. (13.24) becomes $ x1 + v20 x1 = - A 30 sin31v0t + f2 = - A 30 C 34 sin 1v0t + f2 - 1 4 sin 31v0t + f2 D (13.26) The particular solution of Eq. (13.26) is (and can be verified by substitution) x11t2 = A 30 3 tA 30 cos1v0t + f2 sin 31v0t + f2 8v0 32 v20 x1t2 = x01t2 + ax11t2 (13.27) Thus the approximate solution of Eq. (13.20) becomes = A 0 sin1v0t + f2 + A 30 a 3a t 3 sin 31v0t + f2 (13.28) A 0 cos1v0t + f2 8 v0 32 v20 The initial conditions on x(t) can be used to evaluate the constants A 0 and f. Notes 1. It can be seen that even a weak nonlinearity (i.e., small value of a) leads to a nonperiodic solution since Eq. (13.28) is not periodic due to the second term on the right-hand side of Eq. (13.28). In general, the solution given by Eq. (13.21) will not be periodic if we retain only a finite number of terms. 2. In Eq. (13.28), the second term, and hence the total solution, can be seen to approach infinity as t tends to infinity. However, the exact solution of Eq. (13.20) is known to be bounded for all values of t. The reason for the unboundedness of the solution, Eq. (13.28), is that only two terms are considered in Eq. (13.21). The second term in Eq. (13.28) is called a secular term. The infinite series in Eq. (13.21) leads to a bounded solution of Eq. (13.20) because the process is a convergent one. To illustrate this point, consider the Taylor’s series expansion of the function sin1vt + at2: sin1v + a2t = sin vt + at cos vt - a3t3 a2t2 sin vt cos vt + Á 2! 3! (13.29) If only two terms are considered on the right-hand side of Eq. (13.29), the solution approaches infinity as t : q . However, the function itself and hence its infinite series expansion can be seen to be a bounded one. 13.4.2 Lindstedt’s Perturbation Method This method assumes that the angular frequency along with the solution varies as a function of the amplitude A 0. This method eliminates the secular terms in each step of the approximation [13.5] by requiring the solution to be periodic in each step. The solution and the angular frequency are assumed as x1t2 = x01t2 + ax11t2 + a2x21t2 + Á (13.30) RaoCh13ff.qxd 10.06.08 13:57 Page 909 13.4 APPROXIMATE ANALYTICAL METHODS v2 = v20 + av11A 02 + a2v21A 02 + Á 909 (13.31) We consider the solution of the pendulum equation, Eq. (13.20), to illustrate the perturbation method. We use only linear terms in a in Eqs. (13.30) and (13.31): x1t2 = x01t2 + ax11t2 v2 = v20 + av11A 02 v20 = v2 - av11A 02 or (13.32) (13.33) $ $ x 0 + ax1 + [v2 - av11A 02][x0 + ax1] + a[x0 + ax1]3 = 0 Substituting Eqs. (13.32) and (13.33) into Eq. (13.20), we get that is, $ $ x 0 + v20 x0 + a1v2x1 + x30 - v1x0 + x12 + a213x1x20 - v1x12 + a313x21 x02 + a41x312 = 0 (13.34) Setting the coefficients of various powers of a to zero and neglecting the terms involving a2, a3, and a4 in Eq. (13.34), we obtain $ x0 + v2x0 = 0 $ x1 + v2x1 = - x30 + v1x0 x01t2 = A 0 sin1vt + f2 (13.35) (13.36) Using the solution of Eq. (13.35), into Eq. (13.36), we obtain (13.37) $ x 1 + v2x1 = - [A 0 sin1vt + f2]3 + v1[A 0 sin1vt + f2] = - 43 A 30 sin1vt + f2 + 14 A 30 sin 31vt + f2 + v1A 0 sin1vt + f2 (13.38) It can be seen that the first and the last terms on the right-hand side of Eq. (13.38) lead to secular terms. They can be eliminated by taking v1 as v1 = 43 A 20, A0 Z 0 (13.39) With this, Eq. (13.38) becomes $ x 1 + v2x1 = 41 A 30 sin 31vt + f2 (13.40) RaoCh13ff.qxd 910 10.06.08 13:57 Page 910 CHAPTER 13 NONLINEAR VIBRATION The solution of Eq. (13.40) is x11t2 = A 1 sin1v t + f12 - A 30 sin 31v t + f2 (13.41) # Let the initial conditions be x1t = 02 = A and x1t = 02 = 0. Using Lindstedt’s method, we force the solution x01t2 given by Eq. (13.37) to satisfy the initial conditions so that # x102 = A = A 0 sin f, x102 = 0 = A 0v cos f 32 v2 or A0 = A f = and p 2 Since the initial conditions are satisfied by x01t2 itself, the solution x11t2 given by Eq. (13.41) must satisfy zero initial conditions.2 Thus x1102 = 0 = A 1 sin f1 - A 30 32 v2 sin 3f A 30 $ x 1102 = 0 = A 1v cos f1 13v2 cos 3f 32 v2 In view of the known relations A 0 = A and f = p/2, the above equations yield A1 = - a A3 b 32 v2 and f1 = p . 2 Thus the total solution of Eq. (13.20) becomes x1t2 = A 0 sin1vt + f2 - aA 30 32 v2 sin 31vt + f2 (13.42) with v2 = v20 + a 34 A 20 (13.43) For the solution obtained by considering three terms in the expansion of Eq. (13.30), see Problem 13.13. It is to be noted that the Lindstedt’s method gives only the periodic solutions of Eq. (13.20); it cannot give any nonperiodic solutions, even if they exist. 13.4.3 Iterative Method In the basic iterative method, first the equation is solved by neglecting certain terms. The resulting solution is then inserted in the terms that were neglected at first to obtain a second, improved, solution. We shall illustrate the iterative method to find the solution of Duffing’s equation, which represents the equation of motion of a damped, harmonically excited, single degree of freedom system with a nonlinear spring. We begin with the solution of the undamped equation. If x01t2 satisfies the initial conditions, each of the solutions x11t2, x21t2, Á appearing in Eq. (13.30) must satisfy zero initial conditions. 2 RaoCh13ff.qxd 10.06.08 13:57 Page 911 13.4 APPROXIMATE ANALYTICAL METHODS 911 Solution of the Undamped Equation. If damping is disregarded, Duffing’s equation becomes $ x + v20 ⫾ ax3 = F cos vt or $ x = - v20x ⫿ ax3 + F cos vt (13.44) x11t2 = A cos vt (13.45) As a first approximation, we assume the solution to be where A is an unknown. By substituting Eq. (13.45) into Eq. (13.44), we obtain the differential equation for the second approximation: $ x 2 = - A v20 cos v t ⫿ A3a cos3v t + F cos v t (13.46) By using the identity cos3 vt = 3 4 cos vt + 1 4 cos 3 vt $ x2 = - 1Av20 ⫾ 43A3a - F2 cos vt ⫿ 14A3a cos 3vt (13.47) Eq. (13.46) can be expressed as (13.48) By integrating this equation and setting the constants of integration to zero (so as to make the solution harmonic with period t = 2p/v), we obtain the second approximation: x21t2 = A3a 1 1Av20 ⫾ 34 A3a - F2cos vt ⫾ cos 3 vt 2 v 36 v2 (13.49) Duffing [13.7] reasoned at this point that if x11t2 and x21t2 are good approximations to the solution x(t), the coefficients of cos vt in the two equations (13.45) and (13.49) should not be very different. Thus by equating these coefficients, we obtain A = 1 3 aAv20 ⫾ A3a - F b 2 4 v or v2 = v20 ⫾ 3 2 F Aa 4 A (13.50) For present purposes, we will stop the procedure with the second approximation. It can be verified that this procedure yields the exact solution for the case of a linear spring (with a = 0) A = v20 F - v2 where A denotes the amplitude of the harmonic response of the linear system. (13.51) RaoCh13ff.qxd 912 10.06.08 13:57 Page 912 CHAPTER 13 NONLINEAR VIBRATION For a nonlinear system (with a Z 0), Eq. (13.50) shows that the frequency v is a function of a, A, and F. Note that the quantity A, in the case of a nonlinear system, is not the amplitude of the harmonic response but only the coefficient of the first term of its solution. However, it is commonly taken as the amplitude of the harmonic response of the system.3 For the free vibration of the nonlinear system, F = 0 and Eq. (13.50) reduces to v2 = v20 ⫾ 43 A2a (13.52) This equation shows that the frequency of the response increases with the amplitude A for the hardening spring and decreases for the softening spring. The solution, Eq. (13.52), can also be seen to be same as the one given by Lindstedt’s method, Eq. (13.43). For both linear and nonlinear systems, when F Z 0 (forced vibration), there are two values of the frequency v for any given amplitude ƒ A ƒ . One of these values of v is smaller and the other larger than the corresponding frequency of free vibration at that amplitude. For the smaller value of v, A 7 0 and the harmonic response of the system is in phase with the external force. For the larger value of v, A 6 0 and the response is 180° out of phase with the external force. Note that only the harmonic solutions of Duffing’s equation—that is, solutions for which the frequency is the same as that of the external force F cos vt—have been considered in the present analysis. It has been observed [13.2] that oscillations whose frequency is a fraction, such as 21, 13, Á , n1 , of that of the applied force, are also possible for Duffing’s equation. Such oscillations, known as subharmonic oscillations, are considered in Section 13.5. Solution of the Damped Equation. If we consider viscous damping, we obtain Duffing’s equation: $ # x + cx + v20 x ⫾ ax3 = F cos vt (13.53) For a damped system, it was observed in earlier chapters that there is a phase difference between the applied force and the response or solution. The usual procedure is to prescribe the applied force first and then determine the phase of the solution. In the present case, however, it is more convenient to fix the phase of the solution and keep the phase of the applied force as a quantity to be determined. We take the differential equation, Eq. (13.53), in the form $ # x + cx + v20 x ⫾ a x3 = F cos1vt + f2 = A 1 cos vt - A 2 sin vt (13.54) in which the amplitude F = 1A 21 + A 2221/2 of the applied force is considered fixed, but the ratio A 1/A 2 = tan -1 f is left to be determined. We assume that c, A 1, and A 2 are all small, of order a. As with Eq. (13.44), we assume the first approximation to the solution to be x1 = A cos vt (13.55) The first approximate solution, Eq. (13.45), can be seen to satisfy the initial conditions x102 = A and # x102 = 0. 3 RaoCh13ff.qxd 10.06.08 13:57 Page 913 13.4 APPROXIMATE ANALYTICAL METHODS 913 where A is assumed fixed and v to be determined. By substituting Eq. (13.55) into Eq. (13.54) and making use of the relation (13.47), we obtain c1v20 - v22A ⫾ 3 3 aA3 aA d cos vt - c vA sin vt ⫾ cos 3 vt 4 4 = A 1 cos vt - A 2 sin vt (13.56) By disregarding the term involving cos 3 vt and equating the coefficients of cos vt and sin vt on both sides of Eq. (13.56), we obtain the following relations: 1v20 - v22A ⫾ 43 aA3 = A 1 c vA = A 2 (13.57) The relation between the amplitude of the applied force and the quantities A and v can be obtained by squaring and adding the equations (13.57): C 1v20 - v22A ⫾ 34 aA3 D 2 + 1c vA22 = A 21 + A 22 = F 2 S21v, A2 + c2v2A2 = F 2 (13.58) Equation (13.58) can be rewritten as (13.59) S1v, A2 = 1v20 - v22A ⫾ 43 aA3 where (13.60) It can be seen that for c = 0, Eq. (13.59) reduces to S1v, A2 = F, which is the same as Eq. (13.50). The response curves given by Eq. (13.59) are shown in Fig. 13.5. Jump Phenomenon. As mentioned earlier, nonlinear systems exhibit phenomena that cannot occur in linear systems. For example, the amplitude of vibration of the system 冟A冟 冟A冟 F0 F2 0 (a)   0 F1 F2 F0 F1 F2 F1 O 冟A冟 F0 F1 F2  O 0 (b)   0 FIGURE 13.5 Response curves of Duffing’s equation.  O 0 (c)   0  RaoCh13ff.qxd 914 10.06.08 13:57 Page 914 CHAPTER 13 NONLINEAR VIBRATION 兩A兩 兩A兩 3 4 6 7 2 1 6 3 4 O 5  (a)   0 (hard spring) FIGURE 13.6 5 12 7  O (b)   0 (soft spring) Jump phenomenon. described by Eq. (13.54) has been found to increase or decrease suddenly as the excitation frequency v is increased or decreased, as shown in Fig. 13.6. For a constant magnitude of F, the amplitude of vibration will increase along the points 1, 2, 3, 4, 5 on the curve when the excitation frequency v is slowly increased. The amplitude of vibration jumps from point 3 to 4 on the curve. Similarly, when the forcing frequency v is slowly decreased, the amplitude of vibration follows the curve along the points 5, 4, 6, 7, 2, 1 and makes a jump from point 6 to 7. This behavior is known as the jump phenomenon. It is evident that two amplitudes of vibration exist for a given forcing frequency, as shown in the shaded regions of the curves of Fig. 13.6. The shaded region can be thought of as unstable in some sense. Thus an understanding of the jump phenomena requires a knowledge of the mathematically involved stability analysis of periodic solutions [13.24, 13.25]. The jump phenomena were also observed experimentally by several investigators [13.26, 13.27]. 13.4.4 Ritz-Galerkin Method In the Ritz-Galerkin method, an approximate solution of the problem is found by satisfying the governing nonlinear equation in the average. To see how the method works, let the nonlinear differential equation be represented as E[x] = 0 (13.61) An approximate solution of Eq. (13.61) is assumed as x' 1t2 = a1f11t2 + a2f21t2 + Á + anfn1t2 (13.62) where f11t2, f21t2, Á , fn1t2 are prescribed functions of time and a1, a2, Á , an are weighting factors to be determined. If Eq. (13.62) is substituted in Eq. (13.61), we get a function E[x' 1t2]. Since x' 1t2 is not, in general, the exact solution of Eq. (13.61), E ' 1t2 = E[x ' 1t2] will not be zero. However, the value of E ' [t] will serve as a measure of the accuracy of the approximation; in fact, E ' [t] : 0 as x ' : x. RaoCh13ff.qxd 10.06.08 13:57 Page 915 915 13.4 APPROXIMATE ANALYTICAL METHODS The weighting factors a i are determined by minimizing the integral L0 t 2 E ' [t] dt (13.63) where t denotes the period of the motion. The minimization of the function of Eq. (13.63) requires 0E 0 ' [t] 2 a dt = 0, E E ' [t] dtb = 2 ' [t] 0ai L0 0ai L0 t t i = 1, 2, Á , n (13.64) Equation (13.64) represents a system of n algebraic equations that can be solved simultaneously to find the values of a1, a2, Á , an. The procedure is illustrated with the following example. Solution of Pendulum Equation Using the Ritz-Galerkin Method EXAMPLE 13.1 Using a one-term approximation, find the solution of the pendulum equation v20 3 $ E[x] = x + v20x x = 0 6 (E.1) by the Ritz-Galerkin method. Solution: By using a one-term approximation for x(t) as x' 1t2 = A 0 sin vt Eqs. (E.1) and (E.2) lead to E[x' 1t2] = - v2A 0 sin vt + v20 cA 0 sin vt = a v20 - v2 - (E.2) 1 3 sin vt d 6 v20 3 1 2 2 v0 A 0 b A 0 sin vt + A sin 3 vt 8 24 0 (E.3) The Ritz-Galerkin method requires the minimization of L0 E 2[x' 1t2] dt t for finding A 0. The application of Eq. (13.64) gives L0 t E ' 0E ' 0A 0 dt = L0 t c a v20 - v2 - * c av20 - v2 - v20 3 1 2 2 v0 A 0 b A 0 sin vt + A sin 3 vt d 8 24 0 1 3 2 2 v A b sin vt + v20A 20 sin 3 vt d dt = 0 8 0 0 8 (E.4) RaoCh13ff.qxd 916 10.06.08 13:57 Page 916 CHAPTER 13 NONLINEAR VIBRATION that is, A 0 a v20 - v2 + 3 1 2 2 v A 0 b a v20 - v2 - v20 A 20 b sin2 vt dt 8 0 8 L0 t t v20A 30 2 3 a v0 - v2 - v20 A 20 b sin vt sin 3 vt dt 24 8 L0 1 1 2 2 2 v A a v - v2 - v20 A 20 b sin vt sin 3 vt dt 8 0 0 0 8 L0 t + + v40A 50 t 2 sin 3 vt dt = 0 192 L0 that is, A 0 c a v20 - v2 - v40A 40 3 1 2 2 v0A 0 b a v20 - v2 - v20A 20 b + d = 0 8 8 192 For a nontrivial solution, A 0 Z 0, and Eq. (E.5) leads to v4 + v2v20 a 12 A 20 - 2b + v40 a 1 - 1 2 A 20 + 5 4 96 A 0 b = 0 (E.5) (E.6) The roots of the quadratic equation in v2, Eq. (E.5), can be found as v2 = v2011 - 0.147938 A 202 v2 = v2011 - 0.352062 A 202 (E.7) (E.8) It can be verified that v2 given by Eq. (E.7) minimizes the quantity of (E.4), while the one given by Eq. (E.8) maximizes it. Thus the solution of Eq. (E.1) is given by Eq. (E.2) with v2 = v2011 - 0.147938 A 202 (E.9) v2 = v2011 - 0.125 A 202 (E.10) This expression can be compared with Lindstedt’s solution and the iteration methods (Eqs. 13.43 and 13.52): The solution can be improved by using a two-term approximation for x(t) as x' 1t2 = A 0 sin vt + A 3 sin 3 vt (E.11) The application of Eq. (13.64) with the solution of Eq. (E.11) leads to two simultaneous algebraic equations that must be numerically solved for A 0 and A 3. ■ RaoCh13ff.qxd 10.06.08 13:57 Page 917 13.5 SUBHARMONIC AND SUPERHARMONIC OSCILLATIONS 917 Other approximate methods, such as the equivalent linearization scheme and the harmonic balance procedure, are also available for solving nonlinear vibration problems [13.10–13.12]. Specific solutions found using these techniques include the free vibration response of single degree of freedom oscillators [13.13, 13.14], two degree of freedom systems [13.15], and elastic beams [13.16, 13.17], and the transient response of forced systems [13.18, 13.19]. Several nonlinear problems of structural dynamics have been discussed by Crandall [13.30]. 13.5 Subharmonic and Superharmonic Oscillations We noted in Chapter 3 that for a linear system, when the applied force has a certain frequency of oscillation, the steady-state response will have the same frequency of oscillation. However, a nonlinear system will exhibit subharmonic and superharmonic oscillations. Subharmonic response involves oscillations whose frequencies 1vn2 are related to the forcing frequency 1v2 as v n (13.65) vn = nv (13.66) vn = where n is an integer 1n = 2, 3, 4, Á 2. Similarly, superharmonic response involves oscillations whose frequencies 1vn2 are related to the forcing frequency 1v2 as where n = 2, 3, 4, Á . 13.5.1 Subharmonic Oscillations In this section, we consider the subharmonic oscillations of order 13 of an undamped pendulum whose equation of motion is given by (undamped Duffing’s equation): $ x + v20x + ax3 = F cos 3 vt (13.67) where a is assumed to be small. We find the response using the perturbation method [13.4, 13.6]. Accordingly, we seek a solution of the form 2 v = v20 x1t2 = x01t2 + ax11t2 + av1 or v20 (13.68) 2 = v - av1 (13.69) where v denotes the fundamental frequency of the solution (equal to the third subharmonic frequency of the forcing frequency). Substituting Eqs. (13.68) and (13.69) into Eq. (13.67) gives $ $ x0 + ax1 + v2x0 + v2ax1 - av1x0 - a2x1v1 + a1x0 + ax123 = F cos 3vt If terms involving a2, a3, and a4 are neglected, Eq. (13.70) reduces to $ $ x0 + v2x0 + ax1 + av2x1 - av1x0 + ax30 = F cos 3 vt (13.70) (13.71) RaoCh13ff.qxd 918 10.06.08 13:57 Page 918 CHAPTER 13 NONLINEAR VIBRATION We first consider the linear equation (by setting a = 0): $ x0 + v2x0 = F cos 3 vt (13.72) (13.73) x01t2 = A 1 cos vt + B1 sin vt + C cos 3 vt # If the initial conditions are assumed as x1t = 02 = A and x1t = 02 = 0, we obtain A 1 = A and B1 = 0 so that Eq. (13.73) reduces to The solution of Eq. (13.72) can be expressed as x01t2 = A cos vt + C cos 3 vt (13.74) where C denotes the amplitude of the forced vibration. The value of C can be determined by substituting Eq. (13.74) into Eq. (13.72) and equating the coefficients of cos 3 vt on both sides of the resulting equation, which yields C = - F 8v2 (13.75) Now we consider the terms involving a in Eq. (13.71) and set them equal to zero: $ a1x1 + v2x1 - v1x0 + x302 = 0 or $ x1 + v2x1 = v1x0 - x30 (13.76) The substitution of Eq. (13.74) into Eq. (13.76) results in $ x 1 + v2x1 = v1A cos vt + v1C cos 3 vt - A3 cos3vt - C 3 cos3 3 vt - 3A2C cos2 vt cos 3 vt - 3AC 2 cos vt cos2 3 vt (13.77) By using the trigonometric relations cos2 u = 3 cos u = cos u cos f = 1 2 3 4 1 2 ∂ cos u + 41 cos 3u 1 cos1u - f2 + 2 cos1u + f2 + 1 2 cos 2u (13.78) Eq. (13.77) can be expressed as 3 3 3 $ x 1 + v2x1 = Aav1 - A2 - C 2 - AC b cos vt 4 2 4 + a v1C - - 3 3 A3 - C 3 - A2Cb cos 3 vt 4 4 2 3 3AC 2 C3 AC1A + C2 cos 5 vt cos 7 vt cos 9 vt 4 4 4 (13.79) RaoCh13ff.qxd 10.06.08 13:57 Page 919 13.5 SUBHARMONIC AND SUPERHARMONIC OSCILLATIONS 919 The condition to avoid a secular term in the solution is that the coefficient of cos vt in Eq. (13.79) must be zero. Since A Z 0 in order to have a subharmonic response, v1 = 43 1A2 + AC + 2C 22 (13.80) Equations (13.80) and (13.75) give v1 = 3 2 AF 2F 2 aA + b 4 8v2 64v4 (13.81) Substituting Eq. (13.81) into Eq. (13.69) and rearranging the terms, we obtain the equation to be satisfied by A and v in order to have subharmonic oscillation as v6 - v20 v4 - 3a 164A2v4 - 8AFv2 + 2F 22 = 0 256 (13.82) Equation (13.82) can be seen to be a cubic equation in v2 and a quadratic in A. The relationship between the amplitude (A) and the subharmonic frequency 1v2, given by Eq. (13.82), is shown graphically in Fig. 13.7. It has been observed that the curve PQ, where the slope is positive, represents stable solutions while the curve QR, where the slope is negative, denotes unstable solutions [13.4, 13.6]. The minimum value of amplitude for the FIGURE 13.7 Subharmonic oscillations. RaoCh13ff.qxd 920 10.06.08 13:57 Page 920 CHAPTER 13 NONLINEAR VIBRATION 13.5.2 Superharmonic Oscillations existence of stable subharmonic oscillations can be found by setting dv2/dA = 0 as A = 1F/16v22.4 Consider the undamped Duffing’s equation $ x + v20x + ax3 = F cos vt (13.83) The solution of this equation is assumed as x1t2 = A cos vt + C cos 3 vt (13.84) where the amplitudes of the harmonic and superharmonic components, A and C, are to be determined. The substitution of Eq. (13.84) into Eq. (13.83) gives, with the use of the trigonometric relations of Eq. (13.78), cos vt C - v2A + v20A + 34 aA3 + 43 aA2C + 32 aAC 2 D + cos 3 vt C - 9v2C + v20C + 14 aA3 + 34 aC 3 + 32 aA2C D + cos 5 vt + cos 9 vt C 34 aA2C + 34 aAC2 D + cos 7 vt C 34 aAC2 D C 14 aC3 D = F cos vt (13.85) Neglecting the terms involving cos 5vt, cos 7vt, and cos 9vt, and equating the coefficients of cos vt and cos 3 vt on both sides of Eq. (13.85), we obtain v20A - v2A + 43 aA3 + 34 aA2C + 32 aAC 2 = F (13.86) v20 C - 9v2 C + 41 aA3 + 34 a C 3 + 23 aA2 C = 0 (13.87) Equations (13.86) and (13.87) represent a set of simultaneous nonlinear equations that can be solved numerically for A and C. As a particular case, if C is assumed to be small compared to A, the terms involving C 2 and C 3 can be neglected and Eq. (13.87) gives C L and Eq. (13.86) gives C L 4 132 aA2 + v20 - 9v22 - 41 aA3 F - v20A + v2A - 34 aA3 3 2 4 aA Equation (13.82) can be rewritten as 1v223 - v201v222 - which, upon differentiation, gives 31v222 dv2 - 2v20v2 dv2 - 3a 2 2 2 3aF 3aF 2 A 1v 2 + A1v22 = 0 4 32 128 3a 3a 2 2 2 3aF 2 3aF 12A dA21v222 A v dv + v dA + A dv2 = 0 4 2 32 32 By setting dv2/dA = 0, we obtain A = 1F/16v22. (13.88) (13.89) RaoCh13ff.qxd 10.06.08 13:57 Page 921 13.6 SYSTEMS WITH TIME-DEPENDENT COEFFICIENTS (MATHIEU EQUATION) 1- 14 aA32143 aA22 = 132 aA2 + v20 - 9 v22 921 Equating C from Eqs. (13.88) and (13.89) leads to * 1F - v20A + v2A - 34 aA32 9 2 3 33 2 2 3 2 - A5115 16 a 2 + A 1 4 av - 4 av02 + A 12 aF2 which can be rewritten as + A110 v2v20 - 9v4 - v402 + 1v20F - 9v2F2 = 0 13.6 (13.90) (13.91) Equation (13.88), in conjunction with Eq. (13.91), gives the relationship between the amplitude of superharmonic oscillations (C) and the corresponding frequency 13v2. Systems with Time-Dependent Coefficients (Mathieu Equation) Consider the simple pendulum shown in Fig. 13.8(a). The pivot point of the pendulum is made to vibrate in the vertical direction as y1t2 = Y cos vt (13.92) where Y is the amplitude and v is the frequency of oscillation. Since the entire pendulum $ accelerates in the vertical direction, the net acceleration is given by g - y1t2 = g - v2Y cos vt. The equation of motion of the pendulum can be derived as $ $ (13.93) ml2 u + m1g - y 2 l sin u = 0 For small deflections near u = 0, sin u L u and Eq. (13.93) reduces to $ g v2Y u + a cos vt bu = 0 l l y(t) Pivot m 0 y  Y cos  t mg   m Pivot 0 y  Y cos  t mg y(t) (a) (b) FIGURE 13.8 Simple pendulum with oscillations of pivot. (13.94) RaoCh13ff.qxd 922 10.06.08 13:57 Page 922 CHAPTER 13 NONLINEAR VIBRATION If the pendulum is inverted as shown in Fig. 13.8(b), the equation of motion becomes $ ml2 u - mgl sin u = 0 or $ g u - sin u = 0 l (13.95) where u is the angle measured from the vertical (unstable equilibrium) point. If the pivot point 0 vibrates as y1t2 = Y cos vt, the equation of motion becomes $ g v2Y u + a- + cos vt b sin u = 0 l l (13.96) For small angular displacements around u = 0, Eq. (13.96) reduces to $ g v2Y u + a- + cos vt b u = 0 l l (13.97) Equations (13.94) and (13.97) are particular forms of an equation called the Mathieu equation for which the coefficient of u in the differential equation varies with time to form a nonautonomous equation. We shall study the periodic solutions and their stability characteristics of the system for small values of Y in this section. Periodic Solutions Using Lindstedt’s Perturbation Method [13.7]. Consider the Mathieu equation in the form d2y dt2 + 1a + P cos t2y = 0 (13.98) where P is assumed to be small. We approximate the solution of Eq. (13.98) as y1t2 = y01t2 + Py11t2 + P2y21t2 + Á a = a0 + Pa1 + P2a2 + Á (13.99) (13.100) where a0, a1, Á are constants. Since the periodic coefficient cos t in Eq. (13.98) varies with a period of 2p, it was found that only two types of solutions exist—one with period 2p and the other with period 4p [13.7, 13.28]. Thus we seek the functions y01t2, y11t2, Á in Eq. (13.99) in such a way that y(t) is a solution of Eq. (13.98) with period 2p or 4p. Substituting Eqs. (13.99) and (13.100) into Eq. (13.98) results in $ $ 1y 0 + a0y02 + P1y 1 + a1y0 + y0 cos t + a0y12 $ + P2(y2 + a2y0 + a1y1 + y1 cos t + a0y2) + Á = 0 (13.101) RaoCh13ff.qxd 10.06.08 13:57 Page 923 13.6 SYSTEMS WITH TIME-DEPENDENT COEFFICIENTS (MATHIEU EQUATION) 923 $ where y i = d2yi /dt2, i = 0, 1, 2, Á Setting the coefficients of various powers of P in Eq. (13.101) equal to zero, we obtain $ y 0 + a0y0 = 0 (13.102) $ y 1 + a0y1 + a1y0 + y0 cos t = 0 (13.103) P0: P1: P2: $ y 2 + a0y2 + a2y0 + a1y1 + y1 cos t = 0 . . . (13.104) where each of the functions yi is required to have a period of 2p or 4p. The solution of Eq. (13.102) can be expressed as cos 2a0 t n t 2 K d , y01t2 = d n sin 2a0 t sin t 2 cos n = 0, 1, 2, Á (13.105) and n2 , 4 a0 = n = 0, 1, 2, Á Now, we consider the following specific values of n. When n = 0: Equation (13.105) gives a0 = 0 and y0 = 1 and Eq. (13.103) yields $ $ (13.106) y 1 + a1 + cos t = 0 or y1 = - a1 - cos t In order to have y1 as a periodic function, a1 must be zero. When Eq. (13.106) is integrated twice, the resulting periodic solution can be expressed as y11t2 = cos t + a (13.107) where a is a constant. With the known values of a0 = 0, a1 = 0, y0 = 1 and y1 = cos t + a, Eq. (13.104) can be rewritten as $ y 2 + a2 + 1cos t + a2cos t = 0 or $ y2 = - 1 2 - a2 - a cos t - 1 2 In order to have y2 as a periodic function, 1 - 21 - a22 must be zero (i.e., a2 = - 12 ). Thus, for n = 0, Eq. (13.100) gives a = - 21 P2 + Á cos 2t (13.108) (13.109) RaoCh13ff.qxd 924 10.06.08 13:57 Page 924 CHAPTER 13 NONLINEAR VIBRATION When n = 1: For this case, Eq. (13.105) gives a0 = With y0 = cos1t/22, Eq. (13.103) gives 1 4 and y0 = cos1t/22 or sin(t/2). 1 1 t 1 3t $ y 1 + y1 = a - a1 - b cos - cos 4 2 2 2 2 (13.110) The homogeneous solution of Eq. (13.110) is given by y11t2 = A 1 cos t t + A 2 sin 2 2 where A 1 and A 2 are constants of integration. Since the term involving cos(t/2) appears in the homogeneous solution as well as the forcing function, the particular solution will contain a term of the form t cos(t/2), which is not periodic. Thus the coefficient of cos(t/2)— namely 1- a1 - 1/22, must be zero in the forcing function to ensure periodicity of y11t2. This gives a 1 = - 1/2, and Eq. (13.110) becomes 1 1 3t $ y 1 + y1 = - cos 4 2 2 (13.111) By substituting the particular solution y11t2 = A 2 cos13t/22 into Eq. (13.111), we obtain A 2 = 41, and hence y11t2 = 41 cos13t/22. Using a 0 = 14, a1 = - 21 and y1 = 41 cos13t/22, Eq. (13.104) can be expressed as 1 t 1 1 3t 1 3t $ y 2 + y2 = - a2 cos + a cos b - a cos b cos t 4 2 2 4 2 4 2 t 1 3t 1 5t 1 = a - a2 - b cos + cos - cos 8 2 8 2 8 2 (13.112) Again, since the homogeneous solution of Eq. (13.112) contains the term cos(t/2), the coefficient of the term cos(t/2) on the right-hand side of Eq. (13.112) must be zero. This leads to a2 = - 18 and hence Eq. (13.100) becomes a = P P2 1 - + Á 4 2 8 (13.113a) Similarly, by starting with the solution y0 = sin1t/22, we obtain the relation (see Problem 13.17) a = P P2 1 + + Á 4 2 8 (13.113b) When n = 2: Equation (13.105) gives a 0 = 1 and y0 = cos t or sin t. With a0 = 1 and y0 = cos t, Eq. (13.103) can be written as 1 1 $ y1 + y1 = - a1 cos t - - cos 2t (13.114) 2 2 Since cos t is a solution of the homogeneous equation, the term involving cos t in Eq. (13.114) gives rise to t cos t in the solution of y1. Thus, to impose periodicity of y1, we RaoCh13ff.qxd 10.06.08 13:57 Page 925 13.6 SYSTEMS WITH TIME-DEPENDENT COEFFICIENTS (MATHIEU EQUATION) 925 set a1 = 0. With this, the particular solution of y11t2 can be assumed as y11t2 = A 3 + B3 cos 2t. When this solution is substituted into Eq. (13.114), we obtain A 3 = - 21 and B3 = 61. Thus Eq. (13.104) becomes $ y 2 + y2 + a2 cos t + y1 cos t = 0 or $ y 2 + y2 = - a2 cos t - 1- 12 + 1 2 1 6 1 12 2 cos 2t2cos t 1 2 For periodicity of y21t2, we set the coefficient of cos t equal to zero in the forcing function 5 of Eq. (13.115). This gives a 2 = 12 , and hence = cos t1 - a1 + a = 1 + 5 2 12 P - + cos 3t + Á (13.115) (13.116a) Similarly, by proceeding with y0 = sin t, we obtain (see Problem 13.17) a = 1 - 1 2 12 P + Á (13.116b) To observe the stability characteristics of the system, Eqs. (13.109), (13.113), and (13.116) are plotted in the 1a, P2 plane as indicated in Fig. 13.9. These equations represent curves FIGURE 13.9 Stability of periodic solutions. RaoCh13ff.qxd 926 13.7 10.06.08 13:57 Page 926 CHAPTER 13 NONLINEAR VIBRATION that are known as the boundary or transition curves that divide the 1a, P2 plane into regions of stability and instability. These boundary curves are such that a point belonging to any one curve represents a periodic solution of Eq. (13.98). The stability of these periodic solutions can be investigated [13.7, 13.25, 13.28]. In Fig. 13.9, the points inside the shaded region denote unstable motion. It can be noticed from this figure that stability is also possible for negative values of a, which correspond to the equilibrium position u = 180°. Thus with the right choice of the parameters, the pendulum can be made to be stable in the upright position by moving its support harmonically. Graphical Methods 13.7.1 Phase Plane Representation Graphical methods can be used to obtain qualitative information about the behavior of the nonlinear system and also to integrate the equations of motion. We shall first consider a basic concept known as the phase plane. For a single degree of freedom system, two parameters are needed to describe the state of motion completely. These parameters are usually taken as the displacement and velocity of the system. When the parameters are used as coordinate axes, the resulting graphical representation of the motion is called the phase plane representation. Thus each point in the phase plane represents a possible state of the system. As time changes, the state of the system changes. A typical or representative point in the phase plane (such as the point representing the state of the system at time t = 0) moves and traces a curve known as the trajectory. The trajectory shows how the solution of the system varies with time. Trajectories of a Simple Harmonic Oscillator EXAMPLE 13.2 Find the trajectories of a simple harmonic oscillator. Solution: The equation of motion of an undamped linear system is given by $ x + v2nx = 0 (E.1) # By setting y = x, Eq. (E.1) can be written as dy = - v2nx dt dx = y dt (E.2) dy v2nx = dx y (E.3) y2 + v2nx2 = c2 (E.4) from which we can obtain Integration of Eq. (E.3) leads to RaoCh13ff.qxd 10.06.08 13:57 Page 927 13.7 GRAPHICAL METHODS 927 yx x FIGURE 13.10 Trajectories of a simple harmonic oscillator. where c is a constant. The value of c is determined by the initial conditions of the system. Equation (E.4) shows that the trajectories of the system in the phase plane (x-y plane) are a family of ellipses, as shown in Fig. 13.10. It can be observed that the point 1x = 0, y = 02 is surrounded by closed trajectories. Such a point is called a center. The direction of motion of the trajectories can be determined from Eq. (E.2). For instance, if x 7 0 and y 7 0, Eq. (E.2) implies that dx/dt 7 0 and dy/dt 6 0; therefore, the motion is clockwise. ■ Phase-plane of an Undamped Pendulum EXAMPLE 13.3 Find the trajectories of an undamped pendulum. Solution: The equation of motion is given by Eq. (13.1): $ u = - v20 sin u (E.1) # # where v20 = g/l. Introducing x = u and y = x = u, Eq. (E.1) can be rewritten as dx = y, dt dy = - v20 sin x dt or dy v20 sin x = dx y or y dy = - v20 sin x dx (E.2) RaoCh13ff.qxd 928 10.06.08 13:57 Page 928 CHAPTER 13 NONLINEAR VIBRATION z  y/0  /0 3 2 2 0 x 3 FIGURE 13.11 Trajectories of an undamped pendulum. # Integrating Eq. (E.2) and using the condition that x = 0 when x = x0 (at the end of the swing), we obtain y2 = 2v201cos x - cos x02 (E.3) Introducing z = y/v0 Eq. (E.3) can be expressed as z2 = 21cos x - cos x02 (E.4) The trajectories given by Eq. (E.4) are shown in Fig. 13.11. ■ Phase-plane of an Undamped Nonlinear System EXAMPLE 13.4 Find the trajectories of a nonlinear spring-mass system governed by the equation $ x + v201x - 2ax32 = 0 (E.1) Solution: The nonlinear pendulum equation can be considered as a special case of Eq. (E.1). To see this, we use the approximation sin u M u - u3/6 in the neighborhood of u = 0 in Eq. (E.1) of Example 13.3 to obtain $ u3 u + v20 a u - b = 0 6 which can be seen to be a special case of Eq. (E.1). Equation (E.1) can be rewritten as dx = y, dt dy = - v201x - 2ax32 dt RaoCh13ff.qxd 10.06.08 13:57 Page 929 13.7 GRAPHICAL METHODS 929 or v201x - 2ax32 dy = dx y or y dy = - v201x - 2ax32 dx (E.2) z2 + x2 - ax4 = A2 (E.3) # Integration of Eq. (E.2), with the condition x = 0 when x = x0 (at the end of the swing in the case of a pendulum), gives where z = y/v0 and A2 = x2011 - a x202 is a constant. The trajectories, in the phase-plane, given by Eq. (E.3), are shown in Fig. 13.12 for several values of a. It can be observed that for a = 0, Eq. (E.3) denotes a circle of radius A and corresponds to a simple harmonic motion. When a 6 0, Eq. (E.3) represents ovals within the circle given by a = 0 and the ovals touch the circle at the points 10, ⫾A2. When a = 11>42 A2, Eq. (E.3) becomes y2 + x2 - x2 x4 x2 b d c y + aA bd = 0 - A2 = cy - a A 2 2A 2A 4A (E.4) Equation (E.4) indicates that the trajectories are given by the parabolas y = ⫾ aA - x2 b 2A FIGURE 13.12 Trajectories of a nonlinear system. (E.5) RaoCh13ff.qxd 930 10.06.08 13:57 Page 930 CHAPTER 13 NONLINEAR VIBRATION These two parabolas intersect at points 1x = ⫾12A, y = 02, which correspond to the points of unstable equilibrium. When 11>42 A2 Ú a Ú 0, the trajectories given by Eq. (E.3) will be closed ovals lying between the circle given by a = 0 and the two parabolas given by a = 11>42A2. Since these trajectories are closed curves, they represent periodic vibrations. When a 7 11>42A2, the trajectories given by Eq. (E.3) lie outside the region between the parabolas and extend to infinity. These trajectories correspond to the conditions that permit the body to escape from the center of force. ■ To see some of the characteristics of trajectories, consider a single degree of freedom nonlinear oscillatory system whose governing equation is of the form $ # (13.117) x + f1x, x2 = 0 By defining dx # = x = y dt (13.118) dy # = y = - f1x, y2 dt (13.119) f1x, y2 dy/dt dy = f1x, y2, say. = = y dx dx/dt (13.120) and we obtain Thus there is a unique slope of the trajectory at every point (x, y) in the phase plane, provided that f1x, y2 is not indeterminate. If y = 0 and f Z 0 (that is, if the point lies on the x axis), the slope of the trajectory is infinite. This means that all trajectories must cross the x axis at right angles. If y = 0 and f = 0, the point is called a singular point, and the slope is indeterminate at such points. A singular point corresponds to a state of equilibrium of # $ the system—the velocity y = x and the force x = - f are zero at a singular point. Further investigation is necessary to establish whether the equilibrium represented by a singular point is stable or unstable. 13.7.2 Phase Velocity ! The velocity v with which a representative point moves along a trajectory is called the phase velocity. The components of phase velocity parallel to the x and y axes are # vx = x, # vy = y (13.121) ! and the magnitude of v is given by dy 2 dx 2 ! ƒ v ƒ = 4v 2x + v 2y = a b + a b C dt dt (13.122) RaoCh13ff.qxd 10.06.08 13:57 Page 931 13.7 GRAPHICAL METHODS 931 We can note that if the system has a periodic motion, its trajectory in the phase plane is a closed curve. This follows from the fact that the representative point, having started its motion along a closed trajectory at an arbitrary point (x, y), will return to the same point after one period. The time required to go around the closed trajectory (the period of oscillation of the system) is finite because the phase velocity is bounded away from zero at all points of the trajectory. 13.7.3 Method of Constructing Trajectories We shall now consider a method known as the method of isoclines for constructing the trajectories of a dynamical system with one degree of freedom. By writing the equations of motion of the system as dx = f11x, y2 dt dy = f21x, y2 dt (13.123) f21x, y2 dy = = f1x, y2, say. dx f11x, y2 (13.124) f1x, y2 = c (13.125) # where f1 and f2 are nonlinear functions of x and y = x, the equation for the integral curves can be obtained as The curve for a fixed value of c is called an isocline. An isocline can be defined as the locus of points at which the trajectories passing through them have the constant slope c. In the method of isoclines we fix the slope (dy)/(dx) by giving it a definite number c1 and solve Eq. (13.125) for the trajectory. The curve f1x, y2 - c1 = 0 thus represents an isocline in the phase plane. We plot several isoclines by giving different values c1, c2, Á to the slope f = 1dy2/1dx2. Let h1, h2, Á denote these isoclines in Fig. 13.13(a). Suppose that we y y (x0, y0) h1 R1 h2 R2 R2 R3 R4 R3 R2 h3 R3 R4 x O h4 R4 x (a) FIGURE 13.13 Method of isoclines. (b) RaoCh13ff.qxd 932 10.06.08 13:57 Page 932 CHAPTER 13 NONLINEAR VIBRATION are interested in constructing the trajectory passing through the point R1 on the isocline h1. We draw two straight line segments from R1: one with a slope c1, meeting h2 at R 2œ , and the other with a slope c2 meeting h2 at R 2fl. The middle point between R 2œ and R 2fl lying on h2 is denoted as R2. Starting at R2, this construction is repeated, and the point R3 is determined on h3. This procedure is continued until the polygonal trajectory with sides R1R2, R2R3, R3R4, Á is taken as an approximation to the actual trajectory passing through the point R1. Obviously, the larger the number of isoclines, the better is the approximation obtained by this graphical method. A typical final trajectory looks like the one shown in Fig. 13.13(b). Trajectories Using the Method of Isoclines EXAMPLE 13.5 Construct the trajectories of a simple harmonic oscillator by the method of isoclines. Solution: The differential equation defining the trajectories of a simple harmonic oscillator is given by Eq. (E.3) of Example 13.2. Hence the family of isoclines is given by v2nx y c = - or y = - v2n x c (E.1) This equation represents a family of straight lines passing through the origin, with c representing the slope of the trajectories on each isocline. The isoclines given by Eq. (E.1) are shown in Fig. 13.14. Once the isoclines are known, the trajectory can be plotted as indicated above. ■ 13.7.4 Obtaining Time Solution from Phase Plane Trajectories # The trajectory plotted in the phase plane is a plot of x as a function of x, and time (t) does not appear explicitly in the plot. For the qualitative analysis of the system, the trajectories are enough, but in some cases we may need the variation of x with time t. In such cases, it is possible to obtain the time solution x(t) from the phase plane diagram, although the original y Trajectory c0 c 1 4 c O x c c1 1 c 4 FIGURE 13.14 Isoclines of a simple harmonic oscillator. 1 RaoCh13ff.qxd 10.06.08 13:57 Page 933 13.8 STABILITY OF EQUILIBRIUM STATES 933 FIGURE 13.15 Obtaining time solution from a phase plane plot. # differential equation cannot be solved for x and x as functions of time. The method of obtaining a time solution is essentially a step-by-step procedure; several different schemes may be used for this purpose. In this section, we shall present a method based on the rela# tion x = 1¢x2/1¢t2. For small increments of displacement and time ( ¢x and ¢t), the average velocity can # be taken as xav = 1¢x2/1¢t2, so that ¢x ¢t = # xav (13.126) In the phase plane trajectory shown in Fig. 13.15, the incremental time needed for the rep# resentative point to traverse the incremental displacement ¢xAB is shown as ¢tAB. If xAB # denotes the average velocity during ¢tAB, we have ¢tAB = ¢xAB /xAB. Similarly, # ¢tBC = ¢xBC/xBC, etc. Once ¢tAB, ¢tBC, Á are known, the time solution x(t) can be plotted easily, as shown in Fig. 13.15(b). It is evident that for good accuracy, the incremental displacements ¢xAB, ¢xBC, Á must be chosen small enough that the correspon# ding incremental changes in x and t are reasonably small. Note that ¢x need not be constant; it can be changed depending on the nature of the trajectories. 13.8 13.8.1 Stability Analysis Stability of Equilibrium States Consider a single degree of freedom nonlinear vibratory system described by two firstorder differential equations dx = f11x, y2 dt dy = f21x, y2 dt (13.127) RaoCh13ff.qxd 934 10.06.08 13:57 Page 934 CHAPTER 13 NONLINEAR VIBRATION # where f1 and f2 are nonlinear functions of x and y = x = dx/dt. For this system, the slope of the trajectories in the phase plane is given by # f21x, y2 y dy = # = (13.128) dx x f11x, y2 Let 1x0, y02 be a singular point or an equilibrium point so that (dy)/(dx) has the form 0/0: f11x0, y02 = f21x0, y02 = 0 (13.129) A study of Eqs. (13.127) in the neighborhood of the singular point provides us with answers as to the stability of equilibrium. We first note that there is no loss of generality if we assume that the singular point is located at the origin (0, 0). This is because the slope (dy)/(dx) of the trajectories does not vary with a translation of the coordinate axes x and y to x¿ and y¿: x¿ = x - x0 y¿ = y - y0 dy dy¿ = dx dx¿ (13.130) Thus we assume that 1x = 0, y = 02 is a singular point, so that f110, 02 = f210, 02 = 0 If f1 and f2 are expanded in terms of Taylor’s series about the singular point (0, 0), we obtain # x = f11x, y2 = a11x + a12y + Higher order terms # y = f21x, y2 = a21x + a22y + Higher order terms (13.131) where a11 = 0f1 , ` 0x 10, 02 a12 = 0f1 , ` 0y 10, 02 a21 = 0f2 , ` 0x 10, 02 a22 = 0f2 ` 0y 10, 02 In the neighborhood of (0, 0), x and y are small; f1 and f2 can be approximated by linear terms only, so that Eqs. (13.131) can be written as # x a a12 x e # f = c 11 de f (13.132) y a21 a22 y The solutions of Eq. (13.132) are expected to be geometrically similar to those of Eq. (13.127). We assume the solution of Eq. (13.132) in the form x X e f = e f elt y Y (13.133) RaoCh13ff.qxd 10.06.08 13:57 Page 935 13.8 STABILITY OF EQUILIBRIUM STATES 935 where X, Y, and l are constants. Substitution of Eq. (13.133) into Eq. (13.132) leads to the eigenvalue problem c a11 - l a21 ` a12 X 0 de f = e f a22 - l Y 0 (13.134) a12 ` =0 a22 - l The eigenvalues l1 and l2 can be found by solving the characteristic equation a11 - l a21 l1, l2 = 121p ⫾ 4p2 - 4q2 as where p = a11 + a22 and q = a11a22 - a12a21. If e X1 f Y1 and e (13.135) X2 f Y2 denote the eigenvectors corresponding to l1 and l2, respectively, the general solution of Eqs. (13.127) can be expressed as (assuming l1 Z 0, l2 Z 0, and l1 Z l2): x X X e f = C1 e 1 fel1t + C2 e 2 fel2t y Y1 Y2 (13.136) If 1p2 - 4q2 6 0, the motion is oscillatory. where C1 and C2 are arbitrary constants. We can note the following: If 1p2 - 4q2 7 0, the motion is aperiodic. If p 7 0, the system is unstable. If p 6 0, the system is stable. If we use the transformation x X e f = c 1 y Y1 X2 a a d e f = [T]e f Y2 b b where [T] is the modal matrix and a and b are the generalized coordinates, Eqs. (13.132) will be uncoupled: # # a l 0 a a = l 1a # (13.137) e #f = c 1 de f or 0 l2 b b = l2 b b The solution of Eqs. (13.137) can be expressed as a1t2 = el1t b1t2 = el2t (13.138) RaoCh13ff.qxd 936 10.06.08 13:57 Page 936 CHAPTER 13 NONLINEAR VIBRATION 13.8.2 Classification of Singular Points Depending on the values of l1 and l2 in Eq. (13.135), the singular or equilibrium points can be classified as follows [13.20, 13.23]. Case (i) — L1 and L2 Are Real and Distinct (p 2>4q). In this case, Eq. (13.138) gives a1t2 = a0el1t b1t2 = b 0el2t and (13.139) where a0 and b 0 are the initial values of a and b, respectively. The type of motion depends on whether l1 and l2 are of the same sign or of opposite signs. If l1 and l2 have the same sign 1q 7 02, the equilibrium point is called a node. The phase plane diagram for the case l2 6 l1 6 0 (when l1 and l2 are real and negative or p 6 0) is shown in Fig. 13.16(a). In this case, Eq. (13.139) shows that all the trajectories tend to the origin as t : q and hence the origin is called a stable node. On the other hand, if l2 7 l1 7 0 1p 7 02, the arrow heads change in direction and the origin is called an unstable node (see Fig. 13.16b). If l1 and l2 are real but of opposite signs (q 6 0 irrespective of the sign of p), one solution tends to the origin while the other tends to infinity. In this case the origin is called a saddle point and it corresponds to unstable equilibrium (see Fig. 13.16d). y yx y x (a) Stable node y (b) Unstable node y x (d) Saddle point x x y x (e) Stable focus (c) Stable node y x (f) Unstable focus FIGURE 13.16 Types of equilibrium points. x (g) Center RaoCh13ff.qxd 10.06.08 13:57 Page 937 13.9 LIMIT CYCLES Case (ii) — L1 and L2 Are Real and Equal (p 2 ⴝ 4q). a1t2 = a0el1t and 937 In this case Eq. (13.138) gives b1t2 = b 0el1t (13.140) The trajectories will be straight lines passing through the origin and the equilibrium point (origin) will be a stable node if l1 6 0 (see Fig. 13.16c) and an unstable node if l1 7 0. Case (iii) — L1 and L2 Are Complex Conjugates (p 2<4q). Let l1 = u1 + iu2 and l2 = u1 - iu2, where u1 and u2 are real. Then Eq. (13.137) gives # # a = 1u1 + iu22a and b = 1u1 - iu22b (13.141) a1t2 = 1a0eu1t2eiu2t, b1t2 = 1b 0eu1t2e -iu2t This shows that a and b must also be complex conjugates. We can rewrite Eq. (13.138) as (13.142) These equations represent logarithmic spirals. In this case the equilibrium point (origin) is called a focus or a spiral point. Since the factor eiu2t in a1t2 represents a vector of unit magnitude rotating with angular velocity u2 in the complex plane, the magnitude of the complex vector a1t2, and hence the stability of motion, is determined by eu1t. If u1 6 0, the motion will be asymptotically stable and the focal point will be stable (p 6 0 and q 7 0). If u1 6 0, the focal point will be unstable (p 7 0 and q 7 0). The sign of u2 merely gives the direction of rotation of the complex vector, counterclockwise for u2 7 0 and clockwise for u2 6 0. If u1 = 0 1p = 02, the magnitude of the complex radius vector a1t2 will be constant and the trajectories reduce to circles with the center as the equilibrium point (origin). The motion will be periodic and hence will be stable. The equilibrium point in this case is called a center or vertex point and the motion is simply stable and not asymptotically stable. The stable focus, unstable focus, and center are shown in Figs. 13.16(e) to (g). 13.9 Limit Cycles In certain vibration problems involving nonlinear damping, the trajectories, starting either very close to the origin or far away from the origin, tend to a single closed curve, which corresponds to a periodic solution of the system. This means that every solution of the system tends to a periodic solution as t : q . The closed curve to which all the solutions approach is called a limit cycle. For illustration, we consider the following equation, known as the van der Pol equation: $ # x - a11 - x22x + x = 0, a 7 0 (13.143) This equation exhibits, mathematically, the essential features of some vibratory systems like certain electrical feedback circuits controlled by valves where there is a source of power that increases with the amplitude of vibration. Van der Pol invented Eq. (13.143) by introducing a type of damping that is negative for small amplitudes but becomes positive for large # amplitudes. In this equation, he assumed the damping term to be a multiple of - 11 - x 22x in order to make the magnitude of the damping term independent of the sign of x. RaoCh13ff.qxd 938 10.06.08 13:57 Page 938 CHAPTER 13 NONLINEAR VIBRATION The qualitative nature of the solution depends upon the value of the parameter a. Although the analytical solution of the equation is not known, it can be represented using the method of isoclines, in the phase plane. Equation (13.143) can be rewritten as dx # y = x = dt dy # y = = a11 - x22y - x dt (13.144) (13.145) a11 - x22y - x dy/dt dy = c = = y dx dx/dt Thus the isocline corresponding to a specified value of the slope dy/dx = c is given by or x d = 0 - a11 - x22 + c y + c (13.146) By drawing the curves represented by Eq. (13.146) for a set of values of c, as shown in Fig. 13.17, the trajectories can be sketched in with fair accuracy, as shown in the same figure. The isoclines will be curves since the equation, Eq. (13.146), is nonlinear. An infinity of isoclines pass through the origin, which is a singularity. An interesting property of the solution can be observed from Fig. 13.17. All the trajectories, irrespective of the initial conditions, approach asymptotically a particular closed y c 2 0 0.5 1 0.5 0 2 3 2 1 c 3 2 1 0 1 2 3 1 2 1 3 c 2 0 0.5 1 0.5 0 2 FIGURE 13.17 Trajectories and limit cycle for van der Pol’s equation. x RaoCh13ff.qxd 10.06.08 13:57 Page 939 13.10 CHAOS 939 curve, known as the limit cycle, which represents a steady-state periodic (but not harmonic) oscillation. This is a phenomenon that can be observed only with certain nonlinear vibration problems and not in any linear problem. If the initial point is inside the limit cycle, the ensuing solution curve spirals outward. On the other hand, if the initial point falls outside the limit cycle, the solution curve spirals inward. As stated above, the limit cycle in both the cases is attained finally. An important characteristic of the limit cycle is that the maximum value of x is always close to 2 irrespective of the value of a. This result can be seen by solving Eq. (13.143) using the perturbation method (see Problem 13.34). 13.10 Chaos Chaos represents the behavior of a system that is inherently unpredictable. In other words, chaos refers to the dynamic behavior of a system whose response, although described by a deterministic equation, becomes unpredictable because the nonlinearities in the equation enormously amplify the errors in the initial conditions of the system. Attractor. Consider a pendulum whose amplitude of oscillation decreases gradually due to friction, which means that the system loses part of its energy in each cycle and eventually comes to a rest position. This is indicated in the phase plane shown in Fig. 13.18(a). The rest position is called an attractor. Thus the pendulum has just one attractor. If the pendulum is given a push at the end of each swing to replenish the energy lost due to friction, the resulting motion can be indicated as a closed loop in the phase plane (see Fig. 13.18b). In general, for a dynamic system, an attractor is a point (or object) toward which all nearby solutions move as time progresses. Poincaré Section. A pendulum is said to have two degrees of freedom—namely, x and # x. In general, a phase space of a system can be defined with as many axes as there are degrees of freedom. Thus, for a system with three degrees of freedom, the phase space might appear (as a spiral converging in the z-direction) as shown in Fig. 13.19(a). Since the points are displaced from one another and never coincide in Fig. 13.19(a), the system does not have a periodic motion. The intersection of the phase space with the yz plane · Velocity (x) · Velocity (x) Amplitude (x) (a) FIGURE 13.18 Attractor. Amplitude (x) (b) RaoCh13ff.qxd 940 10.06.08 13:57 Page 940 CHAPTER 13 NONLINEAR VIBRATION x z 6 7 4 2 5 3 1 y (a) z 7 6 5 4 3 2 1 y (b) FIGURE 13.19 Phase space for a three degree of freedom system. appears as shown in Fig. 13.19(b). This diagram is known as the Poincaré section or Poincaré map, and it represents points that occur at equal intervals of time, nT1n = 1, 2, Á 2, where T denotes the fundamental period of the forcing function. Note that, if the system is periodic, all the dots would be at the same location in the Poincaré section. 13.10.1 Functions with Stable Orbits Consider the sequence of numbers generated by the following equation: xn + 1 = 2xn; n = 1, 2, Á (13.147) For any two initial values of x1, which differ by a small amount, the values of xn + 1 converge to 1. For example, with x1 = 10.0 and x1 = 10.2, the sequences of numbers given by Eq. (13.147) are shown below: 10.0 : 3.1623 : 1.7783 : 1.3335 : 1.1548 : 1.0746 : 1.0366 : 1.0182 : 1.0090 : 1.0045 : 1.0023 : 1.0011 : 1.0006 : 1.0003 : 1.0001 : 1.0001 : 1.0000 RaoCh13ff.qxd 10.06.08 13:57 Page 941 941 13.10 CHAOS and 10.2 : 3.1937 : 1.7871 : 1.3368 : 1.1562 : 1.0753 : 1.0370 : 1.0183 : 1.0091 : 1.0045 : 1.0023 : 1.0011 : 1.0006 : 1.0003 : 1.0001 : 1.0001 : 1.0000 Note that the influence of a change in the initial value of x1 (by 0.2) is lost very soon and the pattern converges to 1. Any starting value between 0 and 1 also iterates to 1. Thus the functional relation, Eq. (13.147), is said to have a stable orbit at x = 1. 13.10.2 Functions with Unstable Orbits xn + 1 = a xn 11 - xn2; n = 1, 2, Á Consider the sequence of numbers generated by the following equation: (13.148) where a is a constant. Equation (13.148) has been used as a simple model for population growth with no predators such as fish and fowl. In such cases, the constant a denotes the rate of growth of the population, xn represents the population in generation n, and the factor 11 - xn2 acts as a stabilizing factor. It can be seen that Eq. (13.148) has the following limitations [13.31]: 1. x1 has to lie between 0 and 1. If x1 exceeds 1, the iterative process diverges to - q , implying that the population becomes extinct. 2. xn + 1 attains a maximum value of a4 at xn = 21. This implies that a 6 4. 3. If a 6 1, xn + 1 converges to zero. 4. Thus, for a nontrivial dynamic behavior (to avoid population from becoming extinct), a has to satisfy the relation 1 … a … 4. The system will have an equilibrium condition if the birth rate replenishes the loss due to death or migration. The population can be seen to stabilize or reach a definite limiting value (predictable) for some values of a —such as a = 3.0. For some other values of a— such as a = 4.0 with x1 = 0.5—the species can be seen to disappear after only two generations, as indicated below: Equation 113.1482 with a = 4.0 and x1 = 0.5: 0.5 : 1.0 : 0.0 : 0.0 : 0.0 : Á However, for a = 4.0 with x1 = 0.4, the population count can be seen to be completely random, as indicated below: Equation 113.1482 with a = 4.0 and x1 = 0.4: 0.4 : 0.96 : 0.154 : 0.520 : 0.998 : 0.006 : 0.025 : 0.099 : 0.358 : 0.919 : 0.298 : 0.837 : 0.547 : 0.991 : 0.035 : 0.135 : 0.466 : 0.996 : 0.018 : Á This indicates that the system is a chaotic one; even small changes in the deterministic equation, Eq. (13.148), can lead to unpredictable results. Physically, this implies that the system has become chaotic with population varying wildly from year to year. In fact, as will be shown in the following paragraph, the system, Eq. (13.148), has unstable orbits. RaoCh13ff.qxd 942 10.06.08 13:57 Page 942 CHAPTER 13 NONLINEAR VIBRATION Bifurcations. Equation (13.148) exhibits a phenomenon known as bifurcation. To see this, we start with a = 2 and a few different values of x1. With this, xn + 1 can be seen to converge to 0.5. When we start with a = 2.5 and a few different values of x1, the process converges to 0.6. If we use a = 3.0 and x1 = 0.1, the process converges to a single value, but while getting there, oscillates between two separate values (namely, 0.68 Á and 0.65 Á ). If we use a = 3.25 and x1 = 0.5, then the value of xn + 1 oscillates between the two values x112 = 0.4952 and x122 = 0.8124. At that point the system is said to have a periodicity of 2. The solution, in this case, moves into a two-pronged fork-type of state with two equilibrium points. If we use a = 3.5 and x1 = 0.5, the system will have a period 4 —that is, the equilibrium state will oscillate between the four values x132 = 0.3828, x142 = 0.5008, x152 = 0.8269, and x162 = 0.8749. This implies that the stable behavior of each of the previous two solutions has been broken into further bifurcation paths. In fact, the system continues to bifurcate with the range of a needed for each successive birth of bifurcations becoming smaller as a increases, as shown in Fig. 13.20. Figure 13.20 is known as a bifurcation plot or Feigenbaum diagram [13.31, 13.35]. It can be observed that the system has reached a chaotic state through a series of bifurcations, with the number of values assumed by Eq. (13.148) doubling at each stage. Strange Attractors. For several years, it was believed that the attractors toward which physical systems tend are equilibrium or rest points (as in the case of the rest position of a pendulum), limit cycles, or repeating configurations. However, in recent years it has been found that the attractors associated with chaotic systems are more complex than the rest points and limit cycles. The geometric points in state space to which chaotic trajectories are attracted are called strange attractors. xn1 x(6) 1.0 0.9 x(2) x(5) 0.8 0.7 0.6 Chaotic region x(4) 0.5 x(1) 0.4 x(3) 0.3 0.2 0.0 2.5 a8 a4 0.1 a2 3.0 3.5 4.0 ai  Range of a for period i (i  2, 4, 8, . . . ) FIGURE 13.20 Bifurcation plot. a RaoCh13ff.qxd 10.06.08 13:57 Page 943 13.10 CHAOS 13.10.3 Chaotic Behavior of Duffing’s Equation Without the Forcing Term 943 Consider a single degree of freedom system with a nonlinear spring and a harmonic forcing function. The equation of motion (Duffing’s equation) can be expressed as $ # x + m x - a x + b x3 = F0 cos v t (13.149) First, we consider the free, undamped vibration of the system with a = b = 0.5: $ (13.150) x - 0.5 x + 0.5 x3 = 0 The static equilibrium positions of this system can be found by setting the spring force equal to zero, as x = 0, +1, - 1. It can be easily verified that the equilibrium solution x = 0 is unstable (saddle point) with respect to infinitesimal disturbances, while the equilibrium solutions - 1 and + 1 are stable (centers) with respect to infinitesimal disturbances. The stability of the system about the three equilibrium positions can be seen more clearly from the graph of its potential energy. To find the potential energy of the system, we mul# tiply Eq. (13.150) by x and integrate the resulting equation to obtain 1 1 1 # 2 1x2 + x4 - x2 = C 2 8 4 (13.151) # # x1t = 02 = x0, x1t = 02 = x0 (13.153) where C is a constant. The first term on the left-hand side of Eq. (13.151) represents the kinetic energy and the second and third terms denote the potential energy (P) of the system. Equation (13.151) indicates that the sum of the kinetic and potential energies is a constant (conservative system). A plot of the potential energy, P = 18 x4 - 41 x2, is shown in Fig. 13.21. Next, we consider the free, damped vibration of the system. The governing equation is $ # (13.152) x + m x - 0.5 x + 0.5 x3 = 0 Let the boundary conditions be given by It can be expected from Fig. 13.21 that the static equilibria, x = + 1 and x = - 1, are unstable with respect to finite disturbances. Physically, when a finite disturbance is given P Unstable (saddle) point Stable point (center) 1 0 1 Stable point (center) x FIGURE 13.21 Plot of potential energy. RaoCh13ff.qxd 944 10.06.08 13:57 Page 944 CHAPTER 13 NONLINEAR VIBRATION about one static equilibrium point, the system could be driven to the other static equilibrium point. In fact, the steady-state solution can be shown to be extremely sensitive to the initial conditions, exhibiting a form of chaos. It can be verified easily [13.32] that for # x0 = 1 and 0 6 x0 6 0.521799, the steady-state solution is x1t : q 2 = + 1. The # phase plane trajectory for x0 = 0.52 is shown in Fig. 13.22(a). Note that x 7 0 for all # values of t. For x0 = 1 and 0.521799 6 x0 6 0.5572, the steady-state solution is # x1t : q 2 = - 1. Figure 13.22(b) shows the phase plane trajectory for x0 = 0.54, indi# cating a single crossing of the x = 0 axis. For 0.5572 6 x0 6 0.5952, the steady-state # solution is x1t : q 2 = + 1. The phase plane trajectory for x0 = 0.56 is shown in Fig. 13.22(c), which indicates two crossings of the x = 0 axis. # In fact, by giving a series of values to x0, we can construct several phase plane trajectories, from which a composite plot, known as the shell plot, can be constructed as shown in Fig. 13.23 [13.33, 13.34]. Here also, the steady-state solution can be seen to be # x = + 1 or - 1 depending on the initial conditions, x0 and x0. It can be seen that the various regions are identified by the numbers 0, 1, 2, 3, Á First, consider the region labeled “0” with x0 Ú 0. If the initial conditions fall in this region, the solution spirals into x = + 1 as t : q and the solution crosses the axis x = 0 zero times (similar to Fig. 13.22a). Next, consider the region labeled “1.” If the initial conditions fall in this region, the solution moves clockwise, crosses the axis x = 0 once, and settles into x = - 1 as t : q (similar to Fig. 13.22b). Next, consider the region labeled “2.” If the system starts with the initial conditions falling in this region, the phase plane moves clockwise, crosses the axis x = 0, continues to move clockwise, crosses the axis x = 0 again, moves into the region “0” for x 7 0, and spirals into x = + 1 as t : q (similar to Fig. 13.22c). This pattern continues with other regions as well, with the labeled region number indicating the number of crossings of the axis x = 0 by the phase plane trajectory. Figure 13.23 indicates that if there is sufficient uncertainty in the initial conditions x0 # and x0, the final state of the system, x = + 1 or -1, is unpredictable or uncertain. If damping is reduced further, the width of each region in Fig. 13.23 (except the regions labeled “0”) becomes even smaller and vanishes as m : 0. Thus the final steady state of the sys# tem is unpredictable as m : 0 for any finite uncertainty in x0, x0, or both. This means that the system exhibits chaos. 13.10.4 Chaotic Behavior of Duffing’s Equation with the Forcing Term Consider Duffing’s equation with m = 0.2, and a = b = 1 including the forcing term. Within the forcing term, small variations in the frequency 1v2 or the amplitude 1F02 can lead to chaos, as indicated below. When V Is Changed. For a fixed value of F0, the phase plane response of Eq. (13.149) can be periodic or chaotic, depending on the value of v. For example, Figs. 13.24 and 13.25 indicate the two situations that are described by the equations $ # x + 0.2 x - x + x3 = 0.3 cos 1.4 t 1periodic, Fig. 13.242 $ # x + 0.2 x - x + x3 = 0.3 cos 1.29 t 1chaotic, Fig. 13.252 (13.154) (13.155) where F0 = 0.3 has been assumed. Figure 13.24 has been plotted using an approximate analysis, known as the harmonic balance method. On the other hand, Fig. 13.25 represents RaoCh13ff.qxd 10.06.08 13:57 Page 945 13.10 CHAOS . Velocity (x) 0.60 0.48 0.36 0.24 0.12 0.00 0.12 0.24 0.36 0.48 0.60 . x0  1.0, x0  0.52 0 0.24 0.48 0.72 0.96 (a) 1.20 1.44 Displacement (x) . Velocity (x) 0.60 0.48 0.36 0.24 0.12 0.00 0.12 0.24 0.36 0.48 0.60 1.44 . x0  1.0, x0  0.54 0.96 0.48 0 0.48 0.96 1.44 Displacement (x) (b) . Velocity (x) 0.60 0.48 0.36 0.24 0.12 0.00 0.12 0.24 0.36 0.48 0.60 . x0  1.0, x0  0.56 1.5 1.2 0.9 0.6 0.3 0.0 0.3 0.6 0.9 1.2 1.5 (c) Displacement (x) FIGURE 13.22 Phase-plane trajectories for different initial velocities. (From [13.32]; reprinted with permission of Society of Industrial and Applied Mathematics and E. H. Dowell and C. Pierre.) 945 946 10.06.08 13:57 Page 946 CHAPTER 13 NONLINEAR VIBRATION 0.96 0.72   0.0168 0.48 Velocity RaoCh13ff.qxd 0.24 0 0.24 0.48 0.72 0.96 1.8 1.2 0.6 0 0.6 1.2 1.8 Displacement FIGURE 13.23 Shell plot. (From [13.33]; reprinted with permission of American Society of Mechanical Engineers.) 0.80 0.60 0.40 0.20 x· 0.00 0.20 0.40 0.60 0.80 0.40 0.60 0.80 1.00 1.20 1.40 0.50 0.70 0.90 1.10 1.30 x FIGURE 13.24 Phase plane of Eq. (13.154). (From [13.34]; reprinted with permission of Academic Press.) a Poincaré map, which indicates points that occur at equal intervals of time 2p T0, 2T0, 3T0, Á where T0 is the fundamental period of excitation, T0 = 2p v = 1.29 . When F0 Is Changed. Chaos can also be observed when the amplitude of the force changes. To illustrate, we consider the following equation [13.33]: $ # x + 0.168 x - 0.5 x + 0.5 x3 = F0 sin vt K F0 sin t (13.156) # For definiteness, we assume the initial conditions as x0 = 1 and x0 = 0. When F0 is sufficiently small, the response of the system will be a simple harmonic motion (that is, the phase plane will be an ellipse) about its static equilibrium position, x = + 1. When F0 is increased, additional harmonics beyond the fundamental are detected and the phase plane will be distorted from a simple ellipse, as shown in Fig. 13.26(a) for F0 = 0.177. Note that 13:57 Page 947 13.10 CHAOS 947 0.80 0.60 0.40 0.20 x· 0.00 0.20 0.40 0.60 0.80 2.00 1.20 1.60 0.80 0.40 0.40 1.20 0.00 0.80 1.60 x FIGURE 13.25 Poincaré map of Eq. (13.155). (From [13.34]; reprinted with permission of Academic Press.) 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 Velocity Velocity 1.6 1.2 0.8 0.4 0 0.4 0.8 1.2 1.6 Displacement   1, F0  0.177,   0.168 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1.6 1.2 0.8 0.4 0 0.4 0.8 1.2 1.6 Displacement   1, F0  0.178,   0.168 (b) 2 Period (a) 1 Period 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 Velocity 10.06.08 Velocity RaoCh13ff.qxd 1.6 1.2 0.8 0.4 0 0.4 0.8 1.2 1.6 Displacement   1, F0  0.197,   0.168 (c) 4 Period 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1.6 1.2 0.8 0.4 0 0.4 0.8 1.2 1.6 Displacement Trajectory of chaos F0  0.21,   0.168 (d) FIGURE 13.26 Distortion of phase plane. (From [13.33]; reprinted with permission of American Society of Mechanical Engineers.) RaoCh13ff.qxd 948 10.06.08 13:57 Page 948 CHAPTER 13 NONLINEAR VIBRATION the boundary of the region labeled “0” in Fig. 13.23 is also shown in Fig. 13.26(a) for a comparison. The response is called a 1 period motion for 0 … F0 … 0.177, implying that the response oscillates through one period while the force oscillates through one period. For F0 = 0.178, the phase plane is shown in Fig. 13.26(b), which indicates that the response is a 2 period motion. Thus, for the response to oscillate through one period, the force must oscillate through two periods. This change from a 1 to a 2 period response as F0 changes from 0.177 to 0.178 is called a bifurcation. When F0 = 0.197, the response will be a 4 period motion (see Fig. 13.26c). As F0 is increased, 8, 16, Á period motions occur and finally, for F0 Ú 0.205, chaos can be observed with no apparent periodicity, as indicated in Fig. 13.26(d). 13.11 Numerical Methods Most of the numerical methods described in the earlier chapters can be used for finding the response of nonlinear systems. The Runge-Kutta method described in Section 11.4 is directly applicable for nonlinear systems and is illustrated in Section 13.12. The central difference, Houbolt, Wilson, and Newmark methods considered in Chapter 11 can also be used for solving nonlinear multidegree of freedom vibration problems with slight modification. Let a multidegree of freedom system be governed by the equation ! ! ! $ # [m] : x 1t2 + [c] : x 1t2 + P1x(t2) = F1t2 (13.157) ! P are assumed to be nonlinear where the internal set of forces opposing ! the displacements ! ! ! functions of x. For the linear case, P = [k]x. In order to find the displacement vector x that satisfies the nonlinear equilibrium in Eq. (13.157), it is necessary to perform an equilibrium iteration sequence in each time step. In implicit methods (Houbolt, Wilson, and Newmark methods), the equilibrium conditions are considered at the same time for which solution is sought. If the solution is known for time ti and we wish to find the solution for time ti + 1, then the following equilibrium equations are considered ! ! $ # (13.158) x i + 1 + [c] : x i + 1 + Pi + 1 = Fi + 1 [m] : ! ! ! where Fi + 1 = F1t = ti + 12 and Pi + 1 is computed as ! ! ! ! ! ! (13.159) Pi + 1 = Pi + [ki]¢x = Pi + [ki]1xi + 1 - xi2 where [ki] is the linearized or tangent stiffness matrix computed at time ti. Substitution of Eq. (13.159) in Eq. (13.158) gives ! ! $ # ! ! N! [m] : x i + 1 + [c] : x i + 1 + [ki] xi + 1 = Fi + 1 = Fi + 1 - Pi + [ki] xi (13.160) ! Since the right-hand side of Eq. (13.160) is completely known, it can be solved for xi + 1 ! using any of the implicit methods directly. The xi + 1 found is only an approximate vector due to the linearization process used in Eq. (13.159). To improve the accuracy of the solution and to avoid the development of numerical instabilities, an iterative process has to be used within the current time step [13.21]. RaoCh13ff.qxd 10.06.08 13:57 Page 949 13.12 EXAMPLES USING MATLAB 13.12 949 Examples Using MATLAB Solution of the Pendulum Equation EXAMPLE 13.6 Using MATLAB, find the solution of the following forms g = 0.09. v0 = Al $ u + v20 u = (a) $ (b) u + v20 u - 16 v20 u3 = $ u + v20 sin u = (c) Use the following initial conditions: 1i2 u102 = 0.1, p 1ii2 u102 = , 4 p 1iii2 u102 = , 2 of the pendulum equation with 0 (E.1) 0 (E.2) 0 (E.3) # u102 = 0 # u102 = 0 # u102 = 0 (E.4) (E.5) (E.6) # Solution: Using x1 = u and x2 = u, each of the Eqs. (E.1) to (E.3) can be rewritten as a system of two first-order differential equations as follows: (a) (b) (c) # x1 = x2 # x2 = - v20 x1 1Linear equation2 (E.7) 1 # x2 = - v20 x1 + v20 x31 1Nonlinear equation2 6 (E.8) # x1 = x2 # x1 = x2 # x2 = - v20 sin x1 1Nonlinear equation2 (E.9) Equations (E.7) to (E.9) are solved using the MATLAB program ode23 for each of the initial conditions given by Eqs. (E.4) to (E.6). The solutions u1t2 given by Eqs. (E.7), (E.8), and (E.9) for a specific initial conditions, are plotted in the same graph. % Ex13_6.m % This program will use the function dfunc1_a.m, dfunc1_b.m and dfun1_c.m % they should be in the same folder tspan = [0: 1: 250]; RaoCh13ff.qxd 950 10.06.08 13:57 Page 950 CHAPTER 13 NONLINEAR VIBRATION x0 = [0.1; 0.0]; x0_1 = [0.7854; 0.0]; x0_2 = [1.5708; 0.0]; [t, xa] = ode23 ('dfunc1_a', tspan, x0); [t, xb] = ode23 ('dfunc1_b', tspan, x0); [t, xc] = ode23 ('dfunc1_c', tspan, x0); [t, xa1] = ode23 ('dfunc1_a', tspan, x0_1); [t, xb1] = ode23 ('dfunc1_b', tspan, x0_1); [t, xc1] = ode23 ('dfunc1_c', tspan, x0_1); [t, xa2] = ode23 ('dfunc1_a', tspan, x0_2); [t, xb2] = ode23 ('dfunc1_b', tspan, x0_2); [t, xc2] = ode23 ('dfunc1_c', tspan, x0_2); plot (t, xa(:, 1)); ylabel ('Theta(t)'); xlabel ('t'); ylabel ('i.c. = [0.1; 0.0]'); title. . . ('Function a: solid line, Function b: dashed line, Function c: dotted line'); hold on; plot (t, xb(:, 1), '--'); hold on; plot (t, xc(:, 1), ':'); pause; hold off; plot (t, xa1(:, 1)); ylabel ('Theta(t)'); xlabel ('t'); ylabel ('i.c. = [0.7854; 0.0]'); title. . . ('Function a: solid line, Function b: dashed line, Function c: dotted line'); hold on; plot (t, xb1(:, 1), '--'); hold on; plot (t, xc1(:, 1), ':'); pause; hold off; plot(t,xa2(:,1)); hold on; ylabel('Theta(t)'); xlabel('t'); ylabel('i.c. = [1.5708; 0.0]') title. . . ('Function a: solid line, Function b: dashed line, Function c: dotted line'); plot(t, xb2(:,1),'--'); hold on; plot(t,xc2(:,1),':'); % dfunc1_a.m function f = dfunc1_a(t,x); f = zeros(2,1); f(1) = x(2); f(2) = 0.0081 * x(1); % dfunc1_b.m function f = dfunc1_b(t,x); f = zeros(2,1); f(1) = x(2); f(2) = 0.0081 * ((x(1))^3) / 6.0  0.0081 * x(1); % dfunc1_c.m function f = dfunc1_c(t,x); f = zeros(2,1); f(1) = x(2); f(2) = 0.0081 * sin(x(1)); RaoCh13ff.qxd 10.06.08 13:57 Page 951 13.12 EXAMPLES USING MATLAB 951 RaoCh13ff.qxd 952 10.06.08 13:57 Page 952 CHAPTER 13 NONLINEAR VIBRATION Function a: solid line, Function b: dashed line, Function c: dotted line 2 1.5 i.c.  [1.5708; 0.0] 1 0.5 0 0.5 1 1.5 2 0 50 150 100 200 250 t ■ Solution of Nonlinearly Damped System EXAMPLE 13.7 Using MATLAB, find the solution of a single degree of freedom system with velocity squared damping. Governing equation: $ # # m x + c1x22 sign1x2 + k x = F0 sin vt (E.1) $ # m x + ceq x + k x = F0 sin v t (E.2) # Data: m = 10, c = 0.01, k = 4000, F0 = 200, v = 10 and 20, x102 = 0.5, x102 = 1.0 Also find the solution of the system using the equivalent viscous damping constant 1ceq2 where ceq is given by Eq. (E.4) of Example 3.7 as ceq = 8 cv X 3p (E.3) RaoCh13ff.qxd 10.06.08 13:57 Page 953 953 13.12 EXAMPLES USING MATLAB # Solution: By introducing x1 = x and x2 = x, Eqs. (E.1) and (E.2) are written as systems of two first-order differential equations as # 1a2 x1 = x2 F0 c 2 k # x2 = sin v t x sign 1x22 x m m 2 m 1 # 1b2 x1 = x2 ceq F0 k # x2 = sin v t x x m m 2 m 1 1Nonlinear equation2 1Linear equation2 (E.4) (E.5) F and X, in Eq. (E.3), is taken as the steady-state or static deflection of the system as X = k0. The MATLAB solutions given by Eqs. (E.4) and (E.5) are plotted in the same graph for a specific value of v. % Ex13_7.m % This program will use the function dfunc3_a.m, dfunc3_b.m % dfunc3_a1.m, dfunc3_b1.m, they should be in the same folder tspan = [0: 0.005: 10]; x0 = [0.5; 1.0]; [t,xa] = ode23 ('dfunc3_a', tspan, x0); [t,xb] = ode23 ('dfunc3_b', tspan, x0); [t,xa1] = ode23 ('dfunc3_a1', tspan, x0); [t,xb1] = ode23 ('dfunc3_b1', tspan, x0); subplot (211); plot (t,xa (:,1)); title ('Theta(t): function a (Solid line), function b (Dashed line)'); ylabel ('w = 10 '); hold on; plot (t,xb(:,1), '--'); subplot (212); plot (t,xa1(:,1)); ylabel ('w = 20 '); hold on; plot (t,xb1 (:,1), '--'); xlabel ('t'); % dfunc3_a.m function f = dfunc3_a (t,x); f0 = 200; m = 10; a = 0.01; k = 4000; w = 10; f = zeros (2,1); f(1) = x(2); f(2) = f0* sin (w*t) /m  a* x(2)^2 * sign(x(2)) /m  k*x(1) /m; % dfunc3_a1.m function f = dfunc3_a1 (t,x); f0 = 200; m = 10; a = 0.01; k = 4000; w = 20; f = zeros (2,1); RaoCh13ff.qxd 954 10.06.08 13:57 Page 954 CHAPTER 13 NONLINEAR VIBRATION f(1) = x(2); f(2) = f0* sin (w*t) /m  a* x(2)^2 * sign (x(2)) /m  k*x(1) /m; % dfunc3_b.m function f = dfunc3_b (t,x); f0 = 200; m = 10; a = 0.01; k = 4000; ceq = sqrt (8*a*f0 / (3*pi)); w = 10; f = zeros (2,1); f(1) = x(2); f(2) = f0* sin (w*t) /m  ceq * x(2) /m  k*x(1) /m; % dfunc3_b1.m function f = dfunc3_b1 (t,x); f0 = 200; m = 10; a = 0.01; k = 4000; ceq = sqrt (8*a*f0 / (3*pi)); w = 20; f = zeros (2,1); f(1) = x(2); f(2) = f0* sin (w*t) /m  ceq * x(2) /m  k*x(1) /m; ■ RaoCh13ff.qxd 10.06.08 13:57 Page 955 13.12 EXAMPLES USING MATLAB 955 Solution of Nonlinear System Under Pulse Loading 13.8 Using MATLAB, find the solution of a nonlinear single degree of freedom system governed by the equation $ m x + k1 x + k2 x3 = f1t2 (E.1) where F(t) is a rectangular pulse load of magnitude F0 applied over 0 … t … t0. Assume the fol# lowing data: m = 10, k1 = 4000, F0 = 1000, t0 = 1.0, x102 = 0.05, x102 = 5. Solve Eq. (E.1) for two cases: one with k2 = 0 and the other with k2 = 500. # Solution: Using x1 = x and x2 = x, Eq. (E.1) is rewritten, as a set of two first-order differential equations as # x1 = x2 F1t2 k1 k2 3 # x1 x x2 = m m m 1 (E.2) The responses, x(t), found with k2 = 0 (linear system) and k2 = 500 (nonlinear system) are plotted in the same graph. 0.3 0.2 0.1 0 0.1 x(t) EXAMPLE 0.2 0.3 Solid line: k 2  500 0.4 Dashed line: k 2  0 0.5 0.6 0.7 0 0.5 1 1.5 2 2.5 t 3 3.5 4 4.5 5 RaoCh13ff.qxd 956 10.06.08 13:57 Page 956 CHAPTER 13 NONLINEAR VIBRATION % Ex13_8.m % This program will use the function dfunc13_8_1.m and dfunc13_8_2.m % they should be in the same folder tspan = [0: 0.01: 5]; x0 = [0.05; 5]; [t,x] = ode23('dfunc13_8_1', tspan, x0); plot(t,x(:,1)); xlabel('t'); ylabel('x(t)'); hold on; [t,x] = ode23('dfunc13_8_2', tspan, x0); plot(t,x(:,1),'--'); gtext('Solid line: k_2 = 500'); gtext('Dashed line: k_2 = 0') % dfunc13_8_1.m function f = dfunc13_8_1(t,x) f = zeros(2,1); m = 10; k1 = 4000; k2 = 500; F0 = 1000; F = F0 * (stepfun(t, 0) ⴚ stepfun(t, 1)); f(1) = x(2); f(2) = ⴚF/m ⴚ k1 * x(1)/m ⴚ k2 * (x(1))^3/m; % dfunc13_8_2.m function f = dfunc13_8_2(t,x) f = zeros(2,1); m = 10; k1 = 4000; k2 = 0; F0 = 1000; F = F0 * (stepfun(t, 0) ⴚ stepfun(t, 1)); f(1) = x(2); f(2) = ⴚF/m ⴚ k1 * x(1)/m ⴚ k2 * (x(1))^3/m; ■ Solution of Nonlinear Differential Equation EXAMPLE 13.9 Using the fourth-order Runge-Kutta method, develop a general MATLAB program called Program18.m to find the solution of a single degree of freedom equation of the form $ # m x + c x + k x + k* x3 = 0 (E.1) Use the program to solve Eq. (E.1) for the following data: m = 0.01, c = 0.1, k = 2.0, k* = 0.5, # x102 = 7.5, x102 = 0. Solution: Equation (E.1) is rewritten as # x1 = f11x1, x22 = x2 k k* 3 c # x2 = f21x1, x22 = x x x m 2 m 1 m 1 (E.2) 13:57 Page 957 13.12 EXAMPLES USING MATLAB 957 10 5 X(1) 10.06.08 0 5 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 Time 0.6 0.7 0.8 0.9 1 200 100 X(2) RaoCh13ff.qxd 0 100 200 300 Program18.m is developed to accept the values of m, c, k, and k* as YM, YC, YK, and YKS, respectively. The time step 1¢t2 and the number of time steps (NSTEP) are specified as 0.0025 and 400, respectively. A subprogram, called fun; is to be given to define f11x1, x22 and f21x1, x22. The pro# gram gives the values of ti, x1ti2 and x1ti2, i = 1, 2, Á , NSTEP as output. The program also plots # x1t2 = x11t2 and x1t2 = x21t2. Solution of nonlinear vibration problem by fourth order Runge-Kutta method Data: ym = yc = yk = yks = 1.000000e002 1.000000e001 2.00000000e+000 5.00000000e001 Results i 1 6 11 16 21 26 . . . time(i) 2.500000e003 1.500000e002 2.750000e002 4.000000e002 5.250000e002 6.500000e002 x(i,1) 7.430295e+000 5.405670e+000 2.226943e+000 8.046611e001 3.430513e+000 5.296623e+000 x(i,2) 5.528573e+001 2.363166e+002 2.554475e+002 2.280328e+002 1.877713e+002 1.002752e+002 RaoCh13ff.qxd 958 10.06.08 13:57 Page 958 CHAPTER 13 NONLINEAR VIBRATION 371 376 381 386 391 396 9.275000e001 9.400000e001 9.525000e001 9.650000e001 9.775000e001 9.900000e001 1.219287e001 1.209954e001 1.166138e001 1.093188e001 9.966846e002 8.822462e002 7.673075e002 2.194914e001 4.744062e001 6.853283e001 8.512093e001 9.724752e001 ■ 13.13 C++ Program An interactive C++ program called Program18.cpp is given for finding the solution of a single degree of freedom system with a nonlinear spring using the fourth-order RungeKutta method. The input and output of the program are similar to those of the MATLAB program, Program18.m, given in Example 13.9. Solution of Nonlinear Spring-Mass-Damper System EXAMPLE 13.10 Use Program18.cpp to solve the problem described in Example 13.9. Solution: The input data are to be entered interactively. The input and output of the program are given below. SOLUTION OF NONLINEAR VIBRATION PROBLEM BY FOURTH ORDER RUNGE-KUTTA METHOD DATA: YM = 0.010000 YC = 0.100000 YK = 2.000000 YKS = 0.500000 RESULTS: I TIME(I) X(I,1) X(I,2) 1 2 3 4 5 . . . 395 396 397 398 399 400 0.00250000 0.00500000 0.00750000 0.01000000 0.01250000 7.43029541 7.22711455 6.90397627 6.47900312 5.97267903 55.28573134 106.35020086 150.94361827 187.66152131 216.01444338 0.98750000 0.99000000 0.99250000 0.99500000 0.99750000 1.00000000 0.09063023 0.08822462 0.08576929 0.08326851 0.08072652 0.07814748 0.95172604 0.97247522 0.99150470 1.00883370 1.02448309 1.03847532 ■ RaoCh13ff.qxd 10.06.08 13:57 Page 959 REFERENCES 13.14 959 Fortran Program A Fortran program, called PROGRAM18.F, is given for the numerical solution of the nonlinear free vibration problem of a single degree of freedom system using the fourth-order Runge-Kutta method. The input and output of the program are similar to those of the MATLAB program, Program18.m, given in Example 13.9. Solution of Nonlinear Vibration Problem EXAMPLE 13.11 Find the solution of the problem described in Example 13.9 using PROGRAM18.F. Solution: The output of the program is given below. SOLUTION OF NONLINEAR VIBRATION PROBLEM BY FOURTH ORDER RUNGE-KUTTA METHOD DATA: YM = YC = YK = YKS = 0.100000Eⴚ01 0.100000Eⴙ00 0.200000Eⴙ01 0.500000Eⴙ00 RESULTS: I TIME(I) X(I,1) X(I,2) 1 2 3 4 5 . . . 395 396 397 398 399 400 0.002500 0.005000 0.007500 0.010000 0.012500 0.743030E+01 0.722711E+01 0.690398E+01 0.647900E+01 0.597268E+01 0.552857E+02 0.106350E+03 0.150944E+03 0.187662E+03 0.216014E+03 0.987511 0.990011 0.992511 0.995011 0.997511 1.000011 0.906302E01 0.882246E01 0.857693E01 0.832685E01 0.807265E01 0.781475E01 0.951725E+00 0.972474E+00 0.991504E+00 0.100883E+01 0.102448E+01 0.103847E+01 ■ REFERENCES 13.1 C. Hayashi, Nonlinear Oscillations in Physical Systems, McGraw-Hill, New York, 1964. 13.2 A. A. Andronow and C. E. Chaikin, Theory of Oscillations (English language edition), Princeton University Press, Princeton, N.J., 1949. 13.3 N. V. Butenin, Elements of the Theory of Nonlinear Oscillations, Blaisdell Publishing, New York, 1965. 13.4 A. Blaquiere, Nonlinear System Analysis, Academic Press, New York, 1966. 13.5 Y. H. Ku, Analysis and Control of Nonlinear Systems, Ronald Press, New York, 1958. 13.6 J. N. J. Cunningham, Introduction to Nonlinear Analysis, McGraw-Hill, New York, 1958. 13.7 J. J. Stoker, Nonlinear Vibrations in Mechanical and Electrical Systems, Interscience Publishers, New York, 1950. RaoCh13ff.qxd 960 10.06.08 13:57 Page 960 CHAPTER 13 NONLINEAR VIBRATION 13.8 J. P. Den Hartog, Mechanical Vibrations (4th ed.), McGraw-Hill, New York, 1956. 13.9 N. Minorsky, Nonlinear Oscillations, D. Van Nostrand, Princeton, N.J., 1962. 13.10 R. E. Mickens, “Perturbation solution of a highly nonlinear oscillation equation,” Journal of Sound and Vibration, Vol. 68, 1980, pp. 153–155. 13.11 B. V. Dasarathy and P. Srinivasan, “Study of a class of nonlinear systems reducible to equivalent linear systems,” Journal of Sound and Vibration, Vol. 7, 1968, pp. 27–30. 13.12 G. L. Anderson, “A modified perturbation method for treating nonlinear oscillation problems,” Journal of Sound and Vibration, Vol. 38, 1975, pp. 451–464. 13.13 B. L. Ly, “A note on the free vibration of a nonlinear oscillator,” Journal of Sound and Vibration, Vol. 68, 1980, pp. 307–309. 13.14 V. A. Bapat and P. Srinivasan, “Free vibrations of nonlinear cubic spring mass systems in the presence of Coulomb damping,” Journal of Sound and Vibration, Vol. 11, 1970, pp. 121–137. 13.15 H. R. Srirangarajan, P. Srinivasan, and B. V. Dasarathy, “Analysis of two degrees of freedom systems through weighted mean square linearization approach,” Journal of Sound and Vibration, Vol. 36, 1974, pp. 119–131. 13.16 S. R. Woodall, “On the large amplitude oscillations of a thin elastic beam,” International Journal of Nonlinear Mechanics, Vol. 1, 1966, pp. 217–238. 13.17 D. A. Evenson, “Nonlinear vibrations of beams with various boundary conditions,” AIAA Journal, Vol. 6, 1968, pp. 370–372. 13.18 M. E. Beshai and M. A. Dokainish, “The transient response of a forced nonlinear system,” Journal of Sound and Vibration, Vol. 41, 1975, pp. 53–62. 13.19 V. A. Bapat and P. Srinivasan, “Response of undamped nonlinear spring mass systems subjected to constant force excitation,” Journal of Sound and Vibration, Vol. 9, 1969, Part I: pp. 53–58 and Part II: pp. 438–446. 13.20 W. E. Boyce and R. C. DiPrima, Elementary Differential Equations and Boundary Value Problems (4th ed.), Wiley, New York, 1986. 13.21 D. R. J. Owen, “Implicit finite element methods for the dynamic transient analysis of solids with particular reference to nonlinear situations,” in Advanced Structural Dynamics, J. Donéa (ed.), Applied Science Publishers, London, 1980, pp. 123–152. 13.22 B. van der Pol, “Relaxation oscillations,” Philosophical Magazine, Vol. 2, pp. 978–992, 1926. 13.23 L. A. Pipes and L. R. Harvill, Applied Mathematics for Engineers and Physicists (3rd ed.), McGraw-Hill, New York, 1970. 13.24 N. N. Bogoliubov and Y. A. Mitropolsky, Asymptotic Methods in the Theory of Nonlinear Oscillations, Hindustan Publishing, Delhi, 1961. 13.25 A. H. Nayfeh and D. T. Mook, Nonlinear Oscillations, Wiley, New York, 1979. 13.26 G. Duffing, “Erzwungene Schwingungen bei veranderlicher Eigenfrequenz und ihre technische Bedeutung,” Ph.D. thesis (Sammlung Vieweg, Braunschweig, 1918). 13.27 C. A. Ludeke, “An experimental investigation of forced vibrations in a mechanical system having a nonlinear restoring force,” Journal of Applied Physics, Vol. 17, pp. 603–609, 1946. 13.28 D. W. Jordan and P. Smith, Nonlinear Ordinary Differential Equations (2nd ed.), Clarendon Press, Oxford, 1987. 13.29 R. E. Mickens, An Introduction to Nonlinear Oscillations, Cambridge University Press, Cambridge, 1981. RaoCh13ff.qxd 10.06.08 13:57 Page 961 REVIEW QUESTIONS 961 13.30 S. H. Crandall, “Nonlinearities in structural dynamics,” The Shock and Vibration Digest, Vol. 6, No. 8, August 1974, pp. 2–14. 13.31 R. M. May, “Simple mathematical models with very complicated dynamics,” Nature, Vol. 261, June 1976, pp. 459–467. 13.32 E. H. Dowell and C. Pierre, “Chaotic oscillations in mechanical systems,” in Chaos in Nonlinear Dynamical Systems, J. Chandra (ed.), SIAM, Philadelphia, 1984, pp. 176–191. 13.33 E. H. Dowell and C. Pezeshki, “On the understanding of chaos in Duffing’s equation including a comparison with experiment,” ASME Journal of Applied Mechanics, Vol. 53, March 1986, pp. 5–9. 13.34 B. H. Tongue, “Existence of chaos on a one-degree-of-freedom system,” Journal of Sound and Vibration, Vol. 110, No. 1, October, 1986, pp. 69–78. 13.35 M. Cartmell, Introduction to Linear, Parametric, and Nonlinear Vibrations, Chapman and Hall, London, 1990. REVIEW QUESTIONS 13.1 Give brief answers to the following: 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. How do you recognize a nonlinear vibration problem? What are the various sources of nonlinearity in a vibration problem? What is the source of nonlinearity in Duffing’s equation? How is the frequency of the solution of Duffing’s equation affected by the nature of the spring? What are subharmonic oscillations? Explain the jump phenomenon. What principle is used in the Ritz-Galerkin method? Define these terms: phase plane, trajectory, singular point, phase velocity. What is the method of isoclines? What is the difference between a hard spring and a soft spring? Explain the difference between subharmonic and superharmonic oscillations. What is a secular term? Give an example of a system that leads to an equation of motion with time-dependent coefficients. Explain the significance of the following: stable node, unstable node, saddle point, focus, and center. What is a limit cycle? Give two examples of physical phenomena that can be represented by van der Pol’s equation. 13.2 Indicate whether each of the following statements is true or false: 1. Nonlinearity can be introduced into the governing differential equation through mass, spring and/or dampers. 2. Nonlinear analysis of a system can reveal several unexpected phenomena. 3. The Mathieu equation is an autonomous equation. 4. A singular point corresponds to a state of equilibrium of the system. 5. Jump phenomenon is exhibited by both linear and nonlinear systems. RaoCh13ff.qxd 962 10.06.08 13:57 Page 962 CHAPTER 13 NONLINEAR VIBRATION 6. The Ritz-Galerkin method finds the approximate solution by satisfying the nonlinear equation in the average. 7. Dry friction can introduce nonlinearity in the system. 8. Poincaré’s solution of nonlinear equations is in the form of a series. 9. The secular term appears in the solution of free Duffing’s equation. 10. According to Lindstedt’s perturbation method, the angular frequency is assumed to be a function of the amplitude. 11. An isocline is the locus of points at which the trajectories passing through them have a constant slope. 12. Time does not appear explicitly in a trajectory plotted in the phase plane. 13. The time variation of the solution can be found from the phase plane trajectories. 14. A limit cycle denotes a steady-state periodic oscillation. 15. Approximate solutions of nonlinear vibration problems can be found using numerical methods such as Houbolt, Wilson, and Neumark method. 13.3 Fill in each of the following blanks with the appropriate words: 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. When finite amplitudes of motion are involved, _____ analysis becomes necessary. _____ principle is not applicable in nonlinear analysis. _____ equation involves time-dependent coefficients. The governing equation of a simple pendulum whose pivot is subjected to vertical vibration is called _____ equation. The representation of the motion of a system in the displacement-velocity plane is known as _____ plane representation. The curve traced by a typical point in the phase plane is called a _____. The velocity with which a representative point moves along a trajectory is called the _____ velocity. The phenomenon of realizing two amplitudes for the same frequency is known as _____ phenomenon. The forced vibration solution of Duffing’s equation has _____ values of the frequency v for any given amplitude ƒ A ƒ . The Ritz-Galerkin method involves the solution of _____ equations. Mechanical chatter is a _____ vibration. If time does not appear explicitly in the governing equation, the corresponding system is said to be _____. The method of _____ can be used to construct the trajectories of a one degree of freedom dynamical system. Van der Pol’s equation exhibits _____ cycles. 13.4 Select the most appropriate answer out of the choices given: 1. Each term in the equation of motion of a linear system involves displacement, velocity, and acceleration of the (a) first degree (b) second degree (c) zero degree 2. A nonlinear stress-strain curve can lead to nonlinearity of the (a) mass (b) spring (c) damper 3. If the rate of change of force with respect to displacement, df/dx, is an increasing function of x, the spring is called a (a) soft spring (b) hard spring (c) linear spring RaoCh13ff.qxd 10.06.08 13:57 Page 963 REVIEW QUESTIONS 963 4. If the rate of change of force with respect to displacement, df/dx, is a decreasing function of k, the spring is called a (a) soft spring (b) hard spring (c) linear spring 5. The point surrounded by closed trajectories is called a (a) center (b) mid-point (c) focal point 6. For a system with periodic motion, the trajectory in the phase plane is a (a) closed curve (b) open curve (c) point 7. In subharmonic oscillations, the natural frequency 1vn2 and the forcing frequency 1v2 are related as (a) vn = v (b) vn = nv; n = 2, 3, 4, Á v (c) vn = ; n = 2, 3, 4, Á n 8. In superharmonic oscillations, the natural frequency 1vn2 and the forcing frequency 1v2 are related as (a) vn = v (b) vn = nv; n = 2, 3, 4, Á v (c) vn = ; n = 2, 3, 4, Á n 9. If time appears explicitly in the governing equation, the corresponding system is called (a) an autonomous system (b) a nonautonomous system (c) a linear system. 10. Duffing’s equation is given by $ (a) x + v20 x + a x3 = 0 $ (b) x + v20 x = 0 $ (c) x + a x3 = 0 11. Lindstedt’s perturbation method gives (a) periodic and nonperiodic solutions (b) periodic solutions only (c) nonperiodic solutions only 13.5 Match the items in the two columns below for the nature of equilibrium points in the context of the stability analysis of equilibrium states with l1 and l2 as eigenvalues: (1) l1 and l2 with same sign (l1, l2: real and distinct) (2) l1 and l2 6 0 (l1, l2: real and distinct) (3) l1 and l2 7 0 (l1, l2: real and distinct) (4) l1 and l2: real with opposite signs (5) l1 and l2: complex conjugates (a) Unstable node (b) Saddle point (c) Node (d) Focus or spiral point (e) Stable node RaoCh13ff.qxd 964 10.06.08 13:57 Page 964 CHAPTER 13 NONLINEAR VIBRATION 13.6 Match the items in the two columns below: # x $ (1) x + f # + v2n x = 0 (a) Nonlinearity in mass ƒxƒ x3 $ b = 0 (2) x + v20 a x (b) Nonlinearity in damping 6 $ (3) a x x + kx = 0 (c) Linear equation $ # 3 (4) x + cx + k x = a x (d) Nonlinearity in spring PROBLEMS The problem assignments are organized as follows: Problems Section Covered 13.1 13.2–13.8 13.9–13.11 13.12, 13.13 13.14–13.16 13.1 13.2 13.3 13.4 13.5 13.17 13.18–13.22 13.23–13.33 13.34 13.35, 13.36 13.37–13.42 13.43–13.45 13.46–13.49 13.50–13.51 13.6 13.7 13.8 13.9 13.10 13.12 13.13 13.14 — Topic Covered Introduction Examples of nonlinear vibration Exact methods of solution Approximate analytical methods Subharmonic and superharmonic oscillations Mathieu equation Graphical methods Stability of equilibrium states Limit cycles Chaos MATLAB programs C++ program Fortran program Design projects 13.1 The equation of motion of a simple pendulum, subjected to a constant torque, Mt = ml2f, is given by $ u + v20 sin u = f (E.1) $ v20 3 u + v20 u = f + u 6 (E.2) If sin u is replaced by its two-term approximation, u - 1u3/62, the equation of motion becomes RaoCh13ff.qxd 10.06.08 13:57 Page 965 PROBLEMS 965 Let the solution of the linearized equation $ u + v20 u = f (E.3) $ v20 3 u + v20 u = u 6 (E.4) be denoted as u11t2, and the solution of the equation be denoted as u21t2. Discuss the validity of the total solution. u1t2, given by u1t2 = u11t2 + u21t2, for Eq. (E.2). 13.2 Two springs, having different stiffnesses k1 and k2 with k2 7 k1, are placed on either side of a mass m, as shown in Fig. 13.27. When the mass is in its equilibrium position, no spring is in contact with the mass. However, when the mass is displaced from its equilibrium position, # only one spring will be compressed. If the mass is given an initial velocity x0 at t = 0, determine (a) the maximum deflection and (b) the period of vibration of the mass. x(t) k2  k1 k1 m FIGURE 13.27 13.3 Find the equation of motion of the mass shown in Fig. 13.28. Draw the spring force versus x diagram. P sin t x(t) x0 k k x0 m k k FIGURE 13.28 13.4 Two masses m1 and m2 are attached to a stretched wire, as shown in Fig. 13.29. If the initial tension in the wire is P, derive the equations of motion for large transverse displacements of the masses. m1 l1 FIGURE 13.29 m2 l2 l3 RaoCh13ff.qxd 966 10.06.08 13:57 Page 966 CHAPTER 13 NONLINEAR VIBRATION 13.5 A mass m, connected to an elastic rubber band of unstretched length l and stiffness k, is permitted to swing as a pendulum bob, as shown in Fig. 13.30. Derive the nonlinear equations of motion of the system using x and u as coordinates. Linearize the equations of motion and determine the natural frequencies of vibration of the system. Rubber band, stiffness k  lx m FIGURE 13.30 13.6 A uniform bar of length l and mass m is hinged at one end 1x = 02, supported by a spring at x = 2l3 , and acted by a force at x = l, as shown in Fig. 13.31. Derive the nonlinear equation of motion of the system. F(t) 2l 3 l 3  k FIGURE 13.31 13.7 Derive the nonlinear equation of motion of the spring-mass system shown in Fig. 13.32 F(t) x(t) l m k1 k2 FIGURE 13.32 h RaoCh13ff.qxd 10.06.08 13:57 Page 967 PROBLEMS 967 13.8 Derive the nonlinear equations of motion of the system shown in Fig. 13.33. Also, find the linearized equations of motion for small displacements, x(t) and u1t2. F(t) x(t) k m Mt(t) Rigid bar, length l (t) (massless) m FIGURE 13.33 13.9 Find the natural time period of oscillation of the pendulum shown in Fig. 13.1(a) when it oscillates between the limits u = - p/2 and u = p/2, using Eqs. (13.1) and (13.12). 13.10 A simple pendulum of length 30 in. is released from the initial position of 80° from the vertical. How long does it take to reach the position u = 0°? 13.11 Find the exact solution of the nonlinear pendulum equation $ u3 u + v20 a u - b = 0 6 # with u = 0 when u = u0, where u0 denotes the maximum angular displacement. 13.12 Find the solution of Example 13.1 using the following two-term approximation for x(t): x' 1t2 = A 0 sin v t + A 3 sin 3 v t 13.13 Using a three-term expansion in Lindstedt’s perturbation method [Eq. (13.30)], find the solution of the pendulum equation, Eq. (13.20). 13.14 The equation of motion for the forced vibration of a single degree of freedom nonlinear system can be expressed as # $ x + cx + k1x + k2 x3 = a1 cos 3v t - a2 sin 3v t Derive the conditions for the existence of subharmonics of order 3 for this system. 13.15 The equation of motion of a nonlinear system is given by $ # x + cx + k1x + k2x2 = a cos 2 v t Investigate the subharmonic solution of order 2 for this system. 13.16 Prove that, for the system considered in Section 13.5.1, the minimum value of v2 for which the amplitude of subharmonic oscillations A will have a real value is given by vmin = v0 + 21 F 2 2048 v50 RaoCh13ff.qxd 968 10.06.08 13:57 Page 968 CHAPTER 13 NONLINEAR VIBRATION Also, show that the minimum value of the amplitude, for stable subharmonic oscillations, is given by A min = F 16 v2 13.17 Derive Eqs. (13.113b) and (13.116b) for the Mathieu equation. 13.18 The equation of motion of a single degree of freedom system is given by $ # 2x + 0.8x + 1.6x = 0 # with initial conditions x102 = - 1 and x102 = 2. (a) Plot the graph x(t) versus t for 0 … t … 10. (b) Plot a trajectory in the phase plane. 13.19 Find the equilibrium position and plot the trajectories in the neighborhood of the equilibrium position corresponding to the following equation: $ # x + 0.1 1x2 - 12 x + x = 0 13.20 Obtain the phase trajectories for a system governed by the equation # $ x + 0.4x + 0.8x = 0 # with the initial conditions x102 = 2 and x102 = 1 using the method of isoclines. 13.21 Plot the phase-plane trajectories for the following system: # $ x + 0.1x + x = 5 # The initial conditions are x102 = x102 = 0. 13.22 A single degree of freedom system is subjected to Coulomb friction so that the equation of motion is given by # x $ x + f # + v2n x = 0 ƒxƒ Construct the phase plane trajectories of the system using the initial conditions x102 = # 101f/v2n2 and x102 = 0. 13.23 The equation of motion of a simple pendulum subject to viscous damping can be expressed as $ # u + c u + sin u = 0 # If the initial conditions are u102 = u0 and u102 = 0, show that the origin in the phase plane diagram represents (a) a stable focus for c 7 0 and (b) an unstable focus for c 6 0. 13.24 The equation of motion of a simple pendulum, subjected to external force, is given by $ # u + 0.5 u + sin u = 0.8 Find the nature of singularity at u = sin -110.82. RaoCh13ff.qxd 10.06.08 13:57 Page 969 PROBLEMS 969 13.25 The phase plane equation of a single degree of freedom system is given by - cy - 1x - 0.1x32 dy = dx y Investigate the nature of singularity at 1x, y2 = 10, 02 for c 7 0. 13.26 Identify the singularity and find the nature of solution near the singularity for van der Pol’s equation: $ # x - a 11 - x22 x + x = 0 13.27 Identify the singularity and investigate the nature of solution near the singularity for an undamped system with a hard spring: $ x + v2n 11 + k2x22 x = 0 13.28 Solve Problem 13.27 for an undamped system with a soft spring: $ x + v2n 11 - k2x22 x = 0 13.29 Solve Problem 13.27 for a simple pendulum: $ u + v2n sin u = 0 13.30 Determine the eigenvalues and eigenvectors of the following equations: # # y = x + 3y a. x = x - y, # # b. x = x + y, y = 4x + y 13.31 Find the trajectories of the system governed by the equations # x = x - 2y, # y = 4x - 5y 13.32 Find the trajectories of the system governed by the equations # x = x - y, # y = x + 3y 13.33 Find the trajectories of the system governed by the equations # x = 2x + y, # y = - 3x - 2y 13.34 Using Lindstedt’s perturbation method, find the solution of the van der Pol’s equation, Eq. (13.143). 13.35 Verify that the following equation exhibits chaotic behavior: xn + 1 = k xn 11 - xn2 Hint: Give values of 3.25, 3.5 and 3.75 to k and observe the sequence of values generated with x1 = 0.5. RaoCh13ff.qxd 970 10.06.08 13:57 Page 970 CHAPTER 13 NONLINEAR VIBRATION 13.36 Verify that the following equation exhibits chaotic behavior: xn + 1 = 2.0 xn 1xn - 12 Hint: Observe the sequence of values generated using x1 = 1.001, 1.002 and 1.003. 13.37 Using MATLAB, solve the simple pendulum equations, Eqs. (E.1) to (E.3) given in Example 13.6, for the following data: # a. v0 = 0.1, u102 = 0.01, u102 = 0 # b. v0 = 0.1, u102 = 0.01, u102 = 10 13.38 Using MATLAB, find the solution of the nonlinearly damped system, Eq. (E.1) of Example 13.7, for the following data: m = 10, c = 0.1, k = 4000, F0 = 200, v = 20, x102 = 0.5, # x102 = 1.0. 13.39 Using MATLAB, find the solution of a nonlinear single degree of freedom system governed by Eq. (E.1) of Example 13.8 under a pulse load for the following data: m = 10, k1 = 4000, # k2 = 1000, F0 = 1000, t0 = 5, x102 = 0.05, x102 = 5. # $ 13.40 Solve the equation of motion x + 0.5x + x + 1.2x3 = 1.8 cos 0.4t, using the Runge-Kutta # method with ¢t = 0.05, tmax = 5.0, and x0 = x0 = 0. Plot the variation of x with t. Use Program18.m for the solution. 13.41 Find the time variation of the angular displacement of a simple pendulum (i.e., the solution of # Eq. 13.5) for g/l = 0.5, using the initial conditions u0 = 45° and u0 = 0. Use the RungeKutta method given in Program18.m. 13.42 In the static firing test of a rocket, the rocket is anchored to a rigid wall by a nonlinear springdamper system and fuel is burnt to develop a thrust, as shown in Fig. 13.34. The thrust acting on the rocket during the time period 0 … t … t0 is given by F = m0v, where m0 is the constant rate at which fuel is burned and v is the velocity of the jet stream. The initial mass of the rocket is M, so that its mass at any time t is given by m = M - m0t, 0 … t … t0. The data # # is: spring force = 8 * 105x + 6 * 103x3 N, damping force = 10x + 20x2 N, m0 = 10 kg/s, v = 2000 m/s, M = 2000 kg, and t0 = 100 s. (a) Using the Runge-Kutta method of numerical integration, derive the equation of motion of the rocket and (b) Find the variation of the displacement of the rocket. Use Program18.m. x(t) k F c v FIGURE 13.34 13.43 Using Program18.cpp, solve Problem 13.40. 13.44 Using Program18.cpp, solve Problem 13.41. 13.45 Using Program18.cpp, solve Problem 13.42. 13.46 Using PROGRAM18.F, solve Problem 13.40. RaoCh13ff.qxd 10.06.08 13:57 Page 971 DESIGN PROJECTS 971 13.47 Using PROGRAM18.F, solve Problem 13.41. 13.48 Using PROGRAM18.F, solve Problem 13.42. 13.49 Write a computer program for finding the period of vibration corresponding to Eq. (13.14). Use a suitable numerical integration procedure. Using this program, find the solution of Problem 13.41. DESIGN PROJECTS 13.50 In some periodic vibratory systems, external energy is supplied to the system over part of a period and dissipated within the system in another part of the period. Such periodic oscillations are known as relaxation oscillations. Van der Pol [13.22] indicated several instances of occurrence of relaxation oscillations such as a pneumatic hammer, the scratching noise of a knife on a plate, the squeaking of a door, and the fluctuation of populations of animal species. Many relaxation oscillations are governed by van der Pol’s equation: $ # x - a 11 - x22 x + x = 0 (E.1) a. Plot the phase plane trajectories for three values of a: a = 0.1, a = 1, and a = 10. Use # # the initial conditions (i) x102 = 0.5, x102 = 0 and (ii) x102 = 0 and x102 = 5. b. Solve Eq. (E.1) using the fourth-order Runge-Kutta method using the initial conditions stated in (a) for a = 0.1, a = 1, and a = 10. 13.51 A machine tool is mounted on two nonlinear elastic mounts, as shown in Fig. 13.35. The equations of motion, in terms of the coordinates x(t) and u1t2, are given by $ mx + k11 1x - l1 u2 + k12 1x - l1 u23 + k21 1x + l2 u2 + k22 1x + l2 u23 = 0 G x (t) x(t) (t) G l1 k1 FIGURE 13.35 l2 k2 (E.1) RaoCh13ff.qxd 972 10.06.08 13:57 Page 972 CHAPTER 13 NONLINEAR VIBRATION $ J0 u - k11 1x - l1 u2 l1 - k12 1x - l1 u23 l1 + k21 1x + l2 u) l2 + k22 1x + l2 u23 l2 = 0 (E.2) where m is the mass and J0 is the mass moment of inertia about G of the machine tool. Using the Runge-Kutta method, find x(t) and u1t2 for the following data: m = 1000 kg, J0 = 2500 kg-m2, l1 = 1 m, l2 = 1.5 m, k1 = 40x1 + 10x31 kN/m, and k2 = 50x2 + 5x32 kN/m. RaoCh14ff.qxd 10.06.08 14:00 Page 973 Karl Friedrich Gauss (1777–1855) was a German mathematician, astronomer, and physicist. Gauss, Archimedes, and Newton are considered to be in a class by themselves among the great mathematicians. Although Gauss was born in a poor family, his extraordinary intelligence and genius in childhood inspired the Duke of Brunswick to pay all the expenses of Gauss during his entire education. In 1795, Gauss entered the University of Göttingen to study mathematics and began to keep his scientific diary. It was found after his death that his diary contained theories which were rediscovered and published by others. He moved to the University of Helmstedt in 1798 from which he received his doctor’s degree in 1799. He published his most famous work, Disquisitiones Arithmeticae (Arithmetical Researches), in 1801. After that, Gauss was made the director of Göttingen Observatory and also broadened his activities to include the mathematical and practical aspects of astronomy, geodesy, and electromagnetism. The instrument used for measuring magnetic field strength is called the “Gaussmeter.” He invented the method of least squares and the law of normal distribution (Gaussian distribution) which is widely used in probability and random vibration. (Photo courtesy of Dirk J. Struik, A Concise History of Mathematics, 2nd ed., Dover Publications, New York, 1948.) C H A P T E R 1 4 Random Vibration 14.1 Introduction If vibrational response characteristics such as displacement, acceleration, and stress are known precisely as functions of time, the vibration is known as deterministic vibration. This implies a deterministic system (or structure) and a deterministic loading (or excitation); deterministic vibration exists only if there is perfect control over all the variables that influence the structural properties and the loading. In practice, however, there are many processes and phenomena whose parameters cannot be precisely predicted. Such processes are therefore called random processes [14.1–14.4]. An example of a random process is pressure fluctuation at a particular point on the surface of an aircraft flying in air. If several records of these pressure fluctuations are taken under the same flight speed, altitude, and load factor, they might look as indicated in Fig. 14.1. The records are not identical even though the measurements are taken under seemingly identical conditions. Similarly, a building subjected to ground acceleration due to an earthquake, a water tank under wind loading, and a car running on a rough road represent random processes. An elementary treatment of random vibration is presented in this chapter. 973 RaoCh14ff.qxd 974 10.06.08 14:00 Page 974 CHAPTER 14 RANDOM VIBRATION FIGURE 14.1 Ensemble of a random process (x 1i21t2 is the ith sample function of the ensemble). 14.2 14.3 Random Variables and Random Processes Most phenomena in real life are nondeterministic. For example, the tensile strength of steel and the dimensions of a machined part are nondeterministic. If many samples of steel are tested, their tensile strengths will not be the same—they will fluctuate about a mean or average value. Any quantity, like the tensile strength of steel, whose magnitude cannot be precisely predicted, is known as a random variable or a probabilistic quantity. If experiments are conducted to find the value of a random variable x, each experiment will give a number that is not a function of any parameter. For example, if 20 samples of steel are tested, their tensile strengths might be x112 = 284, x122 = 302, x132 = 269, Á , x1202 = 298 N/mm2. Each of these outcomes is called a sample point. If n experiments are conducted, all the n possible outcomes of the random variable constitute what is known as the sample space of the random variable. There are other types of probabilistic phenomena for which the outcome of an experiment is a function of some parameter such as time or a spatial coordinate. Quantities such as the pressure fluctuations shown in Fig. 14.1 are called random processes. Each outcome of an experiment, in the case of a random process, is called a sample function. If n experiments are conducted, all the n possible outcomes of a random process constitute what is known as the ensemble of the process [14.5]. Notice that if the parameter t is fixed at a particular value t1, x1t12 is a random variable whose sample points are given by x1121t12, x1221t12, Á , x1n21t12. Probability Distribution Consider a random variable x such as the tensile strength of steel. If n experimental values of x are available as x1, x2, Á , xn, the probability of realizing the value of x smaller than RaoCh14ff.qxd 10.06.08 14:00 Page 975 14.3 PROBABILITY DISTRIBUTION 975 some specified value x ' can be found as Prob1x … x '2 = n ' n (14.1) where n ' denotes the number of xi values which are less than or equal to x ' . As the number of experiments n : q , Eq. (14.1) defines the probability distribution function of x, P(x): P1x2 = lim n ' n: q n (14.2) The probability distribution function can also be defined for a random time function. For this, we consider the random time function shown in Fig. 14.2. During a fixed time span t, the time intervals for which the value of x(t) is less than x' are denoted as ¢t1, ¢t2, ¢t3, and ¢ 4. Thus the probability of realizing x(t) less than or equal to x ' is given by Prob [x1t2 … x' ] = 1 ¢ti t a i (14.3) As t : q , Eq. (14.3) gives the probability distribution function of x(t): P1x2 = lim t: q 1 ¢ti t a i (14.4) If x(t) denotes a physical quantity, the magnitude of x(t) will always be a finite number, so Prob[x1t2 6 - q ] = P1 - q 2 = 0 (impossible event), and Prob[x1t2 6 q ] = P1 q 2 = 1 (certain event). The typical variation of P(x) with x is shown in Fig. 14.3(a). The function P(x) is called the probability distribution function of x. The derivative of P(x) with respect to x is known as the probability density function and is denoted as p(x). Thus p1x2 = FIGURE 14.2 dP1x2 dx = lim Random time function. ¢x : 0 P1x + ¢x2 - P1x2 ¢x (14.5) RaoCh14ff.qxd 976 10.06.08 14:00 Page 976 CHAPTER 14 RANDOM VIBRATION FIGURE 14.3 Probability distribution and density functions. where the quantity P1x + ¢x2 - P1x2 denotes the probability of realizing x(t) between the values x and x + ¢x. Since p(x) is the derivative of P(x), we have x As P1 q 2 = 1, Eq (14.6) gives P1x2 = P1 q 2 = L- q p1x¿2 dx¿ (14.6) p1x¿2 dx¿ = 1 (14.7) q L- q which means that the total area under the curve of p(x) is unity. 14.4 Mean Value and Standard Deviation If f(x) denotes a function of the random variable x, the expected value of f(x), denoted as mf or E [ f(x)] or f1x2, is defined as q mf = E[f1x2] = f1x2 = L- q f1x2p1x2 dx (14.8) If f1x2 = x, Eq. (14.8) gives the expected value, also known as the mean value of x: q mx = E[x] = x = L- q xp1x2 dx (14.9) RaoCh14ff.qxd 10.06.08 14:00 Page 977 977 14.4 MEAN VALUE AND STANDARD DEVIATION Similarly, if f1x2 = x 2, we get the mean square value of x: q 2 L- q 2 mx = E[x ] = x = 2 x2p1x2 dx (14.10) The variance of x, denoted as s2x, is defined as the mean square value of x about the mean, q s2x 2 = E[1x - x2 ] = L- q 1x - x22 p1x2 dx = 1x22 - 1x22 (14.11) The positive square root of the variance, s1x2, is called the standard deviation of x. Probabilistic Characteristics of Eccentricity of a Rotor EXAMPLE 14.1 The eccentricity of a rotor (x), due to manufacturing errors, is found to have the following distribution: p1x2 = e kx2, 0, 0 … x … 5 mm elsewhere (E.1) where k is a constant. Find (a) the mean, standard deviation, and the mean square value of the eccentricity and (b) the probability of realizing x less than or equal to 2 mm. Solution: The value of k in Eq. (E.1) can be found by normalizing the probability density function: q L- q p1x2 dx = That is, ka That is, L0 5 x3 5 b = 1 3 0 k = 3 125 (a) The mean value of x is given by Eq. (14.9): x = L0 kx2 dx = 1 5 p1x2x dx = ka x4 5 b = 3.75 mm 4 0 The standard deviation of x is given by Eq. (14.11): s2x = = L0 5 1x - x22 p1x2 dx L0 5 1x2 + x2 - 2xx2 p1x2 dx (E.2) (E.3) RaoCh14ff.qxd 978 10.06.08 14:00 Page 978 CHAPTER 14 RANDOM VIBRATION = L0 = ka = ka ‹ 5 kx4 dx - 1x22 x5 5 b - 1x22 5 0 3125 b - 13.7522 = 0.9375 5 sx = 0.9682 mm (E.4) The mean square value of x is x 2 = ka 3125 b = 15 mm2 5 (E.5) (b) Prob [x … 2] = L0 = ka 2 p1x2 dx = k L0 2 x2 dx 8 x3 2 b = = 0.064 3 0 125 (E.6) ■ 14.5 Joint Probability Distribution of Several Random Variables When two or more random variables are being considered simultaneously, their joint behavior is determined by a joint probability distribution function. For example, while testing the tensile strength of steel specimens, we can obtain the values of yield strength and ultimate strength in each experiment. If we are interested in knowing the relation between these two random variables, we must know the joint probability density function of yield strength and ultimate strength. The probability distributions of single random variables are called univariate distributions; the distributions that involve two random variables are called bivariate distributions. In general, if a distribution involves more than one random variable, it is called a multivariate distribution. The bivariate density function of the random variables x1 and x2 is defined by p1x1, x22 dx1 dx2 = Prob [x1 6 x1œ 6 x1 + dx1, x2 6 x2œ 6 x2 + dx2] (14.12) that is, the probability of realizing the value of the first random variable between x1 and x1 + dx1 and the value of the second random variable between x2 and x2 + dx2. The RaoCh14ff.qxd 10.06.08 14:00 Page 979 979 14.5 JOINT PROBABILITY DISTRIBUTION OF SEVERAL RANDOM VARIABLES joint probability density function has the property that q q L- q L- q p1x1, x22 dx1 dx2 = 1 (14.13) P1x1, x22 = Prob [x1œ 6 x1, x 2œ 6 x2] The joint distribution function of x1 and x2 is x1 = x2 L- q L- q p1x1œ , x 2œ 2 dx1œ dx 2œ (14.14) The marginal or individual density functions can be obtained from the joint probability density function as q p1x2 = L- q p1x, y2 dy (14.15) p1x, y2 dx (14.16) q p1y2 = L- q The variances of x and y can be determined as s2x s2y = E[1x - mx2 ] = 2 = E[1y - my2 ] = 2 q L- q q L- q 1x - mx22p1x2 dx 1y - my22p1y2 dy (14.17) (14.18) The covariance of x and y, sxy, is defined as the expected value or average of the product of the deviations from the respective mean values of x and y. It is given by sxy = E[1x - mx21y - my2] = = = q q q q L- q L- q L- q L- q q - mx q q L- q L- q 1x - mx21y - my2p1x, y2 dx dy 1xy - xmy - ymx + mxmy2p1x, y2 dx dy q xyp1x, y2 dx dy - my L- q L- q xp1x, y2 dx dy q q L- q L- q q yp1x, y2 dx dy + mxmy = E[xy] - mxmy q L- q L- q p1x, y2 dx dy (14.19) RaoCh14ff.qxd 980 10.06.08 14:00 Page 980 CHAPTER 14 RANDOM VIBRATION The correlation coefficient between x and y, rxy, is defined as the normalized covariance: rxy = sxy sxsy (14.20) It can be seen that the correlation coefficient satisfies the relation -1 … rxy … 1. 14.6 Correlation Functions of a Random Process If t1, t2, Á are fixed values of t, we use the abbreviations x1, x2, Á to denote the values of x(t) at t1, t2, Á , respectively. Since there are several random variables x1, x2, Á , we form the products of the random variables x1, x2, Á (values of x(t) at different times) and average the products over the set of all possibilities to obtain a sequence of functions: K1t1, t22 = E[x1t12x1t22] = E[x1x2] K1t1, t2, t32 = E[x1t12x1t22x1t32] = E[x1x2x3] (14.21) and so on. These functions describe the statistical connection between the values of x(t) at different times t1, t2, Á and are called correlation functions [14.6, 14.7]. Autocorrelation Function. The mathematical expectation of x1x2—the correlation function K1t1, t22—is also known as the autocorrelation function, designated as R1t1, t22. Thus R1t1, t22 = E[x1x2] (14.22) If the joint probability density function of x1 and x2 is known to be p1x1, x22, the autocorrelation function can be expressed as R1t1, t22 = q q L- q L- q x1x2p1x1, x22 dx1 dx2 (14.23) Experimentally, we can find R1t1, t22 by taking the product of x1i21t12 and x1i21t22 in the ith sample function and averaging over the ensemble: R1t1, t22 = 1 n 1i2 x 1t12x1i21t22 n ia =1 (14.24) where n denotes the number of sample functions in the ensemble (see Fig. 14.4). If t1 and t2 are separated by t (with t1 = t and t2 = t + t), we have R1t + t2 = E[x1t2x1t + t2]. RaoCh14ff.qxd 10.06.08 14:00 Page 981 14.7 STATIONARY RANDOM PROCESS FIGURE 14.4 14.7 981 Ensemble of a random process. Stationary Random Process A stationary random process is one for which the probability distributions remain invariant under a shift of the time scale; the family of probability density functions applicable now also applies five hours from now or 500 hours from now. Thus the probability density function p1x12 becomes a universal density function p(x), independent of time. Similarly, the joint density function p1x1, x22, to be invariant under a shift of the time scale, becomes a function of t = t2 - t1, but not a function of t1 or t2 individually. Thus p1x1, x22 can be written as p1t, t + t2. The expected value of stationary random process x(t) can be written as E[x1t12] = E[x1t1 + t2] for any t (14.25) and the autocorrelation function becomes independent of the absolute time t and will depend only on the separation time t: R1t1, t22 = E[x1x2] = E[x1t2x1t + t2] = R1t2 for any t (14.26) where t = t2 - t1. We shall use subscripts to R to identify the random process when more than one random process is involved. For example, we shall use Rx1t2 and Ry1t2 to denote the autocorrelation functions of the random processes x(t) and y(t), respectively. The autocorrelation function has the following characteristics [14.2, 14.4]: RaoCh14ff.qxd 982 10.06.08 14:00 Page 982 CHAPTER 14 RANDOM VIBRATION 1. If t = 0, R1t2 gives the mean square value of x(t): R102 = E[x2] (14.27) 2. If the process x(t) has a zero mean and is extremely irregular, as shown in Fig. 14.5(a), its autocorrelation function R1t2 will have small values except at t = 0, as indicated in Fig. 14.5(b). 3. If x1t2 M x1t + t2, the autocorrelation function R1t2 will have a constant value as shown in Fig. 14.6. 4. If x(t) is stationary, its mean and standard deviations will be independent of t: E[x1t2] = E[x1t + t2] = m (14.28) sx1t2 = sx1t + t2 = s (14.29) and The correlation coefficient, r, of x(t) and x1t + t2 can be found as r = = = FIGURE 14.5 E[5x1t2 - m65x1t + t2 - m6] s2 E[x1t2x1t + t2] - mE[x1t + t2] - mE[x1t2] + m2 s2 R1t2 - m2 s2 Autocorrelation function. (14.30) RaoCh14ff.qxd 10.06.08 14:00 Page 983 14.7 STATIONARY RANDOM PROCESS FIGURE 14.6 983 Constant values of function. that is, R1t2 = rs2 + m2 (14.31) Since ƒ r ƒ … 1, Eq. (14.31) shows that -s2 + m2 … R1t2 … s2 + m2 (14.32) This shows that the autocorrelation function will not be greater than the mean square value, E[x2] = s2 + m2. 5. Since R1t2 depends only on the separation time t and not on the absolute time t for a stationary process, R1t2 = E[x1t2x1t + t2] = E[x1t2x1t - t2] = R1-t2 (14.33) Thus R1t2 is an even function of t. 6. When t is large 1t : q 2, there will not be a coherent relationship between the two values x(t) and x1t + t2. Hence the correlation coefficient r : 0 and Eq. (14.31) gives R1t : q 2 : m2 (14.34) RaoCh14ff.qxd 984 10.06.08 14:00 Page 984 CHAPTER 14 RANDOM VIBRATION FIGURE 14.7 Autocorrelation function. A typical autocorrelation function is shown in Fig. 14.7. Ergodic Process. An ergodic process is a stationary random process for which we can obtain all the probability information from a single sample function and assume that it is applicable to the entire ensemble. If x1i21t2 represents a typical sample function of duration T, the averages can be computed by averaging with respect to time along x1i21t2. Such averages are called temporal averages. By using the notation 8x1t29 to represent the temporal average of x(t) (the time average of x), we can write E[x] = 8x1t29 = lim 1 x1i21t2 dt T: q T L - T/2 T/2 (14.35) where x1i21t2 has been assumed to be defined from t = - T/2 to t = T/2 with T : q (T very large). Similarly, E[x2] = 8x21t29 = lim 1 [x1i21t2]2 dt q T: T L - T/2 (14.36) 1 x1i21t2x1i21t + t2 dt T: q T L - T/2 (14.37) T/2 and R1t2 = 8x1t2x1t + t29 = lim 14.8 T/2 Gaussian Random Process The most commonly used distribution for modeling physical random processes is called the Gaussian or normal random process. The Gaussian process has a number of RaoCh14ff.qxd 10.06.08 14:00 Page 985 14.8 GAUSSIAN RANDOM PROCESS 985 remarkable properties that permit the computation of the random vibration characteristics in a simple manner. The probability density function of a Gaussian process x(t) is given by p1x2 = 1 x - xq 2 1 e- 2 A s B 12psx (14.38) x where x and sx denote the mean value and standard deviation of x. The mean 1x2 and standard deviation 1sx2 of x(t) vary with t for a nonstationary process but are constants (independent of t) for a stationary process. A very important property of the Gaussian process is that the forms of its probability distributions are invariant with respect to linear operations. This means that if the excitation of a linear system is a Gaussian process, the response is generally a different random process, but still a normal one. The only changes are that the magnitude of the mean and standard deviations of the response are different from those of the excitation. The graph of a Gaussian probability density function is a bell-shaped curve, symmetric about the mean value; its spread is governed by the value of the standard deviation, as shown in Fig. 14.8. By defining a standard normal variable z as z = x - x sx (14.39) 1 2 1 e-2z 12p (14.40) Eq. (14.38) can be expressed as p1x2 = The probability of x(t) lying in the interval -cs and +cs where c is any positive number can be found, assuming x = 0: cs Prob [ -cs … x1t2 … cs] = FIGURE 14.8 density function. Gaussian probability L-cs 1x 1 e - 2 s dx 12ps 2 2 (14.41) RaoCh14ff.qxd 986 10.06.08 14:00 Page 986 CHAPTER 14 RANDOM VIBRATION FIGURE 14.9 Graphical representation of Prob[- cs … x1t2 … cs]. The probability of x(t) lying outside the range ⫿cs is one minus the value given by Eq. (14.41). This can also be expressed as q 1x 2 Prob [ ƒ x1t2 ƒ 7 cs] = e - 2 s dx 12ps Lcs 2 2 (14.42) The integrals in Eqs. (14.41) and (14.42) have been evaluated numerically and tabulated [14.5]; some typical values are indicated in the following table (see Fig. 14.9 also). Value of c Prob[-cs … x1t2 … cs] Prob[ ƒ x1t2 ƒ 7 cs] 14.9 1 2 3 0.6827 0.3173 0.9545 0.0455 0.9973 0.0027 4 0.999937 0.000063 Fourier Analysis 14.9.1 Fourier Series We saw in Chapter 1 that any periodic function x(t), of period t, can be expressed in the form of a complex Fourier series aqcne q x1t2 = inv0t (14.43) n=- where v0 is the fundamental frequency given by v0 = 2p t (14.44) RaoCh14ff.qxd 10.06.08 14:00 Page 987 14.9 FOURIER ANALYSIS 987 and the complex Fourier coefficients cn can be determined by multiplying both sides of Eq. (14.43) with e -imv0t and integrating over one time period: a q t/2 L- t/2 x1t2e -imv0t dt = t/2 - t/2 n=- q L cnei1n - m2v0t dt aqcn q t/2 [cos1n - m2v0t + i sin1n - m2v0t] dt L- t/2 Equation (14.45) can be simplified to obtain (see Problem 14.27) = (14.45) n=- t/2 1 x1t2e -inv0t dt cn = t L- t/2 (14.46) Equation (14.43) shows that the function x(t) of period t can be expressed as a sum of an infinite number of harmonics. The harmonics have amplitudes given by Eq. (14.46) and frequencies which are multiples of the fundamental frequency v0. The difference between any two consecutive frequencies is given by vn + 1 - vn = 1n + 12v0 - nv0 = ¢v = 2p = v0 t (14.47) Thus the larger the period t, the denser the frequency spectrum becomes. Equation (14.46) shows that the Fourier coefficients cn are, in general, complex numbers. However, if x(t) is a real and even function, then cn will be real. If x(t) is real, the integrand of cn in Eq. (14.46) can also be identified as the complex conjugate of that of c- n. Thus cn = c*- n (14.48) The mean square value of x(t)—that is, the time average of the square of the function x(t)—can be determined as 2 1 1 x 1t2 = x21t2 dt = a a cneinv0t b dt t L- t/2 t L- t/2 n = - q t/2 t/2 q 2 2 -1 1 a a cneinv0t + c0 + a cneinv0t b dt t L- t/2 n = - q n=1 q t/2 = 2 1 e a 1cneinv0t + c…n e -inv0t2 + c0 f dt t L- t/2 n = 1 t/2 = 1 2cnc…n + c20 f dt e t L- t/2 na =1 t/2 = q q = c20 + a 2 ƒ cn ƒ 2 = q n=1 2 a ƒ cn ƒ q n=- q (14.49) RaoCh14ff.qxd 988 10.06.08 14:00 Page 988 CHAPTER 14 RANDOM VIBRATION Thus the mean square value of x(t) is given by the sum of the squares of the absolute values of the Fourier coefficients. Equation (14.49) is known as Parseval’s formula for periodic functions [14.1]. Complex Fourier Series Expansion EXAMPLE 14.2 Find the complex Fourier series expansion of the function shown in Fig. 14.10(a). Solution: The given function can be expressed as t t b, - … t … 0 a 2 x1t2 = µ t t Aa 1 - b, 0 … t … a 2 Aa 1 + (E.1) where the period 1t2 and the fundamental frequency 1v02 are given by t = 2a and v0 = p 2p = t a FIGURE 14.10 Complex Fourier series representation. (E.2) RaoCh14ff.qxd 10.06.08 14:00 Page 989 14.9 FOURIER ANALYSIS 989 The Fourier coefficients can be determined as t/2 cn = 1 x1t2e -inv0t dt t L- t/2 t 1 t c A a1 - b e -inv0t dt d Aa 1 + b e -inv0t dt + t L- t/2 a a L0 t/2 0 = Using the relation L tekt dt = cn can be evaluated as cn = ekt 1kt - 12 k2 0 0 A e -inv0t 1 A c + e -inv0t ` e [-inv0t - 1] f ` 2 t - inv0 a 1 -inv02 - t/2 - t/2 + t/2 t/2 A e -inv0t A e -inv0t ` e [-inv0t - 1] f ` d 2 - inv0 a 1 -inv02 0 0 (E.3) (E.4) (E.5) This equation can be reduced to cn = 2A 1 A 1 in p A -in p A 1 -in p 1 A in p c e + e e e t inv0 a n2v20 inv0 a n2v20 a n2v20 + A 1 A 1 1in p2ein p 1in p2e -in p d a n2v20 a n2v20 (E.6) Noting that ein p or 1, e -in p = c - 1, 1, n = 0 n = 1, 3, 5, Á n = 2, 4, 6, Á (E.7) Eq. (E.6) can be simplified to obtain A , 2 cn = e a 4A b atn2v20 0, n = 0 = 2A , n2p2 n = 1, 3, 5, Á (E.8) n = 2, 4, 6, Á The frequency spectrum is shown in Fig. 14.10(b). ■ RaoCh14ff.qxd 990 10.06.08 14:00 Page 990 CHAPTER 14 RANDOM VIBRATION FIGURE 14.11 Nonperiodic function. 14.9.2 Fourier Integral A nonperiodic function, such as the one shown by the solid curve in Fig. 14.11, can be treated as a periodic function having an infinite period 1t : q 2. The Fourier series expansion of a periodic function is given by Eqs. (14.43), (14.44) and (14.46): inv0t a cne q x1t2 = (14.50) n= - q with 2p t (14.51) 1 x1t2e -inv0t dt t L- t/2 (14.52) v0 = and t/2 cn = As t : q , the frequency spectrum becomes continuous and the fundamental frequency becomes infinitesimal. Since the fundamental frequency v0 is very small, we can denote it as ¢v, nv0 as v, and rewrite Eq. (14.52) as q t/2 t: q L - t/2 x1t2e -i vt dt = lim tcn = lim t: q By defining X1v2 as X1v2 = lim 1tcn2 = t: q L- q x1t2e -i vt dt (14.53) q L- q x1t2e -i v t dt (14.54) RaoCh14ff.qxd 10.06.08 14:00 Page 991 14.9 FOURIER ANALYSIS 991 we can express x(t) from Eq. (14.50) as i vt 2pt aqcne 2pt q x1t2 = lim t: q n= - aq1cnt2e q = lim i vt t: q n= q = a 1 X1v2ei vt dv 2p L- q 2p 1 b t 2p (14.55) This equation indicates the frequency decomposition of the nonperiodic function x(t) in a continuous frequency domain, similar to Eq. (14.50) for a periodic function in a discrete frequency domain. The equations q 1 x1t2 = X1v2ei vt dv 2p L- q (14.56) and q X1v2 = L- q x1t2e -i vt dt (14.57) are known as the (integral) Fourier transform pair for a nonperiodic function x(t), similar to Eqs. (14.50) and (14.52) for a periodic function x(t) [14.9, 14.10]. The mean square value of a nonperiodic function x(t) can be determined from Eq. (14.49): 1 x21t2dt = t L- t/2 t/2 = = 2 a ƒ cn ƒ q n= - q tv … 0 aqcncn tv0 = n=q … aqcncn q n=- q v0 1 1tcn21c…n t2 t n =a 2p q - tv0 2p b ta t (14.58) Since tcn : X1v2, tc…n : X*1v2, and v0 : dv as t : q , Eq. (14.58) gives the mean square value of x(t) as t/2 ƒ X1v2 ƒ 2 1 x21t2 dt = dv t: q t L - t/2 L- q 2pt x21t2 = lim q Equation (14.59) is known as Parseval’s formula for nonperiodic functions [14.1]. (14.59) RaoCh14ff.qxd 992 10.06.08 14:00 Page 992 CHAPTER 14 RANDOM VIBRATION Fourier Transform of a Triangular Pulse EXAMPLE 14.3 Find the Fourier transform of the triangular pulse shown in Fig. 14.12(a). Solution: The triangular pulse can be expressed as x1t2 = c A a1 - 0, ƒtƒ b, a ƒtƒ … a otherwise (E.1) The Fourier transform of x(t) can be found, using Eq. (14.57), as q X1v2 = L- q 0 = L- q A a1 Aa 1 + ƒ t ƒ -i vt be dt a t t -i vt Aa 1 - be -i vt dt be dt + a a L0 q FIGURE 14.12 Fourier transform of a triangular pulse. (E.2) RaoCh14ff.qxd 10.06.08 14:00 Page 993 14.10 POWER SPECTRAL DENSITY 993 Since x1t2 = 0 for ƒ t ƒ 7 0, Eq. (E.2) can be expressed as Aa1 + 0 X1v2 = L-a = a t t -i vt be dt + Aa 1 - b e -ivt dt a a L0 a 0 0 A e -i vt A [ivt 1] f ` b e -i vt ` + e -iv a 1 -iv22 -a -a + a a a A e -i vt A e [ivt 1] f ` b e -i vt ` -iv a 1- iv22 0 0 (E.3) Equation (E.3) can be simplified to obtain 2A A A + ei va a b + e -i va a b 2 2 av av av2 2A A A = 1cos va + i sin va2 1cos va - i sin va2 av2 av2 av2 X1v2 = = 2A 4A 2 va 11 - cos va2 = sin a b 2 2 av av2 (E.4) Equation (E.4) is plotted in Fig. 14.12(b). Notice the similarity of this figure with the discrete Fourier spectrum shown in Fig. 14.10(b). ■ 14.10 Power Spectral Density The power spectral density S1v2 of a stationary random process is defined as the Fourier transform of R1t2/2p q 1 S1v2 = R1t2e -i vt dt 2p L- q (14.60) so that q R1t2 = L- q S1v2ei vt dv (14.61) Equations (14.60) and (14.61) are known as the Wiener-Khintchine formulas [14.1]. The power spectral density is more often used in random vibration analysis than the autocorrelation function. The following properties of power spectral density can be observed: 1. From Eqs. (14.27) and (14.61), we obtain q R102 = E[x2] = L- q S1v2 dv (14.62) RaoCh14ff.qxd 994 10.06.08 14:00 Page 994 CHAPTER 14 RANDOM VIBRATION FIGURE 14.13 Typical power spectral density function. If the mean is zero, the variance of x(t) is given by q s2x = R102 = L- q S1v2 dv (14.63) If x(t) denotes the displacement, R(0) represents the average energy. From Eq. (14.62), it is clear that S1v2 represents the energy density associated with the frequency v. Thus S1v2 indicates the spectral distribution of energy in a system. Also, in electrical circuits, if x(t) denotes random current, then the mean square value indicates the power of the system (when the resistance is unity). This is the origin of the term power spectral density. 2. Since R1t2 is an even function of t and real, S1v2 is also an even and real function of v. Thus S1- v2 = S1v2. A typical power spectral density function is shown in Fig. 14.13. 3. From Eq. (14.62), the units of S1v2 can be identified as those of x2/unit of angular frequency. It can be noted that both negative and positive frequencies are counted in Eq. (14.62). In experimental work, for convenience, an equivalent one-sided spectrum Wx1f2 is widely used [14.1, 14.2].1 The spectrum Wx1f2 is defined in terms of linear frequency (i.e., cycles per unit time) and only the positive frequencies are counted. The relationship between Sx1v2 and Wx1f2 can be seen with reference to Fig. 14.14. The differential frequency dv in Fig. 14.14(a) corresponds to the differential frequency df = dv/2p in Fig. 14.14(b). Since Wx1f2 is the equivalent spectrum defined over positive values of f only, we have q 2 E[x ] = L- q Sx1v2 dv K L0 q Wx1f2 df (14.64) When several random processes are involved, a subscript is used to identify the power spectral density function (or simply the spectrum) of a particular random process. Thus Sx1v2 denotes the spectrum of x(t). 1 RaoCh14ff.qxd 10.06.08 14:00 Page 995 14.11 WIDE-BAND AND NARROW-BAND PROCESSES 995 FIGURE 14.14 Two- and one-sided spectrum. In order to have the contributions of the frequency bands dv and df to the mean square value to be same, the shaded areas in both Figs. 14.14(a) and (b) must be the same. Thus which gives 2Sx1v2 dv = Wx1f2 df Wx1f2 = 2Sx1v2 14.11 dv dv = 2Sx1v2 = 4pSx1v2 df dv/2p (14.65) (14.66) Wide-Band and Narrow-Band Processes A wide-band process is a stationary random process whose spectral density function S1v2 has significant values over a range or band of frequencies that is approximately the same order of magnitude as the center frequency of the band. An example of a wide-band random process is shown in Fig. 14.15. The pressure fluctuations on the surface of a rocket due to acoustically transmitted jet noise or due to supersonic boundary layer turbulence are examples of physical processes that are typically wide-band. A narrow-band random process is a stationary process whose spectral density function S1v2 has significant values only in a range or band of frequencies whose width is small compared to the magnitude of the center frequency of the process. Figure 14.16 shows the sample function and the corresponding spectral density and autocorrelation functions of a narrow-band process. A random process whose power spectral density is constant over a frequency range is called white noise, an analogy with the white light that spans the visible spectrum more or less uniformly. It is called ideal white noise if the band of frequencies v2 - v1 is infinitely wide. Ideal white noise is a physically unrealizable concept, since the mean square value of such a random process would be infinite because the area under the spectrum would be infinite. It is called band-limited white noise if the band of frequencies has finite cut-off RaoCh14ff.qxd 996 10.06.08 14:00 Page 996 CHAPTER 14 RANDOM VIBRATION FIGURE 14.15 Wide-band stationary random process. FIGURE 14.16 Narrow-band stationary random process. frequencies v1 and v2 [14.8]. The mean square value of a band-limited white noise is given by the total area under the spectrum—namely, 2S01v2 - v12, where S0 denotes the constant value of the spectral density function. Autocorrelation and Mean Square Value of a Stationary Process EXAMPLE 14.4 The power spectral density of a stationary random process x(t) is shown in Fig. 14.17(a). Find its autocorrelation function and the mean square value. RaoCh14ff.qxd 10.06.08 14:00 Page 997 14.11 WIDE-BAND AND NARROW-BAND PROCESSES 997 FIGURE 14.17 Autocorrelation function of a stationary process. Solution: (a) Since Sx1v2 is real and even in v, Eq. (14.61) can be rewritten as Rx1t2 = 2 L0 q Sx1v2cos vt dv = 2S0 v2 Lv1 cos vt dv v2 2S0 1 1sin v2t - sin v1t2 = 2S0 a sin vt b ` = t t v1 = 4S0 v1 + v2 v2 - v1 cos t sin t t 2 2 This function is shown graphically in Fig. 14.17(b). (b) The mean square value of the random process is given by q E[x2] = L- q Sx1v2 dv = 2S0 v2 Lv1 dv = 2S01v2 - v12 ■ RaoCh14ff.qxd 998 10.06.08 14:00 Page 998 CHAPTER 14 RANDOM VIBRATION 14.12 Response of a Single Degree of Freedom System The equation of motion for the system shown in Fig. 14.18 is $ # y + 2zvny + v2ny = x1t2 (14.67) where x1t2 = F1t2 , m vn = k , Am z = c , cc and cc = 2km The solution of Eq. (14.67) can be obtained by using either the impulse response approach or the frequency response approach. 14.12.1 Impulse Response Approach Here we consider the forcing function x(t) to be made up of a series of impulses of varying magnitude, as shown in Fig. 14.19(a) (see Section 4.5.2). Let the impulse applied at time t be denoted as x1t2 dt. If y1t2 = h1t - t2 denotes the response to the unit impulse2 excitation d1t - t2, it is called the impulse response function. The total response of the system at time t can be found by superposing the responses to impulses of magnitude x1t2 dt applied at different values of t = t. The response to the excitation x1t2 dt will FIGURE 14.18 Single degree of freedom system. The unit impulse applied at t = t is denoted as 2 x1t2 = d1t - t2 where d1t - t2 is the Dirac delta function with (see Fig. 14.19b) d1t - t2 : q d1t - t2 = 0 q L- q as t : t for all t except at t = t d1t - t2 dt = 1 1area under the curve is unity2 RaoCh14ff.qxd 10.06.08 14:00 Page 999 14.12 RESPONSE OF A SINGLE DEGREE OF FREEDOM SYSTEM 999 FIGURE 14.19 Impulse response approach. be [x1t2 dt]h1t - t2, and the response to the total excitation will be given by the superposition or convolution integral: t y1t2 = L- q x1t2h1t - t2 dt (14.68) RaoCh14ff.qxd 1000 10.06.08 14:00 Page 1000 CHAPTER 14 RANDOM VIBRATION 14.12.2 Frequency Response Approach The transient function x(t) can be expressed in terms of its Fourier transform X1v2 as q 1 X1v2ei vt dv x1t2 = 2p Lv = - q (14.69) Thus x(t) can be considered as the superposition of components of different frequencies v. If we consider the forcing function of unit modulus as i vt x ' 1t2 = e (14.70) y' 1t2 = H1v2ei vt its response can be denoted as (14.71) where H1v2 is called the complex frequency response function (see Section 3.5). Since the actual excitation is given by the superposition of components of different frequencies (Eq. 14.69), the total response of the system can also be obtained by superposition as q y1t2 = H1v2x1t2 = L- q H1v2 1 X1v2ei vt dv 2p q = 1 H1v2X1v2ei vtdv 2p L- q (14.72) If Y1v2 denotes the Fourier transform of the response function y(t), we can express y(t) in terms of Y1v2 as q y1t2 = 1 Y1v2ei vt dv 2p L- q (14.73) Comparison of Eqs. (14.72) and (14.73) yields Y1v2 = H1v2X1v2 14.12.3 Characteristics of the Response Function (14.74) The following characteristics of the response function can be noted: 1. Since h1t - t2 = 0 when t 6 t or t 7 t (i.e., the response before the application of the impulse is zero), the upper limit of integration in Eq. (14.68) can be replaced by q so that q y1t2 = L- q x1t2h1t - t2 dt (14.75) RaoCh14ff.qxd 10.06.08 14:00 Page 1001 1001 14.13 RESPONSE DUE TO STATIONARY RANDOM EXCITATIONS 2. By changing the variable from t to u = t - t, Eq. (14.75) can be rewritten as q y1t2 = L- q x1t - u2h1u2 du (14.76) 3. The superposition integral, Eq. (14.68) or (14.75) or (14.76), can be used to find the response of the system y(t) for any arbitrary excitation x(t) once the impulse-response function of the system h(t) is known. The Fourier integral, Eq. (14.72), can also be used to find the response of the system once the complex frequency response of the system, H1v2, is known. Although the two approaches appear to be different, they are intimately related to one another. To see their inter-relationship, consider the excitation of the system to be a unit impulse d1t2 in Eq. (14.72). By definition, the response is h(t) and Eq. (14.72) gives q y1t2 = h1t2 = 1 X1v2H1v2ei vt dv 2p L- q (14.77) where X1v2 is the Fourier transform of x1t2 = d1t2: q q X1v2 = L- q x1t2e -i vt dt = L- q d1t2e -i vt dt = 1 (14.78) since d1t2 = 0 everywhere except at t = 0 where it has a unit area and e -ivt = 1 at t = 0. Equations (14.77) and (14.78) give q 1 h1t2 = H1v2ei vt dv 2p L- q (14.79) which can be recognized as the Fourier integral representation of h(t) in which H1v2 is the Fourier transformation of h(t): q H1v2 = 14.13 L- q h1t2e -i vt dt (14.80) Response Due to Stationary Random Excitations In the previous section, the relationships between excitation and response were derived for arbitrary known excitations x(t). In this section, we consider similar relationships when the excitation is a stationary random process. When the excitation is a stationary random process, the response will also be a stationary random process [14.15, 14.16]. We consider the relation between the excitation and the response using the impulse response (time domain) as well as the frequency response (frequency domain) approaches. RaoCh14ff.qxd 1002 10.06.08 14:00 Page 1002 CHAPTER 14 RANDOM VIBRATION 14.13.1 Impulse Response Approach Mean Values. The response for any particular sample excitation is given by Eq. (14.76): q y1t2 = L- q x1t - u2h1u2 du (14.81) For ensemble average, we write Eq. (14.81) for every (x, y) pair in the ensemble and then take the average to obtain3 q E[y1t2] = Ec L- q q = L- q x1t - u2h1u2 du d E[x1t - u2]h1u2 du (14.82) Since the excitation is assumed to be stationary, E[x1t2] is a constant independent of t, Eq. (14.82) becomes q E[y1t2] = E[x1t2] L- q h1u2 du (14.83) The integral in Eq. (14.83) can be obtained by setting v = 0 in Eq. (14.80) so that q H102 = L- q h1t2 dt (14.84) Thus a knowledge of either the impulse response function h(t) or the frequency response function H1v2 can be used to find the relationship between the mean values of the excitation and the response. It is to be noted that both E[x(t)] and E[y(t)] are independent of t. Autocorrelation. We can use a similar procedure to find the relationship between the autocorrelation functions of the excitation and the response. For this, we first write q L- q y1t2y1t + t2 = # q L- q x1t - u12h1u12 du1 x1t + t - u22h1u22d u2 x1t - u12x1t + t - u22 L- q L- q # h1u12h1u22 du1 du2 q q = 3 (14.85) In deriving Eq. (14.82), the integral is considered as a limiting case of a summation and hence the average of a sum is treated to be same as the sum of the averages, that is, E[x1 + x2 + Á ] = E[x1] + E[x2] + Á RaoCh14ff.qxd 10.06.08 14:00 Page 1003 14.13 RESPONSE DUE TO STATIONARY RANDOM EXCITATIONS 1003 where u1 and u2 are used instead of u to avoid confusion. The autocorrelation function of y(t) can be found as Ry1t2 = E[y1t2y1t + t2] = = 14.13.2 Frequency Response Approach q q q q L- q L- q E[x1t - u12x1t + t - u22]h1u12h1u22du1 du2 L- q L- q Rx1t + u1 - u22h1u12h1u22du1 du2 (14.86) Power Spectral Density. The response of the system can also be described by its power spectral density, which by definition, is (see Eq. 14.60) Sy1v2 = 1 Ry1t2e -i vt dt 2p L- q q (14.87) Substitution of Eq. (14.86) into Eq. (14.87) gives Sy1v2 = q 1 e -i vt dt 2p L- q q q * L- q L- q Introduction of Rx1t + u1 - u22h1u12h1u22du1 du2 ei vu1e -i vu2e -i v1u1 - u22 = 1 (14.88) (14.89) into Eq. (14.88) results in Sy1v2 = q L- q h1u12ei vu1 du1 L- q h1u22e -i vu2 du2 1 Rx1t + u1 - u22e -i v1u1 - u22 dt 2p L- q q * q (14.90) In the third integral on the right-hand side of Eq. (14.90), u1 and u2 are constants and the introduction of a new variable of integration h as h = t + u1 - u2 (14.91) RaoCh14ff.qxd 1004 10.06.08 14:00 Page 1004 CHAPTER 14 RANDOM VIBRATION leads to 1 Rx1t + u1 - u22e -iv1t + u1 - u22 dt 2p L- q q 1 Rx1h2e -ivh dh K Sx1v2 2p L- q q = (14.92) The first and the second integrals on the right-hand side of Eq. (14.90) can be recognized as the complex frequency response functions H1v2 and H1- v2, respectively. Since H1- v2 is the complex conjugate of H1v2, Eq. (14.90) gives Sy1v2 = ƒ H1v2 ƒ 2Sx1v2 (14.93) This equation gives the relationship between the power spectral densities of the excitation and the response. Mean Square Response. The mean square response of the stationary random process y(t) can be determined either from the autocorrelation function Ry1t2 or from the power spectral density Sy1v2: E[y2] = Ry102 = and q 2 E[y ] = L- q q q L- q L- q Rx1u1 - u22h1u12h1u22 du1 du2 Sy1v2 dv = q L- q ƒ H1v2 ƒ 2Sx1v2 dv (14.94) (14.95) Note: Equations (14.93) and (14.95) form the basis for the random vibration analysis of single- and multidegree of freedom systems [14.11, 14.12]. The random vibration analysis of road vehicles is given in Refs. [14.13, 14.14]. Mean Square Value of Response EXAMPLE 14.5 A single degree of freedom system (Fig. 14.20a) is subjected to a force whose spectral density is a white noise Sx1v2 = S0. Find the following: (a) Complex frequency response function of the system (b) Power spectral density of the response (c) Mean square value of the response Solution (a) To find the complex frequency response function H1v2, we substitute the input as eivt and the corresponding response as y1t2 = H1v2eivt in the equation of motion $ # my + cy + ky = x1t2 RaoCh14ff.qxd 10.06.08 14:00 Page 1005 14.13 RESPONSE DUE TO STATIONARY RANDOM EXCITATIONS 1005 FIGURE 14.20 Single degree of freedom system. to obtain and 1 - mv2 + icv + k2H1v2ei vt = ei vt H1v2 = 1 - mv2 + icv + k (E.1) (b) The power spectral density of the output can be found as Sy1v2 = ƒ H1v2 ƒ 2Sx1v2 = S0 ` 2 1 ` - mv2 + icv + k (E.2) RaoCh14ff.qxd 1006 10.06.08 14:00 Page 1006 CHAPTER 14 RANDOM VIBRATION (c) The mean square value of the output is given by4 q E[y2] = L- q Sy1v2 dv q = S0 L- q ` 2 pS0 1 ` dv = kc - mv + k + icv (E.3) 2 which can be seen to be independent of the magnitude of the mass m. The functions H1v2 and Sy1v2 are shown graphically in Fig. 14.20(b). ■ Design of the Columns of a Building EXAMPLE 14.6 A single-story building is modeled by four identical columns of Young’s modulus E and height h and a rigid floor of weight W, as shown in Fig. 14.21(a). The columns act as cantilevers fixed at the ground. The damping in the structure can be approximated by an equivalent viscous damping constant c. The ground acceleration due to an earthquake is approximated by a constant spectrum S0. If each column has a tubular cross section with mean diameter d and wall thickness t = d/10, find the mean diameter of the columns such that the standard deviation of the displacement of the floor relative to the ground does not exceed a specified value d. Solution Approach: Model the building as a single degree of freedom system. Use the relation between the power spectral densities of excitation and output. The building can be modeled as a single degree of freedom system as shown in Fig. 14.21(b) with m = W/g (E.1) FIGURE 14.21 Single-story building. 4 1B20 /A 02A 2 + B21 The values of this and other similar integrals have been found in the literature [14.1]. For example, if H1v2 = ivB1 + B0 q - v A 2 + ivA 1 + A 0 L- q 2 , ƒ H1v2 ƒ 2 dv = p e A 1A 2 f RaoCh14ff.qxd 10.06.08 14:00 Page 1007 14.13 RESPONSE DUE TO STATIONARY RANDOM EXCITATIONS 1007 and k = 4a 3EI b h3 (E.2) since the stiffness of one cantilever beam (column) is equal to 13EI/h32, where E is the Young’s modulus, h is the height, and I is the moment of inertia of the cross section of the columns given by (see Fig. 14.21(c)): I = p 4 1d - d4i 2 64 0 (E.3) Equation (E.3) can be simplified, using d0 = d + t and di = d - t, as I = = = p 2 1d + d2i 21d0 + di21d0 - di2 64 0 p [1d + t22 + 1d - t22][1d + t2 + 1d - t2][1d + t2 - 1d - t2] 64 p dt1d2 + t22 8 (E.4) With t = d/10, Eq. (E.4) becomes I = and hence Eq. (E.2) gives k = 101p 4 d = 0.03966d4 8000 12 E10.03966d42 3 h = 0.47592 Ed4 h3 (E.5) (E.6) When the base of the system moves, the equation of motion is given by (see Section 3.6) $ # $ mz + cz + kz = - mx (E.7) where z = y - x is the displacement of the mass (floor) relative to the ground. Equation (E.7) can be rewritten as k c # $ $ z + z = -x z + m m (E.8) The complex frequency response function H1v2 can be obtained by making the substitution $ x = eivt and z1t2 = H1v2eivt (E.9) RaoCh14ff.qxd 1008 10.06.08 14:00 Page 1008 CHAPTER 14 RANDOM VIBRATION so that c - v2 + iv which gives c k + dH1v2ei vt = - ei vt m m H1v2 = -1 k c a - v + iv + b m m (E.10) 2 The power spectral density of the response z(t) is given by Sz1v2 = ƒ H1v2 ƒ 2Sx$1v2 = S0 ∞ -1 a - v2 + iv k c + b m m ∞ 2 (E.11) The mean square value of the response z(t) can be determined, using Eq. (E.4) of Example 14.5, as q E[z2] = L- q Sz1v2dv q = S0 L- q = S0 a ∞ ∞ dv 2 -1 a - v2 + iv pm2 b kc c k + b m m (E.12) Substitution of the relations (E.1) and (E.6) into Eq. (E.12) gives E[z2] = pS0 W2h3 g c10.47592 Ed42 2 (E.13) Assuming the mean value of z(t) to be zero, the standard deviation of z can be found as sz = 4E[z2] = B 0.47592g2cEd4 pS0W 2h3 (E.14) Since sz … d, we find that pS0W2h3 0.47592g2cEd4 … d2 or d4 Ú pS0W2h3 0.47592g2cEd2 (E.15) RaoCh14ff.qxd 10.06.08 14:00 Page 1009 14.14 RESPONSE OF A MULTIDEGREE OF FREEDOM SYSTEM 1009 Thus the required mean diameter of the columns is given by d Ú e pS0W 2h3 2 0.47592g cEd 2 f 1/4 (E.16) ■ 14.14 Response of a Multidegree of Freedom System The equations of motion of a multidegree of freedom system with proportional damping can be expressed, using the normal mode approach, as (see Eq. 6.128) $ # qi1t2 + 2 zi vi qi1t2 + v2i qi1t2 = Qi1t2; i = 1, 2, Á , n (14.96) where n is the number of degrees of freedom, vi is the ith natural frequency, qi1t2 is the ith generalized coordinate, and Qi1t2 is the ith generalized force. The physical and generalized coordinates are related as ! ! x1t2 = [X]q1t2 or xi 1t2 = a X1j2 i qj1t2 n (14.97) j=1 where [X] is the modal matrix and X1j2 i is the ith component of jth modal vector. The physical and generalized forces are related as ! ! Q1t2 = [X]T F1t2 or Qi1t2 = a X1i2 j Fj1t2 (14.98) Fj1t2 = fj t1t2 (14.99) n where Fj1t2 is the force acting along the coordinate xj1t2. Let the applied forces be expressed as j=1 so that Eq. (14.98) becomes Qi1t2 = a a X1i2 j fj bt1t2 = Ni t1t2 n j=1 (14.100) RaoCh14ff.qxd 1010 10.06.08 14:00 Page 1010 CHAPTER 14 RANDOM VIBRATION where Ni = a X1i2 j fj n (14.101) j=1 By assuming a harmonic force variation t1t2 = ei vt (14.102) the solution of Eq. (14.96) can be expressed as qi1t2 = Ni v2i Hi1v2 t1t2 (14.103) where Hi1v2 denotes the frequency response function Hi1v2 = 1 v 2 v 1 - a b + i 2 zi vi vi (14.104) The mean square value of the physical displacement, xi1t2, can be obtained from Eqs. (14.97) and (14.103) as 1 x2i 1t2 dt t : q 2T L -T x2i 1t2 = lim T n n N 1s2 Nr s = a a X1r2 lim i Xi 2 2 T: q v v r=1 s=1 r s From Eq. (3.56), Hr1v2 can be expressed as Hr1v2 Hs1v2 t21t2 dt L- T T Hr1v2 = ƒ Hr1v2 ƒ e -ifr where the magnitude of Hr1v2, known as the magnification factor, is given by ƒ Hr1v2 ƒ = e c 1 - a v 2 -1/2 v 2 2 b d + a2zr b f vr vr (14.105) (14.106) (14.107) and the phase angle, fr, by fr = tan -1 µ 2zr v vr v 2 1 - a b vr ∂ (14.108) RaoCh14ff.qxd 10.06.08 14:00 Page 1011 14.14 RESPONSE OF A MULTIDEGREE OF FREEDOM SYSTEM 1011 By neglecting the phase angles, the integral on the right-hand side of Eq. (14.105) can be expressed as 1 Hr1v2 Hs1v2 t21t2 dt T : q 2T L -T T lim 1 ƒ Hr1v2 ƒ ƒ Hs1v2 ƒ t21t2 dt T : q 2T L -T (14.109) 1 St1v2 dv t21t2 dt = q T : 2T L -T L- q (14.110) T L lim For a stationary random process, the mean square value of t21t2 can be expressed in terms of its power spectral density function, St1v2, as t21t2 = lim q T Combining Eqs. (14.109) and (14.110) gives 1 Hr1v2 Hs1v2 t21t2 dt T : q 2T L -T T lim q L L- q ƒ Hr1v2 ƒ ƒ Hs1v2 ƒ St1v2 dv (14.111) Substitution of Eq. (14.111) into Eq. (14.105) yields the mean square value of xi1t2 as n n N 1s2 Nr s x2i 1t2 L a a X1r2 i Xi v2 v2 r=1 s=1 r s q L- q ƒ Hr1v2 ƒ ƒ Hs1v2 ƒ St1v2 dv (14.112) The magnification factors, ƒ Hr1v2 ƒ and ƒ Hs1v2 ƒ , are shown in Fig. 14.22. It can be seen that the product, ƒ Hr1v2 ƒ ƒ Hs1v2 ƒ for r Z s, is often negligible compared to ƒ Hr1v2 ƒ 2 and ƒ Hs1v2 ƒ 2. Hence Eq. (14.112) can be rewritten as 2 n N 2r ƒ Hr1v2 ƒ 2St1v2 dv b x 2i 1t2 L a aX1r2 i v4r L- q r=1 q (14.113) For lightly damped systems, the integral in Eq. (14.113) can be evaluated by approximating the graph of St1v2 to be flat with St1v2 = St1vr2 as q L- q ƒ Hr1v2 ƒ 2 St1v2 dv L St1vr2 q L- q ƒ Hr1v2 ƒ 2 dv = pvr St1vr2 2zr (14.114) RaoCh14ff.qxd 1012 10.06.08 14:00 Page 1012 CHAPTER 14 RANDOM VIBRATION FIGURE 14.22 Magnification factors. 2 n pvr St1vr2 2 Nr b a 2 x2i 1t2 = a 1X1r2 i 4 2zr vr r=1 Equations (14.113) and (14.114) yield (14.115) The following example illustrates the computational procedure. Response of a Building Frame Under an Earthquake EXAMPLE 14.7 The three-story building frame shown in Fig. 14.23 is subjected to an earthquake. The ground acceleration during the earthquake can be assumed to be a stationary random process with a power spectral density S1v2 = 0.05 (m2/s4)/(rad/s). Assuming a modal damping ratio of 0.02 in each mode, FIGURE 14.23 Three-story building frame. RaoCh14ff.qxd 10.06.08 14:00 Page 1013 14.14 RESPONSE OF A MULTIDEGREE OF FREEDOM SYSTEM 1013 determine the mean square values of the responses of the various floors of the building frame under the earthquake. Solution: The stiffness and mass matrices of the building frame can be found as 2 [k] = kC - 1 0 1 [m] = m C 0 0 -1 2 -1 0 1 0 0 -1 S 1 0 0S 1 (E.1) (E.2) From Examples 6.10 and 6.11, the eigenvalues and the eigenvectors (normalized with respect to the mass matrix [m]) can be computed using the values k = 106 N/m and m = 1000 kg, as k = 14.0734 rad/s Am (E.3) v2 = 1.2471 k = 39.4368 rad/s Am (E.4) v3 = 1.8025 k = 57.0001 rad/s Am (E.5) v1 = 0.44504 1.0000 0.01037 ! 0.3280 µ 1.8019 ∂ = µ 0.01869 ∂ Z112 = 1m 2.2470 0.02330 (E.6) ! 0.7370 Z122 = µ 1m 1.0000 0.02331 0.4450 ∂ = µ 0.01037 ∂ - 0.8020 - 0.01869 (E.7) 1.0000 0.01869 ! 0.5911 Z132 = µ - 1.2468 ∂ = µ - 0.02330 ∂ 1m 0.5544 0.01036 (E.8) ! ! Note that the notation Z1i2 is used to denote the ith mode shape instead of X1i2 since the relative displacements, zi1t2, will be used instead of the absolute displacements, xi1t2, in the solution. By denoting the ground motion as y(t), the relative displacements of the floors, zi1t2, can be expressed as zi1t2 = xi1t2 - y1t2, i = 1, 2, 3. The equations of motion can be expressed as ! $! #! ! [m]x + [c]z + [k]z = 0 (E.9) RaoCh14ff.qxd 1014 10.06.08 14:00 Page 1014 CHAPTER 14 RANDOM VIBRATION which can be rewritten as $! #! $! ! [m]z + [c]z + [k]z = - [m]y (E.10) $ y $! ! $ where y = µ y ∂ . By expressing the vector z in terms of normal modes, we have $ y ! ! z = [Z]q (E.11) where [Z] denotes the modal matrix. By substituting Eq. (E.11) into Eq. (E.10) and premultiplying the resulting equation by [Z]T, we can derive the uncoupled equations of motion. Assuming a damping ratio zi1zi = 0.022 in mode i, the uncoupled equations of motion are given by where $ qi + 2 zi vi + v2i qi = Qi; i = 1, 2, 3 (E.12) Qi = a Z1i2 j Fj1t2 (E.13) $ $ Fj1t2 = - mj y1t2 = - m y1t2 (E.14) n j=1 and with mj = m denoting the mass of the j th floor. By representing Fj1t2 as we note that and Fj1t2 = fj t1t2 (E.15) fj = - mj = - m (E.16) $ t1t2 = y1t2 (E.17) 3 N 2r p 2 z2i 1t2 = a 1Z1r2 b St1vr2 a i 2 v3r 2zr r=1 (E.18) The mean square values, z2i 1t2, can be determined from Eq. (14.115): ! ! ! Nothing that Z112, Z122, and Z132 are given by Eqs. (E.6), (E.7), and (E.8), respectively, and 112 N1 = a Z112 = - 1000 10.052362 = - 52.36 j fj = - m a Z j 3 3 j=1 j=1 3 3 j=1 j=1 122 = - 1000 10.052372 = - 52.37 N2 = a Z122 j fj = - m a Z j (E.19) (E.20) RaoCh14ff.qxd 10.06.08 14:00 Page 1015 14.15 EXAMPLES USING MATLAB 132 N3 = a Z132 = - 1000 10.052352 = - 52.35 j fj = - m a Z j 3 3 j=1 j=1 1015 (E.21) Eq. (E.18) yields the mean square values of the relative displacements of the various floors of the building frame as z211t2 = 0.00053132 m2 (E.22) z221t2 = 0.00139957 m2 z231t2 = 0.00216455 m2 (E.23) (E.24) ■ Probability of Relative Displacement Exceeding a Specified Value EXAMPLE 14.8 Find the probability of the magnitude of the relative displacement of the various floors exceeding 1, 2, 3, and 4 standard deviations of the corresponding relative displacement for the building frame of Example 14.7. Solution Approach: Assume the ground acceleration to be a normally distributed random process with zero mean and use standard normal tables. $ Since the ground acceleration, y1t2, is assumed to be normally distributed with zero mean value, the relative displacements of the various floors can also be assumed to be normally distributed with zero mean values. Thus the standard deviations of the relative displacements of the floors are given by szi = 3 z2i 1t2; i = 1, 2, 3 The probability of the absolute value of the relative displacement, zi1t2, exceeding a specified number of standard deviations can be found from standard normal tables as (see Section 14.8): 0.31732 for p 0.04550 for p P[ ƒ zi1t2 ƒ 7 pszi] = µ 0.00270 for p 0.00006 for p = = = = 1 2 3 4 ■ 14.15 Examples Using MATLAB Plotting of Autocorrelation Function EXAMPLE 14.9 Using MATLAB, plot the autocorrelation function corresponding to white noise with spectral density S0 for the following cases: (a) band limited white noise with v1 = 0 and v2 = 4 rad/s, 6 rad/s, 8 rad/s (b) band limited white noise with v1 = 2 rad/s and v2 = 4 rad/s, 6 rad/s, 8 rad/s (c) ideal white noise RaoCh14ff.qxd 1016 10.06.08 14:00 Page 1016 CHAPTER 14 RANDOM VIBRATION Solution For (a) and (b), the autocorrelation function, R1t2, can be expressed as (from Example 14.4) R1t2 S0 For (c) it can be expressed as t : 0, R102 = lim e 2 S0 a t:0 = 2 1sin v2 t - sin v1 t2 t (E.1) v1 sin v1 t v2 sin v2 t b - 2 S0 a b f = 2 S01v2 - v12 v2 t v1 t For an ideal white noise, we let v1 = 0 and v2 : q , which yields R = 2S0 d1t2 where d1t2 is the Dirac delta function. The MATLAB program to plot Eq. (E.1) is given below. % Ex14_9.m w1 = 0; w2 = 4; for i = 1:101 t(i) = ⫺5 + 10* (i⫺1)/100; R1(i) = 2 * ( sin(w2 *t(i)) ⫺ sin(w1*t(i)) )/t(i); end w1 = 0; w2 = 6; for i = 1:101 t(i) = ⫺5 + 10*(i⫺1)/100; R2(i) = 2 * ( sin(w2 *t(i)) ⫺ sin(w1*t(i)) )/t(i); end w1 = 0; w2 = 8; for i = 1:101 t(i) = ⫺ 5 + 10*(i⫺1)/100; R3(i) = 2 * ( sin(w2 *t(i)) ⫺ sin(w1*t(i)) )/t(i); end for i = 1:101 t1(i) = 0.0001 + 4.9999*(i⫺1)/100; R3_1(i) = 2 * ( sin(w2 *t(i)) ⫺ sin(w1*t(i)) )/t(i); end xlabel ('t'); ylabel('R/S_0'); plot(t,R1); hold on; gtext('Solid line: w1 = 0, w2 = 4') gtext('Dashed line: w1 = 0, w2 = 6'); plot(t,R2,'--'); gtext('Dotted line: w1 = 0, w2 = 8'); plot(t,R3,':'); w1 = 2; w2 = 4; for i = 1:101 t(i) = ⫺5 + 10*(i⫺1)/100; R4(i) = 2 * ( sin(w2 *t(i)) ⫺ sin(w1*t(i)) )/t(i); end w1 = 2; w2 = 6; for i = 1:101 t(i) = ⫺5 + 10*(i⫺1)/100; R5(i) = 2 * ( sin(w2 *t(i)) ⫺ sin(w1*t(i)) )/t(i); end w1 = 2; w2 = 8; for i = 1 : 101 t (i) = ⫺5 + 10* (i⫺1) / 100; R6 (i) = 2 * ( sin (w2 *t (i) ) ⫺ sin (w1 *t (i) ) ) / t (i); end pause RaoCh14ff.qxd 10.06.08 14:00 Page 1017 14.15 EXAMPLES USING MATLAB 1017 hold off; xlabel ('t'); ylabel ('R/S_0'); plot (t, R4); hold on; gtext ('Solid line: w1 = 2, w2 = 4') gtext ('Dashed line: w1 = 2, w2 = 6'); plot (t, R5, '--'); gtext ('Dotted line: w1 = 2, w2 = 8'); plot (t, R6, ':'); ■ RaoCh14ff.qxd 1018 10.06.08 14:00 Page 1018 CHAPTER 14 RANDOM VIBRATION Evaluation of a Gaussian Probability Distribution Function EXAMPLE 14.10 Using MATLAB, evaluate the following probability for c = 1, 2 and 3: Prob 3 ƒ x1t2 ƒ Ú c s4 = q 1x 2 e - 2 s dx 12ps Lc s 2 2 (E.1) Assume the mean value of x(t) to be zero and standard deviation of x(t) to be one. Solution Equation (E.1) can be rewritten, for s = 1, as 1 2 Prob [ ƒ x1t2 ƒ Ú c] = 2e e -0.5x dx f 12p 3 c (E.2) -q The MATLAB program to evaluate Eq. (E.2) is given below. Ex14_10.m >> q = quad ('normp', ⫺7, 1); >> prob1 = 2 * q prob1 = 1.6827 >> q = quad ('normp', ⫺7, 2); >> prob2 = 2 * q prob2 = 1.9545 >> q = quad ('normp', ⫺7, 3); >> prob3 = 2 * q prob3 = 1.9973 %normp.m function pdf = normp(x) pdf = exp(⫺0.5*x.^2)/sqrt(2.0 * pi); ■ REFERENCES 14.1 S. H. Crandall and W. D. Mark, Random Vibration in Mechanical Systems, Academic Press, New York, 1963. 14.2 D. E. Newland, An Introduction to Random Vibrations and Spectral Analysis, Longman, London, 1975. RaoCh14ff.qxd 10.06.08 14:00 Page 1019 REVIEW QUESTIONS 1019 14.3 J. D. Robson, An Introduction to Random Vibration, Edinburgh University Press, Edinburgh, 1963. 14.4 C. Y. Yang, Random Vibration of Structures, Wiley, New York, 1986. 14.5 A. Papoulis, Probability, Random Variables and Stochastic Processes, McGraw-Hill, New York, 1965. 14.6 J. S. Bendat and A. G. Piersol, Engineering Applications of Correlation and Spectral Analysis, Wiley, New York, 1980. 14.7 P. Z. Peebles, Jr., Probability, Random Variables, and Random Signal Principles, McGraw-Hill, New York, 1980. 14.8 J. B. Roberts, “The response of a simple oscillator to band-limited white noise,” Journal of Sound and Vibration, Vol. 3, 1966, pp. 115–126. 14.9 M. H. Richardson, “Fundamentals of the discrete Fourier transform,” Sound and Vibration, Vol. 12, March 1978, pp. 40–46. 14.10 E. O. Brigham, The Fast Fourier Transform, Prentice Hall, Englewood Cliffs, N.J., 1974. 14.11 J. K. Hammond, “On the response of single and multidegree of freedom systems to non-stationary random excitations,” Journal of Sound and Vibration, Vol. 7, 1968, pp. 393–416. 14.12 S. H. Crandall, G. R. Khabbaz, and J. E. Manning, “Random vibration of an oscillator with nonlinear damping,” Journal of the Acoustical Society of America, Vol. 36, 1964, pp. 1330–1334. 14.13 S. Kaufman, W. Lapinski, and R. C. McCaa, “Response of a single-degree-of-freedom isolator to a random disturbance,” Journal of the Acoustical Society of America, Vol. 33, 1961, pp. 1108–1112. 14.14 C. J. Chisholm, “Random vibration techniques applied to motor vehicle structures,” Journal of Sound and Vibration, Vol. 4, 1966, pp. 129–135. 14.15 Y. K. Lin, Probabilistic Theory of Structural Dynamics, McGraw-Hill, New York, 1967. 14.16 I. Elishakoff, Probabilistic Methods in the Theory of Structures, Wiley, New York, 1983. 14.17 H. W. Liepmann, “On the application of statistical concepts to the buffeting problem,” Journal of the Aeronautical Sciences, Vol. 19, No. 12, 1952, pp. 793–800, 822. 14.18 W. C. Hurty and M. F. Rubinstein, Dynamics of Structures, Prentice Hall, Englewood Cliffs, N.J., 1964. REVIEW QUESTIONS 14.1 Give brief answers to the following: 1. 2. 3. 4. 5. 6. What is the difference between a sample space and an ensemble? Define probability density function and probability distribution function. How are the mean value and variance of a random variable defined? What is a bivariate distribution function? What is the covariance between two random variables X and Y? Define the correlation coefficient, rXY. RaoCh14ff.qxd 1020 10.06.08 14:00 Page 1020 CHAPTER 14 RANDOM VIBRATION 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. What are the bounds on the correlation coefficient? What is a marginal density function? What is autocorrelation function? Explain the difference between a stationary process and a nonstationary process? What are the bounds on the autocorrelation function of a stationary random process? Define an ergodic process. What are temporal averages? What is a Gaussian random process? Why is it frequently used in vibration analysis? What is Parseval’s formula? Define the following terms: power spectral density function, white noise, band-limited white noise, wide-band process, and narrow-band process. How are the mean square value, autocorrelation function, and the power spectral density function of a stationary random process related? What is an impulse response function? Express the response of a single degree of freedom system using the Duhamel integral. What is complex frequency response function? How are the power spectral density functions of input and output of a single degree of freedom system related? What are Wiener-Khintchine formulas? 14.2 Indicate whether each of the following statements is true or false: 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 14.3 A deterministic system requires deterministic system properties and loading. Most phenomena in real life are deterministic. A random variable is a quantity whose magnitude cannot be predicated precisely. The expected value of x, in terms of its probability density function, p(x), is given by q 1- qxp1x2 dx. The correlation coefficient rXY satisfies the relation ƒ rXY ƒ … 1. The autocorrelation function R1t1, t22 is the same as E[x1t12x1t22]. The mean square value of x(t) can be determined as E[x2] = R102. If x(t) is stationary, its mean will be independent of t. The autocorrelation function R1t2 is an even function of t. The Wiener-Khintchine formulas relate the power spectral density to the autocorrelation function. The ideal white noise is a physically realizable concept. Fill in each of the following blanks with the appropriate word: 1. When the vibrational response of a system is known precisely, the vibration is called _____ vibration. 2. If any parameter of a vibrating system is not known precisely, the resulting vibration is called _____ vibration. 3. The pressure fluctuation at a point on the surface of an aircraft flying in the air is a _____ process. 4. In a random process, the outcome of an experiment will be a function of some _____ such as time. 5. The standard deviation is the positive square root of _____. RaoCh14ff.qxd 10.06.08 14:00 Page 1021 REVIEW QUESTIONS 1021 6. The joint behavior of several random variables is described by the _____ probability distribution function. 7. Univariate distributions describe the probability distributions of _____ random variables. 8. The distribution of two random variables is known as _____ distribution. 9. The distribution of several random variables is called _____ distribution. 10. The standard deviation of a stationary random process x(t) will be independent of _____. 11. If all the probability information of a stationary random process can be obtained from a single sample function, the process is said to be _____. 12. The Gaussian density function is a symmetric _____-shaped curve about the mean value. 13. The standard normal variable has mean of _____ and standard deviation of _____. 14. A nonperiodic function can be treated as a periodic function having an _____ period. 15. The _____ spectral density function is an even function of v. 16. If S1v2 has significant values over a wide range of frequencies, the process is called a _____ process. 17. If S1v2 has significant values only over a small range of frequencies, the process is called a _____ process. 18. The power spectral density S1v2 of a stationary random process is defined as the _____ transform of R1t2/2p. 19. If the band of frequencies has finite cut-off frequencies for a white noise, it is called _____ white noise. 14.4 Select the most appropriate answer out of the choices given: 1. Each outcome of an experiment for a random variable is called (a) a sample point (b) a random point (c) an observed value 2. Each outcome of an experiment, in the case of a random process, is called a (a) sample point (b) sample space (c) sample function ' 3. The probability distribution function, P1x2, denotes ' (a) P1x … x2 ' (b) P1x 7 x2 ' ' (c) P1x … x … x + ¢x2 ' 4. The probability density function, p1x2, denotes ' (a) P1x … x2. ' (b) P1x 7 x2 ' ' (c) P1x … x … x + ¢x2 5. Normalization of probability distribution implies (a) P1 q 2 = 1 q L- q 6. The variance of x is given by (b) q p1x2 = 1 (c) L- q p1x2 = 0 (a) x2 (b) 1x22 - 1x22 (c) 1x22 7. The marginal density function of x can be determined form the bivariate density function p(x,y) as q (a) p1x2 = L- q p1x, y2 dy RaoCh14ff.qxd 1022 10.06.08 14:00 Page 1022 CHAPTER 14 RANDOM VIBRATION q (b) p1x2 = L- q q p1x, y2 dx q p1x, y2 dx dy L- q L- q 8. The correlation coefficient of x and y is given by (a) sxy (b) sxy/1sxsy2 (c) sxsy 9. The standard normal variable, z, corresponding to the normal variable x, is defined as x x x - x (a) z = (b) z = (c) z = sx sx sx 10. If the excitation of a linear system is a Gaussian process, the response will be (a) a different random process (b) a Gaussian process (c) an ergodic process 11. For a normal probability density function, Prob[- 3s … x1t2 … 3s] is (a) 0.6827 (b) 0.999937 (c) 0.9973 12. The mean square response of a stationary random process can be determined from the: (a) autocorrelation function only (b) power spectral density only (c) autocorrelation function or power spectral density 14.5 Match the items in the two columns below: (c) p1x2 = 1. All possible outcomes of a random variable (a) Correlation functions in an experiment 2. All possible outcomes of a random process (b) Nonstationary process 3. Statistical connections between the values of x(t) at times t1, t2, Á (c) Sample space 4. A random process invariant under a shift of the time scale (d) White noise 5. Mean and standard deviations of x(t) vary with t (e) Stationary process 6. Power spectral density is constant over a frequency range (f) Ensemble PROBLEMS The problem assignments are organized as follows: Problems 14.1, 14.10 14.3, 14.11 Section Covered 14.3 14.4 Topic Covered Probability distribution Mean value and standard deviation RaoCh14ff.qxd 10.06.08 14:00 Page 1023 PROBLEMS Problems Section Covered 14.2, 14.4 14.7, 14.9 14.5, 14.26 14.12 14.8, 14.13–14.16, 14.27 14.6, 14.17–14.22 14.23–14.25, 14.28–14.30 14.31, 14.32 14.33–14.35 14.36 14.5 14.6 14.7 14.8 14.9 14.10 14.12 14.14 14.15 — 1023 Topic Covered Joint probability distribution Correlation functions Stationary random process Gaussian random process Fourier analysis Power spectral density Response of a single degree of freedom system Response of a multidegree of freedom system MATLAB programs Design project 14.1 The strength of the foundation of a reciprocating machine (x) has been found to vary between 20 and 30 kips/ft 2 according to the probability density function: ka 1 p1x2 = L 0, x b, 30 20 … x … 30 elsewhere What is the probability of the foundation carrying a load greater than 28 kips/ft 2? 14.2 The joint density function of two random variables X and Y is given by xy , 9 pX,Y1x, y2 = L 0, 0 … x … 2, 0 … y … 3 elsewhere (a) Find the marginal density functions of X and Y. (b) Find the means and standard deviations of X and Y. (c) Find the correlation coefficient rX,Y. 14.3 The probability density function of a random variable x is given by 0 p1x2 = c 0.5 0 for x 6 0 for 0 … x … 2 for x 7 2 Determine E[x], E[x2], and sx. 14.4 If x and y are statistically independent, then E[xy] = E[x]E[y]. That is, the expected value of the product xy is equal to the product of the separate mean values. If z = x + y, where x and y are statistically independent, show that E[z2] = E[x2] + E[y2] + 2E[x]E[y]. 14.5 The autocorrelation function of a random process x(t) is given by Rx1t2 = 20 + Find the mean square value of x(t). 5 1 + 3t2 RaoCh14ff.qxd 1024 10.06.08 14:00 Page 1024 CHAPTER 14 RANDOM VIBRATION 14.6 The autocorrelation function of a random process is given by Rx1t2 = A cos vt; - p t … t … 2v 2v where A and v are constants. Find the power spectral density of the random process. 14.7 Find the autocorrelation functions of the periodic functions shown in Fig. 14.24. FIGURE 14.24 14.8 Find the complex form of the Fourier series for the wave shown in Fig. 14.24(b). 14.9 Compute the autocorrelation function of a periodic square wave with zero mean value and compare this result with that of a sinusoidal wave of the same period. Assume the amplitudes to be the same for both waves. 14.10 The life T in hours of a vibration transducer is found to follow exponential distribution pT1t2 = e le -lt, 0, t Ú 0 t 6 0 where l is a constant. Find (a) the probability distribution function of T, (b) mean value of T, and (c) standard deviation of T. 14.11 Find the temporal mean value and the mean square value of the function x1t2 = x0 sin1pt/22. 14.12 An air compressor of mass 100 kg is mounted on an undamped isolator and operates at an angular speed! of 1800 rpm. The stiffness of the isolator is found to be a random variable with mean value k = 2.25 * 106 N/m and standard deviation sk = 0.225 * 106 N/m following normal distribution. Find the probability of the natural frequency of the system exceeding the forcing frequency. 14.13– 14.16 Find the Fourier transform of the functions shown in Figs. 14.25–14.28 and plot the corresponding spectrum. 14.17 A periodic function F(t) is shown in Fig. 14.29. Use the values of the function F(t) at ten equally spaced time stations ti to find (a) the spectrum of F(t) and (b) the mean square value of F(t). Rx1t2 = ae -bƒt ƒ 14.18 The autocorrelation function of a stationary random process x(t) is given by where a and b are constants. Find the power spectral density of x(t). RaoCh14ff.qxd 10.06.08 14:00 Page 1025 PROBLEMS FIGURE 14.25 FIGURE 14.26 FIGURE 14.27 FIGURE 14.28 1025 RaoCh14ff.qxd 1026 10.06.08 14:00 Page 1026 CHAPTER 14 RANDOM VIBRATION FIGURE 14.29 14.19 Find the autocorrelation function of a random process whose power spectral density is given by S1v2 = S0 = constant between the frequencies v1 and v2. 14.20 The autocorrelation function of a Gaussian random process representing the unevenness of a road surface is given by Rx1t2 = s2xe -aƒ vtƒ cos bvt where s2x is the variance of the random process and v is the velocity of the vehicle. The values of sx, a, and b for different types of road are as follows: Type of road Sx A B Asphalt Paved Gravel 1.1 1.6 1.8 0.2 0.3 0.5 0.4 0.6 0.9 Compute the spectral density of the road surface for the different types of road. 14.21 Compute the autocorrelation function corresponding to the ideal white noise spectral density. 14.22 Starting from Eqs. (14.60) and (14.61), derive the relations R1t2 = L0 S1f2 = 4 q L0 S1f2 cos 2pft # df q R1t2 cos 2pft # dt RaoCh14ff.qxd 10.06.08 14:00 Page 1027 PROBLEMS 1027 14.23 Write a computer program to find the mean square value of the response of a single degree of freedom system subjected to a random excitation whose power spectral density function is given as Sx1v2. 14.24 A machine, modeled as a single degree of freedom system, has the following parameters: mg = 2000 lb, k = 4 * 104 lb/in., and c = 1200 lb-in./sec. It is subjected to the force shown in Fig. 14.29. Find the mean square value of the response of the machine (mass). 14.25 A mass, connected to a damper as shown in Fig. 14.30, is subjected to a force F(t). Find the frequency response function H1v2 for the velocity of the mass. FIGURE 14.30 14.26 The spectral density of a random signal is given by S1f2 = e 0.0001 m2/cycle/sec, 0, 10 Hz … f … 1000 Hz elsewhere Find the standard deviation and the root mean square value of the signal by assuming its mean value to be 0.05 m. 14.27 Derive Eq. (14.46) from Eq. (14.45). 14.28 A simplified model of a motor cycle traveling over a rough road is shown in Fig. 14.31. It is assumed that the wheel is rigid, the wheel does not leave the road surface, and the cycle moves at a constant speed v. The cycle has a mass m and the suspension system has a spring constant k and a damping constant c. If the power spectral density of the rough road surface is taken as S0, find the mean square value of the vertical displacement of the motor cycle (mass, m). FIGURE 14.31 RaoCh14ff.qxd 1028 10.06.08 14:00 Page 1028 CHAPTER 14 RANDOM VIBRATION 14.29 The motion of a lifting surface about the steady flight path due to atmospheric turbulence can be represented by the equation 1 $ # x1t2 + 2zvnx1t2 + v2nx1t2 = F1t2 m where vn is the natural frequency, m is the mass, and z is the damping coefficient of the system. The forcing function F(t) denotes the random lift due to the air turbulence and its spectral density is given by [14.17] SF1v2 = ST1v2 a1 + pvc b v where c is the chord length and v is the forward velocity of the lifting surface and ST1v2 is the spectral density of the upward velocity of air due to turbulence, given by ST1v2 = A2 1 + a e1 + a Lv 2 b v Lv 2 2 b f v where A is a constant and L is the scale of turbulence (constant). Find the mean square value of the response x(t) of the lifting surface. 14.30 The wing of an airplane flying in gusty wind has been modeled as a spring-mass-damper system, as shown in Fig. 14.32. The undamped and damped natural frequencies of the wing are found to be v1 and v2, respectively. The mean square value of the displacement of meq (i.e., the wing) is observed to be d under the action of the random wind force whose power spectral density is given by S1v2 = S0. Derive expressions for the system parameters meq, keq, and ceq in terms of v1, v2, d, and S0. FIGURE 14.32 14.31 If the building frame shown in Fig. 14.23 has a structural damping coefficient of 0.01 (instead of the modal damping ratio 0.02), determine the mean square values of the relative displacements of the various floors. RaoCh14ff.qxd 10.06.08 14:00 Page 1029 DESIGN PROJECT 1029 14.32 The building frame shown in Fig. 14.23 is subjected to a ground acceleration whose power spectral density is given by S1v2 = 1 4 + v2 Find the mean square values of the relative displacements of the various floors of the building frame. Assume a modal damping ratio of 0.02 in each mode. 14.33 Using MATLAB, plot the Gaussian probability density function f1x2 = 1 2 e -0.5x 12p over - 7 … x … 7. 14.34 Plot the Fourier transform of a triangular pulse: X1v2 = 4A va , sin2 2 a v2 -7 … va … 7 p (See Fig. 14.12.) 14.35 The mean square value of the response of a machine, E[y2], subject to the force shown in Fig. 14.29, is given by (see Problem 14.24): E[y2] = a N-1 n=0 ƒ cn ƒ 2 a k - mv2n b + c2v2n 2 where cn = 2pnj 2pnj 1 N - i sin f Fj e cos a N j=1 N N with Fj = 0, 20, 40, 60, 80, 100, - 80, -60, - 40, -20, 0 for j = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10; k = 4 * 104, c = 1200, m = 5.1760 and vn = 2 p n. Using MATLAB, find the value of E[y2] with N = 10. DESIGN PROJECT 14.36 The water tank shown in Fig. 14.33 is supported by a hollow circular steel column. The tank, made of steel, is in the form of a thin-walled pressure vessel and has a capacity of 10,000 gallons. Design the column to satisfy the following specifications: (a) The undamped natural frequency of vibration of the tank, either empty or full, must exceed a RaoCh14ff.qxd 1030 10.06.08 14:00 Page 1030 CHAPTER 14 RANDOM VIBRATION value of 1 Hz. (b) The mean square value of the displacement of the tank, either empty or full, must not exceed a value of 16 in2 when subjected to an earthquake ground acceleration whose power spectral density is given by S1v2 = 0.0002 m2/s4 rad/s Assume damping to be 10 percent of the critical value. FIGURE 14.33