JP3746712B2

JP3746712B2 - 拡散物質の拡散状況予測方法

Info

Publication number: JP3746712B2
Application number: JP2001393179A
Authority: JP
Inventors: 良二大場; 清一工藤
Original assignee: Mitsubishi Heavy Industries Ltd
Current assignee: Mitsubishi Heavy Industries Ltd
Priority date: 2001-12-26
Filing date: 2001-12-26
Publication date: 2006-02-15
Anticipated expiration: 2021-12-26
Also published as: TW200424905A; TWI227434B; JP2003196574A

Description

【０００１】
【発明の属する技術分野】
本発明は拡散物質の拡散状況予測方法に関するものである。本発明は、拡散源（例えば放射性物質使用施設や煙突）から大気中に排出された物質（例えば放射性物質や煙）が、大気中にどのように拡散していくかを予測して、各地点で時々刻々変化していく物質の濃度を予測するようにしたものである。
【０００２】
【従来の技術】
放射性物質を扱う施設から、事故により放射性物質が外部に排出された場合には、放射性物質の拡散範囲や各地点での放射性物質の濃度を予測し、放射性物質による危険を受ける恐れがある地域を予測する拡散状況予測方法が開発されつつある。
【０００３】
この拡散状況予測方法は、放射性物質の拡散状況を予測する場合のみならず、例えば工場の煙突から排出されたガス体（煙）が大気中を拡散した場合において、各地点におけるガス体濃度を計算する場合や、環境アセスメントの解析における、拡散物質の拡散状況を解析する場合にも適用することができる。
【０００４】
大気中に排出された物質の拡散状況を、演算により予測するには、次の２つの演算をする必要がある。
（１）気体状況予測演算
（２）拡散状況予測演算
【０００５】
上記（１）の気体状況予測演算とは、気象ＧＰＶ（Grid Point Value）データやＡＭＥＤＡＳ等の気象観測データを基にして、大気現象を解析する偏微分方程式を演算することにより、事象発生（例えば放射性物質の外部排出）時点から所定時間先の時点まで、一定時間刻み毎の時点における、多数の評価地点（格子点位置）の風向・風速を演算により求める、つまり、一定時間刻み毎の風速場データを表す気体状況を求める演算をいう。
【０００６】
また、上記（２）の拡散状況予測演算とは、放出された拡散物質の濃度や性状ならびに前記風速場データを、物質（粒子）の拡散状態を演算する拡散方程式に代入することにより、各時間刻み毎の各格子点位置における拡散物質の濃度を求める演算をいう。
【０００７】
＜気体状況予測演算の説明＞
まず、気体状況予測演算の概略を説明する。気象観測データ、例えば気象ＧＰＶデータは、気象業務支援センターから１２時間ごとに配信される。この気象ＧＰＶデータは、地球の表面を南北方向に沿い伸びると共に東西方向の相互の離間距離が規定距離（２Ｋｍ）となっている複数の緯度仮想線と、地球の表面を東西方向に沿い伸びると共に南北方向の相互の離間距離が規定距離（２Ｋｍ）となっている複数の経度仮想線とが交差する地点（これを「親格子点位置」と称する）における、気象データ（風速ベクトル（風向，風速），気圧，温度，水分量）を示すものである。しかも、気象ＧＰＶデータは、各親格子点位置の気象データとして、配信時点，配信時点から３時間先、６時間先，９時間先，というように３時間間隔の５１時間分のデータが一括して配信される。
【０００８】
上述した気象ＧＰＶデータの親格子点位置の気象データは、空間的には親格子点位置の相互間距離が２Ｋｍと広く、しかも、時間的には３時間間隔と長いため、この親格子点位置の気象データにより示される気体状況（風向，風速）データ即ち風速場データのみでは、拡散物質の拡散濃度を演算することはできない。
【０００９】
このため、空間的に粗く、且つ、時間的にも粗い気象観測データから、空間的にも時間的にも密な気体状況（風向，風速等）を、大気現象を解析する偏微分方程式を演算することにより求める必要がある。
【００１０】
そこで、計算すべき計算領域（地球の表面のなかで予め設定した特定領域）に設定された親格子点位置の間に、子格子点位置を設定する。子格子点位置は、地球の表面を南北方向に沿い伸びると共に東西方向の相互の離間距離が一定距離（５０ｍ）となっている複数の緯度仮想線と、地球の表面を東西方向に沿い伸びると共に南北方向の相互の離間距離が一定距離（５０ｍ）となっている複数の経度仮想線とが交差する地点に配置されている。
【００１１】
そして、演算開始から一定時間刻み毎（例えば２０秒間隔毎）の子格子点位置及び親格子点位置の気象データを、大気現象を解析する偏微分方程式を差分解析演算することにより求める。大気現象を解析する偏微分方程式としては、コロラド州立大学とMission Research社で開発されたＲＡＭＳ（Regional Atmospheric Modeling System）コードで示されている、風速場解析の基本方程式を用いることができる。
【００１２】
このＲＡＭＳコードで示されている風速場解析の基本方程式は、運動方程式，熱エネルギ方程式，水分の拡散方程式及び連続の式からなり、次のような式（１）〜（６）で表される。
【００１３】
【数１】

【００１４】
このようにＲＡＭＳ（Regional Atmospheric Modeling System）コードで示されている風速場解析の基本方程式を演算して、演算開始から一定時間刻み毎（例えば２０秒間隔毎）の、各親格子点位置における気象データと、各子格子点位置における気象データを示す風向ベクトルデータ（風速場データ）が得られる。
【００１５】
＜拡散状況予測演算の概要説明＞
次に拡散状況予測演算について説明する。拡散状況予測演算をするには、コロラド州立大学とMission Research社で開発されたＨＹＰＡＣＴ（Hybrid Particle Concentration Transport Model ）コードに、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた２０秒刻み毎の各親格子点位置及び各子格子点位置の風速場データを次々に代入して、拡散状況の予測演算をする。拡散状況の予測演算の具体例としては、Lagrangian粒子拡散モデルを採用している。
【００１６】
このLagrangian粒子拡散モデルでは、次に示す式（７）〜（９）を用いて粒子の拡散速度（ｕ’，ｖ’，ｗ’）を計算し、各粒子を移動させる。
【００１７】
【数２】

