JP5564074B2 - 流体のシミュレーション方法 - Google Patents

Description

本発明は、粘度の高い流体を対象とした流体のシミュレーション方法に関する。

近年、コンピュータを用いた数値シミュレーションが種々提案されており、例えば、架橋前のゴムや樹脂といった可塑性を有する流体が、バンバリーミキサー等の空間内で混練されるときの流動状態などをナビエストークス方程式を用いて計算するものが提案されている。

特開２０１１−２７５９３号公報 V. Collin, E. Peuvrel-Disdier 、「Numerical and Experimental Study of Dispersive Mixing of Agglomerates」

ところで、空気のように粘度が低い流体の流れを計算する場合、該空気が流れる空間の壁面での空気の速度は０とすることができる。しかし、架橋前のゴムのように粘度の高い流体は、壁面に対してスリップし、壁面でも流れ方向に速度を持つ。従って、このような高粘度の流体についてのシミュレーションでは、壁面での流体のスリップ現象を無視することができない。

従来、壁面での流体の速度（以下、このような速度を「スリップ速度」という。）は、壁面でのせん断応力等の関数で定義されている。Navie's Lawを用いたせん断応力τ_wは、線形のスリップ速度の場合、次式で表される。
τ_w＝Ｆ_slip［ｖ_slip−ｖ_wall］ …（１）

同様に、非線形のスリップ速度の場合、せん断応力τ_wは、次式（２）で表される。
τ_w＝Ｆ_slip［ｖ_slip−ｖ_wall］｜ｖ_slip−ｖ_wall｜^eslip-1 …（２）

ここで、"Ｆ_slip"及び"eslip"はユーザによって決定される定数、"ｖ_slip"は流体の壁面での該壁面と平行な速度、"ｖ_wall"は壁面の速度の該壁面と平行な速度成分である。

定数Ｆ_slipは、スリップのしやすさを示す流体に固有の値であり、通常、実験により決定される。実験は、例えば、上記特許文献１に示される装置を用いて行われる。このような実験により、流体の流れる空間の壁面でのせん断応力τ_w、壁面でのスリップ速度ｖ_slip 、及び、壁面の移動速度ｖ_wallを計測し、これらの関係からＦ_slipが予め同定される。より、具体的には、ｘ軸にスリップ速度ｖ_slip、ｙ軸にせん断応力τ_wをとった両対数のグラフから累乗近似を用いて曲線ｙ＝ａ・ｘ^ｂを求め、このときの係数ａをＦslipとして、また累乗ｂをeslipとしてそれぞれ同定することができる。そして、下記式（３）で壁面におけるスリップ速度(ｖ_slip)を決め、ソルバー側に与えられる。ソルバーは、スリップ速度を用いて計算の収束解を求める。
ｖ_slip＝ｖ_wall＋τ_w／Ｆ_slip …（３）

しかしながら、上記式（３）を用いた場合、定数Ｆ_slipの値が小さいとき又は壁面せん断応力τ_wが異常な値を示しているときには、ソルバー側に与えられる流体のスリップ速度ｖ_slipも異常値を示す。このため、従来の流体シミュレーションにおいて、スリップ速度を考慮すると、計算が収束しないか又は計算が発散することがあり、計算が不安定になるという問題があった。

本発明は、以上のような問題点に鑑み案出なされたもので、粘度の高い流体の流動計算を安定させることができる流体のシミュレーション方法を提供することを主たる目的としている。

本発明のうち請求項１記載の発明は、壁面を有するチャンバー内の流体の状態を、コンピュータを用いて計算する流体のシミュレーション方法であって、前記チャンバーを有限個の要素でモデル化したチャンバーモデルを設定するステップ、前記流体をモデル化した材料モデルを設定するステップ、及び前記チャンバーモデル内に、前記材料モデルを配置しかつ予め定めた条件に基づいて流動計算を行うステップを含み、前記流動計算では、前記材料モデルと前記チャンバーモデルの前記壁面とが接触する接触面において、前記材料モデルの前記接触面と平行な速度である線形のスリップ速度が定義され、前記スリップ速度ｖ_slipを下式（ａ）で計算することを特徴としている。
ｖ_slip＝αｖ_t＋（１−α）・ｖ_wall …（ａ）
ただし、"ｖ_t"は前記接触面から法線方向に距離ｄ_wallを隔てた位置での材料モデルの接触面と平行な速度成分、"ｖ_wall"は前記接触面の速度の前記接触面と平行な速度成分、"α"はスリップ率で下式（ｂ）を満たす０以上１以下の変数である。
α／（１−α）＝μ／（ｄ_wall・Ｆ_slip） …（ｂ）
ただし、"μ"は材料モデルの粘度、"Ｆ_slip"は定数である。

