Nonlinear rheology of entangled wormlike micellar solutions predicted by a micelle-slip-spring model

Surfactant molecules, consisting of a hydrophilic head and a hydrophobic tail, are widely utilized in our daily life [1]. They form self-assembled structures, such as spherical and cylindrical micelles, depending on surfactant and salt compositions. Cylindrical micelles or wormlike micelles (WLMs) with entanglements typically show significant viscoelastic properties and have attracted significant attention from the rheological community.

There have been extensive experimental studies that examine linear rheology [2–4], nonlinear shear rheology [5–12], nonlinear extensional rheology [13–17], and the corresponding microstructures of WLMs [18–20]. In this study, we especially focus on modeling the nonlinear rheological properties of WLM solutions and begin by briefly summarizing the relevant experimental studies. If the WLM solutions remain spatially homogeneous, their rheological properties are similar to those of entangled polymers; that is, the transient shear viscosity growth shows a stress overshoot, and the steady shear viscosity shows shear thinning [5–7]. An important finding is that the Giesekus model, which is a common phenomenological constitutive equation [21], can accurately reproduce the shear rheological properties of WLM solutions [8,9]. An important difference in shear rheological properties between WLM and polymer solutions can be observed, however, at high shear rates, where macroscopic flows can become spatially inhomogeneous. An example of this is shear-banding [22,23], which has been extensively investigated, for example, by Helgeson and co-workers for cetyltrimethylammonium bromide (CTAB) solutions modeled using the Giesekus equation with a stress-diffusion term [10,11]. Extensional flow properties of WLM solutions can be measured by several experimental techniques [24], such as filament-stretching extensional rheometry (FiSER) [13], capillary breakup extensional rheometry (CaBER) [14], and the opposed jet device [15–17]. Rothstein and co-workers reported that the WLM solutions show strain hardening followed by filament rupture before reaching a steady state. This filament rupture was considered to be due to flow-induced scission of micelle chains. Steady extensional viscosities measured by the opposed jet device for several surfactant and salt systems show extension thickening followed by extension thinning [15–17].

To deeply understand such experimental findings, theoretical and numerical approaches are highly desirable. From a theoretical perspective, Cates developed a so-called reptation-reaction model to describe linear “living” polymers [25] that can break at a random point along a chain and whose ends can each randomly fuse with the end of some other chain in a reversible scission scheme. The reptation-reaction model can successfully capture the Maxwell-type relaxation with a single relaxation time as observed in experiments [2]. However, other polymer-like relaxation mechanisms, including chain-length fluctuations and constraint release, which are important to reproduce the rheology of entangled polymers, are not addressed in the original Cates’ model. (In addition to reversible scission, Turner and Cates considered an end-interchange scheme that makes a three-arm star intermediate, and a bond-interchange scheme that makes a four-arm star intermediate [26]. They showed that the relaxation behavior is changed by these mechanisms.)

Based on Cates’ reaction-reptation model, the “pointer algorithm” was developed by Zou and Larson to accurately predict the linear rheology of WLM solutions [27]. In the pointer algorithm, reptation, chain length fluctuations, and constraint release are combined with breakage/rejoining. To compute the relaxation modulus, the boundaries between relaxed and unrelaxed parts of WLM chains are tracked. Analytic formulas for the high-frequency modes are then added to account for the high-frequency region. Zou and Larson showed that this pointer algorithm can reproduce linear rheological properties over a wide frequency range [27]. After their original work, the pointer algorithm was extended in several directions; e.g., the end-interchange and bond-interchange schemes were addressed [28], and relaxation mechanisms of unentangled micelles were included [29]. More recently, Tan and co-workers have established a formula based on pointer algorithm simulations to estimate the micelle length from linear rheological data [30]. Utilizing the pointer algorithm, we can, thus, obtain accurate microscopic micelle parameters from linear rheological data. Compared to these modeling developments in linear rheology, nonlinear rheological models are relatively underdeveloped.

One of the most utilized models for the nonlinear rheology of WLMs is the Vasquez–Cook–McKinley (VCM) model [31]. The VCM model is composed of two species of Hookean dumbbells: long and short dumbbells. A long dumbbell can break to make two short dumbbells and vice versa to mimic the reversible scission scheme. Since the Hookean dumbbell model is employed, the VCM model with a constant breakage rate cannot reproduce nonlinear rheological properties, such as shear thinning observed in WLM solutions. Thus, the VCM model includes the stress-induced micelle breakage under flow [31]. While several nonlinear properties under a homogeneous flow have been successfully reproduced by the VCM model [32], the effects of entanglements (i.e., the constraints they imposed on chain relaxation, and how they might be lost and gained under flow) are not included in the VCM model, which, thus, neglects the most important source of shear thinning in WLM solutions. Shear thinning is predicted by the VCM model, but only through inclusion of shear-induced micelle breakage, whose rate is adjusted arbitrarily to produce shear thinning strong enough to agree with experiments.

Very recently, Peterson and co-workers developed two molecular-based constitutive models for entangled WLM solutions [33–35]. One is simplified tube approximation for the rapid-breaking micelles (STARM) model [33], and the other model is the living Rolie-Poly model [34]. The basic idea of these models is to combine the population balances and polymeric relaxations, i.e., reptation, contour length fluctuation, and constraint release. While these models are physics-based and sophisticated, it is difficult to implement them in macroscopic flow calculations due to their model complexity. It should be noted that there is a simplified version of the living Rolie-Poly model in the fast-breakage limit, where micelle breakage is much faster than reptation, that might be used for macroscopic flow calculations. However, the fast-breakage limit is often not attained in common micellar solutions, and to the best of our knowledge, a simplified physics-based constitutive model beyond the fast-breakage limit has not been developed yet. To avoid complexity appearing in these constitutive models, it may be useful instead to develop an accurate mesoscopic simulation model and fit the predictions of this model with a simpler constitutive equation. This strategy would allow a tractable equation to be developed whose parameters could be correlated with micelle properties, rather than simply being adjusted to fit the experiments, without yielding any connection to microscopic physics.

We developed such a mesoscopic “micelle-slip-spring simulation model” to predict rheological properties of WLM solutions [36]. Our model is based on the slip-spring model proposed by Likhtman [37]. Slip-link and slip-spring models were originally developed for unbreakable entangled polymers and have been found to predict well their linear and nonlinear rheological properties [38]. To address the rheology of WLM solutions, breakage and rejoining events were included in the polymer slip-spring model, and the results were validated through comparison with predictions of the pointer algorithm and with experimental linear viscoelastic data. Moreover, using the micelle-slip-spring model, we examined nonlinear shear properties including shear thinning, which naturally emerge from our micelle-slip-spring model due to the effect of the entanglement constraint [36].

In this article, we further investigate nonlinear rheological properties obtained by the micelle-slip-spring model. For such a purpose, instead of the Hookean springs utilized in our previous study following the original slip-spring model [36], we use the Fraenkel springs to prevent the overstretch of springs so that higher shear rates and extensional flows can be simulated where finite extensibility of the chain becomes important. Additionally, to improve extensional properties under a strong extensional flow, we incorporate stress-induced micelle breakage into the micelle-slip-spring model. Using the extended micelle-slip-spring model, we examine nonlinear shear and extensional properties. After defining and explaining the model in Sec. II, results are presented in Sec. III and summarized in Sec. IV.

Here, we briefly explain the original micelle-slip-spring model [36]. A micellar chain is expressed by a main chain and slip-springs that impose motional constraints on the main chain, as shown in the upper figure in Fig. 1. A micellar chain $i$ $[1≤i≤Nchain(t)]$ consists of $N(i)(t)+1$ beads connected by $N(i)(t)$ springs. Here, the superscript $(i)$ distinguishes different micellar chains, and $Nchain(t)$ is the number of micellar chains in the system at time $t$⁠. Because of breakage and rejoining events, $Nchain(t)$ and $N(i)(t)$ fluctuate in time in our micelle-slip-spring model.

Each main chain spring is modeled by a single Kuhn segment with Kuhn length $b$⁠, and each main chain bead is given a friction coefficient $ξ$⁠. Here, the Kuhn length is related to the persistence length $ℓp$ as $b=2ℓp$ [39]. On length scales less than $ℓp$⁠, the chain is locally rigid.

