JP5209298B2 - 流れのシミュレーション計算方法およびシステム - Google Patents

Description

本発明は、物体の周りの流体の流れのシミュレーション計算方法およびシステムに関する。

物体の周りの流体の流れを、数値流体力学（Computational Fluid Dynamics, CFD）によって解析する方法は、格子の定め方により、構造格子法、非構造格子法および直交格子法に分類される（たとえば、特許文献１の第２段落、図１乃至３）。

表１は、構造格子法、非構造格子法および直交格子法の特徴をまとめた表である。

構造格子法は、高分解能方式の適用に非常に優れているが、格子（メッシュ）生成が煩雑である。また、複雑な形状への適合性や、解適応グリッド制御（格子の形状を解に適応させる方法）の適用にも制約がある。非構造格子法は、解適応グリッド制御の適用に非常に優れており、複雑な形状への適合性も良好である。しかし、高分解能方式の適用や格子生成に制約がある。直交格子法は、格子生成に非常に優れているが、複雑な形状への適合性や、解適応グリッド制御の適用に難がある。また、高分解能方式の適用に制約がある。表１は、上記の３方法の大規模並列演算および大規模可視化への適用性についても示している。

このように、複雑な形状への適合性、解適応グリッド制御の適用、高分解能方式の適用および格子生成に優れた、物体の周りの流体の流れのシミュレーション計算方法は存在しなかった。
特開２００５−７８４１６号公報

したがって、複雑な形状への適合性、解適応グリッド制御の適用、高分解能方式の適用および格子生成に優れた、物体の周りの流体の流れのシミュレーション計算方法およびシステムに対するニーズがある。

本発明による物体の周りの流体の流れのシミュレーション計算方法は、対象領域を、複数の計算単位領域であるキューブに分割し、それぞれのキューブに、等しい数の直交格子要素であるセルを生成し、各計算ステップにおいて、複数のキューブ内の計算を行い、各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する方法である。 前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定値以下であるようにしながら、所望の計算精度が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返して、前記物体の形状に合わせてキューブのサイズを適切に定める。

本発明による物体の周りの流体の流れのシミュレーション計算を行うシステムは、対象領域を、複数の計算単位領域であるキューブに分割し、それぞれのキューブに、等しい数の直交格子要素であるセルを生成する格子生成部と、各計算ステップにおいて、複数のキューブ内の計算を並列に行い、各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する演算部と、シミュレーション計算の制御を行う計算制御部と、を備える。前記格子生成部は、前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定値以下であるようにしながら、所望の計算精度が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返して、前記物体の形状に合わせてキューブのサイズを適切に定める。

本発明においては、物体の形状に応じてサイズを定めることができるキューブを使用するので、複雑な形状への適合性があり、解適応グリッド制御や高分解能スキームにも適している。また、それぞれのキューブには等しい数の直交格子要素であるセルを備えるので、メッシュ生成は容易であり、それぞれのキューブの計算負荷はほぼ等しいので、並列計算における負荷分散が容易である。

本発明の基本的なアイディアは、複雑な形状に対するナビエ・ストークス方程式の時間精度の高い計算用の、高速かつ低記憶容量の方法を開発する際に、単純な直交格子法を適用することである。

本発明の方法においては、計算領域は、複数の３次元サブ領域に分割される。このサブ領域をキューブと呼称する。各々のキューブのサイズは、物体の形状および局所的な特徴的流れ長さに適応するように定められる。各々のキューブにおいて均一に分割された直交格子が生成される。直交格子の要素をセルと呼称する。全てのキューブは、同数のセルを含み、その結果、局所的な計算の分解能は、キューブのサイズによって定まる。全てのキューブの間で格子密度が同じであるので、キューブの計算の負荷はほぼ等しく並列計算が簡単になる。

図３（ａ）は、物体の周りに生成されたキューブを示す図である。図３（ｂ）は、図３（ａ）の部分拡大図である。

図１は、本発明の一実施形態によるシミュレーション計算方法を示す流れ図である。

ステップＳ０１０１０において、キューブが生成される。キューブの生成方法については、後で説明する。

ステップＳ０１０２０において、各々のキューブの流体計算が行われる。キューブの流体計算は、たとえば複数のプロセッサによって並列に行われる。流体計算については後で説明する。

ステップＳ０１０３０において、所定のキューブと該所定のキューブに隣接するキューブとの間で、情報の交換が行われる。キューブ間の情報の交換は、流体計算の１ステップごとに行われる。

