JP6273170B2

JP6273170B2 - 流動解析装置、流動解析方法、及びコンピュータプログラム

Info

Publication number: JP6273170B2
Application number: JP2014115452A
Authority: JP
Inventors: 知和中川; 和英関山; 紗矢香山田
Original assignee: Kobe Steel Ltd
Current assignee: Kobe Steel Ltd
Priority date: 2014-06-04
Filing date: 2014-06-04
Publication date: 2018-01-31
Anticipated expiration: 2034-06-04
Also published as: JP2015230530A; DE102015108687A1; US20150355007A1

Description

本発明は、粒子法を用いて流動性を有する物体の流動解析を行う流動解析装置、流動解析方法、及びコンピュータプログラムに関する。

物体の運動解析の手法は、固定された計算格子を用いる差分法及び有限要素法等と、計算格子を用いずに物体を粒子の集まりとして取り扱う粒子法とに大別される。粒子法では、固定された計算格子を用いる方法では不可能であった、物体の大変形の解析が可能であり、飛沫が生じるような流体の解析に有用とされる。

特許文献１には、粒子法による流体解析方法が開示されている。特許文献１に開示された方法では、あるタイムステップにおける粒子間の作用力から粒子の仮配置を決定し、この仮配置を用いて、粒子の数密度を一定にする計算により、次回タイムステップの圧力を計算し、この圧力を用いて算出した圧力勾配から各粒子の修正後の速度を算出し、この速度を用いて各粒子の位置を修正する。

特開平２００８−１１１６７５号公報

従来の粒子法による流動解析では、解析結果の粒子の分布に密な領域と粗な領域とが発生する。上述した特許文献１では、粒子の数密度が、粒子間距離に応じた重みの和として定義されている。かかる粒子の数密度を一定にするという条件は隣り合う粒子間の距離を一定にすることを保証するものではなく、特許文献１に記載された方法は、仮配置における粒子の分布の粗密を解消するものではない。粒子の分布に粗密が生じてしまうと、微圧縮性流体のような実質的に粗密が生じることがない物体の解析では精度が低下してしまうという問題がある。

本発明は斯かる事情に鑑みてなされたものであり、その主たる目的は、上記課題を解決することができる流動解析装置、流動解析方法、及びコンピュータプログラムを提供することにある。

上述した課題を解決するために、本発明の一の態様の流動解析装置は、粒子から構成された解析対象における粒子の位置の時間変化を推定する流動解析装置であって、前回のタイムステップにおける各粒子の位置及び各粒子の物理的状態に関する粒子情報に基づき、各粒子の速度を決定し、決定された速度に基づいて今回のタイムステップにおける各粒子の位置を決定する位置決定手段と、前記位置決定手段によって決定された位置に配置された各粒子において、隣り合う粒子がバネによって互いに接続されていると想定したときに、前記バネの力による各粒子の位置の変化を推定し、前記解析対象における粒子の分布の粗密を減ずるように、今回のタイムステップにおける各粒子を再配置する再配置手段と、を備える。

この態様において、前記再配置手段は、隣り合う粒子間の距離をそろえるように、今回のタイムステップにおける各粒子を再配置すべく構成されていてもよい。

また、上記態様において、前記再配置手段は、前記バネの力により変化した各粒子の位置を、その変化量が所定の許容値以下に収束するまで繰り返し演算するように構成されていてもよい。

また、上記態様において、前記流動解析装置は、前記再配置手段により再配置された粒子の物理的状態に関する粒子情報を、再配置される前の各粒子の粒子情報に基づいて取得する粒子情報取得手段をさらに備えていてもよい。

また、上記態様において、前記解析対象は、微圧縮性流体であってもよい。

また、本発明の一の態様の流動解析方法は、コンピュータが粒子から構成された解析対象における粒子の位置の時間変化を推定する流動解析方法であって、前記コンピュータが、前回のタイムステップにおける各粒子の位置及び各粒子の物理的状態に関する粒子情報に基づき、各粒子の速度を決定し、決定された速度に基づいて今回のタイムステップにおける各粒子の位置を決定し、前記コンピュータが、決定された位置に配置された各粒子において、隣り合う粒子がバネによって互いに接続されていると想定したときに、前記バネの力による各粒子の位置の変化を推定し、前記解析対象における粒子の分布の粗密を減ずるように、今回のタイムステップにおける各粒子を再配置する。

