A Free energy expressed in terms of the logarithmic elastic strain \boldsymbol{\epsilon}_{e}

Finite-strain formulation and FE implementation of a constitutive model for powder compaction


A finite-strain formulation is developed, implemented and tested for a constitutive model capable of describing the transition from granular to fully dense state during cold forming of ceramic powder. This constitutive model (as well as many others employed for geomaterials) embodies a number of features, such as pressure-sensitive yielding, complex hardening rules and elastoplastic coupling, posing considerable problems in a finite-strain formulation and numerical implementation. A number of strategies are proposed to overcome the related problems, in particular, a neo-Hookean type of modification to the elastic potential and the adoption of the second Piola-Kirchhoff stress referred to the intermediate configuration to describe yielding. An incremental scheme compatible with the formulation for elastoplastic coupling at finite strain is also developed, and the corresponding constitutive update problem is solved by applying a return mapping algorithm.

Keywords: plasticity; elastoplastic coupling; finite element method; automatic differentiation

1 Introduction

The formulation and implementation of elastoplastic constitutive equations for metals at large strain have been thoroughly analyzed in the last thirty years, see for instance [1, 2], so that nowadays they follow accepted strategies. For these materials, pressure-insensitive yielding, -independence, and incompressibility of plastic flow strongly simplify the mechanical behaviour, while frictional-cohesive and rock-like materials (such as granular media, soils, concretes, rocks, ceramics and powders) are characterized by pressure-sensitive, -dependent yielding, dilatant/contractant flow, nonlinear elastic behaviour even at small strain and elastoplastic coupling. There have been several attempts to generalize treatment of metal plasticity at large strain in this context [3, 4, 5, 6, 7, 8, 9], but many problems still remain not completely solved. These include the form of the elastic potential, the stress measure to be employed in the yield function, which has to provide an easy interpretation of experiments, the flow rule and the elastic-plastic coupling laws.

The main difficulty in the practical application of finite-strain elastoplasticity models, as opposed to their small-strain counterparts, is related to development and implementation of incremental (i.e., finite-step) constitutive relationships. The difficulties lie, for instance, in formulation and solution of the highly nonlinear constitutive update problem, consistent treatment of plastic incompressibility (or plastic volume changes), and consistent linearization of the incremental relationships. The last issue is of the utmost importance for overall computational efficiency of the finite element models because consistent linearization (consistent tangent) is needed to achieve the quadratic convergence of the Newton method.

In the present paper, the model for cold forming of ceramic powders proposed by Piccolroaz et al. [10, 11] (called ‘PBG model’ in the following) is developed for large strain analyses, implemented in the finite element method and numerically tested. The need for this large-strain generalization is related to the fact that during ceramic forming the mean strain can easily reach 50%, while peaks can touch 80%. The differences between a small strain and a large strain analysis can be appreciated from Fig. 1, where small-strain and large-strain predictions are reported for the force/displacement relation at the top of a rigid mould containing an alumina ceramic powder. Results (taken from [12] and pertaining to a flat punch and to a punch with a ‘cross-shaped’ groove, respectively reported in Fig. 1a and 1b) clearly show that the large-strain analyses are more consistent and in closer agreement with experimental results than the analyses performed under the small strain hypothesis.

      (a)       (b)

Figure 1: Predictions of the small- and finite-strain versions of the PBG model for compression of alumina powder by: (a) flat punch (uniform uniaxial strain), (b) cross-shaped punch (nonuniform 3D deformation). The results are taken from [12].

The model for powder compaction can be considered as paradigmatic of the difficulties that can be encountered in the implementation of models for geomaterials, since many ‘unconventional’ features of plasticity are simultaneously present to describe the complex transition from a loose granular material (the powder) to a fully dense ceramic (the green body). These difficulties enclose: (i.) the pressure-sensitive, -dependent yield function introduced by Bigoni and Piccolroaz [13] (‘BP yield function’ in the following), which is defined in some regions outside the elastic domain; (ii.) a nonlinear elastic behaviour even at small strain, (iii.) changes in elastic response coupled to plastic deformation (elastoplastic coupling).

In this work, incremental (finite-step) constitutive equations are developed and implemented for the finite-deformation version of the PBG model. In order to improve the computational efficiency, the original model [11] is slightly modified, but its essential features, including the elastoplastic coupling, are preserved. Note that a consistent finite-element implementation of the elastoplastic coupling at finite strain has not been reported in the literature so far. The model is applied to simulate ceramic powder compaction with account for frictional contact interaction.

