1 Introduction
Constitutive models involving softening all predict unlimited localization of strain and damage. This feature generates such undesired phenomena as absence of energy dissipation during crack propagation and mesh size sensitivity in finite element computations. Gurson Gurson (1977)’s famous model for porous ductile materials, which was derived from approximate limitanalysis of some elementary voided cell in a plastic solid, is no exception. In this model, unlimited localization arises from the softening because of the gradual increase of the porosity.
Several proposals have been made to solve this problem. One of these, due to Leblond et al. Leblond et al. (1994) but based on a previous suggestion made by Pijaudier et al. PijaudierCabot and Bazant (1987) in damage of concrete, comprises adopting a nonlocal evolution equation for the porosity involving some spatial convolution of some “local porosity rate” within an otherwise unmodified Gurson model. This simple proposal has attracted the attention of several authors (Tvergaard and Needleman Tvergaard and Needleman (1995), Tvergaard and Needleman Tvergaard and Needleman (1997), Enakoutsa et al. Enakoutsa (2007); Enakoutsa et al. (2007)). It was notably checked by Tvergaard and Needleman Tvergaard and Needleman (1995) that it allows to eliminate mesh size effects. Also, Enakoutsa et al. Enakoutsa (2007); Enakoutsa et al. (2007) showed that with a minor modification, it leads to great numerical reproduction of the results of typical experiments of ductile rupture.
One shortcoming of Leblond et al. Leblond et al. (1994)
’s proposal, however, is that it is purely heuristic and lacks any serious theoretical justification. This was the motivation for a later, more elaborate and physically based proposal of Gologanu
et al. Gologanu et al. (1997). These authors derived an improved variant of Gurson’s model (the GLPD model^{1}^{1}1GLPD: GologanuLeblondPerrinDevaux.) through some refinement of this author’s original homogenization procedure based on Mandel Mandel (1964)’s and Hill Hill (1967)’s classical conditions of homogeneous boundary strain rate. In the approach of Gologanu et al. Gologanu et al. (1997), the boundary velocity is assumed to be a quadratic, rather than linear, function of the coordinates. The physical idea is to account in this way for the possibility of quick variations of the macroscopic strain rate, such as encountered during strain localization, over short distances of the order of the size of the elementary cell considered. The output of the homogenization procedure is a model of “micromorphic” nature involving the second gradient of the macroscopic velocity and generalized macroscopic stresses of “moment” type (homogeneous to the product of a stress and a distance), together with some “microstructural distance” connected to the mean spacing between neighboring voids.The numerical implementation of the GLPD model into a finite element code is quite involved as this task required to introduce extra degrees of freedom representing strains; these extra degrees of freedom will permit the calculation of the spacial derivative of the strains, but their number will increase from 2 to 6 in 2dimensional calculations and from 6 to 9 in the 3dimensional calculations. This will increase the CPU time for the simulations. Another difficulty lies in the necessary operation of “projection” onto the sophisticated yield locus. An implicit algorithm similar in principle to that classically used for the von Mises criterion, although much more complex in detail, is adopted for this purpose. Convergence of the global elastoplastic iterations was difficult. To preserve the quadractic convergence rate of the global Newton iterations stifness tangent moduli are needed. The derivation of tangent moduli is a very difficult task, especially for constitutive models with complex forms. In this paper we derive an exact consistent stiffness matrix for the GLPD model for porous materals. The derivation is based on the small strain formulation, this choice is consistent with the one adopted in many finite element codes, including Abaqus
, LSDyna, Adina, and Systus which have demonstrated their efficiency. In Section 2 we summarize the consitutive equations of the GLPD model. Section 3 presents some aspects of the numerical implementation of the GLPD model into finite element codes. In Section 4 we provide the expressions for the derivatives of the Cauchy stress tensors and the generalized moments stress tensor the model involves with respect to the main field variables. We do not compute all the terms of the tangent stiffness matrix, but only the ones that are critical for the numerical implementation. Finally, Section 5 assesses the performance of the tangent stiffness matrix toward its capacity to reach quadratic convergence in the simulations of small scale ductile fracture problems.2 The GLPD Model
The purpose of this section is to provide a complete description of the GLPD model developed by Gologanu, Leblond, Perrin, and Devaux. The original reference Gologanu et al. (1997) for the GLPD model is not easily accessible, a summary of the equations of this model is given here. A short presentation of the derivation of these equations derived from some homogenization procedure is also provided below; strictly speaking, this presentation is not indispensable, but it is useful to fully grasp the physical foundations of the GLPD model.
2.1 Generalities
In the GLPD model, internal forces are represented through some ordinary secondrank symmetric Cauchy stress tensor plus some additional thirdrank “moment tensor” symmetric in its first two indices only^{2}^{2}2The component is noted in Gologanu et al. (1997)’s original paper. The present notation leads to more naturallooking expressions.. The components of are related through the three conditions
(1) 
(These conditions may be compared to the condition of plane stress in the theory of thin plates or shells).
The virtual power of internal forces is given by the expression
(2) 
where denotes the domain considered, (: material velocity) the Eulerian strain rate, its gradient, the double inner product and the triple inner product .
The virtual power of external forces is given by
(3) 
where represents some surface traction^{3}^{3}3The general equilibrium equations and boundary conditions corresponding to the expressions (2) and (3) of the virtual powers of internal and external forces need not be given since they are not necessary for the numerical implementation..
The corresponding equilibrium equations read, in the absence of body forces and moments:
(4) 
The boundary conditions are complex and will not be given here. (In fact the numerical implementation of the model will require neither the equilibrium equations nor the boundary conditions, but the sole expression of the virtual power of internal forces).
The hypothesis of additivity of elastic and plastic strain rates reads
(5) 
The elastic and plastic parts , of the gradient of the strain rate here do not coincide in general with the gradients , of the elastic and plastic parts of the strain rate.
2.2 Hypoelasticity law
The elastic parts of the strain rate and its gradient are related to the rates of the stress and moment tensors through the following hypoelasticity law:
(6) 
In these expressions and denote the Lamé coefficients and the mean halfspacing between neighboring voids. (In the homogenization procedure, is the radius of the spherical elementary cell considered). Also, and are the Jaumann (objective) timederivatives of and , given by
(7) 
where is the antisymmetric part of the velocity gradient. Finally
is a vector the value of which is fixed by equations
(1) (written in rate form, ):(8) 
(This vector may be compared to the throughthethickness component of the elastic strain rate in the theory of thin plates or shells, the value of which is fixed by the condition of plane stress).
2.3 Yield criterion
The plastic behavior is governed by the following Gursonlike criterion:
(9) 
In this expression:

