« ソルジャー・ブルー祈念日 | トップページ | 眼科に行きました。 »

2010年12月 1日 (水)

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

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

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

 これまでの議論では,未知数Fz,w,Vに対して2つの方程式しか与えられておらず,問題が閉じていません。

 

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

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

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

 そして,プルームの行路長さ:sに対してds≡vsdt,またはdt≡ds/vsによってdtを定義しておきます。

 すると,質量の保存式から明らかにΔ(ρpπr2s)=ρa2πrveΔs,または(d/ds)(ρp2sa)=2rveです。こうした関係があることをエントレインメントの仮定といいます。

 特に,鉛直プルーム:s=z,dz=wdtの場合にはV=ρpπr2w/(πρa)=ρp2w/ρaなのでエントレインメントの仮定はdV/dz=2rveを意味します。

 また,(d/ds)(ρp2sa)=2rveよりρaはsに依らないとして(r2a){d(ρps)/ds}+(2rρpsa)(dr/ds)=2rveです。

さらに,dt=ds/vsよりd/ds=(1/vs)(d/dt)を用いればdr/dt=ρaep-{r/(2ρp)}{d(ρps)/ds}を得ます。

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

また,再びV=ρpπr2w/(πρa)=ρp2w/ρaよりr2=ρaV/(ρpw)~V/wと近似できるという考察から,逆にプルーム半径:rをr≡(V/w)1/2によって定義します。

この定義によれば,V=r2wですが,Δ(ρpπr2s)=ρa2πrveΔs,または(d/ds)(ρp2sa)=2rveの関係は変わらず,ρp~ρaの近似で考えると依然としてエントレインメントの仮定はdV/dz=2rveです。

元々,これも仮定ですから以下ではdV/dz=2rveをエントレインメントの仮定として採用します。この右辺にve=αw,r=(V/w)1/2を代入するとdV/dz=2α(wV)1/2を得ます。

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

(4)dFz/dz=-sV,(12)d(wV)/dz=Fz/w,

および,(13)dV/dz=2α(wV)1/2 です。

一応,未知量V,Fz,wの意味を再び述べておきます。

 

Vは鉛直容積流束,Fzは鉛直方向の浮力(gθ'/T)の流束,wはプルーム面の鉛直速度成分です。

  

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

 

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

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

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

折れ曲がりプルームの場合,積分面をプルームと交わる鉛直面として,この面上ではs=(g/T)(∂θa/∂z),およびa=uは一定であると仮定します。

前と同様,∫pθ'p)dydz=-(∂θa/∂z)∫ρpw'dydzですがプルームの外ではθ'=0,w'=0 とすれば,これは(d/dx)∫Pρpθ'(u+u')dydz=-(∂θa/∂z)∫Pρpw'dydzを意味します。

そこで,浮力流束を(5)Fz≡∫P(gρpθ'w'/Ta)dxdy/(πρa)の代わりに(14)Fz≡∫P{gρpθ'(u+u')/Ta}dydz/(πρa)で定義すれば,(15)dFz/dx=-sMと表わすことができます。

ただし,再びsはs=(g/Ta)(∂θa/∂z)でありMはM≡∫Pρpw'dydz/(πρa)で定義される量です。

