Barber J.R. Intermediate Mechanics of Materials

Подождите немного. Документ загружается.

A.4 Finite element solutions for the bending of beams 585

As we might expect, the total load w

∆

is shared equally between the adjacent nodes,

but less obviously, equal and opposite nodal moments are generated as well. The

distributed loading and the equivalent nodal loading for this case are illustrated in

Figure A.12 (a,b) respectively.

(a)

∆

w per unit length

i+1

(b)

∆

w ∆

Figure A.12: Uniform loading of a single element: (a) continuous loading; (b) equiv-

alent nodal forces and moments

If the whole beam is subjected to a uniformly distributed load, the equivalent

nodal forces and moments of equations (A.84) will be superposed to yield the dis-

crete loading of Figure A.13. Notice that each node receives half of the load on each

of the adjacent elements, whereas the moments from adjacent elements are equal and

opposite and hence cancel each other out except at the two terminal nodes.

Figure A.13: Nodal forces and moments equivalent to uniform loading of the whole

beam

The existence of these terminal moments is perhaps best explained by reference

to the case where the beam is modelled by a single element of length, L, with two

terminal nodes, as in the following example.

∆

w ∆

i = 1

2 3 4

....

w ∆

586 A The Finite Element Method

Example A.5

A beam of length L is simply supported at its ends and subjected to a uniformly

distributed load w

per unit length. Find an approximate solution for the deﬂection,

representing the beam as a single element of length L.

Since there is only one element, the stiffness matrix is the same as K

of equation

(A.73). Furthermore, the beam is simply supported at both nodes, so u

= u

= 0,

leaving the nodal rotations

as the only degrees of freedom. We therefore delete

the ﬁrst and third rows and columns in K

, obtaining

K =

∆

4 2

2 4

, (A.85)

since here the element length

∆

Figure A.14: Nodal forces and moments for the uniformly loaded beam modelled by

a single element of length L

The nodal forces and moments corresponding to the distributed load are given

by equation (A.84) and are illustrated in Figure A.14. The nodal forces act directly

above the supports and hence produce no beam deﬂection, but the nodal moments

−M

= M

do produce deformation, which can be found from the equations

4 2

2 4







−1



, (A.86)

with solution

−

24EI

With only two rotational degrees of freedom, equation (A.65) reduces to

∗

(z) =

(z) +

(z) =

24EI



−

(L −z)

−



= −

z(L −z)

24EI

. (A.87)

w L

A.5 Two and three-dimensional problems 587

This is precisely the same result that would be obtained to this problem by applying

the Rayleigh-Ritz method of §3.6, using a third degree polynomial approximation to

the deformed shape (see Problem A.19).

A.5 Two and three-dimensional problems

The scope of this book has been restricted to systems that can be analyzed using or-

dinary differential equations and hence it is inappropriate to give any extensive dis-

cussion of the application of the ﬁnite element method in two and three dimensions,

which in the continuum formulation leads to partial differential equations. However,

the broad use of the method in engineering is largely due to the fact that it can be

extended to such problems with scarcely any modiﬁcation or added complexity.

Figure A.15: Discretization of a two-dimensional component using triangular ele-

ments

Conceptually, the simplest approach is to use a piecewise linear approximation

function, which leads naturally to a discretization with triangular elements in two

dimensions or tetrahedral elements in three dimensions. Figure A.15 shows a two-

dimensional component and a possible discretization using triangular elements. No-

tice how larger numbers of smaller elements are placed in the re-entrant corner. This

gives a better approximation to the curved boundary, but more importantly, it pro-

vides additional degrees of freedom to approximate the rapidly varying stress and

displacement ﬁelds due to the stress concentration.

Each node in two dimensions will have two degrees of freedom u

, which are

conveniently combined as the vector





Inside each triangular element, the displacement ﬁeld is described by the linear func-

tion

u = Ax + By + C (A.88)

588 A The Finite Element Method

and the conditions

= Ax

+ By

+ C (A.89)

applied at the three corners of the triangle are necessary and sufﬁcient to determine

the three vector constants A,B,C and hence deﬁne the complete piecewise linear

ﬁeld in terms of the nodal displacement vectors u

, which therefore constitute the

degrees of freedom in the formulation.

The two components of the vector shape functions v

(x,y) corresponding to this

discretization are those piecewise linear functions that satisfy equations (A.88, A.89)

with





and





respectively.

