1.31 Adaptive methods

Back

In a plate the residual forces r = p − ph and the jumps tΔh

of the traction

vector at interelement boundaries are assumed to be an appropriate measure

of the quality of an FE solution; see Fig. 1.105. Where these a posteriori error

indicators are large, the size of the elements is diminished or the order of the

polynomial shape functions is increased. Repeating this loop for many cycles,

148 1 What are finite elements?

5.54 24.91

25.41 68.41 5.74 56.1841.3721.60 10.21 5.84 28.08 5.80

10.71 15.46 13.04 22.23 7.47 20.28 23.42

51.02 14.52 22.66 10.14

22.75 21.87 13.35 19 .91 13.98 9.72 44.82

4.62 5.848.9212.73 11.92 14.32 35.03 35.94 19.52 20.08

6.00

Fig. 1.105. Plate analysis with bilinear elements. a) System and load case p (edge

load); b) load case ph; the numbers are the resulting volume forces (kN)

one hopes the solution improves; see Fig. 1.106. This is the idea of adaptive

methods.

Although this whole process seems rather straightforward, it is not quite

clear how to weight the various errors. Would it suffice to look only at the

residual forces ri and neglect the jumps tΔh

, or do the jumps carry the message?

In a beam which is subjected to a uniform load p, the nodal forces

fk = p l/(n+1) (shear force discontinuities) get smaller and smaller the more

elements n are used, but because ph = 0 the element residuals remain the

_ l

0

(p − 0)2 dx = ||p ||20

, (1.420)

same throughout the refinement, see Fig. 1.107,

1.31 Adaptive methods 149

Fig. 1.106. Adaptive refinement

1.107. The residual

the same, although the FE

exact solution for n→∞

although we are convinced that the FE solution wh tends to the exact solution

w. This is an indication that in this example the jump terms—which tend to

zero as the mesh is refined—are a more sensitive error indicator than the

element residuals. In a beam, the jump terms are the discontinuities in the

second- and third-order derivatives, fi = M(x+

i )−M(x

−

i ) and fj = V (x+

j )−

V (x

−

j ), while the element residuals measure the discrepancies in the fourthorder

derivative r = EI wIV − EI wIV

h .

The relevance of the two errors depends on the degree p of the shape

functions. For even order, p = 2, 4, . . ., the residual r dominates while for odd

order, p = 1, 3, . . ., the jump terms dominate (“dichotomy”), see [19] p. 424.

Plate analysis

Following these heuristic remarks, we now concentrate on the adaptive refinement

of plates for the details. The model problem is a plate which is subjected

to volume forces p and to traction boundary conditions t on the part ΓN of

the boundary Γ = ΓD + ΓN. The boundary value problem consists of finding

forces in each element remain

Fig.

solution converges to the

150 1 What are finite elements?

the solution u = [u1, u2]T of the system

− Lu := −μΔu − μ

1 − 2ν

∇div u = p (1.421)

where u is subject to the following boundary conditions

u = 0 on ΓD, Sn = t on ΓN. (1.422)

Let Ω =

_

i Ωi be a suitable partition of Ω into finite elements Ωi with

diameters hi and with edges Γi. The FE solution uh ∈ Vh is the solution of

the variational problem

a(uh, vh) = p(vh) :=

_

Ω

p • vh dΩ +

_

ΓN

￣t

• vh ds for all vh ∈ Vh. (1.423)

Evidently the error e = u − uh satisfies on V the equation

a(e, v) = p(v) − a(uh, v) = p(v) − ph(v) =: r(v) v ∈ V (1.424)

with

r(v) =

_

Ωi

__

Ωi

r · v dΩ +

_

Γi

j · v ds

_

(1.425)

where r are the element residual forces,

r := p − ph = p + Luh on Ωi (1.426)

and j are the (evenly distributed) jumps of the tractions at the element boundaries

j =

⎧⎪⎨

⎪⎩

1

2 (S+

h n+ + S

−

h n−) onΓi        ⊆ Γ

t − Shn on ΓN

0 on ΓD.

(1.427)

Recall that for test functions vh ∈ Vh the fundamental Galerkin orthogonality

a(e, vh) = r(vh) = 0 vh ∈ Vh (1.428)

holds. The single error terms in an element

||ri||20

=

_

Ωi

r • r dΩ ||ji

||20

=

_

Γi

j • j ds (1.429)

