A three-species model for wormlike micellar fluids

Computers and Mathematics with Applications (

)

–

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

A three-species model for wormlike micellar fluids Hunseok Kang a , Young-Ju Lee b,∗ a

Research Center for Nonlinear Ergodic Theory, Sungkyunkwan University, Suwon 440-746, Republic of Korea

b

Department of Mathematics, Texas State University, San Marcos, TX 78666, United States

article

info

Article history: Received 17 September 2015 Received in revised form 4 January 2016 Accepted 10 February 2016 Available online xxxx Keywords: Wormlike micellar fluids Three-species model Global stability Gelation Shear-thickening transition

abstract We present a new population dynamics-based three-species model for modeling wormlike micellar fluids which incorporates the gelation species. Starting at a homogeneous three-species model (Nestor, 2005), we show that the model is globally stable and approaches to the steady state. We then extend the model to the partial differential equation to incorporate the spatial dependence of the species. Their inhomogeneous effect for the flow behavior is modeled as well, by including populations as a polymeric viscosity in the fluid. Using the newly proposed inhomogeneous three-species model, we tackle the challenge to obtain the shear-thickening transition observed in experiments (Liu and Pine, 1996). Guided by Lerouge and Berret (2010), we design parameters in the model and numerically show that they lead to obtain an appropriate shear-thickening transition. The results are shown to be agreeable qualitatively well with experimental results, thereby confirming the conjecture by Lerouge and Berret (2010) affirmatively. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction A micelle consists of a number of amphiphilic molecules or surfactants and each surfactant consists of hydrophobic carbon tail and hydrophilic head group. A large concentration of amphiphiles in the water, can lead to a spherical and cylindrical micellization (or equivalently ‘‘self-organization’’). Such self-organization or the micellization occurs due to the nature of hydrophobicity of the carbon tail and hydrophilicity of the head group of the surfactants. Therefore, wormlike micelles are self-assemblies of amphiphilic surfactant molecules in an aqueous solvent [1,2]. Due to the nature of the self-organization, there can be much variations in size in the formation of micelles, thereby making the modeling of wormlike micellar fluids challenging. Another important aspect in wormlike micellar fluids is the shear-thickening transition. Such transition has been noted to occur above at a critical shear rate and known to be one of the most puzzling phenomena (see [1] and references cited therein). The shear-thickening effect is known to be first observed by Rehage and Hoffmann [3] in 1981. It is only very recent to realize that such a transition is due to a large structure known as so-called the ‘‘gelation’’ or shear induced structure SIS. The first and direct experimental observation for the gelation is made by Liu and Pine [4] using a novel light scattering microscopy technique. In experiments performed by Liu and Pine [4], Taylor–Couette flow is used, in which outer cylinder is rotated with some shear rate, while the inner cylinder is fixed. It is observed that gel bands initially accumulate at the inner boundary after shearing for a few minutes above at a critical shear rate, the gel bands then become unstable, and then finger like gel structures begin to grow from the inner stationary cylinder toward the outer boundary. The emergence of the

∗

Corresponding author. E-mail addresses: [email protected] (H. Kang), [email protected] (Y.-J. Lee).

http://dx.doi.org/10.1016/j.camwa.2016.02.014 0898-1221/© 2016 Elsevier Ltd. All rights reserved.

2

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

Fig. 1.1. Schematics of effective viscosity as a function of shear rate for CTAB/NaSal in water (Left) and the evolution of the shear stress at an applied shear rate exceeding a critical shear rate (Right) in [4].

shear-induced gel phase is observed to dramatically raise the stress in the solution and gives rise to shear-thickening transition [5]. Typical effects of the shear-thickening transition observed in the dilute wormlike micellar solutions are illustrated in Fig. 1.1. Note that the evolution of the effective viscosity is similar to that of the shear stress. There has been a lot of attempt to understand the mechanism of shear-thickening transitions in a shear flow for the last three decades. However, none of them are not fully successful [1]. In particular, earlier models have predicted too large shear rates above which the shear thickening occurs, in comparison with experimentally predicted critical shear rates [6,1]. Our aim in this paper is, therefore to mathematically model the generation and the break of the gelation that can reproduce the following physical phenomenon in a shear flow: 1: Shear banding and the generation of gelation are correlated. 2: Shear thickening is predicted above at a critical shear rate, Γc as seen in Fig. 1.1 (Left). 3: Evolution of effective viscosity as well as the stress agrees qualitatively well with the experiments as seen in Fig. 1.1 (Right). The critical shear rate, Γc will be built into our modeling, which is shown to achieve the shear-thickening transition exactly at the physically observed critical shear rate. In addition, to the best of authors’ knowledge, our work is the first attempt to model the gelation inhomogeneously observed in the wormlike micellar fluids, in the mathematical model. There are in fact two-species model in literatures (see [7] and references cited therein). However, none have taken into account the gelation species in the modelings. We remark that a number of interesting new physical experiments such as the oscillatory falling sphere [8] and jumping bubbles [9] in wormlike micellar fluids have been reported to be conjectured to be correlated with the generation and break of such gelation, SIS. Therefore, our initial attempt to incorporate the gelation due to the shearing in the mathematical model can be an important step toward the modeling of such nontrivial experiments. The rest of this paper is organized as follows. In Section 2, we investigate and generalize the homogeneous three-species model introduced in the work by Handzy [2] and provide stability analysis for the model. In Section 3, an inhomogeneous version of the three-species model is provided. In Section 4, we present the shear-thickening transition obtained by the proposed three-species model. We conclude our paper with several remarks and future directions in Section 5. 2. Homogeneous three-species model In this section, we review and generalize the homogeneous three-species model introduced by Handzy [2] for modeling wormlike micellar fluids. There is a large spectrum of micelle lengths in wormlike micellar fluids. Therefore, it is out of reach to take into account all the variations. In [2], a simple model has been proposed to take into account two sizes of micelles, described as short and long, together with a third species representing shear induced structure (SIS) or the ‘‘gelation’’. We denote the short micelles in concentration by S(t ), long micelles in concentration by L(t ), and ‘‘gelation’’ in concentration by G(t ), respectively. There are two main assumptions in the homogeneous three-species model proposed by Handzy [2]. (1) It is based on mass (or concentrations) exchange property of aforementioned species. Specifically, in the modeling, the long micelles are assumed to be p times as long as short micelles. Namely, p short micelles can join into one long micelle or one long micelle can break into p short micelles. On the other hand, it is assumed that q long micelles can be combined into a gelation, or p × q micelles can be combined to form a gelation. Reversely, the gelation can be broken into a combination of short and long micelles. Note that p and q can be any real numbers that are larger than one. In the treatment of the microscopic dynamics of wormlike micelles by chemical reaction kinetics, an important early study was made by [10]. (2) A scission and reforming of gelation and micelles are crucially dependent on the shear rate and therefore, when considering the three species as chemical reactants, reaction rates are given solely as a function of shear rate. Table in Fig. 2.1 describes the reaction and reaction rates between three species. Remark 2.1. Physically, one can view the small micelle as a spherical micelle and the long micelle as the cylindrical micelle [2,1,7]. In experiments, it is observed that the shear-thickening transition occurs only for surfactants that self-assemble

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

