| A posteriori error estimation and adaptive mesh refinement for reliable finite element solutions |
Nagesh R. IyerEVER since the introduction of computers there has been a remarkable development of versatile and efficient numerical methods for computer-aided structural analysis. Among the different numerical methods, the finite element method (FEM) has been acknowledged as one of the most powerful methods for structural analysis. In this method the structure is visualized/modelled as an assemblage of building block-like elements, interconnected at the nodal points. For each element a standard set of simultaneous linear/nonlinear partial differential equations can be formulated to relate the nodal forces (equivalent to applied loads), and the unknown nodal displacements. The result is a set of simultaneous equations that are amenable for solution by computer. Upon implementing the loading and boundary conditions, the assembled set of equations can be solved for unknown nodal displacements. Substituting these values back to each element formulation provides the distributions of stress and displacement within each element. The method has been employed in many fields of application such as civil engineering, mechanical engineering, aerospace structures, electrical and electronics, bio-engineering/mechanics, etc. including study of interaction problems. Significant advances have been made since 1960s, in the development of new and efficient FE models and techniques which are economically viable and also practical to solve complex and large-size problems. The power and versatility of FEM are being exploited by developing suitable software packages and/or by using commercially available software.Structural Engineering Research Centre (SERC), CSIR Campus, Taramani, Chennai 600 113, India
The paper highlights the importance of reliability of finite element solutions in the context of acceptance of results by designers. The developments made in this direction by introducing the concept of a posteriori error estimation and adaptive mesh refinement for finite element analysis of structures are presented. The paper discusses the scientific basis and details of formulations proposed by researchers along with an illustrative example. The contributions made by SERC, Madras in this area of research are also mentioned. Future directions of research are also highlighted.
For correspondence. (e-mail: nri@cssercm.ren.nic.in)
A posteriori error estimation
Reliability of solutions in finite element analysis (FEA) depends on modelling, analysis procedures, and how well the errors are controlled at different stages of the solution process. The objective of the FE analyst has been to model the geometry, loading, material and other attribute information to closely represent the actual structure and its attribute information in order to compute the structural responses to an acceptable degree of accuracy with least computational effort. Despite significant advances made on the theories and algorithmic tools of discretization, selection of a computational model for a particular problem is largely based on intuition and experience gained from solving similar problems. Much effort is now being directed towards assessing the accuracy of numerical response predictions and improving their quality. In the analysis of complex structures such as ship structures, offshore platforms, tall towers and long span roofs, accurate determination of response is of considerable importance for economical and safe design of such structures. Analysis of any of these structures involves large amount of data to be prepared to model the physical structure. This gives scope for errors in input and inefficient modelling. The quality of results of analysis depends largely on the FE model employed. The discretization error, that is, the difference between the exact and the FE solution, is the result of modelling a continuum with a computational model that has a finite number of degrees of freedom. In FEM, the governing equation, that is the equilibrium conditions can be satisfied only in the weak sense and at a global level. Therefore, the stress solution is discontinuous across element boundaries. And hence, stresses obtained at boundaries are of lower quality than within the domain. Likewise, the imposed natural boundary conditions are not fulfilled in an exact manner. The accuracy or reliability of the results can be improved by better modelling. One can thus observe that errors of this nature are difficult to control straightaway, as it is not easy to estimate or compute the error in the solution obtained. To remove the element of uncertainty involved in modelling, research has been directed towards evolving methodologies for error estimation and adaptive mesh refinements that can be used in FEA. These attempts are intended to provide automatic means to improve mesh design and consequently to offer reliable solution with least computational cost. In a mesh refinement process, the user first selects a starting mesh and estimates the error in the solution corresponding to that mesh. If the error is found to be beyond the prescribed limits, another mesh is designed based on the starting mesh. This is called first level of mesh refinement. This cycle is continued till a satisfactory or acceptable solution (below the prescribed percentage of error in the solution) is found. This process is known as adaptive mesh refinement. Since the process of estimating errors in the solution and mesh refinement technique is based on a starting mesh and its FE solution, this way of estimating errors is known as a posteriori error estimation. A overview of the procedure commonly adopted for error estimation in the case of linear elastic static analysis using FEM is given in the following:
The solution of linear elastic static problems consists of displacements (d ) and stresses (s ). It may be noted that the stresses s are computed at Gauss points, as these are the best sampling points1. However, the displacements d and the stresses s are approximate solutions and differ from the exact values as given here:
ed = d * – d for displacement, (1)
es = s * – s for stresses. (2)
The differences in the corresponding values given in eqs (1) and (2) form the point-wise errors in displacements and stresses, respectively. However, since point-wise errors are difficult to compute, integral measures are conveniently adopted. Among these measures, the ‘energy norm’ is most commonly used. This can be expressed as,
(3)
For specific case of elastic models,
(4)
where B is the strain-displacement matrix, D is the elasticity matrix and ee and es are the error measures in strains and stresses, respectively.
A more direct measure is the so-called L2 norm, which can be associated with the errors in any quantity. This is expressed as,
(5)
and for stresses,
(6)
It can be observed that between the energy norm and L2 norm the expressions of error differ by the elasticity matrix, D. In the energy norm, error in each element will be
(7)
One can compute the total error in the entire domain by summing up the contribution of all elements,
(8)
where nel = number of elements
The relative element-wise error or the element-wise error indicator is defined as
(9)
The error indicator or relative percentage error can therefore be defined as,
(10)
where ||u|| is the exact energy norm.
The following steps may be used to estimate errors:
(1) Compute Gauss point stresses (s ).
(2) Compute smoothed/improved stresses at nodal points (s nd).
(3) Compute projected stresses at Gauss points using shape functions, i.e. project s nd to Gauss points s *nd.
(4) Compute the error estimate for the domain in
energy norm,
(11)
(5) Compute error indicator for the domain and elements using eqs (9) and (10).
The above error norms are similar to those proposed
by Zienkiewicz and Zhu2 and Iyer3. In the case of large and complex
problems, one cannot find the exact solution or obtain solutions by classical methods.
Since the exact values are not known, ‘improved’ ones that are supposed to be
better than the computed values replace the exact values. The procedures or approaches
generally used to determine the ‘improved’ solution are (i) global smoothing
approach, (ii) modified global smoothing approach, (iii) nodal averaging approach, and
(iv) superconvergent patch recovery technique. Iyer3 conducted an extensive
study on these approaches for plane stress/strain problems and found that the modified
global approach is more efficient than the global smoothing approach (especially for
biquadratic elements). The nodal averaging approach is computationally efficient and easy
for implementation. Recently, it was shown4 that if the recovered derivatives
are superconvergent, the Zienkiewicz and Zhu error estimator will always be asymptotically
exact in energy norm. It has been proved4 that the estimated error converges to
the true error if the improved solution is one or more order higher accurate than the FE
solution, a phenomenon known as superconvergence in the literature. Among these, the
methods that recover superconvergent solutions are interesting. The superconvergent patch
recovery (SPR), proposed by Zienkiewicz and Zhu4, is effective and recovers
superconvergent derivatives. The procedure is a least square fit of FE solutions over a
local patch of elements at preselected points, where the rate of convergence is higher
than the global rate or at least more accurate. The idea of SPR has also to some extent
been tested with success to dynamic problems by Wieberg and Li5. Recently, an
enhancement of the
superconvergent patch called superconvergent patch recovery with equilibrium (SPRE) was
achieved6 by adding residual of the equilibrium equation for the improved
solution to the smoothing procedure.
It is observed from the literature that work on development of error analysis procedures has been largely concerned with linear static FEA. For free vibration or dynamic response analysis the errors in the representation of mass as well as stiffness of the structure will significantly influence the quality or reliability of the solution. Errors in natural frequencies and mode shapes will obviously contribute significantly to errors in dynamic response analysis based on frequency domain approach. While considerable successes have been achieved on error estimation and adaptivity for elastostatic problems, the development concerning adaptive FEHs for elastodynamics and plasticity problems is far from complete and is very desirable7. For these problems gradients of the solution (stresses and strains) are often steep at some regions and smooth at others. These steep stress/strain regions may change their locations and shapes both in space and time. It is, therefore, advantageous to vary the spatial mesh and time steps according to the behaviour of the solution, to ensure that the mesh and time step are sufficiently fine in those regions and stages of steep gradients and reasonably coarse in other regions and stages of less steep gradients. In this way numerical computations for a given engineering problem can be carried out economically and efficiently keeping the discretization errors under control and the computational cost nearly minimal. Indeed a posteriori error estimate and adaptive approaches for nonlinear and nonelliptic problems are considered to be one of the most outstanding problems in FEs today. Considerable advances have been made at SERC, Madras8,9 in developing methodologies and algorithms for estimation of errors in FEA of problems in vibration and dynamic response analysis.
Adaptive mesh refinement process
Adaptive improvement of FE solution aims at improving the quality of the solutions by enriching the approximation in some manner that will be enable one to achieve the ‘best’ solution for a given computational effort. Four adaptive strategies have been proposed in the literature. These are:
(1) Successive selection of the meshes by element size refinement (h method).
(2) Moving the nodes (R method).
(3) Successive selection of the order of the poly-
nomial approximation inside some elements (p method).
(4) Simultaneous selection of meshes and local order of the approximation (h–p method).
The maturation of the methods of estimation and
control of discretization errors and the incorporation of these methods into
general-purpose FE systems will
allow the analyst to select only:
(i) the initial discrete model, which is sufficient to resolve the topology of the structure;
(ii) the error measure and the tolerance.
Then the FE software can automatically refine the model until the selected error measure falls below the prescribed tolerance.
Errors estimated based on the results at the current solution stage will serve as a guide to start the adaptive refinement process. In the case of static analysis, the adaptive strategy is concerned with refinement of stiffness representation while in free vibration and dynamic analysis, the adaptive strategy has to deal with refinement of stiffness as well as mass by modifying the FE mesh. The objective of the self-adaptive refinement is to allow the inexperienced user to handle the problem with ease and to achieve reliable solution with required accuracy. It is therefore, important to provide effective and powerful refinement strategy in conjunction with the error analysis that is adopted. One should be able to achieve maximum possible reduction of error in one level of refinement and thereby reduce the number of levels of refinement, which in turn reduces the computational effort.
Initially the structure is idealized by using a
basic/
starting mesh and while carrying out static/vibration analysis error estimation procedures
are used. Based on the elemental error indicators the element needing refinement is
sub-divided equally. Keeping the original node numbering intact, new numbers are assigned
to the nodes as they are being created and element numbers are re-defined appropriately.
During this process, for example, a 4-noded element is divided into four new 4-noded
elements and five new nodes are created for each element refined. Two types of nodes are
created. The first are the nodes that have independent degrees of freedom and the second
are those nodal displacements that are dependent on the neighbouring nodes. It is the
second type which needs special treatment. It may be noted that the FEs are no longer
conforming under such circumstances. These nodes are known as ‘dependent’ nodes.
As far as linear static case is concerned, some schemes are being used to deal with such
situations. These schemes are not readily applicable for vibration and dynamic analysis,
as the problem is now concerned with modal characteristics. After a detailed study, Iyer3
has proposed a scheme of distribution of mass at the constrained nodes to active nodes
equally at each refinement level. An illustration of this is shown in Figure 1. All
elements are initially shaped like a quadrilateral. It is assumed that element number 2
does not satisfy the specified accuracy during the initial analysis and therefore, needs
refinement. Element number 2 in Figure 1 a is the parent element,
which is subdivided into four elements on first level of refinement in Figure 1 b.
Element number 2 has global C0 approximation and is piece-wise linear
along side ‘3-d-4’ in the refined elements. However, the same cannot be
true in the case of element number 1. In the present approach, referring to Figure 1 b,
the dependent node d is constrained in order to maintain inter-element
compatibility which will make the elements conforming. While constraining the node d,
the nodal variables of d are expressed in terms of nodes 3 and 4 as:
(12)
This requires transformation of usual shape functions to constrained shape functions as
(13)
where N = shape function and NC = constrained shape function.
(14)
for 4-noded (bilinear) plane stress/strain element.
Illustrative example
A problem of a square plate with a central square
hole or opening under uniform tension considering quarter symmetry is chosen to
demonstrate the concepts and the procedure. Geometry, loading and material property
information are provided in Figure 2. The analysis has been carried out with uniform
meshes and with h-adaptive refinement procedure. This is done to get an idea of the
effectiveness of the error estimation and adaptive refinement procedures over the
finely-graded FE mesh. To start with, uniform meshes have been used to conduct studies on
proposed procedures of error analysis. The user-defined permissible error in the solution
is taken as 5%. The basic starting mesh used for adaptive refinement is shown in Figure 3.
The uniform meshes that are generated to get an accuracy of about 5% in the solution are
shown in Figure 3 a–d. The resulting mesh obtained
using the proposed h-adaptive
refinement algorithm is shown in Figure 3 e. The performance of the error estimation using the five levels of refinement for uniform meshes and the adaptive mesh refinement (in this case only one level) is shown in Figure 3 f.
A simple example demonstrated here shows that finer uniform meshes need not yield reliable solution. One does not need to have fine mesh covering the full geometry but only in places where the errors are more and need to be controlled. In this example obviously one expects finer mesh at the corner where the stress concentration is high. The adaptive mesh refinement in Figure 3 e shows this. Another point of significance is that the user knows at what stage of mesh refinement the errors are under control and when to stop. With the other approach of using uniform meshes, one keeps increasing the number of equations (or Degrees of Freedom, DOF) and also the computational effort. Further, it gives no assurance of having obtained a reliable solution. From Figure 3 f one can also observe the faster rate of convergence (or reduction) of errors using adaptive mesh refinement compared to uniform meshes.
Concluding remarks
The efficacy of the proposed error estimators and the adaptive refinement has been demonstrated by conducting convergence studies for static analysis of plane stress problems. However, the proposed error estimators and the adaptive refinement strategy are general and can be used appropriately for analysis of any structure and also can be extended to other types of problems. The importance of error analysis is relevant in the current context as a number of decisions are being made on the basis of FEA. In view of the above, it is felt that research effort has to be directed towards developing reliable methods for error analysis and adaptive refinements for FEA of structures.
ACKNOWLEDGEMENTS. I thank Dr T.V.S.R. Appa Rao, Director, SERC and his colleagues Mr J. Rajasankar and Mr G. S. Palani for their valuable inputs in the preparation of this manuscript which is being published with the kind permission of the Director SERC, Chennai.