JP2019101939A - 流体中の繊維状物質の運動状態の解析方法及びその解析装置 - Google Patents

Abstract

【課題】流体中の繊維状物質の運動状態をより短時間でかつより精度高くシミュレーションすることを可能とする流体中の繊維状物質の運動状態の解析方法を提供する。【解決手段】ステップａ：繊維状物質が複数の球体が結合した球体集合体からなるものと仮定する球体集合体モデルを用いる計算方法を用いて、特定の時刻から微小時間経過後における前記繊維状物質の運動を解析する。ステップｂ：繊維状物質の各球体の流体抵抗力の情報に基づいて、流体が繊維状物質から受ける力を計算し、流体が繊維状物質から受ける力と流体中の繊維状物質の体積分率とを考慮して、特定の時刻から微小時間経過後における流体の速度を計算する。ステップｃ：過渡的な流体及び繊維状物質の運動状態の情報を利用して、新たにステップａ〜ｂの処理を、所定条件を満たすまで繰り返し実行して、流体の運動と繊維状物質の運動を連成解析し、流体中の繊維状物質の運動状態を解析する。【選択図】図１

Description

本発明は、流体中の繊維状物質の運動状態の解析方法並びにその解析装置に関する。

従来より、流体中に含まれる粒子がほぼ球状である場合については、その流体情報（流れ場の分布、流体の粘度等）や粒子情報（粒子の径、初期位置等）に基づき、種々の計算方法でコンピュータ・シミュレーションによって粒子の移動や回転を、比較的精度よく解析することが可能であった。これに対して、流体中に含まれる物質が繊維状物質である場合、流体中における繊維状物質の複雑な形状変化や挙動等から、コンピュータ・シミュレーションにより繊維状物質の移動の状態、変形の状態、分布の状態、配向の状態等を精度良くかつ短時間で解析することが困難であった。そのため、近年では、流体中の繊維状物質の運動状態について精度の高いシミュレーションを行うために様々な方法が検討されてきた。

例えば、特開平５−３１４０９１号公報（特許文献１）においては、流動性のある基質に含まれる粒子の移動及び／又は変形及び／又は配向を、前記基質及び粒子の諸特性、基質の流れ場の分布、及び基質の流れ場における粒子の初期の位置及び配向等の情報に基づきＣＡＥ／ＣＡＤ／ＣＡＭシステムを用いて解析する方法であって、前記粒子の形状に擬似した全体形状となるように複数の球体を１次元、２次元又は３次元方向に相互に結合させた球体集合体を前記粒子の解析用モデルとして用い（球体集合体モデルを用い）、この球体集合体を構成する各球体には隣り合う球体との間に粒子の引張り弾性に相当する結合長の自由度及び／又は粒子の曲げ弾性に相当する結合角の自由度を設定したもとで、各球体の並進及び／又は回転を計算し、これら各球体の計算結果の総計である球体集合体全体の移動及び／又は変形及び／又は配向の計算結果から前記粒子の移動及び／又は変形及び／又は配向を解析する方法が開示されている。このような特許文献１に記載のような解析方法は、繊維状物質の運動状態をシミュレーションするための方法として採用可能であるが、解析結果の精度の点では必ずしも十分なものではなかった。

なお、このようなシミュレーションにより得られる流体中の繊維状物質の運動状態の解析結果は、例えば、繊維強化樹脂（ＦＲＰ）や繊維強化金属（ＦＲＭ）等の複合材料の成形条件、成形用金型等の重要な設計指針となり、また、電気粘性流体や磁性流体等における分散粒子の分布や配向の制御のための有効な資料とし得るため、流体中の繊維状物質の運動状態を解析する方法としては、より精度の高いシミュレーションが可能な方法の出現が望まれている。また、このような流体中の繊維状物質の運動状態を解析する方法としては、従来の方法と比較してより短い時間で、より効率よくシミュレーションすることが可能となるような方法の出現も望まれている。

特開平５−３１４０９１号公報

本発明は、前記従来技術の有する課題に鑑みてなされたものであり、流体中の繊維状物質の運動状態をより短時間でかつより精度高くシミュレーションすることを可能とする流体中の繊維状物質の運動状態の解析方法及びそのような解析方法を実行することが可能な流体中の繊維状物質の運動状態の解析装置を提供することを目的とする。

本発明者らは、前記目的を達成すべく鋭意研究を重ねた結果、下記ステップ（ａ）〜（ｃ）を含み、流体の運動と繊維状物質の運動を連成解析して求められる解析結果を順次利用して、特定の時間が経過するまでの間の流体の運動と繊維状物質の運動を連成解析することにより、流体中の繊維状物質の運動状態をより短時間でかつより精度高くシミュレーションすることが可能となることを見出し、本発明を完成するに至った。

すなわち、本発明の流体中の繊維状物質の運動状態の解析方法は、流体の運動方程式と粒子の運動方程式を利用して繊維状物質及び流体の運動を計算する上で必要となる繊維状物質及び流体の情報に基づいて、流体中の繊維状物質の運動状態をシミュレーションするための解析方法であって、下記ステップ（ａ）〜（ｃ）：
ステップ（ａ）：前記繊維状物質が複数の球体が結合した球体集合体からなるものと仮定する球体集合体モデルを用いる計算方法を用いて、特定の時刻から微小時間経過後における前記繊維状物質の運動を解析するためのステップであり、
流体の速度の情報に基づいて、前記繊維状物質の各球体の流体抵抗力及び前記繊維状物質の各球体が前記流体から受けるトルクを計算し、
前記繊維状物質の各球体間の結合力及び前記繊維状物質の各球体間の曲げ変形に対する復元トルクを計算し、
前記繊維状物質の各球体の流体抵抗力、前記繊維状物質の各球体が前記流体から受けるトルク、前記繊維状物質の各球体間の結合力、及び、前記繊維状物質の各球体間の曲げ変形に対する復元トルクの情報に基づいて、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算する、繊維状物質の運動解析ステップ；
ステップ（ｂ）：前記ステップ（ａ）で計算した前記繊維状物質の各球体の流体抵抗力の情報に基づいて、前記流体が前記繊維状物質から受ける力を計算し、前記流体が前記繊維状物質から受ける力と前記流体中の前記繊維状物質の体積分率とを考慮して、特定の時刻から微小時間経過後における流体の速度を計算する、流体の運動と繊維状物質の運動との連成解析ステップ；
ステップ（ｃ）：前記ステップ（ａ）及びステップ（ｂ）で既に求められた過渡的な流体及び繊維状物質の運動状態の情報を利用して、新たに前記ステップ（ａ）〜（ｂ）の計算を順次実行する処理を、所定条件を満たすまで繰り返し実行する反復演算ステップ；
を含み、前記反復演算ステップ中に流体の運動と繊維状物質の運動を解析して求められる微小時間経過後の過渡的な解析結果を順次利用して、流体の運動と繊維状物質の運動を連成解析し、流体中の繊維状物質の運動状態を解析することを特徴とする方法である。

また、上記本発明の流体中の繊維状物質の運動状態の解析方法においては、前記繊維状物質の運動解析ステップにおいて、繊維状物質間の相互作用力に関する流体潤滑力と機械接触力を計算することなく、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算することが好ましい。

また、本発明の流体中の繊維状物質の運動状態の解析装置は、流体の運動方程式と粒子の運動方程式を利用して繊維状物質及び流体の運動を計算する上で必要となる繊維状物質及び流体の情報に基づいて、流体中の繊維状物質の運動状態をシミュレーションするための解析装置であって、
前記繊維状物質が複数の球体が結合した球体集合体からなるものと仮定する球体集合体モデルを用いる計算方法を用いて、特定の時刻から微小時間経過後における前記繊維状物質の運動を解析するための演算手段であり、流体の速度の情報に基づいて、前記繊維状物質の各球体の流体抵抗力及び前記繊維状物質の各球体が前記流体から受けるトルクを計算し、前記繊維状物質の各球体間の結合力及び前記繊維状物質の各球体間の曲げ変形に対する復元トルクを計算し、前記繊維状物質の各球体の流体抵抗力、前記繊維状物質の各球体が前記流体から受けるトルク、前記繊維状物質の各球体間の結合力、及び、前記繊維状物質の各球体間の曲げ変形に対する復元トルクの情報に基づいて、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算する第一の演算手段と；
前記第一の演算手段により計算された前記繊維状物質の各球体の流体抵抗力の情報に基づいて、前記流体が前記繊維状物質から受ける力を計算し、前記流体が前記繊維状物質から受ける力と前記流体中の前記繊維状物質の体積分率とを考慮して、特定の時刻から微小時間経過後における流体の速度を計算する、第二の演算手段と；
前記第一の演算手段及び前記第二の演算手段で既に求められた過渡的な流体及び繊維状物質の運動状態の情報を利用して、新たに前記第一の演算手段及び前記第二の演算手段での計算を順次実行する処理を、所定条件を満たすまで繰り返す反復演算ステップを実行させるために、該条件を満たしたか否かを判定する、判定手段と；
を備え、前記反復演算ステップ中に流体の運動と繊維状物質の運動を解析して求められる微小時間経過後の過渡的な解析結果を順次利用して、流体の運動と繊維状物質の運動を連成解析し、流体中の繊維状物質の運動状態を解析することを特徴とするものである。

さらに、上記本発明の流体中の繊維状物質の運動状態の解析装置においては、前記第一の演算手段が、繊維状物質間の相互作用力に関する流体潤滑力と機械接触力を計算することなく、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算する演算手段であることが好ましい。

本発明によれば、流体中の繊維状物質の運動状態をより短時間でかつより精度高くシミュレーションすることを可能とする流体中の繊維状物質の運動状態の解析方法及びそのような解析方法を実行することが可能な流体中の繊維状物質の運動状態の解析装置を提供することが可能となる。

本発明の流体中の繊維状物質の運動状態の解析装置の好適な一実施形態を模式的に示すブロック図である。 本発明の流体中の繊維状物質の運動状態の解析方法の好適な一実施形態のフローチャートである。 単純せん断流中の繊維状物質の分散系を模式的に示す概略断面図である。 繊維状物質の好適な一実施形態である細長い直線的な繊維状物質と、その解析用モデルである球体集合体モデルとを模式的に示す模式図である。 球体集合体モデル中の隣り合う２つの球体に引張り応力を加えた前後の状態を模式的に示す模式図である。 球体集合体モデル中の隣り合う２つの球体に曲げ応力（曲げる方向に力）を加えた前後の状態を模式的に示す模式図である。 従来技術（流体中の繊維状物質の運動状態の解析方法の一実施形態）において採用されていた解析方法のフォローチャートである。 実施例１及び比較例１の解析結果であって、せん断応力をかけ始めてからの繊維配向テンソルＡｘｘの時間変化に関するグラフである。 実施例２及び比較例２の解析結果であって、せん断応力をかけ始めてからの繊維配向テンソルＡｘｘの時間変化（時間をせん断速度で無次元化した値，すなわち，ひずみ）に関するグラフである。 実施例１及び比較例１で行った解析の計算時間比並びに実施例２及び比較例２で行った解析の計算時間比を示すグラフである。

