# Stress Field Gradient Analysis Technique Using Lower-OrderC0 Elements

Authors: Jianwei Xing; Gangtie Zheng
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

For evaluating the stress gradient, a mathematical technique based on the stress field of lower-orderC0 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 C0 elements. This method can be conveniently introduced into the general finite element analysis programs as a postprocessing module.


---

## Body

## 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 C1 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 C1 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 C0-1 enhanced patch test of a convergence condition. In addition to the reason that only a few types of C1 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 theC1 element, a mixed or coupled formulation using lower-order elements (C0) is implemented. Zhao et al. [13] proposed a mixed quadrilateral element (CQ12+RDKQ) for the couple stress/strain gradient elasticity, in which one satisfies C0 continuity to compute strains and the other one satisfies weak C1 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 C1 continuous Hermitian shape functions for the plastic and damage multipliers and the C0 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 pth 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 C0 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 Section2. 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σx, for example, the gradient vector G(σx) can be written as(1)Gσx=gxσx,gyσx=E01-μ0∑i=1n∂2Ni∂x2ui+μ0∂2Ni∂x∂yvi,00000000000∑i=1n∂2Ni∂y∂xui+μ0∂2Ni∂y2vi,where gx and gy are the components of the stress gradient, Ni is the shape function, n is the number of element nodes, and (ui,vi) is the nodal displacement vector; E0=E, μ0=μ for the plane stress, E0=E/(1-μ2), μ0=μ/(1-μ) for the plane strain; E 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 linearC0 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(2)∂∂ξ∂∂η=∂x∂ξ∂y∂ξ∂x∂η∂y∂η∂∂x∂∂y,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(3)∂2∂ξ2∂2∂η2∂2∂ξ∂η=∂2x∂ξ2∂2y∂ξ2∂2x∂η2∂2y∂η2∂2x∂ξ∂η∂2y∂ξ∂η∂∂x∂∂y+∂x∂ξ2∂y∂ξ22∂x∂ξ∂y∂ξ∂x∂η2∂x∂η22∂x∂η∂y∂η∂x∂ξ∂x∂η∂y∂ξ∂y∂η∂x∂ξ∂y∂η+∂x∂η∂y∂ξ×∂2∂x2∂2∂y2∂2∂x∂y,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.


### 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 Figures1(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.Figure 1
Construction of 2D patches with quadrilateral and triangular elements.
(a)
Q4
(b)
Q8
(c)
T3
(d)
T6With the SCP recovery, the recovered stress gradientgx,y∗ can be expressed as(4)gx,y∗=Pa,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(5)P=1xyx2xyy2,where the recovered stress gradient has the same order of displacement filed [29], and then a=a1a2⋯a6T.The functional integral of the difference between the stress gradient calculated from the SCP recoverygx,y∗ and the stress gradient obtained by the FEM g~x,y can be expressed as(6)Fgx,y∗,g~x,y=∑e=1m∫Segx,y∗-g~x,y2dS=∑e=1m∫SePa-g~x,y2dS,where m is the number of patch elements and Se is the element area.The assumed coefficient vectora can be obtained by minimizing the functional integral in (6), which yields(7)∑e=1m∫SePTPdSa=∑e=1m∫SePTg~x,ydS,where superscript T denotes the transpose, and both sides of (7) can be calculated by the Gauss integral:(8)∑i=1nPTxi,yiPxi,yiwia=∑i=1nPTxi,yig~x,yxi,yiwi,where n=m·ng, ng is the number of Gauss points per element, wi is the weight coefficient of each integration point, and (xi,yi) is the global coordinates of each Gauss point. Finally, there is(9)a=C-1D,in which(10)C=∑i=1nPTxi,yiPxi,yiwi,D=∑e=1nPTxi,yig~x,yxi,yiwi.In [25, 26], utilizing the quadratic completeness polynomial, the SCP recovery for stress gradient has the order of accuracy O(h3) (h is the element characteristic length). And then, it is easy to verify that the error of the stress gradient recovery is of order O(h2) by the same procedure. It has been found that the recovered variable field is more accurate at node points.


## 2.1. Second Derivatives in FEM

