(57)【要約】【課題】スワールが生じている流れ場の数値解析の精度を向上させる。【解決手段】実施の形態に係る乱流のシミュレーション方法では、スワールの生じる流れ場で測定した乱流の特性の測定結果を取得する。また、乱流モデルと前記乱流モデルに対応した第１の乱流モデル係数とを用いて数値解析を行い、前記流れ場における乱流の特性の解析結果を得る。さらに、測定結果と解析結果との誤差を算出する。次に、誤差に対して相関を持つ所定のパラメータを指定し、流れ場全体の領域にわたって前記誤差と前記パラメータとの相関性を数値化する。次に、流れ場全体の領域を、相関性の値が閾値を超える第１の解析領域と相関性の値が閾値以下の第２の解析領域とに分割する。さらに、乱流モデルと第１の乱流モデル係数の値を補正した第２の乱流モデル係数とを用いて、第１の解析領域における乱流の特性を数値解析する。【選択図】図６
2017004103.tif 000002
【特許請求の範囲】
【請求項１】
スワールの生じる流れ場で測定した乱流の特性の測定結果を取得するステップと、
乱流モデルと前記乱流モデルに対応した第１の乱流モデル係数とを用いて数値解析を行い、前記流れ場における乱流の特性の解析結果を得るステップと、
前記取得された測定結果と前記数値解析して得られた解析結果との誤差を算出するステップと、
前記算出された誤差に対して相関を持つ所定のパラメータを指定し、前記流れ場全体の領域にわたって前記誤差と前記パラメータとの相関性を数値化するステップと、
前記流れ場全体の領域を、前記数値化された相関性の値が閾値を超える第１の解析領域と前記相関性の値が閾値以下の第２の解析領域とに分割するステップと、
前記乱流モデルと前記第１の乱流モデル係数の値を補正した第２の乱流モデル係数とを用いて、前記第１の解析領域における乱流の特性を数値解析するステップと、
を有する乱流のシミュレーション方法。
【請求項２】
前記乱流モデルと前記第１の乱流モデル係数とを用いて、前記第２の解析領域における乱流の特性を数値解析するステップ、
をさらに有する請求項１記載の乱流のシミュレーション方法。
【請求項３】
前記乱流モデルと前記第１の乱流モデル係数の値を段階的に増減させた複数の仮の乱流モデル係数のうちの一つとを順次置き換えながら用いて、前記第１の解析領域における乱流の特性の数値解析を繰り返し、この数値解析の解析結果と前記測定結果との誤差が最も小さくなったときの前記仮の乱流モデル係数の値を、前記第２の乱流モデル係数として検出するステップ、
をさらに有する請求項１又は２記載の乱流のシミュレーション方法。
【請求項４】
前記分割するステップでは、前記流れ場全体の領域を前記スワールの径方向に分割して前記第１及び第２の解析領域を設定する、
請求項１ないし３のいずれか１項に記載の乱流のシミュレーション方法。
【請求項５】
前記測定結果を取得するステップ及び前記解析結果を得るステップでは、前記測定結果及び前記解析結果としてそれぞれ流速を取得し、
前記相関性を数値化するステップでは、前記パラメータとして渦度を指定する、
請求項１ないし４のいずれか１項に記載の乱流のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、乱流のシミュレーション方法に関する。
【背景技術】
【０００２】
近年の電子計算機における計算性能の向上に伴い、流体解析技術も急激に進歩を遂げ、最近では流体機械の開発に流体シミュレーションが多用されている。流体機械の周りの流れは、大部分が乱流となっており、乱流をモデル化することなく直接数値シミュレーション（DNS: Direct Numerical Simulation）を行うことが本来望ましいものの、今日の最速の電子計算機を利用したとしても乱流の直接計算は難しくなっている。
【０００３】
また、乱流の長さスケールと時間スケールは、解析対象領域の中で広範囲におよび、これらのスケールは、メッシュ分割に用いる解析格子のサイズよりもはるかに小さいことが一般的である。そこで、乱流のモデル化が必要となり、モデル式と流体の運動方程式とを連立して解析を行う。なお、数値シミュレーションに利用される乱流モデルは、今日まで数多く提案されている。
【先行技術文献】
【特許文献】
【０００４】
【非特許文献１】乱流解析（数値流体力学シリーズ） 数値流体力学編集委員会，東京大学出版会
【非特許文献２】乱流の数値流体力学 モデルと計算法 大宮司久明，三宅裕，吉澤徴，東京大学出版会
【非特許文献３】Sadao KUROSAWA, Turblent Flow Simulation for the Draft Tube of a Kaplan Turbine, 23rd IAHR Symposium, 2006
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、乱流モデルを用いる既存の数値シミュレーション方法は、例えばスワール（渦）が発生する流れ場を対象とした場合、この流れ場全体の領域において、数値シミュレーションによる解析結果が、実際の測定結果を十分に模擬できているとはいえず、解析精度の改善が求められている。
【０００６】
そこで、本発明が解決しようとする課題は、スワールが生じている流れ場の数値解析の精度を高めることができる乱流のシミュレーション方法を提供することである。
【課題を解決するための手段】
【０００７】
実施の形態に係る乱流のシミュレーション方法は、測定結果を取得するステップ、解析結果を得るステップ、誤差を算出するステップ、数値化するステップ、分割するステップ、及び数値解析するステップを有する。測定結果を取得するステップでは、スワールの生じる流れ場で測定された乱流の特性の測定結果を取得する。解析結果を得るステップでは、乱流モデルと前記乱流モデルに対応した第１の乱流モデル係数とを用いて数値解析を行い、前記流れ場における乱流の特性の解析結果を得る。誤差を算出するステップでは、前記取得された測定結果と前記数値解析を行って得られた解析結果との誤差を算出する。相関性を数値化するステップでは、前記算出された誤差に対して相関を持つ所定のパラメータを指定し、前記流れ場全体の領域にわたって前記誤差と前記パラメータとの相関性を数値化する。分割するステップでは、前記流れ場全体の領域を、前記数値化された相関性の値が閾値を超える第１の解析領域と前記相関性の値が閾値以下の第２の解析領域とに分割する。数値解析するステップでは、前記乱流モデルと前記第１の乱流モデル係数の値を補正した第２の乱流モデル係数とを用いて、前記第１の解析領域における乱流の特性を数値解析する。
【図面の簡単な説明】
【０００８】
【図１】実施の形態に係る乱流シミュレーション装置の構成を示す機能ブロック図。
【図２】図１の乱流シミュレーション装置によって流れ場を解析する対象のカプラン水車の構成を概略的に示す断面図。
【図３】図２のカプラン水車内の流れ場に解析格子を配置した状態を示す鳥瞰図。
【図４】図２中の流れ場で測定した乱流の特性の測定結果とこの流れ場の乱流の特性を従来の数値シミュレーションにより解析した解析結果との関係を示す図。
【図５】図１の乱流シミュレーション装置により解析領域を分割する場合の説明図。
【図６】実施の形態に係る乱流のシミュレーション方法を示すフローチャート。
【発明を実施するための形態】
【０００９】
以下、実施の形態を図面に基づき説明する。
本実施形態に係る乱流のシミュレーション方法は、図１に示すように、例えば、乱流シミュレーション装置２０を用いて実現される。この乱流シミュレーション装置２０は、図１に示すように、測定結果取得部２１、数値解析部２２、条件設定部２３、誤差算出部２４、相関性算出部２５、解析領域分割部２６、閾値記憶部２７、乱流モデル係数変更部２８及び乱流モデル係数確定部２９を備えている。
【００１０】
ここで、本実施形態では、図２に示すように、例えばカプラン水車（Kaplan Turbine）３０内のスワール（渦）が生じる流れ場について数値シミュレーションを行う場合を例示する。なお、カプラン水車３０は、図２に示すように、例えば、渦巻状の水管であるケーシング３１、ガイドベーン３２、ランナーベーン３３ａやランナーコーン３３ｂを含むプロペラ状の羽根車（ランナー）３３、吸出し管（ドラフトチューブ）３４などを備えている。
【００１１】
図１、図２に示すように、上述した測定結果取得部２１は、吸出し管３４内のスワールの生じる流れ場において、例えば熱線流速計やレーザ流速計などを用いて測定（実測）された乱流の特性の測定結果（流速や圧力などの測定結果）を取得する。この測定結果取得部２１は、一例として、吸出し管３４内における羽根車３３の例えばランナーコーン３３ｂの直下の領域で測定された軸流速度（羽根車３３の軸方向に沿った流速）などを取得する。
【００１２】
一方、数値解析部２２は、所定の乱流モデルとこの乱流モデルに対応した乱流モデル係数（第１の乱流モデル係数）とを少なくとも用いて数値解析を行い、スワールが生じる流れ場における乱流の特性の解析結果を得る。数値解析部２２は、一例として、上記ランナーコーン３３ｂの底部近傍の領域について数値解析した軸流速度などを取得する。
【００１３】
ここで、乱流モデルを適用する数値シミュレーションは、一般に、渦粘性モデルや応力方程式モデル（RSM: Reynolds Stress Model）などの時間平均モデルを用いる定常解析向きのＲＡＮＳ（Reynolds Averaged Navier-Stokes Simulation）と、空間平均モデルを用いる非定常解析も可能なＬＥＳ（Large-Eddy Simulation）などに分類できる。後者のＬＥＳは、例えば剥離を伴う流れ解析などについては、ＲＡＮＳよりも良好な解析結果が得られるものの、解の格子依存性や計算負荷が高いことが知られている。
【００１４】
一方、前者のＲＡＮＳは、渦粘性モデルとしては、例えばBaldwin-Lomaxモデルなどの０方程式モデル、Spalart-Allmarasモデルなどの１方程式モデル、ｋ−εモデルやｋ−ωモデルなどの２方程式モデルを例示できる。また、レイノルズ応力の輸送方程式を解く上記した応力方程式モデルとしては、Gibson-LaunderモデルやSpeziale-Sarkar-Gatskiモデルなどを挙げることができる。
【００１５】
上述した０方程式モデルは、輸送方程式を解かないため計算負荷が低いものの、その一方で剥離を伴う流れ解析などが過小評価となる。また、１方程式モデルは、渦動粘性の輸送方程式のみを解くものであり、剥離を伴う流れ解析などが過大評価となる。さらに、２方程式モデルは、運動エネルギｋと散逸率ε又は比散逸率ω（ε／ｋ）の輸送方程式を解くものであり、剥離、旋回、二次流れの解析に不向きである。一方、応力方程式モデルは、渦粘性モデルよりも剥離、旋回、二次流れの解析結果が良好に得られる。
【００１６】
なお、標準ｋ−εモデルの場合、運動エネルギｋと散逸率εの輸送方程式は以下の式で定義される。
【００１７】
【数１】
2017004103.tif 000003
【００１８】
【数２】
2017004103.tif 000004
【００１９】
【数３】
2017004103.tif 000005
【００２０】
ここで、上記の式１〜式３において、ｕは速度（流速）、ρは流体の密度、Ｅは変形速度、μは渦粘係数、μtは渦粘性を示している。また、Ｃμ、σk、σε、Ｃ1ε、及びＣ2εは、標準ｋ−εモデルに対応する乱流モデル係数（第１の乱流モデル係数）を示している。乱流シミュレーション装置２０の数値解析部２２は、例えば、この標準ｋ−εモデルとそれぞれ固定値である複数の乱流モデル係数（複数の第１の乱流モデル係数）とを少なくとも用いて数値解析を行う機能を備えている。
【００２１】
【数４】
2017004103.tif 000006
【００２２】
図３は、カプラン水車３０の吸出し管３４内の領域を対象として、数値解析部２２によって解析格子（メッシュ）３５が形成された状態を例示している。数値解析部２２は、解析結果に悪影響を及ぼすといわれる捩れ度が極力少なくなるように直交性を高くして、格子点数が例えば約１００万点の解析格子３５を形成する。
【００２３】
図１、図３に示すように、条件設定部２３は、入口境界Ｅ１、出口境界Ｅ２、壁境界Ｅ３それぞれについて条件を与える。つまり、条件設定部２３は、入口境界Ｅ１の境界条件には、測定結果取得部２１により取得された解析領域の入口部分における流体の速度分布（図２、図３中のＥ１部分における流体の軸流速度［Axial Velocity］及びＥ１部分における流れの接線方向成分の速度［Tangential Velocity］の分布など）を境界条件として設定する。また、条件設定部２３は、出口境界Ｅ２の境界条件には、測定結果取得部２１により取得された解析領域の出口部分における流体の圧力の計測結果を設定し、さらに、壁境界Ｅ３には、壁関数境界条件を設定する。
【００２４】
また、数値解析部２２は、解析手法として、離散化手法には例えば有限体積法を用い、対流項には例えばＱＵＩＣＫ法を用い、差分解法には圧力と速度とを統合した例えばSimple法などを適用する。さらに、数値解析部２２は、乱流モデルには、例えばＲＳＭ（Reynolds Stress Model）や、標準ｋ−εモデルに属する例えばCubic ｋ−εモデルなどを適用した数値解析を行う機能を有する。
【００２５】
図４は、羽根車３３（ランナーコーン３３ｂ）の直下の流れ領域で測定した乱流の特性の測定結果とこの流れ領域の乱流の特性を従来の数値シミュレーション（乱流モデルと固定値の乱流モデル係数とを用いた従来の数値解析）により解析した解析結果との関係を示す図である。具体的には、図４は、羽根車３３の直下の上記流れ領域における羽根車３３の径方向の位置と、羽根車３３の直下の上記流れ領域における軸流速度の測定結果並びにCubic ｋ−εモデル及びＲＳＭによる当該軸流速度の解析結果と、の関係を示している。
【００２６】
図４に示すように、羽根車３３（ランナー）の回転中心から外径側へ約７０ｍｍまでの位置に対応するルート部（回転機械の基端部分）と、この回転中心を基準として羽根車３３の外径側約１５０ｍｍ〜２６０ｍｍまでの範囲に対応するチップ部（回転機械の先端部分）とでは、Cubic ｋ−εモデル及びＲＳＭのいずれの解析結果も、測定結果（実測した計測結果）を適切に模擬できていないことが分かる。
【００２７】
つまり、羽根車の直下（回転翼の下流）のようなスワールが発生する流れ場において、所定の乱流モデルとこの乱流モデル対応する固定値の乱流モデル係数とを用いる既存の数値シミュレーション方法は、定性的に測定結果を模擬できたとしても、定量的には測定結果を模擬することが難しくなっている。カプラン水車に限らず、回転機械におけるスワールの生じる流れ場の解析では、同様な結果になることが想定される。
【００２８】
そこで、本実施形態の乱流シミュレーション装置２０は、図１に示すように、上述した誤差算出部２４、相関性算出部２５、解析領域分割部２６、閾値記憶部２７、乱流モデル係数変更部２８及び乱流モデル係数確定部２９をさらに備えている。
【００２９】
誤差算出部２４は、測定結果取得部２１により取得された流速などの測定結果と数値解析部２２により数値解析して得られた流速などの解析結果（所定の乱流モデルとこの乱流モデルに対応する固定値の第１の乱流モデル係数とを用いて数値解析された結果）との誤差を算出する。相関パラメータ設定部としての機能を有する相関性算出部２５は、誤差算出部２４により算出された誤差に対して相関（Correlation）を持つ例えば渦度（ζ）などの所定のパラメータを指定し（相関パラメータを設定し）、誤差の算出対象となった流れ場全体の領域にわたって誤差とパラメータとの相関性（羽根車３３の径方向の位置に応じた相関性）を数値化する。
【００３０】
解析領域分割部２６は、図５に示すように、誤差算出部２４による誤差の算出対象となった流れ場全体の領域を、スワール（羽根車３３）の径方向に分割する。具体的には、解析領域分割部２６は、誤差の算出対象となった流れ場全体の領域を、相関性算出部２５により数値化された相関性の値が閾値を超える第１の解析領域（相関の大きい領域Ａ［ルート部］、領域Ｃ［チップ部］）と、当該数値化された相関性の値が閾値以下の第２の解析領域（相関の小さい領域Ｂ）と、に分割する。閾値記憶部２７は、解析領域分割部２６が領域の分割の際に参照する閾値（相関性の値に対応した閾値）を予め記憶している。
【００３１】
また、前述した数値解析部２２は、Cubic ｋ−εモデル又はＲＳＭなどの所定の乱流モデルと、この所定の乱流モデルに対応した乱流モデル係数（第１の乱流モデル係数）の値を補正した第２の乱流モデル係数と、を用いて、第１の解析領域（図５に示す領域Ａ、領域Ｃ）における乱流の特性を数値解析（数値シミュレーション）する。さらに、数値解析部２２は、所定の乱流モデルとこの所定の乱流モデルに対応した乱流モデル係数（第１の乱流モデル係数）とを用いて、第２の解析領域（図５に示す領域Ｂ）における乱流の特性を数値解析する。
【００３２】
ここで、上述した第２の乱流モデル係数は、次のようにして求められる。すなわち、乱流モデル係数変更部２８は、第１の乱流モデル係数の値を段階的に増減させ、この第１の乱流モデル係数を、複数の仮の乱流モデル係数に順次置き換える。一方、数値解析部２２は、上記した所定の乱流モデルと複数の仮の乱流モデル係数のうちの一つとを順次用いて、第１の解析領域における乱流の特性の数値解析を繰り返す（乱流モデル係数変化解析を行う）。さらに、乱流モデル係数確定部２９は、この繰り返し行われる数値解析の解析結果と、測定結果取得部２１により取得された測定結果と、の誤差が最も小さくなったときの前記仮の乱流モデル係数の値を、第２の乱流モデル係数として検出する。
【００３３】
次に、本実施形態に係る乱流のシミュレーション方法を、主に図６に示すフローチャートに基づいて説明する。まず、測定結果取得部２１は、スワールの生じる流れ場（図２に示すランナーコーン３３ｂの直下）で測定された乱流の特性（軸流速度など）の測定結果（Vefd）を取得する。一方、数値解析部２２は、所定の乱流モデル（例えば標準ｋ−εモデル）と、この乱流モデルに対応した第１の乱流モデル係数（例えば式４に例示にしたＣμ、σk、σε、Ｃ1ε、Ｃ2εなど）と、を少なくとも用いて数値解析を行い、スワールが生じる流れ場における乱流の特性（軸流速度など）の解析結果（Vcfd）を得る。
【００３４】
次に、誤差算出部２４は、図６に示すように、流速などの測定結果（Vefd）と数値解析して得られた流速などの解析結果（Vcfd）との誤差（Error）を算出する（Ｓ１）。さらに、相関性算出部２５は、この誤差に対して相関（Correlation）を持つ例えば渦度（ζ）などの所定のパラメータを指定し、誤差の算出対象となった流れ場全体の領域（ランナーコーン３３ｂの直下の領域）にわたって誤差とパラメータとの相関性を数値化する（Ｓ２）。
【００３５】
続いて、図６に示すように、数値化された相関性が閾値を超える場合（Ｓ３のＹｅｓ）、解析領域分割部２６は、図５に示すように、誤差の算出対象となった流れ場全体の領域を、相関性の値が閾値を超えるルート部やチップ部などの第１の解析領域（相関大の領域Ａ、領域Ｃ）と、数値化された相関性の値が閾値以下の第２の解析領域（相関小の領域Ｂ）と、に分割する（Ｓ４）。
【００３６】
次いで、数値解析部２２は、乱流モデル係数変更部２８及び乱流モデル係数確定部２９により第１の乱流モデル係数の値を補正した第２の乱流モデル係数と、乱流モデルと、を用いて（乱流モデル係数変化解析にて）、第１の解析領域（図５に示す領域Ａ、領域Ｃ）における乱流の特性を数値解析する（Ｓ５）。また、数値解析部２２は、第２の解析領域（図５に示す領域Ｂ）については、第１の乱流モデル係数と乱流モデルとを用いて数値解析を行う。
【００３７】
既述したように、本実施形態に係る乱流のシミュレーション方法によれば、スワールの生じる流れ場における乱流の特性の測定結果と初期の数値解析による解析結果とを効果的に融合することによって、大幅なコストの増加などを招くことなく、最終的な数値解析（数値シミュレーション）の精度を高めることができる。これにより、スワールを発生させ得る流体機械の性能を向上させることができる。
【００３８】
以上、本発明の実施形態を説明したが、この実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。この新規な実施形態は、その他の様々な形態で実施することが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。この実施形態やその変形例は、発明の範囲や要旨に含まれると共に、特許請求の範囲に記載された発明とその均等の範囲に含まれる。
【００３９】
例えば、本実施形態に係る乱流のシミュレーション方法や乱流シミュレーション装置２０は、カプラン水車の他、フランシス水車や、また、スワールを発生させ得るあらゆる流体機械に対して適用することが可能である。
【００４０】
なお、本実施形態の乱流シミュレーション装置２０が備えていた各構成要素（誤差算出部２４、解析領域分割部２６、乱流モデル係数変更部２８などを始めとする図１中に機能ブロックとして示した各構成要素）を、コンピュータに実行させることのできる乱流シミュレーションプログラムを構成してもよいし、さらには、この乱流シミュレーションプログラムを、光ディスクや半導体メモリなどの記憶媒体に格納して頒布することなども可能である。
【符号の説明】
【００４１】
２０…乱流シミュレーション装置、２１…測定結果取得部、２２…数値解析部、２４…誤差算出部、２５…相関性算出部、２６…解析領域分割部、２７…閾値記憶部、２８…乱流モデル係数変更部、２９…乱流モデル係数確定部。
【図１】
2017004103.tif 000007
【図２】
2017004103.tif 000008
【図３】
2017004103.tif 000009
【図４】
2017004103.tif 000010
【図５】
2017004103.tif 000011
【図６】
2017004103.tif 000012