JP6014457B2 - 粘弾性流体のシミュレーション方法、粘弾性流体のシミュレーション装置及び粘弾性流体のシミュレーションプログラム - Google Patents

Description

本発明は、未加硫ゴムなどの粘弾性流体の流路における流動をシミュレートする粘弾性流体のシミュレーション方法、シミュレーション装置及びシミュレーションプログラムに関する。

従来から、例えば特許文献１に例示されるように、押出機によるゴムの押出成形などの分野においては、安定した成形を可能にするといった目的で、押出ダイなどの金型の設計および成形条件などの最適な設定を支援する各種解析システムが利用されている。そのような解析の一つの手段として、押出機などの流路内における流体のシミュレーションがあり、ゴムなどの成形においては、粘弾性流体の圧力、速度及び応力を含む物理量を算出するシミュレーションが挙げられる。

粘弾性流体に限らず、ニュートン流体などの流体の流れの数値解析には、有限差分法が用いられるのが通例である。有限差分法では、図３に例示するように、流路等の解析領域をメッシュ状に分割した、格子配列の位置が規則正しく並んでいる構造格子と呼ばれる計算格子を用いる。流体の挙動を支配する支配方程式であるナビエ・ストークス方程式（偏微分方程式）を有限差分法を用いて離散化することで、ナビエ・ストークス方程式をコンピュータで取り扱うことが可能な代数方程式に近似でき、近似した代数方程式（差分方程式）を全ての格子点について数値計算を行うことでコンピュータを用いて近似解を求解する。一例として非特許文献１等が挙げられる。この数値計算法は、図４に例示すように、算出対象となる格子点を評価点ＡＰ_{（ｉ，ｊ）}とする場合に、評価点ＡＰ_{（ｉ，ｊ）}の微分値（ｘ方向）を、評価点ＡＰ_{（ｉ，ｊ）}に隣接する２点の格子点GP1，GP2の値に基づき算出する。ｘ方向とｙ方向の微分値をそれぞれ算出するために、周囲４つの格子点GP1〜GP4を用いる。

ところが、ワイゼンベルグ数が高いゴムなどの粘弾性流体の縮流管内の流れにおいては、応力の勾配が大きい領域（急勾配となる領域）が存在する。流体を支配する応力テンソルの式を上記非特許文献１等のアルゴリズムを用いて数値計算するにあたり、応力の勾配が急勾配であれば、方程式を満たさずに発散してしまうことが知られている。ここでは、ワイゼンベルグ数が高い場合について説明したが、流路の形状により速度勾配や圧力勾配が急勾配の場合も同様に発散してしまう。

このような場合であっても粘弾性流体の数値解析を可能にするべく、応力の勾配を緩やかにするような数値計算上の工夫を施した技術が非特許文献２に開示されている。この数値計算法は、図５Ｇに示すように、評価点ＡＰ_{（ｉ，ｊ）}の微分値を、評価点ＡＰ_{（ｉ，ｊ）}の周囲にある上記２点よりも多い数の格子点GP1〜GP20の値に基づき算出するアルゴリズムを用いる。すなわち、評価点の微分値を計算するために用いる格子点の数を広範囲にするものである。

特許第４８００７７６号公報

"M.F.Tome ,A.Castelo ,V.G.Ferreira ,S.McKee"、"A finite difference technique for solving the Oldroyd-B model for 3D-unsteady free surface flows."、Department of Mathematics, University of Strathclyde, Glasgow, UK, Journal of Non-Newtonian Fluid Mechanics. 「壁瀬航平，福井智宏，森西晃嗣」、「弾性体まわりの流れに対する仮想流束法の検討」、日本機械学会 関西支部 第86期定時総会講演会、No.114-1, pp.13-1、京都、2011年3月

しかしながら、上記非特許文献２のアルゴリズムを用いる場合には、数値計算が発散状態にならず的確に収束に向かうという利点があるものの、上記非特許文献１のアルゴリズムを用いる場合に比べて計算量が格段に増大してしまう。

一方、上記非特許文献１のアルゴリズムを用いる場合には、非特許文献２のアルゴリズムを用いる場合に比べて計算量が格別に少ないという利点があるものの、数値計算が発散してしまうおそれがある。

本発明は、このような課題に着目してなされたものであって、その目的は、高い収束性と低計算コストとを両立する粘弾性流体のシミュレーション方法、シミュレーション装置及びシミュレーションプログラムを提供することである。

本発明は、上記目的を達成するために、次のような手段を講じている。

すなわち、本発明の粘弾性流体のシミュレーション方法は、粘弾性流体の圧力、速度及び応力を含む物理量に関する微分方程式に基づき定常状態における前記物理量を算出するために、コンピュータが実行する粘弾性流体のシミュレーション方法であって、
前記粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する前記物理量の初期値の設定、及び前記流路の境界条件の設定を実行するステップと、
算出対象となる格子点を評価点とする場合に、前記評価点に隣接する所定数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第一のアルゴリズムを用い、全ての格子点における前記物理量を算出し、算出値に基づき所定の終了条件が成立するかを判定し、当該所定終了条件の成立まで前記物理量の算出を繰り返すステップと、
前記評価点の周囲にある前記所定数よりも多い数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第二のアルゴリズムを用い、全ての格子点における前記物理量の算出を繰り返すステップと、を含み、
前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、算出値の前回値に対する変動量を表す残差に基づき発散状態であると判定される場合、又は、前記第一のアルゴリズムによる物理量の算出を所定回数以上実行したと判定される場合の少なくともいずれかに該当するときに前記所定の終了条件が成立すると判定し、
当該終了条件が成立すると判定した場合に、前記第一のアルゴリズムによる複数の算出結果のうち前記残差が最も小さい時点の算出値を初期値とし、前記第二のアルゴリズムによる算出を開始することを特徴とする。