以下、図面を参照しながら、本発明の流体中の繊維状物質の運動状態の解析方法及び解析装置の好適な実施形態について詳細に説明する。なお、以下の説明及び図面中、同一又は相当する要素には同一の符号を付し、重複する説明は省略する。

本発明の流体中の繊維状物質の運動状態の解析方法及び本発明の流体中の繊維状物質の運動状態の解析装置は上述の通りである。図１に本発明の流体中の繊維状物質の運動状態の解析装置の好適な一実施形態のブロック図を示す。

図１に示す実施形態の解析装置は、初期条件等を入力するための入力部１と、入力部１によって入力されたデータに基づいて演算を実行するための演算処理部２と、解析結果を出力するための出力部３とを備えるものである。

このような入力部１は、特に制限されず、所望の設定条件等を入力することが可能なものであればよく、例えば、キーボードやマウスなどの他、各種データが予め記憶されたハードディスクやＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）等の記憶媒体であってもよい。

また、出力部３も特に制限されず、演算処理部２による演算結果（解析結果）をデータとして出力できるものであればよく、例えば、データを記憶させる記憶媒体、紙に印刷させるプリンタ、画面に表示させるモニタ等のいずれの形態のものであってもよい。

演算処理部２は、後述する各ステップを実行させるコンピュータプログラムがインストールされたコンピュータによって主に構成される。すなわち、演算処理部２は、後述する各ステップの計算（演算）を実行するためのＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）及びメモリ等からなるハードと、各ステップ（かかるステップについては数式と共に後述する）を実行させるためにインストールされたコンピュータプログラム（ソフト）とを備えるものである。

このようなＣＰＵとしては、例えば、中央処理装置、処理装置、演算装置、プロセッサ、マイクロプロセッサ、マイクロコンピュータ、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）等が挙げられる。

また、このようなメモリとしては、例えば、一時記憶領域としてのＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ Ｐｒｏｇｒａｍｍａｂｌｅ ＲＯＭ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ ＥＰＲＯＭ）等の不揮発性または揮発性の半導体メモリであってもよいし、ハードディスク、フレキシブルディスク等の磁気ディスクであってもよいし、ミニディスク、ＣＤ（Ｃｏｍｐａｃｔ Ｄｉｓｃ）、ＤＶＤ（Ｄｉｇｉｔａｌ Ｖｅｒｓａｔｉｌｅ Ｄｉｓｃ）等の光ディスクであってもよい。なお、このような演算処理部２を構成するメモリには、入力部１により入力される初期条件（基礎情報：例えば、繊維状物質の球体集合体モデルの各球体の初期位置や各球体の初期回転角等の情報等）を記憶する領域が形成されていることが好ましい。

また、演算処理部２は、第一の演算手段２１と、第二の演算手段２２と、判定手段２３とを備える。

このような第一の演算手段２１は、前記繊維状物質が複数の球体が結合した球体集合体からなるものと仮定する球体集合体モデルを用いる計算方法を用いて、特定の時刻から微小時間経過後における前記繊維状物質の運動を解析するための演算手段であり、流体の速度の情報に基づいて、前記繊維状物質の各球体（前記繊維状物質の球体集合体モデルの各球体）の流体抵抗力及び前記繊維状物質の各球体が前記流体から受けるトルクを計算し、前記繊維状物質の各球体間の結合力及び前記繊維状物質の各球体間の曲げ変形に対する復元トルクを計算し、前記繊維状物質の各球体の流体抵抗力、前記繊維状物質の各球体が前記流体から受けるトルク、前記繊維状物質の各球体間の結合力、及び、前記繊維状物質の各球体間の曲げ変形に対する復元トルクの情報に基づいて、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角（角度）を計算するための演算手段（例えば、メモリに格納した、前記計算処理を実行するためのプログラムによりＣＰＵに演算（計算）処理を実行させるための手段）である。

なお、このような演算に利用する流体の速度の情報としては、既に第二の演算手段により計算された計算結果の情報（データ）がある場合にはその計算結果を利用し、また、未だ第二の演算手段により流体の速度を計算していない場合（１回目の計算時）においては、初期条件として入力された流体の速度の情報を利用してもよいし、あるいは、初期条件として入力された情報に基づいて計算して求めてもよい。また、このような球体集合体モデルを用いる計算方法により、流体の速度の情報に基づいて、前記繊維状物質の各球体の流体抵抗力及び前記繊維状物質の各球体が前記流体から受けるトルクを計算する方法、並びに、前記繊維状物質の各球体間の結合力及び前記繊維状物質の各球体間の曲げ変形に対する復元トルクを計算する方法の具体的な方法は後述する。なお、このような繊維状物質の各球体間の結合力及び前記繊維状物質の各球体間の曲げ変形に対する復元トルクを計算する方法としては公知の方法（例えば特開平５−３１４０９１号公報に記載のような方法）をそのまま応用してもよい。

また、第一の演算手段２１においては、特定の時刻から微小時間経過後（時刻：ｔ＋Δｔ）における前記繊維状物質の運動を解析する。本発明においては、このような解析において、上記計算により求められる前記繊維状物質の各球体の流体抵抗力、前記繊維状物質の各球体が前記流体から受けるトルク、前記繊維状物質の各球体間の結合力、及び、前記繊維状物質の各球体間の曲げ変形に対する復元トルクの情報をそのまま用いて、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角（曲げ角度）を計算する。そして、このような計算の際には、繊維状物質間の相互作用力に関する流体潤滑力と機械接触力を計算することなく、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算することが好ましい。このように、第一の演算手段２１は、繊維状物質間の相互作用力に関する流体潤滑力と機械接触力を計算することなく、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算する演算手段であることが好ましい。このように、本発明においては、第一の演算手段２１において、繊維状物質同士の２体相互作用（繊維状物質間の相互作用力に関する流体潤滑力と機械接触力）を計算することなく、特定の時刻から微小時間経過後における前記繊維状物質の各球体（前記繊維状物質の球体集合体モデルの各球体）の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体（前記繊維状物質の球体集合体モデルを構成する球体であって隣り合う球体）間の曲げ角を計算することができる。そのため、本発明によれば、より短時間で運動状態をシミュレーションすることが可能となる。なお、具体的な計算方法等については後述する。

また、第二の演算手段２２は、第一の演算手段２１により、計算された前記繊維状物質の各球体の流体抵抗力の情報に基づいて、前記流体が前記繊維状物質から受ける力を計算し、前記流体が前記繊維状物質から受ける力と前記流体中の前記繊維状物質の体積分率とを考慮して、特定の時刻から微小時間経過後における流体の速度を計算するための演算手段（例えば、メモリに格納した、前記計算処理を実行するためのプログラムによりＣＰＵに演算（計算）処理を実行させるための手段）である。このような計算により、繊維状物質が流体から受ける力と、流体が前記繊維状物質から受ける力とを連成解析して、流体の速度を計算することが可能である。そして、本発明においては、後述の判定手段２３により、所定条件を満たすものと判定されるまで、微小時間経過後の過渡的な解析結果を順次利用して、新たに第一の演算手段２１による計算及び第二の演算手段２２による計算を実行する処理を繰り返す。これにより、流体の運動と繊維状物質の運動を連成解析して求められる微小時間経過後の過渡的な解析結果（例えば、流体の速度、流体中の繊維状物質の各球体の位置、速度及び角速度、前記繊維状物質を構成する隣り合う球体間の曲げ角等の変位の情報等）を順次利用して、流体中の繊維状物質の運動状態を解析することが可能となる。なお、このような解析方法においては、繊維状物質同士の相互作用を計算することがないため、より短時間で運動状態をシミュレーションすることが可能となる。また、このような解析方法によれば、流体の運動と繊維状物質の運動を連成解析しているため、より精度の高いシミュレーションを行うことが可能となる。

判定手段２３は、前記第一の演算手段及び前記第二の演算手段で既に求められた過渡的な流体及び繊維状物質の運動状態の情報を用いて（例えば第一の演算手段２１で用いる前記流体の速度として第二の演算手段で求められる流体の速度を用いて）、第一の演算手段２１及び第二の演算手段２２での計算を順次実行する処理を、所定条件を満たすまで繰り返し実行させるために、該条件を満たしたか否かを判定するための手段である。このような条件の判定方法としては、例えば、前記所定条件として繰り返し計算する回数を設定して、繰り返し計算する回数が設定した回数を満たしたか否かを判定してもよく、あるいは、第一の演算手段２１及び第二の演算手段２２での計算が、所定の時刻ｔから微小時間（Δｔ）追加した時間における計算結果を求めるものであるため、前記所定条件として最終的な時間（Ｔ_ｅｎｄ）を設定して、第一の演算手段２１及び第二の演算手段２２での計算を順次実行する処理が終了する度に、計算結果を求めた時間（ｔ＋Δｔ）が設定時間（Ｔ_ｅｎｄ）よりも小さいか否かを判断することで、所定条件を満たすか否かを判断してもよい。なお、このような判定は、例えば、メモリに格納した判定プログラムに条件を入力してＣＰＵに実行させればよい。また、本発明においては、判定手段２３で条件を判定した後に、所定条件を満たしていなかった場合には、第一の演算手段２１及び第二の演算手段２２での計算を順次実行する処理を、所定条件を満たすまで繰り返し実行させる。このような観点から、判定手段２３は、条件に応じて反復演算を実行させるための手段としても機能させることが好ましい。なお、本実施形態における演算処理部２には、判定手段２３において所定条件を満たさないものと判断した場合に、第一の演算手段２１及び第二の演算手段２２による計算結果を利用して、新たに第一の演算手段２１及び第二の演算手段２２での計算を順次実行させるためのコンピュータプログラム；判定手段２３において所定条件を満たすものと判断した場合に、それまでの計算結果を出力するためのコンピュータプログラム；等をメモリに格納してもよい。

以下、図２に示すフローチャート（本発明の流体中の繊維状物質の運動状態の解析方法の好適な一実施形態）に基づいて、本発明の流体中の繊維状物質の運動状態の解析方法の好適な実施形態、及び、上記演算手段による計算の具体的な計算方法を説明する。なお、このようなフローチャートに記載するような解析処理方法は、ユーザーが入力部１を介して解析処理を行うよう指示した時に、演算処理部２のメモリ中に記憶された解析処理プログラムを実行させることで開始させてもよい。また、以下に説明する本発明の流体中の繊維状物質の運動状態の解析方法の好適な実施形態は、せん断力が加えられた場合（せん断場中）の流体中の繊維状物質の運動状態（例えば、流体中における繊維状物質の配向状態）を解析する方法である。

このような解析に際しては、先ず、ステップＳ１において、入力部１を介して演算処理部２に初期条件（基礎情報）を入力する。

このような初期条件（基礎情報）としては、流体の運動方程式と粒子の運動方程式を利用して繊維状物質及び流体の運動を計算する上で必要となる繊維状物質及び流体の情報であり、例えば、流体の計算に必要となる計算メッシュ（計算格子：計算セル）の情報、繊維状物質を後述の球体集合体モデルにモデル化するために必要な繊維状物質の物性等の情報、流体の運動を解析するために必要な流体の物性等の流体の境界条件の情報、等が挙げられる。なお、このような「流体の運動方程式」は、連続の式（質量保存）およびナビエ・ストークス方程式（Navier-Stokes方程式）の他、流体の運動を計算するために、及び、繊維状物質の運動及び流体の運動を連成して計算するために、必要となる後述の各種計算式等を含む概念である。また、「粒子の運動方程式」は、式（７）又は式（７’）、及び、式（８）の方程式など、繊維状物質の運動を計算するために、及び、繊維状物質の運動及び流体の運動を連成して計算するために、必要となる後述の各種計算式等を含む概念である。

