JP2000322407A - 流体中における媒体の運動解析方法 - Google Patents

Abstract

(57)【要約】 【課題】流体中の媒体の挙動を効率良く解析すること。 【解決手段】流動性のある基質に含まれる媒体の運動解
析方法であって、媒体１０を密球体モデルの構成個数Ｎ
より少ないｎ個の球体１１で構成する。即ち、ｎ個の球
体１１が間隙を有して所定の自由度で結合された集合体
とする。それに伴い、運動方程式の並進力および回転力
に個数Ｎ，ｎに基づく補正係数を乗算する。これにより
密球体モデルと同等の結果を得るとともに計算量を大幅
に低減させる。よって、流体中の媒体の移動、変形、破
断、配向等が効率良く解析できる。

Description

【発明の詳細な説明】

【０００１】

【産業上の利用分野】本発明は、流体中に混入された微
小媒体の運動解析手法に関する。特に、流体中に混入さ
れた微小媒体を間隙を有した球体でモデル化し、その移
動、変形、破断、配向等の運動状態を効率よく解析する
方法に関する。本発明は、樹脂成形において樹脂に混入
されたファイバ等の挙動を解析する解析手法に適用でき
る。

【０００２】

【従来の技術】従来より、流体中の媒体の運動を解析す
る手法として、特開平５−３１４０９１に開示された流
動基質中の媒体の運動解析方法がある。それは、流体中
に混入された媒体形状に対して複数の球体を１次元、２
次元又は３次元方向に結合させてモデル化し、各球体に
複数の自由度を持たせて、各球体に運動方程式を立てて
全体の動きを解析する手法である。ここで、流動性のあ
る基質とは、その内部に媒体を含み得る流動可能な分散
媒質であって、液体、気体、半流動体の他、生産時また
は使用時においてー時的に流動性を持つもの、例えば成
形時における溶融状態のプラスチックや金属等である。
以下、簡単のため流動性媒質を単に流体と記す。また、
媒体とは例えばプラスチックを強化するガラス繊維等で
ある。

【０００３】例えば、図１１（ａ）に示すように流体中
の媒体がファイバ１００である時には、１次元的に配列
された球体集合体モデル１１０が与えられる。その球体
集合体を構成する各球体１１１には隣り合う球体との間
に媒体の引張り弾性に相当する結合長の自由度、及び媒
体の曲げ弾性に相当する結合角の自由度が設定されてい
る。また、媒体が平板状物体２００であれば、各球体１
１１が２次元的に配列された球体集合体モデル２１０が
与えられる。

【０００４】上記結合長の自由度とは、モデル化の対象
となった媒体の引張り弾性に相当するものであり、隣り
合う球体間の距離が伸びたとき、これを元に戻そうとす
る力として設定される。また、結合角の自由度とは、モ
デル化の対象となった媒体の曲げ弾性に相当するもので
あり、隣り合う球体間の結合角が変化したとき、これを
元に戻そうとするトルクとして設定される。

【０００５】上記結合長の自由度と結合角の自由度は、
実際の媒体の引張り弾性と曲げ弾牲に対応して任意に与
えることができる。一般的に、剛直な媒体の場合ほどこ
れらの自由度が小さく設定され、柔軟な媒体の場合ほど
これらの自由度が大きく設定される。上記手法は、それ
らの設定に基づいて各球体の並進、回転を計算し、その
計算結果を総計して全体の移動、変形、配向等を解析す
る解析方法である。

【０００６】実際には、さらに解析の基礎情報として、
基質の粘度、流れ場の分布等の基質情報および基質の流
れ場における初期位置、初期速度等の初期情報を与え
て、流体中で媒体に引き起こされる伸びや曲げ等の変形
を経時的に解析している。

【０００７】

【発明が解決しようとする課題】しかしながら、この方
法は媒体を例えばＮ個の密接して連続した球体の集合体
でモデル化し、ある時刻における各球体の位置や姿勢に
基づいて、各球体に作用する総合力、総合トルクを求
め、それらの力と並進運動、回転運動に関する運動方程
式から次の時刻における各球体の位置や姿勢を求めると
いう方法を逐次繰り返すことで、媒体の時間経過に伴う
運動状態を求める方法である。例えば、直径Ｄ、長さＬ
の繊維状媒体は、直径Ｄの球体をＮ個並べてモデル化
し、各球体について上記位置、姿勢等を求める。従っ
て、その計算に必要な時間と記憶容量はＮ² に比例して
増大する。さらに、流体中には、多数の媒体が存在する
ため、その個々の媒体に関する動作解析を行うには膨大
な計算量を必要し、現実的には不可能である。

【０００８】従って、本発明の目的は、流体中に混入さ
れた混入媒体を間隙を有した球体でモデル化し、媒体の
移動、変形、破断、配向、又は衝突等の運動状態を計算
量を大幅に低減させて容易に解析することである。

【０００９】

【課題を解決するための手段】上記課題を解決するため
に、請求項１の流体中における媒体の運動解析方法は、
流体中に含まれる媒体をＮ個の球体が自由度をもって密
接して結合された密球体モデルで表現し、ある時刻にお
ける各球体に作用する総合力、総合トルクを求め、この
総合力、総合トルクから各球体の並進運動、回転運動に
関する運動方程式を立て、これらの運動方程式を数値解
析して各球体の位置及び回転角を計算し、その各球体の
計算結果を総計して、その媒体の移動、変形、破断、配
向、又は衝突等の運動状態を逐次解析する方法であっ
て、その媒体をＮ＞ｎであるｎ個の球体が間隙を有して
所定の条件で結合された集合体でモデル化し、その運動
方程式の係数に個数Ｎ，ｎに基づく補正を与えて密球体
モデルと同等の結果を得るとともにその計算量を低減さ
せたことを特徴とする。

【００１０】また、請求項２の流体中における媒体の運
動解析方法は、その補正が運動方程式において流体から
受ける並進力、回転力およびせん断力に対して行われる
ことを特徴とする。

【００１１】また、請求項３の流体中における媒体の運
動解析方法は、その補正が運動方程式において媒体の復
元力に関する弾性係数に対して行われることを特徴とす
る。

【００１２】また、請求項４の流体中における媒体の運
動解析方法は、その補正が運動方程式において媒体の相
互作用による外力に対して行われることを特徴とする。

【００１３】また、請求項５の流体中における媒体の運
動解析方法によれば、補正は媒体を各球体間に仮想直線
を結んだセグメント集合体で表現し、衝突が予想される
セグメントの箇所に仮想球体を設けて、その仮想球体を
用いて計算した仮想球体間の衝突による並進力に対して
行い、並進力をそのセグメントに属する各球体に所定の
比率で分配することを特徴とする。

【００１４】また、請求項６の流体中における媒体の運
動解析方法によれば、前記補正は、前記各球体を前記個
数Ｎ，ｎに基づいて隣接球が相互に接触するように拡大
された仮想球体で表現し、該仮想球体間の相互作用によ
って発生する全ての外力に対して行われることを特徴と
する。全ての外力とは、例えば並進力、回転力、せん断
力、摩擦力、弾性力、衝突等を含む。ここで、完全に接
触した仮想球は、摩擦力と摩擦力に基づく外力および弾
性力を演算する場合に適する。完全に接触せず相互の離
間距離が直径以下である仮想球は、衝突および衝突に基
づく外力を演算する場合に適する。

【００１５】

【作用及び発明の効果】請求項１の流体中における媒体
の運動解析方法によれば、その媒体をＮ＞ｎであるｎ個
の球体が間隙を有して所定の条件で結合された集合体で
モデル化している。流体中における媒体をＮ個の球体の
集合で表現し、その各球体に対して運動方程式をたて
て、各球体の位置、回転角等を計算すると、その計算量
はＮ² に比例する。本請求項では、所定の条件のもと
で、上記媒体をＮ＞ｎであるｎ個の球体で間隙を有して
表現している。よって、その計算量を大幅に低減させる
ことができる。ここに、所定の条件とは、各球体間には
例えば球体の直径以上の間隙が存在しても、互いにある
媒体は他の媒体の間隙を通過できないという条件であ
る。また、球体間に摩擦を考慮する時は、上記構成個数
Ｎ、ｎに基づいて、その半径をＮ−１／ｎ−１倍して摩
擦を計算することである。さらに、例えば流体からの外
力を計算する場合には、半径を拡大しないという条件で
ある。

【００１６】また、例えば媒体に作用する従来の外力を
使用して、Ｎ個の球体をｎ個で置き換えると、球体１個
当たりに作用する外力が大きくなる。これにより、密球
体モデルでの解析解と誤差が生じる。本請求項では、個
数Ｎ，ｎに基づく数値実験によりその誤差がなくなる様
に運動方程式の各外力および係数に補正を与えている。
よって、密球体モデルと同等の結果を得ることができ
る。これにより、計算量を低減させて密球体モデルと同
等の結果を得る効率のよい解析方法となる。

【００１７】また、請求項２の流体中における媒体の運
動解析方法によれば、その補正が運動方程式において流
体から受ける並進力、回転力およびせん断力に対して行
われる。上記並進力、回転力およびせん断力は流体の基
質から上記各球体に作用する全ての外力である。その全
ての外力について、上記個数Ｎ、ｎに基づく補正を行っ
ているので、正確にその作用を各球体の位置、角度に反
映させることができる。これにより、正確に流体の作用
をその媒体の位置、配向に反映させて解析することがで
きる。

【００１８】また、請求項３の流体中における媒体の運
動解析方法によれば、その補正が運動方程式において媒
体の復元力に関連する弾性係数に対して行われる。上記
弾性係数は、その媒体の伸縮に関する弾性係数と、その
媒体の曲げ変形に関する曲げ弾性係数、および捻れに関
する全ての弾性係数を含む。全ての弾性係数について、
上記個数Ｎ、ｎに基づく補正を行っているので、正確に
媒体の弾性による作用を各球体の位置、角度に反映させ
て解析することができる。

【００１９】また、請求項４の流体中における媒体の運
動解析方法によれば、その補正は運動方程式において媒
体の相互作用によって発生する外力に対して行われる。
ここに、相互作用による外力とは、媒体と媒体の衝突・
摩擦等を意味する。流体中には多数の媒体が混入され、
移動過程においてそれらの媒体間には衝突、摩擦などの
相互作用がある。この相互作用によっても、その媒体の
位置、形状、配向が決定される。本請求項の運動解析方
法は、その相互作用に対して上記個数Ｎ、ｎに基づく補
正を行っているので、正確に相互作用を各球体の位置、
角度に反映させて解析することができる。これにより、
媒体間の相互作用による媒質の位置、形状、配向、変形
を正確に求めることができる。

【００２０】また、請求項５の流体中における媒体の運
動解析方法は、媒体を各球体間を仮想直線で結んだセグ
メント集合体で表現し、仮想直線の最も接近した２点を
求め、その位置に仮想球体を設けて、その仮想球体間の
衝突を並進力とする。Ｎ個の密球体モデルで表現した媒
体をＮ＞ｎであるｎ個の球体で表現した場合、各球体間
には間隙が生じる。媒体と媒体の衝突を計算する場合に
は、計算量は低減されるものの、この間隙を他の媒体が
通り抜ける事象が生じる。これは、現実の衝突を表すも
のではない。

【００２１】本請求項の運動解析方法では、上記のよう
に衝突時にはセグメントで表現し、そのセグメント上に
仮想球体を設けて衝突を計算している。よって、密球体
モデルでのシミュレーションと同様、媒体と媒体の衝突
を表現することができる。また、仮想球体で得られた上
記外力を仮想球体の位置に応じた比率で、両側の各球体
に分配している。これにより、密球体モデルでのシミュ
レーションと同等の衝突による媒体の位置、形状、配向
の結果を得ることができる。

【００２２】また、請求項６の流体中における媒体の運
動解析方法によれば、各球体に個数Ｎ，ｎに基づいて拡
大され、相互に接触された仮想球体を想定し、その仮想
球体間の摩擦（相互作用）によって発生する外力を各球
体に働く力としている。

【００２３】Ｎ個の密球体モデルで表現した媒体をＮ＞
ｎであるｎ個の球体で表現した場合、各球体間には間隙
が生じる。そして、媒体と媒体の接触を計算する場合に
は、その計算量は低減されるものの、この間隙があるた
め接触のない事象が生ずる。これは、密球体モデルと一
致するものではない。本請求項では、個数Ｎ，ｎに基づ
いて各球体を拡大して仮想球体とし、その仮想球体で作
製した仮想の密球体モデルを用いている。よって、確実
に媒体と媒体の摩擦力を表現し、その摩擦力から並進力
と回転力を得ることができる。

【００２４】

【発明の実施の形態】以下、本発明の実施の形態につい
て図面を参照して説明する。尚、本発明は下記実施例に
限定されるものではない。図１に、本実施例の流体中に
おける媒体１０を示す。媒体１０は、簡単化するため長
さＬ、直径２ａの１次元繊維とする。図１（ａ）に繊維
形状を、図１（ｂ）にＮ個の密接して連続した球体で示
した密球体モデル（以降、従来モデル）を、図１（ｃ）
に本実施例のモデルを示す。尚、媒体が平板状であれば
球体１１が２次元的に配列された集合体で、媒体が立体
であれば球体を３次元に配列した集合体でモデル化す
る。

【００２５】本実施例の特徴は、従来、Ｎ個の球体１１
が密に結合された媒体１０を、Ｎ＞ｎであるｎ個の球体
１１で表現したことである。また、各球体１１間に６自
由度を有した仮想的なバネ１２を想定し、そのバネ１２
によって相互に結合させたことである。これにより、媒
体１０の位置、配向等に対する計算量を大幅に低減する
とともに媒体１０の変形を正確に算出することができ
る。本実施例では、上記ｎ個の球体モデルを粗視化モデ
ルと記し、この粗視化モデルでの媒体１０の変形の解析
方法について説明する。尚、上記６自由度とは、例えば
球体１１の位置及び姿勢を表すｒ（＝（ｘ，ｙ，ｚ），
座標による位置）、オイラー角Θ（＝（θ，φ，ψ））
等の各パラメータである。各球体は、この６つの方向に
所定の条件で変位可能とする。

【００２６】このような条件の媒体１０（粗視化モデ
ル）を流動場中に想定し、各球体１１の回転、及び並進
に対する運動方程式を次式で表現する。尚、下記の式は
全てベクトルである。

【数１】球体の回転運動に関して、次式が成立する。 Ｉ・ｄ² Θ_i ／ｄｔ² ＝ Ａ・Ｔ_i ^h ＋Ｂ・ΣＴ_i ^b ＋ΣＴ_i ^t −Σａ・ｆ_ij×ｎ_ij …（１） 球体間ですべりのない条件式として次式が成立する。 ｄｒ_i ／ｄｔ＋Ｃ・ａ・ｄΘ_i ／ｄｔ×ｎ_ij ＝ｄｒ_j ／ｄｔ＋Ｃ・ａ・ｄΘ_j ／ｄｔ×ｎ_ji …（２） ここに、Ａ，Ｂ，Ｃ：補正係数 Θ_i ：球体の姿勢に関するオイラー角のベクトル量で、
（α，φ，ψ）である。 Ｉ：球体の慣性モーメント Ｔ_i ^h ：流体から受けるトルク Ｔ_i ^b ：バネから受けるトルク Ｔ_i ^t ：バネから受けるねじりトルク ａ ：球体の半径 ｆ_ij：球体ｉが球体ｊから受ける摩擦力 ｎ_ij：球体ｉからｊへの単位ベクトル ｉ，ｊ：球体を示すインデックス 尚、Ｔ_i ^h は詳細には、Ｔ_i ^h ＝−８πη₀ ａ³ （ｄΘ
_i ／ｄｔ−ω（ｒ_i ））である。ω（ｒ_i ）は位置ｒ_i
における流動場の角速度である。

【００２７】又、並進力に対する運動方程式として次式
が成立する。

【数２】 ｍ・ｄ² ｒ_i ／ｄｔ² ＝ Ｐ・Ｆ_ij ^p ＋Ｑ・Ｆ_i ^h ＋ΣＦ_i ^s ＋Σｆ_ij …（３） ここに、 Ｐ，Ｑ，：補正係数 ｒ_i ：位置（ｘ，ｙ，ｚ） ｍ：質量 Ｆ_i ^h ：流体から受ける力 Ｆ_i ^s ：バネから受ける力 Ｆ_ij ^p ：球体ｉの球体ｊからの相互作用 ｆ_ij：球体ｉの球体ｊからの摩擦力 ｉ，ｊ：球体を示すインデックス である。

【００２８】また、Ｆ_i ^h は詳細には、Ｆ_i ^h ＝−６π
η₀ ａ（ｄｒ_i ／ｄｔ−ｖ（ｒ_i ））である。ｖ
（ｒ_i ）は、位置ｒ_i における流動場の速度である。

【００２９】上記の（１）、（３）を（２）式の条件の
基に解く。（１）、（２）、（３）式は、ベクトル量に
関する式で、実際には、各成分毎に式を立てて、解くこ
とになる。

【００３０】上記補正係数Ａ，Ｂ，Ｃは、粗視化モデル
のための補正係数である。これは、本来Ｎ個で表現すべ
き集合体を、ｎ個、即ち間隙を有して表現した場合、流
体から受けるトルクＴ_i ^h およびバネから受けるトルク
Ｔ_i ^b に誤差が生じるからである。そして、その誤差は
その間隙距離に係わると推定されるからである。ここで
は、この係数Ａ，Ｂが乗算される第１項および第２項に
ついて説明する。他の項の作用は、従来モデルと同等で
ある。また、Ｃは後述する仮想球体による半径の拡大係
数であり、上記トルクに関係する。

【００３１】先ず第１項の補正係数ＡおよびＣについて
説明する。例えば、各球体１１間が間隙を有したままで
あると、両球体１１間の摩擦力ｆ_{i j} が０になる（図
２）。従って、図２に示すように球体１１の外周に互い
に接触する直径２λａの仮想球１３を考える。またこの
時、終端条件として、図３に示す様に必ず媒体の終端に
球体１１の中心と、仮想球体１３の中心が位置する様に
設定する。両端の球体間の中心距離Ｌ’は、従来モデル
では２ａ（Ｎ−１）であり、粗視化モデルでは２λａ
（ｎ−１）である。よって、λ＝Ｎ−１／ｎ−１とな
る。即ち、仮想球体１３の直径は従来モデルの構成個数
Ｎと粗視化モデルの構成個数ｎで関連づけられる。ま
た、この２λａの仮想球体１３を想定すれば、従来モデ
ルとその運動方程式およびその関係式が使用できる。ま
た、これにより上記（２）式の係数Ｃはλに設定され
る。

【００３２】この時、仮想球体１３の単位面積当たりの
外力をそのままに設定しておくと、表面積が増大してい
るので、受ける外力が大きくなると予想される。従っ
て、外力の係数を補正して従来モデルと粗視化モデルの
結果を数値実験で一致させる必要がある。

【００３３】数値実験とは、せん断速度場で従来モデル
の媒体１０と粗視化モデルの媒体１０を回転させる実験
である。その回転周期を一致させて、λを変化させた時
の従来モデルと粗視化モデルでの外力比を求めるもので
ある。この外力比を補正係数αとする。その数値実験の
結果を図４のグラフに示す。曲線ａとｂが、周期が一致
した時の係数λと補正係数αの関係である。示される曲
線の数値解析により、次式による補正係数が得られた。

【数３】 α（λ）＝ａ₁ ・ａｔａｎ（ｂ₁ λ）＋ｃ₁ （λ＜λ_th） …（４）

【数４】 α（λ）＝ａ₂ ・ａｔａｎ（ｂ₂ λ）＋ｃ₂ （λ＞λ_th） …（５）

【００３４】ここに、ａ₁ ，ａ₂ ，ｂ₁ ，ｂ₂ ，ｃ₁ ，
ｃ₂ は定数であり、周期運動について、数値実験より理
論解との比較から決定する。また、λ_thは曲線ａとｂの
交点であり、λ＜λ_th時には（４）式を使用し、λ＞λ
_th時には（５）式を使用する。また、このαに基づいた
両モデルでの数値実験によるトルクの比較より、式
（１）の第１項の流体から受けるトルクＴ_i ^h の補正係
数Ａはα³ であることを得た。

【００３５】次に第２項の補正係数Ｂについて説明す
る。第２項は、バネから受けるトルクの項である。これ
は、従来モデルと粗視化モデルの弾性係数の比較によっ
て説明できる（図５）。図５（ａ）が従来モデルであ
り、図５（ｂ）が粗視化モデルである。ともに外力によ
る変形量は同じであるとする。この時の変形角をそれぞ
れθ₁ 、θ₂ とすれば、図５よりθ₁ ＝２ａ／Ｒ、θ₂
＝２λａ／Ｒである。粗視化モデルでは１個当たりの変
形角が通常モデルのλ倍であるので、１個当たりのバネ
定数は通常モデルの１／λ倍となる。換言すれば、粗視
化モデルでは、モーメントを１／λ倍の個数のバネで受
けるので、その全体のバネ定数は通常モデルの１／λ倍
である。以上から、粗視化モデルのバネ定数ｋ_RMは従来
モデルの１／λ² にする必要がある。よって、各球体が
上式（１）の第２項のバネから受けるトルクＴ_i ^b の補
正係数Ｂは１／λ² であることを得る。

