Lecture 4

## Introduction to the Finite Element Method

Objective: To introduce the concept of using displacement interpolation to obtain approximate stiffness matrices for 2 or 3 dimensional elements.  Even modelling a beam with a stiffness matrix is approximate if distributed or intermediate loads are present.

The stiffness matrix of a beam can be found exactly (for some bending theory). To generalize to a 2 or 3 dimensional piece of a structure, however, requires approximation to describe the boundary conditions on the element.

eg. A piece of plate deforming inplane:

Given displacements u1 and u2, how will displacements vary in between? We could assume a linear variation. What forces it takes to cause the corner displacements will depend on what is assumed about how displacement varies. In fact, by assuming a displacement variation along an edge we implicitly distribute loads along the edge to make it happen.

In practice, we typically interpolate displacements between nodal values over the whole area or volume of an element, so that deformation (and hence strain) anywhere is described by nodal values.  We then need to integrate over the element to estimate its stiffness from its strain energy.

### Meshing

If enough small elements are used, a two or three-dimensional model of an object to be analysed can be built, that will successfully approximate the deformation and stressing of the real object.

In practice, placing a mesh on an object of arbitrary shape, requires solid modeling CAD facilities.  An arbitrary area can be meshed automatically with triangular elements, or an arbitrary volume with tetrahedra.

Meshing can also be done by operations that copy a pattern of elements.  For instance "extrusion" can be used to make a solid model from a planar model.

CAD geometry may not be what is needed for meshing, unless a solid model is needed.  There is often a need to create more suitable geometry, which is why finite element software often allows geometry to be created in a preprocessor.  The FEMAP program used with NASTRAN allows regular geometry to be created, and also can estimate the mid-surface of a solid, so plate elements can be put there.  In STRAND7, the geometry must presently be either imported or captured by the mesh itself.

### The Constant Strain Triangle

This is the simplest finite element possible, as constant strain in the simplest state of strain possible that can approximate the behavior of a small region of plate deforming inplane. Displacement anywhere in the element is expressed in terms of x and y displacements at 3 corner nodes. That is, it depends on 6 nodal values. Writing x displacement as u and y displacement as v:

u = A + Bx + Cy                    v = D + Ex + Fy

The six parameters A to F need to be expressed in terms of the nodal displacements shown below. Note this is an arbitrary linear variation of displacement in the x and y directions.

### Strain-Displacement Relations

Direct strain is change in length/length = change in displacement over a short length. Shear strain is change in a right angle. This change is due to the rotation of a line initially in the x direction , which for small angles is Dv/Dx and the rotation of a line initially in the y direction, which for small angles is Du/Dy, as shown below.

Strains are derivatives of displacements.

The definitions of strains in the xy plane, in terms of x and y displacements u and v are

Using u = A + Bx + Cy and v = D + Ex + Fy the strains are

ex = B       ey = F   and    gxy = E + C

i.e. the strains are constants for given nodal displacement values. For elements other than the constant strain triangle, the strains vary over the element.

Why we need expressions for Strain

(a)In order to find a stiffness matrix to represent the element.  If strain is known as f(nodal displacements) then we can find strain energy per volume as f(nodal displacements). Integrating over an element its total strain energy can be found. Strain energy expressed in terms of nodal displacements enables stiffness terms to be found

(b)Once the full set of equations for a structural model has been solved and the displacements at the nodes are known, then we wish to find the stresses in the element. This is done by first finding the strains from the nodal displacements.

### Types of Finite Elements

For 2D problems small triangular or quadrilateral regions are commonly used, with linear or quadratic displacement variations assumed. The quadratic elements have mid-side nodes (a quadratically curving displacement variation can be fitted through 3 nodal values along a side).

In 3D, a small region of a solid is described by tetrahedra, pentahedra (wedges) or hexahedra (bricks).

### h versus p FINITE ELEMENT METHODS

In a typical finite element program, the user selects to use linear of quadratic elements. To obtain greater accuracy, in any region where stresses are changing rapidly, a finer mesh of smaller elements is used, with the same order of polynomials used to express the displacement variation in each element. This method of improving accuracy is called the h finite element method. (Think of h as the side length of an element, which is being reduced.)

An alternative, which is gaining popularity, is the p finite element method, in which the element size is kept the same, but higher order polynomials are used to represent the displacement variation in each element, in order to improve the accuracy. (p = polynomial order)   To enable this to occur, the elements must define the shape (eg the curvature of a surface) accurately.  This requires more data.  As p is increased, element stiffness matrices become larger, and the computation to solve the equations increases, as more equations are written in terms of extra unknowns.

