特許 7276768 液体シミュレーション装置、液体シミュレーション方法および液体シミュレーションプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2023-05-10

(45)【発行日】2023-05-18

(54)【発明の名称】液体シミュレーション装置、液体シミュレーション方法および液体シミュレーションプログラム

(51)【国際特許分類】

   G01N 11/00 20060101AFI20230511BHJP

   B22D 46/00 20060101ALI20230511BHJP

【ＦＩ】

G01N11/00 C

B22D46/00

【請求項の数】 4

(21)【出願番号】P 2019127244

(22)【出願日】2019-07-08

(65)【公開番号】P2021012142

(43)【公開日】2021-02-04

【審査請求日】2022-03-31

(73)【特許権者】

【識別番号】000005223

【氏名又は名称】富士通株式会社

(73)【特許権者】

【識別番号】391002487

【氏名又は名称】学校法人大同学園

(74)【代理人】

【識別番号】110002147

【氏名又は名称】弁理士法人酒井国際特許事務所

(72)【発明者】

【氏名】風間 正喜

(72)【発明者】

【氏名】諏訪 多聞

(72)【発明者】

【氏名】畠中 耕平

(72)【発明者】

【氏名】前田 安郭

【審査官】山口 剛

(56)【参考文献】

【文献】特開２０１８－１８５７４３（ＪＰ，Ａ）

【文献】特開２０１７－０２１５９９（ＪＰ，Ａ）

【文献】特開２０１４－１４６３０２（ＪＰ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｎ １１／００ － １１／１６

Ｇ０１Ｎ ３３／２０５

Ｂ２２Ｄ ４６／００

Ｇ１６Ｚ ９９／００

(57)【特許請求の範囲】

【請求項1】

粒子法を用いて溶解金属の鋳造状況をシミュレーションする際に、前記鋳造状況のシミュレーション対象である溶解金属内の各粒子が前記溶解金属の内部に位置する内部粒子であるか前記溶解金属の表面に位置する表面粒子であるかを判定する判定部と、
前記内部粒子と判定された各粒子に対しては、酸化膜が生じていない状態の粘度である第１粘性係数を算出し、前記表面粒子と判定された各粒子に対しては、前記酸化膜ができてから指数関数に比例して減少する値の関係式を用いて前記酸化膜の回復タイミングとして算出された所定時間が経過するまでは前記第１粘性係数を算出するとともに、前記所定時間の経過後は前記第１粘性係数よりも大きい値の第２粘性係数を算出する算出部と、
前記各粒子に対して算出された前記第１粘性係数または前記第２粘性係数を用いて、前記鋳造状況のシミュレーションを実行する実行部と
を有することを特徴とする液体シミュレーション装置。

【請求項2】

前記判定部は、前記表面粒子と判定された各粒子が、前記表面粒子かつ鋳造先の壁に接している粒子、または、前記表面粒子かつ粒子圧力が閾値以上の粒子か否かを判定し、
前記算出部は、前記表面粒子と判定された各粒子のうち、前記表面粒子かつ鋳造先の壁に接している粒子および前記表面粒子かつ粒子圧力が閾値以上の粒子については、前記所定時間が経過するまでは前記第１粘性係数を算出するとともに、前記所定時間の経過後は前記第２粘性係数を算出する、
ことを特徴とする請求項１に記載の液体シミュレーション装置。

【請求項3】

コンピュータが、
粒子法を用いて溶解金属の鋳造状況をシミュレーションする際に、前記鋳造状況のシミュレーション対象である溶解金属内の各粒子が前記溶解金属の内部に位置する内部粒子であるか前記溶解金属の表面に位置する表面粒子であるかを判定し、
前記内部粒子と判定された各粒子に対しては、酸化膜が生じていない状態の粘度である第１粘性係数を算出し、前記表面粒子と判定された各粒子に対しては、前記酸化膜ができてから指数関数に比例して減少する値の関係式を用いて前記酸化膜の回復タイミングとして算出された所定時間が経過するまでは前記第１粘性係数を算出するとともに、前記所定時間の経過後は前記第１粘性係数よりも大きい値の第２粘性係数を算出し、
前記各粒子に対して算出された前記第１粘性係数または前記第２粘性係数を用いて、前記鋳造状況のシミュレーションを実行する
処理を実行することを特徴とする液体シミュレーション方法。

