PHYSICAL REVIEW E 84, 051404 (2011)

Reentrant transition in the shear viscosity of dilute rigid-rod dispersions
Hideki Kobayashi* and Ryoichi Yamamoto†Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan and
CREST, Japan Science and Technology Agency, Kawaguchi 332-0012, Japan(Received 15 February 2011; revised manuscript received 7 September 2011; published 16 November 2011)
The intrinsic viscosity of a dilute dispersion of rigid rods is studied using a recently developed direct numericalsimulation (DNS) method for particle dispersions. A reentrant transition from shear-thinning behavior to thesecond Newtonian regime is successfully reproduced in the present DNS results around a Peclet number Pe = 150,
which is in good agreement with our theoretical prediction of Pe = 143, at which the dynamic crossover from
Brownian to non-Brownian behavior occurs in the rotational motion of the rotating rod. For Pe > Pec at which
the dynamic crossover occurs, the rigid rod undergoes non-Brownian rotational motion, which means that therotational frequency is proportional to the shear rate. This crossover occurs only when the effect of the shear rateis dominant over the effect of the thermal fluctuation for all areas. Otherwise, the rotational motion is dominatedby the thermal fluctuation. The viscosity undershoot is observed in our simulations before reaching the second
Newtonian regime. The physical mechanisms behind these behaviors are analyzed in detail.
DOI: 10.1103/PhysRevE.84.051404

PACS number(s): 83.80.Hj, 83.50.Ax, 83.60.Fg, 83.80.Rs

I. INTRODUCTION

The viscous properties of dilute dispersions of rigid rodschange drastically as the rate of applied shear flows γ̇
increases. Although many previous studies have investigatedthis phenomenon, the mechanism of this viscosity change is not
yet completely clear. The aim of this paper is to contribute to
the understanding of the detailed mechanism of the viscosity
change by performing direct numerical simulations (DNSs)
for a dilute dispersion of rigid rods that are subject to thermal
fluctuations in a Newtonian host fluid.
The relationship between the measurable bulk rheologicalproperties and the microscale description of dispersions ofrodlike particles has been previously investigated [1–9].
Giesekus obtained the expression for the bulk stress tensor
of diluted spheroidal dispersions under shear flow by taking
into account the effects of the rotational Brownian motionof the spheroids caused by thermal fluctuations [3]. Leal and
Hinch reported that the viscosity behavior is characterized by
the aspect ratio l of the rod and the dimensionless shear rateγ̇ /Dr , where Dr is the rotational diffusion constant [5–8].
In the case of weak-shear flow, γ̇ /Dr  1, the dilute rigid
rod dispersions exhibit the first Newtonian behavior, in which
the viscosity η of the dispersion is constant and equal to the first
Newtonian (zero-shear limiting) value η0 . For an intermediate
regime, 1  γ̇ /Dr  l 3 + l −3 , the dispersions exhibit shearthinning behavior, in which η ∝ (γ̇ /Dr )−1/3 . In the case ofstrong-shear flow, l 3 + l −3  γ̇ /Dr , the dispersions reenter
the second Newtonian regime, in which η becomes constantagain and is equal to the second Newtonian (high-shearlimiting) value η∞ . Similar results have also been obtainedin numerical [10,11] and experimental [12,13] studies. In
this paper the phrase viscosity transition is used to express
the changes in viscosity from the first Newtonian to shearthinning behavior and from shear-thinning to the second
Newtonian behavior. Similar results have been observed for*†

dilute dispersions of flexible chains, both experimentally [14]and theoretically [15].
Hinch and Leal [5,6] proposed a theoretical model for theviscosity transitions. They considered that the viscosity η ofthe dispersion is determined by the ensemble average of thetemporal viscosity η̂(θ,ϕ) using the probability distributionfunction (PDF) Pγ̇ (θ,ϕ) of the two orientation angles θ and ϕof the rod, i.e.,η(γ̇ ) = η̂(θ,ϕ)Pγ̇ (θ,ϕ) dθ dϕ.
(1)
Here the form of Pγ̇ (θ,ϕ) is shear rate dependent, and the shear
rate dependence of the dispersion viscosity η(γ̇ ) is introduced
mainly through this function.
The rigid rod undergoes a random rotational Brownian
motion at low shear rates in the first Newtonian regime, inwhich the effect of thermal fluctuations is dominant over theeffect of shear flow. Therefore,Pγ̇ (θ,ϕ) = const

holds over the entire phase space of θ and ϕ. The viscosity
is thus constant with respect to the shear rate change in thisregime, i.e., η(γ̇ ) = η0 .
In contrast, the rigid rod undergoes a deterministic tumbling
motion attributable to strong-shear flow in the second Newtonian regime. In this case, the tumbling motion is perfectly
described by Jeffrey’s equation [1]. Therefore, the PDF approaches the high-shear limiting (non-Brownian) asymptotic
form with increasing γ̇ :Pγ̇ (θ,ϕ) = PJ (θ,ϕ),(3)

