A unified approach for treating linear multibody systems involving flexible beams

Mechanism and Machine Theory 107 (2017) 197–209

Contents lists available at ScienceDirect

Mechanism and Machine Theory journal homepage: www.elsevier.com/locate/mechmachtheory

A unified approach for treating linear multibody systems involving flexible beams

cross

⁎

Laith K. Abbasa, , Qinbo Zhoua, Dieter Bestleb, Xiaoting Ruia a b

Institute of Launch Dynamics, Nanjing University of Science and Technology, Nanjing 210094, China Engineering Mechanics and Vehicle Dynamics, BTU-Cottbus-Senftenberg, P.O. Box 10 13 44, 03013 Cottbus, Germany

A R T I C L E I N F O

ABSTRACT

Keywords: Linear multibody systems Flexible beams Transfer matrix method for multibody systems Linear vibration

Over the last decades, various strategies have been proposed to exactly treat linear multibody systems composed of multibody system elements and flexible beams. Mostly they involve advanced calculus and are devoted to a specific system topology. On the other hand, general multibody system approaches have been developed where, however, elastic bodies can be treated only approximately. The transfer matrix method for multibody systems (MSTMM) combines both advantages: it may account for any general multibody system structure and treat beams exactly. This is achieved by introducing new connection elements in the transfer matrix method which allow to set up the overall algebraic transfer equations easily and to determine vibration frequencies by solving minimization problems. The proposed approach will be demonstrated for several beam/body examples taken from specialized literature.

1. Introduction A multibody system may be composed of several bodies (rigid or flexible), ideal joints, and force elements in an arbitrary fashion. Typically, dynamics of such multibody systems is deduced from classical mechanics and used for many engineering applications, such as manipulators, gyroscopes, satellites, robots, launch systems, biomechanical systems, mechanisms, etc. [1–4]. For small vibrations, the equations of motion may be linearized and especially the natural vibration frequencies and modes are the most important dynamics characteristics of these kind of systems. This information is typically required for assessing dynamic system performance or designing system control. Classical multibody system dynamics considers rigid bodies only. However, modern engineering technology also requires consideration of flexible bodies originally described by partial differential equations. In order to fit them to the concept of multibody dynamics, typically approximation methods associated with the Finite Element Method are utilized [5,6]. Exact solutions are reported only for simple beam problems where very specific solution procedures have been developed within full papers [7–16]. The transfer matrix method for multibody systems (MSTMM) may resolve this problem of finding exact solutions for beam problems of any topology within a few lines of analytical calculus. Transfer matrix methods have been developed for a long time and used widely in engineering applications [17–19]. The general idea of MSTMM is to firstly break up a complicated multibody system into elements with simple dynamic properties, which can be expressed by algebraic transfer equations (TE) instead of differential equations. These element TEs relating output states to input states may then be easily combined to overall TEs relating boundary states of the multibody system. The element transfer matrices need to be determined only once and may then be collected in a transfer matrix library [20] providing matrices for lumped mass,

⁎

Corresponding author. E-mail addresses: [email protected] (L.K. Abbas), [email protected] (Q. Zhou), [email protected] (D. Bestle), [email protected] (X. Rui).

http://dx.doi.org/10.1016/j.mechmachtheory.2016.09.022 Received 7 July 2016; Received in revised form 26 September 2016; Accepted 26 September 2016 0094-114X/ © 2016 Elsevier Ltd. All rights reserved.

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

rigid body, beam, spring, damper and joints. These matrices are assembled according to the topology of the flexible multibody system. In case of chain topology this is simply matrix multiplication, else additional compliance equations need to be added. A problem is that element matrices are not only related to the type of element, but also to the model dimension. Thus, if a 1Delement like a translational spring is treated in a 2D-model, artificial stiffness has to be added in the missing dimension. In this paper, specific connection elements will be introduced in Section 3 to avoid such extensions and to ease the build-up of the overall transfer equation of the system. Various examples with increasing complexity will be taken from the literature to demonstrate the concept, to introduce the new connection elements step-by-step, and to show the generality of the approach. For sake of simplicity only planar systems will be investigated, although spatial systems may be handled the same way. For completeness, basic notations and transfer matrices for some elements will be firstly given in Section 2 before they are applied to the examples from literature in Section 4. Readers already familiar with the transfer matrix method may skip Section 2 and proceed with the new connection elements in Section 3, or others may proceed with Section 4 and come back to the element and connection transfer equations later when they are used for the first time.

2. MSTMM modeling of elements within the scope of finding vibration characteristics Bodies and coupling elements (both called just elements later on) of a multibody system are connected at nodal points transferring forces/moments and displacements/angles. In the view of MSTMM, one of the boundary points is called root, the others are called tips, and the directions from tips to the root are called transfer directions. Along these transfer directions, the nodes entering into elements are called inputs denoted by I , and the nodes leaving from elements are called outputs O . Elements may also have several inputs and outputs counted by In and Ol , respectively. Typically, there are several elements along a transfer path. A connection point between two elements is both output of the former element and input of the latter element. The state of a connection point is given by kinematic and kinetic quantities in physical coordinates: in case of linear multibody systems, vibrations are described by small displacements x, y, z along the Cartesian axes and small angular rotations θx , θy, θz about these axes; cutting forces and moments are given by qx , qy, qz and mx , m y, mz , respectively. Positive directions at input points are shown in Fig. 1a, whereas positive directions of forces and moments at output points are opposite due to the principle that action equals reaction, Fig. 1b. In 3D-case, the number of state variables equals ns = 12 and the physical states at a connection point may be summarized in a state vector (SV)

z = [x, y , z, θx , θy, θz, mx , m y, mz , qx , qy, qz ]T

(1)

