A slip-spring simulation model for predicting linear and nonlinear rheology of entangled wormlike micellar solutions

Surfactant solutions form multiple self-assembled structures depending on composition, as have been intensively studied [1]. One of these structures, that of wormlike micelles, shows highly viscoelastic rheology, similar to that of polymer solutions. These properties are important in industrial and consumer product applications and have been the subject of many experimental and theoretical studies to elucidate the relationship between microscopic micellar parameters, including micelle length, breakage time, and stiffness, and macroscopic properties, primarily rheology [2,3].

One of the important findings is that semidilute entangled wormlike micelles show a single Maxwell exponential stress relaxation time over a range of surfactant and salt concentrations and temperatures [4–8]. In the nonlinear viscoelastic region under the shear flow, wormlike micellar solutions also have distinctive rheological properties, such as shear thickening, particularly, at low surfactant concentrations [9] and shear banding at higher concentrations in the semidilute and concentrated regions [10]. These nonlinear phenomena have been investigated by various experimental techniques, including flow birefringence [11–13] and neutron or light scattering [14–16]. In fast extensional flows, breakage of wormlike micelles leads to filament rupture and extension thinning in the steady extensional viscosity [17,18]. While both shear and extensional rheological phenomena of wormlike micelles are understood qualitatively to some extent, quantitative predictive ability lags far behind that available for flexible polymers.

This gap in understanding is much less severe for the linear rheological properties of wormlike micellar solutions, thanks to the intensive work over the past 30 years [19]. In particular, a pioneering theory of the linear rheology of wormlike micellar solutions by Cates [20] combines the reptation theory for entangled polymer melts [21] with breakage and rejoining mechanisms. In the fast breakage limit, this reptation-reaction model can successfully reproduce the single-relaxation-time Maxwell behavior seen in experiments. Zou et al. extended this model outside of the limitations of the Cates analysis by developing a coarse-grained simulation model for linear rheology that they called the “pointer algorithm.” [22–24] In the pointer algorithm, breakage and rejoining events are incorporated into the reptation theory, and boundaries between relaxed and unrelaxed parts of micellar chains are tracked to calculate the relaxation modulus. Additionally, analytical contributions from high-frequency modes including the bending and Rouse modes are added to give the total relaxation modulus. The resulting pointer algorithm can fit the linear rheological properties of well-entangled wormlike micellar solutions across a wide range of frequencies [22,23]. However, it cannot be used to predict nonlinear rheological properties.

For nonlinear rheology, constitutive equations that were originally developed for polymer melts and solutions are frequently employed to predict nonlinear rheological properties of wormlike micellar solutions [25–27]. Although shear banding is qualitatively captured using constitutive equations for polymers, these predictions cannot be quantitatively linked to the microstructure of wormlike micelles, since polymer theories lack the breakage and rejoining processes. To address this limitation, Vasquez et al. developed the Vasquez–Cook–McKinley (VCM) model that represents the solution by two types of dumbbells, long ones that can break and short ones that can rejoin [28]. To reproduce the nonlinear behavior of wormlike micellar solutions, stress-induced micelle breakage is introduced, allowing the VCM model to qualitatively represent the nonlinear rheology, in particular, shear thinning, in homogeneous flows [29]. Inhomogeneous flows have also been investigated using the VCM model [30,31]. This very oversimplified VCM model could, in principle, be made more accurate by extending it to multiple species. Some basic problems with microphysics in the VCM model, including its thermodynamics, were brought to light by Adams et al., who rederived the VCM model using the kinetic theory, calling the newly derived model the “VCM-R model” [32]. Another serious problem in the VCM model, which cannot be easily repaired, is that it neglects the effect of entanglements. While entanglements were included in a theory for the nonlinear rheology of wormlike micelles by Cates et al. [33,34], this model is restricted to the “fast breakage” regime with a single relaxation time and neglects chain stretch. Although it can reproduce well the shear thinning and shear banding behavior of entangled wormlike micelles, it cannot deal with fast flows, where micelles are stretched, nor with solutions outside of the “fast breakage” limit. To the best of our knowledge, the only other model that addresses the effect of entanglements was developed by Padding et al. based on the “Twentanglement technique” [35]. In the Twentanglement technique, entanglements are detected when two bonds cross and this bond crossing is then forbidden. While rheological properties of wormlike micellar solutions are successfully predicted by this technique, only a few entanglements per chain could be considered due to the computational cost. Thus, a coarser-grained simulation model that can access larger time and length scales is desirable [36]. As has been shown for polymeric fluids, a promising way to take the effect of entanglements into account is to use recent slip-link and slip-spring models [37–42]. In this study, we account for the effect of entanglements by modifying the slip-spring model developed by Likhtman [40] to allow chain breakage and rejoining and apply it to both linear and nonlinear rheology of wormlike micellar solutions.

The content of this paper is organized as follows. In Sec. II, we describe the single-chain slip-spring model (Sec. II A) and our addition to it of a breakage/rejoining algorithm (Sec. II B). In Sec. III, we present the results of our simulations for both linear and nonlinear viscoelastic regimes. In particular, in the linear viscoelastic regime, we compare our simulation results to predictions of the pointer algorithm and to the experimental data. Finally, we present conclusions in Sec. IV.

In this subsection, we briefly explain the single-chain slip-spring model originally proposed by Likhtman [40], which is schematically illustrated in Fig. 1. Although the model is described in Likhtman’s work, we review it briefly here. In the slip-spring model, a polymer chain $i$ $(1≤i≤Nchain(t))$ consists of $N(i)(t)+1$ beads connected by $N(i)(t)$ springs, where $N(i)(t)$ is the number of springs in the polymer chain $i$⁠, $Nchain(t)$ is the number of polymer chains in the system, and $t$ is time. While in the original model, the number of springs per chain is constant, in our modified model, the number of springs of each chain $i$ changes in time due to breakage and rejoining events, as explained in Subsection II B. Therefore, we use the superscript $(i)$ to distinguish different polymer chains. The conformation of a chain $i$ is expressed by the positions of each bead $j$ given by $rj(i)$ $(0≤j≤N(i)(t))$⁠. We use the Hookean law for the springs making these Rouse chains. The stiffness of each spring is obtained by taking the spring to represent a single Kuhn segment with equilibrium root mean square length $b$⁠, and each bead is given a friction coefficient $ξ$⁠. In the wormlike micelle solutions, the characteristic length over which the thermal motion of micelles becomes dominant compared with the local rigidity is called the persistence length (⁠$ℓp$⁠). When the length of micelles becomes longer than the persistence length, wormlike micelles behave like a flexible polymer chain [19]. Then, based on the wormlike chain model [43], the Kuhn length of wormlike micelles is twice the persistence length (⁠$b=2ℓp$⁠). In the micelle slip-spring model, each spring is taken to have a single Kuhn segment with equilibrium root mean square length $b$⁠, which is the typical choice made in the slip-spring model.

