Microstructure and Rheology of Rigid Rod Suspensions - Industrial

Article pubs.acs.org/IECR

Microstructure and Rheology of Rigid Rod Suspensions YoChan Kim,† André Bénard,‡ and Charles A. Petty*,† †

Department of Chemical Engineering and Materials Science, and ‡Department of Mechanical Engineering, Michigan State University, East Lansing, Michigan 48824, United States ABSTRACT: Nonspherical particles dispersed in a fluid have a tendency to align because of particle−fluid interactions. At high particle concentrations and in the absence of fluid dynamic couples, spontaneous self-alignment can occur due to excluded volume constraints on rotary Brownian motion. This phenomenon prevents a return-to-isotropy from an anisotropic state. In this communication, the low-order moments of the rotary Smoluchowski (S) equation for a rigid rod suspension are used to explore the effect of the Péclet number on the shear viscosity of a suspension. A closure model for the fourth-order orientation tetradic uses an algebraic fully symmetric quadratic (FSQ) mapping of the second-order orientation moment into the fourth-order orientation moment. The algebraic mapping preserves the 6-fold symmetry and the 6-fold contraction properties of the exact fourth-order orientation moment. The theory predicts that, if the orientation director is in the tumbling regime, shear thickening may occur and, if the orientation director is in the wagging regime, shear thinning may occur.

1. INTRODUCTION Flow induced alignment of rigid rod suspensions provides a practical means to understand the microstructure and rheology of structured fluids. The microstructure may be either isotropic or anisotropic. Under extreme conditions, a nematic phase may occur wherein rod-like particles spontaneously align. The tendency for some suspensions to develop nematic-like microstructures has a significant and practical impact on the rheological, optical, and material properties of structured fluids (see the works of Doi;1 Doi and Edwards;2,3 and Larson4). In sections 2−4 below, a closed equation-of-change for the second-order orientation moment is developed by using an algebraic preclosure for the fourth-order orientation moment in terms of the second-order orientation moment. In section 5, the theory is used to illustrate how the rotary Péclet number influences the shear viscosity. Shear thickening may occur if the following three conditions hold: (1) the velocity gradient is large compared with the rotary diffusion coefficient (i.e., ∥∇u∥ ≫ D°R); (2) the excluded volume potential is large compared with the local energetic state of the suspension (i.e., ΔÛ ≫ [ckB*T]); and (3) the director of the microstructure tumbles. At high Péclet numbers, a Newtonian plateau occurs as the microstructure transitions from a director tumbling regime to a steady alignment regime. For a comprehensive development of the theory employed in this communication, see the works of Kim5 and Petty et al.6

local velocity of the suspension (see p 50 in the work of Doi and Edwards3): ⎡⎛ ⎞ ⎤ ⎡ ∂Ψ ⎤ ∂ ⎢⎣ ̂ ⎥⎦ = −⎢⎜⎜ ⎟⎟ ·( p̲ ̇Ψ)⎥ ⎢⎣⎝ ∂ p̲ ⎠ ⎥⎦ ∂t X̲ ̂ t̂ 2π

∫0 ∫0

