(57)【要約】【課題】固体状態の不連続性と溶融状態の連続性の双方を考慮した物質の運動量を予測することができる演算プログラム、演算方法を提供する。【解決手段】粒状の物質の物性を含むパラメータと該物質に与えられる外部エネルギーとに基づいて時間の経過に伴う前記物質の温度変化を算出する第１の算出手段と、前記パラメータと不連続体の運動量を演算するための第１の粒子法とに基づいて前記物質の運動量を算出し、前記算出された物質の温度が予め設定された融点を超えた場合、前記第１の粒子法を、連続体の運動量を演算するための第２の粒子法に切り替えて前記物質の溶融状態の運動量を算出する第２の算出手段として機能させる。【選択図】図２
2014104696.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを、
粒状の物質の物性を含むパラメータと該物質に与えられる外部エネルギーとに基づいて時間の経過に伴う前記物質の温度変化を算出する第１の算出手段と、
前記パラメータと不連続体の運動量を演算するための第１の粒子法とに基づいて前記物質の運動量を算出し、前記算出された物質の温度が予め設定された融点を超えた場合、前記第１の粒子法を、連続体の運動量を演算するための第２の粒子法に切り替えて前記物質の溶融状態の運動量を算出する第２の算出手段
として機能させるための演算プログラム。
【請求項２】
前記第２の算出手段は、
前記算出された物質の温度が、前記融点より所定温度低い予め設定された遷移開始点から前記融点までの温度範囲内に含まれる場合、前記第１の粒子法を、前記第１及び第２の粒子法に基づく演算手法に切り替えて前記物質の固液混在状態の運動量を算出する
請求項１記載の演算プログラム。
【請求項３】
前記所定温度は、前記物質の温度変化による比熱の変化に基づく値である
請求項２記載の演算プログラム。
【請求項４】
前記所定温度は、前記物質の融点が高いほど、大きい値が用いられる
請求項２または請求項３記載の演算プログラム。
【請求項５】
前記所定温度は、１℃以上、２０℃以下である
請求項２乃至請求項４の何れか１項に記載の演算プログラム。
【請求項６】
前記演算手法においては、前記固液混在状態の運動量を算出するための式（１１）が用いられる
請求項２乃至請求項５の何れか１項に記載の演算プログラム。
【数１１】
2014104696.tif 000017
【請求項７】
前記第１の粒子法、前記第２の粒子法、及び前記演算手法においては、それぞれ互いに異なる動粘性係数が用いられる、
請求項２乃至請求項６の何れか１項に記載の演算プログラム。
【請求項８】
前記演算手法においては、動粘性係数を算出するための式（１２）が用いられる
請求項２乃至請求項７の何れか１項に記載の演算プログラム。
【数１２】
2014104696.tif 000018
【請求項９】
前記第１の算出手段は、粘性係数を算出するための式（７）を少なくとも用いて前記物質の温度を算出する
請求項１乃至請求項８の何れか１項に記載の演算プログラム。
【数７】
2014104696.tif 000019
【請求項１０】
前記第１の粒子法はＤＥＭ法であり、前記第２の粒子法はＭＰＳ法である
請求項１乃至請求項９の何れか１項に記載の演算プログラム。
【請求項１１】
前記第１の粒子法はＤＥＭ法であり、前記第２の粒子法はＳＰＨ法である
請求項１乃至請求項９の何れか１項に記載の演算プログラム。
【請求項１２】
粒状の物質の物性を含むパラメータと該物質に与えられる外部エネルギーとに基づいて時間の経過に伴う前記物質の温度変化を算出し、
前記パラメータと不連続体の運動量を演算するための第１の粒子法とに基づいて前記物質の運動量を算出し、前記算出された物質の温度が予め設定された融点を超えた場合、前記第１の粒子法を、連続体の運動量を演算するための第２の粒子法に切り替えて前記物質の溶融状態の運動量を算出する
ことをコンピュータが実行する演算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物質の溶融可塑化工程を解析する演算プログラム及び演算方法に関するものである。
【背景技術】
【０００２】
数値流体力学分野において、演算を行う際に流体が存在する領域を格子状に要素分割し、要素内の計算点の物理量を連立方程式にて求める手段が一般的であり、有限要素法（Ｆｉｎｉｔｅ  Ｅｌｅｍｅｎｔ  Ｍｅｔｈｏｄ：ＦＥＭ）、有限体積法（Ｆｉｎｉｔｅ  Ｖｏｌｕｍｅ  Ｍｅｔｈｏｄ：ＦＶＭ）や有限差分法（Ｆｉｎｉｔｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｍｅｔｈｏｄ：ＦＤＭ）が代表的である。中でも樹脂流動解析分野において、連続体である溶融体について押出成形機のスクリュ混練解析や樹脂配管、ダイスの流動性評価解析に対しこれら格子法が多く採用され、現在ではレオロジー特性評価技術は成熟の段階に至っている。
【０００３】
しかしながら、粉末あるいはペレット状の樹脂原料（その状態により、適宜固体樹脂、溶融樹脂と称する）は、個別に運動を行う不連続体であることから、格子メッシュ法でその運動量を予測することが困難とされている。一部でＶＯＦ（Ｖｏｌｕｍｅ  ｏｆ  Ｆｌｕｉｄ）法を用いた粉体原料の解析例がみられるものの、その不連続性予測が困難なことや自由表面状態の予測精度について、その定量性の確保が難しい等の課題を有する。このため、不連続性を考慮した固体樹脂を供給し装置内で溶融可塑化を行い、溶融状態で吐出される押出成形機や射出成形機などのスクリュ挿入型のプラスチック溶融可塑化装置内の全領域を格子メッシュ法で一貫解析することは困難とされている。
【０００４】
一方で、押出成形機や射出成形機内の溶融可塑化工程を含む一貫解析が可能なものに、ＦＡＮ（Ｆｌｏｗ  Ａｎａｌｙｓｉｓ  Ｎｅｔｗｏｒｋ）法を解法に用いた技術が知られている（下記特許文献１参照）。この演算においては、入力したスクリュ形状から、そのスクリュが有する輸送量や粘性発熱特性を算出し、入力したスクリュ回転数や押出量、シリンダ設定温度などの操作条件から樹脂原料が受けるせん断発熱エネルギーや伝熱エネルギーを求め、そのエネルギーを樹脂の比熱や熱伝導率で換算することで樹脂温度の上昇率を演算する。この演算により結果的に固体樹脂が溶融可塑化される領域を求めることができ、押出機内の充満率や滞留時間なども合わせた物理量分布を一次元的に予測することが可能となる。
【０００５】
ただし、このＦＡＮ法による演算においても、固体樹脂から溶融樹脂までを一貫して連続体と取り扱うため、固体樹脂の不連続性はもとより、ペレット個々の物理量履歴などの詳細予測は不可能であり、結果としてＦＡＮ法は、押出機内のマクロ的な樹脂挙動を予測する手段として実用化されている。即ち、比較的ミクロ的視点にて樹脂の挙動解析を行う格子メッシュ法では固体樹脂の不連続性を適切に予測が行えないことから、溶融可塑化工程の演算が極めて困難である。また、溶融可塑化の演算が可能なＦＡＮ法解析ではマクロ的視点の解析手段であるため、固体樹脂個々のミクロ的挙動の追跡予測が行えず、結果としてミクロ的視点による樹脂の溶融可塑化挙動を予測する演算手法が確立されていない。
【０００６】
一方、流体のミクロ的演算手法として、流体要素を格子状ではなく個々の粒子にモデル化して運動状態を予測する粒子法がある。格子メッシュ法にＦＥＭやＦＶＭ、ＶＯＦ法など数種の演算手法があるのと同様、粒子法においてもいくつかの演算手法が存在する。この粒子法の中で、固体樹脂の粒状サイズ（直径が数十μｍから数ｍｍ）のモデル化に適した手法として、ＭＰＳ（Ｍｏｖｉｎｇ  Ｐａｒｔｉｃｌｅ  Ｓｅｍｉ−ｉｍｐｌｉｃｉｔ）法（下記特許文献２、３参照）、ＳＰＨ（Ｓｍｏｏｔｈｅｄ  Ｐａｒｔｉｃｌｅ  Ｈｙｄｒｏｄｙｎａｍｉｃｓ）法、ＤＥＭ（Ｄｉｓｔｉｎｃｔ  Ｅｌｅｍｅｎｔ  Ｍｅｔｈｏｄ）法がある。
【０００７】
ＭＰＳ法とＳＰＨ法は連続体の数値解析を前提として開発された演算手法であり、連続体の運動量の保存を表すナビエ・ストークス方程式を粒子個々に適用して計算することを特徴とする。そのため、その演算結果は格子メッシュ法と同様の精度で流体の運動量を示すことができる。これらＭＰＳ法とＳＰＨ法を比較すると、ＭＰＳ法が非圧縮性流体を対象として開発されたのに対し、ＳＰＨ法は圧縮性流体を対象とした開発が行われた部分が根本的に大きな違いであり、それゆえ粒子の物理量演算手法が異なる（下記非特許文献１参照）。
【０００８】
ＤＥＭ法は不連続体の運動を解く演算手法として開発されたものであり、粉体の運動状態の予測など、固体粒子の集合体を対象とした予測解析に実績を多く有する（下記非特許文献２、３参照）。
【０００９】
物質が固−液−気体状態の相変化を起こした場合、物質の運動量が極端に変化する。液体が気化する蒸発挙動の予測は、前記の格子メッシュ法で解析することが可能な汎用ソフトウェア（例えば、アドバンスソフト（株）「Ａｄｖａｎｃｅ／ＦｒｏｎｔＦｌｏｗ／ＭＰ」）にて解析することが可能である。また、一部のソフトウェア（例えば、（株）ソフトウェアクレイドル「ＳＴＲＥＡＭ」）では、固−液間の溶融可塑化あるいは凝固挙動を予測することも可能である。また、粒子法を用いた相変化の予測例は、ＭＰＳ法による液体から気体へ変化する挙動予測手法が知られている（下記非特許文献４参照）。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】特許第４４１４４０８号公報
【特許文献２】特開平０７−３３４４８４号公報
【特許文献３】特開２００９−１２９１９３号公報
【非特許文献】
【００１１】
【非特許文献１】越塚  誠一  著「数値流体力学」培風館  １９９７．４
【非特許文献２】酒井  幹夫  著「粉体の数値シミュレーション」丸善出版  ２０１２．８
【非特許文献３】Ｓ．Ｋｏｓｈｉｚｕｋａ，Ａ．Ｎｏｂｅ  ａｎｄ  Ｙ．Ｏｋａ  “Ｎｕｍｅｒｉｃａｌ  Ａｎａｌｙｓｉｓ  ｏｆ  Ｂｒｅａｋｉｎｇ  Ｗａｖｅｓ  Ｕｓｉｎｇ  ｔｈｅ  Ｍｏｖｉｎｇ  Ｐａｒｔｉｃｌｅ  Ｓｅｍｉ−ｉｍｐｌｉｃｉｔ  Ｍｅｔｈｏｄ”Ｉｎｔ．  Ｊ．Ｎｕｍｅｒ．Ｍｅｔｈ．Ｆｌｕｉｄｓ  ２６，７５１−７６９(１９９８)
【非特許文献４】越塚  誠一  著「計算力学レクチャーシリーズ５  粒子法」丸善出版  ２００５．２
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかしながら、蒸発挙動の予測や、固−液間の溶融可塑化あるいは凝固挙動の予測は、前述のとおり全ての物質を連続体として解いているため、粒状の固体樹脂の不連続性を予測できるものではない。前述した粒子法を用いた相変化の予測においても、物質の状態を全て連続体とした計算である。固形塊から液体への相変化をＭＰＳ法にて演算した研究事例もあるが、これも固体樹脂が連続体と近似できる形態から演算を開始しており、固体樹脂の不連続性を予測できるものではない。また、非特許文献２に記載のＤＥＭ法にて固体粒子をモデル化して物質の不連続挙動を演算し、同時にＭＰＳ法にて液体を連続体として演算を行う、ＤＥＭ−ＭＰＳ法による固液二相流の解析事例が記載されているが、この演算ではＤＥＭ粒子は固体状態を維持したままであり、溶融可塑化の演算が行えることを示したものではない。
【００１３】
以上のことから、粒状の固体樹脂の不連続性を考慮した運動量演算を行うと同時に外部エネルギーを与えることでその温度上昇を予測し、融点を超えた段階で液体状態へ相変化を行い、それ以上の温度領域では液体の連続性を考慮した運動量を演算することが困難であることが問題となっていた。
【００１４】
本発明は上述した問題点を解決するためになされたものであり、固体状態の不連続性と溶融状態の連続性の双方を考慮した物質の運動量を予測することができる演算プログラム及び演算方法を提供することを目的とする。
【課題を解決するための手段】
【００１５】
上述した課題を解決するため、本発明の一態様は、コンピュータを、粒状の物質の物性を含むパラメータと該物質に与えられる外部エネルギーとに基づいて時間の経過に伴う前記物質の温度変化を算出する第１の算出手段と、前記パラメータと不連続体の運動量を演算するための第１の粒子法とに基づいて前記物質の運動量を算出し、前記算出された物質の温度が予め設定された融点を超えた場合、前記第１の粒子法を、連続体の運動量を演算するための第２の粒子法に切り替えて前記物質の溶融状態の運動量を算出する第２の算出手段として機能させる。
【００１６】
また、本発明の一態様は、粒状の物質の物性を含むパラメータと該物質に与えられる外部エネルギーとに基づいて時間の経過に伴う前記物質の温度変化を算出し、前記パラメータと不連続体の運動量を演算するための第１の粒子法とに基づいて前記物質の運動量を算出し、前記算出された物質の温度が予め設定された融点を超えた場合、前記第１の粒子法を、連続体の運動量を演算するための第２の粒子法に切り替えて前記物質の溶融状態の運動量を算出する。
【発明の効果】
【００１７】
本発明によれば、固体状態の不連続性と溶融状態の連続性の双方を考慮した物質の運動量を予測することができる。
【図面の簡単な説明】
【００１８】
【図１】本実施の形態に係るシミュレーション装置のハードウェア構成を示すブロック図である。
【図２】本実施の形態に係るシミュレーション装置の機能構成を示す機能ブロック図である。
【図３】樹脂原料の比熱曲線と動粘性係数の関係を示す模式図である。
【図４】本実施の形態に係るシミュレーション装置の溶融可塑化工程解析処理を示すフローチャートである。
【図５】本実施の形態に係るシミュレーション装置の溶融可塑化工程解析処理を示すフローチャートである。
【図６】オフセット温度を変化させた場合の溶融可塑化工程解析処理結果を示す図である。
【図７】樹脂原料をシリンダに配置させた状態を可視化した図である。
【図８】解析開始後１００秒後における樹脂原料の状態を可視化した解析結果を示す図である。
【図９】解析開始後１５０秒後における樹脂原料の状態を可視化した解析結果を示す図である。
【図１０】解析開始後２００秒後における樹脂原料の状態を可視化した解析結果を示す図である。
【図１１】解析開始後３００秒後における樹脂原料の状態を可視化した解析結果を示す図である。
【図１２】シミュレーションプログラムが情報処理装置に適用される場合を示す模式図である。
【発明を実施するための形態】
【００１９】
本発明の実施の形態は、物質の溶融可塑化工程を演算するシミュレーション手段において、固体状態の樹脂原料を粒子状にモデル化し、ＤＥＭ法にて運動量の演算を行い、演算結果が予め設定された設定融点を超えた場合に、同一粒子の運動量演算をＭＰＳ法あるいはＳＰＨ法に切り替えることで固体状態の樹脂原料の不連続性と溶融状態の樹脂原料の連続性の双方を同時に演算可能とする手段を構築したものである。以下、本実施の形態について図面を参照しつつ、その詳細を説明する。
【００２０】
図１は、本実施の形態に係るシミュレーション装置のハードウェア構成を示すブロック図である。図１に示されるように、シミュレーション装置１は、ＣＰＵ(Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ)１１、記憶部１２、入力部１３、表示部１４、ＨＤＤ（Ｈａｒｄ  ｄｉｓｋ  ｄｒｉｖｅ）１５、を有する。
【００２１】
ＣＰＵ１１は記憶部１２上に展開されるＯＳ(Ｏｐｅｒａｔｉｎｇ  Ｓｙｓｔｅｍ)、ＢＩＯＳ(Ｂａｓｉｃ  Ｉｎｐｕｔ／Ｏｕｔｐｕｔ  Ｓｙｓｔｅｍ)、アプリケーション等の各種プログラムを実行し、シミュレーション装置１の制御を行う。記憶部１２はいわゆるＲＡＭ(Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ)などの揮発性のメモリであり、実行されるプログラムの作業領域として利用される。
【００２２】
入力部１３は、シミュレーション装置１を使用する使用者からの入力を受け付けるものであり、例えば、ディスプレイ上の特定の位置を指定するためのポインティングデバイスであるマウスや、文字または特定の機能等が割り当てられた複数のキーが配列されたキーボードである。表示部１４は、ＯＳおよびＯＳ上で動作するアプリケーションのＧＵＩ（Ｇｒａｐｈｉｃａｌ  Ｕｓｅｒ  Ｉｎｔｅｒｆａｃｅ）や解析結果を表示するディスプレイ等の出力装置である。ＨＤＤ１５は、固体樹脂または溶融樹脂である原料の溶融可塑化工程を予測解析する溶融可塑化工程解析処理を行うためのシュミレーションプログラム１０１が格納されている。
【００２３】
次に、シミュレーション装置１の機能構成を説明する。図２は、本実施の形態に係るシミュレーション装置１の機能構成を示す機能ブロック図である。図２に示されるように、シミュレーション装置１は、取得部１１１と、判定部１１２と、演算部１１３と、出力部１１４とを、機能として有する。これら機能は、ＣＰＵ１１や記憶部１２等のハードウェア資源とシミュレーションプログラム１０１とが協働することにより実現される。
【００２４】
取得部１１１は、解析処理における初期設定条件等が入力部に入力された場合、これを取得する。例えば、粒状の樹脂原料の物性を含むパラメータや該樹脂原料に与えられる外部エネルギー等の初期設定条件、シリンダ内の樹脂の配置情報等の情報や、ＭＰＳ法またはＳＰＨ法の何れを用いるのか等の情報、溶融可塑化工程解析処理の終了条件（規定時間により終了や演算回数により終了等）の情報を取得する。また、後述する設定融点Ｔｍ、オフセット温度Ｔｏｆｆｓｅｔ等の情報を取得する。パラメータは、使用者から解析処理の対象となる固体樹脂の形状、樹脂物性、初期粒子温度、総数等がある。外部エネルギーは、固体樹脂が配置されるシリンダの温度等がある。判定部１１２は、取得部１１１により取得された情報や演算部１１３の演算結果に基づいて、樹脂原料が設定融点Ｔｍを超えたか否か等の溶融可塑化工程解析処理における各種判定を行う。
【００２５】
演算部１１３は、取得部１１１により取得された取得情報に基づいて、樹脂原料を粒子状にモデル化し、溶融可塑化工程解析処理における各種演算処理を行うことにより、樹脂原料の溶融可塑化工程の予測解析を行う。例えば、時間の経過に伴う固体樹脂の温度変化や、取得情報と不連続体の運動量を演算するためのＤＥＭ法とに基づく固体樹脂の運動量の算出がある。また、演算部１１３は、算出された固体樹脂の温度が設定融点Ｔｍを超えたと判定部１１２により判定された場合、ＤＥＭ法を、連続体の運動量を演算するためのＭＰＳ法やＳＰＨ法に切り替えて固体樹脂の溶融状態の運動量を算出する。更に、演算部１１３は、算出された固体樹脂の温度が、設定融点Ｔｍより所定温度低い予め設定された遷移開始温度Ｔ０から、設定融点Ｔｍまでの温度範囲内に含まれる場合、ＤＥＭ法を、ＤＥＭ法及びＭＰＳ法に基づく演算手法に切り替えて固体樹脂の固液混在状態の運動量を算出する。出力部１１４は、演算部１１３の演算結果に基づいて樹脂原料の状態を可視化した解析結果を表示部１４に出力する。
【００２６】
次に、演算部１１３の演算処理について、その詳細を説明する。なお、ここでは重力などの外力が作用しない場合を例に取り、説明を行う。
【００２７】
まず、演算部１１３は、固体である粒状の樹脂原料を粒子（以後、原料粒子と称する）にモデル化し、ＤＥＭ法による運動量の計算を行う。ＤＥＭ法では、粒子どうしが接触・衝突した際にのみ粒子間相互作用が働く。この際の並進運動による運動方程式は（１）式に基づく。なお、この式において、ｍは質量、ｕは速度、ｔは時間、ＦＤＥＭは力、Ｆｇは重力である。ここでＦＤＥＭは接触力の大きさであり、粒子同士が衝突した際にその反発力やエネルギー減衰、摩擦力などから計算される（詳細は非特許文献２を参照）。
【００２８】
【数１】
2014104696.tif 000003
【００２９】
原料粒子の温度場を計算するエネルギー方程式は（２）式で表される。
【００３０】
【数２】
2014104696.tif 000004
【００３１】
演算部１１３は、原料粒子の温度上昇として、伝熱エネルギーｅ１とせん断発熱エネルギーｅ２の双方を考慮し、（３）式によりｋ＋１ステップ目の原料粒子の内部エネルギー増分を求める。ｋは時間を進めた回数、つまり一定時間毎に実行される解析の回数を示す。なお、（３）式におけるｉは着目粒子である。
【００３２】
【数３】
2014104696.tif 000005
【００３３】
伝熱エネルギーの増分Δｅ１は、非定常拡散方程式をＭＰＳ法のラプラシアンモデルで離散化すると（４）、（５）式にて計算できる。なお、（４）、（５）式における、ｊは着目粒子ｉを中心とした演算相互作用半径内にある粒子番号、ｔは時間、ｎ０は初期の粒子数密度、λＴｊｉは熱伝導率、ｒは粒子間距離、ｄは演算の次元、Ｔは温度、ｗは重み関数である。
【００３４】
【数４】
2014104696.tif 000006
【００３５】
【数５】
2014104696.tif 000007
【００３６】
せん断発熱エネルギーの増分Δｅ２は運動エネルギーの粘性消散と同等であるとの仮定の下で（６）式にて計算できる。なお、（６）式における、ｕは速度ベクトル、ηは粘性係数である。
【００３７】
【数６】
2014104696.tif 000008
【００３８】
演算部１１３におけるＤＥＭ法の計算では、固体状態の原料粒子が対象のため、粘性は有しない。そこで、粒子間の反発係数ｅを用いて粘性係数を（７）式にて与える。なお、（７）式における、ｍｓ＊は換算質量、πは円周率である。
【００３９】
【数７】
2014104696.tif 000009
【００４０】
演算部１１３は、これら式（２）〜（７）を用いて求めた原料粒子の内部エネルギー増分から、樹脂原料の比熱で換算することで原料粒子の温度上昇分を演算する。演算部１１３は、樹脂原料の初期設定温度に、この算出された温度上昇分を加算することでステップ数ｋ＋１ステップ目の原料粒子の樹脂温度を算出する。これにより、時間経過による原料粒子の温度変化が求められる。
【００４１】
一方、原料粒子の運動量に基づいて求めた樹脂温度が設定融点Ｔｍに達した場合、その粒子が液体に相変化したとみなし、演算部１１３は、原料粒子を液体とした運動量計算に移行する。この運動量計算にはＭＰＳ法またはＳＰＨ法を用いることが好ましい。液体の運動量は（８）式のナビエ・ストークス方程式で表される。なお、（８）式において、ρは密度、Ｐは圧力である。
【００４２】
【数８】
2014104696.tif 000010
【００４３】
ここで、ＭＰＳ法を用いる場合、演算部１１３は、粒子間相互作用を（９）式に基づいた重み関数ｗｉｊに基づき演算を行う。なお、（９）式において、ｒｅは粒子間相互作用の及ぶ範囲の影響半径、ｒｉｊは粒子間距離である。
【００４４】
【数９】
2014104696.tif 000011
【００４５】
一方、ＳＰＨ法を用いる場合、ＳＰＨ法はＭＰＳ法と同様に流体力学的知見から開発された演算手法であるため、（８）式を用いて運動量を演算する点はＭＰＳ法と共通である。ＭＰＳ法は非圧縮性流体を取り扱うことを前提に開発された手法であるが、ＳＰＨ法は圧縮性流体を対象として開発された点が大きく異なる。この考え方の違いは、重み関数ｗｉｊの定義付けに現れる。演算部１１３は、ＳＰＨ法において、（１０）式に示すカーネル関数ｗｉｊ（ｈ）を用いて計算を行う。なお、（１０）式において、ｈはカーネルの大きさを表すパラメータである。
【００４６】
【数１０】
2014104696.tif 000012
【００４７】
演算部１１３は、これら（９）及び（１０）式を用い、（８）式の勾配とラプラシアンの微分演算子に対して着目粒子ｉの離散化を行う。なお、具体的なＭＰＳ法、ＳＰＨ法での運動量計算は、非特許文献１あるいは非特許文献４に詳述されているため、その説明を省略する。
【００４８】
なお、演算部１１３は、ＭＰＳ法による原料粒子の温度上昇演算手法をＤＥＭ法と同様に、（２）〜（６）式にて算出する。なお、ＭＰＳ法では粘度の物性を有することから（６）式の粘性係数ηは（７）式のよる算出は行わず、予め設定された粘性係数そのものを適用する。この粘性係数は使用者により入力されてもよく、予め樹脂原料の種別毎にＨＤＤ１５に格納され、入力された樹脂原料に基づいて取得されるようにしてもよい。
【００４９】
実際に単一粒子の運動量演算をＤＥＭ法からＭＰＳ法へ転移させる場合、その粒子が固体か液体かにより（１）と（８）式をそれぞれ使い分けるには融点近傍における物理量の計算不安定性を招き解が発散する可能性が高くなる。そのため、演算部１１３は、ＤＥＭ−ＭＰＳ連続相変化演算として、（１）と（８）式を融合させた（１１）式の運動方程式を用いる。（１１）式において、νは動粘性係数である。
【００５０】
【数１１】
2014104696.tif 000013
【００５１】
前述した（１１）式の右辺第二項では、動粘性係数νが０以上、つまり粘性を有する場合の原料粒子にのみ計算を行うとし、第三項は粒子の直径に応じた質量計算が必要となるためＤＥＭ法における原料粒子のみの直径を対象とする計算を行うとした。これにより、固体粒子と流体粒子いずれの場合でも（１１）式にて運動量演算が行える。また、重力による外力の影響を考慮する場合には、ナビエ・ストークス方程式に従い重力の計算項を（１１）式右辺第四項に付与することが好ましい。
【００５２】
固体状態の原料粒子が融解し液体状態へ連続変化する相変化モデルでは、前述した設定融点Ｔｍを境にした急激な演算モデルの移行による物理量の特異点（例えば、動粘性係数が突然有限値を持つことによる不連続性によってもたらされる解の発散）が発生する危惧が生じる。この特異点問題を緩和、解消するために、演算部１１３は、設定融点Ｔｍ以下の規定温度である遷移開始温度Ｔ０から設定融点Ｔｍまでの温度範囲（Ｔｏｆｆｓｅｔ＝Ｔｍ−Ｔ０）で（１１）式における混在状態（ＤＥＭ−ＭＰＳ混在）の動粘性係数を、（１２）式を用いて算出する。また、固体状態（ＤＥＭ演算）の動粘性係数を（１３）式、液体状態（ＭＰＳ演算またはＳＰＨ演算）の動粘性係数を（１４）式を用いて算出する。なお、これら式におけるνｍｐｓはＭＰＳ法にて算出される動粘性係数、Ｔは粒子温度である。
【００５３】
【数１２】
2014104696.tif 000014
【００５４】
【数１３】
2014104696.tif 000015
【００５５】
【数１４】
2014104696.tif 000016
【００５６】
上式を採用すると、粒子温度Ｔが遷移開始温度Ｔ０を越え完全固体状態から混在状態に移行した状態では、（１２）式の右辺第二項の値が小さく、動粘性係数νは微小値をとる。一方、粒子温度Ｔが設定融点Ｔｍに上昇するにつれ、動粘性係数νはνｍｐｓに近くなり、設定融点Ｔｍ以降では真の動粘性係数を有した運動量演算が行われる。つまり、（１１）式及び（１２）式を採用することで相変化における動粘性係数の連続性を保つことが可能となるため、解の安定性を保持できることが可能となる。
【００５７】
オフセット温度Ｔｏｆｆｓｅｔ、設定融点Ｔｍ、遷移開始温度Ｔ０について、図３を用いて具体的に説明する。図３は、樹脂原料の比熱曲線と動粘性係数の関係を示す模式図である。一般的な結晶性樹脂の比熱は図３に示されるように、固体および液体状態では一定値と近似可能な数値をとるが、融点近傍では融点を極大値としたシャープな増大を示す。本実施の形態では、この比熱増大領域を相変化の遷移状態（混在状態）と定義し、遷移開始温度をＴ０と設定する。遷移開始温度Ｔ０より低い温度では、粒子は完全固体状態であるため、ＤＥＭ法で演算される。遷移開始温度Ｔ０から設定融点Ｔｍまでの温度領域ではＤＥＭ法とＭＰＳ法の解法が混在するが、粒子温度が設定融点Ｔｍに近くなるほどＭＰＳ法の重みが増すため動粘性係数が大きくなる。このことから判るように、遷移開始温度Ｔ０は解析対象とする樹脂の比熱曲線によって定められる、即ち、オフセット温度Ｔｏｆｆｓｅｔの範囲は、樹脂原料の温度変化による比熱の変化に基づいて決定することが好ましい。
【００５８】
なお、オフセット温度Ｔｏｆｆｓｅｔを狭く（低く）することにより、実在の現象に近い解析結果が得られるが、演算が不安定化する可能性がある。一方、オフセット温度Ｔｏｆｆｓｅｔを広く（高く）することにより、演算が安定化するが、実在の現象とはかけ離れた解析結果が得られる可能性が高くなる。このことから、オフセット温度Ｔｏｆｆｓｅｔは１〜２０℃が好適である。また、解析の対象となる樹脂原料の融点が高いものほど、オフセット温度Ｔｏｆｆｓｅｔを広くとることにより演算が安定化する。よって、オフセット温度Ｔｏｆｆｓｅｔは、対象とする樹脂原料により異なることは言うまでもなく、樹脂原料の種別毎に設定すべきオフセット温度Ｔｏｆｆｓｅｔが対応付けられたテーブルをＨＤＤ１５に格納するようにしてもよい。この場合、初期設定条件の項を設けて更に細分化してオフセット温度Ｔｏｆｆｓｅｔを設定することが好ましい。また、樹脂原料の比熱曲線を取得し、当該比熱曲線からその増大領域を割り出してオフセット温度Ｔｏｆｆｓｅｔおよび遷移開始温度Ｔ０を決定するようにしてもよい。また、設定融点Ｔｍとして必ずしも比熱が極大値を示す温度を与える必要はなく、その極大値以上、相変化の遷移状態終了以下の温度範囲で適宜設定すればよい。
【００５９】
次に、図面を用いて、シミュレーション装置１の動作の詳細を説明する。図４及び図５は、本実施の形態に係るシミュレーション装置の溶融可塑化工程解析処理を示すフローチャートである。図４に示されるように、先ず、シミュレーション装置１の溶融可塑化工程解析処理が実行され、入力部１３により使用者が前述した初期設定条件、シリンダ内の樹脂の配置情報等の情報を入力すると、取得部１１１はこれを取得する（Ｓ１０１）。取得後、判定部１１２は、原料粒子のうちの着目粒子iの粒子温度が設定融点Ｔｍ以上か否かを判定する（Ｓ１０２）。なお、ここで判定される粒子温度は、一定時間前の粒子温度が対象となる。例えば初回（解析処理を行っていない。ｋ＝０）であれば初期粒子温度が判定対象となり、２回目以降（ｋ＝ｋ＋１）であれば、前回のステップｋにおいて着目粒子ｉに対応付けられた、後述する温度上昇算出処理に基づいて算出される粒子温度が判定対象となる。
【００６０】
粒子温度が設定融点Ｔｍ以上でない場合（Ｓ１０２，ＮＯ）、演算部１１３は第１温度上昇算出処理を行う（Ｓ１０３）。この第１温度上昇算出処理は前述した（２）式〜（７）式が用いられる。これら式により着目粒子iの内部エネルギー増分が求められ、この内部エネルギー増分から樹脂の比熱で換算し、着目粒子ｉの温度上昇分が求められる。また、演算部１１３は、この温度上昇分と初期粒子温度とを足し合わせ、現在の着目粒子の温度を算出し、ｉと対応付ける。一方、粒子温度が設定融点Ｔｍ以上である場合（Ｓ１０２，ＹＥＳ）、演算部１１３は第２温度上昇算出処理を行う（Ｓ１０４）。この処理においては、前述した（２）式〜（６）式が用いられる。この処理結果から現在の着目粒子の温度を算出し、ｉと対応付ける処理はステップＳ１０３の第１温度上昇算出処理と同様である。
【００６１】
粒子温度算出後、判定部１１２は、着目粒子の粒子温度が遷移開始温度Ｔ０未満か、遷移開始温度Ｔ０〜設定融点Ｔｍ未満の範囲内か、設定融点Ｔｍ以上かの何れであるかの判定を行う（Ｓ１０５）。着目粒子の粒子温度が遷移開始温度Ｔ０未満である場合（Ｓ１０５，Ｔ０未満）、演算部１１３は、動粘性係数の算出式として（１３）式を選択、即ち動粘性係数を０として算出し（Ｓ１０６）、これを用いて第１運動量算出処理を行い、運動量を算出すると共に、着目粒子ｉの速度を算出する（Ｓ１０７）。この処理においては、前述した（１）式等が用いられる。
【００６２】
一方、着目粒子の粒子温度が遷移開始温度Ｔ０〜設定融点Ｔｍ未満の範囲内である場合（Ｓ１０５，Ｔ０以上〜Ｔｍ未満）、演算部１１３は、動粘性係数の算出式として（１２）式を選択し、この式に基づいて動粘性係数を算出し（Ｓ１０８）、これを用いて第２運動量算出処理を行い、運動量を算出すると共に、着目粒子ｉの速度を算出する（Ｓ１０９）。この処理においては、前述した（１１）式等が用いられる。
【００６３】
また、着目粒子の粒子温度が設定融点Ｔｍ以上である場合（Ｓ１０５，Ｔｍ以上）、演算部１１３は、動粘性係数の算出式として（１４）式を選択し、この式に基づいて動粘性係数を算出し（Ｓ１１０）、これを用いて第３運動量算出処理を行い、運動量を算出すると共に、着目粒子ｉの速度を算出する（Ｓ１１１）。この処理においては、初期設定条件の入力の際にＭＰＳ法が選択されていれば、前述した（８）式、（９）式等が用いられ、ＳＰＨ法が選択されていれば、（８）式、（１０）式等が用いられる。
【００６４】
運動量の算出後、演算部１１３は、着目粒子ｉの移動度（移動のしかた）及び圧力ポワソン方程式の計算を行う（Ｓ１１２）。この演算については従来の解析処理と同様、例えばＭＰＳ法における手法を用いた解析処理と同様であるため、その説明は省略する。図５に移り、演算後、判定部１１２は、着目粒子ｉの粒子温度が設定融点Ｔｍ以上であるか否かの判定を行う（Ｓ１１３）。粒子温度の粒子温度が設定融点Ｔｍ以上である場合（Ｓ１１３，ＹＥＳ）、演算部１１３は、圧力勾配項計算及び着目粒子ｉの速度、位置の修正を行い（Ｓ１１４）。この演算や修正処理についても従来の解析処理と同様、例えばＭＰＳ法における手法を用いた解析処理と同様であるため、その説明は省略する。
【００６５】
ステップＳ１１４の処理後または粒子温度が設定融点Ｔｍ以上でない場合（Ｓ１１３，ＮＯ）、判定部１１２は、全ての原料粒子を着目粒子として選択したか否かの判定を行う（Ｓ１１５）。全ての原料粒子を着目粒子として選択していない場合（Ｓ１１５，ＮＯ）、判定部１１２は、着目粒子ｉに１を加算した結果をｉとし、図４に示されるようにステップＳ１０２の着目粒子ｉの粒子温度の判定処理を行う。一方、全ての原料粒子を着目粒子として選択した場合（Ｓ１１５，ＹＥＳ）、判定部１１２は、前述した各演算に基づく解析結果を出力するか否かを判定する（Ｓ１１７）。この判定は、例えば、出力部１１４が表示部１４に解析結果を出力するか否かを示すＧＵＩを表示し、使用者に入力を促すことによる入力結果に基づいてもよく、規定の回数解析結果を得られたか否かにより判定してもよい。解析結果を出力する場合（Ｓ１１７，ＹＥＳ）、出力部１１４は、解析結果を表示部１４に出力し、表示させる（Ｓ１１８）。
【００６６】
ステップＳ１１８の処理後または解析結果を出力しない場合（Ｓ１１７，ＮＯ）、判定部１１２は、入力された終了条件を満たすか否かを判定する（Ｓ１１９）。この判定は、例えば、終了条件が一定時間の経過（３日間等）であれば、溶融可塑化工程解析処理の開始時間と現在の時間とから経過時間を算出し、一定時間を過ぎているか否かを判定することがある。また、規定回数解析を終えたか否か（解析結果が規定数求められたか否か）等、終了条件については適宜使用者が設定できる。終了条件を満たす場合（Ｓ１１９，ＹＥＳ）、本フローは終了となる。一方、終了条件を満たさない場合（Ｓ１１９，ＮＯ）、判定部１１２は、時間ｔに所定の値（例えば１０−４ｓ）を加算すると共に、ステップ数ｋに１を加算して（ｋ＝ｋ＋１）とし、図４に示されるようにステップＳ１０２の着目粒子ｉの粒子温度の判定処理を行う。
【００６７】
本実施の形態においては、原料粒子が固体状態であればＤＥＭ法、液体状態であればＭＰＳ法またはＳＰＨ法を用いると説明したが、これに限定するものではない。固体状態にある原料粒子の運動量及び速度を演算できる粒子法、液対状態にある原料粒子の運動量及び速度を演算できる粒子法であればよい。
【００６８】
また、本実施の形態においては、着目粒子ｉの温度上昇分、運動量及び速度を演算した後に、次の粒子に移行して演算を行うと説明したが、これに限定されるものではなく、各種演算毎に着目粒子を変えて演算するようにしてもよい。
【００６９】
本実施の形態によれば、原料粒子が固体状態であれば固体状態を考慮したＤＥＭ法等の粒子法を用いて運動量を演算し、原料粒子が液体状態であれば液体状態を考慮したＭＰＳ法やＳＰＨ法等の粒子法に切り替えて運動量を演算することにより、演算が安定化すると共に良好に樹脂原料の溶融可塑化工程を予測解析することができる。また、固体状態から液体状態へと遷移する遷移領域（比熱増大領域）においてはＤＥＭ−ＭＰＳ連続相変化演算として、（１）と（８）式を融合させた（１１）式の運動方程式を用いることにより、より演算が安定化し、より良好な樹脂原料の溶融可塑化工程を予測解析することができる。
【００７０】
（実施例）
溶融可塑化工程解析処理における初期設定条件として、粒子直径３ｍｍの高密度ポリエチレン（ＨＤＰＥ）を対象とし、初期温度５０℃、総数１４０８個のＨＤＰＥの樹脂原料をスクリュとシリンダで構成された加熱装置に配置させた。スクリュとシリンダは静止状態にてそれぞれ１５０℃に加熱された状態で、オフセット温度Ｔｏｆｆｓｅｔを、０℃、１０℃、３０℃に設定した３パターンの解析を開始し、時間経過による原料粒子の温度上昇の予測を実施した。ＨＤＰＥの物性は、融点１３０℃、比熱と熱伝導率は固体・溶融体ともにそれぞれ２８００Ｊ／ｋｇＫと０．３３５Ｗ／ｍＫとした。
【００７１】
図６は、オフセット温度を変化させた場合の溶融可塑化工程解析処理結果を示す図である。図６に示されるように、オフセット温度Ｔｏｆｆｓｅｔが０℃では安定した演算が行えず、演算途中で終了し、１０℃、３０℃では演算は安定する結果となった。即ち、オフセット温度Ｔｏｆｆｓｅｔの値を大きくすることで演算の安定化は可能となった。しかしながら、オフセット温度Ｔｏｆｆｓｅｔを３０℃と広くすると溶融可塑化が遅延する傾向が目立ち、４５０秒後の長時間後でも完全に溶融しなかった。一方、オフセット温度Ｔｏｆｆｓｅｔを１０℃に設定することで、定性的に妥当な解を得ることができた。
【００７２】
以下、オフセット温度Ｔｏｆｆｓｅｔを１０℃に設定した場合の解析結果を図７〜図１１を用いて説明する。図７は樹脂原料をシリンダに配置させた状態を可視化した図であり、図８は解析開始後１００秒後における樹脂原料の状態を可視化した解析結果を示す図である。これら図に示される３１はスクリュ、３２はシリンダ、４は樹脂原料を示す。樹脂原料４の原料粒子は全てＤＥＭ法で演算が開始され、スクリュ３１とシリンダ３２表面の流体温度が上昇することでＤＥＭ−ＭＰＳ混在状態をとっている。また、スクリュ３１とシリンダ３２間の流路寸法が狭い領域に存在する原料粒子は加熱エネルギーを多く受けるため融点以上となりＭＰＳモデルに移行している。
【００７３】
図９は解析開始後１５０秒後における樹脂原料の状態を可視化した解析結果を示す図であり、図１０は解析開始後２００秒後における樹脂原料の状態を可視化した解析結果を示す図である。スクリュ３１、シリンダ３２表面から距離のある流路は伝熱が遅延される傾向を良好に予測できており、時間経過とともに徐々に加熱エネルギーがその粒子へも伝わり、次第に溶融状態に移行する妥当な予測が行えている。図１１は、解析開始後３００秒後における樹脂原料の状態を可視化した解析結果を示す図である。全ての粒子がＭＰＳモデルへと移行しており、樹脂原料４は全て溶融状態に相変化した。これら演算結果はＨＤＰＥの溶融可塑化工程を定性的に妥当な予測が行えており、シミュレーション装置１の溶融可塑化工程解析処理が不連続体の固体ペレットである樹脂原料の溶融可塑化工程を良好に予測できることを示したものといえる。
【００７４】
本発明は、その要旨または主要な特徴から逸脱することなく、他の様々な形で実施することができる。そのため、前述の実施の形態は、あらゆる点で単なる例示に過ぎず、限定的に解釈してはならない。本発明の範囲は、特許請求の範囲によって示すものであって、明細書本文には、何ら拘束されない。更に、特許請求の範囲の均等範囲に属する全ての変形、様々な改良、代替および改質は、全て本発明の範囲内のものである。
【００７５】
また、シミュレーションプログラム１０１を、図１２に示されるような、コンピュータにより読み取り可能な可搬型の記録媒体８に記憶させ、当該記録媒体８を情報処理装置９等に読み込ませることにより、前述した機能を情報処理装置９等に実現させることができる。記録媒体８としては、例えば、光ディスク（ＣＤ−ＲＯＭ、ＤＶＤディスク等）、磁気ディスク（ハードディスクドライブ等）、フラッシュメモリ、ＩＣカード、更にネットワークを介することで伝送可能な媒体等、コンピュータで読み取りや実行が可能な全ての媒体が含まれる。
【００７６】
なお、特許請求の範囲に記載の第１の算出手段、第２の算出手段は、例えば前述の実施の形態における演算部１１３に対応する。また、第１の粒子法は、例えばＤＥＭ法に対応し、第２の粒子法は、例えばＭＰＳ法、ＳＰＨ法に対応する。予め設定された融点は、例えば設定融点Ｔｍに対応し、遷移開始点は遷移開始温度Ｔ０に対応する。固液混在状態は、例えば遷移状態、混在状態に対応する。
【符号の説明】
【００７７】
１  シミュレーション装置、４  樹脂原料、８  記録媒体、９  情報処理装置、１１  ＣＰＵ、１２  記憶部、１３  入力部、１４  表示部、１５  ＨＤＤ、３１  スクリュ、３２  シリンダ、１０１  シミュレーションプログラム、１１１  取得部、１１２  判定部、１１３  演算部、１１４  出力部。
【図１】
2014104696.tif 000020
【図２】
2014104696.tif 000021
【図３】
2014104696.tif 000022
【図４】
2014104696.tif 000023
【図５】
2014104696.tif 000024
【図６】
2014104696.tif 000025
【図１２】
2014104696.tif 000026
【図７】
2014104696.tif 000027
【図８】
2014104696.tif 000028
【図９】
2014104696.tif 000029
【図１０】
2014104696.tif 000030
【図１１】
2014104696.tif 000031