特許 6711332 繊維挙動計算装置、方法、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6711332

(24)【登録日】2020年6月1日

(45)【発行日】2020年6月17日

(54)【発明の名称】繊維挙動計算装置、方法、及びプログラム

(51)【国際特許分類】

   G16Z 99/00 20190101AFI20200608BHJP

【ＦＩ】

   G16Z99/00

【請求項の数】4

【全頁数】19

(21)【出願番号】特願2017-156839(P2017-156839)

(22)【出願日】2017年8月15日

(65)【公開番号】特開2019-36119(P2019-36119A)

(43)【公開日】2019年3月7日

【審査請求日】2019年2月28日

(73)【特許権者】

【識別番号】000003609

【氏名又は名称】株式会社豊田中央研究所

(74)【代理人】

【識別番号】100079049

【弁理士】

【氏名又は名称】中島 淳

(74)【代理人】

【識別番号】100084995

【弁理士】

【氏名又は名称】加藤 和詳

(72)【発明者】

【氏名】笹山 俊貴

(72)【発明者】

【氏名】稲垣 昌英

【審査官】 宮地 匡人

(56)【参考文献】

【文献】 特開２０００−３２２４０７（ＪＰ，Ａ）

【文献】 特開平０５−３１４０９１（ＪＰ，Ａ）

【文献】 山本 智，繊維分散系のシミュレーション法の開発と応用(<小特集>数値計算と現実)，応用数理，２００４年１２月２２日，Vol.14 No.4，pp.323-333

【文献】 大野 将広，流体中の繊維状物質のモデリングに関する基礎的研究，日本機械学会関西支部講演会講演論文集，２００２年，Vo.77，pp.9-19 - 9-20

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ１６Ｚ ９９／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

流体中に含まれる繊維の挙動を計算する繊維挙動計算装置であって、
前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、所定時間間隔で、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部と、
予め設定され、又は計算される流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する作用力算出部と、
毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する挙動計算部と、を含み、
前記補正係数算出部は、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する
繊維挙動計算装置。

【請求項2】

前記粗視化繊維モデルは、前記ｎ個の球と、球と球との間を結ぶセグメントとを含み、
前記球に作用する力は、前記セグメントの両端の球に作用する力と、前記セグメントに渡って分布して作用する力を含み、
前記補正係数算出部は、各球の位置に基づいて、所定時間間隔で、球毎に、前記セグメントの両端の球に作用する力、及び前記セグメントに渡って分布して作用する力の各々に対する補正係数を計算し、
前記挙動計算部は、毎時刻、球毎に、前記補正係数を乗じた前記セグメントの両端の球に作用する力と、前記補正係数を乗じた前記セグメントに渡って分布して作用する力との和に基づいて、前記球の挙動を計算する請求項１記載の繊維挙動計算装置。

【請求項3】

流体中に含まれる繊維の挙動を計算する繊維挙動計算装置における繊維挙動計算方法であって、
補正係数算出部が、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、所定時間間隔で、球毎に、前記球に作用する力の補正係数を計算し、
作用力算出部が、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算し、
挙動計算部が、毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算することを含み、
前記補正係数算出部によって計算することでは、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する
繊維挙動計算方法。

【請求項4】

流体中に含まれる繊維の挙動を計算するためのプログラムであって、
コンピュータを、
前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、所定時間間隔で、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部、
予め設定され、又は計算される流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する作用力算出部、及び
毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する挙動計算部として機能させるためのプログラムであって、
前記補正係数算出部は、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する
プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、繊維挙動計算装置、方法、及びプログラムに係り、特に、流体中の繊維の挙動を計算する繊維挙動計算装置、方法、及びプログラムに関する。

【背景技術】

【0002】

従来より、繊維をＮ個未満の球の連結体としてモデル化し、計算量を低減する運動解析方法が知られている（特許文献１）。この運動解析方法では、球を低減したことにより繊維の運動に誤差が生じるため、球に働く力およびトルクに対して、経験的に決定した補正係数を乗じている。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特開２０００−３２２４０７号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

しかしながら、上記特許文献１に記載の運動解析方法では、以下の３つの問題点がある。

【0005】