【００１８】
ここで、ＨＹＰＡＣＴ（Hybrid Particle Concentration Transport Model ）コードに、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた２０秒刻み毎の各親格子点位置及び各子格子点位置の風速場データを次々に代入して、拡散状況の予測演算をした具体例を説明する。
【００１９】
この演算をするためには、排出源から大気中に排出される物質を多数の粒子Ｐに置換し、排出源の位置から演算周期Δｔ（ここではΔｔ＝２０秒）毎にＮ個（ここでは２０個）の粒子Ｐが発生すると設定する。
【００２０】
つまり、演算開始時点で２０個の粒子Ｐを発生させ、演算開始時点から２０秒後に２０個の粒子を発生させ、演算開始時点から４０秒後に２０個の粒子を発生させるというように、演算周期Δｔ（２０秒）毎に２０個の粒子を発生させていく。そして演算周期Δｔ（２０秒）毎に、各粒子Ｐの位置（空間座標）を演算により求める。
【００２１】
なお、演算開始時点（時刻０秒）において発生させた２０個の粒子Ｐを、Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰として示し、
演算開始時点から２０秒後において発生させた２０個の粒子Ｐを、Ｐ₂₀ ⁰¹，Ｐ₂₀ ⁰²，Ｐ₂₀ ⁰³，Ｐ₂₀ ⁰⁴，Ｐ₂₀ ⁰⁵，Ｐ₂₀ ⁰⁶，Ｐ₂₀ ⁰⁷，Ｐ₂₀ ⁰⁸，Ｐ₂₀ ⁰⁹，Ｐ₂₀ ¹⁰，Ｐ₂₀ ¹¹，Ｐ₂₀ ¹²，Ｐ₂₀ ¹³，Ｐ₂₀ ¹⁴，Ｐ₂₀ ¹⁵，Ｐ₂₀ ¹⁶，Ｐ₂₀ ¹⁷，Ｐ₂₀ ¹⁸，Ｐ₂₀ ¹⁹，Ｐ₂₀ ²⁰として示し、
演算開始時点から４０秒後において発生させた２０個の粒子Ｐを、Ｐ₄₀ ⁰¹，Ｐ₄₀ ⁰²，Ｐ₄₀ ⁰³，Ｐ₄₀ ⁰⁴，Ｐ₄₀ ⁰⁵，Ｐ₄₀ ⁰⁶，Ｐ₄₀ ⁰⁷，Ｐ₄₀ ⁰⁸，Ｐ₄₀ ⁰⁹，Ｐ₄₀ ¹⁰，Ｐ₄₀ ¹¹，Ｐ₄₀ ¹²，Ｐ₄₀ ¹³，Ｐ₄₀ ¹⁴，Ｐ₄₀ ¹⁵，Ｐ₄₀ ¹⁶，Ｐ₄₀ ¹⁷，Ｐ₄₀ ¹⁸，Ｐ₄₀ ¹⁹，Ｐ₄₀ ²⁰として示す。
つまり、符号「Ｐ」の後の下段に示した数字が、演算開始時点からの時間であり、符号「Ｐ」の後の上段に示した数字が、その時点において発生させた２０個の粒子を区別するものである。他の時点において発生させた粒子も同様に表記する。
【００２２】
まず、演算開始時点においては、排出源Ｓから２０個の粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰が発生する。
【００２３】
演算開始時点から２０秒後においては、図１９に示す排出源Ｓから新たに２０個の粒子Ｐ₂₀ ⁰¹，Ｐ₂₀ ⁰²，Ｐ₂₀ ⁰³，Ｐ₂₀ ⁰⁴，Ｐ₂₀ ⁰⁵，Ｐ₂₀ ⁰⁶，Ｐ₂₀ ⁰⁷，Ｐ₂₀ ⁰⁸，Ｐ₂₀ ⁰⁹，Ｐ₂₀ ¹⁰，Ｐ₂₀ ¹¹，Ｐ₂₀ ¹²，Ｐ₂₀ ¹³，Ｐ₂₀ ¹⁴，Ｐ₂₀ ¹⁵，Ｐ₂₀ ¹⁶，Ｐ₂₀ ¹⁷，Ｐ₂₀ ¹⁸，Ｐ₂₀ ¹⁹，Ｐ₂₀ ²⁰が発生する。
このとき、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰は、排出源Ｓから離れた位置にまで達すると共に拡散している。
各粒子Ｐの位置は、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた２０秒刻み毎の風速場データを用いて、Lagrangian粒子拡散モデルにおける各粒子Ｐの拡散速度（ｕ’，ｖ’，ｗ’）を計算し、各粒子を移動させることにより求める。
【００２４】
演算開始時点から４０秒後においては、図２０に示す排出源Ｓから新たに２０個の粒子Ｐ₄₀ ⁰¹，Ｐ₄₀ ⁰²，Ｐ₄₀ ⁰³，Ｐ₄₀ ⁰⁴，Ｐ₄₀ ⁰⁵，Ｐ₄₀ ⁰⁶，Ｐ₄₀ ⁰⁷，Ｐ₄₀ ⁰⁸，Ｐ₄₀ ⁰⁹，Ｐ₄₀ ¹⁰，Ｐ₄₀ ¹¹，Ｐ₄₀ ¹²，Ｐ₄₀ ¹³，Ｐ₄₀ ¹⁴，Ｐ₄₀ ¹⁵，Ｐ₄₀ ¹⁶，Ｐ₄₀ ¹⁷，Ｐ₄₀ ¹⁸，Ｐ₄₀ ¹⁹，Ｐ₄₀ ²⁰が発生する。
このとき、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰は、排出源Ｓから更に離れた位置にまで達すると共に更に拡散している。
また、演算開始時点から２０秒後において発生した２０個の粒子Ｐ₂₀ ⁰¹，Ｐ₂₀ ⁰²，Ｐ₂₀ ⁰³，Ｐ₂₀ ⁰⁴，Ｐ₂₀ ⁰⁵，Ｐ₂₀ ⁰⁶，Ｐ₂₀ ⁰⁷，Ｐ₂₀ ⁰⁸，Ｐ₂₀ ⁰⁹，Ｐ₂₀ ¹⁰，Ｐ₂₀ ¹¹，Ｐ₂₀ ¹²，Ｐ₂₀ ¹³，Ｐ₂₀ ¹⁴，Ｐ₂₀ ¹⁵，Ｐ₂₀ ¹⁶，Ｐ₂₀ ¹⁷，Ｐ₂₀ ¹⁸，Ｐ₂₀ ¹⁹，Ｐ₂₀ ²⁰は、排出源Ｓから離れた位置にまで達すると共に拡散している。
各粒子Ｐの位置は、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた２０秒刻み毎の風速場データを用いて、Lagrangian粒子拡散モデルにおける各粒子Ｐの拡散速度（ｕ’，ｖ’，ｗ’）を計算し、各粒子を移動させることにより求める。
【００２５】
演算開始時点から６０秒後においては、図２１に示す排出源Ｓから新たに２０個の粒子Ｐ₆₀ ⁰¹，Ｐ₆₀ ⁰²，Ｐ₆₀ ⁰³，Ｐ₆₀ ⁰⁴，Ｐ₆₀ ⁰⁵，Ｐ₆₀ ⁰⁶，Ｐ₆₀ ⁰⁷，Ｐ₆₀ ⁰⁸，Ｐ₆₀ ⁰⁹，Ｐ₆₀ ¹⁰，Ｐ₆₀ ¹¹，Ｐ₆₀ ¹²，Ｐ₆₀ ¹³，Ｐ₆₀ ¹⁴，Ｐ₆₀ ¹⁵，Ｐ₆₀ ¹⁶，Ｐ₆₀ ¹⁷，Ｐ₆₀ ¹⁸，Ｐ₆₀ ¹⁹，Ｐ₆₀ ²⁰が発生する。
このとき、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰は、排出源Ｓから更に離れた位置にまで達すると共に更に拡散している。
また、演算開始時点から２０秒後において発生した２０個の粒子Ｐ₂₀ ⁰¹，Ｐ₂₀ ⁰²，Ｐ₂₀ ⁰³，Ｐ₂₀ ⁰⁴，Ｐ₂₀ ⁰⁵，Ｐ₂₀ ⁰⁶，Ｐ₂₀ ⁰⁷，Ｐ₂₀ ⁰⁸，Ｐ₂₀ ⁰⁹，Ｐ₂₀ ¹⁰，Ｐ₂₀ ¹¹，Ｐ₂₀ ¹²，Ｐ₂₀ ¹³，Ｐ₂₀ ¹⁴，Ｐ₂₀ ¹⁵，Ｐ₂₀ ¹⁶，Ｐ₂₀ ¹⁷，Ｐ₂₀ ¹⁸，Ｐ₂₀ ¹⁹，Ｐ₂₀ ²⁰は、排出源Ｓから更に離れた位置にまで達すると共に更に拡散している。
また、演算開始時点から４０秒後において発生した２０個の粒子Ｐ₄₀ ⁰¹，Ｐ₄₀ ⁰²，Ｐ₄₀ ⁰³，Ｐ₄₀ ⁰⁴，Ｐ₄₀ ⁰⁵，Ｐ₄₀ ⁰⁶，Ｐ₄₀ ⁰⁷，Ｐ₄₀ ⁰⁸，Ｐ₄₀ ⁰⁹，Ｐ₄₀ ¹⁰，Ｐ₄₀ ¹¹，Ｐ₄₀ ¹²，Ｐ₄₀ ¹³，Ｐ₄₀ ¹⁴，Ｐ₄₀ ¹⁵，Ｐ₄₀ ¹⁶，Ｐ₄₀ ¹⁷，Ｐ₄₀ ¹⁸，Ｐ₄₀ ¹⁹，Ｐ₄₀ ²⁰は、排出源Ｓから離れた位置に達すると共に拡散している。
各粒子Ｐの位置は、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた２０秒刻み毎の風速場データを用いて、Lagrangian粒子拡散モデルにおける各粒子Ｐの拡散速度（ｕ’，ｖ’，ｗ’）を計算し、各粒子を移動させることにより求める。
【００２６】
上述したように、演算周期Δｔ（２０秒）毎に２０個の粒子を次々と発生させていくと共に、各演算周期Δｔ（２０秒）毎における粒子の位置つまり空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））を求めていく。
【００２７】
そして、演算開始から所定時間経過したときにおいて、排出源Ｓから所定距離離れた単位空間（予測地域の単位体積）に、図２２に示すように、粒子Ｐが存在していた場合、この粒子の数から、この単位空間における物質の濃度を計算することができる。
【００２８】
即ち、排出源Ｓにおいて、１秒間にＱ（ｍ³）の物質が排出されているとすると、粒子Ｐは２０秒間に２０個（換算すると１秒間に１個）発生しているので、各粒子Ｐは、１個につきＱ／１（ｍ³）の排出源強度を有していることになる。そこで、この単位空間に存在する粒子Ｐの数に排出源強度Ｑ／１（ｍ³）を掛けることにより、この単位空間における物質の濃度を求めることができる。
【００２９】
上述した具体例を一般的に示すと次のようになる。排出源から排出されるガス体などの物質を多数の粒子で置換する。そして、排出源から毎秒Ｎ個の粒子を放出する。この場合、計算上での粒子の排出量はＮ／sec である。実際の排出源から排出される物質の排出量がＱ（ｍ³／sec ）である場合、各粒子はＱ／Ｎ（ｍ³）の排出源強度を有することになる。
【００３０】
各粒子ごとに運動方程式を非定常に数値計算することによって、即ち、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた風速場データを、粒子の運動方程式であるＨＹＰＡＣＴ（Hybrid Particle Concentration Transport Model ）コードに代入し、Lagrangian粒子拡散モデルを用いて各粒子Ｐの拡散速度（ｕ’，ｖ’，ｗ’）を計算し、各粒子を移動させることによって、各粒子の座標を非定常的に決定することができる。つまり、各粒子の空間座標を、演算周期Δｔ毎に決定することができる。なお、Lagrangian粒子モデルにより求めてデータ記録装置に記録される各粒子のデータは、各粒子の空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））だけである。
【００３１】
粒子（物質）の運動方程式であるＨＹＰＡＣＴコードは、粒子の移流、拡散、重力沈降現象を表現するものである。ここにおいて、粒子の移流現象は、大気の時間平均速度に依存し、拡散現象は、大気の乱流速度に依存し、重力沈降は、粒子の質量、重力加速度、空気の粘性係数などに依存する（図２３参照）
【００３２】
空気中の単位体積中の粒子個数がｎ個である場合、この空間中のガス濃度（物質濃度）はｎ×Ｑ／Ｎ（ガスｍ³／空気ｍ³）となる。つまり、この単位空間に存在する粒子数ｎに、各粒子が有している排出源強度Ｑ／Ｎを掛けたものとなる。
【００３３】
【発明が解決しようとする課題】
この環境濃度（単位体積における物質濃度）は、排出される物質の排出量の時間変化に依存する。このため、排出量が時間と共に変化する条件では、拡散計算は、各排出条件毎に実施する必要がある。したがって、想定される排出条件が多い場合には、排出量ケース分の拡散計算を実施する必要があり、結果として、膨大な計算時間が必要である。
【００３４】
即ち、図２４に示すように、例えば排出源Ｓ（例えば煙突）からガス（物質）が排出されているときに、風下のある地点Ｆにおけるガス濃度の時間変化は、排出源Ｓから排出される物質の時間変化に応じて変化する。
【００３５】
つまり、図２５（ａ）のように物質の排出量が時間と共に変化する場合には、地点Ｆにおける物質の濃度は、図２５（ｂ）のように時間と共に変化し、図２６（ａ）のように物質の排出量が一定の場合には、地点Ｆにおける物質の濃度は、図２６（ｂ）のように一定値にまで上昇してから一定濃度を維持し、図２７（ａ）のように物質が瞬間的に排出される場合には、地点Ｆにおける物質の濃度は、図２７（ｂ）のように一時的に上昇してから零になる。
【００３６】
このように、物質の排出量が時間と共に変化する場合には、粒子の発生個数を物質の排出量に合わせて時間と共に変化させる必要がある。そして、このように時間の経過と共に発生個数を変化させた粒子の移動位置を求め、この粒子の移動位置から物質の濃度計算をすることになる。したがって、排出量の変化が異なる各ケースごとに、拡散計算をしておかなければならず、膨大な計算結果が必要になるのである。
【００３７】
例えば放射性物質を扱う施設にて放射性物質が外部に排出された事故が発生した場合には、極めて多数の物質（例えば１００種類程度の物質）が排出される。しかも、各物質毎にその排出量が時間に応じてそれぞれ異なっている。したがって、各物質毎に、粒子の発生個数を物質の排出量に合わせて時間と共に変化させ、このように発生個数を変化させた粒子の移動位置を求め、この粒子の移動位置から物質の濃度計算をすることになる。したがって、この場合には、例えば１００種類の物質に対応した１００種類の拡散計算をしておく必要がある。
【００３８】
本発明は、上記従来技術に鑑み、多種類の物質が排出されると共に、各物質の排出量が時間と共に変化する場合であっても、物質の拡散状況を短時間で予測演算することができる拡散物質の拡散状況予測方法及び拡散物質の拡散状況予測システムを提供することを目的とする。
【００３９】
【課題を解決するための手段】
上記課題を解決する本発明の拡散物質の拡散状況予測方法は、拡散物質の拡散状況をコンピュータによる演算処理により求める拡散物質の拡散状況予測方法であって、
前記コンピュータは、
排出源から大気中に排出された物質が大気中を拡散していく状況を予測するため、前記物質を多数の粒子に置換して、排出源の位置から演算周期毎に予め設定した個数の粒子が発生すると設定すると共に、
排出源の位置を含む領域内の多数の地点における、時間の経過に沿って変化する風向・風速を示す風速場データを、粒子の拡散状態を演算する拡散方程式に代入することにより、各粒子の拡散速度を求め、この拡散速度から各粒子が存在する空間位置を示す空間座標を各演算周期毎に求めると共に、前記粒子を最初に発生させた時点からの経過時間である排出後経過時間を計測し、各演算周期における各粒子の空間座標と各粒子の排出後経過時間を対応させてデータ記録装置に記録しておき、
排出された物質の排出後経過時間の時間経過に伴う排出量の変化に比例させて、排出後経過時間の時間経過に沿う粒子に対する排出源強度データを設定しておき、
前記データ記録装置に記録していた、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間を読み出すと共に、読み出した排出後経過時間を参照して各粒子が発生した時点を求め、この時点における各粒子の排出源強度を前記排出源強度データから求め、前記データ記録装置に、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と排出源強度を対応させて再記録させ、
所定の演算周期における所定の領域の前記物質の濃度は、当該所定の演算周期における当該所定の領域に存在する全ての粒子の排出源強度を積算することにより求めることを特徴とする。
【００４０】
また本発明の拡散物質の拡散状況予測方法は、拡散物質の拡散状況をコンピュータによる演算処理により求める拡散物質の拡散状況予測方法であって、
前記コンピュータは、
複数の排出源から大気中に排出された物質が大気中を拡散していく状況を予測するため、前記物質を多数の粒子に置換して、各排出源の位置から演算周期毎に予め設定した個数の粒子がそれぞれ発生すると設定すると共に、
排出源の位置を含む領域内の多数の地点における、時間の経過に沿って変化する風向・風速を示す風速場データを、粒子の拡散状態を演算する拡散方程式に代入することにより、各粒子の拡散速度を求め、この拡散速度から各粒子が存在する空間位置を示す空間座標を各演算周期毎に求めると共に、前記粒子を最初に発生させた時点からの経過時間である排出後経過時間を計測し、各演算周期における各粒子の空間座標と各粒子の排出後経過時間と排出源を識別する排出源識別情報とを対応させてデータ記録装置に記録しておき、
各排出源から排出された物質の排出後経過時間の時間経過に伴う排出量の変化に比例させて、排出後経過時間の時間経過に沿う粒子に対する排出源強度データを各排出源毎にそれぞれ設定しておき、
前記データ記録装置に記録していた、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と各粒子の排出源識別情報を読み出すと共に、読み出した排出後経過時間を参照して各粒子が発生した時点を求め、読み出した排出源識別情報を参照して粒子が発生した時点における各粒子の排出源強度をその粒子が発生した排出源に対応した前記排出源強度データから求め、前記データ記録装置に、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と排出源強度を対応させて再記録させ、
所定の演算周期における所定の領域の前記物質の濃度は、当該所定の演算周期における当該所定の領域に存在する全ての粒子の排出源強度を積算することにより求めることを特徴とする。
【００４１】
また本発明の拡散物質の拡散状況予測方法では、
前記排出源強度データは、前記排出源から実際に排出される物質の濃度を実測することにより求めて設定されたり、
前記排出源強度データは、前記排出源の周囲の観測点で実測した物質の濃度の時間変化を基に設定されることを特徴とする。
【００４３】
【発明の実施の形態】
以下に、本発明の実施の形態を図面に基づき詳細に説明する。
【００４４】
＜第１の実施の形態（排出源が１つの場合）＞
本発明の第１の実施の形態にかかる拡散物質の拡散状況予測方法を、図１〜図１０を参照しつつ説明する。
【００４５】
第１の実施の形態での第１ステップ（計算フロー図である図１参照）では次のような処理をする。即ち、排出源Ｓ（粒子の拡散状況を示す図２〜図４参照）から実際に排出される物質の排出量が一定であっても、時間経過に沿い排出量が変化する場合であっても、まずは物質の排出量Ｑ（ｍ³／sec ）を一定値（＝１．０）として、従来のLagrangian粒子拡散モデルを使用して、粒子の挙動を数値計算する。更に、各粒子が有する情報である空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））の他に、粒子を最初に発生させた時点からの経過時間である排出後経過時間Ｔｉ（ｔ）を、各演算周期Δｔ毎に、データ記録装置１に記録する。
【００４６】
この第１ステップの処理を、具体的に説明すると次のようになる。この演算では、演算周期Δｔ（ここではΔｔ＝２０秒）毎に、２０個の粒子を発生させると共に、演算周期Δｔ（２０秒）毎に粒子Ｐの位置（空間座標）を演算する。
【００４７】
まず、演算開始時点においては、排出源Ｓから２０個の粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰が発生する。
【００４８】
演算開始時点から２０秒後においては、図２に示す排出源Ｓから新たに２０個の粒子Ｐ₂₀ ⁰¹，Ｐ₂₀ ⁰²，Ｐ₂₀ ⁰³，Ｐ₂₀ ⁰⁴，Ｐ₂₀ ⁰⁵，Ｐ₂₀ ⁰⁶，Ｐ₂₀ ⁰⁷，Ｐ₂₀ ⁰⁸，Ｐ₂₀ ⁰⁹，Ｐ₂₀ ¹⁰，Ｐ₂₀ ¹¹，Ｐ₂₀ ¹²，Ｐ₂₀ ¹³，Ｐ₂₀ ¹⁴，Ｐ₂₀ ¹⁵，Ｐ₂₀ ¹⁶，Ｐ₂₀ ¹⁷，Ｐ₂₀ ¹⁸，Ｐ₂₀ ¹⁹，Ｐ₂₀ ²⁰が発生する。
このとき、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰は、排出源Ｓから離れた位置にまで達すると共に拡散している。
各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の位置は、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた２０秒刻み毎の風速場データを用いて、Lagrangian粒子拡散モデルにおける各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の拡散速度（ｕ’，ｖ’，ｗ’）を計算し、各粒子を移動させることにより求める。
【００４９】
更に、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰は、演算開始時点（粒子を最初に発生させた時点）から２０秒が経過している。そこで、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝２０），ｙｉ（ｔ＝２０），ｚｉ（ｔ＝２０））に、排出後経過時間Ｔｉ（ｔ）＝２０秒をそれぞれ対応させてデータ記録装置１に記録する（図１，図２参照）。
【００５０】
演算開始時点から４０秒後においては、図３に示す排出源Ｓから新たに２０個の粒子Ｐ₄₀ ⁰¹，Ｐ₄₀ ⁰²，Ｐ₄₀ ⁰³，Ｐ₄₀ ⁰⁴，Ｐ₄₀ ⁰⁵，Ｐ₄₀ ⁰⁶，Ｐ₄₀ ⁰⁷，Ｐ₄₀ ⁰⁸，Ｐ₄₀ ⁰⁹，Ｐ₄₀ ¹⁰，Ｐ₄₀ ¹¹，Ｐ₄₀ ¹²，Ｐ₄₀ ¹³，Ｐ₄₀ ¹⁴，Ｐ₄₀ ¹⁵，Ｐ₄₀ ¹⁶，Ｐ₄₀ ¹⁷，Ｐ₄₀ ¹⁸，Ｐ₄₀ ¹⁹，Ｐ₄₀ ²⁰が発生する。
このとき、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰は、排出源Ｓから更に離れた位置にまで達すると共に更に拡散している。
また、演算開始時点から２０秒後において発生した２０個の粒子Ｐ₂₀ ⁰¹，Ｐ₂₀ ⁰²，Ｐ₂₀ ⁰³，Ｐ₂₀ ⁰⁴，Ｐ₂₀ ⁰⁵，Ｐ₂₀ ⁰⁶，Ｐ₂₀ ⁰⁷，Ｐ₂₀ ⁰⁸，Ｐ₂₀ ⁰⁹，Ｐ₂₀ ¹⁰，Ｐ₂₀ ¹¹，Ｐ₂₀ ¹²，Ｐ₂₀ ¹³，Ｐ₂₀ ¹⁴，Ｐ₂₀ ¹⁵，Ｐ₂₀ ¹⁶，Ｐ₂₀ ¹⁷，Ｐ₂₀ ¹⁸，Ｐ₂₀ ¹⁹，Ｐ₂₀ ²⁰は、排出源Ｓから離れた位置にまで達すると共に拡散している。
各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰，Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の位置は、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた２０秒刻み毎の風速場データを用いて、Lagrangian粒子拡散モデルにおける各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰，Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の拡散速度（ｕ’，ｖ’，ｗ’）を計算し、各粒子を移動させることにより求める。
【００５１】
更に、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰は、演算開始時点（粒子を最初に発生させた時点）から４０秒が経過している。そこで、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝４０），ｙｉ（ｔ＝４０），ｚｉ（ｔ＝４０））に、排出後経過時間Ｔｉ（ｔ）＝４０秒をそれぞれ対応させてデータ記録装置１に記録する（図１，図３参照）。
また、演算開始時点から２０秒後に発生した粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰は、演算開始時点（粒子を最初に発生させた時点）から２０秒が経過している。そこで、各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の各空間座標（ｘｉ（ｔ＝４０），ｙｉ（ｔ＝４０），ｚｉ（ｔ＝４０））に、排出後経過時間Ｔｉ（ｔ）＝２０秒をそれぞれ対応させてデータ記録装置１に記録する（図１，図３参照）。
【００５２】
演算開始時点から６０秒後においては、図４に示す排出源Ｓから新たに２０個の粒子Ｐ₆₀ ⁰¹，Ｐ₆₀ ⁰²，Ｐ₆₀ ⁰³，Ｐ₆₀ ⁰⁴，Ｐ₆₀ ⁰⁵，Ｐ₆₀ ⁰⁶，Ｐ₆₀ ⁰⁷，Ｐ₆₀ ⁰⁸，Ｐ₆₀ ⁰⁹，Ｐ₆₀ ¹⁰，Ｐ₆₀ ¹¹，Ｐ₆₀ ¹²，Ｐ₆₀ ¹³，Ｐ₆₀ ¹⁴，Ｐ₆₀ ¹⁵，Ｐ₆₀ ¹⁶，Ｐ₆₀ ¹⁷，Ｐ₆₀ ¹⁸，Ｐ₆₀ ¹⁹，Ｐ₆₀ ²⁰が発生する。
このとき、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰²，Ｐ₀₀ ⁰³，Ｐ₀₀ ⁰⁴，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ⁰⁶，Ｐ₀₀ ⁰⁷，Ｐ₀₀ ⁰⁸，Ｐ₀₀ ⁰⁹，Ｐ₀₀ ¹⁰，Ｐ₀₀ ¹¹，Ｐ₀₀ ¹²，Ｐ₀₀ ¹³，Ｐ₀₀ ¹⁴，Ｐ₀₀ ¹⁵，Ｐ₀₀ ¹⁶，Ｐ₀₀ ¹⁷，Ｐ₀₀ ¹⁸，Ｐ₀₀ ¹⁹，Ｐ₀₀ ²⁰は、排出源Ｓから更に離れた位置にまで達すると共に更に拡散している。
また、演算開始時点から２０秒後において発生した２０個の粒子Ｐ₂₀ ⁰¹，Ｐ₂₀ ⁰²，Ｐ₂₀ ⁰³，Ｐ₂₀ ⁰⁴，Ｐ₂₀ ⁰⁵，Ｐ₂₀ ⁰⁶，Ｐ₂₀ ⁰⁷，Ｐ₂₀ ⁰⁸，Ｐ₂₀ ⁰⁹，Ｐ₂₀ ¹⁰，Ｐ₂₀ ¹¹，Ｐ₂₀ ¹²，Ｐ₂₀ ¹³，Ｐ₂₀ ¹⁴，Ｐ₂₀ ¹⁵，Ｐ₂₀ ¹⁶，Ｐ₂₀ ¹⁷，Ｐ₂₀ ¹⁸，Ｐ₂₀ ¹⁹，Ｐ₂₀ ²⁰は、排出源Ｓから更に離れた位置にまで達すると共に更に拡散している。
また、演算開始時点から４０秒後において発生した２０個の粒子Ｐ₄₀ ⁰¹，Ｐ₄₀ ⁰²，Ｐ₄₀ ⁰³，Ｐ₄₀ ⁰⁴，Ｐ₄₀ ⁰⁵，Ｐ₄₀ ⁰⁶，Ｐ₄₀ ⁰⁷，Ｐ₄₀ ⁰⁸，Ｐ₄₀ ⁰⁹，Ｐ₄₀ ¹⁰，Ｐ₄₀ ¹¹，Ｐ₄₀ ¹²，Ｐ₄₀ ¹³，Ｐ₄₀ ¹⁴，Ｐ₄₀ ¹⁵，Ｐ₄₀ ¹⁶，Ｐ₄₀ ¹⁷，Ｐ₄₀ ¹⁸，Ｐ₄₀ ¹⁹，Ｐ₄₀ ²⁰は、排出源Ｓから離れた位置に達すると共に拡散している。
各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰，Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰，Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰の位置は、ＲＡＭＳ（Regional Atmospheric Modeling System）コードにより求めた２０秒刻み毎の風速場データを用いて、Lagrangian粒子拡散モデルにおける各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰，Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰，Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰の拡散速度（ｕ’，ｖ’，ｗ’）を計算し、各粒子を移動させることにより求める。
【００５３】
更に、演算開始時点にて発生した粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰は、演算開始時点（粒子を最初に発生させた時点）から６０秒が経過している。そこで、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））に、排出後経過時間Ｔｉ（ｔ）＝６０秒をそれぞれ対応させてデータ記録装置１に記録する（図１，図４参照）。
また、演算開始時点から２０秒後に発生した粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰は、演算開始時点（粒子を最初に発生させた時点）から４０秒が経過している。そこで、各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））に、排出後経過時間Ｔｉ（ｔ）＝４０秒をそれぞれ対応させてデータ記録装置１に記録する（図１，図４参照）。
また、演算開始時点から４０秒後に発生した粒子Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰は、演算開始時点（粒子を最初に発生させた時点）から２０秒が経過している。そこで、各粒子Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））に、排出後経過時間Ｔｉ（ｔ）＝２０秒をそれぞれ対応させてデータ記録装置１に記録する（図１，図４参照）。
【００５４】
上述したように、演算周期Δｔ（２０秒）毎に２０個の粒子を次々と発生させていくと共に、各演算周期Δｔ（２０秒）毎における粒子の位置つまり空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））を求めていく。また、各演算周期Δｔにおける排出後経過時間Ｔｉ（ｔ）を計測しておき、各演算周期における各粒子の空間座標と各粒子の排出後経過時間を対応させて、データ記録装置１に次々と記録していく。
【００５５】
次に第２ステップ（図１参照）の処理を、具体的に説明する。前述した第１ステップでは、物質の排出量Ｑ（ｍ³／sec ）を一定値（＝１．０）として数値計算を進めていた。しかし、実際の排出源Ｓから排出される物質の排出量は、図５に示すように排出後経過時間Ｔｉ（ｔ）の経過と共に変化することが多い。そこで、このように排出量が時間と共に変化する場合には、この図５に示す物質の排出量変化に応じた、図６に示すような、排出後経過時間Ｔｉ（ｔ）の時間経過に沿う粒子に対する排出源強度を示すデータを設定する。
【００５６】
図６における排出源強度データでは、例えば、排出後経過時間Ｔｉ（ｔ）が０秒，２０秒，６０秒では、排出源強度はそれぞれ０．３，０．９，０．６となる。
【００５７】
次に、データ記録装置１に記録している、各演算周期毎の、各粒子の空間座標と各粒子の排出後経過時間Ｔｉ（ｔ）を読み出すと共に、各粒子毎にその排出後経過時間Ｔｉ（ｔ）を参照してその粒子が発生した時点を求め、この時点における各粒子の排出源強度を図６に示す排出源強度データから求める。更に、各演算周期毎に、各粒子の空間座標と各粒子の排出後経過時間Ｔｉ（ｔ）と排出源強度とを対応させて、データ記憶装置１に再記録する。
【００５８】
具体的に説明すると、排出後経過時間Ｔｉ（ｔ）が２０秒のとき（第１回目の演算周期）のときのデータとしては、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝２０），ｙｉ（ｔ＝２０），ｚｉ（ｔ＝２０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒が対応してデータ記録装置１に記録されている（図２参照）。
【００５９】
そこで、この各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝２０），ｙｉ（ｔ＝２０），ｚｉ（ｔ＝２０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒を読み出し、現在の時刻ｔ＝２０秒から排出後経過時間Ｔｉ（ｔ）＝２０秒を減算することにより、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰が発生した時点である排出後経過時間Ｔｉ（ｔ）＝０秒を求める。そして、図６に示す排出強度データから、粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰が発生したときの排出後経過時間Ｔｉ（ｔ）＝０秒のときの排出源強度０．３を求める。
【００６０】
そして、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝２０），ｙｉ（ｔ＝２０），ｚｉ（ｔ＝２０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒と、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の排出源強度０．３とを対応させて、データ記録装置１に再記録する。
【００６１】
また、排出後経過時間Ｔｉ（ｔ）が４０秒のとき（第２回目の演算周期）のときのデータとしては、
各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝４０），ｙｉ（ｔ＝４０），ｚｉ（ｔ＝４０））と、排出後経過時間Ｔｉ（ｔ）＝４０秒、並びに、
各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の各空間座標（ｘｉ（ｔ＝４０），ｙｉ（ｔ＝４０），ｚｉ（ｔ＝４０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒が対応してデータ記録装置１に記録されている（図３参照）。
【００６２】
そこで、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝４０），ｙｉ（ｔ＝４０），ｚｉ（ｔ＝４０））と、排出後経過時間Ｔｉ（ｔ）＝４０秒を読み出し、現在の時刻ｔ＝４０秒から排出後経過時間Ｔｉ（ｔ）＝４０秒を減算することにより、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰が発生した時点である排出後経過時間Ｔｉ（ｔ）＝０秒を求める。そして、図６に示す排出強度データから、粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰が発生したときの排出後経過時間Ｔｉ（ｔ）＝０秒のときの排出源強度０．３を求める。
同様に、各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の各空間座標（ｘｉ（ｔ＝４０），ｙｉ（ｔ＝４０），ｚｉ（ｔ＝４０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒を読み出し、現在の時刻ｔ＝４０秒から排出後経過時間Ｔｉ（ｔ）＝２０秒を減算することにより、各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰が発生した時点である排出後経過時間Ｔｉ（ｔ）＝２０秒を求める。そして、図６に示す排出強度データから、粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰が発生したときの排出後経過時間Ｔｉ（ｔ）＝２０秒のときの排出源強度０．５を求める。
【００６３】
そして、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝４０），ｙｉ（ｔ＝４０），ｚｉ（ｔ＝４０））と、排出後経過時間Ｔｉ（ｔ）＝４０秒と、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の排出源強度０．３とを対応させて、データ記録装置１に再記録する。
また各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の各空間座標（ｘｉ（ｔ＝４０），ｙｉ（ｔ＝４０），ｚｉ（ｔ＝４０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒と、各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の排出源強度０．５とを対応させて、データ記録装置１に再記録する。
【００６４】
また、排出後経過時間Ｔｉ（ｔ）が６０秒のとき（第３回目の演算周期）のときのデータとしては、
各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝６０秒、及び、
各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝４０秒、並びに、
各粒子Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒、が対応してデータ記録装置１に記録されている（図４参照）。
【００６５】
そこで、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝６０秒を読み出し、現在の時刻ｔ＝６０秒から排出後経過時間Ｔｉ（ｔ）＝６０秒を減算することにより、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰が発生した時点である排出後経過時間Ｔｉ（ｔ）＝０秒を求める。そして、図６に示す排出強度データから、粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰が発生したときの排出後経過時間Ｔｉ（ｔ）＝０秒のときの排出源強度０．３を求める。
同様に、各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝４０秒を読み出し、現在の時刻ｔ＝６０秒から排出後経過時間Ｔｉ（ｔ）＝４０秒を減算することにより、各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰が発生した時点である排出後経過時間Ｔｉ（ｔ）＝２０秒を求める。そして、図６に示す排出強度データから、粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰が発生したときの排出後経過時間Ｔｉ（ｔ）＝２０秒のときの排出源強度０．５を求める。
同様に、各粒子Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒を読み出し、現在の時刻ｔ＝６０秒から排出後経過時間Ｔｉ（ｔ）＝２０秒を減算することにより、各粒子Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰が発生した時点である排出後経過時間Ｔｉ（ｔ）＝４０秒を求める。そして、図６に示す排出強度データから、粒子Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰が発生したときの排出後経過時間Ｔｉ（ｔ）＝４０秒のときの排出源強度０．９を求める。
【００６６】
そして、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝６０秒と、各粒子Ｐ₀₀ ⁰¹〜Ｐ₀₀ ²⁰の排出源強度０．３とを対応させて、データ記録装置１に再記録する。
また各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝４０秒と、各粒子Ｐ₂₀ ⁰¹〜Ｐ₂₀ ²⁰の排出源強度０．５とを対応させて、データ記録装置１に再記録する。
また各粒子Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰の各空間座標（ｘｉ（ｔ＝６０），ｙｉ（ｔ＝６０），ｚｉ（ｔ＝６０））と、排出後経過時間Ｔｉ（ｔ）＝２０秒と、各粒子Ｐ₄₀ ⁰¹〜Ｐ₄₀ ²⁰の排出源強度０．９とを対応させて、データ記録装置１に再記録する。
【００６７】
以降の演算周期においても、同様な処理演算をして、各粒子の各空間座標と、排出後経過時間Ｔｉ（ｔ）と、各粒子の排出源強度を対応させて、データ記録装置に再記録する。
【００６８】
次に第３ステップ（図１参照）の処理を、具体的に説明する。例えば、排出後経過時間Ｔｉ（ｔ）＝１２０秒のときにおいて、図７に示すような排出源Ｓから所定距離離れた所定の格子領域（単位体積となっている単位空間）Ｉ，Ｊ，Ｋでの、物質の濃度を計算するには、データ記録装置１から、排出後経過時間Ｔｉ（ｔ）＝１２０秒における、この格子領域に存在する粒子を読み出す。読み出したところ、図７に示すように粒子が存在していた場合、これら各粒子が有している排出源強度を積算することにより、この単位空間における物質の濃度を計算することができる。
【００６９】
つまり、図７に示す格子領域には、
強度が０．３となっている４つの粒子Ｐ₀₀ ⁰¹，Ｐ₀₀ ⁰⁵，Ｐ₀₀ ¹⁰，Ｐ₀₀ ²⁰と、
強度が０．５となっている３つの粒子Ｐ₂₀ ⁰¹，Ｐ₂₀ ⁰⁷，Ｐ₂₀ ¹⁷と
強度が０．９となっている２つの粒子Ｐ₄₀ ⁰⁸，Ｐ₄₀ ¹⁰と、
強度が０．６となっている１つの粒子Ｐ₆₀ ¹⁷とが存在している。
このため、これら粒子の排出源強度を次のように積算することにより、この単位空間における物質の濃度を５．１と計算することができる。
（０．３×４）＋（０．５×３）＋（０．９×２）＋（０．６×１）＝５．１
【００７０】
前述した第１の実施の形態を一般的（数学的）に説明すると次のようになる。第１の実施の形態では、図８に示すように、排出源Ｓから物質（ガス体等）が排出されているときに、排出源Ｓの風下のある格子領域Ｉ，Ｊ，Ｋにおける物質濃度（ガス濃度）を、時間変化に応じて予測するものである。しかも、図９（ａ）に示すように物質の排出量Ｑが一定の場合において、図９（ｂ）に示すように格子領域Ｉ，Ｊ，Ｋの濃度時間変化を予測演算することができることは勿論、図１０（ａ）に示すように物質の排出量が時間変化する排出量ｑ（ｔ）であっても、図１０（ｂ）に示すように格子領域の濃度時間変化を予測演算することができるようにしたものである。
【００７１】
第１の実施の形態ではまず、排出源Ｓから実際に排出される物質の排出量が一定であっても、時間経過に沿い排出量が変化する場合であっても、まずは物質の排出量Ｑ（ｍ³／sec ）を一定値（＝１．０）として、従来のLagrangian粒子拡散モデルを使用して、物質を粒子に置換し、排出源Ｓから毎秒Ｎ個の粒子を発生させ、各粒子の挙動を数値計算して粒子の位置を示す空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））を求める。更に、各粒子が有する情報である空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））の他に、粒子を最初に発生させた時点からの経過時間である排出後経過時間Ｔｉ（ｔ）を、各演算周期Δｔ毎に、データ記録装置に記録する。これにより、従来のLagrangian粒子拡散モデルを利用して、時間変化するあらゆる排出量ｑ（ｔ）に対応した濃度分布の時間変化を計測することができる。
【００７２】
そのためには、時間変化する物質の排出量ｑ（ｔ）に比例した、排出後経過時間Ｔｉ（ｔ）の時間経過に沿う粒子に対する排出源強度を示すデータを設定する。そして、ある時刻（ｔ）において、各粒子の位置を示す空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））と排出後経過時間Ｔｉ（ｔ）をデータ記憶装置から読み出す。
【００７３】
一定排出量（Ｑ＝１．０）の場合には、各粒子の排出源強度はＱ／Ｎ（ｍ³）＝１／Ｎであるが、時間変化する排出量ｑ（ｔ）の場合には、各粒子の排出源強度はｑ（ｔ−Ｔｉ）／Ｎ（ｍ³）となる。
【００７４】
再度、各粒子が有する情報である空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））の他に、排出後経過時間Ｔｉ（ｔ）および各粒子の排出源強度ｑｉ（ｔ−Ｔｉ）／Ｎ（ｍ³）を、各時刻（ｔ）毎に、データ記録装置に再記録させる。
【００７５】
Lagrangian粒子モデルと、排出後経過時間と、時間変化する排出量ｑ（ｔ）に対応させた排出源強度ｑ（ｔ−Ｔｉ）／Ｎ（ｍ³）を利用した本実施例では、空間中の単位体積中（＝空気１ｍ³）の各粒子が有する排出源強度が異なるので、各粒子の排出源強度ｑｉ（ｔ−Ｔｉ）／Ｎ（ｍ³）を積算したΣｑｉ（ｔ−Ｔｉ）／Ｎ（ｍ³）がこの単位体積中に存在するガス量となる。したがって、この空間中ガス濃度はΣｑｉ（ｔ−Ｔｉ）／Ｎ（ｍ³）／Ｎ（ガスｍ³／空気ｍ³）となる。
【００７６】
＜第２の実施の形態（排出源が複数の場合）＞
排出源が複数（ｊ個）存在し、それぞれの排出源から、時間変化する異なる排出量（ｑｊ（ｔ））で物質が排出されている場合、第１の実施の形態で使用したLagrangian粒子拡散モデルの粒子情報（位置，放出後経過時間）に加えて、各粒子に排出源識別情報（ｓｉ）を持たせることにより、第１の実施の形態と同じ機能を発揮することが可能となる。
【００７７】
例えば図１１に示すように、２つの排出源Ｓ１，Ｓ２がある場合には、排出源Ｓ１から排出される粒子には排出源識別情報ｓ１を持たせ、排出源Ｓ２から排出される粒子には排出源識別情報ｓ２を持たせる。そして、Lagrangian粒子拡散モデルを使用して、物質を粒子に置換し、両排出源Ｓ１，Ｓ２からそれぞれ毎秒Ｎ個の粒子を発生させ、各粒子の挙動を数値計算して粒子の位置を示す空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））を求める。更に、各粒子が有する情報である空間座標（ｘｉ（ｔ），ｙｉ（ｔ），ｚｉ（ｔ））の他に、粒子を最初に発生させた時点からの経過時間である排出後経過時間Ｔｉ（ｔ）と、排出源識別情報ｓ１またはｓ２を持たせる。
【００７８】
図１２（ａ）に示すように、排出源Ｓ１，Ｓ２から排出される排出量Ｑがそれぞれ一定（＝１）として、従来のLagrangian粒子拡散モデルを利用して、図１２（ｂ）に示すような格子領域Ｉ，Ｊ，Ｋにおける物質の濃度時間変化を求めて記録しておく。
【００７９】
次に、排出源Ｓ１から排出される物質の排出量が、図１３（ａ）に点線で示すように、時間変化するｑ１（ｔ）であるとすると、この排出源Ｓ１から放出された粒子については、第１の実施の形態と同じ方法で、粒子の排出源強度をｑ１（ｔ）とする。そして、各粒子毎に、その粒子が発生した時点の排出源強度ｑ１（ｔ）を参照して、排出源強度を設定する。この結果、求める格子領域Ｉ，Ｊ，Ｋにおいて、排出源識別情報ｓ１を有する粒子の排出源強度を積算することにより、排出源Ｓ１から排出された物質の濃度（求める格子領域Ｉ，Ｊ，Ｋにおける物質濃度）の時間変化を求めることができる。
【００８０】
同様に、排出源Ｓ２から排出される物質の排出量がｑ２（ｔ）であるとすると、この排出源Ｓ２から放出された粒子については、第１の実施の形態と同じ方法で、粒子の排出源強度をｑ２（ｔ）とする。そして、各粒子毎に、その粒子が発生した時点の排出源強度ｑ２（ｔ）を参照して、排出源強度を設定する。この結果、求める格子領域Ｉ，Ｊ，Ｋにおいて、排出源識別情報ｓ２を有する粒子の排出源強度を積算することにより、排出源Ｓ２から排出された物質の濃度（求める格子領域Ｉ，Ｊ，Ｋにおける物質濃度）の時間変化を求めることができる。
【００８１】
このようにして、排出源Ｓ１から排出された物質の濃度の時間変化と、排出源Ｓ２から排出された物質の濃度の時間変化とを加算することにより、求める格子領域Ｉ，Ｊ，Ｋにおける物質の濃度を求めることができる。
【００８２】
＜第３の実施の形態＞
第３の実施の形態では、例えば、ガス漏洩事故が発生した後、ガスの排出量の実測結果に基づいて、短時間で現状および将来の濃度分布を計算する方法である。
【００８３】
環境アセスメントなどの拡散計算は、緊急に計算結果が必要ではないので、時間の経過に沿い変化するガスの排出量ｑ（ｔ））を決定してから、数日から数ケ月の時間をかけてガスの拡散状況の予測演算を実施する。しかし、ガス漏洩事故などのような緊急時には、周辺住民の非難対策を緊急にとる必要があるので、できるだけ短時間に拡散計算結果を出力する必要がある。
【００８４】
このような場合、各時刻の３次元空間風速分布に基づいて、図１４（ａ）（ｂ）（ｃ）に示すように、第１の実施の形態と同じLagrangian粒子拡散モデルで、一定放出率（Ｑ＝１）の拡散計算を２４時間連続して実施しておく。
【００８５】
ガス漏洩事故発生後、排出源である煙突の出口に設置されたガス濃度測定装置などから、実測ガス排出量ｑ（ｔ）を測定する。この実測ガス排出量ｑ（ｔ）に応じて、各粒子の排出源強度ｑ（ｔ）を設定することにより、第１の実施の形態と同じ方法で、このガス排出量に対応した濃度時間変化を補正計算することができる（図１４（ｄ）（ｅ）参照）。つまり、各粒子毎に、その粒子が発生した時点の排出源強度ｑ（ｔ）を参照して、排出源強度を設定する（図１４（ｄ））。そして、求める格子領域Ｉ，Ｊ，Ｋにおいて、その領域に存在する粒子の排出源強度を積算することにより、排出源Ｓから排出された物質の濃度（求める格子領域Ｉ，Ｊ，Ｋにおける物質濃度）の時間変化を求めることができる。
【００８６】
なお、将来のガス排出量は、現状の実測ガス排出量ｑ（ｔ）と等価とするか、または、別途、設定される予測式に基づいて計算する。
【００８７】
現在及び将来の３次元空間風速分布を計算する方法としては、気象庁などから定期的に配信される広域格子（２０ｋｍ）の気象データ（ＧＰＶ：Grid Point Value）を利用して、広域気象予測モデル（ＲＡＭＳ，ＭＭ５など）で詳細格子（数ｋｍ〜数１０ｍ）の気象データ（風向，風速，気温）の時間変化を計算する方法がある。
【００８８】
＜第４の実施の形態＞
第４の実施の形態は、例えば、ガス漏洩事故が発生した後、漏洩箇所周辺の濃度実測結果に基づいて、短時間で現状のガスの排出量を計算する方法である。
【００８９】
ガス漏洩事故が発生した場合、図１５（ａ）に示す排出源であるガス漏洩源Ｓからの排出量ｑ（ｔ）が測定困難な場合が多い。このような場合、漏洩箇所周辺観測点（ｘｋ）で実測された濃度（ｃｋ（ｔ）観測）の時間変化から、第１の実施の形態のLagrangian粒子拡散モデルを用いて、短時間で排出量ｑ（ｔ）を推定することができる。
【００９０】
漏洩事故発生前に、第１の実施の形態のLagrangian粒子拡散モデルを用いて、一定排出量（Ｑ＝１）で（図１５（ｂ）参照）、ある周辺観測点（ｘｋ）の濃度（Ｃｋ（ｔ）計算）を計算しておく（図１５（ｃ）参照）。
【００９１】
この濃度（Ｃｋ（ｔ）計算）は、観測点（ｘｋ）を中心とする３次元空間体積（Δｘ×Δｙ×Δｚ）中に存在する粒子（ｎ個）から、次の式（１０）で計算されたものである。
【００９２】
Ｃｋ（ｔ）計算＝ｎ×ΣＱ（ｔ−Ｔｉ）／Ｎ／（Δｘ×Δｙ×Δｚ）・・(10)
この場合、排出量Ｑ（ｔ−Ｔｉ）は一定値（＝１）である。
【００９３】
次に、観測点（ｘｋ）の濃度時間変化（Ｃｋ（ｔ）観測）が測定されると、（１０）式を利用して（１１）式が求まる。
【００９４】

