特許 7471966 乱流数値解析方法、乱流数値解析プログラムおよび乱流数値解析装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-04-12

(45)【発行日】2024-04-22

(54)【発明の名称】乱流数値解析方法、乱流数値解析プログラムおよび乱流数値解析装置

(51)【国際特許分類】

   G06F 30/20 20200101AFI20240415BHJP

   G06F 113/08 20200101ALN20240415BHJP

【ＦＩ】

G06F30/20

G06F113:08

【請求項の数】 13

(21)【出願番号】P 2020148421

(22)【出願日】2020-09-03

(65)【公開番号】P2022042813

(43)【公開日】2022-03-15

【審査請求日】2023-02-20

(73)【特許権者】

【識別番号】000003078

【氏名又は名称】株式会社東芝

(73)【特許権者】

【識別番号】317015294

【氏名又は名称】東芝エネルギーシステムズ株式会社

(74)【代理人】

【識別番号】100091487

【弁理士】

【氏名又は名称】中村 行孝

(74)【代理人】

【識別番号】100105153

【弁理士】

【氏名又は名称】朝倉 悟

(74)【代理人】

【識別番号】100107582

【弁理士】

【氏名又は名称】関根 毅

(74)【代理人】

【識別番号】100150717

【弁理士】

【氏名又は名称】山下 和也

(72)【発明者】

【氏名】中島 峻浩

(72)【発明者】

【氏名】黒澤 貞男

(72)【発明者】

【氏名】内田 竜朗

【審査官】堀井 啓明

(56)【参考文献】

【文献】特開２００４－０１２２４８（ＪＰ，Ａ）

【文献】米国特許第０５９２６３９９（ＵＳ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０６Ｆ ３０／００－３０／２８

(57)【特許請求の範囲】

【請求項1】

翼の乱流モデルの第１の流れ場解析を、解析領域を複数のメッシュに分割して行い、各々の前記メッシュにおける物理量を算出するステップと、
算出された前記物理量に基づいて補正関数を作成するステップと、
作成された前記補正関数を用いて前記物理量に対応する補正値を算出し、前記補正値を用いて渦粘性係数を算出するステップと、
算出された前記渦粘性係数を用いて、前記翼の前記乱流モデルの第２の流れ場解析を行うステップと、を備えた、乱流数値解析方法。

【請求項2】

前記補正関数を作成するステップにおいて、前記補正関数は、前記物理量としての圧力勾配に基づいて作成され、
前記渦粘性係数を算出するステップにおいて、作成された前記補正関数を用いて前記圧力勾配に対応する前記補正値が算出される、請求項１に記載の乱流数値解析方法。

【請求項3】

前記補正関数を作成するステップにおいて、前記補正関数は、前記物理量としての流線の曲率または前記流線の曲率半径に基づいて作成され、
前記渦粘性係数を算出するステップにおいて、前記補正関数が前記曲率に基づいて作成された場合、作成された前記補正関数を用いて前記曲率に対応する前記補正値が算出され、前記補正関数が前記曲率半径に基づいて作成された場合、作成された前記補正関数を用いて前記曲率半径に対応する前記補正値が算出される、請求項１に記載の乱流数値解析方法。

【請求項4】

前記補正関数を作成するステップにおいて、前記補正関数は、前記物理量としての流速に基づいて作成され、
前記渦粘性係数を算出するステップにおいて、作成された前記補正関数を用いて前記流速に対応する前記補正値が算出される、請求項１に記載の乱流数値解析方法。

【請求項5】

前記補正関数を作成するステップにおいて、前記補正関数は、前記物理量としての圧力に基づいて作成され、
前記渦粘性係数を算出するステップにおいて、作成された前記補正関数を用いて前記圧力に対応する前記補正値が算出される、請求項１に記載の乱流数値解析方法。

【請求項6】

前記補正関数を作成するステップにおいて、前記補正関数は、前記物理量としての流れのせん断応力に基づいて作成され、
前記渦粘性係数を算出するステップにおいて、作成された前記補正関数を用いて前記せん断応力に対応する前記補正値が算出される、請求項１に記載の乱流数値解析方法。

【請求項7】

前記補正関数を作成するステップにおいて、前記補正関数は、前記物理量としての乱流運動エネルギに基づいて作成され、
前記渦粘性係数を算出するステップにおいて、作成された前記補正関数を用いて前記乱流運動エネルギに対応する前記補正値が算出される、請求項１に記載の乱流数値解析方法。

【請求項8】

前記補正関数を作成するステップにおいて、前記補正関数は、前記物理量としての前記メッシュの座標に基づいて作成され、
前記渦粘性係数を算出するステップにおいて、作成された前記補正関数を用いて前記座標に対応する前記補正値が算出される、請求項１に記載の乱流数値解析方法。

【請求項9】

前記補正関数は、前記物理量が第１の範囲に含まれるときに、前記補正値が第１補正値であり、前記物理量が前記第１の範囲とは異なる第２の範囲に含まれるときに、前記補正値が前記第１補正値とは異なる第２補正値であるように、構成され、
前記第１補正値および前記第２補正値の少なくとも一方は、一定値である、請求項１～８のいずれか一項に記載の乱流数値解析方法。

【請求項10】

