プルーム上昇のモデル式(2)

　今日もプルーム上昇式の続きです。

§2.理論の単純化,エントレインメントの仮定

　これまでの議論では,未知数Ｆ_z,ｗ,Ｖに対して２つの方程式しか与えられておらず,問題が閉じていません。

　これは,プルームの明確な定義が与えられていないため,プルームの大きさを与える式が無いことに起因すると考えられます。

　そこで,連続の方程式∇(ρ_pｖ_p)＝0 を用いてプルームの大きさを規定する,ある方程式を与えます。

　プルーム半径をｒとして,改めてプルームの平均進行速度の大きさをｖ_sとします。さらに,エントレインメント(entrainment)速度,つまり,プルーム境界から乱流により流入する外気の平均速さをｖ_eとします。

　そして,プルームの行路長さ:ｓに対してｄｓ≡ｖ_sｄｔ,またはｄｔ≡ｄｓ/ｖ_sによってｄｔを定義しておきます。

　すると,質量の保存式から明らかにΔ(ρ_pπｒ²ｖ_s)＝ρ_a2πｒｖ_eΔｓ,または(ｄ/ｄｓ)(ρ_pｒ²ｖ_s/ρ_a)＝2ｒｖ_eです。こうした関係があることをエントレインメントの仮定といいます。

　特に,鉛直プルーム:ｓ＝ｚ,ｄｚ＝ｗｄｔの場合にはＶ＝ρ_pπｒ²ｗ/(πρ_a)＝ρ_pｒ²ｗ/ρ_aなのでエントレインメントの仮定はｄＶ/ｄｚ＝2ｒｖ_eを意味します。

　また,(ｄ/ｄｓ)(ρ_pｒ²ｖ_s/ρ_a)＝2ｒｖ_eよりρ_aはｓに依らないとして(ｒ²/ρ_a){ｄ(ρ_pｖ_s)/ｄｓ}＋(2ｒρ_pｖ_s/ρ_a)(ｄｒ/ｄｓ)＝2ｒｖ_eです。

さらに,ｄｔ＝ｄｓ/ｖ_sよりｄ/ｄｓ＝(1/ｖ_s)(ｄ/ｄｔ)を用いればｄｒ/ｄｔ＝ρ_aｖ_e/ρ_p－{ｒ/(2ρ_p)}{ｄ(ρ_pｖ_s)/ｄｓ}を得ます。

ここでTaylorによる自然な仮定:ｖ_e∝ｗを採用し,鉛直プルームに対するエントレインメント係数αを導入してｖ_e＝αｗとします。

また,再びＶ＝ρ_pπｒ²ｗ/(πρ_a)＝ρ_pｒ²ｗ/ρ_aよりｒ²＝ρ_aＶ/(ρ_pｗ)～Ｖ/ｗと近似できるという考察から,逆にプルーム半径:ｒをｒ≡(Ｖ/ｗ)^1/2によって定義します。

この定義によれば,Ｖ＝ｒ²ｗですが,Δ(ρ_pπｒ²ｖ_s)＝ρ_a2πｒｖ_eΔｓ,または(ｄ/ｄｓ)(ρ_pｒ²ｖ_s/ρ_a)＝2ｒｖ_eの関係は変わらず,ρ_p～ρ_aの近似で考えると依然としてエントレインメントの仮定はｄＶ/ｄｚ＝2ｒｖ_eです。

元々,これも仮定ですから以下ではｄＶ/ｄｚ＝2ｒｖ_eをエントレインメントの仮定として採用します。この右辺にｖ_e＝αｗ,ｒ＝(Ｖ/ｗ)^1/2を代入するとｄＶ/ｄｚ＝2α(ｗＶ)^1/2を得ます。

以上から,(抵抗力)を無視すると次のような簡単な基本方程式系が得られました。

(4)ｄＦ_z/ｄｚ＝－ｓＶ,(12)ｄ(ｗＶ)/ｄｚ＝Ｆ_z/ｗ,

