風車の重大事故防止に資する乱気流診断システムの開発
九州大学 応用力学研究所 風工学分野 准教授 内田 孝紀
１．はじめに太鼓山風力発電所(京都府与謝郡伊根町字野村小字太鼓山地内)において，平成25年3月12日19時30分頃，3号機の発電機付きブレード（翼）部分（高さ約50.0ｍ，重さ約45.0ｔ）が地上に落下する重大事故が発生した．破断したのは支柱とブレードを結合する根元の部分で，事故当時は風速15.0ｍ/ｓ程度と設計基準内に収まっており，「金属疲労」が事故の主原因との見方で調査が進められた．
太鼓山風力発電所には，オランダのLAGERWAY社製の750ｋW風車が6基設置されており，最大出力は4,500ｋWである．発電所計画時の年間発電電力量は8,549MWｈであった．京都府の公益外郭団体により，平成13年11月から運用が開始されている．総事業費は約15億円である．風況精査時における年間平均風速は地上高20.0ｍで5.4ｍ/ｓ(べき指数1/7で高度補正後：年間平均風速6.2ｍ/ｓ＠50.0ｍ)であった．太鼓山風力発電所の位置，風車外形図，事故当時の写真をFig.1～Fig.3に示す．
本研究では，風車周辺に発生している「地形乱流」の影響に注目して太鼓山風力発電所を対象にしたLESによる非定常・数値風況診断を実施した．得られた数値シミュレー
ション結果に基づき，各風車位置における気流性状の相対比較を行った．特に，気流の3次元的な鉛直構造と時系列データから評価した標準偏差の3成分(主流方向，主流直交方向，鉛直方向)の値に注目した．本論文では，「地形乱流」を「地形の凹凸により発生する気流の時間的・空間的な変動である」と定義する．

Fig.1 Location of Taikoyama wind farm in Kyoto Prefecture－ 219 －Fig.2 Picture of wind turbine accident, quoted from RKB Mainichi Broadcasting
CorporationFig.3 Schematic outline of wind turbine２．数値計算手法の概要
本研究では，数値不安定を回避し，複雑地形上の局地的な風の流れを高精度に数値予測するため，一般曲線座標系のコロケート格子に基づいた実地形版RIAM-COMPACT®
を用いた．ここでコロケート格子とは，計算格子のセル中心に物理速度成分と圧力を定義し，セル界面に反変速度成分にヤコビアンを乗じた変数を定義する格子系である．数値計算法は差分法(FDM；Finite-Difference Method)に基づき，乱流モデルとしてLES(Large-Eddy Simulation)を採用する．LESでは流れ場に空間フィルタを施し，大小様々なスケールの乱流渦を計算格子よりも大きなGS(Grid Scale)成分の渦と，それよりも小さなSGS(Sub-Grid Scale)成分の渦に分離する．GS成分の大規模渦はモデルに頼らず直接数値シミュレーションを行う．一方で，SGS成分の小規模渦が担う，主として－ 220 －エネルギー消散作用はSGS応力を物理的考察に基づいてモデル化する．
流れの支配方程式は，フィルタ操作を施した非圧縮流体の連続の式(式(1))とナビエ・ストークス方程式(式(2))である．本研究では，平均風速5.0～6.0ｍ/ｓ以上の強風場を対象にしているので，大気が有する高度方向の温度成層(大気安定度)の効果は省略した．
地表面粗度の影響は，地表面の凹凸を高解像度に再現することで取り入れた．計算アルゴリズムは部分段階法(F-S法)に準じ，時間進行法はオイラー陽解法に基づく．圧力に関するポアッソン方程式は逐次過緩和法(SOR法)により解く．空間項の離散化は式(2)の対流項を除いてすべて2次精度中心差分とし，対流項は3次精度風上差分とする．ここで，対流項を構成する4次精度中心差分は，梶島による4点差分と4点補間に基づいた補間法を用いる．3次精度風上差分の数値拡散項の重みは，通常使用されている河村-桑原
スキームタイプのα=3.0に対して，α=0.5とし，その影響は十分に小さくする．LESのサブグリッドスケールモデルには標準スマゴリンスキーモデルを用いる．壁面減衰関数を併用し，モデル係数は0.1とした．
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
2 SGS   Cs f s   S

(4)

(5)

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

Sij

(3)

1/ 2

fs 
1  exp  z  / 25

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

Fig.4 Computational grids near wind turbines

