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

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

Info

Publication number
JP5842992B2
JP5842992B2 JP2014503325A JP2014503325A JP5842992B2 JP 5842992 B2 JP5842992 B2 JP 5842992B2 JP 2014503325 A JP2014503325 A JP 2014503325A JP 2014503325 A JP2014503325 A JP 2014503325A JP 5842992 B2 JP5842992 B2 JP 5842992B2
Authority
JP
Japan
Prior art keywords
particle
term
fluid
velocity
particles
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
JP2014503325A
Other languages
English (en)
Other versions
JPWO2013132597A1 (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
Publication of JPWO2013132597A1 publication Critical patent/JPWO2013132597A1/ja
Application granted granted Critical
Publication of JP5842992B2 publication Critical patent/JP5842992B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N13/00Investigating surface or boundary effects, e.g. wetting power; Investigating diffusion effects; Analysing materials by determining surface, boundary, or diffusion effects
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/25Design optimisation, verification or simulation using particle-based methods
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N13/00Investigating surface or boundary effects, e.g. wetting power; Investigating diffusion effects; Analysing materials by determining surface, boundary, or diffusion effects
    • G01N13/02Investigating surface tension of liquids
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N19/00Investigating materials by mechanical methods
    • G01N19/02Measuring coefficient of friction between materials
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Analytical Chemistry (AREA)
  • Pathology (AREA)
  • Immunology (AREA)
  • General Health & Medical Sciences (AREA)
  • Biochemistry (AREA)
  • Chemical & Material Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Description

本発明は、シミュレーションプログラム、シミュレーション方法及びシミュレーション装置に関する。
近年の計算機の演算能力の向上につれて、計算機シミュレーションの手法も発展してきた。その結果として、様々な応用分野に計算機シミュレーションが用いられるようになってきた。
流体や弾性体等の連続体の問題を解く数値計算手法としては、格子をベースにして微分方程式の近似解を求解する有限差分法、有限要素法、又は有限体積法等が多く用いられてきた。近年では数値計算をCAE(Computer Aided Engineering)などの応用分野で活用するため、これらの数値計算手法も発展し、流体と構造物が相互作用する問題が解かれるようになってきた。
しかしながら、数値計算に格子を用いる手法では、自由表面などの界面の存在する問題や、流体・構造連成問題などの移動境界が発生する場合には、取り扱いが複雑なため、プログラム作成が困難である場合が多い。
これに対して数値計算に格子を用いない粒子法(MPS(Moving Particle Semi-implicit)法、SPH(Smoothed Particle Hydrodynamics)法等)では、移動境界の取り扱いに特別な処置を必要としない。それゆえ、近年粒子法は広く用いられるようになっている。
図1は、粒子法の問題点を説明するための図である。
自由表面などの、変形し移動する境界を容易に扱いたいがために開発された粒子法であるが、離散化方法から、連続体が、粒子11から構成されるただの粒子群となってしまっているため、図1に示すように、連続体の境界面がどこかはっきりしなくなる。そのため、粒子法では表面張力などの境界を陽に扱う必要がある問題に対しては、統一的手法は出来ていない。
差分法や有限要素法では、CSF(Continuum Surface Force)モデルと呼ばれる表面張力の計算方法がある(例えば、非特許文献1参照。)。CSFモデルの粒子法への応用がX.Y. Hu, N.A. Adamsにより行われ、粒子法による表面張力の計算が可能となった(例えば、非特許文献2参照。)。また、他にも粒子間の相互作用として、表面張力を決定する手法(例えば、特許文献1参照。)や、近傍粒子の重心位置からのずれから表面張力を求める手法がある(例えば、非特許文献3、4参照。)。
ここで、粒子法に適用されているCSFモデルによる計算方法の概略を説明する。
図2は、影響半径、およびカーネル関数の重ねあわせによる連続関数を示す図である。
まず、表面張力の算出に用いるcolor関数<Ci>を以下のように定義する。
Figure 0005842992
ここで、下付き添え字iはi番目の粒子11の値を表わし、jはj番目の粒子11の値を表わす。
Figure 0005842992
はそれぞれ粒子11の位置ベクトル、粒子11の質量および密度であり、各粒子11に値が与えられている。
Figure 0005842992
はSPH法で用いられるカーネル関数22の重ね合わせで、3次のスプライン関数23などがよく用いられる。
Figure 0005842992
はカーネル関数22の値が非ゼロとなる球形領域の半径を表し、図2に示したように影響半径21と呼ばれる。
color関数の勾配(▽Ci)はSPH法の標準的な手法であり、以下のように計算される。
Figure 0005842992
上記(式2)を用いて粒子iにかかる表面張力を以下のように定義する。
Figure 0005842992
ここで、
Figure 0005842992
である。
Figure 0005842992
は表面張力係数、
Figure 0005842992
は次元、
Figure 0005842992
は単位行列である。
CSFモデルによる計算方法は、上記(式3)により得られた表面張力を用いて時間発展させる計算方法である。
図3乃至図7は、従来のCSFモデルによる計算方法を用いた結果の経時的変化を示す図である。
図3乃至図7に示すように、CSFモデルによる計算方法は、流体の塊が表面張力の影響を受けながら重力により落下し、平面に衝突する計算方法である。図3乃至図7に示した挙動の計算は、非特許文献2に開示されている計算方法を用いた。図3乃至図7においては、y=0が壁面を表すように鉛直方向をy軸に設定しており、y<0となる位置には流体31が侵入できないようになっている。
図3に示すように、流体31の初期形状は、0.05[m]×0.05[m]の正方形である。また、
Figure 0005842992
[N/m2]、
Figure 0005842992
[kg/m3]、
Figure 0005842992
[m]、
Figure 0005842992
[m/s]、
Figure 0005842992
[s]、粘性率
Figure 0005842992
は0.001[m2/s]、重力加速度
Figure 0005842992
[m/s2]、カーネル関数は3次のスプライン関数を用いた。初期状態の粒子位置は均等格子状におき(一辺当たり33個で総数は1089個)、速度はすべての粒子で0を与える。
計算結果を見ると、図4に示すように、流体31は時刻0.04[s]では滑らかな円形分布になり、図5に示すように、時刻0.14[s]で壁と衝突し大きく変形する。
特開2008−111675号公報
B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti, "Modelling Merging and Fragmentation in Multiphase Flows with SuRFER", Journal of Computaional Physics, 113, 134, 147, (1994) X.Y. Hu, N.A. Adams, "A multi-phase SPH method for macroscopic and mesoscopic flows", Journal of Computaional Physics, 213, pp. 844-861 (2006) T. Hongo, M. Shigeta, S. Izawa, and Y. Fukunishi, "3次元非圧縮SPH法における気液界面に作用する表面張力モデル", 第23回数値流体力学シンポジウム, A8-5 (2009) M. Agawa, M. Shigeta, S. Izawa, and Y. Fukunishi, "斜面上を流下する液体の非圧縮SPHシミュレーション", 第23回数値流体力学シンポジウム, A9-4 (2009)
しかしながら、上述のような従来の表面張力計算手法においては、衝突した流体31が大きく変形した場合に、流体31において不自然な速度分布が発生するなどの挙動を示す場合がある、という問題点がある。
上述の例では、図6に示すように、時刻0.7[s]では形状が崩れ、図7に示すように、不自然な分布の速度分布32が発生している。
1つの側面では、本発明は、不自然な分布の速度分布が発生するなどの挙動を抑制したシミュレーション結果を出力することが可能なシミュレーションプログラム、シミュレーション方法及びシミュレーション装置を提供することを目的とする。
1つの案では、流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーションプログラムにおいて、
コンピュータに、
格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
更新した該粒子の速度に基づいて、前記流体の形状変化を算出し、
算出された前記形状変化を出力部に出力する
処理を実行させるシミュレーションプログラムであって、
前記該各粒子の速度の更新においては、
Figure 0005842992
を用い、前記各粒子に係る力は、上記(式12)の右辺第2項で示される圧力勾配項、右辺第3項で示される粘性応力項、右辺第4項で示される重力項、右辺第5項は粒子間距離が近づきすぎないようにするための力であり、右辺第6項は表面張力であり、ここで、
Figure 0005842992
Figure 0005842992
Figure 0005842992
また、下付き添え字は粒子の番号を表し、
Figure 0005842992
:i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
Figure 0005842992
:j番目の質量であり、
Figure 0005842992
:粘性項を計算するために導入した定数であり、
Figure 0005842992
であり
Figure 0005842992
Figure 0005842992
y=0の平面が固体の壁面の場合、
Figure 0005842992
であり、
g:重力加速度
y:y方向の単位ベクトル
h:カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
W:SPH法で用いられるカーネル関数の重ね合わせ
Figure 0005842992

Figure 0005842992
の再標準化行列であり、
Figure 0005842992
とおくと、
Figure 0005842992
Figure 0005842992
であり、
Figure 0005842992
は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
まず、粒子i,jに関して、以下の特性量を定め、
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
さらに、各勾配を以下のよう計算し、
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
ここで、
Figure 0005842992
Figure 0005842992
とし、上記の量を用いて、
Figure 0005842992
は以下のように決定する
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
というシミュレーションプログラムが提供される。
不自然な分布の速度分布が発生するなどの挙動を抑制したシミュレーション結果を出力することができる。
粒子法の問題点を説明するための図である。 影響半径、およびカーネル関数の重ねあわせによる連続関数を示す図である。 従来のCSFモデルによる計算方法を用いた結果の経時的変化(時刻=0[s])を示す図である。 従来のCSFモデルによる計算方法を用いた結果の経時的変化(時刻=0.04[s])を示す図である。 従来のCSFモデルによる計算方法を用いた結果の経時的変化(時刻=0.14[s])を示す図である。 従来のCSFモデルによる計算方法を用いた結果の経時的変化(時刻=0.7[s])を示す図である。 従来のCSFモデルによる計算方法を用いた結果の経時的変化(時刻=0.7[s]:不自然な速度分布発生)を示す図である。 表面張力係数を説明するための図である。 本発明を適用したシミュレーション装置の構成例を示す図である。 表面エネルギーΧ(カイ)を幅εで滑らかに近似した関数を説明するための図である。 2次元流体の例を示す図である。 粒子が壁面に衝突した状態を説明するための図である。 本発明を適用したシミュレーション処理の流れを示すフローチャートである。 流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化(時刻=0[s])を示す図である。 流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化(時刻=0.04[s])を示す図である。 流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化(時刻=0.14[s])を示す図である。 流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化(時刻=0.7[s])を示す図である。 流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化(時刻=0.7[s]:不自然な速度分布の発生無し)を示す図である。 第2の実施の形態における近似関数の効果を示す図である。 情報処理装置の構成図である。
以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。
本発明を適用した実施の形態は、コンピュータに実行させるシミュレーションプログラム、シミュレーション方法及びシミュレーション装置であって、シミュレーションの対象である流体を粒子の集合として捉え、変分を基にした粒子間の相互作用力である表面張力を、所定の近似関数を用いて算出し、その算出した表面張力を用いることにより、流体の形状変化をシミュレーションする。
(第1の実施の形態)
図8は、表面張力係数を説明するための図である。
一般に流体の表面エネルギーは以下のように定義できる。
Figure 0005842992
そして、上記の表面エネルギーを粒子法により離散化して計算し、第一変分より流体粒子にかかる表面張力を求めることができる。
ここで、上記(式4)の右辺の積分領域は、流体の境界面全域であり、
Figure 0005842992
は流体が空気に触れている部分では1、それ以外の部分では0の値をとる。
Figure 0005842992
はそれぞれ、図8に示すように、空気に触れている部分の表面張力係数、固体に触れている部分の表面張力係数である。
図9は、本発明を適用したシミュレーション装置の構成例を示す図である。
図9において、シミュレーション装置900は、処理部901、格納部902及び出力部903を備える。そして、シミュレーション装置900は、後述する近似間数を用い、流体等の連続体を粒子の集まりとして表現する粒子法による数値計算を行って表面張力を算出する。その計算結果を用いて連続体の形状変化を算出し、シミュレーション結果912を出力する。
格納部902は、本発明を適用したシミュレーションプログラムを実行する場合に用いられる各計算式の情報を格納する。
処理部901は、第1および第2の実施の形態として後述する、本発明を適用したシミュレーション処理を実行する。
出力部903は、処理部901が実行したシミュレーション結果912を出力する。
まず、近似関数を用いて表面張力を算出する手法について説明する。
上述の(式4)の積分を、下記(式5)の近似関数のように定義する。
Figure 0005842992
ここで、
Figure 0005842992
は定数であり、表面張力波の振動周期と合うように決められる。
Figure 0005842992
は空間の次元、
Figure 0005842992
はそれぞれ、液体が空気と触れているときの表面張力と液体が固体と触れているときの表面張力である。
Figure 0005842992
は粒子の位置から決まる関数で、
Figure 0005842992
Figure 0005842992
Figure 0005842992
とする。
図10は、表面エネルギーΧ(カイ)を幅εで滑らかに近似した関数を説明するための図である。
(式5)中の
Figure 0005842992
は(式4)の
Figure 0005842992
を幅
Figure 0005842992
で、図10に示したように滑らかに近似した関数であり、
Figure 0005842992
は粒子iと固体境界との距離である。
Figure 0005842992
が面積の次元をもち、さらにSPH法では
Figure 0005842992
は流体の内部の粒子では1に近く、表面近くになると値が小さくなり、影響半径h内に1つも粒子がなくなると、
Figure 0005842992
という値となる。
(式5)のように定義した近似関数は、流体表面の粒子の長さと同じ次元を持つ量を足し合わせたものであり、(式4)の表面積分を近似すると見なす。
(式5)で粒子法により離散化された表面エネルギーの勾配を計算し、各粒子にかかる正味の表面張力を導出すると以下のようになる。
Figure 0005842992
以上、近似関数を用いた表面張力の算出について説明した。
次に、上述のようにして算出した表面張力を用いて、流体の形状変化の算出について説明する。
図11は、2次元流体の例を示す図であり、壁面となる平面から鉛直方向にy軸を設定するとともに、壁面に沿ってx軸を設定している。
図11のような状況下、すなわちy=0の固体111の平面の上に置かれた流体112での、非圧縮性粘性流体の運動方程式を考える。
Figure 0005842992
Figure 0005842992
Figure 0005842992
において、(式7)は質量保存則、(式8)は運動量保存則、(式9)は状態方程式である。
Figure 0005842992
はそれぞれ、流体112の密度場、速度場、圧力場、音速である。
Figure 0005842992
は粘性応力テンソルで、流体112の粘性係数を
Figure 0005842992
(定数)とすると、
Figure 0005842992
である。
(式7)は、流体112が集まってくるような速度場の場合に密度が上昇する効果を表わし、逆に流体112が離れていくような速度場の場合に密度が低下する効果を表わす。(式8)の右辺第1項は圧力勾配項で、流体112が圧力の大きい部分から圧力の小さい部分へ向かって力が発生する効果を表す。右辺第2項は粘性応力項で、流体の流れにブレーキがかかるような効果をあらわす。右辺第3項は重力項である。
Figure 0005842992
は重力加速度、
Figure 0005842992
はy方向の単位ベクトルである。
本第1の実施の形態では(式9)のような密度と圧力のみの関係を用いているが、一般の温度や内部エネルギー、エントロピー等を用いた状態方程式を用いても良い。
また、空気と接触部には
Figure 0005842992
の表面張力係数(定数)、固体111との接触部には
Figure 0005842992
の表面張力係数(定数)がかかっているものとする。
そして、(式7)乃至(式9)を、SPH法を用いて以下のように離散化し、上述した表面張力項を導入する。
Figure 0005842992
ここで、(式10)で計算した
Figure 0005842992
のx成分、y成分をそれぞれ
Figure 0005842992

Figure 0005842992
とし、i番目の速度の成分を
Figure 0005842992
として、以下の計算を行う。
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
ここで、下付き添え字は粒子の番号を表す。すなわちi番目の粒子の位置ベクトル、速度ベクトル、密度、圧力はそれぞれ
Figure 0005842992
である。
Figure 0005842992
はj番目の質量である。
Figure 0005842992
は粘性項を計算するために導入した定数である。
(式12)は(式8)の運動方程式を粒子法により離散化したもので、右辺第2項は圧力勾配項、第3項は粘性応力項を表したものである。(式12)の右辺第5項は粒子間距離が近づきすぎないようにするための力で、
Figure 0005842992
を定数とすれば、
Figure 0005842992
となる。
右辺第6項
Figure 0005842992
は、上述したような表面張力による力で、(式6)が適用されるが、上述のシミュレーションの場合y=0の平面が固体111の壁面であるため、
Figure 0005842992
である。
図12は、粒子が壁面に衝突した状態を説明するための図である。
上述の(式12)および(式14)は、粒子の位置を移動させ、y=0の平面の位置にある壁面を通過しようとする場合に、図12の上に示したように壁面をそのまま通過するのではなく、図12の下に示したように壁面に衝突した瞬間に止まるという現実に即した修正を表す。
ここで、
Figure 0005842992

Figure 0005842992
の再標準化行列である。さらに、
Figure 0005842992
とおく。
Figure 0005842992

Figure 0005842992
である。
Figure 0005842992
は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値である。具体的には以下のように決定する。
まず、粒子i,jに関して、以下の特性量を定める。
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
さらに、各勾配を以下のよう計算する。
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
ここで、
Figure 0005842992
Figure 0005842992
とし、上記の量を用いて、
Figure 0005842992
は以下のように決定する。
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
Figure 0005842992
図13は、本発明を適用したシミュレーション処理の流れを示すフローチャートである。
まず、ステップS1301において、シミュレーションの対象とする流体の物理量のデータとして、各粒子の位置ベクトル、速度ベクトルを取得し、ステップS1302において、上述の(式9)を用いて、現在の速度を用いて粒子の位置を時間dt/2更新する。
そして、(式11)におけるyの値による場合分けにおいて、y<0となる結果、粒子が壁面を通り抜けると判断された場合(ステップS1303:Yes)は、ステップS1304において、上述の(式11)を用いて、粒子が壁面を通り抜けないよう、すなわち、y≧0となるように粒子の位置と速度を補正する。
一方、(式11)におけるyの値による場合分けにおいて、y≧0となる結果、粒子が壁面を通り抜けていないと判断された場合(ステップS1303:No)は、粒子の位置と速度は補正しない。
次に、ステップS1305において、上述の(式12)右辺の第2項乃至第5項を算出することにより、物理モデルにより流体粒子にかかる力を算出する。
そして、ステップS1306において、上述の(式12)右辺の第6項を算出することにより、各流体粒子に係る表面張力を算出する。
各境界粒子に対する表面張力の計算が終了したら、ステップS1307において、流体粒子の速度を更新する。
次に、ステップS1308において、上述の(式13)を用いて、ステップS1307で更新した速度を用いて流体粒子の位置を時間dt/2経過後の位置に更新する。
そして、(式14)におけるyの値による場合分けにおいて、y<0となる結果、粒子が壁面を通り抜けると判断された場合(ステップS1309:Yes)は、ステップS1310において、上述の(式14)を用いて、粒子が壁面を通り抜けないよう、すなわち、y≧0となるように粒子の位置と速度を補正する。
一方、(式14)におけるyの値による場合分けにおいて、y≧0となる結果、粒子が壁面を通り抜けていないと判断された場合(ステップS1303:No)は、粒子の位置と速度は補正しない。
その後、ステップS1311において、流体粒子の速度を更新し、上記ステップS1302乃至ステップS1311を繰り返し実行する。
以上、説明したシミュレーション処理を実行すると、流体の形状は図14乃至図18に示したようになる。
図14乃至図18は、流体の形状について、本発明を適用したシミュレーション処理を用いた結果の経時的変化を示す図である。なお、図3乃至図7の従来技術での計算と同様の計算条件である。
図14に示すように、流体141の初期形状は、0.05[m]×0.05[m]の正方形である。また、それぞれの初期条件を、
Figure 0005842992
[N/m2]、
Figure 0005842992
[kg/m3]、
Figure 0005842992
[m]、
Figure 0005842992
[m/s]、
Figure 0005842992
[s]、
Figure 0005842992
、粘性率
Figure 0005842992
は0.001[m2/s]、重力加速度
Figure 0005842992
[m/s2]、
Figure 0005842992

Figure 0005842992
[m]とした。
初期状態の粒子は均等格子状に配置し(一辺当たり33個で総数は1089個)、速度はすべての粒子で0と密度は1000[kg/m3]を与える。
計算結果を見ると、図15に示すように、流体141は時刻0.04[s]では図4と同様に均一な円形分布になる。そして、図16に示すように、時刻0.14[s]で壁と衝突し大きく変形する。
図17に示したように、時刻0.7[s]では、図6に示した従来手法の計算結果とは異なり、流体141の形状は崩れず、図18に示すように、速度分布も落ち着き定常状態が維持されている。
(第2の実施の形態)
本発明を適用した第2の実施の形態では、上述の第1の実施の形態で、(式4)の積分を(式5)の近似関数のように近似した代わりに、(式4)の積分を下記(式31)の近似関数のように近似する。第1の実施の形態中の圧力「p1」を「p2」に置き換えることにより、第1の実施の形態中の説明を第2の実施の形態に適用できる。
Figure 0005842992
ここで
Figure 0005842992
Figure 0005842992
Figure 0005842992
である。
Figure 0005842992
も面積と同様の次元を持つため、(式6)は表面エネルギーを近似したものと見なす。
Figure 0005842992
は再標準化行列であり、粒子の非一様な分布における勾配計算を補正する。
図19は、第2の実施の形態における近似関数の効果を示す図である。
上記(式31)中の
Figure 0005842992
192は、
Figure 0005842992
191を原点付近で滑らかにした近似関数で、
Figure 0005842992
極小値の大きさが、
Figure 0005842992
とほぼ同じであり、「0.1」となっている。
(式6)のエネルギーを基にして粒子kにかかる正味の力を導出すると以下のようになる。
Figure 0005842992
Figure 0005842992
はそれぞれ、x,y,z方向の単位ベクトルである。
上述したように、本発明は粒子法の計算で、表面張力が効果的な流体運動の計算に対して適用可能である。特に溶けた金属の流し込みのシミュレーションや、樹脂の流し込みのシミュレーション等には効果的である。また、壁面等の固体に衝突するような流体の形状変化をシミュレーションする際に、従来のように壁面等に粒子を置く設定を行う必要が無くなる。
図9のシミュレーション装置は、例えば、図20に示すような情報処理装置(コンピュータ)を用いて実現することが可能である。図20の情報処理装置は、CPU(Central Processing Unit)2001、メモリ2002、入力装置2003、出力装置2004、外部記憶装置2005、媒体駆動装置2006及びネットワーク接続装置2007を備える。これらはバス2008により互いに接続されている。
メモリ2002は、例えば、ROM(Read Only Memory)、RAM(Random Access Memory)、フラッシュメモリ等の記憶装置であり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、CPU2001は、メモリ2002を利用してプログラムを実行することにより、上述のシミュレーション処理を行う。メモリ2002は、図9の格納部902としても使用できる。
入力装置2003は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置2004は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置2004は、図9の出力部903としても使用できる。
外部記憶装置2005は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記憶装置2005には、ハードディスクドライブも含まれる。情報処理装置は、この外部記憶装置2005にプログラム及びデータを格納しておき、それらをメモリ2002にロードして使用することができる。
媒体駆動装置2006は、可搬型記録媒体2009を駆動し、その記録内容にアクセスする。可搬型記録媒体2009は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体2009には、Compact Disk Read Only Memory(CD−ROM)、Digital Versatile Disk(DVD)、Universal Serial Bus(USB)メモリ等も含まれる。オペレータは、この可搬型記録媒体2009にプログラム及びデータを格納しておき、それらをメモリ2002にロードして使用することができる。
このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ2002、外部記憶装置2005、及び可搬型記録媒体2009のような、物理的な(非一時的な)記録媒体が含まれる。
ネットワーク接続装置2007は、通信ネットワーク2010に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置2007を介して受け取り、それらをメモリ2002にロードして使用することができる。ネットワーク接続装置2007は、図9の出力部903としても使用できる。
開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。
例えば、(式5)または(式31)のように定義した近似関数の代わりに、適切な近似関数を用いることもできる。
すなわち、本発明は、以上に述べた実施の形態に限定されるものではなく、本発明の要旨を逸脱しない範囲内で種々の構成または形状を取ることができる。
そして、開示されたシミュレーションプログラム、シミュレーション方法及びシミュレーション装置によれば、シミュレーション対象である流体について、不自然な速度分布が発生するなどの挙動を示すことなく、適切なシミュレーション結果を出力する、という効果を奏する。

Claims (5)

  1. 流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーションプログラムにおいて、
    コンピュータに、
    格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
    前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
    前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
    更新した該粒子の速度に基づいて、前記流体の形状変化を算出し、
    算出された前記形状変化を出力部に出力する
    処理を実行させるシミュレーションプログラムであって、
    前記該各粒子の速度の更新においては、
    Figure 0005842992
    を用い、前記各粒子に係る力は、上記(式12)の右辺第2項で示される圧力勾配項、右辺第3項で示される粘性応力項、右辺第4項で示される重力項、右辺第5項は粒子間距離が近づきすぎないようにするための力であり、右辺第6項は表面張力であり、ここで、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    また、下付き添え字は粒子の番号を表し、
    Figure 0005842992
    :i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
    Figure 0005842992
    :j番目の質量であり、
    Figure 0005842992
    :粘性項を計算するために導入した定数であり、
    Figure 0005842992
    であり
    Figure 0005842992
    Figure 0005842992
    y=0の平面が固体の壁面の場合、
    Figure 0005842992
    であり、
    g:重力加速度
    y:y方向の単位ベクトル
    h:カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
    W:SPH法で用いられるカーネル関数の重ね合わせ
    Figure 0005842992

    Figure 0005842992
    の再標準化行列であり、
    Figure 0005842992
    とおくと、
    Figure 0005842992
    Figure 0005842992
    であり、
    Figure 0005842992
    は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
    まず、粒子i,jに関して、以下の特性量を定め、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    さらに、各勾配を以下のよう計算し、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    ここで、
    Figure 0005842992
    Figure 0005842992
    とし、上記の量を用いて、
    Figure 0005842992
    は以下のように決定する
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    ことを特徴とするシミュレーションプログラム。
  2. 前記表面張力の算出は、空気との界面の表面エネルギーを算出した結果と前記空気以外との界面の表面エネルギーを算出した結果の和を算出することを特徴とする請求項1記載のシミュレーションプログラム。
  3. 前記形状変化の算出は、前記流体が他相を通り抜けようとする場合、前記表面張力の算出において用いられた前記粒子の位置及び速度を補正することを特徴とする請求項1記載のシミュレーションプログラム。
  4. 流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーション方法において、
    コンピュータが、
    格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
    前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
    前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
    更新した該粒子の速度に基づいて、前記流体の形状変化を算出し、
    算出された前記形状変化を出力部に出力する
    シミュレーション方法であって、
    前記該各粒子の速度の更新においては、
    Figure 0005842992
    を用い、前記各粒子に係る力は、上記(式12)の右辺第2項で示される圧力勾配項、右辺第3項で示される粘性応力項、右辺第4項で示される重力項、右辺第5項は粒子間距離が近づきすぎないようにするための力であり、右辺第6項は表面張力であり、ここで、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    また、下付き添え字は粒子の番号を表し、
    Figure 0005842992
    :i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
    Figure 0005842992
    :j番目の質量であり、
    Figure 0005842992
    :粘性項を計算するために導入した定数であり、
    Figure 0005842992
    であり
    Figure 0005842992
    Figure 0005842992
    y=0の平面が固体の壁面の場合、
    Figure 0005842992
    であり、
    g:重力加速度
    y:y方向の単位ベクトル
    h:カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
    W:SPH法で用いられるカーネル関数の重ね合わせ
    Figure 0005842992

    Figure 0005842992
    の再標準化行列であり、
    Figure 0005842992
    とおくと、
    Figure 0005842992
    Figure 0005842992
    であり、
    Figure 0005842992
    は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
    まず、粒子i,jに関して、以下の特性量を定め、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    さらに、各勾配を以下のよう計算し、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    ここで、
    Figure 0005842992
    Figure 0005842992
    とし、上記の量を用いて、
    Figure 0005842992
    は以下のように決定する
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    ことを特徴とするシミュレーション方法。
  5. 流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーション装置において、
    前記流体を形成する粒子の位置ベクトル及び速度ベクトルを格納する格納部と、
    前記流体を粒子の集合で表現した流体モデルについて、前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出
    する表面張力算出部と、
    前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新する更新部と、
    更新した該粒子の速度に基づいて、前記流体の形状変化を算出する形状変化算出部と、
    算出された前記形状変化を出力する出力部と、
    を備えるシミュレーション装置であって、
    前記該各粒子の速度の更新においては、
    Figure 0005842992
    を用い、前記各粒子に係る力は、上記(式12)の右辺第2項で示される圧力勾配項、右辺第3項で示される粘性応力項、右辺第4項で示される重力項、右辺第5項は粒子間距離が近づきすぎないようにするための力であり、右辺第6項は表面張力であり、ここで、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    また、下付き添え字は粒子の番号を表し、
    Figure 0005842992
    :i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
    Figure 0005842992
    :j番目の質量であり、
    Figure 0005842992
    :粘性項を計算するために導入した定数であり、
    Figure 0005842992
    であり
    Figure 0005842992
    Figure 0005842992
    y=0の平面が固体の壁面の場合、
    Figure 0005842992
    であり、
    g:重力加速度
    y:y方向の単位ベクトル
    h:カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
    W:SPH法で用いられるカーネル関数の重ね合わせ
    Figure 0005842992

    Figure 0005842992
    の再標準化行列であり、
    Figure 0005842992
    とおくと、
    Figure 0005842992
    Figure 0005842992
    であり、
    Figure 0005842992
    は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値であり、以下のよ
    うに決定し、
    まず、粒子i,jに関して、以下の特性量を定め、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    さらに、各勾配を以下のよう計算し、
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    ここで、
    Figure 0005842992
    Figure 0005842992
    とし、上記の量を用いて、
    Figure 0005842992
    は以下のように決定する
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    Figure 0005842992
    ことを特徴とするシミュレーション装置。
JP2014503325A 2012-03-06 2012-03-06 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置 Active JP5842992B2 (ja)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2012/055697 WO2013132597A1 (ja) 2012-03-06 2012-03-06 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置

Publications (2)

Publication Number Publication Date
JPWO2013132597A1 JPWO2013132597A1 (ja) 2015-07-30
JP5842992B2 true JP5842992B2 (ja) 2016-01-13

Family

ID=49116113

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2014503325A Active JP5842992B2 (ja) 2012-03-06 2012-03-06 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置

Country Status (4)

Country Link
US (1) US9170185B2 (ja)
EP (1) EP2824598A1 (ja)
JP (1) JP5842992B2 (ja)
WO (1) WO2013132597A1 (ja)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109960841B (zh) * 2017-12-26 2022-11-01 中国科学院深圳先进技术研究院 一种流体表面张力的仿真方法、终端设备及存储介质
CN117951973B (zh) * 2024-03-26 2024-06-07 中国科学技术大学 一种稳定的弹塑性实体光滑粒子动力学模拟方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH1015479A (ja) * 1996-07-09 1998-01-20 Matsushita Electric Ind Co Ltd 自由表面を有する液体のシミュレーションの実行方法及びそれを用いた塗布方法
JP3566009B2 (ja) * 1996-12-27 2004-09-15 松下電器産業株式会社 粒子型流体シミュレーション方法及びその装置
US7204252B2 (en) * 2001-12-21 2007-04-17 Eidon, Llc Surface energy assisted fluid transport system
CN1764917A (zh) 2003-02-05 2006-04-26 莫尔德弗洛爱尔兰有限公司 采用混合模型进行过程模拟的设备和方法
US7647214B2 (en) * 2004-12-27 2010-01-12 Seoul National University Industry Foundation Method for simulating stable but non-dissipative water
JP2007122269A (ja) * 2005-10-26 2007-05-17 Sony Corp 流体−構造体の連成数値シミュレーション方法及び流体−構造体の連成数値シミュレーション用記憶装置のプログラム
JP4798661B2 (ja) 2006-10-27 2011-10-19 みずほ情報総研株式会社 流体解析装置、流体解析方法及び流体解析プログラム

Also Published As

Publication number Publication date
EP2824598A4 (en) 2015-01-14
EP2824598A1 (en) 2015-01-14
US20140365145A1 (en) 2014-12-11
US9170185B2 (en) 2015-10-27
WO2013132597A1 (ja) 2013-09-12
JPWO2013132597A1 (ja) 2015-07-30

Similar Documents

Publication Publication Date Title
Bender et al. Divergence-free SPH for incompressible and viscous fluids
JP6657359B2 (ja) ハイブリッド熱格子ボルツマン法のための温度結合アルゴリズム
Wang et al. Algorithms for interface treatment and load computation in embedded boundary methods for fluid and fluid–structure interaction problems
Vincent et al. Facilitating the adoption of unstructured high-order methods amongst a wider community of fluid dynamicists
EP3525120A1 (en) Lattice boltzmann collision operators enforcing isotropy and galilean invariance
JP6249912B2 (ja) 解析装置
JP5644872B2 (ja) シミュレーション装置、シミュレーション方法、及びプログラム
US10031984B2 (en) Method and device for simulating surface tension
EP2919140A1 (en) Apparatus, method, and computer readable storage medium storing a program for simulating injection molding
JP6098190B2 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
Manteaux et al. Adaptive physically based models in computer graphics
JP5704246B2 (ja) 物体運動解析装置、物体運動解析方法、及び物体運動解析プログラム
Haji Mohammadi et al. Moving least squares reconstruction for sharp interface immersed boundary methods
Karakus et al. An adaptive fully discontinuous Galerkin level set method for incompressible multiphase flows
Nguyen et al. A discontinuous Galerkin front tracking method for two-phase flows with surface tension
JP5842992B2 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
JP5839473B2 (ja) 解析装置
Lygidakis et al. Numerical analysis of flow over the NASA Common Research Model using the academic Computational Fluid Dynamics code Galatea
JP2017194884A (ja) 解析装置および解析方法
Nestor et al. Moving boundary problems in the finite volume particle method
Li An arbitrary Lagrangian Eulerian method for three-phase flows with triple junction points
Zhou et al. Extension of local domain-free discretization method to simulate 3D flows with complex moving boundaries
Wang et al. Extended Eulerian SPH and its realization of FVM
Pita et al. A fluid–structure interaction method for highly deformable solids
Zhang et al. An integrated coupling framework for highly nonlinear fluid-structure problems

Legal Events

Date Code Title Description
A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20150526

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20150727

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20150818

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20150928

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

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20151102

R150 Certificate of patent or registration of utility model

Ref document number: 5842992

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150