Interpolation Type Stress Recovery Technique Based Error Estimator for Elasticity Problems

Finite Element Method (FEM) is a powerful technique for simulating real life problems, but it is only able to provide an approximated solution that necessitates reliable error control of computed solutions. The various error estimation techniques have been developing for last several years. A critical review of different error estimation techniques to get the practical finite element solution of linear, non-linear and transient problems is given in [1]. Error estimators can either be classified based on its convergence through the global effectivity index, which is defined as ratio of the estimated error to the true error or can be classified based on the procedures to obtain the estimates. The classification of error estimators based on convergence are asymptotically exact, asymptotically (upper) bounded and asymptotically not bounded. The classification of error estimators as per the estimation procedure are the residualbased error estimators [2], constitutive relation error (CRE) based error estimators [3] and the recovery-type error estimators [4]. The explicit-type recovery error estimator in energy norm is proposed in [5] for the linear elasticity problem using smooth solution. In [6], a posteriori error estimate has been developed and improved convergence is shown using non-coinciding meshes for problems in linear elasticity. An effective error estimator presented in [7] is based on continuous estimated stress field that is achieved by interpolating from nodal stresses over the element with the shape functions. In [8], it is shown that the technique of duality error majorants can serve as an effective local indicator of modeling errors of nonlinear problems. A posteriori error estimate based on an equilibrated stress reconstruction that is obtained from mixed finite element solutions of local Neumann linear elasticity problems is presented in [9]. The quality of recovery based error estimation depends on the way to obtain the smoothed or postprocessed continuous stress field. The recovery techniques are based on the least square fitting of velocity (or the displacement) field or their derivatives (stress field) by a higher order polynomial over a patch of elements or nodes. The various authors have proposed procedures for recovery of postprocessed stress field. The super-convergent patch recovery scheme is proposed in [10] in which postprocessed stresses are obtained by interpolating from a stress surface fitted to the superconvergent stress points surrounding the node of interest. In [11], a patch recovery scheme is employed in which recovery is performed for all components simultaneously. The coupling of the stress components is made through the equations of equilibrium. In [12], computational aspect of Goal-oriented error estimates (GOEE) based on enhanced Superconvergent Patch Recovery is presented in controlling the local error in quantities of interest (QoI). A method of extrapolation based on patch recovery is proposed in [13] for obtaining continuous stress field in a local manner. A recovery technique based upon the least square fitting of velocity field over an element patch is proposed in [14]. In [15], a moving least squares (MLS) recovery-based procedure to obtain postprocessed smoothed stresses field is presented in which the continuity of the recovered field is provided by the shape functions of the underlying mesh. Investigations are reported in [16] for getting improved recovery of stress field using domain decomposition method in heterogeneous structures. In [17], an improvement of the SPR technique, called SPR-C technique (Constrained SPR), is presented and uses the appropriate constraint equations in order to obtain stress interpolation polynomials in the patch. Investigators have also shown interest to study the effectivity of error estimators and propose different techniques to enhance the effectivity indices of error estimators. In [18], the formulations for recovery based error estimators, relying upon the Zienckiewicz-Zhu recovered gradient, are presented and their results in terms of effectivity indices are compared. In [19]., it is analysed the effectivity of the Zienkiewicz-Zhu error estimator considering the depending parameters such as smoothing procedure, the type of refinement (uniform or adaptive), the type of triangular element (linear or quadratic), and the number of integration points used in the numerical integration of the error estimation. The upper and lower bounds for the effectivity index to judge the quality of error estimator for linear and unstructured triangular meshes is computed in [20]. The present study aimed to present a recovery of higher derivative of field variable utilizing moving least squares (MLS) technique over a patch of mesh independent nodes for formulating a posteriori error estimator. The moving least squares (MLS) procedure has inherent completeness, robustness and continuity [21] and, proved to interpolate random data with higher accuracy. The performance of proposed error estimator is compared on benchmark elastic problems with ZZ error estimator [10], which evaluates the error using the recovery technique based upon the least square fitting of stress field over a mesh dependent nodes patch, in term of effectivity and rate of convergence.