The above-mentioned implementation difficulties are efficiently handled by using an advanced hybrid symbolic-numeric approach implemented in AceGen, a symbolic code generation system [14, 15]. AceGen combines symbolic and algebraic capabilities of Mathematica, automatic differentiation (AD) technique, simultaneous optimization of expressions and automatic generation of computer codes, and it is an efficient tool for rapid prototyping of numerical procedures as well as for generation of highly optimized compiled codes (such as finite element subroutines). Finite element computations have been carried out using AceFEM, a highly flexible finite element code that is closely integrated with AceGen.

Selected results of 2D and 3D simulations of powder compaction processes have already been reported in [12], and the model predictions have been compared to experimental data showing satisfactory agreement, see for instance Fig. 1. However, the finite-strain formulation and the numerical strategies adopted for its implementation have not been presented in [12], as that paper was aimed at providing an overview of elastoplastic coupling in powder compaction processes. In the present paper, we provide the details of the formulation and implementation, and, as an application, we study the effect of friction and initial aspect ratio on compaction of alumina powder in a cylindrical die.

2 PBG model at small strain

The small-strain PBG model [10] is briefly described below as a reference for its finite-strain version introduced in the next section, with a slight modification to the notation to make it more convenient for the subsequent extension to the finite-strain framework. The model is fully defined by specifying the free energy, the yield condition, and the plastic flow rule, and these are provided below. For the details, including justification of the specific constitutive assumptions and calibration of the model for alumina powder, refer to Piccolroaz et al. [10].

2.1 Free energy

The total strain is decomposed into the elastic and plastic parts,


and the free energy is assumed in the following form,


where the plastic strain and the forming pressure are adopted as internal variables, and . The elastoplastic coupling is here introduced through the dependence of cohesion , parameter and shear modulus on the forming pressure , namely


where denotes the Macauley brackets operator, and , , , , , , , , and are material parameters. Note that the elastoplastic coupling is related to the variation in , so that, if remains constant and equal to one, the elastic properties of the material remain unchanged during plastic flow.

The forming pressure is assumed to depend on the volumetric part of the plastic strain through the following relationship


where and are material parameters. In view of the above dependence, the free energy could formally be expressed solely in terms of the total strain and the plastic strain . However, the dependence of on is implicit, i.e., cannot be expressed as an explicit function of . It is thus convenient to keep as an internal variable with the additional constraint introduced by Eq. (6), see Section 4.3.

Note that the forming pressure and the cohesion are also used to define the yield surface (see below). The corresponding governing equations (6) and (3) are thus called ‘hardening laws’ in [10].

2.2 Inelastic strain rate

The elastoplastic coupling, introduced above through the dependence of the free energy on the plastic strain and the forming pressure , is a crucial feature of the model. As a result, the stress depends not only on the elastic strain, but also on the internal variables. Indeed, the stress is defined by


and its rate involves the contributions due to the evolution of the internal variables,


Here, is the elastic fourth-order tensor, and P describe the elastoplastic coupling, and the inelastic strain rate is defined as


The inelastic strain rate is thus not equal to the plastic strain rate , and the former will be used in the plastic flow rule, which is crucial for a consistent treatment of the elastoplastic coupling, see Bigoni [16]. The model is rate-independent, hence by the time we understand here a time-like load parameter.

2.3 Yield condition

The yield condition is defined using the Bigoni–Piccolroaz (BP) yield function [13]




and , and are the usual invariants of the stress tensor,


The forming pressure and the cohesion , which depends on through Eq. (3), define the size of the yield surface and its position along the hydrostatic axis. Parameters , , , and define the shape of the yield surface and are assumed constant.

It is seen from Eq. (11) that the BP yield function is defined infinity for , so it cannot be evaluated numerically for an arbitrary stress state, and incremental schemes employing, for instance, the return mapping algorithm cannot be applied directly. Therefore, following Stupkiewicz et al. [17], an alternative implicit yield function is used in practice, which has the same zero level set as the original yield function (i.e., ) but behaves well for arbitrary stress states, see Section 4.1.

2.4 Plastic flow rule

The flow rule is expressed in terms of the inelastic strain rate rather than the plastic strain rate , see Bigoni [16],