where PJ (θ,ϕ) is the theoretical result [6] derived fromJeffrey’s equation [1]. The viscosity, therefore, tends to beconstant again in this regime, i.e., η = η∞ .
The viscosity exhibits strong shear-thinning behavior in theintermediate regime. The PDF is approximately given byh-kobayashi@aist.go.jp
ryoichi@cheme.kyoto-u.ac.jp1539-3755/2011/84(5)/051404(9)(2)

Pγ̇ (θ,ϕ)  PJ (θ,ϕ) + (Dr /γ̇ )P1 (θ,ϕ),051404-1(4)©2011 American Physical Society

HIDEKI KOBAYASHI AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 84, 051404 (2011)

Intrinsic viscosity(a)
1st Newtonian[η] / [ η] 0

1

Shear thinning

2nd Newtonian

10 -1 -2
10

l 3+ l -3.

1

FIG. 1. A schematic illustration of the viscosity transition. (a) A typical behavior of the
intrinsic viscosity [η] as a function of the dimensionless shear rate γ̇ /Dr . [η]0 is the intrinsic
viscosity in the limit γ̇ /Dr → 0. Here P (ϕ) ≡ π2− π cos θP (θ,ϕ) dθ is normalized so that S ≡ π2 2− π2 P (ϕ) dϕ = 1. (b) P (ϕ) in the weak-shear
regime where the rod undergoes random tumbling. (c) P (ϕ) in the strong-shear regime where
the rod undergoes periodic tumbling.

/Dr

Probability
Fluctuation dominant

Shear flow dominant

(b)

(c)

3

2

P(ϕ)

P(ϕ)

2

1

4/π
1

S=1
S=1
-π/2

0
ϕ

π/2

-π/2

π/2

0
ϕ

where P1 represents the leading term of the perturbation
expansion of the thermal effects. It is clearly seen thatthe contribution from the thermal effects decreases as thedimensionless shear rate γ̇ /Dr increases in this regime, which
gives rise to drastic shear-thinning behavior. The solid line
shown in three different flow regimes in Fig. 1 represents aschematic illustration of the viscosity transition based on theabove considerations.
Consistent with the theoretical model of Hinch and Leal
[5,6], the viscosity transition from the first Newtonian tothe shear-thinning regime has already been successfully
reproduced in various numerical studies [11,16]. However,the viscosity transition from the shear-thinning to the second
Newtonian regime has never been successfully reproducedby numerical simulations. For rigid rod dispersions, we
did not find any previous studies that were performed atshear rates high enough to approach the second Newtonian
regime. Several numerical simulations have been conducted
for flexible chain dispersions at high shear rates that are
expected to be in the second Newtonian regime. However,the viscosity transition from shear-thinning to the second
Newtonian behavior has never been correctly reproduced,
not even when the hydrodynamic interactions are taken into
account using the Rotne-Prager-Yamakawa (RPY) tensor [16].
In the present study, we used a different type of approach,called the smoothed profile method (SPM) [17–21] categorized

as a fluid mechanics approach [22], which can accurately take
into account the thermal fluctuations and the hydrodynamic
coupling between bead particles with a finite radius a and a
Newtonian host fluid, based on DNS of particle dispersions.
The viscosity of a rigid rod dispersion has been calculatedusing SPM to reproduce the viscosity transition from shearthinning to the second Newtonian regime and to examine
carefully the validity of the theoretical model proposed byHinch and Leal [5,6].
II. METHODS
A. ModelWe solve the dynamics of a single rigid rod in a Newtonian
solvent using SPM [17–20]. In this method, the boundaries
between solid particles and solvents are replaced with acontinuous interface by assuming a smoothed profile. This
simple modification enables us to calculate the hydrodynamic
interactions both efficiently and accurately without neglecting
many-body interactions. The equation governing a solventwith density ρf and shear viscosity ηf is a modified NavierStokes equation:

051404-2


ρf


∂u
+ (u · ∇)u = ∇ · σ + ρf φfp + fshear∂t(5)

REENTRANT TRANSITION IN THE SHEAR VISCOSITY . . .

PHYSICAL REVIEW E 84, 051404 (2011)

with the incompressible condition ∇ · u = 0, whereσ = −pI + ηf [∇u + (∇u)T ](6)is the Newtonian stress tensor with a solvent viscosity of ηf ,and u(r,t) and p(r,t) are the velocity and pressure fieldsof the solvent, respectively. A smoothed profile function
0  φ(r,t)  1 distinguishes between the fluid and particle
domains, yielding φ = 1 in the particle domain and φ = 0
in the fluid domain. These domains are separated by thininterstitial regions, the thicknesses of which are given by ξ . The
body force φfp is introduced so that the total velocity field u ofthe dispersion satisfies u(r) = (1 − φ)uf (r) + φup (r), where
uf is the fluid velocity and up represents the rigid motions of
the particles. The incompressible condition ∇ · u thus ensures∇φ · (up − uf ) because both uf and up satisfy incompressible
conditions. The gradient of φ is proportional to the surfacenormal vector and has a support on the interfacial domains.
Therefore, the body force φfp introduced to satisfy therigidity of the particles ensures the appropriate impermeability
boundary conditions at the fluid-particle interface, whereas
the nonslip boundary conditions are imposed automatically by
the viscous stress term in the Navier-Stokes equation. Moredetailed explanations and the mathematical expressions for φ
and φfp were also detailed in our previous papers [17,18]. The
external force fshear is introduced to maintain a linear shearwith a shear rate of γ̇ . This force is applied with the obliquecoordinate transformation based on tensor analysis [21,23].
In the present study, we use a bead-spring model consistingof N freely rotating beads in a single rigid rod. The bead
diameter is σ = 2a. The motion of the ith bead is governed by
the following Newton-Euler equations of motion with thermal
fluctuations:
d
d
P
C
V
Vi = FH
Ri = Vi ,i + Fi + Fi + Gi ,dt
dtdIi · i = NHi + Gi ,
dt
Mi

