JP7102741B2 - 流体解析装置、流体解析方法、及び流体解析プログラム - Google Patents

Description

本発明は、流体解析装置、流体解析方法、及び流体解析プログラムに関する。

従来よりコンピュータによって流体の流れを数値解析する数値流体力学（ＣＦＤ：Computational Fluid Dynamics）を用いた技術が知られている。例えば、固体モデルと、有限要素でモデル化した流体モデルとを用いて、流体の流れを解析するにあたって、流体モデルのせん断応力を求める技術が提案されている（例えば、特許文献１参照）。この技術では、直交格子法の一例であるＩｍｍｅｒｓｅｄ ｂｏｕｎｄａｒｙ（ＩＢ）法に基づき、ナビエ・ストークス（Ｎａｖｉｅｒ・Ｓｔｏｋｅｓ）の式により流体モデルの運動を表している。そして、固体モデルと流体モデルとの境界に作用する流体モデルのせん断応力を求めている。この境界上のせん断応力を算出する際には、流体領域の要素の重心点の情報を用いて、境界上の点から流体領域へ延ばした垂線上の点に速度を補間し、その値を用いて境界内側（物体内部）の計算格子に仮想的な速度を外挿し、その値を参照して流れ場を計算する。

また、自動車や鉄道車両などの対象物に対する遮蔽板の流体力学的な影響を解析するため、計算格子中で物体が占める固体体積率を用いて流体と境界の速度差が０となるような仮想流体力を求め、求めた仮想流体力を用いて物体形状を表現する技術が提案されている（例えば、特許文献２参照）。この技術では、境界近傍の計算格子では速度を補間して求めている。

さらに流体と剛体の連成シミュレーションを簡易的に行う技術が提案されている（例えば、特許文献３参照）。この技術では、剛体の形状データをもとに計算格子に含まれる固体体積率を計算し、計算された固体体積率の値に基づいて剛体領域の計算格子に剛体運動の速度を付加し、流れ場を計算する。

特開２０１４－１０２５７４号公報 特開２０１６－１６４７０６号公報 特開２００４－２１８８３号公報

しかしながら、特許文献１の技術では、境界の内側の計算格子に流速を補間し、境界に隣接した計算格子でも境界から離れた計算格子と同じように式を離散化するが、補間操作に三次元の高次式を用いなければならず、外挿アルゴリズムが複雑になり、計算負荷が増大する。また、せん断応力の算出に流れ場の情報を用いた境界上の速度勾配を用いるが、得られる速度勾配は法線方向の勾配のみであり、境界上の三次元方向の速度勾配を求めることは困難である。特許文献２の技術は、計算格子中で物体が占める固体体積率を用いて仮想流体力を算出するため、実際の境界近傍における流れ場を扱うものではなく、境界近傍の流れ場の計算精度は低下する。特許文献３の技術も、固体体積率から剛体領域の計算格子を判定し、その計算格子に剛体運動から求まる速度を付加するため、実際の境界近傍における流れ場を扱うものではなく、境界近傍の流れ場の計算精度は低下する。従って、計算負荷を抑制しつつ簡単な構成で、高精度に流体の流れを解析する流体解析装置を提供することには改善の余地がある。

本発明は、上記事実を考慮してなされたもので、計算負荷を抑制しつつ簡単な構成で、高精度に流体の流れを解析することができる流体解析装置、流体解析方法、及び流体解析プログラムを提供することを目的とする。

上記目的を達成するために、本発明の流体解析装置は、
物体と、前記物体の周囲に流れ、かつ粘度が空間的に変化する流体と、を含む空間を前記物体と前記流体との境界に一致せず、かつ直交する複数の格子線により複数の要素で分割した数値計算対象モデルを定義するモデル定義部と、
前記モデル定義部で定義された数値計算対象モデルを用いて予め定めた方向に連続する前記境界を含む複数の要素各々における情報に基づいて、前記境界の予め定めた方向の流体の流れに関係する前記流体の速度勾配を示す第１の物理量を、前記複数の格子線によりなる格子の中心点又は重心点から前記境界までの距離を用いて演算する第１演算部と、
前記第１演算部で演算された前記第１の物理量に基づいて、前記境界の前記予め定めた方向と交差する方向の流体の流れに関係する前記流体の速度勾配を示す第２の物理量を、前記第１の物理量と前記境界で成立する予め定めた速度勾配に関する関係式とを用いて演算する第２演算部と、
前記第１演算部で演算された前記第１の物理量、及び前記第２演算部で演算された前記第２の物理量に基づいて、前記境界を含む予め定めた注目領域を流れる前記流体の挙動を導出する導出部と、
を備えている。

本発明の流体解析装置によれば、第１演算部は、空間を複数の要素で分割して定義された数値計算対象モデルを用いて予め定めた方向に連続する物体と流体との境界を含む複数の要素各々における情報に基づいて、境界の予め定めた方向の流体の流れに関係する第１の物理量を演算する。第２演算部は、第１演算部で演算された第１の物理量に基づいて、境界の予め定めた方向と交差する方向の流体の流れに関係する第２の物理量を演算する。導出部は、第１の物理量、及び第２の物理量に基づいて、境界を含む予め定めた注目領域を流れる流体の挙動を導出する。このように、境界の予め定めた方向の流体の流れに関係する第１の物理量に基づき交差する方向の流体の流れに関係する第２の物理量を演算するので、計算負荷を抑制しつつ簡単な構成で、高精度に流体の流れを解析することができる。

