An improved stress recovery technique for low-order 3D finite elements

In this paper, we propose a stress recovery procedure for low-order finite elements in 3D. For each finite element, the recovered stress field is obtained by satisfying equilibrium in an average sense and by projecting the directly calculated stress field onto a conveniently chosen space. Compared with existing recovery techniques, the current procedure gives more accurate stress fields, is simpler to implement, and can be applied to different types of elements without further modification. We demonstrate, through a set of examples in linear elasticity, that the recovered stresses converge at a higher rate than that of directly calculated stresses and that, in some cases, the rate of convergence is the same as that of the displacement field.


INTRODUCTION
In a displacement-based finite element (FE) analysis, the stress is usually computed from the displacement field in a postprocessing step. This directly calculated stress is obtained from the derivative of the displacement field, leading inevitably to a lower-order less-accurate field as compared to that of the displacement. In addition, the usual C 0 -continuity at interelement boundaries leads to a discontinuous stress field. As a result, a significant amount of research has been devoted to retrieve the lost accuracy and continuity of the stress field, leading to methods known as "stress recovery techniques." Improved stress fields obtained through these procedures not only provide a better stress representation for reduced computational cost but also are widely used in a posteriori error estimates.
Initial techniques were based on global and local smoothing using the least-squares method. For instance, Hinton and Campbell 1 proposed a procedure in which global or local nodal averaging is used to minimize the error between the improved nodal stress and the FE stress. Similar to global smoothing, Oden and Brauchli 2 proposed a procedure that, assuming a C 0 continuous interpolation of the stress, computes nodal stresses by a least-square fit to the directly calculated stress. Global smoothing techniques are computationally expensive, difficult to implement, and for quadratic elements, results were discouraging. 3 The superconvergent patch recovery (SPR) of Zienkiewicz and Zhu, 3 which outperformed the earlier smoothing methods, became widely adopted and was further explored as an a posteriori error estimate. 4 In this recovery technique, the stress at mesh nodes is obtained by using a complete polynomial of the same order as that used for approximating the displacement field. The polynomial is then fitted in a least-square sense to the directly calculated stresses computed at superconvergent points over a patch of elements. Indeed, at these points, the accuracy of the derivatives is one order higher than at any other point in the element, resulting in a substantial improvement in the corresponding stresses. Superconvergent patch recovery was later further developed to improve its accuracy by adding an equilibrium constraint to the least-squares functional. 5 This enabled the use of higher-order polynomials and thus resulted in more accurate recovered stresses. A similar method by Wiberg et al. 6 proposed an additional improvement by adding the equilibrium equation and the prescribed tractions as constraints on the stress field. However, because of the additional constraints to the least-square functional and the use of higher-order polynomials for the recovered stresses, these improvements required more computational resources. All in all, SPR-based methods have an important drawback, ie, they rely on superconvergent points that are not always guaranteed to exist, eg, in non-regular isoparametric elements, elements with varying material properties, 7 and in linear triangular elements. 8 Several other methods that do not require superconvergent points were later proposed. In the recovery by equilibrium in patches (REP) by Boroomand and Zienkiewicz,8,9 the improved stresses are obtained by minimizing, over a patch of elements and in a least-square sense, the error between the improved and the FE stresses projected onto the space of FE strains. To deal with instability and singularity issues of small patches, an equilibrium term is added to the original minimization formulation of each element. An improved REP method was later proposed by the same authors, 8 in which the additional equilibrium term for each element was no longer needed. This leads to a simpler and efficient recovery technique as compared with the original REP formulation. Ubertini then proposed the recovery by compatibility in patches (RCP), 10 which exceeds the performance of both REP and SPR in terms of accuracy and sensitivity to mesh geometry changes. In RCP, the stresses are improved over a patch of elements by minimizing the error between improved strains and the strains from the FE solution. RCP requires an a priori problem-dependent particular solution to the continuum equation of equilibrium. Payen and Bathe introduced then the nodal-point force (NPF) method. 11,12 The improved stresses in NPF are calculated from the element nodal point forces using two virtual work statements, which express that the virtual work of element nodal point forces is equal to the virtual work of element boundary tractions and to the internal virtual work, respectively. This method was later extended to 3D for 4-node tetrahedral elements. 12 The procedure is not only straightforward to implement but also more accurate because of the higher-order stress field adopted. However, NPF also has drawbacks, which include the need of stress averaging and an element-type specific patch. As a result, Payen and Bathe later proposed another stress improvement procedure (SIP) that alleviates these shortcomings. 13 SIP uses a mixed formulation based on the Hu-Washizu principle 14 where the stress-strain relationship is relaxed pointwise but the fulfilment of equilibrium is enhanced. In their paper, they showed that, under certain assumptions, REP, RCP, and NPF can be derived from this mixed formulation. To date, their stress recovery procedure outperforms SPR, REP, RCP, and NPF. Nonetheless, the authors provided the interpolation matrices needed to implement the stress recovery procedure only in 1D and 2D.
In this paper, we extend SIP and investigate its performance on 3D linear FE analysis. We focus on 4-node tetrahedral and 8-node hexahedral elements and provide the required interpolation matrices needed for the computation of the stress coefficients that describe the improved stress field. With the help of numerical examples, we also investigate the convergence rate of the resulting approach. Even though the examples provided here are for problems in linear elastostatics, the method is general and thus could also be applied in a similar manner to nonlinear and dynamic problems.