【請求項4】

コンピュータに、
粒子法を用いて溶解金属の鋳造状況をシミュレーションする際に、前記鋳造状況のシミュレーション対象である溶解金属内の各粒子が前記溶解金属の内部に位置する内部粒子であるか前記溶解金属の表面に位置する表面粒子であるかを判定し、
前記内部粒子と判定された各粒子に対しては、酸化膜が生じていない状態の粘度である第１粘性係数を算出し、前記表面粒子と判定された各粒子に対しては、前記酸化膜ができてから指数関数に比例して減少する値の関係式を用いて前記酸化膜の回復タイミングとして算出された所定時間が経過するまでは前記第１粘性係数を算出するとともに、前記所定時間の経過後は前記第１粘性係数よりも大きい値の第２粘性係数を算出し、
前記各粒子に対して算出された前記第１粘性係数または前記第２粘性係数を用いて、前記鋳造状況のシミュレーションを実行する
処理を実行させることを特徴とする液体シミュレーションプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

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

【背景技術】

【0002】

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

【0003】

そして、鋳造過程の基本的なシミュレーション技術である、液体が冷えて固まる過程（固化過程）の計算手法には、Ｃｌｅａｒｙの手法が知られている。Ｃｌｅａｒｙの手法では、粒子法の１つであるＳＰＨ法を用いて、各液体粒子の内部エネルギーの時間発展を計算し、温度、密度、粘性係数を内部エネルギーの関数として計算する。

【0004】

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

【0005】

近年では、粒子の位置により、粘性係数を変更してシミュレーションを行う技術が知られている。例えば、シミュレーションを行う際に、各粒子が物質の内部にある内部粒子か、物質の表面にある表面粒子かを判定する。そして、表面粒子だけ係数が高い粘性係数を用いてシミュレーションすることで、粒子データの粒子系を大きく設定しても正確なシミュレーションを行うとともに、処理の高速化を図る。

【先行技術文献】

【特許文献】

【0006】

【文献】特開２０１８－１８５７４３号公報

【発明の概要】

【発明が解決しようとする課題】

【0007】

しかしながら、上記技術では、酸化膜の破れた後、酸化膜の回復による物質の一時的な流動停止と物質の蛇行がシミュレーションできないので、流体シミュレーションの精度が低下する。

【0008】

例えば、アルミニウム合金などの液体金属の流体シミュレーションを行う場合、ラドルから金型に注入された液体金属は、金型の壁にぶつかった部分の酸化膜が破れ、それ以外の部分は酸化膜ができるので、酸化膜が破れた壁側に蛇行するように流れる。その後、破れた酸化膜が回復し、一時的な流動停止が発生する。ところが、上記技術では、酸化膜は表面に瞬時に形成されると仮定されていて、酸化膜の破れと回復の位置やタイミングを表現することができないので、流体シミュレーションに誤差が生じ、正確性が低下する。

【0009】

一つの側面では、流体シミュレーションの精度低下を抑制することができる液体シミュレーション装置、液体シミュレーション方法および液体シミュレーションプログラムを提供することを目的とする。

【課題を解決するための手段】

【0010】

第１の案では、液体シミュレーション装置は、鋳造状況のシミュレーション対象である物質内の各粒子について、予め定めた所定条件を満たすか否かを判定する判定部を有する。液体シミュレーション装置は、前記所定条件を満たさない粒子に対しては、酸化膜が生じていない状態の粘度である第１粘性係数を算出し、前記所定条件を満たす粒子に対しては、所定時間が経過するまでは前記第１粘性係数を算出するとともに、前記所定時間の経過後は前記第１粘性係数よりも大きい値の第２粘性係数を算出する算出部を有する。液体シミュレーション装置は、前記各粒子に対して算出された前記第１粘性係数または前記第２粘性係数を用いて、前記鋳造状況のシミュレーションを実行する実行部を有する。

