 (Parte 1 de 3)

Professor Peter Hunter p.hunter@auckland.ac.nz

Associate Professor Andrew Pullan a.pullan@auckland.ac.nz

Department of Engineering Science

The University of Auckland New Zealand c Copyright 1997-2003

Department of Engineering Science The University of Auckland

Contents

 1.1 Representing a One-Dimensional Field 1 1.2 Linear Basis Functions 2 1.3 Basis Functions as Weighting Functions 4 1.4 Quadratic Basis Functions 7 1.5 Two- and Three-Dimensional Elements 7 1.6 Higher Order Continuity 10 1.7 Triangular Elements 14 1.8 Curvilinear Coordinate Systems 16 1.9 CMISS Examples 20

1 Finite Element Basis Functions 1

 2.1 One-Dimensional Steady-State Heat Conduction 23 2.1.1 Integral equation 24 2.1.2 Integration by parts 24 2.1.3 Finite element approximation 25 2.1.4 Element integrals 26 2.1.5 Assembly 27 2.1.6 Boundary conditions 29 2.1.7 Solution 29 2.1.8 Fluxes 29 2.2 An -Dependent Source Term 30 2.3 The Galerkin Weight Function Revisited 31 2.4 Two and Three-Dimensional Steady-State Heat Conduction 32 2.5 Basis Functions - Element Discretisation 34 2.6 Integration 36 2.7 Assemble Global Equations 37 2.8 Gaussian Quadrature 39 2.9 CMISS Examples 42

 3.1 Introduction 43 3.2 The Dirac-Delta Function and Fundamental Solutions 43 3.2.1 Dirac-Delta function 43
 3.3 The Two-Dimensional Boundary Element Method 48 3.4 Numerical Solution Procedures for the Boundary Integral Equation 53 3.5 Numerical Evaluation of Coefficient Integrals 5 3.6 The Three-Dimensional Boundary Element Method 57 3.7 A Comparison of the FE and BE Methods 58 3.8 More on Numerical Integration 60 3.8.1 Logarithmic quadrature and other special schemes 60 3.8.2 Special solutions 61 3.9 The Boundary Element Method Applied to other Elliptic PDEs 61 3.10 Solution of Matrix Equations 61 3.1 Coupling the FE and BE techniques 62 3.12 Other BEM techniques 64 3.12.1 Trefftz method 64 3.12.2 Regular BEM 64 3.13 Symmetry 65 3.14 Axisymmetric Problems 67 3.15 Infinite Regions 69 3.16 Appendix: Common Fundamental Solutions 72 3.16.1 Two-Dimensional equations 72 3.16.2 Three-Dimensional equations 72 3.16.3 Axisymmetric problems 73 3.17 CMISS Examples 73

i CONTENTS

 4.1 Introduction 75 4.2 Truss Elements 76 4.3 Beam Elements 79 4.4 Plane Stress Elements 81 4.4.1 Notes on calculating nodal loads 83 4.5 Three-Dimensional Elasticity 84 4.5.1 Weighted Residual Integral Equation 85 4.5.2 The Principle of Virtual Work 86 4.5.3 The Finite Element Approximation 87 4.6 Linear Elasticity with Boundary Elements 89 4.7 Fundamental Solutions 91 4.8 Boundary Integral Equation 93 4.9 Body Forces (and Domain Integrals in General) 96 4.10 CMISS Examples 98

4 Linear Elasticity 75

 5.1 Introduction 9 5.2 Finite Differences 9 5.2.1 Explicit Transient Finite Differences 9 5.2.2 Von Neumann Stability Analysis 101 5.3 The Transient Advection-Diffusion Equation 103 5.4 Mass lumping 106 5.5 CMISS Examples 108

CONTENTS i

 6.1 Introduction 1 6.2 Free Vibration Modes 1 6.3 An Analytic Example 113 6.4 Proportional Damping 114 6.5 CMISS Examples 115

6 Modal Analysis 1

 7.1 Achieving a Boundary Integral Formulation 117 7.2 Removing Domain Integrals due to Inhomogeneous Terms 118 7.2.1 The Galerkin Vector technique 118 7.2.2 The Monte Carlo method 119 7.2.3 Complementary Function-Particular Integral method 120 7.3 Domain Integrals Involving the Dependent Variable 120 7.3.1 The Perturbation Boundary Element Method 121 7.3.2 The Multiple Reciprocity Method 122 7.3.3 The Dual Reciprocity Boundary Element Method 124

