Elastoplastic Large Deformation Using Meshless Integral Method

In this paper, the meshless integral method based on the regularized boundary integral equation [1] has been extended to analyze the large deformation of elastoplastic materials. The updated Lagrangian governing integral equation is obtained from the weak form of elastoplasticity based on Green-Naghdi’s theory over a local sub-domain, and the moving least-squares approximation is used for meshless function approximation. Green-Naghdi’s theory starts with the additive decomposition of the Green-Lagrange strain into elastic and plastic parts and considers a J2 elastoplastic constitutive law that relates the Green-Lagrange strain to the second Piola-Kirchhoff stress. A simple, generalized collocation method is proposed to enforce essential boundary conditions straightforwardly and accurately, while natural boundary conditions are incorporated in the system governing equations and require no special handling. The solution algorithm for large deformation analysis is discussed in detail. Numerical examples show that meshless integral method with large deformation is accurate and robust.


Introduction
Over the past two decades the meshless methods have attracted much attention owing to their advantages in adaptivity, higher degree of continuity in the solution field, and the capability to handle moving boundaries and changing geometry.In the meshless method, the concept of an element is eliminated.The model geometry consists of a distribution of nodes over the problem domain, and the approximate solution is constructed entirely based on these nodes.Most development in meshless methods to date has been focused mostly on linear elasticity.Research in large deformation plasticity using meshless methods is only currently gaining attention.
For the formulation of elastoplastic theory at large deformation, Green-Naghdi's theory [2,3] begins with the additive decomposition of the Green-Lagrange strain into elastic and plastic parts as . On the other hand, Lee's theory [3,4] considers the multiplicative decomposition of deformation gradient F into an elastic F e and plastic part F p as e p F F F  .Both theories are formulated on the basis of fundamental laws of continuum mechanics.According to Chiou et al. [3], Green-Naghdi's theory is more flexible because it can be applied to either isotropic or anisotropic materials and the computing procedure involved is relatively straight-forward.In general, theories of large deformation plasticity have to be implemented numerically since it is difficult to obtain closed form solutions to practical problems involving large deformation.Researchers such as Chiou et al. [5], Lee [6], and Hu [7], have used FEM to solve large deformation plasticity problems.Some researchers have used meshless methods to analyze large deformation with elastoplasticity.Belytschko et al. [8] proposed a 3D Element-Free Galerkin (EFG) method intended for dynamic problems with geometric and material nonlinearities solved with explicit time integration.Rossi and Alves [9] applied a modified Element-Free Galerkin method to large deformation processes.The proposed EFG method enables direct enforcement of essential boundary conditions, and the plasticity model assumes a multiplicative decomposition of the deformation gradient into an elastic and plastic part and allows for nonlinear isotropic hardening.Chen et al. [10] and Eskandarian et al. [11] formulated the governing equations for rate-independent large strain plasticity in the framework of meshless method and used the method to simulate high-speed impact.Xiong et al. [12] used meshless local Petrov-Galerkin method to model large deformation of micro-electronic mechanical devices (MEMS) by using local radial point interpolation (RPIM) approximation to construct shape function.The system equation is obtained based on the Total Lagrangian (TL) approach.Hu et al. [13] used meshless local Petrov-Galerkin method for large deformation contact analysis of elastomeric components.A nonlinear formulation of meshless local Petrov-Galerkin finite-volume mixed method was developed by Han et al. to analyze static and dynamic large deformation problems [14].Li et al. [15] developed a coupled finite element and meshless local Petrov-Galerkin method to analyze large deformation problems.Gu et al. [16,17] employed meshless local Kriging method to analyze large deformation problems.A comparison between Total Lagrangian and Updated Lagrangian approaches with an adaptive local meshless formulation was given by Gu in [18].Gun et al. [19] used a meshless formulation of the Euler-Nernouli theory with non-linear strain-displacement relation which covers both axial and bending deformations to predict the large deformation behavior of fabrics.Li and Lee [20,21] employed an adaptive meshless method to solve mechanical contact problems which incorporated a sliding line algorithm with penalty method to handle contact constrains.In [22], Zhu et al. used meshless local natural neighbor interpolation method to analyze 2D elastoplastic large deformation in conjunction with Updated Lagrangian formulation.A Galerkin SPH method was utilized by Wang to solve large deformation fracture problems [23].Reproducing kernel particle methods were used to analyze large deformation problems in [24][25][26].In [27], Li et al. used an element-free Galerkin method to solve die forging problems.Quak et al. [28] applied both SPH and EFG meshless methods to analyze extraction problems.
The authors have developed a meshless integral method for linear elasticity [1] and later extended it to elastoplasticity for small deformation [29].This meshless integral method is truly meshless and does not require a background mesh for integration.In the present work, we extend the meshless integral method to large deformation elastoplasticity.The governing integral equation is obtained from the weak form based on Green-Naghdi's large deformation elastoplasticity theory over a local subdomain, and moving least-squares approximation is used for meshless function approximation.The constitutive law uses A J 2 elastoplastic flow theory based on von Mises yield criterion with isotropic hardening in which the Green-Lagrange strain is related to the second Piola-Kirchhoff stress.The fixed point iteration, because of its simplicity in implementation, is used to solve the nonlinear equations arising from large deformation elastoplasticity.A number of numerical tests were carried out for model verification.It was found that the meshless results are in excellent or satisfactory agreement with either closed form solutions or FEM results, indicating that the large deformation elastoplasticmeshless integral method based on the regularized integral equation is accurate and robust.An application of the current method to the metal forming process has been reported separately [30].
This paper is organized as follows: In Section 2, the regularized local boundary integral equation is derived, and the subtraction method is used to remove the strong singularity that is present in the initial local boundary integral equation.In Section 3, large deformation elastoplastic constitutive equations are presented.Section 4 describes the moving least-squares approximation (MLSA) for approximating the displacement, strain, and stress fields over the problem domain.The meshless implementation of the regularized boundary integral equation using MLSA and the treatment of weak singularity are presented in Section 5. Section 6 discusses the enforcement of essential and natural boundary conditions.The solution algorithm for solving the nonlinear system equations is discussed in Section 7. Numerical examples are presented in Section 8 to assess the accuracy and effectiveness of the method.Discussion and conclusions from this study are given in Section 9.