The main chain $i$ has $Zss(i)(t)$ slip-springs at time $t$⁠. A slip-spring illustrated in Fig. 1 consists of a slip-link shown with the red open circle connected to an anchoring point shown with the red filled diamond by a red slip-spring. The slip-link $k$ on the chain $i$⁠, $sk(i)$⁠, is constrained to move along the main chain $i$⁠. Different from current “discrete” slip-spring models, in which the positions of the slip-links are restricted to jump main chain beads [40,41], in our model, as in the original slip-spring model of Likhtman, the slip-links can move anywhere along the chain [37]. The time evolution of $sk(i)$ is defined by a 1D continuous dummy variable $xk(i)$⁠, which spans between $0$ and $N(i)(t)$⁠. In this study, slip-links are allowed to pass through each other since this has little effect on the results, as noted by Likhtman [37]. While the slip-link can travel along the main chain, the corresponding anchoring point, $ak(i)$⁠, is pinned in space unless the flow is applied. When the shear or extensional flow is applied, the anchoring point moves affinely, as in our previous study [36].

Each slip-spring has $Ns$ Kuhn segments with the Kuhn length $bss$⁠, and each slip-link has a friction coefficient $ξs$⁠. As in our previous study [36] that follows the original study by Likhtman [37], we assume that the Kuhn length of main chain springs is the same as that of slip-springs (i.e., $b=bss$⁠). The slip-link friction coefficient $ξs$ is introduced to control the dynamics of the slip-springs. The initial number of slip-springs per chain is determined by an additional parameter, the average number of beads between slip-links, $Ness$⁠. The slip-links are initially distributed randomly along the chains, and the anchoring points are placed around the corresponding slip-links according to a Gaussian potential of a Hookean spring as in the original model of Likhtman [37] even though we employ the Fraenkel potential. We note that the microscopic detail of the model should not appear in long-time rheological behavior, and thus, the effect of this approximation would be expected to vanish in the long-time region.

The remaining relaxation mechanism needed to reproduce the dynamics of entangled polymers is constraint release (CR). We implement CR as in the original Likhtman’s slip-spring model [37], in which a slip-link is paired with a partner slip-link on a different chain. The CR occurs when the slip-link slides off the end of the chain [i.e., $xk(i)(t)<0$ or $N(i)(t)<xk(i)(t)$]. When CR occurs, this slip-link and the partner slip-link are removed from the chain and the partner chain, respectively. At the same time, two new slip-links are placed on two different randomly chosen chains. One slip-link is placed near the end of a chain, and the other is placed at a random position on a different chain. The corresponding anchoring points are distributed using the same procedure explained previously. Thus, while the number of slip-springs on a given chain $i$⁠, $Zss(i)(t)$⁠, fluctuates in time due to CR or breakage/rejoining events, the total number of slip-springs in the system is unchanged even under flow in most of our simulations. However, the effects of flow-induced reduction in entanglement density are quite possibly important [42] and are, therefore, considered in Sec. III D, where we have examined a limiting case in which the total number of entanglements continuously decreases and eventually becomes zero.

In the original micelle-slip-spring model, we utilized a Hookean spring law following the original polymer slip-spring model. However, in this study, we focus on nonlinear shear and extensional properties of entangled WLMs. Thus, it is important to address the effect of finite extensibility. For springs that each correspond to a single Kuhn length, it is typical to employ the Fraenkel spring law [43]. The total potential energy of the chain $i$ is then expressed as

where $K$ and $Kss$ are the Fraenkel parameters for main-chain springs and slip-springs, respectively; $kB$ is the Boltzmann constant, $T$ is the temperature, $Rj(i)$ is the vector between two adjacent beads defined as $Rj(i)=rj+1(i)−rj(i)$⁠, with $rj(i)$ being the position of the $j$th bead on the $i$th main chain, and $qk(i)$ is the vector between the slip-link and the anchoring point defined as $qk(i)=ak(i)−sk(i)$⁠. Using this potential, we can express the spring force originating from spring $j$ of the main chain, $Fsp,j(i)$⁠, as

Furthermore, the spring force originating from slip-spring $k$⁠, $Fss,k(i)$⁠, can be written as

From these forces, bead and slip-link positions are updated over time. For more details regarding the dynamical equations, please see our previous study [36].

From the forces and conformations of main chains and slip-springs, the $αβ$ [$α,β∈(x,y,z)$] components of the stress tensor of the main chain (i.e., real chain) and of the virtual-spring are expressed as

respectively, where $cchain$ is the chain concentration and $⟨⋯⟩$ is the ensemble average. To compute the relaxation modulus $G(t)$ accurately from an equilibrium simulation using the Green–Kubo formula, Ramírez and co-workers [44] pointed out that the virtual-spring stress tensor $σαβV$ needs to be included so that

where $x$ is the velocity direction, $y$ is the velocity gradient direction, $γ$ is the step strain magnitude, and $V$ is the volume of the system. It should be noted that only the stress from main-chain springs is utilized to obtain $G(t)$ in the original work by Likhtman [i.e., $G(t)=(V/kBT)⟨σxyR(t)σxyR(0)⟩$] [37]. However, this $G(t)$ deviates from relaxation modulus computed by the relaxation after a small step strain, $σxyR(t)/γ$⁠, as explained by Ramírez and co-workers [44]. We note that the virtual stress contribution is not more than around $25%$ of the main-chain contribution and shows similar dependence on time as the main chain contribution under shear flow and is negligibly small under an extensional flow (for more detail, see Sec. S1 in the supplementary material [45]). Thus, because of the rather long runs with small $Δt$ values needed to adequately resolve the virtual-spring contribution and its rather weak effect, we use only the main-chain stress $σαβR$ to compute the nonlinear rheological properties.

Breakage and rejoining events are included in our model as done in our previous study [36]; that is, we use the reversible scission scheme and do not consider the end-interchange and the bond-interchange schemes. To incorporate these events, we introduce the probability $Pbreak$ that a breakage event occurs for any micelle among all micelles during each time step,

where $Δt$ is the simulation time step size and $Δtbreakeq$ is the average time between breakage events at equilibrium. In this study, a rejoining event occurs during each time step with probability $Prejoin$ that is the same as $Pbreakeq$ over all time steps even under flow. At each simulation time step, breakage and rejoining events independently occur with the probabilities of $Pbreakeq$ and $Prejoin$⁠, respectively. Thus, the number of micelles $Nchain(t)$ fluctuates in time. (It should be noted that the breakage probability is increased under the strong extensional flow due to stress-induced micelle breakage, as explained in Sec. II B.) The average breakage time per chain with average length at equilibrium or under weak flow, $τ¯breakeq$⁠, can be defined using $Δtbreakeq$ as

where $⟨Nchain⟩$ is the average number of WLM chains in the system. The dimensionless breakage time, $ζ$⁠, is defined using $τ¯breakeq$ as

where $τ¯rep≡3⟨Ztube⟩3τe$ is the pure reptation time in the absence of other relaxation mechanisms of a micelle of average number of entanglements $⟨Ztube⟩$ and $τe≡(Netube)2τ0/(3π2)$ is the Rouse time of the chain segment between entanglements [46].

The assumptions regarding breakage/rejoining events are as follows. When a breakage event occurs at a certain simulation time step, a randomly selected main-chain bead is chosen to be a breakage point. At the breakage point, the chosen bead is duplicated, and these two beads are made the new ends of the resulting chains, as can be seen in Fig. 1. Since beads are inertialess, there is no fundamental problem with adding a bead. For a rejoining event, any two micelles can fuse with equal probability. That is, we assume the rejoining probability does not depend on micelle length and employ a mean-field approach that does not track the relative positions of micelles. Before and after breakage and rejoining events, configurations of main chains and slip-springs and pairing information of slip-links remain the same. To reduce computational costs, we impose upper and lower bounds on $Nchain(t)$⁠, $Nchainmax$ and $Nchainmin$⁠, respectively. Additionally, we maintain upper and lower bounds on $N(i)(t)$ by introducing the maximum and minimum numbers of springs per chain, $Nmax$ and $Nmin$⁠, respectively.

