Simple, Accurate and User-Friendly Differential Constitutive Model for the Rheology of Entangled Polymer Melts and Solutions from Nonequilibrium Thermodynamics

In a recent reformulation of the Marrucci-Ianniruberto constitutive equation for the rheology of entangled polymer melts in the context of nonequilibrium thermodynamics, rather large values of the convective constraint release parameter βccr had to be used in order for the model not to violate the second law of thermodynamics. In this work, we present an appropriate modification of the model, which avoids the splitting of the evolution equation for the conformation tensor into an orientation and a stretching part. Then, thermodynamic admissibility simply dictates that βccr ≥ 0, thus allowing for more realistic values of βccr to be chosen. Moreover, and in view of recent experimental evidence for a transient stress undershoot (following the overshoot) at high shear rates, whose origin may be traced back to molecular tumbling, we have incorporated additional terms into the model accounting, at least in an approximate way, for non-affine deformation through a slip parameter ξ. Use of the new model to describe available experimental data for the transient and steady-state shear and elongational rheology of entangled polystyrene melts and concentrated solutions shows close agreement. Overall, the modified model proposed here combines simplicity with accuracy, which renders it an excellent choice for managing complex viscoelastic fluid flows in large-scale numerical calculations.


Introduction
Using accurate constitutive models to describe the complex rheological response of high molecular weight (MW) polymers (melts or solutions) to an applied flow field allows for the more economic, rational design of new polymer-based products as well as for the improved description of processing operations as these models are key ingredients in large-scale numerical calculations of viscoelastic fluid flows. A successful model should be able to accurately predict changes in the rheology of the material during processing as a function of chemical structure and composition, molecular architecture, and morphology (possible co-existence of crystalline and amorphous domains).
Since its introduction and after several modifications and improvements, the tube model [1,2] has offered a remarkable conceptual framework for formulating theories that can explain but also 2 of 18 quantitatively describe the behavior of entangled, high-MW polymer melts and concentrated polymer solutions, both under equilibrium and nonequilibrium (flow) conditions. It relies on the idea that due to chain uncrossability and the development of entanglements with other chains, the lateral motion of a reference chain is effectively constrained within a curvilinear tube-like region (a mean-field tube) around the primitive path of the chain with a diameter whose magnitude reflects the strength of these topological interactions. In the simplest case of a linear chain that is restricted within such a tube-like region, relaxation is achieved through backward and forward diffusion (i.e., through reptation) along the axis of the confining tube [1,2]. Today, chain reptation within the mean-field tube is not simply a mathematical invention; it has been directly observed in molecular simulations [3]. However, it is also well-known that for the tube model to be able to quantitatively describe the equilibrium dynamics of actual polymeric systems, it has to account for additional mechanisms, such as contour length fluctuations (CLFs) due to the "breathing" motion of chain ends, and constraint release (CR) due to the destruction of topological constraints as a result of the simultaneous (collective) relaxation of the other chains surrounding the reference chain.
A mechanism that the original theory failed to consider under strong flows is that of convective constraint release (CCR). In this case, when the shear rate exceeds 1/τ d where τ d is the inverse chain reptation or disentanglement time, the simultaneous motion of the surrounding chains comprising the mean-field tube leads to the release of entanglements, thus inducing a faster relaxation [4,5]. Self-consistently, this means that the number of entanglements should be a decreasing function of the strength of the applied flow, a notion that was originally missed by polymer physicists and rheologists. This has been directly demonstrated in the nonequilibrium molecular dynamics (NEMD) simulations of Baig et al. [6] with an entangled polyethylene melt under shear flow, as well as in several other NEMD simulation studies afterwards [7][8][9][10]. The direct observation of flow-induced chain disentanglement was incorporated a few years later into a novel and very successful constitutive rheological model by Ianniruberto and Marrucci accounting explicitly for the CCR mechanism [11]. Chain stretching was also considered [12] by proposing two separate evolution equations: one for the orientation tensor S describing the average orientation of tube segments, and one for the scalar variable λ accounting for the average chain stretch. However, Wapperom et al. [13] argued that such a model leads to anomalous behavior when employed in numerical simulations of complex flows. This motivated Ianniruberto and Marrucci [14] to propose a suitable modification that unified the two different equations into a single one through the use of a single conformation tensor.
Another important mechanism that was originally neglected but was again pointed out by NEMD simulations to be of relevance is rotation or tumbling even of entangled chains under shear [6][7][8]15]. Shear-induced tumbling of entangled polymer chains was incorporated into a constitutive model by Costanzo et al. [16] by introducing a "tumbling" function. In fact, the authors showed that this tumbling behavior is responsible for the appearance of a transient stress undershoot (following the overshoot) at high shear rates in the case of polymer solutions but not in the case of polymer melts. A similar conclusion was reached by Stephanou et al. [17,18] who revisited the kinetic theory of Curtiss-Bird [19][20][21] and solved the corresponding model equations [22] in the presence of rotational diffusion. They showed that the tumbling-snake model for polymeric systems, as they coined the kinetic theory of Curtiss-Bird, predicts pronounced undershoots in the time-dependent shear viscosity (arising from the rotational Brownian diffusion term in the equation of motion for polymer segments) but not in the transient normal stress coefficients. Obviously, such undershoots should be absent in the transient extensional viscosities, since elongational flows are by default free of any rotational contributions to the velocity gradient tensor, and this was indeed verified to be the case [16,23,24]. Additional results from NEMD simulations demonstrated that the overshoot and undershoot in the shear viscosity upon flow startup are correlated to tube orientation and chain tumbling, respectively [9,10].
In a recent contribution [25], we formally re-derived the Ianniruberto-Marrucci [14] and Rolie-Poly [26] models in the context of the generalized bracket [27] and GENERIC [28][29][30] formalisms of nonequilibrium thermodynamics (NET), and appropriately extended them to account, in a self-consistent manner, for a second normal stress difference. However, by applying the 2nd law of thermodynamics and requiring the evolution equation to preserve the positive -definite nature of the conformation tensor between successive entanglement points along a chain for all times and all flow fields was found to constrain the so-called CCR parameter β ccr to rather large values (larger than one), which seems to be too restrictive. In the present article, we provide a reformulation of that model that allows the use of more realistic values of the CCR parameter. We also introduce in the model terms that account for non-affine motion of entanglement strands (allowing for chain tumbling) through a slip parameter ξ. The use of nonequilibrium thermodynamics guarantees that this extra mechanism is self-consistently incorporated in the constitutive model, which also ensures the internal consistency of the final set of transport equations.
The structure of the rest of the paper is as follows: the nonequilibrium thermodynamics formulation of the new (refined) model in the context of the generalized bracket formalism [27] is presented in the next section; the proof of its thermodynamic admissibility is detailed in Appendix A. We note here that at the level of description of interest in this work, the generalized bracket, and the GENERIC [28][29][30] formalisms are equivalent. In Section 3, we derive asymptotic solutions of the model in the cases of steady-state and transient shear and elongational flows, in the limit of weak flows. We also discuss how closely the new model can describe experimental rheological data already reported in the literature over the years both for entangled polymer melts and concentrated polymer solutions. We conclude with Section 4 presenting a summary of the most important results of our study and a brief discussion of future plans.