前記第１の物理量及び第２の物理量は、流体の速度勾配である。流体の挙動を正しく捉えるためには、境界の近傍においても流体の運動方程式に基づいて、流体の速度を特定することが好ましい。そのためには、境界近傍の流体の速度を求める際に、境界上における流体の速度勾配を用いることが好ましい。

前記数値計算対象モデルは、直交する複数の格子線により複数の要素で分割されて定義される。数値計算対象モデルは数値計算による解析を容易に行うことができ、かつ物体及び流体の領域を的確に表現することが好ましい。そこで、数値計算対象モデルを、直交する複数の格子線により複数の要素で分割して定義することで、各領域を容易に扱うことができる。
なお、前記第１演算部は、前記第１の物理量を、
（Ｉ、Ｊ）を複数の格子線により分割される要素の前記予め定めた方向の位置及び前記予め定めた方向と交差する方向の位置とし、予め定めた方向をｘ方向とし、ｂを（Ｉ、Ｊ）の位置の要素における中心点又は重心点を通る前記予め定めた方向の直線と前記境界との交点であることを示す情報とし、前記ｘ方向の流体の速度勾配を（∂ｕ _ｉ ／∂ｘ） _ｂ とし、流速のｘ成分をｕとし、Δｘを要素のｘ方向の幅とし、ε _x を格子の中心点又は重心点から前記境界までの幅をΔｘで正規化した情報とするとき、

で表される式を用いて演算することができる。
また、前記第２演算部は、前記第２の物理量を、
予め定めた方向と交差する方向をｙ方向とｚ方向とし、前記ｙ方向と前記ｚ方向の流体の速度勾配を（∂ｕ／∂ｙ） _ｂ と（∂ｕ／∂ｚ） _ｂ とし、前記境界における接線方向の単位ベクトルをｓとｔ、前記物体が運動する際の回転角速度ベクトルをωとするとき、

で表される式を用いて演算することができる。

前記予め定めた方向と交差する方向は、相互に直交し、かつ各々前記予め定めた方向と直交する２方向である。
流体の流れを把握するには流体の流れに関する３次元方向の速度成分を特定することが好ましい。そこで、予め定めた方向の速度勾配に加えて、相互に直交し、かつ各々前記予め定めた方向と直交する２方向の速度勾配を利用することにより、流体の流れの３次元方向の速度成分を的確に把握することができる。

本発明の流体解析方法は、
コンピュータが、
物体と、前記物体の周囲に流れ、かつ粘度が空間的に変化する流体と、を含む空間を前記物体と前記流体との境界に非適合で、かつ直交する複数の格子線により複数の要素で分割して定義された数値計算対象モデルを用いて予め定めた方向に連続する前記境界を含む複数の要素各々における情報に基づいて、前記境界の予め定めた方向の流体の流れに関係する前記流体の速度勾配を示す第１の物理量を、前記複数の格子線によりなる格子の中心点又は重心点から前記境界までの距離を用いて演算し、
演算された前記第１の物理量に基づいて、前記境界の前記予め定めた方向と交差する方向の流体の流れに関係する前記流体の速度勾配を示す第２の物理量を、前記第１の物理量と前記境界で成立する予め定めた速度勾配に関する関係式とを用いて演算し、
演算された前記第１の物理量、及び演算された前記第２の物理量に基づいて、前記境界を含む予め定めた注目領域を流れる前記流体の挙動を導出する。

本発明の流体解析プログラムは、コンピュータを、請求項１～請求項４の何れか１項に記載された流体解析装置として機能させる。

以上説明したように本発明によれば、計算負荷を抑制しつつ簡単な構成で、高精度に流体の流れを解析することができる、という効果が得られる。

実施形態に係る流体解析装置の構成の一例を示すブロック図である。 直交格子を用いて物体領域と流体領域との境界を含む領域を分割した一例を示すイメージ図である。 剛体運動する二次元物体の一例を示すイメージ図である。 二方向に境界を持つ場合の一例を示すイメージ図である。 比較例の外挿法の説明図である。 実施形態に係るコンピュータによる流体解析装置により流体解析方法を実行するための処理の流れの一例を示すフローチャートである。 実施形態に係る流体解析装置で流体の流れを検証する第１の検証で検証対象とした流体が流れる物体の概略構造の一例を示すイメージ図である。 第１の検証の検証結果を示すイメージ図である。 比較法を含めた第１の検証の検証結果を示すイメージ図である。 実施形態に係る流体解析装置で流体の流れを検証する第２の検証で検証対象とした流体が流れる物体の概略構造の一例を示すイメージ図である。

以下、図面を参照して本発明に係る実施形態を説明する。
本実施形態は、剛体の移動によって作用する流体の流れを解析する場合の一例を説明する。

