Ordering and Crystallization of Entangled Polyethylene Melts under

Nov 21, 2018 - U.S. Army Research Laboratory, Aberdeen Proving Ground , Maryland 21005 , United States. ‡ SURVICE Engineering Company, Aberdeen ...
0 downloads 0 Views 9MB Size
Article Cite This: Macromolecules XXXX, XXX, XXX−XXX

pubs.acs.org/Macromolecules

Ordering and Crystallization of Entangled Polyethylene Melts under Uniaxial Tension: A Molecular Dynamics Study Yelena R. Sliozberg,†,‡ In-Chul Yeh,*,† Martin Kröger,§ Kevin A. Masser,† Joseph L. Lenhart,† and Jan W. Andzelm*,† †

U.S. Army Research Laboratory, Aberdeen Proving Ground, Maryland 21005, United States SURVICE Engineering Company, Aberdeen Proving Ground, Maryland 21005, United States § Polymer Physics, Department of Materials, ETH Zürich, Leopold-Ruzicka-Weg 4, CH-8093 Zürich, Switzerland Downloaded via KAOHSIUNG MEDICAL UNIV on November 22, 2018 at 06:41:40 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



ABSTRACT: Morphological and mechanical properties of semicrystalline polymers are strongly influenced by flow-induced crystallization during processing. We perform extensive molecular dynamics simulations with more than 1 million atoms to describe orientation, drawability, and crystallization of entangled polyethylene melts under uniaxial tensions at three different strain rates and after a subsequent cooling. During tensile deformation at the lowest strain rate of 107 s−1, the polyethylene melt experiences entanglement loss and crystal nucleation. At higher strain rates of 108 and 109 s−1, we observe a higher degree of chain alignment and void formation in addition to disentanglement and crystal nucleation. Chain segments make sharp turns relative to the neighboring chain orientations at the entanglement points, which manifests as a bimodal distribution of the local order parameter. Upon cooling below the melting temperature, semicrystalline polyethylene with a crystallinity close to 50% is formed. The entanglements are located in the amorphous regions of the semicrystalline polyethylene, with some located in the crystal/amorphous interface region. The chain ends of the semicrystalline polyethylene are preferentially localized at the crystal/amorphous interface, which is in agreement with recent experimental results.

1. INTRODUCTION The properties of high-performance semicrystalline polymeric fibers are greatly affected by crystal growth at the molecular level. For instance, molecular orientation and extension of polymer chains in polyethylene (PE) fibers directly influence the morphology and mechanical performance of the final material through flow-induced crystallization (FIC).1−3 The molecular level structure of these polymeric fibers strongly correlates with complex processing conditions and mechanisms of deformation and relaxation over the course of spinning and drawing.4 The degree of polymer chain extension during the fiber drawing is related to topological constraints or entanglements, which restrict the motion of the polymers due to the slow reptation dynamics associated with them. The drawability of the polymer solutions and melts significantly depends on their entanglement density. It has been shown that the attained drawability of flexible polymers crystallized from the melt is limited by the presence of molecular entanglements.4,5 While polymer chains undergo a transition from a random coil to a fully extended chain in dilute solutions,6 only sections of a chain experience the coil−stretch transition in an entangled melt during draw.1 There are challenges in the use of molecular dynamics (MD) simulations associated with the limitation in the time scale, where simulation times are typically not sufficient to observe primary nucleation, and flexible PE chains retain a meltlike configuration without imposed deformations.7 This makes it © XXXX American Chemical Society

difficult to compare directly computational results with experimental data. Nevertheless, MD simulations can provide meaningful insight in molecular mechanisms of the polymer chain orientation and crystallization from entangled melts during FIC to complement experimental efforts by elucidating mechanisms of orientation, elongation, and disentanglement of chains in isolation of other competing molecular and larger scale factors. Recent coarse-grained MD simulations on the crystallization of long polymers and blends in concentrated solutions and melts8−12 showed that entanglements control the thickness of crystalline lamellar structures. Although highly coarse-grained MD simulations could provide faster dynamics,13 modeling of the detailed structure−property relationships requires the retention of atomistic details to provide an accurate description of structure at meltlike densities.7 Theoretical pictures of the PE crystallization have been corroborated and extended successfully with molecular simulations using the united atom (UA) representation of the PE, where the carbons and bonded hydrogens are clustered together into a single coarse-grained particle.7,14−17 For instance, the free energy barrier of nucleation and its dependence on temperature were obtained, and a quantitative model of homogeneous nucleation of a polymer melt from first-principles was proposed.15,16 The crystal nucleation and Received: July 17, 2018 Revised: October 9, 2018

A

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

with an equilibrium bond angle θ0 of 110° and an angle bending constant ka of 60.0 kcal/(mol rad2). The dihedral angle rotation potential is expressed in the following form:

chain extension in crystalline domains happening during the uniaxial deformation have been also investigated previously with relatively small model systems.7,14 Ko et al.14 characterized the crystallization of PE by performing MD simulations at multiple temperatures below the melting temperature at zero load, starting from an oriented PE melt represented by 20 chains, each consisting of 400 UA sites, prepared by uniaxial deformation at 425 K. They observed the temperature dependence of lamellar thickness and tilted chain lamellae in accord with experimental data. Despite the significant impact that these computational studies have made to our understanding of the crystallization in flexible polymers, previous works considered relatively small model systems, which may influence the morphological properties such as number and sizes of crystallites through the limit imposed by periodic boundaries of the simulation cell. In this work, we study ordering and crystallization of well-entangled PE melts during uniaxial tension and subsequent cooling with a significantly larger model system represented by thousands of PE chains, each consisting of 1000 UA sites. The larger model system enabled us to perform a more detailed investigations of structural and morphological properties of PE such as entanglement, chain orientation, void formation, chain-end distributions, and nucleation of multiple crystallites during elongation and crystallization.

3

Ud(ϕi) =

n=1

∑ UvdW + ∑ Ub + ∑ Ua + ∑ Ud

(1)

where UvdW, Ub, Ua and Ud are nonbonded, bond stretching, bond angle bending, and dihedral angle rotation energy terms, respectively. The interactions between i and j UA sites of different chains and between those separated by four or more bonds are expressed by a 12−6 Lennard-Jones potential UvdW(rij) = 4εij[(σ /rij)12 − (σ /rij)6 ]

(2)

where rij represents the distance between the sites i and j. The well depth εij between CH2−CH2 and CH3−CH3 UA sites is 0.09344 and 0.22644 kcal/mol, respectively, while the van der Waals parameter σ was set to 4.009 Å. A cutoff of rc = 2.5σ and a long-range van der Waals tail correction26 were also utilized. Bonded interactions are described by potentials that account for bond stretching, angle, and torsion. A harmonic bond stretching potential for bond i is expressed by Ub(li) = k b(li − l0)2

(3)

where li is the length of the bond i, the equilibrium bond length l0 is set to 1.53 Å, and the bond stretching constant kb equals 317 kcal/(mol Å2). The harmonic bond angle bending potential with respect to the bond angle θi for triplet i is described by Ua(θi) = ka(θi − θ0)2

(5)

