(57)【要約】【課題】格子ボルツマン法による移動体周囲の流れ解析において、移動体の位置や姿勢が変更されても、計算効率の高い移動体周囲の流れ解析方法を提案する。【解決手段】格子ボルツマン法による移動体周囲の流れ解析において、計算空間を分割する複数の格子系１０１，１０２を用い、移動体である回転体１０５を内包する格子系１０２では回転体１０５とともに動く座標系１０４を用い、格子系間の仮想粒子速度分布の座標変換に際しては、転送元格子系マクロ量の座標変換結果から転送先速度分布を構築する。【選択図】図２
2013037434.tif 000002
【特許請求の範囲】
【請求項１】
格子ボルツマン法による移動体周囲の流れ解析を行う流体解析方法であって、計算空間を分割する複数の格子系を用い、前記複数の格子系のうち、移動体を内包する格子系では移動体とともに動く座標系を用い、前記複数の格子系間の仮想粒子速度分布の座標変換に際しては、転送元の格子系マクロ量の座標変換結果から転送先の仮想粒子速度分布を構築することを特徴とする流体解析方法。
【請求項２】
請求項１記載の流体解析方法において、前記転送元の格子系マクロ量の座標変換結果から転送先の仮想粒子速度分布を構築するステップは、転送元の座標系において仮想粒子の速度モーメントから密度ρ、速度ｕ、運動量流束テンソルの非平衡成分Πneqを求め、次に、u, Πneqを転送先の座標系に座標変換してu´,Π´neqを得た後、転送先座標系での仮想粒子速度分布giの平衡成分を式９の平衡分布関数から算出し、非平衡成分を式１１から算出するようにしたことを特徴とする流体解析方法。
【数１】
2013037434.tif 000019
ここで、wiは重み係数、ciは仮想粒子の並進方向。
【数２】
2013037434.tif 000020
ここで、csは音速、β１、β２は回転座標系の成分。
【請求項３】
格子ボルツマン法による移動体周囲の流れ解析を行う流体解析方法であって、
入力装置を用いて移動体の形状や境界条件の種類を含む解析条件を設定するステップと、
前記解析条件に基づき、コンピュータにおいて、静止座標系が定められた格子系の仮想粒子に対し格子ボルツマン法による演算処理を行うステップと、
前記解析条件に基づき、前記コンピュータにおいて、前記移動体とともに動く動座標系が定められ前記移動体を内包する格子系の仮想粒子に対し格子ボルツマン法による演算処理を行うステップと、
前記コンピュータにおいて、前記静止座標系が定められた格子系と前記動座標系が定められた格子系が重なる領域の仮想粒子速度分布を他の座標系への変換処理を行うステップと、
前記コンピュータにおいて、前記変換処理で得られた仮想粒子速度分布を前記他の座標系の格子系に反映するステップを有することを特徴とする流体解析方法。
【請求項４】
入力装置、記憶装置、演算装置、出力装置を有し、格子ボルツマン法による移動体周囲の流れ解析を行う流体解析装置であって、
前記記憶装置は、前記入力装置によって入力された移動体の形状や境界条件の種類を含む解析条件と、流体解析を前記演算装置に実行させるプログラムと、前記演算装置での演算結果を格納し、
前記演算装置は、前記プログラムに基づき、前記記憶装置から移動体の形状や境界条件の種類を含む解析条件を読み込み、前記解析条件に基づき、静止座標系が定められた格子系の仮想粒子に対し格子ボルツマン法による演算処理を行い、前記解析条件に基づき、前記移動体とともに動く動座標系が定められ前記移動体を内包する格子系の仮想粒子に対し格子ボルツマン法による演算処理を行い、前記静止座標系が定められた格子系と前記動座標系が定められた格子系が重なる領域の仮想粒子速度分布を他の座標系への変換処理を行い、前記変換処理で得られた仮想粒子速度分布を前記他の座標系の格子系に反映し、これらの処理を所定時間繰り返し行うことにより流体解析を行い、前記出力装置に流体解析結果を出力することを特徴とする流体解析装置。
【請求項５】
格子ボルツマン法による移動体周囲の流れ解析を行う流体解析方法をコンピュータに実行させるための流体解析プログラムにおいて、
移動体の形状や境界条件の種類を含む解析条件を読み込むステップと、
前記解析条件に基づき、静止座標系が定められた格子系の仮想粒子に対し格子ボルツマン法による演算処理を行うステップと、
前記解析条件に基づき、前記移動体とともに動く動座標系が定められ前記移動体を内包する格子系の仮想粒子に対し格子ボルツマン法による演算処理を行うステップと、
前記静止座標系が定められた格子系と前記動座標系が定められた格子系が重なる領域の仮想粒子速度分布を他の座標系への変換処理を行うステップと、
前記変換処理で得られた仮想粒子速度分布を前記他の座標系の格子系に反映するステップを有する流体解析方法をコンピュータに実行させるための流体解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、流体解析方法、流体解析装置、及び流体解析プログラムに係り、特に、格子ボルツマン法による流体解析方法、流体解析装置、及び流体解析プログラムに関する。
【背景技術】
【０００２】
格子ボルツマン法は、流体を仮想的な粒子の集合と捉え、計算格子上を動く仮想粒子の挙動を計算する流体解析方法であり、特に移動体周囲の流れ場を求める方法である。格子ボルツマン法に関しては、例えば、特表2010-500654号公報（特許文献１）に開示されたものがある。
【０００３】
格子ボルツマン法による移動体周囲の流れ解析には、移動体を静止した格子の中で移動させる移動境界法が従来から用いられている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特表2010-500654号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
移動体を静止した格子の中で移動させる移動境界法では、移動体の位置や姿勢が変更されるたびに、その変更を境界条件に反映させる幾何計算が必要となる。そのため、移動体が大きな表面積を有する場合、幾何計算の処理負荷が増大し計算効率が低下する。
【０００６】
本発明は、計算効率の高い移動体周囲の流体解析方法、流体解析装置、及び流体解析プログラムを提供することを課題とする。
【課題を解決するための手段】
【０００７】
本発明は、上記課題を解決するために、格子ボルツマン法による移動体周囲の流れ解析において、計算空間を分割する複数の格子系を用い、移動体を内包する格子系では移動体とともに動く座標系を用い、格子系間の仮想粒子速度分布の座標変換に際しては、転送元格子系マクロ量の座標変換結果から転送先速度分布を構築することを特徴とする。
【発明の効果】
【０００８】
本発明によれば、移動物体が静止して見える動座標系を用いることで位置や姿勢の変更を境界条件に反映させる必要がなくなり、移動体周囲の流れを効率的に計算可能となる。
【図面の簡単な説明】
【０００９】
【図１】本発明の一実施例にかかる流体解析方法を示すフローチャートである。
【図２】静止座標系の格子と回転座標系の格子および回転体の配置を示した図である。
【図３】接合部仮想粒子の回転座標系への変換処理を示したフローチャートである。
【図４】格子の接合部を示す図である。
【図５】本発明の一実施例で用いられる解析装置（計算処理装置）である。
【発明を実施するための形態】
【００１０】
以下、図面を用いて本発明の一実施例を説明する。
【００１１】
本実施例では、二次元計算空間中に配された移動体の周囲の流れを格子ボルツマン法によって流体解析する例を説明する。具体的には、図２に示すように二次元計算空間中に配された回転移動する物体である回転体１０５の周囲の流れを格子ボルツマン法よって計算する流体解析方法を説明する。
【００１２】
尚、本実施例の流体解析には、例えば、図５に示すように、入力装置１、演算装置（コンピュータ）２、記憶装置３、出力装置４を備えた解析装置（計算処理装置）が用いられる。入力装置１から解析に必要なデータ（後述）が入力され、入力されたデータは記憶装置３に格納されている。また、記憶装置３には後述の流体解析方法を内容とする流体解析プログラムが格納されている。そして、演算装置２が記憶装置３に格納された流体解析プログラムを読み込み、必要なデータを記憶装置３から抽出し流体解析プログラムに従って演算を行う。演算装置２で得られた解析結果は記憶装置３に記憶され、また、出力装置４に出力される。
【００１３】
初めに本実施例における流体解析の計算体系について説明する。
【００１４】
本発明では計算空間を分割する複数の格子系を用いる。本実施例における流体解析では、図２に示したように、二次元計算空間を格子系１０１と回転体を内包する格子系１０２により分割する。格子系１０１には静止座標系１０３が定められ、格子系１０２には回転体１０５と一体となって角速度ωで回転する回転座標系（動座標系）１０４が定められている。即ち、移動体である回転体１０５を内包する格子系では回転体１０５とともに動く座標系を用いている。また、格子系１０１と格子系１０２はともに格子間隔Δxの直交格子とし、時刻tにおける静止座標系１０３と回転座標系１０４との傾きの大きさは角度θで表す。
【００１５】
本実施例における流体解析では、格子系１０１の外周部境界１０７および回転体１０５の境界１０６に滑りなし壁面条件を課す。ただし、境界条件の種類は、本発明の流体解析方法と直接関係するものではない。即ち、他の境界条件を課しても良い。また、格子系１０１と格子系１０２の情報の受け渡しのために、格子系１０２と格子系１０２とが互いに重なる領域である接合部１０８を設けている。
【００１６】
次に、図１を参照して、上述の計算体系に基づく流体解析を説明する。図１は本実施例における流体解析方法を示すフローチャートである。本実施例では、演算処理Ｓ１〜Ｓ１１にしたがって流れ場を求めていく。
【００１７】
初めのＳ１では、入力装置１を用いて解析開始にあたって必要なデータを入力する。主なデータは、回転体１０５の形状、初期の流れ場、境界条件の種類（ここでは、回転体境界１０６、外周部境界１０７に滑りなし壁面条件を設定）、解析パラメータ、および解析終了時刻である。
【００１８】
次に、Ｓ２で格子系１０１の仮想粒子に対し演算処理を行う。同様にＳ３で格子系１０２の仮想粒子に対し演算処理を行う。Ｓ２はサブステップＳ４、Ｓ５、Ｓ６から構成され、Ｓ３はサブステップＳ７、Ｓ８、Ｓ９から構成される。
【００１９】
サブステップＳ４では、静止座標系１０２での仮想格子の並進と衝突の過程を記述した次の式１で示される格子ボルツマン方程式から１タイムステップ(Δt)後の仮想粒子の速度分布fiの値を計算する。さらに、外周部境界１０７での速度分布fiを境界条件に基づいて更新する。
【００２０】
【数１】
2013037434.tif 000003
【００２１】
ここで、ｘは空間位置で、ciは仮想粒子の並進方向であり、c=Δx/Δtとして以下の式２で定められる。
【００２２】
【数２】
2013037434.tif 000004
【００２３】
また、fi(0)は平衡状態における速度分布であり以下の式３で定められる。
【００２４】
【数３】
2013037434.tif 000005
【００２５】
ここで、ρは密度、uは静止座標系での速度である。ρ、uは速度モーメントの式４を用いて算出する。また重み係数wiは式５で与えられる。
【００２６】
【数４】
2013037434.tif 000006
【００２７】
【数５】
2013037434.tif 000007
【００２８】
サブステップＳ５では、図３の手順に従って、静止座標系格子と回転座標系格子の接合部１０８の速度分布を回転座標系へ座標変換する。図３は接合部仮想粒子の回転座標系への変換処理のフローチャートを示す。格子系間の仮想粒子速度分布の座標変換に際しては、転送元格子系マクロ量の座標変換結果から転送先の仮想粒子速度分布を構築するようにしている。
【００２９】
図３で示した演算処理は、接合部１０８における静止座標系マクロ量（密度ρ、速度u、運動量流束テンソルΠの非平衡成分Πneq（粘性応力に相当））の算出（手順２０１）、前記マクロ量の回転座標系１０４への変換（手順２０２）、回転座標系マクロ量から速度分布の平衡成分gi(0)の算出（手順２０３）、非平衡成分gi_neqの算出（手順２０４）、および回転座標系の速度分布関数giの算出（手順２０５）からなる。
【００３０】
手順２０１における静止系マクロ量（ρ、u、Πneq）の算出においては、ρ、uに関しては前記式４で求めた結果を用い、Πneqを以下の式６で計算する。
【００３１】
【数６】
2013037434.tif 000008
【００３２】
ここで、α１およびα２は静止座標系１０３の座標成分を示す。
【００３３】
手順２０２では、ベクトル量uと二階の対称テンソル量Πneqを次の式７、式８にてそれぞれ回転座標系１０４の変数u´、Π´neqに変換する。なお、密度ρはスカラー量であり、座標変換に対して不変であるため、回転座標系１０４においても同じ値を用いる。
【００３４】
【数７】
2013037434.tif 000009
【００３５】
【数８】
2013037434.tif 000010
【００３６】
手順２０３においては、ρ、u´を、前記式３と同様の平衡状態の式９に代入し、回転座標系での速度分布g_iの平衡成分gi(0)を算出する。
【００３７】
【数９】
2013037434.tif 000011
【００３８】
手順２０４においては、クヌーセン数に相当する微少量εを用いて分布関数giを平衡分布gi(0)の周りに次式１０のようにチャップマン・エンスコッグ展開したときに、１次の項fi(1)とρ、u´との間に一般に成り立つ関係式１１、および式１１から導かれる式１２を利用する。
【００３９】
【数１０】
2013037434.tif 000012
【００４０】
【数１１】
2013037434.tif 000013
【００４１】
【数１２】
2013037434.tif 000014
【００４２】
ここで、β１およびβ２は回転座標系の座標成分を表す。δはデツタ関数である。また、csは音速であり、式２中のcとは「csの二乗=c/3」の関係がある。ci、wiは前記式２、式５で定められる。
【００４３】
式１１、式１２において非平衡成分gi_neqとgi(1)が等しいと近似すると、式１３のgi_neqを得る（即ち、実質的には、非平衡成分gi_neqを式１１から算出）。
【００４４】
【数１３】
2013037434.tif 000015
【００４５】
手順２０５においては、前記手順２０３、２０４で求めた速度分布の平衡成分gi(0)および非平衡成分gi_neqを式１４のように加算しgiを算出する。
【００４６】
【数１４】
2013037434.tif 000016
【００４７】
手順２０１〜２０５は、マクロ量を正確に再現する座標変換の構成となっている。即ち、算出したgiから再び、流体の密度、速度、および運動量流束テンソルを計算すると、回転座標系でのそれぞれの値と一致する。また密度と次式１５の関係がある圧力、および密度と速度の積である運動量についても同様のことが言える。
【００４８】
【数１５】
2013037434.tif 000017
【００４９】
以上がサブステップＳ５において実行される座標変換処理の内容である。
【００５０】
次に、サブステップＳ６においては、接合部１０８で求めたgiを格子系２の格子点上に補間する。例えば、図４における格子系１０２の格子点３０１に関しては、周囲にある格子点３０２などの格子系１０１の格子点から格子点３０１のgi値を線形補間する。
【００５１】
以上がＳ２（サブステップＳ４、Ｓ５、Ｓ６）の説明であり、これによって格子系１０１の１タイムステップの処理が終了する。
【００５２】
次に、格子系１０２の仮想粒子に対する演算処理Ｓ３（サブステップＳ７、Ｓ８、Ｓ９）を説明する。
【００５３】
格子系１０２のダイナミクスは、回転体が静止しているように見える回転座標系１０４で記述される。そのため、サブステップＳ７、Ｓ８、Ｓ９には、格子系１０１と同様の演算Ｓ４〜Ｓ６がそれぞれ適用される。ただし、回転座標系であることから、前記式１の格子ボルツマン方程式が、遠心力とコリオリ力といった見かけの力Forceを加えた式１６に置き換わる。また、静止座標系への座標変換、即ち、サブステップＳ８で実行される座標変換処理は式７、式８の逆変換となる。
【００５４】
【数１６】
2013037434.tif 000018
【００５５】
このように前記Ｓ３（サブステップＳ７、Ｓ８、Ｓ９）によって格子系１０２の仮想粒子の状態が１タイムステップ進んだ状態に更新される。また、接合部１０８において仮想粒子の情報交換が適切に行われた状態となる。そして、前記ステップＳ２とＳ３を、Ｓ１０で解析終了時刻と判断されるまで反復実行することで、終了時刻までの流れ場を得ることができる。
【００５６】
その後、Ｓ１１において、解析結果を出力し解析終了の処理を実行して全体の処理を終える。以上の処理が図５に示す解析装置によって実行される。
【００５７】
本実施例の流体解析方法と、従来の解析方法、即ち、物体を静止した格子系のなかで移動させる移動境界法との違いは、従来方法では物体移動に伴って移動する回転体境界１０６で幾何処理が生じるが、回転座標系１０４を用いる本実施例では、静止座標系格子と回転座標系格子の接合部１０８の処理が加わるものの、回転体境界１０６での幾何処理が不要なる点にある。この特徴より回転体表面積が大きな場合ほど従来方法に比べて計算負荷が低減される。
【００５８】
なお、本実施例では、二次元計算空間の特定の格子モデルを対象としたが、三次元の格子モデルを用いる場合も式１〜１６と同様の式が成り立ち、本実施例記載の流体解析方法を三次元計算空間へ容易に適用可能である。
【符号の説明】
【００５９】
１０１  格子系（静止座標系格子）
１０２  格子系（回転座標系格子）
１０３  静止座標系
１０４  回転座標系（動座標系）
１０５  回転体（移動体）
１０６  回転体境界
１０７  静止座標系格子の外周部境界
１０８  静止座標系格子と回転座標系格子の接合部
２０１  静止座標系マクロ量取得処理
２０２  マクロ量変換処理
２０３  回転座標系速度分布の平衡成分算出処理
２０４  回転座標系速度分布の非平衡成分算出処理
２０５  回転座標系速度分布算出処理
３０１  回転座標系格子の格子点
３０２  静止座標系格子の格子点
【図１】
2013037434.tif 000021
【図２】
2013037434.tif 000022
【図３】
2013037434.tif 000023
【図４】
2013037434.tif 000024
【図５】
2013037434.tif 000025