PUBLICATIONS DE L’INSTITUT MATHÉMATIQUE Nouvelle série, tome 91(105) (2012), 83–93 DOI: 10.2298/PIM1205083S

THE EXACT ANALYTICAL SOLUTION FOR THE GAS LUBRICATED BEARING IN THE SLIP AND CONTINUUM FLOW REGIME

Nevena D. Stevanović and Vladan D. Djordjević

Abstract. The exact analytical solution for the compressible two-dimensional gas flow in the microbearing is presented. The general slip-corrected Reynolds lubrication equation is derived and it is shown that it possesses an exact analytical solution. It is obtained by a suitable transformation of the independent variable, and it provides the pressure distribution in the bearing as well as the mass flow rate through it. By neglecting the rarefaction effect, this solution is also applicable to the continuum gas flow in the bearing, which also does not exist in the open literature. The obtained analytical solution can be usefully applied for testing the other, experimental or numerical results.

1. Introduction

The existing technology enables production of micro devices, which have wide applications in everyday life and in scientific investigations. Gas flow is a part of the most micro-electro-mechanical systems (MEMS), such as micro pumps, micro turbines, sensors for pressure, velocity and temperature measurements, systems for electric circuits cooling, magnetic disk storages etc. (Gad-El-Hak 2002). In these small systems the ratio between the mean free path of the molecules and the characteristic length of the microchannel, which is defined as the Knudsen number (Kn), is not negligible even at atmospheric pressure. As Kn increases rarefaction effects become more important and when the Knudsen number value comes over 10^{−2} , the continuum approach breaks down. In the range 10^{−2} < Kn < 10^{−1} , known as the slip flow regime, gas flow still obeys the continuum i.e., Navier– Stokes equations, but now with slip and temperature jump boundary conditions at the walls of the flow boundaries. In the range 10^{−1} < Kn < 10 (transitional flow regime) the Navier–Stokes equations break down and more complex Burnett equations of which the accuracy is of the order O(Kn^{2}) together with boundary conditions of the same, second-order accuracy, are used or the individual particlebased direct simulation Monte Carlo (DSMC) approach is to be employed. Finally, for Kn > 10 the gas flow is considered as a free molecular flow amenable to the methods of kinetic theory of gases (Vinsenti and Kruger 1986).

2010 Mathematics Subject Classification: Primary 76P05.

Key words and phrases: microbearing, Reynolds lubrication equation, analytical solution, second order slip velocity.

Accurate theoretical models and analytical solutions of micro flows are very useful tools in the design and analyzes of MEMS operation. Analytical solutions enable quick qualitative analyzes of MEMS operation, as well as efficient design procedures, without a need for experimental work. Also, reliable analytical solutions are benchmark tests for numerical methods that are usually applied to problems with complex flow channel geometry.

In the MEMS devices gas flow mostly goes on in the slip regime. The slip boundary condition at the channel wall was first defined by Maxwell (1879) as

u_{slip} = u_{gas} − u_{wall} = ±\frac{2 − σ_{ν}}{σ_{ν}} λ\frac{∂u}{∂n} |_{wall} + \frac{3}{4} \frac{μ}{ρT} \frac{∂T}{∂s} |_{wall}

while Smoluchowski (1898) proposed temperature jump at the wall as

T_{jump} = T_{gas} − T_{wall} = ±\frac{2 − σ_{T}}{σ_{T}} ( \frac{2γ}{γ + 1} ) \frac{λ}{Pr} \frac{∂T}{∂n} |_{wall}

where u, T, ρ, μ, γ are respectively velocity, temperature, density, dynamic viscosity, and specific heat ratio, λ is the free path of the molecules, σ_{ν} and σ_{T} are tangential momentum and thermal accommodation coefficients, while Pr is the Prandtl number. The operators ∂/∂n and ∂/∂s denote the normal and tangential derivatives at the wall surface. The sign in front of the first term on the right-hand side in velocity and temperature boundary conditions depends on the orientation of the axis perpendicular to the wall. The sign is plus when the normal axis is orientated from the wall, while it is minus when the normal axis is orientated towards the wall.

Comparisons between results obtained by experimental studies and analytical studies which are based on the first-order Maxwell–Smoluchowski boundary conditions, have shown some discrepancies. This has encouraged attempts in the literature to modify the existing slip-boundary conditions and define higher accuracy second order boundary conditions. By them a higher accuracy in the slip flow regime is achieved, as well as extension of the slip solutions application to the part of the transitional regime. The generalized second order boundary condition for isothermal flow is:

u_{slip} = ±A_{1}λ\frac{∂u}{∂n} |_{wall} − A_{2}λ^{2} \frac{∂^{2}u}{∂n^{2}} |_{wall} ,

where A_{1} and A_{2} are the first and second order slip coefficients which are differently defined by several authors in the literature. Table 1 present some values of A_{1} and A_{2} proposed by different authors and it can be seen that consensus about second order slip coefficient has not been reached (Barber and Emerson 2006, Lockerby et al. 2004). Hence, the value for A_{2} must be tested and chosen carefully for each analyzed flow condition.

In MEMS devices two general gas flow problems exist. These are pressure driven and shear driven flows through a channel or a pipe. In the slow (low Mach number) the pressure driven gas flow in the continuum regime pressure distribution

Author:  A_{1} A_{2}
Schamberg (1947)  1 5π/12
Deissler (1964)  1 9/8
Hsia and Domoto (1983)  1  0.5
Maxwell (1879)  1  0
Beskok and Karniadakis (1996) 1 −0.5

Table 1. First and second order slip coefficients proposed in the literature

is linear, which means that the gas flow could be treated as incompressible. Differently, in the slip regime, despite the low Mach number flow condition, experiments show that the pressure distribution is nonlinear, so that gas flow must be treated as compressible. The consequence is that the pressure distribution derived from the basic flow equations (inertia is neglected) for the slip pressure driven gas flow in a channel or a pipe is governed by a nonlinear first order differential equation. An exact analytical solution of this equation for the second-order boundary condition was readily derived for both the pressure driven gas flow between parallel plates and the pipe flow (Karniadakis and Beskok 2002).

The simplest form of a shear driven flow is the flow between parallel infinite plates caused by moving one of them, i.e., the Couette flow. Since the pressure is uniform, for isothermal conditions the flow can be treated as incompressible, which simplifies the problem and enables obtaining analytical solution for the velocity field in the slip regime. In gas lubrication theory, which is widely used in MEMS technologies, the gas flow is also caused by one moving plate and is considered as a shear driven flow. But now, in order to produce a load capacity, one wall is inclined and of finite length, which leads to the remarkable pressure variation in the stream-wise direction and necessitates the treatment of this flow as compressible. The pressure distribution in a gas lubricated bearing is governed by the so-called Reynolds equation. Under certain conditions (Szeri 1998) it can be readily derived from the Navier–Stokes equations, for both no-slip and slip boundary conditions. As in the case of pressure driven flow, this equation is nonlinear, but until now, and as far as we are informed, an exact analytical solution of it has not been found.

We show in this paper that such a solution exists. It is found by suitably transforming the independent variable (microbearing channel height) in the slipcorrected Reynolds equation and by obtaining the solution by quadratures. The validity of the solution is proved by comparison with numerical results available in the literature.

2. The exact analytical solution of the slip-corrected Reynolds lubrication equation

In this paper the exact analytical solution for the microbearing gas flow is presented. The lubrication problem analysed in the paper is depicted in Fig. 1. The isothermal two-dimensional compressible slip gas flow is considered. The low Mach number flow condition in the bearing is assumed. Hence, inertia effect can

Figure 1. Microbearing geometry

be neglected. Moreover, the channel cross section is slowly varying, thus the crosswise velocity component being much smaller than the stream-wise component. So, within this well known approximations (Szeri 1998) made at the derivation of the Reynolds equation, the continuum and extremely simplified Navier–Stokes equations read:

Ṁ = ∫_{0}^{h} ρu dy = const (2.1)

\frac{∂^{2}u}{∂y^{2}} = \frac{1}{μ} \frac{dp}{dx}, (2.2)

where ρ is the variable density, p is the pressure, uis the stream-wise velocity component, μ = const is the dynamic viscosity, Ṁ is the mass flow rate per unit width, while the other denotations are clearly seen in Fig. 1.

The equation (2.2) for the slip flow regime should be solved with the mentioned second order boundary conditions:

(2.3)  \begin{aligned}y = 0 : u = u_{0}(x) & = u_{w} + A_{1}λ\frac{∂u}{∂y} − A_{2}λ^{2} \frac{∂^{2}u}{∂y^{2}} \\ y = h : u = u_{1}(x) & = −A_{1}λ\frac{∂u}{∂y} − A_{2}λ^{2} \frac{∂^{2}u}{∂y^{2}}, \end{aligned}