where is the plastic multiplier satisfying the usual complementarity conditions,


Here, is a parameter that governs non-associativity of the flow rule (), and corresponds to the associated flow rule.

3 PBG model at finite strain

The PBG model [10] has been extended to the finite-strain framework by the same group of authors in [11]. In that model, the usual multiplicative decomposition of the deformation gradient has been adopted, the free energy has been expressed in terms of the logarithmic elastic strain while keeping the same form (2) of the free energy function, and the BP yield condition has been expressed in terms of the Biot stress tensor referred to the initial configuration. With regard to the elastoplastic coupling and plastic flow rule, the Biot stress and its conjugate strain measure have been used to define the inelastic strain rate (using the general framework developed by Bigoni [16]), and that inelastic strain rate has subsequently been used in the plastic flow rule. Finally, a complete set of rate equations has been derived; however, incremental formulation and it finite element implementation have not been attempted.

In this section, a finite-strain formulation of the PBG model is introduced, which is more convenient for the finite-element implementation than the formulation of Piccolroaz et al. [11]. At the same time, the essential features of that model are preserved, namely the specific form of the free energy function, the elastoplastic coupling framework of Bigoni [16], the BP yield condition [13], and the plastic flow rule. The main difference is in the selection of the internal variables and in the choice of the stress and strain measures used to define the inelastic strain rate and to formulate the flow rule. Also, the yield condition is here expressed in terms of the second Piola-Kirchhoff stress referred to the intermediate configuration rather than in terms of the Biot stress referred to the initial reference configuration which seems more consistent with respect to the experimental testing procedures that are typically used to calibrate the model.

3.1 Free energy

The deformation gradient F is multiplicatively split into elastic and plastic parts,


and the following standard kinematic quantities are introduced,


respectively, the total and plastic right Cauchy–Green tensors, and the elastic left Cauchy–Green tensor. Furthermore, we have


In order to conveniently treat the volumetric strains, which are essential in modeling of powder compaction, the logarithmic elastic and plastic strain tensors are introduced,


where , , , and . The well-known benefit of using the logarithmic strain measure is that the volumetric strain is simply obtained as a trace of the corresponding strain tensor, and the total volumetric strain is additively decomposed into elastic and plastic contributions.

Following Piccolroaz et al. [11], the free energy can be assumed in the same functional form as in the small-strain model, Eq. (2), with the infinitesimal elastic strain simply replaced by the logarithmic strain . However, this form is not efficient in numerical implementation, and a modified free energy function is adopted in this work. For completeness, application of the original free energy of Piccolroaz et al. [11] is discussed in Appendix A.

In the modified free energy function, the volumetric behavior is described in terms of the logarithmic elastic strain , just like in the original model [11], while the shear behavior is described by the term of the neo-Hookean type formulated for the isochoric part of the elastic left Cauchy–Green tensor , namely




The right Cauchy–Green tensor C is adopted as a relevant measure of the total strain in view of the standard objectivity argument, and the plastic right Cauchy–Green tensor is adopted as an internal variable. Since elastic strains are here relatively small, the present modification of the free energy function with respect to that of [11] does not noticeably affect the actual elastic response.

The free energy (21) involves two invariants characterizing the elastic strain, and , that can be easily expressed in terms of C and . Indeed, in view of (18) and (20), we have


so that


Parameters , and in the free energy (21) are assumed to depend on the forming pressure through Eqs. (3)–(5), exactly as in the small-strain model, while the forming pressure is related to by


The above relationship is a consistent generalization of Eq. (6) to the finite deformation regime, where function is specified by Eq. (6).

3.2 Inelastic strain rate

The inelastic strain rate and subsequently the flow rule are introduced using the Green strain tensor and its conjugate stress tensor, the second Piola–Kirchhoff stress . This is a particularly convenient choice because the second Piola–Kirchhoff stress is directly obtained as the derivative of the free energy with respect to C using the following standard relationship:


Clearly, the material response is invariant with respect to the choice of a pair of conjugate strain and stress measures, see [18, 16]. Evaluation of the rate of defines the inelastic strain rate according to


where , , , and


3.3 Yield condition

The yield condition is assumed to be defined by the BP yield function (10) expressed in terms of the second Piola-Kirchhoff stress referred to the intermediate configuration,