3

Fig. 2.1. Scheme of micelle reactions in the three-species model (Left) and Chemical-like reactions of the three species (Right).

into cylindrical micelles [1]. The reaction rates given in Table in Fig. 2.1 show that the gelation G can occur through L, not directly from S. Therefore, it is physically relevant. The chemical-like reactions of three species and their rates given in Fig. 2.1 lead to the mathematical model as a system of ordinary differential equations as follows: L˙ =

˙ = G ˙= S

dL dt dG dt dS dt

= −c1 (γ˙ )L − qc2 (γ˙ )Lq + c3 (γ˙ )Sp + α c4 (γ˙ )G

(2.1a)

= c2 (γ˙ )Lq − c4 (γ˙ )G

(2.1b)

= pc1 (γ˙ )L − pc3 (γ˙ )Sp + β c4 (γ˙ )G,

(2.1c)

where 1. ci = ci (γ˙ ) ≥ 0 for i = 1, 2, 3, 4; 2. α and β are positive numbers satisfying the equality pα + β = pq.

(2.2)

The three-species model (2.1) is conservative and the conservation equation is given through: pL + pqG + S = M,

(2.3)

where M is a constant which represents the total population (in concentrations) of all species. Remark 2.2. Eqs. (2.1) can be considered as a generalization of those proposed in [2] in that the model in [2] specifically assumes that p = 2 and q = 3 in the chemical-like reactions of the three species and takes a specific expression for ci = ci (γ˙ ) as ci = ai + bi γ˙ 2 with ai ≥ 0 and bi ≥ 0 constants for i ≥ 1. Such restrictions are removed in our model. 2.1. Asymptotic analysis of homogeneous three-species model In this section, we present the stability result for the three-species model. More precisely, we shall prove that ∀M > 0, each species converges to a unique steady state solution of the system (2.1) if initial conditions are given legitimately. The legitimacy conditions for the initial condition for L, G and S can be stated as follows: pL + pqG + S = M

and L, G, S ≥ 0.

(2.4)

Let (L∞ , G∞ , S∞ ) be the equilibrium state. We begin our analysis with the local behavior of solutions near the equilibrium state (L∞ , G∞ , S∞ ). 2.1.1. Local analysis of three-species model ˙ = 0 and S˙ = 0 and obtain the For locally analyzing the three-species model (2.1) at the steady state, we set L˙ = 0, G following identities: G=

c2 (γ˙ ) c4 (γ˙ )

 Lp

and S =

pc1 (γ˙ )L + β c2 (γ˙ )Lq pc3 (γ˙ )

.

(2.5)

Plugging the aforementioned relations (2.5) into the conservation equation (2.3), we obtain the critical point (L∞ , G∞ , S∞ ) given as follows: G∞ =

c2 (γ˙ ) c4 (γ˙ )

 Lq∞

and

S∞ =

c1 (γ˙ ) c3 (γ˙ )

L∞ +

β c2 (γ˙ ) q L∞ , pc3 (γ˙ )

(2.6)

4

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

where L∞ ∈ (0, M) is the unique positive root of the equation for x given as follows: px +

pqc2 (γ˙ ) c4 (γ˙ )

 q

x +

c1 (γ˙ ) c3 (γ˙ )

β c2 (γ˙ ) q x = M. pc3 (γ˙ )

x+

(2.7)

By a simple calculation, the Jacobian matrix J = J (γ˙ , L∞ , G∞ , S∞ ) of the system (2.1) at the critical points (L∞ , G∞ , S∞ ) can be shown to be of the following form:

−c1 (γ˙ ) − q2 c2 (γ˙ )Lq∞−1 J = qc2 (γ˙ )Lq∞−1 pc1 (γ˙ ) 

α c4 (γ˙ ) −c4 (γ˙ ) β c4 (γ˙ )

pc3 (γ˙ )Sp∞−1  0 −p2 c3 (γ˙ )Sp∞−1

 (2.8)

and the Jacobian matrix J has three distinct real eigenvalues λi for i = 1, 2 and 3 given as follows:

λ1 = 0,

λ2 = −

  1 X + X 2 − 4Y , 2

and λ3 = −

  1 X − X 2 − 4Y , 2

(2.9)

where X = c1 (γ˙ ) + c4 (γ˙ ) + q2 c2 (γ˙ )Lq∞−1 + p2 c3 (γ˙ )Sp∞−1 , q −1 Y = c1 (γ˙ )c4 (γ˙ ) + qc2 (γ˙ )c4 (γ˙ )(q − α)Lq∞−1 + p2 c3 (γ˙ )(c4 (γ˙ ) + q2 c2 (γ˙ )L∞ )Sp∞−1 .

It is apparent that two eigenvalues λ2 and λ3 are negative, from which we conclude that for any given M, the solution to the system (2.1) is convergent to the critical points (L∞ , G∞ , S∞ ) locally. 2.1.2. Global analysis of three-species model In this section, we provide a global analysis, i.e., the three-species model admits a unique steady state for any legitimate initial conditions for L, G, and S, (2.4). First, we provide a simple, but important lemma: Lemma 2.1. The homogeneous three-species model guarantees the non-negativity of all three populations, L, G and S for all time t ≥ 0. Proof. We prove the result by showing that whenever one of three species becomes zero, it becomes positive immediately in the later time. Let L = 0 at some time level, then, we have that dL dt

= c3 (γ˙ )Sp + α c4 (γ˙ )G > 0.

(2.10)

Similarly, let G = 0, then, we have dG dt

= c2 (γ˙ )Lq > 0.

(2.11)

Lastly, let S = 0, then, we have dS dt

= pc1 (γ˙ )L + β c4 (γ˙ )G > 0.

(2.12)