For special cases with reduced dimension used later in the paper, this state vector may be condensed to the essential coordinates of the specific element type as e.g. ns = 6, z = [x, y, θz, mz , qx , qy]T in 2D or as ns = 4, z = [ y, θz, mz , qy]T , ns = 2, z = [ y, qy]T and ns = 2, z = [θz, mz ]T for 1D elements. Details will become clear later on. Due to the sign convention in Fig. 1, connecting an output point of an element j with an input point of another element k then may be written simply as zIk ≡ zOj . For modal vibrations, solutions may be written as z = Z e λt where Z = [X , Y , Z , Θx , Θy , Θz , Mx , My, Mz , Qx , Qy , Qz ]T denotes the vibration amplitudes in modal coordinates and λ = λr ± iλi , are the eigenvalues. The real part λr is typically negative and related to the magnitude of damping, whereas the imaginary part λi is related to the vibration frequency of the damped system. For undamped systems, λr = 0 and λi = ω . In the context of MSTMM, for a general element numbered as j with single input and single output point, the TEs can be obtained as ZO = Uj ZI where Uj is the transfer matrix of the jth-element expressing the relationship between the element SVs related to its output end ZO and its input end ZI [20]. Some typical elements of the multibody systems, which are considered later on, are shown in Fig. 2. Each element is characterized by its specific transfer matrix, which may be easily derived and saved in a MSTMM library. The transfer matrices of the elements can then be regarded as building blocks and the mechanical behavior of the global system is obtained by assembling these blocks according to the topology by predetermined rules [21]. In the following, transfer matrices of planar elements for beam/body systems vibrating in a plane are introduced briefly which is sufficient for the intended comparison with results from the literature.

Fig. 1. Sign convention at (a) input end and (b) output end.

198

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

Fig. 2. Some typical elements in TMM library: (a) linear spring, (b) torsional spring, (c) lumped mass, (d) viscous damper, (e) rigid body, and (f) beam.

2.1. Linear translational spring For the one-dimensional vibration of a massless linear translational spring shown in Fig. 2a, the SV (1) may be reduced to the essential coordinates zI = [ yI , qy I ]T and zO = [ yO , qy O]T , respectively. From the equilibrium condition and the force law of linear springs with stiffness Ky we find

⎧ qy O = qy I ⎧ yO = yI − qy I / Ky ⎨ or ⎨ q = K ( y − y ) y I O ⎩ yI ⎩ qy O = qy I By applying Laplace transformation y = Y

e λt ,

(2)

qy = Qy

e λt

we obtain the TE

⎡1 −1/ Ky ⎤ ⎡1 −1/ Ky ⎤ ⎡Y ⎤ ⎡Y ⎤ ⎥ ⎥ ⎢ ⎥ or ZO = ULS ZI where ULS = ⎢ ⎢Q ⎥ = ⎢ ⎣ y ⎦O ⎣0 ⎣0 1 ⎦ ⎣Qy ⎦ I 1 ⎦

(3)

is the transfer matrix of a linear translational spring element. Transfer matrices of linear translational springs in other directions are identical, only state variables have to be adapted. 2.2. Torsional spring A torsional spring with stiffness K′z may vibrate about the z − axis as shown in Fig. 2b. The moments at the input and output ends are equal and given by the rotation angles as mz O = mz I = K ′z (θz O − θz I ). Applying a similar transformation θz = Θz e λt , mz = Mz e λt as above results in the TE

⎡ Θz ⎤ ⎡ 1 1/ K ′z ⎤ ⎡ Θz ⎤ ⎢ ⎥ =⎢ ⎥ ⎢ ⎥ or ⎣ Mz ⎦O ⎣ 0 1 ⎦ ⎣ Mz ⎦ I

⎡ 1 1/ K ′z ⎤ ZO = UTS ZI where UTS = ⎢ ⎥ ⎣0 1 ⎦

(4)

2.3. Lumped mass Fig. 2c shows a lumped mass m (point mass without moment effects) vibrating in y − direction. The displacements at the lower and upper ends are equal, i.e., yO = yI . Referring to Newton's second law m y¨I = qy I − qy O , and using the differentiation rule y¨I = YI λ2e λt results in the TE

⎡Y ⎤ ⎡ 1 0⎤ ⎡ Y ⎤ ⎢Q ⎥ = ⎢ ⎥ ⎢ ⎥ or ⎣ y ⎦O ⎣− mλ2 1 ⎦ ⎣Qy ⎦ I

⎡ 1 0⎤ ZO = UM ZI where UM = ⎢ ⎥ ⎣− mλ2 1 ⎦ 199

(5)

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

2.4. Viscous damper The massless dashpot in Fig. 2d with damping coefficient C has equal forces at its two ends given as qy O = qy I = C ( yİ − yȮ ). With y ̇ = Yλ e λt we get the TE

⎡Y ⎤ ⎡ 1 − 1/ Cλ ⎤ ⎡ Y ⎤ ⎢Q ⎥ = ⎢ ⎥ ⎢ ⎥ or 1 ⎦ ⎣Qy ⎦ I ⎣ y ⎦O ⎣ 0

⎡ ⎤ ZO = UD ZI where UD = ⎢ 1 − 1/ Cλ ⎥ 1 ⎦ ⎣0

(6)

2.5. Dynamics of a planar rigid body with single input and single output ends The rigid body in Fig. 2e with mass m and moment of inertia JZ I about its input end I may vibrate in the x, y -plane. In a bodyfixed frame {I , xI yI zI }, the coordinates of mass center C are lIC = [c1, c2 , 0]T and the coordinates of output end O are lIO = [b1, b 2 , 0]T . Due to planar motion, i.e. z = θx = θy = qz = mx = m y = 0 , the SV (1) reduces to Z = [X , Y , Θz , Mz , Qx , Qy]T . For rigid body motion, all rotation angles are equal, i.e.,

θ: =θI ≡ θO ⇒ θ¨ O = θ¨ I

(7)

Further, rigid body kinematics (8)

rO = rI + rIO, vO = vI + ω × rIO, aO = aI + ω̇ × rIO + ω × (ω × rIO) may be formulated in the inertial frame with

⎡ cos θ − sin θ 0 ⎤ rO = [ xO yO 0]T , rI = [ xI yI 0]T , rIO = A lIO, A= ⎢ sin θ cos θ 0 ⎥ , ω= [ 0 0 θİ ]T ⎢⎣ ⎥ 0 0 1⎦

(9)

Linearization of acceleration relation (8) with respect to θ < < 1 then yields for the x, y –coordinates

x¨O = x¨I − θ¨I b 2 ,

y¨O = y¨I + θ¨I b1

(10)

Newton's law

