数理解析研究所講究録
第 824 巻 1993 年 64-93

64

Numerical analysis of a uniform flow of a rarefied gas past asphere on the basis of the Boltzmann equation
for hard-sphere molecules
Shigeru Takata, Yoshio Sone, and Kazuo Aoki
京大・工
青木一生
高田滋 曽根良夫
Department of Aeronautical Engineering, Kyoto University,
Kyoto 606-01, Japan

A slow uniform flow of a rarefied gas past a sphere is investigated on the basis of the linearized Boltzmann equation for hard-sphere molecules and the diffuse reflection condition. With
the aid of a similarity solution, the Boltzmann equation is reduced to two simultaneous integrodifferential equations with three independent variables, which are solved numericaUy. The collision integral is computed efficiently by the use of a numerical collision kernel [Phys. Fluids A
1, 363 (1989)]. The velocity distribution function of the gas molecules, which has discontinuity
in the gas, the density, flow velocity, and temperature fields of the gas, and the drag on thesphere are obtained accurately for the whole range of the Knudsen number. In spite of slow
flow, the temperature is nonuniform (thermal polarization). From the behavior of the velocity
distribution function, the kinetic transition region is clearly seen to separate into the Knudsen
layer and the layer for small Knudsen numbers.
$S$

I.

INTRODUCTION

A slow uniform flow of a rarefied gas past a sphere, or motion of a very small particle in agas, is one of the simplest and most fundamental problems in rarefied gas dynamics, and aerosol
sciences, and has interested many scientists and engineers.1-16 Millikan’s $experiments^{1,3}$ arefamous among the works on this subject. Cercignani et al.8 analyzed the problem by their variational method17 on the basis of the Boltzmann-Krook-Welander (BKW or BGK) $equation^{18,19}$and obtained the drag on the sphere for the whole range of the Knudsen number. $2$ Their results
agree very well with recent more detailed computations although a very simple test function
was used. The behavior of the gas around the sphere is studied in Refs. 10,11,12,14, and 15.
In Refs. 10 and 15, the asymptotic behavior of the gas for small Knudsen numbers is obtained
analytically up to the second order of the Knudsen number with the aid of the asymptotic$theory^{21,22}$ of the Boltzmann system. According to Refs. 10 and 15, in contrast to the solutionof classical fluid dynamics, the temperature around the sphere is nonuniform even in a slow flow,which is called thermal polarization.i3 The exphcit expression of the nonuniform temperaturefield is given in Ref. 15. In Refs. 12 and 14, Loyalka et al. carried out numerical computation of the simultaneous integral equations for macroscopic variables derived from the BKW
equation and diffuse reflection condition. There is some discrepancy between the results of thetemperature field in Refs. 14 and 15. The information obtained so far is for the BKW model
equation, and very few have been done on this subject on the basis of the standard Boltzmann
equation; recently, hmited results on the drag have been reported in Refs. 23 and 24. In thispaper we will carry out accurate and detailed numerical study of the problem on the basis of the
Boltzmann equation for hard-sphere molecules and clarify the comprehensive behavior of the65flow field for the whole range of the Knudsen number. It is shown in Refs. 25 and 26 that thevelocity distribution function of gas molecules has discontinuity in a gas around a convex body.
The discontinuity not only is interesting to be investigated accurately but also should be treatedcarefully to obtain the accurate information of the macroscopic variables. From the behavior of
the velocity distribution function, we will see separation of the kinetic transition region into the
Knudsen layer and the $S1ayer^{26,27}$ for small Knudsen numbers (Sec. V C).

II.

PROBLEM AND PRELIMINARY ANALYSIS

A.