Nonequilibrium Thermodynamics Refinement of the Stephanou et al. Model
Following Stephanou et al. [25], in the modified version of our model, we also make use of only one structural parameter, the entanglement strand conformation tensor c, which is made dimensionless throughc = K c/k B T where K denotes the spring constant of the Hookean dumbbells representing entanglement strands, k B the Boltzmann constant, and T the absolute temperature. The conformation tensorc refers to one entanglement strand and at equilibrium (zero flow field applied) coincides by definition with the unit tensor I. Contrary to our previous model, however, we consider an extra contribution to the dissipative bracket, in addition to relaxation phenomena at the level of the structural variable, namely strand non-affine motion, introduced through a fourth-rank tensor L. Then, using the well-known procedure of the generalized bracket formalism [27], the resulting evolution equation for c reads: . c αβ, [1] where, on the left hand-side, we have defined the upper-convected Maxwell time derivative: .
In Equation (1), A denotes the Helmholtz free energy of the system, H m the mechanical part of the system Hamiltonian, Λ the relaxation fourth-rank matrix, and L the fourth-rank matrix associated with the non-affine motion mechanism. The corresponding expression for the stress tensor is [27].
In Equations (1)-(3), we have assumed Einstein's implicit summation convention for any repeated Greek indices (and this will be followed throughout this work). We will also restrict our analysis to incompressible and isothermal fluids. The mechanical part of the Hamiltonian is simply equal to the kinetic energy plus the Helmholtz free energy: The latter is given as [25,31,32] and includes two terms: the dimensionless potential Φ tr(c − I) accounting for chain stretching, and an entropic contribution involving the logarithm of the determinant of the conformation tensor.
In Equation (4b), G e = n e k B T is the entanglement modulus with n e being the entanglement density (assumed to be constant, independent of the flow). The partial derivative of the potential Φ tr(c − I) with respect to the trace ofc defines the (dimensionless) effective spring constant h [27,30,33]: where with A = A/G e we have denoted the dimensionless free energy. The corresponding Volterra derivative of the dimensionless free energy with respect toc is In Ref. [25], we had employed the FENE-P (Cohen) approximation [34] rather than the FENE-P (Warner) one, as it is a better approximation to the inverse Langevin function [33,34]. However, Varchanis et al. [35] found that the FENE-P (Cohen) force law leads to two different solutions, both stable and both preserving the positive definite nature of the conformation tensor, with one, however, being aphysical. Motivated by this finding, in the present study we chose to work with the FENE-P (Warner) approximation for which [33]: Note that Equation (7a) differs slightly from Equation (23a) in Stephanou et al. [33]: although both equations lead to the same expression for the force law, Equation (7a) vanishes at equilibrium as it should, whereas Equation (23a) in Stephanou et al. [33] does not. In the above expressions, b e ≡ 3L 2 e / a 2 e is the extensibility parameter with L e denoting the maximum length of the entanglement strand. The length of an entanglement strand is equal to L/(Z + 1), where L is the primitive path contour length and Z the number of entanglements per chain. For given polymer chemistry, the extensibility parameter should not be considered as an adjustable (free) parameter but rather as a known constant, equal to b e ≡ 3(0.82) 2 /C ∞ (M e /M 0 ) where C ∞ is the polymer characteristic ratio at infinite chain length, M e the entanglement molecular weight, and M 0 the average molar mass per backbone bond, which is equal to half the monomer molecular weight (note the correction in Appendix A of Ref. [25]). For example, for polystyrene (PS): M 0 = 52 g/mol and C ∞ = 9.6 (Ref. [36], Table 3.3), hence b e = M e /298 g/mol. Consequently, for PS melts for which M e ≈ 13,300 g/mol (Ref. [36], Table 3.3), we find b e = 54, whereas for entangled polymer solutions the value of b e will depend on the polymer concentration of the solution (see Section 3.2.2). In order to be able to predict a transient stress undershoot (following the overshoot) at high shear rates, the new model accounts for non-affine deformation effects in Equation (1) through the tensor L, for which we consider the following form [27]: where ξ denotes the non-affine or slip parameter. The same form has been used by Beris et al. [37] to capture the thixotropic behavior of concentrated star polymer suspensions. Finally, the stress tensor is obtained by substituting Equations (6) and (8) into Equation (3), and reads: Defining the relaxation matrix Λ is the only remaining task. To this, we consider two relaxation times, the disentanglement or reptation time τ d , and the Rouse time τ R , implying that two relaxation matrices should be defined: 2τ R G e c αγ β βε + c αε β βγ + c βγ β αε + c βε β αγ (10) In Equation (10),β is the (dimensionless) mobility tensor (see Ref. [25]) and f τ d and f τ R scalar functions of the trace ofc defined as: The (total) relaxation matrix Λ needed in Equation (1) is then Λ = Λ τ d + Λ τ R . It is of interest to discuss how these expressions for the relaxation matrices differ from the expression used in the previous work, Equation (10a) in Ref. [25]. In the previous version of our model, the characteristic relaxation time associated with changes in the stretching of the entanglement segment (i.e., of the trace of the conformation tensor) was the Rouse time, which necessitated that Λ τ d be selected such that Λ τ d ααγε = 0. To avoid imposing strict restrictions on the CCR parameter in the present work (to guarantee the thermodynamic admissibility of the model, see below), the characteristic relaxation time associated with stretching is taken to be the CCR relaxation time, Equation (12c).
Upon inserting Equations (4a), (6), (8), and (10) in Equation (1), we obtain the following evolution equation for the dimensionless conformation tensor: where, on the left hand-side, we have defined the Johnson-Segalman time-derivative: while the CCR relaxation time is given as: with β ccr being the effective CCR parameter. Note that in order for the model to account for both reptation and constraint release (at least in a crude way), we have invoked the double reptation Materials 2020, 13, 2867 6 of 18 approximation and considered the equilibrium CCR relaxation time to be half the corresponding reptation time when only reptation is considered, i.e., 1 2 τ d [14]. When a non-zero value of the CCR parameter is considered, the CCR relaxation time τ ccr will exhibit a flow-induced reduction from its equilibrium value (i.e., 1 2 τ d ) towards its corresponding value at large flow rates (i.e., the Rouse time). Finally, for the mobility tensor, we invoke Giesekus' postulate that β = I + ασ [25,33], where α is the anisotropic mobility parameter andσ = σ/G e denotes the dimensionless stress tensor. Thus, the evolution equation, Equation (12a), becomes (in the following, for the purpose of simplifying notation, we will drop the arguments in the functions h tr(c − I) and τ ccr trc ): .
As far as the conditions that should be obeyed in order for the model to be thermodynamically admissible and preserve the positive-definite character of the conformation tensor, these are elaborated in Appendix A and require that: According to these, the proposed revised model allows for any non-negative value of the CCR parameter to be used without violating the laws of thermodynamics, which should be contrasted with the earlier version [25], wherein β ccr had to be strictly larger than unity. Finally, note that in deriving the above constitutive model, we have assumed that the contribution of the solvent to the total stress is too small compared to the contribution of polymer (although it can be easily incorporated if needed).

