JP2014081900A - シミュレーションプログラム、シミュレーション方法及びシミュレーション装置 - Google Patents

Abstract

【課題】大きな範囲における連続体の移動を精度良くシミュレートするシミュレーションプログラム、シミュレーション方法及びシミュレーション装置を提供する。
【解決手段】入力されたデータを基に、第１領域の２次元平面の各点における連続体の流量及び水位を算出し、前記第１領域に含まれる第２領域における前記連続体を粒子の集合として表現し、算出した前記流量及び前記水位を基に、各前記粒子の状態を３次元解析する処理をコンピュータに実行させる。
【選択図】図３

Description

本発明は、シミュレーションプログラム、シミュレーション方法及びシミュレーション装置に関する。

従来、数値計算を用いて水や空気といった流体の流れを調べる連続体の運動の解析、すなわち、流体解析において、粒子法と呼ばれる技術が提案されている。具体的には、粒子法とは、連続体の運動を有限の数の粒子の運動として解析する手法である。現在提案されている代表的な粒子法としては、ＳＰＨ（Smoothed Particles Hydrodynamics）法やＭＰＳ（Moving Particles Semi-implicit）法といったものがある。以下では、水や空気といった流体を「連続体」と呼ぶ場合がある（例えば、非特許文献１参照）。

粒子法における標準的な手法としては、ある粒子に対して予め領域（以下では、「影響領域」と呼ぶ。）を設定しておき、その影響領域内に存在する相手粒子のみからの相互作用を計算することで、その粒子に加えられる力を求める方法が従来知られている。

ＳＰＨ法を用いて連続体を表現する際の特徴は、複数の粒子の物理量をカーネル関数と呼ばれる重み関数を用いて平滑化近似して、基礎方程式を離散化することにある。この平滑化処理によって、メッシュ点上で物理量を評価するといったきわめて煩雑なメッシュ操作の計算処理を行わなくてよくなる。そのため、ＳＰＨ法は、自由表面問題やマルチフィジックス問題を扱うのに適しているといえる。具体的には、ＳＰＨ法は、例えば、海の波が護岸に衝突や越波するときの流量や衝撃圧を推定するための手法として適しているといえる。

ただし、ＳＰＨ法による計算を用いて広範囲の計算を行った場合、計算量が大きくなりすぎる。そのため、ＳＰＨ法を用いた場合には、大きくとも数ｋｍ四方程度の限られた領域でのみの計算を行うことになる。そのため、津波の伝播を考えた場合、津波計算に用いる震源直上の津波波源から沿岸部まで波が伝播する数１００ｋｍ以上に広がる過程をＳＰＨ法で扱うことは困難である。

また、地形データ及び構造データと組み合わせて流体のシミュレーションを行う従来技術（例えば、特許文献１参照）があるが、この方法を用いても数１００ｋｍ以上といった広域を扱うことは困難である。

一方、地震による地殻変動から、津波発生を経て、震源から海岸に至る津波の伝播までの過程を２次元の伝播シミュレータによって計算する従来技術も提案されている（例えば、特許文献２、３、非特許文献２参照）。これらの手法は、２次元である地図上の各点で波高及び速度等の情報が設定されており、津波が発生した場合の沿岸における波高及び浸水域等を数値計算によって予測することが可能である。

なお、低い次元で計算した結果を高い次元の計算の初期条件及び境界条件として与えるという従来技術が提案されている。例えば、加震源と受振構造物を三次元有限要素法で扱い、その間の地盤を介して伝わる振動を簡易的なモデルで求めるといった従来技術がある（例えば、特許文献４参照）。また、１次元のBoussinesq波動モデルと２次元のＳＰＨ法の計算を組み合わせる従来技術がある（例えば、非特許文献３参照）。１次元のBoussinesq波動モデルと２次元のＳＰＨ法の計算を組み合わせる従来技術では、例えば、以下のような処理が行われる。ＳＰＨ法で計算を行う領域とBoussinesq法で計算を行う領域とが重なり合う領域であるバッファ領域を設定する。バッファ領域の中では、Boussinesq法の計算結果として得られた速度を強制的にＳＰＨ法の各粒子に設定する。また、バッファ領域の外においては、ＳＰＨ粒子は流体の運動方程式を解いた結果として得られる加速度に基づいて、速度の更新が行われる。ただし、バッファ領域は予め決められており、ＳＰＨ粒子の生成や消滅は行われない。

特開２００５−１１５７０１号公報 特開２００８−０８９３１６号公報 特開２０１０−０５４４６０号公報 特開２００４−２１９２３７号公報

J.J.Monaghan 著 「Smoothed Particle Hydrodynamics」 Annu. Rev. Astron. Astrophys., Vol.30, pp.543-574 Coto, C., Ogawa, Y., Shuto, N., and Imamura, F.: Numerical method of tsunami simulation with the leap-frog scheme (IUGG/IOC Time Project), IOC Manual, UNESCO, No.35, 1997, pp.6-16 Kassiotis et al. 2011, 「Coupling SPH with a 1-D Boussinesq-type wave model」 http://hal.inria.fr/docs/00/60/13/52/PDF/KassiotisEtAl＿Spheric6.pdf

しかしながら、２次元シミュレーションのみを用いた従来技術では、３次元性が重要になる港湾部などの複雑形状を扱うことが困難となることや、建築物にかかる波圧及び波力を算出することが困難であるといった問題がある。

また、加震源と受振構造物を三次元有限要素法で扱い、その間の地盤を介して伝わる振動を簡易的なモデルで求めるといった従来技術は、地盤中の振動に適用されるものであり、そのまま津波シミュレーションに用いることは困難である。

また、Boussinesq波動モデルをＳＰＨ法の計算に用いる従来技術では、津波のように大量の流体が移動を伴う場合には、バッファ領域が狭いと、バッファ領域の粒子が全てＳＰＨ領域に移動してしまい、バッファ領域での計算ができなくなるおそれがある。一方、津波の流体移動を考慮しバッファ領域の幅を決めると、バッファ領域が非常に大きくなり、計算量及びメモリ量が増大しすぎてしまい、実際に処理することが困難となる。さらに、この従来技術では、１次元と２次元との組合せを対象としており、２次元と３次元とを組み合わせることは考慮されておらず、津波の伝播計算にそのまま用いることは困難である。

開示の技術は、上記に鑑みてなされたものであって、連続体の移動を精度良くシミュレートするシミュレーションプログラム、シミュレーション方法及びシミュレーション装置を提供することを目的とする。

本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置は、一つの態様において、入力されたデータを基に、第１領域の２次元平面の各点における連続体の流量及び水位を算出し、前記第１領域に含まれる第２領域における前記連続体を粒子の集合として表現し、算出した前記水位及び前記流量を基に、各前記粒子の状態を３次元解析する処理をコンピュータに実行させる。

本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置の一つの態様によれば、連続体の移動を精度良くシミュレートすることができるという効果を奏する。

図１は、実施例１に係るシミュレーション装置のブロック図である。 図２は、２次元伝播計算及び３次元粒子法計算を行う領域の例を示す図である。 図３は、実施例１に係るシミュレーション装置による解析処理のフローチャートである。 図４は、実施例２に係る２次元伝播計算処理部のブロック図である。 図５は、実施例２に係る３次元ＳＰＨ計算処理部のブロック図である。 図６は、初期値の設定を模式的に表した図である。 図７は、影響範囲及び近傍粒子について説明するための図である。 図８は、粒子の流出を模式的に表した図である。 図９は、境界条件算出における補間を説明するための図である。 図１０は、境界からの反発力を説明するための図である。 図１１は、粒子の流入を模式的に表した図である。 図１２は、実施例２に係る２次元伝播計算の処理のフローチャートである。 図１３は、実施例２に係る３次元ＳＰＨ計算の処理のフローチャートである。 図１４は、実施例２に係る３次元ＳＰＨ計算の処理のフローチャートである。 図１５は、実施例３に係るシミュレーション装置による解析処理のフローチャートである。 図１６は、シミュレーションプログラムを実行するコンピュータを示す図である。

以下に、本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置の実施例を図面に基づいて詳細に説明する。なお、以下の実施例により本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置が限定されるものではない。また、以下では、津波発生時における海水の動きを解析するシミュレーションについて説明するが、連続体の動きを解析するシミュレーションであれば各実施例を適応できる。特に、津波のような連続体の移動量が大きい場合のシミュレーションに適応することが好ましい。

図１は、実施例１に係るシミュレーション装置のブロック図である。図１に示すように、本実施例に係るシミュレーション装置は、ユーザインタフェース１、２次元伝播計算処理部２、３次元ＳＰＨ計算処理部３及び結果記憶部４を有している。

ユーザインタフェース１は、シミュレーション装置を使用する操作者からの数値の入力、操作者へのシミュレーション結果の出力を行なうインタフェース装置であり、具体的にはキーボード等の入力装置及び表示装置等の出力装置である。

操作者は、ユーザインタフェース１を用いて、例えば、図２に示すように領域１００及び領域１０１を設定する。図２は、２次元伝播計算及び３次元粒子法計算を行う領域を示す図である。領域１０１は、後述する２次元伝播計算を行う領域である。領域１０１が、「第１領域」の一例にあたる。領域１００は、後述する３次元粒子法計算を行う領域である。領域１００が、「第２領域」の一例にあたる。このように、２次元平面の各点における連続体の流量及び水位を算出する２次元伝播計算は、大きな領域に対して行い、３次元ＳＰＨ計算は、沿岸部などの３次元的に入り組んだ構造を含む小さな領域に対して行う。ここで、領域１００と領域１０１との境界はある程度水深がある場所であることが好ましい。これは、水深が浅くなると、海水の流れが海底や沿岸の影響を受け易くなるため２次元伝播計算の精度が落ちるからである。２次元伝播計算の計算結果を３次元粒子法計算に用いるためには、２次元伝播計算の精度が、３次元ＳＰＨ計算で計算を行った場合の精度と同程度であることが好ましい。そこで、ある程度の深さがある地点を境界として設定することが好ましい。具体的には、例えば、３０ｍ〜５０ｍの水深を有する地点を境界として設定することが好ましい。

