1.45 Generalized finite element methods

Back

In the past decade the FE method has been extended in various directions.

Most of these extensions, as the Element Free Galerkin (EFG) method or the

X-FEM method, can be characterized in the framework of the Generalized

Finite Element Method (GFEM) or as it was called previously the Partition

of Unity Method (PUM).

According to the GFEM common to all these techniques is that the domain

Ω is divided into different regions ωj and that to each region belongs a local

space Vj of functions, not necessarily polynomials, that match the (assumed)

character of the solution and thus ensure good local approximation. Then a

partition of unity is used to “paste” these spaces together to form the trial

space Vh. The partition of unity of the domain Ω may be based on a simple

triangulation and so it offers more freedom when compared to standard FE

methods, [20], where the shape functions and the mesh are closely linked.

214 1 What are finite elements?

Fig. 1.156. X-FEM and horizontal displacement in a cracked bar: a) regular mesh

and discontinuous shape function; b) point load on the left side and c) on the right

side of the crack; d) half-sided load p e) linear load p(x) = p0 ・ x/l

1.45 Generalized finite element methods 215

Meshless methods

A significant stimulus in the development of meshless methods comes from

the problem of moving boundaries as in crack growth. The Element Free

Galerkin (EFG) method and the Particle or the Finite Point method attempt

to overcome this problem.

The EFG is based on a moving least squares scheme that is a method

of reconstructing continuous functions from a set of “arbitrarily” distributed

point samples via the calculation of a weighted least squares measure which is

biased towards the region around the point at which the reconstructed value

is requested.

The advantages of meshfree methods compared to finite elements are

the higher order of the shape functions, simpler incorporation of h- and padaptivity

and certain advantages in crack problems. The most important

drawback of meshfree methods though is their higher computational cost.

The calculation of the strain energy products kij of “free floating” meshless

shape functions is anything but trivial.

X-FEM

The idea of the so-called extended FE method (X-FEM) is to model cracks and

other discontinuities by locally enriching the trial space Vh with discontinuous

shape functions through a partition of unity method.

To fix ideas we consider the bar in Fig. 1.156 a which is subdivided into

four linear elements. We want to model a crack—or rather a fracture—of the

bar at node 2. To this end we introduce at this node an additional shape

function

ψ2(x) = J(x2) · ϕ2(x) J(x2) :=

_

0.5 x > x2

−0.5 x < x2

(1.605)

which is the product of a step function and the unit displacement ϕ2(x) of

the node x2.

The five standard linear shape functions ϕi(x), i = 0, 1, 2, 3, 4, form a partition

of unity

_

i

ϕi(x) = 1 0 ≤ x ≤ l , (1.606)

that is an assemblage of nodal influence functions, and by multiplying the

step function J(x2) with ϕ2(x) the step function is restricted to the immediate

neighborhood of the node x2. That is any additional function which is

introduced to enrich the trial space Vh lives under the umbrella of one of the

nodal shape functions. In the X-FEM terminology this is written as

uh(x) =

_

i

ϕ________i(x) · (ui +

_

j

ψj(x) ui

j)

=

_

i

ϕi(x) · ui +

_

i

ϕi(x) · (ψ1(x) ui

1 + ψ2(x) ui

2)) . (1.607)

216 1 What are finite elements?

Fig. 1.157. Shape functions and derivatives of the shape functions

The extra functions are the functions ψj(x) and the coefficients ui

j are the

additional degrees of freedom. At nodes i where no functions are added the

ui

j are zero. In (1.607) we assumed that two extra functions are added at each

node. Of course the ψj(x) can vary from node to node and some of the ui

j can

be zero.

Our numbering scheme is simply the following:

uh(x) = u1 · ϕ1(x) + u2 · ϕ2(x) + u3 · ψ2(x) + u4 · ϕ3(x) . (1.608)

1.45 Generalized finite element methods 217

The problem we face is that the new function ψ2(x) is discontinuous at x2

which means that the strain energy product

k33 =

_ l

0

EA(ψ

_

2(x))2 dx = ∞ (1.609)

is infinite. Now we could simply ignore the infinite slope and split the integral

into two parts

k33 =