m aC ≡ m aI + m ω̇ × rIC + m ω × (ω × rIC) =

∑F

(11)

may be handled in the same way resulting in

mx¨I − mc2 θ¨I = qx I − qx O, my¨I + m c1 θ¨I = qy I − qy O

(12)

Finally, the z − coordinate of Euler's law

II ω̇ + ω × II ω + mrIC × aI =

∑ MI

(13)

After linearization yields the relation

JzI θ¨I + m (c1 y¨I − c2 x¨I ) = −mzI + mz O − b1 qyO + b 2 qxO

(14)

where qxO, qyO may be substituted by Eq. (12). Modal transformation of Eqs. (7), (10), (12) and (14) and combination in a TE finally results in ZO = URB ZI where

URB

⎡ 1 0 − b2 ⎢ 0 1 b1 ⎢ 0 0 1 ⎢ = ⎢ mλ2 (b − c ) − mλ2 (b − c ) λ2 [J − m (b c + b c )] 2 2 1 1 zI 2 2 1 1 ⎢ 0 mλ2c2 ⎢ − mλ2 ⎢⎣ 0 − mλ2 − mλ2c1

0 0 0 0 0 0 1 − b2 0 1 0 0

0⎤ ⎥ 0⎥ 0⎥ b1⎥ ⎥ 0⎥ 1 ⎥⎦

(15)

If the longitudinal displacement in x − direction can be neglected, a further reduction of the SV and Eq. (15) to Z = [Y , Θz , Mz , Qy]T and

URB

⎡ 1 b1 ⎢ 0 1 ⎢ = ⎢− mλ2 (b1 − c1) λ2 [JzI − m (b 2 c2 + b1 c1)] ⎢⎣ − mλ2 − mλ2c1

0 0 1 0

0⎤ ⎥ 0⎥ b1⎥ 1 ⎥⎦

(16)

may be achieved. 200

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

2.6. Euler-Bernoulli beam The full derivation of transfer matrices for Timoshenko- and Euler-Bernoulli beams may be found in [22]. For simplicity, only the transfer matrix for the Euler-Bernoulli beam in Fig. 2f will be presented here. The differential equation of an Euler-Bernoulli beam with length l , bending stiffness EI , mass per unit length m = ρA, material density ρ , and cross section area A and its modal transformation are given as

EI

y (x, t ) = Y (x ) e λt ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→

∂ 4y ∂ 2y +m 2 =0 ∂x 4 ∂t

∂ 4Y (x ) mλ2 + Y (x ) = 0 ∂x 4 EI

(17)

The general solution is Y (x ) = A1 cosh βx + A2 sinh βx + A3 cos βx + A4 sin βx with the eigenvalue-dependent abbreviation β = 4 −mλ2 / (EI ) and arbitrary constants A1 , ⋯A4 . For the Euler-Bernoulli beam, the linearized relations Θz = Y ′, Mz = EIY ″ and Qy = M ′z may be added to end up with the modal solution functions Z (x ) = B (x ) a or

⎡ cosh βx ⎤⎡ A ⎤ ⎡ Y (x ) ⎤ sinh βx cos βx sin βx ⎢ ⎥⎢ 1⎥ ⎢ ⎥ β βx β βx β βx β βx sinh cosh − sin cos Θ x ( ) ⎥ ⎢ A2 ⎥ ⎢ z ⎥ =⎢ ⎢ EIβ 2 cosh βx EIβ 2 sinh βx − EIβ 2 cos βx − EIβ 2 sin βx ⎥ ⎢ A3 ⎥ ⎢ Mz (x )⎥ ⎢ ⎥ ⎢ Q (x ) ⎥ ⎣ y ⎦ x ⎣ EIβ 3 sinh βx EIβ 3 cosh βx EIβ 3 sin βx − EIβ 3 cos βx ⎦ ⎢⎣ A4 ⎥⎦

(18)

]T

The coefficient vector a = [A1 , A2 , A3 , A4 of unknowns may be adapted to boundary condition ZI at input end as ZI = [B (0)] a or a = [B (0)]−1 ZI . For the beam output end at x = l we then get from Eq. (18) ZO = [B (l )] a = [B (l )] [B (0)]−1ZI = UBZI where

⎡S ⎢ ⎢ βV UB = ⎢ EIβ 2U ⎢ ⎢⎣ EIβ 3T S=

ch + c , 2

T /β S EIβV EIβ 2U T=

U / EIβ 2 T /EIβ S βV

sh + s , 2

U=

V / EIβ 3 ⎤ ⎥ U / EIβ 2 ⎥ ⎥, T /β ⎥ ⎥⎦ S ch − c , 2

V=

sh − s , 2

ch = cosh(βl ), sh = sinh(βl ), c = cos(βl ), s = sin(βl )

(19)

According to the definition given above, the coefficient β and thus the transfer matrix eigenvalue λ .

UB

is also an implicit function of the

3. Connection elements The main key of MSTMM is the simplicity of transferring a SV from one element to another. However, in some systems the dimensions of SVs may be inconsistent, like connecting the one-dimensional spring (3) to a planar rigid body (15) or (16). This makes it difficult to transfer from one SV to another of different dimension, since some state variables in a high-dimensional SV do not appear in a low-dimensional SV; mathematically, this prevents usual matrix multiplication of transfer matrices with nonmatching size. Another problem arises for branched systems where several transfer paths merge. Both problems can be resolved rather elegantly by introducing proper connection elements resulting in both transfer equations (TE) and geometrical equations (GE) between different input ends of that connection element according to the fundamentals of dynamics and kinematics. Generally

Fig. 3. Connection element with a) consistent and b) inconsistent SVs.

201

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

speaking, connection elements are massless dummy bodies with no size and no inertia. Let us first consider a situation like in Fig. 6b where two transfer paths with consistent SVs merge. Fig. 3a shows the associated connection element with two inputs I1 and I2 and states according to the sign convention in Fig. 1. In order to obtain the transfer equation for this element, we have to formulate kinematics and kinetics relations. Kinematical consistency and force/moment equilibrium for this massless element without size yields

YO = YI1;

Θz O = Θz I1;

With Z = [Y , Θz , Mz , Qy

]T