【発明の効果】

【0011】

一実施形態によれば、流体シミュレーションの精度低下を抑制することができる。

【図面の簡単な説明】

【0012】

【図1】図１は、実施例１にかかる流体シミュレーション装置の機能構成を示す図である。

【図2】図２は、流体シミュレーションでは表現できない流動挙動を説明する図である。

【図3】図３は、表面粒子の粘性係数を説明する図である。

【図4】図４は、粒子の近傍を表す空間領域の分類を示す図である。

【図5A】図５Ａは、実施例に係る粒子判定の模式図である。

【図5B】図５Ｂは、実施例に係る粒子判定の模式図である。

【図6】図６は、フラグを利用した粘性係数の算出を説明する図である。

【図7】図７は、全体的な処理の流れを示すフローチャートである。

【図8】図８は、酸化膜判定処理の流れを示すフローチャートである。

【図9】図９は、シミュレーション結果を説明する図である。

【図10】図１０は、シミュレーション結果を説明する図である。

【図11】図１１は、ハードウェア構成例を説明する図である。

【発明を実施するための形態】

【0013】

以下に、本願の開示する液体シミュレーション装置、液体シミュレーション方法および液体シミュレーションプログラムの実施例を図面に基づいて詳細に説明する。なお、この実施例によりこの発明が限定されるものではない。また、各実施例は、矛盾のない範囲内で適宜組み合わせることができる。

【実施例1】

【0014】

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

【0015】

ここで、以前から行われている一般的な流体シミュレーションでは、注湯実験で現れる流動挙動を表現できないことがある。図２は、流体シミュレーションでは表現できない流動挙動を説明する図である。図２に示すように、流体シミュレーション装置１０は、ラドルから金型に液体アルミニウムを流し込み、その時の流動挙動をシミュレーションする。ところが、流体シミュレーションでは、図２の（ａ）や（ｂ）に示す流動挙動を表現できない。

【0016】

図２の（ａ）に示す注湯実験は、金型確度が０度、傾動速度が１６度／秒で液体アルミニウムを流し込んだ実験である。また、図２の（ｂ）に示す注湯実験は、金型確度が４０度、傾動速度が２１度／秒で液体アルミニウムを流し込んだ実験である。図２の（ａ）に示す注湯実験では、鋳造金型に流し込まれた液体アルミニウムが薄くくびれてから液滴上に膨らむ形状の流動挙動が発生する。また、図２の（ｂ）に示す注湯実験では、鋳造金型に流し込まれた液体アルミニウムが蛇行して最初の反射で細く飛び出る流動挙動が発生する。

【0017】

これらは、物理的には酸化膜が発生したことで起こると考えられる。したがって、酸化膜の破れることを考慮していない一般的な流体シミュレーションでは、注湯実験とは異なる流動挙動となる。そこで、実施例１にかかる流体シミュレーション装置１０は、酸化膜の破れを計算で扱うことにより、図２の（ａ）や（ｂ）に示す注湯実験と同等の流動挙動のシミュレーションを実現する。

【0018】

図１に戻り、流体シミュレーション装置１０は、通信部１１、記憶部１２、制御部２０を有する。通信部１１は、他の装置との間の通信を制御する処理部であり、例えば通信インタフェースなどである。例えば、通信部１１は、管理者の端末からシミュレーションの開始指示や各種データを受信し、管理者の端末にシミュレーションのシミュレーション結果を送信する。

【0019】

記憶部１２は、データや制御部２０が実行する各種プログラムを記憶する記憶装置の一例であり、例えばメモリやハードディスクなどである。この記憶部１２は、物性値情報ＤＢ１３と逐次結果情報ＤＢ１４を記憶する。