Introduction
Finite Element Method (FEM) is a powerful technique for simulating real life problems, but it is only able to provide an approximated solution that necessitates reliable error control of computed solutions.The various error estimation techniques have been developing for last several years.A critical review of different error estimation techniques to get the practical finite element solution of linear, non-linear and transient problems is given in [1].Error estimators can either be classified based on its convergence through the global effectivity index, which is defined as ratio of the estimated error to the true error or can be classified based on the procedures to obtain the estimates.The classification of error estimators based on convergence are asymptotically exact, asymptotically (upper) bounded and asymptotically not bounded.The classification of error estimators as per the estimation procedure are the residualbased error estimators [2], constitutive relation error (CRE) based error estimators [3] and the recovery-type error estimators [4].The explicit-type recovery error estimator in energy norm is proposed in [5] for the linear elasticity problem using smooth solution.In [6], a posteriori error estimate has been developed and improved convergence is shown using non-coinciding meshes for problems in linear elasticity.An effective error estimator presented in [7] is based on continuous estimated stress field that is achieved by interpolating from nodal stresses over the element with the shape functions.In [8], it is shown that the technique of duality error majorants can serve as an effective local indicator of modeling errors of nonlinear problems.A posteriori error estimate based on an equilibrated stress reconstruction that is obtained from mixed finite element solutions of local Neumann linear elasticity problems is presented in [9].The quality of recovery based error estimation depends on the way to obtain the smoothed or post-processed continuous stress field.The recovery techniques are based on the least square fitting of velocity (or the displacement) field or their derivatives (stress field) by a higher order polynomial over a patch of elements or nodes.The various authors have proposed procedures for recovery of postprocessed stress field.The super-convergent patch recovery scheme is proposed in [10] in which post-processed stresses are obtained by interpolating from a stress surface fitted to the superconvergent stress points surrounding the node of interest.In [11], a patch recovery scheme is employed in which recovery is performed for all components simultaneously.The coupling of the stress components is made through the equations of equilibrium.In [12], computational aspect of Goal-oriented error estimates (GOEE) based on enhanced Superconvergent Patch Recovery is presented in controlling the local error in quantities of interest (QoI).A method of extrapolation based on patch recovery is proposed in [13] for obtaining continuous stress field in a local manner.A recovery technique based upon the least square fitting of velocity field over an element patch is proposed in [14].In [15], a moving least squares (MLS) recovery-based procedure to obtain postprocessed smoothed stresses field is presented in which the continuity of the recovered field is provided by the shape functions of the underlying mesh.Investigations are reported in [16] for getting improved recovery of stress field using domain decomposition method in heterogeneous structures.In [17], an improvement of the SPR technique, called SPR-C technique (Constrained SPR), is presented and uses the appropriate constraint equations in order to obtain stress interpolation polynomials in the patch.
Investigators have also shown interest to study the effectivity of error estimators and propose different techniques to enhance the effectivity indices of error estimators.In [18], the formulations for recovery based error estimators, relying upon the Zienckiewicz-Zhu recovered gradient, are presented and their results in terms of effectivity indices are compared.In [19]., it is analysed the effectivity of the Zienkiewicz-Zhu error estimator considering the depending parameters such as smoothing procedure, the type of refinement (uniform or adaptive), the type of triangular element (linear or quadratic), and the number of integration points used in the numerical integration of the error estimation.The upper and lower bounds for the effectivity index to judge the quality of error estimator for linear and unstructured triangular meshes is computed in [20].The present study aimed to present a recovery of higher derivative of field variable utilizing moving least squares (MLS) technique over a patch of mesh independent nodes for formulating a posteriori error estimator.The moving least squares (MLS) procedure has inherent completeness, robustness and continuity [21] and, proved to interpolate random data with higher accuracy.The performance of proposed error estimator is compared on benchmark elastic problems with ZZ error estimator [10], which evaluates the error using the recovery technique based upon the least square fitting of stress field over a mesh dependent nodes patch, in term of effectivity and rate of convergence.