(7)
(8)

where Ri , Vi , and i are the position, translational velocity,and rotational velocity of the beads, respectively. The variables
Mi and Ii are the mass and moment of inertia, and FH
i and
NH
are
the
hydrodynamic
force
and
torqueexerted
by
the
i
and
G
solvent on the beads, respectively [17,18]. GVii are
the random force and torque, respectively, caused by thermal
fluctuations. The temperature of the system is defined suchthat the long-time diffusive motion of dispersed particles
reproduces the correct behavior [19,20]. FPi represents thepotential force attributable to direct interbead interactions.
From our previous study [21], we can evaluate the average
stress tensor of the dispersion σ dis as1σ dis = σ  −(9)drrρφfpVtwith the volume V = Lx Ly Lz , where Li is the system size inthe i direction, · · · denotes averaging over space and time,and · · · t denotes time averaging over the steady state.
We use a bead-spring model as a model of rodlike
objects with a truncated Lennard-Jones potential and a finitely

extensible nonlinear elastic (FENE) potential. The truncatedLennard-Jones interaction is expressed in terms of ULJ :
1
12
6
4 rσij+rij < 2 6 σ− rσij(10)
ULJ (rij ) =1
0
rij > 2 6 σ ,where rij = |Ri − Rj |. The parametercharacterizes thestrength of the interactions, and σ represents the diameterof the beads. Consecutive beads on a chain are connected by a
FENE potential of the form
 

1
r 2
2
,
(11)
UFENE (r) = − kc R0 ln 1 −2
R0
where r = |Ri+1 − Ri |, kc = 30 /σ 2 , and R0 = 1.5σ . Let
us denote FCi as the constraint force acting on the ith
bead attributable to the bond-angle constraints that cause theconnected beads to form a straight rod. Thus,
 N


∂
μ · α ,(12)
FCi =∂Ri α=3 α α = (α − 2)R1 − (α − 1)R2 + Rα ,(13)

where  α = 0 is the constraint condition to be satisfied.
The variable μα is a Lagrange multiplier associated with theconstraints chosen such that the condition  α = 0 is satisfied
at a time t + h, where h is the time increment of a single
simulation step.
The numerical simulations are performed in three dimensions with periodic boundary conditions. The lattice spacing is taken to be the unit of length. The unit of time isgiven by ρf 2 /ηf , where ηf = 1 and ρf = 1. The systemsize is Lx × Ly × Lz = 32 × 32 × 32. The other parameters
include the following: σ = 4, ξ = 2, = 1, Mi = 4π a 3 /3,
N = 5, and h = 6.7 × 10−2 . In the presented simulations
under shear flow, the Navier-Stokes equation is discretizedwith a dealiased Fourier spectral scheme in space and witha Euler scheme in time [21]. To follow the motions of
the beads, the positions, velocities and angular velocities ofthe beads are integrated with the Adams-Bashforth scheme.
The bead particles are assumed to be neutrally buoyant, andhence no gravity effects are considered. At t = 0, the rigid rodaligns along the x axis, which is the flow direction. The total
duration τt of each simulation is set such that γ̇ τt  3500.
The range of kB T is 5.0 × 10−4 < kB T < 32 and that of γ̇
is 5.0 × 10−3 < γ̇ < 2.0 × 10−2 , where kB is the Boltzmann
constant, and T is the temperature. From the symmetry of thesystem, we follow the polar angles θ and ϕ defined in Fig. 2 toconsider the motion of a rigid rod. The angle defined between
the rod and the x-y plane is denoted by θ , and the angle defined
between the rod projected on the x-y plane and the x axis isdenoted by ϕ.
B. Effective aspect ratio

In the present study, the rigid rod is represented asconnected beads. Because the beads composing the rod canrotate freely, the effective aspect ratio l differs from the simple
geometrical aspect ratio L/σ , where L  N σ is the rod’s
length. Instead, we evaluate l numerically with the PDF of

051404-3

HIDEKI KOBAYASHI AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 84, 051404 (2011)
2.5

2

PJ(ϕ)

1.5

1

0.5

0
-π/2

-π/4

0

FIG. 2. The geometry of the rod’s orientation in the present
simulations.