【0020】

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

【0021】

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

【0022】

制御部２０は、流体シミュレーション装置１０全体を司る処理部であり、例えばプロセッサなどである。この制御部２０は、実行部３０とフラグ判定部４０を有する。なお、実行部３０とフラグ判定部４０は、プロセッサが有する電子回路の一例やプロセッサが実行するプロセスの一例である。

【0023】

実行部３０は、シミュレーション受付部３１、粘性係数算出部３２、粒子判定部３３、演算部３４を有し、流体シミュレーションを実行する処理部である。例えば、実行部３０は、特開２０１８－１８５７４３号公報と同様の手法により、流体シミュレーションを実行する。

【0024】

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

【0025】

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

【0026】

【数1】

【0027】

ここで、表面粒子の粘性係数を、図３を参照して説明する。図３は、表面粒子の粘性係数を説明する図である。図３に示すように、液体表面にある粒子ｉ０（表面粒子）が表わされている。粒子ｉ０の粒子径は、ｄｘである。そして、粒子ｉ０の凝固膜の大きさは、ｄ_ｍであり、粒子ｉ０の液体部分の大きさは、ｄｘからｄ_ｍを差し引いた値である。したがって、粘性係数算出部３２は、表面粒子である粒子ｉ０の粘性係数を、式（１－１）示すように算出する。すなわち、粘性係数算出部３２は、粒子が表面にある場合で示すように、（ｄｘ－ｄ_ｍ）／ｄｘを液体時の粘性係数μ_ｆとし、ｄ_ｍ／ｄｘを固体時の粘性係数μ_solidとして、液体部分と固体部分の粘性係数を平均化する。なお、粘性係数算出部３２は、後述するフラグ判定部４０により算出されるフラグの値によっても、式（１－１）または式（１－２）の選択を動的に変更することができる。

【0028】

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

【0029】

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

【0030】

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

【0031】

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

【0032】

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

【0033】

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

【0034】

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

【0035】

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

【0036】

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

【0037】

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

【0038】

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

【数2】

【0039】

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

【数3】

【0040】

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

【数4】

【0041】

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

【数5】

【0042】

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

【0043】

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

【0044】

式（３）の右辺第３項の（総乗）Π_ｉｊ ^{ｎ＋１,＊}は粘性応力を表し、以下の式（６）で表わされる。粘性応力は、式（１）の粘性係数を適用して計算される。すなわち、演算部３４は、粒子データの粒子が表面粒子であると判定された場合には、式（１－１）の粒子が表面にある場合の粘性係数を適用する。これに対して、演算部３４は、粒子データの粒子が内部粒子であると判定された場合には、式（１－２）の表面以外にある場合の粘性係数を適用する。式（６）で表わされる粘性応力Π_ｉｊ ^{ｎ＋１,＊}は、未知数であるｖ_ｉ ^ｎ＋１，ｖ_ｊ ^ｎ＋１を含むが、式（３）は連立一次方程式となっているから、共役勾配法などの反復法により計算される。なお、演算部３４は、後述するフラグ判定部４０により算出されるフラグの値によっても、式（１－１）または式（１－２）の選択を動的に変更することができる。

【数6】

【0045】

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

【0046】

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

【0047】

フラグ判定部４０は、鋳造状況のシミュレーション対象である物質内の各粒子について、予め定めた所定条件を満たすか否かを判定し、判定結果に応じて、演算部３４の処理を制御する処理部である。具体的には、フラグ判定部４０は、所定条件を満たさない粒子に対しては、酸化膜が生じていない状態の粘度である第１粘性係数を用いて、シミュレーションを実行させる。一方、フラグ判定部４０は、所定条件を満たす粒子に対しては、所定時間が経過するまでは第１粘性係数を用いてシミュレーションを実行し、所定時間経過後は第１粘性係数よりも大きい値の第２粘性係数を用いてシミュレーションを実行させる。

【0048】