第一に、補正係数の値を単純な数値実験を基に決定している。第二に、補正係数の値は全ての球で同じ値を取る。第三に、補正係数の値は解析実施前に設定され、解析中は同じ値を取り続ける。

【0006】

２つの球からなる曲がらない剛直な繊維の運動を解く場合等、限られた解析条件下ではこの補正方法は有効である。

【0007】

しかし、補正係数は繊維の変形状態や球の個数に応じて変化するため、あらゆる繊維の状態に対して事前に補正係数を求めておくことは非常に困難である。すなわち。球の数が増えた場合や、繊維に曲げ変形が生じる場合には、補正係数の値を球毎に、逐次変更する必要がある。

【0008】

このことから、従来技術で用いている補正係数の決定方法やその有用性は、汎用的であるとは言い難い。

【0009】

本発明は、上記の事情を鑑みてなされたもので、繊維のアスペクト比に対応するＮ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる繊維挙動計算装置、方法、及びプログラムを提供することを目的とする。

【課題を解決するための手段】

【0010】

上記の目的を達成するために本発明に係る繊維挙動計算装置は、流体中に含まれる繊維の挙動を計算する繊維挙動計算装置であって、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、毎時刻、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部と、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する作用力算出部と、毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する挙動計算部と、を含み、前記補正係数算出部は、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する。

【0011】

本発明に係る繊維挙動計算方法は、流体中に含まれる繊維の挙動を計算する繊維挙動計算装置における繊維挙動計算方法であって、補正係数算出部が、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、毎時刻、球毎に、前記球に作用する力の補正係数を計算し、作用力算出部が、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算し、挙動計算部が、毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算することを含み、前記補正係数算出部によって計算することでは、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算する。

【0012】

本発明に係るプログラムは、流体中に含まれる繊維の挙動を計算するためのプログラムであって、コンピュータを、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、毎時刻、球毎に、前記球に作用する力の補正係数を計算する補正係数算出部、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する作用力算出部、及び毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する挙動計算部として機能させるためのプログラムであって、前記補正係数算出部は、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数を計算するプログラムである。

【0013】

本発明によれば、補正係数算出部が、前記繊維のアスペクト比に対応するＮ個より少ないｎ個の球でモデル化された粗視化繊維モデルを用いて、各球の位置に基づいて、毎時刻、球毎に、前記球に作用する力の補正係数を計算する。このとき、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記補正係数が計算される。

【0014】

また、作用力算出部が、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、前記球に作用する力を計算する。

【0015】

そして、挙動計算部が、毎時刻、球毎に、前記補正係数を乗じた前記球に作用する力に基づいて、前記球の挙動を計算する。

【0016】

このように、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記球に作用する力の補正係数を、毎時刻、球毎に計算することにより、Ｎ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる。

【発明の効果】

【0017】

以上説明したように、本発明の繊維挙動計算装置及びプログラムによれば、前記アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、前記球に作用する力の補正係数を、毎時刻、球毎に計算することにより、Ｎ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる、という効果が得られる。

【図面の簡単な説明】

【0018】

【図1】繊維及び粗視化繊維モデルの一例を示す図である。

【図2】粗視化繊維モデル及び非粗視化繊維モデルの一例を示す図である。

【図3】繊維及び非粗視化繊維モデルの一例を示す図である。

【図4】（ａ）結合した３つの球を示す図、及び（ｂ）曲げ変形に対する復元力を説明するための図である。

【図5】（ａ）仮想球における回転角速度の算出方法を説明するための図、及び（ｂ）粘性トルクに起因する力の負荷を説明するための図である。

【図6】（ａ）粗視化繊維モデルの一例を示す図、及び（ｂ）非粗視化繊維モデルの一例を示す図である。

【図7】本発明の実施の形態に係る繊維挙動計算装置を示すブロック図である。

【図8】本発明の実施の形態に係る繊維挙動計算装置を示す機能ブロック図である。

【図9】本発明の実施の形態に係る繊維挙動計算装置の繊維挙動計算処理ルーチンの内容を示すフローチャートである。

【図10】単純せん断流動場での挙動を説明するための図である。

【図11】回転周期の誤差を示すグラフである。

【図12】高速化割合を示すグラフである。

【図13】ｔ＝０．４ｓにおける変形の様子を示す図である。

【図14】端部の球の軌跡を示すグラフである。