前記補正関数は、前記物理量が大きくなるに従って前記補正値が小さくなる、または大きくなるように、構成される、請求項１～８のいずれか一項に記載の乱流数値解析方法。

【請求項11】

前記渦粘性係数を算出するステップにおいて、前記補正値を用いることなく算出される場合の前記渦粘性係数の数式に前記補正値を乗じることにより、前記渦粘性係数が算出される、請求項９または１０に記載の乱流数値解析方法。

【請求項12】

請求項１～１１のいずれか一項に記載された乱流数値解析方法をコンピュータに実行させる乱流数値解析プログラム。

【請求項13】

翼の乱流モデルの第１の流れ場解析を、解析領域を複数のメッシュに分割して行い、各々の前記メッシュにおける物理量を算出する解析部と、
前記解析部で算出された前記物理量に基づいて補正関数を作成する関数作成部と、
前記関数作成部により作成された前記補正関数を用いて前記物理量に対応する補正値を算出し、前記補正値を用いて渦粘性係数を算出する係数算出部と、を備え、
前記解析部は、前記係数算出部で算出された前記渦粘性係数を用いて、前記乱流モデルの第２の流れ場解析を行う、乱流数値解析装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明の実施の形態は、乱流数値解析方法、乱流数値解析プログラムおよび乱流数値解析装置に関する。

【背景技術】

【0002】

フランシス水車等のような流体機械は、非設計点において、翼周りの流れの剥離が発生し得る。この場合、流れが乱れて、効率が低下し得る。このような水車の開発等では、通常、非定常解析手法と比べて解析リソースを少なくすることができる定常解析手法が使用される。しかしながら、定常解析手法では、流れの剥離を精度良く予測することが困難である。このため、非設計点における効率の予測精度が低下し得る。

【0003】

従来の解析手法では、翼表面からの流れの剥離の発生と、剥離領域を適切に予測することが困難である。その理由としては、剥離発生点において渦粘性係数が過大に予測されることが挙げられる。例えば、翼表面の近傍に形成される境界層内において、渦粘性係数が過大に予測され得る。

【0004】

これに対して近年では、補正値（または補正係数）を用いて渦粘性係数を算出する解析手法が知られている。このことにより、渦粘性係数が過大に予測されることの防止を図っている。しかしながら、補正値が、解析領域全体にわたって一定になっている場合がある。この場合、渦粘性係数を過大に予測していない位置においても、渦粘性係数が補正される。その結果、当該位置における渦粘性係数が過小に予測される。

【0005】

これに対して、剥離の予測精度を向上させるための解析手法も知られている。例えば、渦粘性係数の補正値を、翼表面からの距離に応じて変化させる解析手法が知られている。しかしながら、翼表面からの距離に応じて補正値を変化させる場合、予測精度は、翼の形状の影響を受けやすい。この場合、ある翼の形状では、高い予測精度を得ることができるが、他の翼の形状では、予測精度が低下し得る。

【先行技術文献】

【非特許文献】

【0006】

【文献】Suranaree University of Technology, T.Chitsomboon C.Thamthae, 2011, ”Adjustment of k-ω SST turbulence model for an improved prediction of stalls on wind turbine blades”, World Renewable Energy Congress 2011

【発明の概要】

【発明が解決しようとする課題】

【0007】

実施の形態は、このような点を考慮してなされたものであり流れの剥離の予測精度を向上させることができる乱流数値解析方法、乱流数値解析プログラムおよび乱流数値解析装置を提供することを目的とする。

【課題を解決するための手段】

【0008】

実施の形態による乱流数値解析方法は、翼の乱流モデルの第１の流れ場解析を行うステップと、補正関数を作成するステップと、渦粘性係数を算出するステップと、翼の乱流モデルの第２の流れ場解析を行うステップと、を備えている。第１の流れ場解析を行うステップにおいて、各メッシュにおける物理量が算出される。補正関数を作成するステップにおいて、算出された物理量に基づいて補正関数が作成される。渦粘性係数を算出するステップにおいて、作成された補正関数を用いて物理量に対応する補正値が算出され、補正値を用いて渦粘性係数が算出される。第２の流れ場解析を行うステップにおいて、算出された渦粘性係数が用いられる。

【0009】

実施の形態による乱流数値解析プログラムは、上述した乱流数値解析方法をコンピュータに実行させるプログラムである。

【0010】

実施の形態による乱流数値解析装置は、翼の乱流モデルの第１の流れ場解析を行う解析部と、補正関数を作成する関数作成部と、渦粘性係数を算出する係数算出部と、を備えている。解析部は、各メッシュにおける物理量を算出する。関数作成部は、解析部で算出された物理量に基づいて補正関数を作成する。係数算出部は、関数作成部により作成された補正関数を用いて物理量に対応する補正値を算出し、補正値を用いて渦粘性係数を算出する。解析部は、係数算出部で算出された渦粘性係数を用いて、乱流モデルの第２の流れ場解析を行う。

【発明の効果】

【0011】

実施の形態によれば、流れの剥離の予測精度を向上させることができる。

【図面の簡単な説明】

【0012】

【図1】図１は、本実施の形態におけるフランシス水車を示す子午面断面図である。

【図2】図２は、翼の周囲の流れを示す平面断面図である。

