JP2000122521A

JP2000122521A - 拡散現象シミュレータ、拡散現象解析システム、拡散現象のシミュレーション方法および記憶媒体

Info

Publication number: JP2000122521A
Application number: JP29163798A
Authority: JP
Inventors: Makoto Morishita; 信森下
Original assignee: Individual
Current assignee: Individual
Priority date: 1998-10-14
Filing date: 1998-10-14
Publication date: 2000-04-28

Abstract

(57)【要約】【課題】複雑な処理や複雑な方程式の解法を用いるこ
となく、短時間で精度良く拡散現象を解析するシミュレ
ーションシステムを提供する。【解決手段】拡散シミュレータは、入力されたデータ
などを下に、解析対象となる空間を複数に分割してセル
を生成し、各セルに対応して、当該セルに対応する領域
に含まれる粒子の量を示す第１の状態量を定義するセル
状態量設定部３６、セル管理部３２と、単位時間ごと
の、セルの周囲に位置する局所近傍を構成する他のセル
への第１の状態量の移動を定義するシミュレーション条
件設定／記憶部３８と、定義された他のセルへの状態量
の移動に基づき、セルごとに、定義された第１の状態量
を変更することにより、粒子の移動を実現するシミュレ
ーション実行部４０とを備えている。

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、分子や微粒子が熱
エネルギーにより移動する現象である拡散をシミュレー
トする拡散解析シミュレータに関し、特に、海、河川、
湖沼などの水圏に物質が投入されたときの拡散状況をシ
ミュレートするのに好適な拡散解析シミュレーション手
法に関する。

【０００２】

【従来の技術】従来、粒子などの拡散をシミュレートお
よび解析するために数多くの手法が提案されている。こ
れらは、基本的には、拡散方程式の近似解を求めること
により、シミュレーションを実現している。たとえば、
特開平６−２８３４５８号公報には、半導体の製造プロ
セスにおける不純物の拡散をシミュレートする技術が開
示されている。また、特開平９−７９８９９号公報に
は、熱および流体の運動方程式に基づき所定の位置にお
ける流速および圧力を求めるシミュレーション技術が開
示されている。

【０００３】また、拡散現象の一例として、海上でのタ
ンカーの座礁、老朽化による船舶の破損などが原因とな
って流出した油や、河川や湖沼に廃棄された汚染物質の
拡散が知られている。これらの現象において、その拡散
状況、拡散過程を短時間で精度良く予測し、その対策を
立案することは環境保全の観点から重要となっている。
たとえば、流出した油の拡散過程に関する従来の研究と
して、海面拡散モデルに関して、ある時点で放出された
油が同心円状に広がることを前提に次元解析により導出
されたフェイ(Fay)の式、初期の流出油の形状がある半
径を有する円筒と仮定して、エネルギー保存則から拡散
半径をもとめた元良の式、ナビエ・ストークス方程式お
よび拡散方程式を連立させて水中の質量拡散を考慮する
手法、油を平板上に広がる流体とみなして浅水方程式を
用いる手法、ランダムウォークモデルを用いた手法など
が提案されている。

【０００４】近年のコンピュータシミュレーションにお
いては、ランダムウォークモデルを用いて油の拡散過程
を解析する場合が多いが、特に、解析平面上に汚染物質
の量に相当する粒子を配置し、その粒子をすべての方向
に満遍なく、かつ、統計的分布にしたがって連続的な距
離を飛ばす方法が採用されている。

【０００５】

【発明が解決しようとする課題】しかしながら、上記ラ
ンダムウォークモデルを用いた場合には、解析領域の大
きさ或いはその解析精度によっては、計算量が著しく膨
大なものとなり、多大の計算時間が必要となる。したが
って、現実に汚染が生じている場合に、汚染現場（たと
えば、海上で汚染を監視する船舶）と、陸上の施設（た
とえば、研究施設）に設置された大型コンピュータシス
テムとの通信を確保し、情報の授受を行うことが不可避
となる。その一方、ナビエ・ストークス方程式および拡
散方程式を連立させる手法も、巨大なメモリ空間を有す
る大型スーパーコンピュータを用いる必要があり、か
つ、このようなコンピュータを用いても著しく多大な計
算時間を要する。

【０００６】本発明は、複雑な処理や複雑な方程式の解
法を用いることなく、短時間で精度良く拡散現象を解析
するシミュレーションシステムを提供することを目的と
する。特に、小型コンピュータを用いて、短時間でかつ
高精度に油、汚染物質などの拡散シミュレーションを実
行するシステムを提供することを目的とする。

【０００７】

【課題を解決するための手段】本発明の目的は、解析対
象となる空間を複数に分割してセルを生成するセル生成
手段と、各セルに対応して、当該セルに対応する領域に
含まれる粒子の量を示す第１の状態量を定義する状態量
定義手段と、単位時間ごとの、セルの周囲に位置する局
所近傍を構成する他のセルへの第１の状態量の移動を定
義する移動量定義手段と、前記定義された他のセルへの
状態量の移動に基づき、前記セルごとに、前記定義され
た第１の状態量を変更することにより、粒子の移動を実
現するシミュレーション手段とを備えたことを特徴とす
る拡散現象シミュレータにより達成される。本発明によ
れば、局所近傍系を形成するセルに粒子を移動させるこ
とによりシミュレーションが実現できる。したがって、
方程式の解析解を求めるための複雑な処理を行うことな
く、短時間にて解析現象をシミュレートすることが可能
となる。

【０００８】本発明の好ましい実施態様においては、前
記第１の状態量が、対応するセルの領域に含まれる粒子
の数にて表わされる。すなわち、所定の領域に含まれる
粒子の数により、粒子の濃度を表わしている。

【０００９】本発明のさらに好ましい実施態様において
は、前記シミュレーション手段が、単位時間ごとに、セ
ル中の粒子を、当該セルに隣接する局所近傍を構成する
セルに移動させるように構成され、前記セルの大きさお
よび前記単位時間により、拡散係数が決定される。すな
わち、セルの大きさおよび単位時間を変更することによ
り、所望の等価な拡散係数を設定することが可能とな
り、より精度よく拡散現象をシミュレートすることが可
能となる。

