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

Description

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

従来、流体や弾性体等の連続体の運動を解く数値計算の手法として、格子をベースにして微分方程式の近似解を求解する有限差分法や有限要素法、有限体積法などが用いられている。また、近年では、数値計算をＣＡＥ（Computer Aided Engineering）等の応用分野で活用するため、連続体の状態を解く数値計算の手法も発展し、流体と構造物とが相互作用する問題が解かれる。しかしながら、格子を用いる数値計算の手法では、自由表面等の界面の存在する問題や、流体と構造とを連成して解く流体構造連成解析問題などの移動境界が発生する場合には、連続体の取り扱いが複雑になるため、プログラム作成が困難になる場合がある。

ここで、格子を用いない数値計算の手法として、粒子法がある。粒子法とは、連続体の運動を有限の数の粒子の運動として解析する手法である。現在提案されている代表的な粒子法としては、ＳＰＨ（Smoothed Particles Hydrodynamics）法やＭＰＳ（Moving Particles Semi-implicit）法といったものがある。粒子法では、移動境界の取り扱いに特別な処置をせずに連続体の運動を解析することができる。そのため、粒子法は、近年、連続体の運動を解く数値計算の手法として広く用いられるようになってきている。

特に、鋳造・鍛造などの金属加工では、液体金属の中に冷えて固まった金属（固化した金属）が混じる、固化した金属が成長する、及び、固化過程において金属の体積が変化する等の複雑な過程となる。自由表面の扱いの簡単さや並列性能面、固体との連成計算が比較的容易である、などの利点により鋳造・鍛造シミュレーションは粒子法の活躍が期待される分野と考えられる。

鋳造過程の基本的なシミュレーション技術である、液体が冷えて固まる過程(固化過程)の計算手法として、Ｃｌｅａｒｙの方法が知られている。Ｃｌｅａｒｙの方法は、粒子法の一つであるＳＰＨ法を用いて、各液体粒子の内部エネルギーの時間発展を計算し、温度、密度、粘性係数を内部エネルギーの関数として計算を行う。すなわち、Ｃｌｅａｒｙの方法は、内部エネルギーが小さくなり温度が低くなると、液体の粘性係数を大きくして固化を表現し、液体の密度を大きくして固化に伴う体積低下を表現する。

Ｃｌｅａｒｙの手法は、流体の方程式をＳＰＨ法により以下の式（１）〜式（４）に示すように離散化する。

式（１）は質量保存則、式（２）は運動量保存則、式（３）は状態方程式、式（４）はエネルギー保存則を表す。ここで、ｘ_ｉ，ｖ_ｉ，ρ_ｉ，ｍ_ｉ，ｐ_ｉ，Ｕ_ｉは、それぞれ、粒子ｉの位置ベクトル、粒子ｉの速度ベクトル、粒子ｉの密度、粒子ｉの質量、粒子ｉの圧力、粒子ｉの内部エネルギーである。ｘ_ｉｊ，ｖ_ｉｊは、それぞれ、粒子ｉとｊの相対位置ベクトル、相対速度ベクトルであり、ｘ_ｉｊ＝ｘ_ｉ−ｘ_ｊ，ｖ_ｉｊ＝ｖ_ｉ−ｖ_ｊである。κ_ｉ，μ_ｉは、それぞれ、粒子ｉの熱伝導率、粒子ｉの粘性係数である。また、Ｐ_０＝ρ_０ｃ^２であり、ｃは音速である。ρ_ｓ，ｉは、粒子ｉの基準密度で、ρ_ｉ＝ρ_ｓ，ｉのとき圧力が０となる値である。

また、Ｗは、カーネル関数であり、例えば、Ｗとして、以下の式（５）が示すスプライン関数が用いられる。

ｈは、粒子間の影響半径であり、例えば、ｈとして、初期状態の平均粒子の間隔の２倍から３倍程度の値が用いられる。βは、カーネル関数の全空間積分量が１になるように調整された値であり、２次元の場合は、０．７πｈ^２と定められ、３次元の場合は、πｈ^３と定められる。

Ｃｌｅａｒｙの手法は、内部エネルギーが低下し、温度が融点以下に下がった際に粘性係数μ_ｉが大きくなることにより、式（２）の第三項の粒子同士の相対速度を打ち消す効果が大きくなり、変形しにくくなることで固化を表現する。また、Ｃｌｅａｒｙの手法は、基準密度ρ_ｓ，ｉを大きくすることで圧力を低下し、式（２）の第二項の効果により周りの粒子を集めることで、固化に伴う収縮を表現する。

式（１）〜式（４）を一般的な常微分方程式の数値解法である、オイラー法やリープフロッグ法などで、時間発展を計算することで、シミュレーションが可能となる。

Ｃｌｅａｒｙの手法では、固化過程における粘性係数の値が増大するため、計算上の時間刻みが非常に小さくなる。そのため、計算が終了する時刻までの計算回数が多大となる。それゆえ、Ｃｌｅａｒｙの手法では、多大な計算時間がかかる。

ここで、流体と剛体とが相互作用する問題の計算手法の一例として、液体の部分については液体の運動を解くための方程式を用い、固体の部分については剛体の運動方程式を用いて解く手法がある。かかる手法では、固体の部分を剛体の運動法的式を用いて解くため、Ｃｌｅａｒｙの手法よりも計算時間が短くてすむ。

Paul W. Cleary, "Extension of SPH to predict feeding freezing and defect creation in low pressure die casting", Applied Mathematical Modelling, 34(2010), pp. 3189-3201. Koshizuka, S., Nobe A. and Oka Y. "Numerical Analysis of Breaking Waves Using the Moving Particle Semi-implicit Method", Int. J. Numer. Meth. Fluids, 26, 751-769(1998)