STRESS IMPROVEMENT PROCEDURE
Here, we provide only the equations for obtaining the enhanced stress and we give the full derivation in Appendix A. To compute the enhanced stresses we start from a displacement field that was computed by standard procedures. The enhanced stress (m) in the mth element is computed by the following two projection equations: where the summations involve a patch of N p elements, called the stress calculation domain V p , f B is the known body force, (m) and (m) are Lagrangian multipliers, and (m) h is the directly calculated FE stress. In Equation (1), ie, (m) ∈ is the vector-valued function space of first-order polynomials, from which the displacement field is also chosen. The enhanced stress (m) ∈  2 = {T| T i ∈ P 2 (V ); 1 ≤ i, ≤ 3}, ie, (m) is a member of the second-order-tensor-valued function space of second-order polynomials. Similarly, in Equation (2), ie, (m) ∈ 2 , Equation (1) implies that the enhanced stress field is obtained by satisfying equilibrium in a weak sense over the stress calculation domain V p . Yet, this equation alone is not sufficient to compute the enhanced stress field because any constant field added to (m) would also satisfy the equation. This arbitrariness is removed by Equation (2), which ensures that the difference between enhanced and directly calculated stresses is orthogonal to 2 . In other words, Equation (2) projects the directly calculated stress onto 2 , extracting its useful component. 13 For 1D solutions, N p = 1 is sufficient to obtain an accurate stress approximation. However, for 2D and 3D problems, we use multiple elements in the stress calculation domain to improve the accuracy. In this method, for the interior, side, and edge nodes ( Figures 1A-C), the stress calculation domain (a patch) is defined as the first layer of elements attached to it. The patch for the corner node ( Figure 1D) is defined by 2 element layers. For the mth element, the calculation domain is given by the element itself and all its neighbouring elements, as shown in Figures 2A-D. Clearly, the number of elements in the patch is also different whether we are dealing with an interior or a boundary element. In order to obtain the discrete FE counterparts of (1) and (2), we define the following interpolations for the enhanced stress (m) and the Lagrange multipliers (m) and (m) : where E , E ,Ē are interpolation matrices, and̂,̂,̂their corresponding coefficient vectors. Substituting these interpolations into Equations (1) and (2) results in the discrete system of equations that we solve for obtaining the unknown stress coefficients in each element as follows: In Equation (4), is the differential operator on SIP stresses to compute the divergence of the stresses and is defined as Each of the 6 SIP stress components in 3D is interpolated using 10 coefficients, resulting in an interpolation matrix E of size 6 × 60, which is valid for both hexahedral (H8) and tetrahedral (T4) elements as follows: The interpolation matrix for the Lagrange multiplier (m) is given by Finally, the Lagrange multipliers (m) ∈V are interpolated using the interpolation matrixĒ given in Equation (9). The interpolation matrix for the self-equilibrated stresses should satisfy divĒ = 0 as required by the self-equilibrating property of̄( m) . This can be easily verified by applying the operator given in Equation (5). Therefore, we have ⊺ = All matrices defined here are suitable for both tetrahedral (T4) and hexahedral (H8) elements. Furthermore, all elements in the mesh are treated the same way regardless of their position in the mesh, ie, no special treatment is needed for boundary or corner elements, even though they have a smaller stress calculation domain. Using the interpolation matrices defined in Equations (6) to (9) we solve the discrete system of Equation (4) to calculate the SIP coefficient vector̂. After obtaininĝ, the SIP stress anywhere in the element is obtained by (m) = Ê. When a nodal patch is used, the stress i at the ith node is computed as i = Ê. In this case, the stresses inside the elements are interpolated from the recovered nodal values by using the standard FE shape functions. The computational cost associated with the calculation of the recovered stresses for each element in the mesh is as small as the computational effort required to solve a linear system of equations with 60 unknowns ((n 3 ) with n = 60 if a direct solver is used). Furthermore, this cost scales linearly with the number of elements, and therefore becomes increasingly insignificant relative to the cost of obtaining the displacement solution.

