Elementary Finite Element Method

.

Here p t q> 2
,

Eq. (2-11)

.

.

,

.

written as

is

u in

which a

x

=

and

1

into Eq. (2-10),

we


are

l

=

left

(p Q

.

=±* If the

(2-12)

t
approximate solution u

is

R(x)=Lu-f

(2-13)

i

which

is

zero

In the

if

u

=

u*.

method of weighted

residuals, the

aim

to find an approximate

is

solution u for w* such that the residual R{x) in Eq. (2-13) possible or

is

substituted

with a residual

is

made

as small as

minimized. In other words the error between the approximate

and the exact solution u is minimized. A number of schemes aim of minimization of R(x) details of some of the major schemes collocation, subdomain, least squares, and Galerkin methods are given in Appendix 1. Moreover, in the subsequent chapters, Galerkin's method is used to derive finite element equations for a number of simple problems. For the sake of completeness, only brief and general statements for these methods are given below. solution w*

are possible to achieve the

;

—

—

Mathematically the idea of minimization can be expressed as

f

R(x)W (x)dx t

= 0,

I

=

1, 2,

.

.

.

,

n,

(2-14)

where D denotes the domain of a structure or body under consideration. For the one-dimensional column problem, the domain is simply the linear extent of the column.

:

28

Steps in the Finite Element

Method

Chapter 2

W

In Eq. (2-14) the denote weighting functions. Various residual schemes such as collocation, subdomain, and Galerkin use different weighting t

For

functions.

instance, in the case of the collocation

method,

W=

1.

{

The

expression in Eq. (2-14) implies that the weighted value of R(x) over the

domain of a

structure vanishes. Figure 2-7 shows a schematic representation of Eq. (2-14). The shaded areas in Fig. 2-7(b) denote error between the

approximate and the exact solution u*

D

R(x) over

—

u over the domain D.

[Fig. 2-7(a)] is related to the error

D

integral of R(x) over

u*

—

w,

The residual and the sum or the

minimized.

is

u

i

R(x)= Lu-f

(a)

Figure 2-7 Schematic representation of residual, (a) Integration of

R

over D. (b) Error u*

As a simple

—

u over D.

illustration, let us consider the following

second-order

differ-

ential equation that governs the problems of one-dimensional stress deforma-

column and flow

tion in a

in

Chapters 3 and

d 2v dy

where v*

the

is

unknown

4, respectively

3

(2-1 5a)

f,

2

(deformation), y

is

the coordinate axis, c denotes a

material property, and /is the forcing function. For a column, c the applied load, area.

Assume

E is

that

the

= EA, /is

modulus of elasticity, and A is the cross-sectional 10 units of load; then Eq. (2-1 5a) 1 and

/=

EA =

specializes to

d 2 v* dy

An

(2-15b)

10.

:

approximate solution for v* can be written as a special case of

Eq. (2-11) as v

where


t

=

satisfy the

= ai +

cc 2

y

+

(x 3

y

2

= 2 *t9t,

(2

"

16)

2 This function should be chosen so as to 1,

Chapter 2

Steps in the Finite Element

Method

29

given by

R(y)

According to the

MWR, f*

=

W~ 1

I0

(2 " 17a)

'

Eq. (2-14) leads to

Riy)W {y)dy )

= 0,

i

=

(2-17b)

1, 2, 3,

Jo

or

^(^-10)w,(y)dy = or

j^d^y w

y

(y)dy

0,

(2-17C)

=

(2 . 17d)

Here L

is the length of the column which represents the domain D. Now, we can substitute for the second derivative in Eq. (2-17) by differentiating v in Eq. (2-16) twice. The final results will yield three simultaneous equations in

au

cc 2

,

and a 3

as

[[^^-^]w

[[^^y

2

(y)dy

[ -^-10]w d

which are solved for % u


substituted into Eq. (2-16),

At

this stage,

to proceed in the

it is

from Eq.

and a 3

,

we

.

3

y)dy

i(

=

0,

= 0,

(y)dy

When

(2-18)

=

0,

these values of

a„ a 2 and a 3 ,

are

obtain the approximate solution for v*.

not necessary to go into the details of the steps required (2-17) to Eq. (2-18); they are given in

Appendix

1

and

subsequent chapters.

ELEMENT EQUATIONS Use of either of the two foregoing methods will lead to equations ing the behavior of an element, which are commonly expressed as [k]{q]

= {Q},

describ-

(2-19)

=

=

vector of unknowns at the [k] element property matrix, {q} element nodes, and {Q} vector of element nodal forcing parameters. Equation (2-19) is expressed in a general sense; for the specific problem of vector of nodal displacements, stress analysis, [k] stiffness matrix, {q}

where

=

and {Q}

= = = vector of nodal forces. Details of the matrices in Eq. (2-19) will be

developed and described fully in subsequent chapters.

Step

Assemble Element Equations to Obtain Global

5.

or Assemblage Equations and Introduce Boundary

Conditions

Our

final

aim

to obtain equations for the entire

is

body

approximately the behavior of the entire body or structure. In

that define

fact, as will

discussed in various chapters, use of the variational or residual procedure

body

relevant to the entire

;

it is

for simplicity that

we view

be is

the procedure in

Step 4 as having been applied to a single element.

Once element,

the element equations, Eq. (2-19), are established for a generic

we

are ready to generate equations recursively for other elements by

using Eq. (2-19) again and again. Then equations. This assembling process continuity (Chapter

1). It

is

we add them

together to find global

based on the law of compatibility or

requires that the

body remain continuous; that

is,

the neighboring points should remain in the neighborhood of each other after

the load

is

applied (Fig. 2-4). In other words, the displacements of two

adjacent or consecutive points must have identical values [Fig. 2-8(a)].

Depending on the type and nature of the problem, we may need to enforce the continuity conditions more severely. For instance, for deformations occurring in a plane,

it

may be

sufficient to enforce continuity

of the dis-

Figure 2-8 Interelement compatibility, (a) Compatibility for plane

problems, (b) Compatibility for bending-type problems. Slopes

may not

be

equal

Elements

Equal displacements (a)

Slope

Equal slopes or gradients

Equal displacements (b)

30

Chapter 2

Steps

On

placements only.

in the Finite

Element Method

31

the other hand, for bending problems, the physical

properties of the deformed

body under the load

the continuity of displacements

we

requires that in addition to

ensure that the slopes or the

first

deriva-

of displacements are also continuous or compatible at adjacent nodes [Fig. 2-8(b)]. Often it may become necessary to satisfy compatibility of the curvatures or the second derivative also.

tive

Finally,

we

obtain the assemblage equations, which are expressed in

matrix notation as [Ktfr}

= {R},

(2-20)

where [K] = assemblage property matrix, {r} = assemblage vector of nodal unknowns, and {Q} = assemblage vector of nodal forcing parameters. For stress-deformation problems, these quantities are the assemblage stiffness matrix, nodal displacement vector, and nodal load vector, respectively.

BOUNDARY CONDITIONS Until now we have considered Equation (2-20)

[Fig. 2-9(a)].

only the properties of a body or structure

tells

us about the capabilities of the body to

is an engineer. How he depend on the surroundings and the problems he faces; these aspects can be called constraints. In the case of engineering bodies, the surroundings or the constraints are the boundary conditions. Only when we introduce these conditions can we decide how the body will perform. Boundary conditions, then, are the physical constraints or supports [Fig. 2-9(b)] that must exist so that the structure or body can stand in space uniquely. These conditions are commonly specified in terms of known values of the unknowns on a part of the surface or boundary S and/or gradients or derivatives of the unknowns on S 2 Figure 2- 10(a) depicts a beam. In the case of the simply supported beam, the boundary S is the two end points where

withstand applied forces. will

perform

It is just like

saying that one

his engineering duties will

x

.

t

the displacements are given. This type of constraint expressed in terms of

displacements

is

often called the essential, forced, or geometric boundary

conditions.

At

the simple supports of the beam, the bending

moment

zero

that

is,

the second derivative of displacement vanishes. This type of constraint

is

is

;

often called a natural boundary condition.

Figure 2- 10(b) shows a cylinder through which fluid or temperature On the boundary S^ temperature or fluid head is known; this is the essential boundary condition. The right end is impervious to water, or

flows.

insulated against heat; then the

heat flux, which

is

boundary condition

proportional to the

first

is

specified as fluid or

derivative of fluid

head or tem-

perature. This is the natural boundary condition. To reflect the boundary conditions in the finite element approximation of the body represented by Eq. (2-20), it is usually necessary to modify these

[K]{r}= {R}

(a)

[K](7}=

{R]

Constraints

(b)

Figure 2-9 Boundary conditions or constraints, (a) constraints (b)

Body with

Figure 2-10 Examples

boundary conditions,

Body without

constraints.

of boundary conditions, (b) Pipe flow with

(a)

s,

"TT (a)

$

\ Flow (b)

Beam

with

boundary conditions.

S,

Chapter 2

Steps

in the Finite

Element Method

33

equations only for the geometric boundary conditions. Further details and procedure for such modification are given in subsequent chapters. The final

modified assemblage equations are expressed by inserting overbars as [K]{f}

Step

6.

Solve for the Primary

Equation (2-21)

is

a

set

= {R}.

(2-21)

Unknowns

of linear (or nonlinear) simultaneous algebraic

equations, which can be written in standard familiar form as

K K

ll

rl

2l r t

+ K 12 r + + K22 r +

•

2

•

2

•

•

+K +K

ln r n

2n rn

= Ru = Ri> (2-22)

K s + Km2 r + n

2

x

•

• •

+K

nn rn

=R

r

These equations can be solved by using the well-known Gaussian elimination or iterative methods. Detailed coverage of these methods is beyond the scope of this book; however, later when individual problems are considered, some of these methods will be illustrated. Moreover, brief descriptions of these in Appendix 2. At the end of this step, we have unknowns (displacements) r l9 r2 , rn These are called primary unknowns because they appear as the first quantities sought in the basic Eq. (2-22). The designation of the word primary will change depending on the unknown quantity that appears in Eq. (2-22). For instance, if the

methods are given

solution

solved for the

problem

is

,

formulated by using stresses as unknowns, the stresses

called the primary quantities.

be the

fluid or velocity

Step

7.

.

will

be

For the flow problem the primary quantity can

head or potential.

Solve for Derived or Secondary Quantities

Very often additional or secondary quantities must be computed from the primary quantities. In the case of stress-deformation problems such quantities can be

strains, stresses,

moments, and shear forces for the flow problem they ;

can be velocities and discharges. It is relatively straightforward to find the secondary quantities once the primary quantities are known, since we can make use of the relations between the strain and displacement and stress

and

strain that are defined in Step 3.

Step

8.

Interpretation of Results

final and the important aim is to reduce the results from the use of the element procedure to a form that can be readily used for analysis and design. The results are usually obtained in the form of printed output from

The

finite

the computer.

We then

select critical sections

of the body and plot the values

Steps

34

in the Finite

Element Method

of displacement and stresses along them, or often very convenient

and

less

we can

Chapter 2

tabulate the results. It

is

time consuming to use (available) routines and

ask the computer to plot or tabulate the results.

REFERENCES [1]

Turner, M.

J.,

Clough, R. W., Martin, H. C, and Topp, L. C, "Stiffness and Complex Structures," /. Aero. Sci., Vol. 23, No. 9,

Deflection Analysis of Sept. 1956. [2]

Clough, R. W., "The Finite Element Method in Plane Stress Analysis," Proc. 2nd Conf. on Electronic Computation, ASCE, Pittsburgh, Sept. 1960.

[3]

Crandall,

[4]

Argyris,

S.

J.

H., Engineering Analysis, McGraw-Hill,

New

in

York, 1956.

H., Energy Theorems and Structural Analysis, Butterworth's,

London, 1960. [5]

[6]

Desai, C. S., and Abel, J. F., Introduction Nostrand Reinhold, New York, 1972.

Element Method,

Vah

Pian, T. H.*H„ and Tong, P., "Finite Element Methods in Continuum .

Mechanic^" York, 1972. [7]

to the Finite

Finlayson,

Academic

in

B

.

Advances

A

Press,

.

,

in

Applied Mechanics, Vol.

12,

Academic

Press,

New

The Method of Weighted Residuals and Variational Principles, York, 1972.

New

ONE-DIMENSIONAL STRESS DEFORMATION

203SQ01 INTRODUCTION From

we shall consider engineering problems idealThe main motive for treating these simple problems

here through Chapter 10

ized as one-dimensional.

to introduce the reader to the details of various steps so that basic concepts can be understood and assimilated thoroughly without undue complex and is

lengthy derivations. An advantage to this approach is that hand calculations can be performed for two- and three-dimensional problems this can become ;

increasingly difficult.

Although simple problems are treated (Chapters 3-10), we introduce concepts and terms that are general and relevant to advanced theory and applications. These concepts are explained and defined in easy terms, often with intuitive and physical explanations. It may be mentioned that

many

one-dimensional idealizations permit not only simple derivations but often provide satisfactory solutions for

As our

first

many

practical problems.

problem we consider the case of a column,

strut,

or bar of

uniform cross section subjected to purely axial loading [Fig. 3- 1(a)]. Under these conditions we can assume that the deformations will occur only in one, vertical direction. Consequently, we can further assume that the column can be replaced by a [Fig. 3- 1(b)].

line

Now we

with the axial

stiffness

EA lumped

consider derivations of the

the step-by-step procedure as described in Chapter

finite

at the centerline

element method in

2.

35

One- Dimensional Stress Deformation

36

Chapter 3

ic

s Nodes

Finite

elements

S EA zMx^

777Z

(0

(b)

(a)

Figure 3-1 Axially loaded column, (a) Actual column, (b) One-

dimensional idealization,

Step

1.

(c)

Discretization.

Discretization and Choice of Element Configuration

Before we proceed further it is necessary to describe the coordinates or geometry of the column by using a convenient coordinate system. In the onedimensional approximation, it is necessary to use only one coordinate along the vertical direction. is

We

y axis. Because this coordinate system column (or structure), it can be called the

call this the

defined to describe the entire

global coordinate system.

We now

discretize the

column

into

an arbitrary number of smaller units

that are called finite elements [Fig. 3-l(c)].

The

intersections of elements are

called nodes or nodal points.

At

this stage, before Step 2,

it is

useful to introduce the concept of a local

or element coordinate system. There are a

number of advantages

local system for deriving element equations (Steps 2-4)

for multidimensional problems,

makes the required

tions extremely simple to handle. Indeed,

it is

;

its

derivatives

and integra-

possible to obtain

by using the global system, but use of the local system subsequently) facilitates the derivatives and is economical.

tions

to using a

use, particularly

(as

all

we

deriva-

shall see

EXPLANATION OF GLOBAL AND LOCAL COORDINATES For a simple explanation, let us consider an example, Fig. 3-2(a); here we need to define or survey a plot of land and relate it to a standard or global point or benchmark A. Surveying the plot could entail locating each point and establishing its distance from point A. Let us assume that point A away from the plot and that it is difficult to establish a direct relation with A. A local point B is available, and its distance from A is known, howin is

it

far

Chapter 3

One- Dimensional Stress Deformation

37

Plot of land

"^S B

Local station

"Global" station

(a)

1

r > 1

®i

>s =

©• .L-+1-

1

A

A(Dfo

j

7

-1 s

Y2

=

y2 Y3 Vi

Vi I

Global y

v/^W s

Global y

= Local coordinate = y//

^^

77777,

L =

(y-y 3 )/(//2)

(b)

(c)

Figure 3-2 Global and local coordinates, (a) Concept of global and local coordinate systems, (b) Local coordinate

point

ever. Point

1. (c)

B is

measured from node

Local coordinate measured from midnode

accessible

from

all

points in the plot.

distances of each point in the plot

3.

We

can

first

define the

from the point B; then, knowing the

distance x AB from point A to point B, it is possible to define the distances of the points in the plot with respect to point A. For example, if the distance

of any point P from

B is x B

,

then

*a

its

distance from point

= x* +

*ab-

A is (3-la)

Here we can call the measurements with respect to point B local coordiand those from point A, global coordinates. We can see that the mea-

nates,

One- Dimensional Stress Deformation

38

surements with respect to point B, which of the plot

much

is

in the vicinity,

Chapter 3

make

the definition

simpler than those with respect to point A, which

is

more

There can be a number of possible points like B with respect to which we can define the local coordinate system such choices will depend on the nature of the problem and the convenience and ease of difficult to reach.

;

measurements.

The basic idea of the use of local coordinate system(s) ment method is very similar to the foregoing concept.

in the finite ele-

LOCAL AND GLOBAL COORDINATE SYSTEM FOR THE ONE-DIMENSIONAL PROBLEM As mentioned above,

there are a

number of ways

in

which we can define local

coordinate systems for the one-dimensional problem (Fig. here two local coordinate systems, Figs. 3-2(b) and In the

first

case,

we measure

3-1).

We

consider

(c).

from the node point 1 is measured from analogous to point B and that

the local coordinate

of a generic element e [Fig. 3-2(b)]. The global coordinate the base of the column. the base

is

Note

that

A

analogous to point

node

1

is

in Fig. 3-2(a).

y then the global coordinate of any point ;

y=y+y

t

We

call the local

in the element

coordinate

is

(3-lb)

;

hence

=y- yi-

y Often

it is

0- lc)

more convenient to express the local coordinate as a nondi-

mensional number; such a procedure can considerably facilitate the integrations and differentiations involved in the subsequent computation. Here we nondimensionalize by dividing^ by the length of the element; thus o

y -y\ _ yi-yx

_

y

(3-ld)

.

i

where s = nondimensionalized local coordinate, / = length of the element, and Vi and y 2 = global coordinates of nodes 1 and 2, respectively. Note that because we nondimensionalized the coordinate by dividing by the length of the element, I = y 2 — J>i, the value of s varies from zero at node 1 to unity at node 2. In the second alternative, we can attach the origin of the local system at an intermediate point in the element, say, at the midpoint [Fig. 3-2(c)]. Here the local coordinate

is

written as

L The values of L range from

—

1

=

^-

at point

1

to

(3-2)

at point 3 to

1

at point 2.

One- Dimensional Stress Deformation

Chapter 3

An

important property of these local coordinates

39

is

that they are in

nondimensionalized form, and their values are expressed as numbers and often lie

between zero and unity.

property that imparts simplicity to the

It is this

subsequent derivations. Step

2.

Select Approximation for the

Unknown

One of the main grasp at this stage

Model

or Function

(Displacement)

ideas in the finite element is

method

that the reader should

the a priori or in advance selection of mathematical

functions to represent the deformed shape of an element under loading. This implies that since

it is

a closed form or exact solution,

difficult to find

we

guess a solution shape or distribution of displacement by using an appropriate this function, we must follow the laws, and constraints or boundary conditions inherent in the problem. The most common functions used are polynomials. In the initial stages of the finite element method, the polynomials used were expressed in terms of generalized coordinates however, now most finite element work is done by using interpolation functions, which can often be considered as transformed generalized coordinate functions.

mathematical function. In choosing principles,

;

Generalized Coordinates

The

we can use

simplest polynomial that

is

the one that gives linear

variation of displacements within the element [Fig. 3-3(b)],

v

=a +

oc 2

1

(3-3a)

y,

or in matrix notation. l

[1

y]\

\

<;)

(3-3b)

or

W= where o^ and a 2 in the element,

y and y2 To show x

- 3c

)

= generalized coordinates, y = the coordinate of any point and v = displacement at any point in the element. The a's

contain the displacements at nodes v .

(3

[+]{«}.

this,

we

first

x

and v 2 and the coordinates of nodes 1 and 2 by substituting

evaluate v for points

fory, Vl

=«,

+a y

v2

=

+

<%!

7

l9

a 2j 2

,

(3-4a)

(3-4b)

or in matrix notation,

N=P "IN

(3-40

= [A]{«}.

(3-4d)

or {q}

One- Dimensional Stress Deformation

40

©

H

v

©

.

©

y3

Yi

Y2

h

'1

^

V

Chapter 3

(a)

V =

OLy

+

OL

2

Y

(b)

N, =/ 2 //=1-sorN 1 = 1(1

-

L)

Equal

© A ©© A © (e)

Figure 3-3 Linear interpolation functions and interelement compatibility.

Here

r

{q}

=

[v x

v2]

is

the vector of nodal displacements. Second,

we

solve

for {a} as

{a}

= [A]-i{q},

where

(3-5a)

[A]-

-yi~

yi

1

_-i

l/l

~ 1

yi

-

(3-5b)

i

yi

-y\ _-i

-yC l

(3-5c)

.

Chapter 3

Here

One- Dimensional Stress Deformation

= Jacobian

|/|

element.

=y —y =

determinant

2

I

l

41

equals the length of the

Thus

\[y

fai]

/

x

(2

-y

2

(3-5d) 1

1

x

(2

1)

2)

(2

v2 \ x 1)

Therefore

(3-5e)

which shows that

a's are functions of

and made up of y u y 2i v l9 and v 2 Note .

that a's are related to but are not explicit functions of nodal displacements.

This

is

one of the reasons we call them generalized coordinates. can substitute {a} into Eq. (3-3) to express v in terms of the

Now we

nodal displacements: v

=

y 2v

i

-y,v 2

+

-»«+*> (

(3-6a)

,,

)

-(Z ^)*.-(ZL f 2 )^ =

fv,

=

tf,*,

l

+ fv = l

2

(l

-

y.) Vl

(3

+

j-v 2

+Nv = — s, N = y/l = s,

6b >

(3-6c)

(3-6d)

2i

2

"

where iVj = 1 — y/l 1 and y === /, [Figs. 3-2(b) and 2 3-3]. Equation (3-6d) leads us to the concept of interpolation function models. In this equation, N and N2 are called interpolation, shape, or fozs/s functions. The displacement v at any point in the element can now be expressed as x

v=[N,

N ]H

(3-7a)

2

= [N]{q}, where [N]

is

(3-7b)

called the matrix of interpolation, shape, or basis functions.

property of the interpolation functions

example, in the above,

N +N = t

2

is

that their

sum

A

equals unity. For

1

INTERPOLATION FUNCTIONS Since our aim in the finite element analysis

is

to find nodal displacements

v x and v 2 we can see the advantage of using approximation models of the type in Eq. (3-7). In contrast to Eq. (3-3a), here the displacement v is expressed directly in terms of nodal displacements. Also, the use of interpola,

42

One- Dimensional Stress Deformation

tion functions

makes

it

Chapter 3

quite easy to perform the differentiations

and

integra-

tions required in the finite element formulations.

In simple words, an interpolation function

of unity at the node point to which

it

is

pertains

a function that bears a value

and a value of zero

nodes. For example, Figs. 3-3(c) and (d) show distributions of

along the element; the function point

JV,

pertains to point

1,

and

N

2

N

x

at other

N

and

2

pertains to

2.

For a given element it is possible to devise and use different types of local coordinates and interpolation functions N Let us consider the second coordinate system [Fig. 3-2(c)] described in Step 2 and in Eq. (3-2), and t

.

express the v as

v=i(l-L) Vl +i(l+L)v

= N v, + N v 2

x

(3-8a)

2

(3-8b)

2

= [tfi N ]l Vl \

(3-8c)

=

(3-8d)

2

Note

that

N

t

(i

=

1,

[N]{q].

2) in

Eq. (3-8b) are different from

N

t

(i

=

1,

2) in Eq.

same linear variation within the element. Alternative coordinate systems based on local measurement from other points on the element are possible; for example, see Prob. 3-1. (3-6d); however, they yield the

RELATION BETWEEN LOCAL AND GLOBAL COORDINATES

An

important point to consider is that a one-to-one correspondence exists between the local coordinate s or L and the global coordinate y of a point in the element. For example, for L in Eq. (3-2) we have

y

= i(l - L)y + ^(1 + L)y = #iJ>i + N y = [N]{y„], t

2

where

T {y„]

=

[y

t

y2 ]

is

(3-9a)

2

(3-9b)

2

(3-9c)

the vector of nodal coordinates.

In fact an explanation of the concept of isoparametric elements, which

is

can be given at this stage. A comparison of Eq. (3-8) and (3-9) shows that both the displacement v and the coordinate y at a point in the element are expressed by using the same (iso) the most

common

procedure

interpolation functions.

An

now

in use,

element formulation where

similar) functions for describing the deformations in

(or geometry) of an element

This

is

is

we use the same (or and the coordinates

called the isoparametric element concept

rather an elementary example;

we

and more general isoparametric elements.

shall subsequently

[1].

look at other

VARIATION OF ELEMENT PROPERTIES We often tacitly assume that the material properties such as cross-sectional area A and the elastic modulus E are constant within the element. It is not necessary to assume that they are constant. We can introduce required variation, linear or higher order, for these quantities.

For instance, they can

be expressed as linear functions:

E = N,E + N E = [NtfEJ, A = N A + N A = [N]{A„}, E and {A„} r = [A A are the x

T

where [E n } of E and A

= at

[E

x

X

2

2

(3-10)

t

2

2

(3-11)

2]

nodes

1

2]

x

and

vectors of nodal values

2, respectively.

REQUIREMENTS FOR APPROXIMATION FUNCTIONS As we have

stated before, the choice of an approximation function is guided by laws and principles governing a given problem. Thus an approximation function should satisfy certain requirements in order to be acceptable. For general use these requirements are expressed in mathematical language. However, in this introductory treatment, we shall discuss them in rather

simple words.

An

approximation function should be continuous within an element. The

linear function for v [Eqs. (3-7)

and (3-8)]

is

indeed continuous. In other words,

it

does not yield a discontinuous value of v but rather a smooth variation of

v,

and the variation does not involve openings, overlaps, or jumps.

up

The approximation function should provide interelement compatibility by the problem. For instance, for the column prob-

to a degree required

lem involving

axial deformations,

it

is

necessary to ensure interelement

compatibility at least for displacements of adjacent nodes. That

is,

the

approximation function should be such that the nodal displacements between adjacent nodes are the same. This is shown in Fig. 3-3(e). Note that for this case, the higher derivatives such as the first derivatives may not be compatible.

The displacement

at node 2 of element 1 should be equal to the displacement node 1 of element 2. For the case of the one-dimensional element, the linear approximation function satisfies this condition automatically. As indicated in Chapter 2 (Fig. 2-8), satisfaction of displacement compatibility by the linear function does not necessarily fulfill compatibility of first derivative of displacement, that is, slope. For axial deformations, however, if we provide for the compatibility up to only the displacement, we can still expect to obtain reliable and convergent solutions. Often, this condition is at

tied in with the highest order of derivative in the energy function such as the

potential energy. derivative dv/dy

For example,

= €,

is

1

;

in

Eq. (3-21) below, the highest order of

hence, the interelement compatibility should 43

:

One- Dimensional Stress Deformation

44

include order of v at least

up

to

(zero), that

is,

Chapter 3

displacement

v.

In general,

should provide interelement compatibility up to order

the formulation

—

n 1, where n is the highest order of derivative in the energy function. Approximation functions that satisfy the condition of compatibility can be called conformable.

The

other and important requirement

is

that the approximation function

should be complete; fulfillment of this requirement will assure monotonic convergence. Monotonic convergence can be explained in simple terms as a process in which the successive approximate solutions approach the exact solution consistently without changing sign or direction.

For

instance, in

approximate areas approach the exact area in such a way that each successive value of the area is smaller or greater than the previous value of area for upper and lower bound solutions, respectively. Completeness can be defined in a number of ways. One of the ways is to Fig. l-8(c) the

relate

it

to the characteristics of the chosen approximation function. If the

function for displacement approximation allows for rigid body displacements

(motions) and constant states of strains (gradients), then the function can be

considered to be complete.

mode For

A

rigid

body motion represents a displacement

that the element can experience without development of stresses in

it.

instance, consider the general polynomial for v as

v

=a + 2

cc 2

y

In Eqs. (3-3a) and (3-6)

|

+

a 3y 2

+

a 4 j> 3

we have chosen

the general polynomial as

+

•

•

•

+

n cc n+

(3-3d)

iy

a linear polynomial by truncating

shown by the

vertical

dashed

line.

The

linear

approximation contains the constant term a, which allows for the rigid body displacement mode. In other words, during this mode, the element remains

and does not experience any strain or stress, that is, a 2 y = 0. The requirement of constant state of strain (e y) for the one-dimensional column deformation is fulfilled by the linear model, Eq. (3-3a) because of the existence of term a 2 y. This condition implies that as the mesh is refined, that in each element is, the elements become smaller and smaller, the strain e y

rigid

approaches a constant value. In the case of one-dimensional plane deformations in the column, the

condition of constant state includes only e y

may

— the

first

derivative or gradient

and more general problems such as beam and plate bending. In such cases, it will be necessary to satisfy the constant strain state requirement for all such generalized strain or gradients of the unknowns involved, e.g., see Chapters 7, 11-14. In addition to monotonic convergence, we may be interested in the rate of

v.

Additional constant strain states

exist in other

of convergence. This aspect is often tied in with the completeness of the polynomial expansion used for the problem. For instance, for the onedimensional column problem, completeness of the approximation function requires that a linear function, that

is,

a polynomial of order («) equal to one

Chapter 3

is

One- Dimensional Stress Deformation

needed. Completeness of the polynomial expansion requires that

including and up to the satisfied for the linear

from Eq.

first

order should be included. This

model, Eq. (3-3a) since

(3-3d). In the case of

cubic approximation function is

45

all

terms up to n

beam bending, Chapter is

6,

to

and includes the order

It

may happen

not include satisfy the

all

we

all

it

terms

automatically

=

1

are chosen

shall see that a

Then

required to satisfy completeness.

necessary to choose a polynomial expansion such that

up

is

includes

all

it

terms

3.

that an approximation function of a certain order (n)

may

terms from the polynomial expansion, Eq. (3-3d), and

still

requirements of rigid body motion and constant states of

As an example of a complete approximation

satisfying rigid

strain.

body motion and

constant state of strain in a two-dimensional problem, but not complete

of the polynomial expansion, see Chapter 12, Eq. (12-9). For two-dimensional problems, the requirement of completeness of polynomial expansion can be explained through the polynomial expansion represented by using Pascal's triangle, see Chaps. 11-13. Here we have given rather an elementary explanation of the requirements for approximation functions. The subject is wide in scope and the reader interested in advanced analysis of finite element method will encounter the subject quite often. For example, the completeness requirement can be further explained by using the so-called "patch test" developed by Irons [2]; it is discussed in refs. [3] and [4]. Moreover, the approximation model should satisfy the requirements of isotropy or geometric invariance [5]. These topics are considered beyond the scope of this text. in the sense

Step

3.

Define Strain-Displacement and Stress-Strain Relations

For stress-deformation problems, the actions or causes (Chapter 2) are and the effects or responses become strains, deformations, and stresses. The basic parameter is the strain or rate of change of deformation. The link connecting the action and response is the stress-strain or constitutive law of the material. It is necessary to define relations between strains and displaceforces,

ments and stresses and strains for the derivation of element equations in Step 4. Hence, in this step we consider these two relations. We note at this stage that although we use familiar laws from strength of materials and elasticity for the stress-deformation

problem,

the relations relevant to specific topics.

flow problem (Chapter

4),

the relation

in later chapters

we

shall use

For instance, in the case of the fluid between gradient and fluid head and

Darcy's law will be used.

Returning to the axial deformation of the column element, the displacement relation, assuming small strains, can be expressed as ey

= ^> ay

strain-

(3-12a)

One- Dimensional Stress Deformation

46

Chapter 3

where e y = axial strain. Since we have chosen to use the local coordinate L and since our aim is to find dv/dy in the global system, we can use the chain rule of differentiation as

€

Now, from Eq.

(3-2)

dv

,~ 101 v (3 " 12b)

= d_(l^M)=\ dy\

dy

dL

dL dv -ayTL'

we have d_L

and from Eq.

dv = Ty =

>

(3-12c)

1/2

(3-8)

= ^[±(1 -

L)v t

+ i(l +

=

L)v 2 ]

1]N.

±[-1

(3-12d)

Substitution of Eqs. (3-12b)-(3-12d) into Eq. (3-12a) leads to

= -}-[-

e,

1

1]H

(3-13a)

or

=

{€,} (1

where

[B]

=

(1//)

x

x [— 1

1)

x

{q} 2) (2

1)

is

a one-dimensional problem, the strain vector

contains only one term, and the matrix [B]

{e y}

(3-13b)

,

x

can be called the strain-displacement transforma-

1]

tion matrix. Because this

[B] (1

is

only a row vector.

retain this terminology for multidimensional problems

;

We

shall

however, with multi-

dimensional problems these matrices will have higher orders.

The student can

easily see that Eq. (3-1 2d) indicates constant value

strain within the element; this

the displacement.

is

because

we have chosen

of

linear variation for

We can then call this element a constant-strain-line element.

STRESS-STRAIN RELATION For

simplicity,

we assume

elastic (Fig. 3-4).

that the material of the

column element

is

linearly

This assumption permits use of the well-known Hooke's

law,

ay

— Ey € y

(3- 14a)

,

or in matrix notation, {«,} (1

where

[C]

is

X

1)

=

[C] (1

X

{€,}

1)(1

X

,

(3-14b)

1)

the stress-strain matrix. Here, for the one-dimensional case,

matrices in Eq. (3- 14b) consist of simply one scalar term.

One- Dimensional Stress Deformation

Chapter 3

47

i\

A^ = U cl u

-

III

Figure 3-4 Linear elastic constitutive or stress-strain (Hooke's) law.

Substitution of Eq. (3-13) into Eq. (3-14b)

now

allows us to express {a y}

in terms of {q} as

{«y } Step

A

4.

(3-15)

Derive Element Equations

number of procedures

Among

= [C][B]{qJ.

are available for deriving element equations.

and residual methods. Principles based on and complementary energies and hybrid and mixed methods are used within the framework of variational methods. As described in Chapter 2 and in Appendix 1, a number of schemes such as Galerkin, collocation, and least squares fall under the category of residual methods. We shall use some of these methods in this chapter and subsequently in other chapters. these are the variational

potential

PRINCIPLE OF In simple words,

MINIMUM POTENTIAL ENERGY if

a loaded elastic body

is

in equilibrium

under given geo-

metric constraints or boundary conditions, the potential energy of the

deformed body assumes a stationary value. In the case of linear elastic bodies is a minimum; since most problems we consider involve this specialization, for convenience we shall use the term minimum. Figure 3-5 shows a simple axial member represented by a linear spring in equilibrium, the value

with spring constant k(F/L). Under a load P, the spring experiences a displacement equal to v. The potential or the potential energy n p of the spring is composed of two load (see Chapter 2) parts, strain energy U and potential p of the external

n,

The

strain energy

U

W u+ W

:

B.

(3-16)

can be interpreted as the area under the stress-strain when we minimize II,, we differentiate or

curve (Fig. 3-4). Mathematically,

:

One- Dimensional Stress Deformation

48

Chapter 3

IP

I

n

w^% Figure 3-5 Idealized linear spring.

take variations of Tl p with respect to the displacement v. While doing this we assume that the force remains constant, and we can relate variation of work

W by the load and the potential of the load as

done

SW

6Wr

(3-17)

where 8 denotes arbitrary change, variation, or perturbation. For our purpose

we can consider

it

to imply a series of partial differentiations.

sign in Eq. (3-17) occurs because the potential

into is

work by

these loads,

W

p

The negative

of external loads

is

lost

W. Then the principle of minimum potential energy

expressed as

sn B = su + 8wD =

su-sw =

(3-18)

o.

There are two ways that we can determine the minimum of II,: manual and mathematical. Both involve essentially an examination of the function represented by lip until we find a minimum point. For simple understanding, we first consider the manual procedure and write the potential energy for the spring (Fig. 3-5), assuming undeformed state of the spring as the datum for potentials, as

n

-Pv Ikv - Pv, 2

(3-19a)

where kv = force in the spring and %(kv)v denotes strain energy as the area under the load-displacement curve (Fig. 3-5). Since the load in the spring goes from to kv, we have to use average strain energy. The term Pv denotes the potential of load P; since we have assumed P to be constant, this term does not include £. We further assume that P = 10 units and k = 10 units per unit deformation. Then

U p = \\0v 2

Now we

search for the

;

a positive v

is

=

5v 2

-

minimum by examining

various values of deformations table

lOv

assumed

v.

The

(3- 19b)

10v.

values of the potential Il p for

results are

shown

in the following

to act in the direction of the applied load

One- Dimensional Stress Deformation

Chapter 3

V

n,

-2.000 -1.000

+40.000

49

+ 15.000

0.000

0.000

0.125

-1.1719 -2.1875 -3.7500 -5.0000

0.250 0.500 1.000

2.000

0.0000

3.000

15.0000

4.000

40.0000

5.000

75.0000

etc.

Figure 3-6 shows a plot of lip versus value at v will

=

1.

deform by

On

Hence, under 1

known

is

can be seen that lip has a

when

minimum

in equilibrium,

we can perform the procedure of "going to and minimum by using mathematics. It is

11, to find its

that a function assumes a

derivative

v. It

10 units, the spring,

unit.

the other hand,

on the function

P=

zero.

Applying

minimum

\kvbv

Pdv

well

value at a point(s) where

this principle to II, in

SU,

fro"

Eq.

(3- 19a),

we

its

obtain

=0

or (kv

- P)Sv = 0.

(3-20a)

Figure 3-6 Variation of potential energy.

/

i

40

Is** 30 20

\

\ \ 10

\

\ -3

-1

^^" '2

3

Minimum

point

1

-10 -20 -30

4

5

v

One- Dimensional Stress Deformation

50

Since Sv

is

arbitrary, the term in parentheses

Chapter 3

must vanish. Therefore,

kv-P = or

=P,

kv

(3-20b)

which is the equation of equilibrium for the spring. STl p = in Eq. (3-20) is analogous to equating d\\ p \dv = 0, which will result in the same equilibrium equation (3-20b). Substitution of the numerical values gives \0v

Therefore, v

We method

=

1.0 unit, the

=

10.

same answer

as before.

note here that in most problems solved by using the

Hp

is

a function of a large

ments. Consequently,

it

is

number of parameters

most economical and

matical methods because the manual procedure

finite

element

or nodal displace-

direct to use the is

mathe-

cumbersome and

often

impossible.

The mathematical procedure involved minimization of lip. For simplicity, we may view the process as simply taking derivatives of Tl p In general, how.

ever, the minimization will involve calculus of variations. In

treatment in this book

we

Now we return to the column element (Fig.

n> =

most of the

shall use the simple differentiation concept.

3-7)

and write FL

as

[5]

M

c e dv iff ^ ' >

~

iff

fvdv

~

if

T'vds

S P„v„

Figure 3-7 Generic column element with loads.

jp„

© U\®

©n

% \

^

(3-2 la)

One- Dimensional Stress Deformation

Chapter 3

=

where \o y e y

strain energy per unit

Ty =

(weight) per unit volume,

S

area,

= part

x

M — number

maximum

A

value of

vt

level,

M=

Y = body

=

Pa =

applied

displacement corresponding to

the joint forces

P a is now appropriate. Pu can be treated as

Pih

Since a

applied to the total structure,

is

force

2.

contribution of the joint force

we

volume,

of points at which joint forces are applied; here the

comment concerning

joint force

=

surface loading or traction per unit surface

of the surface on which surface loading acts,

nodal (joint) forces at local

and

V

volume,

51

Pig

applied at point

i

a local

in the global sense. Later

more convenient to add contributions of concentrated when we consider potential energy of the entire body; becomes clearer when the total structure is considered.

shall see that

it is

joint or nodal forces their relevance

The terms

Eq. (3-21 a) are essentially similar to those in Eq.

in

now we

except that

assign

volume

(3- 19a),

to the element instead of treating

it

as a

spring.

For the present, we assume that the cross-sectional area A of the element is

constant; then Eq. (3-2 la) reduces to

- A\

Tyvdy -

n,

=

^

Ty

is

the (surface) loading per unit length along the centerline of the

Here

idealized line

o y e y dy

column

Yvdy-

£ Pav

t

(3-21b)

.

(Fig. 3-7).

Equation (3-2 lb) can

now

be expressed

in

terms of the local coordinate

system by using the transformation of Eq. (3-2) as

dy

=

(3-22)

-L-dL.

Therefore, 1

11,

=

^

Next we

ay€ y

dL-^C

YvdL-l-j TyvdL-tP»Vi'

(«lc)

J"'

and a y from Eqs.

substitute for v, €y9

(3-8), (3-13),

and

(3-14),

respectively, in Eq. (3-21c) to obtain, in matrix notation,

n* = T

j\ _1

F

(1

x

1)(1

-yj "'

^WL "t£

[C]

{€'

x

1)(1

[N] (1

x

x

2) (2

x

(1

1)

dL-^

T,

{q} 1)(1

x

[N]

1)

x 2)(2 x

Pu (1

x

f dL

{q} 1)(1

x

1)

(3-21d)

v,

l)(l

x

1)

One- Dimensional Stress Deformation

52

Chapter 3

or

=^

n

{qf

[W

x

x

(1

Al

2) (2

x

-yj

[B]

x 1)(1 x

1)(1

[N] (1

[C]

Y

{q}

x

2) (2

[N]

x

1)(1

x

2) (2

1)

M r

x

x

2) (2

x

1)(1

(3-2 le)

T

1

(1

1)

dL

fy dL-

{q}

dL

{q}

(1

1)

x

1)(1

x

1)

where Y and Ty are assumed to be uniform. Equation (3-21) represents a quadratic function expressed in terms of v and v 2 In matrix notation, transposing in Eq. (3-2 Id) is necessary to make .

x

the matrix multiplication in {e } r [C]{€ >,} consistent so as to yield the scalar 3;

a y e y = Eel m Eq. (3-2 lc). The need for transposing will become when we expand the terms in Eq. (3-2 lc). The last term denotes sum-

(energy) term clear

mation,

P uv

l

+P

M=

2 iV 2 , if

2.

EXPANSION OF TERMS

We now them

consider the

first

three terms in Eq. (3-2 le) one by one

and expand

as follows:

First term

"-*!>
1

—

{v\

2v v 2 l

+

..{:)-

v\)dL.

(3-23a)

Second term:

AIY

W Here

Y

is

assumed

fy

for uniform

[y(l

to be

-

L)v,

+ i-(l +

L)v 2 ~\^dL.

(3-23b)

uniform gravity load per unit volume. Similarly,

the third term yields

W

pi

Finally, the

'

C

¥

sum of

Eq. (3-2 le) gives

J',

[t (1 "

L)Vl

+

T

(1

+ L >*\ dL

-

(3 " 23c >

Eqs. (3-23a), (3-23b), and (3-23c) and the last term in

One- Dimensional Stress Deformation

Chapter 3

¥l',[i Notice that

2) to find its

+

invoke the principle of

minimum

+ L>

2

|
-

C P,^.

if J

minimum

.

potential energy (see Chapter x

\\

2

=

i-(l

value by differentiating E^, with respect to v and v 2 as

Tyl

gig

L) Vl

a quadratic or second-order function in Vj and v 2

II p is

Now we

-

(1

53

T

(1

<" 2 ^

" +

L)dL

"

2v ^ dL

Pl/

=

~ ^f\\

t

(3 " 25a)

°'

T

(1

+

L)dL

X

~

^2

T

\

Although we used

(1

+

L)flfL

~

P2/

=

(3 " 25b)

°-

and partial differentiations to find advanced applications to multidimensional problems, it is more concise and convenient to use calculus of variations; then both Eq. (3-25a) and Eq. (3-25b) are written together as

the

minimum of lip,

differential calculus

for

SU P = where 6

is

the variation notation.

We

(3-26)

0,

shall briefly discuss variational princi-

ples at various stages, but the majority of derivations in this

book

will

be

obtained by using differential calculus.

INTEGRATION For the one-dimensional problem, the integrations in Eq. (3-25) are simple; in fact, the first term in these equations is independent of the coordinate L

;

hence integrations involve only constant terms. Thus,

aeV 41

L

1

(2»,

-

2v 1 )dL

= ^(2t>, -

-1

=
4\

2v 2 )L

4) AIY 2

One- Dimensional Stress Deformation

54

Chapter 3

Integrations of other terms in Eq. (3-25) lead to

^0>, ^fi-Vi

v2)

+v

2

-Mf-TA- P u =

)-^-^- P

2l

0,

(3-27a)

= 0.

(3-27b)

In matrix notation,

AE

1

-1

I

1

1

or

= {Q}.

M{q} Here it is

[k]

=

stiffness

(3-28b)

matrix of the element, and with the linear approximation

identical to the matrix of stiffness influence coefficients in matrix struc-

tural analysis;

{Q}

=

it is composed of body and the joint loads. It is

element nodal load vector, and

force (due to gravity), surface traction or loading,

interesting that the finite element derivations with the linear result in

approximation

AIY

lumping of the applied loads (equally) at the two nodes:

= weight of the element, where Al = volume of the element, and fy = total l

surface load.

We note that if higher-order approximations were used, the load come out

vector would not necessarily

as

lumped loading; we

shall consider

such cases at later stages.

Equation (3-28) provides a general expression that can be used repeatedly to find stiffness relations for

In the foregoing,

we

all

elements in the assemblage.

derived element equations by differentiating the

expanded function for Tl p This is possible because there are only two variables, v and v 2 Advanced problems, however, involve a large number of variables, and we commonly write the results directly in terms of matrix equations. For instance, differentiation of II P in Eq. (3-21e), .

.

x

dn£

_Q

*n*-o. dv 2 leads to

*L

[BY (2

x

[B]

[C] 1)(1

x

1)(1

x

dL{q}

=

Al 2

(2

2)

Y dL

[Nf _1

x

1)(1

x

+

1)

dL

f,

(2x

1)(1

x

1)

+

{P//} (2

/=

x

1)

1,2,

or [k]{q]

= {Q},

(3-28b)

Chapter 3

One- Dimensional Stress Deformation

55

where ai f

1

h=%

[B]

T

E[B]dL

and (Q}

The terms

in

=

y J'

[NTF^L

+ -J-J"

[Nff/L -

{P,}.

Eq. (3-28) are the same as in Eq. (3-25), except that here they The transpose in [N] r arises because after

are arranged in matrix notation.

differentiations [Eq. (3-25)] the force terms yield

order of 2 in the final results. After the process

?.

load vector that has an

on n^, the a set of (linear) simultaneous equations [Eq. (3-28b)]. Often, the stiffness matrix [k] is called an operator, which means that when operated on (nodal) displacements, the results are the (nodal) forces.

end

result

o. differentiation

is

COMMENT The primary

task in finite element formulation can be considered to

involve derivation of element equations. involve assembly, patching

The subsequent steps essentially up or combination of elements, and use of linear

algebra for solution of the resulting simultaneous equations. Step

5.

Assemble Element Equations

to

Obtain Global Equations

POTENTIAL ENERGY APPROACH Although, for simple understanding, we considered the equilibrium of a single element in the foregoing step, it is necessary to emphasize that it is the equilibrium of the entire structure in which we are interested. Consequently, we look at the total potential energy of the assemblage and find its stationary (minimum) value. The procedure of assembling element equations can be interpreted through minimization of total potential energy.

As an example, we

consider the column (Fig. 3-1) divided into three

elements with four nodes (Fig.

3-8).

Here we have numbered nodes starting from the top of the column and have measured the global coordinate y as positive downward. This is just for convenience; for instance, since the loads act downward, they are positive, and displacements are positive. It is, however, possible to measure y as positive upward or downward from any convenient point and to number nodes in other consistent fashions.

By using Eq. (3-24) we can write the potential energies H p for each element in the assemblage and add them together to obtain total potential energy IP, as

One- Dimensional Stress Deformation

56

Global

Local

Nodes

Displacement

Nodes

Properties

Displacement

S ©

©

A

Chapter 3

" y

©

© © © 1 J

v2

A

An to ,

©

©

«3

© ©

1

J

/

»2

v?

A © ©

v4

Figure 3-8 Discretization of column and numbering.

2

+

»/

—

1

^fj"'

W - 2v v 3

+ vl)dL

4

^f

[y

-

i>.

+

yO

+ L ^] rfL

y^ J'

[yd -

L) Vl

+

-1(1

+ L)v^dL

y^

[yd " L> + yd + L)»4 ]rfi

fir

3

3

J'

~

¥

-

%^ J'

/',

- P!,«, -

[t°

~~

L)v '

+ T (I + LKl rfL

[yd " iK + yd +

(PL

+ P*u)v - (P|, + 2

£)»,]<«.

P\,)v 3 (3-29)

One- Dimensional Stress Deformation

Chapter 3

57

Here we assume that A,E, /, Y, and fy are different for different elements and an element number, e denotes an element, the super= total number of elements. Note script on P denotes an element, and that the local joint or point loads at common nodes are added together and 2 yield the global joint or point loads at these points. For instance, PJ, + P U gives global point load at global node point 2. For global equilibrium we minimize 11^ with respect to all four nodal displacement unknowns, v u v 2 v 3 and v 4 Thus, their subscript denotes

M

,

dU

=

dv f

^f -

^

<2

^ -

.

,

2v ^ dL

-

^r f y 1

~

({

L)dL

1

[yO -

J

- Ph =

L)\dL

(3-30)

0,

and so on. In combined variational notation we write

5n p

(3-31a)

which denotes cttPp

dn