Application of finite element method to three-dimensional elasticity problems

Muhsin J. Jweeg , ... Kadhim Kamil Resan , in Energy Methods and Finite Element Techniques, 2022

9.3 Steps of formulation

The procedure of formulation of the stiffness equations is similar to those presented previously. The complexities of obtaining the stiffness and load vectors will be simplified in detail. The following steps are presented in detail. These steps are summarized as follows,

Step 1 , Define the nodal parameters for the case of 3D elasticity problems, the state of displacements at any point in the domain can be defined in terms of three displacement components u, v, and w in the direction of the x, y, and z axes, respectively. For the 8-node hexahedral element using a local numbering system, the following vectors can be defined,

1.

Nodal displacement vector

At each node, the displacements are u i , v i , and w i , in x i , y i , and z i , directions, where i=1, 8 is the number of nodes of the element and can be arranged in a vector as follows,

(9.7) { δ } t = { u 1 v 1 w 1 u 8 v 8 w 8 }

2.

Nodal loading vector

The loading is in x i , y i , and z i , directions where i=1, 8 is the number of nodes. These can be arranged in a vector containing the element load vector as follows.

(9.8) { F } t = { F x 1 F y 1 F z 1 F x 8 F y 8 F z 8 }

3.

Shape function determination

The shape function N i ( x , y , z ) can be found as follows. Fig. 9.3 shows a 3D element with eight nodes, and by using the natural coordinates ( ξ , η , ζ ) with the local numbering sequence should be adhered to be used in obtaining the system stiffness matrix concerning the global numbering. The displacement field can be assumed as follows,

u ( x , y , z ) = i = 1 8 u i N i ( x , y , z )

v ( x , y , z ) = i = 1 8 v i N i ( x , y , z )

(9.9) w ( x , y , z ) = i = 1 8 w i N i ( x , y , z )

Figure 9.3. 8-Node isoparametric hexahedral element.

To obtain the same function expressions for every 8-node hexahedral element, it is essential to transform the element into a uniform one in an intrinsic space. This may be carried out by transforming the element into a cube of unit side length in the intrinsic ξ , η , ζ space as shown in Fig. 9.3.

Step 2,

Consider one of the displacement components, say (u-component in the x-direction), it can be expressed in terms of ξ , η , ζ using an algebraic polynomial as follows,

u ( ξ , η , ζ ) = α 1 + α 2 ξ + α 3 η + α 4 ζ + α 5 ξ η + α 6 ξ ζ + α 7 η ζ + α 8 ξ η ζ

v ( ξ , η , ζ ) = α 9 + α 10 ξ + α 11 η + α 12 ζ + α 13 ξ η + α 14 ξ ζ + α 15 η ζ + α 16 ξ η ζ

(9.10) w ( ξ , η , ζ ) = α 17 + α 18 ξ + α 19 η + α 20 ζ + α 21 ξ η + α 22 ξ ζ + α 23 η ζ + α 24 ξ η ζ

Using the following eight conditions,

u = u 1 , v = v 1 , w = w 1 , at , ξ = 0 , η = 0 , ζ = 0

u = u 2 , v = v 2 , w = w 2 , at , ξ = 1 , η = 0 , ζ = 0

u = u 3 , v = v 3 , w = w 3 , at , ξ = 1 , η = 1 , ζ = 0

u = u 4 , v = v 4 , w = w 4 , at , ξ = 0 , η = 1 , ζ = 0

u = u 5 , v = v 5 , w = w 5 , at , ξ = 0 , η = 0 , ζ = 1

u = u 6 , v = v 6 , w = w 6 , at , ξ = 1 , η = 0 , ζ = 1

u = u 7 , v = v 7 , w = w 7 , at , ξ = 1 , η = 1 , ζ = 1

(9.11) u = u 8 , v = v 8 , w = w 8 , at , ξ = 0 , η = 1 , ζ = 1

The unknowns α 1 , α 2 , , α 8 in Eq. (9.10) can be found as follows, by substitution of boundary conditions of Eq. (9.5) in Eq. (9.4) as follows, using the conditions of Eq. (9.11), we can write,

u 1 = α 1 , v 1 = α 9 , w 1 = α 17

u 2 = α 1 + α 2

v 2 = α 9 + α 10

w 2 = α 17 + α 18

u 3 = α 1 + α 2 + α 3 + α 5

v 3 = α 9 + α 10 + α 11 + α 13

w 3 = α 17 + α 18 + α 19 + α 21

u 4 = α 1 + α 3

v 4 = α 9 + α 11

w 4 = α 17 + α 19

u 5 = α 1 + α 4

v 5 = α 9 + α 12

w 5 = α 17 + α 20

u 6 = α 1 + α 2 + α 4 + α 6

v 6 = α 9 + α 10 + α 12 + α 14

w 6 = α 17 + α 18 + α 20 + α 22

u 7 = α 1 + α 2 + α 3 + α 4 + α 5 + α 6 + α 7 + α 8

v 7 = α 9 + α 10 + α 11 + α 12 + α 13 + α 14 + α 15 + α 16

w 7 = α 17 + α 18 + α 19 + α 20 + α 21 + α 22 + α 23 + α 24

u 8 = α 1 + α 3 + α 4 + α 7

v 8 = α 9 + α 11 + α 12 + α 15

(9.12) w 8 = α 17 + α 19 + α 20 + α 23

Arrange Eq. (9.12) in matrix form,

(9.13) { u 1 v 1 w 1 u 2 v 2 w 2 u 3 v 3 w 3 u 4 v 4 w 4 u 5 v 5 w 5 u 6 v 6 w 6 u 7 v 7 w 7 u 8 v 8 w 8 } = [ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 1 0 ] { α 1 α 2 α 3 α 4 α 5 α 6 α 7 α 8 α 9 α 10 α 11 α 12 α 13 α 14 α 15 α 16 α 17 α 18 α 19 α 20 α 21 α 22 α 23 α 24 }

Writing Eq. (9.13) in an abbreviated form gives,

{ δ } = [ A ] { α }

Or,

(9.14) { α } = [ A ] 1 { δ }

Where,

[ A ] = [ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 1 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 1 0 ]

Having obtained α s ' from Eq. (9.14), Eq. (9.10) becomes,

u i ( ξ , η , ζ ) / field displacement ( 8 × 1 ) = [ N 1 ( ξ , η , ζ ) N 8 ( ξ , η , ζ ) ] / shape function ( 1 × 8 ) { u 1 u 2 u 8 } / nodal displacement vector ( 8 × 1 )

Similar expressions are obtained for v i ( ξ , η , ζ ) and w i ( ξ , η , ζ ) .

v i ( ξ , η , ζ ) / field displacement ( 8 × 1 ) = [ N 1 ( ξ , η , ζ ) N 8 ( ξ , η , ζ ) ] / shape function ( 1 × 8 ) { v 1 v 2 v 8 } / nodal displacement vector ( 8 × 1 )

