Improved associated flow rule for anisotropic viscoplasticity in thermoplastic polymer systems

.


Introduction
Thermoplastic systems are being increasingly used in load-bearing structures because of their attractive properties such as ease of processability, high specific stiffness and strength and, high toughness.Some of the applications in the automotive industry include front-end carriers, clutch pedals and engine mounts (Sonsino and Moosbrugger, 2008;Launay et al., 2010;Stewart, 2010).
Thermoplastics usually exhibit an orientation-dependent (anisotropic) mechanical response.This could be due to processing or the presence of reinforcements such as fibers (short-fiber, long fiber and continuous fiber composites).The anisotropy along with the time-dependent behavior of the thermoplastic result in a complex mechanical behavior, which makes the prediction of the deformation and failure challenging (Pastukhov, 2019;Kanters et al., 2018;Zhai et al., 2018;Ding et al., 2006).
Constitutive modeling of rate-independent plastic deformation of anisotropic materials has been the subject of many studies, however the combination of anisotropy with a time-dependent response has received less attention.The first yield criterion for anisotropic materials was proposed by Hill as an extension to the von Mises yield criterion (Hill, 1948(Hill, , 1971)).Since this pioneering work, many other anisotropic models have been proposed (Barlat and Lian, 1989;Gotoh, 1977;Barlat et al., 2005;Karafillis and Boyce, 1993;Bron and Besson, 2004).Some examples for anisotropic models for composite materials include (Dean et al., 2016;Notta-Cuvier et al., 2014;Mehdipour et al., 2019).
Considerable progress also has been made on constitutive modeling of the rate-dependent response of isotropic solid polymers over the last decades.One of the first approaches to model the strong ratedependence of solid polymers was to consider the polymer as a fluid with a highly stress-dependent viscosity (Haward and Thackray, 1968).These models usually applied Eyring's (1936) or Argon's (1973) theory to describe the viscosity.Different extensions of the initial 1D models to 3D finite strain deformations were proposed by various authors (Boyce et al., 1988;Buckley and Jones, 1995;Govaert et al., 2000).Effects of various factors such as temperature, hydrostatic pressure, and aging on deformation were effectively modeled through these extensions.Some more recent examples of such models are (Clarijs et al., 2019;Yang, 2019;Zhao et al., 2020).
While in material systems, such as thermoplastics, both anisotropic and time-dependent aspects of the material behavior have a significant contribution in the mechanical response and failure, constitutive modeling of the combination of these two types of behavior has been less explored.One reason being the increased difficulty and complexity of the modeling.Some examples of such models, which include a both orientation-dependent and time-dependent response, are (Huh et al., 2013;Khan and Baig, 2011;Hütter et al., 2013;Hütter and Tervoort, 2015;Senden et al., 2013).While both of these factors play a major role in deformation and failure of thermoplastics, in many thermoplastic systems, the effects of these two behaviors occur in a specific form.As will be demonstrated later in more detail, for many thermoplastic systems, the effects of strain rate and anisotropy on the flow kinetics are separable.This type of behavior is not captured by the majority of the previously proposed models, for example Huh et al. (2013).https://doi.org/10.1016/j.mechmat.2021.104087Received 6 May 2021; Received in revised form 24 September 2021; Accepted 24 September 2021 A rate-dependent anisotropic flow rule was proposed by van Erp et al. (2009), where the equivalent plastic strain rate is dependent on the equivalent Hill stress through an Eyring equation, and the derivative of the Hill stress with respect to the stress tensor results in the direction of plastic deformation.Since the plastic potential used to find the direction of the flow is the Hill equivalent stress itself, this model has been classified as an ''associated'' flow rule.This model utilizes an Eyring type equation and a stress-dependent viscosity.Van Erp et al. performed experiments on oriented polypropylene samples at various strain rates and orientations and observed that the effects of strain rate and orientation on the yield can be separated in a multiplicative manner.The developed model was aimed at capturing this observed factorizable effect of strain rate and orientation on the yield behavior.The yield stress measured in uniaxial tests at various off-axis angles  can be expressed as a product of a strain rate factor and an orientation factor: (1) The model proposed by van Erp et al. (2009) was extended in one of our previous works (Amiri-Rad et al., 2019) to include the orthotropic and nonlinear viscoelastic pre-yield regime.This was achieved by using an orthotropic Hooke's law and a multi-mode approach utilizing a spectrum of relaxation times.This model was later successfully applied to an injection-molded short-fiber composite component with a complex spatial distribution of fiber orientations (Amiri-Rad et al., 2020).While the model produced very good predictions for the tested short-fiber composite, for systems with a higher level of anisotropy an inconsistency is observed between the model predictions and the experiments.The equation for the Hill equivalent stress is fitted to the experimentally measured yield stress at different loading angles to obtain the Hill parameters.But when these parameters are used as input to the developed anisotropic elasto-viscoplastic model, the predicted yield stress at different angles deviates from the Hill fit to the experimental data.In this paper, a modification of the associated flow rule is presented, which removes this inconsistency.The proposed model is able to capture the factorizable effects of the strain rate and the orientation on yield, which is experimentally observed for various thermoplastic systems, such as short fiber, long fiber and continuous fiber reinforced, as well as oriented thermoplastics, with low to high levels of anisotropy (Pastukhov, 2019;Pastukhov and Govaert, 2021).
In the following sections, first the details of the anisotropic elastoviscoplastic constitutive model (Amiri-Rad et al., 2019) will be discussed.Then, the inconsistency, which becomes more apparent for higher levels of anisotropy, is investigated and a modification is proposed.Finally, the results of the improved model are presented and discussed.

