2.3 Comparison finite elements—boundary elements

Back

In the BE analysis of a plate at each point the influence that all the different

sources have on that particular point must be calculated separately. Typically

the influence a source exerts depends on the distance between the source point

y and the observation point x. This means that it typically fades away like

ln r (in the near term) or like r−1 or r−2, and the influence depends on the

angle ϕ between the two points. A change in the position, (r, ϕ) → (ˆr, ˆϕ),

makes that the influence coefficients change, and thereby also the stresses.

In the FE method the situation is basically the same. It is only that the

FE method plugs all the information on how the influence changes when a

point (r, ϕ) moves to a new location (ˆr, ˆϕ) into the nodal displacements ui

of an element, and interpolates the element displacement field using linear or

quadratic polynomials. Because of this technique the element is autonomous.

The price the FE method pays for this technique is that it is not as sensitive

as the BE method. If for example a cantilever plate is modeled with bilinear

elements and Poisson’s ratio is zero, ν = 0, then the FE bending moment is

the same in any vertical plane of an element, which is impossible. Obviously

the 2 × 4 data cells ui of a bilinear element are not sufficient to register the

change in the influence coefficients if the observer moves from x to ˆx = x+Δx.

But one can also view things differently: an FE program replaces the original

load case p only with a load case ph that is compatible with its limited

means.

Green’s function

The FE solution, uh(x) = (G0[x], p), is based on a mesh-dependent approximation

of the Green’s function. Does this also hold true for BE solutions? To

answer this question let us study the Poisson equation −Δu = p in a domain

2.3 Comparison finite elements—boundary elements 263

Fig. 2.17. Influence function for qx a) BE solution Gxx

3 = gh,xx

3 + uh

R, b) regular

part uh

R = Gxx

3

− gh,xx

3

Ω with boundary conditions u = 0 on ΓD (fixed edge) and ∂u/∂n = 0 on ΓN

(free edge) (Γ = ΓD ∪ΓN). It is assumed that the solution can be represented

by the formula

u(x) =

_

Ω

G0(y, x) p(y) dΩy =

_

ΓD

g0(y, x) t(y) dsy

−

_

ΓN

∂g0

∂ν

(y, x) u(y) dsy +

_

Ω

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

where t = ∂u/∂ν is the slope (traction) on the boundary at the integration

point y with normal vector ν and Green’s function G0. The function g0 in the

second integral representation is a suitable fundamental solution, for example

g0 = −1/(2 π) lnr. If the BE solution

264 2 What are boundary elements?

uh(x) =

_

ΓD

g0(y, x) th(y) dsy −

_

ΓN

∂g0

∂ν

(y, x) uh(y) dsy

+

_

Ω

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

is compared with the FE solution

uh(x) =

_

i

ui ϕi(x) =

_

Ω

Gh0

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

there seems not to be much agreement. But the BE solution (2.56) is just the

formula

uh(x) =

_

Ω

Gh0

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

u(x) =

_

Ω

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

=

_

Ω

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

_

Ω

uR(y, x) p(y) dΩy (2.59)

from which it follows that the boundary integrals in the influence function

(2.55) are just an equivalent expression for the work done by the distributed

load p acting through the regular part uR:

_

Ω

uR(y, x) p(y) dΩy =

_

ΓD

g0(y, x) t(y) dsy −

_

ΓN

∂g0

∂ν

(y, x) u(y) dsy .

(2.60)

Hence it can be assumed that the boundary integrals in the BE solution (2.56)

play the same role,

_

Ω

uh

R(y, x) p(y) dΩy : =

_

ΓD

g0(y, x) th(y) dsy −

_

ΓN

∂g0

∂ν

(y, x) uh(y) dsy ,

(2.61)

and so we arrive at (2.58), which makes the two methods look alike, even

though the FE Green’s function is only mesh dependent, while the BE Green’s

function is also load-case dependent.

The difference between the two methods is that the FE method must

approximate both parts of the Green’s function G0 = g0 + uR while the BE

method must only approximate the regular part uR because the fundamental

solution g0 is part of the code.

in disguise, of course with a different Gh0

. To see this note that the Green’s

function G0 can be split (see Fig. 2.17) into a fundamental solution g0 and a

regular part uR,

2.3 Comparison finite elements—boundary elements 265

Note that the BE method uses an ingenious approach: it does not approximate

uR(y, x)—this would be simply too laborious, because at every point

x it would have to approximate a new function uR(y, x)—but instead substitutes

for the work integral (uR[x], p) the work done by the Cauchy data u

on ΓN and t on ΓD via the conjugate terms of the fundamental solution g0.

That is, the program knows that it suffices to approximate the “static” data

u and t leaving the effects caused by a change in the observation point x to

the fundamental solution g0(y, x). This is the essence of (2.61).

Note also that the boundary values

˜uh(x) = uh(x) + ε(x) ˜th(x) = th(x) + η(x) (2.62)

of uh(x) in (2.56) differ (slightly) from the values under the integral signs,

because the solution (2.56) satisfies the integral equation

c(x) u(x) =

_

ΓD

g0(y, x) th(y) dsy −

_

ΓN

∂g0

∂ν

(y, x) uh(y) dsy

+

_

Ω

g0(y, x) p(y) dΩy x ∈ Γ (2.63)