These inequalities (2.10), (2.11), and (2.12) imply that the solution to the system of equations should not be negative in the phase space R3 . This completes the proof.  Lemma 2.1 indicates that the model is physically correct. It implies that the solutions L, G and S must be bounded above as well. This is because the system is conservative, i.e., it holds that pL + pqG + S = M. To complete our main stability result, together with the aforementioned boundedness results, we need to invoke Dulac’s criterion and Poincare–Bendixson lemma regarding the dynamical system given as follows: dx dt

= f (x, y) and

dy dt

= g (x, y).

(2.13)

The first result we shall invoke is so-called, the Bendixson–Dulac Criterion [11]: Lemma 2.2 (Bendixson–Dulac Criterion). Let M be a simply connected region of the phase plane R2 . If there exists a continuously differentiable function φ(x, y) such that

∂ ∂ (φ(x, y)f (x, y)) + (φ(x, y)g (x, y)) is of a single sign in M, ∂x ∂y then the dynamical system (2.13) has no closed orbits completely contained in M .

(2.14)

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

5

Fig. 2.2. M = {L ≥ 0, S ≥ 0, pL + S ≤ M} ∈ R2 , P = (L∞ , S∞ ) and C is a closed orbit contained completely in the complement of B , as a subset of M .

In terms of the three-species model, the region M can be understood as follows:

M = {(L, S) ∈ R2 : L ≥ 0, S ≥ 0, and pL + S ≤ M}.

(2.15)

Next, we invoke the Poincare–Bendixson theorem. In literature, there are a number of versions of Poincare–Bendixson theorem. We shall take the following version, which is a bit rephrased from what is given in [11] for our purpose: Lemma 2.3 (Poincare–Bendixson Lemma). We assume that R is a closed, bounded subset of the phase plane R2 , f and g are continuous and differentiable on an open set containing R and R does not contain any fixed point of the system (2.13). Then, it holds true that if R does not contain a closed orbit, then there is no trajectory C that is completely confined in R. One can understand the closed and bounded subset of the phase plane R2 , denoted by R as the subset of the phase plane R2 given as R = M / B , where B is an open neighborhood of the unique critical points of the three-species model (2.1) such that if an initial condition is chosen in B , the system converges to the critical points. Under this setting, Fig. 2.2 is provided to illustrate the result in Lemma 2.3. We now provide the main theorem in this section. Theorem 2.4. The homogeneous three-species model admits a unique steady state solution no matter how initial conditions are given legitimately in the sense of Eq. (2.4). Proof. Using the conservation equation: pL + pqG + S = M, we reduce the three-species model (2.1) as a system of two equations as follows: dL dt

= f (L, S) and

dS dt

= g (L, S),

(2.16)

where f (L, S) = −c1 (γ˙ )L − qc3 (γ˙ )Lq + c2 (γ˙ )Sp + g (L, S) = pc1 (γ˙ )L − pc2 (γ˙ )Sp +

β c4 (γ˙ ) pq

α c4 (γ˙ ) pq

(M − pL − S) ,

(M − pL − S) .

We note that the system (2.16) has two eigenvalues, which are λ2 and λ3 in (2.9). Invoking the well-known local existence and uniqueness theory for autonomous systems [11], the negativity of two eigenvalues λ2 and λ3 guarantees that there is a (open) neighborhood B ⊂ R2 containing the critical point P = (L∞ , S∞ ) such that the solution to the system (2.16) with any initial condition in B converges to (L∞ , S∞ ). We shall now claim that in fact, if we choose any initial condition in the entire domain of phase space M , the solution to the system (2.16) converges to the critical point (L∞ , S∞ ). This will complete the proof. First of all, we can show that there is no closed orbit in M . Note that the space M is simply connected. We also have that

∂ ∂ (f (L, S)) + (g (L, S)) = −c1 (γ˙ ) − q2 c3 (γ˙ )Lq−1 − p2 c2 (γ˙ )Sp−1 − c4 (γ˙ ) < 0. ∂L ∂S

(2.17)

The Bendixson–Dulac’s criterion, Lemma 2.2 indicate that with φ(L, S) = 1, there are no closed orbits completely contained in M . Secondly, invoking Poincare–Bendixson Lemma 2.3, we conclude that if we start at an initial condition contained in R = M \ B . Then, the trajectory should not be contained completely in R. Therefore, the trajectory eventually meets with the neighborhood B and it will converge to the steady state solution P. This completes the proof. 

6

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

Fig. 2.3. Shear thickening followed by shear thinning in the homogeneous flow.

2.2. Model prediction and two-species model Eqs. (2.5) and (2.6) can be used to analytically compute the equilibrium state, L∞ , G∞ and S∞ for any given M. Typically, the homogeneous rheology is also of interest. Some of the sample model predictions obtained by the homogeneous three-species model are presented in Fig. 2.3. The results show that the shear-thickening behavior can be obtained with appropriate parameters. For the result, we have chosen M = 100, α = β = 2, µL = 0.5, µG = 100, µS = 0.01 and c1 (γ˙ ) = 102 + 101.5 γ˙ 2 ,

c2 (γ˙ ) = 10−3 γ˙ 2 ,

c3 (γ˙ ) = 10 + 10

c4 (γ˙ ) = 10 + 10

2

−1

γ˙ , 2

3

(2.18a) 1.5

γ˙ . 2

(2.18b)

Note that the effective viscosity denoted by µeff is computed by the formula: µeff = µL L∞ + µG G∞ + µS S∞ by setting the Newtonian viscosity µn = 0. We remark that in general, the two-species model, i.e., the reduced model of the three-species model (2.1) can only predict a constant or strictly monotonic viscosity as a function of γ˙ . More precisely, we consider two-species models by taking into account the following two independent cases. Case I: one of micelles is absent and the rates to and from itself are set to be 0; Case II: the gelation G is absent and the rates to and from G are set to be 0. Two independent cases Case I and Case II can be written as the following two-species model in its general form as follows: dV dt

  = −k cV (γ˙ )Vk − cW (γ˙ )W

and

dW dt

= cV (γ˙ )Vk − cW (γ˙ )W,

(2.19)

where k > 1, ci = ci (γ˙ ) ≥ 0, with i = V or W, which is coupled with the conservation equation: V + kW = M.

(2.20)