Constitutive model
In this section, the elasto-viscoplastic constitutive model presented by Amiri-Rad et al. (2019) will be briefly discussed.More detail about the formulation and numerical implementation can be found in the original paper.This model is developed to capture the anisotropy and the rate-dependent response of thermoplastic systems and is inspired by the Eindhoven Glassy Polymer (EGP) model (Govaert et al., 2000;Tervoort et al., 1996;Klompen et al., 2005).The total stress is divided into the contributions of the driving stress resulting from intermolecular interactions of secondary bonds, and the hardening stress represented by a neo-Hookean spring which is induced by an entangled network of primary bonds (Haward and Thackray, 1968).Fig. 1 shows the mechanical analog of the constitutive model.To model the nonlinear pre-yield regime, multiple Maxwell branches (modes) with stress-dependent dashpots for viscoplastic contributions are used (Tervoort et al., 1996;van Breemen et al., 2011).The hardening stress contribution is modeled with a parallel neo-Hookean spring with modulus   .The total stress is written in terms of the different contributions as: where the hardening stress is   , and  , represents the driving stress from the th mode.The parameter  represents the total number of Maxwell branches.
In each Maxwell branch, the deformation gradient   is multiplicatively decomposed to an elastic part  , and a plastic part  , : (3) To make the decomposition unique, it is assumed that the plastic deformation includes no rigid body rotation, i.e.  , =  (Lee, 1969;Boyce et al., 1988).According to this decomposition, the deformation can be pictured as follows: First, a plastic deformation brings the material from the initial configuration  0 to the intermediate configuration , and then through elastic deformation the material reaches its current configuration   .The velocity gradient tensor   is then written as: The plastic rate of deformation tensor in the intermediate configuration is defined as the symmetric part of the plastic velocity gradient tensor: The spring in each Maxwell branch obeys an orthotropic Hooke's law and the stress tensor in each mode is found according to: where, the hat sign ̂is used to indicate the intermediate configuration.
Here, Ê, is the elastic Green-Lagrange strain tensor and Ŝ, is a pullback of the Cauchy stress with the following definition: The plastic rate of deformation tensor for each dashpot is found from an associated flow rule (van Erp et al., 2009): where σ, is the equivalent Hill stress and ̇ε , is the equivalent plastic strain rate.The stress tensor σ is obtained when the Cauchy stress tensor   is rotation neutralized to exclude the rigid body rotations (Boyce et al., 1992): where the rotation tensor  , is found from the polar decomposition of  , .The Hill stress is calculated from the following relation (Hill, 1948(Hill, , 1971)): where where σ,, are the stress components in the material coordinate system (aligned with the principal orientations of the orthotropic material).There are 6 material parameters:  11 ,  22 ,  33 ,  12 ,  23 and  13 which can be found by fitting the Hill equation to uniaxial test results at various off-axis angles.Since according to the experimental observations, these parameters do not change as the strain rate changes (factorizable effect of orientation and strain rate on yield), these parameters can be calculated from tests performed at an arbitrary strain rate.The parameters  11 ,  22 and  33 can be inferred as the ratio of yield in uniaxial tests in each material direction to a reference yield value  y, .The parameters  12 ,  23 and  13 are ratios of yield in pure shear in the respective directions to the reference shear yield stress obtained as  y, =  y, √ 3 .
The equivalent plastic strain rate ̇ε , in Eq. ( 8) is obtained from an Eyring flow equation (Eyring, 1936): where   is the viscosity of mode  and depends on the equivalent Hill stress σ and the hydrostatic pressure  as following (van Erp et al., 2009;Amiri-Rad et al., 2020): where  0, is called the zero viscosity,  0 is the characteristic stress, and  is the pressure dependency parameter.In this equation, part I takes into account the effect of stress on the flow and it is called the stress shift factor, and part II captures the influence of hydrostatic pressure on the flow kinetics.It should be noted that the equivalent Hill stress σ and the hydrostatic pressure  are based on the total driving stress and are the same for all the branches as the assumption of thermorheologically simple behavior is applied.This equation implies that plastic flow happens as soon as stress is applied to the material.This model, inspired by the EGP model for isotropic polymers, does not use an explicit yield surface and considers the material as a fluid with a high initial viscosity that decreases (approximately) exponentially with applied stress.Details on obtaining the material parameters such as the spectrum ( 0, , 4 Ĉ ), characteristic stress  0 , pressure dependency parameter , plastic anisotropy parameters  11 ,  22 ,. . .etc. from

