宇宙航空研究開発機構特別資料

宇宙航空研究開発機構特別資料

第 45 回流体力学講演会／

航空宇宙数値シミュレーション技術シンポジウム 2013 論文集

JAXA Special Publication

Proceedings of 45th Fluid Dynamics Conference / Aerospace Numerical Simulation Symposium 2013

2014 年 3 月 March 2014

宇宙航空研究開発機構

Japan Aerospace Exploration Agency

開 催 日：平成 25 年 7 月 4 日（木）～ 5 日（金）

開催場所：タワーホール船堀

4 July ～ 5 July, 2013

Tower Hall Funabori

今年度の ANSS2013 は、例年通り日本航空宙学会空気力学部門委員会による流体力学講 演会との合同開催という形で、タワーホール船堀（東京都・江戸川区）にて開催いたしま した。本年は、ANSS側、学会側からの特別企画8テーマ及び一般講演について122件の発 表があり、5会場を用いて講演を進め、活発な議論を行うことができました。また特別講演 につきましては、JAXA 航空本部・中橋和博氏より「ペタフロップス計算機時代における CFD」、日建設計の小西厚夫氏より「東京スカイツリーの構造設計」、スタンフォード大学 Gianluca Iaccarino教授より「Progress in Quantification of Uncertainties in Fluid Flow Simulations」 との題目でご講演をいただくとともに、新たな試みとして、パネルディスカッション「航 空教育支援フォーラム：空力教育の共通プラットフォーム構築に向けて」を実施し、239名 の参加登録人数を数えることができました。

特別企画としては、EFD/CFD 融合技術、民間超音速機実現のための空力技術、非定常空 力と空力音響技術、先進流体計測技術、デトネーションエンジン、航空宇宙における HPC 利用技術、宇宙輸送及び推進系技術などに関するセッションを設け、関連分野の情報交流 の活性化を試みました。また一般発表では、低レイノルズ数翼型の性能・特性、プラズマ 関連の研究発表が多くみられ、数値シミュレーションのみならず、実験・計測との対応な どについて活発な議論がなされました。「航空教育支援フォーラム：空力教育の共通プラッ トフォーム構築に向けて」と題したパネルディスカッションでは、空力教育にCFDプログ ラムを共通プラットフォームとして持ち込み、CFD を空力教育の新たなツールとする進め 方に関する議論がなされました。私ども実行委員会の責務として、学会活動の一環として このような新しい技術をいち早く教育の現場へ浸透させるための議論の場を提供し続けて いくべきであるとの感を強くしました。

最後に、本シンポジウムの運営に当たり、日本航空宇宙学会空気力学部門委員長の佐宗章弘 名古屋大学教授をはじめ同部門委員の方々のご努力に謝意を表します。

平成25年12月吉日

航空宇宙数値シミュレーション技術シンポジウム（ANSS） 委員長

越岡 康弘

ANSS運営委員会委員

越岡康弘（委員長）、新城淳史（幹事）、相曽秀昭、青山剛史、池田友明、榎本俊治、齋藤 健一、佐藤茂、嶋英志、清水太郎、長谷川進、藤田直行、牧野好和、松尾裕一、松山新吾、

村上桂一、村山光宏、山根敬、吉田正廣

１．回転翼空力弾性解析コードの風車への適用...1 田辺安忠，杉浦正彦（宇宙航空研究開発機構.航空本部.機体システム研究グループ）

菅原瑛明（菱友システムズ）

２．構造格子 CFD コード UPACS-LES による四輪型簡易脚形状 RLG の空力音響予測の検証...7 田中健太郎（菱友システムズ），村山光宏，山本一臣，池田友明，榎本俊治（宇宙航空研究開発機構）

雨宮和久（エイ・エス・アイ総研）

３．Multipole.Analysis による超音速飛翔体の近傍場波形の改善...13 金森正史，橋本敦，青山剛史，牧野好和（宇宙航空研究開発機構）

石川敬掲（三向ソフトウェア開発），山本雅史（計算力学研究センター），飯村拓哉（菱友システムズ）

４．非線形 Tricomi 方程式解析を用いたフォーカスブームにおける低ブーム波形の効果推算...19 金森正史，橋本敦，青山剛史（宇宙航空研究開発機構），山本雅史（計算力学研究センター）

５．リード弁式吸気機構によるマイクロ波ロケットの推力性能改善...25 福成雅史，山口敏和，齋藤翔平，浅井健太，栗田哲志，小紫公也（東京大学），

小田靖久，梶原健，高橋幸司，坂本慶司（日本原子力研究開発機構）

６．ローテーティングデトネーションエンジンの数値解析：厚み方向の影響とスケール効果...29 図斉健太，桜澤歩，朝原誠，山田英助，林光一（青山学院大学大学院.理工学研究科.理工学専攻.機械創造コース）

坪井伸幸（九州工業大学大学院工学研究院）

７．高次精度非構造格子法の比較に関する研究...35 澤木悠太（東北大），芳賀臣紀，保江かな子（JAXA），澤田恵介（東北大）

８．高迎角剥離流の非定常解析に向けて...41 橋本敦，石田崇，石向桂一，青山剛史（宇宙航空研究開発機構）

９．極超音速希薄風洞流れ場の粒子計算解析...47 小澤宇志，鈴木俊之，藤田和央（宇宙航空研究開発機構）

10．磁気シールドに対する印加磁場配位の効果について...53 永田靖典（JAXA/ISAS），里深優，渡辺理成（早稲田大学），山田和彦，安部隆士（JAXA/ISAS）