【００１０】また、本発明のさらに好ましい実施態様に
おいては、前記状態量定義手段が、前記セルに対応し
て、粒子が通過し得るか否かを示す第２の状態量を定義
している。これにより、たとえば、海面や空間と、陸地
や障害物などを定義することが可能となる。前記移動量
定義手段は、前記第２の状態量に基づき、局所近傍を構
成するセルのうち、粒子を移動させることが可能なセル
を選択するように構成されていても良い。また、前記シ
ミュレーション手段が、前記第２の状態量が障害物であ
ることを示すセルにおいて、粒子を蓄積させるように作
動しても良い。

【００１１】さらに、本発明の好ましい実施態様におい
ては、前記移動量定義手段が、前記局所近傍を構成する
セルの各々に粒子を移動させる移動確率を決定する。ま
た、移動確率は、移流速度を考慮して決定されるのが更
に好ましい。たとえば、移動確率は、移流方向に向かう
粒子の移動確率と、その反対方向に向かう粒子の移動確
率との差が、前記移流速度と比例するように決定するこ
とができる。

【００１２】また、本発明の目的は、上記構成の拡散現
象シミュレータと、風力および風向を示すデータに基づ
き、移流速度を算出する速度算出手段とを備え、前記速
度算出手段にて算出された移流速度に基づき、移動確率
が決定されるように構成された拡散現象解析システムに
よっても達成される。この発明によれば、風力および風
向の観測値に基づき移流速度が算出されるため、現場な
どにて、或いは、現場からの情報に基づき、極めて精度
良く、拡散状況を予測ないし推定することが可能とな
る。

【００１３】上記拡散現象解析システムは、さらに、前
記速度算出手段が、前記風力および風向を示すデータ
と、潮流の大きさおよび向きを示す潮流データとに基づ
き、海上に浮遊する粒子の移流速度を算出するように構
成されていても良い。これは、特に、海上を浮遊する粒
子の拡散の状況を予測ないし推定するのに好ましい。

【００１４】なお、上記セル生成手段などは、図１に示
すＣＰＵにより実現される。より詳細には、セル生成手
段および状態量定義手段は、主として、図２に示すセル
状態量設定部およびセル管理部により、移動量定義手段
は、主として、シミュレーション条件設定／記憶部によ
り、シミュレーション手段は、シミュレーション実行部
により実現することが可能である。

【００１５】また、本発明の目的は、上記シミュレータ
にて実現する処理ステップを備えた方法、および、上記
処理ステップからなるプログラムを格納した記憶媒体に
よっても達成される。

【００１６】

【発明の実施の形態】以下、添付図面を参照して、本発
明の実施の形態につき説明を加える。図１は、本発明の
第１の実施の形態にかかる拡散現象シミュレータのハー
ドウェア構成を示すブロックダイヤグラム、図２は、こ
の拡散現象シミュレータの機能を示すブロックダイヤグ
ラムである。

【００１７】図１に示すように、拡散現象シミュレー
タ１０は、ＣＰＵ１２と、メモリ１４と、外部記憶装置
１６と、インタフェース（Ｉ／Ｆ）１８と、マウス、キ
ーボードなどの入力装置２０と、ＣＲＴディスプレイや
ＬＣＤからなる表示装置２２と、プリンタ２４とから構
成されている。ＣＰＵ１２、メモリ１４およびＩ／Ｆ１
８は、データバス２６を介して接続されている。また、
外部記憶装置１６、入力装置２０、表示装置２２および
入力装置２４は、Ｉ／Ｆ１８を介してデータバス２４と
接続でき、これにより、ＣＰＵ１２との間で、データの
授受が可能となっている。

【００１８】メモリ１４は、シミュレーション中のデー
タや、作動中のシミュレーションプログラムを一時的に
記憶する主記憶を有している。また、外部記憶装置１６
には、本実施の形態にかかる拡散現象シミュレータ１０
を作動させるのに必要なプログラムが記憶されている。
上記プログラムは、ＣＤ−ＲＯＭやフロッピーディスク
に格納され、ＣＤ−ＲＯＭ駆動装置やフロッピーディス
ク駆動装置（図示せず）により読み出され、メモリ１４
の主記憶に記憶され、これにより、必要な処理が実行さ
れるように構成しても良い。

【００１９】この拡散現象シミュレータは、図２に示す
ような機能を有すると考えることができる。すなわち、
拡散現象シミュレータ１０は、後述するセルを管理する
セル管理部３２と、セルごとの状態量を記憶するセル状
態量記憶部３４と、セルに対応する状態量を設定する状
態量設定部３６と、シミュレーション条件を設定すると
ともに、この条件を記憶するシミュレーション条件設定
／記憶部３８と、所定のシミュレーション条件に基づい
てシミュレーションを実行するシミュレーション実行部
４０と、当該シミュレーション結果を記憶するシミュレ
ーション結果記憶部４２とを備えている。

【００２０】図２の構成部分のうち、セル管理部３２、
セル状態量設定部３６、シミュレーション実行部４０の
機能は、主として、ＣＰＵ１２により実現され、シミュ
レーション条件設定／記憶部３８の機能は、ＣＰＵ１２
およびメモリ１４により実現される。さらに、セル状態
量記憶部３４およびシミュレーション結果記憶部４２の
機能は、メモリ１４および外部記憶装置１６により実現
される。

【００２１】以下に、本実施の形態にかかる拡散現象シ
ミュレータの動作原理を説明する。本実施の形態におい
ては、セルラオートマトン（Cellular Automata; Ｃ
Ａ）を用いてモデル化を行っている。すなわち、解析領
域をセルと称する区分領域に分割し、各セル上の内部状
態を表わす離散的状態量を定義し、近傍のセル同志の相
互作用を定める局所近傍則を用いて状態量を離散時間の
経過に沿って遷移させている。