Asymptotic Solutions
We proceed next to analyze the asymptotic behavior of the new model in the limit of weak flows for the following three cases: inception of simple shear flow (SSF) described by the kinematics u = To get asymptotic expressions for the conformation tensor and consequently for the material functions in the limit of small strain rates, we invoke a linearization of the evolution equation for the conformation tensor, which in a second step is solved analytically to also provide the stress tensor components. Then, the following results are obtained for the relevant viscometric or material functions: Inception of shear: where η 0 = (1-ξ) 2 G e 1 2 τ d and Ψ 1,0 = η 0 τ d are also the zero-rate steady-state values of the shear viscosity and of the first normal stress coefficient: lim Materials 2020, 13, 2867 Inception of uniaxial elongation: and thus, Trouton's law holds for the steady-state extensional viscosity, i.e.: Small Amplitude Oscillatory flow:

Results and Discussion
In this section, we will use the FENE-P (Warner) approximation for the function h [Equation (7b)]. In Section 3.1, we will also employ the relation τ d = 3Zτ R [2], take b e = 100, and keep the number of entanglements Z equal to 20.

Material Functions in Shear
In Figure 1, we depict the model predictions for the steady-state shear viscosity ( Figure 1a) and the first ( Figure 1b) and second ( Figure 1c) normal stress coefficients (appropriately scaled with the zero-rate viscosity η 0 and first normal stress coefficient Ψ 1,0 ), as a function of the applied dimensionless shear rate Wi ≡ . γ 0 τ d ; results for various values of the model parameters are shown. When a vanishing value of the anisotropic parameter α is considered, then a vanishing second normal stress coefficient is predicted (the corresponding curves are missing from Figure 1c) As the CCR parameter increases, both η and Ψ 1 exhibit a faster shear thinning behavior, which is the expected behavior given that the CCR parameter introduces a shear-induced decrease of the relaxation time.
However, and irrespective of the value of the CCR parameter assumed, the power-law behavior at large shear rates remains unaltered, with the corresponding exponents equal to −2/3 for the shear viscosity and equal to −4/3 for the first normal stress coefficient. Then, while keeping β ccr constant, upon increasing the value of α from 0 to 0.1, the power-law exponent at large shear rates changes to −1 and −5/3 for η and Ψ 1 , respectively. By further increasing the value of α to 0.2, the shear viscosity is affected only slightly, while the first normal stress coefficient is not affected at all. In the case of the second normal stress coefficient, we note a non-vanishing value, which at small shear rates reaches the zero-rate value dictated by Equation (16c). Thus, as α increases, the zero-rate value increases. However, the large shear rate power-law exponent remains equal to −2 irrespective of the values of the parameters. Finally, for non-zero values of the non-affine parameter ξ, we note that until about Wi = 100, the predictions for η and Ψ 1 are the same as when ξ = 0, but after that the power-law behavior differs, the corresponding exponent for both η and Ψ 1 becomes now equal to −2. On the other hand, the prediction for Ψ 2 remains the same. parameter increases, both η and Ψ1 exhibit a faster shear thinning behavior, which is the expected behavior given that the CCR parameter introduces a shear-induced decrease of the relaxation time. However, and irrespective of the value of the CCR parameter assumed, the power-law behavior at large shear rates remains unaltered, with the corresponding exponents equal to −2/3 for the shear viscosity and equal to −4/3 for the first normal stress coefficient. Then, while keeping βccr constant, upon increasing the value of α from 0 to 0.1, the power-law exponent at large shear rates changes to −1 and −5/3 for η and Ψ1, respectively. By further increasing the value of α to 0.2, the shear viscosity is affected only slightly, while the first normal stress coefficient is not affected at all. In the case of the second normal stress coefficient, we note a non-vanishing value, which at small shear rates reaches the zero-rate value dictated by Equation (16c). Thus, as α increases, the zero-rate value increases. However, the large shear rate power-law exponent remains equal to −2 irrespective of the values of the parameters. Finally, for non-zero values of the non-affine parameter ξ, we note that until about Wi = 100, the predictions for η and Ψ1 are the same as when ξ = 0, but after that the power-law behavior differs, the corresponding exponent for both η and Ψ1 becomes now equal to −2. On the other hand, the prediction for Ψ2 remains the same.
In Figures 2 and 3, we present the model predictions for the growth of the shear viscosity and the first normal stress coefficient, respectively, upon inception of the shear flow at three different values of the dimensionless shear rate Wi (= 1, 10, and 100), along with the prediction of the LVE behavior given by Equations (15) (dotted orange line). Irrespective of the shear rate and the values of the parameters assumed, at sufficiently short times, all viscometric curves follow the LVE prediction, which is the expected theoretical prediction. At small shear rates (Wi = 1), both functions are seen to approach their steady-state values monotonically, whereas at larger shear rates (Wi = 10), they go through an overshoot before they reach their steady-state values. Finally, at large shear rates (Wi = 100) and for a non-zero value of the non-affine parameter, the shear viscosity goes through an undershoot right after the overshoot, which is more pronounced in the curves of Figure 2b. Of course, In Figures 2 and 3, we present the model predictions for the growth of the shear viscosity and the first normal stress coefficient, respectively, upon inception of the shear flow at three different values of the dimensionless shear rate Wi (= 1, 10, and 100), along with the prediction of the LVE behavior given by Equations (15) (dotted orange line). Irrespective of the shear rate and the values of the parameters assumed, at sufficiently short times, all viscometric curves follow the LVE prediction, which is the expected theoretical prediction. At small shear rates (Wi = 1), both functions are seen to approach their steady-state values monotonically, whereas at larger shear rates (Wi = 10), they go through an overshoot before they reach their steady-state values. Finally, at large shear rates (Wi = 100) and for a non-zero value of the non-affine parameter, the shear viscosity goes through an undershoot right after the overshoot, which is more pronounced in the curves of Figure 2b. Of course, this is only a rather small undershoot; more pronounced undershoots will be reported in the next section, when we will discuss how the model can fit true experimental data (see Figures 7a and 11b). Furthermore, note that both the intensity and the position of this undershoot are controlled by the parameters of the model. Such a behavior is not noticed for the first normal stress coefficient. Overall, the same trend is noticed also for the second normal stress coefficient (Figure 4) whereas, contrary to the other two viscometric functions, the transient curves go over (instead of below) the LVE predictions. Overall, it is a remarkable attribute of the new, single-mode viscoelastic model that it can predict undershoots in the transient shear viscosity.
Furthermore, note that both the intensity and the position of this undershoot are controlled by the parameters of the model. Such a behavior is not noticed for the first normal stress coefficient. Overall, the same trend is noticed also for the second normal stress coefficient (Figure 4) whereas, contrary to the other two viscometric functions, the transient curves go over (instead of below) the LVE predictions. Overall, it is a remarkable attribute of the new, single-mode viscoelastic model that it can predict undershoots in the transient shear viscosity.    Furthermore, note that both the intensity and the position of this undershoot are controlled by the parameters of the model. Such a behavior is not noticed for the first normal stress coefficient. Overall, the same trend is noticed also for the second normal stress coefficient (Figure 4) whereas, contrary to the other two viscometric functions, the transient curves go over (instead of below) the LVE predictions. Overall, it is a remarkable attribute of the new, single-mode viscoelastic model that it can predict undershoots in the transient shear viscosity.    Furthermore, note that both the intensity and the position of this undershoot are controlled by the parameters of the model. Such a behavior is not noticed for the first normal stress coefficient. Overall, the same trend is noticed also for the second normal stress coefficient (Figure 4) whereas, contrary to the other two viscometric functions, the transient curves go over (instead of below) the LVE predictions. Overall, it is a remarkable attribute of the new, single-mode viscoelastic model that it can predict undershoots in the transient shear viscosity.

