The Incremental Hybrid Natural Element Method for Elastoplasticity Problems

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: 𝑇_𝑖=𝐱∈𝑅²∣𝑑𝐱,𝐱_𝑖<𝑑𝐱,𝐱_𝑗,∀𝑗≠𝑖, (1) where 𝑑(𝐱,𝐱_𝑖) is distance between any point 𝐱 and the discrete point 𝐱_𝑖.

Figure 1 

Discrete model of region Ω and Voronoi cell about 𝐱𝑖.

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

Figure 2 

The 2nd order Voronoi cell about 𝐱.

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: 𝑓(𝐱)=^𝑛_𝑖=1Φ_𝑖(𝐱)𝑓_𝑖, (2) 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 Φ_𝑖=𝐴_𝑖(𝐱)𝐴(𝐱),𝐴(𝐱)=^𝑛_𝑗=1𝐴_𝑗(𝐱), (3)(4) 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: 𝑢^ℎ(𝐱)=^𝑛_𝑖=1Φ_𝑖(𝐱)𝑢_𝑖 (5) 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 𝐋^T𝝈̇+𝐛̇=0,(inΩ), (6) where 𝐋 is the matrix of differential operator for 2D problems: 𝜕𝐋=𝜕𝑥₁00𝜕𝜕𝑥₂𝜕𝜕𝑥₂𝜕𝜕𝑥₁, (7) where 𝝈̇ is the stress increment of point 𝐱; 𝐛̇ is the body force increment of point 𝐱.

The strain-displacement relationship can be written as 𝜺̇=𝐋𝐮̇,(inΩ), (8) where 𝜺̇ and 𝐮̇ are strain and displacement increment of random point 𝐱, respectively.

The stress-strain relationship can be expressed as 𝝈̇=𝐃𝜺̇. (9)

During elasticity stage, 𝐃=[𝐷]_𝑒. (10)

Within plastic stage𝐃=[𝐷]_𝑒𝑝=𝐸𝑄𝜎^′2₂₂+2𝑃_𝑠−𝜎^′₁₁𝜎^′₂₂+2𝜈𝑃_𝑠−𝜎^′₁₁+𝜈𝜎^′₂₂𝜎1+𝜈¹²−𝜎^′₁₁𝜎^′₂₂+2𝜈𝑃_𝑠𝜎^′₁₁+2𝑃_𝑠−𝜎^′₂₂+𝜈𝜎^′₁₁𝜎1+𝜈¹²−𝜎^′₁₁+𝜈𝜎^′₂₂𝜎1+𝜈¹²−𝜎^′₂₂+𝜈𝜎^′₁₁𝜎1+𝜈¹²𝑅_𝑠2()+1+𝜈2𝐻^′()9𝐸1−𝜈𝜎², (11)

where 𝜎 is the equivalent stress, 𝜎=𝜎²₁₁+𝜎²₂₂−𝜎₁₁𝜎₂₂+3𝜎²₁₂, (12) where 𝜎^′ is the deviatoric stress, 𝜎^′₁₁=𝜎₁₁−𝜎₁₁+𝜎₁₂3,𝜎^′₂₂=𝜎₂₂−𝜎₁₁+𝜎₁₂3,𝑃_𝑠=2𝐻^′9𝐸𝜎²+𝜎²₁₂,𝑅1+𝜈_𝑠=𝜎^′2₁₁+2𝜈𝜎^′₁₁𝜎^′₂₂+𝜎^′2₂₂,𝑄=𝑅_𝑠+21−𝜈²𝑃_𝑠, (13) and 𝐻^′ is the hardening parameter [18].

The boundary conditions can be written as 𝐮̇=𝐮̇,onΓ_𝑢,𝐧⋅𝝈̇=𝐭,̇onΓ_𝑡, (14)(15) 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 𝑢̇^ℎ(𝐱)=^𝑛_𝑖=1Φ^𝑢_𝑖(𝐱)𝑢̇_𝑖,𝜎̇^ℎ(𝐱)=^𝑛_𝑖=1Φ^𝜎_𝑖(𝐱)𝜎̇_𝑖, (16) 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: 𝑢̇^ℎ(𝐱)=𝚽^𝐮̇𝜎̇𝐪,_ℎ(𝐱)=𝐏𝜷̇, (17)(18) where𝚽^𝐮=Φ₁0Φ₂0⋯Φ_𝑛00Φ₁0Φ₂⋯0Φ_𝑛,𝐪̇=𝑢̇₁𝐱₁𝑢̇₂𝐱₁𝑢̇₁𝐱₂𝑢̇₂𝐱₂⋯𝑢̇₁𝐱_𝑛𝑢̇₂𝐱_𝑛,Φ𝐏=¹00Φ₂00⋯Φ_𝑛000Φ₁00Φ₂0⋯0Φ_𝑛000Φ₁00Φ₂⋯00Φ_𝑛,𝜷̇=𝜎̇₁₁𝐱₁𝜎̇₂₂𝐱₁𝜎̇₁₂𝐱₁𝜎̇₁₁𝐱₂𝜎̇₂₂𝐱₂𝜎̇₁₂𝐱₂⋯𝜎̇₁₁𝐱_𝑛𝜎̇₂₂𝐱_𝑛𝜎̇₁₂𝐱_𝑛. (19)

