Stress-Structure Relationship of the Reversible Associating Polymer Network under Start-up Shear Flow

We adopt Langevin dynamics to explore the stress-structure relationship of telechelic reversible associating polymer gel during startup shear flow, with shear strengths varying from Wi = 12.6 to Wi = 12640. At weak shear flow Wi = 12.6, the shear stress proportionally increases with shear strain at short times, followed by a strain hardening behavior and then passes through a maximum (σmax, γmax) and finally decreases until it reaches the steady state. During the evolution of stress, the gel network is only slightly broken and essentially maintains its framework, and the strain hardening behavior originates from the excessive stretching of chains. On the other hand, the stress-strain curve at intermediate shear flow Wi = 505.6 shows two differences from that at Wi = 12.6, namely, the absence of strain hardening and a dramatic increase of stress at large strains, which is caused by the rupture of gel network at small strains and the network recovery at large strains, respectively. Finally, at very strong shear flow Wi = 6319.7, the gel network is immediately broken by shear flow and the stress-strain curve exhibits similar behaviors to those of classical polymeric liquids.


INTRODUCTION
Reversible associating network gels, also known as physical gels, have attracted extensive attention in recent years. A transient gel network is composed of crosslinks that are formed by noncovalent bonds, such as hydrogen bonding, ionic bonds, metal complexation, [1] and π-π stacking. [2−7] Because of their unique abilities such as the outstanding processing and recycling resulting from the nature of reversible cross-linking, physical gels exhibit superior properties over traditional polymer materials, with numerous applications in fields such as biological and industrial technology and engineering artificial tissue. For example, the associated polymeric networks can be used in shape memory, [8−12] drug release, [13] response stimulation [14−17] and self-healing materials. [18−24] In addition, many researchers also focused on polymer hydrogel designs that can absorb water or biological fluids to address medical needs. [25−28] A polymer association gel network is composed of reversible cross-linking multifunctional groups which can be classified into two kinds, i.e., the associating sites distributed along the backbones [29−31] and those located at the ends of polymer chain (telechelic polymers). [32,33] For telechelic polymers, the terminal monomer tends to aggregate with others in the same or different chains into a cluster that effectively works as a crosslink of the physical gel network. Thus, four types of polymer chains are produced in the gel network, namely, bridge, loop, dangling and free chains. The gel network exhibits rich structural and rheological properties under various conditions due to its nature of reversible association properties. Generally, it is believed that the bridge chains play important roles in the elasticity of gel networks. [34−36] Extensive researches by experiment, [37,38] simulation [39−42] and theory [37,43−49] have been performed to understand the structural and rheological properties of gel network. Tanaka and Edwards [46] performed pioneering work in developing theories to explain the complex behaviors of associating polymer gels, such as the mean field theory, percolation network and micelle theory, which describe the assembly and deformation of the networks. They noted that the lifetime of the gel depends on its strength and stability, which are affected by the association group type, location and size. Many other studies examined the dependence of mechanical properties on the structure of the polymer gel networks. [21,38,50] It is found that the bridge chains are responsible for the mechanical strength of the gel network, and chains that have only one end attached to the network (dangling) and free chains that do not form a network cannot withstand stress. In recent years, some theories have been developed to depict the rheological properties of transient gel networks. Green and Tobolsky [43] proposed a transient network theory that asserted that physical crosslinks can break and re-associate by the thermal motion of polymers. Tanaka and Edwald [46] developed an affine deformation theory to explain the dynamic properties and the shear rate dependence of nonlinear viscosity. This theory assumes that a bridge chain deforms affinely with macroscopic deformation. In addition, a nonaffine deformation theory [49] was also established, indicating that the physical crosslinks can fluctuate in space and move from the average positions. Based on this theory, stress overshoot and strain hardening during start-up shear have been investigated. [51] They considered that the strain hardening is caused by the stretching of bridge chains beyond the linear region for a high shear rate. In short, the reversible associating gel network is a soft elastomer. The rupture and formation of the association point leads to the complicated structural rearrangement of the gel network, resulting in the corresponding rheological or mechanical responses. The relationship between the gel structure and macroscopic properties is the key question in this area, which will be illustrated in this work.
The outline of this paper is as follows: The simulation method and molecular model are described briefly in the second section. The stress-strain behaviors and the structural properties of gel network and polymer chains are presented in the third section. Finally, the last section summarizes the findings of this work.

SIMULATION METHODS
In our model system, a telechelic polymer is composed of N beads of mass m. Two beads at the end of a polymer are called stickers. The standard Kremer-Grest model [52] is applied to represent the telechelic polymer chain. The excluded-volume interaction between the stickers is represented by the Lennard-Jones (LJ) potential as follows: The repulsive, shifted, and truncated 6-12 Lennard-Jones potential (known as the WCA potential) is taken into account between the nonstickers: where denotes the distance between the beads i and j, is the Heaviside function (  for  and  for ). The WCA potential is also used to model the nonbonded potential between nonsticker and sticker beads. The standard reduced LJ units, in terms of which all physical quantities are expressed, are used in the simulation, namely m, σ, ε, and (where k B denotes the Boltzmann constant and k B =1), as the mass, length, energy, time, and temperature scales, respectively. The interaction between adjacent monomers within a polymer is expressed by a finitely extensible nonlinear elastic potential: where we use the classic diameters and the fully stretched bond length , guaranteeing a certain stiffness of the bonds while avoiding high-frequency modes and chain crossing. [52] ρ=0.2γ In the simulation, we set the chain length to be 120. There are 1.008×10 5 beads that are confined in a sheared rectangular box with Lees-Edwards boundary conditions and the density is . The box size is L x =140σ, and L y =L z =60σ (where coordinates x, y, and z refer to the flow, vorticity, and gradient direction, respectively). A wide range of shear rates , spanning several orders of magnitude (from 0.0001 to 0.1), are used in these simulations. The strength of the telechelic polymer depends on the mutual interaction between the associative sites. [53,54] We choose a parameter to ensure the formation of a perfect gel network. The association strength affects the lifetime of clusters. The reversible association of stickers gives rise to two characteristic time: the lifetime of being attached to a cluster and the time for a sticker to find a partner of a new cluster . For strong association, it is difficult to break and find new partners for a cluster. In addition to the association strength, chain length and concentration also affect the equilibrium constant of association/dissociation. [55,56] In this work, the interaction of the solvent is expressed by random force and friction. Each monomer in the system is subjected to conservative forces, random forces, and friction forces. The following formula shows the force of each monomer: ∇ν s ∇ν s = (0, 0, 0; 0, 0, 0;γ, 0, 0) where r, ζ, and v refer to the position, friction coefficient, and velocity of the monomers, respectively. Here, refers to the velocity gradient tensor , and we have to deduct the speed of the flow field. The friction constant is set to , and the temperature T is 1.0. We choose the time step to be 0.01. The random force is calculated by this function . To improve the efficiency of the calculation, we utilize GALAMOST [57,58] to perform an equilibrium simulation. First, the pure repulsive potential energy between all the monomers is simulated in 20 million steps, and then the repulsive potential (WCA) between the stickers is replaced by the attractive potential (LJ). Then, 80 million steps are performed to make the system fully reach equilibrium, and a gel network structure is formed. After equilibration, the initial configurations for the nonequilibrium Langevin dynamics simulations are produced from the equilibrium simulations every 10 million steps, start-up shear simulations are performed by LAMMPS. [59] γ=0.1 In the start-up shear simulations, we average a sample within a strain , and the calculated maximum strain is 40. The results are obtained based on sample averages from 30 to 50 in our chosen system. For structural analysis of the gel networks, we utilize cluster analysis [60] to obtain the fractions of four polymer types and cluster sizes. We acquired a transient formation of polymers every 0.1γ. When the telechelic polymers are dissolved in solvents, the stickers located at the extremes of the polymer chains can aggregate into a cluster due to attraction. Within an attractive distance , stickers constitute the same cluster. We judge the type of a polymer by the following rules: if two ends of a polymer are attached to the same cluster or different clusters, this polymer is labeled a loop chain or bridge chain, respectively; a polymer with only one end belonging to a cluster is a dangling chain; a free chain has no end in a cluster. Thus, the sum of all types of polymers is the total. Based on the number of stickers in a cluster, we can count the number of clusters of different sizes.