In addition to these “main” Rouse chains, the effect of entanglements is introduced through a set of additional Hookean springs that impose entanglement constraints on the main Rouse chains. Each slip-spring $k$⁠, represented by the jagged red line in Fig. 1, is bounded at one end by a “slip ring,” which is the open red circle, and at the other by an “anchoring point,” represented by the solid red square in Fig. 1. The slip-ring $k$ on the main chain $i$ has a time-dependent location $sk(i)$ (⁠$1≤k≤Zss(i)(t)$⁠) along the main chain, where $sk(i)$ can vary between $r0(i)$ and $rN(i)(i)$⁠, where $N(i)$ is the number of springs on chain $i$⁠, and $Zss(i)(t)$ is the number of slip-springs on the polymer chain $i$⁠. The location of the slip-ring $sk(i)$ is determined through a one-dimensional continuous dummy variable $xk(i)$ as

where $trunc(x)$ is the closest integer to $x$ less than or equal to $x$⁠. Here, $xk(i)$⁠, when rounded up to the nearest integer, gives the spring index number on which the slip-spring resides. The slip-spring $k$ on the main chain $i$ is bounded on the other end by an anchoring point with a location $ak(i)$⁠. We also use the Hookean spring law for the slip-springs. Each slip-spring has $Ns$ Kuhn segments and the same Kuhn length $bss$ (⁠$=b$⁠) as in the main chain. Moreover, the initial number of slip-springs on the chain $i$ is determined using the average number of beads between slip-links $Ness$⁠. The slip-springs are initially placed randomly spaced along the chains. Then, we run the simulations to equilibrate their positions.

In the slip-spring model, the states of the polymer chains are updated as follows. The location of a bead $j$ of a chain $i$⁠, $rj(i)$⁠, is evolved as

where $kB$ is Boltzmann’s constant, $T$ is the absolute temperature, $κ$ is the velocity gradient tensor, and $fj(i)(t)$ is the random force. The dummy variable $xk(i)$ (⁠$0≤xk(i)(t)≤N(i)(t)$⁠) is updated as

where $ξs$ is the slip-link friction coefficient and $gk(i)(t)$ is the random force. In this study, we allow slip-springs to pass through each other based on the argument in the original work of Likhtman that the crossability of slip-springs has little effect on the results [40]. On the other hand, others have noted that “uncrossability” of slip-springs may affect chain dynamics [44]. Since this issue is not yet completely resolved even for entangled polymer melts, we simply follow Likhtman's original argument and leave the resolution of this issue to future work. The random forces $fj(i)(t)$ and $gk(i)(t)$ in Eqs. (2) and (3) satisfy the following fluctuation-dissipation expressions:

where $I$ is the unit tensor. Moreover, we assume that the anchoring point of a slip-spring moves affinely with the flow [45] as

In addition to the dynamics described in Eqs. (2) and (3), constraint release (CR) should be considered to correctly reproduce entangled polymer dynamics. In this study, we implement CR by pairing two slip-links between different chains as done in the original work of Likhtman [40]. That is, when a slip-link passes off the end of a chain (⁠$xk(i)(t)<0$ or $N(i)(t)<xk(i)(t)$⁠), this slip-link and the partner slip-link with which it is paired are removed, respectively, from the chain and the partner chain. Simultaneously, two new slip-links are generated as follows: one slip-link is placed at the end of a randomly chosen chain, and the other one, with which the first one is paired, is placed at a random position on a different, randomly chosen chain. Therefore, although the number of slip-springs $Zss(i)(t)$ of a polymer chain $i$ fluctuates in time, the total number of slip-springs remains constant during a simulation. It should be noted that the conservation of the number of slip-springs is only strictly valid for large systems under equilibrium, and in many other slip-link models, the number of slip-links is allowed to decrease under the flow [39,46–48]. However, to keep our model simple, we follow Likhtman’s approach and keep the number of slip-links constant even when we compute nonlinear rheological properties.

From the conformations of the polymer chains and slip-springs, the $α$ and $β$ (⁠$α,β=x,y,z$⁠) components of the real-chain stress tensor and the virtual-spring stress tensor are written as

respectively, where $c$ is the chain concentration and $⟨⋯⟩$ is the ensemble average of $(⋯)$⁠. As pointed out in [49], the virtual-spring stress tensor in Eq. (8) needs to be included in the total stress to correctly calculate the stress relaxation modulus using linear response theory. The proper inclusion of stress produced by virtual springs was extensively discussed by Schieber et al. [50], who showed how to obtain the stress from an expression for the free energy of the virtual springs. However, in this study, following Likhtman’s original approach, we use only the stress from main chain springs and do not use the virtual-spring stress tensor. We do this because the virtual-spring contribution to stress is generally roughly proportional to the main-chain contribution in both linear and nonlinear flows [41,45,51] and because it is difficult to compute accurately the virtual-spring stress in nonlinear flows. We note that the original slip-spring model of Likhtman also omitted the stress from the virtual springs, and we, therefore, use the formula for the modulus from the original paper, which was chosen to fit the experimental data.

Simulations are made dimensionless by dividing stress, time, and length by $G0≡kBT/b3$⁠, $τ0≡ξb2/kBT$⁠, and $b$⁠, respectively. Note that the unit time $τ0$ used in this study differs from that used in Likhtman’s original work $τ0Likhtman$ by a factor of $3π2$⁠. Basic properties of the slip-spring model are summarized in  Appendix A.

