Modeling non-linear rheology of PLLA: comparison of Giesekus and Rolie-Poly constitutive models

ABSTRACT

Rheological models for biobased plastics can assist in predicting optimum processing parameters in industrial forming processes for biobased plastics and their composites such as film blowing, or injection stretch-blow molding in the packaging industry. Mathematical descriptions of polymer behavior during these forming processes are challenging, as they involve highly nonlinear, time-, temperature-, and strain-dependent physical deformation processes in the material, and have not been sufficiently tested against experimental data in those regimes. Therefore, the predictive capability of two polymer models, a classical Giesekus and a physically-based Rolie-Poly, is compared here for extensional and shear rheology data obtained on a poly(L-lactide) (PLLA) across a wide range of strain rates of relevance to those forming processes. Generally, elongational and shear melt flow behavior of PLLA was predicted to a satisfactory degree by both models across a wide range of strain rates (for strain rates 0.05–10.0 s⁻¹), within the strain window up to 1.0. Both models show a better predictive capability for smaller strain rates, and no significant differences between their predictions were found. Hence, as the Giesekus model generally needs a smaller number of parameters, this class of models is more attractive when considering their use in computationally demanding forming simulations of biobased thermoplastics.

KEYWORDS:

• Biobased plastics
• PLLA
• rheology
• constitutive modeling
• Giesekus and Rolie-Poly models

Introduction

Bio-renewable, biodegradable, and biocompatible thermoplastic polymers, such as poly(L-lactic acid) (PLA) or poly(L-lactide) (PLLA), have been of interest to both academia and industry for more than a decade [1,2]. PLLA is of special interest in the area of sustainable and biodegradable materials (e.g. for packaging applications), and bioresorbable materials (e.g. for the use in biomedical sector). A detailed mathematical description of polymer behavior during these forming processes is a major and on-going challenge, as it involves highly nonlinear, time-, temperature-, and strain-dependent physical deformation processes in the material. An additional hurdle is that models have generally been insufficiently tested against experimental data in these regimes. Also, as the forming processes for thermoplastic polymers can also be used for processing biobased composites (e.g. nanoparticle-reinforced), the experimentally-validated rheological polymer models can also be exploited for optimizing forming process parameters and morphology for their composites [3].

A general constitutive equation for a polymer should be able to capture different types of deformation (e.g. extensional and shear) within a wide time/temperature window across large deformations, and ideally account for molecular features of the polymer (e.g. molecular weight distribution and molecular architecture). Generally, two classes of constitutive models are used to capture the non-linear viscoelastic flows of polymers: (1) the classical, phenomenological constitutive models of integral and differential type containing a relatively small number of adjustable parameters [4], and (2) constitutive equations based on the Doi-Edwards (DE) reptation theory (known as tube-based models) [5] that can contain a larger number of model parameters, especially for multi-mode representations. The single-mode tube models are found to be much more accurate compared to the phenomenological models for monodisperse linear polymers when predicting rheological experiments. On the other hand, it is known that both phenomenological and tube models make similar predictions for polydisperse linear polymers when multiple modes are used. Therefore, in the latter case the simpler classical phenomenological models may be more attractive compared to more sophisticated tube-based models. Moreover, ultimately these models must be implemented numerically e.g. using the finite-element approach, and thus the mathematical structure of the model should enable its robust numerical implementation.

This work aims at comparing those two types of the models in predicting the rheological response of biobased plastics such as PLLA in its melt state, which has not yet been carried out in the literature. Two of the current authors have previously investigated both classes of models to study rheological behavior of thermoplastic polymers in the nonlinear viscoelastic regime [6,7]. In particular, a class of constitutive models that are kinematically structured to capture geometrical non-linearity; before insertion of a description of the physical response was considered for thermoplastic polymers [6]. On the other hand, a hybrid glass-melt constitutive model, containing Rolie-Poly equations [8,9] that captures both conformational entropy elasticity and relaxation by tube diffusion in the melt [7], was developed to investigate the behavior of a linear thermoplastic in the solid- and melt-state. Here, the two approaches are used to model both the extensional and shear rheology of a biobased thermoplastic PLLA.

