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

前記事の続きです。

次に(10),(11),(12)の∂ｐ/∂ｘ,∂ｐ/∂ｙ,∂ｐ/∂ｚの表現を用いてｖ(ｒ)＝∫_r'=a^r'=∞Ｇ(ｒ,ｒ')gradｐ(ｒ')ｄ³ｒ'＋ｉηＵ∫_r'=∞∂Ｇ(ｒ,ｒ')/∂ｒ'ｄＳ' (8) の右辺第一項を求めます。

ｘ成分は,ｖ_x－Ｕ(1－ａ/ｒ)＝{1/(2η)}Σ_l=0^∞{2(ｌ＋1)Ａ_lＰ_l+1(cosθ)/(2ｌ＋3)}[∫₀^r{ｒ'^l+1/ｒ^l+2－ａ^2l+1/(ｒｒ')^l+2}(1/ｒ'^l)ｄｒ'＋∫_r^∞{ｒ^l+1/ｒ'^l+2－ａ^2l+1/(ｒｒ')^l+2}(1/ｒ'^l)ｄｒ']です。

この積分を実行すると,ｖ_x＝Ｕ(1－ａ/ｒ)＋{1/(2η)}Σ_l=0^∞{(ｌ＋1)Ａ_lＰ_l+1(cosθ)/(2ｌ＋1)}{(ｒ²－ａ²)/ｒ^l+2} (15) となります。

同様に,ｖ_y＝{1/(2η)}Σ_l=0^∞{(ｌ＋1)Ａ_lＰ¹_l+1(cosθ)cosφ/(2ｌ＋1)}{(ｒ²－ａ²)/ｒ^l+2}(16),およびｖ_z＝{1/(2η)}Σ_l=0^∞{(ｌ＋1)Ａ_lＰ¹_l+1(cosθ)sinφ/(2ｌ＋1)}{(ｒ²－ａ²)/ｒ^l+2}(17)を得ます。

圧力ｐの一般形:ｐ＝ｐ_∞＋Σ_l=0^∞(Ａ_l/ｒ^l+1)Ｐ_l(cosθ) (9)と,今求めたばかりのｖに対する(15),(16),(17)の表式により,境界条件(4),(5)と△ｐ＝0 (3)を満たす圧力ｐとポアソン方程式:η△ｖ＝gradｐ (1)の解として流速ｖが得られたことになります。

しかし,これらが連続の方程式:divｖ＝∂ｖ_x/∂ｘ＋∂ｖ_y/∂ｙ＋∂ｖ_z/∂ｚ＝0 (2)を満たすかどうかはわかりません。

そこで,具体的に(15),(16),(17)のdivｖを計算し,それをゼロとおくことによって係数Ａ_lの満足すべき条件を求めます。

泥臭い計算によって,∂ｖ_x/∂ｘ＝ａＵcosθ/ｒ²＋{1/(2η)}Σ_l=0^∞(ｌ＋1)Ａ_l[{－ｌｒ²＋(ｌ＋2)ａ²}cosθＰ_l+1(cosθ)＋(ｒ²－ａ²)sin²θ{ｄＰ_l+1(cosθ)/ｄ(cosθ)}]/{(2ｌ＋1)ｒ^l+3},

∂ｖ_y/∂ｙ＝{1/(2η)}Σ_l=0^∞Ａ_l[{－ｌｒ²＋(ｌ＋2)ａ²}sinθcosφＰ¹_l+1(cosθ)＋(ｒ²－ａ²)(cosθsinθcosφ{ｄＰ¹_l+1(cosθ)/ｄ(cosθ)}－(sinφ/sinθ)Ｐ¹_l+1(cosθ))]/{(2ｌ＋1)ｒ^l+3},

∂ｖ_z/∂ｚ＝{1/(2η)}Σ_l=0^∞Ａ_l[{－ｌｒ²＋(ｌ＋2)ａ²}sinθsinφＰ¹_l+1(cosθ)＋(ｒ²－ａ²)(cosθsinθsinφ{ｄＰ¹_l+1(cosθ)/ｄ(cosθ)}－(cosφ/sinθ)Ｐ¹_l+1(cosθ))]/{(2ｌ＋1)ｒ^l+3}

が得られます。

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

すると,divｖ＝ａＵcosθ/ｒ²＋(1/η)Σ_l=0^∞Ａ_l[sin²θ{ｄＰ_l+1(cosθ)/ｄ(cosθ)}＋(ｌ＋1)cosθＰ_l+1(cosθ)]/{(2ｌ＋1)ｒ^l+1}となります。

さらに,公式(1－ｚ²){ｄＰ_l(ｚ)/ｄｚ}＋ｌｚＰ_l(ｚ)＝ｌＰ_l-1(ｚ)によって,sin²θ{ｄＰ_l+1(cosθ)/ｄ(cosθ)}＋(ｌ＋1)cosθＰ_l+1(cosθ)＝(ｌ＋1)Ｐ_l(cosθ)です。

結局,divｖ＝ａＵcosθ/ｒ²＋(1/η)Σ_l=0^∞(ｌ＋1)Ａ_lＰ_l(cosθ)/{(2ｌ＋1)ｒ^l+1}＝0 (18)を得ます。

これが成立するためには,2Ａ_l/(3η)＝－ａＵδ_l1,つまりＡ_l＝－(3/2)ηａＵδ_l1(19)となるのが必要十分です。

したがって,ｐ＝ｐ_∞－(3/2)ηａＵcosθ/ｒ²＝ｐ_∞－(3/2)ηａＵｘ/ｒ³ (20)です。

また,ｕ≡ｖ_x＝Ｕ(1－ａ/ｒ)－{ａＵ/(4ｒ)}(1－ａ²/ｒ²)(3cos²θ－1)＝Ｕ－{ａＵ/(4ｒ)}(3＋ａ²/ｒ²)－{3ａＵｘ²/(4ｒ³)}(1－ａ²/ｒ²)(21),ｖ≡ｖ_y＝－{3ａＵ/(4ｒ)}(1－ａ²/ｒ²)cos²θsinθcosφ＝－{3ａＵｘｙ/(4ｒ³)}(1－ａ²/ｒ²) (22),ｗ≡ｖ_z＝－{3ａＵ/(4ｒ)}(1－ａ²/ｒ²)cos²θsinθsinφ＝－{3ａＵｘｚ/(4ｒ³)}(1－ａ²/ｒ²) (23)も得られます。

かくして,方程式(1)～(5)は完全に解かれ,解は一意的であることも示されました。

方程式が解かれたので,この解を利用して流体が球に及ぼす力を求めることができます。

この力は球面上に働く応力(単位面積当たりに働く力;法線および接線応力)を小球面にわたって積分することで求めることができます。

まず,一般的な応力テンソルＰ_ijを考えます。運動がナビエ・ストークスの方程式に従う流体は,応力テンソルが次のニュートン・ストークスの法則に従う線形ニュートン流体です。

すなわち,応力テンソルＰ_ijはＰ_ij＝(－ｐ＋ζｅ_ll)δ_ij＋2η(ｅ_ij－ｅ_llδ_ij/3) (ｉ,ｊ＝1,2,3はそれぞれｘ,ｙ,ｚ成分に対応) (24)で与えられます。

ここで,ｐは圧力,η,ζはそれぞれ第一粘性率,第二粘性率であり,ｅ_ij≡(∂ｖ_i/∂ｘ_j＋∂ｖ_j/∂ｘ_i)/2は歪み速度テンソルです。

非圧縮性流体では連続の方程式としてdivｖ＝ｅ_lｌ＝0 が成立するので,(24)はＰ_ij＝－ｐδ_ij＋η(∂ｖ_i/∂ｘ_j＋∂ｖ_j/∂ｘ_i) (25)となります。

応力ｆと応力テンソルＰ_ijの関係はｆ_i＝Ｐ_ijｎ_j(ｊはアインシュタインの規約に従うダミー添字)で与えられます。

また,ｎは応力ｆの働くべき面の外向き法線ベクトルで,小球面上ではｎ≡(ｎ_x,ｎ_y, ｎ_z)＝ｒ/ｒ＝(ｘ/ｒ,ｙ/ｒ,ｚ/ｒ)_r=a＝(cosθ,sinθcosφ,sinθsinφ)と表わされます。

小球全体に働く力Ｄは応力ｆによってＤ＝∫ｆｄＳで与えられます。ところが流れはｘ軸に関して対称(ｕはｙ,ｚの偶関数),ｖ,ｗはｙ,ｚの奇関数)ですから,∫ｆ₂ｄＳ＝∫ｆ₃ｄＳ＝0 となるはずです。

