移流拡散方程式を解く

　特殊な条件の下で移流拡散方程式を解いてみたいと思います。これは30年ほど前に就職した会社で配属後に計算したものです。 

　点源から不断にＱ(ｍ³/ｓ)の気体物質が排出され,風速が一様でｖ＝(ｕ,0,0),拡散係数:Ｋ_x,Ｋ_y,Ｋ_zが定数である場合の濃度計算を考えます。まずは,点源が原点にあるとします。 

　このとき,濃度Ｃに対する拡散方程式は∂Ｃ/∂ｔ＋ｕ(∂Ｃ/∂ｘ)＝Ｋ_x(∂²Ｃ/∂ｘ²)＋Ｋ_y(∂²Ｃ/∂ｙ²)＋Ｋ_z(∂²Ｃ/∂ｚ²)＋Ｑδ³(ｒ)と書けます。

　これを解くため,Ｃをフーリエ(Fourier)変換して,Ｃ(ｒ,ｔ)≡{1/(2π)³}∫ｄ³ｋＣ^(ｋ,ｔ)exp(ｉｋｒ)とします。

　ｒ≠0 ではＣ^＝Ｃ^(ｋ,ｔ)について∂Ｃ^/∂ｔ＋(Ｋ_xｋ_x²＋Ｋ_yｋ_y²＋Ｋ_zｋ_z²＋ｉｋ_xｕ)Ｃ^＝0 となります。

　故に,Ｃ^(ｋ,ｔ)＝Ｃ^(ｋ,ｔ₀)exp[－(Ｋ_xｋ_x²＋Ｋ_yｋ_y²＋Ｋ_zｋ_z²＋ｉｋ_xｕ)(ｔ－ｔ₀)]と書くことができます。