および,(13)ｄＶ/ｄｚ＝2α(ｗＶ)^1/2　です。

一応,未知量Ｖ,Ｆ_z,ｗの意味を再び述べておきます。

Ｖは鉛直容積流束,Ｆ_zは鉛直方向の浮力(ｇθ'/Ｔ)の流束,ｗはプルーム面の鉛直速度成分です。

プルームの行路長さと混同するおそれのある同じ変数ｓを用いていますが,(4)式のｓは,もちろんｓ≡(ｇ/Ｔ_a)(∂θ_a/∂ｚ)で定義される安定度パラメータです。

そして,Morton,Taylor,Turner(1956)によれば,エントレインメント係数αの具体的な値の候補はα＝0.093です。一方,Briggs(1968)によればα＝0.075です。また,Ricou,Spalding(1961)によればジェットに対しα＝0.080です。

　§3.折れ曲がりプルーム(bent-over plume)

　進行方向が完全に鉛直ｚ方向でｓ＝(一定)のプルーム断面が水平面内にある鉛直プルームとは異なる状態のプルームを折れ曲がりプルームといいます。

折れ曲がりプルームの場合,積分面をプルームと交わる鉛直面として,この面上ではｓ＝(ｇ/Ｔ)(∂θ_a/∂ｚ),およびｖ_a＝ｕｉは一定であると仮定します。

前と同様,∫_∞∇(ρ_pθ'ｖ_p)ｄｙｄｚ＝－(∂θ_a/∂ｚ)∫_∞ρ_pｗ'ｄｙｄｚですがプルームの外ではθ'＝0,ｗ'＝0 とすれば,これは(ｄ/ｄｘ)∫_Pρ_pθ'(ｕ＋ｕ')ｄｙｄｚ＝－(∂θ_a/∂ｚ)∫_Pρ_pｗ'ｄｙｄｚを意味します。

そこで,浮力流束を(5)Ｆ_z≡∫_P(ｇρ_pθ'ｗ'/Ｔ_a)ｄｘｄｙ/(πρ_a)の代わりに(14)Ｆ_z≡∫_P{ｇρ_pθ'(ｕ＋ｕ')/Ｔ_a}ｄｙｄｚ/(πρ_a)で定義すれば,(15)ｄＦ_z/ｄｘ＝－ｓＭと表わすことができます。

ただし,再びｓはｓ＝(ｇ/Ｔ_a)(∂θ_a/∂ｚ)でありＭはＭ≡∫_Pρ_pｗ'ｄｙｄｚ/(πρ_a)で定義される量です。

この場合も排出源がhot-sourceなら,浮力流束Ｆ_zの初期値Ｆは{(単位時間にプルーム積分面を通過する総排出熱量)×(ｇ/Ｔ_a)/ｃ_p}/(πρ_a)と考えられるので,やはりＦ＝ｇＱ_H/(πρ_aｃ_pＴ_a)で与えられます。

