#### Abstract

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.

#### 1. 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–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–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–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 two-dimensional 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 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.

#### 2. 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 can be expressed as: where is distance between any point and the discrete point .

Suppose as any point within the region; the natural neighbor nodes around are those found according Delaunay empty circle principle. The natural neighbor nodes of are nodes 1, 2, 3, 6, and 7 (see Figure 2).

During the process of problem-solving, the natural neighbor nodes of must be determined at first. Then we can build a local interpolation format using natural neighbor nodes interpolation as follows: where is the interpolation physical quantity for node , is the serial number of natural neighbor nodes around , is the total number of natural neighbor nodes around , is the physical quantity value of node , and is the base function of interpolation of node . 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 is the area of second order natural neighbor nodes around within first order Voronoi element of node and is the total area of second order natural neighbor nodes.

For mechanics problem we can use the following to be displacement interpolation function: and then solve the problem using discrete control equation using Galerkin approximation. Here can be Sibson interpolation or non-Sibsonian interpolation.

Because the shape functions of natural neighbor nodes are continuous on region except nodes and 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.

#### 3. The Incremental Hybrid Natural Element Method for Elastoplasticity Mechanic

##### 3.1. Basic Equations of Elastoplasticity Plane Problem

For two-dimensional problem, the displacement, strain, and stress fields of a body are defined to be , , and , respectively, after random load history. When the load increases persistently, the increments of body force, surface force, and boundary displacement are , , and , respectively, then the increments of displacement, strain, and stress are , , 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 is the matrix of differential operator for 2D problems: where is the stress increment of point ; is the body force increment of point .

The strain-displacement relationship can be written as
where and are strain and displacement increment of random point **,** respectively.

The stress-strain relationship can be expressed as

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 as where is the prescribed displacement increment at an arbitrary point on the essential boundary , is the prescribed surface force increment at an arbitrary point on the natural boundary , is the boundary of the domain , and , .

##### 3.2. 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 as where is the total number of natural neighbor nodes around and and are interpolating functions of displacement increment 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 is the flexibility matrix, , and is shown in (10) and (11); 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 is the unbalance correction term [20, 21]; it is related to load history.

According to stationary value condition of multivariation increment variation principle,

We have

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 ; 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 calculated at the end of the ()th iteration according to (30). These applied loads give rise to displacement increments . The corresponding increments of the strain 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.

#### 4. Numerical Examples

##### 4.1. Cantilever Beam with Concentrated Load at the End

Figure 3 shows a cantilever beam with concentrated load, N, at the end; the length of the beam is m, and the height is m; this is a plane stress problem. In the example we ignore the gravity. The material elastoplasticity model can choose bilinear model and Tresca yield condition. Parameters of the beam are Young’s modulus, Pa, Poisson’s ratio, , and the yield stress, Pa. The linear hardening elastoplastic model is adopted with .

**(a)**

**(b)**

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 and the shear stress 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.

##### 4.2. Cylinder with Uniform Inner Pressure

Figure 6 is a cylinder with uniform pressure on inner surface, the inner radius of the cylinder is m, the outer radius is m, and the uniform pressure is 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, Pa, Poisson’s ratio, , and the yield stress, Pa. The linear hardening elastoplastic model is adopted with .

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 -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.

##### 4.3. 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 m, and the width is m. The radius of the hole is 1 meter, and the tensile load is 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, Pa, Poisson’s ratio, , and the yield stress, Pa. The linear hardening elastoplastic model is adopted with .

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 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.

#### 5. 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. The corresponding formulae of this method are deduced and the relevant programs of this method are compiled. Three numerical examples of two-dimensional elastoplastic problems are given, showing that with the proposed algorithm of the HNEM, the more accurate 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.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The work was supported by the Natural Science Foundation of Shanghai China (no. 13ZR1415900) and Program of Shanghai Science and technology Commission (no. 11231202700).