ステップＳ０１０４０において、流体計算の所定数の計算ステップが行われたかどうか判断される。流体計算の所定数の計算ステップが行われていれば、ステップＳ０１０５０に進む。流体計算の所定数の計算ステップが行われていなければ、ステップＳ０１０２０に戻り、ステップＳ０１０２０およびステップＳ０１０３０を実施する。

ステップＳ０１０５０において、キューブのサイズの調整が必要かどうか判断される。キューブのサイズの調整が必要であれば、ステップＳ０１０６０に進む。キューブのサイズの調整が必要でなければ、ステップＳ０１０７０に進む。

ステップＳ０１０６０において、キューブのサイズが調整される。キューブのサイズ調整方法については後で説明する。

ステップＳ０１０７０において、シミュレーションを終了する時点であるかどうか判断する。シミュレーションを終了する時点であれば、シミュレーションを終了する。シミュレーションを終了する時点でなければ、ステップＳ０１０２０に戻る。

キューブの生成方法
本実施形態において、格子生成は２段階で行われる。第１の段階は、図３に示されるように、流れの場を満たす種々のサイズのキューブを生成することである。第２の段階は、物体の近くに位置するキューブ内にセルを要素とする直交格子を生成することである。

図２は、本実施形態におけるキューブの生成方法を示す流れ図である。図２の流れ図は、図１のステップＳ０１０１０におけるキューブの生成方法の詳細を示すものである。図２の方法は、各方向に均一のメッシュ間隔の直交格子を備えた、流れの場全体を覆う単一のキューブから開始する。物体境界と交差するキューブは、最小サイズが指定された空間分解能に到達するまで、任意の等しい子キューブに繰り返し分割される。

ステップＳ０２０１０において、物体面とキューブとの交差が判定される。

ステップＳ０２０２０において、物体面と交差するキューブは、均等に分割される。本実施形態において、物体面と交差するキューブは、等しい８（＝２^３）個の子キューブに分割される。本実施形態において、キューブの分割の状態は、ツリー構造によって表現される。

図５は、キューブの分割の状態を表現するツリー構造を示す図である。流れの場全体を覆う単一のキューブは、ルート・キューブとして表現される。ルート・キューブのレベルは１であり、分割が１回行われるごとにレベル数が１増加する。所定のキューブを複数のキューブに分割した場合に、該所定のキューブを親キューブ、該複数のキューブを子キューブと呼称する。親キューブは、その子キューブへのポインタを有する。このように、キューブ生成は、分割のレベルによって示され、最下位のレベルのキューブはリーフ・キューブと呼称される。リーフ・キューブのみが流れ計算に使用される。ツリー構造において、子キューブへのポインタを追いかけることによって、ルート・キューブからリーフ・キューブを追跡することができる。メッシュ適応のために、キューブは、最下位のレベルにおいてのみ追加されまたは削除される。

ステップＳ０２０３０において、隣接するキューブのサイズが調整される。隣接するキューブの間のレベル差は、キューブ間のデータ転送エラーを最小とするために１に制限される。したがって、隣接するキューブに比較して２レベルの差を有するキューブは、検出され分割される。

ステップＳ０２０４０において、物体内部のキューブは除去される。流体計算に必要ないからである。

ステップＳ０２０５０において、物体境界と交差するキューブのサイズは、指定された分解能に到達しているかどうか判断される。物体境界と交差するキューブのサイズが、指定された分解能に到達していれば、ステップＳ０２０６０に進む。物体境界と交差するキューブのサイズが、指定された分解能に到達していなければ、ステップＳ０２０１０に戻り分割が繰り返される。

ステップＳ０２０６０において、非等方キューブが生成される。

図６は、非等方に分割されたキューブを示す図である。キューブは、１方向または２方向に分割され、２（＝２^１）個または４（＝２^２）個の子キューブがそれぞれ生成される。キューブが非等方にＮ方向に分割されると、子キューブの数は２^Ｎである。このように、ツリーは、２^Ｎ−ツリーとなる。特定の方向に繰り返し分割されると、キューブは大きなアスペクト比を有する。図３（ａ）は、そのような場合を示している。

ステップＳ０２０７０において、リーフ・キューブに直交格子が生成される。本実施形態において、リーフ・キューブの直交格子は、一例として、１６ｘ１６ｘ１６個の均一サイズのセルからなる。リーフ・キューブ内のセルの個数は、自由に定めることができる。ただし、それぞれのリーフ・キューブ内のセルの個数は等しくなるようにする。

リーフ・キューブは、後で説明する流体計算の単位である。それぞれのリーフ・キューブ内のセルの個数は等しいので、それぞれのリーフ・キューブの計算の負荷はほぼ等しい。

