APPLICATIONS OF THE FINITE POINTS METHOD IN THE STRAIN-GRADIENT PLASTICITY

In this paper, the meshless Finite Points Method (FPM) is used for numerical simulation of one and two-dimensional elastoplastic problems, which develop a softening response after to exceed the elastic range. In such problems, a pathological behaviour of the solutions is induced by the ill-conditioning of the partial differential equations after reaching the yield limit, producing the localization of deformations. Specifically, this work presents the implementation of a regularization technique through the enrichment of the constitutive equations using strain gradients, with a characteristic length based on a non-local plastic strain formulation in order to maintain the ellipticity of the differential equations. The strain localization phenomenon is replicated objectively in the computational simulation considering an isotropic Von Mises (J2) softening model. This localization phenomenon is induced weakening a region of the material and also developing problems with geometrical discontinuities. The FPM performs the domain’s discretization in a finite number of points in each of which the punctual collocation of the equations is performed. For these reasons the FPM applies the strong formulation, allowing the use of high-order differentiability shape functions (C2 class or higher). This features can improve the computational cost because of the same shape functions are used in order to approximate non-local strain and the displacement nodal fields. Both fields are obtained by using the iterative Newton-Raphson method. To ensure convergence of the iterative method, the algorithmic tangent operator is obtained by Perturbation Method. The theory is developed for one-dimensional geometries, being extended to two and three-dimensional problems. In order to validate the proposal, a benchmarking is developed with typical problems extracted from literature.


INTRODUCTION
The solutions of the computationally simulations of the elastoplastic problems with softening behaviour, are highly dependent of the discretization used ( [1,2]), concentrating the localization of deformations on thin zones of typical width of the characteristic element of the geometry when the yield limit has been exceeded.
The fundamental reason for this behaviour is attributed to the ill-posed of the differential equations governing the problem of softening ( [1, 3, and 4]).
There are many jobs where the localization of deformation has been treated, many of them using the approach of the finite element method(FEM)( [32]).To treat this behaviour in problems of mechanical damage or softening, traditionally techniques of fracture energy regularization ( [11,23,27]) or enrichment approach of the constitutive equations with strain gradients ([3, 6, 7, 22, 31]) has used.
For to enrich the governing equations of the problem, the non-local plastic strains approach will be used ( [4]) through of an internal hardening variable.These non-local strain gradients add a characteristic length which acts as regularizer of the localization, independent of the discretization employed.
The Finite Point Method (FPM) ( [15,19,24,25]) requires shape functions of class C2 or higher to solve the equation of balance, this feature being used also for non-local plastic strain.
In general, incremental iterative schemes such as Newton-Raphson are used.Thereby, ensure the convergence order of the scheme is requiered.A scheme of constant stiffness matrix is used for compression problems, while for the tension problems a tangent algorithm scheme is obtained by perturbation techniques.
This paper is organized as follows: In Section 2 the basic aspects of the formulation of the MPF is presented.Section 3 introduces the nonlocal formulation of a general variable and gradients approach from it.Section 4 describes the constitutive and consistency relations of elastoplastic problems with isotropic softening and its formulation in the MPF.Section 5 describes the perturbation method for the Algorithmic Tangent Operator.Finally, in Section 6 examples of the literature are developed, which validate the proposal in one and two-dimensional geometries.

