Probing the nonequilibrium dynamics of stress, orientation, and entanglements in polymer melts with orthogonal interrupted shear simulations

The concept of entanglements is fundamental to the theory of polymer dynamics [1,2]. In equilibrium, tube theories relate the time scales that govern the dynamics of individual polymer chains to the rheological properties of the melt by scaling relations that involve the number of entanglements per chain. However, it is not clear how this picture can be fully extended to the nonequilibrium regime. Changes in the entanglement structure are one of several relaxation mechanisms that play a role in the complex nonlinear viscoelastic behavior of entangled linear polymer melts under shear [3,4]. Dynamic equations for disentanglement have been proposed [5], but one interesting effect that has been overlooked is the coupling of entanglements to orientation and how this affects the relaxation mechanisms. Entanglements are usually built into constitutive models through the number of entanglements per chain, which is a scalar quantity. However, in melts that have been highly oriented by flow, the interaction of any chain with its neighbors becomes anisotropic. Neglecting this property can lead to problems when using constitutive models to calculate the change in components of stress under flow, when starting from an out of equilibrium state. In this paper, we simulate an orthogonal flow protocol to study the nonlinear rheology of startup shear, starting from states that are partially relaxed after being driven to steady state of shear, and compare the simulation results to predictions by widely used constitutive equations. We also analyze the evolution of the number of entanglements per chain and the orientation at multiple scales along the chain during the transient period of startup.

Several models for the nonlinear viscoelasticity of polymer melts have been proposed. The original tube model of Doi and Edwards [2] relates the polymer stress to the average orientation of the tube and integrates over the whole deformation history of the melt. More recent models incorporated other relaxation mechanisms, such as contour-length fluctuations, chain stretch, and convective constraint release (CCR) [3–7]. Some of these are the Doi-Edwards-Marrucci-Grizzuti (DEMG) model [8], which includes the effects of chain stretch on the simpler Doi–Edwards theory; the Mead-Larson-Doi (MLD) model [9], which incorporates CCR; and the Graham-Likhtman-Milner-McLeish (GLaMM) model [10], which starts from a stochastic differential equation describing the microscopic dynamics of a polymer chain that can be used to derive an equation for the stress. Another approach not based on tube theory is the slip-link model [11]. One of the more successful simple models proposed is the Rolie-Poly model [12], which includes reptation, chain stretch, and CCR as relaxation mechanisms, and is a single-mode version of the GLaMM model [10]. It has been shown to describe, at least semiquantitatively, several observations of nonlinear phenomena under many flow conditions, in both experiment and simulation [12,13].

One of the hallmarks of nonlinear rheology of polymer melts is the existence of stress overshoots during startup shear under strong flows [14,15], particularly with Rouse Weissenberg numbers $WiR>1$⁠. Overshoots have been observed extensively in both experiments and simulations [13,16–22], and connecting the overshoot in the stress components to changes in molecular properties has been a subject of debate in the literature [23,24]. The strain $γη∗$ at the peak of the viscosity $η$ varies from $γη∗≈1$ up to $γη∗≈10$⁠. The first normal stress difference $N1$ also goes through an overshoot, and its peak happens at a strain $γN1∗≈2γη∗$⁠. Simulations find the same qualitative features and also show that the overshoot in the stresses is mirrored by the components of the chain orientation tensor for low and moderate shear rates [13,25,26]. This is usually interpreted as a validation of the stress optical rule, which states that information about the stress is encoded in the orientation and can be correlated to birefringence. Simulations also show an overshoot in the stretch and an undershoot in the number of entanglements per chain at very high rates.

One interesting protocol for investigating nonlinear rheology of polymer melts is interrupted startup shear [17,27,28]. Experiments have shown that when one interrupts shear after reaching steady state and resumes it after a short waiting time, the observed transient overshoot is smaller during the second stage of shear. This has been attributed to changes in the entanglement structure of the melt during the first shear [29], but reductions in the overshoot during the second shear have also been observed for unentangled melts [30]. Associated changes in the first normal stress overshoot have also been observed. Interestingly, while the original overshoot for unentangled systems is recovered for waiting times similar to the longest equilibrium relaxation time (the Rouse time), for entangled systems it is only recovered on time scales longer than several disengagement times.

While this protocol is useful for studying the effects of relaxation on the size of the overshoot, it is challenging to perform in both experiments and simulations due to the unusual changes in the flow over time. However, it allows for a new type of experiment that so far has not been realized: restarting shear with a different flow direction after relaxation. Such a setup can probe the changes in relaxation due to starting shear from a configuration far from equilibrium. For an aligned chain, the local molecular environment is different in the directions along the backbone and in the plane orthogonal to the backbone, and this anisotropy in the entanglement structure might lead to a different stress response. A related concept is orthogonal superposition rheometry, where a fluid in steady state shear is subject to oscillatory shear in a direction orthogonal to the steady-state flow [31–39]. A similar set of experiments was conducted by Hurlimann and Kraft [40–42] in a series of papers where they measured the rheological response of low-density polyethylene melts under a sudden change in shear direction using a multidirectional rheometer. Their experiments are essentially orthogonal interrupted shear experiments with no waiting time between the end of the first stage of shear and start of the second stage. Their focus was on the changes observed in the overshoot in the viscosity when the first stage of shear was interrupted at different prestrains. For prestrains smaller than the strain needed to reach the peak of the overshoot, little change was observed in the size of the overshoot during the orthogonal stage of shear. However, for large prestrains (larger than the strain needed to reach steady state), a significant reduction in the overshoot during the orthogonal shear stage was observed. Jeyaseelan and Giacomin [43] showed that transient network theory predictions were in reasonable agreement with these experiments, but the absence of strong overshoots in the stress for the analyzed rates makes agreement during the transient regime easier.

