Motion of a plate moving in a free-molecular gas

自由分子気体中の物体の一定外力による運動
Motion of a plate moving in a free-molecular gas under a constant external force
機械理工学専攻 分子流体力学研究室 修士１回生 辻 徹郎
A thin plate accelerated or decelerated in a free-molecular gas at rest by a constant
external force is considered. The force is in the direction perpendicular to the plate. In
this situation, the plate velocity approaches its final constant velocity as time goes on. It
is shown numerically that, under the diffuse-reflection boundary condition, the difference
between the plate velocity and its final value decreases in proportion to an inverse power of
time. This agrees with the previous theoretical result obtained under the assumption that
the initial plate velocity is sufficiently close to the final one.
I.
はじめに
一様で静止した気体中を，物体が一定外力によって運動する問題を考えます．物体は外力に平行
な方向に初速度 vw0 (> 0) で運動を始めますが，拘束または対称性により，物体は回転することな
く外力と平行な方向に運動します．気体による抵抗が物体には働き，物体の速度は外力と抵抗が釣
り合うような一定値（最終速度）に近づいていきます．物体に働く抵抗が物体の速さに比例すると
仮定すると，その速度は最終速度に指数関数的に収束します．つまり，τ を時間変数とし，vw (τ )
を外力が働いている方向の物体の速度，vw∞ (> 0) を最終速度とすると，十分大きな τ に対して，
|vw∞ − vw (τ )| ≈ C1 e−C2 τ
(1)
が成り立ちます．ここで C1 , C2 は正の定数です．
ところが，気体が分子間の衝突が無視できる程度に高度に希薄な理想気体（自由分子気体）の
場合は，物体の速度の最終速度への近づき方は (1) とは大きく異なることが既往研究で数学的解析
により示されています [1-4]．既往研究では，対象は空間的 d 次元の円盤で (d = 1, 2, 3 )，円盤はそ
の進行方向に対して垂直です（Fig. 1）．この場合，円盤の最終速度への近づき方は遅く，時間 τ の
逆べきになります．つまり，十分大きな τ に対して，
|vw∞ − vw (τ )| ≈ C1′ τ −n
(2)
が成り立ちます．ここで C1′ は正の定数，n は d と気体分子と円盤表面の相互作用によって決定さ
れる自然数です．気体分子が円盤表面で鏡面反射（完全弾性衝突）するとき n = d + 2 で [1,2]，こ
Figure 1: A circular disk in a free-molecular gas. (a) Circular disk (d = 3), (b) infinite plate with
a rectangular cross section (d = 2), (c) infinite plate with a finite thickness (d = 1).
れは物体が任意の凸型物体であるときに対しても成り立ちます [3]．気体分子が円盤表面で拡散反
射するときは n = d + 1 です [4]．拡散反射とは，円盤表面に衝突した気体分子の衝突後の速度が，
円盤の温度や速度によって特徴づけられる等方的なガウス分布に従うというもので，分子気体力学
の分野ではしばしば用いられる代表的な気体分子と物体表面間の相互作用のモデルです．(2) のよ
うな最終速度への遅いアプローチは，円盤表面を過去に反射した気体分子が再び円盤に衝突するこ
とによります（再衝突）．さらに鏡面反射条件下では，初速度と最終速度が 0 ≤ vw∞ < vw0 の関
係にあるときには円盤の速度 vw (τ ) はその最終速度 vw∞ を一度下回り，その後に最終速度に近づ
いていくことが示されています [2]．
これらの結果は，厳密な数学的解析によるものですが，初速度 vw0 と最終速度 vw∞ との差が十
分に小さい，という条件下でのみ成り立つものです．一方，(2) の遅いアプローチは前述の再衝突
によるものなので，同じ結果が任意の初速度 vw0 や最終速度 vw∞ で成り立つかは明らかではあり
ません．本研究では，気体分子と物体表面の相互作用が拡散反射である問題 (n = d + 1) に対して，
数値計算によって (2) が成り立つかを調べました．
II.
A.
問題の定式化
問題設定
数値計算を簡単にするために，本研究では円盤の代わりに厚さの無視できる平板を対象の物体と
します（Fig. 2）．数値計算は d = 1, 2, 3 に対して行いましたが，ここでは d = 2 の二次元問題に
対して定式化を行います [Fig. 2(b)]．
Figure 2: A rectangular plate without thickness in a free-molecular gas. (a) Rectangular plate
(d = 3), (b) infinite plate with a finite width (d = 2), (c) infinite plate (d = 1).
温度 T0 ，密度 ρ0 で一様な，静止平衡状態にある理想気体を考え，その中に幅 L の無限に長い
厚さの無視できる平板が温度 T0 で気体中に固定されているとします．直交座標系を Fig.2(b) のよ
うにとりますが，このとき平板は X1 = 0，−L/2 ≤ X2 ≤ L/2，−∞ < X3 < ∞ に位置していて，
平板には単位質量あたりの一定外力 F(≥ 0) が X1 方向に働いています．時刻 τ = 0 で平板はその
拘束から解放され，初速度 vw0 (≥ 0) で X1 方向に運動を始めます．平板は X1 軸に沿って運動しま
すが，次第に一定速度の定常状態（このときの速度を最終速度と呼ぶことにします）に近づいてい
きます．本研究では，下に示す仮定のもとで，物体の運動を，特にその最終速度へのアプローチの
仕方を調べます．
(i) 気体の振舞はボルツマン方程式 [5] で記述される.
(ii) 気体は高度に希薄なため，気体分子同士の衝突は無視でき (自由分子気体 [5]), 気体分子には
外力は働かない.
(iii) 境界（平板表面）に入射した気体分子は境界と十分に相互作用し，境界面の条件で規定され
る平衡分布で反射され，また境界面に運ばれる正味の質量は零である．
（拡散反射 [5]）
基礎方程式を示す前に，用いる記号をまとめます．有次元の物理量： Xi は空間における直交座
標系，τ は時間変数，ξi は分子速度，Xw (τ ) は時刻 τ における平板の位置 (X1 座標)，vw (τ ) は対
応する X1 方向の平板の速度，vw0 は平板の初速度，vw∞ は平板の最終速度，f˜(Xi , ξi , τ ) は気体分
子の速度分布関数，F は平板に働く単位質量あたりの外力（X1 方向），G(τ ) は平板に働く単位質
量あたりの抵抗（X1 方向），M は平板の単位面積あたりの質量です．それらに対応する無次元の
物理量 xi , t, ζi , xw , uw , uw0 , uw∞ , f, F, G, M は次のように定義されます．
xi = Xi /L,
ζi = ξi /(2kT0 /m∗ )1/2 ,
uw (t) = vw (τ )/(2kT0 /m∗ )1/2 ,
uw∞ = vw∞ /(2kT0 /m∗ )1/2 ,
F = F/Lτ0−2 ,
M = M/ρ0 L.
t = τ /τ0 ,
xw (t) = Xw (τ )/L,
uw0 = vw0 /(2kT0 /m∗ )1/2 ,
f (xi , ζi , t) = f˜(Xi , ξi , τ )/ρ0 (2kT0 /m∗ )−3/2 ,
G(t) = G(τ )/Lτ0−2 ,
(3)
ここで，基準時間は τ0 = L(2kT0 /m∗ )−1/2 で定義され，k はボルツマン定数，m∗ は気体分子の質
量です．
B.
基礎方程式
二次元問題に対する定式化では物理量は x3 に依存しないため，方程式から分子速度の ζ3 を次の
ような周辺分布関数を定義することで消去することができます．
∫ ∞
g(x1 , x2 , ζ1 , ζ2 , t) =
f (x1 , x2 , ζi , t) dζ3 .
(4)
−∞
このとき，自由分子気体に対するボルツマン方程式は
∂g
∂g
∂g
+ ζ1
+ ζ2
=0
∂t
∂x1
∂x2
(5)
と表されます．対応する初期条件は
g0 = π −1 exp(−ζ12 − ζ22 ),
g = g0 ,
(t = 0).
(6)
平板上での境界条件（拡散反射条件）は
g(x1 , x2 , ζ1 , ζ2 , t) = gw± (x2 , ζ1 , ζ2 , t),
]
[
1
1
x1 = xw± (t), − ≤ x2 ≤ , ζ1 − uw (t) ≷ 0 ,
2
2
(
)
gw± (x2 , ζ1 , ζ2 , t) = π −1 ρw± (x2 , t) exp − [ζ1 − uw (t)]2 − ζ22 ,
∫ ∞∫
√
ρw± (x2 , t) = ∓2 π
[ζ1 − uw (t)]
−∞
(7a)
(7b)
ζ1 −uw (t)≶0
× g(xw± (t), x2 , ζ1 , ζ2 , t) dζ1 dζ2 .
(7c)
ここで x1 = xw± は x1 = xw ± 0 を意味していて，{ (x1 , x2 )|x1 = xw± , −1/2 ≤ x2 ≤ 1/2 } は平板
の表面の x1 が正の側の面（＋の記号）と負の側の面（−の記号）を表しています．
平板の運動方程式は
duw (t)/dt = F − G(t).
dxw (t)/dt = uw (t),
(8)
平板に働く抵抗 G(t) は平板上の速度分布関数を用いて
G± (t) = ±
1
M
∫
1/2
−1/2
G(t) = G+ (t) + G− (t),
{∫ ∞ ∫
[ζ1 − uw (t)]2
−∞
∫
∞
ζ1 −uw (t)≶0
∫
× g(xw± , x2 , ζ1 , ζ2 , t) dζ1 dζ2
+
−∞
(9a)
ζ1 −uw (t)≷0
[ζ1 − uw (t)]2
}
× gw± (x2 , ζ1 , ζ2 , t) dζ1 dζ2 dx2
(9b)
と表せます．G± (t) はそれぞれ平板の表面 xw± (t) に働く抵抗です．運動方程式に対応する初期条
件は
xw (0) = 0,
uw (0) = uw0 .
(5)∼(7) と (8)∼ (10) の連立方程式に対して数値計算を行います．
(10)
III.
数値計算の結果
時間が十分経過したとき (t → ∞)，平板の速度 uw (t) はその最終速度 uw∞ ，つまり平板に働く
抵抗 G と外力 F が釣り合うような速度に近づいていきます．ある M に対して，最終速度 uw∞ は
一意に外力 F によって決まるものなので，実際の計算では F の代わりに uw∞ をパラメータとし
ます．よって，問題を特徴付けるパラメータは uw0 ，uw∞ ，M ですが，ここでは既往研究（文献
[1-4]）に合わせるため M = 1 とします．
平板の速度 uw (t) の最終速度 uw∞ への近づき方，言い換えれば |uw∞ − uw (t)| の減衰の仕方に
興味があるので [(2) 式参照]，|uw∞ − uw (t)| の両対数スケール（log は常用対数），および以下に
定義する両対数スケールにおける傾き αu を中心に，一次元問題と三次元問題結果の説明をします．
αu =
A.
d log |uw∞ − uw (t)|
.
d log t
(11)
一次元問題 (d = 1)
一次元問題では平板は x2 x3 平面に無限に広い [Fig.2(c) 参照] ので，物理量は x1 にしか依存しま
せん．このとき (4) の代わりに次の周辺分布関数 ḡ(x1 , ζ1 , t) を用いることで分子速度 ζ2 も消去す
ることができます．
∫ ∞∫ ∞
f (x1 , ζi , t)dζ2 dζ3 .
(12)
g(x1 , ζ1 , t) =
−∞
−∞
計算に用いたパラメータを示します．Case 1∼Case 6 は初速度 uw0 が最終速度 uw∞ より小さ
い場合です (0 ≤ uw0 < uw∞ )．Case 7 は Case 1 と同じですが，再衝突 (I 節参照) の効果を無視し
ています．一方，Case 8∼Case 10 は 0 ≤ uw∞ < uw0 の場合です．Case 10 は Case 8 で再衝突の
効果を無視したものです．
Case
Case
Case
Case
Case
Case
1:
3:
5:
7:
8:
10:
(uw∞ ,
(uw∞ ,
(uw∞ ,
(uw∞ ,
(uw∞ ,
(uw∞ ,
uw0 ) = (1.5, 0),
uw0 ) = (3.55659, 0),
uw0 ) = (1.5, 1.125),
uw0 ) = (1.5, 0) [no recollision],
uw0 ) = (0, 1),
uw0 ) = (0, 1) [no recollision].
Case 2:
Case 4:
Case 6:
(uw∞ , uw0 ) = (2.35815, 0),
(uw∞ , uw0 ) = (1.5, 0.75),
(uw∞ , uw0 ) = (1.5, 1.3125),
Case 9:
(uw∞ , uw0 ) = (1.5, 6),
Figure 3(a) は Case 1∼Case 7 に対する log |uw∞ − uw (t)| vs log t，Fig. 3(b) は対応する αu vs
log t を示したものです．Figure 3(a) より，再衝突を無視した Case 7 を除いて，およそ t > 10 で
log |uw∞ − uw (t)| は log t に比例して減衰しています．Figure 3(b) はその傾きが −2 に近づいてい
ることを示していますが，これは最終速度と初速度が十分に小さいという条件の下で数学的に導か
Figure 3: Long-time behavior for 0 ≤ uw0 < uw∞ (one-dimensional problem). (a) log |uw∞ − uw |
vs log t, (c) αu vs log t.
れた結果，つまり (2) 式 (n = d + 1) と一致しています．一方，再衝突を無視した Case 7 では最終
速度 uw∞ への近づき方はより速いことが分かります．Case 7 を片対数でプロットすれば，このと
き uw∞ − uw (t) は (1) 式で表わされるような指数関数的な収束の仕方になります．また，これら計
算では uw∞ − uw (t) は常に正値で，単調減少です．
続いて Case 8∼Case 10 の計算結果を Fig. 4 に示します．Figure 4(a) は x1 = xw (t) の軌跡を x1 t
平面に Case 8(実線) と Case 10(点線) に対して示したもので，これらの Case では外力が働いていな
いため平板は時間が十分経過した極限 (t → ∞) で静止します．この図から分かるように，Case 8 で
は xw (t) は時間に対して単調増加ではありません．つまり，平板の位置はその最終値 limt→∞ xw (t)
を少し越えて，その後に最終値に戻ってきます．対照的に，Case 10 ではこのようなオーバーシュー
トは起こっていません．Figure 4(b) は Fig. 4(a) の拡大図で，Fig. 4(c) は uw∞ − uw (t) を Case
8∼Case 10 に対して示したものです．平板の速度 uw (t) は最初 uw∞ より大きいですが，一旦 uw∞
を下回って，その後に下から uw∞ に近づいていきます．Case 8 と Case 9 では uw∞ − uw (t) の極
値を与える時刻がそれぞれ t = t∗ = 3.4775 と 2.6125 と得られていて，それに対応して x1 = xw (t)
は t = t∗ で変曲点をむかえます [Fig. 4(b)]. 再衝突の効果を無視した Case 10 ではこのようなオー
バーシュートは起こりません．Figure 4(d) では αu が −2 に近づいていることが分かります．この
平板の位置 xw (t) と平板の速度差 uw∞ − uw (t) のオーバーシュートは，境界条件 (7) が鏡面反射条
件の場合に対して，最終速度と初速度の差が十分小さいという条件の下で示されています [2]．
Figure 4: Time evolution for 0 ≤ uw0 < uw∞ (one-dimensional problem). (a) trajectory x1 =
xw (t), (b) magnified figure of (a), (c) uw∞ − uw (t) vs t, (d) αu vs log t.
B.
三次元問題 (d = 3)
平板はその辺の長さがそれぞれ L, H であるとし，初期に X1 = 0, −L/2 ≤ X2 ≤ L/2, −H/2 ≤
X3 ≤ H/2 に位置しています．また (4) のような周辺分布関数は用いずに，もとの速度分布関数
f (xi , ζi , t) を扱います．
Figure 5 は Case 1 に対して，アスペクト比 H/L が H/L = 1, 2, 4, 8, 16 としたときの計算結果
で，Fig. 5 では二次元問題の結果 H/L → ∞ (点線) も合わせて示してあります．Figure 5(a)，(b)
から αu は −4 に近づいていく傾向があり，これは (2) 式 (n = d + 1) と一致します．また，Fig. 5(a)
から，アスペクト比 H/L を大きくすればするほど二次元問題の結果に近づいていくことが分かり
ますが，最終的に得られる傾き αu は二次元問題と三次元問題では明らかに違います．二次元問題
の αu は示していませんが，この場合は −3 に近づいていて，やはり (2) 式 (n = d + 1) と一致し
ます．
Figure 5: Long-time behavior for Case 1 with various aspect ratios (three-dimensional problem).
(a) log |uw∞ − uw | vs log t, (b) αu vs log t.
しかしながら，一次元問題の結果 Fig. 3(b) と Fig. 5(b) では示している αu と log t の範囲がか
なり違うことに注意すれば，明確な結果を得るためにはより大きい時間 t までの計算が必要である
ことが分かります．三次元問題の計算は時間を要するので，大きい t までの計算は H/L = 1 の場
合に関してのみ行いました．より詳細な結果から (2) 式 (n = d + 1) を確かめましたが，ここでは
結果は省略します．
IV.
おわりに
本研究では，自由分子気体中を一定外力のもとで運動する物体（平板）の最終的な定常状態への
近づき方に特に注目し，その振舞を数値計算を行い調べました．平板ではなく円盤を用いた既往研
究 [1-4] の数学的結果 (2) は平板の初速度と最終速度が十分に小さいという条件に基づいていまし
たが，本研究では任意の初速度，最終速度に対して同じ結果が成り立つことを数値的に示しました
（ただし境界条件は拡散反射条件）．
自由分子気体中の凸型物体の運動に関してはいくつか研究報告がありますが，それらの研究で
は再衝突の効果は考慮されていません．物体の運動が加速，減速，回転を伴う非定常運動の場合は，
一般に再衝突の効果を無視することはできませんので，物体の振舞は再衝突を無視した場合とは異
なることが予想されます．これらの変化を調べることも今後の興味深い課題の一つとして挙げられ
ます．
参考文献
[1] S. Caprino, C. Marchioro, and M. Pulvirenti, “Approach to equilibrium in a microscopic
model of friction,” Commun. Math. Phys. 264, 167–189 (2006).
[2] S. Caprino, G. Cavallaro, and C. Marchioro, “On a microscopic model of viscous friction,”
Math. Models Methods Appl. Sci. 17, 1369–1403 (2007).
[3] G. Cavallaro, “On the motion of a convex body interacting with a perfect gas in the meanfield approximation,” Rendiconti di Matematica. Ser. VII 27, 123–145 (2007).
[4] K. Aoki, G. Cavallaro, C. Marchioro, and M. Pulvirenti, “On the motion of a body in
thermal equilibrium immersed in a perfect gas,” ESAIM: Math. Model. Num. Anal. 42,
263–275 (2008).
[5] Y. Sone, Molecular Gas Dynamics: Theory, Techniques, and Applications (Birkäuser,
Boston, 2007).