1.19 Influence functions

Back

As in many engineering problems (Fig. 1.46) the accuracy of an FE solution

depends on how well the influence functions for the displacements, stresses,

or support reactions can be approximated (Fig. 1.47). The nature of these

Green’s functions is therefore to be discussed next.

All influence functions are displacements or deflections. In a beam the influence

functions for w(x), w_(x),M(x), or V (x) at a point x are the deflection

curves of the beam if a dual load, i.e., a force P = 1, amoment M = 1, a sharp

bend w_(x+) − w_(x−) = 1 or a dislocation w(x+) − w(x−) = 1, is applied at

x (see Fig. 1.34, p. 48).

In an energy method the focus is not on the loads themselves, but rather

on the work done by the loads acting through the virtual displacements. The

characteristic property of a single force P = 1 at x is that it contributes work

δw(x) upon acting through a virtual displacement δw.

A convenient symbol to describe such an action (point load) is the Dirac

delta

δ0(y − x) = 0 for all y   = x ,

_ l

0

δ0(y − x)ϕi(y) dy = ϕi(x) .

(1.231)

70 1 What are finite elements?

+

-

Fig. 1.47. The scalar product of the load p and the approximate Greens function

Gh1

is the normal force Nh of the FE solution

This is a function that is zero almost everywhere—up to the point x—and

that acting through a nodal unit displacement ϕi(y) or any other virtual

displacement contributes work ϕi(x) (= the value of ϕi(y) at x). This is what

a single force looks like in energy methods.

To have symbols for the other dual quantities, higher Dirac deltas are

introduced:

δ1(y − x) moment

_ l

0

δ1(y − x)ϕi(y) dy = ϕ

_

i(x)

δ2(y − x) sharp bend

_ l

0

δ2(y − x)ϕi(y) dy = Mi(x)

δ3(y − x) dislocation

_ l

0

δ3(y − x)ϕi(y) dy = Vi(x) ,

where Mi(x) and Vi(x) are the moment and shear force, respectively, of ϕi at

the point x.

The equivalent nodal forces fi that belong to a Dirac delta δ0 are

1.19 Influence functions 71

fi =

_ l

0

δ0(y − x) ϕi(y) dy = ϕi(x) (1.232)

and the same holds for the higher Dirac deltas. Because each Dirac delta

extracts from the virtual displacement ϕi just the term conjugate to δi, the fi

are just the values of ϕi at x conjugate to δi. Hence if the influence function

for a quantity Q(x) is to be calculated, the equivalent nodal forces fi are just

fi = Q(ϕi)(x).

Table 1.2. Equivalent nodal forces for influence functions in beams and slabs

Beam dual quantity fi Kirchhoff plate fi

w δ0 force ϕi(x) w ϕi(x)

w

_

δ1 moment ϕ

_

i(x) w,x , w,y ϕi,x (x), ϕi,y (x)

M δ2 sharp bend Mi(x) mxx,mxy,myy m(i)

xx(x),m(i)

xy(x),m(i)

yy (x)

V δ3 dislocation Vi(x) qx, qy q(i)

x (x), q(i)

y (x)

Table 1.2 lists the equivalent nodal forces fi for Euler–Bernoulli beams

and Kirchhoff plates; quantities for second-order equations are listed in Table

1.3.

Table 1.3. Equivalent nodal forces for influence functions of bars and plates

Bar dual quantity fi plate fi

u δ0 force ϕi(x) ux, uy ϕi(x)

N δ1 dislocation Ni(x) σxx, σxy, σyy σ(i)

xx(x), σ(i)

xy (x), σ(i)

yy (x)

Deflection of a taut rope

To calculate the influence function for the deflection w of the rope in Fig. 1.39

p. 58 at a node xk, a single force P = 1 is applied at xk. The equivalent nodal

forces fi are

fi =

_ l

0

δ0(xk − y) ϕi(y) dy = ϕi(xk) =

_

1 i = k

0 i        = k

(1.233)

so that f is identical to the unit vector ek, and with u = K(−1)ek the Green’s

function becomes

Gh0

(xk, y) =

_

i

ui ϕi(y) =

_

i

(K(−1) ek)i ϕi =

_

i

k(−1)

ik ϕi(y) . (1.234)

In a load case p the deflection wh(xk) at the node xk is therefore

72 1 What are finite elements?

i

wh k

_ l

0

Gh0(xk, y) p(y) dy =