Constitutive equations based on tube theory can capture the qualitative features of interrupted shear experiments. The Rolie-Poly model can predict not only the monotonic recovery of the overshoot as a function of waiting time when the first shear is stopped after steady state, but also the nonmonotonic recovery seen if one stops the first stage of shear during the overshoot [44,45]. The Doi–Edwards model, which only includes relaxation due to orientation of the chains, is also able to model the recovery of the overshoot, suggesting that while disentanglement might contribute quantitatively to the recovery of the overshoot it is not the only mechanism responsible for it [46].

In this paper, we describe nonequilibrium molecular dynamics simulations of interrupted startup shear flow. We perform simulations in which the second shear is in the same direction as the original or in a plane orthogonal to the plane of shear during the first startup. The stress overshoot is significantly larger for a given waiting time for an orthogonal second shear. We also show measurements of the normal stresses as well as molecular scale quantities such as the end-to-end distance and orientation during both stages of startup shear. All of the relevant measures show the corresponding overshoots and undershoots before reaching steady-state, though the time at which they occur is not the same for all quantities. We contrast these measurements to predictions by constitutive equations based on the tube theory and show that available models are not capable of capturing this behavior even qualitatively.

Molecular dynamics simulations were carried out with the package LAMMPS [47] and used the Kremer–Grest bead-spring model for linear homopolymers [48]. Monomers were modeled as beads with mass $m$ interacting via a purely repulsive truncated Lennard-Jones (LJ) potential,

All simulation results will be in terms of LJ units (distance $a$⁠, energy $u0$⁠, mass $m$⁠, and time $τ=am/u0$⁠). The simulated melts have a total of 184 000 beads, and each melt comprises polymers with length $N=40$ or 500 beads. The polymers are linear chains of beads connected by covalent bonds represented by the FENE potential,

where $K=30u0/a2$ and $R0=1.5a$⁠. All simulations were performed at a density of $ρ=0.85$⁠.

Table I lists parameters for the systems studied. These parameters have been measured extensively in the literature [13,49,50]. Recently, Svaneborg and Everaers have published values of these parameters for many different stiffnesses [51], and values in the table are equal to the values quoted in their paper when these are readily available. All values quoted were consistent with direct measurements from our simulations. For the time scales, we use the values in the literature for the entanglement Rouse time $τe$ and obtain the chain Rouse time $τR$ and disentanglement time $τd$ from scaling relations. The scaling relation for the Rouse time is the standard $τR=τeZ2$⁠, while for $τd$ we use the scaling relationship by Likhtman and McLeish [52], which accounts for contour length fluctuations,

where $Z=N/Ne$ and $Ne$ is the entanglement length.

For all simulations, the equations of motion were integrated with a time step $δt=0.01τ$⁠, and temperature was kept at $T=1.0u0$⁠. Equilibrated melts were generated using the standard double-bridging algorithm [50]. Temperature during equilibration was controlled by a Langevin thermostat with a time constant of 10$τ$⁠.

The simulation protocol involved two stages of simple shear with a stage of relaxation in between. Figure 1 shows the shear plane and rate used during each stage of the simulations. In the first stage, the simulation box is sheared in the $xy$ plane at a desired rate. After the simulation has reached steady state, we stop the flow in a configuration where the simulation box is orthorhombic. The melt is allowed to relax for a certain time $tw$ and then sheared again. During the second stage of shear, we deform the box in either the $xy$ or $yz$ plane (where the first axis refers to the flow direction and the second refers to the gradient direction), leading to different transient stress responses. The rates used for shearing each system were chosen to produce Rouse Weissenberg numbers $WiR=γ˙τR$ ranging from 1 to 250. In this flow regime, nonlinear viscoelastic effects become relevant and the transient overshoot in the viscosity is very prominent [13,20,25,53].

During the relaxation stage, temperature was controlled by a Langevin thermostat with a time constant of 10$τ$ [54]. The temperature is obtained from the velocities of atoms after subtracting an instantaneous linear velocity profile that is calculated by binning the system in the $y$ direction and calculating the average $x$ velocity in each bin. This is done so that the remaining velocity profile after the first stage of shear does not affect the temperature at the beginning of the relaxation stage. Other methods for calculating temperature such as accounting only for the velocity component in the vorticity direction were also used and resulted in no measurable difference in the measured temperature and stresses.

Both stages of shear were simulated by integrating the SLLOD equations of motion and deforming the simulation box with the desired shear rate, and temperature was controlled with a Nose–Hoover thermostat with a time constant of 10$τ$ [55]. In our simulations, the flow homogeneity was ensured by directly probing the linear velocity profile across the system, and peculiar velocities were defined by subtracting the profile at a particle’s position from its velocity. During these stages, the thermostat was biased by subtracting a linear velocity profile consistent with the simulation box deformation rate before thermostatting and then reinstating the linear profile after the thermostatting was performed on the peculiar velocities.

Since the SLLOD equations of motion simulate a system under homogeneous flow, inhomogeneities such as shear banding are not seen in our simulations [55]. Other methods can be used to probe inhomogeneous flows. Recent Dissipative Particle Dynamics simulations with Lees–Edwards boundary conditions have shown that shear banding can be observed during startup shear of entangled polymer melts [56,57]. In that case, the melt separates into fast and slow bands where stress overshoots are still seen, but the values of the normal stress maxima differ from one band to the other.

The different types of thermostatting can lead to artifacts when calculating the temperature immediately after the flow is turned off or on. This is particularly important for the systems with small chains, where the relevant time scales of the dynamics are shorter and short time stress data can be affected by temperature fluctuations. In order to avoid this problem, the shearing stages were initialized with a velocity profile consistent with the desired shear rate, and peculiar velocities were adjusted to match the desired temperature. Other methods such as thermostatting only the direction perpendicular to the flow plane and quickly ramping the shear rate from zero to the desired value at the start of shear resulted in no significant change in the measured values of the stress in all cases tested.

