太鼓山風力発電所のナセル落下事故に対する数値流体力学的アプローチによる一考察
Numerical Investigation of Terrain-Induced Turbulence in Taikoyama Wind Farm内田 孝紀*1
Takanori UCHIDA*1

Associate Professor, Research Institute for Applied Mechanics (RIAM), Kyushu University, 6-1 Kasuga-kouen,
Kasuga, Fukuoka, 816-8580, Japan

Abstract
Because a significant portion of the topography in Japan is characterized by steep, complex terrain,
which results in a complex spatial distribution of wind speed, great care is necessary for selecting a sitefor the construction of Wind Turbine Generators (WTGs). We have developed a Computational FluidDynamics (CFD) model for unsteady flow called Research Institute for Applied Mechanics, Kyushu
University, COMputational Prediction of Airflow over Complex Terrain (RIAM-COMPACT®). TheRIAM-COMPACT® CFD model is based on Large-Eddy Simulation (LES) technique. In this paper, thenumerical wind simulation over the Taikoyama wind farm was executed using the high resolution
elevation data. As a result, the numerical results obtained suggest that all the six WTGs in the Taikoyama
wind farm are subject to significant influence from separated flow (terrain-induced turbulence) which isdue to the topographical irregularity in front of WTGs. The new reproduction proposal was also done.
キーワード：太鼓山風力発電所，複雑地形，地形乱流，LES
Key Words：Taikoyama wind farm, Complex terrain, Terrain-Induced Turbulence, LES
１．はじめに
太鼓山風力発電所(京都府与謝郡伊根町字野村小字太鼓山地内)において，平成25年3月12日19時30分頃，3号機の発電機付きブレード（翼）部分（高さ約50.0ｍ，重さ約45.0ｔ）が地上に落下する重大事故が発生した．
破断したのは支柱とブレードを結合する根元の部分で，事故当時は風速15.0ｍ/ｓ程度と設計基準内に収まっており，「金属疲労」が事故の主原因との見方で調査が進められた．既に報告書1)が公表され，研究論文2)も報告されている．
太鼓山風力発電所には，オランダのLAGERWAY社製の750ｋW風車が6基設置されており，最大出力は4,500ｋWである．
発電所計画時の年間発電電力量は8,549MWｈであった．京都府の公益外郭団体により，平成13年11月から運用が開始されている．総事業費は約15億円である．風況精査時における年間平均風速は地上高20.0ｍで5.4ｍ/ｓ(べき指数1/7で高度補正後：年間平均風速6.2ｍ/ｓ＠50.0ｍ)であった．
太鼓山風力発電所の位置，風車外形図，事故当時の写真を図1～図3に示す．
文献2)には，各風車位置の気流性状の詳細が取り上げられていない．そこで本研究では，風車周辺に発生している「地形乱流」の影響に注目して太鼓山風力発*1 九州大学応用力学研究所 准教授（〒816-8580 春日市
春日公園 6-1）E-mail:takanori@riam.kyushu-u.ac.jp（原稿受付：年月日）

電所を対象にしたLESによる非定常・数値風況診断3, 4)を実施した．得られた数値シミュレーション結果に基づき，各風車位置における気流性状の相対比較を行った．特に，気流の3次元的な鉛直構造と時系列データから評価した標準偏差の3成分(主流方向，主流直交方向，鉛直方向)の値に注目した．本論文では，「地形乱流」を「地形の凹凸により発生する気流の時間的・空間的な変動である」と定義する．
２．数値計算手法の概要
本研究では，数値不安定を回避し，複雑地形上の局Google EarthTaikoyama wind farmNorth

Fig.1 Location of Taikoyama wind farm in Kyoto Prefecture

No.6

No.4
No.5

No.3

No.2
No.1

