Predicting High-Density Polyethylene Melt Rheology Using a Multimode Tube Model Derived Using Non-equilibrium Thermodynamics

1. Introduction

As reported by the Society of Plastics Industries (SPI) in 2000, the plastic industry in the US is positioned, in terms of shipment, in the fourth place among manufacturing industries after motor vehicles and equipment, electronic components and accessories, and petroleum refining [1]. A more recent survey estimates the global plastic packaging market to worth $269.6 billion by 2025 with a 3.9% compound annual growth rate (CAGR).1 This alone highlights the impact of plastic materials in our lives and, thus, the significance of optimizing the polymer processing technology. Future polymer processing will focus not on the machine, but on the product [1]. Several instabilities appear in the polymer industry, which makes life difficult for polymer engineers. For example, under certain circumstances, when a molten plastic is forced through a die then the shark-skin defect appears [2]. To avoid this, it was suggested to slow down the manufacturing rate; however, this decreases the production rates of commercial products leading to an increase in cost. Wang et al. have suggested that this defect may be related to a molecular instability corresponding to an oscillation of the absorbed chains in the die exit area between coiled and stretched states [3]. As such, it seems that the answer needed should be sought by keeping the molecular level of description and performing molecular dynamics simulations. The ultimate goal is to predict the properties of a product via numerical simulations based on first molecular principles and multiple-scale techniques [1].

Due to computational limitations, however, this aim was unachievable until the last few years, when the extended evolution of simulation algorithms, the parallelization of these algorithms and the race, very recently undertaken, to construct accurate coarse-grained potentials (derived directly from the atomistic simulations) have revolutionized the field. For example, by topologically and dynamically mapping atomistic simulation results onto the tube notion of de Gennes-Edwards, we have recently been able [4], [5] to obtain the most fundamental function of the tube (reptation) theory (according to which the polymer motion due to entanglements is confined within a tube-like region whose axis coincides with the primitive path of the chain and its diameter provides a measure of the strength of the topological interactions), namely, the segment survival probability function, compare the atomistic simulations results against the predictions of modern tube models [6] and even propose modifications to improve these models on a molecular level [7], [8].

Accurate continuum simulations (usually using a finite element scheme [9]) require the use of accurate constitutive models, able to provide the necessary molecular physics associated with the rheological behavior of polymeric systems. As such, the use of empirical models or without reference to molecular physics may fail to represent even qualitative features of the material behavior. Furthermore, a rather small set of parameters should be included in said models, to which a physical significance must be assigned, and the models should have the capacity to fit simultaneously all given data with a single set of parameter values [10]. Only then would polymer engineers be able to correctly predict the rheological response in industrial processes, and solved several long-standing problems that the industry faces.

However, polymers exhibit a wide spectrum of relaxation times, which gives polymeric fluids a partial memory [11]. Conformation tensor-based models that include only a single mode cannot describe small-amplitude oscillatory shear (SAOS) where a spectrum of relaxation times is needed. Even for dilute solutions, both theory and experiments suggest that a superposition of several exponential modes is obtained [12]. Over the past two decades or so, several researchers employed multiple modes of well-known models in order to improve their predictive capabilities. For example, the Kaye-Bernstein–Kearsley–Zapas (K-BKZ) integral model [13], [14], the Phan-Thien and Tanner (PTT) model [14,15,16,17], and the Giesekus model [14], [17] have been used to predict the rheological behaviour of industrial polymers, such as low-density polyethylene (LDPE) and high-density polyethylene (HDPE). Although such well-known rheological models are able to reproduce experimentally observable features of the material functions in various flows, they fail, as mentioned above, to capture the correct physics.

