Next Article in Journal
Properties of Resorbable Conduits Based on Poly(L-Lactide) Nanofibers and Chitosan Fibers for Peripheral Nerve Regeneration
Next Article in Special Issue
Wall Slip-Free Viscosity Determination of Filled Rubber Compounds Using Steady-State Shear Measurements
Previous Article in Journal
Behaviour of the Peri-Implant Soft Tissue with Different Rehabilitation Materials on Implants
Previous Article in Special Issue
Molecular Processes Leading to Shear Banding in Entangled Polymeric Solutions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Pavlina C. Konstantinou
and
Pavlos S. Stephanou
*
Department of Chemical Engineering, Cyprus University of Technology, P.O. Box 50329, 3603 Limassol, Cyprus
*
Author to whom correspondence should be addressed.
Polymers 2023, 15(15), 3322; https://doi.org/10.3390/polym15153322
Submission received: 18 July 2023 / Revised: 31 July 2023 / Accepted: 3 August 2023 / Published: 7 August 2023
(This article belongs to the Special Issue Rheological Properties of Polymers and Polymer Composites)

Abstract

:
Based on the Generalized bracket, or Beris–Edwards, formalism of non-equilibrium thermodynamics, we recently proposed a new differential constitutive model for the rheological study of entangled polymer melts and solutions. It amended the shortcomings of a previous model that was too strict regarding the values of the convective constraint release parameter for the model not to violate the second law of thermodynamics, and it has been shown capable of predicting a transient stress undershoot (following the overshoot) at high shear rates. In this study, we wish to further examine this model’s capability to predict the rheological response of industrial polymer systems by extending it to its multiple-mode version. The comparison with industrial rheological data (High-Density Polyethylene resins), which was based on comparison with experimental data available in (a) Small Amplitude Oscillatory shear, (b) start-up shear, and (c) start-up uniaxial elongation, was noted to be good.

Graphical Abstract

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 fourth place among manufacturing industries, following motor vehicles and equipment, electronic components and accessories, and petroleum refining [1]. A more recent survey predicts that the global plastic packaging market will be worth $269.6 billion by 2025, achieving a compound annual growth rate (CAGR) [2]. This figure alone highlights the impact of plastic materials on our lives and, thus, the significance of optimizing 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 molten plastic is forced through a die, the shark-skin defect appears [3]. To avoid this issue, it is suggested to slow down the manufacturing rate; however, this action 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 that corresponds to an oscillation of the absorbed chains in the die exit area between coiled and stretched states [4]. Thus, it seems that the answer needed should be sought by maintaining the molecular level of description and performing molecular dynamics simulations. The ultimate goal of this process is to predict the properties of a product via numerical simulations based on molecular principles and multiple-scale techniques [1].
Due to computational limitations, however, this goal has been unachievable until the last few years, in which period the extended evolution of simulation algorithms, the parallelization of these algorithms and the race, which has been very recently undertaken, to construct accurate coarse-grained potentials (directly derived from the atomistic simulations) have revolutionized the field. For example, by topologically and dynamically mapping atomistic simulation results onto the tube notion of the de Gennes–Edwards model, we have recently been able [5,6] 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, the axis of which 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 [7]; and even propose modifications to improve these models on a molecular level [8,9].
Accurate continuum simulations (usually using a finite element scheme [10]) require the use of accurate constitutive models, which are able to provide the necessary molecular physics associated with the rheological behavior of polymeric systems. Thus, the use of empirical models or models without reference to molecular physics may fail to represent even the qualitative features of the material’s 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 simultaneously fit all given data using a single set of parameter values [11]. Only in this context would polymer engineers be able to correctly predict the rheological response in industrial processes and solve several long-standing problems that the industry faces.
However, polymers exhibit a wide spectrum of relaxation times, which gives polymeric fluids a partial memory [12]. Conformation tensor-based models that include only a single mode cannot describe small-amplitude oscillatory shear (SAOS), in which 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 [13]. In 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 [14,15], the Phan-Thien and Tanner (PTT) model [15,16,17,18], and the Giesekus model [15,18] have been used to predict the rheological behaviors 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 to capture the correct physics.
Polymers with large molecular weights should be described based on the use of the 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 [19]. Under flow, as polymer chains are oriented, a number of entanglements are expected, on average, to be lost, as dictated by the convective constraint release (CCR) mechanism [20,21] and shown to be the case via detailed atomistic non-equilibrium molecular dynamics (NEMD) simulations [22]. More recently, tube models [23,24,25,26] have been used to predict the appearance of a transient stress undershoot (following the overshoot) at high shear rates in start-up shear, which 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 [27], as well as its thermodynamically admissible version, known as the pom-pon [28] model. Also, several works employed multimode versions of well-known tube models to predict the rheological responses of industrial polymer systems (we only mention a few of these works). Inkson et al. [29] used a multimode version of the pom-pom model and found that it can quantitatively address the rheology of LDPE for shear, uniaxial, and planar elongation. Soulages et al. [30] investigated the lubricated flow of a LDPE in a cross-slot geometry and compared the predictions of the extended pom-pom model [31] and the modified extended pom-pom model [16] to a plethora of rheological data: in shear, they compared it to transient and steady-state shear viscosity and the first normal stress coefficient, as well as the steady-state second normal stress difference, and in uniaxial elongation, they compared it to the transient uniaxial elongational viscosity. They noted that both models performed equally well (note that the thermodynamic admissibility of these two models is shown in Ref. [32]). Hoyle et al. [33] evaluated the performance of the multimode pom-pom model to those of both LDPE and branched HDPE melts, whereas, more recently, Konaganti et al. [15] employed the double convected pom-pom [34] model to predict the rheological behavior of a high-molecular-weight HDPE melt. Since multimode versions have been found to be superior to single-mode versions, our aim in this paper is to generalize the tube model of Stephanou et al. [23] to its multimode 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 non-equilibrium thermodynamics (NET). In Section 3, we derive the expressions of 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, in which 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, in which 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 the method of Stephanou et al. [23,35], is considered to be made dimensionless through c ˜ = K c / k B T , where K denotes the spring constant of the Hookean dumbbells that represents the entanglement strands at equilibrium, kB the Boltzmann constant, and T denotes the absolute temperature [36]. The conformation tensor c ˜ refers to one entanglement strand, and at equilibrium (zero flow field applied), it coincides with the unit tensor. To characterize the multiple modes of the polymer chains, N conformation tensors are considered, with one tensor being considered for each mode [36]. Finally, the momentum density M, which is the hydrodynamic variable, is further considered, meaning that, overall, the vector of state variables is expressed as x = { M , c ( 1 ) , c ( 2 ) , , c ( i ) , , c ( N ) } .

