- Email: [email protected]
Numerical solutions of non-linear kinetic equations for a one-dimensional evaporation-condensation problem
Computers & Fluids, Vol. 1, pp. 367-377 Pergamon Press, 1973. Printed in Great Britain
NUMERICAL SOLUTIONS OF NON-LINEAR KINETIC EQUATIONS FOR A ONE-DIMENSIONAL EVAPORATION-CONDENSATION PROBLEM* S. M. YEN Coordinated Science Laboratory University of Illinois at Urbana-Champaign Urbana, Illinois 61801, U.S.A. (Received 11 December 1972)
Abstract--Numerical solutions have been obtained for both the nonlinear Boltz~ann equation (for two collision laws) and the Krook equation for a one-dimensional evaporationcondensation problem for a range of parameters. Our calculations of the evaporation (condensation) rate indicate that the linear prediction underestimates the non-equilibrium hindrance effect due to intermolecular collisions. The Krook evaporation rates are lower than the corresponding Boltzmann results and, thus, differ even more from the linear values. The Boltzmann evaporation rates for a gas of hard spheres are lower than those for Maxweilian molecules at lower values of Knudsen number. The micorscopic and macroscopic properties obtained from the Krook solutions differ appreciably from the corresponding Boltzmann results. 1. I N T R O D U C T I O N The transport phenomena occur in a gas because of the departure from thermal equilibrium. When the departure is small, the transports may be expressed in terms of the transport coefficients and the property gradients. For a given problem, the property gradients are obtained by solving the basic continuum gas dynamic equations with known relation between the transport coefficients and a property such as the temperature. Since the physics of the transport phenomena is considered implicitly through the transport coefficients, we face only one problem, namely, a mathematical one of solving the continuum gas dynamic equations. Under conditions far from thermal equilibrium, the concept of transport coefficients is not valid and the continuum equation is no longer applicable. The gas under such conditions can be described by the Boltzmann equation, the basic equation in kinetic theory. This equation contains the term collision integral which represents the rate of change of velocity distribution function (or molecular population density) due to intermolecular collisions. The collision integral which depends on the local departure from thermal equilibrium is essentially the microscopic transport property. The solution of the Boltzmann equation yields the velocity distribution function from which any transport property can be evaluated. In non-equilibrium gas dynamics, we, therefore, face two problems: (1) physics: to find the effect of intermolecular collisions on the microscopic transport, and (2) mathematics: to solve the nonlinear Boltzmann equation. * This work was supported in part by the Joint Services Electronics Program (U.S. Army, U.S. Navy and U.S. Air Force) under Contract DAAB-07-67-C-0199, in part by the National Oceanic and Atmospheric Administration under Grant N22-1-72-(G), and in part by the Office of Naval Research under Contract N00014-67-A-305-0001. 367
368
S.M. YEN
The advent of the high speed electronic digital computers has spurred the development of methods of solving equations of kinetic theory, including the Boltzmann equation. The major difficulty of solving the Boltzmann equation is the intractable form of the collision integral in it. Straightforward numerical quadrature would require days on the fastest computers today. The development of Nordsieck's Monte Carlo method[I] has made possible the direct solution of the Boltzmann equation for several non-equilibrium gas flow problems, including those of heat transfer[2,3]. His method consists of two distinct steps: (1) to evaluate the Boltzmann collision integral by using a statistical sampling technique closely resembling the real statistical collision phenomena in the gas and (2) to solve the Boltzmann difference equation by using an iterative scheme. The linearized kinetic theory problem of vapor condensation on or evaporation from spherical liquid droplets under non-equilibrium conditions has been tackled by Springer and Sampson[4], Shankar[5], and Edwards and Collins[6] using the moment method and by Pao[7] using the Krook equation. Springer, et al. used Lees' Ansatz; Edwards, et al., Grad's expression. These efforts have also led to the study of a simpler one-dimensional evaporationcondensation problem which, in essence, is the heat transfer problem of rarefied gases with net mass flux between two plates at different temperatures. The significance of studying this conceptually simpler problem is twofold. First, it is possible to solve this simple problem more accurately using basic equations of kinetic theory under a large range of nonequilibrium conditions. Secondly, the solutions for this problem would be relevant to the problem of more complex geometry and the study of their behavior would lead to better understanding of the non-equilibrium gas behavior for evaporation-condensation problems in general. The moment method with Lees' Ansatz has been used to solve this problem by Springer and Patton[8] in the linearized case and Edwards and Collins[6] in the non-linear case. Pao[7] solved the linearized Krook equation. We have solved both the non-linear Boltzmann equation (for hard spheres and Maxwellian molecules) and the non-linear Krook equations for the one-dimensional evaporationcondensation problem under a wide range of conditions. Our solutions yield the velocity distribution function, the collision integral for the Boltzmann solution as well as their moments of physical interest. 2. FORMULATION OF THE PROBLEM We consider the one-dimensional problem of rarefied gases between two parallel plates having two different temperatures, T1 and T2, and at distance d apart. This problem differs from that of heat transfer in that there exists a net mass flux due to evaporation or condensation. The physical system, as shown in Fig. 1, is characterized by three parameters; a length parameter, a gradient parameter, and a mass transfer parameter. The length parameter Kn and the gradient parameter M r , are defined as follows. Kn = x/ ((2kTl)/m )
beo " d x
(1)
in which/'1 = cold wall temperature, beq. = average local equilibrium collision frequency, and d = distance between the two plates. Mr = [(7'2 - 7"1 )/7"1 ]Kn
(2)
in which 7'2 = hot wall temperature. The parameter Kn may be interpreted as the reciprocal
Numerical solutions of non-linear kinetic equations
369
:]
r~-(d r~*Cd
--Tz
// i
T1- -
I i
x:O "~% Knudsen Leyer Y x
d '
Fig. 1. Schematic representation of the one-dimensional evaporation-condensation problem.
of the number of mean free paths between the two plates and the parameter Mr as the gradient in each mean free path. The mass transfer parameter/3 is defined as follows: rh-(d) _ molecular mass flux toward the cold plate at x = d /3 = n~+(d) molecular mass flux toward the hot plate at x = d
(3)
The mass fluxes n~-(x) and rh+(x) are defined as follows: rh-(x) -- f f f v~f- dr)
(4)
rh+(x) = f~f v x f + d~
(5)
in which vx = velocity in x-direction and f = f ( ~ , x) = velocity distribution function. The velocity distribution function f - a n d f + are for molecules moving toward the cold plate and for those toward the hot plate respectively. It should be pointed out that the mass flux ratio is not constant between two plates; therefore, the solutions obtained for a given value of fl maintained at the hot plate (x = d) would be different from those obtained for the same value of fl but defined at the cold plate (x = 0). The selection of the definition of the mass transfer parameter may depend on the nature of the problem. For example, Cheremisin[9] used the definition of rh+(0)/m-(0) for his study of mass transfer problems with high absorption at the cold plate. The Boltzmann equation for one-dimensional flows may be written as
vx(df/dx) = (a - bf)
(6)
in which f = f ( ~ , x) = velocity distribution function and ( a - bf) = collision integral. a(~, x), the gain term in the collision integral, is the rate of scattering of molecules in d~ of velocity space, b(f~, x)f(~, x), the loss term, is the rate of scattering of molecules out of d~ of velocity space. Our objective is to solve the Boltzmann equation for a problem by using Nordsieck's method with the following boundary conditions: (1) At each boundary the distribution function of the reflected molecules is half-range
370
S.M. YEN
Maxwellian, representing either complete or incomplete energy accommodation and isootropic emission. (Other gas surface interaction models such as the Maxwell model have been incorporated in our computer program and could also be used.) The accommodation coefficient a is defined as follows: = (T~ -
T,)/(T~
-
(7)
Tw)
in which T~ and 7', are the temperatures of the incident and reflected molecules respectively, Tw = plate temperature. (2) The ratio of molecular fluxes at x = d, t h - ( d ) / t h +(d), is ft. As indicated above, we have also obtained the solution of the Krook equation. The Krook equation, which replaces the collision integral by a statistical (or relaxation) model is vx(df/dx ) = b(n~ - f)
(8)
where b = collision frequency nO = Maxwell-Boltzmann distribution function with the same density, average velocity, and temperature as those determined from the distribution function f. 3. NUMERICAL METHOD The molecular flow in one-dimensional problems is axially symmetrical, that is, the velocity distribution functionfis not a function of the azimuthal angle. The velocity space considered is therefore two-dimensional; however, the individual collision is still treated as a threedimensional phenomenon. As shown in Fig. 2, the velocity space is quantized by covering the semi-circular region with 226 fixed cells.
i
i
i
i
-z3 -19 -~s -11
i
-7 -3 6
"
I
I
~ ' 12 1'5 1'9'
~
Vxrn
Ftg. 2. Velocityspace. In our iterative scheme, we replace the position variable x in our iterative scheme to integrate the Boltzmann equation (and the Krook equation); therefore, the Boltzmann equation used is in the following form: = (a - b f )
(9)
(dz/dx) = K n ( x / 2 / n ) n T 1/2
(10)
vx(df/dz)(d~/dx)
where T is the integration variable and
in which n = local number density, T = local temperature, and K n = Knudsen number defined by equation (1).
Numerical solutions of non-linear kinetic equations
~t--A--
-A "1-
371
VI
vx Fig. 3. Set-up of velocity spaces for implementation of the numerical method. In order to explain our numerical method it is convenient to show schematically in Fig. 3 the velocity space at the z-position. The velocity space at each z-position is a semi-circle and has two quadrants. The quadrant with vx > 0 is for the molecules moving toward the hot plate and that with v~ < 0, for the molecules moving toward the cold plate. There are 113 velocity points in each quadrant. One Boltzmann equation is associated with each velocity point. The integration of 113 Boltzmann equations of the right-moving molecules is from z = 0 to z = 1, while, for the left-moving molecules, the integration of the 113 Boltzmann equations starts from z = 1. The form of the distribution function at each starting position is determined by the boundary conditions at each plate and is thus known; however, the density is unknown. The distribution function at the terminating position is also unknown. The integration scheme of the 226 Boltzmann equation is also divided into two parts for the two groups of molecules. The computation scheme was designed according to the following conditions peculiar to the problem: (1) The boundary conditions at z = 0 and z = 1 are half-range Maxwellians with 7'1 and 7'2 equal to either the wall temperatures or a fraction of them (for 0t < 1). The ratio rt2/n 1 is unknown, nl and n 2 are twice the number density of reflected molecules at the cold plate and the hot plate respectively. (2) There is a net molecular mass flux. The iterative procedure consists of the following steps: Step 1 Divide the variable z into J equal intervals. Step 2 Assume initial distribution function at each z-position. Step 3 Evaluate the collision integral for each of these distribution functions. Step 4 Integrate the resulting differential equation for each of the 226 velocity bins to get the first iterate o f f ( J , z). The integration of the 113 Boltzmann equations in the right quadrant is from 0 to 1, and in the left quadrant, from 1 to 0.
372
s.M. YEN
Step 5
The values of rh-(1) and rh+(1) are computed according to equations (4) and (5) and the values o f f - ( g , 1) (i.e., at the left quadrant which is for Vx < 0) are adjusted to insure that the net mass flux ratio is fl for the second iteration. Step 6 Repeat steps (3), (4) and (5) to perform successive iterations. We use the rms difference of successive iterates to measure the risidual 6 f a n d to monitor the convergence of the iterative process. We used the same iterative procedure to solve the Krook equation. Since the systematic errors due to integration for the solutions of the two equations are thus comparable, it would be meaningful to make comparative studies. Furthermore, we have obtained directly the distribution function which would be used to make more detailed comparisons with the corresponding solutions of the Boltzmann equation. 4. DISCUSSION OF RESULTS Solutions of the Boltzmann and Krook equations were obtained for the case of complete energy accommodation (~ = 1) and T1/T 2 = 0"5 and for the following ranges of parameters: (1) Kn = 0.1 to 100, (2)/~ = 1.1-1.9. Of particular interest in this problem is the evaporation (or condensation) rate with pressure and temperature differences as the driving forces. The expression of the evaporation rate obtained by Patton and Springer[8] and Shankar[5] for the linear case is
r&.=~ L2--0-~n+ ~ - L~3o-KT79~
(11)
where rhe is the evaporation rate according to the linear prediction and has the unit of (~)P~o. The above expression shows explicitly the effects of the temperature difference AT and the pressure difference AP. In order to show the effect of the departure from equilibrium through the gradient parameter M r, we rewrite equation (11) as
,h,=~ k3-~+
e~ 30K.+9
We thus observe that, according to the linear theory, the major driving force in the evaporation phenomenon is the pressure difference, as indicated by the first term on the right hand side of equation (11). The non-equilibrium effect to hinder evaporation is expressed explicitly through the gradient parameter M r and the length parameter Kn in the second term. Before we make comparison between our nonlinear calculations of the evaporation rate and the linear predictions, it should be pointed out that both equations (11) and (12) are valid only for the linear case when AT is small, including the free molecule case. As shown in Table 1, for T1/T 2 = 0.5, the free molecule evaporation rate according to the linear formula is less than the exact value for small values of fl and larger for larger values of ft. The range of the ratio of free molecule evaporation rate to the linear calculation varies between 0.8871 and 0.7349 for fl = 1'3 to 1.9. Our solutions were obtained for the temperature ratio (T~/T2) of 0.5; therefore, the length parameter Kn is equal to the gradient parameter M r (see equation (2)). The value AP/Poo is not treated as a parameter but is obtained from the solution itself. We have studied in detail the microscopic and macroscopic gas properties from these solutions and have selected the following results for presentation in this paper: (1) Tabulation of the evaporation rate rh-(1), to be referred to simply as rh in the remainder of this paper, and AP/Po~ for both Boltzmann and Krook solutions for hard spheres
Numerical solutions of non-linear kinetic equations
373
and comparison of rh with the linear prediction for the computed AP/Poo and the parameters Kn and M r . (2) Comparison of Boltzmann evaporation rates rh for hard spheres with those for Maxwellian molecules. (3) Comparison of the Boltzmann values of the temperature, the density and the average velocity distributions with the corresponding K r o o k results. (4) Graphical display of typical velocity distribution functions. We shall briefly describe these results. Table 1. Ratio of free molecule evaporation rate for T~/T2 = 0.5 Mass flow parameter, fl
1"1
1"3
1"5
1-7
1'9
(rhrM)/OhFu)e
1"805
0"8871
0"8051
0"7745
0"7349
rhrM = free molecule evaporation rate (theM)e = free molecule evaporation rate according to linearized formula In Table 2, we give the Boltzmann and K r o o k results on the evaporation rate rh together with the parameters Kn and M r and the calculated values of AP/P~. Using equation (12), we compute the evaporation rate the according to linear prediction. In order to make comparative study with the linear results we tabulate both rhe and the ratio rh/rh e. In addition, we give the ratio of the evaporation rate with the respective free molecule value. We note in Table 2 that, for about the same evaporation rate rh, the K r o o k solutions give lower values of the pressure difference AP/P~ than the Boltzmann solutions. This difference implies that, for the same value of pressure difference, the K r o o k evaporation rates are lower. Figure 4 shows the variation of rh/rhrM for M r = I and Kn = 1 for both Boltzmann and
o
•E ~
I
I.C
•E ~
I
I
MT=I, Kn:l • Boltzmann o Krook
n.'~ o o
o
g ~ os~
g~
o o
,.6~_
•
._o:.g ~
o
o
0.9 0
I 0.5
I 1.0
I 1.5
A_~_p P~
Fig. 4. Variation of the evaporation rate ratio, th/rh~M, with AP]P~ofor Mr = 1 and Kn = 1. K r o o k results. The variation of the difference with the pressure ratio was found to depend on the non-equilibrium parameters, M r and Kn. As discussed above, the linear prediction is in error even in the free molecule case; therefore, it is to be expected that, for larger values of Kn (e.g., 10) the ratio rh/m t should be close to (rhrM)/(dlFM)¢ as given in Table 1. The Boltzmann results shown in Table 2 indeed bear out this prediction. As Kn and M r increase, this ratio decreases, indicating that the linear formula underestimates the non-equilibrium hindrance effect.
/ a. j
374
S. M. YEN Table 2. Evaporation rate ih 7'1/7"2 = 0.5, hard spheres ( 1 ) / ~ = 1.1
Boltzmann Kn Mr 3 3 1 1 0"5 0"5 0"1 0"1
APIP~ 0"5209 0"4617 0"4107 0.2494
rh 0"01556 0'01449 0'01419 0"01144
rht 0"00972 0'01039 0"01291 0"01649
Krook 3 1 0"5 0.1
3 1 0"5 0"1
0"4843 0"4020 0"3356 0-1963
0-01525 0"01440 0"01939 0"01232
0"00395 0"00113 0"00146 0-00875
(2) fl = 1-3 Boltzrnann l0 10 3 3 1 1 0"5 0"5
0'8363 0"8084 0.7508 0"7990
0'04734 0-04598 0"04363 0"04173
0.05270 0"04822 0"04342 0'04175
10 3 1 0"5
0"8192 0"7653 0"6747 0"5998
0-04778 0"04703 0"04550 0'04410
(3) fl = 1.5 Boltzmann 10 10 3 3 1 1 0-5 0"5 0"l 0"1
1 "126 1.100 1.040 0"9773 0"7574
rh#hru 0"9778 0'9101 0'8913 0'7188
/hJ(theu) t 1"098 1'173 1"458 1"862
0"9580 0.9047 0"8630 0"7744
0"4465 0"1279 0"1646 0.9879
0'8984 0-9534 1.005 1"000
0"9915 0.9630 0.9138 0'8739
0'9782 0"8952 0"8060 0"7750
0"05541 0"05502 0"05522 0"05687
0"8623 0"8546 0"8240 0"7755
1.000 0.9850 0.9530 0"9237
1"020 1"021 1"025 1"056
0"07992 0"07890 0"07591 0"07384 0"06464
0"1014 0"1010 0-1000 0"9003 0-09060
0-7879 0"7813 0.7590 0"7434 0"7135
1.000 0"9915 0.9539 0'9279 0"8123
1"026 1"021 1.011 1"000 0"9162
1"105 1-05l 0'9562 0"8754 0-7074
0.7913 0.07714 0.07360 0"07063 0-06458
0.09805 0"90325 0"08706 0"08378 0'08330
0-8070 0-8273 0"8453 0"8431 0.7753
0"9944 0'9694 0'9348 0.8876 0.8116
0"9915 0"9430 0"8804 0.8472 0.8424
(4) fl = 1-7 Boltzmann 3 3 1 l
1.396 1.344
0.1112 0.1089
0.1476 0.1472
0.7534 0-7388
0.9983 0.9763
1.026 1.023
Krook 3 l
3 l
1-342 1"247
0"1088 0.1044
0"1391 . 0'1322
0"7823 0"7901
0"9865 0.9373
0"9830 0-9468
(5) fl = 1.9 Boltzmann 3 3 1 1
1.697 1.641
0.1440 0"1401
0'1950 0'1932
0.7385 0-7249
1.000 0.9780
1.032 1"023
Krook 3 1
1.638 1.548
0"1410 0.1362
0.1857 0.1789
0'7590 0"7615
0.9840 0"9508
0"9830 0.9468
Krook 10 3 1 0"5
Krook 10 3 1 0-5 0"l
l0 3 l 0"5 0.1
3 1
rh#h c 1'601 1-394 1"099 0"6937 3.856 12.914 9"426 1"409
T a b l e 3 s h o w s t h e r e s u l t s o n e v a p o r a t i o n r a t e s f o r M a x w e l l i a n m o l e c u l e s f o r fl = 1.5. I n c o m p a r i n g w i t h t h e c o r r e s p o n d i n g r e s u l t s f o r h a r d s p h e r e s , we n o t e t h a t t h e d i f f e r e n c e is s m a l l f o r K n -- 10 a n d 3 a n d t h a t t h e e v a p o r a t i o n r a t e s f o r M a x w e l l i a n m o l e c u l e s a r e h i g h e r f o r K n = 1 a n d 0"5.
Numerical solutions of non-linear kinetic equations Table 3. Evaporation rate th fl = 1.5, Kn
MT
AP~
10 3 1 0"5
10 3 1 0"5
1 "119 1"~5 1"039 0"9934
T1/T2 =
375
0.5, Maxwellian molecules
m
mt
m/mr
m/mrM
me/(m~)e
0"07969 0"07912 0"07727 0"07526
0"1~2 0"1~2 0"~999 0"1018
0"7950 0"7896 0"7738 0"7395
I'001 0"9942 0"9710 0"9458
1"013 1"013 1"~1 1"029
I
l
l
l
l
l
l
I
l
0.7O
,B :1.5 Kn:l ]-1/1"2=05
o &d 0
0.6 o 0 =
=oa
I - I--
0.5
Q 0
0.4
• Boltzmonn (Max, prob. error < 0.21%) o Krook
• •o LO
f
o.31 o
]
[
I
I
0.2
I
0.4
i
I
I
0.6
I
0,8
10
Distance variable, T Fig. 5. V a r i a t i o n o f temperature w i t h the distance variable -r f o r f l = 1-5 and K n = ! . 1
1
l
I
1.05~~Q
I
1
I
I
I
/3:15 Kn:1
rz/T2:OS o
1.00
• o
~
• o
0.95-
o o •
090
Boltzmonn (Max.
prob. error < 0.26%) o Krook
0.85
o
I
l
02
I
I
[
I
I
°ot
t
04 06 0.8 Distance variable, -c
1.0
Fig. 6. Variation of density with the distance variable ~- for fl = 1"5 and
K n = 1.
376
S. M. YEN -042(
~bo I o
~×-043 ~
/3= 1.5, Kn= 1, T I / T 2 = 0.5 J [ I I I I I I I • Boltzmcmn (Mox | o prob error < 0 . 2 3 % ) - ~ •
o•
°oKrOok
~-044 /
< -045 t i
-046 037
I
036
¢, 0.35
I
i
•
+
]
i
I
I
L
!
I
[
I
8oltzmann ( M a x prob error< 0.].4%
•
•J
o o
0.33
o
1£ ~-I I•
032
O~_
O
o <~ o34
I
O O"
o Krook
I
O
•
• II 02
[I
:i It J [ 0.4 0.6 Distance Variable, T
i 0.8
I 10
Fig. 7. Variation of average velocities vx + (for molecules moving toward plate with /'2) and v~- (for molecules moving toward plate with T1) with distance variable ~-. fl = 1.5, K n = 1.
(/~=1.5,
Kn=3,
T1/T2=0.5,
"r=7/8}
Boltzmonn
Krook
(Vertical scale increased by I 2 . 5 % ) Fig. 8. Computer graphical display of distribution functions at ~" = 7/8. ~ = 1.5.
Kn =
3.
Numerical solutions of non-linear kinetic equations
377
Figures 5 and 6 give the variation of the temperature and density distributions as a function of ~. The differences between the Boltzmann and Krook results are similar to those found for the heat transfer problems, i.e., with/3 = 1. [3] Figure 7 shows the average velocity vx+ (for molecules moving toward the plate with temperature 7"2) and vx- (for molecules moving toward the plate with/'1). As the molecular motion advances, the difference between the Boltzmann and Krook values of average velocities increases. Figure 8 shows the computer graphical display of the velocity distribution function at = 7/8 for/3 = 1.5 and Kn = 3. The Krook result is qualitatively similar to the Boltzmann result; however, the quantitative difference is significant. (The vertical scale for the Krook distribution function is increased by 12-5 per cent.) 5. SUMMARY Our nonlinear Boltzmann calculations of the evaporation (condensation) rate indicate that the linear prediction underestimates the non-equilibrium hindrance effect. The Krook evaporation rates are lower than the corresponding Boltzmann results and, thus, deviate more from the linear values. The Boltzmann evaporation rates for the hard spheres are lower than those for the Maxwellian molecules at lower values of Knudsen number. The microscopic and macroscopic properties obtained from the Krook solutions differ appreciably from the corresponding Boltzmann results. REFERENCES Nordsieek A. and Hicks B. L., Proc. 5th Internat. Sym. on Rarefied Gas Dynamics, Vol. 1, p. 675, (1967). Hicks B. L., Yen S. M., and Reilly B. J.,J. Fluid Mech., 53, 85-111, (1972). Yen S. M., Int. J. Heat Mass Transfer, 15, 1865-1869, (1971). Sampson R. E. and Springer G. S., J. Fluid Mech., 36, 577. Shankar P. N., Ph.D. Thesis, Cal. Inst. of Tech., 1968. Edwards R. H. and Collins R. L., Proc. 6th lnternat. Syrup. on Rarefied Gas Dynamics, Vol. 2, p. 1489, 1969. 7. Pao Y. P. Phys. Fluids, 14, 1340 (1971). 8. Patton A. J. and Springer G. S., Proc. Sixth Int. Syrup. on Rarefied Gas Dynamics, p. 1497, 1969. 9. Cheremisin F. G., Izvestiya Akad. Nauk SSSR, Makhamika Zhidkosti i gaze (Bulletin, Acad. of Sciences of the USSR, Mechanics of Liquids and Gases), Vol. 2, p. 176, 1972.
1. 2. 3. 4. 5. 6.
NUMERICAL SOLUTIONS OF NON-LINEAR KINETIC EQUATIONS FOR A ONE-DIMENSIONAL EVAPORATION-CONDENSATION PROBLEM* S. M. YEN Coordinated Science Laboratory University of Illinois at Urbana-Champaign Urbana, Illinois 61801, U.S.A. (Received 11 December 1972)
Abstract--Numerical solutions have been obtained for both the nonlinear Boltz~ann equation (for two collision laws) and the Krook equation for a one-dimensional evaporationcondensation problem for a range of parameters. Our calculations of the evaporation (condensation) rate indicate that the linear prediction underestimates the non-equilibrium hindrance effect due to intermolecular collisions. The Krook evaporation rates are lower than the corresponding Boltzmann results and, thus, differ even more from the linear values. The Boltzmann evaporation rates for a gas of hard spheres are lower than those for Maxweilian molecules at lower values of Knudsen number. The micorscopic and macroscopic properties obtained from the Krook solutions differ appreciably from the corresponding Boltzmann results. 1. I N T R O D U C T I O N The transport phenomena occur in a gas because of the departure from thermal equilibrium. When the departure is small, the transports may be expressed in terms of the transport coefficients and the property gradients. For a given problem, the property gradients are obtained by solving the basic continuum gas dynamic equations with known relation between the transport coefficients and a property such as the temperature. Since the physics of the transport phenomena is considered implicitly through the transport coefficients, we face only one problem, namely, a mathematical one of solving the continuum gas dynamic equations. Under conditions far from thermal equilibrium, the concept of transport coefficients is not valid and the continuum equation is no longer applicable. The gas under such conditions can be described by the Boltzmann equation, the basic equation in kinetic theory. This equation contains the term collision integral which represents the rate of change of velocity distribution function (or molecular population density) due to intermolecular collisions. The collision integral which depends on the local departure from thermal equilibrium is essentially the microscopic transport property. The solution of the Boltzmann equation yields the velocity distribution function from which any transport property can be evaluated. In non-equilibrium gas dynamics, we, therefore, face two problems: (1) physics: to find the effect of intermolecular collisions on the microscopic transport, and (2) mathematics: to solve the nonlinear Boltzmann equation. * This work was supported in part by the Joint Services Electronics Program (U.S. Army, U.S. Navy and U.S. Air Force) under Contract DAAB-07-67-C-0199, in part by the National Oceanic and Atmospheric Administration under Grant N22-1-72-(G), and in part by the Office of Naval Research under Contract N00014-67-A-305-0001. 367
368
S.M. YEN
The advent of the high speed electronic digital computers has spurred the development of methods of solving equations of kinetic theory, including the Boltzmann equation. The major difficulty of solving the Boltzmann equation is the intractable form of the collision integral in it. Straightforward numerical quadrature would require days on the fastest computers today. The development of Nordsieck's Monte Carlo method[I] has made possible the direct solution of the Boltzmann equation for several non-equilibrium gas flow problems, including those of heat transfer[2,3]. His method consists of two distinct steps: (1) to evaluate the Boltzmann collision integral by using a statistical sampling technique closely resembling the real statistical collision phenomena in the gas and (2) to solve the Boltzmann difference equation by using an iterative scheme. The linearized kinetic theory problem of vapor condensation on or evaporation from spherical liquid droplets under non-equilibrium conditions has been tackled by Springer and Sampson[4], Shankar[5], and Edwards and Collins[6] using the moment method and by Pao[7] using the Krook equation. Springer, et al. used Lees' Ansatz; Edwards, et al., Grad's expression. These efforts have also led to the study of a simpler one-dimensional evaporationcondensation problem which, in essence, is the heat transfer problem of rarefied gases with net mass flux between two plates at different temperatures. The significance of studying this conceptually simpler problem is twofold. First, it is possible to solve this simple problem more accurately using basic equations of kinetic theory under a large range of nonequilibrium conditions. Secondly, the solutions for this problem would be relevant to the problem of more complex geometry and the study of their behavior would lead to better understanding of the non-equilibrium gas behavior for evaporation-condensation problems in general. The moment method with Lees' Ansatz has been used to solve this problem by Springer and Patton[8] in the linearized case and Edwards and Collins[6] in the non-linear case. Pao[7] solved the linearized Krook equation. We have solved both the non-linear Boltzmann equation (for hard spheres and Maxwellian molecules) and the non-linear Krook equations for the one-dimensional evaporationcondensation problem under a wide range of conditions. Our solutions yield the velocity distribution function, the collision integral for the Boltzmann solution as well as their moments of physical interest. 2. FORMULATION OF THE PROBLEM We consider the one-dimensional problem of rarefied gases between two parallel plates having two different temperatures, T1 and T2, and at distance d apart. This problem differs from that of heat transfer in that there exists a net mass flux due to evaporation or condensation. The physical system, as shown in Fig. 1, is characterized by three parameters; a length parameter, a gradient parameter, and a mass transfer parameter. The length parameter Kn and the gradient parameter M r , are defined as follows. Kn = x/ ((2kTl)/m )
beo " d x
(1)
in which/'1 = cold wall temperature, beq. = average local equilibrium collision frequency, and d = distance between the two plates. Mr = [(7'2 - 7"1 )/7"1 ]Kn
(2)
in which 7'2 = hot wall temperature. The parameter Kn may be interpreted as the reciprocal
Numerical solutions of non-linear kinetic equations
369
:]
r~-(d r~*Cd
--Tz
// i
T1- -
I i
x:O "~% Knudsen Leyer Y x
d '
Fig. 1. Schematic representation of the one-dimensional evaporation-condensation problem.
of the number of mean free paths between the two plates and the parameter Mr as the gradient in each mean free path. The mass transfer parameter/3 is defined as follows: rh-(d) _ molecular mass flux toward the cold plate at x = d /3 = n~+(d) molecular mass flux toward the hot plate at x = d
(3)
The mass fluxes n~-(x) and rh+(x) are defined as follows: rh-(x) -- f f f v~f- dr)
(4)
rh+(x) = f~f v x f + d~
(5)
in which vx = velocity in x-direction and f = f ( ~ , x) = velocity distribution function. The velocity distribution function f - a n d f + are for molecules moving toward the cold plate and for those toward the hot plate respectively. It should be pointed out that the mass flux ratio is not constant between two plates; therefore, the solutions obtained for a given value of fl maintained at the hot plate (x = d) would be different from those obtained for the same value of fl but defined at the cold plate (x = 0). The selection of the definition of the mass transfer parameter may depend on the nature of the problem. For example, Cheremisin[9] used the definition of rh+(0)/m-(0) for his study of mass transfer problems with high absorption at the cold plate. The Boltzmann equation for one-dimensional flows may be written as
vx(df/dx) = (a - bf)
(6)
in which f = f ( ~ , x) = velocity distribution function and ( a - bf) = collision integral. a(~, x), the gain term in the collision integral, is the rate of scattering of molecules in d~ of velocity space, b(f~, x)f(~, x), the loss term, is the rate of scattering of molecules out of d~ of velocity space. Our objective is to solve the Boltzmann equation for a problem by using Nordsieck's method with the following boundary conditions: (1) At each boundary the distribution function of the reflected molecules is half-range
370
S.M. YEN
Maxwellian, representing either complete or incomplete energy accommodation and isootropic emission. (Other gas surface interaction models such as the Maxwell model have been incorporated in our computer program and could also be used.) The accommodation coefficient a is defined as follows: = (T~ -
T,)/(T~
-
(7)
Tw)
in which T~ and 7', are the temperatures of the incident and reflected molecules respectively, Tw = plate temperature. (2) The ratio of molecular fluxes at x = d, t h - ( d ) / t h +(d), is ft. As indicated above, we have also obtained the solution of the Krook equation. The Krook equation, which replaces the collision integral by a statistical (or relaxation) model is vx(df/dx ) = b(n~ - f)
(8)
where b = collision frequency nO = Maxwell-Boltzmann distribution function with the same density, average velocity, and temperature as those determined from the distribution function f. 3. NUMERICAL METHOD The molecular flow in one-dimensional problems is axially symmetrical, that is, the velocity distribution functionfis not a function of the azimuthal angle. The velocity space considered is therefore two-dimensional; however, the individual collision is still treated as a threedimensional phenomenon. As shown in Fig. 2, the velocity space is quantized by covering the semi-circular region with 226 fixed cells.
i
i
i
i
-z3 -19 -~s -11
i
-7 -3 6
"
I
I
~ ' 12 1'5 1'9'
~
Vxrn
Ftg. 2. Velocityspace. In our iterative scheme, we replace the position variable x in our iterative scheme to integrate the Boltzmann equation (and the Krook equation); therefore, the Boltzmann equation used is in the following form: = (a - b f )
(9)
(dz/dx) = K n ( x / 2 / n ) n T 1/2
(10)
vx(df/dz)(d~/dx)
where T is the integration variable and
in which n = local number density, T = local temperature, and K n = Knudsen number defined by equation (1).
Numerical solutions of non-linear kinetic equations
~t--A--
-A "1-
371
VI
vx Fig. 3. Set-up of velocity spaces for implementation of the numerical method. In order to explain our numerical method it is convenient to show schematically in Fig. 3 the velocity space at the z-position. The velocity space at each z-position is a semi-circle and has two quadrants. The quadrant with vx > 0 is for the molecules moving toward the hot plate and that with v~ < 0, for the molecules moving toward the cold plate. There are 113 velocity points in each quadrant. One Boltzmann equation is associated with each velocity point. The integration of 113 Boltzmann equations of the right-moving molecules is from z = 0 to z = 1, while, for the left-moving molecules, the integration of the 113 Boltzmann equations starts from z = 1. The form of the distribution function at each starting position is determined by the boundary conditions at each plate and is thus known; however, the density is unknown. The distribution function at the terminating position is also unknown. The integration scheme of the 226 Boltzmann equation is also divided into two parts for the two groups of molecules. The computation scheme was designed according to the following conditions peculiar to the problem: (1) The boundary conditions at z = 0 and z = 1 are half-range Maxwellians with 7'1 and 7'2 equal to either the wall temperatures or a fraction of them (for 0t < 1). The ratio rt2/n 1 is unknown, nl and n 2 are twice the number density of reflected molecules at the cold plate and the hot plate respectively. (2) There is a net molecular mass flux. The iterative procedure consists of the following steps: Step 1 Divide the variable z into J equal intervals. Step 2 Assume initial distribution function at each z-position. Step 3 Evaluate the collision integral for each of these distribution functions. Step 4 Integrate the resulting differential equation for each of the 226 velocity bins to get the first iterate o f f ( J , z). The integration of the 113 Boltzmann equations in the right quadrant is from 0 to 1, and in the left quadrant, from 1 to 0.
372
s.M. YEN
Step 5
The values of rh-(1) and rh+(1) are computed according to equations (4) and (5) and the values o f f - ( g , 1) (i.e., at the left quadrant which is for Vx < 0) are adjusted to insure that the net mass flux ratio is fl for the second iteration. Step 6 Repeat steps (3), (4) and (5) to perform successive iterations. We use the rms difference of successive iterates to measure the risidual 6 f a n d to monitor the convergence of the iterative process. We used the same iterative procedure to solve the Krook equation. Since the systematic errors due to integration for the solutions of the two equations are thus comparable, it would be meaningful to make comparative studies. Furthermore, we have obtained directly the distribution function which would be used to make more detailed comparisons with the corresponding solutions of the Boltzmann equation. 4. DISCUSSION OF RESULTS Solutions of the Boltzmann and Krook equations were obtained for the case of complete energy accommodation (~ = 1) and T1/T 2 = 0"5 and for the following ranges of parameters: (1) Kn = 0.1 to 100, (2)/~ = 1.1-1.9. Of particular interest in this problem is the evaporation (or condensation) rate with pressure and temperature differences as the driving forces. The expression of the evaporation rate obtained by Patton and Springer[8] and Shankar[5] for the linear case is
r&.=~ L2--0-~n+ ~ - L~3o-KT79~
(11)
where rhe is the evaporation rate according to the linear prediction and has the unit of (~)P~o. The above expression shows explicitly the effects of the temperature difference AT and the pressure difference AP. In order to show the effect of the departure from equilibrium through the gradient parameter M r, we rewrite equation (11) as
,h,=~ k3-~+
e~ 30K.+9
We thus observe that, according to the linear theory, the major driving force in the evaporation phenomenon is the pressure difference, as indicated by the first term on the right hand side of equation (11). The non-equilibrium effect to hinder evaporation is expressed explicitly through the gradient parameter M r and the length parameter Kn in the second term. Before we make comparison between our nonlinear calculations of the evaporation rate and the linear predictions, it should be pointed out that both equations (11) and (12) are valid only for the linear case when AT is small, including the free molecule case. As shown in Table 1, for T1/T 2 = 0.5, the free molecule evaporation rate according to the linear formula is less than the exact value for small values of fl and larger for larger values of ft. The range of the ratio of free molecule evaporation rate to the linear calculation varies between 0.8871 and 0.7349 for fl = 1'3 to 1.9. Our solutions were obtained for the temperature ratio (T~/T2) of 0.5; therefore, the length parameter Kn is equal to the gradient parameter M r (see equation (2)). The value AP/Poo is not treated as a parameter but is obtained from the solution itself. We have studied in detail the microscopic and macroscopic gas properties from these solutions and have selected the following results for presentation in this paper: (1) Tabulation of the evaporation rate rh-(1), to be referred to simply as rh in the remainder of this paper, and AP/Po~ for both Boltzmann and Krook solutions for hard spheres
Numerical solutions of non-linear kinetic equations
373
and comparison of rh with the linear prediction for the computed AP/Poo and the parameters Kn and M r . (2) Comparison of Boltzmann evaporation rates rh for hard spheres with those for Maxwellian molecules. (3) Comparison of the Boltzmann values of the temperature, the density and the average velocity distributions with the corresponding K r o o k results. (4) Graphical display of typical velocity distribution functions. We shall briefly describe these results. Table 1. Ratio of free molecule evaporation rate for T~/T2 = 0.5 Mass flow parameter, fl
1"1
1"3
1"5
1-7
1'9
(rhrM)/OhFu)e
1"805
0"8871
0"8051
0"7745
0"7349
rhrM = free molecule evaporation rate (theM)e = free molecule evaporation rate according to linearized formula In Table 2, we give the Boltzmann and K r o o k results on the evaporation rate rh together with the parameters Kn and M r and the calculated values of AP/P~. Using equation (12), we compute the evaporation rate the according to linear prediction. In order to make comparative study with the linear results we tabulate both rhe and the ratio rh/rh e. In addition, we give the ratio of the evaporation rate with the respective free molecule value. We note in Table 2 that, for about the same evaporation rate rh, the K r o o k solutions give lower values of the pressure difference AP/P~ than the Boltzmann solutions. This difference implies that, for the same value of pressure difference, the K r o o k evaporation rates are lower. Figure 4 shows the variation of rh/rhrM for M r = I and Kn = 1 for both Boltzmann and
o
•E ~
I
I.C
•E ~
I
I
MT=I, Kn:l • Boltzmann o Krook
n.'~ o o
o
g ~ os~
g~
o o
,.6~_
•
._o:.g ~
o
o
0.9 0
I 0.5
I 1.0
I 1.5
A_~_p P~
Fig. 4. Variation of the evaporation rate ratio, th/rh~M, with AP]P~ofor Mr = 1 and Kn = 1. K r o o k results. The variation of the difference with the pressure ratio was found to depend on the non-equilibrium parameters, M r and Kn. As discussed above, the linear prediction is in error even in the free molecule case; therefore, it is to be expected that, for larger values of Kn (e.g., 10) the ratio rh/m t should be close to (rhrM)/(dlFM)¢ as given in Table 1. The Boltzmann results shown in Table 2 indeed bear out this prediction. As Kn and M r increase, this ratio decreases, indicating that the linear formula underestimates the non-equilibrium hindrance effect.
/ a. j
374
S. M. YEN Table 2. Evaporation rate ih 7'1/7"2 = 0.5, hard spheres ( 1 ) / ~ = 1.1
Boltzmann Kn Mr 3 3 1 1 0"5 0"5 0"1 0"1
APIP~ 0"5209 0"4617 0"4107 0.2494
rh 0"01556 0'01449 0'01419 0"01144
rht 0"00972 0'01039 0"01291 0"01649
Krook 3 1 0"5 0.1
3 1 0"5 0"1
0"4843 0"4020 0"3356 0-1963
0-01525 0"01440 0"01939 0"01232
0"00395 0"00113 0"00146 0-00875
(2) fl = 1-3 Boltzrnann l0 10 3 3 1 1 0"5 0"5
0'8363 0"8084 0.7508 0"7990
0'04734 0-04598 0"04363 0"04173
0.05270 0"04822 0"04342 0'04175
10 3 1 0"5
0"8192 0"7653 0"6747 0"5998
0-04778 0"04703 0"04550 0'04410
(3) fl = 1.5 Boltzmann 10 10 3 3 1 1 0-5 0"5 0"l 0"1
1 "126 1.100 1.040 0"9773 0"7574
rh#hru 0"9778 0'9101 0'8913 0'7188
/hJ(theu) t 1"098 1'173 1"458 1"862
0"9580 0.9047 0"8630 0"7744
0"4465 0"1279 0"1646 0.9879
0'8984 0-9534 1.005 1"000
0"9915 0.9630 0.9138 0'8739
0'9782 0"8952 0"8060 0"7750
0"05541 0"05502 0"05522 0"05687
0"8623 0"8546 0"8240 0"7755
1.000 0.9850 0.9530 0"9237
1"020 1"021 1"025 1"056
0"07992 0"07890 0"07591 0"07384 0"06464
0"1014 0"1010 0-1000 0"9003 0-09060
0-7879 0"7813 0.7590 0"7434 0"7135
1.000 0"9915 0.9539 0'9279 0"8123
1"026 1"021 1.011 1"000 0"9162
1"105 1-05l 0'9562 0"8754 0-7074
0.7913 0.07714 0.07360 0"07063 0-06458
0.09805 0"90325 0"08706 0"08378 0'08330
0-8070 0-8273 0"8453 0"8431 0.7753
0"9944 0'9694 0'9348 0.8876 0.8116
0"9915 0"9430 0"8804 0.8472 0.8424
(4) fl = 1-7 Boltzmann 3 3 1 l
1.396 1.344
0.1112 0.1089
0.1476 0.1472
0.7534 0-7388
0.9983 0.9763
1.026 1.023
Krook 3 l
3 l
1-342 1"247
0"1088 0.1044
0"1391 . 0'1322
0"7823 0"7901
0"9865 0.9373
0"9830 0-9468
(5) fl = 1.9 Boltzmann 3 3 1 1
1.697 1.641
0.1440 0"1401
0'1950 0'1932
0.7385 0-7249
1.000 0.9780
1.032 1"023
Krook 3 1
1.638 1.548
0"1410 0.1362
0.1857 0.1789
0'7590 0"7615
0.9840 0"9508
0"9830 0.9468
Krook 10 3 1 0"5
Krook 10 3 1 0-5 0"l
l0 3 l 0"5 0.1
3 1
rh#h c 1'601 1-394 1"099 0"6937 3.856 12.914 9"426 1"409
T a b l e 3 s h o w s t h e r e s u l t s o n e v a p o r a t i o n r a t e s f o r M a x w e l l i a n m o l e c u l e s f o r fl = 1.5. I n c o m p a r i n g w i t h t h e c o r r e s p o n d i n g r e s u l t s f o r h a r d s p h e r e s , we n o t e t h a t t h e d i f f e r e n c e is s m a l l f o r K n -- 10 a n d 3 a n d t h a t t h e e v a p o r a t i o n r a t e s f o r M a x w e l l i a n m o l e c u l e s a r e h i g h e r f o r K n = 1 a n d 0"5.
Numerical solutions of non-linear kinetic equations Table 3. Evaporation rate th fl = 1.5, Kn
MT
AP~
10 3 1 0"5
10 3 1 0"5
1 "119 1"~5 1"039 0"9934
T1/T2 =
375
0.5, Maxwellian molecules
m
mt
m/mr
m/mrM
me/(m~)e
0"07969 0"07912 0"07727 0"07526
0"1~2 0"1~2 0"~999 0"1018
0"7950 0"7896 0"7738 0"7395
I'001 0"9942 0"9710 0"9458
1"013 1"013 1"~1 1"029
I
l
l
l
l
l
l
I
l
0.7O
,B :1.5 Kn:l ]-1/1"2=05
o &d 0
0.6 o 0 =
=oa
I - I--
0.5
Q 0
0.4
• Boltzmonn (Max, prob. error < 0.21%) o Krook
• •o LO
f
o.31 o
]
[
I
I
0.2
I
0.4
i
I
I
0.6
I
0,8
10
Distance variable, T Fig. 5. V a r i a t i o n o f temperature w i t h the distance variable -r f o r f l = 1-5 and K n = ! . 1
1
l
I
1.05~~Q
I
1
I
I
I
/3:15 Kn:1
rz/T2:OS o
1.00
• o
~
• o
0.95-
o o •
090
Boltzmonn (Max.
prob. error < 0.26%) o Krook
0.85
o
I
l
02
I
I
[
I
I
°ot
t
04 06 0.8 Distance variable, -c
1.0
Fig. 6. Variation of density with the distance variable ~- for fl = 1"5 and
K n = 1.
376
S. M. YEN -042(
~bo I o
~×-043 ~
/3= 1.5, Kn= 1, T I / T 2 = 0.5 J [ I I I I I I I • Boltzmcmn (Mox | o prob error < 0 . 2 3 % ) - ~ •
o•
°oKrOok
~-044 /
< -045 t i
-046 037
I
036
¢, 0.35
I
i
•
+
]
i
I
I
L
!
I
[
I
8oltzmann ( M a x prob error< 0.].4%
•
•J
o o
0.33
o
1£ ~-I I•
032
O~_
O
o <~ o34
I
O O"
o Krook
I
O
•
• II 02
[I
:i It J [ 0.4 0.6 Distance Variable, T
i 0.8
I 10
Fig. 7. Variation of average velocities vx + (for molecules moving toward plate with /'2) and v~- (for molecules moving toward plate with T1) with distance variable ~-. fl = 1.5, K n = 1.
(/~=1.5,
Kn=3,
T1/T2=0.5,
"r=7/8}
Boltzmonn
Krook
(Vertical scale increased by I 2 . 5 % ) Fig. 8. Computer graphical display of distribution functions at ~" = 7/8. ~ = 1.5.
Kn =
3.
Numerical solutions of non-linear kinetic equations
377
Figures 5 and 6 give the variation of the temperature and density distributions as a function of ~. The differences between the Boltzmann and Krook results are similar to those found for the heat transfer problems, i.e., with/3 = 1. [3] Figure 7 shows the average velocity vx+ (for molecules moving toward the plate with temperature 7"2) and vx- (for molecules moving toward the plate with/'1). As the molecular motion advances, the difference between the Boltzmann and Krook values of average velocities increases. Figure 8 shows the computer graphical display of the velocity distribution function at = 7/8 for/3 = 1.5 and Kn = 3. The Krook result is qualitatively similar to the Boltzmann result; however, the quantitative difference is significant. (The vertical scale for the Krook distribution function is increased by 12-5 per cent.) 5. SUMMARY Our nonlinear Boltzmann calculations of the evaporation (condensation) rate indicate that the linear prediction underestimates the non-equilibrium hindrance effect. The Krook evaporation rates are lower than the corresponding Boltzmann results and, thus, deviate more from the linear values. The Boltzmann evaporation rates for the hard spheres are lower than those for the Maxwellian molecules at lower values of Knudsen number. The microscopic and macroscopic properties obtained from the Krook solutions differ appreciably from the corresponding Boltzmann results. REFERENCES Nordsieek A. and Hicks B. L., Proc. 5th Internat. Sym. on Rarefied Gas Dynamics, Vol. 1, p. 675, (1967). Hicks B. L., Yen S. M., and Reilly B. J.,J. Fluid Mech., 53, 85-111, (1972). Yen S. M., Int. J. Heat Mass Transfer, 15, 1865-1869, (1971). Sampson R. E. and Springer G. S., J. Fluid Mech., 36, 577. Shankar P. N., Ph.D. Thesis, Cal. Inst. of Tech., 1968. Edwards R. H. and Collins R. L., Proc. 6th lnternat. Syrup. on Rarefied Gas Dynamics, Vol. 2, p. 1489, 1969. 7. Pao Y. P. Phys. Fluids, 14, 1340 (1971). 8. Patton A. J. and Springer G. S., Proc. Sixth Int. Syrup. on Rarefied Gas Dynamics, p. 1497, 1969. 9. Cheremisin F. G., Izvestiya Akad. Nauk SSSR, Makhamika Zhidkosti i gaze (Bulletin, Acad. of Sciences of the USSR, Mechanics of Liquids and Gases), Vol. 2, p. 176, 1972.
1. 2. 3. 4. 5. 6.