苫前ウィンビラ発電所における複雑地形上での地形及び風車のウェイクの影響を考慮した風況シミュレーション
Simulation of the wind in consideration of wake flow and complex terrain
at Tomamae Winvilla Wind park
○太田俊広
北谷匠子
露木和彦

(CRC ソリューションズ)
谷川亮一(CRC ソリューションズ)(電源開発)
兵藤浩晃 (電源開発)(CRC ソリューションズ)

Toshihiro Oota, CRC Solutions Corp. , 2-7-5, Minamisuna, Koto-Ku, Tokyo
Ryoichi Tanikawa, CRC Solutions Corp.
Shouko Kitaya, CRC Solutions Corp.
Kazuhiko Tsuyuki, ELECTRIC POWER DEVELOPMENT CO.,LTD.
Hiroaki Hyoudou, ELECTRIC POWER DEVELOPMENT CO.,LTD.
When the large scale wind power site is built, it can't get the expected amount of power
generation in many case with an influence of the wind turbine wake by the arrangementof the wind turbine. Then, the wind model which can calculate the wind which took theinfluence of wake into consideration in complex terrain was developed. Wake formation
is modeled by the momentum theory, and included into the k-ε model, and thedecrement of the velocity of the wind is found. As a result of calculating the influence ofwake at Tomamae Winvilla Wind park by using this model, in the part of wake, About10%-20% of the velocity of the wind decreased.
Keywords: Wind turbine, Complex terrain, Wake, k-ε model, Simulation

１．はじめに
２．基礎方程式
ウィンドファームのような大規模風力発電所における効ｋ‐ε型乱流モデルで用いた基礎方程式 1）を以下に示す。
率的な発電機の配列については、過去の調査・研究によっDu i
1 ∂P∂=−+uuて、いくつかの指標が存在する。海外での過去の研究結果
Dtρ ∂x i ∂x j i jから、一般に風車間の間隔は、風上側に 7〜12D、風向と直角に 1.5〜3.0D とされている。また、風速が弱く、風況 ∂u i ∂u j  2 − kδ ij+の乱れが少ないほど後列の発電ロスが大きいとされており、 u i u j = ν t  3∂x∂xi  j
風向偏差からウェイク領域を算出する方法が用いられている。しかし、これらはアメリカ、ヨーロッパなど、比較的k2ν=Ctµ
平坦で広い面積を確保できる場合に適用されるものであり、ε
日本の風力発電所のように十分な面積を得ることが難しく、Dk
狭く起伏の大きい地形に風車を建設する場合には、その指= Pk − ε + D k
Dt
標は必ずしも適用することができない。仮に適用できてもDε ε
予想通りの発電出力が得られない場合が多い。従って、よ= (C P − Cε 2ε ) + Dε
り狭い面積により多くの風車を効果的に配置することで、Dt k ε 1 k
効率的に発電出力を得る工夫が必要となる。それには今ま21  ∂u i ∂u j 
での指標とは異なる、より実情に沿った新しい手法が必要
Pk = ν t+となる。
2  ∂x j ∂x i 
そこで複雑地形でウェイクの影響を考慮した風況計算を行えるように k‐ε型乱流モデルを用いた風況モデルを開∂  ν t ∂k D=k発した。これは風車によるウェイク生成を運動量理論によ∂x j  σ k ∂x j りモデル化して風速の減衰率を求め、これを k‐εモデルに組み込んだものである。このモデルを用いて山岳地形で∂  ν t ∂ε Dε =ある苫前ウィンビラ発電所での実際の風車配置に対してウ∂x j  σ ε ∂x j 
ェイクの影響を計算したのでその結果を報告する。

( )実際にはこれを一般化座標で差分法で離散化して計算した。

３．風車モデルの組み込み
k−εモデルにウェイクの影響を取り入れ、複雑地形上での影響を計算するモデルを作成する。そのため風車によるウェイク生成をモデル化して風速の減衰率を求め、これを k−εモデルに組み込むこととした。モデル化の方法として観測データを用いて減衰率を求める方法と理論式から求める方法が考えられるが、今回のモデルでは理論式から求める方法を用いた。用いる理論式は発電出力と風車前面の風速から運動量理論を用いて風車直後の風速減衰率を求めるものである 2)。風車を吹き抜ける風のエネルギーのうち、発電された分だけ風のエネルギーが減るとみなし、減ったエネルギーの分だけ風車直後での風速が減衰するのを求めるものである。もとめた減衰率にしたがってモデル内部の風車位置における格子点の風速を減衰させることによりウェイクの影響を組み込んだ計算を行うことができる。
用いた理論式について具体的に説明する。風が風車を吹き抜ける際に連続の式と運動量の釣り合いを考慮すると、風車にかかる力 T は以下のようになる。
T=ρAV(U-V 1)
ρ:空気密度
A:風車面積V:風車直後の風速
U:風車前方の風速
V 1:風車から離れた位置でのウェイク風速（３−１）

また、ベルヌイの式からは
T=A(pu-pd)=0.5ρA(U2-V12)
pu:風車前方の圧力
pd:風車後方の圧力（３−２）となり、（３−１）式と（３−２）式より V と V 1 は以下のようになる。
V=0.5(U+V 1)（３−３）

