The Finite Element Method: A Four-Article Series - Part 3
The following four-article series was published in a newsletter
of the American Society of Mechanical Engineers (ASME).
It serves as an introduction to the recent analysis discipline
known as the finite element method. The author is
an engineering consultant and expert witness specializing in
finite element analysis.
FINITE ELEMENT ANALYSIS: Solution
by Steve Roensch, President, Roensch & Associates
Third in a four-part series
While the pre-processing and post-processing phases of the
finite element method are interactive and time-consuming for the
analyst, the solution is often a batch process, and is demanding
of computer resource. The governing equations are assembled into
matrix form and are solved numerically. The assembly process
depends not only on the type of analysis (e.g. static or
dynamic), but also on the model's element types and properties,
material properties and boundary conditions.
In the case of a linear static structural analysis, the
assembled equation is of the form Kd = r, where K
is the system stiffness matrix, d is the nodal degree of
freedom (dof) displacement vector, and r is the applied
nodal load vector. To appreciate this equation, one must begin
with the underlying elasticity theory. The strain-displacement
relation may be introduced into the stress-strain relation to
express stress in terms of displacement. Under the assumption of
compatibility, the differential equations of equilibrium in
concert with the boundary conditions then determine a unique
displacement field solution, which in turn determines the strain
and stress fields. The chances of directly solving these
equations are slim to none for anything but the most trivial
geometries, hence the need for approximate numerical techniques
presents itself.
A finite element mesh is actually a displacement-nodal
displacement relation, which, through the element interpolation
scheme, determines the displacement anywhere in an element given
the values of its nodal dof. Introducing this relation into the
strain-displacement relation, we may express strain in terms of
the nodal displacement, element interpolation scheme and
differential operator matrix. Recalling that the expression for
the potential energy of an elastic body includes an integral for
strain energy stored (dependent upon the strain field) and
integrals for work done by external forces (dependent upon the
displacement field), we can therefore express system potential
energy in terms of nodal displacement.
Applying the principle of minimum potential energy, we may set
the partial derivative of potential energy with respect to the
nodal dof vector to zero, resulting in: a summation of element
stiffness integrals, multiplied by the nodal displacement
vector, equals a summation of load integrals. Each stiffness
integral results in an element stiffness matrix, which sum to
produce the system stiffness matrix, and the summation of load
integrals yields the applied load vector, resulting in Kd =
r. In practice, integration rules are applied to elements,
loads appear in the r vector, and nodal dof boundary
conditions may appear in the d vector or may be
partitioned out of the equation.
Solution methods for finite element matrix equations are
plentiful. In the case of the linear static Kd = r,
inverting K is computationally expensive and numerically
unstable. A better technique is Cholesky factorization, a form
of Gauss elimination, and a minor variation on the "LDU"
factorization theme. The K matrix may be efficiently
factored into LDU, where L is lower triangular,
D is diagonal, and U is upper triangular,
resulting in LDUd = r. Since L and D are
easily inverted, and U is upper triangular, d may
be determined by back-substitution. Another popular approach is
the wavefront method, which assembles and reduces the equations
at the same time. Some of the best modern solution methods
employ sparse matrix techniques. Because node-to-node
stiffnesses are non-zero only for nearby node pairs, the
stiffness matrix has a large number of zero entries. This can be
exploited to reduce solution time and storage by a factor of 10
or more. Improved solution methods are continually being
developed. The key point is that the analyst must understand the
solution technique being applied.
Dynamic analysis for too many analysts means normal modes.
Knowledge of the natural frequencies and mode shapes of a design
may be enough in the case of a single-frequency vibration of an
existing product or prototype, with FEA being used to
investigate the effects of mass, stiffness and damping
modifications. When investigating a future product, or an
existing design with multiple modes excited, forced response
modeling should be used to apply the expected transient or
frequency environment to estimate the displacement and even
dynamic stress at each time step.
This discussion has assumed h-code elements, for which the
order of the interpolation polynomials is fixed. Another
technique, p-code, increases the order iteratively until
convergence, with error estimates available after one analysis.
Finally, the boundary element method places elements only along
the geometrical boundary. These techniques have limitations, but
expect to see more of them in the near future.
Next month's article will discuss the post-processing phase of
the finite element method.