【００２２】まず、一次元の拡散現象のシミュレーショ
ンにつき説明する。一般に、一次元の移流拡散現象を表
わす微分方程式は、一次元座標をｘ、時間をｔとすれ
ば、以下の（１）式にて与えられることが知られてい
る。

【００２３】

【数１】ここで、ｕは濃度、Ｄは拡散係数、ｖは移流速度を表わ
す。

【００２４】このような微分方程式で表わされる一次元
の移流拡散現象を、本実施の形態においては、以下のよ
うな二次元ＣＡを用いて解析している。まず、図３に示
すようにｙ方向に単位幅を有する解析平面を考え、この
解析平面を所定の数のセル（たとえば、２０×２０００
個のセル）に分割する。各座標（ｘ，ｙ）に対応して、
図４に示すように、壁面／水面を示す状態量（以下、便
宜上「第２の状態量」と称する）（符号４０１参照）お
よび粒子数を示す第１の状態量（符号４０２参照）（以
下、便宜上「第１の状態量」と称する。が与えられる。
図３の例では、ｙ＝０および２０の位置に対応する第２
の状態量は壁面を示し、その一方、他の位置に対応する
第２の状態量は水面を示す。また、あるセルに対応する
第２の状態量が水面を示す場合には、このセルに含まれ
る粒子数（第１の状態量）が割り当てられる。この実施
の形態では、粒子数がセルの領域に存在する粒子の濃度
と考えた。図３の例において、各セルに割り当てられた
第２の状態量から、このモデルは、単位幅を有する容器
と考えることができる。すなわち、第２の状態量が壁面
を示すセルが、溶媒（粒子）を入れた容器の内壁と考え
ることができる。

【００２５】さて、図３に示すモデルにおいて、ｘ＝０
の位置では、常に最大濃度（＝１）の溶媒（粒子）が供
給され、ｘ＝２０００の位置では、常に溶媒の濃度が０
となると考える。すなわち、溶媒の濃度ｕは、以下の
（２）式にてｕ（０，ｔ）＝１、ｕ（２０００，ｔ）＝０・・・・（２）本実施の形態では、（２）式に対応するように、平面の
左端（ｘ＝０の位置）においては、各セルに粒子が常に
最大値（１０００）だけ存在し、その一方、平面の右端
（ｘ＝２０００の位置）では、各セルに粒子が存在しな
いように、セルの第２の状態量を設定している。

【００２６】次に、局所近傍則につき説明を加える。こ
の実施の形態においては、一つのセル上の粒子を、単位
時間経過後に他のセルに振り分けるように定義した。こ
こで、図５に示すように、あるセルに定義された１個の
粒子は、第２の状態量が壁でないセル（すなわち、第２
の状態量が水面を示すセル）に、単位時間経過後、すな
わち、次の時間ステップで、等しい確率で必ず移動する
ものとし、また、ある粒子の運動は、他の粒子の運動に
影響されないと考えた。

【００２７】たとえば、図５（ａ）に示すように、ある
セルの粒子５０１に隣接する四方のセルの第２の状態量
が、すべて水面を示す（符号５０２〜５０５参照）であ
る場合には、当該粒子は、次の時間ステップで、１／４
の確率で、符号５０２〜５０５の位置のセルの何れかに
移動する。その一方、図５（ｂ）に示すように、あるセ
ルの粒子５０６に隣接する一つのセルの第２の状態量の
みが壁面を示す（符号５０７参照）場合には、当該粒子
は、次の時間ステップで、１／３の確率で、符号５０８
〜５１０の位置のセルの何れかに移動する。同様に、図
５（ｃ）および（ｄ）に示すように、あるセル５１１、
５１２の隣接する２つおよび３つのセルの第２の移動量
が壁面を示す場合には、当該粒子は、１／２、１／１の
確率で、それぞれ、残りの位置のセルに移動する。

【００２８】また、単位セル幅（セル間の距離）および
時間ステップ（シミュレーションを実行する際の単位時
間）を、実際の現象に対応するような物理量を割り当て
た。この場合、拡散係数Ｄは、（３）式のように表わす
ことができる。Ｄ＝（単位セル幅）^２／（２×（単位時間ステップ））・・・・（３）たとえば、セル幅を０．００１、時間ステップを０．０
００５として、セル上の粒子数を濃度に相当する量とし
て、粒子１個あたり０．００１の濃度が与えられる。

【００２９】このようにセルの状態量の初期値（境界条
件）、局所近傍系を決定してシミュレーションを実行す
ることにより、単位時間ごとの各セルの濃度（粒子数）
を得ることが可能となる。たとえば、後述するような入
力手順および処理手順にて、一次元の拡散現象のシミュ
レーションを実行した結果を図６に示す。図６において
も、セル幅を０．００１、時間ステップを０．０００５
としている。四角にてプロットされた点（たとえば符号
６０１参照）は、５００ステップだけ（つまり、ｔ＝
０．０２５）シミュレーションした結果、白丸にてプロ
ットされた点（たとえば符号６０２参照）は、３０００
ステップだけ（つまりｔ＝０．１５）シミュレーション
した結果、黒丸にてプロットされた点（たとえば符号５
０３参照）は、３００００ステップだけシミュレーショ
ンした結果を示す。また、図中の実線は、（１）式に示
した拡散方程式の解析解を示している。図６から、拡散
方程式の解析解と、本実施の形態にかかるＣＡによるシ
ミュレーション結果とは良く一致していることが理解で
きる。すなわち、ＣＡは一つの物理モデルであり、拡散
方程式で表現される数学モデルがもつ矛盾点を回避しな
がら、巨視的には簡単な規則を用いて現象を再現できる
という特徴を備えている。

【００３０】上述したような動作原理を有する拡散現象
シミュレータ１０の作動につき以下に説明する。図７
は、本実施の形態にかかる拡散現象シミュレータの必要
な設定処理を示すフローチャートである。操作者は、ま
ず、入力装置２０を操作して、解析空間となるセルの配
置を決定する（ステップ７０１）。入力装置２０からの
データが与えられるとセル状態量設定部３６は、セルの
配置を示すデータをセル管理部３２に与える。セル管理
部３２は、これに応答して、セル状態量記憶部３４中に
必要な領域を確保する。たとえば、上述した例では、操
作者が入力装置２０を操作することにより、２０×２０
００のセルの配置が決定される。

