Modelling of rate-dependent inelasticity and damage in semi-crystalline polymers using an Eulerian framework

The present paper concerns the mechanical response of semi-crystalline polymers during cyclic loading, and it includes both modelling and experimental testing. The model is Eulerian in the sense that it is independent of measures of total deformation and plastic/inelastic deformations. It is able to account for such essential phenomena as strain-rate dependence, work hardening, and damage. The model was applied to uniaxial tension tests performed on high-density polyethylene (HDPE), which is a semi-crystalline polymer widely used in the industry. Two types of tests were conducted: monotonic tests, and loading–unloading tests. The model was able to reproduce the experimental results very well. The proposed model was also implemented as a UMAT in Abaqus, including an analytic tangent. The UMAT was used for simulating two 3D geometries. The implementation seems to be robust, and no convergence problems were observed.

The material in Fig. 1 is exposed to four load cycles, and in each load cycle, a higher maximum stress is reached.The stiffness of the material seems to decrease with each cycle, both for loading and unloading.The width of the hysteresis also increases with each load cycle.The apparent weakening of the material during cyclic loading is not necessarily related to actual damage in the material, but can be due to complex rearrangements of the polymer chains and associated residual stresses at the microlevel.However, at a macroscopic, phenomenological level, these phenomena can be modelled by use of damage models.
HDPE has been tested in a number of studies before.O' Connor and Findley (1962) studied the creep properties of HDPE in both tension and compression.Clausen et al. (2011) and Pawlak and Galeski (2005) studied the response of HDPE in both tension and compression.Drozdov and Christiansen (2008) investigated the response of HDPE at different strain rates and also the relaxation behaviour.In uniaxial tension, HDPE and other semi-crystalline polymers may undergo necking and stable growth of the neck.Hong et al. (2004) study the stress-strain response of a number of semi-crystalline polymers well beyond necking.Zeng et al. (2010) Fig. 1.Stress-strain response of Delrin during cyclic loading (Shojaei & Volgers, 2018).
measured the response of HDPE in both uniaxial tension and also in biaxial tension.Also, the influence of hydrostatic pressure on the viscoelastic properties of HDPE was investigated by Spitzig and Richmond (1979).The tensile properties of HDPE were also investigated under different conditions by Nikolov et al. (2006).HDPE can be manufactured with different degrees of crystallinity, and in general, the yield stress increases with the crystallinity (Peacock & Mandelkern, 1990;Popli & Mandelkern, 1987).HDPE also has the ability to undergo crazing, which was illustrated by Michler and Godehardt (2000).
The models listed above can all be said to be Lagrangian in the sense that they depend on measures of total and inelastic deformation.In the present work, an Eulerian framework for the isotropic elasticity and inelasticity of semi-crystalline polymers is applied to new experimental results.All state variables in the model are defined in the current state of the material (cf.Onat, 1968), and the evolution of elastic deformations is prescribed directly.The model builds on previous work (e.g.Eckart, 1948;Kroon et al., 2021;Kroon & Rubin, 2019, 2020;Leonov, 1976;Rubin, 1994;Rubin & Attia, 1996;Rubin & Papes, 2011).One attractive aspect of this Eulerian framework is that all state variables are, in principle, measurable in the current state of the material.Another advantage is that it is relatively straight-forward to implement numerically, as will be demonstrated in a subsequent section.A slightly different version of the present model was published in Kroon and Rubin (2023), where the model was applied to POM (polyoxymethylene).In the present work, the model is applied to experiments performed on HDPE.Two types of tests were carried out: monotonic loading at a low strain-rate, and loading-unloading tests at two higher strain-rates.The proposed model was calibrated by use of these tests.The hardening behaviour of HDPE is quite different from that of POM, and the present work therefore contains a new formulation of the hardening behaviour.Also, a 3D implicit numerical implementation, including an analytical tangent, is proposed.
The paper is organized as follows: Section 2 describes the experimental setup.The Eulerian constitutive model is presented in Section 3. Some of the main features of the proposed model are described in Section 4. In this section, the model is also compared to the experimental results.A detailed description of the numerical implementation of the model is provided in Section 5. Two numerical examples with more complicated 3D geometries are presented in Section 6. Section 7 contains a discussion and some concluding remarks.

