# The Incremental Hybrid Natural Element Method for Elastoplasticity Problems

Authors: Yongqi Ma; Yi Dong; Yankai Zhou; Wei Feng
Journal: Mathematical Problems in Engineering
 (2014)

---

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


---

## Body

## 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 theC0 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 nt 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 xi can be expressed as:
(1)Ti={x∈R2∣d(x,xi)<d(x,xj),∀j≠i},
where d(x,xi) is distance between any point x and the discrete point xi.Figure 1
Discrete model of regionΩ and Voronoi cell about xi.Supposex 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).Figure 2
The 2nd order Voronoi cell aboutx.During the process of problem-solving, the natural neighbor nodes ofx must be determined at first. Then we can build a local interpolation format using natural neighbor nodes interpolation as follows:
(2)f(x)=∑i=1nΦi(x)fi,
where f(x) is the interpolation physical quantity for node x, i is the serial number of natural neighbor nodes around x, n is the total number of natural neighbor nodes around x, fi is the physical quantity value of node i, and Φi(x) is the base function of interpolation of node xi. The base function can be chosen between Sibson or non-Sibsonian base function of interpolation.The interpolating Sibson base function can be written as(3)Φi=Ai(x)A(x),(4)A(x)=∑j=1nAj(x),
where Ai(x) is the area of second order natural neighbor nodes around x within first order Voronoi element of node xi and A(x) is the total area of second order natural neighbor nodes.For mechanics problem we can use the following to be displacement interpolation function:(5)uh(x)=∑i=1nΦi(x)ui
and then solve the problem using discrete control equation using Galerkin approximation. Here Φi(x) can be Sibson interpolation or non-Sibsonian interpolation.Because the shape functions of natural neighbor nodes areC∞ continuous on region except nodes and C0 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 beu0, ε0, and σ0, respectively, after random load history. When the load increases persistently, the increments of body force, surface force, and boundary displacement are b˙, t-, and u-˙, respectively, then the increments of displacement, strain, and stress are u˙, ε˙, 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(6)LTσ˙+b˙=0,(inΩ),
where L is the matrix of differential operator for 2D problems:
(7)L=[∂∂x100∂∂x2∂∂x2∂∂x1],
where σ˙ is the stress increment of point x; b˙ is the body force increment of point x.The strain-displacement relationship can be written as(8)ε˙=Lu˙,(inΩ),
where ε˙ and u˙ are strain and displacement increment of random point x, respectively.The stress-strain relationship can be expressed as(9)σ˙=Dε˙.During elasticity stage,(10)D=[D]e.Within plastic stage(11)D=[D]ep=EQ[σ22′2+2Ps-σ11′σ22′+2νPs-σ11′+νσ22′1+νσ12-σ11′σ22′+2νPsσ11′+2Ps-σ22′+νσ11′1+νσ12-σ11′+νσ22′1+νσ12-σ22′+νσ11′1+νσ12Rs2(1+ν)+2H′9E(1-ν)σ-2],whereσ- is the equivalent stress,
(12)σ-=σ112+σ222-σ11σ22+3σ122,
where σ′ is the deviatoric stress,
(13)σ11′=σ11-σ11+σ123,σ22′=σ22-σ11+σ123,Ps=2H′9Eσ-2+σ1221+ν,Rs=σ11′2+2νσ11′σ22′+σ22′2,Q=Rs+2(1-ν2)Ps,
and H′ is the hardening parameter [18].The boundary conditions can be written as(14)u˙=u-˙,(on⁡Γu),(15)n·σ˙=t-˙,(on⁡Γt),
where u˙ is the prescribed displacement increment at an arbitrary point x on the essential boundary Γu, t-˙ is the prescribed surface force increment at an arbitrary point x on the natural boundary Γt, Γ is the boundary of the domain Ω, and Γ=Γu∪Γt, Γu∩Γt=∅.