Thus far, we have explained the case where the breakage and rejoining times are assumed to be constant, which is valid at equilibrium or under a weak flow. Extensional flow experiments by FiSER for the CTAB/NaSal system [13] suggest that a micelle will break when the extensional stress per micelle reaches around $3.6×105$ Pa [47]. This stress is about $104$ times larger than the experimental plateau modulus of $GN0∼10$ Pa for the CTAB/NaSal system [13]. Under a strong extensional flow, the extensional stress obtained from the micelle-slip-spring model is expected to become higher than the micelle breakage stress and stress-induced micelle breakage (SIMB) might become evident. To accurately model the nonlinear rheology under strong extensional flow, we address the effect of SIMB in this study.

There are several ways to address SIMB in our micelle-slip-spring model. A simple way is to set the threshold stress for micelle breakage over which the micelles break. The threshold stress can be estimated from experimental data of extensional stress before rupture [13]. However, the experimental threshold stress might be so high that all micelles break very fast. More realistically, under high stress, the micelles can be allowed to break faster than the time scale $τ¯breakeq$ for breakage at lower stress than the threshold stress. Mandal and Larson developed the activation formula that can address the acceleration of breakage at any stress [47]. In this article, we use this formula to obtain the extensional rheological properties under more reasonable conditions. The average breakage time of the chain $i$ can be expressed as

where $τ°$ is the time constant and $Esciss$ is the scission free energy. The molecular dynamics simulations by Mandal and Larson revealed that the dimensionless scission energy $E~sciss=Esciss/kBT$ decreases linearly with stress. Thus, $E~sciss$ of the chain $i$ can be expressed as

where $E~scisseq$ is the scission energy under no external force, $σ~E≡σE/GN0$ is the normalized extensional stress per chain, and $σ~Ebreak≡σEbreak/GN0$ is the normalized characteristic extensional stress of micelle breakage, which can be estimated from molecular dynamics simulations. It should be noted that $σEbreak$ is not identical to the experimental breakage stress per micelle (⁠$∼3.6×105$ Pa for the CTAB/NaSal system [13]) and is expected to be smaller than the experimental breakage stress per micelle for the reason mentioned above. Combining Eqs. (10) and (11), we can obtain the breakage time under the extensional flow as

where the average breakage time at equilibrium, $τ¯breakeq$⁠, is expressed as $τ¯breakeq=τ°exp⁡(E~scisseq)$⁠.

We incorporate SIMB into the slip-spring model as follows. At each simulation time step, we compute the extensional stress of individual chains, which means that each chain has different extensional stress. This is used to determine whether or not breakage occurs using Eq. (12) and a random number. Namely, the breakage probability shown in Eq. (7) is increased under the strong extensional flow. In general, the tension varies along the chain with the highest tension being near the center of the chain. In this study, however, for simplicity, we use the averaged stress over the whole chain and assume that every bead of the breaking chain has an equal probability to be the breakage point. Thus, the bead at which breakage occurs is randomly selected, as done in our original micelle-slip-spring model. For rejoining events, we assume that the rejoining rate retains its equilibrium value even under the strong extensional flow.

Simulations are carried out using dimensionless variables. Units of stress, time, and length are $Gchain≡cchainkBT$⁠, with $cchain$ being the chain density, $τ0≡ξb2/kBT$⁠, and the Kuhn length $b$⁠, respectively. Here, $τ0$ in this study differs from the Rouse time of a single Kuhn segment by a factor of $3π2$⁠.

Basic parameter values used in this article are the same as those in our previous study [36]. For slip-spring parameters, we use $Ness=4$⁠, $Ns=0.5$⁠, and $ξs/ξ=0.1$⁠, which are the same as in Likhtman’s original work [37]. Here, a comment on the number of “entanglements” needs to be made. The average number of slip-springs introduced to reproduce the entangled polymer dynamics, $⟨Zss⟩=⟨N⟩/Ness$⁠, is not the same as the number of entanglements, $⟨Ztube⟩=⟨N⟩/Netube$⁠, defined in the tube model. Here, $Netube$ is the number of beads per entanglement. Comparing the slip-spring model with the Likhtman–McLeish tube model [48], Likhtman found that $Netube≃7$ gives the number of beads per entanglement when using the standard slip-spring model parameters, $Ness=4$ and $Ns=0.5$ [37]. (More precisely, $Netube$ with and without CR is $Netube=6.7$ and $5.7$⁠, respectively [49].) Thus, we use $Netube=7$ to determine the number of entanglements per chain $⟨Ztube⟩$⁠.

The Fraenkel parameters for main chain springs and slip-springs are both set to $K=Kss=102$ for all simulations. The Fraenkel parameter is arbitrary as long as it is large enough to prevent spring stretch. We chose the above value because it accomplishes this without excessive computational cost required when employing very large Fraenkel parameter values. Our choice of $K=Kss=102$ is of the same order as used in several other studies [50–52], while larger values were employed in one particular study [53]. The simulation time step $Δt$ is kept less than $10−3τ0$⁠. A $Δt$ value smaller than $10−3τ0$ is required to accurately obtain the stress under the strong shear and extensional flow. For a shear rate range of $γ˙>τ¯R−1$⁠, where $γ˙$ is the shear rate and $τ¯R$ is the average Rouse time defined as $τ¯R≡(⟨N⟩+1)2τ0/(6π2)$⁠, we set the time step less than the value of $10−3τ0/(γ˙τ¯R)$⁠. Here, in typical simulations, the time step size should be set to be inversely proportional to the strain rate so that the strain increment in each time step is constant. Furthermore, even for the same strain rate, we used a smaller time step in the extensional flow than in the shear flow. Due to the computational cost, we examined the time step size range of $Δt≥10−5τ0$⁠. As discussed later in Sec. III C, the smallest time step size $Δt=10−5τ0$ might not be enough to prevent overstretch under extension. Nevertheless, we confirmed that the $Δt$ values used in this study are small enough to that the stress from main chains $σR$ does not show significant divergence. We show the effect of changing the time step size in Fig. S2 in the supplementary material [45].

For breakage and rejoining events, $Δtbreakeq$ is an input parameter for simulation that is obtained from the micelle average breakage time $τ¯breakeq$ at equilibrium from Eq. (8). To compare results from slip-spring simulations with those from other models or experiments, it is convenient to use the dimensionless breakage time $ζ$ that can be computed from Eq. (9) once $τ¯breakeq$ and $τ¯rep$ are obtained.

In this study, we mainly set the average number of chains in the system to $⟨Nchain⟩=2000$⁠. For several shear and extensional simulations with large strain rates, which require a small $Δt$⁠, we set $⟨Nchain⟩=1000$⁠. Under the strong flow, since the fluctuation in the stress becomes small, we can safely use the smaller number of chains. We set the maximum and minimum numbers of springs to $Nmax=4⟨N⟩$ and $Nmin=2$ in all simulations. The total number of micelles $Nchain(t)$ without SIMB is allowed to change within the range of $0.95⟨Nchain⟩≤Nchain(t)≤1.05⟨Nchain⟩$⁠. The above values are the same as in our previous study [36]. We have checked the effect of changing the ranges of $N(i)(t)$ and $Nchain(t)$ and confirmed that the expansion of these ranges has little effect on $G(t)$ [36]. On the other hand, there is no upper limit to $Nchain(t)$ for simulations with SIMB.

Since the experimental shear rheological properties of WLM solutions can be successfully fitted by the Giesekus model [8,9], we will see the Giesekus model can also fit the micelle-slip-spring simulation results. If we can systematically map micelle parameters of the slip-spring model (e.g., the number of entanglements and the breakage rate) onto a nonlinear parameter appearing in the Giesekus model, physical insights could thereby be indirectly incorporated into macroscopic (inhomogeneous) flow simulations, which might give better understanding of inhomogeneous flows of WLM solutions and how they depend on micellar parameters, such as the breakage rate, rather than merely on the phenomenological nonlinear parameter of the Giesekus model.

We briefly explain the Giesekus model and its steady state solutions under shear and extensional flows. The Giesekus model is a well-known phenomenological constitutive equation [21,54], which for mode $p$ (⁠$1≤p≤Nmode$⁠) is expressed as