【発明を実施するための形態】

【0019】

以下、図面を参照して本発明の実施の形態を詳細に説明する。

【0020】

＜本発明の実施の形態の概要＞
本発明の実施の形態では、密に連結したＮ個の球を用いてモデル化される繊維を、図１に示すようにＮ個未満の球の連結体としてモデル化した粗視化繊維モデルを用いる。この粗視化繊維モデルの運動が、Ｎ個の球でモデル化した繊維である非粗視化繊維モデル（図２参照）の運動と同等になるように、粗視化繊維モデル内の各球に働く力に対して、毎時刻、球毎に異なる補正係数を乗じる。この補正係数は、粗視化繊維モデルおよび非粗視化繊維モデルの運動方程式をそれぞれ考え、前者が後者と等しくなるように理論的に導出され、導出した補正係数を乗じた力を用いて、粗視化繊維モデル内の球の運動方程式を解くことにより、位置、変形及び配向といった、繊維の挙動を解析する。

【0021】

＜本発明の実施の形態の原理＞
＜非粗視化繊維の解析手法＞
まず、本発明の実施の形態において解く球の基礎方程式として、非粗視化繊維モデルについて、球に関する並進および回転の運動方程式の内、並進の運動方程式のみを解く解析手法について説明する。

【0022】

図３に示すように、非粗視化繊維モデルでは、繊維を球の連結体としてモデル化する。このとき、球の数は繊維のアスペクト比と同数である。各球の運動方程式を解くことで、繊維の運動をシミュレートする。球iに関する運動方程式は

【0023】

【数1】

【0024】

で与えられる。ただし、mは球の質量、v_iは球の速度、F^hは流体から受ける粘性力、F^Sは引張変形に対する復元力、F^bは曲げ変形に対する復元力、F^htは流体から受ける粘性トルクに起因する力、F^Pは球同士の流体力学的相互作用力、F^Rは壁面反発力である。球同士の流体力学的相互作用力F^Pは例えば上記特許文献１に記載の方法にて算出する。壁面反発力F^Rは例えば壁面への球の貫入量に応じてペナルティ法等により算出する。また、jは球iに隣接して結合する球の番号を表し、nは球iと結合していない球の番号を表す。各力は次式にて算出される。

【数2】

【0025】

ただし、aは球の半径、Eは繊維の引張弾性率、η(r_i)およびV(x_i)は球の位置x_iにおける流体の粘度および速度ベクトル、r_ijは球iから球jへ向かう方向ベクトル、n_ijは球iから球jへ向かう単位方向ベクトル、|r_{ij_0}|は球iとjの初期結合長である。

【0026】

曲げ変形に対する復元力F_ij^bは、次の手順にて算出する。まず、連続して一列に結合した3つの球i、j、kを考え（図4(a)）、球iとjの中間点、および球jとkの中間点にそれぞれ作用する、曲げ変形に対する復元トルクT_ij^b、T_kj^bを次式にて算出する。

【0027】

【数3】

【0028】

ただし、ρは3つの球の中心位置を通る円の曲率半径、θは方向ベクトルr_ijとr_jkのなす角であり、

【0029】

【数4】

【0030】

で与えられる。式(4)にて算出した復元トルクT_ij^b、T_kj^bを基に、球i、j、kに働く曲げ復元力を

【0031】

【数5】

【0032】

と算出する（図4(b)）。この操作を全ての結合した3つの球の組に対して行う。

【0033】

粘性トルクに起因する力F_ij^htの算出手順を以下に示す。まず、互いに結合した2つの球i、jを考え、この2つの球の中心位置に仮想的な球Imの存在を仮定する（図5(a)）。この仮想球の回転角速度ω_Imを次式にて算出する。

【0034】

【数6】

【0035】

算出した回転角速度ω_Imを基に、仮想球に働く粘性トルクT^h_Imを次式により求める。この粘性トルクに基づき、球iおよびjに作用する力F_ij^htおよびF_ji^htを次式にて算出する（図5(b)）。

【0036】

【数7】

【0037】

ただし、Ω(x_Im)は仮想球の位置x_Imにおける流体の渦度ベクトルである。式(10)にて算出した粘性トルクに基づき、球iおよびjに作用する力F_ij^htおよびF_ji^htを次式にて算出する（図5(b)）。

