Bolund Experimentを対象にした
LES乱流モデルの予測精度の検証―風洞実験との比較―

九州大学応用力学研究所
新エネルギー力学部門 風工学分野
内田 孝紀
連絡先：takanori@riam.kyushu-u.ac.jp，092-583-7776
１. はじめに対し，独自に別途，風洞実験も実施した．本報では，現在，世界では空前の勢いで風力発電産業が成 我々の研究グループが開発を進めるLES系の数値長を遂げている．これは再生可能エネルギーの中で 風況予測技術RIAM-COMPACT® の結果と，風洞風力発電が最も脱化石燃料，CO2削減に対してコス 実験の気流計測結果 の比較を報告する．
トパフォーマンスが高いからである．日本においても風力発電が最も有力な再生可能エネルギーであるこ ２．実地形版RIAM
実地形版RIAMRIAM-COMPACT®
COMPACT®ソフトウエアのソフトウエアの概要とは間違いなく，風力発電をより一層普及させること
RIAM-COMPACT®(Research Institute for Appliedが，地球温暖化の克服「グリーン・イノベーション」に Mechanics, Kyushu University, COMputational
Prediction of Airflow over Complex Terrain；リアムコ
世界的な規模で貢献すると確信する．
風力発電分野においても，風車の配置計画や年 ンパクト) に関して，そのコア技術は九州大学応用力間発電電力量の試算などに対して，数値流体力学 学研究所で開発が続けられており，2006年に著者らCFD(Computational Fluid Dynamics)の利用が急速 が起業した九州大学発ベンチャー企業の(株)リアムに進んでいる．こうした状況を受けて，風力発電分野 コンパクト(http://www.riam-compact.com/)が(株)産で利用されているCFDソフトウエアの予測精度の検 学連携機構九州(九大TLO)から独占的ライセンス使証(Validation Testing)が各方面で進められている．そ 用許諾を受けている(2006年にRIAM-COMPACT®の一つに，Bolud Experimentがある(図1を参照)．これ の商標と実用新案を取得)．現在では，九州電力グルはデンマーク工科大学(DTU)のリソ研究所が中心と ープの西日本技術開発(株)，(株)環境GIS研究所となり，海上にある孤立地形を過ぎる流れ場を対象とし 開発コンソーシアムを作り，「実地形版RIAMて実施したプロジェクトである．このプロジェクトでは， COMPACT®ソフトウエア」と名付け，業界標準モデ
風況観測（カップ式，超音波式風速計およびLIDAR） ルの一つとして広く普及に努めている．現在では，国および数値シミュレーション(CFD)が実施され，CFD 内の風力事業者最大手の(株)ユーラスエナジーホーの予測精度が議論されている．ここでは，非定常・非 ルディングス，電源開発(株)，日本風力開発(株)，エ線形風況解析モデルであるLES（空間平均モデル） コ・パワー(株)を含め，多数の導入実績を有する．
よりも，定常・非線形風況解析モデルであるRANS（時非定常な乱流シミュレーションに主眼を置いた間平均モデル）が比較的高い評価を受けている．一 RIAM-COMPACT®では，計算時間の問題が懸念さ方，日本ではLES系のモデルも高い評価を受けてい れてきた．現行の流体計算ソルバーは，Intel Core i7ることから，改めてその評価を実施した．この目的に などのマルチコアCPU(Central Processing Unit)に対応しており，計算時間は劇的に短縮され，実用面での利用において特段の問題は無くなってきた ．
さらに，現在ではGPGPUコンピューティングへの対応 に も 成 功 し た ． GPGPU(General Purposecomputing on GPU：GPUによる汎用計算)のコンセプ
トとは，グラフィック・レンダリングのみならず，
図1 Bolund Experimentにおける対象地形，点線部分 GPU(Graphics Processing Unit)が有する浮動小数点
1)

2)

1)

3)

4, 5)