(: deviator of ) is the von Mises equivalent stress.

is the mean stress.

represents a kind of average value of the yield stress in the heterogeneous metallic matrix, the evolution equation of which is given below.

is a parameter connected to the porosity (void volume fraction) through the relation:
(10) where is Tvergaard’s parameter, the critical porosity at the onset of coalescence of voids, and () a factor describing the accelerated degradation of the material during coalescence Tvergaard (1981); Tvergaard and Needleman (1984),

is a quadratic form of the components of the moment tensor given by
(11) where and are the quadratic invariants of defined by:
(12) and denoting the mean and deviatoric parts of , taken over its first two indices.

Again, is the mean halfspacing between neighboring voids.
2.4 Flow rule
The plastic parts of the strain rate and its gradient are given by the flow rule associated to the criterion (9) through normality:
(13) 
were as
The term in equation (13) represents a rigidbody motion of the elementary cell, which is left unspecified by the flow rule but fixed in practice by conditions (1). (The vector may be compared to the throughthethickness component of the plastic strain rate in the theory of thin plates or shells, the value of which is fixed by the condition of plane stress).
The values of the derivatives of the yield function in equations (13) are readily calculated to be
(14) 
2.5 Evolution of internal parameters
The evolution of the porosity is governed by the classical equation resulting from approximate incompressibility of the metallic matrix:
(15) 
The parameter is given by
(16) 
where is the function which provides the yield stress of the matrix material in terms of the local equivalent cumulated plastic strain , and represents some average value of this equivalent strain in the heterogeneous matrix. The evolution of is governed by the following equation:
(17) 
3 Numerical implementation
The GLPD model described in Section 2 has been incorporated into the Systus FE code developed by ESI Group, in the 2D case. The trickiest features of the numerical implementation, which stands as an extension of those proposed by Aravas Aravas (1987) and Enakoutsa et al. Enakoutsa et al. (2007) for the original Gurson model, are presented here. Emphasis is mainly placed on the complex problem of projection of the (supposedly known) elastic stress predictor onto the yield locus defined by the yield function (9). (This problem will be called the projection problem for shortness in the sequel).
3.0.1 The GLPD model and the class of generalized standard materials
The class of generalized standard materials, as defined by Halphen and Nguyen Halphen and Nguyen (1975), consists of elasticplastic materials for which the plastic strain plus the internal parameters collectively obey some “extended normality rule”. This class is remarkable in that as shown by Nguyen Nguyen (1977), for such materials, provided that the flow rule is discretized in time with an implicit (backward Euler) scheme, the projection problem is equivalent to minimizing some strictly convex function, which warrants existence and uniqueness of its solution.
It so happens that the GLPD model fits into the framework of generalized standard materials for a fixed porosity. This property is tied to the special evolution equation (17) obeyed by the hardening parameter . The proof is provided in Enakoutsa’s Enakoutsa (2007)’s thesis and is in fact a straightforward extension of that given by Enakoutsa et al. Enakoutsa et al. (2007) for the original Gurson model.
This property strongly suggests adopting an implicit algorithm to solve the projection problem, to take advantage of the guaranteed existence and uniqueness of the solution. However, since the porosity must not be allowed to vary for the GLPD model to be “generalized standard”, it appears necessary, to benefit from this property, to use an explicit scheme regarding this specific parameter. Then will be fixed during the whole calculation of the values of field quantities at time from their values at time , and updated (using a discretized version of equation (15)) only at the end upon the convergence; the projection algorithm will then be exactly the same as if the porosity were a constant.
We shall therefore use the explicit estimate of the porosity at time
given by(18) 
and the explicit estimate of the parameter resulting from there, during the whole “transition from time to time ”, but the projection algorithm developed will otherwise be fully implicit with respect to all other parameters, that is the components of the plastic strain and the plastic strain gradient and the hardening parameter . From now on, all quantities will conventionally be denoted with a lower index if considered at time , and without any special symbol if considered at time . From now on, all quantities will implicitly be considered at time .
3.0.2 Parametrization of the yield locus
One key point of the procedure of solution of the projection problem, aimed at reducing the number of unknowns, lies in a suitable partial parametrization of the yield locus defined by the yield function (9). This parametrization is inspired from the classical one for an ellipse and obtained by looking for the maximum possible value of the quantity , namely , and then writing this quantity in the form for some angle and solving the equation with respect to . One thus gets
(19) 
The sign of the parameter is introduced into equation (19) in order to allow for negative as well as positive values of .
3.1 Solution of the projection problem for a fixed hardening parameter
Momentarily assuming the value of the current yield stress to be known, we shall now show how the projection problem can be solved through combination of the yield criterion and the flow rule. This problem will be reduced to a system of two coupled equations on the unknowns and , which are solved numerically to get
(20) 
and
(21) 
where , , , and are defined as in Enakoutsa and Leblond (2009) .
This is the second equation of the system looked for on the unknowns and . The lefthand side depends only on and the righthand side only on .
The simplest way to solve the system of equations (20, 21) on and may comprise using equation (20) to express as a function of , and inserting its expression into equation (21) to get an equation on the single unknown , to be solved by Newton’s method. But numerical experience reveal that the convergence of the Newton iterations is then often problematic. An alternative method comprises solving equation (20) on through Newton iterations, being calculated as a function of at each step by solving equation (21) through Newton subiterations.
3.1.1 Iterations on the hardening parameter
The value of the current yield stress has been assumed to be known up to now. In reality, it is not and must be determined iteratively. This is done using a fixed point algorithm, starting from the value at time , solving the projection problem with this value, updating it using a discretized form of the evolution equation (17), resolving the projection problem with the new value, etc. up to convergence.
3.1.2 Other features of the numerical implementation
The numerical implementation of the GLPD model, just like that of all secondgradient models, raises a difficulty tied to the clear necessary use of the second derivatives of the shape functions. This seems to require elements of class which are never available in standard finite element codes. This difficulty is circumvented through some trick suggested by Gologanu et al. Gologanu et al. (1997) themselves and used since in several works (see, for instance, Shu et al. Shu et al. (1999), Forest et al. Forest et al. (2000), Matsushima et al. Matsushima et al. (2000) ) for the numerical implementation of various secondgradient models. This trick comprises introducing a new nodal variable in the form of a symmetric secondrank tensor , replacing the gradient of the strain rate by the gradient of in all equations and imposing the approximate coincidence of and at the Gauss points through some penalty method.
The advantage is that the components of are then got from the nodal values of the new variable and the sole first derivatives of the shape functions; thus, classical elements of class are sufficient. The price to pay is an increased number of nodal degrees of freedom: six () instead of two () in 2D. Also, imposing the internal constraints by a penalty method may give rise to locking phenomena, for which sub integration is there natural remedy. In practice, 8node quadratic elements are used with 4Gauss points integration. Numerical experience reveals that this suffices to prevent locking.
4 Derivation of the tangent stiffness matrix of the GLPD model
In this section we derive the equations of the tangent stiffness matrix to circumvent global Newton iterations convergence problems the numerical simulations with the GLPD model have revealed. We do not compute all the terms of this matrix but only the most important ones. Thus, we do not take into account the variations of the stresses due to the variations of the temperature; this is strictly licit because the correction of the stresses is perfomed with explicit scheme, using the stresses at time and not at time ; and as a conseqence, the correction is independent upon the displacement increment between these instants.
Also, we do not take into account the variations of the stresses due to the objective derivation in the law of hypoelasticity, which does indeed depend on and therefore generates theoretically a contribution in the stiffness matrix. Similarly, the influence of geometry on the residual forces will not be taken into account. We can summarize all this by saying that the calculation of the tangent stifness matrix will be carried out by neglecting the effects of large deformations. This choice is in conformity with the standard one used in many finite element codes for the calculation of the tangent stiffness matrix for usual elastoplasticity models (without damage) which practical numerical simulations involing these models has demonstrated the robustness.
We assume that
Hence, we can rewrite the equations giving the expresions , and . With these notations, we have:
and
The value of depends on through its arguments , , and . Moreover, the equations to be solved on and are expressed as
In the subsequent, we introduce the following expressions:
where the terms G and F depend on , , , and , , .
4.1 Calculation of the derivatives
First, for practical purposes, we find the following intermediate expressions:
4.1.1 Derivatives of the
The terms depends on , et ; as a consequence, we have
4.1.2 The derivatives of the terms and
The terms and also depend on ; thus we get
and
4.1.3 The derivatives of and
The terms and depend on the variables . Taking the derivatives we get