また請求項２記載の発明は、壁面を有するチャンバー内の流体の状態を、コンピュータを用いて計算する流体のシミュレーション方法であって、前記チャンバーを有限個の要素でモデル化したチャンバーモデルを設定するステップ、前記流体をモデル化した材料モデルを設定するステップ、及び前記チャンバーモデル内に、前記材料モデルを配置しかつ予め定めた条件に基づいて流動計算を行うステップを含み、前記流動計算では、前記材料モデルと前記チャンバーモデルの前記壁面とが接触する接触面において、前記材料モデルの前記接触面と平行な速度である非線形のスリップ速度が定義され、前記スリップ速度ｖ_slipを下式（ａ）で計算することを特徴としている。
ｖ_slip＝αｖ_t＋（１−α）・ｖ_wall …（ａ）
ただし、"ｖ_t"は前記接触面から法線方向に距離ｄ_wallを隔てた位置での材料モデルの接触面と平行な速度成分、"ｖ_wall"は前記接触面の速度の前記接触面と平行な速度成分、"α"はスリップ率で下式（ｃ）を満たす０以上１以下の変数である。
α／（１−α）＝μ／｛（ｄ_wall・Ｆ_slip）｜ｖ_slip−ｖ_wall｜^eslip-1｝ …（ｃ）
ただし、"μ"は材料モデルの粘度、"Ｆ_slip"及び"eslip"は定数である。

また請求項３記載の発明は、前記チャンバーは、バンバリーミキサーのケーシングと、該ケーシング内を回転する少なくとも１本のロータとの間で形成されるバンバリーミキサーの混練空間であり、前記流体は、前記ロータによって混練される架橋前のゴム又は樹脂材料である請求項１又は２記載の流体のシミュレーション方法である。

本発明のシミュレーション方法は、流動計算において、材料モデルは、チャンバーモデルの壁面との接触面において、スリップを考慮した計算が可能になる。従って、本発明では、粘度の高い架橋前のゴムや樹脂といった可塑性を有する流体について精度の良いシミュレーション方法を提供できる。

また、材料モデルのスリップ速度ｖ_slipは、上記式（ａ）で計算される。

式（ａ）において、スリップ率αは、０以上１以下の値として定められる。従って、式（ａ）においてスリップ率α＝０の場合、スリップ速度ｖ_slipは、接触面の速度ｖ_wallとなる。接触面の速度ｖ_wallが０の場合、材料モデルの速度は０である。これは、接触面において、材料モデルにスリップが発生しない状態を示す。

一方、式（ａ）において、スリップ率α＝１の場合、スリップ速度ｖ_slipは、接触面から法線方向に距離ｄ_wallを隔てた位置での材料モデルの接触面と平行な速度ｖ_t となる。これは、接触面において、材料モデルがフルスリップ（摩擦なし）の状態に相当している。

以上のように、本発明では、スリップ速度ｖ_slipを安定して計算することができ、ひいては、シミュレーションでの流動計算を収束させることができる。

可塑性材料を混練するバンバリーミキサーの部分断面図である。 本実施形態の処理手順の一例を示すフローチャートである。 チャンバーモデルの斜視図である。 （ａ）はチャンバーモデルの断面図、（ｂ）はそのＸ部拡大図である。 チャンバーモデルを分解して示す断面図である。 チャンバーモデル内に材料モデルと気相モデルとを混在して配置した状態を示す断面図である。 流動計算の処理手順の一例を示すフローチャートである。 材料モデルの速度を説明する接触面近傍の断面図である。 スリップ率と式（ｂ）又は（ｃ）の右辺の値との関係を示すグラフである。 混練シミュレーションをした後の状態を視覚化した線図である。