the rotating rigid rods without thermal fluctuations in the x-yplane, i.e., θ = 0, as represented by
PJ (ϕ) = l 2 −1
l 2 +1

C0
1
sin2 ϕ + l 2 +1

(14)

,

where C0 is determined from the normalization condition π2− π2 PJ (ϕ) dϕ = 1 [1]. Equation (14) was obtained in thefollowing manner. The projection of the PDF PJ (ϕ) of arotating rigid rod on the x-y plane is governed by a FokkerPlanck equation of the form
∂PJ (ϕ)∂PJ (ϕ)∂[ωPJ (ϕ)]
,
=
+ 2Dr
∂t
∂ϕ
∂ϕ 2

(15)

where ω = ϕ̇ is the angular velocity of the tumbling rod.
When the rigid rod rotates in the x-y plane without thermal
fluctuations in steady state, the Fokker-Planck equation ismodified to∂[ωPJ (ϕ)]
= 0.
∂ϕ

(16)

In this case, ω is represented as
ω = γ̇

l2 − 1 2
1
sin ϕ + 2
l2 + 1
l +1


(17)

from Jeffrey’s equation [1]. Equation (14) is obtained because
PJ (ϕ) is inversely proportional
 to ω. Figure 3 shows that ournumerical results of P (ϕ) = cos θ Pγ̇ (θ,ϕ) dθ of the strongshear regime agree well with PJ (ϕ) with l = 7.1. We thus use
l = 7.1 for the present rigid rod, which is composed of freely
rotating beads.
C. Analytic formula for the viscosity

Hinch and Leal [5,6] studied the rheological properties of adilute dispersion of rigid nonspherical particles in steady shear
flow. They obtained an analytical formula for the dispersion
disviscosity η ≡ σ12/γ̇ . We analyze our numerical results using
their formula. The dispersion viscosity η is given by theensemble average of the temporal viscosity η̂(θ,ϕ) using the

π/4

π/2

ϕ

FIG. 3. The PDF for a rotational rigid rod as a function of ϕwithout thermal fluctuations. Numerical results (open circle) andEq. (14) with l = 7.1 (solid line).

PDF of the two angles for the rotating rigid rods Pγ̇ (θ,ϕ),which satisfies the normalization condition,
 π
 π
2
2
cos θ dθdϕPγ̇ (θ,ϕ) = 1.
(18)
− π2

− π2

The temporal shear viscosity is found to be
2η̂(θ,ϕ) = ηf 1 +  A cos4 θ sin2 2ϕ + 2B cos2 θ +
I3

Dr 1
+ F cos2 θ sin 2ϕ ,
(19)
γ̇ 2
where  is the volume fraction of suspended particles, Dr is
the rotational diffusion constant, and A,B,F,I3 are the shapefunctions given in previous studies [4–6]. In the case of rigid
rods, A, B, F , and I3 are dependent only on the aspect ratio l
at l = 7.1, A = 8.44, B = 0.06, and I3 = 0.99.
The shear viscosity of the dispersion is obtained bysubstituting Eq. (19) into Eq. (1). When we consider thestrong-shear case Dr  γ̇ , we can safely neglect the final termin Eq. (19). The dynamics of the angle ϕ become decoupledfrom the angle θ because the angle θ is sufficiently small for alarge l [6]. Thus, we obtain the following formula:η(γ̇ ) = η̂(θ,ϕ)P (θ )P (ϕ) dθ dϕ(20)= ηf 1 +  Acos4 θ θ sin2 2ϕ ϕ + 2Bcos2 θ θ
2+ + E ,(21)I3πg(ϕ) ϕ =where
f (θ ) θ = −2 π cos θ dθf (θ )P (θ ),
2
 π2
 π2
− π2 dϕg(ϕ)P (ϕ), P (θ ) ≡ − π2 Pγ̇ (θ,ϕ) dϕ, and E isthe error arising from the separation of integrals over θ and ϕ.
We can neglect E safely because E is sufficiently smallin comparison with the other terms.

051404-4

REENTRANT TRANSITION IN THE SHEAR VISCOSITY . . .

PHYSICAL REVIEW E 84, 051404 (2011)0τp /122τp /123τp /12

FIG. 4. (Color online) The time evolution of the fluid velocity field u at kB T = 0. The times range from 0 to τp /4, where τp is the rotational
period. The velocity fields are plotted only when |ud |/|u| > 0.05, where ud is the deviation velocity field u − γ̇ (y − Ly /2)ex .

III. RESULTS

number that represents the strength of the shear flow normalized by the strength of diffusion caused by thermal fluctuations.
In our work, Pe is defined as