Experiments
The material tested was a high-density polyethylene (HDPE) material, which is a semi-crystalline polymer.The HDPE used had a melt flow index of 26 g/10 min (190 • C, 2.16 kg) and a density of 0.953 g/cm 3 .Uniaxial tensile tests were performed.The test specimens were punched out from plates with a thickness of  = 0.6 mm.The geometry of the test specimen is shown in Fig. 2. The width of the test specimen in the test region was  = 6.0 mm.Thus, the initial cross-sectional area of the specimen in the test region was  0 =   = 3.6 mm 2 .
The HDPE plates had been created through an injection-moulding process.The manufacturing process renders these plates slightly orthotropic, but in the present study the material was only tested in the moulding direction.
In all tests, digital image correlation (DIC) was used to determine the local strains in the test region during the tests.A speckle pattern was sprayed onto the test specimens, and the tests were video-recorded for image analysis after the tests.During the tests with the lowest loading-rate, the thickness of the test specimen was also measured manually with a caliper.
Two types of tests were performed: (i) monotonic loading at low strain rate and up to a relatively high strain (DIC and manual thickness measurements); (ii) loading up to a prescribed strain followed by unloading to zero stress (only DIC).During the tests, force vs. time was also recorded using a load cell and a computer.For each load case, 3-8 tests were performed.Based on the measurements, the true stress-strain relation could be established.

Kinematics
The current position of a material point is denoted by , and the velocity is  = ẋ, where (•) denotes differentiation with respect to time.The velocity gradient  and the rate of deformation tensor  are defined by The dilatation,  , which is taken to be elastic, satisfies the equation In other words, inelastic deformations are taken to be incompressible.The elastic distortional deformation is determined by the symmetric, positive definite, unimodular tensor Be , which is determined by the evolution equation where  ∶  = tr( T ), and  is the second order unit tensor.The first invariant of Be is  1 = Be ∶ .The entities Be and  are similar to the elastic left Cauchy-Green deformation tensor and the elastic dilatation in the Lagrangian theory.The right hand side of (3) can be understood as follows: the first two terms are required for objectivity, the third term is needed to ensure that Be remains unimodular, and the fourth term is where the actual constitutive behaviour is prescribed.In the fourth term,  ≥ 0 is a function that determines the rate of distortional inelasticity, and the direction of inelasticity is determined by the tensor Āp , defined by This tensor satisfies the condition which ensures that Be remains unimodular (i.e.det( Be ) = 1).
Two equivalent deformation measures are also introduced.The equivalent elastic distortional deformation is given by where is the deviatoric part of Be .Henceforth, (•) ′ denotes the deviatoric part of (•), in accordance with the definition in (7).Also, the equivalent rate of total distortional deformation is defined as where  ′ is the deviatoric part of .

Strain energy and stress
The strain energy of the material, , defined per unit mass, is taken to be where  0 is the zero-stress density of the material and A. Abelen and M. Kroon Above,  ∈ [0, 1] is a damage variable, and  > 0 and  > 0 are the shear and bulk moduli of the material.Within the context of a purely mechanical theory, the rate of material dissipation requires that where  is the symmetric Cauchy stress, and  is the current density.The Cauchy stress is specified by and the von Mises stress  e is defined by For later use in the numerical algorithm, the Kirchhoff stress, , is also introduced as

Rate of inelastic deformations and damage
The inelastic deformations relax the elastic deformations and thereby the stresses.The function  is expressed as where  0 ,  0 , and  0 are constants (or functions), and  is a yield function, defined as where  is a hardening variable.The Macaulay brackets ⟨⟩ are defined by The hardening variable  is determined by the evolution equations where  is a state variable that governs the evolution of ,  ≥ 0 governs the rate of , and  s is the saturated value of the hardening rate.
The evolution of damage is taken to be where  is a positive material parameter.Although  0 is zero in a zero-stress state, it is positive when  > 0, i.e. when (20) is active.
Insertion of the expressions above into (11) yields the internal dissipation The fact that Āp ∶  ≥ 0 and the properties of the parameters and entities involved guarantee that the internal dissipation is always non-negative.