しかしながら、固体の部分を剛体の運動方程式を用いて解く手法では、液体から固体が新たに生成される状況では、計算結果の精度が良好ではないという問題がある。液体の金属を型に流し込んで冷却する場合を例に挙げて説明する。この場合、液体の金属は、冷却の状況に応じて、複数の部分において固化が始まり、固化された複数の部分のそれぞれの固化された体積が時間に伴って大きくなる。そして、液体の金属は、全体が固化される。図１４は、従来の手法の問題点の一例を説明するための図である。図１４の例は、冷却によって固化された金属の固体の部分９０の粒子９０ａと、金属の固体の部分９１の粒子９１ａと、金属の液体の部分９２の粒子９２ａとが型に存在する場合を示す。この場合に、冷却によって、固体の部分９０の固化された体積が大きくなり、液体の部分９２が固化されて、固化された部分９２と固体の部分９０とが同一の固体となっても、上述した手法は、別々の固体として扱う。すなわち、上述した手法は、固化された部分９２の運動と固体の部分９０の運動とを別々に剛体の運動方程式を用いて解く。そのため、上述した手法では、本来１つの固体であるにも関わらず、複数の固体のそれぞれについて運動を解くので、計算結果の精度が良好ではない。

１つの側面では、液体の部分については液体の運動を解くための方程式を用い、固体の部分については剛体の運動方程式を用いて解く手法を用いる場合においても、精度良く計算を行うことを目的とする。

本願の開示するシミュレーションプログラムは、１つの態様において、コンピュータに、次の処理を実行させる。すなわち、シミュレーションプログラムは、コンピュータに、液体及び固体を含む連続体を複数の粒子で表した場合における複数の粒子のそれぞれのうち、第１の時間における液体の粒子について、第１の時間よりも後の第２の時間における状態を計算させる。シミュレーションプログラムは、コンピュータに、第２の時間における状態に基づいて、第１の時間における液体の粒子が、第２の時間に第１の固体の粒子となったか否かを判定させる。シミュレーションプログラムは、コンピュータに、第１の時間における液体の粒子が、第２の時間に第１の固体の粒子となったと判定された場合、次の処理を実行させる。すなわち、シミュレーションプログラムは、コンピュータに、第１の固体の粒子と、第１の固体の粒子から所定範囲内に存在する第２の固体の粒子を含む固体に属する全ての粒子とを同一の固体に属する粒子と定義させる。シミュレーションプログラムは、コンピュータに、剛体の運動方程式を用いて、同一の固体に属する粒子のそれぞれの状態を計算させる。

液体の部分については液体の運動を解くための方程式を用い、固体の部分については剛体の運動方程式を用いて解く手法を用いる場合においても、精度良く計算を行うことができる。

図１は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。 図２は、実施例に係るシミュレーション装置の機能構成の一例を示す図である。 図３は、金属モデルデータが示す金属モデルの一例を説明するための図である。 図４は、金属モデルデータが示す金属モデルの一例を説明するための図である。 図５は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。 図６は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。 図７は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。 図８は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。 図９は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。 図１０は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。 図１１は、実施例に係るシミュレーション処理の手順を示すフローチャートである。 図１２は、実施例に係る固体番号定義処理の手順を示すフローチャートである。 図１３は、シミュレーションプログラムを実行するコンピュータを示す図である。 図１４は、従来の手法の問題点の一例を説明するための図である。

以下に、本願の開示するシミュレーションプログラム、シミュレーション方法及びシミュレーション装置の実施例を図面に基づいて詳細に説明する。なお、実施例は開示の技術を限定するものではない。

［シミュレーション装置の構成］
実施例に係るシミュレーション装置について説明する。本実施例に係るシミュレーション装置は、金属の液体を型に流しこんで冷却させるシナリオを用いて、金属の液体及び固体を含む連続体の各粒子の状態を粒子法によりタイプステップｔ_ｔｓ毎に計算する。図１は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。図１に示すように、シミュレーション装置は、金属の液体２０を型２１に流しこんで冷却させて、金属の液体２０、及び、金属の液体２０が固化した金属の固体２１が混ざった場合の各粒子の状態を粒子法により計算する。

図２は、実施例に係るシミュレーション装置の機能構成の一例を示す図である。図２に示すように、シミュレーション装置１０は、入力部１１と、表示部１２と、記憶部１３と、制御部１４とを有する。

入力部１１は、制御部１４に情報を入力する。例えば、入力部１１は、ユーザから後述のシミュレーション処理を実行する指示であるシミュレーション実行指示を受け付けて、受け付けたシミュレーション実行指示を制御部１４に入力する。また、入力部１１は、ユーザから、初期状態における各粒子の初期値を受け付けて、受け付けた各粒子の初期値を制御部１４に入力する。ここで、初期状態における各粒子の初期値には、各粒子の位置、密度、速度、内部エネルギー、状態、固体番号が含まれる。ここで、粒子の状態には、内部エネルギーが所定値より高く、温度が融点より高い場合には、粒子が液体であるため、液体であることを示す情報が設定される。また、粒子の状態には、内部エネルギーが所定値以下であり、温度が凝固点より低い場合には、粒子が固体であるため、固体であることを示す情報が設定される。また、固体番号には、粒子が属する固体の識別番号が設定される。入力部１１のデバイスの一例としては、キーボードやマウスなどが挙げられる。

表示部１２は、各種の情報を表示する。例えば、表示部１２は、後述の表示制御部１４ｅの制御によりシミュレーション結果を表示する。表示部１２のデバイスの一例としては、液晶ディスプレイなどが挙げられる。

記憶部１３は、制御部１４で実行される各種プログラムを記憶する。また、記憶部１３は、金属モデルデータ１３ａを記憶する。金属モデルデータ１３ａは、金属の液体及び固体を含む連続体を複数の粒子として表した金属のモデルを示す。図３及び図４は、金属モデルデータが示す金属モデルの一例を説明するための図である。図３の例は、型２１に流し込まれた金属の液体２０の一部が固化されて、金属の液体２０と金属の固体２１とが混ざっている場合を示す。図３に示す部分２３の金属のモデルの詳細について説明する。図４は、図３に示す部分２３の金属のモデルの詳細を示す図である。図４の例に示すように、金属の液体２０及び金属の固体２２のモデルには、複数の粒子２５が含まれる。本実施例では、各粒子２５の位置、密度、速度、内部エネルギー、状態、固体番号がタイプステップｔ_ｔｓ毎に計算される。