Taking the isotropic normal stressσx, for example, the gradient vector G(σx) can be written as(1)Gσx=gxσx,gyσx=E01-μ0∑i=1n∂2Ni∂x2ui+μ0∂2Ni∂x∂yvi,00000000000∑i=1n∂2Ni∂y∂xui+μ0∂2Ni∂y2vi,where gx and gy are the components of the stress gradient, Ni is the shape function, n is the number of element nodes, and (ui,vi) is the nodal displacement vector; E0=E, μ0=μ for the plane stress, E0=E/(1-μ2), μ0=μ/(1-μ) for the plane strain; E 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 linearC0 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(2)∂∂ξ∂∂η=∂x∂ξ∂y∂ξ∂x∂η∂y∂η∂∂x∂∂y,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(3)∂2∂ξ2∂2∂η2∂2∂ξ∂η=∂2x∂ξ2∂2y∂ξ2∂2x∂η2∂2y∂η2∂2x∂ξ∂η∂2y∂ξ∂η∂∂x∂∂y+∂x∂ξ2∂y∂ξ22∂x∂ξ∂y∂ξ∂x∂η2∂x∂η22∂x∂η∂y∂η∂x∂ξ∂x∂η∂y∂ξ∂y∂η∂x∂ξ∂y∂η+∂x∂η∂y∂ξ×∂2∂x2∂2∂y2∂2∂x∂y,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.


## 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 Figures1(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.Figure 1
Construction of 2D patches with quadrilateral and triangular elements.
(a)
Q4
(b)
Q8
(c)
T3
(d)
T6With the SCP recovery, the recovered stress gradientgx,y∗ can be expressed as(4)gx,y∗=Pa,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(5)P=1xyx2xyy2,where the recovered stress gradient has the same order of displacement filed [29], and then a=a1a2⋯a6T.The functional integral of the difference between the stress gradient calculated from the SCP recoverygx,y∗ and the stress gradient obtained by the FEM g~x,y can be expressed as(6)Fgx,y∗,g~x,y=∑e=1m∫Segx,y∗-g~x,y2dS=∑e=1m∫SePa-g~x,y2dS,where m is the number of patch elements and Se is the element area.The assumed coefficient vectora can be obtained by minimizing the functional integral in (6), which yields(7)∑e=1m∫SePTPdSa=∑e=1m∫SePTg~x,ydS,where superscript T denotes the transpose, and both sides of (7) can be calculated by the Gauss integral:(8)∑i=1nPTxi,yiPxi,yiwia=∑i=1nPTxi,yig~x,yxi,yiwi,where n=m·ng, ng is the number of Gauss points per element, wi is the weight coefficient of each integration point, and (xi,yi) is the global coordinates of each Gauss point. Finally, there is(9)a=C-1D,in which(10)C=∑i=1nPTxi,yiPxi,yiwi,D=∑e=1nPTxi,yig~x,yxi,yiwi.In [25, 26], utilizing the quadratic completeness polynomial, the SCP recovery for stress gradient has the order of accuracy O(h3) (h is the element characteristic length). And then, it is easy to verify that the error of the stress gradient recovery is of order O(h2) 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 Figure1 for 2D problem, the central point Pc(x0,y0) is surrounded by sample points Pjxj,yj, j=1,2,…,ns (ns is number of sample points). The vector from Pc to Pj is ncj=(xj-x0,yj-y0); thus, the unit vector can be written as n-cj=(1/PcPj)(xj-x0,yj-y0), where PcPj=(xj-x0)2+(yj-y0)22 is the vector length.The directional derivative alongncj at the central point Pc can be expressed as(11)∂σ~x∂n-cjPc=limP→PcP∈ncjσ~xP-σ~xPcPcP=1PcPjxPj-xPc∂σ~x∂xPc000000000+yPj-yPc∂σ~x∂yPc,where σ~x is the stress result of commercial finite element software. The limit can be approximately written as (12)limP→PcP∈ncjσ~xP-σ~xPcPcP≈σ~xPj-σ~xPcPcPj.Substituting (12) into (11) yields(13)ncj·GPcσ~xT=σ~xPj-σ~xPc.For all sample points, the stress gradient equation can be written as(14)A·GPcσ~xT=b,where(15)A=nc1nc2⋮ncns,b=σ~xP1-σ~xPcσ~xP2-σ~xPc⋮σ~xPns-σ~xPc.In general,ns>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 Ans×2. The pseudoinverse M2×ns can be calculated by the singular value decomposition of Ans×2; that is,(16)A=UΣVT,where U∈Rns×ns and V∈R2×2 are orthogonal matrices, Σ∈Rns×2 is a diagonal matrix with nonnegative real numbers, and the diagonal entries are known as the singular values. Then, (17)M=VΣ†UT,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 pointPc can be expressed by(18)GPcσ~xT=M·b.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 ncj=(xj-x0,yj-y0,zj-z0) and GPc(σ~x)=(∂σ~x/∂x,∂σ~x/∂y,∂σ~x/∂z), respectively. As an example, the 3D element patch is shown in Figure 2 with 8-node linear hexahedral elements.Figure 2
Construction of 3D patch with hexahedral element.