【００３１】次いで、操作者が入力装置２０を操作する
ことにより、各セルの状態量が与えられる（ステップ７
０２）。これに応答して、セル状態量設定部３６および
セル管理部３２が作動し、入力されたデータにしたがっ
て、各セルの第２の状態量および必要な場合には第１の
状態量が、セル状態量記憶部３４に、各セルに対応する
ように記憶される（ステップ７０２）。上述した例で
は、ｙ＝０およびｙ＝１のセルの第２の状態量が、壁面
を示すものに設定され、或いは、ｘ＝０およびｘ＝２０
００のセルにおける濃度（第１の状態量）も所定の値に
設定される。

【００３２】このようにして状態量が設定／記憶された
後に、シミュレーションに必要な経時的要素が入力され
る（ステップ７０４）。この経時的要素には、シミュレ
ーション時間、単位時間（ステップ）、粒子の供給条件
などを含む境界条件、図５を参照して説明した粒子の移
動確率などが含まれる。さらに、この経時的要素には、
後述するような粒子の速度（移流速度）ｖが含まれる。
この実施の形態において、粒子の移動確率および粒子の
速度（移流速度）などを、時間にしたがって変更するよ
うにしても良いし、固定値としても良い。

【００３３】ステップ７０４にて与えられたデータは、
シミュレーション条件設定／記憶部３８に伝達され、シ
ミュレーション条件設定／記憶部３８中のメモリ（図示
せず）に記憶され（ステップ７０５）、シミュレーショ
ン実行中に必要に応じてシミュレーション実行部４０に
出力される。このようにして、シミュレーション実行に
必要なセルの設定および各種要素の設定が終了する。

【００３４】次に、本実施の形態にかかる拡散現象シミ
ュレータにおけるシミュレーション処理につき、図８を
参照して説明を加える。操作者が入力装置２０を操作し
て、シミュレーションの開始を指示すると、これに応答
してシミュレーション実行部４０が起動されて、シミュ
レーションを開始する。シミュレーション実行部４０
は、セルを選択して、セル管理部３２を介して、セル状
態量記憶部３４から、対応する状態量などを示すデータ
を取得する（ステップ８０１）。次いで、シミュレーシ
ョン実行部４０は、図７のステップ７０４にて与えられ
た移流速度に基づいて、粒子の移動確率を算出する（ス
テップ８０２）。本実施の形態においては、単位時間ス
テップに１セルしか移動しないものと考えている。した
がって、移流速度ｖは、以下の（４）式にて与えられ
る。ｖ＝（単位セル幅）×ε／（単位時間ステップ）・・・・（４）

【００３５】ここで、εは、移動確率の差である。たと
えば、図９に示すように、矢印方向に速度ｖにて粒子が
流れると考える。この場合に、図９（ａ）のように、粒
子９０１が、四方（符号９０２〜９０５）に移動可能で
ある場合には、移動確率の差εは、 ε＝（Ｐ_２＋Ｐ_４）−（Ｐ_１＋Ｐ_３）・・・・（５）となる。図５に示すように、本来は（すなわち、速度が
ない場合には）、Ｐ_１＝Ｐ_２＝Ｐ_３＝Ｐ_４＝１／４であ
ったのが、移流速度ｖを考慮することで、εが算出さ
れ、これにより、Ｐ_１〜Ｐ_４が移流速度を考慮した値に
変更される。

【００３６】たとえば、単位セル幅が０．０５、単位時
間ステップが５×１０^−４であった場合に、移流速度ｖ
が１０と設定されていると、ε＝０．５となり、その結
果、Ｐ_１＝Ｐ_３＝１／８、Ｐ_２＝Ｐ_４＝３／８となる。

【００３７】また、図９（ｂ）に示すように、粒子９０
６が、三方（符号９０７〜９０９）に移動可能である場
合には、 ε＝Ｐ_４−（Ｐ_１＋Ｐ_３）・・・・（６）

【００３８】図９（ｃ）に示すように、粒子９１０が二
方（符号９１１、９１２）に移動可能である場合には、 ε＝Ｐ_２−Ｐ_１・・・・（７）と表わされる。なお、図９（ｄ）に示すように、粒子９
１３が移動可能な方向が一方向（符号９１４）に限定さ
れる場合には、移流速度が考慮されないことは言うまで
もない。

【００３９】このように、本実施の形態においては、移
流速度に基づいて移動確率の差εが算出され、これによ
り、速度を考慮した粒子の各方向への移動確率を得るこ
とができる。

【００４０】次いで、シミュレーション実行部４０は、
算出した粒子の移動確率にしたがって、粒子が移動すべ
きセルを決定する。これにより、このセルの状態量が変
化する。たとえば、上述したように、単位セル幅を０．
０５、単位時間ステップを５×１０^−４、移流速度ｖを
１０とした場合に、図９（ａ）の符号９０１に対応する
位置に、ｎ個の粒子が存在する、すなわち、対応するセ
ルの第２の状態量が、ｎ個の粒子を示すと考える。この
とき、ｎ×１／８（個）の粒子が、符号９０２、９０５
のセルにそれぞれ移動し、ｎ×３／８（個）の粒子が、
符号９０３、９０４のセルに移動する。したがって、こ
れに応じて、セルの状態量が変化することになる（ステ
ップ９０３）。

【００４１】このようにしてセルの状態量を変化させた
後に、シミュレーション実行部４０は、セル管理部３２
に、所定のセルの状態量をセル状態量記憶部３４に記憶
させるように指示するとともに、このセルに関するシミ
ュレーションの結果をシミュレーション結果記憶部４２
に記憶する（ステップ８０４）。さらに、シミュレーシ
ョン実行部４０は、シミュレーション結果を表示装置２
２やプリンタ２４に出力する。