where $σp$ is the stress for mode $p$⁠; $σ▽p$ is the upper convected derivative defined as $σ▽p=σ˙p−σp⋅κ+−κ⋅σp$⁠, with $σ˙p$ being the time derivative of $σp$ and $κ$ being the velocity gradient tensor; $αG$ is the anisotropic parameter, which is the same for all modes [55]; $I$ is the unit tensor; and $τp$ and $Gp$ are the relaxation time and modulus for mode $p$⁠, respectively. Here, $τp$ and $Gp$ are determined by the multimode Maxwell model fitting of the relaxation modulus. The total stress $σ$ is $σ=∑p=1Nmodeσp$⁠.

In the shear flow, the steady shear viscosity for mode $p$⁠, $ηp$⁠, is written as [56]

where $η0,p$$(=Gpτp)$ is the zero shear viscosity for mode $p$⁠, and $fpS[γ˙,αG]$ is the function defined as

with

Here, $γ˙$ is the shear rate. Under the extensional flow, the steady extensional viscosity for mode $p$⁠, $ηE,p$⁠, is [56]

with

and

where $ϵ˙$ is the extension rate. The anisotropic parameter $αG$ is determined by fitting steady shear or extensional viscosity data.

In addition to the comparison of slip-spring simulation results with predictions of the Giesekus model, we will test another common constitutive equation, namely, the Phan–Thien/Tanner (PTT) model [54,57]. The PTT model is expressed as

where $σ□p$ is the Gordon–Schowalter convected derivative defined as $σ□p=σ˙p−σp⋅κ+−κ⋅σp+(1−a)(σp⋅D+D⋅σp)$⁠, with $a$ being the parameter describing a slippage of the strand and $D=(κ+κ+)/2$ being the strain rate tensor, $Y$ is a function of the term in brackets that depends on the average conformation of the polymer chains through the trace of the stress tensor as indicated, and $Gp∗$ is the modulus defined as $Gp∗=a2Gp$⁠. In this study, we set $a=1$ to avoid unphysical oscillations produced by the Gordon–Schowalter convected derivative under the shear flow [54]. Thus, $σ□p$ and $Gp∗$ reduce to $σ▽p$ and $Gp$⁠, respectively. As is the case for the Giesekus model, $τp$ and $Gp$ are determined by the multimode Maxwell model fitting of $G(t)$⁠.

Two functional forms for the dependence of $Y$ on conformation have been suggested, namely, a linear function and an exponential function. The difference between the two forms is evident in the steady extensional viscosity as a function of the extensional rate. For the linear function, the extensional viscosity shows extension hardening followed by a high-extension-rate plateau. For the exponential function, the extensional viscosity shows extension hardening followed by extension thinning. Since typical experimental extensional viscosities obtained by the opposed jet device show extension hardening followed by extension thinning [15–17], we use the exponential function expressed as

where $αPTT$ is the nonlinear parameter in the PTT model. As in the Giesekus model, the nonlinear parameter $αPTT$ is the same for all modes. That is, we use a single nonlinear parameter $αPTT$ to fit the slip-spring data. When we use the exponential function shown in Eq. (21), the steady shear viscosity can be expressed as [58]

where $W[x]$ is the Lambert W function. Since no analytical solution exists for the steady-state solution under the uniaxial extensional flow, we numerically solve Eq. (20) to obtain steady extensional viscosity.

Finally, in this study, we do not use the popular VCM model to fit the slip-spring data because the VCM model is a two species dumbbell model, whose linear viscoelastic behavior is expressed by the superposition of two relaxation times. Thus, it does not give a good fit to the linear viscoelastic data of the slip-spring model. In addition, the VCM model is based on dilute dumbbells with no entanglement effects, making it physically unsuitable to describe the shear thinning produced by entanglements in entangled WLM solutions.

We first show the linear viscoelastic (LVE) properties of the micelle-slip-spring model with the Fraenkel spring potential. Before calculating the LVE data, all WLM chains in the ensemble are equilibrated without flow for time $teq$⁠. The relaxation time needed to equilibrate the micelles $teq$ depends on the dimensionless breakage time $ζ$ as $teq≥2(τ¯repτ¯breakeq)1/2$ for $ζ<1$ or $teq≥2τ¯rep$ for $ζ≥1$⁠. Here, $(τ¯repτ¯breakeq)1/2$ is the theoretically estimated longest relaxation time in the fast-breakage limit by Cates [25].

In our previous paper [36], we presented $G(t)$ data from the slip-spring model using Hookean springs obtained by relaxation after a small step shear strain. We here recompute $G(t)$ for the Fraenkel springs from equilibrium simulations using the Green–Kubo formula. We carry out this recalculation for the following two reasons. First, $G(t)$ of Fraenkel springs is slightly different from that of Hookean springs with the same numbers of main chain springs and slip-springs. (For more detail, see Fig. S3 in the supplementary material [45].) Since $G(t)$ is a small-strain property, predictions of it from the slip-spring model with Fraenkel springs should converge to those with Hookean springs when the Fraenkel spring constant is very large and there are enough short Fraenkel springs to recover Hookean behavior. To simulate with a large Fraenkel spring constant, a very small $Δt$ is required, and thus, it takes a long simulation time to obtain reliable results. Even with this, since the nonlinearity arising from the Fraenkel potential exists even in the small strain region, it is difficult to obtain an accurate $G(t)$ from relaxation after a small step shear strain, as also reported by Lin and Das who compared the relaxation behavior of Hookean and Fraenkel springs [50]. For more detail, see Sec. S4 in the supplementary material [45]. Thus, we use the Green–Kubo formula to obtain accurate $G(t)$ data. To compute $G(t)$⁠, we utilize Eq. (6) and the multiple-tau correlator for efficient calculations [59].

Figure 2 shows the relaxation modulus $G(t)$ for (a) $⟨Ztube⟩=5$ and (b) $7$ at four dimensionless breakage times $ζ=0.01$ (black), $0.1$ (red), $1$ (blue), and $10$ (magenta). Here, $G0$ is the modulus defined as $G0=⟨N⟩Gchain$ utilized as the unit of stress in Likhtman’s work [37]. After the equilibrium simulation for the time period $teq$⁠, we performed the simulation with $κ=0$ for a time duration $tLVE=(ttotal−teq)=1×105τ0$ to compute $G(t)$⁠, where $ttotal$ is the total simulation time. This $tLVE$ is at least 50 times longer than the longest relaxation time (cf. Table S1 and S2 in the supplementary material [45]). It should be noted that $G(t)$ for $⟨Ztube⟩=7$ and $ζ=10$ shown in Fig. 2(b) for $t≳3000τ0$ seems not to be converged due to the limited computational time. Thus, only the $G(t)$ data for $t≲3000τ0$ is utilized in the later analysis. From Fig. 2, the short-time relaxation behavior wherein $G(t)$ exceeds unity for $t≲10−2τ0$ is absent for Hookean chains, for which $G(t)$ normalized as described above approaches unity at short times. The rise of $G(t)$ above unity at short time is due to the Fraenkel potential used in this study, as noted in the literature [50]. Moreover, relaxation behavior in the relatively short time region $t≲τ0$ is almost identical for all breakage rates. This is because the breakage time is much larger than the unit time of the system. The difference between different breakage times can be seen in the long-time region. As expected, the relaxation becomes slower with increasing breakage time (i.e., decreasing breakage rate). The $G(t)$ results obtained from slip-spring simulations with the Fraenkel potential are then fitted by the multimode Maxwell model $G(t)=∑pGpexp⁡(−t/τp)$⁠, where $τp$ and $Gp$⁠, respectively, are the relaxation time and strength for mode $p$⁠. The fitting results are shown in Fig. 2 with solid lines. We will use the set of parameters ${τp,Gp}$ determined by these fits to $G(t)$ in the following analysis of nonlinear rheology under the flow. We tabulate the values of ${τp,Gp}$ in Sec. S5 of the supplementary material [45].

Next, we investigate nonlinear shear rheological properties. In all of the slip-spring simulations in this subsection, which focuses on shear properties, the breakage and rejoining times are assumed to retain constant values (i.e., $τ¯break=τ¯rejoin=τ¯breakeq$⁠) so that the breakage rate does not depend on configuration or stress. Extensional flow experiments for the CTAB/NaSal system [13] suggest that a micelle will break when under a stress of around $3.6×105$ Pa, roughly consistent with molecular dynamics simulations [47]. This scission stress is high enough that micelle scission is likely to occur only under the strong extensional flow and not under the shear flow. The effect of stress-induced micelle breakage at the high-stress level under the shear flow will be addressed in our future work, if necessary.

