遅い粘性流(１)(Stokes近似)

学生時代は素粒子論が専門でしたが,就職難の1977年に流体の数値シミュレーションが主体の会社に入ってからは,必要に迫られて流体力学と数値計算の勉強をするようになり,これが第２の専門となりました。

応用物理学はやったことがなかったし,仕事もひととおりの公式類を知っていれば特に支障はなかったのですが,持ち前の好奇心もあって入社して数年間は流体力学関係の勉強をしました。

もっとも,当時はプログラマー残酷時代といわれる時代で,残業も半端ではなく,暇があると酒を飲んで発散していました。

それでも,寝る前などに細々と勉強して書き綴ったのが,この記事のネタ本である1980年当時の覚書きノートです。

もっとも,私は,昔から勉強が好きだったというわけではありませんでした。

最近,TOKIOの山口達也君が父親役をやっているドラマで,"受験の神様"と呼ばれる女子中学生に「あなたは勉強したいの,それとも受験がしたいの？」と問われてキョトンとしているシーンがありましたが,「勉強＝受験」であった頃は成績は悪くなかったけど,強いられる感じの勉強は好きではありませんでした。 

試験の点数がいいとか速く解けるか?ということで他人との競争に負けるのは大して気にならなかったし,答案が不完全とかケアレスミスなども二の次でしたが,いざ自分が取り掛かった問題が解けないとか理解できないということは我慢できない性格でした。

試験ではそうはいかないでしょうが,わからないときは持ち帰って,ちゃんと解けて,理解できるまでは粘ったという記憶があります。 

未だに「高々知っているか知らないかの違いじゃないか」というような記憶するだけの知識レベルについては,あまり興味はないのですが,自然現象に限らず社会現象についても,なぜそうなるのか？というような合理性を追求することについては,理解できて解決できないと我慢できないという傾向は少なからず残っています。

その他,言わせてもらえば,語学のような学問は,会話で英語などの外国語とかが読めたり話せたりするだけであれば,その国のネイティブな幼児のほうがより優れていると思えますから,バイリンガルであることにこそ,その意味があると思っています。

社会人になってみると,無理に勉強しなくてもよいし,むしろ勉強や研究をする機会が少なくなりましたが,逆に好奇心だけでいろいろとやってみようという気になりました。 

余談はさておき,流体運動の現象というのはこの地球上の表面では薄い層を除いて大部分の領域では粘性が無視できるレイノルズ数(Reynolds number)の大きい流れが主体です。

こうした領域では完全流体の近似が成り立ちます。そして多くは乱流に移行する現象になっています。

しかし,ここでは,1980年当時の覚書きノートから,非圧縮かつ密度が一様な単成分流体のナビエ・ストークスの方程式(Navier-Stokes equation)において,レイノルズ数が比較的小さくストークス(Stokes)近似やオセーン(Oseen)近似が有効なケースについて流体中の球状物体に及ぼされる力＝抵抗力を計算してみます。

今日はまず古典的なストークス近似の扱いからはじめます。

非圧縮,かつ密度が一様な単成分ニュートン流体のナビエ・ストークスの方程式はρ(∂ｖ/∂ｔ)+ρ(ｖ・grad)ｖ＝－gradｐ＋η△ｖで与えられます。

また,この場合の"連続の方程式＝流体質量保存)の方程式"はdivｖ＝0 です。ここでρは流体の密度,ｖ＝ｖ(ｒ,ｔ)は流体速度,ｐは圧力,ηは第一粘性率です。

そして,レイノルズ数をＲ≡ρＵＬ/η＝ＵＬ/νで定義します。ν＝η/ρは動粘性係数と呼ばれ,ＵとＬは,それぞれ流れを特徴付ける代表流速と代表長さです。

レイノルズ数:Ｒが小さい場合は遅い粘性流として,ナビエ・ストークスの方程式の流速の２乗程度の"左辺第２項＝慣性項"は非常に小さいので,"右辺の項＝粘性項"と比較して省略する近似が可能で,これをストークス近似といいます。

