(57)【要約】【課題】繊維状物体の集合の運動を一つ一つ計算する場合と比較して、計算時間を低減させつつ、同等の計算結果が得られるようにする。【解決手段】作用力算出部３６が、繊維状物体の集合の一部の繊維状物体をモデル化した繊維モデルを用いて、予め設定され、又は計算される流体に関する条件と、繊維モデルの各々の位置とに基づいて、毎時刻、繊維モデル間の各々に相互に作用する力、及び繊維モデルの単体に作用する力を計算する。このとき、繊維モデルの繊維半径を所定倍に拡大して、相互に作用する力を計算する。挙動計算部４４が、毎時刻、繊維モデル毎に、他の繊維モデルとの間で相互に作用する力、及び繊維モデルの単体に作用する力に基づいて、繊維モデルの挙動を計算する。【選択図】図４
2020064450.tif 000002
【特許請求の範囲】
【請求項１】
流体中に含まれる繊維状物体の集合の挙動を計算する繊維挙動計算装置であって、
前記繊維状物体の集合の一部の前記繊維状物体をモデル化した繊維モデルを用いて、予め設定され、又は計算される流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル間の各々に相互に作用する力を計算する相互作用力算出部と、
前記流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル毎に、前記繊維モデルの単体に作用する力を計算する単体作用力算出部と、
毎時刻、前記繊維モデル毎に、他の繊維モデルとの間で相互に作用する力、及び前記繊維モデルの単体に作用する力に基づいて、前記繊維モデルの挙動を計算する挙動計算部と、を含み、
前記相互作用力算出部は、前記繊維モデルの繊維半径を所定倍に拡大して、前記相互に作用する力を計算する
繊維挙動計算装置。
【請求項２】
前記所定倍は、前記流体中の、繊維状物体の集合の体積含有率に基づいて定められる請求項１記載の繊維挙動計算装置。
【請求項３】
前記繊維モデルは、前記繊維状物体のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルであり、
各球の位置に基づいて、所定時間間隔で、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部を更に含み、
前記相互作用力算出部は、毎時刻、球毎に、前記球の位置に基づいて、他の繊維モデルとの間で相互に作用する力を算出し、
前記単体作用力算出部は、毎時刻、球毎に、前記繊維モデルの単体に作用する力を計算し、
前記挙動計算部は、
毎時刻、球毎に、前記補正係数を乗じた前記他の繊維モデルとの間で相互に作用する力、及び前記補正係数を乗じた前記繊維モデルの単体に作用する力に基づいて、前記球の挙動を計算する請求項１又は２記載の繊維挙動計算装置。
【請求項４】
前記粗視化繊維モデルは、前記ｎ個の球と、球と球との間を結ぶセグメントとを含み、
前記作用する力は、前記セグメントの両端の球に作用する力と、前記セグメントに渡って分布して作用する力を含み、
前記補正係数算出部は、各球の位置に基づいて、所定時間間隔で、球毎に、前記セグメントの両端の球に作用する力、及び前記セグメントに渡って分布して作用する力の各々に対する補正係数を計算し、
前記挙動計算部は、毎時刻、球毎に、前記補正係数を乗じた前記セグメントの両端の球に作用する力と、前記補正係数を乗じた前記セグメントに渡って分布して作用する力との和に基づいて、前記球の挙動を計算する請求項３記載の繊維挙動計算装置。
【請求項５】
流体中に含まれる繊維状物体の集合の挙動を計算する繊維挙動計算装置における繊維挙動計算方法であって、
相互作用力算出部が、前記繊維状物体の集合の一部の前記繊維状物体をモデル化した繊維モデルを用いて、予め設定され、又は計算される流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル間の各々に相互に作用する力を計算し、
単体作用力算出部が、前記流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル毎に、前記繊維モデルの単体に作用する力を計算し、
挙動計算部が、毎時刻、前記繊維モデル毎に、他の繊維モデルとの間で相互に作用する力、及び前記繊維モデルの単体に作用する力に基づいて、前記繊維モデルの挙動を計算することを含み、
前記相互作用力算出部が計算することでは、前記繊維モデルの繊維半径を所定倍に拡大して、前記相互に作用する力を計算する
繊維挙動計算方法。
【請求項６】
流体中に含まれる繊維状物体の集合の挙動を計算するためのプログラムであって、
コンピュータを、
前記繊維状物体の集合の一部の前記繊維状物体をモデル化した繊維モデルを用いて、予め設定され、又は計算される流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル間の各々に相互に作用する力を計算する相互作用力算出部、
前記流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル毎に、前記繊維モデルの単体に作用する力を計算する単体作用力算出部、及び
毎時刻、前記繊維モデル毎に、他の繊維モデルとの間で相互に作用する力、及び前記繊維モデルの単体に作用する力に基づいて、前記繊維モデルの挙動を計算する挙動計算部として機能させるためのプログラムであって、
前記相互作用力算出部は、前記繊維モデルの繊維半径を所定倍に拡大して、前記相互に作用する力を計算する
プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、繊維挙動計算装置、方法、及びプログラムに係り、特に、流体中の繊維の挙動を計算する繊維挙動計算装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
従来より、流体中を流れる繊維一本一本の運動をシミュレートする技術が知られている（特許文献１、２）。
【０００３】
また、流体中における粉粒体の運動をシミュレートする技術が知られている（特許文献３）。特許文献３では、実粒子よりも大きな径の粒子（粗視化粒子）を用いて、実粒子群の運動を模擬させることで、解析に必要となる粒子数を低減している。このとき、実粒子群と粗視化粒子の運動が一致するように、粒子径の拡大率に応じて粗視化粒子に作用する力を変更する。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１８−８８０８５号公報
【特許文献２】特開２０００−３２２４０７号公報
【特許文献３】特開２０１０−１０７３１４号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
上記特許文献１、２に記載の技術では、例えば、繊維強化樹脂の射出成形過程における繊維運動を解析する場合、膨大な数の繊維を取り扱う必要がある。さらに、繊維同士の衝突を考慮する場合には、さらに計算量が増大する。これは、衝突力を算出するにあたり、周囲に存在する多数の繊維との距離を算出する必要があるためである。
【０００６】
また、上記特許文献３は、粉粒体を対象としたものであり、繊維の運動を取り扱うものではない。
【０００７】
本発明は、上記の事情を鑑みてなされたもので、繊維状物体の集合の運動を一つ一つ計算する場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる繊維挙動計算装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の目的を達成するために本発明に係る繊維挙動計算装置は、流体中に含まれる繊維状物体の集合の挙動を計算する繊維挙動計算装置であって、前記繊維状物体の集合の一部の前記繊維状物体をモデル化した繊維モデルを用いて、予め設定され、又は計算される流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル間の各々に相互に作用する力を計算する相互作用力算出部と、前記流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル毎に、前記繊維モデルの単体に作用する力を計算する単体作用力算出部と、毎時刻、前記繊維モデル毎に、他の繊維モデルとの間で相互に作用する力、及び前記繊維モデルの単体に作用する力に基づいて、前記繊維モデルの挙動を計算する挙動計算部と、を含み、前記相互作用力算出部は、前記繊維モデルの繊維半径を所定倍に拡大して、前記相互に作用する力を計算する。
【０００９】
本発明に係る繊維挙動計算方法は、流体中に含まれる繊維状物体の集合の挙動を計算する繊維挙動計算装置における繊維挙動計算方法であって、相互作用力算出部が、前記繊維状物体の集合の一部の前記繊維状物体をモデル化した繊維モデルを用いて、予め設定され、又は計算される流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル間の各々に相互に作用する力を計算し、単体作用力算出部が、前記流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル毎に、前記繊維モデルの単体に作用する力を計算し、挙動計算部が、毎時刻、前記繊維モデル毎に、他の繊維モデルとの間で相互に作用する力、及び前記繊維モデルの単体に作用する力に基づいて、前記繊維モデルの挙動を計算することを含み、前記相互作用力算出部が計算することでは、前記繊維モデルの繊維半径を所定倍に拡大して、前記相互に作用する力を計算する。
【００１０】
本発明に係るプログラムは、流体中に含まれる繊維状物体の集合の挙動を計算するためのプログラムであって、コンピュータを、前記繊維状物体の集合の一部の前記繊維状物体をモデル化した繊維モデルを用いて、予め設定され、又は計算される流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル間の各々に相互に作用する力を計算する相互作用力算出部、前記流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル毎に、前記繊維モデルの単体に作用する力を計算する単体作用力算出部、及び毎時刻、前記繊維モデル毎に、他の繊維モデルとの間で相互に作用する力、及び前記繊維モデルの単体に作用する力に基づいて、前記繊維モデルの挙動を計算する挙動計算部として機能させるためのプログラムであって、前記相互作用力算出部は、前記繊維モデルの繊維半径を所定倍に拡大して、前記相互に作用する力を計算する。
【００１１】
本発明によれば、相互作用力算出部が、前記繊維状物体の集合の一部の前記繊維状物体をモデル化した繊維モデルを用いて、予め設定され、又は計算される流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル間の各々に相互に作用する力を計算する。このとき、相互作用力算出部は、繊維モデルの繊維半径を所定倍に拡大して、相互に作用する力を計算する。
【００１２】
また、単体作用力算出部が、前記流体に関する条件と、前記繊維モデルの各々の位置とに基づいて、毎時刻、前記繊維モデル毎に、前記繊維モデルの単体に作用する力を計算する。そして、挙動計算部が、毎時刻、前記繊維モデル毎に、他の繊維モデルとの間で相互に作用する力、及び前記繊維モデルの単体に作用する力に基づいて、前記繊維モデルの挙動を計算する。
【００１３】
このように、繊維モデルの繊維半径を所定倍に拡大して、他の繊維モデルとの間で相互に作用する力を計算し、他の繊維モデルとの間で相互に作用する力、及び前記繊維モデルの単体に作用する力に基づいて、前記繊維状物体の集合の一部の前記繊維状物体をモデル化した繊維モデルの各々の挙動を計算することにより、繊維状物体の集合の運動を一つ一つ計算する場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる。
【発明の効果】
【００１４】
以上説明したように、本発明の繊維挙動計算装置、方法、及びプログラムによれば、繊維状物体の集合の運動を一つ一つ計算する場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる、という効果が得られる。
【図面の簡単な説明】
【００１５】
【図１】繊維及び粗視化繊維モデルの一例を示す図である。
【図２】剛体反発力を説明するための図である。
【図３】本発明の実施の形態に係る繊維挙動計算装置を示すブロック図である。
【図４】本発明の実施の形態に係る繊維挙動計算装置を示す機能ブロック図である。
【図５】本発明の実施の形態に係る繊維挙動計算装置の繊維挙動計算処理ルーチンの内容を示すフローチャートである。
【図６】解析対象の形状を示す図である。
【図７】短繊維における繊維配向テンソルのｚ方向分布を示すグラフである。
【図８】計算高速化割合を示すグラフである。
【図９】長繊維における繊維配向テンソルのｚ方向分布を示すグラフである。
【図１０】長繊維の運動解析結果である繊維の様子を示す図である。
【発明を実施するための形態】
【００１６】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１７】
＜本発明の実施の形態の概要＞
本発明の実施の形態では、流体中を流れる繊維状物体（繊維半径a、流体中の繊維体積含有率Vf）の運動シミュレーションを行う際に、実際よりも少ない繊維のみの運動を解く。この繊維同士の相互作用力を計算する場合には、繊維半径aがer倍（er>1）に拡大したものとして算出することで、繊維体積含有率の違いによる相互作用力の違いを補完する。一方で、流体からの力、繊維の引張復元力、曲げ復元力などの、相互作用力以外の力については繊維半径をaのままとして算出し、上記で算出した力を用いて繊維の運動方程式を解き、繊維の位置、変形及び配向といった挙動を解析する。
【００１８】
これにより、計算に必要となる繊維の本数が従来手法と比べて低減されるため、計算量が抑えられる。また、繊維同士の相互作用力計算を行う際に、代表繊維の半径を拡大させることで、従来手法を用いた場合と同じ繊維体積含有率を得ることができる。また、引張復元力や曲げ復元力などを算出する場合に繊維半径を拡大させた場合には、繊維の剛性が大きくなり、従来手法とは異なる変形挙動を示す。このため、繊維半径を拡大せずにこれらの力を算出することで、従来手法と同等の変形挙動が得られる。
【００１９】
＜本発明の実施の形態の原理＞
繊維を球の連結体としてモデル化する。このとき図１のように、連結した球間の距離が球の直径よりも大きくなるようにモデル化する。ここで、この繊維の内、結合した2つの球からなる要素をセグメントと呼称する。例として、図１の繊維のセグメント数は3である。
【００２０】
このような繊維モデルを用いて、繊維体積含有率がVfである系を解くことを考える。本実施の形態では、この繊維体積含有率Vfを得るために必要となる繊維本数よりも少ない繊維本数を用いて解析を行う。このままでは繊維体積含有率が低下するため、元の系の繊維体積含有率Vfと同等になるよう、繊維の半径aをer倍に拡大する。erは半径拡大率である。ここで、繊維半径を単に拡大した場合、繊維の剛性が大きくなり、曲げ変形などの挙動が異なる。このため、繊維同士の相互作用力を算出する場合のみ、繊維半径が拡大したものとして取り扱う。
【００２１】
具体的には、繊維モデルを構成する各球iについて次の運動方程式を解く。
【００２２】
【数１】
2020064450.tif 000003
【００２３】
ここで、jは球iと結合している球の番号、nsはセグメントの番号、iwは壁面の番号、mは球の質量、vは球の速度、Fhは粘性力、FSは引張復元力、Fbは曲げ復元力、Fhtは粘性トルクに起因する力、FRepは繊維間の衝突時に働く剛体反発力、Fwは壁面反発力である。なお、粘性力、引張復元力、曲げ復元力、粘性トルクに起因する力、及び壁面反発力は、繊維モデルの単体に作用する力の一例である。また、繊維間の衝突時に働く剛体反発力は、繊維モデル間に相互に作用する力の一例である。
【００２４】
剛体反発力FRepの算出方法を以下に示す。まず、近接する2つのセグメントisおよびnsを考える（図２）。セグメントisは球iおよびjから構成されているとする。セグメントisがnsより受ける剛体反発力FRsegを次式にて算出する。
【００２５】
【数２】
2020064450.tif 000004
【００２６】
ただし、xcis、 xcnsはそれぞれセグメントis、nsの最近接点の座標、D0、G0は定数である。最後に、セグメントisを構成する球i、jに働く剛体反発力FRepを次式にて算出する。
【００２７】
【数３】
2020064450.tif 000005
(1)
【００２８】
ただし、
2020064450.tif 000006
であり、xi、xjは球iおよびjの座標である。
【００２９】
剛体反発力FRep以外の各力は、以下のように算出される。
【００３０】
まず、壁面反発力Fwは例えば壁面への球の貫入量に応じてペナルティ法等により算出する。また、球iに働く粘性力Fh、球i、jに働く引張復元力FSは、次式にて算出される。
【００３１】
【数４】
2020064450.tif 000007
【００３２】
ただし、aは球の半径、Eは繊維の引張弾性率、η(ri)およびV(xi)は球の位置xiにおける流体の粘度および速度ベクトル、rijは球iから球jへ向かう方向ベクトル、nijは球iから球jへ向かう単位方向ベクトル、|rij_0|は球iとjの初期結合長である。
【００３３】
曲げ変形に対する復元力Fijbは、次の手順にて算出する。まず、連続して一列に結合した3つの球i、j、kを考え、球iとjの中間点、および球jとkの中間点にそれぞれ作用する、曲げ変形に対する復元トルクTijb、Tkjbを次式にて算出する。
【００３４】
【数５】
2020064450.tif 000008
【００３５】
ただし、ρは3つの球の中心位置を通る円の曲率半径、θは方向ベクトルrijとrjkのなす角であり、
【００３６】
【数６】
2020064450.tif 000009
【００３７】
で与えられる。式(4)にて算出した復元トルクTijb、Tkjbを基に、球i、j、kに働く曲げ復元力を
【００３８】
【数７】
2020064450.tif 000010
【００３９】
と算出する。この操作を全ての結合した3つの球の組に対して行う。
【００４０】
粘性トルクに起因する力Fijhtの算出手順を以下に示す。まず、互いに結合した2つの球i、jを考え、この2つの球の中心位置に仮想的な球Imの存在を仮定する。この仮想球の回転角速度ωImを次式にて算出する。
【００４１】
【数８】
2020064450.tif 000011
【００４２】
算出した回転角速度ωImを基に、仮想球に働く粘性トルクThImを次式により求める。この粘性トルクに基づき、球iおよびjに作用する力FijhtおよびFjihtを次式にて算出する。
【００４３】
【数９】
2020064450.tif 000012
【００４４】
ただし、Ω(xIm)は仮想球の位置xImにおける流体の回転角速度ベクトルである。式(10)にて算出した粘性トルクに基づき、球iおよびjに作用する力FijhtおよびFjihtを次式にて算出する）。
【００４５】
【数１０】
2020064450.tif 000013
【００４６】
この操作を全ての結合した2つの球の組に対して行う。
【００４７】
＜粗視化繊維の解析手法＞
本発明の実施の形態では、上述した解析手法を基に、繊維を構成する球の数を低減させ、計算量の低減を図る。
【００４８】
粗視化繊維モデルはNseg個のセグメントから構成されていると考える。また、粗視化によって低減した、1セグメントあたりの球の個数をNRと表記する。
【００４９】
粗視化された繊維は球の数が低減しているため、そのままでは繊維としての質量及び慣性モーメントが実際よりも小さく、運動が正しく表現されない。そこで、これらの影響を考慮することで、非粗視化繊維モデルに対する粗視化繊維モデルの誤差を低減することを試みる。このために、粗視化された繊維の運動が元の繊維の運動と同等となるように定式化を変更することを考える。元の繊維を本発明の実施の形態では「非粗視化繊維モデル」と呼称する。非粗視化繊維モデルは、粗視化繊維モデルと等しい形状を有し、粗視化繊維モデルの各セグメント内において低減した球を考慮した繊維である。また、この繊維は定式化の際に仮想的に考える繊維であり、実際にその運動を解くわけではない。
【００５０】
上記の2種類の繊維について重心の運動方程式を考え、粗視化繊維モデルの運動方程式が、非粗視化繊維モデルの運動方程式と等しくなるよう、粗視化繊維モデルの各球に働く力を補正する。これを行うに当たり、以下の仮定を置く。
【００５１】
(a) セグメントは剛直である。
【００５２】
(b) 一本の繊維内のセグメントは全て等しい長さを持つ。
【００５３】
(c) セグメント内に分布して働く力（例えば粘性力）は、線形分布である。
【００５４】
(d) 一本の繊維内の球は全て等しい質量を持つ。
【００５５】
以上の仮定の下、粗視化繊維モデルおよび非粗視化繊維モデルの重心の運動方程式（並進および回転）を考える。このとき、球に作用する力を次の2つに分類する。
【００５６】
（１）一つ目は、セグメントの両端の球にのみ作用すると考えることのできる力(すなわち、セグメントに渡って分布していない非分布力)である。例えば、引張変形に対する引張復元力FS、曲げ変形に対する曲げ復元力Fbや、流体から受ける粘性トルクに起因する力Fhtである。また、剛体反発力FRepや壁面反発力Fwは、厳密にはセグメント内に分布する場合があるが、本発明の実施の形態では非分布力として取り扱う。
【００５７】
（２）二つ目は、セグメントに渡って分布して作用する力である。例えば、流体から受ける粘性力Fhである。
【００５８】
以降、上記一つ目の力をFND、二つ目の力をFDと表記する。さらに、表記を簡単にするため、球iから球i+1への単位方向ベクトルni i+1を次のように表記する。
【００５９】
【数１１】
2020064450.tif 000014
【００６０】
以上の前提の下で、粗視化繊維モデルに関する並進および回転の運動方程式は次式で与えられる。
【００６１】
ここで、
【００６２】
【数１２】
2020064450.tif 000015
【００６３】
ここで、vfCG、ωfCG、IfCGは粗視化繊維モデルの速度、回転角速度および慣性モーメントである。また、Ni0、N1、N2はそれぞれ次式で与えられる。
【００６４】
【数１３】
2020064450.tif 000016
【００６５】
一方、非粗視化繊維モデルに関する並進および回転の運動方程式は次式で表される。
【００６６】
【数１４】
2020064450.tif 000017
【００６７】
ここで、vf、ωf、Ifはそれぞれ非粗視化繊維モデルの速度、回転角速度および慣性モーメントである。
【００６８】
式(13)が式(17)と等しくなるよう、粗視化繊維における力を
【００６９】
【数１５】
2020064450.tif 000018
【００７０】
と補正することを考える。はじめに、繊維の並進の運動方程式（式(13)第1式および(17)第1式）に着目する。このときの補正係数をαTおよびβTと置くと、粗視化繊維モデルに関する式(13)第1式と非粗視化繊維モデルに関する(17)第1式とを比較することで、補正係数の計算式は次のように一意に導出される。
【００７１】
【数１６】
2020064450.tif 000019
【００７２】
ここで、iは球の番号を表す。同様に、繊維の回転の運動方程式（式(13)第2式および(17)第2式）に着目すると、補正係数αRおよびβRは次のように決定される。
【００７３】
【数１７】
2020064450.tif 000020
2020064450.tif 000021
【００７４】
ここで、
2020064450.tif 000022
である。ただし、riGは繊維の重心位置から球iの位置xiへ向かう方向ベクトルである。
【００７５】
しかしながら、繊維の回転運動を補正する場合と、並進運動を補正する場合とでは、補正係数の値が異なっており、繊維の並進および回転の両運動を同時に補正することは困難である。そこで、少なくとも繊維の回転運動に関しては正しく補正が行われるよう、次のように力の補正方法を変更する。まず、繊維の重心位置から球iの位置xiへ向かう単位方向ベクトルをniGと置き、力FNDおよびFDをそれぞれniGに平行な成分と垂直な成分に分解する。繊維の回転の運動方程式を考えた場合、平行成分は消え、垂直成分のみが残る。このため、垂直成分に対してはαRおよびβRを乗じ、平行成分に対してはαTおよびβTを乗じることとする。即ち、
【００７６】
【数１８】
2020064450.tif 000023
【００７７】
と補正を行う。
【００７８】
以上の手順をまとめると、毎時刻、以下の（Ａ）〜（Ｆ）の処理を繰り返すことになる。
【００７９】
(A) 各球に関して、補正係数を式(18)-(23)により算出する。
【００８０】
(B) 各球に関して、作用する力を式(2)-(11) により算出する。
【００８１】
(C) 各球に関して、(B)にて算出した力をFNDおよびFDに分類する。本発明の実施の形態では、粘性力をFD、その他の力をFNDとする。
【００８２】
(D) 各球に関して、FNDおよびFDを式(24)および(25)により補正する。
【００８３】
(E) 各球に関して、補正後の力
2020064450.tif 000024
の和を取り、補正後の力
2020064450.tif 000025
を算出する。
【００８４】
(F) 各球に関して、次の運動方程式を解いて、次時刻の各球の位置を算出する。
【００８５】
【数１９】
2020064450.tif 000026
【００８６】
ここで、補正係数の計算において、式(19),(21)は球毎に計算する必要があり、式(21),(22)は球の位置xiを用いて計算するため、毎時刻計算する必要がある。
【００８７】
＜本発明の実施の形態の繊維挙動計算装置の構成＞
図３に示すように、本発明の実施の形態に係る繊維挙動計算装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６、ＨＤＤ１８、通信インタフェース２１、及びこれらを相互に接続するためのバス２２を備えている。
【００８８】
ＣＰＵ１２は、各種プログラムを実行する。ＲＯＭ１４には、各種プログラムやパラメータ等が記憶されている。ＲＡＭ１６は、ＣＰＵ１２による各種プログラムの実行時におけるワークエリア等として用いられる。記録媒体としてのＨＤＤ１８には、後述する繊維挙動計算処理ルーチンを実行するためのプログラムを含む各種プログラムや各種データが記憶されている。
【００８９】
本実施の形態における繊維挙動計算装置１０を、繊維挙動計算処理ルーチンを実行するためのプログラムに沿って、機能ブロックで表すと、図４に示すようになる。繊維挙動計算装置１０は、入力部２０、演算部３０、及び出力部５０を備えている。
【００９０】
入力部２０は、解析対象となる粗視化繊維モデルに関する情報、及び流体に関する条件の設定値を入力として受け付ける。
【００９１】
演算部３０は、初期設定部３２、補正係数算出部３４、作用力算出部３６、力補正部４２、及び挙動計算部４４を備えている。なお、作用力算出部３６が、相互作用力算出部及び単体作用力算出部の一例である。
【００９２】
初期設定部３２は、解析対象となる流体中に含まれる繊維の集合の一部の繊維をモデル化した粗視化繊維モデルの各球の位置xiを初期化すると共に、各位置における流体の粘度η(ri)、速度ベクトルV(xi)、回転角速度ベクトルΩ(xIm)を設定する。
【００９３】
補正係数算出部３４は、各球の位置xiに基づいて、毎時刻、球毎に、式(18)-(23)により、当該球に作用する非分布力に対する補正係数αT、αR、分布力に対する補正係数βTi、βRiを計算する。なお、毎時刻が、所定時間間隔の一例である。また、補正係数の計算を毎時刻ではなく、数ステップに一回行うようにしてもよい。
【００９４】
作用力算出部３６は、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、式(1)-(11)より、当該球に作用する非分布力として、引張変形に対する引張復元力FS、曲げ変形に対する曲げ復元力Fb、流体から受ける粘性トルクに起因する力Fht、繊維同士の剛体反発力FRep、壁面反発力Fwを算出し、分布力として、流体から受ける粘性力Fhを算出する。
【００９５】
ここで、繊維同士の剛体反発力FRepを算出する際には、繊維モデルの繊維半径ａをer倍に拡大して算出し、引張変形に対する引張復元力FS、曲げ変形に対する曲げ復元力Fb、流体から受ける粘性トルクに起因する力Fht、壁面反発力Fw、流体から受ける粘性力Fhを算出する際には、繊維モデルの繊維半径ａを拡大せずに算出する。また、拡大倍率erは、予め与えられた、元の繊維の集合の体積含有率と、繊維の集合の一部の繊維をモデル化した粗視化繊維モデルの集合の体積含有率とに基づいて定められる。具体的には、元の繊維の集合の体積含有率Vfと、繊維数を減らした粗視化繊維モデルの集合の体積含有率Vf'とに基づいて、拡大倍率erは、以下の式に従って決定される。
【００９６】
【数２０】
2020064450.tif 000027
【００９７】
なお、流体に関する条件を、予め設定されたものとする場合を例に説明したが、これに限定されるものではなく、繊維の影響を受けて計算された流体に関する条件を用いてもよい。例えば、繊維の影響を受けて流れ場が変化することを考慮した、流体と繊維運動の連成計算を行うようにしてもよい。
【００９８】
力補正部４２は、毎時刻、球毎に、非分布力として扱う、引張変形に対する引張復元力FS、曲げ変形に対する曲げ復元力Fb、流体から受ける粘性トルクに起因する力Fht、繊維同士の剛体反発力FRep、壁面反発力Fwに対して、式（24）により補正係数αT、αRを乗じて補正する。
【００９９】
また、力補正部４２は、毎時刻、球毎に、分布力として扱う、流体から受ける粘性力Fhに対して、式（25）により補正係数αTi、αRiを乗じて補正する。
【０１００】
挙動計算部４４は、毎時刻、球毎に、補正された流体から受ける粘性力Fh、引張変形に対する引張復元力FS、曲げ変形に対する曲げ復元力Fb、流体から受ける粘性トルクに起因する力Fht、繊維同士の剛体反発力FRep、壁面反発力Fwの和を取って、補正後の力
2020064450.tif 000028
を算出し、式(26)により、速度ベクトルviを計算し、次時刻における当該球の位置ｘiを計算する。
【０１０１】
出力部５０は、毎時刻、球毎に計算された位置ｘiの計算結果を、解析対象となる粗視化繊維モデルの挙動として出力する。
【０１０２】
＜繊維挙動計算装置の動作＞
次に、本発明の実施の形態に係る繊維挙動計算装置１０の動作について説明する。
【０１０３】
入力部２０によって解析対象となる流体中に含まれる繊維の集合の一部の繊維をモデル化した粗視化繊維モデルに関する情報、及び流体に関する条件の設定値を受け付けると、繊維挙動計算装置１０によって、図５に示す繊維挙動計算処理ルーチンが実行される。
【０１０４】
まず、ステップＳ１００において、初期設定部３２は、解析対象となる流体中に含まれる繊維の集合の一部の繊維をモデル化した粗視化繊維モデルの各球の位置xiを初期化すると共に、各位置における流体の粘度η(ri)、速度ベクトルV(xi)、回転角速度ベクトルΩ(xIm)を設定する。
【０１０５】
ステップＳ１０２において、補正係数算出部３４は、各球の位置xiに基づいて、現時刻について、球毎に、式(18)-(23)により当該球に作用する非分布力に対する補正係数αT、αR、分布力に対する補正係数βTi、βRiを計算する。
【０１０６】
ステップＳ１０４において、作用力算出部３６は、予め設定された流体に関する条件と、各球の位置とに基づいて、現時刻について、球毎に、式(1)-(11)より、当該球に作用する非分布力として、引張変形に対する引張復元力FS、曲げ変形に対する曲げ復元力Fb、流体から受ける粘性トルクに起因する力Fht、繊維同士の剛体反発力FRep、壁面反発力Fwを算出し、分布力として、流体から受ける粘性力Fhを算出する。
【０１０７】
ステップＳ１０６において、力補正部４２は、現時刻について、球毎に、上記ステップＳ１０４で算出された、非分布力として扱う、引張変形に対する引張復元力FS、曲げ変形に対する曲げ復元力Fb、流体から受ける粘性トルクに起因する力Fht、繊維同士の剛体反発力FRep、壁面反発力Fwに対して、式（24）により補正係数αT、αRを乗じて補正する。
【０１０８】
また、力補正部４２は、現時刻について、球毎に、上記ステップＳ１０４で算出された、分布力として扱う、流体から受ける粘性力Fhに対して、式（25）により補正係数αTi、αRiを乗じて補正する。
【０１０９】
ステップＳ１０８において、挙動計算部４４は、現時刻について、球毎に、上記ステップＳ１０６で補正された流体から受ける粘性力Fh、引張変形に対する復元力FS、曲げ変形に対する復元力Fb、流体から受ける粘性トルクに起因する力Fht、繊維同士の剛体反発力FRep、壁面反発力Fwの和を取って、補正後の力
2020064450.tif 000029
を算出し、式(26)により、速度ベクトルviを計算し、次時刻における当該球の位置ｘiを計算する。
【０１１０】
ステップＳ１１０では、繰り返し処理を終了するか否かを判定する。繰り返し処理を終了しない場合には、上記ステップＳ１０２へ戻り、次時刻について処理を繰り返す。一方、繰り返し処理を終了すると判定された場合には、ステップＳ１１２において、出力部５０は、毎時刻、球毎に計算された位置ｘiの計算結果を、解析対象となる粗視化繊維モデルの挙動として出力し、繊維挙動計算処理ルーチンを終了する。
【０１１１】
＜実験結果＞
本発明の実施の形態で説明した手法を用いて数値解析を行った結果を示す。
【０１１２】
図６に示す金型形状内の流体中の繊維強化樹脂を解析対象とし、繊維強化樹脂の射出成形過程における繊維の運動解析を実施する。はじめに、樹脂流動解析を実施し、樹脂の流れ場を求める。ここで、樹脂流動解析の条件は、流量10cc/s、樹脂温度230℃、金型温度110℃とし、用いる樹脂材料は短ガラス繊維強化ポリプロピレンおよび長ガラス繊維強化ポリプロピレン相当とする。また、保圧は設定しない。
【０１１３】
次に、樹脂流動解析により得られた流れ場を用いて、短繊維の運動解析を行う。解析は、繊維長を200μm、繊維半径を5μm、繊維体積含有率Vfを10%、繊維を構成する球数を2とした条件（即ち、剛直繊維）にて実施する。得られた結果を用いて、上記図６に示す各位置（A、B、C）において算出した繊維配向テンソルを図７に示す。横軸は厚み方向位置を、縦軸は繊維配向テンソルのxおよびy方向成分（AxxおよびAyy）を表す。実線は従来技術（特許文献１、２に記載の、流体中を流れる繊維一本一本の運動をシミュレートする技術）による結果を表し、破線は本実施の形態の方法において半径拡大率をer=3とした結果である。両者は良く一致していることがわかる。また、半径拡大率erを変化させた場合の、従来技術に対する計算高速化割合を図８に示す。本実施の形態の方法の適用により、大幅な計算高速化が得られている。
【０１１４】
最後に、長繊維の運動解析を行う。解析は、繊維長を1000μm、繊維半径を10μm、繊維体積含有率Vfを10%とした条件にて実施する。また、繊維の曲げ変形が生じるよう、繊維を構成する球数を13とする。er=2とした場合の、位置A、B、Cにおける繊維配向テンソルを図９に、位置Cにおける繊維の様子を図１０に示す。図１０に示されるような、繊維の曲げ変形が生じる条件においても、本実施の形態の方法は従来技術と同等の結果を示すことがわかる。
【０１１５】
以上説明したように、本発明の実施の形態に係る繊維挙動計算装置によれば、繊維モデルの繊維半径を所定倍に拡大して、他の繊維モデルとの間で相互に作用する力を計算し、他の繊維モデルとの間で相互に作用する力、及び繊維モデルの単体に作用する力に基づいて、繊維の集合の一部の繊維をモデル化した繊維モデルの各々の挙動を計算することにより、繊維の集合の運動を一つ一つ計算する場合と比較して、計算に必要となる繊維の本数が低減されるため、計算時間を低減させつつ、同等の計算結果が得られる。これによって、高繊維体積含有率の条件での繊維運動解析が可能となる。
【０１１６】
また、アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、球に作用する力の補正係数を、毎時刻、球毎に計算することにより、Ｎ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる。
【０１１７】
また、粗視化繊維モデルを用いることにより、球の数が低減するため、それに伴い計算メモリおよび計算量が低減される。また、球が低減したことにより、非粗視化繊維と比べて質量や慣性モーメントが異なるため、本来は繊維の運動に誤差が生じる。本発明の実施の形態では、理論的に導出した補正係数を用いて粗視化繊維内の各球に働く力を補正するため、このような誤差が低減される。
【０１１８】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１１９】
例えば、上記の実施の形態では、Nseg個のセグメントから構成され、１セグメントあたりの球の個数をNR個に低減した粗視化繊維モデルを用いて、運動を解析する場合を例に説明したが、これに限定されるものではなく、別の粗視化繊維モデルを用いて、繊維の運動を解析するようにしてもよい。
【０１２０】
また、本発明のプログラムは、記憶媒体に格納して提供するようにしてもよい。
【符号の説明】
【０１２１】
１０ 繊維挙動計算装置
２０ 入力部
３０ 演算部
３２ 初期設定部
３４ 補正係数算出部
３６ 作用力算出部
４２ 力補正部
４４ 挙動計算部
５０ 出力部
【図１】
2020064450.tif 000030
【図２】
2020064450.tif 000031
【図３】
2020064450.tif 000032
【図４】
2020064450.tif 000033
【図５】
2020064450.tif 000034
【図６】
2020064450.tif 000035
【図７】
2020064450.tif 000036
【図８】
2020064450.tif 000037
【図９】
2020064450.tif 000038
【図１０】
2020064450.tif 000039