ステップＳ０２０８０において、セルの物体情報が構築される。すなわち、各々のセルは、物体の内側に位置するか外側に位置するかマークを付される。たとえば、物体の内側のセルは０とされ、外側のセルは１とされる。物体の境界は、異なるマークを有する隣接セルの間のセル境界によって定められる。

図４は、物体（固体）と流体との境界を含むキューブおよび該キューブ内のセルを示す図である。図４（ａ）は、物体の周囲のキューブを示す図である。図４（ａ）は、キューブ内のセルを表示していない。図４（ｂ）は、図４（ａ）において太い枠で示したキューブを示す図である。図４（ｂ）は、キューブ内のセルを表示している。図４（ｂ）において、物体の内側のセルには濃い陰影を付し、物体の外側のセルには陰影を付していない。セル境界には薄い陰影を付している。

本実施形態において、キューブ内のセルの数は、１６ｘ１６ｘ１６個であり、大きいので、キューブ内のセルの物体情報を圧縮するのが好ましい。本実施形態においては、ラン・レングス法によってキューブ内のセルの物体情報を圧縮する。

図７は、２次元のラン・レングス・データ・ラインを説明するための図である。図７に示すように、キューブ内の全てのセルは、０または１のフラグを有する。したがって、最小の情報は、セルがフラグの値を変える場所である。図７に示すように、（ｊ，ｋ）によって定義される２次元の場における位置は、以下の式による１次元配列ｉによって記述される。
ここで、ｎ_ｊは、ｊ方向のキューブ内セル数であり、Ｘ、Ｙは整数として、ｍｏｄ（Ｘ，Ｙ）は、Ｘ／Ｙの余りである。

上記の式とともに、１次元アレイJK_runlength(i)が、以下のように定義される。

１次元アレイをデコードすることにより、以下が得られる。
ここで、INTEGER(Z)は数Ｚの整数部である。

図８は、３次元のラン・レングス・データ・ラインを説明するための図である。図８に示すように、（ｊ，ｋ，ｌ）によって定義される３次元の場に対しては、以下のとおりである。

上記のラン・レングス法によって、キューブ内のセルの物体情報は、大幅に圧縮される。

ステップＳ０２０９０において、キューブの境界情報が構築される。本実施形態において、各々のキューブは、以下のデータを有する。
１）親キューブおよび子キューブへのポインタ
２）キューブの系図
３）キューブの分割レベル
４）セルについての物体情報
４）のセルについての物体情報は、既に説明した。１）乃至３）の情報がキューブの境界情報に対応する。キューブの境界情報は、図１のステップＳ０１０３０において、所定のキューブと該キューブに隣接するキューブとの間で、情報の交換を行う際に、隣接するキューブを見つけるために使用される。

図９は、キューブの系図の定義を示すための図である。キューブの系図は、隣接するキューブをすばやくかつ簡単に検索するのに利用される。キューブの位置は、親キューブ内の位置にしたがって、それぞれの方向に０または１によって示される。キューブの系図は、ルート・キューブからのそのファミリの全てのキューブを含むテーブルである。

図１０は、キューブ位置と系図との関係を示す図である。黒のキューブの位置は、それぞれの方向において２進ツリーで表現される。２進ツリー内の斜線で示した領域は、ルート・キューブに対応する。系図自体は、ルート・キューブにおけるキューブ位置である。このように、隣接キューブの系図は、黒いキューブからの相対的な位置によって得られる。ツリーを遡ることなく、一方向に横切って、系図によって隣接キューブを見出すことができる。

ステップＳ０２１００において、物体と流体との境界を含むキューブにおいてサブ・グリッドが生成される。サブ・グリッドの生成については後で説明する。

ステップＳ０２１１０において、キューブの計算の負荷を複数のプロセッサに分散させる。全てのキューブは、同数のセルを含むので、全てのキューブの計算の負荷はほぼ等しい。したがって、たとえば、複数の同一性能のプロセッサに、同一数のキューブを分配することにより、均等な負荷分散を行うことができる。ただし、キューブの計算の負荷を１台のプロセッサで実行してもよい。

サブ・グリッドの生成方法
図１１は、本実施形態におけるサブ・グリッドの生成方法を示す流れ図である。図１１の流れ図は、図２のステップＳ０２１００におけるサブ・グリッドの生成方法の詳細を示すものである。

ステップＳ１１０１０において、リーフ・キューブのインターフェース・セルが抽出される。

図１２は、固体（物体）と流体との境界を含むキューブを示す図である。図１２において、固体中のセルは黒い丸で示し、流体中のセルは四角で示す。固体に隣接するセルは、インターフェース・セルと呼称され、ドットを付した丸で示す。