FINITE POINTS METHOD
Let Ω I be the interpolation sub-domain called cloud, and let s j with j = 1, 2, . . ., n be a collection of n points with coordinates x j ∈ Ω I .The unknown function u(x) may be approximated within Ω I by where α = [α 1 , α 2 , . . ., α m ] T and the vector p(x), called base interpolating functions, contains typically monomial in the space coordinates ensuring that basis is complete.For a 2D problem, and p = [1, x, y, xy, x 2 , y 2 ] T for m = 6 (3) can be used.
The function u(x) can now be sampled at the n points belonging to Ω I giving where u h j = u(xj) are the unknown values of u at x j point, u j = u(x j ) are the approximate values, and p j = p(x j ).In the case of FEM, the number of points is chosen so that m = n.In this case C is a square matrix [32].If n > m ,C is nor a square matrix and the approximation cannot fit all the u h j values.This problem can be solved determining the u values by minimizing with respect to the α parameters the sum of the square of the error at each point weighted with a function ω(x) as: This approximation is called weighted least square (WLS) interpolation.Note that for ω(x) = 1 the standard least square (LSQ) method is reproduced.Function ω(x) is usually built in such way that it takes a unit vale in the vicinity of the point I typically called star node where the function (or its derivatives) are to be computed and vanishes outside a region Ω I surrounding the point.The region Ω I can be used to define the number of sampling points n in the interpolation region.In this work, he normalized Gaussian weight function ω(x) is used, thus n ≥ m is always required.
Several possibilities for selecting the weighting function ω(x) can be found in [13, 14, 18 and 30].Minimization of eq. ( 5) with respect to α gives where The final approximation is obtained by substituting α from eq. ( 6) into eq.( 1) giving where are the shape functions.
It must be noted that because of least square approximation u(x j ) u(x j ) = u h j .That is, the local values of the approximating function do not fit the nodal unknown values.The WLS approximation described above depends strongly on the shape and the way in which the weighting function is applied.The simplest way is to define a fixed function ω(x) for each of the Ω I interpolation sub-domains [21].
Let ω(x) be a weighting function satisfying The Gaussian weight function used in this work is given by where The support of this function is isotropic, circular and spherical in two and three-spatial, respectively.A detailed description of the effects of the parameters q and β on the numerical approximation and some guidelines for their setting have been presented in [18] and [21].The approximate function û(x) is defined in each interpolation domain Ω I .In fact, different interpolation subdomains can yield different shape function Φ I j .As a consequence of this, a point belonging to two or more overlapping clouds has different values of the shape functions which means that The interpolation is now multivalued within Ω I and, therefore for the approximation to be useful, a decision must be made in order to limit the choice to a single value.Indeed, the approximate function û(x) will be typically used to provide the value of the unknown function u(x) and its derivatives only in specific regions within each interpolation sub-domain.For instance by using point collocation the validity of the interpolation is limited to a single point x I .

NON-LOCAL PLASTIC STRAIN GRADIENTS
The classical theory isn't sufficient for known the local stress state in the particular node by the heterogeneous structure of the material.For this reason a non-local approach is used.([2, 8, 29]) A non-local approach for an arbitrary function is defined as the weighted average about a specifically point on the characteristic domain.This approach looks as follow: The weight function α must introduce a characteristic length has relation with the material microstructural adding information in the neighbour of material ( [2]).In this work, the nonlocal approach is applicated to the plastic invariant κ with the weight function: resulting: Using the described in [17] by the local and non-local relation, here was used: For close the problem is necessary to add the boundary condition ε p (x) • n = 0, where n is the exterior normal at the geometry, this condition is related with thermodynamics subjects( [26,28]).
This work is developed about the isotropic theory of Von Mises plasticity and the flux rule of Prandtl-Reuss ( [5]) defined by: s s s ||s s s|| (17) where n n n is the plastic-flow vector, f the yield function, σ σ σ the stress tensor and s s s is the deviatoric component of σ σ σ.The evolution of the plastic strain and of the invariant plastic are relationated ( [10]) by: So, replacing ( 18) into (19) and considering ( 16) results the non-local form:

NUMERICAL IMPLEMENTATION
The general elastoplastic problems with the non-local plastic gradients approach is gobernated by the relations ( [5,10,16]): where σ σ σ, ε, λ and u are the time-evolution of stress, strain, plastic multiplicator and displacements respectively and the consistence condition: wish at the plastic state: Applicating the gradient approach, as used in [3, 4, 6 and 22] and differentiating the yield function gradient dependent (27) we obtain: Is defined: where h is the modulus of linear hardening.
The punctual collocation of the MPF allows a strong formulation ( [19, 24, and 25]), and according to what was done in [22,23] and using the Einstein notation for each star node with E h the nodal value of E, (22) in ( 21) and ( 28) look as: To give an elastoplastic incremental form at the problem, a displacements control is used on the boundary conditions.Applying the MPF to the boundary conditions is obtained:

ALGORITHMIC TANGENT OPERATOR BY PERTURBATION METHOD
There are problems where obtaining analytical tangent operator is very complicated or even impossible, so the numerical derivation methods are the only option.In this work, a numerical derivation using a perturbation method is performed ( [12,20]).The perturbation method consists of applying a small perturbation to the each component of the strain vector and, using the constitutive equation of the material, determines the variation that will be obtained in the stress tensor due to this perturbation.The j column of the tangent constitutive tensor can be computed as: where dε j is the perturbated strain vector at the j−component, dσ j = σ − σ j , σ is the stress updated and σ j is the stress obtained by the application of the constitutive equation on the dε j .For each dε j , the only component not equal to 0 is the j − th, with a small value typically dε j = ∝ √ EPS where EPS is the EPS-machine.

