Modelling of Elongational Flow of HDPE Melts by Hierarchical Multi-Mode Molecular Stress Function Model

The transient elongational data set obtained by filament-stretching rheometry of four commercial high-density polyethylene (HDPE) melts with different molecular characteristics was reported by Morelly and Alvarez [Rheologica Acta 59, 797–807 (2020)]. We use the Hierarchical Multi-mode Molecular Stress Function (HMMSF) model of Narimissa and Wagner [Rheol. Acta 54, 779–791 (2015), and J. Rheology 60, 625–636 (2016)] for linear and long-chain branched (LCB) polymer melts to analyze the extensional rheological behavior of the four HDPEs with different polydispersity and long-chain branching content. Model predictions based solely on the linear-viscoelastic spectrum and a single nonlinear parameter, the dilution modulus GD for extensional flows reveals good agreement with elongational stress growth data. The relationship of dilution modulus GD to molecular characteristics (e.g., polydispersity index (PDI), long-chain branching index (LCBI), disengagement time τd) of the high-density polyethylene melts are presented in this paper. A new measure of the maximum strain hardening factor (MSHF) is proposed, which allows separation of the effects of orientation and chain stretching.


Introduction
Polyethylene is the most broadly used polymer, with applications ranging from plastics packaging found in daily lives to engineering plastics. Processability of commercial thermoplastics has always been a key area of research as we continue to work towards more efficient and greener manufacturing. During polymer processing, polymers are subjected to shear and extensional flow [1,2] simultaneously with extensional flow playing a dominant role in processes such as blow molding [3], fiber spinning [4], film blowing [5] and film casting [6,7]. With processing behaviors of polymers being influenced by their molar mass (M w ), polydispersity (PDI), and long-chain branching (LCB) [8], studying the effects of molecular characteristics and chain architecture on shear and extensional rheological behavior of polyethylene is important in optimizing its processability.
It is has been known that low concentrations of LCB can be beneficial to the processing characteristics of linear polyethylene [9]. However, conventional LCB detection methods that use size exclusion chromatography (SEC) and nuclear magnetic resonance (NMR) may face a challenge in detecting low levels of LCB [10][11][12]. Therefore, over the decades, various methods have been put forth to detect and determine low traces of LCB in linear polyethylene. Dilution rheology [13,14] can provide quantitative signs of LCB but requires the determination of the molecular weight distribution by SEC, and the method is restricted to long-chain branched polymers with well-defined topography, such as metallocene HDPE. The increase in the traces of LCB can be observed in the enhancement of zero-shear viscosity (η 0 ) [8,15]. However, due to the long relaxation times of branched polymers, it is not easy to get an accurate measurement of the zero-shear viscosity [16,17]. Moreover, the effects of molar mass, polydispersity, and LCB need to be separated, which is hard as the increasing contents of LCB and broadening of molecular weight distribution have similar effects on the shear rheological behavior [16][17][18]. The van Gurp-Palmen plot [19] has been used to determine molecular information of a given polymer in terms of its molecular weight, polydispersity and LCB [20][21][22][23]. Shroff and Mavridis [24] proposed a long chain branching index (LCBI) based on the ratio between the zero-shear viscosity and the intrinsic viscosity of linear polyethylene. Using similar principles, Garcia-Franco et al. [21] showed that there was a linear relationship between the LCBI to the loss angle calculated from the van Gurp-Palman plots at complex modulus of |G * | = 10 kPa for truly linear polyethylenes (hydrogenated anionically synthesized polybutadienes). Using the rule-of-thumb relation proposed by Garcia-Franco et al. [21], Morelly and Alvarez [23] proposed a simple qualitative index for LCB.
It is important to be able to directly relate the effects of molecular characteristics of a polymer to its processability, thus, to determine the extent of any improvement. This gap can be bridged by constitutive modeling of the shear and extensional rheological behaviors. A robust constitutive model capable of reaching quantitative agreement with the rheological properties of polymers can be applied in polymer processing simulation for testing the processability of polymers according to the processing conditions. This would also reduce the time and cost in physical testing of polymers to determine the suitable molecular structure of specific polymer processes and plastics applications. The Hierarchical Multimode Molecular Stress Function (HMMSF) model of Narimissa and Wagner [25][26][27][28][29] is based on the concepts of hierarchical relaxation and dynamic dilution, and takes into account the interchain tube pressure. The model accurately predicts the uniaxial and multiaxial extensional and shear rheological behaviors of linear and LCB polymer melts based on their linear-viscoelastic characterization. The HMMSF model features only a single non-linear material parameter, i.e., dilution modulus G D , for extensional flows, and an additional constraint release (CR) parameter for shear flow. The HMMSF model was extended by Narimissa and Wagner [27,30] to encompass monodisperse, bidisperse, and polydisperse linear polymer melts by means of relating the relaxation times to the Rouse stretch-relaxation times for each relaxation mode, which makes the HMMSF model suitable for all commercial polymers. The HMMSF model is a great choice to be developed into numerical simulations of polymer processes due to its requiring of a maximum of two material parameters, its mathematical simplicity, and the quantitative prediction of rheological properties of polymeric melts [29,31]. Moreover, the HMMSF model is now available for linear and nonlinear rheological modeling in the IRIS software [32].
Morelly and Alvarez [23] studied four commercial grades of high density polyethylene (HDPE) with different molecular weight, polydispersity and level of LCB to analyze the effects that molecular architecture and distribution of molecular weight have on the extensional rheology of linear polyethylene. Morelly and Alvarez [23] concluded that while molecular weight and polydispersity affect the nonlinear rheological behavior of HDPE, even low traces of LCB have a significant influence on the extensional rheology. In the same study, Morelly and Alvarez proposed a simplified LCBI that is based on the loss angle (δ) calculated from the van Gurp-Palmen plot and a maximum strain hardening factor (MSHF), which was originally developed as the ratio between zero-shear viscosity of long-chain branched and perfectly linear polyethylene melts [33].
The objective of this paper is to systematically analyze the aforementioned commercial grades of HDPEs (investigated by Morelly and Alvarez [23]) by the HMMSF model. This study will attempt to relate these molecular characteristics with the dilution modulus G D of the HMMSF model. We also propose a new method to determine the extent of the strain hardening in polymers using molecular modeling.