where is the Cauchy stress, so that we have


and the yield function is now defined by Eq. (10) through the invariants of ,


The above invariants, and thus the yield function, can be explicitly expressed in terms of the second Piola-Kirchhoff stress and the plastic right Cauchy–Green tensor , the latter playing the role of a hardening variable. Indeed, we have


which is easily verified in view of the following identity holding for ,


and a similar identity holding for the respective deviators.

Note that the yield function was expressed in [11] in terms of the Biot stress tensor referred to the initial configuration. The present choice of is motivated by the typical model calibration procedures which are based on the Cauchy stress or the nominal stress referred to the intermediate configuration, see, for instance, [4]. Considering that the elastic strains are relatively small, the stress tensor is, in a sense ‘close’ to the Cauchy stress tensor , hence provides a physically sound description of the yield surface. At the same time, as shown above, the yield function depending on can be equivalently expressed solely in terms of and which is not possible if the Cauchy stress is used instead.2

3.4 Plastic flow rule

The flow rule is expressed in terms of the inelastic strain rate defined in Eq. (28) and the plastic flow direction tensor corresponding to the second Piola–Kirchhoff stress . The finite-strain counterpart of the flow rule (15) is thus the following,


with the usual complementarity conditions (16). The term in the formula for corresponds to the volumetric part of the gradient evaluated with respect to in the intermediate configuration. Note that, as in the small-strain model, the associated flow rule and thus the symmetry of the tangent operator are recovered for .

4 Finite element implementation

The essential steps involved in derivation and implementation of incremental constitutive relationships are presented in this section. The algorithmic treatment employs the commonly-used backward-Euler integration scheme and the classical return mapping algorithm, see, e.g., [1, 2]. The present computer implementation has been carried out using a symbolic code generation system AceGen [15], and the related automation is also briefly discussed below.

4.1 Implicit BP yield surface

The BP yield surface specified by Eqs. (10)–(12) is highly flexible considering the shape of its meridian and deviatoric sections. However, this comes at the cost that the original yield function is not continuous and, to enforce convexity, is defined infinity for . As a result, the BP yield function cannot be effectively evaluated for an arbitrary stress state so that the classical return mapping algorithms cannot be directly applied. A general strategy to overcome this problem, see Stupkiewicz et al. [17], is to introduce an implicitly defined yield function that has the same zero level set as the original yield function, , i.e., the same yield surface, but behaves well for arbitrary stress states. The implicit yield function formulation is followed in this work and is briefly described below. Alternative, less general approaches have been proposed in [19, 20].

Construction of a convex yield function generated by a convex yield surface is illustrated in Fig. 2. Consider the -space corresponding to a fixed Lode angle , and introduce a reference point inside the yield surface . Further, denote by the distance between the reference point and the current stress point and by the distance between the reference point and the image point that lies on the yield surface ,


and we have . The yield function is then defined by


By construction, the yield function is convex and generates a family of self-similar iso-surfaces .

Figure 2: Implicit yield function generated by a convex yield surface .

In order to evaluate the yield function for an arbitrary stress , a nonlinear equation must be solved to determine . That equation corresponds to the condition that the image point lies on the yield surface. The yield function is thus an implicit function. Consequently, its derivatives, for instance, the gradient used in the flow rule, involve the derivatives of the implicit dependence of on the stress and, possibly, also on hardening variables. The details can be found in [17].

The present implementation of the PBG model employs the above implicit formulation of the BP yield function. Accordingly, the actual use of the implicit yield function and its gradient is denoted below by a ’’ in the superscript.

4.2 Incremental flow rule

With regard to finite element implementation, the constitutive equations specified in Section 3 must be cast in an incremental form, i.e., an appropriate time integration scheme must be applied to the evolution equations for internal variables.

Using the flow rule (34), the rate constitutive equation (27) is rewritten as


which upon application of the implicit backward-Euler integration scheme yields


where and in the subscript denote, respectively, the current time and the previous time , at which the corresponding quantities are evaluated. The incremental flow rule (38) is accompanied by the complementarity conditions,


that are enforced at the end of the time increment, consistently with the backward-Euler scheme applied to integrate the rate equation (37).

Considering that arbitrary stress states can be encountered during iterative solution of the return mapping algorithm, the plastic flow direction is modified according to


so that , just like for the stresses satisfying . Of course, the flow rule is unaltered for .

