1.27 How to predict changes

Back

How will cracks in a beam, that is changes in the stiffness of single elements,

affect the stress distribution in a structure? How will the Dirac energy change?

This is the topic of this section.

Betti would say that the effects of the cracks can only be predicted by

meticulously evaluating at each integration point y (that is effectively the

whole structure if p = gravity load) the changes in the influence function for,

say, the bending moment M(x)

Mc(x) −M(x) =

_ l

0

[Gc

2(y, x) − G2(y, x)] p(y) dy . (1.319)

Here G2 and Gc

2 are the influence functions for M(x) in the uncracked and

Mc(x) in the cracked structure respectively.

But we want to show that there is a better approach which effectively

restricts the analysis to those members which change. This approach may be

114 1 What are finite elements?

Fig. 1.83. Bridge structure

summarized as follows: To determine the change in a quantity u(x), σ(x), ...

due to a change EI → EI +ΔEI in a certain member it suffices to evaluate

the increase/decrease in the strain energy product in that member only

Mc(x) −M(x) = −

_ x2

x1

ΔEI w

__ (Gc

2)__(y, x) dy = −

_ x2

x1

ΔEI

EI

MMc

2

EI

dy .

(1.320)

Here Gc

2 is the Green’s function for Mc(x) in the structure after the beam has

cracked and w is the deflection of the element before the beam cracked. Note

that x can be any point of the structure, that is we can effectively predict

the changes in any member by integrating over the cracked element only. Of

course the trick is that the distant member is present under the integral sign

via the Green’s function.

Remark 1.13. The bending moments of a Green’s function are very large—

even when they are “small”—but because we divide by the stiffness EI, see

(1.320), the effects cancel.

Localization

In the following we want to do the localization more systematically.

Assume the stiffness in the center span of a bridge deviates from the stiffness

EI in the outer spans by a term ΔEI, see Fig. 1.83.

The weak formulation of the original problem (uniform EI)

_ l

0

EI w

__

v

__

dx

_ _ _

a(w,v)

=

_ l

0

p v dx

_ _ _

(p,v)

v ∈ V (1.321)

and the weak formulation of the changed problem

_ l

0

EI w

__

c v

__

dx

_ _ _

a(wc,v)

+

_ x2

x1

ΔEI w

__

c v

__

dx

_ _ _

d(wc,v)

=

_ l

0

p v dx

_ _ _

(p,v)

v ∈ V (1.322)

differ only by one additional term, the integral from x1 to x2.

Hence in an abstract setting modifications of the stiffness of a structure

lead to weak formulations with an additional symmetric term d(u, v)

1.27 How to predict changes 115

uc ∈ V : a(uc, v) + d(uc, v) = (p, v) v ∈ V . (1.323)

The notation we adopted here is short for: find uc ∈ V such that ... for all

v ∈ V . In FE methods we restrict the search to the FE functions in Vh ⊂ V .

How then does the solution u of the original problem

u ∈ V : a(u, v) = (p, v) v ∈ V (1.324)

differ from the solution uc of the changed/cracked model (1.323)? Or more to

the point: how does J(uc) differ fromJ(u) where J(.) is any output functional,

that is the displacement or the stress or ...

J(u) = u(x) J(u) = σ(x) J(w) = M(x) etc. (1.325)

at a specific point?

Let eu = uc − u the difference between the two solutions. Recall that—

where applicable—

J(u) = a(G, u) ≡ w(x) =

_ l

0

EI G

__(y, x)w

__(y) dy (1.326)

and hence

J(eu) = a(G, eu) ≡ wc(x) − w(x) =

_ l

0

EI G

__(y, x) (w

__

c (y) − w

__(y)) dy

(1.327)

as well.

Obviously we have, subtract (1.324) from (1.323),

a(eu, v) = −d(uc, v) v ∈ V (1.328)

and therefore also, choose v = G,

J(eu) = a(eu,G) = −d(uc,G) (1.329)

or in terms of the beam

wc(x) − w(x)

_ _ _

J(ew)

= −

_ x2

x1

ΔEI G

__(y, x)w

__