Figure 3 shows the steady shear viscosities $η≡σxyR/γ˙$ for three dimensionless breakage times (a) $ζ=0.01$⁠, (b) $0.1$⁠, and (c) $1$⁠. The slip-spring simulation results with $⟨N⟩=35$ springs and $⟨Ztube⟩=5$ entanglements are plotted with symbols. As can be seen in Fig. 3, the micelle-slip-spring model can reproduce shear thinning due to the effect of the entanglement constraint. Here, it should be noted that the shear-rate dependencies of the shear viscosity of WLM and polymer solutions are characterized by a power law, i.e., $η∝γ˙n$⁠, but with different exponents. The exponent $n$ for WLM solutions is typically around $n≃−1$ [1], whereas for nearly monodisperse polymer solutions, it is typically around $n≃−0.8$ [60]. This is an important difference since the power law of $η∝γ˙−1$ in WLM solutions, which corresponds to a stress plateau in the steady-state flow curve, is responsible for shear banding. In Fig. 3, however, the number of entanglements and the shear rate ranges are too small to observe this characteristic exponent for WLMs (⁠$n≃−1$⁠). To investigate the higher shear rate region, the missing mechanisms here (e.g., the entanglement loss under flow) should be incorporated into the micelle-slip-spring model. Although the limiting case in which there is flow-induced loss of entanglements but no regeneration of them is examined in Sec. III D, consideration of a better slip-spring regeneration algorithm under flow is deferred to a future work. Red solid and blue dotted lines in Fig. 3 are the fitting results for $γ˙τ¯R≲1$⁠, where $τ¯R$ is the average Rouse relaxation time, by the Giesekus model [cf. Eqs. (14)$∼$(16)] and the PTT model [cf. Eq. (22)], respectively. The nonlinear parameter appearing in each constitutive equation is the same for all modes. Thus, a single nonlinear parameter is used in each constitutive model to fit the slip-spring data. We observe that the steady shear viscosities for $γ˙τ¯R≲1$ obtained by the slip-spring model can be successfully reproduced by both the Giesekus and the PTT models with nonlinear parameter values summarized in Table I. On the other hand, the slip-spring simulation data for $γ˙τ¯R>1$ are not well reproduced by the constitutive models. Here, we note that there is no physical reason that a near perfect agreement should be obtained since the constitutive models used here are phenomenological and are not microstructurally based. Almost the same trend can be observed for the case of $⟨N⟩=49$ springs and $⟨Ztube⟩=7$ entanglements, as shown in Fig. S5(a) in the supplementary material [45].

Figure 4 shows transient shear viscosities $η+(t)≡σxyR(t)/γ˙$ for (a) $ζ=0.01$⁠, (b) $0.1$⁠, and (c) $1$⁠. Slip-spring simulation results with $⟨N⟩=35$ (i.e., $⟨Ztube⟩=5$⁠) are plotted with symbols. Here, the same color in each figure means the same strain rate. The linear viscosity growth function $η0+(t)$ computed from the LVE data is shown with the black bold dotted line. The stress overshoot shown by the slip-spring model predictions for $γ˙τ¯R≲1$ qualitatively matches that seen in experiments [5,6]. For the high strain rate region (⁠$γ˙τ¯R>1$⁠), the shear viscosity exceeds the linear viscosity growth function, thus showing a “hardening” prior to the stress overshoot. This may be related to the nonlinearity of the chain extension at high strain rates because of the small number of Kuhn steps per spring in our model. This hardening behavior under strong shear, which is attributed to the finite extensibility of springs, is observed in experiments for CTAB/NaSal systems, as reported by Inoue and co-workers [12]. Similarly, stiff entangled polymers are known to produce strain hardening in shear [61]. This strain hardening in our model will likely weaken if strain-induced entanglement loss at high shear rates is included. As can be seen in the solid lines of Fig. 4, the transient shear viscosities, including the stress overshoots seen for $γ˙τ¯R≲1$⁠, are in good agreement with those predicted by the Giesekus model with the nonlinear parameter determined by the steady state data in Fig. 3. While the steady-state slip-spring simulation data for $γ˙τ¯R≲1$ are reasonably captured by the PTT model as shown in Fig. 3, the transient shear viscosity predictions are poorer than those of the Giesekus model in that the stress overshoot in the PTT model is less prominent than in the slip-spring model. The hardening behavior seen in $γ˙τ¯R>1$ is not reproduced by these constitutive models, where the effect of finite extensibility is not addressed. Almost the same trend can be observed for the case of $⟨N⟩=49$ springs and $⟨Ztube⟩=7$ entanglements, as shown in Fig. S5(b) in the supplementary material [45]. Better predictions by the PTT model might be obtained by setting the slip-parameter in the Gordon–Schowalter convected derivative to $a≠1$⁠. However, doing so produces unphysical stress oscillations under the shear flow, which are not seen in current slip-spring simulations.

Figure 5 shows normal stress coefficients $Ψ1+(t)$ $≡N1(t)/γ˙2={σxxR(t)−σyyR(t)}/γ˙2$ for (a) $ζ=0.01$⁠, (b) $0.1$⁠, and (c) $1$ with the same parameters as those in Fig. 4. The first normal stress growth function $Ψ1,0+(t)$ [54] computed from the LVE data is also shown with the black bold dotted line. The Giesekus model (solid lines) can reproduce both the transient and the steady state data predicted by the slip-spring model for $γ˙τ¯R≲1$⁠. However, as is the case for the transient shear viscosity, the PTT model predictions (dotted lines) for both the transient and steady normal stress coefficients deviate from the slip-spring simulation data for $γ˙τ¯R≲1$⁠. Moreover, both constitutive models cannot reproduce $Ψ1+(t)$ for $γ˙τ¯R>1$⁠. Almost the same trend can be observed for the case of $⟨N⟩=49$ springs and $⟨Ztube⟩=7$ entanglements, as shown in Fig. S5(c) in the supplementary material [45].

From the comparison of slip-spring simulation data with the Giesekus and PTT models in Figs. 4 and 5, the Giesekus model matches better the slip-spring data than does the PTT model under a moderate shear flow where the chains do not show significant stretch. Since experimental studies show that the Giesekus model can also reproduce the shear rheological data of WLM solutions [8,9], the micelle-slip-spring model will also be in qualitative agreement with the experimental shear data. We note that the strain hardening behavior seen in both the slip-spring model and the experiments for the CTAB/NaSal system [12] is not reproduced by the Giesekus model. Thus, we can conclude that a constitutive model that can predict rheological properties under a strong shear flow has not yet been developed.

In experiments of shear-rate regions showing shear thinning, flows of typical WLM solutions show shear banding [22,23]. To predict shear banding, spatially dependent simulations are required. Such simulations with the micelle-slip-spring model incur excessive computational costs. Since the slip-spring model predictions for $γ˙τ¯R≲1$ are well reproduced by the Giesekus model, the latter can be used for such spatially dependent simulations. This will be the goal of our future research whose aim is to use slip-spring simulations to choose and tune a constitutive model to produce similar responses in flows such as shear, extensional, and mixed flows, thus validating it for use in more complex flows.

Next, we investigate uniaxial extensional rheological properties obtained by the micelle-slip-spring model again with Fraenkel springs. As for the shear flow, we first omit stress-induced micelle breakage (i.e., $τ¯break=τ¯breakeq$⁠) and show the basic properties of the micelle-slip-spring model under the extensional flow. However, as discussed above, experiments [13] and recent molecular dynamics simulations [47] reveal that stress-induced micelle breakage (SIMB) is expected when the stress becomes several orders of magnitude larger than the plateau modulus. This might be the case at the high extensional rates examined in this study. Thus, we will also test the effect of SIMB for several conditions.

