JP6740879B2

JP6740879B2 - 流体中における構造体の簡易運動解析方法、装置、及びプログラム

Info

Publication number: JP6740879B2
Application number: JP2016230398A
Authority: JP
Inventors: 俊貴笹山; 昌英稲垣
Original assignee: Toyota Central R&D Labs Inc
Current assignee: Toyota Central R&D Labs Inc
Priority date: 2016-11-28
Filing date: 2016-11-28
Publication date: 2020-08-19
Anticipated expiration: 2036-11-28
Also published as: JP2018088085A

Description

本発明は、流体中における構造体の簡易運動解析方法、装置、及びプログラムに関する。

従来、流体中の媒体の運動解析方法として、特許文献１、２、及び非特許文献１に開示された方法がある。これらの方法は、図３に示すように、任意形状の構造体５０を球５２の集合体としてモデル化する方法であり、各球５２について並進および回転の運動方程式を数値解析によって解くことにより、構造体５０の運動をシミュレートする。

並進の運動方程式は下記（１）式で与えられ、回転の運動方程式は下記（２）式で与えられる。

・・・（１）

・・・（２）

ここで、ｍは球の質量、ｒ_ｉは球ｉの位置ベクトル、Ｆ^Ｓ _ｉｊは引っ張り復元力、ｆ_ｉｊは摩擦力、Ｆ^ｈ _ｉは粘性力、Ｆ^Ｐ _ｉｋは球同士の流体力学的相互作用力、Ｆ^Ｒ _ｉｌは壁面反発力、ａは球の半径、θ_ｉは球ｉの配向ベクトル、Ｔ^ｂ _ｉｊは曲げトルク、Ｔ^ｔ _ｉｊはねじりトルク、ｎ_ｉｊは球ｉから球ｊへ向かう単位方向ベクトル、Ｔ^ｈ _ｉは粘性トルクである。

また、上記（１）、（２）式に加えて、結合した球同士がすべることなく運動するための拘束条件（すべり無し条件）を下記（３）式によって与える。

・・・（３）

ここで、ｖ_ｘは球ｘ（ｘはｉまたはｊ）の速度ベクトル、ω_ｘは球ｘ（ｘはｉまたはｊ）の角速度ベクトル、ｎ_ｊｉは球ｊから球ｉへ向かう単位方向ベクトルである。

また、上記（１）式のＦ^Ｓ _ｉｊ、Ｆ^ｈ _ｉは下記（４）、（５）式により算出される。

・・・（４）

・・・（５）

ここで、Ｅは弾性係数、ｒは図４（ａ）に示すように、球ｉの中心と球ｊの中心を結ぶ線の長さである結合長、ｒ_０は初期結合長、ηは流体の粘度、Ｖ（ｒ_ｉ）は位置ベクトルｒ_ｉで表される位置における流体の速度ベクトルである。

また、上記（２）式の、Ｔ^ｂ _ｉｊ、Ｔ^ｔ _ｉｊ、Ｔ^ｈ _ｉは下記（６）〜（８）式により算出される。

・・・（６）
ここで、ｎ^⊥ _ｉｊは、ｎ_ｉｊに垂直な単位ベクトルである。

・・・（７）
ここで、Ｇはせん断弾性係数である。

・・・（８）

ここで、θ_ｂは図４（ｂ）に示すように、球ｉと球ｊとの結合角の変位、θ_ｔは図４（ｃ）に示すように、球ｉと球ｊとのねじり角の変位、Ω（ｒ_ｉ）は位置ベクトルｒ_ｉで表される位置における流体の渦度ベクトルである。

特開２０００−３２２４０７号公報特許平１１−２５０６９号公報 S. Yamamoto, T. Matsuoka, "A method for dynamic simulation of rigid and flexible fibers in a flow field", J. Chem. Phys. 98, 644 (1993)

しかしながら、上記従来技術では、並進の運動方程式である上記（１）式、回転の運動方程式である上記（２）式の両方を解く必要があり、計算時間が増大する、という問題があった。