図１に戻り、記憶部１３は、例えば、フラッシュメモリなどの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。なお、記憶部１３は、上記の種類の記憶装置に限定されるものではなく、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）であってもよい。

制御部１４は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。図２に示すように、制御部１４は、計算部１４ａと、更新部１４ｂと、判定部１４ｃと、定義部１４ｄと、表示制御部１４ｅとを有する。

計算部１４ａは、各種の情報を計算する。例えば、計算部１４ａは、タイムステップｔ_ｔｓ毎に、状態が液体の粒子である液体粒子に対して、流体の運動方程式を用いた時間発展の計算を行うことにより、液体粒子の位置、速度、密度、内部エネルギーを計算する。例えば、計算部１４ａは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての液体粒子に対して、上記の式（１）〜式（４）を用いて、タイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーを計算する。

また、計算部１４ａは、タイムステップｔ_ｔｓ毎に、状態が固体の粒子である固体粒子が属する固体に対して、剛体の運動方程式を用いた時間発展の計算を行うことにより、固体粒子の位置、速度、密度、内部エネルギーを計算する。例えば、計算部１４ａは、タイムステップ（ｔ_ｔｓ−１）の場合における固体粒子に対して、タイムステップｔ_ｔｓでの固体粒子の重心の並進運動と重心周りの回転運動とを計算する。これにより、計算部１４ａは、タイムステップｔ_ｔｓでの固体粒子が属する固体の時間発展を計算する。すなわち、計算部１４ａは、固体毎に、各固体に属する固体粒子にかかる力から、タイムステップｔ_ｔｓでの固体の重心の並進運動と重心周りの回転運動とを計算する。そして、計算部１４ａは、計算された固体の重心の並進運動と重心周りの回転運動とから、剛体の運動方程式を用いて、タイムステップｔ_ｔｓでの各固体に属する固体粒子の位置、速度、密度を計算する。また、計算部１４ａは、上記の式（４）を用いて、タイムステップｔ_ｔｓでの固体粒子の内部エネルギーを計算する。

計算部１４ａの一態様について説明する。例えば、計算部１４ａは、入力部１１からシミュレーション実行指示が入力された場合には、まず、タイムステップｔ_ｔｓの値を０に設定する。そして、計算部１４ａは、入力部１１から各粒子の初期値が入力されたか否かを判定する。初期値が入力された場合には、計算部１４ａは、タイムステップｔ_ｔｓの値を１だけインクリメントする。また、計算部１４ａは、タイムステップｔ_ｔｓの値がシミュレーションの最後のタイムステップＮ_Ｌ以下であると表示制御部１４ｅにより判定された場合にも、タイムステップｔ_ｔｓの値を１だけインクリメントする。

続いて、計算部１４ａは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての液体粒子に対して、上記の式（１）〜式（４）を用いて、タイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーを計算する。

また、計算部１４ａは、タイムステップ（ｔ_ｔｓ−１）の場合における固体粒子が属する全ての固体に対して、剛体の運動方程式を用いた時間発展の計算を行うことにより、次の処理を行う。すなわち、計算部１４ａは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての固体粒子のタイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーを計算する。

更新部１４ｂは、各種の情報を更新する。更新部１４ｂの一態様について説明する。例えば、更新部１４ｂは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての液体粒子の位置、速度、密度、内部エネルギーのそれぞれを、計算したタイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーのそれぞれに更新する。

また、更新部１４ｂは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての固体粒子の位置、速度、密度、内部エネルギーのそれぞれを、計算したタイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーのそれぞれに更新する。

そして、更新部１４ｂは、全ての粒子に対して、計算した内部エネルギーに基づいて、状態を更新する。例えば、更新部１４ｂは、計算した内部エネルギーが所定値より高い場合には、液体であることを示す情報を粒子の状態に設定し、状態を更新する。また、更新部１４ｂは、内部エネルギーが所定値以下である場合には、固体であることを示す情報を粒子の状態に設定し、状態を更新する。

続いて、更新部１４ｂは、全ての粒子の更新結果（位置、速度、密度、内部エネルギー、状態）をタイムステップｔ_ｔｓと対応付けて記憶部１３の所定領域に格納する。

続いて、判定部１４ｃは、各種の判定を行う。判定部１４ｃの一態様について説明する。例えば、判定部１４ｃは、タイムステップｔ_ｔｓにおける全ての粒子の更新前後の状態に基づいて、新たに固化した固体粒子（状態が液体から固体に変化した粒子）があるか否かを判定する。

新たに固化した固体粒子がある場合には、判定部１４ｃは、新たに固化した固体粒子を特定する。そして、判定部１４ｃは、特定した固体粒子の固体番号を未定義の状態にする。

そして、判定部１４ｃは、タイムステップｔ_ｔｓにおける全ての粒子の更新前後の状態に基づいて、新たに溶融した液体粒子（状態が固体から液体に変化した粒子）があるか否かを判定する。

新たに溶融した液体粒子がある場合には、判定部１４ｃは、新たに溶融した液体粒子を特定する。そして、判定部１４ｃは、特定した液体粒子がタイムステップ（ｔ_ｔｓ−１）において固体粒子であったときに属していた固体に属する全ての固体粒子の固体番号を未定義の状態にする。

そして、判定部１４ｃは、タイムステップｔ_ｔｓにおける全ての粒子の更新前後の状態に基づいて、再度、新たに固化した固体粒子があるか否かを判定するとともに、新たに溶融した液体粒子があるか否かを判定する。

定義部１４ｄは、各種の情報を定義する。定義部１４ｄの一態様について説明する。例えば、定義部１４ｄは、判定部１４ｃにより、新たに固化した固体粒子があると判定された場合、又は、新たに溶融した液体粒子があると判定された場合に、次の処理を行う。すなわち、定義部１４ｄは、固体番号が未定義の固体粒子のうち、未選択の固体粒子があるか否かを判定する。

