7.9 Shifted Green’s functions

Back

In more abstract terms most properties of an FE solution are based on a

Shifted Green’s function theorem, that we want to formulate in the following.

The model boundary value problem is the Poisson equation:

Δu = p inΩ u= 0 on Γ . (7.250)

The associated identities are

G(u, ˆu) =

_

Ω

Δu ˆudΩ +

_

Γ

∂u

∂n

ˆuds −

_

Ω

∇u •∇ˆudΩ = 0 (7.251)

and

B(u, ˆu) =

_

Ω

Δu ˆudΩ

_

Γ

∂u

∂n

ˆuds −

_

Γ

u

∂ˆu

∂n

ds

_

Ω

u (−Δˆu) dΩ = 0. (7.252)

In the following G(u/p, ˆu) = 0 denotes the formulation of Green’s first identity

if in G(u, ˆu) the term −Δu is replaced by p and the trace of u (= boundary

value) on Γ by 0, i.e. in a first step for the left-hand side the data on the righthand

side of the boundary value problem (7.250) are substituted—wherever

possible—and in a second step the remaining slots are filled with the function

u or its derivatives, where u is the argument to the left of the slash in u/p. In

542 7 Theoretical details

Fig. 7.10. The essence of the FE methodthe substitute load ph is work equivalent

to p with respect to all virtual displacements ϕi Vh: a) original load case, b)

equivalent load case

the same sense a formulation like B(u/p, ˆu/ˆp) = 0 is understood where ˆu is

the solution of a problem

Δˆu = ˆp in Ω ˆu = 0 on Γ . (7.253)

With these substitutions, the identities become instances of the principle of

virtual displacements:

G(u/p, δu) =

_

Ω

p δudΩ

_ _ _

δWe(p,δu)

_

Ω

∇u •∇δu dΩ

_ _ _

δWi(u,δu)

= 0, (7.254)

7.9 Shifted Greens functions 543

Fig. 7.11. On Vh the approximate Dirac deltas δh

0 are perfect replacements or

proxies for the original Dirac deltas δi

and Betti’s theorem:

B(u/p, ˆu/ˆp) =

_

Ω

p ˆudΩ

_

Ω

u ˆpdΩ = 0. (7.255)

If u is the exact solution of (7.250), the expression G(u/p, δu) is the same as

G(u, δu). Things will become interesting if the exact data are mixed with the

FE solution uh as in G(uh/p, ϕi). This is what is done in FE methods.

Principle of virtual displacements

The essential property of the exact solution u is that it satisfies the principle

of virtual displacements (7.254) for any (sufficiently regular) virtual displacement

δu.

The characteristic property of the FE solution uh is that it satisfies (7.254)

only with regard to the trial functions ϕi ∈ Vh (see Fig. 7.10):

G(uh/p, ϕi) =

_

Ω

pϕi dΩ

_

Ω

∇uh •∇ϕi dΩ = 0 ϕi ∈ Vh (7.256)

but not for arbitrary admissible virtual displacements δu (admissible means

δu = 0 on Γ)

G(uh/p, δu) =

_

Ω

p δudΩ

_

Ω

∇uh •∇δu dΩ = 0, (7.257)

because −Δuh = p.

Equation (7.256) essentially means that G(u/p, ϕi) = 0 and G(uh/p, ϕi) =

0 are the “same”. In the expression G(u/p, ϕi) = 0 the exact solution u can

be replaced by the FE solution uh.

544 7 Theoretical details

An immediate consequence of the equivalence between p and its FE counterpart

ph on Vh is that the approximate Dirac deltas δh

0 are perfect substitutes

for the exact deltas on Vh; see Fig. 7.11.

Betti’s theorem

By similar reasoning Betti’s theorem can be extended to FE solutions. Given

the solutions u1 and u2 of two model problems

Δu1 = p1 u1 = 0 on Γ , −Δu2 = p2 u2 = 0 on Γ (7.258)

it follows that

B(u1, u2) = 0 or

_

Ω

p1 u2 dΩ =

_

Ω

p2 u1 dΩ , (7.259)

and the message is that in the last equation the exact solutions u1 and u2 can

be replaced by their FE counterparts uh1

and uh2

,

_

Ω

p1 uh2

dΩ =

_

Ω

p2 uh1

dΩ , (7.260)

which means that B(u1/p1, u2/p2) = 0 and B(uh1

/p1, uh2

/p2) = 0 are on Vh

the “same”. The proof rests on the Equivalence Theorem (see Eq. (7.402)),

_

Ω

p1 uh2

dΩ =

_

Ω

ph1

uh2

dΩ

_

Ω

p2 uh1

dΩ =

_

Ω

ph2