収束状態ともいえる定常状態では、時間変化に応じて物理量が変わらないので、算出値の前回値に対する変化量を表す残差が０又は限りなく小さくなり、逆に発散状態では、残差が徐々に大きくなるので、残差に基づき発散状態であるか定常状態に収束したかを判定できる。これを前提とし、算出対象となる評価点に隣接する所定数の格子点の値に基づいて評価点の微分値を差分法で算出する第一のアルゴリズムの実行時に、収束が見込めない発散状態であるとき、又は、計算を所定回数以上実行したときのいずれかに該当するときに繰り返し計算を終了するようにしているので、計算が終了しないことを回避しつつ、算出値を或る程度真値に向かわせることが可能となる。それでいて、第一のアルゴリズムによる複数の算出結果のうち残差が最も小さい時点の算出値を初期値とし、第二のアルゴリズムによる算出を開始するので、計算負荷の大きい第二のアルゴリズム単独で計算する場合に比べて全体の計算コストを削減できるとともに、確実に収束させることが可能となる。
したがって、第一及び第二のアルゴリズムの良い所取りが可能となり、高い収束性と低計算コストとを両立することが可能となる。

計算コストを低減するためには、前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、前記残差のうち最小残差を記録しておき、前記最小残差が更新されるときに、全ての格子点について前記最小残差となる時点の物理量を記憶することが好ましい。

複雑な計算を要することなく、発散状態であるかの判定を可能にするためには、前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、前記物理量の算出結果に基づき算出する残差と予め記録されている最小残差との比率が所定閾値を超えない場合に発散状態でないと判定し、前記比率が前記所定閾値を超える場合に発散状態であると判定することが好ましい。

高い収束性を確保したまま、計算コストをより一層削減するためには、前記第二のアルゴリズムにより物理量の算出を繰り返すステップにおいて、定常状態であると規定する残差の規定値よりも大きな値の基準値を予め設定しておき、前記残差と前記基準値とを比較し、前記残差が前記基準値よりも小さいと判定した場合に、前記第二のアルゴリズムによる繰り返し計算を停止し、その後、前記第二のアルゴリズムによる算出結果を初期値とし、再度、前記第一のアルゴリズムによる繰り返し計算を開始し、定常状態であると規定する残差の規定値より残差が小さくなることを以て、繰り返し計算を終了して定常状態における物理量を算出することが好ましい。

本発明は、上記シミュレーション方法を実行する装置として特定可能である。
すなわち、本発明の粘弾性流体のシミュレーション装置は、粘弾性流体の圧力、速度及び応力を含む物理量に関する微分方程式に基づき定常状態における前記物理量を算出するために用いられる粘弾性流体のシミュレーション装置であって、前記粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する前記物理量の初期値の設定、及び前記流路の境界条件の設定を実行する初期設定部と、算出対象となる格子点を評価点とする場合に、前記評価点に隣接する所定数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第一のアルゴリズムを用い、全ての格子点における前記物理量の算出を繰り返し、その後、前記評価点の周囲にある前記所定数よりも多い数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第二のアルゴリズムを用い、全ての格子点における前記物理量の算出を繰り返す物理量算出部と、前記物理量算出部による繰り返し計算の終了条件が成立したか否かを判定する判定部と、を備え、前記判定部は、前記物理量算出部が前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、算出値の前回値に対する変動量を表す残差に基づき発散状態であると判定される場合、又は、前記第一のアルゴリズムによる物理量の算出を所定回数以上実行したと判定される場合の少なくともいずれかに該当するときに終了条件が成立すると判定し、前記物理量算出部は、第一のアルゴリズムによる繰り返し計算の終了条件が成立すると判定された場合に、前記第一のアルゴリズムによる複数の算出結果のうち前記残差が最も小さい時点の算出値を初期値とし、前記第二のアルゴリズムによる算出を開始することを特徴とする。この装置によっても、上記シミュレーション方法が奏する作用効果を得ることができる。

本発明は、上記シミュレーション方法を構成するステップを、プログラムの観点から特定することも可能である。

すなわち、本発明の粘弾性流体のシミュレーションプログラムは、粘弾性流体の圧力、速度及び応力を含む物理量に関する微分方程式に基づき定常状態における前記物理量を算出するために、下記ステップをコンピュータに実行させる粘弾性流体のシミュレーションプログラムであって、
前記粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する前記物理量の初期値の設定、及び前記流路の境界条件の設定を実行するステップと、算出対象となる格子点を評価点とする場合に、前記評価点に隣接する所定数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第一のアルゴリズムを用い、全ての格子点における前記物理量を算出し、算出値に基づき所定の終了条件が成立するかを判定し、当該所定終了条件の成立まで前記物理量の算出を繰り返すステップと、
前記評価点の周囲にある前記所定数よりも多い数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第二のアルゴリズムを用い、全ての格子点における前記物理量の算出を繰り返すステップと、を含み、
前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、算出値の前回値に対する変動量を表す残差に基づき発散状態であると判定される場合、又は、前記第一のアルゴリズムによる物理量の算出を所定回数以上実行したと判定される場合の少なくともいずれかに該当するときに前記所定の終了条件が成立すると判定し、
当該終了条件が成立すると判定した場合に、前記第一のアルゴリズムによる複数の算出結果のうち前記残差が最も小さい時点の算出値を初期値とし、前記第二のアルゴリズムによる算出を開始することを特徴とする。
このプログラムを実行することによっても、上記シミュレーション方法が奏する作用効果を得ることができる。

本発明に係る粘弾性流体のシミュレーション方法を示すフローチャート。 本発明に係る粘弾性流体のシミュレーション装置を模式的に示すブロック図。 本実施形態で用いる解析対象となる流路及び計算格子に関する説明図。 第一のアルゴリズムによる微分値の算出方法に関する説明図。 第二のアルゴリズムによる微分値の算出方法に関する説明図。 第二のアルゴリズムによる微分値の算出方法に関する説明図。 第二のアルゴリズムによる微分値の算出方法に関する説明図。 第二のアルゴリズムによる微分値の算出方法に関する説明図。 第二のアルゴリズムによる微分値の算出方法に関する説明図。 第二のアルゴリズムによる微分値の算出方法に関する説明図。 第二のアルゴリズムによる微分値の算出方法に関する説明図。 数値計算における反復回数と残差の一例を示す説明図。 数値計算における反復回数と残差の他の一例を示す説明図。 数値計算における反復回数と残差の関係を一部拡大して示す説明図。 本発明の他の実施形態に係る粘弾性流体のシミュレーション方法を示すフローチャート。