特に,流速ｖが時刻ｔに依らない定常流に対する定常ストークス近似の問題は,η△ｖ＝gradｐ (1)とdivｖ＝0 (2)を連立させて解く問題になります。

まず,(1)式の発散を取り(2)を適用すれば,△ｐ＝0 (3)を得ますからｐは調和関数であることがわかります。

無限遠で,ｖが一様速度ｕに,圧力ｐがｐ_∞(一定)となる非圧縮性流体中に小半径ａの球を静止固定させた系において,Ｕ≡|ｕ|が小さい低レイノルズ数の線型粘性流に対し定常ストークス近似を適用します。

ｕの向きをｘ軸の正の向きに取り,また小球の中心を原点にとります。さらに極座標(ｒ,θ,φ)に対してはｘ軸を極軸に取ります。

境界条件は次のように与えられます。ｖ＝0 at ｒ＝|ｒ|＝ａ (4),ｖ＝Ｕi,ｐ＝ｐ_∞　at ｒ＝∞ (5)です。

△ｐ＝0 (3)の解ｐ(ｒ)を(1)に代入すると,(4),(5)の境界条件に対して(1),(2)の解を求める問題は,右辺が既知となったポアソン(Poisson)方程式に対する条件付きのディリクレ問題となります。

ここで,η△Ｇ(ｒ,ｒ')＝δ(ｒ－ｒ') (6),Ｇ(ｒ,ｒ')＝0 at ｒ＝ａ, ∞ (7) によって演算子η△に対する第一種のグリーン関数Ｇ(ｒ,ｒ')を定義すると,(1)の解は次式で与えられます。

ｖ(ｒ)＝∫_r'=a^r'=∞Ｇ(ｒ,ｒ')grad'ｐ(ｒ')ｄ³ｒ'＋ｉηＵ∫_r'=∞{∂Ｇ(ｒ,ｒ')/∂ｒ'}ｄＳ' (8)です。

なぜなら,グリーンの公式:∫_Ｖ{ｖ(ｒ)△Ｇ(ｒ,ｒ')－Ｇ(ｒ,ｒ')△ｖ(ｒ)}ｄ³ｒ＝∫_Ｓ{ｖ(ｒ)∂Ｇ(ｒ,ｒ')/∂ｎ－Ｇ(ｒ,ｒ')∂ｖ(ｒ)/∂ｎ}ｄＳにおいて,両辺にηを乗じη△ｖ(ｒ)＝gradｐ(ｒ),η△Ｇ(ｒ,ｒ')＝δ(ｒ－ｒ')を代入すれば,Ｓ上ではＧ(ｒ,ｒ')＝0 なのでｖ(ｒ')＝∫_ＶＧ(ｒ,ｒ')gradｐ(ｒ)ｄ³ｒ＋η∫_Ｓ{∂Ｇ(ｒ,ｒ')/∂ｎ}ｄＳとなるからです。

ここで,積分変数をｒ'としてｒとｒ'を交換すればｖ(ｒ)＝∫_Ｖ'Ｇ(ｒ',ｒ)gradｐ(ｒ')ｄ³ｒ'＋η∫_Ｓ'{∂Ｇ(ｒ',ｒ)/∂ｎ'}ｄＳ'です。