The work starts with a brief description of the experimental protocols for obtaining the rheological data for PLLA. This is then followed by a description of the main components of the two constitutive models, along with their parameterization using the available data for PLLA. The results of numerical simulations of extensional and shear rheological tests are then compared for the two models, and subsequently discussed in the context of their use in forming simulations of PLLA and other biobased plastics.

Experimental methods

Materials

Commercial grade polylactide (PLLA 4032D) was supplied by Nature Works, with a melt index of 19.3 g/min (210 °C, 2.16 kg) and a specific gravity of 1.24. Prior to compression molding the polymer was dried under a reduced pressure at 80 °C for 4 hours. Then, PLLA pellets were melt-compounded using a twin-screw extruder – the compounding procedure was applied to produce a reference system for organoclay-based PLLA material compounded via a similar route [3]. The molecular weight after extrusion was determined as M _n = 78,800 g mol⁻¹, with a polydispersity index M _w/M _n = 1.49 as measured by a size exclusion chromatography (SEC) with a MALLS detector in dichloromethane [3]. Samples for melt rheology were later prepared by compression molding.

Differential scanning calorimetry

Differential scanning calorimetry (DSC) measurements were performed with a TA Instruments DSC Q10 in a nitrogen atmosphere to determine the PLLA phase transition temperatures and to identify an appropriate temperature range for rheological measurements.

Isothermal crystallization for temperatures below the PLLA melting temperature was investigated following the temperature profile shown in Figure 1. The initial heat/cool/heat cycle, from 20 °C up to 190 °C, was applied to define a thermal history and was subsequently employed to determine the glass transition temperature T _g, melting temperature T _m and crystallization enthalpy ∆H. T _g was determined from the peak of the derivative heat flow with respect to temperature, whereas T _m was determined using the endothermic peak of heat flow as a function of temperature. The specimens were heated and cooled at a fixed rate of 10 °C min⁻¹ through the procedure.

xAfter the initial heat/cool/heat cycle, the temperature was reduced to an annealing temperature and held for an isothermal dwell of 20 minutes, as shown in Figure 1, matching the time required for a typical rheometric isothermal frequency sweep. The annealing temperature ranged between 160 °C and 120 °C, in 10° C steps. Following the annealing, the material was cooled to 20 °C followed by a reheating cycle to 190 °C. This reheat cycle with a defined thermal history was used to determine the degree of crystallinity formed during the isothermal dwell. Because of the significant recrystallization during the reheating cycle, this was done by integrating over the entire recrystallization and melting region. The difference between the enthalpic heat of crystallization, ∆H _c and of melting, ∆H _m after an annealing period follows ∆H = ∆H _m – ∆H _c. ∆H was used to determine the crystallinity developed during annealing, and a value of the specific heat of fusion of ∆H _fusion = 91 Jg⁻¹ was used to determine the degree of crystallinity x = ∆H/∆H _fusion [10]. A final heat/cool/heat cycle was performed to ensure that no significant degradation of the material had occurred due to the extended DSC scan period (~360 minutes).

Figure 1. Temperature profile used for DSC with isothermal annealing dwells intended to replicate the temperature history during rheometric measurements

Rheometric measurements

Linear viscoelastic shear rheology was carried out in the temperature range of T = 130–180 °C using an ARES rheometer with 25 mm diameter parallel plates. A strain sweep in the range of ω = 0.1–100 rad s⁻¹ was performed to determine the linear viscoelastic regime. A strain of 1.5% ensured linear viscoelastic response at all temperatures. Measurements of the storage (G’) and loss (G”) moduli were performed from high frequency to low frequency in logarithmically spaced steps. The obtained curves of G’ and G” versus log ω were shifted using time-temperature superposition (TTS) to create a single master curve at the reference temperature of 170 °C.