Polymers with large molecular weights should be described via the use of tube theory mentioned above, which introduces terms such as reptation, chain contour length fluctuations and constraint release (CR) due to the motion of surrounding chains [18]. Under flow, as polymer chains are oriented, a number of entanglements are expected to be on average lost as dictated by the convective constraint release (CCR) mechanism [19], [20] and shown to be the case by detailed atomistic non-equilibrium molecular dynamics (NEMD) simulations [21]. More recently, tube models [22,23,24,25] have been used to predict the appearance of a transient stress undershoot (following the overshoot) at high shear rates in start-up shear, that originated from the molecular tumbling of polymer chains in simple shear. Tube models have also been generalized to account for branches, such as the pom-pom model [26] and, its thermodynamically-admissible version, the Pom-pon [27] model. Also, several works employed multiple-mode versions of well-known tube models to predict the rheological response of industrial polymer systems; we mention here only a few of these works. Inkson et al. [28] used a multimode version of the pom-pom model and found that it can address quantitatively the rheology of LDPE for shear, uniaxial elongation and planar extension. Soulages et al. [29] investigated the lubricated flow of a LDPE in a cross-slot geometry and compared the predictions of the extended Pom-Pom model [30,14] and the modified extended Pom-Pom model [16] versus a plethora of rheological data: in shear they compared against transient and steady-state shear viscosity and first normal stress coefficient, and the steady-state second normal stress difference, and in uniaxial extension they compared against the transient uniaxial extensional viscosity. They noted that both models perform equally well (note that the thermodynamic admissibility of these two models is shown in Ref. [31]). Hoyle et al. [32] evaluated the performance of the multimode pom-pom model against both LDPE and branched HDPE melts, whereas more recently Konaganti et al. [14] employed the double convected pom-pom [26] model to predict the rheological behavior of a high-molecular-weight HDPE melt. Since multimode versions have been illustrated to be superior to single-mode ones, our aim in this paper is to generalize the tube model of Stephanou et al. [22] to its multiple-mode version and use it to predict the rheological response of a HDPE melt.

The structure of the paper is as follows: in Section 2 the new model is briefly derived using NET. In Section 3 we derive the expressions for the relevant rheological material functions obtained by analyzing the asymptotic behavior of the model in the limits of small shear rates. The results obtained with the new model are then presented in Section 4: we first discuss its parameterization and then show how accurately and reliably it can describe the viscoelasticity of HDPE polymer melts. The paper concludes with Section 5 where the most important aspects of our work are summarized, and future plans are highlighted and discussed.

2. The constitutive model

2.1. State variables

This work considers an isothermal and incompressible flow, meaning that the total mass density, ρ, and the entropy density (or temperature) are excluded from the vector of state variables. To characterize the polymer chains, the entanglement strand conformation tensor c, following Stephanou et al. [22], [33], is considered which is made dimensionless through $\tilde{c}=Kc/k_{B}T$, where K denotes the spring constant of the Hookean dumbbells representing the entanglement strands, k_B the Boltzmann constant, and T the absolute temperature [34]. The conformation tensor $\tilde{c}$ refers to one entanglement strand and at equilibrium (zero flow field applied) coincides with the unit tensor. To characterize the multiple modes of the polymer chains, N conformation tensors are considered, one for each mode [34]. Finally, the momentum density M as the hydrodynamic variable is further considered, so that overall, the vector of state variables is expressed as$x=\{M,c^{(1)},c^{(2)},\dots ,c^{(i)},\dots ,c^{(N)}\}$.

2.2. System Hamiltonian

The mechanical part of the system’s Hamiltonian is given by [34]

$$H_m\left(x\right)=K_{en}\left(x\right)+A\left(x\right),$$

(1)

where

$$K_{en}(x)={\displaystyle \int \frac{M^2}{2\rho }}dV,$$

(2)

represents the kinetic energy of the system, whereas [22], [34], [35]:

$$A(x)={\displaystyle \int a(x)}dV=\frac{1}{2}{\displaystyle \sum _{i=1}^N{\displaystyle \int G_e^{(i)}\left[\Phi \left(\mathrm{tr}\left({\tilde{c}}^{(i)}-I\right)\right)-\ln \det {\tilde{c}}^{(i)}\right]}dV},$$

(3)

represents the system’s Helmholtz free energy (with $a(x)$ the Helmholtz free energy density) that includes two contributions: the dimensionless potential $\Phi \left(\mathrm{tr}\left({\tilde{c}}^{(i)}-I\right)\right)$ accounting for chain stretching, and an entropic contribution involving the logarithm of the determinant of the conformation tensor of each mode. Here, $G_e^{(i)}$ is the entanglement modulus of the i^th mode. The partial derivative of the potential with respect to the trace of the conformation tensor defines the (dimensionless) effective spring constant [35], [36] for the i^th mode

$$h\left(\mathrm{tr}{\tilde{c}}^{(i)}\right)\equiv \frac{\partial \Phi \left(\mathrm{tr}{\tilde{c}}^{(i)}\right)}{\partial \mathrm{tr}{\tilde{c}}^{(i)}},$$

(4a)

so that the corresponding Volterra derivative of the free energy with respect to the conformation tensor is

