2.2 Structural analysis with boundary elements

Back

Some examples may serve to highlight the application of BE methods in structural

analysis, and the potentials of the method.

Shear wall

The shear wall in Fig. 2.6 carries its own weight plus line loads coming from

different structural elements that are connected to the wall. The wall rests on

a series of narrow supports, which in FE analysis would be modeled as point

supports. Here they are modeled as short boundary elements.

The outer boundary plus the edges of the openings were subdivided into

a total of 123 quadratic boundary elements. Some minor inconsistencies arise

at the supports, because at the transition point between the free edge and a

support, the physical tractions are discontinuous, while C0 shape functions

cannot model this behavior—but these are very minor details.

The unknowns are the boundary displacements ux, uy along the free edges,

while at the supports, some (at a roller support) or all of the tractions tx, ty

are unknown. After the support reactions and the shape of the deformed

structure are determined by solving the system Hu = Gt + d, the influence

function (2.15) makes it possible to calculate the displacements and stresses

in the interior.

248 2 What are boundary elements?

Fig. 2.7. Wall on point supports: a) system and loading, b) principal stresses

Two-span wall

The two-span beam in Fig. 2.7 was modeled as a plate. To compare the results

with a beam solution, the plate was placed on point supports. Theoretically

the support reactions should be zero, because the influence of a support reaction

P on the displacement of the support itself is −P ln r = −P ln 0 = ∞,

since the distance is zero (r = 0). This implies that the support reaction must

be zero, P = 0, since in a set of equations such as

_

∞ b

c d

__

P

x

_

=

_

e

f

_

(2.20)

P must be zero.

But in reality the BE program replaces the function call ln r for values of

r < 10−3, (= one millimeter) with ln 10−3 = −6.9077 to avoid overflow. This

is equivalent to replacing the point support with a very short line support,

which suffices to make P equal to the beam solution.

Slabs

Engineering plate-bending problems are much more complex than the biharmonic

problems usually discussed in the mathematical literature. Add to this

that the correct modeling of a complex floor plate depends on so many parameters,

and so many assumptions must be made by the analyst, that these

assumptions tend to have much more influence on the results than, say, the

order of the boundary elements or other mathematical “subtleties”. Hence

programming the biharmonic equation KΔΔw = p for the analysis of slabs

is a real challenge, but the results are rewarding; see Fig. 2.8.

2.2 Structural analysis with boundary elements 249

Fig. 2.8. Slab: a) system and 274 boundary elements, b) principal moments under

gravity load, c) 3-D plot of the deflection surface

250 2 What are boundary elements?

The choice of boundary elements for interpolating the boundary values

w,wn,mn, vn is dictated by the regularity assumptions made in the derivation

of the integral equations [116]. The deflection w must be C1 on any smooth

part of the boundary, while for the other functions C0-continuity suffices.

Support conditions are usually of mixed type, not jut w = mn = 0 (hinged)

or w = wn = 0 (clamped). Most supports are modeled as elastic supports, so

that a hinged support becomes cw + vn = mn = 0 and a clamped support

becomes cw + vn = cϕwn + mn = 0, where c and cϕ are elastic constants.

Add to this that at corner points, because of the continuity of the gradients

∇w, ∇∇w etc., boundary conditions on both sides of a corner point cannot

be formulated independently of each other. If in addition the thickness of

a plate changes at the corner point, then things can become really complicated.

In a BE program, all the modeling is done on the boundary, and the

boundary of a Kirchhoff plate is a very thin layer. But luckily the biharmonic

equation—being of elliptic type—is a very patient equation, so that for engineering

purposes the accuracy attainable with BE methods is more than

sufficient, and at least on the same level as FE results.

Walls and T beams

Internal supports are subdivided into boundary elements (or rather line elements)

to model the distribution of the support reaction s with piecewise

linear functions; see Fig. 2.9. The nodal values are so determined that the

deflection w of the plate is zero at the nodes of the wall (or that cw+s = 0 if

the walls are elastic). As in FE analysis, it can be assumed that the calculated

