6.4 Finite Difference Method for a Continuous System

Back

6.4.1 Bar

A thin uniform rod is the first example of a continuous system that will be examined. The equation of

motion for such a system is given by

›2u=›x2 ¼ 1=c2 ›2u=›t2

where c ¼

ffiffiffiffiffi

pE=r: E is the modulus of elasticity and r is the density of the rod. A complete derivation of

this equation can be found in many books (see Thomson and Dahleh 1998). Separation of variables can

be used to find a solution of the form uðx; tÞ ¼ U ðxÞGðtÞ: This substitution gives the equation.

d2U

dx2 GðtÞ ¼

1

c2 UðxÞ

d2GðtÞ

dt2

Upon rearrangement, this becomes

1

UðxÞ

d2U

dx2 ¼

1

c2GðtÞ

d2GðtÞ

dt2

The left-hand side of the equation is independent of t and the right-hand side is independent of x:

Therefore, each side must equal a constant. Let this constant be denoted by a and a ¼ 2ðv=cÞ2: There are

now two differential equations, the displacement equation

d2U =dx2 þ a2U ¼ 0

Numerical Techniques 6-11

© 2005 by Taylor & Francis Group, LLC

and the temporal equation

d2G=dt2 þ v2G ¼ 0

The displacement equation can be solved using

a finite difference approximation. The finite

difference approximation requires dividing the

bar of length L into n 2 1 pieces each of length

h ¼ L=ðn 2 1Þ: The finite difference mesh consists of n points, 1; 2; 3; …; n: The first node (labeled 1)

is the left-hand endpoint of the bar, and the last node, labeled n; is the right-hand endpoint of the bar

(see Figure 6.5). The value of U at each point i is denoted by Ui: Using the centered difference

approximation for the second derivative, the discretized version of the displacement equation becomes

ðUiþ1 2 2Ui þ Ui21Þ=h2 þ aUi ¼ 0

Collecting like terms, the equation becomes

Uiþ1 2 ð2 2 lÞUi þ Ui21 ¼ 0

where l ¼ h2a2: As is typical for the finite difference approximation, this equation holds at the

interior mesh points ði ¼ 2; 3; …; n 2 1Þ and not at the endpoints. To see that this equation does not

hold at the endpoints, substitute i ¼ 1 into the recurrence relationship. The resulting equation is

U2 2 ð2 2 lÞU1 þ U0 ¼ 0

The problem is that there is no value U0: Recall that U1 represents the displacement of the lefthand

endpoint. Since U0 represents a point even further to the left, it is not on the bar. Similarly, the

recurrence relation does not work for i ¼ n: When the relationship is evaluated for i ¼ n; the value

Unþ1 is needed and it does not exist. The displacement equation evaluated for i ¼ n is

Unþ1 2 ð2 2 lÞUn þ Un21 ¼ 0

When the relationship is evaluated for ði ¼ 2; 3; …; n 2 1Þ, the resulting matrix problem is tridiagonal.

This is characteristic of the centered difference approximation. It is given by the following

matrix:

21 2 2 l 21 0 … 0

0 21 22 l 21 0 … 0

0 0 21 22 l 21 … 0

0 21 2 2 l 21

2

66666666666666664

3

77777777777777775

U1

U2

U3

Un

2

66666666666666664

3

77777777777777775

¼

0

0

0

0

2

66666666666666664

3

77777777777777775

This matrix problem does not have a unique solution since there are n unknowns and there are

only n 2 2 equations. The boundary conditions for the bar provide the needed two additional

constraints. For a bar that is fixed at both ends, the deflection is zero at the ends (x ¼ 0 and x ¼ L).

Under these conditions, U1 ¼ Un ¼ 0: If one of the ends is free then the stress is zero at that end,

which means that dU =dx ¼ 0: For the purposes of illustration, let us assume that the end at x ¼ 0 is

fixed and the other end is free. This means that U1 ¼ 0 and dU =dxln ¼ 0: If we use the centered

difference approximation for the derivative, we get

dU=dxln ¼ ðUnþ1 2 Un21Þ=2h ¼ 0

which means that Unþ1 ¼ Un21:

FIGURE 6.5 Discretized bar.

6-12 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

Unþ1 is the displacement of a point off the bar. We never compute Unþ1; however, as one can see

from a previous calculation, this value is used when the recurrence relationship is evaluated for i ¼ n;

Unþ1 2 ð2 2 lÞUn þ Un21 ¼ 0: Under the assumption that the stress is zero at the right-hand endpoint,

the remaining constraint becomes

ð2 2 lÞUn 2 2Un21 ¼ 0

For a bar that has one fixed left-hand endpoint and a free right-hand end, the matrix problem is the

following:

2 2 l 21 0 0 … 0

21 22 l 21 0 … 0

0 21 22 l 21 … 0

0 21 2 2 l 21

22 22 l

2

6666666666664

3

7777777777775

U2

U3

U4

Un21

Un

2

66666666666666664

3

