Explicit, fully implicit and forward gradient numerical integration of a hyperelasto-viscoplastic constitutive model for amorphous polymers undergoing finite deformation

Following the growing use of amorphous polymers in an expanding range of applications, interest for numerical modeling of polymer behavior has greatly increased. Together with reliable constitutive models, stable, accurate and rapid integration algorithms valid for large deformations need to be developed. Here, in the framework of hyperelasto-viscoplasticity and multiplicative split formulation, three integration algorithms (explicit, fully implicit and forward gradient) are generated for a constitutive polymer model and respective stability is investigated. The algorithms are furthermore implemented in a commercial Finite Element code and simulation of a full field tensile test is shown to capture the actual deformation behavior of polymers.


Introduction
Three main groups of polymers have emerged since the discovery of natural rubber: thermoplastics (also called glassy polymers), thermosettings and elastomers. Thermoplastics are made of long covalent chains only linked together by weak Hydrogen and Van der Waals bonds. With increasing temperature, these bonds gradually weaken and eventually break, resulting in the existence of a so-called glass transition temperature T g . Thermoplastics behave therefore like solids below T g and like viscous fluids above T g . Thermosettings on the other hand can be seen as three-dimensional networks where numerous covalent bonds tie polymer chains together; thermosettings are therefore by nature stiff. Finally, elastomers lie in between the two previous groups as usually being made of a soft thermosetting with very few covalent bonds or a blend of two thermoplastics, one being used above its glass transition temperature and the other one below. As a result, elastomers show a pure elastic response up to very high strains (∼ 400%). Thermoplastics can be further split into two sub-groups: amorphous and semi-crystalline. Crystallinity depends on the chains ability to align themselves and therefore polymers with small and regular monomers are more prone to crystallinity than large monomers with bulky side groups. Only amorphous thermoplastics are considered here.
Amorphous thermoplastics being relatively new and increasingly used, their mechanical behavior is currently a field of active research. Noticeable progress has been made by Boyce and coworkers in [1][2][3][4][5][6][7] as well as in [8], where a new mechanical model is presented and the two distinct phenomena occurring during deformation are identified, namely: chain segment rotation and molecular alignment.
Intermolecular resistance to chain segment rotation by surrounding chains occurs first and is responsible for the typical hardening-softening sequence observed at relatively low strains, giving rise to a local stress maximum. Segment rotation causes the rupture of previously mentioned Hydrogen and Van der Waals bonds as well as an increase of the local free volume which together result in micro-shear banding (comparable to Piobert-Lüders bands in metals) and softening. Being non-directional by nature this phenomenon is modeled as isotropic hardening-softening by an Eyring dashpot.
Intramolecular resistance to molecular alignment occurs at larger strains and comes from the polymer chain itself.
In principle, elastic stretching of an elastomer is similar to permanent orientation of molecular chains of glassy polymers (where entanglements replace chemical cross-links). Intramolecular resistance is therefore naturally modeled as rubber elasticity by a Langevin spring leading to kinematic hardening.
Whether a polymer behaves in a ductile or brittle manner is influenced by stability issues in the hardening-softeningrehardening sequence and is coupled to the relationship between the magnitude of isotropic and kinematic hardening. A high isotropic hardening peak combined to a weak kinematic hardening imply that large strains must develop in order to recover the peak stress level. This leads to localization and possible formation of crazes and failure in a brittle manner. Conversely, a low isotropic hardening peak combined to a strong kinematic hardening result in stabilization, as only small additional strains suffice to retrieve the peak stress level, see Fig. 1.
Subsequent to and, to some extent, in parallel with Boyce and coworkers, Anand and coworkers [9][10][11][12][13][14][15][16] have also made notable contributions to the development of the mechanical modeling of polymers. In an early paper [9], benefits of using the Hencky strain for large deformations are presented. In a series of publications devoted to metals subjected to deformations at high temperatures [10][11][12], modeling aspects subsequently applicable to polymers are developed; in particular, identification of state variables for large deformations is carefully studied in [11]. The original polymer model from Boyce is generalized in [13] to include both compressibility and thermal expansion and resulting model predictions are compared to a large amount of experimental data. In [14] the original Boyce model is reformulated within a firm thermodynamical framework, where the consequences of frame invariance on the principal of virtual power as well as consequences of the dissipation inequality on the flow rule are discussed. Note that both Boyce and Anand switch between two back-stress definitions, one based on the plastic deformation gradient (multiplicative split) and one based on the total deformation gradient, see the two one-dimensional rheological models shown in Fig. 2. The model presented in [14] is extended in [15] to a fully coupled thermomechan- ical theory, where also the response at large deformations and unloading is addressed differently. A model based on an elastic-viscoplastic mechanism and a hyperelastic network mechanism is developed in [16] to specifically address the response of thermoplastic polyurethanes. A large deformation viscoelastic-viscoplastic constitutive framework that also allows for a recovery strain at zero stress as observed in experiments is presented in [17]. In addition, constitutive models accounting for damage in thermoplastic polymers is emerging as proposed in [18,19]. Numerical algorithms are also provided in [19], as well as in [20]. In this work, numerical aspects of the constitutive model for glassy polymers developed by Boyce and co-workers and modified by Anand and Gurtin [14] are studied. Specifically, three different integration schemes: explicit, fully implicit and forward gradient are developed and their robustness examined. The resulting algorithms are then implemented into an explicit finite element code as user material routines. The outline of the paper is therefore as follows: first, a summary of the constitutive model for glassy polymers according to [14] is presented. A discussion to highlight key features of the different integration schemes is then given. Finally, a three dimensional model simulation illustrates the capability of capturing typical polymer mechanical behavior.

