唐突ですが,過去の資料を整理していたら1983年(33歳当時)の仕事関係の覚書きノートを見つけたので,ブログに覚書きを残す一環としてこれも書いておこうと思います。
内容は,古典的なBriggsによるプルーム(plume)上昇式の紹介です。(←ノートにはプリューム上昇式と書いてありましたが,今はplumeはプルームと読むべきだと思っているので直しました。)
当時は,環境アセス会社の社員で丁度,冷却塔(cooling-tower)からの温排気のシミュレーションについて調べていた頃と思います。
※(この後,最終的にはANLモデル(Argonne国立研究所モデル)と格闘することになったのですが。。この会社をやめて10年以上後にはSACTIというプログラムとも出会いましたね。illinois大学関連かな？)
ただし,当時は常識だった気象力学の知識は今ではかなり忘れているので,この記事内容だけで私自身がわかるように補足しています。
Ⅰ.プルーム上昇式の導出
§1.基礎理論(鉛直プルーム)
プルーム内の局所密度ρpは不変:∂ρp/∂ｔ＝0 とします。
すると質量保存の連続の方程式は(1)∇(ρpｖp)＝0 となります。ただしｖp,およびρpはプルームガスの局所速度,および密度です。
また,局所プルーム要素の運動は断熱的とします。
断熱条件はプルームの温位をθpとして(2)ｄθp/ｄｔ＝0 です。
温度がＴ,圧力がｐの気体に対して温位はθ≡Ｔ(ｐ00/ｐ)(Ｒ/Ｃp)で定義されます。ただしＲは気体定数:Ｒ≒8.31Ｊ/(molＫ),Ｃvはモル定積比熱,Ｃpはモル定圧比熱です。
ＣvとＣpの間にはMayerの関係式:Ｃp＝Ｃv＋Ｒが成立します。
ｄθp/ｄｔ＝0 は中立大気中の運動では浮力が保存すること,移流以外には熱の放射吸収がないことを意味します。
※(補足):ガス(気体)が希薄で理想気体近似ができるとすると,その状態方程式は気体分子の分子量をｍとしてｐ＝ρＲＴ/ｍです。ｐ,ρ,Ｔは,それぞれ気体の圧力,密度,絶対温度です。
断熱変化の場合には,ｃ1を定数としてPoissonの式:ｐ＝ｃ1ργが成立します。γは比熱比でγ≡Ｃp/Ｃvです。
ｐ＝ρＲＴ/ｍの両辺の対数微分を取ると,ｄlnｐ＝ｄinρ＋ｄlnＴですからｄlnＴ＝ｄlnｐ－ｄlnρですが,断熱変化ではｐ＝ｃ1ργよりlnｐ＝γlnρ＋lnｃ1が成立するため,ｄlnρ＝ｄlnｐ/γです。
それ故,ｄlnＴ＝(1－1/γ)ｄlnｐ＝(Ｒ/Ｃp)ｄlnｐ,すなわちｄln{Ｔ/ｐ(Ｒ/Ｃp)}＝0 ,またはＴ/ｐ(Ｒ/Ｃp)＝(一定)です。
そこで,ある定数ｃ2を与えθ≡ｃ2Ｔ/ｐ(Ｒ/Ｃp)によって温位という量θを定義すると,断熱変化はｄθ＝0 を意味するので,断熱変化の場合の温位の時間変化率はｄθ/ｄｔ＝0 です。
特に,圧力ｐの値がｐ00のときθ＝Ｔとなるように定数ｃ2を定めれば,θ＝Ｔ(ｐ00/ｐ)(Ｒ/Ｃp)という表現ができます。
さて,プルームの周辺大気について対流が生じない条件はｄｐ/ｄｚ＝－ρｇ(静力学平衡)ですが,これは対数表現ではｄlnｐ/ｄｚ＝－ｍｇ/(ＲＴ)です。ただし,鉛直上向きをｚ軸の正の向きとしています。
さらに,断熱変化:ｄθ＝0 ならｄlnｐ＝(Ｃp/Ｒ)ｄlnＴなので,熱の放射吸収のない断熱環境で対流が生じない条件はｄlnＴ/ｄｚ＝－ｍｇ/(ＣpＴ),またはｄＴ/ｄｚ＝－ｍｇ/Ｃpと書けます。
そこで,鉛直高さｚと共にこの比率で温度が減衰すれば対流が生じず,しかも高さ方向の変化ｄｚに対する温度変化は断熱変化なので,温位勾配はゼロ:ｄθ/ｄｚ＝0 です。
この勾配:ｄＴ/ｄｚ＝－ｍｇ/Ｃpの温度減率(気温減率)を断熱減率といい,この状態にある大気を中立状態といいます。
ところで,通常の地上付近成層圏の大気の体積組成(分子数組成)は分子量28の窒素Ｎ2が78.08％,分子量32の酸素Ｏ2が20.95％,その他はアルゴン,二酸化炭素,水蒸気などです。
そこで,地上付近の乾燥大気の分子量:ｍ＝ｍdは28.964です。(小倉義光「一般気象学」より)
また,Ｎ2,Ｏ2のような２原子分子では理論的にはモル定積比熱は(5/2)Ｒ,モル定圧比熱はＣp＝(7/2)Ｒです。
そこで,この大気の断熱減率を計算するとΓd≡ｍdｇ/Ｃp＝28.964(g/mol)×9.8(ｍ/ｓ2)×2/{7×8.31Ｊ/(molＫ)}＝9.8(Ｋ/ｋｍ)です。
これは,100ｍ上昇する毎に気温が約１度下がるような状態です。この気温減率:Γd＝ｍdｇ/Ｃp＝9.8(Ｋ/ｋｍ)を乾燥断熱減率といいます。
一般には,大気は湿度ゼロの乾燥空気ではなく水分を含んでいます。
容積Ｖの中に質量Ｍdの乾燥空気と質量Ｍvの水蒸気(vapor)が混在しているとき,この空気の密度はρ＝(Ｍd＋Ｍv)/Ｖ＝ρd＋ρvです。
そして,ｐdを乾燥大気の分圧,ｅを水蒸気圧とすると,乾燥大気,水蒸気の状態方程式はそれぞれｐd＝ρdＲＴ/ｍd,ｅ＝ρvＲＴ/ｍvです。
全大気圧ｐはｐ＝ｐd＋ｅ＝ρdＲＴ/ｍd＋ρvＲＴ/ｍvですから,ｐ＝ρＲＴ{1＋(ｍv/ｍd)(ρv/ρd)}/{1＋(ρv/ρd)},またはｐ＝ρＲＴ*/ｍd;Ｔ*≡Ｔ{1＋(ｍv/ｍd)(ρv/ρd)}/{1＋(ρv/ρd)}と書けます。
ここで,Ｔ*は仮温度(virtual temperature)と呼ばれる量で具体的にｍd＝28.964(ｇ/mol),ｍv＝18(ｇ/mol)を代入すると,水蒸気の混合比:ｑ≡ρv/ρd＝(ｍd/ｍv){ｅ/(ｐ－ｅ)}に対してＴ*≒Ｔ(1＋1.61ｑ)/(1＋ｑ)＝Ｔ(1＋0.61ｒ)と表現されます。
ただし,ｒ≡ｑ/(1＋ｑ)＝ρv/(ρv＋ρd)は比湿(specific humidity)と呼ばれる量で,式から明らかにｒ≦1です。
乾燥大気の温度Ｔに対する安定条件:ｄＴ/ｄｚ＝－ｍdｇ/Ｃp＝－Γdは,Ｔを仮温度Ｔ*に置き換えると近似的にそのまま成立すると考えれば,ｄＴ*/ｄｚ＝－ｍdｇ/Ｃp＝－Γdなる式が安定条件に対応します。
左辺のＴ*にＴ*＝Ｔ(1＋0.61ｒ)を代入すると,(1＋0.61ｒ)(ｄＴ/ｄｚ)＝－Γdです。そこで,不飽和で凝結のない大気の断熱減率ΓはΓ＝Γd/(1＋0.61ｒ)で与えられ,Γ≦Γdです。
さらに,飽和(saturate)して凝結水滴があり潜熱(latent heat)の存在によって温度降下が鈍くなる場合の湿潤断熱減率をΓwとします。
飽和水蒸気圧をｅsと書き,単位質量当たりの液体の水が水蒸気になるのに必要な潜熱(蒸発熱)をＬとすれば,ｅsと温度Ｔの関係は理論的にはClausius-Clapeyronの式:ｄｅs/ｄＴ＝Ｌｅs/(ＲＴ2)で表わされます。
(上式は,lnｅs＝－Ｌ/(ＲＴ)＋const.またはｅs＝Ａexp{－Ｌ/(ＲＴ)}を意味しますが,2006年７/30の記事「気液平衡の統計力学」ではこれより精密な修正式:lnｅs＝lnＴ－Ｌ/(ＲＴ)－1＋const.またはｅs＝ＡＴexp{－Ｌ/(ＲＴ)}を与えています。)
例えばＴ＝100℃ではＬ＝597(cal/ｇ)ですが常温付近でのＬはＴの変動に鈍感なのでＬがＴに無関係と近似してこれを積分すれば,ln(ｅs/ｅs0)＝(ｍvＬ/Ｒ)(1/Ｔ0－1/Ｔ)です。
これはｅsがＴだけに依存して周囲の大気圧ｐには依らない形の関係式です。
さて,乾燥空気と混合比ｑの水蒸気が混合した空気塊を考えます。
まず,単一気体の断熱変化は,0＝ｄlnθ＝(Ｃp/Ｒ)ｄlnＴ－ｄlnｐで表わされます。これはエネルギー収支を示す式と見るなら,0＝ＲＴｄθ/θ＝ＣpｄＴ－ＲＴｄｐ/ｐと表現されます。
そこで,乾燥空気と水蒸気が混合しただけの大気において水蒸気が飽和に達していないとき,各成分の上記左辺の温位変化:ＲＴｄθ/θへの寄与は,それぞれＣpｄＴ－ＲＴ(ｄｐd/ｐd)＝0,ｑＣpvｄＴ－ｑＲＴ(ｄｅ/ｅ)＝0 です。
ここで,ｑはこの大気の乾燥空気に対する水蒸気の混合比でｑ＝ρv/ρd＝(ｍd/ｍv)(ｅ/ｐd)＝(ｍd/ｍv){ｅ/(ｐ－ｅ)}です。また,Ｃpvは水蒸気の定圧モル比熱です。このとき,もちろんｐd＝ｐ－ｅです。
この空気塊が上昇してｑ＝ｑsとなって飽和した後,さらに上昇してｄｑsの水蒸気が凝結したとすれば,空気塊はＬｄｑsだけの潜熱を受けるので,変化が断熱なら水蒸気が飽和した空気塊のエネルギー収支の式は次のようになります。
すなわち,ＣpｄＴ＋ｑsＣpvｄＴ－ＲＴ{ｄ(ｐ－ｅs)/(ｐ－ｅs)＋ｑsｄｅs/ｅs}＋Ｌｄｑs＝0 となります。ｑsは飽和時の混合比です。
これに,ｄｅs/ｄＴ＝Ｌｅs/(ＲＴ2)を代入すれば,－Ｌｄ(ｑs/Ｔ)＝ＣpｄＴ/Ｔ＋ｑsＣpvｄＴ/Ｔ－Ｒｄ(ｐ－ｅs)/(ｐ－ｅs)を得ます。
ｐd＝ρdＲＴ/ｍd,ｅs＝ρvＲＴ/ｍvですから,静力学平衡:ｄｐ/ｄｚ＝－ρｇ, or ｄ(ｐd＋ｅs)/ｄｚ＝－(ρd＋ρV)ｇは,ｄ(ｐd＋ｅs)/ｄｚ＝－ｇ(ｍdｐd＋ｍvｅs)/(ＲＴ),またはｄｐ＝－ｇ{ｍd(ｐ－ｅs)＋ｍVｅs}/(ＲＴ)ｄｚと表現できます。
そこで,ｄ(ｐ－ｅs)＝－ｇ{ｍd(ｐ－ｅs)＋ｍVｅs}/(ＲＴ)ｄｚ－ｄｅsを－Ｌｄ(ｑs/Ｔ)＝(Ｃp＋ｑsＣpv)ｄＴ/Ｔ－Ｒｄ(ｐ－ｅs)/(ｐ－ｅs)に代入すると,－Ｌｄ(ｑs/Ｔ)＝(Ｃp＋ｑsＣpv)ｄＴ/Ｔ＋ｇ{ｍd＋ｍVｅs/(ｐ－ｅs)}ｄｚ/Ｔ＋Ｒｄｅs/(ｐ－ｅs)です。
これから,ｄＴ/ｄｚ＝－ｇ{ｍd＋ｍVｅs(ｐ－ｅs)}/(Ｃp＋ｑsＣpv)－ＲＴ(ｄｅs/ｄｚ)/{(ｐ－ｅs)(Ｃp＋ｑsＣpv)}＋{ＬＴ/(Ｃp＋ｑsＣpv)}{ｄ(ｑs/Ｔ)/ｄｚ}を得ます。
ところが,ｑs＝(ｍd/ｍv){ｅs/(ｐ－ｅs)}よりｅs/(ｐ－ｅs)＝ｍvｑs/ｍdです。
そこで,ｄｑs/ｄｚ＝(ｍd/ｍv)[(ｄｅs/ｄｚ)/(ｐ－ｅs)－{ｄ(ｐ－ｅs)/ｄｚ}{ｅs/(ｐ－ｅs)2}]＝(ｍd/ｍv)[(ｄｅs/ｄｚ)/(ｐ－ｅs)－{ｑs/(ｐ－ｅs)}{ｄ(ｐ－ｅs)/ｄｚ]と書けます。
つまり,(ｄｅs/ｄｚ)/(ｐ－ｅs)＝(ｍv/ｍd)(ｄｑs/ｄｚ)＋{ｑs/(ｐ－ｅs)}{ｄ(ｐ－ｅs)/ｄｚ}です。これらの式を用いて,できるだけｅsを消去します。
以上から,湿潤断熱減率は,Γw≡－ｄＴ/ｄｚ＝Γd/(1＋ｑsＣpv/Ｃp)＋ｍV2ｇ/{ｍdＣp(1＋ｑsＣpv/Ｃp)}＋ｍvＲＴ(ｄｑs/ｄｚ)/{ｍdＣp(1＋ｑsＣpv/Ｃp}＋ＲＴ{ｑs/(ｐ－ｅs)}{ｄ(ｐ－ｅs)/ｄｚ}/{Ｃp(1＋ｑsＣpv/Ｃp－ＬＴ{ｄ(ｑs/Ｔ)/ｄｚ}/{Ｃp(1＋ｑsＣpv/Ｃp)}です。
この式は閉じていないので古い正野重方著の岩波全書の「気象力学」などには近似公式もあるようですがここでは主要課題ではないので割愛します。ただ,一般に湿潤断熱減率は一定値ではありません。
(なお,水蒸気の定圧比熱については,2008年１/６の記事「氷,水,水蒸気の比熱」,および2009年２/９の記事「水蒸気の比熱」などを参照してください。) (補足終わり) ※
さて,浮力以外の気圧傾度力は無視できると考えて外気圧ｐは近似的に静力学平衡:ｄｐ/ｄｚ＝－ρｇを満足しているとします。
また,大気を流体と見るとき,対象とする運動でのReynolds数:Ｒeは非常に大きいため,1/Ｒeに比例する分子粘性項も無視します。
そして,プルームガスの周辺大気(ambient air)の局所速度,圧力,密度,温度,温位を,それぞれｖa,ｐa,ρa,Ｔa,θaとします。
特に,気圧については,∂ｐa/∂ｘ＝∂ｐa/∂ｙ＝0,∂ｐa/∂ｚ＝－ρaｇが成立しており,プルームガス圧については近似的にｐp＝ｐaが成立するとします。
この場合,周辺大気の運動方程式は,ρaｄｖa/ｄｔ＝－∇ｐa－ρaｇｋ＝0 です。ただしｉ,ｊ,ｋは,それぞれｘ,ｙ,ｚ軸の方向単位ベクトルを表わします。
この方程式の解は明らかにｖa＝constで与えられます。
一方,プルームガスの満たす運動方程式はρpｄｖp/ｄｔ＝－∇ｐp－ρpｇｋ～(－∇ｐa－ρaｇｋ)－(ρp－ρa)ｇｋです。
すなわち,この近似での運動方程式はρp(ｄｖp/ｄｔ)＝(ρa－ρp)/ｇｋ,またはｄｖp/ｄｔ＝{(ρa－ρp)/ρp}ｇｋです。これらの右辺が浮力項(buoyancy term)です。
状態方程式:ｐa＝ρaＲＴa/Ｍ,ｐp＝ρpＲＴp/Ｍから,ｐp～ｐaの近似ではρa/ρp ～Ｔp/Ｔa～θp/θaなので温位偏差θ'をθ'≡θp－θaと定義すれば,(ρa－ρp)/ρp～ (θp－θa)/θa＝θ'/θaです。
周辺大気の温位θaの定義:θa≡Ｔa(ｐ00/ｐa)(Ｒ/Ｃp)において特にｐ00～ｐaなるパラメータを採用すればθa～Ｔaなので,結局この近似(Boussinesq近似)での運動方程式は,(3)ｄｖp/ｄｔ＝(ｇθ'/Ｔa)ｋとなります。
連続方程式が質量の保存を意味するのに対して,運動方程式は運動量の保存を意味します。
さらに,ｖp≡ｖa＋ｖ'と書き,特に速度偏差:ｖ'の鉛直ｚ方向の成分(or第３成分)をｗ'とします。また,周辺大気の風速ｖaは鉛直成分を持たないとして,その向きをｘ軸(or１軸)の正の向きに取ります。
すると,ｖp≡ｕｉ,ｖp3＝ｖ'3＝ｗ'と書けます。さらにθa,ρa,ｕはｚだけの関数でありθa＝θa(ｚ),ρa＝ρa(ｚ),ｕ＝ｕ(ｚ)のように表わせるとします。
そして水平プルーム断面上ではｗ'と浮力ｇθ'/Ｔaが比例する:ｇθ'/Ｔa∝ ｗ'と仮定します。比例係数をｃ(ｚ)とするとｇθ'/Ｔa＝ｃ(ｚ)ｗ'です。
定常状態が実現されていると仮定し,面積分することによってプルームを水平面で切った断面上での平均量に対する方程式を求めます。
ただし,プルーム断面の外では明らかにθ'＝0,かつｖ'＝0 で特にｗ'＝0 です。
この結果は,まず浮力の効果についての方程式:(4)ｄＦz/ｄｚ＝－ｓＶを与えます。ただし,Ｆzは浮力流束で,(5)Ｆz≡∫P(ｇρpθ'ｗ'/Ｔa)ｄｘｄｙ/(πρa)で定義されます。
Ｖは鉛直容積流束で,(6)Ｖ≡∫Pρpｗ'ｄｘｄｙ/(πρa)です。
また,ｓは安定度パラメータと呼ばれる量で,(7)ｓ≡(ｇ/Ｔa)(∂θa/∂ｚ)です。
(証明):(2)ｄθp/ｄｔ＝0 より,∂θp/∂ｔ＋ｖp∇θp＝0ですから,定常状態:∂θp/∂ｔ＝0 ではρpｖp∇(θa＋θ')＝ 0です。
そして,∂θa/∂ｘ＝∂θa/∂ｙ＝∂θ'/∂ｘ＝∂θ'/∂ｙ＝0 よりρpｖp∇(θa＋θ')＝－ρpｗ'(∂θa/∂ｚ)です。
これに,(1)∇(ρpｖp)＝0 にθ'を掛けたものを加えると∇(ρpθ'ｖp)＝－ρpｗ'(∂θa/∂ｚ)を得ます。
これを全水平面で積分すると∫∞∇(ρpθ'ｖp)ｄｘｄｙ＝－(∂θa/∂ｚ)∫∞ρpｗ'ｄｘｄｙですが,無限遠では明らかにθ'＝0 で,プルーム断面を含む全水平面領域では全ての量はｚに依らず一定ですから,左辺の(∂/∂ｚ)も積分記号の外に出せます。
プルーム水平断面の外ではρpθ'ｗ'＝ρpｗ'＝0 なので,上式は(ｄ/ｄｚ)∫Pρpθ'ｗ'ｄｘｄｙ＝－(∂θa/∂ｚ)∫Pρpｗ'ｄｘｄｙとなります。
周辺大気の密度ρaは"上昇高度までの層当たり平均密度＝(一定)"と考え,両辺に(ｇ/Ｔa)/(πρa)を掛けると{∂(1/Ｔa)}{(ｄ/ｄｚ)∫Pρpｇθ'ｗ'ｄｘｄｙ}/(πρa)＝－(ｇ/Ｔa)(∂θ/∂ｚ){∫∞ρpｗ'ｄｘｄｙ}/(πρa)です。
よって,ｄＦz /ｄｚ＝－ｓＶ－(1/Ｔa)(ｄＴa/ｄｚ)Ｆzを得ました。
ところが,ｄＴa/ｄｚ～ｄθa/ｄｚ,Ｆz ～ｇＶ(＜θ'＞/θa）なので(1/Ｔa)(ｄＴa/ｄｚ)Ｆz～ (ｓＶ)・(＜θ'＞/θa)＜＜(ｓＶ)ですから,右辺第１項:(ｓＶ)に比べて右辺第２項:(1/Ｔa)(ｄＴa/ｄｚ)Ｆzを無視します。
そこで,結局ｄＦz /ｄｚ＝－ｓＶを得ます。(証明終わり)
次に,運動量に対する方程式は,(8)ｄ(ｖＶ)/ｄｚ＝(Ｆz/ｗ)ｋ＋ｖa(ｄＶ/ｄｚ)＋(抵抗力)です。ただし,ｖＶは(9)ｖＶ≡∫Pｖpρpｗ'ｄｘｄｙ/(πρa)で定義される運動量流束です。
この式で定義されるベクトルｖをプルームの速度と考えます。ｗはｖのｚ成分です。
(証明):(1)∇(ρpｖp)＝0 より,ｖp・∇(ρpｖp)＝0 ,つまりΣj=13ｖpi∂j(ρpｖpj)＝0 です。
また,(3)ｄｖp/ｄｔ＝(ｇθ'/Ｔa)ｋより,定常状態ではρp(ｖp∇)ｖp＝(ｇρpθ'/Ｔa)ｋ,つまりΣj=13ρpｖpj∂jｖpi＝(ｇρpθ'/Ｔa)δi3です。
したがって,∇(ρpｖpｖpi)＝Σj=13ｖpi∂j(ρpｖpj)＋Σj=13ρpｖpj∂jｖpi＝(ｇρpθ'/Ｔa)δi3です。
プルームの進行方向をｎとするとき,プルームの断面と側面で囲まれた領域ＶにGaussの定理を適用すると,∫V∇(ρpｖpｖpi)ｄＶ＝∫P+ρpｖpnｖpiｄＳ－∫P-ρpｖpnｖpiｄＳ－Δｎ{1＋(ｄｒ/ｄｎ)2}1/2(∫Cρpｖpｖpiｄｓ)です。
ただし,ここでのｒはプルーム半径(plume-radius)を意味します。プルーム断面は等方的なので,これは断面の形を円と想定したときの半径です。
－Δｎｄｓはプルーム側面の外向き法線の向きを持ちΔｎは側面ｎ方向の厚みです。それ故,ｄｓは周の微小線素の長さを持つ内向き法線方向ベクトルです。
故に,Δ∫PρpｖpnｖpiｄＳ≡∫P+ρpｖpnｖpiｄＳ－∫P-ρpｖpnｖpiｄＳ＝Δｎ{1＋(ｄｒ/ｄｎ)2}1/2(∫Cρpｖpｖpiｄｓ)＋δi3Δｎ∫P(ｇρpθ'/Ｔa)ｄＳです。
よって,(ｄ/ｄｎ)(∫PρpｖpnｖpiｄＳ)＝{1＋(ｄｒ/ｄｎ)2}1/2(∫Cρpｖpｖpiｄｓ)＋δi3∫P(ｇρpθ'/Ｔa)ｄＳを得ます。
一方,連続の方程式から(ｄ/ｄｎ)(∫PρpｖpnｄＳ)＝{1＋(ｄｒ/ｄｎ)2}1/2(∫Cρpｖpｄｓです。
さらにｖpi＝ｖai＋ｖ'iですから,(ｄ/ｄｎ)(∫PρpｖpnｖpiｄＳ)＝δi3∫P(ｇρpθ'/Ｔa)ｄＳ＋ｖai(ｄ/ｄｎ)(∫PρpｖpnｄＳ)＋{1＋(ｄｒ/ｄｎ)2}1/2(∫Cρpｖpｖ'iｄｓ)を得ます。
特にｎをｚ方向のｋに選べばｖpn＝ｗ'より,(ｄ/ｄｚ)(∫Pρpｗ'ｖpiｄＳ)＝δi3∫P(ｇρpθ'/Ｔa)ｄＳ＋ｖi(ｄ/ｄｚ)(∫Pρpｗ'ｄＳ)＋(抵抗力)です。
ここで,(抵抗力)≡{1＋(ｄｒ/ｄｚ)2}1/2(∫Cρpｖpｖ'iｄｓ)であり,ｄｓ＝{1＋(ｄｒ/ｄｚ)2}-1/2(－ｄｙｉ＋ｄｘｊ)＋{1＋(ｄｒ/ｄｚ)2}-1/2(ｄｒ/ｄｚ)ｄｌｋです。(ｄｌ＝(ｄｘ2＋ｄｙ2)1/2)
一方,∫P(ｇρpθ'/Ｔa)ｄＳ＝(1/ｗ)∫P(ｇ/Ｔa)(ρpθ'ｗ)ｄＳ＝(1/ｗ)∫P(ｇ/Ｔa)(ρpθ'ｗ')ｄＳ＋(1/ｗ)∫P(ｇ/Ｔa){ρpθ'(ｗ－ｗ')}ｄＳです。
ｖＶ≡∫Pｖpρpｗ'ｄＳ/(πρa),Ｖ≡∫Pρpｗ'ｄＳ/(πρa),およびｖp3＝ｗ'よりｗ＝∫Pρpｗ'2ｄＳ/∫Pρpｗ'ｄＳですから,ｗ－ｗ'＝{∫Pρpｗ'2ｄＳ－ｗ'∫Pρpｗ'ｄＳ}/{∫Pρpｗ'ｄＳ}です。
仮定によりｇθ'/Ｔa＝ｃ(ｚ)ｗ'を代入すると,(1/ｗ)∫P(ｇ/Ｔa){ρpθ'(ｗ－ｗ')}ｄＳ＝{ｃ(ｚ)/ｗ}[{∫Pρpｗ'ｄＳ}{∫Pρpｗ'2ｄＳ}－{∫Pρpｗ'2ｄＳ}{∫Pρpｗ'ｄＳ}]/{∫Pρpｗ'ｄＳ}＝0 です。
したがって,∫P(ｇρpθ’/Ｔa)ｄＳ＝(1/ｗ)∫P(ｇ/Ｔa)(ρpθ'ｗ')ｄＳ故,(ｄ/ｄｚ)(∫Pρpｗ'ｖpiｄＳ)＝(1/ｗ)δi3∫P(ｇ/Ｔa)(ρpθ'ｗ')ｄＳ＋ｖi(ｄ/ｄｚ)(∫Pρpｗ'ｄＳ)＋(抵抗力),
または,(ｄ/ｄｚ)(∫Pｖpρpｗ'ｄＳ)＝(1/ｗ){∫P(ｇ/Ｔa)(ρpθ'ｗ')ｄＳ}ｋ＋ｖa(ｄ/ｄｚ)(∫Pρpｗ'ｄＳ)＋(抵抗力)なる式を得ます。
最後の表式は,まさにｄ(ｖＶ)/ｄｚ＝(Ｆz/ｗ)ｋ＋ｖa(ｄＶ/ｄｚ)＋(抵抗力)そのものです。
※特にプルームが軸対称で完全に鉛直向きなら(抵抗力)＝{1＋(ｄｒ/ｄｚ)2}1/2(∫Cρpｖpｖ'iｄｓ)＝0でｄ(ｗＶ)/ｄｚ＝Ｆz/ｗです。※(証明終わり)
さて,ｄＦz/ｄｚ＝－ｓＶ,ｄ(ｗＶ)/ｄｚ＝Ｆz/ｗを解くために次のような初期条件を与えます。
まず,初期速度は吐出速度に等しいとします。すなわち,(10a)ｖ＝ｗ0ｋです。そして初期運動量流束は,(10b)ｖＶ＝(ｗＶ)ｋ＝Ｆmｋ≡(ρ0/ρa)ｗ02ｒ02ｋとします。ρ0は排出ガス密度,ｒ0は煙突口の口径です。
さらに,初期の浮力流束は,(10c)Ｆz＝Ｆ≡(1－ρ0/ρa)ｇｗ0ｒ02で与えられるとします。
つまり,初期容積流束をＶ0≡ｗ0ｒ02とすれば初期の総質量流束はπρ0Ｖ0＝πρ0ｗ0ｒ02ですから,定義により運動量流束はＦm＝πρ0ｗ0Ｖ0/(πρa)＝(ρ0/ρa)ｗ0Ｖ0です。
そして,初期浮力流束はＦz＝Ｆ＝(1－ρ0/ρa)ｇＶ0～{(Ｔ0－Ｔa)/Ｔ0}ｇＶ0となるわけです。
ところで,排出源がhot-sourceの場合の初期浮力流束は,(11)Ｆ＝ｇＱH/(πρaｃpＴa)で与えられます。
ここで,ＱHは排出されるガスによって単位時間に排出される熱量(cal/sec)であり,ｃpは空気の単位質量当たりの熱容量(定圧比熱)(cal/(kg・Ｋ))です。
なぜなら,Ｆz＝Ｆ＝(ｇθ'/Ｔa)(πρ0ｗ0ｒ02)/(πρa)＝{ｇ/(πρaＴa)}(θ'πρ0Ｖ0)ですが,πρ0Ｖ0は総質量流束なので排出ガスによって単位時間に外気に付加される熱量はＱH＝ｃpＴ'πρ0Ｖ0 ～ｃpθ'πρ0Ｖ0に等しいからです。
つまり,θ'πρ0Ｖ0＝ＱH/ｃpによってＦ＝ｇＱH/(πρaｃpＴa)が得られるわけです。
この場合,排熱強度ＱHと排出ガス密度ρ0の情報を得ることにより逆にＶ0＝ｗ0ｒ02が求まります。
そこで,煙突口径ｒ0のデータから吐出速度ｗ0が算定されて,結局全ての初期情報を得ることができます。
ただし,後述するように３つの未知数Ｆz,ｗ,Ｖに対し微分方程式は２つしかないので解を決めるにはもう１つ方程式が必要です。
今日はここまでですがまだ続きます。
参考文献：G.A.Briggs "Plume Rise",U.S.Atomic Energy Commision Division of Technical Information(1969),小倉義光「一般気象学」(東京大学出版会),小倉義光「気象力学通論」(東京大学出版会)
最近のコメント