7.4 Generalized Green’s functions

Back

In an abstract sense the right-hand side of the weak formulation: find a function

u ∈ V such that

a(u, v) = p(v) for all v ∈ V (7.113)

is a functional on V . Hence we are tempted to associate with any functional

J(v) an element z of V in the sense that

a(z, v) = J(v) for all v ∈ V . (7.114)

We have switched names, u → z (= generalized Green’s function) and p(v) →

J(v), to adopt the notation used in the literature.

Because of Green’s first identity—here for a bar

G(v, z) =

_ l

0

−EAv

__

z dx + [N z]l

0

_ _ _

p(z)

_ l

0

EAv

_

z

_

dx

_ _ _

a(v, z)

= 0 (7.115)

the strain energy product can be replaced by p(z) so that

p(z) = a(z, v) = J(v) v ∈ V . (7.116)

That is, the value of the functional J(v) is identical to p(z) (= work done by

the applied load p = {−EAv__,N(l),N(0)} on acting through z where p is the

load case that belongs to v. If, for example, v is the solution of the problem

− EAv

__ =p N(l) =P v(0) = 0 (7.117)

and z is the Green’s function then

p(z) :=

_ l

0

p z dx + P · z(l) . (7.118)

Note that (7.116) is simply Betti’s theorem with the symmetric strain energy

product in between.

An engineer would say that z is the solution of the load case J or rather

δ0 as in the case J(u) = (δ0, u) and z is the Green’s function G0

p(G0) = a(G0, v) = J(v) = (δ0, v) v ∈ V . (7.119)

The advantage of this abstract approach (7.114) is, that we can associate with

any functional—not just the Dirac deltas

J(v) =

_

Ωe

v dΩ J(v) =

_

Ω

δ0 v dΩ J(v) =

_ l

0

σyydx (7.120)

520 7 Theoretical details

a Green’s function z in the sense that J(u) = p(z).

To make the mechanism more transparent, let L be a linear operator, let

L be the adjoint operator, let u be the solution of Lu = p and assume, that

j is some functional, which by pairing it with u yields a result J(u) = (u, j).

Let z be the solution of the adjoint problem Lz = j then

J(u) = (u, j) = (u,L

z) = (Lu, z) = (p, z) . (7.121)

In the case j = δ0 and z = G0, for example, we have

u(x) = (u, δ0) = (u,L

G0) = (Lu,G0) = (p,G0) . (7.122)

Because of Tottenham’s equation (1.210), p. 64, we know that the FE program

evaluates J(uh) by substituting for z the approximate generalized Green’s

function zh

J(uh) = p(zh) . (7.123)

Hence, the more accurate zh the more accurate J(uh). The strategy of the

goal oriented recovery or simply duality technique can then be summarized as

follows:

• Say you are interested in some point value or integral value of the solution.

• Interpret the value as the result of a functional applied to the solution

J(u) = (δ0, u) point value (7.124)

J(u) =

_ l

0

σyy dx integral value (7.125)

which implies that there is a generalized Green’s function z such that

J(u) = p(z).

• Formulate two weak boundary value problems: one for the original solution

u and one for the generalized Green’s function z and determine the

corresponding FE solutions

a(uh, ϕi) = p(ϕi) ϕi ∈ Vh , (7.126)

a(zh, ϕi) = J(ϕi) ϕi ∈ Vh . (7.127)

• Calculate on each element Ωe error indicators η(p)

e and η(z)

e for the two

problems, multiply the two indicators ηe = η(p)

e · η(z)

e and refine the mesh

• By following this procedure the FE result J(uh) is automatically improved.

• Note that it is not necessary to actually calculate p(zh)

J(uh) = p(zh) , (7.128)

because according to Tottenham’s equation a direct evaluation of the FE

solution yields the same result

J(uh) = (δh

0 , u) = uh(x) J(uh) =

_ l

0

σh

yy dx . (7.129)

where ηe ≥ TOL (some tolerance); see Fig. 7.5.

7.4 Generalized Greens functions 521

1828724 kN

2586207 kN

1828724 kN

2586207 kN

2586207 kN

2586207 kN

2586207 kN

2586207 kN

1828724 kN

1828724 kN

10.00

4.00

xy:

a) initial mesh, b)

mesh

If we can associate with each functional a generalized Green’s function, the

basic FE statement

a(u, v) = p(v) v ∈ V (7.130)

