Barber J.R. Intermediate Mechanics of Materials

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

A.2 Axial loading 575

The nodal forces are given by equations (A.39, A.40) as

∆

L/4

p(x)xdx +

∆

L/2

L/4

p(x)



−x



∆

L/2

L/4

p(x)



x −



dx +

∆

3L/4

L/2

p(x)



−x



∆

3L/4

L/2

p(x)



x −



dx +

∆

3L/4

p(x)(L −x)dx

∆

3L/4

p(x)



x −



dx ,

where the distributed load p(x) is due to the weight of the body and is

p(x) =

gA(x) =



1 −



per unit length. Substituting this result into the above expressions and evaluating the

integrals, we obtain

F =



192



The equation system

Ku = F

can then be solved to give

u =

{0.168, 0.295, 0.376, 0.406}

and in particular, the end deﬂection is

u(L) =

0.406

This problem can in fact be solved exactly,

the resulting end deﬂection being

u(L) =

[3 −2 ln(2)]

≈

0.403

Thus, the four element approximation gives better than 1% accuracy for the end

deﬂection.

To obtain this solution, use equilibrium arguments to determine the axial force and hence

the stress as a function of x. Find the strain using Hooke’s law and substitute the resulting

expression into equation (1.3). Integration with respect to x and susbtitution of the boundary

condition u(0)=0 then yields the expression for u(x).

576 A The Finite Element Method

A.2.5 Direct evaluation of the matrix equation

In the preceding examples, we ﬁrst determined the properties of the individual ele-

ments and then assembled the system of simultaneous equations, using the pattern

of equation (A.24). With the Rayleigh-Ritz formulation, it is not necessary to go

through this preliminary step. The complete system of equations can be generated by

applying the Rayleigh-Ritz method directly to the approximation of equation (A.29).

The longitudinal strain is given for all x by

∂

∗

∂

∑

i=1

′

(x) , (A.41)

from equation (A.29) and hence the strain energy density is

E(x)e

E(x)

∑

i=1

′

(x)

E(x)

∑

i=1

∑

j=1

′

(x)v

′

(x) . (A.42)

Notice that to expand the square of the summation, we have to use a different dummy

index i, j in each sum. Integrating over the volume of the bar, we then obtain the total

strain energy

U =

∑

i=1

∑

j=1

E(x)v

′

(x)v

′

(x)A(x)dx . (A.43)

The potential energy of the loads in Figure A.5 is

Ω

= −

p(x)u

∗

(x)dx −F

∗

(L) = −

∑

j=1

p(x)v

(x)dx −F

(A.44)

and hence the total potential energy is

= U +

Ω

∑

i=1

∑

j=1

E(x)A(x)v

′

(x)v

′

(x)dx

−

∑

j=1

p(x)v

(x)dx −F

. (A.45)

The principle of stationary potential energy requires that

∂Π

∂

= 0

and hence

∑

i=1

E(x)A(x)v

′

(x)v

′

(x)dx −

p(x)v

(x)dx −F

= 0 . (A.46)

A.3 Solution of differential equations 577

This equation can be expressed in the matrix form Ku = F , where

E(x)A(x)v

′

(x)v

′

(x)dx (A.47)

p(x)v

(x)dx + F

. (A.48)

Differentiating (A.16), we have

′

(x) = 0 ; x < x

i−1

and x > x

i+1

∆

; x

i−1

< x < x

(A.49)

= −

∆

; x

< x < x

i+1

Substitution in (A.47) and evaluation of the integrals then leads directly to equations

(A.37, A.38).

Also, substitution of (A.16) into (A.48) and evaluation of the resulting integral

gives

∆

j−1

p(x)(x −x

j− 1

)dx +

∆

j+1

p(x)(x

j+1

−x)dx , (A.50)

for j 6= N and

∆

N−1

p(x)(x −x

N−1

)dx + F

, (A.51)

agreeing with (A.39, A.40).

This method of developing the solution is in some respects conceptually more

straightforward than the assembly from elemental values described in §A.2.2. How-

ever, there are practical reasons for preferring the assembly method in more complex

