JP6808195B2 - 流体シミュレーションプログラム、流体シミュレーション装置および流体シミュレーション方法 - Google Patents

Description

本発明は、流体シミュレーションプログラムなどに関する。

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

これに対して、格子を用いない数値計算手法である粒子法では、移動境界の取り扱いに特別な処置を必要としない。したがって、近年、格子を用いない粒子法が広く用いられている（例えば、特許文献１〜３を参照）。なお、格子を用いない粒子法には、例えば、ＭＰＳ（Moving Particle Semi-implicit）法やＳＰＨ（Smoothed Particle Hydrodynamics）法などが挙げられる。

ところで、鋳造、鍛造、接合などの金属加工では、液体金属の中に冷えて固化した金属が混在したり、成長したりする。また、固化の過程で体積が変化したりする。また、溶融金属の表面では、空気中の酸素との結合により酸化膜が生じる。このような液体から固体への変化である凝固現象の取り扱いや、固体と気体との界面である自由表面の取り扱いは非常に重要である。自由表面の取り扱いの容易さ、並列性能面や固体との連成計算の容易さなどの利点により、鋳造、鍛造シミュレーションでは、粒子法の活躍が期待される分野であると考えられる。なお、液体と気体との界面に形成される膜のことを酸化膜や凝固膜と称される。

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

非特許文献１では、ＳＰＨ法を用いて、流体の方程式を以下のように離散化している。

式（１）は、質量保存則を表す。式（２）は、運動量保存則を表す。式（３）は、状態方程式を表す。式（４）は、エネルギー保存則を表す。ここで、ｘ_ｉは、粒子ｉの位置ベクトルである。ｖ_ｉは、粒子ｉの速度ベクトルである。ρ_ｉは、粒子ｉの密度である。ｍ_ｉは、粒子ｉの質量である。ｐ_ｉは、粒子ｉの圧力である。Ｕ_ｉは、粒子ｉの内部エネルギーである。Ｔ_ｉは、粒子ｉの温度である。ｘ_ｉｊは、粒子ｉと粒子ｊとの相対位置ベクトルであり、ｘ_ｉｊ＝ｘ_ｉ−ｘ_ｊである。ｖ_ｉｊは、粒子ｉと粒子ｊとの相対速度ベクトルであり、ｖ_ｉｊ＝ｖ_ｉ−ｖ_ｊである。ｋ_ｉは、粒子ｉの熱伝導率である。μ_ｉは、粒子ｉの粘性係数である。Ｐ_０は、ρ_０ｃ^２であり、ｃは音速である。ρ_ｓ，ｉは、粒子ｉの基準密度で、ρ_ｉ＝ρ_ｓ，ｉのとき、圧力が０になるような値である。Ｗは、カーネル関数（重み関数と呼ばれることもある）で、式（５）のスプライン関数がよく用いられている。

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

非特許文献１では、内部エネルギーＵ_ｉが低下し、温度Ｔ_ｉが融点以下に下がった際に、粘性係数μ_ｉが大きくなると定義することにより、式（２）の第３項の粒子同士の相対速度ｖ_ｉｊを打ち消す効果が大きくなり、粒子ｉが変形しにくくなることで固化を表現する。また、内部エネルギーＵ_ｉが低下し、温度Ｔ_ｉが融点以下に下がった際に、基準密度ρ_ｓ，ｉが大きくなるように定義することにより、圧力ｐ_ｉが低下し、式（２）の第２項の圧力の効果により周りの粒子が集まることで固化に伴う収縮を表現する。

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

特開２０１０−１０８１８３号公報 特開２０１３−１４７００２号公報 特開２０１４−１０４６９６号公報

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.

ここで、液体金属の表面は、瞬時に酸素と反応して酸化膜（凝固膜）を形成する。形成される凝固膜は、０．１ミリメートル（ｍｍ）以下であり、非常に薄い。シミュレーションを行う装置が液体金属の固化過程を計算する際、かかる薄い凝固膜を凝固現象として取り扱い、従来の手法を用いて計算できる。しかしながら、表面の凝固膜の解析に非常に小さい粒子を用いる必要があるため、演算量および計算時間の効率が悪いという問題がある。また、シミュレーションを行う装置が液体金属の固化過程を計算する際、表面の凝固膜の解析に凝固膜を解像できない粒子を用いて計算できるが、正確にシミュレーションすることができないという問題がある。

ここで、表面の凝固膜を解像できない１ｍｍ径の粒子を用いて液体アルミニウムを注湯した場合の注湯結果を、図９を参照して説明する。図９は、表面の凝固膜を解像できない粒子を用いた液体アルミニウムの注湯結果を示す図である。図９に示すように、左図は、実際に水を注湯した場合のｔ１時点とｔ２時点の水の注湯の実験結果を表している。中図は、実際に液体アルミニウムを注湯した場合のｔ１時点とｔ２時点の液体アルミニウムの注湯の実験結果を表している。右図は、液体アルミニウムを注湯する場合のｔ１時点とｔ２時点の液体アルミニウムの注湯のシミュレーション結果を表している。なお、シミュレーションは、凝固膜を表現できない１ｍｍ径の粒子を用いたＣｌｅａｒｙの手法である。

図９右図に示すように、凝固膜を表現できない１ｍｍ径の粒子を用いて液体アルミニウムの注湯計算を行うと、液体の飛び跳ね方や波の進行速度について、実際に液体アルミニウムを注湯した場合と比べて大きな違いがでる。したがって、実際の実験結果とシミュレーション結果が合致しない。つまり、表面の凝固膜の解析に凝固膜を表現できない粒子を用いて計算すると、正確にシミュレーションすることができない。

