Molecular simulation of flow-enhanced nucleation of polyethylene crystallites in biaxial flows

Semicrystalline polymers constitute the largest fraction of industrial plastics. They are used for fabrication of fibers, films, packaging, and composites.¹ Production processes for many of these products, such as injection molding, blown-film extrusion, and heat-sealing, involve complex flows coupled with solidification from the melt. The effect of flow on the final structure and properties, especially crystallinity, is complicated. Crystallization involves rearrangement of molecules into ordered regions; flow aids in this ordering through reorientation of polymer chain segments. The amount of stretching and orientation that the flow introduces dramatically affects the crystallization rate and consequent morphology.² The evidence for rate enhancement in blown-film processes³ and changes in morphology in injection molding⁴ exemplify this effect of flow, commonly referred to as Flow-Induced Crystallization (FIC). The earliest stage of FIC is Flow-Enhanced Nucleation (FEN), which involves small molecular clusters that arise suddenly in the solution or melt⁵ but plays an important role in determining the final morphology of a material.

Many experimental studies have probed the phenomena of FIC and FEN.^6–13 While experiments have provided valuable insights into these phenomena, their observations are typically indirect. By contrast, computational techniques, such as molecular dynamics simulations, are well suited to resolve the fine spatiotemporal scale of a nucleation event.^14–22 However, to date, only the effects on nucleation due to shear or uniaxial extensional flow fields have been studied. Biaxial extensional flow (BEF) fields are industrially relevant for processes such as blown-film extrusion (BFE), but their effects on polymer crystallization have received only sparse attention.^23–30 Bullwinkle et al.²⁹ studied BFE using small angle light scattering. They modeled FIC but did not specifically study the effects of biaxial flow. Van Drongelen et al.³⁰ measured crystallinity and molecular orientation in BFE using wide-angle x-ray diffraction. They found that the process conditions that lead to high molecular orientation also lead to higher crystallinity. Zhang et al.²⁶ and Chen et al.²⁵ studied experimentally the effect of transient biaxial stretching on the crystallization of polyethylene terephthalate (PET) and natural rubber, respectively. Both reports highlighted the importance of stretch ratio and orientation of amorphous chains. Nie et al.²⁷ modeled nucleation kinetics under biaxial stretching using nonequilibrium molecular dynamics (NEMD) simulations and a coarse-grained model of polyvinyl alcohol; they applied a start-up flow followed by quenching below the melting point with cessation of flow and showed that the induction time for nucleation correlates with the average local orientation parameter at the final state of strain. However, the simulation protocol employed there convolves the effects of nucleation with relaxation after cessation of flow.

Nicholson and Rutledge²¹ examined nucleation rates for a united atom model of n-pentacontane (C150) in steady state shear or uniaxial extensional flow using NEMD simulations with lattice reduction techniques for generalized homogeneous flows. They simulated a range of strain rates and employed an iso-Weissenberg number quench to reduce the effect of the quench on chain conformation. They used these results to test several nucleation rate models proposed in the literature. Each model was classified according to the nature of the independent variable in the model: (a) strain rate, (b) stress, (c) conformation, or (d) free energy. Models based on measures of stress or conformation were found to perform significantly better than the others; in some cases, the two were indistinguishable. They argued that the similarity in performance is a consequence of the stress-optical rule and arises from coaxiality of stress and conformation tensors in shear and uniaxial extensional flow.

In this work, we report the effect of steady state biaxial flow on the orientation of C150 chain segments in the melt state, and in turn on nucleation kinetics, using the methods of Nicholson and Rutledge.²¹ Equibiaxial and non-equibiaxial flows are characterized by different ratios of the strain rate in orthogonal directions within the plane of biaxial extension. The orientation of chain segments under these flows is quantified using the conformation tensor and nematic order parameters. Nucleation kinetics under BEF are compared and contrasted with those published previously for uniaxial and shear flows. A model in which FEN is correlated with the P₂ order parameter tensor for Kuhn segments is found to work well under all of the flow conditions studied.

The simulated system consisted of a melt of 100 C150 chains in which each CH₂ or CH₃ unit was modeled as a united atom (UA). The) force field of Paul, Yoon, and Smith (PYS)³¹ with modifications by Waheed et al.^32,33 was used for the simulations; it has been shown previously to provide a realistic description of polyethylene melts and crystal nucleation.^34–37 

For the simulation procedure, a steady state biaxial flow was imposed in the melt state at 440 K. This temperature is above the experimental melting point of 396.4 K for C150.⁴³ The PYS force field has been shown elsewhere to reproduce experimental melting points for n-alkanes within about ±2%.^35,44,45 The melt state flow simulation was run for 90 ns to reach the steady state, as confirmed by fluctuation of the Kuhn segment order parameter (defined below) about a steady average. The simulation was run for an additional 750 ns to analyze the melt flow state. 25 snapshots were sampled at intervals of 30 ns for subsequent nucleation studies. This time interval is slightly greater than the Rouse time (τ_R) of 26 ns³⁶ for C150 at 440 K, ensuring independent, uncorrelated starting configurations. These configurations were then quenched to a crystallization temperature of 280 K, corresponding to a 29% supercooling (⁠$Tm−T/Tm$⁠), to initiate crystal nucleation trajectories. τ_R for C150 at 280 K is estimated to be 340 ns.³⁶ The strain rate of the simulation after quenching was chosen such that the Rouse Weissenberg number, defined as $Wi=γ̇τR$⁠, was the same before and after the quench (an “iso-Weissenberg number quench”). Example strain rates before and after quench are presented in the supplementary material, Table S1. The steady-state NEMD simulations were conducted for strain rate ratios corresponding to Λ = 1, 1.2, 1.5, 3, and 10. Out of these simulations, configurations from simulations having Λ = 1, 1.5, and 3 were used for nucleation simulations, as these were sufficient to draw the conclusions of this work; with further increases in Λ, the deformation approaches planar extensional flow, which is similar to uniaxial extensional flow, for which data already exist.²¹ 

For supercritical crystalline clusters, we characterized the compactness of the cluster using its fractal dimension, D_f. For this purpose, we followed the box counting algorithm described by Gagnepain and Roques-Carmes⁵⁰ and also employed by Nicholson and Rutledge.³⁹ Briefly, starting with a cubic box of side length L that just encloses the cluster, the box was further divided into voxels of side length δ = L/2, L/3, etc. For a given δ, the number of voxels Q(δ) that intersect with any part of a UA in the cluster was counted. D_f was then determined from the slope of −log(Q(δ)) vs log(δ) in the range 1 nm < δ < σ_LJ, where σ_LJ is the Lennard-Jones radius of a UA. D_f was found to be not very sensitive to the range over which δ was averaged (the 1 nm cutoff is chosen for consistency). For a three-dimensional object (such as our nuclei), D_f = 3 for an object that is compact and completely fills space. D_f < 3 for clusters that are not compact and space-filling.

For each simulation, D_f, κ², and n_max, the size of the largest cluster, were each averaged over five consecutive 3 ns windows, and then, <D_f> and <κ²> were plotted against n_max. For each biaxial flow mode (Wi, Λ), the results for 25 simulations were binned, and the mean and standard deviation of each bin were plotted vs <n_max > for that bin.

Nicholson and Rutledge²² showed that two-dimensional distributions of P_2,K,i and ν_K,i split into two distinct peaks above a certain Wi for both uniaxial and shear flows. The peak associated with the lower values of P_2,K,i and ν_K,i was dubbed “flow-oriented,” whereas the peak associated with higher values of P_2,K,i and ν_K,i was dubbed “nematic.” Nematic UAs in the melt are those belonging to the nematic peak in this two-dimensional distribution.

The imposition of flow induces orientation of the polymer chain segments to varying degrees in the melt state at 440 K. The degree of orientation was characterized by the second invariant of the global nematic Kuhn segment order tensor, J₂(P_2,K), which can range between 0 and 1.5. This second invariant is plotted vs Wi in Fig. 1 for the different flow conditions. The uniaxial extensional flow imposes the highest degree of orientation among all flows studied. J₂(P_2,K) is approximately zero at low Wi < 1, indicative of isotropic melts at small strain rates. J₂(P_2,K) begins to increase around Wi = 1, indicative of oriented melts, and then to level off above Wi = 10, as orientation becomes saturated. The same qualitative behavior is observed for the biaxial flows, but with J₂(P_2,K) being lower for the lower strain rate ratio in the range of Λ = 10 to 1, with the last corresponding to equibiaxial flow. The shear flow curve overlaps approximately the curve for biaxial flow at Λ = 1.2. These behaviors show that orientation can vary over a wide range under biaxial flow conditions, from orientations weaker than in shear flow for equibiaxial flow to orientations approaching that in uniaxial extension for strongly non-equibiaxial flows.

