Viscosity prediction from a modified square well intermolecular potential model

=.

. .i ~.

llllglg

•

ELSEVIER

Fluid Phase Equilibria 117 (1996) 378-385

Viscosity prediction from a modified square well intermolecular potential model Wayne D. Monnery, Anil K. Mehrotra and William Y. Svrcek Department of Chemical and Petroleum Engineering, The University of Calgary, Calgary, Alberta T2N IN4, Canada

ABSTRACT A statistical mechanical based model for simultaneously calculating liquid and gas phase viscosities using the square well intermolecular potential has been developed. Theoretical results are utilized along with modifications to correct the model for the molecular chaos assumption and the inadequacy of the square well potential. The radial distribution function values are determined from perturbation theory. Model parameters are obtained from liquid and gas viscosity data and generalized with group contributions. With the resulting model, liquid and gas viscosities of a variety of nonpolar compounds are correlated within 5% and 3%, respectively, and predicted within 8% and 5%, respectively.

1. INTRODUCTION In the chemical and process industries, viscosity is an important property in hydraulics calculations for surface processing facilities, pipeline systems and flow through porous media. With the increased popularity of process and reservoir simulators, there is a need for a reliable and analytical predictive model for viscosity calculations. Also, like some PVT data, viscosity data provide a physical link to the intermolecular pair potential. Recently, an extensive review of viscosity prediction and correlation methods was completed by Monnery et al. [1], which showed that although there is an abundance of literature and many viscosity models are available, amongst the predictive methods, most have only been fit and/or are restricted to a narrow range of compounds or mixtures with only a few empirical models widely generalized. Also, typically, many are applicable over limited ranges of temperature and pressure and relatively few are applicable continuously from dilute gas to liquid phase. Of those, simple semi-theoretical methods based in statistical mechanics (such as corresponding states, hard sphere and square well theories) appear to be the most promising for engineering use. Their theoretical basis gives the correct qualitative trends with temperature and density, resulting in more confident extrapolation. Also, the parameters typically have some physical significance. However, more work is needed to apply these models to a wide range of compounds and mixtures over wide ranges of temperature and pressure, generalizing parameters in order to keep the models predictive. The modified hard sphere model of Chung et al. [2,3] has been generalized and applied to a variety of compounds and mixtures, although for associating fluids, requires a parameter which has been determined for relatively few compounds (e.g. water, 1alkanols). The extended corresponding states method is also predictive but does not give particularly good results for fluids significantly different in structure from the reference fluid(s). In this work, we develop a statistical mechanical based model to simultaneously calculate liquid and gas phase viscosities based on the square well theory of Davis et al. [4,5]. After modifications, the resulting semitheoretical model is applied to 75 nonpolar compounds with excellent results. 0378-3812/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0378- 38 ! 2(95)02975-3

W.D. Monnery et al. / Fluid Phase Equilibria I 17 (1996) 378-385

379

2. SQUARE W E L L MODEL The square well intermolecular potential is the simplest potential model which accounts for both repulsion and attraction, albeit crudely. However, as stated in Reed and Gubbins [6], the resulting viscosity model is an excellent compromise between realism and ease of calculation. As in dense fluid hard sphere theory, the two important assumptions in deriving the model are: (1) only 2-body interactions are important and (2) molecular chaos for velocities (but not positions), i.e. no statistical correlation between succeeding velocities of particles. Also analogous to dense hard sphere theory, a modified Boltzmann equation was derived accounting for positional correlation at higher density through the radial distribution function. By also accounting for the different types of interactions for a square well particle and solving the resulting modified Boltzmann equation, Davis et al. [4] derived the following model: 11

[1 + 0.4bp(g(o,) + R3g(Ro,)~,)] 2

bpg(°,)*R2bpg(R°,) 2+

+ 2_~(bp)[g(o,). R4g(Ro,),2]

-~

where parameters

ll/1

,, = 1 - e't*r +

