JP5842992B2 - シミュレーションプログラム、シミュレーション方法及びシミュレーション装置 - Google Patents
シミュレーションプログラム、シミュレーション方法及びシミュレーション装置 Download PDFInfo
- 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
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N13/00—Investigating surface or boundary effects, e.g. wetting power; Investigating diffusion effects; Analysing materials by determining surface, boundary, or diffusion effects
-
- 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
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/25—Design optimisation, verification or simulation using particle-based methods
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N13/00—Investigating surface or boundary effects, e.g. wetting power; Investigating diffusion effects; Analysing materials by determining surface, boundary, or diffusion effects
- G01N13/02—Investigating surface tension of liquids
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N19/00—Investigating materials by mechanical methods
- G01N19/02—Measuring coefficient of friction between materials
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-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
自由表面などの、変形し移動する境界を容易に扱いたいがために開発された粒子法であるが、離散化方法から、連続体が、粒子11から構成されるただの粒子群となってしまっているため、図1に示すように、連続体の境界面がどこかはっきりしなくなる。そのため、粒子法では表面張力などの境界を陽に扱う必要がある問題に対しては、統一的手法は出来ていない。
図2は、影響半径、およびカーネル関数の重ねあわせによる連続関数を示す図である。
まず、表面張力の算出に用いるcolor関数<Ci>を以下のように定義する。
はそれぞれ粒子11の位置ベクトル、粒子11の質量および密度であり、各粒子11に値が与えられている。
はSPH法で用いられるカーネル関数22の重ね合わせで、3次のスプライン関数23などがよく用いられる。
はカーネル関数22の値が非ゼロとなる球形領域の半径を表し、図2に示したように影響半径21と呼ばれる。
図3乃至図7に示すように、CSFモデルによる計算方法は、流体の塊が表面張力の影響を受けながら重力により落下し、平面に衝突する計算方法である。図3乃至図7に示した挙動の計算は、非特許文献2に開示されている計算方法を用いた。図3乃至図7においては、y=0が壁面を表すように鉛直方向をy軸に設定しており、y<0となる位置には流体31が侵入できないようになっている。
は0.001[m2/s]、重力加速度
コンピュータに、
格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
更新した該各粒子の速度に基づいて、前記流体の形状変化を算出し、
算出された前記形状変化を出力部に出力する
処理を実行させるシミュレーションプログラムであって、
前記該各粒子の速度の更新においては、
を用い、前記各粒子に係る力は、上記(式12)の右辺第2項で示される圧力勾配項、右辺第3項で示される粘性応力項、右辺第4項で示される重力項、右辺第5項は粒子間距離が近づきすぎないようにするための力であり、右辺第6項は表面張力であり、ここで、
また、下付き添え字は粒子の番号を表し、
:i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
:j番目の質量であり、
:粘性項を計算するために導入した定数であり、
であり
y=0の平面が固体の壁面の場合、
であり、
g:重力加速度
y:y方向の単位ベクトル
h:カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
W:SPH法で用いられるカーネル関数の重ね合わせ
は
の再標準化行列であり、
とおくと、
であり、
は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
まず、粒子i,jに関して、以下の特性量を定め、
さらに、各勾配を以下のよう計算し、
ここで、
とし、上記の量を用いて、
は以下のように決定する
というシミュレーションプログラムが提供される。
本発明を適用した実施の形態は、コンピュータに実行させるシミュレーションプログラム、シミュレーション方法及びシミュレーション装置であって、シミュレーションの対象である流体を粒子の集合として捉え、変分を基にした粒子間の相互作用力である表面張力を、所定の近似関数を用いて算出し、その算出した表面張力を用いることにより、流体の形状変化をシミュレーションする。
は流体が空気に触れている部分では1、それ以外の部分では0の値をとる。
はそれぞれ、図8に示すように、空気に触れている部分の表面張力係数、固体に触れている部分の表面張力係数である。
図9において、シミュレーション装置900は、処理部901、格納部902及び出力部903を備える。そして、シミュレーション装置900は、後述する近似間数を用い、流体等の連続体を粒子の集まりとして表現する粒子法による数値計算を行って表面張力を算出する。その計算結果を用いて連続体の形状変化を算出し、シミュレーション結果912を出力する。
出力部903は、処理部901が実行したシミュレーション結果912を出力する。
は定数であり、表面張力波の振動周期と合うように決められる。
は空間の次元、
はそれぞれ、液体が空気と触れているときの表面張力と液体が固体と触れているときの表面張力である。
は粒子の位置から決まる関数で、
は(式4)の
を幅
で、図10に示したように滑らかに近似した関数であり、
は粒子iと固体境界との距離である。
が面積の次元をもち、さらにSPH法では
は流体の内部の粒子では1に近く、表面近くになると値が小さくなり、影響半径h内に1つも粒子がなくなると、
という値となる。
は粘性応力テンソルで、流体112の粘性係数を
(定数)とすると、
は重力加速度、
はy方向の単位ベクトルである。
を定数とすれば、
上述の(式12)および(式14)は、粒子の位置を移動させ、y=0の平面の位置にある壁面を通過しようとする場合に、図12の上に示したように壁面をそのまま通過するのではなく、図12の下に示したように壁面に衝突した瞬間に止まるという現実に即した修正を表す。
は0.001[m2/s]、重力加速度
本発明を適用した第2の実施の形態では、上述の第1の実施の形態で、(式4)の積分を(式5)の近似関数のように近似した代わりに、(式4)の積分を下記(式31)の近似関数のように近似する。第1の実施の形態中の圧力「p1」を「p2」に置き換えることにより、第1の実施の形態中の説明を第2の実施の形態に適用できる。
上記(式31)中の
192は、
191を原点付近で滑らかにした近似関数で、
極小値の大きさが、
とほぼ同じであり、「0.1」となっている。
Claims (5)
- 流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーションプログラムにおいて、
コンピュータに、
格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
更新した該各粒子の速度に基づいて、前記流体の形状変化を算出し、
算出された前記形状変化を出力部に出力する
処理を実行させるシミュレーションプログラムであって、
前記該各粒子の速度の更新においては、
:i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
:j番目の質量であり、
:粘性項を計算するために導入した定数であり、
であり、
g:重力加速度
y:y方向の単位ベクトル
h:カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
W:SPH法で用いられるカーネル関数の重ね合わせ
は
であり、
は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
まず、粒子i,jに関して、以下の特性量を定め、
とし、上記の量を用いて、
は以下のように決定する
- 前記表面張力の算出は、空気との界面の表面エネルギーを算出した結果と前記空気以外との界面の表面エネルギーを算出した結果の和を算出することを特徴とする請求項1記載のシミュレーションプログラム。
- 前記形状変化の算出は、前記流体が他相を通り抜けようとする場合、前記表面張力の算出において用いられた前記粒子の位置及び速度を補正することを特徴とする請求項1記載のシミュレーションプログラム。
- 流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーション方法において、
コンピュータが、
格納部から、前記流体を形成する粒子の位置ベクトル及び速度ベクトルを取得し、
前記流体を粒子の集合で表現した流体モデルについて、取得した前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出し、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新し、
更新した該各粒子の速度に基づいて、前記流体の形状変化を算出し、
算出された前記形状変化を出力部に出力する
シミュレーション方法であって、
前記該各粒子の速度の更新においては、
:i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
:j番目の質量であり、
:粘性項を計算するために導入した定数であり、
であり、
g:重力加速度
y:y方向の単位ベクトル
h:カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
W:SPH法で用いられるカーネル関数の重ね合わせ
は
であり、
は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値であり、以下のように決定し、
まず、粒子i,jに関して、以下の特性量を定め、
とし、上記の量を用いて、
は以下のように決定する
- 流体の経時的な形状変化をシミュレーションする場合に、不自然な分布の速度分布が発生する挙動を抑制するシミュレーション装置において、
前記流体を形成する粒子の位置ベクトル及び速度ベクトルを格納する格納部と、
前記流体を粒子の集合で表現した流体モデルについて、前記位置ベクトル及び速度ベクトルを用いて、前記流体を形成する各粒子に係る力および該各粒子に係る表面張力を算出
する表面張力算出部と、
前記速度ベクトルに算出した各粒子に係る力及び各粒子に係る表面張力を加えて該各粒子の速度を更新する更新部と、
更新した該各粒子の速度に基づいて、前記流体の形状変化を算出する形状変化算出部と、
算出された前記形状変化を出力する出力部と、
を備えるシミュレーション装置であって、
前記該各粒子の速度の更新においては、
:i番目の粒子の位置ベクトル、速度ベクトル、密度、圧力であり、
:j番目の質量であり、
:粘性項を計算するために導入した定数であり、
であり、
g:重力加速度
y:y方向の単位ベクトル
h:カーネル関数の値が非ゼロとなる球形領域の半径である影響半径
W:SPH法で用いられるカーネル関数の重ね合わせ
は
であり、
は粒子i,j間での1次元リーマン問題を解いて求めた時空間の中間値であり、以下のよ
うに決定し、
まず、粒子i,jに関して、以下の特性量を定め、
とし、上記の量を用いて、
は以下のように決定する
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)
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)
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 | みずほ情報総研株式会社 | 流体解析装置、流体解析方法及び流体解析プログラム |
-
2012
- 2012-03-06 WO PCT/JP2012/055697 patent/WO2013132597A1/ja active Application Filing
- 2012-03-06 EP EP12870651.2A patent/EP2824598A1/en not_active Withdrawn
- 2012-03-06 JP JP2014503325A patent/JP5842992B2/ja active Active
-
2014
- 2014-08-25 US US14/467,703 patent/US9170185B2/en active Active
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 |