Problem formulation
Let us consider the two-dimesional linear elastic problem.It is to find a stress field  and a unknown displacement field u, satisfyimg over a domain  that is bounded by . Static equilibrium: where: P is the force vector and T L is the derivative operator.The strain vector () can be calculated as follows: Boundary conditions: where: D contains elasticity coefficient of linear isotropic constitutive law that relates the stress to strain.
According to the Finite Element Method, the equilibrium equations are numerically solved in a domain that is discretized into several finite elements.The displacements of any point within an element are calculated based on the following equation: where:  is the nodal displacement matrix, u is the displacement of a certain point within an element and N is the matrix of the interpolation functions, also known as shape functions.Since the strains can be related to the nodal displacements in the following formula: In which, B is the strain interpolation matrix, and the stresses within an element can also be related to nodal forces, the relationship between nodal forces and nodal displacements can be described as follows: In the above equation, f is the nodal force vector and D is known as the stiffness matrix.The global equa- tions are obtained from assembly of elemental equations.The K and f matrices are found from following equations.
where b is the body force per unit volume.

Error estimation
The error in computed stress, *  e is defined as the difference between the exact or (recovered) values of stress,  , and respective computed values, h  , i.e.
The error can be evaluated in any appropriate norms.Since the finite element solution minimizes the error in the energy norm, the magnitude of the error in energy norm is a good measure of the overall quality of solution.The integral measure of the error in energy norm may be defined as follows: The standard measure of the quality of error estimators i.e., effectivity index,  , is defined as ratio of esti- mated error and true Error.An estimator is asymptotically exact for a particular problem if its effectivity index converges to one when the mesh size approaches to zero.
, e e es =  (13) where: e denotes the exact error in energy norm, and es e represents the evaluated error estimate.