(9.15) w i ( ξ , η , ζ ) / field displacement ( 8 × 1 ) = [ N 1 ( ξ , η , ζ ) N 8 ( ξ , η , ζ ) ] / shape function ( 1 × 8 ) { w 1 w 2 w 8 } / nodal displacement vector ( 8 × 1 )

These can be arranged in matrix form,

(9.16) { u i v i w i } / field displacement vector of a particular position ( 24 × 1 ) = [ N 1 0 0 N 2 0 0 N 8 0 0 0 N 1 0 0 N 2 0 0 N 8 0 0 0 N 1 0 0 N 2 0 0 N 8 ] / ( 24 × 24 ) { u 1 v 1 w 1 u 8 v 8 w 8 } / ( 24 × 1 )

Therefore we can write the expressions for the shape functions as follows,

The boundary conditions in Eq. (9.11) are substituted in Eq. (9.10), eight equations are obtained and can be arranged as follows,

The unknown parameters α 1 , α 2 , , α 8 can be obtained using the above boundary conditions. Hence, it can be shown that,

N 1 ( ξ , η , ζ ) = ( 1 ξ ) ( 1 η ) ( 1 ζ )

N 2 ( ξ , η , ζ ) = ( ξ ) ( 1 η ) ( 1 ζ )

N 3 ( ξ , η , ζ ) = ( ξ ) ( η ) ( 1 ζ )

N 4 ( ξ , η , ζ ) = ( 1 ξ ) ( η ) ( 1 ζ )

N 5 ( ξ , η , ζ ) = ( 1 ξ ) ( 1 η ) ( ζ )

N 6 ( ξ , η , ζ ) = ( ξ ) ( 1 η ) ( ζ )

N 7 ( ξ , η , ζ ) = ( ξ ) ( η ) ( ζ )

(9.17) N 8 ( ξ , η , ζ ) = ( 1 ξ ) ( η ) ( ζ )

The relationship between the Cartesian system (x, y, z) and the intrinsic system ( ξ , η , ζ ) can be obtained by employing the isoparametric transformation,

x ( ξ , η , ζ ) = i = 1 8 x i N i ( ξ , η , ζ )

y ( ξ , η , ζ ) = i = 1 8 y i N i ( ξ , η , ζ )

(9.18) z ( ξ , η , ζ ) = i = 1 8 z i N i ( ξ , η , ζ )

Step 3, Express the strain components at any point in terms of nodal displacements and shape functions, we have, from the theory of elasticity, the strain components in 3D problems are as follows,

ε x = u x = i = 1 8 u i N i x

ε y = v y = i = 1 8 v i N i y

ε z = w z = i = 1 8 w i N i z

γ x y = u y + v x = i = 1 8 ( u i N i y + v i N i x )

γ y z = v z + w y = i = 1 8 ( v i N i z + w i N i y )

(9.19) γ z x = w z + u x = i = 1 8 ( u i N i x + w i N i z )

Defining the strain tensor,

(9.20) { ε } t = { ε x ε y ε z γ x y γ y z γ z x }

It can be shown that,

(9.21) { ε } ( 6 × 1 ) = [ B ] ( 6 × 24 ) { δ } ( 24 × 1 )

Where,

[ B ] t = [ B 1 B 2 B 6 ]

(9.22) [ B ] = [ N i x 0 0 0 N i y 0 0 0 N i z N i y N i x 0 0 N i z N i y N i z 0 N i x ] :

Where [B] is the strain–displacement matrix of size (6×24). It is a function of natural coordinates f ( ξ , η , ζ )

The following procedure can be employed in order to obtain the Cartesian derivative of the intrinsic shape functions.

Consider a function f ( ξ , η , ζ ) . Using the chain rule of partial differentiation, it can be deduced that,

f ξ = f x x ξ + f y y ξ + f z z ξ

f η = f x x η + f y y η + f z z η

(9.23) f ζ = f x x ζ + f y y ζ + f z z ζ

The above equations can be expressed in matrix form as follows,

(9.24) { f ξ f η f ζ } = [ J ] [ f x f y f z ]

Where, [ J ] is the 3D Jacobian matrix defined as follows,

(9.25) [ J ] = [ x ξ y ξ z ξ x η y η z η x ζ y ζ z ζ ]

Applying the above analysis to the intrinsic shape functions of the hexahedral element, it can be deduced that,

(9.26) { N i x N i y N i z } = [ J ] 1 { N i ξ N i η N i ζ } , where , i = ( 1 , 2 , , 8 )

And,

| J | = x 14 ( y 24 z 34 y 34 z 24 ) + y 14 ( z 24 x 34 z 34 x 24 ) + z 14 ( x 24 y 34 x 34 y 24 )

Therefore the strain–displacement matrix [B] is obtained by using Eqs. (9.25) and (9.26) as follows,

{ N i x N i y N i z } = 1 | J | [ ( y η z ζ z η y ζ ) ( z η x ζ x η z ζ ) ( x η y ζ y η x ζ ) ( z ξ y ζ y ξ z ζ ) ( x ξ z ζ z ξ x ζ ) ( y ξ x ζ x ξ y ζ ) ( y ξ z η z ξ y η ) ( z ξ x η x ξ z η ) ( x ξ y η y ξ x η ) ] { N i ξ N i η N i ζ }

We have,

{ ε x ε y ε z γ x y γ y z γ z x } = [ B ( 1 , 1 ) 0 0 B ( 1 , 22 ) 0 0 0 B ( 2 , 2 ) 0 0 B ( 1 , 23 ) 0 0 0 B ( 3 , 3 ) 0 0 B ( 1 , 24 ) B ( 4 , 1 ) B ( 4 , 2 ) 0 B ( 4 , 22 ) B ( 4 , 23 ) 0 0 B ( 5 , 2 ) B ( 5 , 3 ) 0 B ( 5 , 23 ) B ( 5 , 24 ) B ( 6 , 1 ) 0 B ( 6 , 3 ) B ( 6 , 22 ) 0 B ( 6 , 24 ) ] { u 1 v 1 w 1 u 8 v 8 w 8 }

Or,

{ ε } = [ B ] { δ }

Where,