### 3.2. Incremental Hybrid Natural Element Method

Let us suppose that there existnt nodes in domain Ω, as shown in Figure 1; the displacement increment and stress increment of an arbitrary node can be expressed as
(16)u˙h(x)=∑i=1nΦiu(x)u˙i,σ˙h(x)=∑i=1nΦiσ(x)σ˙i,
where n is the total number of natural neighbor nodes around x and Φiu and Φiσ are interpolating functions of displacement increment u˙ and stress increment σ˙, which can be found from (3).Equations of (16) can be written in matrix as follows:
(17)u˙h(x)=𝚽uq˙,(18)σ˙h(x)=Pβ˙,
where(19)𝚽u=[Φ10Φ20⋯Φn00Φ10Φ2⋯0Φn],q˙=[u˙1(x1)u˙2(x1)u˙1(x2)u˙2(x2)⋯u˙1(xn)u˙2(xn)],P=[Φ100Φ200⋯Φn000Φ100Φ20⋯0Φn000Φ100Φ2⋯00Φn],β˙=[σ˙11(x1)σ˙22(x1)σ˙12(x1)σ˙11(x2)σ˙22(x2)σ˙12(x2)⋯σ˙11(xn)σ˙22(xn)σ˙12(xn)].Applying the incremental Hellinger-Reissner variation principle given by [19] to elastoplastic analysis, we have
(20)πHR=∫Ω[(F-0T+F-T)B(σ˙)-(σ0T-σ˙)(Lu˙)+σ˙T(ε0-Lu˙)+(F-0T+F-T)u˙]dΩ+∫Γσ(T-0T+T-˙T)u˙dΓσ+∫ΓuT-˙T(u˙-u-˙+u0-u-0)+T-0Tu˙dΓu.Among aboveB(σ˙) is incremental surplus energy density function, for elastoplastic problem
(21)B(σ˙)=12σ˙TSσ˙,
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; F-˙ is the body force column vector; and T-˙ 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
(22)ΠHR(σ˙,u˙)=∫Ω[12β˙TPTSPβ˙-β˙TPTL𝚽uq˙+F-˙T𝚽uq˙]dΩ+∫ΓσT-˙T𝚽uq˙dΓσ+Rq˙,ΠHR(β˙,q˙)=∑n(12β˙THβ˙+β˙TGq˙-QnTq˙)+Rq˙,
where
(23)H=∫ΩPTSPdΩ,G=∫ΩPTBdΩ,QnT=∫ΩF-𝚽dΩ+∫ΓσT-𝚽dΓσ,B=L𝚽,
where R is the unbalance correction term [20, 21]; it is related to load history.According to stationary value condition of multivariation increment variation principle,(24)∂ΠHR(σ˙,u˙)∂β˙=0.We have(25)β˙=H-1Gq˙.Substituting (21) into (18) yields
(26)ΠHR(q˙)=∑n(12q˙TKq˙-QnTq˙)+Rq˙,
where
(27)K=GTH-1G.Based on the stationary condition of variation,(28)∂ΠHR(q˙)∂q˙=0.The equations for solving generalized displacement can be found:(29)Kq˙=Q-R.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:
(30)Δψ=Kq˙-ΔF.For the solution of (30) and for each load increment consider the situation existing for the rth iteration. The applied loads for the rth iteration are the residual forces Δψr-1 calculated at the end of the (r-1)th iteration according to (30). These applied loads give rise to displacement increments q˙r. The corresponding increments of the strain ε˙r 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.


## 3.1. Basic Equations of Elastoplasticity Plane Problem