また、本発明の一の態様のコンピュータプログラムは、コンピュータに、粒子から構成された解析対象における粒子の位置の時間変化を推定させるためのコンピュータプログラムであって、前記コンピュータを、前回のタイムステップにおける各粒子の位置及び各粒子の物理的状態に関する粒子情報に基づき、各粒子の速度を決定し、決定された速度に基づいて今回のタイムステップにおける各粒子の位置を決定する位置決定手段、及び、前記位置決定手段によって決定された位置に配置された各粒子において、隣り合う粒子がバネによって互いに接続されていると想定したときに、前記バネの力による各粒子の位置の変化を推定し、前記解析対象における粒子の分布の粗密を減ずるように、今回のタイムステップにおける各粒子を再配置する再配置手段、として機能させる。

本発明に係る流動解析装置、流動解析方法、及びコンピュータプログラムによれば、解析結果における粒子の分布の粗密の発生を抑制し、高精度に物体の流動解析を行うことが可能となる。

実施の形態に係る流動解析装置の構成を示すブロック図。実施の形態に係る流動解析装置による流動解析処理の手順を示すフローチャート。粒子の再配置処理を説明するための模式図。粒子の再配置処理を説明するための模式図。粒子の再配置処理を説明するための模式図。粒子の再配置処理を説明するための模式図。粒子情報の取得処理を説明するための模式図。粒子情報の取得処理を説明するための模式図。評価試験における二重円筒問題の解析モデルを説明するための模式図。従来手法による材料Ａの内筒１回転後の粒子位置の解析結果を示す図。本手法によって解析した全粒子についての１回転後の周方向速度と粒子の半径方向座標との関係を示すグラフ。従来手法によって解析した全粒子についての１回転後の周方向速度と粒子の半径方向座標との関係を示すグラフ。評価試験における回転実験装置と解析モデルとを説明するための模式図。実験装置の観察結果及び本手法による解析結果を示す図。実験装置の観察結果及び本手法による解析結果を示す図。実験装置の観察結果及び本手法による解析結果を示す図。実験装置の観察結果及び本手法による解析結果を示す図。実験装置の観察結果及び本手法による解析結果を示す図。実験装置の観察結果、本手法による解析結果、及び従来手法による解析結果を比較する図。

以下、本発明の好ましい実施の形態を、図面を参照しながら説明する。

本実施の形態に係る流動解析装置は、解析対象である物体の運動を粒子法によりシミュレートするものである。つまり、流動解析装置は、解析対象を複数の粒子から構成されたものとして、粒子の運動方程式に基づいて各粒子の位置をタイムステップ毎に算出する。解析対象は、塑性変形又は流動する連続体、又は連続体とみなすことができる物体であり、具体的には、樹脂、ゴム、プラスチック、鋳造時の金属等の高粘性流体、水のような低粘性流体、鍛造時、塑性変形時における金属等の固体、セメント等の粉体を含む。なお、以下の説明では、微圧縮性流体である高粘性流体を解析対象としている。

［流動解析装置の構成］
図１は、本実施の形態に係る流動解析装置の構成を示すブロック図である。流動解析装置１は、コンピュータ１０によって実現される。図１に示すように、コンピュータ１０は、本体１１と、入力部１２と、表示部１３とを備えている。本体１１は、ＣＰＵ１１１、ＲＯＭ１１２、ＲＡＭ１１３、ハードディスク１１５、読出装置１１４、入出力インタフェース１１６、及び画像出力インタフェース１１７を備えており、ＣＰＵ１１１、ＲＯＭ１１２、ＲＡＭ１１３、ハードディスク１１５、読出装置１１４、入出力インタフェース１１６、及び画像出力インタフェース１１７は、バスによって接続されている。

ＣＰＵ１１１は、ＲＡＭ１１３にロードされたコンピュータプログラムを実行することが可能である。そして、流動解析用のコンピュータプログラムである流動解析プログラム１１０を当該ＣＰＵ１１１が実行することにより、コンピュータ１０が流動解析装置１として機能する。

ＲＯＭ１１２は、マスクＲＯＭ、ＰＲＯＭ、ＥＰＲＯＭ、又はＥＥＰＲＯＭ等によって構成されており、ＣＰＵ１１１に実行されるコンピュータプログラム及びこれに用いるデータ等が記録されている。

ＲＡＭ１１３は、ＳＲＡＭまたはＤＲＡＭ等によって構成されている。ＲＡＭ１１３は、ハードディスク１１５に記録されている流動解析プログラム１１０の読み出しに用いられる。また、ＣＰＵ１１１がコンピュータプログラムを実行するときに、ＣＰＵ１１１の作業領域として利用される。

ハードディスク１１５は、オペレーティングシステム及びアプリケーションプログラム等、ＣＰＵ１１１に実行させるための種々のコンピュータプログラム及び当該コンピュータプログラムの実行に用いられるデータがインストールされている。流動解析プログラム１１０も、このハードディスク１１５にインストールされている。

