A deformation tensor model of Brownian suspensions of orientable particles—the nonlinear dynamics of closure models

Non-Ncwtonin Fluid Mechanics ELSEVIER

J. Non-Newtonian Fluid Mech., 64 (1996) 43-69

A deformation tensor model of Brownian suspensions of orientable particles the nonlinear dynamics of closure models A n d r e w J. Szeri*, David J. Lin University of California, Irvine, Department of Mechanical and Aerospace Engineering, Irvine, CA 92717-3975, USA Received 15 May 1995; revised 18 December 1995

Abstract A new model is developed for the evolution of the orientation distribution within a flowing Brownian suspension of orientable particles, e.g. fibers, disks, rods, etc. Rather than solving the full Fokker Planck equation for the orientation distribution function for the suspended phase, and in place of the usual approach of developing a moment closure model, a new approach is taken in which an evolution equation is developed for an approximate, simplified deformation of the orientable particles associated with a material point. The evolution equation for the remaining degrees of freedom in the assumed class of deformations is developed from the Fokker Planck equation; it is as quick to integrate as direct moment tensor evolution equations. Because the deformation is restricted to a special class, one can show a priori that the model always gives physically sensible results, even in complicated flows of practical interest. The nonlinear dynamics of the model in unsteady, three-dimensional flows is considered. It is shown that the model predicts bounded deformations in all flows. The model has a unique global attractor in any steady, three -dimensional flow, provided the rotary Brownian diffusivity is non-zero. In contrast, commonly used moment closure approximations can only be shown to have a unique global attractor when the rotary Brownian diffusivity is large enough (i.e. when the flow is weak enough). This may be the explanation of why a supurious (multiple) attractor was recently observed at high Pdclet number in un~)rm shear.flow of a dilute suspension modeled using the first composite closure of Hinch and Leal, by Chaubal, Leal and Fredrickson (J. Rheol., 39(1) (1995) 73-103).

Keywords: Brownian motion; Fibers; Modelling; Nonlinear dynamics; Suspensions

I. Introduction

Suspensions of orientable particles have physical properties that depend sensitively on the

conjormation of the suspended phase. By this one generally means the departure from an * Corresponding author. 0377-0257/96/$15.00 © 1996 - Elsevier Science B.V. All rights reserved SSDI 0377-0257(95)01427- 5

44

A.J. Szeri, D.J. Lin / J . Non-Newtonian Fluid Mech. 64 (1996) 43-69

equilibrium state of the distribution with respect to concentration, orientation or deformation of the suspended particles. Such a departure from equilibrium microstructure is brought about by the flow of the suspension; hence there is a subtle interplay between the flow and the evolving microstructure, each of which affects the other. I.I. Kinetic theory and moment closure models

An essential ingredient in a mathematical model of suspensions that will faithfully predict both the flow and the flow-induced microstructure is some means to account for the evolution of microstructure in what may be a very complicated flow. Intense effort has been directed toward this problem, as is apparent from Bird et al. [1], or Larson [2]. Past effort has been concentrated into: (i) kinetic theory, aimed at developing an evolution equation for the distribution function of the suspended phase, and (ii) moment closure models, aimed at developing a simplified model for practical work consisting of an evolution equation for the moments of the distribution function. The moments are required to account for the particle contribution to the state of stress at a point in the suspension. The full kinetic theory is too unwieldy for use in the solution of practical problems. Moment closure models give rise to pathological behavior even in simplified theoretical flow fields; this behavior is reviewed by Advani and Tucker [3] and by Szeri and Leal [4,5]. In this paper, we consider dilute Brownian suspensions of orientable particles in a Newtonian solvent. The particles are assumed to be small with respect to length scales over which the flow varies, with the consequence that only the velocity gradient tensor of the flow evaluated at the center of mass of a given particle will influence the evolution of its orientation. The centers of mass of the particles are assumed to follow a known path through the flow. The restriction to dilute suspensions can be overcome. In Szeri [6] we develop a deformation tensor model for liquid crystalline polymers, which are far from dilute. The result of a phase space kinetic approach to the modeling of the dynamics of a fiber suspension is a Fokker-Planck (or Smoluchowski) equation of the form ~--t+~uu'{O(KU--(u'KU)U)}--~UU"

DR~u =0"

(1)

Here, ~,(u, t) is the orientation distribution function of a given material point of the fiber suspension, K is the "equivalent velocity gradient tensor" of the flow, evaluated at the location of the material point in question (K = f~ + GE, where fl is the vorticity tensor, G is the shape factor of particles, and E is the rate-of-strain tensor); u is the (normalized) direction-vector of a rod-like particle, and DR is the rotary Brownian diffusivity. For ellipsoidal particles, G is derivable from the aspect ratio: one finds G = 0 for spherical particles and G = 1 for infinite aspect ratio fibers. It is important to note that K may depend on time either through time-dependence of the Eulerian velocity field or through the motion of the particle through different regions of the flow (or both). The operator ~/~u is the gradient in conformation space, projected onto the sphere of orientations, i.e. ~/~u - ( I - uu)Vu. For a detailed description of the model we refer the reader to Bird et al. [1], or to Larson [2]. Eq. (1) was developed originally as a mathematical description of a dilute suspension; however Folgar and Tucker [7] have developed a mathematical formulation for concentrated suspensions of non-Brownian fibers (G = 1) in which the role of the "Brownian term" is as a model of interaction between particles.

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

45

The particle contribution to the state of stress at a point in the suspension is phrased in terms of weighted orientation space averages (moment tensors) of the distribution function ¢,, associated with a material point. For example the second m o m e n t tensor is

( u u ) = f uug'(u, t) d2u,

(2)

where the integral is taken over orientation space. The second and fourth moments are required by the theory. The closure models developed thus far for Brownian suspensions are all of the direct m o m e n t tensor evolution equation type. If we multiply (1) by uu and integrate over orientation space, we obtain the well-known second m o m e n t tensor evolution equation: d dt

(uu) = ~(uu)

- 2(uuuu)

+

(uu)f~ -r + G ( E ( u u ) + ( u u ) E "r

(

1)

• E) + 6DR ( u u ) -- ~ I .

(3)

A closure model is required as the equation for ( u u ) depends on (uuuu). We discuss specific closure models in Section 3.

1.2. The deformation tensor model In this paper, we propose a new approach to the development of a practical model for the evolution of the microstructure in a flowing suspension. We shall develop a model for the evolution of the deformation of the suspended phase, rather than approximating the evolution of a m o m e n t tensor directly. Just as in the development of m o m e n t closure models, our point of departure is kinetic theory. However, the innovation we propose is based on the idea of restricting the deformation of the suspended phase to an assumed class of deformations appropriate to the material. Using the new model we propose, one can obtain whatever moments are required of the distribution function associated with the approximate deformation. Therefore, there is no m o m e n t closure model; instead, it may be useful to think of our new approach as a closure model for the deformation rather than a conventional m o m e n t closure. In order to develop a deformation tensor model, we begin by splitting the F o k k e r - P l a n c k equation into two equations: ~-

~u -(g, ti) = 0

(4)