2.2. System Hamiltonian

The mechanical part of the system’s Hamiltonian is given as [36]
H m ( x ) = K e n ( x ) + A ( x )
where
K e n ( x ) = M 2 2 ρ d V
represents the kinetic energy of the system, whereas [23,35,36]
A ( x ) = a ( x ) d V = 1 2 i = 1 N G e ( i ) [ Φ ( tr ( c ˜ ( i ) I ) ) ln det c ˜ ( i ) ] d V
represents the system’s Helmholtz free energy (with a ( x ) the Helmholtz free energy density) that includes the following contributions: the dimensionless potential Φ ( tr ( c ˜ ( i ) I ) ) , which accounts for chain stretching, and an entropic contribution, which involves the logarithm of the determinant of the conformation tensor of each mode [35,36]. Here, G e ( i ) is the entanglement modulus of the ith mode, and I is the unit tensor. The partial derivative of the potential with respect to the trace of the conformation tensor defines the (dimensionless) effective spring constant [37,38] for the ith mode
h ( tr c ˜ ( i ) ) Φ ( tr c ˜ ( i ) ) tr c ˜ ( i )
meaning that the corresponding Volterra derivative of the free energy with respect to the conformation tensor is
δ A δ c ˜ ( i ) = G e ( i ) 2 [ h ( tr c ˜ ( i ) ) I ( c ˜ ( i ) ) 1 ]
Here, the following FENE-P(Wagner) expression is used:
h ( tr c ˜ ( i ) ) = b e ( i ) 3 b e ( i ) tr c ˜ ( i )
where b e ( i ) is the finite extensibility (FENE) parameter of the entanglement strand associated with the ith mode. As shown by the study of Stephanou et al. [35], b e = 3 [ ( 0.82 ) 2 / C ] ( M e / M 0 ) (when all FENE parameters are considered equal), where C is the polymer characteristic ratio at infinite chain length, Me is the entanglement molecular weight, and M0 is the average molar mass of a monomer. For example, for PS melts, be = 54 [35].

2.3. The Poisson and Dissipation Brackets

Following the work of Beris and Edwards [36], the Poisson bracket for multiple modes is given as follows:
{ F , G } c = i = 1 N [ δ F δ c a β ( i ) γ ( c a β ( i ) δ G δ M γ ) δ G δ c a β ( i ) γ ( c a β ( i ) δ F δ M γ ) ] d V + i = 1 N c γ a ( i ) [ δ F δ c a β ( i ) γ ( δ G δ M β ) δ G δ c a β ( i ) γ ( δ F δ M β ) ] d V + i = 1 N c γ β ( i ) [ δ F δ c a β ( i ) γ ( δ G δ M α ) δ G δ c a β ( i ) γ ( δ F δ M α ) ] d V
We note both here and throughout this paper that Einstein’s summation convention for repeated Greek indices is employed. The complete Poisson bracket was then simply given as follows:
{ F , G } = [ δ F δ M γ β ( M γ δ G δ M β ) δ G δ M γ β ( M γ δ F δ M β ) ] d V + { F , G } c
Next, the following expression for the dissipation bracket associated with the conformation tensors is used [36].
[ F , G ] nec = i = 1 N δ F δ c a β ( i ) Λ α β γ ε ( i i ) δ G δ c γ ε ( i ) d V + i = 1 N L α β γ ε ( i ) [ δ F δ c γ ε ( i ) α ( δ G δ M β ) δ G δ c γ ε ( i ) α ( δ F δ M β ) ] d V
The first integral on the right-hand side of Equation (6) accounts for the relaxation effects of each conformation tensor, which is proportional to a fourth–rank relaxation tensor, whereas the second integral allows for the non-affine deformation of each conformation tensor. We note that the subscript “nec”, meaning “no entropy production correction”, is added to the dissipation bracket to indicate that this bracket lacks terms that involve Volterra derivatives with respect to the entropy density, which can be omitted when we consider isothermal systems [36]. We further note that although, in general, the dissipation bracket allows explicit coupling between cross modes, as shown in Equations (8.2–25) of Beris and Edwards [36], in this study, they are omitted for simplicity. Then,
c ˙ ˜ α β , [ 1 ] ( i ) = Λ α β γ ε ( i i ) δ A δ c γ ε ( i ) + L α β γ ε ( i ) γ u ε
where we have defined the upper-convected time derivative:
c ˙ ˜ α β , [ 1 ] ( i ) c ˜ α β ( i ) t + u γ γ c α β ( i ) c ˜ α γ ( i ) γ u β c ˜ γ β ( i ) γ u α
Finally, the extra (polymeric) stress tensor is given as follows:
σ α β = i = 1 N ( 2 c a γ ( i ) δ A δ c γ β ( i ) + L α β γ ε ( i ) δ A δ c γ ε ( i ) )

