(57)【要約】多数の非球形粒子を含む懸濁液でモールドキャビティが充填されるプロセスのシミュレーションにおける非球形粒子の統計的な配向分布を記述するための方法および装置。この方法および装置は、繊維強化成形ポリマー部品の分析または繊維強化金属製品を生産するための金属鋳造プロセスの分析に適用され得る。これらの分析の結果を使用して、部品の張力側面および整経側面を決定し、生産プロセスにおいて使用されるプロセス条件を最適化し得る。【選択図】図３
2010531752.tif 000002
【特許請求の範囲】
【請求項１】
シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算するための、コンピュータにより実装される方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、
ａ）境界条件を特定することと、
ｂ）初期条件を設定することと、
ｃ）巨視的レベルで流体フロー、熱フロー、および質量輸送を求めるために、前記シミュレーション領域の前記セルの少なくとも一部分に関する質量・運動量・エネルギーの平衡方程式を解くことと、
ｄ）空間および時間の関数として、巨視的レベルで前記非球形粒子配向における変化を決定するために、前記解かれた前記平衡方程式の結果に少なくとも部分的に基づいて、非球形粒子配向動的方程式を解くことと、
を含む、方法。
【請求項２】
前記ステップｄ）において、前記粒子配向方程式は、懸濁液を含む計算セルについてのみ解かれる、請求項１に記載の方法。
【請求項３】
前記ステップｃ）は、（ｃｃ）前記流体または懸濁液の更新された自由表面または流頭を決定することをさらに含み、自由表面は、前記解かれた平衡方程式の結果に少なくとも部分的に基づいて、前記懸濁液で充填されたセルを、前記キャビティの空セルから分離する、請求項１または２に記載の方法。
【請求項４】
ステップ（ｃｃ）は、前記更新された流頭に従って、前記境界条件を更新することをさらに含む、請求項３に記載の方法。
【請求項５】
ｅ）前記モールドキャビティが、前記懸濁液のシミュレーション射出によって充填されているかを判断することによって、前記シミュレーションされた射出成形プロセスが完了しているかを決定することと、
ｆ）前記シミュレーションされた射出成形プロセスが完了するまでステップｃ），ｃｃ），ｄ）を繰り返すことと、
をさらに含む、請求項３または４に記載の方法。
【請求項６】
多数の非球形粒子を含む溶媒により形成される懸濁液でモールドキャビティが充填されるプロセスのシミュレーションにおいて、非球形粒子の統計的分布配向を記述するための方法であって、
前記キャビティの幾何学的形状を定義する３次元コンピュータモデルを提供することと、
境界条件を特定することと、
複数のセルによるメッシュを形成するために、前記モデルに基づいて解領域を離散化することと、
前記解領域の少なくとも一部分に関するエネルギー方程式およびフロー方程式を解くことと、
各々の前記セルにおけるフロー条件および温度条件を時間の関数として計算することと、
非球形粒子配向の変化を計算することと、
各々の前記セルにおける前記非球形粒子の前記配向の統計的分布を時間の関数として記述することと、
を含む、方法。
【請求項７】
前記流体フローについての方程式を前記粒子配向についての方程式から部分的に分離する原理を使用し、それによって、前記流体フローに関する方程式の計算が、前記粒子配向に関する方程式の計算から独立して実行され、局所的な流速および前記局所的な流速の空間勾配が、前記解かれた平衡方程式の結果から得られ、前記局所的な流速および前記局所的な流速の勾配は、前記粒子配向に関する方程式の入力として使用される、請求項１から６のいずれかに記載の方法。
【請求項８】
前記懸濁液のフローは、流速Ｕ（ｒ，ｔ）およびその空間勾配を含む前記懸濁液の状態変数を与えるシミュレーションモデルによってモデル化され、前記モデル化は、前記解領域の前記充填部分に含まれる離散点ｒｍの周りの空間に位置する全ての計算セルにおける離散瞬間時間ｔｎにおける、質量・運動量・エネルギーの離散化輸送方程式の解によって行われる、請求項７に記載の方法。
【請求項９】
前記モールドキャビティは射出成形装置の一部であり、前記射出成形装置は、前記懸濁液を前記モールドキャビティに供給するためのゲートシステムおよび入口を有し、前記ゲートシステムおよび前記入口の幾何学的形状も前記シミュレーション領域の一部を形成する、請求項請求項１から８のいずれかに記載の方法。
【請求項１０】
前記非球形粒子は、繊維および／またはプレートレット（ｐｌａｔｅｌｅｔ）である、請求項請求項１から９のいずれかに記載の方法。
【請求項１１】
前記溶媒の粘度が十分大きい場合には、前記非球形粒子の移動に対する慣性の効果を無視する、請求項請求項１から１０のいずれかに記載の方法。
【請求項１２】
前記非球形粒子のサイズを、該粒子の近傍における速度勾配が一定であるとで仮定できる程度に十分小さいものと仮定する、請求項１１に記載の方法。
【請求項１３】
前記懸濁液中の前記粒子または繊維配向（ｆｉｂｒｅ  ｏｒｉｅｎｔａｔｉｏｎ；ＦＯ）は、分布関数によって統計的に巨視的なレベルで記述される、請求項請求項１から１２のいずれかに記載の方法。
【請求項１４】
前記分布関数は、増順（ｉｎｃｒｅａｓｉｎｇ  ｏｒｄｅｒ）に並ぶ、粒子又は繊維の複数の配向テンソル（ＦＯテンソル）によって近似され、２次および４次の粒子又は繊維配向テンソルが、前記分布関数の近似に使用される、請求項１３に記載の方法。
【請求項１５】
前記２次および４次粒子又は繊維配向テンソルのみが、前記分布関数の近似に使用される、請求項１４に記載の方法。
【請求項１６】
前記４次配向テンソルに近似解が使用され、前記近似４次配向テンソルは、閉包近似の使用によって、前記２次配向テンソルから計算される、請求項１４または１５に記載の方法。
【請求項１７】
前記閉包近似は、ハイブリッド閉包近似である、請求項１６に記載の方法。
【請求項１８】
前記ハイブリッド閉包近似は、数値的安定性を達成するように安定化されている、請求項１７に記載の方法。
【請求項１９】
前記ハイブリッド閉包近似の安定化は、［０−１］の数値領域に制限されるスカラー配向因子ｆｓｃの使用を含む、請求項１８に記載の方法。
【請求項２０】
前記スカラー配向因子の値は、前記ｆｓｃの計算値が０未満の場合に値を０に設定すること、および、前記ｆｓｃの計算値が０を上回る場合に値を１に設定することによって、前記数値領域［０−１］に制限される、請求項１９に記載の方法。
【請求項２１】
空間座標および時間座標の関数としてその配向ベクトルにより求められる粒子または繊維の過渡配向を、単一線位の決定論的配向動力学的に関するＪｅｆｆｅｒｙのモデルに基づいたものとし、さらに、繊維間相互作用を確率的にモデル化するランダム項をさらに含む、請求項１から２０のいずれかに記載の方法。
【請求項２２】
前記繊維間相互作用をモデル化する回転拡散係数を含む、前記粒子または繊維分布関数のＦｏｋｋｅｒ−Ｐｌａｎｃｋ方程式を、前記粒子または繊維の前記過渡配向動力学的をモデル化するために使用する、請求項２１に記載の方法。
【請求項２３】
ＦｏｌｇａｒおよびＴｕｃｋｅｒの手法に従い、また、前記回転拡散係数が、前記局所速度勾配の有効せん断速度に比例するように設定される、請求項２２に記載の方法。
【請求項２４】
前記分布関数の前記繊維配向（ＦＯ）テンソル近似は、Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式と組み合わせて使用されることにより、Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式（Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ  Ｅｑｕａｔｉｏｎ；ＦＴＥ）と呼ばれる増順（ｉｎｃｒｅａｓｉｎｇ  ｏｒｄｅｒ）のＦＯテンソルの相互結合方程式の階層をもたらし、また、解かれる前記粒子配向方程式は、前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式または修正されたＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式に少なくとも部分的に基づく、請求項２２または２３に従い且つ請求項１４から２１のいずれかに１項に従う方法。
【請求項２５】
解かれる前記粒子配向方程式は、いわゆるＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式（ＦＴＥ）または修正されたＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式に少なくとも部分的に基づく、請求項請求項１から２４のいずれかに記載の方法。
【請求項２６】
前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式は、粒子配向および輸送における変化を計算するために利用される、請求項請求項１から２５のいずれかに記載の方法。
【請求項２７】
前記４次ＦＯテンソルの近似解を計算するために、前記ＦＴＥの前記２次ＦＯテンソル方程式に閉包近似を適用する、請求項２４、２５、または２６に記載の方法。
【請求項２８】
前記懸濁液中の前記非球形粒子の前記配向の統計的分布を記述するために２次モーメントを利用する、請求項２４から２７のいずれかに記載の方法。
【請求項２９】
前記非球形粒子の前記配向の統計的分布は、前記配向分布密度関数Ψ（ｐ）の前記２次配
2010531752.tif 000345
に記載の方法。
【請求項３０】
前記４次ＦＯテンソルの近似解を計算するために、前記２次ＦＯテンソル方程式に安定化ハイブリッド閉包近似を適用するか、または前記４次配向テンソルを、安定化ハイブリッド閉包近似を使用して、前記２次テンソルの関数として計算する、請求項２７から２９のいずれかに記載の方法。
【請求項３１】
前記ハイブリッド閉包近似の安定化は、［０−１］の数値領域に制限されるスカラー配向因子の使用を含む、請求項３０に記載の方法。
【請求項３２】
ペナルティ項の追加によって修正された、前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式によって、前記２次
2010531752.tif 000346
項３０または３１に記載の方法。
【請求項３３】
2010531752.tif 000347
所要値１に近似的に維持されるように特に選択される、請求項３２に記載の方法。
【請求項３４】
単位トレースを有する全対称３×３行列の集合が、前記修正されたＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式の安定積分多様体（ｓｔａｂｌｅ  ｉｎｔｅｇｒａｌ  ｍａｎｉｆｏｌｄ）になるように前記ペナルティ項を選択することによって、前記配向テンソルがその所要値１に近似的に維持される、請求項３３に記載の方法。
【請求項３５】
前記ハイブリッド閉包近似は、前記ハイブリッド閉包近似によって前記４次ＦＯテンソ
2010531752.tif 000348
の間のその理論的に許容可能な領域に括ること（ｂｒａｃｋｅｔｉｎｇ）によって安定化される、請求項３０から４０のいずれかに記載の方法。
【請求項３６】
前記スカラー配向因子の値は、前記ｆｓｃの計算値が０未満の場合に値を０に設定すること、および、前記ｆｓｃの計算値が０を上回る場合に値を１に設定することによって、前記数値領域［０−１］に制限される、請求項３１または３５に記載の方法。
【請求項３７】
前記繊維配向の計算のための前記時間積分方法は、前記局所せん断速度レベルに局所的に適合される、請求項２４から３６のいずれかに記載の方法。
【請求項３８】
前記適合は、
（ａ）せん断速度増加による時間ステップサイズの減少、
（ｂ）せん断速度増加による前記時間積分方法の次数の増加、および／または、
（ｃ）せん断速度増加による、いくつかの小さい時間ステップにおける前記時間ステップの局所的分割、
を含む、請求項３７に記載の方法。
【請求項３９】
前記適合（ａ）（ｂ）（ｃ）のうちの少なくとも１つは、前記フロー計算により求められるグローバルな時間ステップにおける誤差が、全ての充填セルについて近似的に同一であるように使用され、前記同一であることは、前記せん断速度の局所変動が非常に大きい場合であっても、該記せん断速度の局所変動とは独立して同一であるようにされる、請求項３８に記載の方法。
【請求項４０】
前記繊維配向行列の第２の不変量および第３の不変量の数値によって、前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式の前記数値解のその可能な領域または位相空間への制限を制御するために、不変量（ｉｎｖａｒｉａｎｔ）利用される、請求項２４から３９のいずれかに記載の方法。
【請求項４１】
前記懸濁液の圧縮性の効果が、該圧縮性の効果が非常に小さい場合であっても、前記モデルに含まれる、請求項請求項１から４０のいずれかに記載の方法。
【請求項４２】
2010531752.tif 000349
−対流型（ｒｅａｃｔｉｏｎ  ｃｏｎｖｅｃｔｉｏｎ）の双曲型ＰＤＥの非線形結合系であり、また、純粋な配向輸
2010531752.tif 000350
けられる、請求項２４から４１のいずれかに記載の方法。
【請求項４３】
演算子分割（ｏｐｅｒａｔｏｒ  ｓｐｌｉｔｔｉｎｇ）の手法に従い、２つの式を入れ替わりに使用して数値積分を実行する、請求項４２に記載の方法。
【請求項４４】
前記演算子分割は、対称的演算子分割法（ｓｙｍｍｅｔｒｉｃ  ｏｐｅｒａｔｏｒ  ｓｐｌｉｔｔｉｎｇ）によって、または単純演算子分割法（ｓｉｍｐｌｅ  ｏｐｅｒａｔｏｒ  ｓｐｌｉｔｔｉｎｇ）によって実行される、請求項４３に記載の方法。
【請求項４５】
前記シミュレーションの開始時の前記非球形粒子の初期配向は、前記非球形粒子が実際のモールド充填プロセスにおいて有する分布にかかわらず、ランダム化分布である、請求項請求項１から４４のいずれかに記載の方法。
【請求項４６】
前記自由表面における粒子の噴フロー効果は、前記繊維配向の初期化において無視される、請求項請求項１から４５のいずれかに記載の方法。
【請求項４７】
自由表面における前記新しく充填されたセルの前記繊維の前記繊維配向分布は、前記充填シミュレーションの時間ステップ毎おける近傍セルからの前記質量輸送に従って、加重平均によって初期化される、請求項請求項１から４６のいずれかに記載の方法。
【請求項４８】
前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式（１６）は、縮約記法（ｃｏｎｔｒａｃｔｅｄ  ｎｏｔａｔｉｏｎ；ＣＮ）を使用してベクトル形式で再定式化される、請求項２４から４８のいずれかに記載の方法。
【請求項４９】
前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式（１６）の右辺の代数構造は、いくつかの代数演算によって、安定化ハイブリッド閉包を達成するために利用される、請求項２４から４８のいずれかに記載の方法。
【請求項５０】
前記繊維配向を決定するために必要な代数演算の数を、前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式（１６）の右辺において共通部分式を同定することによって最小化する、請求項４９に記載の方法。
【請求項５１】
前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式（１６）は、効率的な評価に最適である項構造を達成するために、前記ハイブリッド閉包の定義公式を、分析計算および結果として生じる項の再体系化の右辺関数に代入することによって再定式化される、請求項４９または５０のいずれかに記載の方法。
【請求項５２】
制御項は、動的トレース安定化のために、前記再定式化されたＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式（１６）に含まれる、請求項４８から５１のいずれかに記載の方法。
【請求項５３】
前記再定式化された右辺関数は、最適効率のためにＣＮを使用して最少数の演算で評価される、請求項４８から５２のいずれかに記載の方法。
【請求項５４】
前期ＦＴＥの前記数値解が、前記第２の不変量および第３の不変量の制御によって、許容可能な領域にあるか否かをチェックすることをさらに含む、請求項４８から５３のいずれかに記載の方法。
【請求項５５】
ＤＴＳにより近似的にのみ単位トレースを有するＦＯ行列を解明するために、「トレース再スケーリング」の利用をさらに含む、請求項４８から５４のいずれかに記載の方法。
【請求項５６】
位相空間射影のアルゴリズム記述を含む、請求項４８から５５のいずれかに記載の方法。
【請求項５７】
前記許容可能な領域に前記数値解を制限するために、前記ＦＴＥの数値積分に関する他の提案された技術と、位相空間射影とのの組み合わせをさらに含む、請求項５６のいずれかに記載の方法。
【請求項５８】
前記Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ方程式（ＦＴＥ）または修正されたＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式の解であるＦＯテンソルの解軌道または集合は、粒子もしくは繊維配向テンソルまたは行列の位相空間集合ＭＦＴに制限され、行列の前記位相空間ＭＦＴ集合は、非負定値であり、かつ単位トレースを有する実対称３×３行列のベクトル空間の有界凸部分集合であり、それによって、前記ＦＴＥは、微分代数系（ｄｉｆｆｅｒｅｎｔｉａｌ  ａｌｇｅｂｒａｉｃ  ｓｙｓｔｅｍ；ＤＡＳ）と考えられる、請求項２７から５７のいずれかに記載の方法。
【請求項５９】
前記ＦＴＥの解法は、ＦＯテンソルまたは行列の提案された解が、前記テンソルまたは行列が非負定値であり、かつ単位トレースを有するという点において許容される解であるかどうかを決定するための、監視ステップを含み、許容される解でない場合は、前記提案された解が前記位相空間ＭＦＴに射影される位相空間射影ステップを実行する、請求項５８に記載の方法。
【請求項６０】
前記提案されたＦＯ解行列の前記射影は、前記提案されたＦＯ行列を、前記位相空間ＭＦＴにおける最近接行列にマッピングすることによって実行される、請求項５９に記載の方法。
【請求項６１】
前記監視ステップは、前記ＦＯテンソルまたは行列の前記第２の不変量および第３の不変量の数値を決定することを含む不変量監視ステップである、請求項５８または５９のいずれかに記載の方法。
【請求項６２】
前記位相空間ＭＦΓにおける前記最近接行列への前記提案されたＦＯ行列の前記マッピングは、前記提案されたＦＯ行列の固有値のトレース再スケーリングの使用を含む、請求項６０または６１に記載の方法。
【請求項６３】
前記ＦＴＥを解くことは、初期化ステップ、回転前ステップ、対流ステップ、および回転後ステップを含む対称的一括手法である対称的演算子分割の使用を含む、請求項７および請求項２４から６２のいずれかに記載の方法。
【請求項６４】
前記ＦＴＥの前記方程式を解くことは、前記ＦＯテンソルまたは行列の全６つの要素の使用を含む、請求項２４から６３のいずれかに記載の方法。
【請求項６５】
時間間隔内において前記ＦＴＥの前記方程式を解くことは、スケーリングされた流速勾配およびスケーリングされた時間座標の使用を含む、請求項２４から６４のいずれかに記載の方法。
【請求項６６】
前記ＦＴＥの前記方程式の解法は、制御項による動的トレース安定化（ｄｙｎａｍｉｃａｌ  ｔｒａｃｅ  ｓｔａｂｉｌｉｚａｔｉｏｎ）と組み合わせられるフロー制御時間積分（ｆｌｏｗ  ｃｏｎｔｒｏｌｌｅｄ  ｔｉｍｅ  ｉｎｔｅｇｒａｔｉｏｎ）の使用を含む、請求項３４から６５のいずれかに記載の方法。
【請求項６７】
前記シミュレーションの結果は、シミュレーションされる前記射出成形プロセスによって製造される製品の特徴を改善するために使用される、請求項請求項１から６６のいずれかに記載の方法。
【請求項６８】
前記シミュレーションの結果は、製造される製品の強度、剛性、または安定性特徴を改善するために使用される、請求項６７に記載の方法。
【請求項６９】
前記シミュレーションの結果は、製造される製品の変形の量を低減させるために、形状安定性を改善するために使用される、請求項６７に記載の方法。
【請求項７０】
前記粒子配向の統計的分布によって、前記フロー計算の巨視的レベルにおいて前記局所繊維配向状態を特徴付けることをさらに含む、請求項請求項１から６９のいずれかに記載の方法。
【請求項７１】
請求項１から７０のいずれかに記載の方法を実行するための、ソフトウェアアルゴリズムを含むコンピュータ可読媒体上のコンピュータ製品。
【請求項７２】
鋳造装置と、それに接続される、請求項７１に記載のコンピュータ製品を含むコンピュータにより実装される最適化システムとを備える、最適な鋳造のためのシステム。
【発明の詳細な説明】
【技術分野】
【０００１】  本発明は、粒子を含む懸濁液（ｓｕｓｐｅｎｓｉｏｎ）のモールドキャビティへのフロー（ｆｌｏｗ）に関する３次元モデリングの分野に関し、より具体的には、多数の非球形粒子を含む懸濁液でモールドキャビティが充填されるプロセスに関するシミュレーションにおいて、非球形粒子の統計的な配向分布を記述するための方法および装置に関する。
【背景技術】
【０００２】  射出成形プロセスまたは金属鋳造プロセスに関する忠実な３Ｄシミュレーションは、最大で数十万個もの数式系を伴う。これらの複雑な計算に対処するために、シミュレーション方法の効率の向上が進められてきている。最適化ソフトウェアおよび最新のワークステーションの処理能力によって、このようなシミュレーションを普通の職場で実行することが可能となり、十分に適切な結果を得られるようになった。その結果は、純粋に科学的な研究分野ではともかく、射出成形品の研究開発部や、鋳物工場、製造業者の技術者によっては十分に実用的である。【０００３】  ポリマーの組成に繊維等の粒子が加えられ、その繊維の配向分布を記述する必要がある場合、シミュレーションおよび付随する式の集合は、さらに大幅に複雑になる。このようなプロセスに関する忠実な３Ｄシミュレーションは、これまで、計算時間が長過ぎるか、またはシミュレーションの精度が不十分となる恐れがあるため、シミュレーションの複雑性により、最新のワークステーションにおいても結果が許容できないことから、一般の職場への導入に成功していない。【０００４】  繊維で強化された部分では、開発技術者が、部品の張力（ｔｅｎｓｉｏｎ）や整経（ｗａｒｐｉｎｇ）について予測できるように、繊維配向分布の記述を有することが極めて重要であることが多い。繊維は、通常、プラスチック部分の力学的特性を改善するために使用される。しかしながら、（熱）力学的特性（熱膨張、強度、および剛性等）は、繊維の向きに依存する。【０００５】  射出成形プラスチック部品の使用は、近年、多くの産業において着実に増加している。電子機器、消費財、医療機器、および自動車部品の製造業者は、かつてない程、プラスチックから作製される製品に使用する製品および部品の生産を増やしている。【０００６】  繊維で強化された射出成形部品は、改良された強度／重量比率、耐久性、部品の一体化、および低コストを提供するため、構造用金属部品に取って代わりつつある。【０００７】  同時に、競争の圧力により、プラスチック射出成形産業の製造業者は、設計を生産プロセスにさらに適合させるために、設計を最適化する新しい方法を探索することを余儀なくされている。部品やモールド構成を修正する必要性が、設計開発プロセスの後期に発見される場合、必要な変更を実行するための遅延コストおよび関連コストは、設計開発段階の初期段階よりも大幅に高くなる。【０００８】  自社の部品が生産可能であること、かつ最適に機能することを確実にしたい会社は、製造プロセスを十分に理解し、かつこの知識を設計段階の初期において部品設計に組み込むために、コンピュータ支援エンジニアリング技術を使用して、射出成形モールドにおける複雑なフローや結果として生じる繊維配向をシミュレーションおよびモデル化することを望む。【０００９】  繊維強化部品のための、およびその中で生産される繊維強化部品のための、射出成形モールドの設計には、考慮すべき多くの因子が存在する。部品全体の幾何学的形状、最小および最大壁厚さ、モールドにおいて、液体ポリマーや繊維懸濁液が射出されるゲートの数および位置、モールドにおいて、キャビティ中の気体を排出するベントの数および位��、ポリマーの組成や特性、繊維の特性や量、鋳引け、余裕部分、繊維配向分布等のパラメータがある。部品やモールドの設計は、最終部品の形式および機能に単に大きく依存するわけではなく、それに密接に相関する製造プロセスの影響も考慮すべきである。【００１０】  便利なことに、コンピュータ支援エンジニアアリングシミュレーション（Ｃｏｍｐｕｔｅｒ  Ａｉｄｅｄ  Ｅｎｇｉｎｅｅｒｉｎｇシミュレーション；ＣＡＥシミュレーション）は、射出成形プロセス中にモールドキャビティの内部で発生し得ることに関して、視覚的および数的なフィードバックを、設計および製造の技術者に提供するように用いられることができる。これによって、技術者は、想定される部品設計の挙動を十分に理解および予測することが可能になり、結果的に、製造に対する従来の高価な試行錯誤の手法を、実質的に排除することが可能となる。コンピュータ支援エンジニアアリングシミュレーションの使用によって、設計段階中の部品設計やモールド設計、製造プロセスパラメータの最適化が容易になり、設計段階中に、最小コストおよびスケジュールに対する最小の影響で、必要な変更を容易に実行することが可能になる。【００１１】  繊維強化部品のエンジニアリングプロセス内におけるＣＡＥシミュレーション技術の適用は、（ｉ）流体フローおよび熱伝導の計算を含む「射出成形」製造プロセスのシミュレーション、ならびに（ｉｉ）応力や強度、場合によっては耐久性の計算を含み、これらの全ては、これらの部品について巨視的なレベルで実行され、外部負荷の存在下におけるその機能的な力学的特性を決定する。両方の種類のシミュレーションは、液体状態および固体状態において浸漬繊維を含むポリマー材料の特性を記述する適切な材料モデルを必要とする。【００１２】  巨視的レベルにおける長さのスケールは、部品の幾何学的形状の線形サイズ（全体のサイズ、壁厚さ等）によって決定され、典型的には、数ｍｍから最大ｃｍまでの範囲において変動する。計算セルのサイズは、十分正確に、巨視的な長さのスケールを分解しなければならないため、最大でも、最小の巨視的サイズより約１桁小さい大きさである。短繊維で強化された部分における、混入された繊維の典型的なサイズは、巨視的計算セルの典型的なサイズよりも１桁または２桁小さいため、巨視的材料挙動のモデリングに関連する繊維特性は、統計的な手法によって記述される。短繊維強化材料では、それに関連する巨視的特性は、（ａ）通常は全体的に（ほぼ）一定である、体積濃度（Ｖｏｌｕｍｅ  Ｃｏｎｃｅｎｔｒａｔｉｏｎ）と、（ｂ）通常は部分的な幾何学形状によって大幅に変化する、各セル内における繊維配向（Ｆｉｂｅｒ  Ｏｒｉｅｎｔａｔｉｏｎ；ＦＯ）の局所的な分布である（本トピックに関するさらなる詳細は、発明を実施するための形態の節１．１および１．２において論じられる）。【００１３】  局所的ＦＯの統計的分布に関する単純化された記述および実際的な目的における適切な記述は、対応する分布関数の低次（すなわち２次および４次）モーメントによって提供される。その数学的構造により、これらのモーメントは、（それぞれ２次および４次の）配向テンソル（ｏｒｉａｔａｔｉｏｎ  ｔｅｎｓｏｒ）で��される。繊維強化部品のＣＡＥシミュレーションのフレームワーク内において、４次テンソルは、巨視的レベルにおける繊維強化材料のレオロジー特性（ｒｈｅｏｌｏｇｉｃａｌ  ｐｒｏｐｅｒｔｙ）ならびに力学的特性を予測するのに必要とされる。それは、これらが４次テンソル特性を有するためである。２次ＦＯテンソルは、単位トレースを有する実数値の対称的な３×３行列であるため、その９個の成分のうち５つのみが独立している。４次ＦＯテンソルの独立成分の数は、（全）対称性（ｔｏｔａｌ  ｓｙｍｍｅｔｒｙ）によって、３４＝８１から１５に減少する。【００１４】  そのモーメントに関してＦＯ分布を記述する数学的モデルは、閉包関係（ｃｌｏｓｕｒｅ  ｒｅｌａｔｉｏｎ）の観点から、２次テンソルに関する４次配向テンソルを近似計算することによって大幅に単純化される。閉包関係は、関数の観点から、このような計算スキームの数学的記述を提供する。関連する計算手続は、特定の仮定下において閉包関係が近似的に有効である場合にのみ、「閉包近似」と言われる。いくつかの閉包近似とともに２次ＦＯテンソルのみを使用する手法は、モールド充填プロセス中の時間および空間における２次ＦＯテンソルの進化をシミュレーションするための、「Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ」型のモデルを導く（詳細は、後の２．２節から２．５節を参照）。【００１５】  非対称性熱可塑性プラスチックの流動に関する非特許文献１は、繊維配向予測アルゴリズムを含む分析を使用して、強化熱可塑性プラスチック製品内における繊維配向分布をモデル化するための方法を開示する。ソフトウェアは、線形３角形要素から成る２次元有限要素メッシュを使用して、３次元モデルを近似している。繊維配向を計算するために、一般的なＨｅｌｅ−Ｓｈａｗ近似およびＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式の変形を使用して、流動場が計算されている。繊維配向、温度、粘度の計算は、３次元解を得るために、各要素の「厚さ」における１９の層において有限差分技術を使用して実行されている。しかしながら、モデルが、面方向の外においては速度成分をシミュレーションできないため（Ｈｅｌｅ−Ｓｈａｗ近似の限界）、このシステムが、忠実な３次元としては記述可能ではないことに留意することが重要である。この文献に説明する方法は、忠実な３次元シミュレーションに適用する場合に、不安定で、要求される処理能力が高過ぎるシミュレーションとなる可能性があり、開発技術者の普通の環境においては使用され得ないシミュレーションとなりえる。
【先行技術文献】
【非特許文献】
【００１６】
【非特許文献１】Ｇｌａｓｓ  ｆｉｂｒｅ  ｏｒｉｅｎｔａｔｉｏｎ  ｗｉｔｈｉｎ  ｉｎｊｅｃｔｉｏｎ  ｍｏｕｌｄｅｄ  ａｕｔｏｍｏｔｉｖｅ  ｐｅｄａｌ−Ｓｉｍｕｌａｔｉｏｎ  ａｎｄ  ｅｘｐｅｒｉｍｅｎｔａｌ  ｓｔｕｄｉｅｓ，Ｂ．Ｒ．Ｗｈｉｔｅｓｉｄｅら，Ｐｌａｓｔｉｃｓ，Ｒｕｂｂｅｒ  ａｎｄ  Ｃｏｍｐｏｓｉｔｅｓ，２０００，Ｖｏｌｕｍｅ  ２９，Ｎｏ．１
【発明の開示】
【００１７】  このような背景から、本発明の目的は、多数の非球形粒子を含む懸濁液でモールドキャビティが充填されるプロセスのシミュレーションにおいて、非球形粒子の配向分布を決定するための方法を提供することにある。この方法は、従来の方法よりも安定し、かつ多くの計算を必要としない。
【００１８】  上記目的は、シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算するための、コンピュータにより実装される方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、（ａ）境界条件を特定することと；（ｂ）初期条件を設定すること；（ｃ）巨視的レベルで流体フロー、熱フロー、および質量輸送を求めるために、前記シミュレーション領域の前記セルの少なくとも一部分に関する質量・運動量・エネルギーの平衡方程式（ｂａｌａｎｃｅ  ｅｑｕａｔｉｏｎ）を解くこと；（ｄ）空間および時間の関数として、巨視的レベルで前記非球形粒子配向における変化を決定するために、前記解かれた平衡方程式の結果に少なくとも部分的に基づいて、非球形粒子配向動的方程式（ｎｏｎｓｐｅｒｉｃａｌ  ｐａｒｔｉｃｌｅ  ｏｒｉｅｎｄａｔｉｏｎ  ｄｙｎａｍｉｃ  ｅｑｕａｔｉｏｎ）を解くことと、を含む方法を提供することによって、請求項１に従って達成される。なお、ステップｄ）において、前記粒子配向方程式は、懸濁液を含む計算セルについてのみ解かれてもよい。
【００１９】  ステップｃ）は、（ｃｃ）前記流体または懸濁液の更新された自由表面または流頭（ｆｌｏｗ  ｆｒｏｎｔ）を決定することをさらに含むことが好ましい。この自由表面は、前記解かれた平衡方程式の結果に少なくとも部分的に基づいて、前記懸濁液で充填されたセルを、前記キャビティの空セルから分離することが好ましい。ステップ（ｃｃ）は、前記更新された流頭に従って、前記境界条件を更新することをさらに含むことが好ましい。また、本発明の方法は、ｅ）前記モールドキャビティが、前記懸濁液のシミュレーション射出によって充填されているかを判断することによって、前記シミュレーションされた射出成形プロセスが完了しているかを決定することと、ｆ）前記シミュレーションされた射出成形プロセスが完了するまでステップｃ），ｃｃ），ｄ）を繰り返すこととをさらに含むことが好ましい。
【００２０】  本発明の上記目的は、多数の非球形粒子を含む溶媒により形成される懸濁液でモールドキャビティが充填されるプロセスのシミュレーションにおいて、非球形粒子の統計的分布配向を記述するための方法であって、（１）前記キャビティの幾何学的形状を定義する３次元コンピュータモデルを提供することと、（２）境界条件を特定することと、（３）複数のセルによるメッシュを形成するために、前記モデルに基づいて解領域を離散化することと、（４）前記解領域の少なくとも一部分に関するエネルギー方程式およびフロー方程式を解くことと、（５）各々の前記セルにおけるフロー条件および温度条件を時間の関数として計算することと、（６）非球形粒子配向の変化を計算することと、（７）各々の前記セルにおける前記非球形粒子の前記配向の統計的分布を時間の関数として記述することとを含む方法を提供することによって、請求項６に従って達成される。
【００２１】  本発明の方法は、部分的な繊維配向分布および熱力学的特性を大幅に低減された計算量で予測するために、懸濁液のフロー（ｆｌｏｗ）を利用する。シミュレーションの結果は、開発技術者によって使用されることができ、それによって、繊維配向に関して製品が最適化され、製品の強度および形状安定性が改善されうるであろう。
【００２２】  多数の非球形粒子を含む懸濁液でモールドキャビティが充填されるプロセスに関する、本発明に係るシミュレーションにおいて、非球形粒子の配向分布を決定するための方法および装置のさらなる目的、特徴、利点、および特性は、次の詳細説明によって、より明らかになるであろう。
【図面の簡単な説明】
【００２３】  本説明の以下の詳細部分において、図面に示す例示的実施形態を参照して、本発明についてより詳細に説明する。
【図１】射出成形装置およびモールドの断面の概略図である。
【図２】本発明のある実施形態に従う、射出成形プロセスの基本プロセスステップを概説するトップレベルのフローチャートである。
【図３】図２のフローチャートのステップ５をさらに詳細に示すフローチャートである。
【図４】繊維強化プラスチック製品の顕微鏡写真である。
【図５】本発明のある実施形態において使用される繊維のモデルを示す。
【図６】特性多項式（ｃｈａｒａｃｔｅｒｉｓｔｉｃ  ｐｏｌｙｎｏｍｉａｌ）の２次項および３次項ならびに線形項および定数項を示すグラフである。
【図７】行列のトリプレットの固有値を示すグラフである。
【図８】「最大限に対称的である」場合に対応する特別な例を示す。
【図９】４面体形状の体積の変形バージョンを示す。
【図１０】本発明のある実施形態に従うＦＯ行列（すなわち２次ＦＯテンソル）の位相空間集合ＭＦＴの全体構造に関する図を示す。
【図１１】本発明のある実施形態に従う単純化形式における演算子分割プロセスを示すフローチャートである。
【図１２】時間ステップ方法の詳細を示すフローチャートである。
【図１３】本発明のある実施形態に従う、トレース再スケーリング（ｔｒａｃｅ  ｒｅｓｃａｌｉｎｇ）を使用する位相空間射影プロセスを示すフローチャートである。
【詳細な説明】
【００２４】  図１は、射出成形装置１を図式的に示す。射出成形装置には、ホッパー（ｈｏｐｐｅｒ）３に配置されるポリマーペレットを備えるネジ２が設けられる。ポリマーペレットは、ネジ２および加熱要素４の作用によって、粘性の塊に変換され、高圧下でモールド６のキャビティ５に注入される。また、成形装置および射出成形製造サイクルは、当技術分野において周知であり、本明細書において詳細には説明しない。成形装置１は、非繊維強化・繊維強化のどちらのプラスチック部品も製造することができる。
【００２５】  コンピュータによる射出成形プロセスの計算シミュレーションは、図２に示すプロセスに従って実行できる。
【００２６】  一般的に表されるシミュレーションの主要なステップは以下の通りである。
ステップ１：シミュレーション領域の幾何学的形状に関するデジタル表現を準備する。
ステップ２：計算領域を多数の小要素への分割するメッシュ化。これら多数の小要素は、（様々な解法アルゴリズムを利用して、）微分方程式を離散化するための基礎である。シミュレーションされる物理的現象に対する解は、このようにして探索される。
ステップ３：異なる材料領域に関して必要な物理的データをシミュレーションモデル（データベースまたはデータバンク）に取り入れる。
ステップ４：シミュレーション計画の境界条件を特定する。
ステップ５：射出成形プロセスをシミュレーションする（このステップは後に詳細に説明される）。
ステップ６：ワークステーション等のコンピュータのディスプレイ上に、グラフ提示または数値的提示として結果を表示する。
【００２７】  繊維強化製品の射出成形プロセスをシミュレーションするステップ５の詳細について、図３のフローチャートに示す。プロセスのこの部分では、熱フロー、流体フロー、応力、歪みに関する微分方程式が、以下の計算アルゴリズムを使用して解かれる。
ステップ１、熱物理材料特性および流頭（ｆｌｏｗ  ｆｒｏｎｔ）の初期条件が設定される。
ステップ２、領域全体の熱方程式および全流体セルに関するフロー方程式が、質量、エネルギー、運動量保存の式を使用して解かれる。
ステップ３、このステップにおいては流頭が移動し、新しい流頭に従って境界条件が選択される。
ステップ４では、繊維配向および繊維輸送（ｆｉｂｅｒ  ｔｒａｎｓｐｏｒｔ）が、先行するステップにおいて得られた流速から計算される。このステップは、後で詳細に説明される。初期条件が、新しく充填されたセルに適用される。流体物質を含むセルのみが考慮される。
ステップ５では、化学反応等の追加の量（ｑｕａｎｔｉｔｙ）が計算され、セルが凝固する場合に検証される。
ステップ６では、モールドにおいて射出成形プロセスが完了しているか否かが検証され、完了していない場合、シミュレーションは次の時間ステップへの進み、プロセスはステップ２に戻る。
ステップ７では、モールドから取り出した成型品の特性が計算される。
モールドの外部の熱物理材料温度、すなわち、温度、鋳引け、反り等が、射出成形シミュレーションから得た情報を使用して計算される。
【００２８】  図４は、繊維強化プラスチック製品の顕微鏡写真であり、射出成形プロセス完了後の繊維の配向が示されている。最終生産物における繊維の向きは、射出成形プロセス中の塊状熱可塑性プラスチックのフローパターンに大きく依存する。繊維の向きは、各々の繊維について厳密に決定されるのではなく、むしろ分布関数によって記述される。
【００２９】  繊維の向きの計算に関する詳細に進む前に、短繊維熱可塑性プラスチック溶融液の基礎的な流体力学的側面に関して概説する。
【００３０】  １．短繊維熱可塑性プラスチック溶融液の基礎的な流体力学的側面
【００３１】  ポリマーの塊（ｐｏｌｙｍｅｒ  ｍａｓｓ）の中に、多数の短繊維（ｓｈｏｒｔ  ｆｉｂｅｒ）を分散させて混入することにより、成型品を短繊維強化熱可塑性プラスチックとすることができる。融解状態では、すなわち熱可塑性プラスチックマトリクスが液体であるのに十分高い温度である場合、プラスチック溶融液と分散繊維との混合液は、流体力学および流体物理学の専門用語において一般的に粒子懸濁液（ｐａｒｔｉｃｌｅ  ｓｕｓｐｅｎｓｉｏｎ）と言われる複合流体を構成する。一般に、このような懸濁液は、２つの異なる相、つまり（ｉ）分散繊維を含まない融解プラスチック材料に対応する溶媒相、および（ｉｉ）溶媒に浸漬される全ての繊維から成る粒子相から構成される。
【００３２】  １．１  繊維の幾何学的形状および材料特性
【００３３】  回転対称（図５に示す）の軸を有するスフェリカルな粒子は、以下において繊維（ｆｉｂｅｒ）と記述される。用語の粒子（ｐａｒｔｉｃｌｅ）および繊維（ｆｉｂｅｒ）は、別途明示的に示す場合を除
2010531752.tif 000003
ラスチック材料に分散される繊維（例えば、炭素またはガラス型の繊維）に関するこれら
2010531752.tif 000004
比は大きくなる。
【００３４】  繊維強化プラスチックに浸漬される大部分の種類の繊維では、繊維材料の質量密度は、懸濁プラスチック溶融液の密度に似ている。プラスチック溶融液自体は、モールド充填プロセス中に生じる典型的な温度においてかなり高い粘度を有する。これらの２つの側面が組み合わされるので、浮力効果や慣性効果は、好適な実施形態に従う方法において、プラスチック溶融液内における繊維の動きを記述する際に無視することができる。
【００３５】  繊維濃度に関し、懸濁液における繊維は空間的に均一の分布を有すると仮定され、このため濃度は一定であると解釈される。しかし、（ｉ）入口に既に存在する不均一な分布の純枠に対流的な輸送、（ｉｉ）モールド充填段階中のせん断誘起粒子移動効果の発生、という作用は、これらを変化させることがある。
【００３６】  （ｉ）または（ｉｉ）の効果が大きい場合、熱可塑性プラスチック溶融液に浸漬される短繊維の懸濁液フローを記述するために使用されるモデルは、必然的に、粒子および溶媒の両方の相のフローを相互結合し、かつ不均一な粒子濃度を可能にする「２相フロー」モデルでなければならない。しかしながら、これはこの例に当てはまらない。現象（ｉｉ）は拡散プロセスの一種として理解可能であり、この場合、拡散定数は、懸濁粒子の実サイズ（すなわち、繊維の長さ、または球形粒子の場合は半径）の、フローキャビティ（例えば、壁厚）の典型的なサイズに対する比率として規定される相対的な粒子サイズの２乗の大きさである。相対的な粒子サイズは、繊維強化プラスチックからの射出成形によって製造される薄い壁部分であっても、典型的には０．１から０．０１のオーダーであるため、この場合、（ｉｉ）を完仝に無視することは全く妥当である。（ｉ）に関し、ネジにおける溶融液の準備プロセスにより、入口における不均一な濃度プロファイルが存在し得ることが想定され得る。しかしながら、典型的な短繊維強化プラスチック部分について実験的に実際に認められる繊維濃度は、その部分内のほぼどこにおいても一定値から極めて小さい偏移しか示さない。従って、繊維分布が均一であるとの仮定は妥当である［２０］。
【００３７】  １．２  レオロジー特性に対する繊維の影響
【００３８】  熱可塑性プラスチック材料は、非ニュートン的なフロー挙動を呈する。繊維を含まない純枠なプラスチック溶融液の場合、材料のレオロジー特性（ｒｈｅｏｌｏｇｉｃａｌ  ｐｒｏｐｅｒｔｙ）は、スカラー粘度関数（ｓｃａｌａｒ  ｖｉｓｃｏｓｉｔｙ  ｆｕｎｃｔｉｏｎ）によってモデル化される。このスカラー粘度関数は、温度ならびにフロー（一般的ニュートン流体）の状態変数に依存し得る。レオロジー特性は、非ニュートン的なフロー特性（例えば、せん断減粘（ｓｈｅａｒ  ｔｈｉｎｎｉｎｇ））の限られた範囲にしか及ばないことが知られているが、射出成形シミュレーションの目的では有用かつ満足しうる正確性を有することが証明されている。この種類のモデルは、本発明の好滴な実施形態に従う方法においても使用される。
【００３９】  熱可塑性プラスチックマトリクス材に繊維を添加することによって、固体状態において、異方性が強く、埋め込まれた繊維の配向の局所的な分布に強く依存する力学的特性がもたらされる。原理上は、材料が液体（すなわち、融解）状態である場合にも、異方的な材料挙動が示される。この異方性を説明するためには、上述のスカラー粘度を粘度テンソルに置き換えなければならない。繊維懸濁液の材料挙動をモデル化する両方の方法を比較する多くの研究が行なわれている。モールド充填中に見られるほとんどの流動状況では、両方の種類のモデル（すなわち、スカラーモデルおよびテンソルモデル）によって予測される充填パターンが、ほとんど差を示さないことが分かる（例えば、［２１］、［２２］を参照）。ゆえに、異方的な粘度効果は無視され、熱可塑性プラスチック溶融液が繊維を含む場合においても、単純な一般的ニュートンモデルが使用される。これは、溶融液のレオロジー特性に対する分散繊維の配向の影響を無視することに等しい。粘度は、繊維濃度に依存するが、これは一定であると仮定される（上記参照）ことから、この側面は、可塑性プラスチック溶融液の効果的な材料特性に寄与する先験的な既知の材料パラメータとしてのみ関与する。
【００４０】  このアプローチのため、好適な実施形態に従う方法は、フローの計算と繊維配向の計算とを部分的に分けて行う。流速は、分散繊維の向きに局所的に影響を及ぼすが、フローに対する繊維配向の影響はごくわずかである。ゆえに、フローの計算は、繊維配向の計算とは独立に実行される。このように、溶融液の局所的な流速は、外部係数の組として繊維配向の計算に使用されるモデルに関与する。
【００４１】  ２．Ｆｏｌｇａｒ−Ｔｕｃｋｅｒモデル
【００４２】  ２．１  Ｊｅｆｆｅｒｙの方程式
【００４３】
2010531752.tif 000005
細長い粒子である。繊維は、溶媒の局所フロー場において繊維に作用する力学的力が、実質的にいかなる変形をも与えることができないほど十分小さい。ゆえに、個々の繊維は、本明細書において、細長く、回転対称の剛体であるようにモデル化され、繊維の向きは、回転対称軸に沿う向きを有する単位ベクトルｐによって求められる。２つのベクトル±ｐの両方は、粒子の同一の配向状態を表す。
【００４４】  図５は、回転楕円形状の剛体粒子を示す。この粒子の動きは、粒子の近傍における速度
2010531752.tif 000006
に加え、配向ベクトルｐが示されている。流速Ｕの方向およびサイズは、対応する矢印の方向および長さによって示される。速度ベクトルの方向は、全く同じであるが、その長さは同じではなく、粒子近傍の流速が一定でないことが示される。しかしながら、底部から上部へのベクトルの長さ方向において増加の程度が一定であるため、流速の変化の局所的な大きさ、すなわち速度勾配は、同一であるように思われる。粒子表面の各点において、粒子は、局所的な流速と同じように移動するはずである（「滑り無し」境界条件）。粒子の周辺の流速が一定であれば、これは、単純な並進運動（ｔｒａｎｓｌａｔｉｏｎａｌ  ｍｏｔｉｏｎ）、つまり、純粋に対流的な輸送（ｃｏｎｖｅｃｔｉｖｅ  ｔｒａｎｓｐｏｒｔ）をもたらすであろう。そうでない場合は、速度勾配の存在下において、粒子は、破線矢印によって示されるように、回転運動（ｒｏｔａｔｉｏｎａｌ  ｍｏｔｉｏｎ）も行う。全体として、剛体粒子によって実行される最も一般的な種類の運動は、粒子の周囲における「平均速度」による並進（ｔｒａｎｓｌａｔｉｏｎ）と、粒子の近傍におけるその平均値から流速の局所偏差を記述する速度勾配によりもたらされる、質量中心の周りの回転（ｒｏｔａｔｉｏｎ）である。
【００４５】  粘性力が優勢であり、そのため慣性力がわずかであり、単一の繊維が動く範囲における局所的フローの変動が小さいことを仮定することによって、上述の定性的記述は、Ｊｅｆｆｅｒｙ［２］による数学的モデルに組み込まれる。上述のように、熱可塑性プラスチック溶融液中に分散される短繊維の場合、これらの全ての仮定は満たされる。繊維の周囲の流速の
2010531752.tif 000007
述するのに十分であることを含意している。このテンソルは、３×３行列であり、その要素は、空間における点ｒ＝（χ１，χ２，χ３）Ｔのデカルト座標χｉに対する流速ベクトルの成分Ｕの偏導
2010531752.tif 000008
子の中心の空間座標である場合、粒子近傍の点ｒおよび時間ｔにおける流速Ｕは、１次テ
2010531752.tif 000009
近似される。また速度勾配は、局所的に一定量としてみなされうる。すなわち、ｒ０近傍の
2010531752.tif 000010
【００４６】  この場合、空間座標および時間座標の関数としてその配向ベクトルｐ（ｒ，ｔ）によって求められる繊維の過渡的な配向状態（ｔｒａｎｓｉｅｎｔ  ｏｒｉｅｎｔａｔｉｏｎ  ｓｔａｔｅ）は、以下のＪｅｆｆｅｒｙの方程式を使用して計算される。
【００４７】
2010531752.tif 000011
【００４８】  この方程式は、コンパクトオイラー形式（ｃｏｍｐａｃｔ  ｅｕｌｅｒｉａｎ  ｆｏｒｍ）で本明細書に記載されている。方程式（１）の左辺（ｌ．ｈ．ｓ．）における対流（または材料）微分
2010531752.tif 000012
は、フローの局所速度場Ｕ（ｒ，ｔ）における純粋な繊維配向、ＦＯ輸送を記述し、右辺（ｒ．ｈ．ｓ．）
2010531752.tif 000013
無限に細長い繊維（すなわちアスペクト比ｒａ→∞）である理想的な特別の場合では、テン
2010531752.tif 000014
する繊維の一般的な場合では、それは、繊維の幾何学的形状パラメータλ＝（ｒａ２−１）／（ｒａ２＋１）を含む以下の項によって定義される。
【００４９】
2010531752.tif 000015
【００５０】  このパラメータは、粒子幾何学的形状が回転対称楕円形粒子である場合に、それが繊維の回転運動にどのように影響を及ぼすかをコードしている。
【００５１】  Ｊｅｆｆｅｒｙの方程式を記述する他の方法は、速度勾配テンソルのその対称部分および反対称部分への一意的に定義された分解を利用することである。対称部分および反対称部分は、以下によって求められるせん断速度（ｓｈｅａｒ  ｒａｔｅ）および渦度テンソル（ｖｏｒｔｉｃｉｔｙ  ｔｅｎｓｏｒ）によって与えられる。
【００５２】
2010531752.tif 000016
【００５３】  せん断速度および渦度テンソルは、以下の式を介して有効速度勾配テンソルに関連付けられる。
【００５４】
2010531752.tif 000017
【００５５】
2010531752.tif 000018
式の代替形式が以下のように得られる。
【００５６】
2010531752.tif 000019
【００５７】  Ｊｅｆｆｅｒｙの方程式のさらに別の形式は、渦度テンソル（ｖｏｒｔｉｃｉｔｙ  ｔｅｎｓｏｒ）を、流速ベクトル
2010531752.tif 000020
おける第１の項を書き直すことによって得ることができる（例えば、［３］参照）。
【００５８】  幾何学的形状が回転対称であるが、楕円形ではない粒子の場合、方程式（１）の形式ならびに有効速度勾配の定義は、同一のままであるが、幾何学的形状パラメータλの公式は、
2010531752.tif 000021
換えられるように変更しなければならず、その適切な値は、実験により決定しなければならない。この意味では、幾何学的形状パラメータλは、材料パラメータと見なすことができる。非対称非球面粒子の場合、方程式（１）は、３つの幾何学的形状パラメータを含む３つのＪｅｆｆｅｒｙ型方程式の結合系によって置き換えられなければならない［３］。
【００５９】  ２．２  繊維配向の巨視的分布
【００６０】  マクロなレベルで実行される短繊維強化プラスチックの射出成形シミュレーションの場合、単一の計算セル（すなわち計算領域の体積要素（Ｖｏｌｕｍｅ  Ｅｌｅｍｅｎｔ））に含まれる繊維の数は、個々の繊維のサイズが小さいため、非常に多い。これを図４に示す顕微鏡写真によって示す。これは、短繊維強化プラスチックから制作された部品から採取される典型的なサンプルである。ゆえに、計算セルに含まれる材料の局所的ＦＯ状態は、巨視的モデル、すなわち、ＦＯベクトルｐの分布関数Ψ（ｐ，ｒ，ｔ）によって記述される（［５］を参照）。繊維濃度は均一であると仮定されるが、分布関数Ψは、依然として局所的な流速場Ｕ（ｒ，ｔ）およびそ
2010531752.tif 000022
加の点として、その配向に影響を及ぼす繊維の相互作用を説明する項を、マクロレベルでモデルに含めなければならない。
【００６１】  個々の繊維が相互作用しない巨視的モデルとしてのＪｅｆｆｅｒｙの方程式から、相互に作用する多くの繊維のＦＯ統計を与える巨視的モデルへの移行は、以下のＦＯ分布関数Ψを従属変数として有する、対応するＦｏｋｋｅｒ−Ｐｌａｎｃｋ方程式を介して達成される。
【００６２】
2010531752.tif 000023
【００６３】
2010531752.tif 000024
演算子を表す。ＦｏｌｇａｒおよびＴｕｃｋｅｒの手法［４］によると、局所速度勾配の有効せん断速度
2010531752.tif 000025
懸濁液における繊維間相互作用の単純化モデルを得ることができる。式中、Ｇｉｊは、式（１ｂ）
2010531752.tif 000026
る。無次元で非負である相互作用係数Ｃｌは、懸濁液の材料パラメータである。典型的には、この係数は、短繊維強化プラスチックについて、１０−３．．．１０−２の範囲における小さい値（正値）を有する。非圧縮性（ならびにほぼ非圧縮性）のフローでは、「確率的」拡散項の相対的な弱さ（「決定論的な」Ｊｅｆｆｅｒｙ力学を表す項と比較して）は、安定性問題の源であり得ることに留意されたい。
【００６４】  ２．３  繊維配向テンソルおよびＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式
【００６５】  Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式（２）によって局所的ＦＯ分布を計算することは、フローシミュレーション領域の計算セル毎の単位球面Ｓ２について定義されるＰＥＤの数値解を必要とし、これは、「工業規模」の３Ｄ問題について極めて高価なタスクである。ゆえに、ＡｄｖａｎｉおよびＴｕｃｋｅｒ［６］は、分布関数のモーメントとして定義される繊維配向テンソルの利用すること、すなわち、Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式を、ＦＯテンソルのモーメント方程式の階層に置き換えることを提案した。方向±ｐが同一の配向状態に対応するという事実を反映する、変数ｐに対するＦＯ分布の反転対称性Ψ（−ｐ，．．．）＝Ψ（ｐ，．．．）により、全ての奇数次モーメントは、一様にゼロになり、結果的に、Ψのモーメント展開は、偶数次の要素のみを含む。ゆえに、この展開の第１の非自明なモーメント（ｎｏｎｔｒｉｖｉａｌ  ｍｏｍｅｎｔ）は、
【００６６】
2010531752.tif 000027
によって求められる第２のモーメントである。
2010531752.tif 000028
【００６７】
2010531752.tif 000029
行列）と呼ばれ、定義により、実対称３×３行列である。ＦＯ分布（これも定義による）は、
2010531752.tif 000030
である。展開階層（ｅｘｐａｎｓｉｏｎ  ｈｉｅｒａｒｃｈｙ）における次の非自明なモーメントは、
【００６８】
2010531752.tif 000031
【００６９】
2010531752.tif 000032
【００７０】  上述のモーメント方程式の階層（ｈｉｅｒａｒｃｈｙ）は、微分演算および積分演算を置き換えることによって、各次数のモーメントについて得られる。すなわち、
2010531752.tif 000033
（Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ  Ｅｑｕａｔｉｏｎ，ＦＴＥ）は、階層的に体系化された一連の以下のモーメント方程式において最初の非自明的であるものとして得られる。
【００７１】
2010531752.tif 000034
【００７２】  Ｊｅｆｆｅｒｙの方程式のように、（５）の左辺における対流微分は、オイラー基準座標系における繊維の並進運動による局所的ＦＯ状態（モーメント展開のこの次数におけるＦＯ行列
2010531752.tif 000035
の右辺の第３の項は、Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式（２）における拡散項の存在によりもたらされ
2010531752.tif 000036
減衰効果（ｄａｍｐｉｎｇ  ｅｆｆｅｃｔ）をもたらす。
【００７３】  オイラー形式（５）では、ＦＴＥは、対流−反応型の１次ＰＥＤの結合系である。ラグランジ
2010531752.tif 000037
【００７４】  ２．４  ＦＯ行列の基礎特性
【００７５】
2010531752.tif 000038
【００７６】
2010531752.tif 000039
（ｕｎｉａｘｉａｌ  ｏｒｉａｔａｔｉｏｎ  ｓｔａｔｅ）を表し、この場合、繊維の１００％は、±ｐ０により与えられる方向に向いている。したがって、ｐ０の符号は、この特別なＦＯ行列を定義するダイアディック
2010531752.tif 000040
であり、残りの２つの固有値は、ｐ０に直交する平面に存在する対応するベクトルを有するゼロであり、そうでない場合は任意である。
【００７７】
2010531752.tif 000041
の項は、固有値により与えられる重み付けにより、ＦＯ行列の固有ベクトルによって形成さ
2010531752.tif 000042
クトルＥｋの方向に沿って配向される繊維の局所分数として固有値μｋを解釈することが可能になる。この意味では、ＦＯ行列のスペクトルデータ｛μｋ，Ｅｋ｝ｋ＝１，２，３は、小さな体積の懸濁液中における繊維の局所的マクロ配向状態を表す。
【００７８】  これは、ＦＯ行列および位相空間に関してＦＴＢを（数式的に）定義することの動機付けとしての役割を果たす［１２］。
【００７９】
2010531752.tif 000043
のトレースが１に等しい場合に（そしてその場合に限り）、ＦＯ行列である。ＦＴＢの位相空間ＭＦＴは、全ＦＯ行列の組である。
【００８０】  適切に正規化される分布関数か、そうでない場合は任意の分布関数を使用して、集合ＭＦＴが、タイプ（３）のモーメント積分からもたらされる全実対称３×３行列の集合に等しいことが分かる。そのトポロジカル特性および幾何学的特性による位相空間ＭＦＴの数学的特徴付けは、著者のうちの１人によって最近提供されている（［１２］参照）。ＦＴＥの数値積
2010531752.tif 000044
間ステップも）中に上述の特別なスペクトル特性を有することが必要とされるため、行列が集合ＭＦＴに属するか否かを厳密に把握することが、実際的に非常に重要である。
【００８１】  ２．５  閉包問題
【００８２】  いわゆる閉包問題（ｃｌｏｓｕｒｅ  ｐｒｏｂｌｅｍ）は、モーメント展開の各次数において、モーメン
2010531752.tif 000045
系を解くことができない。閉包近似をＦＴＥに適用することは、（５）の右辺における厳密
2010531752.tif 000046
【００８３】  既知の例は、以下の式によって定義されるハイブリッド閉包（ｈｙｂｒｉｄ  ｃｌｏｓｕｒｅ）［７］である。
【００８４】
2010531752.tif 000047
この式は、いくつかの周知の不利点にもかかわらず、その代数的簡易性および数値的ロバスト性により、認められた選択である［７］。
【００８５】
2010531752.tif 000048
【００８６】
2010531752.tif 000049
として定義される。
【００８７】  これは、一軸配向分布の特別な場合において厳密な結果をもたらし、線形閉包は次のように与えられる。
【００８８】
2010531752.tif 000050
これは、等方性（ｉｓｏｔｒｏｐｉｃ）の場合に厳密である。これらの極端な場合の間の補間重み（ｉｎｔｅｒｐｏｒａｔｉｏｎ  ｗｅｉｇｈｔ）は、以下のスカラー配向因子によって提供される。
【００８９】
2010531752.tif 000051
【００９０】
2010531752.tif 000052
てはまる。
【００９１】  ハイブリッド閉包近似（６）と組み合わせられるＦＴＥ（５）の特別な変形を、ＦＴＥ−ｈｙｂと表す。この変形は、非常に実際的な利益を有する。本発明の好適な実施形態は、Ｆｏｌｇａｒ−ＴｕｃｋｅｒモデルのＦＴＥ−ｈｙｂ変形にも基づいている。
【００９２】  ２．６  微分代数系としてのＦＴＥ
【００９３】  （５）の右辺は、任意の実対称行列について数式的によく定義されている。（５）の右辺の特性が解軌道を領域ＭＦＴに自動的に制限するか否かに関する質問は、少なくともその厳密な形式（すなわち、閉包近似を含まない）において、ＦＴＥについては肯定的に回答でき
2010531752.tif 000053
Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ方程式（２）の解Ψを使用してモーメント積分（３）を評価することによって得られるＦＯ分布の厳密な２次モーメントに等しい。この議論によって、厳密な４次ＦＯ
2010531752.tif 000054
軌道が必然的に領域ＭＦＴに制限されることが結論付けられる。
【００９４】  しかしながら、この議論は、閉包近似が適用される場合は有効ではないし、完全分布関数に関する先験的な知識を持たずに数値的にＦＴＥを解かなければならない場合には、閉包近似が常に必要である。従って、ＦＴＥの解の位相空間ＭＦＴへの制限に関する問題は、実用性を有する全ての問題において常に発生する。
【００９５】
2010531752.tif 000055
れる。それは、その不変量の代数不等式によって定式化されうる（［１２］参照）。従属変数に対するこのような制約は、ＦＴＥを微分代数系（ＤＡＳ）に変更し、数値積分に使用される手順において留意しなければならない。
【００９６】  ３．ＦＴＥの位相空間の数学的特徴付け
【００９７】  位相空間ＭＦＴの基礎的なトポロジカル特性は、以下に概説されるように、ＦＯ行列の特別な特性から直接推測できる。
【００９８】
2010531752.tif 000056
限される実対称３×３行列のベクトル空間の有界凸部分集合（ｂｏｕｄｅｄ  ｃｏｎｖｅｘ  ｓｕｂｓｅｔ）である。
【００９９】  ＦＴＥに適切な積分方法を構築するために、任意の適切なＯＤＥ積分と組み合わせられうるＭＦＴへの射影マッピングの定義が、ＭＦＴの凸性（ｃｏｎｖｅｘｉｔｙ）によって可能になる（［９］の章ＩＶ．４参照）。ＦＯ行列の不変代数的特徴（ｉｎｖａｒｉａｎｔ  ａｌｇｅｂｒａｉｃ  ｃｈａｒａｃｔｅｒｉｚａｔｉｏｎ）は、次の実対称３×３行列の特性多項式：
2010531752.tif 000057
の分析によって得られうる。
【０１００】
2010531752.tif 000058
特徴は、以下の定理に従ってこれらの不変量の項において定式化され得る。
【０１０１】
2010531752.tif 000059
およびＤａが非負である場合およびその場合に限り、ＦＯ行列である。
【０１０２】  図６は、特性多項式Ｐａ（μ）の２次項および３次項ならびに線形項および定数項の別々の分析によってこの定理を説明しており、正のトレースＳａは、明らかに正の固有値μｋの存在をもたらし、一方２つの不変ＫａおよびＤａの非負値は、負の固有値μｋの存在を阻止する。
2010531752.tif 000060
ることによって、これは、上記定理の証明を完成させる。
【０１０３】  定理１に従い、５Ｄ物体である位相空間ＭＦＴの幾何学的形状の概念は、ＦＯ行列の対角要素および非対角要素を別々に考察することによって得られる。まず、
2010531752.tif 000061
角形集合（ｔｒｉａｎｇｌｅ  ｓｅｔ）に制限される（図７参照）。
【０１０４】
2010531752.tif 000062
【０１０５】  これは、その対角要素に対するＦＯ行列を特徴付ける必要な条件を提供し、位相空間集合ＭＦＴの「対角部分」の不変表現をもたらす。
【０１０６】  領域ＭＦＴの「非対角」部分の記述形式は、実対称３×３行列の対角要素および非対角要素のトリプレットについて記法（χ，ｙ，ｚ）および（ｕ，ν，ｗ）を導入することによって得られる。ただし、任意の（しかし固定の）対角トリプレット（χ，ｙ，ｚ）∈Δａを考慮する。固定対角トリプレットに属する全ての「許容可能な」非対角トリプレットの以下の集合を数式的に定義する。
【０１０７】
2010531752.tif 000063
【０１０８】  代数的には、集合Ｎ（χ，ｙ，ｚ）は、上記の定理２に関して説明したように、以下の不等式のペアを同時に満たす全ての非対角トリプレット（ｕ，ν，ｗ）の集合として特徴付けられることができる。
【０１０９】
2010531752.tif 000064
【０１１０】
2010531752.tif 000065
【０１１１】  「最大限に対称的（ｍａｘｉｍａｌｌｙ  ｓｙｍｍｅｔｒｉｃ）である」場合のｘ＝ｙ＝ｚ＝１／３に対応する特別な例について図８に示し、他の場合（図９参照）は、この４面体形状ボリュームの変形バージョンに対応する。
【０１１２】  不変量ＫａおよびＤａの非負値への制限（９）および（１０）を補足するために、正のトレー
2010531752.tif 000066
【０１１３】  集合ΔａおよびＮ（χ，ｙ，ｚ）を導入することによって得られるファイブレーション（ｆｉｂｒａｔｉｏｎ）
2010531752.tif 000067
理解するために役立つ。つまり、個々の「繊維」は、許容可能な非対角トリプレットの集合Ｎ（χ，ｙ，ｚ）をその「基点」（χ，ｙ，ｚ）∈Δａにローカルに関連付けるという手順によって、可視化される。そして、３角形集合Δａにおけるこの基点の後の変動は、位相空間全体をカバーする。
【０１１４】  この節において概説される結果は、位相空間ＭＦＴが複雑な数学的オブジェクトであることを明確に示す。定理２に記載される厳密な数学的結果によると、対称３×３行列の実ベクトル空間において発展方程式（ｅｖｏｌｕｔｉｏｎ  ｅｑｕａｔｉｏｎ）として必然的に定義されるＦＴＥの任
2010531752.tif 000068
不変不等式（９）および（１０）の実現を必然的に含意する。この事実によって、必然的にＦＴＥはＤＡＳに変わる。単位トレース条件は、４次ＦＯテンソルの如何なる閉包近似も存在しない場合には、それが満足されることが厳密なＦＴＥに「組み込まれ」、様々な閉包近似（第４節参照）が行われるより一般的な状況の下であっても、依然として有効である。しかし、不変不等式（９）および（１０）の有効性は、概して、現在知られている如何なる閉包近似の存在下においても維持されない（また、「厳密な」ＦＴＥの場合の（９）および（１０）の有効性は、ＦＴＥ自体の代数構造によるよりも、むしろ２．６節に提供される間接的な理由によって推論され得る）。これらの数学的事実は、学術的な著作に通常示される単純な例において見落とされることが多い。
【０１１５】  しかしながら、閉包近似は、ＦＴＥを数値的に解く任意の手順の本質的要素として必然的に含まれる。ゆえに、典型的には工業用途において発生する状況等のより複雑な状況に対処するのに適切なあらゆるシミュレーション手続は、解軌道をその理論的に許容可能な領域に制限する適切な制御手続によって、これに対処しなければならない。したがって、ＦＴＥを微分代数系として適切に取り扱うことは、重要な工業用途に必須である。
【０１１６】  ４．トレース保存および安定性に関する問題
【０１１７】  閉包近似をＦＴＥに適用することは、（５）の右辺における厳密な（しかし未知の）４次
2010531752.tif 000069
さらなる詳細については、閉包問題に関する上の節を参照されたい）。閉包の選択に依存
2010531752.tif 000070
【０１１８】
2010531752.tif 000071
要性は、以下の式の組によって求められる、いわゆる直交異方性対称性（ｏｒｔｈｏｔｒｏｐｉｃ  ｓｙｍｍｅｔｒｙ）である。
【０１１９】
2010531752.tif 000072
【０１２０】
2010531752.tif 000073
とを仮定する場合、そのトレースに関して、以下の微分方程式（略称：ＤＥ）を導き出すことできる（［１７］参照）。
【０１２１】
2010531752.tif 000074
【０１２２】  ＦＴＥの解に関するトレースの安定性（ｓｔａｂｉｌｉｔｙ）および保存性（ｃｏｎｓｅｖａｔｉｏｎ）を、この方
2010531752.tif 000075
においてまだ使用されていない。閉包近似がＦＴＥに適用される場合、後者は、次の正規化条件を更に満たし、
【０１２３】
2010531752.tif 000076
【０１２４】  トレースのＤＥ（１２）は以下のように単純化される。
【０１２５】
2010531752.tif 000077
【０１２６】  拡散パラメータＤｒは、定義により非負であるため、トレース条件により定義される超平面（ｈｙｐｅｒｐｌａｎｅ）が、安定積分多様体（ｓｔａｂｌｅ  ｉｎｔｅｇｒａｌ  ｍａｎｉｆｏｌｄ）であること（Ｄｒ＞０の場合）、
2010531752.tif 000078
レースがＦＴＥの第１の積分であること（Ｄｒ＝０の場合）が結論付けられ得る（厳密な４次
2010531752.tif 000079
意されたい）。
【０１２７】  上の考察は、ＦＴＥへの閉包近似の適用がトレース条件の有効性にいかに影響を及ぼし得
2010531752.tif 000080
単純なものである。
【０１２８】  ハイブリッド閉包（６）が対称特性（１１）のみ有し、正規化条件（１３）を満たさないという特別の場合は、トレースの安定性を慎重に考察しなければならない。このような閉包近似の場合、トレースのＤＥは、以下の形式を仮定する（［１７］も参照）。
【０１２９】
2010531752.tif 000081
【０１３０】  これは次の前因子（ｐｒｅｆａｃｔｏｒ）を含む。
【０１３１】
2010531752.tif 000082
【０１３２】  この前因子は、（１３）を満たす全ての閉包について有効である６Ｄｒの単純化された表現から外れたものであり、ハイブリッド閉包が正規化条件（１３）を満たさないという事実を反映する。特定の条件において前因子φａに出現する追加の項が、この項を負にし得るこ
2010531752.tif 000083
を有するＦＴＥの積分多様体であるが、局所的に不安的であり得ることを意味する。
【０１３３】  前因子φａが負になる可能性は、以下の４段階の手順の議論によって推論できる。
2010531752.tif 000084
という特別な場合では、上記前因子は、単純化された形式である
2010531752.tif 000085
される局所的配向状態がほぼ一軸である場合（すなわち、その固有値μｋのうちの１つが１
2010531752.tif 000086
の値もほぼ１である。従って、特定の仮定下におけるφａの正確な値を近似する、単純化
2010531752.tif 000087
得られる。
2010531752.tif 000088
値を有する少なくとも１つの負の固有値を必然的に有する。
【０１３４】  この事実は、前因子φａの近似式の微分においてこれまでに仮定された種々の状況下に
2010531752.tif 000089
を使用して、この近似式は、以下の形式で記載され得る。
2010531752.tif 000090
【０１３５】
2010531752.tif 000091
固有ベクトルＥｋ（これはせん断速度テンソルの最大の負の固有値の固有方向に沿う方向を有する）を有するＦＯ行列で評価される場合に、φａが負になることが分かる。この場合、
2010531752.tif 000092
が得られる。
【０１３６】  上記議論は厳密に定式化できる。つまり、任意の非圧縮性フローおよび相互作用係数の
2010531752.tif 000093
なるため、ハイブリッド閉包を含むＦＴＥの解のトレースが、位相空間ＭＦＴの特定の領域において局所的に不安定になるという、数学的な証拠を提供することができる。
【０１３７】  したがって、ＦＴＥの数値積分中にトレース安定性問題に留意することが明らかに必要である。
【０１３８】  ５．ＦＴＥの数値積分：一般的側面
【０１３９】  （任意の種類の閉包近似を含む）ＦＴＥの数値積分に適切な方法を選択することは、数式のタイプの数学的分類と、ＦＴＥの特定の代数形式に関する側面とに依存する。閉包ＦＴＥの一般的構造は、以下の段落で詳述するように、従属変数としてのＦＯ行列の要素ａｉｊ（２）の対流−反応型（ｃｏｎｖｅｃｔｉｏｎ−ｒｅａｃｔｉｏｎ  ｔｙｐｅ）の双曲型偏微分方程式（ｈｙｐｅｒｂｏｌｉｃ  ｐａｒｔｉａｌ  ｄｉｆｆｅｒｅｎｔｉａｌ  ｅｑｕａｔｉｏｎ；ＰＤＥ）の結合系を構成することを示す。
【０１４０】  式（５）に与えられるＦＴＥの形式から開始する。この式は、左辺に物質導関数
2010531752.tif 000094
て次のような式が得られる。この式は、閉包ＦＴＥと同等である。
【０１４１】
2010531752.tif 000095
【０１４２】  ＦＴＥのこの形式は、その数学的構造の大部分を既に示している。つまり、左辺は、式（１６）の従属変数であるＦＯ行列の要素ａｉｊ（２）の、分離された純粋な対流輸送を支配する局所的な流速Ｕ（ｒ，ｔ）を含む、輸送型（ｔｒａｎｓｐｏｔ  ｔｙｐｅ）または対流型（ｃｏｎｖｅｃｔｉｏｎ  ｔｙｐｅ）の単純な１次偏微分演算子から成る。式（１６）の右辺における種々の項の代数構造は、右辺が、次のよ
2010531752.tif 000096
称テンソル関数を、（ＦＯ行列自体と同様に、また数学的整合性による必要に応じて、）構成することを示す。
【０１４３】
2010531752.tif 000097
【０１４４】
2010531752.tif 000098
2010531752.tif 000099
【０１４５】
2010531752.tif 000100
存性を式（１ａ）を用いて解き、有効せん断速度に関して回転拡散パラメータの定義方程式Ｄｒ
2010531752.tif 000101
【０１４６】
2010531752.tif 000102
【０１４７】  閉包ＦＴＥのこの形式は、一定モデルパラメータλおよびＣ１に対する右辺の依存性を明示的に示す。しかし、ＦＯ行列に対する明示的依存性と黙示的依存性との間を依然とし区別
2010531752.tif 000103
に対する黙示的依存性との間も依然として区別している。文献によっては、式（１７ｂ）は、場合により添字記法を使用する成分形式、すなわち
2010531752.tif 000104
と記載される。
【０１４８】  全パラメータ依存性を表示せず、またＦＯ行列または速度勾配に対する成分関数Ｆｉｊ（２）（．．．）の明示的依存性および黙示的依存性を区別しない、縮約記法（ｃｏｎｄｅｎｓｅｄ  ｎｏｔａｔｉｏｎ）を使用
2010531752.tif 000105
れる。添字記法（ｉｎｄｅｘ  ｎｏｔａｔｉｏｎ）から、式（１７ａ／ｂ）において使用される完全なテンソル記法に戻すことによって、これらの式に関する以下の単純化バージョンが得られる。
【０１４９】
2010531752.tif 000106
【０１５０】  ６．ＦＴＥの数値積分：「演算子分割（ｏｐｅｒａｔｏｒ  ｓｐｌｉｔｔｉｎｇ）」
【０１５１】
2010531752.tif 000107
ついては、多種多様の手法が存在する。特に非線形右辺関数Ｆ（．．．）の場合にロバストで柔
2010531752.tif 000108
により決定される「演算子」とを、別々に扱うことによって進められる。数学的文献において、この手法は、「分別ステップの方法（ｍｅｔｈｏｄ  ｏｆ  ｆｒａｃｔｉｏｎａｌ  ｓｔｅｐｓ）」、「分割方法（ｓｐｌｉｔｔｉｎｇ  ｍｅｔｈｏｄ）」、または単純に「演算子分割（ｏｐｅｒａｔｏｒ  ｓｐｌｉｔｔｉｎｇ）」として知られている（［２３〜３７］参照。代替方法については、［２８〜３１］を参照）。「演算子分割」においては、
2010531752.tif 000109
のうちのいずれかをゼロに設定することによって、元の式から得られる２つの式
2010531752.tif 000110
等しく、第２の式は、一般に結合型および非線形である常微分方程式（ｏｒｄｉｎａｒｙ  ｄｉｆｆｅｒｅｎｔｉａｌ  ｅｑｕａｔｉｏｎ；ＯＤＥ）の系列を構成する。このタイプの手法は好適な実施形態において使用される。式（１７ａ／ｂ）の記法を使用し、分割方法を閉包ＦＴＥに適用することによって、以下の部分方程式が得られる。
【０１５２】
2010531752.tif 000111
【０１５３】
2010531752.tif 000112
【０１５４】  これらの方程式は、「演算子分割」方法のフレームワーク内における別々の部分問題として扱われる。式（１８ａ）によりモデル化される物理的効果は、流速に従う流体質量の純粋な対流輸送による、モールドキャビティの充填領域内におけるＦＯ統計の空間再分布である（ＦＯ行列の要素によってコードされる）。一方、式（１８ｂ）は、対流輸送の効果を完全に無視し、局所速度勾配および繊維間の相互作用によってもたらされる繊維の回転運動によるＦＯ分布の局所変化率を単に説明するだけである。完全式（１６）（または抽象的記法におけるその同等物（１７ａ／ｂ））の近似解を得るために、両方の部分問題の（数値）解法をどのように組み合わせるかを示す、いくつかのバリエーションについて、以下に説明する。
【０１５５】  ６．１  「単純演算子分割（Ｓｉｍｐｌｅ  Ｏｐｅｒａｔｏｒ  Ｓｐｌｉｔｔｉｎｇ）」
【０１５６】  ２つの部分方程式から実際に解かれる方程式の近似解を得るための、最も単純な方式（通常、「単純演算子分割」と示す）を用いる場合、第１の式の中間解を、第２の式の入力（すなわち初期値）として解釈することにより、第１の式に続けて第２の式を解く。これを詳
2010531752.tif 000113
【０１５７】  フローソルバー（ｆｌｏｗ  ｓｏｌｖｅｒ）は、懸濁液のフローをモデル化するソフトウェアである。これは、モールドキャビティからゲートシステム、入口に及ぶ計算領域全体のうち、充填部分Ω（ｎ）に含まれる空間の離散点ｒｍの周囲に位置する全ての計算セルにおける、離散時間ｔｎにおける質量・運動量・エネルギーの離散化輸送方程式の解によって、流体の状態変数、
2010531752.tif 000114
2010531752.tif 000115
計算は、先行する計算ステップによって既に把握している「古い」領域Ω（ｎ）において定義
2010531752.tif 000116
【０１５８】  これは、「単純演算子分割」の利用可能なバリエーションの１つを構成する以下の３つのステップの手順によって達成される。
2010531752.tif 000117
る。
【０１５９】  図１１は、演算子分割の上の変形例を示す。
【０１６０】
2010531752.tif 000118
従って、好適な実施形態に従う方法において使用される。ただし、ステップサイズによっ
2010531752.tif 000119
する詳細は、後続の節において説明される。
【０１６１】  「単純演算子分割」の別の変形例
2010531752.tif 000120
ップによって出発し、次いで、
2010531752.tif 000121
【０１６２】  この変形では、対流ステップと回転ステップとが逆の順番で実行される。理論的な時間離散化誤差は、上述の第１の変形例と同様に同一次数を有するが、この代替変形例は、異なる瞬間時からの流速および速度勾配を使用するため、一貫性が低い傾向にある。これは、次のの別の変形によって回避できる
2010531752.tif 000122
（ｉｉ）  回転ステップを出発し、その後、
2010531752.tif 000123
う。
【０１６３】  ６．２  「対称的演算子分割（Ｓｙｍｍｅｔｒｉｃ  Ｏｐｅｒａｔｏｒ  Ｓｐｌｉｔｔｉｎｇ）」
【０１６４】  理論的には、時間離散化誤差に関する高次の精度は、「対称的演算子分割」手法によって得られる。この手法の基本概念は、ステップサイズが普通である部分方程式のうちの１つに関する１つの中間ステップを、ステップサイズが半分である他の方程式のステップの中に、括ってしまうことである。この手法に関して可能な１つの変形は、以下の４ステップによって行われる。
2010531752.tif 000124
【０１６５】
2010531752.tif 000125
る
【０１６６】  理論的に可能な別の「対称的分割」による変形例は、ステップサイズが半分である２つの対流ステップの間に回転ステップを一括する。この変形は、
2010531752.tif 000126
流前」ステップに進み、次いで、
2010531752.tif 000127
実行し、
2010531752.tif 000128
流後」ステップを完了させる。
【０１６７】
2010531752.tif 000129
純分割」の場合と同様に、初期初期化ステップ（拡張ステップ）および回転前回転ステップを置き換える第３の変形がさらに存在し、後者から開始することによって、Ω（ｎ）におけ
2010531752.tif 000130
【０１６８】  モールド充填に見られる典型的なフロー状態は、せん断フローを頻繁に伴うが、それは繊維の顕著な回転をもたらす。このため「回転演算子」部分は、純粋な対流輸送の効果に比べて、常に重要な（そして多くの場合、最も主要な）役割を果たす。ゆえに、対称的分割変形のうちの第１の変形に関して詳述するような、回転前ステップおよび回転後ステップによる対流ステップの対称的な括りは、「対称的分割」手法［１９］の好適な変形例を構成する。しかしながら、ほとんどの最も実際的な場合において、流体フロー計算により提供される時間ステップのサイズは十分小さいので、６．１節に示される「単純分割」による変形例は、多くの場合、十分正確に機能する。
【０１６９】  ６．３  別の分割変形
【０１７０】  右辺の関数Ｆ（．．．）の代数構造に依存して、多種多様の更なる分割手法を実現することが
2010531752.tif 000131
2010531752.tif 000132
は全部）を残すかもしれない。
【０１７１】  式（１６）の右辺を考察すると、３つの項の和が認識され、第１の項および第３の項が線形であり、中央の項がＦＯ行列に対して非線形依存性を恐らく有する閉包近似を伴う（閉包の選択に依存して、この中央の項は、下位の項の別の和に分割されうる）。上述のように、式（１６）の右辺の関数の項構造は、多種多様の追加の分割スキームを形成するための可能性を開く。式（１８ａ／ｂ）への分割に基づく変形例とは対照的に、これらの更なるの代替法は、相互に分離可能な物理的効果（例えば、「対流輸送」と局所「回転力学」の場合）の明確な区別に対応しない。ゆえに、これらは、単に「人工的数学的分離」と見なされるだけであり、物理的ＦＯ現象の適切なシミュレーションについては有用でない。これらの「人工的」分離変形の可能性は、完全性のためにのみ本明細書において記述され、本明細書においてさらに論じられない。
【０１７２】  ６．４  「対流ステップ」に対する数値解析的アプローチ
【０１７３】  本節では、閉包ＦＴＥに対する「演算子分割」手法内における部分問題として解決しなければならない対流式（１８ａ）の積分に使用する、数値解析的方法について説明する。新しい
2010531752.tif 000133
上法（ｆｉｒｓｔ  ｏｒｄｅｒ  ｕｐｗｉｎｄ  ｓｃｈｅｍｅ）［３２］を使用して組み立てられ、結果として得られる数式系
2010531752.tif 000134
【０１７４】  ６．５  「初期化ステップ」に対する数値解析的アプローチ
【０１７５】
2010531752.tif 000135
な実施形態に従う方法に適用されるアルゴリズムの一部として、修正されたＶＯＦ方法によって実行される。これは、フロー方程式（すなわち、質量・運動量・エネルギーに関する保存方程式）の数値解析のために行われる。ＦＯ行列の値は、両方の領域で共通するセルに
2010531752.tif 000136
填されたセル」を含むので、ＦＯ行列の要素の合理的な値の初期割当が必要になる。好適な実施形態によると、これは、先行して計算された質量フローに従って、その近傍セルからの又は近傍セルへの計算セル内における質量輸送に対応する加重平均によって達成される。
【０１７６】  粒子濃度が均一であると仮定され、かつＦＯ行列が巨視的なレベルにおける体積平均手順からもたらされるため、セルが、その近傍セルから、その中に輸送された流体質量の正味量に比例して、ＦＯへの寄与を得ることが見込まれる。ＦＴＥの位相空間（すなわち可能なＦＯ行列の組）について上に概説した数学的特徴によれば、後者は、（５次元）Ｔｒ［ａｉｊ（２）］＝１の超平面に制限される、実対称３×３行列の６次元ベクトル空間の有界かつ凸の部分集合である。質量輸送に従って重み付けされる平均は、ＦＯ行列の凸の組み合わせ（ｃｏｎｖｅｘ  ｃｏｍｂｉｎａｔｉｏｎ）を構成し、有効なＦＯ行列を常にもたらすことから、初期化の手順は、この手順の重要な特性であるＦＴモデルの位相空間の（理論的に必要とされる）トポロジカルな構造に互換性がある。
【０１７７】  ６．６  「噴フロー（Ｆｏｕｎｔａｉｎ  Ｆｌｏｗ）」効果の取り扱い
【０１７８】  用語の「噴フロー（Ｆｏｕｎｔａｉｎ  Ｆｌｏｗ）」は、多くの粘性非ニュートン流体の場合、流頭における自由表面付近の巨視的フローパターン全体を特徴付ける。「噴フロー」では、流頭付近の上流の粒子は、中心領域（ｃｏｒｅ  ｒｅｇｉｏｎ）から壁境界（ｗａｌｌ  ｂｏｅｎｄａｒｙ）へと移動する。この効果は、流体の材料特性により、実際に「自律的に」発生し、如何なる追加のモデリングをも必要としない。「噴フロー」は、適切な非ニュートン構成則を含むＮａｖｉｅｒ−Ｓｔｏｋｅｓ方程式に基づく３Ｄフロー計算においては、少なくとも突発的な現象ではある。しかしながら、これは、フロー計算が、単純化モデル（例えば、［１］におけるＣｒｏｃｈｅｔの論文の節１１．３において論じられるＨｅｌｅ−Ｓｈａｗ型のモデル）に基づく場合には当てはまらない。好適な実施形態に従う方法によるチャネルフローに関するシミュレーションでは、「噴フロー」パターンは、粒子トレースにおいて明確に認識可能であり、この現象をフローソルバーで取り扱うことができることが示される。
【０１７９】  自由表面に垂直であるＦＯ成分がゼロになるようにするために、いくつかの特別な措置が必要である。これは、繊維配向の「噴フロー」について観測される効果との一貫性の要請から必要である。またこれは、新しく充填されたセルの繊維配向の初期化手続において行われる。「必要とされる」垂直ＦＯ成分の消去は、次の理由により、上述の「初期化ステップ」において（意図的には）実行されない。数学的観点から、ＦＴＥモデルは、その右辺に結合される双曲輸送方程式系（ｓｙｓｔｅｍ  ｏｆ  ｈｙｐｅｒｂｏｌｉｃ  ｔｒａｎｓｐｏｒｔ  ｅｑｕａｔｉｏｎｓ）（すなわち、非線形反応項を有する対流−反応型系、上記参照）である。ゆえに、自由表面セルにおけるＦＯ行列の境界条件として、自由表面法線の方向におけるＦＯ行列の成分の消滅を規定することは、不適切である（すなわち、数学的に言うとできない）。このＦＯ成分の消滅は、実際は、フローについて計算された「噴フロー」挙動から自動的に生じるはずである。
【０１８０】  しかしながら、時間的および／または空間的離散化誤差、ならびにモールド充填中の自由表面の動きの計算に使用する手続によりもたらされる数値的誤差のため、恐らく、新しく充填されたセルの初期化手順が、「噴フロー」現象により必要とされる自由表面におけるＦＯ挙動と一致するように、いくつかの「訂正」を加えなければならない。全ての新しく充填されたセルについて、繊維配向テンソルのトレースが、理論値１との差が１％以下であることが確認される。差が大き過ぎる場合、繊維配向テンソルは、その固有空間における固有ベクトルの許容可能な３角形に対する直交射影によって訂正される。
【０１８１】  ６．７  シミュレーションの開始におけるＦＯ行列の初期化
【０１８２】  （閉包）ＦＴＥの双曲線構造は、モールド充填シミュレーションの時間ステップ毎に、境界条件として、入口における（ならびに入口付近の）計算セルにおけるＦＯ行列の初期化を要求する。ゆえに、この初期ＦＯ状態について適切な選択を行なわなければならない。好適な実施形態に従う方法では、後述のように、この目的のために、（ランダムＦＯ分布
2010531752.tif 000137
【０１８３】  高度に粘性のせん断フローでは、ＦＯ状態は、局所的準平衡のその「最終」状態に迅速に押し進められるため、高せん断レートの影響が強い。ゆえに、インゲート（溶融液が実際にその部分に入る場所）において観測されるＦＯ状態は、ランナーシステムにおけるそのフロー履歴によって完全に決定され、入口において見込まれるその初期状態からは相当に独立している。一方、ＦＴＥの解析的構造ならびにその運動学（すなわち位相空間における可能な状態）の分析によれば、ランダムな配向状態が速度勾配の全成分の全ての結合をもたらすことから、ランナーシステムの入口付近におけるフロー場によって後で適合されることを考慮すれば、入口における配向がランダムであるとの仮定は、最適な選択であることが示される。
【０１８４】  上述の演算子分割プロセスは、図１１を参照して説明される、単純化された形式である。
ステップ１：  新しく充填されたセルに関して初期条件および境界条件を設定し、
ステップ２：  フロー場に従って繊維配向を移動し、
ステップ３：  速度勾配に従って繊維回転を計算する。
【０１８５】  ７．  「回転ステップ」ＯＤＥ系に対する数値解析的アプローチ
【０１８６】  閉包ＦＴＥに対する「演算子分割」法の一部として解かなければならない、ＯＤＥの結合系（１８ｂ）を積分するには、数値解析的方法が使用される。この方法は、具体的には、好適な実施形態に従う方法のＦＯモジュールの実装のために開発されている。この方法は、短繊維強化熱可塑性プラスチック材料の３Ｄ射出成形シミュレーションの典型的な工業用途の枠の中において関連する種々の側面を含んでいる。例えば次の通りである。
・  ＦＯ行列の６つの独立要素の完全集合を使用するＦＴＥを定式化すること。（対角要素のうちの１つを無視し、ＦＴＥの従属変数を１つ減らすためにトレース条件を使用する「標準手順」とは異なる。）
・  トレース条件により定義された超平面（ｈｙｐｅｒｐｌａｎｅ）を、ＦＴＥの安定積分多様体（ｓｔａｂｌｅ  ｉｎｔｅｇｒａｌ  ｍａｎｉｆｏｌｄ）に変換するＦＴＥの右辺へ、「ペナルティ」制御項を追加すること。
・  局所せん断速度のサイズ（速度勾配の最大成分）に従って時間積分スキームを選択する積分法を構築するために、ＦＴＥの右辺関数の特定のスケーリング挙動を、速度勾配の要素に関して利用すること。
・  区間［０，１］に制限されるスカラー配向因子を含む「安定化ハイブリッド閉包」（ｓｔａｂｉｌｚｅｄ  ｈｙｂｒｉｄ  ｃｌｏｓｕｒｅ）を実装すること。
・  ハイブリッド閉包を有するＦＴＥの特別な右辺関数の評価のための効果的な手順を、最小限の演算数で実装すること。
これらの側面について以下に詳述する。
【０１８７】  ７．１  ハイブリッド閉包の安定化
【０１８８】  式（６）および（６ａ／ｂ）によって定義されるハイブリッド閉包近似は、特に、速度勾配が複雑な構造を有し（すなわち、単純なせん断／伸縮型ではなく、複雑な３Ｄフローパターンを反映する）、かつ時間ステップのサイズ（フローソルバーにより決定される）が比較的大きい場合、安全性問題に直面する。これらの不安定性の主な原因は、スカラー配向因子
2010531752.tif 000138
定になり、ＦＯ行列について許容可能に設定される位相空間集合から大きく離れた値に対して指数関数的に発散する。
【０１８９】
2010531752.tif 000139
る。
【０１９０】  ＦＯ行列の行列式は、常に区間［０，１／２７］に制限されるため、式（６ｃ）により定義される配
2010531752.tif 000140
ッド閉包に関する標準定義（６）は、以下のように修正される。
【０１９１】
2010531752.tif 000141
【０１９２】  式（６ａ／ｂ）により与えられる線形項および２次閉包項に関する定義は不変であり、式（６ｃ）
2010531752.tif 000142
値に評価され、そうでない場合は最小値または最大値として０または１をとる場合に、
2010531752.tif 000143
【０１９３】
2010531752.tif 000144
の区間［０，１］への制限が、考えられる試験例における深刻な不安定性の発達を阻止することが示された。安定化ハイブリッド閉包は、多種多様の３Ｄ射出成形シミュレーションの検査において成功している。
【０１９４】  ７．２  トレース条件を介する変数の数の減少
【０１９５】
2010531752.tif 000145
３×３行列である。その要素ａｉｊ（２）＝ａｊｉ（２）は、従属変数において相互に６つの集合を先験的に構成するため、ＦＴＥは、６つのＰＤＥの結合系である。ＦＴＥの代数構造は、任意の対称３×３行列を形式的に認めるが、その行列がＦＯ行列として認められるためには、その行列の不変量における制限として定式化される多数の追加の条件を満たさなければならない。
【０１９６】
2010531752.tif 000146
列の対角要素ａｋｋ（２）のうちの１つを無視（解消）しうる明らかな可能性をもたらし、変数の数を１つ減少させる。学問的な多くのフロー状況（例えば、単純なせん断フローまたは伸縮フロー）の下における繊維懸濁液フローを研究する大部分の（全部ではない）公表された研究論文は、解消される特定の対角要素について固定的選択を行なうことによって、
2010531752.tif 000147
が、式（１４）または（１５ａ）によって、式上、全ての合理的な閉包近似のトレース条件に適
2010531752.tif 000148
式のように）満たされ、ＦＯ行列の残りの集合の５つの成分およびさらにそのＰＤＥの結合系のみを考慮することで十分である。
【０１９７】  この手順は、純粋に形式的（ｐｕｒｅｌｙ  ｆｏｒｍａｌ）である。このように解消される対角要素の選択が良い選択であるか否かは、フローの局所速度場Ｕの空間特性およびその勾配
2010531752.tif 000149
ーでは、速度勾配の特別な形式が直交ｘ３方向に対して「弱い結合」のみをもたらし、繊維
2010531752.tif 000150
の安定性に関する問題）が生じるため、満足できない結果に至ることも考えられる。しか
2010531752.tif 000151
ために、任意の具体的な単純形式を有することを仮定できないことから、複雑なモールドの幾何学的形状内における３Ｄ射出成形シミュレーションでは、状況は大きく異なっている。そして、（上述のような）固定的選択が、適切な選択ではない可能性が高くなる。ゆえに、トレース条件によってＦＯ行列の対角要素のうちの１つを解消することは、工業用懸濁液フローのシミュレーションには不適切であると考えられる。
【０１９８】  単位トレース条件により対角要素のうちの１つを解消することに反対する、（数値的安定性に関する）別の議論が、Ｓｈａｍｐｉｎｅによる一連の論文［３４〜３６］に見られる。この論文では、特別の代数的制約（「保存則」）を有するＯＤＥ系について研究されている。［３４〜３６］
2010531752.tif 000152
辺関数Ｆ（．．．）に代入することによって、数式から濃度のうちの１つを解消する可能性をもたらし、それによって同時に保存則を厳密に満たす。が、一方、Ｓｈａｍｐｉｎｅは、数値積分手順内におけるこの「トリック」を使用することが、単に、保存則を介して代数的に得られた残りの濃度ｃｎ内において（減少した）ＯＤＥ系の数値積分によって計算される濃度
2010531752.tif 000153
常に厳密に満たされるが、このように計算される数値解は、正確であることが決して保証されない。実際、特に硬い（ｓｔｉｆｆ）ＯＤＥ系では、この「直接的な解消手法」が不注意に使用される場合に、数値解が、真の解から逸脱し得ることが知られている。［３４〜３６］による
2010531752.tif 000154
非線形制約が考慮される場合であっても、同一の結論をもたらす。
【０１９９】  線形な「保存則」である「単位トレース」制約により拡張されるＦｏｌｇａｒ−Ｔｕｃｋｅｒ方程式は、［３４〜３６］に提示される数学的方程式の特定の例であるため、Ｓｈａｍｐｉｎｅの数学的分析は、単位トレース条件により対角要素のうちの１つを解消することによって数値的不安定がもたらされる可能性が高いことを示し、ＦＯ行列の対角要素における制御不可能な不正確性がもたらされることを示す。
【０２００】  好適な実施形態に従う方法は、ＦＯの計算に、ＦＯ行列の６つの要素の全てと、ＦＴＥ系の方程式とを使用する。上述の議論のように、これは、単位トレース条件を介して対角要素のうちの１つを解消する「標準的手法」からは逸脱した方法である。
【０２０１】  ７．３  ペナルティ項（ｐｅｎａｌｔｙ  ｔｅｒｍ）による動的なトレース安定化
【０２０２】  閉包ＦＴＥ系およびその「回転ステップ」部分が、式（１６）および（１８ｂ）により与えられると考えられ、かつトレース条件が、従属変数の集合からＦＯ行列の対角要素のうちの１つを解消するために使用されない場合、数値的に計算されたＦＯ行列のトレースを、その望ましい値１に可能な限り近づけて維持するためには、いくつかの追加の方法が必要である。
【０２０３】  論文［３４〜３６］によれば、（前の７．２節に論じられたように、）このような解消は、いかなる種類の手法であっても避けることが推奨され、代わりに、通常のＯＤＥ積分方法と、所要の保存則または代数的制約を（正確にまたは近似的に）満たす、何らかの射影手続（ｐｒｏｊｅｃｔｉｏｎ  ｐｒｏｃｅｄｕｒｅ）との組み合わせによって、完全な系を取り扱うことが推奨される。個々の時間ステップにおいて、まずは代数的制約を考慮せずに数値的解が計算され、その後、制約方程式によって定義される代数多様体（例えば、線形制約の場合、超平面）の最近接点への射影マッピング（ｐｒｏｊｅｃｔｉｏｎ  ｍａｐｐｉｎｇ）により補正される。
【０２０４】  この手法の別の方法は、適切な制御項をモデル方程式の右辺に追加することによって、制約を考慮することを試みており、数値的に計算されたＦＯ行列のトレースを、その所望値１に可能な限り近似に維持することを試みている。制御項の存在によって、更なる追加の測度（ｍｅａｓｕｒｅ）を含まずに、制約または保存則が近似的に満たされる。多くの場合、「ハード」な制御が、ある種の射影を効果的にもたらし、一方、「ソフト」な制御が、所要の制約からの逸脱を小さくすることを可能にすることを示すことができる。
【０２０５】  好適な実施形態に従う方法において、このような制御法は、特に、数値的に計算されたＦＯ行列のトレースを、その所要単位値に近似的に維持するために使用される。安定化は、適切な制御項を式（１６）または（１８ｂ）の右辺関数に加えることによって、トレース条件
2010531752.tif 000155
の安定した定点になる状況に等しい。
【０２０６】  合理的な制御項により満たされるべき一般的な必要性は、以下であるように定められている。
・  トレース条件が既に満たされている場合、制御項は全く無いものと同様でなければならない。
2010531752.tif 000156
なければならない。
【０２０７】  これらの２つの主要な必要性の他に、合理的な制御項は、速度勾配の再スケーリングに関してＦＴＥの右辺関数のスケーリング挙動に適合するべきであり（この点に関するさらなる詳細については、７．４節のフロー制御時間積分を参照）、また、「ソフト」または「ハード」な制御の調節を可能にする調節パラメータを含むべきである。後者は、「無限にハード」な制御の制限において、実際上射影に類似する。
【０２０８】  ハイブリッド閉包近似（標準的かつ安定化バリエーションの両方を含む）という特定の
2010531752.tif 000157
て、トレースのＤＥは、以下の式を介して得られる。
2010531752.tif 000158
【０２０９】  （標準的な）ハイブリッド閉包に関する特別な場合では、右辺関数のトレースは、
2010531752.tif 000159
れる変数前因子φａを有するトレースのＤＥの単純な形式（１５ａ）をもたらす。標準的なハ
2010531752.tif 000160
因子φａを有するトレースからも、完全に同一のＤＥが得られる。
【０２１０】
2010531752.tif 000161
ナルティ項の或る選択肢は、次の通りである。
【０２１１】
2010531752.tif 000162
【０２１２】  パラメータαは正であることが必要とされるが、時間に対して必ずしも一定である必要
2010531752.tif 000163
定化ハイブリッド閉包を有する閉包ＦＴＥの、一般化された次の変形例も考慮しなければならない。
【０２１３】
2010531752.tif 000164
【０２１４】
2010531752.tif 000165
【０２１５】
2010531752.tif 000166
に従って、式（１５ａ／ｂ）を修正するように調節される。
【０２１６】
2010531752.tif 000167
面は、ハイブリッド閉包を有するＦＴＥの修正された（「コントロールされた」）式（２１）の安定積分多様体（ｓｔａｂｌｅ  ｉｎｔｅｇｒａｌ  ｍａｎｉｆｏｌｄ）になる。結果として、式（２１）の任意の数値解のトレースは、制御項（２０）により所要値１に近似的に維持されるが、それは、この項が、
2010531752.tif 000168
し進められるためである。
【０２１７】  制御項（２０）の「強度」は、パラメータαを適切に調節することによって調整できる。つまり、パラメータαの値が小さいと、トレース値を１の方向へ比較的「ソフト」に訂正
2010531752.tif 000169
影（ｉｎｓｔｒａｎｔ  ｐｒｏｊｅｃｔｉｏｎ）の効果に近い）。トレース訂正が大きい場合、超平面の両側に数値解を交互に押し出すことによるスプリアス振動（ｓｐｕｒｉｏｕｓ  ｏｓｃｉｌｌａｔｉｏｎ）が生じるため、実際は、αの値が大きいと、明示的ＯＤＥ積分方法と共に使用される際に問題が生じる場合がある。しかしながら、中間的な値においては、比較的広範囲のαの選択が、制御項によりもたらされる不要なアーチファクトを含まない、適切な数値的結果をもたらすことが示されている［１８］。
【０２１８】
2010531752.tif 000170
制御項の以下の一般的定義が、トレースＤＥ（２２）を一意（ｉｄｅｎｔｉｃａｌｌｙ）にもたらし、（安
2010531752.tif 000171
たらされる。
【０２１９】
2010531752.tif 000172
【０２２０】
2010531752.tif 000173
は、制御項の最も一般的な形式を構成する。
【０２２１】  ７．４  フローが制御された時間積分（Ｆｌｏｗ  ｃｏｎｔｒｏｌｌｅｄ  ｔｉｍｅ  ｉｎｔｅｇｒａｔｉｏｎ）
【０２２２】  「演算子分割」手法を使用することによって、「回転ステップ」ＯＤＥ系（１８ｂ）の数値
2010531752.tif 000174
て（１８ｂ）の右辺に入る。これらが大きい場合、ＦＯ行列要素の値が、時間間隔の間に大きく変動する。一方、速度勾配が非常に小さい場合、先行する計算ステップにより提供される、ＦＯ行列要素の初期値からの変動は、ごく僅かであるだろうということが予測できる。
【０２２３】  射出成形シミュレーションに関する全ての実際的な場合では、モールドキャビティ内の充填領域における速度勾配は、充填中に大きく変動する。せん断速度の値は、典型的にはキャビティ壁付近で高く、壁に挟まれる間隙の中心領域においては、かなり低い値が観測されている。これにより、流速方向において、壁付近の強いせん断フローの影響のため、繊維の配向状態はかなり迅速に変化するが、中心領域の繊維は、その配向状態を変化させるのに著しく長い時間を要する。このため、周知の「サンドイッチ構造」（壁付近の高度に配向した領域の間に、配向の度合いが比較的低い「コア」中心領域を有する）がもたらされる。これは、典型的には、短繊維強化熱可塑性プラスチックから作製される射出成形部分において観測される。
【０２２４】  ＯＤＥ系の積分の数値スキームは、ステップサイズを適切に（好ましくは適応的に）選択することによって、右辺関数の「強度（ｓｔｒｅｎｇｔｈ）」を勘案する。つまり、値の間隔が大きい粗いステップを用いることや、右辺関数を大きく変化させること、また、右辺関数の（絶対）値がかなり小さい領域において小さ過ぎるステップを取ることにより生じる、不正確性を回避することが試みられる。ステップサイズや積分スキームの適切な選択によって、ＯＤＥ積分の数値誤差を、所望の限界まで制御することが可能になる。
【０２２５】
2010531752.tif 000175
が、これは、速度勾配の変動による右辺関数の変動にもかかわらず、計算領域の全セルにおいて、同じように提供されねばならない。この目的のために、右辺関数の評価の回数を最小限に抑えつつこの目標を達成する、特別な数値スキームが考案されている。このスキームでは、
・  「積分法」（オイラー、中間点、または４次ルンゲ・クッタ（Ｒｕｎｇｅ  Ｋｕｔｔａ））
・  （無次元）局所ステップサイズ
・  選択された法則とサイズによる、ステップの局所的な数
を選択する。これらは、局所速度勾配の値によって決定される局所右辺関数の「強度（Ｓｔｒｅｎｇｔｈ）」に適応させられる。
【０２２６】  このスキームは、数値的ＯＤＥ積分の「適応ステップサイズ制御」に関する既知の手法とは異なる。それは、速度勾配の成分に関する式（１６）または（１８ｂ）の右辺関数のスケーリング挙動（ｓｃａｌｉｎｇ  ｂｅｈａｖｉｏｕｒ）を活用し、ＦＴＥ系に具体的に適応させるためである。積分法ならびにスキームに伴う局所時間ステップの数およびサイズの選択に関する手続は、使用する積分法の理論的離散化誤差に基づく発見的手法である。スキームの両側面について、以下に詳述する。
【０２２７】
2010531752.tif 000176
として、式（１６）から得られてもよい。
【０２２８】
2010531752.tif 000177
よびせん断速度テンソルに関する定義（１ａ／ｂ）、すなわち、
2010531752.tif 000178
の右辺に関して既に実行されているように、）右辺関数の代数的表現を、以下の単純化形式で記述することが役立つ。
【０２２９】
2010531752.tif 000179
【０２３０】
2010531752.tif 000180
る依存性のみは残す。右辺関数（２３）を使用するＯＤＥ系（１８ｂ）の（（１７ｃ）に対応する）同等な再定式化は、以下の形式で記述できる。
【０２３１】
2010531752.tif 000181
【０２３２】
2010531752.tif 000182
らす。
【０２３３】
2010531752.tif 000183
【０２３４】  前述のように、一般的にやらなければならないタスクは、時間間隔［ｔｎ，ｔｎ＋ｌ］におけるＯＤＥの結合系（１８ｃ）の数値積分である。この時間間隔において、グローバル数値誤差εＦＯが存在するが、これは充填セル領域の全セルにおいて均一である。このタスクは、速度勾配の
2010531752.tif 000184
利用することによって、以下のように達成できる。
・  最大の絶対値を有する速度勾配の成分の値、つまり、
2010531752.tif 000185
局所計算セルを識別する空間ベクトルｒｍや、速度勾配が右辺関数に入る時間を表す時間座
2010531752.tif 000186
された量および式（２４）を使用して、ＯＤＥ（１８ｃ）を、以下のスケーリングされた形式に変換する。
【０２３５】
2010531752.tif 000187
【０２３６】
2010531752.tif 000188
領域内の各セルは、それぞれ異なるパラメータγｍａｘ．を示す。
【０２３７】
2010531752.tif 000189
に亘り、次の無次元ステップサイズで積分されなければならない。
2010531752.tif 000190
【０２３８】
2010531752.tif 000191
以下の因子により、ほぼ均一の「単位強度（ｕｎｉｔ  ｓｔｒｅｎｇｔｈ）」を有する。
2010531752.tif 000192
に等しいか、１未満である無次元量である。
2010531752.tif 000193
素の値も、１によって制限される。
2010531752.tif 000194
値を有する。
【０２３９】
2010531752.tif 000195
述は、以下の正確な数学的形式に定式化される。
【０２４０】
2010531752.tif 000196
【０２４１】
2010531752.tif 000197
において、
ｉ．    十分小さい離散化誤差を有する「積分法」、ならびに
ｉｉ．  スケーリングされた時間間隔をカバーする適切な数のサブステップ
を選択することによって、達成できる。
【０２４２】  好適な実施形態に従う方法において使用される積分スキームは、ワン・ステップ・メソッド（ｏｎｅ  ｓｔｅｐ  ｍｅｔｈｏｄ）に分類される単純な積分法の組を使用する。（この積分法については［３８］の節１６．１および［３９］の７．２．１〜７．２．３を参照。これらは参照により本明細書に組み込まれる。）好適な実施形態に従う方法において、ＦＯモジュール内で使用される特別の積分法は次の通りである。
・  ステップ毎に右辺関数の評価を必要とする単純な前進オイラー法（ｆｏｒｗａｒｄ  Ｅｕｌｅｒ  ｒｕｌｅ）。１次の精度を有する。
・  ステップ毎に右辺関数に関する２つの評価を必要とする「中間点（ｍｉｄｏｐｏｉｎｔ）」または２次ルンゲ・クッタ（Ｒｕｎｇｅ−Ｋｕｔｔａ；ＲＫ２）法。２次の精度を有する。
・  ステップ毎に右辺関数に関する４つの評価を必要とする４次ルンゲ・クッタ（Ｒｕｎｇｅ−Ｋｕｔｔａ；ＲＫ４）法。４次の精度を有する。
【０２４３】  積分法に関する以下の特性が、本発明のある実施形態に使用される特別なスキームの構造において使用されるが、これらは次数（ｏｒｄｅｒ）ｐに依存する。
（ａ）次数ｐである方法は、サブステップ毎に右辺関数に関するｐ回の評価を必要とする。
（ｂ）ＯＤＥ（系）を、全サイズΔτについて、ｈ＝Δτ／ＮのＮ個の等距離サブステップで数値積分する場合、積分の最終サブステップにおいて累積される全誤差εｔｏｔは、
2010531752.tif 000198
（ｃ）この全誤差が、事前に選択された閾値ｆεｍａｘよりも小さいことが必要とされる場合、
2010531752.tif 000199
（ｄ）積分を終えたときに、全誤差をεｍａｘよりも小さくするには、積分間隔のサイズが
2010531752.tif 000200
ブステップ（すなわち、Ｎ＝１）で十分であることを示唆する。
【０２４４】  これを考慮しすると、
（ｉ）３つの積分法のうちの１つによる単一ステップ、または、
（ｉｉ）ＲＫ４法を使用する適切なサイズの多数のサブステップ、
を選択することによって、全体誤差がεｍａｘ＝εＦＯより小さくなる積分スキームを得ること
2010531752.tif 000201
された時間間隔で式（２５）の数値積分を行う、「ハイブリッド（ｈｙｂｒｉｄ）」積分スキーム
2010531752.tif 000202
に関わらず、均一である。
【０２４５】  提案された方法は、具体的には、上記論点（ａ）から（ｄ）に提示される関係において推定されるように、所要の誤差境界を達成するために必要な右辺関数評価の数が最小である
2010531752.tif 000203
マー溶融液は間隙サイズが小さく高粘度であるので、γｍａｘ、の値はかなり大きくなり得る（約
2010531752.tif 000204
更新された値も、Ｏ（１）のオーダーとなりうるため、所要の誤差境界の合理的な選択は、
2010531752.tif 000205
よるＲＫ４法は、半分のサイズｈ／２の２つのステップによる「中間点」（ＲＫ２）法や、４分の１のサイズｈ／４の４つのステップによるオイラー法よりも、正確であることを示すことができる。しかしが、数値計算にかかる労力（右辺関数の４つの評価）は、全てのバリエーションについて同一である。
【０２４６】
2010531752.tif 000206
する値ｐ∈｛１，２，４｝について計算される。０＜εＦＯ＜１であるので（実際は、通常εＦＯ≪１である）、誤差境界のサイズは、オーダーパラメータｐの値が増加するにつれて、単調に増加する（常
2010531752.tif 000207
視できるほど小さな差しか与えない。
これは、典型的には、速度勾配が小さい領域のいずれのセルにおいても発生する。例えば、２つの隣接する壁の間に比較的大きい間隙が存在する場合、中心中心領域における速度勾配は非常に小さい。例えば、εｍｉｍ＝１０〜６の最小閾値は、典型的なモールド充填用途において、合理的な選択である。
【０２４７】  上述の全ての側面を考慮すると、領域の各セル内において実行するためのアルゴリズムは、以下のように定式化され得る。なお、このアルゴリズムは図１２のフローチャートにも示される。
【０２４８】
2010531752.tif 000208
【０２４９】
2010531752.tif 000209
（ａ）そうである場合、ＯＤＥ積分を飛ばし（すなわち何も実行しない）、セルにおけるＦＯ行列の先行値を、新しい（更新された）値として維持し、手順を終了する。
（ｂ）そうでない場合、ステップ１２．２に進む。
【０２５０】
2010531752.tif 000210
手順を終了する。
（ｂ）そうでない場合、ステップ１２．３．に進む。
【０２５１】
2010531752.tif 000211
ステップによって、セルにおけるＦＯ行列の先行値を更新し、手順を終了する。
（ｂ）そうでない場合、ステップ１２．４に進む。
【０２５２】
2010531752.tif 000212
クッタ」（ＲＫ４）ステップによって、セルにおけるＦＯ行列の先行値を更新し、手順を終了する。
（ｂ）そうでない場合、ステップ１２．５に進む。
【０２５３】
2010531752.tif 000213
有するいくつかのＲＫ４ステップを実行することによって計算される。
（ａ）所要のεＦＯ精度を達成するための最小のステップ数は、不等式
2010531752.tif 000214
数の整数部分を入手し、ｍａｘ（．．．，．．．）を使用して２つの数字のうちの大きい方を入手する
2010531752.tif 000215
2010531752.tif 000216
（ｃ）Ｎおよびｈが計算されると、ステップサイズｈを有するＲＫ４法のＮ回のステッ
2010531752.tif 000217
し、手順を終了する。
【０２５４】
2010531752.tif 000218
ッド」積分法を定義する。この積分法は、誤差が全体的に均一でεＦＯより小さく、式（１８ｂ）の右辺関数の評価の必要数が最小で済む積分法である。（１８ｂ）のスケーリングされたバージョン（２５）への当該積分法の適用は、右辺関数の評価において、スケーリングされた速
2010531752.tif 000219
イズｈによる計算を使用することによって、「同じように（ｉｎ  ｐｌａｃｅ）」実行される。
【０２５５】  このアルゴリズムの典型的な用途においては、誤差境界パラメータεＦＯ＝０．００１およびεｍｉｍ＝１０−６を有する。この場合、領域Ωｎ＋ｌにおけるセルの大部分において（すなわち約８０％）、ＦＯ行列の「回転ステップ」の更新は、単一オイラーステップによって実行される。また比
2010531752.tif 000220
キップする。恐らくモールドキャビティの壁付近に位置するため、高いせん断速度値を有する適当な数のセルだけが、１回または複数回の（典型的には２０以下の）「４次ルンゲ・クッタ」ステップで更新されることが観測される。
【０２５６】  ７．５  「回転ステップ」ＯＤＥ積分における動的トレース安定化
【０２５７】  動的なトレース安定化（ｄｙｎａｍｉｃａｌ  ｔｒａｃｅ  ｓｔａｌｉｚａｔｉｏｎ）を達成すべく、右辺関数が、（ハイブリッド閉包の場合の特別な型（２０）または式（２０ａ）によって与えられる一般的な型のいずれかの）制御項を含む場合に、式（２０）又は（２０ａ）の制御パラメータがα＝α０・γｍａｘ（α０は典型的なサイズがＯ（１）の無次元制御パラメータ）と定義されるならば、右辺関数のスケーリング挙動は、依然として、本節に説明される「ハイブリッド」積分アルゴリズムに適合する。
【０２５８】  これを、制御項をその一般的な形式（２０ａ）で考察することによって示す。任意の閉包近
2010531752.tif 000221
て、そのトレースを書き直してみる。
2010531752.tif 000222
【０２５９】  また、制御因子α＝α０・γｍａｘを使用して、制御項（２０ａ）について以下の等化表現を得る。
2010531752.tif 000223
【０２６０】
2010531752.tif 000224
る。その他は（２０ａ）式と同一である。これは、以下の式によって表現され得る。
【０２６１】
2010531752.tif 000225
【０２６２】  この式はまた、その一般的な形式（２０ａ）の制御項が、制御パラメータがα＝α０・γｍａｘとして選択される場合に、γｍａｘによる速度勾配のスケーリングに対して、「制御されていな
2010531752.tif 000226
【０２６３】  （安定化）ハイブリッド閉包の特別な場合では、項（２０）を以下のように書き直すこと
2010531752.tif 000227
【０２６４】
2010531752.tif 000228
た制御項の特別な形式についても（２８）が有効であることを明示的に示す。
【０２６５】  上の考察は、ａ０が制御パラメータとして使用され、右辺関数に依存する項がスケーリングされた速度勾配で評価される場合に、「回転ステップ」ＯＤＥ系の「ハイブリッド」積分アルゴリズムが、制御項の存在下においても変更することなく適用できることを示す。式（２１）に対応するＯＤＥのスケーリングされたバージョンは、単に、制御項を式（２５）の右辺に追加することによって得られる。すなわち、以下のＯＤＥを、積分アルゴリズムの離散時間ステップの基礎として、（２５）の代わりに考慮しなければならない。
【０２６６】
2010531752.tif 000229
【０２６７】  好適な実施形態に従う方法のＦＯモジュールに実装される時間積分スキームは、本節に定義される「ハイブリッド」積分アルゴリズムを、制御項（２０）によって、動的トレース安定化を含むＯＤＥ系（２９）に適用することによって、「回転ステップ」ＯＤＥ系の数値積分を実行する。
【０２６８】  ７．６  安定化ハイブリッド閉包を含むＦＴＥの効率的な評価
【０２６９】  右辺関数の評価（ｅｖａｌｕａｔｉｏｎ）は、ＦＯ計算手続において最もコストのかかる部分であるため、本方法は、最も効率的なやり方で評価する。
【０２７０】  ７．４節において提案されるように、「フロー制御時間積分」（ｆｌｏｗ  ｃｏｎｔｒｏｌｌｅｄ  ｔｉｍｅ  ｉｎｔｅｇｒａｔｉｏｎ；ＦＣＴＩ）アルゴリズムによるアプローチによって、効率に関してある程度の対処がなされている。このアプローチを用いることにより、ＦＴＥの右辺関数（恐らく制御項を含む）の評価を最小限に止めつつ、均一な精度を達成可能である。
【０２７１】  ＦＯ計算の計算コストに直接影響を及ぼす最重要因子は、閉包近似の選択である。安定化ハイブリッド閉包は、低い計算コストで合理的な精度をもたらすので、本方法で使用される。
【０２７２】ハイブリッド閉包を有するＦＴＥの代数的再定式化
【０２７３】  効率改善における第１のステップは、ハイブリッド閉包の代数的定義（６）を、ＦＴＥの右辺関数（２３）に代入することである。これは、（数回の代数的変換後、）次の特有の代数形式を有する以下の右辺関数をもたらす。
【０２７４】
2010531752.tif 000230
【０２７５】
2010531752.tif 000231
【０２７６】
2010531752.tif 000232
【０２７７】
2010531752.tif 000233
【０２７８】
2010531752.tif 000234
【０２７９】
2010531752.tif 000235
子φａに関して上記に提供される式は、（１５ｂ）に示される式に等しい。
【０２８０】  公式（３２）は、圧縮性流速場または非圧縮性流速場のいずれの場合においても有効である。ここで、フローが非圧縮性であると仮定すると、理論的にｄｉｖＵ＝０が適用できるは
2010531752.tif 000236
において維持されることが強調される。数値誤差（これはフロー計算中に必然的に発生す
2010531752.tif 000237
に比例する項が、前因子公式（３２）においてゼロであると先験的に仮定される場合は、不安定性が発生することが分かる。
【０２８１】  量（３１）および（３２）を有する形式（３０）における右辺関数の使用は、効率的な評価に対する基礎的なステップである。というのも、それによって、右辺関数（３０）の代数構造が、多種多様の計算が一度だけ（前因子の計算時に）実行されなければならないように体系化されるため、多く演算を行わなくて済むからである。例えば（３１）および（３２）における項２Ｄｒまたは１−ｆｓ）／７など、頻繁に出現する共通の式があり、それらを一度だけ計算すること、および後に使用するためにその値を補助変数に格納することによって、計算作業が大幅に低減される。
【０２８２】  さらに、
【０２８３】
2010531752.tif 000238
ことによって、いくつかの演算を節約することができる。
【０２８４】  代数式
2010531752.tif 000239
【０２８５】
2010531752.tif 000240
に従って再定義される。
【０２８６】
2010531752.tif 000241
関数の評価のために式（３０ａ）および（３１ａ）を使用すると、（３０）および（３１）を使用する場合に比べて、演算の総回数を少なくすることができる。
【０２８７】
2010531752.tif 000242
【０２８８】  以下の制御項の追加（式（２０）参照）
2010531752.tif 000243
は、右辺の評価に必要なさらなる演算をもたらす。この制御項の一般構造は、前因子
2010531752.tif 000244
の場合は、以下の項
2010531752.tif 000245
をもたらす。
【０２８９】
2010531752.tif 000246
【０２９０】右辺評価のための効率的なアルゴリズム
【０２９１】  上記ステップを利用して、以下の一連の計算に従って、式（３０ｂ）の評価を最小回数の代数的演算で行うスキームが実現される。
【０２９２】  １．入力：
2010531752.tif 000247
【０２９３】
2010531752.tif 000248
【０２９４】
2010531752.tif 000249
【０２９５】
2010531752.tif 000250
算する。
【０２９６】
2010531752.tif 000251
ｖｅｒｓｉｏｎ）
2010531752.tif 000252
【０２９７】
2010531752.tif 000253
【０２９８】
2010531752.tif 000254
を計算する。
2010531752.tif 000255
【０２９９】
2010531752.tif 000256
【０３００】
2010531752.tif 000257
一連の演算によって実行する。
2010531752.tif 000258
【０３０１】
2010531752.tif 000259
引き続き、右辺関数（３０ｂ）の評価を完成させる。
2010531752.tif 000260
【０３０２】
2010531752.tif 000261
【０３０３】「縮約記法」（ｃｏｎｔｒａｃｔｅｄ  ｎｏｔａｔｉｏｎ；ＣＮ）を使用する、ベクトル形式におけるＦＴＥの再定式化
【０３０４】
2010531752.tif 000262
く。ゆえに、一般３×３行列としてそれらを格納することは非実際的であり、また、行列
2010531752.tif 000263
発明によると、対称３×３行列は６つの成分実ベクトルとして処理される。これは、ベク
2010531752.tif 000264
ば［２２］参照）（構造力学では、これは「フォークト記法」として既知である）として知られる対称行列の指標対（ｉｊ）＝（ｊｉ）を使用して行われる。
【０３０５】
2010531752.tif 000265
【０３０６】
2010531752.tif 000266
つだけである。上記表は、好適な実施形態のＦＯモジュールにおいても導入される通常の選択を示す（代替例は、通常、非対角行列要素のベクトル指標μ∈｛４，５，６｝への割り当てを選択する）。
【０３０７】
2010531752.tif 000267
れる。
【０３０８】
2010531752.tif 000268
【０３０９】
2010531752.tif 000269
成分は、第１および第２の指標対に対してテンソルが対称的である場合（すなわち、
2010531752.tif 000270
独立要素の数は、テンソル成分間のさらなる対称性および代数的関係の存在によって低減され得る。完全対称４次テンソルでは、以下の６つの追加の式が存在する。
2010531752.tif 000271
を、さらに小さい値に減少させる、代数的関係の組である（詳細については［２２］を参照）。
【０３１０】  ＣＮ法の導入により、以下のように、「ベクトル形式」によるＦＴＥがもたらされる。
【０３１１】
2010531752.tif 000272
【０３１２】
2010531752.tif 000273
（３５ａ）
【０３１３】
2010531752.tif 000274
である。
【０３１４】
2010531752.tif 000275
される。
2010531752.tif 000276
【０３１５】
2010531752.tif 000277
【０３１６】  （３５ｂ）
【０３１７】
2010531752.tif 000278
【０３１８】
2010531752.tif 000279
【０３１９】提案されたアルゴリズムの計算コスト
【０３２０】  右辺関数の効率的な評価について提案されたステップ２から１０の手続において実行される計算演算に関する最も効率的な実装は、ＣＮスキームをこの手順に適用することによってもたらされる。以下の表は、ＣＮ法を使用する手続の個々のステップの実行に必要な数
2010531752.tif 000280
概説を提供する。
2010531752.tif 000281
【０３２１】  以下、上記に提示されるアルゴリズムの種々のステップに関する表において提供される実装の側面ならびに演算数に関してコメントする。
【０３２２】
2010531752.tif 000282
【０３２３】
2010531752.tif 000283
最終ステップ１０において６ベクトルとしても現れる。
【０３２４】  ・  格納演算および割り当て演算は、計算コスト全体に対するそれらの寄与が無視できるため、勘定されていない。
【０３２５】  ・  補助変数の数は、その値が必要とされなくると、既存の変数を再利用することによって減少される。これは、上に提案されるアルゴリズムのバージョンにおいては説明されていないが、それは、それを説明すると、アルゴリズムの説明の明確性が悪くなり得るためである。
【０３２６】  ・  対称行列の対の縮約演算（ステップ４および６参照）は、ＣＮを使用すると、乗除算７回および加減算５回の演算で実装される（行列記法では、それは、
2010531752.tif 000284
【０３２７】  ・  実対称３×３行列の行列式（→ステップ５）は、ＣＮを使用して、乗除算１１回および加減算５回の演算で計算される。
【０３２８】
2010531752.tif 000285
およびξ６を計算するために、２回の乗算を追加する。
【０３２９】
2010531752.tif 000286
１２回の演算を必要とする。
【０３３０】
2010531752.tif 000287
のβの３回の加算のみを必要とする。
【０３３１】  表の最終列は、右辺関数の１回の評価が、９６回の乗除算、８４回の加減算、３回の関数呼び出しという計算コストによって得られることを示している。（この関数呼び出しは、（倍精度）浮動小数点数の平方根の計算と、このような数の対の最小値および最大値の計算するために行われる。）現代のコンピュータハードウェアでは、加算演算および乗算演算の計算コストはほぼ同一であり、ＭＩＮ関数またはＭＡＸ関数に対する呼び出しは、約１．５〜２個の演算のコストがかかり、平方根計算の総コストは、約６〜１０個の演算になる（すなわち、所要の精度に応じて平均８個の演算）。
【０３３２】  従って、提案されるアルゴリズムによれば、式（３０ｂ）の１回の評価は、全体で、約１９０個の演算の計算コストによって達成される。これは、トレース安定化を含む（安定化）ハイブリッド閉包を有するＦＴＥの右辺を評価するための手法としては、最も効率的である。
【０３３３】計算効率の評価
【０３３４】  提案された手順の計算効率に関する評価は、式（３５）によりベクトル形式で与えられる、右辺関数を評価する「標準的な」手法との比較により行うことができる。このような手法は、例えば、コードが特定のモジュール構造を有するように意図される場合に用いなければならないだろう。この特定のモジュール構造とは、（ｉ）閉包近似とは独立の全ての演算
2010531752.tif 000288
固有の選択に従って、別々のルーティンにおいて計算される、というものである。このような「標準的な」手法の計算コストは、以下の考察によって評価できる。
【０３３５】
2010531752.tif 000289
ラメータ２Ｄｒ（ステップ２から４参照）の計算から出発する。これらは合計４６回の演算と、追加の平方根計算が必要である。
【０３３６】
2010531752.tif 000290
行されると共に、６つの成分の結果として得られる式のみが数値的に評価されるる場合、４８回の演算を必要とする。
2010531752.tif 000291
算を必要とする。
・  ２Ｄｒ・１の加算および６Ｄｒ・ａの減算（６Ｄｒ＝３／２・２Ｄｒを使用する）には、別に３＋１２＋１＝１６個の演算を必要とする。
【０３３７】  平方根計算について８回の演算を勘定すると、この手法が式（３５）の右辺を評価するための計算コストは、ここまでで既に最大１９０回の演算となる。トレース安定化のために制御項を含めることによって、いくつかのさらなる演算を行わねばならない（ステップ８参
2010531752.tif 000292
３Ｄｒ＝３／２・２Ｄｒを減算する演算に置き換えることによって省かれ得る。
【０３３８】
2010531752.tif 000293
に必要とされる閉包近似とは別である。
【０３３９】
2010531752.tif 000294
計算に使用されるべき以下の組の式がもたらされる。
【０３４０】
2010531752.tif 000295
【０３４１】
2010531752.tif 000296
４次テンソル式の成分を含む。
【０３４２】
2010531752.tif 000297
【０３４３】
2010531752.tif 000298
【０３４４】
2010531752.tif 000299
たがって決定される。
【０３４５】
2010531752.tif 000300
ステップ５によって実行され、２３回の演算を必要とする（ネストされたｍｉｎ（０，ｍａｘ（．．．，１）
2010531752.tif 000301
【０３４６】
2010531752.tif 000302
独立行列要素の各々について２回の乗算を必要とするか、合計で４２回の演算を必要とする（残りの行列要素は、対称性により割り当てられる）。
【０３４７】
2010531752.tif 000303
列要素に対するζ１または３ζ１の多数の割り当てによって計算される。同様に、更新演算
2010531752.tif 000304
総演算数は、１０（すなわち、９減算と、３ζ１の計算のための１演算）に減少する。
【０３４８】
2010531752.tif 000305
よって計算しなければならない（残りは、割り当てにより得られる）。１２個の式のうちの
2010531752.tif 000306
回である。その後、残りの行列要素は、割り当てによって得られる。
【０３４９】
2010531752.tif 000307
めの割り当てのために２１回の演算を必要とする。
【０３５０】
2010531752.tif 000308
要とする、上述のステップを順に実行する（この手順も非常に効率的に設定されているこ
2010531752.tif 000309
１１８＋１８＝１３６である。
【０３５１】  上記に概説した「標準的手法」は、ハイブリッド閉包への適用時に、全体で約３２６回の演算を必要とする。これは、提案された効率的な手法の計算コストを、およそ１．７２倍（すなわち、約７０％）上回る。提案された効率的な手法を使用することによって、「標準的」手法の計算コストの約６０％で、安定化ハイブリッド閉包（動的トレース安定化を含む）を有するＦＴＥの右辺関数の評価の結果を得ることができる。
【０３５２】
2010531752.tif 000310
数の特定の代数構造の賢明な再体系化によって主に達成される。最適な効率レベルは、ＣＮスキームを介して方程式を減少させることによって、最終的に達成される。
【０３５３】  ７．７  トレース再スケーリング、不変量監視、および位相空間射影
【０３５４】  任意のＦＯ行列は、非負固有値および単位トレースを有することを必要とする（→２．４節）。これは、ＦＴＥの解に対する代数的制約の組があることを示唆し、これによって、ＤＡＳになる（→２．６節）。トレースに関する等式Ｓａ＝Ｔｒ（ａ）＝１とは別に、固有値に関する非負制約は、
2010531752.tif 000311
2010531752.tif 000312
化できる。
【０３５５】  トレース条件は、動的トレース安定化（ｄｙｎａｍｉｃａｌ  ｔｒａｃｅ  ｓｔａｂｉｌｉｚａｔｉｏｎ；ＤＴＳ）を介して、「能動的」かつ計算的に低コストな方式で説明される。ここで、ＦＴＥの右辺には、適切なペナルティ項（→７．３節）が追加され、それは数値解のトレースを所要値Ｓａ＝１に近似であるよ
2010531752.tif 000313
レース条件を満たすＦＴＥの数値解を得ることができる。実際的な観点からは、これは、解の正確性にほとんど影響を及ぼさない。というのも、因子１／Ｓａによる乗算により行列要素を再スケールすることが常にできるためである。このトレース再スケーリング演算は、数学的に次のように記述される。
【０３５６】
2010531752.tif 000314
【０３５７】  また、トレース再スケーリング演算は、次のようないくつかの好ましい特性を有する。
2010531752.tif 000315
ように）。
・  対応する固有ベクトルは、再スケーリング演算により影響を受けない。
【０３５８】  ゆえに、トレース再スケーリングの演算は、計算された数値解に対して質的変化をもたらさないが、行列の全ての本質的特性を保ったまま、解をわずかに訂正する。この観点から、非負固定値μｋ（または同等に、非負第２の不変量および第３の不変量ＫａおよびＤａ）
2010531752.tif 000316
【０３５９】  ＦＴＥの不変量制御の一般的側面
【０３６０】  ［１２］の結果（→３節）によると、実対称３×３行列の固有値の非負性は、そのトレースが正であり、第２の不変量および第３の不変量ＫａおよびＤａが非負である場合に保証される。
【０３６１】
2010531752.tif 000317
レース演算を換算するため、トレースに関する発展方程式を直接的に導くことができる（→４節）。ＦＴＥの右辺の特別な代数構造により、この発展方程式は、「厳密な」場合、ならびにハイブリッドまたは一般的な直交異方性閉包近似（式（１４）および（１５ａ／ｂ）参照）のいずれかの場合に、閉形式（ｃｌｏｓｅｄ  ｆｏｒｍ）をとる。ＦＴＥのこのような特別な代数特性により、特別に（意図的に）トレースの発展方程式に影響を及ぼすＦＴＥ自体の右辺にペナルティ項を追加することができる（→ＤＴＳ）。
【０３６２】  しかしながら、トレースとは違って、第２の不変量および第３の不変量ＫａおよびＤａは、行列要素の非線形関数である。ゆえに、これらの不変量の閉形式発展方程式（ｃｌｏｓｅｄ  ｆｏｒｍｅｖｏｌｕｔｉｏｎ  ｅｑｕａｔｉｏｎ）を得るために、同一の（または類似の）手順を適用することはできない（一般には、不変量の発展方程式の右辺は、それらの不変量の組み合わせだけでなく、全行列要素の関数であり、固有ベクトルにおいてコードされる行列の方向性から独立していない）。同様に、任意の特定の項を、既定の方式で非線形不変量の発展に影響を及ぼすＦＴＥの右辺に追加することは、かなり難しい（不可能でないにしても）。ゆえに、不変量を非負領域に維持するために、ＦＴＥの右辺における適切なペナルティ項による第２の不変量および第３の不変量の「能動的」制御は、利用できない。
【０３６３】  代替の手順は、「非制御」（または部分的に制御される）積分手順と、その許容可能領域に数値解をマッピングする射影演算との組み合わせから成る。このような「受動的」手順は、数値積分およびその後の射影の演算が、時間ステップ毎に適用される場合に、「非制御」スキームの精度を保存する有効な積分スキームをもたらす（詳細説明については、射影方法に関する［９］の章ＩＶ．４を参照）。
【０３６４】  ＦＴＥに適用されるこの組み合わせ手順は、以下から成る。
１．第５節、６節、および７節（７．６節まで）に説明する方法を使用する積分ステップ。
２．続いて、前の積分ステップにおいて計算された行列が依然としてＦＯ行列であるかを決定するための不変量監視ステップ。
３．監視ステップの結果が負である場合（すなわち、不変量条件に反する場合）に、位相空間の組ＭＦＴに行列をマッピングする位相空間射影ステップ。
【０３６５】  位相空間ＭＦＴは、凸集合であるため（→第３節における定理１）、実対称３×３行列毎に、「最小距離」（適切な測定基準により測定される）においてＭＦＴに一意的な行列が存在する。これは、上記手順の第３のステップに必要な射影マッピングを定義する。位相空間射影に関するさらなる詳細について後述する。
【０３６６】効率的な不変量監視および固有値計算
【０３６７】  上述の手順によると、不変量の監視は、時間ステップ毎および計算領域の（充填）セル毎に実行しなければならない。ゆえに、不変量（第２の不変量および第３の不変量またはそれ自体が不変量である固有値）の計算は、いかなる不必要な計算オーバーヘッドをも回避して、できるだけ効率的に実行されなければならない。
【０３６８】  第３節で紹介したように、行列の対角要素および非対角要素の記法（χ，ｙ，ｚ）および（ｕ，ν，ｗ）を使用して、不変量は、以下の公式によって、行列要素の多項式として求められる。
【０３６９】
2010531752.tif 000318
【０３７０】  これは、全部で１１回の加算および１３回の乗算で評価できる。トレースＳａは、ＦＴＥの数
2010531752.tif 000319
とによる不変量の監視は、セル毎に時間ステップ当たり追加の２２回の演算により実行できる。
【０３７１】
2010531752.tif 000320
ＫａおよびＤａのチェックは、実対称３×３行列がＦＯ行列の集合ＭＭＴに属するか否かを決定する最も効率的な方式である。
【０３７２】  詳細に分析すると、行列の固有値の数値計算は、その特性多項式Ｐａ（μ）の根を計算することによって、最も効率的に実行され得ることが示される。これは、Ｐａ（μ）の係数として不変量を計算するのに必要な２４回の演算に加えて約１００回の演算のコストがかかるため、約５倍高くなる。これは、ＦＯ行列特性を監視するために、不変量の代わりに固有値を使用
2010531752.tif 000321
性の効率的な監視について最適な方法であることを明確にする。行列の数値的対角化は、固有値を計算する最も高価な方法である。ゆえに、このような手順は、固有ベクトルが利益を有する（例えば、後述の位相空間射影手順の一部として）場合を除いて推奨されない。
【０３７３】位相空間射影
【０３７４】  不変量監視によって、ＦＴＥの数値解が、現在の時間ステップ内で、位相空間集合ＭＦＴの
2010531752.tif 000322
方に反する場合）、解は、位相空間ＭＦＴへの射影により訂正されなければならない。
【０３７５】
2010531752.tif 000323
称行列の６次元部分空間に対するこれらの量の制限は直接的である。［１８］によれば、任意
2010531752.tif 000324
【０３７６】
2010531752.tif 000325
れる（正式な証明に関しては［１８］を参照）。
2010531752.tif 000326
節における式（７）および図７参照）
2010531752.tif 000327
ユークリッド距離を有する。
【０３７７】
2010531752.tif 000328
において求められる。
【０３７８】
2010531752.tif 000329
【０３７９】  また、射影マッピングの値を実際に計算するアルゴリズムは、以下のステップから成る。
【０３８０】
2010531752.tif 000330
を計算する。
【０３８１】
2010531752.tif 000331
この場合、入力行列を出力に割り当て、終了する。そうでない場合（すなわち、Ｋａ＜０またはＤａ＜０である場合）、以下のステップを実行する。
【０３８２】
2010531752.tif 000332
算する。このステップは、数式として以下のように表現され得る。
2010531752.tif 000333
【０３８３】
2010531752.tif 000334
【０３８４】
2010531752.tif 000335
【０３８５】
2010531752.tif 000336
と一致することに留意されたい。この意味において、（３９ｂ／ｃ／ｄ）は、既に、射影マッピングの一般的定義をコードしている。（３９ｂ／ｃ／ｄ）の結果を計算する個々のステップは、
2010531752.tif 000337
合に、定義を数式的に完成させる。
【０３８６】
2010531752.tif 000338
列を取るこのアルゴリズムの実際的な実装は、手順のステップ１においてトレースの追加の評価を省略し得る。さらに、実際的な実装は、種々のアルゴリズムステップの以下の変形および／または仕様を組み込み得る。図１３を参照されたい。
【０３８７】  ・  ステップ１では、入力行列の不変量ＫａおよびＤａを（３８）によって計算する。
【０３８８】
2010531752.tif 000339
る場合には、入力行列をそのまま出力に割り当てる。
【０３８９】  ・  数値的行列対角化によって、例えば、反復周期ヤコビ回転（［３８］、章１１．１または［３９］、章６．５．２）またはギブンズ／ハウスホルダー低減ステップおよび後続のＱＲ繰り返し（［４２］または［３８］、章１１．２を参照）によって、ステップ３の結果を得る。
【０３９０】  ・  実際には、ステップ４は、入力固有値μｋのうちの少なくとも１つが負である場合に
2010531752.tif 000340
含む）に位置する。
2010531752.tif 000341
することは、（３９）に従ってμ＊を決定するために解かなければならない最小化問題を、３角形Δａの３つの端点（１，０，０），（０，１，０），（０，０，１）によって画定される面Ｓａ＝１に制限する。この「面」最小化問題は、以下の射影された３角形
2010531752.tif 000342
例えば、（μ１，μ２）面において同等に解決できる。次いで、μ＊の第３座標は、
2010531752.tif 000343
算するアルゴリズムよりも単純であり、より効率的に実装できる。
【０３９１】
2010531752.tif 000344
【０３９２】  ８．実装
【０３９３】  製品の所定の部分における繊維の向きに関する可能性を記述する分布関数であるプロセスの結果は、コンピュータワークステーション（図示せず）のディスプレイにおいてグラフ形式または数値形式で提示される。これは、その上で計算が実行されるコンピュータまたはワークステーションのディスプレイであってもよく、または、これは、その上でシミュレーションが実行されるコンピュータにネットワークで接続されるコンピュータのディスプレイであってもよい。
【０３９４】  モールドまたは製品の設計者は、射出成形プロセスにおける製品の質を改善するために、シミュレーションの結果を使用するだろう。また、シミュレーション結果は、関連する製品の製造コストを低減させるために技術者によっても使用されうる。典型的にはポリマー材料よりもさらに高い強度特徴を有する繊維の向きが、製品の強度、剛性、および安定性に対して決定的な影響を有することから、繊維情報に関する信頼できる情報の利点によって、技術者は、製品の強度、剛性、および安定性の特性を決定する前に、良く理解し、情報を有することが可能になる。繊維の向きは、中に懸濁される繊維を含む塊状ポリマーを射出成形する際に発生し得る整経効果に影響を及ぼす。繊維の向きを把握することによって、影響する整経や他の応力を十分に予測することが可能になり、また設計変更を回避することができる。
【０３９５】  また、シミュレーションの結果は、ＣＡＤソフトウェア等の他のアプリケーションに渡すこともできる。したがって、シミュレーションの結果は、設計ソフトウェアとシミュレーションソフトウェアとの間の対話型プロセスにおいて使用できる。
【０３９６】  本発明は、多数の利点を有する。異なる実施形態または実装によって、以下の利点のうちの１つ以上がもたらされ得る。これが包括的なリストではなく、本明細書に説明されない他の利点も存在し得ることに留意されたい。本発明の１つの利点は、繊維強化製品における繊維配向の分布が、大幅に低減された計算量で決定できることである。本発明の利点は、繊維強化製品における繊維配向の分布が、高い精度で決定できることである。本発明の別の利点は、繊維強化製品における繊維配向の分布がより高速に決定できることである。
【０３９７】  本発明について例証目的で詳細に説明したが、このような詳細は、単にその目的に関するだけであり、本発明の範囲を逸脱することなく当業者によって変形を加えてもよいことを理解されたい。
【０３９８】  請求項において使用される用語の「備える」は、他の要素またはステップを除外しない。請求項において使用される単数形の用語は、複数形を除外しない。
【０３９９】参考文献
［１］ＪＳ．Ｇ．Ａｄｖａｎｉ（Ｅｄ．）：Ｆｌｏｗ  ａｎｄ  Ｒｈｅｏｌｏｇｙ  ｉｎ  Ｐｏｌｙｍｅｒ  Ｃｏｍｐｏｓｉｔｅｓ  Ｍａｎｕｆａｃｔｕｒｉｎｇ，Ｅｌｓｅｖｉｅｒ，Ａｍｓｔｅｒｄａｍ  （１９９４）
［２］Ｇ．Ｂ．Ｊｅｆｆｅｒｙ：  Ｔｈｅ  ｍｏｔｉｏｎ  ｏｆ  ｅｌｌｉｐｓｏｉｄａｌ  ｐａｒｔｉｃｌｅｓ  ｉｍｍｅｒｓｅｄ  ｉｎ  ａ  ｖｉｓｃｏｕｓ  ｆｌｕｉｄ，Ｐｒｏｃ．Ｒ．Ｓｏｃ．Ａ  １０２，１６１−１７９  （１９２２）
［３］Ｍ．Ｊｕｎｋ  ａｎｄ  Ｒ．Ｉｌｌｎｅｒ：Ａ  ｎｅｗ  ｄｅｒｉｖａｔｉｏｎ  ｏｆ  Ｊｅｆｆｅｒｙ’ｓ  ｅｑｕａｔｉｏｎ，Ｊ．Ｍａｔｈ  Ｆｌｕｉｄ  Ｍｅｃｈ．８，１−３４  （２００６）
［４］Ｆ．Ｆｏｌｇａｒ  ａｎｄ  ＣＬ．Ｔｕｃｋｅｒ  ＩＩＩ：Ｏｒｉｅｎｔａｔｉｏｎ  ｂｅｈａｖｉｏｕｒ  ｏｆ  ｆｉｂｅｒｓ  ｉｎ  ｃｏｎｃｅｎｔｒａｔｅｄ  ｓｕｓｐｅｎｓｉｏｎｓ，Ｊ．Ｒｅｉｎｆ．Ｐｌａｓｔ．Ｃｏｍｐｏｓ．３，９８−１１９  （１９８４）
［５］ＣＬ．Ｔｕｃｋｅｒ  ａｎｄ  Ｓ．Ｇ．Ａｄｖａｎｉ：Ｐｒｏｃｅｓｓｉｎｇ  ｏｆ  ｓｈｏｒｔ−ｆｉｂｅｒ  ｓｙｓｔｅｍｓ，ｃｈ．６，ｐｐ．１４７  ｉｎ  Ｓ．Ｇ．Ａｄｖａｎｉ  （Ｅｄ．）：Ｆｌｏｗ  ａｎｄ  Ｒｈｅｏｌｏｇｙ  ｉｎ  Ｐｏｌｙｍｅｒ  Ｃｏｍｐｏｓｉｔｅｓ  Ｍａｎｕｆａｃｔｕｒｉｎｇ（ｓｅｅ  ｒｅｆ．［１］ａｂｏｖｅ）
［６］Ｓ．Ｇ．Ａｄｖａｎｉ  ａｎｄ  ＣＬ．Ｔｕｃｋｅｒ  ＩＩＩ：Ｔｈｅ  ｕｓｅ  ｏｆ  ｔｅｎｓｏｒｓ  ｔｏ  ｄｅｓｃｒｉｂｅ  ａｎｄ  ｐｒｅｄｉｃｔ  ｆｉｂｅｒ  ｏｒｉｅｎｔａｔｉｏｎ  ｉｎ  ｓｈｏｒｔ  ｆｉｂｅｒ  ｃｏｍｐｏｓｉｔｅｓ，Ｊ．Ｒｈｅｏｌ．７５１−７８４  （１９８７）
［７］Ｓ．Ｇ．Ａｄｖａｎｉ  ａｎｄ  ＣＬ．Ｔｕｃｋｅｒ  ＩＩＩ：Ｃｌｏｓｕｒｅ  ａｐｐｒｏｘｉｍａｔｉｏｎｓ  ｆｏｒ  ｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌ  ｓｔｒｕｃｔｕｒｅ  ｔｅｎｓｏｒｓ，Ｊ，Ｒｈｅｏｌ．，３６７−３８６  （１９９０）
［８］Ｊ．Ｓ．Ｃｉｎｔｒａ  ａｎｄ  ＣＬ．Ｔｕｃｋｅｒ  ＩＩＩ：Ｏｒｔｈｏｔｒｏｐｉｃ  ｃｌｏｓｕｒｅ  ａｐｐｒｏｘｉｍａｔｉｏｎｓ  ｆｏｒ  ｆｌｏｗ−ｉｎｄｕｃｅｄ  ｆｉｂｅｒ  ｏｒｉｅｎｔａｔｉｏｎ，Ｊ，Ｒｈｅｏｌ．３９，１０９５−１１２２  （１９９５）
［９］Ｅ．Ｈａｉｒｅｒ、Ｃ．Ｌｕｂｉｃｈ、  ａｎｄ  Ｇ．Ｗａｎｎｅｒ：Ｇｅｏｍｅｔｒｉｃ  ｎｕｍｅｒｉｃａｌ  ｉｎｔｅｇｒａｔｉｏｎ，Ｓｐｒｉｎｇｅｒ，Ｂｅｒｌｉｎ  （２００２）
［１０］Ｊ．Ｌｉｎｎ，Ｊ．Ｓｔｅｉｎｂａｃｈ，Ａ．Ｒｅｉｎｈａｒｄｔ：Ｃａｌｃｕｌａｔｉｏｎ  ｏｆ  ｔｈｅ  ３Ｄ  ｆｉｂｅｒ  ｏｒｉｅｎｔａｔｉｏｎ  ｉｎ  ｔｈｅ  ｓｉｍｕｌａｔｉｏｎ  ｏｆ  ｔｈｅ  ｉｎｊｅｃｔｉｏｎ  ｍｏｌｄｉｎｇ  ｐｒｏｃｅｓｓ  ｆｏｒ  ｓｈｏｒｔ−ｆｉｂｅｒ  ｒｅｉｎｆｏｒｃｅｄ  ｔｈｅｒｍｏｐｌａｓｔｓ，ＥＣＭＩ  ２０００  Ｃｏｎｆｅｒｅｎｃｅ，Ｐａｌｅｒｍｏ  （２０００）
［１１］Ｊ．Ｌｉｎｎ：Ｅｘｐｌｏｒｉｎｇ  ｔｈｅ  ｐｈａｓｅ  ｓｐａｃｅ  ｏｆ  ｔｈｅ  Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ  ｅｑｕａｔｉｏｎ，ＳＩＡＭ−ＥＭＳ  Ｃｏｎｆｅｒｅｎｃｅ，Ｂｅｒｌｉｎ  （２００１）
［１２］Ｊ．Ｌｉｎｎ：Ｏｎ  ｔｈｅ  ｆｒａｍｅ−ｉｎｖａｒｉａｎｔ  ｄｅｓｃｒｉｐｔｉｏｎ  ｏｆ  ｔｈｅ  ｐｈａｓｅ  ｓｐａｃｅ  ｏｆ  ｔｈｅ  Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ  ｅｑｕａｔｉｏｎ，ｐ．３２７−３３２  ｉｎ  Ａ．Ｂｕｉｋｉｓ，Ｒ．Ｃｉｅｇｉｓ，Ａ．Ｄ．Ｆｉｔｔ  （Ｅｄｓ．）：Ｐｒｏｇｒｅｓｓ  ｉｎ  Ｉｎｄｕｓｔｒｉａｌ  Ｍａｔｈｅｍａｔｉｃｓ  ａｔ  ＥＣＭＩ  ２００２，Ｓｐｒｉｎｇｅｒ  （２００４）
［１３］Ｓ．Ｈｅｓｓ：Ｆｏｋｋｅｒ−Ｐｌａｎｃｋ  ｅｑｕａｔｉｏｎ  ａｐｐｒｏａｃｈ  ｔｏ  ｆｌｏｗ  ａｌｉｇｎｍｅｎｔ  ｉｎ  ｌｉｑｕｉｄ  ｃｒｙｓ−ｔａｌｓ，Ｚ．Ｎａｔｕｒｆｏｒｓｃｈ．３１Ａ，１０３４  ｆｆ．  （１９７６）
［１４］Ｍ．Ｄｏｉ：Ｍｏｌｅｃｕｌａｒ  ｄｙｎａｍｉｃｓ  ａｎｄ  ｒｈｅｏｌｏｇｉｃａｌ  ｐｒｏｐｅｒｔｉｅｓ  ｏｆ  ｃｏｎｃｅｎｔｒａｔｅｄ  ｓｏｌｕｔｉｏｎｓ  ｏｆ  ｒｏｄｌｉｋｅ  ｐｏｌｙｍｅｒｓ  ｉｎ  ｉｓｏｔｒｏｐｉｃ  ａｎｄ  ｌｉｑｕｉｄ  ｃｒｉｓｔａｌｌｉｎｅ  ｐｈａｓｅｓ，Ｊ．Ｐｏｌｙｍ．ＳｃＬ，Ｐｏｌｙｍ．Ｐｈｙｓ．Ｅｄ．１９，２２９−２４３  （１９８１）
［１５］Ｍ．Ｇｒｏｓｓｏ，Ｐ．Ｌ．Ｍａｆｆｅｔｏｎｅ，Ｆ．Ｄｕｐｒｅｔ：Ａ  ｃｌｏｓｕｒｅ  ａｐｐｒｏｘｉｍａｔｉｏｎ  ｆｏｒ  ｎｅｍａｔｉｃ  ｌｉｑｕｉｄ  ｃｒｙｓｔａｌｓ  ｂａｓｅｄ  ｏｎ  ｔｈｅ  ｃａｎｏｎｉｃａｌ  ｄｉｓｔｒｉｂｕｔｉｏｎ  ｓｕｂｓｐａｃｅ  ｔｈｅｏｒｙ，Ｒｈｅｏｌ，Ａｃｔａ３９，３０１−３１０  （２０００）
［１６］Ｍ．Ｋｒｏｇｅｒ：Ｓｉｍｐｌｅ  ｍｏｄｅｌｓ  ｆｏｒ  ｃｏｍｐｌｅｘ  ｎｏｎｅｑｕｉｌｉｂｒｉｕｍ  ｆｌｕｉｄｓ，Ｐｈｙｓ，Ｒｅｐ．３９０，４５３−５５１  （２００４）
［１７］Ｊ．Ｌｉｎｎ：Ｔｈｅ  Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ  Ｍｏｄｅｌ  ａｓ  ａ  Ｄｉｆｆｅｒｅｎｔｉａｌ  Ａｌｇｅｂｒａｉｃ  Ｓｙｓｔｅｍ  ｆｏｒ  Ｆｉｂｅｒ  Ｏｒｉｅｎｔａｔｉｏｎ  Ｃａｌｃｕｌａｔｉｏｎ，ｐｐ．２１５−２２４  ｉｎ  Ｙ．Ｗａｎｇ，Ｋ．Ｈｕｔｔｅｒ：Ｔｒｅｎｄｓ  ｉｎ  Ａｐ−ｐｌｉｃａｔｉｏｎｓ  ｏｆ  Ｍａｔｈｅｍａｔｉｃｓ  ｔｏ  Ｍｅｃｈａｎｉｃｓ，ｐｒｏｃｅｅｄｉｎｇｓ  ｏｆ  ｔｈｅ  ＳＴＡＭＭ  ２００４  ｃｏｎｆｅｒｅｎｃｅ  ｉｎ  Ｓｅｅｈｅｉｍ  （Ｇｅｒｍａｎｙ），Ｓｈａｋｅｒ  （２００５）
［１８］Ｕ．Ｓｔｒａｕｔｉｎｓ：Ｉｎｖｅｓｔｉｇａｔｉｏｎ  ｏｆ  ｆｉｂｅｒ  ｏｒｉｅｎｔａｔｉｏｎ  ｄｙｎａｍｉｃｓ  ｗｉｔｈｉｎ  ｔｈｅ  Ｆｏｌｇａｒ−Ｔｕｃｋｅｒ  ｍｏｄｅｌ  ｗｉｔｈ  ｈｙｂｒｉｄ  ｃｌｏｓｕｒｅ，ｍａｓｔｅｒ  ｔｈｅｓｉｓ，Ｄｅｐｔ．ｏｆ  Ｍａｔｈｅｍａｔｉｃｓ，Ｕｎｉｖｅｒｉｔｙ  ｏｆ  Ｋａｉｓｅｒｓｌａｕｔｅｒｎ  （２００４）
［１９］Ｊ．Ｌｉｎｎ：Ｄｒｅｉｄｉｍｅｎｓｉｏｎａｌｅ  Ｖυｒａｉｉｓｂｅｒｅｃｈｎｕｎｇ  ｄｅｒ  ａｎｉｓｏｔｒｏｐｅｎ  Ｍａｔｅｒｉａｌｅｉｇｅｎ−ｓｃｈａｆｔｅｎ  ｉｎ  ｔｈｅｒｍｏｐｌａｓｔｉｓｃｈｅｎ  Ｓｐｒｉｔｚｇｕｓｓｅｒｚｅｕｇｎｉｓｓｅｎ  （Ａｎｌｓｏ−３Ｄ），Ｐｒｏｊｅｋｔ−ｂｅｒｉｃｈｔ  ｆｕｒ  ｄｉｅ  ＭＡＧＭＡ  ＧｍｂＨ，Ｔｅｉｌ  Ｈａ  （２００２）
［２０Ｊ  Ｊ．Ｌｉｎｎ：Ｄｒｅｉｄｉｍｅｎｓｉｏｎａｌｅ  Ｖｏｒａｕｓｂｅｒｅｃｈｎｕｎｇ  ｄｅｒ  ａｎｉｓｏｔｒｏｐｅｎ  Ｍａｔｅｒｉａｌｅｉｇｅｎ−ｓｃｈａｆｔｅｎ  ｉｎ  ｔｈｅｒｍｏｐｌａｓｔｉｓｃｈｅｎ  Ｓｐｒｉｔｚｇｕｓｓｅｒｚｅｕｇｎｉｓｓｅｎ  （Ａｎｌｓｏ−３Ｄ），Ｐｒｏｊｅｋｔ−ｂｅｒｉｃｈｔ  ｆｕｒ  ｄｉｅ  ＭＡＧＭＡ  ＧｍｂＨ，Ｔｅｉｌ  ｌｉｂ  （２００２）
［２１］Ｂ．Ｅ．ＶｅｒＷｅｙｓｔ，ＣＬ．Ｔｕｃｋｅｒ，Ｐ．Ｈ．Ｆｏｓｓ，Ｊ．Ｆ．Ｏ’Ｇａｒａ：Ｆｉｂｅｒ  ｏｒｉｅｎｔａｔｉｏｎ  ｉｎ  ３Ｄ  ｉｎｊｅｃｔｉｏｎ  ｍｏｕｌｄｅｄ  ｆｅａｔｕｒｅｓ：ｐｒｅｄｉｃｔｉｏｎ  ａｎｄ  ｅｘｐｅｒｉｍｅｎｔ，Ｉｎｔｅｒｎａｔ．Ｐｏｌｙｍｅｒ  Ｐｒｏｃｅｓｓｉｎｇ  １４，４０９−４２０  （１９９９）；
［２２］Ｂ．Ｅ．ＶｅｒＷｅｙｓｔ：Ｎｕｍｅｒｉｃａｌ  ｐｒｅｄｉｃｔｉｏｎｓ  ｏｆ  ｆｌｏｗ−ｉｎｄｕｃｅｄ  ｆｉｂｅｒ  ｏｒｉｅｎｔａｔｉｏｎ  ｉｎ  ｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌ  ｇｅｏｍｅｔｒｉｅｓ，Ｐｈ．Ｄ  ｔｈｅｓｉｓ，Ｕｎｉｖ．ｏｆ  Ｉｌｌｉｎｏｉｓ  ａｔ  Ｕｒｂａｎａ  Ｃｈａｍ−ｐａｉｇｎ  （１９９８）
［２３］Ｇ．Ｉ．Ｍａｒｃｈｕｋ：Ｓｐｌｉｔｔｉｎｇ  ａｎｄ  Ａｌｔｅｒｎａｔｉｎｇ  Ｄｉｒｅｃｔｉｏｎ  Ｍｅｔｈｏｄｓ，ｐｐ．１９７−４６２  ｉｎ  Ｐ．Ｇ．Ｃｉａｒｅｔ  ＆  Ｊ．Ｌ  Ｌｉｏｎｓ  （Ｅｄｓ．）：Ｈａｎｄｂｏｏｋ  ｏｆ  Ｎｕｍｅｒｉｃａｌ  Ａｎａｌｙｓｉｓ，Ｖｏｌｕｍｅ  Ｉ，Ｎｏｒｔｈ−Ｈｏｌｌａｎｄ，Ｅｌｓｅｖｉｅｒ  （１９９０）
［２４］Ｋ．Ｗ．Ｍｏｒｔｏｎ：Ｎｕｍｅｒｉｃａｌ  ｓｏｌｕｔｉｏｎ  ｏｆ  ｃｏｎｖｅｃｔｉｏｎ−ｄｉｆｆｕｓｉｏｎ  ｐｒｏｂｌｅｍｓ，Ｃｈａｐｍａｎ＆Ｈａｌｌ，Ｌｏｎｄｏｎ  （１９９６）
［２５］Ｒ．Ｊ．ＬｅＶｅｑｕｅ，Ｎｕｍｅｒｉｃａｌ  Ｍｅｔｈｏｄｓ  ｆｏｒ  Ｃｏｎｓｅｒｖａｔｉｏｎ  Ｌａｗｓ，Ｂｉｒｋｈａｕｓｅｒ  （１９９２）
［２６］Ｇ．Ｓｔｒａｎｇ：”Ｏｎ  ｔｈｅ  ｃｏｎｓｔｒｕｃｔｉｏｎ  ａｎｄ  ｃｏｍｐａｒｉｓｏｎ  ｏｆ  ｄｉｆｆｅｒｅｎｃｅ  ｓｃｈｅｍｅｓ”，ＳＩＡＭ  Ｊｏｕｒｎ．Ｎｕｍ．Ａｎａｌ．５，５０６−５１７  （１９６８）
［２７］Ｍ．Ｇ．Ｃｒａｎｄａｌｌ  ａｎｄ  Ａ．Ｍａｊｄａ：”Ｔｈｅ  ｍｅｔｈｏｄ  ｏｆ  ｆｒａｃｔｉｏｎａｌ  ｓｔｅｐｓ  ｆｏｒ  ｃｏｎｓｅｒｖａｔｉｏｎ  ｌａｗｓ”，Ｍａｔｈ．Ｃｏｍｐ．３４，２８５−３１４  （１９８０）
［２８］Ｈ．Ｖ．Ｋｏｊｏｕｈａｒｏｖ，Ｂ．Ｍ．Ｃｈｅｎ：”Ｎｏｎｓｔａｎｄａｒｄ  ｍｅｔｈｏｄｓ  ｆｏｒ  ｔｈｅ  ｃｏｎｖｅｃｔｉｖｅ  ｔｒａｎｓｐｏｒｔ  ｅｑｕａｔｉｏｎ  ｗｉｔｈ  ｎｏｎｌｉｎｅａｒ  ｒｅａｃｔｉｏｎｓ”，Ｎｕｍｅｒ．Ｍｅｔｈ  Ｐａｒｔｉａｌ  Ｄｉｊｆ．Ｅｑ．１４，４６７−４８５  （１９９８）；”Ｎｏｎｓｔａｎｄａｒｄ  ｍｅｔｈｏｄｓ  ｆｏｒ  ｔｈｅ  ｃｏｎｖｅｃｔｉｖｅ−ｄｉｓｐｅｒｓｉｖｅ  ｔｒａｎｓｐｏｒｔ  ｅｑｕａｔｉｏｎ  ｗｉｔｈ  ｎｏｎｌｉｎｅａｒ  ｒｅａｃｔｉｏｎｓ”ｉｎ  Ｒ．Ｅ．Ｍｉｃｋｅｎｓ  （ｅｄ．）：Ａｐｐｌｉｃａ−ｔｉｏｎｓ  ｏｆ  ｎｏｎ−ｓｔａｎｄａｒｄ  ｆｉｎｉｔｅ  ｄｉｆｆｅｒｅｎｃｅ  ｓｃｈｅｍｅｓ，ｍｉｎｉｓｙｍｐｏｓｉｕｍ  ｏｎ  ｎｏｎ−ｓｔａｎｄａｒｄ  ｆｉｎｉｔｅ  ｄｉｆｆｅｒｅｎｃｅ  ｓｃｈｅｍｅｓ：ｔｈｅｏｒｙ  ａｎｄ  ａｐｐｌｉｃａｔｉｏｎｓ，ＳＩＡＭ  ａｎｎｕａｌ  ｍｅｅｔｉｎｇ，Ａｔｌａｎｔａ  ＧＡ，ＵＳＡ  １９９９，ｐｕｂｌ．ｂｙ  Ｓｉｎｇａｐｏｒｅ：Ｗｏｒｌｄ  Ｓｃｉｅｎｔｉｆｉｃ  （２０００）
［２９］Ｈ．Ｗａｎｇ、Ｘ．Ｓｈｉ、  ａｎｄ  Ｒ．Ｅ．Ｅｗｉｎｇ：”Ａｎ  ＥＬＬＥＭ  ｓｃｈｅｍｅ  ｆｏｒ  ｍｕｌｔｉｄｉｍｅｎｓｉｏｎａｌ  ａｄｖｅｃｔｉｏｎ−ｒｅａｃｔｉｏｎ  ｅｑｕａｔｉｏｎｓ  ａｎｄ  ｉｔｓ  ｏｐｔｉｍａｌ−ｏｒｄｅｒ  ｅｒｒｏｒ  ｅｓｔｉｍａｔｅ”，ＳＩＡＭ  Ｊ．  Ｎｕｍｅｒ．Ａｎａｌ．３８，１８４６−１８８５  （２００１）
［３０］ＰＪ．ｖａｎ  ｄｅｒ  Ｈｏｕｗｅｎ：”Ｎｏｔｅ  ｏｎ  ｔｈｅ  ｔｉｍｅ  ｉｎｔｅｇｒａｔｉｏｎ  ｏｆ  ３Ｄ  ａｄｖｅｃｔｉｏｎ−ｒｅａｃｔｉｏｎ  ｅｑｕａｔｉｏｎｓ”，Ｊ．Ｃｏｍｐｕｔ．Ａｐｐｌ．Ｍａｔｈ．１１６，２７５−２７８  （２０００）
［３１］Ｗ．  Ｈｕｎｓｄｏｒｆｅｒ，Ｊ．Ｇ．Ｖｅｒｗｅｒ：”Ａ  ｎｏｔｅ  ｏｎ  ｓｐｌｉｔｔｉｎｇ  ｅｒｒｏｒｓ  ｆｏｒ  ａｄｖｅｃｔｉｏｎ−ｒｅａｃｔｉｏｎ  ｅｑｕａｔｉｏｎｓ”，Ａｐｐｌ．Ｎｕｍｅｒ．Ｍａｔｈ．１８，１９１−１９９  （１９９５）
［３２］Ｓ．Ｖ．Ｐａｔａｎｋａｒ：Ｎｕｍｅｒｉｃａｌ  ｈｅａｔ  ｔｒａｎｓｆｅｒ  ａｎｄ  ｆｌｕｉｄ  ｆｌｏｗ，Ｈｅｍｉｓｐｈｅｒｅ  Ｐｕｂｌ．Ｃｏｒｐ．  （１９８０）
［３３］ＣＡＪ．Ｆｌｅｔｃｈｅｒ：Ｃｏｍｐｕｔａｔｉｏｎａｌ  ｔｅｃｈｎｉｑｕｅｓ  ｆｏｒ  Ｆｌｕｉｄ  Ｄｙｎａｍｉｃｓ，Ｖｏｌｕｍｅ  Ｉ：Ｆｕｎｄａｍｅｎｔａｌ  ａｎｄ  Ｇｅｎｅｒａｌ  Ｔｅｃｈｎｉｑｕｅｓ  （２ｎｄ  ｅｄｉｔｉｏｎ），Ｓｐｒｉｎｇｅｒ  （１９９１）
［３４］Ｌ．Ｆ．Ｓｈａｍｐｉｎｅ：”Ｃｏｎｓｅｒｖａｔｉｏｎ  ｌａｗｓ  ａｎｄ  ｔｈｅ  ｎｕｍｅｒｉｃａｌ  ｓｏｌｕｔｉｏｎ  ｏｆ  ＯＤＥｓ”，Ｃｏｍｐ．Ｍａｔｈ．Ａｐｐｌｉｃ．１２Ｂ，１２８７−１２９６  （１９８６）
［３５］Ｌ．Ｆ．Ｓｈａｍｐｉｎｅ：”Ｌｉｎｅａｒ  ｃｏｎｓｅｒｖａｔｉｏｎ  ｌａｗｓ  ｆｏｒ  ＯＤＥｓ”，Ｃｏｍｐ．Ｍａｔｈ．Ａｐｐｌｉｃ．３５，４５−５３  （１９９８）
［３６］Ｌ．Ｆ．Ｓｈａｍｐｉｎｅ：”Ｃｏｎｓｅｒｖａｔｉｏｎ  ｌａｗｓ  ａｎｄ  ｔｈｅ  ｎｕｍｅｒｉｃａｌ  ｓｏｌｕｔｉｏｎ  ｏｆ  ＯＤＥｓ，ｐａｒｔ  ＩＩ”，Ｃｏｍｐ．Ｍａｔｈ．Ａｐｐｌｉｃ．３８，６１−７２  （１９９９）
［３７Ｊ  Ｊ．Ｌｉｎｎ：”Ｅｎｔｗｉｃｋｌｕｎｇ  ｅｉｎｅｓ  Ｓｏｆｔｗａｒｅ−Ｍｏｄｕｌｓ  ｚｕｒ  Ｂｅｒｅｃｈｎｕｎｇ  ｄｅｒ  Ｆａｓｅｒｏｒｉｅｎ−ｔｉｅｒｕｎｇ  ｉｎ  ｄｅｒ  ＳｐｒｉｔｚｇｉｅＢｓｉｍｕｌａｔｉｏｎ  ｍｉｔ  ＳＩＧＭＡＳＯＦＴ”，ｔｅｃｈｎｉｃａｌ  Ｒｅｐｏｒｔ  ｆｏｒ  ｔｈｅ  ＭＡＧＭＡ  ＧｍｂＨ  （２００１）
［３８］Ｗ．Ｈ．Ｐｒｅｓｓ，Ｓ．Ａ．Ｔｅｕｋｏｌｓｋｙ，Ｗ．Ｔ．Ｖｅｔｔｅｒｌｉｎｇ，Ｂ．Ｐ．Ｆｌａｎｎｅｒｙ：Ｎｕｍｅｒｉｃａｌ  Ｒｅｃｉｐｐｅｅｓｓ  ｉｉｎｎ  ＦＦｏｏｒｒｔｔｒｒａａｎｎ  ７７７７：：ＴＴｈｈｅｅ  Ａ  Ａｒｔ  ｏｆ  Ｓｃｉｅｎｔｉｆｉｃ  Ｃｏｍｐｕｔｉｎｇ  （２”Ｅｄｉｔｉｏｎ），Ｃａｍｂｒｉｄｇｅ  Ｕｎｉｖｅｒｓｉｔｙ  Ｐｒｅｓｓ  （１９９２）
［３９］Ｊ．Ｓｔｏｅｒ，Ｒ．Ｂｕｌｉｒｓｃｈ：Ｉｎｔｒｏｄｕｃｔｉｏｎ  ｔｏ  Ｎｕｍｅｒｉｃａｌ  Ａｎａｌｙｓｉｓ（３ｒｄ  Ｅｄｉｔｉｏｎ），Ｓｐｒｉｎｇｅｒ  （２００２）
［４０］Ｄ．Ｈ．Ｃｈｕｎｇ，Ｔ．Ｈ．Ｋｗｏｎ：”Ａｎ  ｉｎｖａｒｉａｎｔ−ｂａｓｅｄ  ｏｐｔｉｍａｌ  ｆｉｔｔｉｎｇ  ｃｌｏｓｕｒｅ  ａｐｐｒｏｘｉｍａｔｉｏｎ  ｆｏｒ  ｔｈｅ  ｎｕｍｅｒｉｃａｌ  ｐｒｅｄｉｃｔｉｏｎ  ｏｆ  ｆｌｏｗ−ｉｎｄｕｃｅｄ  ｆｉｂｒｅ  ｏｒｉｅｎｔａｔｉｏｎ”，Ｊ．Ｒｈｅｏｌ，４６，１６９−１９４  （２００２）
［４１］Ｆ．Ｄｕｐｒｅｔ，Ｖ．Ｖｅｒｌｅｙｅ，Ｂ．Ｌａｎｇｕｉｌｉｅｒ：”Ｎｕｍｅｒｉｃａｌ  ｐｒｅｄｉｃｔｉｏｎ  ｏｆ  ｍｏｕｌｄｉｎｇ  ｏｆ  ｓｈｏｒｔ−ｆｉｂｒｅ  ｃｏｍｐｏｓｉｔｅ  ｐａｒｔｓ”，Ｐｒｏｃ．１ｓｔ  ＥＳＡＦＯＲＭ  Ｃｏｎｆ．，２９１−２９４  （１９９８）
［４２］Ｇ．Ｈ．Ｇｏｌｕｂ，Ｈ．Ａ．ｖａｎ  ｄｅｒ  Ｖｏｒｓｔ：”Ｅｉｇｅｎｖａｌｕｅ  Ｃｏｍｐｕｔａｔｉｏｎ  ｉｎ  ｔｈｅ  ２０ｔｈ  Ｃｅｎ  ｔｕｒｙ”，Ｊ．Ｃｏｍｐ．Ａｐｐｌ．Ｍａｔｈ．１２３，３５−６５  （２０００）
【図１】
2010531752.tif 000351
【図２】
2010531752.tif 000352
【図３】
2010531752.tif 000353
【図４】
2010531752.tif 000354
【図５】
2010531752.tif 000355
【図６】
2010531752.tif 000356
【図７】
2010531752.tif 000357
【図８】
2010531752.tif 000358
【図９】
2010531752.tif 000359
【図１０】
2010531752.tif 000360
【図１１】
2010531752.tif 000361
【図１２】
2010531752.tif 000362
【図１３】
2010531752.tif 000363
【手続補正書】
【提出日】平成22年1月4日(2010.1.4)
【手続補正１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
シミュレーション領域の幾何学的形状の少なくとも一部分を形成するモールドキャビティが、多数の非球形粒子を含む溶媒により形成される懸濁液で充填又は部分的に充填される射出成形プロセスを、シミュレーションするためのシミュレーションモデルを使用することによって、巨視的レベルにおいて非球形粒子の配向統計を計算する方法であって、前記シミュレーション領域の幾何学的形状のデジタル表現またはコンピュータモデルが提供され、前記シミュレーション領域の少なくとも一部分を再分割または離散化することによって複数の計算セルによるメッシュが形成され、
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テンソルまたは行列の提案された解が、前記テンソルまたは行列が非負定値であり、かつ単位トレースを有するという点において許容される解であるどうかかを決定するための監視ステップを含み、許容される解でない場合は、前記提案された解を許容される位相空間に射影される位相空間射影ステップを実行するか、トレース再スケーリングを行うステップを実行する、方法。
【請求項１５】
請求項１から１４のいずれか１項に記載の方法をコンピュータに実行させる、コンピュータ・プログラム。
【請求項１６】
鋳造装置と、請求項１５に記載のコンピュータ・プログラムにより制御される、前記鋳造装置に接続されたコンピュータとを備える、鋳造システム。
【国際調査報告】
2010531752.tif 000364
2010531752.tif 000365
2010531752.tif 000366
2010531752.tif 000367