水の波(8)(有限振幅の波:非線形波３)

　水の波の続きです。

　このシリーズは今日で終わりますが,非線形波やソリトンについての話はまた別の題名でアップする予定です。

　前回の最後では,Ｋ-ｄＶ方程式:∂ｕ/∂ｔ＋ｕ(∂ｕ/∂ｘ)＋μ(∂³ｕ/∂ｘ³)＝0 の定常波の解をｕ(ｘ,ｔ)＝ｕ(ζ),ζ≡ｘ－σｔとして３階常微分方程式：－σ(ｄｕ/ｄζ)＋ｕ(ｄｕ/ｄζ)＋μ(ｄ³ｕ/ｄζ³)＝0 を求め,これの積分を求めました。

すなわち,１回積分すると,μ(ｄ²ｕ/ｄζ²)＝－ｕ²/2＋σｕ＋2Ａ (Ａは積分定数)となり,さらに両辺にｄｕ/ｄζを掛けて,積分すると(μ/2)(ｄｕ/ｄζ)²＝－ｕ³/6＋σｕ²/2＋Ａｕ＋Ｂ(Ａ,Ｂは積分定数)となります。

両辺を６倍して右辺をｆ(ｕ)とすると,3μ(ｄｕ/ｄζ)²＝－ｕ³＋3σｕ²＋6Ａｕ＋6Ｂ≡ｆ(ｕ)です。この方程式が物理的に意味があるｕの実数解を持つのは,明らかにｆ(ｕ)≧0 の範囲でのみです。

実係数の３次方程式ｆ(ｕ)＝0 は１実根,または３実根を持ちますが,この方程式が１実根のみを持つ場合には,ｆ(ｕ)≧0 を値域とする解ｕが有界でないため,これは"有限範囲での振動＝波"を表わす解ではないので,ここでの考察では対象外とします。

ｆ(ｕ)＝0 が３実根を持つとして,それらをｕ₁,ｕ₂,ｕ₃ (ｕ₁≦ｕ₂≦ｕ₃)と書くと,ｆ(ｕ)＝(ｕ－ｕ₁)(ｕ－ｕ₂)(ｕ₃－ｕ)です。

ｆ(ｕ)＝－ｕ³＋3σｕ²＋6Ａｕ＋6Ｂより,σ＝(ｕ₁＋ｕ₂＋ｕ₃)/3,Ａ＝－(ｕ₁ｕ₂＋ｕ₂ｕ₃＋ｕ₃ｕ₁)/6,Ｂ＝ｕ₁ｕ₂ｕ₃ /6です。

ｕ₂≦ｕ≦ｕ₃でｆ(ｕ)≧0ですから,3μ(ｄｕ/ｄζ)²＝ｆ(ｕ)はｕのこの区間を往復する非線形振動を表わすと思われます。

　そこで,3μ(ｄｕ/ｄζ)²＝ｆ(ｕ)＝(ｕ－ｕ₁)(ｕ－ｕ₂)(ｕ₃－ｕ)の解でｕ＝ｕ₃でζ＝0 となるものを求めます。

　つまり,方程式ｄζ/(3μ)^1/2＝ｄｕ/{(ｕ－ｕ₁)(ｕ－ｕ₂)(ｕ₃－ｕ)}^1/2を解きます。

　まず,ｋ²≡(ｕ₃－ｕ₂)/(ｕ₃－ｕ₁)とします。さらに,ｔ＝{(ｕ₃－ｕ)/(ｕ₃－ｕ₂)}^1/2とおくとｄｔ＝(－1/2)(ｕ₃－ｕ₂)^-1/2(ｕ₃－ｕ)^-1/2ｄｕです。

1－ｔ²＝1－(ｕ₃－ｕ)/(ｕ₃－ｕ₂)＝(ｕ－ｕ₂)/(ｕ₃－ｕ₂),1－ｋ²ｔ²＝1－ｋ²(ｕ₃－ｕ)/(ｕ₃－ｕ₂)＝1－(ｕ₃－ｕ)/(ｕ₃－ｕ₁)＝(ｕ－ｕ₁)/(ｕ₃－ｕ₁)です。

故に,ｄｔ/{(1－ｔ²)(1－ｋ²ｔ²)}^1/2＝(－1/2)(ｕ₃－ｕ₂)^-1/2(ｕ₃－ｕ)^-1/2ｄｕ(ｕ₃－ｕ₂)^1/2(ｕ₃－ｕ₁)^1/2/{(ｕ－ｕ₁)(ｕ－ｕ₂)}^1/2＝(－1/2)(ｕ₃－ｕ₁)^1/2ｄｕ/{(ｕ－ｕ₁)(ｕ－ｕ₂)(ｕ₃－ｕ)}^1/2です。

そこで,方程式ｄζ/(3μ)^1/2＝ｄｕ/{(ｕ－ｕ₁)(ｕ－ｕ₂)(ｕ₃－ｕ)}^1/2は,－{(ｕ₃－ｕ₁)/(12μ)}^1/2ｄζ＝ｄｔ/{(1－ｔ²)(1－ｋ²ｔ²)}^1/2を意味します。

したがって,楕円積分Ｆ(ｘ,ｋ)を用いて{(ｕ₃－ｕ₁)/(12μ)}^1/2ζ＝Ｆ({(ｕ₃－ｕ)/(ｕ₃－ｕ₂)}^1/2,ｋ)なる解の表式を得ます。

それ故,sn({(ｕ₃－ｕ₁)/(12μ)}^1/2ζ,ｋ)＝{(ｕ₃－ｕ)/(ｕ₃－ｕ₂)}^1/2です。そこで,cn²({(ｕ₃－ｕ₁)/(12μ)}^1/2ζ,ｋ)＝1－sn²({(ｕ₃－ｕ₁)/(12μ)}^1/2ζ,ｋ)＝(ｕ－ｕ₂)/(ｕ₃－ｕ₂)です。

※(注):snはヤコービ(Jacobi)の楕円関数:sn(ｕ,ｋ)＝Ｆ^-1(ｕ,ｋ)＝ｘ≡sinφです。またはｕ＝Ｆ(ｘ,ｋ)です。

　Ｆ(ｘ,ｋ)は第１種の楕円積分で,Ｆ(ｘ,ｋ)≡∫₀^xｄｔ/{(1－ｔ²)(1－ｋ²ｔ²)}^1/2＝∫₀^φｄθ/(1－ｋ²sin²θ)^1/2で定義されます。

特に,Ｋ(ｋ)≡Ｆ(π/2,ｋ)は第１種の完全楕円積分と呼ばれ,4Ｋ(ｋ)は楕円関数:snの２つの周期のうちの１つです。

さらに,sn(ｕ,ｋ)＝ｘ＝sinφに対応して関数:cnをcn(ｕ,ｋ)≡cosφ＝{1－sn²(ｕ,ｋ)}^1/2で定義します。これもヤコービの楕円関数と呼ばれています。(注終わり)※

さて,振動区間をＵ≡ｕ₃－ｕ₂とおけば,上に得られた解はｕ(ｘ,ｔ)＝ｕ₂＋Ｕcn²({Ｕ/(12μｋ²)}^1/2(ｘ－σｔ),ｋ)と書けます。ここで,σ＝(ｕ₁＋ｕ₂＋ｕ₃)/3＝ｕ₂＋(ｕ₁＋ｕ₃－2ｕ₂)/3＝ｕ₂＋(2－ｋ^-2)Ｕ/3です。

以上から, ｕ(ｘ,ｔ)＝ｕ₂＋Ｕcn²({Ｕ/(12μｋ²)}^1/2[ｘ－{ｕ₂＋(2－ｋ^-2)Ｕ/3}ｔ],ｋ)と書けます。ただし,Ｕ＝ｕ₃－ｕ₂,ｋ²＝(ｕ₃－ｕ₂)/(ｕ₃－ｕ₁)です。

この解は,周期関数cnで表現される定常波形を持つ波列を表わしており,この意味でクノイド波(cnoidal wave)と呼ばれます。

　クノイド波:ｕ(ｘ,ｔ)＝ｕ₂＋Ｕ・cn²({Ｕ/(12μｋ²)}^1/2[ｘ－{ｕ₂＋(2－ｋ^-2)Ｕ/3}ｔ],ｋ)は,速度ｕ₂の一様流の上に振幅がＵの周期波cn²が重なった形をしており,周期波は一様流に相対的に位相速度(2－ｋ^-2)Ｕ/3で進行するという様を表わしています。

　波形cn²(ｓ,ｋ)は,パラメータｋ(0≦ｋ≦1)の値に応じて,cn²(ｓ,0)＝cos²ｓからcn²(ｓ,1)＝sech²ｓまで変化します。

cnの周期は4Ｋですから,これに応じて波長λを{Ｕ/(12μｋ²)}^1/2λ≡4Ｋで定義すると,これはλ＝8Ｋ(3μｋ²/Ｕ)^1/2なる式で与えられます。

　以上では,前提なしで３実根ｕ₁,ｕ₂,ｕ₃ が全て異なる:ｕ₁＜ｕ₂＜ｕ₃と仮定しましたが,特にｕ₂→ｕ₁の極限:ｕ₁＝ｕ₂の重根の場合を想定すると,ｋ²≡(ｕ₃－ｕ₂)/(ｕ₃－ｕ₁)よりｋ＝1です。

そこで,cn(ｓ,ｋ)＝cn(ｓ,1)＝sechｓより,この定常進行波はｕ(ｘ,ｔ)＝ｕ₁＋Ｕsech²[{Ｕ/(12μ)}^1/2{ｘ－(ｕ₁＋Ｕ/3)ｔ}],Ｕ＝ｕ₃－ｕ₁となります。この場合,解は波列でなく,定常波形がsech²の孤立波を表わします。

この孤立波は振幅Ｕの平方根に反比例した拡がりを持ち,一様流ｕ＝ｕ₁に相対的に,位相速度Ｕ/3で進行します。

いま１つの極限ｕ₃→ｕ₂,つまりｕ₂＝ｕ₃の重根の場合はｋ＝0,かつＵ＝0 なので波はｕ＝ｕ₃の一様流に帰着します。　 