例えば、フラグ判定部４０は、粒子判定部３３により内部粒子と判定された各粒子については、フラグ「１．０」を設定して演算部３４に通知することで、式（１－２）を用いた内部粒子のシミュレーションを実行させる。

【0049】

また、フラグ判定部４０は、粒子判定部３３により表面粒子と判定された各粒子については、一定時間の間、フラグの値「１．０」を設定して演算部３４等に通知することで、表面粒子と判定されたにも関わらず、式（１－１）ではなく式（１－２）を用いた、内部粒子と同様のシミュレーションを実行させる。そして、フラグ判定部４０は、時間経過とともにフラグの値を減少させ、一定時間経過後にフラグの値が「０．５」以下となった場合、フラグの値を演算部３４等に通知して、式（１－１）への切り替えを実行させ、表面粒子としてのシミュレーションを実行させる。

【0050】

つまり、表面粒子については、酸化膜が出来るが、金型の壁などに衝突することで圧力がかかって酸化膜が破れ後に回復する。したがって、酸化膜が破れてから回復する期間は、表面粒子を内部粒子として扱い、酸化膜が回復した以降は、表面粒子を表面粒子として扱う。

【0051】

より詳細には、フラグ判定部４０は、式（７）を用いて、酸化膜の回復タイミングを算出する。ここで、ｔは時間であり、ｋは式（８）により定義される。また、Ｔは、酸化膜が破れてから有効になるまでの時間（半減期）であり、酸素濃度などシミュレーション環境により依存することから、予め算出して設定する値である。

【0052】

【数7】

【数8】

【0053】

ここで、式（７）によるフラグの値（ｆ（ｔ））の値の変化と粘性係数の関係について説明する。図６は、フラグを利用した粘性係数の算出を説明する図である。表面粒子以外の粒子（内部粒子）では、図６の（ａ）に示すように、フラグの値が「１．０」となるので、式（１－２）を用いたシミュレーションが実行される。

【0054】

一方、表面粒子では、図６の（ｂ）に示すように、時間ｔ１が経過するまでフラグの値「０．５」より大きな値となるので、式（１－２）を用いたシミュレーションが実行される。そして、時間ｔ１が経過すると、フラグの値「０．５」以下となるので、式（１－１）を用いたシミュレーションが実行される。

【0055】

このように、フラグ判定部４０は、酸化膜が出来てから、指数関数に比例して減少する関係式を用いて、酸化膜の回復タイミングを算出することで、表面粒子に適用する粘性係数を動的に変化させて、図２の（ａ）や（ｂ）に示す注湯実験と同等の流動挙動のシミュレーションを実現する。

【0056】

なお、フラグ判定部４０は、より厳密に判定するために、所定の条件を変更することもできる。例えば、フラグ判定部４０は、実行部３０の処理や粒子法等により取得される情報を用いて、表面粒子かつ壁に接している粒子、または、表面粒子かつ圧力の絶対値が閾値以上の粒子に対して、フラグの値が「０．５」以下になるまでは、式（１－２）を用いたシミュレーションを実行させることもできる。

【0057】

［全体的な処理の流れ］
図７は、全体的な処理の流れを示すフローチャートである。図７に示すように、シミュレーション受付部３１は、シミュレーション依頼を受け付けたか否かを判定する（Ｓ１０１）。シミュレーション依頼を受け付けていないと判定した場合には（Ｓ１０１：Ｎｏ）、シミュレーション受付部３１は、シミュレーション依頼を受け付けるまで、判定処理を繰り返す。

【0058】

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

【0059】

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

【0060】

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

【0061】

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

【0062】

その後、フラグ判定部４０は、酸化膜判定処理を実行する（Ｓ１０６）。例えば、フラグ判定部４０は、表面粒子と判定された各粒子については、式（８）によりフラグの値を計算し、計算結果に基づいて演算部３４に対して各種制御を実行する。

【0063】