Transient non-linear shear rheology was measured using an ARES rheometer with 25 mm diameter parallel plates, and transient non-linear extensional rheology was measured using the same instrument fitted with the Elongational Viscosity Fixture (EVF) [11]. Both measurements were carried out at T = 168 °C to ensure that the material was amorphous during the test. Measurements were carried out at five strain rates, 0.1 s⁻¹, 0.5 s⁻¹, 1 s⁻¹, 5 s⁻¹ and 10 s⁻¹, to reflect typical processing conditions of industrial forming processes such as film blowing.

Constitutive models

The predictive capabilities of the phenomenological (Giesekus) and molecularly-based (Rolie-Poly) models to describe large, nonlinear, rate– and temperature–dependent deformations of polymer melts were compared to the non-linear rheological measurements for PLLA. The two models are briefly summarized below.

Giesekus model

The Giesekus model has been widely used in describing the melt behavior of polymers, where large elastic strain and finite rotations exist. The model becomes either neo-Hookean or Newtonian in the elastic and viscous limits, respectively.

The current version of the Giesekus model is given in terms of the upper convected (Oldroyd) time derivative of the isochoric stress tensor $\hat{\sigma }$

(1) $\frac{d\hat{\sigma }}{dt}=\hat{\mathbf{L}}\cdot \hat{\sigma }+\hat{\sigma }\cdot {\hat{\mathbf{L}}}^T-\frac{\hat{\sigma }}{\tau }\left(1+\frac{\alpha \cdot \hat{\sigma }}{G}\right)+\frac{2\hat{\mathbf{D}}}{G}$(1)

and in compressible form

(2) $\sigma =KlogJ\mathbf{1}+\hat{\sigma }$(2)

where $\tau =\eta /G$ is the mechanical relaxation time (characteristic time constant), defined through the shear modulus, $G$ and Newtonian shear viscosity $\eta$. The isochoric velocity gradient is denoted by $\hat{\mathbf{L}}$, while $J$ is the volume ratio, and $K$ stands for the bulk modulus.

The adjustable parameter α is helpful in fitting the model into the experimental data in the viscoelastic regime [6]. The limit α = 0 corresponds to the Upper Convected Maxwell (UCM) model, which has several weaknesses when compared with experimental data for polymer melts. Some of them can be mitigated when α > 0, such as nonphysical stress growth at high rates in extensional flows. For α = 1, the Giesekus constitutive equation is identical to the Leonov Equation [6](6) $Z=\frac{M_{\mathrm{w}}}{M_{\mathrm{e}}}$(6) .

Rolie-Poly model

The Rolie-Poly constitutive Equation [9](9) $G^{\prime }(\omega )=\sum _{i=1}^N\frac{G_i{\omega }^2{\tau }_i^2}{1+{\omega }^2{\tau }_i^2}$(9) derives from the theory of Doi and Edwards [5] and includes the process of convective constraint release (CCR).

The single-mode Rolie-Poly constitutive equation used in this work is given by:

(3) $\frac{d\hat{\sigma }}{dt}=\hat{\mathbf{L}}\cdot \hat{\sigma }+\hat{\sigma }\cdot {\hat{\mathbf{L}}}^T-\frac{1}{{\tau }_D}\left(\hat{\sigma }-1\right)-\frac{1\left(1-\sqrt{3/tr\hat{\sigma }}\right)}{{\tau }_R}\left(\hat{\sigma }+\beta {\left(\frac{tr\hat{\sigma }}{3}\right)}^{\delta }\left(\hat{\sigma }-1\right)\right)$(3)

where $\beta$ is the CCR coefficient (analogous to $c_{\nu }$ in the full theory, where $\beta =0$ [8]), ${\tau }_R$ is the Rouse time, $1$ denotes the second-order identity tensor, while ${\tau }_D$ is the reptation time. The ratio of the reptation to Rouse times defines the entanglement number [12]

(4) $Z=\frac{1}{3}\frac{{\tau }_D}{{\tau }_R}$(4)