The original slip-spring model explained in Sec. II A can successfully predict the dynamics of lightly to well-entangled linear polymer melts by varying the number of slip-links per chain. Before adding breakage and rejoining, as shown in  Appendix A, we demonstrated that our code for the slip-spring model gives predictions for the stress relaxation of polymers nearly identical to those of the original Likhtman model [40]. (Most recent studies of slip-spring models confine slip-springs to hopping from bead to bead rather than sliding between beads as is done in the original model used here [41,42]. We do not employ the hopping version of the slip-spring model because it raises uncertainty about how to address a slip-link that resides at a bead where breakage occurs. We, therefore, use Likhtman’s continuous version in the work presented here.) When we apply this model to the rheology of wormlike micellar solutions, it is essential to combine polymerlike relaxation mechanisms present in Likhtman’s model, such as reptation and chain-length fluctuations, as well as Rouse relaxation modes, with breakage and rejoining characteristic of wormlike micellar solutions. Breakage/rejoining accelerates relaxation of micellar chains by creating new unrelaxed micelle ends, as seen in the pioneering study by Cates [20].

We, therefore, combine the slip-spring model with the breakage/rejoining mechanism as indicated in Figs. 2 and 3. That is, during each simulation time step, we track the motion of main chains and slip-springs according to Eqs. (2) and (3). Whether or not a single breakage or rejoining event occurs for any micelle among all micelles in the system during a time step is determined by the probability $Pbreak/rejoin$⁠,

where $Δt$ is the simulation time step size, $Δtbreak$⁠, and $Δtrejoin$ are the average times between breakage and rejoining events, respectively, and $Δtbreak/rejoin$ represents either $Δtbreak$ or $Δtrejoin$ depending on which event is occurring. Therefore, the statistics of the breakage/rejoining events are Poisson distributed over time $t$⁠.

Since the simulation time step $Δt$ needs to be small enough to give accurate breakage and rejoining probabilities using Eq. (9), $Δt$ should not be larger than the average time required for a single micellar breakage or rejoining event within the entire simulation box (⁠$Δt≤Δtbreak/rejoin$⁠). Using the time constant for a breakage event to occur to any one micelle among all of them, $Δtbreak$⁠, we introduce the average breakage time $τ¯break$ per micelle of average length as

where $⟨Nchain⟩$ is the average number of micelle chains in the system. Unlike the original slip-spring model, the total number of micellar chains in our micelle slip-spring model is allowed to fluctuate since whether breakage/rejoining occurs is determined stochastically according to the probability given in Eq. (9). For efficient simulations, upper and lower limits on the number of chains $Nchain(t)$ in the system, $Nchainmax$ and $Nchainmin$⁠, respectively, are introduced.

Figure 3 illustrates the slip-spring model with breakage and rejoining, which follows the reversible scission scheme assumed in the pointer algorithm [22] (and in Cates’ original model [20]). It is assumed that any main-chain bead on any chain has the same probability to be a breakage point as any other bead, and that bead is duplicated to form a new end bead on each of the resulting chain fragments, which means that the total number of springs is unchanged in a breakage event (⁠$N(i)=N(i′)+N(i″)$ in Fig. 3), while the number of beads increases by one. The extra bead added after a breakage event causes no fundamental problem because beads are inertialess. Of course, the reverse occurs in a rejoining event. From the assumption of equal breakage probability per spring, regardless of the chain length, longer chains are more likely to break. The breakage event is assumed to leave chain and slip-spring configurations and pairing information of slip-springs unchanged as illustrated in Fig. 3. We use a mean-field approach which does not track the relative positions of the micelle centers of mass; thus any micelle can be joined to any other micelle in the ensemble with equal probability. Hence, the rejoining probability is independent of micelle length and, therefore, any two micelle chains have an equal probability to rejoin. Chain and slip-spring configurations and pairing information of slip-springs are unchanged by rejoining as shown in Fig. 3. In a rejoining event, since there are two ends of each chain, four cases can be considered (the following notation is shown in Fig. 3): (1) bead $0$ of chain $i′$ rejoins with bead $0$ of chain $i″$⁠; (2) bead $0$ of chain $i′$ with bead $N(i″)$ of chain $i″$⁠; (3) bead $N(i′)$ of chain $i′$ with bead $0$ of chain $i″$⁠; and (4) bead $N(i′)$ of chain $i′$ with bead $N(i″)$ of chain $i″$⁠. Note that if we define the chain end labeled $0$ as the “head” and the opposite end with nonzero bead number as the “tail,” then chains that do not break or rejoin have head/tail statistical symmetry, since there is no difference, statistically, between the two ends. This statistical head-tail symmetry of each chain is broken in each of the fragments produced by chain breakage: each new chain end created by the breakage comes from the center of the original chain, which have been stretched and oriented more than the ends have been. However, one of these two new chain ends affected by the flow in Fig. 3 is labeled 0 and is a “head” and the other is nonzero and is a “tail.” Hence, once again, the ensemble of fragments retains statistical head/tail symmetry since equal numbers of “heads” and “tails” have been affected by the flow. Thus, each of the four cases mentioned above is statistically identical to each of the others, and ensemble-average properties resulting from each are identical. Therefore, we can safely limit ourselves to only one of the four cases. Here, however, we allow both cases (2) and (3).

During the slip-spring simulation, we maintain upper and lower limits on the number of beads per chain by not allowing the shortest micelles $Nmin$ to break and preventing the longest micelles $Nmax$ from rejoining with any other micelles. We block any rejoining that would create a micelle length greater than the longest one allowed.

In this study, since we restrict ourselves to entangled linear micelles for simplicity, additional “end-interchange” and “bond-interchange” schemes [52], which involve “three-arm” and “four-arm” intermediates, are not addressed. However, the slip-spring model with breakage/rejoining presented in this study could be modified to include them.