Remark 1.

The simple backward-Euler scheme is usually avoided in finite-strain plasticity, in particular, in the case of plastically incompressible (or nearly incompressible) materials, e.g., in metal plasticity, and integration schemes employing the exponential map are then preferable, cf. [21, 22]. However, the present plastic flow rule, formulated in terms of the inelastic strain rate to consistently introduce the elastoplastic coupling, is not well suited for application of an exponential map integrator.

4.3 Constitutive update problem

In the constitutive update problem, given are the deformation gradient at the current time step and the internal variables and at the previous time step, and the goal is to find the current internal variables and and the plastic multiplier that satisfy the incremental flow rule (38), the complementarity conditions (39) and the constitutive relationship between and specified by Eq. (25).

The second Piola–Kirchhoff stress , which is needed to evaluate the yield function and its gradient in Eqs. (38)–(40), is defined by the free energy according to


The constitutive update problem is solved here using the classical return-mapping algorithm [1, 2]. The trial stress is first computed by assuming that the response is elastic,


for which the trial value of the yield function is evaluated according to


where the necessary invariants of are expressed in terms of and using the formulae in Eq. (32).

If then the step is elastic, and


If then the step is plastic and a set of nonlinear equations must be solved,


where the vector of unknowns comprises the internal variables and the plastic multiplier ,


and denote the components of . The local residual vector is defined as


where are the component-wise residuals corresponding to the incremental flow rule (38),


and is the equation that relates the hardening variable and the plastic strain , cf. Eq. (25),


Equation (45) is solved using the Newton method according to the following iterative scheme:


Once the constitutive update problem is solved, the first Piola–Kirchhoff stress and the consistent constitutive tangent are computed,


that are needed at the global level where the equilibrium equations are solved using the finite element method. It is reminded here that the internal variables implicitly depend on through Eq. (45), and the derivative of this implicit dependence must be accounted for when computing the consistent tangent , cf. [23].

In the present implementation, the standard return mapping algorithm specified above has been actually replaced by a more robust algorithm employing the augmented primal closest-point projection method proposed by Perez-Foguet and Armero [24]. This improved algorithm is described in Appendix B.

The model has been implemented using AceGen [14, 15], a symbolic code generation system that combines the symbolic capabilities of Mathematica (www.wolfram.com) with the automatic differentiation technique and additional tools for optimization and automatic generation of computer codes. The present formulation of incremental elastoplasticity and the structure of the constitutive update problem fit the general formulation introduced by Korelc [15], hence the automation approach developed in [15] can be directly applied to derive the necessary finite element routines. In particular, the incremental constitutive model is fully defined by specifying the local residual in terms of the internal variables , as done above, while the remaining part of the formulation remains unaltered. The details are omitted here; an interested reader is referred to [15], see also Section 2 in [25] for a concise presentation of the present automation approach in finite-strain elastoplasticity.

Application of the automatic differentiation technique implemented in AceGen results in the exact linearization of the incremental constitutive relationships which is highly beneficial for the overall performance of the Newton-based computational scheme applied to solve the global equilibrium equations.

It has been checked numerically that the consistent elastoplastic tangent corresponding to the associated flow rule, i.e., for , is not symmetric for a finite strain increment. However, the symmetry is recovered for the strain increment decreasing to zero. This shows consistency of the present incremental scheme with the rate formulation in which a symmetric elastoplastic tangent is obtained (for associative plasticity) from the Prager’s consistency condition, see [11].

5 Numerical example

An application of the model to finite element computations is presented in this section. Specifically, compaction of an alumina powder in a cylindrical die is considered with account for die-wall friction, and the effect of friction coefficient and initial aspect ratio is studied in detail. Other examples, including comparison to experimental results, can be found in [12].

Consider thus cold compression of alumina powder into a rigid cylindrical mould of radius . The powder specimen has an initial height and is compressed by a rigid punch with a maximum force corresponding to the average compaction pressure of 160 MPa. Three values of the initial aspect ratio of the alumina sample are employed, , and . Coulomb friction is assumed at the powder-mould and powder-punch contact interfaces, and three values of the friction coefficient are used, , and . Material parameters corresponding to alumina powder are provided in Table 1, see [10, 12]. The value of parameter has been increased with respect to that adopted in [10, 12] since it has been noticed that the latter may lead to exceedingly large elastic shear strains. At the same time, in compression-dominated processes, such as those considered in [10, 12] and in the present paper, the overall response is not significantly affected.