77777777777777775

¼

0

0

0

0

0

2

66666666666666664

3

77777777777777775

This matrix problem has a unique solution because there are n 2 1 linearly independent equations and

there are n 2 1 unknowns.

6.4.1.1 Example

Consider a bar of unit length that is divided into four equal pieces. The left-hand end is fixed so that

U1 ¼ 0; and the right-hand end is free. The eigenvalue problem for this system is the following:

2 2 l 21 0 0

21 22 l 21 0

0 21 22 l 21

0 0 22 22 l

2

6666664

3

7777775

U2

U3

U4

U5

2

6666664

3

7777775

¼

0

0

0

0

2

6666664

3

7777775

The eigenvalues, l; are l ¼ ½0:1522; 1:2346; 2:7654; 3:8478􀀉: The natural frequencies can be recovered

from the eigenvalues since l ¼ h2v=c and c ¼

ffiffiffiffiffi

pE=r is determined by the material properties of the bar.

6.4.2 Beam

A second example of finite difference approximation is given by looking at a centered difference

approximation to compute the transverse vibration of a uniform beam. We will consider a special case of

the Euler equation for the beam. The following fourth-order equation results when the flexural rigidity,

ðEIÞ; is constant:

d4W =dx4 2 b4W ¼ 0

where b contains all of the material properties of the beam (i.e., b4 ¼ rAv2=EI).

The centered difference formula for the fourth derivative is given by

f iv ¼ 1=h4ðfiþ2 2 4fiþ1 þ 6fi 2 4fi21 þ fi22Þ

In order to see that this expansion is in fact second-order, combine the Taylor expansions for fiþ2; fiþ1;

fi; fi21; and fi22; and notice that you obtain an approximation for the fourth derivative which exactly

matches the combined Taylor approximations up to order h2: Using this approximation, the transverse

beam equation at the interior mesh point becomes

Wiþ2 2 4WIþ1 þ ð6 2 lÞWi 2 4Wi21 þ Wi22 ¼ 0

Numerical Techniques 6-13

© 2005 by Taylor & Francis Group, LLC

where l ¼ h4b4: This equation is valid for i ¼ 3; n 2 2: As before, the boundary conditions have to be

used in order to have sufficent conditions to determine all n values for the deflection, W : There will be

two additional constraints placed at each end. In this case, all of the types of boundary condition require

fictitious points. In other words, they require that the beam be extended beyond the physical boundaries

to include W21; W0; Wnþ1; and Wnþ2: Let us consider three types of boundary conditions: fixed end, free

end, and simply-supported end. For simplicity, we will derive all of the conditions for the right-hand

endpoint x ¼ L:

Fixed end: When the end of the beam is fixed, both the deflection W ¼ 0 and dW =dx ¼ 0: As above,

the centered difference approximation for the first derivative requires the introduction of a fictitious

point. These conditions translate into

Wn ¼ 0 and dW =dxln ¼ 1\h2ðWnþ1 2 Wn21Þ ¼ 0

These constraints reduce to the following:

Wnþ1 ¼ Wn21

Free end: The bending moment and the shear force are zero at the free end. These are given by a second

and third derivative, respectively. Assuming the free end is located at N; the two fictitious points

introduced are at n þ 1 and n þ 2: The discrete version of the two constraints is:

d2W =dx2ln ¼ 1=h2ðWn21 2 2Wn þ Wnþ1Þ ¼ 0

d3W =dx3ln ¼ 1=ð2h3ÞðWnþ2 2 2Wnþ1 þ 2Wn21 2 Wn22Þ ¼ 0

Simply-supported end: The boundary conditions of this type require that the deflection and the

moment are zero. When these are applied to the right-hand endpoint

Wn ¼ 0

d2W =dx2 ¼ 1=h2ðWnþ1 2 2Wn þ Wn21Þ ¼ 0

Combining these two equations, one sees that the value of W at the fictitious point is equal to the value

on the beam, i.e., Wnþ1 ¼ Wn21:

All three of these boundary condition types allow one to replace the fictitious values in the recurrence

relations with values on the beam. The procedures followed are very similar to those outlined for the

longitudinal vibration of a beam.

6.4.3 Summary of Finite Difference Methods for a Continuous System

1. Equation of motion for the bar: ›2u=›x2 ¼ 1=c2 ›2u=›t2:

2. Displacement equation for the bar: d2U=dx2 þ a2U ¼ 0:

3. Finite difference approximation: Uiþ1 2 ð2 2 lÞUi þ Ui21 ¼ 0:

4. Equation of motion for the beam: d4W =dx4 2 b4W ¼ 0:

5. Centered difference approximation for the fourth derivative: f iv ¼ 1=h4ð fiþ2 2 4fiþ1 þ 6fi 2

4fi21 þ fi22Þ:

6. Finite difference approximation for the beam equation: Wiþ2 2 4Wiþ1 þ ð6 2 lÞWi 2 4Wi21 þ

Wi22 ¼ 0: