www.nature.com/scientificreports
OPEN
Received: 26 May 2017 Accepted: 17 July 2017 Published: xx xx xxxx
Molecular dynamics for linear polymer melts in bulk and confined systems under shear flow Soowon Cho, Sohdam Jeong , Jun Mo Kim & Chunggi Baig In this work, we analyzed the individual chain dynamics for linear polymer melts under shear flow for bulk and confined systems using atomistic nonequilibrium molecular dynamics simulations of unentangled (C50H102) and slightly entangled (C178H358) polyethylene melts. While a certain similarity appears for the bulk and confined systems for the dynamic mechanisms of polymer chains in response to the imposed flow field, the interfacial chain dynamics near the boundary solid walls in the confined system are significantly different from the corresponding bulk chain dynamics. Detailed molecularlevel analysis of the individual chain motions in a wide range of flow strengths are carried out to characterize the intrinsic molecular mechanisms of the bulk and interfacial chains in three flow regimes (weak, intermediate, and strong). These mechanisms essentially underlie various macroscopic structural and rheological properties of polymer systems, such as the mean-square chain end-to-end distance, probability distribution of the chain end-to-end distance, viscosity, and the first normal stress coefficient. Further analysis based on the mesoscopic Brightness method provides additional structural information about the polymer chains in association with their molecular mechanisms. Polymers undergo a variety of processing conditions in practical polymer processes, such as the plastic extrusion process. It is crucial to understand the structural and dynamical behaviors of polymer molecules under various external conditions to economically manufacture high-quality products in such processes. Accordingly, numerous experimental and theoretical research efforts have explored the fundamental aspects behind the macroscopic rheological responses of dense polymeric fluids1, 2, which have enormously advanced our knowledge and enabled us to predict the material properties of polymers under specific conditions in a variety of practical applications. However, many unresolved rheological issues remain (especially from the microscopic viewpoint) for polymeric materials in bulk or confined systems, e.g., fundamental mechanisms underlying stress overshoot, interfacial slip, and melt instability for polymer melts under shear flow2–8. To systematically control such rheological phenomena, it is essential to comprehend the intrinsic molecular dynamics of individual polymer chains separately in bulk and confined situations and how they compare to each other; such an understanding would greatly help to build general knowledge to accurately capture the physical aspects that underlie such complex macroscopic responses of polymer systems and tune the material properties in response to an arbitrary external flow field. In this work, we performed an in-depth analysis on the fundamental molecular mechanisms and dynamic characteristics of bulk and confined polymer melts under shear flow using atomistic nonequilibrium molecular dynamics (NEMD) simulations of unentangled (C50H102) and weakly entangled (C178H358) linear polyethylene (PE) melts. This work is in addition to various advanced experimental9–13 and numerical14–22 studies to reveal the individual chain dynamics in polymer solutions or melts under an external flow field. This molecular-level information attained by directly tracking down individual chain motions is applied to understand the rheological behaviors of representative mesoscopic and macroscopic structural and dynamical properties in response to the applied flow field in a wide range of flow strengths. We analyze the similarities and differences between the bulk and confined melt systems in the characteristic molecular mechanisms and rheological responses under various flow regimes.
School of Energy and Chemical Engineering, Ulsan National Institute of Science and Technology (UNIST), UNISTgil 50, Eonyang-eup, Ulju-gun, Ulsan, 689-798, South Korea. Correspondence and requests for materials should be addressed to J.M.K. (email: [email protected]) or C.B. (email: [email protected]) SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
1
www.nature.com/scientificreports/
Method
Canonical NEMD simulations of monodispersed unentangled (C50H102) and entangled (C178H358) linear PE melts for bulk and confined systems under shear flow were conducted using the p-SLLOD algorithm23 implemented with a Nosé–Hoover thermostat24, 25. For both bulk and confined systems, we employed the Siepmann-Karaboni-Smit united-atom model26 (wherein the original rigid bond was replaced by a flexible bond with a harmonic spring), which has been most commonly applied to the simulations of PE melts14, 15, 17–19. The set of evolution equations for each system14, 19 was numerically integrated with the reversible reference system propagator algorithm (r-RESPA)27 using multiple time scales: a short time scale (0.47 fs) for bonded (bond-stretching, bond-bending, and bond-torsional) interactions and a long time scale (2.35 fs) for the nonbonded Lennard-Jones interactions, thermostat, and external flow field. Both bulk and confined systems were subjected to a simple shear flow for which only the xy-component of the velocity gradient tensor was non-zero, with x, y, and z in Cartesian coordinates representing the flow, velocity gradient, and neutral directions, respectively. The system conditions for all bulk and confined PE melts studied in this work correspond to a constant temperature T = 450 K and pressure P = 1 atm. Specifically, for bulk systems, the simulations were executed at densities of ρ = 0.743 g/cm3 and ρ = 0.782 g/cm3 for the C50 PE melt and C178 PE melt, respectively. The simulation box dimensions for the bulk systems were set as (93.02 Å × 45.00 Å × 45.00 Å) (x × y × z) with a total of 120 molecules for the C50 PE melt, and [(65.89, 131.78, and 263.55) Å × 65.89 Å × 65.89 Å] (enlarged in the flow (x-)direction depending on the applied shear rate to avoid system-size effects) with a total of 54, 108, and 216 molecules, respectively, for the C178 PE melt. For confined systems where PE melts are confined by the two-layered rigid simple-cubic lattice walls in the velocity gradient (y-)direction, the simulations were carried out at ρ = 0.763 g/cm3 and ρ = 0.789 g/cm3 for the C50 PE melt and C178 PE melt, respectively. The walls were composed of 544 atoms with simulation box dimensions of (93.02 Å × 49.51 Å × 45.00 Å) for the C50 PE melt containing 120 molecules, and 676, 1352, and 2704 atoms with the box dimension of [(65.89, 131.78, and 263.55) Å × 70.51 Å × 65.89 Å] containing 54, 108, and 216 molecules, respectively, for the C178 PE melt. The lattice parameter of the simple cubic wall was set to σw = 5.227 Å (=1.33 σCH2). The surface energy parameter of the wall atoms was set as ɛw/kB = 939 K ( = 20 εCH2/kB), corresponding to a mica surface (~200–400 mJ/m2)28, 29. The wall atoms were fixed in their lattice sites during the simulations. For the confined systems, the shear flow field was generated by moving the upper wall at a constant velocity, V, in the flowing direction and fixing the bottom wall (the readers are referred to the Supplementary Information for further methodological details). A wide range of flow strengths, from the linear to the highly nonlinear viscoelastic regimes, was applied to the bulk and confined C50 and C178 PE melts: 0.02 ≤ Wi ≤ 200 for C50 bulk PE, 0.05 ≤ Wi ≤ 500 for C50 confined PE, 0.39 ≤ Wi ≤ 7000 for C178 bulk PE, and 0.68 ≤ Wi ≤ 5600 for C178 confined PE. (The Weissenberg number (Wi) is defined as the product of the longest relaxation time (τ) of the system and the applied shear rate (γ ), Wi ≡ τγ .) The characteristic relaxation time, τ, for each system evaluated by the integral below the stretched-exponential curve describing the decay of the time autocorrelation function of the unit chain end-to-end vector was τ = 0.5 ± 0.05 ns for C50 bulk PE, τ = 1.2 ± 0.1 ns for C50 confined PE, τ = 15.6 ± 1.0 ns for C178 bulk PE, and τ = 26.7 ± 2.0 ns for C178 confined PE work. (Before collecting data, each system was fully equilibrated for a long time (i.e., several times longer than its longest relaxation time τ) at its target state point. After equilibration, the production run of more than 8–10 times longer than τ was carried out to evaluate statistically-reliable physical properties for each system; i.e., 5 ns and 15 ns for the C50 bulk and confined PE melts, respectively, and 100 ns and 200 ns for the C178 bulk and confined PE melts, respectively).
Results and Discussion
Figure 1A illustrates the bulk and confined systems studied in this work. For confined systems, if the center-of-mass position of a chain is located within 2.5 σ (σ = 3.93 Å) from the wall surfaces, it is considered as an interfacial chain. Under equilibrium conditions, in contrast to isotropic random-coil configurations displayed by chains in the bulk system, interfacial chains of confined systems are shown to possess mostly “L” or “U”-shaped (rather extended) configurations on the wall surfaces via the combined effects of the intramolecular chain conformational entropy and the favorable energetic interactions between the polymer and wall. That is, some parts of an interfacial chain are attached to the wall (i.e., adsorbed) and others are detached from the wall (i.e., non-adsorbed). While the adsorbed chain segments experience the wall friction through direct interactions with the wall, the non-adsorbed segments experience intermolecular interactions with nearby surrounding chains in the bulk region. This dual interaction induces distinctive molecular mechanisms for interfacial chains at a given flow strength under shear from those for chains in bulk system or bulk region in confined system. First, for the bulk system (left panel of Fig. 1B), chains are mainly aligned to the flow (x-)direction at low shear rates without a significant structural deformation because the orientation is easier than stretching in response to the applied rotational shear field. At intermediate flow strengths, chains are substantially deformed (stretched) and nearly aligned to the flow direction, and furthermore begin to execute a whole-chain rotation and tumbling motion. In this flow regime, bulk chains mostly exhibit rather symmetrical S-shaped rotations and tumbling behaviors (Fig. 1B). This indicates that the head and tail portions of a chain almost equally (i.e., symmetrically) move in opposite directions relative to each other, as a result of the relative difference in their streaming velocity according to their different (higher and lower) positions in the velocity gradient direction of the applied shear field. In contrast, at high shear rates, chains mainly exhibit tumbling behaviors with hairpin-like configurations, rather than the symmetrical S-shaped configuration. This hairpin-like rotational characteristic occurs because a strong flow field does not allow a sufficient time for the head and tail portions of chains to symmetrically execute their respective movements during the rotational time span; i.e., either head or tail portion alone leads the overall chain tumbling motion quickly without waiting for the other portion to move in the opposite flowing direction. In comparison, for confined system, chains near the wall exhibit distinctive characteristic molecular mechanisms with respect to the flow strength. In the weak flow regime, chains perform a z-to-x rotation (i.e., rotation SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
2
www.nature.com/scientificreports/
Figure 1. Characteristic molecular mechanisms and rheological properties. (A) Schematic illustration for the simulated bulk and confined polyethylene (PE) melt systems at equilibrium. For confined systems, the interfacial region is defined as the region located vertically within 2.5 σ (σ = 3.93 Å) of the wall surfaces. (B) Characteristic molecular dynamic mechanisms of polymer chains for the bulk system and the interfacial region of the confined system under shear in the three respective (weak, intermediate, and strong) flow regimes. Comparisons between the bulk and confined systems for (C) shear viscosity η and (D) the first normal stress coefficient (Ψ1) as a function of Wi for the C50 PE melt under shear flow. The vertical dotted lines distinguish the three characteristic flow regimes for bulk (black) and confined (green) systems. The three characteristic flow regimes are defined, based on the variation of the mean-square end-to-end distance for the bulk systems and the variation of the degree of slip (referring to ref. 19 for details) for the confined systems. The error bars are smaller than the size of the symbols, unless otherwise specified. See Supplementary Fig. 1 for the corresponding plot of the confined system with the Wi number based on the real shear rate accounting for a non-zero slip at the wall.
SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
3
www.nature.com/scientificreports/ from the z-direction to the x-direction) with still residing at the interface through the attractive polymer-wall interactions19. The alignment of interfacial chains with the flow (x-)direction via this in-plane chain rotation reduces the wall friction experienced by the chains during their movement along the flow direction through a decrease in the effective collision area between the chain and wall19. In the intermediate flow regime, the majority of interfacial chains display the out-of-plane wagging mechanism19, wherein outer parts of interfacial chains exhibit a repetitive motion between the detachment from and attachment to the wall (leading to an “L”-shaped configuration in the x-y plane) via the competitive effects of the applied flow field (inducing detachment) and the attractive polymer-wall interaction (enhancing attachment). In this flow regime, chains outside the interfacial region are substantially deformed and aligned to the flow direction, which significantly affects the degree of interaction between interfacial and nearby surrounding bulk chains. With further increasing flow strength, interfacial chains experience strong dynamical collisions with the wall atoms, which gives rise to highly nonlinear rotational dynamics, i.e., irregular (chaotic) rotation and tumbling mechanisms at the wall, leading the interfacial chains to detach from the wall toward the bulk region19. Interfacial chains can only execute a hairpin-like tumbling motion while respecting the geometrical constraint imposed by the wall. All these characteristic molecular mechanisms and dynamics underlie the macroscopic structural and dynamical properties for bulk and confined systems, respectively. Figure 1C presents the viscosity variation with respect to the Wi number for the bulk and confined C50 PE melt systems. Both bulk and confined PE melts exhibit a typical shear-thinning behavior, and the degrees of shear thinning are similar. However, the viscosity of the confined system is larger than that of the bulk system in the entire flow regime, which is consistent with experiments9. This is attributed to the high degree of momentum transfer via collisions (friction) between the interfacial chains and the wall. A similar behavior is shown for the first normal stress coefficient in Fig. 1D; this is again ascribed to the large contribution of interfacial chains to the overall system elasticity, arising from their highly oriented and stretched conformations. In Fig. 2A, we compare the mean-square chain end-to-end distances 〈Rete2〉 as a function of shear rate for the bulk and confined systems. For both the bulk and confined systems, the overall behavior of 〈Rete2〉 can be characterized by three distinct flow regimes. In the weak flow regime, 〈Rete2〉 displays a slight increase with shear rate, indicating that chains in this regime are mostly oriented to the flow direction without a significant structural deformation or stretch. In the intermediate regime, 〈Rete2〉 exhibits a dramatic increase with the applied flow strength and eventually reaches a maximum. With a further increase of the shear rate, 〈Rete2〉 shows a rather decreasing behavior, which was also observed in previous studies14, 15, 18 for bulk PE melts under shear flow. The decrease is ascribed to strong intermolecular collisions together with intense chain rotation and tumbling dynamics at high shear rates. However, we should keep in mind that the characteristic molecular mechanism behind the variation of 〈Rete2〉 for each flow regime is not the same for the bulk and interfacial chains. The 〈Rete2〉 value is quantitatively quite similar for the bulk and confined systems with respect to Wi (a further quantitative similarity is seen for the confined system in terms of Wi number based on the real shear rate, which accounts for a non-zero slip at the wall (Supplementary Fig. 2)). However, the average value of 〈Rete2〉 for only the interfacial chains appears to be somewhat larger than that of the whole confined system, indicative of a larger deformation of the interfacial chains (in association with a higher degree of molecular interactions with the wall) via the favorable polymer-wall interactions. There is also a relatively steeper increase of 〈Rete2〉 in the intermediate flow regime and a larger decrease of 〈Rete2〉 in the strong flow regime for interfacial chains than that for the bulk system. This phenomenon can be understood based on the molecular mechanisms of interfacial chains described in Fig. 1B. At flow strengths greater than that in the weak flow regime (where the characteristic molecular mechanism is the z-to-x in-plane chain rotation without a substantial structural deformation, i.e., only a slight variation of 〈Rete2〉), the chains exhibit the out-of-plane wagging mechanism with highly deformed structures (i.e., a large variation of 〈Rete2〉) in response to the applied flow. The steep increase of 〈Rete2〉 for the interfacial chains in the intermediate flow regime is associated with the repetitive chain attachment-detachment mechanism, because the interfacial chains mostly tend to be aligned and stretched in the flow direction without executing rotational or tumbling dynamics. In the strong flow regime, interfacial chains exhibit a chaotic (irregular) rotation and tumbling mechanism with strong dynamic collisions with the wall atoms; this considerably reduces the stretched chain conformation to a rather compact structure, leading to a significant decrease of 〈Rete2〉 with increasing flow strength. These stronger variations of 〈Rete2〉 for interfacial chains can be further understood by analyzing the probability distribution function of the chain end-to-end distance (Fig. 2B). Compared to the chains in the bulk system, the interfacial chains in the confined system exhibit a more pronounced stretch peak in the intermediate flow regime and a more distinctive rotation peak in the strong flow regime. In addition, the region between the rotation and stretch peaks displays a more pronounced curvature for interfacial chains compared to that for the corresponding bulk system. The interfacial chains exhibit two distinct peaks, even at equilibrium, indicative of a certain amount of chains with a rather extended conformation because of the energetically favorable polymer-wall interaction. Further detailed structural information for the polymer chains under shear flow can be obtained from the Brightness method15, 30, which categorizes the mesoscale chain structures into several configuration classes based on the monomer distribution along the chain. Figure 2C presents the probability distribution function (PDF) for six representative configurations (Coil, Fold, Kink, Dumbbell, Half-dumbbell, and Stretched) for the C50 PE melt system. (It is noted that the Brightness method focuses on the overall chain configuration or shape without regard to the actual chain size). We note that rather short C50 PE chains produce higher portions for the Fold configuration than the Coil configuration at equilibrium. At low shear rates, the portions of Half-dumbbell and Stretched configurations appear somewhat larger for interfacial chains than those for the bulk system, indicative of more extended conformations of interfacial chains via attractive polymer-wall interactions. More interesting features are exposed in the intermediate flow regime where, with increasing shear rate, (i) the Fold portion exhibits a gradual decrease for the bulk system but a small increase for interfacial chains, (ii) the Half-dumbbell portion SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
4
www.nature.com/scientificreports/
Figure 2. Comparison of the mesoscopic and macroscopic structural properties of the bulk and confined systems. (A) Mean-square chain end-to-end distance 〈Rete2〉 as a function of Wi for the bulk and confined C50 PE melt systems. “Interfacial” represents the corresponding result for only the interfacial chains in the confined system. The interfacial regions of both top and bottom walls for the confined systems are found to produce practically identical results (within the statistical uncertainties) for all the microscopic and macroscopic structural and dynamical properties, e.g., 〈Rete2〉 and its PDF, molecular conformations via the Brightness method, interfacial residence time, and so forth. The error bars are smaller than the size of the symbols unless otherwise specified. Comparison of the bulk system and the interfacial region in the confined system for the C50 ∼ PE melt for (B) the probability distribution function, P( Rete ), of the chain end-to-end distance (normalized by ∼ the equilibrium value ( Rete = Rete / Rete eq ) for each system), and (C) the probability distribution function (PDF) of the six representative mesoscale chain configurations (Coil, Fold, Kink, Dumbbell, Half-dumbbell, and Stretched) as a function of Wi computed by the Brightness method. exhibits a small increase for the bulk system but a small decrease for interfacial chains, and (iii) the Stretched portion is quite larger for interfacial chains than for the bulk system. These results can be understood by considering the characteristic molecular dynamics corresponding to the S-shaped tumbling mechanism for the bulk system and the L-shaped repetitive chain detachment-attachment for interfacial chains. In the strong flow regime, the stretched portion appears to decrease more significantly for interfacial chains compared to that for the bulk system, as associated with irregular (chaotic) chain rotation and tumbling mechanisms of interfacial chains. In SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
5
www.nature.com/scientificreports/
Figure 3. Comparisons of the order parameter, chain orientational distribution, and entanglement density distribution for the bulk and confined systems. (A) Chain order parameter λ, which corresponds to the largest eigenvalue of the order tensor S = (3uu − I) /2, where u denotes the unit chain end-to-end vector and I denotes the second-rank unit tensor, as a function of Wi for the bulk and confined C50 PE melt systems. “Interfacial” represents the corresponding result for only the interfacial chains in the confined system. The error bars are smaller than the size of the symbols, unless otherwise specified. See Supplementary Fig. 3 for the corresponding plot of the confined system with Wi number based on the real shear rate accounting for a nonzero slip at the wall. (B) Probability distribution function (PDF) of the chain orientation angle (based on the chain end-to-end vector) with respect to the flow direction as a function of Wi with a schematic illustration of the molecular mechanisms in conjunction with the PDF in the three respective flow regimes. The total space is divided by six angular regions for∼the bulk system and the interfacial region of the confined system. (C) PDF of the local entanglement density (Zes = Zes/σ 3) along the velocity gradient (y-)direction for the entangled C178 PE melts of the bulk and confined systems. Note that Wi = 0 corresponds to the equilibrium condition, Wi = 65 for the bulk system and Wi = 110 for the confined system are in the intermediate flow regime, and Wi = 650 for bulk system and Wi = 570 for the confined system are in the strong flow regime. addition, the dominant hairpin-like chain tumbling mechanism for the bulk system at high shear rates leads to a relatively larger increase in the Fold portion than that for the interfacial chains. Figure 3A presents the variation of the chain order parameter, λ, for the C50 PE melt with respect to the shear rate, which measures the degree of chain alignment in the flow direction in response to the applied field. The interfacial chains already exhibit a large degree of chain ordering in the weak flow regime, with a gradual increase of λ with the shear rate and almost saturation of chain ordering at the end. As such, the λ-value exhibits only a slight variation with shear rate in the intermediate flow regime. In contrast, the corresponding bulk system displays a rather small increase in λ with respect to the shear rate in the weak flow regime and a steep increase in the intermediate flow regime, followed by a plateau region in the strong flow regime. We note that qualitatively similar behavior is observed for the confined system when λ is calculated over all the chains of bulk and interfacial regions; thus, the order parameter averaged over the entire confined system (which might be the case in typical experiments) may provide erroneous structural information for interfacial chains. Furthermore, in contrast to the plateau values of λ for bulk chains in the strong flow regime, the λ-value for the interfacial chains decreases considerably with increasing shear rate in this regime. This phenomenon is directly associated with the irregular (chaotic) chain rotation and tumbling mechanisms via strong molecular collisions of interfacial chains with the wall. Further information related to the chain orientation was obtained by investigating the PDF of the chain orientation angle with respect to the flow direction for the C50 PE system; for this purpose, six spatial orientation regions were chosen, as depicted in Fig. 3B. In the weak flow regime, with increasing shear rate, the PDFs for the bulk system show a decrease for Regions 1, 2, and 6, but a somewhat increase for Regions 4 and 5. As a noticeable
SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
6
www.nature.com/scientificreports/ difference between the bulk and confined systems, the PDFs for interfacial chains exhibit a significant increase for Region 3 and 4 compared to those of the corresponding bulk system; this result directly shows a significant chain alignment to the flow direction for the interfacial chains at low shear rates, which is consistent with the result for the order parameter in Fig. 3A. In the intermediate flow regime, the bulk system shows a steep increase in the PDF of Region 4 but a decrease in the PDF of Region 5 with increasing shear rate. In addition, there is a slight increase in the PDF of Region 3, which is associated with an increase in the degree of the chain tumbling dynamics. In contrast, the interfacial chains exhibit a less steep increase in the PDF of Region 4 and a nearly constant value for the PDF of Region 3 in this flow regime, as can be understood by considering the out-of-wagging (repetitive chain detachment-attachment) mechanism. In the strong flow regime, the bulk system exhibits a slight increasing behavior for the PDFs of Regions 3 and 4 and a slight decreasing (or a plateau) behavior for those of the other regions. In contrast, the interfacial chains exhibit somewhat opposite trends, which are directly related to the intense irregular (chaotic) chain rotation and tumbling mechanisms of interfacial chains on the wall (see Supplementary Fig. 4 for the average tumbling time of chains in the bulk and confined C50 PE melts). Figure 3C presents the spatial distributions along the velocity gradient (y-)direction for the entanglement ∼ density Zes for the entangled C178 PE melts, which were obtained from the topological analysis of the entanglement network of the system using the Z1-code31, 32. (The average number of kinks (intermolecular entanglements) per chain directly obtained from the Z1-code analysis for the C178PE melt system is approximately equal to 6, which is roughly two times larger than the number of entanglements based on the experimental plateau modulus. Our various analyses show that the rheological characteristics of entanglement network for the C178PE melt under shear in a wide range of flow strengths is qualitatively very similar to that of the longer C400 PE melt for both bulk and confined systems. We therefore consider the present result for the rheological aspects of entanglement network for the C178 PE melt to be qualitatively valid for the practical entangled systems). In contrast to the ∼ ∼ homogeneous distributions of Zes for the bulk system, Zes for the confined system exhibits a distinctive shoulder around each of the interfacial regions at equilibrium. This is mainly associated with a relatively higher polymer ∼ density19 in the interfacial regions via the polymer-wall interaction. In addition, Zes is enhanced by entanglements between the detached chain segments of the L or U-shaped interfacial chains and the nearby surrounding bulk chains. This feature can be associated with a relatively higher viscosity (and thus a higher viscous dissipation ∼ under flowing conditions) in the region near the interface. At high shear rates, the Zes-value decreases throughout the confined system due to disentanglement between chains via chain alignment and stretching. Further, the distribution becomes flattened with less pronounced shoulders with their positions appearing somewhat shifted toward the system center. This phenomenon can be understood by considering the frequent movement of interfacial chains into the bulk region and mixing with bulk chains via strong dynamical collisions with the wall at high flow strengths.
Conclusion
Through a detailed analysis of individual chain dynamics using atomistic NEMD simulations for unentangled (C50H102) and weakly entangled (C178H358) linear PE melts under shear flow in bulk and confined systems, we revealed and contrasted the characteristic molecular mechanisms with respect to the applied flow strength between the bulk and interfacial chains. This molecular dynamic information is very useful in understanding the structural and rheological behavior of bulk and confined systems under shear as a function of flow strength, and would be further beneficial for building our general knowledge of predicting and controlling the material properties of the polymer under various flow conditions. The main features identified in this study are summarized below. • Under equilibrium conditions, while polymer chains in the bulk system display random-coil configurations, the interfacial chains of the confined system possess “L” or “U”-shaped configurations on the wall via the combined effects of the intramolecular entropy and the attractive polymer-wall interaction. The detached chain segments of the “L” or “U”-shaped interfacial chains make entanglements with nearby surrounding bulk chains, enhancing the degree of chain entanglement (and thus the viscosity) around the interfacial regions. • In the weak flow regime, while both bulk and interfacial chains are aligned to the flow (x-)direction without a significant structural deformation, the interfacial chains of the confined system perform the z-to-x in-plane rotation at the wall. In addition, compared to the bulk chains, the interfacial chains achieve a much larger degree of chain ordering in the weak flow regime. • In the intermediate flow regime, both bulk and interfacial chains become nearly aligned to the flow direction with a highly deformed (extended) structure; however, in comparison to bulk chains, interfacial chains display a considerably steeper increase of 〈Rete2〉 with respect to the shear rate and a more pronounced stretch ∼ peak in P( Rete ). As regards the characteristic molecular mechanism, the interfacial chains reveal the L-shaped out-of-plane wagging (repetitive chain detachment-attachment) mechanism while chains in the bulk system exhibit a rather symmetrical S-shaped rotation and tumbling behavior. These distinct dynamic characteristics between the bulk and interfacial chains lead to significantly different results for the probability distribution with respect to the mesoscopic chain configurations in the Brightness method and chain orientation angle as a function of the flow strength in this regime. • In the strong flow regime, both the bulk and interfacial chains exhibit chain end-over-end tumbling behaviors. Specifically, bulk chains exhibit a tumbling behavior mainly with a hairpin-like configuration instead of the symmetrical S-shaped one. In comparison, interfacial chains exhibit highly nonlinear rotational dynamics, such as irregular (chaotic) rotation and tumbling (strictly with a hairpin-like configuration) mechanisms at the wall, via strong dynamical collisions with the wall. This interfacial dynamics results in (i) a significant SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
7
www.nature.com/scientificreports/ ∼ decrease of 〈Rete2〉, (ii) distinctive rotation peak in P( Rete ), and (iii) distinct behavior for the probability distribution of the mesoscopic chain configurations in the Brightness method, with respect to the applied flow strength in this regime. In addition, such strong collisions lead the chains to frequently move out of the interfacial region toward the bulk region.
The findings in this study should be carefully taken into account in theoretical modeling. For example, a naïve use of the order parameter and the mean-square chain end-to-end distance averaged over the whole system may lead to erroneous predictions of rheological properties (e.g., stress, anisotropic diffusion, flow birefringence) for confined systems, since the intrinsic structural and dynamical behaviors of the bulk and interfacial chains in response to the applied flow are quite dissimilar. In addition, an adequate description for the flow-induced crystallization of confined systems would require combined information of the distinct characteristic molecular dynamics in each flow regime for both the bulk and interfacial chains. Considering the rapid advance in experimental technique (e.g., fluorescent video microscopy10–13), we also expect the present findings to be potentially useful in the experimental analysis and practical applications of bulk and confined dense polymeric materials undergoing shear flow in the near future.
References
1. Bird, R. B., Armstrong, R. C. & Hassager, O. Dynamics of Polymeric Liquids, vol. 1: Fluid Mechanics, 2nd ed. (Wiley-Interscience, New York, 1987). 2. Doi, M. & Edwards, S. F. The Theory of Polymer Dynamics. (Oxford University Press, New York, 1986). 3. Tordella, J. P. Capillary flow of molten polyethylene-a photographic study of melt fracture. Trans. Soc. Rheol. 1, 203–212 (1957). 4. Ramamurthy, A. V. Wall slip in viscous fluids and influence of materials of construction. J. Rheol. 30, 337–357 (1986). 5. Denn, M. M. Extrusion instabilities and wall slip. Annu. Rev. Fluid Mech. 33, 265–287 (2001). 6. Hatzikiriakos, S. G. Wall slip of molten polymers. Prog. Polym. Sci. 37, 624–643 (2012). 7. Priezjev, N. V. & Troian, S. M. Molecular origin and dynamic behavior of slip in sheared polymer films. Phys. Rev. Lett. 92, 018302 (2004). 8. Jabbarzadeh, A., Atkinson, J. D. & Tanner, R. I. Effect of the wall roughness on slip and rheological properties of hexadecane in molecular dynamics simulation of Couette shear flow between two sinusoidal walls. Phys. Rev. E 61, 690–699 (2000). 9. Van Alsten, J. & Granick, S. Molecular tribometry of ultrathin liquid films. Phys. Rev. Lett. 61, 2570–2573 (1988). 10. Smith, D. E., Babcock, H. P. & Chu, S. Single-polymer dynamics in steady shear flow. Science 283, 1724–1727 (1999). 11. Hur, J. S., Shaqfeh, E. S. G., Babcock, H. P., Smith, D. E. & Chu, S. Dynamics of dilute and semidilute DNA solutions in the start-up of shear flow. J. Rheol. 45, 421–450 (2001). 12. Schroeder, C. M., Teixeira, R. E., Shaqfeh, E. S. G. & Chu, S. Characteristic periodic motion of polymers in shear flow. Phys. Rev. Lett. 95, 018301 (2005). 13. Teixeira, R. E., Babcock, H. P., Shaqfeh, E. S. G. & Chu, S. Shear thinning and tumbling dynamics of single polymers in the flowgradient plane. Macromolecules 38, 581–592 (2005). 14. Baig, C., Mavrantzas, V. G. & Kröger, M. Flow Effects on Melt Structure and Entanglement Network of Linear Polymers: Results from a Nonequilibrium Molecular Dynamics Simulation Study of a Polyethylene Melt in Steady Shear. Macromolecules 43, 6886–6902 (2010). 15. Kim, J. M., Edwards, B. J., Keffer, D. J. & Khomami, B. Dynamics of individual molecules of linear polyethylene liquids under shear: Atomistic simulation and comparison with a free-drainging bead-rod chain. J. Rheol. 54, 283–310 (2010). 16. Dalal, I. S., Hoda, N. & Larson, R. G. Multiple regimes of deformation in shearing flow of isolated polymers. J. Rheol. 56, 305–332 (2012). 17. Kim, J. M. & Baig, C. Precise analysis of polymer rotational dynamics. Sci. Rep. 6, 19127 (2016). 18. Yoon, J., Kim, J. & Baig, C. Nonequilibrium molecular dynamics study of ring polymer melts under shear and elongation flows: A comparison with their linear analogs. J. Rheol. 60, 673–685 (2016). 19. Jeong, S., Cho, S., Kim, J. M. & Baig, C. Molecular mechanisms of interfacial slip for polymer melts under shear flow. J. Rheol. 61, 253–264 (2017). 20. Kim, J. M., Stephanou, P. S., Edwards, B. J. & Khomami, B. A mean-field anisotropic diffusion model for unentangled polymeric liquids and semi-dilute solutions: Model development and comparison with experimental and simulation data. J. Non-Newtonian Fluid Mech. 166, 593–606 (2011). 21. Nafar Sefiddashti, M. H., Edwards, B. J. & Khomami, B. Steady shearing flow of a moderately entangled polyethylene liquid. J. Rheol. 60, 1227–1244 (2016). 22. Jeong, S. H., Kim, J. M. & Baig, C. Rheological influence of short-chain branching for polymeric materials under shear with variable branch density and branching architecture. Macromolecules 50, 4491–4500 (2017). 23. Baig, C., Edwards, B. J., Keffer, D. J. & Cochran, H. D. A proper approach for nonequilibrium molecular dynamics simulations of planar elongational flow. J. Chem. Phys. 122, 114103 (2005). 24. Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 52, 255–268 (1984). 25. Hoover, W. G. Canonical dynamcis: equilibrium phase-space distributions. Phys. Rev. A 31, 1695–1697 (1985). 26. Siepmann, J. I., Karaborni, S. & Smit, B. Simulating the critical behaviour of complex fluids. Nature 365, 330–332 (1993). 27. Tuckerman, M., Berne, B. J. & Martyna, G. J. Reversible multiple time scale molecular dynamics. J. Chem. Phys. 97, 1990–2001 (1992). 28. Granick, S. Motions and Relaxations of Confined Liquids. Science 253, 1374–1379 (1991). 29. Cui, S. T., Cummings, P. T. & Cochran, H. D. Molecular simulation of the transition from liquidlike to solidlike behavior in complex fluids confined to nanoscale gaps. J. Chem. Phys. 114, 7189–7195 (2001). 30. Venkataramani, V., Sureshkumar, R. & Khomami, B. Coarse-grained modeling of macromolecular solutions using a configurationbased approach. J. Rheol. 52, 1143–1177 (2008). 31. Kröger, M. Shortest multiple disconnected path for the analysis of entanglements in two- and three-dimensional polymeric systems. Comput. Phys. Commun. 168, 209–232 (2005). 32. Jeong, S. H., Kim, J. M., Yoon, J., Tzoumanekas, C., Kröger, M. & Baig, C. Influence of molecular architecture on the entanglement network: topological analysis of linear, long- and short-chain branched polyethylene melts via Monte Carlo simulations. Soft Matter 12, 3770–3786 (2016).
Acknowledgements
This work was supported by the National Research Foundation of Korea (2016R1D1A1A09916994, 10067706, 2016R1C1B1014515, and 2016H1A2A1907051). The simulations in this work were carried out on the computational resources of the UNIST Supercomputing Center. SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
8
www.nature.com/scientificreports/
Author Contributions
C.B. planned and designed the numerical simulations. S.C., S.J. and J.M.K. executed the simulations and measured the properties. S.C., S.J., J.M.K., and C.B. analyzed the data and wrote the manuscript.
Additional Information
Supplementary information accompanies this paper at doi:10.1038/s41598-017-08712-5 Competing Interests: The authors declare that they have no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2017
SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
9
OPEN
Received: 26 May 2017 Accepted: 17 July 2017 Published: xx xx xxxx
Molecular dynamics for linear polymer melts in bulk and confined systems under shear flow Soowon Cho, Sohdam Jeong , Jun Mo Kim & Chunggi Baig In this work, we analyzed the individual chain dynamics for linear polymer melts under shear flow for bulk and confined systems using atomistic nonequilibrium molecular dynamics simulations of unentangled (C50H102) and slightly entangled (C178H358) polyethylene melts. While a certain similarity appears for the bulk and confined systems for the dynamic mechanisms of polymer chains in response to the imposed flow field, the interfacial chain dynamics near the boundary solid walls in the confined system are significantly different from the corresponding bulk chain dynamics. Detailed molecularlevel analysis of the individual chain motions in a wide range of flow strengths are carried out to characterize the intrinsic molecular mechanisms of the bulk and interfacial chains in three flow regimes (weak, intermediate, and strong). These mechanisms essentially underlie various macroscopic structural and rheological properties of polymer systems, such as the mean-square chain end-to-end distance, probability distribution of the chain end-to-end distance, viscosity, and the first normal stress coefficient. Further analysis based on the mesoscopic Brightness method provides additional structural information about the polymer chains in association with their molecular mechanisms. Polymers undergo a variety of processing conditions in practical polymer processes, such as the plastic extrusion process. It is crucial to understand the structural and dynamical behaviors of polymer molecules under various external conditions to economically manufacture high-quality products in such processes. Accordingly, numerous experimental and theoretical research efforts have explored the fundamental aspects behind the macroscopic rheological responses of dense polymeric fluids1, 2, which have enormously advanced our knowledge and enabled us to predict the material properties of polymers under specific conditions in a variety of practical applications. However, many unresolved rheological issues remain (especially from the microscopic viewpoint) for polymeric materials in bulk or confined systems, e.g., fundamental mechanisms underlying stress overshoot, interfacial slip, and melt instability for polymer melts under shear flow2–8. To systematically control such rheological phenomena, it is essential to comprehend the intrinsic molecular dynamics of individual polymer chains separately in bulk and confined situations and how they compare to each other; such an understanding would greatly help to build general knowledge to accurately capture the physical aspects that underlie such complex macroscopic responses of polymer systems and tune the material properties in response to an arbitrary external flow field. In this work, we performed an in-depth analysis on the fundamental molecular mechanisms and dynamic characteristics of bulk and confined polymer melts under shear flow using atomistic nonequilibrium molecular dynamics (NEMD) simulations of unentangled (C50H102) and weakly entangled (C178H358) linear polyethylene (PE) melts. This work is in addition to various advanced experimental9–13 and numerical14–22 studies to reveal the individual chain dynamics in polymer solutions or melts under an external flow field. This molecular-level information attained by directly tracking down individual chain motions is applied to understand the rheological behaviors of representative mesoscopic and macroscopic structural and dynamical properties in response to the applied flow field in a wide range of flow strengths. We analyze the similarities and differences between the bulk and confined melt systems in the characteristic molecular mechanisms and rheological responses under various flow regimes.
School of Energy and Chemical Engineering, Ulsan National Institute of Science and Technology (UNIST), UNISTgil 50, Eonyang-eup, Ulju-gun, Ulsan, 689-798, South Korea. Correspondence and requests for materials should be addressed to J.M.K. (email: [email protected]) or C.B. (email: [email protected]) SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
1
www.nature.com/scientificreports/
Method
Canonical NEMD simulations of monodispersed unentangled (C50H102) and entangled (C178H358) linear PE melts for bulk and confined systems under shear flow were conducted using the p-SLLOD algorithm23 implemented with a Nosé–Hoover thermostat24, 25. For both bulk and confined systems, we employed the Siepmann-Karaboni-Smit united-atom model26 (wherein the original rigid bond was replaced by a flexible bond with a harmonic spring), which has been most commonly applied to the simulations of PE melts14, 15, 17–19. The set of evolution equations for each system14, 19 was numerically integrated with the reversible reference system propagator algorithm (r-RESPA)27 using multiple time scales: a short time scale (0.47 fs) for bonded (bond-stretching, bond-bending, and bond-torsional) interactions and a long time scale (2.35 fs) for the nonbonded Lennard-Jones interactions, thermostat, and external flow field. Both bulk and confined systems were subjected to a simple shear flow for which only the xy-component of the velocity gradient tensor was non-zero, with x, y, and z in Cartesian coordinates representing the flow, velocity gradient, and neutral directions, respectively. The system conditions for all bulk and confined PE melts studied in this work correspond to a constant temperature T = 450 K and pressure P = 1 atm. Specifically, for bulk systems, the simulations were executed at densities of ρ = 0.743 g/cm3 and ρ = 0.782 g/cm3 for the C50 PE melt and C178 PE melt, respectively. The simulation box dimensions for the bulk systems were set as (93.02 Å × 45.00 Å × 45.00 Å) (x × y × z) with a total of 120 molecules for the C50 PE melt, and [(65.89, 131.78, and 263.55) Å × 65.89 Å × 65.89 Å] (enlarged in the flow (x-)direction depending on the applied shear rate to avoid system-size effects) with a total of 54, 108, and 216 molecules, respectively, for the C178 PE melt. For confined systems where PE melts are confined by the two-layered rigid simple-cubic lattice walls in the velocity gradient (y-)direction, the simulations were carried out at ρ = 0.763 g/cm3 and ρ = 0.789 g/cm3 for the C50 PE melt and C178 PE melt, respectively. The walls were composed of 544 atoms with simulation box dimensions of (93.02 Å × 49.51 Å × 45.00 Å) for the C50 PE melt containing 120 molecules, and 676, 1352, and 2704 atoms with the box dimension of [(65.89, 131.78, and 263.55) Å × 70.51 Å × 65.89 Å] containing 54, 108, and 216 molecules, respectively, for the C178 PE melt. The lattice parameter of the simple cubic wall was set to σw = 5.227 Å (=1.33 σCH2). The surface energy parameter of the wall atoms was set as ɛw/kB = 939 K ( = 20 εCH2/kB), corresponding to a mica surface (~200–400 mJ/m2)28, 29. The wall atoms were fixed in their lattice sites during the simulations. For the confined systems, the shear flow field was generated by moving the upper wall at a constant velocity, V, in the flowing direction and fixing the bottom wall (the readers are referred to the Supplementary Information for further methodological details). A wide range of flow strengths, from the linear to the highly nonlinear viscoelastic regimes, was applied to the bulk and confined C50 and C178 PE melts: 0.02 ≤ Wi ≤ 200 for C50 bulk PE, 0.05 ≤ Wi ≤ 500 for C50 confined PE, 0.39 ≤ Wi ≤ 7000 for C178 bulk PE, and 0.68 ≤ Wi ≤ 5600 for C178 confined PE. (The Weissenberg number (Wi) is defined as the product of the longest relaxation time (τ) of the system and the applied shear rate (γ ), Wi ≡ τγ .) The characteristic relaxation time, τ, for each system evaluated by the integral below the stretched-exponential curve describing the decay of the time autocorrelation function of the unit chain end-to-end vector was τ = 0.5 ± 0.05 ns for C50 bulk PE, τ = 1.2 ± 0.1 ns for C50 confined PE, τ = 15.6 ± 1.0 ns for C178 bulk PE, and τ = 26.7 ± 2.0 ns for C178 confined PE work. (Before collecting data, each system was fully equilibrated for a long time (i.e., several times longer than its longest relaxation time τ) at its target state point. After equilibration, the production run of more than 8–10 times longer than τ was carried out to evaluate statistically-reliable physical properties for each system; i.e., 5 ns and 15 ns for the C50 bulk and confined PE melts, respectively, and 100 ns and 200 ns for the C178 bulk and confined PE melts, respectively).
Results and Discussion
Figure 1A illustrates the bulk and confined systems studied in this work. For confined systems, if the center-of-mass position of a chain is located within 2.5 σ (σ = 3.93 Å) from the wall surfaces, it is considered as an interfacial chain. Under equilibrium conditions, in contrast to isotropic random-coil configurations displayed by chains in the bulk system, interfacial chains of confined systems are shown to possess mostly “L” or “U”-shaped (rather extended) configurations on the wall surfaces via the combined effects of the intramolecular chain conformational entropy and the favorable energetic interactions between the polymer and wall. That is, some parts of an interfacial chain are attached to the wall (i.e., adsorbed) and others are detached from the wall (i.e., non-adsorbed). While the adsorbed chain segments experience the wall friction through direct interactions with the wall, the non-adsorbed segments experience intermolecular interactions with nearby surrounding chains in the bulk region. This dual interaction induces distinctive molecular mechanisms for interfacial chains at a given flow strength under shear from those for chains in bulk system or bulk region in confined system. First, for the bulk system (left panel of Fig. 1B), chains are mainly aligned to the flow (x-)direction at low shear rates without a significant structural deformation because the orientation is easier than stretching in response to the applied rotational shear field. At intermediate flow strengths, chains are substantially deformed (stretched) and nearly aligned to the flow direction, and furthermore begin to execute a whole-chain rotation and tumbling motion. In this flow regime, bulk chains mostly exhibit rather symmetrical S-shaped rotations and tumbling behaviors (Fig. 1B). This indicates that the head and tail portions of a chain almost equally (i.e., symmetrically) move in opposite directions relative to each other, as a result of the relative difference in their streaming velocity according to their different (higher and lower) positions in the velocity gradient direction of the applied shear field. In contrast, at high shear rates, chains mainly exhibit tumbling behaviors with hairpin-like configurations, rather than the symmetrical S-shaped configuration. This hairpin-like rotational characteristic occurs because a strong flow field does not allow a sufficient time for the head and tail portions of chains to symmetrically execute their respective movements during the rotational time span; i.e., either head or tail portion alone leads the overall chain tumbling motion quickly without waiting for the other portion to move in the opposite flowing direction. In comparison, for confined system, chains near the wall exhibit distinctive characteristic molecular mechanisms with respect to the flow strength. In the weak flow regime, chains perform a z-to-x rotation (i.e., rotation SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
2
www.nature.com/scientificreports/
Figure 1. Characteristic molecular mechanisms and rheological properties. (A) Schematic illustration for the simulated bulk and confined polyethylene (PE) melt systems at equilibrium. For confined systems, the interfacial region is defined as the region located vertically within 2.5 σ (σ = 3.93 Å) of the wall surfaces. (B) Characteristic molecular dynamic mechanisms of polymer chains for the bulk system and the interfacial region of the confined system under shear in the three respective (weak, intermediate, and strong) flow regimes. Comparisons between the bulk and confined systems for (C) shear viscosity η and (D) the first normal stress coefficient (Ψ1) as a function of Wi for the C50 PE melt under shear flow. The vertical dotted lines distinguish the three characteristic flow regimes for bulk (black) and confined (green) systems. The three characteristic flow regimes are defined, based on the variation of the mean-square end-to-end distance for the bulk systems and the variation of the degree of slip (referring to ref. 19 for details) for the confined systems. The error bars are smaller than the size of the symbols, unless otherwise specified. See Supplementary Fig. 1 for the corresponding plot of the confined system with the Wi number based on the real shear rate accounting for a non-zero slip at the wall.
SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
3
www.nature.com/scientificreports/ from the z-direction to the x-direction) with still residing at the interface through the attractive polymer-wall interactions19. The alignment of interfacial chains with the flow (x-)direction via this in-plane chain rotation reduces the wall friction experienced by the chains during their movement along the flow direction through a decrease in the effective collision area between the chain and wall19. In the intermediate flow regime, the majority of interfacial chains display the out-of-plane wagging mechanism19, wherein outer parts of interfacial chains exhibit a repetitive motion between the detachment from and attachment to the wall (leading to an “L”-shaped configuration in the x-y plane) via the competitive effects of the applied flow field (inducing detachment) and the attractive polymer-wall interaction (enhancing attachment). In this flow regime, chains outside the interfacial region are substantially deformed and aligned to the flow direction, which significantly affects the degree of interaction between interfacial and nearby surrounding bulk chains. With further increasing flow strength, interfacial chains experience strong dynamical collisions with the wall atoms, which gives rise to highly nonlinear rotational dynamics, i.e., irregular (chaotic) rotation and tumbling mechanisms at the wall, leading the interfacial chains to detach from the wall toward the bulk region19. Interfacial chains can only execute a hairpin-like tumbling motion while respecting the geometrical constraint imposed by the wall. All these characteristic molecular mechanisms and dynamics underlie the macroscopic structural and dynamical properties for bulk and confined systems, respectively. Figure 1C presents the viscosity variation with respect to the Wi number for the bulk and confined C50 PE melt systems. Both bulk and confined PE melts exhibit a typical shear-thinning behavior, and the degrees of shear thinning are similar. However, the viscosity of the confined system is larger than that of the bulk system in the entire flow regime, which is consistent with experiments9. This is attributed to the high degree of momentum transfer via collisions (friction) between the interfacial chains and the wall. A similar behavior is shown for the first normal stress coefficient in Fig. 1D; this is again ascribed to the large contribution of interfacial chains to the overall system elasticity, arising from their highly oriented and stretched conformations. In Fig. 2A, we compare the mean-square chain end-to-end distances 〈Rete2〉 as a function of shear rate for the bulk and confined systems. For both the bulk and confined systems, the overall behavior of 〈Rete2〉 can be characterized by three distinct flow regimes. In the weak flow regime, 〈Rete2〉 displays a slight increase with shear rate, indicating that chains in this regime are mostly oriented to the flow direction without a significant structural deformation or stretch. In the intermediate regime, 〈Rete2〉 exhibits a dramatic increase with the applied flow strength and eventually reaches a maximum. With a further increase of the shear rate, 〈Rete2〉 shows a rather decreasing behavior, which was also observed in previous studies14, 15, 18 for bulk PE melts under shear flow. The decrease is ascribed to strong intermolecular collisions together with intense chain rotation and tumbling dynamics at high shear rates. However, we should keep in mind that the characteristic molecular mechanism behind the variation of 〈Rete2〉 for each flow regime is not the same for the bulk and interfacial chains. The 〈Rete2〉 value is quantitatively quite similar for the bulk and confined systems with respect to Wi (a further quantitative similarity is seen for the confined system in terms of Wi number based on the real shear rate, which accounts for a non-zero slip at the wall (Supplementary Fig. 2)). However, the average value of 〈Rete2〉 for only the interfacial chains appears to be somewhat larger than that of the whole confined system, indicative of a larger deformation of the interfacial chains (in association with a higher degree of molecular interactions with the wall) via the favorable polymer-wall interactions. There is also a relatively steeper increase of 〈Rete2〉 in the intermediate flow regime and a larger decrease of 〈Rete2〉 in the strong flow regime for interfacial chains than that for the bulk system. This phenomenon can be understood based on the molecular mechanisms of interfacial chains described in Fig. 1B. At flow strengths greater than that in the weak flow regime (where the characteristic molecular mechanism is the z-to-x in-plane chain rotation without a substantial structural deformation, i.e., only a slight variation of 〈Rete2〉), the chains exhibit the out-of-plane wagging mechanism with highly deformed structures (i.e., a large variation of 〈Rete2〉) in response to the applied flow. The steep increase of 〈Rete2〉 for the interfacial chains in the intermediate flow regime is associated with the repetitive chain attachment-detachment mechanism, because the interfacial chains mostly tend to be aligned and stretched in the flow direction without executing rotational or tumbling dynamics. In the strong flow regime, interfacial chains exhibit a chaotic (irregular) rotation and tumbling mechanism with strong dynamic collisions with the wall atoms; this considerably reduces the stretched chain conformation to a rather compact structure, leading to a significant decrease of 〈Rete2〉 with increasing flow strength. These stronger variations of 〈Rete2〉 for interfacial chains can be further understood by analyzing the probability distribution function of the chain end-to-end distance (Fig. 2B). Compared to the chains in the bulk system, the interfacial chains in the confined system exhibit a more pronounced stretch peak in the intermediate flow regime and a more distinctive rotation peak in the strong flow regime. In addition, the region between the rotation and stretch peaks displays a more pronounced curvature for interfacial chains compared to that for the corresponding bulk system. The interfacial chains exhibit two distinct peaks, even at equilibrium, indicative of a certain amount of chains with a rather extended conformation because of the energetically favorable polymer-wall interaction. Further detailed structural information for the polymer chains under shear flow can be obtained from the Brightness method15, 30, which categorizes the mesoscale chain structures into several configuration classes based on the monomer distribution along the chain. Figure 2C presents the probability distribution function (PDF) for six representative configurations (Coil, Fold, Kink, Dumbbell, Half-dumbbell, and Stretched) for the C50 PE melt system. (It is noted that the Brightness method focuses on the overall chain configuration or shape without regard to the actual chain size). We note that rather short C50 PE chains produce higher portions for the Fold configuration than the Coil configuration at equilibrium. At low shear rates, the portions of Half-dumbbell and Stretched configurations appear somewhat larger for interfacial chains than those for the bulk system, indicative of more extended conformations of interfacial chains via attractive polymer-wall interactions. More interesting features are exposed in the intermediate flow regime where, with increasing shear rate, (i) the Fold portion exhibits a gradual decrease for the bulk system but a small increase for interfacial chains, (ii) the Half-dumbbell portion SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
4
www.nature.com/scientificreports/
Figure 2. Comparison of the mesoscopic and macroscopic structural properties of the bulk and confined systems. (A) Mean-square chain end-to-end distance 〈Rete2〉 as a function of Wi for the bulk and confined C50 PE melt systems. “Interfacial” represents the corresponding result for only the interfacial chains in the confined system. The interfacial regions of both top and bottom walls for the confined systems are found to produce practically identical results (within the statistical uncertainties) for all the microscopic and macroscopic structural and dynamical properties, e.g., 〈Rete2〉 and its PDF, molecular conformations via the Brightness method, interfacial residence time, and so forth. The error bars are smaller than the size of the symbols unless otherwise specified. Comparison of the bulk system and the interfacial region in the confined system for the C50 ∼ PE melt for (B) the probability distribution function, P( Rete ), of the chain end-to-end distance (normalized by ∼ the equilibrium value ( Rete = Rete / Rete eq ) for each system), and (C) the probability distribution function (PDF) of the six representative mesoscale chain configurations (Coil, Fold, Kink, Dumbbell, Half-dumbbell, and Stretched) as a function of Wi computed by the Brightness method. exhibits a small increase for the bulk system but a small decrease for interfacial chains, and (iii) the Stretched portion is quite larger for interfacial chains than for the bulk system. These results can be understood by considering the characteristic molecular dynamics corresponding to the S-shaped tumbling mechanism for the bulk system and the L-shaped repetitive chain detachment-attachment for interfacial chains. In the strong flow regime, the stretched portion appears to decrease more significantly for interfacial chains compared to that for the bulk system, as associated with irregular (chaotic) chain rotation and tumbling mechanisms of interfacial chains. In SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
5
www.nature.com/scientificreports/
Figure 3. Comparisons of the order parameter, chain orientational distribution, and entanglement density distribution for the bulk and confined systems. (A) Chain order parameter λ, which corresponds to the largest eigenvalue of the order tensor S = (3uu − I) /2, where u denotes the unit chain end-to-end vector and I denotes the second-rank unit tensor, as a function of Wi for the bulk and confined C50 PE melt systems. “Interfacial” represents the corresponding result for only the interfacial chains in the confined system. The error bars are smaller than the size of the symbols, unless otherwise specified. See Supplementary Fig. 3 for the corresponding plot of the confined system with Wi number based on the real shear rate accounting for a nonzero slip at the wall. (B) Probability distribution function (PDF) of the chain orientation angle (based on the chain end-to-end vector) with respect to the flow direction as a function of Wi with a schematic illustration of the molecular mechanisms in conjunction with the PDF in the three respective flow regimes. The total space is divided by six angular regions for∼the bulk system and the interfacial region of the confined system. (C) PDF of the local entanglement density (Zes = Zes/σ 3) along the velocity gradient (y-)direction for the entangled C178 PE melts of the bulk and confined systems. Note that Wi = 0 corresponds to the equilibrium condition, Wi = 65 for the bulk system and Wi = 110 for the confined system are in the intermediate flow regime, and Wi = 650 for bulk system and Wi = 570 for the confined system are in the strong flow regime. addition, the dominant hairpin-like chain tumbling mechanism for the bulk system at high shear rates leads to a relatively larger increase in the Fold portion than that for the interfacial chains. Figure 3A presents the variation of the chain order parameter, λ, for the C50 PE melt with respect to the shear rate, which measures the degree of chain alignment in the flow direction in response to the applied field. The interfacial chains already exhibit a large degree of chain ordering in the weak flow regime, with a gradual increase of λ with the shear rate and almost saturation of chain ordering at the end. As such, the λ-value exhibits only a slight variation with shear rate in the intermediate flow regime. In contrast, the corresponding bulk system displays a rather small increase in λ with respect to the shear rate in the weak flow regime and a steep increase in the intermediate flow regime, followed by a plateau region in the strong flow regime. We note that qualitatively similar behavior is observed for the confined system when λ is calculated over all the chains of bulk and interfacial regions; thus, the order parameter averaged over the entire confined system (which might be the case in typical experiments) may provide erroneous structural information for interfacial chains. Furthermore, in contrast to the plateau values of λ for bulk chains in the strong flow regime, the λ-value for the interfacial chains decreases considerably with increasing shear rate in this regime. This phenomenon is directly associated with the irregular (chaotic) chain rotation and tumbling mechanisms via strong molecular collisions of interfacial chains with the wall. Further information related to the chain orientation was obtained by investigating the PDF of the chain orientation angle with respect to the flow direction for the C50 PE system; for this purpose, six spatial orientation regions were chosen, as depicted in Fig. 3B. In the weak flow regime, with increasing shear rate, the PDFs for the bulk system show a decrease for Regions 1, 2, and 6, but a somewhat increase for Regions 4 and 5. As a noticeable
SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
6
www.nature.com/scientificreports/ difference between the bulk and confined systems, the PDFs for interfacial chains exhibit a significant increase for Region 3 and 4 compared to those of the corresponding bulk system; this result directly shows a significant chain alignment to the flow direction for the interfacial chains at low shear rates, which is consistent with the result for the order parameter in Fig. 3A. In the intermediate flow regime, the bulk system shows a steep increase in the PDF of Region 4 but a decrease in the PDF of Region 5 with increasing shear rate. In addition, there is a slight increase in the PDF of Region 3, which is associated with an increase in the degree of the chain tumbling dynamics. In contrast, the interfacial chains exhibit a less steep increase in the PDF of Region 4 and a nearly constant value for the PDF of Region 3 in this flow regime, as can be understood by considering the out-of-wagging (repetitive chain detachment-attachment) mechanism. In the strong flow regime, the bulk system exhibits a slight increasing behavior for the PDFs of Regions 3 and 4 and a slight decreasing (or a plateau) behavior for those of the other regions. In contrast, the interfacial chains exhibit somewhat opposite trends, which are directly related to the intense irregular (chaotic) chain rotation and tumbling mechanisms of interfacial chains on the wall (see Supplementary Fig. 4 for the average tumbling time of chains in the bulk and confined C50 PE melts). Figure 3C presents the spatial distributions along the velocity gradient (y-)direction for the entanglement ∼ density Zes for the entangled C178 PE melts, which were obtained from the topological analysis of the entanglement network of the system using the Z1-code31, 32. (The average number of kinks (intermolecular entanglements) per chain directly obtained from the Z1-code analysis for the C178PE melt system is approximately equal to 6, which is roughly two times larger than the number of entanglements based on the experimental plateau modulus. Our various analyses show that the rheological characteristics of entanglement network for the C178PE melt under shear in a wide range of flow strengths is qualitatively very similar to that of the longer C400 PE melt for both bulk and confined systems. We therefore consider the present result for the rheological aspects of entanglement network for the C178 PE melt to be qualitatively valid for the practical entangled systems). In contrast to the ∼ ∼ homogeneous distributions of Zes for the bulk system, Zes for the confined system exhibits a distinctive shoulder around each of the interfacial regions at equilibrium. This is mainly associated with a relatively higher polymer ∼ density19 in the interfacial regions via the polymer-wall interaction. In addition, Zes is enhanced by entanglements between the detached chain segments of the L or U-shaped interfacial chains and the nearby surrounding bulk chains. This feature can be associated with a relatively higher viscosity (and thus a higher viscous dissipation ∼ under flowing conditions) in the region near the interface. At high shear rates, the Zes-value decreases throughout the confined system due to disentanglement between chains via chain alignment and stretching. Further, the distribution becomes flattened with less pronounced shoulders with their positions appearing somewhat shifted toward the system center. This phenomenon can be understood by considering the frequent movement of interfacial chains into the bulk region and mixing with bulk chains via strong dynamical collisions with the wall at high flow strengths.
Conclusion
Through a detailed analysis of individual chain dynamics using atomistic NEMD simulations for unentangled (C50H102) and weakly entangled (C178H358) linear PE melts under shear flow in bulk and confined systems, we revealed and contrasted the characteristic molecular mechanisms with respect to the applied flow strength between the bulk and interfacial chains. This molecular dynamic information is very useful in understanding the structural and rheological behavior of bulk and confined systems under shear as a function of flow strength, and would be further beneficial for building our general knowledge of predicting and controlling the material properties of the polymer under various flow conditions. The main features identified in this study are summarized below. • Under equilibrium conditions, while polymer chains in the bulk system display random-coil configurations, the interfacial chains of the confined system possess “L” or “U”-shaped configurations on the wall via the combined effects of the intramolecular entropy and the attractive polymer-wall interaction. The detached chain segments of the “L” or “U”-shaped interfacial chains make entanglements with nearby surrounding bulk chains, enhancing the degree of chain entanglement (and thus the viscosity) around the interfacial regions. • In the weak flow regime, while both bulk and interfacial chains are aligned to the flow (x-)direction without a significant structural deformation, the interfacial chains of the confined system perform the z-to-x in-plane rotation at the wall. In addition, compared to the bulk chains, the interfacial chains achieve a much larger degree of chain ordering in the weak flow regime. • In the intermediate flow regime, both bulk and interfacial chains become nearly aligned to the flow direction with a highly deformed (extended) structure; however, in comparison to bulk chains, interfacial chains display a considerably steeper increase of 〈Rete2〉 with respect to the shear rate and a more pronounced stretch ∼ peak in P( Rete ). As regards the characteristic molecular mechanism, the interfacial chains reveal the L-shaped out-of-plane wagging (repetitive chain detachment-attachment) mechanism while chains in the bulk system exhibit a rather symmetrical S-shaped rotation and tumbling behavior. These distinct dynamic characteristics between the bulk and interfacial chains lead to significantly different results for the probability distribution with respect to the mesoscopic chain configurations in the Brightness method and chain orientation angle as a function of the flow strength in this regime. • In the strong flow regime, both the bulk and interfacial chains exhibit chain end-over-end tumbling behaviors. Specifically, bulk chains exhibit a tumbling behavior mainly with a hairpin-like configuration instead of the symmetrical S-shaped one. In comparison, interfacial chains exhibit highly nonlinear rotational dynamics, such as irregular (chaotic) rotation and tumbling (strictly with a hairpin-like configuration) mechanisms at the wall, via strong dynamical collisions with the wall. This interfacial dynamics results in (i) a significant SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
7
www.nature.com/scientificreports/ ∼ decrease of 〈Rete2〉, (ii) distinctive rotation peak in P( Rete ), and (iii) distinct behavior for the probability distribution of the mesoscopic chain configurations in the Brightness method, with respect to the applied flow strength in this regime. In addition, such strong collisions lead the chains to frequently move out of the interfacial region toward the bulk region.
The findings in this study should be carefully taken into account in theoretical modeling. For example, a naïve use of the order parameter and the mean-square chain end-to-end distance averaged over the whole system may lead to erroneous predictions of rheological properties (e.g., stress, anisotropic diffusion, flow birefringence) for confined systems, since the intrinsic structural and dynamical behaviors of the bulk and interfacial chains in response to the applied flow are quite dissimilar. In addition, an adequate description for the flow-induced crystallization of confined systems would require combined information of the distinct characteristic molecular dynamics in each flow regime for both the bulk and interfacial chains. Considering the rapid advance in experimental technique (e.g., fluorescent video microscopy10–13), we also expect the present findings to be potentially useful in the experimental analysis and practical applications of bulk and confined dense polymeric materials undergoing shear flow in the near future.
References
1. Bird, R. B., Armstrong, R. C. & Hassager, O. Dynamics of Polymeric Liquids, vol. 1: Fluid Mechanics, 2nd ed. (Wiley-Interscience, New York, 1987). 2. Doi, M. & Edwards, S. F. The Theory of Polymer Dynamics. (Oxford University Press, New York, 1986). 3. Tordella, J. P. Capillary flow of molten polyethylene-a photographic study of melt fracture. Trans. Soc. Rheol. 1, 203–212 (1957). 4. Ramamurthy, A. V. Wall slip in viscous fluids and influence of materials of construction. J. Rheol. 30, 337–357 (1986). 5. Denn, M. M. Extrusion instabilities and wall slip. Annu. Rev. Fluid Mech. 33, 265–287 (2001). 6. Hatzikiriakos, S. G. Wall slip of molten polymers. Prog. Polym. Sci. 37, 624–643 (2012). 7. Priezjev, N. V. & Troian, S. M. Molecular origin and dynamic behavior of slip in sheared polymer films. Phys. Rev. Lett. 92, 018302 (2004). 8. Jabbarzadeh, A., Atkinson, J. D. & Tanner, R. I. Effect of the wall roughness on slip and rheological properties of hexadecane in molecular dynamics simulation of Couette shear flow between two sinusoidal walls. Phys. Rev. E 61, 690–699 (2000). 9. Van Alsten, J. & Granick, S. Molecular tribometry of ultrathin liquid films. Phys. Rev. Lett. 61, 2570–2573 (1988). 10. Smith, D. E., Babcock, H. P. & Chu, S. Single-polymer dynamics in steady shear flow. Science 283, 1724–1727 (1999). 11. Hur, J. S., Shaqfeh, E. S. G., Babcock, H. P., Smith, D. E. & Chu, S. Dynamics of dilute and semidilute DNA solutions in the start-up of shear flow. J. Rheol. 45, 421–450 (2001). 12. Schroeder, C. M., Teixeira, R. E., Shaqfeh, E. S. G. & Chu, S. Characteristic periodic motion of polymers in shear flow. Phys. Rev. Lett. 95, 018301 (2005). 13. Teixeira, R. E., Babcock, H. P., Shaqfeh, E. S. G. & Chu, S. Shear thinning and tumbling dynamics of single polymers in the flowgradient plane. Macromolecules 38, 581–592 (2005). 14. Baig, C., Mavrantzas, V. G. & Kröger, M. Flow Effects on Melt Structure and Entanglement Network of Linear Polymers: Results from a Nonequilibrium Molecular Dynamics Simulation Study of a Polyethylene Melt in Steady Shear. Macromolecules 43, 6886–6902 (2010). 15. Kim, J. M., Edwards, B. J., Keffer, D. J. & Khomami, B. Dynamics of individual molecules of linear polyethylene liquids under shear: Atomistic simulation and comparison with a free-drainging bead-rod chain. J. Rheol. 54, 283–310 (2010). 16. Dalal, I. S., Hoda, N. & Larson, R. G. Multiple regimes of deformation in shearing flow of isolated polymers. J. Rheol. 56, 305–332 (2012). 17. Kim, J. M. & Baig, C. Precise analysis of polymer rotational dynamics. Sci. Rep. 6, 19127 (2016). 18. Yoon, J., Kim, J. & Baig, C. Nonequilibrium molecular dynamics study of ring polymer melts under shear and elongation flows: A comparison with their linear analogs. J. Rheol. 60, 673–685 (2016). 19. Jeong, S., Cho, S., Kim, J. M. & Baig, C. Molecular mechanisms of interfacial slip for polymer melts under shear flow. J. Rheol. 61, 253–264 (2017). 20. Kim, J. M., Stephanou, P. S., Edwards, B. J. & Khomami, B. A mean-field anisotropic diffusion model for unentangled polymeric liquids and semi-dilute solutions: Model development and comparison with experimental and simulation data. J. Non-Newtonian Fluid Mech. 166, 593–606 (2011). 21. Nafar Sefiddashti, M. H., Edwards, B. J. & Khomami, B. Steady shearing flow of a moderately entangled polyethylene liquid. J. Rheol. 60, 1227–1244 (2016). 22. Jeong, S. H., Kim, J. M. & Baig, C. Rheological influence of short-chain branching for polymeric materials under shear with variable branch density and branching architecture. Macromolecules 50, 4491–4500 (2017). 23. Baig, C., Edwards, B. J., Keffer, D. J. & Cochran, H. D. A proper approach for nonequilibrium molecular dynamics simulations of planar elongational flow. J. Chem. Phys. 122, 114103 (2005). 24. Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 52, 255–268 (1984). 25. Hoover, W. G. Canonical dynamcis: equilibrium phase-space distributions. Phys. Rev. A 31, 1695–1697 (1985). 26. Siepmann, J. I., Karaborni, S. & Smit, B. Simulating the critical behaviour of complex fluids. Nature 365, 330–332 (1993). 27. Tuckerman, M., Berne, B. J. & Martyna, G. J. Reversible multiple time scale molecular dynamics. J. Chem. Phys. 97, 1990–2001 (1992). 28. Granick, S. Motions and Relaxations of Confined Liquids. Science 253, 1374–1379 (1991). 29. Cui, S. T., Cummings, P. T. & Cochran, H. D. Molecular simulation of the transition from liquidlike to solidlike behavior in complex fluids confined to nanoscale gaps. J. Chem. Phys. 114, 7189–7195 (2001). 30. Venkataramani, V., Sureshkumar, R. & Khomami, B. Coarse-grained modeling of macromolecular solutions using a configurationbased approach. J. Rheol. 52, 1143–1177 (2008). 31. Kröger, M. Shortest multiple disconnected path for the analysis of entanglements in two- and three-dimensional polymeric systems. Comput. Phys. Commun. 168, 209–232 (2005). 32. Jeong, S. H., Kim, J. M., Yoon, J., Tzoumanekas, C., Kröger, M. & Baig, C. Influence of molecular architecture on the entanglement network: topological analysis of linear, long- and short-chain branched polyethylene melts via Monte Carlo simulations. Soft Matter 12, 3770–3786 (2016).
Acknowledgements
This work was supported by the National Research Foundation of Korea (2016R1D1A1A09916994, 10067706, 2016R1C1B1014515, and 2016H1A2A1907051). The simulations in this work were carried out on the computational resources of the UNIST Supercomputing Center. SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
8
www.nature.com/scientificreports/
Author Contributions
C.B. planned and designed the numerical simulations. S.C., S.J. and J.M.K. executed the simulations and measured the properties. S.C., S.J., J.M.K., and C.B. analyzed the data and wrote the manuscript.
Additional Information
Supplementary information accompanies this paper at doi:10.1038/s41598-017-08712-5 Competing Interests: The authors declare that they have no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2017
SCIENTIfIC REPOrtS | 7: 9004 | DOI:10.1038/s41598-017-08712-5
9