Rheological data are obtained during the transient regime of startup shear and interrupted shear. Due to the noisy nature of stress measurements at individual time steps, we report stress data for single simulations that have been filtered using a Savitzky–Golay filter [58]. This filtering method works by fitting polynomials to sequential subsets of data points, effectively averaging data over a moving window without distorting the underlying signal.

Individual chain statistics such as length and orientation along different axes at different scales can be obtained by taking suitable averages over the positions of the monomers. Primitive path information can be obtained by one of the several methods published in the literature [59–61], usually by an algorithm that shrinks the backbone of the chains as short as possible without allowing chains to pass through each other.

Even in equilibrium, measurements of the number of entanglements by different estimators can vary; equivalently, one must determine the entanglement length $Ne$⁠, i.e., the number of bonds per entanglement. Hoy et al. [59] showed how different ways of calculating the entanglement length lead to different values. One possible measure, based on random walk statistics, is the rheological entanglement length $Nerheol≡Nb⟨R2⟩/Lpp2$⁠, where $Nb$ is the number of bonds, $⟨R2⟩$ is the average squared end-to-end distance, and $Lpp$ is the primitive path length of a chain. One can also calculate a topological entanglement length, given by $Netopol=Nb/Zk$⁠, where $Zk$ is the number of topological constraints (TCs) (points of contact between chains as obtained by an algorithm that reduces chains to their primitive paths) on a chain. It is usually found that the rheological entanglement length is larger than the topological entanglement length by a factor of around 2. To avoid confusion, we will use $Ne$ to refer to the rheological entanglement length in equilibrium, and this is the definition used for values in Table I and in scaling relations. Since this method does not work away from equilibrium, we use the Z1 code [62] to obtain primitive path statistics as well as information about contact between the primitive paths of different chains. The Z1 code works by shrinking the polymer chains to their primitive paths through a series of geometrical moves and identifying points where a move would force two chains to come in contact or cross. We call these contact points topological constraints (TCs) and analyze the evolution of the average number of TCs per chain $⟨Zk⟩$ under flow. Other methods such as Primitive Path Analysis [61] can generate similar information, and a direct comparison between TC statistics generated by different algorithms in aligned melts is the subject of undergoing research.

We use the Rolie-Poly model as well as a differential version of the Doi–Edwards model [63] to calculate the stress tensor during both stages of startup shear and during relaxation. The differential version of the Doi–Edwards model is

where $σ$ is the polymer stress and $κ$ is the velocity gradient tensor given by $κij=∂jvi$⁠. The operator $d/dt$ is the material derivative. This differential approximation to the original integral Doi–Edwards model is more suitable to numerical analysis. It has been shown to exhibit much of the same behavior as the integral version, except that it predicts a zero second normal stress difference. Like its integral counterpart, this model does not show an overshoot on the first normal stress difference, predicting a monotonically increasing $N1$ that saturates at the steady-state value [15].

The Rolie-Poly [12] equation is a single mode version of the GLaMM [10] model that includes contributions to relaxation from reptation, retraction, and CCR. The Rolie-Poly constitutive equation is

where the parameter $β$ sets the relative strength of the CCR relaxation. The disentanglement and Rouse times are related by Eq. (3). To replicate the flow conditions in our repeated startup shear MD simulations, we use different velocity gradient tensors to represent the different stages of the MD simulation when calculating the polymer stress using the constitutive equations. The chosen $κ$ controls which components of the stress tensor couple to flow in the advective terms of Eqs. (4) and (5). The most general velocity gradient tensor needed for our Rolie-Poly calculations is given by

where we remind the reader that the first index defines the velocity direction $v$ and the second index defines the gradient direction $∇$⁠, i.e., $κab≡κv∇$⁠. The two shear protocols are shown in Fig. 1. The first stage of shear always has nonzero $κxy$ (⁠$κyz=0$⁠), while the shear after relaxation can have either $κxy$ or $κyz$ as the nonzero component of the velocity gradient tensor. The choices $κxy$ or $κyz$ couple different components of the stress tensor to the flow during the second stage of shear while the uncoupled components continue to relax.

For both the parallel and orthogonal interrupted shear protocols, the viscosity $η$⁠, normal stress differences $N1$ and $N2$⁠, and chain statistics were measured for different values of shear rate $γ˙$⁠, waiting time $tw$ between the two shear stages and parameters in Table I. We contrast the transient stresses during the second shear to the transient stresses during the first shear and show that the main features apparent in the transient period vary monotonically with waiting time $tw$⁠. The simulation data show that the magnitude of the transient peak in all relevant properties is larger during the second shear when the flow is resumed in a direction orthogonal to the previous shear. We then show that both the Doi–Edwards and Rolie-Poly models predict the same magnitude of the overshoot in viscosity regardless of the relative orientation of the flow during the second stage, contradicting the simulations.

The measured values of $η$ and $N1$ during the first stage of startup shear are summarized in Fig. 2. Also plotted are components of the bond orientation tensor $Sij=⟨uiuj⟩$⁠, where $u$ is the bond vector connecting neighboring beads along the chain and $i$ and $j$ are Cartesian indices. Measurements of the orientation are performed on a single snapshot of the simulation, leading to noisier data. We do not show data for $Syy−Szz$ or $N2$⁠, but their values are always negative and very small in magnitude, with an undershoot consistent with the observed stress behavior. Under startup shear, according to the stress-optical rule, the stress components are proportional to the components of the bond orientation tensor,