Figure 6 shows steady extensional viscosities $ηE=σE/ϵ˙$ normalized by the zero shear viscosity $η0$ for three dimensionless breakage times (a) $ζ=0.01$⁠, (b) $0.1$⁠, and (c) $1$ with the steady extensional stress $σE$ defined by $σE=σxxR−(σyyR+σzzR)/2$⁠. Slip-spring simulation results without SIMB are plotted with open symbols. To compare simulation results with the experimental data, it is convenient to introduce the Weissenberg number defined as $Wiext≡ϵ˙⟨τ⟩$⁠, where $⟨τ⟩$ is an average relaxation time defined as $⟨τ⟩=(∑p=1NmodeGpτp2)/(∑p=1NmodeGpτp)$⁠. In Fig. 6, the simulation results showing significant spring stretch beyond the theoretical maximum value are plotted with open squares. Under high extensional flows, the average squared end-to-end vector for the micellar chains, $⟨REE2⟩$⁠, becomes larger than the theoretical maximum value $⟨REE,max2⟩=(⟨N⟩b)2$ since longer chains in the ensemble show significant stretch (for more details, see Fig. S6 in the supplementary material [45]). This is evident for slower breakage cases where longer chains less frequently break. There are two possible reasons for this overstretch: one is the Fraenkel spring constant that is not large enough; the other is $Δt$ values (⁠$Δt∼10−5τ0$ for high extensional rates) that are not small enough. However, it is not practical to employ larger Fraenkel parameters or smaller $Δt$ values due to the large computational cost. In the current case, however, since the number of these chains is small, the extensional stress is expected not to be much affected by these chains. This is partly supported by the fact that the extensional stress does not show significant divergence (cf. Fig. 7). In Fig. 6, we also plot with gray solid and dotted lines the predictions of the Giesekus model and the PTT model, respectively, with the same nonlinear parameter as that under the shear flow. From Fig. 6, we can observe that extensional viscosities from slip-spring simulations without SIMB monotonically increase at $Wiext≳0.5$⁠. This is a typical behavior obtained by the bead-spring models [51]. On the other hand, extensional viscosities predicted by the Giesekus and PTT models show extension thickening followed by a high Weissenberg number plateau and extension thinning, respectively. Experimentally, Walker and co-workers reported that the extensional viscosities for CPyCl/NaSal solutions show extension thickening followed by extension thinning, where the increase in $ηE$ occurs at $Wiext(expt)≃0.7$ and the maximum in $ηE$ is seen at around $Wiext(expt)$ $≃4∼7$ [17]. Here, $Wiext(expt)$ is obtained by multiplying experimental apparent strain rates and a terminal relaxation time $τexpt$⁠, which is evaluated by $τexpt=η0/GN0$⁠. While the slip-spring simulations can reasonably capture the onset of the increase in $ηE$⁠, the original slip-spring simulations without SIMB cannot reproduce extension thinning seen in experiments. Thus, there should be a missing mechanism(s) in the original micelle-slip-spring model under a strong extensional flow. It should be noted that the PTT model can qualitatively capture experimental extensional viscosities, i.e., extension thickening followed by extension thinning. However, the maximum in $ηE$ becomes less prominent with increasing dimensionless breakage time.

As discussed earlier, we incorporate stress-induced micelle breakage (SIMB) to improve slip-spring predictions under a strong extensional flow. To implement SIMB, the characteristic extensional stress, $σEbreak$⁠, appearing in Eq. (12) should be determined. This can be done with the assistance of the molecular dynamics simulation [47]. For example, for the cetyltrimethylammonium chloride (CTAC)/NaSal system, utilizing the scission energy vs stress data shown in Fig. 6(b) of the work by Mandal and Larson [47] and performing a linear fit [cf. Eq. (11)], we can estimate $σEbreak$ to be $σEbreak≃3×105$ Pa for the CTAC/NaSal system at the salt concentration giving the lowest breakage rate. This value is much higher than the experimental plateau modulus of CTAC/NaSal systems [62]. It should be noted that the characteristic stress might be system dependent and has not been systematically examined. To the best of our knowledge, for the CPyCl/NaSal system, whose extensional properties are reported by Walker and co-workers [17], no information to estimate the characteristic stress has been reported. Here, based on the finding that $σEbreak$ is much higher than the plateau modulus, we assume a value much lower than that estimated by Mandal and Larson, that is we take $σEbreak=2×104$ Pa, which allows us to find a significant effect of SIMB on rheological properties. (Comparison to the larger value of $σEbreak$⁠, i.e., $σEbreak=1×105$ Pa closer to the value of Mandal and Larson [47], giving much less extension thinning, is shown in Sec. S8 in the supplementary material [45].) If we assume the experimental plateau modulus as $GN0≃10$ Pa, the dimensionless characteristic extensional stress is evaluated as $σ~Ebreak≃2×103$⁠. In the slip-spring model, the plateau modulus can be approximated as $GN,ss0=0.18×0.8G0=0.144G0$ [49]. During simulations, we rescale the extensional stress by $GN,ss0$ to compute $σ~E$ and evaluate the breakage time by Eq. (12).

In Fig. 6, slip-spring simulation results with SIMB are plotted with filled symbols. We can observe that the extensional viscosities show the maximum at around $Wiext$ $≃2∼3$⁠, which is almost the same as in the experimental results by Walker and co-workers [17]. Here, $⟨REE2⟩$ with SIMB is smaller than or at most slightly larger than the theoretical maximum value. As compared to the simulations without SIMB, the degree of overstretch is, thus, significantly suppressed (for more details, see Fig. S6 in the supplementary material [45]). The constitutive models using the nonlinear parameter derived from fitting the shear flow (gray lines in Fig. 6) underestimate $ηE$ obtained by the slip-spring model with SIMB (filled symbols in Fig. 6). To improve the predictions of the constitutive models, we separately fit the Giesekus and PTT models to the extensional data with SIMB of $Wiext≳1$⁠. While we use Eqs. (17)–(19) for the Giesekus model fit, an optimal PTT fit to the extensional data can be obtained by minimizing the value of $∑j{(ηE,jPTT−ηE,j)/ηE,j}2$⁠, where $j$ is the individual data points with different strain rates, $ηE,jPTT$ is the steady extensional viscosity of the PTT model, and $ηE,j$ is the steady extensional viscosity of the slip-spring model. The parameter values under the extensional flow are summarized in Table I. The fitting results are shown in Fig. 6 with red lines. As can be seen in Table I, the Giesekus and PTT models with smaller nonlinear parameter values than under shear flow give better predictions under a uniaxial extensional flow. We can see that the Giesekus model can reproduce the slip-spring simulation results better than the PTT model. However, the Giesekus model cannot predict extension thinning observed in experiments and slip-spring simulations with SIMB. On the other hand, the PTT model can reproduce extension thinning but provides poorer predictions at the onset of extension thickening. It should be noted that the Giesekus and PTT models with the nonlinear parameter determined under the extensional flow give poorer predictions under the shear flow than those shown in Figs. 4 and 5 (for more details, see Fig. S8 in the supplementary material [45]).

Figure 7 shows transient uniaxial extensional viscosities $ηE+(t)$ for (a) $ζ=0.01$⁠, (b) $0.1$⁠, and (c) $1$⁠. For each $ζ$ value, $ηE+(t)$ curves of the slip-spring model with the smallest strain rate (⁠$ϵ˙τ0=1×10−2$ for $ζ=0.01$⁠, $ϵ˙τ0=3×10−3$ for $ζ=0.1$⁠, and $ϵ˙τ0=1×10−3$ for $ζ=1$⁠) are almost identical to $ηE0+$⁠, which indicates that these strain rates are in the linear response regions. In the small strain rate region (⁠$ϵ˙τ0≤3×10−2$ for $ζ=0.01$⁠, and $ϵ˙τ0≤1×10−2$ for $ζ=0.1$ and $ζ=1$⁠), $ηE+(t)$ curves with and without SIMB are almost identical, which indicates that SIMB does not affect $ηE+(t)$ in the small strain rate region. In the larger strain rate region, the effect of SIMB can be detected for all $ζ$ cases. While transient extension hardening behavior with SIMB is almost identical to that without SIMB, $ηE+(t)$ with SIMB shows a maximum before reaching a steady state, which is the qualitatively different behavior from that without SIMB. As shown in Fig. 6 with filled symbols, the steady extensional viscosities show weak extension thinning.