Materials
Four high-density polyethylene (HDPE) samples, PE-191-8, PE-127-9, PE-124-10, and PE-236-12 were studied by Morelly and Alvarez [23] regarding the effects of polymer chain architecture and molar mass distribution on their elongational flow. The samples are named in the form of PE-A-B, where A denotes the molecular weight of the polyethylene (in kg/mol) and B denotes the polydispersity index (PDI). Details of the rheological measurements of the polymers can be found in Morelly and Alvarez.

Modelling Approach
The Hierarchical Multi-mode Molecular Stress Function (HMMSF) model of Narimissa and Wagner for linear and long-chain branched (LCB) polymer melts implements the concepts of (i) hierarchical relaxation, (ii) dynamic dilution, (iii) interchain tube pressure, and (iv) convective constraint release [25][26][27][28]30]. The model consists of a well-defined set of constitutive relations based on clear and physically justified assumptions and comprises the rheology of both linear and LCB melts.
The extra stress tensor of the Hierarchical Multi-mode MSF (HMMSF) model is given as, S I A DE is the Doi and Edwards orientation tensor which is based on the assumption of independent alignment (IA) of tube segments [34], and which is five times the second order orientation tensor S, u represents the length of the deformed unit vector u at time t, and the bracket denotes an average over an isotropic distribution of unit vectors at time t , u(t ), which can be expressed as a surface integral over the unit sphere. The molecular stress functions f i = f i (t, t ) in Equation (1) are defined as the inverse of the relative tube diameters a i of each mode i, In the same way as the orientation tensor S, the molecular stress functions are functions of both the observation time t and the past time time t , when tube segments are created by reptation. The relaxation modulus G(t) of the melt can be expressed as a sum of discrete Maxwell modes with partial relaxation moduli g i and relaxation times τ i , By considering the ratio of the relaxation modulus at time t = τ i to the dilution modulus G D = G(t = τ D ), we determine the mass fractions w i of dynamically diluted linear or LCB polymer segments with relaxation time τ i > τ D , We assume that the value of w i obtained at t = τ i can be attributed to the chain segments with relaxation time τ i , and we consider segments with τ i < τ D to be permanently diluted, i.e., we fix their weight fractions at w i = 1. Although this may seem to be a very rough estimate, we could demonstrate that this is a sufficiently robust assumption to model the rheology of broadly distributed polymers, and that the modeling results are largely independent of the number of discrete Maxwell modes used to represent the relaxation modulus G(t) [25]. The evolution equation for the molecular stress functions f i = f i (t, t ) of each mode can be expressed as [25], with the initial conditions f i (t = t , t) = 1. The first term on the right hand side expresses affine deformation by stretch rate K : S with K being the velocity gradient tensor, the second term takes into account Rouse relaxation in the longitudinal direction of the tube, and the third term limits molecular stretch due to the interchain tube pressure in the lateral direction of a tube segment [35]. The topological parameter α depends on the chain architecture of the melt, and is given by, CR represents a dissipative Constraint Release (CR) term in shear flow which is zero in extensional flow. For more details about the HMMSF model, we refer to the review article [29].
The HMMSF model has also shown promising results in predicting the crystallization rate and morphology of HDPE [36] as well as extensional and shear rheology of lowviscosity polymer melts [37].