以下、本発明の一実施形態を、図面を参照して説明する。

［粘弾性流体のシミュレーション装置］
本実施形態の粘弾性流体のシミュレーション装置は、押出機の流路にゴムなどの粘弾性流体が流れる場合に、定常状態における粘弾性流体の圧力、速度及び応力を含む物理量を算出するために用いられる。説明を簡素にするために、本明細書では、粘弾性流体の圧力、速度及び応力をまとめて物理量と表記する。

図２に示すように、粘弾性流体のシミュレーション装置１は、初期設定部１０と、物理値及び物理量の算出に必要な各種情報を記憶するワークメモリ１１と、物理量算出部１２と、残差算出部１３と、判定部１４と、最良物理データ記憶部１５とを有する。これら各部１０〜１５は、ＣＰＵ、メモリ、各種インターフェイス等を備えたパソコン等の情報処理装置において予め記憶されている図示しない処理ルーチンをＣＰＵが実行することによりソフトウェア及びハードウェアが協働して実現される。

図２に示す初期設定部１０は、キーボードやマウス等の既知の操作部を介してユーザからの操作を受け付け、解析対象となる流路に関する情報の設定、前記流路に対する計算格子の設定、各々の格子点に対する物理量の初期値の設定、流路の境界条件の設定等、シミュレーションに必要な各種設定を実行し、これら設定値をワークメモリ１１に記憶する。具体的には、本実施形態では、解析対象となる流路として、図３に示す流路が設定される。流路の解析領域に対して所定間隔で格子が配置される。全ての格子点に対して物理量の初期値が設定される。ここでは、流れが一様であるとして、全ての格子点に対して同じ速度を設定している。また、流路の境界条件が設定される。流入境界条件としては、例えば流路の入口部ｉｎが全て同じ流速とする条件や、流路の入口中央部の流速が速いいわゆる放物線状の流速とする条件などが挙げられる。流出境界条件としては、例えば出口ｏｕｔから大気圧に流体を開放する境界条件などが挙げられる。壁面ｈｗの境界条件としては、速度０となる完全固着条件や、或る程度の滑りを許容するスリッド条件などが挙げられる。これらの各種設定は、解析の目的に応じて適宜選択可能である。

図２に示す物理量算出部１２は、予め設定されている粘弾性流体の支配する微分方程式に基づき物理量を算出するためのもので、予め設定されている代数方程式（差分方程式）を用いて、全ての格子点について差分方程式の数値計算を行う。この代数方程式は、上述の通り、粘弾性流体の挙動を支配するナビエ・ストークス方程式（偏微分方程式）を有限差分法により離散化することで、ナビエ・ストークス方程式を近似した式である。物理量算出部１２は、格子点の物理量の微分値を差分法で算出するための第一のアルゴリズム１２ａ及び第二のアルゴリズム１２ｂを有する。

以下、第一のアルゴリズム１２ａと第二のアルゴリズム１２ｂの違いについて説明する。二次元問題の場合、粘弾性流体を支配する方程式は６つとなり、未知の物理量が、圧力Ｐ，速度ｕ，ｖ，応力Ｓ^ＸＸ，Ｓ^ＸＹ，Ｓ^ＹＹの６つになることが知られている。これら支配方程式のうち例えば応力テンソルの式を二次元問題として成分で記載すると下記（１）〜（３）の３つの式になる。
ここで、Ｒｅはレイノルズ数を表し、Ｗｅはワイゼンベルグ数を表す。

図２に示す第一のアルゴリズム１２ａは、例えば、上記式（１）の右辺の第一項及び第二項の微分値（∂ｕ／∂ｘ）、（∂ｕ／∂ｙ）を下記式（４）のように差分法で離散化し、数値計算を行うアルゴリズムである。第一のアルゴリズム１２ａにおける差分法を図示すれば図４のように表現される。例えば算出対象となる格子点を評価点と呼ぶ場合に、評価点ＡＰ_(ｉ，ｊ)の（∂ｕ／∂ｘ）は、格子点GP1の値と格子点GP2の値の差分に基づき計算され、評価点ＡＰ_(ｉ，ｊ)の（∂ｕ／∂ｙ）は、格子点GP3の値と格子点GP4の値の差分に基づき計算される。すなわち、図４に示すように、第一のアルゴリズム１２ａは、算出対象となる格子点（ｘ＝ｉ，ｙ＝ｊ）を評価点ＡＰ_(ｉ，ｊ)とする場合に、評価点ＡＰ_(ｉ，ｊ)に隣接する所定数（４つ）の格子点GP1、GP2、GP3、GP4の値（物理量）に基づいて評価点ＡＰ_(ｉ，ｊ)の微分値を差分法で算出する。なお、各々の格子点GP1、GP2、GP3、GP4の物理量（Ｐ，ｕ，ｖ，Ｓ^ＸＸ，Ｓ^ＸＹ，Ｓ^ＹＹ）は、図２に示すワークメモリ１１に記憶されている。

なお、ここでは、説明のために式の一部のみを挙げているが、他の項についても適宜離散化しているのはいうまでもない。

図２に示す第二のアルゴリズム１２ｂは、例えば、上記式（１）の右辺の第二項の微分値（∂ｕ／∂ｙ）を下記式（５）のように差分法で離散化し、数値計算を行うアルゴリズムである。第二のアルゴリズムの式（５）を図示すれば、図５Ａのように表現され、評価点ＡＰ_(ｉ，ｊ)の微分値は、仮想格子点vGP1、vGP2、vGP3、vGP4の値の平均となる。これら仮想格子点vGP1、vGP2、vGP3、vGP4は、式（５）の第１項〜第４項にそれぞれ対応する。

