JP5132768B2

JP5132768B2 - モールド充填プロセスのシミュレーションにおいて、粒子の統計的な配向分布を記述するための方法および装置。

Info

Publication number: JP5132768B2
Application number: JP2010513774A
Authority: JP
Inventors: ヨアキムリン; マティアスモーグ
Original assignee: MAGMA Giessereitechnologie GmbH
Current assignee: MAGMA Giessereitechnologie GmbH
Priority date: 2007-07-02
Filing date: 2008-07-01
Publication date: 2013-01-30
Anticipated expiration: 2028-07-01
Also published as: KR20100036226A; EP2481550A1; ES2389098T3; JP2010531752A; US20100169062A1; EP2164694B1; EP2481550B1; CN101754845A; US8364453B2; KR101524260B1; CN101754845B; EP2164694A1; WO2009003677A1

Description

本発明は、粒子を含む懸濁液（ｓｕｓｐｅｎｓｉｏｎ）のモールドキャビティへのフロー（ｆｌｏｗ）に関する３次元モデリングの分野に関し、より具体的には、多数の非球形粒子を含む懸濁液でモールドキャビティが充填されるプロセスに関するシミュレーションにおいて、非球形粒子の統計的な配向分布を記述するための方法および装置に関する。

射出成形プロセスまたは金属鋳造プロセスに関する忠実な３Ｄシミュレーションは、最大で数十万個もの数式系を伴う。これらの複雑な計算に対処するために、シミュレーション方法の効率の向上が進められてきている。最適化ソフトウェアおよび最新のワークステーションの処理能力によって、このようなシミュレーションを普通の職場で実行することが可能となり、十分に適切な結果を得られるようになった。その結果は、純粋に科学的な研究分野ではともかく、射出成形品の研究開発部や、鋳物工場、製造業者の技術者によっては十分に実用的である。

ポリマーの組成に繊維等の粒子が加えられ、その繊維の配向分布を記述する必要がある場合、シミュレーションおよび付随する式の集合は、さらに大幅に複雑になる。このようなプロセスに関する忠実な３Ｄシミュレーションは、これまで、計算時間が長過ぎるか、またはシミュレーションの精度が不十分となる恐れがあるため、シミュレーションの複雑性により、最新のワークステーションにおいても結果が許容できないことから、一般の職場への導入に成功していない。

繊維で強化された部分では、開発技術者が、部品の張力（ｔｅｎｓｉｏｎ）や整経（ｗａｒｐｉｎｇ）について予測できるように、繊維配向分布の記述を有することが極めて重要であることが多い。繊維は、通常、プラスチック部分の力学的特性を改善するために使用される。しかしながら、（熱）力学的特性（熱膨張、強度、および剛性等）は、繊維の向きに依存する。

射出成形プラスチック部品の使用は、近年、多くの産業において着実に増加している。電子機器、消費財、医療機器、および自動車部品の製造業者は、かつてない程、プラスチックから作製される製品に使用する製品および部品の生産を増やしている。

繊維で強化された射出成形部品は、改良された強度／重量比率、耐久性、部品の一体化、および低コストを提供するため、構造用金属部品に取って代わりつつある。

同時に、競争の圧力により、プラスチック射出成形産業の製造業者は、設計を生産プロセスにさらに適合させるために、設計を最適化する新しい方法を探索することを余儀なくされている。部品やモールド構成を修正する必要性が、設計開発プロセスの後期に発見される場合、必要な変更を実行するための遅延コストおよび関連コストは、設計開発段階の初期段階よりも大幅に高くなる。

自社の部品が生産可能であること、かつ最適に機能することを確実にしたい会社は、製造プロセスを十分に理解し、かつこの知識を設計段階の初期において部品設計に組み込むために、コンピュータ支援エンジニアリング技術を使用して、射出成形モールドにおける複雑なフローや結果として生じる繊維配向をシミュレーションおよびモデル化することを望む。

繊維強化部品のための、およびその中で生産される繊維強化部品のための、射出成形モールドの設計には、考慮すべき多くの因子が存在する。部品全体の幾何学的形状、最小および最大壁厚さ、モールドにおいて、液体ポリマーや繊維懸濁液が射出されるゲートの数および位置、モールドにおいて、キャビティ中の気体を排出するベントの数および位置、ポリマーの組成や特性、繊維の特性や量、鋳引け、余裕部分、繊維配向分布等のパラメータがある。部品やモールドの設計は、最終部品の形式および機能に単に大きく依存するわけではなく、それに密接に相関する製造プロセスの影響も考慮すべきである。

便利なことに、コンピュータ支援エンジニアアリングシミュレーション（ＣｏｍｐｕｔｅｒＡｉｄｅｄＥｎｇｉｎｅｅｒｉｎｇシミュレーション；ＣＡＥシミュレーション）は、射出成形プロセス中にモールドキャビティの内部で発生し得ることに関して、視覚的および数的なフィードバックを、設計および製造の技術者に提供するように用いられることができる。これによって、技術者は、想定される部品設計の挙動を十分に理解および予測することが可能になり、結果的に、製造に対する従来の高価な試行錯誤の手法を、実質的に排除することが可能となる。コンピュータ支援エンジニアアリングシミュレーションの使用によって、設計段階中の部品設計やモールド設計、製造プロセスパラメータの最適化が容易になり、設計段階中に、最小コストおよびスケジュールに対する最小の影響で、必要な変更を容易に実行することが可能になる。

繊維強化部品のエンジニアリングプロセス内におけるＣＡＥシミュレーション技術の適用は、（ｉ）流体フローおよび熱伝導の計算を含む「射出成形」製造プロセスのシミュレーション、ならびに（ｉｉ）応力や強度、場合によっては耐久性の計算を含み、これらの全ては、これらの部品について巨視的なレベルで実行され、外部負荷の存在下におけるその機能的な力学的特性を決定する。両方の種類のシミュレーションは、液体状態および固体状態において浸漬繊維を含むポリマー材料の特性を記述する適切な材料モデルを必要とする。

巨視的レベルにおける長さのスケールは、部品の幾何学的形状の線形サイズ（全体のサイズ、壁厚さ等）によって決定され、典型的には、数ｍｍから最大ｃｍまでの範囲において変動する。計算セルのサイズは、十分正確に、巨視的な長さのスケールを分解しなければならないため、最大でも、最小の巨視的サイズより約１桁小さい大きさである。短繊維で強化された部分における、混入された繊維の典型的なサイズは、巨視的計算セルの典型的なサイズよりも１桁または２桁小さいため、巨視的材料挙動のモデリングに関連する繊維特性は、統計的な手法によって記述される。短繊維強化材料では、それに関連する巨視的特性は、（ａ）通常は全体的に（ほぼ）一定である、体積濃度（ＶｏｌｕｍｅＣｏｎｃｅｎｔｒａｔｉｏｎ）と、（ｂ）通常は部分的な幾何学形状によって大幅に変化する、各セル内における繊維配向（ＦｉｂｅｒＯｒｉｅｎｔａｔｉｏｎ；ＦＯ）の局所的な分布である（本トピックに関するさらなる詳細は、発明を実施するための形態の節１．１および１．２において論じられる）。

局所的ＦＯの統計的分布に関する単純化された記述および実際的な目的における適切な記述は、対応する分布関数の低次（すなわち２次および４次）モーメントによって提供される。その数学的構造により、これらのモーメントは、（それぞれ２次および４次の）配向テンソル（ｏｒｉａｔａｔｉｏｎｔｅｎｓｏｒ）で示される。繊維強化部品のＣＡＥシミュレーションのフレームワーク内において、４次テンソルは、巨視的レベルにおける繊維強化材料のレオロジー特性（ｒｈｅｏｌｏｇｉｃａｌｐｒｏｐｅｒｔｙ）ならびに力学的特性を予測するのに必要とされる。それは、これらが４次テンソル特性を有するためである。２次ＦＯテンソルは、単位トレースを有する実数値の対称的な３×３行列であるため、その９個の成分のうち５つのみが独立している。４次ＦＯテンソルの独立成分の数は、（全）対称性（ｔｏｔａｌｓｙｍｍｅｔｒｙ）によって、３^４＝８１から１５に減少する。

そのモーメントに関してＦＯ分布を記述する数学的モデルは、閉包関係（ｃｌｏｓｕｒｅｒｅｌａｔｉｏｎ）の観点から、２次テンソルに関する４次配向テンソルを近似計算することによって大幅に単純化される。閉包関係は、関数の観点から、このような計算スキームの数学的記述を提供する。関連する計算手続は、特定の仮定下において閉包関係が近似的に有効である場合にのみ、「閉包近似」と言われる。いくつかの閉包近似とともに２次ＦＯテンソルのみを使用する手法は、モールド充填プロセス中の時間および空間における２次ＦＯテンソルの進化をシミュレーションするための、「Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ」型のモデルを導く（詳細は、後の２．２節から２．５節を参照）。

非対称性熱可塑性プラスチックの流動に関する非特許文献１は、繊維配向予測アルゴリズムを含む分析を使用して、強化熱可塑性プラスチック製品内における繊維配向分布をモデル化するための方法を開示する。ソフトウェアは、線形３角形要素から成る２次元有限要素メッシュを使用して、３次元モデルを近似している。繊維配向を計算するために、一般的なＨｅｌｅ−Ｓｈａｗ近似およびＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式の変形を使用して、流動場が計算されている。繊維配向、温度、粘度の計算は、３次元解を得るために、各要素の「厚さ」における１９の層において有限差分技術を使用して実行されている。しかしながら、モデルが、面方向の外においては速度成分をシミュレーションできないため（Ｈｅｌｅ−Ｓｈａｗ近似の限界）、このシステムが、忠実な３次元としては記述可能ではないことに留意することが重要である。この文献に説明する方法は、忠実な３次元シミュレーションに適用する場合に、不安定で、要求される処理能力が高過ぎるシミュレーションとなる可能性があり、開発技術者の普通の環境においては使用され得ないシミュレーションとなりえる。

Ｇｌａｓｓｆｉｂｒｅｏｒｉｅｎｔａｔｉｏｎｗｉｔｈｉｎｉｎｊｅｃｔｉｏｎｍｏｕｌｄｅｄａｕｔｏｍｏｔｉｖｅｐｅｄａｌ−Ｓｉｍｕｌａｔｉｏｎａｎｄｅｘｐｅｒｉｍｅｎｔａｌｓｔｕｄｉｅｓ，Ｂ．Ｒ．Ｗｈｉｔｅｓｉｄｅら，Ｐｌａｓｔｉｃｓ，ＲｕｂｂｅｒａｎｄＣｏｍｐｏｓｉｔｅｓ，２０００，Ｖｏｌｕｍｅ２９，Ｎｏ．１

このような背景から、本発明の目的は、多数の非球形粒子を含む懸濁液でモールドキャビティが充填されるプロセスのシミュレーションにおいて、非球形粒子の配向分布を決定するための方法を提供することにある。この方法は、従来の方法よりも安定し、かつ多くの計算を必要としない。

上記目的は、シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算するための、コンピュータにより実装される方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、（ａ）境界条件を特定することと；（ｂ）初期条件を設定すること；（ｃ）巨視的レベルで流体フロー、熱フロー、および質量輸送を求めるために、前記シミュレーション領域の前記セルの少なくとも一部分に関する質量・運動量・エネルギーの平衡方程式（ｂａｌａｎｃｅｅｑｕａｔｉｏｎ）を解くこと；（ｄ）空間および時間の関数として、巨視的レベルで前記非球形粒子配向における変化を決定するために、前記解かれた平衡方程式の結果に少なくとも部分的に基づいて、非球形粒子配向動的方程式（ｎｏｎｓｐｅｒｉｃａｌｐａｒｔｉｃｌｅｏｒｉｅｎｄａｔｉｏｎｄｙｎａｍｉｃｅｑｕａｔｉｏｎ）を解くことと、を含む方法を提供することによって、請求項１に従って達成される。なお、ステップｄ）において、前記粒子配向方程式は、懸濁液を含む計算セルについてのみ解かれてもよい。

ステップｃ）は、（ｃｃ）前記流体または懸濁液の更新された自由表面または流頭（ｆｌｏｗｆｒｏｎｔ）を決定することをさらに含むことが好ましい。この自由表面は、前記解かれた平衡方程式の結果に少なくとも部分的に基づいて、前記懸濁液で充填されたセルを、前記キャビティの空セルから分離することが好ましい。ステップ（ｃｃ）は、前記更新された流頭に従って、前記境界条件を更新することをさらに含むことが好ましい。また、本発明の方法は、ｅ）前記モールドキャビティが、前記懸濁液のシミュレーション射出によって充填されているかを判断することによって、前記シミュレーションされた射出成形プロセスが完了しているかを決定することと、ｆ）前記シミュレーションされた射出成形プロセスが完了するまでステップｃ），ｃｃ），ｄ）を繰り返すこととをさらに含むことが好ましい。

本発明の上記目的は、多数の非球形粒子を含む溶媒により形成される懸濁液でモールドキャビティが充填されるプロセスのシミュレーションにおいて、非球形粒子の統計的分布配向を記述するための方法であって、（１）前記キャビティの幾何学的形状を定義する３次元コンピュータモデルを提供することと、（２）境界条件を特定することと、（３）複数のセルによるメッシュを形成するために、前記モデルに基づいて解領域を離散化することと、（４）前記解領域の少なくとも一部分に関するエネルギー方程式およびフロー方程式を解くことと、（５）各々の前記セルにおけるフロー条件および温度条件を時間の関数として計算することと、（６）非球形粒子配向の変化を計算することと、（７）各々の前記セルにおける前記非球形粒子の前記配向の統計的分布を時間の関数として記述することとを含む方法を提供することによって、請求項６に従って達成される。

本発明の方法は、部分的な繊維配向分布および熱力学的特性を大幅に低減された計算量で予測するために、懸濁液のフロー（ｆｌｏｗ）を利用する。シミュレーションの結果は、開発技術者によって使用されることができ、それによって、繊維配向に関して製品が最適化され、製品の強度および形状安定性が改善されうるであろう。

多数の非球形粒子を含む懸濁液でモールドキャビティが充填されるプロセスに関する、本発明に係るシミュレーションにおいて、非球形粒子の配向分布を決定するための方法および装置のさらなる目的、特徴、利点、および特性は、次の詳細説明によって、より明らかになるであろう。

本説明の以下の詳細部分において、図面に示す例示的実施形態を参照して、本発明についてより詳細に説明する。
射出成形装置およびモールドの断面の概略図である。本発明のある実施形態に従う、射出成形プロセスの基本プロセスステップを概説するトップレベルのフローチャートである。図２のフローチャートのステップ５をさらに詳細に示すフローチャートである。繊維強化プラスチック製品の顕微鏡写真である。本発明のある実施形態において使用される繊維のモデルを示す。特性多項式（ｃｈａｒａｃｔｅｒｉｓｔｉｃｐｏｌｙｎｏｍｉａｌ）の２次項および３次項ならびに線形項および定数項を示すグラフである。行列のトリプレットの固有値を示すグラフである。「最大限に対称的である」場合に対応する特別な例を示す。４面体形状の体積の変形バージョンを示す。本発明のある実施形態に従うＦＯ行列（すなわち２次ＦＯテンソル）の位相空間集合Ｍ_ＦＴの全体構造に関する図を示す。本発明のある実施形態に従う単純化形式における演算子分割プロセスを示すフローチャートである。時間ステップ方法の詳細を示すフローチャートである。本発明のある実施形態に従う、トレース再スケーリング（ｔｒａｃｅｒｅｓｃａｌｉｎｇ）を使用する位相空間射影プロセスを示すフローチャートである。

詳細な説明

図１は、射出成形装置１を図式的に示す。射出成形装置には、ホッパー（ｈｏｐｐｅｒ）３に配置されるポリマーペレットを備えるネジ２が設けられる。ポリマーペレットは、ネジ２および加熱要素４の作用によって、粘性の塊に変換され、高圧下でモールド６のキャビティ５に注入される。また、成形装置および射出成形製造サイクルは、当技術分野において周知であり、本明細書において詳細には説明しない。成形装置１は、非繊維強化・繊維強化のどちらのプラスチック部品も製造することができる。

コンピュータによる射出成形プロセスの計算シミュレーションは、図２に示すプロセスに従って実行できる。

一般的に表されるシミュレーションの主要なステップは以下の通りである。
ステップ１：シミュレーション領域の幾何学的形状に関するデジタル表現を準備する。
ステップ２：計算領域を多数の小要素への分割するメッシュ化。これら多数の小要素は、（様々な解法アルゴリズムを利用して、）微分方程式を離散化するための基礎である。シミュレーションされる物理的現象に対する解は、このようにして探索される。
ステップ３：異なる材料領域に関して必要な物理的データをシミュレーションモデル（データベースまたはデータバンク）に取り入れる。
ステップ４：シミュレーション計画の境界条件を特定する。
ステップ５：射出成形プロセスをシミュレーションする（このステップは後に詳細に説明される）。
ステップ６：ワークステーション等のコンピュータのディスプレイ上に、グラフ提示または数値的提示として結果を表示する。

繊維強化製品の射出成形プロセスをシミュレーションするステップ５の詳細について、図３のフローチャートに示す。プロセスのこの部分では、熱フロー、流体フロー、応力、歪みに関する微分方程式が、以下の計算アルゴリズムを使用して解かれる。
ステップ１、熱物理材料特性および流頭（ｆｌｏｗｆｒｏｎｔ）の初期条件が設定される。
ステップ２、領域全体の熱方程式および全流体セルに関するフロー方程式が、質量、エネルギー、運動量保存の式を使用して解かれる。
ステップ３、このステップにおいては流頭が移動し、新しい流頭に従って境界条件が選択される。
ステップ４では、繊維配向および繊維輸送（ｆｉｂｅｒｔｒａｎｓｐｏｒｔ）が、先行するステップにおいて得られた流速から計算される。このステップは、後で詳細に説明される。初期条件が、新しく充填されたセルに適用される。流体物質を含むセルのみが考慮される。
ステップ５では、化学反応等の追加の量（ｑｕａｎｔｉｔｙ）が計算され、セルが凝固する場合に検証される。
ステップ６では、モールドにおいて射出成形プロセスが完了しているか否かが検証され、完了していない場合、シミュレーションは次の時間ステップへの進み、プロセスはステップ２に戻る。
ステップ７では、モールドから取り出した成型品の特性が計算される。
モールドの外部の熱物理材料温度、すなわち、温度、鋳引け、反り等が、射出成形シミュレーションから得た情報を使用して計算される。

