The Incremental Hybrid Natural Element Method for Elastoplasticity Problems

An incremental hybrid natural element method (HNEM) is proposed to solve the two-dimensional elasto-plastic problems in the paper. The corresponding formulae of this method are obtained by consolidating the hybrid stress element and the incremental Hellinger-Reissner variational principle into the NEM. Using this method, the stress and displacement variables at each node can be directly obtained after the stress and displacement interpolation functions are properly constructed. The numerical examples are given to show the advantages of the proposed algorithm of the HNEM, and the solutions for the elasto-plastic problems are better than those of the NEM. In addition, the performance of the proposed algorithm is better than the recover stress method using moving least square interpolation.


Introduction
The natural element method (NEM) [1,2] is an algorithm for building the approximate functions or trial functions over the whole region by natural neighbor nodes interpolation method. The NEM interpolation is constructed on the basis of the underlying Voronoi tessellation. The Delaunay triangles which are the dual of the Voronoi diagram are used in the numerical computation of the NEM interpolation. Unlike the finite element method where angle restrictions are imposed on the triangles for the convergence of the method, there are no such constraints on the shape, size, and angles of the triangles in NEM [3]. The NEM is widely applied to numerically solve partial differential equations (PDEs) since it has advantages over the conventional FEM such as the convenience in preprocess, the simple shape function, and it satisfies the Delta interpolation on the boundary [4][5][6].
The elastoplasticity problem is an important nonlinear problem in solid mechanics. Based on the NEM, some research work to solve the elastoplastic problems has been implemented over the last twenty years. Toi and Kang [7][8][9] deduced NEM for elastoplasticity problems by combining natural neighbor nodes interpolation method and incremental shape Galerkin method. Ding et al. [10] formatted elastoplasticity natural element method and calculated actual engineering problems. Hehua et al. [11] utilized the meshless local natural neighbor interpolation method (MLNNIM) to solve problems of two-dimensional elastoplastic large deformation. In recent years, with the development of the numerical methods, other meshless methods [12][13][14] have been used to solve the elastoplasticity problems. On the basis of reproducing kernel particle method (RKPM), using complex variable theory, the complex variable reproducing kernel particle method (CVRKPM) for elastoplasticity problems is discussed by Chen and Cheng [12], and based on the complex variable moving least-squares (CVMLS) approximation and element-free Galerkin (EFG) method, the complex variable element-free Galerkin (CVEFG) method for twodimensional elastoplasticity problems is presented by Peng et al. [13]. The formulation and numerical implementation of the improved complex variable element-free Galerkin (ICVEFG) method is presented for two-dimensional large deformation problems of elastoplasticity in total Lagrangian description by Li et al. [14].  Due to the 0 continuity of traditional NEM shape function on nodes, we cannot obtain smooth and continuous differential functions. Thus we cannot directly solve the stresses on nodes for solid structure analysis. Iteration is necessary for elastoplasticity mechanical problems; the precision of last step stresses influences the evaluation of next iterative step; the accumulated errors may yield an incorrect solution. To improve the accuracy of stresses, we can adopt methods such as intensive nodes, smoothing shape functions [15], and recovering stresses using high order meshless interpolation [16]. However, these methods in line space will cost more calculation time and negatively influence the performance.
By introducing the hybrid stress element into the natural element method and combining the incremental Hellinger-Reissner variational principle [17], an algorithm suitable for solving the elastoplastic problems, the incremental hybrid natural element method (HNEM), is proposed in the paper. The effectiveness and the advantages of the HNEM are verified through the presented numerical examples.

Natural Neighbor Nodes Interpolation Method
For the region Ω of any problem, it can be represented by discrete points as shown in Figure 1. According to Delaunay empty circle principle, the structure can be gridded automatically using triangulation, and then the Voronoi structure of the region Ω can be built based on the information of triangle grids. The Voronoi structure for a discrete point x can be expressed as: where (x, x ) is distance between any point x and the discrete point x . Suppose x as any point within the region; the natural neighbor nodes around x are those found according Delaunay empty circle principle. The natural neighbor nodes of x are nodes 1, 2, 3, 6, and 7 (see Figure 2).
During the process of problem-solving, the natural neighbor nodes of x must be determined at first. Then we can build a local interpolation format using natural neighbor nodes interpolation as follows: where (x) is the interpolation physical quantity for node x, is the serial number of natural neighbor nodes around x, is the total number of natural neighbor nodes around x, is the physical quantity value of node , and Φ (x) is the base function of interpolation of node x . The base function can be chosen between Sibson or non-Sibsonian base function of interpolation.
The interpolating Sibson base function can be written as where (x) is the area of second order natural neighbor nodes around x within first order Voronoi element of node x and (x) is the total area of second order natural neighbor nodes. Engineering   3 For mechanics problem we can use the following to be displacement interpolation function:

Mathematical Problems in
and then solve the problem using discrete control equation using Galerkin approximation. Here Φ (x) can be Sibson interpolation or non-Sibsonian interpolation.
Because the shape functions of natural neighbor nodes are ∞ continuous on region except nodes and 0 continuous on nodes, we can get high precise displacement using displacement interpolation only. But we cannot obtain the stresses on nodes and high precise stresses. Iteration is necessary for elastoplasticity mechanical problems; the precision of last step stresses influences the evaluation of next iterative step; the accumulated errors may yield an incorrect solution. So the incremental hybrid natural element method is an excellent candidate for elastoplasticity mechanic problems when using displacement interpolation and stress interpolation, respectively, and combining incremental multivariable Hellinger-Reissner variational principle in this paper.

Basic Equations of Elastoplasticity Plane Problem.
For two-dimensional problem, the displacement, strain, and stress fields of a body are defined to be u 0 , 0 , and 0 , respectively, after random load history. When the load increases persistently, the increments of body force, surface force, and boundary displacement areḃ , t, andu , respectively, then the increments of displacement, strain, and stress areu ,, and, where the symbol "⋅" means incremental quantity. Thus the basic incremental equations of elastoplasticity mechanic can be expressed as follows.
The equilibrium equation can be written as where L is the matrix of differential operator for 2D problems: wherėis the stress increment of point x;ḃ is the body force increment of point x.
The strain-displacement relationship can be written aṡ wherėandu are strain and displacement increment of random point x, respectively. The stress-strain relationship can be expressed aṡ = D.
(9) During elasticity stage, Within plastic stage where is the equivalent stress, where is the deviatoric stress, and is the hardening parameter [18]. The boundary conditions can be written aṡ 4 Mathematical Problems in Engineering whereu is the prescribed displacement increment at an arbitrary point x on the essential boundary Γ ,ṫ is the prescribed surface force increment at an arbitrary point x on the natural boundary Γ , Γ is the boundary of the domain Ω, and Γ = Γ ∪ Γ , Γ ∩ Γ = 0.

Incremental Hybrid Natural Element
Method. Let us suppose that there exist nodes in domain Ω, as shown in Figure 1; the displacement increment and stress increment of an arbitrary node can be expressed aṡ where is the total number of natural neighbor nodes around x and Φ and Φ are interpolating functions of displacement incrementu and stress increment, which can be found from (3). Equations of (16) can be written in matrix as follows: where Applying the incremental Hellinger-Reissner variation principle given by [19] to elastoplastic analysis, we have Among above () is incremental surplus energy density function, for elastoplastic problem where S is the flexibility matrix, S = D −1 , and D is shown in (10) and (11); (⋅) 0 give the values of every variation in the initial state;Ḟ is the body force column vector; andṪ is the surface force column vector.
Suppose the displacement boundary condition (14) is satisfied, and substitute (17) and (18) into multivariation Hellinger-Reissner variation principle; we can obtain where where R is the unbalance correction term [20,21]; it is related to load history.