【００４２】すべてのセルに関してステップ８０１〜８
０４の処理が終了すると（ステップ８０５でイエス(Ye
s)）、シミュレーション実行部４０は、シミュレーショ
ン時間が経過しているか否かを判断する（ステップ８０
６）。このステップにおいてノー(No)と判断された場合
には、シミュレーション実行部４０は、シミュレーショ
ン時間を単位時間ステップだけインクリメントし（ステ
ップ８０７）、必要な場合（たとえば、時間にしたがっ
て、経時的要素を変更する必要が有る場合）には、移流
速度や移動確率などを変更する。このようにして、シミ
ュレーション時間が経過するまで、ステップ８０１〜８
０５の処理が繰り返される。

【００４３】図１０に、流れの場の中にある濃度の高い
領域が時間とともに拡散する現象をシミュレートした例
を示す。図１０（ａ）に示すように、この例では、流れ
の下流側に固定構造物（壁面）（符号１０００参照）を
設置し、構造物の後方に流れが回り込む現象の再現を試
みている。この構造物は、油が拡散する際に設置される
オイルフェンスに相当すると考えることができる。

【００４４】この例においては、解析領域を１２０×１
１０のセルに分割し、初期状態として、１０００個の粒
子を、解析領域の略中央に位置する１２×１５のセル
（符号１００１参照）に配置した。また、単位セル幅ｄ
ｘ＝ｄｙ＝０．０５、単位時間ステップを５×１０^―４
とした。これにより、（３）式にしたがった等価な拡散
係数Ｄは２．５となる。さらに、矢印方向の移流速度ｖ
を１０としている。これにより、（４）式にしたがっ
て、移動確率の差εが０．５に設定されている。図１０
（ｂ）、（ｃ）、（ｄ）は、それぞれ、５０ステップ、
１００ステップ、２００ステップ経過後の拡散の状態を
示している。なお、本例において、３００ステップのシ
ミュレーションを実行するのに要した時間は、パーソナ
ルコンピュータ（ＣＰＵ：Pentium Pro (登録商標) 200
MHz、ＯＳ：Linux (登録商標)）を用いて５分程度であ
った。

【００４５】さらに、本発明においては、このシミュレ
ータを用いて、平成９年１月２日に日本海で発生したナ
ホトカ号の船体折損事故に伴う流出原油の拡散シミュレ
ーションを実行した。このシミュレーションでは、図１
１に示すように、日本海の事故現場周辺約６００ｋｍ四
方の範囲（符号１１０１参照）を解析平面とした。この
平面を、一辺約５ｋｍの正方形のセルに分割した。これ
により、解析平面は、１２０セル（東西方向）×１１０
セル（南北方向）に分割されている。また、流出原油の
量は、粒子数にて表現し、陸地を壁（すなわち、セルの
第２の状態量が壁面を示す）と考えた。本例において
は、重油に相当する粒子は、壁の前面にて堆積され、堆
積された粒子数は、海岸に打ち上げられた重油の量に相
当すると考えた。すなわち、セルの第２の状態量が壁面
を示す場合に、そのセルに対しても粒子が移動できるも
のとして、かつ、移動した粒子が、当該セルに堆積され
るように設定した。

【００４６】移流として、潮流および吹送流の両者を考
慮し、（４）式に関して説明したように、粒子が隣接す
るセルに移動する確率に差を持たせることで、移流拡散
現象を計算の中で表現した。なお、潮流のデータは、海
上保安庁の日本近海海流統計図中の１〜３月の平均潮流
速度に基づき算出した。また、風向および風速は、気象
庁が日本海に設置している観測ブイ（東経１３４度、北
緯３７度）により記録された３時間ごとのデータに基づ
き与えた。折損事故現場は、鳥取県隠岐島の沖北北東約
１０６ｋｍの海上（東経１３３度５２分、北緯３７度１
０分）であったので、気象庁の観測ブイは、事故現場の
比較的近くに位置する。吹送流を計算するために用いた
風力係数は資料を参考にして１０^−２程度とした。な
お、事故当日以降の気象庁の天気図より事故現場付近の
風向および風速を計算することは困難なため、第１次近
似として、発生現場の海域の４セルに各々１０００個の
粒子を配置し、１〜１０時間ごとの拡散状況のシミュレ
ーションを実行した。今回の事故で流出した重油の拡散
係数は明らかではなく、また、拡散係数は、気温、海水
温による大きく異なることが推察されるので、シミュレ
ーションでは、風力係数と拡散係数に相当するパラメー
タを変化させて計算を行って、計算結果と公表されてい
る重油漂着のデータとを比較することとした。

【００４７】図１２（ａ）〜（ｄ）は、風力係数を０と
して、拡散係数を、それぞれ、１．２５、２．５、６．
２５、１２．５（ｋｍ^２／ｈ）として、事故から１２日
後の重油の拡散をシミュレーションした結果を示す図で
ある。図１３（ａ）〜（ｄ）は、風力係数を０．０１と
して、拡散係数を図１２と同様に設定したものであり、
図１４（ａ）〜（ｄ）は、風力係数を０．０２とした場
合のものである。図１２〜図１４から、風力係数を大き
くするのにしたがって、油の濃度が高い領域は、北東方
向に大きく移動し、能登半島を巻き込むように佐渡、さ
らには、新潟まで到達する可能性が有ることが理解でき
る。

【００４８】図１５（ａ）は、一般に公表された重油漂
着の範囲およびその日付を示す図である。図１２ないし
図１４と、図１５（ａ）とを比較して、拡散係数を２．
５ｋｍ^２／ｈ、風力係数を０．０１とすると、図１５
（ｂ）のようなシミュレーション結果を得ることができ
る。図１５（ａ）と図１５（ｂ）とから、新潟方面の重
油の漂着をかなり正確に予測できたことが理解できる。
その一方、図１５（ｂ）からは、若狭湾から福井県にか
けての沿岸への重油の漂着は予測できなかったが、これ
は、風力および風向を解析領域内にて一定にしているこ
と、重油の流出箇所が事故現場一個所に限定されたこと
などによるものと思われる。なお、風の影響を考慮しな
い場合には、図１５（ｃ）に示すようなシミュレーショ
ン結果となり、若狭湾への実際の漂着と略同等の結果が
得られていることが理解できる。