Constitutive model
The constitutive model to be explored here, pertinent for glassy polymers, is thoroughly described in [14]. In particular, a version suitable for large elastic stretches similar to the one proposed in [21,22] is scrutinized in the following. The material model is based on the multiplicative split of the deformation gradient F = F e F p into an elastic part, F e and a plastic part, F p (Kröner-Lee decomposition [23,24]) and it contains two internal variables (s, η) representing an intermolecular resistance to plastic flow and the local free volume, respectively. A summary of the constitutive model is given next.
The constitutive relations are expressed in the relaxed configuration in terms of the elastic Hencky strain E e and the co-rotated Kirchhoff stress T e (conjugate to the elastic Hencky strain). The relaxed configuration is taken to be rotation free, i.e. F e = R e U e = RU e , and F p = U p , as proposed in [14], where R is the rotation tensor and U e and U p is the elastic and plastic right stretch tensor, respectively. The Hencky strain is defined as E e = 1 2 log(C e ), where C e = F eT F e = (U e ) 2 is the elastic right Cauchy-Green deformation tensor. The co-rotated Kirchhoff stress tensor T e is related to the Cauchy stress tensor by the pull-back operation The equation of stress derives from an elastic free energy and can be expressed by use of the standard isotropic elasticity tensor as where E e 0 = E e − (trE e )I/3 and G and K are the shear and bulk moduli, respectively.
As discussed above, resistance to plastic flow is initially governed by isotropic hardening/softening and at higher stretches by orientation hardening. The orientation hardening is taken to be associated with plastic stretch and modeled as kinematic hardening with a back-stress S b defined as follows. From the left plastic Cauchy-Green deformation tensor B p = F p F pT , an effective plastic stretch governing the orientation hardening (and thus the back-stress modulus) is introduced as The back stress modulus μ is expressed by where L −1 is the inverse of the Langevin function defined by L(x) = coth(x) − 1/x for x > 0, μ R is a material parameter called the rubbery modulus and λ L a material parameter representing the network locking stretch. Finally, with B p 0 = B p − 1 3 tr(B p )I being the deviatoric part of B p , the kinematic hardening or back stress S b is defined as By assuming that the plastic flow is irrotational (zero plastic spin) and incompressible on the relaxed configuration, the evolution law for the plastic part of the deformation gradient can be written aṡ with D p being deviatoric. Then, D p directly corresponds to the plastic part of the velocity gradient on the relaxed configuration, and is given by a standard Mises-type flow rule that accounts for a back stress, where S 0 = dev(C e det(F e )F e−1 TF e−T ) = 2GE e 0 is the deviatoric part of the Mandel stress which is work conjugate to D p on the relaxed configuration, andτ is the equivalent shear stress defined as Here, |A| = √ A : A denotes the magnitude of a 2nd order tensor A. Thus, the magnitude of D p defined in (6) is proportional to the equivalent plastic shear strain rate, which is defined by a viscosity law as where υ 0 is a reference plastic shear strain-rate, m ∈]0; 1] is a strain-rate sensitivity parameter, s represents the intermolecular resistance to plastic flow and α is a pressure sensitivity parameter governing the influence of the mean normal stress π = tr(T)/3. The evolution laws for the two internal variables governing the isotropic hardening/softening are given by the coupled differential equations are additional material parameters. The initial conditions are given by s(0) = s 0 and η(0) = 0.
In summary, this model contains 13 model parameters: two associated with elasticity {G, K }, three associated with the rate of plastic flow {υ 0 , m, α}; two associated with kinematic (orientation) hardening {μ R , λ L }; and six associated with isotropic hardening/softening {h 0 , g 0 , s 0 , 3 Geometric interpretation in the stress space of a the generalized trapezoidal rule and b the generalized midpoint rule for inviscid plasticity. Subscripts n and n + 1 refer to two consecutive time steps and n + α to an intermediate state. r is the plastic flow direction tensor, σ the stress tensor and λ the plastic strain increment

Methods for integration of finite strain rate constitutive equations: explicit, fully implicit and forward gradient
The constitutive models presented in the previous section is expressed in terms of rates and derivatives. The resulting system of coupled differential equations is non-linear and must be integrated numerically. Such an integration inevitably raises the questions of stability, accuracy and efficiency. All integration strategies are indeed a compromise between these three issues. Two well-established integration methods: the generalized trapezoidal rule and the generalized midpoint rule are presented in [25] for inviscid plasticity. Both methods make use of known quantities at step n and a priori unknown quantities at step n + 1, where a scalar coefficient, say θ ∈ [0; 1], weights the influence of the two contributions. An illustration of the generalized trapezoidal and the midpoint rules is shown in Fig. 3 for elasto-plasticity. Depending on the value taken by θ : θ = 0, θ = 1, or 0 < θ < 1, the method is referred to as an explicit formulation, a fully implicit formulation or a gradient forward formulation. The fully implicit formulation reduces to the well-known radial return method for J2 plasticity. Note that for the extreme cases to be considered here, i.e. θ = 0 or 1, the generalized midpoint rule becomes the same as the generalized trapezoidal rule.
The two integration methods can be shown to be firstorder accurate for all θ -values, except for θ = 1/2, where second-order accuracy prevails. When it comes to stability, the midpoint rule is unconditionally stable for θ ≥ 1/2, whereas the generalized trapezoidal rule is unconditionally stable for values of θ depending on the shape of the yield surface but still always greater than 1/2. As opposed to unconditionally stable, conditionally stable means that stability depends on the incremental step size. For conditionally stable methods applied to viscoplasticity and creep, stability issues are discussed in [26], where a critical time step is put forward. This criterion is illustrated by several examples in [27] and is shown to be very much case specific: the time step increment limit depends on the viscoplastic properties of the material, the current elastic state and the current viscoplastic rate of deformation. To partially circumvent this time step limitation, a forward gradient time integration scheme is developed in [28], where the plastic strain directions are taken from step n as in an explicit method but the effective plastic strain increment is calculated from both step n and a forward gradient estimation of step n + 1. This method, which was applied in [29] to integrate the constitutive equations of the Boyce model, is explored in Sect. 4. Note that this time step limitation is purely material related and has nothing to do with the limitation imposed by the Courant criterion derived in the context of explicit time integration within finite element analysis.
Another issue of considerable importance when integrating constitutive equations in the context of finite deformation is rotation neutralization. Roughly two main strategies have been employed: one illustrated in [30,31] where quantities are rotated by the small rotation increment between two consecutive steps and one presented in [32] where all quantities are systematically rotated to an additional rotation-free configuration (pull-back), integrated and rotated back to the current configuration (push-forward). The second strategy is used in the following.
Iterative integrations schemes for hyperelasto-viscoplasticityicity also aim at, from quantities at step n, calculating quantities at step n + 1. Kinematic quantities known at the beginning of the time step are the deformation gradients: F n , F e n , F p n and F n+1 ; in addition, the (Cauchy) stress T n and the internal variables are known. As for small deformations, it seems natural to define trial quantities and, to begin with, a trial elastic deformation gradient: F e * = F n+1 F p n −1 . Carrying out the polar decomposition F e * = R e * U e * and invoking isotropy leads to two interesting results: R e n+1 = R e * and ) the plastic strain directions, D p reads D p = υ p N. The step at which D p is evaluated determines the nature of the algorithm (explicit, fully implicit or forward gradient). Using the definition of the elastic Hencky strain E e = log(U e ) and the isotropic elasticity L tensor leads directly to: T e n+1 = T e * − L : ( tD p ). Note that thanks to the use of the Hencky strain this expression is exact (no approximation was made in the derivation) and that it is precisely the same expression as for small strain plasticity. The Newton scheme is outlined below in Box 1.
The explicit formulation is by far the simplest of the three alternatives to implement. It makes straight-forward use of quantities at step n; it does not require internal Newton-Raphson iterations and it is therefore relatively fast to execute. The fully implicit formulation is as expected of much higher complexity as it contains an iterative Newton-Raphson algorithm built to converge towards the sought solution; this obviously increases the running time substantially. Finally, the forward gradient formulation is of intermediate complexity: it does not require additional internal Newton-Raphson iterations but still the expression for υ p n+1 is quite tedious and involves quite intricate algebra. Note that in this formulation, according to [28], υ p n+1 is only an approximation of the effective value at n + 1. The three algorithms are presented in detail in "Appendices A, B and C".
At first glance, the explicit formulation of "Appendix A" does not seem to follow the update scheme presented above mainly because no trial stress is calculated. The explanation is as follows: instead of first removing the rotation from F e * and then removing the plastic strain from the total strain, it is also possible to first remove the plastic contribution from F n+1 before removing rotation effects in F e n+1 . When it comes to practical implementation, two important issues should be mentioned. Firstly, both the fully implicit and the forward gradient algorithms make use of numerous 4th order tensors and of their inverse. As inversion of 4th order tensors is very time consuming and can lead to numerical instability, it was here chosen to rewrite all four order tensors as 9 by 9 and subsequently 6 by 6 matrices which are much more easily invertible. Secondly, care should be taken when approximating derivatives of the exponential of a matrix: use of less than four terms in the Taylor expansion has proven to be inaccurate enough to jeopardize the whole solution (see [33] for details).

