(57)【要約】【課題】複数の球の集合体としてモデル化した構造体の流体中の挙動を、並進及び回転の２つの運動方程式を解くことで解析する場合と比較して、計算時間を低減させることができる。【解決手段】簡易運動解析装置１０は、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するモデル化部１２と、球の回転角速度を算出する回転角速度算出部１８と、回転角速度に基づいて流体から球に働く粘性トルクを算出する粘性トルク算出部２０と、粘性トルクを等価な第１の力に変換する粘性トルク変換部２２と、球の相対位置関係に基づいて曲げトルクを算出する曲げトルク算出部２４と、曲げトルクを等価な第２の力に変換する曲げトルク変換部２４と、第１の力、第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、構造体の挙動を解析する解析部３２と、を含む。【選択図】図１
2018088085.tif 000002
【特許請求の範囲】
【請求項１】
流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するステップと、
前記球の回転角速度を算出するステップと、
前記回転角速度に基づいて、流体から前記球に働く粘性トルクを算出するステップと、
前記粘性トルクを等価な第１の力に変換するステップと、
前記球の相対位置関係に基づいて曲げトルクを算出するステップと、
前記曲げトルクを等価な第２の力に変換するステップと、
前記第１の力、前記第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、前記構造体の挙動を解析するステップと、
を含む流体中における構造体の簡易運動解析方法。
【請求項２】
前記球の回転角速度を算出するステップは、連結された２個の球の相対位置および前記球の速度に基づいて算出する
請求項１記載の流体中における構造体の簡易運動解析方法。
【請求項３】
前記球の回転角速度を算出するステップは、連結された２個の球の間に設定した仮想球の回転角速度を算出し、
前記粘性トルクを算出するステップは、前記仮想球の位置における渦度ベクトルを算出し、算出した前記渦度ベクトル及び前記仮想球の回転角速度に基づいて、前記粘性トルクを算出する
請求項２記載の流体中における構造体の簡易運動解析方法。
【請求項４】
流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するモデル化部と、
前記球の回転角速度を算出する回転角速度算出部と、
前記回転角速度に基づいて流体から前記球に働く粘性トルクを算出する粘性トルク算出部と、
前記粘性トルクを等価な第１の力に変換する粘性トルク変換部と、
前記球の相対位置関係に基づいて曲げトルクを算出する曲げトルク算出部と、
前記曲げトルクを等価な第２の力に変換する曲げトルク変換部と、
前記第１の力、前記第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、前記構造体の挙動を解析する解析部と、
を含む流体中における構造体の簡易運動解析装置。
【請求項５】
コンピュータに、請求項１〜３の何れか１項に記載の流体中における構造体の簡易運動解析方法の各ステップを実行させるための流体中における構造体の簡易運動解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、流体中における構造体の簡易運動解析方法、装置、及びプログラムに関する。
【背景技術】
【０００２】
従来、流体中の媒体の運動解析方法として、特許文献１、２、及び非特許文献１に開示された方法がある。これらの方法は、図３に示すように、任意形状の構造体５０を球５２の集合体としてモデル化する方法であり、各球５２について並進および回転の運動方程式を数値解析によって解くことにより、構造体５０の運動をシミュレートする。
【０００３】
並進の運動方程式は下記（１）式で与えられ、回転の運動方程式は下記（２）式で与えられる。
【０００４】
【数１】
2018088085.tif 000003
・・・（１）
【０００５】
【数２】
2018088085.tif 000004
・・・（２）
【０００６】
ここで、ｍは球の質量、ｒｉは球ｉの位置ベクトル、ＦＳｉｊは引っ張り復元力、ｆｉｊは摩擦力、Ｆｈｉは粘性力、ＦＰｉｋは球同士の流体力学的相互作用力、ＦＲｉｌは壁面反発力、ａは球の半径、θｉは球ｉの配向ベクトル、Ｔｂｉｊは曲げトルク、Ｔｔｉｊはねじりトルク、ｎｉｊは球ｉから球ｊへ向かう単位方向ベクトル、Ｔｈｉは粘性トルクである。
【０００７】
また、上記（１）、（２）式に加えて、結合した球同士がすべることなく運動するための拘束条件（すべり無し条件）を下記（３）式によって与える。
【０００８】
【数３】
2018088085.tif 000005
・・・（３）
【０００９】
ここで、ｖｘは球ｘ（ｘはｉまたはｊ）の速度ベクトル、ωｘは球ｘ（ｘはｉまたはｊ）の角速度ベクトル、ｎｊｉは球ｊから球ｉへ向かう単位方向ベクトルである。
【００１０】
また、上記（１）式のＦＳｉｊ、Ｆｈｉは下記（４）、（５）式により算出される。
【００１１】
【数４】
2018088085.tif 000006
・・・（４）
【００１２】
【数５】
2018088085.tif 000007
・・・（５）
【００１３】
ここで、Ｅは弾性係数、ｒは図４（ａ）に示すように、球ｉの中心と球ｊの中心を結ぶ線の長さである結合長、ｒ０は初期結合長、ηは流体の粘度、Ｖ（ｒｉ）は位置ベクトルｒｉで表される位置における流体の速度ベクトルである。
【００１４】
また、上記（２）式の、Ｔｂｉｊ、Ｔｔｉｊ、Ｔｈｉは下記（６）〜（８）式により算出される。
【００１５】
【数６】
2018088085.tif 000008
・・・（６）
ここで、ｎ⊥ｉｊは、ｎｉｊに垂直な単位ベクトルである。
【００１６】
【数７】
2018088085.tif 000009
・・・（７）
ここで、Ｇはせん断弾性係数である。
【００１７】
【数８】
2018088085.tif 000010
・・・（８）
【００１８】
ここで、θｂは図４（ｂ）に示すように、球ｉと球ｊとの結合角の変位、θｔは図４（ｃ）に示すように、球ｉと球ｊとのねじり角の変位、Ω（ｒｉ）は位置ベクトルｒｉで表される位置における流体の渦度ベクトルである。
【先行技術文献】
【特許文献】
【００１９】
【特許文献１】特開２０００−３２２４０７号公報
【特許文献２】特許平１１−２５０６９号公報
【非特許文献１】S. Yamamoto, T. Matsuoka, “A method for dynamic simulation of rigid and flexible fibers in a flow field”, J. Chem. Phys. 98, 644 (1993)
【発明の概要】
【発明が解決しようとする課題】
【００２０】
しかしながら、上記従来技術では、並進の運動方程式である上記（１）式、回転の運動方程式である上記（２）式の両方を解く必要があり、計算時間が増大する、という問題があった。
【００２１】
本発明は上記事情に鑑み成されたものであり、本発明の目的は、複数の球の集合体としてモデル化した構造体の流体中の挙動を、並進及び回転の２つの運動方程式を解くことで解析する場合と比較して、計算時間を低減させることができる流体中における構造体の簡易運動解析方法、装置、及びプログラムを提供することにある。
【課題を解決するための手段】
【００２２】
上記目的を達成するために、請求項１に係る発明の流体中における構造体の簡易運動解析方法は、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するステップと、前記球の回転角速度を算出するステップと、前記回転角速度に基づいて、流体から前記球に働く粘性トルクを算出するステップと、前記粘性トルクを等価な第１の力に変換するステップと、前記球の相対位置関係に基づいて曲げトルクを算出するステップと、前記曲げトルクを等価な第２の力に変換するステップと、前記第１の力、前記第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、前記構造体の挙動を解析するステップと、を含む。
【００２３】
なお、請求項２に記載したように、前記球の回転角速度を算出するステップは、連結された２個の球の相対位置および前記球の速度に基づいて算出するようにしてもよい。
【００２４】
また、請求項３に記載したように、前記球の回転角速度を算出するステップは、連結された２個の球の間に設定した仮想球の回転角速度を算出し、前記粘性トルクを算出するステップは、前記仮想球の位置における渦度ベクトルを算出し、算出した前記渦度ベクトル及び前記仮想球の回転角速度に基づいて、前記粘性トルクを算出するようにしてもよい。
【００２５】
請求項４記載の発明の流体中における構造体の簡易運動解析装置は、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するモデル化部と、前記球の回転角速度を算出する回転角速度算出部と、前記回転角速度に基づいて流体から前記球に働く粘性トルクを算出する粘性トルク算出部と、前記粘性トルクを等価な第１の力に変換する粘性トルク変換部と、前記球の相対位置関係に基づいて曲げトルクを算出する曲げトルク算出部と、前記曲げトルクを等価な第２の力に変換する曲げトルク変換部と、前記第１の力、前記第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、前記構造体の挙動を解析する解析部と、を含む。
【００２６】
請求項５記載の発明の流体中における構造体の簡易運動解析プログラムは、コンピュータに、請求項１〜３の何れか１項に記載の流体中における構造体の簡易運動解析方法の各ステップを実行させるための流体中における構造体の簡易運動解析プログラムである。
【発明の効果】
【００２７】
本発明によれば、複数の球の集合体としてモデル化した構造体の流体中の挙動を、並進及び回転の２つの運動方程式を解くことで解析する場合と比較して、計算時間を低減させることができる。
【図面の簡単な説明】
【００２８】
【図１】簡易運動解析装置の構成を示すブロック図である。
【図２】簡易運動解析処理ルーチンのフローチャートである。
【図３】構造体のモデル化について説明するための図である。
【図４】並進及び運動の方程式のパラメータについて説明するための図である。
【図５】回転角速度の算出について説明するための図である。
【図６】粘性トルクの変換について説明するための図である。
【図７】曲げトルクの算出について説明するための図である。
【図８】曲げトルクの算出について説明するための図である。
【図９】曲げトルクの変換について説明するための図である。
【図１０】実施例に係る構造体について説明するための図である。
【図１１】実施例に係る繊維の回転周期に関する、従来技術と本発明との誤差について説明するための図である。
【図１２】実施例に係る従来技術と本発明の計算時間の高速化割合について説明するための図である。
【発明を実施するための形態】
【００２９】
以下、図面を参照して本発明の実施形態を詳細に説明する。
【００３０】
＜本発明の概要＞
【００３１】
まず、本発明の概要について説明する。
【００３２】
本実施形態では、前述したように、図３に示すように任意形状の構造体５０を球５２の集合体としてモデル化する。そして、各球５２について、前述した回転の運動方程式である上記（２）式は解かずに、下記（９）式の並進の運動方程式のみを解くことによって構造体５０の挙動を解析する。
【００３３】
【数９】
2018088085.tif 000011
・・・（９）
【００３４】
上記（９）式の右辺第５項のＦｂｉｊは、前述した回転の運動方程式である上記（２）式の曲げトルクＴｂｉｊを等価な力（第２の力）に変換したものである。また、上記（９）式の右辺第６項のＦh_torqueｉｊは、前述した回転の運動方程式である上記（２）式の粘性トルクＴｈｉを等価な力（第１の力）に変換したものである。なお、本実施形態では、第３の力として、球の引っ張り復元力ＦＳｉｊ、球の粘性力Ｆｈｉ、壁面反発力ＦＲｉｌ、及び流体力学的相互作用力ＦＰｉｋの全てを含む場合について説明するが、これに限らず、第３の力は、球の引っ張り復元力ＦＳｉｊ、球の粘性力Ｆｈｉ、壁面反発力ＦＲｉｌ、及び流体力学的相互作用力ＦＰｉｋの少なくとも１つを含むようにしてもよい。
【００３５】
なお、曲げトルクＴｂｉｊの等価な力Ｆｂｉｊ、粘性トルクＴｈｉの等価な力Ｆh_torqueｉｊに摩擦力の影響が含まれるため、上記（２）式の摩擦力ｆｉｊの算出が不要となる。
【００３６】
＜本発明の実施の形態に係る簡易運動解析装置の構成＞
【００３７】
本発明の実施の形態に係る簡易運動解析装置の構成について説明する。図１に示すように、本実施の形態に係る簡易運動解析装置１０は、ＣＰＵと、ＲＡＭと、後述する発話文抽出処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この簡易運動解析装置１０は、機能的には図１に示すように、モデル化部１２と、引っ張り復元力算出部１４、粘性力算出部１６、回転角速度算出部１８と、粘性トルク算出部２０と、粘性トルク変換部２２と、曲げトルク算出部２４と、曲げトルク変換部２６と、壁面反発力算出部２８と、流体力学的相互作用力算出部３０と、解析部３２と、を含んで構成されている。
【００３８】
モデル化部１２は、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化する。例えば解析対象が前述した図３に示すような構造体５０であった場合、構造体５０を複数の球５２の集合体としてモデル化する。
【００３９】
引っ張り復元力算出部１４は、引っ張り復元力ＦＳｉｊを算出する。
【００４０】
粘性力算出部１６は、粘性力Ｆｈｉを算出する。
【００４１】
回転角速度算出部１８は、球の回転角速度を算出する。具体的には、図５に示すように、結合した２つの球ｉ，ｊの中心を結ぶ線分の中心部に破線で示す仮想球Ｉｍが存在すると仮定し、現時刻における仮想球Ｉｍの回転角速度ωＩｍを下記（１０）式により近似する。
【００４２】
【数１０】
2018088085.tif 000012
・・・（１０）
【００４３】
ここで、ｒｉｊは、球ｉから球ｊへ向かう方向ベクトルである。
【００４４】
粘性トルク算出部２０は、回転角速度算出部１８が算出した回転角速度に基づいて流体から球に働く粘性トルクを算出する。具体的には、まず仮想球Ｉｍの位置における流体の渦度ベクトルΩ（ｒＩｍ）を算出する。渦度ベクトルの求め方は任意であり、例えば球ｉの渦度ベクトルΩ（ｒｉ）および球ｊの渦度ベクトルΩ（ｒｊ）の平均値を渦度ベクトルΩ（ｒＩｍ）としてもよい。
【００４５】
そして、算出した渦度ベクトルΩ（ｒＩｍ）と上記（１０）式で算出した球Ｉｍの回転角速度を用いて、上記（８）式により球Ｉｍの粘性トルクＴｈＩｍを算出する。
【００４６】
粘性トルク変換部２２は、粘性トルク算出部２０が算出した粘性トルクを等価な第１の力に変換する。
【００４７】
まず、仮想球Ｉｍにおける粘性トルク粘性トルクＴｈＩｍは、両隣の球ｉ，ｊに対して、等価な力を与えると仮定する。すなわち、図６に示すように、球ｉ，ｊに対して、粘性トルク粘性トルクＴｈＩｍに相当する力Ｆｈ＿ｔｏｒｑｕｅｉｊ、Ｆｈ＿ｔｏｒｑｕｅｊｉを下記（１１）式によって算出する。
【００４８】
【数１１】
2018088085.tif 000013
・・・（１１）
【００４９】
なお、仮想球Ｉｍの回転角速度を近似する場合に限らず、球ｉ、ｊの回転角速度を近似するようにしてもよい。また、上述の説明において、仮想球Ｉｍの存在を仮定したが、必ずしも球を想定する必要はなく、連結された２つの球の中心間距離と等しい長さを有する円柱などの他の形状の存在を仮定してもよい。
【００５０】
曲げトルク算出部２４は、球の相対位置関係に基づいて曲げトルクを算出する。具体的には、図７に示すように、３個の連結した球ｉ，ｊ，ｋの相対位置関係に基づいて、３個の球ｉ，ｊ，ｋがなす曲率半径ρを下記（１２）式により算出する。
【００５１】
【数１２】
2018088085.tif 000014
・・・（１２）
【００５２】
ここで、φは、図７に示すように、球ｉの中心と球ｊの中心を結ぶ線と、球ｊの中心と球ｋの中心を結ぶ線が成す角度である。
【００５３】
そして、上記（１２）式により算出した曲率半径ρを用いて、図８に示すように、球ｉおよびｋに働く曲げトルクＴｂｉｊ、Ｔｂｋｊを下記（１３）式により算出する。
【００５４】
【数１３】
2018088085.tif 000015
・・・（１３）
【００５５】
ここで、Ｉは断面２次モーメントである。
【００５６】
曲げトルク変換部２６は、曲げトルク算出部２４により算出された曲げトルクを等価な第２の力に変換する。
【００５７】
球ｉとｊを考えた場合、図９に示すように、球ｉは、曲げトルクＴｂｉｊによって球ｊから下記（１４）式で表される力Ｆｂｉｊを受ける。
【００５８】
【数１４】
2018088085.tif 000016
・・・（１４）
【００５９】
ここで、作用反作用の法則から、球ｊは球ｉから下記（１５）式で表される力Ｆｂｊｉを受ける。
【００６０】
【数１５】
2018088085.tif 000017
・・・（１５）
【００６１】
図９に示すように、球ｊと球ｋについても同様である。
【００６２】
壁面反発力算出部２８は、壁面反発力ＦＲｉｌを算出する。
流体力学的相互作用力算出部３０は、流体力学的相互作用力ＦＰｉｋを算出する。
解析部３２は、粘性トルクの等価な力、曲げトルクの等価な力、球の引っ張り復元力、球の粘性力、壁面反発力、及び流体力学的相互作用力を含む並進の運動方程式、すなわち上記（９）式を解くことにより、構造体５０の挙動、例えば球の位置、変形、及び配向の少なくとも１つを解析する。
【００６３】
＜本発明の実施形態に係る簡易運動解析装置の作用＞
【００６４】
次に、本発明の実施形態に係る簡易運動解析装置１０の作用について説明する。簡易運動解析装置１０は、簡易運動解析処理の実行が指示されると、図２に示す簡易運動解析処理ルーチンを実行する。
【００６５】
まず、ステップＳ１００では、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化する。
【００６６】
ステップＳ１０２では、粘性力を算出する。
【００６７】
ステップＳ１０４では、引っ張り復元力を算出する。
【００６８】
ステップＳ１０６では、球の回転角速度を算出する。
【００６９】
ステップＳ１０８では、ステップＳ１０６で算出した回転角速度に基づいて流体から球に働く粘性トルクを算出する。
【００７０】
ステップＳ１１０では、ステップＳ１０８で算出した粘性トルクを等価な第１の力に変換する。
【００７１】
ステップＳ１１２では、球の相対位置関係に基づいて曲げトルクを算出する。
【００７２】
ステップＳ１１４では、ステップＳ１１２で算出された曲げトルクを等価な第２の力に変換する。
ステップＳ１１６では、壁面反発力を算出する。
ステップＳ１１８では、流体力学的相互作用力を算出する。
【００７３】
ステップＳ１２０では、並進の運動方程式である上記（９）式を解くことにより、構造体５０の挙動、例えば球の位置、変形、及び配向の少なくとも１つを解析する。
【００７４】
なお、ステップＳ１０２〜Ｓ１２０は、各球について実行される。従って、ステップＳ１０４〜Ｓ１１０の処理は、結合する２個の球の組み合わせの全てに対して実行され、ステップＳ１１２〜Ｓ１１４の処理は、結合する３個の球の組み合わせの全てに対して行われる。
【００７５】
このように、本実施形態では、並進の運動方程式のみを解くため、上記（２）、（３）式の計算を省略することができ、上記（１）式の摩擦力ｆｉｊの算出も不要となる。従って、計算用のメモリ容量及び計算量を低減することができる。
【００７６】
＜実施例＞
【００７７】
次に、本発明の実施例について説明する。
【００７８】
まず、図１０に示すように、長さＬ［μｍ］の繊維を、半径ａ＝５［μｍ］の球の連結体としてモデル化した。このとき、繊維のアスペクト比ｒは、ｒ＝Ｌ／２ａである。また、１本の繊維を構成する球の個数Ｎは、Ｎ＝Ｌ／２ａである。この繊維の中心座標を（ｘ，ｙ，ｚ）＝（０，０，０）とし、時刻ｔ＝０においてｘ軸に平行に配向しているものとする。そして、この繊維に対して、
【００７９】
【数１６】
2018088085.tif 000018
なる単純せん断流動場を与える。ただし、ｕ，ｖ，ｗはそれぞれｘ，ｙ，ｚ方向の流速である。流動場の粘度はη＝１００［Ｐａ・ｓ］、せん断速度は
【数１７】
2018088085.tif 000019
とする。ただし、壁面反発力及び流体力学的相互作用力の算出は省略する。
【００８０】
はじめに、繊維の引張弾性率をＥ＝７０［ＧＰａ］とし、アスペクト比ｒを変化させた場合の繊維の回転挙動を計算した。本計算条件下において、繊維の曲げ変形は小さく、剛体回転挙動を示した。
【００８１】
図１１には、繊維の回転周期に関する、従来技術と本発明との誤差を示した。さらに、計算時間の高速化割合を図１２に示した。ただし、高速化割合とは、（従来技術を用いた場合の計算時間）÷（本発明を用いた場合の計算時間）で与えられる。また、ここでいう従来技術とは、並進及び回転の運動方程式を両方解くことで構造体の挙動を解析する技術である。
【００８２】
図１１より、本発明を用いた場合、従来技術との誤差は６％以下である。さらに、アスペクト比ｒが増加するに伴って誤差は低下していることがわかる。また、図１２より、アスペクト比ｒの増加に比例して計算高速化割合が上昇しており、アスペクト比１００においては１０倍以上の高速化割合が得られている。
【００８３】
以上より、従来技術と同等の精度を維持しつつ、大幅な高速化が達成できていることが判った。
【００８４】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【００８５】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能であるし、ネットワークを介して提供することも可能である。
【符号の説明】
【００８６】
１０   簡易運動解析装置
１２   モデル化部
１４   復元力算出部
１６   粘性力算出部
１８   回転角速度算出部
２０   粘性トルク算出部
２２   粘性トルク変換部
２４   トルク算出部
２６   トルク変換部
２８   壁面反発力算出部
３０   流体力学的相互作用力算出部
３２   解析部
５０   構造体
５２   球
【図１】
2018088085.tif 000020
【図２】
2018088085.tif 000021
【図３】
2018088085.tif 000022
【図４】
2018088085.tif 000023
【図５】
2018088085.tif 000024
【図６】
2018088085.tif 000025
【図７】
2018088085.tif 000026
【図８】
2018088085.tif 000027
【図９】
2018088085.tif 000028
【図１０】
2018088085.tif 000029
【図１１】
2018088085.tif 000030
【図１２】
2018088085.tif 000031