Figure 4 shows the time evolution of the fluid velocity field ufor kB T = 0. Display of the velocity field, depicted by arrows,
is limited to regions where the quantity |ud |/|u| exceeds 0.05,
where ud is the deviation velocity field defined as u − γ̇ (y −Ly /2)ex . The time interval displayed is from 0 to τp /4, wheret = 0 is the time at which ϕ = 0 and τp is the rotational period.
The rigid rod aligns almost fully with the shear direction inmost of period. This motion is in good agreement with Jeffrey’s
equation [1]. The velocity u is approximately equal to γ̇ (y −Ly /2)ex except in the vicinity of the rod. To examine thedeviation velocity, we plot the time evolution of the deviation
velocity field ud in Fig. 5 in regions where |ud |/|u| > 0.05.
The amplitude of ud is significant only near the rod. This isbecause the deviation velocity is attributable to an exchange ofmomentum between the particle and fluid. The amplitude of
ud at time 2τp /12 is larger than at t = 0,τp /4. This behavior isin qualitative agreement with that of η̂(θ,ϕ) which, as Eq. (19)reveals, has minima at ϕ = 0,π/2 and a maximum at ϕ = π/4.
The momentum exchange between particle and fluid is clearly
evident in η̂(θ,ϕ).
In Fig. 6 we plotted the intrinsic viscosityη − ηf
(22)
[η] ≡
ηf 

We find that the intrinsic viscosity [η] gradually changes
from non-Newtonian (shear-thinning) to Newtonian behavior
with increasing Peclet number, as shown in Fig. 6. The
present simulation data for [η] show shear-thinning behavior
for Pe < 102 and second Newtonian behavior for 104 < Pe.
These results are in good agreement with previous theoreticalstudies [5,6,9].
To quantitatively compare our results with those of Hinchand Leal [5,6], we obtain the relationship between Pe in our
definition and γ̇ /Dr , which is used in Hinch and Leal’s work
[5] instead of Pe. On the basis of the shell model [24,25], therotational diffusion constant Dr for a rigid rod is calculated asof the dispersion obtained from the present simulations as afunction of the Peclet number Pe. Pe is the dimensionlessIn the shell model mentioned above, the contour of themacromolecules of arbitrary shape is represented by a shell

Pe =

Dr =

6π ηf σ 3 γ̇.
kB T

3[ln l + d(l)]kB T,
π η f L3d(l) = −0.662 +

051404-5

0.917 0.05
− 2 .
l
l

(23)

(24)
(25)

HIDEKI KOBAYASHI AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 84, 051404 (2011)0τp /122τp /123τp /12

FIG. 5. (Color online) The time evolution of the deviation velocity field is ud ≡ u − γ̇ (y − Ly /2)ex at kB T = 0. The times range from 0
to τp /4, where τp is the rotational period. The velocity fields are plotted only when |ud |/|u| > 0.05.

composed of many identical small beads. The shell model canbe adequately modeled by decreasing the size of the beads.

20

R1

R2

R3

7.35x10-2

2.64x101

γ̇
l3
Pe.
=
Dr
18[ln l + d(l)]

R4
104

10
9
8
7
6
5
4

l2
4
,
15 ln l + ϒ(l)1.38 8.87 8.82
+ 2 − 3 .
ϒ(l) = −0.90 −
l
l
l

3

2
10-2

[η] =

10-1

100

101

102

103

104

(26)

The theoretical model of Hinch and Leal is also plotted in Fig. 6
with the solid lines in the three different regimes, namely, theweak (R 1 ), intermediate (R 2 ), and strong (R 3 + R 4 ) shear
regimes.
According to the work of Hinch and Leal [5], for the weakshear regime γ̇ /Dr  1, namely, Pe  7.35 × 10−2 , which isdenoted by R 1 in Fig. 6, [η] is constant. On the basis of Ortega’s
work [26], the intrinsic viscosity [η] of the weak-shear flow
regime for a rigid rod with a short aspect ratio is calculated as

14

[η]

From Eqs. (23) and (24), the relationship between γ̇ /Dr and
Pe is expressed as

105

Pe

FIG. 6. The intrinsic viscosity as a function of Pe. [η] (circle)
and [ηθ ] (square). The three solid lines correspond to the theoreticalresults of Hinch and Leal [5,6], with the first Newtonian regime
denoted by R 1 , the shear-thinning regime denoted by R 2 , and thesecond Newtonian regime denoted by R 3 + R 4 . In our simulation, theviscosity shows an undershoot before reaching the second Newtonian
regime R 3 .

(27)
(28)

This expression is identical to Hinch and Leal’s result in thelimit of l → ∞. The lowest shear rate that we consider inthe present simulations is still not in the weak-shear regime
because of the extremely long simulation time needed to obtainreliable data.
For the intermediate shear regime 1  γ̇ /Dr  l 3 + l −3 ,namely, 7.35 × 10−2  Pe  26.4, which is denoted by R 2

051404-6

REENTRANT TRANSITION IN THE SHEAR VISCOSITY . . .

PHYSICAL REVIEW E 84, 051404 (2011)

in Fig. 6, the intrinsic viscosity [η] shows shear-thinning asderived from Eq. (21),12[η] = C1 Pe−1/3 + 2B + 2/I3 ,

10

(29)

− π2