ただし,その重根の極限の近傍のｋ～0,Ｕ～0では,Ｕ/ｋ²＝ｕ₃－ｕ₁より,定常解はｕ(ｘ,ｔ)～ｕ₂＋(ｕ₃－ｕ₂)cos²[{Ｕ/(12μｋ²)}^1/2(ｘ－ｕ₁ｔ)]＝ｕ₃＋2ｕ₂sin[{(ｕ₃－ｕ₁)/(3μ)}^1/2(ｘ－ｕ₁ｔ)]となって微小振幅の正弦波となります。

　こうして,"Korteweg-deVries方程式＝Ｋ-ｄＶ方程式"で記述される有限振幅の長波のうちの定常進行波は,クノイド波,または孤立波になることがわかりました。

　しかし,一般に任意の初期条件から出発した非定常の波が,これらの定常解に漸近するとは限りません。

Ｋ-ｄＶ方程式,および類似の非線形発展方程式の研究は近年急速な発展を遂げ,特にＫ-ｄＶ方程式については,かなり多くのことがわかっているようです。

ここでは孤立波に関する２,３の結果を紹介します。

初期条件が三角関数:ｕ(ｘ,0)＝cosπｘのＫ-ｄＶ方程式∂ｕ/∂ｔ＋ｕ(∂ｕ/∂ｘ)＋μ(∂³ｕ/∂ｘ³)＝0 の初期値問題はZabuskyとKruskal(1965)によって数値的に解かれました。

特にμ＝0 なら,数値計算に頼らずとも,解はｕ＝cos[π(ｘ－ｕｔ)]となることがわかります。

※(注)ｕ(ｘ,ｔ)≡cos[π{ｘ－σ(ｘ,ｔ)ｔ}]とおけば,これはｕ(ｘ,0)＝cosπｘを満たします。

そして,∂ｕ/∂ｔ＝π{σ＋(∂σ/∂ｔ)ｔ}sin[π{ｘ－σ(ｘ,ｔ)ｔ}],∂ｕ/∂ｘ＝π{－1＋(∂σ/∂ｘ)ｔ}sin[π{ｘ－σ(ｘ,ｔ)ｔ}]ですから,∂ｕ/∂ｔ＋ｕ(∂ｕ/∂ｘ)＝0 はσ＋(∂σ/∂ｔ)ｔ＋{－1＋(∂σ/∂ｘ)ｔ}cos[π{ｘ－σ(ｘ,ｔ)ｔ}]＝0 となります。

これから,ｔ＝0 ではσ＝cosπｘ＝ｕ(ｘ,0)を得ます。すなわち,σ(ｘ,0)＝ｕ(ｘ,0)です。

そこで,σ＝ｕ＝cos[π(ｘ－σｔ)]とおいてみると,σ＋(∂σ/∂ｔ)ｔ＋{－1＋(∂σ/∂ｘ)ｔ}cos[π{ｘ－σ(ｘ,ｔ)ｔ}]＝σ＋(∂σ/∂ｔ)ｔ＋{－1＋(∂σ/∂ｘ)ｔ}σ＝{∂σ/∂ｔ＋σ(∂σ/∂ｘ)}ｔ＝0 となります。

よって,解の一意性から解はｕ＝cos[π(ｘ－ｕｔ)]です。(注終わり)※

μ＝0 の解ｕ＝cos[π(ｘ－ｕｔ)]では,∂ｕ/∂ｘ＝π{－1＋(∂ｕ/∂ｘ)ｔ}sin[π(ｘ－ｕｔ)]より,[1－πｔsin[π(ｘ－ｕｔ)]](∂ｕ/∂ｘ)＝－πsin[π(ｘ－ｕｔ)]となります。

この解は,ｔ＝ｔ_B＝1/πに点ｘ＝1/2の付近で突っ立ち:|∂ｕ/∂ｘ|＝∞という現象を呈し,その点以後のｘでは解は多価となって物理的意味を持たなくなります。

ZabuskyとKruskal(1965)の数値計算は,μ^1/2＝0.022 において実行され,その結果得られた数値解は,やはりｔ＝ｔ_B付近で突っ立ちに近い状態を示します。

しかし,この場合はゼロでない分散項の効果により,波の重なりは起こらず,逆に連続的な波が,いくつかの孤立波に分散します。

ｔ＝3.6ｔ_Bでは,波形はほぼ完全に分離した孤立波の連なりとなりますが,これらの孤立波は,先に得られた定常波ｕ(ｘ,ｔ)＝ｕ₁＋Ｕsech²[{Ｕ/(12μ)}^1/2{ｘ－(ｕ₁＋Ｕ/3)ｔ}],Ｕ＝ｕ₃－ｕ₁と同じものであることが,その振幅,幅,進行速度の関係から確かめられます。

Ｋ-ｄＶ方程式に従う有限振幅波がある場合に,いくつかの孤立波の集まりになってしまうとすれば,これらの孤立波の間には,どのような相互作用が働くでしょうか？

孤立波の位相速度は振幅に比例するので,ある時刻に大振幅の波を左に,小振幅の波を右にして,互いに十分遠く離しておいたとすれば,大振幅波が小振幅波に近づき,ついには追いつくと予想されます。

Zabusky(1967)は初期条件として,ｕ_j(ｘ,0)＝Ｕ_jsech²(ｘ/Ｄ_j),Ｄ_j≡(12μ/Ｕ_j)^1/2の形を持つ２つの孤立波ｕ_j(ｘ,ｔ)(ｊ＝1,2)を取り,それぞれの振幅をＵ₁＝180,Ｕ₂＝80として,初期に距離を12Ｄ₁(Ｄ₁＝(12μ/Ｕ₁)^1/2＝(μ/15)^1/2)だけ隔てて置きました。

そして,小振幅の方の波が静止するように,さらにｕ₁＝－26＋2/3の一様流を加えました。

数値計算の結果として得られた２つの波は,重なる以前は互いに独立に運動しますが,驚くべきことに衝突以後に再び分離して衝突前と同じ大きさと形と位相速度で運動を続けます。

衝突の影響は,僅かに両者の位相の変化となって残るだけです。

２つの波が初期と同じ間隔12Ｄ₁だけ離れたときの振幅の変化は,僅かにΔＵ₁/Ｕ₁＜0.07％,ΔＵ₂/Ｕ₂＜0.5％なることが見出されました。

　このように,Ｋ-ｄＶ方程式の孤立波解が,あたかも独立の粒子であるかのように挙動するという結果から,この孤立波をソリトン(soliton)と名付けました。

　非線形波動の現象は,水の波に限らずプラズマ振動,さらに素粒子のモデルなどと関連して研究者の興味を集めており,非線形現象の解明と数学的理論の開発の両面にわたって,その発展が期待されます。

　水の波の話から自然にソリトンを紹介するという所期の目的が達せられたので,水の波の記事シリーズをこれで終わります。

参考文献:巽友正著「流体力学」(培風館)

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

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

