(57)【要約】【課題】  計算格子を移動させる際に生じるつぶれや歪みを抑え、高品質な計算格子を形成する計算格子形成装置、計算格子形成方法、及びプログラムを提供する。【解決手段】  格子を移動させる際、移動ベクトル算出部２１０は、移動対象となる格子点が以下の要請を満たすような移動ベクトルΔＰcを決定する。すなわち、格子線がなるべく直交し（直交化）、格子線上に格子点が所定の規則に基づいて滑らかに分布し（平滑化）、初期状態の計算格子の形状を維持し（初期格子の再現性）、格子点間隔が所定の最小間隔以上となり（最小間隔確保）、格子点間隔が最大間隔以下となる（最大間隔確保）ように、移動ベクトルΔＰcを算出し、算出された移動ベクトルΔＰcに従って、格子を移動させる。【選択図】図９
2010026579.tif 000002
【特許請求の範囲】
【請求項１】
計算空間に複数の格子点からなる計算格子を形成し、各格子点を経時的に移動させる計算格子形成装置であって、
移動対象となる格子点の移動ベクトルを算出する移動ベクトル算出手段と、
前記移動ベクトル算出手段によって算出された移動ベクトルに従って前記格子点の移動後の位置を逐次計算する格子移動手段と、を備え、
前記移動ベクトル算出手段は、
前記計算格子の格子線の交差角が垂直または略垂直となるように前記移動ベクトルを定義する直交化要請と、
前記格子点が所定の規則に基づいて前記計算空間に滑らかに分布するように前記移動ベクトルを定義する平滑化要請と、
移動後の計算格子の格子分布が初期状態の該当する計算格子の格子分布に近づくように前記移動ベクトルを定義する初期格子再現要請と、
前記移動対象となる格子点と隣接する格子点との間隔が所定の最小間隔以上となるように前記移動ベクトルを定義する最小間隔確保要請と、
前記移動対象となる格子点と隣接する格子点との間隔が所定の最大間隔以下となるように前記移動ベクトルを定義する最大間隔確保要請と、
に基づき、前記移動ベクトルを算出することを特徴とする計算格子形成装置。
【請求項２】
前記平滑化要請、前記初期格子再現要請、前記最小間隔確保要請、及び前記最大間隔確保要請では、前記初期状態の計算格子の位置情報と前の格子移動ステップにおける各格子の位置情報とに基づき、それぞれ移動ベクトルを定義することを特徴とする請求項１記載の計算格子形成装置。
【請求項３】
前記移動ベクトル算出手段は、前記直交化要請、前記平滑化要請、前記初期格子再現要請、前記最小間隔確保要請、及び前記最大間隔確保要請を満足するため、それぞれ周囲の格子点との相対的な位置関係から所定の関係式を用いて各移動ベクトルを算出し、算出された各移動ベクトルを所定の重みを用いて重ね合わせた線形和を算出し、
前記格子移動手段は、算出された線形和に従って、各格子点の位置を移動することを特徴とする請求項１記載の計算格子形成装置。
【請求項４】
前の時刻における格子点の位置情報に基づいて、前記移動ベクトル算出手段に対して初期値を外挿する初期値外挿手段を更に備えることを特徴とする請求項１記載の計算格子形成装置。
【請求項５】
計算空間に複数の格子点からなる計算格子を形成し、各格子点を経時的に移動させる計算格子形成方法であって、
移動対象となる格子点の移動ベクトルを算出する移動ベクトル算出ステップと、
前記移動ベクトル算出ステップによって算出された移動ベクトルに従って前記格子点の移動後の位置を逐次計算する格子移動ステップと、を備え、
前記移動ベクトル算出ステップは、
前記計算格子の格子線の交差角が垂直または略垂直となるように前記移動ベクトル定義する直交化要請と、
前記格子点が所定の規則に基づいて前記計算空間に滑らかに分布するように前記移動ベクトルを定義する平滑化要請と、
移動後の計算格子の格子分布が初期状態の該当する計算格子の格子分布に近づくように前記移動ベクトルを定義する初期格子再現要請と、
前記移動対象となる格子点と隣接する格子点との間隔が所定の最小間隔以上となるように前記移動ベクトルを定義する最小間隔確保要請と、
前記移動対象となる格子点と隣接する格子点との間隔が所定の最大間隔以下となるように前記移動ベクトルを定義する最大間隔確保要請と、
に基づき、前記移動ベクトルを算出することを特徴とする計算格子形成方法。
【請求項６】
前記平滑化要請、前記初期格子再現要請、前記最小間隔確保要請、及び前記最大間隔確保要請では、前記初期状態の計算格子の位置情報と前の格子移動ステップにおける各格子の位置情報とに基づき、それぞれ移動ベクトルを定義することを特徴とする請求項５記載の計算格子形成方法。
【請求項７】
前記移動ベクトル算出ステップは、前記直交化要請、前記平滑化要請、前記初期格子再現要請、前記最小間隔確保要請、及び前記最大間隔確保要請を満足するため、それぞれ周囲の格子点との相対的な位置関係から所定の関係式を用いて各移動ベクトルを算出し、算出された各移動ベクトルを所定の重みを用いて重ね合わせた線形和を算出し、
前記格子移動ステップは、算出された線形和に従って、各格子点の位置を移動することを特徴とする請求項５記載の計算格子形成方法。
【請求項８】
前の時刻における格子点の位置情報に基づいて、前記移動ベクトル算出ステップに対して初期値を外挿する初期値外挿ステップを更に備えることを特徴とする請求項５記載の計算格子形成方法。
【請求項９】
コンピュータを請求項１記載の計算格子形成装置として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、計算格子形成装置、計算格子形成方法、及びプログラムに関し、詳細には、流体解析、及び構造解析等に用いられる計算格子を形成する計算格子形成装置、計算格子形成方法、及びプログラムに関する。
【背景技術】
【０００２】
近年コンピュータ等を用いて、物理現象を数値解析する様々な手法が研究、開発されている。例えば、移動物体周囲の流体解析手法の一つとして、ＡＬＥ（Arbitrary Lagrangian Eulerian formulations）法が知られている。ＡＬＥ法では、解析対象となる物理空間を多数の格子点で離散化し、物理空間をメッシュ状の計算格子（以下、格子という）からなる計算空間に写像する。そして、この格子を経時的に移動、変形させる等して移動物体周囲の流れを解析する。ここで形成される格子の質は流体解析の精度や安定性に大きく寄与する。
【０００３】
初期状態の格子（以下、初期格子）の形成法としては、（１）Geometrical Method、（２）代数式による格子形成法（Transfinite法）、（３）座標に対する楕円型偏微分方程式による格子形成法等の手法が知られている。また、格子を初期状態から移動させる格子移動法としては、（４）スプリングアナロジーによる格子移動法、（５）弾性力学的アプローチによる格子移動法、（６）速度場の偏微分方程式による平滑化、等の手法が知られている。
【０００４】
また、（５）の弾性力学的な手法を利用して、複雑な対象に対して形状が一部変わったときの格子を作成する際、格子を弾性体と見立て、既に存在する格子の外形を変えて構造解析を行うことで、内部の格子の移動ベクトルを求め、格子の作成を簡略化する手法が開示されている（特許文献１）。
また、構造解析を対象としては、大変形を伴う場合に物体の変形が過度になりすぎないように各座標方向に対して、幾何学的な関係に基づいて変形の程度を補正する手法が開示されている（特許文献２）。
また、直交座標系において、物理量の変化の大きな箇所の格子間隔が密になるよう代数的に格子を移動させる手法が開示されている（特許文献３）。
【０００５】
【特許文献１】特開平９−１２８４３６号公報
【特許文献２】特開平７−３０２３５４号公報
【特許文献３】特開平７−０６３６４３号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、上述の特許文献１の手法は、直前の格子に対するつりあい方程式をとくため、初期の格子情報が反映されず、空間解像度の制御が難しい。すなわち格子が一度つぶれてしまうと境界の位置が元の位置に戻ってもその格子は元の大きさに戻らないという不具合が起きる。また一つの剛性パラメータで格子の直交化と集中化を調節しているため、複雑な形状では対象や運動によってこのパラメータを試行錯誤的に調節しなければならないという問題がある。
【０００７】
また、特許文献２の手法は、各格子線方向に対して過度に格子が変形しないよう補正を独立に求めるため、２次元もしくは３次元的な歪みを抑えることができないという問題がある。また、特許文献３の手法は、各格子線方向に対して補正を独立に求めるため、２次元もしくは３次元的な歪みを抑えることができない。さらに直交座標系を前提としているため、複雑な形状には対応できないという問題がある。
【０００８】
図１３は、従来手法を用いて角柱状物体６０の回転させた場合における流れ解析の結果出力例であり、（ａ）は初期状態の全体図、（ｂ）は（ａ）の物体６０の左上角部分を含む計算空間を拡大した図、（ｃ）は最大変位時における全体図、（ｄ）は（ｃ）の角柱状物体６０の上角部分の計算空間を拡大した図である。
なお、ここでは、上述のスプリングアナロジーを構造格子に拡張する手法を用いて格子を移動している。角柱の回転角度は最大４５度とする。
図１３に示す従来手法では、特に壁近傍の角部で格子がつぶれてしまっている。
【課題を解決するための手段】
【０００９】
本発明は、前述した問題点に鑑みてなされたもので、その目的は、計算格子を移動させる際に生じるつぶれや歪みを抑え、高品質な計算格子を形成する計算格子形成装置、計算格子形成方法、及びプログラムを提供することである。
【００１０】
前述した目的を達成するために第１の発明は、計算空間に複数の格子点からなる計算格子を形成し、各格子点を経時的に移動させる計算格子形成装置であって、移動対象となる格子点の移動ベクトルを算出する移動ベクトル算出手段と、前記移動ベクトル算出手段によって算出された移動ベクトルに従って前記格子点の移動後の位置を逐次計算する格子移動手段と、を備え、前記移動ベクトル算出手段は、前記計算格子の格子線の交差角が垂直または略垂直となるように前記移動ベクトルを定義する直交化要請と、前記格子点が所定の規則に基づいて前記計算空間に滑らかに分布するように前記移動ベクトルを定義する平滑化要請と、移動後の計算格子の格子分布が初期状態の該当する計算格子の格子分布に近づくように前記移動ベクトルを定義する初期格子再現要請と、前記移動対象となる格子点と隣接する格子点との間隔が所定の最小間隔以上となるように前記移動ベクトルを定義する最小間隔確保要請と、前記移動対象となる格子点と隣接する格子点との間隔が所定の最大間隔以下となるように前記移動ベクトルを定義する最大間隔確保要請と、に基づき、前記移動ベクトルを算出することを特徴とする計算格子形成装置である。
【００１１】
また、前記平滑化要請、前記初期格子再現要請、前記最小間隔確保要請、及び前記最大間隔確保要請では、前記初期状態の計算格子の位置情報と前の格子移動ステップにおける各格子の位置情報とに基づき、それぞれ移動ベクトルを定義する。
【００１２】
前記移動ベクトル算出手段は、前記直交化要請、前記平滑化要請、前記初期格子再現要請、前記最小間隔確保要請、及び前記最大間隔確保要請を満足するため、それぞれ周囲の格子点との相対的な位置関係から所定の関係式を用いて各移動ベクトルを算出し、算出された各移動ベクトルを所定の重みを用いて重ね合わせた線形和を算出し、前記格子移動手段は、算出された線形和に従って、各格子点の位置を移動する。
【００１３】
また、前の時刻における格子点の位置情報に基づいて、前記移動ベクトル算出手段に対して初期値を外挿する初期値外挿手段を更に備えるようにしてもよい。
【００１４】
第２の発明は、計算空間に複数の格子点からなる計算格子を形成し、各格子点を経時的に移動させる計算格子形成方法であって、移動対象となる格子点の移動ベクトルを算出する移動ベクトル算出ステップと、前記移動ベクトル算出ステップによって算出された移動ベクトルに従って前記格子点の移動後の位置を逐次計算する格子移動ステップと、を備え、前記移動ベクトル算出ステップは、前記計算格子の格子線の交差角が垂直または略垂直となるように前記移動ベクトル定義する直交化要請と、前記格子点が所定の規則に基づいて前記計算空間に滑らかに分布するように前記移動ベクトルを定義する平滑化要請と、移動後の計算格子の格子分布が初期状態の該当する計算格子の格子分布に近づくように前記移動ベクトルを定義する初期格子再現要請と、前記移動対象となる格子点と隣接する格子点との間隔が所定の最小間隔以上となるように前記移動ベクトルを定義する最小間隔確保要請と、前記移動対象となる格子点と隣接する格子点との間隔が所定の最大間隔以下となるように前記移動ベクトルを定義する最大間隔確保要請と、に基づき、前記移動ベクトルを算出することを特徴とする計算格子形成方法である。
【００１５】
また、第３の発明は、コンピュータを第１の発明の計算格子形成装置として機能させるためのプログラムである。
【発明の効果】
【００１６】
本発明により、計算格子を移動させる際に生じるつぶれや歪みを抑え、高品質な計算格子を形成する計算格子形成装置、計算格子形成方法、及びプログラムを提供することができる。
【発明を実施するための最良の形態】
【００１７】
以下図面に基づいて、本発明の実施形態を詳細に説明する。
【００１８】
図１は、本実施の形態の計算格子形成装置１により形成される計算格子について説明する図である。
図２は、本実施の形態に係る計算格子形成装置１を実現するコンピュータのハードウェア構成図である。
【００１９】
計算格子形成装置１は、解析対象となる物理空間を計算空間５０に写像するために、計算空間５０を複数の格子点で離散化し、計算格子（以下、格子という）を形成する。図１に示す計算空間５０は縦方向及び横方向に交差する複数の格子線によって表現される。格子線の交差する点が格子点である。
【００２０】
図１は、計算空間５０の中央に置かれた角柱状の物体６０を上面から見た様子を表す。図１の計算格子は、物体６０の各辺の延長線上で、格子線間隔が狭くなるように形成されている。
【００２１】
なお、図１では以下の説明の簡略化のために２次元の計算空間５０を例示したが、３次元としてもよい。格子の形状は四角形に限定さるものではなく六面体としてもよい。また、計算空間５０ではその座標を、２次元の場合は（ξ、η）で表し、３次元の場合は（ξ、η、ζ）で表すものとする。
また、格子ベクトルとは各格子点の位置から決定されるベクトルを表し、移動ベクトルとは、移動対象となる格子点の移動すべき向きと長さを定義するベクトルを表すものとする。
【００２２】
図２に示すように、計算格子形成装置１は、制御部３、記憶部５、メディア入出力部７、通信制御部９、入力部１１、表示部１３、周辺機器Ｉ／Ｆ部１５等が、バス１７と介して接続される。
【００２３】
制御部３は、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等で構成される。
【００２４】
ＣＰＵは、記憶部５、ＲＯＭ、記録媒体等に格納されるプログラムをＲＡＭ上のワークメモリ領域に呼び出して実行し、バス１７を介して接続された各装置を駆動制御し、計算格子形成装置１が行う後述する処理（図９、図１１、図１２参照）を実現する。
ＲＯＭは、不揮発性メモリであり、コンピュータのブートプログラムやＢＩＯＳ等のプログラム、データ等を恒久的に保持している。
ＲＡＭは、揮発性メモリであり、記憶部５、ＲＯＭ、記録媒体等からロードしたプログラム、データ等を一時的に保持するとともに、制御部３が各種処理を行う為に使用するワークエリアを備える。
【００２５】
記憶部５は、ＨＤＤ（ハードディスクドライブ）であり、制御部３が実行するプログラム、プログラム実行に必要なデータ、ＯＳ（オペレーティングシステム）等が格納される。プログラムに関しては、ＯＳ（オペレーティングシステム）に相当する制御プログラムや、アプリケーションプログラムが格納されている。
これらの各プログラムコードは、制御部３により必要に応じて読み出されてＲＡＭに移され、ＣＰＵに読み出されて各種の手段として実行される。
【００２６】
メディア入出力部７（ドライブ装置）は、データの入出力を行い、例えば、フロッピー（登録商標）ディスクドライブ、ＰＤドライブ、ＣＤドライブ（−ＲＯＭ、−Ｒ、ＲＷ等）、ＤＶＤドライブ（−ＲＯＭ、−Ｒ、−ＲＷ等）、ＭＯドライブ等のメディア入出力装置を有する。
【００２７】
通信制御部９は、通信制御装置、通信ポート等を有し、コンピュータとネットワーク１９間の通信を媒介する通信インタフェースであり、ネットワーク１９を介して、他のコンピュータ間との通信制御を行う。
【００２８】
入力部１１は、例えば、キーボード、マウス等のポインティング・デバイス、テンキー等の入力装置であり、データの入力を行う。
【００２９】
表示部１３は、液晶パネル、ＣＲＴモニタ等のディスプレイ装置、及びディスプレイ装置と連携してコンピュータのビデオ機能を実現するための論理回路等（ビデオアダプタ等）を有する。
【００３０】
周辺機器Ｉ／Ｆ（インタフェース）部１５は、コンピュータに周辺機器を接続させるためのポートであり、周辺機器Ｉ／Ｆ部１５を介してコンピュータは周辺機器とのデータの送受信を行う。周辺機器Ｉ／Ｆ部１５は、ＵＳＢやＩＥＥＥ１３９４やＲＳ−２３２Ｃ等で構成されており、通常複数の周辺機器Ｉ／Ｆを有する。周辺機器との接続形態は有線、無線を問わない。
【００３１】
バス１７は、各装置間の制御信号、データ信号等の授受を媒介する経路である。
【００３２】
次に図３を参照しながら、計算格子形成装置１の機能構成について説明する。
図３（ａ）は、計算格子形成装置１全体の機能概略構成図であり、図３（ｂ）は、格子移動演算部２００についてのブロック図である。
【００３３】
計算格子形成装置１は、図３に示すように、初期格子形成部１００、格子移動演算部２００、流体解析部３００、及び結果出力部４００を備えて構成される。
初期格子形成部１００は、初期状態の格子（以下、初期格子という）を形成する。初期格子形成にはいかなる手法を用いてもよい。例えばGeometrical Method、代数式による格子形成法（Transfinite法）、座標に対する楕円型偏微分方程式による格子形成法等といった公知の手法を用いてもよい。また、予め初期格子に関するデータを記憶部５に記憶しておき、解析実行時に初期格子形成部１００により記憶部５から読み出すようにしてもよい。初期格子形成部１００は形成した初期格子を格子移動演算部２００に出力する。
【００３４】
格子移動演算部２００は、初期格子形成部１００から入力される初期格子と、入力部１１から入力されるか、或いは記憶部５から読み込んだ解析対象となる流れ場の情報及び格子情報ｄ１０等と、に基づいて格子を移動させるための演算を実行する。格子の移動とは、時刻ステップ毎に計算空間５０内の各格子点の位置を変化させることであり、格子移動演算部２００により逐次計算される。格子移動演算部２００の行う処理については後述する。
【００３５】
流体解析部３００は、格子移動演算部２００から入力される移動後の格子について、ヤコビ行列式、メトリック関係式等を算出し、流体計算を行う。例えば、流れ場や、各種統計量の算出、格子の経時変化等の解析を行う。そして、解析結果を結果出力部４００に出力する。流体解析部３００の計算後、再び格子移動演算部２００に戻り、次の時刻での格子位置を求め、流体計算を行う。結果出力部４００は流体解析部３００から入力される解析結果を、数値やグラフ或いはグラフィック表示にて表示部１３に表示する。
【００３６】
次に、格子移動演算部２００について、図３（ｂ）、図４、図５、図６、図７、図８を参照して説明する。
格子移動演算部２００は、図３（ｂ）に示すように、直交化部２１１、平滑化部２１２、初期格子再現部２１３、最小間隔確保部２１４、最大間隔確保部２１５、及び線形和算出部２１６を有する移動ベクトル算出部２１０と、格子移動部２２０とを備えて構成される。
【００３７】
移動ベクトル算出部２１０は、各ステップ毎の各格子点の移動ベクトルを算出する。移動ベクトルはその方向及び長さにより定義される。
更に本発明の計算格子生成装置１では、移動ベクトルは、直交化部２１１、平滑化部２１２、初期格子再現部２１３、最小間隔確保部２１４、及び最大間隔確保部２１５の各部によるそれぞれの要請に従って決定される。
すなわち、直交化部２１１により各格子線がなるべく直交（直交または略直交）に交わり、格子の歪みや捩れを抑える。また、平滑化部２１２により格子線上に格子点が滑らかに分布し、格子間隔や位置の急激な変化を少なくする。また、初期格子再現部２１３により、初期の格子分布をなるべく保つ。また、最小間隔確保部２１４により格子間隔が狭くなりすぎないようにする。また、最大間隔確保部２１５により格子間隔が広くなりすぎないようにする。
【００３８】
このように、移動ベクトルの算出の際に直交化等の要請を直接入れることにより、格子がつぶれにくいという効果を得る。また、移動ベクトルの算出の際に格子点の位置に関する要請を直接与えているため、各要請の重みをユーザが直感的にコントロールしやすくなる。
【００３９】
また、移動ベクトル算出部２１０は、直交化以外の各要請（平滑化部２１２、初期格子再現化部２１３、最小間隔確保部２１４、最大間隔確保部２１５）に関しては、初期格子の位置情報と、直前の格子移動ステップによって再構成された格子位置情報とを反映させて、移動ベクトルを算出するものとする。
また、各要請に対応する移動ベクトルの算出には、次に説明するように、周囲の格子点との相対的な位置関係から所定の幾何学的な関係式（後述の式（１）〜式（９））を用いて算出する。
【００４０】
図４は、直交化部２１１により決定される移動ベクトルΔＰcｏｒｔｈについて説明する図であり、図５は、平滑化部２１２により決定される移動ベクトルΔＰcｓｍｏについて説明する図であり、図６は、初期格子再現部２１３により決定される移動ベクトルΔＰcｉｎｉについて説明する図であり、図７は、最小間隔確保部２１４により決定される移動ベクトルΔＰcｍｉｎについて説明する図であり、図８は、最大間隔確保部２１５により決定される移動ベクトルΔＰcｍａｘについて説明する図である。
図４〜図８において、Ｐnw、Ｐn、Ｐne、Ｐe、Ｐse、Ｐs、Ｐsw、Ｐw、Ｐcは、それぞれ格子点を表す。これらの各格子点のうち隣接する格子点を結ぶ線が格子線である。
なお、図４〜図８ではベクトルを表す符合として矢印（→）を付しているが、明細書内には、矢印を付さない。代わりに、例えば点ＰとベクトルＰとを区別するため、ベクトルを意味する場合には「ベクトルＰ」と記載することとする。
【００４１】
図４に示すように、直交化部２１１は、各格子線の交差角が垂直または略垂直に交差するように、移動対象となる格子点（中央の格子点Ｐc）の移動ベクトルを算出する。具体的には、格子点Ｐcは、まず、点Ｐs側からの直交性を満足するために、点Ｐsから立てた垂線上（例えば点Ｐsを通り、点Ｐw、Ｐe、Ｐsw、Ｐseから算出されるベクトルｂに垂直な線上）に移動することが望ましい。この移動ベクトルをΔＰc（ｓ）とおく。
最終的に求める直交性を向上するための移動ベクトルΔＰcｏｒｔｈは、２次元の場合では周囲４方向、３次元の場合では次の式（１）のように周囲６方向の移動ベクトルを重ね合わせて得られる。
【００４２】
【数１】
2010026579.tif 000003
【００４３】
ここで、ｗｗ、ｗｅ、ｗｓ、ｗｎ、ｗｈ、ｗｌは、それぞれｗ方向、ｅ方向、ｓ方向、ｎ方向、ｈ方向、ｌ方向の重みであり、すべて足し合わせると１となるように算出する。
また、ｗｗ、ｗｅ、ｗｓ、ｗｎ、ｗｈ、ｗｌは、各方向の格子の長さ（格子点間隔）や直角からのずれ角を加味してその重みのバランスを変化させるものとする。
【００４４】
平滑化部２１２は、格子点Ｐcが格子線上に所定の規則に従って滑らかに分布するように移動ベクトルΔＰcを定義する。ここでいう「滑らかさ」及び「平滑性」とは、格子点が一部に集中したり、格子間隔比がばらついたりせずに、所定の規則に従って分布することをいう。例えば、間隔や間隔比が等しい場合、或いは等比級数に従って分布する場合等である。
具体的には、図５に示すような格子配置を考えた場合、η方向の平滑性を高める移動ベクトルΔＰc（η）は以下のように与えられる。
【００４５】
【数２】
2010026579.tif 000004
【００４６】
式（２）において、上付き添え字「０」は初期の格子位置を意味する。また、式（２）は、対象とする格子点Ｐcが、両隣の格子Ｐs、Pnを結ぶ直線上に移動するとともに初期の格子間隔比を維持する位置に移動することを示している。
３次元での平滑化の格子移動ベクトルΔＰcｓｍｏは、ξ、η、ζの各方向の移動ベクトルを重ね合わせて次の式（３）のように算出される。
【００４７】
【数３】
2010026579.tif 000005
【００４８】
ここで、ｗξ、ｗη、ｗζは、それぞれξ、η、ζの各方向の重みであり、すべて足し合わせると１となるように算出する。
また、この重みｗξ、ｗη、ｗζは各方向の格子の初期長さに応じて変化するように定義する。
【００４９】
初期格子再現部２１３は、移動後の格子分布を初期状態の格子分布に近づけるように移動ベクトルを定義する。ここで、格子分布とは、各格子線上に存在する格子点の分布を意味するものである。「格子分布を初期状態の格子分布に近づける」とは、格子線全体が伸縮した場合でもそれに伴い格子点も移動し、格子点の相対的な位置関係（分布）は変化しないという意味である。
具体的には、η方向に対して図６のような格子配置を考えた場合、格子点Ｐcに対する点Ｐs側からの初期格子の再現性に関する移動ベクトルΔＰc（ｓ）は以下の式で与えられる。
【００５０】
【数４】
2010026579.tif 000006
【００５１】
式（４）は、点Ｐcと点Ｐsとの格子間隔が、初期の同点間の格子間隔と等しくなるように、点Ｐcを、両隣の格子Ｐs及びＰnを結ぶベクトルＰsＰnの方向に移動することを表している。最終的な初期格子の再現性のための移動ベクトルΔＰcｉｎｉは、３次元では次の式（５）のように周囲６方向の移動ベクトルを重ね合わせて得られる。
【００５２】
【数５】
2010026579.tif 000007
【００５３】
ここで、ｗｗ、ｗｅ、ｗｓ、ｗｎ、ｗｈ、ｗｌは、それぞれｗ方向、ｅ方向、ｓ方向、ｎ方向、ｈ方向、ｌ方向の重みであり、すべて足し合わせると１となるように算出する。
また、この重みｗｗ、ｗｅ、ｗｓ、ｗｎ、ｗｈ、ｗｌは、格子の初期長さに応じて変化するように定義する。
【００５４】
最小間隔確保部２１４は、各格子点の最小間隔を確保するように移動ベクトルを定義する。具体的には、η方向に対して、図７のような格子配置を考えると、Ｐs側の最小格子間隔を確保するための移動ベクトルΔＰc（ｓ）は次の式（６）で与えられる。
【００５５】
【数６】
2010026579.tif 000008
【００５６】
式（６）は、格子点Ｐcと格子点Ｐsとの距離が許容最小格子間隔Δｍｉｎ以下になった場合に、点Ｐcが点Ｐs側の許容最小格子間隔と等しくなるように、両隣の格子点ＰsとＰnとを結ぶベクトルＰsＰnの方向に移動することを表している。
【００５７】
ここで許容最小格子間隔には一定値を用いてもよいし、格子毎に違う値を用いてもよい。
反対側（Ｐn側）の格子間隔も許容最小格子間隔以下となる場合には、両者の移動ベクトルを格子長さに関する重み付けのもとで重ね合わせることでη方向の移動ベクトルを求める。
最終的に最小格子間隔の確保に関する移動ベクトルΔＰｃｍｉｎは、ξ、η、ζ方向の各移動ベクトルを重ね合わせることにより、以下の式（７）で求められる。
【００５８】
【数７】
2010026579.tif 000009
【００５９】
最大間隔確保部２１５は、格子点の最大間隔を確保するように移動ベクトルを定義する。具体的には、η方向に対して、図８のような格子配置を考えると、Ｐs側の最大格子間隔の確保のための移動ベクトルΔＰc（ｓ）は以下の式（８）で与えられる。
【００６０】
【数８】
2010026579.tif 000010
【００６１】
式（８）は点Ｐcと点Ｐsとの距離が許容最大格子間隔Δｍａｘ以上になった場合に、点Ｐcが点Ｐs側の許容最大格子間隔と等しくなるように、両隣の格子ＰsとＰnとを結ぶベクトルＰsＰnの方向に移動することを表している。
ここで、許容最大格子間隔は一定値を用いてもよいし、格子毎に違う値を用いてもよい。反対側（Ｐn側）の格子間隔も許容最大格子間隔以上となる場合には、両者の移動ベクトルを格子長さに関する重み付けのもとで重ね合わせることにより、η方向の移動ベクトルを求める。
最終的に、最大格子間隔の確保に関する移動ベクトルΔＰcｍａｘはξ、η、ζ方向の寄与を重ね合わせることにより、以下の式（９）で求められる。
【００６２】
【数９】
2010026579.tif 000011
【００６３】
線形和算出部２１６は、移動ベクトル算出部２１０の直交化部２１１、平滑化部２１２、初期格子再現部２１３、最小間隔確保部２１４、最大間隔確保部２１５によって算出された各格子点の移動ベクトルの線形和を算出して、次の時刻の格子移動ベクトルΔＰcとする。線形和は次式（１０）で表される。
【００６４】
【数１０】
2010026579.tif 000012
【００６５】
ここで、Ｃｏｒｔｈ、Ｃｓｍｏ、Ｃｉｎｉ、Ｃｍｉｎ、Ｃｍａｘは、それぞれ直交性、平滑性、初期格子再現性、最小格子間隔の確保、最大格子間隔の確保に関する各重みである。
また、ｄｔｓｕｂは、後述する逐次計算に用いるサブステップ（図９参照）の擬似時間刻みである。
【００６６】
格子移動部２２０は、線形和算出部２１６により算出された移動ベクトルの線形和ΔＰcに従って格子点の移動後の位置を逐次計算する。
【００６７】
次に、図９及び図１０を参照して、計算格子形成装置１における動作について説明する。
【００６８】
図９は、計算格子形成装置１の処理手順を示すフローチャートである。
図１０は、本実施の形態の計算格子形成装置１によって得られる結果例を示す図である。
【００６９】
計算格子形成装置１は、図９のフローチャートの手順に従って、格子を形成し、移動させ、流体解析計算を実行する。
図９において、左側が通常の流体計算の部分、右側が本実施の形態の計算格子形成装置１の実行する格子移動処理に関する部分である。格子移動処理は流体計算の時間進行の一部に組み込まれる。
【００７０】
まず、初期格子形成部１００は、初期格子に関する情報を記憶部５から読み込み（ステップＳ１０１）、格子移動演算部２００へ渡す。格子移動演算部２００は、格子を移動する上で必要となる初期格子情報（格子ベクトル、格子間隔）を算出し（ステップＳ１０２）、ＲＡＭに保存する。次に、初期格子形成部１００は、流れ場及び格子に関する情報を有するファイルを記憶部５から読み込み（ステップＳ１０３）、格子移動演算部２００へ渡す。
【００７１】
格子移動演算部２００は、読み込まれた初期格子、流れ場及び格子に関する情報に基づいて、格子を移動する処理を行う。以降の処理では、各格子の移動ベクトルを、流体計算の擬似時間刻み（以下、時刻ステップという）毎に算出する。その際、各時刻ステップ内にサブステップを設け、移動ベクトルを陽解法的に逐次計算する。
【００７２】
まず、ステップＳ１０４において、物体６０との境界を次の時刻における位置に移動させる。ここで、サブステップの初期値を設定する。初期値には、直前の時刻の格子情報を用いる。或いは、初期値には直前の時刻以前の時刻の格子情報から外挿して与えるようにしてもよい。例えば、直前の２時刻の格子情報を初期値として与えるようにする。（ステップＳ１０５）。ただし、第１番目の時刻ステップについてはステップＳ１０５をスキップし、ステップＳ１０６へ移行する。初期値を外挿して与えた場合には、解の収束が高速化し、壁近傍の格子品質が向上することが期待できる。
【００７３】
ステップＳ１０６において、移動ベクトル算出部２１０は、現在の格子点の位置情報から格子ベクトルと格子間隔を算出する。そしてステップＳ１０７において、移動ベクトル算出部２１０は、各格子の移動ベクトルを算出する。
【００７４】
移動ベクトルは、上述したようにまず直交化部２１１、平滑化部２１２、初期格子再現部２１３、最小間隔確保部２１４、及び最大間隔確保部２１５の各部によるそれぞれの要請が満たされるように各移動ベクトルΔＰcｏｒｔｈ、ΔＰcｓｍｏ、ΔＰcｉｎｉ、ΔＰcｍｉｎ、ΔＰcｍａｘが算出される。そして、算出された各部の要請による各移動ベクトルΔＰcｏｒｔｈ、ΔＰcｓｍｏ、ΔＰcｉｎｉ、ΔＰcｍｉｎ、ΔＰcｍａｘの線形和ΔＰcが線形和算出部２２０によって算出される。すなわち、直交化、平滑化、初期格子の再現性、最小格子間隔の確保、及び最大格子間隔の確保のすべての要請を満たす移動ベクトルΔＰcが決定される。格子移動部２２０は決定された移動ベクトルΔＰcに従って、格子の位置を移動させる（ステップＳ１０８）。
【００７５】
移動ベクトル算出部２１０はステップＳ１０６〜ステップＳ１０８の処理をＫ回繰り返し、上限回数（Ｋmax）に達すると（ステップＳ１０９；Ｙｅｓ）、その時の移動ベクトルによって次のステップの格子位置が最終的に決定する。
その後、現ステップ及び前ステップでの各格子位置から格子移動速度を算出して（ステップＳ１１０）、流体計算部３００へ移行する。
【００７６】
流体解析部３００では、移動後の格子位置情報を用いて、公知の流体計算処理と同様にヤコビアン（ヤコビ行列式）、メトリック関係式を算出し（ステップＳ１１１）、当該時刻ステップでの流体計算を実行する（ステップＳ１１２）。流体解析部３００は、ステップS１０４〜ステップＳ１１２の処理をＮ回の時刻ステップ分だけ繰り返し、上限回数（Ｎmax）に達すると（ステップＳ１１３；Ｙｅｓ）、最終的な流れ場の状態や各種統計量の時間履歴、そのときの格子位置等を含む解析結果を結果出力部４００に渡す。結果出力部４００は、流体解析結果を数値やグラフ化し、或いはグラフィック表示により出力する（ステップＳ１１４）。
【００７７】
図１０は、本実施の形態の計算格子形成装置１によって得られる、角柱状物体６０の回転問題の結果出力例であり、（ａ）は初期状態の全体図、（ｂ）は（ａ）の物体６０の左上角部分を含む計算空間を拡大した図、（ｃ）は最大変位時の全体図、（ｄ）は（ｃ）の角柱状物体６０の上角部分を含む計算空間を拡大した図である。
角柱状物体６０の運動は以下の式（１１）により与えられるものとする。
【００７８】
【数１１】
2010026579.tif 000013
【００７９】
式（１１）において、θは回転角度を表す。
また、最大回転角θmaxは４５度とする。
【００８０】
図１０（ｄ）に示すように、本実施の形態の計算格子形成装置１によって形成された格子は最大変位時においてもつぶれや歪みが少ない。特に壁（物体）近傍において、その結果が顕著に表れている。
【００８１】
以上説明したように、本実施の形態の計算格子生成装置１は、格子移動演算部２００によって格子を移動させる際に、移動ベクトル算出部２１０によって、移動対象となる格子点の移動ベクトルを算出する。そのとき算出される移動ベクトルは、直交性、平滑性、初期格子再現性、最小格子間隔の確保、最大格子間隔の確保の各要請を満足するように定義される。その結果、つぶれや歪みが少なく、高品質な格子を形成することが可能となり、流体解析の精度を向上させることができる。
また、初期格子と直前の格子位置情報とを参照しながら移動ベクトルを算出するため、初期の格子の格子分布（空間解像度）を反映しつつ、かつ時間の進行に伴った歪みの少ない格子を構成できる。
【００８２】
なお、上述の図９のフローチャートでは、境界の移動が既知の場合の処理動作を示しているが、例えば、構造解析や機構解析といった他の解析手法と連成させてもよい。その場合は、図９のフローチャートに代えて、図１１または図１２のフローチャートに基づく処理を行う。また、図３（ａ）の機能ブロック図において、「流体解析部３００」に「構造解析または機構解析との連成解析部」が設けられる。
【００８３】
図１１及び図１２は、本実施の形態の移動ベクトル算出手法を用いて格子を形成した際に、構造解析または機構解析と連成させる場合の処理手順を示し、図１１は弱連成の場合、図１２は強連成の場合を示すフローチャートである。
【００８４】
図１１の弱連成の例において、ステップＳ２０１〜ステップＳ２１２は、図９のステップＳ１０１〜ステップＳ１１２と同様の手法で格子を形成し、移動させ、流体計算を実行する。そして、ステップＳ２１３として、流体計算の後に構造計算または機構計算の処理を追加する。その後、ステップＳ２０４〜ステップＳ２１３までの処理を上限回数（Ｎmax）に達するまで繰り返し（ステップＳ２１４）、その後結果を出力して（ステップＳ２１５）、処理を終了する。
【００８５】
また、図１２の強連成の例においても、ステップＳ３０１〜ステップＳ３１２は、図９のステップＳ１０１〜ステップＳ１１２と同様の手法で格子を形成し、移動させ、流体計算を実行する。そして、ステップＳ３１３として、流体計算の後に構造計算または機構計算の処理を追加する（ステップＳ３１３）。強連成の場合は、流体計算と構造計算（または機構計算）の間で、力や変位の値をやり取りし、両者の釣り合いをとる。そのため、時刻ステップＮ内にサブステップＭを設けている。サブステップを上限回数（Ｍmax）に達するまで繰り返しつつ(ステップＳ３１４)、非定常計算における実時間進行に相当する時間ステップを上限回数（Ｎmax）に達するまで繰り返す（ステップＳ３１５）。
その後、結果を出力して（ステップＳ３１６）、処理を終了する。
【００８６】
なお、本実施の形態の計算格子の移動方法は、上述の流体解析（図９）または流体解析と構造解析との連成（図１１、図１２）に適用するだけでなく、構造解析に適用することも可能である。
構造解析に適用する場合には、図９のフローチャートにおいて、ステップＳ１１２における「流体計算」処理を、「構造計算」処理に置き換える。
また、図３（ａ）の機能ブロック図については、流体解析部３００は構造解析部に置き換えられる。
【００８７】
なお、上述の図９、図１１、図１２に示す処理を行うプログラムは、ＣＤ−ＲＯＭ等の記録媒体に保持させて流通させても良いし、通信回線を介して送受信することもできる。
【００８８】
以上、添付図面を参照しながら、本発明に係る計算格子形成装置等の好適な実施形態について説明したが、本発明はかかる例に限定されない。当業者であれば、本願で開示した技術的思想の範疇内において、各種の変更例又は修正例に想到し得ることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
【図面の簡単な説明】
【００８９】
【図１】計算格子形成装置１により形成される計算格子について説明する図
【図２】計算格子形成装置１を実現するコンピュータのハードウェア構成図
【図３】（ａ）計算格子形成装置１全体の機能概略構成図、（ｂ）格子移動部２００についてのブロック図
【図４】直交化部２１１により決定される移動ベクトルについて説明する図
【図５】平滑化部２１２により決定される移動ベクトルについて説明する図
【図６】初期格子再現部２１３により決定される移動ベクトルについて説明する図
【図７】最小間隔確保部２１４により決定される移動ベクトルについて説明する図
【図８】最大間隔確保部２１５により決定される移動ベクトルについて説明する図
【図９】計算格子形成装置１の処理手順を示すフローチャート
【図１０】計算格子形成装置１によって得られる結果出力例を示す図
【図１１】構造解析または機構解析と弱連成させる場合の処理手順を示すフローチャート
【図１２】構造解析または機構解析と強連成させる場合の処理手順を示すフローチャート
【図１３】従来手法を用いた角柱状物体６０の回転問題の結果出力例
【符号の説明】
【００９０】
１………計算格子形成装置
５０………計算格子
６０………物体（角柱状物体）
１００………初期格子形成部
２００………格子移動演算部
２１０………移動ベクトル算出部
２１１………直交化部
２１２………平滑化部
２１３………初期格子再現部
２１４………最小間隔確保部
２１５………最大間隔確保部
２１６………線形和算出部
２２０………格子移動部
３００………流体解析部
４００………結果出力部
【図１】
2010026579.tif 000014
【図２】
2010026579.tif 000015
【図３】
2010026579.tif 000016
【図４】
2010026579.tif 000017
【図５】
2010026579.tif 000018
【図６】
2010026579.tif 000019
【図７】
2010026579.tif 000020
【図８】
2010026579.tif 000021
【図９】
2010026579.tif 000022
【図１０】
2010026579.tif 000023
【図１１】
2010026579.tif 000024
【図１２】
2010026579.tif 000025
【図１３】
2010026579.tif 000026