are weighted with h2iand hi respectively to produce a resulting error term

η2

i := h2i

||ri||20

+ hi ||ji

||20

. (1.430)

Under suitable assumptions this measure η2

i is an upper bound for the strain

energy of the field e in a single element Ωi (see Sect. 7.11, p. 565)

1.31 Adaptive methods 151

||ei||2

E := a(e, e)Ωi

≤ c η2

i (1.431)

where c is an unknown constant. Evidently

||e||2

E = a(e, e) ≤ c η2 := c

_

i

η2

i (1.432)

and so the ratio between η and the energy norm of the exact solution u defines

a global error indicator

ηrel := η

||u||E

. (1.433)

Because of u = uh +e and the Galerkin orthogonality, a(uh, e) = 0, we have,

(see Sect. 7.12, p. 568, and the “Pythagorean theorem”)

||u||2

E = a(uh + e,uh + e) = a(uh,uh) + 2a(uh, e) + a(e, e)

= a(uh,uh) + a(e, e) = ||uh||2

E + ||e||2

E (1.434)

so that

ηrel = η _

||uh||2

E + ||e||2

E

. (1.435)

The relative error in a single element is defined to be the ratio between ηi and

the energy norm ||ui||E := a(u,u)1/2

Ωi

of the exact solution in that element:

ηrel

i := ηi

||ui||E

. (1.436)

Locally the Galerkin orthogonality must not be true, a(e,uh)Ωi

            = 0, so that

￣ηrel

i = _ ηi

||uhi||2

E + ||ei||2

E

(1.437)

is not the same as (1.436), but it is an obvious approximation and to make

it computable ||ei||2

E is replaced by η2

i . Hence the idea is to refine all those

elements where the (so-modified) relative error

˜ηrel

i = _ ηi

||uhi

||2

E + η2

i

(1.438)

exceeds a certain threshold value η0 [145]. Normally only the first 30% of all

elements on the list are refined; otherwise the problematic zones will not crystallize

so well. Often the element residual forces r are also neglected because

the contributions of the jump terms are more significant.

How the energy error is typically distributed can be seen in Fig. 1.108,

where dark regions symbolize high errors and light regions small errors in the

strain energy. The critical regions are simply those with stress concentrations,

so that eventually an engineer can easily predict where the program will refine

the mesh.

152 1 What are finite elements?

Fig. 1.108. Distribution of error

indicators in a plate (gravity load)

L2-projections

Another popular way to quantify the error is to compare the raw FE stresses

with averaged stresses. This is known as the method of Zienkiewicz and Zhu

or simply as Z2 [261].

The averaged stresses are written in terms of the shape functions ψi,

￣σ

h(x) =

_n

i

￣σ

i ψi(x) or

⎡

⎣

￣σh

xx(x)

￣σh

yy(x)

￣σh

xy(x)

⎤

⎦ =

_n

i

⎡

⎣

￣σxxi

￣σyyi

￣σxyi

⎤

⎦ ψi(x) (1.439)

and the nodal values ￣σxxi (analogously for ￣σyyi , ￣σxyi ) are determined in such

a way that the square of the error ||σh

xx

− ￣σxx||20

is minimized (L2-projection),

which leads to the system of equations

_

j

_

Ω

ψi ψj dΩ ￣σxxj =

_

Ω

ψi σh

xx dΩ, i = 1, 2, . . . n . (1.440)

Often this system is only solved approximately, by diagonalizing the mass

matrix. The energy error in each element is then

1.31 Adaptive methods 153

Fig. 1.109. Oscillating longitudinal

displacements produce

oscillating stresses

η2

i =

_

Ωi

(￣σhi

− σhi ) • (￣εhi

− εhi ) dΩ . (1.441)

The method can be improved by interpolating not at the nodes, but at

the superconvergent points [262], [263]. In many cases this method is very

successful.

The weak point of the Z2-method though is the implicit assumption that

• oscillations indicate errors

• smooth stresses mean accurate stresses

and so the method breaks down if the false stresses are smooth as in Fig. 1.103,

p. 136, where pollution produced a large but smooth error in the stresses.

In the following problem [2]

− EAu

__(x) = μ sin(2m π x) u(0) = u(1) = 0, m>0 , (1.442)

a bar is stretched and compressed by rapidly oscillating longitudinal forces

having an amplitude μ. If this problem is solved on a uniform mesh of 2n, n ≤