仮想格子点vGP1、vGP2、vGP3、vGP4には、当然ながら値がないため、仮想格子点vGP1〜4の値は次のように算出される。ここでは、仮想格子点vGP1のみを例に挙げる。仮想格子点vGP1は、下記式（６）及び図５Ｂに示すように、仮想格子点vGP5、vGP6、vGP7、vGP8の値の平均となる。仮想格子点vGP5、vGP6、vGP7、vGP8は、式（６）の第１項〜第４項にそれぞれ対応する。

さらに、仮想格子点vGP5は、下記式（７）及び図５Ｃに示すように、格子点GP1の値と格子点GP5の値の差分に基づき計算される。

同様に、仮想格子点vGP6は、下記式（８）及び図５Ｄに示すように、格子点GP4の値と格子点ＡＰ_(ｉ，ｊ)の値の差分に基づき計算される。

同様に、仮想格子点vGP7は、下記式（９）及び図５Ｅに示すように、格子点GP6の値と格子点GP1の値の差分、及び、格子点GP7の値と格子点ＡＰ_(ｉ，ｊ)の値の差分に基づき計算される。

同様に、仮想格子点vGP8は、下記式（１０）及び図５Ｆに示すように、格子点GP8の値と格子点GP5の値の差分、及び、格子点GP3の値と格子点GP4の値の差分に基づき計算される。

上記式（５）を全て展開し、ｙ方向微分だけでなく、x方向微分にも適用すれば、図５Ｇに示すように、第二のアルゴリズム１２ｂは、算出対象となる格子点（ｘ＝ｉ，ｙ＝ｊ）を評価点ＡＰ_(ｉ，ｊ)とする場合に、評価点ＡＰ_(ｉ，ｊ)の周囲にある上記４つよりも多い数（２０つ）の格子点GP1〜GP20の値（物理量）に基づいて評価点ＡＰ_(ｉ，ｊ)の微分値を差分法で算出する。すなわち、第一のアルゴリズム１２ａでは、第二のアルゴリズム１２ｂは、評価点ＡＰ_(ｉ，ｊ)の微分値を計算するために用いる格子点の数を第一のアルゴリズム１２ａに比べて広範囲にするものである。

図２に示すように、残差算出部１３は、物理量算出部１２が算出した算出値が収束状態であるか否かを判定するための指標を算出する。具体的には、残差算出部１３は、算出値Ｐ^ｎ＋１の前回値Ｐ^ｎに対する変化量を表す残差Ｒを下記式により算出する。本実施形態では、収束状態ともいえる定常状態では、時間変化に応じて物理量が変わらないので、残差Ｒが０又は限りなく小さくなり、逆に発散状態では、残差Ｒが徐々に大きくなるので、残差Ｒに基づき発散状態であるか定常状態に収束したかを判定できる。本実施形態では、残差Ｒを圧力Ｐに基づき算出しているが、物理量（圧力、速度及び応力）のいずれを採用してもよい。残差算出部１３は、物理量算出部１２により全ての格子点の物理量が算出される度に残差Ｒを算出する。

ここで、ｉ_maxは、直交座標におけるｘ方向の格子点数を表す。ｊ_maxは、直交座標におけるｙ方向の格子点数を表す。ｎは、計算の反復回数を表す。Ｐは、各格子点（ｉ，ｊ）における圧力値を表す。また、境界となる格子点の物理量（圧力、速度、応力）は固定値又は外挿値であるため、物理量の変化を観測する観点から境界となる格子点を除く内点のみを参照するようにしている。これにより、全ての格子点を参照する場合に比べて計算コストを削減できるので、好ましい。

最良物理データ記憶部１５は、物理量算出部１２が算出した複数の算出結果のうち残差Ｒが最小となる時点の物理量を記憶する。

判定部１４は、物理量算出部１２における繰り返し計算を終了するか否かを判定する。本実施形態においては、判定部１４は、第一のアルゴリズム１２ａの繰り返し計算時において、残差Ｒに基づき発散状態であると判定される場合、又は、第一のアルゴリズム１２ａによる算出回数ｎが所定回数Ｎ_max以上実行したと判定される場合の少なくともいずれかに該当するときに、終了条件が成立すると判定する。判定部１４は、残差Ｒに基づき発散状態であるか否かを判定するために、残差算出部１３が残差Ｒを算出する度に、予め記録している最小残差Ｒ_minと現残差Ｒ_ｎの比較を行い、最小残差Ｒ_min＞現残差Ｒ_ｎの場合のみ、最小残差Ｒ_min＝現残差Ｒ_ｎとして最小残差を更新する。

また、判定部１４は、現残差Ｒ_ｎと予め記録されている最小残差Ｒ_minとの比率（Ｒ_ｎ／Ｒ_min）が所定閾値ε_errを超えない場合に発散状態でないと判定し、前記比率（Ｒ_ｎ／Ｒ_min）が前記所定閾値ε_errを超える場合に発散状態であると判定する。ここで、計算が収束に向かう状態では、上記比率（Ｒ_ｎ／Ｒ_min）が１又は１に近い値となり、逆に発散する状態では、上記比率（Ｒ_ｎ／Ｒ_min）が１から遠ざかるのを利用している。所定閾値ε_errとしては、例えば１よりも十分に大きな値（例えば１．０×１０^１）などが考えられる。なお、本実施形態では、上記比率を（Ｒ_ｎ／Ｒ_min）としているが、（Ｒ_min／Ｒ_ｎ）としても良い。この場合、所定閾値ε_errは、１よりも十分に小さな値（例えば１．０×１０^−１）などに設定することが挙げられる。

また、判定部１４は、第二のアルゴリズム１２ｂの繰り返し計算時において、残差Ｒが規定値ε_０よりも小さい場合に、停止条件が成立したと判定する。規定値ε_０は、定常状態であると規定するために予め設定される値であり、例えば１．０×１０^−３が挙げられる。

