Universal viscometric functions for dilute polymer solutions

J. Non-Newtonian Fluid Mech., 71 (1997) 245 – 272

Universal viscometric functions for dilute polymer solutions J. Ravi Prakash a, Hans Christian O8 ttinger b,* a

b

Department of Chemical Engineering Indian Institute of Technology, Madras 600 036, India ETH Zu¨rich, Department of Materials, Institute of Polymers, Swiss FIT Rheocenter, CH-8092 Zu¨rich, Switzerland Received 13 January 1997; accepted 27 January 1997

Abstract The Gaussian approximation has been shown previously to be excellent for the treatment of hydrodynamic effects in dilute polymer solutions. However, the computational time required to find the viscometric functions in simple shear flow is prohibitively long for chains with a large number of beads. Here we introduce a new approximation which retains the accuracy of the Gaussian approximation but is significantly less computationally intensive. Thus the rheological behavior of long chains may be explored. Extrapolation of results obtained numerically for long chains to the infinite chain length limit is shown to lead to predictions independent of model parameters. As a result, within the context of the approximation introduced here, universal viscometric functions for dilute polymer solutions in simple shear flow under theta conditions are obtained. © 1997 Elsevier Science B.V. Keywords: Polymer solutions; Viscometric functions; Gaussian approximation

1. Introduction In order to predict the behavior of polymer solutions in complex flows, it is essential to supplement the conservation laws of continuum mechanics with a constitutive equation for the stress tensor. Constitutive equations are frequently obtained by developing a kinetic theory for the polymer solution. Such theories typically represent the polymer molecule by a coarse-grained model and replace the rapidly fluctuating motions of the solvent molecules surrounding the polymer molecule by a randomly varying force field. The usefulness of a constitutive equation developed in this manner is verified by comparing the prediction of viscometric functions in certain simple flows, such as simple shear flows or elongational flows, with experimental results. It is a striking feature of polymer solutions that many properties, including viscometric functions, exhibit universal behavior when the experimental data is plotted in terms of * Corresponding author. Fax: 41 1 6321076 0377-0257/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII S 0 3 7 - 0 2 5 7 ( 9 7 ) 0 0 0 1 2 - 8

246

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

appropriately normalized coordinates. For example, when the viscosity of a dilute solution is plotted as a function of shear rate, the various curves for different types of monomers can be superimposed into a single curve by normalizing the viscosity and shear rate for each type of monomer. Such universal behavior is believed to arize because the macroscopic properties of the polymer solution are determined by a few large scale properties of the polymer molecule. Thus solutions that differ from each other with regard to the chemical structure or molecular weight of the dissolved polymers, the temperature etc., still behave similarly so long as a few parameters specifying the molecular characteristics are the same. Molecular theories developed within the framework of polymer kinetic theory usually cite this universal behavior of different polymer systems in order to justify the introduction of crude mechanical models to represent real polymer molecules. For instance, a frequently used model to represent a flexible polymer molecule is a linear chain of identical beads connected by Hookean springs. While such models do indeed represent certain large scale features of the polymer molecule, such as its stretching and orientation by the solvent flow field, it turns out that apart from a basic length and time scale, there occur other parameters that need to be specified, for example, the number of beads N in the chain, the strength of hydrodynamic interaction h*, the finite spring extensibility parameter b etc. Though one can argue that the values of some of these parameters lie within a fairly restricted range, their ultimate choice is essentially arbitrary, and this choice crucially influences the model’s predictions. In other words, the predictions of these models are by no means universal. It is the purpose of this paper to introduce an approximate treatment of the dynamic behavior of a polymer solution within whose context it is possible to transcend the limitations of a mechanical polymer model and verify that universal predictions of viscometric functions are indeed obtained. An early attempt at obtaining a constitutive equation with kinetic theory was the model of Rouse [1] which depicts a dilute polymer solution as a Newtonian solvent with non-interacting Hookean bead-spring chains suspended in it. The solvent is believed to influence the motion of the beads by exerting a drag force and a Brownian force. This model leads to a constitutive equation that is a multimode generalization of the convected Jeffreys or Oldroyd B model. As is well known, while this constitutive equation accounts for the presence of viscoelasticity through the prediction of a nonzero first normal stress difference in simple shear flow, there are several other features of dilute solution behavior that it does not predict. For instance, the nonvanishing of the second normal stress difference, the shear rate dependence of the viscometric functions, and the molecular weight dependence of the viscosity, the first normal stress difference and the diffusion coefficient. Significant progress was made when Zimm [2] incorporated the effect of hydrodynamic interaction in an equilibrium a6eraged form into the theory of Rouse. Hydrodynamic interaction is a long range interaction between the beads which arizes because of the solvent’s propensity to propagate one bead’s motion to another through perturbations in its velocity field. While it is reasonably straight forward to generalize the Rouse theory to include hydrodynamic interaction, it renders the resultant equations analytically intractable. Zimm introduced the assumption of equilibrium averaging in order to solve them approximately. The predictions by the Zimm theory of the molecular weight dependence of the viscometric functions, the diffusion coefficient and the relaxation time, and of small amplitude oscillatory shear flow material functions in excellent agreement with experiment. However, as in the Rouse case, it too fails to predict the

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

247

occurrence of a nonvanishing second normal stress difference and the shear rate dependence of the viscometric functions. A nonzero positi6e second normal stress difference and shear rate dependent viscometric functions were predicted when O8 ttinger [3,4] accounted for the presence of hydrodynamic interaction through a self-consistent averaging procedure. It must be noted that the sign of the second normal stress difference has not been experimentally conclusively established. However, an exact solution obtained by Brownian dynamics simulations leads, at low shear rates, to a negative value for the second normal stress difference. Furthermore, while the predictions of the shear rate dependence of the viscosity and first normal stress difference by the self-consistent averaging procedure are in qualitative agreement with the Brownian dynamics simulations, they do not agree quantitatively. The introduction of the Gaussian approximation for the hydrodynamic interaction [5–8] led to predictions of viscometric functions which are in close agreement with the results of Brownian dynamics simulations. Unfortunately, the computational time required to examine chains with a large number of beads is prohibitively large. Consequently, it is impossible to obtain results for very long chains with what is by far the best approximation to date. There are two principal reasons why it is important to examine long chain behavior. Firstly, the representation of more degrees of freedom permits the exploration of more length (and time) scales. As a consequence, aspects of polymer solution behavior that might otherwise be hidden when only short chains are considered are revealed. Magda et al. [9] and Kishbaugh and McHugh [10] cite this as their reason for introducing a decoupling approximation which while retaining the accuracy of the self-consistent averaging procedure is much more computationally efficient. This enables them to examine the behavior of chains with a large number of beads. There is a second and more fundamental reason for considering long polymer chains. O8 ttinger has shown, both in the generalized Zimm model [3] and in the zero shear rate limit of the Gaussian approximation [5], that in the limit N “ , results are obtained which are parameter free. Thus predictions of rheological properties in this limit represent general consequences of the assumptions regarding hydrodynamic interaction and are independent of the details of the mechanical model. In this paper we introduce an approximation, which we call the twofold normal approximation, that retains the accuracy of the Gaussian approximation but is not so computationally intensive. By this means, it becomes possible to first accumulate data for chains with as many as 100 beads, and then to obtain universal predictions of viscometric functions in shear flow by extrapolating these finite chain length results to the limit N“ . The basic equations that govern the Gaussian and twofold normal approximation are introduced in Section 2. Section 3 is devoted to examining linear viscoelastic properties that may be obtained with the twofold normal approximation. After deriving a codeformational memory integral expansion for the stress tensor in Section 3.1, small amplitude oscillatory shear flow is considered in Section 3.2. In Section 3.3, expressions for the viscosity and first normal stress difference in the zero shear rate limit of steady shear flow are obtained. Viscometric functions for steady shear flows obtained by numerical calculations are the subject of Section 4. Details of the numerical calculations are given in Section 4.1, and the results of these calculations for chains of finite length are presented in Section 4.2. In each of these flow situations, a systematic comparison with the results of the Gaussian approximation is carried out. The universal viscometric functions predicted the twofold normal approximation are presented in Section 4.3. Section 5 summarizes the conclusions of this paper.