【図3】図３は、本実施の形態における乱流数値解析装置の構成を示すブロック図である。

【図4】図４は、本実施の形態による乱流数値解析方法において、翼の負圧面における圧力勾配の分布を示すグラフである。

【図5】図５は、本実施の形態による乱流数値解析方法において、補正関数の一例を示すグラフである。

【図6】図６は、本実施の形態による乱流数値解析方法を示すフローチャートである。

【図7】図７は、本実施の形態による乱流数値解析方法において、渦粘性係数の分布を示すグラフである。

【図8】図８は、本実施の形態による乱流数値解析方法において得られる、迎え角と揚力係数との関係を示すグラフである。

【図9】図９は、本実施の形態による変形例として、翼の負圧面近傍における流線の曲率の分布を示すグラフである。

【図10】図１０は、本実施の形態による変形例として、翼の負圧面近傍における流速の分布を示すグラフである。

【図11】図１１は、本実施の形態による変形例として、翼の負圧面における圧力の分布を示すグラフである。

【図12】図１２は、本実施の形態による変形例として、翼の負圧面における流れのせん断応力の分布を示すグラフである。

【図13】図１３は、本実施の形態による変形例として、翼の負圧面における乱流運動エネルギの分布を示すグラフである。

【発明を実施するための形態】

【0013】

以下、図面を参照して、本発明の実施の形態における乱流数値解析方法、乱流数値解析プログラムおよび乱流数値解析装置について説明する。

【0014】

まず、図１～図８を用いて、本実施の形態における乱流数値解析方法、乱流数値解析プログラムおよび乱流数値解析装置について説明する。乱流数値解析は、種々の流体機械の翼を対象とすることができる。流体機械の一例としては、水力機械が挙げられる。ここでは、まず、図１を用いて水力機械の一例であるフランシス水車について説明する。

【0015】

図１に示すように、フランシス水車１は、水車運転時に上池から水圧鉄管（いずれも図示せず）を通って水が流入する渦巻き状のケーシング２と、複数のステーベーン３と、複数のガイドベーン４と、ランナ５と、を備えている。

【0016】

ステーベーン３は、ケーシング２に流入した水をガイドベーン４およびランナ５に導くための部材である。ステーベーン３は、周方向に所定の間隔をあけて配置されている。ステーベーン３の間に水が流れる流路が形成されている。

【0017】

ガイドベーン４は、流入した水をランナ５に導くための部材である。ガイドベーン４は、周方向に所定の間隔をあけて配置されている。ガイドベーン４の間には、水が流れる流路が形成されている。各ガイドベーン４は、回動可能に構成されており、各ガイドベーン４が回動して開度を変えることにより、ランナ５に流入する水の流量が調整可能になっている。このようにして、後述する発電機７の発電量が調整可能になっている。

【0018】

ランナ５は、ケーシング２に対して回転軸線を中心に回転可能に構成されている。ランナ５は、水車運転時にケーシング２から流入する水によって回転駆動される。すなわち、ランナ５は、ランナ５に流入する水の圧力エネルギを回転エネルギへと変換するための部材である。

【0019】

ランナ５は、後述する主軸６に連結されたクラウン５ａと、クラウン５ａの外周側に設けられたバンド５ｂと、クラウン５ａとバンド５ｂとの間に設けられた複数のランナ羽根５ｃと、を有している。このうちランナ羽根５ｃは、周方向に所定の間隔を開けて配置されている。ランナ羽根５ｃの間には、水が流れる流路が形成されている。

【0020】

ランナ５には、主軸６を介して発電機７が連結されている。この発電機７は、水車運転時には、ランナ５の回転エネルギが伝達されて発電を行うように構成されている。

【0021】

なお、発電機７は、電動機としての機能をも有し、電力が供給されることによりランナ５を回転駆動するように構成されていてもよい。この場合、吸出し管８を介して下池の水を吸い上げて上池に放出させることができ、フランシス水車１を、ポンプ水車としてポンプ運転（揚水運転）することが可能になる。この際、ガイドベーン４の開度は、ポンプ揚程に応じて適切な揚水量になるように変えられる。

【0022】

ランナ５の水車運転時の下流側には、吸出し管８が設けられている。この吸出し管８は、図示しない下池または放水路に連結されており、ランナ５を回転駆動した水が、圧力を回復して、下池または放水路に放出されるようになっている。

【0023】

次に、図３を用いて、本実施の形態による乱流数値解析装置１０について説明する。ここでは、図２に示すような翼２０を乱流数値解析の対象とする。

【0024】

図２に示す翼２０には、流入方向Ｄの流れが、翼２０に流入する。翼２０は、流入方向Ｄに対して所定の角度（迎え角という）をなしている。このことにより、翼２０の正圧面２１よりも負圧面２２の近傍では、流れの剥離が生じやすい。このため、負圧面２２の近傍では、乱流流れが形成される。本実施の形態による乱流数値解析装置１０は、このような乱流流れの数値解析を行う装置である。翼２０は、上述したガイドベーン４であってもよいが、これに限られることはなく、任意の流体機械の翼であってもよい。以下の説明では、一般的な翼２０を乱流数値解析の対象とする。

【0025】

