9.4 Equations of Motion and Solution Methods

Back

In this section, an overview of the equation of motion and the various solution methods available are

provided. In most cases, detailed derivation of the equations is not given. The intention is to highlight

various solution schemes and the practical aspects of choosing a particular one for a given system or a

given dynamic analysis. A consideration of damping and its idealization are also discussed in this section.

9.4.1 Equation of Motion

The equations of motion describing the dynamic behavior of a structural system may be obtained from

the extended Hamilton’s principle for elastodynamics:

d

ðtf

to ðU 2 TÞdt 2

ðtf

to

dW dt ¼ 0 ð9:1Þ

where

d: first variation

to; tf : two arbitrary time points at which the first variation vanishes

U: strain energy of the system

T: kinetic energy of the system

dW: virtual work of the external forces acting on the system during virtual displacements

The external forces may be conservative or nonconservative. Nonconcentrative forces, e.g., damping

and follower forces, are deformation dependent, i.e., their magnitudes and/or directions depend on the

deformations. Introducing the finite element discretization assumption and substituting expressions for

the virtual work and the kinetic and strain energy terms into the above equation, we obtain the Lagrange

equation which finally leads to the equations of motion of the system that may be written in the following

form:

Mu€ þ Cu_ þ Ku ¼ pðtÞ ð9:2Þ

where M, C, and K are, respectively, the global mass, global damping, and global stiffness matrices (made

up by proper assembly of the elementmatrices); pðtÞ is the time-dependent applied force vector; and u€ ; u_ ;

and u are the nodal acceleration, nodal velocity, and nodal displacement vectors, respectively.

The global mass and damping matrices are assembled from the element matrices that are given by

(Cook et al., 1989; Bathe, 1996):

MðeÞ ¼

ð

V

rðeÞNTN dV ð9:3Þ

CðeÞ ¼

ð

V

cðeÞ s NTN dV ð9:4Þ

where rðeÞ and cðeÞ s are the mass density and the viscous damping coefficient for an element e; and N is the

matrix of the element shape functions.

As indicated above, there are two general approaches to solving Equation 9.1. One is the direct

integration approach and the other is the modal superposition approach. With the direct integration

approach, the equations of motion are integrated directly in the time domain. In the second approach,

the modal superposition method, the equations of motion are first transformed to modal generalized

displacements through the use of the mode shapes of the structure. The transformation yields a set of

Finite Element Applications in Dynamics 9-27

© 2005 by Taylor & Francis Group, LLC

decoupled second-order differential equations for the system that are easier to solve. In the remaining

part of this section, we describe the two solution approaches in more detail.

9.4.2 Direct Integration Method

With this method, the coupled equations of motion are integrated directly without special

transformation. Direct integration schemes seek to satisfy the equilibrium equations of motion at

discrete time points, 0; Dt; 2Dt; …; t; t þ Dt;… etc. This requires assumptions for the variation of u€ ; u_ ;

and u within the time step. These assumptions will determine the accuracy and convergence

characteristics of a particular time integration scheme. They also govern the variation of these quantities

with time and are different from the shape function assumptions for the element, which govern the

displacement variation within the element, i.e., with spatial coordinates. Direct integration schemes may

be categorized as “explicit” and “implicit” schemes. In explicit schemes, the equilibrium equations at

time t are used to solve the unknowns at time t þ Dt: These schemes are usually fast and efficient and

require no factorization of the effective stiffness matrix provided that the mass and damping matrices can

be formulated as diagonal matrices. The main drawback of these schemes is the requirement that the time

step for the analysis has to be smaller than a critical value. On the other hand, in implicit schemes, the

equilibrium equations at time t þ Dt are used to solve unknowns at time t þ Dt: These schemes require

factorization of the effective stiffness matrix of the structure at each time step but do not require the

condition of a critical time step. In implicit schemes, the number of operations at each time step is

proportional to the matrix order times half of the band width, m ¼ anb; where m is the number of

operations, n is the number of equations, a is a constant $ 2, and b is half the band width.

Many integration schemes are available in commercial FE programs. These include Euler forward,

Euler backward, central difference, Houbolt, Wilson-u, and Newmark schemes (Belytschko and Hughes,

1983; Zeng et al., 1992). In the following, one example of an explicit scheme and one example of an

implicit scheme are discussed.

9.4.2.1 The Central Difference Method

This is an explicit scheme in which the velocity and acceleration assumptions are given by

t u_ ¼

1

2Dt

􀀍

tþDt u 2 t2Dt u

􀀎

t u€ ¼