さらに、操作者は、初期データの入力を行う。本実施例では、図２の領域１０１を平面視した２次元平面を所定の格子に分割する。２次元伝播計算では、この分割した格子における海水の動きを解析する。そして、初期データとは、領域１０１を表す座標をｘｙ座標として、分割した各格子における、水位η、ｘ方向の流量Ｍ（以下では、単に「流量Ｍ」という場合がある。）及びｙ方向の流量Ｎ（以下では、単に「流量Ｎ」という場合がある。）の初期値である。ここで、水位ηとは、静水面からの水位上昇量である。また、流量とは、速度と全水深との積であり、ｘ方向の流量Ｍは、Ｍ＝ｖ_ｘ×Ｄ（ｖ_ｘはｘ方向の速度、Ｄは全水深であり、Ｄ＝水位η＋静水深ｈである。）として表される。また、ｙ方向の流量Ｎは、Ｎ＝ｖ_ｙ×Ｄ（ｖ_ｙはｙ方向の速度である。）として表される。全水深と流量から速度を得ることができ、また、逆に全水深と速度から流量を得ることもできるので、流量と速度のどちらも本質的には同等の情報を含んでいるとみなせる。二次元伝播計算では流量の方が扱いやすいため、流量を用いるが、三次元粒子法計算では速度の方が扱いやすいため、流量の情報を必要に応じて速度に変換することになる。速度ｖ_ｘと速度ｖ_ｙとを合成した速度が各格子における流速になる。

そして、２次元伝播計算処理部２は、入力された水位η、流量Ｍ及び流量Ｎの初期値を用いて、次のステップの各格子における水位η、流量Ｍ及び流量Ｎを算出する。さらに、２次元伝播計算処理部２は、順々に各ステップの各格子における水位η、流量Ｍ及び流量Ｎを算出していく。２次元伝播計算処理部２は、予め設定された最終の状態になる最終ステップまでの各ステップの各格子における水位η、流量Ｍ及び流量Ｎの算出を繰り返す。この予め設定された最終状態になったときのステップが、「第１の所定ステップ」にあたる。本実施例では、津波が陸上の最深部まで達し引き始めた時点を最終状態とする。

そして、２次元伝播計算処理部２は、算出したステップ毎の各格子における水位η、流量Ｍ及び流量Ｎのデータを３次元ＳＰＨ計算処理部３へ出力する。

３次元ＳＰＨ計算処理部３は、ステップ毎の各格子における水位η、流量Ｍ及び流量Ｎのデータの入力を２次元伝播計算処理部２から受ける。

そして、３次元ＳＰＨ計算処理部３は、受信したステップ毎の各格子における水位η、流量Ｍ及び流量Ｎのデータのうち、３次元粒子法計算を開始する条件を満たすデータを有する最も時刻の早いステップを特定する。３次元ＳＰＨ計算処理部３は、例えば、図２の領域１００の境界における水位ηが５０ｃｍを超えた時点を、３次元粒子法計算を開始する条件として記憶している。

そして、３次元ＳＰＨ計算処理部３は、特定した３次元粒子法計算を開始するステップにおける各格子における水位ηを用いて領域１００のその時点での水深に合わせて粒子を配置する。ここで、水深とは水位ηで表される頂点から海底までの距離である。また、粒子は、どのような大きさにすることもできるが、直径が１０ｃｍ〜１ｍであることが好ましい。直径が小さいほど精度の良い解析ができるが計算量が増え、直径が大きいほど計算量が減るが解析の精度が落ちる。そのため、粒子の大きさは、使用できる計算資源と供給される精度に応じて決定することが好ましい。粒子は、流体をモデル化したものであり、粒子を表現するために要求されるデータは、例えば、中心座標、速度、影響半径、密度、質量、粘性度などである。このうち、３次元ＳＰＨ計算処理部３は、影響半径、質量、粘性度などを予め記憶している。

さらに、３次元ＳＰＨ計算処理部３は、３次元粒子法計算を開始するステップの各格子における流量Ｍ及び流量Ｎから各粒子の速度の初期値である速度ｖ_０を求める。ここで、「ｖ」は粒子の速度を表すベクトルである。また、３次元ＳＰＨ計算処理部３は、各格子から各粒子の位置座標の初期値である位置ｒ_０を求める。ここで、「ｒ」は粒子の位置を表す位置ベクトルである。さらに、３次元ＳＰＨ計算処理部３は、圧力及び密度の初期値である圧力Ｐ_０及び密度ρ_０として海水の一般的な圧力及び密度をそれぞれ用いる。

そして、３次元ＳＰＨ計算処理部３は、速度ｖ_０、位置ｒ_０、圧力Ｐ_０及び密度ρ_０を用いて、次のステップの各粒子の速度ｖ、位置ｒ、圧力Ｐ及び密度ρを求める。このとき、３次元ＳＰＨ計算処理部３は、図２の領域１００における境界から一定の距離Δｘの範囲を境界領域として記憶している。そして、３次元ＳＰＨ計算処理部３は、各ステップにおける速度ｖ、位置ｒ、圧力Ｐ及び密度ρを算出する場合に、境界領域における粒子の速度ｖなどの境界条件を、受信したステップ毎の各格子における流量Ｍ及び流量Ｎから求める。ここで、３次元粒子法計算を行うステップは、２次元伝播計算を行うステップよりも細かい。そこで、３次元ＳＰＨ計算処理部３は、境界条件として用いる流量Ｍ及び流量Ｎを補間しながら境界条件を求めていく。

このように、３次元ＳＰＨ計算処理部３は、予め設定された最終の状態になるまでの各ステップの各粒子における速度ｖ、位置ｒ、圧力Ｐ及び密度ρの算出を繰り返す。この予め設定された最終状態になったときのステップが、「第２の所定ステップ」にあたる。本実施例では、領域１００において津波が陸上の最深部まで達し引き始めた時点を最終状態とする。

そして、３次元ＳＰＨ計算処理部３は、算出したステップ毎の各粒子における速度ｖ、位置ｒ、圧力Ｐ及び密度ρのデータを結果記憶部４へ出力する。

結果記憶部４は、３次元ＳＰＨ計算処理部３から受信した解析結果であるステップ毎の各粒子における速度ｖ、位置ｒ、圧力Ｐ及び密度ρをメモリやハードディスク等に格納する。

次に、図３を参照して本実施例に係る解析処理の流れについて説明する。図３は、実施例１に係るシミュレーション装置による解析処理のフローチャートである。

２次元伝播計算処理部２は、ユーザインタフェース１を用いて入力された入力データである各格子における、水位η、ｘ方向の流量Ｍ及びｙ方向の流量Ｎの初期値を取得する（ステップＳ１）。

２次元伝播計算処理部２は、水位η、ｘ方向の流量Ｍ及びｙ方向の流量Ｎの初期値を用いて、各ステップにおける２次元伝播計算を最終ステップまで実行する（ステップＳ２）。

３次元ＳＰＨ計算処理部３は、ステップ毎の各格子における、水位η、ｘ方向の流量Ｍ及びｙ方向の流量Ｎの入力を２次元伝播計算処理部２から受ける。次に、３次元ＳＰＨ計算処理部３は、３次元粒子法計算を開始するステップを特定する。そして、３次元ＳＰＨ計算処理部３は、特定したステップにおける２次元伝播計算の水位η、ｘ方向の流量Ｍ及びｙ方向の流量Ｎのデータから３次元粒子法計算の初期条件を生成する（ステップＳ３）。

３次元ＳＰＨ計算処理部３は、現ステップの速度ｖ、位置ｒ、圧力Ｐ及び密度ρを用いて、１ステップ後の３次元粒子法計算を実行し、次のステップの速度ｖ、位置ｒ、圧力Ｐ及び密度ρを算出する（ステップＳ４）。

３次元ＳＰＨ計算処理部３は、３次元粒子法計算の終了条件を満たすか否かを判定する（ステップＳ５）。終了条件を満たしていない場合（ステップＳ５：否定）、３次元ＳＰＨ計算処理部３は、２次元伝播計算処理部２から受信した各格子における、２次元伝播計算の水位η、ｘ方向の流量Ｍ及びｙ方向の流量Ｎから次のステップの境界条件を生成する（ステップＳ６）。その後、３次元ＳＰＨ計算処理部３は、ステップＳ４に戻る。

これに対して、終了条件を満たしている場合（ステップＳ５：肯定）、３次元ＳＰＨ計算処理部３は、解析結果を結果記憶部４に格納し、解析処理を終了する。

以上に説明したように、本実施例に係るシミュレーション装置は、連続体の広い移動に対するシミュレーションにおいて、広い範囲では２次元伝播計算を行い、港湾部や都市部といった３次元の形状の影響が大きい領域では３次元粒子法計算を行う。これにより、津波のように大量の流体の移動が伴う事象においても、破綻無く計算を行うことができ、計算量を抑えながら、精度の良いシミュレーションを行うことができる。

次に、実施例２について説明する。本実施例に係るシミュレーション装置は、全体的な処理は実施例１と同様であるが、処理方法を具体的にしたものである。そこで、以下では、２次元伝播計算及び３次元粒子法計算の処理について主に説明する。また、以下で特に説明の無い限りは、図１と同じ符号を有する各部は同じ機能を有しているものとする。