$$\frac{\delta A}{\delta {\tilde{c}}^{(i)}}=\frac{G_e^{(i)}}{2}\left[h\left(\mathrm{tr}{\tilde{c}}^{(i)}\right)I-{\left({\tilde{c}}^{(i)}\right)}^{-1}\right],$$

(4b)

where I is the unit tensor. Here, the FENE-P(Wagner) expression is used:

$$h\left(\mathrm{tr}{\tilde{c}}^{(i)}\right)=\frac{b_e^{(i)}-3}{b_e^{(i)}-\mathrm{tr}{\tilde{c}}^{(i)}}.$$

(4c)

where $b_e^{(i)}$ is the finite extensibility (FENE) parameter of the entanglement strand associated with the i^th mode. As shown by Stephanou et al. [33], $b_e=3\left[{\left(0.82\right)}^2/C_{\infty }\right]\left(M_e/M_0\right)$ (when all FENE parameters are considered equal), where $C_{\infty }$ is the polymer characteristic ratio at infinite chain length, M_e the entanglement molecular weight, and M₀ the average molar mass of a monomer,. For example, for PS melts b_e=54 [33].

2.3. The Poisson and dissipation brackets

Following Beris and Edwards [34], the Poisson bracket for multiple modes is given as

$$\begin{array}{ll}{\left\{F,G\right\}}_c=& -{\displaystyle \sum _{i=1}^N{\displaystyle \int \left[\frac{\delta F}{\delta c_{a\beta }^{(i)}}{\nabla }_{\gamma }\left(c_{a\beta }^{(i)}\frac{\delta G}{\delta M_{\gamma }}\right)-\frac{\delta G}{\delta c_{a\beta }^{(i)}}{\nabla }_{\gamma }\left(c_{a\beta }^{(i)}\frac{\delta F}{\delta M_{\gamma }}\right)\right]dV}}\\ & +{\displaystyle \sum _{i=1}^N{\displaystyle \int c_{\gamma a}^{(i)}\left[\frac{\delta F}{\delta c_{a\beta }^{(i)}}{\nabla }_{\gamma }\left(\frac{\delta G}{\delta M_{\beta }}\right)-\frac{\delta G}{\delta c_{a\beta }^{(i)}}{\nabla }_{\gamma }\left(\frac{\delta F}{\delta M_{\beta }}\right)\right]dV}}\\ & +{\displaystyle \sum _{i=1}^N{\displaystyle \int c_{\gamma \beta }^{(i)}\left[\frac{\delta F}{\delta c_{a\beta }^{(i)}}{\nabla }_{\gamma }\left(\frac{\delta G}{\delta M_{\alpha }}\right)-\frac{\delta G}{\delta c_{a\beta }^{(i)}}{\nabla }_{\gamma }\left(\frac{\delta F}{\delta M_{\alpha }}\right)\right]dV}}\end{array}.$$

(5a)

Note that here, and throughout this paper, Einstein’s summation convention for repeated Greek indices is employed. The complete Poisson bracket is then simply given as:

$$\begin{array}{ll}\left\{F,G\right\}=& -{\displaystyle \int \left[\frac{\delta F}{\delta M_{\gamma }}{\nabla }_{\beta }\left(M_{\gamma }\frac{\delta G}{\delta M_{\beta }}\right)-\frac{\delta G}{\delta M_{\gamma }}{\nabla }_{\beta }\left(M_{\gamma }\frac{\delta F}{\delta M_{\beta }}\right)\right]dV}\\ & +{\left\{F,G\right\}}_c\end{array}.$$

(5b)

Next, the following expression for the dissipation bracket associated with the conformation tensors is used [34]

$$\begin{array}{l}{\left[F,G\right]}_{\mathrm{nec}}=-{\displaystyle \sum _{i=1}^N{\displaystyle \int \frac{\delta F}{\delta c_{a\beta }^{(i)}}{\Lambda }_{\alpha \beta \gamma \epsilon }^{(ii)}\frac{\delta G}{\delta c_{\gamma \epsilon }^{(i)}}dV}}\\ +{\displaystyle \sum _{i=1}^N{\displaystyle \int L_{\alpha \beta \gamma \epsilon }^{(i)}\left[\frac{\delta F}{\delta c_{\gamma \epsilon }^{(i)}}{\nabla }_{\alpha }\left(\frac{\delta G}{\delta M_{\beta }}\right)-\frac{\delta G}{\delta c_{\gamma \epsilon }^{(i)}}{\nabla }_{\alpha }\left(\frac{\delta F}{\delta M_{\beta }}\right)\right]dV}}\end{array}.$$