EXAMPLES
In this section, we perform numerical tests on problems modeled with 3D linear FEs. No units are presented, so it is assumed that any consistent unit system can be used to interpret the results. h-convergence is performed to investigate the accuracy of the new approach. To that end, the mesh refinement adopted is not based on just splitting the FEs from the coarsest mesh but uses the exact geometry as reference for creating new elements. This is particularly important for geometries with curved surfaces so that the discretization error is minimal. We use the H 0 -norm of the stress 15 defined as where ex is the exact stress and (m) is the SIP stress in an element. The global error in the mesh is computed by adding the error contribution of every element, which is in turn computed by numerical quadrature on Gauss points. The error in directly calculated stresses (m) h can be computed in a similar way by replacing (m) by (m) h in Equation (10).

Example 1: 3D cube under uniaxial loading
This first problem is used to verify SIP in 3D. We start with a cube of unit volume subjected to a sinusoidal body force f b = 10 5 × x sin(6πx)e 1 . Homogeneous essential boundary conditions (u (m) = 0) are prescribed on its left and right faces (see Figure 3). The cube is composed of a material with Young's modulus E = 110×10 9 and Poisson's ratio = 0. The exact solution to this problem can be obtained by integrating the 1D equilibrium equation for a bar. To check for frame invariance, the block is rotated and analyzed in several rotated configurations (with the body force rotated accordingly), and the results are compared with the initial configuration. Namely, the problem is tested in 6 different orientations in 3D space: aligned with the x, y, and z axes, and then rotated 30 degrees in the xy, yz, and zx planes. We start with a mesh of 8 hexahedral elements of equal size and every new mesh refinement in the sequence is obtained by subdividing each hexahedral element by 8. For tetrahedral meshes, we start by creating tetrahedra within original hexahedral mesh, with a total of 48 tetrahedra. For the following level of refinement, we create tetrahedra within the new refined hexahedral mesh (384 tetrahedra). This mesh refinement is obtained by increasing the number of nodes on each edge of the model geometry.
The mesh size h is taken as the diameter of a sphere circumscribing any element in the FE mesh. For the calculation of the SIP stress, we consider 3-point and 8-point quadrature rules for hexahedral and tetrahedral elements, respectively. Figure 4 shows the global mesh error ||e|| H 0 as a function of the mesh size h. This result is representative for all orientations of the 3D block. We show not only that the recovered stresses are more accurate than that of the directly calculated ones but also that the former achieve superconvergence for this problem.