以下、本発明の実施の一形態が図面に基づき説明される。
本発明は、コンピュータ（図示省略）を使用して、チャンバー内を流れる流体の流動状態をシミュレートするシミュレーション方法である。

前記流体としては、チャンバーの壁面との接触面に対して大きなスリップが生じるような比較的粘度の高い流体が対象とされる。該流体としては、例えば、架橋前のゴム、樹脂又はエラストマー等の可塑性材料が相当する。これらの可塑性材料は、十分に練られて安定的な流動状態（流体）とみなすことができる状態のものが前提となる。このような状態としては、架橋前のゴムの場合、十分に練られて約８０℃程度まで昇温した状態が相当する。

前記チャンバーとは、上述の可塑性材料が配されかつ移動する壁面で囲まれた空間を意味する。チャンバーには、種々の形状のものが含まれる。代表的なチャンバーとしては、バンバリーミキサーが有する混練空間が挙げられる。

図１には、バンバリーミキサー１の断面図が示されている。該バンバリーミキサー１は、ケーシング２と、該ケーシング２内を回転する一対のロータ３、３と、これらの間で形成される断面横向きの略８の字状のチャンバー４（混練空間）とを有する。本実施形態では、このチャンバー４が解析対象とされる。ただし、チャンバー４は、このような形状に限定されるものではない。

図２には、本実施形態の処理手順の一例が示されており、以下、各ステップが順に説明される。

［チャンバーモデルの設定］
本実施形態のシミュレーション方法では、先ず、チャンバー４を有限個の要素でモデル化したチャンバーモデルがコンピュータに入力される（ステップＳ１）。

図３は、チャンバーモデル５を視覚化した斜視図、図４（ａ）はその断面図、図４（ｂ）は図４（ａ）のＸ部拡大図がそれぞれ示されている。チャンバーモデル５は、前記ケーシング２の内周面をなす外周面５ｏと、前記回転する一対のロータ３、３の外周面である内周面５ｉと、ロータ３の軸方向の両端側で前記外周面を閉じる両端面５ｓとで閉じられた三次元空間を有し、この空間が三次元の要素ｅによって分割されている。そして、この空間を定義する面がチャンバーモデル５の壁面である。

チャンバーモデル５は、オイラー要素ｅで分割（離散化）されている。要素分割は、四面体、六面体などの他、多面体セル（ポリヘドラルグリッド）といった三次元要素で行われる。そして、各要素について、可塑性材料（材料モデル）の圧力、温度及び速度等の物理量が計算される。

また、チャンバーモデル５の外周面５ｏ及び両端面５ｓは変形しない。チャンバーモデル５の内周面５ｉは、ロータ３の回転に対応して回転する。このため、チャンバーモデル５の容積形状は変化する。本実施形態では、このようなチャンバーモデル５の形状変化をシミュレーションに取り込むために、図５に分解して示されるように、チャンバーモデル５が、一対の回転部５Ａ、５Ｂと、これらの間を継ぐ継ぎ部５Ｃと、これらが収容される外枠部５Ｄとの４つの部分に分けて構成される。

回転部５Ａ、５Ｂは、各々、円形の外周面５Ａｏ、５Ｂｏと、ロータ３の外周面に等しい前記内周面５ｉとを有する筒状でモデル化されている。そして、回転部５Ａ、５Ｂは、各々、その中心Ｏａ、Ｏｂの周りで回転する。これにより、ロータ３の回転に伴うチャンバー３の容積形状の変化が表現される。

継ぎ部５Ｃは、一対の回転部５Ａ、５Ｂ間に配されて静止している。また、継ぎ部５Ｃは、各回転部５Ａ、５Ｂに沿いかつ接触する凹円弧面ｊを有している。継ぎ部５Ｃの凹円弧面ｊと回転部５Ａ、５Ｂの外周面５Ａｏ、５Ｂｏとは、スライディングサーフェース等の境界条件が定義される。これにより、チャンバーモデル５の回転部５Ａ、５Ｂ内で生じる物理的な作用（力及び熱等）は、この凹円弧面ｊを介して継ぎ部５Ｃの中にある材料モデルに伝達される。