図４は、繊維強化プラスチック製品の顕微鏡写真であり、射出成形プロセス完了後の繊維の配向が示されている。最終生産物における繊維の向きは、射出成形プロセス中の塊状熱可塑性プラスチックのフローパターンに大きく依存する。繊維の向きは、各々の繊維について厳密に決定されるのではなく、むしろ分布関数によって記述される。

繊維の向きの計算に関する詳細に進む前に、短繊維熱可塑性プラスチック溶融液の基礎的な流体力学的側面に関して概説する。

１．短繊維熱可塑性プラスチック溶融液の基礎的な流体力学的側面

ポリマーの塊（ｐｏｌｙｍｅｒｍａｓｓ）の中に、多数の短繊維（ｓｈｏｒｔｆｉｂｅｒ）を分散させて混入することにより、成型品を短繊維強化熱可塑性プラスチックとすることができる。融解状態では、すなわち熱可塑性プラスチックマトリクスが液体であるのに十分高い温度である場合、プラスチック溶融液と分散繊維との混合液は、流体力学および流体物理学の専門用語において一般的に粒子懸濁液（ｐａｒｔｉｃｌｅｓｕｓｐｅｎｓｉｏｎ）と言われる複合流体を構成する。一般に、このような懸濁液は、２つの異なる相、つまり（ｉ）分散繊維を含まない融解プラスチック材料に対応する溶媒相、および（ｉｉ）溶媒に浸漬される全ての繊維から成る粒子相から構成される。

１．１繊維の幾何学的形状および材料特性

回転対称（図５に示す）の軸を有するスフェリカルな粒子は、以下において繊維（ｆｉｂｅｒ）と記述される。用語の粒子（ｐａｒｔｉｃｌｅ）および繊維（ｆｉｂｅｒ）は、別途明示的に示す場合を除

ラスチック材料に分散される繊維（例えば、炭素またはガラス型の繊維）に関するこれら

比は大きくなる。

繊維強化プラスチックに浸漬される大部分の種類の繊維では、繊維材料の質量密度は、懸濁プラスチック溶融液の密度に似ている。プラスチック溶融液自体は、モールド充填プロセス中に生じる典型的な温度においてかなり高い粘度を有する。これらの２つの側面が組み合わされるので、浮力効果や慣性効果は、好適な実施形態に従う方法において、プラスチック溶融液内における繊維の動きを記述する際に無視することができる。

繊維濃度に関し、懸濁液における繊維は空間的に均一の分布を有すると仮定され、このため濃度は一定であると解釈される。しかし、（ｉ）入口に既に存在する不均一な分布の純枠に対流的な輸送、（ｉｉ）モールド充填段階中のせん断誘起粒子移動効果の発生、という作用は、これらを変化させることがある。

（ｉ）または（ｉｉ）の効果が大きい場合、熱可塑性プラスチック溶融液に浸漬される短繊維の懸濁液フローを記述するために使用されるモデルは、必然的に、粒子および溶媒の両方の相のフローを相互結合し、かつ不均一な粒子濃度を可能にする「２相フロー」モデルでなければならない。しかしながら、これはこの例に当てはまらない。現象（ｉｉ）は拡散プロセスの一種として理解可能であり、この場合、拡散定数は、懸濁粒子の実サイズ（すなわち、繊維の長さ、または球形粒子の場合は半径）の、フローキャビティ（例えば、壁厚）の典型的なサイズに対する比率として規定される相対的な粒子サイズの２乗の大きさである。相対的な粒子サイズは、繊維強化プラスチックからの射出成形によって製造される薄い壁部分であっても、典型的には０．１から０．０１のオーダーであるため、この場合、（ｉｉ）を完仝に無視することは全く妥当である。（ｉ）に関し、ネジにおける溶融液の準備プロセスにより、入口における不均一な濃度プロファイルが存在し得ることが想定され得る。しかしながら、典型的な短繊維強化プラスチック部分について実験的に実際に認められる繊維濃度は、その部分内のほぼどこにおいても一定値から極めて小さい偏移しか示さない。従って、繊維分布が均一であるとの仮定は妥当である［２０］。

１．２レオロジー特性に対する繊維の影響

熱可塑性プラスチック材料は、非ニュートン的なフロー挙動を呈する。繊維を含まない純枠なプラスチック溶融液の場合、材料のレオロジー特性（ｒｈｅｏｌｏｇｉｃａｌｐｒｏｐｅｒｔｙ）は、スカラー粘度関数（ｓｃａｌａｒｖｉｓｃｏｓｉｔｙｆｕｎｃｔｉｏｎ）によってモデル化される。このスカラー粘度関数は、温度ならびにフロー（一般的ニュートン流体）の状態変数に依存し得る。レオロジー特性は、非ニュートン的なフロー特性（例えば、せん断減粘（ｓｈｅａｒｔｈｉｎｎｉｎｇ））の限られた範囲にしか及ばないことが知られているが、射出成形シミュレーションの目的では有用かつ満足しうる正確性を有することが証明されている。この種類のモデルは、本発明の好滴な実施形態に従う方法においても使用される。

熱可塑性プラスチックマトリクス材に繊維を添加することによって、固体状態において、異方性が強く、埋め込まれた繊維の配向の局所的な分布に強く依存する力学的特性がもたらされる。原理上は、材料が液体（すなわち、融解）状態である場合にも、異方的な材料挙動が示される。この異方性を説明するためには、上述のスカラー粘度を粘度テンソルに置き換えなければならない。繊維懸濁液の材料挙動をモデル化する両方の方法を比較する多くの研究が行なわれている。モールド充填中に見られるほとんどの流動状況では、両方の種類のモデル（すなわち、スカラーモデルおよびテンソルモデル）によって予測される充填パターンが、ほとんど差を示さないことが分かる（例えば、［２１］、［２２］を参照）。ゆえに、異方的な粘度効果は無視され、熱可塑性プラスチック溶融液が繊維を含む場合においても、単純な一般的ニュートンモデルが使用される。これは、溶融液のレオロジー特性に対する分散繊維の配向の影響を無視することに等しい。粘度は、繊維濃度に依存するが、これは一定であると仮定される（上記参照）ことから、この側面は、可塑性プラスチック溶融液の効果的な材料特性に寄与する先験的な既知の材料パラメータとしてのみ関与する。

このアプローチのため、好適な実施形態に従う方法は、フローの計算と繊維配向の計算とを部分的に分けて行う。流速は、分散繊維の向きに局所的に影響を及ぼすが、フローに対する繊維配向の影響はごくわずかである。ゆえに、フローの計算は、繊維配向の計算とは独立に実行される。このように、溶融液の局所的な流速は、外部係数の組として繊維配向の計算に使用されるモデルに関与する。

２．Ｆｏｌｇａｒ−Ｔｕｃｋｅｒモデル

２．１Ｊｅｆｆｅｒｙの方程式

細長い粒子である。繊維は、溶媒の局所フロー場において繊維に作用する力学的力が、実質的にいかなる変形をも与えることができないほど十分小さい。ゆえに、個々の繊維は、本明細書において、細長く、回転対称の剛体であるようにモデル化され、繊維の向きは、回転対称軸に沿う向きを有する単位ベクトルｐによって求められる。２つのベクトル±ｐの両方は、粒子の同一の配向状態を表す。

図５は、回転楕円形状の剛体粒子を示す。この粒子の動きは、粒子の近傍における速度

に加え、配向ベクトルｐが示されている。流速Ｕの方向およびサイズは、対応する矢印の方向および長さによって示される。速度ベクトルの方向は、全く同じであるが、その長さは同じではなく、粒子近傍の流速が一定でないことが示される。しかしながら、底部から上部へのベクトルの長さ方向において増加の程度が一定であるため、流速の変化の局所的な大きさ、すなわち速度勾配は、同一であるように思われる。粒子表面の各点において、粒子は、局所的な流速と同じように移動するはずである（「滑り無し」境界条件）。粒子の周辺の流速が一定であれば、これは、単純な並進運動（ｔｒａｎｓｌａｔｉｏｎａｌｍｏｔｉｏｎ）、つまり、純粋に対流的な輸送（ｃｏｎｖｅｃｔｉｖｅｔｒａｎｓｐｏｒｔ）をもたらすであろう。そうでない場合は、速度勾配の存在下において、粒子は、破線矢印によって示されるように、回転運動（ｒｏｔａｔｉｏｎａｌｍｏｔｉｏｎ）も行う。全体として、剛体粒子によって実行される最も一般的な種類の運動は、粒子の周囲における「平均速度」による並進（ｔｒａｎｓｌａｔｉｏｎ）と、粒子の近傍におけるその平均値から流速の局所偏差を記述する速度勾配によりもたらされる、質量中心の周りの回転（ｒｏｔａｔｉｏｎ）である。

粘性力が優勢であり、そのため慣性力がわずかであり、単一の繊維が動く範囲における局所的フローの変動が小さいことを仮定することによって、上述の定性的記述は、Ｊｅｆｆｅｒｙ［２］による数学的モデルに組み込まれる。上述のように、熱可塑性プラスチック溶融液中に分散される短繊維の場合、これらの全ての仮定は満たされる。繊維の周囲の流速の

述するのに十分であることを含意している。このテンソルは、３×３行列であり、その要素は、空間における点ｒ＝（χ_１，χ_２，χ_３）^Ｔのデカルト座標χ_ｉに対する流速ベクトルの成分Ｕの偏導

子の中心の空間座標である場合、粒子近傍の点ｒおよび時間ｔにおける流速Ｕは、１次テ

近似される。また速度勾配は、局所的に一定量としてみなされうる。すなわち、ｒ_０近傍の

この場合、空間座標および時間座標の関数としてその配向ベクトルｐ（ｒ，ｔ）によって求められる繊維の過渡的な配向状態（ｔｒａｎｓｉｅｎｔｏｒｉｅｎｔａｔｉｏｎｓｔａｔｅ）は、以下のＪｅｆｆｅｒｙの方程式を使用して計算される。

この方程式は、コンパクトオイラー形式（ｃｏｍｐａｃｔｅｕｌｅｒｉａｎｆｏｒｍ）で本明細書に記載されている。方程式（１）の左辺（ｌ．ｈ．ｓ．）における対流（または材料）微分

は、フローの局所速度場Ｕ（ｒ，ｔ）における純粋な繊維配向、ＦＯ輸送を記述し、右辺（ｒ．ｈ．ｓ．）

無限に細長い繊維（すなわちアスペクト比ｒ_ａ→∞）である理想的な特別の場合では、テン

する繊維の一般的な場合では、それは、繊維の幾何学的形状パラメータλ＝（ｒ_ａ ^２−１）／（ｒ_ａ ^２＋１）を含む以下の項によって定義される。

このパラメータは、粒子幾何学的形状が回転対称楕円形粒子である場合に、それが繊維の回転運動にどのように影響を及ぼすかをコードしている。

Ｊｅｆｆｅｒｙの方程式を記述する他の方法は、速度勾配テンソルのその対称部分および反対称部分への一意的に定義された分解を利用することである。対称部分および反対称部分は、以下によって求められるせん断速度（ｓｈｅａｒｒａｔｅ）および渦度テンソル（ｖｏｒｔｉｃｉｔｙｔｅｎｓｏｒ）によって与えられる。

せん断速度および渦度テンソルは、以下の式を介して有効速度勾配テンソルに関連付けられる。

式の代替形式が以下のように得られる。

Ｊｅｆｆｅｒｙの方程式のさらに別の形式は、渦度テンソル（ｖｏｒｔｉｃｉｔｙｔｅｎｓｏｒ）を、流速ベクトル

おける第１の項を書き直すことによって得ることができる（例えば、［３］参照）。

幾何学的形状が回転対称であるが、楕円形ではない粒子の場合、方程式（１）の形式ならびに有効速度勾配の定義は、同一のままであるが、幾何学的形状パラメータλの公式は、

換えられるように変更しなければならず、その適切な値は、実験により決定しなければならない。この意味では、幾何学的形状パラメータλは、材料パラメータと見なすことができる。非対称非球面粒子の場合、方程式（１）は、３つの幾何学的形状パラメータを含む３つのＪｅｆｆｅｒｙ型方程式の結合系によって置き換えられなければならない［３］。

２．２繊維配向の巨視的分布

マクロなレベルで実行される短繊維強化プラスチックの射出成形シミュレーションの場合、単一の計算セル（すなわち計算領域の体積要素（ＶｏｌｕｍｅＥｌｅｍｅｎｔ））に含まれる繊維の数は、個々の繊維のサイズが小さいため、非常に多い。これを図４に示す顕微鏡写真によって示す。これは、短繊維強化プラスチックから制作された部品から採取される典型的なサンプルである。ゆえに、計算セルに含まれる材料の局所的ＦＯ状態は、巨視的モデル、すなわち、ＦＯベクトルｐの分布関数Ψ（ｐ，ｒ，ｔ）によって記述される（［５］を参照）。繊維濃度は均一であると仮定されるが、分布関数Ψは、依然として局所的な流速場Ｕ（ｒ，ｔ）およびそ

加の点として、その配向に影響を及ぼす繊維の相互作用を説明する項を、マクロレベルでモデルに含めなければならない。

個々の繊維が相互作用しない巨視的モデルとしてのＪｅｆｆｅｒｙの方程式から、相互に作用する多くの繊維のＦＯ統計を与える巨視的モデルへの移行は、以下のＦＯ分布関数Ψを従属変数として有する、対応するＦｏｋｋｅｒ−Ｐｌａｎｃｋ方程式を介して達成される。

演算子を表す。ＦｏｌｇａｒおよびＴｕｃｋｅｒの手法［４］によると、局所速度勾配の有効せん断速度

懸濁液における繊維間相互作用の単純化モデルを得ることができる。式中、Ｇ_ｉｊは、式（１ｂ）

る。無次元で非負である相互作用係数Ｃ_ｌは、懸濁液の材料パラメータである。典型的には、この係数は、短繊維強化プラスチックについて、１０^−３．．．１０^−２の範囲における小さい値（正値）を有する。非圧縮性（ならびにほぼ非圧縮性）のフローでは、「確率的」拡散項の相対的な弱さ（「決定論的な」Ｊｅｆｆｅｒｙ力学を表す項と比較して）は、安定性問題の源であり得ることに留意されたい。

２．３繊維配向テンソルおよびＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式

Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式（２）によって局所的ＦＯ分布を計算することは、フローシミュレーション領域の計算セル毎の単位球面Ｓ^２について定義されるＰＥＤの数値解を必要とし、これは、「工業規模」の３Ｄ問題について極めて高価なタスクである。ゆえに、ＡｄｖａｎｉおよびＴｕｃｋｅｒ［６］は、分布関数のモーメントとして定義される繊維配向テンソルの利用すること、すなわち、Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式を、ＦＯテンソルのモーメント方程式の階層に置き換えることを提案した。方向±ｐが同一の配向状態に対応するという事実を反映する、変数ｐに対するＦＯ分布の反転対称性Ψ（−ｐ，．．．）＝Ψ（ｐ，．．．）により、全ての奇数次モーメントは、一様にゼロになり、結果的に、Ψのモーメント展開は、偶数次の要素のみを含む。ゆえに、この展開の第１の非自明なモーメント（ｎｏｎｔｒｉｖｉａｌｍｏｍｅｎｔ）は、

によって求められる第２のモーメントである。

行列）と呼ばれ、定義により、実対称３×３行列である。ＦＯ分布（これも定義による）は、

である。展開階層（ｅｘｐａｎｓｉｏｎｈｉｅｒａｒｃｈｙ）における次の非自明なモーメントは、

上述のモーメント方程式の階層（ｈｉｅｒａｒｃｈｙ）は、微分演算および積分演算を置き換えることによって、各次数のモーメントについて得られる。すなわち、

（Ｆｏｌｇａｒ−ＴｕｃｋｅｒＥｑｕａｔｉｏｎ，ＦＴＥ）は、階層的に体系化された一連の以下のモーメント方程式において最初の非自明的であるものとして得られる。

Ｊｅｆｆｅｒｙの方程式のように、（５）の左辺における対流微分は、オイラー基準座標系における繊維の並進運動による局所的ＦＯ状態（モーメント展開のこの次数におけるＦＯ行列

の右辺の第３の項は、Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式（２）における拡散項の存在によりもたらされ

減衰効果（ｄａｍｐｉｎｇｅｆｆｅｃｔ）をもたらす。

オイラー形式（５）では、ＦＴＥは、対流−反応型の１次ＰＥＤの結合系である。ラグランジ

２．４ＦＯ行列の基礎特性

（ｕｎｉａｘｉａｌｏｒｉａｔａｔｉｏｎｓｔａｔｅ）を表し、この場合、繊維の１００％は、±ｐ_０により与えられる方向に向いている。したがって、ｐ_０の符号は、この特別なＦＯ行列を定義するダイアディック

であり、残りの２つの固有値は、ｐ_０に直交する平面に存在する対応するベクトルを有するゼロであり、そうでない場合は任意である。

の項は、固有値により与えられる重み付けにより、ＦＯ行列の固有ベクトルによって形成さ