ハードディスク１１５には、例えば米マイクロソフト社が製造販売するWindows（登録商標）等のオペレーティングシステムがインストールされている。以下の説明においては、本実施の形態に係る流動解析プログラム１１０は当該オペレーティングシステム上で動作するものとしている。

読出装置１１４は、フレキシブルディスクドライブ、ＣＤ−ＲＯＭドライブ、またはＤＶＤ−ＲＯＭドライブ等によって構成されており、可搬型記録媒体１２０に記録されたコンピュータプログラムまたはデータを読み出すことができる。また、可搬型記録媒体１２０には、コンピュータを流動解析装置として機能させるための流動解析プログラム１１０が格納されており、コンピュータ１０が当該可搬型記録媒体１２０から流動解析プログラム１２０を読み出し、当該流動解析プログラム１２０をハードディスク１１５にインストールすることが可能である。

入出力インタフェース１１６は、例えばUSB，IEEE1394，又はRS-232C等のシリアルインタフェース、SCSI，IDE，又は IEEE1284等のパラレルインタフェース、及びＤ／Ａ変換器、Ａ／Ｄ変換器等からなるアナログインタフェース等から構成されている。入出力インタフェース１１６には、キーボード及びマウスからなる入力部１２が接続されており、ユーザが当該入力部１２を使用することにより、コンピュータ１０にデータを入力することが可能である。

画像出力インタフェース１１７は、ＬＣＤまたはＣＲＴ等で構成された表示部１３に接続されており、ＣＰＵ１１１から与えられた画像データに応じた映像信号を表示部１３に出力するようになっている。表示部１３は、入力された映像信号にしたがって、画像（画面）を表示する。

［流動解析装置の動作］
以下、本実施の形態に係る流動解析装置１の動作について説明する。

流動解析装置１は、以下に説明するような流動解析処理を実行して、解析対象の運動をシミュレートする。本実施の形態に係る流動解析装置１では、粒子法の一つであるＥＦＧＭ（Element-free Galerkin Method）を利用して流動シミュレーションを行う。

図２は、本実施の形態に係る流動解析装置１による流動解析処理の手順を示すフローチャートである。

流動解析処理において、まずＣＰＵ１１１は、計算条件を設定する（ステップＳ１）。この処理において、各種パラメータが設定される。

次にＣＰＵ１１１は、初期粒子配置を設定する（ステップＳ２）。初期粒子配置では、隣り合う粒子間の距離が等しいものとされる。なお、ここでは初期粒子配置において、隣り合う粒子間の距離を等しくないものとしてもよい。後述するステップＳ７の処理により、隣り合う粒子間の距離がそろえられるためである。

次にＣＰＵ１１１は、注目粒子に近接する粒子を探索する（ステップＳ３）。この処理では、注目粒子の運動に影響を与える粒子が存在する領域である影響領域の半径が予め定められており、注目粒子を中心とした影響領域に含まれる粒子が探索される。１つの注目粒子について近接粒子の探索が終了すると、ＣＰＵ１１１は次の注目粒子を設定し、近接粒子を探索する。このようにして、全ての粒子に関して近接粒子の探索が行われる。

次にＣＰＵ１１１は、各粒子の運動方程式を作成する（ステップＳ４）。以下、粒子の運動方程式について説明する。

＜場の近似＞
本実施の形態では、場の近似のためにmoving least square（ＭＬＳ）補間を用いた。ＭＬＳでは、連続体の任意点ｘにおける量Φ（ｘ）は次式のように表される。
ここで、
であり、Φ_Ｉは粒子点ＩのΦ値、Ｎは粒子総数を示す。ｐ（ｘ）はｍ次多項式の基底ベクトルであり、例えば２次元でｍ＝１の場合、ｐ^Ｔ（ｘ）＝［1，ｘ，ｙ］となる。

また、
である。ｗ_Ｉは重み関数で、ここでは以下の式を用いた。
ここにｒ＝｜ｘ−ｘ_Ｉ｜、αは定数、ｒ_０は影響領域の半径である。例えば、α＝７、ｒ_０を初期粒子間隔の２．６倍とすることができる。

＜支配方程式の離散化＞
（１）支配方程式と汎関数
本実施の形態では、Ｒｅ数１．０以下程度の高粘性体を対象とする。そこで、慣性力を無視した以下の支配方程式を用いる。

式（７）乃至（９）に対して、以下の汎関数を定義する。
式（１０）の右辺第３項は速度境界条件を満たすために導入したペナルティ項であり、κは非常に大きな数（本実施の形態では１０^１０とした）である。

