SLIDEREAD

I NTRODUCTION TO THE F INITEE LEMENT M ETHODG. P. Nikishkov Lecture Notes. University of Aizu, Aizu-Wakamatsu 965-8580,

2Updated -01-19

Contents1234Introduction51.1What is the finite element method . . . . . . . . . . . . . . . . . . . . . . . . . . . .51.2How the FEM works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .51.3Formulation of finite element equations . . . . . . . . . . . . . . . . . . . . . . . . .61.3.1Galerkin method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .61.3.2Variational formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .8Finite Element Equations for Heat Transfer112.1Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .112.2Finite element discretization of heat transfer equations . . . . . . . . . . . . . . . . .122.3Different Type Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .13FEM for Solid Mechanics Problems153.1Problem statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .153.2Finite element equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .163.3Assembly of the global equation system . . . . . . . . . . . . . . . . . . . . . . . . .18Finite Elements214.1Two-dimensional triangular element . . . . . . . . . . . . . . . . . . . . . . . . . . .214.2Two-dimensional isoparametric elements . . . . . . . . . . . . . . . . . . . . . . . .224.2.1Shape functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .224.2.2Strain-displacement matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . .234.2.3Element properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .244.2.4Integration in quadrilateral elements . . . . . . . . . . . . . . . . . . . . . . .254.2.5Calculation of strains and stresses . . . . . . . . . . . . . . . . . . . . . . . .26Three-dimensional isoparametric elements . . . . . . . . . . . . . . . . . . . . . . . .284.3.1Shape functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .284.3.2Strain-displacement matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . .294.3.3Element properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .304.3.4Efficient computation of the stiffness matrix . . . . . . . . . . . . . . . . . . .304.3.5Integration of the stiffness matrix . . . . . . . . . . . . . . . . . . . . . . . .314.3.6Calculation of strains and stresses . . . . . . . . . . . . . . . . . . . . . . . .314.3.7Extrapolation of strains and stresses . . . . . . . . . . . . . . . . . . . . . . .314.33

CONTENTS456Discretization335.1Discrete model of the problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .335.2Mesh generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .345.2.1Mesh generators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .345.2.2Mapping technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .345.2.3Delaunay triangulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .36Assembly and Solution376.1Disassembly and assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .376.2Disassembly algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .386.3Assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .386.3.1Assembly algorithm for vectors . . . . . . . . . . . . . . . . . . . . . . . . .386.3.2Assembly algorithm for matrices . . . . . . . . . . . . . . . . . . . . . . . . .39Displacement boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . .396.4.1Explicit specification of displacement BC . . . . . . . . . . . . . . . . . . . .406.4.2Method of large number . . . . . . . . . . . . . . . . . . . . . . . . . . . . .40Solution of Finite Element Equations . . . . . . . . . . . . . . . . . . . . . . . . . . .406.5.1Solution methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .406.5.2Direct LDU method with profile matrix . . . . . . . . . . . . . . . . . . . . .416.5.3Tuning of the LDU factorization . . . . . . . . . . . . . . . . . . . . . . . . .436.5.4Preconditioned conjugate gradient method . . . . . . . . . . . . . . . . . . . .446.46.5