Then, we can show that the system of Eqs. (2.19) has the following characteristics [12]: 1. The system leads to non-negative solutions V and W, which is globally stable and asymptotically convergent to the unique steady state solution; and 2. The effective viscosity µeff = µV V∞ + µW W∞ is either a constant or strictly monotonic as a function of γ˙ . We note that similar results have been provided by Handzy [2], but only using numerical experiments. There has not been a rigorous mathematical justification. The reduced two-species model cannot be appropriate for the shear-thickening transition for the homogeneous fluids. A further limitation of the reduced models is discussed in Remark 4.1.

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

7

3. Inhomogeneous three-species model In this section, we introduce the inhomogeneous three-species model. By the inhomogeneity, we mean both the spatial dependence of the populations on the flow and their effects toward the flow. We shall, therefore, incorporate the dependence of the populations on the flow and the dependence of the flow on the populations. First of all, we extend the homogeneous three-species model into an inhomogeneous three-species model. For this purpose, we take natural two-step process: (1) To extend the three-species model (2.1) into the partial differential equations using Fick’s law [13] and (2) to change the time derivative as the material or convective derivative along the flow velocity. The Fick’s law is often used for such a purpose, see [14] and the references cited therein. More precisely, this consideration results in the following partial differential equations for the three species, in Ω , Dt L = −c1 (γ˙x ) L − qc2 (γ˙x ) Lq + c3 (γ˙x ) Sp + α c4 (γ˙x ) G + ϵL 1L

(3.1a)

Dt G = c2 (γ˙x ) Lq − c4 (γ˙x ) G + ϵG 1G

(3.1b)

Dt S = pc1 (γ˙x ) L − pc3 (γ˙x ) Sp + β c4 (γ˙x ) G + ϵS 1S,

(3.1c)

where Dt = ∂/∂ t + u · ∇ is the material derivative with u being the fluid velocity, ϵL , ϵG and ϵS are diffusion coefficients, and γ˙x is the spatially dependent shear rate or velocity gradient given as follows:

γ˙x :=

√ 

1

(∇ u + (∇ u)T ) (3.2) 2 and D (u) : D (u) is the Frobenius inner product between matrix D (u). In addition, we have the following corresponding expression for the coefficients ci for i = 1, 2, 3 and 4 2 D (u) : D (u),

with D (u) =

1. ci = ci (γ˙x ) ≥ 0 for i = 1, 2, 3, 4; 2. α and β are positive numbers satisfying the equality pα + β = pq, and we note that the addition of the diffusion terms require boundary condition for L, G and S. These will be given as the pure Neumann conditions, namely,

∇ L · n = ∇ G · n = ∇ S · n = 0,

on ∂ Ω ,

(3.3)

where n is the unit outward normal vector on ∂ Ω . Under this condition, the system (3.1a), (3.1b) and (3.1c) can be shown to be still conservative, i.e., pL + pqG + S = M. Eqs. (3.1) model how the flow affects each species, L, G and S. In the inhomogeneous populations model, γ˙ is replaced by the spatially varying shear rate, denoted by γ˙x . Secondly, we consider how the populations L, G and S affect the flow in turn. We begin by reviewing a micro–macro modeling of generic viscoelastic models, where we are able to identify that the contribution of the populations can affect the polymeric viscosity, which then affects the fluid velocity. More precisely, the occurrence of the gelation will be used to increase the polymer viscosity, which in turn affects the computation of the stress and therefore, affects its flow pattern. We consider the flow of the dilute polymeric fluids governed by conservation of mass:

∇ · u = 0,

in Ω ,

(3.4)

and conservation of momentum:

ρ



 ∂u + u · ∇ u = ∇ · σ, ∂t

in Ω ,

(3.5)

where ρ is a fluid’s density and σ is the total stress, which is decomposed into its polymer and Newtonian (solvent) contributions as follows:

σ = −pδ + 2µn ∇ · D (u) + τ,

in Ω ,

(3.6)

where µn is the Newtonian viscosity. To obtain the polymeric contribution to the stress, denoted by τ , we consider the distribution of bead–spring dumbbells (see [15] for more details). The symbol r denotes the location of the center of mass of the dumbbell and the beads are located at Rν with ν = 1, 2. Note that

(−1)ν

Q, with Q = r2 − r1 = R2 − R1 , 2 where Q is called the connector vector. The probability that a dumbbell has center of mass r and ν -th mass located at Rν is given by the configuration density function rν = r + Rν = r +

ψp (r − Rν , Q, t ).

(3.7)

To find the constitutive equation, we need to know the forces on the dumbbell. Typically, three forces are taken into account: 1. The Spring force: FνS = −(−1)ν HQ, where H is a spring constant.

(3.8)

8

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

2. The Brownian force: FνB = −kT

∂ ln ψp (r − Rν , Q, t ), ∂ rν

(3.9)

where k is the Boltzmann constant and T the temperature. 3. The Hydrodynamic force: FνH = −ζ [[r˙ν ]]p −vν − vν ,





(3.10)

where ζ is the drag coefficient, [[r˙ν ]] is the momentum-weighted average of the velocity of the ν th bead over configuration space as defined in [15], vν is the fluid velocity imposed on the dumbbell as evaluated at the position of the ν th bead and  vν is taken to be p

1

 vν = −ξ D (u) · Rν = − ξ (−1)ν D (u) · Q, 2

(3.11)

where ξ is a slippage parameter. By setting the total force on the beads is zero and by the diffusion equation of ψp , we arrive at the equation for the second moment {QQ} as follows: 4H

ζ

{QQ} + {QQ} −

kT 2ζ

∆{QQ} =

4nkT

ζ

δ,

(3.12)

where (·) is the Gordon–Schowalter derivative given by the following: for any w ∈ Rd×d , with d being the dimension of Ω , w = Dt w − R (t )w − wR (t )T ,

(3.13)

where n is the number density of the polymer and R (t ) =



a+1

2

∇u +

a−1

2

T

∇u



with the symbol a = 1 − ξ .

(3.14)

Let τ = aH {QQ} be the conformation tensor and λH = ζ /4H, the molecular relaxation time. We have the equation for the polymeric stress τ :

τ + λH τ  − λH Dtr 1τ =

µp aλH

δ,

(3.15)

where Dtr is a molecular diffusivity coefficient, and µp = a2 nkT λH is the polymer contribution to the viscosity.1 We remark that the number density is the number of polymer per unit volume. Namely, the population is found to contribute to the number density. Therefore, we assume that each population contributes to µp as follows in our model:

µp = µL L + µG G + µS S,

(3.16)

