An Element Free Galerkin method for an elastoplastic coupled to damage analysis

In this work, a Meshless approach for nonlinear solid mechanics is developed based on the Element Free Galerkin method. Furthermore, Meshless is combined with an elastoplastic model coupled to ductile damage. The efficiency of the proposed methodology is evaluated through various numerical examples. Besides these, twodimensional tensile tests under several boundary conditions were studied and solved by a Dynamic-Explicit resolution scheme. Finally, the results obtained from the numerical simulations are analyzed and critically compared with Finite Element Method results.


Introduction
The Finite Element Method (FEM) is the dominant technique for the simulation of forming processes in past decades. It has made significant contributions in computational mechanics. Nevertheless, this method is not well suited to problems including large deformation of materials and moving discontinuities such as crack propagation requiring considerable meshing and remeshing in structural optimization problems.
Meshless methods have attracted much attention in recent decades. These methods aim at avoiding problems related to the use of a mesh eliminating the well-known drawbacks in the FEM.
While in the Finite Element Method, the shape functions are defined inside the elements, the Meshless methods are based on shape functions constructed the nodes positions.
The goal of this study is to develop a Meshless method known as Element Free Galerkin (EFG) based on Moving Least Square (MLS) approximation for solving structural problems involving large deformations in a more efficient manner than with finite elements.
A code incorporating an implementation of the MLS method has been developed and tested on several simulations.
The first step was to develop and describe the proposed method. Secondly, the used behaviour model was presented. Then, various numerical examples based on elastoplasticity coupled to damage were simulated using Meshless method.
To conclude, the results compare well to finite element simulations.

The Element Free Galerkin Method
The Element Free Galerkin (EFG) is different from the Finite Element method on how the shape functions are constructed. In this work, the shape functions in the EFG method are constructed using the Moving Least Square (MLS) approximation introduced by Lancaster and Salkauskas [1].
In what follows, the use of the symbols [ ], { } and < > means respectively the matrix, the column vector and the row vector.

Moving Least Square approximation
In this section, an MLS approximation, which was introduced by Lancanster and Salkauskaus [1] is described.
In the MLS approximation, the trial function around the Where p i (x), i=1, 2…..m are monomial basis functions, m is the number of terms in the basis and α i (x) are the coefficient of the basis functions In order to precisely obtain the local approximation of the function, the difference between the local approximation {U h (X)} and the function {U(X)} must be minimized by a weighted least-squares method.
Define a functional Where w(X-Xi) is a weight function defined on a domain Ω X called domain of influence and N is the number of particles in Ω X . From we obtain the system of equations of order m to solve: The system can be written in vector form as follows: If we replace α ( X _ ) in the local approximation, we can obtain: Then the expression of the local approximation {U h (X)} can be obtained as: where ɸ(X) is called the shape function, i.e.
The shape function of a particle i in a point X is then:

Boundary conditions
In Meshless methods, imposing the essential (Dirichlet) boundary conditions is not evident as the Finite Element Method because MLS is non-interpolative. A number of techniques for the imposition of essential boundary conditions have been developed, including collocation methods [2,3], Lagrange multiplier method [4], penalty method [4], Nitsche's method [4], coupled meshless-Finite Element Method [6,14] and other methods [7].
The advantages and disadvantages of a variety for methods of imposing essential boundary conditions have been discussed briefly in [5].

Numerical integration
In Meshless methods, the numerical evaluation of integrals is more difficult than in the Finite Element Method (FEM). Indeed, the structure is only known by the description of its boundaries and a distributed set of nodes.
In contrast with the FEM, we do not have the concept of mesh used to build shape functions and which is suited as a basis for realising Gauss integration on each element.
In order to well understand and distinguish the most common integration methods, see [5,14].

Elastoplasticity with damage
To further test our Meshless program, a plasticity model coupled to damage was implemented.
The used behavior model is decomposed into an elastic law and a plastic one. For the elastic law, we consider a linear isotropic behavior (Hooke taking into account damage) which is written as follows: where θ is a flag parameter, D is the variable describing the damage, e  is the elasticity tensor and the symbol : means the double contracted product. f (σ, R, D) is the charge function written as: 1 where σ y is the yield stress, R and σ eq are respectively the hardening function and the equivalent stress. The equivalent stress is written in the case of von-Mises in the following way: where S is the deviatoric part of σ. We note n the derivative of σ eq given as: The hardening function is given as follows: where r is the variable describing the isotropic hardening and Q is a constant of the material. The evolution law of plastic deformation is written as follows: where λ is the plastic multiplier and the point refers to the time derivate. The evolution law of the hardening variable is given by the following equation: where b is a constant of the material. The evolution law of the damage variable is written as follows: with S d , s and β are constants of the material and Y is given by the following equation: The problem is closed by the consistency conditions given by the following relations: 4 Numerical examples and discussion

Equilibrium equation and boundary conditions
The basic equation of dynamic problems is where . indicates the divergence,  is the stress tensor,  is the density and    is the acceleration. With the strain-displacement relationship: In the above, the operator  is the gradient, T is the transpose, the point refers to the time derivate and  is the strain tensor, which is decomposed in elastic and plastic parts as follows: Noting that the large strain formulation is considered and the corotational (or Jaumann) frame well defined and developed in [9] is used.