2.4. The Matrices L and Λ

The relaxation matrix of each mode Λ α β γ ε ( i i ) is split into two contributions, following the work of Stephanou et al. [23], which have different relaxation times:
Λ α β γ ε rept , ( i i ) = f rept ( i ) ( t r c ˜ ( i ) ) 2 G e ( i ) τ CR ( i ) ( c ˜ α γ ( i ) β ˜ β ε ( i ) + c ˜ α ε ( i ) β ˜ β γ ( i ) + c ˜ β γ ( i ) β ˜ α ε ( i ) + c ˜ β ε ( i ) β ˜ α γ ( i ) ) Λ α β γ ε Rouse , ( i i ) = f Rouse ( i ) ( t r c ˜ ( i ) ) 2 G e ( i ) τ R ( i ) ( t r c ˜ ) ( c ˜ α γ ( i ) β ˜ β ε ( i ) + c ˜ α ε ( i ) β ˜ β γ ( i ) + c ˜ β γ ( i ) β ˜ α ε ( i ) + c ˜ β ε ( i ) β ˜ α γ ( i ) )
Here, τ CR ( i ) is the CR relaxation time of the ith mode, which is considered to be half of the corresponding reptation time, τ CR ( i ) = 1 2 τ d ( i ) [39] [we note that this time coincides with the CCR relaxation time at equilibrium, as shown in Equation (10)], and τ R ( i ) ( t r c ˜ ) is the Rouse relaxation time of the ith mode,
τ R ( i ) ( t r c ˜ ) = τ R , eq ( i ) ( t r c ˜ ( i ) 3 ) k ( i )
where τ R , eq ( i ) is the equilibrium Rouse relaxation time of the ith mode, which is given as τ d ( i ) = 3 Z τ R , eq ( i ) [19], and k ( i ) is the Extended White–Metzner (EWM) exponent [36] for the ith mode. We 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 rep ( i ) ( t r c ˜ ( i ) ) and f Rouse ( i ) ( t r c ˜ ( i ) ) are scalar functions of the trace of the conformation tensor, as defined via the following equation [23]:
f Rouse ( i ) ( t r c ˜ ( i ) ) = 1 f rep ( i ) ( t r c ˜ ( i ) ) = β c c r ( i ) h ( t r c ˜ ( i ) ) t r c ˜ ( i ) 3 3 + β c c r ( i ) [ h ( t r c ˜ ( i ) ) t r c ˜ ( i ) 3 ]
where β c c r ( i ) is the CCR parameter of the ith mode. For the (dimensionless) mobility tensor β ˜ ( i ) of the ith mode, the Giesekus’ postulate β ˜ ( i ) = I + α ( i ) σ ˜ ( i ) is used [37] with σ ˜ ( i ) = σ ( i ) / G i , and α ( i ) is the anisotropic mobility (Giesekus) parameter of the ith mode. Then, with Λ α β γ ε ( i i ) = Λ α β γ ε rep , ( i i ) + Λ α β γ ε Rouse , ( i i ) , we obtain
Λ α β γ ε ( i i ) = 1 2 G e ( i ) τ CCR ( i ) ( t r c ˜ ( i ) ) ( c ˜ α γ ( i ) β ˜ β ε ( i ) + c ˜ α ε ( i ) β ˜ β γ ( i ) + c ˜ β γ ( i ) β ˜ α ε ( i ) + c ˜ β ε ( i ) β ˜ α γ ( i ) )
where [23,39] the CCR relaxation time is obtained as follows:
1 τ CCR ( i ) ( t r c ˜ ( i ) ) = 1 τ CR ( i ) + ( 1 τ R ( i ) ( t r c ˜ ( i ) ) 1 τ CR ( i ) ) β c c r ( i ) h ( t r c ˜ ( i ) ) t r c ˜ ( i ) 3 3 + β c c r ( i ) [ h ( t r c ˜ ( i ) ) t r c ˜ ( i ) 3 ]
Finally, the expression of the L ( i ) matrix is given via the following equation [36,37]:
L α β γ ε ( i ) = ξ ( i ) 2 ( c ˜ α γ ( i ) δ β ε + c ˜ α ε ( i ) δ β γ + c ˜ β γ ( i ) δ α ε + c ˜ β ε ( i ) δ α γ )
Here, ξ ( i ) is the non-affine/slip parameter of the ith mode. This parameter is important, as it allows for the prediction of a transient stress undershoot (following the overshoot) at high shear rates [23].

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 the mechanical energy, leading to d H m / d t = [ H m , H m ] 0 [36]. For this aspect to be satisfied in our model, the following equation must hold:
i = 1 N δ F δ c a β ( i ) Λ α β γ ε ( i ι ) δ G δ c γ ε ( i ) = i = 1 N G e ( i ) 2 τ CCR ( i ) ( t r c ˜ ( i ) ) k = 1 3 ( h μ k ( i ) 1 ) 2 μ k ( i ) [ 1 + α ( i ) ( 1 ξ ( i ) ) ( h μ k ( i ) 1 ) ] 0
where μ k ( i ) , k = { 1 , 2 , 3 } , i = { 1 , .. , N } are the three eigenvalues of the conformation tensor of the ith mode. Obviously, since the conditions 0 α ( i ) ( 1 ξ ( i ) ) < 1 , 0 ξ ( i ) < 1 , i and β c c r ( i ) 0 , i [23] guarantee that each term of the summation is positive, the sum as a whole is also positive, meaning that the multimode version of the Stephanou et al. model [23] presented in this work is thermodynamically admissible.

2.6. Conformation Tensor Evolution Equation

