A simple damage-gradient enhanced elastoplastic formulation and its numerical implementation

– This paper presents a ductile damage-gradient based nonlocal and fully coupled elastoplastic constitutive equations by adding a Helmholtz equation to regularize the initial and boundary value problem (IBVP) exhibiting some damage induced softening. First, a thermodynamically-consistent formulation of gradient-regularized plasticity fully coupled with isotropic ductile damage and accounting for mixed non linear isotropic and kinematic hardening is presented. For the sake of simplicity, only a simpliﬁed version of this model based on von Mises isotropic yield function and accounting for the single nonlinear isotropic hardening is studied and implemented numerically using an in house FE code. An additional partial diﬀerential equation governing the evolution of the nonlocal isotropic damage is added to the equilibrium equations and the associated weak forms derived to deﬁne the IBVP (initial and boundary value problem). After the time and space discretization, two algebraic equations: one highly nonlinear associated with the equilibrium equation and the second purely linear associated with the damage non locality equation are obtained. Over a typical load increment, the ﬁrst equation is solved iteratively thanks to the Newton-Raphson scheme and the second equation is solved directly to compute the nonlocal damage ¯ D at each node. All the constitutive equations are “strongly” aﬀected by this nonlocal damage variable transferred to each integration point. Some applications show the ability of the proposed approach to obtain a mesh independent solution for a ﬁxed value of the length scale parameter. Comparisons between fully local and nonlocal solutions are given.