【００９５】
ここで、ｎ１，ｎ２，ｎ３，及びｎｌは、放出開始後、０秒前、ΔＴ秒前、（ΔＴ×２）秒前、（ΔＴ×３）秒前及びｔ秒前に放出された粒子で、観測点（ｘｋ）を中心とする３次元空間体積（Δｘ×Δｙ×Δｚ）中に存在する粒子の総数である。
また、Ｃ０（ｔ）は、バックグラウンド濃度と呼ばれ、計算対象とする排出源以外から放出され、観測地点に関係なく存在する濃度であり、一定値または時間ｔによって変化する。
ｑ０，ｑ１，ｑ２，ｑ３およびｑｌは、ｔ秒時から０秒前、ΔＴ秒前、（ΔＴ×２）秒前、（ΔＴ×３）秒前およびｔ（＝ΔＴ×１）秒前の放出率である。
【００９６】
Ｃｋ（ｔ）観測は、２４時間連続して測定されているので、漏洩開始時刻から０秒後、ΔＴ秒後、（ΔＴ×２）秒後、（ΔＴ×３）秒後およびｔ（＝ΔＴ×１）秒後まで（１＋ｌ）個以上、測定することが可能である。
また、観測点ｘｋが、ｋ個ある場合は、ｋ×（１＋ｌ）個の観測データを求めることができる。
【００９７】
（１１）式の未知数はｑ０、ｑ１、ｑ２、ｑ３及びｑｌの（１＋ｌ）個、既知数Ｃｋ（ｔ）観測は、ｋ×（１＋ｌ）個以上あるので、未知数が既知数より少ない。この場合、最小自乗法を用いて、最も、観測濃度Ｃｋ（ｔ）観測に対する最小自乗誤差を小さくするように、未知数のｑ０、ｑ１、ｑ２、ｑ３及びｑｌを決定することができる。なお、図１６は第４の実施の形態の計算状態を示すフロー図である。
【００９８】
＜第５の実施の形態＞
第５の実施の形態は、例えば、ガス漏洩事故が発生した後、漏洩箇所周辺の濃度実測結果に基づいて、短時間で現状のガスの排出量を推定し、濃度分布を計算する方法である。
【００９９】
第４の実施の形態により、漏洩箇所周辺観測点（ｘｋ）で実測された濃度（Ｃｋ（ｔ）観測）の時間変化から、第１の実施の形態のLagrangian粒子拡散モデルを用いて、短時間で排出量ｑ（ｔ）を推定することができる。
【０１００】
図１７に示すように、漏洩事故発生前に、第１の実施の形態のLagrangian粒子拡散モデルを用いて、一定放出率（Ｑ＝１）で、ある周辺観測点（ｘｋ）の濃度（Ｃｋ（ｔ）計算）を計算しておく。
【０１０１】
この濃度（Ｃｋ（ｔ）計算）は、観測点（ｘｋ）を中心とする３次元空間体積（Δｘ×Δｙ×Δｚ）中に存在する粒子（ｎ個）から、次の式（１２）で計算されたものである。
Ｃｋ（ｔ）計算＝ｎ×ΣＱ（ｔ−Ｔｉ）／Ｎ／（Δｘ×Δｙ×Δｚ）・・(12)
この場合、放出率Ｑ（ｔ−Ｔｉ）は一定値（＝１）である。
【０１０２】
次に、観測点（ｘｋ）の濃度時間変化（（Ｃｋ（ｔ）観測）が測定されると、（１２）式を利用して（１３）式が求まる。
【０１０３】