実際,ｆ₂＝Ｐ_2jｎ_j,ｆ₃＝Ｐ_3jｎ_jなので,解からｆ₂(－ｎ)＝－ｆ₂(ｎ),ｆ₃(－ｎ)＝－ｆ₃(ｎ)あるいはＰ_2j,Ｐ_3jがｎ_jの偶関数であることが示されます。

まず,－ｐδ_2j,－ｐδ_3jですが,ｐ＝ｐ_∞－(3/2)ηａＵｘ/ｒ³は明らかにｙ,ｚの偶関数であり,よって－ｐδ_2j,－ｐδ_3jはｙ,ｚの偶関数です。さらにｎ_jのｊ＝1,つまりｘについてはこれらはゼロなのでもちろん偶関数です。

また,ｖ＝ｖ_y＝ｖ₂,ｗ＝ｖ_z＝ｖ₃はｘ,ｙ,ｚの奇関数ですから明らかに∂ｖ₂/∂ｘ_j,∂ｖ₃/∂ｘ_jはｎ_jの偶関数です。

さらに,∂ｖ_j/∂ｘ₂,∂ｖ_j/∂ｘ₃についてはｊ＝1だけの問題になりますが,ｕはｘの偶関数ですから∂ｕ/∂ｙ,∂ｕ/∂ｚもｎ_jの偶関数です。以上からｆ₂(－ｎ)＝－ｆ₂(ｎ),ｆ₃(－ｎ)＝－ｆ₃(ｎ)が確証されました。