　そして,Ｃ^(ｋ,ｔ₀)＝∫ｄ³ｒ'Ｃ(ｒ'＋ｒ,ｔ₀)exp[－ｉｋ(ｒ'＋ｒ)]ですからＣ(ｒ,ｔ)＝{1/(2π)³}∫ｄ³ｋＣ^(ｋ,ｔ₀)exp[－(Ｋ_xｋ_x²＋Ｋ_yｋ_y²＋Ｋ_zｋ_z²＋ｉｋ_xｕ)(ｔ－ｔ₀)＋ｉｋｒ]＝{1/(2π)³}∫ｄ³ｒ'Ｃ(ｒ'＋ｒ,ｔ₀)∫ｄ³ｋexp[－(Ｋ_xｋ_x²＋Ｋ_yｋ_y²＋Ｋ_zｋ_z²)(ｔ－ｔ₀)－ｉ{ｋｒ'－ｋ_xｕ(ｔ－ｔ₀)}]となります。 

　ところが,∫ｄ³ｋexp[－(Ｋ_xｋ_x²＋Ｋ_yｋ_y²＋Ｋ_zｋ_z²)(ｔ－ｔ₀)－ｉ{ｋｒ'－ｋ_xｕ(ｔ－ｔ₀)}]＝∫exp[－Ｋ_x(ｔ－ｔ₀){ｋ_x＋ｉ[ｘ'－ｕ(ｔ－ｔ₀)]/[2Ｋ_x(ｔ－ｔ₀)]}²]ｄｋ_x∫exp[－Ｋ_y(ｔ－ｔ₀){ｋ_y＋ｉｙ'/[2Ｋ_y(ｔ－ｔ₀)]}²]ｄｋ_y∫exp[－Ｋ_z(ｔ－ｔ₀){ｋ_z＋ｉｚ'/[2Ｋ_z(ｔ－ｔ₀)]}²]ｄｋ_zexp[－{ｘ'－ｕ(ｔ－ｔ₀)}²/{4Ｋ_x(ｔ－ｔ₀)}－ｙ'²/{4Ｋ_y(ｔ－ｔ₀)}－ｚ'²/{4Ｋ_z(ｔ－ｔ₀)}]です。

　ガウス積分を実行することにより,右辺＝[π³/{Ｋ_xＫ_yＫ_z(ｔ－ｔ₀)³}]^1/2exp[－{ｘ'－ｕ(ｔ－ｔ₀)}²/{4Ｋ_x(ｔ－ｔ₀)}－ｙ'²/{4Ｋ_y(ｔ－ｔ₀)}－ｚ'²/{4Ｋ_z(ｔ－ｔ₀)}]となります。

　したがって,Ｃ(ｒ,ｔ)＝(1/8)[1/{π³Ｋ_xＫ_yＫ_z(ｔ－ｔ₀)³}]^1/2∫ｄ³ｒ'Ｃ(ｒ'＋ｒ,ｔ₀)exp[－{ｘ'－ｕ(ｔ－ｔ₀)}²/{4Ｋ_x(ｔ－ｔ₀)}－ｙ'²/{4Ｋ_y(ｔ－ｔ₀)}－ｚ'²/{4Ｋ_z(ｔ－ｔ₀)}]です。

　特に初期瞬時濃度が排出強度Ｑ'そのもの,つまりＣ(ｒ,ｔ₀)＝Ｑ'δ³(ｒ)なら,Ｃ(ｒ,ｔ)＝(Ｑ'/8)[1/{π³Ｋ_xＫ_yＫ_z(ｔ－ｔ₀)³}]^1/2exp[－{ｘ－ｕ(ｔ－ｔ₀)}²/{4Ｋ_x(ｔ－ｔ₀)}－ｙ²/{4Ｋ_y(ｔ－ｔ₀)}－ｚ²/{4Ｋ_z(ｔ－ｔ₀)}]となりますね。

　そこで,排出量が単位時間当たりＱ,つまり,微小時間ｄｔ₀ではＱ'＝Ｑｄｔ₀である場合,これを代入してｔ₀で積分します。

　定常状態の濃度Ｃ(ｒ)はこれで得られるはずです。

　Ｃ(ｒ)＝∫₀^∞ｄｔ₀(Ｑ/8)[1/{π³Ｋ_xＫ_yＫ_z(ｔ－ｔ₀)³}]^1/2exp[－{ｘ－ｕ(ｔ－ｔ₀)}²/{4Ｋ_x(ｔ－ｔ₀)}－ｙ²/{4Ｋ_y(ｔ－ｔ₀)}－ｚ²/{4Ｋ_z(ｔ－ｔ₀)}]と書けます。

　この時間積分を実行すると,結局Ｃ(ｒ)＝(Ｑ/4)[1/{π²Ｋ_xＫ_yＫ_z}]^1/2(ｘ²/Ｋ_ｘ＋ｙ²/Ｋ_y＋ｚ²/Ｋ_z)^-1/2 exp{ｕｘ/(2Ｋ_ｘ)}exp(－ｕ(ｘ²/Ｋ_ｘ＋ｙ²/Ｋ_y＋ｚ²/Ｋ_z)^1/2/(2Ｋ_ｘ^1/2))となります。

　特にｘ,ｙを水平方向,ｚを鉛直上向き方向としたとき,一般に地上付近を考えるならばｘ,ｙ＞＞ｚであり,水平方向の対称性から,Ｋ_ｘ＝Ｋ_y＝Ｋ_H とすると,(ｘ²/Ｋ_ｘ＋ｙ²/Ｋ_y＋ｚ²/Ｋ_z)^1/2～ Ｒ/Ｋ_Hと考えることができます。

　ただし,Ｒ≡(ｘ²＋ｙ²)^1/2は点源からの水平距離です。 

　そこで,Ｃ(ｒ) ～ [Ｑ/(4π)][1/(ＲＫ_HＫ_z^1/2)]exp {ｕｘ/(2Ｋ_H)－ｕ{Ｒ²/(2Ｋ_H)²＋ｚ²/(4Ｋ_zＫ_H)^1/2}]です。

　さらに風下方向ｘがｙよりも風に流される影響が優勢であると考えてｘ＞＞ｙとすると,ｕｘ/(2Ｋ_H)－ｕ{Ｒ²/(2Ｋ_H)²＋ｚ²/(4Ｋ_zＫ_H)^1/2} ～ ｕｘ/(2Ｋ_H)－ｕｘ/(2Ｋ_H)[1＋ｙ²/(2ｘ²)＋Ｋ_Hｚ²/(2Ｋ_zｘ²)]と近似すれば,Ｃ(ｒ) ～ [Ｑ/(4π)][1/(ｘＫ_HＫ_z^1/2)]exp[－ｕｙ²/(4Ｋ_Hｘ)－ｕｚ²/(4Ｋ_zｘ)]となります。

ここで,σ_x²＝σ_y²＝2Ｋ_Hｘ/ｕ,σ_z²＝2Ｋ_zｘ/ｕと書けば,Ｃ(ｒ)～[Ｑ/(2πσ_xσ_zｕ)]exp[－ｙ²/(2σ_y²)－ｚ²/(2σ_z²)]となります。

しかし,一般には地上ｚ＝0 より下のｚ＜0 には拡散しないのでｚ＝0 では境界条件として∂Ｃ/∂ｚ＝0 なる反射条件が成立します。

規格化定数は上半面ｚ＞0 のみを考えるので,先の値の２倍となります。また,拡散式は発生点源より風上のｘ≦0 ではＣ(ｒ)＝0 であるという境界条件も満足する必要があります。

もし,点源が高さｚ＝Ｈ_e＝(有効高さ)を持つなら,地下ｚ＝－Ｈ_eに仮想の鏡像点源があると考えればよく,濃度はＣ(ｒ) ～ [Ｑ/(2πσ_xσ_zｕ)]exp{－ｙ²/(2σ_y²)}[exp|－(ｚ－Ｈ_e)²/(2σ_z²)＋exp|－(ｚ＋Ｈ_e)²/(2σ_z²)]]と変更されます。

行政では,簡易式として,ｙ方向について－∞ から＋∞ まで積分し,その濃度が円周2πＲではなく16等分した風向の風下扇形(π/8)Ｒ部分に全て分配されると仮定して近似した式が採用されているようです。

すなわち,濃度は風向をｘ軸としてＣ(ｒ)～Ｃ(Ｒ,ｚ)～ [Ｑ/{(2π)^1/2(π/8)Ｒσ_zｕ}][exp|－(ｚ－Ｈ_e)²/(2σ_z²)＋exp|－(ｚ＋Ｈ_e)²/(2σ_z²)]]で与えられるという簡易式です。

これらの式は有風時のプルーム式と呼ばれています。

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

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

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

    物理学

大気中の移流拡散方程式

　大気中の気体物質の拡散方程式について若干の考察をしてみたいと思います。 

　ある温度で大気の単位体積中に存在する同じ温度の気体物質の体積の値をその濃度と呼びＣで表わすことにします。

　"(大気＋物質)の密度＝大気全体の密度"をρ,濃度対象としている気体物質のみの密度をρ₁とし,質量濃度をｃとします。ρ₁＝ρｃです。

　そして,大気全体には"主流＝平均流"として風速ｖがあるとし,大気全体は非圧縮,かつ空間的にもρは一定と仮定すると,その全体としての連続の方程式はdivｖ＝∇ｖ＝0 となります。

　一方,混合された気体物質に着目した連続の方程式には湧き出しがあるとして,∂ρ₁/∂ｔ＋∇(ρ₁ｖ₁)＝ｑ₁と書くことができます。湧き出しｑ₁は物質の排出強度と呼ばれます。 

　連続の方程式∂ρ₁/∂ｔ＋∇(ρ₁ｖ₁)＝ｑ₁にρ₁＝ρｃを代入すると,∂(ρｃ)/∂ｔ＋∇(ρｃｖ)＝－∇[ρｃ(ｖ₁－ｖ)]となります。

　結局,Ｄｃ/Ｄｔ＝∂ｃ/∂ｔ＋ｖ∇ｃ＝－∇ｉ＋ｑ₁ /ρと書けます。ただし,Ｄ/Ｄｔ＝∂ｃ/∂ｔ＋ｖ∇はラグランジュ微分で,ｉ≡ｃ(ｖ₁－ｖ)は拡散流束密度と呼ばれる量です。

　ここで,Ｎ₀をアボガドロ数,Ｍを大気全体の分子量,Ｍ₁を濃度対象の気体物質の分子量とすると,大気全体と対象物質の単位体積当たりの分子数は,それぞれ,Ｎ₀ρ/Ｍ,Ｎ₀ρ₁/Ｍ₁です。

　同一温度での気体の体積は分子の数に比例するので,Ｃ＝ρ₁Ｍ/(ρＭ₁)＝(Ｍ/Ｍ₁)ｃとなります。よって,質量濃度ｃと(体積)濃度Ｃは単に定数係数の違いしかありません。

　したがって,質量濃度ｃの方程式Ｄｃ/Ｄｔ＝－∇ｉ＋ｑ₁ /ρは体積濃度Ｃで表わすと,ＤＣ/Ｄｔ＝－∇ｊ＋Ｑとなります。

　ｊ≡(Ｍ/Ｍ₁)ｉ＝Ｃ(ｖ₁－ｖ)は(体積濃度の)拡散流束密度,Ｑ≡ｑ₁Ｍ/(ρ/Ｍ₁)は物質の(体積)排出強度と解釈することができます。

　ここで,フィック型の拡散を想定し拡散流束密度ｊが濃度の高い方から低い方へその勾配に比例して流れるとします。

　その拡散流束密度は拡散係数Ｋをテンソル{Ｋ_ij}としてｊ_i＝－Ｋ_ij(∂Ｃ/∂ｘ_j)となります。Ｋは拡散係数ですから正値です。

　ただし,同じ添字が２度出現するときは１から３まで加える,というアインシュタインの規約を用います。

　拡散係数テンソルＫを主軸変換して対角化し,Ｋ_ij＝Ｋ_iδ_ijの主軸をｘ,ｙ,ｚ方向に取れば,その主軸成分はＫ_x,Ｋ_y,Ｋ_z とすることができます。

　このときには,拡散流束密度はｊ＝{Ｋ_x(∂Ｃ/∂ｘ)}ｅ_x＋{Ｋ_y(∂Ｃ/∂ｙ)}ｅ_y＋{Ｋ_z(∂Ｃ/∂ｚ)}ｅ_zです。

　そこで,拡散方程式:∂Ｃ/∂ｔ＋ｖ∇Ｃ＝(∂/∂ｘ){Ｋ_x(∂Ｃ/∂ｘ)}＋(∂/∂ｙ){Ｋ_y(∂Ｃ/∂ｙ)}＋(∂/∂ｚ){Ｋ_z(∂Ｃ/∂ｚ)}＋Ｑが得られます。

　特に拡散係数Ｋがスカラーで近似的に定数と見なせるときは,Ｋ_ij＝Ｋδ_ijであり,通常の移流拡散方程式である∂Ｃ/∂ｔ＋ｖ∇Ｃ＝Ｋ∇²Ｃ＋Ｑに帰着します。 

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

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

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

    物理学

線形微分方程式系の直接解法(ボックスモデル)

　今日は屋内での汚染物質濃度を求めることを仮定した簡易ボックスモデルを例にとって定数係数の連立線形常微分方程式の数値解を求めるための１つの手法を紹介します。

　室内の濃度分布を計算するには簡易モデルではなくて熱対流などの気流を計算する２方程式モデルなどに移流拡散方程式を付加した非定常な数値流体方程式を解くという大げさな方法もあります。

　しかし,ここでは部屋の平均濃度だけを計算する,という簡易ボックスモデルを考察します。

　例えば容積がＶ(ｍ³)の部屋の中央に石油ストーブなどがあって,それからの一酸化炭素ＣＯの排出強度がＱ(ｍ³/s)であるとします。もし部屋が密閉されていれば時間ｔが経った後にそのＣＯの平均濃度をＣとするとＣ＝(Ｑ/Ｖ)ｔとなると考えていいでしょう。

　しかし一般に部屋には隙間があって換気されるわけで,通常１時間に部屋全体の空気が入れ替わる回数を換気率という量で表わします。

　例えば容積Ｖの部屋の換気率がｎならば,１時間の後には総量ｎＶの空気が排出され物質量の保存則によって同時にｎＶの空気が外部から入ってくるわけですから,１秒当たりの空気の流出入量をｋ(ｍ³/s)とするとｋ＝ｎＶ/3600ですね。

そこで,各時刻の部屋の中のＣＯの平均濃度をＣ(ｔ)で表わすとこの部屋の中の単位体積当たりのＣＯの量は各時刻にＣ(ｔ)(ｍ³/ｍ³)なので時刻Δｔの間に外気に流出するＣＯの量はｋＣ(ｔ)Δｔです。

一方外気の一酸化炭素濃度(一定:constant)をＣ₀とすると,流入量はｋＣ₀Δｔですから,部屋内のＣＯ量の時刻変化はＶ(Ｃ(ｔ＋Δｔ)－Ｃ(ｔ))＝ｋ[Ｃ₀－Ｃ(ｔ)]Δｔという式で与えられます。

　さらに排出強度をＱとすると,右辺には発生量ＱΔｔが加わりますから,結局ＣＯ濃度に対する微分方程式はｄＣ/ｄｔ＝ａ(Ｃ₀－Ｃ)＋Ｑ/Ｖと表わすことができます。ここでａ≡ｋ/Ｖです。

　Ｂ＝ａＣ₀＋Ｑ/ＶとおけばｄＣ/ｄｔ＝－ａＣ＋Ｂとなりますが,この微分方程式は簡単に解けてＣ(ｔ)＝(Ｃ( 0 )－ａ^-1Ｂ)exp(－ａｔ)＋ａ^-1Ｂという解が得られます。

こうしたモデルを室内物質濃度の簡易ボックスモデルと言います。

　これを一般化して各部屋の容積がＶ_i(ｉ＝1,2,...ｍ)のｍ個の部屋を有する住宅があって,これらの部屋の平均濃度が(ｍ＋1)次元の列ベクトルでＣ(ｔ)＝^t(Ｃ₀,Ｃ₁,Ｃ₂,...,Ｃ_m)で表わされるとします。

　ここでＣ_kはｋ番目の部屋における物質濃度です。ただしＣ₀は屋外の平均濃度値を示しています。

そして各部屋には排出強度Ｑ_iの排出源があり"部屋iと部屋ｊの間の換気量＝先の方程式のｋに相当する値"をｋ_ij＝ｋ_ji(i≠ｊ;ｉ,ｊ＝0,1,2,．．ｍ)とすれば,このボックスモデルの微分方程式はｄＣ_i/ｄｔ＝∑_j=0^mａ_ij(Ｃ_j－Ｃ_i)＋Ｑ_i/Ｖ_i(ｉ＝0,1,2,．．ｍ)となります。ただしａ_ij≡(ｋ_ij/Ｖ_i)です。

　(排出強度/容積)の列ベクトルを,Ｂ≡^t(Ｑ₀/Ｖ₀,Ｑ₁/Ｖ₁Ｑ₂/Ｖ₂,．．,Ｑ_m/Ｖ_m)で定義します。ただしＱ₀/Ｖ₀＝0 としておきます。

　また,λ_i≡∑_j=0^mａ_ijと定義し行列Ａを成分(ａ_ij－λ_iδ_ij)によって定義すれば,濃度を求めるボックスモデルの微分方程式,すなわち(ｍ＋1)変数の定数係数連立線形非同次微分方程式は(ｍ＋1)次元ベクトル空間での線形非同次微分方程式としてｄＣ/ｄｔ＝－ＡＣ＋Ｂと表わされます。 　

　この方程式の形式的な解は簡単に求めることができて初期値をＣ(0)とすると,Ｃ(ｔ)＝{Ｃ(0)－Ａ^-1Ｂ}exp(－Ａｔ)＋Ａ^-1Ｂと書くことができます。

　さらに濃度を形式解ではなくて実際に具体的に求めるには従来の慣用的な方法であればＡの最大(ｍ＋1)個の固有値:λ_i(ｉ＝0,1,2,．．ｍ)を全て求めて,Ｃ(ｔ)の各成分を具体的にexp(－λ_iｔ)の１次結合で表わすのです。

　しかし,たかだか10部屋程度の住宅で行列Ａが具体的にわかっているときには,固有値を求めるなどという面倒な手続きを省略してコンピュータで無理矢理,力技で行列の指数関数を計算してしまう,という方法があります。

　すなわち,形式解Ｃ(ｔ)＝{Ｃ(0)－Ａ^-1Ｂ}exp(－Ａｔ)＋Ａ^-1Ｂにおいて逆行列:Ａ^-1はＡからコンピュータで掃き出し法などで計算可能ですし,exp(－Ａｔ)＝∑_ν=0^∞(－Ａｔ)^ν/ν!であって,コンピュータで行列のベキ乗(power):(－Ａｔ)^νを具体的に計算することも可能です。

　そこで,左辺のexp(－Ａｔ)を右辺の∑_ν=0^∞のν＝0,1,2,...の最初の数項で近似評価することによって力技で数値的にＣ(ｔ)の近似解を求めることができるわけです。

　逆にこの近似解と物質濃度の直接測定実験の結果から非線形最小二乗法などによって換気率パラメータ:ｋ_ijを推定することもできます。

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

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

非線形最小二乗法（JEA式の作成過程)

　今日は窒素酸化物濃度などを計算する解析的な低煙源拡散式のパラメータの推定のために用いたことがある非線形最小二乗法について解説してみたいと思います。(直角風のJEA式)

　地上の道路上を走行する自動車からの窒素酸化物(NOx)などの汚染物質排出を,有効煙源高さｈe＝0(ｍ)のｙ軸上で単位長さ当たりの定常的な排出強度がＱ(Ｎｍ³/(ｍｓ)の無限線煙源としてモデル化します。

　そして,道路に直角のｘ方向に風が吹いている設定で,風速も拡散係数も高さｚ(ｍ)のベキ乗則に従って鉛直上方に向かって増加するというモデルを想定します。

　このときの拡散方程式の解である汚染物質濃度を,垂直距離ｘ(ｍ)と高さｚ(ｍ)の関数と考えてＣ(ｘ,ｚ)で表わすと,これはＡをＱに比例するあるパラメータとして,Ｃ(ｘ,ｚ)＝Ａｘ^-sexp(－Ｂｚ^p/ｘ)なる形の式(Robertsの式)で表わされる,ことがわかります。

　様々な環境の下で数回に渡って行なわれた実際の実験時の煙源長さはもちろん有限なのですが,とりあえず,これを無限大長さで近似します。

　実はガウスの誤差関数でもって有限長さの効果を取り入れることもできるのですが,今日の話題では割愛します。

　実験はトレーサーガスを地上にある有限長さの直線状のパイプに均等間隔で開けた穴から拡散させるというものです。

そのパイプ状の発生源をｙ軸としたとき,ｎ個の観測点座標:(ｘ_i,ｚ_i)(ｉ＝1,2,．．．ｎ)において,風向が直角に近い環境のときの実測濃度:ｃ_iと先に設定した低煙源拡散式による計算値:Ｃ_i＝Ｃ(ｘ_i,ｚ_i)とを比較することによって,逆にその式のパラメータＡ,ｓ,Ｂ,ｐを推定します。

　具体的には,Ｃ(ｘ,ｚ)＝Ａｘ^-sexp(－Ｂｚ^p/ｘ)の右辺をｆ(Ａ,ｓ,Ｂ,ｐ；ｘ,ｚ)と書いて,Ｃ_i＝ｆ(Ａ,ｓ,Ｂ,ｐ；ｘ_i,ｚ_i)とし,計算値と実測値の誤差の二乗和:Ｓ(Ａ,ｓ,Ｂ,ｐ)≡∑(Ｃ_i－ｃ_i)²を最小にするパラメータＡ,ｓ,Ｂ,ｐを求めるわけです。

　これは次に述べる非線形最小二乗法のパラメータ数が４個のときの例になっています。

　ここでは,より一般的に未知パラメータがｒ個あるとし,それらを順にＡ₁, Ａ₂,．．Ａ_rとします。

　上のケースではｒ＝4で,Ａ₁＝Ａ,Ａ₂＝ｓ,Ａ₃＝Ｂ,Ａ₄＝ｐです。

　そして,Ｃ_i＝Ｃ(ｘ_i,ｚ_i)＝ｆ(Ａ₁,Ａ₂,．．Ａ_r；ｘ_i,ｚ_i)とし,Ｓ(Ａ₁,Ａ₂,．．Ａ_r)≡∑(Ｃ_i－ｃ_i)²と定義するわけです。

　誤差の二乗和Ｓが最小となる条件は通常の線形回帰の計算の場合と同様,Ｆ_i(Ａ₁,Ａ₂,．．Ａ_r)≡∂Ｓ/∂Ａi＝2∑(Ｃ_k－ｃ_k)∂Ｃ_k/∂Ａi＝0 (ｉ＝1,2,．．ｒ)で与えられますから,この r 個の連立方程式を解くことにより,未知数Ａ₁,Ａ₂,．．Ａ_rを求めることが主目的となります。

　これらの方程式は一般に非線形ですから,こうした非線形回帰によるパラメータ推定の方法を非線形最小二乗法と呼ぶことにします。

　具体的には,Ｆ_i(Ａ₁,Ａ₂,．．Ａ_r)＝0 (ｉ＝1,2,．．ｒ)を線形近似することにより多変数ニュートン法(Newtonian method)を実行します。

　そこで,まず連立方程式を線形近似します。

　すなわち,初期値としてＡ_j＝Ａ_j⁰ を適当に与えた後に,　0 ＝ Ｆ_i(Ａ₁,Ａ₂,．．Ａ_r)＝Ｆ_i(Ａ₁⁰,Ａ₂⁰,．．Ａ_r⁰)＋∑(∂Ｆ_i/∂Ａ_j)|_Ａ＝Ａ0(Ａ_j－Ａ_j⁰)と近似します。

これを行列形式で書くため,Ｄという行列をＤ＝(ｄ_ij)≡(∂Ｆ_i/∂Ａ_j)で定義し,特にＤ₀＝(ｄ_ij⁰)≡(∂Ｆ_i/∂Ａ_j)|_Ａ＝Ａ0とします。

一方,列ベクトルＡをＡ≡^t(Ａ₁,Ａ₂,．．Ａ_r)で定義し,特にＡ₀≡^t(Ａ₁⁰,Ａ₂⁰,．．Ａ_r⁰)とします。さらにＦ_i(Ａ₁⁰,Ａ₂⁰,．．Ａ_r⁰)を成分とする同様な列ベクトルをＦ₀と定義します。

こうすれば,先の線形近似は 0＝Ｆ₀＋Ｄ₀(Ａ－Ａ₀)と表わされます。

これを単純に解くと,Ａ＝Ａ₀－Ｄ₀^-1Ｆ₀と近似されることになります。ここで,Ｄ₀^-1はＤ₀の逆行列です。

得られたＡ＝Ａ₀－Ｄ₀^-1Ｆ₀を,改めてＡ₀として代入してＤ₀^-1Ｆ₀を計算し,Ａ＝Ａ₀－Ｄ₀^-1Ｆ₀が収束するまで繰り返します。

すなわち,Ａ_m+1＝Ａ_m－Ｄ_m^-1Ｆ_mの漸化式において誤差|Ａ_m+1－Ａ_m|の相対値が十分小さくなるまで繰り返し計算します。

こうすれば,ｍ→ ∞ でのＡ_m＝^t(Ａ₁^m,Ａ₂^m,．．Ａ_r^m)の極限値としてＦ_i(Ａ₁,Ａ₂,．．Ａ_r)＝0 (ｉ＝1,2,．．．ｒ)を満たすＡ＝^t(Ａ₁,Ａ₂,．．Ａ_r)が得られるはずです。

さらに,実際には収束を速くするために加速係数ｗを与えてＡ_m+1＝Ａ_m－ｗＤ_m^-1Ｆ_mとする方法を用いたほうがいいと思います。

実際のｒ＝4 の計算では,かつてFortranを使ってプログラミングしましたが,初期値を適切に取ると各環境のケースについて大体十数回の繰り返し計算でパラメータの推定値が得られました。

そして濃度の実測値とその推定パラメータによる計算値との相関係数としては,0.9前後の値が得られ,回帰係数の値も0.8 から 1.2程度となって,仮定した計算式が良い近似になることがわかりました。

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

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

低煙源拡散モデル（JEA式）

　与太話が続いているので,ここらで少しアカデミック(Academic)な話をします。

　アカデミックといえば以前,会社でNHC（練馬変態クラブ）を創設した私より年齢が１つ上のフィクサーのような人がいまして,私のことを「チミはアカデミックというよりは垢デミック（垢だらけ）だね」と評されました。

　確かに当時は学生時代の延長で,私は会社で２番目に不潔だと言われていましたから身に覚えがあることです。

　彼と私は共にシュールレアリストを標榜していまして,彼は「シュールレアリストは結婚しない。」というのが口癖でした。私が40歳で退職した翌年に高血圧なのに毎日,飲酒し,しかも過労であったということが原因か,北区のアパートで亡くなっているのが尋ねた人によって発見されました。

　彼(シイタケマン)こそ本当の破滅型人間でしたね。

   本題に入ります。

　以前,4月中旬にS社での「会社員新人時代の思い出」という内容で記事を書きましたが,当時,配属されてすぐの仕事はコンピュータを使うことではなくて,後に環境庁（現在は環境省）の窒素酸化物総量規制マニュアルでJEA式（日本環境庁式）という名称で掲載されることになった自動車からのNOxの解析的拡散式を開発することでした。

　この拡散式は,通常のガウス型拡散式にパスキル・ギフォード(Pasquill-Gifford)流のシアーを持たせたものとは異なり,水平風速と拡散係数が地面からの鉛直高さｚのベキ乗に比例すると仮定した,より現実に近いパラメータに従う移流拡散方程式の定常解を基準にした模型(モデル)です。

　まず,こうした定常解をコンピュータの数値計算に依るのではなく,頭でつまり紙とペンで解析的に解くというのが配属されて最初の仕事内容でした。

　自動車という線煙源が高架ではなく,高さゼロの道路上にあってｙ軸に沿って無限大の長さを持っており,それに直角に風が吹く場合については,線煙源からの距離ｘと高さ:ｚの関数としてロバーツ(Roberts)の式という解があることは既にわかっていました。

　そして,私の使命は,まずは直角風ではなく無風（calm）のときと無限線煙源に平行な風(平行風)が吹く場合の解を求めることでした。

　特に高架道路のように,必ずしも線煙源の高さがゼロでなくて,有効高さHeをを持つとして,一般化して解きました。結果的に,直角風の場合のロバーツの式も有効高さＨｅがゼロでない場合に拡張することに成功しました。

　まずは,境界条件に合う無風時の解を求めることに挑戦しました。

　原理的には点煙源に対する３次元の解がわかれば,それを線煙源のあるｙ軸に沿って積分することで線煙源に対する２次元の解が得られます。

　そこで,まず3次元の無風の場合の偏微分方程式を解くことから始めました。そのため,解を点煙源からの距離 ｒ と高さ ｚ について変数分離しました。

　結局,変数分離で得られるそれぞれの常微分方程式の解は ｒ については 0 次の第２種変形ベッセル関数(modified Bessel's functio)になり,ｚ についても第１種変形ベッセル関数になることがわかりました。

　そこで,後は発生源条件に合うようにベッセル展開をする,つまり展開係数を求めれば解は得られるわけです。

　実を言うと,ここから最終解に到達するのに,岩波の「数学公式集」と２ヶ月以上にらめっこして,やっと運良く境界条件,発生源条件に合致する解を発見することができました。

　ここからは,点煙源の位置;ｙで－∞から＋∞まで積分すればいいわけですが結果は超幾何関数になります。しかし,特に有効高さＨｅがゼロのときには初等関数に帰着することがわかりました。

　そして,２次元無限線煙源の場合には平行風に対する方程式も無風時のときの方程式と同一なので,平行風の解は無風時と同じということで,これは即座に解決しました。

　しかし,風がゼロではなく有風時には,実はｚのベキ乗に特殊な条件がある場合について,既に３次元の解(Ｙｅｈ－Ｈｕａｎｇの解)が存在している,ことが後にわかり,直角風も平行風も共に単にこの既存の解を煙源に沿って積分することで２次元の解が求まることがわかりました。

　これの積分結果から,直角風の場合の解はＨｅがゼロのときはロバーツの式に一致し,Heがゼロでないときには第１種変形ベッセル関数となることもわかりました。

  こうして無限線源の解は全て求まりましたが,現実の道路は曲がりくねっており,ある軸に沿って無限の長さで一直線に伸びているわけではありません。

　しかし,これら２次元の解を求める過程において３次元解が既知となったので,有限長さの場合の解は,ｙについて－∞から＋∞でなく有限区間で積分すれば得られるはずです。

　幸いにして,直角風のときはｙの有限区間で積分しても単に無限線煙源の解にガウスの誤差関数がかかるだけという結果となり,平行風の場合も不完全γ関数,無風の場合も逆正接関数がかかるだけという式で近似できることがわかりました。

　平行風の場合の不完全γ関数をガウス誤差関数で近似することにし,これで所期の定常拡散方程式の解,または近似解の形は全て初等関数やそれに順ずる解析的関数という形で得られました。

　実際には風が線煙源に直角や平行でなく一般の風向であっても,３次元の点煙源解がわかっているので,積分により解を求めることは可能で,計算してみると風向角θを含む第２種変形ベッセル関数となって解は求まります。

　しかし,実際に法律として運用するモデルとしては煩雑なので,有風時（風速１ｍ/ｓ以上）については直角風解と平行風解のみで対応することにしました。

　また,３次元の定常拡散方程式はｘ,ｙ,ｚのうちを ｙによる微分項を時間 t によるそれに置き換え,その拡散係数であったｚのベキ乗を ｚ のゼロ乗,つまり定数に取ることで,ｘ.,ｚの２次元の非定常方程式になるので,元の３次元定常解が２次元の非定常解として得られるわけです。

　これらは任意風向の場合の解と同じく,敢えて発表しませんでした。

　そして風向と線煙源との鋭角θが40度を境として直角風であるか,平行風であるかを判定することにして,両者の式のどちらかを適用することにしました。また,ベッセル関数の計算の煩雑さを避けて高さHeがゼロのときの適用に限ることにして,このモデル式を以ってJEA式 としました。

　実際のパラメータ類を決める作業は,大阪府が府内の地形環境が異なる各地で各種の気象条件下において,約200ｍのパイプの多くの排出口からトレーサーガス（SF6：6フッ化イオウ）を撒いて,それを採取する実験で得た種々の濃度データを利用しました。

　当時の大型コンピュータによって,地形環境,および気象条件ごとにトレーサーガス濃度観測の実測値と先のJEA式による計算値との差の二乗和が最小となるように,多次元ニュートン法による非線形最小二乗法でJEA式の幾つかの未知パラメータの最適値を計算で決めました。

　そして,各地形環境ごとに,得られた式パラメータと気象条件をグラフにプロットし比較して,地形ごとに気象とパラメータの関係を求めることで最終的なモデルが完成したわけです。

  なお,後に六本木や西新宿での高架道路のデータをもとにHeがゼロでない場合の変形ベッセル関数を含むJEA式を修正して「東京都修正モデル」を作ったという記憶もあります。

　実は,当時の資料は今はどこにあるか見つけられなかったのですが,中国から日本の環境アセスメントを知るための研修にきた中国人の役人と大学(専門学校？)の先生の２人に会社の会議室で説明したときに作った英文の資料が見つかったので当時の記憶がかなり鮮明によみがえったわけです。

  うーん,でもこれは企業秘密の部類の一部で一種の自慢話になるかなあ,まあ,もう時効だろうし,ブログだから備忘録として残すのもいいかな。。。

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

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