ENGINEERING RESEARCH INSTITUTE UNIVERSITY OF MICHIGAN ANN ARBOR THE COUETTE FLOW BETWEEN TWO PARALLEL PLATES AS A FUNCTION OF THE KNUDSEN NUMBER C. S. WANG CHANG G. E. UHLENBECK Project 1999 OFFICE OF NAVAIL RESEARCH, U. S. NAVY DEPARTMENT CONTRACT NO. N6onr-23222 June, 1954

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN TUE COUETTE FLOW BETWEEN TWO PARALLEL PLATES AS A FUNCTION OF THE KNUDSEN NUMBER I. INTRODUCTION In a previous report' a formal solution was obtained for the heat transport through a gas between two parallel plates as a function of the Knudsen number d/k (d is the distance between the plates, X is the mean free path of the molecules). The starting point was the linearized Boltzmann equation, since it was assumed that the magnitude of the disturbance from equilibrium, measured by the ratio AT/T CAT = temperature difference between the plates, T = average temperature), was small. In the following, the same method will be applied to the problem of the Couette flow of a gas between two parallel plates as a function of the Knudsen number d/x. We will again assume that the Mach number which in this case is the ratio of the average flow velocity to the mean molecular velocity, is small, so that the disturbance of the equilibrium due to the moving plate is also small in this case. Since the method and the general features of the solution are quite similar to those of the heat-transport problem, only an outline of the calculations will be presented in sections II, III, and IV. In section V some divergence difficulties will be discussed which also occur in the heat-transport problem, and which are due to the parallel-plate geometry. 1C. S. Wang Chang and G. E. Uhlenbeck, "The Heat Transport between Two Parallel Plates as Functions of the Knudsen Number", Univ. of Mich., Eng. Res. Inst., Project M999, Sept., 1955. 1

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN II. FORMULATION OF THE PROBLEM We take the y-z plane halfway between the plates. The upper plate, situated at x = d/2, is stationary, while the lower plate at x = d/2 is moving in the z direction with a velocity w (measured in units of 4m/2kT). The velocity, w, is assumed to be much smaller than unity so that terms of the order w2 and higher will be neglected. The notation of the previous report will be used. The distribution function is written as f-=F. [lt 0 (I C Xw)0 where C = c - w/2, n is the equilibrium-number density, and T the equilibrium temperature, which are all constants. The Boltzmann equation is again > =3(F;) l) 3X X where J is the collision operator. To formulate the boundary conditions w introduce the distribution functions for the molecules going up and down, i.e., f+ and f, and the corresponding disturbances h+ and h-. In terms of the h's, the boundary conditions are r C,(g + c w- n t-(- e- x XdO (2n and etend = (-tied c) X <iob) where J is the accommodation coefficient and the constants B+ and B- are to be determined by the conditions

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN expressing the conditions that the number of molecules per unit area between the plates is given (equal to nd) and that in the steady state there is no net flow in the x-direction. Up to terms linear in w we have the symmetry condition t (C,,?L) - ht ( ) ), (4) and consequently With the symmetry condition (4), Eq. (3a) is automatically satisfied, and the boundary conditions (2a) and (2b) are equivalent with a+ = - C = B. Thus our problem reduces to the solution of the integral differential equation (1), subjected to the boundary condition (2a) and the auxiliary conditions (3b) and (4). The physical quantities we are interested in are the drag on the upper plate, the velocity distribution as a function of x, and the velocity jumps at the plates. III. GENERAL SOLUTION At first glance, this problem seems to be more complicated than the problem of the heat conduction because of the lack of axial symmetry around the x-axis in the present problem. It will be seen, however, that this is not a real difficulty. It is true, though, that in the development in eigenfunctions, we must now use the spherical harmonics instead of the Legendre polynomials used in Ref. (1). We write

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN where e is the angle between the velocity C and the x-axis and 0 is the azimuthal angle measured from the y-axis, and develop h according to Because of the reality of h one must have arlm = arlm. From the symmetry of the problem with regard to the y direction, it follaws that h must be an even function of cy, which has as a consequence that the arm for odd m must be pure imaginary The first few of the development coefficients arlm are related to quantities with physical interest. For instance, the density and temperature are found from (~) = m [-+ _r_3/z ~ooo -(x)= T I- 3- >3/) and the physical quantities (shearing stress, average velocity, and heat flux, all in the z-direction) in which we are especially interested are given by -42 ' toss (C | CfX _ 73b 3- wow ' ci- -l, (tb) From the conservation laws it follows that cx, P Pxxy, Pxz and qx are constants. We will require that cX = 0 (see eq. (3b)); Pxy is zero by symmetry. This implies that* CIO- CA0 C 010= ~t ozo _ O1 oo0 ~ N hco *ozo 0 m in i ton R.10 = C1yL)~. *There are some misprints. in the corresponding formulas in the previous report, Bef. (1)

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Introducing the development of h into the Boltzmann equation, we find that our problem reduces to the solution of an infinite set of homogeneous linear differential equations: _ = n X _,,/_ ' [ ()l with the boundary conditions: 2. + (4 a)C c AQ _+( CA)i icA csV) Finally, a010 = O, and from the symmetry property eq. (4) it follows that even odd Jarm ein x according as I - m is) l odd even. The square bracket in eq. (7) is different from zero only when m = m', and l - l' is odd. Thus coefficients with different values of m are not coupled. Since forjlm>l the integral in the boundary condition (8) is zero, it is clear that the equations for arlm with ImJ>l are completely homogeneous, so that the only solution will be zero. This is also the case for m = O0 because in eq. (8) the quantity B is also an unknown and the inhomogeneous part is only the term containing wCz. This term is also zero for m = O. Hence one may conclude* C -) t7Y = ~) for m / + 1 We use the same procedure as was employed in Ref. (1). One must distinguish between even and odd values of l. Eliminating ar 2L + 1)1 one obtains an infinite set of second-order homogeneous differentia equations for ar(21)l, which, with the use of the symmetry conditions, has as solutions: *This is in contrast to the case of heat conduction. It is also physically evident that a temperraure gradient will piroduce a change of the density n neara the plate, while this is not the case for a velocity gradient. 5

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where the pi's are the positive and nonzero solutions of the secular determinant and where in the sets of the br(2b)l one set can be taken as arbitrary, for which x(i) = (02 12 )z ' foI ) b02 =0 for all i. The coefficients r(21 + 1)1 are now completely determined and one finds fromPutting eq. (.]) in the boundary conditions, one this one needs first the set of inhomogeneous linear eobquations: I*-9 \ 0 C LI - I I J where Putting eq. (11) in the boundary conditions, one is led to the set of inhomogeneous 1inear e quations where the known constant matrices B are defined by whee te kowncontan marics 8aredefned 6

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN In the same sense as in Ref. (1), eq. (12) is a set of inhomogeneous linear equations with the same number of equations as that of unknowns. Thus the problem is formally solved. The complete solution involves infinite determinants. We are, as yet, unable to discuss the convergence of these determinants nor the convergence of the "breaking off" processes which we will use as in Ref. (1). IV. THE MAXWELL MOLECULES We will use the same "successive approximation" method as in Ref. (1). Furthermore, we will restrict ourselves to the Maxwell molecules so that and t(3) t.,,(jjtl\ 0t0 =o unless both r and I are zero which has as a consequence that all br(21)l are zero except bo l. Since we do not expect to get anything from the zeroth approximation, we start with the first approximation. A. First Approximation As in the heat-conduction and sound-propagation2 problems, we will use in the first approximation the eigenfunctions for which 2r + ~ ' 3, and 1<3. since only the I's with m = 1 enter the problem, we therefore need to take only ~oll, V*ll1, and 'o21. We know that (o) a021 = const. = bo21. ((3) 2C. S. Wang Chang and G. E. Uhblenbeck, "On the Propagation of Sound in Monatomic Gases", Univ. of Mich., Eng. Res. Inst., Proj. M999, Oct., 1952. ii ~ ~ it. e ~ ~7

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN The other two a's are bor\ sx 2-1 0t [ I l j3-) Since p = 0 in this case, the boundary condition (12) reduces to (0) (0) 02S -;-( -)b ~ _ T02 where -6- \A- + _ I K +)()c -+)! Using the values of S and T as given in tables in Appendix II, one finds where. (. Using eqs. (6) and (13), one obtains the fo.llowing results: and - 2_ XC W < I- 0( &\FT - (\4 P)

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where w is still in units 22kT-m. Therefore the velocity slip at the upper plate x = d/2 is ( x 3 \+LX 8H For small X/d and ca = 1, this is in agreement with the result of Maxwell. B. Second Approximation (2r + I x 4 and 1 < 4) The eigenfunctions entering the problem are -r.ol-, 021, 111i, 03os1 and. r121. The 8's with odd values of I - m are *021 and *121. The only nonzero value of p is The constants to be determined by the boundary conditions are bo2)1 and. 2 Equations (12) become a set of two equations: and Solving these equations, one then can calculate the quantities of physical interest. It turns out that in all approximations the drag, the velocity distribution, and the velocity jump have the form Pvx- 2 --- a- ______ _ z\T X8 3 - x 11 3 2) C - L 3 (5 v 9

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN where the sum is over all the roots Pi, and si = sinh (npid/2). The functions A, B, and Ei are all functions of K = d/X. In this approximation, one finds I6\ = IS- C where t = tanh (npd/2), c = cosh (npd/2). C. Third Approximation (2r + l - 5, and ~ < 5) The eigenfunctions to be taken in this case are 11o,, t021, k11, t031, *121, I041, *211, and *131. There are four 8's with odd values of I - m. Hence there are two pairs of p differing from zero. They are Pi = 0.7151 A2 and P2 = 1.342 A', where Al = A2 /kT. We find, further, bo(l -0.2456 b () b(2 2539 bz121 Solving the set of three equations for the boundary conditions, one finds then, A= —1 1 = (~ 2.o ) t, (O. 38 t A o i3>t+) a 0331A t2t - <0 t- c~( (os-3 67t o5o7 t0 ) o. t 1 tt > Em= C ( <d O, el 6g t,) E = | 0_ *. o.3oo0+0o173-t1 ) C710 Zo( lO ~

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN D. Fourth Approximation (2r + _ 1 6, and I < 6) We -roceed as before, with the four values of p's and the corresponding ratios of br(21)1 /bPlq given in the following table. i 1 1 2 3 4 p(i) o.4604A, 1.065A, 1.261A | 1.610A2 041 121 (- 5 525) (.2) ((15.58) b0)/b a 1 12 1.218) 16.204) (1.213) (8.442) 141 121 1 12 12 12l 121 | 46 (1.3j2) | 4 (0320) | (-5349)| 4 (32-51) 4 4 4 111 Taking a = 1 to simplify the computation we arrived at; the following results: = a 3g t a+ 3 t,+ a.73 t+ a.sot t a4t< t3 t,t~t D 0 -t,t. 7 titQ + 1 t~ 23 21 ttt +,6 \ t + 3o3 t, tt t3 t - 3 t03 ttztt t ott3t4t t.14 t23,t6, t 3.10 ttt3t4 ~B = 38+ a -t~t+ 137tft 1 t, s-33t4 + 4U-ttt 2s L -t tj+ t 0.4otst~t -t t+, t atit ttt4 4 Z 4 t, t,,t ~ 2.A *r t&bt, t A t433 ttttt ), 2 -~( lS7 t -R i/ t tl- 3. 33+t, t 37 ot1 f 3 r7Lt~t t D ~oQl it,t,4 t 3 '13t3 ~3t4* 3. s&,tt3t4), 11.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN E- 0331 (0.805o o '7t-t oi,.tj oU!33t4- o~1-qt,t~tf O 1-ttIt -t o. oo-t- - o.t4 t-,~),. o. 'G (:,Sc l-:.Loqt, l z.qtSt.:z.70oS,.t ~. -St'tt 1.6rt,t3 -t + 6t3 t C tt tst t S~), For the discussion of these results we have plotted three sets of curves. Figure 1 contains a set of four curves for Pxz/(Pxz) Knudsen against the Knudsen number, d/X. Curve I is for the Stokes-Navier approximation. Curves II, III, and IV are results from our first, second, and third approximations respectively. We have not plotted the curve for the fourth approximation because we do not expect anything new and the numerical computation is very laborious. The straight line on the left is the initial slope for the exact solution. All the curves from the different approximations have the same value (unity) for K = 0, and for K equal to infinity they all approach the same limit, zero. The initial slope for the four approximations are listed below. Initial slope Exact -1.242 First approximation -0.423 Second t" -0.458 Third it -0.534 Fourth " -0.598 Figure 2 consists of a set of three curves for the first three approximations for the velocity slip at the upper plate as a function of K. The initial slope is given by: Initial slope First approximation -0.423 Second "I -0.494 Third " -0.644 Fourth " -0.755 We have drawn Fig. 3 to show the velocity distributions as functions of x for K = 10. Since the velocity distributions for the different approximations deviate very slightly from the straight-line distribution of the first 12

1.0 0.9 0.8 -; I STOKES NAVIER 0.7 \ FIRST APPROXIMATION \ m SECOND APPROXIMATION & l THIRD APPROXIMATION 0.6 Pk L= 0.5 0.4 INITIAL SLOPE 0.3, 0.2 0.1 I I 0 2 4 6 8 10 12 d/X Fig. 1 The Drag as a Function of the Knudsen Number 13

1.0 0.9 0.8 0.7 W 0.6 0.5 -,,,I FIRST APPROX, SECOND APPROX. < THIRD APPROX. 0.4 0.3 0.2 0.1 I I I 0.0.2.4.6.8 1.0 1.2 1.4 d/X 2AW Fig. 2 Velocity Jump, W as a Function of the Knudsen Number w~~~~1

0.0280 0.0240 0.0200 0.0160 0.0120 - 0.0080 / 3rd Approximation IOO40,/ 0.0040 - 2nd Approximation \- ~2 0.4 0.6 0.8 1.0 2x d -0.0040 Fig. 3 Velocity Distribution in the Second and Third Approximation as Compared with the Linear Velocity Distribution R1 (x) From the StokesNavier Approximation with slip for d/\ = O 15

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN approximation we have drawn instead of cz against x, the curves Ri - R, where R = 1 - 2cz(x)/w and R1 is the value of R from the first approximation, namely, As in the case of heat conduction, the boundary effect is shown by the sharp rise of the curve near the wall. V. THE APPROACH TO THE KNUDSEN LIMIT; A DIVERGENCE DIFFICULTY In Ref. (1) we pointed out that the behavior of the exact solution for the heat-conduction problem is for large K quite different from the behavior at small Knudsen number. The approach to the Clausius gas limit (K>>l) is complicated by the occurrence of the hyperbolic functions, so that a development in inverse powers of K is not possible. We see from section III that the same is true for the Couette flow. On the other hand, it seems that the approach to the Knudsen gas regime (K<<K) is quite regular so that a development of all relevant physical quantities in powers of K should be possible. In fact, in Ref. (1) we gave the first two terms in such a development for the heat flux. Analogous results can be found for the drag*. In the zeroth approximation one finds 0)B t4T ic V ( with w still in units of 2kT/m. In the first approximation pa -T) I( I Cg c,(, C 4$4r, (17) For ac = 1 this reduces to the expression given in Ref. (3).3 There the square bracket which is always negative has been evaluated for elastic spheres and Maxwell molecules. The results are as follows. *These results can be found either by the method described in Chap. III of Ref. (1), or from the general solution. 3C. S. Wang and G. E. Uhlenbeck, "Transport Phenomena in Very Dilute Gases" CM 579, UMH-3-F, Univ. of Mich., Eng. Res. Inst., Proj. M604-6, Nov. 15, 1949. 16

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN For elastic spheres of diameter a, e isO 6FE), (ist) 0) For Maxwell molecules, -f\ = -1 i C so-F(e) I r-vzi O -. +)4 } In this case the result obtained from the general solution is expressed in the eigenvalues kl of the collision operator in the form A =__,_ - )j 49tl) ttQk)Lrt-+A!J (K) The summation can be carried out (for details see Appendix I) and confirms eq. (18b). The integral was evaluated numerically in Ref. (3), with the result where )\ = rrsA; <4. Thus for a = 1, the exact value of the initial slope in Fig. 1 is -1.242. A difficulty appears when one wants to calculate the second approximation of the drag, or when one wants to find analogous expansions for the average velocity distribution or for the velocity jump at the plates. Formally, one finds pxa =P~z~ mZxa i) M lx ' (a 4'Ccz CC r X (2- a and in the first approximation, C j C1 = -3/L Ce The bracket expressions are both divergent integrals because of the factor l/cx. The same difficulty occurs in the heat-conduction problem, although we did not notice it at that time. The second approximation to the heat flux given in 17

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN Ref. (1), p. 8, eq. (19) is also divergent, and the temperature distribution in the first approximation for which one can derive the formal expression T'< =_A'y (Ca; at (- ) 'X, ( ), %, 0L~- O( 1 (1) is divergent just as c The divergence of all these integrals is of a logarithmic nature and is always due to the factor 1/cx. The origin of these divergence difficulties is therefore clearly the parallel-plate geometry. Molecules which are emitted nearly parallel to a plate will have to transverse a very long path before reaching the other plate, so that the Knudsen approximation will not be valid for these molecules. Or one can say, the average Knudsen number for molecules emitted in all directions is logarithmically infinite even if d/X<C<l. We have found that by taking two concentric spheres (radii a and b) instead of two parallel plates all divergences disappear. For d = a - b<<K/2)(a + b) = R, quantities such as the temperature distribution will in first approximation in K = d/X contain terms proportional to ln (R/d), which blows up for the parallel plate geometry. One may conclude, therefore, that these divergence difficulties will not affect any real physical situation, and that the behavior of the solutions near the Knudsen limit can still be considered to be regular, so that a development in powers of K is possible. 18

ENGINEERING RESEARCH INSTITUTE UNIVERSITY OF MICHIGAN APPENDIX I DERIVATION OF EQUATION (18b) FROM EQUATION (19) In terms of the eigenvalues,l- of the collision operator, the first order correction to the pressure Pxz in the Knudsen limit is given by (bAn~~) L- o( C2~~~~~~ 3) The expression for Xr2L is In the braces, we see that the second, third, and fourth terms are obtainable from the first term by replacing 9 by t - 9, O, and it/2. Thus, substituting eq. (24) into eq. (23), we see that we need only to calculate the sum Putting r + ~ = i, we have The last sum can be done with the help of the relation C-X=Xt 2 Iue9) T '( ) 19.

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN leading to = 25{ (\tt6~) 8W\,a ~. IB __ _ r~it 4LIr1+ L. A _ r\+-tL j L 98), '~ 7,ci-j and thus to eq. (18b). 20

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN APPENDIX II TABLE USED FOR NUMERICAL COMPUTATIONS 1. Table of Eigenvalues. (A = A, A4 = R = o.636) ~ 1 2 3 4 - A, 9_Al 7 Al_(1 5R 0 | 0 4 ] - s~z | ~8 A2 -2 1 C 2 | A0o 11 A (l1 5R 2A.(l - ~2 8 8 2 Aa 3A - 9 i R. 2 2 9-Al (1 - l. 2. Normalization constants: g) _ ol (>+ 0 )m) N2~.211 21 31 41 51 3 2 2 22 4 - 2 22 2 _3 9 45 3-5.7 3 *5 '7 2 2 21 2 ~t2 5.57 '2.2 4a 2' 22 2 5'7 34'7

ENGINEERING RESEARCH INSTITUTE * UNIVERSITY OF MICHIGAN L EQ,,YT I W}/ e C 4W CA IIM rL 011 111 001 211 131 051 021 3S2No11No21 0 0 0 0 0 3 3 121 |352NO11N 121 | L5 ItN121 0 0 0 0O 041 o 10n2N0j1%41 10ixNl1N041 10~52No31N041 0 O 0 3 15 3 3 2 jN11-1 1N141 llNO31N141- 35yc2N23N141 4 GtN13 41 r2 | 021 121 041 221 141 021 2 9i 15ir 27?t - Nl1No 21 1021 1TNO21N1N1 o221N021 04..2..9..13 2.5..No4.. |t. 3 2'1 121 y NO21N221 | - -N12N1 y f No41N21 | 1N221N121 5tN141N2 0.2 4'-1o4 16 2 1441 | 3 *3NNl 12N1 | No4N14 3 21N14 2o 4.N1 4 -'021 8 4 31 -- No l252 2 2 22

ENGINEERING RESEARCH INSTITUTE ~ UNIVERSITY OF MICHIGAN /dde - _~_ _,' ~?,~___,_re 021 121 041 221 141 3iPg -3ic 5i -9ig N2 1 5il 141 3-4 No021 -8 N121 4 No041 25 28 N4 23