Problem and notations, and temperature
Consider a uniform flow of a rarefied gas [velocity $(U_{\infty}, 0,0)$ , pressure) at rest. We will investigate the steadypast a spherical body (radius and temperature
behavior of the gas under the following assumptions: (i) The gas molecules are hard spheres ofa uniform size and undergo complete elastic coUision between themselves; (ji) The gas molecules
are reflected diffusely on the sphere; and (i\"u) The flow speed
is so small compared with the
speed of the sound that the basic equation (the Boltzmann equation for hard-sphere molecules)
and the boundary condition can be linearized around a uniform equilibrium state at rest.
; $R$ is the specffic gas
We summarize the main notations used in this paper:
constant (the Boltzmann constant divided by the mass of the molecule);
is the mean freepath of the gas molecules at the equilibrium state at rest with pressure
and temperature, where and $m$ are the diameter and[for a hard-sphere molecular gas,$=l_{\infty}/L$(Knudsen number); $k_{\infty}=\pi^{1/2}Kn/2;Lx$ ; is themass of the molecule, respectively]; Kn
Cartesian coordinate system with its origin at the center of the sphere and with its
axis in the$r=0$
direction of the uniform flow;is the polar coordinate system withat the center ofis the molecular
direction;
the sphere and with $\theta=0$ (the polar direction) in the, and components of
are, respectively, thevelocity; ( $=|(i|=(\zeta_{i^{2}})^{1/2};(_{r},$ ( , and$E(\zeta)=\taur^{-3/22}\exp(-();\rho_{\infty}(2RT_{\infty})^{-3/2}E(1+\phi)$is the velocity
in the polar coordinate system;
is the density of the gas; $T_{\infty}(1+\tau)$ is the
distribution function of the gas molecules;
is the flow velocity;is thetemperature; $p_{\infty}(1+P)$ is the pressure;is Kronecker’s delta. The
stress tensor; andis the heat flow vector, where, and, and,
system are expressed by the subscripts
incomponents of $u;,$(e.g.,). (See Fig. 1.)$p_{\infty}$$L$$\tau_{\infty}]$$\tau_{\infty}$$(U_{\infty})$$\rho_{\infty}=p_{\infty}/RT_{\infty}$$l_{\infty}$$T_{\infty}$$p_{\infty}$$l_{\infty}=(\sqrt{2}\pi\sigma^{2}\rho_{\infty}/m)^{-1}$$\sigma$$x_{1}$$(rL, \theta, \varphi)$$(2RT_{\infty})^{1/2}\zeta_{i}$$x_{1}$$r,$$\zeta_{\varphi}$$\theta$$\theta$$\zeta_{i}$$\varphi$$\rho_{\infty}(1+\omega)$$(2RT_{\infty})^{1/2}u_{i}$$p_{\infty}(\delta_{ij}+P_{ij})$$p_{\infty}(2RT_{\infty})^{1/2}Q_{1}$$P_{ij}$$u_{r},$$Q_{i}$$\delta_{ij}$$(r, \theta, \varphi)$$r,$$\theta$$\varphi$$u_{\theta}$

Basic equation and boundary conditionB.

The linearized Boltzmann equation for a steady state is written as:($i^{\frac{\partial\phi}{\partial x_{i}}}= \frac{1}{k_{\infty}}\mathcal{L}(\phi)$.

(1)

For a hard-sphere molecular gas, the linearized coUision integral:$\mathcal{L}(\phi)$

is expressed in the following$form$$\mathcal{L}(\phi)=\mathcal{L}_{1}(\phi)-\mathcal{L}_{2}(\phi)-\nu(\zeta)\phi$

,

(2)

,

$\mathcal{L}_{1}(\phi)=\frac{1}{\sqrt{2}\pi}\int\frac{1}{|\zeta_{1}-\xi_{i}|}\exp(-\xi_{j}^{2}+\frac{(\epsilon_{1jk}\zeta_{j}\xi_{k})^{2}}{((h-\xi_{h})^{2}})\phi(x;, \xi_{i})d\xi_{1}d\xi_{2}d\xi_{3}$$\mathcal{L}_{2}(\phi)=\frac{1}{2\sqrt{2}\pi}\int|(i-\xi;|\exp(-\xi_{j^{2}})\phi(x_{i}, \xi_{i})d\xi_{1}d\xi_{2}d\xi_{3},$

(3a)
(3b)

66

$\nu(\zeta)=\frac{1}{2\sqrt{2}}[\exp(-\zeta^{2})+(2\zeta+\frac{1}{\zeta})\int_{0}^{(}\exp(-\xi^{2})d\xi]$

.

(3c)

Here
is Eddington’s epsilon (thus,
is the vector product of ( and ); and the domain, andof integration in
the following integrals with respect to the molecular velocity(( or ) is the whole molecular velocity space unless otherwise stated.
, etc. are given by the moments of :The (nondimensional) macroscopic variables
$u;,$

au

$\epsilon_{1jk}$$\mathcal{L}_{1},$$\mathcal{L}_{2}$$\epsilon_{ijk}(j\xi_{k}$$\xi_{i}$$i$$\xi_{i}$$i$$\phi$$\tau$$\omega,$(4a)$\omega=\int\phi Ed\zeta_{1}d\zeta_{2}d(3,$$u_{i}= \int\zeta_{t}\phi Ed\zeta_{1}d\zeta_{2}d\zeta_{3}$

,

(4b)
(4c)

$\tau=\frac{2}{3}\int(\zeta_{j^{2}}-\frac{3}{2})\phi Ed(1d\zeta_{2}d\zeta_{3},$$P=\omega+\tau$

,

(4d)
(4e)

$P_{1j}=2 \int\zeta_{i}\zeta_{j}\phi Ed\zeta_{1}d\zeta_{2}d(3,$.

$Q_{i}= \int j^{2}1$(4f)

The linearized form of the diffuse reflection boundary condition on the sphere $(x_{1}^{2}=1)$ at
rest and with temperatureis given by$\tau_{\infty}$$\phi(x_{1}, \zeta;)=\sigma_{w}$,$(\zeta_{i}n;>0)$,$\sigma_{w}=-2\pi^{1/2}\int_{t:n.<0}\xi_{j}n_{j}\phi Ed\xi_{1}d\xi_{2}d\xi_{3}$

(5)

,

(6)

where is the unit normal vector to the boundary, pointed to the gas. The boundary condition$0,0)$ ], temperature, and velocity
at infinity [uniform equilibrium flow with pressure
is given as$\phi=2(1u_{\infty},$(7)$n_{i}$$(U_{\infty},$$T_{\infty}$$p_{\infty}$

with$u_{\infty}=U_{\infty}(2RT_{\infty})^{-1/2}$.

(8)

Similarity solution

C.

According to Ref. 30, the solution of the boundary-value problem, Eqs. (1), (5), and (7),can be expressed by the following similarity solution:(9)$\phi=\Phi_{c}(r, \zeta, \theta_{(})\cos\theta+(\Phi(r, \zeta, \theta_{(})\sin\theta,$$\theta_{\zeta}=\pi-Arc\cos(\zeta_{r}/\zeta)$

is the angle between the molecular velocity
are determined by the following equations:The$\Phi_{s}$,$\pi-\theta_{(}$$(0\leq\theta_{(}\leq\pi)$$\zeta_{i}$

.

(10)

and the radial direction. Then$\mathcal{D}\Phi_{c}+\frac{(\zeta\sin\theta_{\zeta})^{2}}{r}\Phi_{s}=\frac{1}{k_{\infty}}[\mathcal{L}_{1}^{c}(\Phi_{c})-\mathcal{L}_{2}^{c}(\Phi_{c})-\nu(()\Phi_{c}]$,$\mathcal{D}(\Phi_{s}\zeta\sin\theta_{\zeta})-\frac{\zeta\sin\theta_{\zeta}}{r}\Phi_{c}=\frac{1}{k_{\infty}}[\mathcal{L}_{1}^{s}(\Phi_{s}\xi\sin\theta_{\xi})-\mathcal{L}_{2}^{s}(\Phi_{s}\xi\sin\theta_{\xi})-\nu(\zeta)\Phi_{s}\zeta\sin\theta_{\zeta}]$$\Phi_{c}$

and

(11)

,

(12)

67
where
$D \Phi=-\zeta\cos\theta_{\zeta}\frac{\partial\Phi}{\partial r}+\frac{(\sin\theta_{\zeta}\partial\Phi}{r\partial\theta_{\zeta}}$

,

(13)

$\mathcal{L}_{1}^{c}(\Phi)=\frac{1}{\sqrt{2}\pi}\int_{0}^{\infty}\int_{0}^{\pi}\int_{0}^{2\pi}\frac{\xi^{2}\sin\theta_{\xi}}{\mathcal{F}_{1}}\exp(-\xi^{2}+\frac{\mathcal{F}_{2}}{\mathcal{F}_{1}^{2}})\Phi(r, \xi, \theta_{\xi})d\overline{\psi}d\theta_{\xi}d\xi$,,$\mathcal{L}_{2}^{c}(\Phi)=\frac{1}{2\sqrt{2}\pi}\int_{0}^{\infty}\int_{0}^{\pi}\int_{0}^{2\pi}\mathcal{F}_{1}\xi^{2}\sin\theta_{\xi}\exp(-\xi^{2})\Phi(r, \xi, \theta_{\xi})d\overline{\psi}d\theta_{\xi}d\xi$$\mathcal{L}_{1}^{s}(\Phi)=\frac{1}{\sqrt{2}\pi}\int_{0}^{\infty}\int_{0}^{\pi}\int_{0}^{2\pi}\frac{\xi^{2}\sin\theta_{\xi}\cos\overline{\psi}}{\mathcal{F}_{1}}\exp(-\xi^{2}+\frac{\mathcal{F}_{2}}{\mathcal{F}_{1^{2}}})\Phi(r, \xi, \theta_{\xi})d\overline{\psi}d\theta_{\xi}d\xi$$\mathcal{L}_{2}^{s}(\Phi)=\frac{1}{2\sqrt{2}\pi}\int_{0}^{\infty}\int_{0}^{\pi}\int_{0}^{2\pi}\mathcal{F}_{1}\xi^{2}\sin\theta_{\xi}\cos\overline{\psi}\exp(-\xi^{2})\Phi(r, \xi, \theta_{\xi})d\overline{\psi}d\theta_{\xi}d\xi$-2$\cos\theta_{\zeta}\sin\theta_{\zeta}\cos\theta_{\xi}\sin\theta_{\xi}\cos\overline{\psi}$

(14b)
,

(15a)

,

(15b)

(16a)

$\mathcal{F}_{1}=[(2+\xi^{2}-2(\xi(\cos\theta_{\zeta}\cos\theta_{\xi}+\sin\theta_{\zeta}\sin\theta_{\xi}\cos\overline{\psi})]^{1/2},$$\mathcal{F}_{2}=\zeta^{2}\xi^{2}[\cos^{2}\theta_{\zeta}\sin^{2}\theta_{\xi}+\sin^{2}\theta_{\zeta}\cos^{2}\theta_{\xi}+\sin^{2}\theta_{\zeta}\sin^{2}\theta_{\xi}\sin^{2}$

(14a)

di

].

(16b)

corresponds to the polar coordinate expression of ; in Eqs. (3a) and (3b),The setwith the polar direction in the radial direction.
The boundary conditions forandare:at $r=1$ ,$(\xi, \pi-\theta_{\xi},\overline{\psi})$$\xi$$\Phi_{s}$$\Phi_{c}$$\Phi_{c}$$=$,$2 \pi^{3/2}\int_{0}^{\infty}\int_{0}^{\pi/2}\zeta^{3}\sin 2\theta_{\zeta}\Phi_{c}Ed\theta_{\zeta}d($$(\pi/2<\theta_{\zeta}\leq\pi)$$\Phi_{s}$and as$rarrow\infty$

$=$

$0$

,

(17)

,

,
$arrow$$\Phi_{c}$$arrow$$\Phi_{s}$$-2\zeta\cos\theta_{\zeta}u_{\infty}$$-2u_{\infty}$

,

(18)

.

The boundary-value problem, Eqs. (1), (5), and (7), for a single integro-differential equationof six independent variables is reduced to that, Eqs. (11), (12), (17), and (18), for two simultaneous integro-differential equations of three independent variables, which will be analyzed
numerically in the next section.
, etc. are expressed by
as follows:andThe macroscopic variables$\omega,$$u_{r},$$u_{\theta}$$\Phi_{s}$$\Phi_{c}$

(19a)

$=$

$2 \pi(\int_{0}^{\infty}\int_{0}^{\pi}(^{2}\sin\theta_{\zeta}\Phi_{c}Ed\theta_{\zeta}d\zeta)\cos\theta,$$u_{r}$$=$$- \pi(\int_{0}^{\infty}\int_{0}^{\pi}\zeta^{3}\sin 2\theta_{\zeta}\Phi_{c}Ed\theta_{\zeta}d\zeta)\cos\theta$$u_{\theta}$$=$$\pi(\int_{0}^{\infty}\int_{0}^{\pi}\zeta^{4}\sin^{3}\theta_{\zeta}\Phi_{s}Ed\theta_{\zeta}d\zeta)\sin\theta$$u_{\varphi}$

$=$

$0$

$=$

$\frac{4}{3}\pi[\int_{0}^{\infty}\int_{0}^{\pi}\zeta^{2}(\zeta^{2}-\frac{3}{2})\sin\theta_{\zeta}\Phi_{c}Ed\theta_{\zeta}d(]\cos\theta,$$\omega$$\tau$

,

,

,

(19b)
(19c)
(19d)
(19e)

68

$P_{rr}$$=$$4 \pi(\int_{0}^{\infty}\int_{0}^{\pi}(^{4}\cos^{2}\theta_{\zeta}\sin\theta_{\zeta}\Phi_{c}Ed\theta_{\zeta}d()\cos\theta,$(19f)$P_{\theta\theta}$$=$$2 \pi(\int_{0}^{\infty}\int_{0}^{\pi}(^{4}\sin^{3}\theta_{\zeta}\Phi_{c}Ed\theta_{(}d\zeta)\cos\theta,$(19g)$P_{r\theta}$$=$$-2 \pi(\int_{0}^{\infty}\int_{0}^{\pi}\zeta^{5}\cos\theta_{\zeta}\sin^{3}\theta_{\zeta}\Phi_{s}Ed\theta_{\zeta}d()\sin\theta,$(19h)$P_{\varphi\varphi}$$=$$3(\omega+\tau)-P_{rr}-P_{\theta\theta}$$Q_{r}$$=$$- \pi[\int_{0}^{\infty}\int_{0}^{\pi}\zeta^{3}(\zeta^{2}-\frac{5}{2})\sin 2\theta_{\zeta}\Phi_{c}Ed\theta_{\zeta}d(]\cos\theta,$$Q_{\theta}$$=$$\pi[\int_{0}^{\infty}\int_{0}^{\pi}\zeta^{4}((-\frac{5}{2})\sin\theta_{\zeta}\Phi_{s}Ed\theta_{(}d\zeta]\sin\theta,$$Q_{\varphi}$

$=$

$0$

,

$P_{r\varphi}=P_{\theta\varphi}=0$

,

(19i)
(19j)
(19k)

.

(191)

Discontinuity of the velocity distribution function

D.

The Boltzmann equation (1) defines the variation of in the direction of in the ; space.
, it propagates along the characteristic
If there is a discontinuity of at $x;=x_{i}^{(0)}$ and (, where $t(>0)$ is a parameter. Since theterms of the collision integral
andare continuous, the vaniation of the discontinuity is determined by$\phi$$\zeta_{t}$$x$$i=\zeta_{i}^{(0)}$$\phi$$x_{i}=x_{i}^{(0)}+(_{t}^{(0)}t$$\mathcal{L}_{1}$$\mathcal{L}_{2}$$\mathcal{L}(\phi)$($i \frac{\partial[\phi]}{\partial x_{i}}=-\frac{1}{k_{\infty}}\nu(\zeta)[\phi]$

,

(20)

where
is the difference of on both sides of the discontinuity. The discontinuity decays in
the length scale of the free path of the molecule with velocity ( The position of discontinuity
is independent of molecular speed.
On a boundary the velocity distribution function of the molecules leaving the boundary is
specified. That is,( $(|n;>0)$ .
(21)$[\phi]$$\phi$$i\cdot$$\phi=\sigma_{w}$The feature ofis quite different from with ($in_{1}<0$ , which is established by collisionsof surrounding gas molecules. Thus, on the boundary is generally discontinuous at $(in;=0$ .
When the boundary is convex, this discontinuity propagates into the gas along the characteristicof Eq. (1). When the boundary is concave, on the other hand, the characteristic, which is tangent
to the boundary, does not enter into the gas and therefore no discontinuity propagates into thegas. Thus at a given pointin the gas, is discontinuous at ; whose opposite vector(Fig. 2).
lies on the circular cone along which the edge of the sphere is viewed fromspace or directly from the characteristics of Eqs. (11) and
By transformation into ( ,space is found to lie on the common surface:(12), the discontinuity ofandin$\phi$$\sigma_{w}$$\phi$$\phi$$x_{i}$$(-\zeta_{t})$$\zeta$$x_{i}$$r$$\Phi_{c}$$(, \theta_{\zeta})$$\Phi_{s}$$(r, \zeta, \theta_{\zeta})$$r\sin\theta_{\zeta}=1$,$(\pi/2\leq\theta_{\zeta}\leq\pi)$

.

(22)

In the following numerical analysis ofand , the difference formulae of differentiation that
contain the data on both sides of the discontinuity should be avoided.
More about the discontinuity of the velocity distribution function, especially its relation to
the layer27 at the bottom of the Knudsen layer, are found in Refs. 25, 26 and 31.
$\Phi_{c}$$S$$\Phi_{s}$69
ASYMPTOTIC SOLUTIONS FOR LARGE AND SMALL KNUDSEN NUMBERS

III.

Before proceeding to numerical analysis of the problem, we give the asymptotic solution of
the boundary-value problem, Eqs. (11), (12), (17), and (18), for two limiting cases
and$k_{\infty}arrow\infty$$k_{\infty}arrow 0$.

, where the collision term of the Boltzmann equation
For the free molecular limit
is neglected, the solution is obtained as follows:$(k_{\infty}=\infty)$$\frac{\Phi_{c}}{u_{\infty}}=\{\sqrt{\pi}[\cos\theta^{(}(1-r^{2}\sin^{2}\theta_{\zeta})^{1/2}-r\sin^{2}\theta_{\zeta}]-2\zeta\cos\theta_{\zeta},$$[ \pi-A^{\zeta}rc\sin(\frac{1}{r})<\theta_{\zeta}\leq^{\frac{1}{r}}\pi][0\leq\theta<\pi-Arc\sin()],$

$\}$

(23)

,

$\frac{\Phi_{s}}{u_{\infty}}=\{\sqrt{\pi}^{2}(--1[r\cos\theta_{\zeta}+(1-r^{2}\sin^{2}\theta_{(})^{1/2}],$$\}$$[ \pi-A^{\zeta}rc\sin(\frac{1}{r})<\theta\leq^{\frac{1}{r}}\pi][0\leq\theta<\pi-Arcsi_{\zeta}n()],$

,

(25a)

$\frac{\omega}{u_{\infty}\cos\theta}$$=$$- \frac{\sqrt{\pi}}{6}\{(1+\frac{6}{\pi})r^{-2}+[2-(2+r^{-2})(1-r^{-2})^{1/2}]r\}$$\frac{u_{r}}{u_{\infty}\cos\theta}$$=$$\frac{1}{2}+\frac{1}{2}(1-r^{-2})^{3/2}-\frac{1}{4}r^{-3}-\int_{0}^{1/r}t(1-t^{2})^{1/2}(1-r^{2}t^{2})^{1/2}dt$$\frac{u_{\theta}}{u_{\infty}\sin\theta}$$=$$- \frac{1}{2}-\frac{1}{4}(1-r^{-2})^{1/2}(2+r^{-2})-\frac{1}{8}r^{-3}$,$+ \frac{1}{2}\int_{0}^{1/r}t^{3}(1-t^{2})^{-1/2}(1-r^{2}t^{2})^{1/2}dt$$=$$\frac{\tau}{u_{\infty}\cos\theta}$$- \frac{1}{3\sqrt{\pi}}r^{-2}$

(24)

,

(25b)

(25c)

.

(25d)

[up to] is easily obtained with the aid of the
The asymptotic solution for small
asymptotic theory in Re&. 26, 31, 32 and 33 as follows:$O(k_{\infty})$$k_{\infty}$$\gamma_{1},$$\kappa_{0}$

(26a)

$=$

$- \frac{3}{2}\gamma_{1}k_{\infty}r^{-2}$$\frac{u_{r}}{u_{\infty}\cos\theta}$$=$$1- \frac{3}{2}(1+\kappa_{0}k_{\infty})r^{-1}+\frac{1}{2}(1+3\kappa_{0}k_{\infty})r^{-3}$$\frac{u_{\theta}}{u_{\infty}\sin\theta}$$=$$-1+ \frac{3}{4}(1+\kappa_{0}k_{\infty})r^{-1}+\frac{1}{4}(1+3\kappa_{0}k_{\infty})r^{-3}+\frac{3}{2}k_{\infty}Y_{0}(\eta)$

$=$

$0$

$\tau$

where,$\frac{\omega}{u_{\infty}\cos\theta}$, and$Y_{0}(\eta)$

,

,

(26b)

,

(26c)
(26d)

are determined by the molecular model. That is,$\gamma_{1}=1.270042$$\gamma_{1}=1$,,$\kappa_{0}=-1.2540$,$\kappa_{0}=-1.01619$(hard-sphere molecules),

,

(27)

(BKW model),, called Knudsen-layer function for the shear flow, is a function of the stretched coandordinate defined by $\eta=(r-1)/k_{\infty}$ and tabulated in Ref. 34 for hard-sphere molecules and$Y_{0}(\eta)$$\eta$

70(Ref. 34)).
(here, Refs. 26, 32, $33)=[-S(x_{1}),$in Ref. 35 for the BKW model (andare related to the viscosity and the slip coefficient of the shear flow,
The constants,(Ref. 34)]. The pressure is obtained up to
respectively [Re&. 26, 31, etc.;
except for the Knudsen-layer correction, as$[Y_{0}(\eta),$$\gamma_{1}$$x_{1}]$$\eta]$$\kappa_{0}$$\kappa_{0}=-\beta_{A}$$O(k_{\infty}^{2})$$\frac{P}{u_{\infty}\cos\theta}=-\frac{3}{2}\gamma_{1}k_{\infty}(1+\kappa_{0}k_{\infty})r^{-2}$

.

(28)

IV.

METHOD OF NUMERICAL ANALYSIS

A.

Process of numerical computation by a finite difference method

In the present paper we analyze the boundary-value problem, Eqs. (11), (12), (17), and(18), numerically by a finite difference method. As is pointed out in Sec. II , the velocity
distribution functions
andare discontinuous along $r\sin\theta_{\zeta}=1(\pi/2\leq\theta_{\zeta}\leq\pi)$. In orderto describe the discontinuity accurately and avoid errors owing to it, a special care is required inconstructing the difference scheme.25 The problem is considered over a finite domain $(1\leq r\leq r_{D}$ ,$0\leq(\leq\zeta_{D}, 0\leq\theta_{\zeta}\leq\pi)$ , whereandare chosen properly depending on situations. The
process of reduction to the finite domain will be discussed in Sec. IV C. Let (be$I(r^{(0)}=1,r^{(I)}=r_{D}),$the lattice points in the domain, where $i=0,1,$$j=0,1,$$D$$\Phi_{c}$$\Phi_{s}$$\zeta_{D}$$r_{D}$$r^{(i)},$$(^{(j)}, \theta_{\zeta}^{(k)})$$J$$\ldots,$$\ldots,$, and $k=0,1,$. The$f(r)$
ofof , etc. at a lattice point are denoted by the subscripts
functionscorresponding to the lattice point:$(((0)=0, \zeta^{(J)}=\zeta_{D})$$\Phi(r, \zeta, \theta_{(})$$K(\theta_{\zeta}^{(0)}=0, \theta_{(}^{(\hat{K})}=\pi/2, \theta_{\zeta}^{(K)}=\pi)$$\overline{K},$$\ldots,$$(r, \zeta, \theta_{\zeta}),$$\ldots,$$r$$f_{(i)}=f(r^{(i)})$$\Phi_{(i,j,k)}=\Phi(r((, \theta_{\zeta}^{(k)}),$In the followings we introduce a new function$\hat{\Phi}$

.

(29)

, defined by$\hat{\Phi}_{s}=\Phi_{s}\zeta\sin\theta_{\zeta}$

,

(30)

and , instead ofand carry out analysis in terms ofand . We construct the discreteof Eqs. (11), (12), (17), and (18) as the limit of the sequences
solutionsandandobtained by the iteration process described below.
$\hat{\Phi}$$\Phi_{c}$$\Phi_{c(i,j,k)}$$\Phi_{c(i,j)k)}^{(n)}$$\Phi_{c}$$\Phi_{s}$$\hat{\Phi}_{s(:,j,k)}$$\hat{\Phi}_{s(i,j,k)}^{(n)}(n=0,1,2, \ldots)$

Corresponding to Eqs. (11) and (12), the following finite difference equations for$\hat{\Phi}^{(n)}$$s(i,j,k)$$\Phi_{c(\cdot,j,k)}^{(n)}$andare adopted:$-((J) \cos\theta_{\zeta}^{(k)}\nabla_{1^{1}}^{(,j,k)}\Phi_{c}^{(n)}+\frac{((j)_{\sin\theta_{\zeta}^{(k)}}}{r^{(;)}}\nabla_{2}^{(i,j,k)}\Phi_{c}^{(n)}+\frac{\zeta^{(j)}\sin\theta_{\zeta}^{(k)}}{r^{(i)}}\hat{\Phi}_{s(i,j,k)}^{(n)}$$= \frac{1}{k_{\infty}}(C_{(i,j,k)}^{(n-1)}-\nu_{(j)}\Phi_{4i}^{(n)_{j,k)}},)$

,

(31)

,

(32)

$-((j)_{\cos\theta_{\zeta}^{(k)}\nabla_{1}^{(i,j,k)}\hat{\Phi}_{s}^{(n)}+\frac{(^{(j)}\sin\theta_{\zeta}^{(k)}}{r^{(i)}}\nabla_{2}^{(i,j,k)}\hat{\Phi}_{s}^{(n)}-\frac{\zeta^{(j)}\sin\theta_{\zeta}^{(k)}}{r^{(:)}}\Phi_{c(i,k)}^{(n)_{\theta}}}$$= \frac{1}{k_{\infty}}(\hat{S}_{(i,j,k)}^{(n-1)}-\nu_{(j)}\hat{\Phi}_{s(i,j,k)}^{(n)})$

71
where
below, and$\nabla_{1}^{(i,j,k)}$and$C_{(\dot{\cdot},j,k)}^{(n)}$$\nabla_{2}^{(i,j,k)}$andcorrespond tooperators, respectively, given exphcitly
andare the collision integrals:$\partial/\partial r$$\partial/\partial\theta_{\zeta}$$\hat{S}_{(\dot{\cdot},j,k)}^{(n)}$(33)$C_{(1j,k)}^{(n)}=[\mathcal{L}_{1}^{c}(\Phi_{c}^{(n)})-\mathcal{L}_{2}^{c}(\Phi_{c}^{(n)})]_{(i,j,k)}$(34)$\hat{S}_{(i,j,k)}^{(n)}=[\mathcal{L}_{1}^{s}(\hat{\Phi}_{s}^{(n)})-\mathcal{L}_{2}^{s}(\hat{\Phi}_{s}^{(n)})]_{(i,j,k)}$andAn efficient way of computation ofsubsection (Sec. IV B).
We use the following formulae for theand$C_{(ij,k)}^{(n_{)})}$$\hat{S}_{(\cdot,j,k)}^{(n)}$$\nabla_{1}$$\nabla_{2}$, adopted in the paper, is given in the next
operators:(35a)$(k=1)$ ,(35b)$\nabla_{2}^{(i,j,k)}\Phi=\{\begin{array}{l}(\Phi_{(i)j,1)}-\Phi_{(i,j,0)})/\Delta\theta_{\zeta}^{(1)}[(\lambda_{(k-1)}^{-1}+2)\Phi_{(i,j,k)}-(\lambda_{(k-1)}+\lambda_{(k-1)}^{-1}+2)\Phi_{(i,j,k-1)}+\lambda_{(k-1)}\Phi_{(i_{\dot{J}},k-2)}]/(\Delta\theta_{\zeta}^{(k-1)}+\triangle\theta_{\zeta}^{(k)})\end{array}$$(2\leq k\leq K)$,$\Lambda_{(i)}=\frac{\Delta r^{(i+1)}}{\Delta r^{(i)}}\Delta r^{(i)}=r^{(i)}-r^{(i-1)}$$\Delta\theta_{(}^{(k)}\theta_{\zeta}^{(k-1)}\lambda_{(k)}=^{=}\frac{\Delta\theta^{(k+1)}\theta_{\zeta_{(}}^{(k)}-}{\Delta\theta_{\zeta}^{(k)}}$

(35c)

.

. The
are discontinuous alongand
As explained in Sec. II $D,$difference formulae (35a) and (35b) should be modified when they consist of the data on bothsides of the discontinuity. That is,$\Phi_{c}$V$\hat{\Phi}_{s}$$r\sin\theta_{\zeta}=1(\pi/2\leq\theta_{\zeta}\leq\pi)$$(*,j,k)_{\Phi=}1\{\begin{array}{l}A_{1}\Phi_{(i,j)k)}-A_{2}\Phi_{(i-1,j,k)}+A_{3}\Phi_{(i(k),j,k)+},(r^{(i-1)}>1/sin\theta_{\zeta}^{(k)}\geq r^{(i-2)})(\Phi_{(i_{\dot{J}},k)}-\Phi_{(*(k),j,k)+})/(r^{(\dot{\cdot})}-1/sin\theta_{\zeta}^{(k)}),(r^{(i)}>1/sin\theta_{\zeta}^{(k)}\geq r^{(i-1)})\end{array}$

(36a)

72

(36b)

$\nabla_{2}^{(i,j,k)}\Phi=\{\begin{array}{l}B_{1}\Phi_{(*)j,k)}-B_{2}\Phi_{(|)j,k-1)}+B_{3}\Phi_{(i,j,k(e))-}[\theta_{\zeta}^{(k-1)}>\pi-Arcsin(1/r^{(i)})\geq\theta_{\zeta}^{(k-2)}](\Phi_{(|}\cdot,\dot{J}^{k)}-\Phi_{(|j,k(i))-})/[\theta_{\zeta}^{(k)}+Arcsin(1/r^{(*)})-\pi][\theta_{\zeta}^{(k)}>\pi-A\iota csin(1/r^{(|)})\geq\theta_{\zeta}^{(k-1)}]\end{array}$

where

,

(37a)

$\Phi_{(i(k),j,k)\pm}$$=$$\Phi(1/\sin\theta_{\zeta}^{(k)}\pm 0, \zeta^{(j)}, \theta_{\zeta}^{(k)})$$\Phi_{(i,j,k(i))\pm}$$=$$\Phi(r^{(i)}, \zeta^{(j)}, \pi-Arc\sin(1/r^{(i)})\mp 0)$

,

(37b)

is the
, and
, and
(
) are chosen in such a way thatandaccuracy
at
of
the
for
expression
difference
second-order, anddepend only on and . It is noted that $(i(k), j, k)$.
and $(i, j, k(i))$ are not regular lattice points. The $i(k)$ and $k(i)$ represent the intersection of theand $r=r^{(i)}$ respectively. That ischaracteristicwith the lattice lines$A_{1},$$A_{3}$$A_{2}$$B_{1},$$\nabla_{1^{*}}^{(k)}\theta,\Phi(\nabla_{2}^{(i,j,k)}\Phi)$$B_{3}$$B_{2}$$\partial\Phi/\partial r(\partial\Phi/\partial\theta_{\zeta})$$(r^{(i)}, \zeta^{(j)}, \theta_{\zeta}^{(k)})$$A_{1},$$A_{2},$$A_{3},$$B_{1},$$k$$i$$B_{3}$$B_{2}$$\theta_{\zeta}=\theta_{\zeta}^{(k)}$$r\sin\theta_{\zeta}=1$$r^{(i(k))}$$=$$\theta^{(k(*))}$$=$$1/\sin\theta_{\zeta}^{(k)}$

,

(38)

$\pi-Arc\sin(1/r^{(:)})$.

When a discontinuous function on the characteristic is considered, the two hmiting cases as in, i.e.,
Eqs. (37a) and (37b) should be considered. The formulae contain $\Phi_{(i(k),j,k)+}$ and$\Phi_{(*,j,k(i))-}$, correspondingon either side of the discontinuity. The sequences
and
to
and
at
on either side of the discontinuity, are constructed by the following difference
scheme:$\Phi_{c(:}^{(n)_{j,k(:))\pm}}$$\Phi$$r^{(i)}$$\hat{\Phi}_{s}$$\Phi_{c}$$\hat{\Phi}_{s(:}^{(n)_{j,k(*))\pm}}$$- \zeta^{(j)}\cos\theta_{(}^{(k(i))}\nabla_{1^{*}}^{(,j)\pm}\Phi_{c}^{(n)}+\frac{\zeta^{(j)}\sin\theta_{\zeta}^{(k(i))}}{r^{(:)}}\hat{\Phi}_{s(*}^{(n)_{j,k(:))\pm}}$$= \frac{1}{k_{\infty}}(C_{(i,j,k(:))}^{(n-1)}-\nu_{(j)}\Phi_{c(ij,k(i))\pm}^{(n)_{)}})$

,

(39)

,

(40)

$- \zeta^{(j)}\cos\theta_{\zeta}^{(k(*))}\nabla_{1}^{(i,j)\pm}\hat{\Phi}_{s}^{(n)}-\frac{(^{(j)}sin.\theta_{\zeta}^{(k(i))}}{r^{(*)}}\Phi_{4^{i},)}^{(n)_{jk(i))\pm}}$$= \frac{1}{k_{\infty}}(\hat{S}_{(i,j,k(i))}^{(n-1)}-\nu_{(j)}\hat{\Phi}_{s(i}^{(n)_{j,k(i))\pm}})$where$\nabla_{1}^{(i,j)\pm}\Phi$$\Phi_{(i,j,k(i))\pm}$corresponds to$\nabla_{1}^{(i,j,k)}\Phi(\overline{K}+1\leq k\leq K)$. The absence $of\pm sign$ in$C_{(*j,k(i))}^{(n_{)}-1)}$

andreplaced byin Eq. (35a) withcorresponds to the fact that$\Phi_{(i1j,k)}$$\hat{S}_{(i,j,k(i))}^{(n-1)}$$\mathcal{L}_{1}(\phi)$

are continuous.
The boundary conditions forandat $r=1(i=0)$ are obvious from Eq. (17).
As the condition at $r=r_{D}$ , instead of Eq. (18), we impose
and$\mathcal{L}_{2}(\phi)$$\Phi_{c(1j,k)}^{(n)}$$\Phi^{(n)}$$4I-1,j,k)$$\Phi_{c(I,j,k)}^{(n)}$$\hat{\Phi}_{s(i,j,k)}^{(n)}$

$=$

$=$

$G^{(n)}$$c(j,k)$’$H_{c(j,k)}^{(n)}$

,

(41)

73
$\hat{\Phi}_{s(I-1,j,k)}^{(n)}$$=$$G_{s(j,k)}^{(n)}$$\hat{\Phi}_{s(I,j,k)}^{(n)}$$=$$H_{s(j,k)}^{(n)}$

,

(42)

,

where
and
are determined by the as ymptotic solution for large$s(j,k)$ ’$c(j,k)$ ’.
in&c.
With these preparation of difference formulae, we construct the sequencesand$G^{(n)}$$H^{(n)}$$G^{(n)}$$H_{s(j,k)}^{(n)}$$r$$I\wp_{c’}^{k)}$$\Phi_{c(i,j,k)}^{(n)}$by the following process. Let$H^{(n)}$$s(j,k)$$\Phi_{4^{1,j,k)}}^{(n-1)}$and$\hat{\Phi}_{s(l,j,k)}^{(n-1)}$

be known. Then

$G_{s(j,k)}^{(n)},$$G_{4j}^{(n)_{k)}},$$\hat{\Phi}_{s(*,j,k)}^{(n)}$, and$H_{c(j,k)}^{(n)}$are known (Sec. IV C).

(i) For $0\leq k\leq\overline{K}$ , starting from, compute$\Phi_{4*}^{(n)_{j,k)}}$using Eqs. (31), (32), (41), and (42) in descending order of down to$\Phi_{c(0,j,k)}^{(n)}$and$\Phi_{4I-2,j,k)}^{(n)}$$\hat{\Phi}_{s(I-2,j,k)}^{(n)}$$i$The step $i=i+1$ to $i=i$ is as follows. Let

Starting from$\Phi_{4*}^{(n)_{j,0)}}$and$\hat{\Phi}_{s(i,j,0)}^{(n)}$, compute$\Phi_{c(i,j,k)}^{(n)}$$\Phi_{c(i,j,k)}^{(n)}$

and

and

Eqs. (41) and (42) for $i=I-2,$ $I-3$ ] in ascending order ofCarry out this step for every .

$k$

and

$\hat{\Phi}_{s(i,j,k)}^{(n)}$$\hat{\Phi}_{s(0,j,k)}^{(n)}$.

$(i‘ >i)$ be given.

$\hat{\Phi}_{s(i,j,k)}^{(n)}$$\hat{\Phi}_{s(i,j,k)}^{(n)}$and,using Eqs. (31) $-(35c)$ [and

up to

and

$\Phi^{(n)}$$c(i,j,\hat{K})$$\hat{\Phi}^{(n)}$$s(*,j,\hat{K})$$j$(ii) Compute$\Phi_{c(*}^{(n)_{j,k(i))\pm}}$$\Phi_{4^{i,k(i))+}}^{(n)_{\theta}}$

and

and

along the discontinuity using Eqs. (39) and (40):$\hat{\Phi}_{s(ij,k(i))\pm}^{(n)_{)}}$$\hat{\Phi}_{s(*}^{(n)_{j,k(i))+}}$are constructed from the initial dataand$\Phi_{c(0,j,\hat{K})}^{(n)}$$\hat{\Phi}_{s(0,j,\hat{K})}^{(n)}$;andfrom the initial condition
obtained in the preceding step
given by Eq. (17). The computation is continued until the discontinuity becomes negligibly, obtainbysmall. Then, from a set ofinterpolation. [The data on the discontinuity in Eqs. (36a) and (36b) are now prepared.]$(i)^{36}$$\Phi_{4^{1}}^{(n)_{j,k(i))-}}$$\hat{\Phi}_{s(*}^{(n)_{j,k(i))-}}$$\Phi_{4*}^{(n)_{(k),j,k)+}}$$\Phi_{4^{i}}^{(n)_{j,k(i))+}}(\hat{\Phi}_{s(i}^{(n)_{j,k(i))+}})$$(\hat{\Phi}_{s(*(k),j,k)+}^{(n)})$(i\"u) For $\overline{K}+1\leq k\leq K$ , the computation is carried out in as cending order of , independently
in the two regions $r>1/\sin\theta_{\zeta}$ and $r<1/\sin\theta_{\zeta}$ separated by the discontinuity. The
be given. (a) For
step $i=i-1$ to $i=i$ is as follows. Letand$i$$\hat{\Phi}_{s(i,j,k)}^{(n)}(i’<i)$$\Phi_{c(i,j,k)}^{(n)}$$r>1/\sin\theta_{\zeta}$, starting from$\Phi_{c(*,j,\hat{K}+1)}^{(n)}$and$\hat{\Phi}_{s(i,j,\hat{K}+1)}^{(n)}$, compute$\Phi_{c(1j,k)}^{(n)}$and$\hat{\Phi}_{s(i,j,k)}^{(n)}$using

Eqs. (31), (32), (35a), (35b), (36a), and the data obtained in steps (i) and (ii) in ascending, where. (b) Fororder of up toand$k$$\Phi_{c(i,j,\overline{k})}^{(n)}$, starting from$\hat{\Phi}_{s(l,j,\overline{k})}^{(n)}$$\theta_{\zeta}^{(\overline{k})}<\pi-Arc\sin(1/r^{(*)})\leq\theta_{\zeta}^{(\overline{k}+1)}$, compute
and
and
using
Eqs. (31), (32), (35a), (35b), (36b), the data of Eq. (17), and those obtained in step (ii)in ascending order of up to. Carry out these steps for every . Forand,$i>i_{c}$ , where the discontinuity is neghgibly small, neglecting the discontinuity, computeandfrom $k=\overline{K}+1$ to $K$ only with the standard formulae (35a) and (35b)
as in step (i).
$r<1/\sin\theta_{\zeta}$$k$$\Phi_{c(*,j,k)}^{(n)}$$\Phi_{c(*,j,\overline{k}+1)}^{(n)}$$\Phi_{4^{1K)}}^{(n)_{\dot{J}}},$$\Phi_{c(1,j,k)}^{(n)}$$\hat{\Phi}_{s(i,j,\overline{k}+1)}^{(n)}$$j$$\hat{\Phi}_{s(i,j,K)}^{(n)}$$\hat{\Phi}_{s(i,j,k)}^{(n)}$(iv) Applying the Simpson formula to Eqs. $(19b)-(19e)$ , computeand$\hat{\Phi}_{s(*,j,k)}^{(n)}$$(\tau/\cos\theta)_{(I)}^{(n)}$solution for largeand (42).

.

$(u_{r}/\cos\theta)_{(I)}^{(n)},$$(u_{\theta}/\sin\theta)_{(I)}^{(n)}$,Then by the process of connection of these data to the asymptotic$r$(Sec. IV C), determine$G_{4j,k)}^{(n+1)},$$H_{4j,k)}^{(n+1)},$(v) Repeat steps $(i)-(iv)$ with shift of the superscript ($n$$G_{s(j,k)}^{(n+1)}$, and

to $n+1$ ) until$H_{s(j,k)}^{(n+1)}$$\Phi_{4*}^{(n)_{j,k)}}$,in Eqs. (41)and$\hat{\Phi}_{s(i,j,k)}^{(n)}$74) of the problem, fromandconverge. We take the hmits as the solution (which the macroscopic variables are obtained by integration [cf. Eqs. $(19a)-(191)$ ].
$\hat{\Phi}_{s(i,j,k)}$$\Phi_{4^{i,j,k)}}$The order of computation in the preceding process is consistent with the natural course of
integration of Eqs. (11) and (12) along their characteristic in the direction of molecular velocity.

Computation of collision integralB.

efficiently, which take the majority ofandTo compute the colIision integrals
the computing time, we make use of the numerical kernel method introduced in Ref. 37. Thatis,at the nth step and $r=r^{(i)}$ are expanded in a set of basis functions:and$C_{(\dot{*},j,k)}^{(n)}$$\hat{S}_{(1j,k)}^{(n)}$$\hat{\Phi}_{s}$$\Phi_{c}$$(, \theta_{\zeta})=\sum_{\ell,m}\Phi_{c(i,\ell,m)}^{(n)}\Psi_{tm}(\zeta, \theta_{(}),$

(43)

,

(44)

$\Phi_{c}^{(n)}(r^{(i)},$$\hat{\Phi}_{s}^{(n)}(r^{(i)}, \zeta, \theta_{(})=\sum_{t,m}\hat{\Phi}_{s(1,\ell,m)}^{(n)}\Psi_{tm}(\zeta, \theta_{\zeta})$where the basis functionsin (and , takeandthe quadratic functions of andits four corners.38 Then,$\Phi_{4^{1}}^{(n)_{t,m)}}$$\theta_{\zeta}$

are continuous
are chosen in such a way thatand, respectively, at the lattice points and are expressed byin each rectangle oflattices with even lattice points at$\Phi_{c}^{(n)}$$\Psi_{\mathfrak{l}m}(\zeta, \theta_{(})$$\hat{\Phi}_{s}^{(n)}$$\hat{\Phi}_{s(i,l,m)}^{(n)}$$\zeta$$\theta_{\zeta}$$2\cross 2$

,

(45)

,

(46)

$C_{(1j,k)}^{(n)}= \sum_{t,m}\Phi_{4^{1}}^{(n)_{\ell,m)}}K_{jk\ell m}^{c}$$\hat{S}_{(*,j,k)}^{(n)}=\sum_{\ell,m}\hat{\Phi}_{s(1,t,m)}^{(n)}K_{jktm}^{s}$where$K_{j^{c}k\ell m}$and$K_{j^{s}ktm}$

are, respectively,$C_{(\dot{\iota},j,k)}^{(n)}$with$\Phi_{c}^{(n)}=\Psi_{\ell m}$$K_{jk\ell m}^{c}=[\mathcal{L}_{1}^{c}(\Psi_{lm})-\mathcal{L}_{2}^{c}(\Psi_{\ell m})]_{(j,k)}$$K_{jk\ell m}^{s}=[\mathcal{L}I(\Psi_{\ell m})-\mathcal{L}_{2}^{s}(\Psi_{\ell m})]_{(j,k)}$and$\hat{S}_{(i,j,k)}^{(n)}$with$\hat{\Phi}_{s}^{(n)}=\Psi_{tm}$

:

,

(47)

,

(48)

The
and
are universal constants, which can be computed beforehand and used invarious problems. The computations of Eqs. (47) and (48), whose effective domain of integration, as its cross section, should be carried, whereis a slender ring with the support ofout carefully since the kernels of Eqs. (14a) and (15a) have singularity.39
The discontinuities in ( space ofandcannot be describedby the expansions (43) and (44). The expanded functions are inaccurate over a lattice elements(say,, where tc depends on) for all (. These local inaccuraciesininduce only neghgible errors in the results (45) and (46) of integration over ( space. In our
computation, the discontinuity at(i.e. $r=1$ ) is taken care of by introducing two limitingand. If the discontinuity happens to pass the latticevalues 36 atfor$K_{j^{c}ktm}$$K_{j^{s}u_{m}}$$\Psi_{\ell m}$$\Psi_{\ell m}\neq 0$$\Phi_{c}^{(n)}(r^{(i)}, \zeta, \theta_{\zeta})$$\hat{\Phi}_{s}^{(n)}(r^{(i)}, \zeta, \theta_{(})$$\theta_{\zeta}$$r^{(i)}$$\theta_{\zeta}^{(2\kappa)}\leq\theta_{\zeta}<\theta_{\zeta}^{(2\kappa+2)}$$\theta_{\zeta}$$\theta_{\zeta}$$r^{(0)}$$\theta_{\zeta}=\pi/2$

points(are taken as$\theta_{\zeta^{\wedge}}^{(k)}$$k\wedge=2\kappa$$\Phi_{c(0,\ell,\hat{K})}$$\hat{\Phi}_{s(0,\ell,\hat{K})}$or $2\kappa+1$ ), the hmiting values
andand’ respectively in Eqs. (45) and (46), and therefore the expanded$\Phi_{c(*,t,k)}\wedge$$\Phi_{c}(r^{(i)},$$(^{(l)}, \theta_{\zeta^{\wedge}}^{(k)}-0)$$\hat{\Phi}_{s(,k)};,\ell^{\wedge}$

functions are inaccurate over the lattice element$\theta_{\zeta}^{(2\kappa)}\leq\theta_{\zeta}<\theta_{\zeta}^{(2\kappa+2)}$.

$\hat{\Phi}_{s}(r^{(i)}, \zeta^{(l)}, \theta_{\zeta^{\wedge}}^{(k)}-0)$

75
C.

Asymptotic behavior for large

and the condition at $r=r_{D}$$r$

We have hmited the original infinite domain to the finite one $(1\leq r\leq r_{D},$ $0\leq(\leq\zeta_{D},$. It does not introduce any problem to hmit the range of . As is easily checked, and therefore the accurate
decay very rapidly as (by computation,and. The problem iscomputation of the problem can be carried out with a reasonable size ofnot so simple for the case of the variable. From the asymptotic solutions for smalland, the approach of (or) to the uniform state at infinity is expected to beandvery slow $(\sim r^{-m})$ . Thus, a very large
is required to obtain an accurate result by simple
apphcation of Eq. (7) at $r=r_{D}$ . We therefore introduce a method to make use of the asymptotic
solution for large .
Let the deviation $\delta\phi(=\phi-2(1u_{\infty})$ from the equihibrium state at infinity be $O(r^{-m})$ for large, then. That is, for $r>r_{A}$ the characteristic length scale of variation ofis larger than , and therefore the effective Knudsen number $Kn_{e\#}(=t_{\infty}/Lr_{A}=Kn/r_{A})$ is) for large is obtained with. Thus, the asymptotic form of (orandsmall for large
the aid of the asymptotic theory
for smaJl Knudsen numbers as follows:$0\leq$$\theta_{\zeta}\leq\pi)$$\zeta$$\Phi_{s}E$$\Phi_{c}E$$arrow\infty$$\zeta_{D}$$k_{\infty}$$r$$k_{\infty}=\infty$$\phi$$\Phi_{c}$$\Phi_{s}$$r_{D}$$r$$r$$\partial\delta\phi/\partial r=r^{-1}O(\delta\phi)$$\delta\phi$$r_{A}$$\phi$$r_{A}$$=$$\frac{\Phi_{c}}{u_{\infty}}$$\Phi_{s}$$\Phi_{c}$$r$$-2(1+c_{1}r^{-1}+c_{2}r^{-3}) \zeta\cos\theta_{\zeta}+c_{3}r^{-2}(\zeta^{2}-\frac{5}{2})$$+k_{\infty}r^{-2}[ \gamma_{1}c_{1}-2c_{3}r^{-1}(A(()\cos\theta_{\zeta}-\frac{1}{2}(c_{1}+3c_{2}r^{-2})\zeta^{2}B(()(1-3\cos^{2}\theta_{\zeta})]$$+k_{\infty}^{2}r^{-3}\{2c_{1}\zeta D_{1}(\zeta)\cos\theta_{\zeta}+3c_{3}r^{-1}\zeta^{2}F(\zeta)(1-3\cos^{2}\theta_{\zeta})$$-2[c_{1}(3\cos^{2}\theta_{(}-2)+3c_{2}r^{-2}(5\cos^{2}\theta_{\zeta}-3)]\zeta^{3}D_{2}(\zeta)\cos\theta_{\zeta}\}$$=$$\frac{\Phi_{s}}{u_{\infty}}$

,

(49)

$-2-c_{1}r^{-1}+c_{2}r^{-3}+k_{\infty}r^{-3}[c_{3}A(\zeta)-3c_{2}r^{-1}\zeta B(\zeta)\cos\theta_{\zeta}]$$+k_{\infty}^{2}r^{-3}\{-c_{1}D_{1}(\zeta)+6c_{3}r^{-1}(F(\zeta)\cos\theta_{\zeta}$(50)$- \frac{1}{2}[c_{1}(1-3\cos^{2}\theta_{(})+9c_{2}r^{-2}(1-5\cos^{2}\theta_{\zeta})](2D_{2}(()\},$

,

(51a)

$\frac{\omega}{u_{\infty}\cos\theta}$$=$$(\gamma_{1}c_{1}k_{\infty}-c_{3})r^{-2}$$\frac{u_{r}}{u_{\infty}\cos\theta}$$=$$1+c_{1}r^{-1}+c_{2}r^{-3}$$\frac{u_{\theta}}{u_{\infty}\sin\theta}$$=$$-1- \frac{1}{2}c_{1}r^{-1}+\frac{1}{2}c_{2}r^{-3}$$=$$c_{3}r^{-2}$$\frac{\tau}{u_{\infty}\cos\theta}$

,

(51b)

,

(51c)

,

(51d)

, and
where
are undetermined constants, and $A((), B(\zeta),$ $D_{1}((), D_{2}(\zeta)$ , andthe solutions of the following integral equations, which are tabulated in Ref. 40 [ $A(()$ andare first accurately obtained in Ref. 41] :$c_{1},$$c_{2}$$F(\zeta)$$c_{3}$$\mathcal{L}[\zeta;A(()]=-\zeta_{i}(\zeta^{2}-\frac{5}{2})$subsidiary condition:,$\int_{0}^{\infty}\zeta^{4}A(\zeta)\exp(-\zeta^{2})d\zeta=0$

are$B(\zeta)$

(52)

,

$\mathcal{L}[(\zeta_{i}\zeta_{j}-\frac{1}{3}\zeta^{2}\delta_{*j})B(\zeta)]=-2(\zeta;\zeta_{j}-\frac{1}{3}(2\delta_{*j}),$(53)$\mathcal{L}[(\zeta_{i}\zeta_{J}\cdot-\frac{1}{3}\zeta^{2}\delta_{ij})F(\zeta)]=(\zeta_{i}\zeta_{j}-\frac{1}{3}(2\delta_{ij})A((),$

(54)

76

$\mathcal{L}[(\zeta_{i}\delta_{jk}+\zeta_{j}\delta_{k*}+\zeta_{k}\delta_{*j})D_{1}(\zeta)+\zeta_{*}(j(kD_{2}(()]$(55)$=\gamma_{1}((i\delta_{jk}+\zeta_{j}\delta_{k*}+\zeta_{k}\delta_{ij})-(;(j(kB((),$

subsidiary condition: $\int_{0}^{\infty}[5\zeta^{4}D_{1}(\zeta)+\zeta^{6}D_{2}(\zeta)]\exp(-(2)d\zeta=0$ .

, andare introduced since the information on the side of
The undetermined constants
the sphere is not taken into account.
, Eqs. (49) and (50) can be used to evaluate the boundaryIfis taken larger than, andin Eqs. (51b),
conditions (41) and (42). That is, by connecting, andobtained in the step (iv) of(51c), and (51d) with
(say
, and
determinetheconstants
respectively,
weprocess
in
Sec.
IV
A
of
solutionthe, and). Letandin Eqs. (49) and (50) with these constants be denoted by
as
) at
(or
and (
, respectively. We takeand$c_{1},$$c_{3}$$c_{2}$$r_{A}$$r_{D}$$u_{r}/\cos\theta,$$(u_{r}/\cos\theta)_{(I)}^{(n)},$$\tau/\cos\theta$$u_{\theta}/\sin\theta$$(\tau/\cos\theta)_{(I)}^{(n)}$$(u_{\theta}/\sin\theta)_{(I)}^{(n)}$$c_{1},$$c_{1}^{(n)},$$c_{2}^{(n)}$$c_{3}^{(n)}$$c_{2}$$\Phi_{c}$$\Phi_{cA}^{(n)}$$\hat{\Phi}_{sA}^{(n)}$$c_{3}$$\Phi_{cA}^{(n)}$$\hat{\Phi}_{s}$$(r^{(I-1)}, \zeta^{(j)}, \theta_{\zeta}^{(k)})$$\hat{\Phi}_{sA}^{(n)}$$r^{(I)},$$(^{(j)}, \theta_{\zeta}^{(k)})$$G_{c(j,k)}^{(n+1)}$(or) in Eq. (41) [Eq. (42)] respectively. Thus, the boundary
) and
(or
conditions (41) and (42) for the finite difference method are computed.
$H_{s(j,k)}^{(n+1)}$$H_{c(j,k)}^{(n+1)}$$G_{s(j,k)}^{(n+1)}$

V.

RESULT AND DISCUSSIONA.

Macroscopic variables

is proportional, and have simple dependence on
The macroscopic variables[see Eqs. $(19a)-(19e)$]. The result of numericaltoand the others are proportional toin Fig. $3b$ and Table II,
vs is shown in Fig. $3a$ and Table I,computation
$3d$
$3c$
and Table IV. In these figures
in Fig.
and Table III, andin Fig.
(see Sec. III) are also shown. For
and tables, the two hmiting solutions for $k_{\infty}=0$ andcomparison, the corresponding results fdr the BKW equation are shown in Figs. 4a-4d.
The uniform flow is retarded by the sphere, whose effect is larger for smaller Knudsen numbers. One may think the effect of sphere extends further in the gas for larger Knudsen numbers
since the molecules reflected on the sphere proceed in the gas with less collisions for larger Knudsen numbers, but it is not so in a steady flow. In the free molecular flow, only the moleculesthat hit the sphere get the information of the existence of the sphere. That is, only the part$O(1/r^{2})$ of the molecules is affected by the sphere. When the Knudsen number is very small,on the other hand, the molecules reflected on the sphere soon collide with other molecules, givethem the information of the sphere, return to the sphere, and reflected there again. The collided
molecules soon collide with other molecules and so on. Successively the effect of the sphere
propagates into the gas and finally to all the molecules. As the Knudsen number increases,
the molecules travel longer distance before coUision and some of the molecules that return tothe sphere for smaller Knudsen numbers do not return to the sphere. Thus, the effect of thesphere is weaker for larger Knudsen numbers in the steady flow, although in a transient process
disturbances may propagate faster for larger Knudsen numbers.
The temperature is not uniform, and the nonuniformity is larger for larger Knudsen numbers. This is a typical effect of gas rarefaction called thermal polarization.13 Obviously in thecontinuum limit $(k_{\infty}=0)$ , the temperature is uniform in a slow flow where a quantity of thesecond order of the Mach number can be neglected. For $k_{\infty}=0.1$ (a small Knudsen number),
the variation of the temperature for the BKW model is the opposite sign to that for the
hard-sphere molecular gas, although both are very small (Figs. $3d$ and $4d$ ). To understand thethermal polarization, let us examine the velocity distribution function at forward stagnation$\omega,$$u_{r},$$u_{\theta}$$\theta:u_{\theta}$$\tau$$\cos\theta$$\sin\theta$$u_{r}/u_{\infty}\cos\theta$$u_{\theta}/u_{\infty}\sin\theta$$r$$\tau/u_{\infty}\cos\theta$$\omega/u_{\infty}\cos\theta$$k_{\infty}=\infty$$\tau$

77
point ( $r=1,$
for $U_{\infty}>0$ ) for the free molecular flow. The distribution of the molecules, andleaving the sphere is the half $(\zeta_{1}<0)$ of a stationary Maxwellian with temperaturethat of the incoming molecules is the part ( $(1>0)$ of the MaxweUian with temperature
andvelocity. Thus the width of the distribution at the stagnation point is wider by the orderof $U_{\infty}(>0)$ than that of the distribution with temperature, which means the temperatureby the order ofof the gas at the stagnation point is higher than.
$\theta=\pi$$\tau_{\infty}$$T_{\infty}$$U_{\infty}$$T_{\infty}$$T_{\infty}$B.

$U_{\infty}$

Drag on the sphere

The force acting on the sphere is obtained by integrating the momentum flux, which isequivalent to stress in the linearized problem, over any surface enclosing the sphere. Letbe$0,0)$ ] be the uniform flow. Then,the force on the sphere and ; [instead of a special$F_{i}$$U$$(U_{\infty},$$F_{*}\cdot=p_{\infty}L^{2}(2RT_{\infty})^{-1/2}U_{1}h_{D}(k_{\infty})$where$h_{D}(k_{\infty})$$(k_{\infty}=\infty)^{4}$, a function of$k_{\infty}$

,

(56)

, is given in Fig. 5 and Table V. For the free molecular flow,$h_{D}(\infty)=2\sqrt{\pi}(\pi+8)/3$

,

(57a)

and for small Knudsen numbers,$h_{D}(k_{\infty})=6\pi\gamma_{1}k_{\infty}(1+\kappa_{0}k_{\infty})$

,

(57b)

which is calculated from the asymptotic solution in Sec. III (For the BKW equation, anotheris $obtained^{10}$ ). In Fig. 5 various results by Millikan’s experiment or
higher order term inother methods are also shown for comparison. The result for the BKW equation is our recent
detailed computation by the method similar to the present analysis but the difference from theresults in Refs. 8, 12, and 14 is small (See Table VI). The way of comparing the results ofdifferent molecular models is not unique. The relation between the viscosity and the mean
is given by. Since depends on molecular models
free path[Eq. (27)], the result of comparison of a quantity (e.g.,) differs whether or
is taken as
the common independent variable. When is taken as the common independent variable,
for a hard-sphere molecular gas is related tofor the BKW model as$k_{\infty}$$\mu$$l_{\infty}$$\mu=(\sqrt{\pi}/2)\gamma_{1}p_{\infty}(2RT_{\infty})^{-1/2}l_{\infty}$$\gamma_{1}$$h_{D}$$\mu$$\ell_{\infty}$$k_{\infty}$$\mu$$k_{\infty}$$k_{\infty}(BKW)=1.270042k_{\infty}$(hard sphere).

(58)

In the comparison of the drag in Fig. 5, the conversion formula (58) for the BKW model and asimilar one for Millikan’s result are used.
Recently the drag forhas been reported in Ref. 24. The results seem to he nearly on
the same curve as ours (shghtly larger forand smaller forthan ours) in Fig. 5.
In Ref. 24 the interest was the global quantity, i.e., the drag, and the kernel of the collision
integral is expanded and truncated at the fourth term. The result in Ref. 24 is a numericalsolution of the model Boltzmann equation with the truncated collision kernel. In the numericalanalysis the discontinuity of the velocity distribution function is not taken into account. Bothsave a large amount of computation. Thus, the purpose and approach in Ref. 24 are different
from ours, where we tried to solve the original equation faithfully.
The drags for the BKW model computed by different methods are compared in Table VI.
The difference is not large as a whole. The variational result, though a very simple test function
was used, agrees very well with our result, where the discontinuity is taken into account. In
the integral equation approach in Ref. 12, the difficulty of the discontinuity is bypassed42, but
the temperature field is assumed to be uniform, which means that larger errors are expected$k_{\infty\sim}>1$$k_{\infty}\sim 1$$k_{\infty}\sim 10$

78(cf. Fig. $4d$ ). The, since the variation of the temperature is larger for larger
for larger)
and larger for moderate
approach in Ref. 24 shows the same tendency (smaler for large
as that for a hard-sphere molecular gas with the kernel truncation.
$k_{\infty}$$k_{\infty}$$k_{\infty}$$k_{\infty}$

Velocity distribution function

C.

andat several points in the gas are
The velocity distribution functions
andof the uniform flow atshown for $k_{\infty}=0.1,1$ , and 10 in Figs. 6-11. The,
12
infinity, common to all
are shown in Fig. for reference. The difference due to the Knudsen(the free molecular flow; Sec. III), two classes of molecules arenumber is clear. At
distinct. That is, at a point $(r)$ , the molecules with velocity $0\leq\theta_{\zeta}<\pi-Arc\sin(1/r)$ comedirectly from infinity, whose distribution is given by the part $[0\leq\theta_{\zeta}<\pi-Arc\sin(1/r)]$ of theuniform flow in Fig. 12; the molecules with the other velocities come directly from the sphere,
whose distribution is given in Fig. 13 [Eqs. (23) and (24)]. At $k_{\infty}=10$ (alarge Knudsen number)in Figs. 10 and 11, this feature is well preserved and only local correction is seen. In contrast to
the free molecular flow, however, the discontinuity of the velocity distribution function on thesphere decays in a long distance owing to molecular collisions. At $k_{\infty}=0.1$ (Figs. 6 and 7) the
behavior is quite different. The discontinuity of the velocity distribution function on the spheredecays in a short distance (Fig. $6aarrow 6barrow 6c$ ; Fig. $7aarrow 7barrow 7c$ ), and it transforms completely
into a distribution similar to that of the uniform flow in a fairly short distance (Fig. $6aarrow 6d$ ;Fig. $7aarrow 7c$ ). Further away from the sphere the change in the distribution is mainly in size
but not in shape (Fig. $6earrow 12a$ ; Fig. $7carrow 7darrow 12b$ ). At $k_{\infty}=1$ (an intermediate Knudsen
number), the behavior of the distribution retains more of the feature of $k_{\infty}=10$ . For $k_{\infty}=1$
and 10, the velocity distribution function is, more or less, of similar shape to that of the uniform
flow in a region where its discontinuity is vanishingly small. This shows that the flow is nearlyin equihibrium there. On the other hand, for $k_{\infty}=0.1$ , the decay of the discontinuity is so fastthat the velocity distribution function undergoes thorough transformation in a region where thediscontinuity is vanishingly $smaU$ . This kinetic transition region is the Knudsen layer, and theregion with discontinuity is the layer at the bottom of the Knudsen layer, discussed in Refs. 26
and 27.
In order to show the decay of the discontinuity with distance clearly,andplane at a given are shown in Figs. 14-17. The
inandexpressing the[cf. Eq. (18)].
uniform flow at infinity are given by
andThese sinusoidal surfaces are disturbed by the sphere, which imposes the condition (17). The
discontinuity of decays monotonicaUy with distance [cf. Eq. (20)]. This does not mean thatdecay monotonically because of the interaction
the discontinuities ofandterms between
andin Eqs. (11) and (12). The discontinuity is smaller and decays inshorter distance for smaller Knudsen numbers (Figs. 14-16). Figures 15 and 17 show that thediscontinuity decays more rapidly for smaUer molecular speed. This reflects the fact that thediscontinuity decays in the scale of the free path (but not the mean free path). At
and
10 (Figs. 15 and 16), the discontinuity extends to a point of the order of the mean free pathfrom the sphere, but at $k_{\infty}=0.1$ (Fig. 14), it extends only to a point much nearer than a pointof the order of the mean free path from the sphere. This is because the discontinuity is nearlyparallel near the sphere and decays in a distance of the order of the free path. This very thinlayer with discontinuity for $smaU$ Knudsen numbers corresponds to the layer at the bottomof the Knudsen layer. The figures for the BKW model corresponding to Figs. 14-17 are given
in Ref. 26. They are quite similar to the present results. For example, the figure for the BKW
model corresponding to Fig. 15 is given in Fig. 18.
$\Phi_{c}E/u_{\infty}$$\hat{\Phi}_{s}E/u_{\infty}$$\Phi_{c}E/u_{\infty}$$\hat{\Phi}_{s}E/u_{\infty}$$k_{\infty}$$k_{\infty}=\infty$$S$$\Phi_{c}E/u_{\infty}$$\Phi_{c}/u_{\infty}$$\zeta$$r\theta_{(}$$\Phi_{c}/u_{\infty}=-2(\cos\theta_{\zeta}$$\hat{\Phi}_{s}E/u_{\infty}$$\hat{\Phi}_{s}/u_{\infty}$$\hat{\Phi}_{s}/u_{\infty}=-2(\sin\theta_{\zeta}$$\phi$$\Phi_{c}$$\Phi_{c}$$\hat{\Phi}_{s}(=\zeta\sin\theta_{\zeta}\Phi_{s})$$\Phi_{s}$$k_{\infty}=1$

$S$

79
D.

Lattice data and accuracy test

The data of the lattice system that was used to obtain the results in Secs. V $A,$ , and aresummarized. The typical examples of the size of physical space
and the lattice are shown in
Table VII. The number of section is commonly 300 for all. The molecular speed is considered
in $0\leq(\leq\zeta_{D}$ with ($D=5$ , and the interval is divided into 24 sections, $i.e.,$$\zeta^{(i)}-=(D(i/J)^{2}$ with
$J=24$ . The region
is divided into uniform 176 sections, i.e., $K=88$ and $K=176$ .
In this lattice the molecular velocity spacehas 4425 lattice points, which means that 4425
terms are taken in the senies (43) and (44) The criterion of convergence of the sequencesandis that the variations of the macroscopic variables,, andin 10 iterations are less than.
It is inefficient, and also impossible in the present circumstances, to carry out extensive
direct accuracy tests by various lattice systems on the basis of the Boltzmann equation for hardsphere molecules since the computation is very large. The errors come from two different typeof sources: computation of the coUision integral and numerical integration of the differential
system. For the latter the BKW equation is a guide. We therefore investigate the same problem
by the same method on the basis of the BKW equation for various lattice systems, which can be
carried out by engineering workstations such as MIPS RS 3230 and HP 9000730. The typical
data showing the dependence of the results on lattice systems are given in Table VIII, where, andthe maximum differences of the macroscopic variables
by the use of different lattice systems over the whole flow field are shown. In the, the discontinuity of the velocity distribution function isandoriginal computation of
not taken into account since the error of the integrand hes only between two lattice points(cf. Sec. IV B). This error is tested for the BKW equation. The resultsinare also shown Table VIII. Ftom the data in our previous works (Refs. 34, 37, etc.) on the
Boltzmann equation for hard-sphere molecules, in addition to the above tests on BKW equation,
we have chosen the lattice system mentioned in the first paragraph of this subsection.
Now return to solutions of the Boltzmann equation for hard-sphere molecules. The mag, and,nitude ofand
at $\zeta=3.828$ is less than 6
for any
and therefore the size $\zeta_{D}(=5)$ adopted in the present analysis is sufficiently large. Two tests
for different lattice systems are carried out. In the test 1 (T1), the first 241
lattice points$(i=0-240)$ among 301 are reduced to 161, and in the test 2 (T2), the original 176 divisions ofare reduced to 132 divisions. The results of both computations of the macroscopic, andvariables
are compared with those on theoriginal lattice system. Their maximum differences over the whole flow field are shown in Table IX, where the maximum values of the corresponding variables over the whole flow field arealso given for reference. The lattice forseems to be sufficient. In the intermediate Knudsen
numbers, more
requiredtolattices are
get the result of the same accuracy. This feature isconsistent with that of the BKW equation.
The collision integral vanishes for Maxwellian distributions ( $\phi=1,$ , or ). For correand , therefore,
sponding$B$$C$$r_{D}$$k_{\infty}$$0\leq\theta_{\zeta}\leq\pi$$\zeta\theta_{(}$$\Phi_{c}^{(n)}$$\hat{\Phi}_{s}^{(n)}$$(\omega/u_{\infty}\cos\theta)^{(n)},$$(u_{r}/u_{\infty}\cos\theta)^{(n)}$$5\cross 10^{-7}$$(\tau/u_{\infty}\cos\theta)^{(n)}$$(u_{\theta}/u_{\infty}\sin\theta)^{(n)}$$\omega/u_{\infty}\cos\theta,$$u_{r}/u_{\infty}\cos\theta,$$u_{\theta}/u_{\infty}\sin\theta$$\tau/u_{\infty}\cos\theta$$\mathcal{L}_{1}$$\mathcal{L}_{2}$$\theta_{\zeta}^{(2\kappa)}\leq\theta_{\zeta}<\theta_{\zeta}^{(2\kappa+2)}$$\cross 10^{-7}$$\hat{\Phi}_{s}E/u_{\infty}$$\Phi_{c}E/u_{\infty}$$r,$$\theta_{\zeta}$$k_{\infty}$$r^{(i)}$$0\leq\theta_{\zeta}\leq\pi$$\omega/u_{\infty}\cos\theta,$$u_{r}/u_{\infty}\cos\theta,$$u_{\theta}/u_{\infty}\sin\theta$$\tau/u_{\infty}\cos\theta$$\theta_{\zeta}$$r^{(i)}$$\zeta_{i}$$\Phi_{c}$$\zeta_{j^{2}}$$\hat{\Phi}_{s}$$C_{(1,j,k)}=\nu_{(j)}\Phi_{c(i,j,k)}$$\hat{S}_{(i,j,k)}=\nu_{(j)}\hat{\Phi}_{s(i,j,k)}$

,

(59)

.

The maximum differences between the computed and exact values of $(CE)_{(*,j,k)}$ and, andfor $\Phi_{c}=1,$
are 3.3 $x10^{-5},1.4\cross 10^{-5},5.8x10^{-6}$ , and 5.2 $x10^{-5}$ ,$(\hat{S}E)_{(*,j,k)}$$\zeta\cos\theta_{\zeta},$$(^{2}$$\hat{\Phi}_{s}=(\sin\theta_{\zeta}$

80
respectively. Incidentally, the maximums of $\nu(()E, \nu(()(E\cos\theta_{\zeta}, \nu(()(2E$ , andare about 0.13, 0.06, 0.06, and 0.06 respectively.
The drag on the sphere can be computed by integrating the momentum flux over any surface
enclosing the sphere. In principle this serves a good accuracy test of computation. In practice it
is too severe in a three-dimensional problem in an infinite domain such as the present one to beapphed to an arbitrarily large control surface. For a large control surface, the drag is obtained
as a small quantity multiplied by a large area. Thus, the error in local variables is multiphed by
the factor , where is the characteristic linear dimension of the control surface. The resultsin Table. V are computed on the sphere $r=1$ . In the tests for the control spheres $r=r^{(i)}$ for all
between 1 and 4, the variation (max–min) ofis less than 0.15% of
for $k_{\infty}=0.2-10$$k_{\infty}=0.1$$k_{\infty}=0.05$, 3.1% for). Another test to confirm the result of
was made(0.56% forpoints (601)
for the BKW equation: The drag computed with the lattice system of double
remains constant with width 0.26% allowance for the control surfaces $r=r^{(*)}$ with $r^{(:)}=1-10$ .
The data for $r^{(i)}=1-4$ are essentially the same as those computed with the original system oflattice points.
301
A practical point that was introduced to improve the accuracy of numerical computation isadded here. The solution approaches the uniform flow at infinity fairly slowly, and therefore
for a wide range of $r(<r_{D})$ , the difference between the solution and its asymptotic solution isfairly small. We consider the difference :$\nu(()(E\sin\theta_{\zeta}$$r_{c}^{2}$$r_{c}$$r^{(i)}$$h_{D}$$h_{D}$$h_{D}$$r^{(i)}$$r^{(i)}$$\phi’$$\phi’=\phi-\phi_{a}$

,

(60)

$\phi_{a}/u_{\infty}=2(1+c_{1}r^{-1})\zeta_{1}+c_{1}r^{-1}\zeta_{\theta}\sin\theta$

where$\phi_{a}$,corresponds to the terms up to $O(1/r)$ in Eqs. (49) and (50). Then, from Eq. (1)(61)$(. \frac{\partial\phi’}{\partial x_{\dot{\iota}}}=\frac{1}{k_{\infty}}\mathcal{L}(\phi’)-(:\frac{\partial\phi_{a}}{\partial x_{i}},$is a local Maxwelhan. Over a wide range,
is fairly
whereterm disappears since, and correspondingly the error of.
smaller thanis fairly smaller than that of
Furthermore, by analyzing Eq. (61) instead of Eq. (1), we can avoid the error accumulated byfor large (cf. footnote 39
long range integration to form the slowly varying dominant partin Ref. 25).
The computation was carried out by HP 9000730 and MIPS RS 3230 at our laboratory (forthe BKW equation) and FACOM VP-2600 computer at the Data Processing Center of Kyoto
University (for the Boltzmann equation for hard-sphere molecules).
$\mathcal{L}(\phi_{a})$$|\phi’|$$\phi_{a}$$|\phi|$$\mathcal{L}(\phi’)$$\mathcal{L}(\phi)$$\phi_{a}$

VI.

$r$

CONCLUDING REMARKS

We have investigated the behavior of a rarefied gas around a spherical body set in a uniform
flow on the basis of the hnearized Boltzmann equation for hard-sphere molecules and the diffuse
reflection boundary condition. With the aid of the similarity solution, the original hnearized
Boltzmann equation with six independent variables is reduced to two simultaneous integrodifferential equations with three independent variables. These equations are solved numericallyby a finite-difference method. The computation of the collision integral, which takes the majorityof the computing time, is carried out efficiently by the numerical kernel method. The data of thekernel stored can be applied to other computations. Thus, we have obtained the comprehensive
and accurate information of the flow field, including the velocity distribution function with

81
discontinuity, and the force acting on the sphere for the whole range of the Knudsen number.
The layer at the bottom of the Knudsen layer is discussed by the behavior of the velocity
distribution function.
The velocity distribution function of the gas molecules has a discontinuity in the gas, which
is of considerable size and extends to a wide region for intermediate and large Knudsen numbers.
This is a typical feature of a gas around a convex body. Thus it introduces another difficulty inanalysis of rarefied gas flow problems. Some of the results of the present analysis are comparedwith those of the BKW model equation. In view of the result of comparison and the simphcityof the BKW equation, for which the analysis of the present problem can be carried out byengineering workstations such as MIPS RS 3230, it is a nice equation to be used in engineering
problems.
In the present analysis, the temperature of the sphere is as sumed to be equal to that of thegas at infinity. This is realized without any special device when the thermal conductivity of thesphere is much larger than that of the gas. In this case the temperature of the sphere is uniform
and the energy flow from the sphere vanishes since the sphere cannot be a source or sink ofenergy in a steady state. From the similarity solution, the latter condition is seen to be satisfied
by the present solution, where the temperature of the sphere is the same as that of the gas at
infinity. When the sphere is kept at a different temperature, a spherically symmet $ric$ solution
should be superimposed. This does not influence the velocity field and the drag on the sphere.
$S$

REFERENCES
1. R. A. Millikan, “The isolation of an ion, a precision measurement of its charge, and thecorrection of Stokes’s law,” Phys. Rev. 32, 349 (1911).
2. M. Knudsen and S. Weber, “Luftwiderstand gegen die langsame Bewegung kleiner Kugeln,“Ann. der Phys. 36, 981 (1911).
3. R. A. Millikan, “The general law of fall of a small spherical body through a gas, and its
bearing upon the nature of molecular reflection from surfaces,“ Phys. Rev. 2nd Ser. 22, 1(1923).
4. P. S. Epstein, “On the resistance experienced by spheres in their motion through gases,”
Phys. Rev. 2nd Ser. 23, 710 (1924).
5. R. Goldberg, “The flow of a rarefied perfect gas past a fixed spherical obstacle,” Ph. D.
Thesis, New York Univ. (1954).
6. A. B. Basset, Hydrodynamics (Dover, New York, 1961), Vol. II, p. 271.
7. D. R. Wilhs, “Sphere drag at high Knudsen number and low Mach number,” Phys. Fluids
9, 2522 (1966).
8. C. Cercignani, C. D. Pagani, and P. Bassanini, “Flow of a rarefied gas past an axisymmetricbody. II. Case of a sphere,” Phys. Fluids 11, 1399 (1968).
9. W. Phillips, “Drag force on a small sphere moving through a gas,” Phys. Fluids 18, 1089
(1975).
10. Y. Sone and K. Aoki, “Forces on a spherical particle in a shghtly rarefied gas,” in Rarefied
Gas Dynamics, edited by J. L. Potter (AIAA, New York, 1977), Part 1, p. 417.
11. Y. Sone and K. Aoki, “Slightly rarefied gas flow over a specularly reflecting body,“ Phys.
Fluids 20, 571 (1977).
12. K. C. Lea and S. K. Loyalka, “Motion of a sphere in a rarefied gas,” Phys. Fluids 25, 1550
(1982).
13. S. P. Bakanov, V. V. Vysotskij, B. J. Derjaguin, and V. I. Roldughin, “Thermal polarizationof bodies in the rarefied gas flow,” J. Non-Equihib. Thermodyn. 8, 75 (1983).

82
14. W. S. Law and S. K. Loyalka, “Motion of a sphere in a rarefied gas. II. Role of temperature
variation in the Knudsen layer,” Phys. Fluids 29, 3886 (1986).
15. K. Aoki and Y. Sone, “Temperature field induced around a sphere in a uniform flow of ararefied gas,” Phys. Fluids 30, 2286 (1987).
16. S. A. Beresnev, V. G. Chernyak, and G. A. Fomyagin, “Motion of a spherical particle in a
rarefied gas. Part 2. Drag and thermal polarization,” J. Fluid Mech. 219, 405 (1990).
17. C. Cercignani and C. D. Pagani, “Flow of a rarefied gas past an axisymmetric body. I.
General remarks,” Phys. Fluids 11, 1395 (1968).
18. P. L. Bhatnagar, E. P. Gross, and M. Krook, “A model for collision processes in gases, I,”
Phys. Rev. 94, 511 (1954).
19. P. Welander, “On the temperature jump in a rarefied gas,” Ark. Fys. 7, 507 (1954).
20. In this paragraph, we consider the case where the temperature of the sphere is the same asthat of the uniform flow (or the particle has very high thermal conductivity) and the speedof uniform flow is so low that the quantities of the order of the Mach number squared are
negligible.
21. Y. Sone, “Asymptotic theory of flow of rarefied gas over a smooth boundary I,” in Rarefied
Gas Dynamics, edited by L. Trilling and H. Y. Wachman (Academic, New York, 1969), Vol.
1, p. 243.
22. Y. Sone, “Asymptotic theory of flow of rarefied gas over a smooth boundary II,” in Rarefied
Gas Dynamics, edited by D. Dini (Editrice Tecnico Scientifica, Pisa, 1971), Vol. 2, p. 737.
23. S. Takata, Y. Sone, and K. Aoki, “Drag and thermal force on a spherical particle in a rarefied
gas- Numerical analysis for
Knudsen numbers-,” J. Vac. Soc. Jpn. 35, 143 (1992) (inJapanese).
24. S. K. Loyalka, .“Motion of a sphere in a gas: Numerical solution of the hnearized Boltzmann
equation,” Phys. Fluids A 4, 1049 (1992).
25. H. Sugimoto and Y. Sone, “Numerical analysis of steady flows of a gas evaporating from its
cyhndrical condensed phase on the basis of kinetic theory,” Phys. Fluids A 4, 419 (1992).
26. Y. Sone and S. Takata, “Discontinuity of the velocity distribution function in a rarefied gas
around a convex body and the S layer at the bottom of the Knudsen layer,” Transp. Theory
Stat. Phys. 21, 501 (1992).
27. Y. Sone, “New kind of boundary layer over a convex sohd boundary in a rarefied gas,” Phys.
Fluids 16, 1422 (1973).
28. H. Grad, “Asymptotic theory of the Boltzmann equation, II,” in Rarefied Gas Dynamics,edited by J. A. Laurmann (Academic, New York, 1963), Vol. 1, p. 26.
29. C. Cercignani, The Boltzmann Equation and Its Applications (Springer, Berlin, 1987), Chap.
IV, Sec. 5.
30. Y. Sone and K. Aoki, “Similarity solution of the linearized Boltzmann equation with application to thermophoresis of a spherical particle,” J. Mec. Theor. Appl. 2, 3 (1983).
31. Y. Sone, “Lecture Notes on Molecular Gas Dynamics,” (Department of Aeronautical Engineering, Kyoto University, Kyoto, 1990) Chap. III (in Japanese).
32. Y. Sone and K. Aoki, “Steady gas flows past bodies at small knudsen numbers-Boltzmann
and hydrodynamic systems-,” Transp. Theory Stat. Phys. 16, 189 (1987); “Asymptotic
theory of slightly rarefied gas flow and force on a closed body,” Mem. Fac. Eng. Kyoto Univ.
49, 237 (1987).
33. Y. Sone, “Asymptotic theory of a steady flow of a rarefied gas past bodies for small knudsen
numbers,” in Advances in Kinetic Theory and Continuum Mechanics, edited by R. Gatignol
and Soubbaramayer (Springer, Berhn, 1991), p. 19.
34. T. Ohwada, Y. Sone, and K. Aoki, “Numerical analysis of the shear and thermal creep flows
of a rarefied gas over a plane wall on the basis of the hnearized Boltzmann equation forau83
hard-sphere molecules,” Phys. Fluids A 1, 1588 (1989).
35. Y. Sone and Y. Onishi, “Kinetic theory of evaporation and condensation–Hydrodynamic
equation and shp boundary condition-,” J. Phys. Soc. Jpn. 44, 1981 (1978).
.
andshould be taken as the limiting values for
36. This means that$\Phi^{(n)}$$\hat{\Phi}^{(n)}$$40,j,\hat{K})$$\theta_{\zeta}arrow\pi/2-0$$s(0,j,\hat{K})$The other hmiting values are given by Eq. (17)37. Y. Sone, T. Ohwada, and K. Aoki, “TemperatUre jump and Knudsen layer in a rarefied gas
over a plane wall: Numerical analysis of the hnearized Boltzmann equation for hard-spheremolecules,” Phys. Fluids A 1, 363 (1989).
.
corresponds toin Ref. 37 with replacement ((38.
can be carried out analytically. After a similar
and
39. The integration with respect to intransformation in Ref. 34 and detailed classification of the cases depending on the relative
position of the singularity to the ring, the integrations are carried out by Gauss-Legendre
formula.
40. T. Ohwada and Y. Sone, “Analysis of thermal stress slip flow and negative thermophoresis
using the Boltzmann equation for hard-sphere molecules,” Eur. J. Mech., B/Fluids, 11, 389(1992).
41. C. L. Pekeris and Z. Alterman, “Solution of the Boltzmann-Hilbert integral equation II. thecoefficients of viscosity and heat conduction,” Proc. Natl. Acad. Sci. 43, 998 (1957).
42. For the BKW equation with the diffuse reflection condition, the velocity distribution function, and , which are continuous, arecan be ehminated, and the integral equations for.
derived. They are analyzed in Ref. 12.
, andin Eqs. $(14a)-(15b)$ have a symmetric property:
43. The operators
etc. For a uniform lattice system, thereandfore, $K_{j^{c}u_{m}}=K_{j^{c}(K-k)\ell(K-m)}$ , and $K_{j^{s}k\ell m}=K_{j^{s}(K-k)t(K-m)}$ . Thus only$\Psi_{\ell m}(\zeta_{1}, \zeta_{r})$$\Psi_{\ell m}(\zeta, \theta_{\zeta})$$arrow\zeta_{1}, \theta_{\zeta}arrow\zeta_{r})$$\overline{\psi}$$\mathcal{L}_{2}^{c}$$\mathcal{L}_{2}^{s}$$\omega,$$C_{1}^{c},$$C_{2}^{c},$$C_{1}^{s}$$u_{i}$$\tau$$C_{2}^{s}$$\mathcal{L}_{1}^{c}[\Phi(r, \xi, \theta_{\xi})]_{(r,\zeta,\theta_{(})}=\mathcal{L}_{1}^{c}[\Phi(r, \xi, (\pi-\theta_{\xi}))]_{(r,\zeta,\pi-\theta_{C})}$$\theta_{(}$$K_{j^{c}k\ell m}$$(0\leq m\leq K, 0\leq$

k $\leq K, 0\leq j,$l$\leq$J), 9,845,625 each, need to be prepared.

$K_{j^{s}k\ell m}$

84vs r) for various
TABLE I. The distribution of the radial flow velocity (are not common for differen.
accuracy from those at the lattice points,$u_{r}/u_{\infty}\cos\theta$$wl\iota ich$$k_{\infty}$. The data are interpolated with sufficient$tk_{\infty}$$\overline{\overline{\frac{\backslash rk_{\infty}\frac{u_{r}/u_{\infty}\cos\theta}{0005Q1020.406124610\infty}}{100000000000000000000\acute{0}0000000000000000000000000.00000000000000}}}$

101
1.1
1.2
13
14
15
1.7
2
2.5
3
4
5
7
10
15
20
25
50
100
200

500

0001
0.01200394
0737
0.11080.14810.21940.31250.43200.51856328
0.70400.78720.85050.90010.92510.94000.97000.98500.99259970$0$0015 0.0025
0245
0347
0600 0.0779
0.0994 0.1224. 16580.13930.1782 0.2068
2502 0.2808
0.3420 0.3725
4577 0.4856
5407 0.5656
6498
6696
7177 0.7337
7971 0.8083
8575 0.8651
9048 0.9096
9286 0.9321
9428 0.9456
9714 0.9728
9857 0.9864
9929 0.9932

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

9971

$0$

0040
0500

$0$

$0$

$0$

0060
0702
0.14170.2065.I050

0.15780.20680.25150.32900.42110.53040.60560.70160.76010.8276
8785
9185
9386
0.95089753
9876
0.9938
9975
9973
$0$

$0$

$0$

$0$

$0$

$0$

0086
0973
0.19140.27340.34354034
0.49850.59880.70217641
8339
0.87170.91150.93930099
0.11132170
0.30800106
0831
0.11900.16512314
0.23800.32730.30140.38484078
3562
0.44954751
5507 0.5797
4454
0.54240.6547 0.6858
0.64660.7578 0.7888
0.71228169 0.8464
0.8796
0.78989057
0.83420.91139343
0.9423 0.9606
8831
9187
9621
9761
9458 0.9600 0.9757 0.9857
0.9593 0.9700
9821 0.9898
9673 0.9760 0.9857
9920
9833 0.9878
9929
9961
9981
0.9915 0.9937 0.9963
9957 0.9968 0.9981
9990
9992
99830.9987
9996

$0$

$0$

$0$

$0$

0.26400.3146
0.39880.49320.5987

0072

$0$

$0$

0.6678

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

0.75250.80258591
0.90120.9338

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

0.9501$0$9599
0.9796
0.98970.9948

$0$

$0$

$0$

$0$

$0$

$0$

$0$

$0$

0.9979$0$$0$0.01090.12150.23600.33350.4152
0.48350.58930.69610.79924212
4902
0.59670.70390.80700.8564

$0$

$0$

$0$

4291
4990

0.60660.71450.81740.8736$0$0.9892
0.99250.99430.9974
0.99870.9993

$0$

$0$

$0$

9997

$u_{\theta}/u_{\infty}\sin\theta$

$0$

$0$

$0$

$0$

$w1\iota ic|_{l}$0.01130.12630.24480.34530.86380.9210 0.9294
0.9478 0.9551
9715 0.9773
0.98459889
0.9919 0.9951
0.9946 0.9973
0.9961 0.9983
0.99839996
0.9992 0.9999
0.9996 10000
999910000

9144
0.94200.96689809$0$vs r) for various
TABLE II. The distribution of the component of the flow velocity (are not common for different.
sufficient accuracy from those at the lattice points,$\theta$

0111
0.1236
0.23980.3386$0$$0$$k_{\infty}$. The data are interpolated with

$k_{\infty}$$\overline{\overline{\frac{\backslash rk_{\infty}\frac{u_{\theta}/u_{\infty}\sin\theta}{00050.1020406124610\infty}}{100000-00598-0.1078-0.1777-02662-03206-0.3805-04379-04696-04801-04884-0.5000}}}$OI$-00148$$-00854$$-0.1403$$-02216$$-03272$$-0.3929$$-04658$$-0.5362$$-05751$$-0.5882$$-0.5984$$-0$ 6128

1.1
12$-0.1304$$-0.1996$$-02598$$-0.3512$$-04748$$-0.5539$$-0.6425$$-0.7285$$-0.7763$$-07922$$-08047$$-0$ 8219$-0.2303$$-0.2S90$$-0.344S$$-04317$$-05531$$-0.6322$$-0.7215$$-0$ .SOS4$-0$ .S566$-0$ .S725$-0$ .SS50$-0.9021$$13$$-0.3093$$-0.3592$$-04095$$-04895$$-0.6042$$-0.6804$$-0.7668$$-0.8512$$-0.8979$$-0.9133$$-09253$$-0$ 9417

14
15
17

$-0.3732$$-04163$$-04616$$-05345$$-06418$$-0.7141$$-07965$$-0.8774$$-0.9222$$-0.9370$$-09484$$-0$ 9639$-0.4259$$-0.4636$$-05047$$-05712$$-06710$$-0.7393$$-0.8177$$-0.8948$$-0.9376$$-0.9517$$-09625$$-0$ 9772$-0.5079$$-0.5378$$-0.5721$$-06280$$-07146$$-07755$$-0.8460$$-0.9161$$-0.9551$$-0.9679$$-09777$$-0.9908$$2$$-05938$$-06163$$-0.6437$$-06883$$-07595$$-0.8111$$-0.8717$$-09330$$-0.9672$$-0.9784$$-09870$$-0$ 9983

25$-0.6840$$-06999$$-07203$$-07534$$-08075$$-08482$$-0.8967$$-09470$$-0.9756$$-09849$$-09920$

$-1$

3
4
5

$-0.7407$$-0.7531$$-0.7692$$-07954$$-08389$$-08723$$-0.9124$$-0.9549$$-0.9795$$-0.9876$$-09936$$-1$$-0$ .SOS6$-0$ .S174$-0.8284$$-0$ .S470$-08782$$-0.9027$$-0.9321$$-0.9645$$-0.9837$$-0.9902$$-0.9949$

0011
0014
$-1$ 0009$-0.8480$$-0$ .S550$-0$ .S631$-0$ .S775$-09020$$-0.9213$$-0.9444$$-0.9704$$-0.9863$$-0.9916$$-09956$$-1.0006$$7$$-0.8921$$-08972$$-09023$$-09122$$-09294$$-0.9431$$-0.9592$$-09776$$-0.9894$$-0.9935$$-09965$$-1$

10$-09248$$-0.9283$$-0.9315$$-09382$$-09501$$-0.9597$$-0.9708$$-09836$$-0.9920$$-0.9950$$-0.9973$$-1.0001$

$1$

0003

$15$

$-09499$$-0.9523$$-09544$$-09586$$-09664$$-0.9728$$-0.9802$$-09887$$-0.9943$$-09964$$-09980$$-1.0000$$20$$-0.9625$$-0.9642$$-0.9659$$-09690$$-09745$$-0.9794$$-0.9850$$-0.99I3$$-0.9956$$-0.9972$$-09984$$-1.0000$$25$$-0.9700$$-09714$$-0.9727$$-09752$$-09794$$-09833$$-0.9879$$-0.9930$$-0.9964$$-0.9977$$-09987$-I.OOOO

50$-0.9850$$-09857$$-0.9864$$-09876$$-09895$$-0.9914$$-0.9937$$-0.9964$$-09981$$-0.9988$$-09993$$-1.0000$$100$$-0.9925$$-09929$$-0.9932$$-0.9938$$-09948$$-0.9956$$-0.9967$$-09981$$-0.9990$$-09994$$-0.9996$$-1.0000$$200$$-09962$$-0.9964$$-0.9966$$-09969$$-09974$$-09978$$-09983$$-0.9990$$-0.9995$$-0.9997$$-09998$

$-1$

0000

500

$-099S5$$-09986$$-09986$$-09988$$-09990$$-09991$$-09993$$-0.9996$$-0.9998$$-0.9999$$-09999$

$-1$

0000

85
vs r) for various
TABLE III. The distribution of the density ((.
those at the lattice point , which are not common for different$v/u_{\infty}\cos\theta$$k_{\infty}$. The data are interpolated with sufficient accuracy from$k_{\infty}$$s$$\overline{\overline{\frac{r\backslash k_{\infty}\frac{(v/u_{\infty}\cos\theta}{0050.10.2040.6124610\infty}}{I-00925-0.1842-03522-0.6128-0.7897-09968-1.2055-1.3263-1.3683-1.4016-1.4504}}}$

1.01
$11$

$12$

$-00886-0.1740-0.3272-0.5622-0.7214-0.9076-10949-1.2034-1.2410-1.2708-1.3967$$-0.0730-0.1396-02543-0.4266-05432-0.6797-0.8170-0.8964-09239-0.9457-12014$$-0.0613-0.1158-0.2076-0.3435-0.4354-0.5429-0.6513-0.7137-0.7356-0.7526-1.0501$$1.3$$-00523-00451-00845-00982-01490-0..1743-0.\cdot 2855-02421-03048-0.3605-0.4484-0.3784-0.\cdot 5371-04523-0.\cdot 5881-04951-05101-0.6061-05218-0.\cdot 6200-0.9311-0.8345$
$1.4$

1.5

$-0.0394-0.0735-0.1291-0.2085-0.2617-0.3240-0.3872-0.4235-0.4363-0.4461-0.7545$

$1.7$

2
2.5
3

$-00100-00183-00318-00499-0.0616-0.0749-0.0883-0.0964-00988-0.1010-0.2897$$4$$-00056$$-00103$$-00178$$-0.0279$$-00344$$-00417$$-00489$$-00531$$-00545$$-0.0556$$-0.2003$$5$$-00036$$-00066$$-00114$$-00178$$-00219$$-00266$$-0031I$$-00337$$-0.0345$$-0.0351$$-0.1520$$7$$-0.0019-00033-00058-00091-00111-00135-0.0158-00170-00174-0.0177-0.1017$$10$$-00009$

$15$

$20$
$25$
$50$

$100$

$-00028$$-00016$$-00044$$-00054$$-00066$$-00077$$-0.0083$$-00085$$-00086$$-0.0676$$-00004-00007-0.0012-0.0019-0.0024-00029-00034-0.0037-0.0037-0.0038-0.0432$$-00002-00004-0.0007-0.0011-00013-00016-00019-00021$ $-0.0021$ $-0.0022-0.0317$$-00001$$-00004$$-00003$$-0.0007$$-00008$$-00010$$-0.0012$$-00013$$-00013$$-00014$$-0.0250$$-00000-00001-00001-0.0002-00002-0.0002-00003-00003-0.0003-00003-0.0122$$-00000-00000-0.0000-00000-00001-00001-00001$ $-0.0001-0.0001-0.0060$$\ovalbox{\tt\small REJECT} 200-00000-0.0000-00000-00000-0.0000-0.0030500-00000-0.0000-00000-00012$

TABLE IV. The distribution of he temperature (vs r) for various
from those at the lattice poilrts, which are not common for different.
$\tau/u_{\infty}\cos\theta$$t$$k_{\infty}$. The data are interpolated with sufficient accuracy$k_{\infty}$$\ovalbox{\tt\small REJECT}\backslash _{rk_{\infty^{005}}..20406124610}^{\tau/u.\cos\theta}1-00010-0^{0}0^{1}058-0^{0}0236-00620-00903-0^{\infty}1233-0.1551-0.1724-0.1780-0.1822-0^{\infty_{1881}}\ovalbox{\tt\small REJECT}$$1.01$

1.1

1.2
13
$1.4$

$-00001-0.0019-00098-0.0299-00462-00663-00870-00990-0.1032-0.1065-0.1128$$-00001$$-00016-0.0084-0.0255-00394-0.0567-0.0746-00851-0.osss-00917-0.0960$

$1.5$

17
2
2.5

$-00000$$-00006$$-00026$$-0.0077$$-00118$$-0$ OI70$-0.0226$$-00261$$-0.0274$$-00284$$-0.0301$$3$$-00001$$-00004$$-00019$$-00053$$-00082$$-00117$$-00156$$-0.0180$$-0.0189$$-0.0197$$-0.0209$$4$$-00001$$-00003$$-0$ OOI I$-0.0030$$-00046$$-00065$$-00086$$-0.0100$$-00105$$-0$ OI10$-0.0118$57$-00001-00001$$-00002-00001-00001$$-00007-00004-00002$$-0.0020-00010-00005$$-00008-00016-00030$$-00042-00021-00011$$-00055-00027-00013$$-0.0064-00032-000I5$$-00067-00034-00016$$-0.0070-00035-00017$$-00019-00038-0...0075$

$10$

15

$-00001$$-00002$$-00004$$-00005$$-00006$$-00007$$-00007$$-00007$$-0.0008$$20$$-00001$$-00001$$-00002$$-00003$$-00003$$-0.0004$$-00004$$-0.0004$$-0.0005$$25$$-00000$$-0.0001$$-00001$$-00002$$-00002$$-00002$$-0.0003$$-00003$$-0.0003$$-0.0000$$-00000$$-00000$$-00001$$-00001$$-00001$$-00001$$-0.0001$$-00000$$-00000$$-00000-00000-00000$$-00000-00000-00000$$-00000-00000-00000$$-0...0000-00000-00000$$50$$200500100$$-00000-00000$

86
TABLE V. The drag on the sphere:$h_{D}$

vs$k_{\infty}$[cf. Eq. \langle 56)].

$\overline{\overline{\frac{k_{\infty}0.050.10.20.40.6124610\infty}{h_{D}(k_{\infty})110912.11683.8110622927.79519562511277212.233312.555712.807113.1653}}}$

TABLE VI. Comparison ot$h_{D}$

in the drag formula (56) for the BKW model. The conversion (58) of$k_{\infty}$

is not made here.

$\overline{\overline{\frac{k_{\infty}\frac{h_{D}(k_{\infty})/h_{D}(\infty)}{Ref.8Ref.12Ref.24Presentmethod}}{10..64040..64020..647606403}100952909476095080.9540}}$7900
0.88642
4$0$

TABLE VII.

0.78990.88700.79330.88590.79070.8879

Examples ot lattice of(Zl,Rl, O)(Zl,Rl, O)(Zl,R2, O)(Z2,Rl, O)$(Z2,R1, O)$$(Z2,R1, \cross)$$r$

3 $x10^{-4}$2
3 $xI0^{-5}$Z2:$\zeta_{D}=5,$

6 $x10^{-4}$
4 $x10^{-4}$2 $x10^{-4}$1 $x10^{-3}$3 $x10^{-4}$2 $x10^{-4}$$x10^{-\ell}$Upper bounds $oi$ the differences ofwhole lat tice points.
bZl; $(D=6,$ $J=24,$ $K=300$“.

$\omega/u_{\infty}\cos\theta,$$u_{r}/u_{\infty}\cos\theta,$$u_{\theta}/u_{\infty}\sin\theta$, and$\tau/u_{\infty}\cos\theta$

for the systems A and$B$

over the$J=24,$ $K=176$

Rl: $I=S00,$ $r_{D}=13.77\{k_{\infty}=0.1$ ), 199.26 $(k_{\infty}=1)$ , 928.72 $(k_{\infty}=10)$), 928.72 $(k_{\infty}=10)$
R2: $I=600,$ $r_{D}=13.77(k_{\infty}=0.1)$ , 199.26 \langle
R3: $I=320,$ $r_{D}=20.38(k_{\infty}=0.1)$ , 304.52 $(k_{\infty}=1)$ , 158987 $(k_{\infty}=10)$$O$ ; In the computation oi the macroscopic variables in theterm of collision, the discontinuity of: not taken into account.
$k_{\infty}=1$

$g$

$n$

$\phi$

is taken into account.

$x$

TABLE IX. Dependence of the macroscopic variables on the lattice systems ior the hard-sphere gas.

Upper bounds of the differences of a macroscopic variable for the original lattice system (cf. first paragraph of Sec. V D)
lattice points (cf third paragraph of Sec. V D) over the whole lattice points.
and for the system with reducedb Upper bounds of the differences of a macroscopic variable for the original lattice system and for the system with reduced
lattice points (cf. third paragraph of Sec. V D) over the whole lattice points.
c Maximum of the absolute value.
“$r^{(i)}$$\theta_{\zeta}^{(k)}$87$(u_{\infty}, 0,0)$$x_{3}$

FIG. 1. Geometry and coordinate systems [Sec. II A and Eq. (8)].

FIG. 2. Discontinuity of the velocity distribution function. At
he point $x;$ , tlle velocity distribution function‘is discontinuous on the aded
ilt space.
$l$$s1_{1}$$COI\iota e$$\dot{\zeta}_{j}$vsvs , and (d)vsFIG. 3. The profiles of the macroscopic variables (llard sphere). (a)vs . Here, –indicates the presen numerical result, $—-$ -the asymptotic solution $(26a)-(26d)$ with $k_{\infty}=0$ ,. The indicates the value at $r=1$ .
andthe free molecular ftow$u_{r}/u_{\infty}\cos\theta$$\tau/u_{\infty}\cos\theta$

—-

$t$

$r$

$(k_{\infty}=\infty)$$x$$f(b)u_{\theta}/u_{\infty}\sin\theta$$r,$$(c)\omega/u_{\infty}\cos\theta$

$r$

88

vs
vs
vs , and (d)FIG. 4. The profiles of the macroscopic variables (BKW). (a)vs . Here, –indicates the present numerical result, $—-$ -the asymptotic solution $(26a)-(26d)$ with $k_{\infty}=0$ ,. The indicates the value at $r=1$ . The conversion (58) ofand —-the free molecular flow
is not made here.
$u_{r}/u_{\infty}\cos\theta$$\tau/u_{\infty}\cos\theta$$r,$$(b)u_{\theta}/u_{\infty}\sin\theta$$r,$$(c)\omega/u_{\infty}\cos\theta$

$r$

$r$

$(k_{\infty}=\infty)$$x$$k_{\infty}$.

$=p_{\infty}L^{2}(2RT_{\infty})^{-1/2}U,h_{D}$ [Eq. (56)].
Here,
indicates the. The drag F. is given byvsFIG. 5. The drag on the sphere.
present result for hard-spbere molecules, $+the$ present result for tbe BKW model, $Ot$ he result by kernel expansion method[Eq. $(57b)$ ], $—-$ -the(Ref. 24), —–the empirical formula by Millikan (Ref. 3), –the asymptotic result for smallof Eq. $(57b)$ ], and —the result for the free molecular flow [Eq. $(57a)$ ].
tokes drag [the first-order term in$h_{D}$$F$$k_{\infty}$$\bullet$$k_{\infty}$$S$$k_{\infty}$

89

(a)

$r=1$

(c)

$r=1200$(b)$r=1010$

FIG. 6. The velocity distribution function
at five points
in he gas for $k_{\infty}=0.1$ . (a) $r=1,$ \langle ) $r=1.010$ ,(c) $r=1200,$ $(d)r=1.496$ , and (e) $r=1.998$ . The
surface ofis shown as a function ofand(by lines $\zeta=const$ and. The vertical
stripes sltow the discontinuity. Note the difference oi
the ordinate in. $6e$ .
$\Phi_{c}E$

$b$

$t$

$\Phi_{c}E$$\zeta\cos\theta_{\zeta}$$\theta_{\zeta}=const$$\sin\theta_{\zeta}$$I\cdot ig$(e)$r=199S$

at four points in the gas for $k_{\infty}=0.1$ . $(a)r=1,$ $(b)r=1010,$ $(c)r=1.100$ , and (d)FIG. 7. The velocity distribution function
const and. The verticalby lines
and ($r=1496$ . The surface of , is sllown as a function of (tripes show the discontinuity. Note the difference of the ordinate in Fig. $7d$ .
$\hat{\Phi}_{s}E$$\hat{\Phi}$

$s$

$E$

$\cos\theta_{\zeta}$$\sin\theta_{\zeta}$$\zeta=$$\theta_{\zeta}=const$

90

(a)

$r=1$

(c)

$r\Rightarrow 149S$(b)$r=1100$(d)$r=2.014$

FIG. 8. The velocity distribution function
at five points
in the gas ior. (a) $r=1,$ $(b)r=1.100,$ $(c)$$r=1.49S,$ $(d)r=2014$ , and (e) $r=6$ .IO1. (See thecaption oi Fig. 6.)$\Phi_{c}E$$k_{\infty}=1$(e)$r=6101$

FIG. 9. The velocity distribution function , at four points in the gas for$r=6.101$ . (See the caption of Fig. 7.)$\hat{\Phi}$$E$$k_{\infty}=1$. (a) $r=1,$ $(b)r=1.100,$ $(c)r=2.014$ , and (d)

91

(a)

$r=1$

(c)

$r=6018$(b)$r=1990$(d)$r=1104$

FIG. 10. The velocity distribution function
at five points
in the gas for $k_{\infty}=10$ . $(a)r=1,$ $(b)f=1990,$ $(c)$$r=6018,$ $(d)r=11.04$ , and (e) $r=20.S6$ . (See thecaption of Fig. 6)$\Phi_{c}E$(e)$r=20S6$

FIG. 11. The velocity distribution function , at four points in the gas for $k_{\infty}=10$ . $(a)r=1,$ $(b)r=1990,$ $(c)r=1104$ , and (d)$r=2086$ . (See he caption of Fig. 7.)$\hat{\Phi}$

$t$

$E$

92

(a)

FIG. 12. The velocity distribution functions$\Phi_{c}E$and$\hat{\Phi}$,$E$

in the uniform flow at infinity. (a)$\Phi_{c}E$and (b)$\hat{\Phi}_{s}$E. (common to all$k_{\infty}$)

and , at $\zeta=0.556$ for $k_{\infty}=0.1$ . (a)and (b) , . The surfaces oi
FIG. 14. The velocity distribution functions
and, are shown as functions of andby lines $r=$ const andconst. The vertical lines show the discontinuity. Thelines behind tber parts are shown by dashed lines.
$\hat{\Phi}$$\Phi_{c}E$$\hat{\Phi}$$E$$i_{1}\iota vlsible$

$r$

$E$

$\hat{\Phi}$$\Phi_{c}E$$E$$\Phi_{c}E$$\theta_{(}=$$\theta_{(}$$0$

FIG. 15 The velocity distribution functions$\Phi_{c}E$and$\hat{\Phi}_{s}E$at ( $=0.556$ for $k_{\infty}=1$ . $(a)\Phi_{c}E$ and (b)$\hat{\Phi}$, E. (See the caption of Fig. 14.)

93

FIG. 16. The velocity distribution functions

FIG. 17. The velocity distribution functions$\Phi_{c}E$$\Phi_{c}E$

and

and

$\hat{\Phi}$,$E$$\hat{\Phi}_{s}E$at ( $=0.556$ for $k_{\infty}=10$ . $(a)\Phi_{c}E$ and (b)at ( $=0.139$ for $k_{\infty}=1$ . $(a)\Phi_{c}E$ and (b)FIG. 18. The distribution functions
andfor the BKW equation at ( $=0667$ for
sphere). (a)and (b)E. (See the caption of Fig. 14.) The conversion (58) ofof the value (from that in Fig. 15.
$\Phi_{c}E$$\Phi_{c}E$$\hat{\Phi}_{s}$$\hat{\Phi}_{s}E$, E. (See the caption of Fig. 14.)$\hat{\Phi}_{s}$E. (See the caption of Fig. 14.)

, which corresponds to Fig. 15 (hardis not made here. Note the difference$k_{\infty}=1$$k_{\infty}$$\hat{\Phi}$