problems. For example, in a large problem different types of elements and hence

shape functions might be used in different parts of the structure. Efﬁcient solvers for

the resulting equations also depend on the matrix being assembled piece by piece, as

discussed in §A.6.1 below.

A.3 Solution of differential equations

In §§A.2.4, A.2.5, we developed the ﬁnite element method by applying the principle

of stationary potential energy to a piecewise linear approximation to the displace-

ment. This is still the most widely used method for problems in mechanics of mate-

rials, but an alternative method is to apply the ‘least squares ﬁt’ arguments of §A.1

directly to the governing differential equations of equilibrium. The principle of sta-

tionary potential energy is mathematically equivalent to a variational statement of

the equations of equilibrium, so the two methods lead to the same ﬁnal equations,

578 A The Finite Element Method

but one advantage of the least squares method is that it can be used to obtain an

approximate solution of any problem governed by a differential equation.

To apply the methods of §A.1.3 to differential equations, two minor modiﬁcations

are required:-

1. We need to make allowance for appropriate boundary conditions.

2. It is usually necessary to perform one or more integrations by parts to retain

appropriate continuity in the approximating function.

We shall illustrate this procedure in the context of the axial loading problem of Figure

A.5 above.

F(x) F(x + δx)

p(x) δx

δx

Figure A.8: Equilibrium of an inﬁnitesimal bar segment

We ﬁrst need to develop the differential equation of equilibrium for the bar. Fig-

ure A.8 shows the forces acting on an inﬁnitesimal segment

x of the bar. The axial

force F is assumed to be a function of x. Equilibrium of the beam segment then

requires that F(x+

x)+p(x)

x−F(x)=0 and hence

F(x +

x) −F(x)

+ p(x) = 0 .

Taking the limit as

x→0, this gives the differential equation

+ p(x) = 0 . (A.52)

The corresponding stress is

F(x)

A(x)

= E(x)e

= E(x)

(A.53)

and substitution of (A.53) into (A.52) gives



E(x)A(x)



+ p(x) = 0 . (A.54)

This is a second-order ordinary differential equation requiring two boundary condi-

tions, which in this case are u

(0)= 0 and F(L) = F

. The latter can be expressed in

terms of u

using (A.53), with the result

(0) = 0 ; E(L)A(L)u

′

(L) = F

. (A.55)

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

Equations (A.54, A.55) deﬁne the continuum formulation of the problem.

A piecewise linear approximation satisfying the ﬁrst of (A.55) in a collocation

sense is

∗

(x) =

∑

i=1

(x) , (A.56)

where the shape functions v

(x) are given by (A.16).

A set of conditions for the unknowns u

is then obtained by enforcing





E(x)A(x)u

∗′

(x)



+ p(x)



(x)dx = 0 ; j = (1,N) , (A.57)

where as in §A.1 we use the v

(x) for both shape and weight functions.

The obvious next step is to substitute (A.56) into (A.57) to obtain N simultaneous

equations for the u

, but a difﬁculty is encountered in that u

∗′′

(x) is ill-deﬁned for

the piecewise linear function (A.16). To overcome this difﬁculty, we apply partial

integration to the ﬁrst term in (A.57) obtaining

E(x)A(x)u

∗′

(x)v

(x)



−

E(x)A(x)u

∗′

(x)v

′

(x)dx

p(x)v

(x)dx = 0 ; j = (1,N) .

Substituting for u

∗

(x) from (A.56) and for u

∗′

(L) from the second of the boundary

conditions (A.55), we then have

Ku = F , (A.58)

where

E(x)A(x)v

′

(x)v

′

(x)dx (A.59)

p(x)v

(x)dx + F

jN , (A.60)

agreeing with the results (A.47, A.48) obtained using the Rayleigh-Ritz method in

§A.2.5.

A.4 Finite el ement solutions for the bending of beams

Piecewise linear approximations cannot be used for the displacement of beams in

bending, since they involve discontinuities in slope at the nodes. To avoid this dif-

ﬁculty, we use a discretization in which the nodal displacements u