m, linear elements so that the nodes are located at the points

xk = k

2n k = 0, 1, . . . , 2n , (1.443)

the piecewise linear FE solution is uh = 0, because we know that it interpolates

the exact solution (we let EA = 1)

u(x) = μ

4m π2 sin(2m π x) , sin(2m π xk) = sin(2(m−n) π k) = 0 (1.444)

at the nodes xk, where the exact solution is zero; see Fig. 1.109. Hence the

FE stress σh = 0 is infinitely smooth—so all seems well—while in truth the

exact stress oscillates rapidly

σ(x) = E ε(x) = μ

2m π

cos(2m π x) . (1.445)

154 1 What are finite elements?

Fig. 1.110. Green’s function for σxx at a node after an adaptive refinement

Of course, the fact that all fi = 0 are zero should have warned us, but this

problem is apt to illustrate, that deformations caused by oscillating loads are

hard to catch with finite elements; see Fig. 1.3, p. 3. The contributions to

the fi tend to cancel, and also the influence functions for the jumps in the

stresses at the nodes are symmetric (see Fig. 1.69, p. 98) so that oscillating

loads—paradoxically enough—produce smooth FE stress fields.

Additional methods

In FE adaptivity we essentially distinguish between

• h-methods, which reduce the element size

• p-methods, which increase the degree p of the polynomial shape functions .

A combination of the two methods is the hp-method, where an increase in the

order p is combined with a refinement of the mesh.

There are other methods, such as r-adaptivity or remeshing, where the

The latter two methods acknowledge that it makes no sense to make the

solution error smaller than the modeling error. In transition zones, there might

occur effects which cannot be handled with a reduced model (3D → 2D), but

adaptivity means that the model is changed and the constitutive equations

are amended, while d-adaptivity means that the dimension of the problem is

changed.

In particular, an increase in the order of the polynomials should probably

be reserved for regions where the solution is smooth. In regions where the

model does not fit well, it often makes more sense to shrink the elements and

surmount thus the hurdles. In zones were the thickness changes, it is often

much more effective to increase the number of elements than to increase the

degree of the polynomials.

d-adaptivity (dimension), and m-adaptivity (model), [227].

where a switch back to the original 3-D model is necessary. The term mlayout

of a mesh is changed without changing the number of elements,

1.31 Adaptive methods 155

Fig. 1.111. Conventional adaptive refinement of a plate (bilinear elements), two

load cases wind load and gravity load a) system and wind load; b) optimal mesh

for σxx(x) (wind load); c) for σyy(x) (gravity load)

156 1 What are finite elements?

Fig. 1.112. Goal-oriented refinement (duality techniques), optimal meshes for

σxx(x) under wind load; b) σyy(x) under gravity load

Duality techniques

Let us assume that the displacement at a particular point x is to be calculated

to high accuracy. We know that the distance between the exact and

approximate Green’s function is responsible for the error of the FE solution:

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

_

Ω

[G0(y, x) − Gh0

(y, x) ] p(y) dΩy

= a(G0 − Gh0

, u) . (1.446)

So what is needed is a good approximation of the Green’s function. This is

the strategy of duality techniques: We optimize the mesh in such a way that

the error in the Green’s function G0 becomes small. To achieve this goal first

a point load P = 1 is placed at the point x where u(x) is to be calculated, and

the mesh is refined for this load case by a conventional adaptive procedure.

The solution of this load case is the approximate Green’s function Gh0

for

a)

1.31 Adaptive methods 157

Fig. 1.113. A very uneven mesh but Maxwell’s theorem implies that the displacements

are the same, δh

1,2 = δh

2,1

.

u(x). Then in the second step the original load case is solved on this optimized

mesh. And because the mesh refinement in the first step has reduced the error

G0 − Gh0

(which is an intrinsic property of the mesh) the displacement error

e(x) in the original load case at the point x should be markedly reduced. This

is the basic idea, see Fig. 1.110, Fig. 1.111 and Fig. 1.112.

But in doing so we also incorporate information about the error in the

original problem, u−uh, or—as it is called—the error in the primal problem.

That is the mesh is refined where the error u−uh is large and where the error

G0 − Gh0

is large. Both errors steer the refinement process.

To understand this strategy recall that influence functions are based on

Betti’s principle, W1,2 = W2,1. In Fig. 1.113 a single force P1 = 1 is applied at

the point x1. We want to calculate the horizontal displacement at the point