このような計算メッシュに関して、流体中に繊維状物質が分散された状態を解析するために、適切な範囲に解析領域を有限個に分割すればよく、そのような分割を行うための条件等は解析する対象に応じて公知の方法を採用して適宜設定すればよい。

また、初期条件（基礎情報）として入力する前記繊維状物質の情報としては、例えば、繊維状物質の本数、繊維状物質の直径、繊維状物質のアスペクト比、繊維状物質の密度（ρ^ｐ）、繊維状物質のヤング率の情報が挙げられる。

また、初期条件（基礎情報）として入力する前記流体の情報としては、例えば、流体の密度（ρ^ｆ）、粘度係数（μ）、流体中の圧力（ｐ）、流体の速度（壁面速度や、ｕ^ｆ _ｉ：初期条件としてのｔ_ｎ（ｎ＝０）の時の速度（繊維状物質は未考慮））等の境界条件、計算領域の大きさなどの情報が挙げられる。なお、流体中の圧力（Ｐ）等については、例えば、非圧縮性流体を仮定した計算法の場合は流体の質量が保存されるように圧力が決定されるため、特別な初期条件は必要としない。また、前記基礎情報として利用する流体の情報のうち、流体の運動に関する情報（速度や角速度）を入力する場合には、連続の式（質量保存）およびナビエ・ストークス方程式（Navier-Stokes方程式）に基づいて解析した値を利用してもよい。なお、流体の運動に関する情報（速度や角速度）は、流体の特性の情報や境界条件等の情報に基づいて、演算処理部２のＣＰＵで計算させてもよい。ここで、例えば、図３に示すようなせん断流動場（単純せん断流中の繊維状物質の分散系：なお、図３中のＵ_０は、移動するプレート（板）の速度を示し、２δは計算領域の１辺の長さを示す）を例に挙げると、流体の速度（ｕ^ｆ _ｉ）は、そのＵ_０やδ（delta）の情報（せん断速度（＝Ｕ_０／δ）の情報）に基づいて、連続の式（質量保存）およびナビエ・ストークス方程式（Navier-Stokes方程式）に基づいて解析できる。

このような初期条件としては、少なくとも、繊維状物質の物性値の情報（繊維状物質を球体集合体モデルにモデル化するために必要な情報：繊維状物質の直径、アスペクト比、ヤング率など）、繊維状物質の初期位置及び速度の情報、流体の物性値及び境界条件（壁面上のすべりなし条件（非滑り境界条件）、周期境界条件（計算メッシュの情報）、流体の初期速度）を入力する必要がある。

なお、このような初期条件の入力に際しては、例えば、繊維強化樹脂の射出成形を想定する場合には、ストークス数が、概ね１以下となる条件（１に近似した値（ほぼ１）か又はそれよりも小さな値となる条件）となるように、せん断速度や流体の粘度などを設定することが好ましい。このように設定した値を入力することにより繊維状物質に加わる慣性力よりも繊維状物質に加わる流体の粘性力の方が支配的になり、上記成形過程に即した解析を数値計算上、安定に計算することが可能となる。なお、このようなストークス数の条件は特に制限されず、想定する対象における繊維状物質に加わる慣性力と繊維状物質に加わる粘性力のバランスに応じて、その設定値を変更してもよい。

このように、ステップＳ１は、初期条件（基礎情報）として、流体の運動方程式と粒子の運動方程式を利用して繊維状物質及び流体の運動を計算する上で必要となる繊維状物質及び流体の情報（前記繊維状物質の物性値の情報、前記繊維状物質の初期位置及び速度の情報、前記流体の物性値及び境界条件）を入力するステップである。すなわち、ステップＳ１は、後述の計算（球体集合体へのモデル化を含む）をするために必要となる繊維状物質及び流体の初期条件（基礎情報）を入力するステップである。

なお、このように入力部１により入力する初期条件（基礎情報）に関して、例えば、前述のように例示した情報のうちの一部の情報は、入力した他の情報に基づいて演算処理部２のＣＰＵに別途計算させてもよいため、上記情報の中からその解析プログラムの種類等に応じて、後述の計算を行うために必要となる情報を適宜入力すればよい。なお、このような初期条件（基礎情報）及び／又はＣＰＵに別途計算させた初期の状態は入力後、演算処理部２中のメモリの適切な領域に記憶させる等して、計算時に適宜その情報を利用できるようにすればよい。

また、図２に示すフローチャートに記載のような解析方法を実行する際には、ステップＳ１において入力部１を介して入力された初期条件（基礎情報）を利用して、演算処理部２のＣＰＵにおいて繊維状物質１０の球体の集合体によるモデル化が行われる。以下、このような球体集合体モデルについて図面を参照しながら簡単に説明する。

図４は、繊維状物質の好適な一実施形態である細長い直線的な繊維状物質１０と、その解析用モデルである球体集合体モデル１０Ａとを模式的に示す模式図である。このような図４に示す実施形態の球体集合体モデル１０Ａは、５個の直径ｄ_Ｐの球体が直線状に連なった状態を示すものである。このような球体集合体モデル（Bead model）のモデル化を演算処理部２のＣＰＵにより行う場合、例えば、初期条件（基礎情報）として入力された繊維状物質の直径やアスペクト比の情報等に基づいて、球体の直径や個数等を適宜計算してモデル化してもよく、或いは、予めサイズなどを設定した球体により球体集合体を形成するようにしてモデル化してもよい。なお、図４に示す実施形態では、繊維形状の物質１０の直径ｄと、球体集合体モデル１０Ａの各球体の直径ｄ_Ｐを同じ大きさとして、５個の直径ｄ_Ｐの球体が直線状に連なった球体集合体モデルを繊維状物質１０としてモデル化している。

また、このようなモデル化の演算に際しては、例えば、図４に示すように繊維状物質１０が直線的なものである場合、ＣＰＵにより、球体集合体モデルも複数の球体が直線的に結合された集合体１０Ａとしてモデル化する。そして、かかる球体集合体モデル１０Ａを用いて、初期の球体集合体モデルの各球体の中心Ｏ１〜Ｏ５の位置（ｘ、ｙ、ｚ方向の位置情報等）や各球体間の結合角等といった初期の位置情報等は、シミュレーションする対象等応じて適宜設定できる。ここで、隣り合う球体の中心間距離の変化量、各球体の結合角（曲げ角）といった概念を説明するために、図５〜図６に、球体集合体モデル１０Ａ中（５つの球体中）の隣り合う２つの球体の状態を模式的に示す。

図５は、球体集合体モデル１０Ａ中の隣り合う２つの球体に引張り応力を加えた前後の状態を模式的に示すものであり、図５中の（ａ）は引張り応力を加える前の状態（図４に示すような直線状に並んだ５つの球体のうちの任意の２つの隣り合う球体の状態）を模式的に示し、図５中の（ｂ）は引張り応力が加えられた後（微小時間経過後）の状態を模式的に示す。なお、このような球体集合体モデル１０Ａにおいては、引張り応力が加えれらて球体間の距離が伸びた場合に弾性的に元に戻ろうとする復元力（結合力）があるものと仮定して、繊維状物質の状態を解析する。そして、このような球体集合体モデル１０Ａの各球体間の結合力（球体間の距離が伸びた時に元に戻ろうとする力）としては、球体の中心間距離の変化量（σ_ｉｊ）と、繊維状物質のヤング率（Ｅ）とに基づいて求められる値を採用する（具体的な計算方法は後述する）。ここにおいて、球体の中心間距離の変化量（σ_ｉｊ）は、応力が加えられた後の球体の中心間距離ｒ_１から、応力が加えられる前の球体の中心間距離ｒ_０を引いた値（ｒ_１−ｒ_０）として求めることができる。

また、図６は繊維状物質を曲げる方向に力を加えた場合について、球体集合体モデル１０Ａ中（５つの球体中）の隣り合う２つの球体の力を加える前後の状態を模式的に示すものであり、図６中の（ａ）は繊維状物質を曲げる方向に力が加えられる前の状態（図４に示すような直線状に並んだ５つの球体のうちの任意の２つの隣り合う球体の状態）を模式的に示し、図６中の（ｂ）は繊維状物質を曲げる方向に力が加えられた後の状態を模式的に示す。本発明においては、球体集合体モデル１０Ａを利用して、前記繊維状物質の各球体間の曲げ変形に対する復元トルク（繊維状物質を曲げた時に元に戻ろうとする力）を計算する際に、結合角（曲げ角：隣あう二つの球体間の角度）θ_ｉｊに基づいて計算を行う（具体的な計算方法は後述する）。

このように、本発明においては、初期情報として与えられた情報に基づいて、繊維状物質１０をモデル化する。また、このようなモデル化に際しては、繊維状物質１０を構成する各球体の中心Ｏ１〜Ｏ５の初期の位置（中心初期位置：中心のｘｙｚ座標）と、各球体ごとに初期の回転角θｉｊを演算し、かつ、演算して求められた値をメモリの適切な領域に記憶して利用することができ、そのようなモデル化等の演算をするようにプログラミングしたプログラムを適宜利用してもよい。なお、このようなプログラムによる演算は、その演算を行うためのモデル化のための演算手段（図示せず）を更に備えるものとして演算処理部２を構成させて達成してもよいし、第一の演算手段２１により、各種演算を行う前に併せて計算させてもよい。これにより、初期条件（基礎情報）の入力により、演算処理部２において初期の状態を演算させることが可能となる。なお、このような球体集合体モデルへのモデル化の方法としては、特開平５−３１４０９１号公報に記載の公知の方法で採用している方法を利用してもよい。

また、ステップＳ１において、入力部１を介して演算処理部２に初期条件（基礎情報）を入力した後においては、ステップＳ２において、時刻ＴをΔｔだけ更新する。このようなステップにより、第一の演算手段２１において前記繊維状物質の各球体（前記繊維状物質を構成する各球体）の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角の計算値の時刻Ｔをｔ_ｎ＋Δｔ（なお、ここで、ｎは演算の反復回数であり、初回は０となる）とすることができる。

次に、ステップＳ３においては、第一の演算手段２１を利用し、流体の速度（ｕ^ｆ _ｉ）の情報に基づいて、特定の時刻ｔ_ｎ（なお、時刻ｔ_ｎはＴ−Δｔで表される時刻：Δｔ更新する前の時刻）における繊維状物質の流体抵抗力（ｆ^ｈ _ｉ）及び特定の時刻ｔ_ｎにおける流体から受けるトルク（Ｔ^ｈ _ｉ）を計算する。そして、このような流体抵抗力（ｆ^ｈ _ｉ）及び復元トルク（Ｔ^ｈ _ｉ）は、下記式（１）〜（２）：