未選択の固体粒子がある場合には、定義部１４ｄは、固体番号が未定義の未選択の固体粒子を１つ選択する。続いて、定義部１４ｄは、選択した固体粒子を中心とする半径ｈの球内に、固体に属する固体粒子があるか否かを判定する。なお、固体に属する固体粒子を判定する方法の一例としては、球内に存在する固体粒子の中に、固体粒子の固体番号に固体の識別番号が設定されているか否かを判定する方法がある。かかる方法では、識別番号が設定されている場合には、固体に属すると判定され、識別番号が設定されていない場合には、固体に属さないと判定される。また、半径ｈの値としては、任意の値を採用できるが、例えば、粒子法における粒子の影響半径を採用することができる。

図５〜図７は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。図５の例は、固体番号が未定義の固体粒子３０が定義部１４ｄにより選択された場合を示す。また、図５の例は、固体粒子３０と、固体３１に属する固体粒子３１ａのうち固体粒子３０と最も近い位置に存在する固体粒子３１ａとの距離がｒ_１であることを示す。また、図５の例は、固体粒子３０と、固体３２に属する固体粒子３２ａのうち固体粒子３０と最も近い位置に存在する固体粒子３２ａとの距離がｒ_２であることを示す。また、図５の例は、ｒ_１とｒ_２が、共に、半径ｈよりも長い場合を示す。図５の例の場合では、定義部１４ｄは、選択した固体粒子３０を中心とする半径ｈの球内に、固体に属する固体粒子がないと判定する。

図６の例は、固体番号が未定義の固体粒子３５が定義部１４ｄにより選択された場合を示す。また、図６の例は、固体粒子３５と、固体３６に属する固体粒子３６ａのうち固体粒子３５と最も近い位置に存在する固体粒子３６ａとの距離がｒ_３であることを示す。また、図６の例は、固体粒子３５と、固体３７に属する固体粒子３７ａのうち固体粒子３５と最も近い位置に存在する固体粒子３７ａとの距離がｒ_４であることを示す。また、図６の例は、ｒ_３が、半径ｈ以下であり、ｒ_４が、半径ｈよりも長い場合を示す。図６の例の場合では、定義部１４ｄは、選択した固体粒子３５を中心とする半径ｈの球内に、固体３６に属する固体粒子３６ａがあると判定する。

図７の例は、固体番号が未定義の固体粒子４０が定義部１４ｄにより選択された場合を示す。また、図７の例は、固体粒子４０と、固体４１に属する固体粒子４１ａのうち固体粒子４０と最も近い位置に存在する固体粒子４１ａとの距離がｒ_５であることを示す。また、図７の例は、固体粒子４０と、固体４２に属する固体粒子４２ａのうち固体粒子４０と最も近い位置に存在する固体粒子４２ａとの距離がｒ_６であることを示す。また、図７の例は、ｒ_５及びｒ_６が、共に、半径ｈ以下である場合を示す。図７の例の場合では、定義部１４ｄは、選択した固体粒子４０を中心とする半径ｈの球内に、複数の固体４１及び固体４２のそれぞれに属する固体粒子４１ａ及び固体粒子４２ａのそれぞれがあると判定する。

固体粒子がある場合には、定義部１４ｄは、あると判定された固体粒子が、複数の固体のそれぞれに属するか否かを判定する。例えば、図６の例の場合では、定義部１４ｄは、あると判定された固体粒子が、複数の固体のそれぞれに属さないと判定する。また、図７の例の場合では、定義部１４ｄは、あると判定された固体粒子が、複数の固体のそれぞれに属すると判定する。

定義部１４ｄは、あると判定された固体粒子が、複数の固体のそれぞれに属さない場合には、球内の固体粒子の固体番号の値を、選択した固体粒子の固体番号に設定する。例えば、図６の例の場合では、定義部１４ｄは、球内の固体粒子の固体番号の値「３６」を、選択した固体粒子３５の固体番号に設定する。これにより、１つの固体に属する粒子が１つの固体に属するように定義される。

一方、定義部１４ｄは、あると判定された固体粒子が、複数の固体のそれぞれに属する場合には、球内の複数の固体粒子のうち、選択した固体粒子に最も近い位置に存在する固体粒子の固体番号の値を、選択した固体粒子の固体番号に設定する。例えば、図７の例の場合では、定義部１４ｄは、球内の複数の固体粒子４１ａ、４２ａのうち、選択した固体粒子４０に最も近い位置に存在する固体粒子４１ａの固体番号の値「４１」を、選択した固体粒子４０の固体番号に設定する。これにより、１つの固体に属する粒子が１つの固体に属するように定義される。

そして、定義部１４ｄは、球内の複数の固体粒子のうち、選択した固体番号に設定された固体番号の値以外の値が固体番号に設定された固体粒子を特定する。例えば、図７の例の場合では、定義部１４ｄは、球内の複数の固体粒子４１ａ、４２ａのうち、選択した固体番号に設定された固体番号の値「４１」以外の値「４２」が固体番号に設定された固体粒子４２ａを特定する。

続いて、定義部１４ｄは、特定した固体粒子が属する固体に属する全ての固体粒子の固体番号の値を、選択した固体粒子の固体番号に設定された値に更新する。例えば、図７の例の場合では、定義部１４ｄは、特定した固体粒子４２ａが属する固体４２に属する全ての固体粒子４２ａ（３つの固体粒子４２ａ）の固体番号の値を、選択した固体粒子４０の固体番号に設定された値「４１」に更新する。これにより、選択した固体粒子を介して複数の固体のそれぞれに属する固体粒子が１つの固体に属するように定義される。

また、選択した固体粒子を中心とする半径ｈの球内に、固体に属する固体粒子がない場合には、定義部１４ｄは、選択した固体粒子の固体番号に、他の固体粒子の固体番号の値と重複しない値を設定する。例えば、図５の例の場合では、定義部１４ｄは、選択した固体粒子３０の固体番号に、他の固体粒子の固体番号の値と重複しない値「５０」を設定する。