Inconsistency in model predictions for highly anisotropic systems
In the previous section, the anisotropic elasto-viscoplastic model proposed by Amiri-Rad et al. ( 2019) was discussed.As was mentioned before, this model has been successfully applied to a 3D short-fiber composite beam (Amiri-Rad et al., 2020).In this section, the performance of the model for systems with higher levels of anisotropy than for short-fiber composites will be evaluated.As described in the previous section, the Hill parameter  11 is equal to the ratio of the yield stress in the 1 direction to a reference yield value.In this paper, 1 and 2 material directions are taken as the strongest and weakest material directions, respectively.The reference yield stress  y, can be chosen as any arbitrary value and here is taken equal to the yield stress in the 2 direction,  y,22 .Hence,  11 is equal to the following ratio and is used as an indicator of the level of anisotropy: The model described in the previous section was successfully applied to a short-fiber composite (20 wt% short E-glass fiber reinforced polycarbonate) in a previous work (Amiri-Rad et al., 2019).The  11 value for this system was between 1.27 and 1.4.Experiments on other anisotropic thermoplastics systems such as long fiber composites (LFC), unidirectional laminae (UD), solid-state drawn polypropylene tapes etc. show similar yield behavior as short fiber composites: the yield kinetics is dependent on the applied strain rate and loading orientation in a multiplicative manner (Pastukhov, 2019).This can be clearly seen in Fig. 2, where yield stress vs. strain rate curves remain parallel and only shift vertically on a double logarithmic scale when the loading angle is changed for various thermoplastic systems (Pastukhov, 2019).Hence, the model presented in the previous sections seems a good candidate for describing the behavior of these systems as well.However, when the model is applied to systems that show a higher level of anisotropy than short-fiber composites, some inconsistencies in the model predictions are observed.This will be described in more detail for a uniaxial test on oriented iPP with  11 = 9.16 (van Erp et al., 2009).Model parameters for this material are shown in Table 1.The Eyring parameters  0 and  0 are obtained for the reference angle of  = 90 • .For simplicity, only one mode is considered since the number of modes only affects the pre-yield regime and has no influence on the predicted yield stress.Also, it is assumed that there is no hardening stress and, the hardening modulus   (Fig. 1) is zero.Since the goal in this study is to investigate the inconsistencies in the predicted yield stress, which is controlled by the driving stress (Eq.( 2)), the hardening stress can be ignored.
In the following, we will derive analytically the yield stress for a uniaxial test.The material is assumed to be transversely isotropic.With this assumption, there are only 3 independent Hill parameters, which are shown in Table 1.The values for the remaining Hill parameters are:  33 =  22 ,  13 =  12 ,  23 = 1.We consider the case of uniaxial tension, with only one nonzero stress component   , in the in-plane off-axis orientation .The test is performed with a constant strain rate of   , and the yield stress is defined as the maximum point where the change in stress becomes zero,   ∕  = 0, and therefore the rate of plastic strain in the material is equal to the applied rate,  , =   .The boundary conditions can be expressed as: where  is the strain rate tensor and the symbol * shows an unprescribed value.Therefore, this equation implies that a strain rate of   is applied to the specimen in the -direction while all stress components except   are zero.For this boundary condition, the equivalent Hill stress can be written as: where  () is a factor that only depends on the angle  between the loading direction  and the 1 material direction.This function can be calculated by transforming the stress tensor (Eq.( 16)) to the material coordinate system and substituting it into Eq.( 11): As described in the previous section, the associated flow rule can be written as: where  is found by taking the derivative of the equivalent Hill stress σ with respect to the stress tensor : Therefore, the plastic strain rate in  direction can be re-written using Eqs.( 17) and ( 19): The equivalent plastic strain rate ̇ε  can be written according to the Eyring flow rule (Eq.( 13) as: where ε0 is called the rate constant, and the pressure dependence in Eq. ( 14) has been neglected for simplicity.At yield, where the ratio σ ∕ 0 is significantly larger than unity, this equation can be approximated by: Using Eqs. ( 17) and ( 23), the flow in  direction (Eq.( 21)) can be written as: Considering that at yielding, the rate of plastic flow is equal to the applied deformation rate,  , =   , the yield stress  ,y can be obtained from the above equation as: Using Eq. ( 25) and the material parameters presented in Table 1 for the oriented iPP, at the strain rate of   = 10 −3 s −1 , the yield stress for various off-axis angles  is plotted in Fig. 3(a).Also the experimental yield values and the Hill fit to them are shown in this figure.This Hill fit is obtained with Eq. ( 11) and the same parameters as shown in Table 1.
As it can be seen in Fig. 3(b), there is a noticeable inconsistency between the model prediction and the Hill fit.The maximum deviation of 13.93% is observed at  = 0 • and the deviation decreases to zero as the angle increases to 90 • , the direction of the reference yield stress.As the level of anisotropy  11 increases, the maximum deviation (deviation at  = 0 • ) of the model predictions compared to the corresponding Hill curve also increases.To show this trend, the maximum deviation at   = 0 • for a system with  12 = 1 and  22 = 1 is plotted in Fig. 4, where  11 is varied between 1 and 40.This example shows how higher levels of anisotropy lead to higher deviations.The reason for this trend will be more clear when the cause of the deviation is discussed in the next section.

A modification to the flow rule
In the previous section, it was demonstrated that the previously proposed model shows inconsistency with the Hill fit to experimental results.This inconsistency can be ignored for lower levels of anisotropy but becomes apparent for high  11 values.For the example of the uniaxially loaded oriented iPP system, it was shown that a considerable deviation of around 13.93% is observed.In this section, first the reason for this deviation will be discussed, and then a modification to the associated flow rule used in the model will be presented to remove this inconsistency.A rather complicated modification has been previously proposed by Hütter et al. (2013), to improve the predictions, however the deviations are not completely removed.In this section, we will show that at least for the two cases of uniaxial extension and pure shear, our solution completely removes the deviations with a much simpler formulation.
The reason for the inconsistency is revealed by obtaining the equivalent plastic strain rate ̇ε  from the associated flow rule.Again, we assume a uniaxial test with boundary conditions as shown in Eq. ( 16), with a constant strain rate of   applied at different off-axis angles .These boundary conditions resemble a common procedure used experimentally to obtain Hill parameters (van Erp et al., 2009;Pastukhov, 2019).The equivalent plastic strain rate ̇ε  can be found by multiplying both sides of Eq. ( 19) by the stress tensor : In the above equation the term  ∶  can be simplified.The equivalent Hill stress can be written in the following form (Han and Reddy, 2012): where   is the symmetric 4th order Hill anisotropy tensor.Using this equation,  can be written as: Hence: Substituting the above equation into Eq.( 26) results in the following energy equivalency relation: For uniaxial loading, using Eqs.( 16) and ( 17), this will result in: Hence, the equivalent plastic strain rate ̇ε  can be found as: Yield is defined as the maximum stress, where   ∕  = 0, and therefore the rate of plastic flow will be equal to the applied strain rate   =  , .While the applied deformation rate  , is the same for any off-axis angle , it can be seen from Eq. ( 32) that the equivalent plastic strain rate ̇ε  does not remain the same for all angles.
In the construction of the model, it is assumed that for a uniaxial test performed at a constant strain rate of   , the equivalent Hill stress σ is constant, and using this assumption, Hill constants are found by fitting the Hill equation to experimental data (Fig. 5-(b)).As was previously shown in the Eyring equation (Eq.( 22), also plotted in Fig. 5(a)), the equivalent plastic strain rate ̇ε  is dependent on the equivalent Hill stress σ .Therefore, a non-constant ̇ε  also means the Hill stress σ is not constant and this violates the underlying assumption of the model.As is evident from Eq. ( 32), the function  () creates the deviation in the model predictions.This function is plotted in Fig. 6 for different off-axis angles for the uniaxial test of the oriented iPP.As it is can be seen in this figure, this function has the largest deviation from unity at  = 0 • , and this deviation goes to zero as the angle approaches 90 • .This is the same trend as seen in Fig. 3

(b).
As a remedy to this problem, we propose a pre-factor  in the following form to be applied to the flow rule: where σ is the von Mises stress defined as: where    is the deviatoric part of the stress tensor   .Using this pre-factor, we can write the new flow rule as: With this modification of the flow rule, we obtain the equivalent plastic strain rate ̇ε  , by using the energy equation: Specifically in contrast to Eq. ( 32), for uniaxial loading, the equivalent plastic strain rate ̇ε  becomes: where we have used σ =   on the basis of Eq. ( 16).As it can be seen from this equation, now the equivalent plastic strain rate ̇ε  remains constant for all angles, meaning that the inconsistency that was present in Eq. ( 32) is removed.In the next section, we will apply the new flow rule to the case of uniaxial loading discussed previously.Additionally, we will apply the model also to a pure shear loading and will discuss the results of the proposed pre-factor for this case.

Results and discussion
In the previous section, the reason for the observed deviation in the model predictions was discussed and a modification to the flow rule was proposed.In this section, we will show the results of the improved model when applied to two different loading conditions: uniaxial extension and pure shear.1.

Uniaxial deformation
First, the case of uniaxial tension discussed previously will be revisited applying the improved model.The boundary conditions are the same as Eq. ( 16), where the applied deformation rate is   = 10 −3 s −1 .The material parameters are shown in Table 1.The yield stress  ,y can be obtained similar to the way Eq.( 25) was derived.For the new model with pre-factor,  ,y is obtained as: Using this equation, the yield stress  ,y for different angles  is plotted in Fig. 7, and compared to the results for the associated flow rule without pre-factor and the Hill fit to the experimental data.As can be seen in this figure, the obtained yield stress for the modified model now completely matches the Hill fit and the deviation is completely removed.Hence, with the proposed modification no inconsistency is observed.

Pure shear deformation
As the next test of the improved model, we consider pure shear loading.First, we will show that, as for uniaxial extension, the equivalent plastic strain rate does not remain constant for different loading directions in the original model.Then, we investigate the results with application of the proposed pre-factor.The boundary conditions for pure shear loading at different in-plane angles  with only non-zero stress components   and   can be written as: where the same   =   is applied for all the angles.
We assume the same material and hence the same Hill parameters as in the previous example (Table 1).We again examine whether the associated flow rule predictions match the values obtained from the Hill fit.To obtain the yield values for different angles from the Hill fit, one can write the equivalent Hill stress according to: where () has the following form: 2 () =  sin 2 (2) +  cos 2 (2) + 4 sin 2 (2) + 2 cos 2 (2). (41) Then, the yield stress can be obtained as: Using this equation, the yield stress values for an applied strain rate of   = √ 3 2 × 10 −3 s −1 for different angles are plotted in Fig. 8.This value for   is chosen to have the same value of σ as in the previous uniaxial example, but also any other value could have been chosen.
The energy equation (Eq.( 26)) for pure shear loading can be written as: which will result in the following expression for the equivalent rate: where  , equals the imposed rate   during steady yielding.Similar to the uniaxial extension case, here it is observed that the equivalent strain rate ̇ε  does not remain constant when the off-axis angle  is changed, at constant   .Therefore, again a deviation between the associated flow rule and the Hill fit is expected.To examine this, we calculate the yield stress from the associated flow rule.The yield stress values,  ,y , can be calculated by replacing ̇ε  with Eq. ( 23): in the orthotropic case too.In this paper, two loading cases of uniaxial tension and pure shear were discussed.Investigation of the model for more general loading condition can be the subject of future research.

Conclusion
In this paper, a modification to a previously developed macromechanical constitutive model for the time-dependent response of anisotropic thermoplastic systems was presented.It was shown that while the previous model can successfully capture the mechanical response of short-fiber composites, however, there is an inconsistency with the experimental results which can become considerable as the level of anisotropy increases (e.g. for long-fiber composites, continuous fiber composites and oriented polymers).The reason for this deviation was mathematically investigated and was found to be related to the viscoplastic multiplier of the associated flow rule.A modification to the associated flow rule was presented and it was shown that, for uniaxial extension and pure shear loading, the deviation is completely removed.The proposed pre-factor, while being considerably simpler than a previously proposed modification, shows better improvement.The improved model is able to capture the factorizable contributions of strain rate and material orientation on the yield stress for different thermoplastic systems with various levels of anisotropy.

Fig. 2 .
Fig. 2. Rate dependence of the strength for various iPP systems: yield strength of oriented iPP in different loading directions (data taken from van Erp et al., 2009), tensile strength of UD laminates of continuous glass fiber reinforced iPP in various loading directions (UD, data taken from Erartsin, 2020), tensile strength of oriented long fiber reinforced iPP in various loading directions (LFC, data taken from Pastukhov and Govaert, 2021), and yield strength of neat iPP (data taken from Erartsin, 2020).Lines are guides to the eye.

Fig. 3 .
Fig. 3. (a) comparison of yield stress at different off-axis angles obtained from the Hill equation and the associated flow rule.Markers show experimental data of oriented iPP (data taken from van Erp et al., 2009), (b) deviation of the yield stress obtained from the associated flow rule from the Hill equation.

Fig. 4 .
Fig. 4. Maximum deviation in the yield stress as a function of the  11 value.

Fig. 5 .
Fig. 5. (a) Eyring fit to the yield stress measurements at the reference angle of  = 90 • , (b) Hill fit to the yield stress at a constant strain rate of 10 −3 s −1 .Markers show experimental data (data taken from van Erp et al., 2009).

Fig. 6 .
Fig.6.Orientation factor  as a function of the off-axis angle  for the material parameters of Table1.