本発明は上記事情に鑑み成されたものであり、本発明の目的は、複数の球の集合体としてモデル化した構造体の流体中の挙動を、並進及び回転の２つの運動方程式を解くことで解析する場合と比較して、計算時間を低減させることができる流体中における構造体の簡易運動解析方法、装置、及びプログラムを提供することにある。

上記目的を達成するために、請求項１に係る発明の流体中における構造体の簡易運動解析方法は、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するステップと、前記球の回転角速度を算出するステップと、前記回転角速度に基づいて、流体から前記球に働く粘性トルクを算出するステップと、前記粘性トルクを等価な第１の力に変換するステップと、前記球の相対位置関係に基づいて曲げトルクを算出するステップと、前記曲げトルクを等価な第２の力に変換するステップと、前記第１の力、前記第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、前記構造体の挙動を解析するステップと、をコンピュータが実行する。

なお、請求項２に記載したように、前記球の回転角速度を算出するステップは、連結された２個の球の相対位置および前記球の速度に基づいて算出するようにしてもよい。

また、請求項３に記載したように、前記球の回転角速度を算出するステップは、連結された２個の球の間に設定した仮想球の回転角速度を算出し、前記粘性トルクを算出するステップは、前記仮想球の位置における渦度ベクトルを算出し、算出した前記渦度ベクトル及び前記仮想球の回転角速度に基づいて、前記粘性トルクを算出するようにしてもよい。

請求項４記載の発明の流体中における構造体の簡易運動解析装置は、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するモデル化部と、前記球の回転角速度を算出する回転角速度算出部と、前記回転角速度に基づいて流体から前記球に働く粘性トルクを算出する粘性トルク算出部と、前記粘性トルクを等価な第１の力に変換する粘性トルク変換部と、前記球の相対位置関係に基づいて曲げトルクを算出する曲げトルク算出部と、前記曲げトルクを等価な第２の力に変換する曲げトルク変換部と、前記第１の力、前記第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、前記構造体の挙動を解析する解析部と、を含む。

請求項５記載の発明の流体中における構造体の簡易運動解析プログラムは、コンピュータに、請求項１〜３の何れか１項に記載の流体中における構造体の簡易運動解析方法の各ステップを実行させるための流体中における構造体の簡易運動解析プログラムである。

本発明によれば、複数の球の集合体としてモデル化した構造体の流体中の挙動を、並進及び回転の２つの運動方程式を解くことで解析する場合と比較して、計算時間を低減させることができる。

簡易運動解析装置の構成を示すブロック図である。簡易運動解析処理ルーチンのフローチャートである。構造体のモデル化について説明するための図である。並進及び運動の方程式のパラメータについて説明するための図である。回転角速度の算出について説明するための図である。粘性トルクの変換について説明するための図である。曲げトルクの算出について説明するための図である。曲げトルクの算出について説明するための図である。曲げトルクの変換について説明するための図である。実施例に係る構造体について説明するための図である。実施例に係る繊維の回転周期に関する、従来技術と本発明との誤差について説明するための図である。実施例に係る従来技術と本発明の計算時間の高速化割合について説明するための図である。

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

＜本発明の概要＞

まず、本発明の概要について説明する。

本実施形態では、前述したように、図３に示すように任意形状の構造体５０を球５２の集合体としてモデル化する。そして、各球５２について、前述した回転の運動方程式である上記（２）式は解かずに、下記（９）式の並進の運動方程式のみを解くことによって構造体５０の挙動を解析する。

・・・（９）

上記（９）式の右辺第５項のＦ^ｂ _ｉｊは、前述した回転の運動方程式である上記（２）式の曲げトルクＴ^ｂ _ｉｊを等価な力（第２の力）に変換したものである。また、上記（９）式の右辺第６項のＦ^h_torque _ｉｊは、前述した回転の運動方程式である上記（２）式の粘性トルクＴ^ｈ _ｉを等価な力（第１の力）に変換したものである。なお、本実施形態では、第３の力として、球の引っ張り復元力Ｆ^Ｓ _ｉｊ、球の粘性力Ｆ^ｈ _ｉ、壁面反発力Ｆ^Ｒ _ｉｌ、及び流体力学的相互作用力Ｆ^Ｐ _ｉｋの全てを含む場合について説明するが、これに限らず、第３の力は、球の引っ張り復元力Ｆ^Ｓ _ｉｊ、球の粘性力Ｆ^ｈ _ｉ、壁面反発力Ｆ^Ｒ _ｉｌ、及び流体力学的相互作用力Ｆ^Ｐ _ｉｋの少なくとも１つを含むようにしてもよい。