11．回転翼のホバリング飛行解析への CFD./ 規定後流モデルハイブリッド手法の適用...59 武田茂（首都大学東京），田辺安忠，杉浦正彦，張替正敏（宇宙航空研究開発機構）

菅原瑛明（菱友システムズ）

12．ヘリコプタの BVI 騒音予測に向けた CFD と規定後流モデルのハイブリッド手法の改良...65 杉浦正彦，田辺安忠（宇宙航空研究開発機構），菅原瑛明（菱友システムズ）

頓所和之（菱友システムズ）

14．遷音速流における翼後流 PIV 計測による圧力推定...77 松島紀佐，泉知宏（富山大学・工），加藤裕之（JAXA）

15．縮約モデルと粒子フィルタを用いたリアルタイムデータ同化計算...83 菊地亮太，三坂孝志，大林茂（東北大学.流体科学研究所）

16．超音速領域における PIV 計測データの補正方法に関する研究 –.MTV データとの比較 –...89 三井克仁，半田太郎（九州大学），中野葵（川崎重工業），小池俊輔（宇宙航空研究開発機構）

17．吸収飽和誘起蛍光法による光学的に厚いプラズマ流中の温度分布計測...95 伊藤彦，金子剛，小紫公也（東京大学大学院.新領域創成科学研究科），野村哲史（宇宙航空研究開発機構）

SCHÖNHERR.Tony，小泉宏之（東京大学大学院.工学系研究科）

18．実在気体気流条件での空力係数の計測...99 丹野英幸，佐藤和雄，小室智幸，伊藤勝宏（宇宙航空研究開発機構.角田宇宙センター）

藤田和央（宇宙航空研究開発機構.調布航空宇宙センター）

19．観測ロケット実験における柔軟構造飛翔体の空力解析について... 103 高橋裕介，河東顯（北海道大学.大学院工学研究院 / 工学院）

山田和彦，安部隆士（宇宙航空研究開発機構.宇宙科学研究所）

鈴木宏二郎（東京大学.大学院新領域創成科学研究科）

20．大迎角細長物体の横力制御における DBD プラズマアクチュエータ設置位置の検討... 109 佐藤雅幸，西田浩之，松原暁良（東京農工大学），野々村拓（JAXA）

21．OpenFOAM を用いた NACA0012 翼型まわりの準二次元解析... 115 中谷淳，村澤杏樹（岐阜工業高等専門学校）

22．低 Re 数における Ishii 翼型まわりの流れ場... 121 大竹智久，村松旦典，本橋龍郎（日本大学理工学部），互井梨絵（JAXA），神田翔（ケーヒン）

23．低 Re 数領域における NACA0012 翼面上の圧力分布... 127 山口裕太（日本大学大学院），大竹智久，村松旦典（日本大学理工学部）

24．AMR 法による高揚力装置流れの高解像度数値解析... 131 松尾裕一（JAXA），富塚孝之，中森一郎（アドバンスソフト）

25．マルチブロック構造格子における NURBS.Volume を利用した自動細分化ツールの開発... 139 松村洋祐（みずほ情報総研），堤誠司，高木亮治，山本一臣（宇宙航空研究開発機構）

伊藤浩之，竹川国之（菱友システムズ）

27．スクラムジェットエンジン燃料最適化分布に向けた検討－質量流率との対比... 151 佐藤茂（宇宙航空研究開発機構角田），渡邉孝宏（日立東日本ソリューションズ）

福井正明（スペースサービス），宗像利彦（日立東日本ソリューションズ）

回転翼空力弾性解析コードの風車への適用

田辺安忠、杉浦正彦

宇宙航空研究開発機構 航空本部 機体システム研究グループ 菅原瑛明

菱友システムズ

Application of a Rotorcraft Aeroelastic Analysis Code to Wind Turbines

by

Yasutada Tanabe, Masahiko Sugiura（JAXA）

and Hideaki Sugawara (Ryoyu Systems, Co., Ltd) ABSTRACT

Application of the CFD/CSD coupling analysis code rFlow3D which was originally developed for rotorcraft to the wind turbines is described in this paper. NREL Phase VI wind turbine is selected as the test case to validate the accuracy of the analyses. It is found that even with a relative coarse resolution blade grid, the Euler solver can predict the performance of the wind turbine quite satisfactory when the wind turbine is operated in attached flow or deep stall conditions, but the discrepancies in the incipient stall flow region are remarkable where the highest power is generated. With a Navier-Stokes solver, the prediction accuracy is improved with the refinement of the blade grid and the highest power can be predicted satisfactorily.

１．はじめに

クリーンエネルギの代表として、風力発電が注目されて おり、近年設置台数が急激に増加している。設置環境も陸 上から洋上へと多様化が進んでいる。風力発電の効率化を 目指して、風車の大型化に伴い、ブレードの弾性変形が無 視できないレベルにある。日本のような風環境が急変しや すい中、風車の破損も目立っており、荷重の状況を正しく 予測し、設計に反映することが求められている。

JAXA においては、ヘリコプタの空力弾性騒音を解析す るために、CFD/CSD連成解析技術に基づく統合解析ツール

rFlow3Dを開発してきた[1-3]。主にヘリコプタの騒音低減