本実施の形態による乱流数値解析装置１０は、入力部１１と、記憶部１２と、演算部１３と、出力部１７と、を備えている。乱流数値解析装置１０は、例えば、乱流数値解析プログラムがインストールされているコンピュータであってもよい。

【0026】

入力部１１は、図示しない入力装置から入力されたデータを取得する機能を有している。例えば、入力部１１は、翼２０の形状および解析条件等を取得する。入力装置は、ユーザの入力操作を受け付ける装置であり、例えば、キーボードまたはマウスであってもよい。

【0027】

記憶部１２は、本実施の形態による乱流数値解析プログラムを記憶する機能を有している。乱流数値解析プログラムは、記憶部１２に記憶されている。記憶部１２は、例えば、メモリ、ハードディスク、ＣＤ－ＲＯＭ、ＤＶＤ－ＲＯＭなどであってもよい。また、記憶部１２は、演算部１３により得られた解析結果を記憶する機能を有していてもよい。

【0028】

演算部１３は、例えば、マイクロプロセッサ等で構成された演算装置である。演算部１３は、乱流数値解析装置１０の種々の動作を制御する。演算部１３は、記憶部１２に記憶された乱流数値解析プログラムを読み出して実行することにより、乱流数値解析方法を実行する。なお、乱流数値解析プログラムを記録したコンピュータ読み取り可能な記録媒体（図示せず）から乱流数値解析プログラムを記憶部１２にインストールしてもよい。また、乱流数値解析プログラムをネットワークからダウンロードして記憶部１２にインストールしてもよい。

【0029】

演算部１３は、解析部１４と、関数作成部１５と、係数算出部１６と、を有している。

【0030】

解析部１４は、翼２０の乱流モデルの流れ場解析を行う機能を有している。解析部１４は、第１の流れ場解析および第２の流れ場解析を含む複数回の流れ場解析を、解析が収束するまで行う。

【0031】

乱流モデルとしては種々のモデルが知られているが、本実施の形態においては、ＲＡＮＳ（Reynolds Averaged Navier-Stokes Simulation）を乱流モデルとして用いる。ＲＡＮＳは、時間平均モデルを用いる乱流モデルである。

【0032】

ＲＡＮＳは、ｋ－ωモデル、ｋ－εモデルおよびＳＳＴｋ－ωモデル等のいくつかのモデルに分類される。一例として、ｋ－ωモデルにおける輸送方程式を、以下の式（１）および式（２）に示す。

【数1】

【数2】

ここで、ρは流体の密度を示し、ｋは乱流運動エネルギを示し、ωは乱流の比散逸率を示している。また、Ｇ_ｋは乱流運動エネルギｋの生成を表す項であり、Ｇ_ωは、ωの発生を表す項であり、Γ_ｋとΓ_ωはそれぞれｋとωの有効拡散係数を表す。Ｙ_ｋとＹ_ωはそれぞれ乱流によるｋとωの散逸を表す項であり、Ｄ_ωはクロス拡散項を表し、Ｓ_ｋとＳ_ωはソース項を表している。

【0033】

解析部１４は、解析領域を複数のメッシュに分割し、上記式（１）および式（２）を用いて乱流モデルの流れ場解析を複数回行う。流れ場解析を行うことによって、各メッシュにおける種々の物理量が解析結果として算出される。物理量には、例えば、圧力勾配、流線の曲率、流速、圧力、せん断応力、乱流運動エネルギおよびメッシュの座標等が含まれている。

【0034】

解析部１４は、後述する係数算出部１６により算出された渦粘性係数を用いて、上述した乱流モデルの流れ場解析を行う。そして、解析部１４は、各メッシュにおける種々の物理量を算出する。

【0035】

関数作成部１５は、解析部１４で算出された物理量に基づいて補正関数を作成する。ここでは、物理量としての圧力勾配に基づいて補正関数を作成する例について説明する。

【0036】

例えば、図４に示すように、翼２０の負圧面２２おける圧力勾配の大きさは、翼長さ方向位置で異なっている。ここでは、流れの方向に対して圧力が上昇する逆圧力勾配の例が示されている。この圧力勾配の大きさに基づいて、補正関数が作成される。圧力勾配は、流線に沿う方向における勾配であってもよい。なお、図４においては、横軸に翼長さ方向位置を示している。翼長さ方向位置とは、翼２０の負圧面２２に沿った前縁２３からの位置を示している。図４の左端は前縁２３の位置（以下、前縁位置と記す）を示し、右端は後縁２４の位置を示している。前縁２３とは、流れの上流側の縁であって、ステーベーン３の側の縁に相当する。後縁２４とは、流れの下流側の縁であって、ランナ５の側の縁に相当する。

【0037】

図５に示すように、補正関数は、圧力勾配が第１の範囲Ｘ_１に含まれるときに補正値が第１補正値ｆ_１であるように構成されていてもよい。また、圧力勾配が第２の範囲Ｘ_２に含まれるときに補正値が第２補正値ｆ_２であってもよく、圧力勾配が第３の範囲Ｘ_３に含まれるときに補正値が第３補正値ｆ_３であってもよい。

【0038】