Parameters used in this paper are summarized in Table I. In all simulations, we employ $Ness=4$⁠, $Ns=0.5$⁠, and $ξs/ξ=0.1$⁠, which are the same as in Likhtman’s original work [40]. We know of no physical reason for this choice of $ξs/ξ$⁠, other than the observation by Likhtman that $ξs$ cannot be set to zero, and for this value of $ξs/ξ$⁠, the dependence of viscosity on $ξs/ξ$ is small [40]. Resolving whether this choice represents some physical property that might be polymer-specific, or needs to be chosen to counteract some deficiency of the model, remains to be determined. It should be noted that the number of slip-springs is not identical to the number of “entanglements” defined in the tube model. By comparing the slip-spring model with the tube model [53], it has been found that $Netube=7$ should be employed as the number of beads per entanglement if we use the standard slip-spring parameters, $Ness=4$ and $Ns=0.5$ [40]. (More precisely, the number of beads per entanglement when CR is present is $Netube=6.7$ [54].) Then, the number of entanglements can be defined as $Ztube≡N/Netube$⁠. Therefore, the average numbers of entanglements for $⟨N⟩=21$⁠, $35$⁠, $42$⁠, and $63$ are $⟨Ztube⟩=3$⁠, $5$⁠, $6$⁠, and $9$⁠, respectively. The simulation time step used here is $Δt=0.01τ0$ for linear simulations and is $Δt=0.001τ0$ for nonlinear simulations. In previous studies, $Δt=0.01τ0$ [45] and $Δt=0.05τ0$ [42] were employed to predict rheological properties by a discrete slip-spring model in which the slip-links stay on main-chain beads, hopping from one to a neighboring one. Moreover, the virtual-spring stress was not taken into account in those studies. This suggests that the simulation time steps employed in this study are small enough to reasonably evaluate the main-chain stress. However, it should be noted that a smaller simulation time step is needed to correctly address the virtual-spring stress.

Using the average breakage time $τ¯break$ per micelle from Eq. (10), we define the dimensionless breakage time $ζ$ as

where $τ¯rep$ is the average theoretical value of the reptation time defined as $τ¯rep≡3⟨Ztube⟩3τe$ [21]. Here, $τe$ is the Rouse time of the chain segment between entanglements defined as $τe≡(Netube)2τ0/(3π2)$⁠. In what follows, we need to set appropriately the dimensionless breakage time $ζ$ to compare results from slip-spring simulations with those from the pointer algorithm and from experiments. Once $ζ$ is determined, $Δtbreak$ can be calculated using Eq. (10). In all slip-spring simulations including nonlinear cases, we set $Δtbreak=Δtrejoin$ as required by detailed balance, since the average micelle length is fixed by thermodynamics, and we do not consider stress-induced micelle breakage. Whether or not a breakage/rejoining event occurs is determined by Eq. (9) using $Δtbreak$ and $Δtrejoin$⁠.

In all slip-spring simulations, we set the maximum and minimum number of springs per chain as $Nmax=4⟨N⟩$ and $Nmin=2$⁠, respectively, as listed in Table I. The total number of micelles in the ensemble $Nchain(t)$ is allowed to fluctuate within the range $0.95⟨Nchain⟩≤Nchain(t)≤1.05⟨Nchain⟩$⁠. We have checked the effect of the upper bound on the number of springs per chain $Nmax$ and of the imposed upper and lower bounds on $Nchain(t)$ and confirmed that expansion of the range of these parameter values has little effect on $G(t)$ (for more detail, please see the supplementary material [62]).

We here examine basic properties of the slip-spring model with breakage and rejoining. First, all chains in the ensemble are simulated in the absence of the flow (⁠$κ=0$⁠) for a time $teq$ at least two times longer than the theoretically estimated longest relaxation time by Cates $(τ¯repτ¯br)1/2$ [20] when the dimensionless breakage time $ζ<10$⁠, and at least two times longer than the reptation time $τ¯rep$ of a micelle of average length when $ζ≥10$ (note that at high $ζ$⁠, the micelles rarely break).

Figure 4 shows snapshots of a test chain with breakage/rejoining at three different times after an equilibrium simulation. Here, we set the average number of springs per chain as $⟨N⟩=21$⁠, the average number of beads between slip-links as $Ness=4$⁠, and the dimensionless breakage time as $ζ=1$⁠. We can clearly observe that both the number of the test chain beads $Ntest(t)$ and slip-springs $Zsstest(t)$ fluctuate in time. $Ntest(t)$ changes because of breakage/rejoining events, while $Zsstest(t)$ changes because of both CR and breakage/rejoining events.

With the same parameters as used in Figs. 4, 5(a), and 5(b) show the $xx$ and $xy$ components of the stress tensor obtained using Eq. (7) plotted against dimensionless time $t/τ0$⁠, where we remind the reader that the characteristic time $τ0$ is defined as $τ0≡ξb2/kBT$⁠. As expected from Eq. (7), we observe that the $xx$ and $xy$ components of the stress tensor in the absence of the flow are $⟨σxx⟩=⟨N⟩$ and $⟨σxy⟩=0$⁠, respectively, even with breakage and rejoining events. Figure 5(c) shows that the number of chains in the system fluctuates between the upper and lower limits that we set to be $Nchainmax=1.05⟨Nchain⟩$ and $Nchainmin=0.95⟨Nchain⟩$⁠, respectively. Figure 5(d) shows the time-averaged distribution of chain lengths, measured by the distribution in the number of “main-chain” springs in a chain $P(N)$⁠. The red solid line indicates the fitted exponential distribution, which is expected because the chain breakage is a Poisson process with breakage equally likely anywhere along the chain. A slight deviation from the exponential function can be observed especially in the small-$N$ region because of the finite number of chains in the system $⟨Nchain⟩$ and the imposed upper bound on the number of main chain springs $Nmax$⁠.

We next investigate linear viscoelasticity obtained from the slip-spring simulations with breakage/rejoining. Theoretically, the most reliable method to calculate the relaxation modulus $G(t)$ is to employ linear response theory. While $G(t)$ can be calculated from the autocorrelation function of the shear stress at equilibrium, this method takes much longer than the longest relaxation time to obtain an accurate result even if we employ the multiple-tau correlator [55]. Since we employ $⟨Nchain⟩=2000$ wormlike micellar chains, this method is not computationally affordable. A faster method is to obtain $G(t)$ from the relaxation of stress $σxy$ after applying a small step shear strain $γ$ to the sample as $G(t)=σxy(t)/γ$⁠. In this study, we choose $γ=1$ as the step size, which, although not very small, allows one to obtain a reasonable approximation to the linear rheological data with less statistical noise than at smaller strain. (Comparison to the relaxation modulus obtained from the equilibrium simulation is shown in the supplementary material [62].) Previous literature shows that $G(t)$ obtained from a step strain of $γ=1$ is almost identical to that obtained from smaller $γ$ [49]. In the present micelle slip-spring model, rates of breakage and rejoining are taken to be independent of the step size $γ$⁠. Therefore, the step strain does not disturb the distribution functions, such as the number distribution of main chain beads. To further reduce statistical noise, at least 10 independent simulations with different random seeds were performed and the obtained shear stress relaxation curves were averaged. (We have checked that the results are independent of $⟨Nchain⟩$ within the statistical error; this error is larger at long times when fewer chains are used, as shown in the supplementary material [62].) The storage modulus $G′(ω)$ and loss modulus $G″(ω)$ were calculated from the stress relaxation modulus $G(t)$ as