クトルＥ_ｋの方向に沿って配向される繊維の局所分数として固有値μ_ｋを解釈することが可能になる。この意味では、ＦＯ行列のスペクトルデータ｛μ_ｋ，Ｅ_ｋ｝_{ｋ＝１，２，３}は、小さな体積の懸濁液中における繊維の局所的マクロ配向状態を表す。

これは、ＦＯ行列および位相空間に関してＦＴＢを（数式的に）定義することの動機付けとしての役割を果たす［１２］。

のトレースが１に等しい場合に（そしてその場合に限り）、ＦＯ行列である。ＦＴＢの位相空間Ｍ_ＦＴは、全ＦＯ行列の組である。

適切に正規化される分布関数か、そうでない場合は任意の分布関数を使用して、集合Ｍ_ＦＴが、タイプ（３）のモーメント積分からもたらされる全実対称３×３行列の集合に等しいことが分かる。そのトポロジカル特性および幾何学的特性による位相空間Ｍ_ＦＴの数学的特徴付けは、著者のうちの１人によって最近提供されている（［１２］参照）。ＦＴＥの数値積

間ステップも）中に上述の特別なスペクトル特性を有することが必要とされるため、行列が集合Ｍ_ＦＴに属するか否かを厳密に把握することが、実際的に非常に重要である。

２．５閉包問題

いわゆる閉包問題（ｃｌｏｓｕｒｅｐｒｏｂｌｅｍ）は、モーメント展開の各次数において、モーメン

系を解くことができない。閉包近似をＦＴＥに適用することは、（５）の右辺における厳密

既知の例は、以下の式によって定義されるハイブリッド閉包（ｈｙｂｒｉｄｃｌｏｓｕｒｅ）［７］である。

この式は、いくつかの周知の不利点にもかかわらず、その代数的簡易性および数値的ロバスト性により、認められた選択である［７］。

として定義される。

これは、一軸配向分布の特別な場合において厳密な結果をもたらし、線形閉包は次のように与えられる。

これは、等方性（ｉｓｏｔｒｏｐｉｃ）の場合に厳密である。これらの極端な場合の間の補間重み（ｉｎｔｅｒｐｏｒａｔｉｏｎｗｅｉｇｈｔ）は、以下のスカラー配向因子によって提供される。

てはまる。

ハイブリッド閉包近似（６）と組み合わせられるＦＴＥ（５）の特別な変形を、ＦＴＥ−ｈｙｂと表す。この変形は、非常に実際的な利益を有する。本発明の好適な実施形態は、Ｆｏｌｇａｒ−ＴｕｃｋｅｒモデルのＦＴＥ−ｈｙｂ変形にも基づいている。

２．６微分代数系としてのＦＴＥ

（５）の右辺は、任意の実対称行列について数式的によく定義されている。（５）の右辺の特性が解軌道を領域Ｍ_ＦＴに自動的に制限するか否かに関する質問は、少なくともその厳密な形式（すなわち、閉包近似を含まない）において、ＦＴＥについては肯定的に回答でき

Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式（２）の解Ψを使用してモーメント積分（３）を評価することによって得られるＦＯ分布の厳密な２次モーメントに等しい。この議論によって、厳密な４次ＦＯ

軌道が必然的に領域Ｍ_ＦＴに制限されることが結論付けられる。

しかしながら、この議論は、閉包近似が適用される場合は有効ではないし、完全分布関数に関する先験的な知識を持たずに数値的にＦＴＥを解かなければならない場合には、閉包近似が常に必要である。従って、ＦＴＥの解の位相空間Ｍ_ＦＴへの制限に関する問題は、実用性を有する全ての問題において常に発生する。

れる。それは、その不変量の代数不等式によって定式化されうる（［１２］参照）。従属変数に対するこのような制約は、ＦＴＥを微分代数系（ＤＡＳ）に変更し、数値積分に使用される手順において留意しなければならない。

３．ＦＴＥの位相空間の数学的特徴付け

位相空間Ｍ_ＦＴの基礎的なトポロジカル特性は、以下に概説されるように、ＦＯ行列の特別な特性から直接推測できる。

限される実対称３×３行列のベクトル空間の有界凸部分集合（ｂｏｕｄｅｄｃｏｎｖｅｘｓｕｂｓｅｔ）である。

ＦＴＥに適切な積分方法を構築するために、任意の適切なＯＤＥ積分と組み合わせられうるＭ_ＦＴへの射影マッピングの定義が、Ｍ_ＦＴの凸性（ｃｏｎｖｅｘｉｔｙ）によって可能になる（［９］の章ＩＶ．４参照）。ＦＯ行列の不変代数的特徴（ｉｎｖａｒｉａｎｔａｌｇｅｂｒａｉｃｃｈａｒａｃｔｅｒｉｚａｔｉｏｎ）は、次の実対称３×３行列の特性多項式：

の分析によって得られうる。

特徴は、以下の定理に従ってこれらの不変量の項において定式化され得る。

およびＤ_ａが非負である場合およびその場合に限り、ＦＯ行列である。

図６は、特性多項式Ｐ_ａ（μ）の２次項および３次項ならびに線形項および定数項の別々の分析によってこの定理を説明しており、正のトレースＳ_ａは、明らかに正の固有値μ_ｋの存在をもたらし、一方２つの不変Ｋ_ａおよびＤ_ａの非負値は、負の固有値μ_ｋの存在を阻止する。

ることによって、これは、上記定理の証明を完成させる。

定理１に従い、５Ｄ物体である位相空間Ｍ_ＦＴの幾何学的形状の概念は、ＦＯ行列の対角要素および非対角要素を別々に考察することによって得られる。まず、

角形集合（ｔｒｉａｎｇｌｅｓｅｔ）に制限される（図７参照）。

これは、その対角要素に対するＦＯ行列を特徴付ける必要な条件を提供し、位相空間集合Ｍ_ＦＴの「対角部分」の不変表現をもたらす。

領域Ｍ_ＦＴの「非対角」部分の記述形式は、実対称３×３行列の対角要素および非対角要素のトリプレットについて記法（χ，ｙ，ｚ）および（ｕ，ν，ｗ）を導入することによって得られる。ただし、任意の（しかし固定の）対角トリプレット（χ，ｙ，ｚ）∈Δ_ａを考慮する。固定対角トリプレットに属する全ての「許容可能な」非対角トリプレットの以下の集合を数式的に定義する。

代数的には、集合Ｎ_{（χ，ｙ，ｚ）}は、上記の定理２に関して説明したように、以下の不等式のペアを同時に満たす全ての非対角トリプレット（ｕ，ν，ｗ）の集合として特徴付けられることができる。

「最大限に対称的（ｍａｘｉｍａｌｌｙｓｙｍｍｅｔｒｉｃ）である」場合のｘ＝ｙ＝ｚ＝１／３に対応する特別な例について図８に示し、他の場合（図９参照）は、この４面体形状ボリュームの変形バージョンに対応する。

不変量Ｋ_ａおよびＤ_ａの非負値への制限（９）および（１０）を補足するために、正のトレー

集合ΔａおよびＮ_{（χ，ｙ，ｚ）}を導入することによって得られるファイブレーション（ｆｉｂｒａｔｉｏｎ）

理解するために役立つ。つまり、個々の「繊維」は、許容可能な非対角トリプレットの集合Ｎ_{（χ，ｙ，ｚ）}をその「基点」（χ，ｙ，ｚ）∈Δ_ａにローカルに関連付けるという手順によって、可視化される。そして、３角形集合Δ_ａにおけるこの基点の後の変動は、位相空間全体をカバーする。

この節において概説される結果は、位相空間Ｍ_ＦＴが複雑な数学的オブジェクトであることを明確に示す。定理２に記載される厳密な数学的結果によると、対称３×３行列の実ベクトル空間において発展方程式（ｅｖｏｌｕｔｉｏｎｅｑｕａｔｉｏｎ）として必然的に定義されるＦＴＥの任

不変不等式（９）および（１０）の実現を必然的に含意する。この事実によって、必然的にＦＴＥはＤＡＳに変わる。単位トレース条件は、４次ＦＯテンソルの如何なる閉包近似も存在しない場合には、それが満足されることが厳密なＦＴＥに「組み込まれ」、様々な閉包近似（第４節参照）が行われるより一般的な状況の下であっても、依然として有効である。しかし、不変不等式（９）および（１０）の有効性は、概して、現在知られている如何なる閉包近似の存在下においても維持されない（また、「厳密な」ＦＴＥの場合の（９）および（１０）の有効性は、ＦＴＥ自体の代数構造によるよりも、むしろ２．６節に提供される間接的な理由によって推論され得る）。これらの数学的事実は、学術的な著作に通常示される単純な例において見落とされることが多い。

しかしながら、閉包近似は、ＦＴＥを数値的に解く任意の手順の本質的要素として必然的に含まれる。ゆえに、典型的には工業用途において発生する状況等のより複雑な状況に対処するのに適切なあらゆるシミュレーション手続は、解軌道をその理論的に許容可能な領域に制限する適切な制御手続によって、これに対処しなければならない。したがって、ＦＴＥを微分代数系として適切に取り扱うことは、重要な工業用途に必須である。

４．トレース保存および安定性に関する問題

閉包近似をＦＴＥに適用することは、（５）の右辺における厳密な（しかし未知の）４次

さらなる詳細については、閉包問題に関する上の節を参照されたい）。閉包の選択に依存

要性は、以下の式の組によって求められる、いわゆる直交異方性対称性（ｏｒｔｈｏｔｒｏｐｉｃｓｙｍｍｅｔｒｙ）である。

とを仮定する場合、そのトレースに関して、以下の微分方程式（略称：ＤＥ）を導き出すことできる（［１７］参照）。

ＦＴＥの解に関するトレースの安定性（ｓｔａｂｉｌｉｔｙ）および保存性（ｃｏｎｓｅｖａｔｉｏｎ）を、この方

においてまだ使用されていない。閉包近似がＦＴＥに適用される場合、後者は、次の正規化条件を更に満たし、

トレースのＤＥ（１２）は以下のように単純化される。

拡散パラメータＤ_ｒは、定義により非負であるため、トレース条件により定義される超平面（ｈｙｐｅｒｐｌａｎｅ）が、安定積分多様体（ｓｔａｂｌｅｉｎｔｅｇｒａｌｍａｎｉｆｏｌｄ）であること（Ｄ_ｒ＞０の場合）、

レースがＦＴＥの第１の積分であること（Ｄ_ｒ＝０の場合）が結論付けられ得る（厳密な４次

意されたい）。

上の考察は、ＦＴＥへの閉包近似の適用がトレース条件の有効性にいかに影響を及ぼし得

単純なものである。

ハイブリッド閉包（６）が対称特性（１１）のみ有し、正規化条件（１３）を満たさないという特別の場合は、トレースの安定性を慎重に考察しなければならない。このような閉包近似の場合、トレースのＤＥは、以下の形式を仮定する（［１７］も参照）。

これは次の前因子（ｐｒｅｆａｃｔｏｒ）を含む。

この前因子は、（１３）を満たす全ての閉包について有効である６Ｄ_ｒの単純化された表現から外れたものであり、ハイブリッド閉包が正規化条件（１３）を満たさないという事実を反映する。特定の条件において前因子φ_ａに出現する追加の項が、この項を負にし得るこ

を有するＦＴＥの積分多様体であるが、局所的に不安的であり得ることを意味する。

前因子φ_ａが負になる可能性は、以下の４段階の手順の議論によって推論できる。

という特別な場合では、上記前因子は、単純化された形式である

される局所的配向状態がほぼ一軸である場合（すなわち、その固有値μ_ｋのうちの１つが１

の値もほぼ１である。従って、特定の仮定下におけるφ_ａの正確な値を近似する、単純化

得られる。

値を有する少なくとも１つの負の固有値を必然的に有する。

この事実は、前因子φ_ａの近似式の微分においてこれまでに仮定された種々の状況下に

を使用して、この近似式は、以下の形式で記載され得る。

固有ベクトルＥ_ｋ（これはせん断速度テンソルの最大の負の固有値の固有方向に沿う方向を有する）を有するＦＯ行列で評価される場合に、φ_ａが負になることが分かる。この場合、

が得られる。

上記議論は厳密に定式化できる。つまり、任意の非圧縮性フローおよび相互作用係数の

なるため、ハイブリッド閉包を含むＦＴＥの解のトレースが、位相空間Ｍ_ＦＴの特定の領域において局所的に不安定になるという、数学的な証拠を提供することができる。

したがって、ＦＴＥの数値積分中にトレース安定性問題に留意することが明らかに必要である。

５．ＦＴＥの数値積分：一般的側面

（任意の種類の閉包近似を含む）ＦＴＥの数値積分に適切な方法を選択することは、数式のタイプの数学的分類と、ＦＴＥの特定の代数形式に関する側面とに依存する。閉包ＦＴＥの一般的構造は、以下の段落で詳述するように、従属変数としてのＦＯ行列の要素ａ_ｉｊ ^（２）の対流−反応型（ｃｏｎｖｅｃｔｉｏｎ−ｒｅａｃｔｉｏｎｔｙｐｅ）の双曲型偏微分方程式（ｈｙｐｅｒｂｏｌｉｃｐａｒｔｉａｌｄｉｆｆｅｒｅｎｔｉａｌｅｑｕａｔｉｏｎ；ＰＤＥ）の結合系を構成することを示す。

式（５）に与えられるＦＴＥの形式から開始する。この式は、左辺に物質導関数

て次のような式が得られる。この式は、閉包ＦＴＥと同等である。

ＦＴＥのこの形式は、その数学的構造の大部分を既に示している。つまり、左辺は、式（１６）の従属変数であるＦＯ行列の要素ａ_ｉｊ ^（２）の、分離された純粋な対流輸送を支配する局所的な流速Ｕ（ｒ，ｔ）を含む、輸送型（ｔｒａｎｓｐｏｔｔｙｐｅ）または対流型（ｃｏｎｖｅｃｔｉｏｎｔｙｐｅ）の単純な１次偏微分演算子から成る。式（１６）の右辺における種々の項の代数構造は、右辺が、次のよ

称テンソル関数を、（ＦＯ行列自体と同様に、また数学的整合性による必要に応じて、）構成することを示す。

存性を式（１ａ）を用いて解き、有効せん断速度に関して回転拡散パラメータの定義方程式Ｄ_ｒ

閉包ＦＴＥのこの形式は、一定モデルパラメータλおよびＣ_１に対する右辺の依存性を明示的に示す。しかし、ＦＯ行列に対する明示的依存性と黙示的依存性との間を依然とし区別

に対する黙示的依存性との間も依然として区別している。文献によっては、式（１７ｂ）は、場合により添字記法を使用する成分形式、すなわち

と記載される。

全パラメータ依存性を表示せず、またＦＯ行列または速度勾配に対する成分関数Ｆ_ｉｊ ^（２）（．．．）の明示的依存性および黙示的依存性を区別しない、縮約記法（ｃｏｎｄｅｎｓｅｄｎｏｔａｔｉｏｎ）を使用

れる。添字記法（ｉｎｄｅｘｎｏｔａｔｉｏｎ）から、式（１７ａ／ｂ）において使用される完全なテンソル記法に戻すことによって、これらの式に関する以下の単純化バージョンが得られる。

６．ＦＴＥの数値積分：「演算子分割（ｏｐｅｒａｔｏｒｓｐｌｉｔｔｉｎｇ）」

ついては、多種多様の手法が存在する。特に非線形右辺関数Ｆ（．．．）の場合にロバストで柔

により決定される「演算子」とを、別々に扱うことによって進められる。数学的文献において、この手法は、「分別ステップの方法（ｍｅｔｈｏｄｏｆｆｒａｃｔｉｏｎａｌｓｔｅｐｓ）」、「分割方法（ｓｐｌｉｔｔｉｎｇｍｅｔｈｏｄ）」、または単純に「演算子分割（ｏｐｅｒａｔｏｒｓｐｌｉｔｔｉｎｇ）」として知られている（［２３〜３７］参照。代替方法については、［２８〜３１］を参照）。「演算子分割」においては、

のうちのいずれかをゼロに設定することによって、元の式から得られる２つの式

等しく、第２の式は、一般に結合型および非線形である常微分方程式（ｏｒｄｉｎａｒｙｄｉｆｆｅｒｅｎｔｉａｌｅｑｕａｔｉｏｎ；ＯＤＥ）の系列を構成する。このタイプの手法は好適な実施形態において使用される。式（１７ａ／ｂ）の記法を使用し、分割方法を閉包ＦＴＥに適用することによって、以下の部分方程式が得られる。

これらの方程式は、「演算子分割」方法のフレームワーク内における別々の部分問題として扱われる。式（１８ａ）によりモデル化される物理的効果は、流速に従う流体質量の純粋な対流輸送による、モールドキャビティの充填領域内におけるＦＯ統計の空間再分布である（ＦＯ行列の要素によってコードされる）。一方、式（１８ｂ）は、対流輸送の効果を完全に無視し、局所速度勾配および繊維間の相互作用によってもたらされる繊維の回転運動によるＦＯ分布の局所変化率を単に説明するだけである。完全式（１６）（または抽象的記法におけるその同等物（１７ａ／ｂ））の近似解を得るために、両方の部分問題の（数値）解法をどのように組み合わせるかを示す、いくつかのバリエーションについて、以下に説明する。

６．１「単純演算子分割（ＳｉｍｐｌｅＯｐｅｒａｔｏｒＳｐｌｉｔｔｉｎｇ）」

２つの部分方程式から実際に解かれる方程式の近似解を得るための、最も単純な方式（通常、「単純演算子分割」と示す）を用いる場合、第１の式の中間解を、第２の式の入力（すなわち初期値）として解釈することにより、第１の式に続けて第２の式を解く。これを詳