Fig.2 Picture of wind turbine accident, quoted from RKB Mainichi Broadcasting Corporation
地的な風の流れを高精度に数値予測するため，一般曲線座標系のコロケート格子に基づいた実地形版
RIAM-COMPACT®を用いた3, 4)．ここでコロケート格子とは，計算格子のセル中心に物理速度成分と圧力を定義し，セル界面に反変速度成分にヤコビアンを乗じた変数を定義する格子系である．数値計算法は差分法(FDM；Finite-Difference Method)に基づき，乱流モデルとしてLES(Large-Eddy Simulation)を採用する．
LESでは流れ場に空間フィルタを施し，大小様々なスケールの乱流渦を計算格子よりも大きなGS(GridScale)成分の渦と，それよりも小さなSGS(Sub-GridScale)成分の渦に分離する．GS成分の大規模渦はモデルに頼らず直接数値シミュレーションを行う．一方で，SGS成分の小規模渦が担う，主としてエネルギー消散作用はSGS応力を物理的考察に基づいてモデル化する．
流れの支配方程式は，フィルタ操作を施した非圧縮流体の連続の式(式(1))とナビエ・ストークス方程式(式(2))である．本研究では，平均風速5.0～6.0ｍ/ｓ以上の強風場を対象にしているので，大気が有する高度方向の温度成層(大気安定度)の効果は省略した．地
表面粗度の影響は，文献3, 4)と同様，地表面の凹凸を高解像度に再現することで取り入れた．計算アルゴリズ
ムは部分段階法(F-S法)5)に準じ，時間進行法はオイ
ラー陽解法に基づく．圧力に関するポアッソン方程式
は逐次過緩和法(SOR法)により解く．
空間項の離散化は式(2)の対流項を除いてすべて2次精度中心差分とし，対流項は3次精度風上差分とする．ここで，対流項を構成する4次精度中心差分は，梶島による4点差分と4点補6)間に基づいた補間法 を用いる．3次精度風上差分の数

値拡散項の重みは，通常使用されている河村-桑原ス
キーム7)タイプのα=3.0に対して，α=0.5とし，その影響は十分に小さくする．
LESのサブグリッドスケール
モデルには標準スマゴリンスキーモデル8)を用いる．
壁面減衰関数を併用し，モデル係数は0.1とした．

Fig.3 Schematic outline of wind turbine
u i
0xi(1)

 ij
u i
u i p 1 2 uiuj

t
x j
xi Re x j x j x j

(2)

13

 ij  u'i u' j  u'k u'k  ij  2 SGS S ij

(3)

 SGS   Cs f s   S(4)

2





(5)

1  ui u j 
2  x j xi (6)S  2Sij Sij

Sij 1/ 2

f s  1  exp  z  / 25

(7)

   hx hy hz 1/ 3(8)Fig.5 Characteristic two scales (Uin and h)３．数値計算条件の概要
本研究における計算領域は，主流方向(x)，主流直交方向(y)，鉛直方向(z)方向に10.0(x)×5.0(y)×3.3(z)ｋｍの空間を有する．計算領域中の最大標高は681.0ｍ，最小標高は8.0ｍである．地形標高データは国土地理院の空間解像度10.0ｍのデータを使用した．計算格子数は各方向に401(x)×201(y)×41(z)点の合計約330万点である．x方向およびy方向の格子幅は，風車周辺で密になるように不等間隔に分布させた．z方向の格子幅においても不等間隔とし，地表面付近に滑らかに引き寄せた．水平方向の最小格子幅は10.0ｍ，鉛直方向の最小格子幅は2.0ｍである．図4に風車近傍の計算格子図を示す．資料9)によると，この地域の卓越風向は西南西と報告されている．本数値シミュレーション
でも西南西を対象とした．境界条件に関して，流入境界面にはべき乗則(N値=7)に従う風速分布を与えた．
側方境界面と上部境界面は滑り条件，流出境界面は対流型流出条件とした．地表面には粘着条件を課した．
式(2)の無次元パラメータReはレイノルズ数(=Uinh/ν)であり，文献10)を参考にして本シミュレーションではRe=104とした．代表スケールの取扱いは図5に示す通りである．hは計算領域内の標高差(=673.0ｍ)，Uinは流入境界面の最大標高位置における風速，νは動粘性

WSW

Fig.4 Computational grids near wind turbines係数である．時間刻みはΔt=2×10-3h/Uinとした．
４．計算結果と考察
最初に，1号機～6号機のすべての風車に関して，ハブ高さ位置(地上高50.0ｍ，図3を参照)で算出した風速
3 成 分 に よ る 無 次 元 ス カ ラ ー 風 速 Uscalar/Uin(=(u2+v2+w2)1/2/Uin)の時間変化を図6に示す．ここで横軸は 無 次 元 時 間 (=T/(h/Uin)) を 示 す ． 実 風 速 と し て
Uin=5.0ｍ/ｓを仮定すると，横軸は実スケールで約45分間に相当する．図6を観察すると，特異な流動現象が風車周辺において発生していることが見て取れる．すなわち，1号機～6号機の無次元スカラー風速の時間変化は，すべての風車でほぼ同じ傾向を示している．その波形は低風速から高風速を周期的に繰り返しながら推移している．後ほど詳しく議論するが，これらの周期的に推移する波形は，風車近傍に存在する地形の凹凸から発生している地形乱流が風車を通過していることを示唆している．すなわち，1号機～6号機のすべての風車は，地形の凹凸を発生起源とする地形乱流の影響を日常的に受けていたのではないかと推測される．
偶然に3号機のナセルの落下事故が発生したが，今回の事故は起こるべくして起きたとも考えられ，他の風車のナセルが落下していても不思議ではなかったのではないだろうか．なお，事故後の調査では1号機を除くすべての風車で3号機と同様な亀裂が見つかっている．
図6に示したような非定常な流動現象は，乱流モデルにRANS(Reynolds-Averaged Navier-Stokes)モデルを採用した数値風況シミュレーションでは，支配方程式をレイノルズ平均(アンサンブル平均)しているため，再現することが出来ない．図7には，図6で示した無次
元時間t=100～120の時間平均(フレーム平均)を施した風況場から得られた各風車立地点での風速の鉛直分布を示す．これはRANSモデルの出力結果に準ずる．縦軸のz*は地上高(ｍ)を，横軸は正規化した風速を表わしている．ここで注意すべきは次の点であろう．すなわ
ち，図7のみから判断すると，太鼓山風力発電所におけNon-dimensional
scalar wind speed

Fig.6 Time variation of scalar wind speeds at hub height (50.0m),
Time interval of horizontal axis shows about 45 minutes (Uin=5.0m/s)
る1号機～6号機のすべての風車位置において，平均風速は地形効果によって局所的に増速されるものの大きな速度シアーが観察されていないので，風況の面からは1号機～6号機のすべての風車位置において重大な問題は生じないというような結論になりがちである．よって，地形乱流の影響を調査するためには，非定常な流動現象を踏まえての考察が不可欠である．
これ以降は，今回ナセルの落下事故が発生した3号機風車に注目して考察を進める．図8には，風車ハブ高さ位置(地上高50.0ｍ)における水平断面内(ヨー方向)および鉛直断面内に関して，風車に対する風の進入角度の時間変化を示す．○印が水平断面内(ヨー方向)の角度であり，□印が鉛直断面内の角度である．図8には，それぞれの角度の定義を併せて示す．二種類の角度ともに，図6に対応して周期的な時間変化を示していることが観察される．最初に，□印で表示した鉛直断面内の角度の時間変化に注目すると，±25度に近い吹上風が頻繁に発生しているのが分かる．特に図中のA点

75.25m

Swept area

50.0m

24.75mLocal speed-up bytopographic effects

Fig.7 Profiles of mean wind speeds at each of wind turbines(=+28.5度) およびB点(=+28.2度)に示すように，25度を越えるような大きな角度を有する吹上風も発生している．次に，○印で表示した水平断面内(ヨー方向)の角度の時間変化に注目する．
鉛直断面内の角度と同様，±25度前後の角度が周期的に出現している．ここで，図中のC点(=+28.1度，無次元時間106.0)およびD点(=-40.4度，無次元時間107.5)における風車立地点における速度ベクトル図を図9に示す．また，これに対応した風速の鉛直分布を図10に示す．図9のSide Viewおよび図10を吟味すると，C点およびD点の両時間ともに風車受面内では，図10に黒太線で示すべき乗則(N値=7)に従う流入風速分布と比較して，これを大きく逸脱するような風速の鉛直分布は出現していない．しかしながら，図9のRear Viewに注目すると，C点では速度ベクトルの鉛直分布が大きくねじれていることが分かる．
一方，D点では風車ブレードの上端から地表面において速度ベクトルが大きく傾いている．
速度ベクトルの鉛直分布に見られるねじれや傾きは，地形乱流の3次元構造に起因している．結果として，風車支柱とブレードを結合する根元付近に付加的な荷重が作用し，ボルトの金属疲労をより進行させたと推測できる．
紙面の都合上ここでは示さなかったが，べき乗則(N値=7)に従う流入風速分布を大きく逸脱し，風車ハブ中心からブレード下端に大きな速度シアーが発生する時間帯も複数存在した．
図11には，図6で示した無次元時間t=100～120の時間平均(フレーム平均)を施した風況場から算出した各風車立地点での標準偏差3成分の鉛直分布を示す．
本研究では，流入気流が有する変動成分(風の息)は考慮しておらず，地形の凹凸を発生起源として生成された地形乱流の気流変動のみを評価している．図11から，標準偏差の各成分ともに，すべての風車位置の受風面内において比較的大きな値を示しているのが分かる．
また，図11(a)に示すx方向成分と，図11(b)に示すy方向成分の値はほぼ同程度の値を示していることは特筆すべき点である．これは図8および図9で述べたように，すべての風車位置において水平断面内(ヨー方向)の気流のA