を目指して、種々のアクティブ・デバイスの作動によるロ ータ・ブレードの弾性変形、空力騒音変化の評価に用いら れてきており、実験結果による検証でも良好な結果が得ら れている。この解析ツールの同じ回転翼でもある水平軸風 車への適用は、自然な解析技術の応用と思われる。しかし ながら、ヘリコプタは基本的に流れの剥離が起こらない範 囲内で運航されることに対し、風車は自然の風の風速変化 から完全に流れが剥離している範囲まで考える必要があり、

流体力学的にはより複雑である。そのため、CFD手法につ いても、ヘリコプタでは Euler 方程式に基づくソルバーで もほとんどの場合かなり精度の高い結果が得られるが、風 車のケースでは、流れ場の剥離のシミュレーション精度が そのまま結果に大きく反映されるため、Navier-Stokes 方程 式に基づくソルバーであることが前提であるように思われ る。今回は rFlow3Dの粘性計算を基本に、NRELの Phase VIの風車[4]を検証対象として、ソルバーの種類や境界条件、

計算格子の影響などについて、検証を行ったので、その結 果について報告する。

２．計算手法

ヘリコプタや風車のロータ・ブレードは非常に細長い形 状をしており、荷重に応じて弾性変形が生じやすい。近年 海上での風力発電が増加傾向にあり、６MW 級の風車では 直径が 120m以上にもなる。アスペクト比がヘリコプタよ りも大きくなる場合が多く、弾性変形の影響が在来のロー タ・ブレードよりも更に重要になっている。ブレードの弾 性変形解析を含む回転翼の空力性能や騒音を予測するため には、図１に示す流れで計算を進める必要がある。

rMode rFlow3D /

JANUS rNoise

Blade structural properties Rotating speed Root constraint conditions

Natural frequencies Natural mode shapes (represented with

polynomials)

Flow conditions Rotor blade controls

Trim conditions

Surface pressure on rotating blades

Observer conditions

Sound pressure time history at observation

points

図１ 回転翼統合解析ツールチェーン

JAXAで開発されている回転翼向けの統合解析ツール[1- 3]は、3 つのコードで構成されている：まずはロータ・ブ レードの固有周波数と固有モードを計算する rModeという コ ー ド で あ る 。 ブ レ ー ド の 弾 性 振 動 の 支 配 方 程 式 は Houbolt & Brooksの方程式［5］に基づいており、非均一の ねじれのあるビームを仮定している。ブレードの構造物性 値とロータの回転数に基づいて、ブレードのフラップ、ラ グ、ねじりの 3方向の連成固有振動周波数と振動モードの

形を Holzer-Myklestad 法[6]を用いて計算する。ここでの出

力は次の段階の CFD/CSD 連成解析コードの入力となると 同時に、ロータ回転数と固有周波数の関係を示す Fan-plot はロータブレードの回転時の安定解析にも使われる。

解析 の中心である CFD/CSD/TRIM の 連成解析コード rFlow3Dの計算手法は移動重合格子法である[7]。ロータブ レードの枚数分の内部格子をブレード周りに形成し、ブレ ードの回転や弾性変形に合わせて移動・変形をし続ける。

また、ヘリコプタの胴体は複雑な形状をしている場合が多 く、胴体周りの格子は非構造格子も採用できるように拡張 されている。胴体周りが非構造格子の場合は TAS-Code を ベースにしたソルバーを使用しており、JAXA と東北大学 との共同開発で、JANUSというコード名で呼ばれる[8]。背 景格子は直交格子を 2 層まで用いることができ、ブレード の回転領域や後流領域では密な格子を採用して、翼端渦を 精度よく捉えられるようにしている。各格子間のデータの 受け渡しは 2次精度の Tri-Linear 補間法で行っている。内 部格子では移動格子に対応した完全非定常 Euler/Navier- Stokes方程式を 4次精度の SLAU+FCMTスキーム[9,10]で 離散化し，ヘリコプタのような低速から遷音速領域が共存

CFD で求めた空気力に基づいて、構造解析との弱連成手法 で計算される。さらに目標となるロータ推力やモーメント と一致するように、ブレードの制御入力を変更し、トリム を取りながら、現実の飛行条件を再現できるようにしてい る。

また、回転翼からの騒音計算は FW-H 方程式に基づく

Farassat 定式１[11]に従い一周分のブレード表面圧力から積

分で求めており、rNoise という後処理ツール[2]で行うよう になっている。

以上の解析ツールチェーンはアクティブ制御によるヘリ コプタの騒音低減の実験結果と比較検証[3]をしてきており、

良好な一致を得ている。しかしながら、ヘリコプタの場合 は、ブレード周りの流れ場は付着している場合がほとんど で、非粘性の Euler 解でも良好な結果が期待できるが、風 車への適用に当たっては、ブレード周りの流れ場が非常に 幅広い迎え角範囲を有する場合が多く、粘性効果が無視で きないレベルにある。そのため、CFDの支配方程式は以下

のようなThin-Layer形式[12]の粘性項を追加した。

S G

F E

q ˆ

Re ˆ 1 ˆ

ˆ

ˆ

τ

+ ∂ + ∂ + ∂ = ∂

∂

粘性項の離散化については、FVM法で用いるメトリクスと

差分法のJacobianとの対応関係を用いて、中心差分となる

ように定式化した。また、今回の解析対象の風車のレイノ ルズ数が低く、乱流モデルを採用しなかった。

さらに、固定壁面上の圧力の与え方について、内側のセ ルの値からの外挿の精度について、ゼロ次外挿は壁面に隣 接するセルの圧力をそのまま使い、

1 , , , ,jw i j

i

p

p =

で、2次外挿はさらに一個上のセルの値も入れて、

6 / ) 7

