Primitive Chain Network Simulations for Double Peaks in Shear Stress under Fast Flows of Bidisperse Entangled Polymers

A few experiments have reported that the time development of shear stress under fast-startup shear deformations exhibits double peaks before reaching a steady state for bimodal blends of entangled linear polymers under specific conditions. To understand this phenomenon, multi-chain slip-link simulations, based on the primitive chain network model, were conducted on the literature data of a bimodal polystyrene solution. Owing to reasonable agreement between their data and our simulation results, the stress was decomposed into contributions from long- and short-chain components and decoupled into segment number, stretch, and orientation. The analysis revealed that the first and second peaks correspond to the short-chain orientation and the long-chain stretch, respectively. The results also implied that the peak positions are not affected by the mixing of short and long chains, although the intensity of the second peak depends on mixing conditions in a complicated manner.


INTRODUCTION
Viscosity growth under high shear of entangled polymers has been widely investigated as a typical example of the non-linear response [1][2][3][4][5][6].Due to the elastic response, viscosity linearly grows against applied strain in the initial stage, irrespective of the strain rate.When the shear rate is smaller than the slowest relaxation rate (the reciprocal longest relaxation time), viscosity gradually reaches zero-shear viscosity after the longest relaxation time.Under high shear, the steady state viscosity decreases with increasing the shear rate, exhibiting shear thinning.Before reaching the steady state, viscosity exhibits an overshoot.When the shear rate is lower than the reciprocal Rouse time, the peak is located at the total shear strain of ca.2.3, reflecting segment orientation [7].With higher shear, the total strain at the peak increases due to the contribution of chain stretch with increasing the shear rate.In some literature, viscosity exhibits an undershoot following the overshoot, related to the coherent tumbling of molecules [8][9][10][11].Under extremely high shear, strain hardening without a steady state has been reported [12,13].
Concerning the viscosity overshoot, some literature reports multiple peaks rather than the widely observed single peak.Kinouchi et al. [14] reported viscosity growth curves with two peaks for polystyrene solutions with bidisperse molecular weight distributions in which a small amount of long-chain portion was dispersed in short-chain matrices.Osaki et al. [15] performed an extended study that varied the molecular weight of the short-chain matrix.They reported that double-peak behavior is observed only when the molecular weights are sufficiently separated.Snijkers et al. [16] reported multiple peaks for viscosity growth of a comb polymer melt.These multiple overshoots may be attributable to flow instabilities and edge fracture [17][18][19].However, Snijkers et al. [16] performed their experiments with a cone-partitioned plate rheometry [19] to exclude possible effects of edge fracture.
The molecular origin of multiple peaks has yet to be discussed.Concerning bidisperse blends of linear polymers, Kinouchi et al. [14] and Osaki et al. [15] speculated that each mixed portion exhibits its peak.Isram [20] theoretically supported such an idea using his tube model.In particular, he concluded that the first and second peaks reflect overshoots in the stretch of short and long chain fractions, respectively.However, his theoretical model needs to be further considered, particularly regarding the effects of constraint release.In the model, the coupling between different fractions is taken into account by the change of chain contraction through convective constraint release; the relaxation time is modified only under fast flow via chain stretch.In such an implementation, the contribution of segment orientation may need to be further considered.
Nevertheless, the comparison to the data by Osaki et al. [15] was not attained straightforwardly, and no discussion of the segment orientation was given.Other molecular theories that can predict non-linear rheology of mixtures [21][22][23] have not been applied to this problem.This paper examined the double peaks in viscosity growth under high shear for a bidisperse polystyrene solution reported by Osaki et al. [15] employing a multi-chain slip-link simulation [24,25], where the coupling between different components is directly considered.The results demonstrated that the second peak is due to the long chain stretch, as Isram showed [20].
However, the first peak is attributable to the segment orientation, and the contribution from the short-chain stretch is minor.Details are explained below.

MODEL AND SIMULATIONS
Because the model employed in this study is the same as that used in the previous studies for shear flows [25][26][27][28][29][30][31][32], only a brief description is given below.In the primitive chain network model, a network consisting of nodes, strands, and dangling ends represents an entangled polymeric liquid.Each polymer chain in the system corresponds to a path connecting dangling ends through nodes and strands.This path is conceptually equivalent to the primitive path in the tube model [33].
However, no tube is considered in this model, and the path fluctuates due to force balance at each network node and dangling end.Namely, the positions of nodes and dangling ends obey a Langevin-type equation of motion, in which force balance is considered among the drag force, tension acting on each strand, osmotic force suppressing density fluctuations, and thermal random force.Mimicking entanglements, a slip-link is located at each node to bundle two subchains.The chains slide through slip-links, and the chain sliding is described by the change rate equation of the number of Kuhn segments on each strand according to the same force balance for the node position along the chain backbone.Because of chain sliding, when a dangling end protrudes from the connected slip-link beyond a critical Kuhn segment number, a new node with a slip-link is created by hooking another segment randomly chosen from the surroundings.Vice versa, when a dangling end slides off from the connected slip-link, the slip-link disappears, and the bundled chains are released.Concerning simulations under shear flows, infinitesimal step shear strain is applied affinely to the network nodes and dangling ends, followed by relaxations according to the force balance mentioned above.This set of deformation and relaxation is repeatedly applied to attain a designated shear rate.
The length, energy, and time units are chosen as the average strand length under equilibrium a, thermal energy , and diffusion time of the node  !=  " /6, where  is the friction coefficient of the node.For convenience in conversion to the experimental system, instead of  and , units of molecular weight  ! and modulus  ! are used.Here,  ! is the average molecular weight of strands, and  != / # .As discussed earlier,  ! and  ! are similar but different from the entanglement molecular weight  $ and plateau modulus  % .Nevertheless, the parameters  !,  !, and  ! are determined according to the following procedure. ! is estimated from entanglement molecular weight  $ as  != 2 $ /3 according to the empirical relation established earlier [34][35][36][37].Using  !, the average strand number per chain under equilibrium  ! is determined as  != / !, where  is the molecular weight of the examined polymer.A primitive chain network is created according to this  !, and an equilibrium simulation run is conducted for a sufficiently long period, at least ten times longer than the longest relaxation time of the system.From stress fluctuations under equilibrium, linear relaxation modulus () is calculated with the Green-Kubo formula.The obtained () is converted to ′() and "() through the REPTATE software [38,39].These calculations are done with dimensionless units, and the obtained moduli are converted to the experimental value via  !, which is determined from  ! by  != / ! .Here,  is the polymer density.Finally,  ! is determined by fitting ′() and "() to experimental data.
The examined system was polystyrene tricresyl phosphate solutions at 0℃ with a polymer concentration of 0.11 g/cm 3 . !,  !, and  ! were chosen at 10 5 Da, 2.5e 3 Pa, and 0.12 sec, respectively.The molecular weights of short and long chains are  & = 7.1 × 10 ' Da and  ( = 8.2 × 10 ) Da, which were replaced by the chains with the segment number per chain of  & = 7 and  ( = 82 .Hereafter, subscripts L and H denote the low and high molecular weight components; thus, L does not mean the long chain fraction but the short one.According to the experiment [15], the fractions for the long and short chains were 0.09 and 0.91.
The simulations were performed with a cubic simulation box with periodic boundary conditions.The Lees-Edwards boundary was employed for the shear gradient direction for the cases under shear flows.The box size was 16 3 , which is sufficiently larger than the radius of gyration of examined chains under equilibrium.The segment density was fixed at 10. See Fig 1 for a typical snapshot of the single long chain.Because  ! for the examined system is large, finite chain extensibility [40] was not considered.Friction reduction [41,42] was also neglected because the polymer concentration was low.Eight independent simulation runs were conducted for each shear rate, starting from different initial configurations, and quantities reported below are ensemble averages taken among these simulations.

RESULTS AND DISCUSSION
Figure 2 shows the linear viscoelastic response.The simulation reasonably captures the experimental data with two plateaus in ′() corresponding to entanglements between all the chains and only between the long chains.However, the longest relaxation time is overestimated.
Because the primitive chain network simulation has attained semi-quantitative agreement for bidisperse polystyrene melts in previous studies [43,44], the discrepancy may be due to solvent quality that depends on temperature.Nevertheless, the reason is unknown.

Figure 2
The linear viscoelasticity of the examined bidisperse polystyrene solution.Symbols represent experimental data from the literature [15] plotted against the left and bottom axes.Red curves show the simulation results.
In Figure 3, panels (a) and (b) show the time development of shear stress sigma and normal stress difference  * for some shear rates.For low shear rates, both sigma and  * are quantitatively reproduced by the simulation.As the shear rate increases, the discrepancy between the data and the simulation result worsens.Concerning the shear stress shown in panel (a), the time development up to 10 sec is nicely captured even under high shear.However, after exhibiting the first peak, the second peak in the simulation is weaker and occurs earlier than that observed experimentally.Since the double peaks behavior is not apparent in the simulation even for the highest shear rate reported experimentally, the result for a higher shear is added, as shown by the blue curve, which exhibits double peaks.
After bumping, the shear stress reaches a steady value underestimated in the simulation.This discrepancy in the steady value is summarized in panel (c), which shows the steady-state viscosity plotted against the shear rate.The simulation results (red triangle) are smaller than the data (black circle) in the entire range of examined shear rates.However, the simulation is consistent with the Cox-Merz rule; the steady state viscosity agrees with the complex viscosity (red broken curve).
The viscosity obtained experimentally is located beyond the complex viscosity (black broken curve).Nevertheless, the discrepancy in the steady-state viscosity is not significant in the logarithmic plot in panel (c).when the shear rate is high.In a short period of up to 20 seconds, the simulation captures the data irrespective of the shear rate.However, the simulation reaches the maximum at an  * value significantly lower than the experiment.The steady-state value is also underestimated.Note that the blue curve is the simulation result for γ̇=2.5sec - , for which the experiment was not performed, and the coincidence with the data at γ̇=1.165sec -1 (uppermost symbols) is meaningless.
Nevertheless,  * exhibits a single peak even under high shear, for which the shear stress shows double peaks.See blue curves in panels (a) and (b).and red dotted curves are the complex viscosity from the experiment and simulation.
Although the simulation does not excellently reproduce the data, owing to the reasonable agreement, the origin of the double peaks in shear stress is analyzed below for the case of the blue curve at the shear rate of 2.5sec -1 in Figure 3 (a).Figure 4 shows the contributions from the two components for the shear stress and the normal stress difference compared to the response from the entire system (blue) and the results of simulations without mixing.Concerning the double peaks for shear stress seen in panel (a), the first and second peaks come from short (green) and long (red) components, respectively, as suggested by Kinouchi et al. [14] The short-chain contribution shown by the green solid curve coincides with that obtained for the pure short-chain system displayed by the green broken curve, implying that the mixing with the long chain does not affect the short-chain behavior.In contrast, the mixing suppresses the peak stress for the longchain contribution, whereas the peak position is unchanged.See the red curves.For the first normal stress difference, the peak originates from the long component, whereas the contribution from the short chains appears as a shoulder in the short-time region.The effects of mixing in short and long-chain contributions are similar to that observed for shear stress.Green and red broken curves indicate the results for pure long and short-chain systems without mixing reduced according to their mixing fractions.
Further analysis for the molecular mechanism is shown below according to the decoupling approximation [7,28] , in which the stress is approximated as follows.Note that the non-linear spring constant is not considered since the finite chain extensibility is neglected in this study.
Here,  is the stress tensor,  is the number of segments under deformations,  " is the squared segment stretch, and  is the segment orientation tensor.Figure 5 shows these molecular quantities for each component, compared to the cases for the long and short chains without mixing.
Concerning the short chain, all the examined quantities are insensitive to mixing with the long chain, as seen as the coincidence between solid and broken green curves.The peak in shear stress corresponds to the shear component of the orientation tensor  +, shown in panel (c), and it is not related to the stretch  " in panel (b), according to the peak position.The peak of  +, is located at the applied shear strain of ca.2.3, consistent with the tube theory, and no peak is For the long chain, mixing with the short chain suppresses the magnitudes of  reduction and  " enhancement, as seen in panels (a) and (b).However, the characteristic times of their development are unchanged.For instance, the peak position of  " and the mitigation times of  " and  to the steady state are insensitive to mixing.These results are rationalized by the fact that the Rouse relaxation and chain contraction are insensitive to the change in entanglement network.For the orientation  +, shown in panel (c), the peak position is not affected by mixing and is located at the same strain as the short chain.An interesting feature is the undershoot following the peak enhanced by mixing, reflecting the coherent tumbling of the long chains [8,9] .
However, the undershoot does not appear in shear stress because the increase in the segment stretch conceals it.Indeed, the position of undershoot in  +, is close to that in  " .For the orientation  ++ − ,, , no effect of mixing is seen, and thus, the reduction of  * due to mixing in values.However, since it is not significantly greater than for the short chain, the short chain stretch only weakly affects the first peak, which is dominated by the segment orientation.In contrast, the second peak originated by the long chain is well delayed due to the long chain stretch.As a result, the two well-separated peaks are realized.In addition to this first condition concerning the peak positions, the second necessary condition is that the magnitude of peaks must be comparable; if the second peak from the long chain component is too large, the first peak is concealed, and vice versa.The long-chain stress is reduced by mixing, naively due to the reduced long-chain density.
However, mixing also induces suppression of the long chain stretch by culling out the entanglements between long chains, and the corresponding stress overshoot is diminished.See Fig. 4 (a).For the examined case, due to this suppressed long chain stretch, the second peak is somewhat smeared.Consequently, the balance of molecular weights and volume fractions between two components is difficult.The explanation above is consistent with the earlier conjectures [14,15] and the theoretical analysis [20].Nevertheless, the behaviors of molecular parameters apart from chain stretch have been shown for the first time.

CONCLUSIONS
The molecular origin of double peaks in the time development of shear stress under fast shear of the bidisperse polystyrene solution reported by Osaki et al. [15] was analyzed with the primitive chain network simulations.The experimental data for linear and non-linear viscoelastic responses under shear were qualitatively reproduced, the double peaks.Owing to the nature of multi-chain simulation, the entire stress was decomposed into contributions from different components, and the result demonstrated that the first and second peaks come from the short and long-chain components, respectively.Decoupling analysis of stress for each component revealed that the first and second peaks originate from the short-chain orientation and the long-chain stretch, respectively.Although these results are qualitatively consistent with earlier studies, they also imply that the necessary conditions for the emergence of double peaks are not easily fulfilled.
Specifically, the condition is not trivial in realizing comparable peak intensities from the balance in molecular weights and volume fractions between long and short-chain components.This complication is probably the reason why earlier reports for the double peaks are only a few, apart from the experimental difficulties.Further experiments and simulations to clarify such a condition are worth conducting, and the results from supplemental studies will be published elsewhere.

Figure 1 A
Figure 1 A typical snapshot of the system.For clarity, consecutive blue and red cylinders show a pair of long and short chains.Red and yellow spheres indicate the chain ends.Bold green lines show the entangled segments for the indicated chains.Thin blue lines draw other segments.The simulation box is displayed as the yellow frame cube.
The normal stress difference inFig 3 (b)  demonstrates that the data is not correctly reproduced

Figure 3
Figure 3The time development of shear stress (a) and normal stress difference  * (b) for the shear rate of 1.165, 0.58, 0.28, and 0.138 sec -1 from Osaki et al.[15] shown by symbols compared with the simulation indicated by red curves.Blue curves are the simulation results for the shear rate of 2.5 sec -1 , for which the experimental data are unavailable.Panel (c) shows the steady-state viscosity from Osaki et al.[15] and simulations, indicated by the black and red symbols.Black

Figure 4
Figure 4 Shear stress (a) and the first normal stress difference (b) at γ̇ != 0.3.Blue, green, and red curves show the results for the entire system and contributions from short and long chains.
observed for  " .The change in  is relatively minor, slightly decreasing with time.The first normal stress shown in Fig 4 (b) is consistent with the orientation  ++ − ,, and  " , both monotonically increase and reach a plateau around / !~20.

Fig 4 (
Fig 4 (b) originates from the reduction of  " in panel (b).

Figure 5
Figure 5 Decoupled analysis of stress concerning the normalized segment number per chain / !(a), squared segment stretch (b), the shear component of orientation tensor (c), and the difference between normal orientation components of orientation tensor (d).Green and red curves indicate behaviors of the short and long chains, respectively.Solid and broken curves display the cases with and without mixing of short and long chains.The applied shear rate is γ̇ != 0.3.