Regularized Local Boundary Integral Equation Using Subtraction Method
In this work, the updated Lagrangian coordinate system is used in which the stress directions are referred to the last known equilibrium state.Consider a two-dimensional body for a given load increment as shown in Eulerian coordinates, x k (k = 1, 2), and Lagrangian coordinates, X K (K = 1, 2), are employed.We use upper case letter and upper case subscript to indicate that the variable of interest refers to the reference state, and use lower case letter and lower case subscript to indicate that the variable of interest refers to the current state.Within the updated Lagrangian framework, the following procedures are implied: 1) During each load increment, all field variables are defined with respect to the state at the start of this load increment (reference configuration); 2) At the end of this load increment, field variables are updated with respect to the state at the end of this load increment (current configuration), and the current configuration will be the reference configuration for the next load increment.
In large deformation problems, various stress measures can be defined.The most commonly used stresses are: 1) Cauchy stress (true stress) tensor σ which expresses the stress relative to the current configuration; 2) First Piola-Kirchhoff stress (nominal stress) tensor P which relates forces in the current configuration with areas in the reference configuration; 3) Second Piola-Kirchhoff stress tensor T which relates forces in the reference configuration to areas in the reference configuration, where the force in the reference configuration is obtained via a mapping that preserves the relative relationship between the force direction and the area normal in the current configuration.The relationships between these stresses are where F is the deformation gradient and defined as

F
, and J is the determinant of F. The second Piola-Kirchhoff stress tensor is symmetric, and in the current work is linked to the Green-Lagrange strain in the elastoplastic constitutive law.For a large deformation elastoplastic body represented by a reference domain Ω 0 with boundary Γ 0 , the governing elastoplasticity equations are as follows: where τ (bold face denotes vectors or tensors) is the stress measure defined by T LJ is the second Piola-Kirchhoff stress, G is the transformation tensor between the reference configuration and the current configuration, on on here, u represents the prescribed displacement on u  , T represents the prescribed traction on t  ; n is the outward unit normal to the boundary; where the notation  in this paper is used to denote a node (e.g., a indicates node a, b indicates node b, etc.) in order to reserve the usual subscripts, I, J, etc., for denoting degree of freedom (DOF) components, 0 0 a    is a spherical (circular in 2D) sub-domain related to node a , a Y is the position vector of node a which is also called a source point, X is the integration or field point which may or may not coincide with a node, and I g is the test function.In the following, the functional dependence on X, i.e. "(X)", will be dropped for brevity when no ambiguity is caused.In this work, following [31], we use a special test function defined by where J e represents the J-th component of a unit force vector, JI u   is the special test function, given by [31]: is the shear modulus.Note that the special test function is defined in the reference state.The associated traction, for linear elastic behavior, is The special test function IJ u   has the property that it vanishes on the boundary of a spherical 0 a  .With IJ u   as the test function, application of integration by parts to (4) twice leads to:  Equation (8) in its current form cannot be used directly in numerical calculation because when a Y is a boundary node, the equation contains strong singularity (1/r type in line integral) in the traction term IJ t   .For the local inte- gral equation approach to be a valid numerical method, the strong singularity must be handled appropriately.
The subtraction technique is employed in the present study to remove the strong singularity; the technique for small deformation has been presented in 4 [1].For simplicity in implementation, the local sub-domain is always chosen, in the reference state, as a sphere or part of a sphere centered on a node.the section where the traction is prescribed.
Using subtraction method, we obtain, in the limit of 0   , the following: , d The integration of Here  Substituting the above expressions into (11) leads to: In ( 14), the subtraction method has been used in the last term on the right hand side.When the field point X approaches the source node a Y ,   a I I u x u  tends to zero which removes the strong singularity and makes the integral numerically integrable.All other terms in (14) are regular or weakly singular for which special integration quadrature gives convergent and accurate results.The regularized Equation ( 14) holds for any source node a Y , either inside the domain or on the boundary.In the current meshless integral method, both boundary and interior source nodes are used, and the moving leastsquares approximation is employed for approximating the solution field, as presented later.