Mz O = Mz I1 + Mz I2;

Qy O = Qy I1 + Qy I2

(20)

this can be summarized in the TE

⎡1 ⎢ ZO = UI1 Z I1 + UI2 Z I2 where UI1 = ⎢ 0 ⎢0 ⎣0

0 1 0 0

0 0 1 0

0⎤ 0 ⎥, ⎥ 0⎥ 1⎦

⎡0 ⎢ UI2 = ⎢ 0 ⎢0 ⎣0

0 0 0 0

0 0 1 0

0⎤ 0⎥ ⎥ 0⎥ 1⎦

(21)

Additionally, the kinematics quantities at the two input ends coincide, i.e.,

YI1 = YI2;

Θz I1 = Θz I2

(22)

which may be shortly written as GE

⎡ ⎤ H I1 Z I1 = H I2 Z I2 where H I1 = H I2 = ⎢ 1 0 0 0 ⎥ ⎣0 1 0 0⎦

(23)

The same approach may also be applied to connection elements where SVs are of inconsistent dimension at input end I2 compared to input I1 and output end O as shown in Fig. 3b. For example, the SVs Z I1 and ZO may be beam connections like in Fig. 7b with 4 × 1-state vector Z = [Y , Θz , Mz , Qy]T , while input end I2 is connected either to a linear translational spring with 2 × 1-state vector Z = [Y , Qy]T or a torsional spring with 2 × 1-state vector Z = [Θz , Mz ]T . Although this element has non-uniform SVs, the relations between input and output states may be derived similarly to Eqs. (20)–(23). According to vanishing size and inertia, we obtain in case of a translational element input I2

YO = YI1; Θz O = Θz I1; Mz O = Mz I1; Qy O = Qy I1 + Qy I2

(24)

resulting in

ZO = UI1 Z I1 + UI2 Z I2

⎡1 ⎢ where UI1 = ⎢ 0 ⎢0 ⎣0

0 1 0 0

0 0 1 0

0⎤ 0 ⎥, ⎥ 0⎥ 1⎦

0⎤ 0⎥ ⎥ 0⎥ 1⎦

⎡0 ⎢ UI2 = ⎢ 0 ⎢0 ⎣0

(25)

The kinematical consistency of displacements at the two input ends can be written as YI1 = YI2 or as GE

H I1 Z I1 = H I2 Z I2 where H I1 = [1 0 0 0 ] and H I2 = [1 0 ]

(26)

Similarly, in case of a torsional input Θz I2 and Mz I2 instead of YI2 and Qy I2 at I2 , see right part of Fig. 3b, Eqs. (24)–(26) will take the following forms

YO = YI1; Θz O = Θz I1; Mz O = Mz I1 + Mz I2; Qy O = Qy I1, Θz I1 = Θz I2

(27)

resulting in

⎡1 ⎢ ZO = UI1 Z I1 + UI2 Z I2 where UI1 = ⎢ 0 ⎢0 ⎣0

0 1 0 0

0 0 1 0

⎡0 0⎤ 0 ⎥, U = ⎢ 0 ⎥ I2 ⎢ 0⎥ ⎢0 ⎣0 1⎦

0⎤ 0⎥ ⎥ 1⎥ 0⎦

(28)

and

H I1 Z I1 = H I2 Z I2 where H I1 = [ 0 1 0 0 ] and H I2 = [1 0 ]

(29)

Fig. 4. Cantilever beam with eccentric tip mass (a) and (b) corresponding state vectors and transfer direction.

202

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

4. Application of MSTMM to various beam problems 4.1. Beam with rigid body In order to demonstrate the general concept of MSTMM, let us start with a simple problem with chain topology taken from [10– 12]. Fig. 4a shows the model consisting of a uniform Euler-Bernoulli cantilever beam and a rigidly mounted mass m with moment of inertia Jzc about its centroid with eccentricity e . This vibrating system comprises two elements, beam 1 and rigid body 2, Fig. 4b. Choosing the transfer direction from left to right, the element TEs read as Z2 = U1Z1 and Z3 = U2 Z2 . If we consider only lateral vibrations, the SVs (1) reduce to Zi = [ Y , Θz , Mz , Qy]Ti , and the transfer matrices are U1 = UB according to Eq. (19) and U2 = URB according to Eq. (16) where b1 = d , c1 = e, b 2 = c2 = 0, JzI = Jzc + m e 2 . The overall TE results from Z3 = U2 Z2 = U2 U1 Z1 as

Uall Zall = 0 where Uall = [ U2 U1 − I 4×4 ],

T

Zall = [ Z1T ZT3 ]

(30)

with 4 × 4 -identity matrix I 4×4 . It should be noted that this overall TE involves only the boundary SVs Z1 and Z3, whereas SVs of internal connection points are eliminated. Thus, even for longer chains composed of several elements, the size of the overall transfer matrix would remain the same, which is a clear advantage of MSTMM over classical multibody formalisms where the final dynamics equations depend on all degrees of freedom. For common boundary conditions listed in Table 1, always half of the state variables are zero due to constraints, whereas the others are unknown. In our case of a fixed boundary 1 and free boundary 3, we have Zall = [0, 0, Mz1, Qy1, Y3, Θz3, 0, 0]T . Thus, Eq. (30) may be reduced to Uall Zall = 0 where Zall = [Mz1, Qy1, Y3, Θz3]T is composed of the unknown state variables only, and Uall is a 4 × 4 -square matrix resulting from elimination of all columns of Uall associated with zeros in Zall (i.e., Uall is composed of columns 3– 6 of Uall given in Eq. (30)). According to Eqs. (16) and (19), Uall is a function of the unknown eigenvalues λ of the system, where due to missing damping we get λ = ± i ω or λ2 = − ω 2 . Thus, we finally end up with Uall (ω) being a function of the unknown system eigenfrequencies. For non! 0 has to be fulfilled. The natural frequencies ω of the trivial solutions Zall ≠ 0 , the characteristic equation Δ (ω) = det Uall (ω) = system may now be computed by zero search of the determinant, which is based on sign change of Δ (ω) by scanning an interesting frequency range [0, Ω]. However, this procedure can be cumbersome due to a high amount of enumeration effort, and often zeros may be missed, especially if zeros are very close. In case of multiple roots there may be even no sign change at all. Therefore, a reliable and efficient algorithm called recursive scanning approach was developed in [23] by switching from zero search for Δ (ω) to minimization of its absolute value Δ (ω) . Numerical results are obtained for an aluminum beam ( ρ = 2700 kg/m3, E = 69 GN/m2 ) with length l = 1 m and circular crosssection of diameter 2r = 0.01 m , resulting in m = 0.212058 kg/m and EI = 33.8703 Nm2 . The fMin1D algorithm in [23] with frequency range [1500] rad/s, Nx0 = 100 grid size, and absolute precision tolerance ε = 10−6 is used to find the first five natural frequencies ωi . For comparison with [10–12], the tip mass characteristics are normalized by α = m / m l and β = Jzc / m l 3, and the natural frequencies are normalized as