（２）離散化
任意点ｘの速度を式（１）を用いて次式で表す。

汎関数の停留条件と式（１１）から以下の方程式が得られる。
式（１３）及び（１４）において、添え字Ｉ，Ｊは、粒子Ｉ，Ｊに関する成分であることを表している。

２次元問題に関しては、各マトリックス成分は下式のようになる。
ここにμは粘性係数、λは体積圧縮係数である。

非圧縮性流体の場合λ＝∞であるが、本実施の形態では微圧縮性流体を仮定してλ＝１００μとした。このような微圧縮性問題でしばしば問題になるロッキング現象は、後述の評価試験では発生しなかった。

また、式（６）の重み関数を使用すると粒子点の値が補間値と同一にならない（φ（ｘ_Ｉ）≠φ_Ｉ）が、式（１３）、（１４）のペナルティ項によって、式（１１）で計算される境界上の速度を既定速度に一致させることができる。

（３）積分法
一般にＥＦＧＭでは数値積分用のセルが必要であり、これが計算負荷を上げる要因になっている。そこで、積分点を粒子点と一致させ(nodal integration)、積分セルを不要にした。この場合、式（１３）、（１４）は以下のようになる。
ここに、Ｖ_ａは粒子ａの体積、Ｓ_ａ ^ｕ及びＳ_ａ ^ｔのそれぞれは粒子ａのΓ_ｕ及びΓ_ｔ上の面積を表す。Ｖ_ａは、初期粒子が等間隔で配置されていることを前提に計算した。なお、ここではnodal integrationにともなって発生する圧力振動抑制のための処理は行わないが、当該処理を行うこととしてもよい（Beissel, S. et al.: Nodal integration of the element-free Galerkin method, Comput. Methods Appl. Mech. Engrg. , Vol. 139, pp49-74, 1996.を参照）。

図２の説明に戻る。ステップＳ４の処理を行った後、ＣＰＵ１１１は、粒子速度を算出する（ステップＳ５）。本実施の形態では、粒子速度を完全陰解法により求め、粒子座標の更新を２次のRunge-Kutta法を用いた。ステップＳ５の処理について、詳細に説明する。

まず、粒子速度ベクトルを下式のように規定し、式（１１）より粒子点の補間速度Ｕ_ｔを求める。
ここに、ｘは粒子座標ベクトルを、添字ｔは時刻ｔにおける値を示す。

次に、下式よりΔｔ／２後の粒子座標ベクトルを算出する。
ここに、Δｔは時間刻み幅を示している。つまり、ｔを前回のタイムステップとすると、今回のタイムステップはｔ＋Δｔとなる。

時刻ｔ＋Δｔ／２における粒子速度ベクトルは、下式より与えられる。

式（１１）により、粒子点の補間速度Ｕ_{ｔ＋Δｔ／２}を求める。

以上のようにして粒子速度を算出した後、ＣＰＵ１１１は、今回のタイムステップにおける粒子位置を算出し、粒子位置を更新する（ステップＳ６）。時刻ｔ＋Δｔにおける粒子位置は、次式によって与えられる。

μがひずみ速度依存性の場合には、ステップＳ６の処理を終了後ひずみ速度を計算してμを修正し、タイムステップを進めずに再度ステップＳ５に戻って、μの変化が許容値以下になるまで反復計算すればよい。あるいは、反復計算を実施する代りにΔｔを十分小さく取ってΔｔ前の時刻のひずみ速度に対するμを近似値として用いる方法もある。

以上のステップＳ１乃至Ｓ６の処理は、ＥＦＧＭにより実現した処理である。また、ステップＳ５及びＳ６の処理は、前回のタイムステップｔにおける各粒子の位置Ｘ_ｔ及び各粒子の物理的状態に関する粒子情報（速度、圧力、温度、応力、質量、体積、密度等）から、各粒子の速度を決定し、決定された速度を用いて今回のタイムステップｔ＋Δｔにおける各粒子の位置を決定する位置決定処理Ｓ１００である（図２参照）。

上記の位置決定処理Ｓ１００で得られた各粒子の位置では、解析対象中の粒子の分布に粗密が生じる、つまり、隣り合う粒子の間隔がそろっていない場合がある。そこで、ＣＰＵ１１１は、このような粒子の分布の粗密を減ずるよう、粒子を再配置する（ステップＳ７）。この処理では、再配置の前後において解析対象の体積が変化しないように粒子が再配置される。

ステップＳ７の処理について説明する。図３Ａ乃至図３Ｄは、粒子の再配置処理を説明するための模式図である。図３Ａは、粒子位置の更新前の粒子の位置の一例を示しており、図３Ｂは、図３Ａの状態から更新した粒子位置を示している。また、図３Ｃは、図３Ｂの状態から粒子を再配置する処理の概念を示しており、図３Ｄは、図３Ｂの状態から再配置した粒子の位置を示している。

図３Ａに示すように、粒子位置の更新前においては、隣り合う粒子間の距離がそろっている。つまり、各粒子を中心とした円形の衝突感知領域が設定されており、隣り合う粒子において衝突感知領域が接した状態とされ、２以上の衝突感知領域が重なり合ったり、衝突感知領域が互いに離隔したりしないようになっている。粒子表面から衝突感知領域の外縁までの距離が、衝突感知距離とされている。

図３Ａの状態から、算出した粒子速度で各粒子を移動させることで、粒子位置が更新され、図３Ｂに示す状態になる。図３Ｂでは、各粒子が近づくように移動された結果、複数の粒子の衝突感知領域が重なり合う、つまり、隣り合う粒子の間の距離が衝突感知距離の２倍よりも短くなっている。また、図３Ｂでは、隣り合う粒子の間の距離が衝突感知距離の２倍よりも短くなった場合を示しているが、隣り合う粒子の間の距離が衝突感知距離の２倍よりも大きくなる場合もある。

粒子位置が更新された後、隣り合う粒子をバネで接続することを想定する。このバネは、自然長が衝突感知距離の２倍とされている。つまり、バネが自然長のとき、隣り合う粒子の間の距離は衝突感知距離の２倍となる。図３Ｃに示すように、隣り合う粒子間の距離が衝突感知距離の２倍よりも短いとき、これらの粒子を繋ぐバネは自然長よりも縮短しているので、これらの粒子に反発力を作用させ、粒子を互いに離反させる。また、隣り合う粒子間の距離が衝突感知距離の２倍よりも大きいときには、これらの粒子は衝突していないものとして、粒子間力は作用しないものとする。ただし、粒子間力の吸引力を無視できない流動体の場合には、これらの粒子に吸引力を作用させ、粒子を互いに近接させる。また、上記バネの運動モデルには減衰率を設定しておく。

バネの運動モデルに減衰率を設定したことにより、バネによる粒子間距離の調整を繰り返し実行すると、隣り合う粒子間の距離は衝突感知距離の２倍に近づいていく。隣り合う粒子間の距離と衝突感知距離の２倍との差が許容値以下となったとき、ＣＰＵ１１１は粒子の再配置の処理を終了する。これにより、図３Ｄに示すように、隣り合う粒子の衝突感知領域が重なったり互いに離隔したりせず、接する状態となる。つまり、隣り合う粒子間の距離がそろった状態となり、粒子の分布の粗密が解消される。

なお、ここでいう「距離をそろえる」とは、隣り合う粒子間の距離が一定となるだけでなく、数値計算上の誤差又は上記の許容値等によるバラツキがあっても、実質的に粒子間の距離が等しくなることも含む。

粒子の再配置処理をさらに具体的に説明する。粒子の再配置処理は、以下の手順で行われる。

まず、ＣＰＵ１１１は、下式により粒子速度及び粒子位置を初期設定する。

次に、ＣＰＵ１１１は、粒子に作用する力Ｑを次式で計算する。

次に、ＣＰＵ１１１は、粒子速度及び粒子位置を次式で更新する。

上記の粒子に作用する力Ｑの算出ステップ、並びに粒子速度及び粒子位置の更新ステップを、ｘの変動が許容値以下になるまで繰り返し実行する。

上記の粒子の再配置処理によれば、粒子間距離がｓ_０以下になった場合に反発力が働き、一定の繰り返し後粒子間距離ｓ_０の均等配置となる。なお、ｋ、ｃ、Δｔには物理的意味がないので、繰り返し数を少なくするよう自由に設定することができる。

図２の説明に戻る。上記の粒子の再配置処理を終了した後、ＣＰＵ１１１は、再配置後の粒子情報を取得する（ステップＳ８）。

ステップＳ８の処理について説明する。図４Ａ及び図４Ｂは、粒子情報の取得処理を説明するための模式図である。図４Ａは、再配置された粒子を示しており、図４Ｂは、粒子情報の取得処理の概念を示している。

図４Ａにおいて、実線の円は再配置後の粒子を示しており、破線の円は再配置前の粒子を示している。いま、図中斜線を付した粒子を注目粒子Ｐ０とする。注目粒子Ｐ０を中心として円形の影響領域Ａ０が設定される。この注目粒子Ｐ０の粒子情報（速度、圧力、温度、応力、質量、体積、密度等）は、再配置前に前記影響領域Ａ０内に存在していた各粒子に影響を受けると考える。

図４Ｂでは、再配置前の粒子ＰＰの位置及び再配置後の注目粒子Ｐ０の位置が示されている。影響領域Ａ０内に存在する粒子ＰＰ（図中網掛けにて示す）の粒子情報に基づいて、次式により注目粒子Ｐ０の粒子情報が算出される。
ここで、Ｗは重み関数であり、Φは速度、圧力等の粒子情報を示している。

再配置後の粒子情報の取得処理が終了すると、ＣＰＵ１１１は、粒子のひずみ速度及び応力を計算する（ステップＳ９）。この処理において、粒子のひずみ速度は次式より求められる。
ここで、ε_ｉ，ｊはひずみ速度を、ｕ_ｉ，ｊ、ｕ_ｊ，ｉは速度勾配を示している。また、粒子の応力は、得られたひずみ速度に粘性係数を乗ずることにより求められる。

なお、上述したステップＳ７乃至Ｓ９の処理は、毎タイムステップ実行される必要はない。位置決定処理Ｓ１００において決定された各粒子の位置が、前回のタイムステップにおける位置から所定値未満しか変動しておらず、粒子の移動量が少ないと評価されるような場合には、ステップＳ７乃至Ｓ９の処理を実行しなくてもよい。これにより、解析の精度を損なうことなく、計算時間を抑制することができる。ただし、より高精度な解析が求められるような場合には、ステップＳ７乃至Ｓ９の処理を毎タイムステップ実行する構成とすることも可能である。

次にＣＰＵ１１１は、流動解析処理を終了するか否かを判定し（ステップＳ１０）、流動解析処理を終了しない場合には（ステップＳ１０においてＮＯ）、タイムステップを１つ進めて（ステップＳ１１）、ステップＳ３へと処理を移す。以降、ＣＰＵ１１１は、流動解析処理を終了すると判定するまで、ステップＳ３乃至Ｓ１１の処理を繰り返す。これにより、粒子の位置及び粒子情報の時間的な変化がシミュレートされ、時系列の粒子位置情報及び粒子情報が取得される。

ステップＳ１０において、流動解析処理を終了する場合には（ステップＳ１０においてＹＥＳ）、ＣＰＵ１１１は、解析結果を表示部１３に出力し（ステップＳ１２）、流動解析処理を終了する。表示される解析結果としては、あるタイムスステップにおける粒子位置を座標空間に描画した画像、又は粒子情報の数値データ等がある。また、粒子位置を座標空間に描画した画像を時系列で並べたり、時系列で表示遷移したりして、粒子の位置の時間変化をわかりやすくユーザに提示することも可能である。

以上の如く構成したことにより、本実施の形態に係る流動解析装置１は、粒子の分布の粗密が低減され、より実際の物体の運動に近い解析結果を得ることができる。

また、本実施の形態に係る流動解析装置１は、合成樹脂、ゴム等の混練時の挙動、プラスチック等の射出成形時の挙動、金属の鋳造時又は鍛造時の挙動、金属の塑性変形時の挙動、水のような低粘性流体の挙動等、様々な物体の挙動の解析に利用することが可能である。

（評価試験）
上記の実施の形態において説明した流動解析装置を実際に作成し、その性能を評価した。

（１）試験１
本発明者は、二重円筒内が粘性流体で満たされ、内筒が各速度ω_０で回転している場合の定常解を計算する試験を実施した。

本試験において、粘性係数は、下式で示されるべき乗則に従うものとした。

慣性項を省略した円筒座標系のＮＳ方程式を、流速の周方向成分のみを考慮して境界固着条件で解くと以下の解が得られる。
ここで、ｕは周方向速度を、τはせん断応力を、ｒは半径方向の座標を、Ｒは外筒半径を、ξは内筒半径と外筒半径との比を、Ｔは内筒（外筒）に作用するトルクを示している。

図５は、本試験における二重円筒問題の解析モデルを説明するための模式図である。図５には、初期配置された粒子が示されている。この初期配置では、隣り合う粒子の間隔がそろえられている。なお、図５に解析条件も示している。

本試験では、以下の材料を考えた。
材料Ａ：μ_０＝１０００Ｐａ・ｓ，α＝１．０
材料Ｂ：μ_０＝１１３２４．７６Ｐａ・ｓ，α＝０．３

本問題の場合、粒子の周方向速度は周方向位置又は時刻に依らずに式（２０）に一致せねばならないが、実際には粒子の移動とともに誤差が蓄積する。そこで、誤差を評価するために内筒の１回転後の結果を本手法と従来手法（従来のＥＦＧＭによる解析手法）とについて求め、両者を比較した。

