(57)【要約】【課題】繊維のアスペクト比に対応するＮ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られるようにする。【解決手段】補正係数算出部３４が、アスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、所定時間間隔で、球毎に、前記球に作用する力の補正係数を計算する。作用力算出部３６が、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、球に作用する力を計算する。挙動計算部４４が、毎時刻、球毎に、補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する。【選択図】図８
2019036119.tif 000002
【特許請求の範囲】
【請求項１】
流体中に含まれる繊維の挙動を計算する繊維挙動計算装置であって、
前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、所定時間間隔で、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部と、
予め設定され、又は計算される流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する作用力算出部と、
毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する挙動計算部と、を含み、
前記補正係数算出部は、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する
繊維挙動計算装置。
【請求項２】
前記粗視化繊維モデルは、前記ｎ個の球と、球と球との間を結ぶセグメントとを含み、
前記球に作用する力は、前記セグメントの両端の球に作用する力と、前記セグメントに渡って分布して作用する力を含み、
前記補正係数算出部は、各球の位置に基づいて、所定時間間隔で、球毎に、前記セグメントの両端の球に作用する力、及び前記セグメントに渡って分布して作用する力の各々に対する補正係数を計算し、
前記挙動計算部は、毎時刻、球毎に、前記補正係数を乗じた前記セグメントの両端の球に作用する力と、前記補正係数を乗じた前記セグメントに渡って分布して作用する力との和に基づいて、前記球の挙動を計算する請求項１記載の繊維挙動計算装置。
【請求項３】
流体中に含まれる繊維の挙動を計算する繊維挙動計算装置における繊維挙動計算方法であって、
補正係数算出部が、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、所定時間間隔で、球毎に、前記球に作用する力の補正係数を計算し、
作用力算出部が、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算し、
挙動計算部が、毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算することを含み、
前記補正係数算出部によって計算することでは、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する
繊維挙動計算方法。
【請求項４】
流体中に含まれる繊維の挙動を計算するためのプログラムであって、
コンピュータを、
前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、所定時間間隔で、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部、
予め設定され、又は計算される流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する作用力算出部、及び
毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する挙動計算部として機能させるためのプログラムであって、
前記補正係数算出部は、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する
プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、繊維挙動計算装置、方法、及びプログラムに係り、特に、流体中の繊維の挙動を計算する繊維挙動計算装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
従来より、繊維をＮ個未満の球の連結体としてモデル化し、計算量を低減する運動解析方法が知られている（特許文献１）。この運動解析方法では、球を低減したことにより繊維の運動に誤差が生じるため、球に働く力およびトルクに対して、経験的に決定した補正係数を乗じている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０００−３２２４０７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上記特許文献１に記載の運動解析方法では、以下の３つの問題点がある。
【０００５】
第一に、補正係数の値を単純な数値実験を基に決定している。第二に、補正係数の値は全ての球で同じ値を取る。第三に、補正係数の値は解析実施前に設定され、解析中は同じ値を取り続ける。
【０００６】
２つの球からなる曲がらない剛直な繊維の運動を解く場合等、限られた解析条件下ではこの補正方法は有効である。
【０００７】
しかし、補正係数は繊維の変形状態や球の個数に応じて変化するため、あらゆる繊維の状態に対して事前に補正係数を求めておくことは非常に困難である。すなわち。球の数が増えた場合や、繊維に曲げ変形が生じる場合には、補正係数の値を球毎に、逐次変更する必要がある。
【０００８】
このことから、従来技術で用いている補正係数の決定方法やその有用性は、汎用的であるとは言い難い。
【０００９】
本発明は、上記の事情を鑑みてなされたもので、繊維のアスペクト比に対応するＮ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる繊維挙動計算装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１０】
上記の目的を達成するために本発明に係る繊維挙動計算装置は、流体中に含まれる繊維の挙動を計算する繊維挙動計算装置であって、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、毎時刻、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部と、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する作用力算出部と、毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する挙動計算部と、を含み、前記補正係数算出部は、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する。
【００１１】
本発明に係る繊維挙動計算方法は、流体中に含まれる繊維の挙動を計算する繊維挙動計算装置における繊維挙動計算方法であって、補正係数算出部が、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、毎時刻、球毎に、前記球に作用する力の補正係数を計算し、作用力算出部が、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算し、挙動計算部が、毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算することを含み、前記補正係数算出部によって計算することでは、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する。
【００１２】
本発明に係るプログラムは、流体中に含まれる繊維の挙動を計算するためのプログラムであって、コンピュータを、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、毎時刻、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する作用力算出部、及び毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する挙動計算部として機能させるためのプログラムであって、前記補正係数算出部は、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算するプログラムである。
【００１３】
本発明によれば、補正係数算出部が、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、毎時刻、球毎に、前記球に作用する力の補正係数を計算する。このとき、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数が計算される。
【００１４】
また、作用力算出部が、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する。
【００１５】
そして、挙動計算部が、毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する。
【００１６】
このように、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記球に作用する力の補正係数を、毎時刻、球毎に計算することにより、Ｎ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる。
【発明の効果】
【００１７】
以上説明したように、本発明の繊維挙動計算装置及びプログラムによれば、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記球に作用する力の補正係数を、毎時刻、球毎に計算することにより、Ｎ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる、という効果が得られる。
【図面の簡単な説明】
【００１８】
【図１】繊維及び粗視化繊維モデルの一例を示す図である。
【図２】粗視化繊維モデル及び非粗視化繊維モデルの一例を示す図である。
【図３】繊維及び非粗視化繊維モデルの一例を示す図である。
【図４】（ａ）結合した３つの球を示す図、及び（ｂ）曲げ変形に対する復元力を説明するための図である。
【図５】（ａ）仮想球における回転角速度の算出方法を説明するための図、及び（ｂ）粘性トルクに起因する力の負荷を説明するための図である。
【図６】（ａ）粗視化繊維モデルの一例を示す図、及び（ｂ）非粗視化繊維モデルの一例を示す図である。
【図７】本発明の実施の形態に係る繊維挙動計算装置を示すブロック図である。
【図８】本発明の実施の形態に係る繊維挙動計算装置を示す機能ブロック図である。
【図９】本発明の実施の形態に係る繊維挙動計算装置の繊維挙動計算処理ルーチンの内容を示すフローチャートである。
【図１０】単純せん断流動場での挙動を説明するための図である。
【図１１】回転周期の誤差を示すグラフである。
【図１２】高速化割合を示すグラフである。
【図１３】ｔ＝０．４ｓにおける変形の様子を示す図である。
【図１４】端部の球の軌跡を示すグラフである。
【発明を実施するための形態】
【００１９】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００２０】
＜本発明の実施の形態の概要＞
本発明の実施の形態では、密に連結したＮ個の球を用いてモデル化される繊維を、図１に示すようにＮ個未満の球の連結体としてモデル化した粗視化繊維モデルを用いる。この粗視化繊維モデルの運動が、Ｎ個の球でモデル化した繊維である非粗視化繊維モデル（図２参照）の運動と同等になるように、粗視化繊維モデル内の各球に働く力に対して、毎時刻、球毎に異なる補正係数を乗じる。この補正係数は、粗視化繊維モデルおよび非粗視化繊維モデルの運動方程式をそれぞれ考え、前者が後者と等しくなるように理論的に導出され、導出した補正係数を乗じた力を用いて、粗視化繊維モデル内の球の運動方程式を解くことにより、位置、変形及び配向といった、繊維の挙動を解析する。
【００２１】
＜本発明の実施の形態の原理＞
＜非粗視化繊維の解析手法＞
まず、本発明の実施の形態において解く球の基礎方程式として、非粗視化繊維モデルについて、球に関する並進および回転の運動方程式の内、並進の運動方程式のみを解く解析手法について説明する。
【００２２】
図３に示すように、非粗視化繊維モデルでは、繊維を球の連結体としてモデル化する。このとき、球の数は繊維のアスペクト比と同数である。各球の運動方程式を解くことで、繊維の運動をシミュレートする。球iに関する運動方程式は
【００２３】
【数１】
2019036119.tif 000003
【００２４】
で与えられる。ただし、mは球の質量、viは球の速度、Fhは流体から受ける粘性力、FSは引張変形に対する復元力、Fbは曲げ変形に対する復元力、Fhtは流体から受ける粘性トルクに起因する力、FPは球同士の流体力学的相互作用力、FRは壁面反発力である。球同士の流体力学的相互作用力FPは例えば上記特許文献１に記載の方法にて算出する。壁面反発力FRは例えば壁面への球の貫入量に応じてペナルティ法等により算出する。また、jは球iに隣接して結合する球の番号を表し、nは球iと結合していない球の番号を表す。各力は次式にて算出される。
【数２】
2019036119.tif 000004
【００２５】
ただし、aは球の半径、Eは繊維の引張弾性率、η(ri)およびV(xi)は球の位置xiにおける流体の粘度および速度ベクトル、rijは球iから球jへ向かう方向ベクトル、nijは球iから球jへ向かう単位方向ベクトル、|rij_0|は球iとjの初期結合長である。
【００２６】
曲げ変形に対する復元力Fijbは、次の手順にて算出する。まず、連続して一列に結合した3つの球i、j、kを考え（図4(a)）、球iとjの中間点、および球jとkの中間点にそれぞれ作用する、曲げ変形に対する復元トルクTijb、Tkjbを次式にて算出する。
【００２７】
【数３】
2019036119.tif 000005
【００２８】
ただし、ρは3つの球の中心位置を通る円の曲率半径、θは方向ベクトルrijとrjkのなす角であり、
【００２９】
【数４】
2019036119.tif 000006
【００３０】
で与えられる。式(4)にて算出した復元トルクTijb、Tkjbを基に、球i、j、kに働く曲げ復元力を
【００３１】
【数５】
2019036119.tif 000007
【００３２】
と算出する（図4(b)）。この操作を全ての結合した3つの球の組に対して行う。
【００３３】
粘性トルクに起因する力Fijhtの算出手順を以下に示す。まず、互いに結合した2つの球i、jを考え、この2つの球の中心位置に仮想的な球Imの存在を仮定する（図5(a)）。この仮想球の回転角速度ωImを次式にて算出する。
【００３４】
【数６】
2019036119.tif 000008
【００３５】
算出した回転角速度ωImを基に、仮想球に働く粘性トルクThImを次式により求める。この粘性トルクに基づき、球iおよびjに作用する力FijhtおよびFjihtを次式にて算出する（図5(b)）。
【００３６】
【数７】
2019036119.tif 000009
【００３７】
ただし、Ω(xIm)は仮想球の位置xImにおける流体の渦度ベクトルである。式(10)にて算出した粘性トルクに基づき、球iおよびjに作用する力FijhtおよびFjihtを次式にて算出する（図5(b)）。
【００３８】
【数８】
2019036119.tif 000010
【００３９】
この操作を全ての結合した2つの球の組に対して行う。
【００４０】
＜粗視化繊維の解析手法＞
本発明の実施の形態では、上述した解析手法を基に、繊維を構成する球の数を低減させ、計算量の低減を図る。
【００４１】
粗視化繊維モデルはNseg個のセグメントから構成されていると考え、セグメントは2つの球から構成されるものとする。例として、図6(a)に示す粗視化繊維のセグメント数は3である。また、粗視化によって低減した、1セグメントあたりの球の個数をNRと表記する（図6(b)）。
【００４２】
粗視化された繊維は球の数が低減しているため、そのままでは繊維としての質量及び慣性モーメントが実際よりも小さく、運動が正しく表現されない。そこで、これらの影響を考慮することで、非粗視化繊維モデルに対する粗視化繊維モデルの誤差を低減することを試みる。このために、図6(a)に示す粗視化された繊維の運動が図6(b)に示す繊維の運動と同等となるように定式化を変更することを考える。図6(b)の繊維を本発明の実施の形態では「非粗視化繊維モデル」と呼称する。非粗視化繊維モデルは、粗視化繊維モデルと等しい形状を有し、粗視化繊維モデルの各セグメント内において低減した球を考慮した繊維である。また、この繊維は定式化の際に仮想的に考える繊維であり、実際にその運動を解くわけではない。
【００４３】
上記の2種類の繊維について重心の運動方程式を考え、粗視化繊維モデルの運動方程式が、非粗視化繊維モデルの運動方程式と等しくなるよう、粗視化繊維モデルの各球に働く力を補正する。これを行うに当たり、以下の仮定を置く。
【００４４】
(a) セグメントは剛直である。
【００４５】
(b) 一本の繊維内のセグメントは全て等しい長さを持つ。
【００４６】
(c) セグメント内に分布して働く力（例えば粘性力）は、線形分布である。
【００４７】
(d) 一本の繊維内の球は全て等しい質量を持つ。
【００４８】
以上の仮定の下、粗視化繊維モデルおよび非粗視化繊維モデルの重心の運動方程式（並進および回転）を考える。このとき、球に作用する力を次の2つに分類する。
【００４９】
（１）一つ目は、セグメントの両端の球にのみ作用すると考えることのできる力(すなわち、セグメントに渡って分布していない非分布力)である。例えば、引張変形に対する復元力FS、曲げ変形に対する復元力Fbや、流体から受ける粘性トルクに起因する力Fhtである。また、流体力学的相互作用力FPや壁面反発力FRは、厳密にはセグメント内に分布する場合があるが、本発明の実施の形態では非分布力として取り扱う。
【００５０】
（２）二つ目は、セグメントに渡って分布して作用する力である。例えば、流体から受ける粘性力Fhである。
【００５１】
以降、上記一つ目の力をFND、二つ目の力をFDと表記する。さらに、表記を簡単にするため、球iから球i+1への単位方向ベクトルni i+1を次のように表記する。
【００５２】
【数９】
2019036119.tif 000011
【００５３】
以上の前提の下で、粗視化繊維モデル（図6(a)）に関する並進および回転の運動方程式は次式で与えられる。
【００５４】
ここで、
【００５５】
【数１０】
2019036119.tif 000012
【００５６】
ここで、vfCG、ωfCG、IfCGは粗視化繊維モデルの速度、回転角速度および慣性モーメントである。また、Ni0、N1、N2はそれぞれ次式で与えられる。
【００５７】
【数１１】
2019036119.tif 000013
【００５８】
一方、非粗視化繊維モデル（図6(b)）に関する並進および回転の運動方程式は次式で表される。
【００５９】
【数１２】
2019036119.tif 000014
【００６０】
ここで、vf、ωf、Ifはそれぞれ非粗視化繊維モデルの速度、回転角速度および慣性モーメントである。
【００６１】
式(13)が式(17)と等しくなるよう、粗視化繊維における力を
【００６２】
【数１３】
2019036119.tif 000015
【００６３】
と補正することを考える。はじめに、繊維の並進の運動方程式（式(13)第1式および(17)第1式）に着目する。このときの補正係数をαTおよびβTと置くと、粗視化繊維モデルに関する式(13)第1式と非粗視化繊維モデルに関する(17)第1式とを比較することで、補正係数の計算式は次のように一意に導出される。
【００６４】
【数１４】
2019036119.tif 000016
【００６５】
ここで、iは球の番号を表す。同様に、繊維の回転の運動方程式（式(13)第2式および(17)第2式）に着目すると、補正係数αRおよびβRは次のように決定される。
【００６６】
【数１５】
2019036119.tif 000017
2019036119.tif 000018
ここで、
2019036119.tif 000019
【００６７】
である。ただし、riGは繊維の重心位置から球iの位置xiへ向かう方向ベクトルである。
【００６８】
しかしながら、繊維の回転運動を補正する場合と、並進運動を補正する場合とでは、補正係数の値が異なっており、繊維の並進および回転の両運動を同時に補正することは困難である。そこで、少なくとも繊維の回転運動に関しては正しく補正が行われるよう、次のように力の補正方法を変更する。まず、繊維の重心位置から球iの位置xiへ向かう単位方向ベクトルをniGと置き、力FNDおよびFDをそれぞれniGに平行な成分と垂直な成分に分解する。繊維の回転の運動方程式を考えた場合、平行成分は消え、垂直成分のみが残る。このため、垂直成分に対してはαRおよびβRを乗じ、平行成分に対してはαTおよびβTを乗じることとする。即ち、
【００６９】
【数１６】
2019036119.tif 000020
【００７０】
と補正を行う。
【００７１】
以上の手順をまとめると、毎時刻、以下の（Ａ）〜（Ｆ）の処理を繰り返すことになる。
【００７２】
(A) 各球に関して、補正係数を式(18)-(23)により算出する。
【００７３】
(B) 各球に関して、作用する力を式(2)-(11) により算出する。
【００７４】
(C) 各球に関して、(B)にて算出した力をFNDおよびFDに分類する。本発明の実施の形態では、粘性力をFD、その他の力をFNDとする。
【００７５】
(D) 各球に関して、FNDおよびFDを式(24)および(25)により補正する。
【００７６】
(E) 各球に関して、補正後の力
2019036119.tif 000021
の和を取り、補正後の力
2019036119.tif 000022
を算出する。
【００７７】
(F) 各球に関して、次の運動方程式を解いて、次時刻の各球の位置を算出する。
【００７８】
【数１７】
2019036119.tif 000023
【００７９】
ここで、補正係数の計算において、式(19),(21)は球毎に計算する必要があり、式(21),(22)は球の位置xiを用いて計算するため、毎時刻計算する必要がある。
【００８０】
＜本発明の実施の形態の繊維挙動計算装置の構成＞
図７に示すように、本発明の実施の形態に係る繊維挙動計算装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６、ＨＤＤ１８、通信インタフェース２１、及びこれらを相互に接続するためのバス２２を備えている。
【００８１】
ＣＰＵ１２は、各種プログラムを実行する。ＲＯＭ１４には、各種プログラムやパラメータ等が記憶されている。ＲＡＭ１６は、ＣＰＵ１２による各種プログラムの実行時におけるワークエリア等として用いられる。記録媒体としてのＨＤＤ１８には、後述する繊維挙動計算処理ルーチンを実行するためのプログラムを含む各種プログラムや各種データが記憶されている。
【００８２】
本実施の形態における繊維挙動計算装置１０を、繊維挙動計算処理ルーチンを実行するためのプログラムに沿って、機能ブロックで表すと、図８に示すようになる。繊維挙動計算装置１０は、入力部２０、演算部３０、及び出力部５０を備えている。
【００８３】
入力部２０は、解析対象となる粗視化繊維モデルに関する情報、及び流体に関する条件の設定値を入力として受け付ける。
【００８４】
演算部３０は、初期設定部３２、補正係数算出部３４、作用力算出部３６、力補正部４２、及び挙動計算部４４を備えている。
【００８５】
初期設定部３２は、解析対象となる粗視化繊維モデルの各球の位置xiを初期化すると共に、各位置における流体の粘度η(ri)、速度ベクトルV(xi)、渦度ベクトルΩ(xIm)を設定する。
【００８６】
補正係数算出部３４は、各球の位置xiに基づいて、毎時刻、球毎に、式(18)-(23)により、当該球に作用する非分布力に対する補正係数αT、αR、分布力に対する補正係数βTi、βRiを計算する。なお、毎時刻が、所定時間間隔の一例である。また、補正係数の計算を毎時刻ではなく、数ステップに一回行うようにしてもよい。
【００８７】
作用力算出部３６は、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、式(2)-(11)より、当該球に作用する非分布力として、引張変形に対する復元力FS、曲げ変形に対する復元力Fb、流体から受ける粘性トルクに起因する力Fht、球同士の流体力学的相互作用力FP、壁面反発力FRを算出し、分布力として、流体から受ける粘性力Fhを算出する。なお、流体に関する条件を、予め設定されたものとする場合を例に説明したが、これに限定されるものではなく、繊維の影響を受けて計算された流体に関する条件を用いてもよい。例えば、繊維の影響を受けて流れ場が変化することを考慮した、流体と繊維運動の連成計算を行うようにしてもよい。
【００８８】
力補正部４２は、毎時刻、球毎に、非分布力として扱う、引張変形に対する復元力FS、曲げ変形に対する復元力Fb、流体から受ける粘性トルクに起因する力Fht、球同士の流体力学的相互作用力FP、壁面反発力FRに対して、式（24）により補正係数αT、αRを乗じて補正する。
【００８９】
また、力補正部４２は、毎時刻、球毎に、分布力として扱う、流体から受ける粘性力Fhに対して、式（25）により補正係数αTi、αRiを乗じて補正する。
【００９０】
挙動計算部４４は、毎時刻、球毎に、補正された流体から受ける粘性力Fh、引張変形に対する復元力FS、曲げ変形に対する復元力Fb、流体から受ける粘性トルクに起因する力Fht、球同士の流体力学的相互作用力FP、壁面反発力Frの和を取って、補正後の力
2019036119.tif 000024
を算出し、式(26)により、速度ベクトルviを計算し、次時刻における当該球の位置ｘiを計算する。
【００９１】
出力部５０は、毎時刻、球毎に計算された位置ｘiの計算結果を、解析対象となる粗視化繊維モデルの挙動として出力する。
【００９２】
＜繊維挙動計算装置の動作＞
次に、本発明の実施の形態に係る繊維挙動計算装置１０の動作について説明する。
【００９３】
入力部２０によって解析対象となる粗視化繊維モデルに関する情報、及び流体に関する条件の設定値を受け付けると、繊維挙動計算装置１０によって、図９に示す繊維挙動計算処理ルーチンが実行される。
【００９４】
まず、ステップＳ１００において、初期設定部３２は、解析対象となる粗視化繊維モデルの各球の位置xiを初期化すると共に、各位置における流体の粘度η(ri)、速度ベクトルV(xi)、渦度ベクトルΩ(xIm)を設定する。
【００９５】
ステップＳ１０２において、補正係数算出部３４は、各球の位置xiに基づいて、現時刻について、球毎に、式(18)-(23)により当該球に作用する非分布力に対する補正係数αT、αR、分布力に対する補正係数βTi、βRiを計算する。
【００９６】
ステップＳ１０４において、作用力算出部３６は、予め設定された流体に関する条件と、各球の位置とに基づいて、現時刻について、球毎に、式(2)-(11)より、当該球に作用する非分布力として、引張変形に対する復元力FS、曲げ変形に対する復元力Fb、流体から受ける粘性トルクに起因する力Fht、球同士の流体力学的相互作用力FP、壁面反発力FRを算出し、分布力として、流体から受ける粘性力Fhを算出する。
【００９７】
ステップＳ１０６において、力補正部４２は、現時刻について、球毎に、上記ステップＳ１０４で算出された、非分布力として扱う、引張変形に対する復元力FS、曲げ変形に対する復元力Fb、流体から受ける粘性トルクに起因する力Fht、球同士の流体力学的相互作用力FP、壁面反発力FRに対して、式（24）により補正係数αT、αRを乗じて補正する。
【００９８】
また、力補正部４２は、現時刻について、球毎に、上記ステップＳ１０４で算出された、分布力として扱う、流体から受ける粘性力Fhに対して、式（25）により補正係数αTi、αRiを乗じて補正する。
【００９９】
ステップＳ１０８において、挙動計算部４４は、現時刻について、球毎に、上記ステップＳ１０６で補正された流体から受ける粘性力Fh、引張変形に対する復元力FS、曲げ変形に対する復元力Fb、流体から受ける粘性トルクに起因する力Fht、球同士の流体力学的相互作用力FP、壁面反発力Frの和を取って、補正後の力
2019036119.tif 000025
を算出し、式(26)により、速度ベクトルviを計算し、次時刻における当該球の位置ｘiを計算する。
【０１００】
ステップＳ１１０では、繰り返し処理を終了するか否かを判定する。繰り返し処理を終了しない場合には、上記ステップＳ１０２へ戻り、次時刻について処理を繰り返す。一方、繰り返し処理を終了すると判定された場合には、ステップＳ１１２において、出力部５０は、毎時刻、球毎に計算された位置ｘiの計算結果を、解析対象となる粗視化繊維モデルの挙動として出力し、繊維挙動計算処理ルーチンを終了する。
【０１０１】
＜実験結果＞
本発明の実施の形態で説明した手法を用いて数値解析を行った結果を示す。長さ500 μmの繊維を、半径a=5 μmの球の連結体としてモデル化した。このとき、非粗視化繊維モデルを用いた場合、繊維を構成する球の個数はN=50である。この繊維の中心座標を(x,y,z)=(0,0,0)とし、時刻t=0においてx軸に平行に配向しているものとする(図１０)。この繊維に対して、
2019036119.tif 000026
なる単純せん断流れ場を与える。ただし、u,v,wはそれぞれx,y,z方向の流速である。流れ場の粘度はη=100 Pa・s、せん断速度は
2019036119.tif 000027
とする。本計算条件下において、繊維の曲げ変形は小さく、剛体回転挙動を示す。繊維を構成する球の個数を低減した場合における、繊維の回転周期に関する、上記特許文献１に記載の技術と本実施の形態との誤差を図１１に示す。中実シンボルは補正処理を行わない場合、白抜きシンボルは補正処理を行った場合の結果を表す。補正処理を行わない場合、セグメント数が低下するにつれ、誤差が増加していくことがわかる。一方、補正処理を行うことで、誤差はほぼゼロとなり、非粗視化繊維モデルと同等の回転運動を示す。非粗視化繊維モデルでの計算時間に対する、粗視化モデルでの計算時間の比を図１２に示す。セグメント数が低下するにつれて、計算高速化率が上昇していることが分かる。
【０１０２】
次に、繊維が柔軟である場合での繊維の挙動を検証する。柔軟繊維を表現するため、弾性率E=0.5 GPaとして計算を行った。その他の条件は上記の検証と同じである。本計算条件では、繊維はS字状の変形を示す（図１３）。図１４に、繊維の端部の球のy座標の時間変化を示す。図１３、図１４より、セグメント数を15まで低下させた場合においても、粗視化繊維モデルは非粗視化繊維モデルの挙動を良く再現していることがわかる。
【０１０３】
以上説明したように、本発明の実施の形態に係る繊維挙動計算装置によれば、アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、球に作用する力の補正係数を、毎時刻、球毎に計算することにより、Ｎ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる。
【０１０４】
粗視化繊維モデルを用いることにより、球の数が低減するため、それに伴い計算メモリおよび計算量が低減される。また、球が低減したことにより、非粗視化繊維と比べて質量や慣性モーメントが異なるため、本来は繊維の運動に誤差が生じる。本発明の実施の形態では、理論的に導出した補正係数を用いて粗視化繊維内の各球に働く力を補正するため、このような誤差が低減される。
【０１０５】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１０６】
例えば、本発明のプログラムは、記憶媒体に格納して提供するようにしてもよい。
【符号の説明】
【０１０７】
１０   繊維挙動計算装置
２０   入力部
３０   演算部
３２   初期設定部
３４   補正係数算出部
３６   作用力算出部
４２   力補正部
４４   挙動計算部
５０   出力部
【図１】
2019036119.tif 000028
【図２】
2019036119.tif 000029
【図３】
2019036119.tif 000030
【図４】
2019036119.tif 000031
【図５】
2019036119.tif 000032
【図６】
2019036119.tif 000033
【図７】
2019036119.tif 000034
【図８】
2019036119.tif 000035
【図９】
2019036119.tif 000036
【図１０】
2019036119.tif 000037
【図１１】
2019036119.tif 000038
【図１２】
2019036119.tif 000039
【図１３】
2019036119.tif 000040
【図１４】
2019036119.tif 000041