なお、曲げトルクＴ^ｂ _ｉｊの等価な力Ｆ^ｂ _ｉｊ、粘性トルクＴ^ｈ _ｉの等価な力Ｆ^h_torque _ｉｊに摩擦力の影響が含まれるため、上記（２）式の摩擦力ｆ_ｉｊの算出が不要となる。

＜本発明の実施の形態に係る簡易運動解析装置の構成＞

本発明の実施の形態に係る簡易運動解析装置の構成について説明する。図１に示すように、本実施の形態に係る簡易運動解析装置１０は、ＣＰＵと、ＲＡＭと、後述する発話文抽出処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この簡易運動解析装置１０は、機能的には図１に示すように、モデル化部１２と、引っ張り復元力算出部１４、粘性力算出部１６、回転角速度算出部１８と、粘性トルク算出部２０と、粘性トルク変換部２２と、曲げトルク算出部２４と、曲げトルク変換部２６と、壁面反発力算出部２８と、流体力学的相互作用力算出部３０と、解析部３２と、を含んで構成されている。

モデル化部１２は、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化する。例えば解析対象が前述した図３に示すような構造体５０であった場合、構造体５０を複数の球５２の集合体としてモデル化する。

引っ張り復元力算出部１４は、引っ張り復元力Ｆ^Ｓ _ｉｊを算出する。

粘性力算出部１６は、粘性力Ｆ^ｈ _ｉを算出する。

回転角速度算出部１８は、球の回転角速度を算出する。具体的には、図５に示すように、結合した２つの球ｉ，ｊの中心を結ぶ線分の中心部に破線で示す仮想球Ｉｍが存在すると仮定し、現時刻における仮想球Ｉｍの回転角速度ω_Ｉｍを下記（１０）式により近似する。

・・・（１０）

ここで、ｒ_ｉｊは、球ｉから球ｊへ向かう方向ベクトルである。

粘性トルク算出部２０は、回転角速度算出部１８が算出した回転角速度に基づいて流体から球に働く粘性トルクを算出する。具体的には、まず仮想球Ｉｍの位置における流体の渦度ベクトルΩ（ｒ_Ｉｍ）を算出する。渦度ベクトルの求め方は任意であり、例えば球ｉの渦度ベクトルΩ（ｒ_ｉ）および球ｊの渦度ベクトルΩ（ｒ_ｊ）の平均値を渦度ベクトルΩ（ｒ_Ｉｍ）としてもよい。

そして、算出した渦度ベクトルΩ（ｒ_Ｉｍ）と上記（１０）式で算出した球Ｉｍの回転角速度を用いて、上記（８）式により球Ｉｍの粘性トルクＴ^ｈ _Ｉｍを算出する。

粘性トルク変換部２２は、粘性トルク算出部２０が算出した粘性トルクを等価な第１の力に変換する。

まず、仮想球Ｉｍにおける粘性トルク粘性トルクＴ^ｈ _Ｉｍは、両隣の球ｉ，ｊに対して、等価な力を与えると仮定する。すなわち、図６に示すように、球ｉ，ｊに対して、粘性トルク粘性トルクＴ^ｈ _Ｉｍに相当する力Ｆ^{ｈ＿ｔｏｒｑｕｅ} _ｉｊ、Ｆ^{ｈ＿ｔｏｒｑｕｅ} _ｊｉを下記（１１）式によって算出する。

・・・（１１）

なお、仮想球Ｉｍの回転角速度を近似する場合に限らず、球ｉ、ｊの回転角速度を近似するようにしてもよい。また、上述の説明において、仮想球Ｉｍの存在を仮定したが、必ずしも球を想定する必要はなく、連結された２つの球の中心間距離と等しい長さを有する円柱などの他の形状の存在を仮定してもよい。

