The influence that the order of the element (linear or quadratic), the element formulation, and the level of integration have on the accuracy of a structural simulation will be demonstrated by considering a static analysis of the cantilever beam shown in Figure 1.
This is a classic test used to assess the behavior of a given finite element. Since the beam is relatively slender, we would normally model it with beam elements. However, it is used here to help assess the effectiveness of various solid elements.
The beam is 150 mm long, 2.5 mm wide, and 5 mm deep; built-in at one end; and carries a tip load of 5 N at the free end. The material has Young's modulus, E, of 70 GPa and a Poisson's ratio of 0.0. Using beam theory, the static deflection of the tip of the beam for a load P is given as
where,
l is the length, b is the width, and d is the depth of the beam.
For P = 5 N the tip deflection is 3.09 mm.
Full integration
The expression “full integration” refers to the number of Gauss points required to integrate the polynomial terms in an element's stiffness matrix exactly when the element has a regular shape. For hexahedral and quadrilateral elements a “regular shape” means that the edges are straight and meet at right angles and that any edge nodes are at the midpoint of the edge. Fully integrated, linear elements use two integration points in each direction. Thus, the three-dimensional element C3D8 uses a 2 × 2 × 2 array of integration points in the element. Fully integrated, quadratic elements (available only in Abaqus/Standard) use three integration points in each direction. The locations of the integration points in fully integrated, two-dimensional, quadrilateral elements are shown in Figure 2.
Figure 2 Integration points in fully integrated, two-dimensional, quadrilateral elements.
Several different finite element meshes were used in Abaqus/Standard simulations of the cantilever beam problem, as shown in Figure 3. The simulations use either linear or quadratic, fully integrated elements and illustrate the effects of both the order of the element (first versus second) and the mesh density on the accuracy of the results.
Figure 3 Meshes used for the cantilever beam simulations.
The ratios of the tip displacements for the various simulations to the beam-theory value of 3.09 mm are shown in Table 1.
Table 1 Normalized tip displacements with fully-integrated elements.
Element | Mesh Size (Depth × Length) | |||
1 × 6 | 2 × 12 | 4 × 12 | 8 × 24 | |
CPS4 | 0.074 | 0.242 | 0.242 | 0.561 |
CPS8 | 0.994 | 1.000 | 1.000 | 1.000 |
C3D8 | 0.077 | 0.248 | 0.243 | 0.563 |
C3D20 | 0.994 | 1.000 | 1.000 | 1.000 |
The linear elements CPS4 and C3D8 underpredict the deflection so badly that the results are unusable. The results are the least accurate with coarse meshes, but even a fine mesh (8 × 24) still predicts a tip displacement that is only 56% of the theoretical value. Notice that for the linear, fully integrated elements it makes no difference how many elements there are through the thickness of the beam. The underprediction of tip deflection is caused by shear locking, which is a problem with all fully integrated, first-order, solid elements.
As we have seen, shear locking causes the elements to be too stiff in bending. It is explained as follows. Consider a small piece of material in a structure subject to pure bending. The material will distort as shown in Figure 4.
Figure 4 Deformation of material subjected to bending moment M.
Lines initially parallel to the horizontal axis take on constant curvature, and lines through the thickness remain straight. The angle between the horizontal and vertical lines remains at 90°.
The edges of a linear element are unable to curve; therefore, if the small piece of material is modeled using a single element, its deformed shape is like that shown in Figure 5.
Figure 5 Deformation of a fully integrated, linear element subjected to bending moment M.
For visualization, dotted lines that pass through the integration points are plotted. It is apparent that the upper line has increased in length, indicating that the direct stress in the 1-direction, σ11, is tensile. Similarly, the length of the lower dotted line has decreased, indicating that σ11 is compressive. The length of the vertical dotted lines has not changed (assuming that displacements are small); therefore, σ22 at all integration points is zero. All this is consistent with the expected state of stress of a small piece of material subjected to pure bending. However, at each integration point the angle between the vertical and horizontal lines, which was initially 90°, has changed. This indicates that the shear stress, σ12, at these points is nonzero. This is incorrect: the shear stress in a piece of material under pure bending is zero.
This spurious shear stress arises because the edges of the element are unable to curve. Its presence means that strain energy is creating shearing deformation rather than the intended bending deformation, so the overall deflections are less: the element is too stiff.
Shear locking only affects the performance of fully integrated, linear elements subjected to bending loads. These elements function perfectly well under direct or shear loads. Shear locking is not a problem for quadratic elements since their edges are able to curve (see Figure 6). The predicted tip displacements for the quadratic elements shown in Table 1 are close to the theoretical value. However, quadratic elements will also exhibit some locking if they are distorted or if the bending stress has a gradient, both of which can occur in practical problems.
Figure 4–6 Deformation of a fully integrated, quadratic element subjected to bending moment M.
Fully integrated, linear elements should be used only when you are fairly certain that the loads will produce minimal bending in your model. Use a different element type if you have doubts about the type of deformation the loading will create. Fully integrated, quadratic elements can also lock under complex states of stress; thus, you should check the results carefully if they are used exclusively in your model. However, they are very useful for modeling areas where there are local stress concentrations.
Volumetric locking is another form of overconstrained that occurs in fully integrated elements when the material behavior is (almost) incompressible. It causes overly stiff behavior for deformations that should cause no volume changes.
Reduced integration
Only quadrilateral and hexahedral elements can use a reduced-integration scheme; all wedge, tetrahedral, and triangular solid elements use full integration, although they can be used in the same mesh with reduced-integration hexahedral or quadrilateral elements.
Reduced-integration elements use one fewer integration points in each direction than the fully integrated elements. Reduced integration, linear elements have just a single integration point located at the element's centroid. (Actually, these first-order elements in Abaqus use the more accurate “uniform strain” formulation, where average values of the strain components are computed for the element. This distinction is not important for this discussion.) The locations of the integration points for reduced-integration, quadrilateral elements are shown in Figure 7.
Figure 7 Integration points in two-dimensional elements with reduced integration.
Abaqus simulations of the cantilever beam problem were performed using the reduced-integration versions of the same four elements utilized previously and using the four finite element meshes shown in Figure 7. The results from these simulations are presented in Table 2.
Table 2 Normalized tip displacements with reduced-integration elements.
Element | Mesh Size (Depth × Length) | |||
1 × 6 | 2 × 12 | 4 × 12 | 8 × 24 | |
CPS4R | 20.3* | 1.308 | 1.051 | 1.012 |
CPS8R | 1.000 | 1.000 | 1.000 | 1.000 |
C3D8R | 70.1* | 1.323 | 1.063 | 1.015 |
C3D20R | 0.999** | 1.000 | 1.000 | 1.000 |
* no stiffness to resist the applied load, ** two elements through width |
Linear reduced-integration elements tend to be too flexible because they suffer from their own numerical problem called hourglassing. Again, consider a single reduced-integration element modeling a small piece of material subjected to pure bending (see Figure 8).
Figure 8 Deformation of a linear element with reduced integration subjected to bending moment M.
Neither of the dotted visualization lines has changed in length, and the angle between them is also unchanged, which means that all components of stress at the element's single integration point are zero. This bending mode of deformation is thus a zero-energy mode because no strain energy is generated by this element distortion. The element is unable to resist this type of deformation since it has no stiffness in this mode. In coarse meshes this zero-energy mode can propagate through the mesh, producing meaningless results.
In Abaqus a small amount of artificial “hourglass stiffness” is introduced in first-order reduced-integration elements to limit the propagation of hourglass modes. This stiffness is more effective at limiting the hourglass modes when more elements are used in the model, which means that linear reduced-integration elements can give acceptable results as long as a reasonably fine mesh is used. The errors seen with the finer meshes of linear reduced-integration elements (see Table 2) are within an acceptable range for many applications. The results suggest that at least four elements should be used through the thickness when modeling any structures carrying bending loads with this type of element. When a single linear reduced-integration element is used through the thickness of the beam, all the integration points lie on the neutral axis and the model is unable to resist bending loads. (These cases are marked with a * in Table 2.)
Linear reduced-integration elements are very tolerant of distortion; therefore, use a fine mesh of these elements in any simulation where the distortion levels may be very high.
The quadratic reduced-integration elements available in Abaqus/Standard also have hourglass modes. However, the modes are almost impossible to propagate in a normal mesh and are rarely a problem if the mesh is sufficiently fine. The 1 × 6 mesh of C3D20R elements fails to converge because of hourglassing unless two elements are used through the width, but the more refined meshes do not fail even when only one element is used through the width. Quadratic reduced-integration elements are not susceptible to locking, even when subjected to complicated states of stress. Therefore, these elements are generally the best choice for most general stress/displacement simulations, except in large-displacement simulations involving very large strains and in some types of contact analyses.