implies that the equilibrium position u of a structure is the generalized Green’s

function of the functional p(u) = (p, u) where p is the applied load. That is,

the whole concept of a generalized Green’s function is simply an application of

the Riesz’ representation theorem: for each linear (bounded) functional J(v)

there is an element z ∈ V such that a(z, v) = J(v) for each v.

The engineer’s version would go like this: for each load p there is a strong

solution Lu = p, where L is the differential equation. In FE methods the load

p becomes a functional p(u) and the strong solution becomes a weak solution,

a(u, v) = p(v). Hence for each functional p(v) there is a weak solution u. Now

we extend this approach to just any (linear and continuous) functional J(v),

which not necessarily must be associated with a load case p, and we claim, that

for any such functional there is a weak solution z such that a(z, v) = J(v)

for every v ∈ V . Finally we apply integration by parts, so that the virtual

strain energy a(z, u) becomes virtual external work a(z, u) = (p, z), and so

J(u) = (p, z) where Lu = p.

Essentially it is again Betti’s theorem which allows to proceed from the

symmetric middle term a(z, u) in either direction

(p, z)

←=

a(z, u)

= (J, u) = J(u) , (7.131)

Fig. 7.5. Influence

function for N

adaptively refined

522 7 Theoretical details

Fig. 7.6. Taut rope: a) unit load case pi, b) equivalent nodal forces that generate

the pseudo rigid-body motion wh = 1, c) a distributed load, d) equivalent nodal

forces

where J = Lz is the “load” (as for example J = δ0) in the functional J(u).

That is (J, u) = (δ0, u) with J = Lz is the “load form” of the functional

and J(u) = u(x) is the abstract form. The abstract form is that version where

we directly evaluate the solution u while the load forms—actually there are

two

J(u) = (δ0

, u) = (Lz, u) = (z,L

u) = (z,p↑

) , (7.132)

are those versions where we evaluate J(u) indirectly by forming e.g. the scalar

product between z and p.

7.4 Generalized Greens functions 523

Examples

Consider the equation −Lu = p of 2-D elasticity, see (7.83) p. 513. The

functional which provides the sum of the horizontal forces,

_

H, in a patch

Ωp is J(u) = a(u, e1)Ωp, because

G(u, e1)Ωp = p(e1)Ωp

− a(u, e1)Ωp

=

_

Ωp

p e1 dΩ +

_

Γp

t e1 ds

_ _ _ _

H

−a(u, e1)Ωp = 0, (7.133)

so that the generalized Green’s functions z is the solution of the variational

problem

a(z, v) = J(v) = a(v, e1)Ωp v ∈ V . (7.134)

Clearly the solution is the step function z = e1 (inside Ωp and z = 0 outside

of Ωp). The FE approximation zh solves the variational problem

a(zh,ϕi) = a(ϕi, e1)Ωp ϕi

∈ Vh . (7.135)

It has the shape of the “ramp” in Fig. 1.133 p. 185, i.e., zh = e1 inside Ωp

and then it slowly—and not abruptly—drops to zero outside Ωp.

The equivalent nodal forces fi = a(ϕi, e1)Ωp are zero if the support of

the nodal unit displacement ϕi is contained in Ωp, because the external loads

which constitute the associated unit load case pi are self-equilibrated, fi =

a(ϕi, e1) = p(e1) = 0. This is best seen in the case of a taut rope; see Fig.

7.6. The sum of the three nodal forces, that constitute the typical unit load

case is zero, and therefore fi = −F +2F − F = 0. Only at the edge nodes of

the patch Ωp the balance is disturbed, because the nodal forces, that happen

to lie outside of Ωp are not taken into account, so that fi = −F + 2F = F;

see Fig. 7.6 b.

In the same sense the functional

J(v) = a(v,urot)Ωp urot = tan ϕ e3 × x (7.136)

yields the resulting moment for the right portion Ωp of the plate in Fig. 1.60,

p. 89, for a given displacement field v. Here it is assumed, that the point

x = 0 is the point about which the right portion is rotated by an angle ϕ.

The FE solution zh ∈ Vh of the variational equation

a(zh,ϕi) = a(ϕi,urot)Ωp ϕi

∈ Vh (7.137)

is not the true solution z = tan ϕ ex, and this is why the resulting moment

is wrong

Mh = 1.6 kNm =

_

Ωp

zh • p dΩ         =