図１に、本実施形態に係る流体解析装置１０を、コンピュータにより実現する構成の一例を示す。
図１に示すように、流体解析装置１０として動作するコンピュータは、ＣＰＵ（Central Processing Unit）１２Ａ、ＲＡＭ（Random Access Memory）１２Ｂ、およびＲＯＭ（Read Only Memory）１２Ｃを備えた装置本体１２を含んで構成されている。ＲＯＭ１２Ｃは、剛体の移動によって作用する流体の挙動をシミュレーションする各種機能を実現するための流体解析プログラム１２Ｐを含んでいる。装置本体１２は、入出力インタフェース（Ｉ／Ｏ）１２Ｄを備えており、ＣＰＵ１２Ａ、ＲＡＭ１２Ｂ、ＲＯＭ１２Ｃ、及びＩ／Ｏ１２Ｄは各々コマンド及びデータを授受可能なようにバス１２Ｅを介して接続されている。また、Ｉ／Ｏ１２Ｄには、キーボード及びマウス等の数値入力、文字入力及びコマンド入力を可能とする入力部１２Ｆ、画像表示を可能とする表示部１２Ｇ、外部装置と通信する通信部１２Ｈ、及び各種データを記憶する記憶部１２Ｍが接続されている。

装置本体１２は、流体解析プログラム１２ＰがＲＯＭ１２Ｃから読み出されてＲＡＭ１２Ｂに展開され、ＲＡＭ１２Ｂに展開された流体解析プログラム１２ＰがＣＰＵ１２Ａによって実行されることで、流体解析装置１０として動作する。なお、流体解析プログラム１２Ｐは、流体の流れを解析する各種機能を実現するためのプロセスを含む（詳細は後述）。

本実施形態では、粘性係数が空間的に変化する非ニュートン流体を、本実施形態に係る流体の一例として説明する。

まず、本実施形態に係る流体解析装置１０により流体の流れを数値解析することを可能とするために、流れ場を離散化して物体と流体との境界における速度勾配を導出する方法について説明する。

流体の流れを解析する場合、流れ場の基礎方程式等を用いて流体の挙動を表現することが行われている。次に示す（１）式及び（２）式に本実施形態で用いる流れ場の基礎方程式の一例を示す。（１）式は非圧縮性流体の連続の式を示し、（２）式はナビエ・ストークスの式を示す。

ただし、（１）式中及び（２）式中の記号は次を示している。ｕ_ｉは流速、ｐは圧力、ρは密度を示す。τ_ｉｊは粘性応力テンソルを示す。この粘性応力テンソルτ_ｉｊは、粘度ηとひずみ速度テンソルＤ_ｉｊを用いて次に示す（３）式で表すことができる。

以下の説明では、上記の（３）式の右辺第１項、及び右辺第２項の各々の発散を取ったものを、それぞれ粘性項第１項、及び粘性項第２項と称する。

本実施形態では、境界に適合しない計算格子を用いて流体の流れを計算する解析手法を対象に説明する。本実施形態では、直交格子を用いて流体の流れを解析する。直交格子を用いて流体の流れを解析する技術の一例には、文献１（佐藤他、直交格子法における物体境界近傍の直接離散化法（速度場と圧力場の整合性を考慮した高精度化）、日本機械学会論文集Ｂ編 ７９（２０１３）６０５－６２１）に示す技術が挙げられる。

また、本実施形態では、空間離散化にコロケーション格子法を用いる。そして、境界に隣接した計算格子について基礎式を離散化する。以下の説明では、本実施形態に関連する粘性項に関する離散化について説明する。

図２に、直交格子を用いて物体領域２と流体領域３とを離散化した物体領域２と流体領域３との境界４を含む注目領域１を示す。図２に示す注目領域１は、Ｘ軸とＹ軸を通る二次元平面における物体領域２と流体領域３との一部の領域を示している。

図２に示すように、注目領域１は、物体領域２と流体領域３との境界４に適合しない計算格子を用いて複数の要素で定義される。物体領域２と流体領域３とを、直交格子を用いて分割、すなわち、直交する格子線５により分割した複数（図２の例では縦横に３ｘ３の９個）の計算格子６で定義した場合を示している。計算格子６の各々には、流体の流れを解析する場合において、流体の速度、圧力又は温度等を計算する計算点７が定義される。計算点７は計算格子６内の何れの位置に定義してもよいが、計算格子６を代表する代表点として機能するので、計算格子６の中心点又は重心点を設定することが好ましい。本実施形態では、計算点７を計算格子６の中心点とした場合を説明する。

なお、以下の説明では、複数の計算格子６の各々を区別して説明する場合には、位置を示す識別符号を付した計算格子（Ｉ，Ｊ）と表記して説明する。

また、以下の説明では、図２に示すように、ｘ方向に境界を持つ二次元の計算格子６を対象に、粘性項の離散化について説明する。すなわち、特定の計算格子６である計算格子（Ｉ，Ｊ）は、隣り合う計算格子６のうちｘ方向の計算格子（Ｉ－１，Ｊ）のみに境界４を持ち、ｙ方向には境界４を持たない。

ところで、粘度が空間的に一定であると仮定した場合、粘性項第１項のみを扱えばよいが、本実施形態では、粘度の空間分布を考慮して、粘性項第２項も含めて離散化を行う。図２の計算格子（Ｉ，Ｊ）における粘性項第１項の離散化式は、次の（４）式で表すことができる。

ただし、（４）式中の記号は次を示している。δ／δｘ_ｉは２次精度中心差分であることを示し、下付き添え字ｂは、計算格子（Ｉ，Ｊ）の計算点７を通るｘ軸方向の直線と境界４との交点（図２に三角形記号で示す位置、以下、境界点８という）における値であることを示す。