ステップＳ１１０２０において、物体面近傍にサブ・グリッドが生成される。

図１３は、サブ・グリッドが生成されたキューブを示す図である。図１３において、固体中のセルは黒い丸で示し、流体中のセルは大きな四角で示し、インターフェース・セルは、ドットを付した丸で示す。サブ・グリッド・ポイントは白い丸で示し、境界ポイントは小さな四角で示す。

サブ・グリッドは、直交格子数を増加させずに、壁面近傍の計算精度を維持するように、直交格子と独立に付加される新しい計算ノードである。サブ・グリッドは、インターフェース・セルの中心点の位置から流体境界面に向かって、当該流体境界面に垂直方向へ段階的に近づけた複数の位置に設けられる。この手順は、サブ・グリッドが壁境界に到達し、サブ・グリッド内の各層が壁境界から一定の距離となるまで続けられる。

ステップＳ１１０３０において、サブ・グリッドの層数が十分であるかどうか判断される。層数が十分であれば、ステップＳ１１０４０に進む。層数が十分でなければ、ステップＳ１１０２０に戻り、再びサブ・グリッドを生成する。

ステップＳ１１０４０において、クラウドが生成される。クラウドについては後で説明する。

ステップＳ１１０５０において、形状関数が構築される。形状関数については後で説明する。

流体計算
流体計算について以下に説明する。

圧縮性粘性流れを支配する無次元のナビエ・ストークス方程式は、保存形で以下のように記述される。
ここで、t、xおよびReは、それぞれ、時間、座標およびレイノルズ数である。i、jおよびkは、上記および以下の方程式において、座標インデックスを表し、アインシュタインの総和則がそれぞれのインデックスに適用される。保存変数ベクトルＵ、非粘性フラックス・ベクトルＦ(Ｕ) および粘性フラックス・ベクトルＧ(Ｕ)は、
によって定義される。ここで、ρ, ｐ, ｅ，ＴおよびＫは、密度、圧力、全エネルギ、温度および流体の熱伝導率であり、ｕ_iは、座標方向における流れの速度成分を表す。この方程式のセットは、完全気体の状態方程式によって閉じられる。
ここでγは、比熱の比率である。粘性応力テンソルの成分τ_ijは、
によって与えられる。ここでδ_ijは、クロネッカーのデルタ記号であり、μは、サザランドの関係によって定まる粘性係数である。

乱流シミュレーションに対して、粘性μは、たとえば、スパラート・アラマラス一方程式モデル(Spalart and Allmaras one-equation model)によって評価される、渦粘性によって置き換えられる。

可変密度または一定密度の低マッハ数流れ問題の解法における効率を改善するために、前処理が行われる。前処理された支配方程式を導出するために、式（９）によって表現された元の系は、保存変数Ｕから原始変数Ｑへ以下のように変換される。
ここで、原始変数ベクトルＱおよびヤコビアン
は、
および
によって与えられる。ここで、ＨおよびＣ_ｐは、それぞれ、定圧における全エンタルピーおよび比熱である。元のヤコビアン行列
を前処理行列Γによって置き換えることにより、前処理ナビエ・ストークス方程式が得られる。
本実施形態において、前処理行列Гは以下のように記述される。
ここで、
である。ここで、Ｕ_ｒは、基準速度であり、
によって与えられる。ここで、｜ｖ｜は、ローカル速度の大きさであり、ｃは、音速である。また、基準速度は、局所的な拡散速度よりも小さくないはずである。そこで、
である。ここで、Δｄは、特性セル長さである。

前処理ナビエ・ストークス方程式の固有値は以下のとおりである。
ここで、
である。ここで、ｎは、セル境界の外向きの法線ベクトルを示す。基準速度が局所速度と同じ次数である場合には、全ての固有値は、ｕと同じオーダーを有する。基準速度が音の局所速度と等しければ、前処理系は、非前処理ナビエ・ストークス方程式の元の系へ戻る。

圧縮性ナビエ・ストークス方程式は、直交格子における、セルを中心とする、有限体積の方式によって離散化される。式（１６）は、以下の代数形式で記述される。
ここで、Ｖ_ｉは、セルｉの体積である。総和インデックスｊ(ｉ)は、セルｉの全ての面を意味し、ΔＳ_ijは、セルｉおよびｊの間のインターフェース領域である。ｈ(ｎ_ij)は、検査体積境界に垂直な、非粘性数値フラックス・ベクトルである。
は、検査体積境界の両側のプリミティブ値である。数値フラックスｈは、たとえば、Harten-Lax-van Leer-Einfeldt-Wada (HLLEW, Obayashi, S., and Guruswamy, G. P., “Convergence Acceleration of a Navier-Stokes Solver for Efficient Static Aeroelastic Computations,” AIAA Journal, Vol. 33, No. 6, 1995, pp. 1134, 1141.)の近似リーマン解法を使用して計算される。二重のセル境界におけるプリミティブ変数は、３次精度以上を達成するために、たとえば、MUSCLアプローチ(Yamamoto, S., and Daiguji, H., “Higher-Order-Accurate Upwind Schemes for Solving the Compressible Euler and Navier-Stokes Equations,” Computer and Fluids, Vol. 22, No. 2/3, 1933, pp. 259, 270.)によって評価される。粘性フラックス項Ｇは、２次精度中心差分によって評価される。