_

i

k(−1)

ik

_ l

0

ϕi(y) p(y) dy

=

_

i

k(−1)

ik fi (1.235)

where the fi are now the equivalent nodal forces belonging to the distributed

load p.

b)

Fig. 1.48. a) Influence function for the normal force at the center of the bar;

(x ) =

FE approximation. The equivalent nodal forces f are the normal forces of the

nodal unit displacements at the point x = l/2

1.19 Influence functions 73

Fig. 1.49. Influence function for the FE shear force Vh at the center of the beam

Normal force in a bar

The influence function for the normal force N at the center of the bar in Fig.

1.48 is the longitudinal displacement of the bar if the bar is split at the center,

u(x+) − u(x−) = 1; see Fig. 1.48. The equivalent nodal forces fi that belong

to this load case are

fi =

_ l

0

δ 1( l

2

− y) ϕi(y) dy = EAϕ

_

i( l

2

) = Ni( l

2

) . (1.236)

Evidently the FE solution of this load case, the shape in Fig. 1.48 b, is not

the exact influence function for N(l/2).

But it is the exact influence function for the average value Na of the

normal force N(x) in the center element having end points x2 and x3 and

length le = x3 − x2. To see this, note that

Na =

1

le

_ x3

x2

N(x) dx =

1

le

_ x3

x2

_ l

0

G1(y, x) p(y) dy dx

=

1

le

_ x3

x2

_ l

0

EA

d

dx

G0(y, x) p(y) dy dx

= EA

le

_ l

0

[G0(y, x3) − G0(y, x2)] p(y) dy (1.237)

and this kernel EA/le[G0(y, x3) − G0(y, x2)] is just the shape in Fig. 1.48 b.

Namely to reproduce this kernel in Vh the following equivalent nodal forces

must be applied

fi = EA

le

_ l

0

[δ0(y − x3) − δ0(y − x2)] ϕi(y) dy = EA

le

[ϕi(x3) − ϕi(x2)]

(1.238)

74 1 What are finite elements?

which are just the forces, f2 = −EA/le, f3 = EA/le, that were applied previously

to approximate G1.

Remark 1.8. Because the normal force is defined as N(x) = EAu_(x), the

influence function G1(y, x) for N(x) is

G1(y, x) = EA

d

dx

G0(y, x) . (1.239)

So if σij = Op(u), where Op() is some differential operator and G0 is the

Green’s function for u, then the Green’s function for σij is Op(G0) and differentiation

is carried out with respect to x.

Shear force in a beam

To obtain the influence function for the shear force V (l/2) in the middle of

the beam (see Fig. 1.49) a dislocation w+ − w− = 1 must be applied so that

the equivalent nodal forces

fi =

_ l

0

δ3( l

2

− y)ϕi(y) dy = Vi( l

2

) (1.240)

are the shear forces of the nodal unit displacement ϕi at the point x = l/2.

Obviously the FE influence function is not correct. Could the situation be

saved by claiming that the FE approximation is the exact influence function

for the average value of V (x)? No, this time the previous logic fails by a narrow

margin. The influence function for the average value Va of V (x) in the center

element is

Va =

1

le

_ x3

x2

V (x) dx =

1

le

_ x3

x2

_ l

0

G3(y, x) p(y) dy dx

=

1

le

_ x3

x2

_ l

0

d

dx

G2(y, x) p(y) dy dx =

1

le

_ l

0

[G2(y, x3) − G2(y, x2)] p(y) dy

(1.241)

where 1/le[G2(y, x3) − G2(y, x2)] is almost the figure in Fig. 1.49. “Almost”

because the influence functions for bending moments have a sharp bend

at x that is not to be seen at the end points x3 and x2. But it is obviously

possible to come very close to this shape in Vh. In other words, FE solutions

gain in accuracy if averages are studied, rather than point values.

To approximate the shape of the kernel 1/le[G2(y, x3) − G2(y, x2)] in Vh,

the same equivalent nodal forces must be applied as in Fig. 1.49, because the

third-order difference quotient and the third-order derivative of a third-degree

polynomial ϕi at any point x ∈ (x2, x3) are the same:

1.19 Influence functions 75

Fig. 1.50. How a

commercial FE prothe

shear force V at

the center of a beam

result

fi =

1

le

_ l

0

[δ2(y − x3) − δ2(y − x2)] ϕi(y) dy =

_ l

0