where $α$ is the stress-optical coefficient. For the fully flexible FENE model, Cao and Likhtman [13] showed that $α=0.32a3/u0$⁠. For the entangled system, orientation and stress components track each other closely with minor deviations at the highest Weissenberg numbers. This can be attributed to the higher degrees of chain stretch and disentanglement induced at these stronger flows. For more moderate shear rates, the stress overshoot seems to be completely accounted for by the change in orientation, as predicted by the stress-optical law. For unentangled systems, the shear stress is significantly larger than orientation for all Weissenberg numbers larger than one. One possible explanation for this is that for unentangled systems the polymer contribution to the stress is small and, therefore, on the same order of magnitude as the contribution to the stress from unentangled interchain interactions.

Previous simulations of startup shear flow of entangled polymer melts have also observed weak undershoots on $⟨Zk⟩$⁠, the number of TCs per chain as measured by the Z1 code, during startup [21].

The strain at the peak of overshoot or undershoot is different for the shear and normal stresses. In experiments and simulations of both entangled and unentangled melts, the strain at the peak of the overshoot in $N1$ is usually twice as large as the peak strain in $η$ [20,26,30,64]. This is consistent with our findings, and we also observe even larger strains at the peak for the chain end-to-end length.

Figure 3 shows the viscosity, normal stress differences, and end-to-end length measured during the first and second stages of shear for a melt of chains with N = 500 monomers at a Rouse Weissenberg number $WiR=100$ and second shear in the same plane as the initial shear. Components of the bond orientation tensor during the overshoot are included for comparison. The orientation is correlated to the stress in a manner similar to the one observed during the startup from equilibrium, and the difference between stress and orientation is smaller for shorter waiting times. The size of the stress overshoot during startup shear increases monotonically with the waiting time $tw$⁠. The same behavior is observed for the first normal stress difference and end-to-end distance, as well as for the undershoot in the second normal stress difference.

Similar features can be observed when the simulations are repeated with a melt of short chains, as can be seen in Fig. 4. The dependence of the size of the overshoot on the waiting time between stages of shear remains even in the absence of entanglements between chains. In the unentangled case, full recovery of the overshoot occurs in the Rouse time, while for the entangled melt full recovery happens for waiting times longer than the disentanglement time.

Both sets of simulations show that the strain at the overshoot peak is about twice as large for $N1$ as it is for $η$ and $N2$⁠. The strain at the peak of viscosity seems to be independent of the waiting time, while the strain at the peak of the first normal stress difference decreases for the entangled system, moving from a position very close to the onset of steady state for short waiting times to a smaller value at the largest waiting times. The strain at the peak of the overshoot in the end-to-end distance of the chains also appears to decrease with increasing waiting time.

Both the Doi–Edwards and Rolie-Poly models qualitatively capture the recovery of the overshoot in the viscosity [44,46]. Neither the Rolie-Poly nor the differential Doi–Edwards model used here can predict a nonzero value for the second normal stress difference. The integral Doi–Edwards model predicts a value of $N2=−2N1/7$⁠, but neither form can predict an overshoot in $N1$ [63]. Figure 5 shows the values for the peak of the overshoot in the viscosity, $N1$⁠, $N2$⁠, and in the end-to-end distance obtained in the simulations as a function of waiting time. The maxima/minima for all properties for the unentangled system reach the terminal value (identical to the obtained when starting from equilibrium) after a waiting time around the Rouse time, which is the longest relaxation time for unentangled chains, while the entangled systems only reach the terminal value for the longest waiting times studied, which are over four times larger than the disentanglement time. This is consistent with trends observed in experiments on polystyrene melts [30]. A conservative estimate for the error bars on points in Fig. 5 can be obtained by looking at the variation in the ratio $ηmax/ηss$ and other quantities during startup from several equilibrium configurations prepared independently. By simulating startup shear of five different equilibrium melt configurations, we estimate error bars to be of order $±$0.25 in the viscosity and first normal stress ratios. Using the same method to estimate an error bar for the ratio of the second normal stress leads to an error bar of $±$0.75. This shows that all melts at the longest waiting time are statistically equivalent to samples sheared from equilibrium.

Both the entangled and unentangled systems exhibit a monotonic recovery of the overshoot, with corresponding changes to the bond orientation tensor.

Figure 6 compares predictions of the Doi–Edwards and Rolie-Poly models to the viscosity in the molecular dynamics simulations of the fully flexible entangled system. The best Rolie-Poly model fit parameters were $Z=9$ and $β=1$⁠, while the shown Doi–Edwards fit uses $Z=3$⁠. The Rolie-Poly value of $Z$ is consistent with $Ne=60$⁠, which is similar to the value used to fit the Likhtman–McLeish [52] and GLaMM [10] models to stress data by Cao and Likhtman [13]. Although this value differs from the rheological entanglement length $Ne=86$ reported in the literature, previous work has also shown that constitutive models using fit parameters derived from MD simulations do not necessarily correctly match the measured values of either the steady state or the overshoot [21]. The value of $Z$ used in the Doi–Edwards fit is significantly smaller than the simulation and Rolie-Poly values. One possible explanation is the well-known excessive shear thinning predicted by the Doi–Edwards model. The ratio $ηmax/ηss$ grows very strongly with increasing $WiR$⁠, so a small value of $Z$ is needed to fit the values measured in the simulations for our range of Weissenberg numbers. The Doi–Edwards model more accurately describes the data for shorter waiting times, while predicting smaller overshoots for longer waiting times. The best agreement with the Rolie-Poly model is for the longest waiting times. This might be due in part to our choice of using a single mode model, and better agreement for shorter waiting times can probably be obtained by adding a spectrum of shorter relaxation modes.