Comparison with the Experimental Data Presented in
Stephanou et al. Figure 5 shows how the new model can fit the experimental data for the spectra of the storage and loss moduli discussed in Ref. [17] in conjunction with the earlier version of this model; the model predictions for both model variants are those presented in Equation (19). To carry out the fitting, we first need to specify two of the model parameters: the reptation time τ d , and either the entanglement modulus G e or the zero-rate viscosity η 0 . The values that were found to describe the data best are τ d = 11 s and η 0 = 68 kPa.s. The fitting is particularly good for small frequencies, which is the expected behavior. This is particularly pleasing given the use of only one mode in the model and the neglect of fast (Rouse) mode contributions. The discrepancy between the model and the experimental measurements at larger frequencies can be attributed, as is well known, to the neglect of some other equilibrium relaxation mechanisms, such as CLFs [38,39].

Comparison with the Experimental Data Presented in
Stephanou et al. Figure 5 shows how the new model can fit the experimental data for the spectra of the storage and loss moduli discussed in Ref. [17] in conjunction with the earlier version of this model; the model predictions for both model variants are those presented in Equation (19). To carry out the fitting, we first need to specify two of the model parameters: the reptation time τd, and either the entanglement modulus Ge or the zero-rate viscosity η0. The values that were found to describe the data best are τd = 11 s and η0 = 68 kPa.s. The fitting is particularly good for small frequencies, which is the expected behavior. This is particularly pleasing given the use of only one mode in the model and the neglect of fast (Rouse) mode contributions. The discrepancy between the model and the experimental measurements at larger frequencies can be attributed, as is well known, to the neglect of some other equilibrium relaxation mechanisms, such as CLFs [38,39]. Next, in Figure 6, we present the fittings for the steady-state values of the shear viscosity and of the two normal stress coefficients, whereas in Figures 7 and 8, we present the corresponding fittings in the case of the start-up of shear flow. There, the Rouse time was calculated through τR = τd/3Z and, since the entanglement molecular weight of the PS in the solution is equal to 64.6 kDa (see Table S-1 of Ref. [17]), we take be = 260. Then, and in order to choose the values of the rest of the model parameter, we considered the following: (a) the parameter α effectively controls the transient second normal stress coefficient, and (b) the undershoot observed at large shear rates in the transient shear viscosity is highly sensitive to the value of the non-affine parameter ξ (see Figure 7b). Thus, we came up with the following values: α = 0.4, βccr = 0.06, and ξ = 0.03. The graphs show that the model offers an exceptionally good description of all three steady-state viscometric functions. In addition, it can predict quite accurately the growth of the shear viscosity and of the first normal stress coefficient to flow start-up; even the undershoot time position, tu, in the shear viscosity at large shear rates is predicted quite accurately (Figure 7b). The description of the relative undershoot, du, (undershoot depth divided by the steady-state shear viscosity value), on the other hand, is less satisfactory ( Figure  7c). Finally, the growth of the first normal stress coefficient in the start-up of shear flow is noted to be adequately predicted (Figure 8a), whereas for the growth of the second normal stress coefficient we note that the prediction comes with an offset, even in the LVE envelope (Figure 8b). We hypothesize that this is due to the consideration of only one mode in the present analysis. Next, in Figure 6, we present the fittings for the steady-state values of the shear viscosity and of the two normal stress coefficients, whereas in Figures 7 and 8, we present the corresponding fittings in the case of the start-up of shear flow. There, the Rouse time was calculated through τ R = τ d /3Z and, since the entanglement molecular weight of the PS in the solution is equal to 64.6 kDa (see Table S-1 of Ref. [17]), we take b e = 260. Then, and in order to choose the values of the rest of the model parameter, we considered the following: (a) the parameter α effectively controls the transient second normal stress coefficient, and (b) the undershoot observed at large shear rates in the transient shear viscosity is highly sensitive to the value of the non-affine parameter ξ (see Figure 7b). Thus, we came up with the following values: α = 0.4, β ccr = 0.06, and ξ = 0.03. The graphs show that the model offers an exceptionally good description of all three steady-state viscometric functions. In addition, it can predict quite accurately the growth of the shear viscosity and of the first normal stress coefficient to flow start-up; even the undershoot time position, t u , in the shear viscosity at large shear rates is predicted quite accurately (Figure 7b). The description of the relative undershoot, d u , (undershoot depth divided by the steady-state shear viscosity value), on the other hand, is less satisfactory (Figure 7c). Finally, the growth of the first normal stress coefficient in the start-up of shear flow is noted to be adequately predicted (Figure 8a), whereas for the growth of the second normal stress coefficient we note that the prediction comes with an offset, even in the LVE envelope (Figure 8b). We hypothesize that this is due to the consideration of only one mode in the present analysis.

Comparison with the Experimental Data of Costanzo et al.
We next use the new version of our model to describe the experimental data of Costanzo et al. [16] for the polystyrene melt PS185k and the polystyrene solution PS285k/2k-65, both of which have the same number of entanglements, Z = 13.9. For the reptation time, we consider the values depicted as τ m in the Table 2 of Costanzo et al. [16] (5.10 s for the melt, and 5.32 s for the solution). And the same for the Rouse time and the zero-rate viscosity: τ R = 0.271 s and η 0 = 2.28 × 10 5 Pa.s for the PS185k melt, and τ R = 0.319 s and η 0 = 1.04 × 10 5 Pa.s for the PS285k/2k-65 solution. We then note in Figure 9 the overall agreement between the model predictions and the rheological data for both the melt and the solution despite the use of only one mode; however, we still note deviations at large frequencies as before. Finally, as mentioned before, for the melt b e = b e = 54, whereas for the solution M e = 83 (since 20.5 kDa as calculated through M e solution = M e melt /φ with φ being the polymer volume fraction [16]).
Moreover, for the melt, and since no undershoots are observed at large shear rates in the transient shear viscosity, we consider ξ = 0. For the remaining parameters, both for the melt and the solution, we follow the same strategy as above, the only difference being that now we use the model to also fit the steady-state and transient elongational viscosities for the two systems. Overall, the following model parameters are used: α = 0.15, β ccr = 0.2, and ξ = 0 for the melt; and α = 0.47, β ccr = 0.001, and ξ = 0.0015 for the solution. Again, the model can describe very well the steady-state ( Figure 10) and transient ( Figure 11) shear viscosity, and the steady-state ( Figure 12) and transient ( Figure 13) elongational viscosity for both the PS185k melt, and the PS285k/2k-65 solution. Note also that in the case of the elongational flow, we have used ξ = 0 both for the melt and the solution, since no chain tumbling occurs in extensional flows due to their irrotational character. Again, the description of the growth of the shear viscosity is quite gratifying, especially of the position of the undershoot at large shear rates in the case of the polymer solution ( Figure 11c); the relative undershoot, on the other hand, is described again less accurately (Figure 11d). Overall, the revised model, even with the use of only a single mode, can quantitatively capture the experimental trends both for entangled polymer melts and concentrated polymer solutions.
Materials 2020, 13, x FOR PEER REVIEW 12 of 18 the same number of entanglements, Z = 13.9. For the reptation time, we consider the values depicted as τm in the Table 2 of Costanzo et al. [16] (5.10 s for the melt, and 5.32 s for the solution). And the same for the Rouse time and the zero-rate viscosity: τR = 0.271 s and η0 = 2.28 × 10 5 Pa.s for the PS185k melt, and τR = 0.319 s and η0 = 1.04 × 10 5 Pa.s for the PS285k/2k-65 solution. We then note in Figure 9 the overall agreement between the model predictions and the rheological data for both the melt and the solution despite the use of only one mode; however, we still note deviations at large frequencies as before. Finally, as mentioned before, for the melt be = be = 54, whereas for the solution Me = 83 (since 20.5 kDa as calculated through Me solution = Me melt /ϕ with ϕ being the polymer volume fraction [16]). Moreover, for the melt, and since no undershoots are observed at large shear rates in the transient shear viscosity, we consider ξ = 0. For the remaining parameters, both for the melt and the solution, we follow the same strategy as above, the only difference being that now we use the model to also fit the steady-state and transient elongational viscosities for the two systems. Overall, the following model parameters are used: α = 0.15, βccr = 0.2, and ξ = 0 for the melt; and α = 0.47, βccr = 0.001, and ξ = 0.0015 for the solution. Again, the model can describe very well the steady-state ( Figure 10) and transient ( Figure 11) shear viscosity, and the steady-state ( Figure 12) and transient ( Figure 13) elongational viscosity for both the PS185k melt, and the PS285k/2k-65 solution. Note also that in the case of the elongational flow, we have used ξ = 0 both for the melt and the solution, since no chain tumbling occurs in extensional flows due to their irrotational character. Again, the description of the growth of the shear viscosity is quite gratifying, especially of the position of the undershoot at large shear rates in the case of the polymer solution ( Figure 11c); the relative undershoot, on the other hand, is described again less accurately (Figure 11d). Overall, the revised model, even with the use of only a single mode, can quantitatively capture the experimental trends both for entangled polymer melts and concentrated polymer solutions.   We then note in Figure 9 the overall agreement between the model predictions and the rheological data for both the melt and the solution despite the use of only one mode; however, we still note deviations at large frequencies as before. Finally, as mentioned before, for the melt be = be = 54, whereas for the solution Me = 83 (since 20.5 kDa as calculated through Me solution = Me melt /ϕ with ϕ being the polymer volume fraction [16]). Moreover, for the melt, and since no undershoots are observed at large shear rates in the transient shear viscosity, we consider ξ = 0. For the remaining parameters, both for the melt and the solution, we follow the same strategy as above, the only difference being that now we use the model to also fit the steady-state and transient elongational viscosities for the two systems. Overall, the following model parameters are used: α = 0.15, βccr = 0.2, and ξ = 0 for the melt; and α = 0.47, βccr = 0.001, and ξ = 0.0015 for the solution. Again, the model can describe very well the steady-state ( Figure 10) and transient ( Figure 11) shear viscosity, and the steady-state ( Figure 12) and transient ( Figure 13) elongational viscosity for both the PS185k melt, and the PS285k/2k-65 solution. Note also that in the case of the elongational flow, we have used ξ = 0 both for the melt and the solution, since no chain tumbling occurs in extensional flows due to their irrotational character. Again, the description of the growth of the shear viscosity is quite gratifying, especially of the position of the undershoot at large shear rates in the case of the polymer solution ( Figure 11c); the relative undershoot, on the other hand, is described again less accurately ( Figure 11d). Overall, the revised model, even with the use of only a single mode, can quantitatively capture the experimental trends both for entangled polymer melts and concentrated polymer solutions.      Figure  11a,b depicts the LVE envelope, Equation (15a). Same parameter values as in Figure 10.   Figure 9).  Figure  11a,b depicts the LVE envelope, Equation (15a). Same parameter values as in Figure 10.  Figure 11a,b depicts the LVE envelope, Equation (15a). Same parameter values as in Figure 10.
Materials 2020, 13, x FOR PEER REVIEW 14 of 18 Figure 12. Comparison of the model predictions with the experimental measurements of Costanzo et al. [16] for the steady uniaxial extensional viscosity. Same parameter values as in Figure 11 but in all cases ξ = 0.