境界点８におけるｘ方向の速度勾配（∂ｕ_ｉ／∂ｘ）_ｂは、次の（５）式に示す片側３点差分式（空間２次精度）によって求めることができる。

次に、粘性項第２項を離散化する。離散化された粘性項第２項は、次の（６）式で表すことができる。なお、流速ｕ_ｉのｘ成分がｕであり、ｙ成分がｖである。

上記の（６）式の右辺第１項には，境界点８におけるｘ方向の速度勾配（∂ｕ／∂ｘ）_ｂ、及びｙ方向の速度勾配（∂ｕ／∂ｙ）_ｂが含まれる。従って、粘性項第２項を離散化するには、粘性項第１項の離散化式（（４）式）に現れた速度勾配（∂ｕ／∂ｘ）_ｂの他に、ｘ方向と異なる方向の速度勾配、すなわち速度勾配（∂ｕ／∂ｙ）_ｂが必要になる。ところが、ｙ方向の速度勾配は、境界点８からｙ方向に延ばした直線上に計算点７が存在しないため、速度勾配（∂ｕ／∂ｘ）_ｂの計算に用いた式（（５）式）と同じように求めることができない。
そこで、本実施形態では、境界点８におけるｙ方向の速度勾配を、境界４上における速度勾配の関係式と、計算されたｘ方向の速度勾配とを用いて導出する。以下に、境界を持つ特定方向の速度勾配と異なる方向の速度勾配、すなわち、ｘ方向と異なる方向の速度勾配の導出について説明する。

ここでは、説明を簡単にするため、剛体運動する物体について速度勾配を導出する場合を説明する。剛体運動する物体では、物体表面上の境界点では、次の（７）式及び（８）式の関係が成り立つ。

ただし、（７）式中及び（８）式中の記号は次を示している。∇ｕは速度勾配テンソル、ｓ、ｔは境界点８での接線方向（二方向）の単位ベクトル、ωは剛体運動の回転角速度ベクトルを示す。

ここで、（７）式について図３を参照して説明する。
図３に、剛体運動する二次元物体の一例を示す。

図３に示すように、物体表面上の点Ｐとその近傍点Ｑにおける速度ベクトルｕ_ｐ、ｕ_０はそれぞれ次の（７Ａ）式、及び（７Ｂ）式で表すことができる。

ｕ_Ｇとω_Ｇは剛体運動の速度ベクトルと角速度ベクトルであり、ｒ_Ｐとｒ_Ｑは剛体の回転中心点Ｇから点Ｐ、Ｑへの位置ベクトルである。ｓはＰＱ間の単位ベクトル（物体表面上の接線方向の単位ベクトル）であり、ｈはＰＱ間の距離である。（７Ａ）式及び（７Ｂ）式により、ＰＱ間の速度勾配、すなわち剛体運動する物体表面上のｓ方向の速度勾配は、次に示す（７Ｃ）式で表すことができる。

また、（７Ｃ）式はさらに、次に示す（７Ｄ）式で表すことができる。

なお、（８）式についても同様のため、説明を省略する。

ここで、流速ｕの成分について考えると上記の（７）式、及び（８）式はそれぞれ、（９）式及び（１０）式で表すことができる。

ただし、（９）式中及び（１０）式中の記号ｓ_ｘ、ｓ_ｖ、ｓ_ｚと、ｔ_ｘ、ｔ_ｖ、ｔ_ｚとは、それぞれｓとｔの成分を示す。（９）式及び（１０）式のｘ方向の速度勾配（∂ｕ／∂ｘ）_ｂは（５）式の外挿式によって既知と考えると、それ以外の成分である速度勾配（∂ｕ／∂ｙ）_ｂ、及び速度勾配（∂ｕ／∂ｚ）_ｂは、（９）式、及び（１０）式の連立方程式から、次に示す（１１）式で表すことができる。

同様に、流速ｕと直交する方向の流速ｖ、ｗの速度勾配のｙ、ｚ成分についてもそれぞれ（１２）式及び（１３）式で表すことができる。

以上説明したように、境界点８について、ｘ方向の速度勾配は、（５）式を用いて導出でき、ｙ方向及びｚ方向の速度勾配は、（１１）式～（１３）式を用いて導出できる。

なお、図２に示す例で、二次元（ｘ－ｙ）平面内の運動（回転角速度ω）の場合は、（１２）式及び（１３）式は、次の（１４）式及び（１５）式に示すように簡略化して表すことができる。

また、物体が静止している場合、又は物体運動が並進運動のみの場合には、回転角速度ωがゼロベクトルであり、（９）式及び（１０）式の右辺は０となる。

上記では、説明の便宜上、図２に示すｘ方向に速度勾配を持つ場合を説明したが、ｙ、ｚ方向に境界を持つ計算格子６の場合にも同様に離散化可能であることは勿論である。

また、上記では、ｘ方向のみに速度勾配を持つ場合（図２）を説明したが、複数の方向に境界を持つ場合についても、各方向に関して上記と同様の方法を独立に適用することができる。例えば、図４に示すように、ｘ方向、及びｙ方向の二方向に境界４を持つ場合、ｘ、ｙの各方向に関して上記と同様に適用することができる。