Figure 6 shows the resulting storage (circles) and loss (diamonds) moduli for $⟨N⟩=35$ with (closed symbols) and without (open symbols) slip-springs at three dimensionless breakage times $ζ=0.1$⁠, $ζ=1$⁠, and $ζ=10$ where $⟨Ztube⟩=5$ when slip-springs are present. Here, the dimensionless breakage time $ζ$ for the simulations without slip-springs is defined as $ζRouse≡τ¯break/τ¯R$⁠, where $τ¯R$ is the Rouse time of a micelle of average length defined as $τ¯R≡(⟨N⟩+1)2τ0/(6π2)$⁠. From this figure, we can clearly observe the effect of entanglements and breakage/rejoining events. That is, the presence of entanglements slows the relaxation relative to the same chain without slip-springs, as expected. Moreover, the relaxation becomes faster with decreasing dimensionless breakage time both with and without slip-springs.

Figure 7 shows the linear viscoelastic moduli obtained from the slip-spring simulations for $⟨N⟩=21$ [i.e., $⟨Ztube⟩=3$] at four dimensionless breakage times (a) $ζ=0.1$⁠, (b) $ζ=1$⁠, (c) $ζ=10$⁠, and (d) $ζ=100$ compared to pointer algorithm [22] results for the same parameter values. The pointer algorithm also requires specifying the chain semiflexibility factor $αe$ whose value spans between the loosely entangled regime and the tightly entangled regime and includes the crossover between these regimes. To be precise, $αe$ is the ratio of the entanglement length $ℓe$ to the persistence length $ℓp$⁠: $αe≡ℓe/ℓp$⁠. In the slip-spring model, if we employ the standard parameter set $Ness=4$ and $Ns=0.5$⁠, the semiflexibility factor $αe$ can be determined as follows. Since the number of beads per entanglement length is $Netube=7$⁠, the entanglement length in the slip-spring model $ℓess$ is $ℓess=7b$⁠. Moreover, the persistence length can be determined using the Kuhn length $b$ as $ℓpss=b/2$⁠. Therefore, $αe$ in the slip-spring model is $αe=14$ and, therefore, this value of $αe$ is used in the pointer algorithm results in Fig. 7. In addition, the bending modes usually included in the pointer algorithm are here left out since these modes are left out of the chains in the slip-spring simulations. It should be noted that experimentally most wormlike micelle solutions have $αe$ values smaller than this (typically, $αe$ is close to unity for surfactant concentrations of around $10%$ or so). The relatively high value of $αe$ in the slip-spring model is a consequence of the absence of a bending potential in the bead-spring model. Inclusion of a bending potential in the slip-spring model, while possible, would make the model much more computationally expensive, and is not pursued here. Therefore, when we compare the slip-spring simulation results with the experimental data for typical wormlike micellar solutions, we will need to have a way of adjusting results to correct for the unrealistically high $αe$ in the slip-spring model. Such an adjustment procedure will be discussed later. The unit time and stress of the slip-spring model can be theoretically determined a priori without any adjustments as $τ0=2.6×10−5$ s and $G0=410$ Pa (as described in  Appendix B). These values are used in all linear rheological predictions except for Fig. 10, where comparisons are made to the experimental data. From Fig. 7, the slip-spring simulations with breakage and rejoining and with $αe=14$⁠, which is the same value as in the pointer simulations, are in reasonable agreement with those obtained by the pointer algorithm in the terminal region [except for $ζ=0.1$ shown in Fig. 7(a)]. The difference between the slip-spring and pointer algorithm results for $ζ=0.1$ may be because the number of entanglements is too small for the tube model, on which the pointer algorithm is based, to be accurate and because this inaccuracy becomes more severe at lower $ζ$⁠, where reptation becomes more important. We also observe that the moduli predicted by the slip-spring model are larger than those obtained from the pointer algorithm in the transition region at frequencies above the low-frequency crossover. This difference in the intermediate frequency region is thus related to Rouse modes, as discussed later in Fig. 9.

Figure 8 shows the linear viscoelastic moduli for a larger $⟨N⟩=35$ [i.e., $⟨Ztube⟩=5$] at the same four dimensionless breakage times (a) $ζ=0.1$⁠, (b) $ζ=1$⁠, (c) $ζ=10$⁠, and (d) $ζ=100$⁠, again compared to results from the pointer algorithm at the same parameter values without bending modes. The results are similar to those in Fig. 7 for $⟨Ztube⟩=3$⁠; that is, the slip-spring results are in reasonable agreement with those from the pointer algorithm in the terminal region, but differences are observed in the intermediate frequency region.

The relaxation in this intermediate frequency region should be dominated by the Rouse modes. In the original pointer algorithm, the high-frequency Rouse modes are analytically included as

where $GN$ is the plateau modulus (whose value is discussed in  Appendix B), $ϕi$ is the volume fraction of micelles with length $Li$⁠, $Zi$ is the number of entanglements per chain in the pointer algorithm (equivalent to $Ztube$ in this paper), and $p$ is the mode number of the Rouse motions. In this expression, long-range and longitudinal Rouse modes are assumed to be reduced in magnitude by entanglements so that their contribution is neglected as is evident in the lower bound of the sums in Eqs. (14) and (15), which cut off the small mode numbers $p$ that are less than $Zi$⁠. Instead of including just the high-frequency Rouse modes, in Fig. 9, the solid lines show the predictions that include all the Rouse modes along with the predictions of the slip-spring simulations for (a) $⟨Ztube⟩=3$ and (b) $⟨Ztube⟩=5$ at two dimensionless breakage times (i) $ζ=1$ and (ii) $ζ=100$ while the dashed lines omit the low-frequency Rouse modes. Note that for the “complete” Rouse modes, we do not attenuate the slower “longitudinal” Rouse modes, as is done in many models of entangled polymers [53]. From Fig. 9, the pointer simulation results with all Rouse modes at their full amplitudes gives good agreement with the slip-spring model simulations for both $⟨Ztube⟩=3$ and $⟨Ztube⟩=5$⁠. For the majority of wormlike micellar solutions, which have small $αe$ values and values of $⟨Ztube⟩$ of 10 or greater, the pointer algorithm can give reasonably good rheological predictions in the intermediate frequency regime because the number of Rouse modes is very small. However, when the number of entanglements is relatively small [e.g., $⟨Ztube⟩=3$ and $5$], wormlike micellar chains are weakly entangled so that including only “fast” Rouse modes from tube theory, which describe motion at length scales shorter than a tube segment, are evidently not sufficient. Because the “tube” is not very well formed for small numbers of entanglements, the Rouse modes contribute significantly to relaxation at all frequencies, and one should not remove the low-frequency Rouse modes from the theory nor suppress their magnitude.