すると発電量Ｐは（３−２）式と（３−３）式から
P=TV=0.5ρA(U2-V12)V=2ρA(U-V)V 2（３−４）となる。したがってある時点での風車前方の風速と発電出力がわかれば（３−４）式を用いて風車直後での風速を求めることができ、風速の減衰率を知ることができる。
この方法を用いたときの具体的な計算手順は以下のようになる。
① 風車の前面(風上側１メッシュの点)での風速と発電性能曲線を用いて発電出力を求める。
② 風車前面の風速と発電出力から運動量理論により風車背後での風速の減衰率を求める。
③ k−εモデルの風車直後の風速を②で求めた風速に設定し、次の時間ステップに進む。
実際には計算領域内で風車位置に相当する格子点に対して上記の手順を行うことになる。対象となる格子点は風車のローターの大きさと計算メッシュの大きさにもよるが、ローターの占める面積に対応する複数の格子点について行う。

４．計算結果
今回のモデルは苫前ウィンビラ発電所を対象として計算を行った。苫前ウィンビラ発電所は北海道苫前町の海岸近くにあり、出力 1650kw が１４基、1500kw が５基の計１９基の風車が設置された大規模な風力発電所である。計算領域は東西５km、南北５km、鉛直 1ｋｍであり、水平格子間隔は 50ｍ、鉛直格子間隔は地表付近を細かくした不等間隔、格子数は１０１×１０１×２０である。苫前ウィンビラ発電所の風車はローター直径が約 66ｍであり、水平方向のメッシュ間隔を 50ｍとしているので、風車 1 機あたり水平方向で 2 点、鉛直方向はそれよりも細かいため 3〜4点の格子点に対して３節で述べた風車モデルを適用した。
風車の仕様を表１に示す。
table 1 Wind turbine technical data

Nominal output
Units
Nominal wind speed
Cut inCut out
Rotor diameter
Number of blades
Hub height

1650kW
1417m/s
4m/s
25m/s
66m
360.4m

1500kW
513m/s
2.5m/s
25m/s
66m
360m

次に、計算に用いた地形および風車位置を図１に示す。
苫前ウィンビラ発電所の西側には別の大規模な風力発電所があるが、これも計算領域内に入るので計算対象とした。
ただしこの発電所の影響をあまり受けないように北西風の日に対して計算を行うこととし、日時は２００１年１２月
１７日午前 0 時の計算を行った。この日は強い北西風が吹いているケースである。計算の初期値は一定風向・風速とはせず、まず CRC ソリューションズで開発された局地気象モデル LOCALS 3)の 3km メッシュでの計算を行い、その結果を 50ｍメッシュの k−εモデルにネスティングして初期値とした。
以上の条件で計算した結果について、地上高 60ｍでの計算結果を風車無しの場合の風速を図２に、風車ありの場合を図３に示す。図 2 では、海上で１５m/s 程度の風が陸に向かって吹いており、地形の影響を受けて谷筋では風が弱く、海岸近くや標高の高いところでは風が強くなっているのが分かる。また、図３に示す風車ありの場合では基本的な風況は図 2 と同じような傾向を示しながら、風車の風下側では風速の小さい領域が現れており、風車によるウェイクが再現できていることが分かる。また、全体を見ても振動などの計算異常は現れておらず、流れ場全体に矛盾を引き起こすことなく風車モデルを組み込むことができている事がわかる。ウェイクを見やすくするために、風車ありの場合と無しの場合の風速の差をとったものを図４に示す。
図４を見ると風車のウェイクにより風車風下側で風速が減少している様子がとらえられている。また、風上側風車よりも風下側風車のウェイクのほうがより風速が下がっている。また、ウェイクの部分が途切れているところがあるが、これは実際に途切れ途切れになっているのではなく、図４
では地上高６０ｍの断面をとっているため、谷の部分ではウェイクよりも下側を表示しているためである。つまり図５の鉛直断面図に示すようにウェイクは地形に沿って上下にうねるのではなく、山の部分では地形に沿うが、谷の部分では水平に近い状態で形成されることが分かる。
これらの計算結果から、ウェイクの影響を考慮した気流計算が行われていることがわかり、k−εモデルへの風車モデルの組み込みは実現できているといえる。

[m]
fig.1 topographical map with wind turbine position

fig.3 contour of windspeed(with wind turbine)[m/s]
fig.2 contour of windspeed(no wind turbine)

[m/s]

[m/s]
fig.4 difference of windspeed between fig2 and fig3

fig.5 vertical section of A-A lineAA

５．まとめk−εモデルに風車モデルを組み込むことで風車ウェイクの計算機能の開発を行った。計算の結果、ウェイクの水平分布や鉛直断面を知ることが出来た。また、今後は観測データから統計手法により減衰率を求めるなど、より現実に近い減衰率を求めることや、観測データとの比較による精度検証などにより精度の向上をおこなう予定である。

参考文献
1）小林敏雄ら：数値流体力学シリーズ 乱流解析：東京大学出版会、1995
2）Spera,D.A.: WIND TURBINE TECHNOLOGY,ASME PRESS,1995
3）Kikuchi et al., 1981: Numerical Study on the Effectsof Mountains on the Land and Sea Breeze
Circulation in the Kanto District : Journal of theMeteorological Society of Japan Vol. 59, No. 5