フローソルバー（ｆｌｏｗｓｏｌｖｅｒ）は、懸濁液のフローをモデル化するソフトウェアである。これは、モールドキャビティからゲートシステム、入口に及ぶ計算領域全体のうち、充填部分Ω^（ｎ）に含まれる空間の離散点ｒ_ｍの周囲に位置する全ての計算セルにおける、離散時間ｔ_ｎにおける質量・運動量・エネルギーの離散化輸送方程式の解によって、流体の状態変数、

計算は、先行する計算ステップによって既に把握している「古い」領域Ω^（ｎ）において定義

これは、「単純演算子分割」の利用可能なバリエーションの１つを構成する以下の３つのステップの手順によって達成される。

る。

図１１は、演算子分割の上の変形例を示す。

従って、好適な実施形態に従う方法において使用される。ただし、ステップサイズによっ

する詳細は、後続の節において説明される。

「単純演算子分割」の別の変形例

ップによって出発し、次いで、

この変形では、対流ステップと回転ステップとが逆の順番で実行される。理論的な時間離散化誤差は、上述の第１の変形例と同様に同一次数を有するが、この代替変形例は、異なる瞬間時からの流速および速度勾配を使用するため、一貫性が低い傾向にある。これは、次のの別の変形によって回避できる

（ｉｉ）回転ステップを出発し、その後、

う。

６．２「対称的演算子分割（ＳｙｍｍｅｔｒｉｃＯｐｅｒａｔｏｒＳｐｌｉｔｔｉｎｇ）」

理論的には、時間離散化誤差に関する高次の精度は、「対称的演算子分割」手法によって得られる。この手法の基本概念は、ステップサイズが普通である部分方程式のうちの１つに関する１つの中間ステップを、ステップサイズが半分である他の方程式のステップの中に、括ってしまうことである。この手法に関して可能な１つの変形は、以下の４ステップによって行われる。

る

理論的に可能な別の「対称的分割」による変形例は、ステップサイズが半分である２つの対流ステップの間に回転ステップを一括する。この変形は、

流前」ステップに進み、次いで、

実行し、

流後」ステップを完了させる。

純分割」の場合と同様に、初期初期化ステップ（拡張ステップ）および回転前回転ステップを置き換える第３の変形がさらに存在し、後者から開始することによって、Ω^（ｎ）におけ

モールド充填に見られる典型的なフロー状態は、せん断フローを頻繁に伴うが、それは繊維の顕著な回転をもたらす。このため「回転演算子」部分は、純粋な対流輸送の効果に比べて、常に重要な（そして多くの場合、最も主要な）役割を果たす。ゆえに、対称的分割変形のうちの第１の変形に関して詳述するような、回転前ステップおよび回転後ステップによる対流ステップの対称的な括りは、「対称的分割」手法［１９］の好適な変形例を構成する。しかしながら、ほとんどの最も実際的な場合において、流体フロー計算により提供される時間ステップのサイズは十分小さいので、６．１節に示される「単純分割」による変形例は、多くの場合、十分正確に機能する。

６．３別の分割変形

右辺の関数Ｆ（．．．）の代数構造に依存して、多種多様の更なる分割手法を実現することが

は全部）を残すかもしれない。

式（１６）の右辺を考察すると、３つの項の和が認識され、第１の項および第３の項が線形であり、中央の項がＦＯ行列に対して非線形依存性を恐らく有する閉包近似を伴う（閉包の選択に依存して、この中央の項は、下位の項の別の和に分割されうる）。上述のように、式（１６）の右辺の関数の項構造は、多種多様の追加の分割スキームを形成するための可能性を開く。式（１８ａ／ｂ）への分割に基づく変形例とは対照的に、これらの更なるの代替法は、相互に分離可能な物理的効果（例えば、「対流輸送」と局所「回転力学」の場合）の明確な区別に対応しない。ゆえに、これらは、単に「人工的数学的分離」と見なされるだけであり、物理的ＦＯ現象の適切なシミュレーションについては有用でない。これらの「人工的」分離変形の可能性は、完全性のためにのみ本明細書において記述され、本明細書においてさらに論じられない。

６．４「対流ステップ」に対する数値解析的アプローチ

本節では、閉包ＦＴＥに対する「演算子分割」手法内における部分問題として解決しなければならない対流式（１８ａ）の積分に使用する、数値解析的方法について説明する。新しい

上法（ｆｉｒｓｔｏｒｄｅｒｕｐｗｉｎｄｓｃｈｅｍｅ）［３２］を使用して組み立てられ、結果として得られる数式系

６．５「初期化ステップ」に対する数値解析的アプローチ

な実施形態に従う方法に適用されるアルゴリズムの一部として、修正されたＶＯＦ方法によって実行される。これは、フロー方程式（すなわち、質量・運動量・エネルギーに関する保存方程式）の数値解析のために行われる。ＦＯ行列の値は、両方の領域で共通するセルに

填されたセル」を含むので、ＦＯ行列の要素の合理的な値の初期割当が必要になる。好適な実施形態によると、これは、先行して計算された質量フローに従って、その近傍セルからの又は近傍セルへの計算セル内における質量輸送に対応する加重平均によって達成される。

粒子濃度が均一であると仮定され、かつＦＯ行列が巨視的なレベルにおける体積平均手順からもたらされるため、セルが、その近傍セルから、その中に輸送された流体質量の正味量に比例して、ＦＯへの寄与を得ることが見込まれる。ＦＴＥの位相空間（すなわち可能なＦＯ行列の組）について上に概説した数学的特徴によれば、後者は、（５次元）Ｔｒ［ａ_ｉｊ ^（２）］＝１の超平面に制限される、実対称３×３行列の６次元ベクトル空間の有界かつ凸の部分集合である。質量輸送に従って重み付けされる平均は、ＦＯ行列の凸の組み合わせ（ｃｏｎｖｅｘｃｏｍｂｉｎａｔｉｏｎ）を構成し、有効なＦＯ行列を常にもたらすことから、初期化の手順は、この手順の重要な特性であるＦＴモデルの位相空間の（理論的に必要とされる）トポロジカルな構造に互換性がある。

６．６「噴フロー（ＦｏｕｎｔａｉｎＦｌｏｗ）」効果の取り扱い

用語の「噴フロー（ＦｏｕｎｔａｉｎＦｌｏｗ）」は、多くの粘性非ニュートン流体の場合、流頭における自由表面付近の巨視的フローパターン全体を特徴付ける。「噴フロー」では、流頭付近の上流の粒子は、中心領域（ｃｏｒｅｒｅｇｉｏｎ）から壁境界（ｗａｌｌｂｏｅｎｄａｒｙ）へと移動する。この効果は、流体の材料特性により、実際に「自律的に」発生し、如何なる追加のモデリングをも必要としない。「噴フロー」は、適切な非ニュートン構成則を含むＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式に基づく３Ｄフロー計算においては、少なくとも突発的な現象ではある。しかしながら、これは、フロー計算が、単純化モデル（例えば、［１］におけるＣｒｏｃｈｅｔの論文の節１１．３において論じられるＨｅｌｅ−Ｓｈａｗ型のモデル）に基づく場合には当てはまらない。好適な実施形態に従う方法によるチャネルフローに関するシミュレーションでは、「噴フロー」パターンは、粒子トレースにおいて明確に認識可能であり、この現象をフローソルバーで取り扱うことができることが示される。

自由表面に垂直であるＦＯ成分がゼロになるようにするために、いくつかの特別な措置が必要である。これは、繊維配向の「噴フロー」について観測される効果との一貫性の要請から必要である。またこれは、新しく充填されたセルの繊維配向の初期化手続において行われる。「必要とされる」垂直ＦＯ成分の消去は、次の理由により、上述の「初期化ステップ」において（意図的には）実行されない。数学的観点から、ＦＴＥモデルは、その右辺に結合される双曲輸送方程式系（ｓｙｓｔｅｍｏｆｈｙｐｅｒｂｏｌｉｃｔｒａｎｓｐｏｒｔｅｑｕａｔｉｏｎｓ）（すなわち、非線形反応項を有する対流−反応型系、上記参照）である。ゆえに、自由表面セルにおけるＦＯ行列の境界条件として、自由表面法線の方向におけるＦＯ行列の成分の消滅を規定することは、不適切である（すなわち、数学的に言うとできない）。このＦＯ成分の消滅は、実際は、フローについて計算された「噴フロー」挙動から自動的に生じるはずである。

しかしながら、時間的および／または空間的離散化誤差、ならびにモールド充填中の自由表面の動きの計算に使用する手続によりもたらされる数値的誤差のため、恐らく、新しく充填されたセルの初期化手順が、「噴フロー」現象により必要とされる自由表面におけるＦＯ挙動と一致するように、いくつかの「訂正」を加えなければならない。全ての新しく充填されたセルについて、繊維配向テンソルのトレースが、理論値１との差が１％以下であることが確認される。差が大き過ぎる場合、繊維配向テンソルは、その固有空間における固有ベクトルの許容可能な３角形に対する直交射影によって訂正される。

６．７シミュレーションの開始におけるＦＯ行列の初期化

（閉包）ＦＴＥの双曲線構造は、モールド充填シミュレーションの時間ステップ毎に、境界条件として、入口における（ならびに入口付近の）計算セルにおけるＦＯ行列の初期化を要求する。ゆえに、この初期ＦＯ状態について適切な選択を行なわなければならない。好適な実施形態に従う方法では、後述のように、この目的のために、（ランダムＦＯ分布

高度に粘性のせん断フローでは、ＦＯ状態は、局所的準平衡のその「最終」状態に迅速に押し進められるため、高せん断レートの影響が強い。ゆえに、インゲート（溶融液が実際にその部分に入る場所）において観測されるＦＯ状態は、ランナーシステムにおけるそのフロー履歴によって完全に決定され、入口において見込まれるその初期状態からは相当に独立している。一方、ＦＴＥの解析的構造ならびにその運動学（すなわち位相空間における可能な状態）の分析によれば、ランダムな配向状態が速度勾配の全成分の全ての結合をもたらすことから、ランナーシステムの入口付近におけるフロー場によって後で適合されることを考慮すれば、入口における配向がランダムであるとの仮定は、最適な選択であることが示される。

上述の演算子分割プロセスは、図１１を参照して説明される、単純化された形式である。
ステップ１：新しく充填されたセルに関して初期条件および境界条件を設定し、
ステップ２：フロー場に従って繊維配向を移動し、
ステップ３：速度勾配に従って繊維回転を計算する。

７．「回転ステップ」ＯＤＥ系に対する数値解析的アプローチ

閉包ＦＴＥに対する「演算子分割」法の一部として解かなければならない、ＯＤＥの結合系（１８ｂ）を積分するには、数値解析的方法が使用される。この方法は、具体的には、好適な実施形態に従う方法のＦＯモジュールの実装のために開発されている。この方法は、短繊維強化熱可塑性プラスチック材料の３Ｄ射出成形シミュレーションの典型的な工業用途の枠の中において関連する種々の側面を含んでいる。例えば次の通りである。
・ＦＯ行列の６つの独立要素の完全集合を使用するＦＴＥを定式化すること。（対角要素のうちの１つを無視し、ＦＴＥの従属変数を１つ減らすためにトレース条件を使用する「標準手順」とは異なる。）
・トレース条件により定義された超平面（ｈｙｐｅｒｐｌａｎｅ）を、ＦＴＥの安定積分多様体（ｓｔａｂｌｅｉｎｔｅｇｒａｌｍａｎｉｆｏｌｄ）に変換するＦＴＥの右辺へ、「ペナルティ」制御項を追加すること。
・局所せん断速度のサイズ（速度勾配の最大成分）に従って時間積分スキームを選択する積分法を構築するために、ＦＴＥの右辺関数の特定のスケーリング挙動を、速度勾配の要素に関して利用すること。
・区間［０，１］に制限されるスカラー配向因子を含む「安定化ハイブリッド閉包」（ｓｔａｂｉｌｚｅｄｈｙｂｒｉｄｃｌｏｓｕｒｅ）を実装すること。
・ハイブリッド閉包を有するＦＴＥの特別な右辺関数の評価のための効果的な手順を、最小限の演算数で実装すること。
これらの側面について以下に詳述する。

７．１ハイブリッド閉包の安定化

式（６）および（６ａ／ｂ）によって定義されるハイブリッド閉包近似は、特に、速度勾配が複雑な構造を有し（すなわち、単純なせん断／伸縮型ではなく、複雑な３Ｄフローパターンを反映する）、かつ時間ステップのサイズ（フローソルバーにより決定される）が比較的大きい場合、安全性問題に直面する。これらの不安定性の主な原因は、スカラー配向因子

定になり、ＦＯ行列について許容可能に設定される位相空間集合から大きく離れた値に対して指数関数的に発散する。

る。

ＦＯ行列の行列式は、常に区間［０，１／２７］に制限されるため、式（６ｃ）により定義される配

ッド閉包に関する標準定義（６）は、以下のように修正される。

式（６ａ／ｂ）により与えられる線形項および２次閉包項に関する定義は不変であり、式（６ｃ）

値に評価され、そうでない場合は最小値または最大値として０または１をとる場合に、

の区間［０，１］への制限が、考えられる試験例における深刻な不安定性の発達を阻止することが示された。安定化ハイブリッド閉包は、多種多様の３Ｄ射出成形シミュレーションの検査において成功している。

７．２トレース条件を介する変数の数の減少

３×３行列である。その要素ａ_ｉｊ ^（２）＝ａ_ｊｉ ^（２）は、従属変数において相互に６つの集合を先験的に構成するため、ＦＴＥは、６つのＰＤＥの結合系である。ＦＴＥの代数構造は、任意の対称３×３行列を形式的に認めるが、その行列がＦＯ行列として認められるためには、その行列の不変量における制限として定式化される多数の追加の条件を満たさなければならない。

列の対角要素ａ_ｋｋ ^（２）のうちの１つを無視（解消）しうる明らかな可能性をもたらし、変数の数を１つ減少させる。学問的な多くのフロー状況（例えば、単純なせん断フローまたは伸縮フロー）の下における繊維懸濁液フローを研究する大部分の（全部ではない）公表された研究論文は、解消される特定の対角要素について固定的選択を行なうことによって、

が、式（１４）または（１５ａ）によって、式上、全ての合理的な閉包近似のトレース条件に適

式のように）満たされ、ＦＯ行列の残りの集合の５つの成分およびさらにそのＰＤＥの結合系のみを考慮することで十分である。

この手順は、純粋に形式的（ｐｕｒｅｌｙｆｏｒｍａｌ）である。このように解消される対角要素の選択が良い選択であるか否かは、フローの局所速度場Ｕの空間特性およびその勾配

ーでは、速度勾配の特別な形式が直交ｘ_３方向に対して「弱い結合」のみをもたらし、繊維

の安定性に関する問題）が生じるため、満足できない結果に至ることも考えられる。しか

ために、任意の具体的な単純形式を有することを仮定できないことから、複雑なモールドの幾何学的形状内における３Ｄ射出成形シミュレーションでは、状況は大きく異なっている。そして、（上述のような）固定的選択が、適切な選択ではない可能性が高くなる。ゆえに、トレース条件によってＦＯ行列の対角要素のうちの１つを解消することは、工業用懸濁液フローのシミュレーションには不適切であると考えられる。

単位トレース条件により対角要素のうちの１つを解消することに反対する、（数値的安定性に関する）別の議論が、Ｓｈａｍｐｉｎｅによる一連の論文［３４〜３６］に見られる。この論文では、特別の代数的制約（「保存則」）を有するＯＤＥ系について研究されている。［３４〜３６］

辺関数Ｆ（．．．）に代入することによって、数式から濃度のうちの１つを解消する可能性をもたらし、それによって同時に保存則を厳密に満たす。が、一方、Ｓｈａｍｐｉｎｅは、数値積分手順内におけるこの「トリック」を使用することが、単に、保存則を介して代数的に得られた残りの濃度ｃ_ｎ内において（減少した）ＯＤＥ系の数値積分によって計算される濃度

常に厳密に満たされるが、このように計算される数値解は、正確であることが決して保証されない。実際、特に硬い（ｓｔｉｆｆ）ＯＤＥ系では、この「直接的な解消手法」が不注意に使用される場合に、数値解が、真の解から逸脱し得ることが知られている。［３４〜３６］による

非線形制約が考慮される場合であっても、同一の結論をもたらす。

線形な「保存則」である「単位トレース」制約により拡張されるＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式は、［３４〜３６］に提示される数学的方程式の特定の例であるため、Ｓｈａｍｐｉｎｅの数学的分析は、単位トレース条件により対角要素のうちの１つを解消することによって数値的不安定がもたらされる可能性が高いことを示し、ＦＯ行列の対角要素における制御不可能な不正確性がもたらされることを示す。

好適な実施形態に従う方法は、ＦＯの計算に、ＦＯ行列の６つの要素の全てと、ＦＴＥ系の方程式とを使用する。上述の議論のように、これは、単位トレース条件を介して対角要素のうちの１つを解消する「標準的手法」からは逸脱した方法である。

７．３ペナルティ項（ｐｅｎａｌｔｙｔｅｒｍ）による動的なトレース安定化

閉包ＦＴＥ系およびその「回転ステップ」部分が、式（１６）および（１８ｂ）により与えられると考えられ、かつトレース条件が、従属変数の集合からＦＯ行列の対角要素のうちの１つを解消するために使用されない場合、数値的に計算されたＦＯ行列のトレースを、その望ましい値１に可能な限り近づけて維持するためには、いくつかの追加の方法が必要である。

論文［３４〜３６］によれば、（前の７．２節に論じられたように、）このような解消は、いかなる種類の手法であっても避けることが推奨され、代わりに、通常のＯＤＥ積分方法と、所要の保存則または代数的制約を（正確にまたは近似的に）満たす、何らかの射影手続（ｐｒｏｊｅｃｔｉｏｎｐｒｏｃｅｄｕｒｅ）との組み合わせによって、完全な系を取り扱うことが推奨される。個々の時間ステップにおいて、まずは代数的制約を考慮せずに数値的解が計算され、その後、制約方程式によって定義される代数多様体（例えば、線形制約の場合、超平面）の最近接点への射影マッピング（ｐｒｏｊｅｃｔｉｏｎｍａｐｐｉｎｇ）により補正される。