演算能力を，他の数値演算にも幅広く利用すること∆ = (h h h )
-(8)
である．
RIAM-COMPACT®の数値計算手法を説明する．
数値不安定を回避し，複雑地形上の局所的な風の計算アルゴリズムは部分段階法(F-S法) に準じ，流れを高精度に数値予測するため，変数配置には 時間進行法はオイラー陽解法に基づく．圧力に関す
一般曲線座標系のコロケート格子を採用している．コ るポアッソン方程式は逐次過緩和法(SOR法)により解
ロケート格子では，計算格子のセル中心に物理速度 く．空間項の離散化は式(2)の対流項を除いて全て2成分と圧力を定義し，セル界面に反変速度成分にヤ 次精度中心差分とし，対流項は3次精度風上差分とコビアンを乗じた変数を定義する．数値計算法は(有 する．ここで，対流項を構成する4次精度中心差分は，限)差分法(FDM；Finite-Difference Method)に基づき， 梶島による4点差分と4点補間に基づいた補間法 を乱流モデルにはLES(Large-Eddy Simulation)を採用 用いる．3次精度風上差分の数値拡散項の重みは，する．LESでは流れ場に空間フィルタを施し，大小 通常使用される河村-桑原スキーム タイプのα=3に様々なスケールの乱流渦を，計算格子よりも大きな 対して，α=0.5とし，その影響は十分に小さくする．
GS(Grid Scale) 成 分 の 渦 と ， そ れ よ り も 小 さ な LESのサブグリッドスケールモデルには標準スマゴリSGS(Sub-Grid Scale)成分の渦に分離する．GS成分 ンスキーモデル を用いる．壁面減衰関数を併用し，の大規模渦はモデルに頼らず直接数値シミュレーシ モデル係数は0.1とした．
ョンを行う．一方で，SGS成分の小規模渦が担う，主としてエネルギー消散作用はSGS応力を物理的考察に ３．地形標高データ
地形標高データおよびデータおよび計算条件および計算条件など計算条件など基づいてモデル化される．
ここでは，本研究で使用した地形標高データおよ流れ場の支配方程式は，空間フィルタ操作を施さ び計算条件などについて説明する．図2には，地形れた非圧縮流体の連続の式(式(1))とナビエ・ストーク 標高データに対するコンターマップを示す．図中にス方程式(式(2))である．本研究では，約5m/s以上の強風場を対象にしているので，大気が有する高度方向の温度成層(密度成層)の効果は省略した．
1/3x y z6)

7)

8)9)

-(1)

∂ui=0∂x i

-(2)

∂τij
∂ui
∂ ui
∂p
1 ∂2 ui+ uj=−+−∂t∂x j∂x i Re ∂x j∂x j ∂x j

1 ' 'u k u k δij − 2νSGS Sij
3

τij ≈ u 'i u 'j ≈-(4)

2

νSGS = ( Cs fs ∆ ) S(S = 2Sij Sij

Sij =-(6)1  ∂u i ∂u j +
2  ∂x j ∂x i (

図2 本研究で使用した地形標高データなど

-(5)

1/2)

fs = 1 − exp −z + / 25-(3))-(7)

図3 気流の鉛直プロファイルを評価した地点など