【００４９】折損した船首部の漂流およびそこからの流
出油を考慮したシミュレーション結果の例を図１６〜図
１８に示す。図１６では、沈没地点において油が初期的
に流出し（すなわち、初期流出のみが生じ）、漂流する
船首部からの油の流出が次第に減少するという設定の下
でシミュレーションを実行している。図１７では、沈没
地点においては油が一定量流出し続け、漂流する船首部
からの油の流出が次第に減少するという設定がなされ、
その一方、図１８では、沈没地点および漂流する船首部
の双方から一定量の油が流出し続けるという設定がなさ
れている。また、各図において、それぞれ（ａ）は３日
後、（ｂ）は６日後、（ｃ）は１２日後の油の拡散状況
を示している。図１６〜図１８に示すように、油の流出
箇所を変更（すなわち、シミュレーションにおいて油の
粒子を供給すべきセルおよび供給する油の量を経時的に
変更）することにより、重油の漂着をより正確に推定す
ることが可能となることが理解できる。

【００５０】たとえば、上述したように一定の風向およ
び風速で、かつ、一定の拡散係数の下でのシミュレーシ
ョン計算は、通常のパーソナルコンピュータを用いて１
〜１０分程度で実行することができる。したがって、重
油や汚染物質の流出など緊急を要する事故後の対策を立
てるために、迅速にシミュレーションを実行することが
可能である。また、本発明においては、局所近傍則のみ
を用いて現象を再現しているため、従来の支配方程式を
解く手法のように境界条件に拘束されることなく、事故
周辺海域など必要な領域を取り出して解析空間として、
これを解析すれば良く、必要な領域の現象のみを再現す
ることが可能となる。

【００５１】このように、本実施の形態によれば、拡散
現象を解析すべき空間を複数のセルに分割し、かつ、セ
ルに存在する一以上の粒子を局所近傍系を構成する他の
セルに移動させることにより、粒子の拡散をシミュレー
トしている。したがって、複雑な方程式の解析解を求め
ることなく、適切に、拡散現象をシミュレートすること
が可能となる。

【００５２】次に、本発明の第２の実施の形態につき説
明を加える。第２の実施の形態においては、第１の実施
の形態にかかる拡散現象シミュレータに、風速、風向
や、潮流を検出する手段を付加し、重油や汚染物質の流
出現場などで、より正確に、重油や汚染物質の拡散を可
能としている。

【００５３】図１９は、第２の実施の形態にかかる拡散
現象解析システムの機能を示すブロックダイヤグラムで
ある。図１９において、図２に示すシミュレータと同一
の構成部分には同一の符号を付し、また、この拡散現象
解析システムは、図１に示すハードウェア構成にて実現
可能である。図１９に示すように、このシステム５０
は、拡散現象シミュレーションを構成するセル管理部３
２、シミュレーション実行部４０などの他、風向および
風速を示すデータに基づき、水面上に発生する吹送流を
算出する吹送流算出部５２と、統計的な潮流の速度およ
び方向を示すデータを予め記憶した潮流データベース５
４と、吹送流および潮流に基づき、粒子の移動に偏りを
持たせるための移動確率の差を算出する速度算出部５６
とを備えている。

【００５４】吹送流算出部５２には、操作者が入力装置
２０を操作することにより、或いは、吹送流算出部５２
に接続されたセンサにより、風向および風速を示す観測
値が与えられる。吹送流算出部５２は、この観測値に風
力係数を乗じて流れの速度（ベクトル）を推定する。な
お、風力係数は実験値を入力装置２０に入力すれば良
い。

【００５５】速度算出部５６は、日時や現場の位置に基
づき、潮流データベースから、必要な潮流の方向および
大きさを示すデータを読み出すとともに、吹送流算出部
からの流れの速度ベクトルを受け入れ、これらのベクト
ル和を求めて、速度ｖの大きさおよび方向を決定する。
このようにして、速度算出部５６にて算出された速度ｖ
は、シミュレーション条件設定／記憶部３８に一時的に
記憶され、実際のシミュレーション中に、シミュレーシ
ョン実行部４０にて、移動確率の差ε（（４）式参照）
を求めるために用いられ、その結果、各セルの粒子の移
動確率を決定するために用いられる。なお、実際のシミ
ュレーションの処理は、図８のフローチャートに示すも
のと略同様であるため、その説明は省略する。

【００５６】本実施の形態によれば、実際の風力の観測
値や潮流を加味して速度を算出し、この速度に基づき、
粒子の移動に偏りを持たせることができる。また、パー
ソナルコンピュータにて、システムを実現することがで
きる。したがって、重油や汚染物質流出の現場などで、
シミュレーションを実行することで、リアルタイムで拡
散状況を、高精度で予測ないし推定することが可能とな
る。また、現場での拡散状況の調査結果に基づき拡散の
度合いを簡単に変更することができるため、拡散状況
を、さらに高精度で予測なしし推定することが可能とな
る。

【００５７】本発明は、以上の実施の形態に限定される
ことなく、特許請求の範囲に記載された発明の範囲内
で、種々の変更が可能であり、それらも本発明の範囲内
に包含されるものであることは言うまでもない。たとえ
ば、前記実施の形態においては、図４に示すように、各
セル（座標）に対して、壁面／水面を示す第２の状態量
と、粒子数を示す第の状態量とを定義したが、これに限
定されるものではなく、他の状態量を割り当てても良
い。たとえば、図２０に示すように、各セルに対して、
粒子の移動確率を定義しておいても良い。この移動確率
は、セルの位置に応じて、操作者が入力装置２０を操作
して設定しても良いし、図１９に示す速度算出部５６に
て算出された速度などに基づき、自動的に設定しても良
い。

