(57)【要約】【課題】複雑な形状への適合性、解適応グリッド制御の適用、高分解能方式の適用および格子生成に優れた、物体の周りの流体の流れのシミュレーション計算方法を提供する。【解決手段】本発明による計算方法は、対象領域をキューブと呼称される計算単位領域に分割し、それぞれのキューブに等しい数の直交格子要素であるセルを生成し、各計算ステップにおいて、複数のキューブ内の計算を並列に行い、各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する。前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定値以下であるようにしながら、所望の計算精度が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返す。その結果、前記物体の形状に合わせてキューブのサイズが適切に定められる。【選択図】図４
2008165804.tif 000002
【特許請求の範囲】
【請求項１】
物体の周りの流体の流れのシミュレーション計算方法であり、
対象領域を、キューブと呼称される複数の計算単位領域に分割し、
それぞれのキューブに、等しい数のセルと呼称される直交格子要素を生成し、
各計算ステップにおいて、複数のキューブ内の計算を行い、
各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する方法であって、
前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定の範囲であるようにしながら、所望の分解能が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返して、前記物体の形状に合わせてキューブのサイズを適切に定めるシミュレーション計算方法。
【請求項２】
キューブの分割の回数をツリー構造のレベルと対応させて、キューブの集合を該ツリー構造の最も下位のレベルの要素であるリーフによって表現し、該ツリー構造を使用して、所定のキューブに隣接するキューブを見出して、該所定のキューブと該所定のキューブに隣接するキューブとの間で情報を交換するように構成した請求項１に記載のシミュレーション計算方法。
【請求項３】
隣接するキューブは、ツリー構造において同じレベルか隣接レベルであるように構成した請求項２に記載のシミュレーション計算方法。
【請求項４】
キューブが、前記ツリー構造に関し、親および子へのポインタと、系図と、レベルとのデータを有し、該キューブのセルに関し、前記物体の有無のデータを有するように構成した請求項２に記載のシミュレーション計算方法。
【請求項５】
前記物体の有無のデータを、ラン・レングス・データとして記憶するように構成した請求項４に記載のシミュレーション計算方法。
【請求項６】
分割の際にキューブが、Ｎを１、２または３として２Ｎ個の均等サイズのキューブに分割される請求項１に記載のシミュレーション計算方法。
【請求項７】
前記対象領域が３次元の領域であり、キューブが立方体または直方体である請求項１に記載のシミュレーション計算方法。
【請求項８】
所定数の計算ステップを実行した後に、キューブの適応パラメータに基づいて、該キューブのサイズが適切であるかどうか評価し、適切でない場合には該キューブのサイズを変更するように構成した請求項１に記載のシミュレーション計算方法。
【請求項９】
前記物体と前記流体との境界を含むキューブ内の前記境界の近傍に、前記境界の近傍の計算精度を向上させるように、セルとは別に、セルよりも細かな間隔でサブ・グリッドを設けるように構成した請求項１に記載のシミュレーション計算方法。
【請求項１０】
計算の負荷を、単一のキューブの計算を単位として複数のプロセッサに分散させるように構成した請求項１に記載のシミュレーション計算方法。
【請求項１１】
物体の周りの流体の流れのシミュレーション計算を行うシステムであって、
対象領域を、キューブと呼称される複数の計算単位領域に分割し、それぞれのキューブに、等しい数のセルと呼称される直交格子要素を生成する格子生成部と、
各計算ステップにおいて、複数のキューブ内の計算を並列に行い、各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する演算部と、
シミュレーション計算の制御を行う計算制御部と、を備え、
前記格子生成部が、前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定値以下であるようにしながら、所望の計算精度が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返して、前記物体の形状に合わせてキューブのサイズを適切に定める物体の周りの流体の流れのシミュレーション計算を行うシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物体の周りの流体の流れのシミュレーション計算方法およびシステムに関する。
【背景技術】
【０００２】
物体の周りの流体の流れを、数値流体力学（Computational Fluid Dynamics, CFD）によって解析する方法は、格子の定め方により、構造格子法、非構造格子法および直交格子法に分類される（たとえば、特許文献１の第２段落、図１乃至３）。
【０００３】
表１は、構造格子法、非構造格子法および直交格子法の特徴をまとめた表である。
【表１】
2008165804.tif 000003
【０００４】
構造格子法は、高分解能方式の適用に非常に優れているが、格子（メッシュ）生成が煩雑である。また、複雑な形状への適合性や、解適応グリッド制御（格子の形状を解に適応させる方法）の適用にも制約がある。非構造格子法は、解適応グリッド制御の適用に非常に優れており、複雑な形状への適合性も良好である。しかし、高分解能方式の適用や格子生成に制約がある。直交格子法は、格子生成に非常に優れているが、複雑な形状への適合性や、解適応グリッド制御の適用に難がある。また、高分解能方式の適用に制約がある。表１は、上記の３方法の大規模並列演算および大規模可視化への適用性についても示している。
【０００５】
このように、複雑な形状への適合性、解適応グリッド制御の適用、高分解能方式の適用および格子生成に優れた、物体の周りの流体の流れのシミュレーション計算方法は存在しなかった。
【特許文献１】特開２００５−７８４１６号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
したがって、複雑な形状への適合性、解適応グリッド制御の適用、高分解能方式の適用および格子生成に優れた、物体の周りの流体の流れのシミュレーション計算方法およびシステムに対するニーズがある。
【課題を解決するための手段】
【０００７】
本発明による物体の周りの流体の流れのシミュレーション計算方法は、対象領域を、複数の計算単位領域であるキューブに分割し、それぞれのキューブに、等しい数の直交格子要素であるセルを生成し、各計算ステップにおいて、複数のキューブ内の計算を行い、各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する方法である。  前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定値以下であるようにしながら、所望の計算精度が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返して、前記物体の形状に合わせてキューブのサイズを適切に定める。
【０００８】
本発明による物体の周りの流体の流れのシミュレーション計算を行うシステムは、対象領域を、複数の計算単位領域であるキューブに分割し、それぞれのキューブに、等しい数の直交格子要素であるセルを生成する格子生成部と、各計算ステップにおいて、複数のキューブ内の計算を並列に行い、各計算ステップ終了後に、隣接するキューブの間で境界情報を交換する演算部と、シミュレーション計算の制御を行う計算制御部と、を備える。前記格子生成部は、前記対象領域をキューブに分割する際に、隣接するキューブのサイズの比が所定値以下であるようにしながら、所望の計算精度が得られるように、前記物体と前記流体との境界を含むキューブが十分に小さくなるまで分割を繰り返して、前記物体の形状に合わせてキューブのサイズを適切に定める。
【０００９】
本発明においては、物体の形状に応じてサイズを定めることができるキューブを使用するので、複雑な形状への適合性があり、解適応グリッド制御や高分解能スキームにも適している。また、それぞれのキューブには等しい数の直交格子要素であるセルを備えるので、メッシュ生成は容易であり、それぞれのキューブの計算負荷はほぼ等しいので、並列計算における負荷分散が容易である。
【発明を実施するための最良の形態】
【００１０】
本発明の基本的なアイディアは、複雑な形状に対するナビエ・ストークス方程式の時間精度の高い計算用の、高速かつ低記憶容量の方法を開発する際に、単純な直交格子法を適用することである。
【００１１】
本発明の方法においては、計算領域は、複数の３次元サブ領域に分割される。このサブ領域をキューブと呼称する。各々のキューブのサイズは、物体の形状および局所的な特徴的流れ長さに適応するように定められる。各々のキューブにおいて均一に分割された直交格子が生成される。直交格子の要素をセルと呼称する。全てのキューブは、同数のセルを含み、その結果、局所的な計算の分解能は、キューブのサイズによって定まる。全てのキューブの間で格子密度が同じであるので、キューブの計算の負荷はほぼ等しく並列計算が簡単になる。
【００１２】
図３（ａ）は、物体の周りに生成されたキューブを示す図である。図３（ｂ）は、図３（ａ）の部分拡大図である。
【００１３】
図１は、本発明の一実施形態によるシミュレーション計算方法を示す流れ図である。
【００１４】
ステップＳ０１０１０において、キューブが生成される。キューブの生成方法については、後で説明する。
【００１５】
ステップＳ０１０２０において、各々のキューブの流体計算が行われる。キューブの流体計算は、たとえば複数のプロセッサによって並列に行われる。流体計算については後で説明する。
【００１６】
ステップＳ０１０３０において、所定のキューブと該所定のキューブに隣接するキューブとの間で、情報の交換が行われる。キューブ間の情報の交換は、流体計算の１ステップごとに行われる。
【００１７】
ステップＳ０１０４０において、流体計算の所定数の計算ステップが行われたかどうか判断される。流体計算の所定数の計算ステップが行われていれば、ステップＳ０１０５０に進む。流体計算の所定数の計算ステップが行われていなければ、ステップＳ０１０２０に戻り、ステップＳ０１０２０およびステップＳ０１０３０を実施する。
【００１８】
ステップＳ０１０５０において、キューブのサイズの調整が必要かどうか判断される。キューブのサイズの調整が必要であれば、ステップＳ０１０６０に進む。キューブのサイズの調整が必要でなければ、ステップＳ０１０７０に進む。
【００１９】
ステップＳ０１０６０において、キューブのサイズが調整される。キューブのサイズ調整方法については後で説明する。
【００２０】
ステップＳ０１０７０において、シミュレーションを終了する時点であるかどうか判断する。シミュレーションを終了する時点であれば、シミュレーションを終了する。シミュレーションを終了する時点でなければ、ステップＳ０１０２０に戻る。
【００２１】
キューブの生成方法
本実施形態において、格子生成は２段階で行われる。第１の段階は、図３に示されるように、流れの場を満たす種々のサイズのキューブを生成することである。第２の段階は、物体の近くに位置するキューブ内にセルを要素とする直交格子を生成することである。
【００２２】
図２は、本実施形態におけるキューブの生成方法を示す流れ図である。図２の流れ図は、図１のステップＳ０１０１０におけるキューブの生成方法の詳細を示すものである。図２の方法は、各方向に均一のメッシュ間隔の直交格子を備えた、流れの場全体を覆う単一のキューブから開始する。物体境界と交差するキューブは、最小サイズが指定された空間分解能に到達するまで、任意の等しい子キューブに繰り返し分割される。
【００２３】
ステップＳ０２０１０において、物体面とキューブとの交差が判定される。
【００２４】
ステップＳ０２０２０において、物体面と交差するキューブは、均等に分割される。本実施形態において、物体面と交差するキューブは、等しい８（＝２３）個の子キューブに分割される。本実施形態において、キューブの分割の状態は、ツリー構造によって表現される。
【００２５】
図５は、キューブの分割の状態を表現するツリー構造を示す図である。流れの場全体を覆う単一のキューブは、ルート・キューブとして表現される。ルート・キューブのレベルは１であり、分割が１回行われるごとにレベル数が１増加する。所定のキューブを複数のキューブに分割した場合に、該所定のキューブを親キューブ、該複数のキューブを子キューブと呼称する。親キューブは、その子キューブへのポインタを有する。このように、キューブ生成は、分割のレベルによって示され、最下位のレベルのキューブはリーフ・キューブと呼称される。リーフ・キューブのみが流れ計算に使用される。ツリー構造において、子キューブへのポインタを追いかけることによって、ルート・キューブからリーフ・キューブを追跡することができる。メッシュ適応のために、キューブは、最下位のレベルにおいてのみ追加されまたは削除される。
【００２６】
ステップＳ０２０３０において、隣接するキューブのサイズが調整される。隣接するキューブの間のレベル差は、キューブ間のデータ転送エラーを最小とするために１に制限される。したがって、隣接するキューブに比較して２レベルの差を有するキューブは、検出され分割される。
【００２７】
ステップＳ０２０４０において、物体内部のキューブは除去される。流体計算に必要ないからである。
【００２８】
ステップＳ０２０５０において、物体境界と交差するキューブのサイズは、指定された分解能に到達しているかどうか判断される。物体境界と交差するキューブのサイズが、指定された分解能に到達していれば、ステップＳ０２０６０に進む。物体境界と交差するキューブのサイズが、指定された分解能に到達していなければ、ステップＳ０２０１０に戻り分割が繰り返される。
【００２９】
ステップＳ０２０６０において、非等方キューブが生成される。
【００３０】
図６は、非等方に分割されたキューブを示す図である。キューブは、１方向または２方向に分割され、２（＝２１）個または４（＝２２）個の子キューブがそれぞれ生成される。キューブが非等方にＮ方向に分割されると、子キューブの数は２Ｎである。このように、ツリーは、２Ｎ−ツリーとなる。特定の方向に繰り返し分割されると、キューブは大きなアスペクト比を有する。図３（ａ）は、そのような場合を示している。
【００３１】
ステップＳ０２０７０において、リーフ・キューブに直交格子が生成される。本実施形態において、リーフ・キューブの直交格子は、一例として、１６ｘ１６ｘ１６個の均一サイズのセルからなる。リーフ・キューブ内のセルの個数は、自由に定めることができる。ただし、それぞれのリーフ・キューブ内のセルの個数は等しくなるようにする。
【００３２】
リーフ・キューブは、後で説明する流体計算の単位である。それぞれのリーフ・キューブ内のセルの個数は等しいので、それぞれのリーフ・キューブの計算の負荷はほぼ等しい。
【００３３】
ステップＳ０２０８０において、セルの物体情報が構築される。すなわち、各々のセルは、物体の内側に位置するか外側に位置するかマークを付される。たとえば、物体の内側のセルは０とされ、外側のセルは１とされる。物体の境界は、異なるマークを有する隣接セルの間のセル境界によって定められる。
【００３４】
図４は、物体（固体）と流体との境界を含むキューブおよび該キューブ内のセルを示す図である。図４（ａ）は、物体の周囲のキューブを示す図である。図４（ａ）は、キューブ内のセルを表示していない。図４（ｂ）は、図４（ａ）において太い枠で示したキューブを示す図である。図４（ｂ）は、キューブ内のセルを表示している。図４（ｂ）において、物体の内側のセルには濃い陰影を付し、物体の外側のセルには陰影を付していない。セル境界には薄い陰影を付している。
【００３５】
本実施形態において、キューブ内のセルの数は、１６ｘ１６ｘ１６個であり、大きいので、キューブ内のセルの物体情報を圧縮するのが好ましい。本実施形態においては、ラン・レングス法によってキューブ内のセルの物体情報を圧縮する。
【００３６】
図７は、２次元のラン・レングス・データ・ラインを説明するための図である。図７に示すように、キューブ内の全てのセルは、０または１のフラグを有する。したがって、最小の情報は、セルがフラグの値を変える場所である。図７に示すように、（ｊ，ｋ）によって定義される２次元の場における位置は、以下の式による１次元配列ｉによって記述される。
【数１】
2008165804.tif 000004
ここで、ｎｊは、ｊ方向のキューブ内セル数であり、Ｘ、Ｙは整数として、ｍｏｄ（Ｘ，Ｙ）は、Ｘ／Ｙの余りである。
【００３７】
上記の式とともに、１次元アレイJK_runlength(i)が、以下のように定義される。
【数２】
2008165804.tif 000005
【００３８】
１次元アレイをデコードすることにより、以下が得られる。
【数３】
2008165804.tif 000006
ここで、INTEGER(Z)は数Ｚの整数部である。
【００３９】
図８は、３次元のラン・レングス・データ・ラインを説明するための図である。図８に示すように、（ｊ，ｋ，ｌ）によって定義される３次元の場に対しては、以下のとおりである。
【数４】
2008165804.tif 000007
【００４０】
上記のラン・レングス法によって、キューブ内のセルの物体情報は、大幅に圧縮される。
【００４１】
ステップＳ０２０９０において、キューブの境界情報が構築される。本実施形態において、各々のキューブは、以下のデータを有する。
１）親キューブおよび子キューブへのポインタ
２）キューブの系図
３）キューブの分割レベル
４）セルについての物体情報
４）のセルについての物体情報は、既に説明した。１）乃至３）の情報がキューブの境界情報に対応する。キューブの境界情報は、図１のステップＳ０１０３０において、所定のキューブと該キューブに隣接するキューブとの間で、情報の交換を行う際に、隣接するキューブを見つけるために使用される。
【００４２】
図９は、キューブの系図の定義を示すための図である。キューブの系図は、隣接するキューブをすばやくかつ簡単に検索するのに利用される。キューブの位置は、親キューブ内の位置にしたがって、それぞれの方向に０または１によって示される。キューブの系図は、ルート・キューブからのそのファミリの全てのキューブを含むテーブルである。
【００４３】
図１０は、キューブ位置と系図との関係を示す図である。黒のキューブの位置は、それぞれの方向において２進ツリーで表現される。２進ツリー内の斜線で示した領域は、ルート・キューブに対応する。系図自体は、ルート・キューブにおけるキューブ位置である。このように、隣接キューブの系図は、黒いキューブからの相対的な位置によって得られる。ツリーを遡ることなく、一方向に横切って、系図によって隣接キューブを見出すことができる。
【００４４】
ステップＳ０２１００において、物体と流体との境界を含むキューブにおいてサブ・グリッドが生成される。サブ・グリッドの生成については後で説明する。
【００４５】
ステップＳ０２１１０において、キューブの計算の負荷を複数のプロセッサに分散させる。全てのキューブは、同数のセルを含むので、全てのキューブの計算の負荷はほぼ等しい。したがって、たとえば、複数の同一性能のプロセッサに、同一数のキューブを分配することにより、均等な負荷分散を行うことができる。ただし、キューブの計算の負荷を１台のプロセッサで実行してもよい。
【００４６】
サブ・グリッドの生成方法
図１１は、本実施形態におけるサブ・グリッドの生成方法を示す流れ図である。図１１の流れ図は、図２のステップＳ０２１００におけるサブ・グリッドの生成方法の詳細を示すものである。
【００４７】
ステップＳ１１０１０において、リーフ・キューブのインターフェース・セルが抽出される。
【００４８】
図１２は、固体（物体）と流体との境界を含むキューブを示す図である。図１２において、固体中のセルは黒い丸で示し、流体中のセルは四角で示す。固体に隣接するセルは、インターフェース・セルと呼称され、ドットを付した丸で示す。
【００４９】
ステップＳ１１０２０において、物体面近傍にサブ・グリッドが生成される。
【００５０】
図１３は、サブ・グリッドが生成されたキューブを示す図である。図１３において、固体中のセルは黒い丸で示し、流体中のセルは大きな四角で示し、インターフェース・セルは、ドットを付した丸で示す。サブ・グリッド・ポイントは白い丸で示し、境界ポイントは小さな四角で示す。
【００５１】
サブ・グリッドは、直交格子数を増加させずに、壁面近傍の計算精度を維持するように、直交格子と独立に付加される新しい計算ノードである。サブ・グリッドは、インターフェース・セルのポイントを、流体境界面に垂直方向に数段階伸ばすことによって定められる。この手順は、サブ・グリッドが壁境界に到達し、サブ・グリッド内の各層が壁境界から一定の距離となるまで続けられる。
【００５２】
ステップＳ１１０３０において、サブ・グリッドの層数が十分であるかどうか判断される。層数が十分であれば、ステップＳ１１０４０に進む。層数が十分でなければ、ステップＳ１１０２０に戻り、再びサブ・グリッドを生成する。
【００５３】
ステップＳ１１０４０において、クラウドが生成される。クラウドについては後で説明する。
【００５４】
ステップＳ１１０５０において、形状関数が構築される。形状関数については後で説明する。
【００５５】
流体計算
流体計算について以下に説明する。
【００５６】
圧縮性粘性流れを支配する無次元のナビエ・ストークス方程式は、保存形で以下のように記述される。
【数５】
2008165804.tif 000008
ここで、t、xおよびReは、それぞれ、時間、座標およびレイノルズ数である。i、jおよびkは、上記および以下の方程式において、座標インデックスを表し、アインシュタインの総和則がそれぞれのインデックスに適用される。保存変数ベクトルＵ、非粘性フラックス・ベクトルＦ(Ｕ) および粘性フラックス・ベクトルＧ(Ｕ)は、
【数６】
2008165804.tif 000009
によって定義される。ここで、ρ, ｐ, ｅ，ＴおよびＫは、密度、圧力、全エネルギ、温度および流体の熱伝導率であり、ｕiは、座標方向における流れの速度成分を表す。この方程式のセットは、完全気体の状態方程式によって閉じられる。
【数７】
2008165804.tif 000010
ここでγは、比熱の比率である。粘性応力テンソルの成分τijは、
【数８】
2008165804.tif 000011
によって与えられる。ここでδijは、クロネッカーのデルタ記号であり、μは、サザランドの関係によって定まる粘性係数である。
【００５７】
乱流シミュレーションに対して、粘性μは、たとえば、スパラート・アラマラス一方程式モデル(Spalart and Allmaras one-equation model)によって評価される、渦粘性によって置き換えられる。
【００５８】
可変密度または一定密度の低マッハ数流れ問題の解法における効率を改善するために、前処理が行われる。前処理された支配方程式を導出するために、式（９）によって表現された元の系は、保存変数Ｕから原始変数Ｑへ以下のように変換される。
【数９】
2008165804.tif 000012
ここで、原始変数ベクトルＱおよびヤコビアン
【数１０】
2008165804.tif 000013
は、
【数１１】
2008165804.tif 000014
および
【数１２】
2008165804.tif 000015
によって与えられる。ここで、ＨおよびＣｐは、それぞれ、定圧における全エンタルピーおよび比熱である。元のヤコビアン行列
【数１３】
2008165804.tif 000016
を前処理行列Γによって置き換えることにより、前処理ナビエ・ストークス方程式が得られる。
【数１４】
2008165804.tif 000017
本実施形態において、前処理行列Гは以下のように記述される。
【数１５】
2008165804.tif 000018
ここで、
【数１６】
2008165804.tif 000019
である。ここで、Ｕｒは、基準速度であり、
【数１７】
2008165804.tif 000020
によって与えられる。ここで、｜ｖ｜は、ローカル速度の大きさであり、ｃは、音速である。また、基準速度は、局所的な拡散速度よりも小さくないはずである。そこで、
【数１８】
2008165804.tif 000021
である。ここで、Δｄは、特性セル長さである。
【００５９】
前処理ナビエ・ストークス方程式の固有値は以下のとおりである。
【数１９】
2008165804.tif 000022
ここで、
【数２０】
2008165804.tif 000023
である。ここで、ｎは、セル境界の外向きの法線ベクトルを示す。基準速度が局所速度と同じ次数である場合には、全ての固有値は、ｕと同じオーダーを有する。基準速度が音の局所速度と等しければ、前処理系は、非前処理ナビエ・ストークス方程式の元の系へ戻る。
【００６０】
圧縮性ナビエ・ストークス方程式は、直交格子における、セルを中心とする、有限体積の方式によって離散化される。式（１６）は、以下の代数形式で記述される。
【数２１】
2008165804.tif 000024
ここで、Ｖｉは、セルｉの体積である。総和インデックスｊ(ｉ)は、セルｉの全ての面を意味し、ΔＳijは、セルｉおよびｊの間のインターフェース領域である。ｈ(ｎij)は、検査体積境界に垂直な、非粘性数値フラックス・ベクトルである。
【数２２】
2008165804.tif 000025
は、検査体積境界の両側のプリミティブ値である。数値フラックスｈは、たとえば、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法は、式（２７）における非粘性フラックスの最初の総和に対して、ノード・ポイントｊ(ｉ)を下三角行列のグループｊ∈Ｌ(ｉ)と上三角行列のグループｊ∈Ｕ(ｉ)とに分けるテクニックに基づいている。時間レベルｎの関係
【数２３】
2008165804.tif 000026
において、式（２７）は、対称Gauss-Seidel方に含まれる２個のステップ
【数２４】
2008165804.tif 000027
および
【数２５】
2008165804.tif 000028
にしたがって解かれる。ここで、Ｒｉは、式（２７）の右手側の残差である。行列Ｄは、
【数２６】
2008165804.tif 000029
によって与えられる。ここで、Δｔは時間ステップを示し、Ｉは単位行列である。Jameson およびTurkelによって近似されたフラックス・ヤコビアンＪ±を使用することによって、Ｄ対角化形を以下のようにする。
【数２７】
2008165804.tif 000030
ここで、λijは、空間半径
【数２８】
2008165804.tif 000031
である。LU-SGSは、結局２回の掃き出し(sweep)によって解かれる。第１は、低位の（順方向の）掃き出し
【数２９】
2008165804.tif 000032
であり、第２は、上位の（逆方向の）掃き出し
【数３０】
2008165804.tif 000033
である。ここで、
【数３１】
2008165804.tif 000034
である。
【００６３】
境界条件を課すために、インターフェース・セルｉには体積力が導入される。式（２）の元のナビエ・ストークス方程式は以下のように記述される。
【数３２】
2008165804.tif 000035
ＲＨＳｉは、非粘性フラックスおよび粘性フラックスを含み、
【数３３】
2008165804.tif 000036
によって与えられる。体積力Ｓｉは、
【数３４】
2008165804.tif 000037
によって示される。ここで、(・)IBは、ＩＢ法(Immersed Boundary Method)によって課された体積力である。
【００６４】
式（３９）の第２項は、グリッドレス(gridless)法によって評価される。該方法は、クラウドおよびノード位置の形状関数の概念を必要とする。各々のインターフェース・セルのクラウドは、隣接流体セル、インターフェース・セル、境界ポイントおよびサブ・グリッドを含む。インターフェース・セルｉにおける流れ変数は、そのクラウド内の周囲のポイントの流れ変数および形状関数の和によって評価される。K番目の座標方向における任意の特定の関数  の空間導関数は、形状関数の導関数φによって以下のように記述される。
【数３５】
2008165804.tif 000038
ここで、Ｎは、考察対象のクラウド内のノードの数、添え字ｎは、クラウドに属するポイントのインデックスを示す。式（４０）によって、インターフェース・セルとクラウドのメンバー・ポイントの間の全ての中間点において、非粘性フラックスおよび粘性フラックスが計算される。
【００６５】
キューブのサイズ調整方法
図１４は、本実施形態におけるキューブのサイズ調整方法を示す流れ図である。図１４の流れ図は、図１のステップＳ０１０５０およびステップＳ０１０６０におけるキューブのサイズ調整方法の詳細を示すものである。
【００６６】
ステップＳ１４０１０において、全てのキューブの適応パラメータが計算される。適応パラメータは、たとえば、画像処理において輪郭を捉えるために使用されるラプラシアン・フィルタに基づく。流れ解ベクトルの選択されたコンポーネントに対するラプラシアン・フィルタは、以下の式によって与えられる。コンポーネントは、たとえば、マッハ数であってもよい。
【数３６】
2008165804.tif 000039
ここで、Δlは、キューブ・サイズであり、Ｌは、長さスケールである。Ｌの値を増加させると、大きなキューブが細分化される。しかし、増加させすぎると、重要な流れの特性を解くために、十分なキューブのクラスタリングが行われることなく、滑らかな流れ領域において多すぎる数のキューブが生成される。他方、Ｌの値を減少させると、小さなキューブが細分化される。しかし、減少させすぎると、キューブのサイズの分布における空間的な変化が大きくなり、解の精度を劣化させる。本実施形態においては、Ｌの適切な値が、０と１の間で経験的に定められる。
【００６７】
キューブ内の全てのセルのラプラシアンの絶対値の総和を、該キューブの適応パラメータとする。
【００６８】
ステップＳ１４０２０において、全てのキューブの適応パラメータの平均値および標準偏差が計算される。
【００６９】
ステップＳ１４０３０において、全てのキューブの適応パラメータの平均値および標準偏差に基づいて、適応パラメータの細分化のしきい値Ghighおよび粗大化のしきい値Glowが定められる。
【００７０】
ステップＳ１４０４０において、キューブの適応パラメータが細分化のしきい値よりも大きいかどうか判断される。キューブの適応パラメータが細分化のしきい値よりも大きければ、ステップＳ１４０５０に進む。キューブの適応パラメータが細分化のしきい値以下であれば、ステップＳ１４０６０に進む。
【００７１】
ステップＳ１４０５０において、キューブが細分化される。
【００７２】
ステップＳ１４０６０において、キューブの適応パラメータが粗大化のしきい値よりも小さいかどうか判断される。キューブの適応パラメータが粗大化のしきい値よりも小さければ、ステップＳ１４０７０に進む。キューブの適応パラメータが粗大化のしきい値以上であれば、ステップＳ１４０８０に進む。
【００７３】
ステップＳ１４０７０において、キューブが粗大化される。
【００７４】
ステップＳ１４０８０において、隣接するキューブのサイズ調整が行われる。
【００７５】
ステップＳ１４０９０において、キューブの境界情報が構築される。
【００７６】
ステップＳ１４１００において、キューブの計算の負荷を複数のプロセッサに分散させる。
【００７７】
システム構成
図１５は、本発明の一実施形態によるシミュレーション計算システムの構成を示す図である。
【００７８】
シミュレーション計算システムは、計算制御部１０１、格子生成部１０３、入出力部１０５、システム・バス１０７および演算部１０９を含む。演算部１０９は、ｎ台のプロセッサ１０９１乃至１０９ｎを含む。
【００７９】
格子生成部１０３は、キューブの生成やキューブのサイズ調整を行う。演算部１０９は、各キューブの流体計算およびキューブ間情報交換を行う。計算制御部１０１は、複数のプロセッサへのキューブの計算の分散、キューブのサイズ調整の開始時期の判断、計算終了時期の判断などを行う。入出力部１０５は、シミュレーション対象のデータの入力およびシミュレーション結果の出力を行う。
【００８０】
図１６Ａは、本実施形態によるシミュレーション計算システムの性能を示す図である。
【００８１】
図１６Ｂは、図１６Ａに示したシミュレーション計算システムの性能を測定するために、その周りにおいてシミュレーション計算が行われる物体の形状を示す図である。物体はアメッド・ボディ(Ahrmed body)と呼称されるものである。グリッド数は、４百万である。
【００８２】
図１６Ａ（ａ）の横軸および縦軸は、それぞれ、プロセッサ（ＣＰＵ）の台数および該台数のシステムの性能を示す。図１６Ａ（ｂ）の横軸および縦軸は、それぞれ、プロセッサ（ＣＰＵ）の台数および該台数のシステムの性能比率を示す。複数のプロセッサを含むシステムは、プロセッサの台数に比例した性能を備えるのが理想である。ここで、プロセッサの台数に比例した性能に対する実際の性能の比率を性能比率と定義する。理想的には、複数のプロセッサを含むシステムの性能比率は１である。しかし、実際には、プロセッサの台数が増加するにしたがって性能比率は大きく低下する。図１６Ａ（ｂ）に示すように、本実施形態によるシステムにおいては、プロセッサの台数が６４台となっても性能比率は０．８（８０％）以上である。その理由は、本実施形態によるシステムにおいては、計算負荷のほぼ等しいキューブによって、計算負荷をプロセッサ間にほぼ均等に分散できるからである。
【００８３】
図１７乃至２０は、本実施形態のシミュレーション計算方法による計算結果を示す図である。
【００８４】
図１７および図１８は、翼の周りのマッハ数および圧力Ｐを示す図である。流入マッハ数は、０．７３である。レイノルズ数は、６．５６ｘ１０６であり、全流れ領域に対して、乱流を仮定している。迎角は、２．７９度に設定される。衝撃波は、翼の上面において、先端から翼弦の約５５％に形成される。
【００８５】
図１９は、面圧の計算値と実験値との比較を示す図である。図１９の縦軸は圧力係数を示す。図１９の横軸は位置座標を示す。翼の前縁の座標が０、翼の後縁の座標が１である。２つの線のうち上方の線の値が翼の上面の圧力係数値であり、下方の線の値が翼の下面の圧力係数値である。計算値は、全体として実験値によく一致している。
【００８６】
図２０は、計算によって求められた、車両モデル近くの流線の軌跡を示す図である。図２０は、流線における圧力Ｐも示す。
【００８７】
上述のように本発明のシミュレーション計算方法は、複雑な形状への適合性、解適応グリッド制御の適用、高分解能方式の適用および格子生成に優れている。
【００８８】
本発明の実施形態の特徴は以下のとおりである。
【００８９】
本発明の一実施形態において、キューブの分割の回数をツリー構造のレベルと対応させて、キューブの集合を該ツリー構造の最も下位のレベルの要素であるリーフによって表現し、該ツリー構造を使用して、所定のキューブに隣接するキューブを見出して、境界情報を交換するように構成している。
【００９０】
本実施形態によれば、ツリー構造を使用することによって、所定のキューブに隣接するキューブを容易に見出すことができる。
【００９１】
本発明の実施形態において、隣接するキューブは、ツリー構造において同じレベルか隣接レベルであるように構成している。
【００９２】
本実施形態によれば、隣接するキューブは、ツリー構造において同じレベルか隣接レベルであるように構成することによって、隣接するキューブのサイズの比が所定値以内となり、スムージングが行われる。
【００９３】
本発明の実施形態において、キューブが、前記ツリー構造に関し、親および子へのポインタと、系図と、レベルとのデータを有し、該キューブのセルに関し、前記物体の有無のデータを有するように構成している。
【００９４】
本実施形態によれば、キューブの有するデータの量が少なくなる。
【００９５】
本発明の実施形態において、前記物体の有無のデータを、ラン・レングス・データとして記憶するように構成している。
【００９６】
本実施形態によれば、キューブの有するデータの量がさらに少なくなる。
【００９７】
本発明の実施形態において、前記分割が、Ｎを１、２または３として２Ｎ個への均等分割である。
【００９８】
本実施形態によれば、Ｎを変えることにより、キューブの等方分割または非等方分割を行うことができる。
【００９９】
本発明の実施形態において、前記対象領域が３次元の領域であり、キューブが立方体または直方体である。
【０１００】
本発明の実施形態において、所定数の計算ステップを実行した後に、キューブの適応パラメータに基づいて、該キューブのサイズが適切であるかどうか評価し、適切でない場合には該キューブのサイズを変更するように構成している。
【０１０１】
本実施形態によれば、所定数の計算ステップを実行した後にキューブの適応パラメータによって、該キューブのサイズを評価することによってキューブのサイズを調整することができる。その結果、シミュレーション計算の精度が向上する。
【０１０２】
本発明の実施形態において、前記物体と前記流体との境界を含むキューブ内の前記境界の近傍に、前記境界の近傍の分解能を向上させるように、セルとは別に、セルよりも細かな間隔でサブ・グリッドを設けるように構成している。
【０１０３】
本実施形態によれば、セルの数を一定値以内としながら、物体と流体との境界近傍の分解能を向上させることができる。
【０１０４】
本発明の実施形態において、計算の負荷を、単一のキューブの計算を単位として複数のプロセッサに分散させるように構成している。
【０１０５】
複数のキューブは、サイズは異なっても等しい数のセルを備えるので、複数のキューブの計算の負荷はほぼ等しい。したがって、単一のキューブの計算を単位として、所望の大きさの負荷を複数のプロセッサに分散させることができる。
【図面の簡単な説明】
【０１０６】
【図１】本発明の一実施形態によるシミュレーション計算方法を示す流れ図である。
【図２】本実施形態におけるキューブの生成方法を示す流れ図である。
【図３】物体の周りに生成されたキューブを示す図である。
【図４】物体（固体）と流体との境界を含むキューブ及び該キューブ内のセルを示す図である。
【図５】キューブの分割の状態を表現するツリー構造を示す図である。
【図６】非等方に分割されたキューブを示す図である。
【図７】２次元のラン・レングス・データ・ラインを説明するための図である。
【図８】３次元のラン・レングス・データ・ラインを説明するための図である。
【図９】キューブの系図の定義を示すための図である。
【図１０】キューブ位置と系図との関係を示す図である。
【図１１】本実施形態におけるサブ・グリッドの生成方法を示す流れ図である。
【図１２】固体（物体）と流体との境界を含むキューブを示す図である。
【図１３】サブ・グリッドが生成されたインターフェース・セルを示す図である。
【図１４】本実施形態におけるキューブのサイズ調整方法を示す流れ図である。
【図１５】本発明の一実施形態によるシミュレーション計算システムの構成を示す図である。
【図１６Ａ】本実施形態によるシミュレーション計算システムの性能を示す図である。
【図１６Ｂ】図１６Ａに示したシミュレーション計算システムの性能を測定するために、その周りにおいてシミュレーション計算が行われる物体の形状を示す図である。
【図１７】翼の周りのマッハ数を示す図である。
【図１８】翼の周りの圧力を示す図である。
【図１９】面圧の計算値と実験値との比較を示す図である。
【図２０】計算によって求められた、車両近くの流線の軌跡を示す図である。
【符号の説明】
【０１０７】
１０１…計算制御部、１０３…格子生成部、１０５…入出力部、１０９…演算部
【図１】
2008165804.tif 000040
【図２】
2008165804.tif 000041
【図３】
2008165804.tif 000042
【図４】
2008165804.tif 000043
【図５】
2008165804.tif 000044
【図６】
2008165804.tif 000045
【図７】
2008165804.tif 000046
【図８】
2008165804.tif 000047
【図９】
2008165804.tif 000048
【図１０】
2008165804.tif 000049
【図１１】
2008165804.tif 000050
【図１２】
2008165804.tif 000051
【図１３】
2008165804.tif 000052
【図１４】
2008165804.tif 000053
【図１５】
2008165804.tif 000054
【図１６Ａ】
2008165804.tif 000055
【図１６Ｂ】
2008165804.tif 000056
【図１９】
2008165804.tif 000057
【図２０】
2008165804.tif 000058
【図１７】
2008165804.tif 000059
【図１８】
2008165804.tif 000060