外枠部５Ｄは、回転部５Ａ、５Ｂ及び継ぎ部５Ｃを囲む筒状をなす。外枠部５Ｄの軸方向両端は、前記両端面５ｓによって閉じられている。また、外枠部５Ｄも、回転部５Ａ、５Ｂ及び継ぎ部５Ｃの要素との接触面に、スライディングサーフェース等の境界条件が定義される。これにより、外枠部５Ｄ内で回転する回転部５Ａ、５Ｂ内で生じる物理的な作用（力及び熱等）が、両者の接触面を介して外枠部５Ｄに伝えられる。なお、外枠部５Ｄには、ロータの作用によって大きなせん断が生じやすい。材料の速度等をより詳細に計算するために、外枠部５Ｄは、回転部５Ａ、５Ｂ及び継ぎ部５Ｃよりも小さい要素で構成されるのが望ましい。これによって、チャンバーモデル５の外周面付近の材料モデルの速度プロファイルなどがより詳細に計算され得る。

［材料モデルの設定］
次に、本実施形態では、材料モデルが設定され、コンピュータに入力される（ステップＳ２）。材料モデルは、チャンバー４内を流動する流体（可塑性材料）に相当する。材料モデルは、解析対象となる可塑性材料に応じてせん断粘度、比熱、熱伝導率及び比重等の物性値が定義される。

材料モデルのせん断粘度は、例えば、解析対象となる流体から粘弾性特性（Ｇ'及びＧ”）が複数の温度条件で測定され、Cox-Merz則などを用いてせん断粘度に変換することで得られる。このようにして得られたせん断粘度は、例えば下記式（４）のべき乗法則で近似することができる。
μ＝ｍγ'^n-1 …（４）
ここで、μ：せん断粘度、ｍ：係数、γ'：せん断速度、ｎ：係数である。

材料モデルの比熱は、解析対象の可塑性材料から、例えば断熱型連続法（＠２５℃）にて測定され、その値がコンピュータに入力される。

材料モデルの熱伝導率は、解析対象の可塑性材料から、例えば熱線法（＠２５℃）にて測定され、その値が前記コンピュータに入力される。

［気相モデルの設定］
次に、本実施形態では、気相モデルが設定され、コンピュータに入力される（ステップＳ３）。気相モデルは、チャンバー内に存在する第２の流体としての空気をモデル化したものである。このような気相モデルをシミュレーションに取り入れることにより、チャンバーへの可塑性材料の充填率が１００％未満の状態の流動計算が可能になる。

本実施形態において、気相モデルには、粘度及び比重が定義される。気相モデルの比重は、空気の特定の温度での比重に基づいて決定される。気相モデルの粘度は、実際の空気の粘度と同じ値が定義されてもよいし、異なっていても良い。好ましくは、気相モデルの比重には、実際の値が定義される。

粘度が異なる気相モデル（空気）と材料モデル（可塑性材料）との混相流解析では、両者の粘度差が非常に大きくなる。このため、後述する流動計算において、気相モデルと材料モデルとの界面では、両者の速度差が著しく大きくなる。この結果、速度差に伴うせん断発熱を考慮すると、気相モデルと材料モデルとの界面で、せん断発熱が異常に大きくなり、計算が不安定化しやすい。一方、可塑性材料の充填率が１００％の状態は現実的ではない。このため、充填率１００％未満のシミュレーションを行うためには、チャンバーモデル５内に気相モデルを含ませることが必要になる。この点に鑑み、本実施形態では、計算結果に悪影響を及ぼさない範囲で、気相モデルの粘度を極力大きくする。これにより、計算が安定する。

発明者らは、気相モデルの粘度のみを変えて複数種類の流動計算（シミュレーション）を行い、チャンバーモデル５内の圧力場を比較した。その結果、気相モデルの粘度が実際の空気の粘度の１０倍を超えると、チャンバーモデル５内の圧力値が大きくなり、圧力場に悪影響が出ることが判明した。他方、気相モデルの粘度が、実際の空気の粘度の５倍未満であると、計算を安定させるのが困難であることも分かった。従って、「計算に悪影響を及ぼさない範囲」として、気相モデルの粘度を、実際の空気の粘度の５〜１０倍の範囲で定めるのが望ましい。

［境界条件等の設定］
次に、本実施形態では、流動計算のシミュレーションに必要な境界条件等、各種の条件が定義される（ステップＳ４）。設定される境界条件としては、チャンバーモデル５の壁面での流速境界条件及び温度境界条件が挙げられる。