One-dimensional example
The example consists of a rod subjected to uniaxial tensile load by displacement control [3] as show Fig. 1.The central tenth of the bar is weakened 10% by a reduction in Young's modulus and yield limit to force the localization .The dimensionless geometrical and mechanical parameters are summarized in Table 1 Figure 2 shows the stress-displacement profile for different reasons l / L between the characteristic longitude and total longitude.To validate the method used, the l/L rate used hereinafter in this one-dimensional example will be 0.05.Is verified that from the discretization of 161 points, the solution is practically independent( [11]) of the discretization used, this can be checked by the stress-displacement behaviour(Fig.3) and for the solution obtained for the nonlocal plastic strain (Fig. 4) and total strain (Fig. 5).Simple 2D tension test.This example is based on [2].This corresponds to a flat plate subjected to tension condition under plane stress approach with a degraded area in the center, where the yield limit is reduced to 90% of the value used in the full geometry.The parameters used are presented in the Table 2.
Five different discretizations are tested: 8x10, 15x19, 22x28, 29x37 and 36x46 nodes.Figure 6 shows the dimensions and the deformed geometry under the displacement boundary conditions assigned for the 36x46 nodes discretization.
Figure 7 shows force-displacement behaviour of the simple 2D tension test for the five discretizations used.Tension test with diagonal degradation.This example is based on [9].This corresponds to a flat plate subjected to tension condition under plane stress approach with a diagonal degraded area, where the yield limit is reduced to 80% of the value used in the full geometry.
The parameters used are presented in the Table 3.Two different discretizations are tested: 17x17 and 33x33 nodes.Figure 8 shows the dimensions and the deformed geometry under the displacement boundary conditions assigned for the 33x33 nodes discretization.
Figure 9 shows force-displacement behaviour of the simple 2D tension test with diagonal degradation for the two discretizations used.
2D compression test.This example is based on [27].This corresponds to a flat plate subjected to compression condition under plane strain approach with a degraded area in the lower  left corner, where the yield limit is reduced to 90% of the value used in the full geometry.The parameters used are presented in the Table 4. Six different discretizations are tested: 25x13, 37x19, 49x25, 61x31, 73x37 and 85x43 nodes.Figure 10 shows the dimensions and the deformed geometry under the displacement boundary conditions assigned for the 49x25 nodes discretization.2D compression test with central degradation.To test the insensitivity of the results due at the size of imperfection, a 2D compression test with central degradation is tested.This example is based on [27], which corresponds to a flat plate subjected to compression condition under plane strain approach with a degraded area in the center, where the yield limit is reduced to 90% of the value used in the full geometry.The parameters used are presented in the Table 5.Two different size of imperfections are tested for the 49x25 nodes discretization.Figure 12 shows the dimensions and the deformed geometry under the displacement boundary conditions assigned for the both cases.It also shows that the solution delivered by the proposal presented in this paper is relatively independent of the size of the area artificially degraded.
The value = 1.4901 • 10 −8 for the perturbation method in solving the traction problems is used.Indeed the number of iterations was lower overall, but the time required was considerably higher.

CONCLUTIONS
In this paper the implementation of non-local plastic strain gradients has developed to treat the isotropic softening phenomenon using the meshless Finite Point Method.The regularization technique used produces objective and consistent results for the physical problem in one-dimensional and two-dimensional examples.The numerical results prove the independence of the solution respect to the discretization used and validate the proposal.The perturbation parameter of the method to obtain the algorithmic tangent operator must be studied best way to apply to all kinds of problems.Present and future investigations of the authors are related to the application of the FPM on three-dimensional problems, to improve the perturbation method and complex geometries with irregular discretizations.

Figure 2 .
Figure 2. Stress v/s displacement for different l/L ratios.

Figure 6 .
Figure 6.Dimensions and deformation of the tension problem with central degradation.

Figure 7 .
Figure 7. Force-displacement behaviour of the simple 2-D tension test.

Figure 8 .
Figure 8. Dimensions and deformation of the tension problem with diagonal degradation.

Figure 9 .
Figure 9. Force-displacement behaviour of the 2-D tension test with diagonal degradation.

Figure 10 .
Figure 10.Dimensions and deformation of the compression problem.

Figure 11
Figure 11 shows force-displacement behaviour of the 2D compression test for the six discretizations used.

Figure 11 .
Figure 11.Force-displacement behaviour of the 2-D compression test.

Figure 12 .
Figure 12.Geometries, dimensions and deformations with small(up) and big(down) degradations in two-dimensional compression problem.

Figure 13
Figure13shows force-displacement behaviour of the 2D compression test with central degradation for the small and big central degradation.

Figure 13 .
Figure 13.Force-displacement behaviour of the 2-D compression test with central degradation.

Table 2 .
Dimensionless parameters for simple 2-D tension test.

Table 3 .
Dimensionless parameters for 2-D tension test with diagonal degradation.

Table 5 .
Dimensionless parameters for 2-D compression test with central degradation.