(6)

The first integral on the right-hand side of Eq. (6) accounts for relaxation effects of each conformation tensor, which is proportional to a fourth–rank relaxation tensor, whereas the second integral allows for non-affine deformation for each conformation tensor. Note that the subscript “nec”, meaning “no entropy production correction,” is added to the dissipation bracket to indicate that this bracket is without terms involving Volterra derivatives with respect to the entropy density, which can be omitted when one considers (as we do here) isothermal systems [34]. Note further that although, in general, the dissipation bracket allows explicit coupling between cross modes, see Eq. (8.2-25) of Beris and Edwards,[34] here are omitted for simplicity. Then,

$${\tilde{\dot{c}}}_{\alpha \beta ,[1]}^{(i)}=-{\Lambda }_{\alpha \beta \gamma \epsilon }^{(ii)}\frac{\delta A}{\delta c_{\gamma \epsilon }^{(i)}}+L_{\alpha \beta \gamma \epsilon }^{(i)}{\nabla }_{\gamma }{u}_{\epsilon },$$

(7a)

where we have defined the upper-convected time derivatives:

$${\tilde{\dot{c}}}_{\alpha \beta ,[1]}^{(i)}\equiv \frac{\partial {\tilde{c}}_{\alpha \beta }^{(i)}}{\partial t}+u_{\gamma }{\nabla }_{\gamma }{c}_{\alpha \beta }^{(i)}-{\tilde{c}}_{\alpha \gamma }^{(i)}{\nabla }_{\gamma }{u}_{\beta }-{\tilde{c}}_{\gamma \beta }^{(i)}{\nabla }_{\gamma }{u}_{\alpha }.$$

(7b)

Finally, the stress tensor is given as,

$${\sigma }_{\alpha \beta }={\displaystyle \sum _{i=1}^N\left(2c_{a\gamma }^{(i)}\frac{\delta A}{\delta c_{\gamma \beta }^{(i)}}+L_{\alpha \beta \gamma \epsilon }^{\left(i\right)}\frac{\delta A}{\delta c_{\gamma \epsilon }^{(i)}}\right)}.$$

(8)

2.4. The Matrices L and Λ

The relaxation matrix for each mode ${\Lambda }_{\alpha \beta \gamma \epsilon }^{(ii)}$ is split into two contributions, following Stephanou et al. [22], and with different relaxation times:

$$\begin{array}{l}{\Lambda }_{\alpha \beta \gamma \epsilon }^{\mathrm{rept},(ii)}=\frac{f_{\mathrm{rept}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)}{2G_e^{(i)}{\tau }_{\mathrm{CR}}^{(i)}}\left({\tilde{c}}_{\alpha \gamma }^{(i)}{\tilde{\beta }}_{\beta \epsilon }^{(i)}+{\tilde{c}}_{\alpha \epsilon }^{(i)}{\tilde{\beta }}_{\beta \gamma }^{(i)}+{\tilde{c}}_{\beta \gamma }^{(i)}{\tilde{\beta }}_{\alpha \epsilon }^{(i)}+{\tilde{c}}_{\beta \epsilon }^{(i)}{\tilde{\beta }}_{\alpha \gamma }^{(i)}\right)\\ {\Lambda }_{\alpha \beta \gamma \epsilon }^{\mathrm{Rouse},(ii)}=\frac{f_{\mathrm{Rouse}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)}{2G_e^{(i)}{\tau }_R^{(i)}\left(tr\tilde{c}\right)}\left({\tilde{c}}_{\alpha \gamma }^{(i)}{\tilde{\beta }}_{\beta \epsilon }^{(i)}+{\tilde{c}}_{\alpha \epsilon }^{(i)}{\tilde{\beta }}_{\beta \gamma }^{(i)}+{\tilde{c}}_{\beta \gamma }^{(i)}{\tilde{\beta }}_{\alpha \epsilon }^{(i)}+{\tilde{c}}_{\beta \epsilon }^{(i)}{\tilde{\beta }}_{\alpha \gamma }^{(i)}\right)\end{array}.$$

(9a)

Here, ${\tau }_{\mathrm{CR}}^{(i)}$ is the CR relaxation time of the i^th mode, which is considered to be half the corresponding reptation time, ${\tau }_{\mathrm{CR}}^{(i)}=\frac{1}{2}{\tau }_d^{(i)}$[37] [note that this coincides with the equilibrium CCR relaxation time, see Eq. (9e)], ${\tau }_R^{(i)}\left(tr\tilde{c}\right)$ is the Rouse relaxation time of the i^th mode,

