唐突ですが,過去の資料を整理していたら1983年(33歳当時)の仕事関係の覚書きノートを見つけたので,ブログに覚書きを残す一環としてこれも書いておこうと思います。
内容は,古典的なBriggsによるプルーム(plume)上昇式の紹介です。(←ノートにはプリューム上昇式と書いてありましたが,今はplumeはプルームと読むべきだと思っているので直しました。)
当時は,環境アセス会社の社員で丁度,冷却塔(cooling-tower)からの温排気のシミュレーションについて調べていた頃と思います。
※(この後,最終的にはANLモデル(Argonne国立研究所モデル)と格闘することになったのですが。。この会社をやめて10年以上後にはSACTIというプログラムとも出会いましたね。illinois大学関連かな?)
ただし,当時は常識だった気象力学の知識は今ではかなり忘れているので,この記事内容だけで私自身がわかるように補足しています。
Ⅰ.プルーム上昇式の導出
§1.基礎理論(鉛直プルーム)
プルーム内の局所密度ρpは不変:∂ρp/∂t=0 とします。
すると質量保存の連続の方程式は(1)∇(ρpvp)=0 となります。ただしvp,およびρpはプルームガスの局所速度,および密度です。
また,局所プルーム要素の運動は断熱的とします。
断熱条件はプルームの温位をθpとして(2)dθp/dt=0 です。
温度がT,圧力がpの気体に対して温位はθ≡T(p00/p)(R/Cp)で定義されます。ただしRは気体定数:R≒8.31J/(molK),Cvはモル定積比熱,Cpはモル定圧比熱です。
CvとCpの間にはMayerの関係式:Cp=Cv+Rが成立します。
dθp/dt=0 は中立大気中の運動では浮力が保存すること,移流以外には熱の放射吸収がないことを意味します。
※(補足):ガス(気体)が希薄で理想気体近似ができるとすると,その状態方程式は気体分子の分子量をmとしてp=ρRT/mです。p,ρ,Tは,それぞれ気体の圧力,密度,絶対温度です。
断熱変化の場合には,c1を定数としてPoissonの式:p=c1ργが成立します。γは比熱比でγ≡Cp/Cvです。
p=ρRT/mの両辺の対数微分を取ると,dlnp=dinρ+dlnTですからdlnT=dlnp-dlnρですが,断熱変化ではp=c1ργよりlnp=γlnρ+lnc1が成立するため,dlnρ=dlnp/γです。
それ故,dlnT=(1-1/γ)dlnp=(R/Cp)dlnp,すなわちdln{T/p(R/Cp)}=0 ,またはT/p(R/Cp)=(一定)です。
そこで,ある定数c2を与えθ≡c2T/p(R/Cp)によって温位という量θを定義すると,断熱変化はdθ=0 を意味するので,断熱変化の場合の温位の時間変化率はdθ/dt=0 です。
特に,圧力pの値がp00のときθ=Tとなるように定数c2を定めれば,θ=T(p00/p)(R/Cp)という表現ができます。
さて,プルームの周辺大気について対流が生じない条件はdp/dz=-ρg(静力学平衡)ですが,これは対数表現ではdlnp/dz=-mg/(RT)です。ただし,鉛直上向きをz軸の正の向きとしています。
さらに,断熱変化:dθ=0 ならdlnp=(Cp/R)dlnTなので,熱の放射吸収のない断熱環境で対流が生じない条件はdlnT/dz=-mg/(CpT),またはdT/dz=-mg/Cpと書けます。
そこで,鉛直高さzと共にこの比率で温度が減衰すれば対流が生じず,しかも高さ方向の変化dzに対する温度変化は断熱変化なので,温位勾配はゼロ:dθ/dz=0 です。
この勾配:dT/dz=-mg/Cpの温度減率(気温減率)を断熱減率といい,この状態にある大気を中立状態といいます。
ところで,通常の地上付近成層圏の大気の体積組成(分子数組成)は分子量28の窒素N2が78.08%,分子量32の酸素O2が20.95%,その他はアルゴン,二酸化炭素,水蒸気などです。
そこで,地上付近の乾燥大気の分子量:m=mdは28.964です。(小倉義光「一般気象学」より)
また,N2,O2のような2原子分子では理論的にはモル定積比熱は(5/2)R,モル定圧比熱はCp=(7/2)Rです。
そこで,この大気の断熱減率を計算するとΓd≡mdg/Cp=28.964(g/mol)×9.8(m/s2)×2/{7×8.31J/(molK)}=9.8(K/km)です。
これは,100m上昇する毎に気温が約1度下がるような状態です。この気温減率:Γd=mdg/Cp=9.8(K/km)を乾燥断熱減率といいます。
一般には,大気は湿度ゼロの乾燥空気ではなく水分を含んでいます。
容積Vの中に質量Mdの乾燥空気と質量Mvの水蒸気(vapor)が混在しているとき,この空気の密度はρ=(Md+Mv)/V=ρd+ρvです。
そして,pdを乾燥大気の分圧,eを水蒸気圧とすると,乾燥大気,水蒸気の状態方程式はそれぞれpd=ρdRT/md,e=ρvRT/mvです。
全大気圧pはp=pd+e=ρdRT/md+ρvRT/mvですから,p=ρRT{1+(mv/md)(ρv/ρd)}/{1+(ρv/ρd)},またはp=ρRT*/md;T*≡T{1+(mv/md)(ρv/ρd)}/{1+(ρv/ρd)}と書けます。
ここで,T*は仮温度(virtual temperature)と呼ばれる量で具体的にmd=28.964(g/mol),mv=18(g/mol)を代入すると,水蒸気の混合比:q≡ρv/ρd=(md/mv){e/(p-e)}に対してT*≒T(1+1.61q)/(1+q)=T(1+0.61r)と表現されます。
ただし,r≡q/(1+q)=ρv/(ρv+ρd)は比湿(specific humidity)と呼ばれる量で,式から明らかにr≦1です。
乾燥大気の温度Tに対する安定条件:dT/dz=-mdg/Cp=-Γdは,Tを仮温度T*に置き換えると近似的にそのまま成立すると考えれば,dT*/dz=-mdg/Cp=-Γdなる式が安定条件に対応します。
左辺のT*にT*=T(1+0.61r)を代入すると,(1+0.61r)(dT/dz)=-Γdです。そこで,不飽和で凝結のない大気の断熱減率ΓはΓ=Γd/(1+0.61r)で与えられ,Γ≦Γdです。
さらに,飽和(saturate)して凝結水滴があり潜熱(latent heat)の存在によって温度降下が鈍くなる場合の湿潤断熱減率をΓwとします。
飽和水蒸気圧をesと書き,単位質量当たりの液体の水が水蒸気になるのに必要な潜熱(蒸発熱)をLとすれば,esと温度Tの関係は理論的にはClausius-Clapeyronの式:des/dT=Les/(RT2)で表わされます。
(上式は,lnes=-L/(RT)+const.またはes=Aexp{-L/(RT)}を意味しますが,2006年7/30の記事「気液平衡の統計力学」ではこれより精密な修正式:lnes=lnT-L/(RT)-1+const.またはes=ATexp{-L/(RT)}を与えています。)
例えばT=100℃ではL=597(cal/g)ですが常温付近でのLはTの変動に鈍感なのでLがTに無関係と近似してこれを積分すれば,ln(es/es0)=(mvL/R)(1/T0-1/T)です。
これはesがTだけに依存して周囲の大気圧pには依らない形の関係式です。
さて,乾燥空気と混合比qの水蒸気が混合した空気塊を考えます。
まず,単一気体の断熱変化は,0=dlnθ=(Cp/R)dlnT-dlnpで表わされます。これはエネルギー収支を示す式と見るなら,0=RTdθ/θ=CpdT-RTdp/pと表現されます。
そこで,乾燥空気と水蒸気が混合しただけの大気において水蒸気が飽和に達していないとき,各成分の上記左辺の温位変化:RTdθ/θへの寄与は,それぞれCpdT-RT(dpd/pd)=0,qCpvdT-qRT(de/e)=0 です。
ここで,qはこの大気の乾燥空気に対する水蒸気の混合比でq=ρv/ρd=(md/mv)(e/pd)=(md/mv){e/(p-e)}です。また,Cpvは水蒸気の定圧モル比熱です。このとき,もちろんpd=p-eです。
この空気塊が上昇してq=qsとなって飽和した後,さらに上昇してdqsの水蒸気が凝結したとすれば,空気塊はLdqsだけの潜熱を受けるので,変化が断熱なら水蒸気が飽和した空気塊のエネルギー収支の式は次のようになります。
すなわち,CpdT+qsCpvdT-RT{d(p-es)/(p-es)+qsdes/es}+Ldqs=0 となります。qsは飽和時の混合比です。
これに,des/dT=Les/(RT2)を代入すれば,-Ld(qs/T)=CpdT/T+qsCpvdT/T-Rd(p-es)/(p-es)を得ます。
pd=ρdRT/md,es=ρvRT/mvですから,静力学平衡:dp/dz=-ρg, or d(pd+es)/dz=-(ρd+ρV)gは,d(pd+es)/dz=-g(mdpd+mves)/(RT),またはdp=-g{md(p-es)+mVes}/(RT)dzと表現できます。
そこで,d(p-es)=-g{md(p-es)+mVes}/(RT)dz-desを-Ld(qs/T)=(Cp+qsCpv)dT/T-Rd(p-es)/(p-es)に代入すると,-Ld(qs/T)=(Cp+qsCpv)dT/T+g{md+mVes/(p-es)}dz/T+Rdes/(p-es)です。
これから,dT/dz=-g{md+mVes(p-es)}/(Cp+qsCpv)-RT(des/dz)/{(p-es)(Cp+qsCpv)}+{LT/(Cp+qsCpv)}{d(qs/T)/dz}を得ます。
ところが,qs=(md/mv){es/(p-es)}よりes/(p-es)=mvqs/mdです。
そこで,dqs/dz=(md/mv)[(des/dz)/(p-es)-{d(p-es)/dz}{es/(p-es)2}]=(md/mv)[(des/dz)/(p-es)-{qs/(p-es)}{d(p-es)/dz]と書けます。
つまり,(des/dz)/(p-es)=(mv/md)(dqs/dz)+{qs/(p-es)}{d(p-es)/dz}です。これらの式を用いて,できるだけesを消去します。
以上から,湿潤断熱減率は,Γw≡-dT/dz=Γd/(1+qsCpv/Cp)+mV2g/{mdCp(1+qsCpv/Cp)}+mvRT(dqs/dz)/{mdCp(1+qsCpv/Cp}+RT{qs/(p-es)}{d(p-es)/dz}/{Cp(1+qsCpv/Cp-LT{d(qs/T)/dz}/{Cp(1+qsCpv/Cp)}です。
この式は閉じていないので古い正野重方著の岩波全書の「気象力学」などには近似公式もあるようですがここでは主要課題ではないので割愛します。ただ,一般に湿潤断熱減率は一定値ではありません。
(なお,水蒸気の定圧比熱については,2008年1/6の記事「氷,水,水蒸気の比熱」,および2009年2/9の記事「水蒸気の比熱」などを参照してください。) (補足終わり) ※
さて,浮力以外の気圧傾度力は無視できると考えて外気圧pは近似的に静力学平衡:dp/dz=-ρgを満足しているとします。
また,大気を流体と見るとき,対象とする運動でのReynolds数:Reは非常に大きいため,1/Reに比例する分子粘性項も無視します。
そして,プルームガスの周辺大気(ambient air)の局所速度,圧力,密度,温度,温位を,それぞれva,pa,ρa,Ta,θaとします。
特に,気圧については,∂pa/∂x=∂pa/∂y=0,∂pa/∂z=-ρagが成立しており,プルームガス圧については近似的にpp=paが成立するとします。
この場合,周辺大気の運動方程式は,ρadva/dt=-∇pa-ρagk=0 です。ただしi,j,kは,それぞれx,y,z軸の方向単位ベクトルを表わします。
この方程式の解は明らかにva=constで与えられます。
一方,プルームガスの満たす運動方程式はρpdvp/dt=-∇pp-ρpgk~(-∇pa-ρagk)-(ρp-ρa)gkです。
すなわち,この近似での運動方程式はρp(dvp/dt)=(ρa-ρp)/gk,またはdvp/dt={(ρa-ρp)/ρp}gkです。これらの右辺が浮力項(buoyancy term)です。
状態方程式:pa=ρaRTa/M,pp=ρpRTp/Mから,pp~paの近似ではρa/ρp ~Tp/Ta~θp/θaなので温位偏差θ'をθ'≡θp-θaと定義すれば,(ρa-ρp)/ρp~ (θp-θa)/θa=θ'/θaです。
周辺大気の温位θaの定義:θa≡Ta(p00/pa)(R/Cp)において特にp00~paなるパラメータを採用すればθa~Taなので,結局この近似(Boussinesq近似)での運動方程式は,(3)dvp/dt=(gθ'/Ta)kとなります。
連続方程式が質量の保存を意味するのに対して,運動方程式は運動量の保存を意味します。
さらに,vp≡va+v'と書き,特に速度偏差:v'の鉛直z方向の成分(or第3成分)をw'とします。また,周辺大気の風速vaは鉛直成分を持たないとして,その向きをx軸(or1軸)の正の向きに取ります。
すると,vp≡ui,vp3=v'3=w'と書けます。さらにθa,ρa,uはzだけの関数でありθa=θa(z),ρa=ρa(z),u=u(z)のように表わせるとします。
そして水平プルーム断面上ではw'と浮力gθ'/Taが比例する:gθ'/Ta∝ w'と仮定します。比例係数をc(z)とするとgθ'/Ta=c(z)w'です。
定常状態が実現されていると仮定し,面積分することによってプルームを水平面で切った断面上での平均量に対する方程式を求めます。
ただし,プルーム断面の外では明らかにθ'=0,かつv'=0 で特にw'=0 です。
この結果は,まず浮力の効果についての方程式:(4)dFz/dz=-sVを与えます。ただし,Fzは浮力流束で,(5)Fz≡∫P(gρpθ'w'/Ta)dxdy/(πρa)で定義されます。
Vは鉛直容積流束で,(6)V≡∫Pρpw'dxdy/(πρa)です。
また,sは安定度パラメータと呼ばれる量で,(7)s≡(g/Ta)(∂θa/∂z)です。
(証明):(2)dθp/dt=0 より,∂θp/∂t+vp∇θp=0ですから,定常状態:∂θp/∂t=0 ではρpvp∇(θa+θ')= 0です。
そして,∂θa/∂x=∂θa/∂y=∂θ'/∂x=∂θ'/∂y=0 よりρpvp∇(θa+θ')=-ρpw'(∂θa/∂z)です。
これに,(1)∇(ρpvp)=0 にθ'を掛けたものを加えると∇(ρpθ'vp)=-ρpw'(∂θa/∂z)を得ます。
これを全水平面で積分すると∫∞∇(ρpθ'vp)dxdy=-(∂θa/∂z)∫∞ρpw'dxdyですが,無限遠では明らかにθ'=0 で,プルーム断面を含む全水平面領域では全ての量はzに依らず一定ですから,左辺の(∂/∂z)も積分記号の外に出せます。
プルーム水平断面の外ではρpθ'w'=ρpw'=0 なので,上式は(d/dz)∫Pρpθ'w'dxdy=-(∂θa/∂z)∫Pρpw'dxdyとなります。
周辺大気の密度ρaは"上昇高度までの層当たり平均密度=(一定)"と考え,両辺に(g/Ta)/(πρa)を掛けると{∂(1/Ta)}{(d/dz)∫Pρpgθ'w'dxdy}/(πρa)=-(g/Ta)(∂θ/∂z){∫∞ρpw'dxdy}/(πρa)です。
よって,dFz /dz=-sV-(1/Ta)(dTa/dz)Fzを得ました。
ところが,dTa/dz~dθa/dz,Fz ~gV(<θ'>/θa)なので(1/Ta)(dTa/dz)Fz~ (sV)・(<θ'>/θa)<<(sV)ですから,右辺第1項:(sV)に比べて右辺第2項:(1/Ta)(dTa/dz)Fzを無視します。
そこで,結局dFz /dz=-sVを得ます。(証明終わり)
次に,運動量に対する方程式は,(8)d(vV)/dz=(Fz/w)k+va(dV/dz)+(抵抗力)です。ただし,vVは(9)vV≡∫Pvpρpw'dxdy/(πρa)で定義される運動量流束です。
この式で定義されるベクトルvをプルームの速度と考えます。wはvのz成分です。
(証明):(1)∇(ρpvp)=0 より,vp・∇(ρpvp)=0 ,つまりΣj=13vpi∂j(ρpvpj)=0 です。
また,(3)dvp/dt=(gθ'/Ta)kより,定常状態ではρp(vp∇)vp=(gρpθ'/Ta)k,つまりΣj=13ρpvpj∂jvpi=(gρpθ'/Ta)δi3です。
したがって,∇(ρpvpvpi)=Σj=13vpi∂j(ρpvpj)+Σj=13ρpvpj∂jvpi=(gρpθ'/Ta)δi3です。
プルームの進行方向をnとするとき,プルームの断面と側面で囲まれた領域VにGaussの定理を適用すると,∫V∇(ρpvpvpi)dV=∫P+ρpvpnvpidS-∫P-ρpvpnvpidS-Δn{1+(dr/dn)2}1/2(∫Cρpvpvpids)です。

ただし,ここでのrはプルーム半径(plume-radius)を意味します。プルーム断面は等方的なので,これは断面の形を円と想定したときの半径です。
-Δndsはプルーム側面の外向き法線の向きを持ちΔnは側面n方向の厚みです。それ故,dsは周の微小線素の長さを持つ内向き法線方向ベクトルです。
故に,Δ∫PρpvpnvpidS≡∫P+ρpvpnvpidS-∫P-ρpvpnvpidS=Δn{1+(dr/dn)2}1/2(∫Cρpvpvpids)+δi3Δn∫P(gρpθ'/Ta)dSです。
よって,(d/dn)(∫PρpvpnvpidS)={1+(dr/dn)2}1/2(∫Cρpvpvpids)+δi3∫P(gρpθ'/Ta)dSを得ます。
一方,連続の方程式から(d/dn)(∫PρpvpndS)={1+(dr/dn)2}1/2(∫Cρpvpdsです。
さらにvpi=vai+v'iですから,(d/dn)(∫PρpvpnvpidS)=δi3∫P(gρpθ'/Ta)dS+vai(d/dn)(∫PρpvpndS)+{1+(dr/dn)2}1/2(∫Cρpvpv'ids)を得ます。
特にnをz方向のkに選べばvpn=w'より,(d/dz)(∫Pρpw'vpidS)=δi3∫P(gρpθ'/Ta)dS+vi(d/dz)(∫Pρpw'dS)+(抵抗力)です。
ここで,(抵抗力)≡{1+(dr/dz)2}1/2(∫Cρpvpv'ids)であり,ds={1+(dr/dz)2}-1/2(-dyi+dxj)+{1+(dr/dz)2}-1/2(dr/dz)dlkです。(dl=(dx2+dy2)1/2)
一方,∫P(gρpθ'/Ta)dS=(1/w)∫P(g/Ta)(ρpθ'w)dS=(1/w)∫P(g/Ta)(ρpθ'w')dS+(1/w)∫P(g/Ta){ρpθ'(w-w')}dSです。
vV≡∫Pvpρpw'dS/(πρa),V≡∫Pρpw'dS/(πρa),およびvp3=w'よりw=∫Pρpw'2dS/∫Pρpw'dSですから,w-w'={∫Pρpw'2dS-w'∫Pρpw'dS}/{∫Pρpw'dS}です。
仮定によりgθ'/Ta=c(z)w'を代入すると,(1/w)∫P(g/Ta){ρpθ'(w-w')}dS={c(z)/w}[{∫Pρpw'dS}{∫Pρpw'2dS}-{∫Pρpw'2dS}{∫Pρpw'dS}]/{∫Pρpw'dS}=0 です。
したがって,∫P(gρpθ’/Ta)dS=(1/w)∫P(g/Ta)(ρpθ'w')dS故,(d/dz)(∫Pρpw'vpidS)=(1/w)δi3∫P(g/Ta)(ρpθ'w')dS+vi(d/dz)(∫Pρpw'dS)+(抵抗力),
または,(d/dz)(∫Pvpρpw'dS)=(1/w){∫P(g/Ta)(ρpθ'w')dS}k+va(d/dz)(∫Pρpw'dS)+(抵抗力)なる式を得ます。
最後の表式は,まさにd(vV)/dz=(Fz/w)k+va(dV/dz)+(抵抗力)そのものです。
※特にプルームが軸対称で完全に鉛直向きなら(抵抗力)={1+(dr/dz)2}1/2(∫Cρpvpv'ids)=0でd(wV)/dz=Fz/wです。※(証明終わり)
さて,dFz/dz=-sV,d(wV)/dz=Fz/wを解くために次のような初期条件を与えます。
まず,初期速度は吐出速度に等しいとします。すなわち,(10a)v=w0kです。そして初期運動量流束は,(10b)vV=(wV)k=Fmk≡(ρ0/ρa)w02r02kとします。ρ0は排出ガス密度,r0は煙突口の口径です。
さらに,初期の浮力流束は,(10c)Fz=F≡(1-ρ0/ρa)gw0r02で与えられるとします。
つまり,初期容積流束をV0≡w0r02とすれば初期の総質量流束はπρ0V0=πρ0w0r02ですから,定義により運動量流束はFm=πρ0w0V0/(πρa)=(ρ0/ρa)w0V0です。
そして,初期浮力流束はFz=F=(1-ρ0/ρa)gV0~{(T0-Ta)/T0}gV0となるわけです。
ところで,排出源がhot-sourceの場合の初期浮力流束は,(11)F=gQH/(πρacpTa)で与えられます。
ここで,QHは排出されるガスによって単位時間に排出される熱量(cal/sec)であり,cpは空気の単位質量当たりの熱容量(定圧比熱)(cal/(kg・K))です。
なぜなら,Fz=F=(gθ'/Ta)(πρ0w0r02)/(πρa)={g/(πρaTa)}(θ'πρ0V0)ですが,πρ0V0は総質量流束なので排出ガスによって単位時間に外気に付加される熱量はQH=cpT'πρ0V0 ~cpθ'πρ0V0に等しいからです。
つまり,θ'πρ0V0=QH/cpによってF=gQH/(πρacpTa)が得られるわけです。
この場合,排熱強度QHと排出ガス密度ρ0の情報を得ることにより逆にV0=w0r02が求まります。
そこで,煙突口径r0のデータから吐出速度w0が算定されて,結局全ての初期情報を得ることができます。
ただし,後述するように3つの未知数Fz,w,Vに対し微分方程式は2つしかないので解を決めるにはもう1つ方程式が必要です。
今日はここまでですがまだ続きます。
参考文献:G.A.Briggs "Plume Rise",U.S.Atomic Energy Commision Division of Technical Information(1969),小倉義光「一般気象学」(東京大学出版会),小倉義光「気象力学通論」(東京大学出版会)
最近のコメント