【００５８】また、前記実施の形態においては、図５に
示すように、あるセルの粒子は、最大で隣接する４方向
に移動可能なように設定されているが、粒子の移動方向
はこれに限定されるものではなく、たとえば、図２１
（ａ）、（ｂ）に示すように、平面上で粒子が３方向、
６方向などに移動できるようにしても良い。なお、図２
１（ｃ）に示すように、正方形のセルにおいて、粒子が
８方向に移動できるように設定しても良いが、粒子の移
動距離を等しい方が好ましいため、図５や図２１
（ａ）、（ｂ）のような設定の方がより好ましい。

【００５９】さらに、前記実施の形態においては、二次
元の拡散現象をシミュレートするために、平面上の解析
空間をセルにて分割しているが、本発明は、三次元の拡
散現象をシミュレートすることも可能である。たとえ
ば、この場合には、図２２に示すように、三次元の解析
空間を立方体のセルにて分割し、セル（符号２２０１参
照）の粒子が、矢印に示す８方向に移動できるように設
定すれば良い。この場合も、セルの形状は立方体に限定
されるものではなく、また、粒子の移動方向を８方向以
外に設定しても良い。

【００６０】また、前記実施の形態においては、重油や
汚染物質など液体の拡散現象を解析したが、これに限定
されるものではなく、たとえば、工場や焼却炉から排出
される煙などの拡散や、気体の拡散などにも本発明を適
用できることは言うまでもない。

【００６１】さらに、本明細書において、手段とは必ず
しも物理的手段を意味するものではなく、各手段の機能
が、ソフトウェアによって実現される場合も包含する。
さらに、一つの手段の機能が、二つ以上の物理的手段に
より実現されても、若しくは、二つ以上の手段の機能
が、一つの物理的手段により実現されてもよい。

【００６２】

【発明の効果】本発明によれば、複雑な処理や複雑な方
程式の解法を用いることなく、短時間で精度良く拡散現
象を解析するシミュレーションシステムを提供すること
が可能となる。特に、小型コンピュータを用いて、短時
間でかつ高精度に油、汚染物質などの拡散シミュレーシ
ョンを実行することができ、現場などにて観測値を入力
して、迅速に拡散状況の予測、推定をなすことが可能と
なる。

【図面の簡単な説明】

【図１】図１は、本発明の第１の実施の形態にかかる
拡散現象シミュレータのハードウェア構成を示すブロッ
クダイヤグラムである。

【図２】図２は、第１の実施の形態にかかる拡散現象
シミュレータの機能を示すブロックダイヤグラムであ
る。

【図３】図３は、第１の実施の形態にかかる解析平面
中のセルの構成を説明するための図である。

【図４】図４は、第１の実施の形態にかかるセルに割
り当てられる種々のデータを説明するための図である。

【図５】図５は、セル中の粒子の移動を説明するため
の図である。

【図６】図６は、第１の実施の形態において、一次元
の拡散現象を、二次元のＣＡを用いてシミュレートした
結果を示すグラフである。

【図７】図７は、第１の実施の形態にかかる拡散現象
シミュレータにおいて、シミュレーションに必要な設定
の処理を示すフローチャートである。

【図８】図７は、第１の実施の形態にかかる拡散現象
シミュレータのシミュレーション処理を示すフローチャ
ートである。

【図９】図９は、速度を考慮した粒子の移動確率を説
明するための図である。

【図１０】図１０は、流れの場の中にある濃度の高い
領域が時間とともに拡散する現象をシミュレートした例
を示す図である。

【図１１】図１１は、シミュレーションのために用い
た解析空間を示す図である。

【図１２】図１２は、図１１に示す解析空間におい
て、風力係数を０としてシミュレーションを実行した結
果を示す図である。

【図１３】図１３は、図１１に示す解析空間におい
て、風力係数を０．０１としてシミュレーションを実行
した結果を示す図である。

【図１４】図１４は、図１１に示す解析空間におい
て、風力係数を０．０２としてシミュレーションを実行
した結果を示す図である。

【図１５】図１５は、実際の重油の漂着と、シミュレ
ーションによる重油の移動および漂着の予測とを比較し
た図である。

【図１６】図１６は、油の流出位置および流出量の変
動を考慮してシミュレーションを実行した結果を示す図
である。

【図１７】図１７は、油の流出位置および流出量の変
動を考慮してシミュレーションを実行した結果を示す図
である。

【図１８】図１８は、油の流出位置および流出量の変
動を考慮してシミュレーションを実行した結果を示す図
である。

【図１９】図１９は、第２の実施の形態にかかる拡散
現象解析システムの機能を示すブロックダイヤグラムで
ある。

【図２０】図２０は、状態量としてセルに割り当てら
れる種々のデータの他の例を説明するための図である。

【図２１】図２１は、セルおよび粒子の移動方向の他
の例を示す図である。

【図２２】図２２は、セルおよび粒子の移動方向のさ
らに他の例を示す図である。

【符号の説明】

１０拡散現象シミュレータ１２ＣＰＵ１４メモリ１６固定記憶装置１８インタフェース２０入力装置２２表示装置２４プリンタ３２セル管理部３４セル状態量記憶部３６セル状態量設定部３８シミュレーション条件設定／記憶部４０シミュレーション実行部４２シミュレーション結果記憶部

Claims

【特許請求の範囲】