We also compare the slip-spring results with the experimental data for a solution of $1.5$ wt. % sodium lauryl ethylene glycol sulfate (SLE1S) and cocamidopropyl betaine (CAPB) with various sodium chloride concentrations at $25$ $°$C [24]. In [24], the linear rheological properties of these wormlike micellar solutions ranging from partially entangled to well-entangled were measured, and we here compare these measurements with the predictions of the slip-spring model with breakage/rejoining. Parameters used in the slip-spring simulations are determined with the assistance of the pointer algorithm as follows. From fits of the pointer algorithm predictions to the experimental data in Figs. 6(b) and 6(c) of [24], the average number of entanglements $⟨Ztube⟩$ and the dimensionless breakage time $ζ$ can be determined as $⟨Ztube⟩=4$ and $ζ=97$ for $[Na+]=0.65$ $mol/L$⁠, $⟨Ztube⟩=6$ and $ζ=20$ for $[Na+]=0.70$ $mol/L$⁠, and $⟨Ztube⟩=9$ and $ζ=3.5$ for $[Na+]=0.75$ $mol/L$⁠, where $[Na+]$ is the salt concentration. It should be noted that the chain semiflexibility factor $αe$ for each salt concentration is around $αe∼9$⁠, somewhat smaller than that in the slip-spring model (⁠$αe=14$⁠), although significantly higher than for more typical concentrations of surfactant near 10%. However, we have found that variations in the value of $αe$ mostly shifts the rheology curves horizontally and vertically without changing very much their shapes. (For more detail, please see  Appendix C.) Therefore, we use following empirical fitting to compare simulation results with experimental ones. Using the parameters determined from the pointer algorithm (⁠$⟨Ztube⟩$ and $ζ$⁠), the slip-spring simulations are performed. By fitting the simulation results with $⟨Ztube⟩=6$ and $ζ=20$ with the corresponding experimental data at the intermediate salt concentration in Fig. 10(b), the two parameters (⁠$τ0$ and $G0$⁠) are determined as $τ0=2.9×10−4$ $s$⁠, and $G0=66$ $Pa$⁠. We show in  Appendix C that these values are not too far from the values (⁠$τ0=1.1×10−4$ s and $G0=72$ Pa) that we estimate by micelle solution properties directly determined from the experimental data, both shifted by around a factor of $2$⁠, as discussed in  Appendix C, to account for the deviation of the true value (⁠$αe∼9$⁠) from the one assumed in the calculations (⁠$αe=14$⁠). For the other sodium ion concentrations (⁠$[Na+]=0.65$ $mol/L$ and $[Na+]=0.75$ $mol/L$⁠) in Figs. 10(a) and 10(c), the same values of $τ0$ and $G0$ are employed. We observe that the slip-spring simulation results (lines) are in reasonable agreement with the experimental data (symbols), including the changing shapes of the curves, except in the high frequency region. In this high-frequency region, the high-frequency relaxation modes, such as bending modes, are important but are not taken into account in the slip-spring model used here. Addressing these high-frequency relaxation modes is one of our future tasks. Another future task is the development of a more rigorous and systematic way of either accounting for micelle stiffness in the slip-spring model, or of shifting results obtained at $αe=14$ to obtain predictions at smaller $αe$ values.

Next, we examine nonlinear rheological properties obtained from the slip-spring simulations under shear flows. In all of the following simulations, the breakage rate and the rejoining rate are set to be constant (⁠$τ¯break=τ¯rejoin$⁠). That is, the rate of breakage does not depend on configuration or stress. It should be noted that stress-induced micelle breakage is modeled by mimicking the idea of the “partially extending convected derivative” [56] in the VCM model [28]. Since the VCM model is based on a binary mixture of Hookean dumbbells with no entanglement effects, the VCM model without stress-induced micelle breakage cannot produce shear thinning. Thus, stress-induced micelle breakage in the VCM model is introduced phenomenologically to produce nonlinear rheological properties, such as the shear thinning. However, the functional form of stress-induced micelle breakage in the VCM model is phenomenological and not physically well founded. In the VCM model with typical parameter choices, accelerated micelle breakage becomes important when the stress reaches the magnitude of the plateau modulus, which is necessary to predict the shear thinning of micellar solutions, which becomes significant at this stress level. Recent molecular dynamics simulations, however, suggest instead that the stress must be more than an order of magnitude higher than the plateau modulus to produce significant acceleration of micelle breakage [57]. In experimental micellar solutions, shear thinning begins at stresses comparable to the plateau modulus of the solution, which scales with the solution concentration to a power of 2. If shear thinning were produced by stress-induced chain scission, as in the VCM model, shear thinning would begin at a stress that scales linearly with concentration, since this would yield a constant tension per micelle, and there would be no relationship between the onset of shear thinning and the entanglement plateau modulus. Thus, the experimental rheological data strongly support the conclusion that shear thinning in micelle solutions is controlled primarily by entanglements and not by stress-induced micelle breakage. Micelle breakage can become important at much higher stresses, $∼0.1$ bar, at $0.1$M surfactant concentration [57]. Thus, the shear thinning in the VCM model is incorrectly attributed to micelle breakage, while in actual micelle solutions, it is likely due to micelle entanglements, which are left out of the VCM model. Since our slip-spring simulations incorporate entanglements, we have no need to introduce micelle breakage to produce shear thinning at the correct stress level, and any consideration of the effect of stress-induced micelle breakage at much higher stresses is, therefore, deferred to future works.