【０１０４】
ここで、ｎ１，ｎ２，ｎ３，及びｎｌは、放出開始後、０秒前、ΔＴ秒前、（ΔＴ×２）秒前、（ΔＴ×３）秒前及びｔ秒前に放出された粒子で、観測点（ｘｋ）を中心とする３次元空間体積（Δｘ×Δｙ×Δｚ）中に存在する粒子の総数である。
また、Ｃ０（ｔ）は、バックグラウンド濃度と呼ばれ、計算対象とする排出源以外空放出され、観測地点に関係なく存在する濃度であり、一定値または時間ｔによって変化する。
ｑ０，ｑ１，ｑ２，ｑ３およびｑｌは、ｔ秒時から０秒前、ΔＴ秒前、（ΔＴ×２）秒前、（ΔＴ×３）秒前およびｔ（＝ΔＴ×１）秒前の放出率である。
【０１０５】
Ｃｋ（ｔ）観測は、２４時間連続して測定されているので、漏洩開始時刻から０秒後、ΔＴ秒後、（ΔＴ×２）秒後、（ΔＴ×３）秒後およびｔ（＝ΔＴ×１）秒後まで（１＋ｌ）個以上、測定することが可能である。
また、観測点ｘｋが、ｋ個ある場合は、ｋ×（１＋ｌ）個の観測データを求めることができる。
【０１０６】
（１３）式の未知数はｑ０、ｑ１、ｑ２、ｑ３及びｑｌの（１＋ｌ）個、既知数Ｃｋ（ｔ）観測は、ｋ×（１＋ｌ）個以上あるので、未知数が既知数より少ない。この場合、最小自乗法を用いて、最も、観測濃度Ｃｋ（ｔ）観測に対する最小自乗誤差を小さくするように、未知数のｑ０、ｑ１、ｑ２、ｑ３及びｑｌを決定することができる。
【０１０７】
この推定値ｑ０、ｑ１、ｑ２、ｑ３及びｑｌを用いて、第１の実施の形態の方法で、濃度時間変化の補正計算を行えば、各経過時間毎の濃度分布を計算することができる。
【０１０８】
＜第６の実施の形態＞
第６の実施の形態は、ガス漏洩事故が発生した後、放出量実測結果に基づいて、短時間で現状および将来の濃度分布計算結果をインターネットで提供するシステムである。
【０１０９】
ガス漏洩事故などのような緊急時には、周辺住民の避難対策を緊急にとる必要があるので、できるだけ、短時間に拡散計算結果を出力する必要がある。しかし、このような漏洩事故はいつ発生するか分からないため、監督官庁である消防、警察及び自治体と事業者である各工場では、このための計算機システムを２４時間体制で管理運用する必要がある。
【０１１０】
この管理運用には、多大の労力と費用および高度な計算機運用技術を必要とするので、大規模な危機管理システムを常設する組織以外では、管理運用が困難である。
【０１１１】
この問題点を補う方法として、最近のインターネットを利用した情報提供システムを案出した。
【０１１２】
この第６の実施の形態では、図１８に示すように、監督官庁である消防、警察及び自治体などの監督官庁１０や、事業者である各工場１１とは別の場所に、安全解析センター１２を設置する。安全解析センター１２では、インターネット等の情報伝達手段を介して、気象庁１３などの気象データ配信施設から配信される気象観測データを受信し、平常時は大型計算機を用いて計算をすることにより、各時刻の３次元空間風速分布に基づいて、第１の実施の形態と同じLagrangian粒子拡散モデルで、一定放出率（Ｑ＝１）の拡散計算を２４時間連続して実施しておく。
【０１１３】
ガス漏洩事故発生後、事業者１１の煙突出口に設置されたガス濃度測定装置などから、実測したガス放出量ｑ（ｔ）がわかれば、このガス放出量ｑ（ｔ）をインターネット等の情報伝達手段を介して、安全解析センター１２に伝送する。そうすると、安全解析センター１２では、第１の実施の形態と同じ方法で、このガス放出率に対応した濃度時間変化を補正計算することができる。
【０１１４】
安全解析センター１２は、計算した濃度計算結果を消防、警察及び自治体などの監督官庁１０に送信する。監督官庁１０は、濃度に応じて、事業者１１や工場周辺住民１４に避難勧告をする。
【０１１５】
なお、将来のガス放出率は、現状の実測したガス放出率ｑ（ｔ）と等価とするか、または、別途、設定される予測式に基づいて計算する。なお、複数の事業者１１が存在する場合には、安全解析センター１２では、第２の実施の形態に示す演算をして濃度時間変化を予測する。
【０１１６】
現在及び将来の３次元空間風速分布を計算する方法としては、気象庁１３などから定期的に配信される広域格子（２０ｋｍ）の気象観測データ（ＧＰＶ：Grid Point Value ）を利用して、広域気象予測モデル（ＲＡＭＳ，ＭＭ５など）で詳細格子（数ｋｍ〜数１０ｍ）の気象データ（風向，風速，気温）の時間変化を計算する方法がある。
【０１１７】
【発明の効果】
以上実施の形態と共に具体的に説明したように、本発明にかかる拡散物質の拡散状況予測方法では、拡散物質の拡散状況をコンピュータによる演算処理により求める拡散物質の拡散状況予測方法であって、
前記コンピュータは、
排出源から大気中に排出された物質が大気中を拡散していく状況を予測するため、前記物質を多数の粒子に置換して、排出源の位置から演算周期毎に予め設定した個数の粒子が発生すると設定すると共に、
排出源の位置を含む領域内の多数の地点における、時間の経過に沿って変化する風向・風速を示す風速場データを、粒子の拡散状態を演算する拡散方程式に代入することにより、各粒子の拡散速度を求め、この拡散速度から各粒子が存在する空間位置を示す空間座標を各演算周期毎に求めると共に、前記粒子を最初に発生させた時点からの経過時間である排出後経過時間を計測し、各演算周期における各粒子の空間座標と各粒子の排出後経過時間を対応させてデータ記録装置に記録しておき、
排出された物質の排出後経過時間の時間経過に伴う排出量の変化に比例させて、排出後経過時間の時間経過に沿う粒子に対する排出源強度データを設定しておき、
前記データ記録装置に記録していた、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間を読み出すと共に、読み出した排出後経過時間を参照して各粒子が発生した時点を求め、この時点における各粒子の排出源強度を前記排出源強度データから求め、前記データ記録装置に、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と排出源強度を対応させて再記録させ、
所定の演算周期における所定の領域の前記物質の濃度は、当該所定の演算周期における当該所定の領域に存在する全ての粒子の排出源強度を積算することにより求めるようした。
このため、排出源から排出される量が異なっていても、特定領域における物質の濃度を短時間で演算することができ。
【０１１８】
また本発明の拡散物質の拡散状況予測方法は、拡散物質の拡散状況をコンピュータによる演算処理により求める拡散物質の拡散状況予測方法であって、
前記コンピュータは、
複数の排出源から大気中に排出された物質が大気中を拡散していく状況を予測するため、前記物質を多数の粒子に置換して、各排出源の位置から演算周期毎に予め設定した個数の粒子がそれぞれ発生すると設定すると共に、
排出源の位置を含む領域内の多数の地点における、時間の経過に沿って変化する風向・風速を示す風速場データを、粒子の拡散状態を演算する拡散方程式に代入することにより、各粒子の拡散速度を求め、この拡散速度から各粒子が存在する空間位置を示す空間座標を各演算周期毎に求めると共に、前記粒子を最初に発生させた時点からの経過時間である排出後経過時間を計測し、各演算周期における各粒子の空間座標と各粒子の排出後経過時間と排出源を識別する排出源識別情報とを対応させてデータ記録装置に記録しておき、
各排出源から排出された物質の排出後経過時間の時間経過に伴う排出量の変化に比例させて、排出後経過時間の時間経過に沿う粒子に対する排出源強度データを各排出源毎にそれぞれ設定しておき、
前記データ記録装置に記録していた、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と各粒子の排出源識別情報を読み出すと共に、読み出した排出後経過時間を参照して各粒子が発生した時点を求め、読み出した排出源識別情報を参照して粒子が発生した時点における各粒子の排出源強度をその粒子が発生した排出源に対応した前記排出源強度データから求め、前記データ記録装置に、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と排出源強度を対応させて再記録させ、
所定の演算周期における所定の領域の前記物質の濃度は、当該所定の演算周期における当該所定の領域に存在する全ての粒子の排出源強度を積算することにより求めるようにした。
このため、排出源から排出される量が異なっていても、特定領域における物質の濃度を短時間で演算することができる。また、複数の排出源から物質が排出されていても、正確に物質の濃度計算をすることができる。
【０１１９】
また本発明の拡散物質の拡散状況予測方法では、
前記排出源強度データは、前記排出源から実際に排出される物質の濃度を実測することにより求めて設定されたり、
前記排出源強度データは、前記排出源の周囲の観測点で実測した物質の濃度の時間変化を基に設定される。
このため、排出源から排出される物質の濃度データが予め求まっていなくても、実測データを用いて濃度計算をすることができる。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態における計算フローを示すフローチャートである。
【図２】本発明の第１の実施の形態における粒子の拡散状態を示す説明図である。
【図３】本発明の第１の実施の形態における粒子の拡散状態を示す説明図である。
【図４】本発明の第１の実施の形態における粒子の拡散状態を示す説明図である。
【図５】物質の排出量の時間変化を一例を示す特性図である。
【図６】物質の排出量の時間変化に対応する排出源強度の一例を示す特性図である。
【図７】所定の格子領域における粒子分布を示す説明図である。
【図８】排出源と格子領域を示す説明図である。
【図９】排出量が一定の場合の排出量と濃度との関係を示す特性図である。
【図１０】排出量が時間変化する場合の排出量と濃度との関係を示す特性図である。
【図１１】２つの排出源と格子領域を示す説明図である。
【図１２】排出量が一定の場合の排出量と濃度との関係を示す特性図である。
【図１３】排出量が時間変化する場合の排出量と濃度との関係を示す特性図である。
【図１４】第３の実施の形態を示す説明図である。
【図１５】第４の実施の形態を示す説明図である。
【図１６】第４の実施の形態における計算フローを示すフローチャートである。
【図１７】第５の実施の形態における計算フローを示すフローチャートである。
【図１８】第６の実施の形態にかかるシステムを示すシステム構成図である。
【図１９】従来技術における粒子の拡散状態を示す説明図である。
【図２０】従来技術における粒子の拡散状態を示す説明図である。
【図２１】従来技術における粒子の拡散状態を示す説明図である。
【図２２】所定の格子領域における粒子分布を示す説明図である。
【図２３】粒子拡散モデルの機能を示す説明図である。
【図２４】排出源と格子領域を示す説明図である。
【図２５】排出量が時間変化する場合の排出量と濃度との関係を示す特性図である。
【図２６】排出量が一定の場合の排出量と濃度との関係を示す特性図である。
【図２７】排出量が瞬間的な場合の排出量と濃度との関係を示す特性図である。
【符号の説明】
１データ記録装置
１０監督官庁
１１事業者
１２安全解析センター
１３気象庁
１４工場周辺住民
Ｓ，Ｓ１，Ｓ２排出源