RESULTS AND DISCUSSION
Stress-Strain Response σ = Gγ When exposed to start up shear, classic polymeric liquids usually behave similarly to elastic solids in short time frames; [61,62] that is, shear stress proportionally increases with shear strain: [63] . However, polymers are viscoelastic, meaning that they exhibit both elastic solid and simple liquid properties. The liquid feature makes the shear stress drop under large strains, which eventually presents a small but nonzero plateau value. Evidently, between the increasing and decreasing regions, a stress peak will emerge and has been widely reported in various polymeric liquids, including solutions, [64,65] melts, [66−68] and gels. [51,69] Typically, this stress overshoot can usually be observed when the imposed shear rate is greater than the reciprocal of relaxation time. [64,67−69] A polymer gel has a backbone connected by telechelic polymer chains, exhibiting a dynamic network under imposed shear. Clearly, its structural responses under shear are more complicated than those in polymer solutions or melts, including both the network and polymer chain levels, which may lead to various stress responses.γ τ pγ γ = 0.0001, 0.004, 0.05 The strength of applied shear flow is characterized by the dimensionless Weissenberg number W i = , where is shear rate and τ p is the longest equilibrium relaxation time of chains by calculating the autocorrelation function of the terminal vector. Here the relaxation time τ p is about 126393τ. To reveal the stress-structure relationship in polymer gels, we first evaluate the strain dependent shear stresses at various shear strengths (W i =12.6, 505.6, 6319.7) corresponding to the shear rates ( ), as shown in Figs. 1(a), 1(b), and 1(c), respectively. σ xz is obtained by the tensor version of the virial theorem: . At shear strength W i =12.6, the shear stress-strain curve shows a similar shape to that of polymer melts, i.e., σ first increases rapidly, then passes through a maximum (σ max , γ max ) and finally decreases until it reaches the steady state (σ steady , γ steady ), with γ steady representing the position where shear stress reaches σ steady for the first time. Correspondingly, we divide the whole curve into three regions I, II, and III. Thus, there exist two characteristic strains: the first characteristic point (σ max =0.028, γ max =5.7) between I and II and the second point (σ steady =0.004, γ steady =22) between II and III. Specifically, we fit the stress-strain curve by an elastic deformation σ = Gγ in re-  gion I, [63] with G representing the linear moduli. As we can see, in the first half of region I, shear stress can be well described by the expression, similar to the behaviors in polymer melts. In addition, in the second half of region I, σ shows an upward deviation from the reference expression, which is usually called strain hardening. [36,70−72] This behavior has not been reported in traditional polymer systems, and its occurrence indicates that the system shows a stronger resistance to shear, i.e., a larger elasticity. This behavior probably originates from the fact that the gel network remains almost unruptured in region I, which causes excessive stretching of the polymer chains. [51,71] We further conjecture that the stress decrease in Region II is caused by network rupture, and the stress plateau in region III denotes network stabilization.
At the intermediate shear strength W i =505.6, the stressstrain curve also shows the three-region behaviors but with two main differences from that at W i =12.6. First, the strain hardening disappears, and the expression σ = Gγ remains valid in almost the whole range of region I. This outcome indicates that the stress increment induced by the nonlinear stretching of the polymer chains no longer dominates but may be counteracted by the stress loss induced by network rupture. The second difference is that the plateau in region III is replaced by a slow increase, which is also found in experiments. [38,73] Consequently, the critical point (σ min =0.01, γ min =7.8) between region II and region III here is defined by the minimum after (σ max =0.028, γ max =2.2), rather than (σ steady , γ steady ). The fact that stress rising in this range does not exist in polymer melts or solutions means its corresponding structural mechanism cannot be rooted in the dynamic responses of polymer chains. Further analysis of the sheared gel network may settle this issue. Note that the special type of stress-strain curve in Fig. 1(b) first appears at a smaller shear strength W i =252.8 and disappears until a shear strength of 1263.9.
Finally, we also show the stress-strain behaviors at very strong shear flow W i =6319.7, as shown in Fig. 1(c). The definition of regions I, II and III is similar to that at W i =12.6. As we can see, the shear stress shows a downward rather than upward deviation from the reference expression in region I, indicating that the rupture to the gel network leads to a dramatic decrease in the network strength, which cannot be compensated by chain stretching. In addition, the stress drop in region II from (σ max =0.053, γ max =1.6) to (σ steady =0.01, γ steady =9.0) is also much faster than that at W i =12.6, from (σ max =0.028, γ max =5.7) to (σ steady =0.004, γ steady =22), probably indicating more serious network rupture at this shear rate. To observe the changes in the viscoelastic properties of the gel network, we can obtain the curves of transient viscosity versus time from formula in Fig. 1(d). These curve shapes are similar to the stress-strain curves. η changes by orders of magnitude depending on the shear strength. In traditional polymer systems, viscosity is independent of shear rates in short time frames; [68] however, in our system, the viscosity is divergent at shear initiation. When the time is long enough, the steady-state viscosity of the system is also different. The viscosity increase was also observed, which echoed the stress-strain curve at the intermediate shear strength.
To further understand the overshoot behaviors, we plot σ max and γ max as a function of W i in Figs. 2(a) and 2(b), respect-ively. As shown in Fig. 2(a), σ max remains almost unchanged when the shear strength W i is below 1263.9. In addition, when the shear strength is larger than 1263.9, σ max experiences a fast increase. Note that a shear strength of 1263.9 is also the critical point between the different types of curves in Figs. 1(b) and 1(c), indicating that the fast increase in σ max appears simultaneously with serious rupture to the gel network. In Fig.  2(b), we show the position of the overshoot point γ max . As seen, γ max dramatically decreases with the increasing shear strength when the W i is below the critical 1263.9 again, followed by a rising trend at large W i . It is well known that in typical polymeric liquids, γ max is located at approximately 2.2 and slowly increases in the range of large shear strengths. [64,66,67,74−76] Clearly, when W i >1263.9, the behavior of γ max in Fig. 2(b) is similar to those of polymeric liquids, once more indicating that the gel network is seriously ruptured and the polymer chain dynamics dominate.