Applying the incremental Hellinger-Reissner variation principle given by [19] to elastoplastic analysis, we have 𝜋_HR=_Ω𝐵(𝜎̇)−𝝈^0𝑇−𝝈̇(𝐋𝐮̇)+𝝈̇^𝑇𝜺^𝟎−𝐋𝐮̇+𝐅^0𝑇+𝐅^𝑇𝐮̇+𝑑Ω_Γ𝜎𝐓^0𝑇+𝐓̇^𝑇𝐮̇𝑑Γ_𝜎+_Γ𝑢𝐓̇‌^𝑇𝐮̇−𝐮̇+𝐮⁰−𝐮⁰+𝐓^0𝑇𝐮̇𝑑Γ_𝑢. (20)

Among above 𝐵(𝜎̇) is incremental surplus energy density function, for elastoplastic problem 𝐵(𝜎̇)=12𝝈̇^𝑇𝐒𝝈̇, (21) 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 Π_HR(𝜎̇,𝑢̇)=_Ω12𝜷̇^𝑇𝐏^𝑇𝐒𝐏𝜷̇−𝜷̇^𝑇𝐏^𝑇𝐋𝚽^𝐮𝐪̇+𝐅̇^𝑇𝚽^𝐮𝐪̇+𝑑Ω_Γ𝜎𝐓̇‌^𝑇𝚽^𝐮𝐪̇𝑑Γ_𝜎̇Π+𝐑𝐪,_HR𝛽̇,𝑞̇=_𝑛12𝜷̇^𝑇𝐇𝜷̇+𝜷̇^𝑇𝐆𝐪̇−𝐐^𝑇_𝑛𝐪̇+𝐑𝐪̇, (22) where 𝐇=_Ω𝐏^𝑇𝐒𝐏𝑑Ω,𝐆=_Ω𝐏^𝑇𝐐𝐁𝑑Ω,_𝑇𝑛=_Ω𝐅𝚽𝑑Ω+_Γ𝜎𝐓𝚽𝑑Γ_𝜎,𝐁=𝐋𝚽, (23) where 𝐑 is the unbalance correction term [20, 21]; it is related to load history.

According to stationary value condition of multivariation increment variation principle, 𝜕Π_HR(𝜎̇,𝑢̇)𝜕𝜷̇=0. (24)

We have 𝜷̇=𝐇⁻¹𝐆𝐪̇. (25)

Substituting (21) into (18) yields Π_HR𝑞̇=_𝑛12𝐪̇^𝑇𝐊𝐪̇−𝐐^𝑇_𝑛𝐪̇+𝐑𝐪̇, (26) where 𝐊=𝐆^𝑇𝐇⁻¹𝐆. (27)

Based on the stationary condition of variation, 𝜕Π_HR𝑞̇𝜕𝐪̇=0. (28)

The equations for solving generalized displacement can be found: 𝐊𝐪̇=𝐐−𝐑. (29)

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: Δ𝜓=𝐊𝐪̇−Δ𝐅. (30)

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 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, 𝑞=1 N, at the end; the length of the beam is 𝐿=8 m, and the height is ℎ=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, 𝐸=10⁵ Pa, Poisson’s ratio, 𝜈=0.25, and the yield stress, 𝜎_𝑠=25 Pa. The linear hardening elastoplastic model is adopted with 𝐸^′=0.2𝐸.

(a)

(b)

Figure 3 

Cantilever beam subjected to a concentrated force.

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.

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.

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.

Table 1 

The normal stress 𝜎₁₁ of Gaussian quadrature points (unit: Pa).

Table 2 

The shear stress 𝜎₁₂ of Gaussian quadrature points (unit: Pa).

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 𝑎=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⁵ Pa, Poisson’s ratio, 𝜈=0.25, and the yield stress, 𝜎_𝑠=150 Pa. The linear hardening elastoplastic model is adopted with 𝐸^′=0.2𝐸.

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

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 at 𝑥2=𝟎.

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 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⁵ Pa, Poisson’s ratio, 𝜈=0.25, and the yield stress, 𝜎_𝑠=250 Pa. The linear hardening elastoplastic model is adopted with 𝐸^′=0.2𝐸.

Figure 11 

The plate with a central hole under axial uniform tension.

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 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 at 𝑥2=2.

Figure 14 

The vertical displacement at 𝑥2=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.

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