第 28 回数値流体力学シンポジウム
講演番号 B08-2

3 次元 CFD と BEM による風車解析の比較
Comparison between Three Dimensional CFD and BEM for Wind Turbine Analysis○ 喜屋武 拓麻, 琉球大院, 沖縄県中頭郡西原町字千原 1 番地, E-mail:can.tama259@gmail.com
鈴木 正己, 琉球大, 沖縄県中頭郡西原町字千原 1 番地, E-mail:m-suzuki@tec.u-ryukyu.ac.jp
天久 和正, 琉球大, 沖縄県中頭郡西原町字千原 1 番地, E-mail:kazumasa@tec.u-ryukyu.ac.jp
Takuma KYAN, University of Ryukyu, 1 Senbaru, Nishihara-cho. Nakagami-gun, Okinawa, JapanMasami SUZUKI, University of Ryukyu, 1 Senbaru, Nishihara-cho. Nakagami-gun, Okinawa, JapanKazumasa AMEKU, University of Ryukyu, 1 Senbaru, Nishihara-cho. Nakagami-gun, Okinawa, JapanThe Blade Element Momentum theory (BEM) had been used for a design and the performance prediction of the windturbine, while we have come to be able to perform the three-dimensional CFD analysis easily and rapidly by thedevelopment of the computer. In this study, analysis of the wind turbine performance by BEM and three-dimensional
CFD is performed. The results are investigated for the difference between BEM and CFD, the possibility as design tooland the prediction precision of BEM by comparing with an in-house CFD code.
１．緒言
これまで風車の設計や性能予測には翼素運動量理論（BEM）が広く用いられてきているが，近年のコンピュータの発展により 3次元の CFD 解析も容易かつ高速に行えるようになってきており，より詳細で高精度な結果が CFD に期待されている．また，BEM
は 2 次元翼の揚力係数と抗力係数を必要とするが，翼根から翼端にかけてレイノルズ数は変化し，多くのレイノルズ数に対するデータが必要となる．しかし，実験データとして多くのレイノルズ数に対する揚力係数や抗力係数は完備しておらず，現実的には代表的なレイノルズ数のみのデータを用いて解析せざるを得ない．
BEM は簡易的で容易というメリットはあるものの，上述のような課題も存在している．より詳細な設計や解析には，流れ場を把握することが重要であるが，BEM や実験では難しい点である．CFDは実験に比較し，一般的に結果の信頼性においては劣るものの，詳細な流れ場の取得が容易であり，CFD の高精度化にともない利用価値が高まり，広く利用されつつある．
本研究では風車性能の解析をBEM と3 次元CFD によって行い，実験結果と比較することによって BEM と CFD の差異，BEM の予測精度や設計ツールとしての可能性を検討する．と回転力を算出する方法であり，翼素の空力性能を予め実験等から取得する必要がある．運動量理論ではロータ平面を含む流管の運動量変化からスラストと回転力を推定する．

Fig. 1 Schematic view of experimental apparatus ofDelft University of Technology.
翼素運動量理論では両理論から導出されるスラスト，回転力を等しく置くことで風車を通り抜ける軸流速度𝑣𝑎 や風車後流に誘起される旋回速度𝑣𝑡 を求める．翼素理論のみではタービンへの流入風速とブレードの周速による幾何学的な迎え角がわかるのみである．実際の流入速度はタービンを避けて流れるため，タービンを通過する流速は十分上流の流入風速よりも遅くなる．
このため，運動量理論によりタービンを通過する流速を予測し，より妥当な迎え角を算出する．両理論を単独で用いる場合よりも高い性能予測が可能である．一方であくまで簡易的な計算手法であるが故の課題もある．最も大きな負荷（スラスト）となるべき極限状態となる回転円板では運動量理論のスラストは零となる．これを修正し現実的な値となるように，大きなスラストでは経験的な実験式と切り換えて用いられている．これは運動量理論の適用が不十分なために生じたものである．
BEM の計算結果は翼素理論で用いる翼性能データによって大きく左右される．本研究では 2 次元の数値流体解析により翼性能を算出し利用した．また，比較のため Sandia National Laboratoriesの予測データを用いた結果も併せて示す．