前記流速境界条件としては、材料モデルとチャンバーモデル５の壁面とが接触する接触面には、スリップが定義される。これにより、流動計算において、材料モデルは、前記接触面で速度を持つ。これは、スリップ速度と呼ばれる。これについては、後述する。

前記温度境界条件としては、次の２種類が設定可能である。
ａ）断熱条件：温度の計算を単純化するため、チャンバーモデル５の表面において、熱が外に逃げない条件。
ｂ）全てのチャンバーモデル５の表面温度が温調温度（例えば５０℃）に設定される条件。

他の条件としては、材料モデルの初期温度がある。本実施形態の初期温度は、２０℃（２９３Ｋ）である。また、バンバリーミキサーのロータの回転数に対応するものとして、チャンバーモデルの回転部５Ａ、５Ｂの回転数が定義される。さらに、チャンバーモデル５の容積に対する材料モデルの充填率が定義される。

さらに、他の条件としては、流動計算（シミュレーション）の初期状態、タイムステップ、内部処理でのイタレーションの反復回数、計算終了時刻などがある。初期状態は、例えば、図６に示されるように、チャンバーモデル５を横切る水平な境界面Ｓを基準として、それよりも上部を気相モデルの領域Ａとし、それよりも下部を材料モデルの領域Ｍとして定義される。境界面Ｓのレベルを変えることにより、材料モデル（可塑性材料）の充填率が調節される。これらの条件は、シミュレーションの目的等に応じて任意に定められる。

［流動計算］
本実施形態では、図６に示したように、チャンバーモデル５内に、材料モデル（領域Ｍ）と気相モデル（領域Ａ）とを混在して配置しかつ予め定められた上記各条件に基づいて流動計算が行われる（ステップＳ５）。流動計算では、材料モデルの運動状態を特定する３方向（ｘ，ｙ，ｚ）の速度成分、材料モデルの内部状態を特定する未知量である圧力ｐ及び温度Ｔが計算される。つまり、解くべき未知数はこの５つである。なお、本実施形態では、非圧縮性流れの場合のナビエストークス方程式とし、気相モデル及び材料モデルの各密度は、一定として取り扱われる。

本実施形態の流動計算において、材料モデルは、全温度領域で流体として扱われる。このため、流体の方程式（ナビエストークス方程式、質量保存式、エネルギー方程式の連立）を解くことになる。

本実施形態の流動計算では、空気と可塑性材料との２つの流体が一度に扱われるため、自由界面の流れの計算で用いられるＶＯＦ（Volume of Fluid）法が用いられる。ＶＯＦ法では、二つの流体の界面の移動を直接計算するのではなく、チャンバーモデル５の各要素（「セル」ということもある。）の体積中の材料モデルの充填率（体積分率）を定義して自由界面が表現される。なお、支配方程式は次の通りである。

［運動方程式］
本実施形態では、チャンバーモデル５の中に、気相モデル及び材料モデルの２相が配置される混相流モデルが一流体モデルとして扱われる。しかし、運動方程式は、下記式（５）の１組をそれぞれｘ，ｙ，ｚの３方向について解くだけで良い。ＶＯＦ法により、気相と液相とが平均化され、一つの流体として扱われるためである。
上記式において、符号は次の通りである。
ｕ：混相流モデルの速度
ｐ：混相流モデルの圧力
ρ：混相流モデルの密度
ｇ：重力加速度
Ｔ：混相流モデルの絶対温度
Ｆ：外力

また、材料モデルと気相モデルとを含む要素の密度ρ及び粘性係数は、下記式（６）のように、材料モデルの相及び気相モデルの相の各相がそれぞれ占める体積で重み付けされて平均化されたものが使用される。
上記式において、符号は次の通りである。
Δ_ｑ：体積分率
ρ_ｑ：各セルにおける各相の密度
μ_ｑ：各セルにおける各相の粘度

［質量保存式］
同様に、質量保存（連続の式）及び圧力方程式も一組のみ解くことで足りる。つまり、本実施形態のシミュレーション方法では、混相でありながら、流れ場の計算としては単相と同じであり、場所（体積分率）によって、物性値が異なる流れを解いていることになる。各相の位置は、計算結果として得られる体積分率の分布（Volume Fluxion）より判断できる。