第２の範囲Ｘ_２は、第１の範囲Ｘ_１とは異なっていてもよい。本実施の形態においては、第２の範囲Ｘ_２は、第１の範囲Ｘ_１よりも小さくなっている。第２補正値ｆ_２は第１補正値ｆ_１とは異なっていてもよい。本実施の形態においては、第２補正値ｆ_２は第１補正値ｆ_１よりも大きくなっている。第３の範囲Ｘ_３は、第２の範囲Ｘ_２とは異なっていてもよい。本実施の形態においては、第３の範囲Ｘ_３は、第２の範囲Ｘ_２よりも小さくなっている。第３補正値ｆ_３は第２補正値ｆ_２とは異なっていてもよい。本実施の形態においては、第３補正値ｆ_３は第２補正値ｆ_２よりも大きくなっている。

【0039】

第１補正値ｆ_１、第２補正値ｆ_２および第３補正値ｆ_３の少なくとも一方は、一定値であってもよい。図５においては、第１補正値ｆ_１、第２補正値ｆ_２および第３補正値ｆ_３がそれぞれ一定値になっている。すなわち、図４に示すように、各メッシュにおける圧力勾配の最小値から最大値の範囲内で２つの閾値を設定し、最小値から最大値の範囲を３つの範囲に区分けして、補正値を設定している。ここでは、第１の範囲Ｘ_１と第２の範囲Ｘ_２と第３の範囲Ｘ_３を、不均等に区分けしているが、均等に区分けしてもよい。

【0040】

このような補正関数ｆは、例えば以下の式（３）で求めることができる。

【数3】

ここで、Ｘは物理量であり、本実施の形態においては圧力勾配に相当している。ｃ_１～ｃ_４およびｃ’_１～ｃ’_４は定数である。ｃ_１～ｃ_４およびｃ’_１～ｃ’_４は、図５に示す補正関数が得られるように設定されてもよい。また、ｃ_１～ｃ_４およびｃ’_１～ｃ’_４は、物理量の種別に応じて設定されてもよい。

【0041】

係数算出部１６は、関数作成部１５により作成された補正関数と、上述した圧力勾配と、に基づいて、渦粘性係数を算出する機能を有している。より具体的には、係数算出部１６は、補正関数を用いて圧力勾配に対応する補正値を算出し、この補正値を用いて渦粘性係数を算出する。本実施の形態における係数算出部１６は、補正値を用いることなく算出される場合の渦粘性係数の数式に、補正値を乗じることにより、渦粘性係数を算出する。以下、補正値を用いることなく算出される場合の渦粘性係数を、非補正渦粘性係数と称する。

【0042】

ｋ－ωモデルの場合には、非補正渦粘性係数μ_ｔｏの数式は、以下の式（４）で表される。

【数4】

【0043】

本実施の形態による係数算出部１６は、以下の式（５）により、補正値を用いて渦粘性係数μ_ｔを算出する。渦粘性係数μ_ｔを補正渦粘性係数と称してもよい。

【数5】

式（５）に示されているように、非補正渦粘性係数に補正値（または補正関数）を乗じて、渦粘性係数が算出される。渦粘性係数は、上述した式（１）および（２）におけるΓ_ｋとΓ_ω、Ｇ_ｋ、およびＧ_ωで用いられる。

【0044】

ｋ－εモデルの場合には、渦粘性係数の算出式は、以下の式（６）で与えられる。

【数6】

ここで、Ｃ_μはモデル定数、εは乱流の散逸率を示す。

【0045】

ＳＳＴｋ－ωモデルの場合には、渦粘性係数の算出式は、以下の式（７）で与えられる。

【数7】

ここで、ａ_１は調整定数、Ｓはひずみ速度、Ｆ_２はブレンド関数を示す。

【0046】

このように、本実施の形態による係数算出部１６は、ｋ－ωモデル以外の他のモデルにも適用することができる。

【0047】

出力部１７は、図３に示すように、演算部１３により得られた解析結果を、図示しない表示装置に出力する機能を有していてもよい。解析結果としては、例えば、流れ場解析により得られた種々の物理量および渦粘性係数であってもよい。表示装置は、例えば、ディスプレイであってもよい。また、出力部１７は、解析結果を、図示しない外部記憶装置に出力する機能を有していてもよい。

【0048】

次に、本実施の形態による乱流数値解析方法について、図６を用いて説明する。

【0049】

まず、ステップＳ１として、対象となる翼２０の形状と解析条件等が入力部１１に入力される。

【0050】

ステップＳ１の後、ステップＳ２として、入力された翼２０の形状と解析条件等を用いて、解析部１４により、乱流モデルの第１の流れ場解析が行われる。この場合、上述した式（１）および式（２）の輸送方程式が解かれる。このことにより、各メッシュにおける種々の物理量が算出される。この物理量には圧力勾配も含まれる。例えば、翼２０の負圧面２２においては、図４に示すように圧力勾配の分布が存在している。

【0051】

ステップＳ２の後、ステップＳ３として、算出された圧力勾配に基づいて補正関数が作成される。例えば、図５に示すような構成を有する補正関数が作成される。この場合、上述した式（３）を用いて補正関数が作成される。

【0052】