c (y) dy

_ _ _

d(wc,G)

. (1.330)

Eq. (1.329) is the central equation.

To express the same result with u and Gc, the Green’s function of the

linear functional J(u) in the cracked model, we note that

Gc ∈ V : a(Gc, v) + d(Gc, v) = J(v) v ∈ V . (1.331)

116 1 What are finite elements?

A. Not

decay!

Hence the error in the output functional is as well

J(eu) = a(eu,Gc) + d(eu,Gc) = −d(uc,Gc) + d(eu,Gc) = −d(u,Gc)

(1.332)

so we can express J(eu) = −d(uc,G) = −d(u,Gc) both ways. Either combination

uc × G or u × Gc will do.

Remark 1.14. By applying the same arguments to the FE equations the results

can be extended to eh

u = uhc

− uh

J(eh

u) = −d(uhc

,Gh) = −d(uh,Ghc

) (1.333)

Summary

The change in any output value

J(eu) := J(uc) − J(u) (1.334)

Fig. 1.84. Influence

function for support

all Green’s functions

reaction R

1.27 How to predict changes 117

can be expressed by an energy integral, the d-scalar product between u and

Gc or—vice versa—between uc and G. Because d(Gc, u) only extends over the

elements affected by the change, see (1.320), this formula is superior to Betti’s

formula

J(eu) = − d(Gc, u)

_ _ _

local analysis

=

_ l

0

[Gc − G] p(y) dy

_ _ _

global analysis=Betti

(1.335)

or

wc(x) − w(x) = −

_ x2

x1

ΔEI G

__

c w

__

dy =

_ l

0

[Gc − G] p(y) dy

(1.336)

where the second formulation requires that we trace the deviations between

Gc and G over the whole structure—at least in the load case p = dead weight.

Recall that all quantities in linear mechanics are energies, u(x)×1, σ(x)×1,

see Sect. 1.26 p. 110, and so when the stresses and displacements change then

the Dirac energy, (u(x)+Δu)×1—that is the work done by the exterior load

on acting through the Green’s function—changes and this change in energy,

Δu × 1, is just

− d(u,Gc) = ante × post = −d(uc,G) = post × ante (1.337)

so we only need to look at the single spring, the single member, the single

plate or slab element where the change occurs

−J(eu) = d(uc,G) = Δk G(l, x)wc(l) spring k (1.338)

= ΔEI

_ x2

x1

G

__

w

__

c dy beam (1.339)

= Δt

_

Ωe

σG

ij εc

ij dΩy plate (1.340)

= t

_

Ωe

ΔCijklεG

kl εc

ij dΩy plate (1.341)

= ΔK

_

Ωe

mG

ij κc

ij dΩy slab (1.342)

to assess the change in u(x) or σ(x) at an arbitrary point x of the structure.

The stresses in the element induced by the Dirac delta—which may be

located at some very distant point x—act like weights. That is if the concrete

cracks under tension in a slab element, K → K + ΔK, then typically the

effects will be scaled by some negative power r−1, r−2, . . . of the distance r

between the element and the point x.

But please note that not all Green’s functions decay. If rigid-body motions

are involved then the opposite may be true; see Fig. 1.84 and also Fig. 1.60

p. 89.

118 1 What are finite elements?

Simplification

To apply the formula (1.332) for J(eu) we must know the solution u of the

original problem and the Green’s function Gc of the changed problem (or vice

versa). This is not very practical because once we have set up the equations

for both systems we could compare the two solutions u and uc directly.

So let us try a different approach. With the Green’s function G of the

original problem

G ∈ V : a(G, v) = J(v) v ∈ V (1.343)

we obtain the formula

J(eu) = −d(u,Gc) = −d(u,G) − d(u,Gc − G) (1.344)

or if we drop the second term

J(eu) _ −d(u,G) . (1.345)

This approximate formula has the advantage that all terms come from the

original model.

Force terms

In Sect. 7.7 p. 535 we will see that if the Green’s function G is the influence

function for a force term (a stress or any other internal action) at a point x

then the formula

J(w) = a(G,w) ≡ M(x) ?=