(MPa) (MPa) 1.1 2 0.1 0.19 0.9 0.5 0.383 0.124 1.8 40 (MPa) (MPa) (MPa) (MPa) (MPa) (MPa) 2.3 0.026 3.2 0.18 6 20 64 0.04 2.129 0.063

Table 1: Material parameters for alumina powder.

In the finite element implementation, an axisymmetric under-integrated four-node element employing the volumetric-deviatoric split and Taylor expansion of shape functions [26] is used for the solid, and the augmented Lagrangian method used to enforce the frictional contact constraints [27, 28].

The effect of friction coefficient and initial aspect ratio is illustrated in Fig. 3 which shows the average compacting pressure as a function of the height reduction . In each case, in addition to the compacting pressure indicated by a solid line, the corresponding average pressure at the bottom part of the mould is also shown using a dashed line of the same color. Clearly, for frictionless contact at the die wall, the two pressures would be equal one to the other, and the difference increases with increasing friction coefficient and with increasing initial aspect ratio .

    (a)     (b)     (c)

Figure 3: The average compacting pressure (solid lines) and the average pressure at the bottom part of the mould (dashed lines) as a function of the height reduction for: (a) , (b) , (c) .

The finite element mesh representing one half of the cross-section of the (axisymmetric) sample is shown in Fig. 4 for and . The undeformed mesh is shown in the left column, and the deformed meshes corresponding to the maximum compression force and different friction coefficients are shown aside. The color map, identical for all figures, indicates the resulting density of alumina powder. Due to friction, the density is nonuniform within the cross-section of the specimen, and the substantial effect of friction coefficient and initial aspect ratio on the density distribution is clearly seen in Fig. 4.

Figure 4: Axisymmetric compression into a rigid cylindrical mould: the undeformed mesh (left) and the deformed mesh corresponding to a prescribed maximum force and different friction coefficients for (upper row) and (lower row). The color map, identical for all figures, indicates the density (in g/cm).

The effect of friction, which is more pronounced for higher aspect ratios, results in reduced overall compaction so that the final height corresponding to the same prescribed maximum compression force depends on the friction coefficient. This is particularly visible for . Also, the deformation pattern is affected by the shear stresses due to friction at the die wall, which is seen in the distortion of the initially rectangular mesh.

In the PBG model, plastic hardening is governed by the volumetric plastic deformation. Nonuniform distribution of density results thus in nonuniform hardening within the sample. This is illustrated in Fig. 5 which shows the distribution of the cohesion , again indicated by the same color map for all figures. The pattern of inhomogeneity of is qualitatively similar to that shown in Fig. 4, and the same applies to the distribution of the forming pressure (not shown for brevity).

Figure 5: Axisymmetric compression into a rigid cylindrical mould: the undeformed mesh (left) and the deformed mesh corresponding to a prescribed maximum force and different friction coefficients for (upper row) and (lower row). The color map, identical for all figures, indicates the cohesion (in MPa).

6 Conclusion

A finite strain model for powder compaction has been developed through the derivation of an incremental scheme allowing the successful FE implementation of a series of ‘non-standard’ constitutive features, including: (i.) nonlinear elastic behaviour even at small strain; (ii.) coupling between elastic and plastic deformation; (iii.) pressure- and -dependent yielding; (iv.) non-isochoric flow. The extension of concepts established at small strain to the large deformation context has required a proper selection of the stress variable to be employed for the yield function and the definition of a neo-Hookean ‘correction’ to the small strain elastic potential. An incremental scheme compatible with the finite-strain formulation of elastoplastic coupling has also been developed. The numerical tests performed on the model have indicated a correct and robust behaviour of the developed code and have demonstrated the possibility of accurate simulation of the transition from a granular material to a fully dense body occurring during the cold forming of ceramic powders.


The authors gratefully acknowledge financial support from European Union FP7 project under contract number PIAP-GA-2011-286110.

Appendix A Free energy expressed in terms of the logarithmic elastic strain

In this appendix, we provide the formulation corresponding to the original form of the free energy function as proposed in [11]. Consider thus the small-strain free energy function (2) with the infinitesimal elastic strain simply replaced by the logarithmic elastic strain , viz.