δ3(y − x) ϕi(y) dy

= Mi(x3) −Mi(x2)

le

= Vi(x) . (1.242)

Remark 1.9. A commercial FE program calculates the influence function for

V (l/2) as follows. First it keeps all nodes fixed and it applies the dislocation

(= ϕ1) (see Fig. 1.50 c) to the second element. Next it applies the fixed end

forces (×(−1)) of this load case to the structure and adds to the resulting

deflection curve which is just Gh3

—the deflection caused by the dislocation in

element 2. The result is the exact solution (Fig. 1.50).

The deflection in the first element agrees with the exact solution

w(1)

h = −0.5 x = wexact = G3 (1.243)

but to the deflection in the second element the program adds the beam solution

ϕ1 (= dislocation at the left end of the fixed beam)

w(2)

h = −0.5 − 0.5 x + 3x2 − 2 x3

+ 1.0 − 3 x2 + 2x3 (= ϕ1)

wexact = +0.5 − 0.5 x (1.244)

influence function for

it obtains the correct

gram calculates the

(EI = 1), and how

76 1 What are finite elements?

8.00

6.00

8.00

6.00

1 2

3 4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

52 53 54 55 56 57 58 59 60 61 62 63 64 6566 67 68

69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85

86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102

103 104 105 106 107 108 109 110 111 11 2 113 114 115 116 117 118 119

120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136

137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153

154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170

171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187

188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204

205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221

Fig. 1.51. a) FE influence function for the shear force qh

x at node 99; b) FE mesh

and deflection at node 65, E = 3.0 E7 kN/m2, ν = 0.16, thickness d = 0.20 m,

element size = 0.5 m × 0.5 m

In 2-D and 3-D problems the displacement field cannot be split into a homogeneous

and a particular displacement field, so this technique is not applicable

to such problems.

Note that in the FE model of the beam we must distinguish between the

shear force V h on the left-hand side (V h

l ) and on the right-hand side (V h

r )

of the center node. Gh3

is the influence function for V h

r . This explains the

asymmetry of Gh3

(see Fig. 1.50 b).

Shear force in a slab

The influence function Gh3

,x for the shear force qx at node 99 of the slab in

Fig. 1.51 was calculated with conforming Kirchhoff elements by applying a

unit dislocation at that node. The deflection at node 65 was −0.041618 m,

which is exactly the value of qh

x at node 99 if a unit force P = 1 is placed at

node 65.

1.19 Influence functions 77

and d), and Greens function for σxx, c) and b). The scalar product of the nodal

vectors a) fi (gravity load) and b) uGi

(Greens function) gives the stress σh

xx of the

FE solution oralternativelyof the nodal vectors c) fG

i (Greens function) and

d) ui (gravity load)

Nodal influence functions

Nodal influence functions is short for nodal form of influence functions and

by this we mean that in FE analysis the evaluation of the two equivalent

influence functions

uh(x) =

_ l

0

Gh0

(y, x) p(y) dy =

_ l

0

uh(y) δ0(y − x) dy (1.245)

can be done by summing over the nodes. This is a key point of the FE method.

We have

Fig. 1.52. Nodal forces and nodal displacements of two load cases: gravity load, a)

78 1 What are finite elements?

uh(x) =

_ l

0

Gh0

(y, x) p(y) dy =

_ l

0

_

i

ϕi(y) uGi

(x) p(y) dy =

_

i

uGi

(x) fi

= uT

G f = uT

GKu = uT

GKT u = fT

G u =

_

i

fG

i (x) ui

=

_

i

ϕi(x) ui =

_ l

0

_

i

ui ϕi(y) δ(y − x) dy =

_ l

0

uh(y) δ0(y − x) dy .

(1.246)

So the displacement uh(x) is the scalar product between the nodal displacements

of the Green’s function and the equivalent nodal forces of the load

case p or—vice versa—between the nodal displacements of the FE solution

uh and the nodal forces fG

i of the Green’s function

uh(x) =

⎧⎪⎪⎪⎨

⎪⎪⎪⎩

_

i

ϕi(x) ui =

_ l

0

uh(y) δ0(y − x) dy = fT

G u

_ l

0

Gh0

(y, x) p(y) dy = uT

G f .

(1.247)

So when we evaluate uh(x) or σh(x) then we apply the corresponding Dirac

delta to each shape function ϕi and multiply the result (= fG

i ) with the weight