Ωi = ωi / EI / m l 4

(31)

= 44.435 rad / s of the free cantilever beam without tip mass, respectively. Results for or w.r.t. the first eigenfrequency α = 1, β = 1, e = 0.4l , d = 0.8l according to [10] and [11] (see also Fig. 5a for the first two frequencies) and for α = 10, β = 0.5, e = 0.1l , d = 0.3l according to [12] (Fig. 5b) are shown in Table 2 with very good agreement with values reported in the literature. Once the natural frequencies are determined, one may generate corresponding free vibration mode shapes starting from Zall representing boundary displacements, rotations, moments and shear forces (see Ref. [22] for more details). ω1B

4.2. Beam with elastic restraints Next we consider a model investigated in [13], Fig. 6a, requiring a connection element with consistent SVs. It consists of a uniform Bernoulli-Euler beam with elastically restrained ends, where restraints are provided by translational and rotational springs with stiffnesses K y1, K y2 and K ′z1 , K ′z2 , respectively. Fig. 6b shows a system representation comprised of 5 elements with 4 boundaries, namely Z1, Z3, Z 7 (tips) and Z8 (root), where 2 and 4 are connection elements with consistent SVs analogously to Fig. 3a and described by Eqs. (21) and (23). Since only translation in y − direction and rotation are of interest, the SVs are chosen as Table 1 Common boundary conditions for systems with transversal beam vibrations. Support type

Zero terms

Unknown states

Fixed

Y , Θz Y , Mz

Θz , Qy

Pinned

Mz , Qy

Free

Mz , Qy

Guided

Θz , Qy

203

Y , Θz Y , Mz

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

Fig. 5. Determinant of Uall (ω) for uniform cantilever beam with tip mass: (a) α = 1, β = 1, e = 0.4l, d = 0.8l and (b)α = 10, β = 0.5, e = 0.1l, d = 0.3l .

Table 2 Normalized natural frequencies of a cantilever beam with eccentric tip mass obtained from MSTMM in comparison with [10–12]. i

α = 1, β = 1, e = 0.4l, d = 0.8l

Ωi Table 2 [10]

1 2 3 4 5

0.85068 1.98013 4.94508 7.99216 11.0965

MSTMM

0.8506782 1.9801295 4.9450790 7.9921612 11.096517

α = 10, β = 0.5, e = 0.1l, d = 0.3l

Ωi Table 2 [11]

0.72365 3.92091 24.45381 63.87463 123.13274

MSTMM

0.7236534 3.9209128 24.453806 63.874640 123.13269

ωi /ω1B

ωi /ω1B

Table 1 [12]

MSTMM

0.128191 0.938152 6.5233 17.655 34.486

0.12819089 0.93815182 6.52324547 17.6546646 34.4837954

Zi = [Y , Θz , Mz , Qy]Ti , i = 1, …, 8. The transfer matrix U3 = UB of beam element 3 is given by Eq. (19). Elements 1 and 5 are

Fig. 6. Free-free beam with elastic restraints (a) and (b) associated state vectors and transfer directions.

composed of translational and torsional springs. By composing Eqs. (3) and (4) according to the elements of the SV, we get

⎡1 ⎢ Uj = ⎢ 0 ⎢ ⎢0 ⎣0

− 1/ K yj ⎤ ⎥ 1 1/ K z′j 0 ⎥, ⎥ 0 1 0 ⎥ 0 0 1 ⎦ 0

0

j = 1, 5 (32)

Now, we are in position to set up the system's global TE and GEs. Following all transfer paths in Fig. 6b in opposite direction starting from the root, the main TE can be obtained as 204

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

Z8 = U4, I1 Z5 + U4, I2 Z6

= U4, I1 (U3Z4) + U4, I2 (U5Z 7)

= U4, I1 U3 (U2, I1 Z1 + U2, I2 Z2) + U4, I2 U5Z 7

= U4, I1 U3U2, I1 Z1 + U4, I1 U3U2, I2 U1Z3 + U4, I2 U5Z 7

(33)

For connection elements 2 and 4, the GEs are

H2, I1 Z1 = H2, I2 (U1Z3)

(34a)

and

H 4, I1 U3 (U2, I1 Z1 + U2, I2 (U1Z3)) = H 4, I2 U5Z 7 or

H 4, I1 U3U2, I1 Z1 + H 4, I1 U3U2, I2 U1Z3 = H 4, I2 U5Z 7

(34b)

Combining these equations yields the overall TE Uall Zall = 0 where

⎡ T14×4 T2 4×4 T34×4 −I 4×4 ⎤ ⎢ ⎥ T , ZTall = [ Z1T ZT3 ZT7 ZT8 ] , Uall = ⎢ G12×4 G 22×4 O2×4 O2×4 ⎥ 16×1 ⎢⎣ G 3 ⎥⎦ G G O 4 5 2×4 2×4 2×4 2×4 8×16 T1 = U4, I1 U3U2, I1, T2 = U4, I1 U3U2, I2 U1, T3 = U4, I2 U5, G1 = −H2, I1, G 2 = H2, I2 U1 , G 3 = −H 4, I1 U3U2, I1, G4 = −H 4, I1 U3U2, I2 U1, G5 = H 4, I2 U5

(35)