For two-dimensional problem, the displacement, strain, and stress fields of a body are defined to beu0, ε0, and σ0, respectively, after random load history. When the load increases persistently, the increments of body force, surface force, and boundary displacement are b˙, t-, and u-˙, respectively, then the increments of displacement, strain, and stress are u˙, ε˙, 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(6)LTσ˙+b˙=0,(inΩ),
where L is the matrix of differential operator for 2D problems:
(7)L=[∂∂x100∂∂x2∂∂x2∂∂x1],
where σ˙ is the stress increment of point x; b˙ is the body force increment of point x.The strain-displacement relationship can be written as(8)ε˙=Lu˙,(inΩ),
where ε˙ and u˙ are strain and displacement increment of random point x, respectively.The stress-strain relationship can be expressed as(9)σ˙=Dε˙.During elasticity stage,(10)D=[D]e.Within plastic stage(11)D=[D]ep=EQ[σ22′2+2Ps-σ11′σ22′+2νPs-σ11′+νσ22′1+νσ12-σ11′σ22′+2νPsσ11′+2Ps-σ22′+νσ11′1+νσ12-σ11′+νσ22′1+νσ12-σ22′+νσ11′1+νσ12Rs2(1+ν)+2H′9E(1-ν)σ-2],whereσ- is the equivalent stress,
(12)σ-=σ112+σ222-σ11σ22+3σ122,
where σ′ is the deviatoric stress,
(13)σ11′=σ11-σ11+σ123,σ22′=σ22-σ11+σ123,Ps=2H′9Eσ-2+σ1221+ν,Rs=σ11′2+2νσ11′σ22′+σ22′2,Q=Rs+2(1-ν2)Ps,
and H′ is the hardening parameter [18].The boundary conditions can be written as(14)u˙=u-˙,(on⁡Γu),(15)n·σ˙=t-˙,(on⁡Γt),
where u˙ is the prescribed displacement increment at an arbitrary point x on the essential boundary Γu, t-˙ is the prescribed surface force increment at an arbitrary point x on the natural boundary Γt, Γ is the boundary of the domain Ω, and Γ=Γu∪Γt, Γu∩Γt=∅.


## 3.2. Incremental Hybrid Natural Element Method

Let us suppose that there existnt nodes in domain Ω, as shown in Figure 1; the displacement increment and stress increment of an arbitrary node can be expressed as
(16)u˙h(x)=∑i=1nΦiu(x)u˙i,σ˙h(x)=∑i=1nΦiσ(x)σ˙i,
where n is the total number of natural neighbor nodes around x and Φiu and Φiσ are interpolating functions of displacement increment u˙ and stress increment σ˙, which can be found from (3).Equations of (16) can be written in matrix as follows:
(17)u˙h(x)=𝚽uq˙,(18)σ˙h(x)=Pβ˙,
where(19)𝚽u=[Φ10Φ20⋯Φn00Φ10Φ2⋯0Φn],q˙=[u˙1(x1)u˙2(x1)u˙1(x2)u˙2(x2)⋯u˙1(xn)u˙2(xn)],P=[Φ100Φ200⋯Φn000Φ100Φ20⋯0Φn000Φ100Φ2⋯00Φn],β˙=[σ˙11(x1)σ˙22(x1)σ˙12(x1)σ˙11(x2)σ˙22(x2)σ˙12(x2)⋯σ˙11(xn)σ˙22(xn)σ˙12(xn)].Applying the incremental Hellinger-Reissner variation principle given by [19] to elastoplastic analysis, we have
(20)πHR=∫Ω[(F-0T+F-T)B(σ˙)-(σ0T-σ˙)(Lu˙)+σ˙T(ε0-Lu˙)+(F-0T+F-T)u˙]dΩ+∫Γσ(T-0T+T-˙T)u˙dΓσ+∫ΓuT-˙T(u˙-u-˙+u0-u-0)+T-0Tu˙dΓu.Among aboveB(σ˙) is incremental surplus energy density function, for elastoplastic problem
(21)B(σ˙)=12σ˙TSσ˙,
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; F-˙ is the body force column vector; and T-˙ 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
(22)ΠHR(σ˙,u˙)=∫Ω[12β˙TPTSPβ˙-β˙TPTL𝚽uq˙+F-˙T𝚽uq˙]dΩ+∫ΓσT-˙T𝚽uq˙dΓσ+Rq˙,ΠHR(β˙,q˙)=∑n(12β˙THβ˙+β˙TGq˙-QnTq˙)+Rq˙,
where
(23)H=∫ΩPTSPdΩ,G=∫ΩPTBdΩ,QnT=∫ΩF-𝚽dΩ+∫ΓσT-𝚽dΓσ,B=L𝚽,
where R is the unbalance correction term [20, 21]; it is related to load history.According to stationary value condition of multivariation increment variation principle,(24)∂ΠHR(σ˙,u˙)∂β˙=0.We have(25)β˙=H-1Gq˙.Substituting (21) into (18) yields
(26)ΠHR(q˙)=∑n(12q˙TKq˙-QnTq˙)+Rq˙,
where
(27)K=GTH-1G.Based on the stationary condition of variation,(28)∂ΠHR(q˙)∂q˙=0.The equations for solving generalized displacement can be found:(29)Kq˙=Q-R.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:
(30)Δψ=Kq˙-ΔF.For the solution of (30) and for each load increment consider the situation existing for the rth iteration. The applied loads for the rth iteration are the residual forces Δψr-1 calculated at the end of the (r-1)th iteration according to (30). These applied loads give rise to displacement increments q˙r. The corresponding increments of the strain ε˙r 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

