白滝山ウインドファームの風車ブレード損傷事故の原因解明 - コンピュータシミュレーションによるアプローチ - 内田孝紀 ( 九州大学応用力学研究所 ) 92-583-7776, takanori@riam.kyushu-u.ac.jp 丸山敬, 石川裕彦 ( 京都大学防災研究所 ) 座古勝 ( 大阪大学大学院 ) 出口啓 ( 株式会社きんでん ) 1. はじめに 29 年 1 月 7 日から 8 日にかけて本州南部を台風 918 号 (MELOR) が通過した際, 山口県下関市豊北町の白滝山周辺の尾根上では, 非常に大きな乱れを伴う強風が観測された. この強風により, 株式会社きんでんが所有する白滝山ウインドファームでは, 風車のブレード破損事故が発生した ( 図 1, 図 2 を参照 ). 本研究ではその事故原因を解明するため, まず事故当時の気流性状をメソスケール気象モデル WRF-ARW 1) および LES(Large-Eddy Simulation) 乱流モデルに基づいた RIAM-COMPACT 2-4) を併用して詳細に再現した 5, 6). 次に, 得られた LES の時間平均場の結果を境界条件として, 風車ブレード近傍の気流解析を RANS(Reynolds Averaged Navier-Stokes equation) 乱流モデルを用いて別途行い, ブレードに発生する風圧力を求めた. 最後に,RANS 解析結果 を境界条件として有限要素法 (FEM,Finite Element Method) に基づいて, ブレードに作用する応力を調査した. 本研究の解析フローを表 1 に示す. 1. 気象モデル WRF-ARW によるウインドファームを含む広域スケールの気流解析 2. RIAM-COMPACT (LES 乱流モデル ) によるウインドファーム周辺の局所域スケールの気流解析 3. RANS 乱流モデルによる風車ブレード周辺の気流解析 4. 有限要素法 (FEM) による風車ブレードの応力解析 表 1 本研究における解析のフロー 数百 km 数 m 空間スケール風車の障害発生回数 メソスケール気象モデル WRF-ARW の解析時間帯 (29/1/7 の 9 時 ~24 時 ) 損傷部分 台風 18 号の通過径路 損傷部分の拡大図 図 1 29 年台風 18 号 (MELOR) の通過径路, 風車の障害発生回数の時刻歴, 風車ブレードの損傷状況
(a) 白滝山ウインドファームの位置図 GE 社製 2,5kW (b) 風車スペック図 2 白滝山ウインドファームの位置図および風車スペック 2. メソスケール気象モデルによる気流場の再現ウインドファームを含む広域スケールの気流場は, 米国大気科学研究センター (NCAR) を中心に開発された WRF-ARW(http://wrf-model.org/index.php) を用いて再現した. 図 3 に示すように, 計算領域は 4 段階のネスティングを用いた.4 種類の解析領域の高度, 鉛直空間解像度, 格子数は全て同じである. 高度は 14.5km, 空間解像度は不等間隔で格子数は 28 層である. 最も狭い領域 ( 図 3(d) の領域 4) の水平空間解像度は 1m である. 最も外側の領域 ( 図 3(a) の領域 1) は, 台風を含む広域スケールの気流場を再現できるように設定した. 地形標高データは, 宇宙航空研究開発機構 (JAXA) が作成した ASTERGDEM(http://www. ersdac.or.jp/gdem/j/index.html) を用いた. 境界条件に関して, 最も外側の領域 ( 図 3(a) の領域 1) には, 気象庁のメソ客観解析データを 1 日 8 回補間して与えた.29 年 1 月 7 日の午前 9 時に計算を開始し, 同日 24 時までの時間積分を実行した. 解析領域間の各種情報のやり取りは, 外側の格子情報を内側の格子に与えるだけの 1-way ネスティングで行った.
1km油谷アメダス観測点 83km高さ 14.5km, 28 層 1白滝山ウインドファーム (a) 領域 1(31 31 点 ) (b) 領域 2(1 1 点 ) (c) 領域 3(1 1 点 ) (d) 領域 4(13 13 点 ) 水平格子解像度 :2,7m 9m 3m 1m 図 3 メソスケール気象モデルWRF-ARWによる気象場の再現計算で用いた解析領域と4 重ネスティング 事故発生時間帯最大 (m/s) 25 平均 : 観測値 2 最大 : 観測値 : 領域 2の計算値 15 : 領域 3の計算値 1 平均 : 観測値最大 : 観測値 5 : 領域 2の計算値 3 6 9 12 15 18 21 時間 (h) : 領域 3 の計算値 図 4 29 年 1 月 7 日の油谷アメダス観測点における の変化, 観測値と計算値 (WRF-ARW) の比較 N A 3km (a) 午後 2 時頃 (b) 午後 1 時頃 観測値 図 5 29 年 1 月 7 日午後の分布, 計算値 (WRF-ARW, 領域 3) 計算値 45 (m/s) [m/s] 事故発生時間帯 事故発生時間帯 事故発生時間帯 35 25 15 降水を伴う気象擾乱 降水を伴う気象擾乱 降水を伴う気象擾乱 N W S E N 5 12 15 18 21 2412 15 18 21 2412 15 18 21 24 風車 1 号機 time 時間 (h) 風車 7 号機 time 時間 (h) 風車 17 号機 time 時間 (h) 図 6 風車ナセル上の に関する観測値 (1 分平均値 ) と計算値 (WRF-ARW, 領域 3) との比較, 観測高さは地上高約 85m, 計算結果は地上高約 87m,29 年 1 月 7 日
WRF-ARW の領域 3 白滝山ウインドファーム 流入風 ( 北東の風 ) 21km 8km z 2km x y N 白滝山ウインドファーム (a) 解析領域 ( 速度ベクトルは間引いて表示 ) (b) 地形性状図 8 LES(RIAM-COMPACT ) で用いた図 7 LES(RIAM-COMPACT ) で用いた解析領域と地形性状計算格子,211(x) 81(y) 41(z) 点 29 年 1 月 7 日, ウインドファーム周辺では午後からが増加し始め, ウインドファームに最も近い油谷アメダス観測点 ( 図 3(c) を参照 ) では, 図 4 に示すように,21 時頃に最大の発生が確認され ( 図中に表示 ), の計算値はアメダス平均と最大の間にある. 一方, 変化はアメダス観測値とほぼ一致している. 風車のブレード破損事故が発生し始める 1 月 7 日の午後の WRF-ARW の領域 3 における分布 ( 図 5) を吟味すると, 午後 2 時頃 ( 図 5(a)) では, 地上付近において北東の風が吹いており, 上空に行くに従っては増加している. 但し, が最大値を取った後 ( 図中の矢印 A), 上空 3km 付近でははゼロに近づいている. 図には表示していないが,3km より上空ではが西向きに転じ, 高度とともに速くなっている. これに対し, 午後 1 時頃 ( 図 5(b)) では, 地上付近のは午後 2 時頃よりも北に向き, 上空 3km においても北東の風が吹き, その大きさも図 5(a) と比較して増加している. 風車のナセル ( 地上高約 85m) に取り付けた 計の記録 ( 図 6 の実線 ) を観察すると, が北東よりも北寄りに振れた 19 時以降, ともに変動が大きくなっている. この変動は, 降水を伴う気象的な擾乱の通過に対応している. これに対応する計算結果 ( 図 6 の点線 ) も はほぼ同じような変化を示している. しかしながら, その値は計算値と観測値では少し違いがある. に関しては, 計算値の方が観測値よりも高めになっている. また, 計算結果の変動は, 観測値に比べてゆっくりとしており, 実際のに見られるような高周波成分の乱れを再現できていない. よって, 高解像度の計算格子を用いて周囲の地形起伏を詳細に再現した上で, 気流の時間変動を再現可能な LES 乱流モデルを用いた解析を行うことにした. 3.LES 乱流モデルによる気流場の再現図 7 に,LES 乱流モデル (RIAM-COMPACT ) の解析範囲を示す. メソスケール気象モデル WRF-ARW の解析領域 3 の内側でウインドファームを含み, 北東方向を主軸 (x 軸 ) とする長方形の領域を設定した. 図 8 に示すように,x 軸方向に 21km,y 軸方向に 8km,z 軸方向に 2km の空間を有し, 格子数は x 軸方向に 211 点,y 軸方向に 81 点,z 軸方向に 41 点とした. 格子幅はウインドファーム周辺で細かく設定し, 各方向の最小格子幅は x 軸方向では 35m,y 軸方向では 38m,z 軸 (m/s) 681m LES の基準スケール 図 9 LES(RIAM-COMPACT ) で用いた流入プロファイル (WRF-ARW の結果から作成 ) 6 5 4 3 2 水平成分鉛直成分 降水を伴う気象擾乱の通過 < NE < NNE 1 < N 事故発生時間帯 -1 9 12 15 18 21 24時間 (h) 図 1 高度 681mにおける の変化, 29 年 1 月 7 日,WRF-ARWの結果, 領域 3
ウインドファームの上流の地形から形成される変化を伴う空気塊の列 白滝山ウインドファーム 図 11 LES(RIAM-COMPACT ) の計算結果の一例, 北東の風の場合, 風車ハブ高さ ( 地上約 85m) における主流方向 (x) の速度成分コンター図 方向では 2m である. 計算領域の外側に向かって格子幅は徐々に大きくなっている. LES(RIAM-COMPACT ) の計算に関して, 上流側の流入境界面には, 図 9 に示すメソスケール気象モデル WRF-ARW による計算結果を内挿して与えた. また, 流入境界面における高度 681m( 解析領域中の最高標高を基準空間スケールとしている ) における の変化を図 1(WRF-ARW の結果 ) に示す. 風車ブレードの破損事故が発生した 1 月 7 日 14 時から 22 時にかけて, は上空で北東から北に向って変化し, 平均は 25~27m/s 程度の値を示す. ともに,19 時以降にその変動が大きくなっている. 既に述べたように, これらの変動は降水を伴う気象擾乱の通過に対応している. 以下では, 風車ブレードの破損事故が発生した時刻ので計算を行い, ウインドファーム周辺の気流性状の変化を調べた結果を示す. 図 11 には,LES(RIAM-COMPACT ) の計算結果の一例を示す. これは北東の風の場合で, 風車ハブ高さ ( 地上約 85m) における主流方向 (x) の速度成分コンター図である. 全体的な傾向は以下の通りである. (m/s) x y ウインドファームの風上側で地形の影響により発生した変動の大きい領域が, 主流方向 (x) に移動し, ウインドファームに断続的に流れ込む. 結果として, ウインドファーム周辺で大きな変動が生じる. この断続的な変化の大きい領域の発生場所と流下方向はにより変化するので, 各風車ナセル位置における変化もにより大きく異なる. 風車ナセルに取り付けられた 計の観測記録を精査すると, タワーの振動やブレードの被害が見られた時間帯の の変化は大きく, 同時に乱れの強さ (σ) も平均の 5% 程度に達することが分った ( 図 12(a)). これに対応する LES (RIAM-COMPACT ) の計算結果 ( 図 12(b)) に注目する. 変動の高周波成分が再現できていないので, 乱れの強さは観測値よりも小さくなっている. しかしながら, その値は 3% を超え, の変化は観測結果と類似な傾向が明確に再現されている. 一般的に平均が小さい場合には, 発電効率を高めるために大きな風力が発生するようにブレードのピッチ角を制御する ( 後述する図 13 および図 14 を参照 ). この際, ピッチ角が追従できないような急激な気流の変化が生じると, ブレードに損傷を及ぼすほどの大きな風荷重が加わる. このような気流変化が最も危険である. 図 12 に示されるような気流変化がこれに対応している. このとき, 乱れの強さおよび突風率は非常に大きな値を示す. 以上のように,LES(RIAM-COMPACT ) によって得られたウインドファーム周辺の気流性状の計算結果は, 上流側に位置する地形の影響を受けた の急激な変動を定性的に再現可能であることが示された. また, 乱れの強さは LES による計算値は観測値と比較して約 1% から 2% 程度小さくなることも明らかになった. 4.RANS 乱流モデルによる風車ブレード周辺の気流解析と有限要素法 (FEM) による応力解析前節の風況シミュレーションの結果, 風車は激しい 変動の中で運転され, さらに定格出力を発生させるよう制御されていることが確認された. < NE < N < NW 時間 (s) (a) 観測値, 平均 8.4m/s, 乱れ強さ 54% (b)les 計算値, 平均 12.6m/s, 乱れ強さ 35% 図 12 風車ハブ高さ ( 地上高約 85m) における 1 分間の の時系列データの例
内部桁 (TE 部 ) 後縁 :Trailing-Edge(TE) 2 内部桁 ( 中央部 ) 内部桁 (TE 部 ) 左図の A 断面 内部桁 ( 中央部 ) 前縁 :Leading-Edge(LE) 1 抗力側 :Pressure-Side(PS) 4 2 1 A 揚力側 :Suction-Side(SS) 3 図 13 風車ブレードの全体図 ( 左 ) および断面図 ( 右 ) 3 Flow Flow 偏差 + 偏差 ± 偏差 - (a) ピッチ角 : ゼロ度 (b) ピッチ角 :3 度 翼回転方向 図 14 ブレードピッチ角 ( 左 ) および偏差 ( 右 ) の定義 内部桁 ( 中央部 ) 内部桁 (TE 部 ) 図 15 ブレード構造の写真 ( 左 ) および FEM 解析モデル ( 右 ) 表 2 RANS および FEM 解析条件 解析条件 (m/s) 偏差 ( ) 回転速度 (rpm) ブレードピッチ角度 ( ) 模擬状況 1 3 16.5 暴風 - ピッチ非追従 2 3-4 16.5 暴風 - ピッチ ヨー非追従 この様な状況では, 風車のピッチ角制御が変動に追従しない状態や, ヨー制御も変動に追従できない状態が発生することが推測された. そこで, 風車の運転記録から実際に発生したと思われる条件を, 偏差 回転速度 ピッチ角度で分類し, その条件ごとにブレード周りの流動解析 (RANS 解析 ) を行い, ブレードに発生する風圧力を求め, ブレードが破損する条件を有限要素法 (FEM) 解析により調査した ( 図 13, 図 14, 図 15, 表 2 を参照 ).
大 応力 : 大 隔壁から 3,mm 前縁 :Leading-Edge(LE) の損傷 応力 小 図 16 ブレード軸方向応力図 ( 表 2 の解析条件 2) と実際の写真 解析結果によると,LE( 前縁 :Leading-Edge) 接合部に高い応力が発生していることが判明し, その発生場所も損傷事故が発生した場所と同位置であることが示された. ブレードがピッチ角ゼロ度の状態において, ピッチ制御が追従できず,3m/s の突風を受ける解析条件 1 や, 偏差が加わる解析条件 2( 図 16 を参照 ) では, 応力の増大が確認された. LE 接合部に使用される接着剤の静的強度はかなり大きいものの, 同様の荷重を繰り返し受けると低サイクル疲労で破壊する可能性が高い. 今回の 918 号台風の通過では, 強風条件下で数多くの偏差が観測されている. よって, 低サイクル疲労により, 損傷 ( 亀裂など ) が LE 接合部に発生し, その後の運転でその損傷が進展したものと推測された. LE 接合部に高い応力が発生し, 接着剤に損傷が発生することが判明したため, 風車ブレードには幾つかの事故防止対策を講じた. 具体的には, 低サイクル疲労に対応するため,1LE 接合部のオーバーラミネート ( 図 17),2LE 部の内部補強リブの増設,3 内部桁部の補強などを行った. これらの効果を検証するため, 補強したモデルに対して, 表 2 に示した解析条件 2 の下, 再び FEM 解析を行った. その結果を図 18 に示す.LE 部に発生していた応力は約 3 割に軽減され,918 号台風と同等の気象条件下において損傷が発生しないことを確認した. 5. おわりに 29 年台風 18 号は本州の南を北東進したため, 白滝山ウインドファーム一帯は台風中心から遠く離れ, しかも一般に影響が小さいとされる台風進行方向左側に位置していた. しかし, 北側にあった高気圧の影響で気圧勾配が強まり, これによる台風北側の東寄りの強風がウインドファームに吹き込んでいた. 白滝山における最多は北である. 今回のように強い風が, 東および北東からウインドファームに吹 LE 部の補強 揚力側 :Suction-Side(SS) 抗力側 :Pressure-Side(PS) 図 17 LE 部補強 ( オーバーラミネート ) 図 17 に比べて, 応力は 3 割減少 図 18 補強後のブレード軸方向応力図 ( 表 2 の解析条件 2) き込むことは稀な事象である.LES 乱流モデル (RIAM-COMPACT ) による結果から, 北東方向から強い風が流入する場合は, ウインドファームの風上に位置する地形の影響により, の急激な変化を伴う風が風車に流入することが判明した. このような気象条件では, ピッチ制御が変動に追従できない状態が起こり, ブレードには設計条件 7, 8) を超える風圧力が発生する可能性がある.FEM 解析の結果からも, ピッチ制御が変動に追従できない場合には,LE 接合部に大きな応力が発生し, その発生場所も実際の破損箇所と同じ位置であることが明らかとなった. 警報履歴から, 非常に多くの強風下での偏差が観測されており, このような応力が繰り返し LE 接合部に作用したため, 低サイクル疲労により破壊して損傷が入り, その後の運転で進展したものと推定される. 今回の一連の解析結果を総合的に精査 吟味す
ることにより, 我々のグループでは, 風車ブレードに大きな風荷重を加えることなく, 発電量を稼ぐための風車の安全運転制御上の指針を見出し, それを実際の運転制御に活用することを検討している. すなわち, 風車の稼働率を低下させず, 高い水準に保つ方策である. 例えば, 風車ブレードのピッチ角制御が, の変化に追従できなくなるような急激で, かつ大きな の変化が生じるような気象条件を各ナセル位置で事前に予測しておく. そのような気流条件を満たす場合には, 発電を中止するなどの対応を取る ( 図 19 を参照 ). 先に述べたように, このような危険な変動を伴う気流は, 具体的には図 12 に示すような大きなスケールの変動を有し, 急激で, 大きな 変動を伴うものを意味する. 図 19 運転制御の概念図 制御を行う範囲地形の影響による の大きな変動が発生する また, 風車ブレードには幾つかの事故防止対策を講じた. 具体的には, 低サイクル疲労に対応するため,1LE 接合部のオーバーラミネート,2LE 部の内部補強リブの増設,3 内部桁部の補強などを行った. 補強後のモデルを用いた FEM 解析では,LE 部に発 生していた応力は約 3 割に軽減され,918 号台風と同等の気象条件下において損傷が発生しないことを確認した. 参考文献 1) Skamarock, W.C. et al., 28. A Description of the Advanced Research WRF Version 3, NCAR TECHNICAL NOTE, National Center for Atmospheric Research Boulder, Colorado, USA. (http://wrf-model.org/ index.php) 2) T. Uchida and Y. Ohya, Micro-siting Technique for Wind Turbine Generators by Using Large-Eddy Simulation, Journal of Wind Engineering & Industrial Aerodynamics, Vol.96, pp.2121-2138, 28 3) 内田孝紀, 大屋裕二, 李貫行, 風車立地点近傍に発生する地形乱流の高解像度 LES, 風力エネルギー協会誌,Vol.34, 通巻 93, pp.121-126,21 4) 内田孝紀, 大屋裕二,RIAM-COMPACT による数値風況予測の最前線 設計評価とウインドリスク ( 地形乱流 ) 診断, 風力エネルギー協会誌,Vol.34, 通巻 93,pp.69-74,21 5) 丸山敬, 内田孝紀, 石川裕彦, 出口啓, ウインドファーム周辺の変動場の数値解析, 第 24 回数値流体力学シンポジウム論文集, 現在印刷中,21 6) 丸山敬, 石川裕彦, 内田孝紀, 出口啓, メソスケールモデルと LES を用いたウインドファーム周辺の気流解析, 第 21 回風工学シンポジウム論文集, 現在印刷中,21 7) 内田孝紀, 丸山敬, 竹見哲也, 奥勇一郎, 大屋裕二, 李貫行 : 気象モデルと流体工学モデルを用いた風車設置地点における設計評価手法の提案, 風力エネルギー協会誌, Vol.34, 通巻 94,pp.118-124,21 8) 内田孝紀, 丸山敬, 大屋裕二, 連続的な変化を考慮した非定常数値風況予測による風車設置地点における設計評価手法の提案, 風力エネルギー協会誌, 現在印刷中,21