Constitutive Equation for Elastoplasticity with Large Deformation
In this research, the Green-Naghdi's theory with von Mises yield criterion, associated flow rule, and isotropic strain hardening is used to model the material behavior.The Green-Naghdi's theory is attractive because it can be applied to either isotropic or anisotropic materials and the computing procedure involved is relatively straightforward.
Green-Naghdi's theory [2,3] begins with the additive decomposition of the Green-Lagrange strain increment into the elastic and plastic parts as The von Mises yield criterion has the following form: f is the yield function constructed in the space of second Piola-Kirchhoff stress, IJ T is the component of deviatoric stress of the second Piola-Kirchhoff stress, and characterizes the hardening of the material.For linear hardening, where 0 Y is the initial yield stress, p H is the hardening modulus, and p  is the effective plastic strain.p H is given by the following equation where T E is elastoplastic tangent modulus and Y E is the elastic Young's modulus from uniaxial tension test.
The plastic strain increment p IJ dE is expressed by an associated flow rule as: where L is the loading criteria and given by the following equation: [C] is elastic stiffness matrix which works for both plane strain and plane stress.
The positive scalar H is expressed as: Using Equations ( 18) and ( 19), we obtain the plastic strain increment: With p IJ dE , the increment of effective plastic strain is computed as follows: where 3 2 is the equivalent stress of the current stress state.
The constitutive equations are summarized as follows: Where  and  are Lame constants.

Moving Least-Squares Approximation
In the finite element method, the coupling between the nodes is accomplished through the use of shape functions, defined locally over each element, which interpolate the solution field from nodal values.For a meshless method, the absence of elements excludes the use of such shape functions and therefore, and a different local approximation scheme based on nodal values but independent of any elements needs to be devised.In this work, we have chosen to exploit the non-interpolative moving least-squares approximation (MLSA) scheme because of its high accuracy and the ease with which it can be extended to n-dimensional problems.
Consider a reference domain 0  that contains n nodes: Following [29], we define a support domain for node As introduced previously, the sub-domain  The moving least-squares approximant h u to a function u is defined by: The two vectors p and c are both functions of the spatial coordinates: . p is a complete monomial basis of m terms (e.g., in 2D, m = 3 for a linear basis, and 6 for a quadratic basis), c is a coefficient vector which is determined by minimizing a weighted discrete L 2 -norm: where ˆa u is the fictitious nodal displacement that ap- proximates the value of u at node a Y , and the upper limit of summation, x N , is the total number of nodes in the domain of definition of point X.The matrices P, W and û are defined by Minimization of (28) leads to: (32) where: The MLSA for a function exists only when

 
A X is non-singular.A necessary condition for a well-defined MLSA is that for each sample point The influence of the choices of the basis functions and the weight functions on the behavior and the quality of the shape function has been discussed in [1].In this work, we use linear, quadratic monomial basis, and spline weight functions, defined as follows: l , which makes its use simple.It is noted that MLSA is non-interpolative, and there is a difference between the nodal value of the MLSA approximant h u and the fictitious nodal displacement ˆa u .For brevity, the subscript h in h u will be omitted in the remainder of this paper.

Meshless Implementation
We now apply the MLSA to the integral Equation ( 14) to establish the meshless implementation.The shape function, as we have defined it, gives: where x N is the total number of nodes in the domain of definition of point X.
The related traction term J T is where  