［式（１）〜（２）中、ｆ^ｈ _ｉは球状集合体モデルの球体の流体抵抗力を示し、Ｔ^ｈ _ｉは球状集合体モデルの球体の復元トルク（各球体にかかる流体からのトルク）を示し、μは流体の粘度を示し、ｄ_Ｐは前記球状集合体モデルの球体の直径を示し、ｕ^ｆ _ｉは流体の速度を示し、ｕ^Ｐ _ｉは前記球状集合体モデルの球体の速度を示し、ω^ｆ _ｉは流体の角速度を示し、ω^Ｐ _ｉは球状集合体モデルの球体の角速度を示す。なお、球体に関する各記号に関して、下付きの文字ｉは各球体に付した番号であり、計算領域内にＮ個の球体が存在する場合、かかるｉはｉ≦Ｎの条件を満たす整数となる。また、流体に関する各記号に関して、下付きの文字ｉは、デカルト座標の各方向の成分（例えばｘｙｚ等）であることを意味する。］
で表される計算式により求めることができる。かかる計算を行うステップＳ３は、前記ステップ（ａ）中の一ステップに相当し得るステップであり、流体の速度の情報に基づいて、前記繊維状物質の各球体の流体抵抗力（例えば、Stokes抵抗力など）及び前記繊維状物質の各球体が前記流体から受けるトルクを計算するステップである。このような計算は、特開平５−３１４０９１号公報の段落［００２４］〜［００２５］に記載の計算と基本的に同様のものである。

なお、このような計算に利用する流体の速度（ｕ^ｆ _ｉ）に関して、１回目の計算時（第二の演算手段２２による計算結果が得られる前の段階）においては、前記初期条件として速度（予め別途計算した値）の情報が入力されている場合（例えば、予めナビエストークス方程式に基づいて計算した値などが入力されている場合）にはそれを利用してもよい。また、このような計算に利用する流体の速度としては、１回目の計算時（第二の演算手段２２による計算結果がでる前の段階）において、前記初期条件として予め求めた速度（予め別途計算した値）の情報（データ）がない場合には、入力されている初期条件（流体に関する情報）に基づいて、流体中に前記繊維状物質が存在しないものと仮定して初期時刻（ｔ_ｎ：ｎ＝０）における前記流体の速度を別途計算して利用してもよい。このような計算を行う場合、ナビエ・ストークス方程式に基づいて計算する方法を採用できる。なお、１回目の計算時（第二の演算手段２２による計算結果が得られる前の段階）においては、後述の式（１１１）〜（１１２）を利用して計算した流体の速度を利用してもよい。

また、このような計算に利用する球状集合体モデルの球体の速度（ｕ^Ｐ _ｉ）は、１回目の計算時（第二の演算手段２２による計算結果が得られる前の段階）においては、ｕ^Ｐ _ｉ＝０等の初期値を採用してもよく、または、上記のように流体中に繊維状物質が存在しないと仮定して別途計算した流体速度を用いて、その流体速度を各球体の位置において空間補正した値を採用してもよく、２回目以降の計算時（反復実行ステップにおける計算時）には、後述のステップＳ５において計算した値を利用できる。

また、このような計算に利用する流体の角速度（ω^ｆ _ｉ）は、１回目の計算時（第一の演算手段２１による計算結果がでる前の段階）においては、上述のようにして用いた１回目の計算時の流体の速度（ｕ^ｆ _ｉ）の情報から計算できる値を利用でき、また、２回目以降の計算時（反復実行ステップにおける計算時）においては、後述のステップＳ５において計算した流体の速度の情報から計算できる値を利用できる。なお、このような角速度の具体的な算出方法としては、例えば、下記式（Ｉ）：

［式（Ｉ）中、ω^ｆ _ｉは角速度であり、ｘ_ｊ及びｘ_ｋはデカルト座標を示し、ｕ^ｆ _ｋ及びｕ^ｆ _ｊは、デカルト座標の各方向の速度であることを意味する。］
で表される計算式を計算する方法を挙げることができる。なお、上記式（Ｉ）を各方向の成分で検討すると、下記式：

で表される計算式となる。

また、このような計算に利用する球状集合体モデルの球体の角速度（ω^Ｐ _ｉ）は、１回目の計算時においては、初期条件として入力した値又はその入力値から計算した値を採用する。なお、かかる角速度の値（ω^Ｐ _ｉ）は、１回目の計算時においては、ω^Ｐ _ｉ＝０などの初期値を採用してもよく、あるいは、上記のように流体中に繊維状物質が存在しないと仮定して別途計算した流体速度から算出できる流体角速度を用いて、その流体角速度を各球体の位置において空間補正した値を利用してもよく、２回目以降の計算時（反復実行ステップにおける計算時）には、後述のステップＳ５において計算した値を利用できる。このような値を利用して計算することにより、繊維状物質の各球体の流体抵抗力（ｆ^ｈ _ｉ）、繊維状物質の各球体が流体から受けるトルク（Ｔ^ｈ _ｉ）を求めることができる。

次に、ステップＳ４においては、第一の演算手段２１を利用し、時刻ｔ_ｎ（なお、時刻ｔ_ｎはＴ−Δｔで表される時刻：Δｔ更新する前の時刻）における繊維状物質の各球体間の結合力（Ｆ^ｓ _ｉ）及び時刻ｔ_ｎにおける繊維状物質の各球体間の曲げ変形に対する復元トルク（Ｔ^ｂ _ｉ）を計算する。このような各球体間の結合力（Ｆ^ｓ _ｉ）及び復元トルク（Ｔ^ｂ _ｉ）は、下記式（３）〜（６）：

［式（３）〜（６）中、Ｆ^ｓ _ｉは各球体間の結合力（球体集合体モデル中の隣り合う２つの球体間に生じる引張り変形や圧縮変形に対する復元力）を示し、ｋ_ｓは引張り・圧縮変形に対する定数を示し、Ｅはヤング率（縦弾性係数）を示し、ｄ_ｐは球体集合体モデル中の球体の直径を示し、Ｔ^ｂ _ｉは繊維状物質の各球体間の曲げ変形に対する復元トルク（球体集合体モデル中の隣り合う２つの球体間に生じる曲げ変形に対する復元トルク）を示し、ｋ_ｂは曲げ変形に対する定数を示し、σ_ｉｊは球体の中心間距離の変化量を示し（図５参照）、θ_ｉｊは曲げ角（結合角：球体集合体モデル中の隣り合う２つの球体間の角度）を示す（図６参照）。なお、各記号に関して、下付きの文字ｉ、ｊは、例えば、図４に示すような球体集合体モデルで考慮すると、５個のうちのｉ番目の球体及びその隣のｊ番目の球体に関する値であることを意味する（図４に示すモデルを利用する場合にはｉ及びｊはそれぞれ５以下の任意の整数となる）。］
で表される計算式により求めることができる。なお、このような計算式（３）〜（６）は、各球体について、縦弾性係数（Ｅ）を利用して、その変形量に応じた復元力及び復元トルクを求めるための方程式である。かかるステップは、前記ステップ（ａ）中の一ステップに相当し得るステップであり、前記繊維状物質の各球体間の結合力及び前記繊維状物質の各球体間の曲げ変形に対する復元トルクを計算するステップである。

また、σ_ｉｊとしては、時刻ｔ_ｎにおける各球体の位置情報（座標）から求められる値（１回目の計算時には初期条件と入力した値）を利用でき、θ_ｉｊとしては、時刻ｔ_ｎにおける各球体間の曲げ角の値（１回目の計算時には各球体の初期条件から求められる値を採用できる）を利用できる。

次いで、ステップＳ５において、第一の演算手段２１を利用し、前記繊維状物質の各球体の位置（例えばｘｙｚ座標等）、速度（ｕ^Ｐ _ｉ）及び角速度（ω^Ｐ _ｉ）、並びに、前記繊維状物質を構成する隣り合う球体間の曲げ角（θ_ｉｊ）を計算する。このようなステップＳ５における計算は、ステップＳ３〜４において計算した、前記繊維状物質の各球体の流体抵抗力、前記繊維状物質の各球体が流体から受けるトルク、前記繊維状物質の各球体間の結合力、及び、前記繊維状物質の各球体間の曲げ変形に対する復元トルクの情報に基づいて行う。このような各球体の速度（ｕ^Ｐ _ｉ）及び角速度（ω^Ｐ _ｉ）は、下記式（７）〜（１０）：

［式（７）〜（１０）中、ｍ^Ｐは球体の質量を示し、ｕ^Ｐ _ｉは球体ｉの速度を示し、ｕ^P _ｊは球体ｊの速度を示し、ｆ^ｈ _ｉは球状集合体モデルの球体の流体抵抗力を示し、Ｆ^ｓ _ｉは繊維状物質の各球体間の結合力を示し、Ｔ^ｂ _ｉは繊維状物質の各球体間の曲げ変形に対する復元トルクを示し、Ｔ^ｈ _ｉは繊維状物質の各球体が流体から受けるトルクを示し、ｆ_ｉｊは隣り合う２つの粒子間に生じる摩擦力を示し、ｄ_ｐは球体集合体モデル中の球体の直径を示し、ω^Ｐ _ｉは球体ｉの角速度を示し、ω^Ｐ _ｊは球体ｊの角速度を示し、θ_ｉｊは球状集合体モデルの球体ｉと球体ｊの間の曲げ角（結合角）を示し、ｔは時間を示し、ｎ_ｉｊは球状集合体モデルの球体ｉ及び球体ｊのそれぞれの中心間を結ぶ球体ｉから球体ｊに向かう方向の単位ベクトルを示し、ｎ_ｊｉは球状集合体モデルの球体ｉ及び球体ｊのそれぞれの中心間を結ぶ球体ｊから球体ｉに向かう方向の単位ベクトルを示す。なお、ここにいう球体ｉ及び球体ｊにおけるｉ及びｊは、例えば、図４に示すような球体集合体モデルで考慮すると、５個の球体のうちのｉ番目の球体及びその隣のｊ番目の球体であることを意味する（図４に示すモデルを利用する場合にはｉ及びｊはそれぞれ５以下の任意の整数となる）。］
で表される計算式により求めることができる。

ここで、隣り合う２つの球体間に生じる摩擦力（ｆ_ｉｊ）の算出方法は、特に制限されず、公知の方法を適宜採用でき、例えば、上記式（１０）の両辺を時間微分し、上記式（７）及び（８）を代入して摩擦力（ｆ_ｉｊ）について整理した式を、ＳＯＲ（Successive Over-Relaxation）法などの収束解法を用いて解くことにより算出する方法を採用してもよい。すなわち、このような球体間摩擦力（ｆ_ｉｊ）の算出方法としては、特に制限されないが、例えば、上記式（１０）の両辺を時間微分し、上記式（７）及び（８）を代入して求められる下記式（ＩＩ）：

［なお、式（ＩＩ）中の左辺及び右辺の各成分はそれぞれ下記式：

（式中のｉ、ｊ、ｋ、ｌは、繊維状物質を構成する隣り合う球体の番号を示す。）
で表されるものである。］
を、ＳＯＲ（Successive Over-Relaxation）法などの収束解法を用いて解くことにより算出する方法を採用してもよい。なお、このような隣り合う２つの球体間に生じる摩擦力（ｆ_ｉｊ）の算出方法は、上記した算出方法（例示した方法）に限定されるものではなく、公知の他の算出方法を適宜採用してもよい。

