(57)【特許請求の範囲】
【請求項１】
演算手段の演算処理能力に見合う格子数を入力手段によって入力することで、前記演算手段によって、計算空間を有限個数の格子点で分割し、航空機や翼に働く揚力、横力、抵抗およびモーメントの流体運動を支配する偏微分方程式を離散化した差分方程式の格子収束解を算出し、その算出結果を出力手段によって出力する格子収束解の算出システムにおいて、
前記入力手段、演算手段および出力手段を含むコンピュータを、
単一の格子を用いた高次精度の計算手法による計算値と、前記単一の格子を用いた低次精度の計算手法による計算値とを求める手段、および
高次精度の計算値と低次精度の計算値と格子収束解との比が、いずれの格子数においても略一定として、前記単一の格子における高次精度の計算値と低次精度の計算値との関係から格子収束解を算出して表示する手段として機能させるための格子収束解の算出システム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、格子収束解の算出システムに関し、たとえば、回転翼航空機のロータ・ブレードなどの翼断面を決定するためなどに好適に用いることができる格子収束解の算出システムに関する。
【０００２】
【従来の技術】
近年、コンピュータを利用して流体力学の支配方程式を数値的に解き、様々な流体力学の問題に対する答えを出す数値流体力学（Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｆｌｕｉｄ  Ｄｙｎａｍｉｃｓ，ＣＦＤ）が飛躍的に進歩している。たとえば、回転翼航空機のロータ・ブレードなどの翼断面を決定するためには、風洞試験や実機試験を行う前に、前記ＣＦＤを用いて、所望の翼断面をある程度解析的に推定している。
【０００３】
図８は、ＣＦＤを用いて解析した、格子数と、翼の中央付近の代表的２次元断面における抗力係数との関係を示す図表である。ラインＬ１に示すように、抗力係数ＣＤの値は、格子数の増加にともなって予測される格子収束解（実線Ｌ１で図示）に漸近する。すなわち、格子が極限的に細かくなったときの解が格子収束解であり、計算格子を細かくすればするほど算出する収束解の誤差が減少する。しかし、現実的には、コンピュータの演算処理能力の限界があり、必要十分な数の計算格子を設定することができず、満足する計算精度を得られないことが多い。
【０００４】
図９は、翼形流れ計算用の粗い格子をディスプレイ１に表示した場合の説明図であり、図１０は、翼形流れ計算用の細かい格子をディスプレイ１に表示した場合の説明図である。図１１は、従来の格子収束解および誤差の計算法のフローチャートである。ここで、Ｓｉ（ｉ＝１，２，３，…）はステップを示す。従来、粗い格子と細かい格子の二種類の格子を用意し、各格子の計算値から格子収束解を推定する算出システムが実用化されている。
【０００５】
図９〜１１に示すように、一方において、ステップ１では粗い格子を生成し、次にステップ２において、この粗い格子を用いた計算を行い、ステップ３においてその計算結果Ａｃを出力し、その後ステップ７に移行する。他方、ステップ４では細かい格子を生成し、次にステップ５においてこの細かい格子を用いた計算を行い、ステップ６においてその計算結果Ａｆを出力し、その後ステップ７に移行する。ステップ７において、計算結果Ａｃと計算結果Ａｆの差の絶対値が十分に小さければ、ステップ８において計算結果Ａｃは格子収束解に十分に近似していると判定できる。ステップ７で計算結果Ａｃと計算結果Ａｆの差の絶対値が十分に小さくないと判定されたとき、計算結果Ａｃと計算結果Ａｆの差の絶対値は、誤差予測値であり、ステップ９においてこの誤差予測値にもとづいて、計算結果Ａｃと計算結果Ａｆからの格子収束解を予測する。
【０００６】
【発明が解決しようとする課題】
このような従来技術では、複雑な形状の計算空間に対して、複数の格子（粗い格子と細かい格子）を生成するのは、余分な作業時間が必要になる。また、細かい格子を用いるにしても、コンピュータの演算処理能力の制限からあまり細かい格子を用いることができないし、粗い格子を用いるにしても、粗すぎる格子では、計算精度が極端におちるので、現実的には計算格子のレンジをそれほど大きく設定できない。それ故、前記従来の算出システムを用いて格子収束解を予測するのは、精度が低く実質的に困難な場合があり、複雑な形状の計算空間に対しては限界がある。
【０００７】
したがって本発明の目的は、高精度の格子収束解を、従来よりも少ない計算量で得ることができる格子収束解の算出システムを提供することである。
【０００８】
【課題を解決するための手段】
請求項１記載の本発明は、演算手段の演算処理能力に見合う格子数を入力手段によって入力することで、前記演算手段によって、計算空間を有限個数の格子点で分割し、航空機や翼に働く揚力、横力、抵抗およびモーメントの流体運動を支配する偏微分方程式を離散化した差分方程式の格子収束解を算出し、その算出結果を出力手段によって出力する格子収束解の算出システムにおいて、
前記入力手段、演算手段および出力手段を含むコンピュータを、
単一の格子を用いた高次精度の計算手法による計算値と、前記単一の格子を用いた低次精度の計算手法による計算値とを求める手段、および
高次精度の計算値と低次精度の計算値と格子収束解との比が、いずれの格子数においても略一定として、前記単一の格子における高次精度の計算値と低次精度の計算値との関係から格子収束解を算出して表示する手段として機能させるための格子収束解の算出システムである。
【０００９】
本発明に従えば、演算手段の演算処理能力に見合う格子数を入力手段によって入力することで、前記演算手段によって、計算空間を有限個数の格子点で分割し、航空機や翼に働く揚力、横力、抵抗およびモーメントの流体運動を支配する偏微分方程式を離散化した差分方程式の格子収束解を算出し、その算出結果を出力手段によって出力する格子収束解の算出システムにおいて、ある単一の格子における高次精度の計算手法による計算値を求めるとともに、前記と同一の格子における低次精度の計算手法による計算値を求めておき、これら高次精度の計算値と低次精度の計算値と格子収束解との比が、いずれの格子数においても略一定として、前記高次精度の計算値と低次精度の計算値との関係から格子収束解を算出して表示することができる。
【００１４】
【発明の実施の形態】
本実施形態は、格子収束解の算出システムを、回転翼航空機のロータ・ブレードなどの翼断面を決定するために用いるコンピュータ２に適用した場合の一例である。図１は本発明の一実施形態に係るコンピュータ２のブロック図である。コンピュータ２は、主に、制御装置３と、キーボード４およびマウス５と、ディスプレイ６およびプリンタ７などで構成されている。制御装置３は、ＣＰＵ８とＲＯＭ９とＲＡＭ１０とから成るマイクロコンピュータと、バス１１と、入出力インタフェース１２と、駆動回路１３，１４などから成り、ＣＰＵ８とＲＯＭ９とＲＡＭ１０とはバス１１を介して入出力インタフェース１２に接続されている。
【００１５】
入出力インタフェース１２には、キーボード４およびマウス５がそれぞれ接続されている。また、入出力インタフェース１２には、駆動回路１３，１４を介してディスプレイ６、プリンタ７がそれぞれ接続されている。後述する格子収束解の算出プログラム１５は、図示外の磁気媒体などを介してＲＯＭ９に格納されている（図２参照）。キーボード４およびマウス５からの入力データの入力操作などによって、算出プログラム１５を実行（演算処理）し、ディスプレイ６またはプリンタ７に出力表示するようになっている。
【００１６】
図２は、格子収束解の算出プログラム１５の構成を示す説明図である。図３は、格子生成プログラム１６をコンピュータ２によって実行して生成した格子Ｋをディスプレイ６に表示した場合の説明図である。算出プログラム１５は、格子生成プログラム１６と流れ解析プログラム１７とを有し、格子生成プログラム１６においては、ＣＰＵ８の演算処理能力に見合う実用的な格子数と、翼型の入力データを、キーボード４およびマウス５から入力した後、この格子生成プログラム１６を実行し所望の格子Ｋを生成する。格子Ｋは、複数の放射線Ｌｈと複数の環状線Ｌｋとで略矩形状に形成され、複数の放射線Ｌｈは、翼の断面の外形線から略垂直方向外方向きに放射状に離隔し、複数の環状線Ｌｋは、前記外形線から径方向外方向きに所定間隔おきに形成され、複数の放射線Ｌｈに交差する。流れ解析プログラム１７においては、流れの諸条件として揚力係数などの入力データを、キーボード４およびマウス５から入力した後、前記格子生成プログラム１６で生成した格子を用いて、たとえば翼の中央付近の代表的２次元断面の抗力係数ＣＤの解析を行う。
【００１７】
この流れ解析プログラム１７は、非圧縮性粘性流の流れを偏微分方程式によって表すナビエ・ストークスの方程式を、差分法によって離散化し、この離散化した差分方程式の格子収束解（たとえば所望の抗力係数ＣＤ）を算出する機能を有する。コンピュータ２の計算空間においては二次元翼として推定計算する。
【００１８】
図４は、格子収束解Ａｃｏｎｖと誤差予測値ＥＲＲの計算法を示すフローチャートである。ここで、Ｓｉ（ｉ＝１０，１１，１２，…）はステップを示す。図５は、高次精度の計算手法による計算値と、低次精度の計算手法による計算値と、格子収束解との関係を示す図表である。本実施形態における高次精度の計算手法は、たとえば、微分方程式の微分をテーラー展開の３次まで正確にとった３次精度法であり、低次精度の計算手法は、同じく２次精度法である。
【００１９】
図４に示すように、ステップ１０において、単一の格子である前記実用的な格子を生成した後、ステップ１１とステップ１３が並列処理される。一方のステップ１１において低次精度の計算手法が実行され、次にステップ１２に移行して低次精度の計算値Ａｌｏｗを求めた後、ステップ１５に移行する。他方、ステップ１３において高次精度の計算手法が実行され、次にステップ１４に移行して高次精度の計算値Ａｈｉｇｈを求めた後、ステップ１５に合流する。ステップ１５において、格子収束解Ａｃｏｎｖおよび誤差予測値ＥＲＲを以下の式（１）および式（２）を用いて算出する。
【００２０】
Ａｃｏｎｖ＝Ａｈｉｇｈ＋α（Ａｈｉｇｈ−Ａｌｏｗ）      …（１）
ただし、０．１≦α≦１０．０とする。このαは、計算手法に応じて予め定められる第１の定数αであり、代表的な値はα＝１．０である。
【００２１】
ＥＲＲ＝β｜Ａｈｉｇｈ−Ａｌｏｗ｜                      …（２）
ただし、０．１≦β≦１０．０とする。このβは、計算手法に応じて予め定められる第２の定数βであり、代表的な値はβ＝１．０である。これら第１，第２の定数α，βは、いくつかの格子を用いた計算結果から最小自乗法および標準偏差の統計操作により予め算出しておく。具体的には、第１の定数αは、以下の標準的な最小自乗法の式（３）を用いて表すことができる。ここでＮ種類の異なる細かさの格子を、モデル問題について用意し、ｉ番目（１≦ｉ≦Ｎ）の格子に対して、高次精度と低次精度の計算を実施する。高次精度の計算結果をＡｈｉｇｈ（ｉ）とし、低次精度の計算結果をＡｌｏｗ（ｉ）とする。
【００２２】
【数１】
0003621378.tif 000002
【００２３】
式（３）についてＡｆと第１の定数αを変数とし、この式（３）についてｄε／ｄα＝０とｄε／ｄＡｆ＝０をそれぞれ計算すると、第１の定数αと変数Ａｆについての線形連立方程式が得られ、それらを解いて第１の定数αと変数Ａｆとを求める。第２の定数βは標準偏差であり、以下の式（４）により求める。この第２の定数βは誤差の許容範囲に応じて、０．１〜１０の範囲で定数倍する。
【００２４】
【数２】
0003621378.tif 000003
【００２５】
ところで、格子収束解Ａｃｏｎｖを求める式（１）においては、図５に示すように、高次精度の計算値Ａｈｉｇｈと低次精度の計算値Ａｌｏｗと格子収束解Ａｃｏｎｖとの関係が、格子の粗密に拘わらず｛ある任意の同一格子（単一の格子）に対して｝略一定である。したがって、Ａｈｉｇｈ：Ａｌｏｗ：Ａｃｏｎｖの比は、いずれの格子数においても略一定である。また、高次精度の計算値Ａｈｉｇｈが低次精度の計算値Ａｌｏｗに比べて格子収束解Ａｃｏｎｖに早く漸近するので、式（１）の関係が必然的に成立する。すなわち、高次精度の計算値Ａｈｉｇｈから低次精度の計算値Ａｌｏｗを減じた差を求め、この差に前記第１の定数αを乗じた値に、高次精度の計算値Ａｈｉｇｈを加えることによって、格子収束解Ａｃｏｎｖを算出することができる。この求められた格子収束解Ａｃｏｎｖは、入出力インタフェース１２および駆動回路１３，１４を介してディスプレイ６またはプリンタ７に出力表示される。またこの格子収束解ＡｃｏｎｖはＲＡＭ１０の図示外のワークメモリに記憶され、必要に応じてディスプレイ６またはプリンタ７に出力表示される。
【００２６】
誤差予測値ＥＲＲを求める式においては、高次精度の計算値Ａｈｉｇｈと、低次精度の計算値Ａｌｏｗの差を求め、この差の絶対値を算出し、この算出した絶対値に前記第２の定数βを乗じて、格子収束解Ａｃｏｎｖの誤差予測値ＥＲＲを算出することができる。この求められた誤差予測値ＥＲＲは、前記と同様に、入出力インタフェース１２および駆動回路１３，１４を介してディスプレイ６またはプリンタ７に出力表示される。またこの誤差予測値ＥＲＲはＲＡＭ１０のワークメモリに記憶され、必要に応じてディスプレイ６またはプリンタ７に出力表示される。
【００２７】
図６は、前記算出プログラム１５を用いて算出した、格子数と、翼の中央部付近の代表的２次元断面における抗力係数ＣＤとの関係を示す図表である。入力データとしては、翼型ＮＡＣＡ００１２と、揚力係数ＣＬ＝０．２およびマッハ数Ｍ＝０．７５などである。この抗力係数ＣＤの格子収束解を算出して表示するには、第１の定数αおよび第２の定数βを予め定めておき、先ず、上述した格子生成プログラム１６において、ＣＰＵ８の演算処理能力に見合う実用的な格子数と、翼型ＮＡＣＡ００１２の入力データを、キーボード４およびマウス５から入力した後、この格子生成プログラム１６を実行し所望の格子を生成する。
【００２８】
流れ解析プログラム１７においては、揚力係数ＣＬ＝０．２およびマッハ数Ｍ＝０．７５などの入力データを、キーボード４およびマウス５から入力した後、前記格子生成プログラム１６で生成した格子を用いて抗力係数ＣＤの格子収束解を算出する。この求められた抗力係数ＣＤの格子収束解は、入出力インタフェース１２および駆動回路１３，１４を介してディスプレイ６またはプリンタ７に出力表示される。また、この格子収束解はＲＡＭ１０のワークメモリに記憶され、必要に応じてディスプレイ６またはプリンタ７に出力表示される。
【００２９】
その後、前記出力表示された抗力係数ＣＤの格子収束解およびその他の格子収束解に基いて、たとえば風洞試験用モデルを形成し、前記格子収束解と風洞試験結果とを比較検討し、この算出システムによって得た格子収束解と、風洞試験による実験値とが略一致することを確認する。
【００３０】
以上説明した格子収束解Ａｃｏｎｖの算出システムによれば、ある単一の格子における高次精度の計算手法による計算値Ａｈｉｇｈを求めるとともに、前記と同一の格子における低次精度の計算手法による計算値Ａｌｏｗを求めておき、これら高次精度の計算値Ａｈｉｇｈと低次精度の計算値Ａｈｉｇｈと格子収束解Ａｃｏｎｖとの関係が、格子の粗密に拘わらず略一定として、高次精度の計算値Ａｈｉｇｈと低次精度の計算値Ａｌｏｗとの関係から格子収束解Ａｃｏｎｖを算出して表示することができる。
【００３１】
それ故、前記従来の算出システムで用いる計算格子よりも格段に少ない計算格子でもって高精度の格子収束解Ａｃｏｎｖを得ることができるうえ、複雑な形状の計算空間に対して、粗さの異なる複数種類の格子を用意する必要もなくなり、余分な作業時間を削減することができる。具体的には、従来の算出システムに比べて数十分の一の計算量で高精度の格子収束解Ａｃｏｎｖを得ることができる。
【００３２】
したがって、設計開発期間を大幅に短縮することができる。
格子収束解Ａｃｏｎｖを実際に算出して表示するには、第１の定数αを予め求めておき、高次精度の計算値Ａｈｉｇｈから低次精度の計算値Ａｌｏｗを減じた差を求め、この差に前記第１の定数αを乗じた値に高次精度の計算値Ａｈｉｇｈを加え、格子収束解Ａｃｏｎｖを算出し、その後、ディスプレイ６またはプリンタ７に出力表示することができる。
【００３３】
格子収束解Ａｃｏｎｖの誤差予測値ＥＲＲを算出して表示するには、第２の定数βを予め求めておき、高次精度の計算値Ａｈｉｇｈと、低次精度の計算値Ａｌｏｗの差を求め、この差の絶対値を算出し、この絶対値に第２の定数βを乗じ、誤差予測値ＥＲＲを算出し、その後、ディスプレイ６またはプリンタ７に出力表示することができる。
【００３４】
図７は、航空機の翼胴の揚抗特性を示し、格子収束解Ａｃｏｎｖと実験値とを比較する図表である。本発明の実施の他の形態として、格子収束解の算出システムを、航空機の翼胴の揚抗特性予測に適用する場合もある。これら計算値と実験値とでは、層流／乱流遷移の条件が異なるため、この実施形態においては、実験値に遷移補正を施してある。中程度の揚力範囲（揚力係数ＣＬが約０．２〜０．６の範囲）では、計算が二次精度から三次精度になると実験値に近づき、本発明の格子収束解の算出システムを適用すると、前記中程度の揚力範囲においては、実験値と略一致する格子収束解が得られることが判る。
【００３５】
前記第１，第２の定数α，βの代りに、これら定数と実質的に同等な±１０％以内の変動を含む関数を適用してもよい。本実施形態においては、コンピュータ２の計算空間において二次元翼として推定計算しているが、三次元翼として推定計算することも可能である。前記格子収束解の算出プログラム１５をスーパーコンピュータの制御装置に格納するとともに、このスーパーコンピュータを複数のコンピュータに接続線を介して接続し、各コンピュータから入力した入力データの入力信号に基いて、前記算出プログラム１５を実行し、格子収束解を、入力したコンピュータ側のディスプレイまたはプリンタで出力表示することも可能である。その他、前記実施形態に、特許請求の範囲を逸脱しない範囲において種々の部分的変更を行う場合もある。
【００３６】
【発明の効果】
請求項１記載の本発明によれば、ある単一の格子における高次精度の計算手法による計算値を求めるとともに、前記と同一の格子における低次精度の計算手法による計算値を求めておき、これら高次精度の計算値と低次精度の計算値と格子収束解との比が、いずれの格子数においても略一定として、前記高次精度の計算値と低次精度の計算値との関係から格子収束解を算出して表示することができる。
【００３７】
それ故、前記従来の算出システムで用いる計算格子よりも格段に少ない計算格子でもって高精度の格子収束解を得ることができるうえ、複雑な形状の計算空間に対して、粗さの異なる複数種類の格子を用意する必要もなくなり、余分な作業時間を削減することができる。
【図面の簡単な説明】
【図１】本発明の一実施形態に係るコンピュータのブロック図である。
【図２】格子収束解の算出プログラムの構成を示す説明図である。
【図３】格子生成プログラムで生成した格子をディスプレイ６に表示した場合の説明図である。
【図４】格子収束解と誤差予測値の計算法を示すフローチャートである。
【図５】高次精度の計算値と、低次精度の計算値と、格子収束解との関係を示す図表である。
【図６】格子数と、翼の中央付近の代表的２次元断面における抗力係数との関係を示す図表である。
【図７】翼胴の揚抗特性を示し、格子収束解と実験値とを比較する図表である。
【図８】従来の格子数と翼の中央付近の代表的２次元断面における抗力係数との関係を示す図表である。
【図９】翼形流れ計算用の粗い格子をディスプレイに表示した場合の説明図である。
【図１０】翼形流れ計算用の細かい格子をディスプレイに表示した場合の説明図である。
【図１１】従来の格子収束解および誤差の計算法のフローチャートである。
【符号の説明】
２  コンピュータ
３  制御装置
６  ディスプレイ
７  プリンタ
８  ＣＰＵ
９  ＲＯＭ
１０  ＲＡＭ
１５  算出プログラム
１６  格子生成プログラム
１７  流れ解析プログラム
Ａｃｏｎｖ  格子収束解
ＥＲＲ  誤差予測値
Ａｈｉｇｈ  高次精度の計算値
Ａｌｏｗ  低次精度の計算値
α  第１の定数
β  第２の定数
【図１】
0003621378.tif 000004
【図２】
0003621378.tif 000005
【図３】
0003621378.tif 000006
【図４】
0003621378.tif 000007
【図５】
0003621378.tif 000008
【図６】
0003621378.tif 000009
【図７】
0003621378.tif 000010
【図８】
0003621378.tif 000011
【図９】
0003621378.tif 000012
【図１０】
0003621378.tif 000013
【図１１】
0003621378.tif 000014