そして、定義部１４ｄは、固体番号が未定義の固体粒子のうち未選択の固体粒子があるか否かを判定する上述した処理以降の処理を再び行う。これにより、定義部１４ｄは、固体番号が未定義の全ての固体粒子の固体番号を設定することができる。そして、定義部１４ｄは、固体番号が未定義の固体粒子のうち未選択の固体粒子がない場合には、タイムステップｔ_ｔｓに対応付けて、全ての粒子の固体番号を記憶部１３に格納する。

図８〜図１０は、実施例に係るシミュレーション装置が実行する処理の一例を説明するための図である。図８に示すように、定義部１４ｄは、上述した処理によって、半径ｈ内に存在する粒子６０を辿っていくことで、同一の個体に属する粒子６０を認識することができる。図９は、ある粒子６０の半径ｈ内に存在する粒子を連結させたときの粒子間の接続関係を示すリンクドリストを示す。

図１０の例は、タイムステップｔ_ｔｓにおいて新たに溶融した液体粒子７０がタイムステップ（ｔ_ｔｓ−１）において固体粒子であったときに属していた固体に属する全ての固体粒子７１の固体番号が、未定義の状態にされ、新たに固体番号が定義された場合を示す。図１０の例では、固体粒子７１のそれぞれの固体番号に対して、固体８１の固体番号「８１」または固体８２の固体番号「８２」が定義された場合を示す。ここで、図１０の例では、固体８１の粒子７１と、固体８２の粒子７１との最短の距離がｒ（＞ｈ）であるため、粒子７１は、１つの固体ではなく、２つの固体の何れかに属するように定義された場合を示す。

表示制御部１４ｅは、各種の情報の表示を制御する。表示制御部１４ｅの一態様について説明する。例えば、表示制御部１４ｅは、タイムステップｔ_ｔｓに対応付けて、全ての粒子の固体番号が定義部１４ｄにより記憶部１３に格納された場合に、タイムステップｔ_ｔｓの値が、シミュレーションの最後のタイムステップＮ_Ｌ以下であるか否かを判定する。タイムステップｔ_ｔｓの値が、シミュレーションの最後のタイムステップＮ_Ｌ以下でない場合には、表示制御部１４ｅは、次の処理を行う。すなわち、表示制御部１４ｅは、記憶部１３に格納されたタイムステップごとの全ての粒子の位置、速度、密度、内部エネルギー、状態、固体番号を全てのタイムステップについて取得する。そして、表示制御部１４ｅは、シミュレーション結果（全てのタイムステップにおける全ての粒子の位置、速度、密度、内部エネルギー、状態、固体番号）を表示するように表示部１２の表示を制御する。

制御部１４は、ＡＳＩＣ（Application Specific Integrated Circuit）やＦＰＧＡ（Field Programmable Gate Array）などのハードワイヤードロジックである。または、制御部１４は、ＣＰＵ（Central Processing Unit）やＭＰＵ（Micro Processing Unit）などにプログラムを実行させることにより実現される。

［処理の流れ］
次に、本実施例に係るシミュレーション装置１０の処理の流れを説明する。図１１は、実施例に係るシミュレーション処理の手順を示すフローチャートである。このシミュレーション処理の実行タイミングとしては様々なタイミングが考えられる。例えば、シミュレーション処理は、シミュレーション処理を実行するシミュレーション実行指示が入力部１１から入力された場合に、制御部１４により実行される。

図１１に示すように、計算部１４ａは、タイムステップｔ_ｔｓの値を０に設定する（Ｓ１０１）。そして、計算部１４ａは、入力部１１から各粒子の初期値が入力されたか否かを判定する（Ｓ１０２）。初期値が入力されていない場合（Ｓ１０２；Ｎｏ）には、計算部１４ａは、再び、Ｓ１０２での判定を行う。一方、初期値が入力された場合（Ｓ１０２；Ｙｅｓ）には、計算部１４ａは、タイムステップｔ_ｔｓの値を１だけインクリメントする（Ｓ１０３）。

続いて、計算部１４ａは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての液体粒子に対して、上記の式（１）〜式（４）を用いて、タイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーを計算する（Ｓ１０４）。続いて、更新部１４ｂは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての液体粒子の位置、速度、密度、内部エネルギーのそれぞれを、計算したタイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーのそれぞれに更新する（Ｓ１０５）。

そして、計算部１４ａは、タイムステップ（ｔ_ｔｓ−１）の場合における固体粒子が属する全ての固体に対して、剛体の運動方程式を用いた時間発展の計算を行うことにより、次の処理を行う。すなわち、計算部１４ａは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての固体粒子のタイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーを計算する（Ｓ１０６）。

続いて、更新部１４ｂは、タイムステップ（ｔ_ｔｓ−１）の場合における全ての固体粒子の位置、速度、密度、内部エネルギーのそれぞれを、計算したタイムステップｔ_ｔｓでの位置、速度、密度、内部エネルギーのそれぞれに更新する（Ｓ１０７）。

そして、更新部１４ｂは、全ての粒子に対して、計算した内部エネルギーに基づいて、状態を更新する（Ｓ１０８）。続いて、更新部１４ｂは、全ての粒子の更新結果（位置、速度、密度、内部エネルギー、状態）をタイムステップｔ_ｔｓと対応付けて記憶部１３の所定領域に格納する（Ｓ１０９）。

判定部１４ｃは、タイムステップｔ_ｔｓにおける全ての粒子の更新前後の状態に基づいて、新たに固化した固体粒子があるか否かを判定する（Ｓ１１０）。新たに固化した固体粒子がない場合（Ｓ１１０；Ｎｏ）には、後述するＳ１１３へ進む。

一方、新たに固化した固体粒子がある場合（Ｓ１１０；Ｙｅｓ）には、判定部１４ｃは、新たに固化した固体粒子を特定する（Ｓ１１１）。そして、判定部１４ｃは、特定した固体粒子の固体番号を未定義の状態にする（Ｓ１１２）。