BAngle in the x-z plane (β)

CAngle in the x-y plane (α)

DFig.8 Time variation of angles in the x-y and x-z plane at hub
height (50.0m), Time interval of horizontal axis shows about
45 minutes (Uin=5.0m/s)

Side View

Side View

Rear View

Rear View

75.25m

Swept area

50.0m(a) Time=106.0, C-point in Fig.8

24.75m(b) Time=107.5, D-point in Fig.8

Fig.9 Instantaneous velocity vectors at the wind turbine of No.3
時間的・空間的な変動が大きいことを表している．以上から，x方向成分の比較的大きな標準偏差が出現していることに加えて，z方向成分にも大きな値が示されていること，さらにy方向成分の標準偏差がx方向成分のFig.10 Profiles of instantaneous wind speeds
corresponding to Fig.9

75.25m(a) x-component

Swept area

50.0m

24.75m(b) y-component(c) z-component

Fig.11 Profiles of standard deviations at each of wind turbines

15 (％)

GPV1

GPV2

Taikoyama wind farm
GPV1

North15 (％)

Occurrence frequency ofwind direction (％)

Taikoyam (683.1m)

GPV2(135.25deg, 35.7deg)(135.1875deg, 35.7deg)
5km

5km

GPV3

GPV4(135.1875deg, 35.65deg)(135.25deg, 35.65deg)

15 (％)

GPV315 (％)

Fig.12 Wind characteristics around Taikoyama windfarm based on weather GPV data
(MSM-S data obtained in 2012)

値と同程度の値を有していることに起因して，上下左右方向の励振力が発生したのではないかと考えられる．
その結果，風車本体と支柱接合部のボルトの破損，それに伴ってさらなる励振力の増進で風車の支柱上部が疲労破壊したという解釈も可能ではないだろうか．
さらに本研究では，気象庁配信のGPVデータ(MSM-S，地上高10.0ｍ)を用いて太鼓山風力発電所の周辺における風況特性を解析した．その結果，先に記述した西寄りの風と同様，平成22年～平成24年の3年間では年間を通じて南寄りの風が頻繁に発生しているという結果を得た．一例として平成24年の解析結果を図12に示す．いずれの風配図においても，南寄りの風の頻度が高くなっている．
上記の解析結果を踏まえて，南風を対象にしたRIAMCOMPACT®による非定常・数値風況診断を実施したところ(図13を参照)，すべての風車が図中に矢印Aで示す付近から発生した地形乱流の影響を強く受け，風速と風向のGPV4両者が時間とともに著しく変化する気流中に埋没して稼動しているように推察された．さらに，1号機から5号機の風車は南北にほぼ直線的に並んでいるので，太鼓山風力発電所の上空に南風が発生すると，風車ウエイクの相互干渉（前方に位置する風車ブレードの回転により発生した乱気流が後方の風車に影響を及ぼし，風車の故障や発電量の低下を引き起こす現象）も懸念されることが見出された(図14に示すイメージ図を参照)．
以上より，太鼓山風力発電所に南風が発生した場合には，地形乱流と風車ブレードが回転することに起因した乱気流(風車ウエイクの相互干渉問題)の両者の影響で，風車支柱とブレードを結合する根元付近に付加的な荷重が作用し，結果としてボルトの金属疲労をより進行させたのではないかと推測される．
５．おわりに
本研究では，太鼓山風力発電所を対象にLES乱流モ

Taikoyama (683.1m)
S-windAWind speed
-0.35
1.50

North

Computed using RIAM-COMPACT LES Model

Fig.13 Distribution of streamwise velocity components along vertical cross-sections
which include the wind turbine of No.1, S-wind

Wind turbine

5DWind turbine(D : Diameter of swept area)
FlowWind speed
-0.35
1.50