We obtained [η∞ ] = 3.82 from our numerical data in the highshear regime, which is denoted by R 4 in Fig. 6. The error fromthe theoretical value 3.97 is within 4.26%.
For the regime 102 < Pe < 103 , which is denoted by R 3 inFig. 6, the behavior of [η] shows a notable undershoot before
reaching the high-shear limiting second Newtonian viscosity.
This effect is attributable to the fluctuations in θ and gives
rise to the deviations of P (θ ) from its high-shear limitingπform P ∗ (θ ) ≡ −2 π PJ (θ,ϕ) dϕ and increases with decreasing
2
Pe. The rods tend to align in the flow direction with increasing
Pe. Therefore, cos4 θ and cos2 θ monotonically increase up
to their high-shear limiting values with increasing shear rate.
Increasing cos4 θ and cos2 θ leads to an increase in [η]up to [η∞ ] through Eq. (21). To examine the role of thermal
fluctuations in θ in more detail, let us define[ηθ ] = Asin2 ϕ J cos4 θ θ + 2Bcos2 θ θ +2
I3= 1.81cos4 θ θ + 0.12cos2 θ θ + 2.02

(31)

8

[Δη]

where C1 is an arbitrary constant. Figure 6 shows good
agreement between the data from the present simulation data
and those of Hinch and Leal, where C1 was determined to fit thesimulation data. When l is sufficiently large, the contributions
from the final two terms in Eq. (29) become negligible, and[η] ∝ Pe−1/3 .
For the strong-shear regime l 3 + l −3  γ̇ /Dr , namely,
26.4  Pe, which is denoted by R3 and R4 in Fig. 6, thetheory predicts that the intrinsic viscosity [η] is constantat Asin2 2ϕ J + 2B + 2/I3 = 3.99 = [η∞ ] from Eq. (21).
In this case, · · · J denotes the ensemble average, which is
calculated as π
2
f (ϕ) J =dϕf (ϕ)PJ (ϕ).
(30)

= Acos4 θ θ (sin2 2ϕ ϕ − sin2 2ϕ J )

2
0
10-1

100

101

102
Pe

103

104

105

FIG. 7. The behavior of [η] as a function of Pe. [η] goes tozero at approximately Pe ≈ 150.

predicted by Hinch and Leal [5] for the viscosity transition
from shear-thinning to the second Newtonian behavior but
agrees well with our theoretical prediction of Pec = 143, at
which the dynamic crossover from Brownian to non-Brownian
behavior occurs in the rotational motion of the rotating rod at
l = 7.1 [27].
dis
disFigure 8 shows the first normal stress, N1 = (σ11− σ33),normalized with respect to the viscous stress of the solvent,ηγ̇ , as a function of P e. According to Hinch and Leal [5],in the intermediate shear regime, denoted by R2 in Fig. 8,
1/3
the first normal stress behaves as N1 ∝ γ̇ 2/3 Dr . Therefore,
1/3
N1 /ηγ̇ ∝ γ̇ −1/3 Dr ∝ Pe−1/3 in regime R2 . This behavior isin good agreement with our numerical results.
In the strong-shear regime, denoted by R3 and R4 in Fig. 8,
the theory predicts a first normal stress of N1 = 0. This is as
given by Jeffrey’s equation [1], which perfectly describes the

(32)

10-1

R3

R2

R4

γ = 0.005
γ = 0.01
0.0008 (Pe=infinity,γ = 0.005)
0.0017 (Pe=infinity,γ = 0.01)
Pe-1/3

N1/ηγ

10-2

-3

10

-4

10

(33)
(34)

to eliminate the contribution of θ fluctuations from the total
intrinsic viscosity of the dispersion. Figure 7 shows the behavior of [η] as a function of Pe. [η] decreases with increasing
Pe, and finally [η] goes to zero at approximately Pe ≈ 150.
This value is considerably different from the value Pe = 26.4

Pec=143

4

to estimate the contribution of θ fluctuations on the total
intrinsic viscosity of the dispersion. In this case, P (ϕ) =PJ (ϕ) is assumed in Eq. (21), and cos4 θ θ and cos2 θ θare evaluated numerically from the present simulations. The
results are plotted in Fig. 6 with the square symbols. One cansee that the data of [ηθ ] almost perfectly collapse onto thoseof [η] for 102 < Pe.
However, the shear-thinning behavior observed for 102 >Pe is attributable to the effect of the thermal fluctuations in
ϕ. This behavior gives rise to the deviations of P (ϕ) from its
high-shear limiting form PJ (ϕ) and increases with decreasing
Pe. To examine this effect quantitatively, we introduce[η] ≡ [η] − [ηθ ]

6

100

101

102

103

104

105

Pe

FIG. 8. Normalized first normal stress N1 /ηγ̇ as a function of Pe.
γ̇ = 0.01 (closed circle) and γ̇ = 0.005 (open circle). The solid line
corresponds to Pe−1/3 , which is the theoretical result of Hinch andLeal [5,6]. The dashed (dotted) line corresponds to N1 /ηγ̇ in the limitthat Pe → ∞ for γ̇ = 0.01 (γ̇ = 0.005).

051404-7