したがって,Ｄはｘ成分のみ,つまり抵抗成分のみを有すると考えられます。

これをＤで表わすと,Ｄ＝∫Ｐ_1jｎ_jｄＳ＝∫Ｐ_rxｄＳ；Ｐ_rx≡Ｐ_1jｎ_j＝Ｐ₁₁(ｘ/ｒ)＋Ｐ₁₂(ｙ/ｒ)＋Ｐ₁₃(ｚ/ｒ)＝－ｐcosθ＋η[(ｘ_j/ｒ)(∂ｖ_x/∂ｘ_j)＋(ｘ_j/ｒ)(∂ｖ_j/∂ｘ₁)] (26)となります。

そして,∂ｕ/∂ｒ＝{3ａＵ/(4ｒ²)}[(1＋cos²θ)＋ａ²(1－3cos²θ)/ｒ²],∂ｕ/∂θ＝{3ａＵcosθsinθ/(2ｒ)}(1－ａ²/ｒ²)です。

一方,∂/∂ｘ＝cosθ(∂/∂ｒ)－(sinθ/ｒ)(∂/∂θ),∂/∂ｙ＝sinθcosφ(∂/∂ｒ)＋(cosθcosφ/ｒ)(∂/∂θ)－{sinφ/(ｒsinθ)}(∂/∂φ),∂/∂ｚ＝sinθsinφ(∂/∂ｒ)＋(cosθsinφ/ｒ)(∂/∂θ)＋{cosφ/(ｒsinθ)}(∂/∂φ)です。

そこで,∂ｕ/∂ｘ＝{3ａＵcosθ/(4ｒ²)}[(3cos²θ－1)＋ａ²(3－5cos²θ)/ｒ²],∂ｕ/∂ｙ＝{3ａＵsinθcosφ/(4ｒ²)}[(3cos²θ＋1)＋ａ²(1－5cos²θ)/ｒ²],∂ｕ/∂ｚ＝{3ａＵsinθsinφ/(4ｒ²)}[(3cos²θ＋1)＋ａ²(1－5cos²θ)/ｒ²]となります。