The transfer and geometric matrices of elements 2 and 4, namely U2, I1, U4, I1, U2, I2 , U4, I2 , H2, I1, H 4, I1 and H2, I2, H 4, I2 are given by Eqs. (21) and (23), respectively. Applying the boundary conditions according to Fig. 6b and Table 1 yields Uall8×8 (ω) Zall8×1 = 0 . According to [13], the beam parameters are chosen as l = 1 m , EI = 63476.1250 Nm2 and m = 15.31526 kg/m , the stiffnesses are K y1 = K y2 = 100 EI / l 3 and K ′z1 = K ′z2 = EI / l . After applying the fMin1D algorithm on search interval ω ∈ [1, 3000] rad/s with Nx0 = 100 and ε = 10−6 , results (31) according to Table 3 are obtained with very good agreement to [13] where authors suggest a method based on Fourier series for accurate vibration frequency computation of beams. 4.3. Beam with attached lumped mass and elastic support The pinned-pinned beam in Fig. 7a with additional linear spring of stiffness Ky and attached mass m according to [14] is used as an example for connecting inconsistent SVs. The beam is considered as a three-part system shown in Fig. 7b; namely the span left to the spring 0 ≤ x ≤ xs , between spring and lumped mass, xs ≤ x ≤ xm , and right to the lumped mass, xm ≤ x ≤ l . The boundary ends are Z1, Z4 , Z8 (tips) and Z10 (root). Both elements 2 and 5 are connection elements with SVs Z I1 = [Y , Θz , Mz , Qy]T and Z I2 = [Y , Qy]T of inconsistent dimension at their two input ends and SV ZO = [Y , Θz , Mz , Qy]T at the output end. These two connection elements are similar to Fig. 3b. A deduction of the transfer equation along Fig. 7b analogously to Eqs. (33)–(35) results in the system overall dynamics

Uall