ui of the shape function. That is the equivalent nodal forces fG

i of the Green’s

functions are simply the displacements or stresses, etc., of the shape functions

at the point x. This is clearly seen in the middle of the long chain of equations

(1.246) because the identity

_

i

fG

i (x) ui =

_

i

ϕi(x) ui ⇒ fG

i (x) = ϕi(x) (1.248)

just means that.

For example the stress σh

xx(x) of the plate in Fig. 1.52 is equal to the

work done by the equivalent nodal forces fi on acting through the nodal

displacements uGi

of the Green’s function (a and b in Fig. 1.52)

σh

xx(x) =

_

Ω

Gh1

(y, x) • p(y) dΩy =

_

i

uGi

fi

≡ b × a . (1.249)

Because of Betti’s theorem (K = KT) this result is equivalent to

σh

xx(x) =

_

i

f G

i

ui =

_

j

σxx(ϕj)(x) uj ≡ c × d . (1.250)

The first sum extends over all nodes i = 1, 2, . . . N of the structure while the

second sum extends over all degrees of freedom j = 1, 2, . . . 2×N of the nodes.

This result means that in FE methods we calculate stresses as in finite

difference methods, see Sect. 7.6 p. 533.

1.19 Influence functions 79

The inverse of K

The entries in the inverse K

−1 of the stiffness matrix of a structure are the coefficients

gij of the projections (= FE approximations) of the n nodal Green’s

functions G0[xi], i = 1, 2, . . . n onto Vh:

Gh0

(xi, y) =

_n

j=1

gij ϕj(y) K(−1) = [gij ] . (1.251)

That is to each node i belongs a Green’s function G0(xi, y) (= influence

function for the displacement ui = u(xi) at xi) and the entries gij in row i

of K

−1 describe the expansion of the FE Green’s function in terms of the

ϕj . In 1-D problems the gij would be just the nodal values of the Green’s

function G0(xi, y) that is gij = G0(xi, yj ), j = 1, 2, . . . n, where y1, y2, . . . are

the coordinates of the nodes.

Equation (1.251) is easily verified if the analytic result

ui = uh(xi) =

_

Ω

Gh0

(xi, y) p(y) dΩy =

_

j

_

Ω

gij ϕj(y) p(y) dΩy

=

_

j

gij fj (1.252)

is compared with the computer output

ui =

_

j

k(−1)

ij fj (1.253)

and if the n unit vectors, f = ej are substituted consecutively for f. This

establishes k(−1)

ij = gij .

Linear algebra provides the same result: let u and ˆu any two vectors; then

ˆu

T Ku = uK ˆu which implies that if Ku = f, ˆuT f = uK ˆu. Next let gi

be the solution of Kgi = ei; then it follows that ui = gTi

f. If this is compared

with (1.252), it follows that the coefficients gij are the solutions of Kgi = ei,

which is equivalent to saying that K(−1) = [gij ].

Commercial codes

Commercial codes normally provide no routines for to calculate influence functions

but a diligent user can circumvent this problem. The issue is to determine

the equivalent nodal forces fi which will produce the influence function.

Let us assume that the influence function for the shear stress σxy at the

center of a bilinear element which is part of a larger structure is to be calculated.

The analysis is done on a single bilinear element which has the same shape

and size as the original element. In this case the element has eight degrees of

80 1 What are finite elements?

Fig. 1.53. Influence functions in a plate for a) the horizontal displacement ux,

b) xx

freedom. In the first load case we let u1 = 1 and all other ui = 0. The shear

stress σxy at the center of the element in this load case is the equivalent nodal

force f1, etc. So by solving eight different load cases u = ei, i = 1, 2, . . . 8, we

can calculate the eight stresses σxy(ϕi)(x) = fi which—as equivalent nodal

forces fi—produce the FE influence function for σh

yx(x).

Because the equivalent nodal forces depend only on the element which

contains the point x the formulas in Sect. 4.8 p. 357 would have provided

the stresses more easily but if the elements which are implemented are of an

unknown type or if the elements are curved or not affine to the master element

then this technique can help.

The influence functions for nodal stresses, which are usually average values

σh

xy = σ(1)

xy + σ(2)

xy + σ(3)

xy + σ(4)

xy

4

(1.254)

are obtained in the same way: it is only that the 8 × 4 nodal forces—for each

of the four element stresses σ(i)

xy at the node—must be applied simultaneously

and must be weighted with 1/4.