support reactions are relatively accurate.

The support reactions s of T beams are determined such that the deflection

of the plate and the T beams are the same at the nodes of the T beams. The

T beams are modeled with finite elements, and the reduced stiffness matrices

K are inverted to provide flexibility matrices F = K(−1), which better fit a

boundary element scheme, because F s = w.

Columns

The BE method can assess the magnitude of the bending moments near

columns very accurately, because the correct singularity is built into a BE

program; see Fig. 2.10. If p = P/Ωc is the pressure in the contact zone Ωc,

where P is the support reaction, the influence of the support reaction P on

the bending moment mxx(x) is the integral

_

Ωc

[−K(g0,x1x1 +ν g0,x2x2 )] pdΩy, g0,x1x1 = ∂2g0

∂2x1

(y, x) , (2.21)

and these kernel functions (g0(y, x)),xixj are part of the BE code, and must

not be approximated with piecewise polynomials as in the FE method. The

2.2 Structural analysis with boundary elements 251

Fig. 2.9. Slab: a) principal moments, and b) contour plot of the deflection

surface integral in (2.21) also ensures that the bending moments are automatically

rounded out; see Fig. 2.11.

The same can be said about the shear forces. The contribution of the

support reaction in the column is the integral

_

Ωc

[−K(g0,x1x1x1 +g0,x2x2x1 )] pdΩy , (2.22)

where

252 2 What are boundary elements?

Fig. 2.10. Slab on

point supports

(g0(y, x)),x1x1x1 =

2

8πK r

_

−2 r,22

r,1 −r,31

−r,1 r,22

_

, (2.23)

(g0(y, x)),x2x2x1 =

1

8πK r

_

−r,32

−2 r,22

r,1 −r,31

_

. (2.24)

The resolution attainable with these influence functions depends solely on how

tightly the observation points x are packed; see Fig. 2.12.

The internal actions are smooth continuous functions, because each individual

value is calculated by evaluating the proper influence function. Theoretically

there is a loss of accuracy close to the boundary, but for engineering

purposes this loss is negligible; see Fig. 2.13.

115.74 kN

294.88 kN

343.92 kN

181.14 kN

180.78 kN

590.37 kN

739.76 kN

345.09 kN

127.58 kN

304.22 kN

123.73 kN

2.2 Structural analysis with boundary elements 253

Fig. 2.11. Slab 40 m × 30 m on a grid of columns 5 m × 5 m, p = 10 kN/m2 ;

comparison of bending moments with building code

Fig. 2.12. Distribution of the shear forces qx between the columns. The resolution

can be made so fine that the variation of qx across the column head is detailed out

0.10 0.61 1.11 1.62 2.12 2.63 3.13 3.64 4.14 4.65 5.15 5.66 6.16 6.67 7.17 7.68 8.18 8.69 9.19 9.7010.00

-49.63

-44.12

-37.80

-31.49

-25.17

-18.86

-12.54

-6.22

0.09

6.41

12.73

- 36,72

- 49,63 k N m / m

- 40,00

- 28,00

13,00 12,73

5,00 5,00

0,5 x 0,5

x

y

code

(code)

(code)

(BE)

(BE)

0.10 0.61 1.11 1.62 2.12 2.63 3.13 3.64 4.14 4.65 5.15 5.66 6.16 6.67 7.17 7.68 8.18 8.69 9.19 9.7010.00

-137.29

-109.83

-82.37

-54.92

-27.46

27.46

54.92

82.37

109.83

137.29

0.01 0.05 0.09 0.13 0.17 0.21 0.25 0.29 0.33 0.37 0.41 0.45 0.49 0.53 0.57 0.61 0.65 0.69 0.74 0.780.80

-136.57

-109.25

-81.94

-54.63

-27.31

0.00

27.31

54.63

81.94

109.25

136.57

254 2 What are boundary elements?

gravity load 5 kN/m2

Fig. 2.13. Slab with a large opening a) system and loads, b) principal moments;