The evolution equation used for each of the dimensionless conformation tensors is obtained by substituting Equations (4b), (9d), and (11) into Equation (7a)
c ˙ ˜ [ 1 ] ( i ) = 1 τ CCR ( i ) ( tr c ˜ ( i ) ) { α ( i ) ( 1 ξ ( i ) ) h 2 ( tr c ˜ ( i ) ) c ˜ ( i ) c ˜ ( i )         + [ 1 2 α ( i ) ( 1 ξ ( i ) ) ] h ( tr c ˜ ( i ) ) c ˜ ( i ) [ 1 α ( i ) ( 1 ξ ( i ) ) ] I } , i [ 1 , N ]
The CCR relaxation time for the ith mode is given in Equation (10), and the (dimensionless) effective spring constant is given in Equation (4c). Finally, the expression for the polymeric stress tensor is obtained by substituting Equations (4b) and (11) into Equation (8) as follows:
σ = i = 1 N G e ( i ) [ h ( tr c ˜ ( i ) ) c ˜ ( i ) I ]

3. Asymptotic Behavior of the Model in Steady State Shear

In this section, we provide analytical expressions that describe the asymptotic behavior of the multimode version of the Stephanou et al. [23] model in the limit of weak flows for the following cases: inception of simple shear flow (SSF) described by the kinematics u = ( γ ˙ y , 0 , 0 ) , in which γ ˙ is the (constant) shear rate; inception of uniaxial elongation flow (UEF) described by the kinematics u = ( ε ˙ x , 1 2 ε ˙ y , 1 2 ε ˙ z ) , in which ε ˙ is the (constant) elongation rate; and small amplitude oscillatory shear (SAOS) described by the kinematics u = ( γ ˙ cos ( ω t ) y , 0 , 0 ) , in which ω is the oscillation frequency. The material functions to be analyzed are as follows: (a) the transient shear viscosity η + ( t ) (= σ yx ( t ) / γ ˙ ) and the first, Ψ 1 + ( t ) (= ( σ xx ( t ) σ yy ( t ) ) / γ ˙ 2 ), and second, Ψ 2 + ( t ) (= ( σ yy ( t ) σ zz ( t ) ) / γ ˙ 2 ), normal stress coefficients in the case of shear; (b) the transient elongational viscosity, η E + ( t ) (= ( σ xx ( t ) σ yy ( t ) ) / ε ˙ ), in the case of uniaxial elongation; and (c) the storage, G ( ω ) , and loss, G ( ω ) , 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 analytically solve the corresponding ordinary differential equations. After this stage, we obtain the non-zero stress tensor components via Equation (14). Finally, we obtain the following results for the relevant material functions:
Inception of shear:
η + ( t ) = i = 1 N τ CR ( i ) G ˜ e ( i ) [ 1 exp ( t τ CR ( i ) ) ]
Ψ 1 + ( t ) = 2 i = 1 N ( τ CR ( i ) ) 2 G ˜ e ( i ) [ 1 ( 1 + t τ CR ( i ) ) exp ( t τ CR ( i ) ) ]
Ψ 2 + ( t ) = i = 1 N ( τ CR ( i ) ) 2 G ˜ e ( i ) { [ ξ ( i ) + α ( i ) ( 1 ξ ( i ) ) ] [ 1 ( 1 + t τ CR ( i ) ) exp ( t τ CR ( i ) ) ]    α ( i ) ( 1 ξ ( i ) ) 2 exp ( t τ CR ( i ) ) [ 1 t τ CR ( i ) exp ( t τ CR ( i ) ) ] }
Inception of uniaxial elongation:
η E + ( t ) = 3 i = 1 N τ CR ( i ) G ˜ e ( i ) [ 1 exp ( t τ CR ( i ) ) ]
meaning that Trouton’s law is true for the steady-state extensional viscosity.
Small Amplitude Oscillatory Shear:
G ( ω ) = i = 1 N G ˜ e ( i ) ( ω τ CR ( i ) ) 2 1 + ( ω τ CR ( i ) ) 2 G ( ω ) = i = 1 N G ˜ e ( i ) ω τ CR ( i ) 1 + ( ω τ CR ( i ) ) 2
In Equations (15)–(17), we have defined G ˜ e ( i ) ( 1 ξ ( i ) ) 2 G e ( i ) .

4. Results and Discussion

The FENE parameter can be easily calculated via b e = 3 ( 0.82 ) 2 M e / ( M 0 C ) , as mentioned above. In this study, PE M0 = 14 g/mol, whereas C = 7.3 and Me = 828 g/mol (see Table 3.3, p. 151 of Ref. [40]). These values yield be = 16.34. We will compare against the experimental data of Konaganti et al. [15] that have performed rheological measurements of the sample HDPE-1 (reported by the same group [41]), for which Mw = 206 kg/mol; thus, the number of entanglements is equal to Z ≈ 249 >> 1. The relaxation spectrum is the same as the one used by Konaganti et al. [15] (see their paper’s Table 2 for T = 200 °C, though is also provided in Table 1), and it was obtained by fitting the expressions of the storage and loss moduli, which are shown in Equation (17), with the corresponding experimental data. The comparison against the experimental storage and loss moduli is shown in Figure 1.
All of the remaining parameter values are obtained by fitting the model predictions with the experimental data; 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., Konaganti et al. [15]). The following values of the model parameter are chosen to provide a good comparison with the start-up shear flow experimental data provided in Figure 2: ξ = 0.02, α = 0.3, βccr = 4 × 10−4, and k = −3.5. For comparison, we also depict, in the following figures, the predictions of the multimode Giesekus model [which is a special case of our model in which βccr = ξ = k = 0 and be infinite or the function h = 1 in Equation (4c)] with α = 0.3 and the relaxation spectrum of Table 1.

4.1. Comparison with Start-Up Shear Flow Data