本発明は、１つの側面では、液体金属の固化過程を計算する際に、凝固膜を表現できない粒子径であっても、正確なシミュレーションを可能にするとともに、シミュレーションの演算量および計算時間を効率化することを目的とする。

１つの態様では、流体シミュレーションプログラムは、コンピュータに、物性値の第１粘性係数から、物質が凝固する際の表面膜の第２粘性係数であって前記第１粘性係数より大きい値の第２粘性係数を算出し、特定の時間の複数の粒子データについて、各粒子データから、各粒子データの粒子が前記物質の内部にある粒子を示す内部粒子であるか、前記物質の表面にある粒子を示す表面粒子であるかを判定し、前記粒子データの粒子が前記内部粒子であると判定した場合には、前記第１粘性係数を適用し、前記粒子データの粒子が前記表面粒子であると判定した場合には、前記第２粘性係数を適用し、適用した粘性係数を用いて、前記複数の粒子データそれぞれの次の時間の位置データと密度データとを算出し、前記複数の粒子データについて、算出した位置データと密度データとを時間と対応づけて記憶部に格納する、処理を実行させる。

１実施態様によれば、液体金属の固化過程を計算する際に、凝固膜を表現できない粒子径であっても、正確なシミュレーションを可能にするとともに、シミュレーションの演算量および計算時間を効率化することができる。

図１は、実施例に係る流体シミュレーション装置の構成を示す機能ブロック図である。 図２は、表面粒子の粘性係数を説明する図である。 図３は、粒子の近傍を表す空間領域の分類を示す図である。 図４Ａは、実施例に係る粒子判定の模式図（１）である。 図４Ｂは、実施例に係る粒子判定の模式図（２）である。 図５は、実施例に係る粒子判定処理のフローチャートの一例を示す図である。 図６は、実施例に係る流体シミュレーション処理のフローチャートの一例を示す図である。 図７は、実施例に係る流体シミュレーションを用いた液体アルミニウムの注湯結果を示す図である。 図８は、流体シミュレーションプログラムを実行するコンピュータの一例を示す図である。 図９は、表面の凝固膜を解像できない粒子を用いた液体アルミニウムの注湯結果を示す図である。

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

図１は、実施例に係る流体シミュレーション装置の構成を示す機能ブロック図である。図１に示す流体シミュレーション装置１は、格子を用いない粒子法を用いて液体金属の流動をシミュレーションする。かかる流体シミュレーション装置１は、粒子法により液体が固体へ変化する固化過程を計算する計算手法に対して、液体表面にある粒子（「表面粒子」という）の粘性係数を大きくして計算する。なお、実施例では、粒子法の１つであるＳＰＨ法により液体が固体へ変換する固化過程を計算する計算手法としてＣｌｅａｒｙの手法を用いるものとする。また、ここでいう液体金属とは、例えば、液体アルミニウムであるが、これに限定されるものではない。

流体シミュレーション装置１は、制御部１０と、記憶部２０とを有する。

制御部１０は、ＣＰＵ（Central Processing Unit）などの電子回路に対応する。そして、制御部１０は、各種の処理手順を規定したプログラムや制御データを格納するための内部メモリを有し、これらによって種々の処理を実行する。制御部１０は、シミュレーション受付部１１と、粘性係数算出部１２と、粒子判定部１３と、演算部１４と、出力部１５とを有する。

記憶部２０は、例えば、ＲＡＭ、フラッシュメモリ（Flash Memory）などの半導体メモリ素子、または、ハードディスク、光ディスクなどの記憶装置である。記憶部２０は、物性値情報２１と、逐次結果情報２２とを有する。

物性値情報２１は、シミュレーション対象を示す液体金属の物性値を示す情報である。物性値には、例えば、シミュレーション対象の液体金属について、基準密度、液体時および固体時のそれぞれの粘性係数、粒子径および表面凝固膜の大きさが挙げられる。粒子径は、粒子が球体である場合の直径である。粒子径は、一例として１ｍｍであるが、これに限定されず、液体金属の種類によって異なる。表面凝固膜の大きさは、表面が冷えたり酸化したりして凝固する場合の凝固膜の大きさである。表面凝固膜の大きさは、一例として０．００１ｍｍであるが、０．０１ｍｍであっても０．１ｍｍであっても良く、液体金属の種類によって異なる。

逐次結果情報２２は、流体シミュレーションにおいて、粒子群の時間ステップごとのシミュレーション結果の情報である。シミュレーション結果には、例えば、粒子を識別する識別子に対して、時間ステップごとの速度、密度、位置が設定される。

シミュレーション受付部１１は、流体シミュレーションの実行要求を受け付ける。例えば、シミュレーション受付部１１は、各種の物性値を含む流体シミュレーションの実行要求を受け付ける。そして、シミュレーション受付部１１は、実行要求に含まれた各種の物性値を物性値情報２１に設定する。なお、各種の物性値は、あらかじめ物性値情報２１に設定されていても良い。