［エネルギー方程式］
エネルギー方程式については、下記式（７）で表現される。また、材料モデルの温度はこの式によって求めることができる。
上記式において、符号は次の通りである。
Ｅ：エンタルピ
ｋ：熱伝導率
Ｓ：ソース項

［体積分率輸送方程式］
体積分率の分布は、気相モデルと材料モデルとの２相間の界面位置を決定するものである。この体積分率δ_ｑは、下記式（８）を精度良く解くことで得られる。

ここで、チャンバーモデル５の任意の要素ｅにおいて、体積分率δ_ｑ＝０の場合、当該要素にはｑ相（ｑ相：材料モデルの相）は存在しないことを示す。また、体積分率δ_ｑ＝１の場合、当該要素はｑ相（材料モデルの相）で全部が満たされていることを示す。さらに、０＜δ_ｑ＜１の場合、当該要素は、材料モデルと気相モデルとで満たされていること、つまり混相であることを示す。この方程式は、例えば、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 …（９）

次に、上記ステップで得られた圧力補正方程式が、反復法で計算され、圧力補正量ｐ’が計算される（ステップＳ５５）。反復法には、例えば、ＡＭＧソルバー、ＣＧ又はＢｉ−ＣＧ等が採用される。

次に、計算された解を元に、圧力場が修正される（ステップＳ５６）。圧力場の補正は、次式（１０）が用いられる。
ｐ^ｎ＋１＝ｐ^ｎ＋ωｐ’ …（１０）
ここで、ｐは圧力、ｎは現在のステップ数、ωは緩和係数であり、本実施形態では０．３が採用されている。緩和係数には、他の値が採用されても良い。

次に、境界面の境界条件が修正（アップデート）される。即ち、修正された圧力場から圧力勾配を求め、その圧力勾配を境界条件として与える（ステップＳ５７）。上記「境界面」とは、材料モデル及び気相モデルを含む流体モデルと、それ以外の固体（ロータ及びチャンバー）を表すモデルとの境界面である。

次に、セル表面の質量流量が修正される（ステップＳ５８）。修正は次式（１１）が用いられる。
ｍ_ｆ ^ｎ＋１＝ｍ_ｆ ^＊＋ｍ'_ｆ …（１１）
ここで、ｍ_ｆ ^ｎ＋１は修正後のセル表面の質量流量、ｍ_ｆ ^＊は修正前のセル表面の質量流量、ｍ'_ｆは質量流量の修正値である。

次に、速度場が修正される（ステップＳ５９）。速度場の修正は、例えば、次式（１２）が用いられる。
ｖ^ｎ＋１＝ｖ^＊−（Ｖ▽ｐ'／∂_ｐＶ） …（１２）
ここで、Ｖはセルの体積、ｖ^＊は中間的な速度場（運動方程式から求めた補正前のもの）、∂_ｐＶは運動方程式のマトリクスの対角成分、▽ｐ'は圧力補正量の勾配である。

次に、式（７）のエネルギー方程式を解くことにより、材料モデルの温度及び粘度が算出される（ステップＳ６０）。

さらに、チャンバーモデル５の壁面スリップ速度ｖ_slipが下記式（ａ）で計算される（ステップＳ６１）。

図８に示されるように、式（ａ）の符号ｖ_ｔは、チャンバーモデルの接触面ＣＰから法線方向に距離ｄ_wallを隔てた位置での材料モデルの接触面ＣＰと平行な速度成分である。また、符号ｖ_wallは接触面ＣＰの速度の接触面ＣＰと平行な速度成分、符号αはスリップ率で０以上１以下の変数である。さらに、スリップ率αは、図８のような線形のスリップ速度が与えられる場合には下式（ｂ）を満たし、非線形のスリップ速度（図示省略）が与えられる場合には下式（ｃ）をそれぞれ満たす。

式（ｂ）及び（ｃ）において、符号μは材料モデルの粘度、符号Ｆ_slip及びeslipは定数である。

線形のスリップ速度の場合、接触面ＣＰでのせん断応力τwは、「Computational Methods For Fluid Dynamics」, p256-2567 式(8.74), (8.75)によれば、一般的に、下式（１３）で表すことができる。