Figure 2a illustrates the comparison between the experimental data for the growth of the shear viscosity upon inception of shear flow and the simulated results obtained using the new model. The experimental data (symbols) were collected at three different shear rates—0.05 1/s, 0.5 1/s, and 1 1/s—whereas the 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, doing so much more successfully than the Giesekus model (Figure 2b). We should, however, note that the overshoots noted at the two larger shear rates (0.5 1/s and 1 1/s) are higher than the experimental data. The overshoot predicted using our model is controlled by two parameters: the FENE parameter be and the anisotropic mobility (Giesekus) parameter α. As mentioned above, the former parameter is not a free parameter, as it is directly dictated by structural parameters. The latter parameter is a free parameter, and by increasing its value, the start-up shear viscosity overshoots noted at the two larger shear rates (0.5 1/s and 1 1/s) shift downwards and broaden (results not shown), thus more closely agreeing with the experimental data; however, the good comparison identified at the smaller shear rate (0.05 1/s) is reduced. This result might hint that the parameter α should not be a constant, but should increase with the applied strain rate. Similar arguments have been put forth and resulted in a variable non-affine/slip parameter proposed by Nikiforidis et al. [42] and a variable link tension coefficient proposed by Stephanou and Kröger [26]. We note that although a non-zero value of ξ is employed, the undershoots produced are too small to be noted via the scale used in Figure 2. Although no experimental data are provided, we provide the corresponding prediction of the growth of the first and second normal stress coefficients in Figure 3, as well as 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 between the model predictions and the experimental measurements obtained by Konaganti et al. [15] for the growth of the elongational viscosity as a function of time upon inception of uniaxial elongation at the following stretch rates: 0.05 1/s, 0.5 1/s, and 5 1/s. The comparison employs the same parameter values as those used in Figure 2, except ξ = 0, as informed by Stephanou et al. [37], 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 appropriate. We note that the Giesekus model (Figure 5b) fails to reproduce correctly the trend of the experimental data. It should be noted that the uniaxial elongation data obtained by Konaganti et al. [15] do not seem to reach a steady state. This outcome is customary in the literature [11,12,43] due to experimental difficulties, as the polymer samples used become too slender and often break. This issue has led some researchers to consider the highest value as the steady-state elongational viscosity and omit all following data [44]. Thus, the data that follow the overshoot may not be accurate. 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 developed the multiple-mode version of a generalized, conformation-tensor based viscoelastic model for polymer melts, which was proposed in [23], by making use of the generalized bracket formalism of Beris and Edwards [36]. Like its forerunner, i.e., the single-mode version [23], it accounts for the most significant effects that can be realized in entangled polymer 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 have a very good predictive capability with regard to the industrial experimental data for HDPE obtained by Konaganti et al. [15]. Obviously, a better comparison could have been obtained if we had simultaneously fitted all available data, as Konaganti et al. [15] did (we note that some parameters must have different values in different flows, such as the non-affine/slip parameter, which must be explicitly considered in the fitting process). Furthermore, different values of the model parameters for each mode could also have been considered, following the work of Konaganti et al. [15], which would certainly provide more flexibility to the model and, thus, improve its capacity to fit the experimental data. However, in our present work, we mostly focused on deriving the multimode version of the model proposed by Stephanou et al. [23] using non-equilibrium thermodynamics, with less focus devoted to its predictive capacity to almost perfectly fit all available experimental data.
The model, in its present form, only considers 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 the presence of even low levels of long-chain branching [45]. We, therefore, need to generalize it and allow the explicit consideration of branches, following the guidelines provided by the pom-pom [27] model and its thermodynamically admissible version, known as the pom-pon [28] model. Furthermore, the multimode version does not explicitly consider the molecular weight (MW) distribution of industrial samples, such as the log-normal or gamma distributions, that are able to describe experimentally measured distributions [46]. As such, we should also generalize it to handle molecular weight distribution, following the work of Schieber [47,48]. This generalized constitutive model will allow a more accurate prediction of the rheological responses of industrially used polymeric systems that possess an extensive spectrum of MW. Finally, we used an ambiguous model-fitting process wherein the values of the model parameters were obtained to best compare them to the experimental data. However, atomistic non-equilibrium molecular dynamics simulations can be used to directly obtain some of the parameters from the simulations [37,49]. For example, the EWM exponent can be obtained by directly comparing the prediction of Equation (9b) to the relaxation time as a function of shear rate at large shear rates obtained via the NEMD simulations, and the CCR parameter can be obtained by directly comparing the average number of entanglements as a function of shear rate which, due to CCR, is noted to decrease [22] (we note, however, that in our present work we omitted this approach and considered a constant number of entanglements). Only then, polymer engineers would be able to accurately use the predictions of the multimode constitutive model for comparison with the rheological response noted in actual industrial processes. Our findings provide a foundation for future research that aims to enhance the properties of high MW polymers for diverse applications.

Author Contributions