このような計算式（７）及び（８）はそれぞれ球体集合体モデルの各球体の並進運動と回転運動の方程式であり、計算式（９）は繊維状物質を構成する隣り合う球体ｉと球体ｊ間の曲げ角の方程式であり、計算式（１０）は球体が連結点において滑らないことを示す条件（拘束条件）を与える方程式であり、これらを連立させて計算することで、各球体の速度（ｕ^Ｐ _ｉ）及び角速度（ω^Ｐ _ｉ）を求めることができる。そして、このような各球体の速度（ｕ^Ｐ _ｉ）及び角速度（ω^Ｐ _ｉ）に基づいて、球体集合体モデルの各球体の中心Ｏ１〜Ｏ５の移動する位置を算出できる。このようにして、特定の時刻（ｔ_ｎ）から微小時間経過後（時刻：Ｔ＝ｔ_ｎ＋Δｔ）における前記繊維状物質の各球体の位置、速度及び角速度、並びに、特定の時刻（ｔ_ｎ）から微小時間経過後（時刻：Ｔ＝ｔ_ｎ＋Δｔ）における前記繊維状物質を構成する隣り合う球体間の曲げ角を求めることができる。また、このような球体集合体モデルの各球体の運動を求めることで、各球体の中心Ｏ１〜Ｏ５の初期の位置から移動する位置を求めることができ、これにより繊維状物質の配向状態を求めることができる。ただし、曲げ角に関する計算方法は、上記式（９）に示す計算式に限定されるものではなく、四元数などの演算を用いて算出してもよい。

また、このような特定の時刻（ｔ_ｎ）から微小時間経過後（時刻：Ｔ＝ｔ_ｎ＋Δｔ）における前記繊維状物質の各球体の位置、速度及び角速度、並びに、前記繊維状物質を構成する隣り合う球体間の曲げ角を計算する際には、上記式（７）〜（１０）を採用して、繊維状物質間の相互作用力に関する流体潤滑力と機械接触力を計算することなく、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算することが好ましい。本発明においては、上述のように、繊維状物質間の相互作用力に関する流体潤滑力と機械接触力の計算を省略することが可能であり、これにより計算時間を大幅に短縮することが可能となる。なお、本発明において、計算を省略することが可能な前記繊維状物質間の相互作用力（流体潤滑力及び機械接触力）の計算式等については、後述の従来の流体中の繊維状物質の運動状態の解析方法と併せて説明する。

このようなステップＳ５は、前記ステップ（ａ）中の一ステップに相当し得るステップであり、前記繊維状物質の各球体の流体抵抗力、前記繊維状物質の各球体が前記流体から受けるトルク、前記繊維状物質の各球体間の結合力、及び、前記繊維状物質の各球体間の曲げ変形に対する復元トルクの情報に基づいて、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度、並びに、前記繊維状物質を構成する隣り合う球体間の曲げ角を計算するステップである。このように、式（７）〜（１０）中のｆ^ｈ _ｉ、Ｆ^ｓ _ｉ、Ｔ^ｂ _ｉ、Ｔ^ｈ _ｉの値としては、前記ステップＳ３及びＳ４で求められた値を利用する。

このように、ステップＳ２〜Ｓ５は、上記計算内容を考慮すれば、繊維状物質の運動解析を行うステップ（繊維状物質の運動解析ステップ）であるといえ、かかるステップにより、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置（座標）、速度、角速度、並びに、前記繊維状物質を構成する隣り合う球体間の曲げ角を解析することができる。

次に、ステップＳ６において、第二の演算手段２２を利用し、前記流体が前記繊維状物質から受ける力（Ｆ^ｆ _ｉ）を計算する。このような繊維状物質から受ける力（Ｆ^ｆ _ｉ）の計算は、前記ステップＳ３で計算した前記繊維状物質の各球体の流体抵抗力（ｆ^ｈ _ｉ）の情報に基づいて行う。このような流体が前記繊維状物質から受ける力（Ｆ^ｆ _ｉ）は、下記式（１１）：

［式（１１）中、Ｆ^ｆ _ｉは流体が前記繊維状物質から受ける力を示し、Ｖ_ｐは計算セル１つの体積を示し、ｆ^ｈ _ｉは繊維状物質の各球体の流体抵抗力（流体に対する抵抗力）を示し、式（１１）中の下記部分：

は、Ｖｐ内に存在する全球体について総和計算することを示す。］
で表される計算式により求めることができる。このような計算において利用するＶ_ｐは計算セル１つの体積であり、計算領域の大きさと計算セルの分割数から求めることができる。また、このような計算式（１１）は、流体が前記繊維状物質から受ける反力であり、１つの計算メッシュ（計算セル）中のｆ^ｈ _ｉの総和を求める式である。このようにして求められる流体が前記繊維状物質から受ける力（Ｆ^ｆ _ｉ）は、いわゆる繊維状物質が流体から受けるストークス抵抗力に対する反力であるといえる。このようなステップＳ６は、前記ステップ（ｂ）中の一ステップに相当し得るステップであり、前記ステップ（ａ）で計算した前記繊維状物質の各球体の流体抵抗力の情報（ステップＳ３で計算した繊維状物質の各球体の流体抵抗力（ｆ^ｈ _ｉ）の情報）に基づいて、前記流体が前記繊維状物質から受ける力を計算するステップに相当し得る。

次いで、ステップＳ７において、第二の演算手段２２を利用し、特定の時刻から微小時間経過後における流体の速度（ｕ^ｆ _ｉ）計算する。このような流体の速度（ｕ^ｆ _ｉ）の計算は、前記流体が繊維状物質から受ける力（Ｆ^ｆ _ｉ）と、前記流体中の前記繊維状物質の体積分率（ε）とを考慮した下記式（１２）〜（１３）：

［式（１２）〜（１３）中、ｕ^ｆ _ｉは時刻ｔ_ｎから微小時間（Δｔ）経過後（時刻：Ｔ＝ｔ_ｎ＋Δｔ）における流体の速度を示し、ｕ^ｆ _ｊは時刻ｔ_ｎにおける流体の速度を示し、εは流体中の繊維状物質の体積分率を示し、ｔは時間を示し、ｘ_ｉ及びｘ_ｊはデカルト座標（ｘ，ｙ，ｚ等）を示し、ρ^ｆは流体の密度（流体密度）を示し、Ｐは流体中の圧力を示し、Ｆ^ｆ _ｉは流体が繊維状物質から受ける力を示し、μは流体の粘性係数を示す。なお、下付きの文字ｉ及びｊはアインシュタインの総和規約にしたがうものとする。］
で表される計算式により求めることができる。このような計算式（１２）及び（１３）は、流体が繊維状物質から受ける力（反力）と、流体中の繊維状物質の体積分率とを考慮し、流体の運動と繊維状物質の運動の相互作用（繊維が流体から受ける力と流体が繊維状物質から受ける力（反力））を考慮した解析（連成解析）によって、流体の速度を求める計算式である。

なお、このような流体の体積分率（ε）は時刻ｔ_ｎにおける各計算セルの体積と、それぞれのセル内に存在する球体（繊維状物質を構成する各球体）の体積の総和の比（［球体の体積の総和］／［計算セルの体積］）から求めることができる。

このようなステップＳ７は、前記ステップ（ｂ）中の一ステップに相当し得るステップであり、前記流体が前記繊維状物質から受ける力（ステップＳ６での計算値）と、前記流体中の前記繊維状物質の体積分率とを考慮して、特定の時刻から微小時間経過後における流体の速度を計算するステップに相当し得る。なお、このようなステップＳ６及びＳ７は、上記計算内容を考慮すれば、流体の運動と、繊維状物質の運動とを連成して解析するステップ（流体の運動と繊維状物質の運動との連成解析ステップ）であるといえ、これにより、繊維が流体から受ける力と流体が繊維状物質から受ける力（反力））を双方向（２ｗａｙ）から考慮した解析（連成解析）が可能となる。また、このような計算に際しては、前記流体中の前記繊維状物質の体積分率（ε）も考慮して解析を行うことから、本発明においては、様々な濃度の対象物（前記繊維状物質を含む流体）に対して、繊維状物質の運動状態（配向状態など）をシミュレーションできるものと本発明者らは考えている。

次に、ステップＳ８において、判定手段２３により、時刻Ｔが計算の終了時間として設定した時間（Ｔ_ｅｎｄ)よりも大きな値となっているか否かを判定する。そして、終了時間Ｔ_ｅｎｄよりも時刻Ｔが小さな値である場合には、時刻ＴをステップＳ２の時刻ｔ_ｎに変換して、新たにステップＳ２〜７を実行する。そして、ステップＳ８において、Ｔ_ｅｎｄ＜Ｔの条件を満たすまで、流体の運動と繊維状物質の運動を解析して求められる微小時間経過後の過渡的な解析結果（過渡的な流体及び繊維状物質の運動状態の情報）を順次利用して、新たにステップＳ２〜７を順次実行する処理を繰り返し行う。このようにしてステップＳ２〜７を順次繰り返し実行することによって、ステップＳ３〜７に記載の計算を所定条件を満たすまで繰り返し実行させること（反復演算させること）が可能となり、前記反復演算ステップ中に、流体の運動と繊維状物質の運動を解析して求められる微小時間経過後の過渡的な解析結果（過渡的な流体及び繊維状物質の運動状態の情報）を順次利用して、流体中の繊維状物質の運動を解析することができる。また、このような反復演算した解析は、流体の運動と繊維状物質の運動の相互作用（繊維が流体から受ける力と流体が繊維状物質から受ける力（反力））を双方向（２ｗａｙ）から考慮した解析（連成解析）であるため、これにより非常に精度の高いシミュレーションを行うことも可能となる。

そして、ステップＳ８において、時刻Ｔが終了時間（Ｔ_ｅｎｄ）よりも大きくなるという条件（Ｔ_ｅｎｄ＜Ｔという条件）を満たした場合には演算を終了し、ステップＳ９において演算結果を出力する。なお、このような演算により、初期の時刻から所定の時刻Ｔ_ｅｎｄまでの前記繊維状物質（球体集合体モデル）の各球体の速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角の情報とともに前記繊維状物質（球体集合体モデル）の各球体の中心の位置の情報を出力でき、これに基づいて、初期の時刻から所定の時刻Ｔ_ｅｎｄまでの間の流体中における前記繊維状物質の動き（移動位置）や前記繊維状物質の配向状態も出力できる。また、このような演算により、流体の速度の情報も併せて出力できる。このように、流体中の繊維状物質の運動状態を解析することで、繊維状物質及び流体の運動を解析でき、前記流体の速度、前記繊維状物質（球体集合体モデル）の各球体の速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角の情報の他、その解析情報に基づいて、前記繊維状物質の移動位置や配向状態をも演算して出力できる。なお、このような演算結果の出力は、外部記録媒体（出力部３の一態様）等にデータを記録することで行ってもよいし、ディスプレー（出力部３の一態様）等に解析結果を視覚的な情報（例えば、数値データ自体、数値データをグラフ化したもの、計算結果として求められる時間ごとの運動状態の情報に基づいて作成される３次元アニメーション等の画像、等といった視覚的な情報が挙げられる）として出力することで行ってもよい。