２．解析対象解析・比較の対象には Vermeer らによるデルフト工科大学での風洞実験(1)を選択した．Fig.1 に実験装置の概略を示す．吹出型風洞の前方に 2 枚翼風車が設置され，風速，回転数，トルク，スラストが計測されている．設置した風車は NACA0012 で構成され，ロータ直径 1.2[m]，翼弦長は一定で 0.08[m]，翼端から翼根にかけて捩じりがつけられている．ピッチ角𝜃 [deg]は無次元半径𝑟/𝑅 = 0.3~0.9 で 𝜃𝑡𝑖𝑝 + 4° か ら 𝜃𝑡𝑖𝑝 ま で 線 形 に 減 少 し ，𝑟/𝑅 = 0.9~1.0で翼端ピッチ角𝜃𝑡𝑖𝑝 に固定されている．本研究では𝜃𝑡𝑖𝑝 = 2°での解析結果を示す．実験は風速 5[m/s]で実施され，ロータ半径𝑅と風速𝑈を代表長さ，代表速度とするレイノルズ数
は𝑅𝑒 = 𝑈𝑅/𝜈 = 2 × 105 [-]である．ここで𝜈は空気の動粘度である．
３．翼素運動量理論
３．１．概要翼素運動量理論(2)は翼素理論と運動量理論を組み合わせてタービン性能を解析する方法である．風車の出力係数や推力係数の推定，翼弦長と捩じり角の分布の最適設計に広く用いられている．
翼素理論はブレード翼素の揚力係数と抗力係数からスラスト

３．２．翼素理論
翼素理論は揚力と抗力を用いて翼素に働くスラストと回転力（または接線力）を求める方法である．Fig.2 に翼素に働く流体力（揚力𝐿，抗力𝐷，軸力𝐹𝑎 ，接線力𝐹𝑡 ）とブレード間の平均風速𝑣𝑎 ，
1Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 B08-2
次に接線力𝑑𝐹𝑡 ，接線力係数𝐶𝑡 は𝑑𝐹𝑡 = 2𝜌(1 − 𝑎)𝑎𝑡 𝑉𝑎 𝑟Ω(2𝜋𝑟𝑑𝑟)(1.8)= 2𝜌𝑎𝑡 (1 − 𝑎)𝜆𝑟 𝑉𝑎2 (2𝜋𝑟𝑑𝑟)

周速𝑟𝛺，後流側旋回速度𝑣𝑡 = 2𝑎𝑡 𝑟Ω，相対速度𝑊，迎え角𝛼，ブレードピッチ角𝜃，流入角𝜙との関係を表している．

𝐶𝑡=𝑑𝐹𝑡1 2𝜌𝑉 (2𝜋𝑑𝑟)
2 𝑎(1.9)
2𝑎𝑡 (1 − 𝑎)𝑉𝑎2 𝜆𝑟 𝜌2𝜋𝑟𝑑𝑟
1 2𝜌𝑉 (2𝜋𝑟𝑑𝑟)
2 𝑎= 4𝑎𝑡 (1 − 𝑎)𝜆𝑟となる．ここで𝑎𝑡 は角誘導係数で，旋回速度の誘導係数を表し，次式に示すように翼前後の旋回速度差𝑣𝑡 の半分を動翼周速𝑟𝛺で除したものである．
𝑣𝑡(1.10)𝑎𝑡 =2𝑟Ω=Fig. 2 Fluid forces acting a blade element.
半径𝑟において微小区間𝑑𝑟の翼素に作用するスラスト𝑑𝐹𝑎 は1(1.1)
𝑑𝐹𝑎 = (𝐶𝐿 cos 𝜙 + 𝐶𝐷 sin 𝜙)𝜌𝑊 2 𝑐𝑁𝑑𝑟2となる．ここで𝑊は相対速度で𝑊sin𝜙 = (1 − 𝑎)𝑉𝑎 ，𝑐は翼弦長，𝑁は翼枚数である．スラストを無次元化したスラスト係数は
𝐶𝑎𝑑𝐹𝑎
1 2𝜌𝑉 (2𝜋𝑑𝑟)
2 𝑎1−𝑎 2=() (𝐶𝐿 cos 𝜙 + 𝐶𝐷 sin 𝜙)𝜎sin 𝜙