Notice the similarities between this formulation and the beam bending formula-

tion of § A.4, where the discretization also involves two degrees of freedom {u

}

at each node.

The algebraic equations deﬁning the u

can be obtained as before by developing

expressions for the strain energy U and the potential energy of the applied forces

Ω

and then applying the principle of stationary potential energy. Alternatively we can

multiply the governing partial differential equation by the piecewise linear weight

functions v

(x,y) and integrate over the domain in which the unknown function is

deﬁned — here the two dimensional domain occupied by the body, so this will lead

to a double integral. Generally it is necessary to integrate by parts to preserve appro-

priate continuity conditions and, in two and three dimensions, this requires the use

of the divergence theorem, but further discussion of this topic is beyond the scope of

this book.

A.6 Computational consi derations

The ﬁnite element method is widely used to solve problems for engineering sys-

tems of considerable complexity and it is common for the resulting discrete model

to involve very large numbers of elements and hence of unknowns in the result-

ing equations. Various techniques are avaliable for reducing the computational time

required to solve these equations, most of which take advantage of the fact that the

stiffness matrix (i.e. the coefﬁcient matrix for the corresponding system of equations)

is banded.

When the matrix is banded, it follows that each unknown (nodal value) appears in

only a few equations. It is therefore possible to eliminate any given variable from the

system with a modest number of operations. Sequential use of this technique leads to

a solution that is much more efﬁcient than matrix inversion or other techniques that

would typically be needed when the coefﬁcient matrix is full (i.e. not banded).

Example A.6

Solve the system of 5 equations

A.6 Computational considerations 589







2 1 0 0 0

1 4 1 0 0

0 1 4 1 0

0 0 1 4 1

0 0 0 1 2















































. (A.90)

The equations can be expanded as

= 4 (A.91)

+ 4C

= 2 (A.92)

+ 4C

= 3 (A.93)

+ 4C

= 4 (A.94)

+ 2C

= 2 . (A.95)

The ﬁrst unknown, C

appears only in the ﬁrst two equations, so it can be eliminated

by multiplying (A.92) by 2 and subtracting (A.91) from it, with the result

+ 2C

= 0 . (A.96)

We next eliminate C

by multiplying (A.93) by 7 and subtracting (A.96) from it, with

the result

26C

+ 7C

= 21 . (A.97)

Continued use of the same technique yields

98C

+ 26C

= 83 (A.98)

and ﬁnally

170C

= 113 , (A.99)

or C

= 0.6765. Equation (A.98) then gives

83 −26 ×0.6765

= 0.6675

and the remaining unknowns are recovered one by one from equations (A.97, A.96,

A.91) respectively as C

=0.6280, C

=−0.1794, C

=2.0897.

Anyone who has attempted to solve a system of ﬁve simultaneous equations by

hand will realize that this procedure, which only works when the matrix is banded, is

a great deal quicker than other methods. For a system of N simultaneous equations, it

involves only N matrix row operations to ﬁnd the last unknown, followed by another

N operations to recover the full set of unknowns.

This is a particularly simple example where the bandwidth of the matrix is only

three — i.e. there are only three terms in each of the algebraic equations. In more

complex problems, the band will be wider. For example it is of width 6 for the beam

bending formulation of equation (A.74). However, the computational effort involved

in this elimination procedure is still very much less than that required for direct

inversion of the global stiffness matrix.

590 A The Finite Element Method

A.6.1 Da ta storage considerations

When the system being analyzed is very complex, the amount of numerical data to be

transmitted and stored is very large and this in itself requires special consideration.

For example, the stiffness matrix contains a large number of components, but most of

these components are zeros and considerable savings can be achieved by identifying

the bandwidth of the matrix and storing only the non-zero values. We can also take

advantage of the symmetry of the matrix (K

= K

i j

) to store only half of the off-

diagonal components.

In the solution procedure for banded matrices used in Example A.6, only a few

components of the stiffness matrix are needed at each stage and this permits some

added efﬁciency in data handling. Even more efﬁciency is obtained by only assem-

bling parts of the stiffness matrix as they are needed in the solution procedure. In

this way, only a small part of the large stiffness matrix is stored in the computer at

any given stage. The resulting algorithm is known as a front solver

and it makes it

possible to solve problems involving millions of degrees of freedom.

A.7 Use of the ﬁnite element method in design

The ﬁnite element method is typically used at a later stage in the design process when

