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

シミュレーションプログラム、シミュレーション方法及びシミュレーション装置 Download PDF

Info

Publication number
JP6048062B2
JP6048062B2 JP2012231224A JP2012231224A JP6048062B2 JP 6048062 B2 JP6048062 B2 JP 6048062B2 JP 2012231224 A JP2012231224 A JP 2012231224A JP 2012231224 A JP2012231224 A JP 2012231224A JP 6048062 B2 JP6048062 B2 JP 6048062B2
Authority
JP
Japan
Prior art keywords
particle
calculation
flow rate
dimensional
region
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2012231224A
Other languages
English (en)
Other versions
JP2014081900A (ja
Inventor
多聞 諏訪
多聞 諏訪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujitsu Ltd
Original Assignee
Fujitsu Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fujitsu Ltd filed Critical Fujitsu Ltd
Priority to JP2012231224A priority Critical patent/JP6048062B2/ja
Priority to US13/967,468 priority patent/US9557203B2/en
Publication of JP2014081900A publication Critical patent/JP2014081900A/ja
Application granted granted Critical
Publication of JP6048062B2 publication Critical patent/JP6048062B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01FMEASURING VOLUME, VOLUME FLOW, MASS FLOW OR LIQUID LEVEL; METERING BY VOLUME
    • G01F1/00Measuring the volume flow or mass flow of fluid or fluent solid material wherein the fluid passes through a meter in a continuous flow
    • G01F1/76Devices for measuring mass flow of a fluid or a fluent solid material

Landscapes

  • Physics & Mathematics (AREA)
  • Fluid Mechanics (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Description

本発明は、シミュレーションプログラム、シミュレーション方法及びシミュレーション装置に関する。
従来、数値計算を用いて水や空気といった流体の流れを調べる連続体の運動の解析、すなわち、流体解析において、粒子法と呼ばれる技術が提案されている。具体的には、粒子法とは、連続体の運動を有限の数の粒子の運動として解析する手法である。現在提案されている代表的な粒子法としては、SPH(Smoothed Particles Hydrodynamics)法やMPS(Moving Particles Semi-implicit)法といったものがある。以下では、水や空気といった流体を「連続体」と呼ぶ場合がある(例えば、非特許文献1参照)。
粒子法における標準的な手法としては、ある粒子に対して予め領域(以下では、「影響領域」と呼ぶ。)を設定しておき、その影響領域内に存在する相手粒子のみからの相互作用を計算することで、その粒子に加えられる力を求める方法が従来知られている。
SPH法を用いて連続体を表現する際の特徴は、複数の粒子の物理量をカーネル関数と呼ばれる重み関数を用いて平滑化近似して、基礎方程式を離散化することにある。この平滑化処理によって、メッシュ点上で物理量を評価するといったきわめて煩雑なメッシュ操作の計算処理を行わなくてよくなる。そのため、SPH法は、自由表面問題やマルチフィジックス問題を扱うのに適しているといえる。具体的には、SPH法は、例えば、海の波が護岸に衝突や越波するときの流量や衝撃圧を推定するための手法として適しているといえる。
ただし、SPH法による計算を用いて広範囲の計算を行った場合、計算量が大きくなりすぎる。そのため、SPH法を用いた場合には、大きくとも数km四方程度の限られた領域でのみの計算を行うことになる。そのため、津波の伝播を考えた場合、津波計算に用いる震源直上の津波波源から沿岸部まで波が伝播する数100km以上に広がる過程をSPH法で扱うことは困難である。
また、地形データ及び構造データと組み合わせて流体のシミュレーションを行う従来技術(例えば、特許文献1参照)があるが、この方法を用いても数100km以上といった広域を扱うことは困難である。
一方、地震による地殻変動から、津波発生を経て、震源から海岸に至る津波の伝播までの過程を2次元の伝播シミュレータによって計算する従来技術も提案されている(例えば、特許文献2、3、非特許文献2参照)。これらの手法は、2次元である地図上の各点で波高及び速度等の情報が設定されており、津波が発生した場合の沿岸における波高及び浸水域等を数値計算によって予測することが可能である。
なお、低い次元で計算した結果を高い次元の計算の初期条件及び境界条件として与えるという従来技術が提案されている。例えば、加震源と受振構造物を三次元有限要素法で扱い、その間の地盤を介して伝わる振動を簡易的なモデルで求めるといった従来技術がある(例えば、特許文献4参照)。また、1次元のBoussinesq波動モデルと2次元のSPH法の計算を組み合わせる従来技術がある(例えば、非特許文献3参照)。1次元のBoussinesq波動モデルと2次元のSPH法の計算を組み合わせる従来技術では、例えば、以下のような処理が行われる。SPH法で計算を行う領域とBoussinesq法で計算を行う領域とが重なり合う領域であるバッファ領域を設定する。バッファ領域の中では、Boussinesq法の計算結果として得られた速度を強制的にSPH法の各粒子に設定する。また、バッファ領域の外においては、SPH粒子は流体の運動方程式を解いた結果として得られる加速度に基づいて、速度の更新が行われる。ただし、バッファ領域は予め決められており、SPH粒子の生成や消滅は行われない。
特開2005−115701号公報 特開2008−089316号公報 特開2010−054460号公報 特開2004−219237号公報
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
しかしながら、2次元シミュレーションのみを用いた従来技術では、3次元性が重要になる港湾部などの複雑形状を扱うことが困難となることや、建築物にかかる波圧及び波力を算出することが困難であるといった問題がある。
また、加震源と受振構造物を三次元有限要素法で扱い、その間の地盤を介して伝わる振動を簡易的なモデルで求めるといった従来技術は、地盤中の振動に適用されるものであり、そのまま津波シミュレーションに用いることは困難である。
また、Boussinesq波動モデルをSPH法の計算に用いる従来技術では、津波のように大量の流体が移動を伴う場合には、バッファ領域が狭いと、バッファ領域の粒子が全てSPH領域に移動してしまい、バッファ領域での計算ができなくなるおそれがある。一方、津波の流体移動を考慮しバッファ領域の幅を決めると、バッファ領域が非常に大きくなり、計算量及びメモリ量が増大しすぎてしまい、実際に処理することが困難となる。さらに、この従来技術では、1次元と2次元との組合せを対象としており、2次元と3次元とを組み合わせることは考慮されておらず、津波の伝播計算にそのまま用いることは困難である。
開示の技術は、上記に鑑みてなされたものであって、連続体の移動を精度良くシミュレートするシミュレーションプログラム、シミュレーション方法及びシミュレーション装置を提供することを目的とする。
本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置は、一つの態様において、入力されたデータを基に、第1領域の2次元平面の各点における連続体の流量及び水位を算出し、前記第1領域に含まれる第2領域における前記連続体を粒子の集合として表現し、算出した前記水位及び前記流量を基に、前記第2領域の境界から一定距離にある領域境界を越えて流出した粒子を消去し、前記連続体が前記領域境界に流入した体積が所定体積を超えた場合、粒子を新たに生成して、各前記粒子の状態を3次元解析する処理をコンピュータに実行させる。
本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置の一つの態様によれば、連続体の移動を精度良くシミュレートすることができるという効果を奏する。
図1は、実施例1に係るシミュレーション装置のブロック図である。 図2は、2次元伝播計算及び3次元粒子法計算を行う領域の例を示す図である。 図3は、実施例1に係るシミュレーション装置による解析処理のフローチャートである。 図4は、実施例2に係る2次元伝播計算処理部のブロック図である。 図5は、実施例2に係る3次元SPH計算処理部のブロック図である。 図6は、初期値の設定を模式的に表した図である。 図7は、影響範囲及び近傍粒子について説明するための図である。 図8は、粒子の流出を模式的に表した図である。 図9は、境界条件算出における補間を説明するための図である。 図10は、境界からの反発力を説明するための図である。 図11は、粒子の流入を模式的に表した図である。 図12は、実施例2に係る2次元伝播計算の処理のフローチャートである。 図13は、実施例2に係る3次元SPH計算の処理のフローチャートである。 図14は、実施例2に係る3次元SPH計算の処理のフローチャートである。 図15は、実施例3に係るシミュレーション装置による解析処理のフローチャートである。 図16は、シミュレーションプログラムを実行するコンピュータを示す図である。
以下に、本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置の実施例を図面に基づいて詳細に説明する。なお、以下の実施例により本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置が限定されるものではない。また、以下では、津波発生時における海水の動きを解析するシミュレーションについて説明するが、連続体の動きを解析するシミュレーションであれば各実施例を適応できる。特に、津波のような連続体の移動量が大きい場合のシミュレーションに適応することが好ましい。
図1は、実施例1に係るシミュレーション装置のブロック図である。図1に示すように、本実施例に係るシミュレーション装置は、ユーザインタフェース1、2次元伝播計算処理部2、3次元SPH計算処理部3及び結果記憶部4を有している。
ユーザインタフェース1は、シミュレーション装置を使用する操作者からの数値の入力、操作者へのシミュレーション結果の出力を行なうインタフェース装置であり、具体的にはキーボード等の入力装置及び表示装置等の出力装置である。
操作者は、ユーザインタフェース1を用いて、例えば、図2に示すように領域100及び領域101を設定する。図2は、2次元伝播計算及び3次元粒子法計算を行う領域を示す図である。領域101は、後述する2次元伝播計算を行う領域である。領域101が、「第1領域」の一例にあたる。領域100は、後述する3次元粒子法計算を行う領域である。領域100が、「第2領域」の一例にあたる。このように、2次元平面の各点における連続体の流量及び水位を算出する2次元伝播計算は、大きな領域に対して行い、3次元SPH計算は、沿岸部などの3次元的に入り組んだ構造を含む小さな領域に対して行う。ここで、領域100と領域101との境界はある程度水深がある場所であることが好ましい。これは、水深が浅くなると、海水の流れが海底や沿岸の影響を受け易くなるため2次元伝播計算の精度が落ちるからである。2次元伝播計算の計算結果を3次元粒子法計算に用いるためには、2次元伝播計算の精度が、3次元SPH計算で計算を行った場合の精度と同程度であることが好ましい。そこで、ある程度の深さがある地点を境界として設定することが好ましい。具体的には、例えば、30m〜50mの水深を有する地点を境界として設定することが好ましい。
さらに、操作者は、初期データの入力を行う。本実施例では、図2の領域101を平面視した2次元平面を所定の格子に分割する。2次元伝播計算では、この分割した格子における海水の動きを解析する。そして、初期データとは、領域101を表す座標をxy座標として、分割した各格子における、水位η、x方向の流量M(以下では、単に「流量M」という場合がある。)及びy方向の流量N(以下では、単に「流量N」という場合がある。)の初期値である。ここで、水位ηとは、静水面からの水位上昇量である。また、流量とは、速度と全水深との積であり、x方向の流量Mは、M=v×D(vはx方向の速度、Dは全水深であり、D=水位η+静水深hである。)として表される。また、y方向の流量Nは、N=v×D(vはy方向の速度である。)として表される。全水深と流量から速度を得ることができ、また、逆に全水深と速度から流量を得ることもできるので、流量と速度のどちらも本質的には同等の情報を含んでいるとみなせる。二次元伝播計算では流量の方が扱いやすいため、流量を用いるが、三次元粒子法計算では速度の方が扱いやすいため、流量の情報を必要に応じて速度に変換することになる。速度vと速度vとを合成した速度が各格子における流速になる。
そして、2次元伝播計算処理部2は、入力された水位η、流量M及び流量Nの初期値を用いて、次のステップの各格子における水位η、流量M及び流量Nを算出する。さらに、2次元伝播計算処理部2は、順々に各ステップの各格子における水位η、流量M及び流量Nを算出していく。2次元伝播計算処理部2は、予め設定された最終の状態になる最終ステップまでの各ステップの各格子における水位η、流量M及び流量Nの算出を繰り返す。この予め設定された最終状態になったときのステップが、「第1の所定ステップ」にあたる。本実施例では、津波が陸上の最深部まで達し引き始めた時点を最終状態とする。
そして、2次元伝播計算処理部2は、算出したステップ毎の各格子における水位η、流量M及び流量Nのデータを3次元SPH計算処理部3へ出力する。
3次元SPH計算処理部3は、ステップ毎の各格子における水位η、流量M及び流量Nのデータの入力を2次元伝播計算処理部2から受ける。
そして、3次元SPH計算処理部3は、受信したステップ毎の各格子における水位η、流量M及び流量Nのデータのうち、3次元粒子法計算を開始する条件を満たすデータを有する最も時刻の早いステップを特定する。3次元SPH計算処理部3は、例えば、図2の領域100の境界における水位ηが50cmを超えた時点を、3次元粒子法計算を開始する条件として記憶している。
そして、3次元SPH計算処理部3は、特定した3次元粒子法計算を開始するステップにおける各格子における水位ηを用いて領域100のその時点での水深に合わせて粒子を配置する。ここで、水深とは水位ηで表される頂点から海底までの距離である。また、粒子は、どのような大きさにすることもできるが、直径が10cm〜1mであることが好ましい。直径が小さいほど精度の良い解析ができるが計算量が増え、直径が大きいほど計算量が減るが解析の精度が落ちる。そのため、粒子の大きさは、使用できる計算資源と供給される精度に応じて決定することが好ましい。粒子は、流体をモデル化したものであり、粒子を表現するために要求されるデータは、例えば、中心座標、速度、影響半径、密度、質量、粘性度などである。このうち、3次元SPH計算処理部3は、影響半径、質量、粘性度などを予め記憶している。
さらに、3次元SPH計算処理部3は、3次元粒子法計算を開始するステップの各格子における流量M及び流量Nから各粒子の速度の初期値である速度vを求める。ここで、「v」は粒子の速度を表すベクトルである。また、3次元SPH計算処理部3は、各格子から各粒子の位置座標の初期値である位置rを求める。ここで、「r」は粒子の位置を表す位置ベクトルである。さらに、3次元SPH計算処理部3は、圧力及び密度の初期値である圧力P及び密度ρとして海水の一般的な圧力及び密度をそれぞれ用いる。
そして、3次元SPH計算処理部3は、速度v、位置r、圧力P及び密度ρを用いて、次のステップの各粒子の速度v、位置r、圧力P及び密度ρを求める。このとき、3次元SPH計算処理部3は、図2の領域100における境界から一定の距離Δxの範囲を境界領域として記憶している。そして、3次元SPH計算処理部3は、各ステップにおける速度v、位置r、圧力P及び密度ρを算出する場合に、境界領域における粒子の速度vなどの境界条件を、受信したステップ毎の各格子における流量M及び流量Nから求める。ここで、3次元粒子法計算を行うステップは、2次元伝播計算を行うステップよりも細かい。そこで、3次元SPH計算処理部3は、境界条件として用いる流量M及び流量Nを補間しながら境界条件を求めていく。
このように、3次元SPH計算処理部3は、予め設定された最終の状態になるまでの各ステップの各粒子における速度v、位置r、圧力P及び密度ρの算出を繰り返す。この予め設定された最終状態になったときのステップが、「第2の所定ステップ」にあたる。本実施例では、領域100において津波が陸上の最深部まで達し引き始めた時点を最終状態とする。
そして、3次元SPH計算処理部3は、算出したステップ毎の各粒子における速度v、位置r、圧力P及び密度ρのデータを結果記憶部4へ出力する。
結果記憶部4は、3次元SPH計算処理部3から受信した解析結果であるステップ毎の各粒子における速度v、位置r、圧力P及び密度ρをメモリやハードディスク等に格納する。
次に、図3を参照して本実施例に係る解析処理の流れについて説明する。図3は、実施例1に係るシミュレーション装置による解析処理のフローチャートである。
2次元伝播計算処理部2は、ユーザインタフェース1を用いて入力された入力データである各格子における、水位η、x方向の流量M及びy方向の流量Nの初期値を取得する(ステップS1)。
2次元伝播計算処理部2は、水位η、x方向の流量M及びy方向の流量Nの初期値を用いて、各ステップにおける2次元伝播計算を最終ステップまで実行する(ステップS2)。
3次元SPH計算処理部3は、ステップ毎の各格子における、水位η、x方向の流量M及びy方向の流量Nの入力を2次元伝播計算処理部2から受ける。次に、3次元SPH計算処理部3は、3次元粒子法計算を開始するステップを特定する。そして、3次元SPH計算処理部3は、特定したステップにおける2次元伝播計算の水位η、x方向の流量M及びy方向の流量Nのデータから3次元粒子法計算の初期条件を生成する(ステップS3)。
3次元SPH計算処理部3は、現ステップの速度v、位置r、圧力P及び密度ρを用いて、1ステップ後の3次元粒子法計算を実行し、次のステップの速度v、位置r、圧力P及び密度ρを算出する(ステップS4)。
3次元SPH計算処理部3は、3次元粒子法計算の終了条件を満たすか否かを判定する(ステップS5)。終了条件を満たしていない場合(ステップS5:否定)、3次元SPH計算処理部3は、2次元伝播計算処理部2から受信した各格子における、2次元伝播計算の水位η、x方向の流量M及びy方向の流量Nから次のステップの境界条件を生成する(ステップS6)。その後、3次元SPH計算処理部3は、ステップS4に戻る。
これに対して、終了条件を満たしている場合(ステップS5:肯定)、3次元SPH計算処理部3は、解析結果を結果記憶部4に格納し、解析処理を終了する。
以上に説明したように、本実施例に係るシミュレーション装置は、連続体の広い移動に対するシミュレーションにおいて、広い範囲では2次元伝播計算を行い、港湾部や都市部といった3次元の形状の影響が大きい領域では3次元粒子法計算を行う。これにより、津波のように大量の流体の移動が伴う事象においても、破綻無く計算を行うことができ、計算量を抑えながら、精度の良いシミュレーションを行うことができる。
次に、実施例2について説明する。本実施例に係るシミュレーション装置は、全体的な処理は実施例1と同様であるが、処理方法を具体的にしたものである。そこで、以下では、2次元伝播計算及び3次元粒子法計算の処理について主に説明する。また、以下で特に説明の無い限りは、図1と同じ符号を有する各部は同じ機能を有しているものとする。
図4は、実施例2に係る2次元伝播計算処理部のブロック図である。図4に示すように、2次元伝播計算処理部2は、制御部21、水位算出部22及び流量算出部23を有している。
制御部21は、領域101を平面視した2次元平面を分割した格子を予め記憶している。さらに、制御部21は、領域101の点を表す座標をxy座標としたときの、分割した各格子における、水位η、x方向の流量M及びy方向の速度Yの初期値である水位η、流量M及び流量Nの入力をユーザインタフェース1から受ける。
制御部21は、各格子における水位η、流量M及び流量Nを水位算出部22及び流量算出部23に出力する。その後、制御部21は、ステップ1での各格子における水位ηの入力を水位算出部22から受ける。また、制御部21は、ステップ1での各格子における流量M及び流量Nの入力を流量算出部23から受ける。以下では、初期値である水位η、流量M及び流量Nで表される初期状態をステップ0とし、その後の各ステップをステップ1,2,3,・・・とする。そして、各ステップをまとめて、ステップk(k=0,1,2,・・・)と表す。
その後、制御部21は、受信したステップk(kは自然数)での水位η、流量M及び流量Nを水位算出部22及び流量算出部23に出力し、ステップk+1での水位η、流量M及び流量Nの入力を水位算出部22及び流量算出部23から受けることを繰り返す。
また、制御部21は、最終状態の条件の入力をユーザインタフェース1から受ける。本実施例では、制御部21が、最終状態の条件として、津波が陸上の最深部まで達し引き始めた時点という条件の入力を受ける。そして、制御部21は、ステップkまでの計算結果を受信して、ステップ1〜kまでの各格子における水位η、流量M及び流量Nを用いて、最終状態の条件を満たしているか否かを判定する。最終状態の条件を満たしている場合、制御部21は、現ステップであるステップkを最終ステップとする。その場合、制御部21は、受信したステップkでの各格子における水位η、流量M及び流量Nの水位算出部22及び流量算出部23への出力を停止する。
その後、制御部21は、ステップ1〜kの各ステップでの各格子における水位η、流量M及び流量Nを3次元SPH計算処理部3へ出力する。また、制御部21は、各格子の位置の情報を3次元SPH計算処理部3へ出力する。
水位算出部22は、次に示す式(1)で表される連続方程式を記憶している。
Figure 0006048062
水位算出部22は、ステップk+1では、ステップkの各格子における水位η、流量M及び流量Nの入力を制御部21から受ける。ここで、x方向の流量Mは、M=v×D(vはx方向の速度の関数、Dは全水深であり、水位η+静水深hである。)として表される。また、y方向の流量Nは、M=v×D(vはy方向の速度の関数である。)として表される。
そこで、水位算出部22は、ステップk+1の計算において、ステップkにおける流量M及び流量Nを式(1)に代入して、ステップkにおける水位の時間微分である∂η/∂tを求める。
さらに、水位算出部22は、次に示す式(2)で表される時間積分の式を記憶している。ここで、Δtは、各ステップ間の時間を表している。
Figure 0006048062
水位算出部22は、ステップk+1の計算においては、ステップkにおける水位の時間微分である∂η/∂tを式(2)に代入し、ステップkの各格子における水位ηをステップnにおける水位の時間微分である∂η/∂tに基づいて時間積分する。これにより、水位算出部22は、ステップk+1における水位ηを算出する。
そして、水位算出部22は、算出したステップkにおける水位ηを制御部21へ出力する。
流量算出部23は、次に示す式(3)及び(4)で表される運動方程式を記憶している。なお、gは重力加速度、Dは全水深(D=h+η、ただし、hは静水深)、nは海底摩擦を表す粗度係数である。
Figure 0006048062
Figure 0006048062
流量算出部23は、ステップk+1では、ステップkの各格子における水位η、流量M及び流量Nの入力を制御部21から受ける。
流量算出部23は、ステップk+1の計算においては、ステップkにおける水位η、流量M及び流量Nを式(3)に代入して、ステップkにおけるx方向の速度の時間微分である∂M/∂tを求める。
また、流量算出部23は、ステップk+1の計算においては、ステップkにおける水位η、流量M及び流量Nを式(4)に代入して、ステップkにおけるy方向の速度の時間微分である∂N/∂tを求める。
さらに、水位算出部22は、次に示す式(5)及び(6)で表される時間積分の式を記憶している。ここで、Δtは、各ステップ間の時間を表している。
Figure 0006048062
Figure 0006048062
流量算出部23は、ステップk+1の計算においては、ステップkにおけるx方向の速度の時間微分である∂M/∂tを式(5)に代入し、ステップkの各格子における流量Mをステップkにおけるx方向の速度の時間微分である∂M/∂tに基づいて時間積分する。これにより、流量算出部23は、ステップk+1における流量Mを算出する。
また、流量算出部23は、ステップk+1の計算においては、ステップkにおけるy方向の速度の時間微分である∂N/∂tを式(6)に代入し、ステップkの各格子における流量Nをステップkにおけるy方向の速度の時間微分である∂N/∂tに基づいて時間積分する。これにより、流量算出部23は、ステップk+1における流量Nを算出する。
そして、流量算出部23は、算出したステップk+1における流量M及び流量Nを制御部21へ出力する。
次に、図5を参照して本実施例に係る3次元SPH計算処理部について説明する。図5は、実施例2に係る3次元SPH計算処理部のブロック図である。図5に示すように、3次元SPH計算処理部3は、制御部31、速度算出部32、密度算出部33、圧力算出部34、位置算出部35及び粒子生成部36を有している。ここで、3次元粒子計算におけるステップは、2次元伝播計算におけるステップよりも細かい。そこで、以下では、3次元粒子計算におけるステップを2次元伝播計算におけるステップと区別するためステップ#nと記載する。
制御部31は、例えば、図2の領域100の境界における水位ηが50cmを超えた時点を、3次元粒子法計算を開始する条件として記憶している。また、制御部31は、粒子と粒子の運動の境界条件を設定するための固定境界要素を予め記憶している。固定境界要素は、海底面、地面及び建造物の表面などを微小部分に分割した面要素をモデル化したものである。固定境界要素は、例えば、微小な円盤の集合として境界全体をモデル化して、各境界要素に中心座標、面要素の法線ベクトル、面要素の面積などを含んでも良い。また、固定境界要素として、ポリゴンの集合として境界全体を表現して、各境界要素に複数の頂点の位置座標を持たせても良い。
制御部31は、2次元伝播計算の計算結果であるステップ毎の各格子における水位η、流量M及び流量Nの入力を2次元伝播計算処理部2から受ける。また、制御部31は、各格子の位置の情報の入力を2次元伝播計算処理部2から受ける。
そして、制御部31は、各格子の位置の情報及びステップ毎の各格子における水位ηの中から、境界における水位ηが50cmを超えた時点のステップをステップ#0として特定する。
そして、制御部31は、ステップ#0での各格子における水位η、流量M及び流量Nを抽出する。
ここで、3次元粒子法計算は連続体の動作をより詳細に解析するために、2次元伝播計算で用いられた格子よりも細かい精度で計算を行う。そこで、制御部31は、3次元粒子法計算で用いる点を決めるための格子をより細かく分割する数量を予め記憶している。そして、制御部31は、ステップ#0での各格子における水位ηの内挿を行い、3次元粒子法計算で用いる位置座標における水位ηを求める。次に、制御部31は、3次元粒子法計算で用いる各位置座標における水位ηに静水深hを加算して、各位置座標における全水深Dを求める。また、制御部31は、3次元粒子法計算で用いる位置座標におけるステップ#0での流量M及び流量Nを内挿により求める。そして、制御部31は、各位置座標でのステップ#0での流量M及び流量Nを全水深で除算し、各位置座標での速度の初期値vを求める。
そして、制御部31は、各位置座標における全水深Dを用いて、各位置座標に予め決められたサイズの粒子を配置する。例えば、粒子の直径が1mと決められており、ある位置座標における全水深が50mの場合、制御部31は、その位置座標に50個の粒子を積み重ねて配置する。これにより、各位置座標における粒子の位置の初期値が決定する。
さらに、制御部31は、各位置座標に配置した各粒子の速度の初期値として、算出したその位置座標における速度vを割り当てる。図6は、初期値の設定を模式的に表した図である。図6は、図2に示す領域100の断面の一部を模式的に現している。図6の縦軸は水深を表し、横軸は静水面と平行な特定の方向を表している。図6に示すように、制御部31は、各位置座標において粒子200を水深方向に積み上げて配置する。そして、制御部31は、矢印201で示されるように、矢印201の大きさに応じた速度を各粒子に割り当てる。図6では、例えば、横軸が増える方向を沖に向かう方向とすると、横軸が減る方向に向かう矢印201は、粒子が陸に向かっていることを表している。そして、矢印201が大きいほど粒子の速度が速いことを表している。
各粒子の位置及び速度の初期値の設定が終わると、制御部31は、各粒子の影響範囲内にある他の粒子のリストを作成する。以下では、各粒子の影響範囲内にある他の粒子を「近傍粒子」と呼ぶ。また、近傍粒子のリストを「近傍粒子リスト」と呼ぶ。
図7は、影響範囲及び近傍粒子について説明するための図である。ここで、図7を参照して、各粒子の影響範囲及び近傍粒子について説明する。ここでは、着目粒子と他の粒子との距離が2hSPH以内の場合に、相互に影響を与え合うものとして説明する。ここで、影響とは、例えば、粒子が移動する場合に、他の粒子へ力を加えるなどのことを言う。
ここでは、図7の粒子301を例に説明する。粒子301の影響範囲302は、粒子301を中心として半径2hSPHの円の中に含まれる範囲である。この半径2hSPHは、粒子301の影響半径である。そして、影響範囲302の中に含まれる他の粒子が、粒子301の近傍粒子である。図7では、例えば、粒子303のように、斜線で表されている粒子が近傍粒子である。すなわち、粒子301が影響を与えることになる他の粒子は、斜線で表される近傍粒子ということになる。逆に、斜線で表される近傍粒子からのみ、粒子301は影響を受けることになる。ここで、近傍粒子が数十個程度含まれるように半径2hSPHを設定すれば、着目した粒子と他の粒子間の相互作用を解析するには十分であることが経験的に知られている。
また、粒子301の位置を表す位置ベクトルをrとし、粒子303の位置を表す位置ベクトルをrとすると、粒子301と粒子303との間の距離は|r−r|となる。すなわち、|r−r|≦2hを満たす場合に、粒子303は粒子301の近傍粒子とされる。
制御部31は、例えば図7の粒子301に対する近傍リストを作る場合、まず、影響範囲302に含まれる粒子303のような斜線で表される各粒子を粒子301の近傍粒子として抽出する。例えば、制御部31は、着目粒子である粒子301の位置と他の粒子の位置から粒子301と各他の粒子との距離を求め、その距離が2hSPH以下となる粒子を近傍粒子として抽出する。
図5に戻って、制御部31は、近傍粒子の抽出を全ての粒子に対して行う。そして、制御部31は、各着目粒子と抽出した各着目粒子の近傍粒子とが対応する近傍粒子リストを作成する。そして、制御部31は、作成した近傍粒子リストを速度算出部32、密度算出部33、圧力算出部34、位置算出部35及び粒子生成部36へ出力する。
その後、制御部31は、ステップ#nにおける各着目粒子の速度vの入力を速度算出部32から受ける。また、制御部31は、ステップ#nにおける各着目粒子の密度ρの入力を密度算出部33から受ける。また、制御部31は、ステップ#nにおける各着目粒子の圧力Pの入力を圧力算出部34から受ける。さらに、制御部2は、ステップ#nにおける各着目粒子の位置rの入力を位置算出部35から受ける。ここで、1ステップとは、ΔtSPH時間経過した後の状態である。例えば、ステップ#n+1における速度とは、ステップ#nにおける速度からΔtSPH時間後の速度である。
また、制御部31は、図2の領域100における境界から一定の距離Δxの範囲を境界領域として記憶している。ここで、Δxは粒子の大きさから設定されることが好ましい。具体的には、粒子の大きさ程度の幅を持たせることが好ましい。例えば、粒子の直径が1mの場合、1mより少し大きい幅をΔxとすることが好ましい。そして、制御部31は、受信したステップ#nにおける着目粒子の位置rから、図2の領域100における境界の外に着目粒子が出たか否かを判定する。そして、制御部31は、着目粒子が境界の外に出たと判定した場合、ステップ#n+1以降の計算においてその着目粒子を除外する。図8は、粒子の流出を模式的に表した図である。図8に示すように、粒子202が境界の外に出て粒子203となった場合、制御部31は、粒子203を計算から除外する。
さらに、制御部31は、ステップ#nにおいて生成された粒子の位置rの入力を粒子生成部36から受ける。そして、制御部31は、ステップ#n+1以降の計算において受信した新たに生成された粒子を組み入れる。
また、制御部31は、ステップ#nにおける境界領域内に存在する粒子の速度vを、2次元伝播計算処理部2から受信した流量M及び流量Nから求め、速度算出部32及び位置算出部35から受信したステップ#nにおける速度vと置き換える。また、制御部31は、ステップ#nにおける境界領域内に存在する粒子の密度ρ及び圧力Pは、海水の一般的な密度及び圧力に置き換える。このようにして、制御部31は、ステップ#nにおける境界条件を都度求め、境界領域内に存在する粒子に割り当てる。これは、本実施例での3次元粒子法計算においては境界外には粒子が存在しないため、境界領域に存在する粒子の動きは正確に求めることができない。そこで、2次元伝播計算によって求められた正確性の高い情報を基に境界領域に存在する粒子の動きを推定することで、計算の正確性を向上させるためである。
ただし、上述したように2次元伝播計算と3次元粒子法計算ではステップの精度が異なる。すなわち、3次元粒子法計算の各ステップ#nに一致する、2次元伝播計算のステップkがあるとは限らない。そこで、制御部31は、3次元粒子法計算の各ステップ#nに一致するタイミングの水位η、流量M及び流量Nを算出する。ここで、図9は、境界条件算出における補間を説明するための図である。図9では、例えば、2次元伝播計算のステップkにおいて3次元粒子法計算が開始され、2次元伝播計算の次のステップであるステップk+1と3次元粒子法計算のステップ#nのタイミングが一致する場合を表している。この場合、2次元伝播計算のステップ間の時間はΔtであり、3次元粒子法計算のステップ間の時間はΔtSPHである。この場合、ステップ#1では、対応401のように、制御部31は、ステップkにおける水位η、流量M及び流量Nを用いて、境界条件を算出することができる。また、ステップ#nでは、制御部31は、対応402のように、ステップk+1における水位η、流量M及び流量Nを用いて、境界条件を算出することができる。しかし、ステップ#2〜ステップ#n−1に対応する、2次元伝播計算のステップは存在しない。そこで、例えばステップ#2であれば、制御部31は、ステップk及びk+1における水位η、流量M及び水位Nを用いて、内挿を行いステップ#2に対応する内挿ステップ403における水位η、流量M及び水位Nを求める。そして、制御部31は、対応404のように、求めた内挿ステップ403における水位η、流量M及び流量Nを用いて、ステップ#2の境界条件を算出する。
例えば、内挿ステップにおける水位は、次の式(7)のように求めることができる。ここでは、粒子法計算において境界条件を得たいステップの時刻をtSPHとしている。そして、tSPHに最も近い2次元伝播計算のステップであるステップn及びステップn+1の時刻をnΔt及び(n+1)Δtとしている。
Figure 0006048062
そして、制御部31は、ステップ#nにおける各着目粒子の密度ρ、圧力P及び位置rを、ステップ#n+1の計算に用いる各着目粒子の現ステップの密度、圧力及び位置として、速度算出部32へ出力する。
また、制御部31は、ステップ#nにおける各着目粒子の密度ρをステップ#n+1の計算に用いる現ステップの密度として密度算出部33へ出力する。
また、密度算出部33は、ステップ#n+1における各着目粒子の密度ρをステップ#n+1の圧力計算に用いるために圧力算出部34へ出力する。
また、制御部31は、ステップ#nにおける各着目粒子の速度vを、ステップ#n+1の計算に用いる現ステップの速度として位置算出部35へ出力する。
以上の動作をまとめると、制御部31は、ステップ#nからΔtSPH時間経過した後のステップ#n+1における、各着目粒子の速度v、密度ρ、圧力P及び位置rの入力を受ける。次に、制御部31は、粒子の位置から消滅させる粒子を決定し、さらに粒子生成部36からの情報を用いて粒子を生成する。そして、制御部31は、粒子の消滅及び生成を行った後、2次元伝播計算により算出された水位η、流量M及び流量Nから境界領域にある粒子の速度v、密度ρ、圧力P及び位置rを求めて境界条件を設定する。その後、制御部31は、それらの情報を各部に出力し、ステップ#nからΔtSPH時間経過した後のステップ#n+1における、各着目粒子の速度v、密度ρ、圧力P及び位置rの入力を受ける。このように、制御部31は、既にステップ#nにおける速度v、密度ρ、圧力P及び位置rを受信している場合には、次にステップ#n+1における各着目粒子の速度v、密度ρ、圧力P及び位置rを受信する。そして、制御部31は、指定されたシミュレーションの終了のタイミングまで、各ステップにおける各着目粒子の速度v、密度ρ、圧力P及び位置rを順次受信する。
さらに、制御部31は、ステップ#nにおける各着目粒子の最新の速度v、密度ρ、圧力P及び位置rの各情報の入力を受けると、受信したステップ#nにおける各情報を結果記憶部4へ記憶させる。
そして、制御部31は、受信したステップ#nにおける各着目粒子の位置を用いてステップ#nにおける各着目粒子の近傍粒子を求め、近傍粒子リストを作成する。そして、制御部31は、作成した近傍粒子リストを速度算出部32、密度算出部33、圧力算出部34、位置算出部35及び粒子生成部36へ出力する。このように、制御部31は、ステップ毎の速度v、密度ρ、圧力P及び位置r、並びに近傍粒子リストを各部へ出力する。そして、制御部31は、現在のデータを用いて1ステップ後のデータを算出するよう、速度算出部32、密度算出部33、圧力算出部34、位置算出部35及び粒子生成部36に対して指示し、最終ステップまで繰り返させる。
速度算出部32は、ステップ#n+1の計算に用いる各着目粒子の現ステップの密度、圧力及び位置として、ステップ#nにおける各着目粒子の密度ρ、圧力P及び位置rの入力を制御部31から受ける。
速度算出部32は、次の式(8)で表される運動方程式を記憶している。
Figure 0006048062
ここで、式(8)における、変数及び係数の添え字は各粒子を表す。例えば、r,v,ρ,P,mがあった場合、それらは、粒子aの位置ベクトル、速度ベクトル、密度、圧力、質量を表している。Wは、粒子の分布から連続場を構成するために用いるカーネル関数である。カーネル関数としては、例えば、次の式(9)に示すような5次のスプライン関数などが用いられる。ただし、ここでq=r/hSPHである。
Figure 0006048062
また、式(8)におけるΠabは粒子a,b間に働く人工粘性項である。例えば、次の式(10)に示されるような形式が採用される。
Figure 0006048062
ここで、式(10)におけるμは粘性係数であり、次の式(11)で表される。ただし、ここでηSPHは分母が0になることを避けるためのパラメータであり、考えている空間スケールに対して十分小さい値を設定する。例えば、ηSPH=0.01hSPH等でよい。
Figure 0006048062
そして、速度算出部32は、式(8)に受信したステップ#nにおける着目粒子a及び近傍粒子bの密度ρ、圧力P及び位置rを代入し、ステップ#nにおける加速度dv/dtを求める。
さらに、速度算出部32は、境界からの反発力を次の式(12)を用いて求める。
Figure 0006048062
ここで、図10を参照して、式(12)について説明する。図10は、境界からの反発力を説明するための図である。ここでは、図10における粒子500に加えられる境界からの反発力について説明する。
図10に示すように、dは境界面501からの所定の距離である。そして、境界層502、境界面からdの幅の層を示している。また、sは粒子500の境界層502への侵入距離である。
本実施例では、境界面501から幅dの境界層502を考え、境界層502内に粒子500が突入した場合に、粒子500に対して境界面501から粒子500の向きに反発加速度が加わるものとする。反発力fは、弾性ばねによるものと同じ形式でf=ksとして与える。ばね剛性係数kは、音速cの粒子が有する運動エネルギーmc /2と、粒子が境界層502へdの長さだけ進入したときの弾性歪エネルギーkd/2とが同程度となるように定める。すなわち、連続体の運動の解析を行なうシミュレーションにおいては、粒子500は音速を超えることはないと考えられるので、音速で粒子500が進入した場合に境界面501においてエネルギーがなくなるように設定しておけば、粒子500は境界面を突き抜けることはなくなる。このため、βを1程度のパラメータとする。このようにβを設定することで、式(12)が導かれる。
式(12)で表される反発力を用いることで、境界壁からの粒子のすり抜けを防止し、導入する摩擦やfree slip条件の正確性を向上でき、さらに、境界で物理的に無意味なエネルギー消費の発生を抑えることができる。また、式(12)は簡単な式で表されているので、式(12)で表される反発力を用いた場合には、計算負荷を抑えることができる。
ここで、本実施例では式(12)を用いたが、この反発力は他の式を用いてもよい。例えば、反発力が過大になる可能性を考慮して、粒子が境界面に近づいている場合には、弾性ばねを考慮し、逆に粒子が境界面から遠ざかる場合には、弾性ばねを考慮しないように反発力を設定してもよい。このような反発力は一方向弾性ばね法と呼ばれることがある。エネルギー保存の観点からは粒子と境界面の相対速度に関係なく弾性ばねを考慮する式(12)のような反発力が好ましいが、境界での不自然な反発を避けるために一方向弾性ばね法をもちいてもよい。
そして、速度算出部32は、運動方程式である式(8)を用いて算出した加速度dv/dtと式(12)を用いて求めた反発力とを加えてステップ#nにおける加速度v´を算出する。ここでは、加速度v´は、加速度ベクトルを表す。
さらに、速度算出部32は、次の時間積分の式(13)を記憶している。ここで、vはステップ#nにおける速度を表し、v´はステップ#nにおける加速度を表している。
Figure 0006048062
そして、速度算出部32は、求めたステップ#nの速度v及び加速度v´を式(13)に代入して、ステップ#nの速度vをステップ#nにおける加速度v´で時間積分する。これにより、速度算出部32は、ステップ#n+1の速度vを求める。
速度算出部32は、算出したステップ#n+1の各粒子の速度vを制御部31へ出力する。また、速度算出部32は、算出したステップ#n+1の各粒子の速度vを密度算出部33へ出力する。さらに、速度算出部32は、算出したステップ#n+1の各粒子の速度vを位置算出部35へ出力する。
密度算出部33は、ステップ#n+1の計算に用いる現ステップの密度として、ステップ#nにおける各着目粒子の密度ρの入力を制御部31から受ける。さらに、密度算出部33は、ステップ#n+1の各粒子の速度vの入力を速度算出部32から受ける。
密度算出部33は、次の式(14)で表される連続方程式を記憶している。
Figure 0006048062
密度算出部33は、速度算出部32から入力されたステップ#n+1における各粒子の速度vを式(14)に代入し、ステップ#n+1における各粒子の密度時間微分dρ/dtを算出する。
さらに、密度算出部33は、次の式(15)で表される時間積分の式を記憶している。
Figure 0006048062
密度算出部33は、算出したステップ#nにおける各粒子の密度及び各粒子の密度時間微分dρ/dtの値を式(15)に代入して、ステップ#nの密度ρをステップ#nにおける密度時間微分dρ/dtで時間積分する。これにより、密度算出部33は、ステップ#n+1の密度ρを求める。
密度算出部33は、算出したステップ#n+1の密度ρを、制御部31及び圧力算出部34へ出力する。
圧力算出部34は、ステップ#n+1の密度ρの入力を密度算出部33から受ける。
圧力算出部34は、次の式(16)で表される状態方程式を記憶している。ここで、c及びρは、それぞれ流体の音速と基準密度である。
Figure 0006048062
圧力算出部34は、密度算出部33から入力されたステップ#n+1の密度ρを式(16)に代入して、ステップ#n+1における圧力Pを算出する。
圧力算出部34は、算出したステップ#n+1における圧力Pを制御部31へ出力する。
位置算出部35は、ステップ#n+1における速度vの入力を速度算出部32から受ける。
位置算出部35は、次の式(17)で表される時間積分の式を記憶している。
Figure 0006048062
位置算出部35は、速度算出部32から入力されたステップ#n+1における速度vを式(16)に代入して、ステップ#nの位置rをステップ#n+1における速度vで時間積分する。これにより、位置算出部35は、ステップ#n+1の位置rを求める。
粒子生成部36は、境界面を1つの粒子に相当する所定の領域で分割した分割領域を記憶している。ここで、所定の領域は、密に配置されていればどのような形状でも良く、例えば、矩形でも円形でもよい。以下では、この所定の領域を「粒子生成子」と呼ぶ。粒子生成子は、境界において流入があった時に粒子を生成することで、粒子法計算を破たんなく実行させるために設定されるものである。粒子生成子は、領域100の境界面を微小部分に分割した面要素として、例えば、中心座標、面要素の法線ベクトル、面要素の面積、流入出体積、基準体積を持つ。ここで、基準体積とは、粒子一つに対応する体積である。そして、以下で説明する粒子生成子における流出入体積は、計算開始の時点では0に設定され、流入があれば増加し、流出があれば減少する。粒子生成子を通じて海水が流出している場合には流出入体積は負の値を取る。
粒子生成部36は、ステップ#nにおける各粒子生成子の位置における2次元伝播計算により求められた速度の入力を制御部31から受ける。
そして、粒子生成部36は、受信した速度から各粒子生成子の位置における海水の流速を表す流速ベクトルを求める。さらに、粒子生成部36は、粒子生成子が形成する所定の領域の法線方向を向いた大きさが面積に比例する面積ベクトルを生成する。粒子生成子の面積をSとし、法線ベクトルをuとした場合、面積ベクトルは、例えば、uSとして表すことができる。
そして、粒子生成部36は、各粒子生成子における流速ベクトルと面積ベクトルとの内積を求め、ステップ#nにおける各粒子生成子における海水の流出入体積の増減分を算出する。
そして、粒子生成部36は、各粒子生成子におけるステップ#nの海水の流出入体積に流出入体積の増減分を加算して、各粒子生成子における新しいステップ#nの海水の流出入体積を求める。例えば、古い流出入体積をVoldとして、新しい流出入体積をVnewとした場合、粒子生成部36は、次の式(18)により新しい海水の流出入体積が算出できる。
Figure 0006048062
次に、粒子生成部36は、各粒子生成子におけるステップ#n+1の海水の流出入体積が1つの粒子の体積と等しい基準体積以上か否かを判定する。海水の流出入体積が基準体積以上となっている粒子生成子があれば、粒子生成部36は、その粒子生成子に対応する決められた位置への粒子の生成を決定する。ここで、粒子生成子に対応する決められた位置とは、例えば、粒子生成子に接するように粒子を配置する位置でも良いし、粒子生成子から法線方向に所定の距離離れた位置に粒子が配置される位置でも良い。そして、粒子生成部36は、粒子の生成を決定した粒子生成子におけるステップ#n+1の海水の流出入体積から基準体積を減算して、ステップ#n+1の海水の流出入体積とする。
図11は、粒子の流入を模式的に表した図である。境界の外部にある粒子204の体積を有する海水が境界を越えて内部に流入した場合に、粒子生成部36は、流入した粒子204に対応する境界の内部の位置に新しく粒子205を生成することを決定する。
また、ステップ#n+1の海水の流出入体積がマイナス方向の基準体積を下回った粒子生成子がある場合、粒子生成部36は、その粒子生成子におけるステップ#n+1の海水の流出入体積に基準体積を加算して、ステップ#n+1の海水の流出入体積とする。
そして、粒子生成部36は、生成を決定した粒子の位置を制御部31に通知する。生成した粒子の位置とは、例えば、生成した粒子の中心の位置などでよい。
粒子生成部36は、ステップが終了するまで、上記の粒子生成の処理を繰り返す。
次に、図12を参照して、本実施例に係る2次元伝播計算の処理について説明する。図12は、実施例2に係る2次元伝播計算の処理のフローチャートである。
制御部21は、ユーザインタフェース1からの入力により、領域101の分割した各格子における、水位、x方向の速度及びy方向の速度の初期値である水位η、流量M及び流量Nを取得する(ステップS101)。また、制御部21は、各格子の中から解析対象とする格子を選択格子として一つ選択する。そして、制御部21は、選択格子における、水位η、流量M及び流量Nを水位算出部22及び流量算出部23へ出力する。
水位算出部22は、選択格子における水位の時間微分∂η/∂tを式(1)の連続方程式を用いて算出する(ステップS102)。
また、流量算出部23は、選択格子における流量M及び流量Nの時間微分である∂M/∂t及び∂N/∂tを式(3)及び式(4)の運動方程式を用いて算出する(ステップS103)。
制御部21は、全ての格子について計算が完了したか否かを判定する(ステップS104)。全ての格子について計算が完了していない場合(ステップS104:否定)、制御部21は、新たな選択格子を選択し、ステップS102へ戻る。
これに対して、全ての格子について計算が完了している場合(ステップS104:肯定)、制御部21は、選択格子を一つ選択し、水位算出部22及び流量算出部23に次のステップの水位及び速度の算出を指示する。
水位算出部22は、選択格子の水位ηを選択格子の水位の時間微分∂η/∂tを基に時間積分し(ステップS105)、次のステップの水位を求める。そして、水位算出部22は、算出した次のステップの水位ηを制御部21へ出力する。
また、流量算出部23は、選択格子の流量M及び流量Nを選択格子の速度の時間微分∂M/∂t及び∂N/∂tを基に時間積分し(ステップS106)、次のステップの流量M及び流量Nを求める。そして、流量算出部23は、算出した次のステップの流量M及び流量Nを制御部21へ出力する。
制御部21は、全ての格子について計算が完了したか否かを判定する(ステップS107)。全ての格子について計算が完了していない場合(ステップS107:否定)、制御部21は、新たな選択格子を選択し、ステップS105へ戻る。
これに対して、全ての格子について計算が完了している場合(ステップS107:肯定)、制御部21は、解析段階を1ステップ進める(ステップS108)。
そして、制御部21は、現在の解析段階における計算結果を3次元SPH計算処理部3へ出力する(ステップS109)。ここで、このフローでは、制御部21は、ステップ毎に計算が終わった時点で計算結果を出力しているが、これに限らず、例えば、全てのステップが完了してから全てのステップの計算結果をまとめて出力しても良い。
制御部21は、最終状態の条件を満たしたか否かにより、全ステップが完了したか否かを判定する(ステップS110)。全ステップが完了していない場合(ステップS110:否定)、制御部21は、各格子の中から解析対象とする格子を選択格子として一つ選択する。そして、制御部21は、選択格子における、現ステップの水位η、流量M及び流量Nを水位算出部22及び流量算出部23へ出力し、次のステップの計算を指示し、ステップS102へ戻る。
これに対して、全ステップが完了している場合(ステップS110:肯定)、制御部21は2次元伝播計算の処理を終了する。
次に、図13を参照して、本実施例に係る3次元SPH計算の処理について説明する。図13は、実施例2に係る3次元SPH計算の処理のフローチャートである。
制御部31は、ステップ毎の各格子における水位η、x方向の流量M及びy方向の流量Nのデータを2次元伝播計算処理部2から取得する(ステップS201)。そして、制御部31は、3次元粒子法計算の開始条件を満たしたステップを特定し、特定したステップの水位η、x方向の流量M及びy方向の流量Nを取得する。さらに、制御部31は、内挿により、3次元粒子法計算に用いる位置座標を内挿により求める。そして、制御部31は、求めた位置座標の全水深に応じて粒子を配置し、各粒子の速度、密度、圧力及び位置の初期値を設定する。
また、制御部31は、粒子分布に基づいて近傍粒子リストを作成する(ステップS202)。そして、制御部31は、一つ着目粒子を選択し、着目粒子の密度ρ、圧力P及び位置rを速度算出部32へ出力する。また、制御部31は、着目粒子の密度ρを密度算出部33へ出力する。また、制御部31は、着目粒子の速度vを位置算出部35へ出力する。
速度算出部32は、着目粒子の密度ρ、圧力P及び位置rの入力を制御部31から受ける。そして、速度算出部32は、式(8)の運動方程式に着目粒子の密度ρ、圧力P及び位置rを代入し着目粒子の加速度dv/dtを算出する(ステップS203)。
さらに、速度算出部32は、着目粒子に対する境界からの反発力を式(12)を用いて算出する。着目粒子が境界近傍に存在しない場合は、反発力は0になる。そして、速度算出部32は、算出した反発力を加速度dv/dtに加算して(ステップS204)、加速度v´を算出する。
速度算出部32は、式(13)を用いて、現ステップにおける着目粒子の速度vを現ステップにおける加速度v´で時間積分して、次のステップの速度vを算出する(ステップS205)。そして、速度算出部32は、算出した次のステップの着目粒子の速度vを制御部31へ出力する。
制御部31は、着目粒子が境界領域に存在する場合、2次元伝播計算結果から着目粒子の速度vを算出し、速度算出部32から入力された速度vに変えて、算出した速度vを着目粒子の速度として設定する(ステップS206)。
制御部31は、領域100内に存在する全ての粒子について速度の計算が完了したか否かを判定する(ステップS207)。全ての粒子について計算が完了していない場合(ステップS207:否定)、制御部31は、計算を行っていない粒子から着目粒子を選択し、ステップS203へ戻る。
全ての粒子について計算が完了している場合(ステップS207:肯定)、制御部31は、着目粒子を1つ選択して、密度算出部33、圧力算出部34及び位置算出部35に対して計算実行を指示し、ステップS208へ進む。
密度算出部33は、速度算出部32から受信した次のステップの着目粒子の速度vを着目粒子の密度時間微分を式(14)の連続方程式に代入して密度時間微分dρ/dtを算出する(ステップS208)。
次に、密度算出部33は、制御部31から受信した現ステップの着目粒子の密度ρ及び算出した密度時間微分dρ/dtの値を式(15)に代入する。そして、密度算出部33は、現ステップにおける着目粒子の密度ρを密度時間微分dρ/dtで積分して、次のステップの密度ρを算出する(ステップS209)。密度算出部33は、算出した次のステップの密度ρを制御部31へ出力する。
圧力算出部34は、密度算出部33から受信した次のステップの密度ρを、式(16)の状態方程式に代入して、次のステップの圧力Pを算出する(ステップS210)。圧力算出部34は、算出した次のステップの圧力Pを制御部31へ出力する。
位置算出部35は、式(17)を用いて、現ステップの位置rを速度算出部32から受信した次のステップにおける速度vで時間積分し次のステップの位置rを求める(ステップS211)。位置算出部35は、算出した次のステップの位置rを制御部31へ出力する。
制御部31は、次のステップの位置rを用いて領域境界を越えて流出した粒子を特定し、特定した粒子を領域100の中から消去する(ステップS212)。
制御部31は、領域100内に存在する全ての粒子について密度、圧力及び位置の計算が完了したか否かを判定する(ステップS213)。全ての粒子について計算が完了していない場合(ステップS213:否定)、制御部31は、計算を行っていない粒子から着目粒子を選択し、ステップS208へ戻る。
これに対して、全ての粒子について計算が完了している場合(ステップS213:肯定)、制御部31は、粒子生成を粒子生成部36に指示する。粒子生成部36は、各粒子生成子の位置における2次元伝播計算により求められた速度から、現在の海水の流出入体積を求め粒子の生成処理を行う(ステップS214)。この粒子の生成処理の流れについては次に詳細に説明する。そして、粒子生成部36は、生成を決定した粒子の位置を制御部31へ通知する。
制御部31は、解析段階を1ステップ進める(ステップS215)。
そして、制御部21は、現在の解析段階における計算結果を結果記憶部4へ出力する(ステップS216)。ここで、このフローチャートでは、制御部31は、ステップ毎に計算が終わった時点で計算結果を出力しているが、これに限らず、例えば、全てのステップが完了してから全てのステップの計算結果をまとめて出力しても良い。
制御部31は、最終状態の条件を満たしたか否かにより、全ステップが完了したか否かを判定する(ステップS217)。全ステップが完了していない場合(ステップS217:否定)、制御部31は、ステップS202へ戻る。
これに対して、全ステップが完了している場合(ステップS217:肯定)、制御部31は3次元SPH計算の処理を終了する。
次に、図14を参照して、本実施例に係る3次元SPH計算の処理について説明する。図14は、実施例2に係る3次元SPH計算の処理のフローチャートである。図14のフローチャートは、図13のフローチャートにおけるステップS215にあたる。
粒子生成部36は、粒子生成子を一つ選択する。そして、粒子生成部36は、選択した粒子生成子における現ステップの海水の流速を計算する(ステップS301)。
次に、粒子生成部36は、海水の流速と粒子生成子の面積から、選択した粒子生成子における海水の流出入体積の増減を算出する(ステップS302)。
次に、粒子生成部36は、現ステップの流出入体積Vに増減分を加算し、現ステップの新しい流出入体積Vを算出する(ステップS303)。
次に、粒子生成部36は、新しい流出入体積Vが基準体積V以上か否かを判定する(ステップS304)。基準体積以上の場合(ステップS304:肯定)、粒子生成部36は、選択している粒子生成子に対応する位置に粒子の生成を決定する(ステップS305)。
次に、粒子生成部36は、現ステップの新しい流出入体積Vから基準体積Vを減算し、次のステップの流出入体積Vを算出する(ステップS306)。
これに対して、新しい流出入体積Vが基準体積V以下の場合(ステップS304:否定)、粒子生成部36は、新しい流出入体積Vが基準体積Vのマイナス値、すなわち−Vより小さいか否かを判定する(ステップS307)。−Vより小さい場合(ステップS307:肯定)、粒子生成部36は、現ステップの新しい流出入体積Vに基準体積Vを加算し、次のステップの流出入体積Vを算出する(ステップS308)。これに対して、新しい流出入体積Vが−V以上の場合(ステップS307:否定)、粒子生成部36は、ステップS309へ進む。
そして、粒子生成部36は、全粒子生成子について粒子の生成処理が終了したか否かを判定する(ステップS309)。粒子の生成処理が終了していない粒子生成子がある場合(ステップS309:否定)、粒子生成部36は、ステップS301へ戻る。
これに対して、全粒子生成子について粒子の生成処理が終了した場合(ステップS309:肯定)、粒子生成部36は、粒子の生成処理を終了する。
以上に説明したように、本実施例に係るシミュレーション装置は、連続体の広い移動に対するシミュレーションにおいて、広い範囲では2次元伝播計算を行い、港湾部や都市部といった3次元の形状の影響が大きい領域では3次元粒子方計算を行う。そして、本実施例に係るシミュレーション装置は、上述した各式などを用いて、2次元伝播計算及び3次元粒子法計算を行う。これにより、津波のように大量の流体の移動が伴う事象においても、破綻無く計算を行うことができ、計算量を抑えながら、精度の良いシミュレーションを行うことができる。
また、本実施例に係るシミュレーション装置は、バッファ領域を用いた従来技術に比べて、境界領域の幅を粒子の影響半径と同程度に抑えることができ、計算量を確実に抑えることができる。
さらに、粒子生成子を3次元空間内で向きを持つものとして設定できるので、2次元伝播計算と3次元粒子法計算の境界の形状を任意に設定でき、シミュレーションの条件の自由度を向上させることができる。
次に、実施例3について説明する。本実施例に係るシミュレーション装置は、2次元伝播計算を最後まで終わらせずに途中の状態で、並行して3次元粒子法計算を行うことが実施例1及び実施例2と異なるものである。以下では、同じ機能を有する各部については説明を省略する。
2次元伝播計算処理部2は、ステップ毎にそのステップにおける計算結果を3次元SPH計算処理部3へ出力する。
3次元SPH計算処理部3は、計算結果を受信するとそれまでに受信した計算結果から3次元粒子法計算を開始する条件が満たされたか否かを判定する。
3次元粒子法計算を開始する条件が満たされている場合、3次元SPH計算処理部3は、受信した計算結果の水位及び流量から各粒子の位置及び速度の初期条件を生成する。さらに、3次元SPH計算処理部3は、各粒子の密度及び圧力の初期条件は海水の一般的な密度圧力を用いる。
そして、3次元SPH計算処理部3は、ステップ#nの速度、密度、圧力及び位置を用いて3次元粒子法計算を行い、ステップ#n+1の速度、密度、圧力及び位置を算出する。
3次元SPH計算処理部3は、ステップ#n+1の速度、密度、圧力及び位置の算出後、ステップ#n+2における境界条件、すなわち、境界領域に存在する粒子の速度が求められるか否かを判定する。
境界条件が求められるのであれば、3次元SPH計算処理部3は、ステップ#n+1の速度、密度、圧力及び位置及び境界条件を用いて、ステップ#n+2の速度、密度、圧力及び位置の算出を行う。
これに対して、境界条件が求められない場合、3次元SPH計算処理部3は、2次元伝播計算における次のステップの計算結果を2次元伝播計算処理部2から取得する。その後、3次元SPH計算処理部3は、取得した2次元伝播計算における次のステップの計算結果から境界条件を求める。そして、3次元SPH計算処理部3は、ステップ#n+1の速度、密度、圧力及び位置及び境界条件を用いて、ステップ#n+2の速度、密度、圧力及び位置の算出を行う。
以下では、図15を参照して、2次元伝播計算及び3次元粒子法計算の処理の全体的な流れについて説明する。図15は、実施例3に係るシミュレーション装置による解析処理のフローチャートである。
2次元伝播計算処理部2は、2次元伝播計算を1ステップ進める(ステップS401)。そして、2次元伝播計算処理部2は、現在のステップの計算結果を3次元SPH計算処理部3へ出力する。
3次元SPH計算処理部3は、現在までに受信した2次元伝播計算の計算結果から、3次元粒子法計算の開始条件が満たされたか否かを判定する(ステップS402)。例えば、3次元SPH計算処理部3は、領域100の境界における水位が50cmを超えることを、3次元粒子法計算を開始する条件として記憶している。3次元粒子法計算の開始条件が満たされていない場合(ステップS402:否定)、2次元伝播計算処理部2は、ステップS401に戻る。
これに対して、3次元粒子法計算の開始条件が満たされた場合(ステップS402:肯定)、3次元SPH計算処理部3は、2次元伝播計算による水位及び流量から3次元粒子法計算における位置及び速度の初期条件を生成する(ステップS403)。
そして、3次元SPH計算処理部3は、1ステップ後の3次元粒子法計算を実施し、各粒子の速度、密度、圧力及び位置を算出する(ステップS404)。
その後、3次元SPH計算処理部3は、現ステップまでの計算結果から、3次元粒子法計算の終了条件を満たすか否かを判定する(ステップS405)。
3次元粒子法計算の終了条件を満たしていない場合(ステップS405:否定)、3次元SPH計算処理部3は、次のステップの境界条件を設定するための2次元伝播計算による水位及び流量のデータが存在するか否かを判定する(ステップS406)。
次のステップの境界条件を設定するための2次元伝播計算による水位及び流量のデータが存在しない場合(ステップS406:否定)、2次元伝播計算処理部2は、2次元伝播計算を1ステップ進める(ステップS407)。そして、2次元伝播計算処理部2は、現在のステップの計算結果を3次元SPH計算処理部3へ出力する。
これに対して、次のステップの境界条件を設定するための水位及び流量のデータが存在する場合(ステップS406:肯定)、3次元SPH計算処理部3は、ステップS408へ進む。
そして、3次元SPH計算処理部3は、2次元伝播計算による水位及び流量から3次元粒子法計算における境界条件を生成する(ステップS408)。その後、3次元SPH計算処理部3は、ステップS404へ戻る。
一方、3次元粒子法計算の終了条件を満たされた場合(ステップS405:肯定)、3次元SPH計算処理部3は、3次元粒子法計算を終了する。
ここで、図15では、3次元粒子法計算において2次元伝播計算の計算結果が要求されてから、2次元伝播計算を行うように説明したが、2次元伝播計算処理部2は、3次元粒子法計算の進行状況に関らず2次元伝播計算を進めてよい。
以上に説明したように、本実施例に係るシミュレーション装置は、2次元伝播計算の終了を待たずに3次元粒子法計算を開始し、それらを並行して実行する。この場合、例えば、2次元伝播計算を行う計算機と3次元粒子法計算を行う計算機を分けることができる。そして、3次元粒子法計算を行う計算機は、2次元伝播計算を行う計算機が出力した2次元伝播計算の結果のファイルを用いて、3次元粒子法計算を実行できる。これにより、全体の計算時間を減らすことができる。
また、上記実施例で説明した各種の処理は、あらかじめ用意されたプログラムをコンピュータで実行することによって実現することができる。そこで、以下では、図16を用いて、図1に示したシミュレーション装置と同様の機能を有するシミュレーションプログラムを実行するコンピュータの一例を説明する。
図16は、シミュレーションプログラムを実行するコンピュータを示す図である。図16に示すように、コンピュータ1000は、キャッシュ1010と、CPU(Central Processing Unit)1020と、HDD(Hard Disk Drive)1030と、RAM(Random Access Memory)1040及びROM(Read Only Memory)1050を有する。キャッシュ1010、CPU1020、HDD1030、RAM1040及びROM1050は、それぞれバスによって接続されている。
HDD1030には、図1に示したシミュレーション装置と同様の機能を発揮する各種のシミュレーションプログラム1031が予め記憶されている。
そして、CPU1020は、これらのシミュレーションプログラム1031を読み出して実行する。これにより、図16に示すように、シミュレーションプログラム1031は、シミュレーションプロセス1021になる。
なお、上記したシミュレーションプログラム1031については、必ずしもHDD1030に記憶させなくてもよい。例えば、コンピュータ1000に挿入されるフレキシブルディスク、CD(Compact Disk)−ROM、DVD(Digital Versatile Disc)、光磁気ディスク、IC(Integrated Circuit)カードなどの「可搬用の物理媒体」にシミュレーションプログラム1031を記憶させてもよい。または、コンピュータ1000の内外に備えられるハードディスクドライブ(HDD)などの「固定用の物理媒体」にシミュレーションプログラム1031を記憶させてもよい。または、公衆回線、インターネット、LAN(Local Area Network)、WAN(Wide Area Network)などを介してコンピュータ1000に接続される「他のコンピュータ(またはサーバ)」にシミュレーションプログラム1031を記憶させてもよい。そして、コンピュータ1000は、上述したフレキシブルディスク等から各プログラムを読み出して実行するようにしてもよい。
1 ユーザインタフェース
2 2次元伝播計算処理部
3 3次元SPH計算処理部
4 結果記憶部
21 制御部
22 水位算出部
23 流量算出部
31 制御部
32 速度算出部
33 密度算出部
34 圧力算出部
35 位置算出部
36 粒子生成部

Claims (7)

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

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2012231224A JP6048062B2 (ja) 2012-10-18 2012-10-18 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
US13/967,468 US9557203B2 (en) 2012-10-18 2013-08-15 Simulation method and simulation apparatus

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2012231224A JP6048062B2 (ja) 2012-10-18 2012-10-18 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置

Publications (2)

Publication Number Publication Date
JP2014081900A JP2014081900A (ja) 2014-05-08
JP6048062B2 true JP6048062B2 (ja) 2016-12-21

Family

ID=50486097

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2012231224A Active JP6048062B2 (ja) 2012-10-18 2012-10-18 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置

Country Status (2)

Country Link
US (1) US9557203B2 (ja)
JP (1) JP6048062B2 (ja)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6521777B2 (ja) * 2014-10-22 2019-05-29 三菱電機株式会社 津波監視システム
KR102459848B1 (ko) 2014-10-24 2022-10-27 삼성전자주식회사 입자에 기반하여 대상 객체를 고속으로 모델링하는 방법 및 장치
JP6458501B2 (ja) 2015-01-06 2019-01-30 富士通株式会社 シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置
JP6437883B2 (ja) * 2015-06-02 2018-12-12 公益財団法人鉄道総合技術研究所 津波伝播特性を利用した沿岸の早期津波予測方法
KR101718755B1 (ko) * 2015-09-23 2017-03-23 동남이엔씨(주) 홍수발생시 도로 및 하천으로 유입되는 빗물의 흐름을 모의하기 위한 3차원 유체 시뮬레이션 시스템 및 이를 이용한 시뮬레이션 방법
KR101820946B1 (ko) * 2015-11-19 2018-03-09 이에이트 주식회사 파티클 기반의 sph에 의한 강우 시뮬레이션 방법 및 이를 실행하도록 구성되어진 강우 시뮬레이션 플랫폼
JP2018018164A (ja) 2016-07-25 2018-02-01 富士通株式会社 粒子シミュレーションプログラム、粒子シミュレーション方法、及び情報処理装置
JP6897477B2 (ja) * 2017-10-10 2021-06-30 富士通株式会社 流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置
CN108984900B (zh) * 2018-07-12 2023-05-26 交通运输部天津水运工程科学研究所 一种基于深水区风浪条件的港内波高变化分析方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3993826B2 (ja) 2003-01-15 2007-10-17 大成建設株式会社 振動予測方法及び振動予測システム
JP2005115701A (ja) 2003-10-08 2005-04-28 Port & Airport Research Institute 地形と構造物の電子データを用いた数値解析条件の設定方法
US7349832B2 (en) * 2004-02-17 2008-03-25 Pixar Water particle manipulation
JP5007391B2 (ja) 2006-09-29 2012-08-22 国立大学法人京都大学 津波波源推定方法及び津波波高予測方法並びにその関連技術
TW200945249A (en) * 2008-04-28 2009-11-01 Inst Information Industry Method for rendering fluid
JP5433913B2 (ja) 2008-08-29 2014-03-05 国立大学法人京都大学 波浪予測システム
US8878856B1 (en) * 2010-06-29 2014-11-04 Nvidia Corporation System, method, and computer program product for depicting a body of water utilizing a height field and particles
JP5576309B2 (ja) * 2011-01-31 2014-08-20 インターナショナル・ビジネス・マシーンズ・コーポレーション 粒子法における自由表面に位置する粒子の正確な判定

Also Published As

Publication number Publication date
US9557203B2 (en) 2017-01-31
JP2014081900A (ja) 2014-05-08
US20140114588A1 (en) 2014-04-24

Similar Documents

Publication Publication Date Title
JP6048062B2 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
Issakhov et al. Numerical simulation of the movement of water surface of dam break flow by VOF methods for various obstacles
Cremonesi et al. A Lagrangian finite element approach for the simulation of water-waves induced by landslides
JP6061499B2 (ja) 壊食予測方法および壊食予測システム、この予測に用いる壊食特性データベースおよびその構築方法
US9053261B2 (en) Simulation method and simulation apparatus
Fraga Filho et al. Smoothed Particle Hydrodynamics
Liang Evaluating shallow water assumptions in dam-break flows
US10031984B2 (en) Method and device for simulating surface tension
JPWO2013042234A1 (ja) 物体運動解析装置、物体運動解析方法、及び物体運動解析プログラム
Dao et al. Numerical modelling of extreme waves by Smoothed Particle Hydrodynamics
Zhou et al. MLPG method based on Rankine source solution for modelling 3D breaking waves
KR101106548B1 (ko) 텐서형 와점성계수를 가진 2차원 하천흐름모형을 이용하여 천수흐름을 해석하는 방법
Ma et al. Numerical study of 2-D vertical Water-entry problems using Two-phase SPH method
Zhang et al. Numerical study of nonlinear shallow water waves produced by a submerged moving disturbance in viscous flow
Bašic et al. A Lagrangian finite difference method for sloshing: simulations and comparison with experiments
Oliaei et al. Some numerical issues using element‐free Galerkin mesh‐less method for coupled hydro‐mechanical problems
Janßen et al. Efficient simulations of long wave propagation and runup using a LBM approach on GPGPU hardware
Rostami Varnousfaaderani et al. Numerical simulation of solitary wave breaking and impact on seawall using a modified turbulence SPH method with Riemann solvers
Torres et al. Three dimensional numerical modelling of full-scale hydraulic structures
Yue et al. A well-balanced and fully coupled noncapacity model for dam-break flooding
Khanpour et al. Numerical simulation of the flow under sluice gates by SPH model
KR101135307B1 (ko) 유한체적법에 의한 보존변수 구성 방법 및 이를 이용한 2차원 유체 흐름 해석 방법
Grave et al. Residual-based variational multiscale 2d simulation of sediment transport with morphological changes
KR101820946B1 (ko) 파티클 기반의 sph에 의한 강우 시뮬레이션 방법 및 이를 실행하도록 구성되어진 강우 시뮬레이션 플랫폼
Sriram et al. Simulation of 2D breaking waves by using improved MLPG_R method

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20150604

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20160809

A521 Written amendment

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20161006

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20161025

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20161107

R150 Certificate of patent or registration of utility model

Ref document number: 6048062

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150