Microstructural Dynamics of Polymer Melts during Stretching: Radial Size Distribution

Abstract

The transient elongational viscosity ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ of the polymer melt is known to exhibit strain hardening, which depends on the strain rate $\dot{\mathsf{\epsilon }}$. This phenomenon was elucidated by the difference of chain stretching in the entanglement network between extension and shear. However, to date, the microscopic evolution of polymer melt has not been fully statistically analyzed. In this work, the radial size distributions P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) of linear polymers are explored by dissipative particle dynamics during the stretching processes. In uniaxial extensional flow, it is observed that the mean radius of gyration ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$ and standard deviation $\mathsf{\sigma }\left(\mathrm{t}\right)$ remain unchanged until the onset of strain hardening, corresponding to linear viscoelasticity. Both ${\overline{\mathrm{R}}}_{\mathrm{g}}$ and $\mathsf{\sigma }$ rise rapidly in the non-linear regime, and bimodal size distribution can emerge. Moreover, the onset of strain hardening is found to be insensitive to the Hencky strain (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}\mathrm{t}$) and chain length (N).

1. Introduction

The mechanical responses to the extension of entangled polymer melts are crucial in rheology for guiding polymer processing [1]. For example, in the processes of film manufacturing, fiber spinning, and blow molding, extensional deformation is the dominant mode. Moreover, as the melt flows through a sudden contraction, it accelerates, and the extensional flow becomes predominant. Since those processes are often conducted under highly non-linear conditions, the extensional rheology cannot be derived from the shear rheological behavior. The extensional viscosity describing the resistance of a fluid to the extensional flow is frequently used to characterize the molten polymer, which is viscoelastic. The elongational viscosity ${(\mathsf{\eta }}_{\mathrm{e}})$ associated with the uniaxial extension is defined as the ratio of normal stress difference ($\mathsf{\tau }$) to the rate of strain ($\dot{\mathsf{\epsilon }}$), and the ratio between the elongation and shear viscosities, known as Trouton’s ratio, is equal to three for Newtonian fluids.

A lot of experimental set-ups have been developed for measuring extensional viscosity, such as Göttfert rheotens tester [2] and CaBER [3,4,5,6,7]. However, it is still a challenge to acquire ${\mathsf{\eta }}_{\mathrm{e}}$ accurately. The factors limiting the improved characterization of ${\mathsf{\eta }}_{\mathrm{e}}$ include flow instabilities and the rupture of the sample [8]. Moreover, spatial uniformity is difficult to achieve in producing purely uniaxial extension, and the steady-state elongation flow is seldom reached [9]. Measurements at a constant elongational rate ($\dot{\mathsf{\epsilon }}$) corresponding to stressing experiments are easier to perform. As a result, the stress has to increase with time $\mathsf{\tau }\left(\mathrm{t}\right)$ in order to maintain constant $\dot{\mathsf{\epsilon }}$ in uniaxial extensional flow, leading to the transient elongational viscosity ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$. It is known that the elongational viscosity ascends with the time of deformation and follows the linear viscoelastic response in the zero-strain rate limit, ${\mathsf{\eta }}_0\left(\mathrm{t}\right)$. However, at finite strain rates, an exciting feature of the elongational flow occurs, named strain hardening. During the stretching process, the non-linear response emerges, and the transient extensional viscosity starts to deviate upward from ${\mathsf{\eta }}_0\left(\mathrm{t}\right)$, recognized as the onset of strain hardening [1,10]. The behavior of ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ has been reported to be greatly affected by the strain rate, which is mainly manifested in strain hardening. As $\dot{\mathsf{\epsilon }}$ is increased in stressing experiments, the onset of strain hardening takes place in a shorter stretching time. That is, the non-linear viscoelasticity appears earlier.

It was thought that strain hardening was a feature indicative of specific polymer architecture and closely associated with long-chain branching [11]. However, even ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ of fairly monodisperse melts was found to display an upward deviation at sufficiently high strain rates [8]. In elongational flows, chain stretching is commonly believed to be responsible for strain hardening. Recently, the origin of strain hardening has been elucidated based on the comparison of the responses between uniaxial extension and simple shear [1]. Regardless of polymer architectures, well-entangled polymer melts tend to exhibit strain hardening because the chain disentanglement effect is partially suppressed in extension. Moreover, the geometric condensation of the load-bearing chains accompanied by the shrinking cross-sectional area strengthens the entanglement network effectively [1]. In general, the effect of the strain rate on the onset of strain hardening is simply manifested by the onset time. Nonetheless, some experimental studies from a few decades ago reported that the onset of the rapid stress growth took place at an approximately constant value of the total Hencky strain [12,13].

The dynamics of entangled polymer melts have been extensively studied by the reptation (or tube) model [14] and the molecular stress function (MSF) model [15,16]. In the former, the dynamical behavior of a chain due to its interactions with surrounding chains is accounted for by the imposed topological constraints. To estimate ${\mathsf{\eta }}_{\mathrm{e}}$ quantitatively, the MSF model is introduced by incorporating the interchain tube pressure effect into the reptation model, and it contains certain input parameters that are typically acquired from experiments [16,17]. Indeed, the molecular mechanism responsible for the elongational behavior of entangled polymers can be scrutinized via microscopic simulations, such as the primitive chain network (PCN) model [18,19] and molecular dynamics (MD) simulation [20,21]. In the PCN model, the rubber-like network is constructed by a sequence of subchains connecting consecutive entanglements (nodes) [18]. Network rearrangement in real space is facilitated by various chain dynamical mechanisms, including reptation, contour length fluctuation, and constraint release, based on the force balance on entanglements [19]. However, to agree with the experiments, appropriate molecular mechanisms such as stretch/orientation-dependent monomeric friction [18] have to be incorporated into the PCN model. In contrast, MD can provide direct information about chain dynamics in polymer processing [22,23] without prior knowledge of molecular mechanisms. Nonetheless, the high computational costs make MD difficult to simulate larger systems for longer time scales. As a result, the MD approach is seldom employed to study the elongational flow of polymer melt [20].

Coarse-grained molecular dynamics is able to offer microscopic information about stress, chain conformation, and dynamics [24] while requiring a more reasonable computational cost. Complications in experiments, including molecular weight distribution, can be eliminated simply in simulations. By representing a polymer as beads and springs, the steady-state elongational viscosity of weakly entangled polymer melt has been observed after strain hardening [25]. By mimicking the behavior of extensional flows in entangled polymer melts, observed patterns in extensional viscosity regarding time, rate, and molecular weight are replicated and chain conformations are related to extensional stress [26]. Recently, the viscosity overshoot of a blend of the ring and linear polymers has been reported for biaxial extensional flows [27]. In addition, in this paper, a new simulation methodology that uses a hierarchical triple scale approach is demonstrated to predict the dynamic and rheological characteristics of entangled polymer melts with a high molecular weight [28]. However, until now, the microstructural dynamics of polymer melts during stretching have not been analyzed directly. Particularly, polymer conformations associated with strain hardening have not been microscopically examined. In this work, the elongational flow of linear polymer melts is explored by coarse-grained dissipative particle dynamics (DPD). Both weakly and strongly entangled polymer melts are considered by varying the chain length. The transient elongational viscosity is acquired for various strain rates, and strain hardening is identified. The microstructural dynamics are investigated by monitoring the evolution of the radial size distribution of polymers. The relationship between strain hardening and the mean radius of gyration (and standard deviation) is then established. The influences of the strain rate and molecular weight on the onset of strain hardening are studied.

2. Methods

Similar to coarse-grained MD, DPD is a particle-based mesoscale simulation method [29] that can be used to investigate larger length and time scales than conventional molecular dynamics (MD) [29,30,31]. The DPD bead with mass m is comprised of a few molecules or atoms, and the time evolution of the beads is governed by Newton’s equation of motion [32,33]. The forces acting on the DPD bead are generally classified into three types: conservative force (${\mathrm{f}}_{\mathrm{ij}}^{\mathrm{C}}$), dissipative force (${\mathrm{f}}_{\mathrm{ij}}^{\mathrm{D}}$), and random force (${\mathrm{f}}_{\mathrm{ij}}^{\mathrm{R}}$). These forces are soft and repulsive, in addition to being short-ranged and pairwise-additive [33]. As the distance between any two beads (${\mathrm{r}}_{\mathrm{ij}}$) is larger than the cut-off radius (${\mathrm{r}}_{\mathrm{c}}$), these forces vanish. The total force exerting on bead i is then ${\mathrm{F}}_{\mathrm{i}}={\sum }_{\mathrm{j}(\ne \mathrm{i})}{(\mathrm{f}}_{\mathrm{ij}}^{\mathrm{C}}{+\mathrm{f}}_{\mathrm{ij}}^{\mathrm{D}}{+\mathrm{f}}_{\mathrm{ij}}^{\mathrm{R}})$. The conservative force decreases linearly with interparticle distance, ${\mathrm{f}}_{\mathrm{ij}}^{\mathrm{C}}={\mathrm{a}}_{\mathrm{ij}}(1-{\mathrm{r}}_{\mathrm{ij}}{/\mathrm{r}}_{\mathrm{c}}{)\hat{\mathrm{r}}}_{\mathrm{ij}}$, where the interaction parameter ${\mathrm{a}}_{\mathrm{ij}}$ denotes the interaction strength between two beads and ${\hat{\mathrm{r}}}_{\mathrm{ij}}$ is the unit vector of the inter-particle distance. Because the polymer melt is made of the same type of DPD bead, the interaction parameter between all pairs of beads is always set as ${\mathrm{a}}_{\mathrm{mm}}=25$ [34]. ${\mathrm{f}}_{\mathrm{ij}}^{\mathrm{D}}$ is proportional to the relative velocity between two DPD beads and ${\mathrm{f}}_{\mathrm{ij}}^{\mathrm{R}}$ is introduced to satisfy the fluctuation–dissipation theorem [35,36,37,38]. The conservation of momentum in DPD simulations is automatically fulfilled by all pairwise interacting forces between every two DPD particles [31,39]. As a result, the hydrodynamic behavior of the system can be observed more easily compared to MD with too many details of the molecular motion [31]. All the units in our simulations are scaled by the bead mass (m), cut-off distance (${\mathrm{r}}_{\mathrm{c}}$), and thermal energy (${\mathrm{k}}_{\mathrm{B}}\mathrm{T}$). Therefore, the time (t) is nondimensionalized by ${{(m\mathrm{r}}_{\mathrm{c}}^2{/\mathrm{k}}_{\mathrm{B}}\mathrm{T})}^{1/2}$ and stress (${\mathsf{\tau }}_{\mathrm{ij}}$) by ${\mathrm{k}}_{\mathrm{B}}{\mathrm{T}/\mathrm{r}}_{\mathrm{c}}^3$. In this work, the open software Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [40] was used for all simulations.

A linear polymer consists of repeating units linked only to two others [41]. It is made of a string of DPD beads, and each polymer chain possesses N beads. For simplicity, the molecular weight or chain length is represented by N. The string of DPD beads is connected with the harmonic spring (bond), ${\mathrm{F}}_{\mathrm{ij}}^{\mathrm{S}}={\mathrm{k}}_{\mathrm{s}}{(\mathrm{r}}_{\mathrm{ij}}-{\mathrm{r}}_{\mathrm{eq}}{)\hat{\mathrm{r}}}_{\mathrm{ij}}$, where the spring constant ${\mathrm{k}}_{\mathrm{s}}=100$ and equilibrium length ${\mathrm{r}}_{\mathrm{eq}}=0.4$. Note that the original DPD parameters proposed by Groot and Warren [33] are too soft and provide no specific restrictions for steric interactions. Thus, chain crossings may occur. However, our choice of the DPD parameters (${\mathrm{r}}_{\mathrm{eq}}{,\mathrm{k}}_{\mathrm{s}}{, \mathrm{and} \mathrm{a}}_{\mathrm{ij}}$) can effectively maintain the uncrossability of polymeric chains, and at the same time, the large integration step can still be employed. The verification of chain uncrossability is illustrated in the Supplementary Materials. Videos S1 and S2 also demonstrate that our choice of DPD parameters can indeed prevent chain-crossing events. The chain stiffness can be tuned by implementing additional bending forces between two consecutive bonds. The bending force, ${\mathrm{F}}^{\mathsf{\theta }}=-\nabla {\mathrm{U}}^{\mathsf{\theta }}$, is obtained by the bending potential, ${\mathrm{U}}^{\mathsf{\theta }}={\mathrm{k}}_{\mathsf{\theta }}{(\mathsf{\theta }-{\mathsf{\theta }}_{\mathrm{eq}})}^2$. Here, the bending constant is ${\mathrm{k}}_{\mathsf{\theta }}=2,$ and the equilibrium angle is ${\mathsf{\theta }}_{\mathrm{eq}}=\mathsf{\pi }$. In this work, the number density of the linear polymer system is set as $\mathsf{\rho }=3$. There are about ${6\times 10}^5$ DPD beads in the cuboid with a periodic boundary in all three spatial directions. Initially, the system is a cubic box with approximately $58.5\times 58.5\times 58.5$. The time step is chosen as $\Delta \mathrm{t}=0.01$ to avoid the divergence in simulations. Before conducting the elongation study, it is necessary for the simulation to reach equilibrium. At equilibrium, the properties of the system, such as the internal energy and radius of gyration, should remain constant over time.

In the uniaxial extension experiment, melts are stretched uniaxially at a constant strain rate ($\dot{\mathsf{\epsilon }}$) and then the extensional stress ($\mathsf{\tau }$) is measured. To maintain the constant deformation rate, the uniaxial stress has to be adjusted with time, leading to the time-varying elongational viscosity ${\mathsf{\eta }}_{\mathrm{e}}(\mathrm{t})=\mathsf{\tau }(\mathrm{t})/\dot{\mathsf{\epsilon }}$. The extensional stress is ${\mathsf{\tau }=\mathsf{\tau }}_{\mathrm{xx}}-{(\mathsf{\tau }}_{\mathrm{yy}}{+\mathsf{\tau }}_{\mathrm{zz}})/2$ [42], where the average virial stress ${\mathsf{\tau }}_{\mathsf{\alpha }\mathsf{\alpha }}$ is the negative value of the diagonal component ${\mathrm{P}}_{\alpha \alpha }$ of pressure tensor in the $\mathsf{\alpha }$ direction [43]. ${\mathrm{P}}_{\mathsf{\alpha }\mathsf{\alpha }}$ can be computed by ${\mathrm{P}}_{\mathsf{\alpha }\mathsf{\alpha }}=({\sum }_{\mathrm{i}}{\mathrm{m}}^{\mathrm{i}}{\mathrm{v}}_{\mathsf{\alpha }}^{\mathrm{i}}{\mathrm{v}}_{\mathsf{\alpha }}^{\mathrm{i}}/\mathrm{V})+0.5({\sum }_{(\mathrm{j}\ne \mathrm{i})}{\mathrm{F}}_{\mathsf{\alpha }}^{\mathrm{ij}}{\mathrm{r}}_{\mathsf{\alpha }}^{\mathrm{ij}}/\mathrm{V})$, where V is the volume of the simulation box. The first term relates to the components of the kinetic energy tensor, where ${\mathrm{m}}^{\mathrm{i}}$ and ${\mathrm{v}}_{\mathsf{\alpha }}^{\mathrm{i}}$ represent the mass and $\mathsf{\alpha }$-component of the velocity of bead i, respectively. The second term uses the components of the virial tensor, where ${\mathrm{F}}_{\mathsf{\alpha }}^{\mathrm{ij}}$ and ${\mathrm{r}}_{\mathsf{\alpha }}^{\mathrm{ij}}$ denote the $\mathsf{\alpha }$-component of the force and distance between beads i and j, respectively [44]. It is calculated for all pairwise interactions. In this work, the simulation cubic box ($l_0{}^3$) is stretched at a constant “true” strain rate (i.e., Hencky strain rate, ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$) in the x-direction; thus, the x-dimension box size grows non-linearly with time, ${l\left(\mathrm{t}\right)=l}_0{\exp (\dot{\mathsf{\epsilon }}}_{\mathrm{H}}\Delta \mathrm{t})$. Note that it is different from the elongational process based on the constant “engineering” strain rate (${\dot{\mathsf{\epsilon }}}_{\mathrm{e}}$). In order to keep the box volume constant, both of the box sizes in the y- and z-dimensions shrink equally over time.

The polymer melt has a lot of polymer conformations which can be simply described by the radial size distribution P(${\mathrm{R}}_{\mathrm{g}}$), where the radial size of a chain is represented by the radius of gyration (${\mathrm{R}}_{\mathrm{g}}$). The polymer size along the x-axis is depicted by ${\mathrm{R}}_{\mathrm{gxx}}=[{\sum }_{\mathrm{i}}{{(\mathrm{x}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{cm}})}^2{/\mathrm{N}]}^{1/2}$, where ${\mathrm{x}}_{\mathrm{cm}}$ is the center-of-mass position. The radius of gyration of the polymer is ${\mathrm{R}}_{\mathrm{g}}={(\mathrm{R}}_{\mathrm{gxx}}{}^2{+\mathrm{R}}_{\mathrm{gyy}}{}^2{+\mathrm{R}}_{\mathrm{gzz}}{}^2{)}^{1/2}$. The aspect ratio of the polymer is defined as ${\mathrm{s}=\mathrm{R}}_{\mathrm{gxx}}{/\left(\right(\mathrm{R}}_{\mathrm{gyy}}{+\mathrm{R}}_{\mathrm{gzz}})/2)$. Certainly, the mean value of $\mathrm{s}$ at equilibrium is unity. During the stretching process, both the ${\mathrm{R}}_{\mathrm{g}}$ and $\mathrm{s}$ of each chain can vary with time and the radial size distribution of the polymer melt P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) also changes as a function of time. The first moment of P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) is ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$, which combines the second moment to give the standard deviation $\mathsf{\sigma }$ of the distribution.

3. Results

The elongational viscosity (${\mathsf{\eta }}_{\mathrm{e}}$) is commonly determined by extensional flow experiments, such as filament-stretching via an extensional rheometer [45] and capillary break-up via an extensional rheometer [4,5,6,7]. In those experiments, melts are stretched uniaxially at a constant strain rate (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$), and then the extensional stress ${\mathsf{\tau }=\mathsf{\tau }}_{\mathrm{xx}}-{(\mathsf{\tau }}_{\mathrm{yy}}{+\mathsf{\tau }}_{\mathrm{zz}})/2$ is measured. To maintain the constant deformation rate, the uniaxial stress has to be adjusted with time, leading to the time-varying elongational viscosity ${\mathsf{\eta }}_{\mathrm{e}}{(\mathrm{t})=\mathsf{\tau }(\mathrm{t})/\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$. The transient elongational viscosity depends on the strain rate ${\mathsf{\eta }}_{\mathrm{e}}{(\mathrm{t},\dot{\mathsf{\epsilon }}}_{\mathrm{H}})$, and the characteristic stretching time corresponding to the onset of strain hardening decreases with increasing ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$ [1,46]. It is important to note that stretching and shearing are two distinct types of deformation processes, and the trends observed in our results for the stretching process may not necessarily apply to the shearing process.

3.1. Evolution of Radial Size Distribution of Polymers during Stretching

Figure 1a shows the simulation results of the time-varying elongational viscosities for linear polymers with $\mathrm{N}=130$ subject to various strain rates. For the purpose of comparison, the elongational viscosity is non-dimensionalized by the shear viscosity of monomers ($\mathsf{\mu }=0.85$) [47]. For sufficiently low strain rates (e.g., ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-4}{ \mathrm{and} 10}^{-3}$), their curves of ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ coincide with each other essentially. This consequence of the strain rate-independent behavior seems to imply that these two curves can represent the linear viscoelastic response ${\mathsf{\eta }}_0{\left(\mathrm{t}\right)=\mathsf{\eta }}_{\mathrm{e}}{(\mathrm{t},\dot{\mathsf{\epsilon }}}_{\mathrm{H}}\to 0)$. Evidently, other curves of ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ at different ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$ (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}\ge {5\times 10}^{-3}$) overlap with ${\mathsf{\eta }}_0\left(\mathrm{t}\right)$ initially but deviate at different stretching times, indicating the occurrence of strain hardening. The stress in the non-linear regime associated with strain hardening increases with time faster than that in the linear regime. Moreover, the onset point of strain hardening decreases with increasing strain rate. Obviously, the qualitative characteristics of the transient elongational viscosity acquired from DPD simulations agree well with those observed in elongational rheology experiments. The general explanation for the effect of the strain rate on strain hardening is that, as the strain rate is increased, the polymer chains will be stretched more rapidly and reach the non-linear viscoelastic regime earlier.

It is imperative to realize the strain rate effect in terms of the onset time (abscissa) of strain hardening. To further understand the influence of the strain rate on the elongational viscosity, the data of ${\mathsf{\eta }}_{\mathrm{e}}{(\mathrm{t},\dot{\mathsf{\epsilon }}}_{\mathrm{H}})$ can be transformed into ${\mathsf{\eta }}_{\mathrm{e}}{(\mathsf{\epsilon }}_{\mathrm{H}}{,\dot{\mathsf{\epsilon }}}_{\mathrm{H}})$ based on the Hencky strain ${\mathsf{\epsilon }}_{\mathrm{H}}={\ln (l/l}_0)$, where $l$ and $l_0$ denote the instantaneous and initial length of the sample [10]. That is, the transient elongational viscosity is plotted against the elongational deformation, as shown in the inset of Figure 1a. It can be observed that the elongational viscosity decreases with the strain rate subject to the same ${\mathsf{\epsilon }}_{\mathrm{H}}$. The simulation outcome is qualitatively consistent with the experimental result [48,49]. It is known that, at very small strain rates, the polymer deformation is allowed to relax back to the shape close to equilibrium. As a result, the resistance to further deformation is weak, and the stress $\mathsf{\tau }$ increases slowly. However, for larger strain rates, the relaxation process is relatively slow compared to the deformation rate [50]. Therefore, the nonequilibrium entanglement structure is strengthened to impede further deformation, leading to an increment of $\mathsf{\tau }$ to maintain the constant strain rate. Although $\mathsf{\tau }$ grows with ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$, their ratio ${\mathsf{\eta }}_{\mathrm{e}}{(\mathsf{\epsilon }}_{\mathrm{H}}{,\dot{\mathsf{\epsilon }}}_{\mathrm{H}}{)=\mathsf{\tau }/\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$ actually decreases with increasing ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$ at a fixed ${\mathsf{\epsilon }}_{\mathrm{H}}$, indicating that the increasing rate of $\mathsf{\tau }$ is less than that of ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$. However, as the polymer melt is stretched to a specified length, more work is still required for higher strain rates based on the plot of $\mathsf{\tau }$ versus ${\mathsf{\epsilon }}_{\mathrm{H}}$.

Obviously, the steady state elongational viscosity, which is independent of time ${\mathsf{\eta }}_{\mathrm{e}}{\left(\mathrm{t}\right)=\mathsf{\eta }}^{\ast }$ [48], is not attained for $\mathrm{N}=130$, as shown in Figure 1a. Presumably, a very large simulation box is required to sustain the long-time extension for reaching the steady state ${\mathsf{\eta }}_{\mathrm{e}}$ [25]. To demonstrate the possibility of obtaining ${\mathsf{\eta }}^{\ast }$, the polymer melt with short chains $\mathrm{N}=10$ is considered. Figure 1b shows the evolution of the elongational viscosity with time at different strain rates. The comparison between Figure 1a,b reveals significant differences between short and long chains. Although the initial values of ${\mathsf{\eta }}_{\mathrm{e}}$ of both $\mathrm{N}=10 \mathrm{and} 130$ have small differences, the latter becomes much greater than the former after stretching the polymer melts at the same ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$. This is not surprising because of the fewer entanglements in polymer melts and lower resistances to stretching for $\mathrm{N}=10$. Moreover, for sufficiently low strain rates (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={10}^{-2}{ \mathrm{and} 5\times 10}^{-3}$), ${\mathsf{\eta }}_{\mathrm{e}}$ of $\mathrm{N}=10$ grows with time initially but reaches a plateau eventually. The time-invariant characteristic of ${\mathsf{\eta }}_{\mathrm{e}}$ of $\mathrm{N}=10$ is also revealed by the ${\mathsf{\epsilon }}_{\mathrm{H}}$-independent feature illustrated in the inset of Figure 1b. Similar to $\mathrm{N}=130$, ${\mathsf{\eta }}_{\mathrm{e}}{(\mathsf{\epsilon }}_{\mathrm{H}}{,\dot{\mathsf{\epsilon }}}_{\mathrm{H}})$ decreases with increasing ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$ at a given ${\mathsf{\epsilon }}_{\mathrm{H}}$. Note that the steady state ${\mathsf{\eta }}_{\mathrm{e}}$ is also observed for $\mathrm{N}=20,$ but a longer time period is required.

Although the stress and viscosity can be obtained in experiments, the microstructure evolution of polymer melts is difficult to observe. In contrast, the time-varying microstructure, such as polymer conformations subjected to stretching, can be captured by coarse-grained molecular simulations. To understand the microscopic evolution of polymer conformations during extension, the changes in the distribution of the radius of gyration with time P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) are demonstrated in Figure 2. The average radius of gyration at equilibrium is ${\overline{\mathrm{R}}}_{\mathrm{g}0}$. Under stretching for $\mathrm{N}=130$, the distribution becomes widened and shifts toward larger values of ${\mathrm{R}}_{\mathrm{g}}$, as shown in Figure 2a at the strain rate ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-2}.$ Evidently, the increment of ${\mathrm{R}}_{\mathrm{g}}$ parallel to stretching (x-direction) is significantly more than the decrement of ${\mathrm{R}}_{\mathrm{g}}$ perpendicular to stretching (y- and z-direction). That is, the polymer conformation is elongated from a spherical shape. Two interesting findings are observed from P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$). It is found that the distribution P($\mathrm{t}=0.16$) before the onset of strain hardening is essentially the same as that at equilibrium. The deviation of P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) from the equilibrium distribution becomes significant after strain hardening. Moreover, a bimodal distribution emerges after stretching for a long time (e.g., $\mathrm{t}=38.4 \mathrm{and} 40$), corresponding to large Hencky strains. One peak is close to $1.5{\overline{\mathrm{R}}}_{\mathrm{g}0}$, while the other peak is located at about 3.5 ${\overline{\mathrm{R}}}_{\mathrm{g}0}$. In addition to the Hencky strain, the bimodal distribution also depends on the strain rate. As depicted in the inset of Figure 2a for ${\mathsf{\epsilon }}_{\mathrm{H}}=2$, the bimodal feature is more prominent for lower strain rates. The peak close to 2.5 ${\overline{\mathrm{R}}}_{\mathrm{g}0}$ is significantly higher than that close to 1.5 ${\overline{\mathrm{R}}}_{\mathrm{g}0}$. The appearance of the peak with a smaller radial size can be realized by the representative evolutions of polymer conformations shown in Figure S1 of the Supplementary Materials. Some polymers are reoriented to the stretching direction with small elongation, some polymers are extended but then arrested in a specific stretched state, and some polymers return to the weakly stretched state via local relaxation. Figure S2 of the Supplementary Materials depicts the typical evolutions of polymer conformations near the peak with a larger radial size. However, as the strain rate increases, the distribution of ${\mathrm{R}}_{\mathrm{g}}$ becomes more non-uniform because easily stretched polymers can be elongated further, and the number of polymers in a weakly stretched state decreases slightly.

In contrast to the unsteady stretching of long chains ($\mathrm{N}=130$), steady stretching is reached for short chains ($\mathrm{N}=10$). Figure 2b shows the time evolution of the distribution of ${\mathrm{R}}_{\mathrm{g}}$ for $\mathrm{N}=10$ at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-2}$. By comparing Figure 2b with Figure 2a, it is found that the polymer size distributions of $\mathrm{N}=10 \mathrm{and} 130$ are distinctly different during stretching. While the bimodal distribution appears for long chains, the size distribution of the short chains is weakly perturbed by the equilibrium distribution. Moreover, P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) eventually becomes time-invariant for $\mathrm{N}=10,$ but it still evolves toward the bimodal shape for $\mathrm{N}=130$. The consequence of this is that subject to the same strain rate, the relaxation allows short chains to recover their shapes close to the equilibrium structure, but it fails to reduce the significant elongation of long chains. Furthermore, the effect of the strain rate on the size distribution of short chains at ${\mathsf{\epsilon }}_{\mathrm{H}}=2$ is demonstrated in the inset of Figure 2b. At relatively low strain rates (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={10}^{-2}{,5\times 10}^{-3}{, \mathrm{and} 5\times 10}^{-4}$), it is found that their polymer size distributions are essentially the same as the equilibrium distribution. This result indicates that the relaxation is fast enough compared to the stretching for $\mathrm{N}=10$. Their stretching dynamics always follow the linear viscoelastic response and ultimately reach the steady state. On the contrary, for high strain rates (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={10}^{-1}{ \mathrm{and} 5\times 10}^{-2}$), the deviation from the equilibrium distribution is observed, revealing weak strain hardening, which cannot be clearly identified from ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$.

3.2. Polymer Conformations Associated with Strain Hardening

The onset of strain hardening is generally identified as the point at which the elongational viscosity deviates from the linear viscoelastic response. However, ${\mathsf{\eta }}_0\left(\mathrm{t}\right)$ is difficult to acquire accurately because the strain rate needs to be as small as possible. At ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}\to 0$, the thermal noise causes significant fluctuations of stress and ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ [25], as illustrated in Figure 1. It is desirable to understand the change in the microstructure of polymers near the onset point of strain hardening. Instead of the radial size distribution, the mean radius of gyration of polymers ${\overline{\mathrm{R}}}_{\mathrm{g}}$ is easier to use for the monitoring of the microstructural evolution during stretching. Figure 3a shows the variation of ${\overline{\mathrm{R}}}_{\mathrm{g}}{/\overline{\mathrm{R}}}_{\mathrm{g}0}$ with time at various strain rates for $\mathrm{N}=130$. It is found that, regardless of ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$, all curves of ${\overline{\mathrm{R}}}_{\mathrm{g}}{/\overline{\mathrm{R}}}_{\mathrm{g}0}$ remain at unity for a while and then upturn rapidly at some points. The turning point occurs earlier when the strain rate is higher. The comparison of the onset points between Figure 1a and Figure 3a indicates that they coincided with each other for a given strain rate. This reveals that the upturn of the elongational viscosity is actually accompanied by the drastic change in the microstructure of the polymer melt. That is, strain hardening takes place as the average radial size of polymers starts to grow from the equilibrium state. Although the stretching process stops at ${\mathsf{\epsilon }}_{\mathrm{H}}=2$ for all strain rates, the radial size of polymers is larger for higher ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$. To maintain a higher ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$, higher stress is required, which causes larger polymer deformation.

According to ${\overline{\mathrm{R}}}_{\mathrm{g}}{/\overline{\mathrm{R}}}_{\mathrm{g}0}$, the microstructure of the polymer melt is essentially the same as that at the equilibrium state before strain hardening. This result can be further illustrated by the plot of the standard deviation of ${\overline{\mathrm{R}}}_{\mathrm{g}}{/\overline{\mathrm{R}}}_{\mathrm{g}0}$ against time, as depicted in the inset of Figure 3a for $\mathrm{N}=130$. Again, before the occurrence of strain hardening, the standard deviation is a constant and the same as that at equilibrium, $\mathsf{\sigma }=0.235$. Similar to ${\overline{\mathrm{R}}}_{\mathrm{g}}{/\overline{\mathrm{R}}}_{\mathrm{g}0}$, the onset of strain hardening can be recognized at the point of the uprising of $\mathsf{\sigma }$. The equilibrium state can be maintained for longer at lower strain rates. $\mathsf{\sigma }$ grows rapidly after the onset of strain hardening, signifying the incremental increase in the width of the radial size distribution. This occurs simply because the upper bound is elevated due to stretching. Since the higher stress for maintaining the higher ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$ yields more polymers with strong deformation, the distribution of ${\mathrm{R}}_{\mathrm{g}}{/\overline{\mathrm{R}}}_{\mathrm{g}0}$ is wider for the larger strain rate at ${\mathsf{\epsilon }}_{\mathrm{H}}=2$.

According to P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) and ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$, the linear viscoelastic response is observed as they are maintained at those of the equilibrium state. The onset of strain hardening corresponds to their deviation from the equilibrium state. Apparently, the initial elongational perturbation to the polymer conformation is not clearly displayed in ${\mathrm{R}}_{\mathrm{g}}$. To capture the response of the conformation to the stretching process, the aspect ratio ($\mathrm{s}$) of the polymer is calculated. It is defined as the ratio of ${\overline{\mathrm{R}}}_{\mathrm{gxx}}$ (parallel to the stretching direction) to ${(\overline{\mathrm{R}}}_{\mathrm{gyy}}+{\overline{\mathrm{R}}}_{\mathrm{gzz}})/2$ (perpendicular to the stretching direction) of the chain. Figure 3b shows the variation of the aspect ratio with time at different ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$ for $\mathrm{N}=130$. Obviously, the aspect ratio is unity ($\mathrm{s}=1$) at equilibrium, and it ascends rapidly at a certain point. At the end of the stretching (${\mathsf{\epsilon }}_{\mathrm{H}}=2$), ${\overline{\mathrm{R}}}_{\mathrm{gxx}}$ can be as large as approximately eight to ten times (${\overline{\mathrm{R}}}_{\mathrm{gyy}}+{\overline{\mathrm{R}}}_{\mathrm{gzz}})/2$. The representative polymer conformations at different times are shown in the inset of Figure 3b, at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-3}$. Note that the significant deviation of the aspect ratio from unity can be identified before the onset of strain hardening. For example, one has ${\overline{\mathrm{R}}}_{\mathrm{g}}/{\overline{\mathrm{R}}}_{\mathrm{g}0}=1$ but $\mathrm{s}=1.1$ at $\mathrm{t}=15.6$ before strain hardening at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-3}$. This result reveals that the initial elongation tends to align the deformed polymers at the equilibrium state (with an aspect ratio greater than unity) with the stretching direction. That is, the equilibrium distribution P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) is not altered actually but the mean orientation emerges because of the alignment induced by elongation.

It is imperative to realize the strain rate effect on strain hardening from the plot of ${\overline{\mathrm{R}}}_{\mathrm{g}}$ against the stretching time in terms of the characteristic onset time ${\mathrm{t}}^{\ast }{(\dot{\mathsf{\epsilon }}}_{\mathrm{H}})$. To understand the relationship between ${\mathrm{t}}^{\ast }$ and ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$, the plot of ${\overline{\mathrm{R}}}_{\mathrm{g}}$ (t, ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$) can be redrawn against a dimensionless time ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}\mathrm{t}$, which is the Hencky strain ${\mathsf{\epsilon }}_{\mathrm{H}}$. Figure 3c shows the variation of the mean radius of gyration with the elongation deformation, i.e., ${\overline{\mathrm{R}}}_{\mathrm{g}}/{\overline{\mathrm{R}}}_{\mathrm{g}0}$ vs. ${\mathsf{\epsilon }}_{\mathrm{H}}$ for N = 130. For comparison, the plot of the standard deviation against ${\mathsf{\epsilon }}_{\mathrm{H}}$ is depicted in the inset. It is somewhat surprising to find that the deviation of ${\overline{\mathrm{R}}}_{\mathrm{g}}/{\overline{\mathrm{R}}}_{\mathrm{g}0}$ from unity occurs essentially at the same Hencky strain, i.e., ${\mathsf{\epsilon }}_{\mathrm{H}}\approx 0.1$, regardless of the strain rate. That is, despite differences in ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$, the onset of strain hardening always takes place at an approximately constant ${\mathsf{\epsilon }}_{\mathrm{H}}$, which is consistent with the experiments [12,13]. In addition to the mean of the radial size distribution, similar behavior is observed for the variation of the standard deviation with ${\mathsf{\epsilon }}_{\mathrm{H}}.$ The upturn of $\mathsf{\sigma }$ from the equilibrium value ${\mathsf{\sigma }}_0$ takes place essentially at the same point (${\mathsf{\epsilon }}_{\mathrm{H}}\approx 0.1$) for all strain rates. This finding indicates that the linear viscoelasticity is always obeyed by very small Hencky strains, but the non-linear viscoelasticity emerges beyond the critical value (${\mathsf{\epsilon }}_{\mathrm{H}}>0.1$). Evidently, the radial size distribution P (${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) of the linear viscoelastic regime is independent of ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$, but that of the non-linear regime depends on ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$. Moreover, this finding also reveals that the characteristic onset time is inversely proportional to the strain rate, ${\mathrm{t}}^{\ast }{~\dot{\mathsf{\epsilon }}}_{\mathrm{H}}{}^{-1}$.

For long chains ($\mathrm{N}=130$), strain hardening is easy to observe, but the steady elongational viscosity is difficult to reach. On the contrary, for short chains ($\mathrm{N}=10$), the steady ${\mathsf{\eta }}_{\mathrm{e}}$ can be seen, but strain hardening cannot be recognized easily. Following the concept of Figure 3 for $\mathrm{N}=130$, Figure 4 shows the time evolution of polymer conformations for $\mathrm{N}=10$ at different strain rates. As demonstrated in Figure 4a, strain hardening can be clearly identified by ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$ at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={10}^{-1}{ \mathrm{and} 5\times 10}^{-2}$. However, at lower strain rates (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={10}^{-2}{ \mathrm{and} 5\times 10}^{-3}$), the curve of ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$ does upturn at some point (weak strain hardening) but arrive at a plateau eventually (steady ${\mathsf{\eta }}_{\mathrm{e}}$). At the lowest strain rate (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-4}$), the ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$ curve remains at unity all the time, indicating that the process follows the linear viscoelastic response and finally reaches the steady ${\mathsf{\eta }}_{\mathrm{e}}$. Note that the curve of ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-4}$ is not shown in Figure 1b because it coincided closely with ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-3}$.

Another dynamic feature of the radial size distribution P(${\mathrm{R}}_{\mathrm{g}}/{\overline{\mathrm{R}}}_{\mathrm{g}0},\mathrm{t}$) is the standard deviation $\mathsf{\sigma }\left(\mathrm{t}\right)$, which is illustrated in the inset of Figure 4a for $\mathrm{N}=10$. It is found that the behavior of $\mathsf{\sigma }\left(\mathrm{t}\right)$ of $\mathrm{N}=10$ is distinctly different from that of $\mathrm{N}=130$. After the onset of strain hardening for ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={10}^{-1}{ \mathrm{and} 5\times 10}^{-2}$, the standard deviation $\mathsf{\sigma }\left(\mathrm{t}\right)$ starts to grow fast because short chains are more easily stretched due to their weaker entanglement effects. Unlike $\mathsf{\sigma }\left(\mathrm{t}\right)$ of $\mathrm{N}=130$, the curves reach their peak and descend afterward, revealing that the relaxation of polymer conformations comes into play. Nonetheless, the steady state has not been reached, and hence the standard deviation still decays toward the equilibrium value. In contrast, at low strain rates (${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}\le {10}^{-2}$), $\mathsf{\sigma }\left(\mathrm{t}\right)$ always fluctuates around the equilibrium value of $0.109$ during the stretching process. This shows that the relaxation of polymer conformations dominates over the stretching of short chains at low strain rates.

The response of the polymer conformations along the stretching direction for $\mathrm{N}=10$ is also demonstrated in terms of the aspect ratio ($\mathrm{s}$) in Figure 4b. Again, a significant rise in the aspect ratio is accompanied by strain hardening at high strain rates. The value of $\mathrm{s}$ can be as high as about three. However, at lower strain rates, the rising aspect ratio reaches a plateau, which is consistent with the steady viscosity. The representative conformations at different times are depicted in the inset of Figure 4b at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-3}$. At the lowest strain rate, the value of $\mathrm{s}$ remains essentially at unity with a very slight deviation. This result discloses the fact that the equilibrium microstructure (radial size distribution) of short chains is always maintained at sufficiently low strain rates during the stretching process because of their rapid relaxation. Although the radial size distributions P (${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) are very similar at ${\mathsf{\epsilon }}_{\mathrm{H}}=2$ for the lowest three strain rates (see the inset of Figure 2b), their aspect ratios can still be distinguished because the orientation effect is not revealed in P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$). The incremental increase in the elongational viscosity is attributed to the gradual alignment of deformed chains associated with the equilibrium state, which is achieved during stretching.

3.3. Effects of Molecular Weight

In addition to the strain rate, it is well-known that molecular weight is the property that affects the resistance to unentangling under strain as well [51]. That is, the elongational viscosity is expected to depend on the chain length, ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{N}\right)$. In general, the melt strength of elongational viscosity is enhanced as the molecular weight increases [52]. Figure 5a shows the growth of ${\mathsf{\eta }}_{\mathrm{e}}$ with time at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-3}$ for various chain lengths, from $\mathrm{N}=30 \mathrm{to} 390$. The typical behavior of ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ is observed for each $\mathrm{N},$ and the steady state is not reached. However, two regimes can be distinguished: the short and long time periods. At a short stretching time, all the elongational viscosities behave similarly regardless of the molecular weight. Moreover, they are essentially collapsed together even at different strain rates, as illustrated in the inset of Figure 5a. In contrast, at a sufficiently long time, ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ seems to separate from the main curve (asymptotic behavior) represented by the very long chain length, such as $\mathrm{N}=390$. Obviously, the extent of the deviation from the main curve is more significant as the chain length is shorter. Because the linear viscoelastic response ${\mathsf{\eta }}_0\left(\mathrm{t}\right)$ varies with the molecular weight, it is difficult to identify the onset point of strain hardening for a given strain rate.

The characteristics of the two regimes displayed in ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ can be further understood by examining the radial size distribution. According to the identified onset point of strain hardening at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-3}$ for $\mathrm{N}=130$, P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) at $\mathrm{t}=20$ (close to the onset point) remains the same as that at its equilibrium state. Figure 5b shows the radial size distributions at $\mathrm{t}=20$ for $\mathrm{N}=30,60, \mathrm{and} 260.$ Evidently, their P(${\mathrm{R}}_{\mathrm{g}}$) at $\mathrm{t}=20$ still coincides with the corresponding equilibrium distributions. This outcome reveals that P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) is unchanged in the short-time regime, and the onset point of strain hardening is not sensitive to the molecular weight. On the contrary, the radial size distribution deviates from the equilibrium distribution in the long-time regime. As shown in the inset of Figure 5b, the bimodal distribution appears at ${\mathsf{\epsilon }}_{\mathrm{H}}=2$ only for $\mathrm{N}=60$, but not for $\mathrm{N}=30 \mathrm{and} 260$. Similar to $\mathrm{N}=130$ (see Figure 2a), the presence of the peak with a larger radial size implies that some polymers of $\mathrm{N}=60$ have already elongated to their strongly stretched state. On the contrary, as the chain length becomes significantly longer than N = 130, most of the polymers of N = 260 fail to be extended to the strongly stretched state; thus, the bimodal distribution is not obvious. Certainly, for short chains (N = 30), the relaxation is much faster than that of N = 130 and P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) shifts simply toward larger radial sizes upon stretching.

To understand the effect of molecular weight, the microscopic evolution of polymer conformations in terms of ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$ is illustrated in Figure 6a at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-3}$. In the initial period, all ${\overline{\mathrm{R}}}_{\mathrm{g}}{/\overline{\mathrm{R}}}_{\mathrm{g}0}$ remain in unity regardless of the molecular weight, indicating the regime of linear viscoelastic response. It is interesting to find that the upturn of all ${\overline{\mathrm{R}}}_{\mathrm{g}}/{\overline{\mathrm{R}}}_{\mathrm{g}0}$ curves seem to take place essentially at the same moment, revealing that the onset point of strain hardening is insensitive to molecular weight. Nonetheless, the growth rate of ${\overline{\mathrm{R}}}_{\mathrm{g}}$ varies with the chain length and ascends with increasing $\mathrm{N}$. When $\mathrm{N}$ is large enough, the curve of ${\overline{\mathrm{R}}}_{\mathrm{g}}/{\overline{\mathrm{R}}}_{\mathrm{g}0}$ reaches an asymptotic behavior, similar to the behavior of ${\mathsf{\eta }}_{\mathrm{e}}/{\mathsf{\eta }}_0$ (see Figure 5a). In addition to the mean (${\overline{\mathrm{R}}}_{\mathrm{g}}$), the characteristic of the radial size distribution can also be elucidated by the standard deviation $\mathsf{\sigma }$. The inset of Figure 6a illustrates the variation of ${\mathsf{\sigma }/\mathsf{\sigma }}_0$ over time for different molecular weights, where ${\mathsf{\sigma }}_0$ is the standard deviation of P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) at the equilibrium state. Again, the onset of strain hardening (upturn from unity) seems to occur at about the same time. A maximum is observed for short chains (N = 30 and 60) due to their fast relaxation. The asymptotic behavior of $\mathsf{\sigma }\left(\mathrm{t}\right)$ for long chains is observed as well.

The influence of the molecular weight on polymer conformations upon stretching can be illustrated by the aspect ratio (s) as well. The analysis of the strain rate effect indicates that the aspect ratio is more sensitive to the stretching process than the radius of gyration. In fact, the rapid growth of the aspect ratio occurs earlier than that of ${\mathsf{\eta }}_{\mathrm{e}}$ and ${\overline{\mathrm{R}}}_{\mathrm{g}}$, corresponding to the transition from the orientation to the stretching of polymers. Figure 6b shows the variation of s over time at ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}={5\times 10}^{-3}$ for different molecular weights. For comparison, the plot of the radius of gyration along the stretching direction (${\overline{\mathrm{R}}}_{\mathrm{gxx}}/{\overline{\mathrm{R}}}_{\mathrm{gxx}0}$) against time is depicted in the inset. The rapid growth behavior of all of the curves seems to happen at about the same time for both s(t) and ${\overline{\mathrm{R}}}_{\mathrm{gxx}}\left(\mathrm{t}\right)$, regardless of the molecular weight. Moreover, sufficiently long chains behave alike, and their curves coincide with each other, revealing the asymptotic behavior of large N. According to ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$, $\mathsf{\sigma }\left(\mathrm{t}\right)$, and s(t), at a given strain rate, the onset of strain hardening is insensitive to the molecular weight but very sensitive to the strain rate. For sufficiently long chains, the local environments (entanglement networks) look alike and cannot be distinguished easily. Therefore, the onset of strain hardening (weakly stretched) and the non-linear viscoelastic response (intermediately stretched) are similar. However, the degree of stretching, which varies with the molecular weight, is gradually revealed as the process continues. As a result, the onset of strain hardening collapses at short times, while the deviation from the asymptotic behavior of infinitely long chains appears at long times.

4. Conclusions

In this work, the microstructural evolutions of the linear polymer melt during the stretching processes are explored by dissipative particle dynamics. Consistent with experiments, the elongational viscosity ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ at different ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$ initially overlaps with the linear response ${\mathsf{\eta }}_0\left(\mathrm{t}\right)$, but it upturns at a shorter stretching time for a higher ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$. The onset of strain hardening is microscopically examined by the radial size distribution P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$). It is intuitively expected that P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) grows widely and shifts toward larger values of ${\mathrm{R}}_{\mathrm{g}}$ upon stretching. However, it is interesting to find that P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) remains essentially the same as that at equilibrium before strain hardening but is strongly deformed after that. A bimodal distribution can emerge after prolonged stretching, suggesting that not all the polymers are elongated. Some polymers have relaxed back into a weakly stretched state, while some are arrested in a specific stretched state. In contrast to the unsteady stretching of long chains, the steady elongational viscosity can be acquired for short chains. The P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) of short chains is weakly perturbed from their equilibrium distribution, revealing that the relaxation is fast enough to recover the near-equilibrium conformations. Nonetheless, for a sufficiently small ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$, P(${\mathrm{R}}_{\mathrm{g}},\mathrm{t}$) is always the same as the equilibrium distribution, indicating the absence of strain hardening.

The dynamics of the size distribution can be handily captured by the mean ${\overline{\mathrm{R}}}_{\mathrm{g}}\left(\mathrm{t}\right)$ and standard deviation $\mathsf{\sigma }\left(\mathrm{t}\right)$. The upturns of both ${\overline{\mathrm{R}}}_{\mathrm{g}}$ and $\mathsf{\sigma }$ from their equilibrium values signify the onset of strain hardening. However, the deviation of the aspect ratio ${\overline{\mathrm{R}}}_{\mathrm{gxx}}{/\left[\right(\overline{\mathrm{R}}}_{\mathrm{gyy}}{+\overline{\mathrm{R}}}_{\mathrm{gzz}})/2]$ from unity is observed before strain hardening, indicating that the initial elongation tends to align the slightly stretched polymers with the stretching direction without altering the equilibrium distribution. In other words, the increment of the elongational viscosity associated with the linear response corresponds to the alignment of stretched polymers by maintaining the equilibrium size distribution. The strain rate dependence of strain hardening can be realized as ${\overline{\mathrm{R}}}_{\mathrm{g}}/{\overline{\mathrm{R}}}_{\mathrm{g}0}$, and ${\mathsf{\sigma }/\mathsf{\sigma }}_0$ are plotted against the dimensionless time ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}\mathrm{t}$. The close coincidence of the onset points suggests that the linear viscoelasticity always takes place at a sufficiently small Hencky strain. In addition to ${\dot{\mathsf{\epsilon }}}_{\mathrm{H}}$, the viscosity depends on the chain length N, but the ${\mathsf{\eta }}_{\mathrm{e}}\left(\mathrm{t}\right)$ of different N behave similarly in the initial period. On the basis of ${\overline{\mathrm{R}}}_{\mathrm{g}}/{\overline{\mathrm{R}}}_{\mathrm{g}0}$ and ${\mathsf{\sigma }/\mathsf{\sigma }}_0$, it is somewhat surprising to find that the onset of strain hardening is insensitive to the chain length. This is because the local environments (entanglement network) look alike for sufficiently long chains, and the weakly stretched state associated with the onset of strain hardening cannot be distinguished easily. Our study of the microstructure dynamics of polymer melt during stretching can serve as a crucial guide for mastering the art of polymer processing.