uh1dΩ (7.261)

and Betti’s theorem,

B(uh1

, uh2

) =

_

Ω

ph1

uh2

dΩ

_

Ω

ph2

uh1

dΩ = W1,2 −W2,1 = 0. (7.262)

Hence (7.260) means that the FE solutions can serve on Vh as “proxies” for

the exact solutions; see Fig. 7.12.

Principle of virtual forces

Here the sequence of functions is reversed2. The auxiliary state ˆu and the

virtual forces ˆp (the right-hand side of ˆu) comes first

G(ˆu/ˆp, u) =

_

Ω

ˆpudΩ +

_

Γ

∂ˆu

∂n

uds −

_

Ω

∇ˆu •∇udΩ = 0 (7.263)

and the rule is now the following: (i) if ˆuh ∈ Vh is the FE solution of a load

case ˆp, and if (ii) the second argument u lies in Vh (therefore we write uh

2 We mention this principle only to be complete. The result essentially is contained

in the previous formulations

7.9 Shifted Greens functions 545

Fig. 7.12. Bettis theorem and its extension to FE solutions. Tottenhams equation

is the most prominent application of this extension

instead of u), then for the virtual forces ˆp may be substituted the FE forces

ˆph

G(ˆu/ˆph, uh) =

_

Ω

ˆph uh dΩ +

_

Γ

∂ˆuh

∂n

uh ds −

_

Ω

∇ˆuh •∇uh dΩ = 0.

(7.264)

To prove this rule, we calculate the deflection ϕi ∈ Vh of a beam by applying

a “virtual force” P = 1 at a point x:

G(G0/δ0, ϕi) =

_ l

0

δ0 ϕi dy −

_ l

0

M0Mi

EI

dy = 0. (7.265)

Next recall that

a(G0 − Gh0

, ϕi) = 0 (7.266)

and that

a(G0, ϕi) = (δ0, ϕi) a(Gh0

, ϕi) = (δh

0, ϕi) . (7.267)

Hence it follows that

546 7 Theoretical details

G(G0/δh

0, ϕi) =

_ l

0

δh

0 ϕi dy −

_ l

0

M0Mi

EI

dy = 0, (7.268)

which means that

_ l

0

δ0 wh dy =

_ l

0

δh

0 wh dy wh ∈ Vh . (7.269)

Hence any point value wh(x) of a function wh ∈ Vh is equal to the work

done by the approximate load δh

0 (y − x) acting through wh(y), or stated

otherwise, on Vh the kernel δh

0 is a perfect replacement for the kernel δ0.

This is a truly remarkable result, which of course also holds for the higher

Dirac deltas.

Note that the notation G(G0/δh

0 , uh) implies that in the extension of the

principle of virtual forces the FE load case ˆp is substituted for ˆp but that G0

is left untouched! This is just the reverse of the previous substitutions.

Summary

The practical importance of these three extensions, (7.256), (7.260), and

(7.264) is that probably all post-processing in mechanics is applied duality,

is based on Green’s first or second identity:

G(u/p, ˆu) = 0 (principle of virtual displacements) (7.270)

G(ˆu/ˆp, u) = 0 (principle of virtual forces) (7.271)

B(u/p, ˆu/ˆp) = 0 (Betti’s theorem) . (7.272)

Thus to extract information from the solution for u, the exact solution is

substituted (and therewith the right-hand sides p and 0, etc.), and the place

of ˆu is taken by appropriate auxiliary functions. The function ˆu can be a rigid

body motion ˆu = 1 so that

G(u, 1) =

_

Ω

Δu · 1 dΩ +

_

Γ

∂u

∂n

· 1 ds = 0, (7.273)

provides the sum of the vertical forces (ˆu = 1 would then be called a generalized

Green’s function) or it can be a genuine Green’s function if, say, the stress

σ(x) = ∇u •n in the membrane at a specific point in a particular direction

(n)

σ(x) =

_

Ω

G1(y, x) p(y) dΩy (7.274)

is to be calculated. This equation is identical to B(u,G1[x]) = 0. Similarly, the

unit-dummy-load method of structural mechanics, which is used to calculate

the deflection of a beam at a specific point x

7.9 Shifted Greens functions 547

Fig. 7.13. The two approaches to FE analysis

w(x) =

_ l

0

M(y)M2(y, x)

EI

dy M2 = −EI

d2

dy2 G2(y, x) , (7.275)

is identical to G(G2[x], w) = 0 (principle of virtual forces).

Hence the (original) Green’s function Gi and the generalized Green’s function

z allow us to extract information from the exact solution via Green’s

identities.

If we simply speak of Green’s function—and drop the artificial distinction

between original Green’s functions and generalized Green’s functions—we can