それ故,(ｘ_j/ｒ)(∂ｖ_x/∂ｘ_j)＝{3ａＵ/(4ｒ²)}[(1＋cos²θ)＋ａ²(1－3 cos²θ)/ｒ²]＝∂ｕ/∂ｒ (27)が得られます。

同様に,∂ｖ/∂ｘ＝{3ａＵsinθcosφ/(4ｒ²)}[(3cos²θ－1)＋ａ²(1－5cos²θ)/ｒ²],∂ｗ/∂ｘ＝{3ａＵsinθsinφ/(4ｒ²)}[(3cos²θ－1)＋ａ²(1－5cos²θ)/ｒ²] (28)より,(ｘ_j/ｒ)(∂ｖ_j/∂ｘ₁)＝{3ａＵ/(4ｒ²)}[(3cos²θ－1)＋ａ²(1－5cos²θ)/ｒ²]です。

(20),(25),(27),(28)から,ｆ₁＝ｐ_1jｎ_j＝－ｐ_∞cosθ＋{(3/2)ηａＵ/ｒ²)}[3cos²θ＋ａ²(1－3cos²θ)/ｒ²] (29)です。特に,球面上(ｒ＝ａ)では,ｆ₁|r=a＝－ｐ_∞cosθ＋(3/2)ηＵ/ａ＝(1/ａ){－ｐ_∞ｘ＋(3/2)ηＵ}r=a (30)です。

結局,Ｄ＝∫r=aｆ₁ｄＳ＝2πａ²∫_-1¹ｄ(cosθ)ｆ₁＝6πηａＵ (31)を得ました。すなわち,球は流れの向きに大きさ6πηａＵの力を受けるという有名なストークスの法則が得られました。

実は,これほど苦労して偏微分方程式の境界値問題としてマニュアル通りに解かなくても,Lambや今井功氏の流体力学の教科書にもあるように,さまざまな物理的考察や手法を交えた工夫を用いることよって,より簡単に解を求めることができます。

しかし,当時は物理的な意味などをあまり考えることなく,なるべく式に忠実に数理的に解いても同じ答えがちゃんと得られるか？に興味があって,あえて泥臭い計算をしたのでした。

こうした偏微分方程式の境界値問題をマニュアル通りに解くという泥臭い計算は,その後確か,気象学関係で雲などによる光の散乱の散乱断面積を計算する必要があって,実行しこともありました。

すなわちり,ボルン,ウォルフ著の「光学の原理」にもある"光＝電磁波"の古典的なミイ(Mie)散乱の散乱断面積の式を確認する際にもやったという記憶があります。

ところで,最後にストークスの法則のガリレイの相対性原理との関係を記述しておきます。

すなわち,先の"流体＋小球"の系を一様な相対速度ｕで動く並進座標系から見た場合を考えます。

元の座標系での位置ベクトルをｒ,新しい系での同じ点の位置ベクトルをｒ'で表わします。

便宜上ｔ＝0 で両系が一致すると仮定するとｒとｒ'の関係はガリレイ変換:ｒ'＝ｒ－ｕｔ (32)で与えられます。

位置ｒに存在している流体素片の速度ｖ(ｒ)はラグランジュ微分Ｄｒ/ｄｔに等しいので,ｖ'(ｒ')≡Ｄｒ'/ｄｔ＝Ｄｒ/ｄｔ－Ｄ(ｕｔ)/ｄｔ＝ｖ(ｒ)－ｕ,すなわち,ｖ'(ｒ')＝ｖ(ｒ)－ｕ (33)です。

一方,小球は元の系では静止していたので,新しい系では－ｕなる速度を有します。

また,ｘ'_j＝ｘ_j－ｕ_jｔより∂/∂ｘ_j＝(∂ｘ'_k/∂ｘ_j)∂/∂ｘ'_k ＋(∂ｔ/∂ｘ_j)∂/∂ｔ＝∂/∂ｘ_jです。それ故,∂²/∂ｘ_k∂ｘ_j＝∂²/∂ｘ'_k∂ｘ'_jです。特に△＝△'かつgrad＝grad'です。