248

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

2. Basic equations

2.1. Bead-spring chains with hydrodynamic interaction The instantaneous configuration of the bead-spring chain can be specified by N bead position vectors r1, r2,…, rN with respect to a laboratory-fixed frame of reference. The solvent in which the beads are suspended is assumed to have a homogeneous velocity field, i.e., of the form 7 = 70 + k(t) · r, where 70 is a constant vector and k(t) is a traceless tensor. For homogeneous flows, it is reasonable to expect that the configurational distribution function c is independent of the position of the centre of mass, and depends only on the N −1 connector vectors Qk = rk + 1 − rk. The kinetic theory of such a system, with hydrodynamic interaction taken into account, requires that c satisfy the ‘diffusion’ equation [11]: N−1 ( (c · =− % (t j = 1 (Qj

!

k · Qj −

"

H N−1 kBT N − 1 ( ( % A0 jk · Qk c + · A0 jk · c % z k=1 (Qk z j,k = 1 (Qj

(1)

where H is the spring constant, z is the bead friction coefficient (for spherical beads with radius a, in a solvent with viscosity hs, z is given by the Stokes expression: z=6phsa), kB is Boltzmann’s constant, and T is the absolute temperature. The tensor A0 jk accounts for the presence of hydrodynamic interaction, and is defined by A0 jk = Ajk 1 + z(Vj,k + Vj + 1,k + 1 − Vj,k + 1 −Vj + 1,k )

(2)

where 1 is the unit tensor, Ajk is the Rouse matrix, given by

Á 2 for j− k = 0 Ajk = Í −1 for j− k = 1 Ä 0 otherwise

(3)

and the tensor Vmn, which describes the hydrodynamic interaction between beads m and n, is here assumed to be the given by the Oseen-Burgers expression:





Á rmnrmn à 1 1+ 2 , r mn = rm −rn for m"n Vmn = Í8ph srmn r mn Ã Ä 0 for m=n

(4)

The rheological properties of a dilute polymer solution may be obtained by calculating the polymer contribution to the stress tensor, t p, which for bead-spring chains is given by the Kramers expression [11]: N−1

t p = − nH % ŽQjQj  + (N− 1) nkBT1.

(5)

j=1

Here, n is the number density of polymers, and angular brackets represent averaging with respect to the distribution function c.

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

249

It is clear from Eq. (5) that our purpose is essentially served once the second moments ŽQjQj  are calculated. The time evolution of the average of any arbitrary quantity can be obtained from Eq. (1). In particular, by multiplying Eq. (1) by QjQk and integrating over all configurations, the following evolution equation for the connector vector covariances is obtained: 2kBT d ŽA0 jk  ŽQjQk  =k · ŽQjQk  + ŽQjQk  · k T + z dt −

H N−1 % [ŽQjQm · A0 mk  +ŽA0 jm · QmQk ] z m=1

(6)

The central problem of molecular theories that incorporate hydrodynamic interaction, and which attempt to predict the rheological properties of dilute polymer solutions, is to find the solution of Eq. (6). It is not a closed equation for the second moments since it involves more complicated moments on the right hand side. Various treatments of hydrodynamic interaction essentially boil down to finding a suitable closure approximation. The best approximation so far is the Gaussian approximation, which is discussed in the section below.

2.2. The Gaussian approximation The pre-averaging assumption of Zimm, and the self consistent averaging method of O8 ttinger solve the closure problem by replacing the tensor A0 jk with an average. This leads to the neglect of fluctuations in the hydrodynamic interaction. The Gaussian approximation [5–7] makes no assumption with regard to the hydrodynamic interaction. It only assumes that the solution of the diffusion Eq. (1) may be approximated by a Gaussian distribution. This special choice for the configurational distribution function c, makes Eq. (6) a closed equation for the second moments since all complicated averages on the right hand side can be reduced to functions of the second moment with the help of the Gaussian distribution. Before stating the result of this reduction process, it is worth recalling that a Gaussian distribution is completely characterized by the covariances sjk, defined by sjk ŽQjQk 

(7)

The tensors sjk are not symmetric, but satisfy the relation sjk = s Tkj

(8)

After a great deal of algebra, important details of which may be found in [5], one can rewrite Eq. (6) as the following evolution equation for the covariances sjk : d 2kBT H N−1 sjk = k · sjk + sjk · k T + % [sjm · A( mk +A( jm · smk ] A( jk − dt z m=1 z −

H H N−1 % [sjl · Glp,mk :spm +smp :Glp, jm · slk ] z kBT m,l,p = 1

where, the (N− 1)× (N− 1) matrix with tensor components, A( jk, is defined by

(9)

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

250

A( jk = Ajk 1 + 2h*



H(sˆ j,k )+ H(sˆ j + 1,k + 1)

H(sˆ j,k + 1)

H(sˆ j + 1,k )

n

(10)

j − k

j −k−1 j −k+1 and the (N − 1)2 × (N− 1)2 matrix with fourth rank tensor components, Glp, jk, is defined by Glp, jk =

−

−



3 2h* u( j, l, p, k)K(sˆ j,k )+ u( j +1, l, p, k +1)K(sˆ j + 1,k + 1) 4

j −k 3 u( j, l, p, k + 1)K(sˆ j,k + 1)

n

u( j+1 , l, p, k)K(sˆ j + 1,k )

(11)

j− k −1 3

j−k +1 3 In Eq. (10) and Eq. (11), h*= a(H/pkBT)1/2 is the hydrodynamic interaction parameter, the tensors sˆ mn are defined by −

−

1 H max(m,n) − 1 % sjk, m− n kBT j,k = min(m,n)

sˆ mn = sˆ Tmn = sˆ nm =

&

(12)





(13)





(14)

and the two functions of the second moments, H(s) and K(s) are given by H(s)=

3 2(2p)3/2

and K(s) =

−2 (2p)3/2

&

dk

dk

1 kk 1 exp(− k · s · k) 2 1− 2 k k 2

1 kk 1 k exp(− k · s · k). 2 k 1− 2 k k 2

The function u( j, l, p, k) is unity if l and p lie between j and k, and zero otherwise: u( j, l, p, k) =

!

1 0

if j5 l,pBk otherwise

or

k 5l,p Bj

(15)

Note that the conventions H(sˆ jj )/0= 0 and K(sˆ jj )/0=0, have been adopted in Eq. (10) and Eq. (11), respectively. Both the hydrodynamic interaction functions H(s) and K(s) can be evaluated analytically in terms of elliptic integrals. These and many other properties of these functions are discussed extensively in [4,5,12,8,13]. Eq. (9) constitutes a system of (N− 1)2 coupled ordinary differential equations for the components of the covariance matrix sjk. Their solution permits the evaluation of the ‘Gaussian approximation’ to the stress tensor: N−1

t p = − nH % sjj + (N− 1) nkBT1

(16)

j=1

As a consequence, all the material functions of dilute polymer solutions in arbitrary homogeneous flows may be found. Small amplitude oscillatory shear flows and steady shear flow in the limit of zero shear rate have been examined by O8 ttinger [5] for chains with N530 beads. Zylka [8] has obtained the material functions in steady shear flow for chains with N 515 beads and compared his results with those of Brownian dynamics simulations (the comparison was made for chains with N = 12 beads). Of all the approximate treatments of hydrodynamic interaction introduced so far, the Gaussian approximation compares best with simulation results. However, for both the flow situations just discussed, examining chains with larger values of N becomes increasingly difficult because of computational intensity. The precise reasons for this computational intensity will be discussed in greater detail in the subsequent sections.

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

251

It is appropriate at this point to note that the equations that govern the self consistently averaged theory can be obtained by dropping the last term in Eq. (9). In other words, only the matrix A( jk arises due to the presence of hydrodynamic interaction, and the function H(s) is sufficient for its evaluation. The solution of this truncated version of Eq. (9) can be coupled with the expression for the stress tensor Eq. (16) to evaluate all the rheological properties predicted by the consistent averaging approximation. Viscometric functions in simple shear flows were obtained by O8 ttinger [3] for chains with N 525 beads; longer chains required excessive computational time. The decoupling approximation [9,10] overcomes the computational intensity of the consistent averaging approximation and is consequently capable of predicting the behavior of long chains. However, as has been mentioned earlier, results obtained with the consistent averaging approximation are not in agreement with Brownian dynamics simulations. The question then naturally arises whether it is possible to develop an approximation scheme that is not as computationally intensive as the Gaussian approximation, but is just as accurate. This is the central question that this paper addresses. We find that it is indeed possible to develop such an approximation, and exploit the far reaching consequences of its development. Before describing the new approximation, we briefly consider the decoupling approximation in the next section since the twofold normal approximation was originally inspired by it, and shares some of its features. As we shall see shortly, the twofold normal approximation is more general than the decoupling approximation, and does not suffer from some of its shortcomings.

2.3. The decoupling approximation The tensor A0 jk is responsible for coupling the connector vectors to one another in the diffusion equation (1), and for its nonlinear character. Linearity is achieved in the Rouse model by neglecting hydrodynamic interaction and setting A0 jk =Ajk, the constant Rouse matrix. The Zimm model replaces A0 jk by its equilibrium average, A0 jk, which is called the modified Rouse matrix, and is given by A0 jk = Ajk + 2h*



2

j −k

−

1

j− k−1

−

1



j −k+1

(17)

The self consistent averaging method replaces A0 jk with A( jk (defined in Eq. (10)), which is an average carried out with the non-equilibrium distribution function. All these linearizing procedures render the diffusion equation solvable, indeed the solution c is a Gaussian distribution in each case. However, they do not remove the coupling between the connector vectors. Decoupling can be achieved rigorously both in the Rouse and in the Zimm model by mapping the connector vectors, Q1,…, QN − 1 to a new set of ‘normal’ coordinates, Q %1,…, Q %N − 1 with the transformation: N−1

Qj = % PjkQ %k,

(18)

k=1

where, Pjk are the elements of an orthogonal matrix with the property (P − 1)jk = Pkj,

(19)

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

252

such that N−1

% PmjPmk = djk.

(20)

m=1

Here, djk is the Kronecker delta. In the Rouse model, the Rouse orthogonal matrix diagonalizes the Rouse matrix Ajk (see Eq. (3)). The elements of the Rouse orthogonal matrix can be evaluated analytically by the expression, Pjk =

'  

2 jkp . sin N N

(21)

On the other hand, the Zimm orthogonal matrix which diagonolizes the modified Rouse matrix: N−1

% PjiA0 jkPkl = a˜ldil

(22)

j,k = 1

must be found numerically for N\ 4. Here, a˜l are the so called Zimm eigenvalues. The diffusion equation for the Rouse and the Zimm model, in terms of these normal coordinates, admits a solution for the configurational distribution function of the form N−1

c(Q %1,…, Q %N − 1)= 5 ck (Q %k ).

(23)

k=1

As a consequence, it becomes uncoupled and leads to (N −1) diffusion equations, one for each of the ck (Q %k). Since the Q %k are independent of each other, all the covariances ŽQ %Q j %k with j "k %Q % survive. Differential equations that govern are zero, and only the (N− 1) variances s % ŽQ j j j the variances s %j are then derived in a manner entirely analogous to the derivation of Eq. (6) for the covariances ŽQjQk . The stress tensor and consequently all the rheological properties may be obtained once the evolution equations that govern the variances s %j are solved. The decoupling approximation introduced by Magda et al. [9] and Kishbaugh and McHugh [10] (who use FENE springs in place of the Hookean springs of Magda et al.) consists of carrying over this ‘diagonalize and decouple’ procedure to the case of the self consistently averaged theory. Firstly, they map the connector vectors Qj to a new set of coordinates Q %j using the time-invariant Rouse orthogonal matrix Pjk (see Eq. (21)). (Kishbaugh and McHugh also use the Zimm orthogonal matrix). Secondly, and this is the key assumption, they assume that the same orthogonal matrix Pjk diagonalizes the matrix of tensor components A( jk. The result of this assumption is that even in the self consistently averaged theory, the diffusion equation can be solved by the method of separation of variables (see Eq. (23)). Thus, as in the Rouse and Zimm case, only the (N − 1) transformed coordinate variances s %j are non-zero, and differential equations governing these variances can be derived by manipulating the uncoupled diffusion equations. Finally, they solve these equations numerically in order to obtain the appropriate material functions. Kishbaugh and McHugh [10] claim that the decrease in the number of differential equations to be solved (from (N− 1)2 for the covariances sjk to (N−1) for the variances s %) j is responsible for the great reduction in computational time achieved by the decoupling approximation. As we shall see in Section 4.1, this argument is incomplete; the real reason for the reduction lies elsewhere. Furthermore, it turns out that it is unnecessary to obtain evolution equations for the

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

253

variances s %j through a multistage procedure that centres around the assumption that Pjk diagonalizes the matrix A( jk. The same results can be obtained with a less stringent assumption and fewer inconsistencies as demonstrated by the twofold normal approximation below.

2.4. The twofold normal approximation The approximation introduced in this paper consists of the following three steps: 1. The configurational distribution function c is assumed to be a ‘normal’ distribution, as in the Gaussian approximation. Since the hydrodynamic interaction tensor A0 jk in the diffusion equation is not replaced by an average, fluctuations in the hydrodynamic interaction are accounted for. 2. The Rouse or the Zimm orthogonal matrix Pjk is used to map the bead connector vectors Qj to a new set of ‘normal’ coordinates Q %.j 3. The chosen orthogonal matrix is assumed to diagonalize the covariance matrix sjk ; i.e., N−1

% Pjpsjk Pkq = ŽQ %pQ %q = s %pdpq.

(24)

j,k = 1

The twofold normal approximation derives its name from the two very different assumptions of ‘normality’ made above. By making use of Eq. (24), the time evolution equations for the variances s %j can be obtained from Eq. (9). They have the following form: 2kBT H H H N−1 d T % [s %j · Djk :s %k +s %k:Djk · s %]j k · s % +s % · k + − [s % · L +L · s %] − L s %= j j j j j j j j z z kBT k = 1 z dt (25) where Lj L0 jj are the diagonal tensor components of the matrix L0 jk given by N−1

L0 jk = % PljA( lp Ppk.

(26)

l,p = 1

The matrix Djk in Eq. (25) is given by Djk =

N−1

%

Plj Ppk Glp,mn Pmj Pnk.

(27)

l,m,n,p = 1

Note that the tensors Djk have the useful symmetry property, Djk =Dkj. this follows from the fact that Glp,mn is unchanged when the indices l and p are interchanged, and when m and n are interchanged. In Eq. (26) and Eq. (27), the tensors A( jk and Glp,mn are given by Eq. (10) and Eq. (11), respectively. However, the argument of the hydrodynamic interaction functions is now given by sˆ mn =

1 H max(m,n) − 1 N − 1 % % Pjm Pkms %m. m −n kBT j,k = min(m,n) m = 1

(28)

One may recall that the ultimate purpose of the assumption that Pjk diagonalizes A( jk made in the decoupling approximation, was to obtain (N−1) equations for the variances s %.j Here, the same result has achieved without making any attempt to uncouple the diffusion equation.

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

254

Indeed the approximation of Magda et al. and Kishbaugh and McHugh (for the case of Hookean springs) is recovered from the twofold normal approximation when the last term in Eq. (25), which accounts for fluctuations in hydrodynamic interaction, is dropped. In other words, the two different routes for finding governing equations for the quantities s %j lead to the same result. As a consequence, material functions in shear flow found from the twofold normal approximation in the absence of fluctuations are identical to those found by Magda et al. and Kishbaugh and McHugh. The reason both approximations have the same governing equations is because only the diagonal components Lj appear in Eq. (25). While the decoupling approximation assumes that off-diagonal elements of the matrix L0 jk are zero, no such assumptions is made in the twofold normal approximation. As a result, only material properties that also require off-diagonal components of the tensor L0 ik would show up the difference between the two approximations. Computations in simple shear flow reveal that the off-diagonal elements are not close to zero; infact, they are of the same order of magnitude as the diagonal elements. The decoupling approximation is therefore inconsistent in making the assumption hat the matrix L0 jk is diagonal, since calculations within its framework show that this is not true. It is also interesting to note that the implementation of the decoupling approximation by Kishbaugh and McHugh does not fully exploit the diagonalization assumption. In Section 4.1 we discuss how proper exploitation would actually lead to a further optimisation of their code (but to probably rather inconsistent results). Eq. (25) constitutes (N− 1) equations for the variances s %.j Their solution permits the calculation of the ‘twofold normal approximation’ to the stress tensor: N−1

t p = − nH % s %j + (N− 1) nkBT1,

(29)

j=1

and, consequently, all the material functions in arbitrary homogeneous flows.

3. Linear viscoelastic properties

3.1. Codeformational memory integral expansion In this section, a first order codeformational memory integral expansion for the stress tensor will be derived with the help of Eq. (25) and Eq. (29). Linear viscoelastic properties, such as the material functions for small amplitude oscillatory shear flow, and the steady shear viscosity and first normal stress difference in the limit of zero shear rate, can be obtained from such an expression. We begin by expanding the tensor s %j in terms of deviations from its isotropic equilibrium solution up to first order in velocity gradient: s %j =

kBT [1+ ej + …] H

(30)

Substituting this expression into Eq. (25) leads to the following system of ordinary differential equations for the first order term ej :

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

1 N−1 d % bjkek, ej = k +k T − lH k = 1 dt

255

(31)

where lH = (z/4H), the matrix bjk is given by 1 3 2h* ajk bjk = ljdjk − 2 20

(32)

with N−1

lj = % PmjA0 mn Pnj

(33)

m,n = 1

(note that lj = a˜j (see Eq. (22)) only if the Zimm orthogonal matrix is used to perform the diagonalization of sjk ) and the matrix ajk defined by ajk =

N−1

%

Plj Ppk Pmj Pnk

l,m,n,p = 1

−

u(m+ 1, l, p, n)

m −n+ 1

3



u(m, l, p, n)+u(m +1, l, p, n+1)

n

m −n 3

−

u(m, l, p, n+1)

m−n−1 3

.

(34)

In order to derive Eq. (31) from Eq. (25), the equilibrium values of the hydrodynamic interaction functions H and K are sufficient. Detailed expressions for expansions around an isotropic tensor can be found in [5]. Upon integrating Eq. (31) with respect to time, and substituting the solution in Eq. (29), we obtain the following expression for the first order codeformational memory integral expansion: t p = − nkBT

&

t

−

ds G(t −s)g[1](t, s),

(35)

where g[1] is the codeformational rate-of-strain tensor [14] and the memory function G(t) is given by N−1



G(t) = % E − j,k = 1

t b lH

n

.

(36)

jk

The exponential operator E[M] maps one matrix into another according to E[M]jk = djk + Mjk +

1 N−1 % MjmMmk +…, 2! m = 1

(37)

and has the useful property: N−1 N−1 d E[Mt]jk = % Mjm E[Mt]mk = % E[Mt]jmMmk dt m=1 m=1

(38)

The memory function G(t) can also be written as a superposition of exponential functions: N−1

G(t)= % Hj e − t/tj, j=1

(39)

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

256

where the set {(ti, Hi ), weights Hi given by

i= 1…(N −1)} is the discrete relaxation spectrum, with the relaxation

N−1

Hi = % YjiYki,

(40)

j,k = 1

and the relaxation times ti by ti =

lH . bi

(41)

Here, bi are the eigenvalues of the matrix bmn, and the orthogonal matrix Yjk diagonalizes the matrix bmn. In the Zimm model, the relaxation weights Hj are all equal to unity. In the twofold normal approximation, they are not equal to unity, though their sum is equal to N−1. O8 ttinger and Zylka [15] have suggested that in the presence of nonlinear effects such as hydrodynamic interaction, it is impossible to interpret the relaxation spectrum in terms of normal modes, since the number of relaxation times in a bead-spring chain is much larger than the number of springs. In the context of the twofold normal approximation, in spite of the presence of fluctuating hydrodynamic interaction, there are only (N −1) relaxation times. However, the polymer relaxation process cannot be thought of as a superposition of (N−1) ‘independent’ elementary relaxation processes since the ‘normal’ variables Q %j are not independent of each other.

3.2. Small amplitude oscillatory shear flow The memory integral expansion (35) can be used to derive the exact material functions for small amplitude oscillatory shear flow. This flow is characterized by the following matrix representation for the tensor k(t) in the laboratory-fixed coordinate system:

Á0 k(t) =g; 0 cos vt Ã0 Ä0

1 0 0

0Â 0Ã. 0Å

(42)

Here, g; 0 is the ‘small’ amplitude, and v is the frequency. The material function that represents such flows is the complex viscosity, h*= h% −ih¦, since the stresses depend on time through the relation [14]: t pyx = − h%(v) g; 0 cos vt− h¦(v)g; 0 sin vt.

(43)

Expressions for the real and imaginary parts of the complex viscosity can be found from the − ivs ds [14]. Using Eq. (39) for the memory function G(s) leads to expression h*= 0 G(s) e N−1

h%(v) =nkBT %

j=1 N−1

h¦(v) =nkBT %

j=1

tj Hj , 1+ (tjv)2

(44)

vt 2j Hj . 1+ (tjv)2

(45)

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

257

In the Gaussian approximation, the calculation of the complex viscosity involves the inversion of a (N − 1)2 × (N− 1)2 matrix for each frequency v [5]. Here, a single calculation of the eigenvalues and eigenvectors of the (N− 1)×(N −1) matrix bmn (see Eq. (32)) is adequate to obtain h* at all frequencies. Fig. 1 compares the storage modulus G% =vh¦/(nkBT) and the loss modulus G¦=vh%/(nkBT) predicted by the twofold normal approximation with those predicted by the Gaussian approximation, for N = 30 and h*= 0.15. Within the resolution of the figure, the two approximations give identical results, regardless of whether the Rouse or Zimm orthogonal matrix is used for the purpose of diagonalization. The same is found to be true for the loss tangent tan d, defined as tan d= G¦/G%. It is very difficult to use the Gaussian approximation for chains with N \30 beads since the necessity to store and invert a (N− 1)2 × (N−1)2 matrix makes the computer memory required scale with chain length as N 4 and the CPU time required scale as N 6 (for each value of v). In the linear viscoelastic limit of the twofold normal approximation, as mentioned above, only the eigenvalues and eigenvectors of the (N −1)×(N −1) matrix bmn are necessary. Consequently, the memory and CPU requirement are less severe in this case. However, there is a difference in the memory required depending on whether the Rouse or Zimm orthogonal matrix is used for the purpose of diagonalization, for the following reason. When the Rouse orthogonal matrix is used, sums of the following kind can be performed analytically [16]: %mnj =

max(m,n) − 1

%

Plj.

(46)

l = min(m,n)

As a result, in order to calculate ajk for every index j and k (see Eq. (34)), one only has to sum numerically over the indices m and n. This leads to the CPU time required scaling with chain

Fig. 1. Dimensionless storage modulus G%=vh¦/(nkBT) and loss modulus G¦= vh%/(nkBT) vs. dimensionless frequency lHv for chains with N=30 and h*=0.15. The continuous lines represent the Gaussian approximation, the symbol () represents the twofold normal Rouse approximation, and the symbol (x) represents the twofold normal Zimm approximation.

258

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

length as N 4 while the memory required scales as N 2. When the Zimm orthogonal matrix is used, the quantities mnj must be evaluated numerically, and in order to maintain N 4 scaling for the CPU time, it is necessary to store their values. This makes the memory required scale as N 3. In a sense, the assumptions of the Zimm model are adequate for small amplitude oscillatory shear flows since there is remarkable agreement between its predictions and experimental data [17]. However, there is much to be gained by deriving expressions for the quantities h% and h¦ with the more sophisticated approximations. This is because in the limit of v“0, they can be used to obtain exact expressions for the zero shear rate viscosity and first normal stress difference in steady shear flow. It is also in this limit that one can first fruitfully exploit the reduced CPU time of the twofold normal approximation. These statements are elaborated in the section below.

3.3. Zero shear rate limit of steady shear flow Steady shear flows are described by a k tensor which has the following matrix representation in the laboratory-fixed coordinate system:

Á0 1 0Â k =g; Ã0 0 0Ã. Ä0 0 0Å

(47)

Here g; is the constant shear rate. The three independent material functions used to characterize such flows are the viscosity, hp, and the normal stress differences, C1 and C2. These functions are defined by the following relations [14]: t pxy = − g; hp,

(48a)

t pxx − t pyy = − g; 2C1,

(48b)

t pyy − t pzz = − g; 2C2.

(48c)

As mentioned above, the zero shear rate viscosity and first normal stress difference can be obtained from Eq. (44) and Eq. (45) in the limit of vanishing frequency: N−1

h%(v)= nkBT % tj Hj, hp,0 = vlim “0

(49)

j=1

C1,0 = vlim “0

N−1 2h¦(v) = nkBT % 2t 2j Hj. v j=1

(50)

The calculation of these zero shear rate properties requires the same matrix operations as the calculation of the small amplitude oscillatory shear flow material functions. As a result, the scaling with chain length N, of CPU time and memory required, are identical in both the flow situations. Chains with N5 400 beads were examined with the twofold normal Rouse approximation, while the increased memory requirement of the twofold normal Zimm approximation enabled only examination of chains with N5150 beads.

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

259

Fig. 2. The ratio, R0(h*, N)=hp,0/h Zp,0, of the zero shear rate viscosity predicted by the Gaussian approximation and the twofold normal approximation to that predicted by the Zimm model, vs. N − 1/2, for h*= 0.15 and h*=0.25. Continuous lines represent the Gaussian approximation, while the chain-dashed and dotted lines represent the twofold normal Zimm and twofold Rouse approximations, respectively.

The ratios R0 and R1, defined by R0 = hp,0/h Zp,0 and R1 =C1,0/CZ1,0 (where h Zp,0 and CZ1,0 are the zero shear rate properties calculated with the Zimm model), are useful quantities for the comparison of the Gaussian approximation with the twofold normal approximation. Figs. 2 and 3 display the prediction of R0 and R1, respectively made by the two approximations for two

Fig. 3. The ratio, R1(h*, N)=C1,0/CZ1,0, of the zero shear rate first normal stress difference predicted by the Gaussian approximation and the twofold normal approxiamtion to that predicted by the Zimm model, vs. N − 1/2, for h*=0.15 and h*= 0.25. Continuous lines represent the approximation, while the chain-dashed and dotted lines represent the twofold normal Zimm and twofold Rouse approximations, respectively.

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

260

Table 1 Extrapolation of the ratios R0(h*, N), and R1(h*, N) to the limit of infinite chain length R0 1

R0 0

GA TNR TNZ

h*= 0.15

h*= 0.25

h*= 0.15

h*=0.25

0.7289 0.005 0.7149 0.002 0.7319 0.005

0.72990.002 0.71229 0.0002 0.7271 90.0006

0.72190.004 0.676690.0001 0.71769 0.0003

0.72190.003 0.67669 0.0001 0.71759 0.0001

The acronyms GA, TNR and TNZ stand for Gaussian approximation, twofold normal Rouse and twofold normal Zimm, respectively. The values for the Gaussian approximation are taken from O8 ttinger [5].

values of h* and several values of N. Clearly, using the Zimm orthogonal matrix leads to better agreement with the Gaussian approximation than using the Rouse orthogonal matrix. The agreement in both cases seems to get poorer with larger N and h*. The difference in the accuracy of prediction between the twofold normal Rouse and twofold normal Zimm approximations (which is larger in the case of R1 than for R0), may be due to the fact that the Zimm modes capture at least some aspects of hydrodynamic interaction. In the generalized Zimm model, the value of the ratios R0 and R1 is unity [3]. In the case of the Gaussian approximation, O8 ttinger [5] showed by plotting their values vs. N − 1/2, that these ratios become independent of h* when the data is extrapolated to the limit N “ . The reason for using N − 1/2 as the independent variable was motivated by the fact that for the Zimm model, the leading corrections to the N“ limit are of order N − 1/2. In Figs. 2 and 3, we use N − 1/2 as the abscissa with the expectation that the leading order corrections are of this order even in the case of the twofold normal approximation. The results of extrapolating the data in Figs. 2 and 3 to the limit N “ are displayed in Table 1, where, R0 0 =limN “ R0(N, h*) and R0 1 = limN “ R1(N, h*). A rational function extrapolation algorithm was used to carry out the extrapolation of the data [18]. The values for the Gaussian approximation are taken from [5]. Table 1 confirms that the ratios are independent of h* in the present instance also, and are thus universal ratios. The difference between the values of R0 1 for the two values of h* is less than the difference between the values of R0 0 for the same values of h*. Except for the value of R0 0 at h* =0.15, the use of the twofold normal approximation leads to one significant figure of accuracy more than the Gaussian approximation. This is obviously due to the fact that one can examine chains with larger N with the twofold normal approximation. In the generalized Zimm model, there is a special value of h* =0.2424… for which the leading corrections are of order N − 1. Here too one notices that the error bars for h* = 0.25 are consistently lower than for h* =0.15. From the values of R0 0 and R0 1 in Table 1 it is clear that the maximum difference between the twofold normal (Rouse) approximation and the Gaussian approximation is about 6%. Table 2 presents four more ratios of zero shear rate properties, defined as Uhl =

hp,0 ; nkBTl1

UhR = nlim “0

hp,0 ; nhs(4pR 3g/3)

UCh =

nkBTC1,0 ; h 2p,0

2,0

UCC =C C1,0

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

261

where, l1 is the longest relaxation time, and Rg is the root-mean-square radius of gyration at equilibrium. The exact Zimm values of all the ratios, and the Gaussian approximation values of UhR and UCh have been reported by O8 ttinger [19]. The values of C2,0 for various values of N (which cannot be obtained from a first order codeformational memory integral expansion) required to evaluate UCC were obtained by extrapolation of finite shear rate data to zero shear rate. The Gaussian approximation values of Uhl and UCC, and all the twofold normal Zimm values were obtained (with the same extrapolation algorithm used above) for two values of h*. In each case, as expected, these two values were close to each other. The uncertainty in the value of the last digit has been reported in parenthesis. As in all the comparisons made so far, the twofold normal values are close to those predicted by the Gaussian approximation. Miyaki et al. [20] have experimentally obtained a value of UhR =1.49 Eq. (6) for polystyrene in cyclohexane at the theta temperature.

4. Viscometric functions in steady shear flow

4.1. Solution procedure The viscometric functions that characterize steady shear flow have already been introduced in Eq. (48a). It is the purpose of this section to discuss the numerical details involved in the calculation of h, C1 and C2 for various finite values of the shear rate g; . In the consistent averaging approximation [21], and in the Gaussian approximation [8], steady state viscometric functions are obtained by numerically integrating the system of differential equations (9) for the independent components of the tensor sjk, starting from equilibrium, with respect to time, until steady state is reached. (As mentioned earlier, the equations of the former approximation are obtained by dropping the last term in Eq. (9)). The CPU time required scales with chain length as N 4.5 in both these cases. The origin of the exponent 4.5 can be understood by means of the following argument. The right hand side of Eq. (9) can be found, for every index j and k, by summing over the index m from 1 to N −1. The time required for this procedure therefore scales as N 3. (In the case of the Gaussian approximation, the sums over the indices l and p must be performed, and the results stored, prior to the summation over the index m). The additional factor of N 1.5 is proportional to the time required to reach the steady state Table 2 Universal ratios of prefactors in the limit of zero shear rate

Zimm GA TNZ

Uhl

UhR

UCh

UCC

2.39 1.835(1) 1.835(1)

1.66425 1.213(3) 1.210(2)

0.413865 0.560(3) 0.5615(3)

0.0 −0.0226(5) −0.0232(1)

The exact Zimm values and the Gaussian approximation values for UhR and UCh are reproduced from O8ttinger [19]. Numbers in parentheses indicate the uncertainity in the last figure.

262

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

solution. This is related to the fact that the longest relaxation time scales with chain length as N 1.5 lH [19], where lH is the local relaxation time. Any scheme that purports to obtain the viscometric functions with a reduced computational time must consequently scale with chain length with an exponent of N less than 4.5. It is not difficult to see, arguing as above, that evaluating the right hand side of the governing Eq. (25) for the decoupling approximation and the twofold normal approximation will also require a CPU time that scales as N 3. In the case of the decoupling approximation, the right hand side (with the last term dropped), can be found for every index j, by summing over two indices (l and p in Eq. (26)) from 1 to N−1. In the case of the twofold normal approximation, the sums over the indices k, l, m, n and p (see Eq. (27)) required to find the last term on the right hand side can be reduced to a sum over only the indices m and n by evaluating and suitably storing the sums over the remaining indices. Once the right hand side of Eq. (25) is evaluated, subsequent integration in time until steady state is reached, results in a scheme whose CPU time requirement scales as N 4.5, just as in the case of the consistent averaging and Gaussian approximations. It is very clear therefore that the reduction in the number of equations to be solved from (N −1)2 to (N −1) is not the significant factor responsible for a reduction in the computational time required, as claimed by Kishbaugh and McHugh. Indeed, the possibility of reducing the computational time seems entirely related to reducing the factor N 1.5 required for an explicit time integration. The problem with the time integration is that the time step should be small compared to the local time scale lH. Since we are only interested in the steady state solution and not in the fine time resolution, this problem can be solved, at least in principle, by adopting an implicit time integration scheme, whereby a time step much larger than lH can be used; one that scales with N x, with 0 5 x51.5. This attractive proposal is however not implementable for the Gaussian approximation since typical implicit schemes require the evaluation and storage of the Jacobian, a process whose memory requirement scales as N 4. In the case of the twofold normal approximation, this scheme does lead to a reduction in computational time compared to explicit time integration. However, the reduction is not as significant as an alternative scheme discussed below. In the decoupling approximation, Kishbaugh and McHugh [10] avoid the additional factor of N 1.5 required to integrate to steady state by solving the system of non-linear algebraic equations in Eq. (25) to zero. The equations, which may be that result on setting the derivatives ds %/dt j written in the following form, are solved using the method of successive substitution, and the steady state rheological properties are obtained directly: s % j n + 1 = Fj ({s %1, s %2, s %3,…, s %N − 1} n )

(51)

Here, n indicates the solution found after n iterations. Kishbaugh and McHugh indicate that the CPU time for their scheme scales as N 3. The significant reduction in computational time as compared with the consistent averaging approximation enables them to examine chains with as many as 200 beads. It is worth noting, as mentioned earlier in Section 2.4, that proper exploitation of the diagonalization assumption would actually have resulted in the CPU time scaling N 2. By rewriting Eq. (26) as

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

263

N−1

Lj Pmj = % PljA( lm, l=1

one could evaluate Lj in N 2 steps by picking a fixed value of the index m. Since the decoupling approximation is formulated incorrectly, the result would then depend on the choice of m and consequently lead to inconsistent results. Following Kishbaugh and McHugh [10], we have adopted a successive substitution scheme =0 in Eq. (25), and for the twofold normal approximation by setting the derivatives ds %/dt j rewriting it in the form of Eq. (51). As a consequence, it becomes possible to obtain the viscometric functions for shear flow predicted by it. For low shear rates, the equilibrium solution is used as the initial guess for the solution. For low shear rates, the equilibrium solution is used as the initial guess for the solution. The solution is considered to have converged when successive values of the normalized material functions differ by less than a specified tolerance, typically 10 − 8. The CPU time required to obtain the steady state solution with this scheme scales as N x, with 3.155x 5 3.3. This significant reduction enables us to examine chains with as many as 100 beads. Summarizing the discussion above, one might say that diagonalizing the covariance matrix sjk reduces the number of equations to be solved by roughly a factor of N. While this does not reduce the number of function evaluations, it becomes relatively easy to develop a successive substitution algorithm to obtain the steady state solution of these equations. We have unfortunately not been able to device an algorithm to directly obtain the steady state solution of the Eq. (9) that govern the Gaussian approximation. Even if this were possible, the fact that there are a factor of N more equations than in the twofold case may only lead to a marginal improvement in computational time as compared to the option of integrating the differential equations in time. Before concluding this section, it is worthwhile to compare the CPU time requirement of the twofold normal approximation in the linear viscoelastic case and in the general case. In the former case, we have seen in Section 3 that the CPU time required to find the shear flow material functions at any time scaled as N 4 (see Eq. (35)). This scaling arose from a need to evaluate the elements of the matrix bjk. At first sight it appears paradoxical that this requirement is greater than in the general case, where we have just seen that CPU time scales as N 3.3. However, it is important to note that the N 3.3 scaling holds when the steady state rheological properties are obtained directly, rather than through an explicit time integration—in which case the CPU time would scale as N 4.5. Even in the linear viscoelastic limit, the evaluation of the right hand side of Eq. (31) would require a CPU time that scales as N 3 (in this case, it becomes unnecessary to explicitly evaluate bjk for every index j and k). The steady state values of ej could then be obtained with a successive substitution procedure which one anticipates would require a CPU time proportional to N 3.3. Presumably, the material functions at zero shear rate could be obtained by extrapolating the results at finite shear rates to the limit of zero shear. Note however that the saving in CPU time is made at the cost of losing information about the evolution of the material functions in time. In the next section, we present predictions of the viscometric functions for simple shear flows. O8 ttinger [5] has pointed out that the symmetry implied by the relation in Eq. (8), and the symmetry between the two ends of the bead-spring chain, contribute to a reduction in the number of independent components of the tensor sjk by a factor of 4 for large chains. In

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

264

addition, since shear flows remain invariant when the direction of the z-axis is reversed [14], the xz, yz, zx and zy components of these tensors must be zero. In the twofold normal approximation, the symmetry of the flow field, and Eq. (8) together imply that the variances s %j must have the following representation in the laboratory-fixed coordinate system:

Áa 1j a 0j 0 Â s = Ãa 0j a 2j 0 Ã. Ä 0 0 a 3j Å % j

(52)

Remaining on the note of symmetries, the following point is worth noting. As mentioned before, the hydrodynamic interaction functions H(s) and K(s) can be evaluated in terms of elliptic integrals. Very explicit detail on the method of evaluation can be found in the Appendix of [8]. Here we point out in addition that the symmetry imposed by the flow field can be used to show that only 14 out of the 81 components of the fourth rank tensor K in the laboratory-fixed coordinate system are required in simple shear flow.

4.2. Numerical results for chains of finite length Viscometric functions predicted by the twofold normal approximation for Rouse chains with hydrodynamic interaction are now presented. In all the figures below, the viscometric functions h, C1 and C2 are plotted against the reduced shear rate, b, a dimensionless quantity which is defined as b=

[h]0 hsMg; , NAkBT

(53)

where M is the molecular weight of the polymer, NA is Avagadros number, and [h]0 is the intrinsic viscosity in the limit of zero shear rate. It is easy to show, using the definition of the intrinsic viscosity, that for dilute solutions b =hp,0g; /nkBT. For both the Gaussian approximation and the twofold normal approximation, hp,0 can be evaluated explicitly. Fig. 4 compares the prediction of hp/lHnkBT made by the Gaussian approximation with the prediction of the twofold normal approximation, for parameter values N =45 and h* =0.15 at low shear rates, and for N=15 and h*= 0.15 at high shear rates (see inset to Fig. 4). Results at high shear rates were not obtained for N =45 with the Gaussian approximation due to the excessive requirement of computational time. Purely for the purpose of comparison, the value of hp,0 given by the Gaussian approximation is used to calculate the reduced shear rate in all the cases. Using the Zimm orthogonal matrix leads to results indistinguishable (within the resolution of the figure) from the Gaussian approximation at low to moderate shear rates. However, as can be seen from the inset, at high shear rates the twofold normal Zimm approximation predicts a slightly lower value of hp. Use of the Rouse orthogonal matrix causes a small error at low to moderate shear rates, which becomes negligible at high shear rates. This discrepancy between the predictions of the twofold normal Rouse and twofold normal Zimm approximations was also present in the predictions of the zero shear rate viscosity and first normal stress difference (see Figs. 2 and 3). As remarked earlier, the Zimm orthogonal matrix captures at some aspects of hydrodynamic interaction, and therefore leads to more accurate results at low shear rates

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

265

Fig. 4. Dimensionless polymer contribution to the viscosity vs. the reduced shear rate b of the Gaussian approximation, for N =45 at low shear rates and N= 15 at high shear rates (inset). In both cases h*=0.15. Continuous lines represent the Gaussian approximation, while the chain-dashed and dotted lines represent the twofold normal Zimm and Rouse approximations, respectively.

where hydrodynamic interaction is important. At high shear rates, the beads of the chain are far apart, and hydrodynamic interaction is nearly switched off. As a result, the viscometric functions approach their Rouse values. Such a situation does mean of course that the model has become poor. In order to properly incorporate sufficiently local details at high shear rates, one would have to consider chains with a larger number of beads. The good performance of Rouse modes therefore indicates the breakdown of the model (at that number of beads). The curves comparing the prediction of the Gaussian and twofold normal approximations for the first normal stress difference are very similar in shape to those shown in Fig. 4 for the viscosity; the Zimm modes being more accurate at low shear rates and the Rouse modes being more accurate at high shear rates where, as discussed above, the effects of hydrodynamic interaction are less pronounced and the limit of the model’s validity is approached. The second normal stress difference predictions are compared in Fig. 5 for parameter values N= 45 and h* =0.15. As mentioned before, the Gaussian approximation predicts a negative value at low shear rates, in agreement with simulations. As the shear rate is increased, the function goes through a small positive maximum before going to zero at still higher shear rates. The twofold normal approximation also displays all these features, and differs marginally from the Gaussian approximation at low shear rates. As with the other two material functions, the Zimm modes are more accurate at low shear rates. The value N = 15 is the largest that Zylka [8] was able to implement while solving the equations of the Gaussian approximation with the help of a super computer. We were limited to N = 45 using an IBM RISC 6000/560 work station. It took 679 h of CPU time to generate the Gaussian approximation data in Fig. 4. On the other hand the normal approximation data was produced in roughly 17 min. Because of this great reduction in computational

266

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

time, we were able to examine chains with N5100 beads using the twofold normal approximation. The results for these chains are discussed below, and the extrapolation of all the accumulated finite chain results to the infinite chain limit is considered in the next section. Fig. 6(a,b) display the relative values of the viscosity and first normal stress difference for N= 100 and two values of the hydrodynamic interaction parameter, h*=0.15 and h* =0.25 obtained with the twofold normal Zimm approximation. In both the figures, the value of hp,0 given by the twofold normal Zimm approximation is used to calculate the reduced shear rate b. The shape of these functions at low to moderate shear rates predicted by bead-spring chain models with hydrodynamic interaction is now well established. The curve for the viscosity, Fig. 6(a), reveals an initial shear thinning as the shear rate is increased, reaching a minimum at approximately ln b=1, and then shear thickening at higher shear rates. Notice that while for low shear rates the curves for the two values of h* are coincident, they begin to diverge for reduced shear rates greater than 10. Thus, as was remarked earlier, the results of the model are dependent on the values chosen for the parameters N and h*. However, as we shall see below, it is possible to transcend this limitation and predict universal viscometric functions.

4.3. Uni6ersal 6iscometric functions Bead-spring chain models account for the capacity of macromolecules to stretch and orient themselves under the influence of a flow field. Kinetic theories that use such models to predict the rheological behavior of polymer solutions are based on the expectation that these features of the macromolecule are the primary cause for the solution’s viscoelasticity. The successful prediction of various qualitative aspects of the solution’s behavior justifies this expectation. However, as seen above, the predicted material functions are not universal since they depend on parameters related to the mechanical model. In the limit of N “ , however, O8 ttinger [3] has

Fig. 5. Dimensionless polymer contribution to the second normal stress difference vs. the reduced shear rate b of the Gaussian approximation, for N=45 and h*= 0.15. The continuous line represents Gaussian approximation, while the chain-dashed and dotted lines represent the twofold normal Zimm and Rouse approximations, respectively.

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

267

Fig. 6. Normalized polymer contribution to (a) the viscosity and (b) the first normal stress difference, vs. the reduced shear rate b predicted by the twofold normal Zimm approximation for chains with N= 100 beads for two values of the hydrodynamic interaction parameter h*.

verified analytically that for the generalized Zimm model, the predictions for various values of h* coincide and are consequently parameter free. (The limit N “ is carried out by letting H“ 0 while simultaneously holding the end-to-end distance fixed.) It is the purpose of this section to show numerically that in the limit of N“ , the results of the twofold normal approximation are also parameter free and consequently universal. In the generalized Zimm model, it was shown that leading corrections to the infinite chain limit are of order N − 0.5. We expect that this is also the case in the twofold normal approximation. Consequently, the universal dependence of relative viscosity on the reduced shear rate was obtained by extrapolating values of hp/hp,0 as a function of N − 0.5 (accumulated for chains with N5 100) to the N“ limit at values of the reduced shear rate b. The scheme used for

268

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

this purpose is identical to the one used in Section 3.3. The Zimm orthogonal matrix was used to diagonalize sjk as it leads to more accurate results at low to moderate shear rates. Fig. 7 displays the results of the extrapolation, and clearly these are universal since the values of relative viscosity for the two values of relative viscosity for the two values of h* lie within 1% of each other up until a reduced shear rate of roughly 25. Beyond this shear rate, they begin to diverge indicating that the present data accumulated for chains with N5 100 is insufficient to carry out an accurate extrapolation at shear rates greater than 25. At these high values of shear rate, the 100 bead chain is stretched to such an extent that one begins to probe length scales, lets say of order a, that are model dependent. However, one expects that for chains much longer than N = 100, there will exist enough length scales smaller than a such that at a shear rate b\ 25, model dependent features are absent at the length scale a. In the same figure, the result of a renormalization group calculation of the relative viscosity [13] is also presented. The renormalization group method is essentially a means for refining the results of a first order perturbation theory in the hydrodynamic interaction parameter by adding an infinite number of terms of higher order using an exponentiation procedure [12]. The Gaussian and twofold normal approximations are non-pertubative theories that contain an infinite number of higher orders, but unlike the renormalization group method they are ‘uncontrolled’ approximations. Note though that both the Gaussian approximation and the renormalization group method are exact to first order in the hydrodynamic interaction parameter h*. It is clear from the figure that the two methods lead to significantly different results at moderate to high shear rates. The renormalization group does not predict a minimum, while the twofold normal approximation predicts one, albeit a small departure from the zero shear rate value. The steep upturn in the relative viscosity begins to occur in the twofold normal approximation at a shear rate that is higher than that at which the more gradual increase predicted by the renormalization group begins to occur.

Fig. 7. Universal predictions of normalized polymer contribution to the viscosity vs. the reduced shear rate b. Results obtained by extrapolating the twofold normal Zimm approximation data for h*= 0.15 are represented by the continuous line, while the symbol () represents results for h* =0.25. The dotted line is the result of a tion group (RG) calculation obtained by Zylka and O8 ttinger [13].

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

269

Fig. 8. Universal predictions of normalized polymer contribution to the first normal stress difference vs. the reduced shear rate b. Parameter values and symbols are as in the caption to Fig. 7.

The universal dependence of the relative first normal stress difference on the reduced shear rate was obtained by extrapolating values of C1/C1,0 as a function of N − 0.5 (accumulated for chains with N 5 100) to the N“ limit at various values of the reduced shear rate b. The scheme used for this purpose is identical to the one used in Section 3.3. Fig. 8 reveals that extrapolation is unreliable for b \ 20. Once again, renormalization group results are significantly different from the twofold normal approximation predictions. The result of extrapolating the ratio C2/C1 to the N“ limit is displayed in Fig. 9. This becomes necessary because the zero shear rate expression for the second normal stress difference is not known. The curves for the two values of h* are nearly coincident until the highest value of shear rate examined here. Note that the value of this ratio predicted by the renormalization group method is higher than that predicted by the twofold normal approximation at all values of the shear rate. 5. Conclusions A new approximate method is introduced for treating the effects of hydrodynamic interaction in a bead-spring chain model for dilute polymer solutions that makes it possible to predict the viscometric functions in simple shear flow in a computationally efficient manner. Like the most recent such approximation, namely, the Gaussian approximation, it assumes that the configurational distribution function is Gaussian. In addition, it assumes that the Rouse or Zimm orthogonal matrix diagonalizes the matrix of bead connector vector covariances. This leads to a significant reduction in computational time without significantly affecting the accuracy of the Gaussian approximation. As a result, long chain results may be obtained. Explicit expressions for the frequency dependence of the loss and storage moduli G% and G¦ in small amplitude oscillatory shear flow can be obtained with the twofold normal approxima-

270

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

tion. The values of these modulii found at several frequencies are indistinguishable from those obtained with the Gaussian approximation for chains with N 530 beads, regardless of whether the Rouse or the Zimm orthogonal matrix is used to perform the diagonalization (Fig. 1). Expressions for the zero shear rate viscosity and first normal stress difference in steady shear flow are obtained as the zero frequency limit of the material functions in oscillatory shear flow (Eq. (49) and Eq. (50)). The ratios R0 = hp,0/h Zp,0 and R1 =C1,0/CZ1,0, are displayed in Figs. 2 and 3. Values of these ratios calculated with the twofold normal approximation using the Zimm orthogonal matrix are closer to the Gaussian approximation results than those that are obtained when the Rouse orthogonal matrix is used. This follows from the fact that the Zimm orthogonal matrix reflects to some extent the presence of hydrodynamic interaction. Universal (independent of model parameters) values of the ratios R0 and R1 are found by extrapolating finite chain length results to the limit N “ . Use of the twofold normal approximation leads to one significant figure of accuracy more than the Gaussian approximation for most cases (Table 1). The values of the universal ratios Uhl, UhR, UCh, and UCC predicted by the twofold normal Zimm approximation are close to the predictions of the Gaussian approximation. These are compared to the exact predictions of the Zimm model in Table 2. The requirements in terms of computer memory and processor time for the twofold normal approximation in the linear viscoelastic limit as compared to that for the Gaussian approximation, are shown in Table 3. Reasons for the various scaling exponents are explored in Section 3.3. The prediction of the shear rate dependence of the viscosity, the first normal stress difference and the second normal stress difference in steady shear flow by the twofold normal Zimm approximation agrees excellently with the Gaussian approximation for low to moderate shear rates, but departs from it slightly at high shear rates. When the Rouse orthogonal matrix is used

Fig. 9. Universal predictions of normalized polymer contribution to the second normal stress difference vs. the reduced shear rate b. Parameter values and symbols are as in the caption to Fig. 7.

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

271

Table 3 Scaling of memory and CPU time requirement with chain length N for the Gaussian and twofold normal approximations in linear viscoelastic and steady shear flows Linear viscoelastic

GA TNZ TNR

Steady shear flow

Memory

Time

Memory

Time

N4 N3 N2

N6 N4 N4

N2 N2 N2

N 4.5 N 3.3 N 3.3

GA, TNR and TNZ stand for the Gaussian approximation, the twofold normal Rouse and the twofold normal Zimm, respectively.

there is some error at low shear rates which decreases with increasing shear rate (Figs. 4 and 5). However, the convergence to Rouse values at high shear rates indicates that the model predictions are inadequate. The requirements in terms of computer memory and processor time for the numerical evaluation of the viscometric functions using the twofold normal approximation as compared with that for the Gaussian approximation are displayed in Table 3. A detailed examination of the origin of the different exponents of the chain length N is carried out in Section 4.1. The relative viscosity and the relative first normal stress difference as a function of the reduced shear rate for chains with a hundred beads predicted by the twofold normal Zimm approximation are presented in Fig. 6 for two values of the hydrodynamic interaction parameter h*. The divergence of the curves for reduced shear rates greater than ten indicates the model dependent nature of the predictions. The most important consequence of the introduction of the twofold normal approximation is that universal predictions for the dependence of hp/hp,0, C1/C1,0 and C2/C1 on the reduced shear rate b can be obtained by extrapolating finite chain length results to the limit N “ . Comparison with similar predictions made by using a renormalization group method reveals that the two procedures agree qualitatively but differ significantly quantitatively (Figs. 7–9).

Acknowledgements We gratefully acknowledge financial support provided by the Swiss National Foundation for Scientific Research under Grant No. 20-36073.92.

References [1] [2] [3] [4]

P.E. Rouse, J. Chem. Phys. 21 (1953) 1272–1280. B.H. Zimm, J. Chem. Phys. 24 (1956) 269–281. H.C. O8 ttinger, J. Chem. Phys. 86 (1987) 3731– 3749. H.C. O8 ttinger, J. Non-Newtonian Fluid Mech. 26 (1987) 207 – 246. H.C. O8 ttinger, J. Chem. Phys. 90 (1989) 463–473.

J.R. Prakash, H.C. O8 ttinger / J. Non-Newtonian Fluid Mech. 71 (1997) 245–272

272 [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

W. Zylka, H.C. O8 ttinger, J. Chem. Phys. 90 (1989) 474 – 480. L.E. Wedgewood, J. Non-Newtonian Fluid Mech. 31 (1989) 127 – 142. W. Zylka, J. Chem. Phys. 94 (1991) 4628–4636. J.J. Magda, R.G. Larson, M.E. Mackay, J. Chem. Phys. 89 (1988) 2504 – 2513. A.J. Kishbaugh, A.J. McHugh, J. Non-Newtonian Fluid Mech. 34 (1990) 181 – 206. R.B. Bird, C.F. Curtiss, R.C. Armstrong and O. Hassager, Dynamics of Polymeric Liquids, Vol. 2, Kinetic Theory, 2nd ed., Wiley, 1987. H.C. O8 ttinger, Y. Rabin, J. Non-Newtonian Fluid Mech. 33 (1989) 53 – 93. W. Zylka, H.C. O8 ttinger, Macromolecules 24 (1991) 484 – 494. R.B. Bird, R.C. Armstrong and O. Hassager, Dynamics of Polymeric Liquids, Vol. 1, Fluid Mechanics, 2nd ed., Wiley, 1987. H.C. O8 ttinger, W. Zylka, J. Rheol. 36 (1992) 885 – 910. I.S. Gradshteyn and I.M. Ryzhik, Tables of Integrals, Series and Products, Academic, 1980. K. Osaki, J.L. Schrag, J.D. Ferry, Macromolecules 5 (1972) 144 – 147. W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in FORTRAN, 2nd ed., Cambridge University Press, 1992. H.C. O8 ttinger, Stochastic Processes in Polymeric Fluids, Springer-Verlag, 1996. Y. Miyaki, Y. Einaga, H. Fujita, M. Fukuda, Macromolecules 13 (1980) 588 – 592. L.E. Wedgewood, H.C. O8 ttinger, J. Non-Newtonian Fluid Mech. 27 (1988) 245 – 264.

.