the loss of accuracy close to the boundary is negligible

2.2 Structural analysis with boundary elements 255

Loads

Formally it is of no concern which load is applied and where, because every

load makes its own contribution to the influence function; see Fig. 2.13. The

influence of a force P is the term g0(yP , x) · P + regular parts. If it is a line

load, the influence is a line integral, and if it is a surface load on a patch Ωp,

the influence function is a surface integral:

g0(yP , x) · P

_ _ _

single force

_

ΓL

g0(y, x) l(y) dsy

_ _ _

line load

_

Ωp

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

_ _ _

surface load

.(2.25)

In particular, a BE program is well-qualified to calculate influence functions,

because it must only approximate the regular part of the solution.

Foundation slabs

In the Winkler-model the soil is treated as a grid of springs that can move

independently, and that exert a force cw, where c (kN/m3) is the modulus of

subgrade reaction, so that

KΔΔw + cw = p . (2.26)

The simplest approach to solving this equation is it to actually place the slab

on a grid of springs and determine the unknown spring forces by extending

the collocation equations to the springs. In many cases the accuracy attainable

with this technique is sufficient. Its advantage is that it can be directly

incorporated into an existing BE code with very slight modifications.

A more scientific approach is it to employ the two Fourier–Bessel integrals

g0(y, x) = a2

2πK

_ ∞

0

t

t4 + k

K a4

J0(t

r

a

) dt (2.27)

and

g1(y, x) = a rn

2πK

_ ∞

0

t

t4 + k

K a4

J1(t

r

a

) dt rn = ∇xr •n (2.28)

which are fundamental solutions of (2.26) [126]. Formally the influence function

for the deflection w(x) is nearly identical to the influence function (2.8),

c(x)w(x) =

_

Γ

[ g0(y, x) vν(y) − g0ν(y, x)mν(y)

−vν(g0)(y, x)w(y) + mν(g0)(y, x)wν(y)] dsy

+

_

Ω

g0(y, x) p(y) dΩy +

1

8

_

c ˆK

F(w)(x)

+

_

c

[ g0(yc, x) F(w)(yc) − w(yc) F(g0)(yc, x)] , (2.29)

256 2 What are boundary elements?

except for an extra free term (F/8

_

c ˆK ) where ˆK = a4 c/K.

The Fourier–Bessel functions

Ti,j :=

1

2 π

_ ∞

0

ti

t4 + ˆK

Jj(t

r

a

) dt (2.30)

can be represented in the form of zero-order Kelvin or Thomson functions

kei(x) and ker(x) and their respective derivatives

T1,0 = − 1

2 π

kei( r

a

) T3,0 =

1

2 π

ker( r

a

) , (2.31)

which can be calculated using rapidly convergent expansions [127].

Half-space model

In the year 1885 Boussinesq [49] found the solution for an elastic half-space

loaded at is surface with a vertical force P:

uB := u3(x, y) =

1 + ν

2πE

_

(x3 − y3)2

r3 + 2

1 − ν

r

_

P . (2.32)

If the interface between the foundation slab and the soil is subdivided into

rectangular elements, the soil pressure can be expanded with regard to the n

nodal shape functions giving the expression

wS(x) =

_

i

_

Ω

uB(y, x) ϕi(y) dΩy · pi (2.33)

for the deflection of the soil. The nodal values p i are found by requiring that

the deflection of the slab and the soil be the same at each node:

w(xi) = wS(xi) i = 1, 2, . . . n . (2.34)

By this very simple technique, a coupled analysis of a foundation slab on

top of an elastic half-space can be modeled; see Fig. 2.14. What is striking is

how different the deformation patterns of a foundation slab are, depending on

which model—the Winkler model or the half-space model—is used; see also

Sect. 5.17, p. 469.

If the soil consists of different layers with moduli Ei and Poisson ratios νi,

the total settlement beneath a point force P is the sum of all the individual

contributions si within the single layers:

uΣB

(x) =

_

i

si(x)P . (2.35)