図６は、従来手法による材料Ａの内筒１回転後の粒子位置の解析結果を示す図である。図６に示すように、従来手法による解析結果では、粒子配置が初期と比べて大きく乱れており、クラスタリングの発生も見られる。一方本手法においては、内筒１回転後も図５と同様粒子分布は概ね一様となった。

図７Ａ及び図７Ｂに、解析した全粒子についての１回転後の周方向速度と粒子の半径方向座標との関係を示す。図７Ａは、本手法による結果であり、図７Ｂは、従来手法による結果である。図７Ａ及び図７Ｂより、本手法による結果は、厳密な数値計算により得られた解析結果（以下、「解析解」という。）とよく一致しているが、従来手法による結果は大きなばらつきを生じていることがわかる。特に材料Ｂ（非Newton流体）ではその傾向が強い。

次表は、１回転後のトルクの比較結果を示している。なお、ここでは解析解に対する誤差によって本手法と従来手法との比較を行っている。

トルクは、次式の境界粒子力から求めた。
ここに、Ｒ_Ｉは境界上の粒子Ｉに作用する力を、Ｕ_Ｉは式（１１）で計算される補間速度を示している。

上表から、本手法においては解析解に対するトルクの誤差が１%以下になることがわかる。また下表は、本手法において初期粒子間隔のトルクに及ぼす影響を調べた結果を示している。下表から、本手法では、粒子間隔４ｍｍ（流路幅方向に粒子が４体程度）と、かなり疎な配置でも精度を維持できることがわかる。

（２）試験２
本発明者は、図８に示すような円筒バレルと円柱ロータから構成された回転実験装置内に、シリコンオイルを充填した後、円柱ロータを回転させてその挙動をビデオ観察した。また、同様の解析モデルについて本手法及び従来手法によって解析を行い、ビデオ観察した実際の流体運動と、解析結果とを比較した。なお、図８には、解析条件も示している。

シリコンオイルの粘性係数は、室温でせん断ひずみ速度が５０ｓ^−１以下では概ね１００Ｐａ・ｓである。なお、実験では流体内に多数の気泡が見られるが、気泡の少ない実験との比較から、気泡の流体自由表面形状に及ぼす影響はほとんど無いことを確認している。

図９Ａ乃至図９Ｅは、実験装置の観察結果及び本手法による解析結果を示す図である。図９Ａには１／４回転後の結果を、図９Ｂには１／２回転後の結果を、図９Ｃには３／４回転後の結果を、図９Ｄには１回転後の結果を、図９Ｅには３回転後の結果を示している。

なお、図９Ａ乃至図９Ｅにおける本手法の解析結果において、実験装置のビデオ観察結果から得られた自由表面形状を黒色の実線で示している。同図からわかるように、ロータに巻き取られた流体が１回転後に下部流体と合流する近傍において、本手法の解析結果の方が観察結果より流体量がやや多くなるなど若干の相違が見られるが、解析結果と観察結果とは概ね一致している。なお、本解析に要した計算時間は1 回転当り528秒であった（i5-3210M，2.5GHz，1コア）。

図１０は、実験装置の観察結果、本手法による解析結果、及び従来手法による解析結果を比較する図である。図１０では、ロータを１回転させたときの結果を示している。図１０では、本手法による解析結果では粒子の分布が均一であり、自由表面形状が観察結果とよく合致していることがわかる。一方、従来手法では、粒子の分布に粗密が生じた結果、ロータの周囲において流体の欠落が生じており、自由表面形状が観察結果と合致していない。

（その他の実施の形態）
なお、上述した実施の形態においては、粒子の位置が更新された後、隣り合う粒子間の間隔がそろうように粒子の再配置を行う構成について述べたが、これに限定されるものではない。隣り合う粒子間の間隔がそろうとはいえなくても、粒子の分布の粗密が減少するように粒子を再配置する構成とすることも可能である。例えば、粒子の再配置処理において、隣り合う粒子間の距離と衝突感知距離の２倍との差の許容値を比較的大きくして、隣り合う粒子間の距離の制限を緩和することができる。このようにすることによって、粒子間距離の調整の繰り返し回数が抑制され、解析結果で十分な精度を確保しつつ、流動解析処理の計算時間を抑えることができる。また、これにより、粉体のような圧縮性の高い材料についても解析することが可能となる。

また、上述した実施の形態においては、粒子の再配置処理においてバネの運動モデルを利用する構成について述べたが、これに限定されるものではない。隣り合う粒子の間隔が一定とされた粒子配置のテンプレートを予め設けておき、粒子位置が更新された後、更新された粒子が存在する領域に前記テンプレートを適用して粒子を置き換えることで、粒子を再配置することも可能である。なお、ここでいう「粒子の再配置」とは、処理の前後における個々の粒子の対応関係を維持したまま、各粒子の位置関係を変更することだけでなく、処理前の粒子と対応させることなく新たな粒子に置き換えて、各粒子の位置関係を変更することも含む。