B ( 1 , 1 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( 1 ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( 1 η ) )

B ( 1 , 4 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( 1 ζ ) ( z η x ζ x η z ζ ) ( ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( ξ ) ( 1 η ) )

B ( 1 , 7 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( ξ ) ( η ) )

B ( 1 , 10 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( 1 ξ ) ( η ) )

B ( 1 , 13 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( ζ ) ( z η x ζ x η z ζ ) ( 1 ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( 1 η ) )

B ( 1 , 16 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( ζ ) ( z η x ζ x η z ζ ) ( ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( ξ ) ( 1 η ) )

B ( 1 , 19 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( ζ ) + ( z η x ζ x η z ζ ) ( ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( ξ ) ( η ) )

B ( 1 , 22 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( η ) )

B ( 1 , 2 ) = B ( 1 , 3 ) = B ( 1 , 5 ) = B ( 1 , 6 ) = B ( 1 , 8 ) = B ( 1 , 9 ) = B ( 1 , 11 ) = B ( 1 , 12 ) = B ( 1 , 14 ) = B ( 1 , 15 ) = B ( 1 , 17 ) = B ( 1 , 18 ) = B ( 1 , 20 ) = B ( 1 , 21 ) = B ( 1 , 23 ) = B ( 1 , 24 ) = 0

B ( 2 , 2 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( 1 ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( 1 η ) )

B ( 2 , 5 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( 1 ζ ) ( x ξ z ζ z ξ x ζ ) ( ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( ξ ) ( 1 η ) )

B ( 2 , 8 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( ξ ) ( η ) )

B ( 2 , 11 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( η ) )

B ( 2 , 14 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( ζ ) ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( 1 η ) )

B ( 2 , 17 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( ζ ) ( x ξ z ζ z ξ x ζ ) ( ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( ξ ) ( 1 η ) )

B ( 2 , 20 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( ζ ) + ( x ξ z ζ z ξ x ζ ) ( ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( ξ ) ( η ) )

B ( 2 , 23 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( η ) )

B ( 2 , 1 ) = B ( 2 , 3 ) = B ( 2 , 4 ) = B ( 2 , 6 ) = B ( 2 , 7 ) = B ( 2 , 9 ) = B ( 2 , 10 ) = B ( 2 , 12 ) = B ( 2 , 13 ) = B ( 2 , 15 ) = B ( 2 , 16 ) = B ( 2 , 18 ) = B ( 2 , 19 ) = B ( 2 , 21 ) = B ( 2 , 22 ) = B ( 2 , 24 ) = 0

B ( 3 , 3 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( 1 ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( 1 η ) )

B ( 3 , 6 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( 1 ζ ) ( z ξ x η x ξ z η ) ( ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( ξ ) ( 1 η ) )

B ( 3 , 9 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( ξ ) ( η ) )

B ( 3 , 12 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( 1 ξ ) ( η ) )

B ( 3 , 15 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( ζ ) ( z ξ x η x ξ z η ) ( 1 ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( 1 η ) )

B ( 3 , 18 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( ζ ) ( z ξ x η x ξ z η ) ( ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( ξ ) ( 1 η ) )

B ( 3 , 21 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( ζ ) + ( z ξ x η x ξ z η ) ( ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( ξ ) ( η ) )

B ( 3 , 24 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( η ) )

B ( 3 , 1 ) = B ( 3 , 2 ) = B ( 3 , 4 ) = B ( 3 , 5 ) = B ( 3 , 7 ) = B ( 3 , 8 ) = B ( 3 , 10 ) = B ( 3 , 11 ) = B ( 3 , 13 ) = B ( 3 , 14 ) = B ( 3 , 16 ) = B ( 3 , 17 ) = B ( 3 , 19 ) = B ( 3 , 20 ) = B ( 3 , 22 ) = B ( 3 , 23 ) = 0

B ( 4 , 1 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( 1 ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( 1 η ) )

B ( 4 , 2 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( 1 ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( 1 η ) )

B ( 4 , 4 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( 1 ζ ) ( x ξ z ζ z ξ x ζ ) ( ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( ξ ) ( 1 η ) )

B ( 4 , 5 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( 1 ζ ) ( z η x ζ x η z ζ ) ( ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( ξ ) ( 1 η ) )

B ( 4 , 7 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( ξ ) ( η ) )

B ( 4 , 8 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( ξ ) ( η ) )

B ( 4 , 10 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( η ) )

B ( 4 , 11 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( 1 ξ ) ( η ) )

B ( 4 , 13 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( ζ ) ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( 1 η ) )

B ( 4 , 14 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( ζ ) ( z η x ζ x η z ζ ) ( 1 ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( 1 η ) )

B ( 4 , 16 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( ζ ) ( x ξ z ζ z ξ x ζ ) ( ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( ξ ) ( 1 η ) )

B ( 4 , 17 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( ζ ) ( z η x ζ x η z ζ ) ( ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( ξ ) ( 1 η ) )

B ( 4 , 19 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( ζ ) + ( x ξ z ζ z ξ x ζ ) ( ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( ξ ) ( η ) )

B ( 4 , 20 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( ζ ) + ( z η x ζ x η z ζ ) ( ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( ξ ) ( η ) )

B ( 4 , 22 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( η ) )

B ( 4 , 23 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( η ) )

B ( 4 , 3 ) = B ( 4 , 6 ) = B ( 4 , 9 ) = B ( 4 , 12 ) = B ( 4 , 15 ) = B ( 4 , 18 ) = B ( 4 , 21 ) = B ( 4 , 24 ) = 0

B ( 5 , 2 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( 1 ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( 1 η ) )

B ( 5 , 3 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( 1 ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( 1 η ) )

B ( 5 , 5 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( 1 ζ ) ( z ξ x η x ξ z η ) ( ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( ξ ) ( 1 η ) )

B ( 5 , 6 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( 1 ζ ) ( x ξ z ζ z ξ x ζ ) ( ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( ξ ) ( 1 η ) )

B ( 5 , 8 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( ξ ) ( η ) )

B ( 5 , 9 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( ξ ) ( η ) )

B ( 5 , 11 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( 1 ξ ) ( η ) )

B ( 5 , 12 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( 1 ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( 1 ζ ) ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( η ) )

B ( 5 , 14 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( ζ ) ( z ξ x η x ξ z η ) ( 1 ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( 1 η ) )

B ( 5 , 15 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( ζ ) ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( 1 η ) )

B ( 5 , 17 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( ζ ) ( z ξ x η x ξ z η ) ( ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( ξ ) ( 1 η ) )

B ( 5 , 18 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( 1 η ) ( ζ ) ( x ξ z ζ z ξ x ζ ) ( ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( ξ ) ( 1 η ) )

B ( 5 , 20 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( ζ ) + ( z ξ x η x ξ z η ) ( ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( ξ ) ( η ) )

B ( 5 , 21 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( ζ ) + ( x ξ z ζ z ξ x ζ ) ( ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( ξ ) ( η ) )

B ( 5 , 23 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( η ) )

B ( 5 , 24 ) = 1 | J | ( ( z ξ y ζ y ξ z ζ ) ( η ) ( ζ ) + ( x ξ z ζ z ξ x ζ ) ( 1 ξ ) ( ζ ) + ( y ξ x ζ x ξ y ζ ) ( 1 ξ ) ( η ) )

B ( 5 , 1 ) = B ( 5 , 4 ) = B ( 5 , 7 ) = B ( 5 , 10 ) = B ( 5 , 13 ) = B ( 5 , 16 ) = B ( 5 , 19 ) = B ( 5 , 22 ) = 0

B ( 6 , 1 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( 1 ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( 1 η ) )

B ( 6 , 3 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( 1 ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( 1 η ) )

B ( 6 , 4 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( 1 ζ ) ( z ξ x η x ξ z η ) ( ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( ξ ) ( 1 η ) )

B ( 6 , 6 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( 1 ζ ) ( z η x ζ x η z ζ ) ( ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( ξ ) ( 1 η ) )

B ( 6 , 7 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( ξ ) ( η ) )

B ( 6 , 9 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( ξ ) ( η ) )

B ( 6 , 10 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( 1 ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( 1 ζ ) ( x ξ y η y ξ x η ) ( 1 ξ ) ( η ) )

B ( 6 , 12 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( 1 ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( 1 ζ ) ( x η y ζ y η x ζ ) ( 1 ξ ) ( η ) )

B ( 6 , 13 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( ζ ) ( z ξ x η x ξ z η ) ( 1 ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( 1 η ) )

B ( 6 , 15 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( ζ ) ( z η x ζ x η z ζ ) ( 1 ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( 1 η ) )

B ( 6 , 16 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( 1 η ) ( ζ ) ( z ξ x η x ξ z η ) ( ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( ξ ) ( 1 η ) )

B ( 6 , 18 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( 1 η ) ( ζ ) ( z η x ζ x η z ζ ) ( ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( ξ ) ( 1 η ) )

B ( 6 , 19 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( ζ ) + ( z ξ x η x ξ z η ) ( ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( ξ ) ( η ) )

B ( 6 , 21 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( ζ ) + ( z η x ζ x η z ζ ) ( ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( ξ ) ( η ) )

B ( 6 , 22 ) = 1 | J | ( ( y ξ z η z ξ y η ) ( η ) ( ζ ) + ( z ξ x η x ξ z η ) ( 1 ξ ) ( ζ ) + ( x ξ y η y ξ x η ) ( 1 ξ ) ( η ) )

B ( 6 , 24 ) = 1 | J | ( ( y η z ζ z η y ζ ) ( η ) ( ζ ) + ( z η x ζ x η z ζ ) ( 1 ξ ) ( ζ ) + ( x η y ζ y η x ζ ) ( 1 ξ ) ( η ) )

(9.27) B ( 6 , 2 ) = B ( 6 , 5 ) = B ( 6 , 8 ) = B ( 6 , 11 ) = B ( 6 , 14 ) = B ( 6 , 17 ) = B ( 6 , 20 ) = B ( 6 , 23 ) = 0

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780323886666000083

NUMERICAL MODELLING OF THE NON-LINEAR BEHAVIOUR OF MATERIALS WITH APPLICATIONS TO METALS

S.K. Choi , ... P.F. Thomson , in Advances in Engineering Plasticity and its Applications, 1993

3.2 Finite Element Discretisation

The finite element method uses piecewise continuous interpolation (or shape) functions to describe the distribution of the spatially discretised unknown nodal parameters [e.g. 6] ϕI(t):

(6) ϕ ( x , t ) = ϕ i ( t ) N I ( x )

where NI is the interpolation function associated with node I. The repeated index I indicates summation with respect to all the nodes associated with an element. The velocity gradient is simply given by

(7) v i , j = N I x j v iI

Selective integration [e.g. 7] was again employed to compute the deviatoric and hydrostatic components of the strain tensor.

In essence, the finite element method, through semidiscretisation in space, results in a weak form of the equations of motion with time as the only independent variable which, by applying numerical integration, gives

(8) M v ˙ + Cv = f

where M is the mass matrix, C is the damping matrix, v is the nodal velocity vector, f is the nodal total force vector, and the superposed dot denotes differentiation with respect to time.

By diagonalising the mass matrix, the same solution strategy described in this paper can be used for both the finite element and the contour integral finite difference methods.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780444899910500949

Overview of the ATILA finite element method (FEM) software code

K. Uchino , in Applications of ATILA FEM Software to Smart Materials, 2013

1.3.2 Shape functions

The finite element method, being an approximation process, will result in the determination of an approximate solution of the form

[1.17] w w ^ = N i a i

where the N i are shape functions prescribed in terms of independent variables (such as coordinates) and the a i are nodal parameters, known or unknown. The shape functions must guarantee the continuity of the geometry between elements. Moreover, to ensure convergence, it is necessary that the shape functions be at least C m continuous if derivatives of the mth degree exist in the integral form. This condition is automatically met if the shape functions are polynomials complete to the mth order.

The construction of shape functions for an element Ωe, defined by n nodes, usually requires that if N i is the shape function for node i, then N i  =   1 at node i, and N i  =   0 at the other nodes. Also, for any point p in Ωe we must have

[1.18] 1 n N i p = 1

Polynomials are commonly used to construct shape functions. For instance a Lagrange polynomial

[1.19] N i ξ = k = 1 k 1 n ξ k ξ ξ k ξ i

can be used at node i of a one-dimensional element containing n nodes (see Fig. 1.4). It verifies the following conditions.

1.4. A generalized n-node linear element.

[1.20] { N i ξ i = 1 N i ξ k 1 = 0 1 n N i ξ = 1

For the four-node rectangular element of Fig. 1.5, we can write the four shape functions as:

1.5. An example of a four-node quadrilateral element

[1.21] N ξ η = N i ξ N j η

where i and j indicate the row and column of the node in the element and i, j   =   1, 2. The conditions described by Equation [1.20] can be verified for these functions. Finally, the position of any point p with coordinates (x, h) in the element is given by:

[1.22] ξ η = 1 4 1 + ξ 1 + η 1 4 1 + ξ 1 + η 1 4 1 + ξ 1 + η 1 4 1 + ξ 1 + η 1 1 1 1 1 1 1 1

In most cases, the same shape functions are used to describe the element geometry and to represent the solution w ^ .

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780857090652500011

Frames

Dimitrios G. Pavlou PhD , in Essentials of the Finite Element Method, 2015

7.1 Framed Structures

Frames are structures consisting of beams (Figure 7.1 ). The beams composing a frame generally have arbitrary directions. Therefore, all nodal parameters (displacements, rotations, forces, moments) should be expressed with respect to a common system of coordinates, the global coordinate system.

Figure 7.1. Framed structure.

Since vertical deflections of node 2 of a horizontal beam 2-3 can cause axial deflections to the vertical beam 1-2 connected to the same node (Figure 7.2), the axial component of deflections should be taken into account to the beam element equation of a frame member.

Figure 7.2. The displacement e of the node 2 causes axial deflection on the element 1 and vertical deflection on the element 2.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128023860000074

Galerkin Method of Approximation: Irreducible and Mixed Forms

O.C. Zienkiewicz , ... David Fox , in The Finite Element Method for Solid and Structural Mechanics (Seventh Edition), 2014

2.5.3 Mixed displacement/traction condition

The treatment of a mixed condition in which some displacement components are specified together with some traction components often requires a change in the nodal parameters. For example, a shaft with axis in the x 3 direction and radius R that rotates inside a bearing (without friction or gaps) requires u n = u r ( R ) = 0 and t θ ( R ) = t z ( R ) = 0 (where the coordinate origin is placed at the center of the shaft). In this case it is necessary to transform the degrees of freedom at each node on the boundary of the shaft such that

(2.41) u a = ( u ̃ 1 ) a ( u ̃ 2 ) a ( u ̃ 3 ) a = cos θ a - sin θ a 0 sin θ a cos θ a 0 0 0 1 ( u ̃ r ) a ( u ̃ θ ) a ( u ̃ z ) a = L a u ̃ a

This transformation is then applied to a residual as

(2.42a) R a = L a T R a

and to the mass and stiffness as

(2.42b) M a b = L a T M ab L b , M a c = L a T M ac , M cb = M cb L b K a b = L a T K ab L b , K a c = L a T K ac , K cb = K cb L b

where a and b belong to transformed nodes and c to a node which retains its original orientation. It is usually convenient to perform these transformations on each individual element; however, if desired they can be applied to the assembled arrays.

Once the transformation is performed, each individual displacement and traction condition may be imposed as described above.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9781856176347000028

Elasticity: Two- and Three-Dimensional Finite Elements

O.C. Zienkiewicz , ... J.Z. Zhu , in The Finite Element Method: its Basis and Fundamentals (Seventh Edition), 2013

7.10.2 Direct minimization

The fact that the finite element approximation reduces to the problem of minimizing the total potential energy Π defined in terms of a finite number of nodal parameters led us to the formulation of the simultaneous set of equations given symbolically by Eq. (7.54). This is the most usual and convenient approach, especially in linear solutions, but other search procedures, now well developed in the field of optimization, could be used to estimate the lowest value of Π . In this text we shall continue with the simultaneous equation process but the interested reader could well bear the alternative possibilities in mind [42, 43].

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9781856176330000071

Introduction

J.E. Akin , in Finite Element Analysis with Error Estimators, 2005

1.3 Outline of finite element procedures

From the mathematical point of view the finite element method is an integral formulation. Modern finite element integral formulations are usually obtained by either of two different procedures: weighted residual or variational formulations. The following sections briefly outline the common procedures for establishing finite element models. It is fortunate that all these techniques use the same bookkeeping operations to generate the final assembly of algebraic equations that must be solved for the unknowns.

The generation of finite element models by the utilization of weighted residual techniques is increasingly important in the solution of differential equations for nonstructural applications. The weighted residual method starts with the governing differential equation to be satisfied in a domain Ω:

(1.1) L ( ϕ ) = Q ,

where L denotes a differential operator acting on the primary unknown, ϕ, and Q is a source term. Generally we assume an approximate solution, say ϕ*, for the spatial distribution of the unknown, say

(1.2) ϕ ( x ) ϕ * = i n h i ( x ) Φ i * ,

where the hi (x) are spatial distributions associated with the coefficient Φ* i . That assumption leads to a corresponding assumption for the spatial gradient of the assumed behavior. Next we substitute these assumptions on spatial distributions into the differential equation. Since the assumption is approximate, this operation defines a residual error term, R, in the differential equation

(1.3) L ( ϕ * ) Q = R 0.

Although we cannot force the residual term to vanish, it is possible to force a weighted integral, over the solution domain, of the residual to vanish. That is, the integral of the product of the residual term and some weighting function is set equal to zero, so that

(1.4) I i Ω R w i d Ω = 0

leads to the same number of equations as there are unknown Φ* i values. Most of the time we will find it very useful to employ integration by parts on this governing integral.

Substituting an assumed spatial behavior for the approximate solution, ϕ*, and the weighting function, w, results in a set of algebraic equations that can be solved for the unknown nodal coefficients in the approximate solution. This is because the unknown coefficients can be pulled out of the spatial integrals involved in the assembly process.

The choice of weighting function defines the type of weighted residual technique being utilized. The Galerkin criterion selects

(1.5) w i = h i ( x ) ,

to make the residual error orthogonal to the approximate solution. Use of integration by parts with the Galerkin procedure (i.e., the Divergence Theorem) reduces the continuity requirements of the approximating functions. If a Euler variational procedure exists, the Galerkin criterion will lead to the same element matrices.

A spatial interpolation, or blending function is assumed for the purpose of relating the quantity of interest within the element in terms of the values of the nodal parameters at the nodes connected to that particular element. For both weighted residual and variational formulations, the following restrictions are accepted for establishing convergence of the finite element model as the mesh refinement increases:

1.

The element interpolation functions must be capable of modeling any constant values of the dependent variable or its derivatives, to the order present in the defining integral statement, in the limit as the element size decreases.

2.

The element interpolation functions should be chosen so that at element interfaces the dependent variable and its derivatives, of one order less than those occurring in the defining integral statement, are continuous.

Through the assumption of the spatial interpolations, the variables of interest and their derivatives are uniquely specified throughout the solution domain by the nodal parameters associated with the nodal points of the system. The parameters at a particular node directly influence only the elements connected to that particular node. The domain will be split into a mesh. That will require that we establish some bookkeeping processes to keep up with data going to, or coming from a node or element. Those processes are commonly called gather and scatter, respectively. Figure 1.3 shows some of these processes for a simple mesh with one generalized scalar unknown per node, ng = 1, in a one-dimensional physical space. There the system node numbers are shown numbered in an arbitrary fashion. To establish the local element space domain we must usually gather the coordinates of each of its nodes. For example, for element b(= 5) it gathers the data for system node numbers 6 and 4, respectively, so that the element length, L (5) = x 6x 4, can be computed. Usually we also have to gather some data on the coefficients in the differential equation (material properties). If the coefficients vary over space they may be supplied as data at the nodes that must also be gathered to form the element matrices.

Figure 1.3. Gather and scatter concepts for finite elements

After the element behavior has been described by spatial assumptions, then the derivatives of the space functions are used to approximate the spatial derivatives required in the integral form. The remaining fundamental problem is to establish the element matrices, S e and C e . This involves substituting the approximation space functions and their derivatives into the governing integral form and moving the unknown coefficients, D e , outside the integrals. Historically, the resulting matrices have been called the element stiffness matrix and load vector, respectively.

Once the element equations have been established the contribution of each element is added, using its topology (or connectivity), to form the system equations. The system of algebraic equations resulting from FEA (of a linear system) will be of the form SD = C. The vector D contains the unknown nodal parameters, and the matrices S and C are obtained by assembling the known element matrices, S e and C e , respectively. Figure 1.3 shows how the local coefficients of the element source vector, C e , are scattered and added into the resultant system source, C. That illustration shows a conversion of local row numbers to the corresponding system row numbers (by using the element connectivity data). An identical conversion is used to convert the local and system column numbers needed in assembling each S e into S. In the majority of problems S e , and thus, S, will be symmetric. Also, the system square matrix, S, is usually banded about the diagonal or at least sparse. If S is unsymmetric its upper and lower triangles have the same sparsity.

After the system equations have been assembled, it is necessary to apply the essential boundary constraints before solving for the unknown nodal parameters. The most common types of essential boundary conditions (EBC) are (1) defining explicit values of the parameter at a node and (2) defining constraint equations that are linear combinations of the unknown nodal quantities. The latter constraints are often referred to in the literature as multi-point constraints (MPC). An essential boundary condition should not be confused with a forcing condition of the type that involves a flux or traction on the boundary of one or more elements. These element boundary source, or forcing, terms contribute additional terms to the governing integral form and thus to the element square and/or column matrices for the elements on which the sources were applied. Thus, although these (Neumann-type, and Robin or mixed-type) conditions do enter into the system equations, their presence may not be obvious at the system level. Wherever essential boundary conditions do not act on part of the boundary, then at such locations, source terms from a lower order differential equation automatically apply. If one does not supply data for the source terms, then they default to zero. Such portions of the boundary are said to be subject to natural boundary conditions (NBC). The natural boundary condition varies with the integral form, and typical examples will appear later.

The initial sparseness (the relative percentage of zero entries) of the square matrix, S, is an important consideration since we only want to store non-zero terms. If we employ a direct solver then many initially zero terms will become non-zero during the solution process and the assigned storage must allow for that. The 'fill-in' depends on the numbering of the nodes. If the FEA system being employed does not have an automatic renumbering system to increase sparseness, then the user must learn how to number nodes (or elements) efficiently. After the system algebraic equations have been solved for the unknown nodal parameters, it is usually necessary to output the parameters, D. For every essential boundary condition on D, there is a corresponding unknown reaction term in C that can be computed after D is known. These usually have physical meanings and should be output to help check the results.

In rare cases the problem would be considered completed at this point, but in most cases it is necessary to use the calculated values of the nodal parameters to calculate other quantities of interest. For example, in stress analysis we use the calculated nodal displacements to solve for the strains and stresses. All adaptive programs must do a very large amount of post-processing to be sure that the solution, D, has been obtained to the level of accuracy specified by the analyst. Figure 1.3 also shows that the gather operation is needed again for extracting the local results, D e , from the total results, D, so they can be employed in special element post-processing and/or error estimates.

Usually the post-processing calculations involve determining the spatial derivatives of the solution throughout the mesh. Those gradients are continuous within each element domain, but are discontinuous across the inter-element boundaries. The true solution usually has continuous derivatives so it is desirable to somehow average the individual element gradient estimates to create continuous gradient estimate values that can be reported at the nodes. Fortunately, this addition gradient averaging process also provides new information that allows the estimate of the problem error norm to be calculated. That gradient averaging process will be presented in Chapters 2 and 6.

In the next chapter we will review the historical approach of the method of weighted residuals and its extension to finite element analysis. The earliest formulations for finite element models were based on variational techniques. This is especially true in the areas of structural mechanics and stress analysis. Modern analysis in these areas has come to rely on FEA almost exclusively. Variational models find the nodal parameters that yield a minimum (or stationary) value of an integral known as a functional. In most cases it is possible to assign a physical meaning to the integral. For example, in solid mechanics the integral represents the total potential energy, whereas in a fluid mechanics problem it may correspond to the rate of entropy production. Most physical problems with variational formulations result in quadratic forms that yield algebraic equations for the system which are symmetric and positive definite. The solution that yields a minimum value of the integral functional and satisfies the essential boundary conditions is equivalent to the solution of an associated differential equation. This is known as the Euler theorem.

Compared to the method of weighted residuals, where we start with the differential equation, it may seem strange to start a variational formulation with an integral form and then check to see if it corresponds to the differential equation we want. However, from Euler's work more than two centuries ago we know the variational forms of most even order differential equations that appear in science, engineering, and applied mathematics. This is especially true for elliptical equations. Euler's Theorem of Variational Calculus states that the solution, u, that satisfies the essential boundary conditions and renders stationary the functional

(1.6) I = Ω f ( x , y , z , ϕ , ϕ x , ϕ y , ϕ z ) d Ω + Γ ( q ϕ + a ϕ 2 / 2 ) d Γ

also satisfies the partial differential equation

(1.7) f ϕ x f ( ϕ / x ) y f ( ϕ / y ) z f ( ϕ / z ) = 0

in Ω, and satisfies the natural boundary condition that

(1.8) n x f ( ϕ / x ) + n y f ( ϕ / y ) + n z f ( ϕ / z ) + q + a ϕ = 0

on Γ that is not subject to an essential boundary. Here nx, ny, nz are the components of the normal vector on the boundary, Γ. Note that this theorem also defines the natural boundary condition, as well as the corresponding differential equation. In Chapter 7 we will examine some common Euler variational forms for finite element analysis.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780750667227500321

Numerical and Computational Methods

H. Askes , ... R. de Borst , in Comprehensive Structural Integrity, 2003

3.07.2.1.5 Imposition of essential boundary conditions

One aspect of the EFG method that has been paid much attention in the literature is the enforcement of constraint equations. Especially the imposition of essential (Dirichlet) boundary conditions has been focus of interest. Since EFG shape functions are in general not interpolating, the nodal parameters (contained in the vector u in Equation (54)) are not the nodal values of the approximant function u h (x). Therefore, the imposition of constraint equations is not trivial. Below, a number of strategies are summarized that have been proposed to deal with essential boundary conditions in the EFG method. It is noted that the same problems are encountered for point loads (Zuohui, 2000) and tyings between different degrees of freedom (Pannachet and Askes, 2000).

1.

Collocation methods. In collocation methods, the essential boundary conditions are prescribed at the nodes. As has been pointed out above, a direct collocation may lead to problems. Therefore, several modified collocation techniques have been developed.

Direct collocation. From Figure 6 it can be seen that the EFG shape functions of boundary nodes are almost equal to 1 in their nodes. From this it could be concluded that direct collocation can give reasonable results. This somewhat opportunistic approach has been studied in the original articles on the EFG method, and it has been found that severe errors can be encountered (Belytschko et al., 1994; Lu et al., 1994). Therefore, this approach should be avoided.

Modified collocation. Several attempts have been made to reformulate the collocation principle (Zhu and Atluri, 1998; Atluri et al., 1999; Wagner and Liu, 2000; Chen and Wang, 2000; Wu and Plesha, 2002). The nodal parameter u i is expressed in terms of the value of the approximant function at the node u h (x i ). It can be written that

(60) u h ( x i ) = j N j ( x i ) u j

so that

(61) u j = j [ N j ( x i ) ] 1 u h ( x i )

When Equation (61) is substituted in the global system of equations, essential boundary conditions can be imposed directly. Disadvantages are that the inverse of (part of) the shape function matrix must be determined, and that the bandwidth of the global stiffness matrix may increase.

2.

Weak imposition. Rather than prescribing the essential boundary conditions at the nodes, they can be prescribed in a weak sense. Extra integrals are introduced in the Galerkin formulation. As a result, the boundary conditions are not prescribed pointwise, but in an averaged sense.

Lagrange multipliers. In the first contribution on the EFG method, the Lagrange multiplier technique was proposed to circumvent the problems encountered with direct collocation (Belytschko et al., 1994). Although this approach leads to an increased size of the stiffness matrix, loss of positive-definiteness and bandedness, and the application to explicit dynamic problems is not direct, it has been used in various static applications (Belytschko et al., 1995; Krysl and Belytschko, 1996a, 1996b; Hegen, 1996; Askes et al., 1999; Huerta and Fernández-Méndez, 2001). A special form of the Lagrange multiplier technique has been reported in which the system of equations is re-partitioned, and which bears strong similarities with the modified collocation approach (Wu and Plesha, 2002).

Modified variational principle. To overcome the disadvantages of the Lagrange multiplier approach, it has been proposed to replace the Lagrange multipliers by their physical counterparts, i.e., the tractions on the boundary (Lu et al., 1994; Belytschko et al., 1995; Mukherjee and Mukherjee, 1997; Gavete et al., 2000). This leads to a modified variational approach, similar to mixed-field interpolations. Bandedness and positive-definiteness are maintained, while the size of the stiffness matrix is unaffected. Thus, this approach fits well in a standard numerical solution, although accuracy is lower compared to the Lagrange multiplier approach (Lu et al., 1994; Belytschko et al., 1995).

Better results with the same methodology have been reported when the nodal parameters u i in Equation (48) are replaced by u h (x i ) (Y. X. Mukherjee and S. Mukherjee, 1997). This does not affect the EFG shape functions, but the discretized system of the equations takes a somewhat different format.

Penalty functions. Prescribing the essential boundary conditions by means of penalty formulations leads to a banded, positive-definite stiffness matrix. A disadvantage is that a (large) value for the penalty parameter must be selected by the user, and the stiffness matrix can become poor conditioned. However, high accuracies are reported for a large range of penalty parameters (Zhu and Atluri, 1998; Gavete et al., 2000).

Augmented Lagrangian. The combination of penalty terms and Lagrange multiplier terms into the energy functional results in a so-called augmented Lagrangian approach (Ventura 2002). It eliminates the shortcomings of both Lagrange multipliers and penalty functions, whereas the advantages of either method are maintained. Good results have been reported for the augmented Lagrangian approach (Ventura, 2002). Irrespective of the method that is chosen, a weak imposition of the constraint equation may lead to inaccuracies in case the boundary conditions change abruptly (Pannachet and Askes, 2000). The reason is that EFG shape functions have a relatively large domain of influence (compared to FE shape functions) and they are smooth. Therefore, they are not well suited to describe C –1 or C 0 fields (Jirásek, 1998). It has been found that inaccuracies encountered with Lagrange multipliers or with penalty formulations can be largely overcome by applying nodal integration to the constraint integrals (Pannachet and Askes, 2000).

3.

Changing the shape functions. The problems arising from noninterpolating shape functions can be avoided by modifying the shape functions a priori.

Singular weight functions. Instead of solving the problem to prescribe essential boundary conditions a posteriori (i.e., after the shape functions have been computed), a priori modifications of the original MLS formulation have been proposed. The origin of the noninterpolating character of EFG shape functions lies in the (assumed) finite amplitude of the MLS weight function. It has been shown that an interpolating MLS scheme is obtained when the weight functions are chosen to be infinitely large at their associated nodes (Lancaster and Salkauskas, 1981). Thus, using singular weight functions, essential boundary conditions can be enforced straightforwardly (Kaljević and Saigal, 1997; Chen and Wang, 2000).

Although this strategy is elegant and sound from a theoretical point of view, some practical difficulties remain. The resulting EFG shape functions are not as smooth as with regular weight functions (Hegen, 1996, 1997), therefore, this feature of EFG shape functions cannot be exploited anymore. Furthermore, the EFG shape functions have a highly irregular shape when singular weight functions are applied, so that a high-order integration rule must be used (Hegen, 1996). However, interpolating shape functions (and, thus, singular weight functions) are only needed on the Dirichlet boundary; therefore, the effect of these drawbacks is limited (Chen and Wang, 2000).

Coupling between FE and EFG methods. A radically different approach is to split the solution into an interpolating part on the Dirichlet boundary and a part on which no restrictions are put for the remainder of the domain. In the context of EFG methods, it is convenient to use a row of FEs along the Dirichlet boundary, while the rest of the domain is discretized with EFG shape functions (Krongauz, 1996; Hegen, 1996; Krongauz and Belytschko, 1996; Huerta and Fernández-Méndez, 2000; Wagner and Liu, 2001). A gradual transition from the FE shape functions to the EFG shape functions ensures overall continuity and restricts the influence of EFG nodes close to the boundary.

This concept of coupling the FE method and the EFG method should not be confused with the case that EFG shape functions are added to a FE discretization in order to enrich the spatial resolution. In reference (Huerta and Fernández-Méndez, 2000) both concepts are compared and treated in a unified manner.

If interpolating shape functions are constructed, at least at the Dirichlet boundary, essential boundary conditions can be imposed by means of standard collocation.

As can be seen from the attention focused in (recent) literature, a definitive technique to enforce constraint equations in the EFG method has not yet been established. The enforcement of essential boundary conditions (or, more generally, of constraint equations) is subject of ongoing research.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B0080437494030081

Plate Bending Approximation: Thin and Thick Plates

O.C. Zienkiewicz , ... J.Z. Zhu , in The Finite Element Method: its Basis and Fundamentals (Seventh Edition), 2013

13.3.4 Continuity requirement for shape functions ( C 1 continuity)

The weak form given above contains all the second derivatives of w , thus, it is clear that at least all the first derivative of w must be continuous in Ω . As noted above this class of functions is called C 1 continuous. If the first derivatives are not continuous the plate has "kinks" resulting in singular second derivatives (similar to a Kronnecker delta along a line) and, thus, the square of a singular function is not integrable in Ω .

We have shown above that using Hermite interpolation it is not difficult to construct piecewise polynomials along a line that are C 1 continuous. We also noted that Hermite interpolation is not isoparametric and this causes difficulty when we attempt to extend the use of Hermite-based shape functions to general shaped elements for a thin plate formulation.

To ensure the continuity of both w and its normal slope across an interface we must have both w and w / n uniquely defined by values of nodal parameters along such an interface. Consider Fig. 13.7 depicting the side 1–2 of a rectangular element. The normal direction n is in fact that of y and we desire w and w / y to be uniquely determined by values of w , w / x , w / y at the nodes lying along this line.

Figure 13.7. Continuity requirement for normal slopes.

Figure 13.8. A rectangular plate element.

To show the C 1 continuity along side 1–2, we write

(13.37) w = A 1 + A 2 x + A 3 y +

and

(13.38) w y = B 1 + B 2 x + B 3 y +

with a number of constants in each expression just sufficient to determine a unique solution for the nodal parameters associated with the line.

Thus, for instance, if only two nodes are present a cubic variation of w should be permissible noting that w / x and w are specified at each node. Similarly, for the two nodes only a linear, or two term, variation of w / y would be permissible. Note, however, that a similar exercise could be performed along the side placed in the y direction preserving continuity of w / x along this. Along side 1–2 we thus have w / y , depending on nodal parameters of line 1–2 only, and along side 1–3 we have w / x , depending on nodal parameters of line 1–3 only. Differentiating the first with respect to x , on line 1–2 we have 2 w / x y , depending on nodal parameters of line 1–2 only, and similarly, on line 1–3 we have 2 w / y x , depending on nodal parameters of line 1–3 only.

At the common point, 1, an inconsistency arises immediately as we cannot automatically have there the necessary identity for continuous functions

(13.39) 2 w x y 2 w y x

for arbitrary values of the parameters at nodes 2 and 3. It is thus impossible to specify simple polynomial expressions for shape functions ensuring full compatibility when only w and its slopes are prescribed at corner nodes [16].

Thus if any functions satisfying the compatibility are found with the three nodal variables, they must be such that at corner nodes these functions are not continuously differentiable and the cross-derivative is not unique. The above proof has been given for a rectangular element. Clearly, the arguments can be extended for any two arbitrary directions of interface at the corner node 1. Some elements for use in thin plates with corner discontinuous functions have been developed by Clough et al. [17,18], Fraeijs de Veubeke [19,20], Bazeley et al. [21], and Irons [22].

A way out of this difficulty appears to be obvious. We could specify the cross-derivative as one of the nodal parameters. For an assembly of rectangular elements, this is convenient and indeed permissible. Simple functions of this type have been developed by Bogner et al. [23] and used with success in problems where square elements are adequate. For the geometry of a rectangular element as shown in Fig. 13.8, the shape functions may be taken as combinations of the Hermite interpolations expressed as

(13.40a) w ˆ ( ξ , η ) = a = 1 4 N a ( ξ , η ) u ̃ a

where

(13.40b) N a = H a w ( ξ ) H a w ( η ) , H a θ ( ξ ) H a w ( η ) , H a w ( ξ ) H a θ ( η ) , H a θ ( ξ ) H a θ ( η )

and

(13.40c) u ̃ a = w ̃ a , w ̃ a x w ̃ a y 2 w ̃ a x y T

and the geometry is given by

(13.41) x = a = 1 4 N a ( ξ , η ) x a

with node and function ordering as shown in Table 13.2.

Table 13.2. Node and Function for Rectangular Element

Nodes
a 1 2 3 4
Function
H ( ξ ) 1 2 2 1
H ( η ) 1 1 2 2

Indeed this element is the most accurate rectangular element available for thin problems. A development of this type of element to include continuity of higher derivatives is simple using Hermite interpolation and is outlined in Ref. [24].

In their undistorted form the element is, as for all rectangles, of very limited applicability. Unfortunately, the extension to nodes at which a number of element interfaces meet with different angles (Fig. 13.9) is not, in general, permissible since the element is not of isoparametric type. Here, the continuity of cross-derivatives in several sets of orthogonal directions implies, in fact, a specification of all second derivatives at a node.

Figure 13.9. Nodes where elements meet in arbitrary directions.

This, however, violates physical requirements if the plate stiffness varies abruptly from element to element, for then equality of moments normal to the interfaces cannot be maintained. This process has been used with some success in homogeneous plate situations [24–31] although Smith and Duncan [24] comment adversely on the effect of imposing such excessive continuities on several orders of higher derivatives. For cases where thickness or material changes occur, it is possible to transform the nodal parameters to be normal and tangential to the interface and relax continuity on the parameters 2 w a / n 2 .

The difficulties of finding compatible displacement functions have led to many attempts at ignoring the complete slope continuity while still continuing with the other necessary criteria. Proceeding perhaps from a naive but intuitive idea that the imposition of slope continuity at nodes only must, in the limit, lead to a complete slope continuity, several successful, "nonconforming," elements have been developed [21,32–46].

The convergence of such elements is not obvious but can be proved either by application of the patch test or by comparison with finite difference algorithms. We have discussed the importance of the patch test extensively in Chapter 8, and for plate problems, its importance in both design and testing of elements is paramount and this test should never be omitted. Indeed, we shall show how a successful triangular element is developed via this analytical interpretation [47–52].

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9781856176330000137

Mathematical preliminaries

J.E. Akin , in Finite Element Analysis with Error Estimators, 2005

2.14 General matrix partitions

The above small example has led to the most general form of the algebraic system that results from satisfying the required integral form: a singular matrix system with more unknowns that equations. That is because we chose to apply the essential boundary conditions last and there is not a unique solution until that is done. The algebraic system can be written in a general matrix form that more clearly defines what must be done to reduce the system to a solvable form by utilizing essential boundary condition values. Note that the system degrees of freedom, D, and the full equations could always be rearranged in the following partitioned matrix form

(2.48) [ S u u S u k S k u S k k ] { D u D k } = { C u C k + P k }

where D u represents the unknown nodal parameters, and D k represents the known essential boundary values of the other parameters. The sub-matrices S uu and S kk are square, whereas Suk and Sku are rectangular, in general. In a finite element formulation all of the coefficients in the S and C matrices are known. The P k term represents that there are usually unknown generalized reactions associated with essential boundary conditions. This means that in general after the essential boundary conditions (D k ) are prescribed the remaining unknowns are D u and P k . Then the net number of unknowns corresponds to the number of equations, but they must be re-arranged before all the remaining unknowns can be computed.

Here, for simplicity, it has been assumed that the equations have been numbered in a manner that places the prescribed parameters (essential boundary conditions) at the end of the system equations. The above matrix relations can be rewritten as

S u u D u + S u k D k = C u S k u D u + S k k D k = C k + P k

so that the unknown nodal parameters are obtained by inverting the non-singular square matrix Suu in the top partitioned rows. That is,

D u = S u u 1 ( C u S u k D k ) ,

Most books on numerical analysis assume that you have reduced the system to the non-singular form given above where the essential conditions, D u , have already been moved to the right hand side. Many authors use examples with null conditions (D k is zero) so the solution is the simplest form, D u = S −1 uu C u . If desired, the values of the necessary reactions, P k , can now be determined from

P k = S k u D u + S k k D k C k

In nonlinear and time dependent applications the reactions can be found from similar calculations. In most applications the reaction data have physical meanings that are important in their own right, or useful in validating the solution. However, this part of the calculations is optional. If one formulates a finite element model that satisfies the essential boundary conditions in advance then the second row of the partitioned system S matrix is usually not generated and one can not recover reaction data.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780750667227500333