### 3.2. Error Estimation

The Moore-Penrose pseudoinverseMn0×ns (n0=2 or 3 for 2D or 3D case) of Ans×n0∈Rns×n0 is unique and satisfies all of the following four criteria:(19)AMA=A,MAM=M,AMT=AM,MAT=MA,which yields(20)AMAT=ATAMT=ATAM=AT,MAMT=MTMAT=MTMA=MT.For all g∈Rn0×1, there is(21)Ag-b=Ag-b+AMb-AMb=AMb-b+Aw=AM-Ib+Aw,where w=g-Mb∈Rn0×1. With the first equation of (20), AT(AM-I)=O and thus there is a relationship. Consider(22)AM-Ib,Aw=0,in which 〈,〉 denotes the inner product. With this relationship, we can have(23)AM-Ib+Aw=AM-Ib+Aw≥AM-Ib,or(24)AMb-b≤Ag-b,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(25)g^=M^b+I-M^Az,∀z∈Rn0×1,M^∈Rn0×nsin which AM^A=A and (AM^)T=AM^.Similarly, it can be verified that ifM^AM^=M^,(M^A)T=M^A, then(26)M^b≤M^b+I-M^Az.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 O(h) in (12) and the other one is the stress field σ~x. 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: (27)σ~x=w1σ-x1+w2σ-x2+⋯+wmσ-xmw1+w2+⋯+wm,in which m is the number of surrounding elements for one node, wi(i=1,2,…,m) is the weight (element area or volume), andσ-xi(i=1,2,…,m) is the stress calculated by the node displacements of each element. The stress field σ~x evaluated by (27) has the order of accuracy O(h) at least. Therefore, the order of the minimum norm least squares solution of (14) is of O(h) 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.


## 3.1. General Formulation

As shown in Figure1 for 2D problem, the central point Pc(x0,y0) is surrounded by sample points Pjxj,yj, j=1,2,…,ns (ns is number of sample points). The vector from Pc to Pj is ncj=(xj-x0,yj-y0); thus, the unit vector can be written as n-cj=(1/PcPj)(xj-x0,yj-y0), where PcPj=(xj-x0)2+(yj-y0)22 is the vector length.The directional derivative alongncj at the central point Pc can be expressed as(11)∂σ~x∂n-cjPc=limP→PcP∈ncjσ~xP-σ~xPcPcP=1PcPjxPj-xPc∂σ~x∂xPc000000000+yPj-yPc∂σ~x∂yPc,where σ~x is the stress result of commercial finite element software. The limit can be approximately written as (12)limP→PcP∈ncjσ~xP-σ~xPcPcP≈σ~xPj-σ~xPcPcPj.Substituting (12) into (11) yields(13)ncj·GPcσ~xT=σ~xPj-σ~xPc.For all sample points, the stress gradient equation can be written as(14)A·GPcσ~xT=b,where(15)A=nc1nc2⋮ncns,b=σ~xP1-σ~xPcσ~xP2-σ~xPc⋮σ~xPns-σ~xPc.In general,ns>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 Ans×2. The pseudoinverse M2×ns can be calculated by the singular value decomposition of Ans×2; that is,(16)A=UΣVT,where U∈Rns×ns and V∈R2×2 are orthogonal matrices, Σ∈Rns×2 is a diagonal matrix with nonnegative real numbers, and the diagonal entries are known as the singular values. Then, (17)M=VΣ†UT,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 pointPc can be expressed by(18)GPcσ~xT=M·b.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 ncj=(xj-x0,yj-y0,zj-z0) and GPc(σ~x)=(∂σ~x/∂x,∂σ~x/∂y,∂σ~x/∂z), respectively. As an example, the 3D element patch is shown in Figure 2 with 8-node linear hexahedral elements.Figure 2
Construction of 3D patch with hexahedral element.