where u_{0}(x) and u_{1}(x) are slip velocities at the bottom and top wall. The solution for the velocity field is easily found from momentum equation (2.2) and boundary conditions (2.3) as:

(2.4)  u = \frac{h^{2}}{2μ} \frac{dp}{dx} \frac{y^{2}}{h^{2}} ( u_{0} − u_{1} + \frac{h^{2}}{2μ} \frac{dp}{dx}) \frac{y}{h} + u_{0}

This solution is further used in continuity equation (2.1) that expresses the constancy of the mass flow rate Ṁ through the bearing to get a governing equation for the pressure distribution in the bearing. Then, when the independent variable x is replaced by h(x) (Fig. 1) and the equation of state for an ideal gas ρ = p/(RT) is utilized (R is the gas constant), it reads

(2.5)  − \frac{h^{3}}{12μRT} \frac{dh}{dx}p\frac{dp}{dh} + \frac{u_{0} + u_{1}}{2RT} hp = Ṁ

Equation (2.5) is transformed into the non-dimensional form in the following way: X = x/l, H = h/h_{e}, P = p/p_{e}, U_{0} = u_{0}/u_{w}, U_{1} = u_{1}/u_{w}, where l is the microbearing length, h_{e} the microbearing height at the exit, p_{e} pressure at the channel outlet and u_{w} is the velocity of the infinite plate positioned at y = 0. Then equation (2.5) becomes:

(2.6)  − \frac{1}{Λ} \frac{dH}{dX} H^{3}P \frac{dP}{dH} + (U_{0} + U_{1})HP = \frac{Ṁ}{Ṁ_{c}} ,

where Λ = 6μu_{w}l/p_{e}h_{e}^{2} is the bearing number and Ṁ_{c} = p_{e}h_{e}u_{w}/2RT is the mass flow rate in a Couette flow.

Before proceeding further, the sum of the slip velocities U_{0} + U_{1} was evaluated by using the general velocity field (2.4) for the case of the second order boundary conditions (2.3):

(2.7)  U_{0} + U_{1} = 1 − \frac{6}{Λ} \frac{dH}{dX} H^{2} \frac{dP}{dH} (A_{1} Kn+2A_{2} Kn^{2} ),

where Kn = λ/h is the local value of the Knudsen number. Since the mean free path of molecules for isothermal flow conditions is inversely proportional to the pressure, the local Knudsen number would be expressed as:

(2.8)  Kn = Kn_{e} /P H,

where Kn_{e} = λ_{e}/h_{e} is the the Knudsen number at the microbearing exit. Furthermore, for the linearly varying channel cross section defined as H = H_{i} − (H_{i} − 1)X, where H_{i} is non-dimensional parameter defined as ratio of the inlet and outlet microbearing height H_{i} = h_{i}/h_{e}, equation (2.6) finally attains the form of the general slip-corrected Reynolds lubrication equation:

(2.9)  \frac{H_{i} − 1}{Λ} H^{3}P \frac{dP}{dH} [1 + 6 Kn(A_{1} + 2A_{2} Kn)] + P H = \frac{Ṁ}{Ṁ_{c}}

Integration of this equation (analytical or numerical) by employing two boundary conditions (see (2.13)) yields as a result not only the already mentioned pressure distribution, but also the mass flow rate, which is not known beforehand.

The appropriate transformation of equation (2.9) by introducing new independent variable Z:

(2.10)  Z = \frac{1}{P (H)H}

enables us to get analytically the exact solution. The new independent variable Z has a simple physical meaning. It follows from (2.8) that Z = Kn/ Kn_{e}, i.e., it is the ratio between the local value of the Knudsen number and its exit value. The solution obtained in the similar way, but by introducing a dependent variable Z = \frac{1}{P H(P)} which has the same physical meaning as the variable Z is presented in monograph (Stevanovic 2010) and another paper submitted for publication (Djordjević and Stevanović 2012). Now the equation (2.9) is transformed into:

(2.11)  \frac{H_{i} − 1}{Λ} [1 + F (Kn_{e} Z)]\frac{d}{dZ} (lnP) = Z^{2}(1 − mZ)[1 + Z \frac{d}{dZ} (lnP)],

where

(2.12)  F (Kn_{e} Z) = 6 Kn_{e} Z(A_{1} + 2A_{2} Kn_{e} Z)

and m = Ṁ/Ṁ_{c}. For the solution of this equation two boundary conditions are available (see Fig. 1):