Identification of model parameters

Validation of experimental procedure

Table 1 reports T _g, T _m, the change in heat capacity across T _g, ∆C _p, ∆H and x for PLLA. The magnitude of ∆C _p is proportional to the fraction of amorphous polymer involved in the glass transition. As the annealing temperature decreased from 160 °C to 140 °C, ∆H remains approximately constant. In contrast, ∆H, and hence x, increased a little at 130 °C and more significantly at 120 °C. This suggests that significant cold crystallization begins to occur at 120 °C for the temperature range investigated. The results indicate that the lower temperature boundary for rheometry measurements should be no less than 130 °C to ensure that PLLA crystallization does not significantly influence the rheological behavior during the measurement. The similarity of T _g recorded throughout the DSC study indicates that the thermal profile did not induce significant PLA degradation, although this could be corrected for if needed [13]. For all DSC cycles, T _m is approximately 169 °C.

Table 1. T _g, T _m, ∆C _p and the difference between the enthalpy of crystallization on heating and melting enthalpy, ∆H for PLLA for isothermal dwell temperatures of 160 °C to 120 °C

Temperature history	As received	1^st reheat	160 °C	150 °C	140 °C	130 °C	120 °C	Final reheat
T g (°C)	63.4	63.4	62.9	62.9	63.1	63.1	62.5	62.9
T m (°C)	171.4	169.7	169.0	168.8	168.8	168.8	165.5	168.8
∆C p (Jg^−1 °C^−1)	0.57	0.61	0.60	0.59	0.59	0.57	0.20	0.59
∆H (Jg^−1)	0.19	1.81	2.14	2.31	2.63	5.55	50.14	2.95
x	0.2%	2.0%	2.4%	2.5%	2.9%	6.1%	55.1%	3.2%

Linear rheology

Parameters for both models are determined with the help of the linear theory of viscoelasticity using oscillatory shear data at different temperatures. TTS was carried out first and then complemented by determining the Maxwell modes.

Time-temperature superposition

The master curve of storage and loss moduli (at the reference temperature of T _ref = 170 °C) was constructed from oscillatory shear tests within the temperature range from 130 °C to 180 °C shown in Figure 2, to produce a mastercurve presented in Figure 3. No vertical shifting was necessary. To quantify the temperature dependence of the viscoelastic properties of PLLA, horizontal shift factors, a_T , at different temperatures were fitted to the Williams-Landel-Ferry equation [14]:

(5) $\omega (T_{\mathrm{r}\mathrm{e}\mathrm{f}})=a_T\omega (T)\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}log{a}_T=\frac{-C_1\left(T-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)}{T+C_2}$(5)

The constants C ₁ and C ₂ were determined using the optimizer Reptate [15] (see Table 2 for their values), and then used to produce the viscoelastic master curve.

Figure 2. Storage and loss moduli curves for different temperatures from 130–180°C (double logarithmic scale)

Figure 3. Master curve constructed with the WLF equation (double logarithmic scale) at T _ref = 170 °C

Table 2. Parameters of the WLF equation

Parameter	Value
C 1	2.559
C 2	−39.184
T ref [°C]	170

Linear viscoelasticity

The parameters related to the theory of linear viscoelasticity are as follows: Rouse time of one entanglement segment (of length M _e), τ _e, plateau modulus, G _e, entanglement molecular weight M _e, constraint release parameter c _ν, connected with the number of chains that are needed to create one entanglement. For a known molecular weight M _w, the entanglement number can be defined as

(6) $Z=\frac{M_{\mathrm{w}}}{M_{\mathrm{e}}}$(6)

The Rouse time can then be found from

(7) ${\tau }_{\mathrm{R}}={\tau }_{\mathrm{e}}{Z}^2$(7)

while the reptation time is found from

(8) ${\tau }_{\mathrm{d}}=3{\tau }_{\mathrm{e}}{Z}^3\left(1-\frac{2C_I}{\sqrt{Z}}+\frac{C_{II}}{Z}+\frac{C_{III}}{Z^{3/2}}\right)$(8)