上述のように、フローチャート中のステップＳ２〜Ｓ５（特にＳ３〜Ｓ５）は、上記ステップ（ａ）として好適なステップであり、フローチャート中のステップＳ６〜Ｓ７は上記ステップ（ｂ）として好適なステップであり、フローチャート中のステップＳ８の判定結果に基づいてステップＳ２〜Ｓ７を繰り返し実行させるステップは上記ステップ（ｃ）として好適なステップである。そして、このようなステップ（ａ）〜（ｃ）を備える本発明の流体中の繊維状物質の運動状態の解析方法は、上述のように、コンピュータ（電子計算機）を利用した演算により、流体速度、繊維状物質の速度、位置情報、配向状態等を出力できる。このように、本発明の流体中の繊維状物質の運動状態の解析方法は、流体及び繊維状物質の運動状態（繊維状物質の移動位置や配向状態を含む。）を解析する方法であるともいえる。

なお、このようなステップ（ａ）〜（ｃ）を備える流体中の繊維状物質の運動状態の解析方法において、前記繊維状物質の形態としては、いわゆる繊維のような状態のものであればよく、その繊維長、繊維径、材料、形状は特に制限されず、繊維状（短繊維状や長繊維状等）、針状、棒状、錐状、管状等を含む概念である。また、このようなステップ（ａ）〜（ｃ）を備える流体中の繊維状物質の運動状態の解析方法において、そのシミュレーションの精度の観点から、前記繊維状物質としては、直径ｄが約１μｍ以上（ブラウン運動が支配的とならない大きさ）で、かつ、長さＬが計算領域を超えない程度である繊維状のものをシミュレーションの対象物とすることが好ましい。したがって、この範囲で定義されるアスペクト比（Ｌ／ｄ）であれば、原理的に計算可能である。このような繊維状物質としては、例えば、繊維強化樹脂等で用いられる炭素繊維やガラス繊維等が挙げられる。また、このような運動解析を行う流体としては、流動性を有する液状の物質であることが好ましく、例えば、一般的な樹脂や水等が挙げられる。

上述のような解析により、例えば、繊維強化樹脂（ＦＲＰ）などの成形プロセスの理解や、モデル化のために重要な対象である、図３に模式的に示すような単純せん断中の繊維状物質の配向挙動の３次元数値シミュレーションを実施することが可能となる。なお、ＦＲＰでは成形後の繊維配向状態が成形品の強度や剛性に対して重要なパラメーターであるため、このような数値シミュレーションにより、プロセス中の配向挙動を、より精度高く予測することが可能となる。このように、本発明によれば、流体中の繊維状物質の配向状態をより高精度にかつ現実的な計算時間で予測できるシミュレーションでき、各種成形プロセス最適化への応用が見込まれる。

以上、図１〜２等に基づいて本発明の流体中の繊維状物質の運動状態の解析方法及び解析装置の好適な実施形態について説明したが、本発明の流体中の繊維状物質の運動状態の解析方法及び本発明の流体中の繊維状物質の運動状態の解析装置は、上記実施形態に限定されるものではない。例えば、上記実施形態においてはステップＳ８において時間を条件として、第一の演算手段２１及び第二の演算手段２２での計算を順次実行する処理（ステップＳ３〜７に記載の計算を順次実行する処理：上記ステップ（ａ）及び（ｂ）に記載の計算を順次実行する処理）を繰り返し実行させているが、本発明の流体中の繊維状物質の運動状態の解析方法及び解析装置においては、かかる演算処理を反復実行させるための条件は時間に限定されるものではなく、例えば、Δｔの大きさを特定の大きさに設定しつつ繰り返し演算する回数を条件として設定して、上記ステップ（ａ）及び（ｂ）に記載の計算（第一の演算手段２１及び第二の演算手段２２での計算）が所定の繰り返し回数実行されるまで反復実行させてもよい。

なお、本発明の流体中の繊維状物質の運動状態の解析方法及び解析装置によれば、せん断力を加えられた流体中の繊維状物質の運動状態をより短時間でかつより精度高くシミュレーションすることが可能となる。このように、本発明の流体中の繊維状物質の運動状態の解析方法及び解析装置によって、流体中の繊維状物質の運動状態をより短時間でかつより精度高くシミュレーションすることが可能となる理由は必ずしも定かではないが、本発明者らは以下のように推察する。

ここで、先ず、本発明と対比するため、上述の特許文献１（特開平５−３１４０９１号公報）に記載のような、従来の流体中の繊維状物質の運動状態の解析方法について簡単に説明する。このような従来の流体中の繊維状物質の運動状態の解析方法においては、図７のフローチャートに示すようなステップＳ２１〜Ｓ２９の手順で解析が行われていた。このような解析について簡単に説明すると、先ず、ステップＳ２１及びＳ２２は、基本的に、上述の本発明の解析方法の好適な実施形態において説明したステップＳ１及びＳ２と同様のステップである。次いで、ステップＳ２３においては、下記式（１１１）〜（１１２）：

［式（１１１）〜（１１２）中の記号は、式（１２）及び（１３）中の記号と同様のものである。］
により、最初に、繊維状物質の存在や運動を考慮せずに流体の速度を計算する。このように、従来の流体中の繊維状物質の運動状態の解析方法においては、繊維状物質の存在を無視して、流体の速度を計算する。そして、そのような流体の速度の計算値を利用して、ステップＳ２４において繊維状物質の流体抵抗力及び繊維状物質の各球体が前記流体から受けるトルクを計算する。このようなステップＳ２４における計算は、流体の速度（ｕ^ｆ _ｉ）として、ステップＳ２３で計算した値を利用する以外は、上述の本発明の解析方法の好適な実施形態において説明したステップＳ３で説明した計算と同様の計算方法を採用する。そして、ステップＳ２５において繊維状物質の各球体間の結合力及び曲げ変形に対する復元トルクを計算する。このようなステップＳ２５における計算は、上述の本発明の解析方法の好適な実施形態において説明したステップＳ４で説明した計算と同様の計算方法を採用する。次に、ステップＳ２６において、相互に作用し合う繊維状物質を探索して、下記式（１１３）〜（１１５）：

［式（１１３）〜（１１５）中、ｆ^ｌｕｂ _ｉｊ及びｆ^ｍｅｃ _ｉｊは繊維状物質間の相互作用力を示し、ｄ_ｐは球体集合体モデル中の球体の直径を示し、ｍ^Ｐは球体の質量を示し、ｕ^Ｐ _ｉは球状集合体モデルの球体ｉの速度を示し、ｕ^P _ｊは球状集合体モデルの球体ｊの速度を示し、ｎ_ｉｊは球状集合体モデルの球体ｉ及び球体ｊのそれぞれの中心間を結ぶ球体ｉから球体ｊに向かう方向の単位ベクトルを示し、ｒ_ｉは球体ｉの重心の座標を示し、ｒ_ｊは球体ｊの重心の座標を示し、ｋ_ｓは引張り・圧縮変形に対する定数を示し、ｅ_ｄは球体間の反発係数を示し、ｆ^Ｐ _ｉｊは隣り合う２つの粒子間に生じる摩擦力を示す。］
により、繊維状物質間の相互作用力（２体相互作用）を計算する。なお、このようなｆ^ｌｕｂ _ｉｊ及びｆ^ｍｅｃ _ｉｊの概念は、Yamamoto, s. et al.，J.Chem Phys., 102, 2254(1995)やTsuji, Y. et al., Powder Tech., 77, 79(1993)に記載されているものを利用できる。次いで、ステップ２７において、このようなｆ^ｌｕｂ _ｉｊ及びｆ^ｍｅｃ _ｉｊの計算結果及びｆ^Ｐ _ｉｊの計算結果を利用して、下記式（７’）：

［式（７’）中、ｆ^ｌｕｂ _ｉｊ及びｆ^ｍｅｃ _ｉｊはステップＳ２６において算出した繊維状物質間の相互作用力を示し、それ以外の記号の意味は上記式（７）中のものと同様である。］
で表される計算式と上記式（８）〜（１０）で表される計算式とにより、各球体の速度（ｕ^Ｐ _ｉ）及び角速度（ω^Ｐ _ｉ）等を求めて、繊維状物質の位置を求める。次に、ステップＳ２８において、時刻Ｔが計算の終了時間として設定した時間（Ｔ_ｅｎｄ)よりも大きな値となっているか否かを判定し、終了時間Ｔ_ｅｎｄよりも時刻Ｔが小さな値である場合には、時刻ＴをステップＳ２の時刻ｔ_ｎに変換して、ステップＳ２２〜２７を再度実行し、Ｔ_ｅｎｄ＜Ｔの条件を満たした場合には演算を終了し、ステップＳ２９において演算結果を出力する。

このように、図７に示すような従来の解析方法では、初期条件を入力後（ステップＳ２１）、時間を更新し（ステップＳ２２）、繊維状物質の存在やその運動を考慮せずに流体の速度を計算し（ステップＳ２３）、次に、繊維状物質の各球体の流体抵抗力及び繊維状物質の各球体が前記流体から受けるトルクを計算（ステップＳ２４）し、維状物質の各球体間の結合力及び曲げ変形に対する復元トルクを計算する（ステップＳ２５）。次いで、繊維状物質間の２体相互作用を計算するために相互作用し合う相手を探索し、その結果に基づいて繊維状物質間の２体相互作用を計算する（ステップＳ２６）。なお、このようなステップＳ２６において２体相互作用として考慮する力（繊維状物質間の相互作用力）は、前述の流体潤滑力ｆ^ｌｕｂ _ｉｊと機械接触力ｆ^ｍｅｃ _ｉｊである。そして、このようなステップＳ２３〜Ｓ２７までの計算結果に基づいて、繊維の位置、速度及び角速度等を計算する（ステップＳ２７）。その後、指定の条件に到達していない場合は，ステップＳ２２に戻り同様の計算を繰り返す。

これに対して、図２に示す実施形態（本発明の解析方法の好適な一実施形態）では、初期条件を入力後（ステップＳ１）、時間を更新し（ステップＳ２）、繊維状物質の各球体の流体抵抗力及び繊維状物質の各球体が前記流体から受けるトルクを計算（ステップＳ３）し、繊維状物質の各球体間の結合力及び曲げ変形に対する復元トルクを計算（ステップＳ４）し、その計算結果に基づいて、繊維状物質の位置、速度及び角速度等を計算する（ステップＳ５）。次に、ステップＳ３の計算結果に基づいて流体が繊維状物質から受ける力（計算セル中に存在する各球体にかかる流体抵抗力ｆ^ｈ _ｉの反力を１つの計算セル内で総和した値）を計算し（ステップＳ６）、その計算結果に基づいて、繊維状物質から受ける反力を考慮して、流体の速度を計算する（ステップＳ７）。そして、指定の条件に到達していない場合は，ステップＳ２に戻り同様の計算を繰り返す。