曲げトルク算出部２４は、球の相対位置関係に基づいて曲げトルクを算出する。具体的には、図７に示すように、３個の連結した球ｉ，ｊ，ｋの相対位置関係に基づいて、３個の球ｉ，ｊ，ｋがなす曲率半径ρを下記（１２）式により算出する。

・・・（１２）

ここで、φは、図７に示すように、球ｉの中心と球ｊの中心を結ぶ線と、球ｊの中心と球ｋの中心を結ぶ線が成す角度である。

そして、上記（１２）式により算出した曲率半径ρを用いて、図８に示すように、球ｉおよびｋに働く曲げトルクＴ^ｂ _ｉｊ、Ｔ^ｂ _ｋｊを下記（１３）式により算出する。

・・・（１３）

ここで、Ｉは断面２次モーメントである。

曲げトルク変換部２６は、曲げトルク算出部２４により算出された曲げトルクを等価な第２の力に変換する。

球ｉとｊを考えた場合、図９に示すように、球ｉは、曲げトルクＴ^ｂ _ｉｊによって球ｊから下記（１４）式で表される力Ｆ^ｂ _ｉｊを受ける。

・・・（１４）

ここで、作用反作用の法則から、球ｊは球ｉから下記（１５）式で表される力Ｆ^ｂ _ｊｉを受ける。

・・・（１５）

図９に示すように、球ｊと球ｋについても同様である。

壁面反発力算出部２８は、壁面反発力Ｆ^Ｒ _ｉｌを算出する。
流体力学的相互作用力算出部３０は、流体力学的相互作用力Ｆ^Ｐ _ｉｋを算出する。
解析部３２は、粘性トルクの等価な力、曲げトルクの等価な力、球の引っ張り復元力、球の粘性力、壁面反発力、及び流体力学的相互作用力を含む並進の運動方程式、すなわち上記（９）式を解くことにより、構造体５０の挙動、例えば球の位置、変形、及び配向の少なくとも１つを解析する。

＜本発明の実施形態に係る簡易運動解析装置の作用＞

次に、本発明の実施形態に係る簡易運動解析装置１０の作用について説明する。簡易運動解析装置１０は、簡易運動解析処理の実行が指示されると、図２に示す簡易運動解析処理ルーチンを実行する。

まず、ステップＳ１００では、流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化する。

ステップＳ１０２では、粘性力を算出する。

ステップＳ１０４では、引っ張り復元力を算出する。

ステップＳ１０６では、球の回転角速度を算出する。

ステップＳ１０８では、ステップＳ１０６で算出した回転角速度に基づいて流体から球に働く粘性トルクを算出する。

ステップＳ１１０では、ステップＳ１０８で算出した粘性トルクを等価な第１の力に変換する。

ステップＳ１１２では、球の相対位置関係に基づいて曲げトルクを算出する。

ステップＳ１１４では、ステップＳ１１２で算出された曲げトルクを等価な第２の力に変換する。
ステップＳ１１６では、壁面反発力を算出する。
ステップＳ１１８では、流体力学的相互作用力を算出する。

ステップＳ１２０では、並進の運動方程式である上記（９）式を解くことにより、構造体５０の挙動、例えば球の位置、変形、及び配向の少なくとも１つを解析する。

なお、ステップＳ１０２〜Ｓ１２０は、各球について実行される。従って、ステップＳ１０４〜Ｓ１１０の処理は、結合する２個の球の組み合わせの全てに対して実行され、ステップＳ１１２〜Ｓ１１４の処理は、結合する３個の球の組み合わせの全てに対して行われる。

このように、本実施形態では、並進の運動方程式のみを解くため、上記（２）、（３）式の計算を省略することができ、上記（１）式の摩擦力ｆ_ｉｊの算出も不要となる。従って、計算用のメモリ容量及び計算量を低減することができる。

＜実施例＞

次に、本発明の実施例について説明する。