formulate the following theorem:

Shifted Green’s function theorem: the FE solution satisfies all identities

or tests with regard to the projections of the Green’s functions

G(uh/p,Ghi

) = 0 G(Gi/δh

i , uh) = 0 B(uh/p,Ghi/δi) = 0, (7.276)

where the projections Ghi

∈ Vh are the FE approximations of the Green’s

functions.

To appreciate this theorem the reader must understand the importance of

the Green’s identities and Green’s functions for structural mechanics. When

we say that the support reactions maintain the equilibrium with the applied

load then this actually means that

G(u/p, ˆu) = 0 ˆu = a + bx = rigid-body motion . (7.277)

Any property that we are used to attribute to the exact solution as the satisfaction

of the equilibrium conditions or the fact, that the deflection of a

cantilever beam carrying a point load P is w(l) = P l3/(3 EI), is a result

that can be reproduced by substituting for u the exact solution and for ˆu an

appropriate Green’s function into Green’s identities. And all what we do in

FE analysis is that we create a shadow world Vh where anything which is true

in the real world V is true as well if only we consequently substitute for the

exact Green’s function the projections Ghi

.

548 7 Theoretical details

Fig. 7.14. Four load cases

That is with any mesh (or trial space Vh) we can associate a shift-operator

which maps the exact Green’s functions Gi[x] onto functionals Ghi

[x] in the

dual space H_

m(Ω) and the “only” problem with the FE method is that the

algorithm considers these shifted Green’s functions to be the real Green’s

functions.

Just as a change in the elastic parameters or a change in the stiffness of a

support effects a shift in the Green’s functions so an FE mesh produces a shift

of the Green’s functions. This seems to be the whole point of FE analysis.

Two approaches

From the standpoint of a reviewing engineer the FE method basically can be

classified in two ways. Both are depicted in Fig. 7.13.

• In the first approach the FE solution is identified with the solution of an

equivalent loadcase ph and the displacements and the stresses

uh(x) =

_ l

0

G0(y, x) ph dy σh(x) =

_ l

0

G1(y, x) ph dy (7.278)

are the scalar product between the exact Green’s functions and the equivalent

load.

• In the second approach the same displacements and stresses

uh(x) =

_ l

0

Gh0

(y, x)pdy σh(x) =

_ l

0

Gh1

(y, x)pdy (7.279)

are the scalar product between the approximate Green’s functions (the

shifted Green’s functions) and the original load.

7.9 Shifted Greens functions 549

In the first approach the load p is replaced by a work-equivalent load ph and

in the second approach the Green’s functions Gi are replaced by equivalent

Green’s functions Ghi

, where equivalent means that on Ph (= the set of all

FE load cases) the two coincide; (Gi, pj) = (Ghi

, pj) for every unit load case

pj . Note that a change in the Green’s function is the same as a change in the

governing equation.

Proxies

As mentioned earlier, the kernel δh

0 is a perfect substitute on Vh for the kernel

δ0, namely

wh(x) = (δ0, wh) = (δh

0, wh) = (δh

0 ,

_

i

ϕi ui)

=

_

i

(δ0, ϕi) ui =

_

i

ϕi(x) ui . (7.280)

The inner workings of this result are best understood by studying a two-span

beam, which is subdivided into two elements; see Fig. 7.14 a. The influence

function for the deflection at the center x of the first span is the deflection

curve G0[x] if a single force P = 1 is applied at x. To solve this load case on Vh,

the deflections of the nodal unit displacements at the point x must be applied

as equivalent nodal forces fi; see Fig. 7.14 b. This strange rule—deflections

become equivalent nodal forces—is easily understood if this substitute Dirac

delta

δh

0 (x − y) = {f1, f2, f3, f4, f5, f6}

= {ϕ1(x), ϕ2(x), ϕ3(x), ϕ4(x), 0, 0} (7.281)

is applied to a function wh ∈ Vh. Then indeed the value of wh at x is recovered,

_ 2l

0

δh

0 (x − y)wh(y) dy

= f1 wh(0) − f2 w

_

h(0) + f3 wh(l) − f4 w

_

h(l) + 0wh(2 l) + 0w

_

h(2 l)

= −ϕ2(x)w

_

h(0) − ϕ4(x)w

_

h(l) = wh(x) , (7.282)

(a positive f2 contributes negative work upon acting through a positive rotation

w_

h(0)), because wh lies in Vh and therefore

wh(x) = [−ϕ2(x)w

_

h(0) − ϕ4(x)w

_

h(l) − ϕ6(x)w

_

h(2 l)]x=Їx

= [−ϕ2(x)w

_

h(0) − ϕ4(x)w

_

h(l)] . (7.283)

