Abstract
For evaluating the stress gradient, a mathematical technique based on the stress field of lower-order 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 and Hessian matrices. The accuracy and validity of the presented method are demonstrated by two examples, especially its merit of achieving high accuracy with lower-order linear elements. This method can be conveniently introduced into the general finite element analysis programs as a postprocessing module.
1. 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–4]. To overcome the problem of the conventional elasticity theory in describing the effect of size, the couple stress [5–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 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 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 enhanced patch test of a convergence condition. In addition to the reason that only a few types of 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 element, a mixed or coupled formulation using lower-order elements () is implemented. Zhao et al. [13] proposed a mixed quadrilateral element () for the couple stress/strain gradient elasticity, in which one satisfies continuity to compute strains and the other one satisfies weak 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 continuous Hermitian shape functions for the plastic and damage multipliers and the 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 higher-order derivatives such as stress gradients also need higher-order 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 lower-order linear 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.
2. SCP Stress Gradient Recovery
2.1. Second Derivatives in FEM
Taking the isotropic normal stress , for example, the gradient vector can be written aswhere and are the components of the stress gradient, is the shape function, is the number of element nodes, and is the nodal displacement vector; , for the plane stress, , 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 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 arewhere 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 aswhere 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.
2.2. 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.
(a) Q4
(b) Q8
(c) T3
(d) T6
With the SCP recovery, the recovered stress gradient can be expressed aswhere is the polynomial vector and is the coefficient column vector to be determined.
For quadratic elements like Q8 and T6, is assumed aswhere the recovered stress gradient has the same order of displacement filed [29], and then .
The functional integral of the difference between the stress gradient calculated from the SCP recovery and the stress gradient obtained by the FEM can be expressed aswhere is the number of patch elements and is the element area.
The assumed coefficient vector can be obtained by minimizing the functional integral in (6), which yieldswhere 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 isin which
In [25, 26], utilizing the quadratic completeness polynomial, the SCP recovery for stress gradient has the order of accuracy ( is the element characteristic length). And then, it is easy to verify that the error of the stress gradient recovery is of order by the same procedure. It has been found that the recovered variable field is more accurate at node points.
3. 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.
3.1. General Formulation
As shown in Figure 1 for 2D problem, the central point is surrounded by sample points , ( is number of sample points). The vector from to is ; thus, the unit vector can be written as , where is the vector length.
The directional derivative along at the central point can be expressed aswhere is the stress result of commercial finite element software. The limit can be approximately written as Substituting (12) into (11) yields
For all sample points, the stress gradient equation can be written aswhere
In general, ; hence, the coefficient matrix is not square and (14) is overdetermined and inconsistent. An effective way to solve the equation is to employ the Moore-Penrose pseudoinverse of . The pseudoinverse can be calculated by the singular value decomposition of ; that is,where and are orthogonal matrices, 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 and , respectively. As an example, the 3D element patch is shown in Figure 2 with 8-node linear hexahedral elements.
3.2. Error Estimation
The Moore-Penrose pseudoinverse ( or 3 for 2D or 3D case) of is unique and satisfies all of the following four criteria:which yieldsFor all , there iswhere . With the first equation of (20), and thus there is a relationship. Considerin which denotes the inner product. With this relationship, we can haveorwhich 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 asin which and .
Similarly, it can be verified that if , , thenIt 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 must be the Moore-Penrose pseudoinverse . 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, is the weight (element area or volume), and 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.
4. Numerical Example and Discussions
4.1. Example 1: 2D Problem
For existing exact solution and exhibiting strong stress concentration phenomenon, the infinite plate with a center hole under unidirectional tensile stress as shown in Figure 3 becomes a classic and popular two-dimensional example used in previous papers [13–15, 24, 29].
(a)
(b)
In this case, the normal stress is dominant and the exact solution is given aswhere and are the polar coordinates.
Due to the symmetry in the geometry, only a quarter area is modeled using quadrilateral and triangular elements, such as T3, Q4, T6, and Q8 elements. Three different meshes are shown in Figure 4, including the numbers of nodes and elements.
(a) Quadrilateral coarse mesh 400 elements: 441 nodes (Q4), 1281 nodes (Q8)
(b) Quadrilateral fine mesh 1600 elements: 1681 nodes (Q4), 4961 nodes (Q8)
(c) Quadrilateral double-fine mesh 6400 elements: 6561 nodes (Q4), 19521 nodes (Q8)
(d) Triangle coarse mesh 564 elements: 323 nodes (T3), 1209 nodes (T6)
(e) Triangle fine mesh 2220 elements: 1191 nodes (T3), 4601 nodes (T6)
(f) Triangle double-fine mesh 9040 elements: 4681 nodes (T3), 18401 nodes (T6)
Firstly, a checking point around is used to compare the results of different methods as shown in Figure 4. The exact locations are , , and for coarse, fine, and double-fine meshes, respectively. Results from the proposed stress field gradient (SFG) analysis technique and the SCP recovery method with various element types and meshes are compared in Table 1. The average error is the root sum square of -direction gradient () error and -direction gradient () error. The accuracy of these methods is all improved with the mesh refinement no matter which element is employed, especially for quadrilateral elements. Except fine mesh T6 element (F-T6) only, the presented method is always more accurate than the SCP recovery technique. For lower-order linear elements Q4, it is clear that the SFG technique can still provide a better result than the SCP method in all the mesh types. On the other hand, using T3 element, results of the SFG get closer to those of the SCP with the mesh refinement, and finally a higher accurate stress gradient can be obtained in the double-fine mesh. In the stress field gradient (SFG) analysis, higher-order elements do not always correspond to higher accuracy because of the treatment of minimum norm least squares solution. In other words, a lower-order element, which generally means low computation cost, can provide almost the same accuracy as a higher-order element. This is particularly important to practical applications.
Furthermore, results from the two methods along two checking lines, and , are shown in Figures 5 and 6, respectively. Both figures indicate that they have the uniform accuracy in most areas, including the SFG procedure employing linear T3 and Q4 elements. Although the results at the edge of the hole have the lowest precision for all cases, the accuracy is very high once not at the edge. Specifically, at the edge of the hole, the SFG method using linear elements provides the worst stress gradient for line (Line-) as well as the SCP method using quadratic Q8 element for line (Line-). Utilizing quadrilateral elements, the two approaches have almost identical stability for each mesh. On the other hand, the SCP recovery method with triangular elements works better in the stability and convergence comparing with the SFG method.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
4.2. Example 2: 3D Problem
A three-dimensional cantilever beam with a groove under vertical end load is modeled using 8-node hexahedral elements, and its geometry, loading condition, and mesh are shown in Figure 7. All DOFs on the surface are fixed.
(a)
(b)
Based on the software platform ANSYS, the flowchart of postprocessing with the proposed method is shown in Figure 8. The code of the stress field gradient analysis is developed by Matlab language and the contour figures are generated using software Tecplot. The general finite element analysis program (FEAP), which provides the preliminary analysis of the stress field, is not just limited to ANSYS, and 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 mm and mm, in which the total gradient . 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.
(a)
(b)
(c)
(d)
(e)
(f)
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 to . 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 to 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.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
5. 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 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.