本実施形態においては、時間積分に、たとえば、Lower-Upper Symmetric Gauss-Seidel (LU-SGS)インプリシット法が使用される（Sharov, D., and Nakahashi, K., “Reordering of Hybrid Unstructured Grids for Lower-Upper Symmetric Gauss-Seidel Computations,” AIAA Journal, No. 36, 1998, pp. 484, 486.）。この近似要素分解法の最も魅力的な利点は、インプリシット演算子にある近似を行うことによって、LU-SGS法の元の公式において固有のヤコビアン行列の評価および記憶を完全に除去することができることである。結果としてのLU-SGS法は、時間ごとのイクスプリシット法よりも安価である。

LU-SGS法は、式（２７）における非粘性フラックスの最初の総和に対して、ノード・ポイントｊ(ｉ)を下三角行列のグループｊ∈Ｌ(ｉ)と上三角行列のグループｊ∈Ｕ(ｉ)とに分けるテクニックに基づいている。時間レベルｎの関係
において、式（２７）は、対称Gauss-Seidel方に含まれる２個のステップ
および
にしたがって解かれる。ここで、Ｒ_ｉは、式（２７）の右手側の残差である。行列Ｄは、
によって与えられる。ここで、Δｔは時間ステップを示し、Ｉは単位行列である。Jameson およびTurkelによって近似されたフラックス・ヤコビアンＪ^±を使用することによって、Ｄ対角化形を以下のようにする。
ここで、λ_ijは、空間半径
である。LU-SGSは、結局２回の掃き出し(sweep)によって解かれる。第１は、低位の（順方向の）掃き出し
であり、第２は、上位の（逆方向の）掃き出し
である。ここで、
である。

境界条件を課すために、インターフェース・セルｉには体積力が導入される。式（２）の元のナビエ・ストークス方程式は以下のように記述される。
ＲＨＳ_ｉは、非粘性フラックスおよび粘性フラックスを含み、
によって与えられる。体積力Ｓ_ｉは、
によって示される。ここで、(・)_IBは、ＩＢ法(Immersed Boundary Method)によって課された体積力である。

式（３９）の第２項は、グリッドレス(gridless)法によって評価される。該方法は、クラウドおよびノード位置の形状関数の概念を必要とする。各々のインターフェース・セルのクラウドは、隣接流体セル、インターフェース・セル、境界ポイントおよびサブ・グリッドを含む。インターフェース・セルｉにおける流れ変数は、そのクラウド内の周囲のポイントの流れ変数および形状関数の和によって評価される。K番目の座標方向における任意の特定の関数 の空間導関数は、形状関数の導関数φによって以下のように記述される。
ここで、Ｎは、考察対象のクラウド内のノードの数、添え字ｎは、クラウドに属するポイントのインデックスを示す。式（４０）によって、インターフェース・セルとクラウドのメンバー・ポイントの間の全ての中間点において、非粘性フラックスおよび粘性フラックスが計算される。

キューブのサイズ調整方法
図１４は、本実施形態におけるキューブのサイズ調整方法を示す流れ図である。図１４の流れ図は、図１のステップＳ０１０５０およびステップＳ０１０６０におけるキューブのサイズ調整方法の詳細を示すものである。

ステップＳ１４０１０において、全てのキューブの適応パラメータが計算される。適応パラメータは、たとえば、画像処理において輪郭を捉えるために使用されるラプラシアン・フィルタに基づく。流れ解ベクトルの選択されたコンポーネントに対するラプラシアン・フィルタは、以下の式によって与えられる。コンポーネントは、たとえば、マッハ数であってもよい。
ここで、Δlは、キューブ・サイズであり、Ｌは、長さスケールである。Ｌの値を増加させると、大きなキューブが細分化される。しかし、増加させすぎると、重要な流れの特性を解くために、十分なキューブのクラスタリングが行われることなく、滑らかな流れ領域において多すぎる数のキューブが生成される。他方、Ｌの値を減少させると、小さなキューブが細分化される。しかし、減少させすぎると、キューブのサイズの分布における空間的な変化が大きくなり、解の精度を劣化させる。本実施形態においては、Ｌの適切な値が、０と１の間で経験的に定められる。