_ l

0

M2M

EI

dy = 0 (1.346)

makes no sense—at least not in the naive sense. We cannot calculate M(x) or

V (x) at a point x with Mohr’s integral. Rather what (1.346) means is this: if

there is a sequence of Green’s functions {Gh} which converges to G then we

have

lim

h→0

a(Gh, w) = a(G,w) + J(w) = 0+J(w) (1.347)

that is in the limit out of a(Gh, w) pops J(w) but a(G,w) itself is zero, that is

a computer cannot calculate J(w) by evaluating a(G,w) a posteriori that is

when all is done. Rather the computer must follow the action from the start

only then will it have a chance to catch J(w)—as the limit of all expressions

a(Gh, w).

But luckily in FE analysis we have Tottenham’s equation which guarantees

that the weak formulation and Betti are identical

J(wh) =

_ l

0

Gh(y, x) p(y) dy = a(Gh, wh) (1.348)

1.27 How to predict changes 119

Fig. 1.85. Cracks in a continuous beam. The integral over the shaded areas,

(M,M2) × (EI)

−1 = −13.3 kNm, is approximately ΔM = −11.84kNm

whatever the type of Green’s function G—even if J(wh) is a force term. So in

FE analysis we must not be concerned with such “subtile” distinctions. The

equation

Mh(x) =

_ l

0

Mh

2 Mh

EI

dy (1.349)

simply works—because we make it work, see (7.221) p. 535. This guarantees

also that

J(eh

u) = −d(uhc

,Gh) = −d(uh,Ghc

) (1.350)

for any functional J(.).

Example—cracks in a beam

Under the action of the point load the concrete cracks (see Fig. 1.85) so that

the bending stiffness drops from EI = 90, 625 kNm2 to EI +ΔEI = 46, 400

kNm2 which is a drop of nearly 50 % (ΔEI = −44, 225).

To predict the change in the hogging moment at support B we use the

approximation

J(ew) = −d(w,Gc

2) _ −d(w,G2) (1.351)

120 1 What are finite elements?

that is we substitute for the exact Green’s function Gc

2 (cracked beam) the

original Green’s function G2 (sans cracks) so that the change in M is approximately

the d-scalar product between the two deflections, w and G2, of the

uncracked beam

Mc −M = −

_ x2

x1

ΔEI

EI

MM2

EI

dy = −13.3 kNm (1.352)

while the true value is −11.84 kNm.

Remark 1.15. The study of the equation

J(ew) = Mc(x) −M(x) = −

_ x2

x1

ΔEI

EI

MMc

2

EI

dy = −d(w,Gc

2)

(1.353)

makes for an interesting subject.

If in a continuous beam the change, EI → EI + ΔEI, is the same in all

cross sections then [x1, x2] is the whole beam [0, l] and so d(w,Gc

2) essentially

coincides with a(w,Gc

2)

J(ew) = ΔEI

EI

_ l

0

MMc

2

EI

dy = ΔEI

EI

a(w,Gc

2) = 0 (1.354)

but because a(w,Gc

2) is zero, see Sect. 7.7 p. 535, the change in the bending

moment must be zero too. Hence the bending moment distribution in a

continuous beam—with a uniform EI—is independent of the magnitude of

EI.

Vice versa, if the bending stiffness EI in a continuous beam varies locally

then the bending moment distribution is sensitive to such variations in EI.

In a statically determinate beam the Green’s function forM(x) is piecewise

linear so that M2 = 0 and so J(ew) = 0 for any ΔEI, that is in a statically

determinate beam M(x) does not depend on EI.

Nodal form of J(eu)

The nodal form of a functional J(.) is

J(uh) = uT

GKu (1.355)

where uG is the nodal vector of the Green’s function and the nodal form of

the change J(eu) = J(uhc

− uh) in a functional is—as we will show

J(eu) _ −uT

GKΔ u. (1.356)

To start we observe that the vector-and-matrix form of the two equations

(1.323) and (1.324) is

1.27 How to predict changes 121

Fig. 1.86. Contributions of element Ω1 to the stress σxx a) stresses from the edge