Introduction and state of the art
In the framework of conventional (local) plasticity, materials with induced strain softening lead to IBVP's exhibiting a solution highly sensitive to the mesh size.A typical example is given by the damage-induced plastic flow localization inside narrow shear bands prior to the ductile fracture [1,2].Many works have shown that for this kind of problems the solution of the associated IBVP exhibits a strong dependency to the space and time discretization mainly in the softening stage [3][4][5].
In fact, as flow localization takes place and the cracks being initiated in the material, the mechanical fields as well as some material properties become highly inhomogeneous leading to a non-simple material as defined by Noll [3].For this kind of non-simple materials, the knowledge of the transformation gradient (or the first displacement gradient) is not enough to determine the kinematic and state variables of the body.Accordingly, higher gradients of the displacement and/or of some state variables are required to define the mechanical state of each material point of the body.
When the fully coupled constitutive equations are used and because the damage field is obtained numerically by FEA, there is no explicit expression of damage distribution available for the calculation of the damage first and second gradients.A typical way currently followed in dealing with these problems is to consider the first damage gradient as a state variable in the state and/or dissipation potentials, as can be found in [4][5][6][7][8].
Nedjar [5] defined a damage variable under the form of continuity function β(x, t) having the value 1 when the material is undamaged (or fully continuous) and the value 0 when it is completely damaged (fully discontinuous).The damage gradient is introduced to account for the influence of the neighbouring points.The basic idea of his model is to take the damage gradient as new independent internal variable that has additional evolution equations.Consequently, the damage gradient has to be treated as an independent dissipative mechanism together with the damage dissipation.Nedjar [5] presented the following set of equations of balance: where k is a material parameter, which measures the influence of the neighbourhood on the damage of the point under concern and c is a material parameter [5].Equation ( 1) is the conventional equilibrium equation, and Equation ( 2) is the damage evolution equation, which has the same form as the classical heat equation, consequently it can be solved with the same mathematical tools.Results presented by Nedjar [5] indicate that this model has capability to simulate the damage propagation, but its localization phenomenon is not as apparent as that shown in experimental data.
Liebe et al. [9] proposed another gradient-dependent model.Instead of the gradient of plastic hardening variable, the gradient of damage is adopted in Liebe et al. [9] model.The Helmholtz free energy of the process consists of two additive terms according to: Following the standard thermodynamics of irreversible processes Lemaitre [10], and after some specific mathematical transformation the authors Liebe et al. [9] the state relations and the following local dissipation inequality: where Y is the state variable associated with D and Y is the state variable associated with ∇D.Note that the present format of the dissipation inequality suggests independent evolution equations for D and ∇D.Thus due to the compatibility between D and ∇D, these evolution equations would have to be designed such that compatibility is satisfied.To overcome this nontrivial task, Pamin et al. [11] propose introducing a bilinear form of the dissipation power as follows: Comparing the Equation (4) with Equation ( 5) it is obtained: where P is the residual term for a material point of inhomogeneous material and is also regarded as the influence of mechanical variables at neighboring points.
In nonlocal theories, it was defined by Maugin [12] that the integral of P over a representative volume element Ω should be zero i.e.: Substituting the Equation ( 6) in the Equation ( 7) and using the divergence theorem yields to: where n is the vector representing the outward normal to the boundary ∂Ω.From Equation (8) the following local condition results: Geers et al. [6] propose more widespread method, based on the postulation of some appropriate Helmholtz equation for some selected internal variable concerned by the nonlocal effect.For any classical local variable V one associates a nonlocal variable V solution of a partial differential equation of type (9) rewritten here under the form [6,13,14]: where is a kind of internal length scale related to the material under concern [6,7,13,15,16].It has been shown in [6] that, the Equation ( 10) can be obtained as a particular case of the integral form proposed for example in [15,16] to regularize some damage related internal variables.
Voyiadjis et al. [17] introduced a more general and more complex nonlocal formulation.It consists to introducing the first gradients of the overall strain-like state variables into the state potential to have: ρψ = ρψ ε e , ∇ε e , α, ∇α, r, ∇r, D, ∇D from which the state relations are obtained: where σ and X are third-rank tensors.However, the complexity of this formulation limits its application to may be some nonlocal variables.Recently, Forest [18] by using the generalized continuum mechanics theories pioneered since 1900 by the Cosserat brothers [19] and extensively developed from theoretical point of view during the years 1960 [20][21][22][23][24], has shown that all the recent nonlocal or higher gradients formulations derive from an appropriate micromorphic theory.Particularly, in [8] the application of the micromorphic theory, allows to obtain a generalized thermo-elasto-viscoplastic-damage constitutive equations including the damage effect, which extend various nonlocal and gradient enhanced formulations recently purposed.Following this idea, Saanouni et al. [25][26][27] proposed fully coupled micromorphic constitutive equations and their implementation in dynamic explicit FE code showing their ability to ensure the mesh independency of the IBVP's solution.
This study aims at presenting a gradient-dependent nonlocal constitutive model for isotropic damage coupled with plasticity for the simulation of the failure phenomena in the metal forming processes.The evolution equations are obtained from the thermodynamics of irreversible processes.A simple version of this model based on von Mises yield function with non linear isotropic hardening fully coupled with nonlocal damage variable is numerically studied.
Note that, for the sake of simplicity, the present nonlocal formulation is limited to the small elastoplastic strain assumption and does not consider the damage induced volume variation.The extension to more realistic constitutive equations with kinematic hardening, induced volume variation under finite plastic strains and using the framework of generalized continuum theory (micromorphic theory) is developed in [25,26].