まず、図１０に示すように、長さＬ［μｍ］の繊維を、半径ａ＝５［μｍ］の球の連結体としてモデル化した。このとき、繊維のアスペクト比ｒは、ｒ＝Ｌ／２ａである。また、１本の繊維を構成する球の個数Ｎは、Ｎ＝Ｌ／２ａである。この繊維の中心座標を（ｘ，ｙ，ｚ）＝（０，０，０）とし、時刻ｔ＝０においてｘ軸に平行に配向しているものとする。そして、この繊維に対して、

なる単純せん断流動場を与える。ただし、ｕ，ｖ，ｗはそれぞれｘ，ｙ，ｚ方向の流速である。流動場の粘度はη＝１００［Ｐａ・ｓ］、せん断速度は

とする。ただし、壁面反発力及び流体力学的相互作用力の算出は省略する。

はじめに、繊維の引張弾性率をＥ＝７０［ＧＰａ］とし、アスペクト比ｒを変化させた場合の繊維の回転挙動を計算した。本計算条件下において、繊維の曲げ変形は小さく、剛体回転挙動を示した。

図１１には、繊維の回転周期に関する、従来技術と本発明との誤差を示した。さらに、計算時間の高速化割合を図１２に示した。ただし、高速化割合とは、（従来技術を用いた場合の計算時間）÷（本発明を用いた場合の計算時間）で与えられる。また、ここでいう従来技術とは、並進及び回転の運動方程式を両方解くことで構造体の挙動を解析する技術である。

図１１より、本発明を用いた場合、従来技術との誤差は６％以下である。さらに、アスペクト比ｒが増加するに伴って誤差は低下していることがわかる。また、図１２より、アスペクト比ｒの増加に比例して計算高速化割合が上昇しており、アスペクト比１００においては１０倍以上の高速化割合が得られている。

以上より、従来技術と同等の精度を維持しつつ、大幅な高速化が達成できていることが判った。

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

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能であるし、ネットワークを介して提供することも可能である。

１０簡易運動解析装置
１２モデル化部
１４復元力算出部
１６粘性力算出部
１８回転角速度算出部
２０粘性トルク算出部
２２粘性トルク変換部
２４トルク算出部
２６トルク変換部
２８壁面反発力算出部
３０流体力学的相互作用力算出部
３２解析部
５０構造体
５２球

Claims

流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するステップと、
前記球の回転角速度を算出するステップと、
前記回転角速度に基づいて、流体から前記球に働く粘性トルクを算出するステップと、
前記粘性トルクを等価な第１の力に変換するステップと、
前記球の相対位置関係に基づいて曲げトルクを算出するステップと、
前記曲げトルクを等価な第２の力に変換するステップと、
前記第１の力、前記第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、前記構造体の挙動を解析するステップと、
をコンピュータが実行する流体中における構造体の簡易運動解析方法。
前記球の回転角速度を算出するステップは、連結された２個の球の相対位置および前記球の速度に基づいて算出する
請求項１記載の流体中における構造体の簡易運動解析方法。
前記球の回転角速度を算出するステップは、連結された２個の球の間に設定した仮想球の回転角速度を算出し、
前記粘性トルクを算出するステップは、前記仮想球の位置における渦度ベクトルを算出し、算出した前記渦度ベクトル及び前記仮想球の回転角速度に基づいて、前記粘性トルクを算出する
請求項２記載の流体中における構造体の簡易運動解析方法。
流体中に含まれる構造体を、連結した複数の球の集合体としてモデル化するモデル化部と、
前記球の回転角速度を算出する回転角速度算出部と、
前記回転角速度に基づいて流体から前記球に働く粘性トルクを算出する粘性トルク算出部と、
前記粘性トルクを等価な第１の力に変換する粘性トルク変換部と、
前記球の相対位置関係に基づいて曲げトルクを算出する曲げトルク算出部と、
前記曲げトルクを等価な第２の力に変換する曲げトルク変換部と、
前記第１の力、前記第２の力、並びに、前記第１の力及び前記第２の力と異なる第３の力を含む並進の運動方程式を解くことにより、前記構造体の挙動を解析する解析部と、
を含む流体中における構造体の簡易運動解析装置。
コンピュータに、請求項１〜３の何れか１項に記載の流体中における構造体の簡易運動解析方法の各ステップを実行させるための流体中における構造体の簡易運動解析プログラム。