図４は、実施例２に係る２次元伝播計算処理部のブロック図である。図４に示すように、２次元伝播計算処理部２は、制御部２１、水位算出部２２及び流量算出部２３を有している。

制御部２１は、領域１０１を平面視した２次元平面を分割した格子を予め記憶している。さらに、制御部２１は、領域１０１の点を表す座標をｘｙ座標としたときの、分割した各格子における、水位η、ｘ方向の流量Ｍ及びｙ方向の速度Ｙの初期値である水位η_０、流量Ｍ_０及び流量Ｎ_０の入力をユーザインタフェース１から受ける。

制御部２１は、各格子における水位η_０、流量Ｍ_０及び流量Ｎ_０を水位算出部２２及び流量算出部２３に出力する。その後、制御部２１は、ステップ１での各格子における水位ηの入力を水位算出部２２から受ける。また、制御部２１は、ステップ１での各格子における流量Ｍ及び流量Ｎの入力を流量算出部２３から受ける。以下では、初期値である水位η_０、流量Ｍ_０及び流量Ｎ_０で表される初期状態をステップ０とし、その後の各ステップをステップ１，２，３，・・・とする。そして、各ステップをまとめて、ステップｋ（ｋ＝０，１，２，・・・）と表す。

その後、制御部２１は、受信したステップｋ（ｋは自然数）での水位η、流量Ｍ及び流量Ｎを水位算出部２２及び流量算出部２３に出力し、ステップｋ＋１での水位η、流量Ｍ及び流量Ｎの入力を水位算出部２２及び流量算出部２３から受けることを繰り返す。

また、制御部２１は、最終状態の条件の入力をユーザインタフェース１から受ける。本実施例では、制御部２１が、最終状態の条件として、津波が陸上の最深部まで達し引き始めた時点という条件の入力を受ける。そして、制御部２１は、ステップｋまでの計算結果を受信して、ステップ１〜ｋまでの各格子における水位η、流量Ｍ及び流量Ｎを用いて、最終状態の条件を満たしているか否かを判定する。最終状態の条件を満たしている場合、制御部２１は、現ステップであるステップｋを最終ステップとする。その場合、制御部２１は、受信したステップｋでの各格子における水位η、流量Ｍ及び流量Ｎの水位算出部２２及び流量算出部２３への出力を停止する。

その後、制御部２１は、ステップ１〜ｋの各ステップでの各格子における水位η、流量Ｍ及び流量Ｎを３次元ＳＰＨ計算処理部３へ出力する。また、制御部２１は、各格子の位置の情報を３次元ＳＰＨ計算処理部３へ出力する。

水位算出部２２は、次に示す式（１）で表される連続方程式を記憶している。

水位算出部２２は、ステップｋ＋１では、ステップｋの各格子における水位η、流量Ｍ及び流量Ｎの入力を制御部２１から受ける。ここで、ｘ方向の流量Ｍは、Ｍ＝ｖ_ｘ×Ｄ（ｖ_ｘはｘ方向の速度の関数、Ｄは全水深であり、水位η＋静水深ｈである。）として表される。また、ｙ方向の流量Ｎは、Ｍ＝ｖ_ｙ×Ｄ（ｖ_ｙはｙ方向の速度の関数である。）として表される。

そこで、水位算出部２２は、ステップｋ＋１の計算において、ステップｋにおける流量Ｍ及び流量Ｎを式（１）に代入して、ステップｋにおける水位の時間微分である∂η／∂ｔを求める。

さらに、水位算出部２２は、次に示す式（２）で表される時間積分の式を記憶している。ここで、Δｔは、各ステップ間の時間を表している。

水位算出部２２は、ステップｋ＋１の計算においては、ステップｋにおける水位の時間微分である∂η／∂ｔを式（２）に代入し、ステップｋの各格子における水位ηをステップｎにおける水位の時間微分である∂η／∂ｔに基づいて時間積分する。これにより、水位算出部２２は、ステップｋ＋１における水位ηを算出する。

そして、水位算出部２２は、算出したステップｋにおける水位ηを制御部２１へ出力する。

流量算出部２３は、次に示す式（３）及び（４）で表される運動方程式を記憶している。なお、ｇは重力加速度、Ｄは全水深（Ｄ＝ｈ＋η、ただし、ｈは静水深）、ｎは海底摩擦を表す粗度係数である。

流量算出部２３は、ステップｋ＋１では、ステップｋの各格子における水位η、流量Ｍ及び流量Ｎの入力を制御部２１から受ける。

流量算出部２３は、ステップｋ＋１の計算においては、ステップｋにおける水位η、流量Ｍ及び流量Ｎを式（３）に代入して、ステップｋにおけるｘ方向の速度の時間微分である∂Ｍ／∂ｔを求める。

また、流量算出部２３は、ステップｋ＋１の計算においては、ステップｋにおける水位η、流量Ｍ及び流量Ｎを式（４）に代入して、ステップｋにおけるｙ方向の速度の時間微分である∂Ｎ／∂ｔを求める。

さらに、水位算出部２２は、次に示す式（５）及び（６）で表される時間積分の式を記憶している。ここで、Δｔは、各ステップ間の時間を表している。

流量算出部２３は、ステップｋ＋１の計算においては、ステップｋにおけるｘ方向の速度の時間微分である∂Ｍ／∂ｔを式（５）に代入し、ステップｋの各格子における流量Ｍをステップｋにおけるｘ方向の速度の時間微分である∂Ｍ／∂ｔに基づいて時間積分する。これにより、流量算出部２３は、ステップｋ＋１における流量Ｍを算出する。

また、流量算出部２３は、ステップｋ＋１の計算においては、ステップｋにおけるｙ方向の速度の時間微分である∂Ｎ／∂ｔを式（６）に代入し、ステップｋの各格子における流量Ｎをステップｋにおけるｙ方向の速度の時間微分である∂Ｎ／∂ｔに基づいて時間積分する。これにより、流量算出部２３は、ステップｋ＋１における流量Ｎを算出する。

そして、流量算出部２３は、算出したステップｋ＋１における流量Ｍ及び流量Ｎを制御部２１へ出力する。

次に、図５を参照して本実施例に係る３次元ＳＰＨ計算処理部について説明する。図５は、実施例２に係る３次元ＳＰＨ計算処理部のブロック図である。図５に示すように、３次元ＳＰＨ計算処理部３は、制御部３１、速度算出部３２、密度算出部３３、圧力算出部３４、位置算出部３５及び粒子生成部３６を有している。ここで、３次元粒子計算におけるステップは、２次元伝播計算におけるステップよりも細かい。そこで、以下では、３次元粒子計算におけるステップを２次元伝播計算におけるステップと区別するためステップ＃ｎと記載する。

制御部３１は、例えば、図２の領域１００の境界における水位ηが５０ｃｍを超えた時点を、３次元粒子法計算を開始する条件として記憶している。また、制御部３１は、粒子と粒子の運動の境界条件を設定するための固定境界要素を予め記憶している。固定境界要素は、海底面、地面及び建造物の表面などを微小部分に分割した面要素をモデル化したものである。固定境界要素は、例えば、微小な円盤の集合として境界全体をモデル化して、各境界要素に中心座標、面要素の法線ベクトル、面要素の面積などを含んでも良い。また、固定境界要素として、ポリゴンの集合として境界全体を表現して、各境界要素に複数の頂点の位置座標を持たせても良い。

制御部３１は、２次元伝播計算の計算結果であるステップ毎の各格子における水位η、流量Ｍ及び流量Ｎの入力を２次元伝播計算処理部２から受ける。また、制御部３１は、各格子の位置の情報の入力を２次元伝播計算処理部２から受ける。

そして、制御部３１は、各格子の位置の情報及びステップ毎の各格子における水位ηの中から、境界における水位ηが５０ｃｍを超えた時点のステップをステップ＃０として特定する。

そして、制御部３１は、ステップ＃０での各格子における水位η、流量Ｍ及び流量Ｎを抽出する。

ここで、３次元粒子法計算は連続体の動作をより詳細に解析するために、２次元伝播計算で用いられた格子よりも細かい精度で計算を行う。そこで、制御部３１は、３次元粒子法計算で用いる点を決めるための格子をより細かく分割する数量を予め記憶している。そして、制御部３１は、ステップ＃０での各格子における水位ηの内挿を行い、３次元粒子法計算で用いる位置座標における水位ηを求める。次に、制御部３１は、３次元粒子法計算で用いる各位置座標における水位ηに静水深ｈを加算して、各位置座標における全水深Ｄを求める。また、制御部３１は、３次元粒子法計算で用いる位置座標におけるステップ＃０での流量Ｍ及び流量Ｎを内挿により求める。そして、制御部３１は、各位置座標でのステップ＃０での流量Ｍ及び流量Ｎを全水深で除算し、各位置座標での速度の初期値ｖ_０を求める。

そして、制御部３１は、各位置座標における全水深Ｄを用いて、各位置座標に予め決められたサイズの粒子を配置する。例えば、粒子の直径が１ｍと決められており、ある位置座標における全水深が５０ｍの場合、制御部３１は、その位置座標に５０個の粒子を積み重ねて配置する。これにより、各位置座標における粒子の位置の初期値が決定する。

さらに、制御部３１は、各位置座標に配置した各粒子の速度の初期値として、算出したその位置座標における速度ｖ_０を割り当てる。図６は、初期値の設定を模式的に表した図である。図６は、図２に示す領域１００の断面の一部を模式的に現している。図６の縦軸は水深を表し、横軸は静水面と平行な特定の方向を表している。図６に示すように、制御部３１は、各位置座標において粒子２００を水深方向に積み上げて配置する。そして、制御部３１は、矢印２０１で示されるように、矢印２０１の大きさに応じた速度を各粒子に割り当てる。図６では、例えば、横軸が増える方向を沖に向かう方向とすると、横軸が減る方向に向かう矢印２０１は、粒子が陸に向かっていることを表している。そして、矢印２０１が大きいほど粒子の速度が速いことを表している。