B
matrices as: Combining ( 42) with ( 39) and ( 41), we can express the traction in terms of the shape functions as follows: here ep C is the elastoplastic stiffness matrix.With the above discretization and the boundary conditions that  , Equation ( 14) becomes (there is a summation on b and J but not on a and I): , (47) and the upper limit of summation, y N , is the total num- ber of nodes in the domain of influence of node a Y .The third to sixth terms of Equation ( 47) are the nonlinear terms because of the large deformation and elastoplasticity.In all numerical examples tested in this work, the body force term is zero.
Owing to the subtraction technique, the singularity in the third integral on the right hand side of (45) as X ap- is cancelled by the term in (46), and similarly the seventh integral on the right hand side of (47) is also regularized.Even though the subtraction technique removes the strong singularity, the integrands in the first integral on the right hand side of (45) and the second integral on the right hand side of (47) still contain the weakly singular ln (r) term.The logarithmic singularity is integrable, but the accuracy of ordinary Legendre-Gauss integration is poor.We found that the special integration scheme for the logarithmic singularity [22], which is reproduced in Table 1 for completeness, achieves excellent numerical accuracy [1].The remaining integrals of Equation (45) and Equation ( 47) are regular for which standard quadrature can be used with good accuracy.
In our numerical examples, the numbers of integration points were as follows: 8 integration points for any integral along a straight line, and 8 × 8 integration points for any integral over a sub-domain.For regularized integrals, the usual Legendre-Gauss integration was used.For integrals containing logarithmic singularity, the special integration of 8 integration points [32] as listed in Table 1 was used.

Enforcement of Boundary Conditions
Appropriate boundary conditions need to be enforced in order to solve the simultaneous Equations (44).In mesh-Table 1. Special numerical integration for functions containing logarithmic singularity (8 integration points).
A number of collocation methods have been developed.The direct collocation method [33] used the condition to replace the row of the discretized weak form equation corresponding to the degree of freedom with prescribed displacement a I u .This is actually inconsistent with the assumption of MLSA since the fictitious nodal displacement ˆa I u is generally not equal to the approximated displacement value.
A generalized collocation method uses as the collocation condition which was shown to yield more accurate results [44].Wagner and Liu [45] developed a corrected collocation method which restores the consistency of the weak form and enhances convergence.Wu and Plesha [46] proposed a boundary flux collocation method to enforce the boundary conditions exactly.
Generally, there are two types of discretization in meshless methods: 1) Galerkin based method over the global domain, which is used in EFG [33,42,47,48], clouds [41], RKPM [49][50][51][52]; and 2) local weak form over multiple local domains, which is employed in [31,[53][54][55].For Galerkin based method over the global domain, the n system equations are obtained by applying the weak form over the global domain once, hence all equations must hold simultaneously in order to maintain consistency of the weak form.Replacing a row in the matrix equation by (49), which contains a linear combination of DOFs rather than dictating the single value of a constrained DOF, sacrifices the consistency of the weak form and compromises the solution accuracy.
For local weak form over multiple local domains, each equation is obtained by applying the weak form over a particular local domain, and the weak form needs to be applied n times for a problem with n DOFs.Consequently, the generalized collocation method with (49) can be directly used to enforce essential boundary conditions.Because each system equation is independent of the rest, replacing the equation corresponding to a constrained DOF by (49) will not cause any inconsistency in the weak form.
As the current meshless integral method uses local weak form over multiple sub-domains, the generalized collocation method (49) is applicable.For a DOF with essential boundary condition, we simply use the condition I I u u  rather than applying the integral Equation (14).In numerical implementation, this means replacing the governing equation corresponding to the prescribed DOF with generalized collocation condition (49).Note that the governing Equation ( 14) holds for interior nodes as well as boundary (including corner) nodes, therefore Equation ( 49) is applied to nodes on smooth boundaries as well as at corners.Our numerical tests show that the generalized collocation method (49) for enforcing essential boundary conditions works well for the meshless integral method.
For the natural boundary condition I I t t  , no special treatment is needed.The prescribed traction is directly used in the second integral in Equation (47).
Although theoretically essential boundary conditions can be converted into the equivalent natural boundary conditions and the solution remains unchanged, we observed that the meshless method is always more accurate when essential boundary conditions are used.
After the boundary conditions are enforced, the governing equations can be written as where and the upper limit of summation, y N , is the total num- ber of nodes in the domain of influence of node a Y .Detailed expressions of various terms are given in Equations (45) to (47).

Solution Algorithm for Elastoplasticity with Large Deformation
In this work, we will use the fixed point iteration to solve the governing equations which has the following advantages: the derivative of the stiffness matrix is not needed, and the implementation is relatively easy.With this algorithm, for each load increment, the stiffness matrix is computed only once in the first iteration.
The material considered in this study is rate independent, and therefore the actual loading rate has no effect on the solution.It is, however, convenient to introduce a time parameter to describe the loading history.A piecewise linear function is used in the program to describe the history of applied loads and is implemented in a table form in the input data.Each linear segment is referred to as a load step.For a loading history with N load steps, the corresponding data block takes the following form: Since the elastoplastic integral equations are strongly nonlinear, loads have to be applied in small increments, and iterations are needed to achieve convergence within each load increment.For the first load step 1 P at a level sufficient to cause yielding in the material, the solution is first obtained by assuming elastic behavior.The load and solution are then scaled linearly to initial yielding by a scaling factor r where 1 r  .The remaining portion of load from 1 rP to   1 1 r P  is then divided into M increments.Within each load increment, the solution is iterated until convergence.For subsequent load steps, the load is incremented in a similar manner.
Throughout this section, we use       m var i to denote a variable   var (which may be displacement, stress or strain) for m-th load increment at i-th iteration, and use     var m to denote the converged solution for m-th load increment.

Solution Algorithm
1) Divide the load into N load steps, set the load increment number M for each load step and maximum iteration number Imax; initialize the load step index n, load increment index m, and iteration index i to 1.
2) For load step index n and load increment index m, the convergent state at (m -1)-th load increment will be used as the reference configuration.Set the initial displacement relative to the reference configuration , and compute the stiffness matrix [K] based on the reference configuration.Set the second Piola- , where    is the convergent Cauchy stress at (m -1)-th load increment.
where   ΔR is the load increment from (m -1)-th load increment to m-th load increment.
the nonlinear terms because of large deformation and elastoplasticity corresponding to the third to sixth terms of the Equation ( 47) for m-th load increment at (i -1)-th iteration.
4) The following convergence criterion is used to terminate the iteration for each load step: f  is a predefined tolerance, usually in the order of where J is the determinant of the deformation gradient F which is defined as follows: then go to step 5. b) Otherwise, increase the iteration index by one and check if i is greater than Imax.If yes, the program cannot converge for the preset Imax, exit program execution; otherwise, compute the second Piola-Kirchhoff stresses at all Gaussian points for all nodes and go to step 3.
5) Increase the load increment index m by one and check if m is greater than load increment number M. If yes, exit the program; otherwise, go to step 2. The integrals for next load increment are to be performed over the deformed geometry of converged solution of current load increment, and the converged current configuration will be the reference configuration for the next load increment, hence the updated Lagrangian formulation.
As shown in the previous section, the nonlinear terms (the third to sixth terms of the Equation ( 47)) for m-th load increment at (i -1)-th iteration are needed in order to calculate the internal force       The integration is usually performed using Gaussian numerical integration.Therefore, the stress state, along with the stress correction term, is computed at all Gaussian points in each iteration step.The following section describes the algorithm to compute the stress at each Gaussian point.