Comparison between Model Predictions and Data
All the shear and extensional data presented in this section are digitized from figures 1, 3b and 4 of the study by Morelly and Alvarez [23], and shall henceforth be referred to as M&A. The extraction of shear and extensional data was manually done using the WebPlotDigitizer [38] due to overlapping data points. The gaps seen in the presented data are due to obstruction from other data sets. The extracted extensional stress growth coefficient data are validated by comparisons between data digitized from figure 3b,c with figure 4 of [23]. The comparisons showed no change in data trend and strain hardening behavior, but a time offset was observed between the extracted data which will be further discussed in Section 3.2 of this paper.

Linear-Viscoelastic Characterization
The linear viscoelastic data (storage G and loss modulus G ) of four HDPEs were digitized from figure 1 of [23]. Figure 1 shows the best fit (green lines) of G and G (symbols) of PE-191-8, PE-127-9, PE-124-10 and PE-236-12 at T = 160 • C using the IRIS software [32]. We used the IRIS software [32] to extract the relaxation spectrum from the data and to predict the linear viscoelastic behavior. From the relaxation spectrum, we also calculated the zero-shear viscosity η 0 = ∑ i g i τ i of the HDPEs (Table 1), which is representative of the experimental frequency window, but we note that the true zero-shear viscosity may be significantly larger. We also note that the extracted relaxation spectrum using IRIS software has longer relaxation modes than the spectrum reported by M&A, and the underestimation of the weight of long relaxation modes in M&A results in poor prediction of the elongational viscosity at low strain rates [25]. From the linear viscoelastic data (storage modulus G and loss modulus G ), we used the IRIS software [32] to construct van Gurp-Palman plots (vGP), i.e., loss angle δ as a function of complex modulus G*, of the samples (Figure 2). The converging vGP plots show that neither polydispersity nor the molar mass plays a significant role in the rheological behaviours of the samples as δ → 0 (or G * → G 0 N ). Bear in mind that as the HDPE samples are semi-crystalline, their plateau moduli cannot be easily determined. As shown by Trinckle and Friedrich [20], polydispersity has an inverse relationship with the magnitude of the complex viscosity at δ = 60 • (chosen as somewhere between the terminal flow plateau, δ = 90 • , and the minima of the curves). Here, this inverse relationship is evident for all HDPE samples except PE-236-12 ( Figure 2). Furthermore, vGP curves can correlate the slope of G* to the molar mass of the polymer in the low loss angle region [20,23]. The long-chain branching index (LCBI) is another parameter that can be determined from the vGP curves at G * = 10 4 Pa [21]. This will be explained in Section 4.1.     [23]. Continuous green lines represent fit by discrete relaxation spectrum (Equation (1) and Table 1) using the IRIS software [32].  [32] using storage (G ) and loss (G ) moduli data from Morelly and Alvarez [23]. Horizontal (dashed) and vertical (dashed-dotted) lines denote the complex modulus at δ = 60 (for PDI), and loss angle at G* = 10 kPa (for LCBI, Equation (8)).

