eg. A piece of plate deforming inplane:
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.
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.
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
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.
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
= 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
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
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
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
Evaluating Accuracy of Answers
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. ½
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] ue
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
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
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.