【0038】

【数8】

【0039】

この操作を全ての結合した2つの球の組に対して行う。

【0040】

＜粗視化繊維の解析手法＞
本発明の実施の形態では、上述した解析手法を基に、繊維を構成する球の数を低減させ、計算量の低減を図る。

【0041】

粗視化繊維モデルはN^seg個のセグメントから構成されていると考え、セグメントは2つの球から構成されるものとする。例として、図6(a)に示す粗視化繊維のセグメント数は3である。また、粗視化によって低減した、1セグメントあたりの球の個数をN^Rと表記する（図6(b)）。

【0042】

粗視化された繊維は球の数が低減しているため、そのままでは繊維としての質量及び慣性モーメントが実際よりも小さく、運動が正しく表現されない。そこで、これらの影響を考慮することで、非粗視化繊維モデルに対する粗視化繊維モデルの誤差を低減することを試みる。このために、図6(a)に示す粗視化された繊維の運動が図6(b)に示す繊維の運動と同等となるように定式化を変更することを考える。図6(b)の繊維を本発明の実施の形態では「非粗視化繊維モデル」と呼称する。非粗視化繊維モデルは、粗視化繊維モデルと等しい形状を有し、粗視化繊維モデルの各セグメント内において低減した球を考慮した繊維である。また、この繊維は定式化の際に仮想的に考える繊維であり、実際にその運動を解くわけではない。

【0043】

上記の2種類の繊維について重心の運動方程式を考え、粗視化繊維モデルの運動方程式が、非粗視化繊維モデルの運動方程式と等しくなるよう、粗視化繊維モデルの各球に働く力を補正する。これを行うに当たり、以下の仮定を置く。

【0044】

(a) セグメントは剛直である。

【0045】

(b) 一本の繊維内のセグメントは全て等しい長さを持つ。

【0046】

(c) セグメント内に分布して働く力（例えば粘性力）は、線形分布である。

【0047】

(d) 一本の繊維内の球は全て等しい質量を持つ。

【0048】

以上の仮定の下、粗視化繊維モデルおよび非粗視化繊維モデルの重心の運動方程式（並進および回転）を考える。このとき、球に作用する力を次の2つに分類する。

【0049】

（１）一つ目は、セグメントの両端の球にのみ作用すると考えることのできる力(すなわち、セグメントに渡って分布していない非分布力)である。例えば、引張変形に対する復元力F^S、曲げ変形に対する復元力F^bや、流体から受ける粘性トルクに起因する力F^htである。また、流体力学的相互作用力F^Pや壁面反発力F^Rは、厳密にはセグメント内に分布する場合があるが、本発明の実施の形態では非分布力として取り扱う。

【0050】

（２）二つ目は、セグメントに渡って分布して作用する力である。例えば、流体から受ける粘性力F^hである。

【0051】

以降、上記一つ目の力をF^ND、二つ目の力をF^Dと表記する。さらに、表記を簡単にするため、球iから球i+1への単位方向ベクトルn_{i i+1}を次のように表記する。

【0052】

【数9】

【0053】

以上の前提の下で、粗視化繊維モデル（図6(a)）に関する並進および回転の運動方程式は次式で与えられる。

【0054】

ここで、

【0055】

【数10】

【0056】

ここで、v_f^CG、ω_f^CG、I_f^CGは粗視化繊維モデルの速度、回転角速度および慣性モーメントである。また、N_i⁰、N¹、N²はそれぞれ次式で与えられる。

【0057】

【数11】

【0058】

一方、非粗視化繊維モデル（図6(b)）に関する並進および回転の運動方程式は次式で表される。

【0059】

【数12】

【0060】

ここで、v_f、ω_f、I_fはそれぞれ非粗視化繊維モデルの速度、回転角速度および慣性モーメントである。

【0061】

式(13)が式(17)と等しくなるよう、粗視化繊維における力を

【0062】

【数13】

【0063】

と補正することを考える。はじめに、繊維の並進の運動方程式（式(13)第1式および(17)第1式）に着目する。このときの補正係数をα^Tおよびβ^Tと置くと、粗視化繊維モデルに関する式(13)第1式と非粗視化繊維モデルに関する(17)第1式とを比較することで、補正係数の計算式は次のように一意に導出される。