Chapter 1Introduction1.1What is the finite element methodThe finite element method (FEM) is a numerical technique for solving problems which are describedby partial differential equations or can be formulated as functional minimization. A domain of interestis represented as an assembly of finite elements. Approximating functions in finite elements are determined in terms of nodal values of a physical field which is sought. A continuous physical problem istransformed into a discretized finite element problem with unknown nodal values. For a linear problema system of linear algebraic equations should be solved. Values inside finite elements can be recoveredusing nodal values.Two features of the FEM are worth to be mentioned:1) Piece-wise approximation of physical fields on finite elements provides good precision even withsimple approximating functions (increasing the number of elements we can achieve any precision).2) Locality of approximation leads to sparse equation systems for a discretized problem. This helps tosolve problems with very large number of nodal unknowns.1.2How the FEM worksTo summarize in general terms how the finite element method works we list main steps of the finiteelement solution procedure below.1. Discretize the continuum. The first step is to divide a solution region into finite elements. Thefinite element mesh is typically generated by a preprocessor program. The description of mesh consistsof several arrays main of which are nodal coordinates and element connectivities.2. Select interpolation functions. Interpolation functions are used to interpolate the field variables over the element. Often, polynomials are selected as interpolation functions. The degree of thepolynomial depends on the number of nodes assigned to the element.3. Find the element properties. The matrix equation for the finite element should be establishedwhich relates the nodal values of the unknown function to other parameters. For this task differentapproaches can be used; the most convenient are: the variational approach and the Galerkin method.4. Assemble the element equations. To find the global equation system for the whole solutionregion we must assemble all the element equations. In other words we must combine local elementequations for all elements used for discretization. Element connectivities are used for the assemblyprocess. Before solution, boundary conditions (which are not accounted in element equations) shouldbe imposed.5

CHAPTER 1. INTRODUCTION612u13u2xx0L2Lx1x2Figure 1.1: Two one-dimensional linear elements and function interpolation inside element.5. Solve the global equation system. The finite element global equation sytem is typically sparse,symmetric and positive definite. Direct and iterative methods can be used for solution. The nodal valuesof the sought function are produced as a result of the solution.6. Compute additional results. In many cases we need to calculate additional parameters. Forexample, in mechanical problems strains and stresses are of interest in addition to displacements, whichare obtained after solution of the global equation system.1.3Formulation of finite element equationsSeveral approaches can be used to transform the physical formulation of the problem to its finite elementdiscrete analogue. If the physical formulation of the problem is known as a differential equation then themost popular method of its finite element formulation is the Galerkin method. If the physical problemcan be formulated as minimization of a functional then variational formulation of the finite elementequations is usually used.1.3.1Galerkin methodLet us use simple one-dimensional example for the explanation of finite element formulation using theGalerkin method. Suppose that we need to solve numerically the following differential equation:ad2 u b 0,dx20 x 2L(1.1)with boundary conditionsu x 0 0dua x 2L Rdx(1.2)where u is an unknown solution. We are going to solve the problem using two linear one-dimensionalfinite elements as shown in Fig. 1.1.Fist, consider a finite element presented on the right of Figure. The element has two nodes andapproximation of the function u(x) can be done as follows:u N1 u1 N2 u2 [N ]{u}[N ] [N1 N2 ]{u} {u1 u2 }(1.3)

1.3. FORMULATION OF FINITE ELEMENT EQUATIONS7where Ni are the so called shape functionsx x1N1 1 x2 x1x x1N2 x2 x1(1.4)which are used for interpolation of u(x) using its nodal values. Nodal values u1 and u2 are unknownswhich should be determined from the discrete global equation system.After substituting u expressed through its nodal values and shape functions, in the differentialequation, it has the following approximate form:ad2[N ]{u} b ψdx2(1.5)where ψ is a nonzero residual because of approximate representation of a function inside a finite element. The Galerkin method provides residual minimization by multiplying terms of the above equationby shape functions, integrating over the element and equating to zero:Z x2x1[N ]T aZ x2d2[N ]{u}dx dx2x1[N ]T bdx 0(1.6)Use of integration by parts leads to the following discrete form of the differential equation for the finiteelement: · Z x2 ·dN T dNx1adxdxZ x2dx{u} x1(T[N ] bdx 01)dua x x2 dx(10)adu x x1 0 (1.7)dxUsually such relation for a finite element is presented as:[k]{u} {f } · Z x2 ·dN T dN[k] dxadxdx( )x1()Z x2dudu10T{f } a x x1[N ] bdx a x x2 01dxdxx1(1.8)In solid mechanics [k] is called stiffness matrix and {f } is called load vector. In the considered simplecase for two finite elements of length L stiffness matrices and the load vectors can be easily calculated:"[k1 ] [k2 ] ({f1 } bL2aL11)1 1 11#(, {f2 } bL211)( 0R)(1.9)The above relations provide finite element equations for the two separate finite elements. A globalequation system for the domain with 2 elements and 3 nodes can be obtained by an assembly of elementequations. In our simple case it is clear that elements interact with each other at the node with globalnumber 2. The assembled global equation system is: 11 10 u1 bL a 2 1 u22 1 L2 10 11 u3 0 0 R (1.10)