where µL , µG and µS are viscosities of long micelle, gelation and short micelle, respectively. We remark that considering three different species, it is more natural to consider the multimode models by taking into account different relaxation times for different species, λL , λG and λS . However, for the sake of simplicity, we have adopted only a single mode constitutive equation. For the future research, we would consider the incorporation of multimode models of the constitutive equations. In summary, we arrive at the newly proposed three-species model as follows: Dt L = −c1 (γ˙x ) L − qc2 (γ˙x ) Lq + c3 (γ˙x ) Sp + α c4 (γ˙x ) G + ϵL 1L,

(3.17a)

Dt G = c2 (γ˙x ) Lq − c4 (γ˙x ) G + ϵG 1G,

(3.17b)

Dt S = pc1 (γ˙x ) L − pc3 (γ˙x ) Sp + β c4 (γ˙x ) G + ϵS 1S,

(3.17c)

∇ · u = 0,

(3.17d)

ρ Dt u = ∇ · σ, σ = −pδ + 2µn ∇ · D (u) + τ, µp τ + λH τ  − λH Dtr 1τ = δ, and µp = µL L + µG G + µS S, aλH

(3.17e)

subject to appropriate boundary conditions for both u and τ and pure Neumann conditions for L, S and G. Note that there can be variety of choices for ci , α, β, p, q and ϵL , ϵG , ϵS . Such choices should be designed to incorporate the relevant physics. In Section 4, we provide some special choices of these parameters to study the shear-thickening transitions in a shear flow. We remark that the design of ci for i = 1, 2, 3, 4 is particularly crucial.

1 In [16], it was called the characteristic polymer viscosity scale.

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

9

Fig. 4.1. Numerical modeling setting: a simple shear flow between two vertical plates. Transition from Homogeneous Simple Shear Flow to Inhomogeneous Shear Flow—Shear Banding Transition. The flow rate γ˙ = V /dist is provided on the left plate while the right plate is fixed, where V is the target velocity of the left plate while dist denotes the distance between two plates.

4. Three-species model: A shear thickening transition in shear flows of wormlike micellar fluids As discussed in a number of literatures [4,1], the main characteristics of a shear-thickening transition in Taylor–Couette flow of a wormlike micellar fluids can be summarized as the following three regimes (see Fig. 1.1): Regime I: At low shear rate γ˙ < Γc , the viscosity is about the constant as a function of shear rate, indicating the Newtonian behavior. Regime II: At the critical shear rate γ˙ ≈ Γc , the viscosity increases and deviates progressively from the Newtonian behavior. Regime III: Viscosity passes through its maximum, at a level several times that of the solvent viscosity, and then shear thinning is observed. The main result in this section is to identify parameters ci for i = 1, 2, 3, 4 for the three-species model that can predict all three regimes above and present the relevant results. We note that even if the actual experiment has been performed for the Taylor–Couette flow, our modeling is restricted to a simple shear flow of a wormlike micellar fluids. Therefore, throughout this section, we consider a unidirectional simple shear flow inside two vertical plates. We denote the distance between two plates by dist. Note that dist will be assumed to be one throughout this section. The shear flow will be generated by moving the left plate upward (namely, providing a nonzero γ˙ ). The right plate is always assumed to be fixed (see Fig. 4.1 for illustration). The coordinate of the left plate is x = 0 and that of the right plate is given as x = dist. A number of flow rates γ˙ have been applied upward at the left plate. We note that the shear flow will be generated by moving the left plate slowly until it reaches the desired shear rate in time similar to the start-up experiment considered in [7]. Namely, let γ˙ = V /dist be the target shear rate. We provide the flow velocity, at the left plate as a boundary condition from the formula: V tanh(ν t ) where the quantity 1/ν represents roughly the time for the velocity applied at the left plate to reach its final (imposed) shear rate. Note that ν = 10 is chosen in our computations, which is known to be the similar order with that for the experiments [17,7]. The initial conditions for L and G are set to zero while S = M (see (2.3)). A numerical method for solving the model (3.17) is based on the classical finite element method [18,19]. The piecewise linear finite element has been applied for approximating the velocity fields and piecewise constant fields have been used for approximating both polymeric stress fields and three species L, G and S. In fact, we have attempted the finite difference method for both equations as well, which resulted in the identical results. In particular, positivity preserving schemes for the nonlinear solver for the three-species model have been designed so as to obtain physically relevant results [20]. The mesh size we have used is h = 0.01 and 0.005 and the temporal step size is dt = 0.001. The mesh convergence has been tested and validated. For the steady state computation, we have computed solutions up to the time level t = 3000, at which the steady state is definitely reached. Justification of our numerical schemes has also been made in comparison with the analytic solutions (see Fig. 4.4 and Section 4.2 for more details). Detailed numerical analysis for the employed numerical techniques together with their parallel implementations will be the topic of the future publication. In our modeling of the shear-thickening transitions, we assume that the critical shear rate Γc is one. We keep most of parameters fixed in the mathematical equations (3.17), which are 1. p = 2, q = 3, ϵL = ϵG = ϵS = 10−3 , α = β = 2, 2. µn = 0.05, ρ = 1, Dtr = 0, λH = 1 and ξ = 0, or equivalently a = 1 and vary only few parameters, which are µL , µG , µS and γ˙ . 4.1. Conjecture by Lerouge and Berret [1] and Shear Rheology of the three-species model In this section, we discuss the minimal scenario of the shear-thickening transition proposed by Lerouge and Berret [1]. By investigating the shear rheology of the three-species model, we design and present the parameters that can predict such proposed scenario. Lerouge and Berret [1] proposed a schematic diagram accounting for a possible scenario of the shear-thickening transition as described in Fig. 4.2 (Left). The phase A denotes the flow curve for the shear-thickening fluid, while the phase B denotes the

10

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

Fig. 4.2. A schematic diagram of phase A and phase B (Left) provided by Lerouge and Berret [1] and an indicator function centered at Γc = 1 (Right).

flow curve for the shear-thinning fluid. In the proposed scenario, it was conjectured that starting at the Newtonian regime below a certain shear rate, at the critical shear rate, the fluid jumps from the phase A to the phase B. Basically, this means that the fluid remains as Newtonian at the low shear rate (phase A) and at the shear rate exceeding the critical shear rate, the fluid becomes shear-thinning (phase B). Note that shear thinning fluids are known to induce shear banding behaviors [21,22]. Since the gelation observed to occur at Regime II starts at the fixed inner cylinder and apparently induces a shear banding, the phase B can be correlated with the shear banding and the gelation and therefore, the increasing viscosity. In our modeling of the shear-thickening transition, we shall realize the scenario proposed by Lerouge and Berret [1] and therefore, we design model parameters that exhibit the phases A and B as well as their transitions. To this end, we have to identify the relation between the shear stress, strain rate and the effective viscosity. Following the argument in [23,24], we can consider the reformulation of the constitutive equation as an integral equivalent model as follows:

τ=



t

−∞

  µp −(t − s) exp F (s, t )F (s, t )T ds, λH λ2H

(4.1)

where F is the deformation tensor. This is because we set ξ = 0. At the steady simple shear flow, the deformation gradient is given by the following formula:

 dv

F (s, t ) =

dx

1

0

(t − s)



1

,

(4.2)

where v is the vertical velocity. Note that the simple shear flow of our concern is obtained by providing a flow velocity for the left plate while the right plate is fixed. We note that a simple calculation of the shear stress is given as follows:

σ xy = µn γ˙ + τ xy      t 1 (t − s) (t − s) ds γ˙ = µn + lim 2 µp (s) exp − t →∞ λ λH −∞ H = (µn + µp )γ˙ ,

(4.3)

which identifies the relation between the shear stress and the shear rate. In addition, the effective viscosity is found to be given by the following formula:

µeff = µn + µp = µn + µG G∞ + µL L∞ + µS S∞ .

(4.4)

We remark that in our modeling, the zero shear rate µ0 is defined to be the solvent viscosity plus the viscosity due to the small micelles as given below:

µ0 = µn + µS M .

(4.5)

We set M = 100, µL = 0.1, µG = 50 and µS = 0.01 in this section. By careful investigation, we have identified the parameters ci (γ˙ ) for i = 1, 2, 3, 4 that can induce the phases A and B of the fluid, respectively as follows: c1 (γ˙ ) = 104 + 106 γ˙ 2 , c3 (γ˙ ) = γ˙ , 2

c2 (γ˙ ) = 102 + 10γ˙ 2 ,

c4 (γ˙ ) = 10 + 10 2

−3

γ˙ , 2

(4.6a) (4.6b)

and c1 (γ˙ ) = 104 + 106 γ˙ 2 , c3 (γ˙ ) = 10 + γ˙ , 3

2

c2 (γ˙ ) = 102 + 10γ˙ 2 , c4 (γ˙ ) = 10 + 10 2

−3

γ˙ . 2

(4.7a) (4.7b)

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

11

Fig. 4.3. The phase A and phase B; Parameters used are described in (4.6) and (4.7). The critical shear rate Γc is set to be 1.

Fig. 4.4. Shear flow under low shear rate and populations from inhomogeneous model (Left and Middle) and analytic expressions of populations for the homogeneous model as a function of shear rate for the phase B parameters. µL = 0.1, µG = 50, µS = 0.01 and M = 100.

As shown in Fig. 4.3 (Left), the parameter (4.6) predicts the phase A while the parameter (4.7) predicts the phase B. We note that Fig. 4.3 (Right) shows that the graph of Gelation population is very similar to that of the effective viscosity, µeff as well. An outstanding issue is how to obtain a transition from phase A to phase B. We note that the unique difference in the parameter designs for each phase A and B can be found at the parameter c3 (γ˙ ). As discussed in Remark 2.1, it is related to the generation of the cylindrical micelles, thereby producing the gelation. The phase A is obtained with c3 (γ˙ ) = γ˙ 2 . This is because it takes small values at γ˙ < Γc = 1 unlike the case of phase B. Therefore, we conclude that the transition between two phases A and B can be attained by providing the transition in the definition of c3 (γ˙ ). Namely, we shall define c3 (γ˙ ) as follows: c3 (γ˙ ) =



γ˙ 2 ,

103 + γ˙ 2 ,

if γ˙ ≤ Γc , if γ˙ > Γc .

(4.8)

On the other hand, Eq. (4.8) defines c3 (γ˙ ) discontinuously as a function of γ˙ , which can be avoided by introducing a smoothed Heaviside function as follows: c3 (γ˙ ) = 103 × phase(γ˙ , Γc ) + γ˙ 2 ,

(4.9)

where phase(γ˙ , Γc ) :=

1 2

(1 + tanh (κ (γ˙ − Γc ))) .

(4.10)

Note that the symbol κ controls the varying sharpness between 0 and 1 (see Fig. 4.2). The shear rate γ˙ considered for the homogeneous model in (4.9), is the applied shear rate on the left plate as a boundary condition in the shear flow, i.e., γ˙ = V /dist. The parameters ci (γ˙ ) should be replaced by ci (γ˙x ) in the inhomogeneous model for i = 1, 2, 3, 4 and they will be given in the following manner: c1 (γ˙x ) = 104 + 106 γ˙x2 ,

c2 (γ˙x ) = 102 + 10γ˙x2 ,

c3 (γ˙x ) = 10 × phase(γ˙ , Γc ) + γ˙ , 3

2 x

(4.11a)

c4 (γ˙x ) = 10 + 10 2

−3

γ˙ . 2 x

(4.11b)

12

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

Fig. 4.5. Homogeneous three-species model: Gelation and Shear Diagram with varying Gelation Viscosity µG with the phase B parameters. M = 100 and

µL = 0.1, µG = 50–300, µS = 0.01.

The choice of parameter (4.11) is based on the scenario that for γ˙ larger than Γc , the fluid takes the phase B while it takes the phase A for γ˙ smaller than Γc . Physically, the new definition of c3 (4.11) indicates that the populations can feel the neighboring flow rate. Namely, we provide a non-locality in the model. Having identified appropriate parameters and their definitions according to the conjecture by Lerouge and Berret [1], we are in a position to apply them to simulate the inhomogeneous Three-Species model in an attempt to reproduce all regimes observed in the shear-thickening transition, in an order of Regime I and Regimes II, III. 4.2. Newtonian Regimes: Regime I In the regime where γ˙ is smaller than Γc , the flow behaves like Newtonian. Namely, the shear banding does not occur and the gelation population remains negligible. We observe that the shear banding does not occur for the phase B fluid at the lower shear rate regime either. Fig. 4.4 (Left) shows that for γ˙ < Γc , the viscosity remains as constant. However, an undesired generation of gelation is predicted in the phase B. Fig. 4.4 (Right) shows analytic steady solutions to the homogeneous three-species model, L∞ , G∞ and S∞ . At γ˙ = 0.05, it predicts relatively large gelation, i.e., G = 14.900. Therefore, the phase B is not appropriate and the phase A should be taken into account in this regime to obtain Newtonian flow with negligible gelation populations at lower shear rate regime. Note that Fig. 4.4 (Middle) plots the steady state three species obtained by solving the inhomogeneous three-species model. The results are very close to the analytic prediction of the steady populations of three species. The maximum relative difference between the analytic steady state and the numerical steady state from the inhomogeneous model is around 10−2 . Note that the inhomogeneous three-species model incorporates the nonzero diffusion terms (ϵL = ϵG = ϵS = 10−3 ) unlike the homogeneous case. This confirms the reliability of our numerical simulations. 4.3. Shear thickening and thinning regime: Regimes II, III The shear-thickening occurs at γ˙ , exceeding certain critical shear rate Γc [4]. When shear-thickening transition occurs, the gelation is observed to occur starting at near the fixed right plate and expands toward the moving left plate as time progresses. Namely, in the second regime, it is expected that the gelation occurs, with a correlation with the shear banding. The shear banding is often attributed to the constitutive instability as observed in the Johnson–Segalman model, namely, the non-monotone stress–strain relation for some parameter ranges [22,7]. The Johnson–Segalman model can be obtained from the three-species model (3.17) by setting µp as constant. In fact, as illustrated in Fig. 4.5 (Right), the phase B induces the non-monotone shear stress–strain relation. A sample result shows that the three-species model predicts the shear banding (see Fig. 4.6) as well. The observed shear banding is clearly correlated with the gelation. Experimentally, the gelation is observed largely at low shear rate region in the shear banded flow. Fig. 4.6 (Middle and Right) clearly shows such a well-defined correlation. Remark 4.1. Typically, the shear rheology (see Fig. 4.5) exhibits a Newtonian behavior in the regime where γ˙ is small. Such a result is due to the plateau of the gelation population. The constant gelation followed by the decreasing gelation can be obtained by the three-species model. Any of the reduced two-species model of the three-species model cannot predict such results (see (2.19)) [12]. Fig. 4.5 demonstrates that the non-monotone stress–strain relation can be significantly affected by the gelation viscosity. Namely, the larger the gelation viscosity, the larger the region where the shear stress decreases as a function of shear rate. More precisely, we denote the local minimum and maximum of the shear stress by τm and τM , respectively and denote the

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

13

Fig. 4.6. Shear banding obtained by the inhomogeneous three-species model for two different flow velocities given as 5 and 7, respectively. The phase B parameters are used together with M = 100 and µL = 0.1, µG = 50, µS = 0.01.

Table 4.1 For each µG , γ˙τm , is the shear rate at which shear stress takes local min, γ˙τM is the shear rate at which shear stress takes local max, and γ˙bm and γ˙b are the minimum and maximum shear M rate, respectively, between which the computational shear banding is observed. The phase B parameters are used together with M = 100 and µL = 0.1, µS = 0.01.

µG

γ˙τm

γ˙τM

γ˙bm

γ˙bM

50 150

≈1.33 ≈1.31

≈11.51 ≈15.18

≈1.3 ≈1.1

≈13.1 ≈21.8

shear rates at which the shear stress takes τm and τM by γ˙τm and γ˙τM , respectively. Then, it is observed that the larger the µG , the smaller the γ˙τm and the larger the γ˙τM . Therefore, let µGi for i = 1, 2 be two gelation viscosities. Then it holds:

µG1 ≤ µG2 ⇒ [γ˙τm1 , γ˙τM1 ] ⊆ [γ˙τm2 , γ˙τM2 ],

(4.12)

where (γ˙τmi , τmi ) and (γ˙τM , τMi ) denote the local minimum and local maximum of the shear stress for µGi , for i = 1, 2 as a i function of the shear rate. Therefore, the gelation viscosity µG can be used as a control parameter for shear banding transition in the shear thinning fluid. Such a control parameter can be found at the slippage parameter ξ in the Johnson–Segalman model. Unlike the role played by µG , the larger the ξ , the smaller are both the γ˙τm and γ˙τM . Therefore, the relation (4.12) does not hold for the Johnson–Segalman model. We note that in fact, the shear banding is observed in the region that includes the interval [γ˙τm , γ˙τM ] [25] for the Johnson–Segalman model. This is also the case for the three-species model. Let γ˙bm and γ˙bM denote the minimum shear rate and maximum shear rate at which the banding occurs. Then, in our model, it is shown to hold that for a given µG (see Table 4.1),

[γ˙τm , γ˙τM ] ⊆ [γ˙bm , γ˙bM ].

(4.13)

Next, we investigate the dynamic behavior of the shear stress as well as the viscosity in the phase B at the prescribed intermediate shear rate in comparisons with experiments presented in [26,1] (see also Fig. 1.1). We note that typically, the kinetics of shear banding has been considered as a two-stage process. The first stage is known as tilting, at which stage, local shear rate changes smoothly across the gap with no discernible shear band. The tilt proceeds then toward high shear rate region at the left plate [4]. The second stage is the shear banding, at which stage, the low shear rate band starts and develops at the right plate. This then expands toward the left plate [27,5,28,1]. The shear banding mechanism described here is known to be rather general for entangled shear thinning fluids. Such two-step transition features are qualitatively captured in our three-species model as shown in Fig. 4.7. From Fig. 4.7, we can easily expect that the evolution of the gelation progresses toward the left plate, starting at the right plate as time proceeds. Fig. 4.8 shows the dynamics of the stress and viscosity evolutions obtained from the three-species model. Lastly, the phase B could still be used to obtain Regime III, which is in fact a relatively easy regime. From shear-gram shown in Fig. 4.6, we expect that the flow behavior will remain as Newtonian at the larger shear rate. Namely, at the larger γ˙ , the shear stress is linear as a function of the shear rate and the gelation population will be negligible as well. In fact, the series of simulations confirm that the minimal scenario of the shear-thickening transition proposed by Lerouge and Berret [1] is acceptable. Fig. 4.9 demonstrates the prediction of the shear-thickening transition as a function of the shear rate from the parameter designed and proposed in (4.11). The results agree qualitatively well with the experimental data [4].

14

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

Fig. 4.7. The kinetics of shear banding; the First Stage (Left) and the Second Stage (Right). The phase B parameters are used together with M = 100, µL = 0.1, µG = 50, µS = 0.01 and γ˙ = 5.

Fig. 4.8. Evolution of average effective viscosity (Left) and average shear stress (Right). The phase B parameters are used together with µL = 0.1, µG = 50, µS = 0.01, M = 100 and γ˙ = 5.