各粒子の位置及び速度の初期値の設定が終わると、制御部３１は、各粒子の影響範囲内にある他の粒子のリストを作成する。以下では、各粒子の影響範囲内にある他の粒子を「近傍粒子」と呼ぶ。また、近傍粒子のリストを「近傍粒子リスト」と呼ぶ。

図７は、影響範囲及び近傍粒子について説明するための図である。ここで、図７を参照して、各粒子の影響範囲及び近傍粒子について説明する。ここでは、着目粒子と他の粒子との距離が２ｈ_ＳＰＨ以内の場合に、相互に影響を与え合うものとして説明する。ここで、影響とは、例えば、粒子が移動する場合に、他の粒子へ力を加えるなどのことを言う。

ここでは、図７の粒子３０１を例に説明する。粒子３０１の影響範囲３０２は、粒子３０１を中心として半径２ｈ_ＳＰＨの円の中に含まれる範囲である。この半径２ｈ_ＳＰＨは、粒子３０１の影響半径である。そして、影響範囲３０２の中に含まれる他の粒子が、粒子３０１の近傍粒子である。図７では、例えば、粒子３０３のように、斜線で表されている粒子が近傍粒子である。すなわち、粒子３０１が影響を与えることになる他の粒子は、斜線で表される近傍粒子ということになる。逆に、斜線で表される近傍粒子からのみ、粒子３０１は影響を受けることになる。ここで、近傍粒子が数十個程度含まれるように半径２ｈ_ＳＰＨを設定すれば、着目した粒子と他の粒子間の相互作用を解析するには十分であることが経験的に知られている。

また、粒子３０１の位置を表す位置ベクトルをｒ_ａとし、粒子３０３の位置を表す位置ベクトルをｒ_ｂとすると、粒子３０１と粒子３０３との間の距離は｜ｒ_ａ−ｒ_ｂ｜となる。すなわち、｜ｒ_ａ−ｒ_ｂ｜≦２ｈを満たす場合に、粒子３０３は粒子３０１の近傍粒子とされる。

制御部３１は、例えば図７の粒子３０１に対する近傍リストを作る場合、まず、影響範囲３０２に含まれる粒子３０３のような斜線で表される各粒子を粒子３０１の近傍粒子として抽出する。例えば、制御部３１は、着目粒子である粒子３０１の位置と他の粒子の位置から粒子３０１と各他の粒子との距離を求め、その距離が２ｈ_ＳＰＨ以下となる粒子を近傍粒子として抽出する。

図５に戻って、制御部３１は、近傍粒子の抽出を全ての粒子に対して行う。そして、制御部３１は、各着目粒子と抽出した各着目粒子の近傍粒子とが対応する近傍粒子リストを作成する。そして、制御部３１は、作成した近傍粒子リストを速度算出部３２、密度算出部３３、圧力算出部３４、位置算出部３５及び粒子生成部３６へ出力する。

その後、制御部３１は、ステップ＃ｎにおける各着目粒子の速度ｖの入力を速度算出部３２から受ける。また、制御部３１は、ステップ＃ｎにおける各着目粒子の密度ρの入力を密度算出部３３から受ける。また、制御部３１は、ステップ＃ｎにおける各着目粒子の圧力Ｐの入力を圧力算出部３４から受ける。さらに、制御部２は、ステップ＃ｎにおける各着目粒子の位置ｒの入力を位置算出部３５から受ける。ここで、１ステップとは、Δｔ_ＳＰＨ時間経過した後の状態である。例えば、ステップ＃ｎ＋１における速度とは、ステップ＃ｎにおける速度からΔｔ_ＳＰＨ時間後の速度である。

また、制御部３１は、図２の領域１００における境界から一定の距離Δｘの範囲を境界領域として記憶している。ここで、Δｘは粒子の大きさから設定されることが好ましい。具体的には、粒子の大きさ程度の幅を持たせることが好ましい。例えば、粒子の直径が１ｍの場合、１ｍより少し大きい幅をΔｘとすることが好ましい。そして、制御部３１は、受信したステップ＃ｎにおける着目粒子の位置ｒから、図２の領域１００における境界の外に着目粒子が出たか否かを判定する。そして、制御部３１は、着目粒子が境界の外に出たと判定した場合、ステップ＃ｎ＋１以降の計算においてその着目粒子を除外する。図８は、粒子の流出を模式的に表した図である。図８に示すように、粒子２０２が境界の外に出て粒子２０３となった場合、制御部３１は、粒子２０３を計算から除外する。

さらに、制御部３１は、ステップ＃ｎにおいて生成された粒子の位置ｒの入力を粒子生成部３６から受ける。そして、制御部３１は、ステップ＃ｎ＋１以降の計算において受信した新たに生成された粒子を組み入れる。

また、制御部３１は、ステップ＃ｎにおける境界領域内に存在する粒子の速度ｖを、２次元伝播計算処理部２から受信した流量Ｍ及び流量Ｎから求め、速度算出部３２及び位置算出部３５から受信したステップ＃ｎにおける速度ｖと置き換える。また、制御部３１は、ステップ＃ｎにおける境界領域内に存在する粒子の密度ρ及び圧力Ｐは、海水の一般的な密度及び圧力に置き換える。このようにして、制御部３１は、ステップ＃ｎにおける境界条件を都度求め、境界領域内に存在する粒子に割り当てる。これは、本実施例での３次元粒子法計算においては境界外には粒子が存在しないため、境界領域に存在する粒子の動きは正確に求めることができない。そこで、２次元伝播計算によって求められた正確性の高い情報を基に境界領域に存在する粒子の動きを推定することで、計算の正確性を向上させるためである。

ただし、上述したように２次元伝播計算と３次元粒子法計算ではステップの精度が異なる。すなわち、３次元粒子法計算の各ステップ＃ｎに一致する、２次元伝播計算のステップｋがあるとは限らない。そこで、制御部３１は、３次元粒子法計算の各ステップ＃ｎに一致するタイミングの水位η、流量Ｍ及び流量Ｎを算出する。ここで、図９は、境界条件算出における補間を説明するための図である。図９では、例えば、２次元伝播計算のステップｋにおいて３次元粒子法計算が開始され、２次元伝播計算の次のステップであるステップｋ＋１と３次元粒子法計算のステップ＃ｎのタイミングが一致する場合を表している。この場合、２次元伝播計算のステップ間の時間はΔｔであり、３次元粒子法計算のステップ間の時間はΔｔ_ＳＰＨである。この場合、ステップ＃１では、対応４０１のように、制御部３１は、ステップｋにおける水位η、流量Ｍ及び流量Ｎを用いて、境界条件を算出することができる。また、ステップ＃ｎでは、制御部３１は、対応４０２のように、ステップｋ＋１における水位η、流量Ｍ及び流量Ｎを用いて、境界条件を算出することができる。しかし、ステップ＃２〜ステップ＃ｎ−１に対応する、２次元伝播計算のステップは存在しない。そこで、例えばステップ＃２であれば、制御部３１は、ステップｋ及びｋ＋１における水位η、流量Ｍ及び水位Ｎを用いて、内挿を行いステップ＃２に対応する内挿ステップ４０３における水位η、流量Ｍ及び水位Ｎを求める。そして、制御部３１は、対応４０４のように、求めた内挿ステップ４０３における水位η、流量Ｍ及び流量Ｎを用いて、ステップ＃２の境界条件を算出する。

例えば、内挿ステップにおける水位は、次の式（７）のように求めることができる。ここでは、粒子法計算において境界条件を得たいステップの時刻をｔ_ＳＰＨとしている。そして、ｔ_ＳＰＨに最も近い２次元伝播計算のステップであるステップｎ及びステップｎ＋１の時刻をｎΔｔ及び（ｎ＋１）Δｔとしている。

そして、制御部３１は、ステップ＃ｎにおける各着目粒子の密度ρ、圧力Ｐ及び位置ｒを、ステップ＃ｎ＋１の計算に用いる各着目粒子の現ステップの密度、圧力及び位置として、速度算出部３２へ出力する。

また、制御部３１は、ステップ＃ｎにおける各着目粒子の密度ρをステップ＃ｎ＋１の計算に用いる現ステップの密度として密度算出部３３へ出力する。

また、密度算出部３３は、ステップ＃ｎ＋１における各着目粒子の密度ρをステップ＃ｎ＋１の圧力計算に用いるために圧力算出部３４へ出力する。

また、制御部３１は、ステップ＃ｎにおける各着目粒子の速度ｖを、ステップ＃ｎ＋１の計算に用いる現ステップの速度として位置算出部３５へ出力する。

以上の動作をまとめると、制御部３１は、ステップ＃ｎからΔｔ_ＳＰＨ時間経過した後のステップ＃ｎ＋１における、各着目粒子の速度ｖ、密度ρ、圧力Ｐ及び位置ｒの入力を受ける。次に、制御部３１は、粒子の位置から消滅させる粒子を決定し、さらに粒子生成部３６からの情報を用いて粒子を生成する。そして、制御部３１は、粒子の消滅及び生成を行った後、２次元伝播計算により算出された水位η、流量Ｍ及び流量Ｎから境界領域にある粒子の速度ｖ、密度ρ、圧力Ｐ及び位置ｒを求めて境界条件を設定する。その後、制御部３１は、それらの情報を各部に出力し、ステップ＃ｎからΔｔ_ＳＰＨ時間経過した後のステップ＃ｎ＋１における、各着目粒子の速度ｖ、密度ρ、圧力Ｐ及び位置ｒの入力を受ける。このように、制御部３１は、既にステップ＃ｎにおける速度ｖ、密度ρ、圧力Ｐ及び位置ｒを受信している場合には、次にステップ＃ｎ＋１における各着目粒子の速度ｖ、密度ρ、圧力Ｐ及び位置ｒを受信する。そして、制御部３１は、指定されたシミュレーションの終了のタイミングまで、各ステップにおける各着目粒子の速度ｖ、密度ρ、圧力Ｐ及び位置ｒを順次受信する。