$${\tau }_R^{(i)}\left(tr\tilde{c}\right)={\tau }_{\mathrm{R},\mathrm{eq}}^{(i)}{\left(\frac{tr{\tilde{c}}^{(i)}}{3}\right)}^{k^{(i)}},$$

(9b)

where ${\tau }_{R,\mathrm{eq}}^{(i)}$ is the equilibrium Rouse relaxation time of the i^th mode given as ${\tau }_d^{(i)}=3Z{\tau }_{R,\mathrm{eq}}^{(i)}$ [18], and $k^{(i)}$ is the Extended-White Metzner (EWM) exponent [34] for the i^th mode. Note that for the Rouse time a shear-rate-dependency through the use of the trace of the conformation tensor of each mode is considered. The functions$f_{\mathrm{rep}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)$ and $f_{\mathrm{Rouse}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)$ are scalar functions of the trace of the conformation tensor defined via:[22]

$$f_{\mathrm{Rouse}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)=1-f_{\mathrm{rep}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)={\beta }_{ccr}^{(i)}\frac{h\left(tr{\tilde{c}}^{(i)}\right)tr{\tilde{c}}^{(i)}-3}{3+{\beta }_{ccr}^{(i)}\left[h\left(tr{\tilde{c}}^{(i)}\right)tr{\tilde{c}}^{(i)}-3\right]},$$

(9c)

where ${\beta }_{ccr}^{(i)}$ is the CCR parameter of the i^th mode. For the (dimensionless) mobility tensor ${\tilde{\beta }}^{(i)}$ of the i^th mode, the Giesekus’ postulate ${\tilde{\beta }}^{(i)}=I+{\alpha }^{(i)}{\tilde{\sigma }}^{(i)}$ is used [35] with ${\tilde{\sigma }}^{(i)}={\sigma }^{(i)}/G_i$ and ${\alpha }^{(i)}$ is the anisotropic mobility (Giesekus) parameter [14], [17] of the i^th mode. Then, with ${\Lambda }_{\alpha \beta \gamma \epsilon }^{(ii)}={\Lambda }_{\alpha \beta \gamma \epsilon }^{\mathrm{rep},(ii)}+{\Lambda }_{\alpha \beta \gamma \epsilon }^{\mathrm{Rouse},(ii)}$ we obtain

$${\Lambda }_{\alpha \beta \gamma \epsilon }^{\left(ii\right)}=\frac{1}{2G_e^{(i)}{\tau }_{\mathrm{CCR}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)}\left({\tilde{c}}_{\alpha \gamma }^{(i)}{\tilde{\beta }}_{\beta \epsilon }^{(i)}+{\tilde{c}}_{\alpha \epsilon }^{(i)}{\tilde{\beta }}_{\beta \gamma }^{(i)}+{\tilde{c}}_{\beta \gamma }^{(i)}{\tilde{\beta }}_{\alpha \epsilon }^{(i)}+{\tilde{c}}_{\beta \epsilon }^{(i)}{\tilde{\beta }}_{\alpha \gamma }^{(i)}\right),$$

(9d)

where [22], [37] the CCR relaxation time is obtained as

$$\frac{1}{{\tau }_{\mathrm{CCR}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)}=\frac{1}{{\tau }_{\mathrm{CR}}^{(i)}}+\left(\frac{1}{{\tau }_R^{(i)}\left(tr{\tilde{c}}^{(i)}\right)}-\frac{1}{{\tau }_{\mathrm{CR}}^{(i)}}\right){\beta }_{ccr}^{(i)}\frac{h\left(tr{\tilde{c}}^{(i)}\right)tr{\tilde{c}}^{(i)}-3}{3+{\beta }_{ccr}^{(i)}\left[h\left(tr{\tilde{c}}^{(i)}\right)tr{\tilde{c}}^{(i)}-3\right]}.$$

(9e)

Finally, the expression for the $L_{}^{\left(i\right)}$ matrix is given via [34], [35]

$$L_{\alpha \beta \gamma \epsilon }^{\left(i\right)}=-\frac{{\xi }^{(i)}}{2}\left({\tilde{c}}_{\alpha \gamma }^{(i)}{\delta }_{\beta \epsilon }+{\tilde{c}}_{\alpha \epsilon }^{(i)}{\delta }_{\beta \gamma }+{\tilde{c}}_{\beta \gamma }^{(i)}{\delta }_{\alpha \epsilon }+{\tilde{c}}_{\beta \epsilon }^{(i)}{\delta }_{\alpha \gamma }\right).$$