We also analyzed how the peak strain changes with waiting time. We see a systematic trend of a decrease by a factor of 2 for the strain at the peak of $N1$ with increasing waiting time, while there is an increase by a factor of 2 for the viscosity peak strain. The latter is, however, hard to estimate since the shear stress curves tend to be flat near the peak, particularly for short waiting times, making it hard to measure the position of the peak precisely. Comparison with our chosen constitutive models is difficult, since the Doi–Edwards model predicts a strain at peak that does not vary with Weissenberg number and the Rolie-Poly model [65] (as well as the more complicated GLaMM model [22]) has been shown to strongly overpredict $γmax$ for flows in the regime $WiR>10$⁠. This also happens in our case, where our MD simulations show peak strains ranging from 2 to 11 for our different Weissenberg numbers, while Rolie-Poly predicts a peak strain of around 200 for our largest $WiR$⁠. Rolie-Poly does predict a decrease of $γmax$ for $N1$ with increasing $tw$⁠, but also predicts a decrease for the peak strain of viscosity, which goes against what we observe in simulations.

We repeated the MD simulations detailed above with a second shear flow in a plane orthogonal to the first shear. Most of the same features remain, with overshoots and undershoots similar to those seen in Sec. III B. However, the size of the stress overshoot for the same waiting time is significantly larger when the direction of the second shear is changed.

In this orthogonal protocol, $y$ is now the flow direction and $z$ is the gradient direction (previously the vorticity). Hence, both the flow and gradient directions are orthogonal. Later, we will refer to this as “fully orthogonal.” Here, the viscosity during the second stage of shear is given by $η=σyz/γ˙$⁠, and the normal stress differences are given by $N1=σyy−σzz$ and $N2=σzz−σxx$⁠.

Figure 7 shows data equivalent to that shown in Fig. 3 for the orthogonal case. While the qualitative features are the same, for the same waiting time the relative magnitude of the overshoots in $η,N1,N2,$ and $Ree$ is larger than in the parallel case. The most striking difference is for short waiting times: in the parallel case, the overshoot was almost completely gone, while for orthogonal shear there is still a significant overshoot even for waiting times smaller than the Rouse time, a relative increase of around 50%.

Another significant difference can be seen in the values of the second normal stress difference $N2$ for short waiting, which are significantly more negative than during the initial shear. This is mainly due to our change in the definition of $N2$ to accommodate the different flow direction: the inclusion of a contribution from $σxx$⁠, which is still not fully relaxed for the smallest waiting times, leads to a more negative initial value of $N2$⁠.

The increase in the overshoot of the end-to-end distance comes exclusively from the behavior of $Ry2$ along the new flow direction. The $x$ component, which was the largest one in the first stage of shear, decreases monotonically.

The relative magnitudes of the overshoots and undershoot are larger than for parallel interrupted shear for all Weissenberg numbers. Figure 8 shows the values of the overshoot at the peak, normalized by the steady-state values, for both entangled and unentangled systems at all Weissenberg numbers as a function of waiting time. Comparison with Fig. 5 shows that for the shortest waiting times the magnitude of the overshoot is significantly higher when the second stage of shear is orthogonal to the original shear plane.

Figure 9 compares the viscosity overshoot of the entangled system to predictions by both constitutive models, with the Rolie-Poly parameters $Z=9$ and $β=1$⁠. There is significant disagreement for the shortest waiting times at all Weissenberg numbers larger than one. In fact, the model predictions are essentially the same as those for the parallel interrupted shear. An analysis of the form of Eqs. (4) and (5) explains this behavior: the coefficients of the relaxation mechanisms in both models involve only scalar quantities such as $trσ$⁠, which have no directionality and, therefore, cannot distinguish stresses related by rotations. The steady-state properties predicted by the models are the same as for parallel interrupted shear.

Figures 6 and 9 together suggest that there is no physical mechanism contained in the analyzed constitutive models that can explain the increase in the relative size of the overshoot due to the change of shear direction. Including extra faster relaxation modes can increase agreement in the short waiting time regime for either parallel or orthogonal interrupted shear, but the same spectrum of relaxation times cannot fit the different values of the overshoot measured in these regimes in the MD simulations. In Sec. IV, we will highlight measurements of chain orientation and entanglements that we believe can help elucidate possible physical mechanisms driving the observed differences in overshoot size. In particular, we propose that the anisotropic nature of the chains in the aligned state leads to remarkably different dynamics of entanglement in the transient regime when the flow is restarted in a direction orthogonal to the previous alignment direction, generating a contribution to the stress that arises from the larger number of constraints that need to be overcome in order for the melt to flow into steady state. We argue that constitutive models can be improved upon by utilizing a tensorial definition of entanglements instead of a scalar value of Z.

An important metric in polymer rheology is the number of entanglements per chain. Recent computational studies have focused on how the average number of topological entanglements per chain $⟨Zk⟩$ changes under flow in a United Atom model for polyethylene (UA-PE) [18–21,53], with this number usually measured by the Z1 code. They observe a decrease in the average $⟨Zk⟩$ as a function of Weissenberg number in steady-state shear, and a small undershoot in $⟨Zk⟩$ during the transient under startup shear roughly corresponding to the overshoot in the first normal stress difference. The UA-PE model is much more detailed than the simpler coarse-grained FENE model used in our simulations, and differences in quantities such as the friction coefficient between chains might change the degree of disentanglement at comparable Weissenberg numbers.