(2.13)  H = H_{i}, P = 1 ⇒ Z = 1/H_{1},  H = 1, P = 1 ⇒ Z = 1,

where indices i and e refer to the inlet and exit bearing cross sections respectively.

For some other slip velocities models presented in the literature (Barber and Emerson 2006), expressions for F (Kn_{e} Z) = F (Kn) are different, but the form of equation (2.11) remains the same. However, for all proposed boundary conditions F (0) = 0, and, as expected, equations (2.9) and (2.11) reduce to their well known form for a no-slip compressible flow through a bearing.

Equation (2.11) can be written in the form in which the variables P and Z are separated:

(2.14)  d(lnP) = \frac{(1 − mZ)dZ}{\frac{H_{i} −1}{Λ} + (6A_{1} Kn_{e} \frac{H_{i} −1}{Λ} − 1)Z + (12A_{2} Kn_{e}^{2} \frac{H_{i} −1}{Λ} + m)Z^{2}}

and its solution can be obtained by quadratures. For example, when applying the second of boundary conditions (2.13) we get:

(2.15)  lnP = ∫_{Z}^{1} \frac{(1 − mt)}{a + bt + ct^{2}} dt,

where

A = \frac{H_{i} − 1}{Λ} , b = 6A_{1} Kn_{e} \frac{H_{i} − 1}{Λ} − 1, c = 12A_{2} Kn_{e}^{2} \frac{H_{i} − 1}{Λ} + m

Application of the first of boundary conditions (2.13) in (2.15) leads to:

(2.16)  ∫_{1/H_{i}}^{1} \frac{(1 − mt)}{a + bt + ct^{2}} dt = 0

Provided the bearing number Λ, the reference Knudsen number Kn_{e}, and the ratio of the inlet and exit microbearing height are known, parameter m could be calculated by iteration from equation (2.16), which enable the determination of P and Z from (2.15). Integral (2.15) has two possible solutions depending on whether the discriminant of the a + bt + ct^{2} = 0 is positive or negative (Gradshteyn and Ryzhik 1965). If D := b^{2} − 4ac ⩾ 0 the solution is: (2.17)

lnP = \frac{m}{2c} ln \frac{a + bZ + cZ^{2}}{a + b + c} + ( 1 + \frac{mb}{2c}) \frac{1}{\sqrt{D}} ln \frac{(b + 2c − \sqrt{D} )(b + 2cZ + \sqrt{D} )}{(b + 2c + \sqrt{D} )(b + 2cZ − \sqrt{D} )}

Figure 2. The pressure distribution in the bearing obtained with the presented analytical solution and by the Boltzmann equation (Fukui and Kaneko 1988) for continuum flow conditions (Kn = 0), H_{i} = 2 and: Λ = 1, Λ = 10, Λ = 100

The expression for determining the parameter m is found by putting the first of boundary conditions (2.13), which corresponds to the microbearing exit (P = 1, Z = 1), into (2.17) (2.18)

m = \frac{2c}{\sqrt{D}} ( 1 + \frac{mb}{2c}) ln^{−1} \frac{H_{i}^{2}(a + b + c)}{aH_{i}^{2} + bH_{i} + c} ln \frac{(b + 2c − \sqrt{D} )(bH_{i} + 2c + H_{i} \sqrt{D} )}{(b + 2c + \sqrt{D} )(bH_{i} + 2c − H_{i} \sqrt{D} )}

For the case when D = b^{2} − 4ac < 0 the solution of (2.15) is

(2.19) lnP = \frac{m}{2c} ln \frac{a + bZ + cZ^{2}}{a + b + c} + ( 1 + \frac{mb}{2c}) \frac{2}{\sqrt{−D}} arctg \frac{(1 − Z)\sqrt{−D}}{2a + b + Z(b + 2c)}

Now the parameter m is found in the same way:

(2.20)  m = \frac{4c + 2mb}{\sqrt{−D}} ln^{−1} \frac{H_{i}^{2}(a + b + c)}{aH_{i}^{2} + bH_{i} + c} arctg \frac{(H_{i} − 1)\sqrt{−D}}{H_{i}(2a + b) + (b + 2c)}

Precisely speaking, the parameter m is determined from (2.18) or (2.20) iteratively by supposing the initial value for m taking into account whether b^{2} − 4ac is positive or negative during the calculation. Then, for the value of Z between Z = 1/H_{i} at the bearing inlet and Z = 1 at the bearing outlet, the pressure variation is found from (2.17) or (2.19).

3. Results and discussion