Procedure for Computing Stresses at Each Gaussian Point
1) After the solution is converged for (m -1)-th load increment, the corresponding Cauchy stress using the second Piola-Kirchhoff stresses which are basic stress measure obtained from solution is computed.The Cauchy stress obtained is then used as the initial value of the 2) Determine the loading state 2.1) If EPF = 1, the Gaussian point is in an elastoplastic state in previous load step.Compute the loading criterion function L using Equation (18).Use parameter r to denote the scaling factor such that If L > 0, let r = 0, plastic loading occurs; if L ≤ 0, let r = 1, EPF = 0, compute the yield function after the trial stress increment is applied r = 1 and EPF = 0, the point remains in the elastic state; if 0 f  , let EPF = 1, the point enters into elastoplastic state, determine scaling factor r using (63).Update the stress at this point by 2) If EPF = 0, the Gaussian point is in an elastic state in the previous load step.Compute the yield function after the trial stress increment is applied r = 1 and EPF = 0, the point remains in the elastic state; if 0 f  , let EPF = 1, the point enters into elastoplastic state.Determine scaling factor r using (63).Update the stress at this Gaussian point using Equation (64).
3) Compute the sub-increment of strain   where N is an integer.Integrate numerically to compute sub-increment of stress   ij T   with n looping from 1 to is used to denote the second Piola-Kirchhoff stress increment at the end of n-th iteration and let 4) Update the variables:

The Computation of Nonlinear Terms
For the third term in Equation ( 47), for each Gaussian point, and the third term in Equation ( 47) is computed accordingly.
and the fourth term of Equation ( 47) is computed accordingly.
As for the fifth and sixth terms of Equation ( 47), 0 LK T is the initial value of the second Piola-Kirchhoff stress of the current load increment.
is the increment of the second Piola-Kirchhoff stress of the current iteration with respect to 0 LK T .These two terms can be computed easily.