1

ðDtÞ2

􀀍

tþDt u 2 2 t u þ t2Dt u

􀀎

9>>=

>>;

ð9:5Þ

Using Equation 9.5 in the equilibrium equation at time (t) to obtain a solution at time ðt þ DtÞ gives

1

ðDtÞ2 M þ

1

2ðDtÞ

C

􀀒 􀀓

tþDt u ¼ t p 2 K 2

2

ðDtÞ2 M

􀀒 􀀓

t u 2

1

ðDtÞ2 M 2

1

2ðDtÞ

C

􀀒 􀀓

t2Dt u ð9:6Þ

If the mass and damping matrices are diagonal, Equation 9.6 requires no assembly and no factorization of

the effective stiffness matrix of the structure is required. The scheme is, however, conditionally stable and

requires a time step smaller than a critical value given by

Dtcr # Tn=p ð9:7Þ

REMARKS

* Equation 9.2 is the general equation of motion for a dynamic system. Two solution methods

are available: the direct integration approach and the modal superposition approach.

9-28 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

where Tn is the smallest period in the FE assembly with n DoF. This is an important condition and should

be examined carefully. If the mass matrix has a zero diagonal value, the corresponding period will be zero

and Equation 9.7 cannot be satisfied. Introducing very small diagonal values may not provide a practical

solution since the time step will be very small and the number of steps will be excessive. The scheme is

generally applied when a lumped mass assumption with no zero diagonal terms and velocity-dependent

damping are assumed. The initial conditions required for this scheme are the displacements, velocities,

and accelerations at time t ¼ 0:

9.4.2.2 The Newmark Method

This is a commonly used implicit scheme with the following assumptions:

tþDtu_ ¼ t u_ þ

h

ð1 2 aÞt u€ þatþDt u€

i

Dt ð9:8Þ

tþDt u ¼ t u þ t u_ Dt þ

1

2

2 b

􀀏 􀀐

t u€ þbtþDt u€

􀀒 􀀓

ðDtÞ2 ð9:9Þ

where a and b are user input parameters. Equation 9.8 and Equation 9.9 are solved for tþDtu_ and tþDtu€ ;

then substituted into the equilibrium equations at time ðt þ DtÞ to obtain:

K þ

a

bDt

C þ

1

bðDtÞ2 M

􀀒 􀀓

tþDt u ¼ tþDt p􀀊 ð9:10Þ

where the right-hand side is a function of the parameters a; b; Dt; K; M; C; tu; tu_ ; and tu€ : Equation 9.10

indicates that the scheme requires assembly and factorization of the effective stiffness matrix of the

structure. The choice of parameters a and b will determine the stability of the scheme and may reduce

the scheme to be equivalent to other known ones as follows:

* For 2b $ a $ 1=2; the scheme is unconditionally stable. This does not guarantee accuracy and

only ensures that the results will not grow out of bounds (Belytschko and Hughes, 1983).

* A slightly more stringent criterion for unconditional stability that provides artificial damping in

higher modes is given by using (Belytschko and Hughes, 1983):

a $

1

2

and b $

1

4

a þ

1

2

􀀏 􀀐2

* A commonly used option gives a ¼ 1=2 and b ¼ 1=4 which corresponds to a trapezoidal rule with

constant average acceleration.

* Another option that corresponds to a linear acceleration assumption is given by a ¼ 1=6 and

b ¼ 1=2: This method is very good for small Dt but tends to be unstable for large Dt:

Table 9.3 summarizes the above choices for a and b:

TABLE 9.3 Typical Choices of a and b Parameters for the Newmark Scheme

a Value b Value Comments

2b $ a $ 1=2 2b $ a $ 1=2 No guarantee of accuracy

a $

1

2

b $

1

4

􀀄

a þ

1

2

􀀅2 Improved accuracy

Provides artificial damping in higher modes

a $

1

2

b $

1

4

􀀄

a þ

1

2

􀀅2 Improved accuracy

Corresponds to a trapezoidal rule with constant average acceleration

a $

1

2

b $

1

4

Improved accuracy

Provides artificial damping in higher modes

Good for small Dt

a ¼ 1=6 b ¼ 1=2 Improved accuracy

Corresponds to a linear acceleration

Very good for small Dt but tends to be unstable for large Dt

Finite Element Applications in Dynamics 9-29

© 2005 by Taylor & Francis Group, LLC

9.4.3 Modal Superposition Method

In this method, the nodal displacement response is expressed in terms of the normal modes that may be

found in an eigenvalue analysis. The coupled equations of motion, Equation 9.2, are first transformed