To study nucleation kinetics, 25 independent configurations from the melt flow state were quenched to a crystallization temperature T_c = 280 K at constant Wi (Table S1 in the supplementary material shows typical strain rates before and after the quench). Figure 2 shows the size of the largest crystalline cluster as a function of time after quench for a typical trajectory. The first passage time (FPT) of the largest cluster in the system was tracked for each of the 25 nucleation simulations to create a first-passage time distribution, FPTD. Figure 3 shows the mean and standard deviation, τ_FPT and σ_FPT, of the FPT distribution as a function of the size, n, of the largest cluster for both small and large Wi. At small Wi [Fig. 3(a)], σ_FPT levels off at the end of the stochastic nucleation process, whereas τ_FPT continues to increase as growth processes become important for supercritical clusters. For sufficiently large Wi [Fig. 3(b)], neither τ_FPT nor σ_FPT levels off, indicative of a breakdown of classical nucleation theory, which is consistent with previously reported observations.²² For this reason, σ_FPT was evaluated for a fixed range of cluster sizes that was deemed to be only slightly supercritical (150 < n < 300), as described in Sec. II. This range was also chosen for consistency with prior work.²¹ 

Figure 4 shows the nucleation rate as a function of Wi for different flow conditions. The nucleation rate, I_S, increases with increasing Wi for all flows, indicative of FEN. At low Wi < 1, the increase in the nucleation rate is approximately independent of the flow conditions, whereas for Wi > 1, the nucleation rate is strongly dependent on the flow conditions. As was the case for J₂(P_2,K) in Fig. 1, I_S decreases with the decreasing strain rate ratio between 3 and 1 in BEF at any given Wi, from the highest values for strongly non-equibiaxial flow to the lowest values for equibiaxial flow. For the case of equibiaxial flow, Λ = 1, the nucleation rate appears to level off around Wi = 10. This behavior is a preliminary indication that nucleation rate correlates with chain segment orientation. This correlation is explored quantitatively in the discussion of FEN models that follows.

Figure 5 shows the correlation between the second invariants of the extra stress tensor τ = σ − P_extI, where σ is the stress tensor and P_ext is the externally applied pressure (1 atm in our case), for shear, uniaxial, and biaxial flows over a wide range of Wi. The superposition of data for shear and uniaxial flow is indicative of the stress-optical rule and the coaxiality mentioned in the Introduction. Also shown in Fig. 5 are the dependence of these second invariants for the more general case of biaxial extensional flows for three different strain rate ratios. The biaxial flow data do not superpose. At low Wi (where the second invariants are small), the biaxial data are shifted upward relative to the shear and uniaxial data, indicative of elevated extra stresses in biaxial flows. At high Wi, the results for non-equibiaxial flows parallel the uniaxial flow case, whereas the equibiaxial data deviate toward still higher levels of extra stress. Evidently, significant deviations from the stress-optical rule are manifested in biaxial flows at large Wi.

Figure 6(a) shows the nucleation rate for shear, uniaxial, and biaxial flow conditions at three different strain rate ratios, plotted as a function of J₂(P_2,K). Also shown is the conformation-based model for nucleation rate fitted to all the data points shown in the plot. The model parameters found in this study, I_s,q = 32.46 × 10³⁰ m⁻³ s⁻¹ and θ_K = 6.01, are similar to the values reported earlier based on uniaxial and shear data alone (I_s,q = 38.5 × 10³⁰ m⁻³ s⁻¹ and θ_K = 5.79).²¹ The model based on the nematic order tensor, Eq. (18), provides a good description of all the data from NEMD simulations for shear, uniaxial, and biaxial flow conditions.

Figure 6(b) shows the same data for the nucleation rate plotted as a function of J₂(τ). Clearly, the data do not collapse onto a single curve, and no single parameterization can describe all the data. The model based on the extra stress tensor, Eq. (19), using parameters reported by Nicholson and Rutledge,²¹ is also shown. The nucleation rates for the three biaxial flows follow a similar trend as the model at low Wi but deviate at high Wi, indicating a reduced sensitivity of the nucleation rate to extra stress at large applied stresses. This deviation reveals that the extra stress tensor is not as robust as the Kuhn segment order parameter to correlate FEN.