and slopes

are treated as independent variables with separate shape functions of higher order

polynomial form.

For the nodal displacements u

, we require shape functions v

(z) satisfying the

conditions

580 A The Finite Element Method

) = 1 ; j = i

= 0 ; j 6= i (A.61)

) ≡ (v

)

′

) = 0 ; all j .

This corresponds to the situation where all the nodes except j = i are restrained

against displacement and rotation, whilst node j =i is subjected to unit displacement

without rotation, as shown in Figure A.9 (a).

(a)

v (z)

i-1

i+1

(b)

i-1

i+1

v (z)

Figure A.9: Shape functions for (a) nodal displacements u

and (b) nodal rotations

For the nodal slopes

, we require shape functions v

(z) satisfying

) ≡ (v

)

′

) = 1 ; j = i

= 0 ; j 6= i (A.62)

) = 0 ; all j ,

which correspond to the deformation of Figure A.9 (b), where node i is rotated with-

out displacement and all other nodes are restrained against both displacement and

rotation.

The lowest order piecewise polynomial functions satisfying these conditions are

(z) =

3(z −z

i−1

)

∆

−

2(z −z

i−1

)

∆

; z

i−1

< z < z

3(z

i+1

−z)

∆

−

2(z

i+1

−z)

∆

; z

< z < z

i+1

(A.63)

= 0 ; z < z

i−1

and z > z

i+1

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

(z) = −

(z −z

i−1

)

∆

(z −z

i−1

)

∆

; z

i−1

< z < z

i+1

−z)

∆

−

i+1

−z)

∆

; z

< z < z

i+1

(A.64)

= 0 ; z < z

i−1

and z > z

i+1

The reader can easily verify by substitution that these expressions

satisfy the condi-

tions (A.61, A.62).

If there are N elements and N +1 nodes, the discrete approximation to the beam

displacement can then be written

∗

(z) =

N+1

∑

i=1

(z) +

N+1

∑

i=1

(z) . (A.65)

Consider the i-th element of length

∆

between nodes i and i+1, in which the

displacement is approximated as

∗

(z) = u

(z) +

(z) + u

i+1

(z) +

i+1

(z) . (A.66)

The element strain energy is

i+1



∗

(z)



dz (A.67)

and if the ﬂexural rigidity EI is independent of z, we can substitute for the shape

functions v

(z),v

i+1

(z),v

i+1

(z) from (A.63, A.64) and perform the integra-

tions, obtaining

∆

(12u

+ 4

∆

+ 12u

i+1

+ 4

∆

i+1

+ 12

∆

−24u

i+1

+12

∆

i+1

−12

∆θ

i+1

+ 4

∆

i+1

−12

∆

i+1

) . (A.68)

The four columns of the 4×4 element stiffness matrix are now obtained by differen-

tiating U

with respect to the four nodal variables u

i+1

. We obtain

∂

∆

(12u

+ 6

∆θ

−12u

i+1

+ 6

∆θ

i+1

) (A.69)

∂

∂θ

∆

+ 4

∆

−6

∆

i+1

+ 2

∆

i+1

) (A.70)

∂

i+1

∆

(−12u

−6

∆θ

+ 12u

i+1

−6

∆θ

i+1

) (A.71)

∂

∂θ

i+1

∆

+ 2

∆

−6

∆

i+1

+ 4

∆

i+1

) (A.72)

More generally, the shape functions of equations (A.63, A.64) deﬁne a set of cubic inter-

polation functions that can be used in the sense of §A.1 to deﬁne a discrete approximation

to any function that preserves continuity in both value and derivative at the nodes. This is

known as a cubic spline approximation.

582 A The Finite Element Method

and hence the element stiffness matrix is

∆







12 6

∆

−12 6

∆

−6

∆

−12 −6

∆

12 −6

∆

−6

∆







. (A.73)

The global stiffness matrix can be assembled by superposition, as in §A.2.2. For

example, for a beam with 4 elements and 5 nodes, i= (1,5), we obtain

K =

∆







12 6

∆

−12 6

∆

0 0 0 0 0 0