そして、判定部１４ｃは、タイムステップｔ_ｔｓにおける全ての粒子の更新前後の状態に基づいて、新たに溶融した液体粒子があるか否かを判定する（Ｓ１１３）。新たに溶融した液体粒子がない場合（Ｓ１１３；Ｎｏ）には、後述するＳ１１６へ進む。一方、新たに溶融した液体粒子がある場合（Ｓ１１３；Ｙｅｓ）には、判定部１４ｃは、新たに溶融した液体粒子を特定する（Ｓ１１４）。そして、判定部１４ｃは、特定した液体粒子がタイムステップ（ｔ_ｔｓ−１）において固体粒子であったときに属していた固体に属する全ての固体粒子の固体番号を未定義の状態にする（Ｓ１１５）。

そして、判定部１４ｃは、タイムステップｔ_ｔｓにおける全ての粒子の更新前後の状態に基づいて、再度、新たに固化した固体粒子があるか否かを判定するとともに、新たに溶融した液体粒子があるか否かを判定する（Ｓ１１６）。

新たに固化した固体粒子がなく、かつ、新たに溶融した液体粒子がない場合（Ｓ１１６；Ｎｏ）には、Ｓ１１９に進む。一方、新たに固化した固体粒子がある場合、又は、新たに溶融した液体粒子がある場合（Ｓ１１６；Ｙｅｓ）に、次の処理を行う。すなわち、定義部１４ｄは、固体番号定義処理を実行する（Ｓ１１７）。そして、定義部１４ｄは、タイムステップｔ_ｔｓに対応付けて、全ての粒子の固体番号を記憶部１３に格納する（Ｓ１１８）。

そして、表示制御部１４ｅは、タイムステップｔ_ｔｓの値が、シミュレーションの最後のタイムステップＮ_Ｌ以下であるか否かを判定する（Ｓ１１９）。タイムステップｔ_ｔｓの値が、シミュレーションの最後のタイムステップＮ_Ｌ以下である場合（Ｓ１１９；Ｙｅｓ）には、Ｓ１０３に戻る。一方、タイムステップｔ_ｔｓの値が、シミュレーションの最後のタイムステップＮ_Ｌ以下でない場合（Ｓ１１９；Ｎｏ）には、表示制御部１４ｅは、次の処理を行う。すなわち、表示制御部１４ｅは、記憶部１３に格納されたタイムステップごとの全ての粒子の位置、速度、密度、内部エネルギー、状態、固体番号を全てのタイムステップについて取得する。そして、表示制御部１４ｅは、シミュレーション結果（全てのタイムステップにおける全ての粒子の位置、速度、密度、内部エネルギー、状態、固体番号）を表示するように表示部１２の表示を制御し（Ｓ１２０）、処理を終了する。

図１２は、実施例に係る固体番号定義処理の手順を示すフローチャートである。図１２に示すように、定義部１４ｄは、固体番号が未定義の固体粒子のうち、未選択の固体粒子があるか否かを判定する（Ｓ２０１）。

未選択の固体粒子がある場合（Ｓ２０１；Ｙｅｓ）には、定義部１４ｄは、固体番号が未定義の未選択の固体粒子を１つ選択する（Ｓ２０２）。続いて、定義部１４ｄは、選択した固体粒子を中心とする半径ｈの球内に、固体に属する固体粒子があるか否かを判定する（Ｓ２０３）。

固体粒子がある場合（Ｓ２０３；Ｙｅｓ）には、定義部１４ｄは、あると判定された固体粒子が、複数の固体のそれぞれに属するか否かを判定する（Ｓ２０４）。

あると判定された固体粒子が、複数の固体のそれぞれに属さない場合（Ｓ２０４；Ｎｏ）には、定義部１４ｄは、球内の固体粒子の固体番号の値を、選択した固体粒子の固体番号に設定し（Ｓ２０５）、Ｓ２０１に戻る。

一方、あると判定された固体粒子が、複数の固体のそれぞれに属する場合（Ｓ２０４；Ｙｅｓ）には、定義部１４ｄは、次の処理を行う。すなわち、定義部１４ｄは、球内の複数の固体粒子のうち、選択した固体粒子に最も近い位置に存在する固体粒子の固体番号の値を、選択した固体粒子の固体番号に設定する（Ｓ２０６）。

そして、定義部１４ｄは、球内の複数の固体粒子のうち、選択した固体番号に設定された固体番号の値以外の値が固体番号に設定された固体粒子を特定する（Ｓ２０７）。続いて、定義部１４ｄは、特定した固体粒子が属する固体に属する全ての固体粒子の固体番号の値を、選択した固体粒子の固体番号に設定された値に更新し（Ｓ２０８）、Ｓ２０１に戻る。

また、選択した固体粒子を中心とする半径ｈの球内に、固体に属する固体粒子がない場合（Ｓ２０３；Ｎｏ）には、定義部１４ｄは、選択した固体粒子の固体番号に、他の固体粒子の固体番号の値と重複しない値を設定し（Ｓ２０９）、Ｓ２０１に戻る。

固体番号が未定義の固体粒子のうち未選択の固体粒子がない場合（Ｓ２０１；Ｎｏ）には、定義部１４ｄは、処理結果を内部メモリに格納し、リターンする。

