Finite element method in structural mechanics

The finite element method (FEM) is a powerful technique originally developed for numerical solution of complex problems in structural mechanics, and it remains the method of choice for complex systems. In the FEM, the structural system is modeled by a set of appropriate finite elements interconnected at points called nodes. Elements may have physical properties such as thickness, coefficient of thermal expansion, density, Young's modulus, shear modulus and Poisson's ratio.

History

The origin of finite method can be traced to the matrix analysis of structures ^[1][2] where the concept of displacement or stiffness matrix approach was introduced. Finite element concepts were developed based on engineering methods in 1950s. The finite element method obtained its real impetus in the 1960s and 1970s by the developments of John Argyris with co-workers at the University of Stuttgart, Ray W. Clough with co-workers at University of California, Berkeley, Olgierd Zienkiewicz with co-workers Ernest Hinton, Bruce Irons ^[3] and others at the University of Swansea, Philippe G. Ciarlet at the University of Paris and Richard Gallagher with co-workers at Cornell University. The original works such as those by Argyris ^[4] and Clough ^[5] became foundation for today’s finite element structural analysis methods. Earlier books such as by Zienkiewicz ^[6] and more recent books such as by Yang ^[7] give comprehensive summary of developments in finite-element structural analysis. Implementing the method in software is described in the classic text by Smith, Griffiths and Margetts. ^[8]

Element properties

• Straight or curved one-dimensional elements with physical properties such as axial, bending, and torsional stiffnesses. This type of element is suitable for modeling cables, braces, trusses, beams, stiffeners, grids and frames. Straight elements usually have two nodes, one at each end, while curved elements will need at least three nodes including the end-nodes. The elements are positioned at the centroidal axis of the actual members.
• Two-dimensional elements for membrane action (plane stress, plane strain) and/or bending action (plates and shells). They may have a variety of shapes such as flat or curved triangles and quadrilaterals. Nodes are usually placed at the element corners and, if needed for higher accuracy, additional nodes can be placed along the element edges or even inside the element. The elements are positioned at the mid-surface of the actual layer thickness.
• Torus-shaped elements for axisymmetric problems such as thin, thick plates, shells, and solids. The cross-section of the elements are similar to the previously described types: one-dimensional for thin plates and shells, and two-dimensional for solids, and thick plates and shells.
• Three-dimensional elements for modeling 3-D solids such as machine components, dams, embankments or soil masses. Common element shapes include tetrahedrals and hexahedrals. Nodes are placed at the vertexes and possibly in the element faces or within the element.

Element interconnection and displacement

The elements are interconnected only at the exterior nodes, and altogether they should cover the entire domain as accurately as possible. Nodes will have nodal (vector) displacements or degrees of freedom which may include translations, rotations, and for special applications, higher order derivatives of displacements. When the nodes displace, they will drag the elements along in a certain manner dictated by the element formulation. In other words, displacements of any points in the element will be interpolated from the nodal displacements, and this is the main reason for the approximate nature of the solution.

Practical considerations

From the application point of view, it is important to model the system such that:

• Symmetry or anti-symmetry conditions are exploited in order to reduce the size of the domain.
• Displacement compatibility, including any required discontinuity, is ensured at the nodes, and preferably, along the element edges as well, particularly when adjacent elements are of different types, material or thickness. Compatibility of displacements of many nodes can usually be imposed via constraint relations—When such a feature is not available in the software package, a physical model that imposes the constraints may be used instead.
• Elements' behaviours capture the dominant actions of the actual system, both locally and globally.
• The element mesh is sufficiently fine in order to have acceptable accuracy. To assess accuracy, the mesh is refined until the important results shows little change. For higher accuracy, the aspect ratio of the elements should be as close to unity as possible, and smaller elements are used over the parts of higher stress gradient.
• Proper support constraints are imposed with special attention paid to nodes on symmetry axes.

Large scale commercial software packages often provide facilities for generating the mesh, graphical display of input and output, which greatly facilitate the verification of both input data and interpretation of the results.

Theoretical overview of FEM-Displacement Formulation: From elements to system to solution