The difference with respect to the free energy function (21) concerns the shear response introduced by the last two terms in Eq. (52). We also note that, in the original form [11], the logarithmic plastic strain has been adopted as the internal variable rather than the plastic right Cauchy–Green tensor . The two tensors are related by Eq. (20), or by the inverse relationship , hence both formulations are equivalent. As shown below, explicit formulae involving are readily available, hence using as the internal variable would introduce an additional and unnecessary complexity related to the tensor exponential function relating and .

The free energy (52) is expressed in terms of two invariants of the logarithmic elastic strain , namely and . Evaluation of the former does not pose any difficulties, see Eqs. (23)–(24). In order to compute the second invariant, , the logarithmic elastic strain must be computed explicitly (this has been avoided in the case of the first invariant ). In view of Eqs. (18) and (20), we have


However, this expression involves the deformation gradient F and not the right Cauchy–Green tensor C. Note that an explicit dependence on C (or equivalently on ) is needed in the present elastoplasticity framework, see Eqs. (26)–(27). Due to objectivity, the free energy function is invariant to a rigid-body rotation, hence we have


where refers to a special configuration rotated by with respect to the current configuration. It is seen that the resulting formula (54) for is rather complex as it involves the tensor logarithm function and the square root of C. Even more importantly, if the above formulation was adopted, then solution of the constitutive update problem would involve the third derivative the tensor logarithm function, which would be associated with a prohibitively high computational cost.

Considering that the elastic strains are relatively small in the materials of interest, the second invariant can be approximated with a high accuracy by exploiting the following approximation of the logarithmic strain, see [29],


which leads to


Concluding, the free energy function (52) can be directly expressed in terms of C and using Eq. (24) for and Eq. (54) or (56) for .

Appendix B Augmented primal closest-point projection

The augmented primal closest-point projection method proposed in [24] proved to significantly improve the convergence of the Newton method used for the solution of the return-mapping equations in the constitutive update problem defined in Section 4.3. The idea is to apply the augmented Lagrangian method to enforce the inequality constraints corresponding to the incremental complementarity conditions, cf. Eq. (39),


The simple treatment proposed in [24] amounts to replacing the plastic multiplier in the incremental flow rule (38) by the augmented one, ,


where . The condition is modified accordingly,


which now enforces when and otherwise. As a result, the local residual is redefined as


with given by