計算が収束し、正確な流れ場が得られている状態では、式（１）及び式（１３）により、次式（１４）が成り立つ。

式（１４）に、式（ａ）を代入して変形すると、式（１５）が得られる。

式（１５）をさらに整理すると、式（ｂ）が得られる。式（１５）において、符号μは材料モデルの粘度であり、固有の値である。また、式（ｂ）の符号ｄ_wallは、例えば、接触面ＣＰに接する要素の中心点（代表点）が利用でき、これはチャンバーモデル５の要素分割の方法で決まる既知の値である。さらに、式（ｂ）の符号Ｆ_slipは、予めユーザによって決定される定数であり、これも既知の値である。従って、スリップ率αは、式（ｂ）からを簡単に同定することができる。そして、スリップ速度ｖ_slipは、式（ｂ）で求められたスリップ率αを式（ａ）に代入することによって、容易に求めることができる。

このような式を用いることにより、材料モデルのスリップ速度ｖ_slipが異常値を示すことがない。従って、流体シミュレーションのロバスト性が向上し、速度のシュートアップによる精度悪化も防ぐことができる。また、計算の収束性も向上させることができ、計算時間の短縮も図ることができる。

なお、非線形のスリップ速度については、式（２）及び式（１３）より、次式（１６）が成り立つ。そして、この式（１６）に、式（ａ）を代入して整理すると、上記式（ｃ）が得られる。

上記式（ｃ）では、スリップ速度ｖ_slipと接触面の速度ｖ_wallとの差が小さいとき、非常に小さい値が分母に入る。この場合、式（ｃ）の計算中、いわゆるゼロ割りが生じるおそれがある。そこで、非線形のスリップ速度が定義される場合でも、計算を安定させるために、式（ｃ）の速度差ｖ_slip−ｖ_wall の値が、予め定められた閾値（例えば、１．０×１０^−８程度）よりも小さくなった場合には、式（ｂ）を適用することが望ましい。これにより、非線形のスリップ速度を基調とした安定した計算が行える。

図９は、式（ｂ）、式（ｃ）の右辺の値と、スリップ率αとの関係を示すグラフである。該グラフの縦軸はスリップ率αの値、横軸は式（ｂ）又は（ｃ）の右辺の値を示している。式（ｂ）及び式（ｃ）において、右辺の値は正の値をとる。右辺の値を正の数として、非常に小さな値（１．０×１０^−６）から非常に大きな値（１．０×１０^６）まで変化させても、スリップ率αの値は、０〜１の範囲で変化するに過ぎない。従って、本発明によれば、式（ｂ）を用いる限り、スリップ率αが異常値をとることが無い。

次に、コンピュータは、計算の解が収束したか否を判断する（ステップＳ６２）。この収束したか否かは、例えば、修正された質量流量の総量が、予め定めた誤差の範囲内に留まったか否かを調べることにより行いうる。「計算の解」とは、本流体シミュレーションの支配方程式であるナビエ・ストークス方程式（式（５）＋連続の式）の解である。

そして、収束したと判断された場合（ステップＳ６２でＹ）、１タイムステップ分、次の時刻に移動する（ステップＳ６３）。次に、その時刻が指定された時刻に至っていない場合（ステップＳ６４でＮ）、再び、ステップＳ５１以降が繰り返される。他方、ステップＳ６２で収束していないと判断された場合には、ステップＳ５１以降が繰り返される。

なお、本実施形態の流動計算では、せん断発熱については、材料モデルについてのみ計算が行われ、気相モデルのせん断発熱の計算は行われない。材料モデルについては、可塑性材料であるため、せん断発熱を考慮する必要があるが、材料モデルと気相モデルとの界面では、せん断速度が大きくなるケースが多い。ここで、せん断発熱は、「せん断粘度」に「せん断速度」の二乗を乗じたものであるから、せん断速度の異常に大きな値は、温度の計算（エネルギー方程式）に悪影響を与えるおそれがある。そこで、本実施形態では、せん断発熱の計算を材料モデルのみに限定、より具体的には、チャンバーモデル５の要素において、材料モデルの体積分率δqの値が一定値以上（本実施形態では０．９０以上）の要素に限定することで計算の安定化が図られている。