into a set of independent or decoupled differential equations cast in modal generalized coordinates

(Mirovitch, 1980). The dynamic response of the original system is then obtained by superimposing the

responses of the uncoupled modal equations. The generalized coordinates (modal coordinates) are

introduced by the following coordinate transformation

u ¼ Fq and ui ¼

Xm

j¼1

wijqj ð9:11Þ

where F is an n £ m matrix called the eigenvector or mode shape matrix, m is the number of eigenvectors

ðm # nÞ; n is the number of DoF of the system, and q is a vector of size m £ 1 representing the number of

mode-amplitude generalized or normal coordinates. The transformation given by Equation 9.11

represents a change of basis from nodal displacements (u) to modal generalized coordinates (q). In order

for the equations to be decoupled in the modal generalized coordinate system, the triple product FTCF

has to be a diagonal matrix and the following orthogonality property is assumed for damping:

FTCF ¼ diagð2jiviÞ ð9:12Þ

where ji is the damping ratio for mode i and diagð2jiviÞ indicates a diagonal matrix with the ith diagonal

component of 2jivi: This form has been assumed by generalization of damping for a single DoF system.

It is practically more convenient to define the damping by the damping ratios of each mode than it is to

evaluate the damping matrix explicitly. Introducing Equation 9.11 and Equation 9.12 into the equation

of motion, the following decoupled equations of motion are obtained:

q€ þ diagð2jiviÞq_ þ diagðv2i

Þq ¼ FTpðtÞ ð9:13Þ

If the components of the applied force vector have the same time function, so that pðtÞ can be expressed as

pðtÞ ¼ p􀀊 gðtÞ ð9:14Þ

where p􀀊 is a constant vector and gðtÞ is a function giving the time change of the load vector, then the

modal load can be written as

fr ðtÞ ¼ ðwT

r p􀀊 ÞgðtÞ ¼ GrgðtÞ ð9:15Þ

REMARKS

* In the direct integration approach, the equations of motion are integrated directly without

special transformation. Such schemes may be categorized as explicit and implicit schemes.

* Explicit schemes are usually fast, efficient, and require no factorization of the effective

stiffness matrix provided that the mass and damping are diagonal (refer to Equation 9.6 for

an example). However, explicit schemes require the time step to be smaller than a critical

value (Equation 9.7) that normally leads to very large number of steps.

* Implicit schemes require factorization of the effective stiffness matrix of the structure at

each time step (refer to Equation 9.10 for an example) but do not require the condition of a

critical time step.

* Although implicit schemes are unconditionally stable, accuracy is not always guaranteed

without careful consideration of the step size (refer to Table 9.3).

9-30 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

where wT

r is a row modal vector corresponding to mode r: The quantity Gr is referred to as the modal

participation factor (NISA User’s Manual, 1992). A particular use of this quantity is in ground motion

types of excitation.

Once the generalized coordinates, q, are evaluated, the physical response of the original system in

terms of nodal displacements, velocities, accelerations, and stresses are recovered from Equation 9.11. For

the stress recovery, it is noted that (NISA User’s Manual, 1992; Kim et al., 1994):

sðtÞ ¼ DBu􀀊 ðeÞðtÞ ¼ DB FðeÞq ¼ ½DB FðeÞ􀀉 q ¼ Lq or siðtÞ ¼

Xm

r¼1

Lir qr ð9:16Þ

where siðtÞ represents a stress component in an element (i.e., sxx ; syy ; sxy ; etc.), Lir lists the stress

components corresponding to the rth mode at a typical point. That is, Lir ; r ¼ 1; 2; …; m are modal

stresses which should be interpreted as stress shapes rather than absolute values of stress and the matrix D

is the material stress – strain matrix.

The above derivation of the mode superposition method shows the decoupling advantage of the

normal coordinates, whereby the change of basis from nodal displacements to normal modes yields a set

of uncoupled modal equations, with each equation cast in the form of a single DoF oscillator. As

indicated above, another major advantage of the normal mode method is that a good approximation to

the response may be obtained using a drastically reduced number of coordinates ðm p nÞ: For most types

of loading, with the exclusion of shock and impact loading, the contributions of the lowest few modes are

generally more pronounced than the higher modes. Furthermore, practical finite element idealization

approximates the lower modes better and tends to be less reliable for higher modes of vibration.

9.4.4 Damping Formulation

Damping is a source of energy dissipation in the structure and it leads to decay of the free vibration

amplitude with time. Damping sources in a structure include internal friction in the material, Coulomb