Extensional Stress Growth Coefficient
Using the data of M&A, figure 3 of [23] shows discrepancies between the elongational stress growth coefficient as a function of time for PE 124-10 reported by M&A (figure 3b of [23]) on one side, and the stress growth coefficient calculated on the other side from the extensional stress as a function of Hencky strain (data from figure 4 of [23]) by dividing with the corresponding strain rate ( . ε = 0.1, 0.5, 1/s). The same discrepancy between stress growth coefficients reported and calculated from the elongational stress was also observed for the other HDPE investigated by M&A. The discrepancy between the elongation stress growth coefficients amounts to a time offset of approximately 0.3 s (Figure 3). The same time offset of 0.3 s was observed in the data of the other three samples reported by M&A and measured by a filament-stretching rheometer (VADER1000). It is important to note that a similar time offset has also been observed in the filament-stretching rheometer (VADER1000) in our lab. Therefore, as a result of this time offset, all the elongational stress growth coefficient data presented in this paper will be determined from M&A's extensional stress as a function of Hencky strain (figure 4 of [23]).   [39] and HMMSF model (Equations (4)-(6)) using the IRIS software with dilution modulus G D of 2200, 2000, 1700 and 1200 Pa for PE 191-8, PE 127-9, PE 124-10 and PE 236-12, respectively, and using the relaxation spectra given in Table 1. Despite the limited linear-viscoelastic frequency window (a single SAOS measurement at 160 • C), excellent agreement between HMMSF model predictions and experimental data is achieved for elongational stress growth coefficient of PE 191-8, PE 124-10 and PE 236-12, while the observed deviation for PE 236-12 at Hencky strain rate ( . ε = 0.1/s) is most likely due to a measurement issue as it sits below the linear-viscoelastic envelope (LVE) as compared to the other two Hencky strain rates ( . ε = 0.5 and 1/s). We observe fair quantitative agreement between the HMMSF model predictions and the experimental data even at the lowest strain rate ( . ε = 0.1/s), which would be most affected by the underestimation of the weight of longest relaxation modes. Generally, it is important to extend the experimentally accessible frequency window by using time-temperature superposition of multiple SAOS measurements at higher temperatures and extract the relaxation spectra from the mastercurve. Nonetheless, the HMMSF model can provide satisfying quantitative agreement between predictions and experimental data despite the limited linear viscoelastic frequency window. The overshoots in several elongational stress growth coefficients seen in Figure 4 are most likely not true maxima and they are rather signs of inhomogeneous deformation which will be further discussed in Section 4.2.    [39] and HMMSF model (Equations (4)-(6)) using again the IRIS software with dilution moduli G D of 2200, 2000, 1700 and 1200 Pa for PE 191-8, PE 127-9, PE 124-10 and PE 236-12, respectively, and using the relaxation spectra given in Table 1. Despite the mentioned limited linearviscoelastic frequency window, excellent agreement between HMMSF model predictions and experimental data is achieved for the extensional stress of all four HDPEs. Once again, the observed deviation for PE 236-12 at Hencky strain rate ( . ε = 0.1/s) is most likely due to a measurement issue as it seems to fit closer to the tube model of Doi-Edwards which is for non-stretching monodisperse linear melts.  Overall, the predictions of the HMMSF model are in quantitative agreement with the experimental data of M&A despite the limited linear-viscoelastic frequency window.