On the larger space V this substitute Dirac delta δh

0 will not work. Consider

for example the function w(x) = sin x. The result

550 7 Theoretical details

Fig. 7.15. Maxwells theorem can be extended to the substitute Dirac deltas δh

0

_ 2l

0

δh

0 (x − y) sin y dy = −ϕ2(x) cos(0) − ϕ4(x) cos(π)

= −(−π

8

) · 1 − π

8

· (−1) = π

4

= 0.785 (7.284)

does not fit, because sin x = sin(π/2) = 1. Note that in the first span

ϕ2(x) = −x + 2 x2

π

− x3

π2 ϕ2(x) = −π/8 (7.285)

ϕ4(x) = x2

π

− x3

π2 ϕ4(x) = +π/8 . (7.286)

Maxwell’s theorem

Maxwell’s theorem, δij = δji, seems to make no sense in FE analysis, because

we cannot study the effects of true point loads with an FE program since the

program replaces any point load by work-equivalent surface loads and line

loads. But because these substitute loads δh

0 are (on Vh) a perfect substitute

for the original Dirac deltas, it follows that

wh

2 (x1) = (δ(1,h)

0 , wh

2) = (δ(2,h)

0 , wh

1) = wh

1 (x2) . (7.287)

7.9 Shifted Greens functions 551

Here δ1,h

0 is the assemblage of loads that simulate the action of P = 1 at x1

(see Fig. 7.15 a; simply the nodal forces fi in this 1-D problem), and δ2,h

0 has

the equivalent meaning (see Fig. 7.15 b).

Of course linear algebra provides the same result: let

Kui = ei Kuj = ej (7.288)

then

δij = uTj

Kui = uTi

Kuj = δji . (7.289)

And also the extension to arbitrary other pairs of Dirac deltas (that’s what

the point loads are after all) is evident.

Classical Maxwell

δ0

1,2 =

_

Ω

G0(y, x1) δ0(y x2) dΩy =

_

Ω

G0(y, x2) δ0(y x1) dΩy = δ0

2,1

(7.290)

or, to keep it short,

δ0

1,2 = (G0[x1], δ0[x2]) = (G0[x2], δ0[x1]) = δ0

2,1 (7.291)

—the superscript 0 stands for displacement—is simply Betti (L is the selfadjoint

differential operator)

(G0[x1], δ0[x2]) = (G0[x1],LG0[x2]) = (LG0[x1],G0[x2]) = (δ0[x1],G0[x2]) .

(7.292)

The extension to arbitrary pairs {i, j} of Green’s functions amounts to

δi1

,2 = (Gi[x1], δj [x2]) = (Gj [x2], δi[x1]) = δj

2,1 . (7.293)

Let for example i = 3 and j = 2 then the equation

δ3

1,2 = (G3[x1], δ2[x2]) = (G2[x2], δ3[x1]) = δ2

2,1 (7.294)

means that the shear force at the point x1 of the influence function for the

bending moment at the point x2 in a slab is the same as the bending moment

at the point x2 of the influence function for the shear force at the point x1.

Equ. (7.294) must be read as: G3[x1] picks from the “load”δ2 the shear

force (3) of the field G2[x2] (=Lδ2) at point x1, etc..

The extension of this result to finite elements is more or less obvious (we

skip the details)

δi,h

1,2 = (Ghi

[x1], δh

j [x2]) = (Ghj

[x2], δh

i [x1]) = δj,h

2,1 . (7.295)

552 7 Theoretical details

Are all results of point loads adjoint?

Consider the following question: in load case 1 a single force P1 acts at a point

x1 and in load case 2 a force P2 at a point x2. Are the stresses caused by P1

at the foot of P2, say σ(1)

xx (x2), the same as the stress σ(2)

xx (x1) caused by P2

at the foot of P1? No, this is not true

σ(1)

xx (x2)             = σ(2)

xx (x1) . (7.296)

To see this letKu1 = f1 andKu2 = f2 two FE solutions then the symmetry

of K implies that

uT2

f1 = uT1

f2 (7.297)

which is Betti. If the vectors f are the equivalent nodal forces of Green’s

functions then we have Maxwell

(uG2

)T fG1

= (uG1

)T fG2

. (7.298)

So Betti is “two” and Maxwell is “two” but in (7.296) we operate with four

states

σ(1)

xx (x2) = (u(P1))T fGσ2

            = (u(P2))T fGσ1

= σ(2)

xx (x1) (7.299)

and so we cannot use the symmetry of K to switch sides. Why, after all,

should the value of a functional J1(u2) at u2 be the same as the value of a

second functional J2(u1) at u1? The numbering scheme, 1,2, alone is no proof.