粘性係数算出部１２は、物性値の１つである粘性係数から対象の液体金属の表面にある粒子の粘性係数を算出する。例えば、粘性係数算出部１２は、物性値情報２１から、液体時および固体時のそれぞれの粘性係数、粒子径および表面凝固膜の大きさを取得する。粘性係数算出部１２は、粒子が内部粒子である場合の粘性係数を液体時の粘性係数とする。ここでいう「内部粒子」とは、液体金属の内部、すなわち液体金属内にある粒子のことをいう。そして、粘性係数算出部１２は、粒子が表面粒子である場合の粘性係数を、粒子径のうち表面凝固膜の大きさだけを固体時の粘性係数とし、表面凝固膜の大きさを除外した大きさだけを液体時の粘性係数とするように算出する。ここでいう「表面粒子」とは、液体金属の表面にある粒子のことをいう。一例として、粘性係数算出部１２は、粘性係数を以下の式（６）のように算出する。なお、式（６）で用いられる記号は、以下のとおりである。ｄｘは、粒子径である。ｄ_ｍは、表面凝固膜の大きさである。μ_ｆは、液体時の粘性係数である。μ_solidは、固体時の粘性係数である。

ここで、表面粒子の粘性係数を、図２を参照して説明する。図２は、表面粒子の粘性係数を説明する図である。図２に示すように、液体表面にある粒子ｉ０（表面粒子）が表わされている。粒子ｉ０の粒子径は、ｄｘである。そして、粒子ｉ０の凝固膜の大きさは、ｄ_ｍであり、粒子ｉ０の液体部分の大きさは、ｄｘからｄ_ｍを差し引いた値である。したがって、粘性係数算出部１２は、表面粒子である粒子ｉ０の粘性係数を、式（６）の上部に示すように算出する。すなわち、粘性係数算出部１２は、粒子が表面にある場合で示すように、（ｄｘ−ｄ_ｍ）／ｄｘを液体時の粘性係数μ_ｆとし、ｄ_ｍ／ｄｘを固体時の粘性係数μ_solidとして、液体部分と固体部分の粘性係数を平均化する。

図１に戻って、粒子判定部１３は、特定の時間の複数の粒子データについて、各粒子データから、各粒子データの粒子が内部粒子か表面粒子かを判定する。例えば、粒子判定部１３は、特定の時間の複数の粒子データの中から、順次、判定対象の粒子として粒子データを選択する。粒子判定部１３は、判定対象の粒子の近傍に位置する複数の近傍粒子それぞれが、判定対象を中心とした複数の部分領域のうちいずれの部分領域に属するか否かを判定する。ここでいう近傍粒子とは、判定対象の粒子としてある粒子に着目したとき、当該粒子からの距離が所定半径より小さいところに存在する粒子のことをいう。所定半径は、判定に影響を与えるという意味で「影響半径」というものとする。そして、粒子判定部１３は、複数の部分領域の中に近傍粒子を１つも含まない部分領域が存在する場合には、判定対象の粒子が表面粒子であると判定する。粒子判定部１３は、複数の部分領域の中に近傍粒子を１つも含まない部分領域が存在しない場合には、判定対象の粒子が内部粒子であると判定する。すなわち、粒子判定部１３は、判定対象を中心とした影響半径内の近傍粒子がどの部分領域に属するかを判定し、近傍粒子が１つも所属しない部分領域がある場合には、判定対象の粒子が表面粒子であると判定する。

ここで、表面粒子の判定方法を、図３、図４Ａおよび図４Ｂに基づいて説明する。なお、図３、図４Ａおよび図４Ｂでは、空間が２次元の場合について説明する。

図３は、粒子の近傍を表す空間領域の分類を示す図である。図３に示すように、粒子ｉ０が判定対象の粒子であるとする。判定対象の粒子ｉ０を中心とした影響半径内の空間領域がｘ＋、ｘ−、ｙ＋、ｙ−の４つの部分領域に分類される。

図４Ａおよび図４Ｂは、実施例に係る粒子判定の模式図である。図４Ａに示すように、粒子判定部１３は、判定対象の粒子ｉ０を中心とした影響半径内にある近傍粒子が図３で示したいずれの部分領域に属するかを判定し、近傍粒子が１つも所属しない部分領域がある場合には、判定対象の粒子ｉ０が表面粒子であると判定する。

図４Ａ左図では、粒子判定部１３が、判定対象の粒子ｉａを中心とした影響半径内にある近傍粒子がいずれの部分領域に属するかを判定した結果である。ここでは、ｙ＋の部分領域に、１つも近傍粒子が所属しない。このため、粒子判定部１３は、判定対象の粒子ｉａが表面粒子であると判定する。

図４Ａ右図では、粒子判定部１３が、判定対象の粒子ｉｂを中心とした影響半径内にある近傍粒子がいずれの部分領域に属するかを判定した結果である。ここでは、全ての部分領域に近傍粒子が存在する。このため、粒子判定部１３は、判定対象の粒子ｉｂが内部粒子であると判定する。

ただし、かかる表面粒子の判定方法では、判定対象の粒子が表面粒子であるか否かの判定が座標軸の向きに依存するという問題がある。かかる問題について説明する。図４Ｂ左図に示すように、例えば、座標軸に沿って近傍粒子が並んでいる場合には、全ての部分領域に近傍粒子が存在するので、粒子判定部１３は、判定対象の粒子ｉｃを内部粒子として判定する。ところが、かかる近傍粒子の並びで、座標軸を４５度回転させ、再度表面粒子の判定を行うと、判定結果が異なる。