Relationship between Dilution Modulus, Polydispersity, and Long Chain Branching Index
The relationship between the dilution modulus of the HDPE samples obtained by fitting the single nonlinear parameter of the HMMSF model (G D ) to the experimental data, and the polydispersity index (PDI) of the samples is illustrated in Figure 6. It is evident that there is an inverse relation between PDI and G D , meaning that increasing PDI leads to the increase in dynamic dilution, i.e., a decreasing dilution modulus. This is in line with the basic assumption of the HMMSF model: A broader molecular weight distribution with an increased concentration of short chains leads to enhanced permanent dilution, which in turn reduces the value of the dilution modulus.  Figure 6 shows that except for the HDPE with PDI = 12 (i.e., PE-236-12), the LCBI of the samples increases as their PDI increase. This observation may indicate the dependence of long-chain branching (LCB) on the polydispersity. This will be investigated in more detail in the following parts.
The maximum strain hardening factor (MSHF), as used by M&A, follows the zeroshear viscosity enhancement factor of Stadler and Münstedt [33], and represents the ratio between the elongational stress growth factor σ E (t) of the polymer and the time-dependent linear-viscoelastic envelope of the stress growth factor η LVE (t) = 3η 0 (t), Using the values of figure s3 of [23], Figure 7 shows the relation between LCBI, MSHF, and PDI of the samples at strain rates of 0.1, and 1 s −1 . Except for PE-236-12 @ . ε = 0.1/s, a direct relationship is seen between MSHF and LCBI which emphasizes the role of traces of LCB in the strength of the strain hardening behavior of linear PE during elongational flow. The mentioned exception is in line with the erroneous stress growth coefficient data of the same sample reported in Figure 4 of Section 3.2. The effect of LCBs, when there are only traces of LCBs and these are predominantly of a star shape, is represented in the LVE spectrum by broadening the spectrum, but this has no additional effect on the elongational viscosity (in addition to the effect that longer relaxation times lead to higher stretch). Therefore, it may be of interest to calculate the second moment of the relaxation spectrum, i.e., the disengagement (or longest) relaxation time τ d . The disengagement time (τ d ), and the linear-viscoelastic zero-deformation viscosity in elongational flow (η LVE ) can be calculated from the relaxation spectrum as, We note again that these quantities are based on the relaxation spectra obtained from SAOS measurements in the experimental window. Figure 8 shows the relationship between the PDI, LCBI (Equation (8)), and τ d and η 0 (Equation (10)). There is a direct relation between all those rheological parameters (except the mentioned case for the PDI of PE-236-12). As both disengagement time and LCBI are related directly, we can confirm the validity of Equation (8) as estimation of the amount of LCB in linear melts. Henceforth, we will use the disengagement time as the measure of LCB. Figure 8. Relationship between PDI, LCBI (Equation (8)), and τ d and η 0 (Equation (10)).
As already mentioned, the minor maxima observed in the extensional data ( Figure 4) are not true maxima and they are rather signs of inhomogeneous deformation (i.e., necking) at high strains as opposed to being a material property of HDPE (for more detail, see [40][41][42]). Using the HMMSF model, and considering the fact that the occurrence of stress overshoots (maxima) is unlikely in linear melt with only traces of LCB, we can define M&A's MSHF (Equation (9)) as the ratio between the steady-state elongational viscosity and the LVE, η E ( . ε)/η LVE , as both values are readily accessible in the IRIS software ( Figure 4). In order to determine the steady-state elongational viscosity and LVE, we can use the corresponding values of the stress growth coefficient and LVE at Hencky strain ε H = 5 which is experimentally accessible by filament stretching rheometers, Figure 9 shows the comparison between M&A's MSHF (figure s3 of [23]) and those obtained from the HMMSF model, both as functions of the disengagement time and strain rate.  (9)) and those generated by the HMMSF model (circle, Equation (11)).
Ignoring the mentioned error (PE 236-12 @ . ε = 0.1/s), both methods show a direct relationship between LCB and the MSHF independent from the PDI. In contrast to M&A's values, the HMMSF model (Equation (11)) shows that the MSHF is rather strain rate insensitive. The observed rate-insensitivity of the MSHF will be discussed in more detail next.