キューブ内の全てのセルのラプラシアンの絶対値の総和を、該キューブの適応パラメータとする。

ステップＳ１４０２０において、全てのキューブの適応パラメータの平均値および標準偏差が計算される。

ステップＳ１４０３０において、全てのキューブの適応パラメータの平均値および標準偏差に基づいて、適応パラメータの細分化のしきい値Ghighおよび粗大化のしきい値Glowが定められる。

ステップＳ１４０４０において、キューブの適応パラメータが細分化のしきい値よりも大きいかどうか判断される。キューブの適応パラメータが細分化のしきい値よりも大きければ、ステップＳ１４０５０に進む。キューブの適応パラメータが細分化のしきい値以下であれば、ステップＳ１４０６０に進む。

ステップＳ１４０５０において、キューブが細分化される。

ステップＳ１４０６０において、キューブの適応パラメータが粗大化のしきい値よりも小さいかどうか判断される。キューブの適応パラメータが粗大化のしきい値よりも小さければ、ステップＳ１４０７０に進む。キューブの適応パラメータが粗大化のしきい値以上であれば、ステップＳ１４０８０に進む。

ステップＳ１４０７０において、キューブが粗大化される。

ステップＳ１４０８０において、隣接するキューブのサイズ調整が行われる。

ステップＳ１４０９０において、キューブの境界情報が構築される。

ステップＳ１４１００において、キューブの計算の負荷を複数のプロセッサに分散させる。

システム構成
図１５は、本発明の一実施形態によるシミュレーション計算システムの構成を示す図である。

シミュレーション計算システムは、計算制御部１０１、格子生成部１０３、入出力部１０５、システム・バス１０７および演算部１０９を含む。演算部１０９は、ｎ台のプロセッサ１０９１乃至１０９ｎを含む。

格子生成部１０３は、キューブの生成やキューブのサイズ調整を行う。演算部１０９は、各キューブの流体計算およびキューブ間情報交換を行う。計算制御部１０１は、複数のプロセッサへのキューブの計算の分散、キューブのサイズ調整の開始時期の判断、計算終了時期の判断などを行う。入出力部１０５は、シミュレーション対象のデータの入力およびシミュレーション結果の出力を行う。

図１６Ａは、本実施形態によるシミュレーション計算システムの性能を示す図である。

図１６Ｂは、図１６Ａに示したシミュレーション計算システムの性能を測定するために、その周りにおいてシミュレーション計算が行われる物体の形状を示す図である。物体はアメッド・ボディ(Ahrmed body)と呼称されるものである。グリッド数は、４百万である。

図１６Ａ（ａ）の横軸および縦軸は、それぞれ、プロセッサ（ＣＰＵ）の台数および該台数のシステムの性能を示す。図１６Ａ（ｂ）の横軸および縦軸は、それぞれ、プロセッサ（ＣＰＵ）の台数および該台数のシステムの性能比率を示す。複数のプロセッサを含むシステムは、プロセッサの台数に比例した性能を備えるのが理想である。ここで、プロセッサの台数に比例した性能に対する実際の性能の比率を性能比率と定義する。理想的には、複数のプロセッサを含むシステムの性能比率は１である。しかし、実際には、プロセッサの台数が増加するにしたがって性能比率は大きく低下する。図１６Ａ（ｂ）に示すように、本実施形態によるシステムにおいては、プロセッサの台数が６４台となっても性能比率は０．８（８０％）以上である。その理由は、本実施形態によるシステムにおいては、計算負荷のほぼ等しいキューブによって、計算負荷をプロセッサ間にほぼ均等に分散できるからである。

図１７乃至２０は、本実施形態のシミュレーション計算方法による計算結果を示す図である。

図１７および図１８は、翼の周りのマッハ数および圧力Ｐを示す図である。流入マッハ数は、０．７３である。レイノルズ数は、６．５６ｘ１０^６であり、全流れ領域に対して、乱流を仮定している。迎角は、２．７９度に設定される。衝撃波は、翼の上面において、先端から翼弦の約５５％に形成される。

図１９は、面圧の計算値と実験値との比較を示す図である。図１９の縦軸は圧力係数を示す。図１９の横軸は位置座標を示す。翼の前縁の座標が０、翼の後縁の座標が１である。２つの線のうち上方の線の値が翼の上面の圧力係数値であり、下方の線の値が翼の下面の圧力係数値である。計算値は、全体として実験値によく一致している。

図２０は、計算によって求められた、車両モデル近くの流線の軌跡を示す図である。図２０は、流線における圧力Ｐも示す。

上述のように本発明のシミュレーション計算方法は、複雑な形状への適合性、解適応グリッド制御の適用、高分解能方式の適用および格子生成に優れている。

本発明の実施形態の特徴は以下のとおりである。

本発明の一実施形態において、キューブの分割の回数をツリー構造のレベルと対応させて、キューブの集合を該ツリー構造の最も下位のレベルの要素であるリーフによって表現し、該ツリー構造を使用して、所定のキューブに隣接するキューブを見出して、境界情報を交換するように構成している。

本実施形態によれば、ツリー構造を使用することによって、所定のキューブに隣接するキューブを容易に見出すことができる。

本発明の実施形態において、隣接するキューブは、ツリー構造において同じレベルか隣接レベルであるように構成している。

本実施形態によれば、隣接するキューブは、ツリー構造において同じレベルか隣接レベルであるように構成することによって、隣接するキューブのサイズの比が所定値以内となり、スムージングが行われる。

本発明の実施形態において、キューブが、前記ツリー構造に関し、親および子へのポインタと、系図と、レベルとのデータを有し、該キューブのセルに関し、前記物体の有無のデータを有するように構成している。

本実施形態によれば、キューブの有するデータの量が少なくなる。

本発明の実施形態において、前記物体の有無のデータを、ラン・レングス・データとして記憶するように構成している。

本実施形態によれば、キューブの有するデータの量がさらに少なくなる。

本発明の実施形態において、前記分割が、Ｎを１、２または３として２^Ｎ個への均等分割である。

本実施形態によれば、Ｎを変えることにより、キューブの等方分割または非等方分割を行うことができる。

本発明の実施形態において、前記対象領域が３次元の領域であり、キューブが立方体または直方体である。

本発明の実施形態において、所定数の計算ステップを実行した後に、キューブの適応パラメータに基づいて、該キューブのサイズが適切であるかどうか評価し、適切でない場合には該キューブのサイズを変更するように構成している。

本実施形態によれば、所定数の計算ステップを実行した後にキューブの適応パラメータによって、該キューブのサイズを評価することによってキューブのサイズを調整することができる。その結果、シミュレーション計算の精度が向上する。

本発明の実施形態において、前記物体と前記流体との境界を含むキューブ内の前記境界の近傍に、前記境界の近傍の分解能を向上させるように、セルとは別に、セルよりも細かな間隔でサブ・グリッドを設けるように構成している。

本実施形態によれば、セルの数を一定値以内としながら、物体と流体との境界近傍の分解能を向上させることができる。

本発明の実施形態において、計算の負荷を、単一のキューブの計算を単位として複数のプロセッサに分散させるように構成している。

複数のキューブは、サイズは異なっても等しい数のセルを備えるので、複数のキューブの計算の負荷はほぼ等しい。したがって、単一のキューブの計算を単位として、所望の大きさの負荷を複数のプロセッサに分散させることができる。

本発明の一実施形態によるシミュレーション計算方法を示す流れ図である。 本実施形態におけるキューブの生成方法を示す流れ図である。 物体の周りに生成されたキューブを示す図である。 物体（固体）と流体との境界を含むキューブ及び該キューブ内のセルを示す図である。 キューブの分割の状態を表現するツリー構造を示す図である。 非等方に分割されたキューブを示す図である。 ２次元のラン・レングス・データ・ラインを説明するための図である。 ３次元のラン・レングス・データ・ラインを説明するための図である。 キューブの系図の定義を示すための図である。 キューブ位置と系図との関係を示す図である。 本実施形態におけるサブ・グリッドの生成方法を示す流れ図である。 固体（物体）と流体との境界を含むキューブを示す図である。 サブ・グリッドが生成されたインターフェース・セルを示す図である。 本実施形態におけるキューブのサイズ調整方法を示す流れ図である。 本発明の一実施形態によるシミュレーション計算システムの構成を示す図である。 本実施形態によるシミュレーション計算システムの性能を示す図である。 図１６Ａに示したシミュレーション計算システムの性能を測定するために、その周りにおいてシミュレーション計算が行われる物体の形状を示す図である。 翼の周りのマッハ数を示す図である。 翼の周りの圧力を示す図である。 面圧の計算値と実験値との比較を示す図である。 計算によって求められた、車両近くの流線の軌跡を示す図である。

符号の説明

１０１…計算制御部、１０３…格子生成部、１０５…入出力部、１０９…演算部

Claims (10)