While it is still unclear how to couple the dynamics of entanglements as observed in simulations to constitutive equations, Ianniruberto and Marrucci proposed a model for the entanglement dynamics based on CCR [5], which has been used along with constitutive equations to model disentanglement during flow [66,67]. A direct comparison of predictions of this equation to the evolution of the number of entanglements in simulation as measured by the Z1 code is hard since a direct relationship between the concept of entanglements in tube models and those measured by simulation techniques is still not known [68]. In order to determine how $⟨Zk⟩$ is affected by our flow protocol, we analyzed the number of entanglements per chain during both parallel and orthogonal interrupted startup shear as a function of waiting time.

Entanglement data during interrupted shear for the entangled system are shown in Fig. 10. During parallel interrupted shear, we observe an undershoot in the number of entanglements per chain at around the same time as the overshoot in the end-to-end distance of the chains. The recovery of the undershoot as a function of waiting time is monotonic, similar to the transient behavior of other measured quantities. We also observe an overshoot during the transient before the undershoot happens, corresponding to a strain similar to the strain at the viscosity overshoot. One possible origin of the brief overshoot is that disentanglement is triggered by the nonaffine deformation of the chains in the nonlinear regime by mechanisms such as CCR, while deformations in the linear regime either preserve the number of entanglements or drive their formation. A more precise mechanism for forming entanglements can be formulated in terms of the distribution of entangled chains. A melt has a distribution of values of $Zk$⁠; less entangled chains could have a faster response to the imposed flow that is different from the more entangled chains, which may also depend on the direction of alignment of different chains. In our Z1 analysis, we observe that this brief rise in the average number of entanglements is mostly due to chains with a below average number of entanglements becoming more entangled with their neighbors, while strong disentanglement of the more entangled chains does not happen until a time scale of roughly $10−1τR$ (as seen in the sharp drop of $⟨ZK⟩$ in the purple curve in Fig. 10).

During the orthogonal interrupted shear protocol, the evolution of the average number of entanglements $⟨Zk⟩$ differs remarkably from that observed when shearing from equilibrium. Instead of going through an undershoot into steady state, $⟨Zk⟩$ undershoots and then rises again to a value significantly above the expected steady-state value, finally decreasing on a time scale larger than $τR$⁠. The initial overshoot before the first undershoot is also significantly larger than in the parallel case even for the shortest waiting times. We will demonstrate below that this interesting behavior can be understood as arising as a superposition of different mechanisms that affect the entanglement structure and annotate on the bottom panel of Fig. 10 the regimes where each mechanism dominates the dynamics of entanglements,

• the large orientation along an axis orthogonal to the new flow (because the chains have not relaxed after the first shear pulse) enables the formation of new entanglements under shear in a mechanism by which chains with nearly parallel backbones are pulled over each other; this leads to a large rapid increase in the entanglement number; this can be seen in the bottom panel of Fig. 10, where the first and second (from left to right) pictures of the chain “cloud” have a similar width, indicating a similar alignment along the original flow direction, but the second picture shows a cloud with a larger height, indicating that chains have started to align along the new flow direction;

• disentanglement is triggered by the nonlinear deformation of the chain along the new flow axis; in Fig. 10, this can be seen in the difference between the second and third pictures of the “cloud,” where the width of the distribution decreases as $⟨Zk⟩$ goes down;

• reorientation and stretch relaxation in the new vorticity axis lead to reentanglement by forcing chains to collide when aligning along the flow direction; this can be seen in the plot by noticing that the fourth and third pictures of the cloud have a similar aspect ratio despite a large change in $⟨Zk⟩$⁠, because most of the changes in the chain orientations are along the axis pointing into the page; and

• in the final stage, there is disentanglement to the steady-state value after alignment with the new flow axis. Notice how the width of the chain cloud decreases from the fourth to the final picture shown in the plot.

This result suggests that TCs between chains, as measured by the Z1 code, are inextricably linked to the chain conformation, and a successful model for the dynamics of entanglements should not only describe changes in the average number of entanglements per chain, but also how the created/lost entanglements couple to orientation and flow.

One possible way of modeling the evolution of the number of entanglements per chain (or entanglement density) under flow is to assume that the rate of change of the number of entanglements is proportional to the rate of change of the tube length, which can be expressed in terms of components of the bond orientation tensor $S$ and the chain conformation tensor $A=⟨RR⟩/2Rg2$ (where $Rg$ is the radius of gyration and $R$ is a chain’s end-to-end vector), possibly coupled to the velocity gradient tensor $κ$⁠. This approach was used by Ianniruberto and Marrucci [5] to obtain an equation for the change in $ν=Z/Zeq$ under flow which includes a term proportional to $κ:S$⁠. In order to identify relevant changes in the components of $S$ and $A$ that might be responsible for the curious behavior of the entanglement density under orthogonal interrupted shear, we plot the individual components of these tensors during parallel and orthogonal interrupted shear after our shortest waiting time of $0.27τR$ in Figs. 11 and 12. The diagonal components are shown in the insets of the plots and exhibit little change for parallel interrupted shear while showing a reversal of the relevant components for the orthogonal case (relaxing the $xx$ component while increasing the $yy$ component), similar to the behavior observed for the stresses.

All components of orientation shown in Fig. 11 show the same qualitative behavior as the stress curves, as predicted by the stress-optical law. No components under orthogonal shear show the “double peaks” observed for $⟨Zk⟩$⁠.

The conformation tensor $A$ combines the orientation and stretch at the chain scale into a single measure. The diagonal components of $A$ under flow show very little change in the steady-state values of the individual components in parallel interrupted shear but an inversion between $Axx$ and $Ayy$ in orthogonal interrupted shear, where $Axx$ relaxes in the same time scale that $Ayy$ rises to become the dominant component in $A$⁠. The trace of $A$ is the same as the average end-to-end distance in Fig. 3(d), just squared and normalized by the square of the radius of gyration.