【0064】

【数14】

【0065】

ここで、iは球の番号を表す。同様に、繊維の回転の運動方程式（式(13)第2式および(17)第2式）に着目すると、補正係数α^Rおよびβ^Rは次のように決定される。

【0066】

【数15】

ここで、

【0067】

である。ただし、r_i^Gは繊維の重心位置から球iの位置x_iへ向かう方向ベクトルである。

【0068】

しかしながら、繊維の回転運動を補正する場合と、並進運動を補正する場合とでは、補正係数の値が異なっており、繊維の並進および回転の両運動を同時に補正することは困難である。そこで、少なくとも繊維の回転運動に関しては正しく補正が行われるよう、次のように力の補正方法を変更する。まず、繊維の重心位置から球iの位置x_iへ向かう単位方向ベクトルをn_i^Gと置き、力F^NDおよびF^Dをそれぞれn_i^Gに平行な成分と垂直な成分に分解する。繊維の回転の運動方程式を考えた場合、平行成分は消え、垂直成分のみが残る。このため、垂直成分に対してはα^Rおよびβ^Rを乗じ、平行成分に対してはα^Tおよびβ^Tを乗じることとする。即ち、

【0069】

【数16】

【0070】

と補正を行う。

【0071】

以上の手順をまとめると、毎時刻、以下の（Ａ）〜（Ｆ）の処理を繰り返すことになる。

【0072】

(A) 各球に関して、補正係数を式(18)-(23)により算出する。

【0073】

(B) 各球に関して、作用する力を式(2)-(11) により算出する。

【0074】

(C) 各球に関して、(B)にて算出した力をF^NDおよびF^Dに分類する。本発明の実施の形態では、粘性力をF^D、その他の力をF^NDとする。

【0075】

(D) 各球に関して、F^NDおよびF^Dを式(24)および(25)により補正する。

【0076】

(E) 各球に関して、補正後の力

の和を取り、補正後の力

を算出する。

【0077】

(F) 各球に関して、次の運動方程式を解いて、次時刻の各球の位置を算出する。

【0078】

【数17】

【0079】

ここで、補正係数の計算において、式(19),(21)は球毎に計算する必要があり、式(21),(22)は球の位置x_iを用いて計算するため、毎時刻計算する必要がある。

【0080】

＜本発明の実施の形態の繊維挙動計算装置の構成＞
図７に示すように、本発明の実施の形態に係る繊維挙動計算装置１０は、ＣＰＵ１２、ＲＯＭ１４、ＲＡＭ１６、ＨＤＤ１８、通信インタフェース２１、及びこれらを相互に接続するためのバス２２を備えている。

【0081】

ＣＰＵ１２は、各種プログラムを実行する。ＲＯＭ１４には、各種プログラムやパラメータ等が記憶されている。ＲＡＭ１６は、ＣＰＵ１２による各種プログラムの実行時におけるワークエリア等として用いられる。記録媒体としてのＨＤＤ１８には、後述する繊維挙動計算処理ルーチンを実行するためのプログラムを含む各種プログラムや各種データが記憶されている。

【0082】

本実施の形態における繊維挙動計算装置１０を、繊維挙動計算処理ルーチンを実行するためのプログラムに沿って、機能ブロックで表すと、図８に示すようになる。繊維挙動計算装置１０は、入力部２０、演算部３０、及び出力部５０を備えている。

【0083】

入力部２０は、解析対象となる粗視化繊維モデルに関する情報、及び流体に関する条件の設定値を入力として受け付ける。

【0084】

演算部３０は、初期設定部３２、補正係数算出部３４、作用力算出部３６、力補正部４２、及び挙動計算部４４を備えている。

【0085】

初期設定部３２は、解析対象となる粗視化繊維モデルの各球の位置x_iを初期化すると共に、各位置における流体の粘度η(r_i)、速度ベクトルV(x_i)、渦度ベクトルΩ(x_Im)を設定する。

【0086】

補正係数算出部３４は、各球の位置x_iに基づいて、毎時刻、球毎に、式(18)-(23)により、当該球に作用する非分布力に対する補正係数α^T、α^R、分布力に対する補正係数β^T_i、β^R_iを計算する。なお、毎時刻が、所定時間間隔の一例である。また、補正係数の計算を毎時刻ではなく、数ステップに一回行うようにしてもよい。