x2 caused by this force

ux(x2) =

_

Ω

G0(y, x2) • p(y) dΩy = P1 · Gx0

(x1, x2) . (1.447)

Gx0

is the horizontal component of the displacement field G0 = [Gx0

,Gy

0]T ; in

the following we write simply G0 and also u instead of ux.

The Green’s function takes only samples where the load p is applied and if

the load happens to be a point load then it samples at only one point! Hence

to minimize the error in the displacement at the point x2

e(x2) = P1 · [G0(x1, x2) − Gh0

(x1, x2)] = P1 · [δ1,2 − δh

1,2] (1.448)

we have to minimize—so it seems—only the error [. . .] in the Green’s function

at the point x1 where the load is applied.

158 1 What are finite elements?

Now the Green’s function G0(y, x2) takes off from the point x2 where the

Dirac delta, the single force P2 = 1 in Fig. 1.113, presses against the plate

to generate the influence function for u(x2). So would it not make sense to

refine the mesh near this point too so that the effect felt at the distant point

x1 becomes more accurate? The answer is yes! The two points (or rather the

displacements at the two points) are like twins, are adjoint, because according

to Maxwell’s theorem (see p. 550) we have

δh

1,2 = δh

2,1 . (1.449)

Assume we turn the whole situation around: the point load acts at the center

of the “Gargantuan” element Ω2 and we inquire about the displacement u(x1)

at the point x1. According to the previous logic the result should not be too

good because we failed to refine the neighborhood of the point load P2. But

then δh

1,2 = δh

2,1 and so if the second problem is poorly solved then also the

first!

Hence to carry the information from a point A (where the action is) to a

point B where the observer is standing the neighborhood of both points must

be refined. This is an inherent requirement in self-adjoint problems.

But note that the refinement at the two points must not be of the same

order, see Fig. 1.114. It depends on the strength of the singularity at A and

at B, that is what the action is at A (a point force, a single moment, a bend,

a dislocation, a twist) and what we measure at B (the deflection, the slope,

the bending moment, etc.).

If the load is not a point load but spread over a region A then the integration

∂i u(xB) =

_

A

Gi(y, xB) • p(y) dΩy (1.450)

will lower the order of the singularity and then only the refinement at the

point B will be noticeable. (∂i u is short for u,σ, . . .).

The incorporation of the defect p − ph, the error in the primal problem,

is easy, because due to the Galerkin orthogonality (1.446) is equivalent to13

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

_

Ω

[G0(y, x) − Gh0

(y, x) ] • (p(y) − ph(y)) dΩy

= a(G0 − Gh0

,u − uh) . (1.451)

We essentially have added a zero to (1.446), because the error in the Green’s

function is orthogonal to each vh ∈ Vh and therefore also to uh

0 = ph(G0[x]) − ph(Gh0

[x]) = a(G0[x] − Gh0

[x],uh) . (1.452)

13 Our tacit assumption is that the load case ph consists of volume forces ph only.

We shall make this assumption whenever appropriate to simplify the notation.

1.31 Adaptive methods 159

Fig. 1.114. Adaptive refinement a) standard refinement ηp ≤ εTOL b) goal oriented

refinement ηp × ηG ≤ εTOL

160 1 What are finite elements?

After applying the Cauchy-Schwarz inequality to (1.451) we obtain the fundamental

estimate

|u(x) − uh(x)| ≤ ||G0[x] − Gh0

[x]||E ||u − uh||E (1.453)

where ||.||E is the standard energy norm, ||u||2

E := a(u,u).

This result is motivation to minimize the error in the Green’s function and

the error in the solution uh simultaneously. That is, at each refinement step

an error indicator ηeG

for the Green’s function and an error indicator ηe

p for

the original problem is calculated on each element and the combined error

indicator on each element is then ηe = ηeG

· ηe

p.

This last step is based on

|u(x) − uh(x)| = |

_

e

a(G0 − Gh0

,u − uh)Ωe

|

≤

_

e

|a(G0 − Gh0

,u − uh)Ωe

|

≤

_

e

a(G0 − Gh0

,G0 − Gh0

)1/2

_ _ Ωe_

ηeG

· a(u − uh,u − uh)1/2

_ _ Ωe_

ηep

(1.454)

so that the sum of the local errors ηe = ηeG

· ηe

p provides an upper bound for

the error

|u(x) − uh(x)| ≤