The solution is at first verified for the continuum flow regime which is independent of different slip boundary conditions. It follows from the general solution for the pressure distribution in the microbearing under the slip flow regime, i.e., from equations (2.17), (2.18), (2.19) and (2.20) with the expressions for b and c obtained for Kn_{e} = 0. The results are presented in Fig. 2 for bearing geometry defined with H_{i} = 2 and three values of the bearing number Λ = 1, Λ = 10 and Λ = 100, and they are completely in agreement with the numerical solution of the Boltzmann equation (Fukui and Kaneko 1988).

Figure 3. The pressure distribution in the microbearing obtained with the presented analytical solution and different slip coefficients in the boundary conditions and with the Boltzmann equation (Fukui and Kaneko 1988) for Λ = 1, H_{i} = 2 and: Kn_{e} = 0.1, Kn_{e} = 0.2, Kn_{e} = 0.5

In Figs. 3 and 4 the pressure distribution is presented respectively for bearing number Λ = 1, Λ = 10. All results are obtained for the same ratio of the inlet and outlet microbearing height H_{i} = 2 and three Kn_{e}values (Kn_{e} = 0.1, Kn_{e} = 0.2 and Kn_{e} = 0.5). The values of parameter m that correspond to pressure distributions presented in Figs 3 and 4, are calculated by use of equations (2.18) and (2.20), and given respectively in Tables 2 and 3.

m  Λ = 1
Kn_{e}=0.1 Kn_{e}=0.2 Kn_{e}=0.5
Schamberg (1947): A_{1} =1, A_{2} =5π/12  1.379  1.387  1.404
Deissler (1964): A_{1} =1, A_{2} =9/8  1.378  1.386  1.402
Hsia and Domoto(1983): A_{1} =1, A_{2} =0.5  1.376  1.381  1.393
Maxwell (1879): A_{1} =1, A_{2} =0  1.375  1.377  1.381
Beskok, Karniadakis (1996): A_{1} =1, A_{2} =−0.5  1.374  1.373  1.362

Table 2. Parameter m values for different slip coefficients and flow conditions presented in Fig. 3 (defined with Λ = 1, H_{i} = 2 and three Knudsen number values: Kn_{e} = 0.1, Kn_{e} = 0.2, Kn_{e} = 0.5)

Figure 4. The pressure distribution in the microbearing obtained with the presented analytical solution and different slip coefficients in the boundary conditions and with the Boltzmann equation (Fukui and Kaneko 1988) for Λ = 10, H_{i} = 2 and: Kn_{e} = 0.1, Kn_{e} = 0.2, Kn_{e} = 0.5

The presented results show the reliability of obtained analytical solution for the slip flow regime (Kn_{e} = 0.1), as well as for the part of the transitional regime (Kn_{e} = 0.2, Kn_{e} = 0.5). The second order boundary condition defined by Schamberg (1947) leads to the best fit of the analytical solution with the numerical solution of the Boltzmann equation obtained by Fukui and Kaneko (1988) in the slip regime

m  Λ = 10
Kn_{e}=0.1 Kn_{e}=0.2 Kn_{e}=0.5
Schamberg (1947): A_{1} =1, A_{2} =5π/12  1.589  1.543  1.475
Deissler (1964): A_{1} =1, A_{2} =9/8  1.591  1.545  1.476
Hsia and Domoto(1983): A_{1} =1, A_{2} =0.5  1.593  1.551  1.484
Maxwell (1879): A_{1} =1, A_{2} =0  1.595  1.556  1.493
Beskok, Karniadakis (1996): A_{1} =1, A_{2} =−0.5  1.598  1.562  1.506

Table 3. Parameter m values for different slip coefficients and flow conditions presented in Fig. 4 (defined with Λ = 10, H_{i} = 2 and three Knudsen number values: Kn_{e} = 0.1, Kn_{e} = 0.2, Kn_{e} = 0.5)

(Kn_{e} = 0.1), while for the beginning of the transition flow regime (Kn_{e} = 0.2) the Deissler (1964) boundary condition is the most appropriate. For the higher Knudsen number value (Kn_{e} = 0.5) the analytical solution obtained with Hsia and Domoto’s (1983) slip coefficients is in a good agreement with the numerical solution of the Boltzmann equation. Thus, it is confirmed that the analytical solution is valid even for a higher Knudsen number value up to Kn_{e} = 0.5.