この手法の別の方法は、適切な制御項をモデル方程式の右辺に追加することによって、制約を考慮することを試みており、数値的に計算されたＦＯ行列のトレースを、その所望値１に可能な限り近似に維持することを試みている。制御項の存在によって、更なる追加の測度（ｍｅａｓｕｒｅ）を含まずに、制約または保存則が近似的に満たされる。多くの場合、「ハード」な制御が、ある種の射影を効果的にもたらし、一方、「ソフト」な制御が、所要の制約からの逸脱を小さくすることを可能にすることを示すことができる。

好適な実施形態に従う方法において、このような制御法は、特に、数値的に計算されたＦＯ行列のトレースを、その所要単位値に近似的に維持するために使用される。安定化は、適切な制御項を式（１６）または（１８ｂ）の右辺関数に加えることによって、トレース条件

の安定した定点になる状況に等しい。

合理的な制御項により満たされるべき一般的な必要性は、以下であるように定められている。
・トレース条件が既に満たされている場合、制御項は全く無いものと同様でなければならない。

なければならない。

これらの２つの主要な必要性の他に、合理的な制御項は、速度勾配の再スケーリングに関してＦＴＥの右辺関数のスケーリング挙動に適合するべきであり（この点に関するさらなる詳細については、７．４節のフロー制御時間積分を参照）、また、「ソフト」または「ハード」な制御の調節を可能にする調節パラメータを含むべきである。後者は、「無限にハード」な制御の制限において、実際上射影に類似する。

ハイブリッド閉包近似（標準的かつ安定化バリエーションの両方を含む）という特定の

て、トレースのＤＥは、以下の式を介して得られる。

（標準的な）ハイブリッド閉包に関する特別な場合では、右辺関数のトレースは、

れる変数前因子φ_ａを有するトレースのＤＥの単純な形式（１５ａ）をもたらす。標準的なハ

因子φ_ａを有するトレースからも、完全に同一のＤＥが得られる。

ナルティ項の或る選択肢は、次の通りである。

パラメータαは正であることが必要とされるが、時間に対して必ずしも一定である必要

定化ハイブリッド閉包を有する閉包ＦＴＥの、一般化された次の変形例も考慮しなければならない。

に従って、式（１５ａ／ｂ）を修正するように調節される。

面は、ハイブリッド閉包を有するＦＴＥの修正された（「コントロールされた」）式（２１）の安定積分多様体（ｓｔａｂｌｅｉｎｔｅｇｒａｌｍａｎｉｆｏｌｄ）になる。結果として、式（２１）の任意の数値解のトレースは、制御項（２０）により所要値１に近似的に維持されるが、それは、この項が、

し進められるためである。

制御項（２０）の「強度」は、パラメータαを適切に調節することによって調整できる。つまり、パラメータαの値が小さいと、トレース値を１の方向へ比較的「ソフト」に訂正

影（ｉｎｓｔｒａｎｔｐｒｏｊｅｃｔｉｏｎ）の効果に近い）。トレース訂正が大きい場合、超平面の両側に数値解を交互に押し出すことによるスプリアス振動（ｓｐｕｒｉｏｕｓｏｓｃｉｌｌａｔｉｏｎ）が生じるため、実際は、αの値が大きいと、明示的ＯＤＥ積分方法と共に使用される際に問題が生じる場合がある。しかしながら、中間的な値においては、比較的広範囲のαの選択が、制御項によりもたらされる不要なアーチファクトを含まない、適切な数値的結果をもたらすことが示されている［１８］。

制御項の以下の一般的定義が、トレースＤＥ（２２）を一意（ｉｄｅｎｔｉｃａｌｌｙ）にもたらし、（安

たらされる。

は、制御項の最も一般的な形式を構成する。

７．４フローが制御された時間積分（Ｆｌｏｗｃｏｎｔｒｏｌｌｅｄｔｉｍｅｉｎｔｅｇｒａｔｉｏｎ）

「演算子分割」手法を使用することによって、「回転ステップ」ＯＤＥ系（１８ｂ）の数値

て（１８ｂ）の右辺に入る。これらが大きい場合、ＦＯ行列要素の値が、時間間隔の間に大きく変動する。一方、速度勾配が非常に小さい場合、先行する計算ステップにより提供される、ＦＯ行列要素の初期値からの変動は、ごく僅かであるだろうということが予測できる。

射出成形シミュレーションに関する全ての実際的な場合では、モールドキャビティ内の充填領域における速度勾配は、充填中に大きく変動する。せん断速度の値は、典型的にはキャビティ壁付近で高く、壁に挟まれる間隙の中心領域においては、かなり低い値が観測されている。これにより、流速方向において、壁付近の強いせん断フローの影響のため、繊維の配向状態はかなり迅速に変化するが、中心領域の繊維は、その配向状態を変化させるのに著しく長い時間を要する。このため、周知の「サンドイッチ構造」（壁付近の高度に配向した領域の間に、配向の度合いが比較的低い「コア」中心領域を有する）がもたらされる。これは、典型的には、短繊維強化熱可塑性プラスチックから作製される射出成形部分において観測される。

ＯＤＥ系の積分の数値スキームは、ステップサイズを適切に（好ましくは適応的に）選択することによって、右辺関数の「強度（ｓｔｒｅｎｇｔｈ）」を勘案する。つまり、値の間隔が大きい粗いステップを用いることや、右辺関数を大きく変化させること、また、右辺関数の（絶対）値がかなり小さい領域において小さ過ぎるステップを取ることにより生じる、不正確性を回避することが試みられる。ステップサイズや積分スキームの適切な選択によって、ＯＤＥ積分の数値誤差を、所望の限界まで制御することが可能になる。

が、これは、速度勾配の変動による右辺関数の変動にもかかわらず、計算領域の全セルにおいて、同じように提供されねばならない。この目的のために、右辺関数の評価の回数を最小限に抑えつつこの目標を達成する、特別な数値スキームが考案されている。このスキームでは、
・「積分法」（オイラー、中間点、または４次ルンゲ・クッタ（ＲｕｎｇｅＫｕｔｔａ））
・（無次元）局所ステップサイズ
・選択された法則とサイズによる、ステップの局所的な数
を選択する。これらは、局所速度勾配の値によって決定される局所右辺関数の「強度（Ｓｔｒｅｎｇｔｈ）」に適応させられる。

このスキームは、数値的ＯＤＥ積分の「適応ステップサイズ制御」に関する既知の手法とは異なる。それは、速度勾配の成分に関する式（１６）または（１８ｂ）の右辺関数のスケーリング挙動（ｓｃａｌｉｎｇｂｅｈａｖｉｏｕｒ）を活用し、ＦＴＥ系に具体的に適応させるためである。積分法ならびにスキームに伴う局所時間ステップの数およびサイズの選択に関する手続は、使用する積分法の理論的離散化誤差に基づく発見的手法である。スキームの両側面について、以下に詳述する。

として、式（１６）から得られてもよい。

よびせん断速度テンソルに関する定義（１ａ／ｂ）、すなわち、

の右辺に関して既に実行されているように、）右辺関数の代数的表現を、以下の単純化形式で記述することが役立つ。

る依存性のみは残す。右辺関数（２３）を使用するＯＤＥ系（１８ｂ）の（（１７ｃ）に対応する）同等な再定式化は、以下の形式で記述できる。

らす。

前述のように、一般的にやらなければならないタスクは、時間間隔［ｔ_ｎ，ｔ_ｎ＋ｌ］におけるＯＤＥの結合系（１８ｃ）の数値積分である。この時間間隔において、グローバル数値誤差ε_ＦＯが存在するが、これは充填セル領域の全セルにおいて均一である。このタスクは、速度勾配の

利用することによって、以下のように達成できる。
・最大の絶対値を有する速度勾配の成分の値、つまり、

局所計算セルを識別する空間ベクトルｒ_ｍや、速度勾配が右辺関数に入る時間を表す時間座

された量および式（２４）を使用して、ＯＤＥ（１８ｃ）を、以下のスケーリングされた形式に変換する。

領域内の各セルは、それぞれ異なるパラメータγ_ｍａｘ．を示す。

に亘り、次の無次元ステップサイズで積分されなければならない。

以下の因子により、ほぼ均一の「単位強度（ｕｎｉｔｓｔｒｅｎｇｔｈ）」を有する。

に等しいか、１未満である無次元量である。

素の値も、１によって制限される。

値を有する。

述は、以下の正確な数学的形式に定式化される。

において、
ｉ．十分小さい離散化誤差を有する「積分法」、ならびに
ｉｉ．スケーリングされた時間間隔をカバーする適切な数のサブステップ
を選択することによって、達成できる。

好適な実施形態に従う方法において使用される積分スキームは、ワン・ステップ・メソッド（ｏｎｅｓｔｅｐｍｅｔｈｏｄ）に分類される単純な積分法の組を使用する。（この積分法については［３８］の節１６．１および［３９］の７．２．１〜７．２．３を参照。これらは参照により本明細書に組み込まれる。）好適な実施形態に従う方法において、ＦＯモジュール内で使用される特別の積分法は次の通りである。
・ステップ毎に右辺関数の評価を必要とする単純な前進オイラー法（ｆｏｒｗａｒｄＥｕｌｅｒｒｕｌｅ）。１次の精度を有する。
・ステップ毎に右辺関数に関する２つの評価を必要とする「中間点（ｍｉｄｏｐｏｉｎｔ）」または２次ルンゲ・クッタ（Ｒｕｎｇｅ−Ｋｕｔｔａ；ＲＫ２）法。２次の精度を有する。
・ステップ毎に右辺関数に関する４つの評価を必要とする４次ルンゲ・クッタ（Ｒｕｎｇｅ−Ｋｕｔｔａ；ＲＫ４）法。４次の精度を有する。

積分法に関する以下の特性が、本発明のある実施形態に使用される特別なスキームの構造において使用されるが、これらは次数（ｏｒｄｅｒ）ｐに依存する。
（ａ）次数ｐである方法は、サブステップ毎に右辺関数に関するｐ回の評価を必要とする。
（ｂ）ＯＤＥ（系）を、全サイズΔτについて、ｈ＝Δτ／ＮのＮ個の等距離サブステップで数値積分する場合、積分の最終サブステップにおいて累積される全誤差ε_ｔｏｔは、

（ｃ）この全誤差が、事前に選択された閾値ｆε_ｍａｘよりも小さいことが必要とされる場合、

（ｄ）積分を終えたときに、全誤差をε_ｍａｘよりも小さくするには、積分間隔のサイズが

ブステップ（すなわち、Ｎ＝１）で十分であることを示唆する。

これを考慮しすると、
（ｉ）３つの積分法のうちの１つによる単一ステップ、または、
（ｉｉ）ＲＫ４法を使用する適切なサイズの多数のサブステップ、
を選択することによって、全体誤差がε_ｍａｘ＝ε_ＦＯより小さくなる積分スキームを得ること

された時間間隔で式（２５）の数値積分を行う、「ハイブリッド（ｈｙｂｒｉｄ）」積分スキーム

に関わらず、均一である。

提案された方法は、具体的には、上記論点（ａ）から（ｄ）に提示される関係において推定されるように、所要の誤差境界を達成するために必要な右辺関数評価の数が最小である

マー溶融液は間隙サイズが小さく高粘度であるので、γ_ｍａｘ、の値はかなり大きくなり得る（約

更新された値も、Ｏ（１）のオーダーとなりうるため、所要の誤差境界の合理的な選択は、

よるＲＫ４法は、半分のサイズｈ／２の２つのステップによる「中間点」（ＲＫ２）法や、４分の１のサイズｈ／４の４つのステップによるオイラー法よりも、正確であることを示すことができる。しかしが、数値計算にかかる労力（右辺関数の４つの評価）は、全てのバリエーションについて同一である。

する値ｐ∈｛１，２，４｝について計算される。０＜ε_ＦＯ＜１であるので（実際は、通常ε_ＦＯ≪１である）、誤差境界のサイズは、オーダーパラメータｐの値が増加するにつれて、単調に増加する（常

視できるほど小さな差しか与えない。
これは、典型的には、速度勾配が小さい領域のいずれのセルにおいても発生する。例えば、２つの隣接する壁の間に比較的大きい間隙が存在する場合、中心中心領域における速度勾配は非常に小さい。例えば、ε_ｍｉｍ＝１０^〜６の最小閾値は、典型的なモールド充填用途において、合理的な選択である。

上述の全ての側面を考慮すると、領域の各セル内において実行するためのアルゴリズムは、以下のように定式化され得る。なお、このアルゴリズムは図１２のフローチャートにも示される。

（ａ）そうである場合、ＯＤＥ積分を飛ばし（すなわち何も実行しない）、セルにおけるＦＯ行列の先行値を、新しい（更新された）値として維持し、手順を終了する。
（ｂ）そうでない場合、ステップ１２．２に進む。

手順を終了する。
（ｂ）そうでない場合、ステップ１２．３．に進む。

ステップによって、セルにおけるＦＯ行列の先行値を更新し、手順を終了する。
（ｂ）そうでない場合、ステップ１２．４に進む。

クッタ」（ＲＫ４）ステップによって、セルにおけるＦＯ行列の先行値を更新し、手順を終了する。
（ｂ）そうでない場合、ステップ１２．５に進む。

有するいくつかのＲＫ４ステップを実行することによって計算される。
（ａ）所要のε_ＦＯ精度を達成するための最小のステップ数は、不等式

数の整数部分を入手し、ｍａｘ（．．．，．．．）を使用して２つの数字のうちの大きい方を入手する

（ｃ）Ｎおよびｈが計算されると、ステップサイズｈを有するＲＫ４法のＮ回のステッ

し、手順を終了する。

ッド」積分法を定義する。この積分法は、誤差が全体的に均一でε_ＦＯより小さく、式（１８ｂ）の右辺関数の評価の必要数が最小で済む積分法である。（１８ｂ）のスケーリングされたバージョン（２５）への当該積分法の適用は、右辺関数の評価において、スケーリングされた速

イズｈによる計算を使用することによって、「同じように（ｉｎｐｌａｃｅ）」実行される。

このアルゴリズムの典型的な用途においては、誤差境界パラメータε_ＦＯ＝０．００１およびε_ｍｉｍ＝１０^−６を有する。この場合、領域Ω_ｎ＋ｌにおけるセルの大部分において（すなわち約８０％）、ＦＯ行列の「回転ステップ」の更新は、単一オイラーステップによって実行される。また比

キップする。恐らくモールドキャビティの壁付近に位置するため、高いせん断速度値を有する適当な数のセルだけが、１回または複数回の（典型的には２０以下の）「４次ルンゲ・クッタ」ステップで更新されることが観測される。

７．５「回転ステップ」ＯＤＥ積分における動的トレース安定化

動的なトレース安定化（ｄｙｎａｍｉｃａｌｔｒａｃｅｓｔａｌｉｚａｔｉｏｎ）を達成すべく、右辺関数が、（ハイブリッド閉包の場合の特別な型（２０）または式（２０ａ）によって与えられる一般的な型のいずれかの）制御項を含む場合に、式（２０）又は（２０ａ）の制御パラメータがα＝α_０・γ_ｍａｘ（α_０は典型的なサイズがＯ（１）の無次元制御パラメータ）と定義されるならば、右辺関数のスケーリング挙動は、依然として、本節に説明される「ハイブリッド」積分アルゴリズムに適合する。

これを、制御項をその一般的な形式（２０ａ）で考察することによって示す。任意の閉包近

て、そのトレースを書き直してみる。

また、制御因子α＝α_０・γ_ｍａｘを使用して、制御項（２０ａ）について以下の等化表現を得る。

る。その他は（２０ａ）式と同一である。これは、以下の式によって表現され得る。

この式はまた、その一般的な形式（２０ａ）の制御項が、制御パラメータがα＝α_０・γ_ｍａｘとして選択される場合に、γ_ｍａｘによる速度勾配のスケーリングに対して、「制御されていな

（安定化）ハイブリッド閉包の特別な場合では、項（２０）を以下のように書き直すこと

た制御項の特別な形式についても（２８）が有効であることを明示的に示す。

上の考察は、ａ_０が制御パラメータとして使用され、右辺関数に依存する項がスケーリングされた速度勾配で評価される場合に、「回転ステップ」ＯＤＥ系の「ハイブリッド」積分アルゴリズムが、制御項の存在下においても変更することなく適用できることを示す。式（２１）に対応するＯＤＥのスケーリングされたバージョンは、単に、制御項を式（２５）の右辺に追加することによって得られる。すなわち、以下のＯＤＥを、積分アルゴリズムの離散時間ステップの基礎として、（２５）の代わりに考慮しなければならない。

好適な実施形態に従う方法のＦＯモジュールに実装される時間積分スキームは、本節に定義される「ハイブリッド」積分アルゴリズムを、制御項（２０）によって、動的トレース安定化を含むＯＤＥ系（２９）に適用することによって、「回転ステップ」ＯＤＥ系の数値積分を実行する。

７．６安定化ハイブリッド閉包を含むＦＴＥの効率的な評価

右辺関数の評価（ｅｖａｌｕａｔｉｏｎ）は、ＦＯ計算手続において最もコストのかかる部分であるため、本方法は、最も効率的なやり方で評価する。

７．４節において提案されるように、「フロー制御時間積分」（ｆｌｏｗｃｏｎｔｒｏｌｌｅｄｔｉｍｅｉｎｔｅｇｒａｔｉｏｎ；ＦＣＴＩ）アルゴリズムによるアプローチによって、効率に関してある程度の対処がなされている。このアプローチを用いることにより、ＦＴＥの右辺関数（恐らく制御項を含む）の評価を最小限に止めつつ、均一な精度を達成可能である。