ステップＳ３の後、ステップＳ４として、作成された補正関数と、圧力勾配と、に基づいて、渦粘性係数が算出される。ここでは、作成された補正関数を用いて圧力勾配に対応する補正値を算出し、この補正値を用いて渦粘性係数が算出される。より具体的には、補正関数から、各メッシュにおける圧力勾配に対応する補正値が算出される。例えば、図４に示すように、負圧面２２における翼長さ方向位置ｐ_０における圧力勾配がＸ_０とすると、図５に示すように、圧力勾配Ｘ_０に対応する補正値はｆ_２となる。そして、上述した式（５）により、各メッシュにおいて、対応する補正値を用いて渦粘性係数が算出される。

【0053】

ステップＳ４の後、ステップＳ５として、算出された渦粘性係数を用いて、解析部１４により翼２０の乱流モデルの第２の流れ場解析が行われる。そして、解析部１４は、第１の流れ場解析と同様に、各メッシュにおける種々の物理量を算出する。

【0054】

ステップＳ５の後、ステップＳ６として、解析部１４により、解析が収束しているか否かが判断される。

【0055】

解析が収束していると判断された場合には、数値解析を終了する。

【0056】

解析が収束していないと判断された場合には、上述したステップＳ３に戻る。この場合、ステップＳ５の第２の流れ場解析によって得られた圧力勾配を用いて補正関数が再び作成される。そして、ステップＳ４～ステップＳ６が行われる。解析が収束するまで、ステップＳ３～ステップＳ６が繰り返される。

【0057】

ここで、上述のステップＳ４で算出された渦粘性係数を、図７の破線で示す。図７は、横軸に翼長さ方向位置を示しており、縦軸に渦粘性係数を示している。図７の実線は、上述した式（４）で算出される非補正渦粘性係数を示している。図７の破線は、上述した式（５）で算出される渦粘性係数である。

【0058】

図７に示されているように、破線で示される渦粘性係数は、実線で示される非補正渦粘性係数よりも全体的に低減されている。そして、前縁位置において渦粘性係数は比較的大きくなっており、この比較的大きい渦粘性係数が、高い低減率となるように補正される。このことにより、前縁位置における渦粘性係数が、過大に評価されることが抑制されている。一方、比較的小さい渦粘性係数は、小さくなり過ぎないように補正される。言い換えると、比較的小さい渦粘性係数の低減率は低くなっている。

【0059】

通常、流れ方向に対する翼２０の迎え角が大きくなると、翼２０の負圧面２２において、前縁位置から流れの剥離が発生する。このため、図７に実線で示すように、前縁位置における非補正渦粘性係数が、過大に予測される傾向にあり、流れの剥離が過小評価されるという問題がある。

【0060】

これに対して本実施の形態においては、物理量としての圧力勾配を用いて補正関数を作成し、この補正関数を用いて渦粘性係数を補正している。この補正により、図７に破線で示すように、比較的大きい渦粘性係数の低減率を高くして、渦粘性係数が小さくなるように補正することができる。このため、渦粘性係数が過大に予測されることを抑制できる。とりわけ、図４に示すように、翼２０の前縁位置において圧力勾配が急激に大きくなっている。圧力勾配が大きくなると、流れの剥離が生じ得る。流れの剥離が生じた領域において、渦粘性係数が過大に評価される傾向にある。このため、この圧力勾配が大きい領域で、渦粘性係数の低減率を高めることにより、渦粘性係数を効果的に補正することができる。

【0061】

また、本実施の形態においては、翼２０の形状が異なる場合においても、渦粘性係数を適切に補正することができる。すなわち、本実施の形態では、各メッシュにおける渦粘性係数を、対応するメッシュにおける物理量としての圧力勾配を用いて補正している。このため、翼２０の形状が渦粘性係数の補正に影響を及ぼすことを抑制できる。

【0062】

上述した本実施の形態による乱流数値解析方法によって得られたデータの一例を図８に示す。図８に、迎え角と揚力係数との関係が示されている。図８の実線は、非補正渦粘性係数を用いて得られる揚力係数を示している。図８の破線は、本実施の形態による乱流数値解析方法によって渦粘性係数を用いて得られる揚力係数を示している。図８の一点鎖線は、実験により得られたデータを結んだ線である。

【0063】

図８に示すように、抑え角が比較的小さい範囲では、非補正渦粘性係数を用いる場合の揚力係数と、本実施の形態による渦粘性係数を用いる場合の揚力係数はいずれも、実験結果と良好に一致している。しかしながら、抑え角が比較的大きい範囲では、非補正渦粘性係数を用いる場合の揚力係数は、実験結果よりも大きくなっている。これは、上述したように、剥離現象を過小に評価したことにより、剥離による揚力の低減が十分に反映されずに、その結果として揚力係数が大きくなっていることを意味している。これに対して、本実施の形態による渦粘性係数を用いる場合の揚力係数は、実験結果の揚力係数の低下傾向と同様の傾向を示している。すなわち、この渦粘性係数を用いることにより、揚力係数を、実験結果の揚力係数に近づけることができる。このため、剥離の予測精度が向上していることがわかる。なお、図８は、ＳＳＴｋ－ωモデルに本実施の形態による乱流数値解析方法を適用した例であるが、ｋ－ωモデル等のＲＡＮＳの他のモデルでも、同様の結果を得ることができる。

【0064】

