移流拡散方程式を解く

　特殊な条件の下で移流拡散方程式を解いてみたいと思います。

　これは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(ｔ－ｔ₀)}]

です。

　Gauss積分を実行することにより,

　右辺＝[π³/{Ｋ_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にも仮想の鏡像点源があるとすれば∂Ｃ/∂ｚ＝0 を満たす解となります。

濃度は,Ｃ(ｒ) ～ [Ｑ/(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²)]] で与えられる,という簡易式を多用するようです。

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

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

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

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

    物理学