Results
The algorithms presented above were implemented and tested in two different software: Matlab [34]-for simple homogeneous deformation fields, and Abaqus Explicit [35]-for full field solutions. The material parameters used for Polycarbonate are listed in Table 1.

Stability analysis
To analyze the characteristics of stability of the different algorithms, a homogeneous solid subjected to uniform stretching is considered. The directions of principal stretches are held constant in time. Specifically, a deformation history of relevance for uniaxial tension is examined, i.e. one stretch amplitude increases, here denoted λ 1 , whereas the other two decreases as λ 2 = λ 3 = 1/ √ λ 1 . During loading λ 1 is increased from 1 to 2. The total loading time is set to 1 s and it is divided into equally long time steps. This implies that the tensile stretch rate is constant and equal to 1/s whereas the corresponding logarithmic strain rate decreases from 1/s to 0.5/s. For figure clarity only the tensile stress is plotted in the figures. For this first problem, direct integration of the constitutive equation system is carried out for all three algorithms: explicit, fully implicit and forward gradient. This is accomplished by implementing the algorithms in Matlab. Results for the explicit algorithm are given in Fig. 5 and show a gradual loss of stability between 1800 and 800 . Results for the fully implicit algorithm are given in Fig. 6 and as expected, complete stability can be observed irrespective of the number of steps. Finally results for the forward gradient algorithm are given in Fig. 7 and show no observable loss of stability between 1800 and 1000 steps. Some instability can be spotted at 900 steps and fatal instability eventually occurs around 800 steps. Note that compared to the explicit algorithm, the forward gradient algorithm improves the solution quality up the fatal instability but it does not substantially change the limit number of steps at which fatal instability occurs.