_

e

ηeG

· ηe

p . (1.455)

The energy norms of the fields eG = G0−Gh0

and eu = u−uh are calculated

by measuring the eigenwork done by the element residual forces and jump

terms on the element edges as on page 149, so that for example

ηp = a(eu, eu)1/2

Ωe

= r(eu)1/2 . (1.456)

Because this technique consists in applying a quantity (a Dirac delta)

conjugate to the original value, we speak of duality techniques. In particular

the method is well adapted to problems where the focus is on one or two

values that are approximated to a high degree of accuracy [22], [67], [132].

Which is why the technique is also called goal-oriented recovery.

The interesting point about the duality technique is that the error in the

Green’s function serves as a weight to the error in the primal problem. The

consequences of this can be seen in Fig. 1.114. The first mesh, the mesh in

Fig. 1.114 a, is the result of a standard adaptive refinement

ηe

p

≤ εTOL . (1.457)

To push the error below the preset error margin the program has to refine all

those elements—in practice only the first, say, 30%—where the error ηe

p of the

primal problem exceeds this margin. It has no other choice.

1.31 Adaptive methods 161

The mesh in Fig. 1.114 b is based on weighting the primal error ηp with

the error ηG of the Green’s function elementwise

ηeG

· ηe

p

≤ εTOL . (1.458)

In most parts of the mesh the error ηeG

of the numerical Green’s function

is so low that this inequality is automatically satisfied. That is many of the

refinements in Fig. 1.114 a are not necessary if we are only interested in the

stress σyy at x.

For a more detailed analysis of duality techniques and the concept of generalized

Green’s functions see Sect. 7.4, p. 519.

Remark 1.17. Clearly the duality technique is based on Tottenham’s equation.

But historically the duality technique was developed independently, see [19]

and [22]. Only latter it was discovered that the fundamental equation had

already been published by Tottenham, [67]. Probably there are also other

precursors as, for example, the L∗L-method of Kato-Fujita where to obtain

bounds for a solution at a point a similar technique is applied [174].

Nonlinear problems

Although the duality technique is motivated by Betti’s theorem, it can be

applied successfully to nonlinear problems as well, [22], [67], [200]. Near an

equilibrium point a nonlinear structure essentially exhibits a linear behavior

with regard to load increments and this suffices to establish a connection

between the error in the functional, J(u)−J(uh), and the error in the Green’s

function, z − zh, at the linearization point.

By doing a “Taylor expansion” of Green’s first identity it follows, see Sect.

4.21 p. 403, that the displacement increment uΔ and the load increment pΔ

satisfy—to first order—the variational statement

aT (u;uΔ, v) = (pΔ, v) v ∈ V (1.459)

where aT is the Gateaux derivative of the strain energy product, that is the

gradient of a(u, v) at u in the direction of uΔ.

The Newton-Raphson algorithm for the solution of the nonlinear system

of equations k(u) = f is based on this expansion

KT (ui)ui+1 = f − k(ui) (1.460)

where the tangential stiffness matrix KT corresponds to aT .

If the Gateaux derivative aT is linear in the second and third argument

the associated Euler equation

LΔ uΔ = pΔ (1.461)

is linear which means that near an equilibrium point the response of the

structure to load increments is linear and so for each functional J(uΔ) there

162 1 What are finite elements?

is a (generalized) Green’s function z such that the effect of a load increment

on any functional value of the displacement increment can be predicted as in

the linear theory, J(uΔ) = (pΔ, z).

Now the key point is—we have discussed this before in Sect. 1.27—that

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

J(u) − J(uh) =

1

2ρ(uh)(z − zh) +

1

2ρ

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

h (1.462)

where the first term

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

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

ρ

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

,u − uh) (1.464)

h and R(3)

h

So the error z−zh in the Green’s function at the linearization point is one

part of the error J(u) − J(uh) or else: if the error z − zh can be made small

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

small—at least that is what we hope. This is the idea; see also the Remark

on p. 124.

To repeat: the Green’s function z at the linearization point can not be

used to calculate

J(u) =

_

Ω

z(y, x) • p(y) dΩy . (1.465)

This does not work. But the study of the error z − zh allows to adapt the

mesh so that eventually the error J(u) − J(uh) will become small.

Hence basically what is done is, see Box 1.1, that at the actual equilibrium

point the tangential stiffness matrix is solved for the generalized Green’s

function,

KT (u) z = j (1.466)

corresponding to

aT (uh; zh,ϕi) = J

_(uh;ϕi) ϕi

∈ Vh (1.467)

and the error in the primal problem is weighted with the error in this dual

problem. For an application see the example in Sect. 4.21, p. 412 and for a more

detailed analysis of goal-oriented refinement applied to nonlinear problems see

Sect. 7.5 p. 526.

remainder.

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

1.31 Adaptive methods 163

Start with an initial mesh Tk and let k = 0

(A) For t = 1 to M (loop over all load increments M)

λ = t ・ 1

M

actual load parameter

1. Newton-Raphson: Let i = 0, uΔk = 0, u(0)

k = uk−1 and solve the equations

KT (u(i)

k ) u(i+1)

Δk

= λ ・ f − k(u(i)

k )

u(i+1)

k = u(i+1)

Δk

+ u(i)

k

Let i = i + 1 and repeat iteration till convergence is achieved.

2. Calculate the error indicators of this primal problem η(p)

e .

3. Formulate and solve the dual problem at the actual equilibrium point uk:

KT (uk) zk = jk

4. Calculate the error indicators of the dual problem η(z)

e .

5. Mesh refinement:

• Determine the weighted error indicator

ηe = η(p)

e

・ η(z)

e .

• Calculate the error estimator

J(e) ≈ η =

_

Ωe

ηe .

• IF |η| ≤ TOL (global error margin) → t = t + 1, GOTO (A).

• IF ηe > TOLe (local error margin) → refine element Ωe.

• Generate a new mesh Tk+1, transfer data, let k = k + 1.

• GOTO 1.

Box 1.1: Goal oriented refinement of nonlinear problems

Model adaptivity

Often the FE model of a structure is based on simplifying assumptions and

so we would like to have estimates for the modeling error and also indicators

which could steer an adaptive process by which the model can be updated if

necessary14. The equation

J(uc) − J(u) = −d(uc,G) _ −d(u,G) , (1.468)

14 ’

One should make a model as simple as possible but not simpler’ (A. Einstein).

164 1 What are finite elements?

Fig. 1.115. Effect of the modeling error on the bending moment M(x) at the center

of the first span a) original design b) simplified model c) M(x) d) bending moment

of the influence function G2 for M(x) e) elements

which allows to asses changes in any output functional J(u) due to changes in

the stiffness of a structure serves just this purpose. That is Green’s functions

and the d-scalar product allow to estimate how sensitive the results in a part

ΩA of a structure are to the choice of the elastic parameters in a part ΩB of

the structure.

An introductory problem, the bridge in Fig. 1.115 a, may illustrate the

basic idea. The simplified FE model consists of three beam elements with a

uniform stiffness EI. The goal is the evaluation of the bending moment M(x)

at the center of the first span. The error in M(x) is approximately

Mex(x) −M(x) _ −

_ l

0

ΔEI

EI

MM2

EI

dy = −(0.4 − 1.4 + 0.2) (1.469)

1.31 Adaptive methods 165

Fig. 1.116. Change of Young’s modulus, E = 1 ・ 107 → 1 ・ 105 kN/m2 a) stress

distribution under gravity load, uniform E, b) change in σxx along section A−A if

E changes in element # 181 c) change of σxx in element # 184 if E changes—one

element at a time—in section A − A

166 1 What are finite elements?

Fig. 1.117. Shear wall, 10 m × 3 m. Change of Young’s modulus, E = 1・107 → 1・105

kN/m2 a) stress distribution under gravity load, uniform E, b) change in σyy in

element # 112 if E changes elementwise in section A−A c) change of σxx in element

# 12 if E changes elementwise in section B − B

where the three numbers are the contributions of the three spans to the error.

Because the error is largest in the second span in the next step this span

should be subdivided into, say, three separate elements with more realistic

values of EI and then the analysis should be repeated.

Of course the modeling error

ηm =

_

e

|d(uc,G)Ωe

| _

_

e

|d(u,G)Ωe

| (1.470)

1.31 Adaptive methods 167

should be combined with the discretization errors ηp and ηG of the primal and

η = ηp · ηG + ηm (1.471)

may serve as an engineering error indicator.

In Figure 1.116 this technique was applied to a cantilever plate. In the

element # 181, which gets stressed the most, the modulus of elasticity was

reduced from E = 1 · 107 to E = 1 · 105 kN/m2. Plotted in Fig. 1.116 b are

the true change of the horizontal stress σxx in cross section A − A and the

change as predicted by the d-scalar product

σc

xx

− σxx _ −d(Gh1

,uh) (1.472)

where both fields, Gh1

and uh, were taken from the original unmodified model.

Next the modulus of elasticity was changed in cross-section A − A—one

element at a time—and the change in the stress σxx at the center of element

# 184 was studied. The effects on σxx are plotted in Fig. 1.116 c.

Similar modifications were done to the plate in Fig. 1.117 and also the

results are similar.

It is clearly visible (1) that the influence of local changes, E → E +ΔE,

on the stress field is limited and that (2) in most of the cases the d-scalar

product d(Gh1

,uh) with both fields taken from the unmodified structure is

capable of predicting the effects of the change.

Nonlinear problems

Model adaptivity can be applied to nonlinear problems as well because the

change in any output functional can be approximately predicted by applying

Equ. 1.375 in Sect. 1.27, p. 123.

Fundamental solutions

The formula

uh(x) =

_

Ω

Gh0

(y, x) p(y) dΩy (1.473)

suggests, that all that is needed is a good approximation to G0(y, x). To ease

the burden for the FE program, to approximate the function G0 it can be

split into a fundamental solution g0 and a regular part uR,

G0(y, x) = g0(y, x) + uR(y, x) (1.474)

and then only the regular part uR needs to be approximated with finite elements.

The regular part uR is a homogeneous solution of the differential

dual problem so that—eventually after some proper scaling—the expression

168 1 What are finite elements?

Fig. 1.118. Improved goal-oriented refinement with fundamental solutions, optimal

meshes for a) σxx(x) under wind load; b) σyy(x) under gravity load

equation, and its boundary values are such that g0 + uR satisfies the same

boundary conditions as G0.

In the following, this split will be exemplified with the kernels for the

stresses σij in a plate:

Gij

1 (y, x) = uij

R(y, x) + gij

1 (y, x) . (1.475)

The three fundamental solutions form a matrix

D(y, x) = [g11

1 (y, x), g12

1 (y, x), g22

1 (y, x)] (1.476)

with the elements [116]

Dkij(y, x) =

1

r

{(1 − 2 ν)[δkir,j +δkjr,i −δij r,k ] + 2r,i r,j r,k } 1

4π(1 − ν) .

(1.477)

The index k is the vector component, and the indices ij ≡ {11, 12, 22} denote

the stresses. The traction vectors of these fields have the components (k)

1.31 Adaptive methods 169

Skij(y, x) =

2μ

r2

{2 ∂r

∂n

[(1 − 2 ν)δij r,k +ν(δikr,j +δjkr,i ) − 4r,i r,j r,k ]

+ 2ν (nir,j r,k +nj r,i r,k ) + (1 − 2 ν)(2nkr,i r,j +njδik

+ niδjk) − (1 − 4 ν)nkδij} 1

4π(1 − ν) . (1.478)

If the boundary conditions are assumed to be homogeneous on Dirichlet (ΓD)

and Neumann (ΓN) boundaries, the regular solution must satisfy the equations

Gij

1 = uij

R + gij

1 = 0 on ΓD, (1.479)

t(Gij

1 ) = tij

R + tij

1 = 0 on ΓN . (1.480)

The approximate regular part is the solution of the variational problem

a(uij

R,h,ϕi) = −

_

ΓN

tij

1

•ϕi ds ϕi

∈ Vh . (1.481)

The improved Green’s function is

Gij+

1,h (y, x) = uij

R,h(y, x) + gij

1 (y, x) , (1.482)

and the stresses are

σ+

ij,h(x) =

_

Ω

Gij+

1,h (y, x) • p(y) dΩy . (1.483)

Note that the error in the Green’s function is replaced in (1.453) by the error

in the regular solution uR

|u(x) − uh(x)| ≤ || uR[x] − uh

R[x]||E ||u − uh||E , (1.484)

which will speed up the convergence, because at points not too close to the

boundary the error in the regular part will be much smaller than the error

in the Green’s function. As before two error indicators ηR (indicating the

error in the regular part) and ηp are calculated and the combination of both

indicators, η = ηR · ηp, steers the adaptive refinement.

Although there is one subtle difference between the new approach and the

previous technique. In the new approach we actually construct for, say, the