さらに、制御部３１は、ステップ＃ｎにおける各着目粒子の最新の速度ｖ、密度ρ、圧力Ｐ及び位置ｒから１ステップ後の情報の入力を受けると、受信したステップ＃ｎにおける各情報を結果記憶部４へ記憶させる。

そして、制御部３１は、受信したステップ＃ｎにおける各着目粒子の位置を用いてステップ＃ｎにおける各着目粒子の近傍粒子を求め、近傍粒子リストを作成する。そして、制御部３１は、作成した近傍粒子リストを速度算出部３２、密度算出部３３、圧力算出部３４、位置算出部３５及び粒子生成部３６へ出力する。このように、制御部３１は、ステップ毎の速度ｖ、密度ρ、圧力Ｐ及び位置ｒ、並びに近傍粒子リストを各部へ出力する。そして、制御部３１は、現在のデータを用いて１ステップ後のデータを算出するよう、速度算出部３２、密度算出部３３、圧力算出部３４、位置算出部３５及び粒子生成部３６に対して指示し、最終ステップまで繰り返させる。

速度算出部３２は、ステップ＃ｎ＋１の計算に用いる各着目粒子の現ステップの密度、圧力及び位置として、ステップ＃ｎにおける各着目粒子の密度ρ、圧力Ｐ及び位置ｒの入力を制御部３１から受ける。

速度算出部３２は、次の式（８）で表される運動方程式を記憶している。

ここで、式（８）における、変数及び係数の添え字は各粒子を表す。例えば、ｒ_ａ，ｖ_ａ，ρ_ａ，Ｐ_ａ，ｍ_ａがあった場合、それらは、粒子ａの位置ベクトル、速度ベクトル、密度、圧力、質量を表している。Ｗは、粒子の分布から連続場を構成するために用いるカーネル関数である。カーネル関数としては、例えば、次の式（９）に示すような５次のスプライン関数などが用いられる。ただし、ここでｑ＝ｒ／ｈ_ＳＰＨである。

また、式（８）におけるΠ_ａｂは粒子ａ，ｂ間に働く人工粘性項である。例えば、次の式（１０）に示されるような形式が採用される。

ここで、式（１０）におけるμは粘性係数であり、次の式（１１）で表される。ただし、ここでη_ＳＰＨは分母が０になることを避けるためのパラメータであり、考えている空間スケールに対して十分小さい値を設定する。例えば、η_ＳＰＨ＝０．０１ｈ_ＳＰＨ等でよい。

そして、速度算出部３２は、式（８）に受信したステップ＃ｎにおける着目粒子ａ及び近傍粒子ｂの密度ρ、圧力Ｐ及び位置ｒを代入し、ステップ＃ｎにおける加速度ｄｖ／ｄｔを求める。

さらに、速度算出部３２は、境界からの反発力を次の式（１２）を用いて求める。

ここで、図１０を参照して、式（１２）について説明する。図１０は、境界からの反発力を説明するための図である。ここでは、図１０における粒子５００に加えられる境界からの反発力について説明する。

図１０に示すように、ｄは境界面５０１からの所定の距離である。そして、境界層５０２、境界面からｄの幅の層を示している。また、ｓは粒子５００の境界層５０２への侵入距離である。

本実施例では、境界面５０１から幅ｄの境界層５０２を考え、境界層５０２内に粒子５００が突入した場合に、粒子５００に対して境界面５０１から粒子５００の向きに反発加速度が加わるものとする。反発力ｆ_ｎは、弾性ばねによるものと同じ形式でｆ_ｎ＝ｋｓとして与える。ばね剛性係数ｋは、音速ｃ_０の粒子が有する運動エネルギーｍｃ_０ ^２／２と、粒子が境界層５０２へｄの長さだけ進入したときの弾性歪エネルギーｋｄ^２／２とが同程度となるように定める。すなわち、連続体の運動の解析を行なうシミュレーションにおいては、粒子５００は音速を超えることはないと考えられるので、音速で粒子５００が進入した場合に境界面５０１においてエネルギーがなくなるように設定しておけば、粒子５００は境界面を突き抜けることはなくなる。このため、βを１程度のパラメータとする。このようにβを設定することで、式（１２）が導かれる。

式（１２）で表される反発力を用いることで、境界壁からの粒子のすり抜けを防止し、導入する摩擦やfree slip条件の正確性を向上でき、さらに、境界で物理的に無意味なエネルギー消費の発生を抑えることができる。また、式（１２）は簡単な式で表されているので、式（１２）で表される反発力を用いた場合には、計算負荷を抑えることができる。

ここで、本実施例では式（１２）を用いたが、この反発力は他の式を用いてもよい。例えば、反発力が過大になる可能性を考慮して、粒子が境界面に近づいている場合には、弾性ばねを考慮し、逆に粒子が境界面から遠ざかる場合には、弾性ばねを考慮しないように反発力を設定してもよい。このような反発力は一方向弾性ばね法と呼ばれることがある。エネルギー保存の観点からは粒子と境界面の相対速度に関係なく弾性ばねを考慮する式（１２）のような反発力が好ましいが、境界での不自然な反発を避けるために一方向弾性ばね法をもちいてもよい。

そして、速度算出部３２は、運動方程式である式（８）を用いて算出した加速度ｄｖ／ｄｔと式（１２）を用いて求めた反発力とを加えてステップ＃ｎにおける加速度ｖ´を算出する。ここでは、加速度ｖ´は、加速度ベクトルを表す。

さらに、速度算出部３２は、次の時間積分の式（１３）を記憶している。ここで、ｖ_ｎはステップ＃ｎにおける速度を表し、ｖ_ｎ´はステップ＃ｎにおける加速度を表している。

そして、速度算出部３２は、求めたステップ＃ｎの速度ｖ及び加速度ｖ´を式（１３）に代入して、ステップ＃ｎの速度ｖをステップ＃ｎにおける加速度ｖ´で時間積分する。これにより、速度算出部３２は、ステップ＃ｎ＋１の速度ｖを求める。

速度算出部３２は、算出したステップ＃ｎ＋１の各粒子の速度ｖを制御部３１へ出力する。また、速度算出部３２は、算出したステップ＃ｎ＋１の各粒子の速度ｖを密度算出部３３へ出力する。さらに、速度算出部３２は、算出したステップ＃ｎ＋１の各粒子の速度ｖを位置算出部３５へ出力する。

密度算出部３３は、ステップ＃ｎ＋１の計算に用いる現ステップの密度として、ステップ＃ｎにおける各着目粒子の密度ρの入力を制御部３１から受ける。さらに、密度算出部３３は、ステップ＃ｎ＋１の各粒子の速度ｖの入力を速度算出部３２から受ける。

密度算出部３３は、次の式（１４）で表される連続方程式を記憶している。

密度算出部３３は、速度算出部３２から入力されたステップ＃ｎ＋１における各粒子の速度ｖを式（１４）に代入し、ステップ＃ｎ＋１における各粒子の密度時間微分ｄρ／ｄｔを算出する。

さらに、密度算出部３３は、次の式（１５）で表される時間積分の式を記憶している。

密度算出部３３は、算出したステップ＃ｎにおける各粒子の密度及び各粒子の密度時間微分ｄρ_ａ／ｄｔの値を式（１５）に代入して、ステップ＃ｎの密度ρをステップ＃ｎにおける密度時間微分ｄρ／ｄｔで時間積分する。これにより、密度算出部３３は、ステップ＃ｎ＋１の密度ρを求める。

密度算出部３３は、算出したステップ＃ｎ＋１の密度ρを、制御部３１及び圧力算出部３４へ出力する。

圧力算出部３４は、ステップ＃ｎ＋１の密度ρの入力を密度算出部３３から受ける。

圧力算出部３４は、次の式（１６）で表される状態方程式を記憶している。ここで、ｃ_ｓ及びρ_０は、それぞれ流体の音速と基準密度である。

圧力算出部３４は、密度算出部３３から入力されたステップ＃ｎ＋１の密度ρを式（１６）に代入して、ステップ＃ｎ＋１における圧力Ｐを算出する。

圧力算出部３４は、算出したステップ＃ｎ＋１における圧力Ｐを制御部３１へ出力する。

位置算出部３５は、ステップ＃ｎ＋１における速度ｖの入力を速度算出部３２から受ける。

位置算出部３５は、次の式（１７）で表される時間積分の式を記憶している。

位置算出部３５は、速度算出部３２から入力されたステップ＃ｎ＋１における速度ｖを式（１６）に代入して、ステップ＃ｎの位置ｒをステップ＃ｎ＋１における速度ｖで時間積分する。これにより、位置算出部３５は、ステップ＃ｎ＋１の位置ｒを求める。

粒子生成部３６は、境界面を１つの粒子に相当する所定の領域で分割した分割領域を記憶している。ここで、所定の領域は、密に配置されていればどのような形状でも良く、例えば、矩形でも円形でもよい。以下では、この所定の領域を「粒子生成子」と呼ぶ。粒子生成子は、境界において流入があった時に粒子を生成することで、粒子法計算を破たんなく実行させるために設定されるものである。粒子生成子は、領域１００の境界面を微小部分に分割した面要素として、例えば、中心座標、面要素の法線ベクトル、面要素の面積、流入出体積、基準体積を持つ。ここで、基準体積とは、粒子一つに対応する体積である。そして、以下で説明する粒子生成子における流出入体積は、計算開始の時点では０に設定され、流入があれば増加し、流出があれば減少する。粒子生成子を通じて海水が流出している場合には流出入体積は負の値を取る。