［粘弾性流体のシミュレーション方法］
上記シミュレーション装置１を用いた粘弾性流体のシミュレーション方法について説明する。

まず、ステップＳＴ１（図１参照）において、図２に示す初期設定部１０は、図３に示す粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する物理量の初期値の設定、及び流路の境界条件の設定を実行する。

次のステップＳＴ２（図１参照）において、図２に示す物理量算出部１２は、第一のアルゴリズム１２ａを用いて、全ての格子点における物理量を算出する。次のステップＳＴ３において、残差算出部１３は、物理量算出部１２が計算した算出値に基づき現残差Ｒ_ｎを算出する。次のステップＳＴ４において、判定部１４は、現残差Ｒ_ｎと予め記録されている最小残差Ｒ_minとを比較し、最小残差Ｒ_min＞現残差Ｒ_ｎの場合のみ（ステップＳＴ４：ＹＥＳ）、現残差Ｒ_ｎを最小残差Ｒ_minとして記録する。その際に、ワークメモリ１１に記憶されている現時点の各物理量を最良物理データ記憶部１５に記録する（ステップＳＴ５参照）。

次のステップＳＴ６（図１参照）において、現残差Ｒ_ｎに基づき発散状態であるか否かを判定する。具体的には、比率（Ｒ_ｎ／Ｒ_min）≧所定閾値ε_errの場合（ＳＴ６：ＹＥＳ）には、終了条件が成立するとして、第一のアルゴリズム１２ａによる繰り返し計算を終了し、ステップＳＴ８に移行する。

一方、ステップＳＴ６において、比率（Ｒ_ｎ／Ｒ_min）＜所定閾値ε_errの場合（ＳＴ６：ＮＯ）には、次のステップＳＴ７において、第一のアルゴリズム１２ａによる物理量の算出を所定回数Ｎ_max以上実行した否かを判定する。具体的には、第一のアルゴリズムによる計算の反復回数ｎ≧所定回数Ｎ_maxの場合（ＳＴ７：ＹＥＳ）には、終了条件が成立するとして、第一のアルゴリズムによる繰り返し計算を終了し、ステップＳＴ８に移行する。一方、ステップＳＴ７において、反復回数Ｎ＜所定回数Ｎ_maxの場合（ＳＴ７：ＮＯ）には、ステップＳＴ２に戻り、第一のアルゴリズム１２ａによる算出を繰り返し実行する。

第一のアルゴリズム１２ａによる繰り返し計算が終了すれば、最良物理データ記憶部１５に記録されている物理量をワークメモリ１１にセットし（ステップＳＴ８参照）、第二のアルゴリズム１２ｂによる算出を開始する（ステップＳＴ９参照）。すなわち、第一のアルゴリズム１２ａによる複数の算出結果のうち残差が最も小さい時点の算出値を初期値として、第二のアルゴリズムによる算出を開始する。

次のステップＳＴ１０において、残差算出部１３は、物理量算出部１２が計算した算出値に基づき現残差Ｒ_ｎを算出する。次のステップＳＴ１１において、判定部１４は、残差Ｒ_ｎに基づき収束状態であるか否かを判定する。具体的には、残差Ｒ_ｎ≧規定値ε_０の場合（ＳＴ１１：ＮＯ）に、収束状態ではないと判定し、ステップＳＴ９に戻り、第二のアルゴリズムによる算出を繰り返し実行する。一方、判定部１４は、残差Ｒ_ｎ＜規定値ε_０の場合（ＳＴ１１：ＹＥＳ）に、収束状態である判定し、第二のアルゴリズムによる繰り返し計算を終了し、定常解を得る（ＳＴ１２参照）。

第一のアルゴリズムによる数値計算の場合、反復計算回数ｎと残差Ｒ_ｎの関係は、図６Ａ及び図６Ｂに示す２つのパターンが考えられる。第一のパターンでは、図６Ａに示すように、計算が所定回数Ｎ_maxに至るまでに発散しない場合である。この場合には、反復回数Ｎ≧所定回数Ｎ_maxの終了条件が成立し、計算回数Ｎ_maxの時点の物理量が第二のアルゴリズムの初期値として用いられる。第二のパターンでは、図６Ｂに示すように、計算が所定回数Ｎ_maxに至るまでに発散する場合である。この場合、比率（Ｒ_ｎ／Ｒ_min）≧所定閾値ε_errの終了条件が成立し、計算回数Ｎ_ｔの時点の物理量が第二のアルゴリズムの初期値として用いられる。

以上のように、本実施形態の粘弾性流体のシミュレーション方法は、粘弾性流体の圧力、速度及び応力を含む物理量に関する微分方程式に基づき定常状態における物理量を算出するために、コンピュータが実行する粘弾性流体のシミュレーション方法であって、
粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する物理量の初期値の設定、及び流路の境界条件の設定を実行するステップＳＴ１と、
算出対象となる格子点を評価点ＡＰ_(ｉ，ｊ)とする場合に、評価点ＡＰ_(ｉ，ｊ)に隣接する所定数（４つ）の格子点GP1〜GP4の値に基づいて評価点ＡＰ_(ｉ，ｊ)の微分値を差分法で算出する第一のアルゴリズム１２ａを用い、全ての格子点における物理量を算出し、算出値に基づき所定の終了条件が成立するかを判定し、所定終了条件の成立まで物理量の算出を繰り返すステップＳＴ２〜７と、
評価点ＡＰ_(ｉ，ｊ)の周囲にある所定数よりも多い数の格子点GP1〜GP20の値に基づいて評価点ＡＰ_(ｉ，ｊ)の微分値を差分法で算出する第二のアルゴリズム１２ｂを用い、全ての格子点における物理量の算出を繰り返すステップＳＴ９〜１１と、を含み、
第一のアルゴリズム１２ａにより物理量の算出を繰り返すステップＳＴ２〜７において、算出値Ｐ^ｎ＋１の前回値Ｐ^ｎに対する変動量を表す残差Ｒ_ｎに基づき発散状態であると判定される場合（ＳＴ６：ＹＥＳ）、又は、第一のアルゴリズム１２ａによる物理量の算出を所定回数Ｎ_max以上実行したと判定される場合（ＳＴ７：ＹＥＳ）の少なくともいずれかに該当するときに所定の終了条件が成立すると判定し、
終了条件が成立すると判定した場合に、第一のアルゴリズム１２ａによる複数の算出結果のうち残差Ｒ_ｎが最も小さい時点の算出値を初期値とし（ＳＴ８）、第二のアルゴリズムによる算出を開始する（ＳＴ９）。