Conclusions
Some time ago, we presented a reformulation of the Marrucci-Ianniruberto constitutive equation [14] for the rheology of entangled polymer melts and concentrated solutions in the context of nonequilibrium thermodynamics [25]. However, rather large values of the CCR parameter had to be invoked for the model not to violate the second law of thermodynamics. Here, we have presented an appropriate modification that avoids the splitting of the evolution equation for the conformation tensor into an orientation and a stretching part; as a result, thermodynamic admissibility allows for a friendlier range of values of the CCR parameter to be used (they must only be non-negative). In addition, and in view of recent experimental evidence for the appearance of undershoots in the transient stress (following the overshoot) at high shear rates whose origin may be traced back to molecular tumbling at high shear rates, the new model contains additional terms accounting for nonaffine chain deformation through a slip parameter.
The capability of the new model to describe very satisfactorily rheological data was demonstrated by using it to fit available experimental data for the transient and steady-state shear and elongational rheology of entangled PS melts and solutions. Remarkable success was observed for all viscometric functions checked, considering in particular the use of only a single mode. Given its simplicity and thermodynamic rigor (proven thermodynamic admissibility), this is a remarkable feature of the new model. At the same time, it can also predict quite accurately the elongational rheological data. It can thus be used as an excellent alternative to other successful constitutive models, Figure 12. Comparison of the model predictions with the experimental measurements of Costanzo et al. [16] for the steady uniaxial extensional viscosity. Same parameter values as in Figure 11 but in all cases ξ= 0.
Materials 2020, 13, x FOR PEER REVIEW 14 of 18 Figure 12. Comparison of the model predictions with the experimental measurements of Costanzo et al. [16] for the steady uniaxial extensional viscosity. Same parameter values as in Figure 11 but in all cases ξ = 0.

