(57)【特許請求の範囲】
【請求項１】
バンバリーミキサーのケーシングと、前記ケーシング内を回転する少なくとも１本のロータとの間で形成されるバンバリーミキサーの混練空間であるチャンバー内で、前記ロータによって混練されるゴム又は樹脂材料である可塑性材料が混練されるときの流動状態を、コンピュータを用いて計算する可塑性材料のシミュレーション方法であって、
前記チャンバーを有限個の要素でモデル化したチャンバーモデルを設定するステップ、
前記可塑性材料をモデル化した材料モデルを設定するステップ、
空気をモデル化した気相モデルを設定するステップ、及び
前記チャンバーモデル内に、前記材料モデルと前記気相モデルとを混在して配置しかつ予め定めた条件に基づいて流動計算を行うステップを含み、
前記材料モデルは、全温度領域で流体として扱われ、しかも、前記材料モデルには、前記可塑性材料の粘度が定義され、
前記気相モデルには、実際の空気の粘度の５〜１０倍の粘度が定義され、
前記流動計算を行うステップは、前記チャンバーモデルの中で、前記気相モデル及び前記材料モデルの２相を平均化して一つの流体モデルとして扱うＶＯＦ法が用いられ、
前記流体モデルの密度及び粘度は、前記材料モデルの相及び前記気相モデルの相の各相がそれぞれ占める体積で重み付けされて平均化されたものが使用され、
前記流動計算では、せん断発熱については、体積分率が０．９０以上の材料モデルの要素についてのみ行い、前記気相モデルのせん断発熱の計算を行わないことを特徴とする可塑性材料のシミュレーション方法。
【請求項２】
前記可塑性材料モデルと前記チャンバーモデルとの接触面には、スリップが定義される請求項１記載の可塑性材料のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、架橋前のゴムや樹脂等の可塑性材料の流動状態を解析するのに役立つシミュレーション方法に関する。
【背景技術】
【０００２】
近年、コンピュータを用いた数値シミュレーションが種々提案されており、例えば、架橋前のゴムや樹脂といった可塑性材料が、バンバリーミキサー内で混練されるときの流動状態などを計算するものが提案されている。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】V. Collin (1)*, E. Peuvrel-Disdier (1) ら「Numerical and Experimental Study of Dispersive Mixing of Agglomerates」
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ところで、上記非特許文献１のシミュレーションでは、可塑性材料が移動する空間であるチャンバー内が１００％可塑性材料で満たされる状態、つまりチャンバーの容積に対する可塑性材料の充填率が１００％であることが前提となる。
【０００５】
しかしながら、バンバリーミキサー等において、充填率が１００％であることは通常あり得ないし、充填率が変わることにより可塑性材料の練りの状態が変わる事も経験的に知られている。そうであるにも関わらず、上述の非特許文献１の方法では、可塑性材料の充填率を１００％以外として計算することができないという問題があった。特に、実機を使用しても、充填率１００％の実験は困難であり、シミュレーション結果を検証することもできない。
【０００６】
他方、可塑性材料と空気とは、粘度の値が著しく相違するため、このような２相を混在した流動計算のシミュレーションを行うと、計算が不安定化するという問題がある。
【０００７】
本発明は、以上のような問題点に鑑み案出なされたもので、チャンバーモデル内に、可塑性材料をモデル化した材料モデルと、空気をモデル化した気相モデルとを混在して配置して流動計算を行うとともに、材料モデルには、前記可塑性材料の粘度を定義する一方、前記気相モデルには、実際の空気よりも大きい粘度を定義することを基本として、計算時間の著しい増加や計算の不安定化を抑えつつ、可塑性材料の充填率を変えたシミュレーションを可能としうる可塑性材料のシミュレーション方法を提供することを目的としている。
【課題を解決するための手段】
【０００８】
本発明のうち請求項１記載の発明は、バンバリーミキサーのケーシングと、前記ケーシング内を回転する少なくとも１本のロータとの間で形成されるバンバリーミキサーの混練空間であるチャンバー内で、前記ロータによって混練されるゴム又は樹脂材料である可塑性材料が混練されるときの流動状態を、コンピュータを用いて計算する可塑性材料のシミュレーション方法であって、前記チャンバーを有限個の要素でモデル化したチャンバーモデルを設定するステップ、前記可塑性材料をモデル化した材料モデルを設定するステップ、空気をモデル化した気相モデルを設定するステップ、及び前記チャンバーモデル内に、前記材料モデルと前記気相モデルとを混在して配置しかつ予め定めた条件に基づいて流動計算を行うステップを含み、前記材料モデルは、全温度領域で流体として扱われ、しかも、前記材料モデルには、前記可塑性材料の粘度が定義され、前記気相モデルには、実際の空気の粘度の５〜１０倍の粘度が定義され、前記流動計算を行うステップは、前記チャンバーモデルの中で、前記気相モデル及び前記材料モデルの２相を平均化して一つの流体モデルとして扱うＶＯＦ法が用いられ、前記流体モデルの密度及び粘度は、前記材料モデルの相及び前記気相モデルの相の各相がそれぞれ占める体積で重み付けされて平均化されたものが使用され、前記流動計算では、せん断発熱については、体積分率が０．９０以上の材料モデルの要素についてのみ行い、前記気相モデルのせん断発熱の計算を行わないことを特徴とする。
【００１１】
また請求項２記載の発明は、前記可塑性材料モデルと前記チャンバーモデルとの接触面には、スリップが定義される請求項１記載の可塑性材料のシミュレーション方法である。
【発明の効果】
【００１２】
本発明のシミュレーション方法は、チャンバー内の可塑性材料の流動状態を、コンピュータを用いて計算する可塑性材料のシミュレーション方法であって、前記チャンバーを有限個の要素でモデル化したチャンバーモデルを設定するステップ、前記可塑性材料をモデル化した材料モデルを設定するステップ、空気をモデル化した気相モデルを設定するステップ、及び前記チャンバーモデル内に、前記材料モデルと前記気相モデルとを混在して配置して流動計算を行うステップを含み、前記材料モデルには、前記可塑性材料の粘度が定義される一方、前記気相モデルには、実際の空気よりも大きい粘度が定義されることを特徴としている。
【００１３】
このようなシミュレーション方法では、チャンバーモデル内に、材料モデルと気相モデルとを混在させているため、チャンバー内への可塑性材料の充填率が１００％を下回る状態についても流動計算のシミュレーションを行うことができる。
【００１４】
また、材料モデルには、前記可塑性材料の粘度が定義される一方、前記気相モデルには、実際の空気よりも大きい粘度が定義される。空気を、その本来の粘度を用いて通常の計算を行った場合、非常に多くの計算コストを必要とする他、計算が不安定となり易いが、気相モデルの粘度を、本来の値よりも大きくすること、とりわけ計算結果に影響を与えない範囲内で大きくすることによって、高速にかつ安定して計算することができる。
【図面の簡単な説明】
【００１５】
【図１】可塑性材料を混練するバンバリーミキサーの部分断面図である。
【図２】本実施形態の処理手順の一例を示すフローチャートである。
【図３】チャンバーモデルの斜視図である。
【図４】チャンバーモデルの断面図である。
【図５】チャンバーモデルを分解して示す断面図である。
【図６】チャンバーモデル内に材料モデルと気相モデルとを混在して配置した状態を示す断面図である。
【図７】流動計算の処理手順の一例を示すフローチャートである。
【図８】混練シミュレーションをした後の状態を視覚化した線図である。
【発明を実施するための形態】
【００１６】
以下、本発明の実施の一形態が図面に基づき説明される。
本発明は、コンピュータ（図示省略）を使用して、チャンバー内の可塑性材料の流動状態をシミュレートする可塑性材料のシミュレーション方法である。
【００１７】
前記可塑性材料としては、架橋前のゴム、樹脂又はエラストマー等の可塑性材料が挙げられる。本発明でいう可塑性材料は、粘度を有するものであって、十分に練られて安定的な流動状態（流体）とみなすことができる状態のものを前提としている。このような状態としては、ゴムの場合、十分に練られて約８０℃程度まで昇温した状態が相当する。
【００１８】
前記チャンバーとは、上述の可塑性材料が配されかつ移動する閉じた空間を意味する。チャンバーには、種々の形状のものが含まれる。代表的なチャンバーとして、バンバリーミキサーが有する混練空間などを挙げることができる。図１に示されるように、本実施形態のチャンバー４は、バンバリーミキサー１のケーシング２と、該ケーシング２内を回転する一対のロータ３、３との間で形成される断面横向きの略８の字状の混練空間とされる。ただし、チャンバーは、このような形状に限定して解釈されるものではない。
【００１９】
図２には、本実施形態の処理手順の一例が示されており、以下、各ステップが順に説明される。
【００２０】
［チャンバーモデル］
本実施形態のシミュレーション方法では、先ず、前記チャンバー４を有限個の要素でモデル化したチャンバーモデル５が前記コンピュータに入力、設定される（ステップＳ１）。
【００２１】
図３には、チャンバーモデル５を視覚化した斜視図が、図４には、その断面図がそれぞれ示されている。チャンバーモデル５は、前記ケーシング２の内周面をなす外周面５ｏと、前記回転する一対のロータ３、３の外周面である内周面５ｉと、ロータ３の軸方向の両端側で前記外周面を閉じる両端面５ｓとで閉じられた三次元空間を有し、この空間が三次元の要素ｅによって分割されている。
【００２２】
チャンバーモデル５は、オイラー要素ｅで分割（離散化）されている。要素分割は、四面体、六面体などの他、多面体セル（ポリヘドラルグリッド）といった三次元要素で行われる。そして、各要素について、可塑性材料（材料モデル）の圧力、温度及び／又は速度等の物理量が計算される。
【００２３】
また、チャンバーモデル５は、前記外周面５ｏ及び両端面５ｓは変形しないが、内周面５ｉは、ロータ３の回転によって回転するため、容積形状が変化する。本実施形態では、このような形状変化をなすチャンバーモデル５は、図５に分離させて示されるように、一対の回転部５Ａ、５Ｂと、これらの間を継ぐ継ぎ部５Ｃと、これらが収容される外枠部５Ｄとの４つの部分に分けて構成される。
【００２４】
前記回転部５Ａ、５Ｂは、各々、円形の外周面５Ａｏ、５Ｂｏと、ロータ３の外周面に等しい前記内周面５ｉを有する筒状でモデル化されている。そして、回転部５Ａ、５Ｂは、各々、その中心Ｏａ、Ｏｂの周りで回転することにより、ロータ３の回転に伴うチャンバー３の容積形状の変化を表現することができる。
【００２５】
前記継ぎ部５Ｃは、一対の回転部５Ａ、５Ｂ間に配されて静止している。また、継ぎ部５Ｃは、各回転部５Ａ、５Ｂと接触する凹円弧面ｊを有している。継ぎ部５Ｃの凹円弧面ｊと回転部５Ａ、５Ｂの外周面５Ａｏ、５Ｂｏとは、スライディングサーフェース等の境界条件が定義され、これにより、チャンバーモデル５の回転部５Ａ、５Ｂ内で生じる物理的な作用（力及び熱等）が、この凹円弧面ｊを介して継ぎ部５Ｃの中にある材料モデルとｂへと伝達される。
【００２６】
前記外枠部５Ｄは、前記回転部５Ａ、５Ｂ及び継ぎ部５Ｃを囲む筒状をなし、その軸方向両端は、前記両端面５ｓによって閉じられている。また、外枠部５Ｄも、回転部５Ａ、５Ｂ及び継ぎ部５Ｃの要素との接触面に、スライディングサーフェース等の境界条件が定義される。これにより、外枠部５Ｄ内で回転する回転部５Ａ、５Ｂ内で生じる物理的な作用（力及び熱等）が、両者の接触面を介して外枠部５Ｄにも伝えられる。なお、外枠部５Ｄは、ロータの作用によって大きなせん断が生じるため、材料の速度等をより詳細に計算すべく、回転部５Ａ、５Ｂ及び継ぎ部５Ｃよりも小さい要素で構成される。これによって、チャンバーモデル５の外周面付近の材料モデルの速度プロファイルなどをより詳細に計算するのに役立つ。
【００２７】
［材料モデル］
次に、本実施形態では、材料モデルが設定され、コンピュータに入力される（ステップＳ２）。該材料モデルは、前記チャンバー４内を流動する可塑性材料をモデル化したもので、解析対象となる可塑性材料に応じてせん断粘度、比熱、熱伝導率及び比重といった物性が定義されコンピュータに入力される。
【００２８】
前記せん断粘度は、例えば、解析対象となる可塑性材料から粘弾性特性（Ｇ'及びＧ”）が複数の温度条件で測定され、Cox-Merz則などを用いてせん断粘度に変換することで得られる。このようにして得られたせん断粘度は、例えば下記式（１）のべき乗法則で近似することができる。
η＝ｍγ'n-1  …（１）
ここで、η：せん断粘度、ｍ：係数、γ'：せん断速度、ｎ：係数である。
【００２９】
前記比熱は、解析対象の可塑性材料から、例えば断熱型連続法（＠２５℃）にて測定され、その値が前記コンピュータに入力される。
【００３０】
前記熱伝導率は、解析対象の可塑性材料から、例えば熱線法（＠２５℃）にて測定され、その値が前記コンピュータに入力される。
【００３１】
［気相モデル］
次に、気相モデルが設定され、コンピュータに入力される（ステップＳ３）。気相モデルは、チャンバー内に存在する空気をモデル化したものである。このような気相モデルをシミュレーションに取り入れることで、可塑性材料の充填率が１００％未満の状態の流動計算が可能になる。
【００３２】
本実施形態において、気相モデルには、粘度及び比重が定義される。気相モデルの比重は、空気の比重に基づいてそれと同一の値が決定される。しかしながら、本発明において、気相モデルの粘度は、実際の空気の粘度よりも大きい値が定義され、コンピュータに入力される。
【００３３】
気相モデル（空気）と材料モデル（可塑性材料）との混相流解析においては、通常では、両者の粘度差が非常に大きい。このため、後述する流動計算において、気相モデルと材料モデルとの界面で両者の速度差が著しく大きくなる。この結果、速度差に伴うせん断発熱を考慮した計算では、気相モデルと材料モデルとの界面におけるせん断発熱が異常に大きくなって不安定化しやすい。一方、可塑性材料の充填率を１００％未満とするためには、チャンバーモデル５内に可塑性材料と空気とを混在させた混相流解析が必要になる。
【００３４】
そこで、本発明では、計算結果に悪影響を及ぼさない範囲で、気相モデルの粘度を極力大きくすることによって、計算を安定させる。また、例えば、シミュレーションでのタイムステップ（時間増分）を大きく設定することで、計算速度の向上を図ることもできる。
【００３５】
発明者らは、気相モデルの粘度のみを変えて複数種類のシミュレーションを行い、チャンバーモデル５内の圧力場を比較したところ、気相モデルの粘度が実際の空気の粘度の１０倍を超えると、チャンバーモデル５内の圧力値が大きくなり、圧力場に悪影響が出ることを知見した。他方、気相モデルの粘度が、実際の空気の粘度の５倍未満であると、計算を安定させることができない。従って、「計算に悪影響を及ぼさない範囲」として、気相モデルの粘度は、実際の空気の粘度の５〜１０倍の範囲で定められるのが望ましい。
【００３６】
［境界条件］
次に、本実施形態では、流動計算のシミュレーションに必要な境界条件等、各種の条件が定義される（ステップＳ４）。設定される境界条件としては、チャンバーモデル５の表面での流速境界条件及び温度境界条件が挙げられる。
【００３７】
前記流速境界条件としては、下記の２通りが設定可能であり、シミュレーションの用途や精度等に応じて、いずれかが採用される。
ａ）壁面ノースリップ条件
ｂ）壁面スリップ条件
【００３８】
前記壁面ノースリップ条件では、材料モデルは、チャンバーモデル５の表面での流速は常に零とされる。一方、前記壁面スリップ条件では、材料モデルは、チャンバーモデル５の表面において流速を持つ。この場合、材料モデルとチャンバーモデル５との接触面のスリップ現象は、以下のモデル式（２）（Navier's Law）を適用することができる。
【数１】
0005514236.tif 000002
【００３９】
ここで、上記式の符号は次の通りである。
ｆs：チャンバーモデル５の表面での材料モデルのせん断応力
ｖslip：チャンバーモデル５の表面での材料モデルのスリップ速度（接線速度）
ｖwall：チャンバーモデル５の表面での材料モデルの接線速度
Ｆslip、eslip：物性値
（ただし、eslip＝１のときは線形、０＜eslip＜１のときはべき乗法則になる）。
【００４０】
また、前記チャンバーモデル５の表面での材料モデルのスリップ速度ｖslipを設定するには、先ず、初期の速度場を求め、次に、チャンバーモデル５の表面での材料モデルのせん断応力をその速度場（又は修正後の速度場）から求め、しかる後、チャンバーモデル５の表面での材料モデルのスリップ速度が上記式（２）で求められる。この際、急激なスリップ速度の増加を抑制するために、必要に応じて、１未満の緩和係数などを前記スリップ速度に乗じることで、チャンバーモデル５の各表面での壁面スリップ速度を、徐々に従来のSIMPLE法のループの中で真の値に収束させることが望ましい（後述）。
【００４１】
前記温度境界条件としては、次の２種類が設定可能である。
ａ）断熱条件：温度の計算を単純化するため、チャンバーモデル５の表面において、熱が外に逃げない条件。
ｂ）全てのチャンバーモデル５の表面温度が温調温度（例えば５０℃）に設定される条件。
【００４２】
他の条件としては、材料モデルの初期温度（本実施形態では２０℃（２９３Ｋ））、ロータの回転数（これは、チャンバーモデルの回転部５Ａ、５Ｂの回転数によって表現される）、チャンバーモデル５の表面のスリップ率、チャンバーモデル５の容積に対する可塑性材料の充填率などが挙げられる。さらに、他の条件としては、流動計算（シミュレーション）の初期状態、タイムステップ、内部処理でのイタレーションの反復回数、計算終了時刻などがある。初期状態は、例えば、図６に示されるように、チャンバーモデル５を横切る水平な境界面Ｓを基準として、それよりも上部を気相モデルの領域Ａとし、それよりも下部を材料モデルの領域Ｍとして混在配置できる。また、境界面Ｓのレベルを変えることにより、可塑性材料（材料モデル）の充填率が調節される。これらの条件は、シミュレーションの目的等に応じて任意に定められる。
【００４３】
［流動計算］
次に、本実施形態では、チャンバーモデル５内に、材料モデルと気相モデルとを混在して配置しかつ上記にて定めた条件に基づいて流動計算が行われる（ステップＳ５）。流動計算では、流体の運動状態を特定する３方向（ｘ，ｙ，ｚ）の速度成分と、流体の内部状態を特定する未知量である圧力ｐ及び温度Tとが計算される。つまり、解くべき未知数はこの５つである。なお、本実施形態では、非圧縮性流れの場合のNavier-Stoks方程式とし、気相モデル及び材料モデルの各密度を一定としている。
【００４４】
本実施形態において、材料モデルは、全温度領域で流体として扱われる。このため、流体の方程式（Navier-Stoks 方程式、質量保存式、エネルギー方程式の連立）を解くことになる。また、本実施形態では、空気と可塑性材料との２つの流体を一度に扱う必要があるため、自由界面の流れの計算で用いられるＶＯＦ（Volume of Fluid）法が用いられる。ＶＯＦ法では、二流体の界面の移動を直接計算するのではなく、各要素（「セル」ということもある。）の体積中の流体の充填率（体積分率）を定義して自由界面を表現するものである。なお、支配方程式は次の通りである。
【００４５】
［運動方程式］
本実施形態では、チャンバーモデル５の中に、気相モデル及び材料モデルの２相が配置される混相流モデルが一流体モデルとして扱われる。この場合、運動方程式は、下記式（３）の１組をそれぞれｘ，ｙ，ｚの３方向について解くだけで良い。これは、ＶＯＦ法により、気相と液相とを平均化して、一つの流体として扱うことで実現されている。
【数２】
0005514236.tif 000003
上記式において、符号は次の通りである。
ｕ：混相流モデルの速度
ｐ：混相流モデルの圧力
ρ：混相流モデルの密度
ｇ：重力加速度
Ｔ：混相流モデルの絶対温度
Ｆ：外力
【００４６】
また、材料モデルと気相モデルとの界面付近の密度ρ及び粘性係数μ等は、式（４）のように、材料モデルの相及び気相モデルの相の各相がそれぞれ占める体積で重み付けされて平均化されたものが使用される。
【数３】
0005514236.tif 000004
上記式において、符号は次の通りである。
αq:体積分率
ρq：各セルにおける各相の密度
μq：各セルにおける各相の粘度
【００４７】
［質量保存式］
同様に、質量保存（連続の式）及び圧力方程式も一組のみ解くことで足りる。つまり、本実施形態のシミュレーション方法では、混相でありながら、流れ場の計算としては単相と同じであり、場所（体積分率）によって、物性値が異なる流れを解いていることになる。各相の位置は、計算結果として得られる体積分率の分布（Volume Fluxion）より判断できる。
【００４８】
［エネルギー方程式］
エネルギー方程式については、下記式（５）で表現される。また、材料モデルの温度はこの式によって求めることができる。
【数４】
0005514236.tif 000005
上記式において、符号は次の通りである。
Ｅ：エンタルピ
ｋ：熱伝導率
Ｓ：ソース項
【００４９】
［体積分率輸送方程式］
体積分率の分布は、気相モデルと材料モデルとの２相間の界面位置を決定するものである。この体積分率αqは、下記式（６）を精度良く解くことで得られる。
【数５】
0005514236.tif 000006
ここで、チャンバーモデル５の任意の要素ｅにおいて、体積分率αq＝０の場合、当該要素にはｑ相（ｑ相：材料モデルの相）は存在しないことを示す。また、体積分率αq＝１の場合、当該要素はｑ相（材料モデルの相）で全部が満たされていることを示す。さらに、０＜αq＜１の場合、当該要素は、ｑ相と他の相（気相モデルの相）とで満たされていること、つまり混相であって両者の界面が存在することを示す。この方程式は、例えば、Modified-HRIC (Implicit)  で解くことができる。詳細は、ANSYS Fluent User’s Manual, 26.2.9 Modified HRIC Schemeに詳しく述べられている。
【００５０】
また、本実施形態において、上記各方程式は、圧力ベースの分離型解法で解かれる。圧力方程式と運動方程式とのカップリングには、例えばSIMPLE（Semi-Implicit Method for Pressure-Linked Equations）アルゴリズムが用いられるのが望ましい。
【００５１】
図７には、コンピュータが行う流動計算の具体的な処理手順の一例が示されている。
本実施形態では、先ず、速度と圧力勾配の下限及び上限であるリミットを設定し、速度勾配及び圧力勾配が計算される（ステップＳ５１）。
【００５２】
次に、最新の圧力場から離散化した前記運動方程式をセットアップし、反復法ソルバーを用いて、その解、即ち、材料モデル（又は混相）の３方向の速度が計算される（ステップＳ５２）。反復法としては、例えばガウス・ザイデル法等が挙げられる。
【００５３】
次に、上記速度が質量保存式を満たすか否かをチエックするために、先ず、チャンバーモデル５のセル表面における材料モデルの修正前の質量流量を求める（ステップＳ５３）。質量流量は、質量保存則を満たすか否かを判定するために必要である。また、「修正前の質量流量」とあるのは、SIMPLEアルゴリズムのループの中で、最初に取り敢えず使用される質量流量であり、誤差が大きいので、「修正前」と表現されている。
【００５４】
次に、速度場と圧力場とを、前記SIMPLEアルゴリズムを用いてカップリングし、圧力場を修正するための下記のような圧力補正方程式が構築される（ステップＳ５４）。
▽［ｋ▽φ］＝src
【００５５】
次に、上記ステップで得られた圧力補正方程式が、例えばＡＭＧソルバー、ＣＧ又はＢｉ−ＣＧ等の反復法で計算され、圧力補正量ｐ’が計算される（ステップＳ５５）。
【００５６】
次に、計算された解を元に、圧力場が修正される（ステップＳ５６）。圧力場の補正は、次のように行われる。
ｐn+1＝ｐn＋ωｐ’
ここで、ｐは圧力、ｎは現在のステップ数、ωは緩和係数であり、本実施形態では０．３が採用されている。
【００５７】
次に、境界面の境界条件が修正（アップデート）される。即ち、修正された圧力場から圧力勾配を求め、その圧力勾配を境界条件として与える（ステップＳ５７）。
【００５８】
次に、セル表面の質量流量が修正される（ステップＳ５８）。修正は次のように行われる。
ｍfn+1＝ｍf*＋ｍ'f
ここで、ｍf n+1は修正後のセル表面の質量流量、ｍf*は修正前のセル表面の質量流量、ｍ'fは質量流量の修正値である。
【００５９】
次に、速度場が修正される（ステップＳ５９）。速度場の修正は、例えば、次のようにして行われる。
ｖn+1＝ｖ*−（Ｖ▽ｐ'／∂pV）
ここで、Vはセルの体積、ｖ*は中間的な速度場（運動方程式から求めた補正前のもの）、∂pVは運動方程式のマトリクスの対角成分、▽ｐ'は圧力補正量の勾配である。
【００６０】
次に、式（５）のエネルギー方程式を解くことにより、材料モデルの温度及び粘度が算出される（ステップＳ６０）。
【００６１】
さらに、チャンバーモデル５の壁面スリップ速度ｖslipが下記式（７）で計算される（ステップＳ６１）。そして、これらの解が収束したか否かが判断される（ステップＳ６２）。なお、この収束判定条件は、予め定義されているものとする。
【数６】
0005514236.tif 000007
【００６２】
ここで、ｖslipは、先に述べたように、チャンバーモデル５の表面での材料モデルのスリップ速度、τwは壁面せん断応力、ｖwallはチャンバーモデル５の表面での材料モデルの接線速度、Ｆslipは、スリップのしやすさを示す物性値（係数）であり、値が大きいと壁面ノースリップ条件となる。本実施形態ではｖwall＝零であるが、ロータは回転しているので、ｖwall＝ｖ＝ｒ・ωで計算される速度を持っている。
【００６３】
そして、収束したと判断された場合（ステップＳ６２でＹ）、１タイムステップ分、次の時刻に移動する（ステップＳ６３）。前記収束条件は、例えば、修正された質量流量の総量が、予め定めた誤差の範囲内に留まったか否かで行いうる。そして、その時刻が指定された時刻に至っていない場合（ステップＳ６４でＮ）、再び、ステップＳ５１以降が繰り返される。他方、ステップＳ６２で収束していないと判断された場合には、ステップＳ５１以降が繰り返される。
【００６４】
また、本実施形態の流動計算では、せん断発熱については、材料モデルについてのみ計算が行われ、気相モデルのせん断発熱の計算は行われない。材料モデルについては、可塑性材料であるため、せん断発熱を考慮する必要があるが、材料モデルと気相モデルとの界面では、せん断速度が大きくなるケースが多い。ここで、せん断発熱は、「せん断粘度」に「せん断速度」の二乗を乗じたものであるから、せん断速度の異常に大きな値は、温度の計算（エネルギー方程式）に悪影響を与えるおそれがある。そこで、本実施形態では、せん断発熱の計算を材料モデルのみに限定、より具体的には、チャンバーモデル５の要素において、材料モデルの体積分率αqの値が一定値以上（本実施形態では０．９０以上）の要素に限定することで計算の安定化が図られている。
【００６５】
また、本実施形態では、前記材料モデルと前記チャンバーモデル５の表面との接触面に、壁面スリップ条件が定義されている。本実施形態では、材料モデルに少しずつスリップ速度が与えられる。この手順としては、コンピュータは、最新の圧力場を元にセットアップされた運動方程式の解として得られるせん断応力から、与えるべき壁面スリップ速度ｖslipが前記式（７）で計算される（上記ステップＳ６１）。
【００６６】
次に、ここで、計算されたスリップ速度ｖslipをいきなり接触面に与えると、式（２）’で表される壁面せん断応力τwが異常値を取っている場合等、壁面に異常なスリップ速度を与えるおそれがある。好ましくは、下記式（８）のように、緩和係数βを予め決めておき、少しずつスリップ速度を目標値まで上げることが望ましい。なお、ｎは繰り返しの計算回数である。
【数７】
0005514236.tif 000008
【００６７】
そして、計算されたスリップ速度ｖslipを与えた速度場と圧力場とをSIMPLEアルゴリズム等でカップリングし、圧力場を修正するための方程式が構築される（上記ステップＳ５２）。
【実施例】
【００６８】
本発明の効果を確認するために、図３のチャンバーモデルと表１の条件とに基づいて、可塑性材料のシミュレーションが行われた。評価は、計算時間及び計算精度とした。また、シミュレーションのタイムステップは１．９７３×10-3秒、可塑性材料の充填率は７０％とした。また、計算精度の評価は次の通りとし、○：良好、△：やや不良、×：不良の３段階で評価した。
分配：チャンバーモデル内でどれだけ均一に練られているかを示す指標
分散：粉体を入れて混ぜるとき、その粉体がいかに小さくなっているかを示す指標
ロータートルク：任意の回転速度でロータを回転させるのに必要なトルク
可塑性材料の温度：シミュレーションの場合は、着目している相の温度の平均値
テストの結果は、表１に示される。
【００６９】
【表１】
0005514236.tif 000009
【００７０】
テストの結果、実施例の方法では、計算時間を短縮しつつ計算精度を維持していることが確認できた。また、図８には、図６の状態から混練シミュレーションを開始した後の状態を視覚化した断面図が示される。黒い部分が材料モデルである。このように、本実施形態によれば、可塑性材料の混練状態を視覚化して再現することが出来る。
【符号の説明】
【００７１】
１  バンバリーミキサー
２  ケーシング
３  ロータ
４  チャンバー
５  チャンバーモデル
【図１】
0005514236.tif 000010
【図２】
0005514236.tif 000011
【図５】
0005514236.tif 000012
【図６】
0005514236.tif 000013
【図７】
0005514236.tif 000014
【図８】
0005514236.tif 000015
【図３】
0005514236.tif 000016
【図４】
0005514236.tif 000017