Simulation of simple load-cases
In the present section, the model responses for some uniform, simple load-cases -uniaxial tension, biaxial tension, and simple shear -are considered.For these three cases, the components of the deformation gradient take the forms Table 1 Standard model parameters.
respectively, where a Cartesian coordinate system  1 −  2 −  3 has been assumed.The components of the associated velocity gradients are For the cases of uniaxial and biaxial tension,  1 ≥ 1 is the stretch in the direction(s) of tension, and  2 ≤ 1 denotes the transverse stretch(es) in the perpendicular direction(s).In simple shear, the deformation is parameterized by the shearing .The active components of the associated Be tensors are Time is discretized using a constant time step , and discrete time steps are denoted by   , where  = 0, 1, 2, 3, … and  0 = 0.The state variable (•) at time   is denoted by (•)(  ).Eqs (3), (2), ( 18), ( 19), and (20) are discretized using an implicit scheme.Together with the boundary condition  33 = 0 these equations constitute the problem to be solved.For a given time step  +1 and a given value  1 ( +1 ) or ( +1 ), the problem is solved for the remaining unknowns using a Newton-Raphson scheme.

Model calibration
The model is calibrated by use of the results from the uniaxial tension tests.The optimal set of parameters used in the following simulations is listed in Table 1.The optimization was carried out manually.Unless otherwise stated, these are the parameters used in the simulations below.
With regard to the elastic behaviour of the material, the shear modulus and the Poisson's ratio are treated as free parameters, and the bulk modulus is then given as a function of these two entities.The shear modulus  = 350 MPa and the Poisson's ratio  = 0.4 enabled a good fit to the experimental data.This yields a bulk modulus of  = 2(1 + )∕3(1 − 2) = 1633 MPa.When plotting the results, the engineering strain is utilized, i.e.  11 =  11 − 1 and  22 =  22 − 1.
In Fig. 3(a)-(c), the model is compared to the results from the experimental study.Fig. 3(a) shows the results for monotonic loading at low loading rate ( ε11 = 0.0004∕s).The blue lines represent the experimental tests, and 9 tests were performed.As can be seen, the 9 tests are in very good agreement with each other, and the dispersion is very low.Some deviations can be observed just before failure ( 11 > 0.4), but overall the results agree very well.The green line denotes the calibrated model response.The model is able to capture the main features of the experimental curve very well.In particular, the introduction of the state variable  enables a nice prediction of the smooth transition from the elastic to the inelastic regime.Fig. 3(b) shows the lateral strains (width and thickness) as a function of the longitudinal strain for the same tests as in Fig. 3(a).The width strain ( 22 , blue lines) was measured using DIC whereas the thickness strain ( 33 , red lines) was measured by use of a caliper.The width strain shows a very consistent and well collected behaviour, whereas the thickness strain shows a significant dispersion.The 'stair-like' graphs for the thickness strain stem from the caliper measurements.Even though the caliper measurements are a bit crude, the consistent tendencies in the strain paths still make them reliable.It is clear, however, that the dispersion in the thickness strains is much larger than for the width strain.This should not be taken to mean that the measurements as such are unreliable, but rather that the material behaviour in the thickness direction is more unpredictable compared to the width direction.The response of the thickness direction is something that should be investigated further.The proposed model is isotropic, and it can be discussed if the material that the model is applied to is isotropic.Up to longitudinal strains of about 0.1, this seems to be a fair assumption, based on the results in Fig. 3(b).Beyond that is seems more questionable, even though the average thickness strain behaviour can be said to follow the width strains fairly well (but with a significant scatter).The green line is the calibrated model response, and it agrees well with the width strain, at least up to a longitudinal strain of about 0.4.The initial slope in Fig. 3(b) is governed by Poisson's ratio, and as indicated above,  = 0.4 enabled a nice fit in Fig. 3(b).
In Fig. 3(c), the outcome from loading-unloading tests at two different loading rates is shown.Results for 5 tests with ε11 = 0.029∕s (blue lines) and 3 tests with ε11 = 0.22∕s (red lines) are shown.For the lower loading rate, the loading part of the curves agree well with each other, whereas the dispersion during unloading is larger.This is mainly due to the fact that in the tests, the machine was run to a specified machine displacement after which the unloading started.Since, specimens have a slightly different prestretch at the start of the test, different specimens will experience slightly different longitudinal strain at a given machine displacement.Hence, the unloading will start at slightly different strains in different tests.
The calibrated model agrees very well with the experimental results at the lower loading rate (solid green line).At the higher loading rate, the model captures the tendency and the rate effect well, even though the model seems to underestimate the initial stiffness of the material and also the width of the hysteresis in the experimental data.
The slope of the curve during loading is higher than the slope of the curve during unloading.This is an indication of damage.This illustrates why some kind of damage component is needed in the model.
The model behaviour depends on the model parameters in a relatively complicated and inter-connected way, and the model parameters cannot be fitted in a simple sequence.But some general principles can still be discerned.The shear modulus, , is adjusted so that the initial stiffness of the material in Fig. 3(a) and (c) is reproduced, and the Poisson's ratio, , is adjusted so that the evolution of the lateral strains in Fig. 3(b) are reproduced.The parameters  0 ,  0 ,  s , and  are used for capturing the initiation and evolution of plastic hardening, and  0 and  0 model the influence of strain-rate.Finally,  enables modelling of the decreased stiffness during unloading as compared to loading.