上述してきたように、本実施例に係るシミュレーション装置１０は、液体及び固体の状態となる金属を複数の粒子として表した場合の複数の粒子のそれぞれのうち、タイムステップ（ｔ_ｔｓ−１）における液体の粒子について、次の処理を行う。すなわち、シミュレーション装置１０は、タイムステップ（ｔ_ｔｓ−１）よりも後のタイムステップｔ_ｔｓにおける状態を計算する。そして、シミュレーション装置１０は、計算した状態に基づいて、タイムステップ（ｔ_ｔｓ−１）における液体の粒子が、タイムステップｔ_ｔｓに固体の粒子となったか否かを判定する。続いて、シミュレーション装置１０は、タイムステップ（ｔ_ｔｓ−１）における液体の粒子がタイムステップｔ_ｔｓに固体の粒子となったと判定した場合には、次の処理を行う。すなわち、シミュレーション装置１０は、かかる固体の粒子と、かかる固体の粒子から球内に存在する他の固体の粒子が属する固体に属する全ての粒子とを同一の固体に属する粒子と定義する。その後、シミュレーション装置１０は、同一の固体に対して剛体の運動方程式を用いて同一の固体に属する粒子のそれぞれの状態を計算する。それゆえ、本実施例に係るシミュレーション装置１０によれば、新たに液体から固体となった粒子を中心とする半径ｈの球内に存在する他の固体の粒子が属する固体に属する全ての粒子とを同一の固体に属する粒子と定義して、次の処理を行う。すなわち、本実施例に係るシミュレーション装置１０によれば、同一の固体に対して剛体の運動方程式を用いて同一の固体に属する粒子のそれぞれの状態を計算する。よって、本実施例に係るシミュレーション装置１０によれば、粒子を中心とする半径ｈの球内に存在するような、同一の個体に属する粒子を同一の個体として定義して、粒子のそれぞれの状態を計算する。したがって、本実施例に係るシミュレーション装置１０によれば、液体の部分については液体の運動を解くための方程式を用い、固体の部分については剛体の運動方程式を用いて解く手法を用いる場合においても、精度良く計算を行うことができる。

また、本実施例に係るシミュレーション装置１０によれば、固体の粒子については、剛体の運動方程式を用いて時間発展を計算するため、固体の粒子の粘性係数を増大させて時間発展を計算する場合よりも、粘性係数の増大を抑制することができる。そのため、本実施例に係るシミュレーション装置１０によれば、計算上の時間刻みを小さくせずにすむ。それゆえ、本実施例に係るシミュレーション装置１０によれば、計算が終了する時刻までの計算回数が多大となることを抑制することができる。したがって本実施例に係るシミュレーション装置１０によれば、計算時間が多大となることを抑制することができる。

また、シミュレーション装置１０によれば、新たに固化された固体の粒子から所定範囲内に存在する他の固体の粒子が属する固体が複数である場合には、新たに固化された固体の粒子と、複数の固体に属する全ての粒子とを同一の固体に属する粒子と定義する。したがって、新たに固化された固体粒子を介して複数の固体のそれぞれに属する全ての固体粒子が１つの固体に属するように定義することができる。

また、シミュレーション装置１０は、複数の粒子のそれぞれのうち、タイムステップ（ｔ_ｔｓ−１）における固体の粒子について、タイムステップ（ｔ_ｔｓ−１）よりも後のタイムステップｔ_ｔｓにおける状態を計算する。そして、シミュレーション装置１０は、タイムステップｔ_ｔｓにおける状態に基づいて、タイムステップ（ｔ_ｔｓ−１）における固体の粒子が、タイムステップｔ_ｔｓに液体の粒子となったか否かを判定する。続いて、シミュレーション装置１０は、タイムステップ（ｔ_ｔｓ−１）における固体の粒子がタイムステップｔ_ｔｓに液体の粒子となったと判定した場合には、次の処理を行う。すなわち、シミュレーション装置１０は、タイムステップ（ｔ_ｔｓ−１）における固体の粒子が属した固体に属する全ての粒子のそれぞれについて、粒子と、粒子を中心とする半径ｈの球内に存在する他の粒子とを同一の固体に属する粒子と定義する。そして、シミュレーション装置１０は、同一の固体に対して剛体の運動方程式を用いて同一の固体に属する粒子のそれぞれの状態を計算する。よって、本実施例に係るシミュレーション装置１０によれば、新たに溶融した粒子が属した固体に属する粒子を中心とする半径ｈの球内に存在するような、同一の個体に属する粒子を同一の個体として定義して、粒子のそれぞれの状態を計算することができる。したがって、本実施例に係るシミュレーション装置１０によれば、新たに粒子が溶融した場合には、新たに溶融した粒子が属した固体に属する固体の粒子について、属する固体の定義が行われる。そのため、本実施例に係るシミュレーション装置１０によれば、新たに溶融した粒子が属した固体に属した固体の粒子が属する固体の定義の精度が良好なものとなる。それゆえ、本実施例に係るシミュレーション装置１０によれば、精度良く粒子のそれぞれの状態を計算することができる。

さて、これまで開示の装置に関する実施例について説明したが、本発明は上述した実施例以外にも、種々の異なる形態にて実施されてよい。例えば、実施例などにおいて説明した各処理のうち、自動的に行われるものとして説明した処理の全部または一部を手動的に行うこともできる。

また、各種の負荷や使用状況などに応じて、実施例などにおいて説明した各処理の各ステップでの処理を任意に細かくわけたり、あるいはまとめたりすることができる。また、ステップを省略することもできる。

また、各種の負荷や使用状況などに応じて、実施例などにおいて説明した各処理の各ステップでの処理の順番を変更できる。

また、図示した装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、装置の分散・統合の具体的状態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。

［シミュレーションプログラム］
また、上記のシミュレーション装置１０のシミュレーション処理は、あらかじめ用意されたプログラムをパーソナルコンピュータやワークステーションなどのコンピュータシステムで実行することによって実現することもできる。そこで、以下では、図１３を用いて、上記のシミュレーション装置１０と同様の機能を有するシミュレーションプログラムを実行するコンピュータの一例を説明する。

図１３は、シミュレーションプログラムを実行するコンピュータを示す図である。図１３に示すように、コンピュータ３００は、ＣＰＵ（Central Processing Unit）３１０、ＲＯＭ（Read Only Memory）３２０、ＨＤＤ（Hard Disk Drive）３３０、ＲＡＭ（Random Access Memory）３４０を有する。これら３００〜３４０の各部は、バス３５０を介して接続される。

ＨＤＤ３３０には、上記の実施例で示す計算部１４ａ、更新部１４ｂ、判定部１４ｃ、定義部１４ｄ、表示制御部１４ｅと同様の機能を発揮するシミュレーションプログラム３３０ａが予め記憶される。なお、シミュレーションプログラム３３０ａについては、適宜分離しても良い。