具体的には、図４Ｂ右図に示すように、座標軸を４５度回転させた結果、判定対象の粒子ｉｃを中心としたｘ＋、ｘ−、ｙ＋、ｙ−の４つの部分領域がｘ´−、ｘ´＋、ｙ´＋、ｙ´−の４つの部分領域に分類される。座標軸を４５度回転させた後に、粒子判定部１３は、再度、判定対象の粒子ｉｃの判定を行う。すると、近傍粒子が１つも所属しない部分領域ｙ´＋が存在するので、粒子判定部１３は、判定対象の粒子ｉｃを表面粒子として判定する。座標軸の４５度回転前後で判定対象の粒子ｉｃの判定結果が異なる。つまり、判定対象の粒子ｉｃが表面粒子であるか否かの判定が座標軸の向きに依存する場合がある。

そこで、粒子判定部１３は、回転前の部分領域（ｘ＋、ｘ−、ｙ＋、ｙ−）と回転後の部分領域（ｘ´＋、ｘ´−、ｙ´＋、ｙ´−）の全ての部分領域に近傍粒子が存在する場合に、判定対象の粒子が内部粒子であると判定する。また、粒子判定部１３は、いずれかの部分領域に近傍粒子が存在しない場合に、判定対象の粒子が表面粒子であると判定する。なお、図４Ｂでは、粒子判定部１３は、座標軸を４５度回転させてから判定対象の粒子を判定すると説明した。しかしながら、粒子判定部１３は、これに限定されず、判定対象の粒子に対する近傍粒子の位置（相対位置ベクトル）を４５度回転させてから判定対象の粒子を判定しても良い。

図１に戻って、演算部１４は、粒子ごとに適用される粘性係数を用いて、以下の方程式（７）〜（１０）により時間発展を計算することで、流動凝固の計算を行う。なお、ｘ_ｉは、粒子ｉの位置ベクトルである。ｖ_ｉは、粒子ｉの速度ベクトルである。ρ_ｉは、粒子ｉの密度である。ｍ_ｉは、粒子ｉの質量である。ｐ_ｉは、粒子ｉの圧力である。ｎは、時間ステップである。

方程式（７）は、粒子ｉが速度で移動することを表す式である。

方程式（８）は、運動量保存則を示す式（２）に対応する。

方程式（９）は、質量保存則を示す式（１）に対応する。

方程式（１０）は、粒子ｉが速度で移動し、位置を更新することを表す式である。

演算部１４は、時間ステップ０での、液体金属を表す粒子群の各粒子ｉの位置、速度、密度および質量を初期条件として入力する。すると、演算部１４は、時間ステップｎごとに、式（７）〜式（１０）を計算し、各粒子ｉの位置ｘ_ｉ、速度ｖ_ｉ、密度ρ_ｉを取得する。そして、演算部１４は、各粒子ｉについて、時間ステップｎごとの位置ｘ_ｉ、速度ｖ_ｉ、密度ρ_ｉを逐次結果情報２２に設定する。

式（８）のｐ_ｉｊ ^ｎは、時間ステップｎについて、粒子ｉの圧力と粒子ｊの圧力との平均である。すなわち、ｐ_ｉｊ ^ｎは、基準密度をρ_ｓとして、音速をｃとした状態方程式ｐ_ｉ ^ｎ＝ｃ^２（ρ_ｉ ^ｎ−ρ_ｓ）からｐ_ｉｊ ^ｎ＝（ｐ_ｉ ^ｎ＋ｐ_ｊ ^ｎ）／２と計算される。

式（８）の右辺第３項の（総乗）Π_ｉｊ ^{ｎ＋１,＊}は粘性応力を表し、以下の式（１１）で表わされる。粘性応力は、式（６）の粘性係数を適用して計算される。すなわち、演算部１４は、粒子データの粒子が表面粒子であると判定された場合には、式（６）の粒子が表面にある場合の粘性係数を適用する。これに対して、演算部１４は、粒子データの粒子が内部粒子であると判定された場合には、式（６）の表面以外にある場合の粘性係数を適用する。

式（１１）で表わされる粘性応力Π_ｉｊ ^{ｎ＋１,＊}は、未知数であるｖ_ｉ ^ｎ＋１，ｖ_ｊ ^ｎ＋１を含むが、式（８）は連立一次方程式となっているから、共役勾配法などの反復法により計算される。

演算部１４は、方程式（７）〜（１０）を繰り返し計算することにより、各粒子ｉの時間発展を行い、流動凝固のシミュレーションを実行する。

このように、演算部１４は、表面粒子であれば、内部粒子の粘性係数より大きい粘性係数に代えて計算することで、表面凝固膜による変形のし難さを表現することができる。

［粒子判定処理のフローチャート］
図５は、実施例に係る粒子判定処理のフローチャートの一例を示す図である。なお、図５では、空間が３次元の場合について説明する。

図５に示すように、粒子判定部１３は、各粒子の位置を入力する（ステップＳ１１）。そして、粒子判定部１３は、１つの粒子を選択し、選択した粒子を粒子ｉとする（ステップＳ１２）。粒子判定部１３は、粒子ｉの近傍粒子ｊに対して相対位置ベクトルｘ_ｉｊを計算する（ステップＳ１３）。例えば、粒子判定部１３は、粒子ｉの近傍粒子から１つの近傍粒子を選択し、粒子ｊとする。粒子判定部１３は、粒子ｉと粒子ｊとの相対位置ベクトルｘ_ｉｊを算出する。一例として、ｘ_ｉｊの成分は、（ｒ_ｘｉｊ，ｒ_ｙｉｊ，ｒ_ｚｉｊ）であるとする。