さらに流速のオイラー微分は,∂ｖ(ｒ,ｔ)/∂ｔ＝∂ｖ(ｒ'＋ｕｔ,ｔ)/∂ｔ－(ｕ・grad')ｖ(ｒ'＋ｕｔ,ｔ)となります。

そこで,ナビエ・ストークスの方程式:ρ(∂ｖ/∂ｔ)＋ρ(ｖ・grad)ｖ＝－gradｐ＋η△ｖを,新しい座標系で見ると,ρ{∂(ｖ－ｕ)/∂ｔ}＋ρ{(ｖ－ｕ)・grad'}(ｖ－ｕ)＝－grad'ｐ＋η△'(ｖ－ｕ)です。

ところで,先に見たようにｖ'＝ｖ－ｕです。また,圧力ｐ(ｒ,ｔ)は運動学的な量ではないのでガリレイ変換のスカラーです。すなわちｐ'(ｒ',ｔ)＝ ｐ(ｒ,ｔ)です。

結局,ρ(∂ｖ'/∂ｔ)＋ρ(ｖ'・grad')ｖ'＝－grad'ｐ'＋η△'ｖ'となります。

よって,新しい系でもナビエ・ストークスの方程式の形は保存されます。また,連続の方程式が保存されること:div'ｖ'＝0 も自明です。

したがって,ガリレイ変換(32),(33)に対してナビエ・ストークスの方程式が保存されるので,ニュートン流体に対してはガリレイの相対性原理が成立します。

また,境界条件が与えられた場合にもガリレイの相対性原理が無矛盾であることは自明です。

それ故,ストークス近似の流体方程式の組:η△ｖ＝gradｐ (1),divｖ＝0 (2) において,無限遠で静止した流体中を一様速度－ｕ＝－Ｕｉで運動する小球がある系に対する流体速度の解は,既に求めた(21)～(23)式のｖに－ｕを加えたもので与えられることがわかります。

以上の考察から,静止流体中を速度Ｕｉで運動する小球があるような系に対して,境界条件:ｖ＝Ｕｉ at ｒ＝ａ (34),ｖ＝0 ,ｐ＝ｐ_∞ at ｒ＝∞ (35)を満足する(1),(2)の解は次式で与えられます。

すなわち,ｐ＝ｐ_∞－(3/2)ηａＵｘ/ｒ³ (36),また,ｕ≡ｖ_x＝{ａＵ/(4ｒ)}(3＋ａ²/ｒ²)＋{3ａＵｘ²/(4ｒ³)}(1－ａ²/ｒ²) (37),ｖ≡ｖ_y＝{3ａＵｘｙ/(4ｒ³)}(1－ａ²/ｒ²) (38),ｗ≡ｖ_z＝{3ａＵｘｚ/(4ｒ³)}(1－ａ²/ｒ²) (39)です。また,Ｄ＝－6πηａＵ (40)も自明です。

つまり小球は進行方向に正反対の向きに大きさ6πηａＵの抵抗を受けることになります。

ここで,レイノルズ数:Ｒ＝ρＵＬ/η＝ＵＬ/νにおいて通常は代表長さＬとして小球の直径ｄ＝2ａを選びます。

このとき無次元化された抵抗係数Ｃ_D≡Ｄ/(ρＵ²Ｓ/2)＝Ｄ/(ρＵ²πａ²/2)をストークスの抵抗法則で表わすとＣ_D＝24/Ｒとなります。

しかし,こうしたストークス近似が可能で抵抗が運動速度に比例するというストークスの抵抗法則が有効なレイノルズ数Ｒの範囲は非常に狭く,この法則はせいぜいＲ～100前後でしか成立しません。

通常の雨滴落下や冷却塔の温排気の凝結した水滴落下の終末速度の計算などにおいても,経験上はほとんど適用不可能であることがわかっています。

次は,もう少しはましな近似であるオセーン近似(Oseen近似)の計算をする予定です。

http://folomy.jp/heart/「folomy 物理フォーラム」サブマネージャーです。

人気blogランキングへ　← クリックして投票してください。(１クリック＝１投票です。１人１日１投票しかできません。）

http://homepage2.nifty.com/toshis-kaiga-auction/健康商品の店 「TRS健康ランド」

← クリックして投票してください。(ブログ村科学ブログランキング）

    物理学