Fig.5 Characteristic two scales (Uin and h)－ 222 －を参照)で算出した風速3成分による無次元スカラー風速Uscalar/Uin(= (u2+v2+w2)1/2/Uin)の時間変化をFig.6に示す．ここで横軸は無次元時間(=T/(h/Uin))を示す．実風速としてUin=5.0ｍ/ｓを仮定すると，横軸は実スケールで約45分間に相当する．Fig.6を観察すると，特異な流動現象が風車周辺において発生していることが見て取れる．すなわち，1号機～6号機の無次元スカラー風速の時間変化は，すべての風車でほぼ同じ傾向を示している．その波形は低風速から高風速を周期的に繰り返しながら推移している．後ほど詳しく議論するが，これらの周期的に推移する波形は，風車近傍に存在する地形の凹凸から発生している地形乱流が風車を通過していることを示唆している．すなわち，1号
機～6号機のすべての風車は，地形の凹凸を発生起源とする地形乱流の影響を日常的に受けていたのではないかと推測される．偶然に3号機のナセルの落下事故が発生したが，今回の事故は起こるべくして起きたとも考えられ，他の風車のナセルが落下していても不思議ではなかったのではないだろうか．なお，事故後の調査では1号機を除くすべての風車で3号機と同様な亀裂が見つかっている．

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)
Fig.6に示したような非定常な流動現象は，乱流モデルにRANS(Reynolds-AveragedNavier-Stokes)モデルを採用した数値風況シミュレーションでは，支配方程式をレイノルズ平均(アンサンブル平均)しているため，再現することが出来ない．Fig.7には，Fig.6
で示した無次元時間t=100～120の時間平均(フレーム平均)を施した風況場から得られた各風車立地点での風速の鉛直分布を示す．これはRANSモデルの出力結果に準ずる．
縦軸のz*は地上高(ｍ)を，横軸は正規化した風速を表わしている．ここで注意すべきは次の点であろう．すなわち，Fig.7のみから判断すると，太鼓山風力発電所における1号機～6号機のすべての風車位置において，平均風速は地形効果によって局所的に増速されるものの大きな速度シアーが観察されていないので，風況の面からは1号機～6号機のすべての風車位置において重大な問題は生じないというような結論になりがちである．
よって，地形乱流の影響を調査するためには，非定常な流動現象を踏まえての考察が不可欠である．
これ以降は，今回ナセルの落下事故が発生した3号機風車に注目して考察を進める．
Fig.8には，風車ハブ高さ位置(地上高50.0ｍ)における水平断面内(ヨー方向)および鉛直－ 223 －Fig.7 Profiles of mean wind speeds at each of wind turbines

Fig.8 Time variation of angles in the x-y and x-z plane at hub height(50.0m),
Time interval of horizontal axis show about 45 minutes(Uin =5.0m/s)
断面内に関して，風車に対する風の進入角度の時間変化を示す．○印が水平断面内(ヨー方向)の角度であり，□印が鉛直断面内の角度である．Fig.8には，それぞれの角度の定義を併せて示す．二種類の角度ともに，Fig.6に対応して周期的な時間変化を示していることが観察される．最初に，□印で表示した鉛直断面内の角度の時間変化に注目す
ると，±25度に近い吹上風が頻繁に発生しているのが分かる．特に図中のA点(=+28.5度) およびB点(=+28.2度)に示すように，25度を越えるような大きな角度を有する吹上風も発生している．次に，○印で表示した水平断面内(ヨー方向)の角度の時間変化に注－ 224 －目する．鉛直断面内の角度と同様，±25度前後の角度が周期的に出現している．ここで，図中のC点(=+28.1度，無次元時間106.0)およびD点(=-40.4度，無次元時間107.5)における風車立地点における速度ベクトル図をFig.9に示す．また，これに対応した風速の鉛

Fig.9 Instantaneous velocity vectors at the wind turbine of No.3

Fig.10 Profiles of instantaneous wind speeds corresponding to Fig.9
直分布をFig.10に示す．Fig.9のSide ViewおよびFig.10を吟味すると，C点およびD点の両時間ともに風車受面内では，Fig.10に黒太線で示すべき乗則(N値=7)に従う流入風速分布と比較して，これを大きく逸脱するような風速の鉛直分布は出現していない．しかしながら，Fig.9のRear Viewに注目すると，C点では速度ベクトルの鉛直分布が大きくねじれていることが分かる．一方，D点では風車ブレードの上端から地表面において速度ベクトルが大きく傾いている．

－ 225 －

速度ベクトルの鉛直分布に見られるねじれや傾きは，地形乱流の3次元構造に起因している．結果として，風車支柱とブレードを結合する根元付近に付加的な荷重が作用し，ボルトの金属疲労をより進行させたと推測できる．紙面の都合上ここでは示さなかったが，べき乗則(N値=7)に従う流入風速分布を大きく逸脱し，風車ハブ中心からブレード下端に大きな速度シアーが発生する時間帯も複数存在した．
５．おわりに
本研究では，太鼓山風力発電所を対象にLES乱流モデルによる非定常・数値風況診断を実施した．得られたシミュレーション結果に基づいて，風車周辺に発生している「地形乱流」の影響に注目して，数値風況面からの一考察を行った．西南西の風の場合，ナセル落下事故の起きた3号機を含むすべての風車位置で，べき乗則(N値=7)に従う風速分布から逸脱した速度シアーが存在すること，また大きな風向変化が風車のヨー方向で頻繁に発生しているという知見を得た．特に，主流方向(x)の比較的大きな標準偏差に加えて，鉛直方向(z)へも大きな値をとっていること，さらに主流直交方向(y)の標準偏差がx方向成分の値とほぼ同程度の値を有していることも認められた．
以上から，地形乱流の影響により励振力が増大し，結果として風車支柱とブレードを結合する根元付近に付加的な荷重が発生し，ボルトの金属疲労がより進行したのではな
いかと推測される．

＜謝辞＞
本研究は「公益財団法人JFE21世紀財団」の研究助成により 実施されたものであります．ここに記して感謝の意を表します．

－ 226 －