7.2 Formulation

Back

7.2.1 Differential Formulation

In a majority of engineering vibration problems,

the amplitude of vibrations is very small, so that

the following assumptions hold: (1) a linear form

of strain – displacement relationships, and (2) a

linear form of stress – strain relationships (Hooke’s

Law). If the energy losses are negligible, it is

straightforward to apply Newton’s (second) law

and Hooke’s Law to derive the equations of

motion, which appear as differential equations.

Consider a single-DoF spring – mass system, as

shown in Figure 7.1. The two laws are given by

mu€ ðtÞ ¼ 2f ; Newton’s law

f ¼ kuðtÞ; Hooke’s Law

(

ð7:1Þ

The first equation describes the inertia force, and the second equation describes the elastic force. These

two forces are essential for mechanical vibration to exist (see Chapter 1).

In a similar way, the differential equations are given directly when Newton’s law plus Hooke’s Law is

applied to a multiple-DoF spring – mass system, shown in Figure 7.2

Mu€ ðtÞ ¼ 2F; Newton’s law

F ¼ KuðtÞ; Hooke’s Law

(

ð7:2Þ

where M is the (diagonal) mass matrix, and K is the stiffness matrix.

k

m

u(t)

FIGURE 7.1 Single-DoF spring – mass system.

7-2 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

In the case of continua, the differential equations of motion can be derived by means of Newton’s law

and Hooke’s Law in the same way as above. But in this case the boundary conditions have to be specified

in order to make the problem statement complete (see Chapter 4). For example, as a direct consequence

of Newton’s law and Hooke’s Law, the differential equation of bending vibration of a clamped – clamped

Euler beam may be given as (Chapter 4)

r

›2uðx; tÞ

dt2 ¼ 2f ; Newton’s law

f ¼ EI

›4uðx; tÞ

›x4

; Hooke’s Law

uð0; tÞ ¼ uðl; tÞ ¼

›uð0; tÞ

›x ¼

›uðl; tÞ

›x ¼ 0; Boundary conditions

8>>>>>>><

>>>>>>>:

ð7:3Þ

where r represents the mass per unit length, l the beam length, and EI the bending stiffness (flexural

rigidity).

7.2.2 Integral Formulation and Rayleigh – Ritz Discretization

Besides the approach in which Newton’s law and Hooke’s Law are directly used to establish equations of

motion, there are alternatives: Hamilton’s principle, the minimum potential energy principle, and the

virtual work principle; all of which appear in integral form. From a mathematical standpoint, the

differential equations and the integral equations are equivalent in that one can be derived from another.

However, they are very different in that the integral equations facilitate the application of the discretization

schemes such as the finite element method, an element-wise application of Rayleigh – Ritz method.

Therefore, Hamilton’s principle, as one of the integral formulations, and its Rayleigh – Ritz discretization

are briefly introduced here in order to provide a better understanding of the finite element method.

Denote T as the system kinetic energy, V the system potential energy, and dW the virtual work done by

nonconservative forces. Hamilton’s principle [11] states that the variation of the Lagrangian ðT 2 V Þ

Standard terminology plus the line integral of the virtual work done by the nonconservative forces

during any time interval must be equal to zero. If the time interval is denoted by ½t1; t2􀀉; then Hamilton’s

principle can be expressed as

d

ðt2

t1 ðT 2 V Þdt þ

ðt2

t1

dW dt ¼ 0 ð7:4Þ

In the case of a continuum, we look for an approximate solution uðx; y; z; tÞ in the form of

uðx; y; z; tÞ ¼

Xn

i¼1

wiðx; y; zÞqiðtÞ ð7:5Þ

where wiðx; y; zÞ is called a Rayleigh – Ritz shape function and qiðtÞ is called a generalized coordinate.

In this way, the system kinetic energy and the system potential energy can be, respectively, expressed

k1

m1

u1(t) u2(t) un(t)

m2

k2

mn

kn

FIGURE 7.2 Multiple-DoF spring – mass system.

Vibration Modeling and Software Tools 7-3

© 2005 by Taylor & Francis Group, LLC

as follows:

T ¼

1

2

Xn

i¼1

Xn

j¼1

mijq_iq_j ; 1

2

u_ TðtÞMu_ ðtÞ ð7:6Þ

and

V ¼

1

2

Xn

i¼1

Xn

j¼1

kijqiqj ; 1

2

uTðtÞKuðtÞ ð7:7Þ

where uTðtÞ ; ½q1; q2; …; qn􀀉; M ; ½mij􀀉; K ; ½kij􀀉:

The virtual work done by the generalized forces is

dW ¼

Xn

i¼1

fiðtÞdqi ¼ FT duðtÞ ð7:8Þ

where FT ; ½ f1ðtÞ; f2ðtÞ; …; fnðtÞ􀀉 and fiðtÞ is the generalized force corresponding to the nonconservative

force f ðx; y; z; tÞ

fiðtÞ ¼

ð

wiðx; y; zÞf ðx; y; z; tÞdv ð7:9Þ

Substituting Equation 7.6, Equation 7.7 and Equation 7.8 into Hamilton’s principle (Equation 7.4) and

conducting a routine variation operation, one has

ðt2

t1 ðu_ TM du_ 2 uTK du þ FT duÞdt ¼ 0 ð7:10Þ

Applying the separation integration to the first term of the above equation and noting that the variations

of the generalized coordinate du at times t1 and t2 equal zero, Equation 7.10 is rewritten as

ðt2

t1 ð2Mu€ ðtÞ 2 KuðtÞ þ FÞdu dt ¼ 0 ð7:11Þ

Because du; the variation of the generalized coordinate vector, is arbitrary and independent, from the

above equation one obtains

Mu€ ðtÞ þ Ku ¼ F ð7:12Þ

which is the vibration equation resulting from a Rayleigh – Ritz discretization.

7.2.3 Finite Element Method

In the finite element method (FEM) [7,10,12], a continuum is divided into a number of relatively small

regions called elements that are interconnected at selected nodes. This procedure is called discretization.

The deformation within each element is expressed by interpolating polynomials. The coefficients of

these polynomials are defined in terms of the element nodal DoF that describe the displacements and

slopes of selected nodes on the element. By using the connectivity between elements, the assumed

displacement field can then be written in terms of the nodal DoF by means of the element shape

function. Using the assumed displacement field, the kinetic energy and the strain energy of each

element are expressed in the form of the element mass and stiffness matrices. The energy expressions

for the entire continua can be obtained by adding the energy expressions of its elements. This leads to

the assembled mass matrix and the assembled stiffness matrix, and finally to the finite element

vibration equation.

The displacement in the interior of an element e is determined by a polynomial

uðx; y; z; tÞ ¼ Nue ð7:13Þ

7-4 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

where the matrix N is called the shape function matrix of the element e; and ue the vector of the

nodal DoF.

Based on the element displacement expression Equation 7.13, one can obtain the strain and the stress

in the element e and finally the strain energy.

The strain and the stress in the element e are

1 ¼ ›uðx; y; z; tÞ ¼ ›Nue ¼ Bue ð7:14Þ

and

s ¼ D1 ¼ DBue ¼ Sue ð7:15Þ

respectively, where › is the differential operator matrix, B ¼ ›N is the element strain matrix, D is the

elastic matrix, and S ¼ DB is called the element stress matrix.

The strain energy in the element e is given by the element strain 1 and stress s as

V e ¼

1

2

ð

1Ts dv ¼

1

2 ðueÞTKeue ð7:16Þ

where

Ke ¼

ð

BTDB dv ð7:17Þ

is called the element stiffness matrix.

The velocity at a point ðx; y; zÞ in the element e can be obtained from Equation 7.13 as

u_ ðx; y; z; tÞ ¼ Nu_ e ð7:18Þ

So the kinetic energy of the element e is

Te ¼

1

2

ð

ru_ T u_ dv ¼

1

2 ð_ ueÞTMe u_ e ð7:19Þ

where

Me ¼

ð

rNTN dv ð7:20Þ

is called the element mass matrix.

The equivalent nodal force Fe corresponding to the force f e applied onto the element e is determined

by equaling the work done by Fe to the work done force by f e along any virtual displacement. This leads to

the following:

ðdueÞTFe ¼

ð

duTf e dv ¼ ðdueÞT

ð

NTfe dv

􀀏 􀀐

ð7:21Þ

Note that as the variation of the nodal displacement is arbitrary, one can obtain the expression of the

equivalent nodal force Fe from Equation 7.21 as

Fe ¼

ð

NTfe dv ð7:22Þ

Now we have the kinetic energy, the strain energy, and the equivalent nodal force of the element e: But

these quantities are expressed in the local coordinate system ðXe; Y e; ZeÞ of the element e; not in the global

coordinate system ðX; Y ; ZÞ: In order to calculate the corresponding counterparts for the whole structure,

it is necessary to transform the expressions of the kinetic energy, the strain energy, and the equivalent

nodal force of the element e from the local coordinate system into the global one.

Let L be the transformation matrix from the global coordinate system to the local coordinate system.

Then the nodal displacement vector ue in the local coordinate system is related to the nodal displacement

Vibration Modeling and Software Tools 7-5

© 2005 by Taylor & Francis Group, LLC

vector u􀀊 e in the global coordinate system by the following:

ue ¼ Lu􀀊 e ð7:23Þ

Similarly, the equivalent nodal force vector Fe in the local coordinate system is related to the equivalent

nodal force vector F􀀊 e in the global coordinate system by

Fe ¼ LF􀀊 e ð7:24Þ

Substituting Equation 7.23 into Equation 7.16 and Equation 7.19, and noting that L is a normal matrix

ðLT ¼ L21Þ; the element stiffness and mass matrices in the global coordinate system can be, respectively,

expressed as

K􀀊 e ¼ LTKeL ð7:25Þ

and

M􀀊 e ¼ LTMeL ð7:26Þ

The equivalent nodal force vector in the global coordinate system is solved from Equation 7.24

F􀀊 e ¼ LTFe ð7:27Þ

In this way, we can obtain the total strain energy of the structure as

V ¼

X

e

V􀀊 e ¼

1

2

X

e ðu􀀊 eÞTK􀀊 e u􀀊 e ¼

1

2

uTKu ð7:28Þ

where the matrix

K ¼

X

e

K􀀊 e ð7:29Þ

is called the global stiffness matrix of the structure. The vector u is the global nodal displacement vector

of the structure.

Similarly, the total kinetic energy of all of the elements can be written as

T ¼

X

e

T􀀊 e ¼

1

2

X

e ð_􀀊ueÞTM􀀊 e_􀀊

ue ¼

1

2

u_ TMu_ ð7:30Þ

where the matrix

M ¼

X

e

M􀀊 e ð7:31Þ

is called the global mass matrix. The vector u_ is the global nodal velocity vector.

The total virtual work done by the external forces is

dW ¼

X

e

dW􀀊 e ¼

X

e ðdu􀀊 eÞTF􀀊 e ¼ ðduÞTF ð7:32Þ

where the vector

F ¼

X

e

F􀀊 e ð7:33Þ

is a global generalized force vector.

7-6 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

Substituting Equation 7.28, Equation 7.30, and Equation 7.32 into Hamilton’s principle (Equation 7.4)

and conducting the routine variation operation, one has

Mu€ þ Ku ¼ F ð7:34Þ

which is the vibration equation resulting from the finite element discretization.

7.2.4 Lumped Mass Matrix

The element mass matrix given by Equation 7.20 is normally a full symmetric matrix, because the

element shape functions are not orthogonal with each other. It is desirable to reduce this full matrix into a

diagonal matrix. In practice, this is achieved by lumping the element mass at its nodes. For example, the

consistent element mass matrix of a beam element is

Me ¼

rAl

420

156 22l 54 213l

22l 4l2 13l 23l2

54 13l 156 222l

213l 23l2 222l 4l2

2

6666664

3

7777775

ð7:35Þ

When the inertia effect associated with the rotational DoF is negligible, the element lumped mass matrix

can be obtained by lumping one half of the total beam element mass at each of the two nodes along the

translation DoF:

Me ¼

rAl

2

1 0 0 0

0 0 0 0

0 0 1 0

0 0 0 0

2

6666664

3

7777775

ð7:36Þ

When the inertia effect associated with the rotational DoF is not negligible, the mass moment of inertia of

one half of the beam element about each node can be computed and included at the diagonal locations

corresponding to the rotational DoF:

Me ¼

rAl

2

1 0 0 0

0 l2=12 0 0

0 0 1 0

0 0 0 l2=12

2

6666664

3

7777775

ð7:37Þ

7.2.5 Model Reduction

The finite element discretization of an engineering vibration problem usually generates a very large

number of DoF. In particular, when automatic meshing schemes are not properly applied, or threedimensional

elements must be used, the number of elements created could become too great to be costeffectively

handled with limited computer capabilities. To solve this problem, modelers have to pay

close attention to how the meshing is done in commercial software packages. Very often, simplification

and idealization based on the nature of the problem of concern can tremendously reduce the number

Vibration Modeling and Software Tools 7-7

© 2005 by Taylor & Francis Group, LLC

of elements. For example, there could be two ways of generating the finite elements of a clamped-free

steel beam with a metal block attached to its free end. One way is to mesh both the beam and the block

using three-dimensional elements; the other way is to mesh the beam with one-dimensional beam

elements and treat the block as a lumped mass, zero-dimensional element. It is obvious that the first

approach will result in many more elements than the second approach. However, both approaches will

give very similar results for the first several natural frequencies and the associated mode shapes. Another

technique for reducing the number of elements comes from deleting the detailed features. The detailed

features here imply those geometrical details, such as filets, chamfers, small holes, and so on, which do

not have significant contributions to the vibration behavior of the entire structure, but increase the

number of elements. Generally these detailed features can be deleted without any visible effect on the

results, if the global behavior of the vibration problem is of concern. Note that such detailed features may

have to be kept if the localized behavior such as fatigue (stress) induced by vibration is to be evaluated.

When further model reduction is necessary, Guyan reduction [3] is considered. It was proposed two

decades ago when computer capabilities were much more limited than today. In fact, Guyan reduction is

still in use today and has been cast into many commercial software packages. In Guyan reduction, the

model scale is reduced by removing those DoF (called slave DoF) that can be approximately expressed by

the rest of the DoF (called master DoF) through a static relation. The DoF associated with zero mass or

relatively small mass are likely candidates for slave DoF.

By rearranging the DoF u so that those to be removed, denoted by u2; appear last in the vector, and

partitioning the mass and the stiffness matrices accordingly, one obtains

M11 M12

M21 M22

" #

u€ 1

u€ 2

( )

þ

K11 K12

K21 K22

" #

u1

u2

( )

¼

F

0

( )

ð7:38Þ

If we assume M22 ¼ 0; and M21 ¼ 0; then the second equation in Equation 7.38 can be written as

u2 ¼ 2K21

22 K21u1 ð7:39Þ

Define the transformation

u ¼ Qu1 ð7:40Þ

where the transformation matrix Q is

Q ¼

I

2K21

22 K21

" #

ð7:41Þ

and I is the unit (identity) matrix.

Substituting Equation 7.40 into Equation 7.38 and premultiplying the resulting equation by QT; one

obtains a new reduced-order model

QTMQu€ 1 þ QTKQu1 ¼ F ð7:42Þ

where

QTMQ ¼ M11 2 M12K21

22 K21 þ KT

21K21

22 M22K21

22 K21 ð7:43Þ

and

QTKQ ¼ K11 2 K12K21

22 K21 ð7:44Þ

7-8 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC