JP3587827B2 - 翼形性能の推定方法および翼形性能の推定プログラム - Google Patents

翼形性能の推定方法および翼形性能の推定プログラム Download PDF

Info

Publication number
JP3587827B2
JP3587827B2 JP2002164867A JP2002164867A JP3587827B2 JP 3587827 B2 JP3587827 B2 JP 3587827B2 JP 2002164867 A JP2002164867 A JP 2002164867A JP 2002164867 A JP2002164867 A JP 2002164867A JP 3587827 B2 JP3587827 B2 JP 3587827B2
Authority
JP
Japan
Prior art keywords
value
wing
velocity
velocity gradient
constant
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP2002164867A
Other languages
English (en)
Other versions
JP2004012248A (ja
Inventor
英志 嶋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Kawasaki Motors Ltd
Original Assignee
Kawasaki Jukogyo KK
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Kawasaki Jukogyo KK filed Critical Kawasaki Jukogyo KK
Priority to JP2002164867A priority Critical patent/JP3587827B2/ja
Publication of JP2004012248A publication Critical patent/JP2004012248A/ja
Application granted granted Critical
Publication of JP3587827B2 publication Critical patent/JP3587827B2/ja
Anticipated expiration legal-status Critical
Expired - Fee Related legal-status Critical Current

Links

Images

Landscapes

  • Complex Calculations (AREA)

Description

【0001】
【発明の属する技術分野】
本発明は、回転翼航空機のロータブレードの翼形の決定および翼形決定に関与する最大揚力係数などの翼特性の確認などを行うための翼形性能を、高精度で推定することができる翼形性能の推定方法およびコンピュータによって実行可能な翼形性能の推定プログラムに関する。
【0002】
本発明において、「翼型」とは、米国航空評議委員会(National Advisory
Committee for Aeronautics;略称NACA)で規定・分類された2次元の翼断面形状をいうものとし、この「翼型」を含む任意の翼断面形状をいう場合には、「翼形」と記す(以下同じ)。
【0003】
【従来の技術】
航空機の翼(3次元形状をいう)および翼形(2次元断面形状をいう)は、航空機の性能に決定的な影響をもっており、翼形を含む翼については次のようなことが知られている。すなわち、翼は気流に対する角度である迎角が増大するにつれて、揚力が徐々に大きくなるが、ある角度で最大揚力を発生し、その後、揚力の低下および抗力の急増を引き起こす失速現象を生じる。この最大揚力およびその発生角度は、翼の重要な特性のひとつである。
【0004】
このような翼の特性は、風洞実験によって測定されることが多いが、この風洞実験では実の飛行状態を模擬することが難しく、多くの時間およびコストを要するという問題があり、近年のコンピュータの急速な進歩に応じて、数値シミュレーションによる流れ解析によって、翼特性の測定などが多く行われている。この流れの数値シミュレーションは、数値流体力学(略称CFD;Computational Fluid Dynamics)が用いられ、航空機の空力設計における開発経費の削減および期間の短縮に大きく貢献でき、ナビエ・ストークス解析によって、0°〜13°程度の低迎角では遷音速特性、および乱流遷移の問題を除く抗力の推算も、高い精度で可能であるが、レイノルズ数Reが10にも達する流れ場における最大揚力性能の推定を行うには、現在のコンピュータの処理能力では、乱流現象を近似した有限体積法3次元レイノルズ平均Navier−Stokes(Reynolds Averaged Navier
−Stokes;略称RANS)解析を用いるしかなく、乱流モデルの限界によって、最大揚力の推定精度は制限されている。
【0005】
図5は従来の技術の翼の性能を推定するための概略的なアルゴリズムを示すフローチャートである。たとえばヘリコプタのロータブレードの翼断面の性能を、コンピュータを用いてシミュレーションするにあたっては、ステップa1で、コンピュータによって翼特性解析プログラムが実行され、ステップa2で、翼形および流れなどに関する諸条件を入力して初期設定すると、ステップa2で、入力された諸条件に基づいて乱流粘性係数の計算が開始される。
【0006】
この乱流粘性係数は、これを数値的に解くために、連続的な偏微分方程式を離散的な格子点上で定義される差分方程式として書き直す必用がある。流れ場の計算格子は、流速uの方向をX方向とし、動粘性係数νの方向をY方向とするXY直交座標系上で、流れ場をX方向およびY方向の長さがΔx,Δyの格子で分割し、格子点番号(i,j)(ここに、i,jは正の整数)で表される各格子点で囲まれた格子中央に圧力pを定義する。時間ステップをΔtとしたとき、時刻t=t+nΔt(ここに、nは正の整数)における各格子点位置の渦度ω(x,y;t)と流れ関数ψ(x,y;t)とを導入することによって、流れを規定する2次元ナビエ・ストークス方程式と連続の式とを、渦度輸送方程式とポアソン方程式とに変形し、それらを差分化する。
【0007】
次に、ステップa4で、流れ場の計算が開始され、前記差分化によって得られた差分方程式に対して、時刻tに関する計算を全格子点について行い、ステップa5で流れ関数ψが収束するまで上記ステップa3,a4を繰り返し、ステップa5で定常解に収束すると、次のステップa6に移り、迎角に対する揚力係数などの翼形状特性が求められ、ステップa7で計算が終了する。
【0008】
前記ステップa3において、渦度|Ω|は、|Uy−Vx|によって計算される。ここで、UyはY方向の速度勾配を示す偏微分であり、VxはX方向の速度勾配を示す偏微分であり、U,VはX,Y方向の速度である。Uy,Vxは差分近似によって計算される。このような各速度勾配Uy,Vxの近似法は、従来から多くの手法が周知である。流れの計算方法には、上記のように渦度ωと流れ関数ψとを用いるω‐ψ法以外に、流速u,vと圧力pとを未知変数として用いるu−p法がある。さらに翼面に平行な速度勾配を次式、
{2Ux+Uy+2Uy・Vx+Vx−2/3・(Ux+Vy)1/2
によって求めることもできる。
【0009】
このような各種の手法のいずれによって渦度Ω等を求めても、翼面近傍の境界層では、翼面に平行な速度成分だけが残り、翼面に垂直な方向の変化量が非常に少なくなるので、結局は同じ値に収束して、壁面がX軸に平行な場合、Y方向の速度勾配と同等になり、したがって渦度Ωも同等になることが知られているが、約13°〜15°の高迎角における最大揚力の定常的予測は困難であり、翼の空力設計を高精度でシミュレートすることができないという問題がある。
【0010】
【発明が解決しようとする課題】
本発明の目的は、翼の最大揚力付近の特性を高精度で予測することができる翼形性能の推定方法および翼形性能の推定プログラムを提供することである。
【0011】
【課題を解決するための手段】
請求項1記載の本発明は、翼面近傍の渦度に漸近する流れ方向の速度勾配絶対値Sを生成項とする乱流モデルにおいて、2つの変数a,bの大なる方を値とする関数をMax(a,b)とし、翼面に沿う速度の翼面に沿う方向の速度勾配をUxとし、壁面に垂直な方向の速度勾配をUyとし、Uyに漸近する渦度値等をΩとし、定数をαとし、座標軸Xは翼に対する気流の相対的な速度Uの方向を正としたとき、
S=Max(|Ω|−α|Ux|,0)
または、
S=Max(|Ω|−αMax(0,Ux),0)
によって修正された速度勾配絶対値Sを用いて、翼の迎角に対する揚力を求めることを特徴とする翼形性能の推定方法である。
【0012】
本発明に従えば、翼面近傍の渦度に対して、壁面に沿う方向Xの速度勾配Sを、S=|Ω|−α|Ux|またはS=|Ω|−αMax(0,Ux)によって修正するので、高迎角においても速度勾配絶対値Sを定常解に収束させることができ、最大揚力および失速角などの翼特性を正確に予測することが可能となり、翼の最大揚力付近の特性を、たとえばコンピュータを用いた数値シミュレーションによって高精度で模擬することが可能となる。
【0013】
請求項2記載の本発明は、翼近傍では流体の性質から、壁面に対して垂直方向速度の垂直方向勾配をVyとしたとき、
Ux+Vy=0
が近似的に成立することを用いて、前記翼面に沿う方向の速度勾配Uxを翼面に垂直な方向の速度勾配−Vyで置換えることを特徴とする。
【0014】
本発明に従えば、2次元断面では壁面に沿う方向の速度の壁面に沿う方向の速度勾配Uxの計算は容易であるが、3次元の翼で流れが翼を斜めに横切る場合には、この算出は容易ではない。一方、壁面に対して垂直方向速度の垂直方向勾配Vyは、この場合でも算出が容易であるという利点がある。
【0015】
請求項3記載の本発明は、前記定数αは、複数の異なる既知の翼形の最大揚力の計算値と実験値との相関が最良になる値に選ばれることを特徴とする。
【0016】
本発明に従えば、定数αが複数の異なる翼形の最大揚力係数が実験値と最も近くなる値に選ばれるので、高い信頼性で高迎角における速度勾配絶対値Sが求められ、最大揚力および失速角などの翼特性の予測精度が向上され、翼設計の高精度化が図られるとともに、翼断面決定に要する時間の短縮を図ることができる。
【0017】
請求項4記載の本発明は、前記定数αは、気流の流れのレイノルズ数に応じて異なる値に選ばれることを特徴とする。
【0018】
本発明に従えば、流れの状態は、流れの代表速度をU、流れの代表長さをL、動粘度をνとしたとき、レイノルズ数Re=UL/ν、すなわちRe=(慣性力)/(粘性力)で表され、このレイノルズ数Reに応じて定数αの値を変化させることによって、流れの慣性力と粘性力との相対的影響度を反映させて、前記速度勾配の絶対値Sが求められ、最大揚力および失速角などの翼特性を正確に予測することが可能となる。
【0019】
請求項5記載の本発明は、前記定数αは、複数のレイノルズ数に基づいて求めた値を、簡易式で近似した値に置換えることを特徴とする。
【0020】
本発明に従えば、複数のレイノルズ数から求めた値を簡易式で近似した値に置換えた定数αを用いて、前記速度勾配の絶対値Sを求めるので、高迎角の最大揚力をより高い信頼性で求めることができる。簡易式の例としては、2つのレイノルズ数Re1,Re2において得られた定数をα1,α2とし、任意のレイノルズ数Rnでの定数αを、次式
α=α1+(Rn−Re1)(α2−α1)/(Re2−Re1)
で与える方法がある。
【0021】
請求項6記載の本発明は、翼面近傍の渦度に漸近する速度勾配絶対値Sを生成項とする乱流モデルを用いて翼形性能を求めるコンピュータを、
翼形性能を生成するために必要なデータを入力するための入力手段、
入力手段によって入力されたデータに基づいて、2つの変数a,bの大なる方を値とする関数をMax(a,b)とし、翼面に沿う速度の翼面に沿う方向の速度勾配をUxとし、壁面に垂直な方向の速度勾配をUyとし、Uyに漸近する渦度値等をΩとし、定数をαとし、座標軸Xは翼に対する気流の相対的な速度Uの方向を正としたとき、
S=Max(|Ω|−α|Ux|,0)
または、
S=Max(|Ω|−αMax(0,Ux),0)
によって修正された速度勾配絶対値Sを用いて、翼の迎角に対する揚力を求める演算手段、および
演算手段によって求められた前記翼の迎角に対する揚力を表示する出力手段、として機能させるコンピュータによって実行可能な翼形性能の推定プログラムである。
【0022】
本発明に従えば、コンピュータによって上記翼形性能の推定プログラムが実行され、入力手段によって翼形性能を生成するために必用な、たとえば翼断面形状、流れの条件を特定するための諸量、格子点などの各種のデータが入力されると、演算手段は前記入力手段から入力されたデータに基づいて、翼面近傍の渦度に対して、壁面に沿う方向Xの速度勾配絶対値Sを、S=|Ω|−α|Ux|またはS=|Ω|−αMax(0,Ux)によって修正し、これによって高迎角においても速度勾配絶対値Sを定常解に収束させて、出力手段に翼の迎角に対する揚力を表示させることができる。このようにして出力手段に翼の迎角と揚力とが表示されるので、最大揚力および失速角などの翼特性を正確に予測することが可能となり、翼の最大揚力付近の特性を高精度で模擬することが可能な翼形性能の解析ツールをコンピュータ上に実現することができる。
【0023】
請求項7記載の本発明は、翼近傍では流体の性質から、壁面に対して垂直方向速度の垂直方向勾配をVyとしたとき、
Ux+Vy=0
が近似的に成立することを用いて、前記翼面に沿う方向の速度勾配Uxを翼面に垂直な方向の速度勾配−Vyで置換えることを特徴とする。
【0024】
本発明に従えば、2次元断面では壁面に沿う方向の速度の壁面に沿う方向の速度勾配Uxの計算は容易であるが、3次元の翼で流れが翼を斜めに横切る場合には、この算出は容易ではない。一方、壁面に対して垂直方向速度の垂直方向勾配Vyは、この場合でも算出が容易であるという利点がある。
【0025】
請求項8記載の本発明は、前記定数αは、複数の異なる既知の翼形の最大揚力の計算値と実験値の相関が最良になる値に選ばれることを特徴とする。
【0026】
本発明に従えば、定数αが複数の異なる翼形の最大揚力係数が実験値と最も近くなる値に選ばれるので、高い信頼性で高迎角における速度勾配絶対値Sが求められ、最大揚力および失速角などの翼特性の予測精度が向上され、翼設計の高精度化が図られるとともに、翼断面決定に要する時間の短縮を図ることができる。
【0027】
請求項9記載の本発明は、前記定数αは、気流の流れのレイノルズ数に応じて異なる値に選ばれることを特徴とする。
【0028】
本発明に従えば、流れの状態は、流れの代表速度をU、流れの代表長さをL、動粘度をνとしたとき、レイノルズ数Re=UL/ν、すなわちRe=(慣性力)/(粘性力)で表され、このレイノルズ数Reに応じて定数αの値を変化させることによって、流れの慣性力と粘性力との相対的影響度を反映させて、前記速度勾配の絶対値Sが求められ、最大揚力および失速角などの翼特性を正確に予測することが可能となる。
【0029】
請求項10記載の本発明は、前記定数αは、複数のレイノルズ数に基づいて求めた値を、簡易式で近似した値に置換えることを特徴とする。
【0030】
本発明に従えば、複数のレイノルズ数から求めた値を境界層近似した値に置換えた定数αを用いて、前記速度勾配の絶対値Sを求めるので、高迎角の最大揚力をより高い信頼性で求めることができる。
【0031】
【発明の実施の形態】
図1は本発明の実施の一形態の翼形性能の推定方法を示すフローチャートであり、図2は翼形性能の推定プログラムを実行して翼形性能をシミュレーションするために用いられるコンピュータ1の電気的構成を示すブロック図である。本実施の形態では、ヘリコプタのロータブレードの翼断面の性能を、コンピュータを用いて推定する方法について説明する。ステップb1で、コンピュータ1によって翼形性能の推定プログラムが実行され、翼形および流れなどに関する諸条件が入力手段2によって入力されて、初期設定される。前記翼形性能の推定プログラムは、記憶手段3に記憶され、演算処理装置4からの実行指令によって読み出されて実行され、出力手段である表示手段5に計算結果を表示させ、数値流体解析システムを構築する。
【0032】
このコンピュータ1による解析は、圧縮性有限体積法3次元レイノルズ平均ナビエ・ストークス(Reynolds Averaged Navier−Stokes,略称RANS)方程式を使用する。レイノルズ応力は乱流粘性として導入し、乱流モデルは後述のSpalant−Allmaras(略称SA)乱流モデルを用いる。計算手法は、有限体積法を用い、基礎変数にVan Albadaのリミタを適用したMUSCL(Monotonic Upwind Scheme for
Conservation Low)によって空間2次精度とする。近似Riemann流速としては、SHUS(Simple High resolution Upwind Scheme)を使用する。また、時間積分には、MFGS(Matrix Free Gauss Seidel)陰解法を用い、収束加速のために局所時間法を併用する。
【0033】
この局所時間法によって収束解が得られない場合には、物理時間法による非定常計算に切り換える。この場合の時間精度は一次とし、時間刻みはコード長と音速とを基準にした無次元時間で0.01とする。計算は、最小迎角から開始し、収束解を次の迎角の初期値として使用する。高迎角においては、実験においても流れの初期値依存性を生じることがあるので、風洞試験での迎角変化の模擬を念頭においている。
【0034】
計算格子は、代数的手法によるO型格子で、全4角形格子を使用し、格子点数は周方向に256とし、垂直方向は60とする。また翼弦長を基準として、垂直方向の最小格子間隔は10−5とし、計算領域の半径は20とし、周方向の最小格子間隔は前縁および後縁で約8×10−4とし、6%翼の前縁については1.6×10−4とする。流れの条件は、マッハ数が0.2、レイノルズ数が5.8×10の一様流とする。
【0035】
こうしてステップb1で初期設定すると、ステップb2で、入力された諸条件に基づいて乱流粘性係数の計算が開始される。この乱流粘性係数を求めるにあたっては、これを数値的に解くために、連続的な偏微分方程式を離散的な格子点上で定義される差分方程式として書き直す必用がある。すなわち、ステップb21で、流れ場の計算格子は、流速uの方向をX方向とし、動粘性係数νの方向をY方向とするXY座標上で、流れ場をX方向およびY方向の長さがΔx,Δyの格子で分割し、格子点番号(i,j)(ここに、i,jは正の整数)で表される各格子点で囲まれた格子中央に圧力pを定義する。時間ステップをΔtとしたとき、時刻t=t+nΔt(ここに、nは正の整数)における各格子点位置の渦度ω(x,y;t)と流れ関数ψ(x,y;t)とを導入することによって、流れを規定する2次元ナビエ・ストークス方程式と連続の式とを、渦度輸送方程式とポアソン方程式とに変形し、それらを差分化して、X方向の速度Vxおよびその速度勾配UxならびにY方向の速度Vyおよびその速度勾配Uyを求める。ここで、UyはY方向の速度勾配を示す偏微分であり、VxはX方向の速度勾配を示す偏微分である。
【0036】
上記のナビエ・ストークス方程式をコンピュータ1によって数値シミュレーションするには、何らかの形で近似化する必要があり、このような近似化の手法として、有限差分法、有限体積法、および有限要素法などを用いて、現実の連続な空間を、数多くの離散的な計算点で置換え、ナビエ・ストークス方程式を差分方程式などによって近似し、計算を実行する。最大揚力の予測に関しては、乱流モデルを用いる。
【0037】
この乱流モデルには、多くの種類があるが、共通点として、何らかの形式で速度の空間勾配を用いていることである。x,y,zを3次元の空間座標、u,v,wを各方向X,Y,Zの速度とすると、速度勾配は速度の各方向の座標での偏微分の記号によって、u,u,u,ν,ν,ν,w,w,wで表される。
【0038】
これらの式で表される偏微分は、ステップb21において、数値計算上、差分で近似される。たとえば3次元の直交格子で空間を離散化し、その格子点の3次元座標が格子インデックス(i,j,k)とX,Y,Zの各方向の格子間隔Δx,Δy,Δzとを用いて、次のように表される場合を考える。
(xi,j,k,yi,j,k,zi,j,k)=(iΔx,jΔy,kΔz)
【0039】
この格子上の各点は、インデックス(i,j,k)で区別され、(i,j,k)での速度はui,j,kのように表される。この場合、インデックス(i,j,k)における速度勾配uxi,j,kは、次のように近似することができる。
x,i,j,k=(ui+1,j,k−ui−1,j,k)/2Δx
【0040】
数値流体力学(略称CFD;Computational Fluid Dynamics)における格子およびメッシュの選び方、その上での差分近似の方法は多くの手法が周知であるが、どのような格子、メッシュを用い、どのような差分近似を用いても、格子を細かくしていけば、元の偏微分に限りなく近づく。したがって、乱流モデルの適用において、具体的な差分近似の手法は、その効果に対しては副次的な影響しかもたらさない。そのため、現実の適用にあたっては、それらを差分近似などで離散化したものを用いるが、本実施の形態では、差分近似表示ではなく、元の偏微分の記号を使って説明する。乱流モデルの適用は、このような計算を数万〜数百万点存在するすべての格子点に対して、数百〜数万ステップある全時間ステップについて、繰り返すことによって実現される。
【0041】
乱流モデル中での速度勾配指標は、乱流モデルの効果に直接の影響があり、いくつかの計算手法が知られており、たとえば渦度ωと各速度勾配u,u,ν,ν,wとの関係を次式のように想定する次のような手法がある。
|Ω|={(u−ν+(ν−w+(w−u1/2
【0042】
航空機まわりの流れのような高レイノルズ数の流れでは、物体近傍に境界層と呼ばれるごく薄い、速度の遅い領域が発生する。物体表面では速度が0になるので、境界層の中で、壁に近づくにつれて速度は急速に減少する。このため、速度勾配の成分のうち壁に垂直な方向の変化が平行な方向に比べてはるかに大きくなることが知られている。また、壁の近傍では、流れは明らかに壁に平行な方向であるので、壁の近くでは壁に沿う速度成分の壁に垂直な方向の変化は、他の速度勾配に比べて大きくなることが知られている。仮に、壁がXY平面に平行な場合、UxとWy成分の大きさ、すなわち絶対値は、他のVx,Vy,Vz,Ux,Uy,Uz,Wx,Wzよりも大きくなることが知られている。
【0043】
上記のように、高レイノルズ数の流れでは、大きな速度勾配が境界層の内部に発生し、壁に沿う速度成分の壁に垂直な方向の変化が、他の成分に比べて大きくなるが、翼の最大揚力付近では、翼の前縁付近で壁に沿った方向の速度勾配がかなり大きな値となる。
【0044】
乱流モデルについて説明する。最大揚力係数の予測精度に関して、本件発明者は数種類の乱流モデルをテストした結果、最良であったSpalart−Allmarasモデル(SAモデルと略記する)を採用する。テストした他の乱流モデルは、Baldwin−Lomax代数モデル(略称BL)、Baldwin−Barthモデル、q−ωモデルである。SAモデルは、基本的には乱流粘性係数に対するソース項付移流拡散方程式を解く一方程式モデルであり、次のように表される。なお、下式中の記載「ν」は外部領域で乱流粘性係数に漸近する乱流変数である。
【0045】
Figure 0003587827
ここで、νは分子粘性係数、νは乱流粘性係数、Sは渦度Ωの絶対値(=|Ω|)、各定数の値は、cb1=0.1355,σ=2/3,cb2=0.622,k=0.41,cw1=cb1/k2+(1+cb2)/σ,cw2=0.3,cw3=2,cv1=7.1,ct1=1,ct2=2,ct3=1.2,ct4=0.5である。
【0046】
上記のSA乱流モデルについて、次のような変更を行っている。
変数変換に関しては、乱流方程式について非構造有限体積法で二次精度法を適用するには保存形式の方が扱いやすいため、変数をν ̄からμ ̄≡ρν ̄に変換した等価な次式を用いる。
(μ ̄)+(uμ ̄)+(νμ ̄)=ρ[R.H.S.]
ただし、ρは密度、uはX方向の速度、νはY方向の速度であり、これらのρ,u,νには質量保存側、
ρ+(ρu)+(ρν)=0
が成立しているものとする。
【0047】
空間二次精度化に関しては、乱流方程式には一次精度風上を用いる場合が多く、平板境界層の解析等では精度的に充分であるが、高迎角では格子解像度の影響が顕著になるので、本実施の形態ではナビエ・ストークス方程式に対するのと同様のMUSCLを用いて二次精度化する。上記の高迎角において格子解像度の影響が顕著になる点については、本件発明者の研究によって既に確認されており、航空宇宙学会30周年記念講演会講演集(1999年)、第66頁〜第68頁に「CFDによる最大揚力係数の予測について」として掲載されている。
【0048】
全場乱流に関しては、SAモデルには遷移点を制御するための遷移項(ft1)と、それ以外の遷移を回避するための安定化項(ft2)とが含まれているが、本実施の形態ではこれらの二項を零とし、全場乱流の条件で解析する。ただし、一様流での値は低く保たれており、淀み点から境界層内で急速に増大するので、結果的には数値的な遷移を用いていることになる。
【0049】
次に、ステップb22で、渦度の絶対値|Ω|、
|Ω|=SQRT|Uy−Vx|
を計算し、ステップb23で、翼面である壁に平行方向の速度勾配Vt、
Vt=|UxXt+(Vx+Uy)XtYt+VyYt
を計算する。
【0050】
ステップb24で、渦度ωの修正計算を次式、
|Ω|←Max(|Ω|Uy−αVt,0)
または、
|Ω|←Max(|Ω|−αMax(0,Ux),0)
によって行う。上式において、Max(0,Ux)は、2つの変数0,Uxの大なる方を変数値として選択する関数を意味する。翼面に沿う速度の翼面に沿う方向の速度勾配をUxとし、定数をαとする。この定数αは、0.01<α<100に選ばれ、代表的には3.0に選ばれる。座標軸Xは、翼に対する気流の相対的な速度Uの方向を正とする。
【0051】
ステップb25で、上記の式3または式4によって修正された渦度の絶対値|Ω|を用いて、翼の迎角に対する揚力を求め、次に、ステップb25で、前記差分化によって得られた差分方程式に対して、時刻tに関する計算を行い、ステップa4で流れ関数ψが収束するまで上記ステップa2,a3を繰り返し、ステップa5で、迎角に対する乱流粘性係数などの翼特性が求められる。
【0052】
図4は12%翼形(NACA63−012)の揚力CL、抵抗CD、モーメントCMの3分力特性の本件発明者による実験結果を示すグラフである。同図において、横軸は迎角であり、縦軸は揚力CL、抵抗CD、モーメントCMであり、本発明の実施例のシミュレーション結果は「SA改」と示され、従来手法による比較例のシミュレーション結果は「従来」と示され、従来手法による実験値は「Exp」と示される。迎角13°以上では計算結果と実験結果とはやや異なるが、最大揚力は計算によってほぼ正しく予想されることが確認された。
【0053】
図5は18%翼形(NACA63−018)の揚力CL、抵抗CD、モーメントCMの3分力特性の本件発明者による実験結果を示すグラフである。同図において、横軸は迎角であり、縦軸は揚力CL、抵抗CD、モーメントCMであり、本発明の実施例のシミュレーション結果は「SA改」と示され、従来手法による比較例のシミュレーション結果は「従来」と示され、従来手法による実験値は「Exp」と示される。図4に示される翼(NACA63−012)と同様に、最大揚力係数の計算値が実験値に近づいており、高迎角においても高精度で最大揚力係数を求めることが可能であることが確認された。
【0054】
図6に各種の翼型についての、最大揚力係数の実験値と計算値の比較を示す。斜めに引かれた破線上に点が乗っていれば、実験値と計算値が完全に一致することを示している。ここに示されるように、従来法では実験値に比べて大きく外れるとともに、破線からの散らばりが大きいのに対し、本発明でα=3.0と選んだ場合には、実験値に近づくとともに、散らばりも小さくなることが確認された。
【0055】
【発明の効果】
請求項1記載の本発明によれば、翼面近傍の渦度に対して、壁面に沿う方向Xの速度勾配絶対値Sを、S=Max(Ω−α|Ux|,0)またはS=Max(Ω−αMax(0,Ux),0)によって修正するので、高迎角においても速度勾配絶対値Sを定常解に収束させることができ、最大揚力および失速角などの翼特性を正確に予測することが可能となり、翼の最大揚力付近の特性を、たとえばコンピュータを用いた数値シミュレーションによって高精度で模擬することができる。
【0056】
請求項2記載の本発明によれば、翼近傍では流体の性質から、壁面に対して垂直方向速度の垂直方向勾配をVyとしたとき、
Ux+Vy=0
が近似的に成立することを用いて、Uxを−Vyで置換えるので、2次元断面では容易であるが、3次元の翼で流れが翼を斜めに横切る場合には、困難な壁面方向速度の壁面に沿う方向の速度勾配Uxの計算を、この場合でも算出が容易な壁面に対して垂直方向速度の垂直方向勾配Vyで置換えることができる。
【0057】
請求項3記載の本発明によれば、定数αが複数の異なる翼形の最大揚力係数が計算値と実験値の相関が最良になる値に選ばれるので、高い信頼性で高迎角における速度勾配絶対値Sが求められ、最大揚力および失速角などの翼特性の予測精度が向上され、翼設計の高精度化が図られるとともに、翼断面決定に要する時間の短縮を図ることができる。
【0058】
請求項4記載の本発明によれば、レイノルズ数Reに応じて定数αの値を変化させることによって、流れの慣性力と粘性力との相対的影響度を反映させて、前記速度勾配の絶対値Sが求められ、最大揚力および失速角などの翼特性を正確に予測することができる。
【0059】
請求項5記載の本発明によれば、複数のレイノルズ数から求めた値を境界層近似した値に置換えた定数αを用いて、前記速度勾配の絶対値Sを求めるので、高迎角の最大揚力をより高い信頼性で求めることができる。
【0060】
請求項6記載の本発明によれば、コンピュータによって上記翼形性能の推定プログラムが実行されることによって、翼面近傍の渦度に対して、壁面に沿う方向Xの速度勾配|Ω|を、|Ω|=Max(|Ω|−α|Ux|,0)または|Ω|=Max(|Ω|−αMax(0,Ux),0)によって修正するので、高迎角においても速度勾配絶対値Sを定常解に収束させることができ、最大揚力および失速角などの翼特性を正確に予測することが可能となり、翼の最大揚力付近の特性を高精度で模擬することが可能な翼形性能の解析ツールを、コンピュータ上に実現することができる。
【0061】
請求項7記載の本発明によれば、翼近傍では流体の性質から、壁面に対して垂直方向速度の垂直方向勾配をVyとしたとき、
Ux+Vy=0
が近似的に成立することを用いて、Uxを−Vyで置換えるので、2次元断面では容易であるが、3次元の翼で流れが翼を斜めに横切る場合には、困難な壁面方向速度の壁面に沿う方向の速度勾配Uxの計算を、この場合でも算出が容易な壁面に対して垂直方向速度の垂直方向勾配Vyで置換えることができる。
【0062】
請求項8記載の本発明によれば、定数αが複数の異なる翼形の最大揚力係数の計算値と実験値の相関が最良になる値に選ばれるので、高い信頼性で高迎角における速度勾配絶対値Sが求められ、最大揚力および失速角などの翼特性の予測精度が向上され、翼設計の高精度化が図られるとともに、翼断面決定に要する時間の短縮を図ることができる。
【0063】
請求項9記載の本発明によれば、レイノルズ数Reに応じて定数αの値を変化させることによって、流れの慣性力と粘性力との相対的影響度を反映させて、前記速度勾配の絶対値Sが求められ、最大揚力および失速角などの翼特性を正確に予測することができる。
【0064】
請求項10記載の本発明によれば、複数のレイノルズ数から求めた値を簡易式で近似した値に置換えた定数αを用いて、前記速度勾配の絶対値Sを求めるので、高迎角の最大揚力をより高い信頼性で求めることができる。簡易式の例としては、2つのレイノルズ数Re1,Re2において得られた定数をα1,α2とし、任意のレイノルズ数Rnでの定数αを、次式
α=α1+(Rn−Re1)(α2−α1)/(Re2−Re1)
で与える方法がある。
【図面の簡単な説明】
【図1】本発明の実施の一形態の翼形性能の推定方法を示すフローチャートである。
【図2】翼形性能の推定プログラムを実行して翼形性能をシミュレーションするために用いられるコンピュータ1の電気的構成を示すブロック図である。
【図3】12%翼形(NACA63−012)の揚力CL、抵抗CD、モーメントCMの3分力特性の本件発明者による計算結果と文献による実験結果とを示すグラフである。
【図4】18%翼形(NACA63−018)の揚力CL、抵抗CD、モーメントCMの3分力特性の本件発明者による計算結果と文献による実験結果とを示すグラフである。
【図5】従来の技術の翼の性能を推定するための概略的なアルゴリズムを示すフローチャートである。
【図6】実験と計算による複数の翼型の最大揚力係数の比較について、従来と本発明による方法とを示すグラフである。
【符号の説明】
1 コンピュータ
2 入力手段
3 記憶手段
4 演算処理装置
5 表示手段

Claims (10)

  1. 翼面近傍の渦度に漸近する流れ方向の速度勾配絶対値Sを生成項とする乱流モデルにおいて、2つの変数a,bの大なる方を値とする関数をMax(a,b)とし、翼面に沿う速度の翼面に沿う方向の速度勾配をUxとし、壁面に垂直な方向の速度勾配をUyとし、Uyに漸近する渦度値等をΩとし、定数をαとし、座標軸Xは翼に対する気流の相対的な速度Uの方向を正としたとき、
    S=Max(|Ω|−α|Ux|,0)
    または、
    S=Max(|Ω|−αMax(0,Ux),0)
    によって修正された速度勾配絶対値Sを用いて、翼の迎角に対する揚力を求めることを特徴とする翼形性能の推定方法。
  2. 壁面に垂直な方向の速度をVとしたとき、その垂直方向の勾配はVyで表わし、流体の性質から、壁面近傍では、
    Ux+Vy=0
    が概ね成立することを利用し、前記翼面に沿う方向の速度勾配Uxを翼面に垂直な方向の速度勾配−Vyで置換えることを特徴とする請求項1記載の翼形性能の推定方法。
  3. 前記定数αは、複数の異なる既知の翼形の最大揚力の計算値と実験値との相関が最良になる値に選ばれることを特徴とする請求項1または2記載の翼形性能の推定方法。
  4. 前記定数αは、気流の流れのレイノルズ数に応じて異なる値に選ばれることを特徴とする請求項1,2または3記載の翼形性能の推定方法。
  5. 前記定数αは、複数のレイノルズ数に基づいて求めた値を、簡易式で近似した値に置換えることを特徴とする請求項1,2または3記載の翼形性能の推定方法。
  6. 翼面近傍の渦度に漸近する速度勾配絶対値Sを生成項とする乱流モデルを用いて翼形性能を求めるコンピュータを、
    翼形性能を生成するために必要なデータを入力するための入力手段、
    入力手段によって入力されたデータに基づいて、2つの変数a,bの大なる方を値とする関数をMax(a,b)とし、翼面に沿う速度の翼面に沿う方向の速度勾配をUxとし、壁面に垂直な方向の速度勾配をUyとし、Uyに漸近する渦度値等をΩとし、定数をαとし、座標軸Xは翼に対する気流の相対的な速度Uの方向を正としたとき、
    S=Max(|Ω|−α|Ux|,0)
    または、
    S=Max(|Ω|−αMax(0,Ux),0)
    によって修正された速度勾配絶対値Sを用いて、翼の迎角に対する揚力を求める演算手段、および
    演算手段によって求められた前記翼の迎角に対する揚力を表示する出力手段、として機能させるコンピュータによって実行可能な翼形性能の推定プログラム。
  7. 壁面に垂直な方向の速度をVとしたとき、その垂直方向の勾配はVyで表わし、流体の性質から、壁面近傍では、
    Ux+Vy=0
    が概ね成立することを利用し、前記翼面に沿う方向の速度勾配Uxを翼面に垂直な方向の速度勾配−Vyで置換えることを特徴とする請求項6記載のコンピュータによって実行可能な翼形性能の推定プログラム。
  8. 前記定数αは、複数の異なる既知の翼形の最大揚力の計算値と実験値の相関が最良になる値に選ばれることを特徴とする請求項6または7記載のコンピュータによって実行可能な翼形性能の推定プログラム。
  9. 前記定数αは、気流の流れのレイノルズ数に応じて異なる値に選ばれることを特徴とする請求項6,7または8記載のコンピュータによって実行可能な翼形性能の推定プログラム。
  10. 前記定数αは、複数のレイノルズ数に基づいて求めた値を、簡易式で近似した値に置換えることを特徴とする請求項6,7または8記載のコンピュータによって実行可能な翼形性能の推定プログラム。
JP2002164867A 2002-06-05 2002-06-05 翼形性能の推定方法および翼形性能の推定プログラム Expired - Fee Related JP3587827B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2002164867A JP3587827B2 (ja) 2002-06-05 2002-06-05 翼形性能の推定方法および翼形性能の推定プログラム

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2002164867A JP3587827B2 (ja) 2002-06-05 2002-06-05 翼形性能の推定方法および翼形性能の推定プログラム

Publications (2)

Publication Number Publication Date
JP2004012248A JP2004012248A (ja) 2004-01-15
JP3587827B2 true JP3587827B2 (ja) 2004-11-10

Family

ID=30432896

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2002164867A Expired - Fee Related JP3587827B2 (ja) 2002-06-05 2002-06-05 翼形性能の推定方法および翼形性能の推定プログラム

Country Status (1)

Country Link
JP (1) JP3587827B2 (ja)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4134132B2 (ja) * 2005-09-28 2008-08-13 社団法人日本航空宇宙工業会 ブレード翼型の設計方法
JP5015559B2 (ja) * 2006-11-22 2012-08-29 デンソーテクノ株式会社 管路内流れの圧力損失の評価用プログラム及び装置
US8137066B2 (en) * 2009-04-16 2012-03-20 Frontier Wind, Llc Pressure based load measurement
FR2974062B1 (fr) * 2011-04-13 2013-05-03 Onera (Off Nat Aerospatiale) Rotor de voilure tournante et pale pour un tel rotor
JP7102741B2 (ja) * 2018-01-15 2022-07-20 株式会社豊田中央研究所 流体解析装置、流体解析方法、及び流体解析プログラム
JP7471966B2 (ja) 2020-09-03 2024-04-22 株式会社東芝 乱流数値解析方法、乱流数値解析プログラムおよび乱流数値解析装置
CN114354147B (zh) * 2021-10-22 2023-06-20 中国华能集团清洁能源技术研究院有限公司 一种风力发电机组叶片环境损伤试验系统及其方法和应用

Also Published As

Publication number Publication date
JP2004012248A (ja) 2004-01-15

Similar Documents

Publication Publication Date Title
US8457939B2 (en) Generating inviscid and viscous fluid-flow simulations over an aircraft surface using a fluid-flow mesh
US9348956B2 (en) Generating a simulated fluid flow over a surface using anisotropic diffusion
Biancolini et al. Static aeroelastic analysis of an aircraft wind-tunnel model by means of modal RBF mesh updating
WO2012129440A2 (en) Predicting transition from laminar to turbulent flow over a surface
Abobaker et al. Effect of mesh type on numerical computation of aerodynamic coefficients of NACA 0012 airfoil
Lopez-Morales et al. Verification and validation of HiFiLES: A high-order LES unstructured solver on multi-GPU platforms
López et al. Verification and Validation of HiFiLES: a High-Order LES unstructured solver on multi-GPU platforms
Pandya et al. Accuracy, Scalability, and Efficiency of Mixed-Element USM3D for Benchmark Three-Dimensional Flows
JP3587827B2 (ja) 翼形性能の推定方法および翼形性能の推定プログラム
Ganesh Ram et al. Design Optimization and Analysis of NACA 0012 Airfoil using Computational fluid dynamics and Genetic algorithm
Bird et al. Leading edge vortex formation on finite wings using vortex particles
Mei et al. Implicit numerical simulation of transonic flow through turbine cascades on unstructured grids
Lavoie et al. A penalization method for 2d ice accretion simulations
Lei et al. CFD validation for high-lift devices: two-element airfoil
Aliaga et al. Automatic Mesh Optimization for Wing-Fuselage Juncture Flow Separation Predictions
Nakashima et al. Rans solutions on high lift common research model with scFLOW, a polyhedral finite-volume solver
El Maani et al. CFD Analysis of the Transonic Flow over a NACA 0012 Airfoil
CN113158339A (zh) 一种针对sst湍流模型的湍流长度尺度修正方法
Ruffin et al. Adaptation of a k-epsilon model to a cartesian grid based methodology
Sadrehaghighi Dynamic & Adaptive Meshing
Sugaya et al. Grid metrics modification approach for flow simulation around 3D geometries on Cartesian CFD method
Morton et al. Computational aircraft and armament stability and control techniques applied to the F-16
Delgado-Gutiérrez et al. An efficient implementation of the graphics processing unit-accelerated single-step and simplified lattice Boltzmann method for irregular fluid domains
Emory et al. Componentality-based wall-blocking for RANS models
Sukas et al. HEMLAB Algorithm Applied to 4th AIAA CFD High Lift Prediction Workshop

Legal Events

Date Code Title Description
A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20040617

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20040803

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20040810

R150 Certificate of patent or registration of utility model

Free format text: JAPANESE INTERMEDIATE CODE: R150

Ref document number: 3587827

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20080820

Year of fee payment: 4

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20090820

Year of fee payment: 5

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20100820

Year of fee payment: 6

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110820

Year of fee payment: 7

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110820

Year of fee payment: 7

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120820

Year of fee payment: 8

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120820

Year of fee payment: 8

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130820

Year of fee payment: 9

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130820

Year of fee payment: 9

LAPS Cancellation because of no payment of annual fees