本実施形態は境界４上の境界点８（図２の三角形記号の位置）において、差分参照点（流体計算点）が存在する方向の速度勾配（図２の場合は（∂ｕ／∂ｘ）_ｂ）については流体計算点の情報を使った差分式から求め、その他の方向の速度勾配（図２の場合はｙ方向の速度勾配（∂ｕ／∂ｙ）_ｂ）については境界上で成り立つべき速度勾配の関係式を利用して求める。従って、ｙ方向の速度勾配（∂ｕ／∂ｙ）_ｂは、ｘ方向の速度勾配（∂ｕ／∂ｘ）_ｂと同じ空間精度が保証される。

なお、上記の（９）式及び（１０）式では「差分参照点が存在する方向の速度勾配（∂ｕ／∂ｘ）_ｂ」を２次精度の差分式（（５）式）で求めたが、空間精度を変えることも可能である。例えば、１次精度の場合は（５）式の代わりに次の（５Ａ）式を用いる。

本技術で得られた速度勾配は、流体の非ニュートン性を有する計算対象、粘度に温度依存性がある計算対象、ＬＥＳ（Large Eddy Simulation）など乱流モデルを使う計算対象など、粘性係数が空間的に分布する計算に広く利用可能である。また物体表面上での摩擦応力の算出時にも利用できる。

ところで、境界点で速度勾配（∂ｕ／∂ｙ）_ｂを算出する一般的な方法として、境界近傍の流体領域の計算点で各座標軸方向の速度勾配を求めておき、それを境界点に外挿する方法（外挿法）が挙げられる。
ここで、外挿法により境界点における速度勾配を算出する方法を、比較例として説明する。例えば、図２に示す速度勾配（∂ｕ／∂ｙ）_ｂに関して空間２次精度の外挿式を適用した場合は、（１６）式で表すことができる。

上記の（１６）式の右辺の速度勾配は、次に示す（１７）式から（１９）式を用いて算出することができる。

図５に、比較例の外挿法を適用する場合の物体領域２と流体領域３との境界４を含む注目領域１を示す。

図５に示すように、比較例の外挿法を適用した場合、本実施形態に係る速度勾配を導出する方法に比べて、多い個数の計算点を要する。すなわち、上記の（１７）式から（１９）式に示すように、比較例の外挿法では、速度勾配（∂ｕ／∂ｙ）_ｂの算出に図５に黒丸記号で示す合計７個の計算点の情報が必要となる。一方、本実施形態に係る速度勾配を導出する方法では、（５）式の差分式と解析的な関係式だけで速度勾配（∂ｕ／∂ｙ）_ｂを導出することができるため、図５にＸ記号で示す合計３個の計算点の情報のみで良い。また、本実施形態では、（５）式を用いた場合に、２次外挿法と同じ空間２次精度が得られる。すなわち、本実施形態に係る速度勾配を導出する方法は外挿法に比べて少ない計算点の情報で同等の精度を得ることができる。さらに、補間操作を最小限にすることによって計算の安定性も大幅に向上し、２次外挿法に比べて大きな時間刻みで安定に計算ができ、計算時間の短縮も達成される。

次に、上記説明した速度勾配を導出する方法を基にして、流体解析装置１０において実行される、剛体の移動によって作用する流体の挙動のシミュレーション、すなわち、剛体の移動によって作用する流体の流れを解析する流体解析方法を説明する。コンピュータで実現した流体解析装置１０のＲＯＭ１２Ｃには、本実施形態に係る流体解析方法を実行するための処理の流れ（プロセス）が流体解析プログラム１２Ｐとして記憶されている。

図６には、本実施形態に係る流体解析方法を実行するための処理の流れの一例が示されている。装置本体１２では、流体解析プログラム１２ＰがＲＯＭ１２Ｃから読み出されてＲＡＭ１２Ｂに展開され、ＲＡＭ１２Ｂに展開された流体解析プログラム１２ＰをＣＰＵ１２Ａが実行する。

まず、ステップＳ１００では、計算格子６、計算条件、及び初期値を示す情報の取得処理が実行される。

計算格子６を示す情報は、物体領域２と流体領域３とを分割する直交格子、すなわち、直交する格子線５の間隔が定義された情報である。この計算格子６に関する情報は、ＲＯＭ１２Ｃ又は記憶部１２Ｍに記憶される。また、計算条件を示す情報は、本実施形態で流体の流れを数値解析する場合の条件である。計算条件の一例には、物体モデルの流体モデルに対する位置情報、及び物体モデルの動きを定義する角速度等の情報、及び流体モデルの計算条件の一例には、流体モデルの境界条件、速度場計算の収束条件、及び反復回数等の予め定めた情報が含まれる。計算条件は、ＲＯＭ１２Ｃ又は記憶部１２Ｍに記憶される。また、初期値を示す情報は、本実施形態で流体の流れを数値解析する場合の最初の各種状態を示す情報である。初期値の一例には、物体モデルの解析当初の動きを定義する角速度等の情報、及び、流体の当初の速度、圧力及び温度等の情報が挙げられる。

また、ステップＳ１００では、物体を示す物体モデル、及び流体を示す流体モデルも取得する。

物体モデルは、物体として例えば、剛体（非流体物）の構造を定義したモデルである。なお、物体モデルは、数値解析法により取り扱い可能な有限個の要素でモデル化することで定義してもよい。また、物体モデルには、座標値及び材料特性（例えば、粘度、密度、及び材質等）などの数値データが定義される。これらの数値は、ＲＯＭ１２Ｃ又は記憶部１２Ｍに記憶される。