where $C_I=1.69$, $C_{II}=4.17$ and $C_{III}=-1.55$. The parameters are shown on the mastercurve in Figure 4, and reported in Table 3.

Figure 4. Mastercurve with illustrated molecular parameters obtained from linear viscoelasticity

Table 3. Parameters for the constitutive model from the linear viscoelasticity theory at T _ref = 170^ºC

Parameter	Value
τ e [s]	1.66∙10^−5
G e [Pa]	7.02∙10^5
M e [kDa]	6.65
cν	0.1
τ d [s]	0.111
τ R [s]	5.14∙10^−3

Determination of Maxwell modes

The mastercurve can be described using the discrete Maxwell relaxation time spectra [12]:

(9) $G^{\prime }(\omega )=\sum _{i=1}^N\frac{G_i{\omega }^2{\tau }_i^2}{1+{\omega }^2{\tau }_i^2}$(9)

(10) $G^{{ }^{\prime \prime }}(\omega )=\sum _{i=1}^N\frac{G_i\omega {\tau }_i}{1+{\omega }^2{{\tau }_i}^2}$(10)

where G_i are the relevant moduli corresponding to the relaxation times τ_i.

For the mastercurve obtained from TTS, the Maxwell mode fitting was carried out for either 2 or 5 modes – for both options, relevant G’ and G” curves are shown in Figures 5 and 6, and the modes are reported in Tables 4 and 5. More modes enable better predictions but requires additional parameters.

Figure 5. Mastercurve fitting with two discrete Maxwell modes (double logarithmic scale)

Figure 6. Mastercurve fitting with five Maxwell modes (double logarithmic scale)

Table 4. Parameters for two Maxwell modes (reference temperature 170°C)

i	τi [s]	Gi [Pa]
1	0.06185	41,893
2	0.0035737	3.4452∙10^5

Table 5. Parameters for five Maxwell modes (reference temperature 170°C)

i	τi [s]	Gi [Pa]
1	1.2734	71.634
2	0.21284	3232.3
3	0.035577	50,579
4	0.0059466	1.5936∙10^5
5	0.00099397	2.3831∙10^5

Results and discussion

Predictions of stress-strain curves are produced for non-linear shear and extensional rheology for the different models parameterized using the linear viscoelastic mastercurves. Due to the wide range of strain rates investigated in this work, each case contains plots each split into two subplots, to expose more clearly the relevant strain rate range.

Shear rheology

Experimental results (Figures 8 and 9) show a gradual buildup of stresses with increasing strain, approaching a steady state for the largest strains. Generally, the two models qualitatively capture the trends of experimental data. The predictive capability of the models decreases with increasing strain rate – they capture the data well up to 0.5 s⁻¹, slightly underpredicting the experimental results. Interestingly, for the largest strain rate studied in this work (10 s⁻¹), both models significantly overpredict the experimental data for smaller strains, but ultimately converge toward the experimental data for larger strains. It is expected that a much larger number of modes would be needed to capture the low strain region, and thus the linear viscoelastic response of the system – the response predicted at larger strains is predicted quite satisfactorily by the current models. Predictions coming from both models (for the matching number of modes) are close to each other, which makes it difficult to favor one over the other.

Figure 8. Comparison of Giesekus and Rolie-Poly models (two modes) for transient shear stress-strain curves (continuous lines – Giesekus; dashed lines – Rolie-Poly; discrete points – experimental data)

Figure 9. Comparison of Giesekus and Rolie-Poly models (five modes) for shear stress-strain curves (continuous lines – Giesekus; dashed lines – Rolie-Poly; discrete points – experimental data)

Extensional rheology