　そして,運動量に対する方程式は,(8)ｄ(ｖＶ)/ｄｚ＝(Ｆ_z/ｗ)ｋ＋ｖ_a(ｄＶ/ｄｚ)＋(抵抗力)の証明の際に遭遇した式:(ｄ/ｄｎ)(∫_Pρ_pｖ_pnｖ_piｄＳ)＝{1＋(ｄｒ/ｄｎ)²}^1/2(∫_Cρ_pｖ_pｖ_piｄｓ)＋δ_i3∫_P(ｇρ_pθ'/Ｔ_a)ｄＳにおいて進行方向ｎをｉに置くと得られます。

　ｄ/ｄｎをｄ/ｄｘに変え,ｖ_p＝(ｕ＋ｕ')ｉ＋ｖ'ｊ＋ｗ'ｋ,ｄｓ＝{1＋(ｄｒ/ｄｘ)²}^-1/2(－ｄｚｊ＋ｄｙｋ)＋{1＋(ｄｒ/ｄｘ)²}^-1/2(ｄｒ/ｄｘ)ｄｌｉとして,必要なｉ＝3の場合のみを書き下します。

　それは,(ｄ/ｄｘ){∫_Pρ_p(ｕ＋ｕ')ｗ'ｄＳ}＝∫_P(ｇρ_pθ'/Ｔ_a)ｄＳ＋{1＋(ｄｒ/ｄｘ)²}^1/2(∫_Cρ_pｖ_pｗ'ｄｓ)＝∫_P(ｇρ_pθ'/Ｔ_a)ｄＳ－∫_Cρ_pｖ'ｗ'ｄｚ＋∫_Cρ_pｗ'²ｄｙ＋∫_Cρ_p(ｕ＋ｕ')ｗ'(ｄｒ/ｄｘ)ｄｌです。

　ｕ'＜＜ｕなのでｕ＋ｕ'～ ｕと近似して有効運動量流束をｕＭ_eff;Ｍ_eff≡{∫_Pρ_pｗ'ｄＳ－∫ｄｘ∫_Cρ_pｗ'(ｄｒ/ｄｘ)ｄｌ}/(πρ_a)と定義しｕ',ｖ',ｗ'の２次以上の項を無視します。

　すると,Ｆ_z＝{ｕ∫_P(ｇρ_pθ'/Ｔ_a)ｄＳ}/(πρ_a)より,結局運動方程式として(16)ｕ(ｄＭ_eff/ｄｘ)＝Ｆ_z/ｕが得られます。

　鉛直運動量の水平方向流束ｕＭ_effの初期値はプルーム面を通過する鉛直運動量ρ_pｗπｒ²/(πρ_a)＝ρ_pｗｒ²/ρ_aの総流束(鉛直方向流束)ですから,その初期値はＦ_m＝(ρ₀ｗ₀ｒ₀²/ρ_a)×ｗ₀＝ρ₀ｗ₀²ｒ₀²/ρ_a＝ρ₀ｗ₀Ｖ₀/ρ_aです。

そして,Briggs(1972),Richards(1963)によれば有効運動量流束ｕＭ_effの運動量流束ｕＭ;Ｍ＝∫_Pρ_pｗ'ｄｙｄｚ/(πρ_a)に対する比率の評価値はＭ_eff/Ｍ ～ 2.3です。

ｚをプルーム軸心の高さ,ｄｔ≡ｄｘ/ｕとしてｗ≡ｄｚ/ｄｔ＝ｕｄｚ/ｄｘと定義します。ｗはプルーム軸心の上昇速度を与え,解(ｘ,ｚ),or ｚ＝ｚ(ｘ)はプルームの軌跡を示します。

ここで,折れ曲がりプルームにおけるエントレインメント係数をβとします。すなわちｖ_e＝βｗです。有効プルーム半径をｒとすると,Ｍ_eff＝πｒ²ρ_pｗ/(πρ_a)＝ｗρ_pｒ²/ρ_aです。

折れ曲がりプルームの場合のエントレインメントの仮定はｄ(ρ_pπｒ²ｕ)/ｄｘ＝2πｒρ_aｖ_e,つまりｄ(ρ_pｒ²ｕ/ρ_a)/ｄｘ＝2ｒｖ_eと表わせます。

そこで,ｕｄ(Ｍ_eff/ｗ)/ｄｘ＝2ｒｖ_e＝2ｒβｗです。

ｒ≡(Ｍ_eff/ｗ)^1/2と定義すれば,(17)ｕｄ(Ｍ_eff/ｗ)/ｄｘ＝2β(ｗＭ_eff)^1/2を得ます。

さらにｗＶ≡ｕＭ_effによってＶを定義すれば,Ｖは水平方向の容積流束でありｒ＝(Ｖ/ｕ)^1/2です。これとｕ/ｄｘ＝ｗ/ｄｚなる関係を用いて基本方程式系:(15),(16),(17)を書き直せば次のようになります。

すなわち,(18)ｄＦ_z/ｄｚ＝－ｓ'Ｖ,(19)ｓ'＝ｓ(Ｍ/Ｍ_eff),(20)ｄ(ｗＶ)/ｄｚ＝Ｆ_z/ｗ,(21)ｄＶ/ｄｚ＝2β(ｕＶ)^1/2です。

式(21):ｄＶ/ｄｚ＝2β(ｕＶ)^1/2は,ｚ＝0 でＶ＝0 という点源(point source)の初期条件に対しては容易に解けて,(Ｖ/ｕ)^1/2＝βｚ,つまりｒ＝βｚです。

Richards(1963),TVA-Report(1967)によれば折れ曲がりプルームのエントレインメント係数はβ～ 0.5です。

§4.プルーム上昇高度の算出

　具体的にプルーム上昇式を求めるに際して,さらに次のような２つの仮定を与えます。

(仮定Ⅰ):鉛直プルームを満足する状態を無風状態(静穏:calm)に対応させ,折れ曲がりプルームを有風時の状態に対応させる。

(仮定Ⅱ):初期時には鉛直運動量か浮力流束のいずれか一方が支配的である。前者はジェットと呼んでＦ＝0 とし,後者を浮力プルームと呼んでＦ_m＝0 とする。

　さて,まずｄｘ＝ｕｄｔ,ｄｚ＝ｗｄｔにより,全ての方程式を独立変数ｔの式に変換します。

鉛直プルームについては,(4),(12),(13)より

(22)ｄＦ_z/ｄｔ＝－ｓ(ｗＶ),(23)ｄ(ｗＶ)/ｄｔ＝Ｆ_z,(24)ｄＶ/ｄｔ＝2α(ｗＶ)^3/2/Ｖです。

また,折れ曲がりプルームについては,(18),(20),(21)より

(25)ｄＦ_z/ｄｔ＝－ｓ'(ｗＶ),(26)ｄ(ｗＶ)/ｄｔ＝Ｆ_z,(27)ｄＶ/ｄｔ＝2βｗＶ(ｕ/Ｖ)^1/2です。

　以下,順に安定大気:ｓ＞0,ｓ'＞0 (ｓ,ｓ'は一定値を仮定)の中でのプルーム上昇高さを求めます。

(ⅰ)鉛直プルーム(無風:calm)

　このケースの基本方程式系は(22)ｄＦ_z/ｄｔ＝－ｓ(ｗＶ),(23)ｄ(ｗＶ)/ｄｔ＝Ｆ_z,(24)ｄＶ/ｄｔ＝2α(ｗＶ)^3/2/Ｖです。

(22),(23)よりｄ²(ｗＶ)/ｄｔ²＝－ｓ(ｗＶ)を得ます。これはよく知られている線形振動の方程式です。

そこで,ω≡ｓ^1/2としてこれをｄ²(ｗＶ)/ｄｔ²＝－ω²(ｗＶ)と書き,初期条件:ｔ＝0 でｗＶ＝Ｆ_m,ｄ(ｗＶ)/ｄｔ＝Ｆ_z＝Ｆの解を求めると,ｗＶ＝Ｆ_mcosωｔ＋(Ｆ/ω)sinωｔを得ます。

鉛直運動量流束ｗＶがゼロとなる最初のｔ＝ｔ₁がプルーム上昇高さを与えると考えられます。

上昇高さはΔｈ＝∫₀^t1ｗｄｔですから,仮定Ⅱの浮力プルーム(Ｆ_m＝0)ではｔ₁＝π/ω,よりΔｈ＝ω^-1∫₀^πｗｄλです。

同様にジェットプルーム(Ｆ＝0)ではｔ₁＝π/(2ω)よりΔｈ＝ω^-1∫₀^π/2ｗｄλ(λ＝ωｔ)です。

①     浮力プルーム(Ｆ_m＝0)の場合

　ｄ(Ｖ²/2)/ｄｔ＝Ｖ(ｄＶ/ｄｔ)＝2α(ｗＶ)^3/2の右辺にｗＶ＝(Ｆ/ω)sinωｔを代入してｄ(Ｖ²/2)/ｄｔ＝2α(Ｆ/ω)^3/2sin^3/2ωｔです。

故に,Ｖ²/2＝2α(Ｆ/ω)^3/2ω^-1{∫^ωtsin^3/2λｄλ＋const.},すなわちＶ＝2α^1/2Ｆ^3/4ω^-5/4{∫^ωtsin^3/2νｄν＋const.}^1/2です。

　ｗ＝ｗＶ/Ｖ＝(Ｆ/ω)sinωｔ/Ｖより,Δｈ＝ω^-1∫₀^πｗｄλ＝(1/2)α^-1/2Ｆ^1/4ω^-3/4[∫₀^πsinλｄλ{∫₀^λsin^3/2νｄν＋const.}^-1/2]を得ます。

したがって,ω＝ｓ^1/2より定数ｃをｃ≡(1/2)α^-1/2[∫₀^πsinλｄλ{∫₀^λsin^3/2νｄν＋const.}^-1/2]で定義すれば,Δｈ＝ｃＦ^1/4ｓ^-3/8なる関係式が成立することになります。

　ｔ＝0 のとき,Ｖ＝0 の点源(const.＝0)ではMorton,Taylor,Turner(1956)の数値計算があり,その示唆するところによればα＝0.132としてｃ～ 5.0ですから,浮力プルーム上昇高さの無風状態での評価式として(28)Δｈ＝5.0Ｆ^1/4ｓ^-3/8を得ます。

②     ジェットプルーム(Ｆ＝0)

　ｗＶ＝Ｆ_mcosωｔよりｄ(Ｖ²/2)/ｄｔ＝Ｖ(ｄＶ/ｄｔ)＝2α(ｗＶ)^3/2の右辺にｗＶ＝(Ｆ/ω)sinωｔを代入すれば,ｄ(Ｖ²/2)/ｄｔ＝2αＦ_m^3/2cos^3/2ωｔです。

故に,Ｖ²/2＝2αＦ_m^3/2ω^-1{∫^ωtcos^3/2λｄλ＋const.},すなわちＶ＝2α^1/2Ｆ_m^3/4ω^-1/2{∫^ωtcos^3/2νｄν＋const.}^1/2です。

　そして,ｗ＝ｗＶ/Ｖ＝Ｆ_mcosωｔ/Ｖより,Δｈ＝ω^-1∫₀^π/2ｗｄλ＝(1/2)α^-1/2Ｆ_m^1/4ω^-1・2[∫₀^π/2cosλｄλ{∫₀^λcos^3/2νｄν＋const.}^-1/2]です。

したがって,ω＝ｓ^1/2より定数ｃ₀をｃ₀≡(1/2)α^-1/2[∫₀^π/2cosλｄλ{∫₀^λcos^3/2νｄν＋const.}^-1/2]と定義すれば,Δｈ＝ｃ₀Ｆ_m^1/4ｓ^-1/4なる関係の成立することがわかります。

　ｔ＝0 のとき,Ｖ＝0 の点源(const.＝0)では,μ＝π/2－νと変数変換すれば,cosν＝sinμより∫₀^λcos^3/2νｄν＝－∫_π/2^π/2－λsin^3/2μｄμ,ｃ₀＝(1/2)α^-1/2[∫₀^π/2cosλｄλ{∫_π/2-λ^π/2sin^3/2μｄμ}^-1/2]＝(1/2)α^-1/2[∫₀^π/2sinσｄσ{∫₀^σsin^3/2μｄμ}^-1/2]となります。

　点源のとき,ｃ＝(1/2)α^-1/2[∫₀^πsinλｄλ{∫₀^λsin^3/2νｄν＋const.}^-1/2]に対してｃ₀＝(1/2)α^-1/2[∫₀^π/2sinλｄλ{∫₀^λsin^3/2νｄν＋const.}^-1/2]であり,ほとんどｃと同じです。

ｃを求めたのと同じ数値計算によれば,α＝0.132に対してｃ₀～ 4.0ですから,無風状態でのジェットプルームの評価式として(29)Δｈ＝4.0(Ｆ_m/ｓ)^1/4を得ます。

(ⅱ)折れ曲がりプルーム(有風時)

　このケースの基本方程式系は(25)ｄＦ_z/ｄｔ＝－ｓ'(ｗＶ),(26)ｄ(ｗＶ)/ｄｔ＝Ｆ_z,(27)ｄＶ/ｄｔ＝2βｗＶ(ｕ/Ｖ)^1/2です。

(25),(26)よりｄ²(ｗＶ)/ｄｔ²＝－ｓ(ｗＶ)なのでω'＝ｓ'^1/2としてｗＶ＝Ｆ_mcosω'ｔ＋(Ｆ/ω')sinω'ｔです。

これを(27)から得られるｄ(2Ｖ^3/2/3)/ｄｔ＝2βｕ^1/2ｗＶの右辺に代入すると,ｄ(2Ｖ^3/2/3)/ｄｔ＝2βｕ^1/2{Ｆ_mcosω'ｔ＋(Ｆ/ω')sinω'ｔ}です。

　故に,Ｖ＝[3βｕ^1/2{(Ｆ_m/ω')sinω'ｔ＋(Ｆ/ω'²)(1－cosω'ｔ)}]^2/3＝[(3βｕ^1/2/ｓ'){(ω'Ｆ_msinω'ｔ＋Ｆ(1－cosω'ｔ))}^2/3です。

　一方,(21)ｄＶ/ｄｚ＝2β(ｕＶ)^1/2のｚ＝0 (ｔ＝0)でＶ＝0 の解はＶ＝ｕ(βｚ)²ですから,上式のＶと等置するとｚ＝(Ｖ/ｕ)^1/2/β＝{3/(β²ｕｓ')}^1/3{(ω'Ｆ_msinω'ｔ＋Ｆ(1－cosω'ｔ))}^1/3を得ます。

　三角関数の合成ω'Ｆ_msinω'ｔ－Ｆcosω'ｔ＝(ω'²Ｆ_m²＋Ｆ²)^1/2sin(ω'ｔ＋δ)を用いると,ｚの最大値がｚ_max＝{3/(β²ｕｓ')}^1/3{Ｆ＋(ω'²Ｆ_m²＋Ｆ²)^1/2}^1/3で与えられることがわかります。

　そして,(19)ｓ'＝ｓ(Ｍ/Ｍ_eff)によりｓ'∝ ｓなので,有風時の浮力プルーム(Ｆ_m＝0)の場合は,Δｈ＝ｃ₁{Ｆ/(ｕｓ)}^1/3,ジェットプルーム(Ｆ＝0)の場合は,Δｈ＝ｃ₂(Ｆ_m/ｕ)^1/3ｓ^-1/6と書けます。

　ただし,ｃ₁＝[6/{β²(ｓ'/ｓ)}]^1/3,ｃ₂＝[3/{β²(ｓ'/ｓ)}]^1/3ですから,例えばβ～ 0.5,(ｓ'/ｓ)～ 2.3を代入して計算するとｃ₁～ 2.19,ｃ₂～ 1.73です。

一方,Briggsによれば,有風時の浮力プルームの場合の上昇式は,(30)Δｈ＝2.4{Ｆ/(ｕｓ)}^1/3,有風時のジェットプルームの場合の上昇式は(31)Δｈ＝1.5(Ｆ_m/ｕ)^1/3ｓ^-1/6で与えられるらしいです。

今日はここまでです。この項目はまだ少し続きます。

参考文献：G.A.Briggs "Plume Rise",U.S.Atomic Energy Commision Division of Technical Information(1969)