Damage modelling
Rate-independent as well as rate-dependent inelasticity can be modelled with the present model.The hardening behaviour, the rate-dependent and the rate-independent response of this type of model have been explored in previous works, see e.g.Hollenstein et al. (2013), Kroon et al. (2021), Kroon and Rubin (2019).In this subsection, the primary focus is on the new aspect of this model, i.e. the damage aspect.The influence of damage is therefore demonstrated for some simple load-cases.associated prediction of the evolution of the damage variable, .The damage evolution for uniaxial and biaxial tension is virtually identical.The damage curve for simple shear is slightly below the curves for the other two cases.
Fig. 5(a) shows the stress-strain response during loading and unloading in uniaxial tension for the damage rate  in Table 1.With every load cycle, the accumulated damage increases, and the effective stiffness of the material decreases.In Fig. 5(b), the associated evolution of  is illustrated.The evolution of  is most pronounced right after the onset of inelastic deformations.During unloading,  remains constant, since no inelastic deformations take place.

Prerequisites
The present model was implemented as a User Material (UMAT) in Abaqus (Hibbitt et al., 2016) using implicit integration of the evolution equations.This section provides a description of the stress update algorithm and the computation of the consistent (algorithmic) tangent stiffness.

Numerical integration of the evolution equations
We consider a time step that starts at  =   , ends at  =  +1 , and where the time increment is  =  +1 −   .At  =   , the set of variables { (  ), Be (  ), (  ), (  ), (  )} is known, and the numerical integrator determines the values of these quantities at the end of the time step, i.e. { ( +1 ), Be ( +1 ), ( +1 ), ( +1 ), ( +1 )}.Following the work in Simo (1988), Papes (2013), and Rubin and Papes (2011), it is possible to develop an implicit robust and strongly objective numerical algorithm for integrating the evolution equations.In the following, entities are assumed to be associated with time  =  +1 unless otherwise indicated.
It is convenient to introduce the relative deformation gradient  r () = () −1 (  ), which satisfies the evolution equation and initial condition The relative dilatation,  r = det  r , satisfies the equations The exact solution of ( 2) is then given by The unimodular part of  r , Fr =  It can be shown that the deviatoric part of the evolution Eq. ( 3) becomes Then, the evolution in ( 33) is approximated by which guarantees an exact solution for an elastic response ( = 0).Using a backward Euler approximation of the derivative, the solution of ( 34) is written in the implicit form with Then with the help of ( 6) and ( 35), it follows that Also, the evolution of  in ( 18) is approximated as where is an implicit estimate of  based on (19).Using ( 15), the value of  can be expressed as For the evaluation of ε, an objective estimate of  is needed, and this estimate is given by Also, an implicit estimate of , based on (20), is given by When evaluating  0  0 above,  1 is needed, and  1 can be determined from the third-order polynomial

