3.1. Rheological Properties
Figure 1a–c display the steady-state state viscosity
η as a function of the imposed strain rate
for the simulated PEO-1k, PEO-2k and PEO-5k ring and linear melts. In each Figure, we have also added a second
x axis at the upper part of the graphs indicating the corresponding Weissenberg number
(subscript C coming from Cyclic) defined as
, i.e., with respect to the ring orientational relaxation time
under equilibrium conditions. Qualitatively, for all three different
Mw ring PEO melts simulated, we observe the typical characteristic behavior of the shear viscosity with shear rate already known for linear polymer melts: (a) At low shear rates (
), the viscosity is practically constant defining what we know as the Newtonian plateau. (b) At higher shear rates, the viscosity starts decreasing, exhibiting what we know from the corresponding behavior of linear polymers as shear thinning. (c) At even higher shear rates, we enter the highly nonlinear regime where the viscosity drops rapidly with applied shear rate. It is interesting that for the PEO-1k system, ring and linear melts in the Newtonian regime exhibit similar viscosity values. On the other hand, for the PEO-2k and PEO-5k systems, the viscosity of the linear melt in the Newtonian regime is larger than that of the corresponding ring, which is in full agreement with experimental measurements [
1,
2,
6,
13,
33] and previous simulation studies [
8,
12,
14]. Overall, it appears that differences in the viscosity between ring and linear PEO melts in the Newtonian regime become more important as the
Mw increases.
Given that the calculation of the stress tensor in atomistic NEMD simulations at very low shear rates (approximately below
) suffers from large statistical fluctuations, to compute the corresponding zero shear rate viscosity
η0 we fit the simulation data to the Carreau Model [
34,
35]:
where
p is a characteristic exponent and
λ a fitting parameter with units of time. All results derived from the Carreau model for
η0, along with their standard deviation, are listed in
Table 2. Very rough estimates of
η0 were also obtained from the longest relaxation times for linears and rings
(and
, respectively) as computed directly from the equilibrium MD simulations using the corresponding Rouse equations, namely
and
; they are reported in the fourth and fifth column of
Table 2. In the sixth and seventh column of
Table 2, we report available experimental data for the zero-shear rate viscosity of PEO [
12,
19].
Using Equation (6) to compute
η0,L proved very efficient allowing us to reliably estimate the ratio
of the zero-shear rate viscosity of the linear melt to the corresponding zero shear rate viscosity of the ring. As already mentioned, for the shortest melt examined (PEO-1k), the NEMD prediction for the viscosity ratio between linear and ring melt is
. This observation agrees nicely with the experimental measurements of Nam et al. [
36] for short ring and linear PEO melts at a lower temperature (
T = 329 K) than in our NEMD simulations (
T = 363 K), where the ratio
was equal to 1 for melts with
Mw ≈ 1500 g/mol. We remind the reader that according to the Rouse model the ratio
is equal to 2 [
10]. The reason for the deviation of the NEMD prediction from the Rouse model should be due to the excess free volume phenomena present in the linear melts due to chain ends, which accelerate chain dynamics and, in turn, cause the viscosity of the linear melt to decrease. Of course, as the chain length increases, these excess free volume phenomena become less and less important, thus we expect the ratio of the two viscosities to come closer to the value of 2.0 predicted by the Rouse theory. Indeed, for the PEO-2k melt, the corresponding ratio is
, a result which is fully consistent with the Rouse model [
10]. On the other hand, by further increasing
Mw, entanglements start developing between chains in the linear melt which restrict their dynamics; thus, now, we expect the viscosity of the linear melt to increase faster than the viscosity of the ring. Indeed, according to our NEMD simulations, for PEO-5k (a marginally entangled melt),
. The fact that a change in the relaxation mechanism takes place as we cross over from PEO-2k to PEO-5k that cannot be accommodated by the Rouse model is also reflected in the unrealistically large value of the ratio
(=8.5 ± 1), see data in the fourth and fifth column in
Table 2 for PEO-5k, predicted by naïve application of the Rouse model equations on the basis of the computed chain orientational relaxation times for the ring and linear melt from the equilibrium MD data. Overall, we can say that, according to our NEMD simulations, the ratio
for very short PEO melts (well below
Me) starts from a value close to 1 and increases smoothly with
Mw, reflecting the faster increase of the viscosity of the corresponding linear melt as inter-chain entanglements start playing a role. Our NEMD data for
η0 are also in very favorable agreement with experimentally measured data from Refs. [
19,
33].
In the non-linear regime, the viscosities of ring and linear melts behave qualitatively very similar for all chain lengths, but the degree of shear thinning in the linear melt is sharper compared to that in the corresponding ring melt, especially as the chain length increases. By fitting the viscosity curves with a power law of the form
we can extract the value of the exponent
b quantifying the rate of shear thinning with applied shear rate. For the ring melts, the values obtained are: (a)
b = 0.54 for R-1k, (b)
b = 0.54 for R-2k, and (c)
b = 0.58 for R-5k. For the linear melts, the corresponding values are: (a)
b = 0.54 for L-1k, (b)
b = 0.65 for L-2k, and (c)
b = 0.76 for L-5k. We see that the shear thinning exponents for the linear melts increase more rapidly with chain size than the corresponding ones for the rings. This agrees with recent simulation [
14] and experimental [
13] studies, and is typically attributed (at least to a large extent) to the reduced deformability of ring chains [
13,
14] due to their closed-loop structure. For rings, it has been further argued that ring-specific dynamical mechanisms (including cooperative orientation of parallel segments or tumbling) can lead to molecular individualism, i.e., to a different behavior of different ring molecules in the melt [
13,
14]. We will come back to this issue later in
Section 3.2 where we will examine the effect of flow on the rate of chain orientation of individual chains in the two types of melts (ring and linear).
In
Figure 2 and
Figure 3 we show the NEMD results for the first
and second
normal stress coefficient. We see that Ψ
1 is positive and Ψ
2 negative (this explains why in
Figure 3, we plot −Ψ
2), and that the magnitude of Ψ
2 is much smaller than the magnitude of Ψ
1. For weak-to-moderate shear rates (
), ring melts exhibit smaller Ψ
1 and −Ψ
2 values compared to linear melts. For higher shear rates, both Ψ
1 and −Ψ
2 exhibit large power law regions, similar to those for
, decreasing by several orders of magnitude with shear rate. We also observe that the rate of decline of Ψ
1 and −Ψ
2 with
is greater in the linear melts than in the rings. To quantify this difference, we fit the simulation data for Ψ
1 and −Ψ
2 with
for
with a power law of the form
and
. For PEO-1k the two exponents are similar:
b1 = 1.38 ± 0.05 and
b2 = 1.21 ± 0.05 for the ring, and
b1 = 1.38 ± 0.04 and
b2 = 1.41 ± 0.05 for the linear melt. With increasing chain size, the differences in the power law exponents between ring and linear melts increase. Thus, for PEO-2k,
b1 = 1.48 ± 0.04 and
b2 = 1.44 ± 0.09 for the ring melt, which should be compared to
b1 = 1.61 ± 0.03 and
b2 = 1.63 ± 0.03 for the linear melt. For the marginally entangled PEO-5k,
b1 = 1.55 ± 0.06 and
b2 = 1.53 ± 0.03 for the ring melt, while
b1 = 1.76 ± 0.05 and
b2 = 1.75 ± 0.04 for the linear melt. Our simulation findings are in reasonable agreement with the results of [
14] for marginally entangled (
Z = 6) linear and ring PE melts. In particular, the Ψ
1 and Ψ
2 exponents reported in [
14] are
b1 = 1.45 ± 0.03 and
b2 = 1.49 ± 0.11 for ring PE, which increased to
b1 = 1.61 ± 0.03 and
b2 = 1.67 ± 0.12 for linear PE. The finite extensibility of rings, due to their loopy geometry and more compact arrangement of monomers around their center of mass, should be considered as the main factor responsible for the weaker dependence of Ψ
1 and Ψ
2 on shear rate compared to linear melts.
Figure 4 presents the ratio
for ring and linear melts, for the three PEO systems studied. For all of them, rings exhibit larger
values compared to linear melts, which reflects again the reduced deformability of rings due to their closed loop architecture. As the molecular weight increases, for both linear and ring melts the ratio
decreases. Overall, and for the range of shear rates addressed in the present study, the ratio
for the ring melts ranges between 0.23 and 0.36 for R-1k, between 0.28 and 0.33 for R-2k, and between 0.11 and 0.21 for R-5k. For the linear melts, on the other hand,
ranges between 0.1 and 0.19 for L-1k, between 0.15 and 0.17 for L-2k, and between 0.05 and 0.10 for L-5k.
3.2. Conformational Properties
Figure 5a,b provide representative atomistic snapshots from the NEMD simulations with the marginally entangled PEO L-5k and R-5k melts, respectively, at four different shear rates corresponding to
numbers equal to 1, 10, 100 and 1000. Admittedly, the applied flow has a strong effect on melt conformation; as shear rate increases, strong stretching and significant alignment of ring and linear PEO molecules along the flow (
x-) direction are observed, causing significant deviations of their average shape from the equilibrium, isotropic configuration. To quantify this deviation, we have computed the average radius-of-gyration tensor as a function of applied flow. For a given chain, the
ab-component of the second order radius-of-gyration tensor
G is defined as
where
N is the number of atoms per chain,
and
denote the position vectors of atom
i and of the center-of-mass (com) of the chain this atom belongs to, respectively, and the subscripts
a and
b indicate the space directions
a and
b, respectively.
In
Figure 6,
Figure 7,
Figure 8 and
Figure 9, we compare the
xx,
yy,
xy and
zz components of the radius-of-gyration tensor (
Gxx,
Gyy,
Gzz and
Gxy) between linear and ring melts. We focus our attention first on the
xx component. At low shear rates (
), the magnitude of
Gxx remains practically unaffected by the flow for both types of melts. At intermediate shear rates, chains in the melt deform and at the same time align in the direction of the flow, which causes a considerable increase in the value of
Gxx. At even higher shear rates, the rate of increase of
Gxx with shear rate declines and
Gxx approaches constant values which are, however, different between ring and linear melts. The fact that
Gxx assumes eventually constant values is due to two factors: (a) the finite extensibility of the simulated chains due to bond stretching and bond-bending interactions, and (b) chain rotation and tumbling due to the nature of the applied flow (shear) [
32]. An interesting point to notice in the curves of
Figure 6 is that for all three PEO melts studied, the ratio of the asymptotic
Gxx values between linear and ring melts is very similar and approximately equal to 2, which seems to suggest that, as far as their fully extended conformations are concerned, rings can be considered as linear chains of half the length.
On the other hand, the two other diagonal components Gyy and Gzz display exactly the opposite behavior. At low shear rates, their values remain unaffected by the imposed flow field; however, as the strength of the flow increases, both decrease considerably due to chain alignment in the direction of flow. It is also true that Gyy decreases more rapidly than Gzz, which should have been expected given that z is the neutral axis. As already mentioned, due to the nature of the applied flow, the average melt velocity varies along the y-axis, and this can cause chain rotation and tumbling (see below), which is another reason why Gyy decreases faster than Gzz at higher shear rates. As with many other (rheological and conformational) properties discussed so far, the values of Gyy and Gzz for rings are smaller than for linear melts, which is another manifestation of the more compact structure of cyclic molecules due to the more symmetric arrangement of atoms around their center-of-mass.
The shear rate dependence of the
xy component of the radius-of-gyration tensor is presented in
Figure 9. In all cases,
Gxy initially increases with applied shear rate, goes through a maximum at an intermediate value of shear rate, and then starts decreasing. It is known [
32] that the overall change of
Gxy with respect to applied shear rate can be understood by considering two competing effects: (a) more open molecular conformations due to flow stretching in the
x-direction and the spatial correlations between the
x and
y components of the chain end-to-end vector, leading to an increase in
Gxy, and (b) chain orientation along the flow direction, leading to a decrease in
Gxy. It is the competition between these two effects that gives rise to the maximum. Overall, and despite the obvious similarities in the overall qualitative behavior of
Gxy with shear rate, certain differences are observed. For example, for linear melts,
Gxy increases rapidly with shear rate even for small shear rates whereas for rings the increase is slower and is initiated at higher shear rates. This is another manifestation of stronger resistance to the applied flow, and can be explained again by the more compact structure of rings compared to linear melts which can attain more open and wider conformations. Similar results have been reported by Yoon et al. [
14] for PE.
To check how the intrinsic principal dimensions of the ring molecules change with applied flow, we have also calculated the eigenvalues (
G1 >
G2 >
G3) of the radius-of-gyration tensor for all three different PEO melts simulated. The results found for the PEO-5k are reported in
Figure 10. Also reported in
Figure 10 is the ratio of the largest eigenvalue
G1 with the smallest
G3. Admittedly, the changes of the three eigenvalues with shear rate follow (to a large degree) the corresponding changes of the
xx,
zz and
yy components of the radius-of-gyration tensor
G.
We have also studied the effect of the applied shear flow on chain alignment. To this, we calculated the alignment angle defined as
where
A denotes a second order tensor, e.g., the radius-of-gyration tensor, the stress tensor or the birefringence tensor. In the present work, we choose to work with the radius-of-gyration tensor
G.
Figure 11, then, shows how the alignment angle changes with applied shear rate between rings and linears. At low shear rates,
θ tends to the value of 45° for both types of melts. The limiting value of 45° has already been observed in other simulation studies [
14,
32] and is also explained theoretically. Specifically, at low shear rates, the chains start to align, leading to a nonzero value of <
Axy> but with negligible chain deformation (i.e., negligible normal stresses). At intermediate shear rates, the chains align in the flow direction and, also, stretch considerably; as a result, normal stresses develop in the melt, which in turn cause a steep rise in the value of
θ. At high shear rates, chains reach their maximum allowable stretch, thus
θ reaches a plateau value. Comparing the results between ring and linear melts at sufficiently high shear rates, we notice that rings are characterized by a lower degree of alignment. Once more, this demonstrates the stronger resistance of rings to the applied flow in comparison to linear chains [
14] as a result of their shorter spatial extent due to their closed structure.
3.3. Terminal Relaxation
In the literature, it has been suggested that, similar to their linear counterparts, ring polymers can stretch, collapse, tumble and re-stretch when subjected to shear flow [
37,
38]. Moreover, due to their loopy structure, they can undergo tank-treading (i.e., a rotational motion where chain segments rotate along the deformed contour of the polymer) [
37,
38]. Simulations of dilute ring polymer solutions have shown [
37,
38] that at weak flow rates it is not easy to distinguish tumbling from tank-treading. In contrast, at strong flow rates, tumbling and tank-treading can occur independently. We observed similar tumbling and tank-treading dynamics in our NEMD simulations with the PEO melt systems addressed here. This is clearly discernible in
Figure 12a,b, showing representative atomistic snapshots of two ring chains from the NEMD simulations with the R-5k melt under strong flow conditions (
WiC = 1000) undergoing tumbling and tank-treading motion, respectively. The interested reader can also visualize the two types of motion for the chosen pairs of ring molecules in the two videos that we prepared and uploaded as
Supplementary Material to this manuscript from the NEMD simulation with the R-5k melt at
WiC = 1000.
Although polymers are characterized by a broad spectrum of characteristic times that span a wide range of scales, the most important one is the time associated with the relaxation of length scales on the order of the entire chain, defining what we call orientational or terminal relaxation. For linear melts, terminal relaxation is quantified by calculating the orientational autocorrelation function (OACF) <u(t)·u(0)>L of the unit vector u directed along the chain end-to end vector (Ree), and its time decay in the course of equilibrium MD and NEMD runs. For ring melts, terminal relaxation is quantified by looking at the OACF <u(t)·u(0)>R of the unit vector u directed along the diameter vector (Rd) of the ring, averaged over all possible such vectors for a given ring molecule. The rate with which <u(t)·u(0)> approaches the zero value is a measure of how fast the chain forgets its initial configuration, i.e., of the rate of the overall orientational relaxation of the chain.
Figure 13 shows the spectrum of the individual OACF functions separately for each chain in the simulated PEO-1k ring and linear melts, together with the corresponding average values (thick black lines), under equilibrium conditions and under strong shear flow corresponding to
and
, respectively. The simulations have been carried out until the individual
u(
t)·
u(0) functions for all chains in the melt have dropped to zero. Surprisingly, our MD simulations reveal that, for both molecular architectures studied (ring and linear), the orientational dynamics are highly heterogeneous, since the individual
u(
t)·
u(0) curves deviate significantly from the average <
u(
t)·
u(0)> curve. In fact, this dynamic heterogeneity is more pronounced in the case of linear melts where, in addition, chains require much longer times to relax. For the PEO-1k system examined in
Figure 13, this behavior is counter-intuitive, since it is unentangled and, in addition, excess free volume around chain ends should accelerate chain relaxation in the case of the linear melt [
12]. However, we can explain this peculiar behavior if we recall that, due to their looped structure, ring chains assume conformations that are spatially extended to shorter distances than chains in the linear melt. That is, terminal relaxation for linear chains involves the decorrelation of the unit vector along the end-to-end chain vector
Ree which corresponds to twice the molecular length covered by the diameter
Rd vector used to define the OACF function in the case of rings. Under strong flow conditions (
and
), both types of melts relax faster but the main features of the decorrelation still remain. With increasing molecular weight (see
Figure 14), the effect of excess free volume in the linear melts becomes less and less important, and this is reflected in the enhanced dynamic heterogeneity of linear chains. Similar behavior was observed for the entangled R-5k melt (see
Figure 15).
By integrating the time autocorrelation functions <
u(
t)·
u(0)>
R and <
u(
t)·
u(0)>
L for all chains in the melt, one obtains a measure of the characteristic time constants for relaxation for ring and linear melts, respectively. The corresponding histograms of these times are displayed in
Figure 16 and
Figure 17. For all PEO systems studied, characteristic relaxation times for chains in the linear melt span a much wider range than in the ring melt. With increasing molecular length, the distributions for both types of melts become broader and, of course, their peaks are shifted to larger times. The applied flow has a strong effect on the distributions causing their width and average value to decrease.
3.4. Topological Analysis
As already mentioned several times so far in our article, the linear PEO-5k melt is weakly entangled, characterized by approximately
Z = 2.5 entanglements per chain. On the other hand, several computational works published recently [
4,
5,
39] have provided convincing evidence that in melts of ring molecules, strong ring–ring threading events develop as a ring molecule penetrates through other rings and, simultaneously, itself is threaded by other rings. In fact, the geometric analysis of Tsalikis et al. [
5] demonstrated that, in many cases, these events involve many rings: for example, one ring can be threaded not by one but by several other rings, and vice versa. The same work also demonstrated that the majority of these ring–ring threading events are short-lived and thus have little influence on ring molecule dynamics. However, it also revealed that a non-negligible fraction of the events are long-lived, thus having a strong effect on ring dynamics. Indeed, accounting for these long-lived topological interactions in the
G(
t) of the ring melt in the same way as reptation theory accounts for entanglements in the
G(
t) of linear polymers could explain the long relaxation modes observed experimentally in melts of polymer rings [
5].
It would therefore be of interest to analyze, even in a preliminary way, how the applied flow modifies the above picture. To do this, atomistic configurations from the simulated PEO melts under equilibrium and nonequilibrium conditions were subjected to a detailed topological analysis following the methodology described in earlier works [
5,
40], with the help of the CReTA algorithm of Tzoumanekas and Theodorou [
41] which reduces atomistic configurations to an ensemble of primitive paths (PPs) in which chains are not allowed to cross each other as the algorithm simultaneously reduces the contour length of each chain. For linear melts, at the end of CReTA, kinks along the PP of each polymer chain are identified as points where a rectilinear segment of one chain meets a rectilinear segment of another chain. These kinks are considered as representing entanglement points in the melt, thus their number
Zk per chain is proportional to the number of underlying entanglements
Z per chain. For ring melts, the same concepts apply but the subsequent identification of threading events between rings requires a more elaborate analyis based on 3D vector calculus. All details regarding the implementation of the CReTA algorithm for melts of polymer rings and the geometric analysis for identifying ring–ring threadings can be found elsewhere [
5,
40].
In
Figure 18, we examine the effect of shear rate on kink number density for the L-5k melt. In particular, we report the average number of kinks per chain normalized with the corresponding value,
, under equilibrium conditions. At low shear rates (
), <
Zk> remains practically unaffected by the flow. However, as the shear rate increases, chains orient in the direction of flow (see
Figure 5a), which causes the average number of kink points per chain <
Zk> to drop. This indicates a significant change of the underlying topological network with the applied flow, exactly as has been reported in other computational studies [
32,
42,
43,
44] for linear PE melts.
Figure 19 presents the effect of shear rate on the statistics of threading events for the corresponding ring system, the R-5k melt. At small shear rates (
), the results are practically indistinguishable from the equilibrium ones. Surprisingly enough, the same holds also at higher shear rates (
= 10 or 100), except for some differences in the statistics of multiple threadings involving three or more rings. Given the large decrease in the number of entanglements per chain in the same regime of shear rates for the linear PEO-5k melt shown in
Figure 18, this is a rather unexpected result. Even more surprising is that the same behavior is observed at extremely high shear rates (
) where rings are considerably deformed and aligned: the population of threading events seems to remain practically unaffected by the flow except for the statistics of multiple threadings which are seen to decrease considerably (or even to disappear completely). To explain this behavior, we calculated the survival (or disengagement) times of those threading events that were found to survive the longest under strong shear rates and their penetration length [
4,
5]. Interestingly enough, we found out that practically all of them are short-lived, and (even more importantly) that they are characterized by very short penetration lengths. More precisely, for shear rates corresponding to
, truly long-lived threadings were never observed. We also found that, in distinct contrast to the situation under equilibrium conditions where about 5.8% of threadings correspond to full penetrations [
5], under flow conditions, full penetrations occur only occasionally. To make this message more clear, in
Figure 20 we show a typical atomistic snapshot of a multi-threading event from the NEMD simulations with the R-5k melt at
where the ring with the green color is simultaneously threaded by three other rings shown in yellow, orange and blue. As we can see, in this multi-threading event only a small portion of the surfaces spanned by the three threading rings crosses through the surface of the threaded green chain. As one can understand, the topological analysis can reveal a wealth of information regarding the interplay between flow, topological interactions and rheological response of ring melts under flow which deserves a more detailed and thorough study in the future.
Our main explanation for the transition from the long survival times observed under equilibrium conditions to the short ones observed under flow conditions is ring chain deformation, which prevents other rings from penetrating it. According to our NEMD simulations, polymer rings under flow get significantly elongated in the direction of flow which, in conjunction with their closed loop structure, increases the degree of chain compactness, thus rendering it too difficult or too improbable for other rings to thread them. What we observe under flow conditions is only a very small degree of ring–ring inter-penetration, as a ring can get inside another ring only infinitesimally. As a result, and in clear contrast to the situation under equilibrium conditions, strong ring–ring threadings under shear flow disappear, implying a very minor or secondary role of the corresponding topological interactions on ring dynamics compared to other effects such as local intramolecular re-arrangements or large conformational fluctuations.