In their experimental study of strain-induced crystallization of natural rubber, Chen et al.²⁵ observed that crystallization did not occur under biaxial stretching conditions for stretch ratios less than 1.6. This suppression of crystallization was attributed to frustration in the alignment of amorphous segments under nearly equibiaxial stretch, which is consistent with the low levels of Kuhn segment orientational order observed under equibiaxial flow conditions for all strain rates (Wi) in the current work (Fig. 1). The low levels of orientational order achieved, despite the application of relatively high Wi (and correspondingly high applied stress, as shown in Fig. S2 of the supplementary material) seems to be indicative of this frustration. As the strain rate ratio increases, a dominant direction of extension emerges and the strain rate required to realize a given amount of orientation decreases.

The post-critical nuclei can be observed to understand the effect of flow on the shape and morphology of these nuclei. However, calculating the shape metrics from a single snapshot can be misleading since short-lived appendages appear and disappear with some frequency. Figure 7 shows examples of nuclei showing all UAs appearing as part of the largest cluster in a 3 ns window. The UAs that have low persistence are more transparent and yellow. Within 3 ns, these nuclei can fluctuate about the same shape [Figs. 7(a) and 7(b)], have short-lived appendages [Figs. 7(c) and 7(d)], or consider two nuclei that are close to each other [Figs. 7(e) and 7(f)]. Thus, averaging shape metrics (as detailed in Sec. II) gives us a better estimate of the persistent part of the nuclei. More details on the averaging and rendering of Fig. 7 are given in the supplementary material.

Fractal dimensions and relative shape anisotropy are plotted as functions of cluster size in Figs. 8 and 9, respectively. The fractal dimension values are all above 2.5, which for a three-dimensional object means that they are largely space filling. The shape anisotropy values are also below 0.5, with most of them being in the 0.2–0.3 range. Thus, the nuclei are mostly spherical and not cylindrical or rod-like. For a given Wi, increasing the size of nucleus decreases the fractal dimension, but the anisotropy is largely unaffected. Thus, the larger nuclei are not as completely filled as the smaller ones and become slightly more spherical as they grow. On the other hand, increasing Wi decreases the fractal dimension of nuclei and increases their shape anisotropy. Thus, slower flows allow nuclei to be more completely filled and more spherical in shape. Slower flows lead to slower formation of nuclei, which gives polymer chains more time to organize themselves in more filled and spherical nuclei. Comparing between Λ = 1 and Λ = 3, there are no significant differences in shape metrics. Thus, we can say that the type of biaxial flow does not affect the shape of nuclei for the flow conditions examined in this work.

Flow biases the orientation of Kuhn segments in the melt, so it may be expected that the orientation of crystallites is similarly biased. The crystallite orientation can be conveniently characterized by the direction cosines of the director of S_cluster with respect to the axes of flow. Figure 10(a) shows the relevant section of a triangle plot of the square of the direction cosines of the cluster director, averaged over clusters of size (or n) approximately equal to 200 from 25 trajectories. To show the variation in orientations of these 25 nuclei, the director is projected onto the plane of extension (xy plane). The angle this projection makes with the x axis is called θ_proj and the mean and standard deviation of θ_proj is shown in Fig. 10(b). The data point for crystallite orientation in the absence of flow (quiescent crystallization) lies close to the centroid of the triangle in Fig. 10(a), showing that crystal stems have no orientation preference in this case. With increasing Wi, all data points for biaxial flows are shifted toward the xy plane. The data points for equibiaxial flow (in blue) are all approximately equidistant from the x and y vertices, but with a very broad distribution of values in the angle θ_proj. Figure 10(b) confirms that there is a large variation in orientation in the xy plane for the case of Λ = 1. Taken together, these observations indicate that the clusters are randomly oriented in the xy plane for all Wi when Λ = 1. For strain rate ratios Λ > 1, the crystallites orientation shifts toward the x axis with increasing Wi, and the distribution of orientation becomes rather narrow. Thus, crystallites in strong biaxial flows with Λ > 1 are oriented along the x axis (i.e., the axis of strongest extension). Some example angle distributions are shown in the supplementary material, Fig. S3.

Following the method of Nicholson and Rutledge,²² the formation of nematic domains in the melts under biaxial flow is examined using two-dimensional distributions of the populations of UAs as functions of the Kuhn segment order parameters ν_K and P_2,K. For Λ > 1, these distributions exhibit two peaks above a critical Wi. This behavior is shown in Fig. 11.