HIDEKI KOBAYASHI AND RYOICHI YAMAMOTO

PHYSICAL REVIEW E 84, 051404 (2011)
1〈 sin22ϕ 〉∞tumbling motion in this regime. However, even if we assume
kB T = 0 so that Pe → ∞, the numerical error in the firstnormal stress E(γ̇ ) remains finite and depends on the shear
rate. Nevertheless, because of this error, our results departfrom the theoretical prediction of Hinch and Leal [5] for highPe. We consider E(γ̇ ) to be a numerical error arising from thecompetition between shear force and interbeads interaction, as
described by Eqs. (10) and (11), when the rigid rod rotates inJeffrey’s orbit. Figure 8 shows E(γ̇ ) for γ̇ = 0.01,0.005, when
Pe > 104 , E(γ̇ ) is in good agreement with N1 /ηγ̇ at each shear
rate. This agreement is because the tumbling motion is welldescribed by Jeffrey’s equation [1].
Also of interest are the off-diagonal entries of the stresstensor and the second normal stress. However, these values
are not accurately calculated by this numerical method. In
comparison with theoretical results, the numerical error in thedisstress tensor σ12
is near 4%. The other entries of the stressdis. For example, the second
tensor are less than 0.01 × σ12
dis
dis
dis
normal stress, N2 ≡ σ22 − σ33 ≈ 0.01 × σ12, is buried innumerical error.

0.1

0.01
100

(39)

It is demonstrated in Fig. 9 that the first term in Eq. (38) showssin2 2ϕ ∞ ∝ (γ̇ /Dr )−1/3 for the entire range of Pe. Figure 10
shows that only A is increasing with increasing l, while B
and 2/I3 tend to be decreasing or constant with increasing l.
We estimate A/(2B + 2/I3 ) ∼ l 1.8 for l → ∞. This relational
expression indicates that[η] ∝ sin2 2ϕ ∞ ∝ (γ̇ /Dr )−1/3

(35)

(38)

− π2



γ̇
dψ exp −f (ψ,ϕ) ,P (ϕ) = C1
4Dr
0

2
f (ψ,ϕ) = ψ − 1 − 2sin ψ cos(ψ − 2ϕ),
l +1

106

is rewritten as[η] = Asin2 2ϕ ∞ + 2B + 2/I3 , π2dϕf (ϕ)P∞ (ϕ).
f (ϕ) ∞ =Let us discuss the numerical models with which theviscosity transition to the second Newtonian regime occurs
based on Eq. (21). For the strong-shear regime 1  γ̇ /Dr ,we can estimate cos2 θ θ  1 and cos4 θ θ  1 becausePγ̇ (θ,ϕ)  P ∗ (θ )P (ϕ). In this case P (ϕ) satisfies the FokkerPlanck equation, shown as Eq. (15), for which the formal
solution is given by


102
103
104
105
Dimensionless shear rate

FIG. 9. The behavior of sin2 2ϕ ∞ as a function of γ̇ /Dr .
Numerical results (circle). The solid line corresponds to (γ̇ /Dr )−1/3

IV. DISCUSSION

 π

101

(40)

holds for the entire range of Pe without indicating the occurrence of second Newtonian behavior. The same conclusion(36)lwhere C1 is determined from the normalization condition, π2− π2 P (ϕ) dϕ = 1. When γ̇ /Dr is sufficiently large, P (ϕ)
converges to PJ (ϕ), represented by Eq. (14), and the viscosity
displays second Newtonian behavior.
The above discussion is not valid for the limit of l →∞, which corresponds to an infinitely long or equivalently
infinitely thin rod. In this limit the angular velocity of the
tumbling rod becomes zero at ϕ = 0 from Eq. (17). Thus,
the rod cannot continue its rotational motion without thermal
fluctuations because the hydrodynamic torque acting on therod becomes zero at ϕ = 0 for l → ∞. Therefore, P (ϕ) in
Eq. (35) is modified toP∞ (ϕ) = C2γ̇dψ exp −[ψ − sin ψ cos(ψ − 2ϕ)] ,
4Dr
0
(37)

 π

where
C2 is determined from the normalization condition, π2− π P∞ (ϕ) dϕ = 1. Using Eq. (37), the intrinsic viscosity [η]2lFIG. 10. The behaviors of A, B, and 2/I3 as a function of the
aspect ratio l. A (bold solid line), B (dashed line), 2/I3 (dotted line).
The thin solid line represents l 1.8 .

051404-8

REENTRANT TRANSITION IN THE SHEAR VISCOSITY . . .

PHYSICAL REVIEW E 84, 051404 (2011)