Explicit scheme
The well-known dynamic explicit algorithm based on the finite difference discretization of the derivatives relative to the time variable is used in the resolution of the resultant system [8].

Implicit integration
An Euler implicit integration algorithm has been established for our behavior model. Therefore and according to [9,10,12], knowing the different variables σ n , R n , ε pn and D n at time increment t n , we propose to determine the parameters σ n+1 , R n+1 , εp n+1 and D n+1 at t n+1 .

MLS discretization
For a two-dimensional plane strain dynamic problem, the vectors of displacement, velocity and acceleration as developed in [11] Substituting these expressions in the weak form in Eq.
(34), the final system is written as follows: In this part, The Galerkin weak form was employed to obtain the discretized equations system. Since the shape functions arising from the MLS approximation do not possess the δ -property i.e., the basic unknown, quantities Ui (xi) contained in the array U in Eq. (9) are not real nodal variables.
Therefore, enforcement of essential boundary conditions needs to be dealt differently. A range of methods has been mentioned in the sub-section 2.2.

Examples and discussion
In this Section, the numerical solution strategy adopted in this work to perform the numerical simulation is summarized. Indeed, the results obtained by implementing numerical simulations of two specimens with the previously described constitutive formulations, will be presented and discussed. Firstly, the description of specimens with different geometries and boundary conditions is undertaken. Then, both qualitative and quantitative comparisons are proven in order to show the efficiency of our developed Meshless method. Finally, representative graphs are chosen to describe the evolution of internal variables. First example considered is a plate with hole under concentrated load at the end (see Figure 1). Second part will investigate the efficiency of the method for a butterfly specimen.
All the aforementioned material parameters are properly listed in table1.
In the case of the plate hole specimen, the geometry, the boundary conditions and the nodes distribution are presented respectively in Figure 1 and Figure 2. In order to capture the evolution of internal variables, a relatively condensate and irregular distribution is used in the region surrounding the critical zone of the specimen (see Figure  2). For the Meshless Method, we have used quadratic weight function with a circular support domain having as radius r=1.4*dx, see [5].
Concerning the integration, we have considered a quadratic cell with four Gauss point in each one.
For the Finite Element simulation, we have used Abaqus with a CPE4R element.
Also, it is very interesting to mention that the same number of nodes has been used for both methods. The evolution of representative variables is analyzed. In particular and for different levels of displacement, the isovalues of the equivalent stress, the cumulative plasticity and the damage obtained by Meshless simulation are compared with those issued from the FEM especially in critical zones, knowing that we have used the same mesh for both methods.

Case 1:
Imposed displacement U f =0.025mm.   As we consider the ductile damage, it follows the cumulative plasticity. Thus, it is equal to zero in the whole plate in this first case.

Case 2:
Imposed displacement U f =0.45mm.   In the three cases, it can be observed an agreement between results got from Meshless and those from FEM.
Certainly, the damage distribution for each method, when the critical damage is attained, can be seen on Figure  13.
In addition, the damage variable field obtained in the numerical simulation and illustrated by the contour plots, shown in Figures 10 and 13, improve clearly the ability of our developed and implemented method to predict damage growth under tensile loading conditions.
The curves depicted in Figures 14, 15, 16 and 17 describe respectively the evolution of the equivalent stress and the damage versus the cumulative plasticity. Then the damage and the cumulative plasticity according the time in different zones of the structure.  As expected, the results obtained in the four points, p1, p2, p3 and p4 differently distributed in the plate (see Figure 1) demonstrates clearly the inhomogeneity of the quantities in the structure.
Certainly, the nearest point (p4) of the hole is the most quickly damaged.
In the next part, numerical results for the butterfly specimen, which is illustrated in Figure 18, are presented. The focus here is to study the behavior of different variables. Figure 18. Geometry and boundary conditions of the butterfly specimen. The specimen was taken from [13] With the aim to more push our study and prove that using our developed Meshless method, we can obtain the fracture location.   When a crack propagates, a new frontier is forming. Consequently, the nodes of the both sides of the crack lips must be taken into account. This can explain the difference between FEM and Meshless method when damage reaches Dc (critical damage).

Conclusions
Based on the MLS approximation and the EFG method, the EFG method for two-dimensional elasto-plasticity problems is presented in this paper. The Galerkin weak form is employed to obtain the equations system, the penalty method is used to apply the essential boundary conditions. Some selected numerical examples are given to show that the EFG method has high computational precision and efficiency.
Further, this paper presents the application of Element Free Galerkin (EFG) method to evaluate the damage evolution. Evaluation is mainly on the efficiency of this method to capture the initiation of damage and its propagation for many steps of time.
By using the MLS approximation with consideration of the elastoplastic model coupled to damage, evolution of damage is obtained.
Validation of the method has been achieved through comparison of the simulated results with those of the Finite Element Method.
According to the comparison, both qualitative and qualitative agreement between the two methods can be found, indicating that the used Meshless method can reasonably predict the propagation of damage until the saturation.