このように本実施の形態によれば、第１の流れ場解析によって算出された圧力勾配（または物理量）に基づいて補正関数が作成され、作成された補正関数を用いて圧力勾配に対応する補正値が算出される。そして、この補正値を用いて渦粘性係数が算出される。このことにより、圧力勾配に応じて、渦粘性係数を補正することができる。このため、流れの剥離の予測精度を向上させることができる。とりわけ、本実施の形態によれば、圧力勾配に応じて渦粘性係数が算出される。圧力勾配は、流れの剥離に関連性が高い指標であるため、渦粘性係数の精度を向上させることができる。このため、流れの剥離の予測精度をより一層向上させることができる。

【0065】

また、本実施の形態によれば、圧力勾配が第１の範囲に含まれるときに補正値が第１補正値であり、圧力勾配が第２の範囲に含まれるときに補正値が第２補正値である。そして、第１補正値および第２補正値の少なくとも一方は、一定値である。このことにより、補正関数を簡素化することができ、解析負荷を低減することができる。

【0066】

また、本実施の形態によれば、非補正渦粘性係数の数式に補正値を乗じることにより、渦粘性係数が算出される。このことにより、圧力勾配が比較的大きい場合には、小さい補正値を上述した式（４）の非補正渦粘性係数に乗ずることにより、渦粘性係数を小さくすることができる。一方、圧力勾配が比較的小さい場合には、大きい補正値を上述した式（４）の非補正渦粘性係数に乗ずることにより、渦粘性係数を大きくすることができる。このため、圧力勾配に応じて、渦粘性係数の精度を向上させることができる。

【0067】

なお、上述した本実施の形態においては、補正関数が、物理量としての圧力勾配に基づいて作成される例について説明した。しかしながら、このことに限られることはない。例えば、補正関数は、物理量としての流線の曲率（例えば、曲率の大きさ）に基づいて作成されるようにしてもよい。この場合、作成された補正関数を用いて流線の曲率に対応する補正値が算出される。図９には、一例として、翼２０の負圧面２２の近傍における流線の曲率の分布が示されている。負圧面２２の近傍とは、例えば、負圧面２２に接するメッシュのことを意味している。図９に示すように、翼２０の前縁位置において流線の曲率が急激に大きくなっている。流線の曲率は、流れの剥離に関連性が高い指標であり、流線の曲率が大きくなると、流れの剥離が生じ得る。流れの剥離が生じた領域において、渦粘性係数が過大に評価される傾向にある。このことにより、この流線の曲率が大きい領域で、渦粘性係数の低減率を高めることにより、渦粘性係数を効果的に補正することができる。このため、渦粘性係数の精度を向上させることができる。なお、流線の曲率に代えて、流線の曲率半径を用いてもよい。この場合、曲率半径は、流れの剥離が生じた領域において小さくなることから、物理量として曲率半径を用いる場合には、上述した式（３）に示す「－Ｘ」は「Ｘ」に置き換えてもよい。

【0068】

また、補正関数は、物理量としての流速（例えば、流速の大きさ）に基づいて作成されてもよい。この場合、作成された補正関数を用いて流速に対応する補正値が算出される。図１０には、一例として、翼２０の負圧面２２の近傍における流速の分布が示されている。負圧面２２の近傍とは、例えば、負圧面２２に接するメッシュのことを意味している。図１０に示すように、翼２０の前縁位置において流速が急激に大きくなっている。流速は、流れの剥離に関連性が高い指標であり、流速が大きくなると流れの剥離が生じ得る。流れの剥離が生じた領域において、渦粘性係数が過大に評価される傾向にある。このことにより、この流速が大きい領域で、渦粘性係数の低減率を高めることにより、渦粘性係数を効果的に補正することができる。このため、渦粘性係数の精度を向上させることができる。

【0069】

また、補正関数は、物理量としての圧力（例えば、圧力の大きさ）に基づいて作成されてもよい。この場合、作成された補正関数を用いて圧力に対応する補正値が算出される。図１１には、一例として、翼２０の負圧面２２における圧力の分布が示されている。図１１に示すように、翼２０の前縁位置において圧力が急激に小さくなっている。圧力は、流れの剥離に関連性が高い指標であり、圧力が小さくなると流れの剥離が生じ得る。流れの剥離が生じた領域において、渦粘性係数が過大に評価される傾向にある。このことにより、この圧力が小さい領域で、渦粘性係数の低減率を高めることにより、渦粘性係数を効果的に補正することができる。このため、渦粘性係数の精度を向上させることができる。圧力は、流れの剥離が生じた領域において小さくなっていることから、物理量として圧力の大きさを用いる場合には、上述した式（３）に示す「－Ｘ」は「Ｘ」に置き換えてもよい。

【0070】

また、補正関数は、物理量としての流れのせん断応力（例えば、せん断応力の大きさ）に基づいて作成されてもよい。この場合、作成された補正関数を用いて流れのせん断応力に対応する補正値が算出される。図１２には、一例として、翼２０の負圧面２２における流れのせん断応力の分布が示されている。図１２に示すように、翼２０の前縁位置において流れのせん断応力が急激に大きくなっている。流れのせん断応力は、流れの剥離に関連性が高い指標であり、流れのせん断応力が大きくなると流れの剥離が生じ得る。流れの剥離が生じた領域において、渦粘性係数が過大に評価される傾向にある。このことにより、この流れのせん断応力が大きい領域で、渦粘性係数の低減率を高めることにより、渦粘性係数を効果的に補正することができる。このため、渦粘性係数の精度を向上させることができる。