Example 2: infinite plate with a hole
Consider, in Figure 5, a plate with a hole subjected to a unidirectional tensile stress 0 = 1 at an infinite distance from a cylindrical hole with the radius a = 1. This plane strain elasticity problem has been extensively used in the literature as

FIGURE 6
The discretizations labeled with (A)-(E) use 8-node hexahedral elements those labeled with (F)-(J) use 4-node tetrahedral elements a benchmark test for stress recovery procedures, 3,4,9,10 and its exact solution is given by 16 the following: ) + 3a 4 2r 2 sin 4 0 0  In Figure 7, the global error given by Equation (10) is plotted as a function of mesh size h, which this time is obtained as the diameter of a sphere circumscribing the biggest element in the FE mesh. We observe a convergence rate of 1.32 for hexahedral elements and a convergence rate of 1.35 for tetrahedral elements. We notice throughout that the directly calculated ones and therefore that a much coarser mesh can be used for the same level of accuracy. Contrary to Section 3.1, we did not obtain superconvergence for this example, which we attribute to the irregularity of the mesh. A convergence rate of 1.3 was also obtained in the SPR, REP, and RCP methods. 3,4,9,10 Finally, Figures 8A-D show the x-component of the stress for both SIP and directly calculated FE stresses. From the exact solution, x = 3 at point A in Figure 5, which shows that SIP has higher accuracy for all meshes considered. In fact, the SIP stress for the coarser mesh in the figure is more accurate than the FE stress for the finer one despite of being 64 times coarser. The maximum element errors are concentrated near the points of maximum stress, as shown in Figure 9. Additionally, it can be seen in Figure 9B that SIP gives a lower error for the corner elements on edge AA ′ (see Figure 5) compared with the inner elements on that edge. This is because the stress calculation domain for the corner nodes is larger than that of the inner nodes on the boundary edge.
Starting with a coarse mesh containing 54 hexahedral elements, we again perform uniform refinement by subdividing each element into 8 equal elements. As aforementioned, we emphasize that the refinement strategy is based on the exact geometry. This results in finer meshes containing 432, 3456, and 27 648 hexahedral elements (see Figures 11A-D). A similar approach was used for the tetrahedral meshes where the coarsest mesh has 324 elements, followed by 2592, 20 736, and 165 888 elements in the sequence (see Figures 11E-H). Figure 12 summarizes the convergence results for this problem where convergence rates of 1.44 and 1.49 were obtained for hexahedral and tetrahedral meshes, respectively. As in Section 3.2, the convergence rate is not optimal because of the use of irregular meshes. However, for the meshes considered, the SIP stresses show higher accuracy than those of directly calculated from the FE solution. Figure 13 shows that the SIP stresses are much closer to the exact ones. The higher accuracy of SIP is also shown in Figure 14. Similar to Section 3.2, we notice that the element error is lower in the corner  elements A, B, and C for which we use a bigger patch of elements (see Figure 14B). Nevertheless, the SIP results are a clear improvement over directly calculated FEM stresses in terms of accuracy and convergence rate.

Example 4: spherical cavity in an infinite solid subjected to remote stress
As a final numerical example, we consider a spherical cavity of radius a = 1 in an infinite solid subjected to a remote tensile stress 0 in the z-direction. Taking advantage of the symmetry of the problem, we perform the analysis on a cube with edge length l = 4 with 1/8 of the spherical cavity in one of its corners, as shown in Figure 15. The exact solution of this 3D problem, whose derivation is given in more detail in Appendix B, is given by where = 1 + , and = −7 + 5 . The exact displacement given by Equation (15) is prescribed in the cube surfaces, and a traction-free boundary condition is prescribed on the spherical surface. Moreover, we take E = 1000 and = 0.3. As in the previous examples, we use hexahedral and tetrahedral elements to discretize the problem and we follow the same refinement strategy. The hexahedral meshes shown in Figures 16A-D contain 96, 768, 6144, and 49 152 elements, respectively. For the tetrahedral meshes in Figures 16E-H, the number of elements are 288, 2304, 18432, and 147456, respectively. Figure 17 shows the convergence plots for this problem with convergence rates of 1.40 and 1.46 for hexahedral and tetrahedral meshes, respectively. As in the previous examples, we know that though the SIP stresses convergence rates are not optimal due to irregular meshes, they are still higher than the directly calculated FEM stresses. Besides the faster convergence, the SIP stresses are also more accurate. Figure 18 also shows that SIP stresses are much closer to the exact stresses, and the errors of the SIP stresses in Figure 19 are smaller.