7 Domain Integrals in the BEM 117

 8.1 Time-Stepping Methods 135 8.1.1 Coupled Finite Difference - Boundary Element Method 135 8.1.2 Direct Time-Integration Method 137 8.2 Laplace Transform Method 138 8.3 The DR-BEM For Transient Problems 139 8.4 The MRM for Transient Problems 140

8 The BEM for Parabolic PDES 135 Bibliography 143 Index 147

Chapter 1 Finite Element Basis Functions

1.1 Representing a One-Dimensional Field

Consider the problem of finding a mathematical expression to represent a one-dimensional field e.g., measurements of temperature against distance along a bar, as shown in Figure 1.1a. (b)

FIGURE 1.1: (a) Temperature distribution along a bar. The points are the measured temperatures. (b) A least-squares polynomial fit to the data, showing the unacceptable oscillation between data points.

One approach would be to use a polynomial expression and to estimate the values of the parameters , , and from a least-squares fit to the data. As the degree of the polynomial is increased the data points are fitted with increasing accuracy and polynomials provide a very convenient form of expression because they can be differentiated and integrated readily. For low degree polynomials this is a satisfactory approach, but if the polynomial order is increased further to improve the accuracy of fit a problem arises: the polynomial can be made to fit the data accurately, but it oscillates unacceptably between the data points, as shown in Figure 1.1b.

To circumvent this, while retaining the advantages of low degree polynomials, we divide the bar into three subregions and use low order polynomials over each subregion - called elements. For later generality we also introduce a parameter which is a measure of distance along the bar. is plotted as a function of this arclength in Figure 1.2a. Figure 1.2b shows three linear polynomials in fitted by least-squares separately to the data in each element.

2 FINITE ELEMENT BASIS FUNCTIONS

(b) (a)

FIGURE 1.2: (a) Temperature measurements replotted against arclength parameter . (b) The domain is divided into three subdomains, elements, and linear polynomials are independently fitted to the data in each subdomain.

1.2 Linear Basis Functions A new problem has now arisen in Figure 1.2b: the piecewise linear polynomials are not continuous in across the boundaries between elements. One solution would be to constrain the parameters ,, etc. to ensure continuity of across the element boundaries, but a better solution is to replace the parameters and in the first element with parameters and , which are the values of at the two ends of that element. We then define a linear variation between these two values by where is a normalized measure of distance along the curve.

We define such that and refer to these expressions as the basis functions associated with the nodal parameters and

. The basis functions and are straight lines varying between and as shown in

Figure 1.3.

It is convenient always to associate the nodal quantity with element node and to map the temperature defined at global node onto local node of element by using a connectivity matrix i.e., where = global node number of local node of element . This has the advantage that the

1.2 LINEAR BASIS FUNCTIONS 3 1

holds for any element provided that and are correctly identified with their global counterparts, as shown in Figure 1.4. Thus, in the first element node node element element element node node

FIGURE 1.4: The relationship between global nodes and element nodes.

with and .

In the second element is interpolated by with and , since the parameter is shared between the first and second elements

4 FINITE ELEMENT BASIS FUNCTIONS the temperature field is implicitly continuous. Similarly, in the third element is interpolated by with and , with the parameter being shared between the second and third elements. Figure 1.6 shows the temperature field defined by the three interpolations (1.1)–(1.3).

node node

node + node + element elementelement

FIGURE 1.5: Temperature measurements fitted with nodal parameters and linear basis functions. The fitted temperature field is now continuous across element boundaries.

1.3 Basis Functions as Weighting Functions

It is useful to think of the basis functions as weighting functions on the nodal parameters. Thus, in element 1 at which is the value of at the left hand end of the element and has no dependence on which depends on and , but is weighted more towards than at which depends equally on and

1.3 BASIS FUNCTIONS AS WEIGHTING FUNCTIONS 5 which depends on and but is weighted more towards than which is the value of at the right hand end of the region and has no dependence on .

Moreover, these weighting functions can be considered as global functions, as shown in Fig- ure 1.6, where the weighting function associated with global node is constructed from the basis functions in the elements adjacent to that node.

(b) (c)

(d)