Nonlocal formulation 2.1 Nonlocal state potential
In this section a damage gradient enhanced formulation developed at UTT since 2002 is recalled.Considering only the mechanical phenomena (neglecting the thermal phenomena), the Clausius-Duhem inequality can be written under a general integral form [28]: Assuming small elastic strain, the total strain rate is additively decomposed as: The state potential, convex function of the state variables and damage gradient, is postulated: Accordingly ( 12) can be written as: State relations are classically obtained as: The resulting residual dissipation inequality writes: with ∇D = ∇ Ḋ + V div − − → ∇D where V is the velocity field of the damaged zone.Accordingly (18) becomes: Note that if Y or ∇D is zero, the classical local intrinsic dissipation is recovered.In Equation ( 19), the state relation being defined by ( 16) and ( 17) the flux variables εp , ṙ, α, Ḋ, ∇D should be derived form appropriated dissipation potential together with yield function.It is worth noting, that the thermodynamic force associated with Ḋ is Ȳ = Y + div Y and not only Y as in classical local formulation.
Following the total strain equivalence hypothesis [29], we introduce the effective state variables in the case of isotropic damage: Those effective variables are introduced in dissipation and state potentials in order to perform the strong coupling between the isotropic ductile damage and the elastoplastic behaviour with hardening.Assuming the initial anisotropy of the overall physical phenomena excepting the ductile damage, the state potential can be chosen as in [29,30]: where the fourth rank symmetric tensors Λ and C are the elasticity and kinematic hardening modules and Q is the isotropic hardening modulus all relative to the virgin undamaged material.Form this potential the state relations are easily obtained: In Equations ( 23), ( 24) and ( 25) the local damage variable D has been replaced by its nonlocal counterpart D directly in the state relations, in order to define the "nonlocal" properties of the damaged material by Λ =

Nonlocal dissipation analysis
Following the classical standard material approach, the fluxes variables are deduced from the plastic potential F σ, R, X, D together with the yield function f σ, R, X in the framework of non associative plasticity theory, both closed convex functions, in the effective stress space, of their main arguments.These are chosen as in Saanouni et al. [29,30] but replacing the local damage variable D by its nonlocal counterpart D: is the equivalent stress in the Hill sense, σ y is the yield stress in simple tension, a and b are the nonlinearity parameters relative to the kinematic and isotropic hardening respectively and S, s, β, Y 0 are some material parameters characterizing the ductile damage evolution.Finally the positive definite symmetric fourth rank tensor H is the classical Hill operator characterizing the initial orthotropy of the plastic flow thanks to six material constants namely: F, G, H, L, M, N .Assuming the generalized normality rule, the evolution equations relative to each dissipative phenomenon are easily obtained: The plastic multiplier λp is deduced from the consistency condition applied to the yield function in case of the plastic yielding to obtain: The scalar H pd > 0 is the elastoplastic hardening modulus given by: in which the following notation has been used for convenience: It is worth noting that Equation ( 32) is never used to calculate the plastic multiplier, since this latter is computed numerically as can be found in Section 3.3.here after.The nonlocal damage variable D is governed by the following Helmholtz problem (see Eq. ( 10)): where ω = 2 is the characteristic length scale of the material.This is the strong form related to the ductile damage non locality which will be used to derive a weak form to define the associated IBVP.Finally let's mention that, for the sake of simplicity, a fully isotropic version of the fully coupled nonlocal constitutive equation will be used and implemented into an in house FE code.This is done simply, by neglecting the kinematic hardening variables (α = X = 0) and by assuming the von Mises yield function defined by:

Variational formulation
The IBVP under concern is defined by the two strong forms corresponding to the equilibrium and damage non locality recalled here as: Using the weighted residual method, the weak forms associated to (36) are: Taking into account the natural boundary conditions from (36) leads: Using (38) into (37), allows obtaining the weak forms:

Finite element discretization
Now consider that the domain Ω t is discretized into a finite number of appropriate finite elements noted Ω (e) and the time interval is discretized into a finite number of sub-intervals typically noted [t n , t n+1 = t n + Δt] where Δt is a non constant time increment.Inside each FE Ω (e) the classical Bubnov-Galerkin approximation of the main unknowns namely the displacement field and the nonlocal damage, are used: ), one can obtain the following global algebraic system to be solved with respect to the main nodal unknowns {u} and D : where the following notations have been used (Nte is the total number of elements, Ntpg is the total number of Gauss points of each element, π j is the weight of the jth integration point, J e v and J e s are the determinant of the volume and surface Jacobean matrices mapping the each element to its reference one): In this paper, taking the advantage of linearity of Equations (43b), the following global resolution scheme is adopted: first, the Equations (43a) is resolved iteratively by using Newton-Raphson scheme using the old value of the variables Dn .After convergence u n+1 , σ n+1 , R n+1 , Y n+1 , ε p n+1 and D n+1 are obtained.Secondly, the linear Equations (43b) is directly solved to calculate Dn+1 .