The extensional rheological data shown in Figures 10 and 11 shows again a strong dependence on strain rate. For all strain rates, the stress initially increases with strain, and gradually falls after reaching a maximum value, such phenomenon is known as stress overshoot. The two models again predict the experimental results better at smaller strain rates (up to around 1 s⁻¹). The worst predictions by both models are for the measurements at the largest strain rate of 10 s⁻¹. The predictions of both models improve when the number of modes is increased from two (2) to five (5). Once more, the two models make relatively close predictions to each other.

Figure 10. Comparison of Giesekus and Rolie-Poly models (two modes) for extensional stress-strain curves (continuous lines – Giesekus; dashed lines – Rolie-Poly; discrete points – experimental data)

Figure 11. Comparison of Giesekus and Rolie-Poly models (five modes) for extensional stress-strain curves (continuous lines – Giesekus; dashed lines – Rolie-Poly; discrete points – experimental data)

Some remarks on the application of the models to forming simulations

Although finite element simulations of forming processes for thermoplastic polymers may involve complex tool geometries such as in the injection-stretch blow molding, both models can be used without significantly affecting the simulation times. Some exceptions might arise when a significant amount of modes is required for the Rolie-Poly model, and in turn, the corresponding parameter identification procedure may become quite lengthy.

However, when forming simulations of PLLA (or other biobased thermoplastics) systems reinforced with nanoparticles are of interest, where one is interested in how process conditions may affect morphological changes at the microstructural level of the material, they typically requires application of a multiscale approach based on an RVE concept [16–19]. In those cases, finite element simulations of forming processes for PLLA-based composite materials, might prove extremely costly computationally, especially in 3D. In such cases it is expected that constitutive models with a smaller number of parameters/modes, will be much more efficient computationally. Therefore, in those cases the phenomenological models of Giesekus type may be more favored over multiple-mode and multi-parameter molecularly-based constitutive models.

Conclusions

The predictive capabilities of classical Giesekus and tube-based Rolie-Poly constitutive models in describing non-linear (extensional and shear) rheology data for a commercially available biobased polymer PLLA were compared in this work. Generally, both elongational and shear melt flow behavior of PLLA was predicted to a satisfactory degree by both models across a wide range of strain rates (for strain rate 0.05–10.0 s⁻¹), within a strain window up to 1.0. Typically, the predictions of both models were closer to the experimental data for smaller strain rates. The predictions of both models were relatively close to each other across all strain rates. An increase in the number of modes from 2 to 5 helped to make predictions closer to the data, especially for extensional rheology, but at the expense of a larger number of parameters.

In general, the Rolie-Poly model requires identification of a larger number of model parameters, compared to the Giesekus model, especially for multi-mode representations. Therefore, as the Giesekus model has only four parameters (compressible, single-mode form), this class of simple models is still quite attractive to carry out finite-element forming simulations for biobased plastics when more complex tool geometries are considered and/or polymer behavior is non-uniform as in the case of multi-component polymeric systems such as composites.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research has been financially supported by the National Science Centre (Poland) through Grant No. [DEC-2011/01/B/ST8/06492] within the programme OPUS, and Maria Skłodowska- Curie RISE action (691238).

Notes on contributors

Maja Stępień

Maja Stępień is a former PhD student in the Centre of Molecular and Macromolecular Studies of the Polish Academy of Sciences in Łódź.

Gabriel Y.H. Choong

Gabriel Choong is a Research Fellow with the Composites Research Group at the University of Nottingham. His research interests are polymer processing and characterisation of solid-state and melt-state properties of polymeric materials.

Davide S.A. De Focatiis

Davide De Focatiis is an Associate Professor at the University of Nottingham. He is a member of the Composites Research group and of the Department of Mechanical, Materials and Manufacturing Engineering. His main research interests are in rheological and solid-state properties of polymeric materials.

Łukasz Figiel

Łukasz Figiel is an Associate Professor at the University of Warwick – he is a core staff in the International Institute for Nanocomposites Manufacturing (IINM) in WMG Department of the University of Warwick, and in the Warwick Centre for Predictive Modelling (WCPM). His main research interests are in multiscale modelling of multifunctional materials.