(57)【特許請求の範囲】
【請求項１】
翼面近傍の渦度に漸近する流れ方向の速度勾配絶対値Ｓを生成項とする乱流モデルにおいて、２つの変数ａ，ｂの大なる方を値とする関数をＭａｘ（ａ，ｂ）とし、翼面に沿う速度の翼面に沿う方向の速度勾配をＵｘとし、壁面に垂直な方向の速度勾配をＵｙとし、Ｕｙに漸近する渦度値等をΩとし、定数をαとし、座標軸Ｘは翼に対する気流の相対的な速度Ｕの方向を正としたとき、
Ｓ＝Ｍａｘ（｜Ω｜−α｜Ｕｘ｜，０）
または、
Ｓ＝Ｍａｘ（｜Ω｜−αＭａｘ（０，Ｕｘ），０）
によって修正された速度勾配絶対値Ｓを用いて、翼の迎角に対する揚力を求めることを特徴とする翼形性能の推定方法。
【請求項２】
壁面に垂直な方向の速度をＶとしたとき、その垂直方向の勾配はＶｙで表わし、流体の性質から、壁面近傍では、
Ｕｘ＋Ｖｙ＝０
が概ね成立することを利用し、前記翼面に沿う方向の速度勾配Ｕｘを翼面に垂直な方向の速度勾配−Ｖｙで置換えることを特徴とする請求項１記載の翼形性能の推定方法。
【請求項３】
前記定数αは、複数の異なる既知の翼形の最大揚力の計算値と実験値との相関が最良になる値に選ばれることを特徴とする請求項１または２記載の翼形性能の推定方法。
【請求項４】
前記定数αは、気流の流れのレイノルズ数に応じて異なる値に選ばれることを特徴とする請求項１，２または３記載の翼形性能の推定方法。
【請求項５】
前記定数αは、複数のレイノルズ数に基づいて求めた値を、簡易式で近似した値に置換えることを特徴とする請求項１，２または３記載の翼形性能の推定方法。
【請求項６】
翼面近傍の渦度に漸近する速度勾配絶対値Ｓを生成項とする乱流モデルを用いて翼形性能を求めるコンピュータを、
翼形性能を生成するために必要なデータを入力するための入力手段、
入力手段によって入力されたデータに基づいて、２つの変数ａ，ｂの大なる方を値とする関数をＭａｘ（ａ，ｂ）とし、翼面に沿う速度の翼面に沿う方向の速度勾配をＵｘとし、壁面に垂直な方向の速度勾配をＵｙとし、Ｕｙに漸近する渦度値等をΩとし、定数をαとし、座標軸Ｘは翼に対する気流の相対的な速度Ｕの方向を正としたとき、
Ｓ＝Ｍａｘ（｜Ω｜−α｜Ｕｘ｜，０）
または、
Ｓ＝Ｍａｘ（｜Ω｜−αＭａｘ（０，Ｕｘ），０）
によって修正された速度勾配絶対値Ｓを用いて、翼の迎角に対する揚力を求める演算手段、および
演算手段によって求められた前記翼の迎角に対する揚力を表示する出力手段、として機能させるコンピュータによって実行可能な翼形性能の推定プログラム。
【請求項７】
壁面に垂直な方向の速度をＶとしたとき、その垂直方向の勾配はＶｙで表わし、流体の性質から、壁面近傍では、
Ｕｘ＋Ｖｙ＝０
が概ね成立することを利用し、前記翼面に沿う方向の速度勾配Ｕｘを翼面に垂直な方向の速度勾配−Ｖｙで置換えることを特徴とする請求項６記載のコンピュータによって実行可能な翼形性能の推定プログラム。
【請求項８】
前記定数αは、複数の異なる既知の翼形の最大揚力の計算値と実験値の相関が最良になる値に選ばれることを特徴とする請求項６または７記載のコンピュータによって実行可能な翼形性能の推定プログラム。
【請求項９】
前記定数αは、気流の流れのレイノルズ数に応じて異なる値に選ばれることを特徴とする請求項６，７または８記載のコンピュータによって実行可能な翼形性能の推定プログラム。
【請求項１０】
前記定数αは、複数のレイノルズ数に基づいて求めた値を、簡易式で近似した値に置換えることを特徴とする請求項６，７または８記載のコンピュータによって実行可能な翼形性能の推定プログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、回転翼航空機のロータブレードの翼形の決定および翼形決定に関与する最大揚力係数などの翼特性の確認などを行うための翼形性能を、高精度で推定することができる翼形性能の推定方法およびコンピュータによって実行可能な翼形性能の推定プログラムに関する。
【０００２】
本発明において、「翼型」とは、米国航空評議委員会（Ｎａｔｉｏｎａｌ  Ａｄｖｉｓｏｒｙ
Ｃｏｍｍｉｔｔｅｅ  ｆｏｒ  Ａｅｒｏｎａｕｔｉｃｓ；略称ＮＡＣＡ）で規定・分類された２次元の翼断面形状をいうものとし、この「翼型」を含む任意の翼断面形状をいう場合には、「翼形」と記す（以下同じ）。
【０００３】
【従来の技術】
航空機の翼（３次元形状をいう）および翼形（２次元断面形状をいう）は、航空機の性能に決定的な影響をもっており、翼形を含む翼については次のようなことが知られている。すなわち、翼は気流に対する角度である迎角が増大するにつれて、揚力が徐々に大きくなるが、ある角度で最大揚力を発生し、その後、揚力の低下および抗力の急増を引き起こす失速現象を生じる。この最大揚力およびその発生角度は、翼の重要な特性のひとつである。
【０００４】
このような翼の特性は、風洞実験によって測定されることが多いが、この風洞実験では実の飛行状態を模擬することが難しく、多くの時間およびコストを要するという問題があり、近年のコンピュータの急速な進歩に応じて、数値シミュレーションによる流れ解析によって、翼特性の測定などが多く行われている。この流れの数値シミュレーションは、数値流体力学（略称ＣＦＤ；Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｆｌｕｉｄ  Ｄｙｎａｍｉｃｓ）が用いられ、航空機の空力設計における開発経費の削減および期間の短縮に大きく貢献でき、ナビエ・ストークス解析によって、０°〜１３°程度の低迎角では遷音速特性、および乱流遷移の問題を除く抗力の推算も、高い精度で可能であるが、レイノルズ数Ｒｅが１０８にも達する流れ場における最大揚力性能の推定を行うには、現在のコンピュータの処理能力では、乱流現象を近似した有限体積法３次元レイノルズ平均Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ（Ｒｅｙｎｏｌｄｓ  Ａｖｅｒａｇｅｄ  Ｎａｖｉｅｒ
−Ｓｔｏｋｅｓ；略称ＲＡＮＳ）解析を用いるしかなく、乱流モデルの限界によって、最大揚力の推定精度は制限されている。
【０００５】
図５は従来の技術の翼の性能を推定するための概略的なアルゴリズムを示すフローチャートである。たとえばヘリコプタのロータブレードの翼断面の性能を、コンピュータを用いてシミュレーションするにあたっては、ステップａ１で、コンピュータによって翼特性解析プログラムが実行され、ステップａ２で、翼形および流れなどに関する諸条件を入力して初期設定すると、ステップａ２で、入力された諸条件に基づいて乱流粘性係数の計算が開始される。
【０００６】
この乱流粘性係数は、これを数値的に解くために、連続的な偏微分方程式を離散的な格子点上で定義される差分方程式として書き直す必用がある。流れ場の計算格子は、流速ｕの方向をＸ方向とし、動粘性係数νの方向をＹ方向とするＸＹ直交座標系上で、流れ場をＸ方向およびＹ方向の長さがΔｘ，Δｙの格子で分割し、格子点番号（ｉ，ｊ）（ここに、ｉ，ｊは正の整数）で表される各格子点で囲まれた格子中央に圧力ｐを定義する。時間ステップをΔｔとしたとき、時刻ｔ＝ｔ０＋ｎΔｔ（ここに、ｎは正の整数）における各格子点位置の渦度ω（ｘ，ｙ；ｔ）と流れ関数ψ（ｘ，ｙ；ｔ）とを導入することによって、流れを規定する２次元ナビエ・ストークス方程式と連続の式とを、渦度輸送方程式とポアソン方程式とに変形し、それらを差分化する。
【０００７】
次に、ステップａ４で、流れ場の計算が開始され、前記差分化によって得られた差分方程式に対して、時刻ｔに関する計算を全格子点について行い、ステップａ５で流れ関数ψが収束するまで上記ステップａ３，ａ４を繰り返し、ステップａ５で定常解に収束すると、次のステップａ６に移り、迎角に対する揚力係数などの翼形状特性が求められ、ステップａ７で計算が終了する。
【０００８】
前記ステップａ３において、渦度｜Ω｜は、｜Ｕｙ−Ｖｘ｜によって計算される。ここで、ＵｙはＹ方向の速度勾配を示す偏微分であり、ＶｘはＸ方向の速度勾配を示す偏微分であり、Ｕ，ＶはＸ，Ｙ方向の速度である。Ｕｙ，Ｖｘは差分近似によって計算される。このような各速度勾配Ｕｙ，Ｖｘの近似法は、従来から多くの手法が周知である。流れの計算方法には、上記のように渦度ωと流れ関数ψとを用いるω‐ψ法以外に、流速ｕ，ｖと圧力ｐとを未知変数として用いるｕ−ｐ法がある。さらに翼面に平行な速度勾配を次式、
｛２Ｕｘ２＋Ｕｙ２＋２Ｕｙ・Ｖｘ＋Ｖｘ２−２／３・（Ｕｘ＋Ｖｙ）２｝１／２
によって求めることもできる。
【０００９】
このような各種の手法のいずれによって渦度Ω等を求めても、翼面近傍の境界層では、翼面に平行な速度成分だけが残り、翼面に垂直な方向の変化量が非常に少なくなるので、結局は同じ値に収束して、壁面がＸ軸に平行な場合、Ｙ方向の速度勾配と同等になり、したがって渦度Ωも同等になることが知られているが、約１３°〜１５°の高迎角における最大揚力の定常的予測は困難であり、翼の空力設計を高精度でシミュレートすることができないという問題がある。
【００１０】
【発明が解決しようとする課題】
本発明の目的は、翼の最大揚力付近の特性を高精度で予測することができる翼形性能の推定方法および翼形性能の推定プログラムを提供することである。
【００１１】
【課題を解決するための手段】
請求項１記載の本発明は、翼面近傍の渦度に漸近する流れ方向の速度勾配絶対値Ｓを生成項とする乱流モデルにおいて、２つの変数ａ，ｂの大なる方を値とする関数をＭａｘ（ａ，ｂ）とし、翼面に沿う速度の翼面に沿う方向の速度勾配をＵｘとし、壁面に垂直な方向の速度勾配をＵｙとし、Ｕｙに漸近する渦度値等をΩとし、定数をαとし、座標軸Ｘは翼に対する気流の相対的な速度Ｕの方向を正としたとき、
Ｓ＝Ｍａｘ（｜Ω｜−α｜Ｕｘ｜，０）
または、
Ｓ＝Ｍａｘ（｜Ω｜−αＭａｘ（０，Ｕｘ），０）
によって修正された速度勾配絶対値Ｓを用いて、翼の迎角に対する揚力を求めることを特徴とする翼形性能の推定方法である。
【００１２】
本発明に従えば、翼面近傍の渦度に対して、壁面に沿う方向Ｘの速度勾配Ｓを、Ｓ＝｜Ω｜−α｜Ｕｘ｜またはＳ＝｜Ω｜−αＭａｘ（０，Ｕｘ）によって修正するので、高迎角においても速度勾配絶対値Ｓを定常解に収束させることができ、最大揚力および失速角などの翼特性を正確に予測することが可能となり、翼の最大揚力付近の特性を、たとえばコンピュータを用いた数値シミュレーションによって高精度で模擬することが可能となる。
【００１３】
請求項２記載の本発明は、翼近傍では流体の性質から、壁面に対して垂直方向速度の垂直方向勾配をＶｙとしたとき、
Ｕｘ＋Ｖｙ＝０
が近似的に成立することを用いて、前記翼面に沿う方向の速度勾配Ｕｘを翼面に垂直な方向の速度勾配−Ｖｙで置換えることを特徴とする。
【００１４】
本発明に従えば、２次元断面では壁面に沿う方向の速度の壁面に沿う方向の速度勾配Ｕｘの計算は容易であるが、３次元の翼で流れが翼を斜めに横切る場合には、この算出は容易ではない。一方、壁面に対して垂直方向速度の垂直方向勾配Ｖｙは、この場合でも算出が容易であるという利点がある。
【００１５】
請求項３記載の本発明は、前記定数αは、複数の異なる既知の翼形の最大揚力の計算値と実験値との相関が最良になる値に選ばれることを特徴とする。
【００１６】
本発明に従えば、定数αが複数の異なる翼形の最大揚力係数が実験値と最も近くなる値に選ばれるので、高い信頼性で高迎角における速度勾配絶対値Ｓが求められ、最大揚力および失速角などの翼特性の予測精度が向上され、翼設計の高精度化が図られるとともに、翼断面決定に要する時間の短縮を図ることができる。
【００１７】
請求項４記載の本発明は、前記定数αは、気流の流れのレイノルズ数に応じて異なる値に選ばれることを特徴とする。
【００１８】
本発明に従えば、流れの状態は、流れの代表速度をＵ、流れの代表長さをＬ、動粘度をνとしたとき、レイノルズ数Ｒｅ＝ＵＬ／ν、すなわちＲｅ＝（慣性力）／（粘性力）で表され、このレイノルズ数Ｒｅに応じて定数αの値を変化させることによって、流れの慣性力と粘性力との相対的影響度を反映させて、前記速度勾配の絶対値Ｓが求められ、最大揚力および失速角などの翼特性を正確に予測することが可能となる。
【００１９】
請求項５記載の本発明は、前記定数αは、複数のレイノルズ数に基づいて求めた値を、簡易式で近似した値に置換えることを特徴とする。
【００２０】
本発明に従えば、複数のレイノルズ数から求めた値を簡易式で近似した値に置換えた定数αを用いて、前記速度勾配の絶対値Ｓを求めるので、高迎角の最大揚力をより高い信頼性で求めることができる。簡易式の例としては、２つのレイノルズ数Ｒｅ１，Ｒｅ２において得られた定数をα１，α２とし、任意のレイノルズ数Ｒｎでの定数αを、次式
α＝α１＋（Ｒｎ−Ｒｅ１）（α２−α１）／（Ｒｅ２−Ｒｅ１）
で与える方法がある。
【００２１】
請求項６記載の本発明は、翼面近傍の渦度に漸近する速度勾配絶対値Ｓを生成項とする乱流モデルを用いて翼形性能を求めるコンピュータを、
翼形性能を生成するために必要なデータを入力するための入力手段、
入力手段によって入力されたデータに基づいて、２つの変数ａ，ｂの大なる方を値とする関数をＭａｘ（ａ，ｂ）とし、翼面に沿う速度の翼面に沿う方向の速度勾配をＵｘとし、壁面に垂直な方向の速度勾配をＵｙとし、Ｕｙに漸近する渦度値等をΩとし、定数をαとし、座標軸Ｘは翼に対する気流の相対的な速度Ｕの方向を正としたとき、
Ｓ＝Ｍａｘ（｜Ω｜−α｜Ｕｘ｜，０）
または、
Ｓ＝Ｍａｘ（｜Ω｜−αＭａｘ（０，Ｕｘ），０）
によって修正された速度勾配絶対値Ｓを用いて、翼の迎角に対する揚力を求める演算手段、および
演算手段によって求められた前記翼の迎角に対する揚力を表示する出力手段、として機能させるコンピュータによって実行可能な翼形性能の推定プログラムである。
【００２２】
本発明に従えば、コンピュータによって上記翼形性能の推定プログラムが実行され、入力手段によって翼形性能を生成するために必用な、たとえば翼断面形状、流れの条件を特定するための諸量、格子点などの各種のデータが入力されると、演算手段は前記入力手段から入力されたデータに基づいて、翼面近傍の渦度に対して、壁面に沿う方向Ｘの速度勾配絶対値Ｓを、Ｓ＝｜Ω｜−α｜Ｕｘ｜またはＳ＝｜Ω｜−αＭａｘ（０，Ｕｘ）によって修正し、これによって高迎角においても速度勾配絶対値Ｓを定常解に収束させて、出力手段に翼の迎角に対する揚力を表示させることができる。このようにして出力手段に翼の迎角と揚力とが表示されるので、最大揚力および失速角などの翼特性を正確に予測することが可能となり、翼の最大揚力付近の特性を高精度で模擬することが可能な翼形性能の解析ツールをコンピュータ上に実現することができる。
【００２３】
請求項７記載の本発明は、翼近傍では流体の性質から、壁面に対して垂直方向速度の垂直方向勾配をＶｙとしたとき、
Ｕｘ＋Ｖｙ＝０
が近似的に成立することを用いて、前記翼面に沿う方向の速度勾配Ｕｘを翼面に垂直な方向の速度勾配−Ｖｙで置換えることを特徴とする。
【００２４】
本発明に従えば、２次元断面では壁面に沿う方向の速度の壁面に沿う方向の速度勾配Ｕｘの計算は容易であるが、３次元の翼で流れが翼を斜めに横切る場合には、この算出は容易ではない。一方、壁面に対して垂直方向速度の垂直方向勾配Ｖｙは、この場合でも算出が容易であるという利点がある。
【００２５】
請求項８記載の本発明は、前記定数αは、複数の異なる既知の翼形の最大揚力の計算値と実験値の相関が最良になる値に選ばれることを特徴とする。
【００２６】
本発明に従えば、定数αが複数の異なる翼形の最大揚力係数が実験値と最も近くなる値に選ばれるので、高い信頼性で高迎角における速度勾配絶対値Ｓが求められ、最大揚力および失速角などの翼特性の予測精度が向上され、翼設計の高精度化が図られるとともに、翼断面決定に要する時間の短縮を図ることができる。
【００２７】
請求項９記載の本発明は、前記定数αは、気流の流れのレイノルズ数に応じて異なる値に選ばれることを特徴とする。
【００２８】
本発明に従えば、流れの状態は、流れの代表速度をＵ、流れの代表長さをＬ、動粘度をνとしたとき、レイノルズ数Ｒｅ＝ＵＬ／ν、すなわちＲｅ＝（慣性力）／（粘性力）で表され、このレイノルズ数Ｒｅに応じて定数αの値を変化させることによって、流れの慣性力と粘性力との相対的影響度を反映させて、前記速度勾配の絶対値Ｓが求められ、最大揚力および失速角などの翼特性を正確に予測することが可能となる。
【００２９】
請求項１０記載の本発明は、前記定数αは、複数のレイノルズ数に基づいて求めた値を、簡易式で近似した値に置換えることを特徴とする。
【００３０】
本発明に従えば、複数のレイノルズ数から求めた値を境界層近似した値に置換えた定数αを用いて、前記速度勾配の絶対値Ｓを求めるので、高迎角の最大揚力をより高い信頼性で求めることができる。
【００３１】
【発明の実施の形態】
図１は本発明の実施の一形態の翼形性能の推定方法を示すフローチャートであり、図２は翼形性能の推定プログラムを実行して翼形性能をシミュレーションするために用いられるコンピュータ１の電気的構成を示すブロック図である。本実施の形態では、ヘリコプタのロータブレードの翼断面の性能を、コンピュータを用いて推定する方法について説明する。ステップｂ１で、コンピュータ１によって翼形性能の推定プログラムが実行され、翼形および流れなどに関する諸条件が入力手段２によって入力されて、初期設定される。前記翼形性能の推定プログラムは、記憶手段３に記憶され、演算処理装置４からの実行指令によって読み出されて実行され、出力手段である表示手段５に計算結果を表示させ、数値流体解析システムを構築する。
【００３２】
このコンピュータ１による解析は、圧縮性有限体積法３次元レイノルズ平均ナビエ・ストークス（Ｒｅｙｎｏｌｄｓ  Ａｖｅｒａｇｅｄ  Ｎａｖｉｅｒ−Ｓｔｏｋｅｓ，略称ＲＡＮＳ）方程式を使用する。レイノルズ応力は乱流粘性として導入し、乱流モデルは後述のＳｐａｌａｎｔ−Ａｌｌｍａｒａｓ（略称ＳＡ）乱流モデルを用いる。計算手法は、有限体積法を用い、基礎変数にＶａｎ  Ａｌｂａｄａのリミタを適用したＭＵＳＣＬ（Ｍｏｎｏｔｏｎｉｃ  Ｕｐｗｉｎｄ  Ｓｃｈｅｍｅ  ｆｏｒ
Ｃｏｎｓｅｒｖａｔｉｏｎ  Ｌｏｗ）によって空間２次精度とする。近似Ｒｉｅｍａｎｎ流速としては、ＳＨＵＳ（Ｓｉｍｐｌｅ  Ｈｉｇｈ  ｒｅｓｏｌｕｔｉｏｎ  Ｕｐｗｉｎｄ  Ｓｃｈｅｍｅ）を使用する。また、時間積分には、ＭＦＧＳ（Ｍａｔｒｉｘ  Ｆｒｅｅ  Ｇａｕｓｓ  Ｓｅｉｄｅｌ）陰解法を用い、収束加速のために局所時間法を併用する。
【００３３】
この局所時間法によって収束解が得られない場合には、物理時間法による非定常計算に切り換える。この場合の時間精度は一次とし、時間刻みはコード長と音速とを基準にした無次元時間で０．０１とする。計算は、最小迎角から開始し、収束解を次の迎角の初期値として使用する。高迎角においては、実験においても流れの初期値依存性を生じることがあるので、風洞試験での迎角変化の模擬を念頭においている。
【００３４】
計算格子は、代数的手法によるＯ型格子で、全４角形格子を使用し、格子点数は周方向に２５６とし、垂直方向は６０とする。また翼弦長を基準として、垂直方向の最小格子間隔は１０−５とし、計算領域の半径は２０とし、周方向の最小格子間隔は前縁および後縁で約８×１０−４とし、６％翼の前縁については１．６×１０−４とする。流れの条件は、マッハ数が０．２、レイノルズ数が５．８×１０６の一様流とする。
【００３５】
こうしてステップｂ１で初期設定すると、ステップｂ２で、入力された諸条件に基づいて乱流粘性係数の計算が開始される。この乱流粘性係数を求めるにあたっては、これを数値的に解くために、連続的な偏微分方程式を離散的な格子点上で定義される差分方程式として書き直す必用がある。すなわち、ステップｂ２１で、流れ場の計算格子は、流速ｕの方向をＸ方向とし、動粘性係数νの方向をＹ方向とするＸＹ座標上で、流れ場をＸ方向およびＹ方向の長さがΔｘ，Δｙの格子で分割し、格子点番号（ｉ，ｊ）（ここに、ｉ，ｊは正の整数）で表される各格子点で囲まれた格子中央に圧力ｐを定義する。時間ステップをΔｔとしたとき、時刻ｔ＝ｔ０＋ｎΔｔ（ここに、ｎは正の整数）における各格子点位置の渦度ω（ｘ，ｙ；ｔ）と流れ関数ψ（ｘ，ｙ；ｔ）とを導入することによって、流れを規定する２次元ナビエ・ストークス方程式と連続の式とを、渦度輸送方程式とポアソン方程式とに変形し、それらを差分化して、Ｘ方向の速度Ｖｘおよびその速度勾配ＵｘならびにＹ方向の速度Ｖｙおよびその速度勾配Ｕｙを求める。ここで、ＵｙはＹ方向の速度勾配を示す偏微分であり、ＶｘはＸ方向の速度勾配を示す偏微分である。
【００３６】
上記のナビエ・ストークス方程式をコンピュータ１によって数値シミュレーションするには、何らかの形で近似化する必要があり、このような近似化の手法として、有限差分法、有限体積法、および有限要素法などを用いて、現実の連続な空間を、数多くの離散的な計算点で置換え、ナビエ・ストークス方程式を差分方程式などによって近似し、計算を実行する。最大揚力の予測に関しては、乱流モデルを用いる。
【００３７】
この乱流モデルには、多くの種類があるが、共通点として、何らかの形式で速度の空間勾配を用いていることである。ｘ，ｙ，ｚを３次元の空間座標、ｕ，ｖ，ｗを各方向Ｘ，Ｙ，Ｚの速度とすると、速度勾配は速度の各方向の座標での偏微分の記号によって、ｕｘ，ｕｙ，ｕｚ，νｘ，νｙ，νｚ，ｗｘ，ｗｙ，ｗｚで表される。
【００３８】
これらの式で表される偏微分は、ステップｂ２１において、数値計算上、差分で近似される。たとえば３次元の直交格子で空間を離散化し、その格子点の３次元座標が格子インデックス（ｉ，ｊ，ｋ）とＸ，Ｙ，Ｚの各方向の格子間隔Δｘ，Δｙ，Δｚとを用いて、次のように表される場合を考える。
（ｘｉ，ｊ，ｋ，ｙｉ，ｊ，ｋ，ｚｉ，ｊ，ｋ）＝（ｉΔｘ，ｊΔｙ，ｋΔｚ）
【００３９】
この格子上の各点は、インデックス（ｉ，ｊ，ｋ）で区別され、（ｉ，ｊ，ｋ）での速度はｕｉ，ｊ，ｋのように表される。この場合、インデックス（ｉ，ｊ，ｋ）における速度勾配ｕｘｉ，ｊ，ｋは、次のように近似することができる。
ｕｘ，ｉ，ｊ，ｋ＝（ｕｉ＋１，ｊ，ｋ−ｕｉ−１，ｊ，ｋ）／２Δｘ
【００４０】
数値流体力学（略称ＣＦＤ；Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｆｌｕｉｄ  Ｄｙｎａｍｉｃｓ）における格子およびメッシュの選び方、その上での差分近似の方法は多くの手法が周知であるが、どのような格子、メッシュを用い、どのような差分近似を用いても、格子を細かくしていけば、元の偏微分に限りなく近づく。したがって、乱流モデルの適用において、具体的な差分近似の手法は、その効果に対しては副次的な影響しかもたらさない。そのため、現実の適用にあたっては、それらを差分近似などで離散化したものを用いるが、本実施の形態では、差分近似表示ではなく、元の偏微分の記号を使って説明する。乱流モデルの適用は、このような計算を数万〜数百万点存在するすべての格子点に対して、数百〜数万ステップある全時間ステップについて、繰り返すことによって実現される。
【００４１】
乱流モデル中での速度勾配指標は、乱流モデルの効果に直接の影響があり、いくつかの計算手法が知られており、たとえば渦度ωと各速度勾配ｕｙ，ｕｚ，νｘ，νｚ，ｗｚとの関係を次式のように想定する次のような手法がある。
｜Ω｜＝｛（ｕｙ−νｘ）２＋（νｚ−ｗｘ）２＋（ｗｘ−ｕｚ）２｝１／２
【００４２】
航空機まわりの流れのような高レイノルズ数の流れでは、物体近傍に境界層と呼ばれるごく薄い、速度の遅い領域が発生する。物体表面では速度が０になるので、境界層の中で、壁に近づくにつれて速度は急速に減少する。このため、速度勾配の成分のうち壁に垂直な方向の変化が平行な方向に比べてはるかに大きくなることが知られている。また、壁の近傍では、流れは明らかに壁に平行な方向であるので、壁の近くでは壁に沿う速度成分の壁に垂直な方向の変化は、他の速度勾配に比べて大きくなることが知られている。仮に、壁がＸＹ平面に平行な場合、ＵｘとＷｙ成分の大きさ、すなわち絶対値は、他のＶｘ，Ｖｙ，Ｖｚ，Ｕｘ，Ｕｙ，Ｕｚ，Ｗｘ，Ｗｚよりも大きくなることが知られている。
【００４３】
上記のように、高レイノルズ数の流れでは、大きな速度勾配が境界層の内部に発生し、壁に沿う速度成分の壁に垂直な方向の変化が、他の成分に比べて大きくなるが、翼の最大揚力付近では、翼の前縁付近で壁に沿った方向の速度勾配がかなり大きな値となる。
【００４４】
乱流モデルについて説明する。最大揚力係数の予測精度に関して、本件発明者は数種類の乱流モデルをテストした結果、最良であったＳｐａｌａｒｔ−Ａｌｌｍａｒａｓモデル（ＳＡモデルと略記する）を採用する。テストした他の乱流モデルは、Ｂａｌｄｗｉｎ−Ｌｏｍａｘ代数モデル（略称ＢＬ）、Ｂａｌｄｗｉｎ−Ｂａｒｔｈモデル、ｑ−ωモデルである。ＳＡモデルは、基本的には乱流粘性係数に対するソース項付移流拡散方程式を解く一方程式モデルであり、次のように表される。なお、下式中の記載「ν〜」は外部領域で乱流粘性係数に漸近する乱流変数である。
【００４５】
0003587827.tif 000002
ここで、νは分子粘性係数、νｔは乱流粘性係数、Ｓは渦度Ωの絶対値（＝｜Ω｜）、各定数の値は、ｃｂ１＝０．１３５５，σ＝２／３，ｃｂ２＝０．６２２，ｋ＝０．４１，ｃｗ１＝ｃｂ１／ｋ２＋（１＋ｃｂ２）／σ，ｃｗ２＝０．３，ｃｗ３＝２，ｃｖ１＝７．１，ｃｔ１＝１，ｃｔ２＝２，ｃｔ３＝１．２，ｃｔ４＝０．５である。
【００４６】
上記のＳＡ乱流モデルについて、次のような変更を行っている。
変数変換に関しては、乱流方程式について非構造有限体積法で二次精度法を適用するには保存形式の方が扱いやすいため、変数をν￣からμ￣≡ρν￣に変換した等価な次式を用いる。
（μ￣）ｔ＋（ｕμ￣）ｘ＋（νμ￣）ｙ＝ρ［Ｒ．Ｈ．Ｓ．］
ただし、ρは密度、ｕはＸ方向の速度、νはＹ方向の速度であり、これらのρ，ｕ，νには質量保存側、
ρｔ＋（ρｕ）ｘ＋（ρν）ｙ＝０
が成立しているものとする。
【００４７】
空間二次精度化に関しては、乱流方程式には一次精度風上を用いる場合が多く、平板境界層の解析等では精度的に充分であるが、高迎角では格子解像度の影響が顕著になるので、本実施の形態ではナビエ・ストークス方程式に対するのと同様のＭＵＳＣＬを用いて二次精度化する。上記の高迎角において格子解像度の影響が顕著になる点については、本件発明者の研究によって既に確認されており、航空宇宙学会３０周年記念講演会講演集（１９９９年）、第６６頁〜第６８頁に「ＣＦＤによる最大揚力係数の予測について」として掲載されている。
【００４８】
全場乱流に関しては、ＳＡモデルには遷移点を制御するための遷移項（ｆｔ１）と、それ以外の遷移を回避するための安定化項（ｆｔ２）とが含まれているが、本実施の形態ではこれらの二項を零とし、全場乱流の条件で解析する。ただし、一様流での値は低く保たれており、淀み点から境界層内で急速に増大するので、結果的には数値的な遷移を用いていることになる。
【００４９】
次に、ステップｂ２２で、渦度の絶対値｜Ω｜、
｜Ω｜＝ＳＱＲＴ｜Ｕｙ−Ｖｘ｜
を計算し、ステップｂ２３で、翼面である壁に平行方向の速度勾配Ｖｔ、
Ｖｔ＝｜ＵｘＸｔ２＋（Ｖｘ＋Ｕｙ）ＸｔＹｔ＋ＶｙＹｔ２｜
を計算する。
【００５０】
ステップｂ２４で、渦度ωの修正計算を次式、
｜Ω｜←Ｍａｘ（｜Ω｜Ｕｙ−αＶｔ，０）
または、
｜Ω｜←Ｍａｘ（｜Ω｜−αＭａｘ（０，Ｕｘ），０）
によって行う。上式において、Ｍａｘ（０，Ｕｘ）は、２つの変数０，Ｕｘの大なる方を変数値として選択する関数を意味する。翼面に沿う速度の翼面に沿う方向の速度勾配をＵｘとし、定数をαとする。この定数αは、０．０１＜α＜１００に選ばれ、代表的には３．０に選ばれる。座標軸Ｘは、翼に対する気流の相対的な速度Ｕの方向を正とする。
【００５１】
ステップｂ２５で、上記の式３または式４によって修正された渦度の絶対値｜Ω｜を用いて、翼の迎角に対する揚力を求め、次に、ステップｂ２５で、前記差分化によって得られた差分方程式に対して、時刻ｔに関する計算を行い、ステップａ４で流れ関数ψが収束するまで上記ステップａ２，ａ３を繰り返し、ステップａ５で、迎角に対する乱流粘性係数などの翼特性が求められる。
【００５２】
図４は１２％翼形（ＮＡＣＡ６３１−０１２）の揚力ＣＬ、抵抗ＣＤ、モーメントＣＭの３分力特性の本件発明者による実験結果を示すグラフである。同図において、横軸は迎角であり、縦軸は揚力ＣＬ、抵抗ＣＤ、モーメントＣＭであり、本発明の実施例のシミュレーション結果は「ＳＡ改」と示され、従来手法による比較例のシミュレーション結果は「従来」と示され、従来手法による実験値は「Ｅｘｐ」と示される。迎角１３°以上では計算結果と実験結果とはやや異なるが、最大揚力は計算によってほぼ正しく予想されることが確認された。
【００５３】
図５は１８％翼形（ＮＡＣＡ６３３−０１８）の揚力ＣＬ、抵抗ＣＤ、モーメントＣＭの３分力特性の本件発明者による実験結果を示すグラフである。同図において、横軸は迎角であり、縦軸は揚力ＣＬ、抵抗ＣＤ、モーメントＣＭであり、本発明の実施例のシミュレーション結果は「ＳＡ改」と示され、従来手法による比較例のシミュレーション結果は「従来」と示され、従来手法による実験値は「Ｅｘｐ」と示される。図４に示される翼（ＮＡＣＡ６３１−０１２）と同様に、最大揚力係数の計算値が実験値に近づいており、高迎角においても高精度で最大揚力係数を求めることが可能であることが確認された。
【００５４】
図６に各種の翼型についての、最大揚力係数の実験値と計算値の比較を示す。斜めに引かれた破線上に点が乗っていれば、実験値と計算値が完全に一致することを示している。ここに示されるように、従来法では実験値に比べて大きく外れるとともに、破線からの散らばりが大きいのに対し、本発明でα＝３．０と選んだ場合には、実験値に近づくとともに、散らばりも小さくなることが確認された。
【００５５】
【発明の効果】
請求項１記載の本発明によれば、翼面近傍の渦度に対して、壁面に沿う方向Ｘの速度勾配絶対値Ｓを、Ｓ＝Ｍａｘ（Ω−α｜Ｕｘ｜，０）またはＳ＝Ｍａｘ（Ω−αＭａｘ（０，Ｕｘ），０）によって修正するので、高迎角においても速度勾配絶対値Ｓを定常解に収束させることができ、最大揚力および失速角などの翼特性を正確に予測することが可能となり、翼の最大揚力付近の特性を、たとえばコンピュータを用いた数値シミュレーションによって高精度で模擬することができる。
【００５６】
請求項２記載の本発明によれば、翼近傍では流体の性質から、壁面に対して垂直方向速度の垂直方向勾配をＶｙとしたとき、
Ｕｘ＋Ｖｙ＝０
が近似的に成立することを用いて、Ｕｘを−Ｖｙで置換えるので、２次元断面では容易であるが、３次元の翼で流れが翼を斜めに横切る場合には、困難な壁面方向速度の壁面に沿う方向の速度勾配Ｕｘの計算を、この場合でも算出が容易な壁面に対して垂直方向速度の垂直方向勾配Ｖｙで置換えることができる。
【００５７】
請求項３記載の本発明によれば、定数αが複数の異なる翼形の最大揚力係数が計算値と実験値の相関が最良になる値に選ばれるので、高い信頼性で高迎角における速度勾配絶対値Ｓが求められ、最大揚力および失速角などの翼特性の予測精度が向上され、翼設計の高精度化が図られるとともに、翼断面決定に要する時間の短縮を図ることができる。
【００５８】
請求項４記載の本発明によれば、レイノルズ数Ｒｅに応じて定数αの値を変化させることによって、流れの慣性力と粘性力との相対的影響度を反映させて、前記速度勾配の絶対値Ｓが求められ、最大揚力および失速角などの翼特性を正確に予測することができる。
【００５９】
請求項５記載の本発明によれば、複数のレイノルズ数から求めた値を境界層近似した値に置換えた定数αを用いて、前記速度勾配の絶対値Ｓを求めるので、高迎角の最大揚力をより高い信頼性で求めることができる。
【００６０】
請求項６記載の本発明によれば、コンピュータによって上記翼形性能の推定プログラムが実行されることによって、翼面近傍の渦度に対して、壁面に沿う方向Ｘの速度勾配｜Ω｜を、｜Ω｜＝Ｍａｘ（｜Ω｜−α｜Ｕｘ｜，０）または｜Ω｜＝Ｍａｘ（｜Ω｜−αＭａｘ（０，Ｕｘ），０）によって修正するので、高迎角においても速度勾配絶対値Ｓを定常解に収束させることができ、最大揚力および失速角などの翼特性を正確に予測することが可能となり、翼の最大揚力付近の特性を高精度で模擬することが可能な翼形性能の解析ツールを、コンピュータ上に実現することができる。
【００６１】
請求項７記載の本発明によれば、翼近傍では流体の性質から、壁面に対して垂直方向速度の垂直方向勾配をＶｙとしたとき、
Ｕｘ＋Ｖｙ＝０
が近似的に成立することを用いて、Ｕｘを−Ｖｙで置換えるので、２次元断面では容易であるが、３次元の翼で流れが翼を斜めに横切る場合には、困難な壁面方向速度の壁面に沿う方向の速度勾配Ｕｘの計算を、この場合でも算出が容易な壁面に対して垂直方向速度の垂直方向勾配Ｖｙで置換えることができる。
【００６２】
請求項８記載の本発明によれば、定数αが複数の異なる翼形の最大揚力係数の計算値と実験値の相関が最良になる値に選ばれるので、高い信頼性で高迎角における速度勾配絶対値Ｓが求められ、最大揚力および失速角などの翼特性の予測精度が向上され、翼設計の高精度化が図られるとともに、翼断面決定に要する時間の短縮を図ることができる。
【００６３】
請求項９記載の本発明によれば、レイノルズ数Ｒｅに応じて定数αの値を変化させることによって、流れの慣性力と粘性力との相対的影響度を反映させて、前記速度勾配の絶対値Ｓが求められ、最大揚力および失速角などの翼特性を正確に予測することができる。
【００６４】
請求項１０記載の本発明によれば、複数のレイノルズ数から求めた値を簡易式で近似した値に置換えた定数αを用いて、前記速度勾配の絶対値Ｓを求めるので、高迎角の最大揚力をより高い信頼性で求めることができる。簡易式の例としては、２つのレイノルズ数Ｒｅ１，Ｒｅ２において得られた定数をα１，α２とし、任意のレイノルズ数Ｒｎでの定数αを、次式
α＝α１＋（Ｒｎ−Ｒｅ１）（α２−α１）／（Ｒｅ２−Ｒｅ１）
で与える方法がある。
【図面の簡単な説明】
【図１】本発明の実施の一形態の翼形性能の推定方法を示すフローチャートである。
【図２】翼形性能の推定プログラムを実行して翼形性能をシミュレーションするために用いられるコンピュータ１の電気的構成を示すブロック図である。
【図３】１２％翼形（ＮＡＣＡ６３１−０１２）の揚力ＣＬ、抵抗ＣＤ、モーメントＣＭの３分力特性の本件発明者による計算結果と文献による実験結果とを示すグラフである。
【図４】１８％翼形（ＮＡＣＡ６３３−０１８）の揚力ＣＬ、抵抗ＣＤ、モーメントＣＭの３分力特性の本件発明者による計算結果と文献による実験結果とを示すグラフである。
【図５】従来の技術の翼の性能を推定するための概略的なアルゴリズムを示すフローチャートである。
【図６】実験と計算による複数の翼型の最大揚力係数の比較について、従来と本発明による方法とを示すグラフである。
【符号の説明】
１  コンピュータ
２  入力手段
３  記憶手段
４  演算処理装置
５  表示手段
【図１】
0003587827.tif 000003
【図２】
0003587827.tif 000004
【図３】
0003587827.tif 000005
【図４】
0003587827.tif 000006
【図５】
0003587827.tif 000007
【図６】
0003587827.tif 000008