Conceptualization, P.S.S.; Data curation, P.C.K. and P.S.S.; Formal analysis, P.C.K. and P.S.S.; Investigation, P.C.K. and P.S.S.; Methodology, P.S.S.; Project administration, P.S.S.; Resources, P.C.K. and P.S.S.; Software, P.C.K. and P.S.S.; Supervision, P.S.S.; Validation, P.C.K. and P.S.S.; Visualization, P.C.K. and P.S.S.; Writing—original draft, P.C.K.; Writing—review and editing, P.S.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Tadmor, Z.; Gogos, C.G. Principles of Polymer Processing, 2nd ed.; Wiley: Hoboken, NJ, USA, 2006. [Google Scholar]
  2. PlasticsToday, Global Plastic Packaging Market Worth $269.6 billion by 2025. Available online: https://www.plasticstoday.com/packaging/global-plastic-packaging-market-worth-2696-billion-2025 (accessed on 19 June 2023).
  3. Graham, M.D. The sharkskin instability of polymer melt flows. Chaos 1999, 9, 154–163. [Google Scholar] [CrossRef]
  4. Barone, J.R.; Plucktaveesak, N.; Wang, S.Q. Interfacial molecular instability mechanism for sharkskin phenomenon in capillary extrusion of linear polyethylenes. J. Rheol. 1998, 42, 813–832. [Google Scholar] [CrossRef]
  5. Baig, C.; Stephanou, P.S.; Tsolou, G.; Mavrantzas, V.G.; Kröger, M. Understanding dynamics in binary mixtures of entangled cis- 1,4-polybutadiene melts at the level of primitive path segments by mapping atomistic simulation data onto the tube model. Macromolecules 2010, 43, 8239–8250. [Google Scholar] [CrossRef]
  6. Stephanou, P.S.; Baig, C.; Tsolou, G.; Mavrantzas, V.G.; Kröger, M. Quantifying chain reptation in entangled polymer melts: Topological and dynamical mapping of atomistic simulation results onto the tube model. J. Chem. Phys. 2010, 132, 124904. [Google Scholar] [CrossRef] [Green Version]
  7. Stephanou, P.S.; Baig, C.; Mavrantzas, V.G. Projection of atomistic simulation data for the dynamics of entangled polymers onto the tube theory: Calculation of the segment survival probability function and comparison with modern tube models. Soft Matter 2011, 7, 380–395. [Google Scholar] [CrossRef]
  8. Stephanou, P.S.; Baig, C.; Mavrantzas, V.G. Toward an improved description of constraint release and contour length fluctuations in tube models for entangled polymer melts guided by atomistic simulations. Macromol. Theory Simul. 2011, 20, 752–768. [Google Scholar] [CrossRef]
  9. Stephanou, P.S.; Mavrantzas, V.G. Accurate prediction of the linear viscoelastic properties of highly entangled mono and bidisperse polymer melts. J. Chem. Phys. 2014, 140, 214903. [Google Scholar] [CrossRef]
  10. Nassehi, V. Practical Aspects of Finite Element Modelling of Polymer Processing; Wiley: Hoboken, NJ, USA, 2002. [Google Scholar]
  11. Larson, R.G. Constitutive Equations for Polymer Melts and Solutions, 1st ed.; Butterworth-Heinemann: Oxford, UK, 1988. [Google Scholar]
  12. Bird, R.B.; Armstrong, R.C.; Hassager, O. Dynamics of Polymeric Liquids. Volume 1. Fluid Mechanics, 2nd ed.; Wiley-Interscience: Hoboken, NJ, USA, 1987. [Google Scholar]
  13. Ferry, J.D. Viscoelastic Properties of Polymers, 1st ed.; John Wiley & Sons, Ltd.: Hoboken, NJ, USA, 1980. [Google Scholar]
  14. Luo, X.-L.; Tanner, R.I. Finite element simulation of long and short circular die extrusion experiments using integral models. Int. J. Numer. Methods Eng. 1988, 25, 9–22. [Google Scholar] [CrossRef]
  15. Konaganti, V.K.; Ansari, M.; Mitsoulis, E.; Hatzikiriakos, S.G. Extrudate swell of a high-density polyethylene melt: II. Modeling using integral and differential constitutive equations. J. Nonnewton. Fluid Mech. 2015, 225, 94–105. [Google Scholar] [CrossRef]
  16. Langouche, F.; Debbaut, B. Rheological characterisation of a high-density polyethylene with a multi-mode differential viscoelastic model and numerical simulation of transient elongational recovery experiments. Rheol. Acta 1999, 38, 48–64. [Google Scholar] [CrossRef]
  17. Shiromoto, S.; Masutani, Y.; Tsutsubuchi, M.; Togawa, Y.; Kajiwara, T. The effect of viscoelasticity on the extrusion drawing in film-casting process. Rheol. Acta 2010, 49, 757–767. [Google Scholar] [CrossRef]
  18. Khan, S.A.; Larson, R.G. Comparison of Simple Constitutive Equations for Polymer Melts in Shear and Biaxial and Uniaxial Extensions. J. Rheol. 1987, 31, 207–234. [Google Scholar] [CrossRef]
  19. Doi, M.; Edwards, S.F. The Theory of Polymer Dynamics, 1st ed.; Clarendon Press: Oxford, UK, 1988; Available online: https://global.oup.com/academic/product/the-theory-of-polymer-dynamics-9780198520337 (accessed on 17 May 2023).
  20. Ianniruberto, G.; Marrucci, G. On compatibility of the Cox-Merz rule with the model of Doi and Edwards. J. Nonnewton. Fluid Mech. 1996, 65, 241–246. [Google Scholar] [CrossRef]
  21. Marrucci, G. Dynamics of entanglements: A nonlinear model consistent with the Cox-Merz rule. J. Nonnewton. Fluid Mech. 1996, 62, 279–289. [Google Scholar] [CrossRef]
  22. Baig, C.; Mavrantzas, V.G.; Kröger, M. Flow effects on melt structure and entanglement network of linear polymers: Results from a nonequilibrium molecular dynamics simulation study of a polyethylene melt in steady shear. Macromolecules 2010, 43, 6886–6902. [Google Scholar] [CrossRef]
  23. Stephanou, P.S.; Tsimouri, I.C.; Mavrantzas, V.G. Simple, accurate and user-friendly differential constitutive model for the rheology of entangled polymer melts and solutions from nonequilibrium thermodynamics. Materials 2020, 13, 2867. [Google Scholar] [CrossRef]
  24. Stephanou, P.S.; Schweizer, T.; Kröger, M. Communication: Appearance of undershoots in start-up shear: Experimental findings captured by tumbling-snake dynamics. J. Chem. Phys. 2017, 146, 161101. [Google Scholar] [CrossRef] [Green Version]
  25. Costanzo, S.; Huang, Q.; Ianniruberto, G.; Marrucci, G.; Hassager, O.; Vlassopoulos, D. Shear and Extensional Rheology of Polystyrene Melts and Solutions with the Same Number of Entanglements. Macromolecules 2016, 49, 3925–3935. [Google Scholar] [CrossRef] [Green Version]
  26. Stephanou, P.S.; Kröger, M. Non-constant link tension coefficient in the tumbling-snake model subjected to simple shear. J. Chem. Phys. 2017, 147, 174903. [Google Scholar] [CrossRef]
  27. Mcleish, T.C.B.; Larson, R.G. Molecular constitutive equations for a class of branched polymers: The pom-pom polymer. J. Rheol. 1998, 42, 81–110. [Google Scholar] [CrossRef]
  28. Öttinger, H.C. Thermodynamic admissibility of the pompon model for branched polymers. Rheol. Acta 2001, 40, 317–321. [Google Scholar] [CrossRef]
  29. Inkson, N.J.; McLeish, T.C.B.; Harlen, O.G.; Groves, D.J. Predicting low density polyethylene melt rheology in elongational and shear flows with “pom-pom” constitutive equations. J. Rheol. 1999, 43, 873–896. [Google Scholar] [CrossRef]
  30. Soulages, J.; Schweizer, T.; Venerus, D.C.; Kröger, M.; Öttinger, H.C. Lubricated cross-slot flow of a low density polyethylene melt. J. Nonnewton. Fluid Mech. 2008, 154, 52–64. [Google Scholar] [CrossRef]
  31. Verbeeten, W.M.H.; Peters, G.W.M.; Baaijens, F.P.T. Differential constitutive equations for polymer melts: The extended Pom–Pom model. J. Rheol. 2001, 45, 823–843. [Google Scholar] [CrossRef] [Green Version]
  32. Soulages, J.; Hütter, M.; Öttinger, H.C. Thermodynamic admissibility of the extended Pom-Pom model for branched polymers. J. Nonnewton. Fluid Mech. 2006, 139, 209–213. [Google Scholar] [CrossRef]
  33. Hoyle, D.M.; Harlen, O.G.; Auhl, D.; McLeish, T.C.B. Non-linear step strain of branched polymer melts. J. Rheol. 2009, 53, 917–942. [Google Scholar] [CrossRef]
  34. Clemeur, N.; Rutgers, R.P.; Debbaut, B. On the evaluation of some differential formulations for the pom-pom constitutive model. Rheol. Acta 2003, 42, 217–231. [Google Scholar] [CrossRef]
  35. Stephanou, P.S.; Tsimouri, I.C.; Mavrantzas, V.G. Flow-Induced Orientation and Stretching of Entangled Polymers in the Framework of Nonequilibrium Thermodynamics. Macromolecules 2016, 49, 3161–3173. [Google Scholar] [CrossRef]
  36. Beris, A.N.; Edwards, B.J. Thermodynamics of Flowing Systems: With Internal Microstructure; Oxford University Press: New York, NY, USA, 1994. [Google Scholar]
  37. Stephanou, P.S.; Baig, C.; Mavrantzas, V.G. A generalized differential constitutive equation for polymer melts based on principles of nonequilibrium thermodynamics. J. Rheol. 2009, 53, 309–337. [Google Scholar] [CrossRef]
  38. Öttinger, H.C. Beyond Equilibrium Thermodynamics; John Wiley and Sons: Hoboken, NJ, USA, 2005. [Google Scholar]
  39. Marrucci, G.; Ianniruberto, G. Flow-induced orientation and stretching of entangled polymers. Philos. Trans. R. Soc. London. Ser. A Math. Phys. Eng. Sci. 2003, 361, 677–688. [Google Scholar] [CrossRef]
  40. Larson, R.G. The Structure and Rheology of Complex Fluids; Oxford University Press: New York, NY, USA, 1999. [Google Scholar]
  41. Behzadfar, E.; Ansari, M.; Konaganti, V.K.; Hatzikiriakos, S.G. Extrudate swell of HDPE melts: I. Experimental. J. Nonnewton. Fluid Mech. 2015, 225, 86–93. [Google Scholar] [CrossRef]
  42. Nikiforidis, V.-M.; Tsalikis, D.G.; Stephanou, P.S. On The Use of a Non-Constant Non-Affine or Slip Parameter in Polymer Rheology Constitutive Modeling. Dynamics 2022, 2, 380–398. [Google Scholar] [CrossRef]
  43. Rasmussen, H.K.; Wingstrand, S.L.; Hassager, O. On the universality in the extensional rheology of monodisperse polymer melts and oligomer dilutions thereof. Rheol. Acta 2019, 58, 333–340. [Google Scholar] [CrossRef]
  44. Huang, Q.; Hengeller, L.; Alvarez, N.J.; Hassager, O. Bridging the Gap between Polymer Melts and Solutions in Extensional Rheology. Macromolecules 2015, 48, 4158–4163. [Google Scholar] [CrossRef] [Green Version]
  45. Wood-Adams, P.M. The effect of long chain branches on the shear flow behavior of polyethylene. J. Rheol. 2001, 45, 203–210. [Google Scholar] [CrossRef]
  46. Rogošić, M.; Mencer, H.J.; Gomzi, Z. Polydispersity index and molecular weight distributions of polymers. Eur. Polym. J. 1996, 32, 1337–1344. [Google Scholar] [CrossRef]
  47. Schieber, J.D. Kinetic theory of polymer melts. VIII. Rheological properties of polydisperse mixtures. J. Chem. Phys. 1987, 87, 4917–4927. [Google Scholar] [CrossRef]
  48. Schieber, J.D. Kinetic theory of polymer melts. IX. Comparisons with experimental data. J. Chem. Phys. 1987, 87, 4928–4936. [Google Scholar] [CrossRef]
  49. Stephanou, P.S.; Tsalikis, D.G.; Skountzos, E.N.; Mavrantzas, V.G. Understanding the rheological behavior of polymer nanocomposites: Non-equilibrium thermodynamics modeling coupled with detailed atomistic non-equilibrium molecular dynamics simulations. Mater. Today Proc. 2018, 5, 27589–27598. [Google Scholar] [CrossRef]