ＦＯ計算の計算コストに直接影響を及ぼす最重要因子は、閉包近似の選択である。安定化ハイブリッド閉包は、低い計算コストで合理的な精度をもたらすので、本方法で使用される。

ハイブリッド閉包を有するＦＴＥの代数的再定式化

効率改善における第１のステップは、ハイブリッド閉包の代数的定義（６）を、ＦＴＥの右辺関数（２３）に代入することである。これは、（数回の代数的変換後、）次の特有の代数形式を有する以下の右辺関数をもたらす。

子φ_ａに関して上記に提供される式は、（１５ｂ）に示される式に等しい。

公式（３２）は、圧縮性流速場または非圧縮性流速場のいずれの場合においても有効である。ここで、フローが非圧縮性であると仮定すると、理論的にｄｉｖＵ＝０が適用できるは

において維持されることが強調される。数値誤差（これはフロー計算中に必然的に発生す

に比例する項が、前因子公式（３２）においてゼロであると先験的に仮定される場合は、不安定性が発生することが分かる。

量（３１）および（３２）を有する形式（３０）における右辺関数の使用は、効率的な評価に対する基礎的なステップである。というのも、それによって、右辺関数（３０）の代数構造が、多種多様の計算が一度だけ（前因子の計算時に）実行されなければならないように体系化されるため、多く演算を行わなくて済むからである。例えば（３１）および（３２）における項２Ｄ_ｒまたは１−ｆ_ｓ）／７など、頻繁に出現する共通の式があり、それらを一度だけ計算すること、および後に使用するためにその値を補助変数に格納することによって、計算作業が大幅に低減される。

さらに、

ことによって、いくつかの演算を節約することができる。

代数式

に従って再定義される。

関数の評価のために式（３０ａ）および（３１ａ）を使用すると、（３０）および（３１）を使用する場合に比べて、演算の総回数を少なくすることができる。

以下の制御項の追加（式（２０）参照）

は、右辺の評価に必要なさらなる演算をもたらす。この制御項の一般構造は、前因子

の場合は、以下の項

をもたらす。

右辺評価のための効率的なアルゴリズム

上記ステップを利用して、以下の一連の計算に従って、式（３０ｂ）の評価を最小回数の代数的演算で行うスキームが実現される。

１．入力：

算する。

ｖｅｒｓｉｏｎ）

を計算する。

一連の演算によって実行する。

引き続き、右辺関数（３０ｂ）の評価を完成させる。

「縮約記法」（ｃｏｎｔｒａｃｔｅｄｎｏｔａｔｉｏｎ；ＣＮ）を使用する、ベクトル形式におけるＦＴＥの再定式化

く。ゆえに、一般３×３行列としてそれらを格納することは非実際的であり、また、行列

発明によると、対称３×３行列は６つの成分実ベクトルとして処理される。これは、ベク

ば［２２］参照）（構造力学では、これは「フォークト記法」として既知である）として知られる対称行列の指標対（ｉｊ）＝（ｊｉ）を使用して行われる。

つだけである。上記表は、好適な実施形態のＦＯモジュールにおいても導入される通常の選択を示す（代替例は、通常、非対角行列要素のベクトル指標μ_∈｛４，５，６｝への割り当てを選択する）。

れる。

成分は、第１および第２の指標対に対してテンソルが対称的である場合（すなわち、

独立要素の数は、テンソル成分間のさらなる対称性および代数的関係の存在によって低減され得る。完全対称４次テンソルでは、以下の６つの追加の式が存在する。

を、さらに小さい値に減少させる、代数的関係の組である（詳細については［２２］を参照）。

ＣＮ法の導入により、以下のように、「ベクトル形式」によるＦＴＥがもたらされる。

（３５ａ）

である。

される。

（３５ｂ）

提案されたアルゴリズムの計算コスト

右辺関数の効率的な評価について提案されたステップ２から１０の手続において実行される計算演算に関する最も効率的な実装は、ＣＮスキームをこの手順に適用することによってもたらされる。以下の表は、ＣＮ法を使用する手続の個々のステップの実行に必要な数

概説を提供する。

以下、上記に提示されるアルゴリズムの種々のステップに関する表において提供される実装の側面ならびに演算数に関してコメントする。

最終ステップ１０において６ベクトルとしても現れる。

・格納演算および割り当て演算は、計算コスト全体に対するそれらの寄与が無視できるため、勘定されていない。

・補助変数の数は、その値が必要とされなくると、既存の変数を再利用することによって減少される。これは、上に提案されるアルゴリズムのバージョンにおいては説明されていないが、それは、それを説明すると、アルゴリズムの説明の明確性が悪くなり得るためである。

・対称行列の対の縮約演算（ステップ４および６参照）は、ＣＮを使用すると、乗除算７回および加減算５回の演算で実装される（行列記法では、それは、

・実対称３×３行列の行列式（→ステップ５）は、ＣＮを使用して、乗除算１１回および加減算５回の演算で計算される。

およびξ_６を計算するために、２回の乗算を追加する。

１２回の演算を必要とする。

のβの３回の加算のみを必要とする。

表の最終列は、右辺関数の１回の評価が、９６回の乗除算、８４回の加減算、３回の関数呼び出しという計算コストによって得られることを示している。（この関数呼び出しは、（倍精度）浮動小数点数の平方根の計算と、このような数の対の最小値および最大値の計算するために行われる。）現代のコンピュータハードウェアでは、加算演算および乗算演算の計算コストはほぼ同一であり、ＭＩＮ関数またはＭＡＸ関数に対する呼び出しは、約１．５〜２個の演算のコストがかかり、平方根計算の総コストは、約６〜１０個の演算になる（すなわち、所要の精度に応じて平均８個の演算）。

従って、提案されるアルゴリズムによれば、式（３０ｂ）の１回の評価は、全体で、約１９０個の演算の計算コストによって達成される。これは、トレース安定化を含む（安定化）ハイブリッド閉包を有するＦＴＥの右辺を評価するための手法としては、最も効率的である。

計算効率の評価

提案された手順の計算効率に関する評価は、式（３５）によりベクトル形式で与えられる、右辺関数を評価する「標準的な」手法との比較により行うことができる。このような手法は、例えば、コードが特定のモジュール構造を有するように意図される場合に用いなければならないだろう。この特定のモジュール構造とは、（ｉ）閉包近似とは独立の全ての演算

固有の選択に従って、別々のルーティンにおいて計算される、というものである。このような「標準的な」手法の計算コストは、以下の考察によって評価できる。

ラメータ２Ｄ_ｒ（ステップ２から４参照）の計算から出発する。これらは合計４６回の演算と、追加の平方根計算が必要である。

行されると共に、６つの成分の結果として得られる式のみが数値的に評価されるる場合、４８回の演算を必要とする。

算を必要とする。
・２Ｄ_ｒ・１の加算および６Ｄ_ｒ・ａの減算（６Ｄ_ｒ＝３／２・２Ｄ_ｒを使用する）には、別に３＋１２＋１＝１６個の演算を必要とする。

平方根計算について８回の演算を勘定すると、この手法が式（３５）の右辺を評価するための計算コストは、ここまでで既に最大１９０回の演算となる。トレース安定化のために制御項を含めることによって、いくつかのさらなる演算を行わねばならない（ステップ８参

３Ｄ_ｒ＝３／２・２Ｄ_ｒを減算する演算に置き換えることによって省かれ得る。

に必要とされる閉包近似とは別である。

計算に使用されるべき以下の組の式がもたらされる。

４次テンソル式の成分を含む。

たがって決定される。

ステップ５によって実行され、２３回の演算を必要とする（ネストされたｍｉｎ（０，ｍａｘ（．．．，１）

独立行列要素の各々について２回の乗算を必要とするか、合計で４２回の演算を必要とする（残りの行列要素は、対称性により割り当てられる）。

列要素に対するζ_１または３ζ_１の多数の割り当てによって計算される。同様に、更新演算

総演算数は、１０（すなわち、９減算と、３ζ_１の計算のための１演算）に減少する。

よって計算しなければならない（残りは、割り当てにより得られる）。１２個の式のうちの

回である。その後、残りの行列要素は、割り当てによって得られる。

めの割り当てのために２１回の演算を必要とする。

要とする、上述のステップを順に実行する（この手順も非常に効率的に設定されているこ

１１８＋１８＝１３６である。

上記に概説した「標準的手法」は、ハイブリッド閉包への適用時に、全体で約３２６回の演算を必要とする。これは、提案された効率的な手法の計算コストを、およそ１．７２倍（すなわち、約７０％）上回る。提案された効率的な手法を使用することによって、「標準的」手法の計算コストの約６０％で、安定化ハイブリッド閉包（動的トレース安定化を含む）を有するＦＴＥの右辺関数の評価の結果を得ることができる。

数の特定の代数構造の賢明な再体系化によって主に達成される。最適な効率レベルは、ＣＮスキームを介して方程式を減少させることによって、最終的に達成される。

７．７トレース再スケーリング、不変量監視、および位相空間射影

任意のＦＯ行列は、非負固有値および単位トレースを有することを必要とする（→２．４節）。これは、ＦＴＥの解に対する代数的制約の組があることを示唆し、これによって、ＤＡＳになる（→２．６節）。トレースに関する等式Ｓ_ａ＝Ｔｒ（ａ）＝１とは別に、固有値に関する非負制約は、

化できる。

トレース条件は、動的トレース安定化（ｄｙｎａｍｉｃａｌｔｒａｃｅｓｔａｂｉｌｉｚａｔｉｏｎ；ＤＴＳ）を介して、「能動的」かつ計算的に低コストな方式で説明される。ここで、ＦＴＥの右辺には、適切なペナルティ項（→７．３節）が追加され、それは数値解のトレースを所要値Ｓ_ａ＝１に近似であるよ

レース条件を満たすＦＴＥの数値解を得ることができる。実際的な観点からは、これは、解の正確性にほとんど影響を及ぼさない。というのも、因子１／Ｓ_ａによる乗算により行列要素を再スケールすることが常にできるためである。このトレース再スケーリング演算は、数学的に次のように記述される。

また、トレース再スケーリング演算は、次のようないくつかの好ましい特性を有する。

ように）。
・対応する固有ベクトルは、再スケーリング演算により影響を受けない。

ゆえに、トレース再スケーリングの演算は、計算された数値解に対して質的変化をもたらさないが、行列の全ての本質的特性を保ったまま、解をわずかに訂正する。この観点から、非負固定値μ_ｋ（または同等に、非負第２の不変量および第３の不変量Ｋ_ａおよびＤ_ａ）

ＦＴＥの不変量制御の一般的側面

［１２］の結果（→３節）によると、実対称３×３行列の固有値の非負性は、そのトレースが正であり、第２の不変量および第３の不変量Ｋ_ａおよびＤ_ａが非負である場合に保証される。

レース演算を換算するため、トレースに関する発展方程式を直接的に導くことができる（→４節）。ＦＴＥの右辺の特別な代数構造により、この発展方程式は、「厳密な」場合、ならびにハイブリッドまたは一般的な直交異方性閉包近似（式（１４）および（１５ａ／ｂ）参照）のいずれかの場合に、閉形式（ｃｌｏｓｅｄｆｏｒｍ）をとる。ＦＴＥのこのような特別な代数特性により、特別に（意図的に）トレースの発展方程式に影響を及ぼすＦＴＥ自体の右辺にペナルティ項を追加することができる（→ＤＴＳ）。

しかしながら、トレースとは違って、第２の不変量および第３の不変量Ｋ_ａおよびＤ_ａは、行列要素の非線形関数である。ゆえに、これらの不変量の閉形式発展方程式（ｃｌｏｓｅｄｆｏｒｍｅｖｏｌｕｔｉｏｎｅｑｕａｔｉｏｎ）を得るために、同一の（または類似の）手順を適用することはできない（一般には、不変量の発展方程式の右辺は、それらの不変量の組み合わせだけでなく、全行列要素の関数であり、固有ベクトルにおいてコードされる行列の方向性から独立していない）。同様に、任意の特定の項を、既定の方式で非線形不変量の発展に影響を及ぼすＦＴＥの右辺に追加することは、かなり難しい（不可能でないにしても）。ゆえに、不変量を非負領域に維持するために、ＦＴＥの右辺における適切なペナルティ項による第２の不変量および第３の不変量の「能動的」制御は、利用できない。

代替の手順は、「非制御」（または部分的に制御される）積分手順と、その許容可能領域に数値解をマッピングする射影演算との組み合わせから成る。このような「受動的」手順は、数値積分およびその後の射影の演算が、時間ステップ毎に適用される場合に、「非制御」スキームの精度を保存する有効な積分スキームをもたらす（詳細説明については、射影方法に関する［９］の章ＩＶ．４を参照）。

ＦＴＥに適用されるこの組み合わせ手順は、以下から成る。
１．第５節、６節、および７節（７．６節まで）に説明する方法を使用する積分ステップ。
２．続いて、前の積分ステップにおいて計算された行列が依然としてＦＯ行列であるかを決定するための不変量監視ステップ。
３．監視ステップの結果が負である場合（すなわち、不変量条件に反する場合）に、位相空間の組Ｍ_ＦＴに行列をマッピングする位相空間射影ステップ。

位相空間Ｍ_ＦＴは、凸集合であるため（→第３節における定理１）、実対称３×３行列毎に、「最小距離」（適切な測定基準により測定される）においてＭ_ＦＴに一意的な行列が存在する。これは、上記手順の第３のステップに必要な射影マッピングを定義する。位相空間射影に関するさらなる詳細について後述する。

効率的な不変量監視および固有値計算

上述の手順によると、不変量の監視は、時間ステップ毎および計算領域の（充填）セル毎に実行しなければならない。ゆえに、不変量（第２の不変量および第３の不変量またはそれ自体が不変量である固有値）の計算は、いかなる不必要な計算オーバーヘッドをも回避して、できるだけ効率的に実行されなければならない。

第３節で紹介したように、行列の対角要素および非対角要素の記法（χ，ｙ，ｚ）および（ｕ，ν，ｗ）を使用して、不変量は、以下の公式によって、行列要素の多項式として求められる。

これは、全部で１１回の加算および１３回の乗算で評価できる。トレースＳ_ａは、ＦＴＥの数

とによる不変量の監視は、セル毎に時間ステップ当たり追加の２２回の演算により実行できる。

Ｋ_ａおよびＤ_ａのチェックは、実対称３×３行列がＦＯ行列の集合Ｍ_ＭＴに属するか否かを決定する最も効率的な方式である。

詳細に分析すると、行列の固有値の数値計算は、その特性多項式Ｐ_ａ（μ）の根を計算することによって、最も効率的に実行され得ることが示される。これは、Ｐ_ａ（μ）の係数として不変量を計算するのに必要な２４回の演算に加えて約１００回の演算のコストがかかるため、約５倍高くなる。これは、ＦＯ行列特性を監視するために、不変量の代わりに固有値を使用

性の効率的な監視について最適な方法であることを明確にする。行列の数値的対角化は、固有値を計算する最も高価な方法である。ゆえに、このような手順は、固有ベクトルが利益を有する（例えば、後述の位相空間射影手順の一部として）場合を除いて推奨されない。

位相空間射影

不変量監視によって、ＦＴＥの数値解が、現在の時間ステップ内で、位相空間集合Ｍ_ＦＴの

方に反する場合）、解は、位相空間Ｍ_ＦＴへの射影により訂正されなければならない。

称行列の６次元部分空間に対するこれらの量の制限は直接的である。［１８］によれば、任意

れる（正式な証明に関しては［１８］を参照）。

節における式（７）および図７参照）

ユークリッド距離を有する。

において求められる。

また、射影マッピングの値を実際に計算するアルゴリズムは、以下のステップから成る。

を計算する。

この場合、入力行列を出力に割り当て、終了する。そうでない場合（すなわち、Ｋ_ａ＜０またはＤ_ａ＜０である場合）、以下のステップを実行する。

算する。このステップは、数式として以下のように表現され得る。

と一致することに留意されたい。この意味において、（３９ｂ／ｃ／ｄ）は、既に、射影マッピングの一般的定義をコードしている。（３９ｂ／ｃ／ｄ）の結果を計算する個々のステップは、

合に、定義を数式的に完成させる。

列を取るこのアルゴリズムの実際的な実装は、手順のステップ１においてトレースの追加の評価を省略し得る。さらに、実際的な実装は、種々のアルゴリズムステップの以下の変形および／または仕様を組み込み得る。図１３を参照されたい。

・ステップ１では、入力行列の不変量Ｋ_ａおよびＤ_ａを（３８）によって計算する。

る場合には、入力行列をそのまま出力に割り当てる。

・数値的行列対角化によって、例えば、反復周期ヤコビ回転（［３８］、章１１．１または［３９］、章６．５．２）またはギブンズ／ハウスホルダー低減ステップおよび後続のＱＲ繰り返し（［４２］または［３８］、章１１．２を参照）によって、ステップ３の結果を得る。

・実際には、ステップ４は、入力固有値μ_ｋのうちの少なくとも１つが負である場合に

含む）に位置する。

することは、（３９）に従ってμ^＊を決定するために解かなければならない最小化問題を、３角形Δａの３つの端点（１，０，０），（０，１，０），（０，０，１）によって画定される面Ｓ_ａ＝１に制限する。この「面」最小化問題は、以下の射影された３角形

例えば、（μ_１，μ_２）面において同等に解決できる。次いで、μ^＊の第３座標は、

算するアルゴリズムよりも単純であり、より効率的に実装できる。

８．実装