This latter approach to reanalysis is easier to automate, as remeshing is avoided.  The error in a solution can be estimated (for instance from the size of the jumps in stress at element boundaries), and the solution rerun with higher p values until it converges.

The stresses calculated by a finite element program disagree at element boundaries. (eg for a constant strain triangle CST, they are constant in each element). An estimate of their accuracy is needed in order to decide if a problem should be reanalyzed. The size of the stress discontinuity is one indication of accuracy.  The most accurate stresses are those calculated at element centers, away from the discontinuities at their edges. Programs often average stresses at the nodes between adjacent elements. This may not improve the accuracy, as for instance, on a free surface, both elements meeting at a node may be overestimating.  Options exist in both STRAND7 and NASTRAN to control whether the program does this averaging of stresses at nodes.

Accuracy can be improved by iteratively reanalyzing a problem, with either a finer mesh of elements  (h method) or by increasing the order of polynomial functions used to represent the displacement variation (p method).

### Obtaining An Approximate Element Stiffness matrix

The stiffness coefficients relating forces and displacements at all nodes of an element can be found by estimating its strain energy. This is analogous to the procedure used in dynamics with lumped parameter systems. Equations of motion are written by differentiating a work-energy expression (Lagranges’ equations). Here each nodal displacement serves as a generalized coordinate. In a statics problem, there is no kinetic energy, so if potential energy of one element = V then

V is usually just the elastic strain energy of the element.

ui = one of the terms in the vector of displacements at nodes of the element.

Fi is the force on the element at some node corresponding to ui.

By expressing V in terms of nodal displacements, stiffness terms are obtained.

### Example – Stiffness Matrix Of a Spring

For a simple spring, with displacement u1 at one end and u2 at the other,

potential energy,   V = ½ k (u2 – u1)2
½ k u22 + k u1u2 + 1/2 k u12

Hence the element stiffness matrix is

### General Strain Energy Expression

Any elastic structure resembles a network of springs. If strain energy of one finite element is written in terms of its nodal displacements, two types of terms arise (found by integrating over the element).

Terms involving one node, say node i, ie.  ½ kii ui2

Terms involving an interaction between two nodes i and j ie.  kijuiuj

Hence the force on node i is

These give row i of the element stiffness matrix (kii is the diagonal term). In matrix notation, the strain energy of an element is ½ ueT[Ke] u where ue is a vector of nodal displacements and [Ke] is the element stiffness matrix. Note the analogy to ½ k u2 for a spring. This expression differentiates to give [Ke]ue.

### Summary – Overview of the Finite Element Method

A finite element model divides a 2 or 3 dimensional region into small, but finite elements. In the usual (h) method, may small elements are used, in each of which simple polynomial functions are used to interpolate displacements (or some other field variable) between values at nodes. This enables strains and hence strain energy to be found in terms of nodal displacements. An approximate stiffness matrix can be found for each element, linking all forces and displacements at its nodes, by differentiating an expression for the strain energy of the element. The displacements can therefore be calculated, as can stresses within the element, once the nodal displacements are known. The stresses and strains will be approximate and discontinuous at element boundaries.  Displacements are approximate but continuous at element boundaries.

### Generalization to other Problems

While the finite element method was developed to solve structural problems, it has since been generalized to solve many other problems. One example is the computation of electric and magnetic fields. Another common example is the computation of temperature distributions in solids.

The field variable (pressure, velocity potential, temperature or whatever) is interpolated over each element between nodal values. Typically, integration over the element enables an energy balance to be established approximately and equations are found by differentiating work or energy expressions with respect to each nodal variable.

Different problem types may be combined (eg thermal stress problems involving first finding the temperature distribution).

### Interpolation Functions

Instead of writing displacement in terms of unknown constants, and then relating those constants to nodal displacements, we can write displacement directly in terms of nodal displacements.

ie u = N1(x,y) u1 + N2(x,y)u2 + N3(x,y) u3

The N functions are called interpolation functions. Each such function is equal to 1 at one node and zero at the others, so that at node 1 u = u1 etc.

### Weighted Residual Viewpoint

This is a more general view of the approximation being made in a finite element model.

Consider a partial differential equation.
Eg. A(T,x,y,z) = 0 where A is some operator (eg del-squared). T is temperature (say). We wish to make this equation true approximately over the volume of one element, in a weighted average sense.

ie. Integrate:  over an element, so that

T = SNiTi, Ti being the temperature of node I, and Ni an interpolation function. If we make weighting factor w =Nj (the interpolation function for node j) then we have the Galerkin method. The finite element method is a type of Galerkin method. That is, by choosing w to be each interpolation function in turn, we force the equation being solved to be approximately true in the vicinity of each node, where that interpolation function is approximately equal to one, each time generating an approximate equation, by integrating over the element.