Local integration scheme
At each time t n+1 , the calculation of the internal forces {F e int } (Eq.( 44)) and F e D (Eq. ( 47)) needs the computation of the stress {σ} n+1 and local damage {D} n+1 .This is done through the numerical integration of the fully coupled constitutive equations presented above (Eqs.( 29)-( 31)).The standard elastic prediction and plastic correction algorithm [29,32,33] is used for the computation of the stress tensor together with the other state variables of the model.

Elastic prediction
Elastic prediction stage consists to assume that the given total strain increment is elastic without any dissipation (i.e.Δλ p = 0.).This leads to define the "trial" stress as: where ε * n+1 = ε n+1 + Δε − ε p n is the known trial strain supposed as purely elastic.The "trial" yield function is then defined by: The nonlocal damage " D" at each integration point is obtained by linear interpolation from the nodes and is kept constant equal to " Dn " during the iterative resolution of the equilibrium equations.If f * n+1 < 0 the solution is effectively elastic (during the elastic unloading) and the solution is: If f * n+1 > 0, the trial solution should corrected in order to fulfil the yield criterion.

Plastic correction
For the present model the yield function should be fulfilled at each end of the time step i.e. f n+1 (Δλ p ) = 0.The problem is then to compute σ n+1 , R n+1 , ε p n+1 , D n+1 which are plastically admissible i.e. verifying at t n+1 : By using the time discretization scheme of the constitutive equations one can obtain: From (52), the deviatoric stress is obtained: The normals to the yield surface f n+1 and "trial" yield surface f * are given by: ñn+1 = 3 2 Combining Equations ( 54) and (56) leads to: This equation leads to: and Equation ( 57) means that the normal is still constant during the iterative return to the yield surface (radial return algorithm).Using Equation (58b) together with Equation (51) leads to a single scalar equation with Δλ p as a single unknown: This non-linear equation should be solved iteratively thanks to a classical Newton-Raphson scheme in order to compute the value of Δλ p for the current load increment.Once Δλ p is known, the stresses σ n+1 , R n+1 are computed from Equations ( 52) and (53), while ε p n+1 , D n+1 are computed using: whith

Consistent tangent operator
The linearization of Equations (43a) using the iterative Newton-Raphson method, requires the computation of the tangent stiffness matrix with a manner consistent with time discretization of the stress σ n+1 as discussed above.
From Equation (52) the stress at t n+1 can be rewritten: ) where e n+1 is the deviatoric strain tensor.The consistent tangent matrix is then given by: In which all the derivatives can be easily analytically calculated except the terms dΔλ dε n+1 which is obtained from the derivation of fn+1 Equation (58) performed during the local Newton-Raphson iterative procedure applied to fn+1 .The model presented above has been implemented into in house FE program written using FORTRAN language.It has been written using the same format of the program PLAST2 developed by Owen and Hinton in their classical texts book [34].Also, this program has been adapted for structural damage analysis.The finite-element T3 is implemented with three degrees of freedom per node (u 1 , u 2 , D) to resolve the problem.The nodal variable D is systematically transferred to Gauss points to achieve the coupling plasticity-damage.Once a Gauss point is completely damaged ( D = 1), the corresponding element is removed from the loop calculation.The preparation of data as well the visualization of the results, are achieved with the help of the GID (graphical user interface for geometric modelling, data input, and visualization of results for all types of numerical simulation by CIMNE).The Figure 1 shows the flowchart of the program.

Some applications 4.1 Initially homogeneous plane strain tension test
The case of 2D plane strain test is studied with a nonlocal pattern.The material parameters are: E = 210 000 MPa; ν = 0.   (ω = 0, ω = 0.3 et ω = 1.2).The analysis of the iso-values shows that the range and the shape of the localization zone in local computation are always done inside one row of elements.However, in nonlocal case (ω = 0.3 et ω = 1.2), we see that the range and the shape of the localization zones spread over many element rows as the value of ω increases.The volume of the localisation zone depends, therefore, on the internal length ω.When the value of ω increases, we see that the localization range increases too.