この場合も排出源がhot-sourceなら,浮力流束Fzの初期値Fは{(単位時間にプルーム積分面を通過する総排出熱量)×(g/Ta)/cp}/(πρa)と考えられるので,やはりF=gQH/(πρapa)で与えられます。

 そして,運動量に対する方程式は,(8)d(V)/dz=(Fz/w)a(dV/dz)+(抵抗力)の証明の際に遭遇した式:(d/dn)(∫PρppnpidS)={1+(dr/dn)2}1/2(∫Cρpppi)+δi3P(gρpθ'/Ta)dSにおいて進行方向に置くと得られます。

 d/dnをd/dxに変え,p=(u+u')i+v'+w',d={1+(dr/dx)2}-1/2(-dz+dy)+{1+(dr/dx)2}-1/2(dr/dx)dlとして,必要なi=3の場合のみを書き下します。

 

 それは,(d/dx){∫Pρp(u+u')w'dS}=∫P(gρpθ'/Ta)dS+{1+(dr/dx)2}1/2(∫Cρppw'd)=∫P(gρpθ'/Ta)dS-Cρpv'w'dz+∫Cρpw'2dy+∫Cρp(u+u')w'(dr/dx)dlです。

 u'<<uなのでu+u'~ uと近似して有効運動量流束をuMeff;Meff≡{∫Pρpw'dS-∫dx∫Cρpw'(dr/dx)dl}/(πρa)と定義しu',v',w'の2次以上の項を無視します。

 

 すると,Fz={u∫P(gρpθ'/Ta)dS}/(πρa)より,結局運動方程式として(16)u(dMeff/dx)=Fz/uが得られます。

 鉛直運動量の水平方向流束uMeffの初期値はプルーム面を通過する鉛直運動量ρpwπr2/(πρa)=ρpwr2aの総流束(鉛直方向流束)ですから,その初期値はFm=(ρ0002a)×w0=ρ00202a=ρ000aです。

そして,Briggs(1972),Richards(1963)によれば有効運動量流束uMeffの運動量流束uM;M=∫Pρpw'dydz/(πρa)に対する比率の評価値はMeff/M ~ 2.3です。

zをプルーム軸心の高さ,dt≡dx/uとしてw≡dz/dt=udz/dxと定義します。wはプルーム軸心の上昇速度を与え,解(x,z),or z=z(x)はプルームの軌跡を示します。

ここで,折れ曲がりプルームにおけるエントレインメント係数をβとします。すなわちve=βwです。有効プルーム半径をrとすると,Meff=πr2ρpw/(πρa)=wρp2aです。

折れ曲がりプルームの場合のエントレインメントの仮定はd(ρpπr2u)/dx=2πrρae,つまりd(ρp2u/ρa)/dx=2rveと表わせます。

そこで,ud(Meff/w)/dx=2rve=2rβwです。

 

r≡(Meff/w)1/2と定義すれば,(17)ud(Meff/w)/dx=2β(wMeff)1/2を得ます。

さらにwV≡uMeffによってVを定義すれば,Vは水平方向の容積流束でありr=(V/u)1/2です。これとu/dx=w/dzなる関係を用いて基本方程式系:(15),(16),(17)を書き直せば次のようになります。

すなわち,(18)dFz/dz=-s'V,(19)s'=s(M/Meff),(20)d(wV)/dz=Fz/w,(21)dV/dz=2β(uV)1/2です。

(21):dV/dz=2β(uV)1/2は,z=0 でV=0 という点源(point source)の初期条件に対しては容易に解けて,(V/u)1/2=βz,つまりr=βzです。

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

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

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

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

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

 さて,まずdx=udt,dz=wdtにより,全ての方程式を独立変数tの式に変換します。

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

 

(22)dFz/dt=-s(wV),(23)d(wV)/dt=Fz,(24)dV/dt=2α(wV)3/2/Vです。

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

 

(25)dFz/dt=-s'(wV),(26)d(wV)/dt=Fz,(27)dV/dt=2βwV(u/V)1/2です。

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

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

 このケースの基本方程式系は(22)dFz/dt=-s(wV),(23)d(wV)/dt=Fz,(24)dV/dt=2α(wV)3/2/Vです。

(22),(23)よりd2(wV)/dt2=-s(wV)を得ます。これはよく知られている線形振動の方程式です。

 

そこで,ω≡s1/2としてこれをd2(wV)/dt2=-ω2(wV)と書き,初期条件:t=0 でwV=Fm,d(wV)/dt=Fz=Fの解を求めると,wV=Fmcosωt+(F/ω)sinωtを得ます。

鉛直運動量流束wVがゼロとなる最初のt=t1がプルーム上昇高さを与えると考えられます。

  

上昇高さはΔh=∫0t1wdtですから,仮定Ⅱの浮力プルーム(Fm=0)ではt1=π/ω,よりΔh=ω-10πwdλです。

 

同様にジェットプルーム(F=0)ではt1=π/(2ω)よりΔh=ω-10π/2wdλ(λ=ωt)です。

     浮力プルーム(Fm=0)の場合

 (V2/2)/dt=V(dV/dt)=2α(wV)3/2の右辺にwV=(F/ω)sinωtを代入してd(V2/2)/dt=2α(F/ω)3/2sin3/2ωtです。

故に,V2/2=2α(F/ω)3/2ω-1{∫ωtsin3/2λdλ+const.},すなわちV=2α1/23/4ω-5/4{∫ωtsin3/2νdν+const.}1/2です。

 w=wV/V=(F/ω)sinωt/Vより,Δh=ω-10πwdλ=(1/2)α-1/21/4ω-3/4[∫0πsinλdλ{∫0λsin3/2νdν+const.}-1/2]を得ます。

したがって,ω=s1/2より定数cをc≡(1/2)α-1/2[∫0πsinλdλ{∫0λsin3/2νdν+const.}-1/2]で定義すれば,Δh=cF1/4-3/8なる関係式が成立することになります。

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

     ジェットプルーム(F=0)

 wV=Fmcosωtよりd(V2/2)/dt=V(dV/dt)=2α(wV)3/2の右辺にwV=(F/ω)sinωtを代入すれば,d(V2/2)/dt=2αFm3/2cos3/2ωtです。

故に,V2/2=2αFm3/2ω-1{∫ωtcos3/2λdλ+const.},すなわちV=2α1/2m3/4ω-1/2{∫ωtcos3/2νdν+const.}1/2です。

 そして,w=wV/V=Fmcosωt/Vより,Δh=ω-10π/2wdλ=(1/2)α-1/2m1/4ω-1・2[∫0π/2cosλdλ{∫0λcos3/2νdν+const.}-1/2]です。

したがって,ω=s1/2より定数c0をc0≡(1/2)α-1/2[∫0π/2cosλdλ{∫0λcos3/2νdν+const.}-1/2]と定義すれば,Δh=c0m1/4-1/4なる関係の成立することがわかります。

 t=0 のとき,V=0 の点源(const.=0)では,μ=π/2-νと変数変換すれば,cosν=sinμより∫0λcos3/2νdν=-∫π/2π/2-λsin3/2μdμ,c0=(1/2)α-1/2[∫0π/2cosλdλ{∫π/2-λπ/2sin3/2μdμ}-1/2]=(1/2)α-1/2[∫0π/2sinσdσ{∫0σsin3/2μdμ}-1/2]となります。

 点源のとき,c=(1/2)α-1/2[∫0πsinλdλ{∫0λsin3/2νdν+const.}-1/2]に対してc0=(1/2)α-1/2[∫0π/2sinλdλ{∫0λsin3/2νdν+const.}-1/2]であり,ほとんどcと同じです。

cを求めたのと同じ数値計算によれば,α=0.132に対してc0~ 4.0ですから,無風状態でのジェットプルームの評価式として(29)Δh=4.0(Fm/s)1/4を得ます。

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

 このケースの基本方程式系は(25)dFz/dt=-s'(wV),(26)d(wV)/dt=Fz,(27)dV/dt=2βwV(u/V)1/2です。

(25),(26)よりd2(wV)/dt2=-s(wV)なのでω'=s'1/2としてwV=Fmcosω't+(F/ω')sinω'tです。

 

これを(27)から得られるd(2V3/2/3)/dt=2βu1/2wVの右辺に代入すると,d(2V3/2/3)/dt=2βu1/2{Fmcosω't+(F/ω')sinω't}です。

 故に,V=[3βu1/2{(Fm/ω')sinω't+(F/ω'2)(1-cosω't)}]2/3=[(3βu1/2/s'){(ω'Fmsinω't+F(1-cosω't))}2/3です。

 一方,(21)dV/dz=2β(uV)1/2のz=0 (t=0)でV=0 の解はV=u(βz)2ですから,上式のVと等置するとz=(V/u)1/2/β={3/(β2us')}1/3{(ω'Fmsinω't+F(1-cosω't))}1/3を得ます。

 三角関数の合成ω'Fmsinω't-Fcosω't=(ω'2m2+F2)1/2sin(ω't+δ)を用いると,zの最大値がzmax={3/(β2us')}1/3{F+(ω'2m2+F2)1/2}1/3で与えられることがわかります。

 そして,(19)s'=s(M/Meff)によりs'∝ sなので,有風時の浮力プルーム(Fm=0)の場合は,Δh=c1{F/(us)}1/3,ジェットプルーム(F=0)の場合は,Δh=c2(Fm/u)1/3-1/6と書けます。

 ただし,c1=[6/{β2(s'/s)}]1/3,c2=[3/{β2(s'/s)}]1/3ですから,例えばβ~ 0.5,(s'/s)~ 2.3を代入して計算するとc1~ 2.19,c2~ 1.73です。

一方,Briggsによれば,有風時の浮力プルームの場合の上昇式は,(30)Δh=2.4{F/(us)}1/3,有風時のジェットプルームの場合の上昇式は(31)Δh=1.5(Fm/u)1/3-1/6で与えられるらしいです。

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

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

|

« ソルジャー・ブルー祈念日 | トップページ | 眼科に行きました。 »

202. 気象・地学・環境」カテゴリの記事

コメント

コメントを書く



(ウェブ上には掲載しません)




トラックバック

この記事のトラックバックURL:
http://app.f.cocolog-nifty.com/t/trackback/71281/37921881

この記事へのトラックバック一覧です: プルーム上昇のモデル式(2):

« ソルジャー・ブルー祈念日 | トップページ | 眼科に行きました。 »