An interesting effect of the change in flow direction in the second stage of shear can be seen in the $Axy$ component of the chain conformation tensor during orthogonal interrupted shear. During parallel interrupted shear, both off-diagonal components that do not couple to the flow in $κ:S$ are zero at all times, and one might expect the same would happen to $Axy$ after switching the flow direction. However, our simulations show that $Axy$ goes through a more complicated relaxation process, which exhibits a prolonged undershoot, at time scales comparable to the complex behavior appearing in the Z1 data, and featuring two peaks. The shear component $Ayz$ exhibits an overshoot similar to the stresses as shown in Fig. 7 and is not qualitatively different from what is observed for parallel interrupted shear when taking into account the appropriate flow, shear, and vorticity directions (apart from changes in the relative size of the overshoot for the same waiting time, as discussed in Sec. III C).

At first glance, the qualitative difference observed between $Sxy$ and $Axy$ might seem to come from the contribution of stretch, which is accounted for in $A$ but not in $S$⁠. However, another important difference between the two measurements is the scale over which these measurements are taken. The bond orientation tensor $S=⟨uu⟩$ measures orientation at the monomer scale, since $u$ is the bond length between neighboring beads in a chain, while $A$ measures the stretch and orientation of the whole chain. We define orientation and conformation tensors at different scales as $S(n)=⟨unun⟩$ and $A(n)=⟨RnRn⟩/2Rg,n2$⁠, where $Rn$ is the distance vector connecting two beads n monomers apart along the chain, $un=Rn/|Rn|$ is the normalized distance vectors, $Rg,n2=Rn2/6$ is the equivalent radius of gyration, and averages are taken over all suitable pairs of monomers. In Fig. 13, we show $Sxy(n)$ and normalized $Axy(n)$ measured at different relevant scales along the chain during orthogonal interrupted shear in order to understand the contribution of orientation at different scales to the double peak feature observed in $Axy$⁠. We see that while the orientation in the monomer scale does not undershoot the steady-state value, the orientation at larger scales does show an undershoot at a time scale close to $τR$⁠, and at scales larger than $Ne$ the prolonged undershoot also appears in $Sxy(n)$⁠. Interestingly, the double peak only seems to appear in the conformation tensor and only for scales comparable to the whole length of the chain.

While the relevance of these changes in the $xy$ components of $S$ and $A$ during orthogonal interrupted shear is not clear, a term including off-diagonal components of orientation and conformation that do not couple to flow through terms such as $κ:S$ in a kinetic equation for $dZ/dt$ might be able to explain the difference in the dynamics of $Z$ in our two different flow protocols and can be investigated more thoroughly in future work. Terms that contain a relevant contribution from $Sxy$ and $Axy$⁠, such as those involving tensor invariants like the determinant, might also need to account for the flow history of the melt at different scales along the chain.

We used nonequilibrium molecular dynamics simulations to probe the transient rheological properties of polymer melts under startup shear, both starting from the equilibrium configuration and from a state that has relaxed after steady state of shear for a waiting time $tw$⁠. The observed behavior during startup shear from equilibrium is in line with what has been previously reported in the literature and shows that the stress optical law holds up to large Weissenberg numbers (⁠$WiR>100$⁠).

Our simulations of interrupted startup shear show a monotonic recovery of the transient viscosity overshoot with increasing waiting time, as is seen in experiments. Recovery of associated overshoots and undershoots in the normal stresses, end-to-end distance of the chains, and the number of entanglements per chain are also observed. Constitutive models based on tube theory can qualitatively reproduce the stress measurements, as has been demonstrated extensively before. Semiquantitative agreement can be obtained by carefully tuning model parameters, though these do not necessarily reflect the parameters as measured in simulations (such as the number of entanglements per chain). Better agreement can probably be obtained if one uses multimode versions of the constitutive equations, which we have not used.

By performing simulations of orthogonal interrupted shear, where shearing is resumed in a different plane after the waiting time, we were able to show that these same models fail to accurately predict the recovery of the overshoot even qualitatively. Both the differential Doi–Edwards model and Rolie-Poly predict the same overshoot size regardless of the relative orientation of the shear flow during the second stage of shear. Figure 14 shows the ratio of the peak value of the viscosity overshoot during orthogonal interrupted shear to parallel interrupted shear in our simulations and in predictions by both constitutive models studied.

The maximum viscosity under orthogonal interrupted shear is up to 50% larger than the maximum under parallel interrupted shear for the systems studied, while the constitutive models predict at most differences on the order of 3% for very short waiting times. We note that the maximum viscosity values predicted by our Rolie-Poly model are smaller than that obtained by MD simulations for short waiting times (as can be seen in Figs. 6 and 9) and, as a result, a 3% increase is significantly smaller when compared to the 50% increase seen in MD simulations. It is easy to see that such behavior cannot be modeled by the constitutive models due to the lack of any anisotropic quantities in the incorporated relaxation mechanisms. When integrating the differential equations, the only differences between the second and first stages of shear are the relative magnitudes of the coefficients for the different relaxation mechanisms, which depend on time exclusively through scalar quantities. Small differences between the maximum viscosity in orthogonal and parallel interrupted shear can be seen for short waiting times due to the relatively large difference in initial values of the components of the stress tensor. Better models that include the effects of chain orientation or an anisotropic number of entanglements (representing the very anisotropic nature of neighboring chains in aligned states) in the relative magnitudes of the relaxation mechanisms can plausibly solve this issue. A coarse-grained model for star polymer melts that incorporates a mechanism similar to that described above was moderately successful at capturing some of the nonlinear rheological behavior observed in experiment [69].