FEM simulations of a dog bone specimen
To illustrate the ability of the algorithms developed for the polymer material model to accurately capture mechanical behavior in three-dimensional tensile loading, a standard dog bone specimen is analyzed as a second problem, see Fig. 8 for geometry and dimensions. For this full field analysis, the algorithms are implemented in Abaqus Explicit as user defined material subroutines VUMAT. Full symmetry is assumed allowing for only 1/8 of the geometry to be modelled. The rate of loading is small enough to limit effects of inertia and hence the time step imposed by the Courant condition is significantly smaller than the one required for a stable solution based on the explicit stress integration algorithms. For this reason, no difference is observed between the three different algorithms and only the results from the implicit algorithm are presented here. The force displacement curve is plotted in Fig. 9 and a series of corresponding field results (total strain in the loading direction) are shown in Fig. 10, where the degree of loading increases from left to right. It can be observed that localization sets in already at small strains but instead of catastrophic localization (as would be the case for metals) the region with smaller cross section gradually  Fig. 9. Due to symmetry, only a 1/8 of the specimen is modeled. The scale and colors quantify the total strain in the loading direction extends throughout the specimen in a stable manner since material re-hardening overcomes area reduction causing the deformed region to be stronger than the neighboring undeformed regions.

Conclusion
Three integration algorithms (explicit, fully implicit and forward gradient) are developed and presented in this paper, which all capture the specific behavior of amorphous polymer at large deformations. As expected, higher degree of complexity leads to higher stability: both the explicit and the forward gradient algorithms are only conditionally stable whereas the fully implicit algorithm is unconditionally stable. Implementation in Abaqus of these algorithms is also carried out as user defined subroutines. Due to the large number of time steps inherent to explicit solvers, simulation of a tensile test gave the same results for all subroutines. To make full use of the implicit algorithm, computation of the consistent tangent stiffness matrix would need to be performed.
One first order and two second order tensor products are used in the following.
Inner product of two tensors A and B is noted A : B and the deviatoric part of a tensor A is noted A 0 . The first and second order identity tensors are noted I and I respectively; the fourth order isotropic elasticity tensor is noted L and the Langevin function L is definined according to: L(x) = cot(x) − 1/x for x > 0. Square root, exponential and logarithm of symmetric and real second order tensors M are, as follows, defined through their spectral decomposition.