Structural Properties of Gel Network
As discussed above, under shear, two levels of structural responses exist in polymer gels, i.e., network-level and chainlevel responses. To reveal the relationship between the stress response and network evolution, we investigate the details of the network under shear. Generally, the gel network can be characterized by four types of polymer chains: loop, bridge,  Fig. 3(a). Hence, we characterize the structural changes of the gel network by calculating the proportions of four types of chains, i.e., f l , f b , f d , and f f with their values equal to the number of corresponding type chains normalized by the total number of polymers. Note that at equilibrium, f l , f b , f d , and f f have values of 14%, 86%, 0% and 0%, respectively. The absence of dangling and free chains indicates the formation of a perfect gel network. The structural changes of the gel network under various shear strengths are displayed in Figs. 3(b), 3(c) and 3(d), with the discussed shear strengths (W i =12.6, 505.6, 6319.7) corresponding to those in Fig. 1. In detail, at the weak shear strength W i =12.6 shown in Fig. 3(b), it is clear that up to , f l , f b , f d and f f remain almost unchanged from those at equilibrium, indicating that the network is only deformed, and its backbone is maintained in this region. This situation causes the continuous and excessive stretching of chains, closely consistent with the conjecture suggested from the strain hardening behavior in Fig. 1(a). On the other hand, when the shear strain exceeds , the proportion of bridge chains f b slightly decreases with an amplitude of approximately 3%, i.e., from 86% to 83%, while the dangling chains weakly increase, implying slight rupture to the network; this rupture will surely lead to a decrease in stress, which is also consistent with the behavior of region II. Under further large shear strain, the network remains steady, the chains of the four types dynamically evolve with time, and the shear stress develops a plateau region. It should be noted that during the entire shear flow process, loop f l and free f f chains remain unchanged from the equilibrium state. In short, the gel network is only slightly ruptured and remains unbroken at weak shear.  Fig. 3(c), the bridge chains f b considerably decrease at a relatively small shear strain of approximately 1.8 when compared with the weak shear W i =12.6, while the decreasing amplitude can reach 8.5%. It is clear that the excessive stretching of polymer chains will not occur at such a small strain level; hence, the strain hardening behavior in region I disappears, and moreover, it also leads the stress peak (σ max , γ max ) to move to a smaller strain. As expected, the proportion of dangling chains f d correspondingly increases along with the decrease in bridge chains. With the further increase in shear strain, it is found that the proportion of bridge chains f b remarkably shows a successive recovery process rather than remains steady and gradually reaches the value at equilibrium. This process can enhance the strength of the network and eventually result in the stress rising in region III, which is closely consistent with the stress behaviors. Note that the network recovery means the dangling chains f d will decrease, as shown, and during the shear process, the loop and free chains remain almost unchanged with time. In summary, the gel network is first broken at short time frames and then gradually recovers to its equilibrium state at intermediate shear.
At the very strong shear W i =6319.7 shown in Fig. 3(d), the proportion of bridge chain f b starts to decrease immediately after the onset of shear flow, with a large amplitude ranging from 86% to 72%. This serious rupture to the network causes the stress peak (σ max , γ max ) to emerge at a very small strain.  that of the free chains f f dramatically increases, rather than both of them remaining steady with time. The emergence of free chains indicates that the network is broken. At further large shear strains, the four chain types f l , f b , f d , and f f remain steady, and no obvious network recovery is found. In short, the gel network is seriously ruptured and broken under strong shear.
In the gel network, chain end monomers aggregate into clusters of different sizes. To further elaborate the structural properties of the network, we analyze the transient size distribution of these clusters, as shown in Fig. 4. N A denotes the number of stickers in a cluster, and f represents the number of clusters at this size. Again, the shear strengths W i =12.6, 505.6, and 6319.7 indicate weak, intermediate and strong shear flows, respectively. As we can see in Fig. 4(a), at weak shear, the distribution is nearly Gaussian in the absence of shear flow, which can be described by . At small strains such as in region I, the distribution function f remains unchanged from the equilibrium state, indicating a perfect gel network in this range of shear once again. With increasing shear strain, f becomes more expansive, while its shape remains almost unchanged. According to our results, even in the very large strain , the number of chain end monomers that anchor in the gel network is nearly the same as that in the absence of shear flow (approximately 1680), which means that the network is unbroken all the time. On the other hand, the evolution of the distribution function at the intermediate shear strength is quite different. As shown in Fig. 4(b), the number of clusters that have a small size, such γ = 3.0 γ = 30.0 as N A =1, 2, or 3, first undergoes an increase and then drops at large strains. Correspondingly, the number of end monomers belonging to the gel network drops to 1615 at from the equilibrium value 1680 and then recovers to 1650 at large strains . These results confirm that the network is partly broken and finally recovers, as predicted.
At strong shear flow, our results have already indicated that the network is broken, and dangling and even free chains begin to appear in the system. As expected, the data in Fig. 4(c) show a dramatic growth in the quantity of clusters with small sizes, exhibiting a bimodal distribution. Note that no obvious decay of the first peak is observed, even in the largest strain studied, confirming that the network is broken and can no longer recover.