本実施形態に係る粘弾性流体のシミュレーション装置は、粘弾性流体の圧力、速度及び応力を含む物理量に関する微分方程式に基づき定常状態における物理量を算出するために用いられる粘弾性流体のシミュレーション装置１であって、
粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する物理量の初期値の設定、及び流路の境界条件の設定を実行する初期設定部１０と、
算出対象となる格子点を評価点ＡＰ_(ｉ，ｊ)とする場合に、評価点ＡＰ_(ｉ，ｊ)に隣接する所定数（４つ）の格子点GP1〜GP4の値に基づいて評価点ＡＰ_(ｉ，ｊ)の微分値を差分法で算出する第一のアルゴリズム１２ａを用い、全ての格子点における物理量の算出を繰り返し、その後、評価点ＡＰ_(ｉ，ｊ)の周囲にある所定数（４つ）よりも多い数の格子点GP1〜GP20の値に基づいて評価点ＡＰ_(ｉ，ｊ)の微分値を差分法で算出する第二のアルゴリズム１２ｂを用い、全ての格子点における物理量の算出を繰り返す物理量算出部１２と、
物理量算出部１２による繰り返し計算の終了条件が成立したか否かを判定する判定部１４と、を備え、
判定部１４は、物理量算出部１２が第一のアルゴリズム１２ａにより物理量の算出を繰り返すステップＳＴ２〜７において、算出値Ｐ^ｎ＋１の前回値Ｐ^ｎに対する変動量を表す残差Ｒ_ｎに基づき発散状態であると判定される場合（ＳＴ６：ＹＥＳ）、又は、第一のアルゴリズム１２ａによる物理量の算出を所定回数Ｎ_max以上実行したと判定される場合（ＳＴ７：ＹＥＳ）の少なくともいずれかに該当するときに終了条件が成立すると判定し、
物理量算出部１２は、第一のアルゴリズム１２ａによる繰り返し計算の終了条件が成立すると判定された場合に、第一のアルゴリズム１２ａによる複数の算出結果のうち残差Ｒ_ｎが最も小さい時点の算出値を初期値とし（ＳＴ８）、第二のアルゴリズム１２ｂによる算出を開始する（ＳＴ９）。

収束状態ともいえる定常状態では、時間変化に応じて物理量が変わらないので、算出値Ｐ^ｎ＋１の前回値Ｐ^ｎに対する変化量を表す残差Ｒ_ｎが０又は限りなく小さくなり、逆に発散状態では、残差Ｒ_ｎが徐々に大きくなるので、残差Ｒ_ｎに基づき発散状態であるか定常状態に収束したかを判定できる。これを前提とし、算出対象となる評価点ＡＰ_(ｉ，ｊ)に隣接する所定数の格子点GP1〜GP4の値に基づいて評価点ＡＰ_(ｉ，ｊ)の微分値を差分法で算出する第一のアルゴリズムの実行時に、収束が見込めない発散状態であるとき、又は、計算を所定回数Ｎ_max以上実行したときのいずれかに該当するときに繰り返し計算を終了するようにしているので、計算が終了しないことを回避しつつ、算出値を或る程度真値に向かわせることが可能となる。それでいて、第一のアルゴリズムによる複数の算出結果のうち残差Ｒ_ｎが最も小さい時点の算出値を初期値とし、第二のアルゴリズムによる算出を開始するので、計算負荷の大きい第二のアルゴリズム単独で計算する場合に比べて全体の計算コストを削減できるとともに、確実に収束させることが可能となる。
したがって、第一及び第二のアルゴリズムの良い所取りが可能となり、高い収束性と低計算コストとを両立することが可能となる。

特に、本実施形態の方法において、第一のアルゴリズム１２ａにより物理量の算出を繰り返すステップＳＴ２〜７において、残差Ｒ_ｎのうち最小残差Ｒ_minを記録しておき、最小残差Ｒ_minが更新されるときに、全ての格子点について最小残差Ｒ_minとなる時点の物理量を記憶するステップＳＴ４，ＳＴ５を有する。
本実施形態の装置において、物理量算出部１２が第一のアルゴリズム１２ａにより物理量の算出を繰り返すステップＳＴ２〜７において、最小残差Ｒ_minが更新されるときに、全ての格子点について最小残差Ｒ_minとなる時点の物理量を記憶する最良物理データ記憶部１５を有する。

このようにすれば、終了条件が成立して第一のアルゴリズムによる算出が終了した時点で、全ての格子点について残差Ｒ_ｎが最も小さい時点の物理量が最良物理データ記憶部１５に記憶されているので、これら物理量を第二のアルゴリズムの初期値としてそのまま用いることができ、残差Ｒ_ｎが最小となる時点の物理量を再度算出する必要がなくなり、計算コストを低減することが可能となる。