また、上述した実施の形態においては、粒子の再配置処理の後に粒子情報取得処理を実行し、再配置前に注目粒子の影響領域内に存在していた各粒子の粒子情報に基づいて、注目粒子の粒子情報を算出する構成について述べたが、これに限定されるものではない。粒子情報取得処理を省くことも可能である。この場合、粒子の再配置処理における粒子間距離の調整の繰り返し時間間隔Δｔを十分に小さく設定し、粒子間距離の調整を細かくすることで、１回の調整における粒子の座標変動量を小さくして粒子情報の精度を確保すればよい。

また、上述した実施の形態においては、流動解析処理において、ＥＦＧＭを利用して各粒子の速度を決定し、決定された速度を用いて各粒子の位置を決定する構成について述べたが、これに限定されるものではない。ＥＦＧＭとは異なる粒子法であるＭＰＳ（Moving Particle Semi-implicit）又はＳＰＨ（Smoothed Particle Hydrodynamics）等を利用して、各粒子の速度を決定し、決定された速度を用いて各粒子の位置を決定する構成とすることも可能である。

本発明の流動解析装置、流動解析方法、及びコンピュータプログラムは、粒子法を用いて流動性を有する物体の流動解析を行う流動解析装置、流動解析方法、及びコンピュータプログラムとして有用である。

１流動解析装置
１０コンピュータ
１２入力部
１３表示部
１１０流動解析プログラム
１１１ＣＰＵ
１１５ハードディスク
１１６入出力インタフェース
１１７画像出力インタフェース

Claims

粒子から構成された解析対象における粒子の位置の時間変化を推定する流動解析装置であって、
前回のタイムステップにおける各粒子の位置及び各粒子の物理的状態に関する粒子情報に基づき、各粒子の速度を決定し、決定された速度に基づいて今回のタイムステップにおける各粒子の位置を決定する位置決定手段と、
前記位置決定手段によって決定された位置に配置された各粒子において、隣り合う粒子がバネによって互いに接続されていると想定したときに、前記バネの力による各粒子の位置の変化を推定し、前記解析対象における粒子の分布の粗密を減ずるように、今回のタイムステップにおける各粒子を再配置する再配置手段と、
を備える、
流動解析装置。
前記再配置手段は、隣り合う粒子間の距離をそろえるように、今回のタイムステップにおける各粒子を再配置すべく構成されている、
請求項１に記載の流動解析装置。
前記再配置手段は、前記バネの力により変化した各粒子の位置を、その変化量が所定の許容値以下に収束するまで繰り返し演算するように構成されている、
請求項１又は２に記載の流動解析装置。
前記再配置手段により再配置された粒子の物理的状態に関する粒子情報を、再配置される前の各粒子の粒子情報に基づいて取得する粒子情報取得手段をさらに備える、
請求項１乃至３の何れかに記載の流動解析装置。
前記解析対象は、微圧縮性流体である、
請求項１乃至４の何れかに記載の流動解析装置。
コンピュータが粒子から構成された解析対象における粒子の位置の時間変化を推定する流動解析方法であって、
前記コンピュータが、前回のタイムステップにおける各粒子の位置及び各粒子の物理的状態に関する粒子情報に基づき、各粒子の速度を決定し、決定された速度に基づいて今回のタイムステップにおける各粒子の位置を決定し、
前記コンピュータが、決定された位置に配置された各粒子において、隣り合う粒子がバネによって互いに接続されていると想定したときに、前記バネの力による各粒子の位置の変化を推定し、前記解析対象における粒子の分布の粗密を減ずるように、今回のタイムステップにおける各粒子を再配置する、
流動解析方法。
コンピュータに、粒子から構成された解析対象における粒子の位置の時間変化を推定させるためのコンピュータプログラムであって、
前記コンピュータを、
前回のタイムステップにおける各粒子の位置及び各粒子の物理的状態に関する粒子情報に基づき、各粒子の速度を決定し、決定された速度に基づいて今回のタイムステップにおける各粒子の位置を決定する位置決定手段、及び、
前記位置決定手段によって決定された位置に配置された各粒子において、隣り合う粒子がバネによって互いに接続されていると想定したときに、前記バネの力による各粒子の位置の変化を推定し、前記解析対象における粒子の分布の粗密を減ずるように、今回のタイムステップにおける各粒子を再配置する再配置手段、
として機能させる、コンピュータプログラム。