そして、粘性係数算出部３２は、表面粒子の粘性係数を計算する（Ｓ１０７）。例えば、粘性係数算出部３２は、表面粒子のフラグが立っている粒子については、式（１－１）を用いて粘性係数を算出し、表面粒子のフラグが立っていない粒子については、式（１－２）を用いて粘性係数を算出する。ここで、粘性係数算出部３２は、表面粒子であっても、フラグ判定部４０から通知されたフラグの値が「０．５」より大きい間は、式（１－１）ではなく式（１－２）を用いて粘性係数を算出することもできる。

【0064】

その後、演算部３４は、各粒子の速度、密度、位置を更新し（Ｓ１０８）、各粒子の更新データを逐次結果情報ＤＢ１４に出力する（Ｓ１０９）。

【0065】

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

【0066】

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

【0067】

また、演算部３４は、粒子ｉについて、式（６）を用いて、時間ステップｎ＋１での位置ベクトルｘ_ｉ ^ｎ＋１を計算する。そして、演算部３４は、時間ステップｎ＋１での各粒子ｉの位置ベクトルｘ_ｉ ^ｎ＋１、速度ベクトルｖ_ｉ ^ｎ＋１、密度ρ_ｉ ^ｎ＋１を、逐次結果情報ＤＢ１４に設定する。

【0068】

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

【0069】

［酸化膜判定処理の流れ］
図８は、酸化膜判定処理の流れを示すフローチャートである。なお、この処理は、図７のＳ１０６で実行される。

【0070】

図８に示すように、フラグ判定部４０は、各粒子について酸化膜の破れ条件に該当するか否かを判定する（Ｓ２０１）。例えば、フラグ判定部４０は、表面粒子に該当するか否かを判定する。または、フラグ判定部４０は、表面粒子かつ壁に接している粒子に該当するか否か、または、表面粒子かつ圧力の絶対値が閾値以上である粒子に該当するか否かを判定する。

【0071】

そして、フラグ判定部４０は、酸化膜の破れ条件に該当する場合（Ｓ２０１：Ｙｅｓ）、式（８）に示したフラグ関数を初期化する（Ｓ２０２）。続いて、フラグ判定部４０は、シミュレーションが開始されてからの時間経過によってフラグ関数を更新する（Ｓ２０３）。

【0072】

その後、フラグ判定部４０は、フラグ関数の値が閾値以下の場合（Ｓ２０４：Ｙｅｓ）、表面粒子の条件で粘性係数を計算させる（Ｓ２０５）。一方、フラグ判定部４０は、フラグ関数の値が閾値より大きい（Ｓ２０４：Ｎｏ）、または、Ｓ２０１において酸化膜の破れ条件に該当しない場合（Ｓ２０１：Ｎｏ）、表面粒子以外の条件で粘性係数を計算させる（Ｓ２０６）。

【0073】

なお、フラグ判定部４０は、フラグ関数の値が閾値以下の場合、粘性係数算出部３２に式（１－２）を用いて粘性係数の算出を指示したり、演算部３４に式（１－２）を用いて各種演算を実行するように指示したりすることができる。

【0074】

［効果］
上述したように、流体シミュレーション装置１０は、表面粒子と内部粒子の判定に加えて、酸化膜の破れおよび回復のタイミングを計算により求めることで、表面粒子であっても酸化膜の破れと回復を考慮した粘性係数を用いたシミュレーションを実行することができる。

【0075】

図９と図１０は、シミュレーション結果を説明する図である。図９と図１０は、シミュレーションを開始してからの各時刻のシミュレーション結果と、各時刻に対応する注湯実験とを対比して図示している。なお、酸化膜が破れる圧力を４０００［Pas］、半減期を０．０５５［秒］とする。

【0076】

図９に示すように、注湯実験では、液体アルミニウムの注湯を開始してから２．０［秒］で酸化膜が発生して流動が停止する。その後、２．２［秒］後に、液体アルミニウムが金型の壁に衝突して酸化膜が破れて、蛇行した流動挙動が発生する。そして、２．４［秒］後では、再度酸化膜が発生して流動が弱まった後、２．６秒後には金型全体へ広がるように流し込まれる。