Double notched specimen
The second example is done with a double notched plate of plane strain state, with uniform tensile loading at both ends (top and bottom).
The dimensions of the plate are (1 × 2) with a circular cut of R = 0.3.Figure 7 shows the used three meshes.The model will be compared for three given values of ω (ω = 0, ω = 0.2 and ω = 0.7).
The Figure 8 shows the force-displacement curves for the different values of ω.As expected the local model leads to a belated break point compared to the local model (as given by these figures show a decreasing mesh effects as ω increases.Consequently, one can conclude that the (h=.17   proposed nonlocal formulation ensures the uniqueness of the solution of the IBVP until the final fracture.
Figure 9 shows the damage distribution in the section orthogonal to loading direction near the final fracture.One can say, for a coarse mesh h = 0.1 the concentration of damage is located in one Gauss point in the middle of the test tube.On the other hand, for a refined mesh h = 0.05 and h = 0.025 the width of the band location in the middle of the test tube is the same.This can be interpreted with the fact that the iso-values of damage are better distributed in the middle of the test tube.
The Figures 10,11 and 12 show the repartition of the damage and the accumulated plastic strain for three given values of ω (ω = 0, ω = 0.2 and ω = 0.7).One can see that in the local model, the localization will have a width depending on the size of the elements.On the other hand, in the case of nonlocal model, one can observe that the damage and the strain distributions are as better distributed in the middle of the specimen as the value of ω is high and as the mesh is fine.

Conclusion
A simple formulation of nonlocal elastoplastic-damage constitutive equations is proposed and numerically implemented inside an in house FE code in order to investigate the mesh size independency of the solution of the associated IBVP.
The obtained results indicate clearly the efficiency of the proposed formulation to ensure a solution clearly insensitive to the mesh size even when the damage values reach higher values near the final fracture.This result is in contradiction with the one of [28] which shows that the fully damaged zone (macroscopic crack) is always localized inside a single row of elements even if the damaged zone with lower values of damage is clearly mesh independent.Note that the difference between the two formulations resides only on the element type used in the two studies.
This kind of formulation is recently reformulated in the framework of the generalized micromorphic continuum and will be reinvestigated for various cases including the large plastic strains [25,26].

⎩T
{u e } = [N e u ] {u e n } {δu e } = [N e u ] with respect to the space coordinates, the first gradients of these nodal variables are obtained: and (41) into (39) written for one typical element, leads to:⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ I e u (u e , δu e ) = − Ωe [B e u ] {σ e } .dΩe + Ωe {D} dΩ e (42) where the matrix [ω] is the internal length matrix and D is the nodal vector of nonlocal damage.The computation of the integrals in (42) is made using the well known Gauss method.After the classical assembly over the all structural elements (represented by the symbol Nte A e=1 3; σ y = 400 MPa; Q = 1000 MPa; b = 40; β = 1; S = 1; s = 0.8; Y 0 = 0.The boundary conditions are shown in Figure 2: one side is clamped and on the other side a displacement is imposed.Two regular meshes are used with the following sizes: h = 0.17 mm, h = 0.1 mm.The nonlocal computation is done with two different values of ω: ω = 0.3 and ω = 1.2, finally the local computation is achieved with ω = 0.The plots of Figures 3a, b and c show the forcedisplacement curves for the different values of ω.Clearly, when ω is large the solution is independent from the mesh size.The distribution of the local damage D, the plastic deformation and the von-Mises stress are shown in Figures 4, 5 and 6 respectively for the different values of ω

Fig. 2 .
Fig. 2. Mesh and boundary conditions applied to a part on a tensile test (l = 1 mm, L = 3 mm so that L/l = 3).

Fig. 7 .
Fig. 7. Presentation of the double notched specimen and the three used meshes h = 0.1, h = 0.05 and h = 0.025.