さらにグリーン関数の対称性からＧ(ｒ',ｒ)＝Ｇ(ｒ,ｒ')であり,よって明らかに∂Ｇ(ｒ',ｒ)/∂ｎ'＝∂Ｇ(ｒ,ｒ')/∂ｎ'ですから,結局ｖ(ｒ)＝∫_Ｖ’Ｇ(ｒ,ｒ')grad'ｐ(ｒ')ｄ³ｒ'＋η∫_Ｓ'{∂Ｇ(ｒ,ｒ')/∂ｎ'}ｄＳ'が得られます。

ただし,Ｖ'はｒ'＝ａとｒ'＝∞で囲まれる同心球の領域はｒ'＝ａとｒ'＝∞ の同心球面を示しており,∂/∂ｎ'はｒ'＝ａでは内向き,ｒ'＝ ∞では外向きの法線方向の微分係数を示しています。

したがって(6),(7)を満足するＧ(ｒ,ｒ')を求めれば,解ｖ(ｒ)が得られることになります。

そこで,まず,△ｐ＝0.を満足する解ｐ(ｒ)を求めます。

系の極軸に対する対称性からｐはｒとcosθのみの球関数によって表わされます。

ｒ→ ∞でｐ＝ｐ_∞となる解の一般形はｐ＝ｐ_∞＋Σ_l=0^∞(Ａ_l/ｒ^l+1)Ｐ_l(cosθ) (9)で与えられます。ここでＰ_l(ｚ)はルジャンドル多項式(Legendre多項式)です。

次に,gradｐ＝ｉ(∂ｐ/∂ｘ)＋ｊ(∂ｐ/∂ｙ)＋ｋ(∂ｐ/∂ｚ)を極座標に変換します。

ｘ軸を極軸としているので,∂/∂ｘ＝cosθ(∂/∂ｒ)＋(sin²θ/ｒ){∂/∂(cosθ)},∂/∂ｙ＝sinθcosφ(∂/∂ｒ)－(cosθsinθcosφ/ｒ){∂/∂(cosθ)}－{sinφ/(ｒsinθ)}(∂/∂φ),∂/∂ｚ＝sinθsinφ(∂/∂ｒ)－(cosθsinθsinφ/ｒ){∂/∂(cosθ)}－{cosφ/(ｒsinθ)}(∂/∂φ)です。

故に,∂ｐ/∂ｘ＝－Σ_l=0^∞{(ｌ＋1)Ａ_l/ｒ^l+2}Ｐ_l+1(cosθ) (10),∂ｐ/∂ｙ＝－Σ_l=0^∞(Ａ_l/ｒ^l+2)Ｐ¹_l+1(cosθ)cosφ (11) です。

ただし,Ｐ¹_l+1(cosθ)≡sinθ{ｄＰ_l+1(cosθ)/ｄ(cosθ)}＝{(ｌ＋1)/sinθ}{Ｐ_l(cosθ)－cosθＰ_l+1(cosθ)}を用いています。

また,∂ｐ/∂ｚ＝－Σ_l=0^∞{(ｌ＋1)Ａ_l/ｒ^l+2}Ｐ¹_l+1(cosθ)sinφ (12) です。

(実際のノートには泥臭い計算の詳細があるのですが,ここでは割愛しています。)

次に(6),(7)を満足する第一種グリーン関数Ｇ(ｒ,ｒ')を求めます。

良く知られているように,△(1/|ｒ－ｒ'|)＝－4πδ(ｒ－ｒ')です。

そこで,ｒ≧ａの全領域でｒ＝ｒ'を除いて特異点を持たず有限確定値を取るポアソン方程式の任意のグリーン関数は,－1/(4π|ｒ－ｒ'|)にラプラス方程式の一般解を加えることによって得られます。

ただし,この場合の極軸(対称軸)はｒとｒ'を結ぶ有向直線です。

これはグリーン関数の対称性によるものです。そこで,Ｇ(ｒ,ｒ')＝{－1/(4πη)}{1/|ｒ－ｒ'|－Σ_l=0^∞(ｃ_l/ｒ^l+1)Ｐ_l(cosω)}(ｒ≧ａ)と書けます。

ωはｒとｒ'のなす角です。ｒ＝(ｒ,θ,φ),ｒ'＝(ｒ',θ',φ')なら,cosω＝cosθcosθ'＋sinθsinθ'cos(φ－φ'),cosω＝(ｒ,ｒ')/(|ｒ||ｒ'|)が成立します。

また,式Ｐ_l(cosω)＝Ｐ_l(cosθ)Ｐ_l(cosθ')＋2Σ_m=1^l{(ｌ－ｍ)!/(ｌ＋ｍ)!}Ｐ^m_l(cosθ)Ｐ^m_l(cosθ')cos{ｍ(φ－φ')}も成立します。

1/|ｒ－ｒ'|＝1/(ｒ²＋ｒ'²－2ｒｒ'cosω)^1/2をＰ_l(cosω)で展開すると,ｒ＝ａでＧ(ｒ,ｒ')＝0 となるような展開係数ｃ_lがＰ_l(ｚ)の直交性から一意的に定まります。

これは静電気学における鏡像の原理と同じ結果を与えます。つまりアポロニウスの定理により,Ｇ(ｒ,ｒ')の表現式において右辺第二項の級数和－Σ_l=0^∞(ｃ_l/ｒ^l+1)Ｐ_l(cosω)はｒ'と同一の向きで半径がａ²/ｒ'の位置に大きさ－ａ/ｒ'の源を持つニュートン・ポテンシャルに一致します。

すなわち,Ｇ(ｒ,ｒ')＝{－1/(4πη)}[1/(ｒ²＋ｒ'²－2ｒｒ'cosω)^1/2－(ａ/ｒ'){1/(ｒ²＋ａ⁴/ｒ'²－2ｒａ²/ｒ'cosω)^1/2}]です。

そこで,Ｇ(ｒ,ｒ')＝{－1/(4πη)}Σ_l=0^∞{ｒ^l/ｒ'^l+1－ａ^2l+1/(ｒｒ')^l+1}Ｐ_l(cosω)(ａ≦ｒ＜ｒ'),Ｇ(ｒ,ｒ')＝{－1/(4πη)}Σ_l=0^∞{ｒ'^l/ｒ^l+1－ａ^2l+1/(ｒｒ')^l+1}Ｐ_l(cosω) (ａ≦ｒ'＜ｒ) (13) となります。

(13)によって,(6),(7)を満足するグリーン関数が完全に決定されたのでｖ(ｒ)＝∫_r'=a^r'=∞Ｇ(ｒ,ｒ')gradｐ(ｒ')ｄ³ｒ'＋ｉηＵ∫_r'=∞{∂Ｇ(ｒ,ｒ')/∂ｒ'}ｄＳ' (8)に従って速度ｖ(ｒ)を求めます。

まず,右辺第二項,すなわちｖ_xの剰余的な項を計算します。 

以下の計算では,Ｐ_l(ｚ)の直交性∫_-1¹ｄｚＰ_l(ｚ)Ｐ_l'(ｚ)＝2δ_ll'/(2ｌ＋1)や,∫_-1¹ｄｚ{Ｐ^m_l(ｚ)}²＝{2/(2ｌ＋1)}{(ｌ＋ｍ)!/(ｌ－ｍ)!},∫₀^2πｄφ'cos{ｍ(φ－φ')}＝0 ,∫₀^2πｄφ'cos{ｍ(φ－φ')}cos{ｍ'(φ－φ')}＝∫₀^2πｄφ'sin{ｍ(φ－φ')}sin{ｍ'(φ－φ')}＝πδ_mm'等の公式を用います。

ηＵ∫_r'=∞{∂Ｇ(ｒ,ｒ')/∂ｒ'}ｄＳ'において,∂Ｇ(ｒ,ｒ')/∂ｒ'を考えると(13)の右辺の級数展開で,この積分への寄与がゼロでないものはｌ＝0 に対応する項のみなので次式が得られます。

すなわち,ηＵ∫_r'=∞{∂Ｇ(ｒ,ｒ')/∂ｒ'}ｄＳ'＝－Ｕ・lim_R→∞[Ｒ²{－1/Ｒ²＋ａ/(ｒＲ²)}]＝Ｕ(1－ａ/ｒ)です。

途中ですが疲れたので今日はここまでにします。