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

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

Info

Publication number
JP5670832B2
JP5670832B2 JP2011116344A JP2011116344A JP5670832B2 JP 5670832 B2 JP5670832 B2 JP 5670832B2 JP 2011116344 A JP2011116344 A JP 2011116344A JP 2011116344 A JP2011116344 A JP 2011116344A JP 5670832 B2 JP5670832 B2 JP 5670832B2
Authority
JP
Japan
Prior art keywords
particle
density
equation
pressure
unit
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
JP2011116344A
Other languages
English (en)
Other versions
JP2012243288A (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
Kobe Steel Ltd
Original Assignee
Fujitsu Ltd
Kobe Steel 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, Kobe Steel Ltd filed Critical Fujitsu Ltd
Priority to JP2011116344A priority Critical patent/JP5670832B2/ja
Priority to US13/433,456 priority patent/US9053261B2/en
Publication of JP2012243288A publication Critical patent/JP2012243288A/ja
Application granted granted Critical
Publication of JP5670832B2 publication Critical patent/JP5670832B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Description

本発明は、シミュレーション方法、シミュレーション装置及びシミュレーションプログラムに関する。
従来、数値計算を用いて水や空気といった流体の流れを調べる連続体の運動の解析、すなわち、流体解析において、粒子法と呼ばれる技術が提案されている。具体的には、粒子法とは、連続体の運動を有限の数の粒子の運動として解析する手法である。現在提案されている代表的な粒子法としては、SPH（Smoothed Particles Hydrodynamics）法やMPS（Moving Particles Semi-implicit）法といったものがある。以下では、水や空気といった流体を「連続体」と呼ぶ場合がある。
粒子法における標準的な手法としては、着目したある粒子に対して予め領域（以下では、「影響領域」と呼ぶ。）を設定しておき、その影響領域の内に存在する相手粒子のみからの相互作用を計算することで着目したその粒子に加えられる力を求める方法が従来知られている。
特に、ＳＰＨ法を用いて連続体を表現する際の特徴は、複数の粒子の物理量をカーネル関数と呼ぶ重み関数を用いて平滑化処理による近似を行い、基礎方程式を離散化することにある。この平滑化処理によって、メッシュ点上で物理量を評価するといった煩雑なメッシュ操作の計算処理が無くなるため、ＳＰＨ法は、自由表面における流れを解析する自由表面問題や異なる支配方程式で表される複数の物理現象を解析するマルチフィジックス問題を扱うのに適しているといえる。
このため、ＳＰＨ法は、例えば、海波が護岸に衝突・越波するときの流量や衝撃圧を推定するための手法として適していると考えられている。
J.J.Monaghan 著 「Smoothed Particle Hydrodynamics」Annu. Rec. Astron. Astrophys., Vol.30, pp.543-574 鈴木幸人、越塚誠一、岡芳明 著 「ＨＭＰＳ（Hamiltonian Moving Particle Semi-implicit）法の開発（第２報，シンプレクティックスキームの適用）」 Trans. JSCES, Paper No.20050017(2005) Paul W. Cleary 著 「Modelling confined multi-material heat and mass flows using SPH」 Appl. Math. Modelling, Vol.22,pp.981-993,1998 M.G.Gesteira, B.D.Rogers, R.A.Dalrymplem, A.J.C.Crespo, M.Narayanaswamy 著 「User Guide for the SPHysics Code, September 2010」 pp.9-15 Fang,J., Parriaux,A., Rentschler,M., Ancey,C. 著 「Improved SPH methods for simulating free surface flows of viscous fluids」 Applied Numerical Mathematics, Vol.59,pp.251-271,2009 Morris,J., Fox,P.J., Zhu,Y. 著 「Modeling Low Reynolds number incompressible flows using SPH」 J.comp. Phys., Vol.136,pp.214-226,1997
しかしながら、従来の粒子法を用いたシミュレーションの場合、表面波の伝播を長時間扱うと、表面波が減衰してしまうという問題点がある。これについては、実際に標準的なＳＰＨ法で波の伝播解析を行ったところ、一定の時間を過ぎると波高の大きな減衰や発散が生じてしまい、造波試験を実施できないことが分かった。このように、従来の粒子法を用いたシミュレーションでは、現実との乖離が大きくなってしまうため、現実の水槽内に設置した造波板により造波を行う水理実験を減らすことができなかった。
開示の技術は、上記に鑑みてなされたものであって、波高が減衰しない現実に近い解析を行うことができ、水理実験の実施回数を減らすことができるシミュレーション方法、シミュレーション装置及びシミュレーションプログラムを提供することを目的とする。
本願の開示するシミュレーション方法、シミュレーション装置及びシミュレーションプログラムは、一つの態様において、連続体を粒子の集合として表現する場合の各粒子の状態を解析するシミュレーション方法において、入力した前記各粒子の初期状態の速度、密度、圧力及び位置の初期値に基づき、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを表すカーネル関数により離散化された、ハミルトン方程式を満たす加速度運動方程式を用いて、各粒子の加速度と各粒子が境界面から受ける反発力とを算出し、前記算出した各粒子の加速度と各粒子が境界面から受ける反発力とを基に、現在の各粒子の速度に対して現在の加速度を時間積分した値を加えることで、単位時間経過後の各粒子の速度を算出し、前記カーネル関数を用いて連続体の密度の時間変化を表す離散化した連続方程式を用いて各粒子の密度変動を算出し、前記算出した各粒子の密度変動に対して、前記単位時間経過後の各粒子の速度を用いて時間による積分を行うことで、前記単位時間経過後の各粒子の密度を算出し、前記各粒子の密度を算出する所定回数毎に、前記単位時間経過後の各粒子の密度に対して平滑化を行い、前記単位時間経過後の各粒子の密度を基に、状態方程式を用いて前記単位時間経過後の各粒子の圧力を算出し、前記単位時間経過後の各粒子の速度を基に時間による積分を行い、前記単位時間経過後の各粒子の位置を算出し、前記速度の算出、前記密度の算出、前記圧力の算出及び前記位置の算出を、前記粒子の全てについて前記初期状態から所定時間が経過するまで繰返し、各単位時間経過後における速度、密度、圧力及び位置を取得する。
本願の開示するシミュレーション方法、シミュレーション装置及びシミュレーションプログラムの一つの態様によれば、波高が減衰しない現実に近い解析を行うことができ、水理実験の実施回数を減らすことができるという効果を奏する。
図１は、実施例に係るシミュレーション装置のブロック図である。 図２は、影響範囲及び近傍粒子について説明するための図である。 図３は、実施例で用いた境界からの反発力を説明するための図である。 図４は、実施例に係るシミュレーションの処理のフローチャートである。 図５は、水理実験に用いた造波水槽の概念図である。 図６は、波高計３０２における波高の解析結果と水理実験との比較を表す図である。 図７は、波高計３０３における波高の解析結果と水理実験との比較を表す図である。 図８は、波高計３０４における波高の解析結果と水理実験との比較を表す図である。 図９は、波高計３０５における波高の解析結果と水理実験との比較を表す図である。 図１０は、波高計３０６における波高の解析結果と水理実験との比較を表す図である。 図１１は、実施例に係るシミュレーション方法におけるカーネル関数、運動方程式又は時間積分スキームを変化させた場合の波高を表す図である。 図１２は、従来手法におけるカーネル関数、運動方程式又は時間積分スキームを変化させた場合の波高を表す図である。 図１３は、垂直護岸を設置した実験条件を示す図である。 図１４は、外洋条件を用いた場合の波高履歴を示す図である。 図１５は、瀬戸内条件を用いた場合の波高履歴を示す図である。 図１６は、シミュレーションプログラムを実行するコンピュータを示す図である。
以下に、本願の開示するシミュレーション方法、シミュレーション装置及びシミュレーションプログラムの実施例を図面に基づいて詳細に説明する。なお、以下の実施例により本願の開示するシミュレーション方法、シミュレーション装置及びシミュレーションプログラムが限定されるものではない。
図１は、実施例に係るシミュレーション装置のブロック図である。図１に示すように、本実施例に係るシミュレーション装置は、ユーザインタフェース１、制御部２、結果記憶部３、速度算出部４、密度算出部５、圧力算出部６、位置算出部７及び平滑化部８を有している。ユーザインタフェース１は、シミュレーション装置を使用するユーザからの数値の入力、ユーザへのシミュレーション結果の出力を行なうインタフェース装置であり、具体的にはキーボード等の入力装置及び表示装置等の出力装置である。
制御部２は、初期状態における各粒子の質量、速度、密度、圧力及び位置、すなわち各粒子の質力、並びに、速度、密度、圧力及び位置の初期値の入力をユーザインタフェース１から受ける。そして、制御部２は、受信した各粒子の質量、速度、密度、圧力及び位置の初期値を速度算出部４へ出力する。また、制御部２は、受信した各粒子の質量及び密度の初期値を密度算出部５へ出力する。また、制御部２は、受信した各粒子の圧力の初期値を圧力算出部６へ出力する。また、制御部２は、受信した各粒子の質量を平滑化部８へ出力する。さらに、制御部２は、受信した各粒子の位置の初期値を位置算出部７へ出力する。
また、制御部２は、シミュレーションの終了のタイミングの指定をユーザインタフェース１から受ける。本実施例では、シミュレーションの終了のタイミングの指定は、ステップ単位で行うものとし、ここではシミュレーションがＬステップ後に終了するよう指定されたものとする。ここで、ステップについて説明する。１ステップとは、本実施例ではΔｔ時間経過した後の状態をいう。すなわち、ｎ（１≦ｎ≦Ｌ）ステップ後とは、ｎΔｔ時間経過した後のことを指す。以下では、初期状態を１ステップ目とし、ｎステップ後のステップのことを「ステップｎ」と言う。すなわち、初期状態は１回目のステップといえる。以下では、既に各物理量の算出が終わった最新のステップをステップｎとし、次に物理量を計算するステップをステップｎ＋１として説明する。
ここで、本実施例では、ステップ単位で終了のタイミングを指定したが、これは他の方法でも良く、例えば、初期状態から終了までの時間を指定してもよい。
また、制御部２は、ユーザインタフェース１より入力を受けた各粒子の位置の初期値から、各粒子の影響範囲内にある他の粒子のリストを作成する。以下では、各粒子の影響範囲内にある他の粒子を「近傍粒子」と呼ぶ。また、近傍粒子のリストを「近傍粒子リスト」と呼ぶ。
図２は、影響範囲及び近傍粒子について説明するための図である。ここで、図２を参照して、各粒子の影響範囲及び近傍粒子について説明する。ここでは、着目粒子と他の粒子との距離が２ｈ以内の場合に、相互に影響を与え合うものとして説明する。ここで、影響とは、例えば、粒子が移動する場合に、他の粒子へ力を加えるなどのことを言う。
ここでは、図２の粒子１０１を例に説明する。粒子１０１の影響範囲１０２は、粒子１０１を中心として半径２ｈの円の中に含まれる範囲である。この半径２ｈは、粒子１０１の影響半径である。そして、影響範囲１０２の中に含まれる他の粒子が、粒子１０１の近傍粒子である。図２では、例えば、粒子１０３のように、斜線で表されている粒子が近傍粒子である。すなわち、粒子１０１が影響を与えることになる他の粒子は、斜線で表される近傍粒子ということになる。逆に、斜線で表される近傍粒子からのみ、粒子１０１は影響を受けることになる。ここで、近傍粒子が数十個程度含まれるように半径２ｈを設定すれば、着目した粒子と他の粒子間の相互作用を解析するには十分であることが経験的に知られている。
また、粒子１０１の位置を表す位置ベクトルをｒａとし、粒子１０３の位置を表す位置ベクトルをｒｂとすると、粒子１０１と粒子１０３との間の距離は｜ｒａ−ｒｂ｜となる。すなわち、｜ｒａ−ｒｂ｜≦２ｈを満たす場合に、粒子１０３は粒子１０１の近傍粒子とされる。
制御部２は、例えば図２の粒子１０１に対する近傍リストを作る場合、まず、影響範囲１０２に含まれる粒子１０３のような斜線で表される各粒子を粒子１０１の近傍粒子として抽出する。例えば、制御部２は、着目粒子である粒子１０１の位置と他の粒子の位置から粒子１０１と各他の粒子との距離を求め、その距離が２ｈ以下となる粒子を近傍粒子として抽出する。
図１に戻って、制御部２は、近傍粒子の抽出を全ての粒子に対して行う。そして、制御部２は、各着目粒子と抽出した各着目粒子の近傍粒子とが対応する近傍粒子リストを作成する。そして、制御部２は、作成した近傍粒子リストを速度算出部４、密度算出部５、圧力算出部６、位置算出部７及び平滑化部８へ出力する。
また、制御部２は、ステップｎ＋１における各着目粒子の速度の入力を速度算出部４から受ける。また、制御部２は、ステップｎ＋１における各着目粒子の密度及び圧力の入力を平滑化部８から受ける。さらに、制御部２は、ステップｎ＋１における各着目粒子の位置の入力を位置算出部７から受ける。ここで、１ステップとは、上述したように、Δｔ時間経過した後の状態である。そこで、例えば、ステップｎ＋１における速度とは、ステップｎにおける速度からΔｔ時間後の速度である。
すなわち、制御部２は、初期状態からΔｔ時間経過した後の、各着目粒子の速度、密度、圧力及び位置の入力を受ける。次に、制御部２は、初期状態から２Δｔ時間経過した後の、各着目粒子の速度、密度、圧力及び位置の入力を受ける。このように、制御部２は、既にステップｎにおける速度、密度、圧力及び位置を受信している場合には、次にステップｎ＋１における各着目粒子の速度、密度、圧力及び位置を受信する。そして、制御部２は、指定されたシミュレーションの終了のタイミングまで、各ステップにおける各着目粒子の速度、密度、圧力及び位置を順次受信する。
さらに、制御部２は、ステップｎ＋１における各着目粒子の最新の速度、密度、圧力及び位置から１ステップ後の情報の入力を受けると、受信したステップｎ＋１における各情報を結果記憶部３へ記憶させる。
また、制御部２は、ステップｎ＋１における各着目粒子の密度、圧力及び位置を、次の計算に用いる各着目粒子の密度、圧力及び位置として、速度算出部４へ出力する。
また、制御部２は、ステップｎ＋１における各着目粒子の密度を次の計算に用いる密度として密度算出部５へ出力する。
また、制御部２は、ステップｎ＋１における各着目粒子の圧力を次の計算に用いる現在の圧力として圧力算出部６へ出力する。
また、制御部２は、ステップｎ＋１における各着目粒子の速度を位置算出部７へ出力する。
さらに、制御部２は、受信したステップｎ＋１における各着目粒子の位置を用いてステップｎ＋１における各着目粒子の近傍粒子を求め、近傍粒子リストを作成する。そして、制御部２は、作成した近傍粒子リストを速度算出部４、密度算出部５、圧力算出部６、位置算出部７及び平滑化部８へ出力する。このように、制御部２は、ステップ毎の速度、密度、圧力及び位置、並びに近傍粒子リストを各部へ出力する。そして、制御部２は、現在のデータを用いて１ステップ後のデータを算出するよう、速度算出部４、密度算出部５、圧力算出部６、位置算出部７及び平滑化部８に対して指示し、ステップＬまで繰り返させる。
結果記憶部３は、メモリやハードディスクなどの記憶装置である。
速度算出部４は、例えば、カーネル関数として、次の式（１）で表される５次のスプライン関数を記憶している。
Figure 0005670832
ここで、αＤは、２次元の計算を行なう場合にはαＤ＝７／（４πｈ２）という値であり、３次元の計算を行なう場合にはαＤ＝２１／（１６πｈ３）という値である。そして、ｒは、カーネル関数を用いた計算の対象としている２つの粒子の距離を表し、ｈは影響範囲を表しており、ｑはｑ＝ｒ／ｈと定義される。
式（１）は、ある粒子における、他の粒子へ影響を及ぼす寄与の大きさを表している。そして、式（１）は、粒子同士が離れると一方の粒子が他方の粒子に影響を及ぼす寄与の大きさが小さくなることを表わしている。
そして、式（１）は、（イ）全空間で積分すると１になる、（ロ）影響範囲が０に近づくと、すなわち、一方の粒子と影響を受ける他の粒子との間の距離が０に近づくと極限値がディラックのデルタ関数になる、（ハ）１階以上微分可能である、という３つの条件を満たしている。
さらに、式（１）は、ｑが０より小さい又は２より大きい場合の値が０になることから、有限の領域の外で値が０となる性質を有している。この性質は、コンパクトサポートと呼ばれる場合がある。
ここで、圧縮場では、ｄ２Ｗ／ｄｑ２が負の領域でtensile instabilityと呼ばれる数値不安定が生じやすい。そこで、カーネル関数としては、ｄ２Ｗ／ｄｑ２が負となる領域が比較的狭く、不安定が発生しにくいと考えられる式（１）で表される５次のスプライン関数を用いることが特に好ましい。ただし、数値的不安定さを抑える要求に応じて、カーネル関数は他のものを用いることもできる。例えば、他のカーネル関数としては、次の式（２）〜（４）のようなものを用いてもよい。
Figure 0005670832
ここで、ｑ＝ｒ／ｈである。また、αＤは、２次元の計算を行なう場合にはαＤ＝１／（πｈ２）という値であり、３次元の計算を行なう場合にはαＤ＝１／（π３／２ｈ３）という値である。式（２）は、ガウシアン関数である。
Figure 0005670832
ここで、αＤは、２次元の計算を行なう場合にはαＤ＝２／（πｈ２）という値であり、３次元の計算を行う場合にはαＤ＝５／（４πｈ３）という値である。
Figure 0005670832
ここで、αＤは、２次元の計算を行なう場合にはαＤ＝１０／（７πｈ２）という値であり、３次元の計算を行なう場合にはαＤ＝１／（πｈ３）という値である。
式（２）〜（４）はいずれも、ある粒子における、他の粒子へ影響を及ぼす寄与の大きさを表している。そして、式（２）〜（４）はいずれも、粒子同士が離れると一方の粒子が他方の粒子に影響を及ぼす寄与の大きさが小さくなることを表わしている。そして、式（２）〜（４）はいずれも、（イ）、（ロ）及び（ハ）の条件を満たしている。
特に、式（２）のガウシアン関数は、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを正規分布として表現している。そして、式（２）を用いた場合、後で詳述する、速度、密度、圧力及び位置の算出において、数学的には厳密な値が求められる。しかし、式（２）を用いた場合、Ｗ（ｒ，ｈ）＝０になる点を求めることができない。この場合、影響の範囲がはっきりしなくなることから、影響半径をどの値に取ったらいいのかが不明確になってしまう。そこで、本実施例では、精度とのバランスを考慮してカーネル関数として式（１）を用いている。
さらに、速度算出部４は、例えば、粒子の加速度を表す離散化された加速度運動方程式として次の式（５）を記憶している。
Figure 0005670832
また、添字ａは、粒子ａについての情報であることを表しており、添字ｂは粒子ｂについての情報であることを表している（以下、添字のａ及びｂは同様であるものとする。）。すなわち、式（５）の運動方程式は、粒子ａの加速度を表している。ここで、粒子ａとは、全ての粒子の中から任意に抽出した粒子である。そして、粒子ｂとは、粒子ａの近傍粒子である。すなわち、和の範囲は粒子ａの近傍粒子について亘ることになる。そして、ρａは粒子ａの密度であり、ρｂは粒子ｂの密度である。そして、ρａｂ（ここでは、ρａｂの上部の横棒の記号を省いて記載している。）はρａとρｂとの算術平均を表している。また、ｐａは粒子ａにかかる圧力であり、ｐｂは粒子ｂにかかる圧力である。また、Ｗａｂは粒子ａと粒子ｂにおけるカーネル関数であり、本実施例では式（１）で表される。さらに、∇ａは、粒子ａの位置におけるベクトル微分演算子を表している。また、ｍｂは粒子ｂの質量である。また、ｇの項（ここで、ｇはベクトルを表しているものとする）は外力項を表している。
式（５）の右辺の第１項は圧力勾配を表している。また、式（５）の右辺の第２項は人工粘性を表している。さらに、式（５）の第３項は重力を表している。
さらに、式（５）の右辺の第２項の人工粘性項は、粒子ａと粒子ｂとが互いに近づく場合に限り粒子同士の衝突を妨げる向きに力を及ぼすための項である。人工粘性項におけるβａｖは、人工粘性の強さを調整するパラメータである。また、人工粘性項におけるφは次の式（６）で表される。
Figure 0005670832
ここで、ｒｂａ＝（ｒｂ−ｒａ）（ここで、ｒｂａ、ｒｂ及びｒａはそれぞれベクトルを表しているものとする。）であり、ｒａ、ｒｂはそれぞれ、粒子ａの位置ベクトルと粒子ｂの位置ベクトルとを表している。すなわち、ｒｂａは粒子ａと粒子ｂとの距離である。ｖｂａ＝（ｖｂ−ｖａ）（ここで、ｖｂａ、ｖｂ及びｖａはそれぞれベクトルを表しているものとする。）である。また、ηは、分母が発散しないためのパラメータである。
ここで、式（５）の導出について説明する。流体の運動方程式は次の式（７）で表される。
Figure 0005670832
ここで、ｖ（ここで、ｖはベクトルを表しているものとする。）は速度ベクトルを表し、ρは流体の密度を表しており、ｐは流体の圧力を表しており、ｇの項（ここで、ｇはベクトルを表しているものとする。）は外力項を表している。また、Δは、ラプラス演算子を表している。そして、式（７）の右辺の第１項は圧力勾配を表し、第２項は粘性応力による力を表し、第３項は外力を表している。
ここで、圧力勾配を表す−∇Ｐ／ρは、圧力Ｐの大きいところから小さなところに向かって加速度が働き、さらに密度ρが大きいほど加速度が小さくなることを表している。さらに、∇Ｐ／ρは、次の式（８）のように変形できる。
Figure 0005670832
ここで、式（８）におけるσは、次に行うように式（７）に対して粒子法を用いて離散化を行った場合に、各粒子について対称な形にする場合を一般化するための調整用のパラメータである。
次に、式（７）の圧力勾配を、式（８）を用いて変形したうえで、ＳＰＨ法を用いて離散化を行う。ここで、ＳＰＨ法では、粒子ａの物理量Ａは、次の式（９）のように離散化表現される。物理量Ａは例えば、圧力や密度などである。
Figure 0005670832
そこで、式（７）の圧力勾配を式（８）を用いて変形した式を、式（９）を用いて離散化表現することで、次の式（１０）のように表すことができる。
Figure 0005670832
ここで、式（１０）では、式（７）における粘性応力の項が省かれている。これは、本実施例で解析の対象としている波などの粘性が低い流体においては、粘性応力が無視してよいほどの力であるためであり、離散化後の式から除いている。この式（１０）で分かるように、運動方程式が、粒子ａ及び粒子ｂのそれぞれに対して対称な式として表されている。すなわち、式（１０）を用いることで、粒子ａと粒子ｂとの間で相互に働く力が等しくなる。このため、式（１０）を用いた解析においては、運動量保存側やエネルギー保存則を守ることが容易となる。
そして、調整用のパラメータをσ＝１とすると、式（５）が求まる。ここで、本実施例では、調整用パラメータをσ＝１とした場合としたが、σの値は他の値をとってもよく、例えば、従来粒子法において一般的に用いられてきたσ＝０を用いてもよい。
さらに、後述するように、本実施例では時間積分をするにあたり、エネルギー散逸が比較的生じにくいSymplectic Euler法を用いる。そして、Symplectic Euler法では、運動方程式は次の式（１１）及び式（１２）のHamilton方程式で記述される。
Figure 0005670832
Figure 0005670832
ここで、ｅａ＝ｍａｖａ（ここで、ｅａはベクトルを表しているものとする。）は、粒子ａの運動量ベクトルである。また、∇ｒは位置に関する微分演算子であり、∇ｅは運動量に関する微分演算子である。また、Ｈは全エネルギ（Hamiltonian）である。
Symplectic Euler法は、Hamilton方程式で記述される系を時間積分するスキームであり、次の式（１３）及び式（１４）で表されるステップで実行される。
Figure 0005670832
Figure 0005670832
ここで、添え字ｎは時間ステップ数を表し、Δｔは時間増分を表している。式（１３）及び式（１４）は両辺にｎ＋１ステップの値を含むので、一般には陰的公式である。しかし、HamiltonianがＨ（ｅ，ｒ）＝Ｔ（ｅ）＋Ｕ（ｒ）とｅのみに依存する項と、ｒのみに依存する項の和として表されるときには、式（１３）は次の式（１５）として表せ、式（１４）は次の式（１６）として表される。
Figure 0005670832
Figure 0005670832
式（１５）及び式（１６）は、右辺に未知の値を含まないため、陽的に解くことができる。
これを、粒子法に適用すると、Hamiltonianは、次の式（１７）のように定義できる。
Figure 0005670832
ここで、Ｎは総粒子数である。また、式（１７）の第１項は運動エネルギーを表し、第２項は粒子間ポテンシャルを表し、第３項は外力ポテンシャルを表している。ただし、∇Ｇａ＝ｇ（ここで、ｇはベクトルを表しているものとする。）である。
そして、式（１７）と定義した場合、式（１２）から式（５）の運動方程式を得ることができる。さらに、式（１３）からｄｒａ／ｄｔ＝ｖａとなる。すなわち、式（５）は、Hamilton systemを形成していることがわかり、式（５）は、Symplectic法を用いる運動方程式として適切であるといえる。
ここで、従来は運動方程式として、例えば、次の式（１８）などが用いられてきた。
Figure 0005670832
ここで、Πａｂは、粘性応力を表す粘性項であり、次の式（１９）のように表される。
Figure 0005670832
ここで、μは粘性係数であり、水の場合はμ＝１０−３（Ｐａ・ｓ）である。また、ξは物理粘性を調整するパラメータであり、例えば、ξ＝４．９６３３３などである。また、式（１９）の右辺の第２項は物理量の不連続に起因する数値誤差を是正する人工粘性項である。
そして、式（１８）で表される運動方程式を用いてシミュレーションを行った場合、波高に現実よりも大きな減衰などが生じてしまう。
これに対して、本実施例における運動方程式である式（５）から人工粘性を省いた運動方程式は、Hamilton方程式で表されるエネルギー保存則から導出される。すなわち、式（５）の運動方程式はエネルギー保存則や運動量保存則を的確に反映することができ、波高の減衰などを抑えることができる。したがって、本実施例における運動方程式である式（５）を用いて解析を行うことで、式（１８）のような従来の式を用いた場合と比較して、波動伝播のような安定性が要求される解析をより適切に行うことができる。
図１に戻って、速度算出部４の説明を続ける。速度算出部４は、運動方程式（５）を用いて加速度ｄｖ／ｄｔを算出する。
次に、速度算出部４は、境界からの反発力を次の式（２０）を用いて求める。
Figure 0005670832
ここで、図３を参照して、式（２０）について説明する。図３は、実施例で用いた境界からの反発力を説明するための図である。ここでは、図３における粒子２００に加えられる境界からの反発力について説明する。
図３に示すように、ｄは境界面２０１からの所定の距離である。そして、境界層２０２、境界面からｄの幅の層を示している。また、ｓは粒子２００の境界層２０２への侵入距離である。
本実施例では、境界面２０１から幅ｄの境界層２０２を考え、境界層２０２内に粒子２００が突入した場合に、粒子２００に対して境界面２０１から粒子２００の向きに反発加速度が加わるものとする。反発力ｆｎは、弾性ばねによるものと同じ形式でｆｎ＝ｋｓとして与える。ばね剛性係数ｋは、音速ｃ０の粒子が有する運動エネルギーｍｃ０ ２／２となり、粒子が境界層２０２へｄの長さだけ進入したときの弾性歪エネルギーｋｄ２／２が同程度となるように定める。すなわち、連続体の運動の解析を行なうシミュレーションにおいては、粒子２００は音速を超えることはないと考えられるので、音速で粒子２００が進入した場合に境界面２０１においてエネルギーがなくなるように設定しておけば、粒子２００は境界面を突き抜けることはなくなる。このため、βを１程度のパラメータとする。このようにβを設定することで、式（２０）が導かれる。
式（２０）で表される反発力を用いることで、境界壁からの粒子のすり抜けを防止し、導入する摩擦やfree slip条件の正確性を向上でき、さらに、境界で物理的に無意味なエネルギー消費の発生を抑えることができる。また、式（２０）は簡単な式で表されているので、式（２０）で表される反発力を用いた場合には、計算負荷を抑えることができる。
ここで、波高の現実よりも大きな減衰などの各種問題が発生しにくい式として本実施例では式（２０）を用いたが、この反発力は他の式を用いてもよい。例えば、反発力が過大になる可能性を考慮して、粒子が境界面に近づいている場合には、弾性ばねを考慮し、逆に粒子が境界面から遠ざかる場合には、弾性ばねを考慮しないように反発力を設定してもよい。このような反発力は一方向弾性ばね法と呼ばれることがある。エネルギー保存の観点からは粒子と境界面の相対速度に関係なく弾性ばねを考慮する式（２０）のような反発力が好ましいが、境界での不自然な反発を避けるために一方向弾性ばね法をもちいてもよい。
そして、速度算出部４は、運動方程式（５）を用いて算出した加速度ｄｖ／ｄｔと式（２０）を用いて求めた反発力とを加えてステップｎにおける加速度Ｆｎ ａ（ここで、Ｆｎ ａは加速度ベクトルを表しているものとする。）を算出する。
次に、速度算出部４は、算出した加速度Ｆｎ ａを用いて次の式（２１）で表される時間積分を行い、ステップｎ＋１における速度を求める。
Figure 0005670832
ここで、ｖｎはステップｎにおける速度を表している。
また、時間増分Δｔは、本実施例では次の式（２２）とした。以下の説明での時間増分Δｔは全て式（２２）とした。
Figure 0005670832
ここで、αΔｔは発散防止用の係数で、本実施例では一例として０．１とした。
速度算出部４は、以上に説明した、運動方程式を用いて算出した加速度に境界面からの反発力を加えてステップｎにおける加速度を算出し、この加速度を用いてステップｎ＋１の速度を算出するという処理を、全ての粒子に対して行う。
そして、速度算出部４は、算出したステップｎ＋１における各粒子の速度を制御部２、密度算出部５及び位置算出部７へ出力する。
密度算出部５は、式（１）で表されるカーネル関数を記憶している。さらに、密度算出部５は、次の式（２３）で表される粒子の密度の時間変化を表す密度変動の連続方程式を離散化した式を記憶している。
Figure 0005670832
ここで、式（２３）の導出について説明する。離散化しない場合の流体の密度の時間変化の連続方程式は次の式（２４）で表される。
Figure 0005670832
この式（２４）に対して、上述した式（９）で表されるＳＰＨ法における粒子ａの物理量Ａを用いて離散化することで式（２３）が得られる。この連続方程式は、各粒子の質量を重みとして、重みを付けたカーネル関数で表される各粒子が有する密度を重ね合わせたものである。
図１に戻って、密度算出部５は、各粒子の質量、ステップｎにおける各粒子の位置及びステップｎにおける近傍粒子リストの入力を制御部２から受ける。また、密度算出部５は、ステップｎ＋１における速度を速度算出部４から受ける。そして、密度算出部５は、受信した各粒子の質量、ステップｎ＋１における各粒子の速度、ステップｎにおける各粒子の位置及びステップｎにおける近傍粒子リストを式（２３）に代入する。そして、密度算出部５は、ステップｎ＋１における密度変動である密度時間微分ｄρａ／ｄｔ＝Ｄａ ｎ＋１を算出する。
次に、密度算出部５は、前回のステップで算出したステップｎにおける密度及び算出したステップｎ＋１における密度変動を用いて次の式（２５）で表される時間積分を行い、ステップｎ＋１の密度を求める。
Figure 0005670832
密度算出部５は、密度変動の算出及び算出した密度変動を用いた密度の算出を全ての粒子に対して行う。
そして、密度算出部５は、算出したステップｎ＋１における各粒子の密度を圧力算出部６及び平滑化部８へ出力する。
平滑化部８は、ステップｎ＋１における各粒子の密度の入力を密度算出部５から受ける。
平滑化部８は、カウンタを有している。そして、平滑化部８は、ステップｎ＋１における各粒子の密度及び圧力の入力の回数をカウンタによりカウントし、２０回に一回（例えば、２０の倍数のステップ毎）、受信したステップｎ＋１における各粒子の密度及び圧力に対して以下の平滑化を行う。この２０回に一回のステップが「所定ステップ」の一例にあたる。ここで、本実施例では、２０回に一回の割合としたが、これは他の割合としてもよく、さらに、定期的でなくてもよい。この平滑化を行うステップは、密度及び圧力の算出による表面における粒子のみだれの状態を考慮して決定されることが好ましい。
平滑化部８は、受信したステップｎ＋１における密度に対して、次の式（２６）を用いて圧力の平滑化を行う。
Figure 0005670832
ここで、ρａ ｎｅｗは、平滑化処理後の密度である。
そして、平滑化部８は、平滑化処理後の密度をステップｎ＋１における密度として圧力算出部６へ出力する。また、平滑化部８は、平滑化を行わなかった場合は、密度算出部５から受信したステップｎ＋１における各粒子の密度を制御部２に出力し、平滑化を行った場合には、平滑化後のステップｎ＋１の各粒子の密度を制御部２に出力する。
後述するように、密度の平滑化に伴い圧力に対しても平滑化が行われる。このように、密度及び圧力の平滑化を行うことで、圧力振動を抑制することができる。
ここで、本実施例では、密度算出部５が算出したステップｎ＋１の各粒子の密度は常に平滑化部８に渡されているが、これは次のような方法でもよい。例えば、密度算出部５は、今回処理を行なうステップが２０回に一回のステップに該当するか否かを判定し、今回処理を行なうステップが２０回に一回のステップの場合には、密度算出部５は、平滑化部８へステップｎ＋１の各粒子の密度を出力する。これに対して、今回処理を行うステップが２０回に一回のステップ以外の場合には、密度算出部５は、制御部２へステップｎ＋１の各粒子の密度を出力する。このような方法を取ることで、平滑化部８における処理を軽減できる。
圧力算出部６は、ステップｎ＋１における各粒子の密度の入力を密度算出部５から受ける。また、２０回に一回のステップにおいて、平滑化されたステップｎ＋１における各粒子の密度の入力を平滑化部８から受ける。そして、圧力算出部６は、平滑化部８からステップｎ＋１における各粒子の密度の入力を受けた場合、圧力算出部６から受信したステップｎ＋１における各粒子の密度を平滑化部８からステップｎ＋１における各粒子の密度に更新する。
そして、圧力算出部６は、次の式（２７）で表される状態方程式にステップｎ＋１における密度を用いて、ステップｎ＋１における圧力を算出する。
Figure 0005670832
ここで、ρ０は平均密度を表す。本実施例では、一例として、ρ０＝１（ｇ／ｃｍ３）、ｃ０＝５０（ｍ／ｓｅｃ）として計算を行った。
ここで、状態方程式により、密度の大きな領域は圧力が大きくなり、その周囲へ向かって圧力勾配が働くため、式（２７）を用いた結果として密度の凹凸がならされる。このように、密度の凹凸がならされることで、密度が均一に近くなり、より数値計算の正確性を向上することができる。
また、密度の揺らぎの程度δｐ＝（ρ／ρ０）−１は、流体の典型的な流速ｖ０と音速を用いて表されるｖ０ ２／ｃ０ ２以下と見積もることができる。そして、本実施例における流速は、音速（５０ｍ／ｓｅｃ）に比べて、十分小さい値であるので、密度の揺らぎも十分小さい値となる。そして、本実施例では流体の体積変化が０．１％程度であったことから、密度変化も同様に０．１％程度であり、数値計算上問題となる値ではない。
圧力算出部６は、ステップｎ＋１における圧力の算出を全ての粒子について行う。そして、圧力算出部６は、算出したステップｎ＋１における各粒子の圧力を制御部２へ出力する。
ここで、本実施例では、圧力算出部６は、常に密度算出部５からステップｎ＋１における各粒子の密度を受信し、平滑化部８から平滑化が施された密度を受信した場合に、ステップｎ＋１における各粒子の密度を平滑化後の値に更新している。しかし、次のような方法でもよい。すなわち、密度算出部５が、今回処理を行なうステップが２０回に一回のステップに該当するか否かを判定し、今回処理を行なうステップが２０回に一回のステップであれば、密度算出部５は、圧力算出部６にはステップｎ＋１における各粒子の密度を出力せず、平滑化部８にのみ出力する。そして、圧力算出部６は、平滑化部８からのみステップｎ＋１における各粒子の密度を受信する。一方、密度算出部５が、今回処理を行なうステップが２０回に一回のステップ以外のステップと判定した場合、密度算出部５は、圧力算出部６と平滑化部８の双方にステップｎ＋１における各粒子の密度を出力する。そして、圧力算出部６は、密度算出部５及び平滑化部８の双方からステップｎ＋１における各粒子の密度を受信する。このような方法を用いることで、圧力算出部６の処理を軽減することができる。
位置算出部７は、各粒子の位置の初期値の入力を制御部２から受ける。さらに、位置算出部７は、ｎ＋１ステップにおける各粒子の速度の入力を速度算出部４から受ける。
そして、位置算出部７は、受信したｎ＋１ステップにおける速度を用いて、次の式（２８）で表される時間積分を行い、ステップｎ＋１における粒子の位置を求める。
Figure 0005670832
位置算出部７は、ステップｎ＋１における粒子の位置の算出を全ての粒子に対して行う。そして、位置算出部７は、算出したステップｎ＋１における各粒子の位置を制御部２へ出力する。
このように、時間積分により次のステップの速度をまず求め、次に、求めた次のステップの速度を用いて時間積分を行い、次のステップの密度、圧力及び位置を算出する時間積分法はSymplectic Euler法と呼ばれることがある。シミュレーションにおいて時間積分を行なう手順を時間積分スキームと称することがあり、このSymplectic Euler法は、Hamilton方程式で記述される系に対する時間積分スキームの一種である。
以上に説明したように、本実施例にかかるシミュレーションでは、まず速度を時間微分し、次のステップの速度を求め、その後に次のステップの密度及び圧力を求め、最後に次のステップの位置を求めることになる。すなわち、本実施例に係るシミュレーションにおいては、ステップｎ＋１における密度及びステップｎ＋１の位置のいずれにおいても、それらを求める場合にステップｎの位置を用いている。すなわち、密度及び位置を求めるタイミングを一致することができる。これに対して、従来一般的に用いられている、leap-frog法は、計算の最初に位置座標を半ステップ進め、以降、位置座標と速度を半ステップずつずらしたまま時間積分を行う方法である。この方法では、密度と位置を求めるタイミングがずれてしまう。この点、密度の算出においては、位置が大きく影響する。そのため、密度と位置とは同じタイミングで算出することが数値計算を適正に行ううえで好ましい。すなわち、本実施例で用いたSymplectic Euler法による時間微分を行うことで、エネルギーの保存を数値計算上正確に保存することができる。
次に、図４を参照して、本実施例に係るシミュレーションの処理の流れについて説明する。図４は、実施例に係るシミュレーションの処理のフローチャートである。ここでは、平滑化部８が、平滑化処理を施す２０回に１回のステップとして、２０の倍数のステップ毎に平滑化を施す場合を例に説明する。
制御部２は、各粒子の質量、並びに各粒子の速度、密度、圧力及び位置の初期値を含む入力データをユーザインタフェース１から取得する（ステップＳ１）。そして、制御部２は、受信した各粒子の質量、速度、密度、圧力及び位置の初期値を速度算出部４へ出力する。また、制御部２は、受信した各粒子の質量及び密度の初期値を密度算出部５へ出力する。また、制御部２は、受信した各粒子の圧力の初期値を圧力算出部６へ出力する。また、制御部２は、受信した各粒子の質量を平滑化部８へ出力する。さらに、制御部２は、受信した各粒子の位置の初期値を位置算出部７へ出力する。
制御部２は、受信した各粒子の位置を用いて近傍粒子リストを作成する（ステップＳ２）。
速度算出部４は、式（５）の運動方程式を用いて粒子の現在のステップにおける加速度を算出する（ステップＳ３）。
次に、速度算出部４は、粒子が境界近傍層に含まれる位置にある境界近傍の粒子か否かを判定する（ステップＳ４）。速度算出部４が境界近傍の粒子で無いと判定した場合（ステップＳ４否定）、ステップＳ６へ進む。
これに対して、境界近傍の粒子である場合（ステップＳ４肯定）、速度算出部４は、式（２０）を用いて境界からの反発力を算出する（ステップＳ５）。そして、速度算出部４は、現在のステップにおける加速度に反発力を加えた粒子の加速度を求める（ステップＳ６）。
次に、速度算出部４は、現在のステップにおける速度及び算出した現在のステップにおける加速度を用いて式（２１）の時間積分を行い、次のステップにおける粒子の速度を算出する。
そして、制御部２は、速度算出部４によって次のステップにおける速度の算出が全ての粒子について行われたか否かを判定する（ステップＳ７）。制御部２は、次のステップにおける速度の算出を行っていない粒子があると判定した場合（ステップＳ７否定）、速度算出部４に対してステップＳ３からステップＳ６までを繰り返させる。
これに対して、全ての粒子について次のステップにおける速度の算出が完了している場合（ステップＳ７肯定）、密度算出部５は、式（２３）の連続方程式を用いて次のステップにおける密度時間微分を算出する（ステップＳ８）。そして、速度算出部４は、算出した次のステップにおける各粒子の速度を、制御部２、密度算出部５及び位置算出部７へ出力する。
次に、密度算出部５は、現在のステップにおける密度及び次のステップにおける密度時間微分を用いて式（２５）の時間積分を行い、次のステップにおける密度を算出する（ステップＳ９）。そして、密度算出部５は、算出した密度を圧力算出部６及び平滑化部８へ出力する。
圧力算出部６は、密度算出部５から受信した次のステップにおける密度を式（２７）の状態方程式に用いて次のステップにおける圧力を算出する（ステップＳ１０）。
位置算出部７は、速度算出部４から受信した次のステップにおける速度を用いて式（２８）の時間積分を行い、次のステップの粒子の位置を算出する（ステップＳ１１）。
平滑化部８は、次のステップにおける密度の入力を密度算出部５から受ける。そして、平滑化部８は、現在のステップが２０の倍数のステップか否かを判定する（ステップＳ１２）。平滑化部８は、現在のステップが２０の倍数のステップで無いと判定した場合（ステップＳ１２否定）、密度算出部５から受信した次のステップを制御部２へ出力し、ステップＳ１５へ進む。
これに対して、現在のステップが２０の倍数のステップの場合（ステップＳ１２肯定）、平滑化部８は、受信した次のステップにおける密度に式（２６）を用いて平滑化を行う（ステップＳ１３）。そして、平滑化部８は、平滑化を施した次のステップの密度を制御部２及び圧力算出部６へ出力する。
そして、圧力算出部６は、平滑化部８から受信した平滑化が施された次のステップの密度を式（２７）の状態方程式に用いて次のステップにおける平滑化された圧力を算出する（ステップＳ１４）。そして、圧力算出部６は、平滑化された圧力に次のステップの圧力を更新する。その後、圧力算出部６は、算出した次のステップにおける圧力を制御部２へ出力する。
そして、制御部２は、密度算出部５、圧力算出部６、位置算出部７及び平滑化部８によって次のステップにおける密度、圧力及び位置の算出が全ての粒子について行われたか否かを判定する（ステップＳ１５）。制御部２は、次のステップにおける密度、圧力及び位置の算出を行っていない粒子があると判定した場合（ステップＳ１５否定）、密度算出部５、圧力算出部６、位置算出部７及び平滑化部８に対してステップＳ８からステップＳ１４までを繰り返させる。
これに対して、全ての粒子について次のステップにおける密度、圧力及び位置の算出が完了している場合（ステップＳ１５肯定）、制御部２は、今回のステップの計算を完了し、計算対象の時刻を１ステップ進める（ステップＳ１６）。
そして、制御部２は、今回のステップの各粒子の速度、密度、圧力及び位置を含む計算結果データをユーザインタフェース１に出力する（ステップＳ１７）。さらに、制御部２は、今回のステップの結果データを結果記憶部３に格納する。
そして、制御部２は、指定された最終ステップが完了したか否かを判定する（ステップＳ１８）。制御部２は、全ステップが完了していないと判定した場合（ステップＳ１８否定）、ステップＳ２へ戻る。これに対して、最終ステップが完了している場合（ステップＳ１８肯定）、制御部２は、シミュレーションの処理を終了する。
〔水理実験との比較〕
次に、本実施例に係るシミュレーションの正確さを調べるため、水理実験との比較を行った場合の結果について説明する。
図５は、水理実験に用いた造波水槽の概念図である。この水理実験では、垂直護岸概要条件を用いて実験を行った。そして、この水理実験では、図５に示す位置に、波高計３０１〜３０６を配置し、波高記録を得た。
図５の護岸４０１は模擬的な防波堤であり、斜面４０２は模擬的な水底から護岸４０１に向かう模擬的な斜面である。そして、造波板４０３が振幅運動を行うことによって、造波がなされる。そして、水槽の波高計３０１から護岸４０１までの距離を１８．５ｍとし、水深を０．６７９ｍとした。また、斜面４０２は、１０分の１の勾配を有している。
そして、この検証では、図５と同様の水槽の波高計３０１から護岸４０１に亘る解析モデルを用いて、上述したシミュレーション方法により解析を行い、波高履歴及び越波量を図５で示される水槽を用いた水理実験の結果と比較した。ここで、シミュレーションに用いた造波板は、図５の水槽の造波板と同じ位置にあると想定して解析を行なっている。そして、波高計３０６において、水槽実験によって得られた波高履歴を再現するよう、造波板４０３の移動速度を無反射造波条件式の定常波項を省略した次の式（２９）を用いて計算を行った。
Figure 0005670832
ここで、ηＤは目標波形、η０は造波板前面の水位変動、ωは角周波数であり、ｘｗａｌｌは造波板４０３の位置座標である。また、Ｂは次の式（３０）を満たす。
Figure 0005670832
さらに、Ｈは静水面高さ、Ｋは波数である。Ｋは次の式（３１）から求めることができる。
Figure 0005670832
ここで、ｇは重力加速度である。
表面波伝播の過程のみを調べるためにまず、護岸を設置しない状態での水理実験、本実施例に係るシミュレーション方法の数値解析結果、及び従来手法による数値解析結果を比較する。実験条件は、波高２０ｃｍ、周期２ｓｅｃである。数値解析における粒子数は１２万１９８０体、初期粒子間隔は１ｃｍ、影響半径（＝２ｈ）は３ｃｍである。従来手法による数値解析としては、次の点が本実施例に係るシミュレーション方法と異なるものを採用した。すなわち、従来手法による数値解析では、（Ａ）カーネル関数として、数式（４）で表されるcubic spline関数を用いる。（Ｂ）運動方程式として式（１８）を用いる。ここで、式（１９）において、βａｖ＝０．１、η＝０．０１ｈとした。（Ｃ）時間積分スキームとしてはleap-frog法を用いる。ここで、leap-frog法とは、計算の最初に位置座標を半ステップ進め、以降、位置座標と速度を半ステップずつずらしたまま時間積分を行う方法である。
図６は、波高計３０２の位置における波高の解析結果と水理実験との比較を表す図である。図７は、波高計３０３の位置における波高の解析結果と水理実験との比較を表す図である。図８は、波高計３０４の位置における波高の解析結果と水理実験との比較を表す図である。図９は、波高計３０５の位置における波高の解析結果と水理実験との比較を表す図である。図１０は、波高計３０６の位置における波高の解析結果と水理実験との比較を表す図である。図６〜図１０における実線５１１〜５５１で表されるグラフが本実施例に係るシミュレーション方法を用いた数値解析結果である。図６〜図１０は、いずれも、横軸が時刻であり、縦軸が波高である。また、図６〜図１０における点線５１２〜５５２で表されるグラフが水理実験による結果である。また、図６〜図１０における一点鎖線５１３〜５５３で表されるグラフが本実施例に係るシミュレーション方法を用いた数値解析結果である。
図６〜図１０から分かるように、従来手法においては、造波板４０３に近い波高計３０２では水理実験と波高がほぼ一致しているが、波の伝播距離が増すにつれて波高が大きく減衰している。そして、波高計３０６の位置では、従来手法で求められた波高が水理実験と大きく乖離している。例えば、図６の点線５１２と一点鎖線５１３とは近似している。しかし、図１０の点線５５２と一点鎖線５５３とは大きなズレが生じている。
これに対して、本実施例に係るシミュレーション方法を用いた数値解析の場合には、図６〜図１０に示すように、造波板４０３から遠くなり波の伝播距離が長くなっても、解析結果が水理実験とほぼ一致している。例えば、図６の点線５１２と実線５１１とは近似しており、さらに図１０の点線５５２と実線５５１とも近似している。すなわち、従来手法と比較して本実施例に係るシミュレーション方法は、より水理実験に近い結果を得ることができる。
次に、カーネル関数、運動方程式又は時間積分スキームによる、数値解析の正確性の向上に対する寄与の度合いを測るため、カーネル関数、運動方程式又は時間積分スキームのそれぞれを変更した場合の計算結果について説明する。図１１は、実施例に係るシミュレーション方法におけるカーネル関数、運動方程式又は時間積分スキームを変化させた場合の波高を表す図である。図１１では、カーネル関数、運動方程式又は時間積分スキームをそれぞれ前述した従来手法における（Ａ）〜（Ｃ）に切り替えている。図１１は、横軸が時刻であり、縦軸が波高である。
ここで、太実線６０１が水理実験を示し、細実線６０２が本発明を示し、点線６０３がカーネル関数を切り替えた場合を示し、一点鎖線６０４が運動方程式を切り替えた場合を示し、二点鎖線６０５が時間積分スキームを切り替えた場合を示している。
図１１から分かるように、太実線６０１と細実線６０２とはもっとも近似している。そして、カーネル関数、運動方程式又は時間積分スキームのそれぞれを切り替えた場合、細実線６０２と比べて、太実線６０１から離れていく。その中で、カーネル関数を切り替えた点線６０３が最も大きく太実線６０１から乖離している。
また、図１２は、従来手法におけるカーネル関数、運動方程式又は時間積分スキームを変化させた場合の波高を表す図である。図１２では、従来手法として（Ａ）〜（Ｃ）を満たす式を用いた場合とし、それら（Ａ）〜（Ｃ）を、本実施例で用いた式（１）のカーネル関数、式（５）の運動方程式及びSymplectic Euler法を用いた時間積分に切り替えている。図１２は、横軸が時刻であり、縦軸が波高である。
ここで、太実線６１１が水理実験を示し、細実線６１２が本発明を示し、点線６１３がカーネル関数を切り替えた場合を示し、一点鎖線６１４が運動方程式を切り替えた場合を示し、二点鎖線６１５が時間積分スキームを切り替えた場合を示している。
図１２から分かるように、太実線６１１と細実線６１２とはもっとも乖離している。そして、カーネル関数、運動方程式又は時間積分スキームのそれぞれを切り替えた場合、細実線６１２と比べて、太実線６１１に近づいていく。その中で、カーネル関数を切り替えた点線６１３が最も太実線６０１に近似している。
すなわち、カーネル関数として、式（１）で表される５次のスプライン関数を用いることが、シミュレーションの正確性の向上に最も効果が大きいといえる。
さらに、護岸からの反射波の影響について調べるために直立護岸を設置した状態での水理実験と本実施例に係るシミュレーション方法による数値解析結果との比較を行う。実験条件として、「瀬戸内条件」及び「外洋条件」の２種類の条件を設定した。図１３は、垂直護岸を設置した実験条件を示す図である。ここでは、「瀬戸内条件」及び「外洋条件」の条件として、図１３に示すそれぞれの値を用いた。
図１４は、外洋条件を用いた場合の波高履歴を示す図である。また、図１５は、瀬戸内条件を用いた場合の波高履歴を示す図である。図１４及び図１５のいずれも横軸が時刻であり、縦軸が波高である。そして、実線７０１及び実線７１１が本実施例に係るシミュレーション方法を用いた数値解析結果である。また、点線７０２及び点線７１２が水理実験による結果である。
図１４において、実線７０１と点線７０２とは波高、周期及び波形のいずれも十分に近似している。また、図１５において、実線７１１と点線７１２とは波高、周期及び波形のいずれも十分に近似している。すなわち、本実施例に係るシミュレーション方法を用いた数値解析では、波高、周期及び波形とも、水理実験が十分に再現できているといえる。
以上に説明したように、本実施例に係るシミュレーション方法を用いた場合と従来の粒子法を用いた数値解析とを比較した場合、後者では水理実験の結果と比較して波高が大きく減衰するが、前者では波高の減衰がほとんど無く、水理実験の波高とよく一致する。また、護岸前面の砕波を生じる非線形性の高い領域においても、本実施例に係るシミュレーション方法を用いた場合、波高が実測と良く一致する。これにより、護岸砕波過程における越波量や衝撃圧については造波水槽を用いた水理実験が行われるが、本実施例に係るシミュレーション方法による数値解析を行えば、水理実験をシミュレーションに置き換えることにより水理実験の回数の低減ができ、さらに、実海域の波浪に対する応答も直接評価できる。したがって、護岸設計コストの低減及び合理化に大きく寄与することができる。
また、上記実施例で説明した各種の処理は、あらかじめ用意されたプログラムをコンピュータで実行することによって実現することができる。そこで、以下では、図１６を用いて、図１に示したシミュレーション装置と同様の機能を有するシミュレーションプログラムを実行するコンピュータの一例を説明する。
図１６は、シミュレーションプログラムを実行するコンピュータを示す図である。図１６に示すように、コンピュータ１０００は、キャッシュ１０１０と、ＣＰＵ（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 (10)

  1. 連続体を粒子の集合として表現する場合の各粒子の状態を解析するシミュレーション方法において、
    入力した前記各粒子の初期状態の速度、密度、圧力及び位置の初期値に基づき、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを表すカーネル関数により離散化された、ハミルトン方程式を満たす加速度運動方程式を用いて、各粒子の加速度と各粒子が境界面から受ける反発力とを算出し、
    前記算出した各粒子の加速度と各粒子が境界面から受ける反発力とを基に、現在の各粒子の速度に対して現在の加速度を時間積分した値を加えることで、単位時間経過後の各粒子の速度を算出し、
    前記カーネル関数を用いて連続体の密度の時間変化を表す離散化した連続方程式を用いて各粒子の密度変動を算出し、
    前記算出した各粒子の密度変動に対して、前記単位時間経過後の各粒子の速度を用いて時間による積分を行うことで、前記単位時間経過後の各粒子の密度を算出し、
    前記各粒子の密度を算出する所定回数毎に、前記単位時間経過後の各粒子の密度に対して平滑化を行い、
    前記単位時間経過後の各粒子の密度を基に、状態方程式を用いて前記単位時間経過後の各粒子の圧力を算出し、
    前記単位時間経過後の各粒子の速度を基に時間による積分を行い、前記単位時間経過後の各粒子の位置を算出し、
    前記速度の算出、前記密度の算出、前記圧力の算出及び前記位置の算出を、前記粒子の全てについて前記初期状態から所定時間が経過するまで繰返し、各単位時間経過後における速度、密度、圧力及び位置を取得する
    ことを特徴とするシミュレーション方法。
  2. 前記シミュレーション方法において、
    前記カーネル関数は、全空間で積分すると１になり、前記粒子のうち、一つの粒子と影響を受ける他の粒子との間の距離が０に近づくと極限値がディラックのデルタ関数になり、さらに１階以上微分可能であることを特徴とする請求項１記載のシミュレーション方法。
  3. 前記シミュレーション方法において、
    前記カーネル関数は、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを正規分布として表現するガウシアン関数であることを特徴とする請求項１又は請求項２記載のシミュレーション方法。
  4. 前記シミュレーション方法において、
    前記カーネル関数は、有限の領域外で０になる５次のスプライン関数であることを特徴とする請求項１又は請求項２記載のシミュレーション方法。
  5. 前記シミュレーション方法において、
    前記離散化された加速度運動方程式は、圧力勾配に応じた加速度を表すとともに、密度が大きいほど加速度が小さくなるように加速度を表すことを特徴とする請求項１〜請求項４のいずれか一項に記載のシミュレーション方法。
  6. 前記シミュレーション方法において、
    前記離散化された加速度運動方程式は、各粒子について対称形となるように調整用パラメータを用いて前記加速度運動方程式を変形したことを特徴とする請求項５記載のシミュレーション方法。
  7. 前記シミュレーション方法において、
    前記調整用パラメータが１であることを特徴とする請求項６記載のシミュレーション方法。
  8. 前記シミュレーション方法において、
    前記調整用パラメータが０であることを特徴とする請求項６記載のシミュレーション方法。
  9. 連続体を粒子の集合として表現する場合の各粒子の状態を解析するシミュレーション装置において、
    前記各粒子のうち、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを表すカーネル関数により離散化された、ハミルトン方程式を満たす加速度運動方程式を用いて、各粒子の加速度と境界面からの反発力とを算出し、前記算出した各粒子の加速度と各粒子が境界面から受ける反発力とを基に、現在の各粒子の速度に対して現在の加速度を時間積分した値を加えることで、単位時間経過後の各粒子の速度を算出する速度算出部と、
    前記カーネル関数を用いて連続体の密度の時間変化を表す離散化した連続方程式を用いて求めた各粒子の密度変動に対して、前記単位時間経過後の各粒子の速度を用いて時間による積分を行うことで、前記単位時間経過後の各粒子の密度を算出する密度算出部と、
    前記密度算出部が各粒子の密度を算出する所定回数毎に、前記単位時間経過後の各粒子の密度に対して平滑化を行う前記平滑化部と、
    前記単位時間経過後の各粒子の密度を基に、状態方程式を用いて前記単位時間経過後の各粒子の圧力を算出する圧力算出部と、
    前記単位時間経過後の各粒子の速度を基に時間による積分を行うことで、前記単位時間経過後の各粒子の位置を算出する位置算出部と、
    前記速度算出部、前記密度算出部、前記圧力算出部、前記平滑化部及び前記位置算出部に対して、初期状態の速度、密度、圧力及び位置の初期値を与え、前記初期状態から所定時間経過後までの各単位時間経過後における速度、密度、圧力及び位置の算出を行わせ、算出された各単位時間経過後における速度、密度、圧力及び位置を取得する制御部と、
    を備えたことを特徴とするシミュレーション装置。
  10. 連続体を粒子の集合として表現する場合の各粒子の状態を解析するシミュレーションプログラムにおいて、
    コンピュータに、
    入力した前記各粒子の初期状態の速度、密度、圧力及び位置の初期値に基づき、前記各粒子のうち、一つの粒子が他の粒子に影響を及ぼす寄与の大きさを表すカーネル関数により離散化された、ハミルトン方程式を満たす加速度運動方程式を用いて、各粒子の加速度と各粒子が境界面から受ける反発力とを算出させ、
    前記算出した各粒子の加速度と各粒子が境界面から受ける反発力とを基に、現在の各粒子の速度に対して現在の加速度を時間積分した値を加えることで、単位時間経過後の各粒子の速度を算出させ、
    前記カーネル関数を用いて連続体の密度の時間変化を表す離散化した連続方程式を用いて各粒子の密度変動を算出させ、
    前記算出した各粒子の密度変動に対して、前記単位時間経過後の各粒子の速度を用いて時間による積分を行うことで、前記単位時間経過後の各粒子の密度を算出させ、
    前記各粒子の密度を算出する所定回数毎に、前記単位時間経過後の各粒子の密度に対して平滑化を行わせ、
    前記単位時間経過後の各粒子の密度を基に、状態方程式を用いて前記単位時間経過後の各粒子の圧力を算出させ、
    前記単位時間経過後の各粒子の速度を基に時間による積分を行うことで、前記単位時間経過後の各粒子の位置を算出させ、
    前記速度の算出、前記密度の算出、前記圧力の算出及び前記位置の算出を、前記各粒子の全てについて前記初期状態から所定時間が経過するまで繰り返させ、各単位時間経過後における速度、密度、圧力及び位置を取得させる
    ことを特徴とするシミュレーションプログラム。
JP2011116344A 2011-05-24 2011-05-24 シミュレーション方法、シミュレーション装置及びシミュレーションプログラム Active JP5670832B2 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2011116344A JP5670832B2 (ja) 2011-05-24 2011-05-24 シミュレーション方法、シミュレーション装置及びシミュレーションプログラム
US13/433,456 US9053261B2 (en) 2011-05-24 2012-03-29 Simulation method and simulation apparatus

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2011116344A JP5670832B2 (ja) 2011-05-24 2011-05-24 シミュレーション方法、シミュレーション装置及びシミュレーションプログラム

Publications (2)

Publication Number Publication Date
JP2012243288A JP2012243288A (ja) 2012-12-10
JP5670832B2 true JP5670832B2 (ja) 2015-02-18

Family

ID=47219809

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2011116344A Active JP5670832B2 (ja) 2011-05-24 2011-05-24 シミュレーション方法、シミュレーション装置及びシミュレーションプログラム

Country Status (2)

Country Link
US (1) US9053261B2 (ja)
JP (1) JP5670832B2 (ja)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6458501B2 (ja) * 2015-01-06 2019-01-30 富士通株式会社 シミュレーションプログラム、シミュレーション方法、およびシミュレーション装置
JP6897477B2 (ja) * 2017-10-10 2021-06-30 富士通株式会社 流体シミュレーションプログラム、流体シミュレーション方法および流体シミュレーション装置
CN108897068B (zh) * 2018-04-20 2020-11-13 南京大学 气象台站面密度度量方法
CN108647449B (zh) * 2018-05-15 2022-01-28 长江水利委员会长江科学院 一种基于絮凝动力学的粘性泥沙运动数值模拟方法
CN112446126B (zh) * 2019-09-03 2023-06-16 南京理工大学 超空泡航行体尾拍运动状态的模拟方法
KR102139815B1 (ko) * 2020-04-27 2020-07-30 한밭대학교 산학협력단 두 파장 이상의 소산계수를 이용한 초미세먼지 및 조대입자 총량 및 크기 추출 방법 및 시스템
KR102357256B1 (ko) * 2020-07-21 2022-02-07 한밭대학교 산학협력단 카메라의 임의 풍경 영상을 이용한 방향 의존 시정거리 및 2 차원 공간 미세먼지 분포 측정 방법
CN116882214B (zh) * 2023-09-07 2023-12-26 东北石油大学三亚海洋油气研究院 基于dfl粘弹性方程的瑞雷波数值模拟方法及系统
CN116992796B (zh) * 2023-09-27 2023-12-22 中国科学技术大学 一种自适应低耗散的sph-hllc黎曼求解器耦合方法
CN116992747B (zh) * 2023-09-28 2024-03-22 深圳十沣科技有限公司 一种基于sph流固耦合的冲击式水轮机动力学分析方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH08314983A (ja) * 1995-05-12 1996-11-29 Hitachi Ltd 回路シミュレーション方法
JP5454693B2 (ja) * 2010-08-24 2014-03-26 富士通株式会社 連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置
US8831916B2 (en) * 2011-05-05 2014-09-09 Siemens Aktiengesellschaft Simplified smoothed particle hydrodynamics

Also Published As

Publication number Publication date
JP2012243288A (ja) 2012-12-10
US20120303338A1 (en) 2012-11-29
US9053261B2 (en) 2015-06-09

Similar Documents

Publication Publication Date Title
JP5670832B2 (ja) シミュレーション方法、シミュレーション装置及びシミュレーションプログラム
Liu et al. Three-dimensional liquid sloshing in a tank with baffles
Davidson et al. Linear parametric hydrodynamic models for ocean wave energy converters identified from numerical wave tank experiments
Fang et al. A numerical study of the SPH method for simulating transient viscoelastic free surface flows
Zheng et al. Incompressible SPH method based on Rankine source solution for violent water wave simulation
Raessi et al. A volume-of-fluid interfacial flow solver with advected normals
KR100984048B1 (ko) 파티클 유체 시뮬레이션에서의 강성체 상호작용 처리 방법
US9557203B2 (en) Simulation method and simulation apparatus
NO322925B1 (no) Fremgangsmate og system for a lose element-modeller ved a bruke flerfase-fysikk
CA2912674A1 (en) Mass exchange model for relative permeability simulation
Liang Evaluating shallow water assumptions in dam-break flows
Aristodemo et al. Assessment of dynamic pressures at vertical and perforated breakwaters through diffusive SPH schemes
Salehizadeh et al. A coupled ISPH-TLSPH method for simulating fluid-elastic structure interaction problems
Zhang et al. An immersed boundary method for simulation of inviscid compressible flows
Yuan et al. An immersed-boundary method for compressible viscous flows and its application in the gas-kinetic BGK scheme
WO2014045416A1 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
Zhang et al. Investigations on the hydroelastic slamming of deformable wedges by using the smoothed particle element method
KR101106548B1 (ko) 텐서형 와점성계수를 가진 2차원 하천흐름모형을 이용하여 천수흐름을 해석하는 방법
Cui et al. Solving 2-D highly nonlinear free-surface problems with an improved smoothed particle hydrodynamics method
Ma et al. A study of point moving adaptivity in gridless method
Rijas et al. Numerical modelling of forced heaving of mono hull and twin hull in particle method
Marichal An immersed interface vortex particle-mesh method
Oliaei et al. Some numerical issues using element‐free Galerkin mesh‐less method for coupled hydro‐mechanical problems
US20150186572A1 (en) Analyzing method and analyzing device
Rostami Varnousfaaderani et al. Numerical simulation of solitary wave breaking and impact on seawall using a modified turbulence SPH method with Riemann solvers

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20140306

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20141015

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20141021

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20141120

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: 20141216

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20141218

R150 Certificate of patent or registration of utility model

Ref document number: 5670832

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

S111 Request for change of ownership or part of ownership

Free format text: JAPANESE INTERMEDIATE CODE: R313117

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350