_

Ω

z p dΩ = 2.5 kNm = M .(7.138)

524 7 Theoretical details

Error analysis

So any point value or integral value is obtained by applying a linear functional

J(u) to the solution

J(u) = u(x) J(u) = σxx(x) J(u) =

_ l

0

σyy dx , (7.139)

and the values we see on the computer screen are obtained by applying the

same functionals to the FE solution uh. Hence in linear problems the error of

an FE solution is simply J(e) = J(u − uh) and this error is—as we will show

in the following—“symmetric”.

The Galerkin orthogonality implies

J(e) = a(e, z) = a(e, z − zh) = a(u − uh, z − zh)

= a(u, z − zh) = p(z − zh) , (7.140)

so that, for example, in the case J(v) = (δ0, v)

J(e) = u(x) − uh(x) =

_

Ω

p (G0 − Gh0

) dΩy . (7.141)

Now we have as well

J(e) = a(e, z) = (δ0, u − uh) (7.142)

so that

J(e) =

1

2 p(z − zh) +

1

2

(δ0, u − uh) , (7.143)

that is we can look at the error J(e) in both ways, either we attribute it to

the error in the generalized Green’s function z − zh or to the error in the FE

solution u−uh; see Fig. 7.13, p. 547. According to Betti’s theorem both errors

are the same. In nonlinear problems the two errors are different, see Sect. 7.5,

p. 526.

Error bounds

Next let us derive error bounds. To this aim we consider the equation −Lu = p

of 2-D elasticity with boundary conditions u = 0. Let e = u uh the error

of the FE solution uh then

G(e,ϕk) :=

_

i

G(e,ϕk)Ωi =

_

i

{

_

Ωi

ri •ϕk dΩ +

_

Γi

ji

ϕk ds

−a(e,ϕk)Ωi

} = 0 ϕk

∈ Vh (7.144)

where on each element Ωi

7.4 Generalized Greens functions 525

ri := p + Luh ji :=

1

2

[t+ + t

−] (7.145)

are the element residuals and the evenly split jumps of the tractions between

the elements.

Next let J(v) be any functional and z ∈ V the associated generalized

Green’s function, the value of J(e) is just the virtual work of the right-hand

side of e acting through z or

J(e) =

_

i

{(ri, z)Ωi + (ji, z)Γi

} =

_

i

{(ri, z zh)Ωi + (ji, z zh)Γi

}

(7.146)

where in the last equation we made use of the Galerkin orthogonality. Note

also that because of z = zh = 0 on Γ the boundary integrals are zero on each

portion Γi ⊂ Γ.

By applying the Cauchy-Schwarz inequality if follows [22]

|J(e)| ≤ ηω :=

_

i

η(p)

i η(z)

i , (7.147)

where with the diameter hi of the single elements Ωi

η(p)

i : = (||ri||20

,Ωi +

1

hi

||ji

||20

,Γi )1/2 (7.148)

η(z)

i : = (||z zh||20

,Ωi + hi||z zh||20

,Γi )1/2 . (7.149)

Equation (7.147) is the basis of the duality approach or the goal-oriented

techniques, see Sect. 1.31, p. 156.

The two special functionals

J(v) : = a(v, e)

||e||E

→ J(e) = a(e, e)

||e||E

=

||e||2

E

||e||E

= ||e||E (7.150)

J(v) : =

(v, e)

||e||0

→ J(e) =

(e, e)

||e||0

=

||e||20

||e||0

= ||e||0 (7.151)

coincide at v = e with ||e||E and ||e||0 respectively and consequently the dual

weighted residual error estimate (7.147) can be applied to these global norms

as well.

In a classical L2-error estimate appears a constant cI , that reflects the

interpolation properties of the space Vh and a global stability constant cS

||e||0 ≤ ηL2 := cI cS

3

_

i

h4i

(η(p)

i )2

41/2

cS := ||a(z, z)||E,Ω , (7.152)

which represents the global energy norm of the generalized Green’s function,

while in the duality approach

526 7 Theoretical details

||e||0 ≤ ηωL

2 := cI

_

i

h2i

η(p)

i η(z)

i η(z)

i := ||a(z, z)||E,Ωi (7.153)

the weights η(z)

i reflect the local character of the generalized Green’s function

and “it can be beneficial to keep the dual weights within the error estimator

rather than condensing them into just one global stability constant” [22].