Appendix A: Explicit algorithm
Quantities known at the beginning of the time-step: -F n , T n , R e n , F p n , s n , η n (saved as internal variables) -F n+1 (deformation gradient calculated at the end of the step) (a) Quantities at time step n: Stress pull-back; Cauchy stress, T n , to corotated Kirchhoff, T e n , stress transformation: Left plastic Cauchy-Green deformation tensor: Deviatoric parts: Mean normal pressure and effective plastic stretch: Back stress modulus and back stress tensor: Equivalent shear stress: Equivalent plastic shear strain rate: Plastic strain rate tensor: (b) Update from time-step n to time-step n + 1: New plastic deformation gradient tensor: New elastic deformation gradient tensor: New elastic right Cauchy-Green deformation tensor: New elastic stretch tensor: New elastic Hencky strain tensor: New rotation tensor: New corotated Kirchhoff stress tensor: Stress push-forward; corotated Kirchhoff stress to Cauchy stress transformation: Trial elastic right Cauchy-Green deformation tensor: Trial elastic stretch tensor: Trial elastic Hencky strain tensor: Deviatoric trial corotated Kirchhoff stress tensor: Trial elastic rotation tensor: (c) Iterations: Deviatoric left Cauchy-Green deformation tensor: Mean trial normal pressure and effective plastic stretch: Back stress modulus and back stress tensor: Equivalent shear stress: Plastic strain increment tensor: Equivalent plastic shear strain increment: Deviatoric stress-strain relations residual and viscosity law residual: Residuals Taylor expansion: to zero and solving for δ λ (k) gives: where: with: (d) Update from iteration k to iteration k + 1: New plastic strain increment: New isotropic hardening internal variables: New corotated Kirchhoff stress tensor: New plastic strain increment tensor: New plastic deformation gradient tensor: (e) Final values when convergence criterion is met: (f) Stress push-forward; corotated Kirchhoff stress to Cauchy stress transformation:

Appendix C: Forward gradient algorithm
This algorithm differs from the explicit algorithm in the calculation of equivalent plastic shear strain rate, υ p , which is presented below. Quantities known at the beginning of the time-step: -F n , T n , R e n , F p n , s n , η n (saved as internal variables) -F n+1 (deformation gradient calculated at the end of the step) -θ ∈]0; 1]  New elastic deformation gradient tensor: New elastic right Cauchy-Green deformation tensor: New elastic stretch tensor: New elastic Hencky strain tensor: New rotation tensor: New corotated Kirchhoff stress tensor: Stress push-forward; corotated Kirchhoff stress to Cauchy stress transformation: