(57)【特許請求の範囲】
【請求項１】
流体中を移動する円柱状繊維が、繊維の長軸方向において収縮する方向への流れである収縮流中におかれたときに座屈する条件である座屈条件を導出する方法であって、
円柱状繊維の縦横比ｒ０’に対して一意的に求まる無次元の流体応力分布に、円柱が存在している場所でのμｋ（μは流体粘度、ｋは収縮速度）を掛けることによって収縮流中の円柱表面に働く流体応力分布を得る工程と、
円柱表面での流体応力分布から導かれる座屈する閾値である最小固有値λ０を用いることにより、座屈条件を導出する工程と、
を含み、
前記座屈条件は、
0006675129.tif 000034
で表される（Ｅはヤング率、ｒ０’は円柱の縦横比）ことを特徴とする、座屈条件導出方法。
【請求項２】
請求項１に記載の前記座屈条件を用いて円柱状繊維が破断する条件である破断条件を算出する方法であって、
前記座屈条件において、
等号が成り立つときのμｋの値をμｋｂｕとする工程と、
円柱状繊維が破断する際のμｋの値をμｋｂｒ、μｋｂｒのμｋｂｕに対する比を閾値ｒｂｒ＝μｋｂｒ／μｋｂｕとおく工程と、
閾値ｒｂｒをフィッティングパラメータとみなして実験結果に合うように決定する、又は破断モデルを組み込んだ繊維破断構造解析によって決定する工程と、
前記座屈条件においてμｋをμｋｂｒ／ｒｂｒ、不等号を等号で置き換えて破断する際のμｋｂｒを求める工程と、
を含み、
前記破断条件は、
0006675129.tif 000035
で表されることを特徴とする、破断条件算出方法。
【請求項３】
請求項２に記載の破断条件を用いて円柱状繊維の長さの分布である繊維長分布を予測する方法であって、
適宜設定した計測点に到達するまでのμｋの履歴の中での最大値μｋｍａｘを求める工程と、
前記最大値μｋｍａｘを前記破断条件のμｋｂｒに代入することによって得られる円柱状繊維の縦横比ｒ０’ｍａｘを算出する工程と、
を含み、
円柱状繊維の直径ｄｆを用いて、それぞれｄｆ／ｒ０’ｍａｘおよびｄｆ／（２ｒ０’ｍａｘ）として破断前後の繊維長を求めて、繊維長分布を予測することを特徴とする、繊維長分布予測方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、流体中を移動する繊維の座屈条件導出方法及び破断条件の算出方法並びに繊維長分布予測方法に関する。
【背景技術】
【０００２】
繊維で強化された繊維強化プラスチックは、近年、さまざまな分野で利用されている。繊維強化プラスチックの強度は、繊維の分散、配向、長さ等、繊維のミクロな状態に依存するため、成形過程を経て出てくるプラスチック中の繊維の微視的状態を把握、制御することが重要となる。このような繊維の微視的状態を把握するため、高粘度流体中での繊維の分散、配向状態をシミュレーションで予測する試みが広く行われるようになってきている（非特許文献１、２参照）。一方、繊維長を予測するためには、流れの中での繊維破断過程を解析することが必要になるが、繊維破断を定量的に評価することは難しく、シミュレーションで繊維長を予測するまでには至っていない。【０００３】
流れの中を移動する繊維状物体の配向および破断予測解析は、山本らにより試みられている（非特許文献３〜７参照）。山本らの解析では、繊維を複数の球粒子の結合体として表現し、流れの中で個々の球が受ける流体抵抗を合算することにより、繊維状物体が流れから受ける力を算出している。【０００４】
また、フェルプス（Ｐｈｅｌｐｓ）らは、収縮流中の細長い円柱に働く流体応力に対する近似解を用いることにより、収縮流の中で円柱状繊維が座屈する際の条件を導出した（非特許文献８、９参照）。さらに、導出された座屈条件式は、成形ソフトにも組み込まれている（非特許文献１０参照）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】岡田有司，中野亮：成形加工，27，134(2015)
【非特許文献２】正木炭典，名嘉山祥也，梶原稔尚:成形加工，27，533(2015)
【非特許文献３】Yamamoto, S. and T. Matsuoka: J. Chem. Phys., 98, 644(1993)
【非特許文献４】Yamamoto, S. and T. Matsuoka: Polymer Eng. Sci., 35, 1022(1995)
【非特許文献５】山本智:豊田中央研究所R&Dレビュー，33，63(1998)
【非特許文献６】山本智，松岡孝明:成形加工，11，510(1999)
【非特許文献７】山本智:日本レオロジー学会誌，29，185(2001)
【非特許文献８】Phelps, J. H.: Ph.D. Thesis, University of Illinois at Urbana-Champaign, IL 61801(2009)
【非特許文献９】Phelps, J. H., A.I.A. El-Rahman, V. Kunc and C.L. Tucker III: Composites: Part A, 51, 11(2013)
【非特許文献１０】Nguyen, B.N., X. Jin, J. Wang, J. H. Phelps, C. L. Tucker III, V. Kunc, S.K. Bapanapalli and M.T. Smith: FY 2010 Quarterly Report, PNNL-19185, Pacific Northwest National Laboratory(2010)
【発明の概要】
【発明が解決しようとする課題】
【０００６】
ところで、山本らの解析では、繊維を複数の球粒子の結合体として表現し、個々の球粒子に働く流体抵抗を合算することにより、繊維状物体が流れから受ける力を算出している。しかしながら、流れの中で物体が受ける力は、一般に、物体を模擬した結合粒子を構成する個々の球粒子に働く流体抵抗の和として表すことはできない。たとえば、一様流中におかれた球粒子の結合体で表された物体を考えると、レイノルズ数の低いストークス領域では、個々の球に働く流体抵抗は粒子径に比例するのに対し、物体を構成する球粒子の数は粒子径の３乗に反比例するため、結合粒子全体に働く流体抵抗は粒子径の２乗に反比例することになる。とくに粒子径が０の極限では、結合粒子全体に働く流体抵抗は無限大になる。このような結合粒子を用いる方法では、流れの中での繊維の配向予測は行えるものの、繊維が流れから受ける力を求めることは難しく、繊維破断を予測することまではできないという問題があった。【０００７】
また、フェルプスらの導出過程には、後述のように、理論的な側面でいくつかの疑問が残る。【０００８】
そこで本発明は、このような問題に鑑みてなされたものであり、その目的とするところは、フェルプスらが対象とした収縮流中の円柱周りの流れ場を、より厳密な方法により算出し、その上で、得られた円柱周りの流体応力分布から導かれる座屈に関する固有値問題を解くことによって、収縮流中におかれた円柱に対する座屈条件を導出する。また、座屈条件を利用した破断条件の算出方法、および繊維長分布を予測するための方法についても示す。
【課題を解決するための手段】
【０００９】
上記課題を解決するために、本発明の第１の観点によれば、流体中を移動する円柱状繊維が、繊維の長軸方向において収縮する方向への流れである収縮流中におかれたときに座屈する条件である座屈条件を導出する方法であって、円柱状繊維の縦横比ｒ０’に対して一意的に求まる無次元の流体応力分布に、円柱が存在している場所でのμｋ（μは流体粘度、ｋは収縮速度（繊維が存在する場所での長手方向の流速の勾配））を掛けることによって収縮流中の円柱表面に働く流体応力分布を得る工程と、
円柱表面での流体応力分布から導かれる座屈する閾値である最小固有値λ０を用いることにより、座屈条件を導出する工程と、
を含み、
前記座屈条件は、
0006675129.tif 000002
で表される（Ｅはヤング率、ｒ０’は円柱の縦横比）ことを特徴とする、座屈条件導出方法が提供される。【００１０】
上記課題を解決するために、本発明の第２の観点によれば、本発明の第１の観点にかかる前記座屈条件を用いて円柱状繊維が破断する条件である破断条件を算出する方法であって、前記座屈条件において、等号が成り立つときのμｋの値をμｋｂｕとする工程と、円柱状繊維が破断する際のμｋの値をμｋｂｒ、μｋｂｒのμｋｂｕに対する比を閾値ｒｂｒ＝μｋｂｒ／μｋｂｕとおく工程と、閾値ｒｂｒをフィッティングパラメータとみなして実験結果に合うように決定する、又は破断モデルを組み込んだ繊維破断構造解析によって決定する工程と、前記座屈条件においてμｋをμｋｂｒ／ｒｂｒ、不等号を等号で置き換えて破断する際のμｋｂｒを求める工程と、を含み、
前記破断条件は、
0006675129.tif 000003
で表されることを特徴とする、破断条件算出方法が提供される。【００１１】
上記課題を解決するために、本発明の第３の観点によれば、本発明の第２の観点にかかる前記破断条件を用いて円柱状繊維の長さの分布である繊維長分布を予測する方法であって、適宜設定した計測点に到達するまでのμｋの履歴の中での最大値μｋｍａｘを求める工程と、前記最大値μｋｍａｘを前記破断条件のμｋｂｒに代入することによって得られる円柱状繊維の縦横比ｒ０’ｍａｘを算出する工程と、を含み、円柱状繊維の直径ｄｆを用いて、それぞれｄｆ／ｒ０’ｍａｘおよびｄｆ／（２ｒ０’ｍａｘ）として破断前後の繊維長を求めて、繊維長分布を予測することを特徴とする、繊維長分布予測方法が提供される。
【発明の効果】
【００１２】
本発明によれば、フェルプスらが対象とした収縮流中の円柱周りの流れ場を、より厳密な方法により算出し、その上で、得られた円柱周りの流体応力分布から導かれる座屈に関する固有値問題を解くことによって、収縮流中におかれた円柱に対する座屈条件を導出することができる。これにより、流体中を移動する繊維の座屈条件導出方法及び破断条件の算出方法並びに繊維長分布予測方法を提供することができる。本発明のその他の効果については、以下の発明を実施するための形態の項でも説明する。
【図面の簡単な説明】
【００１３】
【図１】収縮流の中におかれたｒ０’＝１／５０の円柱周りの流速分布を示す図である。軸対称およびｚ’＝０面に対する面対称モデルを用い、ＣＦＤによって算出。左図が無次元の流速分布で、右図がｚ’方向速度成分コンター。円柱が存在しないときの収縮流（主流）を除いた主流からのずれの流速分布を表示する。
【図２】図１の流速分布から算出した円柱表面でのｒ’方向の無次元の速度勾配ｆ（ｚ’）＝δｖｚ’／δｒ’のグラフ（実線）である。破線は、式（１４）で表される速度勾配を無次元化したグラフを示す。
【図３】収縮流中の円柱に対する座屈解析で用いる記号の説明を示す図である。
【図４】数値計算によって求めた最小固有値λ０＝５．０９に対応する固有関数δ（ｚ’）の分布を示す図である。δの最大値が１、δ＝０（ｚ’＝±１）となるように規格化してある。
【図５】上図の流速分布の中を、２連結粒子がクリアランスを通過する様子を解析した事例を示す図である。下図中の粒子の濃淡は、２連結粒子のそれぞれの粒子位置での流速から求まる連結方向の速度勾配（流れの収縮・伸長速度ｋ）を表す。正の値が収縮流、負の値が伸長流に対応する。
【図６】フェルプスらが繊維長分布測定に用いた実験装置の概略を示す図である。
【図７】図６の実験装置の円板部分を対象とした流動解析結果を示す図であり、上段が流速分布で、下段が粘度分布を示す。
【図８】本繊維破断予測理論より求めた図６の装置内のＡ点での繊維長分布（実線）とフェルプスらの実験データ（棒線）との比較を示す図である。
【発明を実施するための形態】
【００１４】
以下に添付図面を参照しながら、本発明の好適な実施の形態について詳細に説明する。なお、本明細書及び図面において、実質的に同一の機能構成を有する構成要素については、同一の符号を付することにより重複説明を省略する。
【００１５】
（１）本発明の一実施形態について、図１〜図８を参照しながら説明する。以下の説明において、上記非特許文献１〜１０のほか、以下の非特許文献１１〜２０を参照する。
【非特許文献１１】巽友正：流体力学，培風館(1982)
【非特許文献１２】吉澤微：流体力学，東京大学出版会(2001)
【非特許文献１３】Forgacs, O.L. and S.G. Mason: J. Colloid Interface Sci., 14, 457 (1959)
【非特許文献１４】Forgacs, O.L. and S.G. Mason: J. Colloid Interface Sci., 14, 473 (1959)
【非特許文献１５】Salinas, A. and J.F.T. Pittman: Polymer Eng. Sci. 21, 23 (1981)
【非特許文献１６】von Turkovich R. and L. Erwin: Polymer Eng. Sci. 23, 743 (1983)
【非特許文献１７】中村喜代治：非ニュートン流体力学，コロナ社(1997)
【非特許文献１８】Dinh, S.M. and R.C. Armstrong: J. Rheology, 28, 207(1984)
【非特許文献１９】Jeffery, G.B.:Proc. R. Soc. London, Ser. A, 102, 161(1922)
【非特許文献２０】Advani, S.G. and C.L. Tucker III: J. Rheology, 31, 751(1987)
【００１６】
（２）対象とする流れ場
高粘度流体中を移動する繊維の周りの流れに対しては、一般的に、ストークス近似が成り立つ。また、繊維近傍では、粘度は近似的に一定とみなすことができる。このような条件下では、繊維周りの流れの支配方程式は線形になり、下記のように記述される（非特許文献１１、１２参照）。
【００１７】
0006675129.tif 000004
【００１８】
ここに、ｖは流速ベクトル、ｐは流体圧力、μは粘度である。なお、繊維とともに移動する座標系で考えるものとし、式（１）では時間微分項についても無視した。ここで、流れ場を、細長い繊維状物体の周りで、繊維中心ｘ０を中心に、以下の（式２）のようにテイラー展開する。ここに、下付き添え字ｉ（ｉ＝１，２，３）、ｊ（ｊ＝１，２，３）はベクトルの（ｘ、ｙ、ｚ）成分を表す。流れの支配方程式（１）が線形であることから、流れ場の重ね合わせを行うことができ、流れ場をテイラー展開の各項ごとに分離して考えることが可能となる。
【００１９】
0006675129.tif 000005
【００２０】
まず、式（２）の右辺第１項は繊維の移動速度ｖ０で移動する一様流を表しており、繊維とともに移動する座標系で考えると無視できる。次に、右辺第２項に注目すると、第２項に含まれる９つの速度勾配テンソル成分の内、せん断変形を表す項は、ｘ０を中心とする回転流を表し、繊維とともに回る回転系で考えると、無視できる。残りの速度勾配テンソル成分は、細長い繊維の長軸方向の伸長もしくは収縮流に対応しており、この内の収縮流が、以前より指摘されているように、繊維破断を引き起こす要因と考えられる（非特許文献１３〜１６参照）。そこで、本実施形態でも、繊維の周りの流れ場から一様流と回転流成分を取り除いた収縮流を考え、その中におかれた繊維を模擬した細長い円柱に働く流体応力について考察する。
【００２１】
収縮（伸長）流は、一般に下記の流速分布によって表される（非特許文献１７参照）。
【００２２】
0006675129.tif 000006
【００２３】
ここに、（ｖｘ、ｖｙ、ｖｚ）は流速の（ｘ、ｙ、ｚ）成分で、流れはｚ＝±∞からｚ＝０に向かって収縮するものとする。また、ｋはｚ方向の速度勾配で、収縮速度（ｋ＜０のときは伸長速度）に対応している。なお、収縮（ｚ）方向に垂直な方向の速度成分（ｖｘ、ｖｙ）はパラメータｍ（０≦ｍ≦１）に依存し、ｍ＝０のときは軸対称流、ｍ＝１のときは平面流れとなる。しかしながら、座屈を引き起こすのは円柱表面に働くｚ方向の流体応力で、（ｖｘ、ｖｙ）の影響は２次的と考えられるため、ここでは、ｍの違いによる（ｖｘ、ｖｙ）の違いは無視する。こうしたことから、次節「（３）収縮流中の円柱に働く流体応力分布」では、式（３）においてｍ＝０としたときの軸対称収縮流（以下の（式４））の中におかれたｚ軸を中心軸とする細長い円柱の周りの流れを対象とする。なお、フェルプスらが行った定式化でも、同様の軸対称収縮流が仮定されている。
【００２４】
0006675129.tif 000007
【００２５】
（３）収縮流中の円柱に働く流体応力分布
式（４）で表される軸対称収縮流の中におかれた半径ｒ０、長さ２ｌの円柱の周りの流れを考える。ただし、円柱の中心軸はｚ軸で、中央がｚ＝０とする。流れの支配方程式は式（１）、円柱表面および無限遠での境界条件は
【００２６】
0006675129.tif 000008
【００２７】
となる。なお、以下の式展開は、収縮流（ｋ＞０）、伸長流（ｋ＜０）のいずれの場合にも適用されるが、以降、繊維破断の直接の要因とされる収縮流（ｋ＞０）を前提に説明する。式（１）および境界条件（５）を以下により無次元化する。
【００２８】
0006675129.tif 000009
【００２９】
すると、以下の式（６）となる。
【００３０】
0006675129.tif 000010
【００３１】
式（６）に含まれるパラメータとしては、境界条件も含めて円柱の縦横比ｒ０’＝ｒ０／ｌのみとなるため、式（６）を解いて得られる無次元の流れ場（流速と圧力）は、ｒ０’のみに依存することになる。
【００３２】
（３−１）収縮流中の円柱周りの流れに対する数値解
式（６）は、ＣＦＤ（Computational Fluid Dynamics：数値流体力学）によって解くことができる。図１に、本件出願人が開発している商用の流体・粉体解析ソフトＲ−ＦＬＯＷを用いて、式（６）を軸対称およびｚ’＝０面に対する面対称条件の下で解いた結果（無次元の流速分布とｚ’方向速度成分分布）を示す。ただし、ｒ０’＝１／５０の場合となっている。
【００３３】
なお、式（６）がｒ０’のみに依存することから、このようにして求まった無次元の流速分布（ｖｒ’’、ｖｚ’’）あるいは（ｖｒ’、ｖｚ’）は、ｒ０’が変わらない限り、不変となる。図１のように円柱周りの流速と圧力分布が得られると、円柱表面での流体応力分布の算出に必要な円柱表面での速度勾配と円柱両端での圧力分布が求まる。
【００３４】
図２は、円柱表面でのｒ’方向の無次元の速度勾配ｆ（ｚ’）＝δｖｚ’／δｒ’をｚ’の関数としてプロットしたグラフである。関数ｆ（ｚ’）はｒ０’が変わらない限り不変な関数となるため、同じｒ０’に対しては、ｆ（ｚ’）は一度求めておくだけよい。同様に、円柱の両端面にかかる無次元の流体圧力ｐ’および粘性応力を求めるための速度勾配についても、ｒ’の関数として求まるため、面全体での中心に向かう無次元の流体応力が、無次元の圧力と速度勾配をｒ’≦ｒ０’の面内で積分することにより算出される。なお、無次元の圧力と速度勾配に対する積分は数値的に行う必要があるが、積分値についてもｒ０’のみに依存することになるため、同じｒ０’に対しては、積分値も一度求めておくだけでよい。
【００３５】
実際の流れの中で円柱側面にかかる粘性応力μδｖｚ／δｒは、図２中の無次元の速度勾配に収縮速度ｋを掛けることによって有次元に変換した後、粘度を掛けて、以下の式（７）によって得られる。
【００３６】
0006675129.tif 000011
【００３７】
同様に、円柱の両端面にかかる流体応力についても、無次元の圧力と速度勾配にμｋを掛けることにより算出することができる。
【００３８】
以上のことから、収縮流中の円柱表面に働く流体応力分布は、与えられた円柱の縦横比ｒ０’に対して一意的に求まる無次元の流体応力分布に、円柱が存在している場所でのμｋを掛けることによって得られることになる。なお、図２に示したような円柱表面上の無次元の速度勾配、あるいは両端での圧力を求める作業は、与えられたｒ０’の値に対しては一度行えばよい。
【００３９】
（３−２）両端近傍を除く円柱周りの流れに対する理論解
円柱が十分に細長い、すなわち縦横比が十分に小さい場合（ｒ０’＝ｒ０／ｌ≪１）、両端近傍を除くと、円柱周りの流れは両端の影響を受けなくなる。両端の影響が無視できる場合の円柱周りの軸対称流を求めるにあたって、支配方程式（１）をストークスの流れの関数φ（ｒ、ｚ）を用いて書き直すと、次のようになる。
【００４０】
0006675129.tif 000012
【００４１】
式（８）を解くにあたり、
【００４２】
0006675129.tif 000013
【００４３】
とおくと、ｆ（ｒ）に関して
【００４４】
0006675129.tif 000014
【００４５】
が得られる。上式を積分するとｆ（ｒ）に対して次式が得られる。
【００４６】
0006675129.tif 000015
【００４７】
ここに、ｃ１、ｃ２、ｃ３、ｃ４は積分定数である。速度成分（ｖｒ、ｖｚ）は、式（８ｂ）より、
【００４８】
0006675129.tif 000016
【００４９】
により得られる（ｒ’＝ｒ／ｌ）。積分定数ｃ１、ｃ２、ｃ３、ｃ４は（ｖｒ、ｖｚ）に対する境界条件（８ｃ）から決定されるが、ストークス近似の下では、無限遠で対数関数的に増大する解が含まれることが知られている（ストークスのパラドックス）（非特許文献１１、１２参照）。そのため、無限遠での境界条件からはｃ１のみを０とし、ｃ２については円柱表面の境界条件から決定するものとする。こうして、積分定数が境界条件から
【００５０】
0006675129.tif 000017
【００５１】
として決まり、流速分布が下記のように求まる。
【００５２】
0006675129.tif 000018
【００５３】
式（１３）で表される流速の内、右辺第１項の（ｋｒ／２、−ｋｚ）は円柱が存在しない場合の収縮流、すなわち主流に対応し、右辺第２項以下が主流からのずれを表している。式（１３ｂ）より、円柱表面での速度勾配は、以下の式（１４）となる。
【００５４】
0006675129.tif 000019
【００５５】
図２に、ｒ０’＝１／５０の場合の式（１４）で表される速度勾配を、無次元化した上で、破線でプロットしてある。図をみると、円柱の先端付近を除くと、速度勾配を式（１４）で近似しても、ずれはそれ程大きくはなく、繊維形状を円柱で表すこと自体が近似であることを考えると、繊維形状がある程度細長い場合には、式（１４）で近似してもよいと考えられる。
【００５６】
（４）収縮流中の円柱に働く流体応力に基づく座屈・破断解析
本節では、前節で求めた収縮流中の円柱表面に働く流体応力分布を用いることにより、収縮流中の円柱に対する座屈条件を導く。また、座屈・破断解析を行う方法についても述べる。
【００５７】
（４−１）座屈解析
円柱の両端に中心に向かう力が作用する場合の座屈条件は、オイラーによって初めて解析され、オイラー座屈として知られている。フェルプスら（非特許文献１８参照）は、ディン（Ｄｉｎｈ）ら（非特許文献１６参照）が近似的に求めた収縮流中の円柱表面に働く流体応力を用いて、オイラー座屈の結果を適用することにより、収縮流中の円柱に対する座屈条件を導出している。しかしながら、ディンらは、円柱表面での速度勾配が円柱表面から遠く離れた場所での（円柱が存在しないときの）収縮流速に比例すると仮定し、その際の比例係数を“抵抗係数”としているが、抵抗係数についての理論的導出は行っておらず、不確定なパラメータとして残っている。そのため、フェルプスらがディンらの近似解に基づいて導出した座屈条件にも、この抵抗係数がそのまま残っており、結果的にフィッティングパラメータになってしまっている。なお、「（３−２）両端近傍を除く円柱周りの流れに対する理論解」で、円柱の両端の影響が無視できる場合の流速分布の導出を行ったが、実際にはこれはディンらの論文における抵抗係数の理論的導出に相当しており、ディンらの論文における“抵抗係数”は、式（１４）の速度勾配と粘度の積として表される。したがって、ディンらの近似解の代わりに「（３−２）両端近傍を除く円柱周りの流れに対する理論解」で得られた解を用いることにより、不確定なパラメータを排除することができる。
【００５８】
フェルプスらが行った座屈条件の導出には、上記以外にもう一つ次のような問題がある。フェルプスらは、オイラー座屈問題をそのままの形で適用するために、ディンらの近似理論から求めた円柱中央部に働く流れによる圧縮力を、オイラー座屈問題での円柱の両端に働く内向きの力に置き換えることにより、座屈条件を導出しているが、そのような置き換えが許されることの根拠は示されていない。このようなことから、フェルプスらの結果の妥当性には疑問が残る。
【００５９】
そこで、本節では、ディンらの近似解の代わりに前節で得られた収縮流中の円柱表面での流体応力分布を用いて、収縮流中の円柱に対する座屈条件を導出する。その際、フェルプスらが行ったようなオイラー座屈解析の結果をそのまま転用するという方法ではなく、前節で求まった円柱表面での流体応力分布そのものを用いることによって導かれる固有値問題を解くことにより、座屈条件を導出する。
【００６０】
収縮流中におかれた円柱の座屈条件を求めるにあたり、図３に示すように、円柱が座屈する前の中心軸（ｚ方向）に沿って長さｄｚの微小要素を取り出し、軸に垂直方向（ｘ方向）の力のつり合いを考えると、次式が得られる。
【００６１】
0006675129.tif 000020
【００６２】
ここに、Ｎ、ｄＮ、Ｔ、ｄＴ、Ｆ、θ１、θ２、θ、ｄｓの各変数の定義は、図３に記載されている通りで、Ｆは円柱表面に働くｚ方向の単位長さあたりの収縮流による力に相当する。ここで、ｘ方向の変位をδとし、微小変形を仮定すると、
【００６３】
0006675129.tif 000021
【００６４】
と近似されるため、式（１５）は
【００６５】
0006675129.tif 000022
【００６６】
となる。式（１６）の両辺をｄｚで割り、さらに、軸(ｚ)方向の力の釣り合いからｄＮ≒−Ｆｄｓ≒−Ｆｄｚが成り立つことを用いると、式（１６）は以下の式（１７）となる。
【００６７】
0006675129.tif 000023
【００６８】
ここで、曲げモーメントＭ、ヤング率Ｅ、断面２次モーメントＩ（＝πｒ０４／４）を導入の上、以下の式（１８）の関係を用いる。
【００６９】
0006675129.tif 000024
【００７０】
すると、式（１７）は以下の式（１９）となる。
【００７１】
0006675129.tif 000025
【００７２】
なお、式（１９）のＮは、円柱内部の軸方向の力の釣り合いから、円柱の端面で受ける流体圧力と粘性応力の総和Ｐと側面で受ける粘性応力Ｆの関数として次のように表される。
【００７３】
0006675129.tif 000026
【００７４】
式（２０）におけるＰとｖｚは、「（３）収縮流中の円柱に働く流体応力分布」で求めた流れ場を用いて算出される。
【００７５】
円柱が十分に細長く（ｒ０’≪１）、円柱表面での速度勾配が式（１４）で近似できる場合には、円柱両端における圧力および粘性応力の影響は無視でき（Ｐ＝０）、式（２０）は次式（２１）で表される。
【００７６】
0006675129.tif 000027
【００７７】
式（２１）を式（１９）に代入すると、式（２２）が得られる。
【００７８】
0006675129.tif 000028
【００７９】
ここでさらに、独立変数をｚからｚ’（＝ｚ／ｌ）に変換すると、式（２２）は以下の式（２３）となる。
【００８０】
0006675129.tif 000029
【００８１】
式（２３）に対する境界条件としては、円柱両端において自由端の条件、すなわち、曲げモーメントおよび円柱両端でのｘ方向のせん断力が０という条件が課せられる。この境界条件の内、曲げモーメントが０という条件は、式（１８ａ）から、以下の式（２４）となる。
【００８２】
0006675129.tif 000030
【００８３】
一方、ｘ方向のせん断力は、一般的にはＴとＰの線形結合として表されるが、円柱が十分に細長い場合には、Ｐ＝０とおいてよいため、Ｔに等しくなる。そのため、ｘ方向のせん断力が０という境界条件はＴ＝０と表され、式（１８ｂ）を用いると、以下の式（２５）となる。
【００８４】
0006675129.tif 000031
【００８５】
微分方程式（２３）は、境界条件（２４）、（２５）とともに、λを固有値とする固有値問題となっており、固有値と固有関数δ（ｚ’）は、式（２３）、（２４）、（２５）を数値的に解くことにより得られる。固有値の内、最小固有値をλ０とするとき、実際に数値計算によって求めた値は、λ０＝５．０９となっている。
【００８６】
図４に、対応する固有関数δ（ｚ’）の分布を示す。ただし、δの最大値が１、δ＝０（ｚ’＝±１）となるように規格化してある。図４の固有関数分布を見ると、円柱中央で変位δが最も大きくなることから、最小固有値を超えるような流体応力が働いたときには、円柱中央部で２つに折れ曲がるように座屈することがわかる。
【００８７】
オイラーの座屈条件に相当する座屈条件は、式（２３ｂ）および最小固有値λ０を用いて、λ≧λ０とおくことにより、以下の式（２６）あるいは（２７）として得られる。
【００８８】
0006675129.tif 000032
【００８９】
フェルプスらがディンらの近似解を利用して得た座屈条件式とは異なり、式（２７）には、抵抗係数のような不確実なパラメータは含まれない。
【００９０】
（４−２）繊維配向解析を利用した繊維座屈予測
前節で得られた収縮流中の円柱に対する座屈条件を、繊維配向解析の結果に適用することにより、流れの中のどの位置で円柱状繊維が座屈するかを予測できる。
【００９１】
繊維配向解析の手法には、楕円あるいは楕円球を用いる方法（非特許文献１９、２０参照）と結合粒子を用いる方法（非特許文献１、３〜７参照）の２通りの方法があり、いずれの方法を用いてもよい。
【００９２】
図５は、２結合粒子を用いて、繊維状粒子がどのようにクリアランスの中を通過するかを解析した事例である。クリアランスの中で、繊維状粒子は流れ方向に揃う様子が解析されている。なお、図５中の粒子の濃淡は、２結合粒子のそれぞれの粒子位置での流速から求まる連結方向の速度勾配、すなわち流れの収縮・伸長速度ｋを表しており、正の値が収縮流、負の値が伸長流に対応している。こうして繊維配向解析により求まった収縮速度ｋと粘度μに対して、繊維が十分に細長いとみなせる場合の座屈条件（２７）を適用することにより、円柱状繊維がどの位置で座屈するかを予測することができる。
【００９３】
（４−３）破断解析
「（４−２）繊維配向解析を利用した繊維座屈予測」では、流れの中での繊維の座屈予測を行う方法を示したが、座屈した繊維が必ず破断するわけではなく、座屈は破断の必要条件でしかない。そのため、繊維長分布を予測するためには、繊維破断そのものを予測するあるいは破断条件と座屈条件の関連付けを行う必要がある。繊維が十分に細長い（ｒ０’≪１）とみなせる場合を考え、座屈条件（２７）において、等号が成り立つときのμｋの値をμｋｂｕとすると、流体中の繊維が存在している場所でのμｋがμｋｂｕを越えた時点で繊維は座屈し始めることになる。その後、繊維がより大きなμｋの場所に移動すると破断にいたると考えられるが、その破断する際のμｋの値をμｋｂｒ、μｋｂｒのμｋｂｕに対する比をｒｂｒ＝μｋｂｒ／μｋｂｕとおく。なお、収縮流の中で繊維表面に働く流体応力μｋに比例するため、ｒｂｒ（＞１）は繊維に働く流体応力が座屈開始時の何倍になると破断するかを表している。以下、破断に関する閾値μｋｂｒあるいはｒｂｒをどのようにして求めるかについて言及する。
【００９４】
破断条件あるいは閾値の算出方法としてまず考えられるのが、円柱状繊維を対象とした破断モデルを組み込んだ構造解析を行うことによって、破断する際の閾値を求める方法である。ただし、実際に破断構造解析を実施する際には、繊維破断モデルが必要になる。
【００９５】
破断条件あるいは閾値を求める上でのもう一つ考えられる方法は、ｒｂｒを定数で近似し、ｒｂｒをフィッティングパラメータとみなして実験結果に合うように決定する方法である。この方法では、繊維長分布に関する実験データが必要になるが、実際には、ｒｂｒは同じ繊維に対しては同じ値を取ると考えられるため、ある特定の実験データに対してｒｂｒをフィッティングにより決定してしまえば、同じ繊維に対しては共通に用いることができる。閾値ｒｂｒが求まると、座屈条件（２７）においてμｋをμｋｂｒ／ｒｂｒ、不等号を等号で置き換えた次式によって、μｋｂｒを求めることができる。
【００９６】
0006675129.tif 000033
【００９７】
破断予測解析の例として、フェルプスらの繊維長分布測定実験（非特許文献８、９参照）を追随することにより、座屈条件を利用した本繊維破断予測理論によってどの程度繊維長分布を予測できるかを検証する。フェルプスらは、図６の装置中央の突起部の上端から直径１８０ｍｍ、厚さ３ｍｍの円板領域内に、直径１７μｍの円柱状繊維を含有した液を流し込んだ上で、図中の円板内のＡ、Ｂ、Ｃの３か所（中心からそれぞれ１５ｍｍ、４５ｍｍ、７５ｍｍの位置）での繊維長分布を測定している。
【００９８】
なお、フェルプスらは、繊維長の測定に加えて、フェルプスらが提案している繊維破断予測理論を用いて得られた繊維長分布と実測データとの比較についても行っている。ただ、フェルプスらが行った解析では、装置中央の突起部を無視して、円板部のみを対象としている他、繊維長の流入条件として、装置内に流入時の繊維長分布ではなく、Ａ点での繊維長分布実測データを流入条件として与えて、Ｂ点、Ｃ点での繊維長分布を解析により予測するということを行っている。
【００９９】
しかしながら、図６の装置では、円板内の外側にいくほど流速は遅くなることから、座屈を引き起こす収縮速度は円板中央部で最大となり、外側ほど弱くなる。そのため、座屈の結果として生じる繊維破断のほとんどは装置中央付近で起きていると考えられる。実際、フェルプスらが得た繊維長分布の実測データでは、Ａ点、Ｂ点、Ｃ点でほぼ同様の繊維長分布になっている。このようなことから、フェルプスらが行ったＡ点での繊維長分布実測データを入力条件として与えて、Ｂ点、Ｃ点での繊維長分布を予測するという方法は、検証としてどれほどの意味を持つか疑問である。
【０１００】
そこで、本解析では、より上流の地点での繊維長分布を与え、Ａ点での繊維長分布を予測することを行う。ただ、フェルプスらの文献には、装置中央の突起物に関する寸法の記載がないため、本解析でもフェルプスらが行った解析と同様、装置中央部の突起物を省略し、円板内部のみを対象とした解析を実施する。なお、流れ場の解析では軸対称性を仮定し、流体は円板中央部上面の半径５ｍｍ（具体的な寸法がわからないため適当に仮定）の領域から一定の速度で下向きに流入すると仮定した。流入流量、粘度モデル等、その他の条件、物性については、フェルプスらの文献値をそのまま用いた。
【０１０１】
一方、繊維の投入については、流入面の下部領域に図５と同様の繊維を模擬した約５０００本の結合粒子を発生し、その際の方向は乱数によって与えた。なお、フェルプスらの文献には、装置内流入時の繊維長分布データが掲載されていないため、本解析では、投入時の繊維長は適当に６ｍｍと仮定した。ただし、フェルプスらが得たＡ点での実測データでは、大半の繊維は繊維長が６ｍｍよりも短いことから、６ｍｍ程度の長い繊維は、Ａ点に到達する前にほとんどが破断し、より短い繊維となっていると考えられる。したがって、流入時の繊維長については、ある程度長く設定しておけば、解析結果にはあまり依存しないものと考えられる。
【０１０２】
このように、フェルプスらの繊維長分布実測結果を追随するにあたっては不明な点があり、厳密な再現はできないものの、検証に必要な詳細な繊維長分布測定データが掲載されている文献は見当たらないため、本解析ではフェルプスらの測定データを用いる。このようなことから、本解析では、本繊維破断予測理論の定量的精度の検証ではなく、本繊維破断モデルによって予測される繊維長分布が、概ね実験結果を再現し得るのか、あるいは全くできないのかを検証することを目的とする。
【０１０３】
図７に、ＣＦＤ解析で得られた円板中央付近での流速および粘度分布を示す。図８には、本繊維破断予測理論に基づいて得られたＡ点における繊維長分布を、フェルプスらの実測データとともにプロットしてある。ただし、ｒｂｒとしては多少の試行錯誤の結果、ｒｂｒ＝２を用いた。
【０１０４】
図８の繊維長の具体的な計算方法としては、Ａ点を通過した各結合粒子に対して、Ａ点に到達するまでのμｋの履歴の中での最大値μｋｍａｘを求めておく。さらに、このμｋｍａｘを式（２８）中のμｋｂｒに代入することによって得られるｒ０’（円柱の縦横比）の値、ｒ０’ｍａｘを算出する。なお、繊維が座屈の結果破断する場合、図４から、繊維の中央で２分割される可能性が最も高いと考えると、破断前後の繊維長は、繊維の直径ｄｆ（＝１７μｍ）を用いて、それぞれｄｆ／ｒ０’ｍａｘおよびｄｆ／（２ｒ０’ｍａｘ）として求まる。ただし、こうして求まった破断前の繊維長ｄｆ／ｒ０’ｍａｘが装置内投入時の繊維長よりも大きい場合には、繊維は破断することなく、投入時の繊維長が維持される。
【０１０５】
なお、本解析では、投入時の繊維長を比較的長めの６ｍｍと仮定したため、ほとんどの繊維はＡ点を通過するまでに破断すると思われるが、実際の実験では、繊維が装置内に流入した時点で繊維長に分布があり、短い繊維も存在していると考えられるため、一部の繊維については、破断することなく、投入時の繊維長を維持したままＡ点を通過していることも考えられる。
【０１０６】
図８をみると、繊維長分布は実測と本解析で違いが見られるものの、平均値を含めた分布は概ね捉えられていることがわかる。先に述べたように、本解析では、繊維投入時の繊維長分布として装置投入時の実測データではなく一定値を仮定しており、投入部の装置形状および投入位置も実装置とは異なっている。このような条件下であっても、繊維長分布が概ね再現できていることを考えると、本繊維破断モデルは、フェルプスらの実験装置内で起きている繊維破断現象の本質は捉えているものと考えられる。
【０１０７】
（５）結言
収縮流の中におかれた細長い円柱状繊維表面に働く流体応力分布をＣＦＤ等によって簡便に求める方法を提案した。また、円柱の両端の影響が無視できる場合の理論解を導出し、得られた流体応力分布に基づいて、固有値問題を解くことにより、座屈条件を導出した。また、座屈条件を繊維配向解析の結果に適用することにより、流れの中で繊維が座屈する位置を予測する方法についても示した。さらに、座屈条件を利用した破断条件を算出することにより、流れの中での繊維破断を予測するための方法についても提案した。また、 フェルプスらの繊維長分布測定実験結果を追随することにより、本繊維破断予測理論が繊維長分布を予測する上で有効であることを示した。
【０１０８】
（本実施形態の効果）
以上説明したように、本実施形態によれば、フェルプスらが対象とした収縮流中の円柱周りの流れ場を、より厳密な方法により算出し、その上で、得られた円柱周りの流体応力分布から導かれる座屈に関する固有値問題を解くことによって、収縮流中におかれた円柱に対する座屈条件を導出することができる。これにより、流体中を移動する繊維の座屈条件導出方法及び破断条件の算出方法並びに繊維長分布予測方法を提供することができる。
【０１０９】
以上、添付図面を参照しながら本発明の好適な実施形態について説明したが、本発明はかかる例に限定されないことは言うまでもない。当業者であれば、特許請求の範囲に記載された範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
【図１】
0006675129.tif 000036
【図２】
0006675129.tif 000037
【図３】
0006675129.tif 000038
【図４】
0006675129.tif 000039
【図５】
0006675129.tif 000040
【図６】
0006675129.tif 000041
【図７】
0006675129.tif 000042
【図８】
0006675129.tif 000043