Fig.14 Image of wake interaction between two wind turbines by visualizing Laplacian of pressure field
デルによる非定常・数値風況診断を実施した．得られたシミュレーション結果に基づいて，風車周辺に発生している「地形乱流」の影響に注目して，数値風況面からの一考察を行った．西南西の風の場合，ナセル落下事故の起きた3号機を含むすべての風車位置で，べき乗則(N値=7)に従う風速分布から逸脱した速度シアーが存在すること，また大きな風向変化が風車のヨー方向で頻繁に発生しているという知見を得た．特に，主流方向(x)の比較的大きな標準偏差に加えて，鉛直方向(z)へも大きな値をとっていること，さらに主流直交方向(y)の標準偏差がx方向成分の値とほぼ同程度の値を有していることも認められた．以上から，地形乱流の影響により励振力が増大し，結果として風車支柱とブレードを結合する根元付近に付加的な荷重が発生し，ボルトの金属疲労がより進行したのではないかと推測される．南風の場合には，地形乱流に加えて，風車ブレードが回転することに起因した乱気流(風車ウエイクの相互干渉)の影響が現れていることが見出された．
謝 辞本研究の一部は，科学研究補助金 挑戦的萌芽研究(課題番号：25560166)の援助を受けました．ここに記して感謝の意を表します．

参考文献

付 録1)

本研究の数値風況シミュレーション結果に基づくな
らば，太鼓山風力発電所の「再生」として，以下の方向性が検討に値すると思われる．既存の風車はすべて撤去し，地形乱流に強い風車を選定して1～2機程度を建設する．その際，風車支柱を出来るだけ高くすることが望ましい．地形乱流は基本的に地表面に近いところで発生・出現するので，支柱が高くなると地形乱流の影響は劇的に減少するからである．次に，図12に示した平成24年のGPV1位置における風速の時系列データ(地上高10.0ｍ)を高度補正し，これを用いて支柱高さ70.0ｍの風車が1機設置された場合の経済性に関わる試算を行った（付図1を参照)．得られた結果は，日本国内の風力発電所との比較から平均的であるとみなされるので，事故報告書1)にも提示されているように適切な維持・管理を実施することで，太鼓山風力発電所の「再生」は十分に実現可能と思われる．

京都府太鼓山風力発電所3号機ナセル落下事故報告書,
http://www.safety-kinki.meti.go.jp/denryoku/2013/
taikoyama5.htm, 2013.
2) Liu, Z.Q. and Ishihara, T., The Sixth International Symposiumon Computational Wind Engineering, CWE2014, CD-ROM,2014.
3) 内田孝紀, 他4名, 風力エネルギー協会誌, Vol.34, No.96,
2011, pp.77-84.
4) 内田孝紀・大屋裕二, 土木学会論文集A2 (応用力学),
Vol.67, No.2, (応用力学論文集Vol.14), 2011, I_757-I_764.
5) Kim, J. and Moin, P., J. Comput. Phys., Vol.59, 1985,
pp.308-323.
6) 梶島岳夫・太田貴士・岡崎和彦・三宅裕, 日本機械学会
論文集, (B編), 63巻, 614号, 1997, pp.3247-3254.
7) Kawamura, T., Takami, H. and Kuwahara, K., Fluid Dyn.
Res., Vol.1, 1986, pp.145-162.
8) Smagorinsky, Mon. Weather Rev., Vol.91, 1963, pp.99-164.
9) 第2回京都府風力発電事業評価委員会資料,
http://www.pref.kyoto.jp/shingikai/koei-02/documents/1291615930487.pdf, 2008.
10) 加藤真規子, 日本風工学会誌, 第59号, 1994, pp.89-92.

about 6.0 m/s @ 70m
Power output (kW)Frequency of occurrence (％)

Annual Average Wind Speed :Wind speed (m/s)

Wind speed (m/s)
Annual Energy Production (kWh)Capacity Factor (％)

Income by Electric Power Selling (Yen)

750kW
29,104,891

1,000kW35,315,003

1,500kW
58,048,175

2,000kW
79,653,460

Remarks) We assume that availability=90 (％) and Feed-in Tariff=22 (Yen/kWh).
Appended Fig.1 Results of trial calculation of economical efficiency in the case of one wind turbine(hub height=70.0m) based on weather GPV data (MSM-S data at a position of GPV1 in Fig.12)