さらに、本実施形態の方法では、第一のアルゴリズムにより物理量の算出を繰り返すステップＳＴ２〜７において、物理量の算出結果に基づき算出する残差Ｒ_ｎと予め記録されている最小残差Ｒ_minとの比率（Ｒ_ｎ／Ｒ_min）が所定閾値ε_errを超えない場合に発散状態でないと判定し、比率（Ｒ_ｎ／Ｒ_min）が所定閾値ε_errを超える場合に発散状態であると判定する。
本実施形態の装置では、物理量算出部１２が第一のアルゴリズムにより物理量の算出を繰り返すステップＳＴ２〜７において、判定部１４は、物理量の算出結果に基づき算出する残差Ｒ_ｎと予め記録されている最小残差Ｒ_minとの比率（Ｒ_ｎ／Ｒ_min）が所定閾値ε_errを超えない場合に発散状態でないと判定し、比率（Ｒ_ｎ／Ｒ_min）が所定閾値ε_errを超える場合に発散状態であると判定するように構成されている。

比率が所定閾値を超えるとは、値が大きくなる方向へ向けて比率が所定閾値を超える場合だけでなく、値が小さくなる方向へ向けて比率が所定閾値を超える場合も含まれる。
このように、残差Ｒ_ｎと最小残差Ｒ_minとの比率（Ｒ_ｎ／Ｒ_min）および所定閾値ε_errに応じて発散状態であるかを判定するので、複雑な計算を要することなく、発散状態であるかの判定が可能となる。それでいて、所定閾値ε_errを１よりも十分に大きい値か、１よりも十分に小さい値に設定すれば、図６Ｃに示すように、残差Ｒ_ｎが脈動しながら収束する場合に、誤って発散状態であると判定することを低減することが可能となる。

本実施形態に係る粘弾性流体のシミュレーションプログラムは、上記粘弾性流体のシミュレーション方法を構成する各ステップをコンピュータに実行させるプログラムである。このプログラムを実行することによっても、上記シミュレーション方法の奏する作用効果を得ることが可能となる。言い換えると、上記シミュレーション方法を使用しているとも言える。

以上、本発明の実施形態について図面に基づいて説明したが、具体的な構成は、これらの実施形態に限定されるものでないと考えられるべきである。本発明の範囲は、上記した実施形態の説明だけではなく特許請求の範囲によって示され、さらに特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれる。

例えば、図２に示す各部１０〜１５は、所定プログラムをコンピュータのＣＰＵで実行することで実現しているが、各部を専用メモリや専用回路で構成してもよい。

［他の実施形態］
また、本実施形態では、第一のアルゴリズムによる計算の後、第二のアルゴリズムの繰り返し計算を行い、定常解を得ているが、第二のアルゴリズムの繰り返し計算により或る程度収束した時点で、第一のアルゴリズムに切り替え、当該第一のアルゴリズムによる計算で定常解を得るように構成してもよい。すなわち、図７に示すように、物理量算出部１２が第二のアルゴリズムにより物理量の算出を繰り返すステップＳＴ９〜１１において、判定部１４は、定常状態であると規定する残差の規定値ε_０（例えば１．０×１０^−３）よりも大きな値の基準値ε_１（例えば１．０×１０^−２）を予め設定しておき、残差Ｒ_ｎと基準値ε_１とを比較する（ステップＳＴ１１）。残差Ｒ_ｎが基準値ε_１よりも小さいと判定部１４が判定した場合（ＳＴ１１：ＹＥＳ）に、物理量算出部１２は第二のアルゴリズム１２ｂによる繰り返し計算を終了する。その後、物理量算出部１２は、第二のアルゴリズムによる算出結果を初期値とし、再度、第一のアルゴリズムによる繰り返し計算を開始する（ＳＴ１３参照）。定常状態であると規定する残差の規定値ε_０より残差が小さくなると判定部１４が判定した場合に（ＳＴ１４：ＹＥＳ）、物理量算出部１２は第一のアルゴリズムによる繰り返し計算を終了し、定常状態における物理量を算出する（ＳＴ１２）。

このようにすれば、計算コストの高い第二のアルゴリズムで定常状態になるまで計算せずに、定常状態になる前に、計算コストの低い第一のアルゴリズムでの計算に移行するので、計算コストを低減させることが可能となる。それでいて、第二のアルゴリズムによる算出結果を初期値とするので、初期値が或る程度収束しており、第一のアルゴリズムを用いても適切に収束させることが可能な場合が多い。
したがって、高い収束性を確保したまま、計算コストをより一層削減することが可能となる。

上記の各実施形態で採用している構造を他の任意の実施形態に採用することは可能である。各部の具体的な構成は、上述した実施形態のみに限定されるものではなく、本発明の趣旨を逸脱しない範囲で種々変形が可能である。

１０…初期設定部
１２…物理量算出部
１２ａ…第一のアルゴリズム
１２ｂ…第二のアルゴリズム
１４…判定部
１５…最良物理データ記憶部
ＡＰ_{（ｉ，ｊ）}…評価点
Ｐ^ｎ＋１…算出値
Ｐ^ｎ…前回値
Ｒ_min…最小残差
Ｒ_ｎ…残差
（Ｒ_ｎ／Ｒ_min）…残差と最小残差の比率
ε_err…所定閾値
ε_０…定常状態であると規定する残差の規定値
ε_１…基準値

Claims (9)

1. 粘弾性流体の圧力、速度及び応力を含む物理量に関する微分方程式に基づき定常状態における前記物理量を算出するために、コンピュータが実行する粘弾性流体のシミュレーション方法であって、
   前記粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する前記物理量の初期値の設定、及び前記流路の境界条件の設定を実行するステップと、
   算出対象となる格子点を評価点とする場合に、前記評価点に隣接する所定数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第一のアルゴリズムを用い、全ての格子点における前記物理量を算出し、算出値に基づき所定の終了条件が成立するかを判定し、当該所定終了条件の成立まで前記物理量の算出を繰り返すステップと、
   前記評価点の周囲にある前記所定数よりも多い数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第二のアルゴリズムを用い、全ての格子点における前記物理量の算出を繰り返すステップと、を含み、
   前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、算出値の前回値に対する変動量を表す残差に基づき発散状態であると判定される場合、又は、前記第一のアルゴリズムによる物理量の算出を所定回数以上実行したと判定される場合の少なくともいずれかに該当するときに前記所定の終了条件が成立すると判定し、
   当該終了条件が成立すると判定した場合に、前記第一のアルゴリズムによる複数の算出結果のうち前記残差が最も小さい時点の算出値を初期値とし、前記第二のアルゴリズムによる算出を開始することを特徴とする粘弾性流体のシミュレーション方法。