製品の所定の部分における繊維の向きに関する可能性を記述する分布関数であるプロセスの結果は、コンピュータワークステーション（図示せず）のディスプレイにおいてグラフ形式または数値形式で提示される。これは、その上で計算が実行されるコンピュータまたはワークステーションのディスプレイであってもよく、または、これは、その上でシミュレーションが実行されるコンピュータにネットワークで接続されるコンピュータのディスプレイであってもよい。

モールドまたは製品の設計者は、射出成形プロセスにおける製品の質を改善するために、シミュレーションの結果を使用するだろう。また、シミュレーション結果は、関連する製品の製造コストを低減させるために技術者によっても使用されうる。典型的にはポリマー材料よりもさらに高い強度特徴を有する繊維の向きが、製品の強度、剛性、および安定性に対して決定的な影響を有することから、繊維情報に関する信頼できる情報の利点によって、技術者は、製品の強度、剛性、および安定性の特性を決定する前に、良く理解し、情報を有することが可能になる。繊維の向きは、中に懸濁される繊維を含む塊状ポリマーを射出成形する際に発生し得る整経効果に影響を及ぼす。繊維の向きを把握することによって、影響する整経や他の応力を十分に予測することが可能になり、また設計変更を回避することができる。

また、シミュレーションの結果は、ＣＡＤソフトウェア等の他のアプリケーションに渡すこともできる。したがって、シミュレーションの結果は、設計ソフトウェアとシミュレーションソフトウェアとの間の対話型プロセスにおいて使用できる。

本発明は、多数の利点を有する。異なる実施形態または実装によって、以下の利点のうちの１つ以上がもたらされ得る。これが包括的なリストではなく、本明細書に説明されない他の利点も存在し得ることに留意されたい。本発明の１つの利点は、繊維強化製品における繊維配向の分布が、大幅に低減された計算量で決定できることである。本発明の利点は、繊維強化製品における繊維配向の分布が、高い精度で決定できることである。本発明の別の利点は、繊維強化製品における繊維配向の分布がより高速に決定できることである。

本発明について例証目的で詳細に説明したが、このような詳細は、単にその目的に関するだけであり、本発明の範囲を逸脱することなく当業者によって変形を加えてもよいことを理解されたい。

請求項において使用される用語の「備える」は、他の要素またはステップを除外しない。請求項において使用される単数形の用語は、複数形を除外しない。

参考文献
［１］ＪＳ．Ｇ．Ａｄｖａｎｉ（Ｅｄ．）：ＦｌｏｗａｎｄＲｈｅｏｌｏｇｙｉｎＰｏｌｙｍｅｒＣｏｍｐｏｓｉｔｅｓＭａｎｕｆａｃｔｕｒｉｎｇ，Ｅｌｓｅｖｉｅｒ，Ａｍｓｔｅｒｄａｍ（１９９４）
［２］Ｇ．Ｂ．Ｊｅｆｆｅｒｙ：Ｔｈｅｍｏｔｉｏｎｏｆｅｌｌｉｐｓｏｉｄａｌｐａｒｔｉｃｌｅｓｉｍｍｅｒｓｅｄｉｎａｖｉｓｃｏｕｓｆｌｕｉｄ，Ｐｒｏｃ．Ｒ．Ｓｏｃ．Ａ１０２，１６１−１７９（１９２２）
［３］Ｍ．ＪｕｎｋａｎｄＲ．Ｉｌｌｎｅｒ：ＡｎｅｗｄｅｒｉｖａｔｉｏｎｏｆＪｅｆｆｅｒｙ’ｓｅｑｕａｔｉｏｎ，Ｊ．ＭａｔｈＦｌｕｉｄＭｅｃｈ．８，１−３４（２００６）
［４］Ｆ．ＦｏｌｇａｒａｎｄＣＬ．ＴｕｃｋｅｒＩＩＩ：Ｏｒｉｅｎｔａｔｉｏｎｂｅｈａｖｉｏｕｒｏｆｆｉｂｅｒｓｉｎｃｏｎｃｅｎｔｒａｔｅｄｓｕｓｐｅｎｓｉｏｎｓ，Ｊ．Ｒｅｉｎｆ．Ｐｌａｓｔ．Ｃｏｍｐｏｓ．３，９８−１１９（１９８４）
［５］ＣＬ．ＴｕｃｋｅｒａｎｄＳ．Ｇ．Ａｄｖａｎｉ：Ｐｒｏｃｅｓｓｉｎｇｏｆｓｈｏｒｔ−ｆｉｂｅｒｓｙｓｔｅｍｓ，ｃｈ．６，ｐｐ．１４７ｉｎＳ．Ｇ．Ａｄｖａｎｉ（Ｅｄ．）：ＦｌｏｗａｎｄＲｈｅｏｌｏｇｙｉｎＰｏｌｙｍｅｒＣｏｍｐｏｓｉｔｅｓＭａｎｕｆａｃｔｕｒｉｎｇ（ｓｅｅｒｅｆ．［１］ａｂｏｖｅ）
［６］Ｓ．Ｇ．ＡｄｖａｎｉａｎｄＣＬ．ＴｕｃｋｅｒＩＩＩ：Ｔｈｅｕｓｅｏｆｔｅｎｓｏｒｓｔｏｄｅｓｃｒｉｂｅａｎｄｐｒｅｄｉｃｔｆｉｂｅｒｏｒｉｅｎｔａｔｉｏｎｉｎｓｈｏｒｔｆｉｂｅｒｃｏｍｐｏｓｉｔｅｓ，Ｊ．Ｒｈｅｏｌ．７５１−７８４（１９８７）
［７］Ｓ．Ｇ．ＡｄｖａｎｉａｎｄＣＬ．ＴｕｃｋｅｒＩＩＩ：Ｃｌｏｓｕｒｅａｐｐｒｏｘｉｍａｔｉｏｎｓｆｏｒｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌｓｔｒｕｃｔｕｒｅｔｅｎｓｏｒｓ，Ｊ，Ｒｈｅｏｌ．，３６７−３８６（１９９０）
［８］Ｊ．Ｓ．ＣｉｎｔｒａａｎｄＣＬ．ＴｕｃｋｅｒＩＩＩ：Ｏｒｔｈｏｔｒｏｐｉｃｃｌｏｓｕｒｅａｐｐｒｏｘｉｍａｔｉｏｎｓｆｏｒｆｌｏｗ−ｉｎｄｕｃｅｄｆｉｂｅｒｏｒｉｅｎｔａｔｉｏｎ，Ｊ，Ｒｈｅｏｌ．３９，１０９５−１１２２（１９９５）
［９］Ｅ．Ｈａｉｒｅｒ、Ｃ．Ｌｕｂｉｃｈ、ａｎｄＧ．Ｗａｎｎｅｒ：Ｇｅｏｍｅｔｒｉｃｎｕｍｅｒｉｃａｌｉｎｔｅｇｒａｔｉｏｎ，Ｓｐｒｉｎｇｅｒ，Ｂｅｒｌｉｎ（２００２）
［１０］Ｊ．Ｌｉｎｎ，Ｊ．Ｓｔｅｉｎｂａｃｈ，Ａ．Ｒｅｉｎｈａｒｄｔ：Ｃａｌｃｕｌａｔｉｏｎｏｆｔｈｅ３Ｄｆｉｂｅｒｏｒｉｅｎｔａｔｉｏｎｉｎｔｈｅｓｉｍｕｌａｔｉｏｎｏｆｔｈｅｉｎｊｅｃｔｉｏｎｍｏｌｄｉｎｇｐｒｏｃｅｓｓｆｏｒｓｈｏｒｔ−ｆｉｂｅｒｒｅｉｎｆｏｒｃｅｄｔｈｅｒｍｏｐｌａｓｔｓ，ＥＣＭＩ２０００Ｃｏｎｆｅｒｅｎｃｅ，Ｐａｌｅｒｍｏ（２０００）
［１１］Ｊ．Ｌｉｎｎ：ＥｘｐｌｏｒｉｎｇｔｈｅｐｈａｓｅｓｐａｃｅｏｆｔｈｅＦｏｌｇａｒ−Ｔｕｃｋｅｒｅｑｕａｔｉｏｎ，ＳＩＡＭ−ＥＭＳＣｏｎｆｅｒｅｎｃｅ，Ｂｅｒｌｉｎ（２００１）
［１２］Ｊ．Ｌｉｎｎ：Ｏｎｔｈｅｆｒａｍｅ−ｉｎｖａｒｉａｎｔｄｅｓｃｒｉｐｔｉｏｎｏｆｔｈｅｐｈａｓｅｓｐａｃｅｏｆｔｈｅＦｏｌｇａｒ−Ｔｕｃｋｅｒｅｑｕａｔｉｏｎ，ｐ．３２７−３３２ｉｎＡ．Ｂｕｉｋｉｓ，Ｒ．Ｃｉｅｇｉｓ，Ａ．Ｄ．Ｆｉｔｔ（Ｅｄｓ．）：ＰｒｏｇｒｅｓｓｉｎＩｎｄｕｓｔｒｉａｌＭａｔｈｅｍａｔｉｃｓａｔＥＣＭＩ２００２，Ｓｐｒｉｎｇｅｒ（２００４）
［１３］Ｓ．Ｈｅｓｓ：Ｆｏｋｋｅｒ−Ｐｌａｎｃｋｅｑｕａｔｉｏｎａｐｐｒｏａｃｈｔｏｆｌｏｗａｌｉｇｎｍｅｎｔｉｎｌｉｑｕｉｄｃｒｙｓ−ｔａｌｓ，Ｚ．Ｎａｔｕｒｆｏｒｓｃｈ．３１Ａ，１０３４ｆｆ．（１９７６）
［１４］Ｍ．Ｄｏｉ：Ｍｏｌｅｃｕｌａｒｄｙｎａｍｉｃｓａｎｄｒｈｅｏｌｏｇｉｃａｌｐｒｏｐｅｒｔｉｅｓｏｆｃｏｎｃｅｎｔｒａｔｅｄｓｏｌｕｔｉｏｎｓｏｆｒｏｄｌｉｋｅｐｏｌｙｍｅｒｓｉｎｉｓｏｔｒｏｐｉｃａｎｄｌｉｑｕｉｄｃｒｉｓｔａｌｌｉｎｅｐｈａｓｅｓ，Ｊ．Ｐｏｌｙｍ．ＳｃＬ，Ｐｏｌｙｍ．Ｐｈｙｓ．Ｅｄ．１９，２２９−２４３（１９８１）
［１５］Ｍ．Ｇｒｏｓｓｏ，Ｐ．Ｌ．Ｍａｆｆｅｔｏｎｅ，Ｆ．Ｄｕｐｒｅｔ：Ａｃｌｏｓｕｒｅａｐｐｒｏｘｉｍａｔｉｏｎｆｏｒｎｅｍａｔｉｃｌｉｑｕｉｄｃｒｙｓｔａｌｓｂａｓｅｄｏｎｔｈｅｃａｎｏｎｉｃａｌｄｉｓｔｒｉｂｕｔｉｏｎｓｕｂｓｐａｃｅｔｈｅｏｒｙ，Ｒｈｅｏｌ，Ａｃｔａ３９，３０１−３１０（２０００）
［１６］Ｍ．Ｋｒｏｇｅｒ：Ｓｉｍｐｌｅｍｏｄｅｌｓｆｏｒｃｏｍｐｌｅｘｎｏｎｅｑｕｉｌｉｂｒｉｕｍｆｌｕｉｄｓ，Ｐｈｙｓ，Ｒｅｐ．３９０，４５３−５５１（２００４）
［１７］Ｊ．Ｌｉｎｎ：ＴｈｅＦｏｌｇａｒ−ＴｕｃｋｅｒＭｏｄｅｌａｓａＤｉｆｆｅｒｅｎｔｉａｌＡｌｇｅｂｒａｉｃＳｙｓｔｅｍｆｏｒＦｉｂｅｒＯｒｉｅｎｔａｔｉｏｎＣａｌｃｕｌａｔｉｏｎ，ｐｐ．２１５−２２４ｉｎＹ．Ｗａｎｇ，Ｋ．Ｈｕｔｔｅｒ：ＴｒｅｎｄｓｉｎＡｐ−ｐｌｉｃａｔｉｏｎｓｏｆＭａｔｈｅｍａｔｉｃｓｔｏＭｅｃｈａｎｉｃｓ，ｐｒｏｃｅｅｄｉｎｇｓｏｆｔｈｅＳＴＡＭＭ２００４ｃｏｎｆｅｒｅｎｃｅｉｎＳｅｅｈｅｉｍ（Ｇｅｒｍａｎｙ），Ｓｈａｋｅｒ（２００５）
［１８］Ｕ．Ｓｔｒａｕｔｉｎｓ：ＩｎｖｅｓｔｉｇａｔｉｏｎｏｆｆｉｂｅｒｏｒｉｅｎｔａｔｉｏｎｄｙｎａｍｉｃｓｗｉｔｈｉｎｔｈｅＦｏｌｇａｒ−Ｔｕｃｋｅｒｍｏｄｅｌｗｉｔｈｈｙｂｒｉｄｃｌｏｓｕｒｅ，ｍａｓｔｅｒｔｈｅｓｉｓ，Ｄｅｐｔ．ｏｆＭａｔｈｅｍａｔｉｃｓ，ＵｎｉｖｅｒｉｔｙｏｆＫａｉｓｅｒｓｌａｕｔｅｒｎ（２００４）
［１９］Ｊ．Ｌｉｎｎ：ＤｒｅｉｄｉｍｅｎｓｉｏｎａｌｅＶυｒａｉｉｓｂｅｒｅｃｈｎｕｎｇｄｅｒａｎｉｓｏｔｒｏｐｅｎＭａｔｅｒｉａｌｅｉｇｅｎ−ｓｃｈａｆｔｅｎｉｎｔｈｅｒｍｏｐｌａｓｔｉｓｃｈｅｎＳｐｒｉｔｚｇｕｓｓｅｒｚｅｕｇｎｉｓｓｅｎ（Ａｎｌｓｏ−３Ｄ），Ｐｒｏｊｅｋｔ−ｂｅｒｉｃｈｔｆｕｒｄｉｅＭＡＧＭＡＧｍｂＨ，ＴｅｉｌＨａ（２００２）
［２０ＪＪ．Ｌｉｎｎ：ＤｒｅｉｄｉｍｅｎｓｉｏｎａｌｅＶｏｒａｕｓｂｅｒｅｃｈｎｕｎｇｄｅｒａｎｉｓｏｔｒｏｐｅｎＭａｔｅｒｉａｌｅｉｇｅｎ−ｓｃｈａｆｔｅｎｉｎｔｈｅｒｍｏｐｌａｓｔｉｓｃｈｅｎＳｐｒｉｔｚｇｕｓｓｅｒｚｅｕｇｎｉｓｓｅｎ（Ａｎｌｓｏ−３Ｄ），Ｐｒｏｊｅｋｔ−ｂｅｒｉｃｈｔｆｕｒｄｉｅＭＡＧＭＡＧｍｂＨ，Ｔｅｉｌｌｉｂ（２００２）
［２１］Ｂ．Ｅ．ＶｅｒＷｅｙｓｔ，ＣＬ．Ｔｕｃｋｅｒ，Ｐ．Ｈ．Ｆｏｓｓ，Ｊ．Ｆ．Ｏ’Ｇａｒａ：Ｆｉｂｅｒｏｒｉｅｎｔａｔｉｏｎｉｎ３Ｄｉｎｊｅｃｔｉｏｎｍｏｕｌｄｅｄｆｅａｔｕｒｅｓ：ｐｒｅｄｉｃｔｉｏｎａｎｄｅｘｐｅｒｉｍｅｎｔ，Ｉｎｔｅｒｎａｔ．ＰｏｌｙｍｅｒＰｒｏｃｅｓｓｉｎｇ１４，４０９−４２０（１９９９）；
［２２］Ｂ．Ｅ．ＶｅｒＷｅｙｓｔ：Ｎｕｍｅｒｉｃａｌｐｒｅｄｉｃｔｉｏｎｓｏｆｆｌｏｗ−ｉｎｄｕｃｅｄｆｉｂｅｒｏｒｉｅｎｔａｔｉｏｎｉｎｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌｇｅｏｍｅｔｒｉｅｓ，Ｐｈ．Ｄｔｈｅｓｉｓ，Ｕｎｉｖ．ｏｆＩｌｌｉｎｏｉｓａｔＵｒｂａｎａＣｈａｍ−ｐａｉｇｎ（１９９８）
［２３］Ｇ．Ｉ．Ｍａｒｃｈｕｋ：ＳｐｌｉｔｔｉｎｇａｎｄＡｌｔｅｒｎａｔｉｎｇＤｉｒｅｃｔｉｏｎＭｅｔｈｏｄｓ，ｐｐ．１９７−４６２ｉｎＰ．Ｇ．Ｃｉａｒｅｔ＆Ｊ．ＬＬｉｏｎｓ（Ｅｄｓ．）：ＨａｎｄｂｏｏｋｏｆＮｕｍｅｒｉｃａｌＡｎａｌｙｓｉｓ，ＶｏｌｕｍｅＩ，Ｎｏｒｔｈ−Ｈｏｌｌａｎｄ，Ｅｌｓｅｖｉｅｒ（１９９０）
［２４］Ｋ．Ｗ．Ｍｏｒｔｏｎ：Ｎｕｍｅｒｉｃａｌｓｏｌｕｔｉｏｎｏｆｃｏｎｖｅｃｔｉｏｎ−ｄｉｆｆｕｓｉｏｎｐｒｏｂｌｅｍｓ，Ｃｈａｐｍａｎ＆Ｈａｌｌ，Ｌｏｎｄｏｎ（１９９６）
［２５］Ｒ．Ｊ．ＬｅＶｅｑｕｅ，ＮｕｍｅｒｉｃａｌＭｅｔｈｏｄｓｆｏｒＣｏｎｓｅｒｖａｔｉｏｎＬａｗｓ，Ｂｉｒｋｈａｕｓｅｒ（１９９２）
［２６］Ｇ．Ｓｔｒａｎｇ：”Ｏｎｔｈｅｃｏｎｓｔｒｕｃｔｉｏｎａｎｄｃｏｍｐａｒｉｓｏｎｏｆｄｉｆｆｅｒｅｎｃｅｓｃｈｅｍｅｓ”，ＳＩＡＭＪｏｕｒｎ．Ｎｕｍ．Ａｎａｌ．５，５０６−５１７（１９６８）
［２７］Ｍ．Ｇ．ＣｒａｎｄａｌｌａｎｄＡ．Ｍａｊｄａ：”Ｔｈｅｍｅｔｈｏｄｏｆｆｒａｃｔｉｏｎａｌｓｔｅｐｓｆｏｒｃｏｎｓｅｒｖａｔｉｏｎｌａｗｓ”，Ｍａｔｈ．Ｃｏｍｐ．３４，２８５−３１４（１９８０）
［２８］Ｈ．Ｖ．Ｋｏｊｏｕｈａｒｏｖ，Ｂ．Ｍ．Ｃｈｅｎ：”Ｎｏｎｓｔａｎｄａｒｄｍｅｔｈｏｄｓｆｏｒｔｈｅｃｏｎｖｅｃｔｉｖｅｔｒａｎｓｐｏｒｔｅｑｕａｔｉｏｎｗｉｔｈｎｏｎｌｉｎｅａｒｒｅａｃｔｉｏｎｓ”，Ｎｕｍｅｒ．ＭｅｔｈＰａｒｔｉａｌＤｉｊｆ．Ｅｑ．１４，４６７−４８５（１９９８）；”Ｎｏｎｓｔａｎｄａｒｄｍｅｔｈｏｄｓｆｏｒｔｈｅｃｏｎｖｅｃｔｉｖｅ−ｄｉｓｐｅｒｓｉｖｅｔｒａｎｓｐｏｒｔｅｑｕａｔｉｏｎｗｉｔｈｎｏｎｌｉｎｅａｒｒｅａｃｔｉｏｎｓ”ｉｎＲ．Ｅ．Ｍｉｃｋｅｎｓ（ｅｄ．）：Ａｐｐｌｉｃａ−ｔｉｏｎｓｏｆｎｏｎ−ｓｔａｎｄａｒｄｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅｓｃｈｅｍｅｓ，ｍｉｎｉｓｙｍｐｏｓｉｕｍｏｎｎｏｎ−ｓｔａｎｄａｒｄｆｉｎｉｔｅｄｉｆｆｅｒｅｎｃｅｓｃｈｅｍｅｓ：ｔｈｅｏｒｙａｎｄａｐｐｌｉｃａｔｉｏｎｓ，ＳＩＡＭａｎｎｕａｌｍｅｅｔｉｎｇ，ＡｔｌａｎｔａＧＡ，ＵＳＡ１９９９，ｐｕｂｌ．ｂｙＳｉｎｇａｐｏｒｅ：ＷｏｒｌｄＳｃｉｅｎｔｉｆｉｃ（２０００）
［２９］Ｈ．Ｗａｎｇ、Ｘ．Ｓｈｉ、ａｎｄＲ．Ｅ．Ｅｗｉｎｇ：”ＡｎＥＬＬＥＭｓｃｈｅｍｅｆｏｒｍｕｌｔｉｄｉｍｅｎｓｉｏｎａｌａｄｖｅｃｔｉｏｎ−ｒｅａｃｔｉｏｎｅｑｕａｔｉｏｎｓａｎｄｉｔｓｏｐｔｉｍａｌ−ｏｒｄｅｒｅｒｒｏｒｅｓｔｉｍａｔｅ”，ＳＩＡＭＪ．Ｎｕｍｅｒ．Ａｎａｌ．３８，１８４６−１８８５（２００１）
［３０］ＰＪ．ｖａｎｄｅｒＨｏｕｗｅｎ：”Ｎｏｔｅｏｎｔｈｅｔｉｍｅｉｎｔｅｇｒａｔｉｏｎｏｆ３Ｄａｄｖｅｃｔｉｏｎ−ｒｅａｃｔｉｏｎｅｑｕａｔｉｏｎｓ”，Ｊ．Ｃｏｍｐｕｔ．Ａｐｐｌ．Ｍａｔｈ．１１６，２７５−２７８（２０００）
［３１］Ｗ．Ｈｕｎｓｄｏｒｆｅｒ，Ｊ．Ｇ．Ｖｅｒｗｅｒ：”Ａｎｏｔｅｏｎｓｐｌｉｔｔｉｎｇｅｒｒｏｒｓｆｏｒａｄｖｅｃｔｉｏｎ−ｒｅａｃｔｉｏｎｅｑｕａｔｉｏｎｓ”，Ａｐｐｌ．Ｎｕｍｅｒ．Ｍａｔｈ．１８，１９１−１９９（１９９５）
［３２］Ｓ．Ｖ．Ｐａｔａｎｋａｒ：Ｎｕｍｅｒｉｃａｌｈｅａｔｔｒａｎｓｆｅｒａｎｄｆｌｕｉｄｆｌｏｗ，ＨｅｍｉｓｐｈｅｒｅＰｕｂｌ．Ｃｏｒｐ．（１９８０）
［３３］ＣＡＪ．Ｆｌｅｔｃｈｅｒ：ＣｏｍｐｕｔａｔｉｏｎａｌｔｅｃｈｎｉｑｕｅｓｆｏｒＦｌｕｉｄＤｙｎａｍｉｃｓ，ＶｏｌｕｍｅＩ：ＦｕｎｄａｍｅｎｔａｌａｎｄＧｅｎｅｒａｌＴｅｃｈｎｉｑｕｅｓ（２^ｎｄｅｄｉｔｉｏｎ），Ｓｐｒｉｎｇｅｒ（１９９１）
［３４］Ｌ．Ｆ．Ｓｈａｍｐｉｎｅ：”ＣｏｎｓｅｒｖａｔｉｏｎｌａｗｓａｎｄｔｈｅｎｕｍｅｒｉｃａｌｓｏｌｕｔｉｏｎｏｆＯＤＥｓ”，Ｃｏｍｐ．Ｍａｔｈ．Ａｐｐｌｉｃ．１２Ｂ，１２８７−１２９６（１９８６）
［３５］Ｌ．Ｆ．Ｓｈａｍｐｉｎｅ：”ＬｉｎｅａｒｃｏｎｓｅｒｖａｔｉｏｎｌａｗｓｆｏｒＯＤＥｓ”，Ｃｏｍｐ．Ｍａｔｈ．Ａｐｐｌｉｃ．３５，４５−５３（１９９８）
［３６］Ｌ．Ｆ．Ｓｈａｍｐｉｎｅ：”ＣｏｎｓｅｒｖａｔｉｏｎｌａｗｓａｎｄｔｈｅｎｕｍｅｒｉｃａｌｓｏｌｕｔｉｏｎｏｆＯＤＥｓ，ｐａｒｔＩＩ”，Ｃｏｍｐ．Ｍａｔｈ．Ａｐｐｌｉｃ．３８，６１−７２（１９９９）
［３７ＪＪ．Ｌｉｎｎ：”ＥｎｔｗｉｃｋｌｕｎｇｅｉｎｅｓＳｏｆｔｗａｒｅ−ＭｏｄｕｌｓｚｕｒＢｅｒｅｃｈｎｕｎｇｄｅｒＦａｓｅｒｏｒｉｅｎ−ｔｉｅｒｕｎｇｉｎｄｅｒＳｐｒｉｔｚｇｉｅＢｓｉｍｕｌａｔｉｏｎｍｉｔＳＩＧＭＡＳＯＦＴ”，ｔｅｃｈｎｉｃａｌＲｅｐｏｒｔｆｏｒｔｈｅＭＡＧＭＡＧｍｂＨ（２００１）
［３８］Ｗ．Ｈ．Ｐｒｅｓｓ，Ｓ．Ａ．Ｔｅｕｋｏｌｓｋｙ，Ｗ．Ｔ．Ｖｅｔｔｅｒｌｉｎｇ，Ｂ．Ｐ．Ｆｌａｎｎｅｒｙ：ＮｕｍｅｒｉｃａｌＲｅｃｉｐｐｅｅｓｓｉｉｎｎＦＦｏｏｒｒｔｔｒｒａａｎｎ７７７７：：ＴＴｈｈｅｅＡＡｒｔｏｆＳｃｉｅｎｔｉｆｉｃＣｏｍｐｕｔｉｎｇ（２”Ｅｄｉｔｉｏｎ），ＣａｍｂｒｉｄｇｅＵｎｉｖｅｒｓｉｔｙＰｒｅｓｓ（１９９２）
［３９］Ｊ．Ｓｔｏｅｒ，Ｒ．Ｂｕｌｉｒｓｃｈ：ＩｎｔｒｏｄｕｃｔｉｏｎｔｏＮｕｍｅｒｉｃａｌＡｎａｌｙｓｉｓ（３^ｒｄＥｄｉｔｉｏｎ），Ｓｐｒｉｎｇｅｒ（２００２）
［４０］Ｄ．Ｈ．Ｃｈｕｎｇ，Ｔ．Ｈ．Ｋｗｏｎ：”Ａｎｉｎｖａｒｉａｎｔ−ｂａｓｅｄｏｐｔｉｍａｌｆｉｔｔｉｎｇｃｌｏｓｕｒｅａｐｐｒｏｘｉｍａｔｉｏｎｆｏｒｔｈｅｎｕｍｅｒｉｃａｌｐｒｅｄｉｃｔｉｏｎｏｆｆｌｏｗ−ｉｎｄｕｃｅｄｆｉｂｒｅｏｒｉｅｎｔａｔｉｏｎ”，Ｊ．Ｒｈｅｏｌ，４６，１６９−１９４（２００２）
［４１］Ｆ．Ｄｕｐｒｅｔ，Ｖ．Ｖｅｒｌｅｙｅ，Ｂ．Ｌａｎｇｕｉｌｉｅｒ：”Ｎｕｍｅｒｉｃａｌｐｒｅｄｉｃｔｉｏｎｏｆｍｏｕｌｄｉｎｇｏｆｓｈｏｒｔ−ｆｉｂｒｅｃｏｍｐｏｓｉｔｅｐａｒｔｓ”，Ｐｒｏｃ．１ｓｔＥＳＡＦＯＲＭＣｏｎｆ．，２９１−２９４（１９９８）
［４２］Ｇ．Ｈ．Ｇｏｌｕｂ，Ｈ．Ａ．ｖａｎｄｅｒＶｏｒｓｔ：”ＥｉｇｅｎｖａｌｕｅＣｏｍｐｕｔａｔｉｏｎｉｎｔｈｅ２０ｔｈＣｅｎｔｕｒｙ”，Ｊ．Ｃｏｍｐ．Ａｐｐｌ．Ｍａｔｈ．１２３，３５−６５（２０００）