Moving Least Squares based recovery technique
In Moving Least Square (MLS) technique, three components are used to express the function , a weight function ) (x w associated to each node, a basis ) (x P , usually consisting of a polynomial, and a set of coefficient, ) (x a which are functions of the coordinates.Let the nodes be defined by 1 x ... n x approximation can be defined as: where: ) (x  is shape function and m is the total number of terms in the basis.In the present study, the basis function is ) (x P has adopted m as six for linear elements.
The vector of coefficients ) (x a can be obtained by minimizing a weighted residual as follows: where: n is the number of nodes i and is a weight function in two-dimension associated to each node (domain of influence of that node) which is usually built in such a way that it takes a unit value in the vicinity of the point where the function and its derivatives are to be computed and vanishes outside a region where: By putting MLS shape function into the wellknown form of shape function equation, we get: where: The derivatives of weight function with respect to i x is computed easily using the chain rule of differentiation.

Adaptive Mesh Improvement Strategy
The refinement strategy depends on the nature of the accuracy criterion that is to be satisfied.A very common requirement is to specify the achievement of a certain minimum percentage error in the energy norm.The error in individual elements in relation to the global error helps us to decide which portions of the mesh need improvement.In the present study, the energy norm of the error is adopted for assessing of the quality of the solution and hrefinement scheme is employed for improving the mesh.The detail of the adopted mesh improvement scheme is given in [22].
The permissible global error is given by: where: p is the order of the approximating polynomial and old h is the old size of the ith element.

Numerical example
The performance and reliability of the error estimators are investigated by2-Dbenchmark example and plate problem in elasticity for which exact solution is available.The linear and higher order elements have been used for the discretization of the domain.
The constants α and β are given as: where: E and  are modulus of elasticity and Poisson's ratio respectively with a value of 1.0 N/mm 2 and 0.3.
The test problem domain is discretized with twodimensional three and six node triangular elements.Uniform subdivision of elements is used for refinement process.The adaptive mesh improvement strategy is also used to study the local error distribution with desired solution accuracy having 3 and 8% global error.The refinement is guided through the error predicted ZZ and proposed interpolation type error estimators.The initial structured mesh with three node and six node triangular elements is shown in Fig. 1 The variation of error in energy norm i.e. convergence rate and effectivity of error estimators for the test problem sequence of refinement with3-node triangular elements are given in Tables 1 to 4.The convergence rate and effectivity of error estimatorswith6-node triangular elements are given in Tables 5 and 6.
Table 1 Convergence rate for test problem using 3-node Triangular Mesh (Structured mesh) This problem involves square portion (4x4) from an infinite elastic plate having a rigid central circular inclusion with a radius (a) of 1 unit and it is the typical problem of a stress concentration in an infinite membrane.The exact stress field for the problem are given in following equations [23]:

Discussion
The finite element analysis numerical results to compare the performance of the proposed error estimation technique and Zienkiewicz-Zhu (ZZ) error estimation technique on two benchmark problems has been presented for linear (3-node) and higher order elements (6-node) structured as well as unstructured triangular mesh in Table 1-8 and Figure 1-6.For both 3-node structured and unstructured mesh, and 6-node unstructured mesh, the convergence obtained with the help of proposed mesh independent node patch (element free Galerkin (EFG) approach) based recovery scheme is found to be better than that for the mesh dependent node patch based ZZ super convergent recovery scheme.Besides, the order of error in comparison to ZZ scheme is smaller, thereby indicating higher effectivity of the propose error estimation scheme.
The effectivity of error estimation is also compared for 3-node and 6-node structured and unstructured meshes through Tables 2, 4, 6 and 8, and seen that element free Galerkin (EFG) approach based recovery technique perform better than the ZZ super convergent recovery technique specially in problems with singularity.In ZZ recovery scheme, especially linear elements, there is a possibility of less accurate recovery for boundary nodes since element dependency provide lesser number of nodes.Such difficulty does not arises with the propose element free recovery scheme.
The adaptively refined meshes plot given in Figs. 3 to 6, obtain through the ZZ and proposed EFG error estimators using triangular elements give a picture of distribution of local error in the finite element solution domain.The meshes depict that the structured mesh adaptively refined with finer element in areas of high local error to get a uniform accuracy throughout the domain.It also concludes from the mesh plots that the initial domain sub division greatly affects the quality of the refined meshing.The coarser initial mesh will result in a uniform error in the domain overlooking the high gradients zones.It also infers from the mesh plots that the propose error estimator effectively and efficiently compute the error in energy norm of the recovered solution both at local and global levels.
Table 9 and Table 12 present actual global error of the FE solution and global error respectively of structured and unstructured mesh predicted by the proposed and ZZ error estimators.The number of element in the mesh and degree of freedom (DOF) is also given.Table 10 and Table 13 depict the number of element and DOF in adaptively refined triangular structured and unstructured mesh respectively at prescribed target error of 3% and 8% using ZZ error estimator.11 and Table 14 depict the number of element and DOF in adaptively refined structured and unstructured mesh at prescribed target error of 3% and 8% using proposed interpolation based error estimator.It is clear from the tables that triangular mesh and mesh improvement as per local error prediction through zz error estimator needs more elements to achieve the target accuracy as compared to proposed estimator.The tables show that more or less similar number of elements is required in uniform mesh or adaptive mesh to get the lower target accuracy.However, required number of elements increased to get the higher accuracy of solution.

Conclusions
In the present work, effective interpolation type gradient recovery based estimation technique for error in finite element solution is proposed.The recovery is based on Moving Least Squares approximation (element free Galerkin approach) of the displacement field or their derivatives over a patch of nodes in a circular boundary.The performance of proposed error estimator has been compared with that of the ZZ error estimator by applying it to benchmark elastic problems using structured as well as unstructured meshes.The proposed mesh independent node patch based recovery scheme is found to be better than that for the mesh dependent node patch based ZZ super convergent recovery scheme.The adaptive mesh improvement strategy, based on guidance of the local error predicted by ZZ and proposed error estimators, has also been used to study the error distribution in the domain.It conclude from the results that propose error estimator effectively and efficiently compute the errors both at local and global levels in FE domains.

INTERPOLATION TYPE STRESS RECOVERY TECHNIQUE BASED ERROR ESTIMATOR FOR ELASTICITY PROBLEMS
S u m m a r y A finite element method coupled with error estimation has gained considerable prominence in industry.However, effective and reliable error control of finite element solution is always a challenging task particularly for incompressible and large deformations problems.Effective interpolation type gradient recovery based error estimation procedure is proposed in the present study.The recovery is based on Moving Least Squares approximation (mesh free approach) of the displacement field or their derivatives by a higher order polynomial over a patch of nodes in a circular boundary.The performance of error estimation scheme in terms of its effectivity and convergence has been compared with that of Zienkiewicz-Zhu (ZZ) super-convergent recovery scheme by applying the scheme to benchmark elastic problems.Error estimators compute the error in energy norm of the recovered solution both at local and global levels.The adaptive meshing based on guidance of the local error predicted by ZZ and proposed interpolation type error estimators, is also used to study the error distribution in the domain.The proposed mesh independent node patch based recovery scheme is found to be better than that for the mesh dependent node patch based ZZ super convergent recovery scheme

id
 surrounding the point i x .The support of the shape function ) ( 1 x  is equivalent to the support of the weight function.The following cubic spline weight function with circular domain of influence is considered in the present study: is the size of influence domain of the point i x .Minimization of weighted residual leads to:

  1 ,
error percentage, e is the global strain energy error, k is a factor lying between 1.0 to 1.5 to prevent oscillation.The following relation gives the permissible error in the ith element:where: N are total number of element.The so-called ele- ment refinement parameter i  guides the refinement: refinement is needed.The new element size new h is found with the help of the following equation:

6. 1 .
1×1 Square test problem The benchmark examples considers an infinite domain problem where it is extracted a 1×1 square domain, centred at the origin of coordinates with body loads ) .Mesh sequences with 3and6 node triangular elements are considered for the analyses.The uniform and unstructured refinements has been used.The exact order polynomials, are given in Equation 27 to 30.The corresponding Neumann boundary condition is imposed.

Fig. 1
Fig. 1 Initial domain discretization with triangular elements quarter of square domain need to be modeled because of symmetry of plate problem.The finite portion of the plate discretized with triangular mesh is shown in Fig.2.The boundary conditions are as follows.Along the circular arc, both displacement components are zero.Along the symmetry line, the normal displacement component and shear stress are zero.Static boundaries conditions are imposed from the traction computed from the above equations.

Table 3 Convergence
rate for test problem using 3-node Triangular Element (Unstructured mesh)

Table 4 Effectivity
for test problem using 3-node Triangular Element (Un-structured mesh)

Table 7
Error in energy norm for Plate problem using 3-node Triangular Element (Unstructured mesh)

Table 9
Number of element and DOF in triangular element structured mesh along with actual and predicted global errors

Table 11 Number
of element and DOF in adaptively refined mesh (Proposed estimator)

Table 12 Number
of element and DOF in unstructured mesh along with actual and predicted global errors Fig.5Adaptive mesh with triangular elements at different initial unstructured mesh (ZZ error estimator, 3% and 8% target error) Table14Number of element and DOF in adaptively refined mesh (Proposed estimator, 3% and 8% target error) Fig.6Adaptive meshing at different initial unstructured discretization (Proposed estimator, 3% and 8% error)