The contribution si = si(Ei, νi, hi) of one layer can be calculated approximately

as follows:

2.2 Structural analysis with boundary elements 257

Fig. 2.14. Foundation plate: a) system and loads, b) deformation according to the

Winkler model, c) deformation according to the half-space model

gravity load g = 7.8 kN/m2

78.0 kN/m

332.0 kN/m

193.8 kN/m 453.8 kN/m

332.3 kN/m

89.1 kN/m

50.7 kN/m

108.9 kN/m

88.5 kN/m

90.5 kN/m

110.9 kN/m

79.3 kN/m

94.9 kN/m

121.9 kN/m

91.9 kN/m

99.8 kN/m

97.7 kN/m

87.3 kN/m 110.3 kN/m

117.6 kN/m

107.6 kN/m

136.6 kN/m

109.2 kN/m

26.6 kN/m

77.4 kN/m 142.9 kN/m

65.3 kN/m

23.3 kN/m

175.1 kN/m

192.7 kN/m

400.3 kN/m

100.1 kN/m

105.0 kN/m

174.5 kN/m

53.5 kN/m

x

y

13.6

15.3

X

Y

Z

X

Y

Z

258 2 What are boundary elements?

Fig. 2.15. Piled raft foundation

si(x) = [uB(x, yi

u) − uB(x, yil

)]|Ei,νi,hi

(2.36)

|yi

u

− yil

| = hi thickness of the layer , (2.37)

where the points yi

u and yil

are the projections of the surface point x onto the

upper and lower cover of the layer i. That is, the contribution of an individual

layer is calculated as if the total half-space had the elastic properties Ei, νi.

Piled raft foundations

In a piled raft foundation the contribution of the piles as well as the raft is

taken into account to transfer the building load. The piles transfer a part of

the loads into deeper and stiffer layers of soil and thereby allow the reduction

To model such foundations it must be possible to predict the soil deformation

if a vertical point force acts in the interior of the elastic half-space.

Mindlin found the following solution for this problem [166]:

uM := u3(x, y) = P

16πE (1 − ν)

_

3 − 4 ν

R3

1

+

8 (1 − ν)2 − (3 − 4 ν)

R2

+

(x3 − y3)2

R3

1

+

(3 − 4 ν)(x3 + y3)2 − 2 y3 x3

R3

2

+

6 y3 x3 (x3 + y3)2

R5

2

_

, (2.38)

where

R1 =

 

x21

+ x22

+ (x3 − y3)2 R2 =

 

x21

+ x22

+ (x3 + y3)2 . (2.39)

The coupled problem of the foundation plate and the elastic half-space can be

stated as follows: on the free surface the tractions must vanish, t = 0, while

settlement and differential settlement in a very economic way; see Fig. 2.15.

2.2 Structural analysis with boundary elements 259

at the interface between the slab and the soil surface the deflection must be

the same, and the vertical stresses must have opposite signs:

u3 =w t3 − pS = 0. (2.40)

Here pS is the soil pressure, which also appears on the right-hand side of the

plate equation

KΔΔw = p − pS . (2.41)

Modeling the slab

The unknowns are the boundary values w and wn of the slab, the soil pressure

distribution pS under the slab and the shear forces s along the piles so that

for example the first integral equation becomes

c(x)w(x) =

_

Γ

[−Vν(g0)w +Mν(g0)∂w

∂ν

] dsy +

_

Ω

g0(p − pS) dΩy

−

_

i

Qi g0[yi] +

_

c

[−w(yc)F(g0)(yc, x)] (2.42)

where the Qi are the pile forces.

Modeling the soil

The soil pressure distribution pS is interpolated with piecewise linear shape

functions

pS(x) =

_

i

pi ϕi(x) , (2.43)

and to simplify the notation it is assumed that the same is done with the

surface load p coming from the building.

When the edge of the slab is divided into boundary elements and the first

and second integral equation (2.42) formulated at the collocation points, the

pertinent system of equations is

Cu = Hu+D(p − pS) + Lq (2.44)

