On the Use of a Non-Constant Non-Afﬁne or Slip Parameter in Polymer Rheology Constitutive Modeling

P.S


Introduction
A fundamental understanding of the rheological and microstructural behavior of polymeric fluids under flow is essential in practical polymer processing operations [1][2][3]. Polymers are fluids that exhibit a non-Newtonian rheological character, which stems from their internal microstructure. The dynamical behavior of the microstructure is usually challenging to obtain experimentally, leaving ample space for computational simulations. Nowadays, due to the growth of computational power and the advent of sophisticated computational algorithms, computationally executed experiments offer the only alternative to experiments that cannot be conducted physically [4], such as quantifying the key role of threading events in linear-ring blends [5][6][7][8]. This is particularly true for polymeric materials, even the simplest of all, i.e., linear polyethylene (PE) chains, since one must span a spectrum of orders of magnitude in both time and space [4]. Today, simulations, particularly atomistic ones (i.e., simulations in which individual atoms of ensembles of molecules are tracked in phase space), offer the best avenue to validate theories without hypotheses. Two such examples are the Rouse theory for cyclic polymers (rings) [9] and the tube notion in entangled polymeric systems (i.e., systems whose molecular weights exceed the entanglement molecular weight, and whose dynamics is dictated by a slithering or snake-like motion within a confining mean-field tube produced by the surrounding chain molecules) [10]. However, similar success has not been accomplished yet in nonequilibrium systems, i.e., systems that are perturbed away from equilibrium, such as those under flow, which is more complicated to simulate [4].
Still, non-equilibrium molecular dynamics simulations (NEMD) have provided invaluable insight into the dynamical behavior of polymer chains under flow. Experimentally, it has been possible to directly visualize individual polymer chains when subjected to flow. For example, we mention the work of Smith et al. [11] and LeDuc et al. [12], who used video microscopy to study the dynamics of individual, tagged chains under shear in dilute ξ = tr(C)/ 1 3 R 2 eq β + tr(C) , where C is the dimensional conformation tensor, tr(C) is the trace of C, and R 2 eq is the average squared end-to-end polymer distance at equilibrium (see the next section) [33]. However, as the polymer aspect ratio increases due to deformation, one would expect that ξ should also be time dependent. Such a consideration has not been cogitated in the past, with the sole exception being the work of Beris et al. [34], that introduced a phenomenological kinetic equation for the slip parameter. They also introduced a limiting value for the non-affine parameter at high shear rates. However, in their model, the slip parameter was only coupled to the shear rate and not the structure itself. To the best of our knowledge, there has not been any other work wherein a (shear-rateor time-) dependent slip parameter was considered.
Although other constitutive models have considered a variety of mechanisms to accommodate the tumbling of polymer chains under shear, such as the use of a tumbling function by Costanzo et al. [24] or the explicit consideration of rotational diffusion by Stephanou et al. [23], the necessity of having a variable slip parameter remains. In this work, we generalize the constitutive model for unentangled polymer melts of Stephanou et al. [29] to accommodate a varying slip parameter by including an additional evolution equation for ξ, following our recent work on the use of a scalar structural variable [35].
This paper is structured as follows: in Section 2, the new model is introduced, whereas Section 3 presents the details concerning the simulated system and technical information regarding the MD and NEMD performed. Then, in Section 4, we offer the model predictions along with a comparison with the accumulated simulation data. The paper concludes with Section 5, where we elaborate on the significance of our work and discuss future plans.