Figure 1. Comparison between the model predictions Equation (17) and the experimental data presented in Ref. [15] for the storage and loss moduli of an HDPE sample at 200 °C. The relaxation spectrum is provided in Table 1.
Figure 1. Comparison between the model predictions Equation (17) and the experimental data presented in Ref. [15] for the storage and loss moduli of an HDPE sample at 200 °C. The relaxation spectrum is provided in Table 1.
Polymers 15 03322 g001
Figure 2. (a) Model prediction of the growth of the shear viscosity prediction (lines) upon the inception of the shear flow at three different dimensionless shear rates, along with comparison with the experimental data (symbols) considered in Ref. [15]. In panel (b), we depict the same comparison for the multimode Giesekus model. The thick black line depicts the LVE envelope, which is shown in Equation (15a).
Figure 2. (a) Model prediction of the growth of the shear viscosity prediction (lines) upon the inception of the shear flow at three different dimensionless shear rates, along with comparison with the experimental data (symbols) considered in Ref. [15]. In panel (b), we depict the same comparison for the multimode Giesekus model. The thick black line depicts the LVE envelope, which is shown in Equation (15a).
Polymers 15 03322 g002
Figure 3. Model predictions (red, blue, and orange lines) of the growth of the first (a) and second (b) normal stress coefficients upon the inception of shear flow. The thick black lines depict the LVE envelope, which is shown in Equations (15b) and (15c). The s parameter values are the same as those used in in Figure 2. The multimode Giesekus model’s predictions are also depicted (red, blue, and orange dashed lines). The thick black dashed lines depict the LVE envelope, which are again shown in Equations (15b) and (15c).
Figure 3. Model predictions (red, blue, and orange lines) of the growth of the first (a) and second (b) normal stress coefficients upon the inception of shear flow. The thick black lines depict the LVE envelope, which is shown in Equations (15b) and (15c). The s parameter values are the same as those used in in Figure 2. The multimode Giesekus model’s predictions are also depicted (red, blue, and orange dashed lines). The thick black dashed lines depict the LVE envelope, which are again shown in Equations (15b) and (15c).
Polymers 15 03322 g003
Figure 4. Model predictions of the steady-state (a) shear viscosity and the (b) first and (c) second normal stress coefficients of the HDPE-1 sample. The parameter values are the same as those used in Figure 2. The multimode Giesekus model prediction is also depicted (black dashed lines).
Figure 4. Model predictions of the steady-state (a) shear viscosity and the (b) first and (c) second normal stress coefficients of the HDPE-1 sample. The parameter values are the same as those used in Figure 2. The multimode Giesekus model prediction is also depicted (black dashed lines).
Polymers 15 03322 g004
Figure 5. (a) Comparison between the model predictions (lines) and the experimental measurements (symbols) of Konaganti et al. [15] for the growth of the elongational viscosity as a function of time for several stretch rates. The thick black line depicts the LVE envelope, as given in Equation (16). The sparameter values are the same as those used in Figure 2 except ξ = 0. In panel (b), we depict the same comparison for the multimode Giesekus model.
Figure 5. (a) Comparison between the model predictions (lines) and the experimental measurements (symbols) of Konaganti et al. [15] for the growth of the elongational viscosity as a function of time for several stretch rates. The thick black line depicts the LVE envelope, as given in Equation (16). The sparameter values are the same as those used in Figure 2 except ξ = 0. In panel (b), we depict the same comparison for the multimode Giesekus model.
Polymers 15 03322 g005
Figure 6. Model prediction of the steady-state uniaxial elongational viscosity. The parameter values are the same as those used in Figure 2 except ξ = 0. The multimode Giesekus model prediction is also depicted (black dashed line).
Figure 6. Model prediction of the steady-state uniaxial elongational viscosity. The parameter values are the same as those used in Figure 2 except ξ = 0. The multimode Giesekus model prediction is also depicted (black dashed line).
Polymers 15 03322 g006
Table 1. Relaxation spectrum [15].
Table 1. Relaxation spectrum [15].
Mode G ˜ e ( i )   ( Pa ) τ CR ( i )   ( s )
1387,8080.00086
2185,3070.0075
393,3380.0548
437,7660.403
512,9342.99
6502530.78
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Konstantinou, P.C.; Stephanou, P.S. Predicting High-Density Polyethylene Melt Rheology Using a Multimode Tube Model Derived Using Non-Equilibrium Thermodynamics. Polymers 2023, 15, 3322. https://doi.org/10.3390/polym15153322

AMA Style

Konstantinou PC, Stephanou PS. Predicting High-Density Polyethylene Melt Rheology Using a Multimode Tube Model Derived Using Non-Equilibrium Thermodynamics. Polymers. 2023; 15(15):3322. https://doi.org/10.3390/polym15153322

Chicago/Turabian Style

Konstantinou, Pavlina C., and Pavlos S. Stephanou. 2023. "Predicting High-Density Polyethylene Melt Rheology Using a Multimode Tube Model Derived Using Non-Equilibrium Thermodynamics" Polymers 15, no. 15: 3322. https://doi.org/10.3390/polym15153322

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop