説明

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

【課題】格子ボルツマン法による移動体周囲の流れ解析において、移動体の位置や姿勢が変更されても、計算効率の高い移動体周囲の流れ解析方法を提案する。
【解決手段】格子ボルツマン法による移動体周囲の流れ解析において、計算空間を分割する複数の格子系101,102を用い、移動体である回転体105を内包する格子系102では回転体105とともに動く座標系104を用い、格子系間の仮想粒子速度分布の座標変換に際しては、転送元格子系マクロ量の座標変換結果から転送先速度分布を構築する。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、流体解析方法、流体解析装置、及び流体解析プログラムに係り、特に、格子ボルツマン法による流体解析方法、流体解析装置、及び流体解析プログラムに関する。
【背景技術】
【0002】
格子ボルツマン法は、流体を仮想的な粒子の集合と捉え、計算格子上を動く仮想粒子の挙動を計算する流体解析方法であり、特に移動体周囲の流れ場を求める方法である。格子ボルツマン法に関しては、例えば、特表2010-500654号公報(特許文献1)に開示されたものがある。
【0003】
格子ボルツマン法による移動体周囲の流れ解析には、移動体を静止した格子の中で移動させる移動境界法が従来から用いられている。
【先行技術文献】
【特許文献】
【0004】
【特許文献1】特表2010-500654号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
移動体を静止した格子の中で移動させる移動境界法では、移動体の位置や姿勢が変更されるたびに、その変更を境界条件に反映させる幾何計算が必要となる。そのため、移動体が大きな表面積を有する場合、幾何計算の処理負荷が増大し計算効率が低下する。
【0006】
本発明は、計算効率の高い移動体周囲の流体解析方法、流体解析装置、及び流体解析プログラムを提供することを課題とする。
【課題を解決するための手段】
【0007】
本発明は、上記課題を解決するために、格子ボルツマン法による移動体周囲の流れ解析において、計算空間を分割する複数の格子系を用い、移動体を内包する格子系では移動体とともに動く座標系を用い、格子系間の仮想粒子速度分布の座標変換に際しては、転送元格子系マクロ量の座標変換結果から転送先速度分布を構築することを特徴とする。
【発明の効果】
【0008】
本発明によれば、移動物体が静止して見える動座標系を用いることで位置や姿勢の変更を境界条件に反映させる必要がなくなり、移動体周囲の流れを効率的に計算可能となる。
【図面の簡単な説明】
【0009】
【図1】本発明の一実施例にかかる流体解析方法を示すフローチャートである。
【図2】静止座標系の格子と回転座標系の格子および回転体の配置を示した図である。
【図3】接合部仮想粒子の回転座標系への変換処理を示したフローチャートである。
【図4】格子の接合部を示す図である。
【図5】本発明の一実施例で用いられる解析装置(計算処理装置)である。
【発明を実施するための形態】
【0010】
以下、図面を用いて本発明の一実施例を説明する。
【0011】
本実施例では、二次元計算空間中に配された移動体の周囲の流れを格子ボルツマン法によって流体解析する例を説明する。具体的には、図2に示すように二次元計算空間中に配された回転移動する物体である回転体105の周囲の流れを格子ボルツマン法よって計算する流体解析方法を説明する。
【0012】
尚、本実施例の流体解析には、例えば、図5に示すように、入力装置1、演算装置(コンピュータ)2、記憶装置3、出力装置4を備えた解析装置(計算処理装置)が用いられる。入力装置1から解析に必要なデータ(後述)が入力され、入力されたデータは記憶装置3に格納されている。また、記憶装置3には後述の流体解析方法を内容とする流体解析プログラムが格納されている。そして、演算装置2が記憶装置3に格納された流体解析プログラムを読み込み、必要なデータを記憶装置3から抽出し流体解析プログラムに従って演算を行う。演算装置2で得られた解析結果は記憶装置3に記憶され、また、出力装置4に出力される。
【0013】
初めに本実施例における流体解析の計算体系について説明する。
【0014】
本発明では計算空間を分割する複数の格子系を用いる。本実施例における流体解析では、図2に示したように、二次元計算空間を格子系101と回転体を内包する格子系102により分割する。格子系101には静止座標系103が定められ、格子系102には回転体105と一体となって角速度ωで回転する回転座標系(動座標系)104が定められている。即ち、移動体である回転体105を内包する格子系では回転体105とともに動く座標系を用いている。また、格子系101と格子系102はともに格子間隔Δxの直交格子とし、時刻tにおける静止座標系103と回転座標系104との傾きの大きさは角度θで表す。
【0015】
本実施例における流体解析では、格子系101の外周部境界107および回転体105の境界106に滑りなし壁面条件を課す。ただし、境界条件の種類は、本発明の流体解析方法と直接関係するものではない。即ち、他の境界条件を課しても良い。また、格子系101と格子系102の情報の受け渡しのために、格子系102と格子系102とが互いに重なる領域である接合部108を設けている。
【0016】
次に、図1を参照して、上述の計算体系に基づく流体解析を説明する。図1は本実施例における流体解析方法を示すフローチャートである。本実施例では、演算処理S1〜S11にしたがって流れ場を求めていく。
【0017】
初めのS1では、入力装置1を用いて解析開始にあたって必要なデータを入力する。主なデータは、回転体105の形状、初期の流れ場、境界条件の種類(ここでは、回転体境界106、外周部境界107に滑りなし壁面条件を設定)、解析パラメータ、および解析終了時刻である。
【0018】
次に、S2で格子系101の仮想粒子に対し演算処理を行う。同様にS3で格子系102の仮想粒子に対し演算処理を行う。S2はサブステップS4、S5、S6から構成され、S3はサブステップS7、S8、S9から構成される。
【0019】
サブステップS4では、静止座標系102での仮想格子の並進と衝突の過程を記述した次の式1で示される格子ボルツマン方程式から1タイムステップ(Δt)後の仮想粒子の速度分布fiの値を計算する。さらに、外周部境界107での速度分布fiを境界条件に基づいて更新する。
【0020】
【数1】

【0021】
ここで、xは空間位置で、ciは仮想粒子の並進方向であり、c=Δx/Δtとして以下の式2で定められる。
【0022】
【数2】

【0023】
また、fi(0)は平衡状態における速度分布であり以下の式3で定められる。
【0024】
【数3】

【0025】
ここで、ρは密度、uは静止座標系での速度である。ρ、uは速度モーメントの式4を用いて算出する。また重み係数wiは式5で与えられる。
【0026】
【数4】

【0027】
【数5】

【0028】
サブステップS5では、図3の手順に従って、静止座標系格子と回転座標系格子の接合部108の速度分布を回転座標系へ座標変換する。図3は接合部仮想粒子の回転座標系への変換処理のフローチャートを示す。格子系間の仮想粒子速度分布の座標変換に際しては、転送元格子系マクロ量の座標変換結果から転送先の仮想粒子速度分布を構築するようにしている。
【0029】
図3で示した演算処理は、接合部108における静止座標系マクロ量(密度ρ、速度u、運動量流束テンソルΠの非平衡成分Πneq(粘性応力に相当))の算出(手順201)、前記マクロ量の回転座標系104への変換(手順202)、回転座標系マクロ量から速度分布の平衡成分gi(0)の算出(手順203)、非平衡成分gi_neqの算出(手順204)、および回転座標系の速度分布関数giの算出(手順205)からなる。
【0030】
手順201における静止系マクロ量(ρ、u、Πneq)の算出においては、ρ、uに関しては前記式4で求めた結果を用い、Πneqを以下の式6で計算する。
【0031】
【数6】

【0032】
ここで、α1およびα2は静止座標系103の座標成分を示す。
【0033】
手順202では、ベクトル量uと二階の対称テンソル量Πneqを次の式7、式8にてそれぞれ回転座標系104の変数u´、Π´neqに変換する。なお、密度ρはスカラー量であり、座標変換に対して不変であるため、回転座標系104においても同じ値を用いる。
【0034】
【数7】

【0035】
【数8】

【0036】
手順203においては、ρ、u´を、前記式3と同様の平衡状態の式9に代入し、回転座標系での速度分布g_iの平衡成分gi(0)を算出する。
【0037】
【数9】

【0038】
手順204においては、クヌーセン数に相当する微少量εを用いて分布関数giを平衡分布gi(0)の周りに次式10のようにチャップマン・エンスコッグ展開したときに、1次の項fi(1)とρ、u´との間に一般に成り立つ関係式11、および式11から導かれる式12を利用する。
【0039】
【数10】

【0040】
【数11】

【0041】
【数12】

【0042】
ここで、β1およびβ2は回転座標系の座標成分を表す。δはデツタ関数である。また、csは音速であり、式2中のcとは「csの二乗=c/3」の関係がある。ci、wiは前記式2、式5で定められる。
【0043】
式11、式12において非平衡成分gi_neqとgi(1)が等しいと近似すると、式13のgi_neqを得る(即ち、実質的には、非平衡成分gi_neqを式11から算出)。
【0044】
【数13】

【0045】
手順205においては、前記手順203、204で求めた速度分布の平衡成分gi(0)および非平衡成分gi_neqを式14のように加算しgiを算出する。
【0046】
【数14】

【0047】
手順201〜205は、マクロ量を正確に再現する座標変換の構成となっている。即ち、算出したgiから再び、流体の密度、速度、および運動量流束テンソルを計算すると、回転座標系でのそれぞれの値と一致する。また密度と次式15の関係がある圧力、および密度と速度の積である運動量についても同様のことが言える。
【0048】
【数15】

【0049】
以上がサブステップS5において実行される座標変換処理の内容である。
【0050】
次に、サブステップS6においては、接合部108で求めたgiを格子系2の格子点上に補間する。例えば、図4における格子系102の格子点301に関しては、周囲にある格子点302などの格子系101の格子点から格子点301のgi値を線形補間する。
【0051】
以上がS2(サブステップS4、S5、S6)の説明であり、これによって格子系101の1タイムステップの処理が終了する。
【0052】
次に、格子系102の仮想粒子に対する演算処理S3(サブステップS7、S8、S9)を説明する。
【0053】
格子系102のダイナミクスは、回転体が静止しているように見える回転座標系104で記述される。そのため、サブステップS7、S8、S9には、格子系101と同様の演算S4〜S6がそれぞれ適用される。ただし、回転座標系であることから、前記式1の格子ボルツマン方程式が、遠心力とコリオリ力といった見かけの力Forceを加えた式16に置き換わる。また、静止座標系への座標変換、即ち、サブステップS8で実行される座標変換処理は式7、式8の逆変換となる。
【0054】
【数16】

【0055】
このように前記S3(サブステップS7、S8、S9)によって格子系102の仮想粒子の状態が1タイムステップ進んだ状態に更新される。また、接合部108において仮想粒子の情報交換が適切に行われた状態となる。そして、前記ステップS2とS3を、S10で解析終了時刻と判断されるまで反復実行することで、終了時刻までの流れ場を得ることができる。
【0056】
その後、S11において、解析結果を出力し解析終了の処理を実行して全体の処理を終える。以上の処理が図5に示す解析装置によって実行される。
【0057】
本実施例の流体解析方法と、従来の解析方法、即ち、物体を静止した格子系のなかで移動させる移動境界法との違いは、従来方法では物体移動に伴って移動する回転体境界106で幾何処理が生じるが、回転座標系104を用いる本実施例では、静止座標系格子と回転座標系格子の接合部108の処理が加わるものの、回転体境界106での幾何処理が不要なる点にある。この特徴より回転体表面積が大きな場合ほど従来方法に比べて計算負荷が低減される。
【0058】
なお、本実施例では、二次元計算空間の特定の格子モデルを対象としたが、三次元の格子モデルを用いる場合も式1〜16と同様の式が成り立ち、本実施例記載の流体解析方法を三次元計算空間へ容易に適用可能である。
【符号の説明】
【0059】
101 格子系(静止座標系格子)
102 格子系(回転座標系格子)
103 静止座標系
104 回転座標系(動座標系)
105 回転体(移動体)
106 回転体境界
107 静止座標系格子の外周部境界
108 静止座標系格子と回転座標系格子の接合部
201 静止座標系マクロ量取得処理
202 マクロ量変換処理
203 回転座標系速度分布の平衡成分算出処理
204 回転座標系速度分布の非平衡成分算出処理
205 回転座標系速度分布算出処理
301 回転座標系格子の格子点
302 静止座標系格子の格子点