load b) strains from G1

.

uc ∈ Rn vT Kuc + vT KΔ uc = vT f for all v ∈ Rn (1.357)

and

u ∈ Rn vT Ku = vT f for all v ∈ Rn (1.358)

where the matrix KΔ encapsulates the change in the stiffness matrix. It corresponds

to the d-scalar product.

By subtracting these two equations we obtain

vT K(uc − u) = −vT KΔ uc v ∈ Rn (1.359)

or if we let v = uG, the nodal vector of a Green’s function which belongs to

a functional J(.)

122 1 What are finite elements?

Fig. 1.87. Change of Young’s modulus in a single element and influence on the

vertical displacement of node N: a) nodal displacements under gravity load;

b) incremental equivalent nodal element forces fΔ, e

G(i) in element Ωe due to the

.

J(uhc

− uh) = uT

GK(uc − u) = −uT

GKΔ uc _ −uT

GKΔ u. (1.360)

In most cases the stiffness of an element Ωe changes, E → E +ΔE, and then

the d-scalar product is

d(uh

G,uh) = ΔE

E

a(uh

G,uh)Ωe (1.361)

so that the stiffness matrix of the modified structure is (E is a common factor

of all entries ke

ij)

Kmod = K +KΔ (KΔ)ij = ΔE

E

ke

ij (1.362)

where the additional matrix KΔ contains only contributions (ke

ij) from the

element Ωe weighted with ΔE/E —in all other elements ΔE = 0—and so

J(uhc

− uh) _ −d(Gh,uh) = −uT

GKΔ u = −uTKΔ uG . (1.363)

Because of the “local” character of KΔ this triple product can be identified

with the work done by the equivalent element nodal forces fΔ, e

G = Ke

Δ uG

of the change of the Green’s function on Ωe on acting through the nodal

displacements u or vice versa, that is the d-scalar product reduces to

− d(Gh,uh) = −uT

(8)(Ke

Δ)(8×8) (uG)(8) (1.364)

(e.g. bilinear element with eight degrees of freedom) or (see Fig. 1.87)

− d(Gh,uh) = −uT

(8)(fΔ, e

G )(8) . (1.365)

Dirac delta at node N

1.27 How to predict changes 123

If only one entry in the stiffness matrix changes, say the entry k11,11 →

k11,11 +Δk of an elastic support, then the change in the output functional is

simply

J(uhc

− uh) = fG

11

· Δk · u11 (1.366)

where fG

11 is the nodal force (= support reaction) in the spring due to the

action of the Dirac delta (located somewhere else) and u11 is the nodal displacement

under load.

Nonlinear problems

In nonlinear problems there are no Green’s functions and so the equation

J(uc) − J(u) = −d(z,uc) (1.367)

is not applicable. But if we substitute for the missing z the Green’s function

z at the linearization point,

aT (u; z, v) = J

_(u; v) v ∈ V (1.368)

—in an FE setting this would be

KT (u) z = j (1.369)

where KT is the tangential stiffness matrix—it follows, [50], that

J(uc) − J(u) = −d(u, z) − 1

2

{d(u, ez) + d

_(u)(eu, z) − R} (1.370)

where d_ is the Gateaux derivative of d and R is a remainder that is cubic in

eu = uc − u and ez = zc − z.

The error in the output functional of a nonlinear problem can be expressed

as, see Sect. 7.5, p. 526,

J(u) − J(uh) =

1

2ρ(uh)(z − zh) +

1

2ρ

∗(uh, zh)(u − uh) + R(3)

h (1.371)

where the first term

ρ(uh)(z − zh) = (p − ph, z − zh) (1.372)

is the work done by the residual forces on acting through the error, z −zh, in

the Green’s function and the second term is—in a somewhat symbolic notation

ρ

∗(uh, zh)(u − uh) = (δ0 − δh0

,u − uh) (1.373)

h

(3)

h

remainder.

the error in the functional J(.) evaluated at u − u and R is a cubic

124 1 What are finite elements?

By combining the discretization error (1.371) with the model error (1.370)

it follows that