̂ ϕ(t )), ̂ χ̲ (X̲ ̂ , t )] ̂ sin(θ ) dθ dϕ = 1 Ψ[ p̲ (θ(t ), (3)

The vector X̲ ̂ represents the spatial position of a material f luid element at some fixed time (i.e., t ̂ = 0); the spatial position of a material f luid element at a later time is x̲ ̂ = χ(̲ X̲ ̂ , t )̂ . The vector field χ(̲ X̲ ̂ , t )̂ represents the motion of the suspension. The differential operator on the left-hand-side of eq 2 represents the material (or substantial) derivative of the orientation density function. The instantaneous rotary flux of orientation states relative to a material frame-of-reference can be separated into a rotary convective f lux and a rotary drif t (dif f usive) f lux: [ p̲ ̇Ψ] ≡ [ p̲ ̇ Ψ] + [( p̲ ̇ − p̲ ̇ )Ψ] c c

(4)

The rigid rods have the same density as the suspending fluid, so gravity is unimportant. Equation 2 neglects spatial diffusion of the suspension relative to the mass average velocity. A balance of angular momentum on an axisymmetric rigid rod yields the following equation-of-change for the instantaneous angular velocity (see the works of Jeffery;7 Batchelor;8,9 Bibbo et al.;10 and Petty et al.6):

2. ROTARY SMOLUCHOWSKI EQUATION FOR RIGID ROD SUSPENSIONS The instantaneous orientation vector of a rigid rod in a suspension can be represented in terms of Euclidean base vectors and the time-dependent Euler angles θ and ϕ:

Special Issue: Scott Fogler Festschrift

p̲ = sin(θ ) cos(ϕ) e̲ x + sin(θ ) cos(ϕ) e̲ y + cos(θ ) e̲ z

Received: October 12, 2014 Revised: January 9, 2015 Accepted: January 20, 2015

(1)

The following rotary Smoluchowski (RS) equation governs the changes in the rotary orientation density function relative to the © XXXX American Chemical Society

π

(2)

A

DOI: 10.1021/ie503995y Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research ⎡ ∂ p̲ ⎤ [ p̲ ̇ ] ≡ ⎢ c ⎥ = −[W̲ ̂ · p̲ ] + [λ T][ I̲ − p̲ p̲ ]·[ S̲ ·̂ p̲ ] c ⎢⎣ ∂t ̂ ⎥⎦ X̲ ̂

Doi

̃ (IIb, IIIb)DRo , ⟨DR ⟩ ≡ FTD (5)

(11)

2

(6)

o ⎧ dilute ⎪C TD , ̃ Doi = ⎨ FTD ⎪ o 2 ⎩C TD/[1 − 3IIb/2] , semidilute

For needlelike particles, [λT] = +1; for disklike particles, [λT] = −1. A typical tumbling parameter for slender rodlike particles is about +0.7, which corresponds to L/d ≅ 2.38 (see p 280 in the work of Larson4). L is a length scale associated with a rigid rod, and d is its diameter. The rotational period of a rigid rod depends on the tumbling parameter. For [λT] = ±1, a prolate spheroidal particle and an oblate particle have infinite rotation periods (i.e., they are not rotating). Experimental evidence for particle tumbling in rigid rod suspensions has been reviewed by Anczurowski and Mason,13 Frattini and Fuller,14 and Larson4 (see, especially, p 280). In eq 5, the operator S̲ ̂ represents the rate-of-strain tensor (dyadic) S̲ ̂ =

1 ̂ ||∇ u̲ || [∇ u̲ ̂ + (∇̂ u̲ )̂ T ] = [∇ u̲ + (∇ u̲ )T ] 2 2

(7)

ΔÛ = −U[ckB*T ] p̲ p̲ : ⟨ p̲ p̲ ⟩

duz ≡ γ ̇, constant dy

(8)

⎛ ∂Ψ ⎞ ∂ ⎜ ⎟ = − Pe ·[( −W̲ · p̲ + λ T[ I̲ − p̲ p̲ ]·[ S̲ · p̲ ])Ψ] ⎝ ∂t ⎠ X̲ ∂ p̲

(9)

The rotary drift flux in eq 4 can be approximated by (see the works of Doi and Edwards;3 Larson;4 and Kuzuu and Doi15):

Doi + FTD

⎛ ∂Ψ ∂ ⎛ ΔÛ ⎞⎞ ( p̲ ̇ − p̲ ̇ )Ψ = −⟨DR ⟩⎜⎜ +Ψ ⎜ ⎟⎟⎟ c ∂ p̲ ⎝ ckB*T ⎠⎠ ⎝ ∂ p̲ ΔÛ [ = ]

(thermal energy) (suspension volume)

⎞ ∂ ⎛ ∂Ψ ∂ ·⎜⎜ − U Ψ ( p̲ p̲ : ⟨ p̲ p̲ ⟩)⎟⎟ ∂ p̲ ⎝ ∂ p̲ ∂ p̲ ⎠

(14)

In the above equation, t ≡ CoTDDoRt;̂ Pe ≡ γ̇/(CoTDDoR); and, FDoi TD o ≡ F̃ Doi TD /CTD.

(excluded volume energy) (suspension volume)

ckB*T[ = ]

(13)

In the above equation, U is a positive dimensionless, scalarvalued coefficient that depends on the volume fraction of particles, the excluded volume (i.e.,VE ≡ πao2L), and the volume of a single rod (i.e.,VR ≡ πd2L/4). The parameter ao (2ao > d) represents the radius of a hypothetical tube of length L that contains a single rigid rod (see p 66 in the work of Larson4). With Jeffery’s model for rotary convection and Doi’s model for rotary diffusion, the RS-equation can be written as

For homogeneous shear, the velocity field is steady and has only one component, u̲ = uz(y)e̲z, and the shear rate also has only one component: ∇̂ u̲ ̂ = γ ̇ e̲ y e̲ z ,

(12)

The excluded volume potential ΔÛ in eq 10 accounts for the interaction of a rigid rod particle with neighboring particles. This phenomenon has important consequences that partly explain the self-alignment and the flow-induced alignment of rigid rod suspensions. The excluded volume potential supports the physical idea that particles cannot occupy the same place at the same time. An instantaneous model for the excluded volume potential has been developed by minimizing the free energy for rigid rod suspensions with the result that (see p 359 in the work of Doi and Edwards;3 Doi;1 Ilg et al.;16 and Onsager17)

and the operator W̲ ̂ represents the vorticity tensor (dyadic) 1 ||∇ u̲ || W̲ ̂ = [∇̂ u̲ ̂ − (∇̂ u̲ )̂ T ] = [∇ u̲ − (∇ u̲ )T ] 2 2

3ckB*T 1 [=] πηS time

In the above equation, ck*B and ηS represent, respectively, the Boltzmann constant for the suspension and the solvent viscosity. The volume fraction of the dispersed phase is represented by c, and T is the absolute temperature. DoR is the rotary diffusion coefficient for a dilute suspension of rigid rods (see p 334 in the work of Doi and Edwards3 and p 281 in the work of Larson4). The Doi tube dilation factor, F̃ Doi TD (IIb, IIIb), accounts for the influence of the local microstructure on rotary diffusion (see section 3 below). The Doi model assumes that the third invariant IIIb of the orientation dyadic does not influence tube dilation:

Note that the instantaneous orientation vector and the instantaneous angular velocity of the rigid rod are mutually orthogonal. In the above equation, [λT] is a dimensionless scalar-valued tumbling coef f icient that depends on the aspect ratio (i.e., L/d) of the rigid rod (see the work of Trevelyn and Mason;11 Bretherton;12 and Larson4): L −1 ( d) −1 ≤ λ T ≡ ≤ +1 2 ( Ld ) + 1

DRo ≡

3. MICROSTRUCTURE The orientation density function satisfies the symmetry condition Ψ(p̲, t) = Ψ(−p̲, t) inasmuch as there is no distinction between the head and the tail of a rod. Thus, for a suspension of axisymmetric particles, the first moment, ⟨p⟩, ̲ and all higher-order odd moments of the orientation density function are null. The second moment of the orientation density function (orientation dyadic) is defined as

(10)

In the above equation, ΔÛ is a specific excluded volume potential and ⟨DR⟩ represents a pre-averaged rotary diffusion coefficient. In general, ⟨DR⟩ depends on the particle aspect ratio, the volume fraction of particles, the temperature, the pressure, and the invariants of the microstructure (p 360 in the work of Doi and Edwards3 and p 520 in the work of Larson4):

2π

⟨ p̲ p̲ ⟩ ≡ B

∫0 ∫0

π

p̲ p̲ Ψ̲(θ , ϕ , t )sin(θ ) dθ dϕ

(15)


Article

Industrial & Engineering Chemistry Research The orientation dyadic ⟨p̲ p⟩̲ is a symmetric, non-negative operator, and its trace is unity tr(⟨ p̲ p̲ ⟩) ≡ ⟨ p̲ · p̲ ⟩ = 1

(16)

The anisotropic orientation dyadic (or, equivalently, the structure tensor) is defined as b̲ ≡ ⟨ p̲ p̲ ⟩− I/3 ̲

(17)

Note that the structure tensor is symmetric and traceless ( b̲ T = b̲ and Ib ≡ tr( b̲ ) = 0). The two nontrivial invariants of b̲ (i.e., IIb and IIIb) can be expressed in terms of the eigenvalues of b̲ : 3

IIb ≡ tr( b̲ · b̲ ) =

∑ λbi 2

(18)

i=1

and 3

IIIb ≡ tr( b̲ · b̲ · b̲ ) =

∑ λbi 3

(19)

i=1

An equation for ⟨p̲ p⟩̲ follows directly from eq 14 by first multiplying the equation by p̲ p̲ and then integrating over the unit sphere (see the work fo Kim5 and Petty et al.6): ⎛ ∂⟨ p̲ p̲ ⟩ ⎞ ⎜ ⎟ + Pe[W̲ T·⟨ p̲ p̲ ⟩+⟨ p̲ p̲ ⟩W̲ ] ⎝ ∂t ⎠ X̲ Figure 1. Invariants of the anisotropic component of the orientation dyadic.

⎡⎛ 1 ⎞ Doi ⎜ = − FTD ⎢⎝⟨ p̲ p̲ ⟩− I̲ ⎟⎠ − U (⟨ p̲ p̲ ·⟨ p̲ p̲ ⟩−⟨ p̲ p̲ p̲ p̲ ⟩ ⎣ 3 ⎤ : ⟨ p̲ p̲ ⟩)⎥ + λ TPe[ S̲ T·⟨ p̲ p̲ ⟩+⟨ p̲ p̲ ⟩· S̲ − 2⟨ p̲ p̲ p̲ p̲ ⟩: S̲ ] ⎦

For example, in the absence of external forces, all stable and unstable anisotropic equilibrium states are either on the prolate boundary (F) or on the oblate boundary (D) of the invariant diagram. The invariants associated with second-order orientation moments with prolate quadratic forms are all on the Fboundary:

(20)

The two terms on the left-hand-side of eq 20 define the Jaumann derivative of ⟨p̲ p̲⟩, which represents the rate-ofchange of ⟨p̲ p⟩̲ relative to a frame rotating with an angular velocity proportional to the vorticity (see the work of Bird et al.18). The first bracket on the right-hand-side of eq 20 represents rotary dif fusion due to Brownian motion and the mitigation of diffusion due to the excluded volume phenomenon. The second bracket on the right-hand-side accounts for rotary convection due to fluid/particle coupling. Equation 20 is unclosed due to the explicit appearance of the orientation tetradic. A fully symmetric quadratic (FSQ) representation for this statistical property is discussed in section 4 below. Figure 1 identifies all possible realizable microstructures of a suspension associated with a non-negative operator ⟨p̲ p̲⟩ (see the works of Bird et al.;19 Lumley;20 Parks et al.;21 and Petty et al.6). The eigenvalues of the orientation dyadic are non-negative and sum to unity (i.e., tr⟨p̲ p̲⟩ = 1. On the other hand, the trace of the structure tensor is zero (i.e., tr⟨ b̲ ⟩ = 0), so the eigenvalues of the structure tensor are positive and negative. Figure 1 indentifies all physical possible anisotropic orientation states parametrized by the two invariants of the structure tensor (see eq 18 and eq 19 above). The local microstructure of the suspension depends on IIb and IIIb. The quadratic form associated with the orientation dyadic defines the microstructure: ̲ ̲ > 0, Q ( z̲ ; IIb, IIIb) ≡ ⟨ p̲ p̲ : zz

IIb = 6(IIIb/6)2/3 ,

(IIb, IIIb)F ⊆ F‐boundary

(22)

and the invariants associated with oblate quadratic forms are all on the D-boundary: IIb = 6( −IIIb/6)2/3 ,

(IIb, IIIb)D ⊆ D‐boundary (23)

Two-dimensional planar anisotropic states are located on the Bboundary. The states on the B-boundary of Figure 1 are associated with second-order orientation operators with one eigenvalue equal to zero and two unequal positive eigenvalues. The two invariants of the structure tensor b̲ on the B-boundary are related: 2 IIb = + 2IIIb, (IIb, IIIb)B ⊆ B‐boundary (24) 9 For two-dimensional planar isotropic states (point C), ⟨p̲ p⟩̲ has one zero eigenvalue and two equal eigenvalues; therefore, (IIb, IIIb)C = (6/36, −1/36). A fully aligned microstructure forms an ideal nematic phase with one nonzero eigenvalue (point A): (IIb, IIIb)A = (24/36, 8/36). All microstructures associated with axisymmetric suspensions must fall either on the boundaries or within the bounded region of the (IIb, IIIb)-plane. Microstructures outside this domain are unphysical and, do not exist. The realizability domain arises

∀ z̲ ∈ E3 , || z̲ || = 1 (21) C


Article

Industrial & Engineering Chemistry Research [AA ̲ ̲ ] ⇒ AijAkl + Aik A jl + Ail A jk

directly from the algebraic properties of real, symmetric, nonnegative operators. This mathematical property places important theoretical constraints on allowable models for the orientation dyadic. It is noteworthy that as a suspension approaches the nematic state (IIb → 2/3), the rotary diffusion coefficient (defined by eq 12 above) becomes very large.

and [AB ̲ ̲ ] ⇒ AijBkl + Aik Bjl + Ail Bjk + Bij Akl + Bik Ajl + Bik Ajk

4. SYMMETRIC QUADRATIC CLOSURE FOR THE ORIENTATION TETRADIC The orientation tetradic is defined as 2π

⟨ p̲ p̲ p̲ p̲ ⟩ ≡

∫0 ∫0

(25)

01(⟨ p̲ p̲ ⟩) ≡ −

This operator has 6-fold symmetry and 6-fold projection properties, which must be retained by any algebraic closure model that uses the following hypothesis: ⟨ p̲ p̲ p̲ p̲ ⟩ = 0(⟨ p̲ p̲ ⟩)

02(⟨ p̲ p̲ ⟩) ≡

(26)

Clearly, a transpose operation on a tetrad (i.e., p̲ p̲ p̲ p̲) can be executed in six different ways. The operator remains exactly the same after each operation. The transpose operation commutes with the ensemble average operator defined by eq 25, so the tetradic operator ⟨p̲ p̲ p̲ p̲⟩ must retain the 6-fold symmetry structure of the tetrad. Consequently, the right-hand-side of eq 26 must also have 6-fold symmetry. In addition, the trace operation (i.e., projection operation) applied to eq 26 implies that tr[⟨p̲ p̲ p̲ p̲⟩] = tr[0 (⟨p̲ p̲⟩] = ⟨p̲ p̲⟩. Clearly, this operation can also be executed in six different ways with the same result. Thus, this 6-fold projection property must hold for the lefthand-side and the right-hand-side of eq 26. Petty et al.6 used the Cayley−Hamilton theorem (see the work of Frazer et al.22) and developed an irreducible representation for 0 (⟨p̲ p̲⟩) in terms of the following six independent tetradic operators: ⎫ ⎪ ⎪ (c) [ I̲ ][⟨ p̲ p̲ ⟩]; (d) [⟨ p̲ p̲ ⟩][⟨ p̲ p̲ ⟩] ⎬ ⎪ (e) [⟨ p̲ p̲ ⟩·⟨ p̲ p̲ ⟩][ I̲ ]; (f) [ I̲ ][⟨ p̲ p̲ ⟩·⟨ p̲ p̲ ⟩]⎪ ⎭

2 ⟨ p̲ p̲ ⟩: ⟨ p̲ p̲ ⟩[ I̲ I̲ ] + [⟨ p̲ p̲ ⟩⟨ p̲ p̲ ⟩] 35 10 − [⟨ p̲ p̲ ⟩·⟨ p̲ p̲ ⟩ I̲ ] 35

(32)

(33)

⟨ p̲ p̲ p̲ p̲ ⟩ = [1 − C2(IIb, IIIb)]01(⟨ p̲ p̲ ⟩) + C2(IIb, IIIb)02(⟨ p̲ p̲ ⟩)

(34)

The above representation for the orientation tetradic satisfies all 6-fold symmetry operations and all 6-fold projection operations associated with the exact orientation tetradic. Kim6 discovered that solutions to eq 20 are realizable for all combinations of the three physical property groups (i.e., 0 ≤ Pe < ∞; 0 ≤ U < ∞; and, −1 ≤ λT ≤ +1) provided the Cayley− Hamilton coefficient, C2(IIb, IIIb), in eq 34 satisfies the following inequality: C2(IIb, IIIb) ≥

8 + 45IIIb 18(1 + 9IIIb)

(35)

Inequality 35 stems from the condition that all initial orientation states on the boundary of the realizable region (see Figure 1) must remain either on the boundary or be attracted by states within the realizable region. The lower bound for C2(IIb, IIIb) is used in section 5 below to determine the influence of the Péclet number on the shear viscosity of a suspension.

(27)

0 (⟨p̲ p⟩) ̲ can be expressed as a linear combination of two tetradic operators with 6-fold symmetry and 6-fold contraction: (28)

In the above representation, 0 1(⟨p̲ p̲⟩) is first-order in ⟨p̲ p̲⟩ and 0 2(⟨p̲ p̲⟩) is second-order in ⟨p̲ p̲⟩. The scalar coefficients C1 and C2 are explicit functions of the eigenvalues of ⟨p̲ p⟩̲ or, equivalently, the invariants of the structure tensor b̲ (see Figure 1). The closure coefficients C1 and C2 are not independent because eq 28 must also satisfy the 6-fold projection conditions. Thus,

5. SUSPENSION RHEOLOGY The deviatoric component of the total stress for a rigid rod suspension consists of three contributions: a solvent stress due to fluid deformation ( τ̲̂S); a viscous stress due to fluid/particle interactions ( τ̲̂V); and, an elastic stress due to particle/particle interactions ( τ̲̂E). A superposition of these three stresses defines the total deviatoric component of the Cauchy stress:

tr[⟨ p̲ p̲ p̲ p̲ ⟩] = tr[0(⟨ p̲ p̲ ⟩)]

τ̲ ̲̂ = τ̲ Ŝ̲ + τ̲ ̲̂ V + τ̲ Ê̲ ,

= C1tr[01(⟨ p̲ p̲ ⟩)] + C2 tr[02(⟨ p̲ p̲ ⟩)] = (C1 + C2)⟨ p̲ p̲ ⟩

1 5 [ I̲ I̲ ] + [⟨ p̲ p̲ ⟩ I̲ ] 35 35

Equations 28, 29, 32, and 33 define a preclosure for the orientation tetradic ⟨p̲ p̲ p̲ p̲⟩:

(a) [ I̲ ][ I̲ ]; (b) [⟨ p̲ p̲ ⟩][ I̲ ]

⟨ p̲ p̲ p̲ p̲ ⟩ = C101(⟨ p̲ p̲ ⟩) + C2 02(⟨ p̲ p̲ ⟩)

(31)

In the above equations, the dyadic-valued operators A̲ and B̲ are symmetric and may be one of the following independent second-order operators: I ̲ , ⟨p̲ p̲⟩, or ⟨p̲ p̲⟩·⟨p̲ p̲⟩. The tetradicvalued operators 0 1(⟨p̲ p̲⟩) and 0 2(⟨p̲ p̲⟩) are (see the work of Petty et al.6)

π

p̲ p̲ p̲ p̲ Ψ̲(θ , ϕ , t)sin(θ ) dθ dϕ

(30)

T̲ ̂ = −P I̲ + τ̲ ̲̂

(36)

If the solvent is a Newtonian fluid, then a scalar-valued solvent viscosity, which depends on the local temperature and pressure, relates the suspension shear rate to the solvent stress:

(29)

The above result implies that C1 + C2 = 1. Explicit algebraic expressions for 0 1(⟨p̲ p̲⟩) and 0 2(⟨p̲ p̲⟩) can be constructed by using the following two 6-fold symmetric tetradic operators:

τ̲ Ŝ̲ = 2ηŜ S̲ ̂ ηŜ ≡ ( τ̲ Ŝ̲ : e̲ z e̲ y)/(2 S̲ :̂ e̲ z e̲ y) D

(37) DOI: 10.1021/ie503995y Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

Figure 2. Effect of the rotary Péclet number on the time averaged suspension viscosity for U = 27 and λT = 0.987.

[ τ̲ Ŝ̲ ]: [ e̲ z e̲ y] η̂ ≡ ̲̂ [ e ̲ z e ̲ y ] 2[S]:

The microstructure of the suspension couples with the strain rate of the suspension to produce an extra viscous stress due to fluid/particle interactions (see p 521 in the work of Larson4):

=

τ̲ ̲̂ V = 2[cζR ]⟨ p̲ p̲ p̲ p̲ ⟩: S̲ ̂ ηV̂ ≡ ( τ̲ ̲̂ V : e̲ z e̲ y)/(2 S̲ :̂ e̲ z e̲ y)

3ckB*T o C TD DRo

(42)

With ηV ≡ η̂V/[cζR] and ηE ≡ η̂E/[cζR], eq 38 and eq 40 imply that the dimensionless extra viscous viscosity and the dimensionless elastic viscosity depend on the Péclet number, the excluded volume coefficient U, and the microstructure parameters λT, IIb, and IIIb: ηV = ⟨py pz pz py ⟩

(43)

and (39)

ηE =

Doi1 developed a model for the elastic stress based on the Onsager free energy with the result that

1 [⟨p p ⟩−U (⟨py p̲ ⟩·⟨ pp ̲ z ⟩−⟨py pz p̲ p̲ ⟩: ⟨ p̲ p̲ ⟩)] Pe y z

(44)

The components of the orientation tetradic operator are defined by eq 34 with a closure coefficient defined as follow (cf. inequality 35 above):

⎡⎛ ⎤ 1 ⎞ τ̲ Ê̲ = 3ckB*T ⎢⎜⟨ p̲ p̲ ⟩− I̲ ⎟ − U (⟨ p̲ p̲ ⟩·⟨ p̲ p̲ ⟩−⟨ p̲ p̲ p̲ p̲ ⟩: ⟨ p̲ p̲ ⟩)⎥ ⎣⎝ ⎦ 3 ⎠

C2 =

ηÊ ≡ ( τ̲ Ê̲ : e̲ z e̲ y)/(2 S̲ :̂ e̲ z e̲ y)

8 + 45IIIb 18(1 + 9IIIb)

(45)

For IIIb = 0 (isotropic), C2 = 4/9; and for IIIb = 6/18 (nematic), C2 = 3/9. With the above closure, eq 20 governs the temporal behavior of the orientation dyadic ⟨p̲ p⟩. ̲ The initial condition must have invariants consistent with the principle of realizability (see Figure 1). The results summarized by Figure 2 were developed from a planar anisotropic state (i.e., point C in Figure 1). For simple shear flows, eq 45 guarantees that solutions to eq 20 will be realizable for 0 ≤ Pe < ∞, 0 ≤ U < ∞, and, −1 ≤ λT ≤ +1. In summary, the underlying microstructure of the suspension is governed by eq 20, which is a dyadic-valued nonlinear ordinary differential equation for the second-order moment of the orientation distribution function. The system of nonlinear ordinary differential equations was solved by using a fourthorder Runga−Kutta algorithm (see the work of Kim,5 Appendix

(40)

The scalar-valued coefficient [3ck*B T] has units of stress and does not depend on the invariants of the microstructure. The parameter “c” represents the volume fraction of the dispersed phase. The first term on the right-hand-side of eq 40 represents the stress induced by rotary Brownian motion; the second term is the stress caused by the excluded volume potential. For homogeneous shear flows 2 S̲ :̂ e̲ z e̲ y ≡ γ ̇ ⇒ ∇ u̲ = [Pe] e̲ y e̲ z

γ̇

= ηŜ + ηV̂ + ηÊ

(38)

For this situation, the relationship between the local kinematics of the suspension is connected to the local extra stress by a fourth-order orientation moment. The scalar-valued coefficient [ cζR] has units of viscosity and does not depend on the invariants of the microstructure: [cζR ] ≡

[ τ̲ Ŝ̲ ]: [ e̲ z e̲ y] + [ τ̲ ̲̂ V ]: [ e̲ z e̲ y] + [ τ̲ Ê̲ ]: [ e̲ z e̲ y]

(41)

Equations 37, 38, and 40 imply that the suspension viscosity is E


Article

Industrial & Engineering Chemistry Research G). The components of ⟨p̲ p⟩̲ and ⟨p̲ p̲ p̲ p⟩̲ were averaged over several periods to obtain time-averaged results.

of the microstructure from an isotropic state to a nematic state in the absence of an external field, they also predict unrealizable behavior for simple shear flows (see the works of Chaubal and Leal29 and Feng et al.30). Interestingly, Larson31 used the Doidecoupling approximation for the excluded volume potential and the HL-closure for the convective flux in eq 20. This approach also predicts unrealizable behavior. Cintra and Tucker32 developed a closure for the orientation tetradic based on an orthotropic operator. This approach assumes that the symmetry properties of a tetradic operator can be represented in terms of the eigenvectors of the orientation dyadic operator. The closure coefficients of the tetradic operator are obtained by fitting the moment equation with the exact solution based on a spherical harmonic expansion or the orientation density function (CT-data set). Although this approach can clearly reproduced the limited bench mark results used to determine the closure coefficients, it is unlikely that predictions for other flows will consistently produce realizable results. It is noteworthy that the CT data set could be used to benchmark the lower bound on the Cayley−Hamilton closure coefficient defined by Ineq 35 above. Equation 20 has been used by numerous researchers as a practical template for more than 60 years to understand the behavior of simple shear flows of suspensions. Although eq 20 is a valuable resource for understanding flows of complex suspensions, it is essential that the closure for the orientation tetradic be consistent with the exact symmetry and projection properties of a tetradic. This is a strong mathematical constraint on the mapping hypothesis expressed by eq 26. Furthermore, it is also essential that the closure for the tetradic operator support the requirement that the orientation dyadic determined by eq 20 is non-negative. It is noteworthy that the closure formulated herein retains the 6-fold symmetry and the 6-fold projection conditions characteristic of tetradic operators and that the closure produces orientation dyadics that are nonnegative for all simple shear flows.

6. DISCUSSION Shear Thickening Prediction. Figure 2 shows how the Péclet number affects the time-averaged extra viscosity coefficient (i.e., the sum of eqs 43 and 44). For these calculations, U = 27 and λT = 0.987. The tumbling parameter correspond to a rigid rod with L/d ≃ 12. The relatively large value of U was selected to produce a microstructure with a periodic director (i.e., the eigenvector associated with the largest eigenvalue of the orientation dyadic). For Pe = 0, the microstructure changes from an isotropic state to an anisotropic state for U = 5, which is consistent with experimental observations and other theoretical models. See the work of Kim2 (Chapter 2) for a comprehensive treatment of the microstructure in the absence of an external field. The Doi tube dilation model for semidilute suspensions (see eq 12 above) supports the existence of a director tumbling regime for rotary Péclet numbers around 100 and a director alignment regime for rotary Péclet numbers above 1000 (see Figure 2). A director wagging regime occurs for Péclet numbers around 500. At high Péclet numbers, the extra viscosity coefficient approaches a Newtonian plateau lower than the low Péclet number plateau. The model predicts shear thickening behavior between the director tumbling regime and the director wagging regime (50 < Pe < 500). Shear thinning occurs in the director wagging regime (i.e., 500 < Pe < 1000). Tube dilation and director tumbling are key physical processes that support the rheological response illustrated by Figure 2. Method-of-Moments. There are several ways to study the RS-equation defined by eq 2. One approach is to develop a solution using spherical harmonics (see the works of Larson23 and Chaubal and Leal24). Another approach is based on the method-of-moments. Developing an explicit expansion for the orientation density function is a complicated process and requires significant computational resources. On the other hand, developing a solution based on low order moments requires a closure approximation that satisfies the realizability condition. Doi1 used a quadratic closure approximation (also called a decoupling approximation) for the orientation tetradic: ⟨ p̲ p̲ p̲ p̲ ⟩ = [⟨ p̲ p̲ ⟩][⟨ p̲ p̲ ⟩]

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was funded by the National Science Foundation Industry/University Cooperative Research Center for Multiphase Transport Phenomena (NSF/ECC 0331977; NSF/I IP 093474). The authors gratefully acknowledge the members of the Center for their support. Y.K. received an NSF summer fellowship in support of his research on complex fluids with Professor Takimoto at Nagoya University, Japan.

(46)

The above closure, albeit appealing, does not retain the 6-fold symmetry and the 6-fold projection properties of a fourth-order moment. Note that the right-hand side of eq 46 is only one of six independent operators needed to represent a tetradic operator based on eq 26. Unfortunately, eq 46 lacks completeness for a theoretical analysis of eq 20. Most significantly, however, it is well-known that the use of the decoupling approximation with eq 20 yields unrealizable solutions. Hand25 used a first-order closure approximation of the orientation tetradic. This closure is equivalent to setting C2(IIb, IIIb) = 0 in eq 26 above. Hand’s closure is limited to microstructures near the isotropic state (see point E in Figure 1) and also produces unrealizable behavior for strongly anisotropic suspensions. Hinch and Leal26,27 as well as Tucker28 developed hybrid approaches by using a superposition of the Hand model and a model near the nematic state (point A in Figure 1). Although these approaches could predict a transition

■

F

NOMENCLATURE ao = radius of dilation tube [=] m b̲ = anisotropic orientation dyadic c = volume fraction of rigid rod phase [cζR] = fluid/particle viscosity coefficient [=] Pa s C1(IIb, IIIb) = closure coefficient C2(IIb, IIIb) = closure coefficient CoTD = tube dilation coefficient d = diameter of rigid rods [=] m ⟨DR⟩ = rotary diffusion coefficient [=] 1/s DOI: 10.1021/ie503995y Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research DoR = dilute rotary diffusion coefficient [=] 1/s (ex, ey, ez) = Cartesian base vectors FDoi TD = normalized tube dilation factor (eq 12) I ̲ = unit operator for a 3D Euclidean vector space k*B = dispersed phase Boltzmann coefficient [=] J/(K m3) L = length of hypothetical dilation tube [=] m p̲ = orientation vector of unit length p̲ ̇ = angular velocity of orientation vector [=] 1/s p̲ ċ = angular velocity of convective orientation vector [=] 1/s p̲ p̲ = instantaneous orientation dyad p̲ p̲ p̲ p̲ = instantaneous orientation tetrad P = Cauchy pressure [=] Pa Pe = Péclet number Q = quadratic form of the orientation dyadic ⟨p̲ p⟩̲ t = time made dimension with CoTDDoR T = absolute temperature [=] K T̲ ̂ = Cauchy stress [=] N/m2 u̲ ̂ = suspension velocity [=] m/s U = excluded volume coefficient (see eq 13) ΔÛ = suspension excluded volume potential [=] J/m3 VE = volume of dilation tube [=] m3 VR = volume of rigid rod [=] m W̲ ̂ = vorticity dyadic made dimensionless with γ̇ X̲ ̂ = reference position of material element [=] m z = arbitrary constant unit vector in 3D vector space

(3) Doi, M.; Edwards, S.F. Theory of Polymer Dynamics; Oxford University Press: London, 1986. (4) Larson, R. G. The Structure of Complex Fluids; Oxford University Press: New York, 1999. (5) Kim, Y.-C. Prediction of Rheological Properties of Structured Fluids in Homogeneous Shear Based on a Realizable Model for the Orientation Dyad. PhD Dissertation, Michigan State University, 2006. (6) Petty, C. A.; Parks, S. M.; Shao, S. M. Flow Induced Alignment of Fibers. In Proceedings of 12th International Conference on Composite Materials, ICCM12, Paris, July 5−9,1999; paper 339. (7) Jeffery, G. B. The Motion of Ellipsoidal Particles Immersed in Viscous Fluid. Proc. R. Soc. London A 1922, 102, 161−179. (8) Batchelor, G. K. Brownian Diffusion of Particles with Hydrodynamic Interaction. J. Fluid Mech. 1976, 74, 1−29. (9) Batchelor, G. K. Sedimentation in a Dilute Polydisperse System of Interacting Spheres, Part I: General Theory. J. Fluid Mech. 1982, 119, 379−408. (10) Bibbo, M. A.; Dinh, S. M.; Armstrong, R. C. Shear Flow Properties of Semiconcentrated Fiber Suspensions. J. Rheol. 1985, 29 (6), 905−929. (11) Trevelyn, B. J.; Mason, S. G. Particle motions in sheared suspensions. I. Rotations. J. Colloid Science 1951, 6, 354−367. (12) Bretherton, F. P. The Motion of Rigid Particles in a Shear Flow at Low Reynolds Number. J. Fluid Mech. 1962, 14, 284−304. (13) Anczurowski, E.; Mason, S. G. The Kinetics of Flowing Dispersions III. Equilibrium Orientations of Rods and Discs (Experimental). J. Colloid Interface Sci. 1967, 23, 533−546. (14) Frattini, P. L.; Fuller, G. G. Rheo-Optical Studies of the Effect of Weak Brownian Rotations in Sheared Suspensions. J. Fluid Mech. 1986, 168, 119−150. (15) Kuzuu, N.; Doi, M. Constitutive Equation for Nematic Liquid Crystals under Weak Velocity Gradient Derived from a Molecular Kinetic Equation. J. Phys. Soc. Jpn. 1983, 52 (10), 3486−3494. (16) Ilg, P.; Karlin, I. V.; Ö ttinger, H. C. Generating moment equation in the Doi model of liquid-crystalline polymers. Phys. Rev. E 1999, 60 (5), 5783−5787. (17) Onsager, L. The Effect of Shape on the Interaction of Colloidal Particles. Ann. N.Y. Acad. Sci. 1949, 51, 627−659. (18) Bird, R. B.; Armstrong, R. C.; Hassager, O. Dynamics of Polymeric Liquids: Vol. 1, Fluid Mechanics, second ed.; WileyInterscience: New York, 1987. (19) Bird, R. B.; Curtiss, C. F.; Armstrong, R. C.; Hassager, O. Dynamics of Polymeric Liquids: Vol. 2, Kinetic Theory, second ed.; WileyInterscience: New York, 1987. (20) Lumley, J. L. Computational Modeling of Turbulent Flows. Adv. Appl. Mech. 1978, 18, 123. (21) Parks, S. M.; Weispfennig, K.; Petty, C. A. An Algebraic Preclosure Theory for the Reynolds Stress. Phys. Fluids 1998, 10 (13), 645−653. (22) Frazer, R. A.; Duncan, W. J.; Collar, A. R. Elementary Matrices and Some Applications to Dynamics and Differential Equations; Cambridge, 1960. (23) Larson, R. G. Constitutive Equations for Polymer Melts and Solutions; Butterworth: Boston, 1988. (24) Chaubal, C. V.; Leal, L. G. Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems”, Part I. Method. J. Non-Newtonian Fluid Mech. 1997, 70, 125−154. (25) Hand, G. L. A Theory of Anisotropic Fluids. J. Fluid Mech. 1962, 13, 33−46. (26) Hinch, E. J.; Leal, L. G. Constitutive Equations in Suspension Mechanics. Part 1. General Formulation. J. Fluid Mech. 1975, 71 (Part 3), 481−495. (27) Hinch, E. J.; Leal, L. G. Constitutive Equation in Suspension Mechanics. Part 2. Approximate forms for a suspension of rigid particles affected by Brownian Rotations. J. Fluid Mech. 1976, 76, 187− 208. (28) Tucker, C. L. Predicting Fiber Orientations in Short Fiber Composites. In Proceedings of Manufacturing Int’l; ASME Publication: New York, 1988; p 95.

Greek Symbols

γ̇ = strain rate [=] 1/s η = suspension viscosity made dimensionless with [cζR] ηE = elastic viscosity made dimensionless with [cζR] ηS = solvent viscosity made dimensionless with [cζR] ηV = extra viscosity made dimensionless with [cζR] θ = Euler angle associated with orientation vector λbi = eigenvalue of anisotropic orientation dyadic (i =1, 2, 3) λT = tumbling coefficient defined by eq 6 τ̲̂ = deviatoric stress of suspension [=] Pa τ̲̂E = elastic stress [=] Pa τ̲̂S = solvent stress [=] Pa τ̲̂V = fluid/particle stress [=] Pa ϕ = Euler angle for the orientation vector χ̲̂ (X̲ ̂ , t )̂ = motion of a material element [=] m Ψ = rotary orientation density function (see eqs 2 and 3) Operators

⟨p̲ p̲⟩ = ensemble average (see eq 15) ⟨p̲ p̲ p̲ p⟩̲ = ensemble average (see eq 25) Ib = first invariant of b̲ ⇒ Ib = tr( b̲ ) = ∑3i=0 λbi = 0 IIb = second invariant of b̲ ⇒ IIb = tr( b̲ · b̲ ) = ∑3i=1 λbi2 IIIb = third invariant of b̲ ⇒ IIIb = tr( b̲ · b̲ · b̲ ) = ∑3i=1 λbi3 [A̲ B̲ ] = fully symmetric 6-fold tetradic operator (see eq 31) 0 (⟨p̲ p⟩) ̲ = preclosure ∇̂ = spatial gradient operator

■

REFERENCES

(1) Doi, M. Molecular Dynamics and Rheological Properties of Concentrated Solutions of Rodlike Polymers in Isotropic and Liquid Crystalline Phases. J. Polym. Sci. 1981, 19, 229−243. (2) Doi, M.; Edwards, S. F. Dynamics of Rod-like Macromolecules in Concentrated Solution. J. Chem. Soc. Faraday Trans. II 1978, 74 (Part 1), 560−570. Part 2, 918−932. G


Article

Industrial & Engineering Chemistry Research (29) Chaubal, C. V.; Leal, L. G. Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems. Part II. The Effect of Flow Perturbations on the simple shear behavior of LCPs. J. Non-Newtonian Fluid Mech. 1997, 82, 25−55. (30) Feng, J.; Chaubal, C. V.; Leal, L. G. Closure Approximations for the Doi Theory: Which to Use in Simulating Complex Flows of Liquid-Crystalline Polymers? J. Rheol. 1998, 42 (5), 1095−1119. (31) Larson, R. G. Arrested Tumbling in Shearing Flows of Liquid Crystal Polymers. Macromolecules 1990, 23, 3983−3992. (32) Cintra, J. S.; Tucker, C. L., III Orthotropic Closure Approximations for Flow-Induced Fiber Orientation. J. Rheol. 1995, 39 (6), 1095−1122.

H


Microstructure and Rheology of Rigid Rod Suspensions - Industrial

Recommend Documents