Figure 11 shows transient shear viscosities $η+$ (⁠$≡σxy/γ˙$⁠) for the slip-spring model with $⟨N⟩=35$ chain springs but with no slip-springs, i.e., $⟨Ztube⟩=0$⁠, at two different breakage times (a) $ζRouse=0.1$ and (b) $ζRouse=10$⁠. Here, $γ˙$ is the shear rate and $τ¯R$ for $⟨N⟩=35$ is $τ¯R≃22$ $τ0$⁠. Thus, the shear rates in Fig. 11 correspond to Rouse–Weissenberg numbers, $WiR≡τ¯Rγ˙$ ranging from $0.22$ to $2.2$⁠. The linear viscosity growth function $η0+(t)$⁠, obtained from the linear viscoelastic data (Fig. 6), is also plotted with the black solid line. As expected, transient shear viscosities are almost identical to the linear viscosity growth function $η0+(t)$ at all shear rates. This is because we employ Rouse chains with Hookean springs that cannot show the shear thinning behavior even when we consider breakage and rejoining events if there are no entanglements. Comparing Figs. 11(a) and 11(b), the steady shear viscosity for $ζRouse=0.1$ is smaller than that for $ζRouse=10$ because the increased rate of chain breakage for $ζRouse=0.1$ relative to $ζRouse=10$⁠, which leads to faster relaxation in the former, which can also be seen in Fig. 6.

Figure 12 shows (i) shear viscosities $η+$ and (ii) normal stress coefficients $Ψ1+$ (⁠$≡N1/γ˙2$⁠) for slip-spring model with $⟨N⟩=35$ chain springs and $⟨Ztube⟩=5$ entanglements at two dimensionless breakage times (a) $ζ=0.1$ and (b) $ζ=10$⁠. The linear viscosity growth function $η0+(t)$ and the first normal stress growth function $Ψ1,0+(t)$ [56] obtained from the linear viscoelastic data are also plotted with the black solid lines in Figs. 12(a-i), 12(b-i), 12(a-ii), and 12(b-ii), respectively. The shear rates in Fig. 12 correspond to Rouse-Weissenberg numbers ranging from $0.022$ to $2.2$⁠. (For reptation-time-based Weissenberg numbers defined as $Wid≡τ¯repγ˙$⁠, the shear rates in Fig. 12 correspond to $Wid$ ranging from $0.62$ to $62$⁠, where the reptation time $τ¯rep$ is $τ¯rep≃620τ0$⁠.) From Figs. 12(a-i) and 12(b-i), in the low shear rate region (⁠$γ˙≤0.01τ0−1$ for $ζ=0.1$ and $γ˙≤0.001τ0−1$ for $ζ=10$⁠), the shear viscosities obtained from slip-spring simulations are almost the same as the linear viscosity growth function $η0+(t)$⁠. On the other hand, in the high-shear rate region, we can clearly observe the shear thinning behavior, which comes entirely from the effect of adding entanglements. In Figs. 12(b-i) and 12(b-ii), stress overshoots can be seen in the strain-rate region $γ˙≥0.01τ0−1$⁠. However, in Figs. 12(a-i) and 12(a-ii), stress overshoots are hardly observed because of the fast relaxation from the breakage/rejoining events. A similar result was also seen in a previous study (see Fig. 10 of [35]). In addition to the transient behavior, we show in Fig. 13 the steady state values of (a) shear viscosities $η$ and (b) first normal stress coefficients $Ψ1$ for the slip-spring model with $⟨N⟩=35$ springs and $⟨Ztube⟩=5$ entanglements for two dimensionless breakage times $ζ=0.1$ and $ζ=10$⁠. From Fig. 13(a), we can clearly see the greater shear thinning at a larger breakage time whose viscosity converges with that for smaller breakage time at large shear rates. The shear rate at which convergence occurs roughly corresponds to the region where the strain rates are larger than the inverse of the average breakage time for $ζ=0.1$ [$γ˙≳τ¯br−1(ζ=0.1)$]. In this high-shear-rate region, the effect of breakage and rejoining is no longer detected.

We note that the steady shear viscosity curve is closely related to shear banding. In the current work, we perform homogeneous shear simulations for which shear banding cannot be observed. To perform inhomogeneous shear simulations, the equations of conservation of mass and momentum should be coupled to the slip-spring model [30]. While possible, this would make the simulations much more computationally expensive. Nevertheless, in the future work, the model might be extended to allow shear banding.

We have presented a mesoscopic simulation model for predicting linear and nonlinear rheology of entangled wormlike micellar solutions. Starting from the slip-spring model originally proposed for entangled polymers by Likhtman [40], breakage and rejoining events are added using the reversible scission scheme developed by Cates [20]. With this new simulation model, both linear and nonlinear properties were examined.

In the linear region, results for slip-spring simulations for $⟨Ztube⟩=3$ and $5$ entanglements per chain were found to be in reasonable agreement with those from the coarser-grained pointer algorithm in the terminal region for the same micelle parameters, with no adjustable parameters. However, in the transition region at frequencies above the low-frequency crossover, the moduli from the slip-spring simulations are higher than those from the original pointer algorithm, which includes high-frequency Rouse modes only, omitting those corresponding to length scales longer than an entanglement spacing. Using instead the complete set of Rouse modes gives good agreement between the pointer algorithm and the slip-spring model simulations even at frequencies above the crossover. We conclude that for the rather weakly entangled micelles simulated here, with only $3$ or $5$ entanglements per molecule, the low-frequency Rouse modes are not suppressed enough by entanglements to be ignored, and hence including them leads to much better agreement with the slip-spring model. To perform the comparisons for the same parameter values, we needed to use a rather large value $αe=ℓe/ℓp=14$ of the ratio of micelle entanglement spacing $ℓe$ to persistence length $ℓp$⁠, since the slip-spring model has only been developed for flexible polymer chains, in the so-called “loose entanglement” limit. Typical micelles are semiflexible, however, with values of $αe$ less than $14$⁠. Therefore, to apply the slip-spring simulations to the experimental data, an empirical method of shifting results for large $αe$ along both the horizontal and vertical axes was developed. Using this method, results from the slip-spring simulations were found to be in reasonable agreement with the experimental data of wormlike micellar solutions (⁠$1.5$ wt. % SLE1S + CAPB with different sodium chloride concentrations) [24], except in the high-frequency region. In this high frequency region, the stress relaxation is dominated by high-frequency modes, such as bending modes, which are not considered in the current slip-spring simulations but could easily be added.