流体モデルは、物体の外側で流体が流れる空間の構造を定義したモデルである。なお、流体モデルは、数値解析法により取り扱い可能な有限個の要素でモデル化することで定義してもよい。また、流体モデルには、座標値及び材料特性（例えば、粘度、密度、及び材質）などの数値データが定義される。これらの数値は、ＲＯＭ１２Ｃ又は記憶部１２Ｍに記憶される。

さらに、ステップＳ１００では、取得した情報を用いて物体と流体とを含む解析対象の全領域を定義する処理を実行する。すなわち、物体モデルによる物体領域２、流体モデルによる流体領域３、及び物体領域２と流体領域３との境界４を含む全領域について、直交格子を用いて分割、すなわち、直交する格子線５により分割した複数の計算格子６で定義する処理が実行される。なお、計算格子６の各々には、計算点７の位置も定義される。
従って、ステップＳ１００は、物体と、物体の周囲に流れる流体とを含む空間を複数の要素で分割した数値計算対象モデルを定義するモデル定義部として動作する機能を有する。

次に、ステップＳ１０２では、境界の移動に関する処理を実行する。このステップＳ１０２の処理は、物体が時々刻々と位置が移動する場合に伴って移動する境界４の位置を定義する処理である。従って、物体の位置が固定の場合は、ステップＳ１０２は省略される。また、本処理ルーチンの実行当初は、ステップＳ１０２では、初期値を用いればよい。

次のステップＳ１０４では、境界４と計算格子６の計算点７との交点である境界点８の抽出処理を実行する。このステップＳ１０４では、境界４を含む計算格子６の各々について境界点８を抽出する演算が行われる。

次のステップＳ１０６では、ステップＳ１０４で抽出された境界点８の各々における境界面上の接線方向の単位ベクトルを算出する処理を実行する。すなわち、上記の（７）式、及び（８）式を用いて、境界点８の各々における境界面上の接線方向の単位ベクトルを各々演算する。

次のステップＳ１０８では、境界点８の各々における一方向の速度勾配を計算する処理を実行する。すなわち、例えば、上記の（５）式を用いた片側３点差分式（空間２次精度）によって、境界点８におけるｕ成分のｘ方向の速度勾配（∂ｕ_ｉ／∂ｘ）_ｂと、ｖ成分のｘ方向の速度勾配（∂ｖ_ｉ／∂ｘ）_ｂと、ｗ成分のｘ方向の速度勾配（∂ｗ_ｉ／∂ｘ）_ｂと、を各々演算する。

次のステップＳ１１０では、境界点８の各々における残存方向の速度勾配を計算する処理を実行する。すなわち、上記の（１１）式から（１３）式）を用いて、ｕ成分の速度勾配（∂ｕ／∂ｙ）_ｂ、及び速度勾配（∂ｕ／∂ｚ）_ｂと、ｖ成分の速度勾配（∂ｖ／∂ｙ）_ｂ、及び速度勾配（∂ｖ／∂ｚ）_ｂと、ｗ成分の速度勾配（∂ｗ／∂ｙ）_ｂ、及び速度勾配（∂ｗ／∂ｚ）_ｂと、を各々演算する。

次のステップＳ１１２では、流体の速度場、及び圧力場を計算する処理を実行する。すなわち、ステップＳ１０８で求めた境界点８の各々における一方向の速度勾配、及びステップＳ１１０で求めたその他の方向の速度勾配を用いて、上記の（１）式から（３）式によって、流体の速度場、及び圧力場を演算する。

次のステップＳ１１４では、終了判断処理を実行し、肯定判断の場合はステップＳ１１６へ処理を移行し、否定判断の場合は、ステップＳ１０２へ処理を戻す。ステップＳ１１４の終了は、速度場計算の収束条件、ステップＳ１０２からステップ１１２の処理の反復回数、及び計算開始からの経過時間等の予め定めた情報に基づいて判断される。

次のステップＳ１１６では、ステップＳ１１２の計算結果を出力する処理を実行する。ステップＳ１１６では、ステップＳ１１２で求めた流体の速度場、及び圧力場の計算結果をそのまま出力してもよく、ステップＳ１１２で求めた流体の速度場、及び圧力場の計算結果に基づいて流体の流れの速度、及び圧力等の物理量を演算し、演算結果を出力してもよい。

このように、本実施形態によれば、剛体の移動によって作用する流体の速度勾配を求めて流体の挙動をシミュレーションすることができる。すなわち、外挿法に比べて少ない計算点の情報で高精度に流体の流れを解析することができる。

ここで、本実施形態に係る流体解析装置１０における流体の流れの解析を検証する。
第１の検証では、粘性係数が空間的に変化する非ニュートン流体を対象として、同軸二重円筒間に流れる流体における流体の流れの精度検証を行う。

図７に、第１の検証で検証対象とした流体が流れる物体（剛体）の概略構造の一例を示す。

図７に示すように、検証対象の領域２０は、内径ｒ_ｉ＝ＫＲ、外径ｒ_０＝Ｒの流路２２を備えた物体領域と、その流路２２を流れる流体領域とを含む。詳細には、流路２２の外周側の半径をＲとしたとき、ｘ方向に０．２Ｒ、ｙ方向及びｚ方向にそれぞれ２．４Ｒとし、ｙ－ｚ面内の計算領域の中央に内径ｒ_ｉ＝ＫＲ、外径ｒ_０＝Ｒの同軸二重円筒間流路が形成された流路２２からなる領域を検証対象の領域２０とする。