(

,

,jw ij i j

i

p p

p = −

である。壁面上の運動量を満足する圧力解[12]は差分法的 に書くと、以下のようになる：

( ) ( )

( )

[ ]

( )

(

)

ξ ξ ξ

τ τ

τ τ

ζ

η ξ

ζ ζ ζ ρ

ζ ζ ζ ρ

ζ ζ

ζ ζ ρ

ζ ζ ζ

ζ η ζ η ζ η ζ

ξ ζ ξ ζ ξ

w v u V

w v u U

w v

u p

p p

z y x

z y x

z y

x t

z y x

z z y y x x z

z y y x x

+ +

−

+ +

−

∂ +

∂ +

∂ +

∂ +

+ +

−

=

+ + + +

+

2 2 2

（４）

式（4）はブレードの周方向については循環型3列対角マト リクス、スパン方向については通常の3列対角マトリクス を解くことで、壁面上の圧力を得ることができる。

３．水平軸風車検証計算例

NRELが NASA-Amesの 40ftx80ft の大型風洞を使って

Phase VIという水平軸風車の風洞試験を実施しており、風

洞試験と同時に、各種解析手法の予測結果の比較が行われ た[13]。この実験データはその後も広く解析コードの検証 例として利用されてきた[14-17]。風車の半径は 5.029mで、

試験回転数は 72RPM である。このとき、ブレード先端の 回転速度マッハ数は 0.11 である。ヘリコプタの場合は通 常 0.6ぐらいで、圧縮性の解析コードをこのような低いマ ッハ数の流れ場に適用できるか、SLAU 全速度スキームを 採用しているとは言え、十分な確認が必要である。

図2にrFlow3Dで使用される重合格子系を示す。今回の

解析では、外部背景格子は3方向に101点の等間隔直交格 子で、風車の回転中心から、3方向±20m の範囲をカバー

界面はゼロ次外挿の流出境界とした。

内部背景格子はブレードの後流を捉えるために細かく分 割しており、各辺 14m の立方体を 151 点で分割している。

V

Outer Background Grid

Inner Background Grid

Blade Grid

図2 風車の重合計算格子

図3にNREL Phase VIのブレードのモデル形状を示す。

回転中心からのブレードの半径は 5.029mで、翼型は s809、

ルート部から25%Rの位置から翼型となっており、そこで

の翼弦は0.737mで、先端での翼弦長は0.356mで、線型的

なテーパ形状となっている。また、この部分から先端に向 けて、捩じりあげ角が約 20度である。ブレードの形状と 風洞試験の詳細は文献[12]を参照されたい。

図3 NREL Phase VIのブレード形状

表1 ブレード計算格子

Blade Grid Points (IxJxK) ∆hmin [m]

EU 81x81x21 5e-3

NS2 81x81x61 5e-5

NS3I 161x81x61 5e-5

NS3J 81x161x61 5e-5

NS4 161x161x61 5e-5

表1に解像度による依存性を調べるために使用した5種 類ブレード格子の格子点やブレード表面に接する層の代表 高さを示す。EU というオイラーソルバー用の格子は図 4 に示すように、比較的に粗くブレード周りを覆うようにし ており、効率的な CFD 解析が可能である。一方、粘性計 算では、特に表面に隣接するセルの高さが粘性境界層を捉 えるのに重要なため、今回は代表高さが 0.05mmとし、こ れは代表レイノルズ数において、y⁺≒5に相当し、層流境 界層を捉えるには十分である。半径方向と周方向の分解能 の影響をみるため、NS2から NS4 まで 4 段階の格子を用 いた。NS４の格子の様子も図4に示す。

EU_grid_81x81x21

NS4_grid_161x161x61

図4 ブレード計算格子の様子

オイラー計算においては、壁面境界条件の影響を調べる ため、表2に示す3種類を行った。粘性計算については、

境界条件を2次外挿による壁面圧力のみにし、格子の影響 について調べた。

風車の運転条件については、実験ではブレードにコーニ ング角を付けたり、ヨー角のある場合の計測も含めている が、今回の検証計算では、コーニング角無し、ヨー角も無 い、ピッチ角を 3度に固定した、風速のみ変化させたケー スと、風速を 15m/s に固定し、ピッチ角を変化させたケー スの 2ケースを選んだ。この2ケースはこれまでに多くの 解析結果との比較に使われてきた[13-17]。

表2 計算条件 Abbr. Solver Blade

Grid Boundary Condition EU00 Euler EU Zeroth Interpolation

EU01 Euler EU 2nd Interpolation

EU10 Euler EU Wall momentum eqn

NS2TL01 TL NS NS2 2nd Interpolation NS3ITL01 TLNS NS3I 2nd Interpolation NS3JTL01 TLNS NS3J 2nd Interpolation NS4TL01 TLNS NS4 2nd Interpolation

風車のピッチを3度に固定し、風速を5，7，10，15，20，

25m/sの 6つに変化させたときの流れ場の様子を図 5に示

す。5m/s と 7m/sの流れ場は非常に滑らかで、後流渦も若 干広がりながら、はっきりとしたらせん状渦を形成してい る。10m/s以上の風速ではブレードまわりの流れが剥離し て、後流がかなり乱れているのが分かる。計算条件によっ

て、特に 10m/s以降では若干の違いが見られるが、定性的

には同じである。

10m/s, NS4TL01 5m/s, NS4TL01 7m/s, NS4TL01

25m/s, NS4TL01 20m/s, NS4TL01

15m/s, NS4TL01

図5 風速変化による風車の後流の変化

オイラーソルバーの解析能力をみるため、敢えて流れが 大きく剥離している条件でもそのまま計算した。図６の推 力変化をみると、若干計算値が大きくなっているが、風速

による推力のほぼ直線的な変化を正しく捉えており、また、

境界条件による差がほとんど見られない。しかしながら、

図 7に示す回転軸周りの空力トルクから計算した発電量は、

壁面圧力を0次または2次外挿で決める境界条件のEU00 とEU01は7m/sまでは実験値とも一致しているが、

7~15m/sの間は発電量を顕著に低く予測していた。偶然か

もしれないが、15m/s以降の流れが大きく剥離する条件で は、また実験値とよく一致した。普段オイラー計算の非粘 性解は流れ場の剥離をかなり遅く予測する傾向にあること が知られており、壁面圧力を単純外挿で決める境界条件で は、却って早期剥離が起こったものと考えられる。一方、

壁面上の運動量方程式を満足するEU10の境界条件では、

風速15m/sまでは非常に実験値と一致する結果が得られ、

理想的なオイラー解であると思われる。それ以降は実験値 よりも大きく発電量を見積もった。

流れ場が付着している状態から大きく剥離する条件まで 全領域にわたって風車の性能を正確に予測できるか、粘性 計算を実施した。5m/sにおける計算格子を変更した場合の 発電量の予測値を図8に示す。総じて実験値より若干大き くなっているが、計算値同士の差がほとんどない。なお、

粘性計算においても、トルクの計算は表面圧力の積分のみ を考慮し、摩擦成分による寄与を無視している。これは実 験値自体がブレード表面の圧力分布から計算されたもので あるからである。なお、参考文献［17］にあるように、こ の場合、摩擦成分の寄与が圧力値成分と比較すると無視で きるレベルである。

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

5 10 15 20 25

Thrust [N]

Wind Speed [m/s]

NREL_VI @ 3deg

Exp. (NREL) EU00 (JAXA) EU01 (JAXA) EU10 (JAXA)

図6 境界条件による推力の変化

0 5 10 15 20

5 10 15 20 25

Power [kW]

Wind Speed [m/s]

NREL_VI @3deg EXP (NREL)

EU00 (JAXA) EU01 (JAXA) EU10 (JAXA)

図7 境界条件による発電量の変化

0 2 4 6 8 10 12

EU01 NS2TL01 NS3ITL01 NS3JTL01 NS4TL01 NREL_exp

Power [kW]

NREL_VI @ 5m/s, 3deg

図8 風速5m/sにおける各種計算手法の予測結果 流れ場が初期剥離状態にある10m/sにおける発電量の予 測値を図9に示す。このときの流れ場は計算格子の分解能 に非常に敏感であり、格子が密になるとともに、計算値が 実験値に近づいた。なお、格子点数ではNS3IとNS3Jは同 じであるが、周方向を細分したNS3Jのほうが実験値に近 かった。

0 2 4 6 8 10 12 14

EU01 NS2TL01 NS3ITL01 NS3JTL01 NS4TL01 NREL_exp

Power [kW]

NREL_VI @ 10m/s, 3deg

図9 風速10m/sにおける各種計算手法の予測結果

大剥離状態にある25m/sにおける発電量の計算手法によ る変化を図10に示すが、粘性計算値が総じて実験とオイラ ー計算値より低かった。さらに格子が密になる改善効果が 見られなかった。おそらくこの速度では、もはや流れが層 流である仮定が成り立っておらず、何らかの乱流モデルの 採用が必要だったかもしれない。乱流モデルの採用による 予測については今後の課題としたい。

0 2 4 6 8 10 12 14

EU01 NS2TL01 NS3ITL01 NS3JTL01 NS4TL01 NREL_exp

Power [kW]

NREL_VI @ 25m/s, 3deg

図10 風速25m/sにおける各種計算手法の予測結果

推力について、NS4の粘性解と同じ境界条件のオイラー 解と計測値との比較を図11に示す。ほとんどの場合、粘性 解のほうが若干推力を大きく予測する傾向にあった。

風車の発電量については、図12に実験値との比較を示す が、剥離が顕著な領域では若干低く予測しているが、性能 に重要なピーク発電の予測値が計測値とよく一致しており、

実用上、十分な精度の性能予測ができると考えている。

きく変化させたときの推力の変化を図13に示す。EU00条 件の計算結果が実験値とよく一致していた。

図14に発電量の変化を示す。全体的な傾向が概ね一致し ているが、ピッチ角 15 度時にピーク発電量があり、EU00 の予測値が 25%も低く見積もっていた。今後は粘性解によ る解析も進めていく予定である。

0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000

5 10 15 20 25

Thrust [N]

Wind Speed [m/s]

NREL_VI @ 3deg

Exp. (NREL) EU01 (JAXA) NS4TL01(JAXA)

図11 風車推力の粘性解

0 5 10 15 20

5 10 15 20 25

Power [kW]

Wind Speed [m/s]

NREL_VI @3deg EXP (NREL)

EU01 (JAXA) NS4TL01 (JAXA)

図12 風車発電量の粘性解

-1500 -10001000150020002500300035004000-5005000

-20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45

Thrust [N]

Pitch [deg]

NREL_VI @15m/s

Exp (NREL) EU00 (JAXA)

図13 風車のピッチ角による推力変化

-25 -20 -15 -10 -50 5 10 15 20 25

-20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45

Power [kW]

Pitch [deg]

NREL_VI @15m/s

EXP (NREL) EU00 (JAXA)

図14 風車のピッチ角による発電量変化

なお、以上の試験計算は風車のブレードを完全に剛とし て、弾性変形しないと仮定している。今後は弾性変形も含 めて解析を進め、風車性能への影響を把握していきたい。

４．まとめ

回転翼航空機用に開発した統合解析ツールrFlow3Dを水 平軸風車に適用し、幅広い試験条件下での空力荷重と発電 量のテスト計算を行った。

NRELのPhase VIの風洞試験データに基づいて、特にピ

ッチ角3度固定で、風速を変化させた試験例に対して、境 界条件の影響や、格子の解像度の影響について詳細に調べ た。

風車の受ける空力荷重となる推力は、格子の解像度の影 響をあまり受けず、どのケースも概ね実験値と良い一致を 示した。

初期剥離状態（風速 10m/s）においては、壁面上の圧力 を単純外挿した Euler解には、流れ場の早期剥離現象が見 られ、最大発電量を低く予測する傾向があった。壁面上の 圧力を運動量保存方程式を満足するように与えた場合は、

Deep Stallとなる風速15m/sまではEuler法でもかなり精度 よく性能を予測できが、それ以降の速度領域では発電量を かなり大きく予測した。

計算格子を密にした粘性計算の結果は全体的に実験値と 一致したが、Deep Stall域では若干発電量を低く見積もる 傾向にあった。

これまでの試験計算の結果、回転翼用の統合解析ツール の CFD ソルバーは風車の荷重や性能予測に十分使えるこ とが分かった。今後はさらに弾性変形や騒音の予測などの 計算機能についても、検証を進めていきたい。また、これ までは全領域を CFD で解いたが、計算の効率化と精度の 両立を図るため、風車の後流を渦モデルで表現し、ブレー ドの周りのみをCFDで解くHybrid Methodの開発を進めて おり、風車の計算例での検証を進めていく予定である。

参考文献

1) 田辺：回転翼の空力弾性計算手法の考察、第 50 回飛行 機シンポジウム、新潟市、2012年11月5~7日。