CHAPTER 1. INTRODUCTION84u32Exact100.0FEM0.5x1.01.52.0Figure 1.2: Comparison of finite element solution and exact solution.After application of the boundary condition u(x 0) 0 the final appearance of the global equationsystem is: 0100 u1 bL a 2 1 u22 0 L2 10 11 u3 0 0 R (1.11)Nodal values ui are obtained as results of solution of linear algebraic equation system. The value ofu at any point inside a finite element can be calculated using the shape functions. The finite elementsolution of the differential equation is shown in Fig. 1.2 for a 1, b 1, L 1 and R 1.Exact solution is a quadratic function. The finite element solution with the use of the simplest element ispiece-wise linear. More precise finite element solution can be obtained increasing the number of simpleelements or with the use of elements with more complicated shape functions. It is worth noting thatat nodes the finite element method provides exact values of u (just for this particular problem). Finiteelements with linear shape functions produce exact nodal values if the sought solution is quadratic.Quadratic elements give exact nodal values for the cubic solution etc.1.3.2Variational formulationThe differential equationd2 u b 0,dx2u x 0 0dua x 2L Rdxa0 x 2L(1.12)with a EA has the following physical meaning in solid mechanics. It describes tension of theone dimensional bar with cross-sectional area A made of material with the elasticity modulus E andsubjected to a distributed load b and a concentrated load R at its right end as shown in Fig 1.3.Such problem can be formulated in terms of minimizing the potential energy functional Π:ZΠ u x 0µ1dua2dxL 0¶2Zdx Lbudx Ru x 2L(1.13)

1.3. FORMULATION OF FINITE ELEMENT EQUATIONSb19R23L2Lx0Figure 1.3: Tension of the one dimensional bar subjected to a distributed load and a concentrated load.Using representation of {u} with shape functions (1.3)-(1.4) we can write the value of potential energyfor the second finite element as:·Z x21dNdx T · dNΠe a{u}{u}dx2dx (x)Z x210TTT {u} [N ] bdx {u}Rx1T(1.14)The condition for the minimum of Π is:δΠ Π Πδu1 . δun 0 u1 un(1.15)which is equivalent to Π 0, uii 1.n(1.16)It is easy to check that differentiation of Π in respect to ui gives the finite element equilibrium equationwhich is coincide with equation obtained by the Galerkin method: Z x2 ·dN Tx1dx· Z x2dNdx{u} EAdxx1([N ]T bdx 0R) 0(1.17)Example. Obtain shape functions for the one-dimensional quadratic element with three nodes. Uselocal coordinate system 1 ξ 1.1-12031xSolution. With shape functions, any field inside element is presented as:u(ξ) XNi ui , i 1, 2, 3At nodes the approximated function should be equal to its nodal value:u( 1) u1u(0) u2u(1) u3

CHAPTER 1. INTRODUCTION10Since the element has three nodes the shape functions can be quadratic polynomials (with three coefficients). The shape function N1 can be written as:N1 α1 α2 ξ α3 ξ 2Unknown coefficients αi are defined from the following system of equations:N1 ( 1) α1 α2 α3 1N1 (0) α1 0N1 (1) α1 α2 α3 0The solution is: α1 0, α2 1/2, α3 1/2. Thus the shape function N1 is equal to:1N1 ξ(1 ξ)2Similarly it is possible to obtain that the shape functions N2 and N3 are equal to:N2 1 ξ 21N3 ξ(1 ξ)2

Chapter 2Finite Element Equations for HeatTransfer2.1Problem StatementLet us consider an isotropic body with temperature dependen