【００３６】図６にこの補正係数Ｂを導入した片持ち梁
の解析結果を示す。片持ち梁は、一端が固定された縦弾
性係数Ｅ＝１０ＭＰａ、直径２ａ＝０．５ｍｍ，長さＬ
＝５．０ｍｍの繊維からなる。図６は、他の一端に荷重
Ｆ＝１．０ｍＮを与えた結果である。横軸が繊維長さｘ
ｍｍ（座標）であり、縦軸がたわみ角θである。折れ線
ｅが球体の構成個数を４とした場合である。従来モデル
の場合、球体個数は１０であるが、個数を４としても理
論曲線ｆに近似できることが示された。

【００３７】このように、流体中の媒体にＮ＞ｎである
ｎ個の球体からなる粗視化モデルを採用し、回転の運動
方程式にＮ、ｎに基づく補正係数を導入すれば、計算量
を大幅に低減させて、正確に媒体の変形を予測すること
ができる。

【００３８】尚、上記係数αを大幅に増大させて、せん
断解析を行うこともできる。せん断流動場で、最大トル
クを受ける様に媒体を配置し回転させる。せん断流動場
がある大きさを越えると媒体は破断する。破断の決定
は、次式のＭａｓｏｎのせん断速度場における繊維の座
屈条件を用いた。

【００３９】

【数５】 η₀ ・γ／Ｅ＝（ｌｎ（２ｄ）−１．７５）／２ｄ⁴ …（６） ここに、η₀ ：係数 １／γ：回転周期 Ｅ：係数 ｄ：繊維の直径 である。尚、破断ルールは、１．破断後も媒体長は変化
しない。２．破断後も媒体の形状を正しく表現する。
３．破断後も１媒体を最低球数の２個で表現できること
とした。

【００４０】図７に媒体のアスペクト比（Ｌ／２ａ）に
対する限界せん断応力曲線を示す。実線ｎがＭａｓｏｎ
の理論値である。プロット点が流体力補正である補正係
数Ａ（＝α³ ）＝３６９％時の限界せん断応力および形
状補正である補正係数Ｂ（１／λ² ）＝５０％の限界せ
ん断応力である。形状補正とは、曲げ変形に対する補正
である。流体力補正係数α³ ＝３６９％とすると、Ｍａ
ｓｏｎの理論式とよく一致する。これにより、粗視化モ
デルでもせん断解析を精度よく行うことができる。

【００４１】また、図８にせん断速度場に置かれた媒体
の周期運動ついて、上記補正を用いて数値実験を行った
結果を示す。図８（ａ）に、せん断速度場と媒体１０の
配置を示す。尚、媒体１０は２球体で表現したものであ
る。図８（ｂ）にアスペクト比２０の、図８（ｃ）にア
スペクト比１００の媒体１０の終端軌跡を示す。両者と
もＪeffreyの理論解にほぼ近似されている。また、図９
に補正の効果をアスペクト比と周期の関係図で示す。横
軸がアスペクト比であり、縦軸が周期である。補正を行
った場合をｐ₁ でプロットし、行わない場合をｐ₂ でプ
ロットした。補正を行った場合は、Ｊefferyの理論解に
一致することが示された。これにより、媒体をＮ＞ｎで
あるｎ個の球体を間隙を有して構成しても、上記補正を
行うことで、ほぼ従来モデルと同等の解析が可能となっ
た。よって大幅に計算時間を低減する優れた解析方法と
なる。

【００４２】以上のようにして、従来モデルに代えてｎ
個の球体モデルからなる粗視化モデルを採用し、回転に
対する運動方程式に係数を乗算することで媒体１０の変
形を精度良く解析することができる。次に、並進運動に
関する方程式（３）を解くことにより、同じく粗視化モ
デルを用いて、媒体の位置を正確に解析する方法につい
て説明する。

【００４３】図１（ｃ）に示す粗視化モデルの媒体１０
を流れ場に設定すると、上記（３）式のうよに並進力に
対する運動方程式が成立する。

【００４４】上記相互作用Ｆ_ij ^p を式（７）〜（８）に
示す。相互作用Ｆ_ij ^p は、接近する２球体の中心の離間
距離｜ｒ_j −ｒ_i ｜の関数であり、媒体同士が非常に接
近した場合に発生し、先ず流体力学的相互作用（引力）
が発生し、さらに接近すると剛体反発力が働く。

【数６】 Ｆ_ij ^p =０ …（７）

【数７】 Ｆ_ij ^p =3/2a² πη(v_j -v_i )n_ij/(｜r _j -r_i ｜-2a)n _ij …（８）