2) 田辺：回転翼機の空力騒音予測ツールの構築について、

第 44 回流体力学講演会／航空宇宙数値シミュレーショ ン技術シンポジウム 2012、富山市、2012年7月5~6日。

3) Tanabe, Y., Saito, S. and Sugawara, H., "Construction and Validation of an Analysis Tool Chain for Rotorcraft Active Noise Reduction," 38th European Rotorcraft Forum, Amsterdam, NL, September 4-6, 2012.

4) Hand, M., Simms, D., Fingersh, L. J., Jager, D., Larwood, S., Cotrell, J., and Schreck, S., "Unsteady Aerodynamics Experiment Phase VI: Wind Tunnel Test Configurations and

5) Houbolt, J.C. and Brooks, G.W., "Differential Equations of Motion for Combined Flapwise Bending, Chordwise Bending, and Torsion of Twisted Nonuniform Rotor Blades," NASA TN-3905, 1957.

6) Isakson, G. and Eisley, G.J., "Natural Frequencies in coupled bending and torsion of twisted rotating and non-rotating blades," NASA CR-65, 1964.

7) Tanabe, Y. and Saito, S., “Significance of All-Speed Scheme in Application to Rotorcraft CFD Simulations,” The 3rd International Basic Research Conference on Rotorcraft Technology, Nanjing, China. October, 2009.