(11)

Here, ${\xi }^{(i)}$ is the non-affine/slip parameter of the i^th mode. This parameter is important as it allows for the prediction of a transient stress undershoot (following the overshoot) at high shear rates [22].

2.5. Thermodynamic admissibility

Any thermodynamic system must obey the restriction of a non-negative rate of total entropy production. When the fluid studied is isothermal and incompressible, the entropy production results from the degradation of mechanical energy leading to $d{H}_m/dt=\left[H_m,H_m\right]\le 0$.[34] For this to be satisfied in our model the following must hold:

$${\displaystyle \sum _{i=1}^N\frac{\delta F}{\delta c_{a\beta }^{(i)}}{\Lambda }_{\alpha \beta \gamma \epsilon }^{(i\iota )}\frac{\delta G}{\delta c_{\gamma \epsilon }^{(i)}}}={\displaystyle \sum _{i=1}^N\frac{G_e^{(i)}}{2{\tau }_{\mathrm{CCR}}^{(i)}\left(tr{\tilde{c}}^{(i)}\right)}{\displaystyle \sum _{k=1}^3\frac{{\left(h{\mu }_k^{(i)}-1\right)}^2}{{\mu }_k^{(i)}}\left[1+{\alpha }^{(i)}\left(1-{\xi }^{(i)}\right)\left(h{\mu }_k^{(i)}-1\right)\right]}}\ge 0,$$

(12)

where ${\mu }_k^{(i)},k=\{1,2,3\},i=\{1,\mathrm{..},N\}$ are the three eigenvalues of the conformation tensor of the i^th mode. Obviously, since the condition$0\le {\alpha }^{(i)}\left(1-{\xi }^{(i)}\right)<1,0\le {\xi }^{(i)}<1,\forall i$and ${\beta }_{ccr}^{(i)}\ge 0,\forall i$ [22] guarantees that each term of the summation is positive, then the sum as a whole is positive as well, meaning that the multi-mode version of the model is thermodynamically admissible.

2.6. Conformation tensor evolution equation

The evolution equation for each of the dimensionless conformation tensors is obtained by substituting Eqs. (4b), (9d) and (11) in Eq. (7a),

$$\begin{array}{l}{\tilde{\dot{c}}}_{[1]}^{(i)}=-\frac{1}{{\tau }_{\mathrm{CCR}}^{(i)}\left(\mathrm{tr}{\tilde{c}}^{(i)}\right)}\{{\alpha }^{(i)}\left(1-{\xi }^{(i)}\right)h^2\left(\mathrm{tr}{\tilde{c}}^{(i)}\right){\tilde{c}}^{(i)}\cdot {\tilde{c}}^{(i)}\\ +\left[1-2{\alpha }^{(i)}\left(1-{\xi }^{(i)}\right)\right]h\left(\mathrm{tr}{\tilde{c}}^{(i)}\right){\tilde{c}}^{(i)}-\left[1-{\alpha }^{(i)}\left(1-{\xi }^{(i)}\right)\right]I\},\forall i\in [1,N]\end{array},$$

(13)

The CCR relaxation time for the i^th mode is given by eq. (9d) and the (dimensionless) effective spring constant is given by Eq. (4c). Finally, the expression for the extra (polymeric) stress tensor is obtained by substituting Eqs. (4b) and (11) in Eq. (8) as

$$\sigma ={\displaystyle \sum _{i=1}^N{G}_e^{(i)}\left[h\left(\mathrm{tr}{\tilde{c}}^{(i)}\right){\tilde{c}}^{(i)}-I\right]}.$$

(14)

3. Asymptotic behavior of the model in steady state shear