このように、フローチャートに記載した実施形態の対比や上記計算式等から、本発明と上記特許文献１に記載のような従来技術との相違点は、従来技術において計算する２体相互作用の計算（ステップＳ２６）を本発明では省略している点、本発明では繊維状物質の位置、速度及び角速度等を計算する際に２体相互作用の計算結果（繊維状物質間の相互作用力：流体潤滑力ｆ^ｌｕｂ _ｉｊと機械接触力ｆ^ｍｅｃ _ｉｊの値）を利用しない点（本発明で利用する計算式（７）は、式（７’）中の流体潤滑力ｆ^ｌｕｂ _ｉｊと機械接触力ｆ^ｍｅｃ _ｉｊの値を０とした場合と同じ式である）、本発明では繊維状物質の存在やその運動を考慮して流体の速度を計算している点（ステップ６及び７：流体運動と繊維運動の相互作用（具体的には繊維が流体から受ける力およびその反力）を双方向（２ｗａｙ）に考慮する２ｗａｙ計算をしている点）にある。

ここで、従来技術において採用していた２体相互作用（粒子間相互作用：繊維状物質間の相互作用力）の計算について検討すると、流体潤滑や機械接触といった２体相互作用の計算では、繊維状物質のある粒子（粒子Ａとする）と別の粒子（粒子Ｂとする）の位置情報や、速度情報を互いに参照し合う必要があり、粒子Ａと粒子Ｂの情報が計算機のメモリ空間上において遠く離れた位置に格納されている場合、それらの通信時間に多くの時間を要することになる。特に、２体相互作用の計算において、並列計算を行う場合、このような通信時間はさらに長くなる傾向にある。このように、２体相互作用の計算は、上記従来の解析方法において最も計算負荷の大きい（計算時間のかかる）部分である。

また、このような従来の解析方法では、流体の速度の計算に繊維状物質の存在やその運動が考慮されていないため、せん断中の繊維状物質の運動において繊維状物質の配向が初期状態から定常値に到達するような基本的な配向挙動をシミュレーションが困難であった。これに対して、本発明においては、このような２体相互作用（繊維状物質間に生じる流体潤滑や機械接触）の計算、すなわち、繊維状物質間の相互作用力（流体潤滑力ｆ^ｌｕｂ _ｉｊ及び機械接触力ｆ^ｍｅｃ _ｉｊ）の計算を省略し、その代替として、流体運動と繊維運動の相互作用（具体的には繊維が流体から受ける力およびその反力）を双方向（２ｗａｙ）に考慮する２ｗａｙ計算をしている。そのため、本発明においては、計算精度を低下させることなく、計算速度が大幅に向上すること（計算速度を大幅に短くすること）が可能となるものと本発明者らは推察する。

また、従来手法では、繊維が流体から受ける力のみを考慮する１方向（１ｗａｙ）の計算手法により、流体中の繊維状物質の運動状態を解析しているが、本発明では、上述のように、繊維状物質が流体から力を受けることによって生じる運動だけではなく、繊維の存在やその運動が流体の運動にも影響を与えることに着目して、繊維状物質と流体の相互作用を双方向（２ｗａｙ）に考慮して計算をしているため、せん断中の繊維状物質の運動において繊維状物質の配向が初期状態から定常値に到達するような基本的な配向挙動についてのシミュレーションを行うことも可能である。なお、このようなシミュレーションは、繊維強化樹脂（ＦＲＰ）などの成形プロセスの理解やモデル化のために好適に利用することが可能である。また、このようなシミュレーションに際して、上述のように、流体と繊維状物質との間の相互作用を２ｗａｙで考慮するため、繊維状物質の配向の予測精度もより向上したものとなる。このように、本発明によれば、流体の運動と繊維状物質の運動の相互作用（具体的には、繊維が流体から受ける力およびその反力）を双方向（２ｗａｙ）に考慮するため、従来と比較して、せん断中の繊維状物質の配向状態の予測をより精度高く行うことが可能となるものと本発明者らは推察している。

以上、説明したように、本発明では、繊維状物質間に生じる流体潤滑や機械接触といった２体相互作用は計算せずに、流体及び繊維状物質の基礎方程式において、流体の運動と繊維状物質の運動の相互作用を双方向に考慮する２ｗａｙ計算を行うため、従来技術と比較してせん断中の繊維配向をより高精度に予測することが可能となる。また、このような本発明の解析方法では、従来の解析方法において最も計算負荷の大きかった２体相互作用の計算を省略して、上述の双方向の相互作用の考慮により（流体の運動と繊維状物質の運動の相互作用の計算により）代替して流体中の繊維状物質の運動状態を解析するため、計算速度がより向上するものと本発明者らは推察する。

以下、実施例及び比較例に基づいて本発明をより具体的に説明するが、本発明は以下の実施例に限定されるものではない。

（実施例１）
実施例１においては、ニュートン流体からなる流体中に下記式：

［式中、ｄ_ｐは繊維状物質の直径（球体集合体モデルとした場合の球体の直径）を示し、γはせん断速度を示し、ρ^ｐは繊維状物質の密度を示し、μは流体の粘性係数を示し、Ｓｔはストークス数を示す］
で表されるストークス数（Ｓｔ）がほぼ１となるような繊維状物質の密度（ρ^ｐ）と直径（ｄ_ｐ）を有する繊維が分散した系について、図２に示すフローチャートに基づいて、上記式（１）〜（１３）に記載の計算を行って、せん断力をかけ始めてから所定の時間までの流体の運動と繊維状物質の運動を解析して、流体中の繊維状物質の配向状態を解析した。

なお、このようにして解析（シミュレーション）する対象の実験は、文献“Cieslinski, M.J. et al., J. Non-Newtonian Fluid Mech., vol.222, 163頁〜170頁(2015年発行)”「以下、かかる文献を場合により単に「参考文献１」と称する）に記載されているスライディング型レオメータを用いた繊維状物質を体積濃度３．８ｖｏｌ％の割合（Ｖｆ＝３．８％）で含む流体に対する実験を模したものであり、上記解析は、具体的には、該参考文献１に記載のスライディング型レオメータ（具体的には大変形、かつ、一様変形を発生させ、時系列データを取得できるスライディング型レオメータ）を用いて、ポリプロピレンからなる流体にガラス繊維を体積濃度（Ｖｆ）２．２ｖｏｌ％の割合で分散させた繊維含有流体（流体（ポリプロピレン）及び繊維状物質（ガラス繊維））に対してせん断力をかけた場合の実験についてのシミュレーションとして行った。

なお、このような解析に際しては、初期条件として、繊維状物質の密度（ρ^ｐ）、繊維状物質の直径（ｄ_ｐ）、繊維状物質の直径に対する繊維長Ｌの比（Ｌ／ｄ_ｐ）であるアスペクト比（Ｇ＝２０）、繊維の本数（６９０本）、繊維状物質を体積濃度（Ｖｆ＝２．２ｖｏｌ％）、流体の粘性係数（μ）、せん断速度（γ）といった情報を入力した。なお、上記参考文献１に記載の実験と同様に、粒子に働く慣性力が支配的とならない条件で、かつ数値的に安定に計算できる条件として、ストークス数がほぼ１（概ね１）以下となるような、ｄ_ｐ、ρ^ｐ、γ、μの組み合わせとして、ｄ_ｐ＝２５［μｍ］、ρ^ｐ＝１×１０^６［ｋｇ／ｍ^３］、γ＝１０００［１／秒］、μ＝０．１［Ｐａ・ｓ］という値を入力した。また、計算上の初期の繊維配向が、Ａｘｘ〜０．５，Ａｙｙ〜０，Ａｚｚ〜０．５に設定されるように、繊維の初期条件を入力した。ここで、Ａｘｘはｘ方向（流れ方向）の配向テンソルを示し、Ａｙｙはｙ方向（壁面垂直方向）の配向テンソルを示し、Ａｚｚはｚ方向（流路の奥行方向）の配向テンソルを示す。

（実施例２）
繊維の本数を６９０本から４１４０本に変更して、初期条件の繊維状物質の体積濃度の値が２．２ｖｏｌ％（Ｖｆ＝２．２％）から１２．７ｖｏｌ％（Ｖｆ＝１２．７％）に変わった以外は、実施例１と同様にして、せん断応力をかけ始めてから所定の時間までの流体の運動と繊維状物質の運動を解析して、流体中の繊維状物質の配向状態を解析した。

なお、このようにして解析（シミュレーション）する対象の実験は、前記参考文献１［“Cieslinski, M.J. et al., J. Non-Newtonian Fluid Mech., vol.222, 163頁〜170頁(2015年発行)］に記載されているスライディング型レオメータを用いた繊維状物質を体積濃度１３．３ｖｏｌ％の割合（Ｖｆ＝１３．３％）で含む流体に対する実験を模したものであり、上記解析は、具体的には、該参考文献１に記載のスライディング型レオメータを用いて、ポリプロピレンからなる流体にガラス繊維を体積濃度（Ｖｆ）１２．７ｖｏｌ％の割合で分散させた繊維含有流体（流体（ポリプロピレン）及び繊維状物質（ガラス繊維））に対してせん断力をかけた場合の実験についてのシミュレーションとして行った。

（比較例１）
図２に示すフローチャートに基づいて上記式（１）〜（１３）に記載の計算を行う代わりに、図７に示すフローチャートに基づいて式（１）〜（６）、式（７’）、式（８）〜（１０）及び式（１１１）〜（１１５）に記載の計算を行う以外は、実施例１と同様にして、せん断応力をかけ始めてから所定の時間までの流体の運動と繊維状物質の運動を解析して、流体中の繊維状物質の配向状態を解析した。

（比較例２）
図２に示すフローチャートに基づいて上記式（１）〜（１３）に記載の計算を行う代わりに、図７に示すフローチャートに基づいて式（１）〜（６）、式（７’）、式（８）〜（１０）及び式（１１１）〜（１１５）に記載の計算を行う以外は、実施例２と同様にして、せん断応力をかけ始めてから所定の時間までの流体の運動と繊維状物質の運動を解析して、流体中の繊維状物質の配向状態を解析した。

〈配向状態の解析結果について〉
実施例１及び比較例１の解析結果として、せん断応力をかけ始めてからの繊維配向テンソルＡｘｘの時間変化（時間をせん断速度で無次元化した値，すなわち，ひずみ）に関するグラフを図８に示す。なお、図８には、繊維状物質の体積分縦率の値が近いことから、参考文献１［“Cieslinski, M.J. et al., J. Non-Newtonian Fluid Mech., vol.222, 163頁〜170頁(2015年発行)”］に記載されているスライディング型レオメータを用いた繊維状物質を体積分率（ε）３．８ｖｏｌ％の割合で含む流体に対する実験の結果も参考例１として併せて示す。

なお、このように、参考のために図８に記載する前記参考文献１に記載されているスライディング型レオメータを用いた実験の結果は、該参考文献１に記載のスライディング型レオメータを用い、直径ｄの平均値が１３．３μｍであり、かつ、アスペクト比Ｇが１０５であるガラス繊維（繊維状物質）を体積濃度（Ｖｆ）の条件が３．８ｖｏｌ％となる割合でポリプロピレン（流体）に分散させ、せん断速度（γ）が１．０［１／秒］となる条件でせん断を加える実験における測定結果である。

また、実施例２及び比較例２の解析結果として、せん断応力をかけ始めてからの繊維配向テンソルＡｘｘの時間変化（時間をせん断速度で無次元化した値，すなわち，ひずみ）に関するグラフを図９に示す。なお、図９には、繊維状物質の体積分立の値が近いことから、前記参考文献１に記載されているスライディング型レオメータを用いた繊維状物質を体積分率（ε）：１３．３ｖｏｌ％の割合で含む流体に対する実験の結果も参考例２として併せて示す。