## 3.2. Error Estimation

The Moore-Penrose pseudoinverseMn0×ns (n0=2 or 3 for 2D or 3D case) of Ans×n0∈Rns×n0 is unique and satisfies all of the following four criteria:(19)AMA=A,MAM=M,AMT=AM,MAT=MA,which yields(20)AMAT=ATAMT=ATAM=AT,MAMT=MTMAT=MTMA=MT.For all g∈Rn0×1, there is(21)Ag-b=Ag-b+AMb-AMb=AMb-b+Aw=AM-Ib+Aw,where w=g-Mb∈Rn0×1. With the first equation of (20), AT(AM-I)=O and thus there is a relationship. Consider(22)AM-Ib,Aw=0,in which 〈,〉 denotes the inner product. With this relationship, we can have(23)AM-Ib+Aw=AM-Ib+Aw≥AM-Ib,or(24)AMb-b≤Ag-b,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(25)g^=M^b+I-M^Az,∀z∈Rn0×1,M^∈Rn0×nsin which AM^A=A and (AM^)T=AM^.Similarly, it can be verified that ifM^AM^=M^,(M^A)T=M^A, then(26)M^b≤M^b+I-M^Az.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 O(h) in (12) and the other one is the stress field σ~x. 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: (27)σ~x=w1σ-x1+w2σ-x2+⋯+wmσ-xmw1+w2+⋯+wm,in which m is the number of surrounding elements for one node, wi(i=1,2,…,m) is the weight (element area or volume), andσ-xi(i=1,2,…,m) is the stress calculated by the node displacements of each element. The stress field σ~x evaluated by (27) has the order of accuracy O(h) at least. Therefore, the order of the minimum norm least squares solution of (14) is of O(h) 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 Figure3 becomes a classic and popular two-dimensional example used in previous papers [13–15, 24, 29].Figure 3
An infinite plate with a central hole subjected to a remotely uniform tension: (a) theory model, (b) finite element model.
(a)
(b)In this case, the normal stressσx is dominant and the exact solution is given as(28)σx=P01-r02r232cos⁡2θ+cos⁡4θ+32r04r4cos⁡4θ,where r 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 Figure4, including the numbers of nodes and elements.Figure 4
Meshes of quadrilateral and triangular 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 pointP around (2.2r0,2.2r0) is used to compare the results of different methods as shown in Figure 4. The exact locations are (2.212r0,2.212r0), (2.225r0,2.225r0), and (2.232r0,2.232r0) 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 x-direction gradient (Sxx) error and y-direction gradient (Sxy) 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.Table 1
Stress gradients at pointP with different meshes and methods.
Mesh-Method-Element
Numerical solution
Exact solution
Error
S
x
x/108
S
x
y/108
S
x
x/108
S
x
y/108
S
x
x
S
x
y
Average
C-SFG-Q4
−10.3667
4.0309
−10.1284
3.7238
2.3524%
8.2482%
5.3003%
C-SCP-Q8
−8.7140
3.3071
−13.9646%
−11.1893%
12.5769%
C-SFG-Q8
−1.0908
4.1159
7.6936%
10.5311%
9.1123%
C-SFG-T3
−9.2718
4.8073
−8.4575%
29.0992%
18.7783%
C-SCP-T6
−10.1010
3.9277
−0.2813%
5.4765%
2.8789%
C-SFG-T6
−10.2858
3.6292
1.5533%
−2.5402%
2.0467%
F-SFG-Q4
−10.1920
3.6648
−9.9690
3.6431
2.2362%
0.5960%
1.4161%
F-SCP-Q8
−9.54347
3.2959
−4.2689%
−9.5305%
6.8997%
F-SFG-Q8
−10.4925
3.8314
5.2506%
5.1685%
5.2096%
F-SFG-T3
−10.1635
3.7822
1.9511%
3.8174%
2.8842%
F-SCP-T6
−10.2152
3.6521
2.4697%
0.2479%
1.3588%
F-SFG-T6
−10.1343
3.5049
1.6575%
−3.7928%
2.7251%
DF-SFG-Q4
−10.1587
3.6648
−9.8932
3.6050
2.6837%
1.6590%
2.1714%
DF-SCP-Q8
−9.8485
3.3578
−0.4519%
−6.8582%
3.6551%
DF-SFG-Q8
−10.2297
3.5939
3.4010%
−0.3094%
1.8552%
DF-SFG-T3
−9.9612
3.7442
0.6875%
3.8603%
2.2739%
DF-SCP-T6
−10.1757
3.5357
2.8555%
−1.9237%
2.3896%
DF-SFG-T6
−10.1554
3.5479
2.6497%
−1.5834%
2.1166%
Meshes: C: coarse; F: fine; DF: double-fine.
Methods: SFG: stress field gradient; SCP: superconvergent patch.Furthermore, results from the two methods along two checking lines,θ=45° and θ=90°, 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 C0 elements provides the worst stress gradient for line θ=45° (Line-45°) as well as the SCP method using quadratic Q8 element for line θ=90° (Line-Y). 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.Figure 5
Stress gradients of sample lines with quadrilateral elements.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Stress gradients of sample lines with triangular elements.
(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 Figure7. All DOFs on the x=0 surface are fixed.Figure 7
Geometry, loading condition, and mesh of 3D problem.
(a)
(b)Based on the software platform ANSYS, the flowchart of postprocessing with the proposed method is shown in Figure8. 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.Figure 8
The flowchart of postprocessing using the presented method.Details of the stress gradients are shown in Figure9 for h=4mm and d=5mm, in which the total gradient G(σ~)Total=(∂σ~/∂x)2+(∂σ~/∂y)2+(∂σ~/∂z)22. 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.Figure 9
Von-Mises stress and its gradients of 3D problem: (a) Von-Mises stress, (b) Von-Mises stress gradient-total, (c) stress gradient slices-total, (d) stress gradient-X, (e) stress gradient-Y, and (f) stress gradient-Z.
(a)
(b)
(c)
(d)
(e)
(f)Results of different groove widths and depths obtained by the SFG technique are compared in Figure10 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.15E+10 to 2.52E+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.0531E+10 to 2.3719E+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 x-direction.Table 2
Maximal Von-Mises stresses and gradients of different groove widths and depths.
h/mm
d/mm
Max-Von-Mises/Pa
Max_Total gradient
5
5
1.4339E+ 08
4.1542E+ 10
10
5
1.2381E+ 08
2.4833E+ 10
20
5
1.1546E+ 08
2.5172E+ 10
4
4
1.3163E+ 08
4.0531E+ 10
4
8
2.3143E+ 08
9.2233E+ 10
4
12
4.6672E+ 08
2.3719E+ 11Figure 10
Von-Mises total stress gradients for different groove widths and depths: (a)h=5 mm, d=5 mm; (b) h=10 mm, d=5 mm; (c) h=20 mm, d=5 mm; (d) h=4 mm, d=4 mm; (e) h=4 mm, d=8 mm; (f) h=4 mm, d=12 mm.
(a)
(b)
(c)
(d)
(e)
(f)Figure 11
Total gradients of middle lines for different groove widths and depths.
(a)
(b)


## 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 Figure3 becomes a classic and popular two-dimensional example used in previous papers [13–15, 24, 29].Figure 3
An infinite plate with a central hole subjected to a remotely uniform tension: (a) theory model, (b) finite element model.
(a)
(b)In this case, the normal stressσx is dominant and the exact solution is given as(28)σx=P01-r02r232cos⁡2θ+cos⁡4θ+32r04r4cos⁡4θ,where r 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 Figure4, including the numbers of nodes and elements.Figure 4
Meshes of quadrilateral and triangular 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 pointP around (2.2r0,2.2r0) is used to compare the results of different methods as shown in Figure 4. The exact locations are (2.212r0,2.212r0), (2.225r0,2.225r0), and (2.232r0,2.232r0) 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 x-direction gradient (Sxx) error and y-direction gradient (Sxy) 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.Table 1
Stress gradients at pointP with different meshes and methods.
Mesh-Method-Element
Numerical solution
Exact solution
Error
S
x
x/108
S
x
y/108
S
x
x/108
S
x
y/108
S
x
x
S
x
y
Average
C-SFG-Q4
−10.3667
4.0309
−10.1284
3.7238
2.3524%
8.2482%
5.3003%
C-SCP-Q8
−8.7140
3.3071
−13.9646%
−11.1893%
12.5769%
C-SFG-Q8
−1.0908
4.1159
7.6936%
10.5311%
9.1123%
C-SFG-T3
−9.2718
4.8073
−8.4575%
29.0992%
18.7783%
C-SCP-T6
−10.1010
3.9277
−0.2813%
5.4765%
2.8789%
C-SFG-T6
−10.2858
3.6292
1.5533%
−2.5402%
2.0467%
F-SFG-Q4
−10.1920
3.6648
−9.9690
3.6431
2.2362%
0.5960%
1.4161%
F-SCP-Q8
−9.54347
3.2959
−4.2689%
−9.5305%
6.8997%
F-SFG-Q8
−10.4925
3.8314
5.2506%
5.1685%
5.2096%
F-SFG-T3
−10.1635
3.7822
1.9511%
3.8174%
2.8842%
F-SCP-T6
−10.2152
3.6521
2.4697%
0.2479%
1.3588%
F-SFG-T6
−10.1343
3.5049
1.6575%
−3.7928%
2.7251%
DF-SFG-Q4
−10.1587
3.6648
−9.8932
3.6050
2.6837%
1.6590%
2.1714%
DF-SCP-Q8
−9.8485
3.3578
−0.4519%
−6.8582%
3.6551%
DF-SFG-Q8
−10.2297
3.5939
3.4010%
−0.3094%
1.8552%
DF-SFG-T3
−9.9612
3.7442
0.6875%
3.8603%
2.2739%
DF-SCP-T6
−10.1757
3.5357
2.8555%
−1.9237%
2.3896%
DF-SFG-T6
−10.1554
3.5479
2.6497%
−1.5834%
2.1166%
Meshes: C: coarse; F: fine; DF: double-fine.
Methods: SFG: stress field gradient; SCP: superconvergent patch.Furthermore, results from the two methods along two checking lines,θ=45° and θ=90°, 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 C0 elements provides the worst stress gradient for line θ=45° (Line-45°) as well as the SCP method using quadratic Q8 element for line θ=90° (Line-Y). 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.Figure 5
Stress gradients of sample lines with quadrilateral elements.
(a)
(b)
(c)
(d)
(e)
(f)Figure 6
Stress gradients of sample lines with triangular elements.
(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 Figure7. All DOFs on the x=0 surface are fixed.Figure 7
Geometry, loading condition, and mesh of 3D problem.
(a)
(b)Based on the software platform ANSYS, the flowchart of postprocessing with the proposed method is shown in Figure8. 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.Figure 8
The flowchart of postprocessing using the presented method.Details of the stress gradients are shown in Figure9 for h=4mm and d=5mm, in which the total gradient G(σ~)Total=(∂σ~/∂x)2+(∂σ~/∂y)2+(∂σ~/∂z)22. 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.Figure 9
Von-Mises stress and its gradients of 3D problem: (a) Von-Mises stress, (b) Von-Mises stress gradient-total, (c) stress gradient slices-total, (d) stress gradient-X, (e) stress gradient-Y, and (f) stress gradient-Z.
(a)
(b)
(c)
(d)
(e)
(f)Results of different groove widths and depths obtained by the SFG technique are compared in Figure10 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.15E+10 to 2.52E+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.0531E+10 to 2.3719E+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 x-direction.Table 2
Maximal Von-Mises stresses and gradients of different groove widths and depths.
h/mm
d/mm
Max-Von-Mises/Pa
Max_Total gradient
5
5
1.4339E+ 08
4.1542E+ 10
10
5
1.2381E+ 08
2.4833E+ 10
20
5
1.1546E+ 08
2.5172E+ 10
4
4
1.3163E+ 08
4.0531E+ 10
4
8
2.3143E+ 08
9.2233E+ 10
4
12
4.6672E+ 08
2.3719E+ 11Figure 10
Von-Mises total stress gradients for different groove widths and depths: (a)h=5 mm, d=5 mm; (b) h=10 mm, d=5 mm; (c) h=20 mm, d=5 mm; (d) h=4 mm, d=4 mm; (e) h=4 mm, d=8 mm; (f) h=4 mm, d=12 mm.
(a)
(b)
(c)
(d)
(e)
(f)Figure 11
Total gradients of middle lines for different groove widths and depths.
(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 linearC0 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.



---

*Source: 457046-2015-05-20.xml*