【0087】

作用力算出部３６は、予め設定された流体に関する条件と、各球の位置とに基づいて、毎時刻、球毎に、式(2)-(11)より、当該球に作用する非分布力として、引張変形に対する復元力F^S、曲げ変形に対する復元力F^b、流体から受ける粘性トルクに起因する力F^ht、球同士の流体力学的相互作用力F^P、壁面反発力F^Rを算出し、分布力として、流体から受ける粘性力F^hを算出する。なお、流体に関する条件を、予め設定されたものとする場合を例に説明したが、これに限定されるものではなく、繊維の影響を受けて計算された流体に関する条件を用いてもよい。例えば、繊維の影響を受けて流れ場が変化することを考慮した、流体と繊維運動の連成計算を行うようにしてもよい。

【0088】

力補正部４２は、毎時刻、球毎に、非分布力として扱う、引張変形に対する復元力F^S、曲げ変形に対する復元力F^b、流体から受ける粘性トルクに起因する力F^ht、球同士の流体力学的相互作用力F^P、壁面反発力F^Rに対して、式（24）により補正係数α^T、α^Rを乗じて補正する。

【0089】

また、力補正部４２は、毎時刻、球毎に、分布力として扱う、流体から受ける粘性力F^hに対して、式（25）により補正係数α^T_i、α^R_iを乗じて補正する。

【0090】

挙動計算部４４は、毎時刻、球毎に、補正された流体から受ける粘性力F^h、引張変形に対する復元力F^S、曲げ変形に対する復元力F^b、流体から受ける粘性トルクに起因する力F^ht、球同士の流体力学的相互作用力F^P、壁面反発力Frの和を取って、補正後の力

を算出し、式(26)により、速度ベクトルv_iを計算し、次時刻における当該球の位置ｘ_iを計算する。

【0091】

出力部５０は、毎時刻、球毎に計算された位置ｘ_iの計算結果を、解析対象となる粗視化繊維モデルの挙動として出力する。

【0092】

＜繊維挙動計算装置の動作＞
次に、本発明の実施の形態に係る繊維挙動計算装置１０の動作について説明する。

【0093】

入力部２０によって解析対象となる粗視化繊維モデルに関する情報、及び流体に関する条件の設定値を受け付けると、繊維挙動計算装置１０によって、図９に示す繊維挙動計算処理ルーチンが実行される。

【0094】

まず、ステップＳ１００において、初期設定部３２は、解析対象となる粗視化繊維モデルの各球の位置x_iを初期化すると共に、各位置における流体の粘度η(r_i)、速度ベクトルV(x_i)、渦度ベクトルΩ(x_Im)を設定する。

【0095】

ステップＳ１０２において、補正係数算出部３４は、各球の位置x_iに基づいて、現時刻について、球毎に、式(18)-(23)により当該球に作用する非分布力に対する補正係数α^T、α^R、分布力に対する補正係数β^T_i、β^R_iを計算する。

【0096】

ステップＳ１０４において、作用力算出部３６は、予め設定された流体に関する条件と、各球の位置とに基づいて、現時刻について、球毎に、式(2)-(11)より、当該球に作用する非分布力として、引張変形に対する復元力F^S、曲げ変形に対する復元力F^b、流体から受ける粘性トルクに起因する力F^ht、球同士の流体力学的相互作用力F^P、壁面反発力F^Rを算出し、分布力として、流体から受ける粘性力F^hを算出する。

【0097】

ステップＳ１０６において、力補正部４２は、現時刻について、球毎に、上記ステップＳ１０４で算出された、非分布力として扱う、引張変形に対する復元力F^S、曲げ変形に対する復元力F^b、流体から受ける粘性トルクに起因する力F^ht、球同士の流体力学的相互作用力F^P、壁面反発力F^Rに対して、式（24）により補正係数α^T、α^Rを乗じて補正する。

【0098】

また、力補正部４２は、現時刻について、球毎に、上記ステップＳ１０４で算出された、分布力として扱う、流体から受ける粘性力F^hに対して、式（25）により補正係数α^T_i、α^R_iを乗じて補正する。