In Fig. 8, we plot the Giesekus and PTT model predictions with modified nonlinear parameters for the extensional flow. In the small strain rate region, the Giesekus model predictions (solid lines) can be brought into reasonable agreement with transient extensional data from the slip-spring model. However, deviations from the slip-spring simulation results can be observed for the larger strain rate region. In particular, the viscosity overshoot observed in the slip-spring simulations is not reproduced by the Giesekus model. Almost the same arguments are true for the PTT model predictions (dotted lines). Thus, the nonlinear shear and extensional slip-spring data are not simultaneously reproduced by the Giesekus and PTT models with a single nonlinear parameter.

Thus, our attempt to fit slip-spring model predictions with those of common constitutive equations under both shear and extensional flows has not yet succeeded, although the Giesekus model can give good predictions under the moderate shear flow. A complete constitutive equation might, however, be developed by replacing the nonlinear term of the Giesekus equation with a functional form sensitive to the flow type or by changing the nonlinear parameter $αG$ to a flow-type parameter that can distinguish shear from the extensional flow. The latter idea would be tricky and would require choosing a suitable metric of flow type that satisfies the proper invariance principles. If such a flow-type parameter is found that allows the model to predict results in both shear and extension, it would need to be checked by comparing its predictions to those of slip-spring simulations in mixed flows. These mixed flows are intermediate between shear and extension, for which the flow switches between shear to extension as a function of time. Such work is beyond the scope of this study and, therefore, deferred to our future work.

Here, we examine the effect of the entanglement loss, which is one of the important nonlinear mechanisms and is not considered so far, on the rheological properties. We implement this effect by considering only a limiting case where the total number of slip-springs (⁠$Ztotal$⁠) continually decreases and eventually becomes zero. Namely, when a slip-link goes through its chain end, this slip-link and the partner slip-link are simultaneously destroyed, and no regeneration of slip-links is considered. This method is similar to that employed by Moghadam and co-workers [63]. From Figs. 9(a) and 9(b), the shear and extensional viscosities without slip-spring regeneration are smaller than those with slip-spring regeneration. This stress reduction corresponds to the decrease in the number of slip-springs seen in the insets of Figs. 9(a) and 9(b) and should contribute to shear or extension thinning behavior. As can be seen in Fig. 9(a), the strain hardening is hardly observed for the case without slip-spring regeneration. This is because the effect of slip-spring constraint is quickly lost due to fast breakage. Nevertheless, the hardening behavior is expected to reappear since the unentangled bead-spring model with the Fraenkel springs itself shows the hardening under fast shear, as shown in Fig. S9 in the supplementary material [45].

Since the extent of entanglement loss in the nonlinear flow is still under investigation, our limiting cases of constant entanglement density and no regeneration of entanglements provide useful bounds on the expected behavior. A correct algorithm is somewhere between two cases since the entanglement regeneration process should exist even under the flow. The correct entanglement regeneration algorithm should be developed first for polymers without breakage and rejoining events.

Finally, we present a comparison between the slip-spring results and experimental data for an aqueous solution of 100 mM cetylpyridinium chloride (CPyCl) with 45 mM sodium salicylate (NaSal) examined by Gaudino and co-workers [7]. It should be noted that the experimental results for 100 mM CPyCl with 45 mM NaSal have a rather broad relaxation time distribution, as can be seen in Fig. 10. However, it is typical that entangled linear WLM solutions with a sufficiently high surfactant or salt concentration show a Maxwell-type stress relaxation with a single relaxation time [2–4]. This relaxation behavior is observed if the number of entanglements is large enough that the reptation dynamics is well separated from Rouse modes, and the breakage rate is fast enough that the relaxation process can be described by a single relaxation time. In this study, the computational cost limits the model to WLM solutions with moderate numbers of entanglements (⁠$≲7$⁠), which is not enough to show the Maxwell-type stress relaxation with a single relaxation time. Nevertheless, in the future, we hope to perform nonlinear simulations with a larger number of entanglements by developing a further computationally efficient code or a constitutive model for well-entangled WLMs based on our micelle-slip-spring model.

To test the ability of the slip-spring model to predict the nonlinear experimental data, the slip-spring parameters (i.e., the number of entanglements $⟨Ztube⟩$ and the dimensionless breakage time $ζ$⁠) should be determined before performing linear and nonlinear rheological simulations. We can estimate $⟨Ztube⟩$ from the linear rheological data with the aid of the following relationship developed from the pointer algorithm by Tan and co-workers [30]:

Here, $Gmin″$ is the local minimum in $G″$⁠, $Gmin′$ is the $G′$ value at the frequency at which $G″$ shows its minimum value $Gmin″$⁠, $⟨L⟩$ is the average micelle length, and $ℓe$ is the entanglement length. However, the target experimental data for the CPyCl/NaSal system does not have a clear minimum in $G″$ [see Fig. 1(a) of [7]]. In this case, it is reasonable to assume that the $Gmin′/Gmin″$ value is a small number, $1−2$⁠, which gives $⟨Ztube⟩≃4−9$⁠. In this study, we, therefore, assume that $Gmin′/Gmin″$ is $1.5$⁠, which gives $⟨Ztube⟩≃7$⁠. (The experimental data showing $Gmin′/Gmin″=1.5$ are at about $ω≃20$$s−1$⁠.) This $⟨Ztube⟩$ value is larger than the value $⟨Ztube⟩=2.5$ that Gaudino and co-workers obtained from an earlier scaling formula lacking a prefactor, namely, $⟨Ztube⟩=GN0/Gmin″$⁠, with $GN0$ being a rough estimate of the plateau modulus since the experimental data do not show a clear plateau from which to obtain $GN0$ [7]. Thus, the value of $⟨Ztube⟩=2.5$ is expected to be smaller than the more realistic value estimated from Eq. (23). For the dimensionless breakage time $ζ$⁠, we use a relatively large value $ζ=10$ since the WLM solution is expected to be well outside the fast breakage limit. This can be inferred from the fact that four modes are required to fit the LVE data (cf. Table I of [7]).

In addition to $⟨Ztube⟩$ and $ζ$⁠, we need to estimate the theoretical values of the unit stress $G0(theory)$ and the unit time $τ0(theory)$ of the slip-spring model. Before computing these values, we estimate the plateau modulus of this WLM solution, $GN0$⁠. We use the following relation established by Tan and co-workers [30]: $GN0/Gmin′=4.25/(Gmin′/Gmin″)+0.625$⁠. This equation is recommended for $Gmin′/Gmin″<10$⁠, which is the case in our study. From our assumption of $Gmin′/Gmin″=1.5$ and $Gmin′=G′(ω≃20$ $s−1)≃10$ Pa, $GN0$ is determined as $GN0≃35$ Pa. Using the same procedure shown in our previous study [36], we can determine $G0(theory)$ as $G0(theory)=GN0/(0.8×0.18)≃2×102$ Pa. In the loosely entangled regime, $GN0$ is expressed as $GN0=9.75kBT/ℓblob3$ where $ℓblob$ is the blob size [27]. Using this equation and setting $T=300$ K, we can compute $ℓblob$ as $1.1×10−7$ m. Taking for the persistence length $ℓp$ the typical value $ℓp=25$ nm and using the relation $ℓblob=ℓe0.6ℓp0.4$⁠, we can determine the entanglement length $ℓe$ as $ℓe≃280$ nm. Thus, the semiflexibility factor $αe$ is $αe=ℓe/ℓp≃11$⁠, which is close to the standard value used in the slip-spring model, $αess=14$⁠. In the slip-spring model, the unit of time $τ0$ is $τ0=ξb2/kBT=ξ′b3/kBT$⁠, with $ξ′$ being the drag coefficient per unit micelle length. The value of $ξ′$ can be computed using the drag coefficient for a cylinder of length $ℓblob$ and micelle diameter $d$ as $ξ′=2πηs/ln⁡(ℓblob/d)$⁠. Here, $ηs$ is the solvent viscosity. From $T=300$ K, $ηs=0.85$ $mPa⋅s$⁠, and the typical value of $d=4$ nm, the theoretical value of the unit time can be computed as $τ0(theory)≃5×10−5$ s.

