Stress Field Gradient Analysis Technique Using Lower-Order C 0 Elements

For evaluating the stress gradient, a mathematical technique based on the stress field of lower-order C0 elements is developed in this paper. With nodal stress results and location information, an overdetermined and inconsistent equation of stress gradient is established and the minimum norm least squares solution is obtained by the Moore-Penrose pseudoinverse. This technique can be applied to any element type in comparison with the superconvergent patch (SCP) recovery for the stress gradient, which requires the quadratic elements at least and has to invert the Jacobi andHessianmatrices.The accuracy and validity of the presentedmethod are demonstrated by two examples, especially itsmerit of achieving high accuracy with lower-order linearC elements.Thismethod can be conveniently introduced into the general finite element analysis programs as a postprocessing module.


Introduction
Effects and evaluations of the stress gradient have always been important to the design of an engineering structure.As nowadays the stress is usually calculated with a finite element method (FEM), the calculation technique of the stress gradient has been the effort of many research works.
The size effect on the stress/strain gradient has already been manifested in many experiments [1][2][3][4].To overcome the problem of the conventional elasticity theory in describing the effect of size, the couple stress [5][6][7] and strain gradient plasticity [8,9] theories have been developed and proved effective in introducing strain gradient terms into constitutive equations and yield functions.However, according to these theories, the element should be at least  1 continuity.This complicates the element construction.Dasgupta and Sengupta [10] developed an 18-DOF (degree-of-freedom) triangular plate bending element, whose displacement function is fifth-order polynomial in the area coordinates.Zervos et al. [11] proposed a  1 continuity elastoplasticity triangular element with assumption that the stress increment consists of both the strain increment and its Laplacian.Chen and Li [12] put forward a 24-DOF quadrilateral spline element for the couple stress/strain gradient elasticity, which can pass the  0-1 enhanced patch test of a convergence condition.In addition to the reason that only a few types of  1 elements are successful and effective, the limitation in the element shape and the number of DOFs all impede the practical applications.
In order to circumvent problems of the  1 element, a mixed or coupled formulation using lower-order elements ( 0 ) is implemented.Zhao et al. [13] proposed a mixed quadrilateral element (CQ12 + RDKQ) for the couple stress/strain gradient elasticity, in which one satisfies  0 continuity to compute strains and the other one satisfies weak  1 continuity to calculate strain gradients.Amanatidou and Aravas [14] developed an alternative mixed finite element formulation and studied the Mindlin [15] strain gradient elasticity theories in detail.The  1 continuous Hermitian shape functions for the plastic and damage multipliers and the  0 continuous interpolation functions for displacement field were employed by Dorgan and Voyiadjis [16].For reasons such as the additional and inevitable computational cost and only a few available element types, the mixed elements are inconvenient to practical applications.
Another motivation for the stress gradient analysis is the fatigue life calculation and the fatigue behavior prediction.Neuber [17] developed a widely adopted method in predicting the notch fatigue life; nevertheless, the method only focuses on the stress of dangerous point but without considering the influence of stress gradient and multiaxial effects.Even if the stresses of dangerous points for various structures are the same, the fatigue life and behavior are different because of different stress gradient distribution [18].Taking into account the stress gradient for assessing the fatigue strength [19], the stress field intensity method [20] and the advanced volumetric method [21] are two of the most popular techniques.Although the stress field intensity method is simple, it needs experiment to provide the size of notch damaged zone.The advanced volumetric method requires the accurate calculation of stress gradients near notch roots; however, it is usually difficult to guarantee the calculation accuracy.
The technique of finite element postprocessing is another approach to calculating continuous stress and its gradient.The superconvergent patch (SCP) recovery is known as the most effective technique [22], which can obtain higher accuracy than any other points at the Gauss integration point within an element.Utilizing the patches of elements in FEM, Zienkiewicz and Zhu [23,24] obtained continuous stress field with data at the superconvergent points.A detailed analysis of the SCP recovery technique and error estimation was presented in [25,26].Since then, Wiberg and Abdulwahab [27], Wiberg et al. [28], and Park and Shin [29] improved the SCP approach using the method of least square for fitting the residuals of the equilibrium equation and applying the principle of virtual work.Recently, the SCP technique has been extended to recover the higher-order derivatives.Using Q8 (8-node quadratic, 9-point quadrature) quadrilateral element, T6 (6-node quadratic, 4-point quadrature) triangular element, and the mixed mesh of Q8 and T6, Gan and Akin [30] developed an element patch based superconvergent second derivative recovery technique for a lower-order strain gradient plasticity model.Bank et al. [31] proposed a superconvergent derivative recovery scheme with the Lagrange triangular element and the superconvergence estimation of the th derivatives.Nevertheless, the higherorder derivatives such as stress gradients also need higherorder element shape functions, which lead to additional DOFs and calculation cost, especially for complex structures or large systems.
Therefore, the effort is still required to enhance techniques for calculating stress gradient in such aspects as reducing computation cost, developing more element types, and providing codes for engineering applications.In the rest of this paper, a mathematical technique based on the stress field results of general finite element analysis programs (FEAP) is proposed to calculate the stress gradient.Because it works with the stress field results, this method can be considered as a stress field gradient (SFG) analysis technique.As no transformation from parametric coordinate to physical coordinate is required, this technique does not involve Jacobi and Hessian matrices and hence does not need to calculate their inverses.With the location and stress information of the central node and its surrounding sample nodes, an overdetermined and inconsistent equation of the stress gradient at the central node is established.
The Moore-Penrose pseudoinverse, which is computed with the singular value decomposition, is employed to obtain the minimum norm least squares solution [32,33] of the stress gradient.A popular 2D example demonstrates its satisfactory accuracy in comparison with the SCP recovery.The presented method especially can also achieve high accuracy with lowerorder linear  0 elements.Another 3D example shows that the method can be easily introduced into general FEAP (such as ANSYS, PATRAN, and ABUQUS) as a postprocessing module.
This paper is organized as follows.First, the superconvergent patch recovery method is briefly introduced in Section 2.Then, the stress filed gradient analysis technique and the error estimation are presented in Section 3. In Section 4, two numerical examples are employed to evaluate the performance of the presented method.The last section of the paper is the conclusions.