only at the collocation points xi. On ΓD the left-hand side is zero, u(x) = 0,

while on ΓN the left-hand side is c(x) uh(x) where uh(x) is the same function

as under the second boundary integral.

Accuracy

Can we expect the same tendency as in the FE method, namely that the

accuracy reflects the nature of the Green’s function? To check this let us apply

Betti’s theorem to the pair {uh,G0}, where uh is the BE solution (2.56) of

the Poisson equation −Δu = p with boundary conditions u = 0 on ΓD and

t = 0 on ΓN and where G0 is the Green’s function for u(x) at a point x. Then

_

Ω

G0 pdΩy +

_

ΓN

˜th G0 dsy = uh(x) · 1 +

_

ΓD

t0 ˜uh dsy (2.64)

and

u(x) =

_

Ω

G0 pdΩy +

_

ΓN

G0 t dsy −

_

ΓD

t0 udsy , (2.65)

or

u(x) − uh(x) =

_

ΓN

G0 (t − ˜th) dsy −

_

ΓD

t0 (u − ˜uh) dsy (2.66)

where t0 are the tractions of the Green’s function on Γ. Obviously the BE

solution will be exact if ˜uh = u on ΓD (which is only true at the collocation

points xi), and if ˜th = t on ΓN.

266 2 What are boundary elements?

In view of our previous interpretation it can be stated as well that the

error is attributable to the approximation error in uh

R:

u(x) − uh(x) =

_

Ω

(uR − uh

R)pdΩy . (2.67)

Compare this with the FE error

u(x) − uh(x) =

_

Ω

(G0 − Gh0

)pdΩy , (2.68)

and the difference between the two methods becomes apparent. In the FE

method the accuracy is linked to the nature of the Green’s function, and in

the BE method the accuracy is linked to the character of the regular solution

uR, which clearly deteriorates the closer the point comes to the boundary.

(Hence “regular solution” may be an euphemism.)

This is confirmed by (2.66), because the propagation of the error is governed

by the Green’s function G0 = O(ln r) and the traction t0 = O(1/r) of

the Green’s function. Both functions “live” on the boundary. Seemingly, the

closer the point x comes to the boundary Γ, the more negative the influence

of the boundary error on the results. This agrees with our observations.

The load case ph

In FE methods p is replaced by a work-equivalent load case ph. In BE analysis

the choice of ph is not so simple to explain. But the BE method bears

some resemblance to the force method. The volume potential (G0, p) of a BE

formulation corresponds to the deflection curve w0 of the statically determinate

structure and the boundary potentials are the deflection curves wi of the

redundant forces Xi.

Because the BE solution satisfies the differential equation in any patch

ΩP—“a truck remains a truck”—local equilibrium is satisfied in any patch

ΩP. But this is not true if ΩP = Ω because on the boundary there are

“parasitic” stresses (the η in (2.62)) that ensure that R = Rh, that is, that

the resultant Rh of the BE load (p + parasitic stresses) deviates from the

true resultant R, and therefore the sum of the BE support reactions deviates

from −R; see Tab. 2.1.

These parasitic stresses on Γ are a result of the deviations between the

tractions of the BE solution and the tractions of the exact solution. Usually

the imbalance R − Rh   = 0 is very small (less than 3%).

Theoretically the BE method is restricted to linear problems, because the

scalar product is distributive

_

Ω

G0(p1 + p2) dΩ =

_

Ω

G0 p1 dΩ +

_

Ω

G0 p2 dΩ (2.69)

2.3 Comparison finite elements—boundary elements 267

Table 2.1. FE versus BE

FE BE

differential equation satisfied no yes

geometric boundary conditions yes (no)

static boundary conditions no no

“global equilibrium” yes no

local equilibrium no yes

but in practice the BE method can also be successfully applied to many such

problems by placing the nonlinear terms on the right-hand side and solving

the equations iteratively [257].

The ideal application for boundary elements are dynamic problems in unbounded

domains. An example of such a problem is the analysis of displacements

of the surface of an elastic half-space if a (very) fast train is speeding

across the surface; see Fig. 2.18. The velocity of propagation for the p- and

s-waves respectively is, cp = 995.1 m/s, cs = 300.0 m/s, and the density

of the half-space is ρ = 2000 kg/m3, from which follows that the speed of

propagation for Rayleigh waves is cR = 284.7 m/s. The train is a point load

simulated by a triangular stress distribution. The train moved at different

speeds v across the surface:

0 < v < cR Figure a cR < v < cs Figure b

cs < v < cp Figure c cp < v < +∞ Figure d .

At v < cR (see Fig. 2.18 a) a symmetric depression appeared (as in the static

case on both sides of the point load the deformation was the same) which

moved at speed v across the surface. At speeds cR < v < cs (see Fig. 2.18 b),

a wave front builds up ahead of the moving point load. At speeds cs < v < cp

the wave front becomes wedge-like (see Fig. 2.18 c) while at the same time

the impression the triangular “point load” makes become smaller. At speeds

cp < v = 1350 m/s < ∞, which exceed the p-wave velocity, deformations only

occur after the point load has passed, i.e., they trail the point load (see Fig.

2.18 d).

268 2 What are boundary elements?

Fig. 2.18. Moving point load and surface waves: a) v < cR, b) cR < v < cs,

c) c s < v < cp, d) cp < v < ∞

Навигация