そして、ＣＰＵ３１０が、シミュレーションプログラム３３０ａをＨＤＤ３３０から読み出して実行する。

そして、ＨＤＤ３３０には、図２の例に示す記憶部１３に記憶された金属モデルデータが設けられる。

そして、ＣＰＵ３１０は、ＨＤＤ３３０からデータを読み出してＲＡＭ３４０に格納する。さらに、ＣＰＵ３１０は、ＲＡＭ３４０に格納された各種のデータを用いて、シミュレーションプログラム３３０ａを実行する。なお、ＲＡＭ３４０に格納される各データは、常に全てのデータがＲＡＭ３４０に格納されなくともよく、全てのデータのうち処理に用いられるデータのみがＲＡＭ３４０に格納されれば良い。

なお、上記したシミュレーションプログラム３３０ａについては、必ずしも最初からＨＤＤ３３０に記憶させなくともよい。

例えば、コンピュータ３００に挿入されるフレキシブルディスク（ＦＤ）、ＣＤ−ＲＯＭ、ＤＶＤディスク、光磁気ディスク、ＩＣカードなどの「可搬用の物理媒体」にプログラムを記憶させておく。そして、コンピュータ３００がこれらからプログラムを読み出して実行するようにしてもよい。

さらには、公衆回線、インターネット、ＬＡＮ、ＷＡＮなどを介してコンピュータ３００に接続される「他のコンピュータ（またはサーバ）」などにプログラムを記憶させておく。そして、コンピュータ３００がこれらからプログラムを読み出して実行するようにしてもよい。

１０ シミュレーション装置
１３ 記憶部
１３ａ 金属モデルデータ
１４ 制御部
１４ａ 計算部
１４ｂ 更新部
１４ｃ 判定部
１４ｄ 定義部
１４ｅ 表示制御部

Claims (5)

1. コンピュータに、
   液体及び固体を含む連続体を複数の粒子で表した場合における前記複数の粒子のそれぞれのうち、第１の時間における液体の粒子について、第１の時間よりも後の第２の時間における状態を計算させ、
   前記第２の時間における状態に基づいて、前記第１の時間における液体の粒子が、前記第２の時間に第１の固体の粒子となったか否かを判定させ、
   前記第１の時間における液体の粒子が、前記第２の時間に前記第１の固体の粒子となったと判定された場合、予め各粒子について設定される、各粒子の位置、物理量、状態および固体番号に基づき、前記第１の固体の粒子と、前記第１の固体の粒子から所定範囲内に存在する第２の固体の粒子を含む固体に属する全ての粒子とに同一の固体番号を設定させ、
   剛体の運動方程式を用いて、前記同一の固体番号を設定された粒子のそれぞれの状態を計算させることを特徴とするシミュレーションプログラム。
2. 前記コンピュータに、
   前記第１の固体の粒子と、前記第２の固体の粒子を含む固体に属する全ての粒子とを同一の固体に属する粒子と定義させる場合、前記第１の固体の粒子から所定範囲内に存在する第２の固体の粒子が属する固体が複数存在する場合、前記第１の固体の粒子と、前記第２の固体の粒子が属する複数の固体に属する全ての粒子とを同一の固体に属する粒子と定義させることを特徴とする請求項１記載のシミュレーションプログラム。
3. 前記シミュレーションプログラムはさらに、
   前記コンピュータに、
   前記複数の粒子のそれぞれのうち、第３の時間における固体の粒子について、第３の時間よりも後の第４の時間における状態を計算させ、
   前記第４の時間における状態に基づいて、前記第３の時間における固体の粒子が、前記第４の時間に第１の液体の粒子となったかを判定させ、
   前記第３の時間における固体の粒子が、前記第４の時間に前記第１の液体の粒子となったと判定された場合、前記第３の時間において固体の粒子を含む固体に属する全ての粒子のそれぞれについて、前記第３の時間において固体の粒子を含む固体に属する全ての粒子と、前記第３の時間において固体の粒子を含む固体に属する全ての粒子のそれぞれの所定範囲内に存在する他の粒子とを同一の固体に属する粒子と定義させ、
   剛体の運動方程式を用いて、前記同一の固体に属する粒子のそれぞれの状態を計算させることを特徴とする請求項１または２記載のシミュレーションプログラム。
4. コンピュータが、
   液体及び固体を含む連続体を複数の粒子で表した場合における前記複数の粒子のそれぞれのうち、第１の時間における液体の粒子について、第１の時間よりも後の第２の時間における状態を計算し、
   前記第２の時間における状態に基づいて、前記第１の時間における液体の粒子が、前記第２の時間に第１の固体の粒子となったか否かを判定し、
   前記第１の時間における液体の粒子が、前記第２の時間に前記第１の固体の粒子となったと判定された場合、予め各粒子について設定される、各粒子の位置、物理量、状態および固体番号に基づき、前記第１の固体の粒子と、前記第１の固体の粒子から所定範囲内に存在する第２の固体の粒子を含む固体に属する全ての粒子とに同一の固体番号を設定し、
   剛体の運動方程式を用いて、前記同一の固体番号を設定された粒子のそれぞれの状態を計算することを特徴とするシミュレーション方法。
5. 液体及び固体を含む連続体を複数の粒子で表した場合における前記複数の粒子のそれぞれのうち、第１の時間における液体の粒子について、第１の時間よりも後の第２の時間における状態を計算するとともに、剛体の運動方程式を用いて、同一の固体番号を設定された粒子のそれぞれの状態を計算する計算部と、
   前記第２の時間における状態に基づいて、前記第１の時間における液体の粒子が、前記第２の時間に第１の固体の粒子となったかを判定する判定部と、
   前記第１の時間における液体の粒子が、前記第２の時間に前記第１の固体の粒子となったと判定された場合、予め各粒子について設定される、各粒子の位置、物理量、状態および固体番号に基づき、前記第１の固体の粒子と、前記第１の固体の粒子から所定範囲内に存在する第２の固体の粒子を含む固体に属する全ての粒子とに同一の固体番号を設定する定義部と
   を有することを特徴とするシミュレーション装置。

Images