【0099】

ステップＳ１０８において、挙動計算部４４は、現時刻について、球毎に、上記ステップＳ１０６で補正された流体から受ける粘性力F^h、引張変形に対する復元力F^S、曲げ変形に対する復元力F^b、流体から受ける粘性トルクに起因する力F^ht、球同士の流体力学的相互作用力F^P、壁面反発力Frの和を取って、補正後の力

を算出し、式(26)により、速度ベクトルv_iを計算し、次時刻における当該球の位置ｘ_iを計算する。

【0100】

ステップＳ１１０では、繰り返し処理を終了するか否かを判定する。繰り返し処理を終了しない場合には、上記ステップＳ１０２へ戻り、次時刻について処理を繰り返す。一方、繰り返し処理を終了すると判定された場合には、ステップＳ１１２において、出力部５０は、毎時刻、球毎に計算された位置ｘ_iの計算結果を、解析対象となる粗視化繊維モデルの挙動として出力し、繊維挙動計算処理ルーチンを終了する。

【0101】

＜実験結果＞
本発明の実施の形態で説明した手法を用いて数値解析を行った結果を示す。長さ500 μmの繊維を、半径a=5 μmの球の連結体としてモデル化した。このとき、非粗視化繊維モデルを用いた場合、繊維を構成する球の個数はN=50である。この繊維の中心座標を(x,y,z)=(0,0,0)とし、時刻t=0においてx軸に平行に配向しているものとする(図１０)。この繊維に対して、

なる単純せん断流れ場を与える。ただし、u,v,wはそれぞれx,y,z方向の流速である。流れ場の粘度はη=100 Pa・s、せん断速度は

とする。本計算条件下において、繊維の曲げ変形は小さく、剛体回転挙動を示す。繊維を構成する球の個数を低減した場合における、繊維の回転周期に関する、上記特許文献１に記載の技術と本実施の形態との誤差を図１１に示す。中実シンボルは補正処理を行わない場合、白抜きシンボルは補正処理を行った場合の結果を表す。補正処理を行わない場合、セグメント数が低下するにつれ、誤差が増加していくことがわかる。一方、補正処理を行うことで、誤差はほぼゼロとなり、非粗視化繊維モデルと同等の回転運動を示す。非粗視化繊維モデルでの計算時間に対する、粗視化モデルでの計算時間の比を図１２に示す。セグメント数が低下するにつれて、計算高速化率が上昇していることが分かる。

【0102】

次に、繊維が柔軟である場合での繊維の挙動を検証する。柔軟繊維を表現するため、弾性率E=0.5 GPaとして計算を行った。その他の条件は上記の検証と同じである。本計算条件では、繊維はS字状の変形を示す（図１３）。図１４に、繊維の端部の球のy座標の時間変化を示す。図１３、図１４より、セグメント数を15まで低下させた場合においても、粗視化繊維モデルは非粗視化繊維モデルの挙動を良く再現していることがわかる。

【0103】

以上説明したように、本発明の実施の形態に係る繊維挙動計算装置によれば、アスペクト比に対応するＮ個の球でモデル化された非粗視化繊維モデルの運動方程式と、粗視化繊維モデルの運動方程式とを比較した結果に基づいて導出された計算式により、球に作用する力の補正係数を、毎時刻、球毎に計算することにより、Ｎ個の球の運動を解く場合と比較して、計算時間を低減させつつ、同等の計算結果が得られる。

【0104】

粗視化繊維モデルを用いることにより、球の数が低減するため、それに伴い計算メモリおよび計算量が低減される。また、球が低減したことにより、非粗視化繊維と比べて質量や慣性モーメントが異なるため、本来は繊維の運動に誤差が生じる。本発明の実施の形態では、理論的に導出した補正係数を用いて粗視化繊維内の各球に働く力を補正するため、このような誤差が低減される。

【0105】

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

【0106】

例えば、本発明のプログラムは、記憶媒体に格納して提供するようにしてもよい。

【符号の説明】

【0107】

１０ 繊維挙動計算装置
２０ 入力部
３０ 演算部
３２ 初期設定部
３４ 補正係数算出部
３６ 作用力算出部
４２ 力補正部
４４ 挙動計算部
５０ 出力部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】