the basic geometry of the system is largely ﬁnalized and we need to establish more

accurate values for the stresses and displacements.

It is extremely widely used in

industry and the chances are very high that you will be called upon to use it at some

time in your career as an engineer. However, commercial ﬁnite element codes are

very ﬂexible and user friendly and they cover the vast majority of design needs. It is

therefore almost never necessary to write a program yourself. The primary purpose of

this appendix is to explain to the reader in a simple context some of the fundamental

processes involved in these codes and to introduce some of the terminology used.

The manuals and other supporting materials for the best known codes contain nu-

merous examples and it is often possible to modify one of these to cover the problem

under consideration, merely by redeﬁning the number of nodes and their location.

The best advice here is to start formulating the problem after the bare minimum of

introductory reading, since the output from the program will often teach you how

to correct errors or improve the model. Make sure you can reproduce the example

solution on your computer system before you make the modiﬁcations.

The error messages you will get in early attempts can be obscure and frustrating.

If you know someone in the company who has used the code before, their advice

B.M. Irons (1970), A frontal solution program for ﬁnite element analyses, International

Journal for Numerical Methods in Engineering, Vol.2, pp.5–32.

It is rather paradoxical that we use an avowedly approximate method when we want more

accuracy. The reason of course is that by using a ﬁner discretization, the accuracy of the ﬁ-

nite element method can be improved as much as we need. By contrast, the ‘exact’ methods

discussed elsewhere in this book make an implicit approximation at the beginning when the

component is idealized, for example as a beam or a shell.

A.8 Summary 591

will be invaluable. Alternatively, if you cut and paste an obscure error message into

Google or another search engine, often it will direct you to an explanation of the

problem and a way of ﬁxing it.

Once the program appears to be working satisfactorily, there is a strong temp-

tation to breathe a sigh of relief and accept the results as accurate. This is very

risky, particularly since most commercial codes have sophisticated postprocessors

that can make any results look very impressive and plausible. If the design is impor-

tant enough for you to perform a ﬁnite element analysis, it is also important enough

to make every effort to be sure that there are no errors in your formulation.

There

are three tests you might perform here:-

(i) Compare the ﬁnite element results with the analytical solution of the simplest

idealized problem you can think of that is somewhat like the real problem. The

results should certainly be of the same order of magnitude and there is reason to

be suspicious if they differ by more than a factor of 2.

(ii) Try to ﬁnd a system of loads and/or a small modiﬁcation of the geometry that

would cause the complex system to have a simple stress ﬁeld such as uniform

stress or a linearly varying bending stress. Make sure the ﬁnite element program

conﬁrms this result.

(iii) Test the convergence of the program by making several runs with different de-

grees of mesh reﬁnement, particularly in critical areas of high stress.

A.8 Summary

In this chapter, the basic principles of the ﬁnite element method were introduced

and applied to problems in one dimension. In mechanics of materials problems, the

method is equivalent to the use of the Rayleigh-Ritz method with a piecewise poly-

nomial approximation to the displacement. The approximation is deﬁned through a

ﬁnite set of nodal values that constitute the degrees of freedom for the problem.

An alternative mathematical formulation can be obtained from the governing

equilibrium equation by choosing the nodal values to minimize the integral of the

square of the error over the domain. This leads to the Galerkin approximation in

which the same functions are used as approximating (shape) functions and weight

functions. The advantage of this method is that it can be applied to any problem

governed by a differential equation.

The ﬁnite element method was applied to problems of axial loading and the bend-

ing of beams. The extension of the method to problems in two and three dimensions

was also brieﬂy discussed.

This argument cuts both ways. If you do not feel high accuracy and reliability are crucial,

it is probably not worth doing the ﬁnite element analysis, unless of course you simply want

to convince your boss of your thoroughness!

592 A The Finite Element Method

Further reading

K.J. Bathe (1982), Finite Element Procedures in Engineering Analysis, Prentice-

Hall, Englewood Cliffs, NJ.

Y.K. Cheung and M.F. Yeo (1979), A Practical Introduction to Finite Element

Method, Pittman, London.

T.J.R. Hughes (1987), The Finite Element Method, Prentice Hall, Englewood Cliffs,

NJ, Chapters 1,2.

H. Kardestuncer and D.H. Norrie, eds. (1987), Finite Element Handbook, McGraw-

Hill, New York.