なお、内径と外径の半径比Ｋは０．６とする。壁面はすべりなし条件とし、流路２２内の流れは周方向に付加した平均圧力勾配ｄｐ／ｄθによって駆動される。検証対象の領域２０のｘ方向の両端面には周期条件を課す。非ニュートンモデルは、次の（２０）式で表されるＯｓｔｗａｌｄ－ｄｅ Ｗａｅｌｅモデル（べき乗則モデル）を用いる。

粘度ηは、ｎ＝１（ニュートン流体）の場合を除き、流れ場の速度勾配に依存して空間的に変化する。

第１の検証における検証対象である同軸二重円筒間の流路２２を流れる流体の挙動（流れ）には、解析解が存在する。解析解は、次に示す（２１）式により求めることができる。

ただし、（２１）式中の記号ｕ^＊ _θは無次元化された周方向流速である。

上記の（２１）式は以下のようにして導出できる。
層流域の同軸二重円筒間流れに関するナビエ・ストークス方程式は、円筒座標系では、次に示す（２１ａ）式で表すことができる。

粘性応力τ_ｒθは、べき乗則モデルを用いた場合には、次に示す（２１ｂ）式で表すことができる。

ただし、（２１ｂ）式中の記号γ_ｒθはせん断速度であり、周方向流速ｕ_θとの間に次の（２１ｃ）式に示す関係が成り立つ。

上記の（２１ａ）式から（２１ｃ）式を解くことにより、無次元化された周方向流速ｕ^＊ _θの解析解は（２１）式に示すように導出できる。この解析解を導出する技術の一例には、文献２（N. Prasanth, U.V. Shenoy, “Poiseuille flow of a power-law fluid between coaxial cylinders”, Journal of Applied Polymer Science, Vol. 46 (1992), pp. 1189-1194.）に記載された技術がある。

上記（２１）式における無次元化された周方向流速ｕ^＊ _θと、有次元の流速ｕ_θとの間は、次の（２２）式に示す関係を有する。

ただし、（２１）式中の記号λはせん断応力が０となる位置の半径距離である。なお、ｒ／Ｒ＝λでの流速分布の連続性は、次に示す（２３）式の関係から求めることができる。

図８に、流体解析装置１０を用いて第１の検証を行った検証結果を示す。図８では、格子間隔０．０２Ｒの条件について、各べき指数ｎでの周方向流速ｕ_θの半径方向の分布を打点（図８に示す黒丸）で示している。また、上記の（２１）式により求めた解析解の分布を実線で示している。

図８に示すように、流体解析装置１０を用いて非ニュートン流体の流れの解析する流体解析方法である本手法により予測された速度分布は、半径方向の何れの位置（何れのｎ）についても解析解と適合している。

図９に、流体解析装置１０を用いて第１の検証を行った検証結果と、外挿により流体の流れの解析する比較法の検証結果とを示す。図９では、ｎ＝０．５、格子間隔０．０８Ｒ（低格子解像度）の条件について、周方向流速ｕ_θの半径方向の分布が示されている。

図９に示す検証結果で用いた外挿により流体の流れの解析する比較法は、一般的な解析方法の一例として、線形補間に基づくＩｍｍｅｒｓｅｄ ｂｏｕｎｄａｒｙ（ＩＢ）法による流体の流れの解析する方法を比較法として用いている。この比較法の技術の一例には、特許文献１の技術、及び文献３（R. Mittal et al., A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries, J. Comput. Phys. 227 (2008) 4827-4825.）に記載された技術がある。文献３の技術では、境界上の流速と周囲の流体計算点の流速を用いて境界の内側（物体内部）の計算格子に仮想的な流速を外挿し、境界に隣接する計算格子ではそれらを用いることで境界から離れた計算格子と同様に流体の基礎式を離散化して計算している。

図９に示すように、比較法では、流速分布が解析解から乖離するのに対し、本手法では解析解に適合した結果が得られている。

次に、本実施形態に係る流体解析装置１０における流体の流れの解析を検証した第２の検証について説明する。
第２の検証では、流体の流れとして、剛体運動する物体周りの流れ場を対象に検証を行った。

図１０に、第２の検証で検証対象とした流体が流れる物体（剛体）の概略構造の一例を示す。

図１０に示すように、検証対象の領域３０は、筒内３２を同じ方向に一定回転する剛体を含む物体領域と、その筒内３２を流れる流体領域とを含む。詳細には、筒内３２で剛体運動する物体は、同じ方向に一定回転するスクリュ（図１０に斜線で示した領域の物体）であり、スクリュの回転によって周囲の流体が流動する。流体は溶融した樹脂相当の粘度特性とし、Ｃｒｏｓｓモデルで樹脂の非ニュートン性を考慮した。計算格子は直交等間隔であり、約０．０６ｍｍとした。スクリュの長径は約１５ｍｍ、回転数は２００ｒｐｍであり、このときの長径端部での最大周速は約０．１５７ｍ／ｓである。

第２の検証では、流体解析装置１０を用いて非ニュートン流体の流れの解析する流体解析方法である本手法による解析において、安定に計算できる計算時間刻みの上限値は、外挿による比較法を用いて解析した場合の約７倍の結果となり、計算時間が大幅に短縮された。

