Numerical simulation of large deformations of amorphous polymer with finite element method : Application to normal impact test

During the last decades, the part of polymeric materials considerably increased in automotive and packaging applications. However, their mechanical behaviour is difficult to predict due to a strong sensitivity to the strain rate and the temperature. Numerous theories and models were developed in order to understand and model their complex mechanical behaviour. The one proposed by Richeton et al. [Int. J. Solids Struct. 44, 7938 (2007)] seems particularly suitable since several material parameters possess a strain rate and temperature sensitivity. The aim of this study is to implement the proposed constitutive model in a commercial finite element software by writing a user material subroutine. The implementation of the model was verified on a compressive test. Next a normal impact test was simulated in order to validate the predictive capabilities of the model. A good agreement is found between the FE predictions and the experimental results taken from the literature.


Introduction
Amorphous polymeric materials exhibit interesting mechanical and physical properties.Thus, polymers are widely used for large-strain applications in the automotive or packaging industries.Below the glass transition temperature, amorphous polymers exhibit a complex behaviour which is strain rate and temperature dependent [1,2].According to the literature, the initial elastic linear response is followed by a viscoelastic response which ends at the yield point.Then, a viscoplastic response takes place characterized by a strain softening followed by a large strain hardening.
Nowadays, finite element (FE) simulations are widely used to improve geometry of structures.However, the mechanical behaviour of polymeric materials needs to be correctly reproduced.During the last decades, numerous models based on physical considerations were developed [3][4][5] in order to simulate the thermomechanical behaviour of amorphous polymer in the glassy region.A model particularly suitable is the one developed by Richeton et al. [4].In the model, the elastic modulus and yield strength are strain rate and temperature dependent.They also proposed a thermally activated flow rule based on the cooperative model initially presented in the work of Fotheringham and Cherry [6].These constitutive equations [4] allow reproducing the stress-strain curves of amorphous polymers over a wide range of strain rate and of temperature.
The aim of this study is to implement the constitutive equations proposed by Richeton et al. [4] in the commercial FE code ABAQUS/Explicit through a VUMAT subroutine.Firstly, the model formulation is detailed.Then, the integration algorithm is presented and validation a Corresponding author: chrystelle.bernard@etu.unistra.frtests at constant strain rate are performed in order to verify the correct implementation of the constitutive equation.Finally, FE simulations of an impact test were carried out.A good agreement is observed between the FE predictions and the experimental results [7].