Claims

シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で充填又は部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算する方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、
a）境界条件を特定することと、
b）初期条件を設定することと、
c）巨視的レベルで流体フロー、熱フロー、および質量輸送を求めるために、前記シミュレーション領域の前記セルの少なくとも一部分に関する質量・運動量・エネルギーの平衡方程式を解くことと、
d）空間および時間の関数として、巨視的レベルで前記非球形粒子配向における変化を決定するために、前記解かれた前記平衡方程式の結果に少なくとも部分的に基づいて、非球形粒子配向動的方程式を解くことと、
を含み、前記非球形粒子配向は、分布関数により統計的に記載され、前記分布関数は、増順（increasing order）に並ぶ複数の非球形粒子配向テンソルにより近似され、また、安定化ハイブリッド閉包近似を用いて、４次の配向テンソルが２次の配向テンソルの関数として計算されることにより近似される、方法。
前記分布関数は、縮約記法（contracted notation; CN）を使用してベクトル形式で再定式化された、反応-対流型（reaction convection）の双曲型偏微分方程式の非線形結合系である、請求項１に記載の方法。
前記反応-対流型の双曲型偏微分方程式の系の右辺の代数構造を、いくつかの代数演算によって、安定化ハイブリッド閉包を達成するために利用する、請求項２に記載の方法。
非球形粒子の配向を決定するために必要な代数演算の数を、前記双曲型偏微分方程式の系の右辺において、共通の部分式を同定することにより最小化する、請求項３に記載の方法。
シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で充填又は部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算する方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、
a）境界条件を特定することと、
b）初期条件を設定することと、
c）巨視的レベルで流体フロー、熱フロー、および質量輸送を求めるために、前記シミュレーション領域の前記セルの少なくとも一部分に関する質量・運動量・エネルギーの平衡方程式を解くことと、
d）空間および時間の関数として、巨視的レベルで前記非球形粒子配向における変化を決定するために、前記解かれた前記平衡方程式の結果に少なくとも部分的に基づいて、非球形粒子配向動的方程式を解くことと、
を含み、前記非球形粒子配向は、分布関数により統計的に記載され、前記分布関数は、反応-対流型（reaction convection）の双曲型偏微分方程式の非線形結合系として、増順（increasing order）に並ぶ複数のテンソルにより近似され、前記双曲型偏微分方程式の非線形結合系は、純粋な配向輸送成分と回転配向動力学的要素とに分けられる、方法。
演算子分割（operator splitting）の手法に従い、２つの式を入れ替わりに使用して数値積分を実行する、請求項５に記載の方法。
前記演算子分割は、対称的演算子分割法（symmetric operator splitting）によって、または単純演算子分割法（simple operator splitting）によって実行される、請求項６に記載の方法。
シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で充填又は部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算する方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、
a）境界条件を特定することと、
b）初期条件を設定することと、
c）巨視的レベルで流体フロー、熱フロー、および質量輸送を求めるために、前記シミュレーション領域の前記セルの少なくとも一部分に関する質量・運動量・エネルギーの平衡方程式を解くことと、
d）空間および時間の関数として、巨視的レベルで前記非球形粒子配向における変化を決定するために、前記解かれた前記平衡方程式の結果に少なくとも部分的に基づいて、非球形粒子配向動的方程式を解くことと、
を含み、前記非球形粒子配向は、分布関数により統計的に記載され、前記分布関数は、反応-対流型（reaction convection）の双曲型偏微分方程式の非線形結合系として、増順（increasing order）に並ぶ複数のテンソルにより近似され、前記双曲型偏微分方程式の非線形結合系の解法は、制御項による動的トレース安定化（dynamical trace stabilization）と組み合わせられるフロー制御時間積分（flow controlled time integration）の使用を含む、方法。
前記制御項は、動的トレース安定化のために、前記双曲型偏微分方程式の非線形結合系に含まれる、請求項８に記載の方法。
ペナルティ項の追加によって修正された、前記双曲型偏微分方程式の非線形結合系によって、２次テンソルのトレースが動的に安定化される、請求項８または９に記載の方法。
前記ペナルティ項の関数形式は、前記２次テンソルの前記トレースがその所要値１に近似的に維持されるように特に選択される、請求項１０に記載の方法。
単位トレースを有する全対称行列の集合が、前記双曲型偏微分方程式の非線形結合系の安定積分多様体（stable integral manifold）になるように前記ペナルティ項を選択することによって、前記２次テンソルがその所要値１に近似的に維持される、請求項１１に記載の方法。
シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で充填又は部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算する方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、
a）境界条件を特定することと、
b）初期条件を設定することと、
c）巨視的レベルで流体フロー、熱フロー、および質量輸送を求めるために、前記シミュレーション領域の前記セルの少なくとも一部分に関する質量・運動量・エネルギーの平衡方程式を解くことと、
d）空間および時間の関数として、巨視的レベルで前記非球形粒子配向における変化を決定するために、前記解かれた前記平衡方程式の結果に少なくとも部分的に基づいて、非球形粒子配向動的方程式を解くことと、
を含み、前記非球形粒子配向は、分布関数により統計的に記載され、前記分布関数は、反応-対流型（reaction convection）の双曲型偏微分方程式の非線形結合系として、増順（increasing order）に並ぶ複数のテンソルにより近似され、
さらに、局所的なせん断速度のサイズに従う時間積分のスキームを選択する積分法を構築すべく、前記双曲型偏微分方程式の非線形結合系の右辺関数の特定のスケーリング挙動を、速度勾配の要素に関して利用することを含む、方法。
シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で充填又は部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算する方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、
a）境界条件を特定することと、
b）初期条件を設定することと、
c）巨視的レベルで流体フロー、熱フロー、および質量輸送を求めるために、前記シミュレーション領域の前記セルの少なくとも一部分に関する質量・運動量・エネルギーの平衡方程式を解くことと、
d）空間および時間の関数として、巨視的レベルで前記非球形粒子配向における変化を決定するために、前記解かれた前記平衡方程式の結果に少なくとも部分的に基づいて、非球形粒子配向動的方程式を解くことと、
を含み、前記非球形粒子配向は、分布関数により統計的に記載され、前記分布関数は、反応-対流型（reaction convection）の双曲型偏微分方程式の非線形結合系として、増順（increasing order）に並ぶ複数のテンソルにより近似され、
前記双曲型偏微分方程式の非線形結合系の解法は、FOテンソルまたは行列の提案された解が、前記テンソルまたは行列が非負定値であり、かつ単位トレースを有するという点において許容される解であるどうかかを決定するための監視ステップを含み、許容される解でない場合は、前記提案された解を許容される位相空間に射影される位相空間射影ステップを実行するか、トレース再スケーリングを行うステップを実行する、方法。
請求項１から１４のいずれか１項に記載の方法をコンピュータに実行させる、コンピュータ・プログラム。
鋳造装置と、請求項１５に記載のコンピュータ・プログラムにより制御される、前記鋳造装置に接続されたコンピュータとを備える、鋳造システム。