Claims

拡散物質の拡散状況をコンピュータによる演算処理により求める拡散物質の拡散状況予測方法であって、
前記コンピュータは、
排出源から大気中に排出された物質が大気中を拡散していく状況を予測するため、前記物質を多数の粒子に置換して、排出源の位置から演算周期毎に予め設定した個数の粒子が発生すると設定すると共に、
排出源の位置を含む領域内の多数の地点における、時間の経過に沿って変化する風向・風速を示す風速場データを、粒子の拡散状態を演算する拡散方程式に代入することにより、各粒子の拡散速度を求め、この拡散速度から各粒子が存在する空間位置を示す空間座標を各演算周期毎に求めると共に、前記粒子を最初に発生させた時点からの経過時間である排出後経過時間を計測し、各演算周期における各粒子の空間座標と各粒子の排出後経過時間を対応させてデータ記録装置に記録しておき、
排出された物質の排出後経過時間の時間経過に伴う排出量の変化に比例させて、排出後経過時間の時間経過に沿う粒子に対する排出源強度データを設定しておき、
前記データ記録装置に記録していた、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間を読み出すと共に、読み出した排出後経過時間を参照して各粒子が発生した時点を求め、この時点における各粒子の排出源強度を前記排出源強度データから求め、前記データ記録装置に、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と排出源強度を対応させて再記録させ、
所定の演算周期における所定の領域の前記物質の濃度は、当該所定の演算周期における当該所定の領域に存在する全ての粒子の排出源強度を積算することにより求めることを特徴とする拡散物質の拡散状況予測方法。
拡散物質の拡散状況をコンピュータによる演算処理により求める拡散物質の拡散状況予測方法であって、
前記コンピュータは、
複数の排出源から大気中に排出された物質が大気中を拡散していく状況を予測するため、前記物質を多数の粒子に置換して、各排出源の位置から演算周期毎に予め設定した個数の粒子がそれぞれ発生すると設定すると共に、
排出源の位置を含む領域内の多数の地点における、時間の経過に沿って変化する風向・風速を示す風速場データを、粒子の拡散状態を演算する拡散方程式に代入することにより、各粒子の拡散速度を求め、この拡散速度から各粒子が存在する空間位置を示す空間座標を各演算周期毎に求めると共に、前記粒子を最初に発生させた時点からの経過時間である排出後経過時間を計測し、各演算周期における各粒子の空間座標と各粒子の排出後経過時間と排出源を識別する排出源識別情報とを対応させてデータ記録装置に記録しておき、
各排出源から排出された物質の排出後経過時間の時間経過に伴う排出量の変化に比例させて、排出後経過時間の時間経過に沿う粒子に対する排出源強度データを各排出源毎にそれぞれ設定しておき、
前記データ記録装置に記録していた、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と各粒子の排出源識別情報を読み出すと共に、読み出した排出後経過時間を参照して各粒子が発生した時点を求め、読み出した排出源識別情報を参照して粒子が発生した時点における各粒子の排出源強度をその粒子が発生した排出源に対応した前記排出源強度データから求め、前記データ記録装置に、各演算周期毎の各粒子の空間座標と各粒子の排出後経過時間と排出源強度を対応させて再記録させ、
所定の演算周期における所定の領域の前記物質の濃度は、当該所定の演算周期における当該所定の領域に存在する全ての粒子の排出源強度を積算することにより求めることを特徴とする拡散物質の拡散状況予測方法。
請求項１または請求項２において、前記排出源強度データは、前記排出源から実際に排出される物質の濃度を実測することにより求めて設定されることを特徴とする拡散物質の拡散状況予測方法。
請求項１または請求項２において、前記排出源強度データは、前記排出源の周囲の観測点で実測した物質の濃度の時間変化を基に設定されることを特徴とする拡散物質の拡散状況予測方法。