Algorithmic tangent modulus
The Kirchhoff stress is given by The updated Kirchhoff stress ultimately depends on  r , and the variation of  is given by The consistent tangent modulus, C, can then be identified as Differentiation of  yields The tangent tensor C can be expressed on the form where The details of the derivation of the tangent modulus and the definitions of D, , and  are provided in Appendix.

Preliminaries
Two numerical 3D examples are provided based on the numerical implementation in the previous section.The first problem consists of an initially bent beam that is opened up.The second problem consists of a plate with a hole that is exposed to combined stretching and shearing.The analyses are carried out as implicit dynamic simulations, and the density of the material is set to 1000 kg/m 3 .

Problem 1: Problem formulation
The implemented model is first applied to the problem illustrated in Fig. 6, i.e. an initially curved and stress-free beam.The dimensions of the beam are indicated in Fig. 6.The lower end of the bent beam is pinned/fixed.The upper, free end of the beam is exposed to a shear traction, , that initially acts upwards.The shear traction  is the true traction acting on the surface and  rotates with the deformation.The magnitude of the traction  increases linearly with time and reaches its maximum value,  0 = 1 MPa, after 1 s.After this, the magnitude is decreased linearly to 0 during 1 s.Hence, this loading-unloading procedure takes 2 s.The geometry of the beam is discretized using second-order wedge elements.

Problem 1: Results
The shear traction  causes the curved beam to open up.After 1s, the shear traction reaches its maximum amplitude, which is associated with the beam having its maximum opening.When the traction is decreased to zero, this causes the beam to spring back somewhat.Fig. 7 shows the distribution of the von Mises stress at maximum opening ( = 1 s) and after the traction has been removed ( = 2 s).When the traction reaches its maximum, the inner curved part of the beam experimences the highest stress, which is slightly above 40 MPa.The upper part of the beam experiences significant rotation.When the traction is removed, the beam springs back, and the stresses decrease.After unloading, there is, however, a residual stress field in the beam.The peak residual stress is about 20 MPa.
Fig. 8 shows the contour of the damage variable, , after unloading.The inner, curved part of the beam experiences most damage, and  reaches a value of almost 0.6 there.Large parts of the beam experience damage to some extent.It is mainly the top part, close to where the shear traction was applied, that does not show any damage at all.

Problem 2: Problem formulation
The second problem to be analysed consists of a plate with a hole in it, see Fig. 9.The dimensions of the plate are indicated in Fig. 9.The lower horizontal surface of the plate is pinned/fixed.The upper surface of the plate is exposed to a traction vector , that varies with time.This traction can be written as the sum of a normal component,  n , and a shear component,  s , i.e.  =  n  +  s , where  is the outward normal to the upper surface, and  is oriented along the upper surface and in the plane of the plate, see Fig. 9.The loading path in terms of these two components is illustrated in Fig. 10.As in problem 1,  is defined as the true traction acting on the surface, and  and  rotate with the deformation, even though the rotations in problem 2 are much smaller compared to problem 1.The ranges of  n and  s are 0 ≤  n ≤ 10 MPa and −5 MPa ≤  s ≤ 5 MPa.The loading rate is 20 MPa/s, so that the whole loading cycle takes 2s.The geometry of the plate is discretized using second-order wedge elements.

Problem 2: Results
Fig. 11 shows the outcome from the analysis in terms of distributions of the von Mises stress.Results are shown for four different stages of loading: after the negative shear traction and the positive normal traction has been applied (a), after the shear traction has become positive (b), when only the positive shear traction is left (c), and the residual stress state after all external tractions have been removed (d).When external tractions are applied, the peak stress is about 35 MPa.The locations of these peaks vary somewhat depending on the load-case, but they are located around the hole, as is to be expected.There is also a peak in the stress field at the lower corners of the plate, associated with the fact that the plate is pinned at the bottom, which causes stress concentrations in the lower corners.In the residual state, there are some remaining stresses, primarily in connection to the hole.The peak residual stress at the hole reaches about 20 MPa.
Fig. 12 shows the distributions of the damage variable, , at the same time instants as in Fig. 11.The damage variable reaches a peak value of just below 0.4, and this occurs to the left and right of the hole.As mentioned before, the stress concentrations at the lower edges also cause damage in these regions.