We used the Z1 code [62], a standard method for analysis of primitive paths of entangled polymers, to determine changes in the number of entanglements per chain during both types of interrupted shear. In parallel interrupted shear, a simple undershoot is observed during the second stage of shear, mirroring the overshoot seen for quantities such as the viscosity, first normal stress difference, and end-to-end distance. Increasing the waiting time between stages of shear controls the relative size of the undershoot in $⟨Zk⟩$⁠, which changes monotonically with $tw$⁠. For orthogonal interrupted shear, the same behavior is observed for long waiting times, but short waiting times lead to peculiar behavior of $⟨Zk⟩$⁠, which develops two successive overshoots during the second stage of shear. The only components of orientation and conformation that show similar behavior are the relaxing off-diagonal components $Sxy$ and $Axy$⁠, which do not couple to flow (since $κxy=0$⁠) in the second stage of shear in standard models of polymer rheology.

Our work suggests that constitutive models that couple stress to polymer orientation and entanglements could be modified to accommodate the directional nature inherent in actual entanglement constraints (rather than a smooth scalar quantity), which is evidently discriminated by complex time-dependent flows such as this example of orthogonal interrupted shearing flows. Future work will focus on the effect of chain stiffness on the phenomena described in this paper. In particular, some of our preliminary results show that stress and $⟨Zk⟩$ data for melts of stiffer chains still exhibit qualitatively the same features observed here.

Mark O. Robbins sadly passed away during the writing of this manuscript. We thank Vicky Nguyen, Ben Dolata, Jon Seppala, Kalman Migler, Gretar Tryggvason, Martin Kröger, and Thomas O’Connor for advice and discussions. This work was funded by NSF (No. DMREF 1628974), the Ives Foundation, and Georgetown University.

The authors have no conflicts to disclose.

The orthogonal interrupted shear protocol discussed in the main text rotates both the velocity and gradient by $90°$⁠: the gradient direction in the first stage becomes the flow direction in the second stage, while the flow in the first stage becomes the neutral vorticity direction in the second stage. Inverting the flow and gradient direction in the second stage of shear such that the gradient direction remains the same during both stages of shear and the initial vorticity direction becomes the flow direction during the second stage can lead to a different stress response. This second flow is partially orthogonal or “orthogonal-parallel” since the gradient direction remains the same, while the first flow is “orthogonal-orthogonal.”

In order to test this, we repeat the simulations detailed in the main text using a velocity gradient tensor given by

During the first stage of shear, $κxy=γ˙$ and $κzy=0$⁠, while during the second stage $κxy=0$ and $κzy=γ˙$⁠. This geometry is similar to that used in experiments by Kraft [41,42].

Figure 15 shows a comparison between the viscosity and first normal stress difference at $WiR=100$ for the three strain protocols (parallel and the two different kinds of orthogonal) for short and long waiting times.

As expected, there is no difference in the stress response after a long waiting time. However, for a short waiting time, the peak value of the viscosity lies between the peaks of parallel shear and the orthogonal-orthogonal shear discussed in the main text, while the peak in the first normal stress difference shows no difference from the one observed for parallel interrupted shear.

The number of entanglements per chain also displays some similarities to the behavior described in the main text for short waiting times (Fig. 16). The number of entanglements rises as the melt is sheared and then decreases when approaching steady state. Although the initial rise and final decrease in entanglements are very similar to that observed in the fully orthogonal shear protocol, at short waiting times it lacks the double peak structure found between $0.1τR$ and $τR$ in that case (compare Figs. 10 and 16). The characteristic two peaks are recovered for intermediate waiting times, while for the longest waiting times we observed the expected behavior similar to a standard startup shear, with a decrease in $⟨Zk⟩$ followed by an undershoot before reaching steady state. The orientation and conformation tensor also lack a clear sign of independent peaks in the off-diagonal components.

The orientation evolution during orthogonal-parallel interrupted shear also shows striking differences at short waiting times, compared to that observed for the orthogonal-orthogonal shear (see Fig. 17). The $σxz$ component of the stress tensor, which is not the primary imposed stress during either the first or second stages of shear (where $xy$ and $zy$ are, respectively, imposed) has a pronounced overshoot in the orthogonal-parallel case, which is not seen in the other flow protocols. Since the $xz$ component couples the flow directions during the first (⁠$x$⁠) and second (⁠$z$⁠) stages of shear, we speculate that this overshoot is due to the realignment of chains when the flow direction is changed, leading to a noticeable bump in stress. This is similar to the small undershoot observed during startup that we associate with the tumbling of chains.

As described by the stress-optical law [Eq. (7)], this change in orientation leads to a change in the corresponding component of stress. This change in $σxz$ is shown in Fig. 17. As a general rule, all flow protocols satisfy the stress-optical law for the entire orientation and stress tensors, with systematic deviations at higher Weissenberg numbers (see Fig. 2) and for unentangled melts. Kraft [41,42] performed experiments with a similar setup and an effective waiting time $tw=0$⁠. However, they only monitored $σxy$ and $σzy$⁠, observing behavior similar to the one we see in the simulations, and did not show measurements of $σxz$⁠. Experimental measurement of a nonzero shear stress in the $xz$ plane during orthogonal interrupted shear with constant gradient direction could then be another way of verifying the results of our simulations for model linear polymer melts.

In the main text, we showed changes in the average number of entanglements per chain $⟨Zk⟩$⁠, as measured by the Z1 code, as well as in the components of orientation and conformation tensors during parallel and orthogonal interrupted shear simulations of our entangled fully flexible melt at $WiR=100$⁠. For completeness, in this appendix, we show the changes in the same quantity for different Weissenberg numbers.

Figure 18 shows data for $⟨Zk⟩$ under parallel and orthogonal interrupted shear of our fully flexible melt for Weissenberg numbers from 25 to 250. While data are available for $WiR=1$⁠, the absolute change in the value of $⟨Zk⟩$ is very small because of the weak imposed flow.