本実施形態に係る流体解析装置１０は、開示の技術の流体解析装置の一例である。また、図６に示すステップＳ１００処理は、モデル定義部で実行される機能の一例である。また、ステップＳ１０８の処理は、第１演算部で実行される機能の一例であり、ステップＳ１１０の処理は、第２演算部で実行される機能の一例であり、ステップＳ１１２の処理は、導出部で実行される機能の一例である。

なお、本実施形態に係る流体解析装置１０に含まれる装置本体１２は、構成する各構成要素を、上記で説明した各機能を有する電子回路等のハードウェアにより構築してもよく、構成する各構成要素の少なくとも一部を、コンピュータにより当該機能を実現するように構築してもよい。

また、本発明を実施の形態を用いて説明したが、本発明の技術的範囲は上記実施の形態に記載の範囲には限定されない。発明の要旨を逸脱しない範囲で上記実施の形態に多様な変更または改良を加えることができ、当該変更または改良を加えた形態も本発明の技術的範囲に含まれる。また、本明細書に記載された全ての文献、特許出願及び技術規格は、個々の文献、特許出願及び技術規格が参照により取り込まれることが具体的かつ個々に記された場合と同程度に、本明細書中に参照により取り込まれる。

１０ 流体解析装置
１２ 装置本体
１２Ａ ＣＰＵ
１２Ｂ ＲＡＭ
１２Ｃ ＲＯＭ
１２Ｐ 流体解析プログラム

Claims (6)

1. 物体と、前記物体の周囲に流れ、かつ粘度が空間的に変化する流体と、を含む空間を前記物体と前記流体との境界に一致せず、かつ直交する複数の格子線により複数の要素で分割した数値計算対象モデルを定義するモデル定義部と、
   前記モデル定義部で定義された数値計算対象モデルを用いて予め定めた方向に連続する前記境界を含む複数の要素各々における情報に基づいて、前記境界の予め定めた方向の流体の流れに関係する前記流体の速度勾配を示す第１の物理量を、前記複数の格子線によりなる格子の中心点又は重心点から前記境界までの距離を用いて演算する第１演算部と、
   前記第１演算部で演算された前記第１の物理量に基づいて、前記境界の前記予め定めた方向と交差する方向の流体の流れに関係する前記流体の速度勾配を示す第２の物理量を、前記第１の物理量と前記境界で成立する予め定めた速度勾配に関する関係式とを用いて演算する第２演算部と、
   前記第１演算部で演算された前記第１の物理量、及び前記第２演算部で演算された前記第２の物理量に基づいて、前記境界を含む予め定めた注目領域を流れる前記流体の挙動を導出する導出部と、
   を備えた流体解析装置。
2. 前記第１演算部は、前記第１の物理量を、
   （Ｉ、Ｊ）を複数の格子線により分割される要素の前記予め定めた方向の位置及び前記予め定めた方向と交差する方向の位置とし、予め定めた方向をｘ方向とし、ｂを（Ｉ、Ｊ）の位置の要素における中心点又は重心点を通る前記予め定めた方向の直線と前記境界との交点であることを示す情報とし、前記ｘ方向の流体の速度勾配を（∂ｕ _ｉ ／∂ｘ） _ｂ とし、流速のｘ成分をｕとし、Δｘを要素のｘ方向の幅とし、ε _x を格子の中心点又は重心点から前記境界までの幅をΔｘで正規化した情報とするとき、
   

   

   
   で表される式を用いて演算する
   請求項１に記載の流体解析装置。
3. 前記第２演算部は、前記第２の物理量を、
   予め定めた方向と交差する方向をｙ方向とｚ方向とし、前記ｙ方向と前記ｚ方向の流体の速度勾配を（∂ｕ／∂ｙ） _ｂ と（∂ｕ／∂ｚ） _ｂ とし、前記境界における接線方向の単位ベクトルをｓとｔ、前記物体が運動する際の回転角速度ベクトルをωとするとき、
   

   

   
   で表される式を用いて演算する
   請求項２に記載の流体解析装置。
4. 前記予め定めた方向と交差する方向は、相互に直交し、かつ各々前記予め定めた方向と直交する２方向である
   請求項３に記載の流体解析装置。
5. コンピュータが、
   物体と、前記物体の周囲に流れ、かつ粘度が空間的に変化する流体と、を含む空間を前記物体と前記流体との境界に非適合で、かつ直交する複数の格子線により複数の要素で分割して定義された数値計算対象モデルを用いて予め定めた方向に連続する前記境界を含む複数の要素各々における情報に基づいて、前記境界の予め定めた方向の流体の流れに関係する前記流体の速度勾配を示す第１の物理量を、前記複数の格子線によりなる格子の中心点又は重心点から前記境界までの距離を用いて演算し、
   演算された前記第１の物理量に基づいて、前記境界の前記予め定めた方向と交差する方向の流体の流れに関係する前記流体の速度勾配を示す第２の物理量を、前記第１の物理量と前記境界で成立する予め定めた速度勾配に関する関係式とを用いて演算し、
   演算された前記第１の物理量、及び演算された前記第２の物理量に基づいて、前記境界を含む予め定めた注目領域を流れる前記流体の挙動を導出する
   流体解析方法。
6. コンピュータを、請求項１から請求項４の何れか１項に記載された流体解析装置として機能させるための流体解析プログラム。
   

Images