近似の数理 第8章 §8.6 ★★★ 粘性係数と熱伝導率が分子の衝突から計算できること、BGK 模型で計算を最後まで実行する
Hilbert 展開との違い
Chapman（1916–17）と Enskog（1917）は、次の二点で展開を組み替えた。
場 n , u , T n,\bm{u},T n , u , T 真の解 f f f で定義し、展開しない。補正 f 1 , f 2 , … f_1,f_2,\dots f 1 , f 2 , … 1 , v , ∣ v ∣ 2 1,\bm{v},\lvert \bm{v}\rvert ^2 1 , v , ∣ v ∣ 2
f f f ε \varepsilon ε ∂ t = ∂ t ( 0 ) + ε ∂ t ( 1 ) + ⋯ \partial_t=\partial_t^{(0)}+\varepsilon \partial_t^{(1)}+\cdots ∂ t = ∂ t ( 0 ) + ε ∂ t ( 1 ) + ⋯ ∂ t ( 0 ) \partial_t^{(0)} ∂ t ( 0 ) ∂ t ( 1 ) \partial_t^{(1)} ∂ t ( 1 )
この組み替えによって、O ( ε ) \mathcal{O}(\varepsilon ) O ( ε ) ε \varepsilon ε
BGK 模型
剛体球の L − 1 \mathcal{L}^{-1} L − 1
∂ t f + v ⋅ ∇ x f = M f − f τ (8.38) \partial_tf+\bm{v}\cdot\nabla_{\bm{x}}f=\frac{M_f-f}{\tau}\tag{8.38} ∂ t f + v ⋅ ∇ x f = τ M f − f ( 8.38 )
を使う。M f M_f M f f f f n , u , T n,\bm{u},T n , u , T τ \tau τ Q Q Q
M f M_f M f f f f 1 , v , ∣ v ∣ 2 1,\bm{v},\lvert \bm{v}\rvert ^2 1 , v , ∣ v ∣ 2 ∫ ( M f − f ) ϕ d v = 0 \int(M_f-f)\phi\,\mathrm{d}\bm{v}=0 ∫ ( M f − f ) ϕ d v = 0 右辺が 0 0 0 f = M f f=M_f f = M f
H 定理が成り立つ：∫ ( M f − f ) ln f d v = ∫ ( M f − f ) ( ln f − ln M f ) d v ≤ 0 \int(M_f-f)\ln f\,\mathrm{d}\bm{v}=\int(M_f-f)(\ln f-\ln M_f)\,\mathrm{d}\bm{v}\le0 ∫ ( M f − f ) ln f d v = ∫ ( M f − f ) ( ln f − ln M f ) d v ≤ 0 ln M f \ln M_f ln M f ( a − b ) ( ln b − ln a ) ≤ 0 (a-b)(\ln b-\ln a)\le0 ( a − b ) ( ln b − ln a ) ≤ 0
BGK は C 型の近似である
BGK 模型は、Boltzmann–Grad 極限のような小さいパラメータの極限から出てくるものではない。衝突演算子という未知の複雑な項を、保存則と H 定理を満たす単純な形で置き換えたモデルであり、第1章の分類では C 型（閉じる）に属する。そのため保存則・平衡・H 定理のような構造的な性質は正しく再現するが、定量的な輸送係数の比（次に見る Prandtl 数）は外れる。
1 次の補正
以下では次元を戻し、τ \tau τ τ → ε τ \tau\to\varepsilon \tau τ → ε τ f = M + f ( 1 ) + ⋯ f=M+f^{(1)}+\cdots f = M + f ( 1 ) + ⋯ O ( τ 0 ) \mathcal{O}(\tau^0) O ( τ 0 )
f ( 1 ) = − τ ( ∂ t ( 0 ) + v ⋅ ∇ x ) M = − τ M D ln M , D : = ∂ t ( 0 ) + v ⋅ ∇ x f^{(1)}=-\tau(\partial_t^{(0)}+\bm{v}\cdot\nabla_{\bm{x}})M=-\tau M\,D\ln M,\qquad D:=\partial_t^{(0)}+\bm{v}\cdot\nabla_{\bm{x}} f ( 1 ) = − τ ( ∂ t ( 0 ) + v ⋅ ∇ x ) M = − τ M D ln M , D := ∂ t ( 0 ) + v ⋅ ∇ x
である（BGK の線型化演算子は ( h − Π h ) / τ (h-\Pi h)/\tau ( h − Π h ) / τ τ \tau τ D ln M D\ln M D ln M ln M = ln n − 3 2 ln T − m ∣ c ∣ 2 2 k B T + 定数 \ln M=\ln n-\frac32\ln T-\frac{m\lvert \bm{c}\rvert ^2}{2 k_{\mathrm B} T}+\text{定数} ln M = ln n − 2 3 ln T − 2 k B T m ∣ c ∣ 2 + 定数 n , u , T n,\bm{u},T n , u , T
( ∂ t + u ⋅ ∇ ) n = − n ∇ ⋅ u , ( ∂ t + u ⋅ ∇ ) u = − ∇ p m n , ( ∂ t + u ⋅ ∇ ) T = − 2 3 T ∇ ⋅ u (\partial_t+\bm{u}\cdot\nabla)n=-n\nabla\cdot\bm{u},\qquad
(\partial_t+\bm{u}\cdot\nabla)\bm{u}=-\frac{\nabla p}{mn},\qquad
(\partial_t+\bm{u}\cdot\nabla)T=-\frac23T\nabla\cdot\bm{u} ( ∂ t + u ⋅ ∇ ) n = − n ∇ ⋅ u , ( ∂ t + u ⋅ ∇ ) u = − mn ∇ p , ( ∂ t + u ⋅ ∇ ) T = − 3 2 T ∇ ⋅ u
で置き換え、v ⋅ ∇ = u ⋅ ∇ + c ⋅ ∇ \bm{v}\cdot\nabla=\bm{u}\cdot\nabla+\bm{c}\cdot\nabla v ⋅ ∇ = u ⋅ ∇ + c ⋅ ∇ ∇ p / ( n k B T ) = ∇ n / n + ∇ T / T \nabla p/(n k_{\mathrm B} T)=\nabla n/n+\nabla T/T ∇ p / ( n k B T ) = ∇ n / n + ∇ T / T ∇ n \nabla n ∇ n ∇ ⋅ u \nabla\cdot\bm{u} ∇ ⋅ u − 1 + 1 − m ∣ c ∣ 2 3 k B T -1+1-\frac{m\lvert \bm{c}\rvert ^2}{3 k_{\mathrm B} T} − 1 + 1 − 3 k B T m ∣ c ∣ 2
D ln M = m k B T ( c i c j − 1 3 ∣ c ∣ 2 δ i j ) ∂ j u i + ( m ∣ c ∣ 2 2 k B T − 5 2 ) c ⋅ ∇ T T (8.41) D\ln M=\frac{m}{ k_{\mathrm B} T}\Bigl(c_ic_j-\frac13\lvert \bm{c}\rvert ^2\delta_{ij}\Bigr)\partial_ju_i+\Bigl(\frac{m\lvert \bm{c}\rvert ^2}{2 k_{\mathrm B} T}-\frac52\Bigr)\frac{\bm{c}\cdot\nabla T}{T}\tag{8.41} D ln M = k B T m ( c i c j − 3 1 ∣ c ∣ 2 δ ij ) ∂ j u i + ( 2 k B T m ∣ c ∣ 2 − 2 5 ) T c ⋅ ∇ T ( 8.41 )
である（検算スクリプト。Euler 方程式で時間微分を消去した式を記号計算で照合した）。第1項は速度勾配のうち対称でトレースのない部分（ずれ変形）、第2項は温度勾配にだけ比例し、密度勾配は残らない。また式 (8.41) は 1 , c , ∣ c ∣ 2 1,\bm{c},\lvert \bm{c}\rvert ^2 1 , c , ∣ c ∣ 2 M M M c \bm{c} c 1 1 1 ∣ c ∣ 2 \lvert \bm{c}\rvert ^2 ∣ c ∣ 2 0 0 0 ∫ ∣ c ∣ 2 c i c j M = 5 n s 2 δ i j \int\lvert \bm{c}\rvert ^2c_ic_jM=5ns^2\delta_{ij} ∫ ∣ c ∣ 2 c i c j M = 5 n s 2 δ ij ∫ ∣ c ∣ 4 M = 15 n s 2 \int\lvert \bm{c}\rvert ^4M=15ns^2 ∫ ∣ c ∣ 4 M = 15 n s 2 1 / 3 1/3 1/3 c \bm{c} c 0 0 0 c \bm{c} c 1 , ∣ c ∣ 2 1,\lvert \bm{c}\rvert ^2 1 , ∣ c ∣ 2 0 0 0 c k c_k c k ( 5 n s 2 2 s − 5 2 n s ) δ k j = 0 \bigl(\frac{5ns^2}{2s}-\frac52ns\bigr)\delta_{kj}=0 ( 2 s 5 n s 2 − 2 5 n s ) δ k j = 0
粘性率・熱伝導率・Prandtl 数
Maxwell 分布についての Gauss 積分
∫ c i c j M d v = n s δ i j , ∫ c i c j c k c l M d v = n s 2 ( δ i j δ k l + δ i k δ j l + δ i l δ j k ) , ∫ ∣ c ∣ 2 c i c j M d v = 5 n s 2 δ i j , ∫ ∣ c ∣ 4 c i c j M d v = 35 n s 3 δ i j ( s = k B T / m ) \begin{gather*}\int c_ic_jM\,\mathrm{d}\bm{v}=ns\,\delta_{ij},\qquad \int c_ic_jc_kc_lM\,\mathrm{d}\bm{v}=ns^2(\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}),\\\int\lvert \bm{c}\rvert ^2c_ic_jM\,\mathrm{d}\bm{v}=5ns^2\delta_{ij},\qquad \int\lvert \bm{c}\rvert ^4c_ic_jM\,\mathrm{d}\bm{v}=35ns^3\delta_{ij}\qquad(s= k_{\mathrm B} T/m)\end{gather*} ∫ c i c j M d v = n s δ ij , ∫ c i c j c k c l M d v = n s 2 ( δ ij δ k l + δ ik δ j l + δ i l δ j k ) , ∫ ∣ c ∣ 2 c i c j M d v = 5 n s 2 δ ij , ∫ ∣ c ∣ 4 c i c j M d v = 35 n s 3 δ ij ( s = k B T / m )
を使うと、f ( 1 ) f^{(1)} f ( 1 ) 8.7 ）。結果は
P i j = p δ i j − μ ( ∂ i u j + ∂ j u i − 2 3 δ i j ∇ ⋅ u ) , q = − κ ∇ T , (8.44) P_{ij}=p\,\delta_{ij}-\mu\Bigl(\partial_iu_j+\partial_ju_i-\frac23\delta_{ij}\nabla\cdot\bm{u}\Bigr),\qquad \bm{q}=-\kappa\nabla T,\tag{8.44} P ij = p δ ij − μ ( ∂ i u j + ∂ j u i − 3 2 δ ij ∇ ⋅ u ) , q = − κ ∇ T , ( 8.44 )
μ = p τ = n k B T τ , κ = 5 2 k B m p τ (8.45) \mu=p\tau=n k_{\mathrm B} T\tau,\qquad \kappa=\frac52\,\frac{ k_{\mathrm B}}{m}\,p\tau\tag{8.45} μ = p τ = n k B T τ , κ = 2 5 m k B p τ ( 8.45 )
である（検算スクリプト）。体積粘性（∇ ⋅ u \nabla\cdot\bm{u} ∇ ⋅ u 0 0 0
定理 8.10（BGK 模型の Chapman–Enskog 1 次 ★★★）
BGK 模型 (8.38) の Chapman–Enskog 展開を O ( τ ) \mathcal{O}(\tau) O ( τ ) (8.22) –(8.24) は構成則 (8.44) 、輸送係数 (8.45) の圧縮性 Navier–Stokes–Fourier 方程式になる。Prandtl 数（運動量の拡散と熱の拡散の比）は
P r = c p μ κ = 5 2 k B m ⋅ p τ 5 2 k B m p τ = 1 \mathrm{Pr}=\frac{c_p\,\mu}{\kappa}=\frac{\frac52\frac{ k_{\mathrm B}}{m}\cdot p\tau}{\frac52\frac{ k_{\mathrm B}}{m}\,p\tau}=1 Pr = κ c p μ = 2 5 m k B p τ 2 5 m k B ⋅ p τ = 1
である（c p = 5 2 k B / m c_p=\frac52 k_{\mathrm B}/m c p = 2 5 k B / m
実際の単原子気体（アルゴンなど）の Prandtl 数は約 2 / 3 2/3 2/3 κ = 15 4 k B m μ \kappa=\frac{15}{4}\frac{ k_{\mathrm B}}{m}\mu κ = 4 15 m k B μ P r = 5 / 2 15 / 4 = 2 3 \mathrm{Pr}=\frac{5/2}{15/4}=\frac23 Pr = 15/4 5/2 = 3 2 τ \tau τ
剛体球の粘性率
剛体球の Chapman–Enskog 計算の結果は
μ H S = 1.016 … × 5 16 d 2 m k B T π \mu_{\mathrm{HS}}=1.016\ldots\times\frac{5}{16\,d^2}\sqrt{\frac{m k_{\mathrm B} T}{\pi}} μ HS = 1.016 … × 16 d 2 5 π m k B T
である。5 16 d 2 m k B T / π \frac{5}{16d^2}\sqrt{m k_{\mathrm B} T/\pi} 16 d 2 5 m k B T / π 1.016 … 1.016\ldots 1.016 … 密度によらない 。圧力を下げて分子が少なくなっても粘性は変わらない（運ぶ分子が減る分だけ平均自由行程が伸びて、遠くまで運動量を運ぶ）。Maxwell がこの性質を予言し、実験で確かめたことは分子運動論の初期の成功として知られる。BGK でも τ ∝ 1 / n \tau\propto1/n τ ∝ 1/ n μ = p τ \mu=p\tau μ = p τ n n n 28.97 28.97 28.97 T = 293 K T=293\,\mathrm{K} T = 293 K d ≈ 0.37 n m d\approx0.37\,\mathrm{nm} d ≈ 0.37 nm μ H S ≈ 1.83 × 10 − 5 P a s \mu_{\mathrm{HS}}\approx1.83\times10^{-5}\,\mathrm{Pa\,s} μ HS ≈ 1.83 × 1 0 − 5 Pa s 1.81 × 10 − 5 P a s 1.81\times10^{-5}\,\mathrm{Pa\,s} 1.81 × 1 0 − 5 Pa s d ≈ 0.37 n m d\approx0.37\,\mathrm{nm} d ≈ 0.37 nm
例題 8.7 BGK 模型から粘性率と熱伝導率を求める ★★★
オリジナル [標準的な題材を自作]
BGK 模型の Chapman–Enskog 1 次の補正 f ( 1 ) = − τ M D ln M f^{(1)}=-\tau M\,D\ln M f ( 1 ) = − τ M D ln M D ln M D\ln M D ln M (8.41) ）について、次を求めよ。s = k B T / m s= k_{\mathrm B} T/m s = k B T / m
温度一様な単純ずれ流れ u = ( γ y , 0 , 0 ) \bm{u}=(\gamma y,0,0) u = ( γ y , 0 , 0 ) γ \gamma γ P x y ( 1 ) = m ∫ c x c y f ( 1 ) d v P^{(1)}_{xy}=m\int c_xc_yf^{(1)}\mathrm{d}\bm{v} P x y ( 1 ) = m ∫ c x c y f ( 1 ) d v μ = p τ \mu=p\tau μ = p τ
静止した気体（u = 0 \bm{u}=\bm{0} u = 0 ∇ T = ( T x , 0 , 0 ) \nabla T=(T_x,0,0) ∇ T = ( T x , 0 , 0 ) q x ( 1 ) = m 2 ∫ ∣ c ∣ 2 c x f ( 1 ) d v q^{(1)}_x=\frac{m}{2}\int\lvert \bm{c}\rvert ^2c_xf^{(1)}\mathrm{d}\bm{v} q x ( 1 ) = 2 m ∫ ∣ c ∣ 2 c x f ( 1 ) d v κ \kappa κ
Prandtl 数 P r = c p μ / κ \mathrm{Pr}=c_p\mu/\kappa Pr = c p μ / κ
解答 (1) ∂ j u i \partial_ju_i ∂ j u i 0 0 0 ∂ y u x = γ \partial_yu_x=\gamma ∂ y u x = γ ∇ ⋅ u = 0 \nabla\cdot\bm{u}=0 ∇ ⋅ u = 0 (8.41) は D ln M = m k B T c x c y γ = γ s c x c y D\ln M=\frac{m}{ k_{\mathrm B} T}c_xc_y\gamma=\frac{\gamma}{s}c_xc_y D ln M = k B T m c x c y γ = s γ c x c y
P x y ( 1 ) = − τ m γ s ∫ c x 2 c y 2 M d v = − τ m γ s n s 2 = − τ m n s γ = − n k B T τ γ . P^{(1)}_{xy}=-\tau m\frac{\gamma}{s}\int c_x^2c_y^2M\,\mathrm{d}\bm{v}=-\tau m\frac{\gamma}{s}\,ns^2=-\tau\,mns\,\gamma=-n k_{\mathrm B} T\tau\,\gamma. P x y ( 1 ) = − τ m s γ ∫ c x 2 c y 2 M d v = − τ m s γ n s 2 = − τ mn s γ = − n k B T τ γ .
ここで ∫ c x 2 c y 2 M = n s 2 ( δ x x δ y y + 0 + 0 ) = n s 2 \int c_x^2c_y^2M=ns^2(\delta_{xx}\delta_{yy}+0+0)=ns^2 ∫ c x 2 c y 2 M = n s 2 ( δ xx δ y y + 0 + 0 ) = n s 2 P x y = − μ ∂ y u x = − μ γ P_{xy}=-\mu\,\partial_yu_x=-\mu\gamma P x y = − μ ∂ y u x = − μ γ μ = n k B T τ = p τ \mu=n k_{\mathrm B} T\tau=p\tau μ = n k B T τ = p τ
(2) u = 0 \bm{u}=\bm{0} u = 0 D ln M = ( ∣ c ∣ 2 2 s − 5 2 ) c x T x T D\ln M=\bigl(\frac{\lvert \bm{c}\rvert ^2}{2s}-\frac52\bigr)\frac{c_xT_x}{T} D ln M = ( 2 s ∣ c ∣ 2 − 2 5 ) T c x T x
q x ( 1 ) = − τ m 2 T x T ∫ ∣ c ∣ 2 c x 2 ( ∣ c ∣ 2 2 s − 5 2 ) M d v = − τ m 2 T x T ( 35 n s 3 2 s − 5 2 ⋅ 5 n s 2 ) = − τ m 2 T x T ⋅ 5 n s 2 = − 5 2 n k B 2 T τ m T x . \begin{aligned}q^{(1)}_x&=-\tau\frac{m}{2}\frac{T_x}{T}\int\lvert \bm{c}\rvert ^2c_x^2\Bigl(\frac{\lvert \bm{c}\rvert ^2}{2s}-\frac52\Bigr)M\,\mathrm{d}\bm{v}
=-\tau\frac{m}{2}\frac{T_x}{T}\Bigl(\frac{35ns^3}{2s}-\frac52\cdot5ns^2\Bigr)\\
&=-\tau\frac{m}{2}\frac{T_x}{T}\cdot5ns^2=-\frac52\,\frac{n k_{\mathrm B}^2T\tau}{m}\,T_x.\end{aligned} q x ( 1 ) = − τ 2 m T T x ∫ ∣ c ∣ 2 c x 2 ( 2 s ∣ c ∣ 2 − 2 5 ) M d v = − τ 2 m T T x ( 2 s 35 n s 3 − 2 5 ⋅ 5 n s 2 ) = − τ 2 m T T x ⋅ 5 n s 2 = − 2 5 m n k B 2 T τ T x .
Fourier の法則 q x = − κ T x q_x=-\kappa T_x q x = − κ T x κ = 5 2 k B m n k B T τ = 5 2 k B m p τ \kappa=\frac52\frac{ k_{\mathrm B}}{m}\,n k_{\mathrm B} T\tau=\frac52\frac{ k_{\mathrm B}}{m}\,p\tau κ = 2 5 m k B n k B T τ = 2 5 m k B p τ
(3) P r = ( 5 k B / 2 m ) p τ ( 5 k B / 2 m ) p τ = 1 \mathrm{Pr}=\frac{(5 k_{\mathrm B}/2m)\,p\tau}{(5 k_{\mathrm B}/2m)\,p\tau}=1 Pr = ( 5 k B /2 m ) p τ ( 5 k B /2 m ) p τ = 1 μ / ρ \mu/\rho μ / ρ κ / ( ρ c p ) \kappa/(\rho c_p) κ / ( ρ c p )
次の次数と展開の限界
Chapman–Enskog 展開を O ( ε 2 ) \mathcal{O}(\varepsilon ^2) O ( ε 2 ) K n \mathrm{Kn} Kn K n \mathrm{Kn} Kn
Navier–Stokes は極限方程式ではない
圧縮性 Navier–Stokes 方程式の粘性項は O ( K n ) \mathcal{O}(\mathrm{Kn}) O ( Kn ) K n → 0 \mathrm{Kn}\to0 Kn → 0
AI（Claude）が書いた教材を、木下が勉強しながら公開しています。式と数値は検算していますが、誤りが残っているかもしれません。誤りのご指摘・削除のご依頼をいただければ直します。 連絡先：support@tecktom.com