1. コンピュータにより実行される、流体力学に基づく物体の周りの流体の流れのシミュレーション計算方法であり、
   対象領域を、キューブと呼称される複数の計算単位領域に分割し、
   それぞれのキューブに、等しい数のセルと呼称される直交格子要素を生成し、
   各計算ステップにおいて、複数のキューブ内について流体力学方程式に基づく数値計算を行い、
   各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する、
   方法であって、
   前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定の範囲であるようにしながら、所望の分解能が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返して、前記物体の形状に合わせてキューブのサイズを定め、かつ、
   前記物体と前記流体との境界面を含むキューブ内の前記境界面の近傍に、前記境界面の近傍の計算精度を向上させるように、前記セルとは別に、前記セルよりも細かな間隔でサブ・グリッドを設けるように構成され、
   前記サブ・グリッドは、前記境界面の近傍の計算精度を向上させるため前記セルとは独立に付加される計算ノードであって、前記境界面に隣接する前記セルであるインターフェース・セルの中心点の位置から当該境界面に向かって、当該境界面に垂直な方向へ段階的に近づけた複数の位置に設けられ、
   前記インターフェース・セルの中心点における流れは、当該中心点における流れの形状関数の空間導関数を、当該インターフェース・セルに隣接する流体中のセルの中心点、前記境界面上の点、及びサブ・グリッドを含む各点の周囲のポイントにおける流れの形状関数の線形結合として表現することにより求められる、
   シミュレーション計算方法。
2. キューブの分割の回数をツリー構造のレベルと対応させて、キューブの集合を該ツリー構造の最も下位のレベルの要素であるリーフによって表現し、該ツリー構造を使用して、所定のキューブに隣接するキューブを見出して、該所定のキューブと該所定のキューブに隣接するキューブとの間で情報を交換するように構成した請求項１に記載のシミュレーション計算方法。
3. 隣接するキューブは、ツリー構造において同じレベルか隣接レベルであるように構成した請求項２に記載のシミュレーション計算方法。
4. キューブが、前記ツリー構造に関し、親および子へのポインタと、系図と、レベルとのデータを有し、該キューブのセルに関し、前記物体の有無のデータを有するように構成した請求項２に記載のシミュレーション計算方法。
5. 前記物体の有無のデータを、ラン・レングス・データとして記憶するように構成した請求項４に記載のシミュレーション計算方法。
6. 分割の際にキューブが、Ｎを１、２または３として２Ｎ個の均等サイズのキューブに分割される請求項１に記載のシミュレーション計算方法。
7. 前記対象領域が３次元の領域であり、キューブが立方体または直方体である請求項１に記載のシミュレーション計算方法。
8. 所定数の計算ステップを実行した後に、キューブの適応パラメータに基づいて、該キューブのサイズが適切であるかどうか評価し、適切でない場合には該キューブのサイズを変更するように構成した請求項１に記載のシミュレーション計算方法。
9. 計算の負荷を、単一のキューブの計算を単位として複数のプロセッサに分散させるように構成した請求項１に記載のシミュレーション計算方法。
10. 流体力学に基づく物体の周りの流体の流れのシミュレーション計算を行うシステムであって、
    対象領域を、キューブと呼称される複数の計算単位領域に分割し、それぞれのキューブに、等しい数のセルと呼称される直交格子要素を生成する格子生成部と、
    各計算ステップにおいて、複数のキューブ内について流体力学方程式に基づく数値計算を並列に行い、各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する演算部と、
    シミュレーション計算の制御を行う計算制御部と、を備え、
    前記格子生成部が、前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定値以下であるようにしながら、所望の計算精度が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返して、前記物体の形状に合わせてキューブのサイズを定め、かつ、
    前記物体と前記流体との境界面を含むキューブ内の前記境界面の近傍に、前記境界面の近傍の計算精度を向上させるように、前記セルとは別に、前記セルよりも細かな間隔でサブ・グリッドを設け、
    前記サブ・グリッドは、前記境界面の近傍の計算精度を向上させるため前記セルとは独立に付加される計算ノードであって、前記境界面に隣接する前記セルであるインターフェース・セルの中心点の位置から当該境界面に向かって、当該境界面に垂直な方向へ段階的に近づけた複数の位置に設けられ、
    前記インターフェース・セルの中心点における流れは、当該中心点における流れの形状関数の空間導関数を、当該インターフェース・セルに隣接する流体中のセルの中心点、前記境界面上の点、及びサブ・グリッドを含む各点の周囲のポイントにおける流れの形状関数の線形結合として表現することにより求められる、
    物体の周りの流体の流れのシミュレーション計算を行うシステム。