For all results presented in Figs. 3 and 4 i.e., for the bearing number values Λ = 1, and Λ = 10 and the Knudsen numbers Kn_{e} = 0.1, Kn_{e} = 0.2 and Kn_{e} = 0.5 the boundary condition defined by Beskok and Karniadakis (1996) gives pronounced deviation from the Fukui and Kaneko (1988) results.

Also, the analytical solutions which correspond to the Maxwell first order boundary condition are presented in Figs. 3 and 4. It is obvious that Schamberg (1947), Deissler (1964) and Hsia and Domoto (1983) second order boundary conditions provide higher accuracy than the Maxwell (1879) first order boundary condition.

4. Conclusion

The appropriate transformation of the slip corrected Reynolds lubrication equation by introducing new independent variable Z leads us to the exact analytical solution. Moreover, it is shown that the variable Z has also a firm physical meaning - it represents the ratio between the local value of the Knudsen number and its value at the microbearing exit. The obtained analytical solution is presented for both the compressible slip corrected Reynolds lubrication equation and the classical compressible Reynolds lubrication equation for continuum flow conditions. The analytical solutions of these equations have not been reported in the open literature.

The slip flow results for a wide range of the Knudsen number and the continuum flow conditions, provided by the general analytical solution from this paper, are in excellent agreement with Fukui and Kaneko’s (1988) numerical solution of the Boltzmann equation.

This result is important since it provides the mathematically exact problem solution and the benchmark for the validation of numerical and experimental methods.

Acknowledgments. This work was supported by the Ministry of Education and Science of the Republic of Serbia (Grant 174014).

References

1. R. E. Barber, D. R. Emerson, Challenges in modelling gas-phase flow in microchannels: from slip to transition, Heat Transfer Engineering 27 (2006), 3–12.

2. A. Beskok, G. Karniadakis, W. Trimmer, Rarefaction and Compressibility Effects in Gas Microflows, J. Fluids Eng. 118 (1996), 448–456.

3. R. G. Deissler, An analysis of second-order slip flow and temperature-jump boundary conditions for rarefied gases, J. Heat Mass Transfer 7 (1964), 681–694.

4. N. D. Stevanovic, Strujanje razređenih gasova u mikrokanalima, Mašinski fakultet Univerziteta u Beogradu, Beograd, (2010).

5. V. D. Djordjevic, N. D. Stevanovic Slip-corrected Reynolds lubrication equation possesses the exact analytical solution, Tribology International (2012) (submitted).

6. S. Fukui, R. Kaneko, Analysis of ultra-thin gas film lubrication based on linearized Boltzmann equation. First report-derivation of a generalized lubrication equation including thermal creep flow, J. Tribol. 110 (1988), 253–262.

7. M. Gad-El-Hak, The MEMS Handbook, CRC Press, Boca Raton, Florida, (2002).

8. Y. Hsia, G. Domoto, An experimental investigation of molecular rarefaction effects in gaslubricated bearings at ultra low clearances, J. Lubr. Technol. 105 (1983), 120–130.

9. G. Karniadakis, A. Beskok, Micro Flows: Fundamentals and Simulation, Springer, New-York, (2002).

10. D. A. Lockerby, J. M. Reese, D. R. Emerson, R. W. Barber, Velocity boundary condition at solid wall in rarefied gas calculation, Phys. Rev. E 70 (2004), 017303.

11. J. C. Maxwell, On stresses in rarefied gases arising from inequalities of temperature, Phil. Trans. R. Soc. London 170 (1879), 231–256.

12. I. S. Gradshteyn, I. M. Ryzhik, Table of integrals, series, and products, Academic Press, New York, 1965.

13. R. Schamberg, The fundamental differential equations and the boundary conditions for high speed slip-flow, and their application to several specific problems, Dissertation, California Institute of Technology, 1947. http://resolver.caltech.edu/CaltechETD:etd-12272004-161736

14. A. Z. Szeri, Fluid Film Lubrication: Theory and Design, Cambridge University Press, Cambridge, 1998.

15. G. W. Vinsenti, H. C. Kruger, Introduction to Physical Gas Dynamics, R. E. Krieger, Malabar, Florida, 1986.

16. M. Von Smoluchowski, Ueber Warmeleitung in verdunnten Gasen, Ann. Phys. Chem. 64 (1898), 101–130.

Faculty of Mechanical Engineering (Received 16 10 2011) University of Belgrade 11120 Belgrade, Serbia

 nevenastev@eunet.rs

Serbian Academy of Sciences and Arts 11000 Belgrade Serbia