A mathematical model for surface shear viscosity measurements

A MATHEMATICAL MODEL FOR SURFACE SHEAR VISCOSITY MEASUREMENTS ignazio Becchi - Laura R e b a u d e n g o L a n d 6 * S O M M A R I O : Reeen/i risultati sperimenlali hanno evidenzialo l'influenza degli effetti della viscosita superficiale in correnti di ]abora/orio e di conseguenza l'importanza di poter conh'ollare tale grandezza lramile tell misuratore. Parlendo dalle leorie esislenti si formula un sistema eli equazioni differenziali /a cui soluzione numerica ~ confrontata con risultati analitici di altri AulorL In ultimo si presenta ]a proposta eli un misuratore di viscosith superficiale che rich&de solo misure di tipo cinematico. S U M M A R Y : A s it* recent experimental work ave have evidenced that the surface shear viscosiO, seems to be an important physical properO, in small scale free surface phenomena, we propose an anat.ytical approach for studying a surface viscometer. Starting from the existing thenO, on surface laA,er behaviom', we formulate a set of differential equation whose numerical solution is compared with a recent anaO,tical result. B3, means of the so calibrated method we anaO,ze the behaviour of a proposed meter .for surface shear viscosiO, which needs onO, kinematic measurements. An interface property like p, was originally introduced to explain the effect of adsorbed molecular layers, as the calming of sea waves by oil, which cannot be obtained by considering the bulk viscosity and the monomoleeular layer thickness only. The mechanical interaction between an interface layer and the surrounding fluids was early introduced by Boussinesq [3] and then extended by Scriven [4]. Following Scriven the dynamic equilibrium may be expressed, in the space (x, y, ~0, for a plane interface Z = cost by A,~= = T.. + (k + ,u.) (v..~ + ,,y.').. + ~,(v..,. - - v...).y Aa.,, = 7'.,, + (k + / s . ) (v..= + v'.')," + ~ . ( v ' . . - - v..')." A~,~ = 0 (1.1) where: T ~- the interface tension [kg m - q k ~- the coefficient of surface dilational viscosity [kg s m - q ZIg = the viscous traction across the interface [kg m -°-] 1. In unidirectional motion, as in uniform flow, equation (1.1) becomes, in.y direction, Foreword. In an experimental study on free surface laminar flows [1] the surface shear viscosity was introduced to justify a kinematic distortion from the theoretical behaviour. In fact the presence of this property induces two effects: a flow rate reduction and a velocity decrease near the surface, similar to those obtained from experimental results. Studying the free surface velocity distribution in turbulent uniform flows [2], we detected a typical behaviour whose significant parameters were influenced by fluid pollution and temperature like a surface property. To find out more about the knowledge of free surface mechanical behaviour we have studied a measurement technique which could be used to estimate the surface shear viscosity for an experimental fluid. The surface shear viscosity F, [kg s m-l], related to a two-dimensional scheme, may be referred to the bulk viscosity /, [kg s m -s] by means of a significant length - - Aa~u = / ~ v~,y. (1.2) Hereafter we will utilize 2 as the parameter for/z~ estimation. supposing no surface dilatation. Eq. (1.2) is practically analogous to the one introduced in reference [1]; following this equivalence is it possible to correlate the surface shear viscosity for free surface liquids with the shear viscosity of a monomolecular adsorbed layer. This observation may suggest that surface viscosity effects are relevant only for highly polluted fluids but experimental results show appreciable effects on surface viscosity even with slight pollution. Furthermore highly (1) polluted fluids may show a non-Newtonlan surface behaviour as a two-dimensional Bingham body or a pseudoplastic interface. This behaviour depends on chemical properties of both, fluid and pollution, as evidence by Joly [6]. Our purpose is to set up a mechanical method for evaluating the surface shear viscosity for a given liquid sample in order to know the surface behaviour in terms of temperature and to foresee surface viscosity effects. In this paper we analyze theoretically the flow field that will be established in a surface shear viscosity meter. * (1) The term "highly" is only a qualitative one because critical concentrations depend on the chemical nature of the solutes. )- =/~d/~ Istituto di Idraulica, Facolt8 di Ingegneria, Genova. MARCH 1975 49 2. on the fixed walls, and Analysis of the viscometric flow. The Couette flow with zero pressure gradient seems to be the best one for evidencing the kinematic effects of viscous forces. I n this kind of stationary flow the NavieroStokes equations are simplified, because the power needed for the motion is supplied only by the viscous forces [7]. Therefore the dynamic equation becomes, in unidirectional motion (Vv = vz = 0: 3p/c~x = 0): 12V-°v~ = 0 op/0y = ap/& o = -- on the moving walls, where Q is the angular velocity of the training body 12, 0 r Or r k - - ~ r ] - - / 2 ~ Q¢ = fluid density [kg m -4 s2] g = acceleration due to g r a v i t y - 9.81 [ m s --°] Though this motion is simple and stable in regard to turbulence inception, it is not possible to evidence the surface viscosity effects and, moreover, its practical realization is very difficult. In a Couette flow rotating about a vertical axis the equations (2.1) become, in cylindrical coordinates: = 0 (2.5) T h e finite differences n u m e r i c a l solution. To solve our integral problem in a bounded domain we have chosen the finite differences method. Following this method it is possible to test the stability and the convergence of the numerical solution of equation (2.2), even if from an analytical point of view it is not possible to prove the existence and the uniqueness of the solution for the complete problem. I n addition, this method enables us to represent any type of boundary geometry. Then, utilizing the finite differences criterion to the first order approximation we can represent the second of the equations (2.2) and the boundary conditions (2.3), (2.4), and (2.5) respectively in the forms: ,,, (4 + 2 Op PO - - - ~ o - Or ~ r l, + i 1 r Ovo - -Or -- l' 0 02PO r--Oyq- ~OZ (2.2) op - - v/.]+l - - vi,j-1 = 0 (3.1) v~,j = 0 (3.2) vm = £2.r,- (3.3) 0 g where t'z = vr = 0. Rotational Couette flows are frequently used for bulk viscosity measurements with a simple apparatus having two concentric rotating cylindrical tanks. This viscometer is useful because it needs a relatively small amount of liquid. Nevertheless, some limits in angular velocity are necessary when the inner cylinder is moving, because instability (Taylor vortices [8]) may occur in presence of three-dimensional disturbances. In order to estimate shear viscosity of adsorbed layers other Authors 15] have used rotating Couette meters, but they have not an exact analytical solution for the calibration of this apparatus. More recently Goodrich [9], [10], [11], [12] proposed some numerical solutions of a particular Couette flow in the presence of surface shear viscosity. Goodrich solved equations (2.2) in the half-space superiorly bounded by the Z = 0 plane, for a moving boundary of unit radius. As the method used by Goodrich cannot solve the dynamic equation in a limited domain, we have used his results only for checking the reliability of our method. The flow field of a Couette motion internal to a cylindrical tanks presents the following boundary conditions: v= 0 50 Ovo on the free surface, from equation (1.1) by neglecting the influence of the upper fluid viscosity. 3. p = isotropic pressure [kg m -2] OZ { Ovo l (2.1) where: cYaro , 0 = 12 ~ r (2.4) Vo = r,Q 0 +%r,)--- vi-1,] 2 (1----~r~)--v;.:-,=0 ,_/-~ (3.4) where : A = the dimension of the square mesh i = the discrete variable in r direction (re j = the discrete variable in Z direction = i • A) 2 / A = P = the dimensionless parameter for the surface shear viscosity. When P = 0 the equation (3.4) becomes: Pi,J-- lli,J-1 = O (3.5) As the equation (3.5) predominates over the others we have modified equation (3.4) in order to obtain for P - + 0 a classical symmetry condition: thus we have utilized a linear combination of equations (3.1), (3.4), and (3.5) obtaining, for the free surface, the form: (2.3) MECCANICA - - v,+,J ( P + + ) --v*-*.J(P + + ) (l+2-~-)-- (1--2-~)--v,.J-*=O (3.6) When P = 0, on the surface we obtain: v,., (4-k-r~-)- v,+,.s (1 + -~r~) -- v,-,., (1 -- ~r~) - - 2v~.j-1 = 0 (3.7) I n an n points discrete domain the equations (3.1), (3.2), (3.3), and (3.6) give an algebraic linear nonhomogeneous system of n equations; hence a numerical solutions may be obtained by computer. As A becomes smaller n increases and the numerical solution converges on the theoretical one, but too large a value of u makes system solution too heavy. To solve the numerical system we have chosen the iterative Gauss-Seidel method, which for a given computer allows the calculation of the maximum points number. In our case we reached several thousands of mesh points. As the presence of a second order condition on the free surface makes the system more difficult to solve, we have introduced the overrelaxing coefficient ~o [13], which improves the convergence of the method. The introduction of co transforms the equations (3.1), (3.2), (3.3), and (3.6) respectively as follows: - _ v,d o~ 1 obtained with 165 iterations utilizing coopt = 1.855, while introducing into the system co = 1.80, it was not possible to reach the solution after 8000 iterations! First we have applied the equations just obtained for the Goodrich problem named "rotating ring of finite thickness" [14]. In this way it is possible to compare our results with the ones obtained by Goodrich. As only a finite number of equations is consistent with a numerical solution we have bounded the integration field as shows Fig. 1. The domain so obtained is limited by: the free surface (OB), the symmetry axis (OD), the (BC) orizontal cut ( D C ) in Z = Z, and a vertical cut in r = r. The boundary conditions are: on OD, DC, BC eq. (3.9) in A(r = r0, Z = 0) eq. (3.10) on OA and AB eq. (3.11) ti - 1 F.--; i -I ' (4 ~--rlg_) {v,+,.j (1 + 2--~[) + v , - z . j ( - - ~ 7 7 r i ) + Fig. 1 - Our approximation to the Goodrich domain. + vt.j+l + vt.j-, t + (1 - - (3.8) co) v,'.j J vi,; = 0 (3.9) vi.j = Qr, (3.10) Since Goodrich domain was undefined both in direction "r" and "-Z", we have checked the influence of the rations r/ro, z/ro on our numerical solution. Keeping: vo(r,O) oo=lim r~0 -0" = co t ,,(" ÷+) (' ~ -'-2-)'v'+z"/P+I'+ \, 1o-' \ \ \ \ lo-2 \ \ \\ \ 1 +-v/b,,,,,- and 1 \ 60 (-Omax ~ " 1 --k.kA.~.,- ~lro ~ 1 The exact value of COoptis very important: for example, for a one thousand points mesh the solution has been MARCH 1975 ~'2 ~k + where the marked variable represents the value obtained from any new iteration. The best value of co may be found by solving the homogenous associated system: setting b=,, and bmax respectively the m i n i m u m and the maximum values for the ratios of a variable in two successive iterations, coopt is found to be the c o m m o n limit towards which tend: (Dmln 0 = - coo I (3.11) 60 and r 10-3 0 I \\fifo e\ 3 4 5 Fig. 2 - The influence of the fixed boundaries. 51 as significant parameters of the flow, we have checked the values of r/ro and I~]/ro rations behind them 0 is practically constant; in Fig. 2 numerical results obtained for ,~, = 0 show that for ;'~to > 4 and I~]/ro > 2 0 value change only by 10 -3 . Utilizing a domain where r / r o = 5 , I.zl/ro=2 and ~/ro = 0.1 we have calculated several solutions for different values o f the surface viscosity parameter P. In Fig. 3 we report the most significant results; the plot o f surface velocity for different values of P, and two maps o f equal velocity lines, for P = 0 and P = 50. The values of 0 obtained have been compared with those of Goodrich. In Fig. 4, where 0 is represented versus 2/ro, the continuous line plots our variable while the dotted lines represent Goodrich solutions for rings of relative thickness s/ro = 0 and s/ro ~ 0,002, that are the extreme values used by Goodrich. As for .~/ro > 0.3 the three plots are practically identical, it can be deduced that disagreements observed for smaller values o f 2/ro are due to the influence o f the dimension o f the rotating ring. In our case this dimension is not evident but it seems to be of the same size o f A. I. v~ r~r~ ,8 .6 .4 .2 1 t ~xk\\ 2 .~ ~ " X X ," . 3 " " x 4 I....... I t / /" 5 r/r~ __/ /I ~ I/ " - / / [÷ "k. -r-- .... i ~ \ - ~" ~ P=50 - r-- I Itltl IIIII ~ Ill I\l\l\l I Y3 V I tt I ~ I~ ~ ' I , - V V I A I il I I Ikl fN r~-"l.Xl I/I ~\H ~'~J I I f I I ~ I~,1 I I CI--F'FI L.LI~.I ]XL I I I I l I L~l\l I IN.I I I I I I ll~l k i I"f-.,LI I I I,_.U~\L..~I I I ITT~ I~ I ]",,,l 1 I 1 I I I IU~Itl I I"~LI I I I I L~ xJ ~ L_L\I l',,ll I I I I I I L.Lk ['KI I I I I I I I I~I\I l'l'-kl I I I I I l\l IXI I I 'r-f---~ i i l I I',,,.II",I..~I I I I I I I I f'~T"~ ' IIIIIIIIIII 1 I 111 t I I I U i [iiiiiit t i I III i I I I I I I I I I I II I I I I I 111 III l Jl 11 I l I I L/rJ t ~ V[ I ! ~ I/] t I ~ ] 1 t tl] I I I IA I I I I /I I I I I /I II / I I I I I/rl I I I vi I I I I I Y[ III VI I ,~ I I I l~rl i i i IA [ i i i i i I/I IJ..-Vq I I I I IA"I I I I I I I X/ I I/I I I/I I I I I I 1 I 1 IXI I I I I I I I VI I I 1,41 X I I I t I I I L,,'r I I i I I I I/V'l I I I ~ I/I I/I I I I I I I 1.1111 I I i I I I I I Jill I I I J/l I .4 L..k-4'q I I,,I/I I I I I IJ/l I I I IXl I I I I I I I I II I I I I l,,'r'll I I I LA"I I I [ I.¥I I I I I I I I I I I I lJ..4"l ILl I IX" I 1 11"I" I I I I I I I I IJ..-4"~l I I ] ]~.l'q I I I l...lll I I t ~ t ~.-..~;-q i I I I I ~ I i i i i J...4-.q i I I I I I I I IJ..-D--I-"4"TI I | L_..L-~ ' ' i'~l I ~ I I I I 1 illllllllllllllllllllllll II II IIII I Fig. 3 - Some typical results for the domain of Fig. 1. 52 MECCANICA I I I --i t I "-I .5 - - t II II II II I1 I 11 IIIII IIIII IIIII IIII tilt L,--I-4-H-VIiII 0 , 10-4 tO-3 I I I ! I IIIIIIII IIIIIIII III IIIIIII1 lt111III /Y Iit l lltl~" 11111 II11t III I111 t 1 1111 1111 1tl 111 10-s I ttNll llltl 10-2 IIIlll I I I I IIII I 1t I I1 1 111 1 1t IIII Illi till III11 I Ii I l I0 IIII1 tllll Io2 ,t/to Fig. 4 - Comparison between our solution and two of Goodrich (dotted lines). Fig. 5 - Effect of ring thickness for different edge shapes. If the contact surface area equals zero, when/z, is zero, no power can be supplied for the rotating ring with finite velocity gradient. Therefore when ~./ro tends to zero the influence of the ring thickness becomes appreciable. For a better knowledge of the relationship between the mesh size zl and the equivalent ring thickness we have carried out some calculations with thickned meshes. As the number of points which can be introduced in our calculation is four thousand at maximum, the minimum value of A Iro that may be obtained is 0.005, while the Goodrich maximum value of s/ro is 0.002. In order to analyze flow behaviour further, when P is zero, varying the form of the moving ring, we have set up a new numerical method utilizing a procedure first introduced by Viviani and Centurioni [15]. Fig. 5 shows the values of 0 versus ro/sfor four different solutions : c) 11 s ro .I l b) U S a) i) Goodrich results (white points) ro I I ii) our square mesh (black points) (cf. Hg. 6a) iii) two solutions obtained by means of the above mentioned method; the upper for a balf-torus shaped moving ring (cf. Fig. 6b), the lower for a Goodrich-like boundary (cf. Fig. 6c). MARCH 1975 Fig. 6 - Geometrical schemes of the different edge shapes. 53 ++)++1'- As appears in Fig. 5, in the case of square mesh the equivalent dimension of ring thickness may be evaluated at about one half of the mesh size. v~+l,j 4. T h e behaviour o f a surface shear viscosity meter. On the basis of the above mentioned results we have worked out a system for estimating surface viscosity by -kinematic measurements. Figure 7 shows a practical design for a viscometer consisting of: i) a cylindrical vessel having internal radius L and liquid height H adjustable between 1/10 and 1/5 of L, ii) a rotating cup whose sharp edge skims the liquid surface at a distance R from the symmetry axis, iii) a disk, of radius r, freely rotating coaxially with the cup over the liquid surface. Therefore to estimate the surface shear viscosity it is enough to measure the angular velocity o f the disk. The behaviour of this apparatus may be interpreted by equations (2.2), only if the flow does not attain Taylor vortex conditions. For any liquid these conditions result in a upper limit for the angular velocity of the rotating cup that should be determined experimentally. *(* + ½) (P + 1 k = 0 (4.2) k + 1 for i = k = r / A , j = H/A. Equations (4.1), and (4.2) with equations (3.1), (3.2), (3.3), and (3.6) form a proper system for solving the flow field. The solution is once again obtained with the iterative Gauss-Seidel method utilizing the overrelaxing coefficient. By this way the set of equations used are: for all the internal points eq. (3.8), - - for all the points on the fixed boundaries, enclosing the symmetry axis, eq. (3.9), - - for the cup sharp edge contact point eq. (3.10), - - for the points on free surface eq. (3.11), - - along the disk the bulk rotation condition (4.1) in the form t¢ //m,j vj.j = r, ~,, - 1 Fm (4.3) on the edge of the disk the equilibrium condition (4.2) in the form vl¢,l -~- 1 + - g - v~,j-1 - - k+l -- 7. r m ( V , . d - - V,.,]-I) k-I 1 1 + (1 - - ,0) V~,j As the disk boundary conditions make numerical stability more difficult, the w value increases and depends on P values. For a two thousand points mesh, we have obtained the approximate relation ~t Fig. 7 - Proposed surface shear viscosity meter. For studying the flow the square mesh method is the only one possible, because the other methods mentioned above cannot represent the boundary conditions on the disk. In addition the Goodrich method is not suitable in presence of immersed fixed boundaries and the finite differences program carried out by Molinari and Viviani [16] cannot calculate the surface conditions when P differs from zero. By introducing the finite differences scheme, the boundary conditions on the disk may be obtained: v,.j/ri = vi+i.~[r,+l for and k-1 7.ra rm(vmd-- rind-x) q" 1 54 0 < i < k, / -----H / A (4.1) = 1.7 + 0.1 e~.sP Fig. 8 shows the results so obtained; the surface velocity distribution for several values of P and two maps of equal velocity fines for the extreme values of P. Then, calling 0o the ratio of the disk to the moving cup angular velocity, in Fig. 9 we have plotted 0n versus 2]L for four different situations, varying H and r in respect of L. Comparing the obtained behaviour with the ones shown in Fig. 4 we may deduce that the closeness o f the vessel bottom shifts the sensitivity range towards higher values of 2[L. An increase of the disk radius balances this effect but at the same time reduces the sensitivity bandwidth. From the calibrations o f point three it is possible to give an equivalent dimension to the sharp edge o f the cup, s = 7.14 10 -a L. Errors in disk level control give relative errors for 0n values smaller then themselves, but their sign cannot be foreseen. In addition equivalent errors in the cup vertical position give for 0o value even smaller differences of the same sign. MECCANICA I, k ve ?ifi .8 .VL: ~ o / / .6 .4 k f ,I .2 i o j J J ,1 .2 IkJlli ii i~ i P,,i i : I I"I..I I ' ' : ~ I -al~:~ " - - T I I i , i. I I ¸it)! ) Ii II\ .4 .5 ~-".-._i i 'iii LII ;IIIII ~ I .3 i---._ I ~-~ --~ ~" . . . . . ~ .6 ~ .... i , ~ - ...... l , .7 - , - . .8 i . ': -'-l--'--i. . . ~ I mi lmmnmmmmmmmm~mmmmmlmmmmmmmmmmmmmmii Immmmmmmmmmllmmmmmmmmmmmmmmmmmmmm, iMmmmmmmmmmmmmmmmmmmmm~mmmmnmmmmm Im~mmmmmmmm)mmmmmmimmlmmmm~mmmmmm Imm~mmmmmmmmlmmmmumm~mmmMmmmmm ImmmlmmmmmmmlmmmmmiimHmmBimmmm Immmmlnimmiimpmmmmmliiml)mmmwami! Immmmmi~mmmmmmummmmwmmmwmmmmismm lmmmmmmmlammmmmlwmmmmtammt~mmmmi~ immmmmmmmmla!mmma~ammmmm)mmw~mmmm Immmmmmmmmmmw~mmmmmmnBmmw~)mi~a) ummmmmmmmmmmmmm~'~,~)mmBww~mm!mm~=~m Imli~e)Mmmmmmmmmmimi.e)~mw.~)mmm -m i ) " r/L i - .~ , _ _ ~' - - ! .9 ,, I - : f~ P- ~b j .-', i + -~--VI I rill II L, , Fig. 8 - Some typical results for estimating meter flow fields. 1, I0 -s 10-2 10-I Fig. 9 - Performance of surface viscosity meter for R/Z. = 0 . 8 r/L=0.4, H/L=O.1; MARCH 1975 1 10 I0 e 3./L I0 a and: I: r/L ~ 0.4, H/L = 0.2; II: r/L = 0.2, I~/L = 0.2; ] l ] : IV: r/L=0.2, H/L=O.1. 55 Conclusions. 1) Our finite differences method gives a numerical solution for a fluid flow having the surface shear viscosity as a boundary condition. 2) By means of several tests the numerical solution has been improved by introducing the overreta:dng coefficient and by checking the mesh size influence against some known solutions. REFERENCES [I] I. BECCHI, and G. SEmN^RA,1971, Udine I ° Congr. Naz. AIMETA, voI. III, pt. I, pag. 75. [2] I. BECCHI, and L. REBAUDENGOL^ND6, 1973, L'Energia F3le#rica, n ° 9, vol. I, pag. 1. [3] J. BousslNESQ,19t3 Ann. de Chim. Phys., vol. 29, pag. 349. [4] L. E, ScRivm,,',1960 Chem. Eng. Science, vol. 12, pag. 98. [5] A. G. BRowr,:, W. C. THUM^N, and J. W. McB^I~, 1951 J. of Colloid Sciences, pag. 491. [6] M. JOLY, 1950J. de Physique et Radium, vol. 11, pag. 171. [7] G. SuvtNo, 1949, Lincei, Rend. Sc. Fis. ~'[at., vol. 6, pag. 708. [8] Ho SCI-ILICHTING,J3otmdary Laver Theory 1968 ~vlcGraw56 3) The finite differences equations has been proved suitable for solving the flow field for the proposed surface viscosity meter, whose efficiency was demonstrated in a wide range of surface viscosity values and for several geometric dimensions. 4) The performance of this meter seems to be acceptable also in the light of any incidental errors that might occur in practice. Received 9 September 1974. Hill, ch. XVII, f, pag. 500. [9] F. C. GOODRICH,1969 Proc. Roy Soc., A, vol. 310, pag. 359. [10] F. C. GOODRICH, and A. K. CHA1"rERJEE, 1970, J. of Coll. and Interf. Sci., voL 34, pag. 36. [11] F. C. GOODRICH, L. H. ALLEN, and A. K. CHATTERJEE, 1971 Proc. Roy. Soc., A, vol. 320, pag. 537. [12J F. C. GOODRICH,and L. H. ALLEN, 1971, J. of Coll. and Interf. Sci., vol. 37, pag. 68. [13] R. S. VARGA,Alatrix llerative Analysis, 1962 Prentice-Hall. [14] F. C. GOODRICH,and L. H. ALLgN, 1972 J. of Coll. and Interf. Sci. vol. 40, pag. 329. [15] L. Cm,~TURIONI,and A. VlVt^N[, 1974 Ist. Elettr. Univ. di Genova (to be published). [16] G. MOLtNAaI, and A. Vlvr^m, 1974 Ist. Elettr. Univ. di Genova (to be published). MECCANICA