そして、粒子判定部１３は、相対位置ベクトルｘ_ｉｊを、ｘ軸中心に−４５度回転させ、ｘ´_ｉｊとする。粒子判定部１３は、相対位置ベクトルｘ_ｉｊを、ｙ軸中心に−４５度回転させ、ｘ´´_ｉｊとする。粒子判定部１３は、相対位置ベクトルｘ_ｉｊを、ｚ軸中心に−４５度回転させ、ｘ´´´_ｉｊとする（ステップＳ１４）。一例として、ｘ´_ｉｊの成分は、（ｒ´_ｘｉｊ，ｒ´_ｙｉｊ，ｒ´_ｚｉｊ）であるとする。ｘ´´_ｉｊの成分は、（ｒ´´_ｘｉｊ，ｒ´´_ｙｉｊ，ｒ´´_ｚｉｊ）であるとする。ｘ´´´_ｉｊの成分は、（ｒ´´´_ｘｉｊ，ｒ´´´_ｙｉｊ，ｒ´´´_ｚｉｊ）であるとする。

ここで、ｘ軸中心の−４５度回転は、角度をラジアン単位で表記して、以下の式（１２）で計算される。

ｙ軸中心の−４５度回転は、角度をラジアン単位で表記して、以下の式（１３）で計算される。

ｚ軸中心の−４５度回転は、角度をラジアン単位で表記して、以下の式（１４）で計算される。

そして、粒子判定部１３は、相対位置ベクトルｘ_ｉｊを用いて粒子ｊが属する部分領域がｘ＋、ｘ−、ｙ＋、ｙ−、ｚ＋、ｚ−のいずれであるかを判定する（ステップＳ１５）。これらの部分領域は、２次元の場合に図３で示したものを、３次元に拡張したものであある。すなわち、ｒ_ｘｉｊ／｜ｘ_ｉｊ｜≧ｃｏｓπ／４ならば、粒子ｊが部分領域ｘ＋に属すると判定される。ｒ_ｘｉｊ／｜ｘ_ｉｊ｜≦-ｃｏｓπ／４ならば、粒子ｊが部分領域ｘ−に属すると判定される。ｒ_ｙｉｊ／｜ｘ_ｉｊ｜≧ｃｏｓπ／４ならば、粒子ｊが部分領域ｙ＋に属すると判定される。ｒ_ｙｉｊ／｜ｘ_ｉｊ｜≦-ｃｏｓπ／４ならば、粒子ｊが部分領域ｙ−に属すると判定される。ｒ_ｚｉｊ／｜ｘ_ｉｊ｜≧ｃｏｓπ／４ならば、粒子ｊが部分領域ｚ＋に属すると判定される。ｒ_ｚｉｊ／｜ｘ_ｉｊ｜≦-ｃｏｓπ／４ならば、粒子ｊが部分領域ｚ−に属すると判定される。

そして、粒子判定部１３は、相対位置ベクトルｘ´_ｉｊを用いて粒子ｊが属する部分領域がｘ＋、ｘ−、ｙ＋、ｙ−、ｚ＋、ｚ−のいずれであるかを判定する（ステップＳ１６）。これらの部分領域は、ステップＳ１５で説明したものと同様である。

そして、粒子判定部１３は、相対位置ベクトルｘ´´_ｉｊを用いて粒子ｊが属する部分領域がｘ＋、ｘ−、ｙ＋、ｙ−、ｚ＋、ｚ−のいずれであるかを判定する（ステップＳ１７）。これらの部分領域は、ステップＳ１５で説明したものと同様である。

そして、粒子判定部１３は、相対位置ベクトルｘ´´´_ｉｊを用いて粒子ｊが属する部分領域がｘ＋、ｘ−、ｙ＋、ｙ−、ｚ＋、ｚ−のいずれであるかを判定する（ステップＳ１８）。これらの部分領域は、ステップＳ１５で説明したものと同様である。

そして、粒子判定部１３は、未処理の近傍粒子ｊが存在するか否かを判定する（ステップＳ１９）。未処理の近傍粒子ｊが存在すると判定した場合には（ステップＳ１９；Ｙｅｓ）、粒子判定部１３は、次の未処理の近傍粒子ｊを処理すべく、ステップＳ１３に移行する。

一方、未処理の近傍粒子ｊが存在しないと判定した場合には（ステップＳ１９；Ｎｏ）、粒子判定部１３は、ステップＳ１５〜ステップＳ１８の２４部分領域の中に、近傍粒子ｊを１つも含まない部分領域が存在するか否かを判定する（ステップＳ２０）。２４部分領域の中に、近傍粒子ｊを１つも含まない部分領域が存在すると判定した場合には（ステップＳ２０；Ｙｅｓ）、粒子判定部１３は、その粒子ｉを表面粒子と判定する（ステップＳ２２）。そして、粒子判定部１３は、ステップＳ２３に移行する。

一方、２４部分領域の中に、近傍粒子ｊを１つも含まない部分領域が存在しないと判定した場合には（ステップＳ２０；Ｎｏ）、粒子判定部１３は、その粒子ｉを内部粒子と判定する（ステップＳ２１）。そして、粒子判定部１３は、ステップＳ２３に移行する。