Fig. 4.9. Shear thickening transition in shear flow for non-local inhomogeneous three-species model: Transition parameters (4.11) are used together with

µL = 0.1, µS = 0.01 and M = 100, while µG = 50 (Left) and µG = 150 (Right), respectively.

5. Conclusion We have presented new population-dynamics based three-species model in this paper, together with some mathematical analysis. Our model has been shown to predict a shear-thickening transition, qualitatively agreeable with experimental results. The model parameter selection takes into account the scenario conjectured by Lerouge and Berret [1]. Therefore, the presented result in this paper confirms the conjecture proposed by Lerouge and Berret [1] affirmatively. More interesting and challenging future works will be to use the proposed three-species model to investigate the oscillatory falling sphere in wormlike micellar fluids [8,29] as well as jumping bubbles [9]. H. Kang and Y.-J. Lee would like to acknowledge the helpful discussion of Dr. Handzy and his friendship. Acknowledgments The authors would like to thank two referees for their constructive comments that lead to significant improvements of the paper. This research is in part supported by NRF 2012R1A1A1014924 and 2015R1A3A2031159. The research of the second author was supported in part by NSF DMS-1358953.

H. Kang, Y.-J. Lee / Computers and Mathematics with Applications (

)

–

15

References [1] Sandra Lerouge, Jean-Francois Berret, Shear induced transitions and instabilities in surfactant wormlike micelles, in: Polymer Characterization, Rheology, Laser Interferometry, Electrooptics, Spring-Verlag, 2010, pp. 1–71. [2] Handzy Nestor, Experimental observations and mathematical description of micellar fluid flow (Ph.D. thesis), The Pennsylvania State University, 2005. [3] H. Rehage, H. Hoffmann, Shear induced phase transitions in highly dilute aqueous detergent solutions, Rheol. Acta 21 (1982) 561. [4] C. Liu, D.J. Pine, Shear-induced gelation and fracture in micellar solutions, Phys. Rev. Lett. 77 (1996) 2121. [5] Y. Hu, E.F. Matthys, Characterization of micellar structure dynamics for a drag-reducing surfactant solution under shear: normal stress studies and flow geometry effects, Rheol. Acta 34 (1995) 450–460. [6] J.-F. Berret, R. Gamez-Corrales, J. Oberdisse, L.M. Walker, P. Linder, Flow-structure relationship of shear-thickening surfactant solutions, Europhys. Lett. (1998). [7] L. Zhou, G.H. McKinely, L.P. Cook, Wormlike micellar solutions: III. VCM model predictions in steady and transient shearing flows, J. Non-Newton. Fluid Mech. 211 (2014) 70–83. [8] A. Jayaraman, A. Belmonte, Oscillations of a solid sphere falling through a wormlike micellar fluid, Phys. Rev. E 67 (2003). [9] Nestor Z. Handzy, Andrew Belmonte, Oscillatory rise of bubbles in wormlike micellar fluids with different microstructure, Phys. Rev. Lett. 92 (2004). [10] J.L. Goveas, P.D. Olmsted, A minimal model for vorticity and gradient banding in complex fluids, Eur. Phys. J. E 6 (2001) 79–89. [11] S.H. Strogatz, Nonlinear Dynamics and Chaos with Applications to Physics, Biology, Chemistry andEngineering, Perseus Books, 1994. [12] Hunseok Kang, Young-Ju Lee, An observation on two-species models, reduced models of the three-species model, in preparation. [13] A. Fick, On liquid diffusion, J. Membr. Sci. 100 (1995) 33–38. [14] R. deForest, A. Belmonte, Spatial pattern dynamics due to the fitness gradient flux in evolutionary games, Phys. Rev. E 87 (2013) 062138. [15] R.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Volume 2: Kinetic Theory, Weiley Interscience, New York, 1987. [16] L.P. Cook, L.F. Rossi, Slippage and migration in models of dilute wormlike micellar solutions and polymeric fluids, J. Non-Newton. Fluid Mech. 116 (2004) 347–369. [17] C.J. Pipe, N.J. Kim, P.A. Vascuez, L.P. Cook, G.H. McKinley, Wormlike micellar solutions. ii:comparison between experimental data and scission model predictions, J. Rheol. 54 (2010) 881–914. [18] D.N. Arnold, A concise introduction to numerical analysis. unpublished, 2001. Available at: http://www.ima.umn.edu/~arnold/597.00-01/notes.html. [19] Susanne C. Brenner, L. Ridgway Scott, The Mathematical Theory of Finite Element Methods, second ed., in: Texts in Applied Mathematics, vol. 15, Springer-Verlag, New York, 2002. [20] Y.-J. Lee, H. Kang, C.-S. Zhang, Scalable computation of the three species models for the wormlike micellar fluids, 2015, in preparations. [21] P.D. Olmsted, Two-state shear diagrams for complex fluids in shear flow, Europhys. Lett. 48 (1999) 339–345. [22] P.D. Olmsted, Perspectives on shear banding in complex fluids, Rheol. Acta 47 (2008) 283–300. [23] Y.-J. Lee, Jinchao Xu, New formulations, positivity preserving discretizations and stability analysis for non-newtonian flow models, Comput. Methods Appl. Mech. Engrg. 195 (2006) 1180–1206. [24] Young-Ju Lee, Jinchao Xu, Chen-Song Zhang, Stable finite element discretizations for viscoelastic flow models, in: Handbook of Numerical Analysis, Vol. XVI, 2011. [25] G.C. Georgiou, D. Vlassopoulos, On the stability of the simple shear flow of a johnson-segalman fluid, J. Non-Newton. Fluid Mech. 75 (1998) 77–97. [26] Y.T. Hu, P. Boltenhagen, E. Matthys, D.J. Pine, Shear thickening in low-concentration solutions of wormlike micelles. II slip, fracture, and stability of the shear-induced phase, J. Rheol. 42 (1998) 1209. [27] Y. Hu, C.V. Rajaram, S.Q. Wang, A.M. Jamieson, Shear thickening behavior of a rheopectic micellar solution:salt effects, Langmuir 10 (1994) 80–85. [28] Y.T. Hu, P. Boltenhagen, D.J. Pine, Shear thickening in low-concentration solutions of wormlike micelles. I direct visualization of transient behavior and phase transitions, J. Rheol. 42 (1998) 1185. [29] Y.-J. Lee, C.-S. Zhang, Self-sustaining oscillations of the falling sphere through a fluids governed by the johnson-segalman model, Soft Matter (2015) submitted for publication.