記載しているように，オリジナルの地形データはリソ研 分布を流入境界面に与えた．ここで，本研究におけ究所から配布されているデータを使用し，これに基づ る風速の代表スケールU は，流入境界面の最大標いて空間解像度1mのDEMデータを作成した．計算 高H(=11m)位置における風速値である．このU とHにを実施した風向は，北をゼロ度として239度である． 基づいたレイノルズ数Re(=U H/ν)は，風洞実験に
図3には，気流の鉛直プロファイルを評価した地点 合わせて2.8×10 とした．などを示す．本研究では，Bolud Experimentで定義されたAライン上の4点(M1, M2, M3, M4)において， ４．風洞実験の風洞実験の概要
RIAM-COMPACT®と風洞実験で得られた平均風速 ここでは，本研究で実施した風洞実験の概要につ分布を比較した．図中に記載しているように，Aライン いて説明する ．図5には，風洞実験設備の概要を示
上のM0地点で得られた風洞実験の平均風速分布を す．本研究で使用した風洞装置は，東京大学生産技数値シミュレーションの流入境界面に与えた．
術研究所所有の境界層型風洞であり，長さ16.5m，図4には，計算領域，計算格子，座標系などを示 幅2.2m，高さ1.8mの測定胴を有する(詳細はhttp://す．計算領域は，計算領域中の最大標高H(=11m)を venus.iis.u-tokyo.ac.jp/equip/を参照して頂きたい)．
基準にし，主流方向(x)および主流直交方向(y)に気流場の計測には，逆流の検出が可能なスプリッ
38.55H四方，鉛直方向(z)に15.45Hの空間を有する． トファイバープローブを使用した(図6を参照)．サンプ
特に鉛直方向は，ブロッケージの影響を避けるため， リング間隔は1kHzとし，計測時間は60sとした．
風洞実験と同じ条件に設定した．計算格子はx方向
図7には，本研究で作成した地形模型を示す．地
およびy方向に213点，z方向に65点の合計約300万 形模型の縮尺は1/150であり，最大標高H(=11m)は，点を設定した．空間解像度に関して，x-y断面は⊿x= 風洞実験では約7cmとなる(図7(a)を参照)．図7(b)に⊿y=2m(=0.18H)の等間隔とし，z方向は地面付近で は地形模型の拡大写真を示す．図8には，気流計測密になるように不等間隔⊿z=0.038m～7.64m(0.003H の様子を示す(Rear View)．地形模型の上流側にスパ～0.7H)とした．流入条件は，図3で既に述べたように， イヤーとラフネスブロックを設置し，1/7のべき乗則にAライン上のM0地点で得られた風洞実験の平均風速 従う境界層流を生成した．
refrefref42)

図4 計算領域，計算格子，座標系など

(a)全体図(ｂ)拡大図図7 地形模型の様子
図5 本研究で使用した境界層型風洞（測定部長さ16.5m，幅2.2m，高さ1.8m），東京大学生産技術研究所所有

図8 気流計測の様子，Rear View
実験を参考にして地表面からz*≦0.0286Hまでの空気層とした(風洞実験スケールで2mm)．ここで，Hは本研究の代表長さスケールである最大標高(=11m，図6 逆流の検出が可能な風洞実験では7cm)である．図9(a)の滑面ケースおよスプリットファイバープローブ
び図9(b)の粗面ケースの両者を観察すると，地形背５．結果および後に渦領域が形成されているのが分かる．図9(b)で結果および考察
および考察
こ こ で は ， 数 値 シ ミ ュ レ ー シ ョ ン (RIAM- は，地表面粗度の影響により，図中に示すように渦領COMPACT® LES乱流モデル)と風洞実験の比較に 域の伸長が明確に示されている ．
ついて説明する．
図10には，図9に示すM1～M4の地点における主図 9 に は ， 数 値 シ ミ ュ レ ー シ ョ ン (RIAM- 流方向(x)の速度成分の鉛直プロファイルの比較を示
COMPACT® LES乱流モデル)に関して，時間平均 す．ここで，ラインが風洞実験(滑面)，シンボル○が風況場(無次元時間t =100～200)の可視化図(主流 LES数値シミュレーション(滑面)，シンボル＋がLES数方向(x)の速度成分の分布図)を示す．比較のため， 値シミュレーション(粗面)をそれぞれ意味する．全て地表面粗度を考慮したケースも合わせて表示する． の地点，高さレベルにおいてLES数値シミュレーショ
地表面粗度は文献 に記載している外力項(キャノピ ンと風洞実験の結果は良好な一致を示しているのがーモデル)を支配方程式に付加することで表現した 分かる．また，LES数値シミュレーションにおけるM1(抵抗係数Cd=10)．粗度モデルの適用範囲は，風洞 ～M3地点の平均風速の鉛直プロファイルに関して，
10)

*10)

M1地点 M2地点

M3地点

M4地点

渦中心

Flow，風向239度(a)地表面粗度無しのケース(滑面)渦領域の伸長渦中心(b)地表面粗度有りのケース(粗面)
図9 数値シミュレーション(RIAM-COMPACT® LES乱流モデル)の結果，
時間平均風況場(無次元時間t =100～200)，主流方向(x)の速度成分の分布図*