Numerical Examples
This section presents the numerical solutions to several large deformation elastoplastic problems using the meshless integral method.The tests include the uniaxial tension tests, shear tests, rotation tests, and punch test.For large deformation elasticity cases, closed form solutions are used for comparison.For large deformation elastoplasticity cases, the finite element results obtained using commercial software, ANSYS (http://www.ansys.com/), is used as the basis for comparison.Since FEM results are approximate solutions, convergence study was carried out for each FEM model in which the mesh was refined progressively until the global L 2 -norm between two successive meshes was within at least 4 10  .With such tight convergence, we regarded the FEM solutions as practically exact.The FEM models used for comparison, therefore, are much finer than the corresponding meshless models, but care was taken to ensure that all nodes in a meshless model coincide exactly with a subset of nodes in the corresponding FEM model in order to facilitate direct comparison.
The guidelines for the selection of the size of subdomain and the size of support domain a w l have been discussed in detail in [1].To summarize briefly, the sizes are set to be proportional to the minimum nodal distance for a node, minDist, defined as the minimum distance between the node and it neighbors, with proportional coefficients being of the order of 1.A global error indicator, the L 2 -norm error in displacement, defined by is used as a measure of the overall performance of the numerical method, where num denotes the numerical result, and ex denotes the exact solution.
In all numerical examples presented below, for elastic material response the properties of AISI 1020 steel with Young's modulus E = 203 GPa and Poisson's ratio 0.3 under plain strain condition were used.If the material behavior was elastoplastic, then after yielding, bilinear stress-strain behavior was assumed with a yield stress of 260 MPa, and an elastoplastic tangent modulus of 1 GPa.

Uniaxial Tension Tests
The uniaxial tension patch test geometry is a 1 × 1 m 2 square plate shown in Figure 7(a).Three meshless models, shown in Figures 7(b) to (d), were simulated in which spline weight function was used.The left edge was constrained from moving in X 1 direction and traction free in X 2 direction; the bottom edge was constrained from moving in X 2 direction and traction free in X 1 direction; the top edge was traction free in both X 1 and X 2 directions; and the right edge was traction free in X 2 direction and had various prescribed displacement in X 1 depending on the test case.The uniaxial tension patch tests were investigated for both elastic only material and elastoplastic material.
For the first elastic only case, a prescribed displacement U 1 of 0.1 m in X 1 direction was applied to the right edgeand the effect of load increments was investigated.For one and two load increments, closed form solutions were obtained for comparison.Meshless results for the same load increments were identical with closed form results, as shown in Tables 2 and 3.For larger load increments, closed form solutions were quite laborious and therefore not pursued.Meshless results for up to 12 load increments are plotted in   result of the large prescribed displacement corresponding to a 10% engineering strain and the enforced elastic only behavior.For comparison, FEM results for the same load increments are σ 11 = 21261 MPa, σ 33 = 6378 MPa, and b = -0.04027m.From this example, a guideline for determining load increment size was adopted for subsequent simulations: the largest engineering strain increment in a load increment should be in the range of 0.005 to 0.01.For the second elastic only case, a prescribed displacement U 1 of 0.5 m in X 1 direction was applied to the right edge using 20, 40, and 60 load increments.Difference in results using 40 and 60 load increments was small.Figures 9 and 10 show the variation of σ 11 and σ 33 , respectively, with applied displacement U 1 from the meshless model using 60 load increments.Corresponding noticeable.
For elastoplastic material behavior, since closed form results are not available, FEM results were used as the basis for comparison.For the first elastoplastic case, a prescribed displacement U 1 of 0.1 m in X 1 direction was applied to the right edge under different number of load increments.Figures 11(a It is noted that the current meshless method is based on Green-Naghdi's theory, while FEM is based on E. H. Lee's theory [2].The differences between meshless results and FEM results show that large deformation elastoplasticity theories may have some influence on the numerical results.
For the second elastoplastic case, a prescribed displacement U 1 of 0.5 m in X 1 direction was applied to the right edge using 60 load increments.Figures 12 and 13 show the variation of σ 11 and σ 33 , respectively, with applied displacement U 1 from the meshless model and FEM model.Good agreement between the meshless results and FEM results is evident.

The Shear Tests
The shear patch test geometry is a 1 × 1 m 2 square plate shown in Figure 14.Three meshless models, shown in Figures 14(b) to (d), with spline weight function were simulated.The motion of the square was described by .Three values of  , 0.1, 0.5 and 1.0, were simulated.For all three models, prescribed displacements were applied to all edges of the plate.For the shear tests, both elastic only material and elastoplastic material were investigated.For elastic only material behavior, analytical results in terms of the Cauchy stress using Jaumann rate for this problem are given as [56]: The constitutive equation in terms of the Jaumann rate is given by: where D is the rate-of-deformation tensor and 0.1

 
W is the spin tensor, and the superscript, J, denotes that the material constants are used with the Jaumann rate.
For the first elastic only shear test with 0.1   , the For elastoplastic material behavior, FEM results were used as the basis for comparison.For 0.1

 
, the L 2 -norms between the meshless and FEM results for the 9 node, 25 node regular, and 25 node irregular models are 3.7 × 10 -13 , 4.6 × 10 -13 , and 6.4 × 10 -14 respectively.The shear stress for all three models is σ 12 = 182.32MPa which is in good agreement with the FEM solution of σ 12 = 182.60MPa.
For elastoplastic behavior with 0.5

 
, the L 2 -norms between the meshless and FEM results for the 9 node, 25 node regular, and 25 node irregular models are 5.2 × 10 -13 , 1.8 × 10 -13 , and 9.4 × 10 -11 respectively.The shear stress for all three models is σ 12 = 316.52MPa which is again in good agreement with the FEM solution of σ 12 = 316.89MPa. Figure 16 shows the variation of shear stress σ 12 with  from the meshless model, together with FEM results.The meshless results match the FEM results almost identically.

The Rigid Body Rotation Tests
The rigid body rotation test geometry is a 1 × 1 m 2 square plate shown in Figure 17(a).Figure 17

The Punch Test
The punch test example is illustrated in Figure 21.The model geometry had 561 nodes, and the spline weight function and linear monomial basis were used.Both the left edge and the right edge were constrained from moving in X 1 direction and were traction free in X 2 direction.The bottom edge was constrained from moving in X 2 direction and was traction free in X 1 direction.To simulate the compression from a punch on top, uniform prescribed displacement was applied on the left half of the top edgein the X 2 direction, and the traction free condition was enforcedon the left half of the top edgein the X 1 direction.The right half of the top edge was traction free.For the punch test, elastoplastic material was investigated.A prescribed displacement U 2 of -0.05 m was first applied to the left half of the top edge.Figure 22 shows the

Concluding Remarks
In this paper, the large deformation elastoplasticmeshless integral method based on regularized boundary integral equation [1] is presented.The updated Lagrangian governing integral equation is obtained from the weak form of elastoplasticity based on Green-Naghdi's theory over a local sub-domain.Green-Naghdi's theory starts with the additive decomposition of the Green-Lagrange strain into the elastic part and plastic part and considers a J 2 elastoplastic constitutive law that relates the Green-Lagrange strain to second Piola-Kirchhoff stress.The generalized collocation method is employed to enforce the essential boundary conditions exactly, which is simple and computationally efficient.The natural boundary conditions are incorporated in the system governing equations and require no special handling.Numerical results show that this method is accurate and robust.
gradient, E is the Green strain tensor associated with the displacement u, ρ 0 is the mass density in the reference configuration Ω 0 , f I is the body force per unit mass, e IJKL C is the linear elastic stiffness matrix, and cor IJ dT is the stress correction term because of nonlinearities in either geometry or material.An index I following a comma designates partial differentiation with respect to X I , and repeated indices indicate summation over the dimensionality of the problem.The essential and the natural boundary conditions on the boundary Γ are respectively:

JT
is the J-th component of the traction at the reference state, I u  is the displacement increment from the reference state to the current state, 0 LK T is the second Piola-Kirchhoff stress at the reference state, LK T  is the second Piola-Kirchhoff stress increment from the reference state to the current state, and   , a  X Y is the Dirac delta function.

Figure 2 .
Figure 2. Schematic diagram showing the sub-domain for an interior or a boundary node a Y .

Figure 3 .
Figure 3. Exclusion of a tiny sphere Ω Δ of radius Δ centered at a node for removing the strong singularity.

2 1 
    is the internal boundary angle subtended by material at a Y on the boundary, as shown in Figure 4.

Figure 4 .
Figure 4. Schematic diagram showing the internal boundary angle θ 2 -θ 1 = θ at node a Y on the boundary.Several special cases are worth noting.For an interior node, 2π   ; for a boundary node where the boundary is smooth, π   ; for a corner node, θ = corner angle.

aY
, which is a sphere (3D) or disk (2D) centered on a Y with a radius a w l .A weight function a w is a continuous function that is positive in the support domain and zero outside, i.e.

Y
, located entirely inside 0  , is a sphere or part of a sphere centered on a Y with a radius a s h .

Figure 5
illustrates the meaning of local sub-domain and support domain.Two other frequently used concepts are the domain of definition and the domain of influence.The domain of definition of a point X is the set of all nodes whose weight functions are non-zero at X, while the domain of influence of a node a Y is the set of all nodes whose weight functions are non-zero in some part of the subdomain of node a Y .The domain of definition and the domain of influence are convenient terms in the description of MLSA and local boundary integrals, and are illustrated schematically in Figure 6.

Figure 5 .
Figure 5. Schematic diagram illustrating the meaning of local sub-domain and support domain for node a Y and node b Y .

Figure 6 .
Figure 6.Schematic diagram showing the domain of influence for node a Y and the domain of definition for point X.
or a quadrature point), at least m weight functions are nonzero and the nodes in the domain of definition of X are not arranged in a degenerate pattern (such as on a straight line).In MLSA, the shape function related to node a Y is   a  X .The size of the support domain should be large enough to ensure the coupling between a minimum set of nodes, but small enough to capture local variations.
of support domain.The weight function has only one parameter, the size of support domain a w

1 2 ,
n n is the normal to the plane passing X over which the traction acts.For a node b Y , we define N and b the program converges in this iteration, update the geometry, the displacement field, and the second Piola-Kirchhoff stress for each Gaussian point (Section 7.2) and obtain the corresponding Cauchy stress by Equation (57): compute the fourth term of Equation(47), are predicted by Equations (61) and (62) with     to zero.The parameter EPF (elastoplastic flag) indicates the stress state at a Gaussian point under consideration with EPF = 0 being elastic and EPF = 1 being elastoplastic.Initial values of EPF for all Gaussian points are set to zero.

Figure 8 .Figure 9 .Figure 10 .
Figure 8.(a) The convergence test of meshless method for Cauchy stress σ 11 for elastic case; (b) The convergence test of meshless method for Cauchy stress σ 33 for elastic case; (c) The convergence test of meshless method for displacement of upper edge for elastic case.FEM simulation was carried out for comparison.Good agreement between the meshless results and FEM results is evident.Nonlinear behavior, although not strong, is ) to (c) show variations of σ 11 , σ 33 , and U 2 , respectively, of the upper edge (denoted as b) with the number of load increments, together with the corresponding FEM results.The figures indicate trend of convergence, and converged values for FEM and meshless method are close to each other.The meshlessresults for 12 load increments are σ 11 = 426 MPa, σ 33 = 212 MPa, and b = -0.09009m.The corresponding results of FEM with 12 load increments are σ 11 = 426 MPa, σ 33 = 212 MPa, and b = -0.08978m.