【数８】 Ｆ_ij ^p =-D₀ ・exp(G₀(1- ｜r _j -r_i ｜/ 2a))n _ij …（９） 但シ、（７）式は｜r _j -r_i ｜>3a が条件である。（８）
式は2.0001a<｜r _j -r_i ｜≦3aが条件である。（９）式
は、｜r _j -r_i ｜≦ 2.0001a が条件である。ここに、
ｖ_i は球体ｉの速度、n _ijはｒ_ijベクトルの単位ベクト
ルであり、D₀，G₀各球体の反発力に関する係数である。
これらの相互作用Ｆ_ijp は、粒子登録法により時間が短
縮されて計算されている。粒子登録法とは、全ての球体
に対しては計算せず、媒体を構成するある特定の粒子を
登録し、その粒子に関してのみ上記相互作用を計算する
ものである。この登録粒子は、適当な時間後他の粒子に
登録される。

【００４５】上記（３）式の補正係数Ｐ，Ｑは、粗視化
モデルのための補正係数である。これは、本来Ｎ個で表
現すべき集合体を、ｎ個、即ち間隙を有して表現したか
らである。即ち、流体から受ける並進力Ｆ_i ^h 、球体間
の相互作用から受ける反発力Ｆ_ij ^pに誤差を生じるから
である。そして、その誤差はその間隙距離に係わると推
定されるからである。ここでは、上記補正係数Ｐ，Ｑが
乗算される第１項および第２項について説明する。他の
項は、従来モデルと同等である。

【００４６】先ず第１項の補正係数Ｐについて説明す
る。第１項は衝突による相互作用である。例えば、各球
体１１間が間隙を有したままであると、現実にはあり得
ない両媒体が互いに通り抜ける事象がある。従って、図
１０（ａ）に示すように互いに隣接した球体１１間を結
ぶ仮想直線を考えセグメント１５、１６とする。その
際、相互の直線距離が所定値以内であるセグメントが選
択され、登録される。そして、登録されたセグメントに
おいて、両セグメントの間隔が最小となる２点ｋ₁ 、ｋ
₂ を設定し、その上に仮想球体１７、１８を想定する。

【００４７】この仮想球体１７、１８の衝突による相互
作用Ｆ_ij ^p は上式（７）〜（９）に従って求められる。
そして、その相互作用Ｆ_ij ^p は、セグメント１５の点ｋ
₁ 上に設けた仮想球体１７の位置による配分比ｓ：（１
−ｓ）に応じて、両側の各球体１１ａ、１１ｂに分配さ
れる（図１０（ｂ））。図１０（ｂ）は式（９）で表せ
られる反発力である。これにより、粗視化モデルにおい
ても従来モデルと同等の衝突が計算される。

【００４８】次に、第２項の補正係数Ｑについて説明す
る。第２項は、流体から受ける並進力である。仮想球体
を想定すると流体から受ける並進力が大きくなる。従っ
て、その流体による外力を補正する必要がある。これ
は、（４）、（５）式の従来モデルと粗視化モデルとの
外力比αが使用できる。よって、補正係数Ｑはαとな
る。これにより、流体からの並進力が補正され、それを
用いた運動方程式から従来モデルと同等の解が得られ
る。

【００４９】このように、Ｎ＞ｎであるｎ個の球体から
なる粗視化モデルを採用して、並進に関する運動方程式
の係数Ｐ，Ｑを補正すれば、計算量を大幅に低減させる
とともに、正確に媒体の位置を予測することができる。

【００５０】（変形例）以上、本発明を表す１実施例を
示したが、他に様々な変形例が考えられる。例えば、係
数αは、関数ａ・ｔａｎ（ｂ₁ λ）から求めたが、他の
関数を近似解としてそれから求めてもよい。例えば、双
曲線関数でもよい。

【００５１】又、１つの１次元線状物体について、近似
する球体の数は任意である。又、両端に１個づつの球を
設けることで、最も簡単で、計算量を最も減少させるこ
とができる。上記の（１）〜（３）式を同時に逐次的に
数値解析することで、媒体の並進、回転、曲げ、せん断
状態を同時に求めることができるが、ある状態情報だけ
必要とするならば、それに関する方程式のみを解けば良
い。上記の数値解析は、コンピュータにより実施される
が、本発明は、媒体の運動状態を解析する方法の他、上
記のように補正された係数を有する方程式を解く各ステ
ップを有するプログラムを記憶媒体に記憶させて提供す
ることもできる。又、そのプログラムを実行するＣＰＵ
を有したコンピュータシステムにより解析装置として提
供することも可能である。また、各球体の位置、角度か
ら媒体の位置、配向、変形を数値で求めたが、それを視
覚化してもよい。また、上記例では各媒体間の相互作用
を求めたが、流路壁との衝突、摩擦に応用し、流路内で
の媒体の偏在等の解析にも使用できる。さらに上記実施
例では、互いに等しい直径の球体、あるいは仮想球体を
想定したがそれらの大きさを異ならせることもできる。