粒子生成部３６は、ステップ＃ｎにおける各粒子生成子の位置における２次元伝播計算により求められた速度の入力を制御部３１から受ける。

そして、粒子生成部３６は、受信した速度から各粒子生成子の位置における海水の流速を表す流速ベクトルを求める。さらに、粒子生成部３６は、粒子生成子が形成する所定の領域の法線方向を向いた大きさが面積に比例する面積ベクトルを生成する。粒子生成子の面積をＳとし、法線ベクトルをｕとした場合、面積ベクトルは、例えば、ｕＳとして表すことができる。

そして、粒子生成部３６は、各粒子生成子における流速ベクトルと面積ベクトルとの内積を求め、ステップ＃ｎにおける各粒子生成子における海水の流出入体積の増減分を算出する。

そして、粒子生成部３６は、各粒子生成子におけるステップ＃ｎの海水の流出入体積に流出入体積の増減分を加算して、各粒子生成子における新しいステップ＃ｎの海水の流出入体積を求める。例えば、古い流出入体積をＶ^ｏｌｄとして、新しい流出入体積をＶ^ｎｅｗとした場合、粒子生成部３６は、次の式（１８）により新しい海水の流出入体積が算出できる。

次に、粒子生成部３６は、各粒子生成子におけるステップ＃ｎ＋１の海水の流出入体積が１つの粒子の体積と等しい基準体積以上か否かを判定する。海水の流出入体積が基準体積以上となっている粒子生成子があれば、粒子生成部３６は、その粒子生成子に対応する決められた位置への粒子の生成を決定する。ここで、粒子生成子に対応する決められた位置とは、例えば、粒子生成子に接するように粒子を配置する位置でも良いし、粒子生成子から法線方向に所定の距離離れた位置に粒子が配置される位置でも良い。そして、粒子生成部３６は、粒子の生成を決定した粒子生成子におけるステップ＃ｎ＋１の海水の流出入体積から基準体積を減算して、ステップ＃ｎ＋１の海水の流出入体積とする。

図１１は、粒子の流入を模式的に表した図である。境界の外部にある粒子２０４の体積を有する海水が境界を越えて内部に流入した場合に、粒子生成部３６は、流入した粒子２０４に対応する境界の内部の位置に新しく粒子２０５を生成することを決定する。

また、ステップ＃ｎ＋１の海水の流出入体積がマイナス方向の基準体積を下回った粒子生成子がある場合、粒子生成部３６は、その粒子生成子におけるステップ＃ｎ＋１の海水の流出入体積に基準体積を加算して、ステップ＃ｎ＋１の海水の流出入体積とする。

そして、粒子生成部３６は、生成を決定した粒子の位置を制御部３１に通知する。生成した粒子の位置とは、例えば、生成した粒子の中心の位置などでよい。

粒子生成部３６は、ステップが終了するまで、上記の粒子生成の処理を繰り返す。

次に、図１２を参照して、本実施例に係る２次元伝播計算の処理について説明する。図１２は、実施例２に係る２次元伝播計算の処理のフローチャートである。

制御部２１は、ユーザインタフェース１からの入力により、領域１０１の分割した各格子における、水位、ｘ方向の速度及びｙ方向の速度の初期値である水位η_０、流量Ｍ_０及び流量Ｎ_０を取得する（ステップＳ１０１）。また、制御部２１は、各格子の中から解析対象とする格子を選択格子として一つ選択する。そして、制御部２１は、選択格子における、水位η_０、流量Ｍ_０及び流量Ｎ_０を水位算出部２２及び流量算出部２３へ出力する。

水位算出部２２は、選択格子における水位の時間微分∂η／∂ｔを式（１）の連続方程式を用いて算出する（ステップＳ１０２）。

また、流量算出部２３は、選択格子における流量Ｍ及び流量Ｎの時間微分である∂Ｍ／∂ｔ及び∂Ｎ／∂ｔを式（３）及び式（４）の運動方程式を用いて算出する（ステップＳ１０３）。

制御部２１は、全ての格子について計算が完了したか否かを判定する（ステップＳ１０４）。全ての格子について計算が完了していない場合（ステップＳ１０４：否定）、制御部２１は、新たな選択格子を選択し、ステップＳ１０２へ戻る。

これに対して、全ての格子について計算が完了している場合（ステップＳ１０４：肯定）、制御部２１は、選択格子を一つ選択し、水位算出部２２及び流量算出部２３に次のステップの水位及び速度の算出を指示する。

水位算出部２２は、選択格子の水位ηを選択格子の水位の時間微分∂η／∂ｔを基に時間積分し（ステップＳ１０５）、次のステップの水位を求める。そして、水位算出部２２は、算出した次のステップの水位ηを制御部２１へ出力する。

また、流量算出部２３は、選択格子の流量Ｍ及び流量Ｎを選択格子の速度の時間微分∂Ｍ／∂ｔ及び∂Ｎ／∂ｔを基に時間積分し（ステップＳ１０６）、次のステップの流量Ｍ及び流量Ｎを求める。そして、流量算出部２３は、算出した次のステップの流量Ｍ及び流量Ｎを制御部２１へ出力する。

制御部２１は、全ての格子について計算が完了したか否かを判定する（ステップＳ１０７）。全ての格子について計算が完了していない場合（ステップＳ１０７：否定）、制御部２１は、新たな選択格子を選択し、ステップＳ１０５へ戻る。

これに対して、全ての格子について計算が完了している場合（ステップＳ１０７：肯定）、制御部２１は、解析段階を１ステップ進める（ステップＳ１０８）。

そして、制御部２１は、現在の解析段階における計算結果を３次元ＳＰＨ計算処理部３へ出力する（ステップＳ１０９）。ここで、このフローでは、制御部２１は、ステップ毎に計算が終わった時点で計算結果を出力しているが、これに限らず、例えば、全てのステップが完了してから全てのステップの計算結果をまとめて出力しても良い。

制御部２１は、最終状態の条件を満たしたか否かにより、全ステップが完了したか否かを判定する（ステップＳ１１０）。全ステップが完了していない場合（ステップＳ１１０：否定）、制御部２１は、各格子の中から解析対象とする格子を選択格子として一つ選択する。そして、制御部２１は、選択格子における、現ステップの水位η、流量Ｍ及び流量Ｎを水位算出部２２及び流量算出部２３へ出力し、次のステップの計算を指示し、ステップＳ１０２へ戻る。

これに対して、全ステップが完了している場合（ステップＳ１１０：肯定）、制御部２１は２次元伝播計算の処理を終了する。

次に、図１３を参照して、本実施例に係る３次元ＳＰＨ計算の処理について説明する。図１３は、実施例２に係る３次元ＳＰＨ計算の処理のフローチャートである。

制御部３１は、ステップ毎の各格子における水位η、ｘ方向の流量Ｍ及びｙ方向の流量Ｎのデータを２次元伝播計算処理部２から取得する（ステップＳ２０１）。そして、制御部３１は、３次元粒子法計算の開始条件を満たしたステップを特定し、特定したステップの水位η、ｘ方向の流量Ｍ及びｙ方向の流量Ｎを取得する。さらに、制御部３１は、内挿により、３次元粒子法計算に用いる位置座標を内挿により求める。そして、制御部３１は、求めた位置座標の全水深に応じて粒子を配置し、各粒子の速度、密度、圧力及び位置の初期値を設定する。

また、制御部３１は、粒子分布に基づいて近傍粒子リストを作成する（ステップＳ２０２）。そして、制御部３１は、一つ着目粒子を選択し、着目粒子の密度ρ、圧力Ｐ及び位置ｒを速度算出部３２へ出力する。また、制御部３１は、着目粒子の密度ρを密度算出部３３へ出力する。また、制御部３１は、着目粒子の速度ｖを位置算出部３５へ出力する。

速度算出部３２は、着目粒子の密度ρ、圧力Ｐ及び位置ｒの入力を制御部３１から受ける。そして、速度算出部３２は、式（８）の運動方程式に着目粒子の密度ρ、圧力Ｐ及び位置ｒを代入し着目粒子の加速度ｄｖ／ｄｔを算出する（ステップＳ２０３）。

さらに、速度算出部３２は、着目粒子に対する境界からの反発力を式（１２）を用いて算出する。着目粒子が境界近傍に存在しない場合は、反発力は０になる。そして、速度算出部３２は、算出した反発力を加速度ｄｖ／ｄｔに加算して（ステップＳ２０４）、加速度ｖ´を算出する。

速度算出部３２は、式（１３）を用いて、現ステップにおける着目粒子の速度ｖを現ステップにおける加速度ｖ´で時間積分して、次のステップの速度ｖを算出する（ステップＳ２０５）。そして、速度算出部３２は、算出した次のステップの着目粒子の速度ｖを制御部３１へ出力する。

制御部３１は、着目粒子が境界領域に存在する場合、２次元伝播計算結果から着目粒子の速度ｖを算出し、速度算出部３２から入力された速度ｖに変えて、算出した速度ｖを着目粒子の速度として設定する（ステップＳ２０６）。

制御部３１は、領域１００内に存在する全ての粒子について速度の計算が完了したか否かを判定する（ステップＳ２０７）。全ての粒子について計算が完了していない場合（ステップＳ２０７：否定）、制御部３１は、計算を行っていない粒子から着目粒子を選択し、ステップＳ２０３へ戻る。