【特許請求の範囲】
【請求項1】
格子ボルツマン法による移動体周囲の流れ解析を行う流体解析方法であって、計算空間を分割する複数の格子系を用い、前記複数の格子系のうち、移動体を内包する格子系では移動体とともに動く座標系を用い、前記複数の格子系間の仮想粒子速度分布の座標変換に際しては、転送元の格子系マクロ量の座標変換結果から転送先の仮想粒子速度分布を構築することを特徴とする流体解析方法。
【請求項2】
請求項1記載の流体解析方法において、前記転送元の格子系マクロ量の座標変換結果から転送先の仮想粒子速度分布を構築するステップは、転送元の座標系において仮想粒子の速度モーメントから密度ρ、速度u、運動量流束テンソルの非平衡成分Πneqを求め、次に、u, Πneqを転送先の座標系に座標変換してu´,Π´neqを得た後、転送先座標系での仮想粒子速度分布giの平衡成分を式9の平衡分布関数から算出し、非平衡成分を式11から算出するようにしたことを特徴とする流体解析方法。
【数1】

ここで、wiは重み係数、ciは仮想粒子の並進方向。
【数2】

ここで、csは音速、β1、β2は回転座標系の成分。
【請求項3】
格子ボルツマン法による移動体周囲の流れ解析を行う流体解析方法であって、
入力装置を用いて移動体の形状や境界条件の種類を含む解析条件を設定するステップと、
前記解析条件に基づき、コンピュータにおいて、静止座標系が定められた格子系の仮想粒子に対し格子ボルツマン法による演算処理を行うステップと、
前記解析条件に基づき、前記コンピュータにおいて、前記移動体とともに動く動座標系が定められ前記移動体を内包する格子系の仮想粒子に対し格子ボルツマン法による演算処理を行うステップと、
前記コンピュータにおいて、前記静止座標系が定められた格子系と前記動座標系が定められた格子系が重なる領域の仮想粒子速度分布を他の座標系への変換処理を行うステップと、
前記コンピュータにおいて、前記変換処理で得られた仮想粒子速度分布を前記他の座標系の格子系に反映するステップを有することを特徴とする流体解析方法。
【請求項4】
入力装置、記憶装置、演算装置、出力装置を有し、格子ボルツマン法による移動体周囲の流れ解析を行う流体解析装置であって、
前記記憶装置は、前記入力装置によって入力された移動体の形状や境界条件の種類を含む解析条件と、流体解析を前記演算装置に実行させるプログラムと、前記演算装置での演算結果を格納し、
前記演算装置は、前記プログラムに基づき、前記記憶装置から移動体の形状や境界条件の種類を含む解析条件を読み込み、前記解析条件に基づき、静止座標系が定められた格子系の仮想粒子に対し格子ボルツマン法による演算処理を行い、前記解析条件に基づき、前記移動体とともに動く動座標系が定められ前記移動体を内包する格子系の仮想粒子に対し格子ボルツマン法による演算処理を行い、前記静止座標系が定められた格子系と前記動座標系が定められた格子系が重なる領域の仮想粒子速度分布を他の座標系への変換処理を行い、前記変換処理で得られた仮想粒子速度分布を前記他の座標系の格子系に反映し、これらの処理を所定時間繰り返し行うことにより流体解析を行い、前記出力装置に流体解析結果を出力することを特徴とする流体解析装置。
【請求項5】
格子ボルツマン法による移動体周囲の流れ解析を行う流体解析方法をコンピュータに実行させるための流体解析プログラムにおいて、
移動体の形状や境界条件の種類を含む解析条件を読み込むステップと、
前記解析条件に基づき、静止座標系が定められた格子系の仮想粒子に対し格子ボルツマン法による演算処理を行うステップと、
前記解析条件に基づき、前記移動体とともに動く動座標系が定められ前記移動体を内包する格子系の仮想粒子に対し格子ボルツマン法による演算処理を行うステップと、
前記静止座標系が定められた格子系と前記動座標系が定められた格子系が重なる領域の仮想粒子速度分布を他の座標系への変換処理を行うステップと、
前記変換処理で得られた仮想粒子速度分布を前記他の座標系の格子系に反映するステップを有する流体解析方法をコンピュータに実行させるための流体解析プログラム。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate


【公開番号】特開2013−37434(P2013−37434A)
【公開日】平成25年2月21日(2013.2.21)
【国際特許分類】
【出願番号】特願2011−171041(P2011−171041)
【出願日】平成23年8月4日(2011.8.4)
【出願人】(000005108)株式会社日立製作所 (27,607)
【Fターム(参考)】