Conclusions
Some time ago, we presented a reformulation of the Marrucci-Ianniruberto constitutive equation [14] for the rheology of entangled polymer melts and concentrated solutions in the context of nonequilibrium thermodynamics [25]. However, rather large values of the CCR parameter had to be invoked for the model not to violate the second law of thermodynamics. Here, we have presented an appropriate modification that avoids the splitting of the evolution equation for the conformation tensor into an orientation and a stretching part; as a result, thermodynamic admissibility allows for a friendlier range of values of the CCR parameter to be used (they must only be non-negative). In addition, and in view of recent experimental evidence for the appearance of undershoots in the transient stress (following the overshoot) at high shear rates whose origin may be traced back to molecular tumbling at high shear rates, the new model contains additional terms accounting for nonaffine chain deformation through a slip parameter.
The capability of the new model to describe very satisfactorily rheological data was demonstrated by using it to fit available experimental data for the transient and steady-state shear and elongational rheology of entangled PS melts and solutions. Remarkable success was observed for all viscometric functions checked, considering in particular the use of only a single mode. Given its simplicity and thermodynamic rigor (proven thermodynamic admissibility), this is a remarkable feature of the new model. At the same time, it can also predict quite accurately the elongational rheological data. It can thus be used as an excellent alternative to other successful constitutive models, Figure 13. Comparison of the model predictions with the experimental measurements of Costanzo et al. [16] for the extensional stress growth coefficient as a function of time for several stretch rates. The thick orange line depicts the LVE envelope, Equation (17). Same parameter values as in Figure 12