全ての粒子について計算が完了している場合（ステップＳ２０７：肯定）、制御部３１は、着目粒子を１つ選択して、密度算出部３３、圧力算出部３４及び位置算出部３５に対して計算実行を指示し、ステップＳ２０８へ進む。

密度算出部３３は、速度算出部３２から受信した次のステップの着目粒子の速度ｖを着目粒子の密度時間微分を式（１４）の連続方程式に代入して密度時間微分ｄρ／ｄｔを算出する（ステップＳ２０８）。

次に、密度算出部３３は、制御部３１から受信した現ステップの着目粒子の密度ρ及び算出した密度時間微分ｄρ／ｄｔの値を式（１５）に代入する。そして、密度算出部３３は、現ステップにおける着目粒子の密度ρを密度時間微分ｄρ／ｄｔで積分して、次のステップの密度ρを算出する（ステップＳ２０９）。密度算出部３３は、算出した次のステップの密度ρを制御部３１へ出力する。

圧力算出部３４は、密度算出部３３から受信した次のステップの密度ρを、式（１６）の状態方程式に代入して、次のステップの圧力Ｐを算出する（ステップＳ２１０）。圧力算出部３４は、算出した次のステップの圧力Ｐを制御部３１へ出力する。

位置算出部３５は、式（１７）を用いて、現ステップの位置ｒを速度算出部３２から受信した次のステップにおける速度ｖで時間積分し次のステップの位置ｒを求める（ステップＳ２１１）。位置算出部３５は、算出した次のステップの位置ｒを制御部３１へ出力する。

制御部３１は、次のステップの位置ｒを用いて領域境界を越えて流出した粒子を特定し、特定した粒子を領域１００の中から消去する（ステップＳ２１２）。

制御部３１は、領域１００内に存在する全ての粒子について密度、圧力及び位置の計算が完了したか否かを判定する（ステップＳ２１３）。全ての粒子について計算が完了していない場合（ステップＳ２１３：否定）、制御部３１は、計算を行っていない粒子から着目粒子を選択し、ステップＳ２０８へ戻る。

これに対して、全ての粒子について計算が完了している場合（ステップＳ２１３：肯定）、制御部３１は、粒子生成を粒子生成部３６に指示する。粒子生成部３６は、各粒子生成子の位置における２次元伝播計算により求められた速度から、現在の海水の流出入体積を求め粒子の生成処理を行う（ステップＳ２１４）。この粒子の生成処理の流れについては次に詳細に説明する。そして、粒子生成部３６は、生成を決定した粒子の位置を制御部３１へ通知する。

制御部３１は、解析段階を１ステップ進める（ステップＳ２１５）。

そして、制御部２１は、現在の解析段階における計算結果を結果記憶部４へ出力する（ステップＳ２１６）。ここで、このフローチャートでは、制御部３１は、ステップ毎に計算が終わった時点で計算結果を出力しているが、これに限らず、例えば、全てのステップが完了してから全てのステップの計算結果をまとめて出力しても良い。

制御部３１は、最終状態の条件を満たしたか否かにより、全ステップが完了したか否かを判定する（ステップＳ２１７）。全ステップが完了していない場合（ステップＳ２１７：否定）、制御部３１は、ステップＳ２０２へ戻る。

これに対して、全ステップが完了している場合（ステップＳ２１７：肯定）、制御部３１は３次元ＳＰＨ計算の処理を終了する。

次に、図１４を参照して、本実施例に係る３次元ＳＰＨ計算の処理について説明する。図１４は、実施例２に係る３次元ＳＰＨ計算の処理のフローチャートである。図１４のフローチャートは、図１３のフローチャートにおけるステップＳ２１５にあたる。

粒子生成部３６は、粒子生成子を一つ選択する。そして、粒子生成部３６は、選択した粒子生成子における現ステップの海水の流速を計算する（ステップＳ３０１）。

次に、粒子生成部３６は、海水の流速と粒子生成子の面積から、選択した粒子生成子における海水の流出入体積の増減を算出する（ステップＳ３０２）。

次に、粒子生成部３６は、現ステップの流出入体積Ｖに増減分を加算し、現ステップの新しい流出入体積Ｖを算出する（ステップＳ３０３）。

次に、粒子生成部３６は、新しい流出入体積Ｖが基準体積Ｖ_０以上か否かを判定する（ステップＳ３０４）。基準体積以上の場合（ステップＳ３０４：肯定）、粒子生成部３６は、選択している粒子生成子に対応する位置に粒子の生成を決定する（ステップＳ３０５）。

次に、粒子生成部３６は、現ステップの新しい流出入体積Ｖから基準体積Ｖ_０を減算し、次のステップの流出入体積Ｖを算出する（ステップＳ３０６）。

これに対して、新しい流出入体積Ｖが基準体積Ｖ_０以下の場合（ステップＳ３０４：否定）、粒子生成部３６は、新しい流出入体積Ｖが基準体積Ｖ_０のマイナス値、すなわち−Ｖ_０より小さいか否かを判定する（ステップＳ３０７）。−Ｖ_０より小さい場合（ステップＳ３０７：肯定）、粒子生成部３６は、現ステップの新しい流出入体積Ｖに基準体積Ｖ_０を加算し、次のステップの流出入体積Ｖを算出する（ステップＳ３０８）。これに対して、新しい流出入体積Ｖが−Ｖ_０以上の場合（ステップＳ３０７：否定）、粒子生成部３６は、ステップＳ３０９へ進む。

そして、粒子生成部３６は、全粒子生成子について粒子の生成処理が終了したか否かを判定する（ステップＳ３０９）。粒子の生成処理が終了していない粒子生成子がある場合（ステップＳ３０９：否定）、粒子生成部３６は、ステップＳ３０１へ戻る。

これに対して、全粒子生成子について粒子の生成処理が終了した場合（ステップＳ３０９：肯定）、粒子生成部３６は、粒子の生成処理を終了する。

以上に説明したように、本実施例に係るシミュレーション装置は、連続体の広い移動に対するシミュレーションにおいて、広い範囲では２次元伝播計算を行い、港湾部や都市部といった３次元の形状の影響が大きい領域では３次元粒子方計算を行う。そして、本実施例に係るシミュレーション装置は、上述した各式などを用いて、２次元伝播計算及び３次元粒子法計算を行う。これにより、津波のように大量の流体の移動が伴う事象においても、破綻無く計算を行うことができ、計算量を抑えながら、精度の良いシミュレーションを行うことができる。

また、本実施例に係るシミュレーション装置は、バッファ領域を用いた従来技術に比べて、境界領域の幅を粒子の影響半径と同程度に抑えることができ、計算量を確実に抑えることができる。

さらに、粒子生成子を３次元空間内で向きを持つものとして設定できるので、２次元伝播計算と３次元粒子法計算の境界の形状を任意に設定でき、シミュレーションの条件の自由度を向上させることができる。

次に、実施例３について説明する。本実施例に係るシミュレーション装置は、２次元伝播計算を最後まで終わらせずに途中の状態で、並行して３次元粒子法計算を行うことが実施例１及び実施例２と異なるものである。以下では、同じ機能を有する各部については説明を省略する。

２次元伝播計算処理部２は、ステップ毎にそのステップにおける計算結果を３次元ＳＰＨ計算処理部３へ出力する。

３次元ＳＰＨ計算処理部３は、計算結果を受信するとそれまでに受信した計算結果から３次元粒子法計算を開始する条件が満たされたか否かを判定する。

３次元粒子法計算を開始する条件が満たされている場合、３次元ＳＰＨ計算処理部３は、受信した計算結果の水位及び流量から各粒子の位置及び速度の初期条件を生成する。さらに、３次元ＳＰＨ計算処理部３は、各粒子の密度及び圧力の初期条件は海水の一般的な密度圧力を用いる。

そして、３次元ＳＰＨ計算処理部３は、ステップ＃ｎの速度、密度、圧力及び位置を用いて３次元粒子法計算を行い、ステップ＃ｎ＋１の速度、密度、圧力及び位置を算出する。

３次元ＳＰＨ計算処理部３は、ステップ＃ｎ＋１の速度、密度、圧力及び位置の算出後、ステップ＃ｎ＋２における境界条件、すなわち、境界領域に存在する粒子の速度が求められるか否かを判定する。

境界条件が求められるのであれば、３次元ＳＰＨ計算処理部３は、ステップ＃ｎ＋１の速度、密度、圧力及び位置及び境界条件を用いて、ステップ＃ｎ＋２の速度、密度、圧力及び位置の算出を行う。

これに対して、境界条件が求められない場合、３次元ＳＰＨ計算処理部３は、２次元伝播計算における次のステップの計算結果を２次元伝播計算処理部２から取得する。その後、３次元ＳＰＨ計算処理部３は、取得した２次元伝播計算における次のステップの計算結果から境界条件を求める。そして、３次元ＳＰＨ計算処理部３は、ステップ＃ｎ＋１の速度、密度、圧力及び位置及び境界条件を用いて、ステップ＃ｎ＋２の速度、密度、圧力及び位置の算出を行う。

以下では、図１５を参照して、２次元伝播計算及び３次元粒子法計算の処理の全体的な流れについて説明する。図１５は、実施例３に係るシミュレーション装置による解析処理のフローチャートである。

２次元伝播計算処理部２は、２次元伝播計算を１ステップ進める（ステップＳ４０１）。そして、２次元伝播計算処理部２は、現在のステップの計算結果を３次元ＳＰＨ計算処理部３へ出力する。

３次元ＳＰＨ計算処理部３は、現在までに受信した２次元伝播計算の計算結果から、３次元粒子法計算の開始条件が満たされたか否かを判定する（ステップＳ４０２）。例えば、３次元ＳＰＨ計算処理部３は、領域１００の境界における水位が５０ｃｍを超えることを、３次元粒子法計算を開始する条件として記憶している。３次元粒子法計算の開始条件が満たされていない場合（ステップＳ４０２：否定）、２次元伝播計算処理部２は、ステップＳ４０１に戻る。