Constitutive equations
The three-dimensional constitutive model used to describe the mechanical behaviour of amorphous polymers under finite strain is based on the work of Richeton et al. [4].According to the kinematics of finite strain, the deformation gradient tensor F is multiplicatively decomposed into an elastic, thermal and plastic components: with F p the plastic deformation gradient defined from the reference configuration to the first relaxed configuration, the plastic deformation gradient defined from the reference configuration to the first relaxed configuration, F th the thermal deformation gradient defined from the first configuration to the second relaxed configuration and F e the elastic deformation gradient defined from the second relaxed configuration to the current configuration.The polar decomposition of F, F e and F p are given by: with R the rotation tensor, U the right stretch tensor and V the left stretch tensor.The velocity gradient tensor L is given by: Incorporating Eq. (1) into Eq.(3), one obtains: EPJ Web of Conferences with The velocity gradient tensor L can be decomposed as follows: where D is the rate of deformation tensor (symmetric part) and W the spin tensor (skew-symmetric part).According to Eq. ( 6), the elastic, thermal and plastic velocity gradients are given by: In the following, two basic assumptions are made.Firstly, we assumed that the plastic flow is incompressible.Thus, We also assumed that the thermoplastic flow is irrotational.By consequence, In this work, the elastic behaviour of amorphous polymers is assumed to be linear isotropic and can be modelled using the classical Hooke's law.The Cauchy stress tensor T is defined by: with ln V e the Hencky strain and C e the fourth order stiffness tensor defined by: (11) with E the Young's modulus, ν the Poisson's ratio and δ the Kronecker-delta symbol.The strain rate and temperature dependence of the Young's modulus is given by Richeton et al. [8]: where θ is the absolute temperature, ε is the effective strain rate, T β is the β-transition temperature, T g is the glass transition temperature, T f is the flow temperature, E i are the instantaneous stiffness of the material measured by DMA at the beginning of each transition region (glassy, glass transition and rubbery region) and m i the Weibull moduli corresponding to the statistics of bond breakage.The different instantaneous moduli and transition temperatures are strain rate dependent and follow the expressions: the instantaneous moduli and the transition temperatures at the chosen reference strain rate εref , s i the strain rate sensitivity parameters, H β the β-activation energy, R the ideal gas constant and c 1 , c 2 the WLF parameters [9].According to literature [10][11][12], the Poisson's ratio of amorphous polymers is strain rate and temperature dependent in the glass transition region.In this region, the Poisson's ratio sharply increases and achieves a value close to the Poisson's ratio of a rubbery material.In order to take into account this phenomenon, the following phenomenological model is proposed by Bernard et al. [13]: where ν 0 and ν c are respectively the Poisson's ratio values at the beginning and at the end of the glass transition region, θ is related to the half of the temperature range across which the glass transition occurs and where 2ω = ω 1 / √ ln 4 with ω 1 the full width at half maximum of the Gaussian function.
The end of the elastic behaviour is characterized by the yield point.The yield strength is modelled with the equations proposed by Richeton et al. [14].Below the glass transition temperature T g , the yield strength τ y is given by: with τ i (θ ) the evolving internal shear stress at the absolute temperature θ , k the Bolzmann constant, V the activation volume, γ the effective shear strain rate, γ0 a preexponential shear strain rate constant and n describes the cooperative character of the yield strength.At the yield point, the internal stress is expressed as: with τ i (0) the internal stress at 0 K, m a material constant, α p a pressure sensitivity parameter and P the hydrostatic pressure.
After the yield point, the mechanical behaviour of amorphous polymer is characterized by a strain softening followed by an important strain hardening due to a reorientation and rearrangement of polymer chain segments.In order to describe the plastic flow, we 04043-p.2DYMAT 2015 introduced the expression of the thermoplastic strain rate tensor D thp as: (17) with β(θ ) the thermal expansion coefficient equal to β(298 K ) if θ < T g (in agreement with [15,16]) and I the identity tensor.The plastic strain rate tensor D p is given by: where T * is the deviatoric part of the driving stress tensor T * expressed in the isothermally unloading configuration, τ is the effective shear stress defined by: and ε p is the effective plastic shear rate.Below T g , the effective plastic shear rate is expressed as: where the material constants are the same as in Eq. ( 15).The strain softening is modelled by the evolution of the internal stress τ i .This stress evolves significantly from its initial value τ i y .In agreement with Boyce et al. [17], the evolution of τ i follows the phenomenological law given by: with h the softening slope and τ ps the stress referring to the preferred structural state of the material.In Eqs. ( 18) and (19), the deviatoric part of the driving stress tensor appears.It is expressed in the isothermally unloading configuration.
In the intermediate configuration, the driving stress tensor is given by: T * = R e T TR e − B (22) with B the backstress tensor which allow to take into account the important strain hardening of the material.In [4], the backstress tensor B is assumed to be coaxial to the plastic stretch tensor V p .In accordance with Arruda and Boyce [18], the principal components of B are given by: with C R the rubbery modulus, N the number of rigid chain links between entanglements, L −1 the inverse Langevin function that is estimated with the Padé approximation proposed by Cohen [19], λ i are the principal components of V p , and λ chain is the plastic stretch on a chain (or called network locking stretch): According to Richeton et al. [4], C R and N are temperature dependent parameters.More details on the model can be found in [4].

Computational procedure
The constitutive model presented above was implemented into the FE code ABAQUS/Explicit through a user material subroutine VUMAT [20].At each step time and for each integration point of the elements, the FE code provides the Cauchy stress tensor T and the user-defined state variables defined at time t, and the time increment t.In addition, the FE code estimates the strain increment tensor ε, the deformation gradient tensor F, the temperature θ and the right stretch tensor U at time t + t.At the end of each time increment, the user needs to return to ABAQUS/Explicit the updated values of Cauchy stress tensor T and of state variables.
At the beginning of the computation, at time t = 0, the material is assumed to have an elastic behaviour (F p = F th = F thp = I) and have an initial temperature θ .In small strain, amorphous polymers exhibit an elastic behaviour.Thus, the elastic deformation gradient tensor F e equals the deformation gradient tensor F. After performing the polar decomposition on F e , the Cauchy stress tensor T is computed at time t + t using Eq.(10).When the effective stress τ reaches the yield strength τ y , the plastic and thermal deformation gradient tensors deviate from the identity tensor.These tensors have to be derived from the viscoplastic flow rules, Eqs. ( 17)-( 24).In order to update the Cauchy stress tensor T t+ t , an implicit time computation scheme was used.More details about the computational procedure could be found in [13].

Computational test
To check the implementation of the constitutive model, preliminary FE simulations were performed.The material is a polycarbonate (PC) studied by Richeton et al. [4].The material was tested in uniaxial compression at 25 • C and at a constant strain rate of 10 s −1 .The material parameters can be found in the work of Richeton et al. [4].
The validation tests were performed on a square firstly meshed with one element CPE4R (continuum plane strain element with reduced integration) and next divided into 10 × 10 elements.The meshed square was subjected to a compression test with a constant strain rate under isothermal and adiabatic conditions.
A good agreement is obtained between the numerical predictions found in [4] and the FE predictions with one element.However, for a large number of elements, a distortion of the mesh is noticed as illustrated by Fig. 1(a).This phenomenon is due to an instability of the effective strain rate, εt+ t that leads to a non-convergence of the problem as shown in Fig. 1(b) for three different elements.Thus, two averaging techniques are investigated in order to stabilize the effective strain rate: geometry homogenization and time homogenization.