where ϕi is the bond torsion angle. We use the following torsion bond constants: k1 = 1.6 kcal/mol, k2 = −0.867 kcal/ mol, and k3 = 3.24 kcal/mol. This force field utilized in our study was originally developed to describe the structural and mechanical properties of PE melts and has been shown to yield realistic values for persistence length and glass transition temperature14 and to reproduce a variety of melt phase structures and dynamics.16 However, this UA force field does not describe the orthorhombic or lower symmetry of the PE crystal lattice, favoring the pseudohexagonal crystal phase.16,19,22 Although experimental research has shown that the majority of the crystalline phase of PE is formed by orthorhombic crystals,4 it was also demonstrated that PE initially crystallizes from amorphous melts into a transient hexagonal phase under extended crystallization conditions and then transforms into a final orthorhombic phase.27,28 The glass transition of C1000 for this force field has been estimated to be 223 K,16 which is in agreement with the range 120−220 K obtained from experimental observations.29 Note that the simulation values tend to slightly overestimate the glass transition temperature due to very high cooling rates used in MD simulations.16 The crystallization temperatures for the force field are estimated to be 280 K for C15016 and 300 K for C192.24 They are both lower than the crystallization temperature of 349 K for C140 obtained from experimental data.30,31 The melting temperature of this force field is also close to the experimental value of 396.4 K.15,16,22,23,32 Therefore, temperatures ranging from 250 to 375 K are expected to be appropriate to study crystallizations and characterizations of semicrystalline PE under thermodynamic equilibrium between amorphous and crystalline phases under this force field as in previous simulation studies.14,16,22,23 We studied the crystallization of PE at 300 K as will be described later. 2.2. Preparation of Initial Configuration. To study large and well-entangled polymer systems, we prepared an initial configuration of PE melt using the fast equilibration procedure for million-atom systems of highly entangled linear PE chains as described in detail previously.25,33 The procedure is implemented with a massively parallel MD code LAMMPS34,35 with a soft-core harmonic potential. The amorphous PE system was built in a cubic simulation cell at T = 400 K, the approximate melting temperature, and had 1500000 UA sites in total, composed of 1500 PE chains, where each chain contains 1000 UA sites (Mw = 14 kg/mol). The resulting mass density at 400 K was 0.82 g/cm3, which is consistent with 0.85 and 0.87 g/cm3 at 350 and 300 K, respectively, estimated from previous MD simulations with the same force field.22,23,25 The box length of the cubic simulation cell of the initial PE melt was 348.77 Å. The average end-to-end distance ⟨R⟩ in a long chain can be estimated from the formula for the mean-square end-to-end distance ⟨Ree2⟩ = C∞Nbl2, where l is the average length of the backbone bond, C∞ is the Flory characteristic ratio, and Nb is the number of bonds of a chain.36 With l = 1.52

2. MODELS AND METHODS 2.1. Force Field. In this work, we employ a UA representation of the PE, where the interactions between the particles are described by the UA force field developed by Paul et al.18 and used in previous computational studies of PE.7,14−16,19−25 In this force field, the interactions are described by the total interaction potential Utotal given by Utotal =

∑ (kn/2)[1 + (−1)n+ 1 cos nϕi]

Å, C∞ = 7.79,25 and Nb = 999, ⟨R⟩ = R ee 2 is estimated to be

(4) B

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules 134 Å. The box length significantly larger than the estimated R enables us to study entanglement, chain orientation, and nucleation of multiple crystallites during a large elongation process with little self-interactions of chains. 2.3. Uniaxial Tension and Quenching. Because the nucleation of the crystallites is too slow to be observed at the time scale of MD, the direct simulation of the crystallization of a well-entangled PE melt using UA description is not possible even with modern computing power.14 At the same time, it was shown previously that crystallization is greatly accelerated by orientation of chains.14 In this study, we employed uniaxial tensions at high strain rates to facilitate the crystallization process via orientation of polymer chains. Three uniaxial tensile deformations along the z-direction were performed on the PE melt at constant true strain rates ė of 107, 108, and 109 s−1. We investigated the effect of the strain rates on structural properties such as chain alignments and entanglement loss by analyzing the configurations obtained at these three different strain rates. As will be described later, we prepared an oriented melt configuration for the subsequent quenching using the uniaxial tensile deformation at the highest strain rate of 109 s−1 because a highly oriented structure was formed more efficiently at 109 s−1 than at lower strain rates. The extension or stretch ratio λ is defined as Lz/L0z , the length of the simulation cell in the z-direction Lz relative to its initial value L0z . The true strain e is related to the extension ratio λ by the equation e = ln(λ). The MD time step is 1 fs. The Langevin thermostat with damping time 100 ps is used to maintain the temperature at T = 400 K, and a Nosé−Hoover barostat with damping time 1 ps is used to maintain atmospheric pressure along the transverse directions. The thermostat imposes an isothermal condition during tensile deformation, which precludes the softening due to adiabatic heating observable experimentally during deformation at high strain rates.37 Because our polymer systems undergo the tensile deformation, their melting point is expected to increase compared with the melting temperature measured without stress due to the chain alignment.38 For example, it was shown that PE would melt at 25 K higher than its melting temperature under elongation.39 Lavine et al.7 observed small nucleus formation at ∼380 K. As the temperature of the simulation is close to the melting temperature, Tm ≈ 400 K,16 some extent of crystal nucleation is expected over the course of the deformation. The crystallization temperature was estimated to be Tc ≈ 300 K for the shorter chains using the same force field at the rates investigated.16,24 We stopped the deformation simulation at λ = 13.5 for all considered systems since lateral dimensions of the simulation cell become too small, i.e., comparable with chain dimension, at large values of strain. As will be shown below, a highly oriented melt configuration was prepared at an intermediate extension ratio λ = 7.4 during the deformation at ė = 109 s−1, and we chose it as the starting configuration of the subsequent crystallization simulation to avoid small lateral dimensions of the simulation cell. After quenching to 300 K, we performed a simulation for 180 ns under the constant pressure condition to study the crystallization. All simulations were performed using the LAMMPS MD simulator.34,35 2.4. Entanglement Analysis. The geometric algorithm implemented in the Z1 code was used to reduce the polymer configurations to their primitive paths (PP) through a procedure described previously.40,41 To obtain the number of interior “kinks” (Z) in the three-dimensional primitive path of each chain, the algorithm uses geometric moves in the limit of

zero primitive chain thickness to reduce the contour lengths of the chains monotonically. In this algorithm, the chain ends are fixed, and the uncrossability of chains is obeyed. The average number of Z, ⟨Z⟩, is proportional to the number of entanglement residing on a single chain, and the code also resolves the spatial location of these “entanglement” points. 2.5. Order Parameter. The orientational order and alignment are described quantitatively using local and global uniaxial orientational order parameters, P 2 , i and P2,z,7,14,19,20,22,23 respectively, both involving the second Legendre polynomial. The local orientational order parameter P2,i42 at each UA site i except at chain ends was defined by the relation P2, i = (3⟨cos2 θij⟩ − 1)/2

(6)

where the average was taken over all UA sites j, excluding chain end sites, within the cutoff distance rc from site i, where θij is the angle formed by the chords that join the UA sites bonded to i and j, respectively. Similarly, to estimate the overall alignment of PE chains along the tensile direction, we calculated a global orientational order parameter defined as P2, z = (3⟨cos2 θi , z⟩ − 1)/2

(7)