8) Tanabe, Y., Saito, S., Takayama, O., Sasaki, D. and Nakahashi, K., A New Hybrid Method of Overlapping Structured Grids Combined with Unstructured Fuselage Grids for Rotorcraft Analysis, 36th European Rotorcraft Forum, Paris, France, September 9-11, 2010.

9) Shima, E., and Kitamura, K., “On New Simple Low- Dissipation Scheme of AUSM-Family for All Speeds,” 47th AIAA Aerospace Sciences Meeting, Orlando, FA, January 5- 8 2009, AIAA Paper 2009-136.

10)Yamamoto, S. & Daiguji, H., “Higher- Order- Accurate Upwind Schemes for Solving the Compressible Euler and Navier-Stokes Equations,” Computers & Fluids, Vol.22, No.2/3, pp.259-270, 1993.

11) Farassat, F., "Derivation of Formulation 1 and 1A of Farassat," NASA TM-2007-214853, 2007.

12) Pulliam, T. H. and Steger, J. L., “Implicit Finite-Difference Simulations of Three-Dimensional Compressible Flow”, AIAA J., Vol. 18, No. 2, Feb. 1980, pp. 159-167.

13) Simms, D., Schreck, S., Hand, M., Fingersh, L. J., "NREL Unsteady Aerodynamics Experiment in the NASA-Ames Wind Tunnel: A Comparison of Prediction to Measurements,"

NREL/TP-500-29494, 2001.

14) Duque, E. P. N., Burklund, M. D., and Johnson, W., "Navier- Stokes and Comprehensive Analysis Performance Predictions of the NREL Phase VI Experiment," AIAA-2003-355, January 6-9, 2003.

15) Schmitz, S., and Chattot, J. J., "Application of a 'Parallelized Coupled Navier-Stokes/Vortex-Panel Solver' to the NREL Phase VI Rotor,", AIAA-2005-593, January 10-13, 2005.

16) Yu, D. O., Kwon, H. I., and Kwon, O. J., "Performance Enhancement of HAWT Rotor Blades by Aerodynamic Shape Optimization," AIAA-2012-1292, January 9-12, 2012.

