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.
copyright 2005 Roensch & Associates. All rights reserved.