Second Derivatives in FEM.
Taking the isotropic normal stress   , for example, the gradient vector G(  ) can be written as where   and   are the components of the stress gradient,   is the shape function,  is the number of element nodes, and (  , V  ) is the nodal displacement vector;  0 = ,  0 =  for the plane stress,  0 = /(1 −  2 ),  0 = /(1 − ) for the plane strain;  and  are Young's modulus and Poisson's ratio, respectively.Firstly, the second derivatives of the shape functions are necessary for evaluating the stress gradient, but they are all equal to zero for lower-order linear  0 elements, such as Q4 (4-node quadrilateral element) and T3 (3-node triangular element).
Assuming using parametric elements, the first derivative transformations from the parametric coordinate to the physical coordinate are where  and  are the parametric coordinates, and the matrix is the so-called Jacobi matrix.The relationship between the second parametric derivatives and the second physical derivatives can be written as where the second square matrix on the right hand side is the Hessian matrix.From (3), the physical second derivatives can be expressed by the first and second parametric derivatives and the inverse of the Jacobi and Hessian matrices.

SCP Recovery and Error Estimation.
The Gauss integration points are also the superconvergent points for the stress gradient, and the superconvergent patches can be built as shown in Figures 1(b) and 1(d).The element shape functions must satisfy quadratic at least in the SCP stress gradient recovery because of the requirement on the existence of the second derivatives.
With the SCP recovery, the recovered stress gradient  * , can be expressed as where P is the polynomial vector and a is the coefficient column vector to be determined.For quadratic elements like Q8 and T6, P is assumed as where the recovered stress gradient has the same order of displacement filed [29], and then a The functional integral of the difference between the stress gradient calculated from the SCP recovery  * , and the stress gradient obtained by the FEM g, can be expressed as where  is the number of patch elements and   is the element area.
The assumed coefficient vector a can be obtained by minimizing the functional integral in (6), which yields where superscript  denotes the transpose, and both sides of ( 7) can be calculated by the Gauss integral: where  =  ⋅   ,   is the number of Gauss points per element,   is the weight coefficient of each integration point, and (  ,   ) is the global coordinates of each Gauss point.Finally, there is In [25,26], utilizing the quadratic completeness polynomial, the SCP recovery for stress gradient has the order of accuracy (ℎ 3 ) (ℎ is the element characteristic length).And then, it is easy to verify that the error of the stress gradient recovery is of order (ℎ 2 ) by the same procedure.It has been found that the recovered variable field is more accurate at node points.

Stress Field Gradient Analysis
It is known that the SCP recovery for the stress gradient needs quadratic elements at least, and the shape functions should be determined beforehand to compute the second derivatives of the displacement field.All of this leads to an additional calculation cost and difficulty, especially to those engineers who are used to work with commercial finite element software, and providing the preliminary analysis by lower-order linear elements, where the SCP technique cannot be implemented yet.In this section, a simple and convenient mathematical technique for calculating the stress gradient will be developed.This method is suitable for any element types and only needs the information of nodal location and stress.

