クロスフロー風車の出力特性と流れ場に関する数値計算 *
Numerical Simulation of the Flow Pattern and the Power Characteristic for a Cross-flow Wind Turbine
Yoichi A NDO **, Hideki A OKI *** and Sumio Y AMAGUCHI **
Key Words : Cross-flow, Wind Turbine, Wind power, Multiple Reference Frame
1. 緒言
現在，CO
2排出量の削減，化石燃料の枯渇や価格の 高騰などの面から無尽蔵でクリーンな自然エネルギーが 注目されている．自然エネルギーを利用した発電では火 力や原子力とは異なり，風力発電では風向と風速の依存 度が高く，大出力の風車の設置は風向が比較的安定し周 囲に障害物のない山間部や海岸などに限られ，そのほと んどがヨー機構を持たない大・中型の水平軸型風車であ る．これに対して，季節を通じて風向が変化し，1 m/s から 3 m/s の低風速である郊外や市街地のビルなどの 構造物隙間の設置に適した風車には，ヨー機構を持つ小 型の水平軸型風車やヨー機構が必要ない垂直軸型風車が ある．垂直軸型風車は水平軸型風車と比べると構造が簡 単で製作が容易であり，低速回転で高トルクを発生し起 動が容易である．さらに，回転速度が小さいので，水平 軸型風車と比べると静かである．
このような観点から，本報では図１に示すような垂直 軸型のクロスフロー風車について，設計の基礎データを 得るために必要な周囲の流れを数値計算により解析する 方法と結果を報告する．
2. 解析手法
クロスフロー風車は回転数が小さいので，風車周囲 の流れは非圧縮性として取り扱う．さらに，一様流中 に置いた風車の回転軸に垂直な断面の流れは，翼の固 定板周辺を除き二次元と考えてよい．また，風車周囲の 流れは非定常であるが，計算時間の短縮のために定常計 算を行い，時間平均流れ場を求めることにする．以上の 条件から，支配方程式には二次元非圧縮定常 Navier - Stokes 方程式 ( N.S. 方程式 )，乱流モデルには標準 k - e 低 Re を使用した．数値計算には汎用流体解析ソフト STAR – CD ( ver.4.06 ) を使用する．
* 平成 23 年 11 月 30 日受付 ** 機械工学科
*** ㈱アイペックス
安 東 洋 一 **
青 木 英 樹 ***
山 口 住 夫 **
図１ クロスフロー風車の形状
2 - 1. 支配方程式
風車はその構造から翼の迎え角が一回転する間に大き く変化する．また，風車全体の幾何形状は風車を構成す る翼枚数を n とすると，1/n 回転毎に元の状態と同じ になる．翼枚数が多い場合，この間の翼の迎え角の変化 は小さい．このような風車の回転運動を模擬した定常計 算を行うために，流れ場全体を 2 つの領域，つまり風 車回転軸を中心とする翼要素含む円形の解析領域 ( A ) とこれ以外の領域 ( B ) に分け，領域 A は回転座標系，
領域 B には静止座標系で表した支配方程式をそれぞれ 適応する．この方法では領域 A 内のメッシュを回転さ せる必要は無い．さらに，回転座標系の支配方程式 ( 連 続の式と N.S. 方程式 ) を式 (1),(2) に示すように変形す ると，角速度の値に応じて静止座標系または回転座標系 の支配方程式を表すことができる．このような解析領域 分割計算手法を Multiple Reference Frame ( MRF ) と いう．
2 – 2. ローター形状
本報のクロスフロー風車のローター形状と翼断面形 状を図２( a ) とその寸法を表１に示す．風車の翼断面 形状は翼弦長 20 mm，最大翼厚 4 mm，前後縁半径 1 mm であり，中心反り線が半径 20 mm の円弧翼であ る．翼はそり線の中点を翼の取付け位置として，直径 d= 100 mm の円上に 12 枚の翼を取り付ける．翼の取 付け点から引いた取付け円の接線と翼弦とのなす角，つ
まり取付け角は 45° である．この場合，ローター外径 d’
は 113.8 mm である．なお，風車の代表直径は翼の取 付け円直径とする．翼の各部の名称は図２( b ) に示す．
中心反り線を構成する円よりも外側の面を背面，内側の 面を腹面と呼び，翼縁がローター内側に位置する方を翼 内縁，ローター外側に位置する方を翼外縁と呼ぶ．
風車を図２( a ) 中に示した方向の一様流中に置くと，
この風車は時計方向に回転する．このとき，図に示す翼 の位置関係となるローターの回転角を基準状態とし，翼 の位置に A ～ L まで名称を付ける．角度θの定義は A 翼の取付け位置を０度として反時計方向を正とする．
2 – 3. 解析領域と境界条件
クロスフロー風車の全解析領域の概要を図３と寸法を 表２に示す．全解析領域の形状は W×L ( = 25d×32d ) の矩形とし，気流は左側面から右側面へ向かって流れる．
ローター回転中心軸から上流側へ向かって B ( = 10d ) の位置に流入境界面，下流側に 15d の位置に流出境界
図２ クロスフロー風車の概要と寸法
表１ クロスフロー風車のローター寸法および翼断面寸法 表２ 全解析領域の寸法 ( d = 100 mm )
面，上下 16d の位置にはスリップ境界面を配置した．
流入境界には U = 4 m/s の一様速度を与え，翼面は壁 面境界として翼に対する相対速度をゼロとする．この解 析領域に対する風車のブロッケージ効果を調べるため に，L = 220d の解析領域についても計算を行う．
全解析領域は，流れ場を MRF で解くためにローター 回転軸を中心とする直径 1.14d の円に囲まれた領域 A とこの領域を除いた領域 B に分ける．
2 – 4. 解析メッシュ
領域 A 内に配置する解析メッシュは翼周りに発生す る渦やはく離を捉えるため，メッシュを翼周りに集中 させる必要がある．したがって，翼弦長 a を代表長さ，
一様流速度 U を代表速度としたレイノルズ数を Re_
a としたとき，翼表面に形成される層流境界層の厚み は a/Re_a
0.5で見積もれる．境界層内には層状のメッシ ュ ( レイヤーメッシュ ) を配置し，乱流モデルに標準 k - e 低 Re を使用するので，壁面に接するメッシュの厚 み y+ が 1 程度になるようにする．メッシュの生成に は STAR - CD の機能の一つである Auto Mesh モジュ ールを使用する．この機能は，STAR - CD の解析に用 いられる様々な形状の解析メッシュを，メッシュ形状，
寸法および境界条件を与えることで自動的に作成する機 能である．図４には作成した領域 A 内の解析メッシュ ( AM0 ) の 1/4 領域を示す．図中の領域①では，メッ シュはレイヤーメッシュで作成しており，翼表面の法線 方向の 2 mm の厚みの中に 0.04 mm から 0.2 mm に 等比数列的 ( 比率 1.2 ) に厚みを増やして 20 層を作成 している．また領域②と③のメッシュについては，寸法 を 0.409 mm としている．図５に E 翼周囲のメッシュ を示す．
領域 B 内に配置するメッシュ ( BM ) はローターから
遠ざかるにつれて流れの変化は少なくなるので，領域 A との境界面のメッシュ寸法 0.409 mm から徐々に大 きくなるように配置する．全解析領域のメッシュ数はお よそ 20 万である．
本報に用いた風車の MRF 解析では，ローターが 1/12 回転する間に主流に対する翼迎え角の変化は 30° であり無視できない．そこで，上述の領域 A に配 置するメッシュを基本 ( AM0 ) として，これをロー ター回転軸の周りに 10° と 20° 回転させたメッシュ ( AM1, AM2 ) を用意する．つぎに，領域 B 内に配置 するメッシュ ( BM ) は領域 A と B の境界面において，
AM0, AM1, AM2 のいずれのメッシュと組み合わせた 場合もメッシュが連続に接続するようにメッシュ間隔を 調節する．流れ場の解析は３種類の解析メッシュの組み 合わせ （AM0 - BM, AM1 - BM, AM2 - BM ） について 行い，これらのメッシュを用いて求まる３つの定常流れ 場は，それぞれのローター回転角における位相平均流れ 場と考え，風車特性のトルク係数や出力係数などはこれ らの流れ場の平均値を用いる．
2 – 5. 緩和係数と収束判定
支配方程式および乱流エネルギー輸送方程式の解法に は Simple 法を使用する．Simple 法は速度と圧力およ
図４ 領域 A 内の解析メッシュ ( AM0 )
図５ E 翼周囲の解析メッシュ 図３ 全解析領域
び乱流エネルギーを反復計算によって求める．反復計算 におけるある時点での速度と圧力および乱流エネルギー をそれぞれ予測値とその修正値の和で表す．修正値は反 復の最終段階で予測値が真の解と一致した場合に 0 に なる．反復計算ごとに求まる修正値は過大に予測される ために発散しやすい．この問題を回避するために修正量 に緩和係数 ( < 1 ) を掛け，つぎの反復計算の速度と圧 力を求める．通常，速度の緩和係数は 0.7，圧力の緩和 係数には 0.3，乱流エネルギーの緩和係数は 0.7 が使用 される．
修正値を収束判定値に用いる場合はこの値を十分小さ くとる必要があり，解く流れ場によって適切な判定値は 異なる．STAR - CD ( Ver. 4.06 ) では，修正値に関す る情報はユーザー側から取得できないが，全メッシュに ついて支配方程式 ( 質量，運動量 ) の残差と乱流エネル ギー輸送方程式の残差の絶対値総和を取得できる
(2)．し たがって，これらを収束判定に使用する．また，渦の放 出などの非定常性の強い流れ場の定常計算では，流れ場 は反復計算の繰り返し回数に対し周期的に変化するの で，流れ場内の渦放出場所近傍の圧力や速度の変化を監 視して収束を判断する．表３に本報に用いた緩和係数を 示す．速度，圧力と乱流エネルギーの緩和係数は一般的 な値よりも低く設定しており，特に圧力の緩和係数は他 に比べて一桁下げている．これは残差が一般的な値を用 いると発散したためで，緩和係数を表２の値から，さら に下げても質量残差は 10
-4，運動量残差は 10
-5，乱流 エネルギー残差は 10
-4以下には下がらなかった．した がって，繰り返し回数に対し，B 翼の下流側 0.2d 付近 の圧力が一定または周期的に変化する場合を収束と判断 した．
3. 解析結果
解析を行うローターの周速比λ = d ω /(2U) を表３ に示す．周速比λに対する風車の出力係数 Cp = T ω / (1/2 ρ U
3A ) の変化を図６の青線で示す．ここで，ω はローターの角速度，U は流入境界の速度，T はロータ ーの軸トルク，ρは流体密度，A はローターの投影面積 である．出力係数 Cp は原点を通りλ = 0.26 で最大値 となる放物線状に変化しており，最大出力係数 Cpmax は 0.088，無拘束周速比は 0.65 である．図中の赤線は 谷野らが行った実験結果
(1)である．実験ではλ = 0.35 で Cpmax = 0.052 となり，無拘束周速比は 0.54 であ
った．
数値計算と実験を比較すると数値計算の Cpmax は 実験結果の 1.7 倍大きい．この原因は二つ考えられる．
一つ目は，ローターのブロッケージ効果である．解析 領域はローターとスリップ境界面の間隔が 16d の場 合，ブロッケージ係数βは 0.968 であり流路幅が不十 分であることが疑われる．そこでこの間隔を 110 d に 広げ，β = 0.995 として，λ = 0.26 の Cp を求めると 0.082 ( 実験結果の 1.6 倍 ) となる．これはβ = 0.968 の Cp と比較し 10 ポイントの低下であるので，数値計 算の Cpmax が実験値と比べて高い原因がローターのブ ロッケージ効果とは考えられない．二つ目は，収束条件 の残差が十分小さくなっていないことが考えられるが，
STAR – CD を使用する限り，これ以上残差を下げるこ とができないので，この確認はできない．しかしながら，
一般に残差が 10
-4程度まで低下すれば，速度場と圧力 場の定性的な評価は可能であるとされている．したがっ て，風車のトルク係数や出力係数を定量的に評価するた めには，翼に作用する力 ( 壁面のせん断応力と圧力 )，
つまり速度場と圧力場を正確に定量評価する必要があ り，残差を 10
-4からさらに小さくする必要があると考 えられる．
数値計算で最大出力係数を示した周速比λ = 0.26 に おけるローター内の各翼周りの詳細な絶対速度ベクトル 分布を図７に示す．図中のカラーコンターは風速の絶対 値を示す．図８には，図７で示した流れ場の圧力分布図 を示す．図７より，ローター付近の流れ場の特徴として は，流体はローターのθ = 120 ～ 240°( 図中の D ～ H 翼に相当 ) の間から流入し，一方θ =-30°( 330°) ～ 60°
( L ～ C 翼 ) の間から流出する．ローター内部には，時 計方向の渦が形成されている．この流れ場において，A
～ L 翼のそれぞれの位置によって翼面にあたる流れが 大きく異なり，翼に発生するトルクも大きく異なる．以 下に詳細を示す．
表３ 反復計算の繰返し回数に対する緩和係数の調整
図６ 数値計算と実験(1)の出力特性の比較
θ = 0°，30°，60° の翼 ( A，B，C 翼 ) では，翼内縁 を曲がる流速の速い流れが背面で剥離しており，このた め背面側では低い負圧となっている．さらに，翼の腹面 では流れが衝突して圧力が上がり，高い正圧となってい る．このように，この３枚の翼は衝動翼のように働くこ とによってかなり高い正のトルクを発生する．θ = 90°
の D 翼からθ = 240° の I 翼まではローターの回転に よる翼移動にともない翼内縁に対する流れの方向が逆と なっている．
また，θ = 90° の翼 ( D 翼 ) の翼外縁に流れが衝突し て圧力が上がり，翼がローターの回転する方向に抗力を 受けることによって正のトルクを発生する．
θ = 120° の翼 ( E 翼 ) では，D 翼周りの流れに加えて，
翼外縁を曲がる流速の速い流れが背面で剥離しており，
このため背面側では負圧となっている．このように翼が ローターの回転する方向により強い抗力と揚力を受ける ことによって正のトルクを発生する．
θ = 150°，180°，210°，240° の翼 （ F，G，H，I 翼）
では，A，B，C 翼と比較してローターの回転による翼 移動にともない翼内縁に対する流れの方向が逆となって いるため， A，B，C 翼の腹面で発生していた翼面に衝
突する流れが F，G，H，I 翼では背面に発生している．
さらに A，B，C 翼では翼内縁を曲がる流速の速い流れ が剥離し低い負圧となる面が背面にあったが，F，G，H，
I 翼では腹面に表れている．このように，この４枚の翼 は翼がローターの回転する方向とは逆に強い抗力を受け ることによってかなり高い負のトルクを発生する．
θ = 240° の位置で流れがローターから剥離しθ = 270°，300°，330° の翼 ( J, K, L 翼 ) 周囲では剥離域と なっている．
θ = 270° の翼 ( J 翼 ) 周囲ではほぼ無風状態となっ ており，トルクはほぼ零である．それに比べ，θ = 300°，
330° の翼 ( K，L 翼 ) では風車内に発生している渦の流 れが翼内縁から腹面付近に衝突して圧力が上がってい る．このようにこの 2 枚の翼は翼がローターの回転す る方向により強い抗力を受けることによって正のトルク を発生する．
ローターにおける風車の各翼に働くトルクを図９に 示す．翼枚数 12 枚の場合には，30° 毎に反時計方向の ( B ～ L ) 翼が対応する．縦軸は翼のトルクであり，正 がローターの回転方向と同じ向きのトルク，負がロータ ーの回転方向と反対向きのトルクである．図９よりロー ターは，θ =-90°( 270°) ～ 130°( 図中の J ～ E 翼に相 当 ) で正のトルクを発生しているが，θ =150° ～ 250°
( F ～ I 翼 ) ではトルクは負となっている．しかしなが ら，このローターは負のトルクよりも正のトルクのほう が高いため，ローター全体 ( θ = 0 ～ 350°) でのトルク が正の値となっている．
以 上 の こ と か ら， ロ ー タ ー に お い て は， 流 体 は 主 に θ =120 ～ 240°( E ～ I 翼 ) の 間 で 流 入 し，
θ =-30°( 330°) ～ 60°( L ～ C 翼 ) の間に流出する．ま た，ローター内部では時計回転の渦が形成されている．
図７ ローター周囲の絶対速度ベクトル
図８ ローター周囲の圧力分布
図９ ローターの１枚の翼が回転中に発生するトルク
図５の圧力分布を見ると流入域のうちθ = 120°( E 翼 ) では翼の揚力によって正トルクが，θ = 150 ～ 240°( F
～ I 翼 ) では抗力によって負トルクが生じる．これに対 し，流出領域であるθ =-30°( 330°) ～ 60°( L ～ C 翼 ) は抗力によって正のトルクが生じていることが分かる．
したがって，翼に働く力はθ = 80° 付近で揚力から抗力 に切り替わる。またこの付近は上流側の翼の後流であり、
翼トルクが小さくなるのでトルクの変化はこの位置で不 連続となる．さらにこのローターは流入側の揚力よりも 流出側の抗力によって生じる正トルクのほうが高い抗力 型の風車であることがわかる．
4. 結言
クロスフロー風車の出力特性と流れ場について数値計 算を行い以下の結論を得た．
(1)STAR-CD ( Ver.4.06 ) を使用した場合，十分な残差 の低下を得られず，トルクと出力係数の正確な予測がで
きなかった．
(2) ローターにおいては，流体は主にθ =120 ～ 240°
の間で流入し，θ =-30°( 330°) ～ 60° の間に流出する．
また，ローター内部では時計回転の渦が形成されている．
(3) 流入域のうちθ = 120° では翼の揚力によって正ト ルクが，θ = 150 ～ 240° では抗力によって負トルク が生じる．これに対し，流出領域であるθ =-30°( 330°)
～ 60° は抗力によって正のトルクが生じる．
(6) このローターは流入側の揚力よりも流出側の抗力に よって生じる正トルクのほうが高い抗力型の風車であ る．
参 考 文 献