Figure3 shows a cantilever beam with concentrated load, q=1 N, at the end; the length of the beam is L=8 m, and the height is h=1 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, E=105 Pa, Poisson’s ratio, ν=0.25, and the yield stress, σs=25 Pa. The linear hardening elastoplastic model is adopted with E′=0.2E.Cantilever beam subjected to a concentrated force.
(a)
(b)Figure4 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.Figure 4
The vertical displacement of the neutral axis of beam.Figure 5
Relationship between displacements of the midpoint at the end of the beam and the load.Tables1 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.Table 1
The normal stressσ11 of Gaussian quadrature points (unit: Pa).
Coordinate ofGaussian points
Analytical
HNEM
NEM
(4.267, −0.292)
−13.067
−12.988
−10.076
(5.333, −0.458)
−14.667
−14.055
−12.288
(5.876, 0.208)
5.333
5.203
4.092
(6.533, 0.042)
0.733
0.779
0.413
(7.467, 0.488)
2.933
2.764
2.111Table 2
The shear stressσ12 of Gaussian quadrature points (unit: Pa).
Coordinate of Gaussian points
Analytical
HNEM
NEM
(4.133, 0.292)
−0.990
−0.959
−4.615
(5.067, 0.208)
−1.240
−1.377
−2.574
(6.533, 0.167)
−1.333
−1.566
−2.488
(6.933, 0.042)
−1.490
−1.376
−0.824
(7.333, 0.167)
−1.333
−1.166
−1.704


### 4.2. Cylinder with Uniform Inner Pressure

Figure6 is a cylinder with uniform pressure on inner surface, the inner radius of the cylinder is a=1 m, the outer radius is b=5 m, and the uniform pressure is p=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, E=105 Pa, Poisson’s ratio, ν=0.25, and the yield stress, σs=150 Pa. The linear hardening elastoplastic model is adopted with E′=0.2E.Figure 6
A cylinder with uniform pressure on inner surface.We can only consider one quarter of the cylinder because of the symmetry, as shown in Figure7. We can arrange 99 nodes totally as shown in Figure 8. Figure 9 gives numerical solution of axial displacement of partial nodes on x1-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.Figure 7
A quarter of the cylinder with boundary conditions.Figure 8
Node distribution and Delaunay triangles of a quarter of the cylinder.Figure 9
The radial displacement atx2=0.Figure 10
Relationship between the radial displacements at the point (5.0, 0.0) and the loading.