なお、このように参考のために図９に記載する前記参考文献１に記載されているスライディング型レオメータを用いた実験の結果は、該参考文献１に記載のスライディング型レオメータを用い、直径ｄの平均値が１３．３μｍであり、かつ、アスペクト比Ｇが１０５であるガラス繊維（繊維状物質）を体積濃度（Ｖｆ）の条件が１３．３ｖｏｌ％となる割合でポリプロピレン（流体）に分散させ、せん断速度（γ）が１．０［１／秒］となる条件でせん断を加える実験における測定結果である。

図８及び図９に示す結果からも明らかなように、比較例１及び２で採用している従来の解析方法（以下、場合により「Int1way」と称する）により得られたシミュレーション結果［体積濃度（Ｖｆ）が２．２％の場合のシミュレーション結果（比較例１：図８）及びＶｆが１２．７％の場合（比較例２：図９）のシミュレーション結果］において、前記参考文献１に記載されている実験結果とは全く異なる挙動を示していることが明らかとなった。これに対して、実施例１及び２で採用している本発明の流体中の繊維状物質の運動状態の解析方法（以下、場合により「Noint2way」と称する）により得られたシミュレーション結果［体積分率（濃度）がＶｆが２．２％の場合のシミュレーション結果（実施例１：図８）及びＶｆが１２．７％のシミュレーション場合（実施例２：図９）の結果］は、初期値から徐々に定常値に達する挙動を示しており、上記参考文献１に記載されている実験結果とグラフの傾向がほぼ一致しており、配向状態の解析結果と実験による配向状態の観察結果がほぼ同様な結果となることが分かった。このような結果から、本発明で採用する解析方法（Noint2way）によれば、従来の解析方法（Int1way）と比較して、せん断応力を印加した場合の流体中の繊維状物質の運動状態（本実施例では配向挙動）を、より精度高くシミュレーションできることが確認された。

〈解析時間について〉
実施例１〜２及び比較例１〜２で行った解析に関して、６４並列（４ＣＰＵ、１６コア／ＣＰＵ）の環境下で、繰り返し演算する回数が５０００となるような計算（時刻：ｔｎのｎが５０００となるまでの計算：５０００ステップの計算）を実施した場合の計算時間を測定した。このような計算結果として、体積濃度Ｖｆが２．２％の場合、体積濃度Ｖｆが１２．７％の場合のそれぞれのケースにおけるＮｏｉｎｔ２ｗａｙ（実施例１又は実施例２）で採用している解析方法の計算時間を１とした場合の計算時間比を図１０に示す。

図１０に示す結果からも明らかなように、Ｖｆが２．２％の場合では、Noint2wayによる計算方法（実施例１で採用している解析方法）は、Int1wayによる計算方法（比較例１で採用している解析方法）に比べて、計算時間が約４倍も速くなることが確認された。また、Ｖｆが１２．７％の場合では、Noint2wayによる計算方法（実施例２で採用している解析方法）は、Int1wayによる計算方法（比較例２で採用している解析方法）に比べて、計算時間が約１２倍も速くなることが確認された。このような結果から、Noint2wayによる計算方法（本発明で採用する解析方法）によれば、より高精度のシミュレーション結果をより短時間で得ることが可能となることが分かった。

なお、体積濃度Ｖｆをさらに増加させた場合（繊維の本数をさらに増加させた場合）、２体相互作用を計算する際に生じるメモリ空間上の離れた位置への通信が増大するため、前記Noint2wayによる計算方法（本発明で採用する解析方法）では、Int1wayによる計算方法（従来の解析方法）と比べて、高速化の効果が更に向上するものと考えられる。

〈考察〉
前記Noint2wayによる計算方法（本発明で採用する解析方法）の適用範囲について検討すると、一般に、粒子含有流体の計算において相互作用に関する議論を行う場合、粒子濃度が重要なパラメーターのひとつとなる。ここで、粒子濃度は、一般に、濃度が低い順で並べると、希薄系、準濃厚系、濃厚系の３つの系に分類されている。そして、実施例において解析したような繊維を含有する流体の場合には、上述のような３つの系は、下記条件式：

で表されるような条件で分類できることが知られている（Doi, M. and Edwards, S.F., The Theory of Polymer Dynamics, Oxford Univ. Press, New York(1988)）。なお、このような条件式においてｎは単位体積当たりの繊維の本数である。そして、希薄系では、本発明において省略した機械接触等の２体相互作用はそれほど重要ではないが（無視できる程度であるが）、濃厚系では２体相互作用は支配的なパラメータとなる可能性がある。そこで、実施例においてシミュレーションしたアスペクト比Ｌ／ｄ＝２０の繊維について、上記条件式にあてはめて分類すると、濃厚系となる条件は、体積濃度Ｖｆが３．９％となる場合であることが分かる。一方で、上述の実施例２においては、図９に示す結果からも明らかなように、体積濃度Ｖｆが１２．７％の場合においても精度高くシミュレーションできることが確認されている。このような実施例１〜２の結果と、上記一般的な分類とを併せ勘案すれば、Noint2wayによる計算方法（本発明で採用する解析方法）は、希薄系〜濃厚系のいずれの濃度の分類の範囲（濃度分類の全範囲）において、十分に精度高く繊維状物質の運動状態（配向状態など）を解析できること明白であり、どのような濃度の系においても適用できるものといえる。

以上説明したように、本発明によれば、流体中の繊維状物質の運動状態をより短時間でかつより精度高くシミュレーションすることを可能とする流体中の繊維状物質の運動状態の解析方法及びそのような解析方法を実行することが可能な流体中の繊維状物質の運動状態の解析装置を提供することが可能となる。

したがって、本発明の流体中の繊維状物質の運動状態の解析方法は、繊維強化樹脂（ＦＲＰ）や繊維強化金属（ＦＲＭ）等の複合材料の成形条件や成形用金型等の設計指針を導き出すためのシミュレーションの方法等として特に有用である。

１…入力部、２…演算処理部、２１…第一の演算手段、２２…第二の演算手段、２３…判定手段、３…出力部、１０…繊維状物質、１０Ａ……繊維状物質の球体集合体モデル。

Claims (4)

1. 流体の運動方程式と粒子の運動方程式を利用して繊維状物質及び流体の運動を計算する上で必要となる繊維状物質及び流体の情報に基づいて、流体中の繊維状物質の運動状態をシミュレーションするための解析方法であって、下記ステップ（ａ）〜（ｃ）：
   ステップ（ａ）：前記繊維状物質が複数の球体が結合した球体集合体からなるものと仮定する球体集合体モデルを用いる計算方法を用いて、特定の時刻から微小時間経過後における前記繊維状物質の運動を解析するためのステップであり、
   流体の速度の情報に基づいて、前記繊維状物質の各球体の流体抵抗力及び前記繊維状物質の各球体が前記流体から受けるトルクを計算し、
   前記繊維状物質の各球体間の結合力及び前記繊維状物質の各球体間の曲げ変形に対する復元トルクを計算し、
   前記繊維状物質の各球体の流体抵抗力、前記繊維状物質の各球体が前記流体から受けるトルク、前記繊維状物質の各球体間の結合力、及び、前記繊維状物質の各球体間の曲げ変形に対する復元トルクの情報に基づいて、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算する、繊維状物質の運動解析ステップ；
   ステップ（ｂ）：前記ステップ（ａ）で計算した前記繊維状物質の各球体の流体抵抗力の情報に基づいて、前記流体が前記繊維状物質から受ける力を計算し、前記流体が前記繊維状物質から受ける力と前記流体中の前記繊維状物質の体積分率とを考慮して、特定の時刻から微小時間経過後における流体の速度を計算する、流体の運動と繊維状物質の運動との連成解析ステップ；
   ステップ（ｃ）：前記ステップ（ａ）及びステップ（ｂ）で既に求められた過渡的な流体及び繊維状物質の運動状態の情報を利用して、新たに前記ステップ（ａ）〜（ｂ）の計算を順次実行する処理を、所定条件を満たすまで繰り返し実行する反復演算ステップ；
   を含み、前記反復演算ステップ中に流体の運動と繊維状物質の運動を解析して求められる微小時間経過後の過渡的な解析結果を順次利用して、流体の運動と繊維状物質の運動を連成解析し、流体中の繊維状物質の運動状態を解析することを特徴とする流体中の繊維状物質の運動状態の解析方法。
2. 前記繊維状物質の運動解析ステップにおいて、繊維状物質間の相互作用力に関する流体潤滑力と機械接触力を計算することなく、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算することを特徴とする請求項１に記載の流体中の繊維状物質の運動状態の解析方法。
3. 流体の運動方程式と粒子の運動方程式を利用して繊維状物質及び流体の運動を計算する上で必要となる繊維状物質及び流体の情報に基づいて、流体中の繊維状物質の運動状態をシミュレーションするための解析装置であって、
   前記繊維状物質が複数の球体が結合した球体集合体からなるものと仮定する球体集合体モデルを用いる計算方法を用いて、特定の時刻から微小時間経過後における前記繊維状物質の運動を解析するための演算手段であり、流体の速度の情報に基づいて、前記繊維状物質の各球体の流体抵抗力及び前記繊維状物質の各球体が前記流体から受けるトルクを計算し、前記繊維状物質の各球体間の結合力及び前記繊維状物質の各球体間の曲げ変形に対する復元トルクを計算し、前記繊維状物質の各球体の流体抵抗力、前記繊維状物質の各球体が前記流体から受けるトルク、前記繊維状物質の各球体間の結合力、及び、前記繊維状物質の各球体間の曲げ変形に対する復元トルクの情報に基づいて、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算する第一の演算手段と；
   前記第一の演算手段により計算された前記繊維状物質の各球体の流体抵抗力の情報に基づいて、前記流体が前記繊維状物質から受ける力を計算し、前記流体が前記繊維状物質から受ける力と前記流体中の前記繊維状物質の体積分率とを考慮して、特定の時刻から微小時間経過後における流体の速度を計算する、第二の演算手段と；
   前記第一の演算手段及び前記第二の演算手段で既に求められた過渡的な流体及び繊維状物質の運動状態の情報を利用して、新たに前記第一の演算手段及び前記第二の演算手段での計算を順次実行する処理を、所定条件を満たすまで繰り返す反復演算ステップを実行させるために、該条件を満たしたか否かを判定する、判定手段と；
   を備え、前記反復演算ステップ中に流体の運動と繊維状物質の運動を解析して求められる微小時間経過後の過渡的な解析結果を順次利用して、流体の運動と繊維状物質の運動を連成解析し、流体中の繊維状物質の運動状態を解析することを特徴とする流体中の繊維状物質の運動状態の解析装置。
4. 前記第一の演算手段が、繊維状物質間の相互作用力に関する流体潤滑力と機械接触力を計算することなく、特定の時刻から微小時間経過後における前記繊維状物質の各球体の位置、速度及び角速度並びに前記繊維状物質を構成する隣り合う球体間の曲げ角を計算する演算手段であることを特徴とする請求項３に記載の流体中の繊維状物質の運動状態の解析装置。