Although the individual order parameters are clearly correlated, projections of the distributions onto the individual order parameter planes reveal that the two distinct peaks are more evident as a function of ν_K than P_2,K. Thus, Fig. 12 shows projections of the UA populations on the ν_K plane for a range of Wi at three different strain rate ratios, Λ = 1.0, 1.5, and 3.0. Figure 12(a) and the first column of Fig. 11 show that equibiaxial flow does not give rise to the nematic order peak, and only a shoulder can be seen at the highest Wi. On the other hand, higher strain rate ratios [Figs. 12(b) and 12(c) and the second and third columns of Fig. 11] show bimodal distributions. The lower valued peak is designated the “flow-oriented” population, and the higher valued peak is designated the “nematic” population.

Figure 13 shows the function $WνK(n)$ as a function of n for several combinations of Wi and Λ. In all cases, one observes $WνK(n)$ > $νK,i,0i$ for all n, with $WνK(n)$ approaching$νK,i,0i$ in the limit that n approaches N, the total number of UAs in the system. This observation confirms that nucleation occurs preferentially in nematic domains as measured by ν_K.

Nie et al.²⁷ performed NEMD simulations of a polyvinyl alcohol (PVA) undergoing biaxial stretch; they also reported a spatial correlation between the more ordered sites in the melt and the sites found in crystal nuclei. Using all-atom MD simulations of polyethylene (C1000) under shear flow, Tang et al.⁵¹ concluded that FEN is a two-step process in which nucleation occurs in regions of higher density. The higher density of these regions was attributed to their being more ordered. Thus, similar to this work, they observe nucleation occurring preferentially from more ordered regions.

Figure 14 shows a phase diagram based on whether a nematic phase was observed in the simulations at 440 K. Each of the biaxial flows (except Λ = 1) shows a minimum Wi for the formation of the nematic phase. Wi at which nematic domains emerge decreases as the strain rate ratio increases, suggesting that such nematic domains are stabilized by more anisotropic, or non-equibiaxial, flow conditions. This idea is supported by looking at extremes of anisotropy in flow. Equibiaxial flow, which is the most isotropic, fails to show nematic phase formation even at the highest Wi studied. On the other hand, uniaxial flow shows nematic phase formation at Wi = 2.31,²² which is lower than for all biaxial flows studied here. Thus, flows with a higher strain rate ratio not only give a higher average Kuhn segment order but also are more efficient in the formation of nematic order domains within the melt.

We performed NEMD simulations to study FEN under biaxial flows. Varying the strain rate ratio (Λ), we observed a variety of orientations in the melt. These orientations range from low orientation in equibiaxial flows to high orientation in anisotropic biaxially flows, asymptotically approaching uniaxial flow for large values of the strain rate ratio. Thus, biaxial extensional flows (BEFs) provide a good candidate to test FEN models. It is found that a model based on the Kuhn segment conformation tensor best explains FEN for all flow conditions. The stress-optical rule fails for BEF, and consequently, the model based on the extra stress tensor does not fit the data. Segmental orientation at the Kuhn segment scale can be obtained by melt birefringence studies, allowing for experimental validation of this model.

An NEMD study allows us to study the post-critical nuclei in detail and observe orientation and shape. We show that flow greatly influences the direction of crystal stems. Equibiaxial flows introduce randomness in orientation within the plane of extension, whereas asymmetric flows tend to orient all nuclei in one direction. We also observe that nuclei are largely spherical even for strong and anisotropic flows. Thus, nuclei do not show evidence of change in shape up to the sizes studied here. Hence, any change in morphology observed in experiments must occur at larger length scales than the ones studied here.

We are also able to show nematic phase formation in melt for biaxial flows that are asymmetric enough. Equibiaxial flows behave differently and fail to show the formation of a nematic phase. This nematic phase is also shown to nucleate preferentially. Thus, we can say that FEN occurs through the formation of more ordered parts in the melt phase. Our observations agree with other FEN and FIC studies conducted in the literature. From these observations, one can say that biaxial flows (especially equibiaxial) behave differently from uniaxial or shear flows previously studied. They orient melt and crystal stems differently and may not form a nematic phase. Thus, these flows must be treated appropriately in FEN experiments, simulations, and models.

Additional information and images are provided as the supplementary material in a separate file.

The authors gratefully acknowledge support from Dow through the University Partnership Initiative. The authors also acknowledge the Lammot du Pont Chair in Chemical Engineering to G.C.R. for funding of this work.

The authors have no conflicts to disclose.

Chinmay S. Gangal: Formal analysis (lead); Investigation (lead); Software (lead). Gregory C. Rutledge: Conceptualization (lead); Supervision (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.