【図面の簡単な説明】

【図１】本発明の具体的な一実施例に係る流体中に置か
れる媒体を間隙を有した球体の集合体でモデル化した模
式図。

【図２】実施例に係る間隙を有した球体間に接触する仮
想球体の設定説明図。

【図３】実施例に係る仮想球体の直径を求める説明図。

【図４】実施例にかかる仮想球体の係数λと補正係数α
の関係図。

【図５】実施例に係る曲げ弾性発生時の球体および仮想
球体の配置図。

【図６】実施例に係る補正係数Ｂを導入した片持ち梁の
解析結果説明図。

【図７】実施例に係る媒体のアスペクト比と限界せん断
応力の関係図。

【図８】実施例に係るせん断速度場における媒体の軌跡
を用いた比較図。

【図９】実施例に係るせん断速度場における媒体のアス
ペクト比と周期の関係図。

【図１０】実施例に係るセグメントによる衝突説明図。

【図１１】従来の流体中における媒体を連続した球体で
モデル化した模式図。

【符号の説明】

１０ 媒体 １１ 球体 １１ａ、１１ｂ 球体 １２ バネ １３ 仮想球体 １５，１６ セグメント １７，１８ 仮想球体

───────────────────────────────────────────────────── フロントページの続き (72)発明者 松岡 孝明 愛知県愛知郡長久手町大字長湫字横道41番 地の１ 株式会社豊田中央研究所内 Ｆターム(参考） 5B056 AA06 BB03 BB83 HH00

Claims (6)

【特許請求の範囲】
1. 【請求項１】流動性のある基質に含まれる媒体の運動解
   析方法であって、該媒体をＮ個の球体が自由度をもって
   密接して結合された密球体モデルで表現し、ある時刻に
   おける各球体に作用する総合力、総合トルクを求め、こ
   の総合力、総合トルクから各球体の並進運動、回転運動
   に関する運動方程式を立て、これらの運動方程式を数値
   解析して各球体の位置及び回転角を計算し、その計算結
   果を総計して、前記媒体の移動、変形、破断、配向、又
   は衝突等の運動状態を逐次解析する方法を簡略化した媒
   体の運動解析方法において、 前記媒体をＮ＞ｎであるｎ個の球体が間隙を有して所定
   の条件で結合された集合体でモデル化し、前記運動方程
   式の並進力および回転力に前記個数Ｎ，ｎに基づく補正
   を与えて前記密球体モデルと同等の結果を得るとともに
   前記計算の計算量を低減させることを特徴とする流体中
   における媒体の運動解析方法。
2. 【請求項２】前記補正は、前記運動方程式において流体
   から受ける並進力、回転力およびせん断力に対して行わ
   れることを特徴とする請求項１に記載の流体中における
   媒体の運動解析方法。
3. 【請求項３】前記補正は、前記運動方程式において前記
   媒体の復元力に関連する弾性係数に対して行われること
   を特徴とする請求項１又は請求項２に記載の流体中にお
   ける媒体の運動解析方法。
4. 【請求項４】前記補正は、前記運動方程式において前記
   媒体の相互作用による外力に対して行われることを特徴
   とする請求項１乃至請求項３のいずれか１項に記載の流
   体中における媒体の運動解析方法。
5. 【請求項５】前記補正は、前記媒体を前記各球体間に仮
   想直線を結んだセグメント集合体で表現し、衝突が予想
   されるセグメントの箇所に仮想球体を設けて、該仮想球
   体を用いて計算した仮想球体間の衝突による並進力に対
   して行い、該並進力を前記セグメントに属する各球体に
   所定の比率で分配することを特徴とする請求項４に記載
   の流体中における媒体の運動解析方法。
6. 【請求項６】前記補正は、前記各球体を前記個数Ｎ，ｎ
   に基づいて隣接球が相互に接触するように拡大された仮
   想球体で表現し、該仮想球体間の相互作用によって発生
   する全ての外力に対して行われることを特徴とする請求
   項１乃至請求項４のいずれか１項に記載の流体中におけ
   る媒体の運動解析方法。