【0077】

一方、シミュレーションでも同様に、液体アルミニウムの注湯を開始してから２．０［秒］で酸化膜が発生して流動が停止した後、２．２［秒］後に酸化膜が破れて蛇行した流動挙動が発生することをシミュレーションできている。その後も、シミュレーションにおいても、注湯実験と同様の流動挙動を示すことができた。

【0078】

また、図１０に示すように、注湯実験では、注湯開始から２．８［秒］後に、液体アルミニウムが注湯箇所から反対の壁に到達し、生成された酸化膜を先頭に、注湯箇所まで戻ってくる。一方、シミュレーションにおいても、実験と同様、液体アルミニウムが注湯箇所から反対の壁に到達し、酸化膜を先頭に注湯箇所まで折り返して戻ってくることがシミュレーションできている。

【0079】

したがって、流体シミュレーション装置１０は、酸化膜の回復による物質の一時的な流動停止と物質の蛇行を正確にシミュレーションすることができるので、シミュレーションの精度低下を抑制することができる。

【実施例2】

【0080】

さて、これまで本発明の実施例について説明したが、本発明は上述した実施例以外にも、種々の異なる形態にて実施されてよいものである。

【0081】

［データや数値等］
上記実施例で用いたデータ例、数値例、閾値、表示例、秒などの時間例等は、あくまで一例であり、任意に変更することができる。

【0082】

［システム］
上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。

【0083】

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

【0084】

さらに、各装置にて行なわれる各処理機能は、その全部または任意の一部が、ＣＰＵおよび当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。

【0085】

［ハードウェア］
図１１は、ハードウェア構成例を説明する図である。図１１に示すように、流体シミュレーション装置１０は、通信装置１０ａ、ＨＤＤ（Hard Disk Drive）１０ｂ、メモリ１０ｃ、プロセッサ１０ｄを有する。また、図１１に示した各部は、バス等で相互に接続される。

【0086】

通信装置１０ａは、ネットワークインタフェースカードなどであり、他の装置との通信を行う。ＨＤＤ１０ｂは、図１に示した機能を動作させるプログラムやＤＢを記憶する。

【0087】

プロセッサ１０ｄは、図１に示した各処理部と同様の処理を実行するプログラムをＨＤＤ１０ｂ等から読み出してメモリ１０ｃに展開することで、図１等で説明した各機能を実行するプロセスを動作させる。例えば、このプロセスは、流体シミュレーション装置１０が有する各処理部と同様の機能を実行する。具体的には、プロセッサ１０ｄは、実行部３０とフラグ判定部４０等と同様の機能を有するプログラムをＨＤＤ１０ｂ等から読み出す。そして、プロセッサ１０ｄは、実行部３０とフラグ判定部４０等と同様の処理を実行するプロセスを実行する。

【0088】

このように、流体シミュレーション装置１０は、プログラムを読み出して実行することで流体シミュレーション方法を実行する情報処理装置として動作する。また、流体シミュレーション装置１０は、媒体読取装置によって記録媒体から上記プログラムを読み出し、読み出された上記プログラムを実行することで上記した実施例と同様の機能を実現することもできる。なお、この他の実施例でいうプログラムは、流体シミュレーション装置１０によって実行されることに限定されるものではない。例えば、他のコンピュータまたはサーバがプログラムを実行する場合や、これらが協働してプログラムを実行するような場合にも、本発明を同様に適用することができる。

【符号の説明】

【0089】

１０ 流体シミュレーション装置
１１ 通信部
１２ 記憶部
１３ 物性値情報ＤＢ
１４ 逐次結果情報ＤＢ
２０ 制御部
３０ 実行部
３１ シミュレーション受付部
３２ 粘性係数算出部
３３ 粒子判定部
３４ 演算部
４０ フラグ判定部

【図1】

【図2】

【図3】

【図4】

【図5A】

【図5B】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】