While the theory of FEM can be presented in different perspectives or emphases, its development for structural analysis follows the more traditional approach via the virtual work principle or the minimum total potential energy principle. The virtual work principle approach is more general as it is applicable to both linear and non-linear material behaviours.

The principle of virtual displacements for the structural system expresses the mathematical identity of external and internal virtual work:

$\mbox{External virtual work} = \int_{V}\delta\boldsymbol{\epsilon}^T \boldsymbol{\sigma} \, dV \qquad \mathrm{(1)}$

The virtual internal work in the right-hand-side of the above equation may be found by summing the virtual work in the individual elements—This is the crucial step where we will need displacement functions written only for the small domain rather than over the entire system. As shown in the subsequent sections, Eq.(1) leads to the following governing equilibrium equation for the system:

$\mathbf{R} = \mathbf{Kr} + \mathbf{R}^o \qquad \qquad \qquad \mathrm{(2)}$

where

$\mathbf{R}$ = vector of nodal forces, representing external forces applied to the system's nodes.

$\mathbf{r}$ = vector of system's nodal displacements, which will, by interpolation, yield displacements at any point of the finite element mesh.

$\mathbf{R}^o$ = vector of equivalent nodal forces, representing all external effects other than the nodal forces which are already included in the preceding nodal force vector R. These external effects may include distributed or concentrated surface forces, body forces, thermal effects, initial stresses and strains.

$\mathbf{K}$ = system stiffness matrix, which will be established by assembling the elements' stiffness matrices :$\mathbf{k}^e$.

Once the supports' constraints are accounted for, the nodal displacements are found by solving the system of linear equations (2), symbolically:

$\mathbf{r} = \mathbf{K}^{-1} (\mathbf{R}-\mathbf{R}^o ) \qquad \qquad \qquad \mathrm{(3)}$

Subsequently, the strains and stresses in individual elements may be found as follows:

$\mathbf{\epsilon} = \mathbf{Bq} \qquad \qquad \qquad \qquad \mathrm{(4)}$

$\mathbf{\sigma} = \mathbf{E}(\mathbf{\epsilon} - \mathbf{\epsilon}^o)+\mathbf{\sigma}^o = \mathbf{E}(\mathbf{Bq} - \mathbf{\epsilon}^o)+\mathbf{\sigma}^o\qquad \qquad \qquad \mathrm{(5)}$

where

$\mathbf{q}$ = vector of element's nodal displacements--a subset of the system displacement vector r that pertains to the element under consideration.

$\mathbf{B}$ = strain-displacement matrix that transforms nodal displacements q to strains at any point in the element.

$\mathbf{E}$ = elasticity matrix that transforms effective strains to stresses at any point in the element.

$\mathbf{\epsilon}^o$ = vector of initial strains in the element.

$\mathbf{\sigma}^o$ = vector of initial stresses in the element.

By applying the virtual work equation (1) to the system, we can establish the element matrices $\mathbf{B}$, $\mathbf{k}^e$ as well as the technique of assembling the system matrices $\mathbf{R}^o$ and $\mathbf{K}$. Other matrices such as $\mathbf{\epsilon}^o$, $\mathbf{\sigma}^o$, $\mathbf{R}$ and $\mathbf{E}$ can be directly set up from data input.

Interpolation or shape functions

Let $\mathbf{q}$ be the vector of nodal displacements of a typical element. The displacements at any point of the element may be found by interpolation functions as, symbolically:

$\mathbf{u} = \mathbf{N} \mathbf{q} \qquad \qquad \qquad \mathrm{(6)}$

where

$\mathbf{u}$ = vector of displacements at any point {x,y,z} of the element.

$\mathbf{N}$ = matrix of shape functions serving as interpolation functions.

Equation (6) gives rise to other quantities of great interest:

• Virtual displacements consistent with virtual nodal displacements: $\delta \mathbf{u} = \mathbf{N} \delta \mathbf{q} \qquad \qquad \qquad \mathrm{(6b)}$
• Strains in the elements:$\mathbf{\epsilon} = \mathbf{Du} = \mathbf{DNq} \qquad \qquad \qquad \qquad \mathrm{(7)}$

where $\mathbf{D}$ = matrix of differential operators that convert displacements to strains using linear elasticity theory. Eq.(7) shows that matrix B in (4) is