Model Modification
Following previous work [29,30], we define the conformation tensor as the second moment of the distribution function Ψ(R, r, t) for the chain end-to-end vector R, with its center-of-mass at position r, i.e., C(r, t) = RR (r, t) = RRΨ(R, r, t)dR with the brackets denoting a configurational average. The evolution equation for the dimensionless conformation tensor c(r, t) = 3C(r, t)/ R 2 eq , as derived by Stephanou et al. [29], is: . where: .
denotes the GS or JS mixed derivative, with . γ ≡ ∇u + (∇u) T being the rate-of-deformation tensor (X T is the transpose of X), and ∇u is the velocity gradient tensor. Additionally, I is the unit tensor, α is the anisotropic Giesekus parameter, is the effective spring constant, accounting for FENE effects, with b the FENE parameter, and the function: with 0 ≤ B(c) ≤ 1, which ensures that the entropy density remains bounded even at high deformation rates [36], with a 0 (c) the dimensionless unbounded free energy [29] given as: and the relaxation, or Rouse, time: with ε the Phan-Thien Tanner parameter, and τ R,eq the equilibrium Rouse time (note that the Rouse time in Stephanou et al. [29] was defined as λ(tr(c))_. Note that other expressions could also be employed, such as the extended White-Metzner expression [30,35,37]. Finally, the corresponding expression relating the stress tensor with the conformation tensor is given as: where G is the elastic modulus.

383
We now proceed to modify the Stephanou et al. [29] model by considering a variable slip parameter. To this end, we simply propose an expression for the evolution equation of ξ bearing in mind our previous work concerning the scalar structural variable λ [35]: where κ = (∇u) T , τ ξ is the characteristic time for the increase in the slip parameter, and ξ 0 is the upper bound of the slip parameter at high shear rates. The first term in Equation (3) is a relaxation term returning the slip parameter to its equilibrium, null, value when the flow is ceased. In contrast, whereas the second is a term that increases the slip parameter as a result of the applied flow. In the following, we define γ = τ ξ /τ R,eq .

Asymptotic Behavior of the Model for Steady-State and Transient Shear Flow
Here, we analyze the asymptotic behavior of the revised model in the limits of low deformation rates for both steady-state and transient simple shear flow (SSF), described by the kinematics u = ε is the elongation rate (x is the flow direction, y is the velocity gradient direction, and z is the neutral direction). The material functions to analyze are the shear viscosity, η ≡ σ yx / . γ, and the two normal stress coefficients, respectively, in the case of shear, and the extensional viscosity, ε, in the case of uniaxial elongation. By expanding the conformation tensor and the slip parameter up to second-order terms in the dimensionless shear rate, Wi = . γτ R,eq , we arrive at the following expressions for the conformation tensor and the slip parameter (when considering a finite value of γ): and for the zero-shear-rate shear viscosity and normal stress coefficients: These are the same as the ones presented by Stephanou et al. [29] by considering ξ = 0 in their Equations (41) and (42). On the other hand, at large shear rates, the slip parameter will approach the upper bound value ξ 0 ; see the next section.
Upon inception of the simple shear flow, the explicit solutions for the time-dependent viscometric functions in the linear viscoelastic (LVE) limit, following the methodology of Stephanou et al. [38], are given as: In steady-state uniaxial elongation, by expanding the conformation tensor and the slip parameter up to first-order terms in the dimensionless elongation rate, Wi = . ετ R,eq , we obtain: and for the zero-elongation-rate elongation viscosity: meaning that Trouton's law holds. Upon inception of elongation flow, the explicit solutions for the time-dependent elongation viscosity in the linear viscoelastic limit, again following the methodology of Stephanou et al. [38], is given as:

Molecular Model and System Studied
In this work, we conducted equilibrium molecular dynamics (MD) simulations, and NEMD simulations of PE oligomer C 48 H 98 melts. Equilibrium MD simulations were performed in the NPT ensemble to fully relax the initial PE configurations at temperature T = 450 K and pressure P = 1 atm using the united-atom potential model of Siepmann et al. [39]. The simulations were carried out using the LAMMPS simulation engine [40], employing the Nosé-Hoover thermostat [41,42] and the Parrinello-Rahman barostat [43] to preserve the temperature and pressure, respectively, at their prescribed values. Subsequently, several fully relaxed configurations from the MD runs were selected as input for the NEMD simulations under shear. The NEMD runs were performed again with LAMMPS in the NVT ensemble at T = 450 K, using the SLLOD algorithm [44], together with the Nosé-Hoover thermostat to control the temperature. The microscopic set of equations of motion was integrated numerically using the reversible Reference System Propagator Algorithm (r-RESPA) [45], with 2 different time steps: (a) a large one (dt = 4 fs) for the integration of the slowest varying forces arising from non-bonded interactions at long interatomic distances, and (b) a small one (dt = 0.5 fs) for the integration of the fast-varying forces corresponding to the bonded (i.e., bonds, angles, and dihedrals) interactions.
All simulations were conducted using large cells containing 16,000 chain molecules of C 48 H 98 and were subjected to periodic boundary conditions in all three directions (x, y, and z). In the course of the NEMD simulations, the xand y-directions were selected as the flow and the shear gradient directions, respectively, whereas z was the neutral direction. The simulation cell had dimensions (462 Å) × (231 Å) × (231 Å) along the x-, y-, and z-directions. The cell was purposefully enlarged in the (x-) flow direction to ensure minimal system size effects, particularly with the NEMD runs at high shear rates, where the polymer chains tend to stretch and orient towards the flow direction. To this end, for the C 48 H 98 chains, the equilibrium root-mean-square of the chain end-to-end vector R 2 eq and the theoretical maximum chain extension of |R| max were calculated to be equal to R 2 eq = 27.2 ± 0.12 and |R| max = 63.1 , respectively. Compared to the simulation cell dimensions, the maximum chain length |R| max is 7.3 times shorter than the dimension in the x-direction and 3.6 times shorter than the dimension in the y-direction. Thus, we can safely expect that the simulation cell is sufficiently large to ensure the absence of system size effects due to chain alignment in the flow direction or tumbling motion in the shear gradient direction.
In the course of the equilibrium MD simulations, the equilibrium orientational relaxation time τ R,eq of the simulated C 48 H 98 PE chains at T = 450 K and P = 1 atm was found to be equal to τ R,eq = 0.6 ± 0.01 ns, as estimated by integrating the stretched-exponential curve [46] over time, describing the time autocorrelation function of the chain end-to-end unit vector. The MD simulations were conducted for a total of 6 ns, which is 10 times larger than the chain relaxation time to ensure that the PE chains were fully equilibrated. The NEMD simulations were executed over a broad range of shear rates spanning the range from the linear up to the highly non-linear viscoelastic regime, corresponding to Weissenberg numbers (with τ R,eq = 0.6 ns) in the interval [0.1, 285].

Model Predictions in Steady-State Shear Flow
In this section, we present the predictions of the new model in the case of steadystate shear flow. The results were obtained numerically by solving the constitutive model, Equation (1), under steady-state conditions using MATLAB [47] and then calculating the stress tensor using Equation (2).
In Figure 1, we depict the model prediction for the slip parameter as a function of Wi and its dependence on the parameters ξ 0 , γ, and ε while keeping the constant α = 0.4, b = 50 (panel (a)), and (b) the parameters α, ε, and b while maintaining the constant ξ 0 = 0.05, γ = 1 (panel (b)). As noted in Figure 1a, at large values of γ, the slip parameter steeply reaches its limiting value, which is unaltered as the shear rate increases further. This situation resembles the previous version of the model [29], where the slip parameter was considered a constant. However, as γ is reduced to unity, the slip parameter increases with the shear rate and reaches its limiting value at about Wi = 10. As the parameter ξ 0 increases, the curve shifts upwards, whereas by increasing ε, we note a slight shift upwards at higher shear rates. This is because when ε = 0 the slip parameter does not fully reach its limiting value since c xy delays its reduction with the shear rate, as shown in Figure 2b. On the other hand, the slip parameter seems insensitive to the precise values of α and b (panel (b)). Note that analytical expression Equation (4e) is accurate in all cases until about Wi = 1.    Next, in Figure 2, we depict the model predictions for the conformation tensor as a function of Wi whilst keeping = 0.4, = 50 constant, whereas, in Figure 3, we depict the same comparison whilst keeping = 0.05, = 1 constant. Irrespective of the values of the parameters (see both Figures 2 and 3), we note that cxx is reported to increase from its equilibrium value after about ≈Wi = 0.3 and eventually reach its limiting value, which, however, differs from the value of the parameter b. On the other hand, cxy is noted to Next, in Figure 2, we depict the model predictions for the conformation tensor as a function of Wi whilst keeping α = 0.4, b = 50 constant, whereas, in Figure 3, we depict the same comparison whilst keeping ξ 0 = 0.05, γ = 1 constant. Irrespective of the values of the parameters (see both Figures 2 and 3), we note that c xx is reported to increase from its equilibrium value after about ≈Wi = 0.3 and eventually reach its limiting value, which, however, differs from the value of the parameter b. On the other hand, c xy is noted to increase linearly with the shear rate at low shear rates, as dictated by Equation (4b), reaching a maximum value, and then decreasing inversely proportional to the shear rate. Finally, the two remaining diagonal elements of the conformation tensor in the shear gradient direction and neutral direction (c yy and c zz , respectively), are observed to be mirrors of the noted behavior of c xx : they initially decrease from their equilibrium value, eventually reaching a finite asymptotic value at high shear rates. We note, in Figure 2, that the value of γ plays only a very modest role for all elements, since the value of ξ 0 is small; note that one would expect such small values to be used since larger values would lead to very intense oscillations in the time-dependent material functions (see the next section). On the other hand, by increasing the parameter ξ 0 , we note the predictions at low shear rates to be insensitive. Still, the limiting asymptotic values at high shear rates are pointed out to decrease for c xx and increase for both c yy and c zz , whereas the c xy curve shifts to lower shear rates at higher shear rates whilst keeping the power-law unaffected. This is a direct result of allowing tumbling to occur sooner, since the slip parameter is larger, thus refraining the flow to deform the chain further. Finally, when we increase the value of ε while keeping the same ξ 0 and γ values (Figure 3), we note that the asymptotic values of the diagonal elements remain the same, but the curves are shifted to the right, a direct result of the steeper decrease in the relaxation time, cf. Equation (1g). As the anisotropic (or Giesekus) parameter, α, is increased, we observe that the curve of c xx shifts downwards and the one of c zz shifts upwards (see Figure 3). In constrast, when the FENE parameter is increased (from 20 to 50, meaning that the chain is now longer or its molecular weight is larger), we observe the reverse behavior, which is the expected outcome. However, note that during these parameter value changes, the predictions of the other elements, c xy (panel (b)) and c yy (panel (c)) are only modestly affected. As in Figure 2, the value of ε controls the rate at which the asymptotic values, in the case of the diagonal elements, are reached, whereas the c xy curve shifts rightwards. other hand, by increasing the parameter , we note the predictions at low shear rates to be insensitive. Still, the limiting asymptotic values at high shear rates are pointed out to decrease for cxx and increase for both cyy and czz, whereas the cxy curve shifts to lower shear rates at higher shear rates whilst keeping the power-law unaffected. This is a direct result of allowing tumbling to occur sooner, since the slip parameter is larger, thus refraining the flow to deform the chain further. Finally, when we increase the value of while keeping the same and values (Figure 3), we note that the asymptotic values of the diagonal elements remain the same, but the curves are shifted to the right, a direct result of the steeper decrease in the relaxation time, cf. Equation (1g). As the anisotropic (or Giesekus) parameter, , is increased, we observe that the curve of cxx shifts downwards and the one of czz shifts upwards (see Figure 3). In constrast,when the FENE parameter is increased (from 20 to 50, meaning that the chain is now longer or its molecular weight is larger), we observe the reverse behavior, which is the expected outcome. However, note that during these parameter value changes, the predictions of the other elements, cxy (panel (b)) and cyy (panel (c)) are only modestly affected. As in Figure 2, the value of controls the rate at which the asymptotic values, in the case of the diagonal elements, are reached, whereas the cxy curve shifts rightwards. In Figures 4 and 5, we depict the same comparison as in Figures 2 and 3, respectively, but for the three dimensionless material functions: the shear viscosity, η (panel (a)), and the first, Ψ 1 (panel (b)), and negative second, −Ψ 2 (panel (c)), normal stress coefficients. We again note, in Figure 4, that the predictions are insensitive to the value of the slip parameter due to the small value of ξ 0 . By increasing the value of ξ 0 , both η and Ψ 1 are noted to shift to lower shear rates, which is much more intense for the shear viscosity, whereas when increasing ε, we note that the shear viscosity curve shifts to higher shear rates, and Ψ 1 remains almost unaltered at large shear rates but shifts rightwards at intermediate Wi.
On the contrary, −Ψ 2 is noted to be almost completely unaffected. Finally, Figure 5 shows that both η and Ψ 1 are invariant to changes in the values of α and b. On the other hand, the zero-shear-rate negative second normal stress coefficient increases as the Giesekus parameter increases, cf. Equation (5c), but again remains invariant at higher shear rates.
We again note, in Figure 4, that the predictions are insensitive to the value of the slip parameter due to the small value of . By increasing the value of , both and Ψ are noted to shift to lower shear rates, which is much more intense for the shear viscosity, whereas when increasing , we note that the shear viscosity curve shifts to higher shear rates, and Ψ remains almost unaltered at large shear rates but shifts rightwards at intermediate Wi. On the contrary, −Ψ is noted to be almost completely unaffected. Finally, Figure 5 shows that both and Ψ are invariant to changes in the values of and b. On the other hand, the zero-shear-rate negative second normal stress coefficient increases as the Giesekus parameter increases, cf. Equation (5c), but again remains invariant at higher shear rates.

Model Predictions in Start-Up Shear Flow
In this section, we present the predictions of the new model in the case of start-up shear flow. The results were obtained numerically by solving the system of differential equations, Equations (1), using MATLAB [47] and then calculating the stress tensor, Equation (2). Here, we consider only the case with γ = 1, since the case of large γ is identical to the predictions presented in Stephanou et al. [29].
In Figure 6, we depict the growth of the slip parameter under start-up shear flow at two different values of the dimensionless shear rate equal to Wi = 1 and 10 and various values of the model parameters. We note that irrespective of the model parameter, it increased exponentially, reaching its steady-state value. A modest dumping behavior is noted wherein the time-dependent prediction is observed to oscillate around the steadystate value before reaching it. The overall behavior of ξ is seen not to depend heavily on the values of the parameters α, ε, and b.

Model Predictions in Start-Up Shear Flow
In this section, we present the predictions of the new model in the case of start-up shear flow. The results were obtained numerically by solving the system of differential equations, Equations (1), using MATLAB [47] and then calculating the stress tensor, Equation (2). Here, we consider only the case with = 1, since the case of large is identical to the predictions presented in Stephanou et al. [29].
Ιn Figure 6, we depict the growth of the slip parameter under start-up shear flow at two different values of the dimensionless shear rate equal to Wi = 1 and 10 and various values of the model parameters. We note that irrespective of the model parameter, it increased exponentially, reaching its steady-state value. A modest dumping behavior is noted wherein the time-dependent prediction is observed to oscillate around the steadystate value before reaching it. The overall behavior of is seen not to depend heavily on the values of the parameters , , and b.
Next, in Figures 7 and 8, we depict the growth of the conformation tensor under startup shear flow at two different values of the dimensionless shear rate equal to Wi = 1 and 10. We observe that the predictions for cxx, cxy, and cyy also present a dumping behavior, which is much more intense relative to the one noted in the slip parameter ( Figure 6). It intensifies as the parameter increases and becomes less intense when the parameter increases. Such a behavior is obtained only when > 0. On the other hand, czz seems not to present such a behavior, although a small saddle point is noted when = 0. The timedependent behavior of the conformation tensor is seen in Figure 8 to be insensitive to the value of the parameters and .  γτ R,eq , and dependence on: (a) the parameters ξ 0 and ε for α = 0.4, b = 50; and (b) the parameters α, ε, and b for ξ 0 = 0.05 and γ = 1. In all cases, γ = 1.
Next, in Figures 7 and 8, we depict the growth of the conformation tensor under startup shear flow at two different values of the dimensionless shear rate equal to Wi = 1 and 10. We observe that the predictions for c xx , c xy , and c yy also present a dumping behavior, which is much more intense relative to the one noted in the slip parameter ( Figure 6). It intensifies as the ξ 0 parameter increases and becomes less intense when the ε parameter increases. Such a behavior is obtained only when ξ 0 > 0. On the other hand, c zz seems not to present such a behavior, although a small saddle point is noted when ε = 0. The time-dependent behavior of the conformation tensor is seen in Figure 8 to be insensitive to the value of the parameters α and ε.   γτ R,eq , and dependence on the parameters ξ 0 , γ, and ε for α = 0.4, b = 50.  Note that the dotted dark yellow and light grey lines in Figure 10(c) in each panel depict the LVE envelope given by Equations (6). We again note that the predictions for the shear viscosity present a dumping behavior for both shear rates when ≠ 0, which intensifies as the shear rates increases; such a behavior can also be noted when is a constant ( ≫ 1) [29,38]. It also intensifies when increases and is less intense when increases (Figure 9(a)), following the behavior noted in the case of the conformation tensor. A similar behavior is also noted in Ψ ( ) (panel (b)), although the dumping behavior seems to be less intense. On the other hand, −Ψ ( ) seems not to exhibit such a behavior (Figure 9(c)). As either the parameter or the FENE parameter b is increased, the time-dependent behavior is unaffected (Figure (10)). However, when increasing , the overshoot is noted to decrease and shift to shorter times for both ( ) (panel (a)) and Ψ ( ) (panel (b)). On the other hand, as the parameter increases, the −Ψ ( ) curve is noted to go over its LVE envelope at short times, which is never the case for both ( ) and Ψ ( ). γτ R,eq , and dependence on the parameters α, ε, and b for of ξ 0 = 0.01 and γ = 1.
Finally, in Figures 9 and 10, we depict the growth of the three dimensionless materials functions: the shear viscosity, η + (t) (panel (a)), and the first, Ψ + 1 (t) (panel (b)), and negative second, −Ψ + 2 (t) (panel (c)), normal stress coefficients under start-up shear flow. Note that the dotted dark yellow and light grey lines in Figure 10c in each panel depict the LVE envelope given by Equations (6). We again note that the predictions for the shear viscosity present a dumping behavior for both shear rates when ξ 0 = 0, which intensifies as the shear rates increases; such a behavior can also be noted when ξ is a constant (γ 1) [29,38]. It also intensifies when ξ 0 increases and is less intense when ε increases (Figure 9a), following the behavior noted in the case of the conformation tensor. A similar behavior is also noted in Ψ + 1 (t) (panel (b)), although the dumping behavior seems to be less intense. On the other hand, −Ψ + 2 (t) seems not to exhibit such a behavior ( Figure 9c). As either the parameter α or the FENE parameter b is increased, the time-dependent behavior is unaffected ( Figure 10). However, when increasing ε, the overshoot is noted to decrease and shift to shorter times for both η + (t) (panel (a)) and Ψ + 1 (t) (panel (b)). On the other hand, as the parameter α increases, the −Ψ + 2 (t) curve is noted to go over its LVE envelope at short times, which is never the case for both η + (t) and Ψ + 1 (t).

Comparison with NEMD Simulation Data for an Unentangled PE Melt
In this section, we aim to compare the predictions of the revised model against rheological data obtained through newly accumulated, fully atomistic NEMD simulations of a short unentangled C48 PE melt over a broad spectrum of shear rates. We only analyze the steady-state data of these NEMD simulations. To fit the simulation data, following Stephanou et al. [29], we first identify the asymptotes of the diagonal elements of the conformation tensor in the limit of high shear rates. Additionally, the equilibrium relaxation (Rouse) time, as mentioned in Section 3, is equal to , = 0.6 ns, whereas the zeroshear-rate viscosity can be obtained from the shear viscosity NEMD data, which is equal to = 5 mPa. s. Note that as small shear rates, the NEMD data of the viscometric functions come with large error bars, and it is difficult, particularly for the normal stress coefficients, to accurately estimate their zero-shear-rate values. Then, the value of ≈ 0.104 is obtained by using the first equation of Equation (45) of Stephanou et al. [29]. Next, the value of the Giesekus parameter ≈ 0.2 can be obtained by fitting the Ψ NEMD data at low shear rates and using Equation (4c) (the value of Ψ , is easily calculated from Equation (4b)). Note that this value differs from the value ≈ 0.06 obtained using the second equation of Equation (45) of Stephanou et al. [29]; however, the fitting of Ψ is much improved when using the former value, and the comparison against the γτ R,eq , and dependence on the parameters α, ε, and b for of ξ 0 = 0.01 and γ = 1. Note that two LVE envelopes are provided for −Ψ + 2 (t), one for each α value.

Comparison with NEMD Simulation Data for an Unentangled PE Melt
In this section, we aim to compare the predictions of the revised model against rheological data obtained through newly accumulated, fully atomistic NEMD simulations of a short unentangled C 48 PE melt over a broad spectrum of shear rates. We only analyze the steady-state data of these NEMD simulations. To fit the simulation data, following Stephanou et al. [29], we first identify the asymptotes c ∞ ii of the diagonal elements of the conformation tensor in the limit of high shear rates. Additionally, the equilibrium relaxation (Rouse) time, as mentioned in Section 3, is equal to τ R,eq = 0.6 ns, whereas the zero-shearrate viscosity can be obtained from the shear viscosity NEMD data, which is equal to η 0 = 5 mPa.s. Note that as small shear rates, the NEMD data of the viscometric functions come with large error bars, and it is difficult, particularly for the normal stress coefficients, to accurately estimate their zero-shear-rate values. Then, the value of ξ 0 ≈ 0.104 is obtained by using the first equation of Equation (45) of Stephanou et al. [29]. Next, the value of the Giesekus parameter α ≈ 0.2 can be obtained by fitting the Ψ 2 NEMD data at low shear rates and using Equation (4c) (the value of Ψ 1,0 is easily calculated from Equation (4b)). Note that this value differs from the value α ≈ 0.06 obtained using the second equation of Equation (45) of Stephanou et al. [29]; however, the fitting of Ψ 2 is much improved when using the former value, and the comparison against the conformation tensor data is only mildly worsened. Next, the value of b eff =5.78 can be obtained from Equation (48) of Stephanou et al. [29], which differs from the value 3 2 / R 2 eq = 16.15. The remaining two parameters, ε and γ, can be obtained by simply fitting the NEMD data, since they do not affect the c ∞ ii values; we obtain ε = 0.4 and γ = 0.01. Figure 11 shows how well the new model can fit the simulation data for the c xx , c xy , c yy , and c zz elements of the dimensionless conformation tensor for the C 48 PE system in steady shear. We observe that the predictions of the revised model are in remarkable agreement with the NEMD extracted simulation results over the entire wide range of shear rates considered, especially for the diagonal elements. To quantify how well the model predicts the simulation data, we calculated the sum of the squares of the residuals (i.e., the residual between the simulation value and the one obtained by the model); this turns out to be ≈0.2474. The corresponding comparison for the material functions η, Ψ 1 , and −Ψ 2 is presented in Figure 12. Contrary to the very good agreement between the refined model predictions and the NEMD simulation data for the dimensionless conformation tensor, the comparison against the viscometric functions is less satisfactory. Deviations from the NEMD data are mainly observed at large shear rates in the case of the shear viscosity (panel (a)) and the second normal stress coefficient (panel (c)). As also mentioned by Stephanou et al. [29], this disaccord should be related to the postulated relation between the stress and conformation tensors, Equation (2), which stems from the assumption of purely entropic elasticity [48]. As such, a more accurate expression for the free energy needs to be invoked in the future [29].
Dynamics 2022, 2, FOR PEER REVIEW 15 conformation tensor data is only mildly worsened. Next, the value of beff=5.78 can be obtained from Equation (48) of Stephanou et al. [29], which differs from the value 3 〈 〉 =16. 15. The remaining two parameters, and , can be obtained by simply fitting the NEMD data, since they do not affect the values; we obtain = 0.4 and = 0.01. Figure 11 shows how well the new model can fit the simulation data for the , , , and elements of the dimensionless conformation tensor for the C48 PE system in steady shear. We observe that the predictions of the revised model are in remarkable agreement with the NEMD extracted simulation results over the entire wide range of shear rates considered, especially for the diagonal elements. To quantify how well the model predicts the simulation data, we calculated the sum of the squares of the residuals (i.e., the residual between the simulation value and the one obtained by the model); this turns out to be ≈0.2474. The corresponding comparison for the material functions , Ψ , and −Ψ is presented in Figure 12. Contrary to the very good agreement between the refined model predictions and the NEMD simulation data for the dimensionless conformation tensor, the comparison against the viscometric functions is less satisfactory. Deviations from the NEMD data are mainly observed at large shear rates in the case of the shear viscosity (panel (a)) and the second normal stress coefficient (panel (c)). As also mentioned by Stephanou et al. [29], this disaccord should be related to the postulated relation between the stress and conformation tensors, Equation (2), which stems from the assumption of purely entropic elasticity [48]. As such, a more accurate expression for the free energy needs to be invoked in the future [29].

Conclusions
Today, it is well established, both experimentally [11,12] and computationally [13][14][15][16][17][18][19][20][21][22], that polymer chains subjected to flow fields that possess a rotational contribution exhibit, in addition to deformation, a tumbling/rotational behavior, which has been shown to be unambiguously responsible for the appearance of a transient stress undershoot (following the overshoot) at high shear rates [23,24,38]. This rotational behavior has been related to the slippage of polymer chains, relative to their surrounding, which in constitutive models is considered, among other methodologies, via the use of a non-affine or slip parameter, [25][26][27]29,38]. Although this parameter has been exclusively considered a constant, evidence suggests it should be a function of the chain's aspect ratio [27]. Probably, with the sole exception of the works of Rallison and Hinch [32,33] and Beris et al. [34], no other constitutive model has considered a shear-rate-(and time-) dependent slip parameter.
In our present work, we modified a constitutive model [29] that has been quite successful in predicting the data (both on the level of the conformation tensor but also the viscometric functions, albeit not accurately enough) obtained from detailed atomistic NEMD simulations of unentangled polymer systems over a wide molecular weight span to accommodate a variable slip parameter. The central idea is that the increase of the slip parameter from its equilibrium (null) value should be both shear-rate-and time-dependent due to the increasing rotational contribution of the imposed shear flow as the shear rate increases. The revised model still accounts for the most significant effects realized in physical systems, such as anisotropic drag, finite extensibility, non-affine motion, variable chain relaxation, and a bounded non-equilibrium free energy, all together as introduced