これに対して、３次元粒子法計算の開始条件が満たされた場合（ステップＳ４０２：肯定）、３次元ＳＰＨ計算処理部３は、２次元伝播計算による水位及び流量から３次元粒子法計算における位置及び速度の初期条件を生成する（ステップＳ４０３）。

そして、３次元ＳＰＨ計算処理部３は、１ステップ後の３次元粒子法計算を実施し、１ステップ後の各粒子の速度、密度、圧力及び位置を算出する（ステップＳ４０４）。

その後、３次元ＳＰＨ計算処理部３は、現ステップまでの計算結果から、３次元粒子法計算の終了条件を満たすか否かを判定する（ステップＳ４０５）。

３次元粒子法計算の終了条件を満たしていない場合（ステップＳ４０５：否定）、３次元ＳＰＨ計算処理部３は、次のステップの境界条件を設定するための２次元伝播計算による水位及び流量のデータが存在するか否かを判定する（ステップＳ４０６）。

次のステップの境界条件を設定するための２次元伝播計算による水位及び流量のデータが存在しない場合（ステップＳ４０６：否定）、２次元伝播計算処理部２は、２次元伝播計算を１ステップ進める（ステップＳ４０７）。そして、２次元伝播計算処理部２は、現在のステップの計算結果を３次元ＳＰＨ計算処理部３へ出力する。

これに対して、次のステップの境界条件を設定するための水位及び流量のデータが存在する場合（ステップＳ４０６：肯定）、３次元ＳＰＨ計算処理部３は、ステップＳ４０８へ進む。

そして、３次元ＳＰＨ計算処理部３は、２次元伝播計算による水位及び流量から３次元粒子法計算における境界条件を生成する（ステップＳ４０８）。その後、３次元ＳＰＨ計算処理部３は、ステップＳ４０４へ戻る。

一方、３次元粒子法計算の終了条件を満たされた場合（ステップＳ４０５：肯定）、３次元ＳＰＨ計算処理部３は、３次元粒子法計算を終了する。

ここで、図１５では、３次元粒子法計算において２次元伝播計算の計算結果が要求されてから、２次元伝播計算を行うように説明したが、２次元伝播計算処理部２は、３次元粒子法計算の進行状況に関らず２次元伝播計算を進めてよい。

以上に説明したように、本実施例に係るシミュレーション装置は、２次元伝播計算の終了を待たずに３次元粒子法計算を開始し、それらを並行して実行する。この場合、例えば、２次元伝播計算を行う計算機と３次元粒子法計算を行う計算機を分けることができる。そして、３次元粒子法計算を行う計算機は、２次元伝播計算を行う計算機が出力した２次元伝播計算の結果のファイルを用いて、３次元粒子法計算を実行できる。これにより、全体の計算時間を減らすことができる。

また、上記実施例で説明した各種の処理は、あらかじめ用意されたプログラムをコンピュータで実行することによって実現することができる。そこで、以下では、図１６を用いて、図１に示したシミュレーション装置と同様の機能を有するシミュレーションプログラムを実行するコンピュータの一例を説明する。

図１６は、シミュレーションプログラムを実行するコンピュータを示す図である。図１６に示すように、コンピュータ１０００は、キャッシュ１０１０と、ＣＰＵ（Central Processing Unit）１０２０と、ＨＤＤ（Hard Disk Drive）１０３０と、ＲＡＭ（Random Access Memory）１０４０及びＲＯＭ（Read Only Memory）１０５０を有する。キャッシュ１０１０、ＣＰＵ１０２０、ＨＤＤ１０３０、ＲＡＭ１０４０及びＲＯＭ１０５０は、それぞれバスによって接続されている。

ＨＤＤ１０３０には、図１に示したシミュレーション装置と同様の機能を発揮する各種のシミュレーションプログラム１０３１が予め記憶されている。

そして、ＣＰＵ１０２０は、これらのシミュレーションプログラム１０３１を読み出して実行する。これにより、図１６に示すように、シミュレーションプログラム１０３１は、シミュレーションプロセス１０２１になる。

なお、上記したシミュレーションプログラム１０３１については、必ずしもＨＤＤ１０３０に記憶させなくてもよい。例えば、コンピュータ１０００に挿入されるフレキシブルディスク、ＣＤ（Compact Disk）−ＲＯＭ、ＤＶＤ（Digital Versatile Disc）、光磁気ディスク、ＩＣ（Integrated Circuit）カードなどの「可搬用の物理媒体」にシミュレーションプログラム１０３１を記憶させてもよい。または、コンピュータ１０００の内外に備えられるハードディスクドライブ（ＨＤＤ）などの「固定用の物理媒体」にシミュレーションプログラム１０３１を記憶させてもよい。または、公衆回線、インターネット、ＬＡＮ（Local Area Network）、ＷＡＮ（Wide Area Network）などを介してコンピュータ１０００に接続される「他のコンピュータ（またはサーバ）」にシミュレーションプログラム１０３１を記憶させてもよい。そして、コンピュータ１０００は、上述したフレキシブルディスク等から各プログラムを読み出して実行するようにしてもよい。

１ ユーザインタフェース
２ ２次元伝播計算処理部
３ ３次元ＳＰＨ計算処理部
４ 結果記憶部
２１ 制御部
２２ 水位算出部
２３ 流量算出部
３１ 制御部
３２ 速度算出部
３３ 密度算出部
３４ 圧力算出部
３５ 位置算出部
３６ 粒子生成部

さらに、制御部３１は、ステップ＃ｎにおける各着目粒子の最新の速度ｖ、密度ρ、圧力Ｐ及び位置ｒの各情報の入力を受けると、受信したステップ＃ｎにおける各情報を結果記憶部４へ記憶させる。

そして、３次元ＳＰＨ計算処理部３は、１ステップ後の３次元粒子法計算を実施し、各粒子の速度、密度、圧力及び位置を算出する（ステップＳ４０４）。

Claims (8)

1. 入力されたデータを基に、第１領域の２次元平面の各点における連続体の流量及び水位を算出し、
   前記第１領域に含まれる第２領域における前記連続体を粒子の集合として表現し、算出した前記流量及び前記水位を基に、各前記粒子の状態を３次元解析する
   処理をコンピュータに実行させることを特徴とするシミュレーションプログラム。
2. 算出した前記流量及び前記水位を基に初期条件及び境界条件を生成し、生成した前記初期条件及び前記境界条件を基に、各前記粒子の状態を解析する処理をコンピュータに実行させることを特徴とする請求項１に記載のシミュレーションプログラム。
3. 前記第１領域の２次元平面の各点における連続体の流量及び水位の初期値の入力を受け、
   前記点毎に、現ステップにおける水位の時間微分を算出し、前記現ステップの水位を、算出した前記水位の時間微分を基に時間積分することで次のステップの水位を算出し、
   前記点毎に、現ステップにおける流量の時間微分を算出し、前記現ステップの流量を、算出した前記流量の時間微分を基に時間積分することで次のステップの流量を算出し、
   前記次のステップの水位及び速度の算出を第１の所定ステップまで繰り返す
   処理をコンピュータに実行させることを特徴とする請求項１又は請求項２に記載のシミュレーションプログラム。
4. 前記粒子毎に、現ステップにおける加速度を運動方程式に基づき算出し、前記現ステップの速度を、算出した前記加速度を基に時間積分することで次のステップの速度を算出し、
   前記粒子毎に、現ステップにおける密度の時間微分を連続方程式に基づき算出し、前記現ステップの密度を、算出した前記密度の時間微分を基に時間積分することで次のステップの密度を算出し、
   前記粒子毎に、算出した次のステップの密度を基に、状態方程式を用いて次のステップの圧力を算出し、
   前記粒子毎に、現ステップにおける位置座標を、算出した次のステップの速度で時間積分して次のステップの位置座標を算出し、
   前記次のステップの速度、密度、圧力及び位置座標の算出を第２の所定ステップまで繰り返すことで各前記粒子の状態を解析する
   処理をコンピュータに実行させることを特徴とする請求項１〜３のいずれか一つに記載のシミュレーションプログラム。
5. 前記第１領域の２次元平面の各点における連続体の速度を、前記第２領域の境界から一定距離にある領域境界に存在する前記粒子の速度として用いる処理をコンピュータに実行させることを特徴とする請求項１〜４のいずれか一つに記載のシミュレーションプログラム。
6. 前記第２領域の境界から一定距離にある領域境界を越えて流出した粒子を消去し、前記連続体が前記領域境界に流入した体積が所定体積を超えた場合、粒子を新たに生成することで、各前記粒子の状態を解析する処理をコンピュータに実行させることを特徴とする請求項１〜５のいずれか一つに記載のシミュレーションプログラム。
7. 入力されたデータを基に、第１領域の２次元平面の各点における連続体の流量及び水位を算出し、
   前記第１領域に含まれる第２領域における前記連続体を粒子の集合として表現し、算出した前記流量及び前記水位を基に、各前記粒子の状態を３次元解析し、
   前記解析の結果をメモリに記憶させる
   処理をコンピュータに実行させることを特徴とするシミュレーション方法。
8. 入力されたデータを基に、第１領域の２次元平面の各点における連続体の流量及び波高を求める２次元解析処理部と、
   前記２次元解析処理部により算出されたデータを基に、前記連続体を粒子の集合として表現する場合の各粒子の状態を解析する３次元解析を、前記第１領域に含まれる第２領域に対して行う３次元解析処理部と、
   を備えたことを特徴とするシミュレーション装置。

Images