Conclusions
Some time ago, we presented a reformulation of the Marrucci-Ianniruberto constitutive equation [14] for the rheology of entangled polymer melts and concentrated solutions in the context of non-equilibrium thermodynamics [25]. However, rather large values of the CCR parameter had to be invoked for the model not to violate the second law of thermodynamics. Here, we have presented an appropriate modification that avoids the splitting of the evolution equation for the conformation tensor into an orientation and a stretching part; as a result, thermodynamic admissibility allows for a friendlier range of values of the CCR parameter to be used (they must only be non-negative). In addition, and in view of recent experimental evidence for the appearance of undershoots in the transient stress (following the overshoot) at high shear rates whose origin may be traced back to molecular tumbling at high shear rates, the new model contains additional terms accounting for non-affine chain deformation through a slip parameter.
The capability of the new model to describe very satisfactorily rheological data was demonstrated by using it to fit available experimental data for the transient and steady-state shear and elongational rheology of entangled PS melts and solutions. Remarkable success was observed for all viscometric functions checked, considering in particular the use of only a single mode. Given its simplicity and thermodynamic rigor (proven thermodynamic admissibility), this is a remarkable feature of the new model. At the same time, it can also predict quite accurately the elongational rheological data. It can thus be used as an excellent alternative to other successful constitutive models, such as the Rolie-Poly [26] for characterizing the rheological response of polymeric fluids or for carrying out large scale numerical calculations of complex viscoelastic flows.
Of course, the new version of the model is far from being considered complete. There is a few more features and mechanisms of relevance to the flow of entangled polymer fluids that must be accounted for in the corresponding constitutive and transport equations. First, we certainly need to include in the model some additional relaxation mechanisms that were neglected in the present treatment (such as CLFs), in combination with a more rigorous treatment of CR effects. Second, we mention the need to implement a configuration dependent friction coefficient as has been proposed by Mead et al. [40], as well as the proper extension to a multi-mode formulation. Third, and given the wealth of information provided by recent NEMD studies [6][7][8]15] at the level of individual chains and their conformations, it would be highly desirable to check how well the model can describe the conformational properties of entangled polymer melts under flow (and not just the viscometric functions). For example, Sefiddashti et al. [9] have already carried out such a first MD test of several constitutive models. In turn, such a direct comparison would help fully parameterize the model on the outcome of these direct NEMD simulations. It is by bringing together simulations and theory that one can achieve a more complete description of polymer flows, which can then be encoded in the form of friendly and easy-to-use, yet powerful and accurate constitutive laws for the subsequent use in numerical calculations of viscoelastic flows through geometries of practical relevance. Finally, it would be of interest to solve the model in the case of a large amplitude oscillatory shear (LAOS) flow to check how well it can address the nonlinear viscoelastic oscillatory response of entangled polymers [41].