3.7 Hierarchical elements

Back

In the p-method the element size is kept fixed while the degree p of the polynomial

expansion is increased to order 2, 3, . . .. If the added shape functions

are orthogonal in the sense of the strain energy product to the previous set of

functions the new stiffness matrix is simply obtained by amending the previous

matrix. Such elements are called hierarchical elements.

Consider a bar element [−1, 1]. To its two linear shape functions

ϕ1(x) =

1

2

(1 − x) ϕ2(x) =

1

2

(1 + x) (3.156)

we add three shape functions which vanish at the end points of the element

[−1, +1] (see Fig. 3.33)

ϕ3(x) =

1 √

6

[−1 +

1

2

(−1 + 3x2)] ϕ4(x) =

1 √

10

[−x +

1

2

(−3 x + 5x3)]

ϕ5(x) =

1 √

14

[

1

2

(1 − 3 x2) +

1

8

(3 − 30 x2 + 35x4)] . (3.157)

Fig. 3.33. Bar ele-

1 to +1. The linear

shape functions

hierarchical functions

310 3 Frames

Fig. 3.34. Approximation of the Greens function for the support reaction at the left

end with linear elements and hierarchical elements (increasing order of p = 2. . . 10)

These new functions are orthogonal to the first two functions ϕ1 and ϕ2

k1j = a(ϕ1, ϕj) = k2j = a(ϕ2, ϕj) = 0 j = 3, 4, 5 (3.158)

and they are also mutually orthogonal

kij = a(ϕi, ϕj) = δij · EA i, j = 3, 4, 5 (3.159)

so that the amended matrix is simply, [26] p. 252,

K = EA

2

⎢⎢⎢⎢⎣

1 −1 0 0 0

−1 1 0 0 0

0 0 2 0 0

0 0 0 2 0

0 0 0 0 2

⎥⎥⎥⎥⎦

. (3.160)

For an application we consider a bar [−1, +1] fixed at the left end and

stretched by uniform forces p = 1

− EAu

__(x) = 1 u(−1) = 0 N(1) = EAu

_(1) = 0 . (3.161)

Because the left end is fixed u1 = 0 but the other four ui are unknown. The

equivalent nodal forces are

3.7 Hierarchical elements 311

Fig. 3.35. In the porder

p of the shape

functions

f2 =

_ +1

−1

1 · ϕ2(x) dx = 1 f3 =

_ +1

−1

1 · ϕ3(x) dx = −

%

2

3

f4 =

_ +1

−1

1 · ϕ4(x) dx = 0 f5 =

_ +1

−1

1 · ϕ5(x) dx = 0. (3.162)

The solution of the system

K = EA

2

⎢⎢⎣

1 0 0 0

0 2 0 0

0 0 2 0

0 0 0 2

⎥⎥⎦

⎢⎢⎣

u2

u3

u4

u5

⎥⎥⎦

=

⎢⎢⎣

f2

f3

f4

f5

⎥⎥⎦

(3.163)

is (we let EA = 1)

u = {2,−

_

2/3, 0, 0}T (3.164)

so that

uh(x) = 2 · ϕ2(x) −

%

2

3

· ϕ3(x) =

3

2

+ x − x2

2

(3.165)

which is the exact solution.

Note that the sum of the equivalent nodal forces

_5

i=1

fi =

_ +1

−1

p · (ϕ1 + ϕ2 + ϕ3 + ϕ4 + ϕ5)

_ _ _

            =1

dx = 1.1835

            = 2.0 =

_ +1

−1

pdx =

_ +1

−1

p · (ϕ1 + ϕ2)

_ _ _

=1

dx = f1 + f2 (3.166)

is not the sum of the applied load (2.0) because the extended set of shape

functions does not form a partition of unity

method the support

reaction oscillates

considerably with the

312 3 Frames

_5

i=1

ϕi(x) = 0.62 − 0.79 x − 0.79 x2 + 0.79 x3 + 1.17 x4 (3.167)

—simply said the sum is not 1 at each point x. So in checking the global equilibrium

condition in the p-method we must restrict the count to the original

fi.

Point loads and hierarchical elements

The p-method will improve the accuracy considerably if the exact solution is

smooth but it can run into difficulties in the presence of point loads, [238] p.

196.

The model problem is a bar [−1, +1] that is fixed at both ends and is

subjected to a horizontal point load P = 1 at the quarter point x = −0.5. We

study the support reaction at the left end of the bar, at x = −1. The exact

Green’s function of the support reaction is the straight line G1(x) = −0.5 x+

0.5, dropping from + 1 at the left end to zero at the right end of the bar. The

triangle in Fig. 3.34 is the best approximation with 10 linear elements. For our

purposes it is perfect because the FE Green’s function Gh1

is exact at x = −0.5

so that the linear model gives the correct answer, Nh(0) = N(0) = 0.75.

If we solve the same problem with the p-method — just one element but

different orders p of shape functions (Pi = Lagrange polynomial)

ϕi(x) = Pi+1(x) − Pi−1(x) √

2

_

2(i + 1) − 1

i = 1, 2, 3, . . . p (3.168)

then the support reaction oscillates considerably (see Fig. 3.35) because in

the p-method the approximate Green’s functions Gh1

tend to wobble; see Fig.

3.34. Note that the exact solution u(x) does not lie in the trial space Vh of

the p-method because all ϕi(x) have continuous first-order derivatives and

also the Green’s function G1(x) is not contained in Vh—linear functions are

excluded because of the boundary conditions—so we must expect an error in

the support reaction.

If we would place a node where the single force is applied the solution

would lie in Vh so the problem could have readily been resolved — in this case

— but evidently care must be taken in the presence of point loads.

Remark 3.2. We add some details. The stiffness matrix in the foregoing problem

is EA × I (the unit matrix), the equivalent nodal forces for Gh1

are

fi = N(ϕi)(−1) = EAϕ_

i(−1) so that the nodal values ui of the Green’s

function Gh1

are the derivatives ϕ_

i of the shape functions at x = −1 and thus

Gh1

(x) =

_p

i=1

ϕ

_

i(−1) ϕi(x) . (3.169)

3.8 Sensitivity analysis 313

Fig. 3.36. Change

in the stiffness of a

spring

If the load p is uniformly distributed all is well. In the influence integral (Gh1

, p)

only the integral of the first term of the series (3.169) is not zero and because

the slope of ϕ1 at x = −1 and the integral of ϕ1 are well tuned the result

Nh(−1) =

_ +1

−1

Gh1

(x)pdx = p ϕ

_

1(−1)

_ +1

−1

ϕ1(x)dx

= p · 0.816 · 1.22474 = p · 1.0 (3.170)

is exact—the total load is p · 2.0.