Figure 10 compares the storage and loss moduli for the slip-spring simulation with $⟨Ztube⟩=7$ and $ζ=10$ (symbols) with the experimental data (lines). As noted in Sec. III A, we can use the slip-spring data for $t/τ0≲3000$ to accurately evaluate $G′$ and $G″$ (cf. Fig. 2). Here, rather than plot the original $G′$ and $G″$ experimental data, we smooth these data using the fits to them by a set of Maxwell modes whose relaxation times and strengths are given in Table I of [7]. The unit time and unit stress determined by shifting the dimensionless $G′$ and $G″$ curves to match the experimental $G′$ and $G″$ are $τ0=5×10−5s$ and $G0=2×102Pa$⁠, respectively. These values are the same as the theoretically determined unit time and stress since the semiflexibility factor of the slip-spring model $αess=14$ is close to the experimental value of $αe≃11$⁠. Unlike our previous study [36], we need not correct $τ0$ and $G0$ to address the different $αe$ values between the slip-spring simulation and the experiment. From Fig. 10, we can see that $G′$ and $G″$ obtained from slip-spring simulation are in reasonable agreement with the experimental data. This indicates that the choices of slip-spring parameters, i.e., $⟨Ztube⟩$ and $ζ$⁠, are reasonable.

Figure 11 compares $η+(t)$ predicted by the slip-spring simulations without SIMB and entanglement loss for $⟨Ztube⟩=7$ and $ζ=10$ (lines) with the experimental data [7] (symbols). Here, we also use the values of $τ0$ and $G0$ determined in Fig. 10. The dimensionless shear rates $γ˙τ0$ in the slip-spring model that correspond to the experimental values of $20$⁠, $40$⁠, and $80s−1$ are, therefore, about $γ˙τ0≃1×10−3$⁠, $2×10−3$⁠, and $4×10−3$⁠, respectively. The slip-spring results in Fig. 11 are shown with solid lines, and $η0+(t)$ computed from the LVE data is presented with the black dotted line. Figure 11 indicates that the slip-spring data are in good agreement with the experimental shear data. Again, this shows that the choices of slip-spring parameters are reasonable.

To the best of our knowledge, experimental data in the extensional flow for WLM solutions with micelles of length accessible by the slip-spring model are not available currently. Extensional measurements for WLM solutions with low viscosities are typically conducted by the Capillary Breakup Extensional Rheometry (CaBER). However, the strain rate of the CaBER measurements is not constant over time, making it difficult to compare the extensional experimental data with the slip-spring data. Additionally, an initial step extensional strain in the CaBER measurements might affect the extensional viscosity [14]. It is difficult to produce the initial state of micellar chains for the slip-spring simulations. The comparison with extensional experiments is an important problem to study in the future.

We have presented linear and nonlinear rheological predictions of a mesoscopic simulation model for entangled WLM solutions, the so-called “micelle-slip-spring model” based on the polymer slip-spring model proposed by Likhtman [37] with breakage and rejoining events added to reproduce WLM dynamics [36]. To capture the nonlinear rheological properties, instead of the Hookean springs employed in our previous study [36], we introduced the Fraenkel potential for main-chain springs and slip-springs. Additionally, to improve extensional properties under the strong extensional flow, we incorporated stress-induced micelle breakage into the micelle-slip-spring model. We used this extended model to examine both linear and nonlinear rheological properties.

In the linear response region, the relaxation modulus $G(t)$ was computed using the Green–Kubo formula from equilibrium simulation. We found that $G(t)$ for Frankel springs differs slightly from that for Hookean springs due to the difference in spring potential. The results for $G(t)$ are fitted using the multimode Maxwell model to determine relaxation times and strengths.

This micelle-slip-spring model with Fraenkel springs can reproduce the characteristic rheological behavior in start-up and steady-state shear of WLM solutions. That is, the transient shear viscosities for moderate strain rates show the experimentally observed overshoot and the steady shear viscosities show the observed shear thinning. The experimentally observed hardening behavior under the strong shear flow can be reproduced by the micelle-slip-spring model with finite extensibility. The slip-spring simulation results for the strain rates where the chains do not show significant chain stretch were fitted by two common phenomenological constitutive equations: the Giesekus model and the PTT model. The steady shear viscosity obtained by the slip-spring model can be well reproduced by both constitutive models. Transient shear predictions by the Giesekus model are also in good agreement with those of the slip-spring model. This result is qualitatively the same as in experimental observations. In contrast, transient predictions by the PTT model (with the non-affine deformation parameter set to zero to avoid spurious oscillations in shear) deviate from those by the slip-spring model. These results indicate that the Giesekus model captures better the shear properties of the slip-spring model for not so high shear rates than does the PTT model.

The micelle-slip-spring model was also applied to predict rheological properties under the uniaxial extensional flow. The Fraenkel potential prevents the overstretch of WLM chains under the extensional flow, and thus, we can reach steady-state extensional viscosities. While the micelle-slip-spring model reasonably reproduces the onset of extension thickening behavior, this model fails to predict extension thinning behavior observed in experiments. Thus, we incorporated stress-induced micelle breakage (SIMB), which allows the simulation results with SIMB to successfully capture extension thinning behavior. We also show the Giesekus and PTT constitutive model predictions with the same nonlinear parameter as under shear flow underestimate the slip-spring results. To reproduce the steady extensional data of the slip-spring model, a smaller value of the nonlinear parameter (corresponding to smaller deviations from the upper-convected Maxwell model) is needed to fit the slip-spring results in the extensional flow. However, the constitutive model predictions even with the modified nonlinear parameter cannot capture the transient extensional viscosities especially under the strong flow. Thus, it is not possible to predict both shear and extensional properties by the Giesekus and PTT models with a single constant nonlinear parameter.

Furthermore, we examined the effect of entanglement loss on rheological properties by considering the limiting case in which we disable the regeneration of slip-springs when previous ones disappear, so that the number of slip-springs continuously decreases and eventually becomes zero. We found that the stress without the slip-spring regeneration is smaller than that with slip-spring regeneration. Thus, the effect of entanglement loss should contribute to shear/extension thinning behavior. A correct algorithm for entanglement loss is somewhere between two cases of constant entanglement density and no regeneration of slip-springs and should be developed first for unbreakable polymers.

The slip-spring simulation results were also compared with experimental linear rheological and nonlinear shear data. The number of entanglements $⟨Ztube⟩$ was estimated from the linear rheological data using the recent relation found by Tan and co-workers [30], and the dimensionless breakage time $ζ$ was set to be large based on the inference that the WLM solution is well outside the fast breakage limit. The simulation results for $G′$ and $G″$ were compared with the experimental results, and the unit time $τ0$ and the stress $G0$ were determined. These $τ0$ and $G0$ are the same as expected theoretical values, indicating the validity of the chosen parameters. Moreover, nonlinear shear viscosities obtained by the slip-spring model are in reasonable agreement with the experimental shear data. These results support the predictive ability of the slip-spring model under the homogeneous shear flow.

Further studies are required to understand more deeply the rheological properties of WLM solutions using the micelle-slip-spring model. For example, strain hardening in nonlinear stress relaxation is experimentally detected in several WLM systems [64,65], which has been linked to shear-induced associations of WLMs. Hence, to reproduce strain hardening in nonlinear stress relaxation, it might be important to incorporate this mechanism into the micelle-slip-spring model. Moreover, it would be interesting to examine the change in strain-rate dependence of shear stress, which is related to shear banding, caused by stress-induced acceleration of micelle breakage. Additionally, the extensional rheology obtained by the slip-spring model should be compared with experimental results. The current simulations are limited to the homogeneous flow. To incorporate inhomogeneous shear such as shear banding, the equations of conservation of mass and momentum should be coupled. While “multiscale” slip-link simulations have been performed for unbreakable polymer melts [66], such simulations are currently too extensive to be carried out for WLM solutions since an ensemble of many WLM chains is required. Thus, the fitting of predictions of the slip-spring model by a phenomenological constitutive equation, which is attempted for the first time in this study, might be effective for such a purpose. This might be done more efficiently in the future with the aid of machine-learning techniques [67]. We hope to continue our study along these directions.

T.S. would like to express gratitude to Professor T. Taniguchi for his comments and fruitful discussion on this work. This work was partially supported by JSPS KAKENHI under Grant No. 21K13893. R.G.L. was funded by 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.

The authors have no conflicts to disclose.

Ethics approval is not required.