(e/2kT)[l + (4/V/~)e °*r f e-':x2dx]

and ~2 are defined as:

012 = e " / k r - e / 2 k T -

2 fx2(x

2 ÷ e/kT)lae:dx

(2)

(3)

0

For a given component at a given density and temperature, this model requires five parameters, o I (or b), e, R, g(ol) and g(ROl), as well as values of the parameters 1111and ~2. 2.1. Parameter reduction The integrals in Equations (2) and (3) have been approximated by the following algebraic equations which reproduce the values of ~1 and 1[/2 within 0.8% over 0~ e / k T ~ 4 and extrapolate smoothly:

I(~1) = 1.0 / [2.2519 + 0.46614(E/k/) + 1.4843 (e/kT) 2 - 0.44117 (e/kT) 3 + 0.18849(e/kT) 4]

(4)

/(1112) = [0.50419 • 0.28304(e/kT)] / [1.0.0.15360(e/kT)]

(5)

Although slightly more complicated, the above expressions are considerably more accurate than those given by Du and Guo [7] who also modified the Davis et al. [4,5] theory. The well width parameter, R, was set to 1.5 which is a typical value shown to provide a good compromise between fitting second virial coefficient data and critical temperature and volume for real fluids [8]. According to Benavides et al. [9], with this value, the phase behavior of square well fluids is similar to a Lennard-Jones 6-12 fluid. The radial distribution function values are determined by the EXP approximation of the WCA perturbation theory [10,1 1], as proposed by Vimalchand et al. [12]: g(p,T;r)

~ gm(P ;r)

exp(e~e/kT)

(6)

380

W.D. Monnery et al. / Fluid Phase Equilibria 117 (1996) 378-385

where gHS(9 ;r) = the hard sphere value of the radial distribution function, determined from the CarnahanStarling [13] equation at o 1 and a correlation fitted to the Monte Carlo simulation results of Barker and Henderson [14] at R a p The empirical coefficients tt at o z and Ro I (-0.2 and -0.5, respectively) were determined by minimizing the deviations between the predicted values from Equation (6) and the Monte Carlo simulation values of Henderson et al. [15], matching these values within 5%. The advantage of the EXP approximation is that it reduces correctly to the theoretical dilute gas limit, although in this case, it is modified somewhat by the values of a. With the well width parameter set and the radial distribution function values determined, the square well model is reduced to two parameters, b (or o 1) and e/k. 2.2. Modifications In the ideal gas limit, the square well model reduces to rl = rlo/d(e/kT) where:

d(e/kT) = g(o,) + R2g(Ro,)fle/kT)

(7)

which is the denominator of the first term of Equation (1) and: fle/kT)

= *2 ÷ (1/6)(e/kT) =

(8)