General Formulation.
As shown in Figure 1 for 2D problem, the central point   ( 0 ,  0 ) is surrounded by sample points   (  ,   ),  = 1, 2, . . .,   (  is number of sample points).The vector from   to   is n  = (  −  0 ,   −  0 ); thus, the unit vector can be written as Substituting ( 12) into (11) yields For all sample points, the stress gradient equation can be written as where In general,   > 2; hence, the coefficient matrix A is not square and ( 14) is overdetermined and inconsistent.An effective way to solve the equation is to employ the Moore-Penrose pseudoinverse of A   ×2 .The pseudoinverse M 2×  can be calculated by the singular value decomposition of A   ×2 ; that is, where U ∈    ×  and V ∈  2×2 are orthogonal matrices, Σ ∈    ×2 is a diagonal matrix with nonnegative real numbers, and the diagonal entries are known as the singular values.Then, where Σ † is the transpose of Σ and its diagonal entries are the inverses of those nonzero diagonal elements of Σ.Finally, the stress gradient at the central point   can be expressed by As the same procedure as 2D problem, the stress gradient equation of 3D problem can also be written as (14), in which the direction vector and stress gradient vector are n  = (  −  0 ,   −  0 ,   −  0 ) and G   (σ  ) = (σ  /, σ  /, σ  /), respectively.As an example, the 3D element patch is shown in Figure 2 with 8-node linear hexahedral elements.

Error Estimation.
The Moore-Penrose pseudoinverse M  0 ×  ( 0 = 2 or 3 for 2D or 3D case) of A   × 0 ∈    × 0 is unique and satisfies all of the following four criteria: which yields or which means that Moore-Penrose pseudoinverse provides the linear least squares solution.
The least squares solution is usually not unique, and the general solution can be expressed as It indicates this solution has the minimum norm.
Apparently, the constraint conditions of ( 25) and ( 26) are the same as the four criteria listed in (19), so that the matrix M must be the Moore-Penrose pseudoinverse M. In this way, it is shown that the Moore-Penrose pseudoinverse can provide the unique minimum norm least squares solution of (14).
The error of ( 14) consists of two parts, one is the limit approximation of order (ℎ) in (12) and the other one is the stress field σ .The stress calculated with FEM is not continuous across elements, and the general FEAP commonly calculates the continuous stress by the method of nodal stress weighted average: in which  is the number of surrounding elements for one node,   ( = 1, 2, . . ., ) is the weight (element area or volume), and   ( = 1, 2, . . ., ) is the stress calculated by the node displacements of each element.The stress field σ evaluated by (27) has the order of accuracy (ℎ) at least.Therefore, the order of the minimum norm least squares solution of ( 14) is of (ℎ) at least.Moreover, it should not be ignored that the least squares processing has the strong ability to improve the accuracy.The following example shows that the accuracy of the suggested method, even employing the lower-order linear elements, is equal to or higher than the SCP recovery using quadratic elements.the proposed method can also be applied to other software platforms as PATRAN and ABUQUS, and so forth.Details of the stress gradients are shown in Figure 9 for ℎ = 4mm and  = 5mm, in which the total gradient (σ) Total = 2 √ (σ/) 2 + (σ/) 2 + (σ/) 2 .It can be seen that the stress gradient in Figure 9(b) is more powerful in characterizing and identifying the stress concentration than the stress itself in Figure 9(a).The gradient slices show that the stress gradient maximum appears around the surfaces and corners of the groove, which are the sensitive or dangerous areas to the fatigue failure and the plastic deformation.
Results of different groove widths and depths obtained by the SFG technique are compared in Figure 10 and Table 2.The high-gradient zone moves to the fixed end with the increase of the width, and the maximum of the stress gradient is reduced from 4.15 + 10 to 2.52 + 10.The distribution of stress gradient becomes more and more uniform and smooth with the width increasing, which is beneficial to prevent structural failure.As the depth increases, the total stress gradient goes up from 4.0531 + 10 to 2.3719 + 11 and the high-gradient zone in the bottom of the groove becomes more and more concentrated.More details are shown in Figure 11, in which the middle lines cross the bottom surfaces of the grooves with different widths and depths along the -direction.

Conclusions
An FEAP-based mathematical technique is developed for accurately extracting stress gradient.Comparing with the SCP recovery method, which needs the quadratic elements at least and must invert the Jacobi and Hessian matrices, this method only requires nodal stress results as well as location information and can be implemented to any element types.
The classic plane example shows that the suggested method can achieve more accurate stress gradient than the SCP recovery technique with the same elements.In addition, its performance in calculating the stress gradient with the linear  0 elements is proved better than the SCP recovery method.It is also observed that with this method the quadrilateral elements can provide better stability than the triangular elements.Besides, based on the proposed technique, a   three-dimensional example has presented the postprocessing analysis with commercial finite element software, and the method should be very convenient to be introduced into FEAP as a code module.

Figure 1 :
Figure 1: Construction of 2D patches with quadrilateral and triangular elements.

Figure 3 :Figure 4 :
Figure 3: An infinite plate with a central hole subjected to a remotely uniform tension: (a) theory model, (b) finite element model.

Figure 6 :
Figure 6: Stress gradients of sample lines with triangular elements.

Figure 11 :
Figure 11: Total gradients of middle lines for different groove widths and depths.

Table 1 :
Stress gradients at point  with different meshes and methods.

Table 2 :
Maximal Von-Mises stresses and gradients of different groove widths and depths.