FIGURE 1.6: (a) (d) The weighting functions associated with the global nodes , respectively. Notice the linear fall off in the elements adjacent to a node. Outside the immediately adjacent elements, the weighting functions are defined to be zero.

For example, weights the global parameter and the influence of falls off linearly in the elements on either side of node 2.

We now have a continuous piecewise parametric description of the temperature field but in order to define we need to define the relationship between and for each element. A convenient way to do this is to define as an interpolation of the nodal values of . For example, in element 1

(1.4) and similarly for the other two elements. The dependence of temperature on , , is therefore

6 FINITE ELEMENT BASIS FUNCTIONS defined by the parametric expressions where summation is taken over all element nodes (in this case only ) and the parameter (the “element coordinate”) links temperature to physical position . provides the mapping between the mathematical space and the physical space , as illustrated in

Figure 1.7. at

FIGURE 1.7: Illustrating how and are related through the normalized element coordinate .

The values of and are obtained from a linear interpolation of the nodal variables and then plotted as . The points at are emphasized.

The essential property of the basis functions defined above is that the basis function associated with a particular node takes the value of when evaluated at that node and is zero at every other node in the element (only one other in the case of linear basis functions). This ensures the linear independence of the basis functions. It is also the key to establishing the form of the basis functions for higher order interpolation. For example, a quadratic variation of over an element requires three nodal parameters , and

The quadratic basis functions are shown, with their mathematical expressions, in Figure 1.8. Notice that since must be zero at (node ), must have a factor and since it is also zero at (node ), another factor is

. Finally, since is at (node ) we have . Similarly for the other two basis functions.

(c) (a) (b)

FIGURE 1.8: One-dimensional quadratic basis functions.

1.5 Two- and Three-Dimensional Elements

Two-dimensional bilinear basis functions are constructed from the products of the above onedimensional linear functions as follows

8 FINITE ELEMENT BASIS FUNCTIONS Let where (1.6)

Note that = where and are the one-dimensional linear basis functions. Similarly, = etc. These four bilinear basis functions are illustrated in Figure 1.9.

node node node node

FIGURE 1.9: Two-dimensional bilinear basis functions.

Notice that is at node and zero at the other three nodes. This ensures that the temperature receives a contribution from each nodal parameter weighted by and that when is evaluated at node it takes on the value . As before the geometry of the element is defined in terms of the node positions ,

1.5 TWO- AND THREE-DIMENSIONAL ELEMENTS 9 by which provide the mapping between the mathematical space (where ) and the physical space .

Higher order 2D basis functions can be similarly constructed from products of the appropriate 1D basis functions. For example, a six-noded (see Figure 1.10) quadratic-linear element (quadratic in and linear in ) would have where

FIGURE 1.10: A -node quadratic-linear element (node numbers circled).

Three-dimensional basis functions are formed similarly, e.g., a trilinear element basis has eight nodes (see Figure 1.1) with basis functions

10 FINITE ELEMENT BASIS FUNCTIONS

1.6 Higher Order Continuity

All the basis functions mentioned so far are Lagrange1 basis functions and provide continuity of across element boundaries but not higher order continuity. Sometimes it is desirable to use basis functions which also preserve continuity of the derivative of with respect to across element boundaries. A convenient way to achieve this is by defining two additional nodal parameters

. The basis functions are chosen to ensure that and and since is shared between adjacent elements derivative continuity is ensured. Since the number of element parameters is 4 the basis functions must be cubic in . To derive these cubic Hermite2 basis functions let

1Joseph-Louis Lagrange (1736-1813). 2Charles Hermite (1822-1901).

1.6 HIGHER ORDER CONTINUITY 1 and impose the constraints

These four equations in the four unknowns , , and are solved to give

Substituting , , and back into the original cubic then gives or, rearranging, (1.14) where the four cubic Hermite basis functions are drawn in Figure 1.12. One further step is required to make cubic Hermite basis functions useful in practice. The derivative defined at node is dependent upon the element -coordinate in the two ad- jacent elements. It is much more useful to define a global node derivative where is arclength and then use where is an element scale factor which scales the arclength derivative of global node to the -coordinate derivative of element node

. Thus is constrained to be continuous across element boundaries rather than . A two- dimensional bicubic Hermite basis requires four derivatives per node and

12 FINITE ELEMENT BASIS FUNCTIONS slope slope