$\mathbf{B} = \mathbf{DN} \qquad \qquad \qquad \qquad \mathrm{(8)}$

• Virtual strains consistent with element's virtual nodal displacements: $\delta \boldsymbol{\epsilon} = \mathbf{B} \delta \mathbf{q} \qquad \qquad \qquad \qquad \mathrm{(9)}$

Internal virtual work in a typical element

For a typical element of volume $V^e$, the internal virtual work due to virtual displacements is obtained by substitution of (5) and (9) into (1):

$\mbox{Internal virtual work} = \int_{V^e}\delta\boldsymbol{\epsilon}^T \boldsymbol{\sigma} \, dV^e = \delta\ \mathbf{q}^T \int_{V^e} \mathbf{B}^T \big\{\mathbf{E}(\mathbf{Bq} - \mathbf{\epsilon}^o)+\mathbf{\sigma}^o\big\} \, dV^e \qquad \mathrm{(10)}$

Element matrices

Primarily for the convenience of reference, the following matrices pertaining to a typical elements may now be defined:

Element stiffness matrix $\mathbf{k}^e = \int_{V^e} \mathbf{B}^T \mathbf{E} \mathbf{B} \, dV^e \qquad \mathrm{(11)}$

Equivalent element load vector $\mathbf{Q}^{oe} = \int_{V^e} - \mathbf{B}^T \big( \mathbf{E}\mathbf{\epsilon}^o - \mathbf{\sigma}^o\big ) \, dV^e \qquad \mathrm{(12)}$

These matrices are usually evaluated numerically using Gaussian quadrature for numerical integration. Their use simplifies (10) to the following:

$\mbox{Internal virtual work} = \delta\ \mathbf{q}^T \big( \mathbf{k}^e \mathbf{q} + \mathbf{Q}^{oe} \big) \qquad \mathrm{(13)}$

Element virtual work in terms of system nodal displacements

Since the nodal displacement vector q is a subset of the system nodal displacements r (for compatibility with adjacent elements), we can replace q with r by expanding the size of the element matrices with new columns and rows of zeros:

$\mbox{Internal virtual work} = \delta\ \mathbf{r}^T \big( \mathbf{k}^e \mathbf{r} + \mathbf{Q}^{oe} \big) \qquad \mathrm{(14)}$

where, for simplicity, we use the same symbols for the element matrices, which now have expanded size as well as suitably rearranged rows and columns.

System virtual work

Summing the internal virtual work (14) for all elements gives the right-hand-side of (1):

$\mbox{System internal virtual work} = \sum_{e} \delta\ \mathbf{r}^T \big( \mathbf{k}^e \mathbf{r} + \mathbf{Q}^{oe} \big) = \delta\ \mathbf{r}^T \big( \sum_{e} \mathbf{k}^e \big)\mathbf{r} + \delta\ \mathbf{r}^T \sum_{e} \mathbf{Q}^{oe} \qquad \mathrm{(15)}$

Considering now the left-hand-side of (1), the system external virtual work consists of:

• The work done by the nodal forces R: $\delta\ \mathbf{r}^T \mathbf{R} \qquad \mathrm{(16)}$
• The work done by external forces $\mathbf{T}^e$ on the part $\mathbf{S}^e$ of the elements' edges or surfaces, and by the body forces $\mathbf{f}^e$

$\sum_{e} \int_{S^e} \delta\ \mathbf{u}^T \mathbf{T}^e \, dS^e + \sum_{e} \int_{V^e} \delta\ \mathbf{u}^T \mathbf{f}^e \, dV^e$

Substitution of (6b) gives:

$\delta\ \mathbf{q}^T \sum_{e} \int_{S^e} \mathbf{N}^T \mathbf{T}^e \, dS^e + \delta\ \mathbf{q}^T \sum_{e} \int_{V^e} \mathbf{N}^T \mathbf{f}^e \, dV^e$

or $-\delta\ \mathbf{q}^T \sum_{e} (\mathbf{Q}^{te} + \mathbf{Q}^{fe}) \qquad \mathrm{(17a)}$

where we have introduced additional element's matrices defined below:

$\mathbf{Q}^{te} = -\int_{S^e} \mathbf{N}^T \mathbf{T}^e \, dS^e \qquad \mathrm{(18a)}$

$\mathbf{Q}^{fe} = -\int_{V^e} \mathbf{N}^T \mathbf{f}^e \, dV^e \qquad \mathrm{(18b)}$

Again, numerical integration is convenient for their evaluation. A similar replacement of q in (17a) with r gives, after rearranging and expanding the vectors $\mathbf{Q}^{te}, \mathbf{Q}^{fe}$:

$-\delta\ \mathbf{r}^T \sum_{e} (\mathbf{Q}^{te} + \mathbf{Q}^{fe}) \qquad \mathrm{(17b)}$

Assembly of system matrices

Adding (16), (17b) and equating the sum to (15) gives: $\delta\ \mathbf{r}^T \mathbf{R} -\delta\ \mathbf{r}^T \sum_{e} (\mathbf{Q}^{te} + \mathbf{Q}^{fe}) = \delta\ \mathbf{r}^T \big( \sum_{e} \mathbf{k}^e \big)\mathbf{r} + \delta\ \mathbf{r}^T \sum_{e} \mathbf{Q}^{oe}$

Since the virtual displacements $\delta\ \mathbf{r}$ are arbitrary, the preceding equality reduces to:

$\mathbf{R} = \big( \sum_{e} \mathbf{k}^e \big)\mathbf{r} + \sum_{e} \big( \mathbf{Q}^{oe} + \mathbf{Q}^{te} + \mathbf{Q}^{fe} \big)$

Comparison with (2) shows that:

• The system stiffness matrix is obtained by summing the elements' stiffness matrices:

$\mathbf{K} = \sum_{e} \mathbf{k}^e$

• The vector of equivalent nodal forces is obtained by summing the elements' load vectors:

$\mathbf{R}^o = \sum_{e} \big( \mathbf{Q}^{oe} + \mathbf{Q}^{te} + \mathbf{Q}^{fe} \big)$

In practice, the element matrices are neither expanded nor rearranged. Instead, the system stiffness matrix $\mathbf{K}$ is assembled by adding individual coefficients ${k}_{ij}^e$ to ${K}_{kl}$ where the subscripts ij, kl mean that the element's nodal displacements ${q}_{i}^e, {q}_{j}^e$ match respectively with the system's nodal displacements ${r}_{k}, {r}_{l}$. Similarly, $\mathbf{R}^o$ is assembled by adding individual coefficients ${Q}_{i}^e$ to ${R}^o_{k}$ where ${q}_{i}^e$ matches ${r}_{k}$. This direct addition of ${k}_{ij}^e$ into ${K}_{kl}$ gives the procedure the name Direct Stiffness Method.

See also

• Finite element method
• Flexibility method
• Matrix stiffness method
• Modal analysis using FEM
• List of finite element software packages
• Structural analysis
• Virtual work
• Interval finite element

References

1. Matrix Analysis Of Framed Structures, 3rd Edition by Jr. William Weaver, James M. Gere, Springer-Verlag New York, LLC, ISBN 978-0-412-07861-3, 1966
2. Theory of Matrix Structural Analysis, J. S. Przemieniecki, McGraw-Hill Book Company, New York, 1968
3. Hinton, Ernest; Irons, Bruce (July 1968). "Least squares smoothing of experimental data using finite elements". Strain 4: 24–27.
4. Argyris, J.H and Kelsey, S. Energy theorems and Structural Analysis Butterworth Scientific publications, London
5. Clough, R.W, “The Finite Element in Plane Stress Analysis.” Proceedings, 2nd ASCE Conference on Electronic Computations, Pittsburgh, Sep 1960
6. The Finite Element Method for Solid and Structural Mechanics, Zienkiewicz O. C and Taylor R L ISBN 978-0-7506-6321-2, 1967, McGraw Hill, New York
7. Finite Element Structural Analysis , T.Y Yang, Prentice-Hall, Inc, Englewood, NJ, 1986
8. "Wiley: Programming the Finite Element Method, 5th Edition - I. M. Smith, D. V. Griffiths, L. Margetts". eu.wiley.com. Retrieved 2015-09-18.