University of Newcastle Upon Tyne s1

Newcastle University

School of Mechanical & Systems Engineering

An Introduction to

THE FINITE ELEMENT METHOD

The Finite Element method for 2D continuum problems

(c) John C. Appleby, 2008, School of Mechanical & Systems Engineering,

Newcastle University, Newcastle upon Tyne NE1 7RU, U.K.

Tel (+44) (0) 191-222 6286 Fax (+44) (0) 191-222 8600

E-mail: Web: http://www.staff.ncl.ac.uk/john.appleby

Contents

23

1. Introduction

The Finite Element Method is a computational method for the analysis of stress, vibration, heat conduction, fluid flow, electrostatics and acoustics; it is suitable for a wide range of 2D and 3D continuum problems, which may be steady or unsteady, linear or non-linear. The method has been developed over about thirty years, starting with traditional matrix methods of structural analysis for frameworks consisting of rod and beam elements. In fact, problems with rods and beams are now also handled in finite element terms, since so many aspects are the same, but some features are necessarily distinct.

The distinguishing feature of the finite element method for continuum problems is that we choose to represent the solution of the problem in a particular way at the outset. We in effect replace the original problem domain and material, which behaves in a smooth, continuous but probably complex way, with a domain, geometry and 'material' which have simpler behaviour, although not as smooth. We then apply the normal equations for the problem (eg equations of plane stress, Laplace's equation etc) to this simplified form of solution. Because our anticipated solution is not smooth, we need the equations in integral form, and these must be derived from the more familiar differential equations for many problems. The equations may themselves involve approximation, such as plane stress or plate bending formulations, or be exact, as for heat conduction.

An example

If we wished to solve the problem of a freely-hanging cable (the catenary), with true solution smooth and continuous, we might choose to represent it as a chain of finite links. For the catenary problem, this is not in fact a practical method of solution, but it is a useful illustration of some of the characteristic features of the finite element method.

Since each link is defined completely by its end points, what was a continuum problem of infinitely many degrees of freedom has been represented by a system (mesh) of 8 elements joined at 9 nodes, each of which has an x-displacement and a y-displacement, giving 18 degrees of freedom.

If we went on to formulate the whole problem in this way, we would obtain a system of 18 simultaneous linear equations in these 18 variables. Now, 4 variables are known, namely zero displacements at the end points, but there are 4 corresponding unknown reactions, so we do indeed have 18 unknowns and a well-defined problem.

After solving the equations for the unknown displacements and forces, we may also solve for derived quantities, such as stresses. Moreover, we can calculate the displacements etc at any point in the structure, since each link is defined entirely in terms of its end-points.

For framework problems, such as pylons, bridges, conveyors etc made from girders, the structure is already in discrete form, and the only approximation is in the choice of equations such as beams with or without shear effects.

This example, though artificial, illustrates the principal stages of a static, linear finite element analysis for a 2D continuum problem:

· Discretise the structure (design a mesh of nodes and elements), and choose a form of solution in each element.

· Apply the equations to obtain a system of simultaneous linear equations.

· Apply boundary conditions and forces (gravity in this example).

· Solve the equations for the primary variables at the nodes.

· Solve for reactions and any derived quantities.

· Calculate values within each element in terms of the nodal values if desired.

A 3D problem follows exactly the same steps. A framework problem is different only in the first step. A non-linear problem requires iteration to solve the main system of equations. An unsteady problem solves the equations at each time step. Thus a vast range of problems are susceptible to essentially the same approach if suitable equations can be derived, and it is possible to construct general purpose finite element packages which use common modules for many different problem types. In these notes we concentrate on continuum problems, in contrast to treatments designed primarily for structural engineers; rod and beam elements are considered in section 7.

Principal advantages of the finite element method

The finite element method was developed initially for structural engineering problems. Because the mesh is designed by the user (though sometimes with computer assistance), it is possible to conform to features of the structure such as:

· fitting the mesh to the geometry, which can involve the use of curved-sided elements,

· having different materials and eg thicknesses within different elements (since the integral equations do not fail at discontinuities),

· elements of varying sizes in different regions to reflect the complexity of the solution in each part.

Disadvantages of finite element analysis are that it is more difficult to understand and to program than eg the finite difference method, and can be more difficult to apply for some kinds of domain and boundary conditions, such as free-surface fluid flows.

Comparison with other methods

The finite element method applies the exact equations (or the normal approximation) to a chosen approximated form of the solution. In contrast, the finite difference method approximates the differential equations for the problem, applied to the solution at points on a grid. For example, .

A third method, the boundary element method, also uses an integral formulation, but discretises the solution only on the boundary, which is more economical, particularly when it is the boundary that is actually of greatest interest.

The first two methods produce large matrices of coefficients that are usually banded, producing economies of storage and solution time; the third method produces a smaller matrix, but more densely populated.

2. Some examples of 2D problems and their chosen solution form

(i) A plane stress problem

A simple triangular bracket, attached to a wall, is loaded at its tip. The bracket is made from thin sheet mild steel, and the load is in-plane, so the equations of plane stress apply. The simplest model (or representation) of this structure is a single element which deforms in a linear way, ie the displacements are linear functions of x and y, and the deformed plate remains a straight-sided triangle.

With this assumption about the way the structure behaves, we know at the outset that the deformed shape must be something like that shown:

Note that we have not so far done any calculation, but simply imposed on the problem a chosen form of solution. Normally, we would choose a form that we expect to be capable of representing the true solution to reasonable accuracy. In this case, the bracket should really bend much more near the tip than we have allowed, and it is not surprising if our model is much stiffer than it should be.

Exercise

To obtain an order of magnitude estimate of the solution, model the domain as a rectangle of comparable dimensions and apply the force as a one-dimensional analogue. (E - Young's modulus for mild steel is )

We can improve things by using four elements, each of which behaves linearly, but which together permit more complex behaviour:

In practice, we often use more than one mesh, of progressively increasing refinement, to help us assess the likely accuracy of the solution. If the increased resolution makes little difference, it suggests that we may have reached convergence to the true solution. It is important to realise that convergence may be very slow, so that we are much further from the true values than it appears, and convergence to the wrong solution is also possible! We must have independent verification, if only from experience or an approximate estimate, before we can rely on any results.

Overall, these meshes produce 6 and 12 variables respectively (3 or 6 nodes x 2 displacements). 3 variables are known (fixed, zero displacements) and also 3 or 9 forces (one imposed, the remainder zero). Therefore we solve for 3 displacements and 3 reactions or 9 displacements and 3 reactions in the two cases.

If the boundary conditions are changed, the matrix of coefficients remains unchanged, and it is economical to re-use this with new boundary conditions, especially as the matrix may be processed independently of the altered information, and it is this stage that takes most computer time. All finite element packages should offer such a re-solution facility.

(ii) A heat conduction problem

A square region (1m along each edge) has the temperature fixed along the edges, and we wish to find the temperature within the region, as well as the heat flow. The thickness and material affect the size of the fluxes, but not the temperature distribution, provided the material is homogeneous. If we model the solution as four elements, joined at five nodes, and we assume the temperature varies linearly within each triangular element, we are in effect assuming that the temperature contours are straight lines.

Exercise

Sketch the expected form of the contours, and give a rough estimate of the central value of temperature.

Since it is obvious that the central value will be below , we can sketch the expected solution as shown:

Once again, the form of solution is clearly not satisfactory, and we can design a better mesh which will permit a smoother solution, but even such a simple mesh as this, or the first one for the bracket example, may be useful to:

· check the operation of the program, data input, boundary conditions etc

· check that we have a well-defined problem

· give an order of magnitude estimate of the solution

· identify areas where the mesh should be refined

This problem has one degree of freedom at any point, namely the temperature, and so there are 5 variables (the nodal temperatures) overall. We obtain a system of 5 equations. With the boundary conditions indicated, we know 4 of these temperatures, and our 5 unknowns are the central nodal temperature, together with 4 values of heat flux, with the results given at the nodes. These must be interpreted as integrated fluxes into or out of the boundary surrounding each node. For this mesh, that is the region extending from mid-way between the nodes on each side. Here the total flux will be zero, with flow into the domain around one node, and out of it around three.

If we change the boundary conditions, eg imposing fluxes rather than temperatures at the boundary, we can use the same set of equations, as described above.

Also note that Laplace's equation, or Poisson's equation if source terms are included, applies to a range of other problems, including electrostatics, torsion and groundwater and other potential flows.

(iii) Curved-sided elements

One of the principal advantages of the finite element method is the ability to match closely the true geometry of the problem. Simple linear elements, like the 3 node triangle, always have straight sides, but elements whose behaviour is described by quadratic or higher functions can also have a geometry described by curves.

Here a single element is able to model a curved side, though in practice many more elements would be used. In general, one quadratic or higher order element of this type has superior performance to the eight linear triangles that could replace it with about the same number of degrees of freedom.

(iv) Mesh grading

As well as modelling the geometry of the problem, element sizes and shapes may be chosen to give uniform accuracy to the solution, by using more nodes and elements in areas of greatest 'activity'; where changes are most rapid and complex, smaller elements should be used.

Here, quadrilateral elements (which might be 4-node, 8-node or more) are graded in size to suit a problem of a point load on a block. Note that 'quadrilateral' does not imply 'rectangular'. In practice, the natural symmetry of this problem would be used to reduce the domain and corresponding calculation by half.

3. Solving the problems

The above examples have indicated how we can model the geometry and the anticipated solution behaviour. Next we shall look at two of them to see what the computer program does to reach a solution, but deferring until later the derivation of the finite element equations.

(i) The simple bracket

The bracket problem above, when modelled with a single triangular element, requires the following data (if the package FINEL is to be used) or very similar data. Some packages offer mouse-driven interactive mesh design, so that the geometry is 'drawn' on to the screen, but the data created will be of the same type.

title

bracket in plane stress optional title

coordinates

1 0.0 0.0 node number, x coordinate, y coordinate

2 0.2 0.0

3 0.0 0.2

elements

1 1 2 3 element number, nodes numbered anticlockwise

dimensions

1 0.01 size code 1 will have thickness 0.01m

constraints

1 1 0.0 node 1, d.o.f. 1, constrained to value 0.0