Although the solution sought in the plastic state () corresponds to a strictly positive , the above simple treatment leads to a significant increase of the radius of convergence of the Newton method (actually, now a semi-smooth Newton method due to the function in the definition of ) so that the computations may proceed with larger time increments.


  1. Corresponding author
  2. If the yield function is directly defined in terms of the Cauchy stress , it has to additionally depend on the total strain, for instance, on C, so that If now Prager’s consistency, , is imposed, the term yields an unsymmetrizing contribution to the tangent constitutive operator. Therefore, the choice of the Cauchy stress in the yield function leads to a model which does not fit the elastoplasticity framework of [18, 16].


  1. J. C. Simo and T. J. R. Hughes. Computational Inelasticity. Springer-Verlag, New York, 1998.
  2. E. A. de Souza Neto, D. Peric, and D. R. J. Owen. Computational Methods for Plasticity: Theory and Applications. Wiley, 2008.
  3. R. I. Borja and C. Tamagnini. Cam-clay plasticity part iii: Extension of the infinitesimal model to include finite strains. Comput. Method. Appl. M., 155:73–95, 1998.
  4. G. Meschke and W. N. Liu. A re-formulation of the exponential algorithm for finite strain plasticity in terms of Cauchy stresses. Comp. Meth. Appl. Mech. Engng., 173:167–187, 1999.
  5. A. Perez-Foguet, A. Rodriguez-Ferran, and A. Huerta. Efficient and accurate approach for powder compaction problems. Comput. Mech., 30:220–234, 2003.
  6. M. Ortiz and A. Pandolfi. A variational cam-clay theory of plasticity. Comput. Method. Appl. M., 193:27–29, 2004.
  7. M. Rouainia and D. M. Wood. Computational aspects in finite strain plasticity analysis of geotechnical materials. Mech. Res. Commun., 33:123–133, 2006.
  8. G. Frenning. Analysis of pharmaceutical powder compaction using multiplicative hyperelasto-plastic theory. Powder Technol., 172:103–112, 2007.
  9. A. Karrech, K. Regenauer-Lieb, and T. Poulet. Frame indifferent elastoplasticity of frictional materials at finite strain. Int. J. Solids Struct., 48:397–407, 2011.
  10. A. Piccolroaz, D. Bigoni, and A. Gajo. An elastoplastic framework for granular materials becoming cohesive through mechanical densification. Part I – small strain formulation. Eur. J. Mech. A/Solids, 25:334–357, 2006.
  11. A. Piccolroaz, D. Bigoni, and A. Gajo. An elastoplastic framework for granular materials becoming cohesive through mechanical densification. Part II – the formulation of elastoplastic coupling at large strain. Eur. J. Mech. A/Solids, 25:358–369, 2006.
  12. S. Stupkiewicz, A. Piccolroaz, and D. Bigoni. Elastoplastic coupling to model cold ceramic powder compaction. J. Eur. Cer. Soc., 2014. doi:10.1016/j.jeurceramsoc.2013.11.017.
  13. D. Bigoni and A. Piccolroaz. Yield criteria for quasibrittle and frictional materials. Int. J. Sol. Struct., 41:2855–2878, 2004.
  14. J. Korelc. Multi-language and multi-environment generation of nonlinear finite element codes. Engineering with Computers, 18:312–327, 2002.
  15. J. Korelc. Automation of primal and sensitivity analysis of transient coupled problems. Comp. Mech., 44:631–649, 2009.
  16. D. Bigoni. Bifurcation and instability of non-associative elastoplastic solids. In H. Petryk, editor, Material Instabilities in Elastic and Plastic Solids, volume 414 of CISM Courses and Lectures, pages 1–52. Springer-Verlag Wien, 2000.
  17. S. Stupkiewicz, R. P. Denzer, A. Piccolroaz, and D. Bigoni. Implicit yield function formulation for granular and rock-like materials. (submitted).
  18. R. Hill and J. R. Rice. Elastic potentials and the structure of inelastic constitutive laws. SIAM J. Appl. Math., 25:448–461, 1973.
  19. R. M. Brannon and S. Leelavanichkul. A multi-stage return algorithm for solving the classical damage component of constitutive models for rocks, ceramics, and other rock-like media. Int. J. Fract., 163:133–149, 2010.
  20. M. Penasa, A. Piccolroaz, L. Argani, and D. Bigoni. Integration algorithms of elastoplasticity for ceramic powder compaction. J. Eur. Cer. Soc. doi:10.1016/j.jeurceramsoc.2014.01.041.
  21. G. Weber and L. Anand. Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic–viscoplastic solids. Comp. Meth. Appl. Mech. Engng., 79:173–202, 1990.
  22. P. Steinmann and E. Stein. On the numerical treatment and analysis of finite deformation ductile single crystal plasticity. Comp. Meth. Appl. Mech. Engng., 129:235–254, 1996.
  23. P. Michaleris, D. A. Tortorelli, and C. A. Vidal. Tangent operators and design sensitivity formulations for transient non-linear coupled problems with applications to elastoplasticity. Int. J. Num. Meth. Engng., 37:2471–2499, 1994.
  24. A. Perez-Foguet and F. Armero. On the formulation of the closest-point projection algorithms in elastoplasticity—part II: Globally convergent schemes. Int. J. Num. Meth. Engng., 53:331–374, 2002.
  25. J. Korelc and S. Stupkiewicz. Closed-form matrix exponential and its application in finite-strain plasticity. Int. J. Num. Meth. Engng., 2014. doi:10.1002/nme.4653.
  26. J. Korelc and P. Wriggers. Improved enhanced strain four-node element with Taylor expansion of the shape functions. Int. J. Num. Meth. Engng., 40:407–421, 1997.
  27. P. Alart and A. Curnier. A mixed formulation for frictional contact problems prone to Newton like solution methods. Comp. Meth. Appl. Mech. Engng., 92:353–375, 1991.
  28. J. Lengiewicz, J. Korelc, and S. Stupkiewicz. Automation of finite element formulations for large deformation contact problems. Int. J. Num. Meth. Engng., 85:1252–1279, 2011.
  29. Z. P. Bazant. Easy-to-compute tensors with symmetric inverse approximating Hencky finite strain and its rate. Trans. ASME J. Eng. Mat. Technol., 120:131–136, 1998.