17) Li, Y., Paik, K. J., Xing, T., and Carrica, P. M., "Dynamic Overset CFD Simulations of Wind Turbine Aerodynamics,"

Renewable Energy 37 (2012), pp.285-298.

構造格子 CFD コード UPACS-LES による四輪型簡易脚形状 RLG の 空力音響予測の検証

田中健太郎¹，村山光宏²，山本一臣²，雨宮和久³，池田友明²，榎本俊治²

1菱友システムズ，²宇宙航空研究開発機構，³エイ・エス・アイ総研

Validation Study of Aeroacoustics Analysis of Rudimentary Landing Gear Using Structured Grid CFD-code UPACS-LES

Kentaro Tanaka, Mitsuhiro Murayama, Kazuomi Yamamoto, Kazuhisa Amemiya, Tomoaki Ikeda and Shunji Enomoto by

ABSTRACT

In this paper, a validation study of aeroacoustics analysis of Rudimentary Landing Gear is conducted using a multi-block structured grid CFD-code, UPACS-LES, developed in JAXA. Sensitivities of the numerical schemes to solve the convective term and the order of the scheme to the flowfield and near-field unsteady pressure are investigated. Improved effects on the result to capture finer vortical structures and to predict far-field noise prediction are shown by the 3rd-order low dissipation SLAU scheme for low-Mach number flow as well as the 6th-order Compact scheme, compared with the 3rd-order Roe scheme. Influence of sound reflection on the floor on the far-field noise prediction is evaluated by Ffowcs Williams and Hawkings code using permeable surface data set in space. In addition, sensitivities of the placement of wake plane of the permeable surface on the far-field noise prediction are investigated with a method to suppress spurious noise due to vortex passing through the wake plane.

１．はじめに

近年，航空機のエンジン騒音の大幅な低減に伴い，特に エンジン出力を抑えた着陸時には高揚力装置や脚を主要な 騒音源とする機体空力騒音がエンジン騒音を上回る例も見 られる⁽¹⁾．そのため今後の旅客機開発においては将来更に 厳しくなると予測される騒音規制値に対応するため，機体 空力騒音の予測技術の精度向上と低騒音化技術の確立は重 要な課題である．

脚騒音は，数多くの単純な断面形をした構造部材の後流 や乱流せん断層，干渉によって構造物周辺で発生する広域 帯騒音が主である⁽²⁾．最近では非定常流体音響解析による 騒音予測が実施される様になってきているが，形状や現象 の複雑性からその結果の信頼性について十分な検証が必要 であり，機体空力騒音の予測技術向上を目指した国際ワー ク シ ョ ッ プ BANC (Workshop on Benchmark problems for Airframe Noise Computations)-I，-II⁽³⁾が 2010年，2012年に 開催されるなど，国際的な取り組みも始められている．

JAXAでは構造格子CFDコードUPACS-LES^(4,5)を用いて，