Mathematical Problems in
According to stationary value condition of multivariation increment variation principle, We havė= Substituting (21) into (18) yields where Based on the stationary condition of variation, The equations for solving generalized displacement can be found: For the elastoplasticity problem with small deformation, the relationship of stress and strain in the plastic field is nonlinear. The Newton-Raphson iterative solver is used to solve nonlinear system. We convert the nonlinear system into a series of linear problems using the method of gradually increasing loading. After the structure enters the yield state, the load is added with load increment method. The structure is added with the load increment ΔF; the residual forces Δ is defined: For the solution of (30) and for each load increment consider the situation existing for the th iteration. The applied loads for the th iteration are the residual forces Δ −1 calculated at the end of the ( − 1)th iteration according to (30). These applied loads give rise to displacement incrementsq . The corresponding increments of the straiṅ are calculated. The incremental stress change assuming elastic behavior is computed. The computed stresses are then brought down to the yield surface and are used to calculate the equivalent nodal forces. These nodal forces can be compared with the externally applied loads to form the residual forces for the next iteration. The system of residual forces is brought sufficiently close to zero through the iterative process, before moving to the next load increment.  Figure 4 gives the results solved with natural element method and hybrid natural element method under different nodes layout. From the comparison we can see that, under the same nodes layout, hybrid natural element method has higher calculation precision. Figure 5 shows the relationship between the deflection of middle point at the right end of the beam and the load. If the load is small, the deformation of the beam is linear elasticity; when the load is bigger than the yield load, the material begins to yield and enters into the plastic state, and the relation between deformation and load is nonlinear. Tables 1 and 2 give the normal stress 11 and the shear stress 12 at some Gaussian quadrature points solved with hybrid natural element method and natural element method. From the comparison with the analytical stresses we can see that, at the same Gaussian quadrature points, hybrid natural element method has higher calculation precision to demonstrate the effectiveness and the advantages of HNEM. Figure 6 is a cylinder with uniform pressure on inner surface, the inner radius of the cylinder is = 1 m, the outer radius is = 5 m, and the uniform pressure is = 1000 N/m. The cylinder can be calculated by plane strain, and the material can choose bilinear strengthening model and von-Mises criterion. The material parameters are Young's modulus, = 10 5 Pa, Poisson's ratio, ] = 0.25, and the yield stress, = 150 Pa. The linear hardening elastoplastic model is adopted with = 0.2 .

Cylinder with Uniform Inner Pressure.
We can only consider one quarter of the cylinder because of the symmetry, as shown in Figure 7. We can arrange 99 nodes totally as shown in Figure 8. Figure 9 gives numerical solution of axial displacement of partial nodes on 1 -axis, and the solution is compared with that solved by ABAQUS. It can be seen that the calculation results of hybrid natural element   method fit well with that of ABAQUS. From the results we can see that, comparing with the results by recovering stress method using moving least square interpolation, the method of this paper can cut down the time cost in calculation and improve computational efficiency on condition of similar calculation precision. Figure 10 shows relationship between radial displacement and the load of node (5.0, 0.0) on the circle. We can see that the material turns into elastoplasticity from elasticity with the load increasing; when the load is bigger than the yield load, the material begins to yield and enters into the plastic state, and the displacement appears nonlinear variation with the increasing of the load.

Plate with a Central Hole under Axial Uniform Tension.
A plate with a central hole is fixed at one end and subjected to axial uniform tension at the other end, as shown in Figure 11.   The length of the plate is = 10 m, and the width is ℎ = 4 m. The radius of the hole is 1 meter, and the tensile load is = 1000 N/m. The plate can be considered as plane stress state, and the material can use bilinear strengthening model and von-Mises criterion. The parameters are Young's modulus, = 10 5 Pa, Poisson's ratio, ] = 0.25, and the yield stress, = 250 Pa. The linear hardening elastoplastic model is adopted with = 0.2 .
The nodes layout of the plate is given in Figure 12; there are 200 nodes in total. The nodes near the hole are intensive so as to improve the calculation precision. Figures 13 and  14 are comparison chart of horizontal displacement and longitudinal displacement on line 2 = 2 solved by FEM and the method of this paper, respectively. We can find that the results calculated using hybrid natural element method are close to those using ABAQUS. Figure 15 shows the curve of   horizontal displacement of node (5, 0) with the increasing of the load. It can be seen that when the load is bigger than the elastic limit of the material, the material begins to yield and enters into the plastic state, and show nonlinear behavior.

Conclusions
Based on the natural element method (NEM), by introducing the hybrid stress element into the NEM and combining the incremental Hellinger-Reissner variational principle, the incremental hybrid natural element method (HNEM) has been proposed for solving two-dimensional elastoplastic problems in the paper.  solutions for the elastoplastic problems can be obtained by the HNEM than those of the NEM. In addition, the efficiency of the proposed algorithm has been improved compared with the recover stress method using moving least square interpolation.
Using the HNEM in this paper, the stress values and displacement values of nodes can be easily got at the same time by structuring suitable node stress function and displacement interpolation function. Then not only the advantages of the convenience in pre-process, the simple shape function, and it satisfies the Delta interpolation on the boundary of the NEM are reserved by the HNEM, but also the problem of not directly solve the stresses on nodes is resolved and the accuracy of stresses is improved.