friction in sliding and pin joints, and other viscous friction forces due to fluid friction (Beards, 1982;

Fretzen, 1986). The overall damping of a system is normally quite difficult to obtain and, in general, must

be determined experimentally (Kareem and Gurley, 1996). In free vibrations and modal analysis, the

damping may be neglected or an overall small value for the system may be assumed. This is a reasonable

assumption since in practice damping is small enough and can usually be assumed to be viscous. In

forced vibrations, however, and when the forcing frequency is close to one of the system’s natural

frequencies, the response of the system is dominated by the specified damping values.

REMARKS

* In the modal superposition approach, the equations of motion are first transformed into a

set of decoupled equations cast in modal generalized coordinates. The response of the

system is then obtained by superimposing the solutions of the decoupled modal equations

(Equation 9.13 and Equation 9.16).

* The modal superposition approach may only be achieved for linear systems with

proportional or directly assumed modal damping (refer to Section 9.4.4).

* In the modal superposition approach, a good approximation to the response may be

obtained using a drastically reduced number of modes. For most types of loading, with the

exclusion of shock and impact loading, the contributions of the lowest few modes are

generally more pronounced than the higher modes.

* Practical FE idealization approximates the lower modes better and tends to be less reliable

for higher modes of vibration.

Finite Element Applications in Dynamics 9-31

© 2005 by Taylor & Francis Group, LLC

Commonly, damping is described in viscous or structural form. Both descriptions are used for their

mathematical convenience but may not truly represent the actual damping behavior and mechanism. For

viscous damping, the damping force is proportional to and opposes the velocity. With structural

damping, the damping force is proportional in magnitude to the internal elastic force (i.e., to the

displacement) and is in the opposite direction to the velocity. In practical FE analysis, most programs

provide the following damping representation that may be used for dynamic analysis:

1. Discrete viscous damper elements (dashpots). As discussed in Section 9.1.2, these elements are

damping counterparts of the spring elements and are discrete idealization of viscous damping in the

structure. The damping matrix resulting from these elements, in general, cannot be decoupled as in

Equation 9.12. Hence, these elements may be used only with direct integration solution methods.

2. Proportional viscous damping (Rayleigh damping). In this type, the following arbitrary form for

the damping matrix is assumed:

C ¼ c1K þ c2M ð9:17Þ

where c1 and c2 are constants to be determined and supplied by the user. A commonly used method in

determining the constants c1 and c2 is to define two damping rations j1 and j2 corresponding to two

unequal natural frequencies v1 and v2; respectively. Since C is proportional to K and/or M, it satisfies the

orthogonality property and we have

c1v2i

þ c2 ¼ 2jivi ð9:18Þ

which leads to

ji ¼

c1vi

2 þ

c2

2vi ð9:19Þ

Equation 9.19 is then used to calculate the coefficient c1 and c2 from the knowledge of the specified

damping rations j1 and j2 and their corresponding natural frequencies v1 and v2: For the case when C is

only proportional to K ðc2 ¼ 0Þ; the damping ratio is proportional to the natural frequency, and thus the

higher modes will be damped more than the lower modes. Similarly, if C is only proportional to M

ðc1 ¼ 0Þ; the damping ratio is inversely proportional to the natural frequency and the higher modes will

have less damping than the lower modes. The relationship of Equation 9.19 is diagrammatically

illustrated in Figure 9.24.

The above method for calculating the coefficients c1 and c2 should not restrict the use of proportional

damping to only modal analysis. Proportional damping may be used in harmonic, modal, and transient

analyses, as well as substructure and reduced types of analyses.

For direct transient dynamic analysis, different sets of constants c1 and c2 may be assigned to different

parts of the model. This may result in a nonorthogonal global damping matrix. The damping

orthogonality condition, however, is not required in direct transient dynamic analysis since the governing

equations are directly integrated.

w

x

combined

x2

w1

x1

w2

stiffness proportional (c2 = 0, xi = )

2

c1wi

mass proportional (c1 = 0, xi = )

2wi

c2

FIGURE 9.24 Proportional damping.

9-32 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

3. Modal viscous damping, in which the damping ratios are directly specified for the

participating modes. This type of damping can be used in modal dynamic analyses and usually

provides good results for systems with small damping. In this type, the damping matrix is assumed

to be diagonal with values of 2jivi; where i is the ith diagonal and ji is the damping ratio for

mode i with frequency vi:

Typical values of overall damping ratios for various systems are given in Table 9.4 (Adams and

Askenazi, 1999).

Навигация