CONCLUSIONS
Building on the work of Payen and Bathe, 13 we have presented a methodology that can be used to recover stresses with high accuracy in 3D elasticity problems. We provide the interpolation matrices required to implement this method and verified it with the aid of several examples that have an analytical solution. The method is straightforward to implement, computationally efficient, and scalable. It is also general in the sense that it can be used for any type of low-order 3D element and mesh configuration. Moreover, no special treatment is required for corner nor boundary elements in the mesh; however, a bigger patch is used for corner nodes when computing stress at those locations. Similar to other recovery techniques, the method fails to reach superconvergence when using irregular meshes. Yet, the results obtained are, to date, the most accurate 3D recovered stresses in problems of linear elasticity; for comparison with other recovery techniques, we refer to the work of Payen and Bathe. 13 Although we only considered and tested linear elastic problems, the method could also be extended and applied to nonlinear and dynamic problems.
Regarding future research, the effectivity robustness tests proposed by Babuška et al. 17 could be performed to explore the feasibility of this method as an a posteriori error estimator. The current recovery technique could also be extended for shell and plate elements and its application could be explored, for example in topology optimization with stress constraints.

DERIVATION OF THE ENHANCED STRESS
Following the procedure outlined by Payen and Bathe, 13 we start with a three-field variational functional. Given the mth element displacements u (m) , the stresses (m) , the strains (m) u , the constitutive matrix C (m) , the body force vector f B , and the surface traction vector f S , a three-field functional can be written as where the independent fields are u (m) , (m) , and (m) . The latter denotes the Lagrangian multiplier field that enforces weakly the stress-strain relationship over the element of volume V (m) . Strains (m) u follow from the kinematic relationship through (m) u = B (m)Û(m) , whereÛ (m) is the element local degree of freedom vector and B (m) is the strain-displacement matrix. It was shown by Payen and Bathe 13 that REP, RCP, and NPF formulations can all be derived from (A1). The difference in these methods lies in the interpolations used for the Lagrange multipliers (m) . Invoking the stationarity of Equation (A1) with respect to the independent fields u (m) , (m) , and (m) yields Substituting (A4) into (A5) yields the conventional weak problem statement that is usually implemented in displacement-based finite element codes and from which the displacement can be readily computed as follows: After computing the displacements, the directly calculated FE stresses can be obtained by using proper kinematic and constitutive relationships. However, the accuracy of these stresses is an order lower than that of the displacements. Thus, instead of computing the stresses directly we use this a posteriori technique to obtain stresses with higher accuracy.

EXACT SOLUTION FOR A SPHERICAL CAVITY SUBJECTED TO REMOTE STRESS
In the following, we use a different notation for the vector coordinate x = (x, y, z) as x = (x 1 , x 2 , x 3 ). Following the procedure outlined by Bower, 18 the displacement field for the problem given in Section 3.4 can be obtained by where i (x) is a vector-valued function and (x) is a scalar function. For this problem, i (x) and (x) are given by 18 the following: respectively. The displacement field is then obtained by replacing Equations (B2) and (B3) into (B1), resulting in Equation (15). Finally, the stress tensor for this problem is given by where A = 1∕(2 R 9 ), B = 3 − 5ν, C = −2 + 5ν, D = 2 + 5ν, F = 18 + 5ν, L = −6 + 5ν, M = 19 + 5ν, N = −19 + 5 , Q = 6 + 5ν, and X ij is the sum of the squares of the ith and jth coordinates, eg, X 23 = x 2 2 + x 2 3 .