← クリックして投票してください。(ブログ村科学ブログランキング投票です。１クリック＝１投票です。１人１日１投票しかできません。クリックするとブログ村の人気ランキング一覧のホ－ムページームに跳びます。）

http://www.mediator.co.jp/category/pages.php?id=115「中古パソコン！メディエーター巣鴨店」

Dell-個人のお客様ページ

(Dellの100円パソコン(Mini9)↓私も注文しました。)

 

ヤーマン　プラチナゲルマローラー　1日3分コロコロエステ！ローラー型プラチナ配合美顔器  

　

お売りください。ブックオフオンラインのインターネット買取 展開へ！ ▼コミック　尾田栄一郎 「ONE PIECE（52）」 ▼コミック 「ONE PIECE」をオトナ買い

三国志特集　▼コミック 横山光輝 「三国志全巻セット」 「三国志（文庫版）全巻セット」   「三国志（ワイド版）全巻セット」 　▼書籍 「三国志」/吉川英治   「三国志」/北方謙三   「三国志」/宮城谷昌光

水の波(2)(浅水波,深水波,表面張力波)

　水の波のつづきです。微小振幅波を近似して分類します。

微小振幅波の位相速度ｃと波数ｋ or 波長λとの関係,つまり分散関係はｃ＝ω/ｋ＝{(ｇ/ｋ)tanh(ｋｈ)}^1/2＝{ｇλtanh(2πｈ/λ)/(2π)}^1/2で与えられます。

これは,一見しただけではよくわからない関数形をしていますが,波長λが水深ｈに対して極めて大きい波,または水深ｈの方がλに対して極めて小さい波に対しては分散関係は簡単になります。

すなわち,ｈ/λ＜＜1では,tanh(ｋｈ)＝tanh(2πｈ/λ)～ ｋｈ＝2πｈ/λなる近似が成り立つので,位相速度ｃもｃ～ (ｇｈ)^1/2と近似できます。この場合,波の速さｃは深さの平方根:ｈ^1/2に比例します。

ｃが波長λに依らないので波は非分散的です。

ｈ/λ＜＜1はｋｈ＜＜1を意味しますから,今の場合には速度ポテンシャルΦも,Φ(ｘ,ｚ,ｔ)＝－{Ａｃ/sinh(ｋｈ)}cosh{ｋ(ｚ＋ｈ)}cos(ｋｘ－ωｔ)の正弦波の係数が,sinh(ｋｈ)～ ｋｈ,cosh{ｋ(ｚ＋ｈ)～ 1と近似できて,Φ(ｘ,ｚ,ｔ)～ {－Ａｃ/(ｋｈ)}cos(ｋｘ－ωｔ)と簡単になります。

これによって,ｄｘ/ｄｔ＝ｕ＝∂Φ/∂ｘ＝(Ａｃ/ｈ)sin(ｋｘ－ωｔ),ｄｚ/ｄｔ＝ｕ＝∂Φ/∂ｚ＝{－Ａｃｋ(ｚ＋ｈ)/ｈ}cos(ｋｘ－ωｔ)が得られます。

そこで前に求めた近似軌道はｘ＝ｘ₀＋[Ａ/(ｋｈ)]cos(ｋｘ₀－ωｔ),ｚ＝ｚ₀＋[Ａ(ｚ₀＋ｈ)/ｈ]sin(ｋｘ₀－ωｔ)と簡単になります。

これは,楕円軌道:(ｘ－ｘ₀)²/ａ²＋(ｚ－ｚ₀)²/ｂ²＝1における長半径ａがａ＝Ａ/(ｋｈ),短半径ｂがｂ＝Ａ(ｚ₀＋ｈ)/ｈとなってａ＞＞ｂを意味します。

そこで,ｚ方向の振幅ｂ＝Ａ(ｚ₀＋ｈ)/ｈは無視できて,水の運動は事実上,水平方向の単振動ｘ＝ｘ₀＋ａcos(ｋｘ₀－ωｔ),ｚ＝ｚ₀と見なせます。 

そして振幅ａ＝Ａ/(ｋｈ)は深さｚ＝ｚ₀に依存しませんから,同じｘ₀の位置にある水の鉛直層は一体になって水平に往復運動をするという描像になります。

その水平運動の速度はｕ＝ｄｘ/ｄｔ＝Ａ(ｇ/ｈ)^1/2sin(ｋｘ－ωｔ)となりますが,これは水面の波形(高さ):η(ｘ,ｔ)＝Ａsin(ｋｘ－ωｔ)と全く同位相なので水面の山では正の向き,谷では負の向きに動くことになります。

このような波は,波長λが長いか,深さｈが小さい場合の波なので,長波,または浅水波(shallow water wave)と呼ばれます。

そして,これを導くために用いた上述の近似を,長波近似,または浅水波近似といいます。

一方,逆に波長λが深さｈに対して極めて短かい波,あるいはｈがλに比べて極めて大きい場合,深い場合の近似を考えることもできます。

この場合はｋｈ＝2πｈ/λ＞＞1なのでtanh(ｋｈ)～ 1と近似できます。そこで,位相速度ｃ＝ω/ｋ＝{(ｇ/ｋ)tanh(ｋｈ)}^1/2＝{ｇλtanh(2πｈ/λ)/(2π)}^1/2の近似はｃ～ (ｇ/ｋ)^1/2＝{ｇλ/(2π)}^1/2となります。

そこで,これの極限での波は分散的です。

このｋｈ＞＞1の場合,sinh(ｋｈ)＝(1/2){exp(ｋｈ)－exp(－ｋｈ)}～ (1/2)exp(ｋｈ),cosh{ｋ(ｚ＋ｈ)}＝(1/2)[exp{ｋ(ｚ＋ｈ)}＋exp{－ｋ(ｚ＋ｈ)}] ～ (1/2)exp{ｋ(ｚ＋ｈ)}と近似できます。

そこで,Φ(ｘ,ｚ,ｔ)＝－{Ａｃ/sinh(ｋｈ)}cosh{ｋ(ｚ＋ｈ)}cos(ｋｘ－ωｔ)も,近似的にΦ(ｘ,ｚ,ｔ)～－Ａｃexp(ｋｚ)cos(ｋｘ－ωｔ)と書けます。

したがって,ｄｘ/ｄｔ＝ｕ＝∂Φ/∂ｘ～ Ａｃｋexp(ｋｚ)sin(ｋｘ－ωｔ),ｄｚ/ｄｔ＝ｗ＝∂Φ/∂ｚ～－Ａｃｋexp(ｋｚ)cos(ｋｘ－ωｔ)ですから,水の運動もｘ～ｘ₀＋Ａexp(ｋｚ₀)cos(ｋｘ₀－ωｔ),ｚ～ ｚ₀＋Ａexp(ｋｚ₀)sin(ｋｘ₀－ωｔ)となります。

前に,ｘ＝ｘ₀＋[Ａcosh{ｋ(ｚ₀＋ｈ)}/sinh(ｋｈ)]cos(ｋｘ₀－ωｔ),ｚ＝ｚ₀＋[Ａsinh{ｋ(ｚ₀＋ｈ)}/sinh(ｋｈ)]sin(ｋｘ₀－ωｔ)によって,楕円軌道:(ｘ－ｘ₀)²/ａ²＋(ｚ－ｚ₀)²/ｂ²＝1,ａ≡Ａcosh{ｋ(ｚ₀＋ｈ)}/sinh(ｋｈ),ｂ≡Ａsinh{ｋ(ｚ₀＋ｈ)}/sinh(ｋｈ)を求めました。

これは,今の近似では,ａ＝ｂ＝Ａexp(ｋｚ₀)なので,水の運動は(ｘ－ｘ₀)²＋(ｚ－ｚ₀)²＝Ａ²exp(2ｋｚ₀)となり,(ｘ₀,ｚ₀)を中心とする半径がＡexp(ｋｚ₀)の円運動になります。

半径Ａexp(ｋｚ₀)は深さ|ｚ₀|の増加に連れて指数関数的に減少するため,水面の波による水の運動は実質的にある有限深さまでに限られ,それより下方には及びません。

その有限深さ|ｚ₀|の値はｋ|ｚ₀|＝2π|ｚ₀|/λのオーダーで決まるため,波の波長λに比例しますから,波長が短くなると共に小さく浅くなることがわかります。

このように,水深に比べて波長が極めて短かい波を短波,または深水波(deep-water wave)と呼び,これを導くために用いた近似を短波近似,または深水波近似といいます。これが成立するのは,tanh(ｋｈ)～ 1の近似が許される場合です。

例えば,1－tanh(ｋｈ)≦0.01となるのはｋｈ＝2πｈ/λ≧2.65のときですから,λ≦2.4ｈ,すなわち水深の約２倍の波長に対して短波近似は相対誤差１％以内で成立します。

さて,これまでの扱いでは,水面に働く表面張力を考慮しませんでしたが表面張力の影響は特に波長の短かい波に対しては無視できません。

表面張力を考慮する場合でも,元々の境界値問題:ｕ＝∇Φ,∇²Φ＝0,∂Φ/∂ｚ＝∂η/∂ｔ＋(∂Φ/∂ｘ)(∂η/∂ｘ)＋(∂Φ/∂ｙ)(∂η/∂ｙ),∂Φ/∂ｔ＋ｐ/ρ＋(∇Φ)²/2＋ｇｚ＝ｆ(ｔ)は,そのまま成立します。

そして,大気と接している水面上の大気圧をｐ₀とするとき,前には∂Φ/∂ｔ＋ｐ/ρ＋(∇Φ)²/2＋ｇｚ＝ｆ(ｔ)なる圧力方程式で,任意関数ｆ(ｔ)をｆ(ｔ)≡ｐ₀/ρと置き,ｚ＝ηでｐ＝ｐ₀として,η＝－ｇ^-1(∂Φ/∂ｔ)－(2ｇ)^-1(∇Φ)²なる式を得ました。

しかし,表面張力による圧力差δｐ＞0 があれば水面での水の圧力ｐは大気圧ｐ₀と同じではなくて,ｐ₀＋δｐに等しくなります。 

そこで,この場合は同じｆ(ｔ)＝ｐ₀/ρに対しｚ＝ηでの圧力としてｐ＝ｐ₀＋δｐを代入すれば,水面での境界条件はη＝－ｇ^-1(∂Φ/∂ｔ)－(2ｇ)^-1(∇Φ)²－δｐ/(ρｇ)となります。

ところで,１つの曲面をその法線を含む平面で切ると,切り口の曲線の曲率半径Ｒは一般に平面の方向によって異なります。

このときのＲの極大値をＲ₁,極小値をＲ₂とするとき,それぞれを与える平面の方向を主方向,κ₁≡1/Ｒ₁,κ₂≡1/Ｒ₂を主曲率,Ｈ≡(κ₁＋κ₂)/2を平均曲率といいます。

γを表面張力とすると,これが働いている曲面を隔てての圧力差は,δｐ＝γ[(1/Ｒ₁)＋(1/Ｒ₂)]＝γ(κ₁＋κ₂)＝2γＨで与えられます。

ところで,水面に限らず一般的な空間曲面の性質について2008年７/30の過去記事「相対論の幾何学(第Ⅰ部-4,空間曲面(1))」http://maldoror-ducasse.cocolog-nifty.com/blog/2008/07/41_5587.html,および2008年８/３の記事「相対論の幾何学(第Ⅰ部-5,空間曲面(2))」http://maldoror-ducasse.cocolog-nifty.com/blog/2008/08/52_65b1.htmlに詳しく論じています。

この中では,滑らかな曲面をある定義域Ｄにある２つのパラメータ(ｕ,ｖ)∈Ｄのベクトル値関数としてｒ(ｕ,ｖ)＝(ｘ(ｕ,ｖ),ｙ(ｕ,ｖ),ｚ(ｕ,ｖ))で定義しています。

そして,この曲面の曲率に関して得られる平均曲率Ｈ≡(κ₁＋κ₂)/2 はｕ,ｖの６つの関数Ｅ,Ｆ,Ｇ,Ｌ,Ｍ,Ｎによって,Ｈ＝(ＥＮ＋ＧＬ－2ＦＭ)/{2(ＥＧ－Ｆ²)}なる式で与えられることが示されています。

ただし,Ｅ,Ｆ,Ｇは,ｒ_u≡∂ｒ/∂ｕ,ｒ_v≡∂ｒ/∂ｖに対しＥ≡ｒ_u²,Ｆ≡ｒ_uｒ_v＝ｒ_vｒ_u,Ｇ≡＝ｒ_v²で定義されます。

また,Ｌ,Ｍ,Ｎはｒ_uu≡∂²ｒ/∂ｕ²,ｒ_uv≡∂²ｒ/∂ｕ∂ｖ,ｒ_vv≡∂²ｒ/∂ｕ²に対し,Ｌ≡(ｒ_uu,ｅ),Ｍ≡(ｒ_uv,ｅ),Ｎ≡(ｒ_vv,ｅ)で定義されます。

ただし,ｅ≡(ｒ_u×ｒ_v)/|ｒ_u×ｒ_v|は曲面ｒ(ｕ,ｖ)の法線単位ベクトルです。 

これらの関数を今の(ｘ,ｙ)をパラメータとして水面の高さを表わす式ｚ＝η(ｘ,ｙ),または,これの曲面ベクトルの表現ｒ(ｘ,ｙ)＝(ｘ,ｙ,η(ｘ,ｙ))で表わすと,Ｅ＝1＋(∂η/∂ｘ)²,Ｆ≡(∂η/∂ｘ)(∂η/∂ｙ),Ｇ＝1＋(∂η/∂ｙ)²です。

それ故,ＥＧ－Ｆ²＝1＋(∂η/∂ｘ)²＋(∂η/∂ｙ)²となります。 

そして,Ｌ＝－∂²η/∂ｘ²,Ｍ＝－∂²η/∂ｘ∂ｙ,Ｎ＝－∂²η/∂ｙ²ですから,ＥＮ＋ＧＬ－2ＦＭ＝－(∂²η/∂ｘ²){1＋(∂η/∂ｘ)²}－(∂²η/∂ｙ²){1＋(∂η/∂ｙ)²}＋2(∂²η/∂ｘ∂ｙ)(∂η/∂ｘ)(∂η/∂ｙ)が得られます。

そこで,微小振幅波近似を採用してηの２次以上の項を無視すると,平均曲率はＨ～(－1/2)(∂²η/∂ｘ²＋∂²η/∂ｙ²)/2となります。

したがって,表面張力による圧力差はδｐ＝γ[(1/Ｒ₁)＋(1/Ｒ₂)]＝2γＨ～－γ(∂²η/∂ｘ²＋∂²η/∂ｙ²)と近似表現されます。

そこで,δｐを考慮した方程式:η＝－ｇ^-1(∂Φ/∂ｔ)－(2ｇ)^-1(∇Φ)²－δｐ/(ρｇ)においてηとΦの２次以上の項を無視した微小振幅波では,η＝－ｇ^-1(∂Φ/∂ｔ)＋{γ/(ρｇ)}(∂²η/∂ｘ²＋∂²η/∂ｙ²)が得られます。

これと,水面ｚ＝ηにおける運動学的境界条件∂Φ/∂ｚ＝∂η/∂ｔ＋(∂Φ/∂ｘ)(∂η/∂ｘ)＋(∂Φ/∂ｙ)(∂η/∂ｙ)の微小振幅波近似∂Φ/∂ｚ＝∂η/∂ｔ(ｚ～0)から,ｚ＝0で成立すべき式:∂²Φ/∂ｔ²＋[ｇ－(γ/ρ)(∂²/∂ｘ²＋∂²/∂ｙ²)](∂Φ/∂ｚ)＝0 を得ます。

これは,前の記事でｚ＝0で∂²Φ/∂ｔ²＋ｇ(∂Φ/∂ｚ)＝0,とｚ＝－ｈで∂Φ/∂ｚ＝0 なる境界条件の下でラプラス方程式∇²Φ＝0 を解いて近似解を得た際のｚ＝0で∂²Φ/∂ｔ²＋ｇ(∂Φ/∂ｚ)＝0 という境界条件に代わる式です。

ラプラス方程式∇²Φ＝∂²Φ/∂ｘ²＋∂²Φ/∂ｚ²＝0 の変数分離解がΦ(ｘ,ｚ,ｔ)＝Ｃcosh{ｋ(ｚ＋ｈ)}cos(ｋｘ－ωｔ),η(ｘ,ｔ)＝Ａsin(ｋｘ－ωｔ) (ただし,Ａ≡－Ｃωcosh(ｋｈ))で与えられるのは前と同じです。

違うのはωやｃをｋ,またはλと関係付ける分散式です。

これは解Φの陽な式を,"ｚ＝0 で∂²Φ/∂ｔ²＋[ｇ－(γ/ρ)(∂²/∂ｘ²)](∂Φ/∂ｚ)＝0 が成り立つ"という境界条件式に代入すれば得られて－ω²cosh(ｋｈ)＋ｋ(ｇ＋γｋ²/ρ)sinh(ｋｈ)＝0 となります。

これから,ω＝{ｋ(ｇ＋γｋ²/ρ)tanh(ｋｈ)}^1/2,あるいはｃ＝ω/ｋ＝{(ｇ/ｋ＋γｋ/ρ)tanh(ｋｈ)}^1/2＝[{ｇλ/(2π)＋2πγ/(ρλ)}tanh(2πｈ/λ)]^1/2を得ます。 

しかし,長波,または浅水波;ｋｈ＝2πλ＜＜1の場合はtanh(ｋｈ)～ｋｈ,(γｋ/ρ)tanh(ｋｈ)～γｋ²ｈ/ρ＜＜ｇｈより,ｃ～(ｇｈ)^1/2と近似されて,表面張力を考慮しない場合と全く同じです。

つまり,長波,または浅水波には表面張力の影響は現われません。

これに対して,短波,または深水波;ｋｈ＝2πλ＞＞1の場合には,tanh(ｋｈ)～1によって,ｃ～(ｇ/ｋ＋γｋ/ρ)^1/2＝{ｇλ/(2π)＋2πγ/(ρλ)}^1/2となります。表面張力γは速度ｃにもろに影響します。

ここで,ｃ²＝ｇλ/(2π)＋2πγ/(ρλ)をλで微分してゼロと置くと,ｄｃ²/ｄλ＝ｇ/(2π)－2πγ/(ρλ²)＝0 です。

そこで,λ_m≡2π{γ/(ρｇ)}^1/2と置けば,位相速度ｃは波長λ＝λ_mにおいて,最小値ｃ_m≡(4γｇ/ρ)^1/4をとることになります。

短波の中でも比較的波長が短かくλ＜λ_mの場合には,ｃ～(ｇ/ｋ＋γｋ/ρ)^1/2＝{ｇλ/(2π)＋2πγ/(ρλ)}^1/2の右辺の括弧の中の第２項の表面張力項が優越します。

このような波を表面張力波,あるいは"さざ波"といいます。表面張力波では波長λが小さいほど位相速度ｃは大きく速い波で,λが増すにつれて遅い波になります。

一方,比較的波長が長くてλ＞λ_mの場合には,第１項ｇ/ｋ＝ｇλ/(2π)が優越します。このような波を重力波といいます。この場合には表面張力がない場合と同じく,波長λの増減と共にｃも増減します。

今日はここまでにします 

参考文献:巽友正著「流体力学」(培風館),小林昭七 著「曲線と曲面の微分幾何」(裳華房)

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

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

← クリックして投票してください。(ブログ村科学ブログランキング投票です。１クリック＝１投票です。１人１日１投票しかできません。クリックするとブログ村の人気ランキング一覧のホ－ムページームに跳びます。）

http://www.mediator.co.jp/category/pages.php?id=115「中古パソコン！メディエーター巣鴨店」

Dell-個人のお客様ページ

(Dellの100円パソコン(Mini9)↓私も注文しました。)

 

ヤーマン　プラチナゲルマローラー　1日3分コロコロエステ！ローラー型プラチナ配合美顔器  

　

お売りください。ブックオフオンラインのインターネット買取 展開へ！ ▼コミック　尾田栄一郎 「ONE PIECE（52）」 ▼コミック 「ONE PIECE」をオトナ買い

三国志特集　▼コミック 横山光輝 「三国志全巻セット」 「三国志（文庫版）全巻セット」   「三国志（ワイド版）全巻セット」 　▼書籍 「三国志」/吉川英治   「三国志」/北方謙三   「三国志」/宮城谷昌光

揺動散逸定理

　非平衡統計熱力学の１過程としてブラウン運動などに関わる揺動散逸定理(fluctuation dissipation theorem)について述べてみます。

　この定理は誰が起源なのかよく知らないのですが,日本では線型応答理論の一環として統計物理学の重鎮であった久保亮五先生などが関わっていたと記憶しています。

　一般に,ある物理量(示量変数の組):ａ＝(ａ₁,ａ₂,..,ａ_n)があって,エントロピーＳがａの関数であるとき,系の時間発展はａに対する１階微分方程式で表現されて,それはｄａ_i/ｄｔ＝(ｄａ_i/ｄｔ)_rev＋(ｄａ_i/ｄｔ)_irrのように,可逆な部分(ｄａ_i/ｄｔ)_revと不可逆な部分(ｄａ_i/ｄｔ)_irrの和として与えられます。

　そして,可逆部分はさらに,(ｄａ_i/ｄｔ)_rev＝Σ_j{ａ_i,ａ_j}(∂Ｓ/∂ａ_j)という構造を持つとします。

　ここで,{Ｘ,Ｙ}はポアソン括弧のように,{Ｙ,Ｘ}＝－{Ｘ,Ｙ}という反対称性を持ち,{Ｘ,ｆ}＝Σ_j{Ｘ,ａ_j}(∂ｆ/∂ａ_j)という性質で規定される量であるとします。

こう定義すると,(ｄＳ/ｄｔ)_rev＝Σ_i(∂Ｓ/∂ａ_i)(ｄａ_i/ｄｔ)_rev＝Σ_iΣ_j(∂Ｓ/∂ａ_i){ａ_i,ａ_j}(∂Ｓ/∂ａ_j)＝{Ｓ,Ｓ}＝0 となって,可逆過程ではエントロピーは生成されないことになります。

実際には,断熱可逆変化か,あるいは可逆であって,かつサイクルである場合以外なら,可逆過程でもエントロピーの変化はありますから,これはどう解釈すべきなんでしょうか？

ｄａ_i/ｄｔ＝(ｄａ_i/ｄｔ)_rev＋(ｄａ_i/ｄｔ)_irrは"可逆な部分と不可逆な部分の和である。"というよりもむしろ,"エントロピー非生成部分とエントロピー生成部分の和である。"と事実のままを述べた方がいいのかもしれません。

一方,不可逆部分については現象論的に(ｄａ_i/ｄｔ)_irr＝Σ_jＬ_ij(∂Ｓ/∂ａ_j)と表わすことにします。

というのは,ａ_jが示量変数のとき∂Ｓ/∂ａ_jは示強パラメータであり,平衡の近傍では不可逆部分は示強パラメータ,あるいはその空間勾配に比例して進行するからです。

実際,問題としている系を局所平衡状態にある部分系の集まりと考えると,それぞれの部分系の物質密度をρ,単位質量当たりのエントロピーをｓとしたとき,エントロピー密度(ρｓ)の変化は,一般に熱力学の関係式により示強パラメータＦ_iと示量変数ａ_iによって,ｄ(ρｓ)＝Σ_iＦ_iｄａ_iと表わされます。

この表式では確かに示強パラメータは,Ｆ_i＝(∂Ｓ/∂ａ_j)/Ｖと表わされています。

そして今,対象としている系が隣り合う２つの部分系Ａ,Ｂだけから成るとし,Ａ,Ｂが等しい体積Ｖを持つとするとき,各部分系のエントロピーＳ＝Ｖρｓは,Ｘ_i＝Ｖａ_iの関数であると考えられます。

今,Ａ,Ｂのエントロピーを,それぞれＳ^A,Ｓ^Bとし,Ｘ_iがＡからＢにΔＸ_i＝ＶΔａ_iだけ移動するとします。

Ｘ_iが系全体では保存する量であって,初めＡにはＸ_iがＸ_i^AだけＢにはＸ_i^Bだけあったとすると,ＡからＢへのΔＸ_i＝ＶΔａ_iの移動による系全体のエントロピーＳの増分は,ΔＳ＝Ｓ^A(Ｘ_i^A－ΔＸ_i)＋Ｓ^B(Ｘ_i^B＋ΔＸ_i)－[Ｓ^A(Ｘ_i^A)＋Ｓ^B(Ｘ_i^B)]～Σ_i(－∂Ｓ/∂Ｘ_i^A＋∂Ｓ/∂Ｘ_i^B)ΔＸ_i＝Σ_i(－Ｆ^A_i＋Ｆ^B_i)ΔＸ_iとなります。

ここでＦ^A_i,Ｆ^Bはそれぞれ部分系Ａ,ＢにおけるＦ_iの値を表わしています。

系全体を孤立系と考えると熱力学第二法則によって,ΔＳ＞0 でなければならないので,１つの示量変数Ｘ_i＝Ｖａ_iのみに着目してＡからＢへと微小量ΔＸ_i＞0 の移動が起こるためには,示強パラメータＦ_iについてはＦ^B_i＞Ｆ^A_iであることが必要になります。

このことから,示量変数Ｘ_i＝Ｖａ_iが保存量のときはＸ_iの輸送を引き起こす駆動力となるのは示強パラメータＦ_iの空間勾配であると考えられます。

また,示量変数Ｘ_i＝Ｖａ_iが非保存量のときはΔＳ＝Σ_iＦ_iΔＸ_iにおいてＦ_i＝(∂Ｓ/∂ａ_i)/Ｖ＞0 ならば,ΔＸ_i＝ＶΔａ_i；Δａ_i＞0 なる変化が不可逆過程として進行し得ます。

つまり,一般にｄＳ＝Σ_iＦ_iｄａ_i,Ｆ_i＝(∂Ｓ/∂ａ_i)と書けますが,量ａ_kが保存量のとき,その保存方程式は∂ａ_k/∂ｔ＋∇Ｊ_k＝0 であり,その流れＪ_kは一般にＪ_k＝Σ_jＬ_kj∇Ｆ_jと,示強的な量Ｆ_j＝(∂Ｓ/∂ａ_j)の勾配を駆動力とする形に表わされます。

そこで,ａ_k'≡ａ_kｒとおくと,Ｊ_k＝ａ_kｖ＝ａ_k(ｄｒ/ｄｔ)ですから,ｄａ_k'/ｄｔ＝ｄ(ａ_kｒ)/ｄｔ～Ｊ_k＝Σ_jＬ_kj∇Ｆ_j＝Σ_jＬ_kj∇(∂Ｓ/∂ａ_j)です。

そこで,もしもＳがａ_jの１次関数なら,∇(∂Ｓ/∂ａ_j)～∂Ｓ/∂ａ_j'となり,示量変数ａ_k'の求める時間発展の形式ｄａ_k'/ｄｔ＝Σ_jＬ_kj(∂Ｓ/∂ａ_j')が得られます。しかし正直なところかなり苦しいです。

こうして,はっきりと証明されたわけではないのですが,現象論的発展方程式はｄａ_i/ｄｔ＝{ａ_i,Ｓ}＋Σ_jＬ_ij(∂Ｓ/∂ａ_j)と表わされるとします。

これは非平衡な初期条件から平衡状態への緩和を表わしています。

そこでａ(ｔ)の時間発展は初期条件ａ(0)＝ａに対してａ_i(ｔ)＝ａ_i＋ｔ[{ａ_i,Ｓ}＋Σ_jＬ_ij(∂Ｓ/∂ａ_j)]＋..で与えられます。

平衡状態においても,一般に巨視的な物理量ａ＝(ａ₁,ａ₂,..,ａ_n)はゆらいでいます。たまたま,ある時刻にａ(ｔ)がａ(ｔ₀)＝ａという値をとったとして,その後の時間発展を観測します。

ａ(ｔ₀＋ｔ)は初期時刻ｔ₀によってさまざまな値を取りますが,それらの平均を取ったものも現象論的発展と同じになるとします。

すなわち,初期条件ａ(ｔ₀)＝ａが与えられると短い時間では平均的に＜ａ_i(ｔ₀＋ｔ)＞_{ａ(ｔ0)=ａ}＝ａ_i＋ｔ[{ａ_i,Ｓ}＋Σ_jＬ_ij(∂Ｓ/∂ａ_j)]となると仮定します。これを線型減衰の仮定と呼びます。

平衡状態におけるゆらぎの時間相関関数＜ａ_i(ｔ₀＋ｔ)ａ_k(ｔ₀)＞を求めるには＜ａ_i(ｔ₀＋ｔ)＞_{ａ(ｔ0)＝ａ}に初期値ａ_k(ｔ₀)＝ａ_kを掛けてａ_kの分布について平均すればいいので,ｔ≧0 に対して＜ａ_i(ｔ₀＋ｔ)ａ_k(ｔ₀)＞_eq＝＜＜ａ_i(ｔ₀＋ｔ)＞_{ａ(ｔ0)=ａ}ａ_k＞_eq ＝＜ａ_iａ_k＞_eq＋ｔ[＜{ａ_i,Ｓ}ａ_k＞_eq＋Σ_jＬ_ij＜(∂Ｓ/∂ａ_j)ａ_k＞_eq]となります。

平衡状態におけるゆらぎに対するボルツマン・アインシュタインの原理,つまり,ボルツマンの原理Ｓ＝k_BlnＷから,逆に微視的状態数ＷがＷ(ａ)＝exp[Ｓ(ａ)/ｋ_B]と書ける,ことを用います。

全状態数をＷとしたときにＷ(ａ)/Ｗが変数ａの状態が実現する確率となりますから,平衡状態での関数ｆ(ａ)の平均値は＜ｆ(ａ)＞_eq＝(1/Ｗ)∫ｄａ₁..ｄａ_nｆ(ａ)exp[Ｓ(ａ)/ｋ_B]で与えられます。

そこで,＜(∂Ｓ/∂ａ_j)ａ_k＞_eq＝∫ｄａ₁..ｄａ_n(∂Ｓ/∂ａ_j)ａ_kexp[Ｓ(ａ)/ｋ_B]＝－ｋ_Bδ_kjとなります。

したがって,＜ａ_i(ｔ₀＋ｔ)ａ_k(ｔ₀)＞_eq＝＜ａ_iａ_k＞_eq＋ｔ[＜{ａ_i,Ｓ}ａ_k＞_eq－ｋ_BＬ_ik]となりますが,右辺の{ａ_i,Ｓ}はｄａ_i/ｄｔの可逆部分です。

ところが,平衡状態では物理量の任意の関数は可逆変化では変化しないので,＜{ｆ(ａ),Ｓ}＞_eq＝0 です。もっともこれはその現象論の範囲で厳密に証明できるわけではないので仮説として導入するわけです。

そして,この式でｆ(ａ)＝ａ_iａ_kを代入すると,＜ａ_j{ａ_iδ_jk,Ｓ}＋ａ_j{ａ_iδ_ijａ_k,Ｓ}＞_eq＝0 :すなわち＜{ａ_i,Ｓ}ａ_k＞_eq＝－＜ａ_i{ａ_k,Ｓ}＞_eqが得られます。

そこで、物理変数ａ_i,ａ_kの時間反転対称性に関して次の２つの場合を考えます。：

(Ⅰ)変数ａ_i,ａ_kが共に時間反転に対して対称である場合

このときは＜ａ_i(ｔ₀＋ｔ)ａ_k(ｔ₀)＞_eq＝＜ａ_i(ｔ₀－ｔ)ａ_k(ｔ₀)＞_eqです。さらに平衡状態の定常性から物理量の時間相関は時間差のみに依存しｔ₀には依存しないので＜ａ_i(ｔ₀－ｔ)ａ_k(ｔ₀)＞_eq＝＜ａ_i(ｔ₀)ａ_k(ｔ₀＋ｔ)＞_eqです。

以上から,＜ａ_i(ｔ₀＋ｔ)ａ_k(ｔ₀)＞_eq＝＜ａ_i(ｔ₀)ａ_k(ｔ₀＋ｔ)＞_eqが得られます。

そこで両辺に線型減衰の仮定を適用すると,＜ａ_iａ_k＞_eq＋ｔ[＜{ａ_i,Ｓ}ａ_k＞_eq－ｋ_BＬ_ik]＝＜ａ_iａ_k＞_eq＋[＜{ａ_k,Ｓ}ａ_i＞_eq－ｋ_BＬ_ki]となるはずです。

これが実際に成り立つためには,＜{ａ_i,Ｓ}ａ_k＞_eq＝＜{ａ_i,Ｓ}ａ_k＞_eqかつＬ_ik＝Ｌ_kiが満たされなければなりません。

前者は＜{ａ_i,Ｓ}ａ_k＞_eq＝－＜ａ_i{ａ_k,Ｓ}＞_eqと組み合わせると＜{ａ_i,Ｓ}ａ_k＞_eq＝＜{ａ_i,Ｓ}ａ_k＞_eq＝0 となります。要するに,ａ_i,ａ_kが共に時間反転に対して対称ならば,その時間微分{ａ_k,Ｓ},{ａ_i,Ｓ}は明らかに時間反転に対して反対称となりますから,時間反転対称なａ_iやａ_kとの積は平衡状態ではゼロとなることを表現しています。

(Ⅱ) 時間反転に対して変数ａ_iが反対称でａ_kが対称の場合

このときは＜ａ_i(ｔ₀＋ｔ)ａ_k(ｔ₀)＞_eq＝－＜ａ_i(ｔ₀－ｔ)ａ_k(ｔ₀)＞_eqです。そこで(Ⅰ)と同様に考えて＜ａ_i(ｔ₀＋ｔ)ａ_k(ｔ₀)＞_eq＝－＜ａ_i(ｔ₀)ａ_k(ｔ₀＋ｔ)＞_eqとなります。

したがって線型減衰の式から＜ａ_iａ_k＞_eq＋ｔ[＜{ａ_i,Ｓ}ａ_k＞_eq－ｋ_BＬ_ik]＝－＜ａ_iａ_k＞_eq－ｔ[＜{ａ_k,Ｓ}ａ_i＞_eq－ｋ_BＬ_ki](ｔ≧0)ですが,これが成り立つためには＜ａ_iａ_k＞＝0 でかつＬ_ik＝－Ｌ_kiが満たされなければなりません。

かくして,変数ａ_i,ａ_kが同じ時間対称性を持つときにはＬ_ik＝Ｌ_ki, 反対の時間対称性を持つときにはＬ_ik＝－Ｌ_kiとなります。

さらに係数Ｌ_ikが時間反転に対して反対称な外部パラメータＡ(Ａは例えば磁場Ｂや速度ｖ)に依存するときには,それぞれＬ_ik(Ａ)＝Ｌ_ki(－Ａ),Ｌ_ik(Ａ)＝－Ｌ_ki(－Ａ)となります。

発展方程式の不可逆部分を熱力学的力∂Ｓ/∂ａで表現する輸送係数Ｌ_ijについての上述の対称性をオンサーガーの相反定理と呼び,この係数Ｌ_ikをオンサーガー係数と呼びます。

さらに発展方程式の不可逆部分はエントロピー生成をするという熱力学第二法則の要請から係数行列(Ｌ_ij)は正値行列です。

なぜなら,ｄＳ/ｄｔ＝(ｄＳ/ｄｔ)_irr＝Σ_i(∂Ｓ/∂ａ_i)(ｄａ_i/ｄｔ)_irr＝Σ_iΣ_j(∂Ｓ/∂ａ_i)Ｌ_ij(∂Ｓ/∂ａ_j)＞0 となるべきことが要求されますが,(∂Ｓ/∂ａ_i)はベクトルとして任意の値を取ると考えてよいからです。

平衡状態での巨視的変数の現象論的な発展方程式がｄａ_i/ｄｔ＝{ａ_i,Ｓ}＋Σ_jＬ_ij(∂Ｓ/∂ａ_j)で与えられるので,必然的に存在する"ゆらぎ＝揺動あるいは雑音"Ｒ_i(ｔ)の存在を考慮すると,ａの時間発展は一般的な確率微分方程式ｄａ_i/ｄｔ＝{ａ_i,Ｓ}＋Σ_jＬ_ij(∂Ｓ/∂ａ_j)＋Ｒ_i(ｔ)で表わされると考えられます。

ここで通常はゆらぎＲ_i(ｔ)は完全にランダムであり＜Ｒ_i(ｔ)＞＝0 ,＜Ｒ_i(ｔ)Ｒ_j(ｔ')＞_eq＝2Ｄ_ijδ(ｔ－ｔ')なる白色雑音(white noise)で与えられます。

こうすると,時刻ｔにおいて状態ａが実現する確率Ｐ(ａ,ｔ)に対して,次のフォッカー・プランク方程式(Fokker-Planck)が得られます。

∂Ｐ(ａ,ｔ)/∂ｔ＝－Σ_i(∂/∂ａ_i){ａ_i,Ｓ}Ｐ(ａ,ｔ)＋Σ_iΣ_j(∂/∂ａ_i)[－Ｌ_ij(∂Ｓ/∂ａ_j)＋Ｄ_ij(∂/∂ａ_j)]Ｐ(ａ,ｔ)です。

一般に,１次元で考えたとき外力Ｆがないときのゆらぎも含めた粒子の運動は,その速度をｕとするとき,次の"運動方程式＝ランジュバン方程式(Langevin)"ｄｕ/ｄｔ＝－γｕ＋Ｒ(ｔ)/ｍ に従います。

そして,ここでもゆらぎＲ(ｔ)は白色雑音,つまり＜Ｒ(ｔ)＞＝0 ,＜Ｒ(ｔ)Ｒ(ｔ')＞＝2Ｄ_uδ(ｔ－ｔ')を満たしているとします。

このとき,時刻ｔ₁に速度ｕ₁を持っていた粒子が時刻ｔに速度ｕを持つ条件付の確率分布Ｔ(ｕ,ｔ|ｕ₁,ｔ₁)は,(∂/∂ｔ)Ｔ(ｕ,ｔ|ｕ₁,ｔ₁)＝γ(∂/∂ｕ)[ｕ＋(ｋ_BＴ/ｍ)∂/∂ｕ＋(Ｄ_u/ｍ²)∂²/∂ｕ²]Ｔ(ｕ,ｔ|ｕ₁,ｔ₁)という方程式に従うことがわかります。

確率分布が従うこの方程式を,フォッカー・プランク方程式と呼ぶわけですね。

そして,先のａについての運動方程式ｄａ_i/ｄｔ＝{ａ_i,Ｓ}＋Σ_jＬ_ij(∂Ｓ/∂ａ_j)＋Ｒ_i(ｔ)を上の１次元速度ｕに対するランジュバン方程式におきかえ,時刻ｔに状態ａが実現する確率分布Ｐ(ａ,ｔ)を上の確率分布Ｔ(ｕ,ｔ|ｕ₁,ｔ₁)におきかえれば,先述のＰ(ａ,ｔ)に対するフォッカー・プランク方程式が得られます。

これが定常解として平衡分布Ｐ_eq(ａ)＝exp[Ｓ(ａ)/ｋ_B](あるいはこの定数倍)を持つためにはｋ_BＬ_ij＝Ｄ_ijとなることが必要条件になります。

つまり,∂Ｐ_eq(ａ)/∂ａ_j＝(1/ｋ_B)(∂Ｓ/∂ａ_j)Ｐ_eq(ａ)なので∂Ｐ(ａ,ｔ)/∂ｔ＝－Σ_i(∂/∂ａ_i){ａ_i,Ｓ}Ｐ(ａ,ｔ)＋Σ_iΣ_j(∂/∂ａ_i)[－Ｌ_ij(∂Ｓ/∂ａ_j)＋Ｄ_ij(∂/∂ａ_j)]Ｐ(ａ,ｔ)においてＰ(ａ,ｔ)＝Ｐ_eq(ａ)とおくと,ｋ_BＬ_ij＝Ｄ_ijが満たされる場合には右辺の第２項はゼロになります。

一方,第１項はΣ_i(∂/∂ａ_i){ａ_i,Ｓ}Ｐ_eq(ａ)＝Ｐ_eq(ａ)[Σ_i(∂/∂ａ_i){ａ_i,Ｓ}＋(1/ｋ_B)Σ_i(∂Ｓ/∂ａ_i){ａ_i,Ｓ}]となりますが,定義によってΣ_i(∂Ｓ/∂ａ_i){ａ_i,Ｓ}＝{Ｓ,Ｓ}＝0 です。

これほど自明ではありませんが,Σ_i(∂/∂ａ_i){ａ_i,Ｓ}＝0 も成立します。

実際,Σ_i(∂/∂ａ_i){ａ_i,Ｓ}＝Σ_i[{1,Ｓ}＋{ａ_i,∂Ｓ/∂ａ_i}]＝Σ_iΣ_i[{1,ａ_j}(∂Ｓ/∂ａ_j)＋{ａ_i,ａ_j}(∂²Ｓ/∂ａ_i∂ａ_j)]＝Σ_iΣ_i{1,ａ_j}(∂Ｓ/∂ａ_j)＝－Σ_iΣ_i,k(∂１/∂ａ_k){ａ_k,ａ_j}(∂Ｓ/∂ａ_j)＝0 となります。１という関数はδ-関数であり,∂１/∂ａ_kは汎関数微分ですね。

そこで,ｋ_BＬ_ij＝Ｄ_ijはＰ_eq(ａ)＝exp[Ｓ(ａ)/ｋ_B]が解になるための必要十分条件であることがわかりました。

それ故,"揺動力＝ゆらぎ"の時間相関関数＜Ｒ_i(ｔ)Ｒ_j(ｔ')＞(白色雑音の2Ｄは揺動力の強さと呼ばれる)と,輸送係数:Ｌ_ijの間には＜Ｒ_i(ｔ)Ｒ_j(ｔ')＞＝2ｋ_BＬ_ijδ(ｔ－ｔ')という関係が成り立ちます。

これはさらに∫₀^∞ｄ(ｔ－ｔ')＜Ｒ_i(ｔ)Ｒ_j(ｔ')＞＝ｋ_BＬ_ij：すなわち∫₀^∞ｄτ＜Ｒ_i(ｔ)Ｒ_j(ｔ＋τ)＞＝ｋ_BＬ_ijと書き直すことができます。

つまり,輸送係数あるいはオンサーガー係数:Ｌ_ijはゆらぎの時間相関関数で与えられます。この法則を揺動散逸定理と呼びます。

こうして,数式的に表現された形の定理が示されても,これが実際の自然現象において物理的にどのような意味を持つのかを理解しなければ,こうした定理の重要性を認識することはできません。

そこで,１例として熱伝導に適用してみます。

平均流速ｖがゼロの１成分流体中の熱伝導方程式はｅを単位質量当たりの内部エネルギーとして∂(ρｅ)/∂ｔ＋∇Ｊ_q＝0 で与えられます。

ここでＪ_qは熱流であり熱拡散の線形近似モデルではＪ_q＝λ∇(1/Ｔ)＝－κ∇Ｔと表わされます。κ＝λ/Ｔ²は熱伝導率と呼ばれています。

この現象論的方程式に対して,これのランジュバン方程式は∂(ρｅ)/∂ｔ＋∇Ｊ_q＝－∇ｑ(ｒ,ｔ)となります。ただしｑ(ｒ,ｔ)は熱流Ｊ_qのゆらぎです。

Ｊ_qは熱流ですから,流体の局所流速をｖ(ｒ,ｔ)とすると,Ｊ_q＝ρｅｖ＝ρｅ(ｄｒ/ｄｔ)です。

示量変数の１つとしてａ＝ρｅｒとすると,ａの不可逆変化部分に対する表式ｄａ_i/ｄｔ＝Σ_jＬ_ij(∂Ｓ/∂ａ_j)＋Ｒ_i(ｔ)は,ｄａ_i/ｄｔ＝ｄ(ρｅｒ_i)/ｄｔ～(Ｊ_q)_i＝Σ_jＬ_ij[∂Ｓ/∂(ρｅｒ_j)] ＋Ｒ_i(ｔ)と書けますが,平衡状態ではρＶｄｅ＝ＴｄＳ－ＰｄＶなので体積一定(ｄＶ＝0 )ならＳ＝ρｅＶ/Ｔ＋(定数)です。

それ故,結局(Ｊ_q)_i～Σ_jＬ_ijＶ(∂/∂ｒ_j)(1/Ｔ)となります。これをＪ_q＝λ∇(1/Ｔ)：すなわち(Ｊ_q)_i＝λ(∂/∂ｒ_j)(1/Ｔ)と比較すると,輸送係数＝オンサーガー係数についてＬ_ij＝(λ/Ｖ)δ_ijという表式が得られます。

そこで揺動散逸定理によると,λは平衡状態における揺動熱流ｑ(ｒ,ｔ)の時間相関関数によって表現されることになります。すなわち,揺動熱流ｑ(ｒ,ｔ)の時間相関関数は,∫₀^∞ｄｔ＜ｑ_α(ｒ,ｔ)ｑ_β(ｒ',0)＞_eq＝ｋ_B(λ/Ｖ)δ_αβδ(ｒ－ｒ')と書けます。

そして,熱流Ｊ_q(ｒ,ｔ)のゆらぎｑ(ｒ,ｔ)を体積Ｖの対象領域全体で空間積分したｑ(ｔ)＝∫ｄｒｑ(ｒ,ｔ)という量を定義して,これを時刻ｔでの"熱流のゆらぎ"と呼ぶことにすれば,熱伝導に対する揺動散逸定理の表現は∫₀^∞ｄｔ＜ｑ_α(ｔ)ｑ_β(0)＞_eq＝ｋ_Bλδ_αβと書くことができます。

今は,平均流速がゼロの流体を考えており,平衡状態では＜Ｊ_q(ｒ,ｔ)＞_eq＝0 なので,熱流のゆらぎｑ(ｒ,ｔ)がエネルギー流そのものになります。そこで対流がない流体においての平衡状態ではｑ(ｔ)はエネルギー流を空間全体で積分したもののゆらぎとなります。

この例では,輸送係数＝オンサーガー係数の一種である熱伝導率がミクロな流れｑ(ｔ)の時間相関関数で与えられるということが揺動散逸定理からの重要な帰結と言えます。

19日木曜日夜から風邪気味で,症状自体は軽いのですがあまり筆が進みません。

北原和夫 著「非平衡系の統計力学」(岩波書店)

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

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

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

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

    物理学

ベナール対流の安定性とレイリー数

　今日は上下にそれぞれ一様な温度の床と天井があるところに閉じ込められた流体があって,元々はその高さに比例した温度変化と圧力変化を持って静止していた流体の,上下についての対流の発生,および安定性とレイリー数(Rayleigh number)との関係について述べたいと思います。

　これは,いわゆるベナール対流(Benard's convection)の問題です。

　そのために,まず流体方程式系の非圧縮近似について考察するところから始めましょう。

　流体の速度が音速に比較して非常に小さいときは非圧縮性流体:Ｄρ/Ｄｔ＝0 としてよいわけですが,特に重力による圧力傾度とつりあった静力学平衡:∂Ｐ/∂ｚ＝－ρｇを満足するような密度ρ＝ρ₀ からの温度分布などによるずれ:Δρだけを圧縮性として考慮して,それ以外は非圧縮性流体として扱う近似が正当化されます。

　これは気象力学などではブジネスク(Boussinesq)近似としてよく利用される近似です。準備として,まずその定式化をすることにします。

　"質量保存の方程式＝連続の方程式"は(∂ρ/∂ｔ)＋div(ρｖ)＝ 0 ですが,これをラグランジュ微分で表すと,Ｄρ/Ｄｔ＋ρdivｖ＝0 です。そして非圧縮性流体Ｄρ/Ｄｔ＝0 の近似を行なうと,連続の方程式はdivｖ＝0 となります。(Ｄ/Ｄｔ≡∂/∂ｔ＋ｖ∇です。)

　次に,運動方程式は初めから圧力項を除く部分では非圧縮 divｖ＝0 を仮定し,外力は重力だけであるとするナビエ・ストークスの方程式(Navier-Stokes equation)とします。

　すなわち,Ｄｖ/Ｄｔ＝－∇Ｐ/ρ－ｇｋ＋ν∇²ｖですね。

　ここで,ρ＝ρ₀＋Δρ,Ｐ＝Ｐ₀＋Ｐ'とし,∇Ｐ₀＝－ρ₀ｇｋとして運動方程式の右辺の－∇Ｐ/ρを展開して２次の微小量を無視します。

　すると,－∇Ｐ/ρ＝－∇(Ｐ₀＋Ｐ')/(ρ₀＋Δρ)＝－(1/ρ₀)[１－(Δρ/ρ₀)]∇(Ｐ₀＋Ｐ')＝[1－(Δρ/ρ₀)]ｇｋ－∇Ｐ'/ρ₀です。

　ここで,ｇは重力の加速度,νは動粘性係数＝η/ρです。

　そこで,運動方程式はＤｖ/Ｄｔ＝－∇Ｐ'/ρ₀－(Δρ/ρ₀)ｇｋ＋ν∇²ｖとなり,これがブジネスク近似です。

　温度Ｔは密度と圧力で決まるので,温度による表現も与えておきましょう。ρ₀,Ｐ₀に対応する温度をＴ₀とし,Ｔ＝Ｔ₀＋ΔＴと置くと,圧縮率をαとしてρ＝ρ₀(1－αΔＴ)となり,Δρ/ρ₀＝－αΔＴですから,運動方程式は,Ｄｖ/Ｄｔ＝－∇Ｐ'/ρ₀＋αｇΔＴｋ＋ν∇²ｖと変形できます。

　最後にエネルギー方程式ですが,これは温度の方程式に直すことができます。粘性による散逸エネルギーを無視すれば,これはＤＴ/Ｄｔ＝κ∇²Ｔで与えられます。

　ここにκ＝ｋ/(ρＣ_p)は温度拡散係数です。ｋは熱伝導係数であり,Ｃ_pは定圧比熱です。

　以下では,ｖ＝(ｕ,ｖ,ｗ)として話を進めます。まず最初の設定としては領域全体でｕ＝ｖ＝ｗ＝0 とします。

　そして,温度は床ｚ＝0 ではＴ＝Ｔ₀,天井ｚ＝ＨではＴ＝Ｔ₁で時間的に変わらないとします。この境界条件を満たす定常熱伝導の方程式∇²Ｔ＝0 の解はＴ(ｚ)＝Ｔ₀－βｚ＝Ｔ₀＋{(Ｔ₁－Ｔ₀)/Ｈ}ｚですね。

　ここで,β＞0 ,すなわち,Ｔ₀＞Ｔ₁とします。なぜなら,そうでないと地球上の重力場の中では対流が起きないからです。

　そして,そればかりでなく|β|＝－βの値が,いわゆる断熱温度減率(これは湿度によって変化しますが)より大きくないと対流は起きません。(例えばランダウ著「流体力学１」(東京図書)参照）

　そして,圧力はＰ(ｚ)＝－ρｇｚ＋const.とします。ｖ＝(ｕ,ｖ,ｗ)を改めてｕ＝ｖ＝ｗ＝0 の平衡状態からの摂動流速とみなし,温度はＴ≡Ｔ(ｚ)＋θ≡Ｔ－βｚ＋θ,圧力もＰ≡Ｐ(ｚ)＋Ｐ'と摂動温度,摂動圧力を取ります。

　連続方程式は∂ｕ/∂ｘ＋∂ｖ/∂ｙ＋∂ｗ/∂ｚ＝0 であり,またブジネスク近似での運動方程式はＤｕ/Ｄｔ＝－(∂Ｐ'/∂ｘ)/ρ₀＋ν∇²ｕ,Ｄｖ/Ｄｔ＝－(∂Ｐ'/∂ｙ)/ρ₀＋ν∇²ｖ,Ｄｗ/Ｄｔ＝－(∂Ｐ'/∂ｚ)/ρ₀＋αｇθ＋ν∇²ｗとなります。

　温度方程式はＤＴ/Ｄｔ＝κ∇²Ｔですが,ＤＴ/Ｄｔの中の移流項ｖ∇Ｔは主流を取って,ｗ∂Ｔ/∂ｚ＝－βｗのみですから,∂θ/∂ｔ＝βｗ＋κ∇²θとなります。

　ここで,運動方程式からｕとｖを消去するために,ｖ＝(ｕ,ｖ,ｗ)の"回転＝渦"：ω＝rotｖを取ります。ζ＝∂ｖ/∂ｘ－∂ｕ/∂ｙとすると,運動方程式から∂ζ/∂ｔ＝ν∇²ζが得られ,また(∂/∂ｔ)∇²ｗ＝αｇ(∇²－∂²/∂ｚ²)θ＋ν∇⁴ｗも得られます。

　こうして基本方程式系としては,ｗ,θ,ζの３つを未知数とする３つの微分方程式が得られたことになります。

　さらに,ｘ,ｙ,ｔの独立変数をｚから分離するため,変数分離形を仮定して,ｗ≡Ｗ(ｚ)exp[i(ｋ_xｘ＋ｋ_yｙ)＋ｐｔ],θ≡Θ(ｚ)exp[i(ｋ_xｘ＋ｋ_yｙ)＋ｐｔ],ζ≡Ζ(ｚ)exp[i(ｋ_xｘ＋ｋ_yｙ)＋ｐｔ]と置きます。

　これは方程式系において∂/∂ｔ＝ｐ,∇²－∂²/∂ｚ²＝－ｋ²≡－(ｋ_x²＋ｋ_y²),∇²＝ｄ²/ｄｚ²－ｋ²と置き換えることに相当します。

　すなわち,ｐ(ｄ²/ｄｚ²－ｋ²)Ｗ＝－αｇｋ²Θ＋ν(ｄ²/ｄｚ²－ｋ²)²Ｗ,ｐΘ＝βＷ＋κ(ｄ²/ｄｚ²－ｋ²)Θ,ｐΖ＝ν(ｄ²/ｄｚ²－ｋ²)Ζとなります。

　これは３つめの方程式とは別に,前の２つのＷとΘだけで独立しているのでこれら２つだけでも解ける形になっています。

　そして,床と天井での境界条件はｚ＝0 とｚ＝ＨでΘ＝Ｗ＝Ζ＝0 ,ｄＷ/ｄｚ＝0 です。

　ここで,これらの方程式を無次元化するために単位を[Ｌ]＝Ｈ,[Ｔ]＝Ｈ²/νに取ることにして,ａ≡ｋＨ,σ≡ｐＨ²/νとします。

　ｘ,ｙ,ｚの単位もＨとした後に,Ｄ≡ｄ/ｄｚと定義し,温度拡散のプラントル数(Prandtl number)をＰｒ≡ν/κとします。

　すると,３つのうちの前の２つの方程式は(Ｄ²－ａ²)(Ｄ²－ａ²－σ)Ｗ＝(αｇＨ²/ν)ａ²Θ,および(Ｄ²－ａ²－Ｐｒ・σ)Θ＝－(βＨ²/κ)Ｗとなります。

　そして境界条件は,ｚ＝0 とｚ＝1 においてΘ＝Ｗ＝Ζ＝0 ,ＤＷ＝0 です。(ただしＷとΘは通常の単位であり無次元化していません。)

　Θを消去すると,(Ｄ²－ａ²)(Ｄ²－ａ²－σ)(Ｄ²－ａ²－Ｐｒ・σ)Ｗ＝－Ｒａ²Ｗとなります。Ｒはレイリー数で,これはＲ≡ｇαβＨ⁴/(κν)で定義される無次元数であって,方程式を支配する特徴的な数です。

　安定性を調べるためにはσ＝0 (つまりｐ＝0)の特別なケースを調べればよいと思われます。なぜなら時間的に減衰するか,増幅するかの限界は,exp(ｐｔ)という因子でのｐ＝0 を臨界点として生じるからです。

　そして,その境界におけるレイリー数に最初に達したときに流れは初めて不安定になると考えられるからです。

そして対称性を利用するために,ｚの原点を移動して,ｚ＝±1/2を境界に取ります。

基本となる微分方程式は(Ｄ²－ａ²)³Ｗ＝－Ｒａ²Ｗであり,境界条件はＷ＝ＤＷ＝(Ｄ²－ａ²)²Ｗ＝0 atｚ＝±1/2です。

この方程式の特性方程式はＷ≡ｅ^qzとして,(ｑ²－ａ²)³＝－Ｒａ²であり,Ｒａ²＝τ³ａ⁶と置いて解くとｑ²＝－ａ²(τ－1),ｑ²＝ａ²[1＋(1/2)τ(1±ｉ√3)],あるいは平方根を取って,±iｑ₀,±ｑ,±ｑ^* (ｑ₀＝ａ(τ－1)^1/2 etc.)ですね。

　あとの詳細は簡単なので本質的な部分以外は省略しますが,解は偶関数解と奇関数解の２種類に分かれ,偶関数解としてＷ＝Ａ₀cos(ｑ₀ｚ)＋Ａcosh(ｑｚ)＋Ａ^*cosh(ｑ^*ｚ),奇関数解としてＷ＝Ａ₀sin(ｑ₀ｚ)＋Ａsinh(ｑｚ)＋Ａ^*sinh(ｑ^*ｚ)が得られます。

　後は,境界ｚ＝±1/2でＷ­＝ＤＷ＝(Ｄ²－ａ²)²Ｗ＝0 である,という条件からＡ₀,Ａ,Ａ^*を決めればいいわけです。

　これらが自明でない解を持つための条件は境界条件を満足する係数ベクトル:Ａ₀,Ａ,Ａ^*の係数行列の行列式がゼロでなければならないという1つの固有値方程式になります。

　例えば偶関数解ではこの行列式＝0 の条件は,Im[(i＋√3)ｑtanh(ｑ/2)]＋ｑ₀tanh(ｑ₀/2)＝0 となります。

　こうした固有値問題を満足するＲを求めると,流れが不安定になり熱伝導だけでは熱を運ぶことが不可能となって物質で熱を運ぶ対流が激しくなる臨界点が求まり,最後には乱流へと移行することになる境界の臨界レイリー数が求められます。

　私がかつて数値的に解いた偶関数の臨界レイリー数Ｒはたとえばａ＝3.0でＲ＝1710 etc.であり文献と一致しています。

　つまり,ａ＝3.0ではこのＲを境として不安定流となり対流が激化していくことになるわけですね。

　参考文献：S.Chandrasekhar 「Hydromagnetic and Hydromagnetic stability」(Dover Books),モーニン,ヤグロム(A.S.Monin&A.M.Yaglom)著(山田豊一訳)「統計流体力学1」(総合図書)

http://fphys.nifty.com/(ニフティ「物理フォーラム」サブマネージャー）                                       TOSHI 

http://blog.with2.net/link.php?269343(ブログ・ランキングの投票）↑ここをクリックすると投票したことになります。

人口増加とロジスティック曲線

　今日は軽い話題を一つ述べましょう。

　現在の人口をN 人としΔｔ 年間に人口 N 人に比例して (ｋΔｔ)Ｎ人だけ人口が増加するとします。今の時刻（年）を ｔ とし増加率を ｋ ＝ｂ－ｄ （ｂは出生率,ｄは死亡率）とします。

　ｋは一定であるとすると, ｔ＋Δｔ の人口:N₁= N＋ΔN人は,ΔN /Ｎ＝ｋΔｔによってＮ₁＝Ｎ (１＋ｋΔｔ）となりますから,時刻（年） ｔ＋nΔｔでの人口:N_n人はN_n＝N（１＋ｋΔｔ)ⁿ となります。ｋ＞0 であれば,まさに人口はネズミ算的に増えてゆきますね。

　そしてΔｔ が無限小ｄｔであるとすると,ΔN /Ｎ＝ｋΔｔはｄＮ/ｄｔ＝ｋＮ となりますから ,ｔ＝0 での人口をＮ＝Ｎ₀ 人とするとＮ＝Ｎ(ｔ)＝N₀ ｅｘｐ ( ｋｔ )ですね。これを見ると,ｋ＞0 なら ｔ → ∞ でＮ → ∞ ですが, ｋ＜0 なら ｔ → ∞ で Ｎ → 0 であり絶滅してしまいます。

　しかし,実際には Δｔ の間にいろいろな災害や環境の変化などで増加率 ｋ は一定ではなく,かなりの変化を受けると考えられます。

　一般に人口 Ｎが増えれば増えるほど個体の増加は妨げられる傾向がありますから,それは増加率がｋから ｋ（１－αＮ） (α＞0 )になるという影響を与えるという効果で表わすことができます。

　このモデルをロジスティックモデルといいます。ｋ（１－αＮ） (α＞0 )において,αＮ＞１なら人口Ｎは増加し,逆なら人口Ｎは減少しますね。

　このモデルでは,Ｎに対する微分方程式という形ではｄＮ/ｄｔ＝ｋＮ（１－αＮ）という非線形微分方程式になります。これを解くとＮ＝Ｎ(ｔ)＝Ｎ₀ /[αＮ₀ { 1－ｅｘｐ ( - ｔ ）}＋exp (－ｔ ) ] です。

　これの描くＮ-ｔ曲線がロジスティック曲線と呼ばれるものです。これを見ると ｔ → ∞ の極限ではＮ → １/αとなって,人口Ｎはある極限値に到達するわけです。それ以上は増加も減少もしません。

　ロジスティックモデルは実際に生態学(ecology)で個体の増加減少の履歴と一致する例が多々あり,人口にもこれが適用できると考えられます。

　正に「増え過ぎた生物は抑制される。」というのが自然の摂理ですが,人類は天敵がいないことや医学の進歩,(そして軍縮などによる戦争の減少？)などによって自然の摂理を破壊し生態系を破壊しつつあります。やがてこの自然の摂理の破壊の報いを受けるかもしれません。

　ところで,ロジスティック微分方程式のΔｔ の刻みを調節して中心差分の差分方程式として離散化すると ,ｋ の値によっては,ｔ が大きいところで不安定な人口増減の振動をするカオスの現象を起こすことが知られています。

　この不安定性は数値解析の目的で「離散化＝差分化」を行なったために生じたものですが,現実の現象のモデルとしては時間刻みが無限小の微分方程式よりも時間刻み有限の差分方程式の方が適切かもしれません。

　カオスの例としては,ロジスティック模型のｘ_n+1＝ａｘ_n(１－ｘ_n)は典型的な例であり,a の値によっては「リー・ヨーク(Li-Yorke)の定理」においてカオスになる条件が満足されます。

（参考文献；山口昌也著「カオス入門」（朝倉書店）、山口昌也編著「数値解析と非線型現象」（日本評論社）　）

http://fphys.nifty.com/(ニフティ「物理フォーラム」サブマネージャー）                                       TOSHI 

http://blog.with2.net/link.php?269343(ブログ・ランキングの投票）

↑ここをクリックして投票をお願いします。

ここをクリックして投票をお願いします。