【0071】

また、補正関数は、物理量としての乱流運動エネルギ（例えば、乱流運動エネルギの大きさ）に基づいて作成されてもよい。この場合、作成された補正関数を用いて乱流運動エネルギに対応する補正値が算出される。乱流運動エネルギは、上述した式（１）等で用いた符号ｋに相当する。図１３には、一例として、翼２０の負圧面２２における乱流運動エネルギの分布が示されている。図１３に示すように、翼２０の前縁位置において乱流運動エネルギが急激に大きくなっている。乱流運動エネルギは、流れの剥離に関連性が高い指標であり、乱流運動エネルギが大きくなると流れの剥離が生じ得る。流れの剥離が生じた領域において、渦粘性係数が過大に評価される傾向にある。このことにより、この乱流運動エネルギが大きい領域で、渦粘性係数の低減率を高めることにより、渦粘性係数を効果的に補正することができる。このため、渦粘性係数の精度を向上させることができる。

【0072】

また、補正関数は、物理量としての座標に基づいて作成されてもよい。この場合、作成された補正関数を用いて座標に対応する補正値が算出される。例えば、解析領域内に原点を任意に設定し、各メッシュにおけるｘ座標値を物理量として用いてもよい。あるいは、各メッシュにおけるｙ座標値またはｚ座標値を物理量として用いてもよい。更には、ｘ座標値、ｙ座標値およびｚ座標値のうちの２つの平均値（または二乗平均値）を物理量としてもよく、またはｘ座標値、ｙ座標値およびｚ座標値の平均値（または二乗平均値）を物理量としてもよい。ｘ座標値、ｙ座標値またはｚ座標値を、ｘ座標値、ｙ座標値およびｚ座標値のうちの少なくとも２つの二乗平均値で除算した値を物理量として用いてもよい。ｘ座標値、ｙ座標値またはｚ座標値を、ｘ座標値、ｙ座標値およびｚ座標値の二乗平均値で除算した値を物理量として用いてもよい。ｘ座標値、ｙ座標値およびｚ座標値の少なくとも１つを用いた値であれば、物理量は任意である。上述したように、翼２０の前縁位置において流れの剥離が発生しやすくなっている。このため、翼２０の前縁位置に比較的近い座標位置で流れの剥離が生じ得る。流れの剥離が生じた座標位置において、渦粘性係数が過大に評価される傾向にある。このことにより、この流れの剥離が生じやすい座標位置で、渦粘性係数の低減率を高めることにより、渦粘性係数を効果的に補正することができる。このため、渦粘性係数の精度を向上させることができる。この流れの剥離が生じやすい座標位置で補正値が大きくなるように、上述した式（５）におけるｃ_１～ｃ_４およびｃ’_１～ｃ’_４を設定してもよい。このことにより、図７に示すような補正関数を作成することができる。なお、流れの流入方向Ｄおよび翼２０の向きは、図２に示す例に限られることはなく、任意である。

【0073】

また、上述した本実施の形態においては、補正関数が、第１補正値と、第２補正値と、第３補正値と、を含むように構成されている例について説明した。しかしながら、このことに限られることはない。例えば、補正関数は、第３補正値を含まなくてもよい。すなわち、圧力勾配の最小値から最大値の範囲を、２つの範囲に区分けして、補正関数を作成してもよい。あるいは、圧力勾配の最小値から最大値の範囲の区分け数は、４つ以上でもよく、任意である。

【0074】

また、上述した本実施の形態においては、補正関数が、圧力勾配が第１の範囲に含まれるときに補正値が第１補正値であり、圧力勾配が第２の範囲に含まれるときに補正値が第２補正値であるよう構成されている例について説明した。しかしながら、このことに限られることはない。例えば、補正関数は、圧力勾配（または物理量）が大きくなるに従って補正値が小さくなるように構成されていてもよい。この場合、圧力勾配が大きくなるに従って補正値が徐々に小さくなるようにしてもよい。補正関数は、線形関数でもよく、または非線形関数でもよい。また、補正関数は、ある圧力勾配の範囲では補正値が一定であるとともに、他の圧力勾配の範囲では圧力勾配が大きくなるに従って補正値が徐々に小さくなる、ように構成されていてもよい。

【0075】

また、上述した本実施の形態においては、翼の乱流モデルとしてＲＡＮＳを用いている例について説明した。しかしながら、このことに限られることはなく、ＲＡＮＳ以外の乱流モデルを、本実施の形態による乱流数値解析のモデルとしてもよい。

【0076】

以上述べた実施の形態によれば、流れの剥離の予測精度を向上させることができる。

【0077】

本発明の実施形態といくつかの変形例を説明したが、これらの実施形態および変形例は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態および変形例は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。また、当然のことながら、本発明の要旨の範囲内で、これらの実施の形態および変形例を、部分的に適宜組み合わせることも可能である。

【符号の説明】

【0078】

１０：乱流数値解析装置、１４：解析部、１５：関数作成部、１６：係数算出部、２０：翼

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】