ステップＳ２３において、粒子判定部１３は、未処理の粒子ｉが存在するか否かを判定する(ステップＳ２３）。未処理の粒子ｉが存在すると判定した場合には（ステップＳ２３；Ｙｅｓ）、粒子判定部１３は、次の未処理の粒子ｉを処理すべく、ステップＳ１２に移行する。

一方、未処理の粒子ｉが存在しないと判定した場合には（ステップＳ２３；Ｎｏ）、粒子判定部１３は、粒子判定処理を終了する。

［流体シミュレーション処理のフローチャート］
図６は、実施例に係る流体シミュレーション処理のフローチャートの一例を示す図である。図６に示すように、シミュレーション受付部１１は、シミュレーション依頼を受け付けたか否かを判定する（ステップＳ３１）。シミュレーション依頼を受け付けていないと判定した場合には（ステップＳ３１；Ｎｏ）、シミュレーション受付部１１は、シミュレーション依頼を受け付けるまで、判定処理を繰り返す。

一方、シミュレーション依頼を受け付けたと判定した場合には（ステップＳ３１；Ｙｅｓ）、シミュレーション受付部１１は、シミュレーションする液体金属の物性値データおよび計算パラメータを入力する（ステップＳ３２）。例えば、シミュレーション受付部１１は、物性値データとして基準密度ρ０、液体時および固体時のそれぞれの粘性係数μ_ｆ，μ_{ｓｏｌｉｄ}を入力する。さらに、シミュレーション受付部１１は、物性値データとして、表面凝固膜の大きさｄ_ｍ、粒子径ｄｘを入力する。シミュレーション受付部１１は、計算パラメータとして音速ｃを入力する。

そして、シミュレーション受付部１１は、各粒子の初期条件を入力する（ステップＳ３３）。例えば、シミュレーション受付部１１は、時間ステップ０での、液体金属を表す粒子群の位置ベクトルｘ_ｉ ^０、速度ベクトルｖ_ｉ ^０、密度ρ_ｉ ^０、質量ｍ_ｉ ^０を入力する。なお、ｉやｊは、粒子群の中の粒子を識別するインデックスを示すものとする。

そして、演算部１４は、時間ステップｎ＋１／２での各粒子の位置を計算する（ステップＳ３４）。例えば、演算部１４は、現在の時間ステップｎでの位置ベクトルｘ_ｉ ^ｎ、速度ベクトルｖ_ｉ ^ｎから、式（７）を用いて、時間ステップｎ＋１／２での粒子の位置を計算する。

そして、粒子判定部１３は、時間ステップｎ＋１／２での各粒子の粒子判定処理を実行する（ステップＳ３５）。例えば、粒子判定部１３は、時間ステップｎ＋１／２での粒子群の位置ｘ_ｉ ^{ｎ＋１／２}を入力して、図５で示した粒子判定処理を実行する。そして、粒子判定部１３は、粒子群の中で表面粒子であると判定された粒子に対して表面粒子のフラグを立てる。

そして、粘性係数算出部１２は、時間ステップｎ＋１／２での表面粒子の粘性係数を計算する（ステップＳ３６）。例えば、粘性係数算出部１２は、表面粒子のフラグを入力して、式（６）を用いて、各粒子の粘性係数を計算する。具体的には、粘性係数算出部１２は、粒子ｉについて、表面粒子のフラグが立っている場合には、式（６）の上部（粒子が表面にある場合）の粘性係数を計算する。粘性係数算出部１２は、粒子ｉについて、表面粒子のフラグが立っていない場合には、式（６）の下部（表面以外にある場合）の粘性係数を計算する。なお、粘性係数算出部１２は、表面粒子の粘性係数を、各時間ステップで計算するが、これに限定されない。粘性係数算出部１２は、表面粒子の粘性係数を、最初の時間ステップで計算しておき、計算されたものを次の時間ステップで用いても良い。

続いて、演算部１４は、時間ステップｎ＋１での各粒子の速度を更新する（ステップＳ３７）。例えば、演算部１４は、粒子ｉについて、圧力ｐ_ｉ ^ｎを（音速をｃとした状態方程式）ｐ_ｉ ^ｎ＝ｃ^２（ρ_ｉ ^ｎ−ρ_０）から取得する。さらに、演算部１４は、粒子ｉについて、取得した圧力ｐ_ｉ ^ｎ、時刻ｎでの速度ベクトルｖ_ｉ ^ｎ、密度ρ_ｉ ^ｎ、時刻ｎ＋１／２での位置ベクトルｘ_ｉ ^{ｎ＋１／２}および計算された粘性係数を用いて、式（８）を共役勾配法により計算する。これにより、演算部１４は、時間ステップｎ＋１の各粒子の速度ベクトルｖ_ｉ ^ｎ＋１を計算する。

続いて、演算部１４は、時間ステップｎ＋１での各粒子の密度を更新する（ステップＳ３８）。例えば、演算部１４は、粒子ｉについて、式（９）を用いて、時間ステップｎ＋１での密度ρｉを計算する。なお、質量ｍ_ｉは、時間ステップｎのデータをそのまま時間ステップｎ＋１のデータとして用いれば良い。すなわち、演算部１４は、ｍ_ｉ ^ｎ＋１をｍ_ｉ ^ｎとして計算すれば良い。

続いて、演算部１４は、時間ステップｎ＋１での各粒子の位置を更新する（ステップＳ３９）。例えば、演算部１４は、粒子ｉについて、式（１０）を用いて、時間ステップｎ＋１での位置ベクトルｘ_ｉ ^ｎ＋１を計算する。

続いて、演算部１４は、時間ステップｎ＋１での各粒子のデータを出力する（ステップＳ４０）。例えば、演算部１４は、時間ステップｎ＋１での各粒子ｉの位置ベクトルｘ_ｉ ^ｎ＋１、速度ベクトルｖ_ｉ ^ｎ＋１、密度ρ_ｉ ^ｎ＋１を、例えば、逐次結果情報２２に設定する。

そして、演算部１４は、目的の時間ステップになったか否かを判定する（ステップＳ４１）。目的の時間ステップになっていないと判定した場合には（ステップＳ４１；Ｎｏ）、演算部１４は、目的の時間ステップまで流体シミュレーション処理を繰り返すべく、ステップＳ３４に移行する。一方、目的の時間ステップになったと判定した場合には（ステップＳ４１；Ｙｅｓ）、演算部１４は、流体シミュレーション処理を終了する。

［流体シミュレーションを用いた液体アルミニウムの注湯結果］
次に、実施例に係る流体シミュレーションを用いた液体アルミニウムの注湯結果を、図７を参照して説明する。図７は、実施例に係る流体シミュレーションを用いた液体アルミニウムの注湯結果を示す図である。図７に示すように、左図は、実際に水を注湯した場合のｔ１時点とｔ２時点の水の注湯の実験結果を表している。中図は、実際に液体アルミニウムを注湯した場合のｔ１時点とｔ２時点の液体アルミニウムの注湯の実験結果を表している。右図は、液体アルミニウムを注湯する場合のｔ１時点とｔ２時点の液体アルミニウムの注湯における実施例に係る流体シミュレーション結果を表している。

図７右図に示すように、流体シミュレーション結果は、図７中図の実際に液体アルミニウムを注湯した場合と一致する結果となっている。これにより、実施例に係る流体シミュレーション処理は、液体表面にある粒子の大きさを表面凝固膜より大きい液体時の大きさで計算しても、流体シミュレーションを正確に行うことができる。また、実施例に係る流体シミュレーション処理は、液体表面にある粒子の大きさを表面凝固膜より大きい液体時の大きさで計算することで、表面凝固膜を解像できる大きさで計算するより演算量を抑制できるとともに、計算時間を高速化することができる。

［実施例の効果］
上記実施例によれば、流体シミュレーション装置１は、物性値の第１粘性係数から、物質が凝固する際の表面膜の第２粘性係数であって第１粘性係数より大きい値の第２粘性係数を算出する。流体シミュレーション装置１は、特定の時間の複数の粒子データについて、各粒子データから、各粒子データの粒子が物質の内部にある粒子を示す内部粒子であるか、物質の表面にある粒子を示す表面粒子であるかを判定する。流体シミュレーション装置１は、粒子データの粒子が内部粒子であると判定した場合には、第１粘性係数を適用し、粒子データの粒子が表面粒子であると判定した場合には、第２粘性係数を適用する。そして、流体シミュレーション装置１は、適用した粘性係数を用いて、複数の粒子データそれぞれの次の時間の位置データと密度データとを算出する。そして、流体シミュレーション装置１は、複数の粒子データについて、算出した位置データと密度データとを時間と対応づけて記憶部２０に格納する。かかる構成によれば、流体シミュレーション装置１は、物資の凝固過程を計算する際に、表面膜を解像できない粒子径であっても、液体表面にある粒子の粘性係数を変えることで、正確なシミュレーションを可能にすることができる。また、流体シミュレーション装置１は、シミュレーションの演算量を抑制するとともに、計算時間を高速化することができる。

また、上記実施例によれば、流体シミュレーション装置１は、特定の大きさの粒子における液体時および固体時の粘性係数を示す第１粘性係数から、粒子の表面が凝固する際の表面膜の大きさに応じた粒子の第２粘性係数を算出する。かかる構成によれば、流体シミュレーション装置１は、大きさが同じ粒子であっても、液体時の場合と表面が凝固する場合とで粘性係数を変えることで、正確なシミュレーションが可能にすることができる。また、流体シミュレーション装置１は、シミュレーションの演算量を抑制するとともに、計算時間を高速化することができる。

また、上記実施例によれば、流体シミュレーション装置１は、判定対象の粒子の近傍に位置する複数の近傍粒子それぞれが、判定対象の粒子を中心とした複数の部分領域のうちいずれの部分領域に属するかを判定する。流体シミュレーション装置１は、複数の部分領域の中に近傍粒子を１つも含まない部分領域が存在する場合には、判定対象の粒子が表面粒子であると判定する。かかる構成によれば、流体シミュレーション装置１は、表面粒子を検出することで、表面粒子だけを物性値の粘性係数より大きい粘性係数に代えてシミュレーションでき、シミュレーションを正確に行うことが可能となる。

また、上記実施例によれば、流体シミュレーション装置１は、複数の部分領域の中に近傍粒子を１つも含まない部分領域が存在しない場合には、複数の部分領域それぞれを同じ方向に特定の角度だけ回転させる。そして、流体シミュレーション装置１は、回転後の複数の部分領域の中に近傍粒子を１つも含まない部分領域が存在する場合には、判定対象の粒子が表面粒子であると判定する。かかる構成によれば、流体シミュレーション装置１は、正確に表面粒子を検出することで、さらに、シミュレーションを正確に行うことが可能となる。

［その他］
なお、図示した流体シミュレーション装置の各構成要素は、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、流体シミュレーション装置１の分散・統合の具体的態様は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。例えば、演算部１４と出力部１５とを１つの部として統合しても良い。また、演算部１４を、各粒子の速度を演算する演算部と、各粒子の密度を演算する演算部と、各粒子の位置を演算する演算部とに分離しても良い。また、記憶部２０を流体シミュレーション装置１の外部装置としてネットワーク経由で接続するようにしても良い。

また、上記実施例で説明した各種の処理は、予め用意されたプログラムをパーソナルコンピュータやワークステーションなどのコンピュータで実行することによって実現することができる。そこで、以下では、図１に示した流体シミュレーション装置１と同様の機能を実現する流体シミュレーションプログラムを実行するコンピュータの一例を説明する。図８は、流体シミュレーションプログラムを実行するコンピュータの一例を示す図である。

図８に示すように、コンピュータ２００は、各種演算処理を実行するＣＰＵ２０３と、ユーザからのデータの入力を受け付ける入力装置２１５と、表示装置２０９を制御する表示制御部２０７とを有する。また、コンピュータ２００は、記憶媒体からプログラムなどを読取るドライブ装置２１３と、ネットワークを介して他のコンピュータとの間でデータの授受を行う通信制御部２１７とを有する。また、コンピュータ２００は、各種情報を一時記憶するメモリ２０１と、ＨＤＤ２０５を有する。そして、メモリ２０１、ＣＰＵ２０３、ＨＤＤ２０５、表示制御部２０７、ドライブ装置２１３、入力装置２１５、通信制御部２１７は、バス２１９で接続されている。

ドライブ装置２１３は、例えばリムーバブルディスク２１０用の装置である。ＨＤＤ２０５は、流体シミュレーションプログラム２０５ａおよび流体シミュレーション関連情報２０５ｂを記憶する。

ＣＰＵ２０３は、流体シミュレーションプログラム２０５ａを読み出して、メモリ２０１に展開し、プロセスとして実行する。かかるプロセスは、流体シミュレーション装置１の各機能部に対応する。流体シミュレーション関連情報２０５ｂは、物性値情報２１および逐次結果情報２２に対応する。そして、例えばリムーバブルディスク２１０が、流体シミュレーションプログラム２０５ａなどの各情報を記憶する。

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

１ 流体シミュレーション装置
１０ 制御部
１１ シミュレーション受付部
１２ 粘性係数算出部
１３ 粒子判定部
１４ 演算部
１５ 出力部
２０ 記憶部
２１ 物性値情報
２２ 逐次結果情報

Claims (5)

1. コンピュータに、
   物質の特定の大きさの粒子における液体時および固体時の粘性係数を示す第１粘性係数から、前記粒子の表面が凝固する際の表面膜の大きさに応じた前記粒子の第２粘性係数を算出し、
   特定の時間の複数の粒子について、各粒子が前記物質の内部にある粒子を示す内部粒子であるか、前記物質の表面にある粒子を示す表面粒子であるかを判定し、
   前記粒子が前記内部粒子であると判定した場合には、前記第１粘性係数を適用し、前記粒子が前記表面粒子であると判定した場合には、前記第２粘性係数を適用し、適用した粘性係数を用いて、前記複数の粒子それぞれの次の時間の位置データと密度データとを算出し、
   前記複数の粒子について、算出した位置データと密度データとを時間と対応づけて記憶部に格納する
   処理を実行させることを特徴とする流体シミュレーションプログラム。
2. 該判定する処理は、判定対象の粒子の近傍に位置する複数の近傍粒子それぞれが、前記判定対象の粒子を中心とした複数の部分領域のうちいずれの部分領域に属するかを判定し、前記複数の部分領域の中に前記近傍粒子を１つも含まない部分領域が存在する場合には、前記判定対象の粒子が前記表面粒子であると判定する
   処理を実行させることを特徴とする請求項１に記載の流体シミュレーションプログラム。
3. 該判定する処理は、前記複数の部分領域の中に前記近傍粒子を１つも含まない部分領域が存在しない場合には、前記複数の部分領域それぞれを同じ方向に特定の角度だけ回転させ、回転後の前記複数の部分領域の中に前記近傍粒子を１つも含まない部分領域が存在する場合には、前記判定対象の粒子が前記表面粒子であると判定する
   処理を実行させることを特徴とする請求項２に記載の流体シミュレーションプログラム。
4. 物質の特定の大きさの粒子における液体時および固体時の粘性係数を示す第１粘性係数から、前記粒子の表面が凝固する際の表面膜の大きさに応じた前記粒子の第２粘性係数を算出する第１の算出部と、
   特定の時間の複数の粒子について、各粒子が前記物質の内部にある粒子を示す内部粒子であるか、前記物質の表面にある粒子を示す表面粒子であるかを判定する判定部と、
   前記粒子が前記内部粒子であると判定した場合には、前記第１粘性係数を適用し、前記粒子が前記表面粒子であると判定した場合には、前記第２粘性係数を適用し、適用した粘性係数を用いて、前記複数の粒子それぞれの次の時間の位置データと密度データとを算出する第２の算出部と、
   前記複数の粒子について、算出した位置データと密度データとを時間と対応づけて記憶部に格納する格納部と、
   を有することを特徴とする流体シミュレーション装置。
5. コンピュータが、
   物質の特定の大きさの粒子における液体時および固体時の粘性係数を示す第１粘性係数から、前記粒子の表面が凝固する際の表面膜の大きさに応じた前記粒子の第２粘性係数を算出し、
   特定の時間の複数の粒子について、各粒子が前記物質の内部にある粒子を示す内部粒子であるか、前記物質の表面にある粒子を示す表面粒子であるかを判定し、
   前記粒子が前記内部粒子であると判定した場合には、前記第１粘性係数を適用し、前記粒子が前記表面粒子であると判定した場合には、前記第２粘性係数を適用し、適用した粘性係数を用いて、前記複数の粒子それぞれの次の時間の位置データと密度データとを算出し、
   前記複数の粒子について、算出した位置データと密度データとを時間と対応づけて記憶部に格納する
   各処理を実行することを特徴とする流体シミュレーション方法。

Images