A New Approach for Evaluating Strain Hardening
We propose a new method for investigating the strength of the strain hardening behavior in polymers based on the molecular modelling. The tube model of Doi-Edwards [39], which was originally developed for monodisperse linear melts, assumes the tension in the deformed chain is equal to its equilibrium value even in the nonlinear viscoelastic regime, which is commensurate with the assumption that the tube diameter remains constant with deformation [31]. According to the model, chains are only oriented, but not stretched. Therefore, the segmental tube orientation S DE is the only contributor to the extra-stress tensor σ, and if expressed as a single integral constitutive equation, it is given by, Therefore, by finding the ratio between the steady-state viscosities of the HMMSF model (η E ( . ε)) and the DE model (η DE ( . ε)) at a Hencky strain well-within the steady-state regime (e.g., ε H = 5 which is also experimentally accessible), we can define a more robust measure of strain hardening since the effects of orientation and stretch are distinctly separated. Hence, our proposed MSHF becomes, As both models are readily available in the IRIS software, calculating the MSHF is an effortless task (see Figure 11d). Figure 10 demonstrates the values of our proposed MSHF (Equation (13)) at different rates. It is obvious that the MSHF is rather rate-insensitive for all three rates, and it becomes rate-independent for the 2 highest elongation rates investigated ( . ε = 0.5, 1/s). The observed rate-insensitivity/independence indicates that the stretch part of the elongational flow becomes saturated in the experimental window (at ε H = 5). To examine this observation, one must investigate the spectral dependence of the stretch (see figure 3c of [25]). To further investigate the stretch saturation observed in Figure 10, the relaxation spectral dependence of the stretch f (obtained from evolution equation Equation (6)) for PE124-10 at . ε = 0.1, 0.5, 1/s is plotted in Figure 11. The spectral dependence demonstrates the importance of the modes with longest relaxation times τ i , which feature the high dynamic dilution and, subsequently, lower mass fraction w i after dilution; hence, these modes have a large effect on the prediction of the elongational viscosity [25].
It is evident (Figure 11a-c) that the stretch is fully saturated at Hencky strain ε H = 5 for all strain rates. Moreover, it explains the minor difference observed between the values measured by our proposed MSHF Equation (13) at . ε = 0.1 , 0.5, 1/s. As shown in Figure 11a-c, only the four largest modes show f i > 1 at . ε = 0.1 /s while f i > 1 is seen in five largest modes at . ε = 0.5, 1/s. Figure 11d compares the predictions of the DE and HMMSF models as functions of Hencky strain for PE 124-10 obtained from the IRIS software. The inverse relation between the steady-state viscosity and strain rate is evident by both HMMSF and DE models.
Overall, it can be concluded that our proposed MSHF method delivers rate-insensitive and LCB-dependent results which are in agreement with the sheer effect of branching on the extent of the strain hardening behavior.   Figure 12 illustrates the mass fraction of dynamically diluted chains w 2 i (Equation (5) (see also figure 3 of [30]) as a function of relaxation modes and PDI of the samples. Results show that (except for PE 236-12) the polydispersity of the HDPE samples has a direct relation with their w 2 i ; i.e., an increase in PDI leads to the growth of the mass fraction. In other words, polydispersity enhances the dynamic dilution through the broadening of the relaxation spectrum. The same trend was observed for disengagement time (or LCBI) as a function of PDI (Figure 6), and G D as a function of PDI ( Figure 6).

Conclusions
In this study, we showed the successful application of the HMMSF model on the elongational data of four sets of high-density polyethylene (PE 191-8, PE 127-9, PE 124-10 and PE 236-12) with traces of LCB. The HMMSF model, with only one nonlinear material parameter in uniaxial and multiaxial extensional flows, namely the dilution modulus G D , shows excellent prediction of rheological properties of the tested polymers.
We propose a new method for investigating the strength of the strain hardening behavior in polymers based on molecular modelling by determining the ratio between the steady-state viscosities of the HMMSF model (η E ( . ε)) and the DE model (η DE ( . ε)) at a Hencky strain well within the steady-state regime.
The polydispersity of the samples has a direct relation with the mass fraction of dynamically diluted chains, w 2 i , i.e., an increase in PDI leads to the growth of the mass fraction. In other words, polydispersity enhances the dynamic dilution through broadening the relaxation spectrum.

Conflicts of Interest:
The authors declare no conflict of interest.