can be derived by considering a characteristic shear rate γ̇ ∗at which the first term in Eq. (38) becomes comparable to theremaining terms. The condition is satisfied atγ̇ ∗∼ l 5.4 .
Dr(41)Equation (41) indicates γ̇ ∗ → ∞ for l → ∞.
In the case of a previous numerical study [16], thehydrodynamic force acting on each bead particle was considered through the RPY tensor. Although the translationalhydrodynamic force was properly considered, the rotationalhydrodynamic torque acting on each bead particle was completely ignored in that study. Therefore, it is suspected thatthe hydrodynamic torque acting on the rod becomes zero at
ϕ = 0. The rod cannot continue rotational motion at a highshear rate, where the effect of thermal fluctuations disappears.
This situation is identical to the case of l → ∞. We expect thatthe second Newtonian regime could be correctly reproducedwith the RPY tensor approach if the hydrodynamic torque isproperly taken into account.
V. CONCLUSION

In the present study, we numerically calculated the intrinsic
viscosity [η] of a dilute dispersion of rigid rods using a DNS
method known as SPM. Simulations were conducted under the
influence of thermal fluctuations and shear flow in the ranges of
5.0 × 10−4 < kB T < 32 and 5.0 × 10−3 < γ̇ < 2.0 × 10−2 ,respectively. We have successfully reproduced the viscosity transition from shear-thinning to the second Newtonian[1] G. B. Jeffrey, Proc. R. Soc. A 102, 161 (1922).
[2] A. Peterlin and H. A. Stuart, Z. Phys. 112, 1 (1939).
[3] V. H. Giesekus, Rheol. Acta. 2, 50 (1962).
[4] G. K. Batchelor, J. Fluid Mech. 41, 545 (1970).
[5] E. J. Hinch and L. G. Leal, J. Fluid Mech. 52, 683 (1972).
[6] L. G. Leal and E. J. Hinch, J. Fluid Mech. 46, 685 (1971).
[7] E. J. Hinch and L. G. Leal, J. Fluid Mech. 71, 481 (1975).
[8] E. J. Hinch and L. G. Leal, J. Fluid Mech. 76, 187 (1976).
[9] L. G. Leal and E. J. Hinch, Rheol. Acta 12, 127 (1973).
[10] H. A. Scheraga, J. Chem. Phys. 23, 1526 (1955).
[11] J. F. Ryder and J. M. Yeomans, J. Chem. Phys. 125, 194906(2006).
[12] H. Chen, Y. Ding, and A. Lapkin, Powder Tech. 194, 132(2009).
[13] G. Chauveteau, J. Rheol. 26, 111 (1982).
[14] I. Noda, Y. Yamada, and M. Nagasawa, J. Phys. Chem. 72, 2890
(1968).
[15] R. L. Christiansen and R. B. Bird, J. Non-Newtonian Fluid Mech.
3, 161 (1977).

regime, as was correctly predicted by the theoretical model ofHinch and Leal [5,6].
However, some discrepancies were noted between the theoretical predictions and the results of the present simulations.
By defining [η] to eliminate the effects of fluctuations in θ ,which is not considered in the theoretical model, we confirmedthat the viscosity transition from shear-thinning to the second
Newtonian behavior occurs at approximately Pe = 150. Thisvalue is considerably larger than the value of 26.4 predicted
by Hinch and Leal [5] but agrees well with our theoreticalprediction of Pec = 143, at which the dynamic crossover from
Brownian to non-Brownian behavior occurs in the rotationalmotion of the rotating rod [27].
We have analyzed the mechanism of the viscosity undershoot observed in our simulation before reaching the second
Newtonian regime. Shear flow suppresses fluctuations in ϕ andθ as its rate is increased. The former contributes to decreasing[η], but the latter contributes to increasing [η]. The undershootoccurs because of the two competing effects.
We also conclude that the viscosity transition to thesecond Newtonian regime can be reproduced correctly only
if the hydrodynamic torque is properly taken into account innumerical models of the dispersions.
ACKNOWLEDGMENTS

The authors would like to express their gratitude to
Dr. T. Murashima and Dr. Y. Nakayama for useful comments
and discussions.

[16] D. Petera and M. Muthukumar, J. Chem. Phys. 111, 7614 (1999).
[17] Y. Nakayama and R. Yamamoto, Phys. Rev. E 71, 036707(2005).
[18] Y. Nakayama, K. Kim, and R. Yamamoto, Eur. Phys. J. E 26,
361 (2008).
[19] T. Iwashita, Y. Nakayama, and R. Yamamoto, J. Phys. Soc. Jpn.
77, 074007 (2008).
[20] T. Iwashita and R. Yamamoto, Phys. Rev. E 79, 031401 (2009).
[21] H. Kobayashi and R. Yamamoto, J. Chem. Phys. 134, 064110(2011).
[22] I. Pagonabarraga, B. Rotenberg, and D. Frenkel, Phys. Chem.
Chem. Phys. 12, 9566 (2010).
[23] A. Onuki, J. Phys. Soc. Jpn. 66, 1836 (1997).
[24] M. M. Tirado and J. G. de la Torre, J. Chem. Phys. 73, 1986
(1980).
[25] M. M. Tirado, J. Chem. Phys. 81, 2047 (1984).
[26] A. Ortega and J. G. de la Torre, J. Chem. Phys. 119, 9914 (2003).
[27] H. Kobayashi and R. Yamamoto, Phys. Rev. E 81, 041807(2010).

051404-9