Discussion and concluding remarks
An isotropic model for the elastic, inelastic and damage behaviour of semi-crystalline polymers is proposed.The model contains a total number of 10 model parameters and was compared to two types of experiments: monotonic loading tests and loadingunloading tests.The experiments were well captured by the model.Based on this comparison, an optimal set of parameters could be determined.
The damage aspect of the model is a way to model the weakening response observed in many polymers during cyclic loading.This observed weakening is not necessarily due to actual damage in the material, but may be due to complex residual stress states at the micro-level in the polymer.But to introduce damage in the model is a way to model this observation phenomenologically.The model was essentially calibrated for the material response during one load cycle.More experiments would be needed to determine the model response during reloading and continued cyclic loading.
The model was implemented in Abaqus as a UMAT.Results from Abaqus and the UMAT routine were compared to Matlab simulations to confirm the correctness of the implementation.An analytic tangent was derived and implemented as well, and the nice convergence rate observed in the simulations indicate that the convergence rate is quadratic.The 3D simulations were performed as implicit dynamic simulations.However, the loading rate was so low that the inertia will not play any significant role in the equilibrium equations.The presence of damage in the model can cause loss of ellipticity of the tangent modulus.This will, in turn, cause convergence problems when solving the equations of equilibrium (or motion).Such  problems were never observed in the present simulations.The reason is probably that the rate parameter  is low enough to keep the total stiffness positive.For higher values of , loss of ellipticity and convergence problems are to be expected.
The present material exhibits significant elasticity compared to most metals.This was seen most clearly when simulating the curved beam, where the spring-back was significant when the external load was removed.
In summary, an Eulerian constitutive model for isotropic, semi-crystalline polymers has been proposed.The model is able to account for such essential phenomena as strain-rate dependence, work hardening, and damage.The model was applied to HDPE, which is a semi-crystalline polymer widely used in the industry.The model was able to reproduce the experimental results well.
The 3D implementation in Abaqus seems to be robust, and no convergence problems were observed.
Furthermore, a number of differential relations need to be established.In order to simplify the notation, a number of help variables, (•) h , are introduced.Differentiation of ( 40

Fig. 3 .
Fig. 3. Comparison between experiments (blue and red lines) and model predictions (green lines) for uniaxial tensile tests: (a) Monotonic stress-strain response with ε11 = 0.0004∕s; (b) evolution of width strain (blue lines) and thickness strain (red lines); (c) stress-strain response during loading-unloading with ε11 = 0.029∕s (blue lines and green solid line) and ε11 = 0.22∕s (red lines and green dashed line).
Fig. 4(a)  and (b) illustrates the model response for the load-cases uniaxial tension, biaxial tension, and simple shear.The model parameters in Table 1 are applied.Fig.4(a) shows the predicted stress response.The case of biaxial tension (dashed line) provides the stiffest response (which is expected), and simple shear (dotted line) is the weakest of the three load cases.Fig.4(b) shows the

Fig. 5 .
Fig. 5. (a) Stress-strain response for repeated loading and unloading in uniaxial tension; (b) predicted evolution of damage variable  for the same load-case.Loading rate: ε11 = 0.0004∕s.

Fig. 11 .
Fig. 11.Contours of the von Mises stress at 4 stages of loading.

Fig. 12 .
Fig. 12. Contours of the damage variable at 4 stages of loading.
r +  r ⊕  ) Be (  ),(55)and the notation ( ⊖ )  =     and ( ⊕ )  =     has been introduced for convenience.It then follows that But as soon as  * e exceeds (  ),  becomes positive and the state variables associated with inelasticity start evolving.In the following, e < (  ), the material response remains elastic, and  = 0.This also means that Be = B * e ,  = (  ), and  = (  ).* e > (  ) is assumed.
) gives d =  h d ε +  h d = =  h d ε +  h  h =  h +  h ε,  h = r is the deviatoric part of  r .Differentiation of (43) gives