Conclusions
Today, it is well established, both experimentally [11,12] and computationally [13][14][15][16][17][18][19][20][21][22], that polymer chains subjected to flow fields that possess a rotational contribution exhibit, in addition to deformation, a tumbling/rotational behavior, which has been shown to be unambiguously responsible for the appearance of a transient stress undershoot (following the overshoot) at high shear rates [23,24,38]. This rotational behavior has been related to the slippage of polymer chains, relative to their surrounding, which in constitutive models is considered, among other methodologies, via the use of a non-affine or slip parameter, ξ [25][26][27]29,38]. Although this parameter has been exclusively considered a constant, evidence suggests it should be a function of the chain's aspect ratio [27]. Probably, with the sole exception of the works of Rallison and Hinch [32,33] and Beris et al. [34], no other constitutive model has considered a shear-rate-(and time-) dependent slip parameter.
In our present work, we modified a constitutive model [29] that has been quite successful in predicting the data (both on the level of the conformation tensor but also the viscometric functions, albeit not accurately enough) obtained from detailed atomistic NEMD simulations of unentangled polymer systems over a wide molecular weight span to accommodate a variable slip parameter. The central idea is that the increase of the slip parameter from its equilibrium (null) value should be both shear-rate-and time-dependent due to the increasing rotational contribution of the imposed shear flow as the shear rate increases. The revised model still accounts for the most significant effects realized in physical systems, such as anisotropic drag, finite extensibility, non-affine motion, variable chain relaxation, and a bounded non-equilibrium free energy, all together as introduced in its predecessor [29]. We compared the predictions of the revised model against newly executed atomistic NEMD simulations of a short unentangled PE melt with a molecular length equal to C 48 . Although the predictions at large shear rates were not significantly modified, the revision amended the problems associated with having ξ-dependent zero-shear-rate viscometric functions [29,38], cf. Equations (4), and linear-viscoelastic properties [38], cf. Equations (5). It should be emphasized that although the revised model was not derived through the use of a non-equilibrium thermodynamics formalism [28,30], its thermodynamic admissibility still holds, since 0 ≤ ξ ≤ 1 (provided 0 ≤ ξ 0 ≤ 1). Additionally, it is a straightforward exercise to extend the model to entangled systems by following our recent work [38]. We expect that the future use of the refined model will allow for more reliable prediction of macroscopic viscoelastic behavior and, therefore, for the development of more reliable computational tools, aiming to tailor-design large-molecular-weight polymeric systems.