Geometry homogenization
In the performed compressive FE tests, the predicted effective strain rate oscillates around a mean value before diverging.From the previous remark, we assume that the effective strain rate values should be equal for all the elements with a same abscissa.Thus, at each time increment and for all the elements with the same abscissa, an average effective strain rate is computed.Thus, the geometry averaged effective strain rate εHG is defined by: with n e the number of element with the same abscissa and εi the effective strain rate for the element i.

Time homogenization
The time increment averaging method is based on the assumption that the effective strain rate evolves smoothly.The time increment homogenization method consists in averaging the effective strain rate on a given interval of time for each element.Thus, the time averaged effective strain rate εHT is defined by: where t N − t i represents the interval of time where the time homogenization is done.

Results and discussion
In Fig. 2, the mechanical response and the strain rate response obtained with 10 × 10 elements using geometry homogenization and time homogenization are compared with the ones obtained with one element.For both homogenization method used, a good agreement is found.We can note a discrepancy, in the stress-strain response (Fig. 2(a)), occurring at 80% strain in the case of the time homogenization method.The analysis of the strain rate response (Fig. 2(b)) shows that the strain rate predicted with the time homogenization geometry is slightly instable 04043-p.4 DYMAT 2015  contrary to the strain rate predicted by the geometry homogenization.The two homogenization techniques gave good results.In spite of some instabilities of the effective strain rate for the time homogenization method, this technique is more suitable since it is easier to use.

Normal impact test
A normal impact test was modelled and simulated with ABAQUS/Explicit.This impact test and the experimental results are taken from the work of Ames et al. [7].In this test a circular plate of PC is clamped and is impacted by a spherical-tipped cylindrical projectile.The mechanical behaviour of the PC was described with the constitutive model presented above.Some material properties were determined by Ames et al. [7].The other parameters were identified from stress-strain curves presented in [7] and are listed in Table 1.ABAQUS/Explicit code was used for the FE simulation using the VUMAT subroutine previously validated.In order to reduce the computation time, the axi-symmetry of the problem was used and thus only one slice of the   geometry was meshed.The PC plate has an initial radius of 101.6 mm and a thickness of 5.334 mm.The projectile, with a mass of 80 kg and a radius of 22.5 mm, impacts the plate at a velocity of 3.6 ms −1 .
The PC plate was meshed with 300 elements CAX4RT (4-node bilinear displacement and temperature, reduced integration with hourglass control element).All the nodal displacements were constrained on the right side of PC plate.The geometry of the test and of the mesh are presented in Fig. 3.The contact between the projectile and the plate is assumed frictionless.
At the end of the impact test, the reaction force on the projectile (Fig. 4) and the numerically predicted profile surface (Fig. 5) are compared with the experimental results [7].For both figures, a good agreement is found between the predicted results and the experimental ones.Moreover, a temperature increase of 44 K has been observed under the projectile.This value is closed with FE predictions found in [7].

04043-p.5
EPJ Web of Conferences However, in Fig. 5, a discrepancy between the profile shapes is observed under the projectile.In order to analyse the effect of projectile geometry, FE simulations with three different projectile diameters were carried out.The predicted results are compared with the experimental ones by Ames et al. [7] in Fig. 6.We can note that the best fit is obtained with a projectile diameter of 17.5 mm.

Conclusion
An elastic-viscoplastic model has been implemented in the commercial FE code ABAQUS/Explicit, via a VUMAT subroutine, in order to simulate the mechanical response of amorphous polymers.The numerical implementation of the constitutive model was first validated for a unique and several elements in plane strain conditions.In order to stabilize the mechanical response, homogenization techniques for averaging the strain rate have been proposed.Both homogenization methods gave good results.A normal impact test on a PC plate was modelled and simulated.A good agreement was found between the numerical prediction and the experimental results taken from the literature [7].

Figure
Figure (a) Deformed shape of a square geometry meshed with 10 × 10 elements tested in uniaxial compression at 25 • C and at 10 s −1 .(b) Strain rate-time curve for three elements of the mesh.

Figure 2 .
Figure 2. FE predictions of uniaxial compressive test at a constant strain of 10 s −1 for an unique element and 10 × 10 elements: (a) stress-strain curve and (b) effective strain rate.

Figure 3 .
Figure 3. FE mesh used in the thermo-mechanically coupled analysis of normal impact test.

Figure 4 .
Figure 4. Reaction force versus time response for the projectile.

Figure 5 .
Figure 5.Comparison between experimental and numerically predicted section profiles.

Figure 6 .
Figure 6.Influence of the projectile diameter size in the section profiles.

Table 1 .
[7]erial parameters of the proposed model for PC identified the experiments given in Ames et al.[7].