### 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 Figure11. The length of the plate is L=10 m, and the width is h=4 m. The radius of the hole is 1 meter, and the tensile load is q=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, E=105 Pa, Poisson’s ratio, ν=0.25, and the yield stress, σs=250 Pa. The linear hardening elastoplastic model is adopted with E′=0.2E.Figure 11
The plate with a central hole under axial uniform tension.The nodes layout of the plate is given in Figure12; 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 x2=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.Figure 12
Node distribution and Delaunay triangles of plate.Figure 13
The horizontal displacement atx2=2.Figure 14
The vertical displacement atx2=2.Figure 15
Relationship between the vertical displacements at the point (5.0, 0.0) and the loading.


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

Figure3 shows a cantilever beam with concentrated load, q=1 N, at the end; the length of the beam is L=8 m, and the height is h=1 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, E=105 Pa, Poisson’s ratio, ν=0.25, and the yield stress, σs=25 Pa. The linear hardening elastoplastic model is adopted with E′=0.2E.Cantilever beam subjected to a concentrated force.
(a)
(b)Figure4 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.Figure 4
The vertical displacement of the neutral axis of beam.Figure 5
Relationship between displacements of the midpoint at the end of the beam and the load.Tables1 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.Table 1
The normal stressσ11 of Gaussian quadrature points (unit: Pa).
Coordinate ofGaussian points
Analytical
HNEM
NEM
(4.267, −0.292)
−13.067
−12.988
−10.076
(5.333, −0.458)
−14.667
−14.055
−12.288
(5.876, 0.208)
5.333
5.203
4.092
(6.533, 0.042)
0.733
0.779
0.413
(7.467, 0.488)
2.933
2.764
2.111Table 2
The shear stressσ12 of Gaussian quadrature points (unit: Pa).
Coordinate of Gaussian points
Analytical
HNEM
NEM
(4.133, 0.292)
−0.990
−0.959
−4.615
(5.067, 0.208)
−1.240
−1.377
−2.574
(6.533, 0.167)
−1.333
−1.566
−2.488
(6.933, 0.042)
−1.490
−1.376
−0.824
(7.333, 0.167)
−1.333
−1.166
−1.704


## 4.2. Cylinder with Uniform Inner Pressure

Figure6 is a cylinder with uniform pressure on inner surface, the inner radius of the cylinder is a=1 m, the outer radius is b=5 m, and the uniform pressure is p=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, E=105 Pa, Poisson’s ratio, ν=0.25, and the yield stress, σs=150 Pa. The linear hardening elastoplastic model is adopted with E′=0.2E.Figure 6
A cylinder with uniform pressure on inner surface.We can only consider one quarter of the cylinder because of the symmetry, as shown in Figure7. We can arrange 99 nodes totally as shown in Figure 8. Figure 9 gives numerical solution of axial displacement of partial nodes on x1-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.Figure 7
A quarter of the cylinder with boundary conditions.Figure 8
Node distribution and Delaunay triangles of a quarter of the cylinder.Figure 9
The radial displacement atx2=0.Figure 10
Relationship between the radial displacements at the point (5.0, 0.0) and the loading.


## 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 Figure11. The length of the plate is L=10 m, and the width is h=4 m. The radius of the hole is 1 meter, and the tensile load is q=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, E=105 Pa, Poisson’s ratio, ν=0.25, and the yield stress, σs=250 Pa. The linear hardening elastoplastic model is adopted with E′=0.2E.Figure 11
The plate with a central hole under axial uniform tension.The nodes layout of the plate is given in Figure12; 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 x2=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.Figure 12
Node distribution and Delaunay triangles of plate.Figure 13
The horizontal displacement atx2=2.Figure 14
The vertical displacement atx2=2.Figure 15
Relationship between the vertical displacements at the point (5.0, 0.0) and the loading.


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



---

*Source: 979216-2014-03-30.xml*