2. 前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、前記残差のうち最小残差を記録しておき、前記最小残差が更新されるときに、全ての格子点について前記最小残差となる時点の物理量を記憶する請求項１に記載の粘弾性流体のシミュレーション方法。
3. 前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、前記物理量の算出結果に基づき算出する残差と予め記録されている最小残差との比率が所定閾値を超えない場合に発散状態でないと判定し、前記比率が前記所定閾値を超える場合に発散状態であると判定する請求項１又は２に記載の粘弾性流体のシミュレーション方法。
4. 前記第二のアルゴリズムにより物理量の算出を繰り返すステップにおいて、定常状態であると規定する残差の規定値よりも大きな値の基準値を予め設定しておき、前記残差と前記基準値とを比較し、前記残差が前記基準値よりも小さいと判定した場合に、前記第二のアルゴリズムによる繰り返し計算を停止し、その後、前記第二のアルゴリズムによる算出結果を初期値とし、再度、前記第一のアルゴリズムによる繰り返し計算を開始し、定常状態であると規定する残差の規定値より残差が小さくなることを以て、繰り返し計算を終了して定常状態における物理量を算出する請求項１〜３のいずれかに記載の粘弾性流体のシミュレーション方法。
5. 粘弾性流体の圧力、速度及び応力を含む物理量に関する微分方程式に基づき定常状態における前記物理量を算出するために用いられる粘弾性流体のシミュレーション装置であって、
   前記粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する前記物理量の初期値の設定、及び前記流路の境界条件の設定を実行する初期設定部と、
   算出対象となる格子点を評価点とする場合に、前記評価点に隣接する所定数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第一のアルゴリズムを用い、全ての格子点における前記物理量の算出を繰り返し、その後、前記評価点の周囲にある前記所定数よりも多い数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第二のアルゴリズムを用い、全ての格子点における前記物理量の算出を繰り返す物理量算出部と、
   前記物理量算出部による繰り返し計算の終了条件が成立したか否かを判定する判定部と、
   を備え、
   前記判定部は、前記物理量算出部が前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、算出値の前回値に対する変動量を表す残差に基づき発散状態であると判定される場合、又は、前記第一のアルゴリズムによる物理量の算出を所定回数以上実行したと判定される場合の少なくともいずれかに該当するときに終了条件が成立すると判定し、
   前記物理量算出部は、第一のアルゴリズムによる繰り返し計算の終了条件が成立すると判定された場合に、前記第一のアルゴリズムによる複数の算出結果のうち前記残差が最も小さい時点の算出値を初期値とし、前記第二のアルゴリズムによる算出を開始することを特徴とする粘弾性流体のシミュレーション装置。
6. 前記物理量算出部が前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、最小残差が更新されるときに、全ての格子点について前記最小残差となる時点の物理量を記憶する最良物理データ記憶部を有する請求項５に記載の粘弾性流体のシミュレーション装置。
7. 前記物理量算出部が前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、前記判定部は、前記物理量の算出結果に基づき算出する残差と予め記録されている最小残差との比率が所定閾値を超えない場合に発散状態でないと判定し、前記比率が前記所定閾値を超える場合に発散状態であると判定するように構成されている請求項５又は６に記載の粘弾性流体のシミュレーション装置。
8. 前記物理量算出部が前記第二のアルゴリズムにより物理量の算出を繰り返すステップにおいて、前記判定部は、定常状態であると規定する残差の規定値よりも大きな値の基準値を予め設定しておき、前記残差と前記基準値とを比較し、
   前記残差が前記基準値よりも小さいと前記判定部が判定した場合に、前記物理量算出部が前記第二のアルゴリズムによる繰り返し計算を停止し、その後、前記第二のアルゴリズムによる算出結果を初期値とし、再度、前記第一のアルゴリズムによる繰り返し計算を開始し、定常状態であると規定する残差の規定値により残差が小さくなると前記判定部が判定する場合に、前記物理量算出部は繰り返し計算を終了して定常状態における物理量を算出する請求項５〜７のいずれかに記載の粘弾性流体のシミュレーション装置。
9. 粘弾性流体の圧力、速度及び応力を含む物理量に関する微分方程式に基づき定常状態における前記物理量を算出するために、下記ステップをコンピュータに実行させる粘弾性流体のシミュレーションプログラムであって、
   前記粘弾性流体が流れる流路に対する計算格子の設定、各々の格子点に対する前記物理量の初期値の設定、及び前記流路の境界条件の設定を実行するステップと、
   算出対象となる格子点を評価点とする場合に、前記評価点に隣接する所定数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第一のアルゴリズムを用い、全ての格子点における前記物理量を算出し、算出値に基づき所定の終了条件が成立するかを判定し、当該所定終了条件の成立まで前記物理量の算出を繰り返すステップと、
   前記評価点の周囲にある前記所定数よりも多い数の格子点の値に基づいて前記評価点の微分値を差分法で算出する第二のアルゴリズムを用い、全ての格子点における前記物理量の算出を繰り返すステップと、を含み、
   前記第一のアルゴリズムにより物理量の算出を繰り返すステップにおいて、算出値の前回値に対する変動量を表す残差に基づき発散状態であると判定される場合、又は、前記第一のアルゴリズムによる物理量の算出を所定回数以上実行したと判定される場合の少なくともいずれかに該当するときに前記所定の終了条件が成立すると判定し、
   当該終了条件が成立すると判定した場合に、前記第一のアルゴリズムによる複数の算出結果のうち前記残差が最も小さい時点の算出値を初期値とし、前記第二のアルゴリズムによる算出を開始することを特徴とする粘弾性流体のシミュレーションプログラム。