In this section, we provide analytical expressions describing the asymptotic behavior of the multiple-mode version of the Stephanou et al. [22] model in the limit of weak flows for the following three cases: inception of simple shear flow (SSF) described by the kinematics $\mathbf{u}=\left(\dot{\gamma }y,0,0\right)$ where $\dot{\gamma }$ is the (constant) shear rate, inception of uniaxial elongation flow (UEF) described by the kinematics $\mathbf{u}=\left(\dot{\epsilon }x,-\frac{1}{2}\dot{\epsilon }y,-\frac{1}{2}\dot{\epsilon }z\right)$ where $\dot{\epsilon }$ is the (constant) elongation rate, and small amplitude oscillatory shear (SAOS) described by the kinematics $u=\left(\dot{\gamma }\cos \left(\omega t\right)y,0,0\right)$ where ω is the oscillation frequency. The material functions to be analyzed are: a) the transient shear viscosity ${\eta }^{+}(t)$ (= ${\sigma }_{\mathrm{yx}}\left(t\right)/\dot{\gamma }$) and the first, ${\Psi }_1^{+}(t)$, (=$\left({\sigma }_{\mathrm{xx}}(t)-{\sigma }_{\mathrm{yy}}(t)\right)/{\dot{\gamma }}^2$) and second, ${\Psi }_2^{+}(t)$, (=$\left({\sigma }_{\mathrm{yy}}(t)-{\sigma }_{\mathrm{zz}}(t)\right)/{\dot{\gamma }}^2$) normal stress coefficients in the case of shear, b) the transient elongational viscosity ${\eta }_E^{+}(t)$, (=$\left({\sigma }_{\mathrm{xx}}(t)-{\sigma }_{\mathrm{yy}}(t)\right)/\dot{\epsilon }$), in the case of uniaxial elongation, and c) the storage, $G^{\prime }\left(\omega \right)$, and loss, $G^{″}\left(\omega \right)$, moduli in the case of SAOS.

To obtain the asymptotic behavior, we need to expand the conformation tensor for each mode in the limit of small strain rates (by invoking a linearization of the evolution equation for each of the conformation tensors) and solve the corresponding ordinary differential equations analytically; after this, we obtain the non-zero stress tensor components via Eq. (14). Finally, we obtain the following results for the relevant material functions:

Inception of shear:

$${\eta }^{+}\left(t\right)={\displaystyle \sum _{i=1}^N{\tau }_{\mathrm{CR}}^{(i)}{\tilde{G}}_e^{(i)}\left[1-\exp \left(-\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}\right)\right]},$$

(15a)

$${\Psi }_1^{+}\left(t\right)=2{\displaystyle \sum _{i=1}^N{\left({\tau }_{\mathrm{CR}}^{(i)}\right)}^2{\tilde{G}}_e^{(i)}\left[1-\left(1+\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}\right)\exp \left(-\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}\right)\right]},$$

(15b)

$$\begin{array}{l}-{\Psi }_2^{+}\left(t\right)={\displaystyle \sum _{i=1}^N{\left({\tau }_{\mathrm{CR}}^{(i)}\right)}^2{\tilde{G}}_e^{(i)}\{\left[{\xi }^{(i)}+{\alpha }^{(i)}\left(1-{\xi }^{(i)}\right)\right]\left[1-\left(1+\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}\right)\exp \left(-\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}\right)\right]}\\ {\alpha }^{(i)}{\left(1-{\xi }^{(i)}\right)}^2\exp \left(-\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}\right)\left[1-\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}-\exp \left(-\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}\right)\right]\}\end{array}.$$

(15c)

Inception of uniaxial elongation:

$${\eta }_E^{+}\left(t\right)=3{\displaystyle \sum _{i=1}^N{\tau }_{\mathrm{CR}}^{(i)}{\tilde{G}}_e^{(i)}\left[1-\exp \left(-\frac{t}{{\tau }_{\mathrm{CR}}^{(i)}}\right)\right]}.$$

(16)

meaning that Trouton’s law holds for the steady-state extensional viscosity. Small Amplitude Oscillatory Shear:

$$\begin{array}{l}{G}^{\prime }\left(\omega \right)={\displaystyle \sum _{i=1}^N{\tilde{G}}_e^{(i)}\frac{\left(\omega {\tau }_{\mathrm{CR}}^{(i)}\right)}{1+{\left(\omega {\tau }_{\mathrm{CR}}^{(i)}\right)}^2}.}\\ G^{″}\left(\omega \right)={\displaystyle \sum _{i=1}^N{\tilde{G}}_e^{(i)}\frac{\omega {\tau }_{\mathrm{CR}}^{(i)}}{1+{\left(\omega {\tau }_{\mathrm{CR}}^{(i)}\right)}^2}}\end{array}$$

(17)

In equations (15)-(17), we have defined ${\tilde{G}}_e^{(i)}\equiv {\left(1-{\xi }^{(i)}\right)}^2{G}_e^{(i)}$.

4. Results and Discussion