J(uc) − J(uh) = −d(uh, zh)

+

1

2

{(p − ph, z − zh) + (δ0 − δh0

,u − uh)}

− 1

2

{d(uh, ez) + d

_(uh)(eu, zh)} +

1

2 R (1.374)

or if we neglect higher order terms

J(uc) − J(uh) = −d(uh, zh)

+

1

2

{(p − ph, z − ih z) + (δ0 − δh0

,u − ih u)} .

(1.375)

Here uc is the exact solution of the modified problem, E → E + ΔE, and uh

is the FE solution of the original (or simplified) problem. The terms ihu and

ih z signal that uh and zh respectively can be replaced by any two functions

that interpolate u and z on Vh—or come close to u and z in what sense

ever—so that they provide a tight upper bound on the error.

To make this formula applicable the unknown functions z and u are approximated

by higher-order interpolations of the FE solutions zh and uh

z − ih z _ i(2)

2h zh − zh (1.376)

u − ih u _ i(2)

2h uh − uh (1.377)

that is if Vh exists for example of piecewise linear functions then a quadratic

interpolation may be used. For additional details see [50].

Remark 1.16. To put these results in a better perspective we add some remarks:

In some sense analysis is all about inequalities, about bounds, about

estimates or—as we might say as well—about errors. The error in a linear

interpolation uI of a regular function u is bounded by the maximum value of

the second derivative u__(x) and the element length h

|u(x) − uI (x)| ≤ h2 · max |u

__| . (1.378)

That is these two parameters control the error. The aim of any numerical

analysis must be the search for such estimates.

With the influence function for u(x)

u(x) − uI (x) = −

_ h

0

G0(y, x) (u

__ − u

__

I ) dy = −

_ h

0

G0(y, x) (u

__ − 0) dy

= −

_ h

0

__ ___ _ u

__

dy . (1.379)

1.27 How to predict changes 125

it should not be too difficult to prove that h2 · max |u__| is a bound for the

interpolation error. The proof never mentions the influence function explicitly

but it relies of course on the properties of the triangle G0(y, x), see Sect. 7.11,

p. 560.

Now in a nonlinear problem there is no Green’s function z such that

J(u) =

_

Ω

z(x, y) · p(y) dΩy . (1.380)

But we know that the error

J(u) − J(uh) (1.381)

of the FE solution is small if the error z − zh in the Green’s function at

the linearization point can be made small. This is the important observation.

The Green’s function z at the linearization point is of no direct use—in the

sense of (1.380)—but it seems reasonable to assume that if the error z−zh is

small (and the error u − uh) then also the error J(u) − J(uh) will be small.

Obviously the error in the Green’s function is—strangely enough—correlated

with the discretization error J(u)−J(uh) and also, as shown above, with the

modeling error.

For more on the discretization error J(u) − J(uh) in nonlinear problems,

see Sect. 7.5 p. 526.

Linear versus nonlinear

The variational formulation of a nonlinear problem, see Sect. 4.21 p. 401,

a(u, v) :=

_

Ω

Eu(v) •S(u) dΩ = (p, v) (1.382)

and a linear problem

a(u, v) :=

_

Ω

E(v) •S(u) dΩ = (p, v) (1.383)

differ only in the term

d(u, v) =

_

Ω

(Eu(v) − E(v)) •S(u) dΩ

=

_

Ω

(∇uT∇v + ∇vT ∇u) •S(u) dΩ (1.384)

and so it should not be too far-fetched to assume that nonlinear effects can

be predicted by the formula

J(uc) − J(u) _ −d(z,uc) = −

_

Ω

(∇uTc

∇z + ∇zT ∇uc) •S(uc) dΩ

(1.385)

126 1 What are finite elements?

c

element Ωe is differc

where the Green’s function z is taken from the linear model and uc is the

solution of the nonlinear problem.

In FE analysis we would substitute also for uc the FE solution of the linear

problem and so

J(uc) − J(uh) ≈ −d(zh,uh) = −

_

Ω

(∇uTh∇zh + ∇zTh

∇uh) •S(uh) dΩ .

(1.386)

Навигация