図10 数値シミュレーション(RIAM-COMPACT® LES乱流モデル)と風洞実験の比較，
図9に示すM1～M4の4地点，主流方向(x)の速度成分の鉛直プロファイル，ライン：風洞実験(滑面)，○：LES数値シミュレーション(滑面)，＋：LES数値シミュレーション(粗面)地表面粗度の有無による有意な差異は見られなかった．地表面粗度の影響は，図9で述べたように，地形背後の渦領域の大きさに顕著に現れている．これに伴い，M4地点の地面近傍において平均風速の鉛直プロファイルの値に違いが確認される．

上の気流性状に関する風洞模型実験，日本建築学会学術講演梗概集（東海），2012年9月
3)Bridge-EP 」 を 用 い た RIAM-COMPACT® のOpenMP並列風況シミュレーション―ここまできた！
デスクトップPC1台による大規模計算とその高速化

６. おわりに本研究では，LES乱流モデル(CFD)に基づいた実
地 形 版 RIAM-COMPACT® を 用 い て ， BolundExperimentを対象に数値風況予測を実施した．また同時に，東京大学生産技術研究所(加藤研究室)の境界層風洞にて，スプリットファイバープローブによる気流計測を実施した．Bolund Experiment観測サイト
周辺の気流性状を明らかにするとともに，LES乱流シミュレーション(CFD)の予測精度を検証した．その結果，CFDによる予測値は，全ての地点，高さレベルにおいて風洞実験結果と良好な一致を示し，LES乱流モデルに基づいたCFDの有効が明らかになった．
謝辞

本研究の一部は，NEDO「風力等自然エネルギー技術研究開発／次世代風力発電技術研究開発／基礎・応用技術研究開発（数値シミュレーション技術を用いた風車性能評価技術等の国際標準化に係る研究開発），2010年度～2012年度」，「科学研究費補助金・基盤研究(B) (課題番号：24310120)，2012度～2014度」による援助を受けました．ここに記して感謝の意を表します．
参考文献
1)

Uchida, T. and Ohya, Y. : Micro-siting Techniquefor Wind Turbine Generators by Using Large-Eddy
Simulation, Journal of Wind Engineering & Industrial
Aerodynamics, Vol.96, pp.2121-2138, 2008.

2)

高橋岳生，加藤信介，中尾圭佑，飯田誠：孤立丘

内 田 孝 紀 ： イ ン テ ル 次 世 代ハ イ エ ン ド CPU「 Sandy―，日本風力エネ ルギー学会論文集，Vol.36，通巻101，pp.82-85，20124)

内田孝紀，大屋裕二：GPGPUコンピューティングによる数値風況予測技術RIAM-COMPACT®の高速化 ， 風 力 エ ネ ル ギ ー 協 会 誌 ， Vol.35 ， 通 巻 98 ，pp.78-84，20115)

内田孝紀：GPGPUコンピューティングによる数値風況予測技術RIAM-COMPACT®の高速化―第2報
マルチGPUによる大規模計算の試み―，日本風力エネルギー学会誌，通巻103，20126)

Kim,J.

andMoin,fractional-step

P. ： Application

methodtoofa

incompressible

Navier-Stokes equations, J. Comput. Phys., Vol.59,pp.308-323, 1985.
7)

梶島岳夫，太田貴士，岡崎和彦，三宅裕：コロケー
ト格子による非圧縮流れの高次差分解析，日本機械学会論文集，(B編)，63巻，614号，pp.3247-3254，
1997.

8)

Kawamura,T., Takami,

H.

andKuwahara,

K.

:

Computation of high Reynolds number flow around acircular cylinder with surface roughness, Fluid Dyn.
Res., Vol.1, pp.145-162, 1986.
9)

Smagorinsky, J. : General circulation experiments
withtheprimitiveexperiments,Mon.

equations,PartWeatherRev.,

1,

Basic
Vol.91,

pp.99-164, 1963.
10) 内田孝紀，大屋裕二：一様流中の2次元崖状地形
まわりの気流性状に関する数値的研究，応用力学
論文集，Vol.8，pp.831-838，2005