⎡ T14×4 T2 4×2 T34×2 −I 4×4 ⎤ ⎢ ⎥ T ]T = ⎢ G11×4 G 21×2 O1×2 O1×4 ⎥ , ZTall = [ Z1T ZT4 Z8T Z10 , 12×1 ⎢⎣ G 3 ⎥ 1×4 G 41×2 G 51×2 O1×4 ⎦6×12

T1 = U7 U5, I1 U4 U2, I1 U1, T2 = U7U5, I1 U4 U2, I2 U3, G1 = −H2, I1 U1, G 2 = H2, I2 U3 , G 3 = −H5, I1 U4 U2, I1 U1, G4 = −H5, I1 U4 U2, I2 U3,

T3 = U7U5, I2 U6 , G5 = H5, I2 U6

(36)

Here in, U1, U4 and U7 represent beam transfer matrices (19) where lengths l are set to xs , xm − xs and l − xm , respectively; U3 is U6 the linear translational spring transfer matrix (3), is the lumped transfer matrix (5), and U2, I1, U5, I1, U2, I2 , U5, I2 , H2, I1, H5, I1, H2, I2 and H5, I2 are given by Eqs. (25) and (26), respectively. Applying the boundary conditions given by Fig. 7a according to Table 1 yields Z1 = [0, Θz1, 0, Qy1]T , Z4 = [0, Qy 4]T , Z8 = [Y8, 0]T , and Z10 = [0, Θz10 , 0, Qy10]T . Thus, the overall transfer matrix Uall reduces to a 6 × 6 matrix. For comparison with Table 2 in [14], it is assumed that the cross-section of the beam is rectangular with width b = 0.01 m and height h = 0.01 m , l = 1 m , ρ = 7.8 × 10 3 kg/m3 and E = 200 GN/m2 resulting in m = ρ × (b × h ) = 0.78 kg/m and Table 3 Non-dimensional natural frequencies

1 2 3

Ωi in comparison with [13] of a free-free beam with elastic restraints. [13] Tables (1–3)

MSTMM

3.031 4.681 6.216

3.02967 4.66396 6.16578

205

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

Fig. 7. Pinned-pinned beam with spring and lumped mass (a) and (b) state vectors and transfer directions.

Table 4 Non-dimensional natural frequencies Ωi in comparison with [14] of a pinned-pinned system with attached spring and lumped mass. [14] Table 2

MSTMM

Analytical method

Galerkin's technique

1 2 3

9.670 38.540 86.343

9.671 38.541 86.353

9.67048 38.54018 86.34348

EI = E × (b × h3 /12) = 166.67 N m2 . The attached mass, spring constant and their distances are Γ = m / m l = 0.1, ξs = xs / l = 0.3, ξm = xm / l = 0.4 and k = Ky l 3 / EI = 10 . By applying the fMin1D algorithm on search interval [12, 000] rad/s with Nx0 = 100 , ε = 10−6 results in the first 3 natural normalized frequencies (31) as tabulated in Table 4, which shows the excellent agreement between MSTMM and [14]. 4.4. Beam with intermediate support and elastically attached mass The system in Fig. 8a requires a combination of connection elements of mixed type. It consists of a uniform two-span pinned– pinned beam with intermediate support at x1 and a spring-mass system at x2 with stiffness Ky and lumped mass m [15]. According to Fig. 8b, the vibration system has 4 boundaries, namely tips Z1, Z3 and root Z10 where Z = [Y , Θz , Mz , Qy]T , and tip Z8 = [Y , Qy]T . The beam is split into three parts with lengths x1, x2 − x1 and l − x2 , respectively, and U1, U3 and U7 are represented by Eq. (19). U5 and U6 are related to the translational spring (3) and lumped mass (5), respectively. Additionally two connection elements are introduced where element 2 has consistent SVs similar to Fig. 3a and element 4 has inconsistent SVs similar to Fig. 3b. Sweeping through Fig. 8b opposite to the transfer directions yields

Z10 = U7 (U4, I1 (U3Z4) + U4, I2 (U5Z 7))

= U7 (U4, I1 (U3 (U2, I1 (U1Z1) + U2, I2 Z3)) + U4, I2 U5U6 Z8)

= U7U4, I1 U3U2, I1 U1 Z1 + U7U4, I1 U3U2, I2 Z3 + U7U4, I2 U5U6 Z8    

 

   

T1 4×4

T2 4×4

(37)

T3 4×2

The GEs for connection elements 2 and 4 are

− H2, I1 U1 Z1 +  

G1 2×4

H2, I2 Z3 = 0

⏟

G 2 2×4

− H 4, I1 U3U2, I1 U1 Z1 −H 4, I1 U3U2, I2 Z3 + H 4, I2 U5U6 Z8 = 0      

 

G3 1×4

G4 1×4

(38)

G5 1×2

where H2, I1 = H2, I2 is given by Eq. (23), where as H 4, I1,H 4, I2 are defined by Eq. (26). Combining Eqs. (37) and (38) results in an overall TE similar to Eq. (36) with appropriate dimensions of the sub-matrices. The eigenvalue search on the 7 × 7-matrix Uall is performed with the fMin1D algorithm on ω ∈ [120, 000] rad/s where Nx0 = 500 , ε = 10−6 . According to [15], the beam has a circular cross-section of diameter 2r = 0.05 m resulting in I = 3.06796 × 10−7m4 . The remaining system parameters are l = 1 m , x1 = 0.4 m , x2 = 0.75 m , E = 2.069 × 1011 N/m2 , m = 15.3875 kg/m , Ky = 190428.375 N/m and m = 3.0775 kg . The lowest five natural frequencies of the loaded beam are compared to [15] in Table 5. It is seen that very good agreement between the corresponding results is achieved. 4.5. Beam with damped oscillator The last example is a damped system presented in [16] and shown in Fig. 9a. It consists of a simply supported, uniform EulerBernoulli beam with grounded linear translational spring K y1 = 5EI / l 3 at position x1 = 0.2l , a damped oscillator with 206

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

Fig. 8. Two-span beam carrying intermediate spring-lumped mass system (a) and (b) corresponding elements.

Table 5 Natural frequencies ωi in comparison with [15] of two-span beam with one intermediate support and attached spring-mass system. [15] Table 3

MSTMM

NAM

FEM

247.6358 2156.8778 4938.1451 8220.1624 15847.8928

247.6360 2156.8780 4938.1520 8220.1982 15848.1531

247.63609 2156.87786 4938.14344 8220.16456 15847.8993

m = 0.12 m l , C = EI m / l 2 and K y2 = 4EI / l 3 at position x2 = 0.6l , and a grounded linear torsional spring K ′z = 10EI / l at position x3 = 0.75l . The beam parameters are l = 1 m , ρ = 7.8 × 10 3 kg/m3, A = 1.96349 × 10−3m2 , E = 2.069 × 1011 N/m2 and I = 3.06796 × 10−7m2 . Fig. 9b shows 5 boundaries with different dimensions and mixed state variables, namely Z1 and Z15 (root) in the form of Z = [Y , Θz , Mz , Qy]T , Z4 and Z9 in the form of Z = [Y , Qy]T , and Z13 in the form of Z = [Θz , Mz ]T . However, with the connection elements in Section 3 they may easily be combined. Elements 1, 4, 8, and 11 are beams with lengths x1, x2 − x1, x3 − x2 and l − x3, and their transfer matrices U1, U4 , U8 and U11 are given by Eq. (19). Element 3 is a linear translational spring with transfer matrix U3 according to Eq. (3), whereas U10 for the torsional spring is given by Eq. (4). Element 6 combines a linear translational spring and a dashpot. Thus, the transfer matrix combines the transfer matrices (3) and (6) as [20]

⎡1 − 1 ⎤ Ky + Cλ ⎥ U6 = ⎢ ⎢⎣ 0 1 ⎥⎦

(39)

representing qy O = qy I = Ky ( yI − yO) + C ( yİ − yȮ ). The lumped mass element 7 has a transfer matrix U7 according to Eq. (5). Elements 2, 5 and 9 are massless connection elements with inconsistent SVs analogously to Fig. 3b. Here, U2, I1, U5, I1, U2, I2 , U5, I2 , H2, I1, H5, I1, H2, I2 and H5, I2 are matrices according to Eqs. (25) and (26), while U9, I1, U9, I2 , H 9, I1 and H 9, I2 are matrices given by Eqs. (28) and (29). Following the transfer directions in Fig. 9b and adding the GEs of connection elements 2, 5 and 9 yields the 7 × 14 overall transfer matrix and SV

Uall

⎡ U11U9, I1 U8U5, I1 U4 U2, I1 U1 U11U9, I1 U8U5, I1 U4 U2, I2 U3 U11U9, I1 U8U5, I2 U6 U7 U11U9, I2 U10 − I 4×4 ⎤ ⎢ ⎥ − H2, I1 U1 H2, I2 U3 01×2 01×2 01×4 ⎥ = ⎢⎢ − H5, I1 U4 U2, I1 U1 − H5, I1 U4 U2, I2 U3 H5, I2 U6 U7 01×2 01×4 ⎥ ⎢ ⎥ 01×4 ⎦ ⎣ − H 9, I1 U8U5, I1 U4 U2, I1 U1 − H 9, I1 U8U5, I1 U4 U2, I2 U3 − H 9, I1 U8U5, I2 U6 U7 H 9, I2 U10

T ZTall = [[Y , Θz , Mz , Qy]1 [Y , Qy]4 [Y , Qy]9 [Θz , Mz ]13 [Y , Θz , Mz , Qy]15]14×1

7×14

(40a)

(40b)

Imposing the boundary conditions in Table 1, Eqs. (40) reduce to a 7 × 7 matrix Uall with complex eigenvalues λ = −δ ± iω due to damping. Therefore, we have to use the fMin2D algorithm from [23] where we may choose the search area ω ∈ [−10, 12000] rad/s, δ ∈ [−5300], grid size Nx 0 × Ny0 = 40 × 20 , and tolerances εx = εy = 10−6 . The first five eigenvalues found by MSTMM are shown in Fig. 10, which agree very well with two approaches described in [16], see Table 6.

207

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

Fig. 9. Simply supported beam with additional translational and rotational spring support and attached damped oscillator (a), as well as (b) corresponding elements.

Fig. 10. Determinant of Uall and first five eigenvalues of the simply supported beam with several attachments.

Table 6 Non-dimensional eigenvalues in comparison of MSTMM with [16] for simply supported beam with several attachments. [16] Table 9

MSTMM

λ / EI /(m l 4 )

Eq. (14)

FEM

λ1

−4.2116 × 10 0 + 4.1765 × 10 0i

−4.2128 × 10 0 + 4.1794 × 10 0i

−4.212833 × 10 0 + 4.179366 × 100i

λ2

−8.1376 × 10−1 + 1.2746 × 101i

−8.1473 × 10−1 + 1.2652 × 101i

−8.147330 × 10−1 + 1.265174 × 101i

λ3

−3.3614 × 10−1 + 3.9540 × 101i

−3.3615 × 10−1 + 3.9540 × 101i

−3.361446 × 10−1 + 3.953997 × 101i

λ4

−2.9362 × 10−1 + 9.2496 × 101i

−2.9602 × 10−1 + 9.2374 × 101i

−2.960221 × 10−1 + 9.237342 × 101i

λ5

−1.0244 × 10 0 + 1.6635 × 10 2i

−1.0223 × 10 0 + 1.6602 × 10 2i

−1.022282 × 10 0 + 1.660219 × 10 2i

5. Conclusions The paper demonstrates that MSTMM can easily handle various coupled problems with beams, dashpots, springs and rigid bodies in a unified way which have been discussed in the literature on all together 94 pages. Within a few lines only, transfer equations can be directly deduced from the system drawing based on the element transfer matrices taken from a MSTMM library. By defining the state vectors properly, the overall transfer matrix can be deduced for both chain and branched system topologies, where in the latter case additional kinematical relations are formulated. We always end up with a system of linear algebraic equations for the modal vibration coordinates of boundary states. A recursive scanning approach for the minima of the absolute determinant of the overall, reduced transfer matrix then results in the desired frequency and damping information which can be computed with any desired accuracy. The broad band of applications to a couple of problems taken from literature demonstrates simplicity and high efficiency of

208

Mechanism and Machine Theory 107 (2017) 197–209

L.K. Abbas et al.

modeling with a uniform MSTMM approved over the specific proposals in the literature. Although these examples have been only planar vibration problems, there is no inherent difficulty in extending the current method to spatial systems. Acknowledgments The research was supported by the Research Fund for the Doctoral Program of Higher Education of China (20113219110025, 20133219110037), the Natural Science Foundation of China Government (11102089, 61304137), the Program for New Century Excellent Talents in University (NCET-10-0075) and the Research Innovation Program 2013 for Graduates in Common Universities of Jiangsu Province (CXLX13-203). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

W. Schiehlen, Multibody Systems Handbook, Springer-Verlag, Berlin, Germany, 1990. J. Wittenburg, Dynamics of Systems of Rigid Bodies, B.G. Teubner, Stuttgart, Germany, 1977. W. Schiehlen, Multibody system dynamics: roots and perspectives, Multibody Syst. Dyn. 1 (1997) 149–188. T.R. Kane, P.W. Likins, D.A. Levinson, Spacecraft Dynamics, McGraw-Hill Book Company, New York, USA, 1983. J. Wittenburg, Dynamics of Multibody Systems, Springer-Verlag, Berlin, Germany, 2008. A.A. Shabana, Dynamics of Multibody Systems, Cambridge University Press, New York, USA, 2010. S. IIanko, Transcendental dynamic stability functions for beams carrying rigid bodies, J. Sound Vib. 279 (2005) 1195–1202. C.A. Rossit, P.A.A. Laura, Free vibration of a cantilever beam with a spring-mass system attached to the free end, Ocean Eng. 28 (2001) 933–939. J.R. Banerjee, Free vibration of beams carrying spring-mass systems - a dynamic stiffness approach, Comput. Struct. 104–105 (2012) 21–26. C.W.S. To, Vibration of a cantilever beam with a base excitation and tip mass, J. Sound Vib. 83 (1982) 445–460. N. Popplewell, D. Chang, Free vibration of a complex Euler-Bernoulli beam, J. Sound Vib. 190 (1996) 852–856. C.F.T. Matt, Simulation of the transverse vibrations of a cantilever beam with an eccentric tip mass in the axial direction using integral transforms, Appl. Math. Model. 37 (2013) 9338–9354. H.K. Kim, M.S. Kim, Vibration of beams with generally restrained boundary conditions using Fourier series, J. Sound Vib. 245 (2001) 771–784. A.D. Mohammad, K. Siavash, H.G. Megen, Free vibration of beam-mass-spring systems: analytical with numeric confirmation, Acta Mech. Sin. 28 (2012) 468–481. Y.L. Hsien, C.T. Ying, Free vibration analysis of a uniform multi-span beam carrying multiple spring–mass systems, J. Sound Vib. 302 (2007) 442–456. D.C. Philip, B.S. Andrew, Eigenvalue sensitivities of a linear structure carrying lumped attachments, J. Am. Inst. Aeronaut. Astronaut. 49 (2011) 2470–2481. E.C. Pestel, F.A. Leckie, Matrix Method in Elastomechanics, McGraw-Hill Book Company, New York, USA, 1963. R. Uhrig, Elastostatik and Elastokinetik in Matrizenschreibweise - Das Verfahren Der Uebertragungsmatrizen, Springer-Verlag, Berlin, Germany, 1973. X.T. Rui, G.P. Wang, Y.Q. Lu, L.F. Yun, Transfer matrix method for linear multibody system, Multibody Syst. Dyn. 19 (2008) 179–207. X.T. Rui, L.F. Yun, Y.Q. Lu, B. He, G.P. Wang, Transfer Matrix Method Of Multibody System And Its Application, Science Press, Beijing, China, 2008 (in Chinese). X.T. Rui, J.S. Zhang, Q.B. Zhou, Automatic deduction theorem of overall transfer equation of multibody system, Adv. Mech. Eng. 2014 (2014) Article ID 378047. L.K. Abbas, M.J. Li, X.T. Rui, Transfer matrix method for the determination of the natural vibration characteristics of realistic thrusting launch vehicle - Part I, Math. Probl. Eng. 2013 (2013) Article ID 764673. D. Bestle, L.K. Abbas, X.T. Rui, Recursive eigenvalue search algorithm for transfer matrix method of linear flexible multibody systems, Multibody Syst. Dyn. 32 (2014) 429–444.

209