averaged over all UA sites i, excluding chain end sites, where θi,z is the angle between the z-axis and the chord that join the UA sites bonded to site i. 2.6. Characterization of Voids. As will be discussed later, voids were observed to form during tensile deformation. We have probed an extent of formation of voids larger than 10 Å (nanometer-size voids) by calculating the fraction of the volume occupied by the voids. First, we partitioned our simulation box into cubic bins, where the linear bin size is equal to 10 Å, which is similar to the cutoff radius of van der Waals forces, rc = 2.5σ ≈ 10 Å. Then, we computed the number of the vacant bins, nv (not occupied by any UA sites). We evaluated the volume of the empty space Vf = nvVb, where Vb is the volume of the bin. The fraction of these nanometersize voids f v is defined as f v = Vf/V, where V is the total volume of the system. Likewise, we computed f low and f high, the fractions of low and high densities, respectively, by counting the number of bins with corresponding densities. Low and high densities were defined as ρlow < 0.650 g/cm3 and ρhigh > 0.925 g/cm3, respectively, based on the lower and upper bounds of the local density distribution of the initial melt, which will be shown below. We estimated λsim max, defined as the maximum stretch ratio λ before formation of large (nanometer-size) voids by analyzing the profile of f v obtained from our simulations. λsim max is analogous to λmax, which is defined as an elongation before the mechanical failure or a maximum draw ratio estimated from experiments. It has been shown that the increase of elongation rate leads to the decrease of λmax.43 To characterize voids smaller than nanometer, we calculated the pore radius distribution of the voids based on Euclidean distance maps, using the efficient method developed by Bhattacharya et al.44 Starting from the spatial coordinates of a configuration, and assuming that each particle is represented by a hard sphere of rs = 2.0 Å, there is a separate pore radius distribution for each chosen radius (rt) of a spherical test particle.44 For given rt a pore radius distribution can be estimated by first dividing the configuration space into a very fine, but finite, grid (the grid rt of the largest sphere, which is C

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

associated with the whole chain motion becoming hindered by surrounding chains.36 tR can be estimated from tR = te(N/Ne)2, where te and Ne are the characteristic entanglement time and the entanglement spacing, respectively, and N is the number of UA sites in the whole chain. te does not depend on the molecular weight of a chain.36 Ne = 44 from Z1 code using the modified “S-kink” estimator,40,47 The characteristic entanglement time te for the employed PE model was estimated to be 7.3 × 10−9 s in a previous MD simulation study.24 Subsequently, the estimated Rouse time is tR = 3.8 × 10−6 s. Then the longest relaxation time td required to the entangled chain to relax fully could be estimated from td = 3te(N/Ne)3.48 td = 2.6 × 10−4 s for our chain length of C1000. With strain rates used in our simulation, the Weissenberg number, the strain rate compared to the inverse Rouse time estimated by Wi = ėtR, is ≥ 38, and significant orientations and stretching of polymer chains are expected during tensile deformations. 3.2. Void Formation and Crystal Nucleation. During tensile deformation, we observe that the polymer undergoes void formation and crystal nucleation. First, the density decrease over the course of deformation for polymer melts at ė = 109 s−1 shown in Figure 2 indicates the existence of voids

centered at the grid point and can be inserted without overlap into the existing structure. If a test particle sphere at a grid point overlaps with the existing structure, the grid point is simply disregarded, as it does not belong to the void space, and the amount of such irrelevant grid points is used to calculate the fraction of volume occupied by the structure, φs, where the void volume fraction equal to φp = 1 − φs. Afterward, one searches for every grid point P within the void space separately the largest radius Rp stored on the whole grid that is within distance Rp from P. This way every relevant grid point receives a value for a new radius, denoted as pore radius rp, which is equal or larger than its Rp by construction, and the distribution of pore radii follows from the set of rp’s. The resulting distribution reflects the radii of largest possible spheres that can be inserted into the configuration without overlap.

3. RESULTS AND DISCUSSION 3.1. Stress−Strain Response during the Tensile Deformation and Relaxation Time Scale. Stress−strain curves for PE melts for all strain rates as a function of the stretch ratio λ are shown in Figure 1. The employed strain

Figure 1. True tensile stress σt vs stretch ratio λ during uniaxial tension at the true strain rates ė = 107−109 s−1 and the temperature Τ = 400 Κ.

Figure 2. Mass density ρ vs stretch ratio λ during uniaxial tension at the true strain rates ė = 107−109 s−1 and the temperature 400 K.

rates ė = 107−109 s−1 are much higher than typical experimental strain rates. Stress−strain relationships depended strongly on strain rates. The overstretching of the bond, which we will describe later, may have caused large stresses developed at higher strain rates of 109 and 108 s−1. The shear modulus G can be estimated from linear rubber elasticity, i.e., the true tensile stress σt during the uniaxial deformation with the relation σt = G(λ2 − 1/λ).45 The estimated values of the shear moduli are 30.04, 8.39, and 2.05 MPa for PE melts at ė = 109, 108, and 107 s−1, respectively. Given that the estimated experimental value of the plateau modulus of PE melt at T = 413 K is 2.6 MPa,46 the estimated shear moduli at higher strain rates deviate significantly from the experimental value. However, our result at the lowest strain rate is in reasonable agreement with the experiment even though a careful analysis with additional simulations at lower strain rates and with different chain lengths may be needed for a more quantitative comparison. Here, we relate our reported results to the underlying relaxation time scales of the systems. One important relaxation time scale is the Rouse time of an entangled strand, tR,

or crazes, which were also observed in previous studies of cross-linked or glassy polymers.49−51 Crazes refer to microcracks bridged by polymer chains, where many microvoids are found within the craze structure. Our results show that the density increases over the course of the tensile deformation at the lowest strain rate ė = 107 s−1 for the entire strain range, while the polymer density initially rises and then drops during the imposed tension at ė = 108 s−1. During the deformation carried out at the fastest rate (ė = 109 s−1), the density decreases and reaches some plateau value at a large strain near λ = 9. Figure 3 shows f v, the fraction of volume occupied by the nanometer-size voids, along with f low and f high, fractions of low and high density regions, respectively. Figure 3a indicates that f v is small and less than 2%. As mentioned earlier, low and high densities were defined as ρlow < 0.650 g/cm3 and ρhigh > 0.925 g/cm3, respectively. The average equilibrium melt density at T = 400 K is 0.82 g/cm3 (Figure 2), and the density of crystalline PE was estimated to be around 1.0 g/cm3.16,22 Therefore, the region with the low density ρlow is associated with formation of the voids smaller than the nanometer (subnanometer-size D

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

7 −1 than λsim becomes even smaller than max = 5.70. f low at ė = 10 s f low of the initial melt over the course of the deformation (cf. Figure 3b), which indicates the absence of voids. In contrast, the denser fraction of the polymer designated by 7 −1 f high increases rapidly until λ reaches λsim (Figure max at ė > 10 s 3c). This rise slows down or approaches some plateau value in f high at ė of 108 and 109 s−1. The fraction f high increases over the entire course of deformation at ė = 107 s−1. While the polymer systems obtained from tensile deformations at rates of 108 and 107 s−1 have a similar amount of the high density fraction (f high ≈ 0.35) at λ = 13.5, the polymer at ė = 108 s−1 has a significant amount of the low density fraction (see Figure 3b). f high values for ė < 109 s−1 are significantly greater than the corresponding value obtained over the course of deformation at 109 s−1, which is around 0.18 (λ = 13.5). Distributions of the local density ρi computed within 10 Å of each UA site i at λ = 13.5 are depicted in Figure 4. The

Figure 4. Probability density distributions of the local density ρi at λ = 13.5 and Τ = 400 Κ for the true strain rates ė = 107−109 s−1. The reference local density distribution of the undeformed melt at λ = 1 is also given. Figure 3. Fractions of (a) nanometer-size voids, f v, (b) low-density regions, f low, and (c) high-density regions, f high, vs stretch ratio λ during uniaxial tension at the true strain rates ė = 107−109 s−1 and T = 400 K.

undisturbed melt (λ = 1) has a narrow density distribution with a mean value close to 0.82 g/cm3, the melt density at T = 400 K for the employed PE model as shown in Figure 4. The densities corresponding to the peaks in density distributions at ė of 107 and 108 s−1 are shifted higher. Both density distributions at higher strain rates, ė > 107 s−1, are broader than the distribution of the undisturbed melt due to formation of both low and high density fractions. We calculated the pore radius distribution at λ = 13.5 and for the undeformed melt at λ = 1 as shown in Figure 5 by using

voids). Similarly, the region with the high density ρhigh is related to the precursor of the crystalline nucleation formed by means of the chain alignment over the course of the simulation. Consequently, the increase in the high and low density fractions indicates the beginning of crystal nucleation for all considered strain rates and void formation for the strain rate ė = 109 s−1, respectively. Figure 3a shows that the nanometer-size void formation starts at λsim max = 4.05 for the deformation at the highest strain rate (ė = 109 s−1). At the lower strain rate of ė = 108 s−1, λsim max is increased to 5.70. Corresponding stress values at the onset of the nanometer-size void formation estimated from the results shown in Figure 1 were 673 and 432 MPa at ė = 109 and 108 s−1, respectively. These stress values are significantly larger than critical stress values estimated from simulations of crosslinked or semicrystalline PE at lower strain rates.22,23,52 We did not observe any void formation for the entire course of the tensile deformation up to λ = 13.5 at the lowest strain rate (ė = 107 s−1). Likewise, the low density fraction of polymer f low starts growing immediately with tension at ė = 109 s−1 and rises faster when the deformation reaches λsim max = 4.05 until it reaches some constant value at λ around 9. At the lower rate of ė = 108 s−1, f low does not change at small strains and increases at λ greater

Figure 5. Probability density distributions of the pore radius rp at λ = 13.5 and Τ = 400 Κ for the true strain rates ė = 107−109 s−1. The distribution of the undeformed melt at λ = 1 is also given. E

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules the method developed by Bhattacharya et al.44 as described earlier in detail. The peak position in the distribution for the undeformed melt is found at a slightly larger pore radius, reflecting larger carbon−carbon radial distances in the melt compared to the corresponding distances in the aligned chain segments obtained during deformations. We defined a void as a pore of the radius rp larger than the half of the van der Waals diameter of the UA site, σ/2 ≈ 2 Å for the employed force field. The distributions shown in Figure 5 confirmed that the samples obtained at ė = 107 s−1 and from the undeformed melt do not contain a significant number of pores to be considered as voids. The PE structures obtained at rates ė > 107 s−1 have broader pore size distributions, where rp could be larger than 2 Å. This verifies the occurrence of voids over the course of deformation at rates ė > 107 s−1. The distributions of the local density in a slice of the PE system with x, y, and z dimensions of 10, 100, and 300 Å, respectively, taken from the simulations at different strain rates at λ = 13.5 are reported in Figure 6. From this picture, one can

Figure 7. Number of entanglements per chain Z as a function of the extension ratio λ at the true strain rates ė = 107−109 s−1 and the temperature of 400 K. The result for the polymer melt at λ = 1 prior to the stretching is also shown.

has three distinct regimes: negligible entanglement loss, fast disentanglement, and slow disentanglement. Initially, within the first regime, the entanglement loss is insignificant for small strains (λ < 1.5), where chains seem to “slide” while most entanglements are preserved. A significant disentanglement follows for the larger strains of 1.5 < λ < 6.5. At strain rates exceeding 107 s−1, the pattern of disentanglement is similar. However, the entanglement loss starts earlier than the loss at the uniaxial tension at ė = 107 s−1 because chains may have less time to slide at faster deformation rates, which are much greater than the inverse Rouse relaxation time. The entanglement loss leads to greater orientation and/or crystal nucleation of the polymer material. At larger λ, the loss of entanglements slows down for all considered strain rates. 3.4. Chain Unfolding and Alignment. Because our results show that the polymer chains experience significant entanglement loss during uniaxial elongations, polymer chains are expected to undergo unfolding. To describe the extent of chain unfolding, we computed ⟨R⟩/Rmax, which is a ratio of the average end-to-end distance ⟨R⟩ to the value of largest possible length of chain Rmax.36 We report our results in Figure 8. All chains undergo stretching under considered strain rate range (Figure 8a). However, polymer systems experience void formation at strain rates ė > 107 s−1, namely, at λsim max = 4.05 and 5.70 at the strain rates ė = 109 and 108 s−1, respectively (Figure 3a). Our results show that the system could be extended much more for the lower rate without formation of voids. For example, in the case of ė = 107 s−1, no void formation was observed, and the final extension of polymer chains at λ = 13.5 is around half of their maximum extension, i.e., ⟨R⟩/Rmax ≈ 0.52. Over the course of the tensile deformation at the high strain rate, the chain segments containing entanglements can be stretched to different degrees depending on the relative locations along the chain. Here, we computed Re/Re,max, the ratio of the end-to-end distance to the largest possible extension of the entangled chain segment, along a polymer chain by dividing a chain into 10 strands, each strand containing 100 UA sites. Strand 1 and strand 10 represent the chain terminal sections. The value of 100 UA sites is about 2.3 times of the length of chain segment between entanglements Ne, which was estimated to be 44 earlier. Our results demonstrate that the polymer chains are more stretched in the middle compared to inherently mobile terminal sections (Figure 8c). While the extension of the entire polymer chain

Figure 6. Snapshots taken from the tensile deformation simulations at T = 400 K and at the true strain rates equal to (a) ė = 109 s−1, (b) ė = 108 s−1, and (c) ė = 107 s−1. The color scheme represents the specimen local density. The snapshots were taken at λ = 13.5 and represent a small slice of the actual simulation box, with x, y, and z dimensions of 10, 100, and 300 Å, viewed along the x-direction.

see the precursor of crystal nucleation marked by high local densities for all strain rates and void formations for the higher strain rates ė > 107 s−1. The polymer undergoes a crazing-like deformation at ė > 107 s−1, which appears as an aligned fibril network with voids.51 It can be seen as an increase of the regions with low density as shown in Figure 3b. 3.3. Entanglement Loss. Because entanglements play a significant role in long polymer chain crystallization under tension, we computed the number of entanglements per chain, ⟨Z⟩, for our polymer system using the Z1 code.40 As shown in Figure 7, ⟨Z⟩ decreased significantly during the tensile deformation at all strain rates. The results also show that the degree of polymer chain disentanglement depends on the time scale of the deformation process with respect to the relaxation time scales as described below. First, we discuss the entanglement loss at the lowest deformation rate ė = 107 s−1. Even though the strain rate of ė = 107 s−1 is 38 times larger than the inverse Rouse time of the polymer chain as described earlier, it is the slowest among three strain rates used in our study to probe the loss of entanglements during the deformation. The change of entanglement number over the course of the deformation F

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

contributed to high stress values shown in Figure 1. Apparently, deformation at low strain rates would not produce any significant stretching, and below some critical strain rate ėc the polymer chains will remain in the coiled state.6,53 As discussed earlier, the Weissenberg number given by Wi = ėtR is ≥38. In this range of Wi, the polymer chains are expected to undergo stretching. Our simulation findings are in agreement with experimental data, where a section of polymer chains experiences the coil−stretch transition in an entangled melt during draw.1 The development of orientational order and alignment can be described quantitatively using the local and global orientational order parameters, P2,i and P2,z defined in eqs 6 and 7, respectively.7,14,19,20,22 We report the average order parameter of the whole sample ⟨P2,i⟩ and P2,z in Figure 9. The

Figure 8. Ratio of the average end-to-end distance to the largest possible extension distance ⟨R⟩/Rmax during the tensile deformation at strain rates ė = 107−109 s−1 and the temperature of 400 K. (a) shows the results for the entire chain while (b) shows the results for an entangled strand in the middle of a polymer chain. Re/Re,max is the ratio of the average end-to-end distance to the largest possible extension of the entangled chain segment. (c) shows Re/Re,max along the chain at λ = 13.5. The strands 1 and 10 represent the chain terminal sections.

Figure 9. (a) Average local order parameter ⟨P2,i⟩ and (b) global order parameter P2,z as a function of stretch ratio λ.

does not significantly depend on the strain rate at the range of considered rates, the middle strands of the chains are more stretched at the higher strain rates (Figure 8a,b). For instance, the middle strands of the chains are reaching their maximum possible extension Re/Re,max = 0.97 at λ = 13.5 for ė = 109 s−1, which indicates that the middle section of the chain is pulled taut at this strain. We calculated the average bond length l at the end of deformation at different strain rates. At λ = 13.5, ⟨l⟩ is estimated to be 1.56 (0.04), 1.53 (0.04), and 1.52 (0.04) Å at strain rates of 109, 108, and 107 s−1, respectively. The reference value for the melt before tension is measured as ⟨l⟩ = 1.52 (0.04) Å. Similarly, the estimated average bond angles are 115(5)°, 111(4)°, and 110(5)° at strain rates of 109, 108, and 107 s−1, respectively, and 110(5)° for the melt. The uncertainties shown in parentheses are the standard deviations. The average bond length and angle increased at higher strain rates while they remained the same at the strain rate of 107 s−1. The significant deviation from the equilibrium bond distance and angle at the highest strain rate of 109 s−1 may have

polymer does not show a linear increase in ⟨P2,i⟩ even for a low strain, which suggests the significant conformational relaxation and shows that polymer behavior is nonglassy but rubbery.7 Our results demonstrate that the orientation of PE chains increases with the increase of the deformation rate. The increase of ⟨P2,i⟩ indicates the increase of relative orientation of the chains, where ⟨P2,i⟩ = 1 represents the ideal parallel mutual alignment of chains, and values close to zero correspond to an isotropic sample.54 Our results are in agreement with the earlier simulation studies on orientation and crystallization of PE during uniaxial extension, which found PE chains aligned toward the deformation.7 We observe a correlation between the degree of alignment and the density, which increases and decreases, respectively, with the increase of the strain rate. As we mentioned earlier, polymer undergoes void formation at rates higher than 107 s−1. Interestingly, the values of ⟨P2,i⟩ at λsim max for the systems under deformation rates of ė = 109 and 108 s−1 are similar and G

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules equal to 0.36 and 0.42, respectively. We computed ⟨P2,i⟩ = 0.51 for the polymer after deformation of λ = 13.5 at ė = 107 s−1. The results for P2,z are similar to ⟨P2,i⟩, showing that chains are aligned in direction of elongation (Figure 9b). Because the polymer could have crazing and void formation, a value of local order parameter alone is not used to designate the site as either crystalline or amorphous for oriented polymer melts. Instead, we defined the material to be crystalline if both conditions of density and order parameter for crystallinity are met, namely, ρi > 0.925 g/cm3 and P2,i > 0.6.14 Figure 10 shows a plot of crystallinity growth over the course of the tensile deformation.

Figure 10. Crystallinity ϕ as a function of the stretch ratio λ at the true strain rates ė = 107−109 s−1 and the temperature 400 K.

Figure 11. P(P2,i), the probability density distributions of the local order parameter P2,i, at different strain rates at (a) λ = 13.5 and (b) λ 7 −1 and λ = 13.5 for ė = 107 s−1. The distribution at = λsim max for ė > 10 s the reference melt at λ = 1 is also shown.

The extent of crystallization of PE at λ = 13.5 is small and measured as ϕ ≈ 0.25 at ė = 107 and 108 s−1. ϕ at ė = 108 s−1 is slightly greater than the value at 107 s−1 although, as we discussed earlier, the polymer material obtained over the course of deformation at ė = 108 s−1 has larger amount of voids. The highest strain rate of ė = 109 s−1 produces highly aligned PE structure with negligible amount of crystallinity due to crazing. The probability density distribution of the local orientation order parameters, P2,i, is shown in Figure 11. The mean value of P2,i is close to zero for the undisturbed melt (λ = 1), indicating a random distribution of chain segments (Figure 11). At λ = 13.5 at the end of the tensile deformation, the polymer segments are significantly aligned at ė > 107 s−1 while the polymer chain at ė = 107 s−1 exhibits less local ordering (Figure 11a). Although alignment of polymer chains is substantial for ė > 107 s−1, the polymer system does not possess a significant fraction of crystalline domains due to the low density and void formation (Figure 3a,b). Thus, we conclude that with increasing strain rate the fraction of the polymer material without voids decreases while the alignment of PE segments increases. Remarkably, the distribution of P2,i shows clear bimodal behavior when it was taken at the onset of the void formation in the polymer, at λ of 4.05 and 5.70, λsim max values for the strain rates ė = 109 and 108 s−1, respectively (Figure 11b). The second peak of this distribution at larger P2,i corresponds to increase of the chain alignment, while the first peak at smaller P2,i is located at the negative values of P2,i. The sign of the orientational parameters is determined by the direction of local chain segments relative to those of the surroundings. According to eq 6, P2,i close to −0.5 indicates that chain segments oriented nearly perpendicular to those of the neighboring chain segments.54 Thus, the negative values of the order parameters are expected if chains segments make sharp turn or fold relative to the orientations of the

neighboring chain segments. We found that this peak is related to the location of interior “kinks” as provided by the Z1 code. Figure 12 shows the distribution of the local order

Figure 12. Entanglement distribution in the oriented melt. The location of the entanglements are shown with black crosses. The color scheme represents the local order parameter P2,i. The snapshot was 8 −1 taken at λ = λsim max = 5.70 during the tensile deformation at ė = 10 s and T = 400 K and represents the fragment of 10 × 100 × 300 Å3 of the actual simulation box.

parameter and entanglement points in the representative snapshots from the simulation, where the entanglements computed by Z1 are located at the domains with the negative order parameter. These results suggest that long entangled chain segments are excluded from straight aligned segments. Similar segregations of long entangled chain segments were obtained in coarsegrained simulations of semicrystalline polymers.9,11 Because polymer chains are entangled in the melt state, individual chain-folding events may interfere with each other due to topological constraints. The first peak at smaller P2,i value is vanishing over the course of deformation due to the entanglement loss. 3.5. Crystallization upon Cooling below Melting Temperature. As described earlier, an oriented melt at λ = H

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules 7.4 formed via the uniaxial deformation at ė = 109 s−1 and T = 400 K was quenched to T = 300 K and subjected to the constant pressure condition for ∼180 ns to study crystallization. The density and the average order parameter ⟨P2,i⟩ of the oriented melt are 0.745 g/cm3 and 0.64 as shown in Figures 2 and 9, respectively. The crystallinity of the oriented melt sample at λ = 7.4, ė = 109 s−1, and T = 400 K, considering both density and order parameter, is equal to 0.013 (Figure 10). The oriented melt configuration composed of bundles of parallel chain segments ordered in the flow direction can be considered as a flow-induced precursor, an intermediate state between random coil and crystalline phase prepared during flow-induced crystallization.2,3 Because of quenching and rearrangement, the density increases from 0.745 to 0.927 g/ cm3 (Figure 13a) during the 10 ns equilibration period,

increase in the number of available conformations for the chain segments, the void closing, and the increase of the density. After the initial decrease of chain alignment, the polymer chains start ordering after 30 ns, and the order parameter slowly increases (Figure 13b). The crystallinity of the semicrystalline PE at T = 300 K is also computed from (1/ρa − 1/ρ)/(1/ρa − 1/ρc), where the densities of crystalline and amorphous PE at T = 300 K are equal to ρc = 1.00 g/cm3 and ρa = 0.87 g/cm3,16,25 respectively. That gives a crystallinity of 0.47. As described earlier, the crystallinity ϕ estimated by considering both density and order parameter was also 0.47. Because the equilibrated polymer at T = 300 K does not possess either voids or low density portions, the crystallinity could be also computed using the orientational order parameter P2,i alone.14 The crystallinity computed by designating all sites with P2,i > 0.6 as crystalline14 was 0.45. Consequently, all three estimates give the similar values of the crystallinity. We have also shown distributions of the local density and the order parameter of the oriented melt at T = 400 K (Figure 14a,c) and semicrystalline polymer at T = 300 K (Figure 14b,d). The polymer at T = 400 K has a characteristic bimodal distribution of P2,i when the large second peak relates to increase of the chain alignment and the smaller first peak is located at the negative values as shown in Figure 14c. As we established earlier, the first peak corresponds to the location of entanglements (Figure 15a). The polymer at T = 300 K represents two interpenetrated phases: crystalline and an entangled amorphous and manifest as a bimodal distribution of the local order parameter as shown in Figure 14d. Here, the first broad peak is distributed around zero value and represents the amorphous portion of the polymer, while the second narrow peak corresponds to the crystalline material (Figure 14d). The phase separation could be observed clearly in Figure 15c. The crystalline phase (P2,i > 0.6) is surrounded with pronounced interphase (0.2 < P2,i < 0.6). The amorphous phase is defined as P2,i < 0.2. Orientational ordering and conformational transition of PE during FIC were studied experimentally by analyzing gauche− trans transition from vibrational spectra using Fourier transform infrared (FTIR) and Raman spectroscopy.3,55−57 We analyzed the dihedral angles along the chains and estimated the fraction of trans conformations of PE configurations at several different stages. The fraction of trans conformations at the initial melt at 400 K was estimated to be 64%, which is consistent with the estimate from the previous simulation with the same force field.25 In the oriented melt configuration in the beginning of the crystallization, the fraction of trans conformations was 95%. In the semicrystalline PE at the end of the crystallization at 300 K, the fraction of trans conformations was 84%, which is smaller than that of the oriented melt but larger than the value at the initial melt, reflecting 47% crystallinity estimated earlier. These are consistent with the trend of the order parameter shown in Figures 9a and 13a and are in qualitative agreement with experimental observations.3,55−57 The loss of entanglements was negligible during the crystallization upon cooling as ⟨Z⟩ = 13.5 at T = 400 K to ⟨Z⟩ = 12.7 at T = 300 K. The entanglements are clustered at the disoriented portion of the polymer at T = 400 K. In particular, the entanglements are located in the chain segments that make a sharp turn relative to the local chain orientation

Figure 13. Progress of crystallization of PE after cooling to the temperature of 300 K. (a) and (b) show changes in the density ρ and the average local order parameter ⟨P2,i⟩, respectively, with respect to time t.

indicating the rapid closing of the voids. The average order parameter ⟨P2,i⟩ rapidly drops from 0.64 to 0.42, and afterward it rises slowly to 0.45 (Figure 13b). Although the average order parameter ⟨P2,i⟩ decreases, the crystallinity considering both density and order parameter increases from 0.013 to 0.47 upon cooling to 300 K. Our results suggest that entangled chain strands have become overextended in the course of the tensile deformation (Figure 8c). As a result, the entropy of a polymer melt is reduced. When the deformation ceases, polymer chains raise their entropy by the decrease of alignment, which leads to the I

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 14. Probability density distributions of the local density and the local order parameter P2,i of the oriented melt and semicrystalline PE. (a) and (c) show the distributions with the oriented melt at λ = 7.4 obtained over the course of deformation with strain rate ė = 109 s−1 at 400 K. (b) and (d) show the distributions with the semicrystalline PE configuration at the end of the subsequent cooling at T = 300 K.

compared to the one at the end of crystallization shown in Figure 15c. The snapshot in Figure 15c indicates multiple crystallites are formed during crystallization, in contrast to lamellar structures observed in experiments and simulations with smaller model systems.7,14 Figure 16 shows the distributions of the local density and the entanglement points during the crystallization process. The

Figure 15. Distribution of the local order parameter P2,i and entanglements in (a) the oriented melt at T = 400 K, (b) the configuration at 10 ns after quenching to T = 300 K, and (c) the semicrystalline PE at T = 300 K. The picture shows the 10 × 100 × 300 Å3 slice of the simulation box. The color scheme represents the local order parameter P2,i while the black crosses denote the entanglement locations. The distributions in (a), (b), and (c) represent typical distributions at three different time points during crystallization but do not cover the same region in the simulation cell. Figure 16. Distribution of the local density and entanglements in (a) the oriented melt at T = 400 K, (b) the configuration at 10 ns after quenching to T = 300 K, and (c) the semicrystalline PE at T = 300 K. The picture shows the same slice of PE shown in Figure 15. The color scheme represents the local density while the black crosses denote the entanglement locations.

(Figure 15a). Upon cooling, most entanglements are located in the amorphous, disordered phase, although some of them are located in the crystal/amorphous interface region (Figure 15b,c). The distribution at the early stage of crystallization shown in Figure 15b display smaller crystalline regions J

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

amorphous interface. The snapshots in Figure 18 illustrate that the majority of the chain ends are located near the crystal−

density distributions shown in Figure 16 are consistent with the density fluctuations observed experimentally during FIC prior to full crystallization.2,3 In the oriented melt, the entanglements are also located in the melt phase (0.65 < ρ < 0.925 g/cm3), and they are absent from the low (ρ < 0.65 g/ cm3) and high density phase (ρ > 0.925 g/cm3), as is seen in Figure 16a. In a similar fashion, entanglements are clustered in the melt density phase (Figure 16b,c) upon cooling to 300 K. Note that the fractional crystallinity of semicrystalline polymer at T = 300 K could be equally described by the orientational order parameter P2,i and the density (Figures 15c and 16c). We have also evaluated the location of chain ends of the polymer. Theoretical studies have demonstrated that the chain ends are unlikely to be found deeply in the crystal, since it is energetically costly and creates network defects.58 In addition, on the basis of sufficiently realistic chain models and simulations, it is concluded that the commonly used models, where the chain ends are found in the amorphous portion of semicrystalline polymers, unintentionally but inevitably contain layers with a higher density than in the crystallites.59 In the modified model, in addition to chain tilt in the crystallites, termination of chains at the crystal surface keeps dangling chain ends out of the crowded interfacial layer, reducing the density at the interface.59 Earlier solid-state NMR studies by VanderHart and Pérez60,61 assumed that the chain ends are distributed uniformly throughout the crystals and estimated that 50−75% of the chain ends reside in the crystalline region, which contradicts the older perception that the chain ends are found in the amorphous portion of semicrystalline PE.62,63 Recent experimental findings have demonstrated that the chain ends are indeed located at the crystal surface near the crystal− amorphous interface.59,64 To assign the location of ends of a polymer chain to phases, we evaluated ⟨P2,ends⟩, the average of the local order parameters around chain end CH3 beads. For this purpose, we computed the average value of the local order parameter of CH2 particles, which are located inside of a sphere with a radius of 20 Å, where the chain end CH3 unit represents a center of the sphere. ⟨P2,ends⟩ would be close to the unity if the ends are buried deeply in the crystal phase while it is expected to be close to zero if the ends are located in the amorphous phase (Figure 14). Figure 17 demonstrates that the mean value of ⟨P2,ends⟩ is ∼0.5, which implies that chain ends are located at the crystal−

Figure 18. Distribution of the chain ends in (a) the oriented melt at T = 400 K, (b) the configuration at 10 ns after quenching to T = 300 K, and (c) the semicrystalline PE at T = 300 K. The picture shows the same slice of semicrystalline PE shown in Figures 15 and 16. The black circles denote the locations of the chain ends. The color schemes represent the local order parameter P2,i.

amorphous interface even though some are found in the amorphous phase. This is in agreement with the experimental results,59,64 which showed that the chain ends are preferentially found near the crystal surface of semicrystalline PE. Next, we investigated the distribution of the order parameter along the polymer chain by estimating the average and the probability distribution of the order parameter from all the chains in the simulation cell as a function of the relative location along the chain as shown in Figure 19. The relative location along the chain was described by the location index ic, ranging from 1 for one end to 1000 for the other end. We find that the order parameters near the chain ends are lower than those in the middle part of the chain in the starting configuration taken from the oriented melt at a higher temperature T = 400 K as shown in Figure 19a, indicating the lower degree of alignment near the chain ends compared with the middle part of the chain. However, at the end of the simulation at T = 300 K, the order parameters near the chain ends became higher than those in the middle, indicating that monomers near the chain ends are preferably located in the crystalline phase of the semicrystalline PE while the chain ends are located at the interface in agreement with the recent experimental studies.64 At the beginning, the average order parameter increases monotonically as the relative location along the chain moves away from the chain end to the middle of the chain, as shown in Figure 19a. However, at the end of the run at T = 300 K, the highest average order parameter was obtained at ic of 14 near the chain ends, as shown in the inset of Figure 19b. Probability distributions of the order parameters near the chain end and in the middle of the chain shown in Figures 19c and 19d, respectively, are consistent with the

Figure 17. Probability density distribution of ⟨P2,ends⟩, the average of order parameters within 20 Å radius of the chain ends, in the semicrystalline PE at T = 300 K. K

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 19. (a) and (b) show the profiles of the average local order parameter ⟨P2,i⟩ along the chain in the oriented melt at 400 K and the semicrystalline PE at T = 300 K, respectively. The relative location in the chain was indicated by the location index ic ranging from 1 for one end to 1000 for the other end. The insets in (a) and (b) show the profiles near the chain end. (c) and (d) show the probability density distributions of local order parameter P2,i in the middle of a chain (ic = 500) and near the chain end (ic = 14) in the oriented melt at 400 K and the semicrystalline PE at T = 300 K, respectively.

profiles of the average order parameter shown in Figures 19a and 19b. The broader distribution around the peak corresponding to the crystalline region in the middle of the chain in Figure 19d indicates that locations of crystalline region are more variable in the middle of the chain than those near the chain ends. The progress of profiles of average order parameter ⟨P2,i⟩ near the chain end, in the middle of the chain, and over the whole chain at T = 300 K is shown in Figure 20a. The highest average order parameter was located near the chain end and increased over the time as shown in Figure 20a,b. The progress of the average order parameter in the middle of the chain follows a similar pattern as that of the average order parameter averaged over the entire chain as shown in Figure 20a. The range of the sites near the chain end with the highest average order parameters also increased over the time as shown Figure 20b, indicating the growth of the crystalline region.

4. CONCLUSIONS The molecular-level structures and subsequent mechanical properties of semicrystalline polymers are greatly affected by chain orientations and disentanglements during processing. In this study, we use MD simulations to study orientation, disentanglements, and crystallization of well-entangled PE melts during uniaxial tensile deformations at three strain rates of 107, 108, and 109 s−1 and after the subsequent quenching using an unparalleled large simulation system, which includes 1500 PE chains each containing 1000 UA sites. During tensile deformations at 400 K near the melting temperature, the PE melt experiences entanglement loss, chain alignment, and crystal nucleation, in agreement with experimental observations. Even though we observe a significant void formation at higher strain rates, we did not observe any void formation at the slowest strain rate of 107 s−1, which indicates that no void formation would be expected at

Figure 20. (a) Average local order parameter ⟨P2,i⟩ with respect to time t after cooling to 300 K for the monomer with the highest value of ⟨P2,i⟩ near the chain end, the monomer in the middle of the chain (ic = 500), and the entire chain. (b) Average local order parameter ⟨P2,i⟩ near the chain end at different times. ic for the monomer with the highest value of ⟨P2,i⟩ near the chain end varies from 8 to14.

L

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Highly Drawn Gel-Spun Ultrahigh Molecular Weight Polyethylene Fibers at the Final Stages of Drawing by Saxs, Waxs, and 1h SolidState Nmr. Macromolecules 2011, 44, 9254−9266. (5) Van Dingenen, J. L. J. 3 - Gel-Spun High-Performance Polyethylene Fibres A2. In High-Performance Fibres; Hearle, J. W. S., Ed.; Woodhead Publishing: 2001; pp 62−92. (6) De Gennes, P. G. Coil-Stretch Transition of Dilute Flexible Polymers under Ultrahigh Velocity Gradients. J. Chem. Phys. 1974, 60, 5030−5042. (7) Lavine, M. S.; Waheed, N.; Rutledge, G. C. Molecular Dynamics Simulation of Orientation and Crystallization of Polyethylene During Uniaxial Extension. Polymer 2003, 44, 1771−1779. (8) Luo, C.; Sommer, J.-U. Frozen Topology: Entanglements Control Nucleation and Crystallization in Polymers. Phys. Rev. Lett. 2014, 112, 195702. (9) Luo, C.; Kröger, M.; Sommer, J.-U. Entanglements and Crystallization of Concentrated Polymer Solutions: Molecular Dynamics Simulations. Macromolecules 2016, 49, 9017−9025. (10) Luo, C.; Kröger, M.; Sommer, J.-U. Molecular Dynamics Simulations of Polymer Crystallization under Confinement: Entanglement Effect. Polymer 2017, 109, 71−84. (11) Triandafilidi, V.; Rottler, J.; Hatzikiriakos, S. G. Molecular Dynamics Simulations of Monodisperse/Bidisperse Polymer Melt Crystallization. J. Polym. Sci., Part B: Polym. Phys. 2016, 54, 2318− 2326. (12) Higuchi, Y.; Kubo, M. Deformation and Fracture Processes of a Lamellar Structure in Polyethylene at the Molecular Level by a Coarse-Grained Molecular Dynamics Simulation. Macromolecules 2017, 50, 3690−3702. (13) Agrawal, V.; Peralta, P.; Li, Y.; Oswald, J. A PressureTransferable Coarse-Grained Potential for Modeling the Shock Hugoniot of Polyethylene. J. Chem. Phys. 2016, 145, 104903. (14) Ko, M. J.; Waheed, N.; Lavine, M. S.; Rutledge, G. C. Characterization of Polyethylene Crystallization from an Oriented Melt by Molecular Dynamics Simulation. J. Chem. Phys. 2004, 121, 2823−2832. (15) Yi, P.; Rutledge, G. C. Molecular Simulation of Crystal Nucleation in N-Octane Melts. J. Chem. Phys. 2009, 131, 134902. (16) Yi, P.; Locker, C. R.; Rutledge, G. C. Molecular Dynamics Simulation of Homogeneous Crystal Nucleation in Polyethylene. Macromolecules 2013, 46, 4723−4733. (17) Bourque, A.; Locker, C. R.; Rutledge, G. C. Molecular Dynamics Simulation of Surface Nucleation During Growth of an Alkane Crystal. Macromolecules 2016, 49, 3619−3629. (18) Paul, W.; Yoon, D. Y.; Smith, G. D. An Optimized United Atom Model for Simulations of Polymethylene Melts. J. Chem. Phys. 1995, 103, 1702−1709. (19) Kim, J. M.; Locker, R.; Rutledge, G. C. Plastic Deformation of Semicrystalline Polyethylene under Extension, Compression, and Shear Using Molecular Dynamics Simulation. Macromolecules 2014, 47, 2515−2528. (20) Lee, S.; Rutledge, G. C. Plastic Deformation of Semicrystalline Polyethylene by Molecular Simulation. Macromolecules 2011, 44, 3096−3108. (21) Lempesis, N.; in ‘t Veld, P. J.; Rutledge, G. C. Atomistic Simulation of the Structure and Mechanics of a Semicrystalline Polyether. Macromolecules 2016, 49, 5714−5726. (22) Yeh, I.-C.; Andzelm, J. W.; Rutledge, G. C. Mechanical and Structural Characterization of Semicrystalline Polyethylene under Tensile Deformation by Molecular Dynamics Simulations. Macromolecules 2015, 48, 4228−4239. (23) Yeh, I.-C.; Lenhart, J. L.; Rutledge, G. C.; Andzelm, J. W. Molecular Dynamics Simulation of the Effects of Layer Thickness and Chain Tilt on Tensile Deformation Mechanisms of Semicrystalline Polyethylene. Macromolecules 2017, 50, 1700−1712. (24) Ramos, J.; Vega, J. F.; Martínez-Salazar, J. Molecular Dynamics Simulations for the Description of Experimental Molecular Conformation, Melt Dynamics, and Phase Transitions in Polyethylene. Macromolecules 2015, 48, 5016−5027.

significantly lower strain rates used in typical experiments. We find that polymer chain segments make sharp turns relative to the local chain orientations at the entanglement junctions, which is reflected as bimodal distributions of the local order parameter. Even though oriented melt configurations observed during uniaxial tensile deformations display high degrees of chain alignment and orientational order and can be considered as flow-induced precursors, they lack truly amorphous and crystalline phases expected for a semicrystalline polymer after crystallization. By quenching an oriented melt configuration to 300 K, we obtain a semicrystalline polymer material, which has a percent crystallinity close to 50%. The loss of entanglements is found to be negligible upon cooling. The semicrystalline polymer has two interpenetrated phases: crystalline and an entangled amorphous. Most of entanglements are located in the amorphous phase of the semicrystalline material. Our simulation shows that the polymer chain ends are preferentially located at the crystal−amorphous interface, which is in agreement with the recent experimental results.



AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (I.-C.Y.). *E-mail: [email protected] (J.W.A.). ORCID

In-Chul Yeh: 0000-0001-9637-1700 Martin Kröger: 0000-0003-1402-6714 Jan W. Andzelm: 0000-0002-2451-3770 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors thank Professors Mark Robbins and Gregory Rutledge and Drs. Tanya Chantawansri, Thomas O’Connor, Robert Elder, Timothy Sirk, and K. Michael Salerno for useful discussions. Calculations were performed using resources at Air Force Research Laboratory and Navy DoD Supercomputing Resource Centers. The research reported in this document was performed in connection with contract/ instrument W911QX-16-D-0014 with the U.S. Army Research Laboratory. The views and conclusions contained in this document are those of Survice Engineering and the U.S. Army Research Laboratory. Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the use thereof. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation hereon. M.K. acknowledges support by the Swiss National Science Foundation through grant SNF 200021_156106.



REFERENCES

(1) Somani, R. H.; Yang, L.; Zhu, L.; Hsiao, B. S. Flow-Induced Shish-Kebab Precursor Structures in Entangled Polymer Melts. Polymer 2005, 46, 8587−8623. (2) Wang, Z.; Ma, Z.; Li, L. Flow-Induced Crystallization of Polymers: Molecular and Thermodynamic Considerations. Macromolecules 2016, 49, 1505−1517. (3) Cui, K.; Ma, Z.; Tian, N.; Su, F.; Liu, D.; Li, L. Multiscale and Multistep Ordering of Flow-Induced Nucleation of Polymers. Chem. Rev. 2018, 118, 1840−1886. (4) Litvinov, V. M.; Xu, J.; Melian, C.; Demco, D. E.; Möller, M.; Simmelink, J. Morphology, Chain Dynamics, and Domain Sizes in M

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

(48) Hou, J.-X.; Svaneborg, C.; Everaers, R.; Grest, G. S. Stress Relaxation in Entangled Polymer Melts. Phys. Rev. Lett. 2010, 105, 068301. (49) Mukherji, D.; Abrams, C. F. Mechanical Behavior of Highly Cross-Linked Polymer Networks and Its Links to Microscopic Structure. Phys. Rev. E 2009, 79, 061802. (50) Mukherji, D.; Abrams, C. F. Microvoid Formation and Strain Hardening in Highly Cross-Linked Polymer Networks. Phys. Rev. E 2008, 78, 050801. (51) Rottler, J.; Robbins, M. O. Growth, Microstructure, and Failure of Crazes in Glassy Polymers. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 011801. (52) Morozinis, A. K.; Tzoumanekas, C.; Anogiannakis, S. D.; Theodorou, D. N. Atomistic Simulations of Cavitation in a Model Polyethylene Network. Polym. Sci., Ser. C 2013, 55, 212−218. (53) Keller, A.; Kolnaar, H. W. H. Flow-Induced Orientation and Structure Formation. In Materials Science and Technology; Wiley-VCH Verlag GmbH & Co. KGaA: 2006. (54) Lafrance, C.-P.; Nabet, A.; Prud’homme, R. E.; Pézolet, M. On the Relationship between the Order Parameter and the Shape of Orientation Distributions. Can. J. Chem. 1995, 73, 1497−1505. (55) Pigeon, M.; Prud’homme, R. E.; Pezolet, M. Characterization of Molecular Orientation in Polyethylene by Raman Spectroscopy. Macromolecules 1991, 24, 5687−5694. (56) Chai, C. K.; Dixon, N. M.; Gerrard, D. L.; Reed, W. RheoRaman Studies of Polyethylene Melts. Polymer 1995, 36, 661−663. (57) Sasaki, S.; Tashiro, K.; Kobayashi, M.; Izumi, Y.; Kobayashi, K. Microscopically Viewed Structural Change of Pe During the Isothermal Crystallization from the Melt: Ii. Conformational Ordering and Lamellar Formation Mechanism Derived from the Coupled Interpretation of Time-Resolved Saxs and Ftir Data. Polymer 1999, 40, 7125−7135. (58) Predecki, P.; Statton, W. O. Dislocations Caused by Chain Ends in Crystalline Polymers. J. Appl. Phys. 1966, 37, 4053−4059. (59) Fritzsching, K. J.; Mao, K.; Schmidt-Rohr, K. Avoidance of Density Anomalies as a Structural Principle for Semicrystalline Polymers: The Importance of Chain Ends and Chain Tilt. Macromolecules 2017, 50, 1521−1540. (60) VanderHart, D. L.; Perez, E. A Carbon-13 Nmr Method for Determining the Partitioning of End Groups and Side Branches between the Crystalline and Noncrystalline Regions in Polyethylene. Macromolecules 1986, 19, 1902−1909. (61) Pérez, E.; Vanderhart, D. L. Morphological Partitioning of Chain Ends and Methyl Branches in Melt-Crystallized Polyethylene by 13c-Nmr. J. Polym. Sci., Part B: Polym. Phys. 1987, 25, 1637−1653. (62) Gedde, U. W. Polymer Physics; Chapman & Hall: London, 1995. (63) Keller, A.; Priest, D. J. Experiments on the Location of Chain Ends in Monolayer Single Crystals of Polyethylene. J. Macromol. Sci., Part B: Phys. 1968, 2, 479−495. (64) Wutz, C.; Tanner, M. J.; Brookhart, M.; Samulski, E. T. Where Are the Chain Ends in Semicrystalline Polyethylene? Macromolecules 2017, 50, 9066−9070.

(25) Sliozberg, Y. R.; Kröger, M.; Chantawansri, T. L. Fast Equilibration Protocol for Million Atom Systems of Highly Entangled Linear Polyethylene Chains. J. Chem. Phys. 2016, 144, 154901. (26) Sun, H. Compass: An Ab Initio Force-Field Optimized for Condensed-Phase Applications - Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102, 7338−7364. (27) Kakiage, M.; Uehara, H.; Yamanobe, T. Novel in Situ Nmr Measurement System for Evaluating Molecular Mobility During Drawing from Highly Entangled Polyethylene Melts. Macromol. Rapid Commun. 2008, 29, 1571−1576. (28) Kakiage, M.; Yamanobe, T.; Komoto, T.; Murakami, S.; Uehara, H. Transient Crystallization During Drawing from Ultra-High Molecular Weight Polyethylene Melts Having Different Entanglement Characteristics. Polymer 2006, 47, 8053−8060. (29) Gaur, U.; Wunderlich, B. The Glass Transition Temperature of Polyethylene. Macromolecules 1980, 13, 445−446. (30) Cormia, R. L.; Price, F. P.; Turnbull, D. Kinetics of Crystal Nucleation in Polyethylene. J. Chem. Phys. 1962, 37, 1333−1340. (31) Gornick, F.; Ross, G. S.; Frolen, L. J. Crystal Nucleation in Polyethylene: The Droplet Experiment. J. Polym. Sci., Part C: Polym. Symp. 1967, 18, 79−91. (32) Ungar, G.; Stejny, J.; Keller, A.; Bidd, I.; Whiting, M. C. The Crystallization of Ultralong Normal Paraffins: The Onset of Chain Folding. Science 1985, 229, 386−389. (33) Sliozberg, Y. R.; Andzelm, J. W. Fast Protocol for Equilibration of Entangled and Branched Polymer Chains. Chem. Phys. Lett. 2012, 523, 139−143. (34) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (35) Placzek, G.; Nijboer, B. R. A.; van Hove, L. Effect of Short Wavelength Interference on Neutron Scattering by Dense Systems of Heavy Nuclei. Phys. Rev. 1951, 82, 392−403. (36) Rubinstein, M.; Colby, R. H. Polymer Physics; Oxford University Press: Oxford, UK, 2003. (37) Furmanski, J.; Cady, C. M.; Brown, E. N. Time−Temperature Equivalence and Adiabatic Heating at Large Strains in High Density Polyethylene and Ultrahigh Molecular Weight Polyethylene. Polymer 2013, 54, 381−390. (38) Uehara, H.; Nakae, M.; Kanamoto, T.; Zachariades, A. E.; Porter, R. S. Melt Drawability of Ultrahigh Molecular Weight Polyethylene. Macromolecules 1999, 32, 2761−2769. (39) Derakhshandeh, M.; Hatzikiriakos, S. G. Flow-Induced Crystallization of High-Density Polyethylene: The Effects of Shear and Uniaxial Extension. Rheol. Acta 2012, 51, 315−327. (40) Kröger, M. Shortest Multiple Disconnected Path for the Analysis of Entanglements in Two- and Three-Dimensional Polymeric Systems. Comput. Phys. Commun. 2005, 168, 209−232. (41) Karayiannis, N. C.; Krö ger, M. Combined Molecular Algorithms for the Generation, Equilibration and Topological Analysis of Entangled Polymers: Methodology and Performance. Int. J. Mol. Sci. 2009, 10, 5054. (42) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Oxford University Press: New York, 1986. (43) Termonia, Y.; Allen, S. R.; Smith, P. Kinetic Model for Tensile Deformation of Polymers. 3. Effects of Deformation Rate and Temperature. Macromolecules 1988, 21, 3485−3489. (44) Bhattacharya, S.; Gubbins, K. E. Fast Method for Computing Pore Size Distributions of Model Materials. Langmuir 2006, 22, 7726−7731. (45) Treloar, L. R. G. The Physics of Rubber Elasticity; Oxford University Press: New York, 1975. (46) Fetters, L. J.; Lohse, D. J.; Colby, R. H. Chain Dimensions and Entanglement Spacings. In Physical Properties of Polymers Handbook; Mark, J. E., Ed.; Springer: New York, 2007; pp 447−454. (47) Hoy, R. S.; Foteinopoulou, K.; Kröger, M. Topological Analysis of Polymeric Melts: Chain-Length Effects and Fast-Converging Estimators for Entanglement Length. Phys. Rev. E 2009, 80, 031803. N

DOI: 10.1021/acs.macromol.8b01538 Macromolecules XXXX, XXX, XXX−XXX