Figure 11 .
Figure 11.(a) The convergence test of FEM results and meshless results for Cauchy stress σ 11 for elastoplastic case; (b) The convergence test of FEM results and meshless results for Cauchy stress σ 33 for elastoplastic case; (c) The convergence test of FEM results and meshless results for displacement of upper edge for elastoplastic case.

Figure 12 .
Figure 12.FEM results versus meshless results of σ 11 for uniaxial tension simulation for elastoplastic case.

Figure 13 .
Figure 13.FEM versus meshless results of Cauchy stress σ 33 for uniaxial tension simulation for elastoplastic case.

L 2 -Figure 15 .
Figure 15.(a) Analytical versus meshless results of Cauchy stress σ 11 of the shear test for elastic case; (b) Analytical versus meshless results of Cauchy stress σ 22 of the shear test for elastic case; (c) Analytical versus meshless results of Cauchy stress σ 12 of the shear test for elastic case.

Figure 16 .
Figure 16.FEM results versus meshless results of the shear test for elastoplastic case.

Figure 17 .
Figure 17.Rotation of a prestressed square with no deformation.(a) Original configuration; (b) Current configuration.densities, virtually identical results were obtained.Fig- ures 18-20 show meshless and analytical results of σ 11 , σ 12 , and σ 22 , respectively, versus different rotation angle θ for the rotation tests.The figures reveal that the meshless results are practically identical with the analytical results.