Structural Properties of Polymer Chain
As we stressed before, under shear, two levels of structural responses exist in the gel solutions, i.e., network level and chain level responses, both of which can deeply determine the macroscopic stress behaviors. Hence, following the network analysis, we investigate the structural properties of the chains, including chain stretching and orientation behaviors. In Fig. 5, we calculate the mean-square end-to-end distances of polymer chains. Fig. 5(a) graphically illustrates the definition of end-toend distances. Because the end-to-end distance of the loop chain is always zero and the proportion of free chains is very low, we only show the results of the other two main kinds of chains, i.e., bridge and dangling chains, denoted as and , respectively. In addition, the average end-to-end distance for all the polymer chains is also presented. At the weak shear shown in Fig. 5(b), and are displayed, while is ignored because of the absence of dangling chains. The evolutions of and with increasing shear strain are essentially synchronous since the major contribution of comes from the bridge chains. It can be found that at the small strains, polymers show affine deformation with shear due to the constraint of the network; in other words, and increase rapidly. Then, and exhibit overshoot behavior and begin to decrease, indicating a chain retraction caused by the rupture of the network. At large strains, and tend to level off. In this region, bridge chains are not permanently anchored in the network but dynamically switch between each other, resulting in a reasonable degree of chain stretching rather than excessive stretching. At intermediate shear strength, and , as well as the , are presented in Fig. 5(c). As shown, first undergoes an overshoot and then begins to decrease due to the rupture of the network. With the further increase of strains, shows an evident increment rather than the plateau behavior at weak shear, originating from the recovery process of the gel network. Because bridge chains dominate in gel solutions, it can be expected that the behavior of is similar to that of . In addition, the dangling chains show a plateau region at large strains after the overshoot behaviors because of the existence of chain free ends. Finally, the chain stretching at strong shear flow is shown in Fig. 5(d). As we can see, reaches its peak value at a very small strain, confirming that the network is broken at the onset of shear flow. In addition, the overshoot point of appears at a much larger strain since the dangling chains can be further stretched after the rupture of the network, while the bridge chains cannot be stretched. For both and , the curve tends to level off at large strains, exhibiting a plateau region. Note that even though the bridge chains predominate, the behavior of is similar to that of rather than , since the former is several orders of magnitude larger than the latter.
The orientation behaviors of polymer chains under shear flow are presented in Fig. 6. Symbol θ refers to the angle between the projection of the unit eigenvector of the chain on the flow-gradient plane and the direction of the flow field, as denoted by the inset. Angle θ can characterize the degree of chain orientation along the flow direction. As we can see, the directions of the polymer are isotropic in the absence of shear flow and rapidly orient into the flow direction with the application of shear flow. At weak shear flow, θ gradually de-   creases from 45° at equilibrium to a plateau value of 8° at large strains. At an intermediate shear strength W i =505.6, θ drops much faster than that at weak shear, undergoes an undershoot and then climbs to a plateau value of approximately 12°. The undershoot and the subsequent rising behaviors indicate a stronger resistance to the applied shear flow, originating from the recovery of the network. At a strong shear strength, θ rapidly drops to its minimum at approximately , and no obvious recovery process is observed.
In this work, we studied a reversible association gel network under start-up shear by adopting Langevin dynamics to explore the relationship between the stress response and the structural properties of the gel network. Various shear rates from 10 −4 to 10 −1 are studied, and three typical shear strengths of W i =12.6, 505.6, and 6319.7 corresponding to the shear rates , 0.004, and 0.05 are chosen to represent weak, intermediate and strong shear flows, respectively. The stress-strain curves, as well as the structural properties that include information of two levels, i.e., network-level and chain-level, are discussed in detail. The main conclusions of this paper are as follows: at weak shear flow, elastic deformation is first observed at the onset of shear flow, i.e., shear stress linearly depends on shear strain, followed by a novel region of strain hardening that originates from the excessive stretching of polymer chains. With a further increase in shear strain, the network begins to break, and the stress reaches its maximum and then drops rapidly. A plateau region of stress appears in the large strain limit due to the steadiness of the gel network and polymer chains. At intermediate shear flow, strain hardening disappears due to the breakage of the gel network at relatively small strains, and the plateau region in large strains is replaced by a rising region that is caused by the recovery process of the gel network. At strong shear flow, the serious rupture of the network leads the stress overshoot point to move to a very small value of strain, and a plateau region similar to that of classic polymer liquids is observed in the large strains.
For polymer gels, the association strength is a key factor that can determine both the structural and stress responses in shear flow. To precisely estimate its effect on those proper-ε = 4.0 ties, another gel solution with an intermediate association strength of is studied. It is found that the association strength only quantitatively affects the stress responses and the structural properties under shear flow. Most of the reported results, including the stress rising at intermediate shear flow, as well as the corresponding evolution of network and polymer chains, such as the rupture or recovery of network and the stretching or retraction of polymer chains, are reproduced. In addition to the association strength, there are many other factors, such as chain length and solution concentration, that can affect the structural, rheological and dynamic properties of polymer gels, which will be discussed in future work.