stress σxx at a point x an approximate Green’s functions G+

1,h and we calculate

the stress with this function, see (1.483). While in the previous method we

simply calculate the FE stress σh

xx(x) =

_

i σxx(ϕi)(x) ui (after the mesh

has been improved) knowing that in so doing we inherently form the scalar

product between the improved kernel Gh1

and the applied load.

Comparative study

The plate in Fig. 1.111 on p. 155 was analyzed with all three methods

170 1 What are finite elements?

• conventional energy-norm refinement

• goal oriented refinement

• goal oriented refinement based on fundamental solutions

The plate was subject to a wind load and gravity load. The aim of the refinement

was to predict the horizontal stress σxx at the point x = (4.8, 2.1)

under wind load and the vertical stress σyy at the same point under gravity

load to high accuracy. Hence each load case required a different refined mesh.

The “exact” values in Tables 1.6 and 1.7 are the results of a BE analysis.

As is typical for goal-oriented methods the adaptive refinement concentrates

on the neighborhood of the point x; see Fig. 1.112. While if the Green’s

function is split into a fundamental solution and a regular part then the refinement

spares the neighborhood of x because the fundamental solution is

already the quasi-optimal solution near this point (see Fig. 1.118), and all

effort goes into minimizing the energy error of the regular part, see (1.484).

Table 1.6. Gravity load, σyy, exact = −87.126 kN/m2

d.o.f. Energy norm d.o.f. Goal-or. d.o.f. Fund. sol.

350 −73.9557 350 −73.9557 350 −85.1865

558 −74.2296 528 −84.0145 576 −85.8080

900 −74.9279 823 −87.7160 870 −86.0676

1462 −74.6239 1206 −85.5208 1235 −86.2783

2325 −85.0523 1619 −85.0156 1694 −86.4994

Table 1.7. Wind load, σxx, exact = −43.634 kN/m2

d.o.f. Energy norm d.o.f. Goal-or. d.o.f. Fund. sol.

350 −45.2505 350 −45.2505 350 −43.7053

538 −45.2789 518 −42.5768 537 −43.7753

858 −45.0542 716 −43.4873 766 −43.6761

1370 −45.0911 982 −44.0888 1065 −43.6685

2079 −45.0155 1254 −43.8960 1378 −43.6596

– – 1565 −43.6094 – –

Changes in the elastic parameters

Even if the material is not homogeneous, the Green’s function can still be

split into a fundamental solution and a regular part. It is merely necessary

that additional forces must be applied at the interface of the different zones.

1.31 Adaptive methods 171

Fig. 1.119. The Green’s function is split into a fundamental solution and a regular

solution, which consists of two parts

These forces account for the discontinuities in the material parameters. The

1-D problem of a simple bar will suffice to illustrate the technique.

The bar in Fig. 1.119 consists of two elements with different cross sections

A1 and A2. The Green’s function for the longitudinal displacement at the

center of the second element is to be calculated. The displacement g0 in Fig.

1.119 a is an appropriate fundamental solution because it corresponds to the

application of a single force P = 1 at the center of the second element. Now the

Green’s function has to have (i) zero displacement at x = 0, (ii) a zero normal

force N(l) = 0 at the free end of the bar and (iii) the jump in the normal

force must be zero at the transition point x1 between the two zones. Note

that the fundamental solution does not and cannot satisfy this last condition

because the slope g_

0 of the fundamental solution is continuous at x1. Its value

is 1/(2EA2), so

N+(x1) − N

−(x1) = EA2 g

_

0(x1) − EA1 g

_

0(x1)

=

1

2

(1 − A1

A2

) = −A1 − A2

2 . (1.485)

172 1 What are finite elements?

Because A1 > A2, the net effect is that a force (A1−A2)/(2A2) pulls the node

x1 to the left:

1

2

A1

A2

← • → 1

2 . (1.486)

(Multiply the right-hand side by A2/A2.) Now to correct this defect and

enforce N(l) = 0 at the end of the bar (currently N(l) = −1/2) the regular

solution must solve the load case in Fig. 1.119 c, and in addition a

constant term 1/EA2 must be added to correct the nonzero displacement,

g0(0, xc) = −1/EA2, of the fundamental solution at x = 0:

G0 = g0 + u(Fig 1.119 c) + 1/EA2 . (1.487)

Навигация

• Главная
• Книги
• Статьи
• Обратная связь
• Поиск