_ x2−ε

0

. . . dx +

_ l

x2+ε

. . . dx =

_

0.5

le

_2

· le +

_

0.5

le

_2

· le = EA · 0.5 le

(1.610)

or we could argue as follows: the discontinuous shape function ψ2(x) can be

considered the limit of the zig-zag function ψε

2(x) in Fig. 1.156 if ε tends to

zero. And if we take the strain energy of the zig-zag function the integral in

the middle

k33 =

_ x2−ε

0

. . . dx +

_ x2+ε

x2−ε

EA[(ψε

2(x))_]2 dx +

_ l

x2+ε

. . . dx (1.611)

explodes if ε → 0

_ x2+ε

x2−ε

EA[(ψε

2(x))_]2 dx _ EA · 0.52

ε2

· 2 ε . (1.612)

To prevent this from happening we let EA = 0 in that part of the bar. The

result is of course the same as before but our treatment of the singularity

seems better justified because that is what we want to model: a bar where

EA = 0 at one point. We will see in the following that the structure follows

exactly this argument.

The stiffness matrix of the bar is calculated as usually, (see Fig. 1.157),

by forming the strain energy products of the shape functions—the fact that

EA = 0 at x2 has no influence on the other values kij

K = EA

le

⎡

⎢⎢⎣

2 −1 0.5 0

−1 2 0 −1

0.5 0 0.5 −0.5

0 −1 −0.5 2

⎤

⎥⎥⎦

. (1.613)

For a first try a point load P = 1 is applied at the node x1 and the FE

program promptly produces the shape in Fig. 1.156 b. Obviously does the

program understand what we want. Simply by adding a discontinuous shape

function the structure develops a crack! And if the point load acts on the

other side of the “crack” the situation is simply reversed, see Fig. 1.156 c.

Also the response to a half-sided constant load p, Fig. 1.156 d, and a load

p(x) = p0 · x/l, Fig. 1.156 e, is correct.

218 1 What are finite elements?

Fig. 1.158. X-FEM and crack between two nodes: a) regular mesh and discontinuous

shape function; b) horizontal displacement due to the point load on the left

side and c) on the right side of the crack

If the crack lies between two nodes (see Fig. 1.158) the discontinuous shape

function extends only from node x1 to x2, (H(ξ) = Heaviside function, which

is zero for ξ < 0 and +1 for ξ > 0),

ψ2(x) = ϕ1(x) · H(x − 1.5) + ϕ2(x) · (H(x − 1.5) − 1) (1.614)

(H(ξ) − 1) = Heaviside on the negative axis) so that

uh(x) = u1 · ϕ1(x) + u2 · ψ2(x) + u3 · ϕ2(x) + u4 · ϕ3(x) . (1.615)

Now the stiffness matrix is

K = EA

le

⎡

⎢⎢⎣

2 1 −1 0

1 1 −1 0

−1 −1 2 −1

0 0 −1 2

⎤

⎥⎥⎦

. (1.616)

And as before do the actions of the point loads stop at the crack, see Fig.

1.158 b and c.

The extension of this technique to 2-D and 3-D problems is straightforward,

see Fig. 1.159,

1.45 Generalized finite element methods 219

Fig. 1.160. Crack-tip functions

Fig. 1.159. Crack between nodes, [169]

220 1 What are finite elements?

Fig. 1.161. Various elements

uh(x) =

_

i

ui ϕi(x) +

_

j∈J

uHj

ϕj(x)H(x) +

_

k∈K

ϕk(x) (

_4

l=1

cl

k Fl(x))

(1.617)

where the Heaviside function H(x) which models the discontinuity at the

crack is multiplied with the shape functions of the circled nodes, set J, and

the asymptotic crack-tip displacement field which consists of four modes (see

Fig. 1.160)

F1(x) =

√

r sin θ

2 F2(x) =

√

r cos θ

2

(1.618)

F3(x) =

√

r sin θ

2

· sinθ F4(x) =

√

r cos θ

2

· sin θ (1.619)

r =

_

x2 + y2 θ = arctan y

x

(1.620)

is multiplied with the shape functions of the set K of squared nodes, [77],

[169].

Навигация