Figure 18 .Figure 19 .
Figure 18.Analytical versus meshless results of Cauchy stress σ 11 for rotation test for different rotation angle.

Figure 20 .
Figure 20.Analytical versus meshless results of Cauchy stress σ 22 of the rotation test for different rotation angle.

Figure 21 .
Figure 21.The square plate for punch test.

Figure 22 .
Figure 22.The undeformed model (solid diamonds), deformed meshless model (solid squares), and FEM model (triangles) for the punch test with U 2 = -0.05m. undeformedmeshless model (solid diamonds), deformed meshless model (solid squares) and FEM model (empty triangles).The meshless results match FEM results very well and the L 2 -norm between meshless and FEM results is 0.0029.Figures 23 and 24 present the distribution of σ 22 and von Mises stress, respectively, along X 2 = 0, indicatingthat the meshless results and FEM solution are comparable with each other.

Figure 25
Figure25 shows the undeformed model (solid diamonds), deformed meshless model (solid squares), and FEM model (triangles) when the prescribed displacement of the left half of the top edge was increased to U 2 = -0.08 m.The L 2 -norm error between meshless results and FEM solution is 0.0056.Figures26 and 27present the distribution of σ 22 and von Mises stress, respectively, along X 2 = 0.The two figures indicate that the meshless results are in reasonable agreement with FEM results.
If node