Y.W. Kwon and H. Bang (1996), The Finite Element Method using MATLAB, CRC

Press, Boca Raton, LA.

G. Strang and G.J. Fix (1973), An Analysis of the Finite Element Method, Prentice-

Hall, Englewood Cliffs, NJ.

O.C. Zienkiewicz (1977), The Finite Element Method, McGraw-Hill, New York.

Problems

Section A.1

A.1. Find a piecewise linear approximation to the parabola y = x

in the range 0 <

x<1, using two equal elements and three nodes, 0,0.5, 1. Plot a graph comparing the

function and its approximation.

A.2. Find a piecewise linear approximation to the function sin(

x/2) in the range

0<x<1, using two equal elements and three nodes, 0,0.5,1. Plot a graph comparing

the function and its approximation.

A.3. Find the best straight line ﬁt to the exponential function e

in the range 0<x<1.

Plot a graph comparing the function and its approximation.

A.4*. Find the straight line ﬁt to the exponential function e

in 0 < x < 1 that mini-

mizes the percentage error,

f (x) − f

∗

(x)

f (x)

rather than the absolute error, f (x)−f

∗

(x).

Section A.2

A.5. A uniform bar of length L, cross-sectional area A and Young’s modulus E, is

supported at x= 0 and subjected to a distributed load p(x) = p

x/L per unit length.

Develop a ﬁnite element solution for the problem, using two equal length ele-

ments and hence estimate the displacement at the free end x=L.

Problems 593

A.6. The bar of Figure PA.6 is 100 mm long and has a cross-sectional area A that

increases linearly from 100 mm

at one end to 200 mm

at the other. The bar is

subjected to a tensile axial force F

=300 kN.

Use the stiffness matrix of equation (A.37) with four equal length elements to

estimate the increase in length of the bar, if Young’s modulus for the material is 80

GPa.

300 kN

100 mm

200 mm

100 mm

Figure PA.6

A.7. The composite bar of Figure PA.7 is supported at A and B and loaded only by

its own weight. The upper segment is made of steel (

g = 75 kN/m

, E = 210 GPa)

and has cross-sectional area A=100 mm

, whilst the lower segment is of aluminium

alloy (

g =27 kN/m

, E = 80 GPa) and cross-sectional area A= 200 mm

Use the ﬁnite element method with two equal length elements in each segment

to estimate the displacement at C.

100 mm

200 mm

100 mm

steel

aluminium

alloy

Figure PA.7

594 A The Finite Element Method

A.8. A uniform bar 0<x<L of length L, cross-sectional area A and Young’s modulus

E is supported at both ends and subjected to a distributed load p(x)= p

x/L per unit

length.

Develop a ﬁnite element solution for the problem, using four equal length ele-

ments, and hence estimate the support reaction at x= L.

Section A.3

A.9. The one-dimensional steady-state conduction of heat is governed by the equa-

tions

q(x) = −K

;

= Q(x) ,

where T is the temperature, q is the heat ﬂux per unit area in the x-direction, Q(x)

is the heat generated rate per unit volume and K is the thermal conductivity of the

material (here assumed constant).

In the body 0 < x < a there is uniform heat generation Q(x)= Q

and the bound-

aries x=0,a are both maintained at zero temperature, T (0)=T (a)=0. Using a piece-

wise linear approximation to the temperature with two elements and three equally-

spaced nodes, estimate the temperature at the mid-point x=a/2.

A.10. Use the shape functions (A.16) to develop an appropriate piecewise linear

approximation for the displacement u governed by equation (12.3) with boundary

conditions (12.6, 12.7), based on the subdivision of the length L into four equal

elements.

Use the method of §A.3 to obtain a set of homogeneous equations for the un-

known constants and hence estimate the buckling load P

. Compare your result with

the exact value (12.13).

A.11. It is proposed to use a piecewise linear approximation for the ordinary differ-

ential equation

′′

(x) = g(x) ; 0 < x < L ,

where g(x) is a known function and the boundary conditions are

f (0) = 0 ; f

′

(L) = 0 .

Deﬁne an appropriate piecewise approximation using N degrees of freedom and de-

velop the resulting set of simultaneous equations.

Section A.4

A.12. Figure PA.12 shows a beam of ﬂexural rigidity EI and length L, built-in at z=0

and loaded by a force F at the free end.

Develop a one element solution to this problem, noting that the displacement u

and the rotation

at the free end will be the only degrees of freedom.