BANC での脚騒音に関する課題の中で最も単純化された直 列 2円柱問題に対して，格子解像度や乱流モデルの影響等 を 調 べ て き た⁽⁶⁾． 本 研 究 で は よ り 複 雑 な 問 題 と し て ， BANC の課題の中から四輪型脚の簡易形状である RLG (Rudimentary Landing Gear）に対して空力音響予測を行う．

低マッハ数流れとなる着陸時条件の解析に圧縮性CFDコ ードを用いる場合，過大な数値散逸による解像度の低下や 移流速度と音速の大きな差によるスティフネスが問題とな るが，嶋らはパラメータ調整なしに低マッハ数流れにおい ても過大な数値散逸を低減する SLAU スキーム⁽⁷⁾を提案し ている．本稿では3次精度 RoeとSLAUスキームおよび6 次精度 Compactスキームを用いた Delayed Detached Eddy Simulation (DDES)解析を行い，対流項スキームの違いによ る流れ場への影響について報告する．

遠方場騒音を分離解法にて予測する際にFfowcs Williams- Hawkings（FW-H）法⁽⁸⁾が多く用いられる^(3,9)．物体から発 生した主要な音源を取り囲む閉曲面（音響面）を定義し，

面上での時間履歴の値を用いて遠方場圧力変動を予測する．

FW-H 法では通例，計算コストの観点から体積分項を省略 しているが，渦が音響面を通過する際に非物理的な音（偽

音）が発生する原因となるため，音響面の取り方には注意 を要する．本稿では FW-H 法における床面での音波の反射 の影響を調べると共に，音響面の取り方による遠方場予測 値への影響と，複数の音響面の平均を取る事で偽音の影響 を緩和させるShurらの手法⁽¹⁰⁾の適用効果ついても報告する．

２．解析対象

航空機の四輪型脚の簡易形状であるRLGは，BANCの課 題の一つとされた形状である．四つの車輪を繋ぐ車軸・支 柱などは全て単純な角柱で構成されており，風洞試験と CFDとの間で境界層遷移の影響差も小さく，基礎的な検証 問題として提供されている．本研究ではインドの NAL に て実施された風洞試験^(11,12)の模型形状を用い，車輪直径

D=0.4064mとした．この模型形状は設計形状に対してわず

かに形状差を有するが，差分は模型スケールで 2mm 程度 である．

３．解析手法

空力音響予測の解析手順として，まず音源となる流れ場 解析を行い，そこで得られた近傍場データを用いて遠方場 騒音を予測する分離解法を用いる．

３．１．流れ場解析と計算格子

流れ場の解析には，JAXA 航空本部で開発されている UPACS-LES^(4-6)を用いる．ベースとなる UPACS⁽¹³⁾はマルチ ブロック構造格子に対応したセル中心型有限体積法による 3次元圧縮性 Navier-Stokes方程式ソルバーで，MPIを用い て並列化されている．UPACS-LES は非定常解析に特化し たバージョンで，LES，DESを扱える様に拡張されている．

乱流モデルとして Spalart-Allmaras 1 方程式モデル(SA-

noft2)を用いた DDES 解析を行う．ただし生成項では渦度

の代わりに歪み速度を用いる．対流項スキームの違いによ る結果の比較を行うため，3次精度 Roeスキームと，数値 粘性の小さい 3次精度 SLAUスキーム⁽⁷⁾，高次精度スキー

ムである 6次精度 Compactスキームの 3種類を用いる．

Compactスキームでは境界層内での計算の安定化のため，

物体近傍領域に Roeスキームを用いるハイブリッド法（以

降 Compact/Roeスキームと呼ぶ）とする．ここで Compact

と の切り替えには簡易的に物体表面からの距離 を閾

値とし，本研究ではd/D=0.01を用いる．Compactスキーム の領域には振動を抑えるため最大14次精度の空間フィルタ ーを用いる．粘性項には 2次精度中心差分を用いる．時間 積分には MFGS陰解法を用い，時間ステップ毎にニュート ン法による内部反復を5回行う事で時間精度を維持する．

図1に示すマルチブロック構造格子は Gridgenを用いて 生成している．格子解像度はSpalartら⁽¹⁴⁾の計算に用いられ た重合型構造格子と同程度となる様にしている．ただし粘 性壁での最小格子幅はやや粗く∆y/D=3.36×10^-5で，これは y⁺=1.4程度に相当する．格子点数は約3370万点，並列計算 時の負荷バランス平滑化のため 592 ブロックに分割してい る．外部境界は音波を減衰させるバッファ域を設け，100D で一様流条件を与える．床面はすべり壁面を想定し対称境 界条件を与える．ブロック間接続のデータ転送幅（以降袖 幅と呼ぶ）は，RoeおよびSLAUスキームでは2セル，高

次精度の Compact/Roeスキームでは 8セルを用いる．また

Compact/Roeで2セルとした時の影響も調査する．

(a) RLG表面格子とワイヤーフレーム

(b) 空間断面格子（z/D=0） 図1 マルチブロック構造格子

３．２．解析条件

解析条件は，一様流流速 U=40m/s（マッハ数 0.115），

車輪直径基準のレイノルズ数は 1×10⁶，一様流静温度は

288.15Kとする．迎角および横滑り角は共に 0°である．空

力係数を算出する際の基準面積はD²を用いる．

流 れ 場 解 析 の 手 順 と し て ， ま ず 局 所 時 間 刻 み に よ る RANS解析にて定常流れ場を求める．続いて時間刻み幅∆t を固定した DDES解析に切り替え，非定常流れ場を発達さ せるため，一様流流速Uと車輪直径Dとで規格化された時 間TU/D=10程度（約0.1secに相当）計算を進める．境界層 付近を除きCFL<1となる様に，時間刻み幅∆tU/D=3.45×10^-4 とする．その後，統計量データ（平均流れおよび変動成 分）と音響解析用データを取得するための計算へ移行する．

統計量および音響解析用データの取得時間は，Roe，

SLAU，Compact/RoeスキームでそれぞれTU/D=44，28，16 程度である．データの取得は 20step毎に行い，サンプリン グ周波数は fD/U(=St)=140 程度（約 14kHz に相当）である．

計算にはノード間通信がInfiniBand化されたPCクラスタ

（Intel Xeon X5690×2CPU×16ノード，全192コア）の 128 コアを用いる．統計量および音響解析用データの取得には，

TU/D=14につきRoeとSLAUスキームでは約4日，袖幅に 8セルを用いた Compact/Roeスキームでは約 5日を要する．

３．３．遠方場音響解析

遠方場騒音の予測にはFW-H法を用いる．RLGおよびそ の後流など主要な音源となる非定常渦が生成する領域を取 り囲む閉曲面（音響面）を定義し，その面上でサンプリン グした時系列の値を用いて遠方場圧力変動を予測する．図 2には本解析に用いる音響面を示す．水色の面が設定した 音響面を表す．解析対象のRLGは床面により支持されてい る形態であるため，閉曲面を半分に割った形となる．音響 面が物体に近過ぎると，音響面を通過する渦に起因した偽 音の影響を大きく受け，また物体から遠過ぎると格子粗さ や数値粘性により音波が減衰してしまう事を考慮して領域 の大きさを決める．偽音の影響の緩和には複数の音響面か ら得た圧力変動値の平均を取る手法（Shurらの手法⁽¹⁰⁾（図 3））を用いる．Shur らの手法を適用出来る様に，後流側 にはx/D=2～8の間に間隔∆x/D=1で面を配置する．

図 4には床面での音波の反射の補正法について示す．無 限遠まで存在する床面で音波は全反射すると仮定し，床面 に対称な鏡像を考える．実像音響面の時系列データから求 める実像観測点と鏡像観測点の圧力変動を時系列に重ね合 わせる事により反射の影響を考慮する．

図2 FW-H法データ取得用音響面

音響面

観測点

Flow

図3 Shurらの手法⁽¹⁰⁾の概略図

Flow

X

Y Z

100D

8D 8D

Flow