３．４．翼素運動量理論
翼素理論と運動量理論から得られるスラスト係数，式(1.2)と(1.7)を等値とすると
1−𝑎 2
𝐶𝑎 = () (𝐶𝐿 cos 𝜙 + 𝐶𝐷 sin 𝜙)𝜎sin 𝜙(1.11)
4𝑎(1 − 𝑎),(𝑎 ≤ 𝑎𝑐 )={4[1 − √𝐶𝑎,𝑎=1 ](1 − 𝑎) + 𝐶𝑎,𝑎=1 , (𝑎 > 𝑎𝑐 )=(1.2)

を得る．なお2 − √𝐶𝑎,𝑎=1𝑎𝑐 =2である．ただし，運動量理論によるスラスト係数は誘導係数が𝑎 = 0.5で最大となり𝑎 = 1で𝐶𝑎 = 0となるが，実際は𝑎 > 0.5でも𝑎とともに大きくなる．このため実験結果と合うように，𝑎 > 𝑎𝑐は一次式で近似している．𝐶𝑎,𝑎=1 は𝑎 = 1における𝐶𝑎 であり，本研究では𝐶𝑎,𝑎=1 = 1.816としている．接線力係数については式(1.4)と(1.9)より

2(1 + 𝑎𝑡 )𝜆𝑟=() (𝐶𝐿 cos 𝜙 + 𝐶𝐷 sin 𝜙)𝜎cos 𝜙で表される．ここで，𝜎 = 𝑐𝑁/(2𝜋𝑟)はソリディティである．また，接線力𝑑𝐹𝑡 は1(1.3)
𝑑𝐹𝑡 = (𝐶𝐿 sin 𝜙 − 𝐶𝐷 cos 𝜙)𝜌𝑊 2 𝑐𝑁𝑑𝑟2であり，接線力を無次元化した接線力係数𝐶𝑡 は次のようになる．
𝐶𝑡𝑑𝐹𝑡
1 2
𝜌𝑉 (2𝜋𝑑𝑟)2 𝑎1−𝑎 2=() (𝐶𝐿 sin 𝜙 − 𝐶𝐷 cos 𝜙)𝜎sin 𝜙
=

2𝐶𝑡
(1.4)

= 4𝑎𝑡 (1 − 𝑎)𝜆𝑟2= 4𝑎𝑡 𝜆𝑟 (1 + 𝑎𝑡 )𝜆𝑟 tan 𝜙0 ， ∵ tan 𝜙0 =(1 + 𝑎𝑡 )𝜆𝑟=() (𝐶𝐿 sin 𝜙 − 𝐶𝐷 cos 𝜙)𝜎cos 𝜙→𝑉 −𝑣𝑎
𝑉𝑎

𝑉 −𝑉𝑎,𝑜𝑢𝑡， 2𝑎 = 𝑎𝑉𝑎=𝐶𝐿𝐶𝐷tan 𝜙
4𝑎𝑡 = (1 + 𝑎𝑡 )(−)𝜎(1.12)
cos 𝜙 sin 𝜙 tan 𝜙0となる．翼端損失がないとき𝜙 = 𝜙0 である．
翼端損失を考慮するために，Prandtl により提案された翼端損失係数𝐹を導入する．
2(𝑁/2)(1 − (𝑟/𝑅))(1.13)
𝐹 = cos −1 [exp {−}]
𝜋(𝑟/𝑅) sin 𝜙0
ここでは翼端渦による誘導速度の影響として，翼端損失係数を迎え角に関連付けて定式化する．流入角𝜙は迎え角𝛼とピッチ角𝜃の和であり，これと誘導係数との関係は
𝛼 = (𝜙0 − 𝛼0 − 𝜃)𝐹 + 𝛼0(1.14)
𝜙 = 𝛼+𝜃で表す．ここで，𝛼0 は揚力係数𝐶𝐿 = 0となる迎え角である．

(1.6)

𝑑𝐹𝑎
1 2𝜌𝑉 (2𝜋𝑑𝑟)
2 𝑎2𝑎(1 − 𝑎)𝑉𝑎2 𝜌2𝜋𝑟𝑑𝑟
1 2𝜌𝑉 (2𝜋𝑟𝑑𝑟)
2 𝑎= 4𝑎(1 − 𝑎)=

1 + 𝑎𝑡𝜎cos2 𝜙∴なお，𝑣𝑎 は翼間の平均軸流速度，𝑣𝑎,𝑜𝑢𝑡 は風車から十分下流での軸流速度である．スラスト係数は𝐶𝑎(𝐶𝐿 sin 𝜙 − 𝐶𝐷 cos 𝜙)1−𝑎(1−𝑎𝑡 )𝜆𝑟= 4𝑎𝑡 tan 𝜙0３．３．運動量理論
運動量理論は翼前後の運動量変化からスラスト，トルク（接線力）を求める方法である．運動量理論によるスラストは(1.5)
𝑑𝐹𝑎 = 2𝜌𝑉𝑎 2 (1 − 𝑎)𝑎(2𝜋𝑟𝑑𝑟)
ここで𝑎は軸流速度の誘導係数を表し，軸流速度間で次の関係が成り立つ．
𝑎= 𝑎(1 + 𝑎𝑡 )𝜆𝑟=() (𝐶𝐿 sin 𝜙 − 𝐶𝐷 cos 𝜙)𝜎cos 𝜙(1.7)

2Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 B08-2
４．3 次元数値解析
数値流体解析は研究室にて開発を行っている 3 次元非圧縮性・粘性解析コードで実施した(3)．この計算コードは Collocated Gridによる有限体積法(4)で構成され，連続の式を満たすため SIMPLEアルゴリズム(5)を用いて圧力補正を行っている．対流行の近似には QUICK スキーム(6)を用い，その他の項は 2 次精度の中心差分で離散化されている．乱流モデルは Kato-Launder k-𝜀 モデルにLaunder-Sharma の低レイノルズ数型 k-𝜀モデル(7)を導入している．
Fig.3 にブレード周りの計算格子を示す．下面は周期境界条件が適用され，1 翼すなわち半球領域で計算が行われている．この半球の半径はロータ半径の 5 倍としている．なお，Fig.1 に示す通り実験に用いられた吹出型風洞の内径はロータ直径の 2 倍程度となっており，十分な大きさの計算領域が設けられている．

５．結果・考察
５．１．解析結果と実験の比較
Fig.4 に周速比ごとのスラスト係数を示す．3 次元 CFD およびBEM のいずれも実験結果とよく一致している．式(1.1)からもわかる通り，スラストは揚力とほぼ一致し，抗力の影響は小さい．単独翼の揚力係数はポテンシャル計算でも実験値を 1 割ほど上回る程度であるため，数値解析による高精度な推定が可能であり，スラストも高精度な予測が同様に期待できる．
一方で，式(1.3)に示すように接線力は揚力とともに抗力の影響も強く受け，揚力による力と抗力の差となるため揚力に比して 1割以下の小さな力となる．このため数値解析，実験のいずれにおいても接線力および出力係数は大きな誤差を生じやすい．このように予測の難しい出力係数ではあるが，Fig.5 に示すように周速比
𝜆 > 6では実験結果と良好な一致が得られている．実験結果は周速比𝜆 = 7.5にて最大出力係数が現れ，𝜆 < 6では失速域に入り，急激に出力が低下している．CFD では失速角が実際よりも大きく（すなわち失速しづらく）なる傾向があり，3 次元 CFD および 2次元CFD を基にしたBEM がこの領域で実験値を上回る結果を示した原因と考えられる．

ここに用いた O-O 型格子は格子を，翼近傍に集中して配置でき，無駄の少ない適切な格子配置が得られる．格子は代数的な格子生成法である Transfinite 補間法(8)を用いて生成している．格子点数は翼周りに 130 点，翼幅方向に 29 点，翼面法線方向に 52 点，総格子点数 196,040 点になっている．ロータ半径を単位長さとした翼面法線方向の最小格子幅は5 × 10−5 であり，𝑦 + はすべて 0.65 以下に収まっている．

５．２．3 次元 CFD と BEM による解析結果の詳細比較
前節に示したように，計算によるスラスト係数や出力係数などのマクロな特性は実験結果と良く一致しており，3 次元 CFD 計算による詳細な流れ場も高精度に予測されていると推察できる．このため，3 次元 CFD の結果を基準に BEM 解析の評価を行う．
Fig.6 は周速比𝜆 = 9.0におけるロータ平面を通り抜ける平均軸流速度の半径方向分布を示している．図より，BEM では翼端と翼根において CFD よりも急激な変化となっている．BEM で用いた翼端損失係数（式(1.13)）は翼枚数，半径，流入角のみを因子とした単純な表現であり，アスペクト比や翼面負荷の分布などの影響が十分に反映されていない．このように，BEM では翼幅方向の影響が十分に評価できるとは言えず，簡易的な解析方法であることが見て取れる．

Fig. 4 Thrust coefficients.

Fig. 6 Span-wise distribution of axial induction factor.

Fig. 3 Computational grid of a half domain.

Fig. 5 Power coefficients.

Fig. 7 Span-wise distribution of thrust coefficients.
3Copyright © 2014 by JSFM

第 28 回数値流体力学シンポジウム
講演番号 B08-2
６．結言
BEM は計算負荷が少なく解析方法が容易であるため，現在でも設計や性能予測に広く利用されている．ただし，単純なモデルに基づく解析であるため，その結果の信頼性や適用限度に関して多くの検証が必要である．本研究では BEM に関して，CFD や実験結果との比較を通した検討結果を示している．以下に本研究で得られた結果のまとめを列挙する．
１） スラスト係数や出力係数といったマクロな量の性能予測において，BEM は 3 次元 CFD に劣らない高い予測精度が得られている．ただし，BEM による解析結果は揚力係数，抗力係数のデータに大きく影響されるため，正確なデータの入手，選定が課題である．
２） 翼幅方向の軸流速度やスラスト係数，出力係数の分布など，BEM では翼幅方向の影響について十分な推定ができるとは言えない．CFD ではこれらについても予測できる
他，詳細な流れ場が得られるなど非常に有力な設計，性能予測ツールである．
３） BEM は揚力係数，抗力係数のデータが存在すれば非常に短時間で解析でき，手軽に活用できる点が利点である．3次元 CFD では多くの計算時間を必要とするが，コンピュータの性能向上や適切な格子配置による格子点数の削減により，実用的な範囲に計算時間も収まってきており，ますます期待できる解析ツールである．
なお，本研究の一部は科学研究費補助金（21500095）の助成により実施されたことを記して謝意を表する．

Fig. 8 Span-wise distribution of power coefficients.
Fig.7 はスラスト係数の半径方向分布を示している．3 次元 CFDと BEM は良好な一致を示している．多少の違いは軸流速度に大きな差が生じていた翼根および翼端付近に表れており，BEM の方が翼端付近で急激な減少を示している．
Fig.8 は出力係数の分布である．スラスト係数の分布と比較して半径位置𝑟/𝑅 = 0.7~0.8で顕著な違いが表れ，BEMの結果はCFDと比較して低い値を示し，分布も緩やかな変化となっている．ただし，Fig.5 に示したようにこれらを積分した出力係数は BEM とCFD で良く一致する結果となっている．
Fig.9 に𝜆 = 9.0における 3 次元 CFD での翼表面上の流線を示す．
翼端にて縦渦により腹面から背面へ流れが回り込む様子が見て取れる．また，腹面においてはく離が生じ，逆流した流れが遠心力により翼端へ押し流されている．Fig.10 は𝑟/𝑅 = 0.9における翼表面の流線と圧力分布を腹面側から表したものである．前縁付近のよどみ点から最大翼厚付近までの圧力降下と，後縁にかけての逆圧力勾配によるはく離域が確認できる．
このように，
3 次元 CFD
では BEM からは得られない翼周りの詳細な流れ場を確認することができる．

参考文献(1) Vermeer, N. J., “Performance measurements on a rotor model with
Mie-vanes in the Delft open jet tunnel”, Institute for Wind Energy,(1991).
(2) J. F. Manwell, J. G. MacGowan, and A. L. Rogers, “Wind energy
explained theory, design and application, 2nd edition,” John Wiley &Sons, Ltd., (2009).
(3) 鈴木，”数値計算による水平軸風車の性能予測，”日本機械学会論文集, B 編 79(808)，(2013)，pp. 2557-2561(4) Rhie, C. M. and Chow, W. L., “Numerical study of the turbulent
flow past an airfoil with trailing edge separation,” AIAA J. Vol.21
No.11,(1983), pp.1525-1532(5) Patankar, S. V., “Numerical heat transfer and fluid flow,”
McGraw-Hill, New York, (1980)(6) Leonard, B. P., “a stable and accurate convective modelling
procedure based on quadratic upstream interpolation,” Computer
Methods in Applied Mechanics and Engineering, (1979), pp.59-98.
(7) Launder, B. E. and Sharma, B. I., Letters in Heat Mass Transfer, 1,(1974), pp.131-138.
(8) Eriksson, L.E., “Generation of boundary conforming grids aroundwing-body configurations using transfinite interpolations,” AIAA J.
Vol.20 No.10, (1982), pp.1313-1320

Fig. 9 Limiting Stream lines on the blade surface.

Fig. 10 Pressure distribution of near the tip.

4Copyright © 2014 by JSFM