or if the unknowns are placed on the left-hand side and an index B is attached

to denote either that the terms are boundary terms, u = uB, or to indicate

that the collocation points lie on the boundary

CB uB −HB uB +DB pS

− LB q = DBp . (2.45)

The vector q = [Q1,Q2, . . .]T contains the pile forces. The matrices H,D,L

are influence matrices that describe the influence of the nodal values ui, pi, pSi

260 2 What are boundary elements?

and the pile forces Qi at the collocation points xi on the boundary. The matrix

CB is nearly identical to 1/2×I. Only at corner points are the entries different.

The soil pressure pS and frictional stress fsi =: fi along the individual

piles are responsible for the deflection of the surface of the elastic half-space

so that the deformation is a superposition of Boussinesq’s solution uB and

Mindlin’s solution uM:

u3(x) = u(x) =

_

Ω

uB(y, x) pS(y) dΩy +

_

i

_

Γi

uM(y, x) fi(y) dsy (2.46)

where Γi is the surface of the individual piles i = 1, 2, . . . and fi is the frictional

stress acting along pile i.

Let uI denote the deflection of the soil surface at the nodes of the grid,

that is the points at which the soil pressure pS is interpolated. The influence

of the soil pressure pS and the frictional stress fi as expressed by (2.46) can

then be written as

usoil

I = BI pS +MI f . (2.47)

At the grid points the deflection of the soil (2.47) and the deflection of the

slab

uslab

I = HI uB +DI (p − pS) + LI q (2.48)

must be the same, uslab

I

− usoil

I = 0, which gives

HI uB −DI pS + LI q − BI pS

−MIf = −DIp . (2.49)

Modeling the piles

At the interface between the piles and the soil, the pile displacement u and the

vertical soil displacement u3 must be the same (slipping is neglected), which

means that the deformations must be the same,

u3(xi) = uB(xi) + uM(xi) = u(xi) , (2.50)

at the k collocation points along the axis of each of the n piles. At one individual

pile i therefore,

BPi pS +MPi f = ui = K

−1

i Nifi = Fi fi (2.51)

where Fi is the flexibility matrix of the pile and the matrix Ni has the

elements

Nkj = π d

_ li

0

ψk ψj dx d = diameter of the pile . (2.52)

2.2 Structural analysis with boundary elements 261

Fig. 2.16. Study of a piled raft foundation a) settlement, b) pile forces

The ψk are the shape functions of the pile. The vector Ni fi = si represents

the equivalent nodal forces si, due to the frictional stress f. Because each pile

has an influence on all the others it follows that

BP pS + (MP −F N) f = 0 (2.53)

where the diagonal entries in the matrix F are the flexibility matrices Fi

of the individual piles and the subscript P on the Boussinesq and Mindlin

influence matrices indicates that the control points, the collocation points, lie

on the piles.

3.67

4.24

4.52

4.80

2.26

1.98

5.09

1.70

1.41

1.13

0.85

0.57

5.37

387.47 kN

422.19 kN

911.16 kN

1046.00 kN

514.76 kN

479.12 kN

377.23 kN

467.07 kN

425.49 kN

896.96 kN

262 2 What are boundary elements?

Finally at each pile the sum (the integral) of the shear forces fi must be

equal to the resultant pile force Qi which, with an abstract matrix Σ whose

somewhat simplified entries are the lengths li of the single bar elements, can

be stated as

ΣNf = q . (2.54)

If the results are collected, the following system of equations is obtained:

⎡

⎢⎢⎣

CB −HB DB −LB 0

HI −DI − BI LI −MI

0 BP 0 MP −F N

0 0 I −ΣN

⎤

⎥⎥⎦

⎡

⎢⎢⎣

uB

pS

q

f

⎤

⎥⎥⎦

=

⎡

⎢⎢⎣

DB p

−DI p

0

0

⎤

⎥⎥⎦

(slab)

(soil)

(piles)

(sum)

The slab in Fig. 2.16 was analyzed with this technique [117].

Навигация