The Finite Element Method: An Introduction
Dr. Garth N. Wells
[email protected]
University of Cambridge and Delft University of Technology
c G. N. Wells 2009
[email protected]
Contents 1 Introduction 1.1 Tensor basics . . . . . . . . . . . . . . . 1.2 Vector calculus . . . . . . . . . . . . . . 1.3 Linear algebra . . . . . . . . . . . . . . 1.4 Essentials from continuum mechanics 1.5 Exercises . . . . . . . . . . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
7 7 13 15 16 19
2 Strong and weak forms of the governing equations 2.1 One-dimensional bar . . . . . . . . . . . . . . . 2.2 Continuum elasticity . . . . . . . . . . . . . . . 2.3 Poisson equation . . . . . . . . . . . . . . . . . 2.4 Minimisation of potential energy . . . . . . . . 2.5 Regularity requirements . . . . . . . . . . . . . 2.6 Definition of trial and weight function spaces . 2.7 Exercises . . . . . . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
21 21 23 25 26 27 28 29
3 The 3.1 3.2 3.3 3.4
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
31 31 32 33 34
4 Formulation of the finite element method 4.1 Finite element method with piecewise linear basis functions 4.2 General finite element basis functions . . . . . . . . . . . . . 4.3 Governing equations . . . . . . . . . . . . . . . . . . . . . . . 4.4 Isoparametric mapping . . . . . . . . . . . . . . . . . . . . . . 4.5 Numerical integration . . . . . . . . . . . . . . . . . . . . . . 4.6 Imposition of Dirichlet boundary conditions . . . . . . . . . 4.7 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
35 35 37 49 50 54 57 60
5 Implementation of the finite element 5.1 Preprocessing . . . . . . . . . . . 5.2 Program flow . . . . . . . . . . . 5.3 Local-global . . . . . . . . . . . . 5.4 Stiffness matrix storage . . . . . . 5.5 Post-processing . . . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
63 63 64 65 70 71
Galerkin method Approximate solution . . . . . . . . Basic error analysis . . . . . . . . . . Convergence of the Galerkin method Exercises . . . . . . . . . . . . . . . .
. . . .
. . . . .
. . . .
. . . . .
. . . .
method . . . . . . . . . . . . . . . . . . . . . . . . .
6 Structural elements for finite element analysis
. . . . .
. . . .
. . . . .
. . . . .
. . . .
. . . . .
. . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
73
5
Contents 6.1 6.2 6.3 6.4 6.5
Rod elements in space Beams . . . . . . . . . . Plate . . . . . . . . . . Shell elements . . . . . Exercises . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. 73 . 76 . 91 . 104 . 104
7 Analysis of the finite element method 7.1 A priori error estimation . . . . . . 7.2 A posteriori error estimation . . . . 7.3 Adaptivity . . . . . . . . . . . . . . 7.4 Exercises . . . . . . . . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
105 105 107 108 109
8 Solvers for large linear systems 111 8.1 LU-decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 8.2 Iterative solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 9 Time-dependent problems 9.1 Elastodynamics . . . . . . . . . . . . . 9.2 Heat equation . . . . . . . . . . . . . . 9.3 Frequency analysis for elastodynamics 9.4 Modal analysis for elastodynamics . . 9.5 Time stepping algorithms . . . . . . . 9.6 Space-time formulations . . . . . . . . 9.7 Exercises . . . . . . . . . . . . . . . . . References
6
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
113 113 116 117 118 119 124 124 125
Version 0.2
January 21, 2011
1 Introduction This course provides an introduction to the finite element method for linear problems. The underlying mathematics of the method are introduced and details of its computer implementation are discussed. Students will develop an understanding of the fundamentals underlying the finite element method and commercial finite element software. This course also provides a foundation to the more advanced courses, such as ‘Finite element methods for nonlinear analysis’ (CT5142). This chapter establishes the notation used in these notes and reviews the basic tools needed to develop and understand the finite element method. The material should serve as a review, and possibly an extension of familiar concepts. It also provides a reference point for later developments in the text. What is the finite element method? The finite element method (FEM) is a numerical method for solving partial differential equations. It is particularly suited for solving partial differential equations on complex geometries, as are commonly encountered in solid and structural mechanics. The procedure can be largely automatised, making it well suited to efficient computer implementation. The application of the finite element method leads to systems (often very large) of linear equations (matrices) which can be solved using a computer. The finite element method differs from the commonly used finite difference method as it is based on a variational formulation. When applied for solving problems in solid and structural mechanics, it is closely related to the concepts of energy and virtual work.
1.1 Tensor basics Throughout these notes, a Cartesian, orthonormal basis is assumed. This simplifies vector and tensor operations considerably as the basis vectors can be ignored and operations expressed in terms of indices only. Orthonormal basis vectors are shown in Figure 1.1. Scalars are denoted by italic type face, a, b, s.
(1.1)
Vectors involve components and a basis, which for simplicity is assumed to be orthonormal, and are generally denoted by lowercase bold characters, a, b, f , u, x.
(1.2)
7
1 Introduction x2 a e2 = [0 1 0]
e1 = [1 0 0] x1
e3 = [0 0 1] x3 Figure 1.1: Example of orthonormal base vectors. A vector a in an n-dimensional space R n can be expressed in terms of its components, ai where i runs from one to n. In the case n = 2, a (1.3) a= 1 . a2 Consider the dot product between two vectors a and b, which is denoted a · b, and is defined through the operation which gives the length k xk of a vector x,
k xk =
√
x · x.
(1.4)
The dot product of two vectors is given by: a · b = a i bi ,
(1.5)
where repeated indices imply summation: n
a·b =
∑ a i bi .
(1.6)
i =1
For n = 3, 3
a·b =
8
∑ a i b i = a 1 b1 + a 2 b2 + a 3 b3 .
(1.7)
i =1
Version 0.2
January 21, 2011
1.1 Tensor basics For an orthonormal basis, ei · e j = δij , allowing the base vectors to be discarded and the more familiar notation a · b = ai b j δij = ai bi ,
(1.8)
can be used. At times, the dot product maybe written as a T b, in the vein of matrix multiplication. The symbol δij is known as the Kronecker delta: ( 1 if i = j, δij = (1.9) 0 if i 6= j. Second-order tensors (which possess some similar properties to matrices) are generally denoted by uppercase bold characters, (1.10)
A, K, M, R
As with vectors, a second-order tensor A can be expressed in terms of components, Aij , A = Aij .
(1.11)
The second-order tensor A in R3 can be expressed as A11 A12 A13 A = Aij = A21 A22 A23 . A31 A32 A33
(1.12)
Second-order tensors are linear operators that act upon vectors. A second-order tensor A takes the vector b and transforms it into the vector a, n
a = Ab = Aij b j =
∑ Aij bj
(1.13)
j =1
where n is the dimension of the vector b. Consider the case a = ABc = Dc.
(1.14)
Second-order tensors can be multiplied, which is expressed as D = AB,
(1.15)
and in index notation as Dij = Aik Bkj .
(1.16)
The transpose of a tensor is denoted by the superscript ‘T’. It is defined by: a · Bc = Bc · a = c · B T a.
January 21, 2011
(1.17)
Version 0.2
9
1 Introduction In terms of indices: Aij T = A ji .
(1.18)
The transpose of a second-order tensor implies interchanging the rows and columns of its components. A useful identity involving the transpose is:
( AB) T = B T A T .
(1.19)
Second-order tensors can be formed through the dyadic production of two vectors, A = a ⊗ b, = ai b j ,
(1.20)
which is equivalent to Aij = ai b j .
(1.21)
The trace of a second-order tensor is essentially the summation of its diagonal components. For example, tr ( a ⊗ b) = ai bi
(1.22)
tr ( A) = Aii .
(1.23)
and
For a second-order tensor A in R3 , this implies: tr ( A) = A11 + A22 + A33 .
(1.24)
Using the trace, an operation similar to the dot product of two vectors, an inner product, is provided for two second-order tensors. It also yields a scalar. It is defined by: A : B = tr A T B . (1.25) This is equivalent to
A : B = Aij Bij .
(1.26)
The repeated i and j indices implies: n
Aij Bij =
n
∑ ∑ Aij Bij .
(1.27)
i =1 j =1
This operation has no equivalent in matrix-vector convention. For two second-order tensors in R2 , A : B = A11 B11 + A12 B12 + A21 B21 + A22 B22 .
10
Version 0.2
(1.28)
January 21, 2011
1.1 Tensor basics The identity tensor I is defined by: a = Ia,
(1.29)
and is equal to I = δij ei ⊗ e j .
(1.30)
A second-order tensor A, multiplied by its inverse A−1 , is equal to the identity tensor, AA−1 = A−1 A = I.
(1.31)
Second-order tensors are often characterised by their special properties. A matrix A is symmetric if: A T = A.
(1.32)
Using index notation, a tensor Aij is symmetric if: Aij = A ji .
(1.33)
A tensor A is orthogonal if: A T = A −1 .
(1.34)
As will be shown in the following section, matrices which rotate the reference frame are orthogonal. Another important operation is the cross, or vector, product. It is defined for a, b ∈ R3 by: a 2 b3 − a 3 b2 (1.35) [ a × b ] = a 3 b1 − a 1 b3 . a 1 b2 − a 2 b1
Such a definition however is not convenient for manipulation. A third-order tensor can be defined which will make manipulations more convenient. A third-order tensor maps a second-order tensor to a first order tensor. A common third-order tensor is the alternating (or permutation) tensor E . The components of E are given by: (1.36) [E ] = Eijk = ei · e j × ek . The expression for E in terms of the basis ei is lengthy and not shown here. The result of the above equation is rather simple. If any of the indices are repeated, Eijk = 0. If the permutation {i, j, k} is even, Eijk = 1, and if the permutation {i, j, k} is odd, Eijk = −1. The cross product of two vectors yields a vector. The cross product can be introduced via the third-order tensor E ,
E : ( a ⊗ b) = a × b. January 21, 2011
(1.37)
Version 0.2
11
1 Introduction x2
e2′ = [− sin θ cos θ ]
e2 = [0 1] e1′ = [cos θ sin θ ] θ
e1 = [1 0] x1
Figure 1.2: Rotation between coordinate systems. This operation is particularly important when considering moments. For an orthonormal basis, in terms of indices,
A : B = Aijk Bjk .
(1.38)
In elasticity, the stress and strain are second-order tensors. They are related to each other via a fourth-order tensor, denoted here as C .
C = Cijkl .
(1.39)
Fortunately, due to physical symmetries, it will not be necessary to manipulate fourthorder tensors directly. Coordinate transformations For many problems, it is convenient to rotate the coordinate frame. Vector and tensor components can be rotated (change of basis) using a rotation matrix Q, which is orthogonal. Consider the orthonormal coordinate systems in Figure 1.2. It may prove convenient to transform vector components relative to the ei system to components in the ei′ system. Consider that for a vector a, the components in ei system are known, and it would be convenient to have the components in the ei′ system. Defining Qij by Qij = ei · e′j ,
(1.40)
the components ai′ are given by ai′ = Q ji a j ,
12
(1.41)
Version 0.2
January 21, 2011
1.2 Vector calculus which can also be expressed as: a′ = Q T a.
(1.42)
Considering the coordinate systems illustrated in Figure 1.2,
e1 · e2′ cos θ = sin θ e2 · e2′
e1 · e1′ Q= e2 · e1′
− sin θ , cos θ
(1.43)
and note that Q is orthogonal, QQ T = Q T Q = I.
(1.44)
Therefore, components in the rotated reference frame can be transformed back to the original reference frame by: a = Qa′ .
(1.45)
At times it is necessary to also rotate the components of a second-order tensor from one coordinate system to another. This is done by the operation: A′ = Q T AQ,
(1.46)
A = QA′ Q T .
(1.47)
and
Often coordinate transforms are used to simplify a problem. An example is a onedimensional rod element in a three-dimensional space. It is convenient to rotate the reference such such the position along the bar is given by [ x]′ = x1′ , 0, 0 , which effectively reduces the three-dimensional problem to a one-dimensional problem.
1.2 Vector calculus Subscripts are commonly used as a short-hand notation for derivatives. Characters after a comma in the subscript denote differentiation with respect to that variable. Consider a function u which is dependent on the position x. Some of its derivatives can be written as: ∂u = u,x , ∂x
∂4 u ∂2 u = u , = u,xxxx , ,xx ∂x2 ∂x4
∂2 u ∂3 u = u,xyy . = u,xy , ∂x∂y ∂x∂y2
(1.48)
This notation is used commonly throughout these notes. An important operator is the divergence. The divergence of a vector field a yields a scalar. It is expressed as:
∇ · a. January 21, 2011
(1.49)
Version 0.2
13
1 Introduction In terms of indices, ∇ · a is equal to
∇·a =
∂ai = ai,i . ∂xi
(1.50)
In a three-dimensional space R3 ,
∇·a =
∂ay ∂a ∂a x ∂a ∂az ∂a1 + 2+ 3 = + + . ∂x1 ∂x2 ∂x3 ∂x ∂y ∂z
(1.51)
The divergence of a second-order tensor is given by:
∇·A =
∂Aij = Aij,j . ∂x j
(1.52)
For a two-dimensional tensor A, ∂A11 ∂A12 + ∂x ∂x2 . ∇ · A = ∂A 1 ∂A22 21 ∂x + ∂x
(1.53)
∂a = a,i , ∂xi
(1.54)
1
2
Another important operator is the gradient. The gradient of a scalar a is given by:
∇a =
which yields a vector. In a three-dimensional space, the components of ∇ a are given by: ∂a ∂x ∂a (1.55) ∇a = ∂y . ∂a ∂z The gradient of a vector is given by:
∇a =
∂ai = ai,j , ∂x j
(1.56)
which yields a second-order tensor. In a three-dimensional space, the components of ∇ a are given by: ∂a1 ∂a1 ∂a1 ∂x ∂y ∂z ∂a2 ∂a2 ∂a2 . ∇a = (1.57) ∂y ∂z ∂x ∂a3 ∂a3 ∂a3 ∂x ∂y ∂z 14
Version 0.2
January 21, 2011
1.3 Linear algebra
n Ω
Γ Figure 1.3: Continuous body Ω ⊂ R n bounded by Γ. An essential theorem in formulating the finite element method is the divergence theorem (also known as Gauss’s theorem). It transforms a volume integral to a surface integral. Consider the body Ω in Figure 1.3. The boundary of the body Ω is denoted Γ = ∂Ω. The outward unit normal to the body is denoted n. The divergence theorem states: Z
Ω
∇ · A dΩ =
Z
∂Ω
An dΓ,
(1.58)
where n is the outward normal to the body Ω. The surface of Ω is given by ∂Ω. For a vector a, Z
Ω
∇ · a dΩ =
Z
∂Ω
a · n dΓ.
(1.59)
Integration by parts is used frequently in the formulation of the finite element method. It implies that: Z
Ω
a · (∇ · B) dΩ = −
Z
Ω
∇ a : B dΩ +
Z
∂Ω
a · Bn dΓ.
(1.60)
The above equation also holds for the case in which a is replaced by a scalar and B is replaced by a vector. The above formula is simple to derive by using the product rule for differentiation to expand Z
Ω
∇ · ( aB ) dΩ =
Z
Ω
∂ a B dΩ, ∂x j i ij
(1.61)
and then applying the divergence theorem.
1.3 Linear algebra The eigenvalues and eigenvectors of a matrix provide information as to its properties. For example, the eigenvalues of the stress tensor are the principal stresses, and the
January 21, 2011
Version 0.2
15
1 Introduction corresponding eigenvectors are the principal stress directions. Also, the eigenvalues of a matrix indicates if a matrix has a unique inverse. For a matrix A, vectors b and scalars λ which satisfy
( A − λI ) b = 0
(1.62)
are known as eigenvectors and eigenvalues, respectively, of A. For a symmetric matrix, all eigenvalues are real. For a positive definite matrix A, all eigenvalues are positive and real, which guarantees that the inverse of A exists, and is unique. Equation (1.62) implies that: det ( A − λI ) = 0.
(1.63)
For a 2 × 2 matrix, finding λ requires finding the roots of a quadratic equation. For a 3 × 3 matrix, the roots of a cubic equation must be found. The polynomial is known as the characteristic equation of A. For large matrices, special numerical algorithms are used to calculate the eigenvalues. A zero eigenvalue implies that the matrix A has a linearly dependent column. Eigenvalue tests are often used in studying the matrices arising in the finite element method. A system of n linear equations can be expressed as: Ku = f ,
(1.64)
where K is an n × n matrix, and u and f are vectors of length n. In the context of the finite element method, K will often be referred to as the ‘stiffness matrix’ and f the ‘right-hand side vector’. The solution u is given by: u = K −1 f .
(1.65)
Typically, the system of equations in a finite element equation will be very large. It is not practical to compute the inverse of the matrix K as it is computationally expensive (both in terms of the number of operations required and the memory required). The system of equations is usually solved using a direct method (such as Gauss elimination or LU decomposition) or an approximate iterative method. Solvers for finite element problems are discussed in Chapter 8.
1.4 Essentials from continuum mechanics The finite element method is used in these notes to solve primarily problems in continuum mechanics. This requires some background in continuum mechanics and linear-elasticity. This section is intended to cover only the important concepts which are used later in these notes. It is not a comprehensive coverage of the topic. The strain ǫ is defined as the symmetric gradient of the displacement vector: ǫij =
16
1 ui,j + u j,i , 2
(1.66)
Version 0.2
January 21, 2011
1.4 Essentials from continuum mechanics which is clearly a second-order tensor. In compact notation, the above expression can be written as: 1 ǫ = ∇s u = (1.67) ∇u + (∇u) T . 2 For convenience, the symmetric gradient using index notation will be expressed as: u(i,j) =
1 ui,j + u j,i , 2
(1.68)
where the bracket around the subscript denotes the symmetric gradient. Obviously, the strain tensor is symmetric. The stress tensor is defined through the traction vector t (force per unit area) on a surface, σn = t
(1.69)
where n is the unit normal vector to the surface. Similar to the strain tensor, the stress tensor is symmetric (which can be proved by considering the balance of angular momentum, as will be shown in the next chapter). A constitutive equation relates the stresses in a material to the strain. For solid, linear-elastic continua, this can be written as: σ = C : ǫ,
(1.70)
or in index notation, σij = Cijkl ǫkl
(1.71)
where C is a fourth-order tensor. For linear elasticity, it is defined by: Cijkl = µ δik δjl + δil δjk + λδij δkl
(1.72)
where
λ=
νE (1 + ν) (1 − 2ν)
(1.73)
µ=
E 2 (1 + ν )
(1.74)
and
with E the Young’s modulus and ν Poisson’s ratio. The constants λ and µ are often referred to as Lam´e constants. In a homogeneous body, Lam´e constants are independent of the position. Due to symmetry, in three-dimensions the stress and strain tensors have only six independent components. In two-dimensions, this reduces to four. For convenience,
January 21, 2011
Version 0.2
17
1 Introduction ‘engineering’ notation is often adopted, in which the stress and strain are represented as vectors, u,x ǫ11 ǫ11 v,y ǫ22 ǫ22 u,z ǫ33 ǫ33 . (1.75) = = ǫ= u,y + v,x γ12 2ǫ12 u,z + w,x γ13 2ǫ13 v,z + w,y γ23 2ǫ23
The stress tensor is expressed in a vector form, although without the factor ‘2’ on the shear terms. σ11 σ22 σ33 (1.76) σ= σ12 σ13 σ23 Using the vector format for stress and strain, the constitutive relationship can be written as: σ = Dǫ
(1.77)
where the matrix D is the constitutive matrix, containing components of the fourthorder tensor C . For isotropic linear-elasticity, D has the following form: λ + 2µ λ λ 0 0 0 λ λ + 2µ λ 0 0 0 λ λ λ + 2µ 0 0 0 (1.78) D= 0 0 0 µ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ
When reducing a three-dimensional problem to two-dimensions, there are three possibilities. The first is known as plane strain. This assumes that the normal strain in the out-of-plane direction is zero (ǫ33 = 0). It is then likely that the stress in the out-of-plane direction is not equal to zero. This is reasonable for problems which are relatively thick in the third direction, such as a dam wall. The second possibility is plane stress (σ33 = 0). In this case, the normal stress out-of-plane is assumed to be zero. This approximation is good for very thin members, such as many beams. The third possibility is axisymmetric. In this case, the strain in the out-of-plane direction is dependent on the distance from an axis of rotation. An axisymmetric formulation is for example commonly used when simulating a circular foundation. In these notes, particular attention will be paid to solving elastic continuum problems. This will require identification of the equations which govern equilibrium of
18
Version 0.2
January 21, 2011
1.5 Exercises elastic bodies. In order to derive the governing partial differential equation, consider first translational equilibrium (balance of linear momentum) of a body Ω ⊂ R n with boundary Γ = ∂Ω, which requires that Z
∂Ω
t dΓ +
Z
Ω
b dΩ =
Z
∂Ω
σn dΓ +
Z
Ω
b dΩ = 0
(1.79)
where b is a body force. Application of the divergence theorem to the integral over ∂Ω gives: Z
Ω
∇ · σ + b dΩ = 0
(1.80)
Balance of linear momentum should also be satisfied for any subdomain of Ω. Therefore, the integral can be removed from equation(1.80), leading to:
∇ · σ + b = 0 in Ω
(1.81)
Consider now moment equilibrium (balance of angular momentum). Moment equilibrium of the body Ω requires that: Z
∂Ω
r × t dΓ +
Z
Ω
r × b dΩ =
Z
∂Ω
r × σn dΓ +
Z
Ω
r × b dΩ = 0
(1.82)
where r is the position vector of a point on ∂Ω. Considering that t = σn, and applying the divergence theorem, Z
Ω
r × (∇ · σ + b) dΩ +
Z
Ω
E : σ T dΩ = 0
(1.83)
If translational equilibrium (equation (1.80)) is satisfied, the first term in equation (1.83) is zero. Then, moment equilibrium is satisfied if the stress tensor is symmetric, since E : σ T = 0, if σ is symmetric. Therefore, the partial differential equation of equilibrium is given by equation (1.80). This will also be referred to as the ‘strong’ governing equation. Note that this equation is general, and not specific to elasticity.
1.5 Exercises 1. For three-dimensions (i = 1 → 3, etc.), expand the terms: a) b) c) d) e) f)
a i bi ai b j Aij Bij Cij bi Cij b j Fik Gkj
2. How many independent components does a symmetric second-order tensor in R n have?
January 21, 2011
Version 0.2
19
1 Introduction 3. What is the trace of the identity tensor I equal to in R n ? 4. Confirm that the rotation matrix R in equation (1.43) is orthogonal. 5. Using the product rule for differentiation and the divergence theorem, derive the integration by parts formula in equation (1.60). 6. A common and important operator in applied mathematics (and mechanics) is the Laplace operator ∆ (sometimes denoted ∇2 ). It is equivalent to ∇ · ∇. Expand ∆u in a three-dimensional space using index notation. 7. Derive the isotropic linear-elastic constitutive matrix D for plane stress and plane strain conditions
20
Version 0.2
January 21, 2011
2 Strong and weak forms of the governing equations An essential step in developing the finite element method is the identification of the governing equation and casting it in its ‘weak form’. To do this, the strong form is multiplied by a weight function and integrated by parts. This has the effect of reducing the order of the derivatives appearing in the equation, and leads to a form which is convenient for later numerical solution. The weak form is also commonly known as a the ‘variational form’ or a ‘variational equation’. The finite element method is a variational method, addressing the weak (variational) form of the relevant governing equation.
2.1 One-dimensional bar Consider the one-dimensional bar in Figure 2.1. The governing equation for a onedimensional rod element with unit cross-sectional area is (see equation (1.81)):
−σ,x = f ,
(2.1)
where σ,x is the normal stress in the rod and f is a distributed load along the rod. This is complemented by a constitutive relationship which gives the stress in the rod in terms of the strain. For a linear-elastic rod, the normal stress is given by: σ = Eǫ = E
du , dx
(2.2)
where E is the Young’s modulus. The governing equation can then be expressed as:
− Eu,xx = f .
(2.3)
To complete the problem, boundary conditions must be specified. A Dirichlet (essential) boundary condition prescribes the displacement at a point. For example, u = 0 at
x=0
(2.4)
0110 111 000 101010 0 1 0 111111111111 00000000000 f
h
x
L
x=0
Figure 2.1: One-dimensional bar.
21
2 Strong and weak forms of the governing equations prescribes the displacement at one end of the bar in Figure 2.1 to be equal to zero. A Neumann (natural) boundary condition prescribes the applied traction, σn = h
at x = L,
(2.5)
where n is the outward normal to the bar (equal to 1 in this case). The above boundary condition sets the applied traction at x = L equal to h. Considering that σ = Eu,x , applying a force at the end of the bar is equivalent to prescribing u,x . The mathematical problem can now be summarised as: find u such that:
− Eu,xx = f u=0 Eu,x n = h
0 < x < L, at x = 0, at x = L.
(2.6) (2.7) (2.8)
To develop the weak form of the governing equation, both sides of equation (2.6) are multiplied by a scalar weight function w ∈ V (w ∈ V means that w comes from the appropriately defined space of functions V , which has an infinite number of members). An important requirement on the space V is that functions w be zero where Dirichlet boundary conditions are applied. Since both sides of the equation are multiplied by the weight function, the equation must still hold for all admissible weight function (this is a technical point which is discussed in section 2.6). Multiplying equation (2.1) by w,
−wσ,x = w f
∀w ∈ V ,
(2.9)
where ‘∀w ∈ V ’ means ‘for all w ∈ V ’. If equation (2.9) must hold point-wise, then
−
Z L 0
wσ,x dx =
Z L 0
w f dx
∀w ∈ V ,
(2.10)
is also true. Taking guidance from equation (1.60), equation (2.10) can be integrated by parts and the Neumann boundary condition from equation (2.8) inserted, which yields:
−
Z L 0
wσ,x dx =
Z L 0
w,x σ dx − wh| x= L
∀w ∈ V .
(2.11)
Note the term wEu,x n| x=0 is not included as by definition the weight function w is equal to zero at x = 0. Inserting now the constitutive relationship σ = Eu,x , solving the weak form involves: find u ∈ S such that Z L 0
w,x Eu,x dx =
Z L 0
w f dx + wh| x= L
∀w ∈ V ,
(2.12)
where S is an appropriately defined space of functions which satisfies the Dirichlet boundary conditions. This is the weak form of the governing equation. Note that in the strong form, after inserting the constitutive relationship (2.3), derivatives of order
22
Version 0.2
January 21, 2011
2.2 Continuum elasticity two with respect to x appear. In the weak form (2.12), only first derivatives with respect to x appear. Since the weak form is derived from the strong form, it follows that a solution of the strong form is also a solution of the weak form. What remains to be shown is that the weak form is equivalent to the strong form (that a solution of the weak form is also a solution of the strong form, at least in a generalised sense). Assuming that E is constant and integrating by parts the first term in (2.12) yields
−
Z L 0
wEu,xx dx + wEu,x n| x= L =
Z L 0
w f dx + wh| x= L .
(2.13)
Note again that no terms at x = 0 appear, as by construction, w(0) = 0. Following Hughes (1987), consider a weight function w ∈ V which is equal to φ ( Eu,xx + f ), where φ is greater than zero over the interval ]0, L[ and zero outside (specifically at x = 0 and x = L). Inserting the weight function into equation (2.13),
−
Z L 0
φ ( Eu,xx + f ) Eu,xx dx =
Z L 0
φ ( Eu,xx + f ) f dx,
(2.14)
which can be rearranged to give Z L 0
φ ( Eu,xx + f )2 dx = 0.
(2.15)
Since φ > 0 and ( Eu,xx + f )2 ≥ 0, the above equation can only hold if − Eu,xx = f , satisfying the strong governing equation (in a distributional sense). What remains now is to ensure that the Neumann boundary condition is satisfied (Dirichlet boundary conditions are satisfied by construction). Returning to equation (2.13), it is now proven that the first term is equal to zero. Hence, for any allowable weight function at x = L, Eu,x n = h,
(2.16)
proving that the Neumann boundary conditions are satisfied. For a more elaborate mathematical treatment, see Brenner and Scott (1994).
2.2 Continuum elasticity A more complicated example is continuum elasticity. The basic form of the equations is the same as the one-dimensional bar in the previous section. The complication arises from the three-dimensional setting. Consider now the continuous body in Figure 2.2. The strong governing equation for this problem was developed in the previous chapter. Recall that the governing equation is given by:
∇ · σ + b = 0 in Ω. January 21, 2011
(2.17)
Version 0.2
23
2 Strong and weak forms of the governing equations
h
Γh Ω
Γg
n
Figure 2.2: Continuous body Ω ⊂ R n with boundary Γ (Γ = Γg ∪ Γh , Γg ∩ Γh = ∅). As in the one-dimensional case, a constitutive relationship is required. For a linearelastic material, σ = C : ǫ,
(2.18)
where C is a fourth-order tensor. Boundary conditions are still required to complete the problem. On part of the boundary, the displacements are prescribed (Dirichlet boundary conditions). In Figure 2.2, the part of the boundary where displacements are prescribed is denoted Γg , u=g
on Γg .
(2.19)
On the remaining boundary, Γh , tractions (potentially zero) are applied (traction ≡ force per unit area on a surface), t=h
on Γh .
(2.20)
The boundary-value problem is now complete. To derive the weak form, the same steps as in the case of the one-dimensional bar are followed. First, the governing equation is multiplied by a weight function w ∈ V , and integrated over the body. Note that in this case, the weight function w is a vector. Multiplying equation (2.17) by w ∈ V , where again V is an appropriately defined function space (w = 0 on Γg ), and integrating over the body Ω yields: Z
Ω
w · (∇ · σ ) dΩ +
Z
Ω
w · b dΩ = 0
∀w ∈ V .
(2.21)
The first term in equation (2.21) can be integrated by parts (see equation (1.60)), yielding:
−
24
Z
Ω
∇w : σ dΩ +
Z
Ω
w · b dΩ +
Z
Γh
w · h dΓ = 0,
Version 0.2
(2.22)
January 21, 2011
2.3 Poisson equation where σn has been replaced by the prescribed traction h (the Neumann boundary condition). Inserting now the constitutive relationship from equation (2.18), a solution to the weak form involves: find u ∈ S such that
−
Z
Ω
∇w : C : ǫ dΩ +
Z
Ω
w · b dΩ +
Z
Γh
w · h dΓ = 0
∀w ∈ V ,
(2.23)
where S is an appropriate space of functions which satisfy the Dirichlet boundary conditions. This is the weak equilibrium equation for a stationary continuous elastic body. Consistency can be proven in the same fashion as was used for the onedimensional problem. Often the term ∇w is written as ∇s w, the symmetric gradient of w. It is simple to show that ∇w : σ = ∇s w : σ if σ is symmetric (see exercise 1). Consider now w as a ‘virtual displacement’ δu. Therefore ∇s w ≡ δǫ. Inserting this into equation (2.23), Z
Ω
δǫ : C : ǫ dΩ =
Z
Ω
δu · b dΩ +
Z
Γh
δu · h dΓ,
(2.24)
which is the equation of virtual work. The equation of virtual work is therefore the weak form of the governing equation. There exists a close link between weak forms and virtual work for a range of problems. However, the variational framework is more general.
2.3 Poisson equation An equation which is often solved using the finite element method is the Poisson equation. A wide range of physical phenomena are governed by this equation, including steady-state heat conduction and flow in permeable media (Darcy’s law). The previous example of continuum elasticity was an example of a system of Poisson equations – one for each spatial direction. Here the Poisson equation is addressed in a more generic fashion. Consider the following equation:
−∇ · q + f = 0 in Ω
(2.25)
where q is a flux vector and f is a source term. In the case of heat conductivity, q is the heat flux vector. In Darcy’s law, q is the flow rate. The constitutive relationship is given by: q = −κ ∇ u
(qi = −κij u,j )
(2.26)
where u is a potential. Inserting the constitutive equation into (2.25) yields a Poisson equation. The constitutive relationship depends on the problem being solved. For linear heat conduction, the scalar u is the temperature. For an isotropic medium, κ = κI, where κ is the thermal conductivity. For Darcy’s law, u is the hydraulic head and κ is known as the hydraulic conductivity.
January 21, 2011
Version 0.2
25
2 Strong and weak forms of the governing equations As for all previous examples, the boundary-value problem requires boundary conditions. The Dirichlet conditions impose: u=g
on Γg
(2.27)
which prescribes the potential on the boundary. For heat conduction, this is the temperature on the boundary Γg . For Darcy flow, it is the hydraulic head. The Neumann boundary condition requires that:
−q · n = h on Γh
(2.28)
where n is the outward normal to the surface Γh (see Figure 2.2). In the case of heat conduction, this imposes the heat flux on Γh , and for Darcy flow is it the inflow/outflow across a boundary. The weak form of the Poisson equation is derived following the same steps as in the previous two examples. Multiplying equation (2.25) by a weight function w and integrating over the body Ω,
−
Z
Ω
w (∇ · q) dΩ +
Z
Ω
w f dΩ = 0
∀w ∈ V .
(2.29)
Integrating the LHS of the above equation by parts yields: Z
Ω
∇w · q dΩ −
Z
Γh
wq · n dΓ +
Z
Ω
w f dΩ = 0
∀w ∈ V .
(2.30)
Inserting now the constitutive relationship and the boundary condition −q · n = h leads to the problem: find u ∈ S such that:
−
Z
Ω
∇w · κ∇u dΩ +
Z
Γh
wh dΓ +
Z
Ω
w f dΩ = 0
∀w ∈ V ,
(2.31)
which is the weak form of the Poisson equation.
2.4 Minimisation of potential energy For particular equations, such as equilibrium of an elastic body, there exists a close link between minimisation of energy and solution of the weak form. For an elastic body Ω, suppose that u is the solution to the weak problem, and consider the potential energy functional I, I (u) =
1 2
Z
Ω
∇s u : C : ∇s u dΩ −
Z
Ω
u · b dΩ −
Z
∂Ω
u · h dΩ.
(2.32)
The potential energy for another displacement field, w = u + v, is given by I (u) ≤ I (u + v) =
1 2
Z
Ω
∇s (u + v) : C : ∇s (u + v) dΩ −
26
Z
Ω
(u + v) · b dΩ −
Version 0.2
Z
∂Ω
(u + v) · h dΩ, (2.33)
January 21, 2011
2.5 Regularity requirements
f
f ,x
x
x
Figure 2.3: Example of a C0 function f . which is equal to I (u) ≤ I (u + v) Z Z Z s s = I (u) + ∇ v : C : ∇ u dΩ − v · b dΩ − Ω
+
1 2
Z
Ω
Ω
∂Ω
v · h dΩ
(2.34)
∇s v : C : ∇s v dΩ,
in which the terms inside the brackets sum to zero for R any v since u is a solution to the weak form (see equation 2.23), and the last term Ω ∇s v : C : ∇s v dΩ ≥ 0, which proves that solution of the weak form corresponds to minimisation of the potential energy.
2.5 Regularity requirements A common topic in finite element analysis is continuity. Functions are often classified as being C n continuous. Continuity refers to whether or not the derivatives of a function are continuous. If a function is C n , this means that its nth derivative is continuous. In finite element analysis, functions which are C0 are most commonly used. This means that the functions are continuous, but their first derivatives are not. A C0 function is shown in Figure 2.3. The function shown in Figure 2.3 is a continuous, piecewise linear function. As will be seen in the following chapters, this is a commonly used function in finite element analysis. From Figure 2.3, it is clear that it is simple to take the first derivative of the function f . However, in a classical sense, its second derivative does not exist. Consider now if the function f is a possible solution u of the strong equation for a bar (equation (2.3)). It is not a possible solution, since for equation (2.3) to make sense, the second derivative of f with respect to x ( f ,xx ) must exist. Examining now the weak form for the rod (equation (2.12)), it is clear that it involves only first derivatives with respect to x of u and w. Therefore a function like f is a possible solution of the weak form and a possible weight function. In summary, the weak form allows for more possible solutions. This is a classic mathematical property of weak forms, and is essential for the finite element method.
January 21, 2011
Version 0.2
27
2 Strong and weak forms of the governing equations The theory of which functions are allowed is discussed on more detail in the following advanced section.
2.6 Definition of trial and weight function spaces There are certain mathematical requirements that must be met by the trial and weight functions for the weak form to ‘make sense’. A function u is said to be square integrable if: Z
Ω
(u)2 dx < ∞.
(2.35)
Functions which are square integrable come from the Hilbert space known as L2 (Ω). Clearly, such functions may be discontinuous. When solving second-order differential equations, functions which have square-integrable derivatives are of interest. Functions for which: Z
Ω
(u)2 + (u,x )2 dΩ < ∞
(2.36)
are known as members of the Sobolev space on degree one, which is denoted H 1 (Ω). The notation u ∈ H 1 means that u comes from (is an element of) H 1 . In physical terms, the requirement that the trial functions be square-integrable implies that energy must be finite. Note that C0 functions belong to H 1 . For multiple spatial dimensions, there do however exist functions from H 1 (Ω) which are not C0 continuous. This can be seen in the well-known Sobolev embedding theorem. Importantly, Sobolev spaces are infinite-dimensional. That is, an infinite number of different functions belong to a given Sobolev space. This is of crucial importance, since a fundamental step in the finite element method will involve finite-dimensional function spaces. It is useful to define function spaces, from which the trial and weight functions come. Functions spaces can be considered the ‘family’ of functions from which u and w can come. Consider first a scalar problem in an n-dimensional domain Ω. Trial functions u come from the space S (u ∈ S ). The function space S is defined by:
S = {u | u ∈ H 1 (Ω) , u| Γg = g}.
(2.37)
This says that the S is a collection of functions from H 1 (Ω) which satisfy the Dirichlet boundary conditions (technically this is not a space if g 6= 0, but this point is not important here). The weight functions come from the space V (w ∈ V ), which is defined by:
V = {w | w ∈ H 1 (Ω) , w| Γg = 0}.
(2.38)
This says that V is a collection of functions which come from H 1 (Ω) which are equal to zero where Dirichlet boundary conditions are applied.
28
Version 0.2
January 21, 2011
2.7 Exercises In multiple dimensions, u ∈ S where
S = {ui | ui ∈ H 1 (Ω) , ui | Γg = gi }
(2.39)
and similarly for weight functions w ∈ V ,
V = {wi | wi ∈ H 1 (Ω) , wi | Γg = 0}.
(2.40)
Further details can be found in texts on functional analysis and the mathematical analysis of the finite element method.
2.7 Exercises 1. For 3 × 3 second-order tensors (matrices), prove that ∇w : σ = ∇s w : σ if σ is symmetric. (Hint: use index notation.) 2. Derive the weak equilibrium equation (2.23) for elasticity using index notation. 3. The Helmholtz equation is given by:
∇ · (∇φ) + kφ = 0 and is often used in wave propagation problems. Derive the weak form. 4. Show whether the following functions belong to H 1 on the given domain: a) b) c) d)
u = 1/r on 1 < r < 2 u = 1/r on 0 < r < 2 √ u = 1/ r on 0 < r < 2 u = ax for x < 0 and u = bx for x > 0 on −1 < x < 1 (a and b are arbitrary constants, a 6= b)
5. Solutions for the√displacement field from linear-elastic fracture mechanics often involve a term r, where r is the distance √ from the crack tip (r ≥ 0). If the displacement field was of the form u = a r (all derivatives of a exist and are well behaved and a is not dependent on r), is u a possible trial solution?
January 21, 2011
Version 0.2
29
3 The Galerkin method The finite element method is one particular Galerkin method, named after the Russian engineer Galerkin. It is a method for finding approximate solutions to partial differential equations. It is closely related to the Rayleigh-Ritz method which involves choosing functions (a basis) for the solution and finding the amplitude of each function by minimising the energy. The Galerkin method is however more general, being able to solve a greater range of problems. The essence of the Galerkin method involves taking the weak form of the governing equation, as developed in the previous chapter, and finding the best solution to a problem given a collection of functions.
3.1 Approximate solution First, consider the approximate solution to some problem, uh ∈ S h , where S h ⊂ S is a finite-dimensional space. This means that there is a limited number of possibilities. For example, uh could be a combination of low order polynomial functions. The Galerkin problem for an elastic bar (equation (2.12)) involves: find uh ∈ S h such that
−
Z L 0
h h dx + w h h| x= L = 0 Eu,x w,x
∀w h ∈ V h .
(3.1)
where V h ⊂ V is a finite dimensional space. Commonly, the above equation is expressed in the abstract format as: find uh ∈ S h such that ∀w h ∈ V h (3.2) B wh , uh = L wh where
and
Z B wh , uh =
L 0
h h w,x Eu,x dx
(3.3)
L wh = wh h| x= L .
(3.4)
This abstract format is introduced to keep the derivation of some later developments compact. Furthermore, it is generic for a range of different problems. In a multidimensional context, the finite-dimensional trial and weight functions are denoted uh and wh , respectively. Considering now the elasticity problem in equation (2.23), the Galerkin problem involves: find uh ∈ S h such that
−
Z
Ω
∇s wh : C : ǫh dΩ +
Z
Ω
wh · b dΩ +
Z
Γh
wh · h dΓ = 0
∀wh ∈ V h
(3.5)
31
3 The Galerkin method where ǫ h = ∇s uh . In the abstract notation of equation 3.2, Z ∇s wh : C : ǫh dΩ B wh , uh = Ω
(3.6)
and
Z Z wh · b dΩ + L wh = Ω
Γh
wh · h dΓ.
(3.7)
The Galerkin method (more specifically, the Bubnov-Galerkin method) requires that the weight and trial functions come from the same finite-dimensional space, taking into account the special requirements on the weight and trial functions where Dirichlet boundary conditions are applied. In Petrov-Galerkin method, the weight functions come from a different function space than the trial functions. This method is used for special applications, often in fluid mechanics. Different Galerkin-based methods are defined by how the unknown field uh is represented. In the finite element method, uh and w h will be simple continuous, piecewise low-order polynomials defined on ‘finite elements’. Spectral Galerkin methods for example use a truncated Fourier series as the basis. A basic question which arises when computing an approximate solution is how uh relates to the exact solution u. Given a finite number of possibilities in S h , which solution does the method seek? Understanding this requires some basic error analysis. The error analysis will tell how the computed solution uh differs from the actual solution u and why the Galerkin method works (or for problems not considered here, why it doesn’t work). This is examined in the following section.
3.2 Basic error analysis It is interesting to check how the solution computed using the Galerkin procedure compares to the exact solution. For the one-dimensional problem, the error in the displacement e at a point is defined by: e = u − uh
(3.8)
uh = u − e
(3.9)
where u is the exact solution and uh is the solution to equation (3.1). The approximate solution is therefore equal to:
For generality, this will be investigated using the abstract notation. Inserting equation (3.9) into equation (3.2), B w h , u − B w h , e − L w h = 0 ∀w h ∈ V h (3.10) Since u is the exact solution, B w h , u − L w h = 0 ∀w h ∈ V h , this implies that: (3.11) B w h , e = B w h , u − B w h , uh = 0 ∀w h ∈ V h 32
Version 0.2
January 21, 2011
3.3 Convergence of the Galerkin method which in mathematical terms means that the error is orthogonal to the function space V h with respect to B (·, ·). This important result is commonly known as Galerkin orthogonality. This means that the approximate solution uh is a projection of the exact solution u onto the space of the weight functions. In the Bubnov-Galerkin, the weight functions w h come form the same space as the trial functions uh , hence the solution uh is the projection of the exact solution onto the finite dimensional space of trial functions. It can be shown that the Galerkin finite element method is optimal in terms of the energy. This error analysis tells something of what the Galerkin method calculates. Given some approximate functions, the Galerkin method will yield the best fit to the exact solution in terms of energy. Consider the following: B u − uh + vh , u − uh + vh = B u − uh , u − uh + vh + B vh , u − uh + vh (3.12) = B u − uh , u − uh + 2B vh , u − uh + B vh , vh for any vh ∈ V h , where u is the exact solution and uh is the solution to the Galerkin problem. From Galerkin orthogonality (equation (3.11)), the term B(u − uh , vh ) is equal to zero. Furthermore, B(vh , vh ) ≥ 0. Denoting now w h = uh − vh , the above result leads to the conclusion that: B u − uh , u − uh ≤ B u − wh , u − wh ∀w h ∈ V h (3.13)
This implies that the solution uh is closer to u than any other element of V h in terms of B (·, ·). Consider the ‘energy’ norm:
kvk2E = 12 B (v, v)
(3.14)
which for an elastic body is the strain energy for a given displacement field v. Equation (3.13) can then be expressed as:
ku − uh k E ≤ ku − wh k E
∀w h ∈ V h
(3.15)
which says that the solution computed using the Galerkin method yields a solution which is optimal in terms of the strain energy. Given a choice of functions, the Galerkin method therefore chooses those which minimise the error in terms of the strain energy. There is a close relationship between the Rayleigh-Ritz method and the Galerkin method. It has been shown that the Galerkin method minimises the error in terms of the energy, which is the principle behind the Rayleigh-Ritz method. For many problems in solid mechanics, the two are equivalent.
3.3 Convergence of the Galerkin method It has been shown here that the Galerkin method works. How well it works requires a priori error estimation. This is reserved for a later section as it relies on some details
January 21, 2011
Version 0.2
33
3 The Galerkin method of the finite element method. Crucially, convergence requires that: lim uh = u
(3.16)
h →0
For finite element analysis, this corresponds to the exact solution being approached upon mesh refinement. A major question which arises is how fast the exact solution is approached as h is reduced. Details of these procedures and more elaborate mathematical analysis of the issue considered in this chapter can be found in a range of books relating to the mathematics of the finite element method (Braess, 2001; Brenner and Scott, 1994; Reddy, 1998; Strang and Fix, 1973).
3.4 Exercises 1. For elasticity, show that a solution to the Galerkin problem corresponds to a minimisation of the potential energy.
34
Version 0.2
January 21, 2011
4 Formulation of the finite element method In this chapter, a key component of the finite element method is developed. This is the discretisation of the governing equations using finite element shape functions. The unknown field uh (the displacement for elasticity problems) will be described using basis functions (the shape functions) and discrete nodal values which represent the amplitude of the basis functions. Discretisation is a step toward computer implementation. The terminology in this chapter relates primarily to elasticity problems. The developments are however general and the same procedure can be used to solve other problems.
4.1 Finite element method with piecewise linear basis functions The simplest finite element shape functions are continuous, piecewise linear functions. Consider the elastic bar in Figure 4.1, which is restrained at both ends and loaded by a distributed force f . Along the bar, a number of ‘nodes’ are located, and the domain between two nodes is known as an element. Associated with each node is a hat-like basis function which has a value of one at the node, and zero at all other
f
1 0 0 1 0 1 0 1 1
2
3
4
5
6
7
8
1 0 0 1 0 1 0 1
L Figure 4.1: One-dimensional bar.
35
4 Formulation of the finite element method nodes. The basis function at node i, is given by x x i −1 − x i −1 < x ≤ x i , x − x x i − x i −1 i −1 i x i +1 x Ni = − x i < x < x i +1 , x i − x i +1 x i − x i +1 0 otherwise,
(4.1)
where xi is the coordinate of node i. In finite element analysis, these basis functions are known as ‘shape functions’. For a one-dimensional bar, discretised with n nodal points, the approximate displacement field uh is given by uh =
n
∑ Ni ( x) ui ,
(4.2)
i =1
where ui is the value of the field uh at node i. Hence, the approximate displacement field is given in terms of the shape functions and the approximate displacement at a finite number of points (nodes). From this expression, the strain field can be computed easily by taking the derivative of the shape functions, h = ǫh = u,x
n
dNi u. dx i i =1
∑
(4.3)
To develop the Galerkin problem, an expression is needed for the weight function w h . It is expressed using the same basis functions as the displacement field, hence wh =
n
∑ Ni wi .
(4.4)
i =1
Recalling the weak form for a one-dimensional elastic rod, Z
Ω
h h Eu,x dΩ = w,x
Z
Ω
w h f dΩ +
Z
Γh
w h h dΓ
∀ wh ∈ V h ,
(4.5)
and inserting the expression for the expressions for the approximate displacement field and the weight function, ! ! Z n
n
Ω
∑ Ni,x wi
E
∑ Nj,x u j
dΩ
j =1
i =1
=
Z
n
Ω
∑ Ni wi i =1
!
f dΩ +
Z
n
Γh
∑ Ni,x wi i =1
!
h dΓ
(4.6)
for all possible values of wi . Given that wi is not a function of spatial position, it can be taken outside of the integrals, ! Z Z Z n n n n ∑ wi Ni,x E ∑ Nj,x u j dΩ = ∑ wi Ni f dΩ + ∑ wi Ni h dΓ. (4.7) i =1
36
Ω
j =1
i =1
Version 0.2
Ω
i =1
Γh
January 21, 2011
4.2 General finite element basis functions Likewise, u j can also be taken outside of the integrals, n
n
∑ wi
∑ uj
i =1
j =1
Z
Ω
Ni,x E Nj,x dΩ
!
n
=
∑ wi i =1
Z
n
Ω
Ni f dΩ + ∑ wi i =1
Z
Γh
Ni h dΓ.
(4.8)
This equation must hold for all possible combinations of wi . Consider the case that all but one wi is equal to zero. Therefore, for each i n
∑ uj j =1
Z
Ω
Ni,x E Nj,x dΩ =
Z
Ω
Ni f dΩ +
Z
Γh
Ni h dΓ
(4.9)
must hold. This yields n equations (one for each wi , and n unknowns in the form of ui ), which can be expressed as the system of equations, Kij u j = f i ,
(4.10)
where Kij =
Z
Ω
Ni,x E Nj,x dΩ
(4.11)
and fi =
Z
Ω
Nj f dΩ +
Z
Γh
Nj h dΓ
(4.12)
The matrix K = Kij is commonly known as the stiffness matrix, and f = f i as the right-hand side (RHS) vector. Solving this linear system of equations provides u j , and hence an expression for the approximate displacement field along the rod. In order the solve the linear system, Dirichlet boundary conditions must be enforced. For the zero Dirichlet conditions, this is relatively simple. If a boundary condition is applied at node k, then the row and the column k is removed from the stiffness matrix, and the kth term is deleted from the vector f . A discussion of more general boundary conditions is delayed until the end of the chapter.
4.2 General finite element basis functions The concept introduced in the previous section can be extended to higher spatial dimensions and more complex shape functions. It is useful to define shape functions on finite elements. Matrices and vectors can be built for each element before being assembled into global matrices and vectors. A mesh divides a body into a number of elements. A typical finite element mesh is shown in Figure 4.2. The problem has been divided into many quadrilateral elements. A finite element mesh consists of nodal points and elements. Nodes are points in space, while elements are lines (1D), surfaces (2D) or volumes (3D) which are constructed by joining nodal points. Elements may not overlap. A simple finite element
January 21, 2011
Version 0.2
37
4 Formulation of the finite element method
Figure 4.2: Typical two-dimensional finite element mesh.
6 7
4 3
8 4
7
5
8
3 1
6 5
2
9 1
2
Figure 4.3: Simple two-dimensional finite element mesh.
38
Version 0.2
January 21, 2011
4.2 General finite element basis functions
1
1
1
2
L x=0 Figure 4.4: One-dimensional linear bar element and associated shape functions.
mesh is shown in Figure 4.3. In Figure 4.3, the mesh is constructed with triangular elements. Both nodes and elements are numbered. Each element has three nodes. As in the one-dimensional case, the unknown (displacement) field uh is discretised by developing a method in which the displacement at any point in a body is determined in terms of a discrete number of values which are stored at the nodes (known as degrees of freedom) and basis functions. The displacement field is given by a linear combination of a finite number of basis functions (a summation of the basis functions, each multiplied by a nodal value, the amplitude). In the finite element method, the basis functions are defined on simple geometric elements, known as finite elements. Finite element shape functions are typically piecewise continuous polynomial functions with a ‘compact support’. Each shape function Ni is associated with a node i. They are characterised by being equal to unity at their node, and zero at all other nodes. Having a compact support, shape functions are non-zero only close to their node. A finite element shape function is only non-zero on elements to which it is attached. Once the displacement field uh and its derivatives ∇uh can be expressed in terms of nodal unknowns, the discretised field can be inserted into the weak governing equations. Shape functions for various elements are discussed in detail in the following sections.
4.2.1 One-dimensional bar elements To begin, linear one dimensional elements, as introduced in Section 4.1, are examined. Consider the one-dimensional bar element in Figure 4.4. The bar has two nodes, denoted by the solid circles. For a node i, the shape function associated with that node, Ni , is equal to one at the node and zero at all other nodes. The displacement field inside the bar in terms of discrete nodal values and shape functions is given by: uh ( x ) = N1 ( x ) a1 + N2 ( x ) a2 ,
January 21, 2011
(4.13)
Version 0.2
39
4 Formulation of the finite element method where ai is the displacement at node i (it is ‘stored’ at the node). The shape functions corresponding to each node are given by: N1 = − N2 =
x + 1, L
(4.14)
x . L
(4.15)
Hence, x x u h = − + 1 a1 + a2 . L L
(4.16)
This type of element is known as a linear element, as its shape functions are linear polynomials. The weak form of the governing equations involves derivatives with respect to the spatial position. Therefore, to solve a problem, the derivative of the displacement field with respect to x (the strain) is required. Taking the derivative of equation (4.16) with respect to x, ǫh =
duh dN1 dN2 1 1 = a + a2 = − a1 + a2 . dx dx 1 dx L L
(4.17)
In the case of a linear displacement interpolation, the derivative of the shape function with respect to x is a constant function within an element. Now, in terms of nodal degrees of freedom ae , the displacement and strain can be calculated in any part of the element. In matrix-vector notation, the displacement field is expressed as: uh = Nae ,
(4.18)
where
and
N = N1
N2 ,
(4.19)
a1 . ae = a2
(4.20)
The strain is given by: ǫh = Bae ,
(4.21)
where the matrix B is equal to: dN1 dN2 B= . dx dx
(4.22)
The shape functions for two linear elements are shown in Figure 4.5. Note also that these shape functions posses C0 continuity – they are continuous but their derivatives involve jumps across element boundaries. The matrix-vector notation for uh and ǫh may seem trivial for one-dimensional problems. It does however make the generalisation to multiple dimensions simple.
40
Version 0.2
January 21, 2011
4.2 General finite element basis functions
1
1
L1
1
L2
Figure 4.5: Shape functions for two one-dimensional linear bar elements. Higher-order one-dimensional elements It is possible to develop one-dimensional elements with higher-order polynomial interpolations. Higher-order elements simply have more nodes. More nodes through which the shape function must interpolate naturally means a higher-order polynomial is required. Figure 4.6 illustrates shape functions for quadratic and cubic elements. For a one-dimensional quadratic element, the displacement field is given by: uh =
3
∑ Ni (x) ai .
(4.23)
i =1
Therefore, the N matrix has the form: N = N1 N2 N3 ,
(4.24)
and the B matrix has a similar form. Since the displacement field is quadratic, obviously the strain field is linear. The interpolation basis for higher-order elements is richer since both constant, linear and quadratic (and even higher) variations in the unknown field can be described within an element. 4.2.2 Two-dimensional continuum elasticity elements In two or more dimensions, each unknown field (such as displacement in each direction) is interpolated using the polynomial shape functions. The displacement field for an element is given by: uh = Nae ,
(4.25)
and the strain field (using engineering notation) is given by: ǫ h = Bae .
(4.26)
The matrix N contains the element shape functions. In two dimensions, it has the form: N1 0 N2 0 . . . Nnn 0 , (4.27) N= 0 N1 0 N2 . . . 0 Nnn
January 21, 2011
Version 0.2
41
4 Formulation of the finite element method
1
1
1
1
2
3
quadratic shape functions
1
1
1
1
2
3
1
4
cubic shape functions Figure 4.6: Higher-order one-dimensional elements. where nn is the numbers of nodes of the element. For two-dimensional problems, the matrix B has the form: ∂N2 ∂Nnn ∂N1 0 0 . . . 0 ∂x ∂x ∂x ∂N1 ∂N2 ∂Nnn 0 ... 0 (4.28) B= 0 . ∂y ∂y ∂y ∂N1 ∂N1 ∂N2 ∂N2 ∂Nnn ∂Nnn ... ∂y ∂x ∂y ∂x ∂y ∂x
The nodal degrees of freedom (normally one for each spatial dimension at each node for elasticity problems) are stored in a vector a. For a two-dimensional problem, a1 x a1 y a 2x a2 y . (4.29) ae = .. . a nn x ann y
The simplest element in two-dimensions is the three-node triangle. It is shown in Figure 4.7. Its shape functions are of the form:
42
Version 0.2
January 21, 2011
4.2 General finite element basis functions
( x3 , y3 ) 3
1 2
( x1 , y1 )
( x2 , y2 )
Figure 4.7: Three-node triangular element.
Figure 4.8: Shape function for a three-node element. N = c1 x + c2 y + c3 .
(4.30)
The shape function for a node is illustrated in Figure 4.8. Given that a shape function should have a value of unity at its node and zero at other nodes, the coefficients c can be found by solving a linear system of equations. Consider node 1 in Figure 4.7. Its shape function must satisfy: N1 ( x1 , y1 ) = c1 x1 + c2 y2 + c3 = 1,
(4.31)
N1 ( x2 , y2 ) = c1 x2 + c2 y2 + c3 = 0, N1 ( x3 , y3 ) = c1 x3 + c2 y3 + c3 = 0,
(4.32) (4.33)
where ( xi , yi ) is the location of the ith node. This can be cast in a matrix form as: x 1 y 1 1 c 1 1 x2 y2 1 c2 = 0 . (4.34) c3 x3 y3 1 0
Solving the above system of equations yields the coefficients for the shape function of node 1. It is clear from the linear form of the shape functions for the three-node triangle that the derivatives are constant, hence the strain in the element is constant. Example 4.1 Consider the element in Figure 4.7. Given the following coordinates of each of
January 21, 2011
Version 0.2
43
4 Formulation of the finite element method the nodes, find the shape function for node 2.
( x1 , y1 ) = (1, 2) ( x2 , y2 ) = (4, 1) ( x3 , y3 ) = (2, 3) To find the coefficients, a system of equations (see equation (4.34)) is formed. Since the shape function for node 2 is being calculated, the only non-zero term on the RHS corresponds to node 2. 1 2 1 c1 0 4 1 1 c 2 = 1 c3 2 3 1 0 Solving this system of equations gives: N2 = 0.25x − 0.25y + 0.25 Once the shape functions have been computed for each node, the N and B matrices can be formed. The N matrix for this element is of the form: N1 0 N2 0 N3 0 N= 0 N1 0 N2 0 N3 and the B matrix is of the form: N1,x 0 N2,x 0 N1,y 0 N2,y B= 0 N1,y N1,x N2,y N2,x
N3,x 0 N3,y
0 N3,y N3,x
A very commonly used element in finite element analysis is the four-node quadrilateral. Its shape functions are of the form: N = c1 xy + c2 x + c3 y + c4 .
(4.35)
It is also know as a bilinear element, due to the presence of the xy terms in the shape functions. A four-node quadrilateral element is shown in Figure 4.9. A shape function for a four-node element is shown in Figure 4.10. Unlike the three-node triangle, the strain in this element is not constant. Higher-order solid elements As in one-dimensional problems, plane and three-dimensional elements can have higher-order interpolations. The shape functions for a six-node triangle are of the form: N = c1 x2 + c2 y2 + c3 xy + c4 x + c5 y + c6
44
Version 0.2
(4.36)
January 21, 2011
4.2 General finite element basis functions 4 3 1
2 Figure 4.9: Four-node quadrilateral element.
Figure 4.10: A shape function for four-node quadrilateral element. As with the three-node triangle, the coefficients can be found by solving a linear system of equations. Solving the below system of equations would yield the coefficients for node 1,
x12
y21
x1 y1
x1
y1
2 x2 2 x3 2 x4 x2 5
y22
x2 y2
x2
y2
y23 y24 y25 y26
x3 y3
x3
y3
x4 y4
x4
y4
x5 y5
x5
y5
x6 y6
x6
y6
x62
c 1 1 1 0 c2 1 0 1 c3 = c4 0 . 1 0 c 5 1 0 c 6 1
(4.37)
Two shape functions for a six-node triangle are shown in Figure 4.11 The polynomial terms needed for triangular elements can be taken from a Pascal triangle (see Figure 4.12). Note that triangular elements are complete – that is the shape functions contain all polynomial terms up to a given order. Each time the order of a triangle is increased, a new ‘row’ of terms from the Pascal triangle is added to the shape functions. Linear, quadratic and cubic triangular elements are shown in Figure 4.13. Higher-order quadrilateral elements can also be constructed. They belong to one of two families: serendipity or Lagrange. Serendipity elements have nodes only on element boundaries, whereas Lagrange elements have nodes on the element interior. Several serendipity and Lagrange finite elements are shown in Figure 4.14. A Pascal triangle for serendipity elements is shown in Figure 4.15. Note that serendipity
January 21, 2011
Version 0.2
45
4 Formulation of the finite element method
Figure 4.11: Two shape functions for a six-node triangular element.
1 x
y
x2 x3 x4 ...
x3 y ...
y2
xy x2 y
xy2 x 2 y2
...
...
y3 xy3
...
y4 ...
...
Figure 4.12: Pascal triangle – polynomial terms for triangular elements of order n.
Figure 4.13: Linear, quadratic and cubic triangular elements.
46
Version 0.2
January 21, 2011
4.2 General finite element basis functions
(a)
(b) Figure 4.14: Serendipity (a) and Lagrange (b) quadrilateral finite elements. 1 x x2 x3 x4 ...
y y2
xy x2 y
xy2
x3 y
y3 xy3
...
y4 ...
...
Figure 4.15: Pascal triangle for serendipity elements.
elements of order higher than three ’miss’ polynomial terms, and are therefore not recommended. Internal nodes must be introduced to ensure all polynomial terms are included. The Pascal triangle for Lagrange elements is shown in Figure 4.16. The elements do not suffer from the limitation of missing polynomial terms for any order interpolation. Eight- and nine-node quadrilateral elements are commonly used in finite element analysis. The shape functions for the eight-node serendipity quadrilateral have the form: N = c1 x2 y + c2 xy2 + c3 x2 + c4 y2 + c5 xy + c6 x + c7 y + c8 ,
(4.38)
and for the nine-node Lagrange quadrilateral the shape functions have one extra term and are of the form: N = c1 x2 y2 + c2 x2 y + c3 xy2 + c4 x2 + c5 y2 + c6 xy + c7 x + c8 y + c9 .
January 21, 2011
Version 0.2
(4.39)
47
4 Formulation of the finite element method 1 x
y
x2 x3 ...
x2 y x3 y
...
y2
xy xy2 x 2 y2 x 3 y2
...
y3 ...
xy3 ...
x 2 y3 x 3 y3
...
Figure 4.16: Pascal triangle for Lagrange elements.
Figure 4.17: Tetrahedral and brick three-dimensional finite elements. 4.2.3 Three-dimensional elements Three-dimensional elements are increasingly used in engineering analysis. They can be formed in the same fashion as one- and two-dimensional elements. Commonly used shapes are tetrahedra and bricks. These two elements are shown in Figure 4.17. The computational effort required for three-dimensional analysis can become high. Not only do the elements tend to have more nodes, they often have more degrees of freedom per node. A four-node tetrahedral element has linear shape functions. They are of the form: N = c1 x + c2 y + c3 z + c4 .
(4.40)
Similar to the three-node triangle, the derivatives of the shape functions are constant within an element. The simplest brick element has eight nodes. Its shape functions are of the form: N = c1 xyz + c2 xy + c3 xz + c4 yz + c5 x + c6 y + c7 z + c8 .
(4.41)
This is known as a ‘trilinear element’. Higher-order version of both tetrahedral and brick elements also exist. The numbers of nodes per element increases rapidly. Also, wedge type elements are often used for generating meshes for complex geometries.
48
Version 0.2
January 21, 2011
4.3 Governing equations x = l1
x=0
x = l2
x=L
x
Figure 4.18: One-dimensional bar showing the element that extends from x = l1 to x = l2 . The three-dimensional equivalent of a Pascal triangle can be used to find the polynomial terms for higher-order three-dimensional elements (Zienkiewicz and Taylor, 1989)
4.3 Governing equations At this point, it is known how the displacement and the displacement gradient at a point in space can be expressed in terms of the nodal variables (degrees of freedom). What remains is the insertion of the discretised field uh into the governing weak equation. The weight (test) functions are discretised in the same fashion as the unknown field uh . Therefore, for the multi-dimensional case within an element, wh = Nbe , s
(4.42)
h
∇ w = Bbe ,
(4.43)
where the notation has been ‘abused’ as ∇s wh is being expressed as a vector (recall that this is possible due to symmetry). Now, the one-dimensional governing equation for an elastic rod (equation (3.1)) is considered. Take a single element (of any order) which extends from l1 to l2 (see Figure 4.18). The discretised governing equation for the element is expressed as: Z l2 l1
( Bbe ) T EBae dΩ =
Z
Γh,e
( Nbe ) T h dΓ,
(4.44)
where the integral over Γh,e is only non-zero if the boundary of the element coincides with the boundary Γh . The discrete nodal values can be removed from the integrals, which after some rearranging (using ( Ac) T = c T A T ) leads to: be T
Z l2 l1
B T EB dΩ ae = be T
Z
Γh,e
N T h dΓ.
(4.45)
The be terms appear on both sides of the equality, and can therefore be eliminated, Z l2 l1
B T EB dΩ ae =
January 21, 2011
Z
Γh,e
N T h dΓ.
Version 0.2
(4.46)
49
4 Formulation of the finite element method The term ke =
Z l2 l1
B T EB dΩ
(4.47)
is known as the ‘element stiffness matrix’. Once the stiffness matrix has been formed for an element, its contribution is added to the ‘global’ stiffness matrix K. For a continuum elasticity problem, the discretised displacement and strain fields are inserted into equation (3.5). This yields: Z
Ωe
( Bbe ) T DBae dΩ =
Z
Γh,e
( Nbe ) T h dΓ +
Z
Ωe
( Nbe ) T b dΩ
(4.48)
(recall that b denotes the body force). Following the same steps as for the onedimensional problem, the element stiffness matrix is expressed as: Z
Ωe
B T DB dΩ ae =
Z
Γh,e
N T h dΓ +
Z
Ωe
N T b dΩ.
(4.49)
Once the stiffness matrix has been formed for each element in the problem, it must be assembled into the global stiffness matrix. This is denoted symbolically by the operation: K = Ane e =1 k e , ne f = A e =1 f e ,
(4.50) (4.51)
where A represents the assembly operation and ne is the number of elements in the mesh. For an element stiffness matrix, the term k ij relates local nodes i and j. In the global mesh, local node i has global node number I, and similarly for J. The component k ij is added to the location K I J . The assembly process is discussed in more detail in Chapter 5.
4.4 Isoparametric mapping In practice, isoparametric elements are used in finite element software. The shape functions are formed for a simple element configuration (typically elements with unit length side and with sides aligned with the coordinate system and with a convenient origin). This approach has a number of advantages. From the point of view of implementation, it requires the programming of only one function to evaluate the shape functions of a type of element, irrespective of the exact shape of the element. It also allows the simple application of numerical integration, as is discussed in the following section. A third advantage is that isoparametric mapping allows higher-order elements to have curved edges. An isoparametric mapping is defined via a mapping using the nodal shape functions. For a two dimensional problem, the point ( x, y) in the ‘physical’ Cartesian
50
Version 0.2
January 21, 2011
4.4 Isoparametric mapping x 3 (−1,1)
(1,1)
4
2
3
4
η ξ 1 (−1,−1)
y 2
x
(1,−1)
1
ξ
Figure 4.19: Isoparametric mapping for a bi-unit square.
configuration is given by: nn
∑ Ni (ξ, η ) xi ,
x=
i =1 nn
(4.52)
∑ Ni (ξ, η ) yi ,
y=
i =1
where (ξ, η ) are known as the ‘natural coordinates’, and nn is number of nodes of the element. The mapping for a four-node quadrilateral element is illustrated in Figure 4.19. The map x takes a point within the bi-unit square in the (ξ, η ) domain and maps it to a point in the ‘real’ element in the Cartesian ( x, y) system. The position on the bi-unit square is given by the natural coordinates, ξ and η. The map ξ is the inverse of x, taking a point in the real Cartesian ( x, y) system and mapping to a point (ξ, η ) in the bi-unit square. Using an isoparametric mapping, shape functions can be defined on simple shapes, such as the bi-unit square. For example the displacement in the x-direction at a point is given by:
uh =
nn
∑ Ni (ξ, η ) aix .
(4.53)
i =1
The difficulty that arises is when computing derivatives with respect to the real coordinates, as the shape functions are defined in terms of the natural coordinates. In the governing weak equations, derivatives of the unknown field with respect to x and y are required. To proceed, consider the application of the product rule for differentia-
January 21, 2011
Version 0.2
51
4 Formulation of the finite element method tion for a function f , ∂f ∂ f ∂x ∂ f ∂y = + , ∂ξ ∂x ∂ξ ∂y ∂ξ ∂f ∂ f ∂x ∂ f ∂y = + . ∂η ∂x ∂η ∂y ∂η
(4.54) (4.55)
This can be written in matrix form as: ∂x ∂y ∂ f ∂f ∂x ∂ξ ∂ξ ∂ξ = , ∂f ∂x ∂y ∂f ∂η ∂η ∂η ∂y | {z }
(4.56)
J
where the matrix J is commonly known as the Jacobian matrix. Taking the inverse of the Jacobian, ∂y ∂ f ∂y ∂f − ∂x ∂ξ 1 ∂η ∂ξ , = (4.57) ∂f ∂f ∂x ∂x j − ∂y ∂η ∂ξ ∂η | {z } J
where
j = det J.
(4.58)
In light of the isoparametric mapping (equation (4.52)), the terms in the matrix J can be computed. Using the shape function, defined in terms of ξ and η, nn ∂x ∂Ni =∑ x, ∂ξ ∂ξ i i =1 nn ∂x ∂Ni =∑ x, ∂η ∂η i i =1
(4.59)
nn ∂Ni ∂y =∑ y, ∂ξ ∂ξ i i =1 nn ∂y ∂Ni =∑ y. ∂η ∂η i i =1
Once terms of the matrix J have been formed, from equation (4.57) the derivatives of the shape function of node i with respect to x and y can be computed, ∂Ni ∂Ni nn nn − ∑ j=1 Nj,ξ y j ∂ξ ∂x 1 ∑ j=1 Nj,η y j = . (4.60) nn ∂N ∂Ni j − ∑nn Nj,η x j N x ∑ j,ξ j i j =1 j =1 ∂y ∂η 52
Version 0.2
January 21, 2011
4.4 Isoparametric mapping s
3
s (0,1)
3
(0,1)
5
6
2 1
(0,0)
r
2 1
(1,0)
(0,0)
4
r
(1,0)
Figure 4.20: Three-node and six-node triangles using area coordinates.
Now, once the derivatives of the shape functions have been calculated on the simple isoparametric domain, they can be calculated for the real configuration. While this procedure may seen complex, it can be implemented in a computer code in a simple and compact fashion.
Quadrilateral elements The shape functions for quadrilateral elements on a bi-unit square is particularly simple. For the bi-unit quadrilateral in Figure 4.19, the shape function of the ith node is given by: Ni (ξ, η ) =
1 (1 + ξ i ξ ) (1 + ηi η ) . 4
(4.61)
Similar formulae can be found for eight- and nine-node elements (Hughes, 1987).
Triangular elements Isoparametric mappings are also applied for triangular elements. A popular approach is to ‘degenerate’ a quadrilateral element to a triangular element (Bathe, 1996; Hughes, 1987). This way, a four-node quadrilateral element can be reduced to a threenode triangular element, and an eight-node quadrilateral element can be reduced to a six-node triangular element. Is is also possible to directly address isoparametric triangular elements. The latter approach is followed here. Triangular elements are often described using area coordinates. Consider the triangle in Figure 4.20. The position of a point inside the triangle is given by r, s and t, where: t = 1 − r − s.
January 21, 2011
(4.62)
Version 0.2
53
4 Formulation of the finite element method The shape functions for the nodes of the element in Figure 4.20 are then: N1 = t, N2 = r, N3 = s.
(4.63) (4.64) (4.65)
Note the numbering of the nodes for the quadratic triangle in Figure 4.20. It is common for higher-order elements that the apexes are first numbered (corner nodes in the case of quadrilateral elements), then the mid-side nodes. For a quadratic (six-node) triangular element, the shape functions are equal to: N1 = 2t2 − t,
N4 = 4rt,
2
N2 = 2r − r,
N5 = 4rs,
N3 = 2s2 − s,
(4.66)
N6 = 4st.
There is a formula for the shape functions of arbitrary order in terms of r, s and t on triangles (Hughes, 1987, p. 166).
4.5 Numerical integration At this stage, the problem is not yet fully discrete. While the displacement field has been discretised, the element stiffness matrices and RHS vectors still involve integration over volumes and surfaces (see section 4.3). To make the formulation fully discrete (and amendable to computer implementation), numerical integration is applied. Before proceeding to the element matrices, numerical integration in one-dimension is considered. For the integration of a function f (ξ ) from −1 to 1, Z 1
−1
nint
f (ξ ) dξ ≃
∑ f ( ξ i ) wi ,
(4.67)
i =1
where ξ i are discrete points on the integral domain, nint is the number of discrete points (integration points) at which the function is evaluated and wi is the weight asn signed to each point ∑i=int1 wi = 2 . Numerical integration is illustrated in Figure 4.21. Different integration schemes specify where the points are located and the weight associated with each point. Two schemes which arise in finite element analysis are Newton-Cotes and Gauss integration. To evaluate the stiffness matrix of an element, it convenient to integrate in the (ξ, η ) domain. To do this, consider in two dimensions that: Z
54
Ωe
B T DB dΩ =
Z 1 Z 1
−1 −1
B T DBj dξ dη.
Version 0.2
(4.68)
January 21, 2011
4.5 Numerical integration f
ξ1
−1
ξ2
ξ3
ξ4
1
ξ
Figure 4.21: Numerical integration of the function f . Note the appearance of the determinant of the Jacobian, j. This is due to dx dy = j dξ dη (the substitution rule for integration has been applied). Application of numerical integration leads to: Z 1 Z 1
−1 −1
B T DBj dξ dη ≃
n
∑ B (ξ i , ηi )T DB (ξ i , ηi ) j (ξ i , ηi )
wi .
(4.69)
i
Note that B (ξ i , ηi ) is the usual B matrix, containing derivatives of the shape functions with respect to x and y, evaluated at the point (ξ i , ηi ). To perform the numerical integration of the element stiffness matrix, the last remaining issue is the selection of an appropriate integration scheme. A Newton-Cotes scheme uses equally spaced integration points. The sampling points and weights are listed in Table 4.1. Importantly, (n + 1) integration points are required when using a Newton-Cotes scheme to integrate an nth-order polynomial. Most commonly in finite element analysis, Gauss integration is used. It is also known as ‘Gauss quadrature’. Gauss quadrature is the optimal numerical integration scheme for polynomials in one dimension. Optimal means that to integrate a polynomial exactly, it requires the least number of points. For a (2n − 1)th-order polynomial, Gauss integration requires n points to integrate the polynomial exactly. The location and weights for Gauss integration in one dimension for up to three points are given in Table 4.2. Three points integrate as 5th order polynomial exactly in one dimension. For multi-dimensional problems, the one-dimensional scheme can be extended in each spatial direction. This is simple when using isoparametric elements. However, in multiple dimensions Gauss integration is not necessarily the most efficient scheme. The question that arises in finite element analysis is: How many integration points should be used? Generally, full integration schemes are recommended for their robustness. Full integration means that the stiffness matrix is integrated exactly (assuming j = constant). The integrand of the stiffness matrix is B T DB, where B contains derivatives of the shape functions. Therefore, if the shape functions are linear, the terms in B are constant and the stiffness matrix is also constant throughout the element. Therefore, one integration point is sufficient (in one, two and three dimensions). If the shape functions are quadratic, B contains linear terms, which when squared give a
January 21, 2011
Version 0.2
55
4 Formulation of the finite element method
n 1
location ξ i 0
weight wi 2
2
1
1
1
1
3
1 3 4 3 1 3 1 4 3 4 3 4 1 4
−1 0 1
4
−1 −
1 3 1 3 1
Table 4.1: Newton-Cotes integration rules for one dimension on the domain [−1, 1].
n 1 2
3
location ξ i 0
weight wi 2
1 −√ 3 1 √ 3 r 3 − 5 0 r
3 5
1 1 5 9 8 9 5 9
Table 4.2: Gauss integration rules for one dimension on the domain [−1, 1].
56
Version 0.2
January 21, 2011
4.6 Imposition of Dirichlet boundary conditions
Figure 4.22: Spurious zero energy mode for the eight-node quadrilateral element when using 2 × 2 integration. quadratic variation. To integrate a quadratic function in two dimensions, two points are required in each direction. For some applications it may be advantageous to use reduced integration schemes in combination with particular elements. The four-node quadrilateral element requires 2 × 2 points for full integration, although it is often integrated with just one point. Similarly, the eight-node quadrilateral requires 3 × 3 points for full integration, although a 2 × 2 scheme is often used. Reduced integration schemes can improve the performance of some elements for special applications. The danger of underintegrating elements is that spurious modes may arise. This means that a deformation mode exists which does not contribute to the energy. This deformation mode can develop in an uncontrolled fashion. Figure 4.22 shows the classic ‘hour-glass’ spurious mode for an eight-node quadrilateral when using reduced 2 × 2 integration. It is often argued that a zero energy mode will be restrained by neighbouring elements, although it is safer to use full integration. A test for the element stiffness matrix is to calculate its eigenvalues. The number of zero eigenvalues indicates the number of rigid body modes - modes which do not contribute to the energy. In one dimension, there should only be one zero eigenvalue which corresponds to translation. In two dimensions, there should be only three zero eigenvalues, two relating the translation (x and y directions) and one for rigid body rotation. If a two-dimensional element has more than three zero eigenvalues, there is likely a deficiency in the formulation. In summary, recommended integration schemes for commonly used one- and twodimensional elements are shown in Table 4.3.
4.6 Imposition of Dirichlet boundary conditions The imposition of Neumann (force) boundary conditions in the finite element method is straightforward. They appear naturally in the weak formulation (hence the common name ‘natural boundary conditions’). Non-zero Dirichlet boundary conditions are more complicated to enforce. In the Galerkin method, the trial functions (uh ) must satisfy the Dirichlet boundary conditions. In solid and structural mechanics,
January 21, 2011
Version 0.2
57
4 Formulation of the finite element method
element
integration scheme
two-node bar
1
three-node bar
2
three-node triangle (T3)
1
six-node triangle (T6)
or
3
four-node quadrilateral (Q4)
2×2
eight-node quadrilateral(Q8)
3×3
Table 4.3: Recommended integration schemes for commonly used elements.
58
Version 0.2
January 21, 2011
4.6 Imposition of Dirichlet boundary conditions this means that the uh must satisfy the prescribed displacements where they are specified. Classical enforcement of Dirichlet boundary conditions in the finite element method relies on the property of the shape functions that they are non-zero at their own node and equal to zero at all other nodes. Two common approached are outlined below. 4.6.1 Elimination of ‘Dirichlet’ degrees of freedom The displacement field can be expressed by: uh = vh + gh ,
(4.70)
where uh is the displacement, vh is equal to zero where Dirichlet (displacement) boundary conditions are applied and gh satisfies the Dirichlet boundary conditions. In discretised form at element level, e kvv kevg fv av e , (4.71) = k a= e fg g k gv kegg where av are the nodal unknowns at nodes where no Dirichlet boundary conditions are applied, and g are the applied Dirichlet boundary conditions. Since the weight function is defined to be zero where Dirichlet boundary conditions are applied,
kevv kegv
kevg kegg
av g
=
fv . 0
(4.72)
Since g is known, it is not included in the global system of equations to be solved. It is possible to write: kevv av + kevg g = f v ,
(4.73)
where av is the only unknown. Since g is known, Dirichlet boundary conditions can be enforced by modifying the RHS vector, kevv av = f e − kevg g.
(4.74)
In terms of indexes, the element RHS vector is modified, f ie,mod = f ie −
nn
∑ keij gej .
(4.75)
j =1
Once the element stiffness matrix has been formed, the element RHS vector can be modified to enforce Dirichlet boundary conditions before assembly into the global RHS vector. The advantages of this approach are that degrees of freedom to which Dirichlet boundary conditions are applied do not appear in the global stiffness matrix and symmetries which may be present in the stiffness matrix are preserved.
January 21, 2011
Version 0.2
59
4 Formulation of the finite element method 4.6.2 Retention of ‘Dirichlet’ degrees of freedom A simpler approach to imposing Dirichlet boundary conditions is to set all terms on the row of the stiffness matrix corresponding to a node where a Dirichlet boundary condition is applied (row k) equal to zero, except the diagonal term which is set equal to one. The value of Dirichlet boundary condition is then inserted into the vector f at position k. The advantage of this scheme is its simplicity, and that the value of the Dirichlet boundary condition is returned in the solution vector which simplifies postprocessing. A disadvantage is that the size of the system of equations to be solved is larger than when the extra degrees of freedom are eliminated (but only marginally so in a typical simulation) and if the stiffness matrix is symmetric before the imposition of Dirichlet boundary conditions, symmetry is lost. However, modern iterative linear solvers can deal with the consequences of boundary conditions being applied in this fashion is a very efficient manner (Ern and Guermond, 2004).
4.7 Exercises 1. For a plane elasticity problem, at which points is the stress from the finite element solution not uniquely defined when using C0 elements? 2. Form the shape functions for a Lagrange nine-node quadrilateral element 3. The three-node triangle element is often referred to as the ‘constant strain triangle’ (CST). Show why this is. 4. Under what conditions the matrix B T DB is symmetric? 5. Consider the Laplace equation, which is given by: ∆u = ∇ · (∇u) = 0 To solve this problem, both the finite element method and the finite difference method require the solution of a system of equations Ka = f . In one dimension and for equally spaced nodal points (nodes are a distance h apart), compare the matrix K for the finite element method with linear elements and for the secondorder central finite difference method. The second-order finite difference equation for the second derivative is: d2 u u − 2ui + ui+1 = i −1 dx2 i h2 For unequally spaced nodes, the difference equation for the second derivative is: 2ui+1 d2 u 2ui−1 2ui − + = 2 h j h j +1 dx i h j h j + h j +1 h j +1 h j + h j +1
Comment on the symmetry of the operators for finite elements and finite differences on unequally spaced grid.
60
Version 0.2
January 21, 2011
4.7 Exercises 6. Calculate and plot the shape functions for the eight-node quadrilateral element on a bi-unit square. 7. Formulate the stiffness matrix for the element in Figure 4.7 using the nodal coordinates in Example 4.1 and for E = 1 and ν = 0.2. Calculate the eigenvalues of the stiffness matrix and comment on their significance. 8. Form the stiffness matrix for a four-node quadrilateral on a bi-unit square. Calculate its eigenvalues using one-point and 2 × 2 integration.
January 21, 2011
Version 0.2
61
5 Implementation of the finite element method This chapter addresses practical aspects of implementing the finite element method. To this point, the underlying theory has been addressed and the discretised formulation for individual elements discussed. One of the strengths of the finite element method is the efficiency with which it can be implemented in a computer code. This chapter describes how the finite element method can be implemented in a computer code. It is illustrated with schematic extracts of computer code. The notation used in this chapter follows as closely as possible the notation used in the accompanying Matlab code. Note: Code examples in the chapter does not reflect the structure of the current version of the Matlab code used in CT5123.
5.1 Preprocessing The first step in finite element analysis is known as preprocessing. This step involves generating a finite element mesh. For simple problems with only a few elements, this can be done by hand. For problems typically analysed using the finite element method, a mesh is too complicated to produce by hand. Mesh generation programs exist which produce meshes suitable for finite element analysis. In large finite element software packages, often a mesh generator is included. The generated mesh files serve as input for a finite element program. A finite element mesh consists of nodal coordinates and a connectivity. Figure 5.1 shows a simple finite element mesh of four-noded quadrilaterals. Both nodes and elements have been numbered. Also shown are the applied boundary conditions (both prescribed displacements and the applied force). A typical input file containing
3
F= 5.3
9
4 3 2
8 5 4
2
7
1 6
1
Figure 5.1: Finite element mesh with elements and nodes numbered.
63
5 Implementation of the finite element method % node 1 2 3 4 5 6 7 8 9
x 0.0 0.0 0.0 0.5 0.6 0.76 2.1 2.1 2.1
y 0.0 0.8 2.0 2.0 1.0 0.7 0.82 1.5 2.0
Figure 5.2: Nodal coordinates. % element 1 2 3 4
node1 2 3 9 5
node 2 1 2 4 6
node 3 6 5 5 7
node 4 5 4 8 8
Figure 5.3: Element connectivity. the nodal coordinates is shown in Figure 5.2. For each node, the x and y coordinate is given. The connectivity list for Figure 5.1 is given in Figure 5.3. For each element, the node numbers are listed. The list can begin with any node of the element, although the numbering must process in an anti-clockwise direction (to avoid calculating a negative volume and hence a negative element stiffness). In addition to nodal coordinates and the connectivity, boundary conditions must be specified. The problem in Figure 5.1 involves both Dirichlet (displacement) and Neumann (force) boundary conditions. The list in Figure 5.4 gives the applied boundary condition. First, the node to which a boundary condition is applied is listed. The second column gives the boundary condition type: Dirichlet (displacement) = 1 and Neumann (force) = 0. The third column gives the degree-of-freedom to which the boundary condition is applied (dof=1 represents the x-direction and dof=2 represents the y-direction). The final column is the value of the applied boundary condition. This is either a prescribed force or the prescribed displacement. The preprocessing is completed by specify the material data. For an isotropic linear-elastic analysis, the is Young’s modulus, the Poisson’s ratio and the density of the material. The density is necessary when taking the self-weight into account.
5.2 Program flow Once the input for a finite element analysis has been prepared, the calculation phase can begin. The outline of a typical finite element program for linear static problems
64
Version 0.2
January 21, 2011
5.3 Local-global % node 1 1 2 3 9
bc type 1 1 1 1 0
dof 1 2 1 1 1
value 0 0 0 0.0 5.3
Figure 5.4: Specification of boundary conditions. is shown in Figure 5.5. It begins with reading the input data from the preprocessing stage. Once the input data has been read, each degree of freedom at which no Dirichlet boundary condition is applied is given an equation number. This is stored in an array ID(node#,dof#) (see Figure 5.6). For degrees of freedom where a Dirichlet boundarycondition is applied, ID=0. Degrees of freedom where Dirichlet boundary conditions are applied do not appear in the stiffness matrix. Hence, for the problem in Figure 5.1 the size of stiffness matrix is 14 × 14. After the number of equations has been determined, it is then possible to allocate the memory required for the calculation. Once the memory has been allocated, it is possible to start a major element of the analysis – the calculation of the element stiffness matrix and element RHS. For each element, the element stiffness matrix ke is formed and added into the global stiffness matrix K. Similarly, the element RHS vector f e for each element is formed and added into the global RHS vector f . If a Dirichlet boundary condition is applied to a degree of freedom (or degrees of freedom) of the element, the RHS vector is modified to impose the boundary condition. The element RHS vector is also assembled into the global RHS vector. After looping over all the elements and forming the stiffness matrix and RHS vector, Neumann boundary conditions are added into the RHS vector. The second major step is the solution of the system of equations: Ka = f . Various methods can be used to do this. The best method depends on the nature of the problem and its size. This aspect is discussed in Chapter 8. After the system of equations has been solved, the displacement can be calculated at any point in the body. For large problems, the amount of data can be immense. This data requires post-processing to interpret the data and to extract the desired information.
5.3 Local-global A key to implementing the finite element method is the local-global relationship between nodes. All nodes in a mesh have a global number. Locally, for each element, each node of the element has a local number. The global node numbers can be seen in Figure 5.1. From the connectivity (Figure 5.3), the relationship between the local and global node numbers can be inferred for each element. Consider element 3
January 21, 2011
Version 0.2
65
5 Implementation of the finite element method
% start program % read input data read_nodes read_connectivity read_boundary_conditions read_material_data % number equations number_equations % setup memory allocation memory_allocate % loop over elements for i=1:number_elements local_arrays % copy global to local arrays element_form % form stiffness matrix k_e and f_e for element assemble_K % add element stiffness matrix into K apply_dirichlet_bc % modify f_e for Dirichlet boundary conditions assemble_f end
% add element RHS into global RHS
% apply Neumann boundary conditions apply_neumann_bc % solve system Ku = f u = K\f % post-process results calculate_stress plot_displacements % end program
Figure 5.5: Finite element program flow for a linear problem.
66
Version 0.2
January 21, 2011
5.3 Local-global ID(1,1) ID(2,1) ID(3,1) ID(4,1) ID(5,1) ID(6,1) ID(7,1) ID(8,1) ID(9,1)
= = = = = = = = =
0 0 0 3 5 7 9 11 13
ID(1,2) ID(2,2) ID(3,2) ID(4,2) ID(5,2) ID(6,2) ID(7,2) ID(8,2) ID(9,2)
= = = = = = = = =
0 1 2 4 6 8 10 12 14
Figure 5.6: Equation numbering. 1(9)
2(4)
4(8) 3 (5)
Figure 5.7: Global/local numbering of element 3. The global node number is shown in brackets. in Figure 5.1. For element three, the local and global node numbers are shown in Figure 5.7. An array (matrix) connect(element#,local node#) contains the connectivity data which is in Figure 5.3. For example, connect(3,1) = 9, connect(3,2) = 4, connect(3,3) = 5, connect(3,4) = 8. Before forming the element i, information specific to that element is copied from the global arrays containing nodal positions, displacements at nodes and any other relevant information to smaller local arrays for the element. An example piece of code is shown in Figure 5.8 for copy elements of the global arrays to the local element arrays. Once the local arrays have been assembled for element i, its stiffness matrix can be formed. A simple function to form the stiffness matrix and RHS vector of an element is shown Figure 5.9. The process is repeated for every element in the mesh. Once the stiffness matrix has been calculated for an element, it is possible to impose Dirichlet boundary conditions by modifying the element RHS vector f e . The process is shown in Figure 5.10. After the stiffness matrix and RHS vector has been formed for an element and the RHS vector has been modified for any Dirichlet boundary conditions, the element stiffness matrix and RHS vector is added into the global stiffness matrix and global RHS vector. This process relies on indexing between local degrees of freedom and their global equation number. An example piece of code is given in Figure 5.11. Note that degrees of freedom where Dirichlet boundary conditions are applied are not assembled into the global arrays.
January 21, 2011
Version 0.2
67
5 Implementation of the finite element method
% setup local arrays for element i for j=1:nodes_per_element % loop over each node x_e(j,:) = x(connect(i,j),:) % copy nodal postions % copy displacements from global to local array for k=1:dof_per_node % loop over each degree of freedom jj = j*dof_per_node - dof_per_node + k if ID(connect(i,j),k) ~= 0 a_e(jj) = a( ID(connect(i,j), k) ) else a_e(jj) = g_e(jj) % displacement is prescribed end end end % return to main program
Figure 5.8: Set-up of local element arrays for element i.
% form element stiffness matrix ip_scheme
% determine integration scheme
for i=1:number_ip_points % loop over integration points ip_location % calculate position of integration point shape_function % evaluate shape functions at integration point form_B_matrix form_D_matrix k_e = k_e + B’*D*B*w(i)*xj(i) % add contribution of ip to k_e f_e = f_e + N’*b*w(i)*xj(i) % add contribution of ip to f_e if postprocess == true % if stress is required stress(i) = D*B*a_e end end % return to main program
Figure 5.9: Formation of element stiffness matrix.
68
Version 0.2
January 21, 2011
5.3 Local-global
% impose Dirichlet boundary conditions for i=1:nodes_per_element*dof_per_node f_e(:) = f_e(:) - k_e(:,i)*g_e(i) end % return to main program
Figure 5.10: Imposing Dirichlet boundary conditions.
% assemble k_e and f_e for element i into K and f for j=1:nodes_per_element for k=1:dof_per_node ii = ID(connect(i,j), k) kk = j*dof_per_node - dof_per_node + k for l=1:nodes_per_element for m=1:dof_per_node jj = ID(connect(i,l), m) ll = l*dof_per_node - dof_per_node + m if ii~=0 & jj~=0 K(ii,jj) = K(ii,jj) + k_e(kk,ll) end end end if ii ~= 0 f(ii) = f(ii) + f_e(kk) end end end % return to main program
Figure 5.11: Assembly of local arrays into global arrays.
January 21, 2011
Version 0.2
69
5 Implementation of the finite element method 0
500
1000
1500 0
500
1000
1500
Figure 5.12: Graphical representation of a stiffness matrix
5.4 Stiffness matrix storage A feature of the finite element method is that the stiffness matrix is sparse. This means that the stiffness matrix contains many zeros. Taking advantage of the sparse nature of the stiffness matrix yields significant computational savings, both in terms of memory and computational time. Figure 5.12 represents graphically the structure of the stiffness matrix for a problem involving 700 four-noded quadrilateral elements. This is relatively small for a finite element analysis. Non-zero terms are marked by a dot. The dimension of the matrix is 1530 × 1530 (giving a total of 2 340 900 entries). However, only 23 638 entries are non-zero (1% of the total matrix elements). Assuming double precision storage, storing the whole matrix would require almost 18 megabytes of memory. Storing only the non zero terms requires on 283 kilobytes of memory! In addition, solvers are implemented such that operations are not performed on non-zero terms which reduces the number of floating point operations required, and hence the computational time, significantly. Yet further savings can be achieved by taking into account symmetry of the stiffness matrix. Many finite element programs store the stiffness matrix in a skyline format. For this method to be efficient, the mesh must be constructed in an ‘optimal’ fashion to minimise the bandwidth of the stiffness matrix. This essentially requires that the node numbers of nodes which are close to each should also be close. Typically, mesh generation programs will label nodes optimally, or have the option to relabel nodes to minimise the stiffness matrix bandwidth. The stiffness matrix in Figure 5.12 is well-numbered, as all non-zero terms are close to the diagonal which minimises the bandwidth.
70
Version 0.2
January 21, 2011
5.5 Post-processing
5.5 Post-processing After a calculation has been completed, post-processing of the data is required. Three quantities of potential interest from a linear-elastic calculation are displacements, stresses and reaction forces. It is important to realise that the stress computed at a point is not always reliable. More reliable are the reaction forces at nodes where Dirichlet boundary conditions are applied. Reaction forces can be evaluates by calculating the so-called internal force vector f int . The internal force vector for an element is given by: f eint =
Z
Ωe
B T σ dΩ +
Z
Ωe
N T b dΩ
(5.1)
Once assembled for all elements (at all nodes, including where Dirichlet boundary conditions are applied), it gives the reaction force at each node. For static problems, the internal force vector should be zero where no Dirichlet boundary conditions are applied. This means that the body is in equilibrium.
January 21, 2011
Version 0.2
71
6 Structural elements for finite element analysis Special finite elements are often used in structural mechanics. Examples are beam, plate and shell elements. They are used in situations where conventional solid elements perform poorly with a practicably allowable number of elements, or in situations where they provide a higher degree of accuracy (in terms of the quantity of interest) than plane elements for the same computational effort. These elements are derived from different governing equations, which are closely related to classical equations from structural mechanics. Their governing equations derive from the governing equation of elastostatics, with simplifications and assumptions which are applicable for common structural problems. The formulation of structural elements tends to be more complex than continuum elements. A particular problem is that classical thin bending problems in structural mechanics are governed by fourth-order equations. That is, the strong governing equations involve fourth-order derivatives, in contrast to continuum elasticity problems which involve only second-order derivatives. This has serious consequences for the shape functions which can be employed.
6.1 Rod elements in space In Chapter 4, a one-dimensional rod element was developed that could only transmit normal forces. It is however possible to assemble this element into a truss structure in two or three dimensions. A simple two-dimensional truss structure is shown in Figure 6.1. Joints (nodes) of the truss can translate in both the x- and y-directions. Hence, a finite element model will require two degrees of freedom at each node. The use of truss elements in a two- or three dimensional structure relies upon a rotation between different coordinate systems. Each element is rotated to a convenient coordinate system. The most convenient system is one in which the axis of the truss element is aligned with the x coordinate axis. Consider the truss element in Figure 6.2, which is oriented in the three-dimensional space. For a truss structure in
F Figure 6.1: Two-dimensional truss structure.
73
6 Structural elements for finite element analysis Fx* x*
y*
y z*
x z
Figure 6.2: Linear truss element in three-dimensional space. three dimensions, a node will have three degrees of freedom, and in two dimensions, two degrees of freedom. For a two dimensional truss structure, the rotation is particularly simple. Denoting the angle between the x-axis and the x ⋆ -axis as θ, the rotation matrix has the form: cos θ sin θ (6.1) QT = − sin θ cos θ The vector ae holds the displacement components at the nodes of an element. Therefore, the vector ae has four components (a component in the x and y direction at each node). The displacement components can be rotated using the matrix Q (see section 1.1) to give the nodal displacements in terms of the coordinate system aligned with the axis of the truss element. Since truss elements are based on the assumption that only normal forces (σx⋆ ) can be transmitted, displacements in the y⋆ -direction are not resisted by an element. The only displacement of any consequence for an element is in the x ⋆ -direction. The displacements at nodes one and two, of a linear bar, in the x ⋆ - and y⋆ -directions are given by: a1x⋆ Q Q21 0 0 a1x 11 ⋆ a1y a1y Q Q 0 0 22 12 = , (6.2) 0 Q11 Q21 a ⋆ 0 a 2x⋆ 2x 0 0 Q12 Q22 a2y a2y where Qij are components of the matrix Q. Given that the only displacements of interest are the in x ⋆ -direction, the above problem can be reduced to give: a1x a1y a1x⋆ Q11 Q21 0 0 . (6.3) = a2x⋆ 0 0 Q11 Q21 a 2x a2y
The following definitions are now adopted: a1x⋆ ⋆ , ae = a2x⋆ 74
Version 0.2
(6.4)
January 21, 2011
6.1 Rod elements in space and
Q11 q= 0
Q21 0
0 Q11
0 . Q21
(6.5)
The displacement in the x ⋆ -direction along the bar is therefore given by: uh⋆ = N ⋆ a⋆e = N ⋆ qae ,
(6.6)
where N ⋆ contains the shape functions for a one-dimensional bar (see equation (4.19) for a linear element). The ‘star’ denotes that the shape functions are constructed relative to the coordinate system on the bar. The strain in the x ⋆ -direction is therefore given by: ⋆ h⋆ ǫh⋆ = u,x ⋆ = B qa e
(6.7)
Similarly, the derivative of a weight function w can be expressed in the local coordinate system in terms of the global degrees of freedom, ⋆ h⋆ w,x ⋆ = B qb e
(6.8)
(recall that be is a variation). Inserting now the discretised fields into the weak governing equation for a onedimensional bar, and eliminating be , Z
Ωe
( Bbe ) T EBae dΩ =
Z
Γh
( Nbe ) T h dΓ,
(6.9)
where the LHS of equation (6.9) is the stiffness matrix. Inserting the relationships in equations (6.7) and (6.8) into the weak governing equation (2.12), the element stiffness matrix is equal to: ke =
Z
Ωe
q T B⋆ T EB⋆ q dΩ.
(6.10)
This is equivalent to ke = q T k⋆e q, where k⋆e is the standard one-dimensional stiffness matrix for a bar. Clearly, the stiffness matrix for a one-dimensional element can be formed, and with the aid of the matrix q it can be transformed into the global, multidimensional coordinate system. It is important to ensure that a truss structure does not form a mechanism. A classic sign of a mechanism is singularity of the global stiffness matrix of the structure.
January 21, 2011
Version 0.2
75
6 Structural elements for finite element analysis
y fy
Fy x T
x = x1
x = x2
Figure 6.3: Plane beam element Ω = ( x1 , x2 ). v,x
v,x y
x
Figure 6.4: Rotation of a line normal to the axis in a beam segment.
6.2 Beams Two types of beams are considered in this section. The first is the classic BernoulliEuler beam, which is valid for relatively slender beams. The second is the Timoshenko beam, which takes into account shear deformations. The development of beam elements follows the same steps as for continuum elements. The governing equation is identified and then the weak form developed. For simplicity, straight, in-plane beams are considered; there are no out-of plane forces or moments. Such a beam is illustrated in Figure 6.3. Consistent with previous sections, the ‘domain’ of the beam, is denoted Ω = ( x1 , x2 ), and its ends (the boundary of Ω) are denoted Γ = x1 ∪ x2 . The outward normal is denoted n. 6.2.1 Kinematics of a beam Beam theory is based upon the assumption that planes which are normal to the beam’s axis remain plane. For a transverse displacement v, Figure 6.4 shows the rotation of the plane due to a rigid body rotation. Consider now a fibre in a beam
76
Version 0.2
January 21, 2011
6.2 Beams
γ
v,x y
x Figure 6.5: Beam segment subjected to pure shear. The fibre is indicated by the heavy line. which is initially perpendicular to the beam axis. Subjecting a segment from a beam to pure shear, the fibre will not rotate, as illustrated in Figure 6.5. Relative to the fibre which has not rotated, a plane which remains perpendicular to the axis undergoes a rotation γ. For a beam segment subjected to both rotation and shear deformation, the rotation θ of a fibre which is initially perpendicular to the beam’s axis is shown in Figure 6.6. It is clear from Figure 6.6 that θ = v,x − γ.
(6.11)
6.2.2 Equilibrium of a beam The equilibrium equations for a beam can be developed in two ways. The first is to elaborate the kinematics of the beam, which together with some assumptions as to the stress in different directions can be inserted into the elasticity equations from Chapter 1. In the second method, which will be followed here, equilibrium of a beam can be considered directly. The bending moment m in a beam is defined as m=
Z h/2
− h/2
σxx y dy,
(6.12)
and the shear stress q is defined by q=
Z h/2
− h/2
σxy dy,
(6.13)
Figure 6.7 show the ‘resultant’ force Q and moment M due to a moment m and a shear force q in a beam. Note that M = m n and Q = q n, where n is the outward unit normal vector. Considering translational equilibrium of a beam segment Ω, it is clear
January 21, 2011
Version 0.2
77
6 Structural elements for finite element analysis
γ
θ
v,x y
x
Figure 6.6: Rotation of fibre in a beam.
+Q n
n y
+M
x
Figure 6.7: Sign conventions for a bending moment and shear force resultant.
78
Version 0.2
January 21, 2011
6.2 Beams that: Z
dΩ
q n dΓ +
Noting that Z
Ω
R
dΩ
Z
Ω
f y dΩ = 0.
(6.14)
q n dΓ = q| x= L2 − q| x= L1 , it is clear that
q,x dΩ +
Z
f y dΩ = 0.
Ω
(6.15)
Since equilibrium must hold for an infinitely small segment of a beam, translational equilibrium requires that q,x + f y = 0.
(6.16)
For rotational equilibrium, it is required that: Z
dΩ
m n dΓ −
Z
dΩ
q nx dΓ −
Z
Ω
f y x dΩ = 0.
(6.17)
This expression can be rearranged such that Z
Ω
m,x dΩ −
Z
Ω
q dΩ −
Z
Ω
q,x x dΩ −
Z
Ω
f y x dΩ = 0.
Satisfaction of the translational equilibrium equation implies that 0, therefore rotational equilibrium requires that m,x − q = 0.
(6.18)
R
Ω (q,x
+ f y ) x dΩ = (6.19)
The boundary Γ of a beam is partitioned such that: Γv ∪ ΓQ = Γ, Γθ ∪ ΓM = Γ,
Γv ∩ ΓQ = ∅, Γθ ∩ ΓM = ∅
(6.20a) (6.20b)
Denoting applied end forces Fy , distributed loads f y and applied moments T (all shown in Figure 6.3), the boundary conditions are v = gv θn = gθ mn = T q n = Fy
on Γv , on Γθ , on ΓM , on ΓQ ,
(6.21a) (6.21b) (6.21c) (6.21d)
These equations, and boundary conditions, are valid for both Bernoulli-Euler and Timoshenko beam theories.
January 21, 2011
Version 0.2
79
6 Structural elements for finite element analysis 6.2.3 Bernoulli-Euler beam A fundamental assumption in the Bernoulli-Euler theory is that a plane which is initially normal to the longitudinal axis remains a plane and normal to the longitudinal axis. This assumption implies that the shear rotation γ is equal to zero, so the rotation can be directly related to the displacement v, θ=
dv , dx
(6.22)
which also implies that κ=
d2 v . dx2
(6.23)
The bending moment in the beam is related to the curvature by: m = − EIκ = − EI
d2 v , dx2
(6.24)
where I is the moment of inertia of the beam. The negative sign is due to the sign convention in which the rotation and the resultant moment are in opposite directions. Taking now the derivative of all terms in equation (6.19) with respect to x, and then inserting equation (6.16) yields: d2 m + f y = 0. dx2
(6.25)
Assuming EI to be constant, and inserting the constitutive relationship from equation (6.24),
− EI
d4 v + f y = 0, dx4
(6.26)
which is the strong equation of equilibrium for a Bernoulli-Euler beam. Being a fourth-order equation, two boundary conditions are required at both ends of the beam. Dirichlet boundary conditions involve the prescription of the displacement or the rotation (equations (6.21a) and (6.21b)), and Neumann involve either the shear force or moment (equations (6.21c) and (6.21d)). With appropriate boundary conditions, the boundary value problem is complete and can be solved. Weak governing equation Following the procedures from Chapter 2, the weak form of equilibrium for a beam ¯ from an approcan be developed. Multiplying equation (6.25) by a weight function v, priately defined space, which is equal to zero where Dirichlet (kinematic) boundary conditions are applied (the bar is used in this section to denote a weight function) and integration over the beam Ω yields: Z
80
Ω
¯ ,xx dΩ + vm
Z
Ω
v¯ f y dΩ = 0.
(6.27)
Version 0.2
January 21, 2011
6.2 Beams Integrating by parts the term involving the moment M once yields:
−
Z
Ω
v¯,x m,x dΩ +
Z
Γ
¯ ,x n dΓ + vm
Z
Ω
v¯ f y dΩ = 0.
(6.28)
Applying integrating by parts again, this time to the term Z
Ω
v¯,xx m dΩ −
Z
Γ
v¯,x mn dΓ +
Z
Γ
¯ ,x n dΓ + vm
Z
Ω
R
Ω v¯,x m,x
dΩ,
v¯ f y dΩ = 0.
(6.29)
Inserting now the Neumann (natural) boundary conditions from equations (6.21c) and (6.21d), and inserting the constitutive relation in equation (6.24), solving the governing weak equation for a beam involves: find v ∈ S such that
−
Z
Ω
v¯,xx EIv,xx dΩ −
Z
ΓM
v¯,x T dΓ +
Z
ΓQ
¯ y dΓ + vF
Z
Ω
v¯ f y dΩ = 0
∀v¯ ∈ V , (6.30)
where S and V are appropriately defined spaces. The integration domain for the boundary integrals has been changed since v¯ = 0 on Γv , and v¯,x = 0 on Γθ . Note that a fourth-order problem, after integrating by parts twice, has second-order derivatives in its weak form. Consistency of the above weak form can be proven following the same procedure as in in Section 2.1. The existence of second-order derivatives in the weak form deserves some special attention. For the weak form to ‘make sense’, the trial functions (v ∈ S ) and the weight (test) functions (v¯ ∈ V ) must possess a higher degree of regularity (roughly speaking, continuity) than for classical continuum problems. Both S and V must be subspaces of the Sobolev space H 2 (Ω), which contains functions satisfying: Z
Ω
v2 + v2,x + v2,xx dx < ∞.
(6.31)
In one dimension this is equivalent to the trial and test functions being at least C1 continuous, which means that both their first and second derivatives exist. This is in contrast to second-order problems, which require C0 continuity only. For completeness, the spaces S and V are formally defined by: o n (6.32a) S = v | v ∈ H 2 (Ω) , v = gv on Γv , v,x = gθ on Γθ , o n V = v¯ | v¯ ∈ H 2 (Ω) , v¯ = 0 on Γv , v¯,x = 0 on Γθ . (6.32b) For the case of an elastic continuum, it was shown that if the weight function could be considered as a ‘virtual displacement’, it implies that the weak form of the governing equation is identical to the equation of virtual work. Consider now: v¯ ≡ δv
(6.33)
Inserting this relationship into equation (6.30) gives: Z
Ω
δv,xx EIv,xx dΩ = −
Z
ΓM
δv,x T dΓ +
Z
ΓQ
δvFy dΓ +
Z
Ω
δv f y dΩ
(6.34)
which is the equation of virtual work for a Bernoulli-Euler beam beam.
January 21, 2011
Version 0.2
81
6 Structural elements for finite element analysis 6.2.4 Timoshenko beam Timoshenko beams are more general Bernoulli-Euler beams as they allow for shear deformation, and are hence suitable for relatively short beams. They also provide a solid introduction to moderately thick plate theories where the advantage of shear deformable theories will become more apparent. Timoshenko beam theory does not rely on the assumption that a plane which is initially normal to the longitudinal axis remains normal to the longitudinal axis. Planes must remain plane, but not necessarily normal. This means that shear deformations can be taken into account. Recall from equation (6.11) that rotation θ of a plane is given by θ=
dv − γ, dx
(6.35)
where γ is the rotation due to shear. In terms of the shear force and properties of the beam, it is related to the shear force by the constitutive relationship: γ=
q , GAS
(6.36)
where G is the shear modulus and As is the effective shear area. Relative to the Bernoulli-Euler theory, an additional unknown, the shear strain γ, has been introduced. The governing equations are given by the equilibrium conditions in equations (6.19) and (6.16) with the constitutive relationships m = − EIκ, q = GAs γ = GAs
(6.37)
dv −θ , dx
(6.38)
and considering v and θ to be the fundamental unknowns. Inserting the constitutive relationships into the equilibrium equations in (6.19) and (6.16) leads to: d2 θ dv − EI 2 − GAs −θ = 0 dx dx 2 dθ d v + fy = 0 − GAs dx dx2
in Ω,
(6.39)
in Ω,
(6.40)
which are two coupled second-order equations, the first for the rotation θ, and the second for the displacement v. Together with the previously defined boundary conditions, the problem is complete.
82
Version 0.2
January 21, 2011
6.2 Beams Weak governing equations Multiplying equations (6.39) and (6.40) by appropriately defined weight functions θ¯ ¯ and v,
−
Z
Z
Ω
Ω
¯ θEIθ ,xx dΩ −
Z
Ω
¯ θGA s ( v,x − θ ) dΩ = 0,
¯ vGA s (v,xx − θ,x ) dΩ +
Z
(6.41)
v¯ f y dΩ = 0.
Ω
(6.42)
Applying integration by parts once, Z
Ω
−
θ¯,x EIθ,x dΩ −
Z
Ω
Z
Γ
¯ θEIθ ,x n dΓ −
v¯,x GAs (v,x − θ ) dΩ +
Z
Γ
Z
Ω
¯ θGA s (v,x − θ ) dΩ = 0,
¯ vGA s (v,x − θ ) n x dΓ +
Z
Ω
(6.43)
v¯ f y dΩ = 0,
(6.44)
and then inserting the Neumann boundary conditions yields the weak problem for a Timoshenko beam of: find v ∈ Sv and θ ∈ Sθ such that Z
Ω
−
θ¯,x EIθ,x dΩ −
Z
Ω
Z
Ω
¯ θGA s (v,x − θ ) dΩ +
v¯,x GAs (v,x − θ ) dΩ +
Z
¯ y dΓ + vF
ΓQ
Z
Z
¯ dΓ = 0, θT
∀θ¯ ∈ Vθ
(6.45)
v¯ f y dΩ = 0
∀v¯ ∈ Vv .
(6.46)
ΓM
Ω
A more precise definition of the necessary functions spaces should be provided to complete the problem. It is however assumed that θ¯ = 0 on Γθ and that θ satisfies the rotation boundary condition, and v¯ = 0 on Γv , and v satisfies the transverse displacement boundary condition. The test and trial functions should be at least C0 continuous (in contrast to the C1 continuity for the Bernoulli-Euler theory). As with all problems, consistency can be proven through the application of integration by parts. ¯ and using γ = v,x − θ, adding the weak As before, considering δθ ≡ θ¯ and δv ≡ v, forms in equations (6.45) and (6.46) is equivalent to the virtual work equation,
−
Z
Ω
δκ m dΩ +
Z
Ω
δγ q dΩ = −
Z
ΓM
δθ T dΓ +
Z
ΓQ
δv Fy dΓ +
Z
Ω
δv f y dΩ (6.47)
For many mechanical problems, the weak form or virtual work expressions appear naturally and can be used as a starting point in developing a finite element model. 6.2.5 Finite element formulations form beams Euler beams A Galerkin problem for a Bernoulli-Euler beam involves: find vh ∈ S h such that: Z
Ω
h h dΩ = − EIv,xx v¯,xx
January 21, 2011
Z
ΓM
h T dΓ + v¯,x
Z
ΓQ
v¯ h Fy dΓ +
Version 0.2
Z
Ω
v¯ h f y dΩ = 0
∀v¯ h ∈ V h (6.48)
83
6 Structural elements for finite element analysis where S h ⊂ S and V h ⊂ V are finite-dimensional spaces. As for continuum elements, we wish to express the displacement field v in terms of shape functions and nodal degrees of freedom. The problem that arises is that simple C0 finite element shape functions are not suitable. The above equation requires the evaluation of the second derivative of the interpolated displacement field. However, the second derivative of a C0 continuous function does not exist in a classical sense. The use of C0 interpolations for fourth-order problems is not mathematically consistent and can lead to unpredictable results. Crucially, convergence of the solution is not assured for interpolations with insufficient continuity. The solution is to use C1 shape functions. Such shape functions can be constructed relatively easily in one dimension (the extension to multiple dimensions is however far from trivial). Hermitian polynomials are C1 functions, and involve both displacement and rotational degrees of freedom. A Hermitian beam element with two nodes has four degrees of freedom (two displacement degrees of freedom and two rotation degrees of freedom). This results in a cubic interpolation of the displacement along the element. The displacement field vh is given by: 2
vh ( x ) = ∑ ( Ni ( x ) vi + Mi ( x ) θi ) ,
(6.49)
i
where vi and θi are degrees of freedom associated with node i, and the shape functions are equal to:
− ( x − x2 )2 (− h + 2 ( x1 − x )) (6.50a) h3 ( x − x1 )2 ( h + 2 ( x2 − x )) N2 = (6.50b) h3 ( x − x1 ) ( x − x2 ) 2 (6.50c) M1 = h2 ( x − x1 )2 ( x − x2 ) M2 = (6.50d) h2 for an element of length h with ends from x1 to x2 (x2 > x1 ). The displacement at a point in the beam is given by: v 1 θ1 (6.51) vh = Nae = N1 M1 N2 M2 v 2 θ2 N1 =
It is necessary to compute both the first and second derivatives of v with respect to x. The first derivative is given by: v 1 θ1 h (6.52) = N,x ae = N1,x M1,x N2,x M2,x v,x v 2 θ2 84
Version 0.2
January 21, 2011
6.2 Beams The second derivative of v is given by:
h = N,xx ae = N1,xx v,xx
M1,xx
N2,xx
v1 θ1 M2,xx v2 θ2
(6.53)
h and v h can be computed given a , they can be inserted into the Now that vh , v,x e ,xx Galerkin problem,
Z
Ω
Z
( N,xx be ) T EI N,xx ae dΩ = −
ΓM
( N,x be ) T T dΓ +
Z
( Nbe ) T Fy dΓ
ΓQ
+
Z
Ω
( Nbe ) T f y dΩ (6.54)
After some rearranging,
Z
Ω
N,xx T EI N,xx dΩ ae = −
Z
ΓM
N,x T T dΓ +
Z
ΓQ
N T Fy dΓ
+
Z
Ω
N T f y dΩ.
(6.55)
The element stiffness matrix is given by:
ke =
Z
Ω
N,xx T EI N,xx dΩ,
(6.56)
and the RHS vector by:
fe =
Z
ΓQ
N T Fy dΓ −
Z
ΓM
N T T dΓ +
Z
Ω
N T f y dΩ.
(6.57)
The operation to form the RHS vector essentially translates the applied loads into equivalent nodal shear forces and moments.
January 21, 2011
Version 0.2
85
6 Structural elements for finite element analysis Taking derivatives of the Hermitian shape functions in equation (6.50), 2 ( x − x2 ) (3x + h − 2x1 − x2 ) , h3 2 (6x + h − 2x1 − 4x2 ) , N1,xx = h3 ( x − x2 ) (3x − 2x1 − x2 ) M1,x = , h2 2 (3x − x1 − 2x2 ) M1,xx = , h2 2 ( x − x1 ) (3x − h − x1 − 2x2 ) N2,x = − , h3 2 (6x − h − 4x1 − 2x2 ) N2,xx = − , h3 ( x − x1 ) (3x − x1 − 2x2 ) , M2,x = h2 2 (3x − 2x1 − x2 ) . M2,xx = h2 N1,x =
(6.58a) (6.58b) (6.58c) (6.58d) (6.58e) (6.58f) (6.58g) (6.58h)
Assuming the centre of the element is at x = 0 (x2 + x1 = 0, x1 = − h/2), the above equations can be simplified significantly. Considering just the second derivatives with respect to x, 12x , h3 6x 1 M1,xx = 2 − , h h 12x N2,xx = − 3 , h 6x 1 M2,xx = 2 + . h h
N1,xx =
(6.59a) (6.59b) (6.59c) (6.59d)
Inserting these terms into equation (6.56), the stiffness matrix is of the form: 12x h3 6x 1 Z h/2 h2 − h EI 12x ke = 12x h3 − h/2 − 3 h 6x 1 + h h2
6x 1 − h h2
12x − 3 h
6x 1 + h h2
dx.
(6.60)
Integrating the terms in the stiffness matrix exactly from − h/2 to h/2 (assuming EI 86
Version 0.2
January 21, 2011
6.2 Beams to be constant), the element stiffness matrix is equal to:
ke =
12EI h3 6EI 2 h 12EI − h3 6EI h2
6EI h2 4EI h 6EI − h2 2EI h
6EI h2 2EI h . − 6EI h2 4EI h
− 12EI h3 − 6EI h2 12EI h3 − 6EI h2
(6.61)
Timoshenko beams The finite element formulation for the Timoshenko beam is relatively simple, due to the second-order nature of the governing equations. This will allow the use of the standard shape functions introduced in Chapter 4. The Galerkin problem for the Timoshenko beam involves: find vh ∈ Svh and θ h ∈ Sθh such that Z Z Z h h h θ¯ h T dΓ ∀θ¯ h ∈ Vθh , (6.62) − θ h dΩ = − θ¯ h GAs v,x θ¯,x EIθ,x dΩ − ΓM Ω Ω Z Z Z h h v¯ h Fy dΓ ∀v¯ h ∈ Vvh . (6.63) v¯ h f y dΩ + v¯,x GAs v,x − θ h dΩ = ΓQ
Ω
Ω
Consider the representation of the fields vh and θ h , and corresponding derivatives, in terms of shape functions and nodal variable, vh = N v ave , h
θ = h
v¯ = θ¯h =
(6.64a)
N θ aθe , N v bve , N θ bθe ,
(6.64b) (6.64c) (6.64d)
A distinction is made between the shape functions for vh and θ h as these may differ. Inserting the expressions for the unknown fields in terms of nodal variables and variations into equations (6.62) and (6.63) leads to: Z
Z
T
Ωe
Ωe
Bθ EIBθ dΩ aθe +
Z
Bv T GAs Bv dΩave −
T
Ωe
Z
N θ GAs N θ dΩ aθe −
Ωe
T
Ωe
N θ GAs Bv dΩ ave
=−
Z
Γe,M
N v T Fy dΓ +
Z
Ωe
T
N θ T dΓ.
(6.65)
N v T f y dΩ,
(6.66)
Bv T GAs N θ dΩaθe
=
January 21, 2011
Z
Z
Γe,Q
Version 0.2
87
6 Structural elements for finite element analysis The stiffness matrix of an element is of the form θθ k kθv aθ f = θ fv kvθ kvv av
(6.67)
where the components of the the element stiffness matrix are given by kθθ =
Z
T
Ωe
kθv = − kvθ = − kvv =
Z
Z
Z
Ωe
T
Bθ EIBθ + N θ GAs N θ dΩ,
(6.68)
T
Ωe Ωe
N θ GAs Bv dΩ,
(6.69)
Bv T GAs N θ dΩ,
(6.70)
Bv T GAs Bv dΩ,
(6.71)
and components of the element RHS vector are given by fθ = − fv =
Z
Z
T
Γe,M
Γe,Q
N θ T dΓ,
N v T Fy dΓ +
(6.72) Z
Ωe
N v T f y dΩ.
(6.73) T
Note that the stiffness matrix is symmetric (kθv = kvθ ). Locking Conceptually, Timoshenko beam theory is superior to the Bernoulli-Euler theory in that it is valid for both short and slender beams. Ideally, only a Timoshenko beam element would be necessary in practice, and both short and slender beams could be analysed. This would be advantageous from the point of view of generality, and for simplicity as the Timoshenko element used C0 shape functions. However, when applying Timoshenko elements for thin bending problems, the result are plagued by shear locking. This is the effect in which the shear contribution to the energy does not vanish, resulting in an overly stiff response. Consider a Timoshenko beam element of length L, with linear shape functions for both the displacement and rotation. The components of the stiffness matrix are given by: GAs L EI GAs L EI − + L + 3 L 6 , kθθ (6.74) e = GAs L GAs L EI EI + − + L 6 L 3 GAs GAs − 2 2 , (6.75) kθv e = GAs GAs − 2 2 88
Version 0.2
January 21, 2011
6.2 Beams and
kvv e
GAs L = GA s − L
GAs − L GAs L
(6.76)
Therefore, the stiffness matrix is of the form
EI GAs L L + 3 EI GAs L + − 6 L ke = GAs 2 GAs − 2
EI GAs L + L 6 EI GAs L + L 3 GAs 2 GAs − 2
−
GAs 2 GAs 2 GAs L GAs − L
GAs 2 GAs − 2 GAs − L GAs L
−
(6.77)
If at one end the beam is clamped, θ1 = v1 = 0, and a shear load P is applied one end,
EI GAs L L + 3 ke = GAs − 2
GAs θ2 0 2 = GAs v2 P L
−
(6.78)
Therefore, the displacement v2 is equal to v2 =
P ( EI/L + GAs L/3) EIGAs /L2 + G2 A2s /12
(6.79)
If L is large, v2 ≈
4PL , GAs
(6.80)
which is independent of I, will yield a very stiff response. As L → 0, v2 ≈
PL , GAs
(6.81)
which is the correct result. Even with a very fine mesh, this element will exhibit a very stiff response. R T If the term Ωe N θT GAs N θ dΩ is evaluated using one-point numerical integration
January 21, 2011
Version 0.2
89
6 Structural elements for finite element analysis (which is reduced for this term), the element stiffness matrix would have the form EI GAs L EI GAs L GAs GAs + − + − L 4 L 4 2 2 EI GAs L GAs L GAs EI GAs + + − − 4 L 4 2 2 L (6.82) ke = GAs GAs GAs GAs − 2 2 L L GAs GAs GAs GAs − − − 2 2 L L Now, for the one element problem clamped at one end, v2 is equal to v2 =
P ( EI/L + GAs L/4) EIGAs /L2
(6.83)
If L is very large, v2 ≈
PL3 4EI
(6.84)
which qualitatively the desired response, and as L → 0, v2 ≈
PL GAs
(6.85)
When using reduced integration, this element performs quite well when the mesh is sufficiently well-refined. The case is similar for quadratic elements, where selective integration (two-point) leads to a dramatic reduction in locking response.
90
Version 0.2
January 21, 2011
6.3 Plate x3
p θ1 x2 t
x1 θ2
Figure 6.8: Coordinate system for a plate.
6.3 Plate Plates are effectively generalisations of beam elements to two spatial dimensions. They are flat two-dimensional entities, are are typically loaded out of plane. Plate finite elements can be extremely complex and are an area of ongoing research. A serious problem stems from the need for C1 interpolations in the classic thin plate and shell bending theories. Given the enormous number of different plate and shell elements, this section does not provide an exhaustive coverage. The basic governing equations of shell and plate elements and the most common finite element approaches are discussed. Also, the relevant dangers in analysing plate and shell structures with the finite element method are presented. More detailed coverage can be found in a number of books (Zienkiewicz and Taylor, 1991; Hughes, 1987; Bathe, 1996). A particularly accessible account can be found in Cook et al. (1989). 6.3.1 Plate kinematics The adopted coordinate system and sign convention for rotations θα for plates is shown in Figure 6.8. For convenience, the convention differs from the right-hand rule. The adopted definitions simplify the formulation as it will avoid the need for a −1 factor on particular terms. The surface of the plate in the x1 –x2 plane is denoted Ω, and the boundary is Γ. The three coordinates are denoted x, y and z, or x1 , x2 and x3 . Similarly, displacement may be denoted u, v and w (displacements in the x, y and z directions, respectively), or u1 , u2 and u3 when using indexes. It is convenient to express many relationships for plates using index notation. T When using Greek subscripts, summation is implied from one to two. The plate has a thickness t, and the coordinate z = 0 corresponds to the mid-surface of the plate. In plate theory, the displacement of a point is given by: uα = −zθα ( x, y) ,
(6.86a)
u3 = u3 ( x, y) ,
(6.86b)
where θα is the rotation of a fibre in the plate (see Figure 6.9). The rotation of the
January 21, 2011
Version 0.2
91
6 Structural elements for finite element analysis
γα θ α
w,α 1
w
x3
xα
Figure 6.9: Plate kinematics with shear deformation included.
92
Version 0.2
January 21, 2011
6.3 Plate mid-surface of the plate is given by w,α , and the shear strain is given by γα , γα = w,α − θα .
(6.87)
Note that if γα = 0, then θα = w,α . From the displacement field in equations (6.86), the strain field can be calculated at any point (recall ǫ = (1/2)(ui,j + u j,i )). −z θα,β + θ β,α 2
ǫαβ =
(6.88a)
−θα + w,α 2
ǫα3 =
(6.88b)
A useful quantity is the curvature, which is a second-order tensor and at a point is given by: καβ =
1 θα,β + θ β,α = θ(α,β) . 2
(6.89)
6.3.2 Equilibrium of a plate The equilibrium equations for a plate are developed here in the same fashion as for plates. The alternative approach would be to insert the plate kinematics into the elasticity equations. The moment tensor mαβ in a plate is defined as: Z t/2
−t/2
σαβ z dz
(6.90)
From symmetry of the stress tensor, it is clear that the moment tensor is symmetric. The shear force vector qα is defined as: Z t/2
−t/2
σα3 dz
(6.91)
The different resultant moments and transverse forces acting on cross-sections of a plate are shown in Figure 6.10 (‘resultants’ are forces which are equivalent to a moment or shear stress). It is useful to define a vector s which is normal to the vector n, and lies in the plane of Ω. This allows the definition of following notation for
January 21, 2011
Version 0.2
93
6 Structural elements for finite element analysis q1
q2 m11
m21
m12 m22
m21
m22 m11 x3
m12 x2
x1
q1
q2
θ1
θ2 Figure 6.10: Moment and transverse shear force resultants acting on cross-sections of a plate (adapted from Hughes (1987)). quantities relative to an arbitrarily oriented plate edge: mnn = mαβ n β nα ,
(6.92a)
mns = mαβ n β sα ,
(6.92b)
∂mnn , ∂s ∂mns , mns,s = ∂s qn = qα nα , θn = θα nα , w,n = w,α nα w,s = w,α sα mnn,s =
(6.92c) (6.92d) (6.92e) (6.92f) (6.92g) (6.92h)
The quantities are illustrated in Figure 6.11. Translation equilibrium of a plate requires that: Z
∂Ω
qα nα dΓ +
Z
Ω
pz dΩ = 0.
(6.93)
Applying the divergence theorem to the term involving qα leads to Z
Ω
qα,α dΩ +
Z
Ω
pz dΩ = 0,
(6.94)
which must also hold for an infinitesimally small piece of the plate, therefore qα,α + pz = 0,
94
(6.95)
Version 0.2
January 21, 2011
PSfrag
6.3 Plate q1
q2 m11
m21
m12 s
m22
mnn qn
mns
x3
x1
x2
n
θ1
θ2 Figure 6.11: Moment and transverse shear force resultants acting on an edge of a plate (adapted from Hughes (1987)). is the equation for translational equilibrium. Considering now rotational equilibrium, Z
∂Ω
mαβ n β dΓ −
Z
∂Ω
q β n β xα dΩ −
Z
Ω
pz xα dΩ = 0,
(6.96)
which leads to: Z
Ω
mαβ,β dΩ −
Z
Ω
qα dΩ −
Z
Ω
q β,β xα dΩ −
Z
Ω
pz xα dΩ = 0,
(6.97)
noting that q β xα ,β = q β xα,β + q β,α xα = qα + q β,α xα . Since translation equilibrium requires that q β,β + pz = 0, rotational equilibrium requires that mαβ,β − qα = 0.
(6.98)
The boundary conditions for a plate are more complicated that for a beam, and are specific to the considered plate theory, hence boundary conditions are elaborated upon in the context of a given theory. To complete the problem, constitutive relationships are required which relate the moment tensor and the shear force to deformations. These constitutive relationships will also be elaborated upon in the context of the different theories. 6.3.3 Thin plates - Kirchhoff theory The Kirchhoff plate model theory is analogous to the Bernoulli-Euler beam. It is applicable for thin plates where shear deformations are negligible. Ironically, the outwardly simplest plate theory presents the most problems when formulating the finite
January 21, 2011
Version 0.2
95
6 Structural elements for finite element analysis element method. Assuming transverse shear deformations are zero, equation (6.87) reduces to: θα = w,α
(6.99)
which corresponds to the classical theory for beams and plates which assumes fibres which are initially normal to the mid-surface remain normal to the mid-surface. The shear force qα is no longer ‘independent’. The goal is now to eliminate qα from the governing equations. This can be done by differentiating equation (6.98) and inserting equation (6.95), mαβ,αβ + pz = 0.
(6.100)
Since shear deformations are zero, the curvature κ is equal to: καβ =
1 w,αβ + w,βα . 2
(6.101)
The moment is related to the curvature by mαβ = −Cαβγδ κγδ ,
(6.102)
where the constitutive tensor is equal to:
Cαβγδ =
t3 ¯ αβ δγδ . µ δαβ δβδ + δαδ δβγ + λδ 12
(6.103)
where λ¯ = 2λµ/ (λ + 2µ). Inserting the constitutive relationship (6.102) into equation (6.100) yields a fourthorder partial differential equation. After some rearranging, the governing equation can be expressed in a particularly convenient form: 12 1 − ν2 w,ααββ = pz . (6.104) Et3 Expanding this equation, 12 1 − ν2 ∂4 w ∂4 w ∂4 w +2 2 2 + 4 = pz . ∂x y Et3 ∂x4 ∂y
(6.105)
This equation is known as the ‘biharmonic equation’. It is often written in a shorthand format using the biharmonic operator ∇4 ,
∇4 w =
∂4 w ∂4 w ∂4 w + + 2 . ∂x2 y2 ∂x4 ∂y4
(6.106)
As a fourth-order equation, four different types of boundary conditions are possible, and the Kirchhoff plate model requires two boundary conditions at all points
96
Version 0.2
January 21, 2011
6.3 Plate on the boundary. Considering the partition of the boundary in equation (6.20), the boundary conditions are given by: w = gw θ n = gθ mnn = T qn + mns,s = Q
on Γw on Γθ on ΓM on ΓQ
(6.107a) (6.107b) (6.107c) (6.107d)
The first two boundary conditions are Dirichlet boundary conditions, and the last two are Neumann boundary conditions. The last boundary condition is known as the ‘modified shear’ boundary condition. Weak form The weak form of the governing equation for a thin plate is derived using the same processes applied for second-order problems. The governing equation is multiplied by a weight function, then integrated by parts. Rather than multiplying the governing ¯ and integrating over the plate surface Ω, it equation (6.105) by a weight function w, is convenient to carry out the process on equation (6.100). Multiplying by a weight ¯ for which w = 0 on Γw and w,n = 0 on Γθ , and integrating over the surface function w, of the plate, Z
Ω
¯ αβ,αβ dΩ + wm
Z
Ω
¯ z dΩ = 0, wp
(6.108)
where Ω is the surface of the plate. Integrating the moment term by parts once,
−
Z
Ω
w¯ ,α mαβ,β dΩ +
Z
Γ
¯ αβ,β nα dΓ + wm
Z
Ω
¯ z dΩ = 0. wp
(6.109)
Using integration by parts again on the first term, Z
Ω
w¯ ,αβ mαβ dΩ −
Z
Γ
w¯ ,α mαβ n β dΓ +
Z
Γ
¯ αβ,β nα dΓ + wm
Z
Ω
¯ z dΩ = 0, (6.110) wp
which is the weak form of the governing equation for thin plate bending. The problem here is that too many boundary terms appear in the first boundary integral. The moment can only be applied normal to the boundary (mαβ n β nα ),
−
Z
Γ
w¯ ,α mαβ n β dΓ = −
=−
Z
ZΓ Γ
w¯ ,n mnn dΓ − w¯ ,n mnn dΓ +
Z
ZΓ Γ
w¯ ,s mns dΓ ¯ ns,s dΓ − wm ¯ ns | BA wm
(6.111)
The last term will be ignored, and leads to the well known ‘corner forces’ in thin plate theory. Inserting the above relationship into equation (6.110) (ignoring the point terms), yields Z
Ω
w¯ ,αβ mαβ dΩ −
January 21, 2011
Z
Γ
w¯ ,n mnn dΓ +
Z
Γ
w¯ (mns,s + qn ) dΓ +
Version 0.2
Z
Ω
¯ z dΩ = 0. (6.112) wp
97
6 Structural elements for finite element analysis Inserting the Neumann boundary conditions, solving the weak problems then involves: find w ∈ S such that Z
Ω
w¯ ,αβ mαβ dΩ −
Z
ΓM
w¯ ,n T dΓ +
Z
ΓQ
¯ dΓ + wQ
Z
Ω
¯ z dΩ = 0 wp
∀w¯ ∈ V , (6.113)
where the functions in S satisfy the Dirichlet boundary conditions. Considering w¯ as a virtual displacement δw, w¯ ,α as a virtual rotation δθα , and w¯ ,αβ as the virtual curvature, δκαβ , the above equation is the equation of virtual work for a plate: Z
Ω
δκαβ mαβ dΩ −
Z
ΓM
δθn T dΓ +
Z
ΓQ
δwQ dΓ +
Z
Ω
δwpz dΩ = 0
(6.114)
Recall that the bending moment mαβ is a function of the curvature, which in turn involves second derivatives of w, hence the highest order derivatives in the weak governing equation are order two. For completeness, the spaces S and V for the weak form of Kirchhoff plate theory are defined as follows: o n (6.115a) S = w | w ∈ H 2 (Ω) , w = gw on Γw , w,n = gθ on Γθ , o n (6.115b) V = w¯ | w¯ ∈ H 2 (Ω) , w¯ = 0 on Γw , w¯ ,n = 0 on Γθ . Consistency can proven through the repeated application of integration by parts. 6.3.4 Moderately thick plates – Reissner-Mindlin theory The Reissner-Mindlin theory for plates includes shear deformation. It is analogous to the Timoshenko beam. As such, it is applicable for both thick and thin plates. However, some serious practical problems arise in finite element analysis when applying Reissner-Mindlin elements for the analysis of thin plates. For this theory, the governing equations for equilibrium are the same as for the Kirchhoff theory. The difference lies in the kinematics and the constitutive model. The rotation of the plate cannot be calculated simply as a derivative of the transverse displacement. The key to the success of this approach in a finite element context is that the moment tensor mαβ is calculated from the curvature, which depends on θα which is no longer calculated directly from w (θα 6= w,α ). The problem now involves two different unknowns; w and θα . The governing equations for a thick plate are those given in Section 6.3.2. Unlike for thin plates, the governing equations cannot be further reduced since θα is not a direct function of w. The equations of equilibrium are: mαβ,β − qα = 0,
(6.116a)
qα,α + pz = 0.
98
(6.116b)
Version 0.2
January 21, 2011
6.3 Plate In addition to the constitutive relationship relating moment and curvature (equation (6.102)), a constitutive relation relating the shear forces to the shear strain is introduced, qα = Cαβ γβ = Cαβ w,β − θ β , (6.117) where
Cαβ = ktµδαβ ,
(6.118)
and k is a correction factor which is equal to 5/6. This set of equations can be reduced by eliminating qα , mαβ,β − Cαβ w,β − θ β = 0 Cαβ w,αβ − θ β,α + pz = 0
(6.119a) (6.119b)
It is these equations which will be solved. This set of equations cannot be reduced further. This means than the governing equations are being solved in their irreducible form. The structure of the governing equations for Mindlin-Reissner plates theory offers a greater degree of flexibility in applying boundary conditions than the Kirchhoff theory. This is due to θα and w being decoupled. The ‘naturally’ occurring boundary conditions involve w = gw θα = gαθ mαn = Tα qn = Q
on Γw , on Γθ , on ΓM , on ΓQ .
(6.120a) (6.120b) (6.120c) (6.120d)
Given that one of the unknowns, θα , is a vector, three boundary conditions must be applied at all points on the boundary. Weak form To derive the weak form, both equations (6.119a) and (6.119b) are addressed. Firstly, equation (6.119a) is multiplied by a weight function θ¯α , Z
Ω
θ¯α mαβ,β dΩ −
Z
Ω
θ¯α Cαβ w,β − θ β dΩ = 0.
(6.121)
Integrating the above equation by parts,
−
Z
Ω
θ¯(α,β) mαβ dΩ +
Z
Γ
θ¯α mαβ n β dΓ −
Z
Ω
θ¯α Cαβ w,β − θ β dΩ = 0.
(6.122)
Inserting now Neumann boundary conditions (θ¯α = 0 where θα is prescribed) yields:
−
Z
Ω
θ¯α,β mαβ dΩ +
January 21, 2011
Z
ΓM
θ¯α Tα dΓ −
Z
Ω
θ¯α Cαβ w,β − θ β dΩ = 0,
Version 0.2
(6.123)
99
6 Structural elements for finite element analysis which is the weak form of equation (6.119a). Multiplying now equation (6.119b) by a ¯ weight function w, Z
Z ¯ z dΩ = 0. ¯ αβ w,αβ − θ( β,α) dΩ + wp wC
(6.124)
Ω
Ω
Integrating the above equation by parts yields:
−
Z
Ω
Z Z ¯ αβ γβ nα dΓ + ¯ z dΩ = 0. w¯ ,α Cαβ w,β − θ β dΩ + wC wp Γ
Ω
(6.125)
Inserting again Neumann boundary conditions (w¯ = 0 where w is prescribed),
−
Z
Ω
Z w¯ ,α Cαβ w,β − θ β dΩ +
ΓQ
¯ dΓ + wQ
Z
Ω
¯ z dΩ = 0, wp
(6.126)
which is the weak form of equation (6.119b). The weak problem therefore involves: find θα ∈ Sθ and w ∈ Sw such that:
−
Z
Ω
θ¯α,β mαβ dΩ +
Z
ΓM
θ¯α Tα dΓ
−
−
Z
Ω
Z
Ω
θ¯α Cαβ w,β − θ β dΩ = 0,
Z w¯ ,α Cαβ w,β − θ β dΩ +
ΓQ
∀θ¯α ∈ Vθ
(6.127a)
∀w¯ ∈ Vw
(6.127b)
¯ dΓ wQ
+
Z
Ω
¯ z dΩ = 0 wp
Note that both weak equations contains derivatives no higher than order one. It is equations (6.127a) and (6.127b) that are solved using the finite element method. This will require the interpolation of both the transverse displacement w and the rotation θα . However, the highest order derivative appearing for both terms is one. Hence, C0 shape functions can be applied. Combining both equations (multiplying equation (6.127b) by minus one), adopting ¯ δκαβ = θ(α,β) and δγα = w¯ ,α − θ¯α , and taking advantage of the the notation δw = w, symmetry of mαβ , yields: Z
Ω
−δκαβ mαβ dΩ +
Z
Ω
δγα qα dΩ +
Z
ΓM
δθα Tα dΓ −
Z
−
ΓQ
Z
Ω
δwQ dΓ δwpz dΩ = 0 (6.128)
which is the equation of virtual work for a Reissner-Mindlin plate. Note that different to the virtual work equation for a KirchhoffR plate (6.114), the above equation includes a contribution to the energy due to shear ( Ω δγα qα dΩ). 100
Version 0.2
January 21, 2011
6.3 Plate 6.3.5 Plate finite elements As in elasticity, symmetry of mαβ and καβ allows the use of simplified ‘engineering’ notation. m = −Db κ
(6.129)
where m = (m11 , m22 , m12 ) T and κ = (κ11 , κ22 , 2κ12 ) T . The matrix D b has the form:
Et3 D = 12 (1 − ν2 ) b
1
ν 0
ν
0
1
0
0
(1 − ν ) 2
.
(6.130)
For the shear forces qα , q = D s γ,
(6.131)
where
tkµ D = 0 s
0 . tkµ
(6.132)
Kirchhoff plate elements There are few reliable, flexible and robust Kirchhoff plate elements. The complications are rooted in the need for C1 continuity. The majority of Kirchhoff plate elements are non-conforming – this means that they do not satisfy C1 continuity. However, some non-conforming elements yield satisfactory results. In place of requiring C1 continuity, a less severe requirement is that elements satisfy a ‘patch test’. A patch test is a numerical test which tests the ability of an element to reproduce particular results which provides and indication as to whether or not a particular element is suitable for analysis. In particular, it tests whether an element can represent rigid body modes and a state of constant strain. In the context of plates, constant strain implies constant curvature. A common Kirchhoff plate element is rectangular and has four nodes and twelve degrees of freedom. It does not however satisfy C1 continuity. The element is still useful as it passes the patch test. Its shape functions are of the form: N = c1 x3 y + c2 xy3 + c3 x3 + c4 y3 + c5 x2 y + c6 xy2 + c7 x2 + c8 y2 + c9 xy + c10 x + c11 y + c12 .
(6.133)
The applicability of this element is however limited by its inability to satisfy the patch test as a quadrilateral – the element can only be used as a rectangle. There exists a range of Kirchhoff plate elements, none of which are of general applicability.
January 21, 2011
Version 0.2
101
6 Structural elements for finite element analysis Reissner-Mindlin plate elements Plate finite elements based on the Reissner-Mindlin principle do not require C1 continuity. Two fields are interpolated (the transverse displacement and the rotation), with derivatives of order no higher than one for both fields appearing in the weak form. This is a significant advantage. It is then possible to form different elements using the same procedures as for plane and three-dimensional continuum elements. Isoparametric mappings can be used to construct arbitrary shapes. The procedures developed in Chapter 4 for continuum elements can be applied to develop plate elements. When using C0 shape functions with full integration, convergence requirements are guaranteed. The unknown fields, w h and θh are interpolated using C0 shape functions in terms of nodal unknowns, w h = N w aw h
(6.134)
θ θ
θ =N a
(6.135)
Taking derivatives, h w,α = Bw aw h
(6.136)
θ θ
κ =B a
(6.137)
Inserting the discretised fields into equations (6.127a) and (6.127b) gives: Z
T
Ωe
Bθ D b Bθ dΩ aθe −
Z
T
Ωe
N θ CBw dΩaw e +
Z
T
Ωe
N θ CN θ dΩaθe
=−
Z
Ωe
Bw T CBw dΩaw e −
Z
Ωe
Z
NθT T dΓ
(6.138a)
N w T pz dΩ
(6.138b)
ΓM
Bw T CN θ dΩaθe
=
Z
ΓQ
N w T Q dΓ +
Z
Ωe
for an element. This can be expressed as: 102
kww e kθw e
kwθ e kθθ e
w f aw e = eθ , fe aθe
(6.139)
Version 0.2
January 21, 2011
6.3 Plate where
= kww e
Z
Ωe
kwθ e =− kθw e =− kθθ e = f ew =
Z
Z
Z
Z
Ωe
ΓQ
f eθ = −
Z
Bw T CBw dΩ
Ωe
(6.140a)
Bw T CN θ dΩ
(6.140b)
T
Ωe
N θ CBw dΩ T
Bθ D b Bθ dΩ + N w T Q dΓ +
Z
Ωe
(6.140c) Z
T
Ωe
N θ CN θ dΩ
N w T pz dΩ
T
ΓM
NθT T dΓ
(6.140d) (6.140e) (6.140f)
Note again that the element stiffness matrix is symmetric. These are the equations that are assembled into a global stiffness matrix, and solved. Shape functions of different orders can be formulated (even different orders for the transverse displacement and the rotation) and numerical integration is applied. While different elements can be used which will result in a consistent formulation, the performance of different elements may vary drastically. It is possible to simulate thin plates using elements based on the Reissner-Mindlin formulation. It is of course attractive if both thick and thin plates can be analysed using the same elements. To do this, the response of Reissner-Mindlin elements in the limit (t → 0) should be studied. The result should converge to the thin plate theory. This is not the case in practice, and is discussed in the following section. Mindlin-Reissner plate elements for thin plate problems Ideally, thick plate formulations could be used for analysing thin plates. The problem which is confronted is shear locking. As t → 0, Reissner-Mindlin elements typically show an overly stiff response, compared to exact thin plate solutions. The reason for this is that the contribution of the transverse shear deformation to the energy does not vanish. The locking effect can be so extreme as to render elements unusable. For this reason, extensive research efforts have been aimed at developing Reissner-Mindlin elements which do not lock in the thin plate limit. The simplest solution to shear locking is to use reduced or selective integration. This excludes the spurious contribution of the transverse shear deformation to the bending energy. For example, a four-node quadrilateral with 2 × 2 integration for bending terms and one-point integration for transverse shear terms yields good results for thin plate problems. Using full integration for this element results in extreme shear locking. Many Reissner-Mindlin elements involve reduced or selective integration to avoid shear locking. Another possibility to avoid shear locking is to use a mixed formulation. The governing equations are not reduced to their primal form in terms of w and θα , but three
January 21, 2011
Version 0.2
103
6 Structural elements for finite element analysis fields are interpolated: w, θα and qα . It can be shown however that most mixed formulations are equivalent to reduced or selective integration (Hughes, 1987). Discrete Kirchhoff elements are based on thick shell formulations. What they require is that the Kirchhoff thin plate requirement of zero transverse shear strain is met at a discrete number of points within an element. These elements do not require reduced integration to avoid shear locking. However, some ambiguities arise when applying distributed loads.
6.4 Shell elements Shells exist in a three-dimensional space. They are essentially curved plates. A simple approximation to a shell structure is to use a collection of flat combined plate/membrane elements to approximate the curved surface. Shell elements can also be formulated directly. The three-dimensional nature of the element makes them more complicated than plates. The direct formulation of shell elements is not addressed here. Details of their formulation can be found in numerous references (Bathe, 1996; Hughes, 1987; Cook et al., 1989; Zienkiewicz and Taylor, 1991).
6.5 Exercises 1. For a truss element with ends located at ( x1 , y1 ) and ( x2 , y2 ), find the rotation matrix R that rotates from the ( x, y) system to the convenient ( x ′ , y′ ) system. 2. For a two-node Hermitian beam element, how does the bending moment and the shear forces vary along the length of the element? 3. Sketch the Hermitian shape functions in equation (6.50). 4. Give the stiffness matrix for a two-node beam element in two-dimensions that allows for axial deformations. 5. Give two examples of when finite element analysis with continuum elements is preferable and two examples where beam, plate or shell elements is preferable. Provide a short reasoning for each case. 6. What is the danger in applying reduced or selective integration for plate elements? 7. If a constant distributed load is applied normal to a plate surface, pz , calculate the equivalent nodal forces for an eight-node Reissner-Mindlin element.
104
Version 0.2
January 21, 2011
7 Analysis of the finite element method It was shown in Chapter 3 that the finite element method calculates a solution which is optimal in terms of the energy, given a collection of basis functions. While the solution may be optimal, this does not tell how large the error is. It is of course important to have some idea of the magnitude of the error and how it can be controlled and reduced. Also, it is useful to know ‘how accurate’ different elements are. This is known as the order of convergence – if the size of the element is halved, how much smaller will the error be?
7.1 A priori error estimation A priori error estimation does not provide quantitative information as to the error in a finite element simulation. What it does provide is the order of convergence. The order of convergence indicates how quickly the error reduces upon mesh refinement. The order depends on several factors. One is with which ‘norm’ the error is being measured. For example, is the error in terms of displacements being measured, or the error in terms of the strain. Another factor is the smoothness of the exact solution. For finite element analysis, the order of convergence depends heavily on the type of element. We want to know if we halve the element size (leading to at least a doubling in computational effort), how much smaller will the error be. Recall from Chapter 3 the error e at a point is defined as: e = u − uh ,
(7.1)
where u is the exact solution and uh is the approximate solution. To give an indication of the size of the error, it is necessary to consider a norm of the error, kekn . The subscript n indicates the type of norm. The Sobolev norm, which is the type of norm which will be used here, is given by:
kukn =
n Z
∑
α =0 Ω
α
2
( D u) dΩ
!1/2
,
(7.2)
where D α denotes the α derivative. For example, setting n = 2, the Sobolev norm is equal to:
k u k2 =
Z
Ω
uu + u,i u,i + u,ij u,ij dΩ
1/2
.
(7.3)
This scalar value is a measure of the ‘magnitude’ of a function. The magnitude of the error is quantified by calculating a norm of e. Two norms of particular interest
105
7 Analysis of the finite element method involve n = 0 and n = 1. For example, setting n = 0 gives:
k e k0 =
Z
Ω
(e · e) dΩ
1/2
,
(7.4)
which is a measure of the error in the displacements. Setting n = 1,
k e k1 =
Z
Ω
(e · e + ∇e : ∇e) dΩ
1/2
,
(7.5)
which now includes the gradient of e, thereby introducing a measure of the error in the strain. A Sobolev semi-norm is given by Z
|u|n =
Ω
( D n u)2 dΩ
1/2
,
(7.6)
which involves only a particular order derivative. An important semi-norm is the energy norm. It is defined by:
kuk2E = |u|2r =
Z
Ω
( Dr u)2 dΩ,
(7.7)
where r is the order of the derivatives appearing in the weak form. Consider the nodal interpolate u I , uI =
∑ Ni aiI .
(7.8)
i
The nodal interpolate is equal to the exact solution at the nodes (u − u I = 0 at the nodes), hence it interpolates the exact solution using the provided shape functions. Interpolation theory provides the inequality
ku − u I km ≤ chk+1−m |u|k+1 ,
(7.9)
where k is the polynomial order, c is an unknown constant, and h is a measure of the element size. Note that to ensure convergence, k + 1 > m. Given that a finite solution is known to be the best possible solution from the finite-dimensional space S h , it holds that
ku − uh k E ≤ ku − u I k E .
(7.10)
This means that the finite element solution is at least as accurate as the nodal interpolate in terms if energy. In this way, the question of determining an error estimate is no longer a finite element question, rather a question which can be answered from interpolation theory. Inserting now equation (7.10) into the above expression,
ku − uh kr ≤ Chk+1−r |u|k+1 , 106
(7.11)
Version 0.2
January 21, 2011
7.2 A posteriori error estimation where C is a constant, independent of h. For an elasticity problem, r = 1, hence
ku − uh k1 ≤ Chk |u|k+1 .
(7.12)
This implies that the solution converges in this norm at a rate equal to the polynomial order. For practical purposes, this norm is dominated by the error in the strain, hence it represents the error in the energy. This result is however conditional upon the seminorm |u|k+1 existing. If the exact solution is not sufficiently smooth, there may be no gain in accuracy through increasing the polynomial order. It would be useful to find an estimate in terms of lower norms, such as the error in the displacement, ku − uh k0 . An estimate can be obtained using ‘Nitsche’s trick’ (see Strang and Fix (1973) for details). It leads to the estimate:
ku − uh ks ≤ Ch β |u|k+1 .
(7.13)
where β = min (k + 1 − s, 2(k + 1 − r )). In elasticity, the two quantities of special interest are the displacements and the strains (and hence the stresses). The error in the displacements is given by kek0 , R 1/2 . For k ≥ 1, the displacement error is: which is equal to Ω e · e dΩ
kek0 ≤ Chk+1 |u|k+1 .
(7.14)
In terms of the strain,
kek1 ≤ Chk |u|k+1 ,
(7.15)
which implies that the approximate strains converge to the exact result slower than the approximate displacements. For linear elements, k = 1. Therefore,
kek0 ≤ Ch2 |u|2 kek1 ≤ Ch|u|2 .
(7.16)
Hence linear elements are known as being O(h2 ) in terms of displacements, and O(h) in terms of strains. The notation O(h a ) means that the error is ‘of the order h a ’. Of course as h becomes small, the greater a the smaller the error. Larger a implies faster convergence. The convergence order for different polynomial orders is given in Table 7.1.
7.2 A posteriori error estimation Once a finite element solution has been calculated, a natural question is how close the solution is to the exact solution. This is error analysis, which tells ‘how good’ the calculated solution is. The order of convergence of a particular element type indicates how quickly the exact solution is approached, but does not tell how far the computed solution is from the exact solution.
January 21, 2011
Version 0.2
107
7 Analysis of the finite element method polynomial order 1 2 3 4
displacements k 2 3 4 5
strains k 1 2 3 4
Table 7.1: Order of convergence O(hk ) for commonly used finite elements in terms of displacements and strain, assuming that the exact solution is sufficiently smooth. A posteriori error estimators are generally either residual-based or rely on stress recovery techniques. Fortunately, some relatively simple error estimators exist for linear problems. The most common is the Zienkiewicz-Zhu (ZZ) (Zienkiewicz and Zhu, 1987) stress recovery error estimator. It is based on the property that the stresses at certain points within an element are super-convergent. That is, they converge faster toward the exact solution than other points. The essence of the error estimator is to take the stresses at the super convergent points and form an interpolation of the stress using these points. Then, by comparing the stresses interpolated from the super convergent points with the stress computed from the derivative of the displacement field, the error can be estimated. Error estimation is a rich field of applied mathematical research. Research is ongoing for more sophisticated error estimators for more complicated problems.
7.3 Adaptivity To improve the accuracy of a finite element solution, adaptivity techniques can be employed in combination with an error estimator. Based on an analysis of the error, a finite element mesh can be adapted to improve the solution. Two forms of adaptivity are commonly used: p-adaptivity and h-adaptivity. padaptivity involves increasing the polynomial order of elements to reduce the error in the calculated solution. h-adaptivity uses another technique to reduce the error. It relies upon reducing the size of the element. Hence, a new mesh must be constructed. It is possible to generate an entirely new mesh, or to devise procedures which sub-divide elements where the error is large. A mathematically appealing concept is hp-adaptivity. It has an extraordinary rate of convergence, but is challenging to implement. It may be desirable to reduce the error to a specified value throughout a mesh. Given a estimation of the error, adaptivity can be employed to reduce the error where needed. Then, the problem is re-analysed, the error is again estimated, and adaptivity is employed. The process can be repeated until the prescribed error tolerance is met. Given a measure of the error, an adaptive scheme is required to reduce the error to the prescribed levels. At this point, knowing the order of convergence of an element may prove useful. To halve the error in a particular problem using h-adaptivity, it
108
Version 0.2
January 21, 2011
7.4 Exercises helps to know how quickly the error reduces for a particular element.
7.4 Exercises 1. Give the order of convergence in terms of the strain for the following element: a) b) c) d) e)
T3 Q4 T6 Q8 Q9
2. When using Q4 elements, how much smaller would you expect the error in the displacements and the stresses to be if all the elements are made 2.7 times smaller? 3. List some advantages and disadvantages of p- and h-adaptivity.
January 21, 2011
Version 0.2
109
8 Solvers for large linear systems Note: This chapter is still being developed. The finite element solution of practical problems often leads to very large systems of equations. Large systems require significant memory to store the matrices and large computational effort to solve the system. To find the finite element solution to a problem, the following system needs to be solved: Ka = f
(8.1)
The stiffness matrix K typically has a number of special properties. Due to the compact support of finite element shape functions, the system of equations is sparse. This means that the stiffness matrix contains many zero elements. Furthermore, if nodes are numbered in an efficient manner, the matrix is also banded. This mean that nonzero terms are located close to the diagonal. Another special property is symmetry. For all problems examined in these notes, the stiffness matrix is symmetric. This is a consequence of using a Galerkin formulation. Storing only non-zero elements in the computer memory leads to enormous savings in memory. Further savings can be achieved for symmetric matrices by storing only the diagonal terms and non-zero terms above the diagonal. However, unlike finite difference methods, the exact structure of the stiffness matrix is not constant. It varies for different meshes. Therefore, a general method is required to solve the system of equations. Gauss elimination – number of operations O n3 In finite element code however, Gauss elimination is rarely used. Typically more efficient direct or iterative solvers are used. Direct solvers yield a result with a known number of steps. Iterative solvers approach the exact solution. The number of steps performed depends on the desired accuracy.
8.1 LU-decomposition A commonly used direct solution technique is known as LU-decomposition. It is particularly effective for solving moderately sized problems where the bandwidth of K is limited (as is the case for a well constructed finite element mesh). The system of equations Ka = f is rewritten as: LUa = f
(8.2)
where L is a lower triangular matrix with ones on the diagonal and U is an upper diagonal matrix. The motivation is that it is relatively simple to solve triangular matrices.
111
8 Solvers for large linear systems The first step is to factorise the matrix K in L and U. In a computer implementation, the factorised matrices can be stored in the memory allocated for K. Then, the system of equations: Ly = f
(8.3)
is solved (the ‘result’ being y). This is step is known as ‘forward substitution. The following step, ‘backward substitution’, is then to solve the system: Ua = y
(8.4)
which yields the solution to the original problem, a. The procedure can be implemented in a computer programming surprisingly simply. The following code extracts for solving a banded system are taken from Golub and Loan (1996). for k=1:n-1 for i=k+1:min(k+p,n) K(i,k) = K(i,k)/A(k,k) end for j=k+1:min(k+q,n) for i=k+1:min(k+p,n) K(i,j) = K(i,j) - K(i,k)*K(k,j) end end end There are several variations on this procedure to improve its performance. ‘Pivoting’ is often used to avoid numerical problems associated with large and small terms on the diagonal. Also, for symmetric systems the algorithm can be modified, yielding considerable computational savings.
8.2 Iterative solvers Direct solvers are relatively simple to program and are robust. However, they can become prohibitively expensive for very large problems. An alternative is the use of iterative solvers. Iterative solvers ‘iterate’ toward the exact solution. The process is stopped when the desired tolerance is reached. For large system of equations, iterative solvers are generally faster than direct solvers, although the time required for convergence depends heavily on the nature of the matrix K. Unlike direct solvers, it is not possible to compute exactly the required computational effort before starting the computation. To be added: Gauss-Siedel Krylov subspace methods Conjugate-gradient The efficiency and robustness of iterative solvers can be improved significantly through the use of preconditioning.
112
Version 0.2
January 21, 2011
9 Time-dependent problems Until now, the finite element method has been applied only for time-independent problems. In this chapter, the unsteady heat equation and elastodynamics are considered in which the time dimension must be addressed.
9.1 Elastodynamics The majority of this chapter deals with elastodynamics problems. Two approaches will be elaborated. The first, modal analysis, relies on linearity of the underlying problem and is suited to vibration analysis. The second, less subtle, approach is more general and suited to wave propagation problems. 9.1.1 Governing equation and the weak form for elastodynamics The equation of motion is given by: ρu¨ = ∇ · σ + b
(9.1)
where ρ is the density and u¨ = ∂2 u/∂t2 is the acceleration vector. In mathematical terms, this equation differs significantly from the static case as the governing equation is hyperbolic. Similar to static problems, its is complemented by boundary conditions. On parts of the boundary, time-dependent Dirichlet boundary conditions are prescribed, u=g
on Γg .
(9.2)
On other parts of the boundary, time-dependent Neumann boundary conditions are applied, σn = h
on Γh
(9.3)
Dynamic problems require some ‘extra’ information compared to their quasi-static counterparts. These are initial conditions. Initial conditions correspond to the conditions at time zero (t = 0). Elastodynamics involves the second derivative of the displacement with respect to time, therefore two initial conditions are required: u ( x, 0) = u0 ( x) u˙ ( x, 0) = v0 ( x)
(9.4a) (9.4b)
where v is the velocity (v = ∂u/∂t). All problems up to this chapter were elliptic. With elliptic problems, information travels through the domain immediately. For example, when a load is applied at the
113
9 Time-dependent problems free end of a restrained rod, a reaction force is felt immediately at the other end. The equations for elastodynamics are hyperbolic. With hyperbolic problems ‘information’ travels at finite speeds. If a bar, restrained at one end, is hit with a hammer at the free end, a stress wave will travel through the bar and a force at the restrained end will not be felt until the stress wave arrives. Deriving the weak form for dynamic problems follows the familiar procedure. Multiplying equation (9.1) by a weight function w, which is independent of time, and integrating gives: Z
Ω
ρw · u¨ dΩ =
Z
Ω
w · (∇ · σ ) dΩ +
Z
Ω
w · b dΩ
(9.5)
Integrating by parts, using the same steps outlined in Chapter 2, and inserting Neumann boundary conditions, solving the weak form involves: for a given t > 0, find u ∈ S such that Z
Ω
ρw · u¨ dΩ = −
Z
Ω
∇s w : σ dΩ +
Z
Ω
w · b dΩ +
Z
Γh
w · h dΓ
∀w ∈ V . (9.6)
9.1.2 Semi-discrete Galerkin form To reach a finite element formulation, the Galerkin form must be developed. From the weak from in the previous section, the Galerkin problem involves: find uh ∈ S h such that Z
Ω
ρwh · u¨ h dΩ = −
Z
Ω
∇s wh : σ h dΩ +
Z
Ω
wh · b dΩ
+
Z
Γh
wh · h dΓ
∀wh ∈ V h , (9.7)
which is known as ‘semi-discrete’. This stems from the fact that the dependence on time t has not been addressed in the same fashion as the spatial dependency x. The time dimension is left continuous. To solve dynamic problems using a computer, like for the static case the governing equation must be expressed in terms of nodal unknowns and basis functions. This is again done using finite element shape functions. The velocity and acceleration fields are represented using the shape functions, which when considering the semidiscrete formulations are constant in time. Therefore, the finite-dimensional velocity and acceleration fields in an element are given by: u˙ h = N a˙ e ,
(9.8a)
h
u¨ = N a¨ e ,
(9.8b)
where a˙ e = dae /dt and a¨ e = d2 ae /dt2 . Inserting the discretised displacement field and the expression for the acceleration field into equation (9.7) and rearranging, for an element, Z
114
Ωe
N T ρN dΩ a¨ e = −
Z
Ωe
B T DB dΩ ae +
Version 0.2
Z
Ωe
N T b dΩ +
Z
Γh,e
N T h dΓ.
(9.9)
January 21, 2011
9.1 Elastodynamics This equation is commonly called the ‘semi-discrete’ equation. This equation can be written is a shorthand format as: Me a¨ e + ke ae = f e
(9.10)
where Me is the element mass matrix. Once the element contributions have been assembled into the global mass and stiffness matrices and the global RHS vector, for a given time tn M a¨ n + Kan = f n .
(9.11)
Before this equation can be solved, a strategy is required to deal with the time derivatives of a. 9.1.3 Mass matrix There are several methods to form the mass matrix. The obvious choice from the variational approach is the consistent mass matrix. This is formed similarly to the stiffness matrix, Mecons =
Z
Ωe
ρN T N dΩ.
(9.12)
This is typically formed by integrating numerically, in the same fashion as for the element stiffness matrix. As with the element stiffness matrix, the mass matrix is symmetric and will have off-diagonal terms. Before the variational basis of the finite element method was properly understood, it was not entirely clear how the mass matrix should be formed. A common approach was to form a lumped mass matrix. This corresponds to having discrete, lumped masses at each node. Since finite element shape functions are equal to unity at their node and zero at all other nodes, the lumped mass matrix has non-zero terms only on the diagonal. This property is highly advantageous in combination with some time integration schemes, as will be discussed in the following sections. Since lumped schemes do not follow naturally from the variational formulation, there is some freedom in determining exactly how the lumped mass matrix should be calculated. Perhaps the simplest procedure to form lumped mass matrices is to form the matrices using numerical integration, but with integration points located only at element nodes. Since shape functions are equal to unity at their node and zero at all other nodes, this will yield a diagonal mass matrix. It is however possible that nodes will have zero or negative masses. Another procedure is based on summing rows of the consistent mass matrix. lump
Mii
n
=
∑ Mijcons
(9.13)
j =1
where n is the dimension of the mass matrix. All off-diagonal terms are equal to zero. This procedure however can lead to negative masses on the diagonal. Other more
January 21, 2011
Version 0.2
115
9 Time-dependent problems sophisticated lumping schemes are possible which will avoid negative masses. There is no rich theory underlying the choices. The performance of lumped mass matrices is guided primarily by experience. They two key points are that zero and negative masses on the diagonal should be avoided. It will become clear when addressing particular time integration schemes why lumped mass matrices are popular. 9.1.4 Damping Problems in structural dynamics often involve damping. The precise source of damping and its mathematical form is often vague. Damping is often included by adding the term C a˙ to the governing equations, M a¨ + C a˙ + Ka = f ,
(9.14)
where the matrix C represents the sources of damping in the structure. A simple method of introducing damping is to say: C = τM + ψK,
(9.15)
where τ and ψ are user chosen parameters. This is known as Rayleigh damping. The contribution of the stiffness matrix to the damping matrix increases damping with increasing frequency. Conversely, the contribution of the mass matrix increases damping with decreasing frequency. The choice of the parameters τ and ψ is certainly arbitrary, representing the uncertainty in the source of damping. A discussion can be found in Cook et al. (1989). Damping due to the material response can also be included. It can be included through the damping matrix C, or through the constitutive model relating stress and strain. This topic is left to other courses which address non-linear material behaviour.
9.2 Heat equation The heat equation is introduced here as a step to the problems in elastodynamics. Unlike elastodynamics, it involves first order derivatives with respect to time. The steady-state version of the equation was introduced in Section 2.3. In strong form, unsteady heat diffusion is governed by: ρcu˙ + ∇ · q = f ,
(9.16)
where in the context of diffusion of heat, c is the ‘capacity’, u is the temperature, q = −κ ∇u is the heat flux, kappa is the conductivity. The problem requires boundary conditions, as outlined in Section 2.3, and initial conditions. Given that only the first derivative with respect to time is involved, only one type of initial condition may be provided, namely the temperature at t = 0, u ( x, t) = u0 ( x) .
116
(9.17)
Version 0.2
January 21, 2011
9.3 Frequency analysis for elastodynamics The heat equation is parabolic. If one end of a rod is heated, this is immediately felt (an infinitesimal amount) and increases with time. Following the usual processes in developing the weak form, the semi-discrete Galerkin problem for the heat equation involves: for a given t > 0, find uh ∈ S h such that Z
Ω
w h ρcu˙ h dΩ +
Z
Ω
∇wh · κ ∇uh dΩ −
Z
Γh
w h h dΓ =
Z
Ω
w h f dΩ
∀wh ∈ V h (9.18)
In matrix form, this can be expressed at time tn as: M a˙ n + Kan = f n .
(9.19)
9.3 Frequency analysis for elastodynamics If an undamped system is allowed to vibrate freely (no forcing terms), its motion is harmonic. The displacements at the nodes can be described by a = a¯ cos (ωt − α) ,
(9.20)
where ω is the frequency (radians per second). The acceleration is therefore given by: a¨ = −ω 2 a¯ cos (ωt − α) .
(9.21)
Inserting these results into equation (9.11) and rearranging, −ω 2 M + K a¯ = 0,
(9.22)
which is a matrix eigenvalue problem. Non-trivial solutions (a¯ 6= 0) require that: det K − ω 2 M = 0, (9.23)
where ω are the natural frequencies of the system. If stiffness and mass matrices are of dimension n, there are n eigenvalues λ (λ = ω 2 ). The eigenvectors are a, ¯ which are also known as natural modes. For explicit time integration procedures, it is important to know the highest natural frequency of the discrete problem. The difficulty is that the size of the problem in equation (9.23) can be very large. An upper bound to the highest frequency can be found by calculating the frequencies for all individual elements. det ke − ωe2 me = 0 (9.24) The highest frequency from all elements provides an upper bound to the highest frequency of the problem. For a one-dimensional linear element, it is relatively simple to calculate the natural frequencies. It turns out that for a lumped mass matrix, the highest frequency is equal to: s 2 E ωe = (9.25) L ρ √ p for an element of length L. For the consistent mass matrix, ωe = 2 3/L E/ρ.
January 21, 2011
Version 0.2
117
9 Time-dependent problems
9.4 Modal analysis for elastodynamics Modal analysis involves the dividing of the response of a structure into a number of modes. The sum of the modes gives the total displacement. For linear problems, it is then possible to treat each mode individually. Consider a solution a = ψi cos (ωi t + θi ) .
(9.26)
Inserting the above equation into the unforced version ( f = 0) of equation (9.11) , K − ωi2 M Ψi = 0, (9.27)
where Ψi are eigenvectors, and are known as vibration modes, and ωi are the natural frequencies. The nodal displacements can be expressed as a summation of eigenfunctions, a (t) =
∑ αi (t) Ψi
(9.28)
i
where αi (t) is the amplitude of the ith eigenmode. A special property of the eigenmodes Ψi is that they are orthogonal. This means that: ( 1 i = j, T Ψi MΨj = (9.29) 0 i 6= j, and T
Ψi KΨj =
(
ωi2 0
i = j, i 6= j.
(9.30)
Orthogonality implies that the vibration modes do not interact with each other. It is then possible to treat each mode separately. Inserting equation (9.28) into equation (9.11) yields: M ∑ α¨ i Ψi + K ∑ αi Ψi = f
(9.31)
i
i
Pre-multiplying both sides of this equation by Ψj T , Ψj T M ∑ α¨ i Ψi + Ψj T K ∑ αi Ψi = Ψj T f i
(9.32)
i
The orthogonality property means that the equation has been decoupled into nm simple, separate equations, α¨ i + ωi2 αi = Ψi T f .
118
(9.33)
Version 0.2
January 21, 2011
9.5 Time stepping algorithms Now, for each mode of interest, the above equation can be solved. Typically, this equation can be solved analytically. For the unforced case, equation (9.33) is a homogeneous ordinary partial differential equation, with the solution αi = c1 sin (ωi t) + c2 cos (ωi t) .
(9.34)
Solutions for forced cases are dependent on the type of loading. In applying modal analysis, the computational effort lies in finding the eigenvectors and eigenvalues. In structural analysis, usually only the low frequency modes have a significant influence on the structural response. Therefore, equation (9.33) is only solved for a few modes. Modal analysis relies on the orthogonality property of the eigenfunctions. Therefore it is limited to linear analysis. For non-linear problems, the orthogonality property does not hold. Further discussion of modal methods can be found in Craig (1981) and Cook et al. (1989).
9.5 Time stepping algorithms Time stepping algorithms involve schemes to increment the time in discrete steps ∆t. Given the values of an , a˙ n and a¨ n at time tn (and perhaps also at earlier times steps), the values an+1 , a˙ n+1 and a¨ n+1 at time tn+1 = tn + ∆t are needed. To do this, a time stepping scheme is needed. Time is then incremented until the desired time is reached. Time stepping schemes are distinguished by being either explicit or implicit. Explicit integrators rely only on information at time step n, and perhaps earlier steps (such as n − 1) to compute the solution at time step n + 1. Implicit schemes rely on information at time step n + 1, and perhaps information at time step n and earlier steps. This means a system of equations must be solved for implicit schemes. Hence, implicit schemes generally require significantly extra computational effort, compared to explicit schemes. Accuracy is not the difference between explicit and implicit schemes. Explicit schemes can be constructed to have a high degree of accuracy (such as the 4th order Runge-Kutta method, which is O(∆t4 ) accurate). Conversely, implicit schemes may have a low degree of accuracy (such as the backward Euler scheme, which is O(∆t) accurate). The difference between explicit and implicit schemes lies in stability. Stability means the range of time steps ∆t for which an integrator will calculate a stable result (small errors do not grow exponentially). If an integrator is unstable, it will ‘blow-up’, rapidly diverging from any sensible result. There are known schemes of exceptional accuracy, but due to their extremely poor stability properties they are useless. There are no known unconditionally stable explicit schemes. 9.5.1 One-step algorithms for the heat equation To introduce time stepping algorithms, the unsteady heat equation is first considered. It provides a step towards the more complicated schemes for elastodynamics.
January 21, 2011
Version 0.2
119
9 Time-dependent problems A popular family of algorithms are ‘generalised trapezoidal’ methods. They involve: an+1 = an + ∆t ((1 − θ ) a˙ n + θ a˙ n+1 ) ,
(9.35)
where 0 ≤ θ ≤ 1. Well known methods are the explicit forward Euler (θ = 0) method, and the implicit trapezoidal (θ = 1/2) and backward Euler (θ = 1) methods. In practice, the forward Euler is not particularly useful as its stability properties are poor. For the problems to be considered here, the trapezoidal method is particularly attractive as it is unconditionally stable, and is second-order accurate. The backward Euler method is also unconditionally stable, but only first-order accurate. In matrix form, the heat equation at time tn+1 must satisfy: M a˙ n+1 + Kan+1 = f n+1 .
(9.36)
Rearranging equation (9.35), a˙ n+1 =
1 1 (1 − θ ) a − an − a˙ n , θ∆t n+1 θ∆t θ
and inserting the above result into equation (9.36) gives, 1 1 (1 − θ ) M + K a n +1 = f n +1 + Man + M a˙ n . θ∆t θ∆t θ
(9.37)
(9.38)
In compact notation, the problem to solve an+1 has become ˆ n+1 = fˆn+1 , Ka
(9.39)
where 1 Kˆ = M+K θ∆t
(9.40a)
1 (1 − θ ) fˆn+1 = f n+1 + Man + M a˙ n , θ∆t θ
(9.40b)
and a˙ n is computed from equation (9.37). 9.5.2 One-step algorithms for elastodynamics The most commonly used time integration scheme in solid and structural mechanics is the generalised Newmark scheme. Newmark time integrators are a family of schemes, both explicit and implicit. Choosing parameters appropriately for a Newmark algorithm yields a number of different, well known integrators. The Newmark algorithm gives an+1 and a˙ n+1 as:
120
∆t2 ((1 − 2β) a¨ n + 2β a¨ n+1 ) 2 = a˙ n + ∆t ((1 − γ) a¨ n + γ a¨ n+1 )
an+1 = an + ∆t a˙ n +
(9.41a)
a˙ n+1
(9.41b)
Version 0.2
January 21, 2011
9.5 Time stepping algorithms where β and γ are parameters which define the nature and properties of the algorithm. Appropriate choices of β and γ yield well known integration methods. The Newmark method is unconditionally stale for: 2β ≥ γ ≥
1 2
(9.42)
For γ = 1/2, a Newmark scheme is second-order accurate. A commonly used time integrator is the central difference method. It is necessary to express the terms a˙ and a¨ in terms of a. Calculating the necessary terms using central differences, a˙ n =
a n +1 − a n −1 2∆t
(9.43a)
a¨ n =
an−1 − 2an + an+1 . ∆t2
(9.43b)
It is possible to rearrange these terms and show that they are equivalent to equation (9.41) for β = 0 and γ = 1/2. It is clear that this algorithm is explicit – it is possible to express an+1 in terms of quantities at n and earlier times (which are known). This scheme is O(∆t2 ) accurate, but it is conditionally stable. It is stable for: ∆t ≤
2 ωh
(9.44)
where ω h is the highest natural frequency (highest value from det K − ω 2 M = 0. The term ω h can be p estimated for linear elements as 2c/L, where c is the dilatational wave speed (c = E/ρ) and L is the length of an element. Therefore, ∆t ≤
L , c
(9.45)
which is the Courant, Friedrichs and Lewy (CFL) stability condition. Another commonly used time integrator is based on the trapezoidal rule, ∆t ( a¨ n + a¨ n+1 ) , 2 ∆t = an + ( a˙ n + a˙ n+1 ) . 2
a˙ n+1 = a˙ n +
(9.46a)
a n +1
(9.46b)
Inserting equation (9.46a) into equation (9.46b), an+1 = an + ∆t a˙ n +
∆t2 ( a¨ n + a¨ n+1 ) 4
(9.47)
This integrator is implicit. It corresponds to β = 1/4 and γ = 1/2 in the Newmark family of integrators. This method is unconditionally stable and the accuracy is O(∆t2 ). Its stability properties are attractive, but as will be shown later, it requires the factorisation of a large matrix. It is also energy conserving.
January 21, 2011
Version 0.2
121
9 Time-dependent problems Explicit time integration – central difference approach The system of equations to be solved, with damping, has the appearance: M a¨ n + C a˙ n + Kan = f n .
(9.48)
It is assumed that a and its time derivatives are known at time t and unknown at time t + ∆t (corresponding to n + 1). Inserting the central difference expressions from equation (9.43), M
an−1 − 2an + an+1 a − a n −1 + Kan = f n . + C n +1 2 2∆t ∆t
(9.49)
Rearranging the above equation such that all terms at t + ∆t are on the LHS and all terms at time t are on the RHS, 1 1 1 1 M + C an+1 = 2 M (2an − an−1 ) + Ca − Kan + f n . (9.50) 2 2∆t 2∆t n−1 ∆t ∆t Recall that the displacement vector a is known at step n and all previous steps, and unknown at step n + 1. Given an and an−1 , it is possible to calculate an+1 . The above equations can be written in a shorthand format as: ˆ n+1 = fˆn , Ma
(9.51)
where ˆ = M
1 1 C M+ 2 2∆t ∆t
(9.52)
and 1 1 fˆn = 2 M (2an − an−1 ) + Ca − Kan + f n . 2∆t n−1 ∆t
(9.53)
Then, the displacement vector a at step n + 1 is given by: ˆ −1 fˆn a n +1 = M
(9.54)
If both M and C are diagonal matrices, an+1 can be calculated without solving a system of equations (hence the attractiveness of lumped mass matrices when using explicit time integrators). In practice, it is not necessary to form the stiffness matrix K. The term Kan is equivalent to: Kan =
Z
Ω
B T σn dΩ.
(9.55)
A difficulty with explicit central-difference time integration is the need for information at n − 1 when starting a calculation. A procedure is required to ‘start’ the 122
Version 0.2
January 21, 2011
9.5 Time stepping algorithms algorithm at t = 0. Combining equations (9.43) it is possible to express an−1 in terms of quantities at n (time t = 0), which are the initial conditions. an−1 = an − ∆t a˙ n +
∆t2 a¨ n 2
(9.56)
The vector an−1 can be used to start the time integration procedure. The central difference procedure is popular in finite element analysis. In combination with a lumped mass matrix, it is very efficient and suitable for very large problems. Furthermore, it is second-order accurate. The significant drawback is the conditional stability. The time step must be chosen carefully for a particular discretisation. Implicit time integration – the Newmark family Formulation of a Newmark scheme starts by inserting the expressions for the velocity and accelerations into: M a¨ n+1 + C a˙ n+1 + Kan+1 = f n+1
(9.57)
Rearranging the expressions in equation (9.41) for β 6= 0 to express the accelerations and velocity at step n + 1 in terms of step n and earlier, a¨ n+1 =
1 1 1 a˙ n − a¨ n + a¨ n ( a n +1 − a n ) − β∆t 2β β∆t2
a˙ n+1 = a˙ n 1 1 1 a ˙ − a ¨ + a ¨ − a − a + ∆t (1 − γ) a¨ n + γ ) ( n n n n n +1 β∆t 2β β∆t2
(9.58a)
(9.58b)
These equations are then inserted into equation (9.57). For simplicity, damping is ignored. This yields, 1 1 1 M a˙ n − a¨ n + a¨ n + Kan+1 = f n+1 (9.59) ( a n +1 − a n ) − β∆t 2β β∆t2 Rearranging such that all unknown quantities appear on the LHS, and all known quantities on the RHS, 1 1 1 1 a˙ n + a˙ n − a¨ n + f (9.60) M + K a n +1 = M an + β∆t 2β β∆t2 β∆t2 which in a short-hand notation can be expressed as: ˆ n+1 = fˆn+1 Ka
January 21, 2011
(9.61)
Version 0.2
123
9 Time-dependent problems Since the stiffness matrix K is not diagonal, calculating an requires the solution of a system of equations. At time step n, the vector fˆ can be calculated from an , the terms in f are prescribed and both M and K can be formed. Solving the resulting system of equations yields an+1 . Note that calculating the ‘new’ a˙ n+1 introduces a dependency on the parameter γ. With the Newmark scheme, numerical damping can be introduced by choosing γ > 1/2. Numerical damping aids in controlling undesirable spurious high frequency modes in the solution. However, choosing γ 6= 1/2 reduces the integration scheme to first order accuracy O(t). 9.5.3 α-Method for elastodynamics In dynamic simulations, it is usually desirable to introduce some numerical damping to suppress spurious high-frequency modes. The problem with the Newmark scheme is that the introduction of any numerical damping reduces the order of accuracy from two to one. A number of algorithms have been developed to preserve second-order accuracy while allowing for the introduction of numerical damping. One such algorithm is the α-Method, which is also known as the Hilber-Hughes-Taylor Method (Hughes, 1987). It involves using a Newmark scheme to solve : M a¨ n+1 + (1 − α) Kan+1 − αKan = f n+1+α ,
(9.62)
where f n+1+α = f n+1+θ∆t . If α = 0, the α-Method reduces to the Newmark scheme. If −1/3 < α < 0, γ = (1 − 2α) /2 and β = (1 − α)2 /4, the algorithm is unconditionally stable and second-order accurate. Numerical dissipation can be controlled through the parameter α. Decreasing α increases the dissipation.
9.6 Space-time formulations Time as a degree of freedom.
9.7 Exercises 1. Calculate the consistent and lumped mass matrices for the four-noded quadrilateral element. 2. Derive the mass matrix for a two-noded beam element. 3. Rearrange equation (9.41) to show for β = 0 and γ = 1/2 the scheme is equivalent to a central-difference procedure.
124
Version 0.2
January 21, 2011
References Bathe, K.-J. (1996). Finite Element Procedures. Prentice-Hall, New Jersey. Braess, D. (2001). Finite Elements: Theory, fast solvers, and applications in solid mechanics. Cambridge University Press, Cambridge. Brenner, S. C. and Scott, L. R. (1994). The Mathematical Theory of Finite Elements Methods. Springer-Verlag, New York. Cook, R. D., Malkus, D. S., and Plesha, M. E. (1989). Concept and Applications of Finite Element Analysis. John Wiley & Sons, New York, third edition. Craig, R. R. (1981). Structural Dynamics. John Wiley and Sons, Inc., New York. Ern, A. and Guermond, J. L. (2004). Theory and Practice of Finite Elements. SpringerVerlag, New York. Golub, G. H. and Loan, C. F. V. (1996). Matrix Computations. The John Hopkins University Press, Baltimore. Hughes, T. J. R. (1987). The Finite Element Method - Linear Static and Dynamic Finite Element Analysis. Prentice-Hall, London, England. Reddy, B. D. (1998). Introductory Functional Analysis. Springer, London. Strang, G. and Fix, G. (1973). Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, New Jersey. Zienkiewicz, O. C. and Taylor, R. L. (1989). The Finite Element Method, volume 1. McGraw-Hill Book Company, Berkshire, fourth edition. Zienkiewicz, O. C. and Taylor, R. L. (1991). The Finite Element Method, volume 2. McGraw-Hill Book Company, Berkshire, fourth edition. Zienkiewicz, O. C. and Zhu, J. Z. (1987). A simple error estimator and adaptive procedure for practical engineering analysis. International Journal for Numerical Methods in Engineering, 24:337–357.
125