∆

−6

∆

0 0 0 0 0 0

−12 −6

∆

24 0 −12 6

∆

0 0 0 0

∆

0 8

∆

−6

∆

0 0 0 0

0 0 −12 −6

∆

24 0 −12 6

∆

0 0

0 0 6

∆

0 8

∆

−6

∆

0 0

0 0 0 0 −12 −6

∆

24 0 −12 6

∆

0 0 0 0 6

∆

0 8

∆

−6

∆

0 0 0 0 0 0 −12 −6

∆

12 −6

∆

0 0 0 0 0 0 6

∆

−6

∆







(A.74)

In this matrix, the internal rectangles indicate the way in which the 4×4 element

matrices are overlaid along the diagonal to construct the global stiffness matrix.

The stiffness matrix of equation (A.74) contains terms for all the nodal displace-

ments that are possible with the given number of nodes, but in practice some of these

degrees of freedom must be constrained if the beam is to be kinematically supported.

If the beam is simply supported at node i, the corresponding nodal displacement

=0 and ceases to be an unknown, whilst the corresponding nodal force F

becomes

an unknown reaction. We therefore remove u

from the equation system by delet-

ing the row and column corresponding to u

from the stiffness matrix. In the same

way, the row and column corresponding to

will be deleted if node i is prevented

from rotating, as at a built in support.

Example A.4

A beam of length 2L is built in at z = 0, simply supported at z = L and loaded by a

force F at the end z = 2L, as shown in Figure A.10. Using a discretization with two

equal elements each of length L, estimate the deﬂection under the force.

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

Figure A.10

Figure A.11 shows the discretization, with 2 elements of length L and three

nodes, for which the full stiffness matrix is

K =







12 6L −12 6L 0 0

6L 4L

−6L 2L

0 0

−12 −6L 24 0 −12 6L

6L 2L

0 8L

−6L 2L

0 0 −12 −6L 12 −6L

0 0 6L 2L

−6L 4L







, (A.75)

by analogy with (A.74).

F , u

M , θ

Figure A.11

However, the boundary conditions at the supports require that u

= u

= 0

and we therefore delete the ﬁrst three rows and columns from K, leaving

K =







−6L 2L

−6L 12 −6L

−6L 4L







. (A.76)

The remaining nodal forces and moments (i.e. the external forces on the beam) are

=0,F

=−F, M

=0 and hence







−6L 2L

−6L 12 −6L

−6L 4L

























−F







, (A.77)

584 A The Finite Element Method

with solution

= −

4EI

; u

= −

7FL

12EI

;

= −

3FL

4EI

. (A.78)

Thus, the end deﬂection is 7FL

/12EI downwards.

The deleted columns of the stiffness matrix could be reinstated to determine the

unknown reactions. For example, the reaction at the simple support is

([0]

−12u

+ 6L

) =

after substituting for

A.4.1 Nodal forces and moments

In Example A.4, the only loads were applied at the nodes, so the load vector was

immediately known. If there is a distributed load, w(z) per unit length acting on the

beam, it must be shared out between adjacent nodes and this can be done using the

principle of stationary potential energy, as in §A.2.4.

The potential energy of the distributed load is

Ω

∗

(z)w(z)dz , (A.79)

where L is the length of the beam, u

∗

(z) is the discrete approximation to the dis-

placement in the entire beam deﬁned by equation (A.65) and the distributed load is

assumed to act downwards.

Differentiating with respect to the nodal displacements, we have

∂Ω

∂

(z)w(z)dz (A.80)

∂Ω

∂θ

(z)w(z)dz (A.81)

and we conclude that the distributed load contributes

= −

(z)w(z)dz (A.82)

= −

(z)w(z)dz (A.83)

to the nodal forces and moments at node i.

If there is a uniformly distributed load w

per unit length in the element z

< z <

i+1

and no load elsewhere, substitution into equations (A.82, A.83) shows that the

only non-zero nodal forces and moments occur at the ends of the loaded element and

are

= F

i+1

= −

∆

; −M

= M

i+1

∆

. (A.84)