【請求項１】解析対象となる空間を複数に分割してセ
ルを生成するセル生成手段と、各セルに対応して、当該セルに対応する領域に含まれる
粒子の量を示す第１の状態量を定義する状態量定義手段
と、単位時間ごとの、セルの周囲に位置する局所近傍を構成
する他のセルへの第１の状態量の移動を定義する移動量
定義手段と、前記定義された他のセルへの状態量の移動に基づき、前
記セルごとに、前記定義された第１の状態量を変更する
ことにより、粒子の移動を実現するシミュレーション手
段とを備えたことを特徴とする拡散現象シミュレータ。
【請求項２】前記第１の状態量が、対応するセルの領
域に含まれる粒子の数にて表わされることを特徴とする
請求項１に記載の拡散現象シミュレータ。
【請求項３】前記シミュレーション手段が、単位時間
ごとに、セル中の粒子を、当該セルに隣接する局所近傍
を構成するセルに移動させるように構成され、前記セル
の大きさおよび前記単位時間により、拡散係数が決定さ
れることを特徴とする請求項１または２に記載の拡散現
象シミュレータ。
【請求項４】前記状態量定義手段が、前記セルに対応
して、粒子が通過し得るか否かを示す第２の状態量を定
義することを特徴とする請求項１ないし３の何れか一項
に記載の拡散現象シミュレータ。
【請求項５】前記第２の状態量が、粒子が通過し得る
領域であること、および、障害物であることの何れか一
方を示すことを特徴とする請求項４に記載の拡散現象シ
ミュレータ。
【請求項６】前記移動量定義手段が、前記第２の状態
量に基づき、局所近傍を構成するセルのうち、粒子を移
動させることが可能なセルを選択することを特徴とする
請求項４または５に記載の拡散現象シミュレータ。
【請求項７】前記シミュレーション手段が、前記第２
の状態量が障害物であることを示すセルにおいて、粒子
を蓄積させるように作動することを特徴とする請求項４
または５に記載の拡散現象シミュレータ。
【請求項８】前記移動量定義手段が、前記局所近傍を
構成するセルの各々に粒子を移動させる移動確率を決定
することを特徴とする請求項１ないし７の何れか一項に
記載の拡散現象シミュレータ。
【請求項９】前記移動量定義手段が、移流速度を考慮
して、前記移動確率を決定することを特徴とする請求項
８に記載の拡散現象シミュレータ。
【請求項１０】前記移動量定義手段が、移流方向に向
かう粒子の移動確率と、その反対方向に向かう粒子の移
動確率との差が、前記移流速度と比例するように、前記
粒子の移動確率を決定することを特徴とする請求項９に
記載の拡散現象シミュレータ。
【請求項１１】請求項９または１０に記載の拡散現象
シミュレータと、風力および風向を示すデータに基づ
き、移流速度を算出する速度算出手段とを備え、前記速
度算出手段にて算出された移流速度に基づき、移動確率
が決定されることを特徴とする拡散現象解析システム。
【請求項１２】前記速度算出手段が、前記風力および
風向を示すデータと、潮流の大きさおよび向きを示す潮
流データとに基づき、海上に浮遊する粒子の移流速度を
算出するように構成されたことを特徴とする請求項１１
に記載の拡散現象解析システム。
【請求項１３】解析対象となる空間を複数に分割して
セルを生成するステップと、各セルに対応して、当該セルに対応する領域に含まれる
粒子の量を示す第１の状態量を定義ステップと、単位時間ごとの、セルの周囲に位置する局所近傍を構成
する他のセルへの第１の状態量の移動を定義するステッ
プと、前記定義された他のセルへの状態量の移動に基づき、前
記セルごとに、前記定義された第１の状態量を変更する
ことにより、粒子の移動を実現するステップとを備えた
ことを特徴とする拡散現象のシミュレーション方法。
【請求項１４】前記第１の状態量が、対応するセルの
領域に含まれる粒子の数にて表わされることを特徴とす
る請求項１３に記載の拡散現象のシミュレーション方
法。
【請求項１５】前記シミュレーションのステップが、
単位時間ごとに、セル中の粒子を、当該セルに隣接する
局所近傍を構成するセルに移動させるように構成され、
前記セルの大きさおよび前記単位時間により、拡散係数
が決定されることを特徴とする請求項１３または１４に
記載の拡散現象のシミュレーション方法。
【請求項１６】前記状態量を定義するステップが、前
記セルに対応して、粒子が通過し得るか否かを示す第２
の状態量を定義するステップを含むことを特徴とする請
求項１３ないし１５の何れか一項に記載の拡散現象のシ
ミュレーション方法。
【請求項１７】前記移動量を定義するステップが、前
記局所近傍を構成するセルの各々に粒子を移動させる移
動確率を決定するステップを含むことを特徴とする請求
項１３ないし１６の何れか一項に記載の拡散現象のシミ
ュレーション方法。
【請求項１８】前記移動確率を決定するステップが、
移流速度を考慮して、移動確率を決定することを特徴と
する請求項１７に記載の拡散現象のシミュレーション方
法。
【請求項１９】拡散現象をシミュレートするプログラ
ムを記憶した、コンピュータにより読出し可能な記憶媒
体であって、解析対象となる空間を複数に分割してセルを生成するス
テップと、各セルに対応して、当該セルに対応する領域に含まれる
粒子の量を示す第１の状態量を定義ステップと、単位時間ごとの、セルの周囲に位置する局所近傍を構成
する他のセルへの第１の状態量の移動を定義するステッ
プと、前記定義された他のセルへの状態量の移動に基づき、前
記セルごとに、前記定義された第１の状態量を変更する
ことにより、粒子の移動を実現するステップとを備えた
プログラムを格納したことを特徴とするコンピュータに
より読み出し可能な記憶媒体。
【請求項２０】前記第１の状態量が、対応するセルの
領域に含まれる粒子の数にて表わされることを特徴とす
る請求項１９に記載のコンピュータにより読み出し可能
な記憶媒体。
【請求項２１】前記シミュレーションのステップが、
単位時間ごとに、セル中の粒子を、当該セルに隣接する
局所近傍を構成するセルに移動させるように構成され、
前記セルの大きさおよび前記単位時間により、拡散係数
が決定されることを特徴とする請求項１９または２０に
記載のコンピュータにより読み出し可能な記憶媒体。
【請求項２２】前記状態量を定義するステップが、前
記セルに対応して、粒子が通過し得るか否かを示す第２
の状態量を定義するステップを含むことを特徴とする請
求項１９ないし２１の何れか一項に記載のコンピュータ
により読み出し可能な記憶媒体。
【請求項２３】前記移動量を定義するステップが、前
記局所近傍を構成するセルの各々に粒子を移動させる移
動確率を決定するステップを含むことを特徴とする請求
項１９ないし２４の何れか一項に記載のコンピュータに
より読み出し可能な記憶媒体。
【請求項２４】前記移動確率を決定するステップが、
移流速度を考慮して、移動確率を決定することを特徴と
する請求項２３に記載のコンピュータにより読み出し可
能な記憶媒体。