The need for the second-order cross-derivative term can be explained as follows; If is cubic in and cubic in , then is quadratic in and cubic in , and is cubic in and quadratic

. Now consider the side 1–3 in Figure 1.13. The cubic variation of with is specified by the four nodal parameters , , and

. But since (the normal derivative) is also cubic in along that side and is entirely independent of these four parameters, four additional parameters are required to specify this cubic. Two of these are specified by and , and the remaining two by and .

1.6 HIGHER ORDER CONTINUITY 13 nodenode node node

FIGURE 1.13: Interpolation of nodal derivative along side 1–3.

The bicubic interpolation of these nodal parameters is given by

14 FINITE ELEMENT BASIS FUNCTIONS where are the one-dimensional cubic Hermite basis functions (see Figure 1.12).

As in the one-dimensional case above, to preserve derivativecontinuity in physical x-coordinate space as well as in -coordinate space the global node derivatives need to be specified with respect to physical arclength. There are now two arclengths to consider: , measuring arclength along the-coordinate, and , measuring arclength along the -coordinate. Thus where and are element scale factors which scale the arclength derivatives of global node to the -coordinate derivatives of element node .

The bicubic Hermite basis is a powerful shape descriptor for curvilinear surfaces. Figure 1.14 shows a four element bicubic Hermite surface in 3D space where each node has the following twelve parameters and

1.7 Triangular Elements

Triangular elements cannot use the and coordinates defined above for tensor product elements

(i.e., two- and three- dimensional elements whose basis functions are formed as the product of onedimensional basis functions). The natural coordinates for triangles are based on area ratios and are called Area Coordinates . Consider the ratio of the area formed from the points , and in Figure 1.15 to the total area of the triangle

Area Area

1.7 TRIANGULAR ELEMENTS 15 12 parameters per node

FIGURE 1.14: A surface formed by four bicubic Hermite elements.

P( , ) Area

FIGURE 1.15: Area coordinates for a triangular element.

16 FINITE ELEMENT BASIS FUNCTIONS where is the area of the triangle with vertices , and

Notice that is linear in and . Similarly, area coordinates for the other two triangles containing and two of the element vertices are

Area Area

Area Area where and .

Notice that .

Area coordinate varies linearly from when lies at node or to when lies at node and can therefore be used directly as the basis function for node for a three node triangle. Thus, interpolation over the triangle is given by where , and .

Six node quadratic triangular elements are constructed as shown in Figure 1.16.

FIGURE 1.16: Basis functions for a six node quadratic triangular element.

1.8 Curvilinear Coordinate Systems

It is sometimes convenient to model the geometry of the region (over which a finite element solution is sought) using an orthogonal curvilinear coordinate system. A 2D circular annulus, for ex-

1.8 CURVILINEAR COORDINATE SYSTEMS 17 ample, can be modelled geometrically using one element with cylindrical polar -coordinates, e.g., the annular plate in Figure 1.17a has two global nodes, the first with and the second with .

(b) (c)(a)

FIGURE 1.17: Defining a circular annulus with one cylindrical polar element. Notice that element vertices and in -space or -space, as shown in (b) and (c), respectively, map onto the single global node in -space in (a). Similarly, element vertices and map onto global node .

Global nodes and , shown in -space in Figure 1.17a, each map to two element vertices in -space, as shown in Figure 1.17b, and in -space, as shown in Figure 1.17c. Thecoordinates at any point are given by a bilinear interpolation of the nodal coordinatesand as where the basis functions are given by (1.6).

Three orthogonal curvilinear coordinate systems are defined here for use in later sections.

Cylindrical polar :

Spherical polar :

18 FINITE ELEMENT BASIS FUNCTIONS

Prolate spheroidal :

y r

The prolate spheroidal coordinates rae illustrated in Figure 1.18 and a single prolate spheroidal element is shown in Figure 1.19. The coordinates are all trilinear in . Only four global nodes are required provided the four global nodes map to eight element nodes as shown in Figure 1.19.

1.8 CURVILINEAR COORDINATE SYSTEMS 19

(d)(c)

(a) (b)

FIGURE 1.19: A single prolate spheroidal element, shown (a) in -coordinates, (c) in-coordinates and (d) in -coordinates, (b) shows the orientation of the

(Parte 1 de 3)