以上本発明の実施形態について詳述したが、本発明は、上記の実施形態に限定されるものではなく、種々の態様に変更して実施することができる。

本発明の効果を確認するために、図３のチャンバーモデルを用いて、未加硫ゴムからなる材料モデルを、実時間で２０秒間混練する流体シミュレーションが行われた。流体シミュレーションの主な条件等は次の通りである。
定数Ｆ_slip：
ロータに対応している回転部の回転数：３０rpm
タイムステップ： １．９７３×１０^-3秒
ゴム材料のチャンバーモデル内への充填率：７０％

本発明の式（ａ）及び（ｂ）を用いてスリップ速度を計算した実施例では、スリップ速度に異常値は一切見られず、最後まで安定して計算が行われた。図１０には、このような混練シミュレーションが視覚化されて示されている。図中、黒色部が、材料モデルである。一方、従来の技術で説明した方法によるスリップ速度の計算方法では、スリップ速度に異常値が頻出し、計算が収束しなかった。

１ バンバリーミキサー
２ ケーシング
３ ロータ
４ チャンバー
５ チャンバーモデル

Claims (3)

1. 壁面を有するチャンバー内の流体の状態を、コンピュータを用いて計算する流体のシミュレーション方法であって、
   前記チャンバーを有限個の要素でモデル化したチャンバーモデルを設定するステップ、
   前記流体をモデル化した材料モデルを設定するステップ、及び
   前記チャンバーモデル内に、前記材料モデルを配置しかつ予め定めた条件に基づいて流動計算を行うステップを含み、
   前記流動計算では、前記材料モデルと前記チャンバーモデルの前記壁面とが接触する接触面において、前記材料モデルの前記接触面と平行な速度である線形のスリップ速度が定義され、
   前記スリップ速度ｖ_slipを下式（ａ）で計算することを特徴とする流体のシミュレーション方法。
   ｖ_slip＝αｖ_t＋（１−α）・ｖ_wall …（ａ）
   ただし、"ｖ_t"は前記接触面から法線方向に距離ｄ_wallを隔てた位置での材料モデルの接触面と平行な速度成分、"ｖ_wall"は前記接触面の速度の前記接触面と平行な速度成分、"α"はスリップ率で下式（ｂ）を満たす０以上１以下の変数である。
   α／（１−α）＝μ／（ｄ_wall・Ｆ_slip） …（ｂ）
   ただし、"μ"は材料モデルの粘度、"Ｆ_slip"は定数である。
   
2. 壁面を有するチャンバー内の流体の状態を、コンピュータを用いて計算する流体のシミュレーション方法であって、
   前記チャンバーを有限個の要素でモデル化したチャンバーモデルを設定するステップ、
   前記流体をモデル化した材料モデルを設定するステップ、及び
   前記チャンバーモデル内に、前記材料モデルを配置しかつ予め定めた条件に基づいて流動計算を行うステップを含み、
   前記流動計算では、前記材料モデルと前記チャンバーモデルの前記壁面とが接触する接触面において、前記材料モデルの前記接触面と平行な速度である非線形のスリップ速度が定義され、
   前記スリップ速度ｖ_slipを下式（ａ）で計算することを特徴とする流体のシミュレーション方法。
   ｖ_slip＝αｖ_t＋（１−α）・ｖ_wall …（ａ）
   ただし、"ｖ_t"は前記接触面から法線方向に距離ｄ_wallを隔てた位置での材料モデルの接触面と平行な速度成分、"ｖ_wall"は前記接触面の速度の前記接触面と平行な速度成分、"α"はスリップ率で下式（ｃ）を満たす０以上１以下の変数である。
   α／（１−α）＝μ／｛（ｄ_wall・Ｆ_slip）｜ｖ_slip−ｖ_wall｜^eslip-1｝ …（ｃ）
   ただし、"μ"は材料モデルの粘度、"Ｆ_slip"及び"eslip"は定数である。
   
3. 前記チャンバーは、バンバリーミキサーのケーシングと、該ケーシング内を回転する少なくとも１本のロータとの間で形成されるバンバリーミキサーの混練空間であり、
   前記流体は、前記ロータによって混練される架橋前のゴム又は樹脂材料である請求項１又は２記載の流体のシミュレーション方法。