and i t = K U - - ( t t "KU)U-- D R

log ~b ~u

-

(5)

Next, we replace (4) by the equivalent statement in the Lagrangian frame in conformation space [4,51:

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

46

~,*(t, uo)

det

O*(0, Uo)

(6)

3

The Lagrangian form of the distribution function ~*(t, uo) is simply the orientational probability attached to an orientation u(t, Uo) that evolves according to the evolution equation (5). The quantity ~u(t, Uo)/OUo is the deformation gradient tensor between reference and deformed configurations of the fiber suspension. Now we state the approximate form of the deformation that we assume in the remainder of the paper: Q(t)Uo

.(t, .0) - [ie(t).oll,

(7)

where Q is a second rank tensor that depends on time. This form of the deformation is motivated by the following considerations. For a dilute non-Brownian suspension of rod-like particles in a flow, the deformation (7) is exact, with

dQ dt

-

K(t)Q

(8)

and initial condition Q ( 0 ) = I ; see Bretherton [8]. If we take G = 1 then Q is simply the deformation gradient tensor for the underlying fluid flow. When G ~: 1, we call Q the equivalent deformation gradient tensor, or simply the deformation tensor. The deformation (7) has the attractive feature that all orientation vectors are forever of unit length. Moreover, as a consequence of fluid incompressibility tr(Vv)(t)= 0, and therefore det[Q(t)] = 1 for all time [9]. In fact, as we shall see, it is possible to replace t¢ by any other deviatoric (i.e. traceless) second rank tensor, and still preserve det[Q(t)] = 1. In the present work, the goal is to retain the elegant and simple form of the deformation (7), but to modify the evolution equation for Q (8) in a way that accounts for the presence of Brownian motions. In this paper, we develop the model

where h is a scalar-valued function at our disposal to fine-tune the model. Note that the additional Brownian term of this equation makes use of the deviatoric part of Q rQ I as a finite strain tensor, in much the same way as (3) makes use of (uu). In both models, the Brownian motions lead to a term proportional to a finite strain measure of the suspended phase. When Brownian motions are considered, the advantage of the formulation (9) over kinetic theory (1) is clear: one solves an ordinary differential equation rather than a partial differential equation at every material point of the suspension. In addition, (9) guarantees a physically sensible solution (as is shown below), something that solution of direct moment evolution equations does not.

1.3. Related work Our work on the development of a deformation tensor model is related to the work of Currie [10] on a simplified constitutive equation for concentrated polymer solutions and melts; which

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

47

was more recently applied to non-Brownian fiber suspensions by Malamataris and Papanastasiou [11]. These authors based their work on the observation that in the absence of Brownian motion, the dynamics of suspended, orientable particles of infinite aspect ratio may be determined from the deformation gradient tensor of the flowing fluid. Very recently, Rao, et al. [12] have extended this notion to non-Brownian suspensions of orientable particles of finite aspect ratio, by noting that the equivalent deformation gradient tensor can play the role of the deformation gradient tensor. We refer the interested reader to Szeri and Leal [9] for a detailed description of the equivalent deformation gradient tensor and a closely related study of particle orientation dynamics. In a sense, the present paper can be viewed as a further development along the line of research initiated by Currie. Our goal is to develop a model to include the deterministic effects of Brownian motion on the evolution of the equivalent deformation gradient tensor, just as Rao et al. realized how to incorporate the effects of finite aspect ratio particles. This will have ramifications for our understanding of the nonlinear dynamics of particles suspended in unsteady, three-dimensional flows, which is based on the equivalent deformation gradient tensor; see Szeri and Leal [9] and Szeri [13].

2. Development of the model In developing the deformation tensor model, there are several important features that should be kept in mind. First, and foremost, one wants to be sure to develop a model which leads to physically reasonable deformations. In other words, there should be a one-to-one correspondence between every particle in the reference configuration and every particle in the deformed configuration. Stated differently, the mapping from reference to deformed configuration should be invertible. Because it is described in terms of a tensor Q, this clearly implies the requirement 0 < det Q < ~ . It is a simple matter to show that this is equivalent to det Q = 1, because the magnitude of Q cancels out of the deformation (7). Hence it is without loss of generality that we permit only det Q(t) = 1 for all time. This is equivalent to restricting our modifications of (8) to deviatoric (i.e. traceless) tensor coefficients of Q on the right-hand side. Finally, we would like to model the action of the Brownian term in (1) in a modified version of (8). What is required is a term that drives the system toward isotropy in the absence of flow. Clearly, the Brownian term added to (8) should depend on the instantaneous state of the suspension, and not on the history of deformations. Thus, we expect that the Brownian term will involve Q but not derivatives of Q with respect to time, or integrals of Q with respect to time. This observation is of critical importance in restricting the possibilities for the modification of (8) that we propose.

2.1. Brownian motions In this section we develop an expression for the Brownian term that we shall add to (8). This involves writing the Brownian term in (5) in terms of the Lagrangian form of the distribution function given by (6).

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

48

For the deformation we have assumed, and assuming also an isotropic reference configuration (and det Q = 1), we show in the Appendix that (6) leads to the result: ~,*(t, u0) = ~(u, t) -

Ila(t)u°l[3 4re

1

4zcl[a-,(t)ull 3"

(10)

This allows us to calculate the Brownian term of (5) directly: log~, _ ( I - uu)3llO(t)Uol[2Q - T Q - 'u . 8u

(11)

This is the composition of a linear term Q - T Q - i u with the projection operator ( I - uu), which projects a vector onto the tangent plane of the sphere of orientations anchored at the orientation u. The appearance of this projection is natural, as the orientable particles are of fixed length--hence there can be no component of the Brownian term in the direction u. An addtional source of nonlinear dependence on orientation is the term [Ie(t)Uoll 2. Before incorporating the Brownian term into (8), we review how the non-Brownian equation (8) is derived from (5) with DR = 0. To begin, we substitute (7) into (5) with DR = 0; this yields d

Q(t)Uo

- ( I - uu)t¢ Q(t)Uo

at Ila(t)Uo[I

IIQ(t)uol[ "

Next, we differentiate the left-hand side and collect terms to obtain

1

(dQ

ileft).o[i (t-.,,)

)

-dT- xQ u0 = 0.

(12)

Because this equation must apply for all initial orientations Uo, we obtain (8). In order to incorporate the new Brownian term, we simply amend this simple derivation with (11). Following precisely the same steps, we find that (12) is modified to read

iiQ(t)uol[ (t-,,,,)

I--d;- KQ-

3DR[[Q(t)Uo]I2Q

- r

1Uo--0.

(13)

However, we cannot proceed as before. It is not possible to conclude that the term inside the square brackets must be zero in (13) because this term depends on u0. This is where we shall be forced to make an approximation in order to continue. We shall replace IIQ(t),,oll 2 by an isotropic scalar function of the deformation h ( Q - T Q ) (which, of course, does not depend on u0). Under this type of approximation, the model Brownian term now reads log ~, Ou

(14)

~ (I- uu)h(Q~Q)Q-'rQ -'u.

The reader will note that this approximation is now linear, except for the projection. Hence, (13) simplifies to

'

E

ilQ(t)Uo[I ( I - ,,,,) ~

- xQ - 3DRh(QrQ)Q -r

1

Uo = 0.

Because this equation must be satisfied for all Uo, this leads to

(15)

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

49

dQ = tcQ + 3DRh(QTQ)Q T. dt

To obtain (9) we subtract the trace of the tensor coefficient of Q on remains is to find a suitable approximation for h. It is worth approximation we have made is to simplify the Brownian term to ramifications of this approximation will emerge in the tests of the

the right-hand side. All that emphasizing that the only be linear in u. The physical model.

2.2. An approximation for the deformation-dependent jactor--fine-tuning the model Now we consider I[Q(t)uoll 2, which we intend to approximate by h(QTQ). The simplest course of action is simply to average the term IlQ(t)Uo[I2 over the sphere of orientations. This yields 1

h(QTQ) = (llQ(t)Uol[ 2> = -~ tr(QTQ).

(16)

As one might suspect, this approximation works well when the orientation distribution function is not far removed from isotropy. However, it is not the best approximation for practical work. It is possible to improve the model, although it is important to note that the model is not very sensitive to the choice of h. We are after a slight improvement of the model in the remainder of the section. We shall make use of an available analytic solution, given in Bird et al. [1], to the Fokker-Planck equation (1) to compute the second moment tensor for the steady state orientation distribution in steady, irrotational flow fields. In an irrotational flow field, the asymptotic orientation distribution is given by ~ ' ( u ) = j exp

u'Ku

,

(17)

where J is a normalization factor. Next, for the same irrotational flow field, we shall insist that the asymptotic (steady) solution of (9) should match the analytic solution for the distribution function. By "match", we mean that eigenvalues of the second moment tensors obtained by the two approaches should agree, insofar as is possible. This will suffice to determine a suitable form for h(QTQ). Note that we shall make use of the results developed in the Appenidx for the analytic calculation of (uu), given Q. This strategy will now be pursued for two different irrotational flow fields, namely uniaxial/biaxial extensional flow, and general irrotational flow. For purposes of exposition only do we pursue this approach for uniaxial/biaxial extensional flow, as the result is a special form of h that is restricted in application to problems where QTQ and ( u u ) have two equal eigenvalues. This is not very general. 2.2.1. Uniaxial/biaxial extensional flow In the case of uniaxial/biaxial extensional flow, the orientation distribution is axisymmetric when viewed in the coordinate system in which

50

0j

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 4 3 - 6 9

~=G

~ 0

0 , - 2d

where ~ is a parameter which describes the type and strength of the flow field; ~ > 0 denotes a biaxial extensional flow and ~ < 0 a uniaxial extensional flow. The asymptotic orientation distribution is given by specialization of (17):

1 [G~ (1 - 3 sin 2 0) 1 ;

~,(u) = ~ exp

this is written in modified spherical polar coordinates in which 0 - - 0 corresponds to the 1 - 2 plane, 0 = _+ re/2 corresponds to the positive and negative parts of the 3-axis, and u - (cos 0 cos o-, cos 0 sin o-, sin 0). The normalization factor J may be obtained by integration; the fully determined form of ~, is 1 (3~*~ 1/2 e x p ( - 3/2i* sin 2 0) 0 = ~-~ \ - ~ - / t er~(3~)~-- ~ ,

(18)

where i* = G~/DR. In the present case of uniaxial/biaxial extensional flow, it is most useful to compute
( 2 ~'/2 e x p ( - 3e*/2) k,~J erf[(3~*/2)'/2] "

(19)

It is crucial to note that this expression for (uu)33 is m o n o t o n i c in the flow parameter ~*. This result, which is exact, will now be c o m p a r e d with the result one obtains from the deformation tensor model (9). In the eigenvector basis for K we have

QTQ =

[10 0] 22

0

0

1/2122

.

(20)

where we have made use of det Q = 1. Note 21 = ),2 -= 2 when there is uniaxial symmetry. In the Appendix, we obtain an exact expression for the second m o m e n t tensor associated with Q. The 3-3 c o m p o n e n t is 23 )3/2 1og[23/2 + ( _ 1 + )`3)1/2] ~U//)33 = 1 + ~1 - 2 + ( - 1 + 23) 3/2 (21) This expression for (uu)33 is monotonically decreasing in the parameter Now, we shall attempt to choose h so that the steady state solution tensor model will have the same m o m e n t tensor c o m p o n e n t (uu)33 as exact solution. At an asymptotic steady state, Q = 0; hence from (9) we to=

-3DRh(QTQ)[Q-YQ

1 -- -~ 1 tr(Q -TQ

l)/]

2. from the deformation one obtains from the have

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

51

The 1-1 component yeilds the equation 2~* h - 23~;

(22)

the other two diagonal components yield the same equation. Therefore it is possible to obtain a single scalar value for h in a given flow by computing it from appropriate values of ~* and 2. Finally, to determine h, we only need to determine the implicit function ~*(2) defined by forcing equality of the exact (uu)33 associated with ~* and the approximate ~uu)33 associated with 2. In the present case we can only detemine the implicit function numerically, in view of the complicated expressions for (uu)33. A useful approximate form of the function ~*(2) is ~*(2) ~ - 1.44974 x 1 0 - 6 2 - 4 _~_0.000801582 3 - 0.07768822 2 _ 0.6525572- i _ 0.995121 + 1.528572 + 0.0070877422 + 0.11471723 - 0.0024000224.

(23)

Over the range - 50 < ~* < 100, this approximate form is very accurate. Now, we simply use this relationship ~*(2) in (22) to obtain a functional form for h(2) that is guaranteed to yield the correct steady state behavior in uniaxial/biaxial extensional flow. This function is plotted in Fig. 1, vs. the flow parameter i* for ease of interpretation. One can see that h is not overly sensitive to the type of flow. Hence the improvements in the model by choosing the "best" h are modest, but worthwhile. A rather subtle benefit of the deformation tensor model with h finely tuned is as follows. Suppose one is computing either a uniaxial or biaxial extensional flow--possibly with ~* time-dependent. Because: (i) the Brownian term of (1) depends only on the instantaneous deformation, (ii) the finely tuned deformation tensor model yields a O_TQ which has uniaxial symmetry at all times in such a flow, and (iii) the model gives an exact asymptotic match for (uu) compared with the full F o k k e r - P l a n c k equation, we can say the following. The transient 2

h

0

50

i00

Fig. 1. A plot of h vs. ~* for an exact match of from the deformation tensor model with the analytic solution. This is valid in flows with uniaxial symmetry only.

52

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43 69

(uu) in such a flow, calculated by the model, exactly matches the transient of (uu) associated with the full F o k k e r - P l a n c k equation! However, as we mentioned above, a finely tuned version of the deformation tensor model so obtained would only apply to flow in which QSQ retained uniaxial symmetry. In general, QrQ will have three unequal (real) eigenvalues. We explore a fine-tuning of the deformation tensor model that is generally applicable next. 2.2.2. General irrotationalflow Now we turn to the more complicated case of general irrotational flow. In this case the tensor x can be written

~=G

42

0

0

-(41 +42)

.

The analytic solution for the orientation distribution after the decay of initial transients is

1El

1,

~, = ~ exp ~ (4"1 cos 2 cr + 4*2 sin 2 o) cos 2 0 - ~ (4 1 -t- 4* 2) sin 2 0

1

The constant J must be obtained numerically in the case where 4"1 - G41/DR=/:4*2 - G42/DR. Likewise, the m o m e n t tensor (uu),~ corresponding to a given pair (4"1,4*2) must be determined numerically. The moment tensor so obtained is diagonal; in other words the coordinates are in the eigenvector basis of x. To summarize, given the rates of strain (4"1,4*2), one can compute numerically a unique, diagonal (uu)K, which is exact. Now we turn to the deformation tensor model, which also achieves an asymptotic steady state in an irrotational steady flow. Our goal will be to find the best approximation h for the deformation-dependent factor; we shall define the best approximation to mean that value of h to which is associated an equilibrium deformation Q that has a second moment tensor that most closely matches the exact moment tensor (uu),~. Before describing the technique for finding the best h, we must first obtain several preliminary results that will shortly be of use. These results concern (i) the relation between Q, h and ~, and (ii) the space of possible deformations Q. Let the eigenvalues of QVQ be {21, 22, 23}. There are four important relationships between these eigenvalues. First, by virtue of the fact that det Q = 1, 23 = 1/(21 22). The other three relationships between {21, 22, 23} can be obtained from the condition that they are the eigenvalues of QVQ at steady state. At steady state, the coefficient of Q on the right-hand side of (9) must vanish. This yields an independent pair of equations to be satisfied by the steady state eigenvalues of QVQ: 2

1

21

22

2 22

1 21

4"1 2122 = - - h

(24)

~'2 h

(25)

and 2122 --

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

53

In what follows, we shall regard (24) and (25) as parametric representations of the rates of strain (4"1 (h), 4*2 (h)) for a given steady state deformation characterized by (21, ,t2). The second preliminary result concerns what is the space of allowable deformations Q. We have supposed that the scalar function h depends only on the deformation Q. In order that this dependence should be coordinate-invariant, we shall suppose that the dependence is on the principle invariants of Q-rQ. The three "usual" principal invariants are il(Q'rQ)= tr (Q-rQ), i2(QXQ) = l{trZ(Q-rQ)- tr [(QVQ)2]} and i3(Q'rQ)= det (Q-rQ)= 1. Clearly only il and i2 can change; moreover, using the invertible nature of Q and the Cayley-Hamilton theorem, one can easily show that i2(QTQ)= tr (Q-TQ-1). Therefore, we shall suppose that h = h ( & 7), where 6 - tr (QTQ) and 7 =- tr ( Q - T Q - t ) . Now, given a pair (6, 2) one can determine the eigenvalues (21, 22) analytically from 1 ¢~ :

1

1

7 ----771 q- -;- q- 2122 ,

21 q- 22 -~'- - -

2122'

z2

although the expressions are too long to repeat here. However, it is not directly evident what are the physically reasonable values that (cS,7) may take. The physically reasonable values may be established by the following, somewhat indirect argument. Consider the special case of uniaxial symmetry where 21 = 22 : 2. In this special case, one can write 1

c~=22+2-- 5 ,

2

7=~+22.

It is possible to eliminate the variable 2 between these two equations and so to obtain an expression independent of 2: ~b(6, 2) ~- 27 - 18c~7 - 6272 + 4(63

-4- 7 3) :

O.

If we evaluate ~b at a general (6, 7), we find

Clearly, ~b ___0 makes physical sense, with ~b = 0 when 2i = 27 for i :~j. In the (6, 7)-plane, the level set is the wedge-shaped region shown in Fig. 2; the physically reasonable values of (6, 7) lie on the wedge boundary and in the interior of the wedge. This same wedge has surfaced elsewhere in the study of microstructured fluids [1,10], but the equation for the boundary and the proof that the physically sensible points must lie in the interior of the wedge seem not to have been described. Therefore, within the wedge shown in Fig. 2, one can choose a point (6, 7). Corresponding to that point is a pair of eigenvalues (21,).2) that one finds analytically. Corresponding to the e i g e n v a l u e s (21,22) is a second moment tensor (uu)o_ that may be computed using the expressions developed in the Appendix. Thus, what remains is to find a consistent set {6, 7, h, 4"1,4*2 } which satisfies (i) Eqs. (24) and (25) corresponding to the steady state in the deformation tensor model, and (ii) which minimize the "distance" between (uu)K and (uu)o. This leads to the following scheme for the determination of the best h. First we choose a (& 7) in the wedge; this determines (21, 22) and (uu)Q analytically. Next, we guess a value for h. Using the guessed value for h, and the fixed

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

54

10

0 0

5

I0

Fig. 2. The wedge of physically reasonable values of c5 - tr(OTQ) and 7 = tr(O x O - ,). The wedge is bounded by the level set ~b = 0 of the function developed in the text.

(21, 22) , Eqs. (24) and (25) yield the associated values of (6#1, {~*2)" These values of (4'1, e*2)

,,which we compute numerically. Next, we compute the error between ,,;this we define to be the square-root of the sum of the squares of the differences between the (corresponding) two larger eigenvalues of (uu)o. and (uu),,. Essentially, we have the error yield a

for a given h. We simply repeat this procedure and optimize the choice of h until we obtain the m i n i m u m error. F o r each of the 1700 points we c o m p u t e d in this way, we found a single sharp m i n i m u m for the error as a function of h - t h i s is the "best for the point (& 2) in the wedge. Repeating this procedure for m a n y initial choices of (~5, 7) allows us to determine discrete points on the surface h(g, 7) over the wedge. In view of the fact that (24) and (25) are a degenerate pair when 21 = 22 (i.e. when qS(a, 7) = 0), this procedure yields an h that has an error of zero when the flow in uniaxial or biaxial extensional flow--i.e, on the b o u n d a r y of the wedge. The points on the surface h(& 7) are fitted as follows. First, we define convenient coordinates for the wedge, with origin at the cusp:

h"

r / = log 7 - log c$, 1

/t = ~ (log 7 + log c5) -- log 3. In terms of (/~, 1/), the wedge boundaries are given by

rl=rl+(tt)3 +ecosh-'(-8e-3~' - ~ = +~

+~ e*l).

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

55

Along the biaxial edge of the wedge (i'/> 0), a good fit for h is developed using standard methods: hb(//)

=

1+ x

- ~ +/~ -0.417058+0.0517998

log3--~+g

+0.00988129

log3-~+/~

;

along the uniaxial edge of the wedge (1,/< 0), a good fit for h is hu(~)=l+

-5+~

x 0.275366-0.145799

log3-~+/~

+0.0321363

log3-~+/~

.

The two expressions h b and hu are simply an alternative specification of h(2) shown in Fig. 1. Within the wedge, the surface is well approximated by two sheets that meet along a "weld-line" near r/+(~) with equation q = r/w(/~) = - 0.648208/~ + 1.62149/~ 3/2 - 0.540019~ 2. Along the weld-line, a good fit for h is hw(/O -- 1 + 0.735151/~ + 0.419544/~ 2. Now, for r/> r/w, the fit is simply an interpolation between hw and h = hw + r / - r/w ( h b i?+ - r/w

_

hb:

hw).

(26)

Now, for r/< r/w, the fit is an interpolation between hw and hu, plus a small correction: h = hu + ~ - r/ (hw - hu) + (r/w - r/)(r/- 11- ) r/w- r/_ × [0.636280 - 0.104022(r/- r/_) - 0.331413/~].

(27)

The fine-tuning of the model is clearly rather complicated, although straight-forward in detail. However, in order to u s e the model with confidence in any flow, one simply requires (9), (26) and (27).

3. T e s t s o f the d e f o r m a t i o n tensor m o d e l

Now, having determined a suitable approximation h based on analytic solutions for the steady state orientation distribution in steady, irrotational flow fields, we move on to test the model equations thus developed. We shall test the model (9) with (26) and (27) against (i) the exact asymptotic solution in irrotational flows, (ii) against various other closure models, and (iii) against numerical solutions of the full F o k k e r - P l a n c k equation by the double-Lagrangian technique [4,5].

A.J. Szeri, O.J. Lin I J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

56

It is well-known that (3) requires a closure approximation of some sort to eliminate the unknown term ( u u u u ) • E. The simplest is the quadratic closure (referred to compactly as Quad, below) ( u u u u ) • E ~ ( u u ) ( u u ) • E. The next two closures share the feature that they attempt to interpolate between limiting forms that are accurate in weak and in strong flow. They are really designed for use in strong flows where there is a single dominant orientation. The first and second composite closures are due to Hinch and Leal [14] (we refer to these as HL-I and HL-II). Also we consider the hybrid closure of Advani and Tucker [3] (A-T), which is an interpolation between Hand's linear closure [15] and the quadratic closure. A careful listing and the most complete set of comparisons of these closure models in three-dimensional flows has been made by Advani and Tucker [3] and more recently by Cintra and Tucker [16]. They found that Quad performed stably but poorly when the flows are rotational or the orientation distributions are isotropic to moderately aligned. HL-II was much better, except that it appeared to lead to a significantly under-damped model that exhibits strong artificial oscillations in uniform shear flow with weak Brownian effects. HL-I was better in uniform shear, but gave physically unrealistic (negative) diagonal components of (uu) in biaxial extensional flow. A-T performed much like the quadratic closure, but it was a little more accurate overall. Our findings agree on all these points. Where possible, we avoid a great many unnecessary figures by reference to those of Advani and Tucker [3]. Very recently, Cintra and Tucker [16] have described a new type of closure they term an orthotropic closure. This is tested along with another, similar closure (the natural closure) of Verleye and Dupret [17]. Both closures are of the form

(uuuu)

, = S((uu),.,.),

where the functional relationship is compactly written in terms of either the eigenvalues of (uu) (Cintra and Tucker) or the invariants of ( u u ) (Verleye and Dupret). These are equivalent approaches. Both groups fit the function f to either numerical solutions of the full F o k k e r Planck equation (Cintra and Tucker), or to analytical solutions of the Fokker-Planck equation with DR = 0 (Verleye and Dupret). The accuracies of the technique are impressive with tests by Cintra and Tucker [16] demonstrating good agreement of the two closure models with solution of the Fokker-Planck equation in a variety of flows. We remark, however, that the orthotropic and natural closures are not constructed with the idea of guaranteeing physically sensible results in all flows. To aid the reader in drawing comparisons between the present work and the recent paper by Cintra and Tucker [16], we shall make use of the error norm they introduce. The average error is T

~= ~

e(t)dt,

o where the scalar error is e(t)=t~e~eiJ ]

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

57

with e~/= < u u ) ; y " -


With this in mind, our goal is to demonstrate that the deformation tensor model (DTM) can perform at least as well as the other closures we have mentioned. By constructing the model so as to avoid nonphysical behavior, we can be sure that the deformation tensor model will be reliable in other flows, besides the ones we test.

3.1. Tests in uniaxial and biaxial extensional flows First we test the D T M in uniaxial/biaxial extensional flow. Recall that the model was " t u n e d " so as to yield the exact asymptotic value for , subject of course to the error associated with the curve fit for h. In Fig. 3. we show the 1-1 and 3-3 components of together with the analytic asymptotes for mild, moderate and strong uniaxial extensional flows with ~*, = ~'2 - ~ * = - 1 , - 10, and - 100. Time has been nondimensionalized by DR ~. The behavior of the other models is shown in Fig. 3 of Advani and Tucker [3], where the tests are conducted at 4 * = - 2 . 8 8 6 8 and - 2 8 . 8 6 8 . At ~ * = - 1, Quad and A-T overpredict, and HL-I underpredicts the degree of anisotropy. The D T M and HL-II are virtually exact. At stronger values of 4", the other models become more accurate. Next we turn to biaxial extensional flow, in which ~*~ = ~'2 - ~ * is positive. In Fig. 4, we show the results for weak, moderate and strong biaxial extensional flows with 4 * = 1, 10 and 100. The other models are shown in Fig. 4 of Advani and Tucker [3], where the tests are conducted at ~ * = 2.8868 and 28.868. All the models are accurate at 4 * = 1. At ~ * = 10, the D T M is the most accurate closure. Quad and A-T tend to overpredict anisotropy. HL-II slightly underpredicts anisotropy, and HL-I yields nonphysical negative diagonal m o m e n t tensor components. This behavior in HL-I is well known. In strong biaxial extensional flow (~* = 100), the situation is worse for HL-I; the other models are virtually exact.

1.0

~

m

<

' 0.75

>33

/ p / iS

v

0.5

/ I

0.25

i " " " " . . . . . . . . .
01s

Fig. 3. Transients o f the 1-1 a n d 3-3 c o m p o n e n t s of
A.J. Szeri, D.J. Lin / J . Non-Newtonian Fluid Mech. 64 (1996) 43 69

58

O. 5 [ - - - A

0.25

~

.....

33

%%, li

-

-

-

""

i

0

,

0.25

0.5

t Fig. 4. Transients of the 1-1 and 3-3 components of (uu> for the DTM in weak, moderate and strong biaxial extensional flows with 4* = 1, 10 and 100, together with the analytic asymptotes.

To summarize our findings in irrotational flows, we can say the following. Quad and A-T tend to overpredict anisotropy. HL-I tends to underpredict anisotropy in either weakly or moderately extensional flows, and to give non-physical results in strongly extensional flows. HL-II and the D T M are consistently the most accurate in all flows. In terms of the average error, the worst performance of the D T M is ~ = 4 x 10 - 3 when 4 * = 10.

3.2. Tests in shear flow Now we leave irrotational flows, where there is an analytic asymptotic solution for
c, =

-

9/2

0

0

0

E=

I0 i/2

0

0

0

,

where ~ is the shear rate. If we define a new dimensionless time r = t / D R , then the dimensionless P6clet number emerges as Pe = ~ / D R . The P6clet number quantifies the relative influence of the flow and of Brownian motions in determining the orientation distribution. Besides fixing Pe we must also fix the shape factor G; we take G = 1 corresponding to infinite aspect ratio fibers. We begin with a weak shear flow, with Pe = 1. In Fig. 5 we show the 1-1, 2-2 and 1-2 components of (uu> vs. the dimensionless time, for both the D T M and for the numerical solution of the full F o k k e r - P l a n c k equation by the double-Lagrangian technique. The D T M performs very well (e = 7 × 10-4), as does A-T (not shown). Quad, characteristically, overpredicts the degree of anisotropy in this weak flow. The two closures due to Hinch and Leal are virtually exact. In a moderate shear flow (Pe-- 10), less faithful agreement by the D T M is observed in Fig. 6 (e = 2 x 10 2). The results for the other closures are shown in Fig. 1 of Advani and Tucker [3]. The DTM, HL-I and HL-II are all quite accurate, whereas A-T and Quad overpredict anisotropy.

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43 69

59


0.25

I

0

0'5

i'.

Fig. 5. Transients of the 1-t, 2-2 and 1-2 components of for the DTM (dashed curve) in weak shear flow with Pe = 1 and G = 1, together with the numerical solution of the full Fokker-Planck equation (solid curve).

Finally, in a strong shear flow with Pe = 100, one observes from Fig. 7 that the D T M overpredicts anisotropy (i.e. the 1-1 component of ) by about 10% with ~ = 8 × 10 2 However the 2-2 and 1-2 components of are very well approximated. The results for the other closures are shown in Fig. 2 of Advani and Tucker. A-T and Quad each overpredict anisotropy by about 20%. HL-II suffers from fairly large-amplitude (artifical) oscillations about the steady state. Asymtotically, HL-II underpredicts anisotropy by about 10%. HL-I appears to be the best overall at Pe = 100, except for a subtle point to which we earlier made reference. As we have reported elsewhere [5], in a slightly stronger shear flow 2 2 c a n become negative for H L - I - - a pathological situation that corresponds to negative orientational probabilities. By the same analysis that was used in Szeri and Leal [5] one can show that this pathological behavior is possible in Quad, A-T and HL-II, also.

iiiiiiii,iii iiiiiiiii ,iljl

v

0.25

___

0. 375

0/75

Fig. 6. Transients of the 1-1, 2-2 and 1-2 components of
A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

60

s~

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0.75

--

<

>

I

u u 11 0.5

A

V

021 -

_

-

r

-

-

-

B

m

m

~

-

-

_

=

_

-

-

-

. - -

0 5

0.25

t Fig. 7. Transients of the 1-1, 2-2 and 1-2 components of <,u> for the DTM (dashed curve) in strong shear flow with 100 and G = 1, together with the numerical solution of the full Fokker-Planck equation (solid curve).

Pe =

To summarize our results in shear flow, we can say the following. The Quad and A-T closures overpredict anisotropy. The HL-I closure gives nonphysical results in strong shear flow. The HL-II closure suffers from large fictitious oscillations in strong shear flows. The performance of the D T M in shear flow seems about the same as that of the orthotropic smooth closure of Cintra and Tucker [16]; it is not as accurate as the orthotropic fitted closure or the natural closure. 3.3. The limit o f high shear N o w we turn to the behavior of the model in shear flows with large Pe. Asymptotically, as t ~ oo in any steady shear flow, we note that Q(t)Qv(t)---,constant,--although Q(t) remains time dependent. Hence, by the Polar Decomposition Theorem, we write O_.(t) - , c P T ( t ) ,

t - , oo,

where P'rp(t) = I and C is a constant tensor. F r o m (9), we can construct an ordinary differential equation for QQT: dt QQ'r =

x + 3DRh(QTQ)

Q - T Q - 1 _ 5 tr(Q

+QQT{tcT+BDRh(QTQ) I Q - T O - I I r ( Q1

T Q - ~)1]}.

As t ~ oo in any steady shear flow, the left-hand side is zero and Q(t)QT(t)--, B = CC T. This yields 0 = x B + BK T + 6 h D R I - 2hDRtr(B- l)B; we regard this as an equation that determines B. Now, we let z - 6 h / P e asymptotic solution to this equation when G = 1, valid as z--, O:

and find an easy

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 4 3 - 6 9 BI2=Z

-1/5

52+

61

zll/5-t-...

1 8/5 B22 = B33 = z 2/5 - ~ z -~- ... B l l = 2 z - 4 / 5 - ~ z1 2/5 + ~ -7~ z8,5 ' +... and BI3 = B23 = 0. We note that there is only a single solution. Now, we define V = B 1/2 in the usual way through the diagonal basis; thereafter a short calculation leads to

(..)

=

1

vnv,

where H is the eigenvector basis of Q'rQ. Next, we develop asymptotic relations for the (diagonal) components of H using the expressions given in the Appendix: Hll = 6.2831% 4/5 -- 6.00075z 7/5 + 3.76991z 2 + "'" 022 =

7.52814z 1/5 _ 6.28319z 4/5 - 4.89329z 7/5 - 2.722712 2 -k- "'"

H33 =

8.94674z 1/5 -- 6.28319z 4/5 + 3.57869z 7/5 -- 1.67552z 2 + -..

and

Finally, the asymptotic result for the components of ( u u ) is

(uu)ll = 1 - 0.65551z 3/5 - 0.29954z 6/5 + ..(uu)22 = 0.29954(z 3/5 + z 6/5) + ' " (uu)33 = 0.71196z 3/5 - 0.5z 6/5 + ... (uu)12 = 0.79954z 3/5 - 0.97752z 6/5 + " and o f course (uu)~3 = (uu)23 = 0. We checked this asymptotic result against our numerical solution of the D T M and calculated a difference of e = 5 x 10-2 at Pe = 100. We remark that in this high shear-rate analysis we have left h u n s p e c i f i e d - - a l t h o u g h we assume it depends very gradually on Pe. This dependence may be determined by comparison with experiments or numerical solutions (when these are available) to create the best model possible.

4. Nonlinear dynamics of the model in more general fluid flows We have carefully constructed the D T M to give physically sensible results in all flows by forcing det Q = 1. The model has been tested successfully now in the usual irrotational and rotational flows. Finally, we consider the nonlinear dynamics of the model in more general flows. The first question to answer is, can Q ever become infinite? We show that the answer is no; for large IlQll, the Brownian terms in the equation force IlQll to decrease.

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

62

To investigate this, we develop a differential equation for the norm of {2, IIQII:

ld

2 dt

I[o[12= o ' o = GQQT E + 9hDR - hDRJlQ- ljl2l[O[[2.

We use the bound

QQT. E <_ IIQII21Ell, where E~ __ 0 is the maximum eigenvaule of E; hence, ld 2 dt IIQll2 -< allQIHEll + 9hDR - hDRiIQ-lll2][Q][2. This may be rewritten ld 2 dt Ilell= - 9hDR +

[IQI[2(GIE,I- hDRIIQ-ill2).

Finally, it is easy to show that [[Q-'I[2 ~ ~ and also lie-'ll2lIQ[I ~-~ ~ as I[Q[[~-~ o0, by working in components. Therefore, as IIQll-~ o~ the right-hand side tends to - o o . Hence, I[QI[2 always remains finite if DR > 0. In fact, any initial condition that begins at a large enough IIQfl will lead to a decrease in IIQII. This is a definition of a dissipative dynamical system. These important systems have the useful property that there is a unique global attractor for the dynamics [18]. It is interesting to check whether or not the m o m e n t closure models have this same dissipative property and hence a unique global attractor for all DR > 0. We cannot show that Q u i d and HL-I posses the dissipative property unless DR is large e n o u g h - - n o t simply positive. To see this, we begin by taking the contraction of (uu > with the second moment equation (3). This yields, after some manipulation, ld

2 dt

II112= 2 D R - 6D~.ll112+ 2GE.(uu

>2_ 2G(uu > . ( ( u u u u )'E).

If we substitute in Quid, we find ld 2 dt II(uu >11= = 2DR -- 6DRII(uu >112+ 2GE" 2

-

-

2Gll(u u>ll2E . (uu >.

Next, we use the spectrum of E, {El, E2, E3}, with El > E2 > E3, to develop upper bounds for the two terms involving E. This yields ld 2dt

IIll=<2DR +l[(uu>ll2[--6OR + 2G(IEll + IE31)].

Hence, for Quid, we can show the dissipative property if

D~ _>G(IE, I + [U,I).

(28)

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43 69

63

Similar arguments for HL-I lead to the inequality 2at II(uu)II2
E

--6DR+

1

(41E,I+3[E3[) .

Hence, for HL-I, we can show the dissipative property if 2G DR ___- ~ (alE, I + 31E3]).

(29)

An example will help to clarify the meaning of the preceeding argument. In uniform shear flow, E1 = ~/2, E2 = 0 and E3 = - 9~/2. Hence, Quad is dissipative in shear flow when G P e < 3. HL-I is dissipative in shear flow when G P e < 15/7. In other words, when these criteria are satisfied, we know there is a unique global attractor for the dynamics of the closure model. In contrast, for any finite Pc, we know there is a unique global attractor for the DTM. Of course, we have not shown that there are multiple attractors for the moment closure models at moderate or large Pe. We simply cannot prove that there is a unique attractor. We remark that this may be a reflection of the use of the bounds in the preceding arguments rather than the inherent properties of the dynamical system. These arguments concerning the dissipative property are not just of mathematical interest; there are practical consequences of the possible non-dissipative nature of the m o m e n t closure models. Such a consequence may have emerged recently in a study by Chaubal, et al. [19] of liquid crystalline polymers. Part of the study concerned the model with the excluded-volume potential set to zero, identical to the dilute suspension model we consider in the present paper. What Chaubal et al. found was a second, spurious attractor for the dynamics of the HL-I system in simple shear flow. This second attractor was an attractor only for the m o m e n t closure model; there was no corresponding solution for the orientation distribution function. The supurious attractor was shown to attract only a limited part of the space of initial conditions, and to vanish with the addition of a slight perturbation to the flow. Nevertheless, Chaubal et al. point out that such spurious solutions might be associated with numerical instability. It is worth noting that the primary (non-spurious) attractor was qualitatively similar to the one we presented in Section 3.3. Note that when D R = 0, the DTM, which becomes exact, has neither bounded solutions nor the dissipative property. This indicates that in flows that lead to ]IQII- ~ when DR = 0, a very small amount of Brownian diffusion has a very large effect on the orientation distribution, asymptotically in time. A second interesting difference between the D T M and the moment closures concerns volume contraction in phase space. Volume contraction in phase space is associated with the existence of attractors [20]. We will now show that the D T M constantly contracts volumes in phase space whenever DR > 0. However, Quad and HL-I contract volumes everywhere in phase space only when DR is large e n o u g h - - n o t simply positive. To determine volume contraction, we compute the divergence of the right-hand side of (9), regarded as a vector field in Q-space. A lengthy but straight-forward calculation yields the divergence -

3hOR(trO) 2 - 7hORtr(O-TQ-l).

64

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

As this is a striclty negative quantity (for D R > 0), volumes in the phase space are always contracting. Now, we apply the same test to the moment closure models. For this calculation, one must regard 5-dimensional (uu)-space as the state space for the dynamical system. The space is only 5-dimensional because (uu) is symmetric, and also t r ( u u ) = 1. The divergence of the vector field that is the right-hand side of the moment equation with Quad is - 3 0 D R - 12G (uu)" E. Although this looks at first glance negative, it is easy to see that if

00!]

0

1/6

0

(30)

0

then we obtain the divergence of the Quad system to be - 3 0 D R + 60G. Hence if DR is small enough, the divergence can be positive in at least a part of the physically sensible region of the phase space. A similar analysis shows that the HL-I system leads to the divergence -- 30D R - - ~ G

( u u ) " E.

(31)

At the same specific values of (uu) and E, we have - 30D R + 28G for the divergence of the HL-I system. A-T and HL-II are much more complicated. We do not write here the general divergences of the systems, only in the specific flow. For A-T, we obtain - 3 0 D R q-(381/14)G, and for HL-II the divergence is - 30DR + [20/27 + 4352/(63e2)]G. Hence, the commonly employed closure models contract volume everywhere in phase space only if the rotary Brownian diffusivity is large enough. The DTM contracts volume everywhere in phase space if the rotary Brownian diffusivity is simply positive. 5. Conclusions

In this paper, we have developed a new approach to the derivation of a practical model of the evolving orientation distribution in a dilute, Brownian suspension of orientable particles. Beginning from kinetic theory, we assumed a restricted class of deformations that the orientable particles associated with a material point can undergo. Thereafter, we developed an evolution equation for the finite number of degrees of freedom of the deformation, working from the Fokker-Planck equation for the full distribution function. The result is a robust model for the evolving microstructure in a flowing Brownian suspension that is as simple to use as a moment closure model, but is guaranteed to give physically sensible results in all flows as a consequence of the restricted class of allowable deformations. This situation contrasts strongly with moment closure models, which are known to misbehave even in very simple flow fields of primarily theoretical interest. In an extensive series of tests of the deformation tensor model against (i) analytic asymptotic solutions of the full Fokker-Planck equation in irrotational flows, (ii) other closure models, and

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

65

(iii) numerical integration of the full Fokker-Planck equation, we found the new model to be among the most accurate. As it is the only closure model that has been developed in such a way that it always gives physically sensible results, it should be a safe model to use in numerical calculations. Finally, we considered the nonlinear dynamics of the deformation tensor model. It was shown that the deformation is bounded in any unsteady, three-dimensional flow of the suspension, provided the rotary Brownian diffusivity is positive. In a steady, three-dimensional flow, there is a unique global attractor for the deformation associated with a material point, provided the rotary Brownian diffusivity is positive. In contrast, we could not show the existence of a unique global attractor for the commonly employed moment closure models, unless the Brownian diffusivity is large enough compared to the strength of the flow.

Acknowledgements A.J.S. would like to thank C.V. Chaubal and L.G. Leal for pointing out the ingenious paper by Currie. Also, we would like to thank M.C. Altan for providing us with a preprint of his paper on non-Brownian spheroidal particle suspensions. We are grateful to C.L. Tucker for a preprint of his recent work with J.S. Cintra. Finally, we thank S. Hank and E. Titi for helpful discussions.

Appendix A

A.I. The distribution function for the deformation We consider the deformation (7) and determine the associated orientation distribution function using (6). This requires the determinant of the mapping Ou/~uo. Normally when working in a three-dimensional space, the determinant of a linear mapping at a point can be related to a ratio of volumes of two parallelipipeds, whose edges are related through the mapping. However in this case, the mapping preserves lengths in the radial direction, and so we must compare areas of parallelograms instead. To begin, we compute

~u(t, no) ~Uo =- d

Q(t) (no + ESUo) ~=o

de IlQ(t) (no +

uo)ll

1

=

IlQ(t)uoll

( I - uu)Q(t)$Uo.

(A1)

Now, consider an orthonormal triple of vectors in the reference configurtion (no, Vo, Wo). By the tensor Q, these three vectors are mapped to the vectors (Quo, Qvo, Qwo). The volume of the associated parallelipiped is det Q. Now we divide through the orthonormal triple by IIQuoll and obtain the vectors (Quo/llQuollu, Q,,o/llQuoll, Qwo/llQuol]). The volume of the associated parallelipiped is simply det Q/llQuoll3. Now, the triple of vectors (no, u0 + Ero, Uo + eWo) is mapped by the deformation (7) to u, u +

(1- ..)Qvo, u + c

(!-

uu)Qwo .

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

66

This triple encloses a parallelipiped with volume E2 det Q/llQuoll 3. The volume of the parallelipiped in the reference configuration is simply e2. Because Uo and u are both unit vectors, this implies that the ratio of the areas of the parallelograms shown in Fig. A1 is det Q/llQuo[[3. This leads directly to (10).

A.2. Practicalformulae for (uu) in terms of Q Now we describe how to obtain the second m o m e n t of the orientation distribution function in terms of the deformation given by Q. The goal is to write (uu) exactly in terms of Q. First, using (6) and (7), it is useful to write (uu) as an integral over the reference configuration, which is here assumed to be isotropic:

1 fQuoQuo f uu ~(u't)d2u= ~ J IQuol 2 d2uo. It is possible to pull Q from the numerator of the integrand, and so to write

(u.)

_1 Q [ ~ UoUoJ2 4re

~'~,T _ ~1QHQT"

~J~aUo)~

NOW the tensor H depends only on H = j-

UoUo

=

QTQ; this is apparent upon rewriting in the form

d2uo .

Uo" QTQuo Because QTQ is a symmetric, positive tensor, it may be written in a basis of orthogonal eigenvectors in the form (20). In the same eigenvector basis of QTQ, the tensor H is also diagonal. Therefore, if P is the orthogonal tensor which diagonalizes QTQ, one easily obtains (in general components)

(uu) = ~1 QPHPTQT. In this expression, H is the diagonal tensor obtained A rather involved calculation yields the diagonal spherical coordinates based on the eigenvectors (cos 0 cos a, cos 0 sin a, sin 0) where - zc/2 _< 0 _< ~/2 nent of H is integrated to yield (for 2~ ¢ 22)

in what follows. components of H. Written in modified of Q~Q, in terms of which u0 = and 0 _< a _< 2zc, the first diagonal compo-

1 i

Reference Configuration

~

(t uu) 0~,

Deformed Configuration

Fig. A1. A diagram of the transformation of a differential area from the reference to the deformed configurations of the sphere of orientations.

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

Hi,1 4zc

67

E(qS/m)Al/222

1 21 -- ~2

(/~1 -- /~2)(fl~12fl~2-- 1) 1/2'

where E is the elliptic integral of the second kind with arguments ~b(#l, ,~2) = arcsin

[( 1 - A#2---~/ ' ~|/2jl

and m(21,22) =

2 2 /~122 -- 21

2 2

The second and third components H2.2 and H3,3 are most easily obtained by solving the following equations: tr
H3,3 = 4~z

and 4zcF(q9 Ira) t r H = H1,1 + H2,2 +//3,3 = (2122 - 2i- ])1/2 • Here F is the elliptic integral of the first kind. The elliptic integrals are very efficiently computed using the methods of Carlson [21]. As an illustration, we offer the following. We time the integration of the DTM, including the determination of
4re HS

4re

/~2

2( - 1 + 23) 2( - 1 + 2 3 ) 3/2 25 27/2 log[23/2 + ( - 1 + 23) 1/2]

3,3 = ,,;[2"Jr-2"-----~ "1

4~

21/2 1og[~ 3/2 + ( __ 1 + 23) 1/2]

1-

(-1+

'

23)3/2

A.3. Practicalformulae for @uuu) in terms of Q Because the constitutive equation for fiber suspensions includes a dependence on the fourth moment of the orientation distribution function, we briefly describe how this may be calculated efficiently. As in computing (uu), we work in the eigenvector basis for Q'rQ. We begin by writing
1[


d2.0

where u ° are the components of Uo. Because QTQ can be written in the form (20), there are in general 21 non-zero components (of the original 8 1 ) - - o f which only 6 are independent:

< HUUU>1111, < H/g/H/> 2222,3333, 1122, < HHH//>1133, < UUHH> 2233-

A.J. Szeri, D.J. Lin /J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

68

The other 15 may be obtained from well-known symmetries, e.g.

("")1212 = Also, we have the trace property

(UUUU)~kk= (UU) a, which yields the system

(,,,,..),,1, + ( . . , , . ) , , 22 + ( ..,,,,),,33 = ( . . ) , , . (

+ (

+

= (

The right-hand sides of these equations are known from the calculation of system solvable, we integrate to obtain

(uu). To make this

(.,,..),,,1, by Gaussian quadrature, for example. By performing only three integrals over the reference configuration, and by use of the trace and symmetry properties we can generate the components of (uuuu) efficiently.

References [1] R.B. Bird, C.F. Curtiss, R.A. Armstrong and O. Hassager, Dynamics of Polymeric Liquids, Vol. 2, Kinetic Theory, John Wiley & Sons, New York, 1987. [2] R.G. Larson, Constitutive Equations for Polymer Melts and Solutions, Butterworth, 1988. [3] S.G. Advani and C.L. Tucker, Closure approximations for three-dimensional structure tensors, J. Rheol., 34 (1990) 367-386. [4] A.J. Szeri and L.G. Leal, A new computational method for microstructured fluids. Part 1. Theory, J. Fluid Mech., 242 (1992) 549-576. [5] A.J. Szeri and L.G. Leal, A new computational method for microstructured fluids. Part 2. lnhomogeneous shear flow of a suspension, J. Fluid Mech., 262 (1994) 171-204. [6] A.J. Szeri, A deformation tensor model for liquid crystalline polymers, J. Rheol., 39 (1995) 873-891. [7] F. Folgar and C.L. Tucker, Orientation behavior of fibers in concentrated suspensions, J. Reinf. Plas. Compos., 3 (1984) 98-119. [8] F.P. Bretherton, The motion of rigid particles in a shear flow at low Reynolds number, J. Fluid Mech., 14 (1962) 284-304. [9] A.J. Szeri and L.G. Leal, Microstructure suspended in three-dimensional flows, J. Fluid Mech., 250 (1993) 143-167. [10] P.K. Currie, Constitutive equations for polymer melts predicted by the Doi-Edwards and Curtiss-Bird kinetic theory models, J. Non-Newtonian Fluid Mech., 11 (1982) 53-68. [11] N. Malamataris and T.C. Papanastasiou, Closed-form material functions for semidilute fiber suspensions, J. Rheol., 35 (1991) 449-464. [12] B.N. Rao, L. Tang and M.C. Altan, Rheological properties of non-Brownian spheroidal particle suspensions, J. Rheol., 38 (1994) 1-17. [13] A.J. Szeri, Pattern formation in recirculating flows of suspensions of orientable particles, Philos. Trans. R. Soc. London, Ser. A, 345 (1993) 477-506. [14] E.J. Hinch and L.G. Leal, Constitutive equations in suspension mechanics. Part 2. Approximate forms for suspension of rigid particles affected by Brownian rotations, J. Fluid Mech., 76 (1976) 187-208. [15] G.L. Hand, A theory of anisotropic fluids, J. Fluid Mech., 13 (1962) 33-46.

A.J. Szeri, D.J. Lin / J. Non-Newtonian Fluid Mech. 64 (1996) 43-69

69

[16] J.S. Cintra and C.L. Tucker, Orthotropic closure approximations for flow-induced fiber orientations, J. Rheol., 39 (1995) 1095-1122. [17] V. Verleye and F. Dupret, Prediction of fiber orientation in complex injection molded parts, Proc. ASME Winter Annual Meeting, 28 Nov.-3 Dec., 1993, New Orleans. [18] J.K. Hale, Asymptotic Behavior of Dissipative Systems, American Mathematical Society, Providence, RI, 1988. [19] C.V. Chaubal, L.G. Leal and G.H. Fredrickson, A comparison of closure approximations for the Doi theory of LCPs, J. Rheol., 39 (1995) 73-103. [20] J.M.T. Thompson and H.B. Stewart, Nonlinear Dynamics and Chaos, J. Wiley & Sons, Chichester, 1986. [21] B.C. Carlson, Computing elliptic integrals by duplication, Numer. Math., 33 (1979) 1-16.