In the nonlinear region, we used the slip-spring model to study nonlinear shear rheology of both unentangled and entangled micellar chains with breakage and rejoining events. We assumed that the breakage time is unchanged under flows, thus neglecting stress-induced micelle breakage, which is reasonable except at very high stress. Since we employ Rouse chains with Hookean springs, shear thinning was not observed when entanglements (i.e., slip-links) are absent, and the steady shear viscosity decreases with decreasing micelle breakage time, as expected. For the same chains, with slip-links added, entanglement-induced shear thinning was observed. Furthermore, although clear stress overshoots could be seen for the wormlike micelle solutions with the large breakage time, stress overshoots were hardly observed for the wormlike micelle solutions with small breakage time. This result is consistent with “Twentanglement” simulations of Padding et al. [35].

To make the slip-spring model with breakage/rejoining an effective tool for analyzing rheology of various wormlike micellar solutions, further studies are obviously required. One of the current problems is that the slip-spring model with the standard parameters (⁠$Ness=4$ and $Ns=0.5$⁠) is too flexible to predict rheology in the tightly entangled regime or at the transition between loosely and tightly entangled micelles, which is typical for wormlike micellar solutions. Adding a bending potential to the main chains is one possible, although very expensive, solution for this problem. Another possibility is to refine and make systematic the crude method, introduced here, of shifting data for a high value of the flexibility parameter to more realistic values. After careful parameter optimization using the linear rheological data, the nonlinear rheological predictions from the micelle slip-spring model should be compared with the nonlinear experimental data. The model could also be extended to branched micelles, whose rheology is difficult to predict with mean-field theories. Furthermore, this model can be extended to predict the rheology of associating polymers, whose linear rheology is not well described by the current theory [58]. Our method also has the potential of including other nonlinear effects including stress- or flow-induced micelle breakage. Additionally, the model could be adapted to a multiscale hybrid simulation technique [59] for predicting inhomogeneous phenomena, such as shear banding. We hope to continue our study on the slip-spring model for wormlike micellar solutions along these directions.

T.S. would like to express the gratitude to Professor T. Taniguchi, Professor J. Takimoto, and Professor S. K. Sukumaran for their comments and fruitful discussion on this work. This work was partially supported by the Grant-in-Aid for JSPS Research Fellow Grant (No. 18J10643) and Mazume Research Encouragement Prize in Kyoto University. R.G.L. acknowledges support from the National Science Foundation (NSF) (under Grant No. CBET-1907517). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of NSF.

To validate the slip-spring simulation code, we calculated the relaxation modulus $G(t)$ of monodisperse polymer melts (without breakage/rejoining) for chains with $N=8$⁠, $16$⁠, $32$⁠, and $64$ springs. Here, other parameters are the same as in the original work of Likhtman [40]. Moreover, CR is modeled in the same manner as in the original work. That is, each slip-link has a partner slip-link, with which it remains paired as long as both slip-links survive. When one slip-link passes through the chain end, this and the partner slip-link are removed. Simultaneously, two new slip-links are placed on two randomly chosen chains: one at the end of one of these chains, and the other at a random location on the other. These two slip-links are then partners until they are destroyed. The relaxation modulus $G(t)$ is obtained from equilibrium simulations using the linear response theory as

To efficiently calculate $G(t)$⁠, we use the multiple-tau correlator of Ramírez et al. [55]. Here, we do not consider the stress contribution from the virtual springs [49] because the objective in this subsection is to validate the slip-spring simulation code against the original calculations of Likhtman [40], shown as open symbols in Fig. 14, which also neglect the virtual-spring contribution. Figure 14 shows that slip-spring results agree with those from the previous study.

We have also investigated the rheological properties of the slip-spring model under steady shear flows. We note that there are a few studies that examine nonlinear rheological properties of the polymer slip-spring model with no breakage/rejoining[45,51,60,61]. Figure 15 shows (a) transient shear viscosities, (b) steady shear viscosities, (c) normal stress coefficients, and (d) shear stress and normal stress for $N=64$ (without breakage/rejoining). In (a), the linear viscosity growth function $η0+(t)$ extracted from the linear viscoelastic predictions is also plotted with a black solid line. In (b), the absolute value of the complex dynamic viscosity $|η∗(ω)|$ is also shown. In (c), the first normal stress growth function $Ψ1,0+(t)$ from linear viscoelasticity is also plotted with a black solid line. In (a) and (c), each color corresponds to the same strain rate given in the figure caption. For $N=64$⁠, the Rouse time $τR$ can be calculated as $τR≃71τ0$⁠, and, therefore, the inverse of the Rouse time is $τR−1≃1.4×10−2τ0−1$⁠. Since for the short chains considered here the reptation time is roughly $Ztube$ times the Rouse time, the reptation time is around $1000τ0$⁠, and the nonlinear regime, therefore, begins at around $γ˙τ0∼10−3$⁠. (We note that for very long chains the reptation time converges to $3Ztube$ times the Rouse time.) From (a) and (b), in the small-shear-rate region (⁠$γ˙τ0≤1×10−3$⁠), the shear viscosity is almost the same as the linear viscosity growth function, while at larger shear rates, we can observe shear thinning and a stress overshoot. From (b), in the shear rate region (⁠$γ˙τ0≤1×10−2$⁠), the shear viscosity is almost the same as the complex dynamic viscosity; i.e., the Cox–Meltz rule holds. However, in the larger shear rate region (⁠$γ˙τ0>1×10−2$⁠), we observe that the shear viscosity falls less steeply than the complex dynamic viscosity and the Cox–Meltz rule does not hold. This shear rate region roughly corresponds to the chain-stretching region (⁠$γ˙≳τR−1$⁠). The deviation in this shear rate region occurs because the effect of entanglements is weakened under the fast shear flows [45]. Moreover, since we use the Hookean law for both main chain springs and slip-springs, chain stretch becomes significant and, therefore, the shear stress becomes larger in this region. From (c) and (d), we can confirm that nonlinear rheological properties under shear flows can qualitatively be captured by the current slip-spring model. We note that the shear rates made dimensionless with $τ0$ are all less than unity, and so the “trapping” of slip-links in folded regions of the chain, discussed in Moghadam et al. [61], is not an issue here.