According to Chapman-Enskog theory, the dilute gas viscosity for a fluid interacting with the square well intermolecular potential is written: rl 0.sw = (5/16 ~ 1/2)(/nkT)lr2/(o2 -sw°{='=~*, = lq0/ -sw°C='2~

(9)

This implies that in the ideal gas limit, d(e/kT) should be equivalent to °°sw c~ (2,2)*, "~ ~o long as the radial distribution functions approach the correct ideal gas limit, exp(-dkT). As shown in Figure 1, this is not the case, although they are close. Furthermore, when the perturbation theory is used for determining the radial distribution function values, d(e/kT) does not match ~s(~2)*at e/kT ~ 0.8, due to the empirical coefficients.

6 5 4

I

O

£.-~(2,2)*

[] o

d(WkT) d(WkT)

'

I

'

I

..._~

actual RDF mod.

.....-" """ t

-~3 2

0

0.0

i

I

0.5

Figure 1. d(e/kT) versus (e/kT)

I

I

1.0 dkT

i

I

1.5

~

I

2.0

W.D. Monnery et al. / Fluid Phase Equilibria 117 (1996) 378-385

381

To ensure the proper ideal gas limit is reached with the radial distribution function values determined from the perturbation theory,f(e/kT) has been modified and is determined from the following expression:

f ( c / k T ) ~ -0.811 + 0.788exp(e/2.055kT)

(10)

As previously stated, the assumptions made in the derivation of the square well model are essentially the same as those made in the derivation of the dense hard sphere model of Enskog. As such, it likely suffers from the same shortcomings. For dense hard spheres, by comparing the results of the Enskog model to the results of molecular dynamics simulations, Alder et al. [16] showed that the assumption of molecular chaos breaks down at sufficiently high densities, resulting in significant underprediction of viscosity and provided approximate corrections for this. Analogously, Michels and Trappeniers [17] have compared viscosity results of the square well model to molecular dynamics simulation results which show significant underprediction of the model at high densities and low temperatures. Based on these results, a multiplicative correction term was formulated as a function of reduced density, b9, such that the model matched the simulation results at higher densities and tended to unity in the ideal gas limit, where no correction is required: C = 1.0

+

ks(bP) k'

(11)

Parameter k5 was found to be a function of e/kT and expressed as: (12)

k s = k, + k2(e/kT) ~

With the values of k I to k4 set to match the square well model viscosity results to the molecular dynamics results of Michels and Trappeniers [ 17], values of b and e/k were regressed in an attempt to simultaneously match liquid and gas phase viscosity data for real fluids. The model could not match gas viscosities well and there was only limited success in matching liquid viscosities. Although both gas and liquid phase data were used, sensitivity analysis of the regressions showed that liquid viscosity tended to control the optimum values of b and that c/k did not have enough "correlating power" to match gas viscosity. This is the well known problem of parameters regressed from a gas phase property giving poor results in the liquid phase and vice versa, which indicates inadequacy of the intermolecular potential and resulting model [ 18,19]. Also, the square well model does not represent the temperature dependence of dilute gas viscosity as well as a Lennard-Jones 6-12 potential. Based on these facts, the dilute gas part of the correction was modified by multiplying by a power function which approximated the ratio of the reduced collision integrals of the square well and Lennard-Jones potentials such that the dilute gas square well model behaved like a Lennard-Jones model. This ultimately resulted in the correction being expressed as: c

= k o ( k r / c ) °2s + [ t ,

+ k~(c/kr)'~](bp)

~'

(13)

Physically, this correction accounts for the breakdown of the model assumptions at high densities and low temperatures and the inadequacy of the square well potential. In addition, parameters k0 to k4 were made adjustable to observe the optimum values and trends, resulting in 7 adjustable parameters in the model. As expected, the model now matched both gas and liquid viscosities very well (_<1% error). Analysis of regressed values of the parameters indicated definite trends for k I , k2 and k4, and in the interest of minimizing adjustable parameters, these were generalized as functions of molar mass and acentric factor, reducing to their square well values as the acentric factor approaches zero: k I = 0.03072 + O.O0128Mt~ TM

(14)

W.,D. Monnery et al. / Fluid Phase Equilibria I 17 (1996) 378-385

382

k2 = 0.6270

exp(-O.OO242Mt~in)

(15) (16)

k4 = 1.0 + exp(-0.03895Mt~ la) Further analysis indicated that c/k could be generalized with critical temperature adjustable parameters, k0, k3 and b, which were then regressed again.

(c/k = 0.65Tc), leaving three

3. DATABASE The data used to regress the parameters were obtained from the DIPPR database [20]. The nonpolar compounds consist of 16 n-paraffins ranging from C 1 to C2o, 14 iso-paraffins up to C 7, 7 1-olefins from C 2 to C8, C 4 and C 5 2-olefins, 2 diolefins, 12 naphthenes from C 3 to Clo, 10 aromatics from C 6 to C12, 3 acetylenes from C 2 to C 4 and 5 simple inorganic fluids. The temperature range of the gas and liquid viscosity and density data was typically from a reduced temperature of 0.48 to 0.80 except for the inorganic fluids for which data ranged up to Tr = 0.90, which was necessary in order to maintain a span of at least 55 K. However, ultimately, bounds were set in accordance with the recommended temperature ranges for the data as stated in the database along with judgement based on consistency checks. Note that, in the database, inevitably there are some gaps in available experimental data which have been filled with predicted values. An attempt was made to minimize the use of predicted data, only using it where necessary to cover the desired range of compounds in a given family. Based on quality indicators that are stated in the database, for these nonpolar compounds, predicted gas and liquid viscosity data are typically considered to be accurate within 5% and 10%, respectively.

4. M O D E L RESULTS Analysis of the regressed optimum values of the remaining parameters suggested that they were dependent on the structure of the compounds; hence, they have been generalized as functions of group contributions, with their numerical values given in Table 1.

g

b,/100 =

g

~ngAb/+

O.06617(~n~Ab) 2

/.l

(17)

/.1

ko~ = ~ noAko~

(18)

/-1

g

k3., = 1.0 . ~

no.Ak3j

(19)

.hi

The group contributions were obtained using one family at a time, starting with n-paraffins to determine CH 3 and CH 2 groups and continuing in this manner to determine the remaining groups. As expected, first members of a homologous series (methane, ethylene and acetylene) did not trend well with higher members and to overcome this problem, each of these compounds is treated as its own group (using regressed values of the

383 Table 1 Square well viscosity model parameter group contributions Group

Ab

A k0

A k3

CH 3 (end)

0.35859

0.35138

0.46613

CH 2

0.28986

-0.00405

0.22505

CH

0.31614

-0.32709

0.50214

C

0.33288

-0.64424

1.03788

CH2--CH

0.50625

0.32162

0.65625

CH--CH

0.42781

-0.02776

0.07048

CH2=C

0.48341

0.01524

0.61742

CH-C

0.45400

0.31160

0.31682

Cyclo-C 5 Ring CH 2

0.27980

0.15240

0.527 exp(-0.20nCH 2)

Cyclo-C 5 Ring CH

0.19405

-0.21698

-0.00301

Cyclo-C6 Ring CH 2

0.30033

0.13283

0.636 exp(-0.15ncn2)

Cyclo-C6 Ring CH

0.13749

-0.27153

-1.07501

ACH

0.22117

0.11850

0.527 exp(-0.11 nEll2)

AC 1 (1 Ring)

0.06539 x (0.613+0.387nAC)

-0.24900

-0.63202/nAC

AC2 (2 Rings)

0.15925

-0.11800

-1.14377

AbcH3(branch) = A bcH3(end) X exp(-0.16riCH 2 - 0.06nCH ) Ak3,cH3(branch) = Ak3,cH3(end) X exp(-0.10(nCH2(nCH + nc) ) di-olefin correction:fdio = 0.80, if(ncH2=CH + nCH=CH + nCH2=C) > 1

parameters). The results of correlated and predicted viscosities are summarized in Table 2. Considering the accuracy of the viscosity data, both correlated and predicted results are excellent. These results compare very favorably with the methods recommended in Reid et al. [21]. For example, for the n-paraffin data used here, the method of Chung et al. [2,3] has an overall AADL and AADG of 10.0% and 3.1%, respectively. For the other families, AADLs typically ranged from 10-35% while AADGs were also about 3%. The model is now being extended to several families of polar and associating fluids. An example to illustrate the viscosity calculation steps, using the group contribution method, is given in Appendix.

NOMENCLATURE

= aromatic (Table 1) = covolume = (2/3)~N A o cm3/mol = correction factor = radial distribution function at repulsion diameter g(Ro 1) = radial distribution function at attraction diameter A b C g(ol)

13,

384

W.D. Monnery et al. / F l u i d Phase Equilibria 117 (1996) 378-385

= = = = = = =

Boltzmann constant parameters in correction term k0..k4 particle mass, g m molar mass, g/mol M Avogadro's Number NA number ofjth groups in compound i nq well width parameter R T = temperature, K = critical temperature rr = reduced temperature WCA = Weeks Chandler Andersen empirical coefficient in radial distribution function perturbation theory APj = group contribution ofjth group for parameter P where P = b, k0 or k 3 E = energy well depth rl = viscosity, mPaos hard sphere dilute gas viscosity = (5/16n l/2)(mkT)l/2/Ol 2 = 3.11630xlO-3(MT)l/2/b 2/3, mPa.s lqO = p = density, mol/cm 3 0 1 = repulsion diameter (Ro I = attraction diameter) 63 = acentric factor ~(2,2)* = reduced collision integral AADL = average absolute deviation (liquid) AADG = average absolute deviation (gas) k

Table 2 Results of square well viscosity model correlation and prediction Family

No. of Compounds

Regressed Parameters

Group Parameters

R/G

% AADL

% AADG

% AADL

% AADG

n-Paraffins

17/16

1.49

1.29

3.83

2.56

Iso-paraffin s

14/14

0.81

0.56

6.18

2.26

1-Olefins

8/7

1.10

0.54

3.72

2.70

2-Olefins

4/4

0.91

0.34

4.31

2.64

Diolefins

2/2

0.86

0.26

7.75

4.42

Naphthenes

12/10

1.60

0.53

5.74

1.30

Aromatics

10/10

1.31

0.51

3.29

1.17

Acetylenes

3/2

1.64

0.46

3.39

1.94

Inorganics

5/0

1.58

1.32

75/65

1.28

0.74

4.66

2.19

Overall

R = regressed, G = group contribution

385

W.D. Monnery et aL / Fluid Phase Equilibria 117 (1996) 378-385

REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18.

Monnery, W.D., W.Y. Svrcek and A.K. Mehrotra, Can. J. Chem. Eng., 73 (1995) 3. Chung, T.H., L.L. Lee and K.E. Starling, Ind. Eng. Chem. Fundam., 23 (1984) 8. Chung, T.H., L.L. Lee and K.E. Starling, Ind. Eng. Chem. Res., 27 (1988) 671. Davis, H.T., S.A. Rice and J.V. Sengers, J. Chem. Phys., 35 (1961) 2210. Davis, H.T. and K.D. Luks, J. Phys. Chem., 69 (1965) 869. Reed, T.M. and K.E. Gubbins, "Applied Statistical Mechanics", McGraw-Hill, 1973. Du, L.G. and T.M. Guo, Chem. Eng. J., 47 (1991) 163. Collings, A.F. and I.L. McLaughlin, J. Chem. Phys., 73 (1980) 3390. Benavides, A.L., Y. Guevara and F. Del Rio, Physica, 202A (1994) 420. Weeks, J.D., D. Chandler and H.C. Andersen, J. Chem. Phys., 54 (1971) 5237. Andersen, H.C., D. Chandler and H.C. Weeks, Adv. Chem. Phys., 34 (1976) 105. Vimalchand, P., A. Thomas, I.G. Economou and M.D. Donohue, Fluid Phase Equil., 73 (1992) 39. Camahan, N.F. and K.E. Starling, J. Chem. Phys., 51 (1969) 635. Barker, J.A. and D. Henderson, Mol. Phys., 21 (1971) 187. Henderson, D., W.G. Madden and D.D, Fitts, J. Chem. Phys., 64 (1976) 5026. Alder, B.J., D.M. Gass and T.E. Wainwright, J. Chem. Phys., 56 (1970) 3813. Michels, J.P.J. and N.J. Trappeniers, Physica, 104A (1980) 243. Mazo, R.M., "Statistical Mechanical Theories of Transport Processes", International Encyclopedia of Physical Chemistry and Chemical Physics, Topic 9, Volume 1, Pergamon Press, 1967. 19. Chapman, S. and T.G. Cowling, "The Mathematical Theory of Non-Uniform Gases", 3rd edition, Cambridge University Press, 1970. 20. Daubert, T.E., R.P. Danner, H.M. Sibul and C.C. Stebbins, DIPPR® Data Compilation of Pure Compound Properties, Version 9.0, Database 11, NIST Standard Reference Data Program, Gaithersburg, MD, 1994. 21. Reid, R.C., J.M. Prausnitz and B.E. Poling, "The Properties of Gases and Liquids", McGraw-Hill, 1987.

APPENDIX Estimate the viscosity of naphthalene at 475 K (experimental "qexp = 0.343 mPaos). Data:

molar volume at 475 K = 146.24 cm3/mol; M = 128.17 g/mol; Groups: 8 × ACH and 2 x AC2

From Equations (14)-(16): From Equations (17)-(19):

k I =0.152, k2 =0.529, k4= 1.064

From Equations (2)- (5):

~1 = -0.4719, ~2 = 0.9002

Tc

= 748.35 K; 60 = 0.3022

k0 = 0.712, k3 = 2.928, b = 237.63 cm3/mol e / k = 0.65Tc = 486.43 K; e/kT = 1.024

rio = 3.1163x10 -3 (128.17x475)1/2/(237.63) 2/3 = 0.02004 mPa°s From Equation (6):

g(ol) = 3.1016 and g(Rol) = 0.4564

From Equation (1):

rlpred = 0.389 mPa.s; error = +13% (from rlexp = 0.343 mPa-s)

Note that the overall error for naphthalene (0.48 ~

T r <_ 0.80)

AADL = 8.3%.