The FENE parameter can be easily calculated via $b_e=3{\left(0.82\right)}^2{M}_e/\left(M_0{C}_{\infty }\right)$ as mentioned above. For PE M₀=14 gr/mole, whereas C_∞=7.3 and M_e = 828 g/mole (see Table 3.3, p. 151 of Ref. [38]). These values yield b_e = 16.34. We will compare against the experimental data of Konaganti et al. [14] that have performed rheological measurements of the sample HDPE-1 reported by the same group [39] for which M_w=206 kg/mole; thus, the number of entanglements is equal to Z≈249. The relaxation spectrum is the same as the one used by Konaganti et al. [14] (see their Table 2 for T=200 ^oC) but is also provided in Table 1, and it was obtained by fitting the storage and loss moduli, Eqs. (7), with the corresponding experimental data. The comparison against the experimental storage and loss moduli is available in Figure 1.

All the remaining parameter values are obtained by fitting the model predictions against the experimental data; note that for simplicity we assume that each parameter has the same value for all modes although, in general, different values for each mode could be considered (e.g., see the work of Konaganti et al. [14]). The following values of the model parameter are chosen to provide a good comparison against the experimental data: ξ=0.02, α=0.3, β_ccr=4×10⁻⁴, and k=-3.5.

4.1. Comparison with start-up shear flow data

Figure 2 illustrates the comparison between the experimental data for the growth of the shear viscosity upon inception of shear flow and the simulated results obtained from the model. The experimental data (blue symbols) were collected at three different shear rates: 0.05 1/s, 0.5 1/s, and 1 1/s, whereas the orange lines represent the simulated shear viscosity values at the corresponding shear rates. It can be observed that the model accurately captures the trends and magnitude of the shear viscosity over time. Note that although a non-zero value of ξ is employed, the undershoots produced are too small to note in the scale used for Figure 2. Although no experimental data are provided, we further provide here the corresponding prediction of the growth of the first and second normal stress coefficients in Figure 3, and the steady-state values of all viscometric material functions in Figure 4.

4.2. Comparison with start-up uniaxial elongational flow data

In Figure 5 we provide the comparison of the model predictions with the experimental measurements of Konaganti et al. [14] for the growth of the elongational viscosity as a function of time upon inception of uniaxial elongational at three different stretch rates: 0.05 1/s, 0.5 1/s, and 5 1/s. The comparison employs the same parameter values as the ones used in Figure 2 except ξ=0, as instructed by Stephanou et al. [35], and the corresponding steady-state prediction is provided in Figure 6. Given that the parameter values were selected based on the start-up shear data (Figure 2), the comparison is noted to be adequately well. Overall, the proposed model demonstrates good agreement with the experimental data, indicating its effectiveness in describing the rheological behavior of polymers of industrial significance.

5. Conclusions

In this work, we have developed the multiple-mode version of a generalized, conformation-tensor based viscoelastic model for polymer melts, that has been proposed recently by one of us [22], by making use of the generalized bracket formalism of Beris and Edwards [34]. As its forerunner, the monodisperse version [22], it accounts for the most significant effects that can be realized in entangled polymers systems: anisotropic drag, finite extensibility, non-affine motion (leading to the exhibition of a transient stress undershoot at large shear rates), and variable chain relaxation due to convective constraint release. The multiple-mode version of the model has been shown to bear a very good predictive capability with regards to the industrial experimental data for HDPE of Konaganti et al. [14].

The model in its present form considers only strictly linear chains. Industrial samples, particularly LDPE, are never strictly linear having either short or long branches distributed along their entire backbone. There is clear evidence that all material functions of PE are considerably affected by even the presence of low levels of long chain branching [40]. We therefore need to generalize it and allow for the explicit consideration of branches, following the guidelines provided by the pom-pom [26] and, its thermodynamically-admissible version, the Pom-pon [27] model. Furthermore, the multi-mode version does not explicitly consider the molecular weight (MW) distribution of industrial samples, such as the log-normal of gamma distributions that are able to describe experimental distributions [41]. As such, we should also generalize it to handle molecular weight distribution following the work of Schieber [42], [43]. This generalized constitutive model will allow for the more accurate prediction of the rheological response of polymeric systems used industrially that do possess an extensive spectrum of MW and not a very narrow distribution as it is customarily assumed when deriving a constitutive model. Only then would polymer engineers be able to accurately use the predictions of the revised constitutive model against the rheological response noted in actual industrial processes. Our findings provide a foundation for future research aimed at enhancing the properties of high MW polymers for diverse applications.