電力中央研究所報告

電 力 輸 送
電力中央研究所報告
発電用風車の動的応答解析コードの開発
研究報告 : Ｎ08015
平成21 年 4 月
財団法人
電力中央研究所
地球工学研究所
発電用風車の動的応答解析コードの開発
幹夫*
清水
Key Words：Wind power generation
Wind turbine
Dynamic response analysis
Finite element method
Finite rotation
キーワード：風力発電
風車
動的応答解析
有限要素法
有限回転
Development of a dynamic response analysis code for wind turbines
by Mikio Shimizu
Abstract
Introduction of wind power generation is growing in recent years. In addition, there are various wind turbine design
codes that can be used to estimate wind turbine behavior. However, the number of degrees of freedom that can be modeled
by many of the design codes is limited, because the codes adopt a conventional modal approach as a dynamic analysis
method in order to execute efficient calculation. In this study, in order to eliminate restriction of the number of degrees of
freedom in modeling, we developed an original code for wind turbine dynamic response simulations. The code is based on
the finite element method which adopts a direct integration method and it can take geometrical nonlinearity including finite
rotation into account for modeling turbine rotation. In this paper, the formulations of the developed code are shown and a few
examples of large deformation problems and wind turbine behavior are analyzed to verify the developed code. The results
calculated by the developed code are in fairly good agreements with the results obtained in terms of a literature or an existing
design code “FAST”, and the usefulness of the developed code is verified.
（Civil Engineering Research Laboratory Rep.No.N08015）
（平成 21 年 2 月 19 日 承認）
*
地球工学研究所 構造工学領域
上席研究員
i
©CRIEPI
背
景
発電用風車については，近年，日本の自然環境を考慮した，設計合理化への国内
独自の取り組みが進み，「日本型風力発電ガイドライン」や「風力発電設備支持物構
造設計指針・同解説」などが発行されている．しかしながら，強風時におけるタワー
の倒壊，ブレードの損傷などの事故に際し，その荷重あるいは挙動評価には，依
然として海外で開発された設計プログラムが用いられる傾向が強い．特に，国際規
格の認証を受けた汎用プログラムについては，設計計算の効率化が図られ，ねじり
モードの無視および主要な振動モードのみの選択により解析自由度が制限されてお
り，事故原因究明のための詳細な損傷評価や高次モードを含む振動現象のシミュレー
ションに適しているとは言い難い．このため，解析自由度を制限することなく発電
用風車の風応答を解析し得るプログラムの開発が必要と考えられる．
目
的
発電用風車を対象として，ロータの回転運動を含む風応答シミュレーションに適
用可能な，動的応答解析プログラムを開発する．
主な成果
（1）コードの開発
既往の解析プログラムの多くが，計算の効率化のためモード重畳法を採用し，自
由度数を制限している現状を鑑み，直接積分法に基づく有限要素解析コードを開発
した．これにより，自由度数を制限することなく，タワー，ブレード各部の応答変
位および断面力の算出を可能とした．また，この開発コードは，ロータの回転を含
む風車の挙動を精度良く模擬するため，時々刻々変化しながら作用する風荷重を準
定常的に近似する機能を有し，はり要素の有限回転を考慮し得るものとした．
（2）有限回転問題の解析機能の検証
開発コードの有限回転に関する解析機能を用い，固定放物線アーチの面外座屈問
題を解析し，既往の文献にみられる解との比較を行った．解析結果は，文献にみら
れる解と概ね一致し，面外座屈を表現した．これにより，ロータの回転および風車
各部の大変形挙動の模擬に対する，開発コードの適用性を確認した．
ii
©CRIEPI
（3）風車の風応答解析機能の検証
開発コードの風応答解析機能を検証するため，一様流が作用する 2 枚および 3 枚
ブレード風車の挙動を解析し，既往の解析プログラム FAST との間で解を比較した．
その結果，開発コードによって算出されたブレード先端の挙動には，ロータの回転
を示す定常的な変位がみられ，定常状態におけるロータの回転数は，FAST による解
と一致した．さらに，両コードの離散化手法などの相違により過渡応答波形が異な
るものの，タワー基礎の風向方向曲げモーメントおよびロータの軸力は開発コード
と FAST との間で概ね一致する結果となった．以上により，開発コードにより，ロー
タの回転を含む風車の風応答挙動の解析が可能であることが示された．
今後の展開
開発コードを，設計時の詳細な動的照査あるいは事故時の挙動分析などに適用し，
その信頼性ならびに解析精度の向上を図る．
iii
©CRIEPI
目
次
１．はじめに ···················································································································· 1
２．プログラムの解析対象と開発方針 ··················································································· 2
2.1 解析対象················································································································ 2
2.2
開発方針と機能概要································································································· 2
３．風車の風応答の定式化 ·································································································· 3
3.1 計算フロー············································································································· 3
3.2
構造系および空気力のモデル化 ·················································································· 4
3.3
大変形問題の定式化································································································· 7
４．開発コードの機能検証 -1- 大変形問題の解析機能の検証 ···················································10
4.1
片持ち梁のエラスティカ··························································································10
4.2
放物線固定アーチの面外座屈 ····················································································10
５．開発コードの機能検証 -2- 風車の風応答解析機能の検証 ···················································11
5.1 FAST の概要··········································································································12
5.2 解析対象···············································································································12
5.3 計算条件···············································································································12
5.4 解析結果···············································································································16
６．試計算による風車の挙動評価例 ·····················································································18
6.1
タワーの抗力係数の影響··························································································18
6.2
風向とロータ回転数との関係 ····················································································18
７．まとめと今後の展開 ····································································································19
謝辞 ·······························································································································20
参考文献 ·························································································································20
iv
©CRIEPI
証を受けており，発電量予測，断面力計算を含め，
1.　はじめに
発電用風車の統合的な性能評価に国内外で広く
近年，地球温暖化が問題視される中，環境負
用いられている 5）．ここに，GL Wind Guideline
担が少ない新エネルギーの一つとして，風力発
は，ドイツ船級協会（Germanischer Lloyd）が
電への期待，関心が高まっている．政府が 2010
制定した発電用風車の規格である．
年の風力発電導入目標を 300 万 k W とする一方，
た だ し， 上 述 の B L A D E D を は じ め， 既 往 の
独立行政法人新エネルギー・産業技術総合開発
design code の多くについては，設計計算の効
機構（以下，N E D O と称する）は，2005 ～ 2007
率化を図って主要な振動モードのみが選択され
年度の 3 カ年にわたる事業により，風力発電の
るなど，風応答計算における自由度が制限され
利用率向上および導入促進を図り，「日本型風
ている 6）．このため，ブレード，タワーなど構
力発電ガイドライン」1）を策定している．また，
造系各部についての詳細な損傷評価，および高
2004 年には，土木学会構造工学委員会に「風力
次モードが含まれる振動現象のシミュレーショ
発電設備耐風設計小委員会」が設置され，3 年
ン，などへの適用を想定した場合，自由度を制
間の活動により「風力発電設備支持物構造設計
限することなく風応答を解析し得るプログラム
指針・同解説」2）が上梓された．
の整備が必要と考えられる．また，運用上の利
こうした事業，活動の背景として，日本特有
便性を勘案すれば，電気事業として国産の解析
の自然環境条件，特に台風により，発電用風車
プログラムを所有する意義は大きく，また，そ
に倒壊，ブレードの損傷などの事故が発生して
の要望も高まりつつある．
3）
こと，および日本の実情を勘案した風力
以上のような現状を鑑み，本研究では，ガイ
発電設備の導入法，設計法が必要とされている
ドラインや設計指針類の策定と同様，構造解析
こと，などが挙げられる．従来，国内の風力発
についても国内独自の取り組みが確立されるこ
電設備の設計においては，先行的に導入が進ん
とを目標として，発電用風車の風応答解析プロ
でいる海外の規格や基準類が適用される傾向が
グラムの開発を試みた．本研究では特に，モノ
いる
みられたが，上記 NEDO のガイドライン
1）
ポール支持式のプロペラ式風車を対象として，
および
土木学会の設計指針 2） は，こうした傾向を改善
ロータの回転運動を大変形問題として模擬し得
する国内独自の取り組みということができる．
る，直接積分法に則った 3 次元有限要素解析プ
しかしながら，上記事故や設計に関連する発
ログラムの開発を目的とした．
電用風車の強度，断面力あるいは挙動評価に必
本報告には，プログラムの解析対象となる風
要不可欠な構造解析手法，プログラムについて
車の概要，プログラムの開発方針と機能概要，
は，国内における開発，整備は積極的にはなさ
風車挙動を大変形問題として扱う場合の定式化，
れていないように思われる．この原因として，
機能の検証として既往の文献，あるいは代表的
海外においては，既に，“d e s i g n c o d e”と称
な風車の design code である「FAST」7）との解
される発電用風車の設計に特化した計算プログ
の比較結果，および試計算による風車挙動の評
ラムが多数存在し，国内でも，これらの中の幾
価例を示す．
つかを使用できることが挙げられる．例えば，
Garrad Hassan 社が開発した「BLADED」4）は，
設計と検定を目的とした発電用風車の断面力計
算プログラムとして，GL Wind Guideline の認
－1－
©CRIEPI
2.　プログラムの解析対象と開発方針
࣮ࣟࢱ⣔
2.1　解析対象
ᨭᣢ࣭ᵓ㐀⣔
ࣈ࣮ࣞࢻЍ
本研究で開発したプログラム（以下，開発コー
ࣁࣈ㸪࣮ࣟࢱ㍈Ў
ドと称する）の解析対象は，発電用風車のうち，
Ўࢼࢭࣝ
これまでの事故件数が多く，今後の導入の可能
性が高い，鋼製モノポール支持式のプロペラ式
風車（以下，風車と称する）であり，そのロー
タ系および支持・構造系モデルとなる．図 1 に
Ћࢱ࣮࣡
一般的な風車の外観を示す．図のとおり，ロー
タ系は回転翼であるブレード，その回転軸とな
るロータ軸，およびブレードの付け根をロータ
軸に連結するハブで構成される．また，支持・
Ћᇶ♏
構造系は，発電機の収納部であるナセル，ロー
タ系およびナセルを指示するタワー，およびタ
図 1　風車の外観
ワーを支える基礎で構成される．
現在，風車のブレードは，ガラス繊維強化プ
ラスチック製で 3 枚とされることが多く，一般
停止あるいは遊転状態にあるカットアウト時を
的に定格出力が 600k W の場合，ロータ直径は 45
基本として，ブレード，タワーなどに高次モー
～ 50m，タワー高さは 40 ～ 50m，定格出力が 1,000
ドの部材振動が発生した場合も想定した，風車
～ 2000k W の場合，ロータ直径，タワー高さとも
の挙動評価，損傷評価への適用性を重視した．
2）
60 ～ 90m となる ．
特に，風車設計への適用を目的とした既往の
なお，プロペラ式以外の水平軸風車にはオラ
プログラムの多くは，表 1 に示すとおり，構造
ンダ式，セルウイング式，多翼式などがあり，
系のモデル化にモード重畳法を採用し，振動モー
8）
鉛直軸風車も存在する ．また，タワーには鋼
ドを取捨選択している 6）．このため，開発コー
製モノポールの他，ラチスタワーやコンクリー
ドは，自由度数と等しい振動モードを考慮し得
トタワーが存在する．
る，直接積分法を採用した非線形有限要素解析
コードとした．開発コードの機能概要を表 2 に
示す．
2.2　開発方針と機能概要
モード重畳法には，主要な振動モードのみを
プログラムの開発については，風の作用を受
考慮し，自由度を減らして効率的な応答解析を
ける風車の応答を時間領域で確定論的にシミュ
実施し得る利点がある．また，ロータの挙動を
レーションすることに主眼を置き，風速の入力
回転モードとして扱うことにより，通常運転時
に対して，風車各部の変位および断面力を任意
の解析の安定性が向上すると考えられる．表 1
に出力し得ることを条件とした．また，ロータ
より，既往のプログラムは風車に特化した設計
が回転状態にあるカットインからカットアウト
への適用を目的として，効率的で安定した計算，
まで（以下，通常運転時と称する），および回転
発電設備としての性能評価を重視していること
－2－
©CRIEPI
表 1　既往の解析プログラムの自由度 5)，6）
項目
開発国
開発元
Garrad Hassan
Oregon St.
Univ.
National
Renewable
Energy Lab.
モードの取捨選択
Edge 方向 *
ブレード 考慮する
Flap 方向 *
モード数 ねじり方向
モードの取捨選択
前後方向 *
考慮する
タワー
左右方向 *
モード数 ねじり方向
有
6
6
0（無視）
有
3
3
0（無視）
モード重畳
法採用
有
1
2
0（無視）
有
2
2
0（無視）
モード重畳
法採用
無
有
有
無
無
無
ユーザ設定
1
2
ユーザ設定 ユーザ設定 ユーザ設定
ユーザ設定
1
2
ユーザ設定 ユーザ設定 ユーザ設定
ユーザ設定 0（無視） 0（無視） ユーザ設定 ユーザ設定 ユーザ設定
無
無
有
無
無
有
ユーザ設定 ユーザ設定
2
ユーザ設定 ユーザ設定
8
ユーザ設定 ユーザ設定
2
ユーザ設定 ユーザ設定
8
ユーザ設定 ユーザ設定
1
ユーザ設定 ユーザ設定
1
モード重畳 モード重畳
モード重畳
法採用
法採用
法採用
備　考
FAST
アメリカ
ADAMS/WT
アメリカ
プログラム名
DUWECS
FLEX5
オランダ デンマーク
BLADED
イギリス
Delft Univ.
GAST
ギリシャ
HAWC
PHATAS-IV
デンマーク オランダ
Dutch Energy
National Tech. Wind Energy
Tech. Univ. of
Dept. of Risø Research
Univ. of
Denmark
National Lab. Found. (ECN)
Athens
注 *：Edge 方向はブレード断面の翼型の弦方向，これと直交する方向が Flap 方向．前後，左右方向は，それぞれ通
常運転時におけるロータ軸方向，これと直交する方向．
表 2　開発コードの機能概要
がわかる．
積
分
法
離 散 化 法
使 用 要 素
未
知
量
参照座標系
大
変
形
収束計算法
空
気
力
構造系の減衰
これに対し，開発コードについては表 2 のと
おり，詳細な挙動評価を目的として，一般的な
構造解析手法を風車の風応答シミュレーション
に適用する方針を採った．すなわち，構造系の
モデル化に際し，自由度および出力点を設定す
る上での柔軟性，拡張性を勘案して，離散化に
直接積分（Newmark - β）法
有限要素法
線形はり要素，ばね要素
変位
Updated - Lagrange 系
幾何学的非線形性，有限回転
Newton - Raphson 法
風速の時刻歴から準定常的に近似
Rayleigh 減衰
は有限要素法を採用した．これにより，構造系
各部の断面力を詳細にみることが可能となる．
設計の目的により異なる場合がある 6）．したがっ
ブレード，ナセルおよびタワーは，はり要素で
て，開発コードについては，現状では，入力さ
モデル化し，ロータの回転を模擬するため，有
れる風速を準定常的に空気力に近似するのみに
限回転および幾何学的非線形性を考慮するとと
止め，ロータの回転により流れ場に生ずる誘導
もに，ロータ軸とタワー頂部とが接合される節
速度は，今後の適用状況に応じて，上記の補正
点には，回転方向に機能するばね要素の適用を
定数の導入によって考慮する．
可能とした．
また，直接積分法は，高次モードの挙動を含
3.　風車の風応答の定式化
めた，非線形性のある振動解析に適していると
いえる．
3.1　計算フロー
なお，空気力については，表 1 に示した既往
のプログラムのすべてが，翼素運動量理論 8） に
開発コードは有限要素法に基づき，離散化さ
基づき，ロータ回転面内の流れ場における誘導
れた各要素に対して局所座標系で成立する運動
速度の影響を含めたモデル化をしている．誘導
方程式を，全体座標系に変換した上で解くこと
速度は，経験的に得られた補正定数を用いるこ
により，未知量すなわち変位を算出する．この
とにより，準定常的な空気力の近似において考
過程を図 2 の全体フローに示す．
慮可能と考えられる．ただし，補正の方法は，
また，開発コードでは，ロータの回転を大変
－3－
©CRIEPI
❶変形前の要素剛性行列読込・
座標変換
①入力データ読込
③要素質量行列
作成
⑦全体節点荷重ベ
クトル作成
④要素線形・幾何
剛性行列作成
⑧ 積 分 (Newmark- β
法適用)
⑤要素減衰行列
作成
⑨変位増分算出
⑥各要 素行列の
座標変換・総和
⑩残差力の計算・変
位補正・収束判定
変位補正ステップ分（残差が収束するまで）繰返し
時間ステ ップ分 繰 返 し
②全体質量・減衰・
剛性行列作成
←図 3 の❶～❽が
含まれる．
❷変形後の要素の反力ベクトル
を全体座標系で作成
❸要素の残差力ベクトルを全体
座標系で作成
❹全要素の反力ベクトル・残差
力ベクトルの総和
❺変形後の要素剛性行列作成・
座標変換・総和
❻残差力ベクトルの二乗和算
出・収束判定
❼積分(Newmark-β法適用)
⑫出力
❽変位補正量算出・変位補正
図 2　開発コードの全体フロー
収束
⑪座標値および構
造データ書換
図 3　残差力の計算・変位補正・収束判定過程
形問題として扱い，有限回転，幾何学的非線形
荷重
性を考慮する．こうした幾何学的非線形解析で
は，ある計算ステップ（例えば，t ～ t+Dt 間，
残差力
仮定値→
D は増分を表す）において，実際には変化する
反力| t+Dt
真の剛性を接線剛性で仮定することによる不釣
←真値（未知）
合力，すなわち残差力が生じる．この残差力を
接線剛性
0 に収束させるため，全体フローの⑩の段階に
反力| t
変位増分
含まれる，図 3 の変位補正過程が必要不可欠と
なる．図 4 には，残差力発生の概念図を示す．
t
次節以下には，フローの中で参照される定式
t+Dt
変位
図 4　残差力発生の概念図
化について示す．
3.2　構造系および空気力のモデル化
風車各部は，図 5 に示す，はり要素によって
離散化される．ある 1 要素について，両端にあ
3.2.1　構造系のモデル化
る節点 i，j の各断面の重心点を a，b，これらを
構造系のモデル化すなわち運動方程式の記述
結ぶ直線すなわち要素軸を変形前後の x 軸とし，
は，図 2 のフローの②～⑥の段階でなされる．
これと直交右手系をなす y，z 軸を定める．また，
－4－
©CRIEPI
vb
ここに，{ } T はベクトルの転置を表す．要素の
qb
va
y
i,a
[ m ]D {u} + [ c ]D {u }
+ ([ k L ] + [ kG ]) D {u} = D { p} (2)
ここに，[m]：要素質量行列，[c]：要素減衰行列，
ua
ya
x
示できる．
wb
j,b zb
qa
運動方程式は，式 (2) のとおり増分形で行列表
ub
yb
wa
[kL]：要素線形剛性行列，[kG]：要素幾何剛性行列，
za
{p}：節点力ベクトルであり，(・) は時間微分を
z
表す．[m]，[k L] および [k G] は式 (3) ～ (5) で表
図 5　はり要素と局所座標系・節点の自由度
される 9）．
式 (3) ～ (5) において，r：要素密度，A：要
節点変位成分を図 5 中のとおり定めれば，節点
素断面積，l：要素長，I y，I z：それぞれ y，z 軸
変位ベクトル {u} は次式で表される．
回りの断面二次モーメント，J x：x 軸回りの極慣
性モーメントである．また，E：ヤング率，G：
{u} = {ua , va , wa , qa , ya , z a , ub , vb , wb , qb , yb , z b }
T
æ1
çç
ççç 3
çç
6I z
13
+
çç 0
2
35 5 Al
çç
çç
çç 0
0
çç
çç
çç 0
0
çç
çç
çç
0
ççç 0
çç
çç
11l
I
+ z
çç 0
210 10 Al
ç
[ m ] = r Al ççç 1
çç
0
çç 6
çç
6I z
9
çç 0
2
çç
70 5 Al
çç
çç
0
çç 0
çç
çç
çç 0
0
çç
çç
çç 0
0
ççç
çç
I
13l
+ z
çççè 0 420 10 Al
13
35
+
6I y
5 Al
11l
210
Jx
3A
-
Iy
10 Al
9
+
0
0
0
0
6I y
5 Al
2
0
0
-
13l
420
0
0
l
+
2
105
+
-
l
2
140
0
2I z
15 A
1
0
13l
420
Iy
10 Al
Iy
30 A
-
l
3
Iz
-
0
6A
Iy
15 A
0
Jx
10 Al
2I y
2
0
0
13l
l
105
0
-
420
0
0
0
70
sym.
2
0
-
せん断弾性係数，F：軸力（張力），f y，f z：せ
………(1)
10 Al
13
35
+
6I z
5 Al
0
0
0
0
0
0
0
0
0
2
140
0
-
Iz
30 A
0
-
11l
210
-
2
13
35
+
6I y
5 Al
Jx
0
11l
210
Iz
10 Al
+
3A
Iy
10 Al
0
－5－
2
0
0
l
2
105
+
0
2I y
15 A
ö÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
2
2 I z ÷÷
l
÷÷
+
105 15 A ø
…………………(3)
©CRIEPI
æ EA
çç
çç l
çç
çç 0
çç
çç
çç
çç 0
çç
çç
çç
çç 0
çç
çç
çç 0
çç
çç
çç
çç 0
ç
[ k L ] = ççç EA
çççç l
çç
ççç 0
çç
ççç
çç 0
ççç
çç
ççç 0
çç
ççç
çç 0
ççç
çç
ççç 0
è
ö÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷÷
÷÷
÷÷
÷÷÷
÷÷
÷÷
÷÷÷
÷÷
÷÷
÷÷÷
÷÷
÷÷
÷÷
÷÷
÷÷÷
÷÷
(4 + fy ) EI z ÷÷
÷÷
l (1 + fy ) ÷ø
12 EI z
3
l (1 + fy )
12 EI y
0
sym.
3
l (1 + fz )
0
GJ x
0
l
-6 EI y
0
6 EI z
2
l (1 + fy )
0
-12 EI z
0
0
0
0
0
0
-12 EI y
0
0
6 EI z
0
2
l (1 + fy )
0
2
l (1 + f y )
0
(2 - fz ) EI y
0
0
0
0
0
0
(2 - fy ) EI z
0
12 EI z
3
l (1 + f y )
0
0
l (1 + fz )
0
l
-6 EI z
6 EI y
0
2
l (1 + fz )
EA
0
2
l
-6 EI y
l (1 + fy )
l (1 + fz )
-GJ x
0
(4 + fy ) EI z
0
0
3
l (1 + fz )
0
l (1 + fz )
0
0
3
l (1 + fy )
(4 + fz ) EI y
0
2
l (1 + fz )
0
l (1 + fy )
0
-6 EI z
2
l (1 + fy )
12 EI y
3
l (1 + fz )
GJ x
0
l
6 EI y
2
l (1 + fz )
0
(4 + fz ) EI y
0
l (1 + fz )
0
0
æ0
çç
çç
6
çç0
5
çç
çç
6
çç0 0
çç
5
çç
0
çç0 0
çç
çç0 0 - l
çç
10
çç
çç
l
0
ç0
Fç
10
[ kG ] = l ççç0 0
0
çç
çç
çç0 - 6
0
çç
5
çç
6
çç
çç0 0
5
çç
çç0 0
0
çç
çç
l
çç0 0 10
çç
çç
ç0 l
0
èçç
10
sym.
0
0
2l
15
0
0
0
0
0
0
0
0
2l
0
-
10
0
-
l
2
15
l
0
0
2
0
l
10
6
0
5
0
0
0
0
0
0
6
5
0
2
0
30
0
-
l
0
2
30
0
l
0
-
10
l
10
0
0
0
0
2l
2
15
0
ö÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷ ÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
÷÷
2 ÷
2l ÷÷
÷÷
15 ø
…………………(4)
…………………………………………… (5)
ん断変形のパラメータであり，それぞれ y，z 方
モード減衰定数をそれぞれ w n，h n としたとき，
向のせん断に対する有効せん断面積 Ay，A z を用
次式より求められる．
い，式 (6) のとおり表される．
fy =
12 EI z
GAy l
2
,
fz =
12 EI y
GAz l
2
(6)
hn =
a
2w n
+
bw n
2
………(8)
要素の減衰は，Rayleigh 減衰として，次式で
3.2.2　空気力のモデル化
表される．
風車各部には準定常的に近似した空気力とし
[ c ] = a [ m ] + b ([ k L ] + [ kG ]) (7)
ここに，a，b は，第 n 次モードの固有円振動数，
て，タワーおよびナセルには風向方向に抗力の
みが，ブレードには抗力，揚力および空力モー
－6－
©CRIEPI
メントが作用すると考える．以下には，図 2 の
モーメント係数である．また，図中の s は空気
フローの⑦の段階に相当する，ブレードに作用
力の作用点である．
する空気力の定式化を示す．
タワーおよびナセルに作用する空気力として
図 6 に示すように，ブレード断面に対して，
は，抗力のみを式 (12) に基づいて考慮する．
はり要素と同様，軸方向（紙面直交方向）に x 軸，
以上の単位長さ当りの空気力は，受風面積を
これと直交右手系をなす y，z 軸を定める．主風速，
考慮し，節点力の増分に換算される．
これと直交する風速成分をそれぞれ W，V とし，
ブレード速度の主風向，これと直交する方向の
3.3　大変形問題の定式化
成分をそれぞれ w ， v とすれば，相対風速 Vr お
よび相対風速と主風向とのなす角 j は，それぞ
図 2 のフローの②～⑦，⑩，および図 3 のフ
れ次式 (9)，(10) で表される．
ローの❶～❺の段階における，各要素行列の作
Vr =
2
2
(V - v ) + (W - w ) æW - w ö÷
çè V - v ÷÷ø 成，座標変換あるいは要素の反力ベクトルの作
(9)
成などにおいては，時間および変位補正ステッ
プごとに変化する局所座標系を決定し，座標変
-1
j = tan çç
………(10)
換行列を構成する必要がある．その際，開発コー
また，翼弦の傾き角 q および迎え角 x を図中の
ドでは，風車の挙動，特にロータの回転を大変
とおり定義すれば，次式の関係となる．
形問題として扱い，前田・林 10）の定式化に倣っ
x = j -q て有限回転を考慮する．以下にその概要を示す．
…………………(11)
単位長さ当りのブレードに作用する空気力の三
成分すなわち抗力 F D，揚力 F L および空力モー
3.3.1　回転行列
メント F M は，上記 Vr と x を用いて次式で表す
図 7 のとおり，ある節点 o が空間ベクトル f
ことができる．
を中心に有限な角 a だけ回転した場合を考える．
FD =
FL =
FM =
1
2
1
2
1
2
r a dC D (x )Vr ………(12)
r a dC L (x )Vr ………(13)
2
2
ベクトル f の全体座標系に対する方向余弦を (fX,
f Y , f Z ) とし，ベクトル f の方向に単位ベクトル
e 3 を，これと直交右手系をなすように 2 つの単
f
r a d CM (x )Vr ………(14)
2
2
a
こ こ に，r a： 空 気 密 度，d： 断 面 の 代 表 長 さ，
CD(x)：抗力係数，CL(x)：揚力係数，CM(x)：空力
a
z
FL
W- w
q
x V
r
j
r
FD
y
s
r*
e3
FM
e1
V- v
図 6　ブレード断面に作用する空気力
o
e2
図 7　ベクトルの回転
－7－
©CRIEPI
位ベクトル e1，e2 を定める．
点は，それぞれ変形前，後とも称する．
このとき，有限な回転に対する回転行列 Ro は，
いま，図 8 のとおり要素の両端 a，b に固定さ
次式で与えられる．
れた 2 組の直交座標系を考える．これらは変形
Ro = E cos a + (1 - cos a) e3 Ä e3 + e3 ´ E sin a
………(15)
前，(xa,t*, ya,t*, za,t*)，(xb,t*, yb,t*, zb,t*) にあるものとし，
t から t+Dt の間の節点 i，j の回転をそれぞれ回
こ こ に，E： 単 位 行 列， Ä ： ダ イ ア ド 積， ×：
転行列 Ri，Rj で表せば，変形後は式 (16) により，
ベクトル積を表す．R o により，節点に固定され
x a ,t +Dt = Ri x a ,t , ya ,t +Dt = Ri ya ,t , z a , t +Dt = Ri z a , t
*
た任意のベクトル r は，次式により r に変換さ
2
R
o
*
*
*
*
*
*
*
*
*
*
となる．
*
r = Ro r …………………(16)
また，Ro の各要素は次式で表される．
X
*
xb ,t +Dt = R j xb ,t , yb ,t +Dt = R j yb ,t , zb ,t +Dt = R j zb ,t
………(18)
れる．
é f
ê
= êf f
ê
êf f
ë
*
a 端において，変形後の x 軸上に単位ベクト
(1 - cos a ) + cos a
ル x a,t+Dt を考え，x a,t+Dt * となす角を a a とする．
Y
X
(1 - cos a ) + fZ sin a
次に，x a,t+Dt* と x a,t+Dt とに直交する軸を回転軸と
Z
X
(1 - cos a ) - fY sin a
して，直交系 (x a,t+Dt*, y a,t+Dt*, z a,t+Dt*) を角 a a だけ
fX fY (1 - cos a ) - fZ sin a
ú
ú
ú
û
fY (1 - cos a ) + cos a
fY fZ (1 - cos a ) - fX sin a ú
fZ fY (1 - cos a ) + fX sin a
fZ (1 - cos a ) + cos a
2
回転させ，x a,t+Dt* を x a,t+Dt に一致させる．このと
fX fZ (1 - cos a ) + fY sin a ù
2
………(17)
き，ya,t+Dt*，za,t+Dt* はそれぞれ ya,t+Dt，za,t+Dt になっ
たものとする．こうした関係は次式で表される．
3.3.2　局所座標系の更新
*
x a ,t +Dt = Ra x a ,t +Dt ,
ya , t +Dt = Ra ya , t +Dt ,
*
*
図 2，図 3 のフローには，それぞれ時刻歴解
z a , t +Dt = Ra z a ,t +Dt ………(19)
上式中の Ra は，式 (15) により次式で表される．
析の時間ステップ，残差を収束させるための変
Ra = (x a ,t +Dt × x a ,t +Dt )E
*
位補正ステップのループが示されている．局所
座標系は，これらステップごとに更新される．
+
(x
*
a , t +Dt
´ x a , t +Dt )(x a , t +Dt ´ x a , t +Dt )
*
1 + (x a , t +Dt × x a , t +Dt )
*
れぞれ t，t+Dt と表し，この間に図 5 に示した
+ (x a ,t +Dt ´ x a ,t +Dt )´ E (20)
ここに，・：スカラー積を表し，上式では，回転
要素が図 8 のように変形する場合を想定して，
軸方向の単位ベクトルを fa とした場合，
ステップ終了時の局所座標系の yt+Dt, zt+Dt 軸の導
出を示す．x t+Dt 軸については，前述のとおり，
x a ,t +Dt ´ x a ,t +Dt = fa sin aa の関係があることを用いている．
両端 a，b を結ぶ直線すなわち要素軸が変形前後
b 端においても同様にして，直交座標系 (xb,t+Dt,
の x 軸となる．なお，ステップの開始，終了時
yb,t+Dt, zb,t+Dt) を得ることができる．
ここでは，1 ステップの開始，終了時点をそ
*
*
(21)
a，b 端において，ベクトル y a,t+Dt と z a,t+Dt およ
び y b,t+Dt と z b,t+Dt とが作る平面は互いに平行であ
xb*
i
a
ya*
*
aa xa
xa
za*
b
yb*
j
ab
るため，y a,t+Dt と y b,t+Dt または z a,t+Dt と z b,t+Dt のな
xb
す角 f t+Dt が，要素の相対ねじり角となる．f t+Dt
は次式の関係によって得られる．
zb*
y
×y
= z a ,t +Dt × zb ,t +Dt = cos ft +Dt (22)
a ,t +Dt b ,t +Dt
- ya ,t +Dt × zb ,t +Dt = z a ,t +Dt × yb ,t +Dt = sin ft +Dt (23)
図 8　はり要素の節点に固定された座標系
－8－
©CRIEPI
a，b 両端のねじり角 qa,t+Dt*，qb,t+Dt* を，絶対値
3.3.4　要素の反力ベクトルの定式化
が等しく，符号が逆になるように，
図 4 に示すとおり，残差力は反力に基づいて
*
*
-qa ,t +Dt = qb ,t +Dt =
ft +Dt
2
得られる．開発コードでは，図 3 に示したフロー
= qt +Dt (24)
の❷の段階において，以下の定式化により要素
と 選 び，y a,t+Dt と y b,t+Dt と に 直 交 す る 軸 回 り に
の反力ベクトルを求めている．
q t+Dt だけ回転させる回転行列 R q を用いれば，変
まず，変形後の局所座標系における節点変位
形後の yt+Dt, zt+Dt 軸は次式で定められる．
ベクトル {u} を式 (1)，
{u} = {ua , va , wa , qa , ya , z a , ub , vb , wb , qb , yb , z b }
T
yt +Dt = Rq ya ,t +Dt , zt +Dt = Rq z a ,t +Dt (25)
以上により，変形後の局所座標系が一義的に
で表す．以下では，常に変形後の局所座標系を
決定される．
参照するため，添え字 t+Dt は省略する．
なお，R q は，式 (17) の a に q t+Dt を，f X , f Y ,
{u} に含まれる各軸方向変位は，
fZ に y a,t+Dt と y b,t+Dt とに直交する軸の全体座標系
ub - u a = Dl ,
に対する方向余弦を代入した行列となる．
va = wa = vb = wb = 0 ………(28)
と表される．ここに，Dl：要素の伸びである．
3.3.3　座標変換行列の更新
また，y，z 軸回りの回転角は次式で定義される．
局所座標系の決定に伴い，全体座標系を局所
ya = -
座標系に変換する座標変換行列 T t+Dt は次式のと
Tt +Dt
¶x
, za =
a
¶v
¶x
, yb = a
おり得られる．
¶w
é xt +Dt T ù
ê
ú
T
= ê yt +Dt ú …………………(26)
ê
ú
êz T ú
ë t +Dt û
¶w
¶x
, zb =
b
¶v
¶x
b
………(29)
したがって，a，b 両端における重心軸の接線ベ
クトル x a *，x b * の，変形後の座標系に対する方
向余弦，
要素ごとに得られた行列は，例えば要素剛性
行列の場合，Tt+Dt を用いた次式 (27) により全体
(la , ma , na ) = T x a , (lb , mb , nb ) = T xb ………(30)
座標系の剛性行列 [K] に変換できる．その他の
を用いれば，
T
要素質量，減衰行列も，式 (27) と同様，各行列
ya = -
の左から座標変換行列の転置，右から座標変換
行列をかければ，全体座標系に変換される．
[ K ]t +Dt
éTt +Dt
ê
ê 0
=ê
ê 0
ê
ê 0
ë
0
0
Tt +Dt
0
0
Tt +Dt
0
0
éTt +Dt
ê
ê 0
([ k L ] + [ kG ])t +Dt ê
ê 0
ê
ê 0
ë
ù
ú
0 ú
ú
0 ú
ú
Tt +Dt úû
0
*
0
Tt +Dt
0
0
Tt +Dt
0
0
la
*
, za =
-qa = qb =
T
0
na
T
f
2
T
*
ma
la
, yb = -
nb
lb
T
, zb =
*
mb
lb
,
…………………(31)
と表される．式 (30) 中の T は，式 (26) の座標
変換行列である．
以上のように得られる {u} および変形後の座
標系における要素剛性行列 ([k L ]+[k G ]) を用い，
ù
ú
0 ú
ú
0 ú
ú
Tt +Dt úû
0
反力ベクトル {f} は次式に基づいて計算される．
{ f } = ([ k L ] + [ kG ]){u} ………(32)
開発コードでは，図 3 の❸の段階において，
変形前の要素剛性行列と変位増分ベクトルとの
………(27)
ここに，0 はすべての成分が 0 の 3 行 3 列の行
積から，式 (32) で得られた反力ベクトルを差し
列を示す．
引いて，復元力による残差を算出する．
－9－
©CRIEPI
なお，減衰力も残差を生ずるが，これについ
P
ては，変形前後の要素減衰行列の差に変形後の
U0
d0=L/1000
速度ベクトルすなわち {u} の時間微分を掛け合
W0
わせたものとなる．
L=200
4.　開発コードの機能検証 -1大変形問題の解析機能の検証
開発コードの検証としては，まず，前節に定
図 9　片持ち梁のエラスティカ
式化を示した大変形問題の解析機能の妥当性を
確認するため，既往の文献 10），11） にみられる片
表 3　片持ち梁のエラスティカ解析諸元 *
持ち梁のエラスティカ問題，および放物線固定
アーチの面外座屈問題を解析した．
寸
法
剛
性
座
要
境
初
4.1　片持ち梁のエラスティカ
屈 荷
素分割
界 条
期 不
重
数
件
整
集 中 荷 重P
載 荷 方 法
図 9 に示す，圧縮方向の集中荷重が作用する
片持ち梁の自由端について，水平および鉛直変
長さ L：200
伸び剛性 EA：1 × 106
曲げ剛性 EI：1 × 106
オイラー座屈荷重 Pcr：61.685
20
固定端の全自由度を固定
水平方向変位 d0：L/1000
1 ～ 5 ステップまで増分 10，
6 ～ 55 ステップまで増分 1．
（最終値：P=100）
注 *：諸量はすべて無次元量とした．
位を算出した．表 3 に解析諸元を，図 10 に解析
結果を文献 11） の理論値と比較して示す．なお，
この問題は面内問題であり，面外およびねじり
方向の自由度については検討対象としていない．
図 10 より，計算値は理論値と良く一致してお
り，開発コードの機能は妥当と判断できる．
4.2　放物線固定アーチの面外座屈
3 次元問題として，図 11 の等分布荷重が作用
する放物線固定アーチについて，3 通りに水平
図 10　片持ち梁のエラスティカ解析結果
方向の初期不整を与えた場合の面外座屈問題を
解析し，アーチ頂部節点の面外および鉛直方向
変位を出力した．表 4 に解析諸元を，図 12 に解
析結果を文献
10）
なお，開発コードの結果には文献の結果との
の結果と比較して示す．
乖離，特に図 12(a) において，座屈荷重を上回っ
図 12 より，解析結果と文献の結果は完全には
た段階で座屈の発生がみられる．その原因とし
一致していないが，座屈の発生を含め，荷重に
て，座屈荷重として文献 10）に示された値を用い
対する変位の傾向は，両者の間で共通するもの
たこと，および載荷方法と剛性の相違が考えら
となっている．したがって，3 次元解析におい
れる．文献には，載荷方法に関する詳細な記述
ても開発コードの機能は妥当と判断できる．
が無かったため，ここでは表 4 に示す載荷方法
－10－
©CRIEPI
p
Y
30
X
L=100
Z
図 11　放物線固定アーチの面外座屈
表 4　放物線固定アーチの面外座屈解析諸元 *1
寸
法
剛
性
座 屈 荷 重
要素分割数
境 界 条 件
初 期 不 整
等分布荷重 p*3
載 荷 方 法
支持点間距離 L：100，高さ：30
伸び剛性 EA：4 × 106
曲げ剛性 EIy：1 × 106，EIz：1 × 105
ねじり剛性 GJ：1 × 105
面外弾性座屈荷重 pcr：3.880*2
20
支持点の全自由度を固定
等分布水平荷重 q：
q/pcr=1/10000，1/100，1/10
増分 0.01pcr を 200 ステップで載荷
（最終値：p=2.0pcr）
(a)q/pcr=1/10000 のケース
注 *1：諸量はすべて無次元量とした．
注 *2：固有値問題に対する解として，文献 10）に示さ
れた値．
注 *3：節点集中荷重 qL'/19 または pL/19 に換算し両端
以外の節点に載荷．ここに，L'：アーチ実長．
(b)q/pcr=1/100 のケース
を考えた．また，幾何学的非線形性について，
開発コードの式 (5) に対し，文献 10）では，さら
に高次の変位の非線形項を考慮している．本研
究では，コード化への適用性を勘案し，式 (5)
を採用した．
5.　開発コードの機能検証 -2風車の風応答解析機能の検証
本章では，開発コードを用いた 2 枚および 3
枚ブレード風車の風応答解析結果を示し，挙動，
特にロータの回転運動の模擬が可能であること
(c)q/pcr=1/10 のケース
図 12　放物線固定アーチの面外座屈解析結果
を，既往の風車解析プログラム「F A S T」とのベ
ンチマークにより検証する．解析の対象，条件は，
FAST に添付された例題に基づいて設定した．
－11－
©CRIEPI
5.1　FAST の概要 6），7）
限回転，幾何学的非線形性を考慮する．さらに，
ロータ軸とタワー頂部との交点を二重節点とし，
FAST（Fatigue, Aerodynamics, Structures
その間を 6 自由度を有するばね要素で連結した．
a n d T u r b u l e n c e）は，風車の設計過程におけ
ばね要素の剛性は，全体座標系の X 方向回りを
る断面力の概算を目的として，オレゴン州立大
0，これ以外の 5 自由度の方向を剛として，ロー
学が National Renewable Energy Laboratory
タの回転を可能とした．
(NREL) との連携の下に開発した設計コードであ
なお，上記回転ばねの自由度は，現状では全
る．FAST は，風車専用の解析機能として，
体座標系に固定された状態で設定されるため，
ブレードのアクティブ制御，
ロータ軸の傾きに伴う自由度方向の変化は厳密
ヨーのアクティブ制御，
には考慮されない．
発電機のトルク制御，
ギヤボックス等の駆動系，
HSS ブレーキ，
5.3　計算条件
などのモデル化機能を有し，疲労評価や騒音予
5.3.1　ブレードの空気力係数
測にも適用可能とされている．ただし，モード
ブレードの空力特性としては，F A S T に添付さ
重畳法による解析の時間短縮のため，表 1 に示
れた例題の空気力係数データファイルを用いた．
したとおり，自由度数に制限がある．
表 6 に，ファイル名と設定条件を示す．ブレー
本研究では，F A S T のソースプログラムおよび
ドの断面形状は軸方向に変化するため，要素に
実行形式がフリーウェアとして公開されている
よって異なるファイルと迎え角の初期値が割り
ことから，ベンチマークに用いることとした．
当てられている．表には，各要素の長さと，迎
え角の初期値を併せて示した．初期値の設定に
ついては後述する．
5.2　解析対象
図 15 には代表例として，各ブレードの両端
図 13 に示す，2 枚および 3 枚ブレードの風車
近 傍 の 空 気 力 係 数， す な わ ち 2 枚 ブ レ ー ド の
モデルを解析の対象とした．ここに，全体座標
AWT27_05.dat，AWT27_95.dat，3 枚ブレードの
系の X，Y，Z 方向を図中のとおり定義する．図
S814_1.dat，S813_15.dat の値を示す．FAST で
14 には，開発コードによる有限要素モデルのプ
は，ねじりの自由度を考慮しないため，空力モー
ロッタ出力図を示す．また，ブレード，タワー
メント係数は入力されない．
およびナセルの断面諸元を，開発コードにおけ
る要素分割に対応させて表 5 に示す．表の値は，
5.3.2　その他の計算条件
FAST に添付された例題，TEST01 および TEST08（そ
解析の時間刻み，減衰定数などの計算条件を
れぞれ 2 枚ブレード，3 枚ブレードモデルに対応）
表 7 に ま と め て 示 す． こ れ ら の 条 件 は， 例 題
を参照し，F A S T への入力と等価となるよう設定
T E S T01 および T E S T08 の内容に倣い，開発コー
したものである．ただし，F A S T については，ね
ドと F A S T とに共通するよう設定したものであ
じりの自由度が存在しないため，表中のねじり
るが，両コードの機能の相違に応じて，一部変
剛性は考慮されない．
更されている．特に，現状の開発コードでは，
ロータの回転については，前述のように，大
F A S T と等しく時間刻みおよびブレードの減衰定
変形問題として扱うため，ブレードの要素の有
数を設定した場合，解析精度，収束性が低下し
－12－
©CRIEPI
E
C
C
D
D
E
A
B
A
B
F
F
Z
Z
Z
X
A（ ハ ブ 高 さ ）
B（ タ ワ ー 高 さ ）
C（ ハ ブ 先 端 位 置 ）
D（ ロ ー タ 軸 長 さ ）
E（ ブ レ ー ド 長 さ ）
F（ブレード軸傾き角）
Y
Z
X
42.7m
42.0m
2.7m
2.3m
12.6m
7.0 度
A（ ハ ブ 高 さ ）
B（ タ ワ ー 高 さ ）
C（ ハ ブ 先 端 位 置 ）
D（ ロ ー タ 軸 長 さ ）
E（ ブ レ ー ド 長 さ ）
F（ブレード軸傾き角）
(a)2 枚ブレード風車
Y
25.0m
24.4m
1.3m
1.3m
7.2m
6.0 度
(b)3 枚ブレード風車
図 13　風応答解析の対象
側面 (XZ 面 )　正面 (YZ 面 )　見取図
側面 (XZ 面 )　正面 (YZ 面 )　見取図
(a)2 枚ブレード風車モデル
(b)3 枚ブレード風車モデル
図 14　開発コードによる解析モデルのプロッタ出力図
－13－
©CRIEPI
表 5　解析モデルの諸元
(a) タワー
モデル
要素 No.
2 枚ブレード
1 ～ 21 共通
1( 基礎側 )
2
3
4
5
6
7
8
9
10
11( 頂部側 )
3 枚ブレード
伸び剛性
(N)
8.792 × 1010
1.511 × 1010
1.395 × 1010
1.347 × 1010
1.319 × 1010
1.307 × 1010
1.303 × 1010
1.303 × 1010
1.299 × 1010
1.286 × 1010
1.259 × 1010
1.212 × 1010
曲げ剛性
(Nm2)
1.564 × 1010
2.350 × 109
2.033 × 109
1.817 × 109
1.608 × 109
1.399 × 109
1.190 × 109
9.782 × 108
7.680 × 108
5.573 × 108
3.464 × 108
1.355 × 108
ねじり剛性
(Nm2)
1.203 × 1010
1.808 × 109
1.564 × 109
1.398 × 109
1.237 × 109
1.076 × 109
9.154 × 108
7.523 × 108
5.908 × 108
4.288 × 108
2.665 × 108
1.042 × 108
断面積
(m2)
8.792
1.511
1.395
1.347
1.319
1.307
1.303
1.303
1.299
1.286
1.259
1.212
単位長質量
(kg/m)
879.2
151.1
139.5
134.7
131.9
130.7
130.3
130.3
129.9
128.6
125.9
121.2
注：2，3 枚ブレードモデルとも，要素は等分割，各要素の断面係数は 1.333（両主軸方向共通）．
(b) ブレード
モデル
要素 No.
1( ハブ側 )
2
3
4
5
2 枚ブレード
6
7
8
9
10( 先端側 )
1( ハブ側 )
2
3
4
5
3 枚ブレード
6
7
8
9
10( 先端側 )
伸び剛性 Flap 方向曲げ
(N)
剛性 (Nm2)
9
2.760 × 107
5.850 × 10
9
1.582 × 107
4.911 × 10
9
1.043 × 107
4.891 × 10
9
6.827 × 106
4.274 × 10
9
4.463 × 106
3.602 × 10
9
2.821 × 106
3.027 × 10
9
1.677 × 106
2.500 × 10
9
8.835 × 105
2.093 × 10
9
3.821 × 105
1.620 × 10
8
1.030 × 105
9.536 × 10
9
7.425 × 106
4.448 × 10
9
5.211 × 106
3.154 × 10
9
2.528 × 106
2.194 × 10
9
1.792 × 106
2.191 × 10
9
1.250 × 106
1.835 × 10
9
7.610 × 105
1.641 × 10
9
4.079 × 105
1.441 × 10
9
1.820 × 105
1.191 × 10
9
7.048 × 104
1.074 × 10
8
4.298 × 104
9.560 × 10
Edge 方向曲げ
剛性 (Nm2)
8.618 × 107
9.772 × 107
1.065 × 108
8.982 × 107
6.640 × 107
4.508 × 107
2.751 × 107
1.819 × 107
9.658 × 106
2.780 × 106
8.530 × 106
1.047 × 107
1.134 × 107
1.226 × 107
9.626 × 106
5.224 × 106
3.285 × 106
1.864 × 106
1.055 × 106
7.780 × 105
ねじり剛性
(Nm2)
4.377 × 107
4.365 × 107
4.496 × 107
3.717 × 107
2.725 × 107
1.842 × 107
1.123 × 107
7.335 × 106
3.862 × 106
1.109 × 106
6.138 × 106
6.031 × 106
5.335 × 106
5.404 × 106
4.185 × 106
2.302 × 106
1.420 × 106
7.869 × 105
4.327 × 105
3.158 × 105
断面積
(m2)
5.850 × 10-1
4.911 × 10-1
4.891 × 10-1
4.274 × 10-1
3.602 × 10-1
3.027 × 10-1
2.500 × 10-1
2.093 × 10-1
1.620 × 10-1
9.536 × 10-2
4.448 × 10-1
3.154 × 10-1
2.194 × 10-1
2.191 × 10-1
1.835 × 10-1
1.641 × 10-1
1.441 × 10-1
1.191 × 10-1
1.074 × 10-1
9.560 × 10-2
単位長質量
(kg/m)
58.50
49.11
48.91
42.74
36.02
30.27
25.00
20.93
16.20
9.536
44.48
31.54
21.94
21.91
18.35
16.41
14.41
11.91
10.74
9.560
注：2，3 枚ブレードモデルとも，各要素の断面係数は 1.333（両主軸方向共通）．各モデルのブレード先端節点に，
11.34k g，5.90k g の集中質量を付加．要素分割は，2 枚ブレードモデルは等分割，3 枚ブレードモデルは，N o .1
要素のみ全長の 6.5%，No.2 要素以降はほぼ等分割（後述の表 6 参照）．
(c) ロータ軸・ナセル
モデル
2 枚ブレード
3 枚ブレード
ロータ軸の要素分割数・集中質量 + 慣性モーメント
2 要素，
1330.000kg+335.340kgm2（ハブ），4604.840kgm2（タワー側）．
1 要素，
247.3kg+100.0kgm2（ハブ），1747.0kg+976.3kgm2（タワー側）．
－14－
ナセルの要素分割数・集中質量
1 要素，
5015.430kg（タワー側）
0（ロータ軸と一体化）
©CRIEPI
表 6　ブレードの要素分割と空気力係数の分布
たことから，表中の値に改めた．これについては，
ブ 要 2 枚ブレードモデル
迎え角
レ
素 長さ
｜
ﾌｧｲﾙ名 初期値
ド No. (m)
(度)
1 1.26 AWT27_05 45.00
2 1.26 AWT27_15 40.34
3 1.26 AWT27_25 36.16
4 1.26 AWT27_35 28.94
5 1.26 AWT27_45 20.48
6 1.26 AWT27_55 12.34
7 1.26 AWT27_65 5.66
8 1.26 AWT27_75 1.78
9 1.26 AWT27_85 0.62
10 1.26 AWT27_95 0.23
Rayleigh 減衰により高次モードの減衰が過大と
3 枚ブレードモデル
迎え角
長さ
ﾌｧｲﾙ名 初期値
(m)
(度)
0.47
0.75
0.75
0.75
0.75
0.74
0.76
0.74
0.76
0.74
S814_1
S814_1
S814_1
S814_1
S814_15
S814_15
S812_15
S812_15
S813_15
S813_15
なること，および減衰力の残差などが影響して
いると考えられる．
70.00
45.88
41.87
38.78
35.05
28.67
22.30
15.93
9.56
3.19
また，開発コードでは考慮不可能なパラメー
タについては，F A S T の入力においてすべて無視
あるいは 0 とし，解に影響がないよう配慮した．
こうしたパラメータの主な例を以下に示す．
注：ファイル名の欄には，F A S T の添付ファイルの名
称を拡張子 .dat を省いて示した．
ヨーおよびピッチ制御を無視，
発電機の動作を無視，
風車の初速度を 0 に設定，
ロータ軸のピンまたは弾性支持を無視，
(a)2 枚ブレード
(b)3 枚ブレード
図 15　ブレードの空力特性の代表例
表 7　計算条件のまとめ
風
速
空
気
空 気 力
係
数
時
間
計 算 ス テ
条
件
密
度
ブレード
タ ワ ー
刻
み
ッ プ 数
ブレード
モ ー ド
減衰定数
タ ワ ー
自
重
の
扱
い
開発コード
FAST
12m/s の水平方向の一様流， ダウンウィンド （図 13 の全体座標系 X 軸方向の風向）．
1.225kg/m3 （2 枚）， 0.9526kg/m3 （3 枚）
FAST に添付された例題の空気力係数ファイル使用 （第 5.3.1 項参照）．
考慮不可．
CD=CL=CM=0 の一定値を指定．
0.00025 秒
0.04 秒 (2 枚 )， 0.005 秒 (3 枚 )
200000
12500(2 枚 )， 10000(3 枚 )
無視 （現状ではブレードに減衰を考慮し Edge 方向 1 次 ： 0.114(2 枚 )， 0.04(3 枚 )
た場合， 解析精度が低下．）
Flap 方向 1 次 ： 0.039(2 枚 )， 0.04(3 枚 )
Flap 方向 2 次 ： 0.120(2 枚 )， 0.04(3 枚 )
X 軸方向 1 次 ： 0.05(2 枚 )， 0.03(3 枚 )
1 次 1.385Hz ： 0.05(2 枚 )
X 軸方向 2 次 ： 0.0958(2 枚 )， 0.03(3 枚 )
2 次 8.701Hz ： 0.0958(2 枚 )
Y 軸方向 1 次 ： 0.05(2 枚 )， 0.03(3 枚 )
1 次　3.42Hz ： 0.03(3 枚 )
Y 軸方向 2 次 ： 0.0962(2 枚 )， 0.03(3 枚 )
2 次 17.47Hz ： 0.03(3 枚 )
自重解析を静的解析として実施後， 引き 風応答解析時に自重を考慮．
続き風応答解析を実施．
注：(2 枚 )，(3 枚 ) は，それぞれ 2 枚，3 枚ブレードモデルに対する値であることを示す．
－15－
©CRIEPI
翼端損失を無視，
向方向曲げモーメントおよびロータ軸力につい
風速の鉛直分布（上空逓増）を無視．
ては，両コードの結果は定常状態，すなわちロー
ここに，ロータ軸のピンまたは弾性支持につい
タ回転数が一定値となる 30 秒以降でほぼ一致
ては，F A S T においては 2 枚ブレード風車に対し
しているが，これ以前，特に解析開始後約 10 秒
てのみ指定可能な機能であり，本研究では，2 枚，
間の過渡応答時では，異なる振幅で変動してい
3 枚ブレード両風車に対し，条件を統一するた
る．この原因は，表 7 に示したとおり，両コー
めに無視した．
ドの間で減衰および自重の扱いが異なり，動的
なお，上記に関連し，発電機の作動を無視し
解析の開始直後に作用する外力が一致していな
た場合，ロータの回転に対する負荷が省かれ，
いためと考えられる．特に，開発コードでは，
空力特性を変更しなければ，ロータの回転数が
Rayleigh 減衰により，FAST では考慮されない自
過大となり，解が発散した．このため，通常運
由度，モードに対する減衰が含まれる．
転時として妥当と考えられる回転数となるよう，
表 6 のとおり迎え角の初期値を調整した．これ
5.4.2　3 枚ブレードの結果
らは，変更前の TEST01，TEST08 の値の，それぞ
図 17( a ) ～ ( d ) に，3 枚ブレードモデルの解
れ 7.76 倍，9.10 倍となっている．
析結果として，開発コードによるブレード先端
変位，開発コードと FAST との間のロータの回転
数，タワー基礎の風向方向曲げモーメントおよ
5.4　解析結果
びロータ軸力の比較を示す．前項に示した 2 枚
解析結果より，ブレードの変位によりロータ
ブレードモデルの結果と同様，以下の傾向，
の回転を確認した上で，その回転数，タワー基
・一定風速の下でブレード先端変位が定常的に
礎の風向方向（X 軸方向）曲げモーメントおよ
変位し，ロータの回転を示している．
びロータ軸力を，開発コードと FAST との間で比
・開発コードと F A S T との間でロータの回転数が
較した．
ほぼ一致している．
なお，開発コードは，表 7 に示すとおり，自
・タワー基礎の風向方向曲げモーメントおよび
重と風の作用を二段階に分けて順に解析する仕
ロータ軸力は，開発コードと FAST との間で，過
様となっており，以降の結果は，これらが足し
渡応答の波形は異なるが，定常状態の値はほぼ
合わされた値となっている．
一致している．
がみられる．また，両コードの結果の差異につ
5.4.1　2 枚ブレードの結果
いては，2 枚ブレードの場合と同様の原因が考
図 16( a ) に開発コードによるブレード先端の
えられる．
水平および鉛直変位を示す．図より，一定風速
の下，ブレード先端は等振幅で定常的に変位し
以上のとおり，開発コードと FAST との間には，
ており，ロータが回転状態にあることを確認で
自重の扱いの相違に起因すると考えられる結果
きる．
の差異がみられたが，ロータおよびタワーの挙
図 16(b) には，開発コードおよび FAST による
動に良好な一致がみられ，開発コードにより風
ロータの回転数を比較して示す．図より，両コー
車の風応答挙動が可能であると判断できる．
ドの結果はほぼ一致していることがわかる．
図 16( c ) および ( d ) に示す，タワー基礎の風
－16－
©CRIEPI
(a) ブレード先端変位：開発コード
(a) ブレード先端変位：開発コード
(b) ロータ回転数
(b) ロータ回転数
(c) タワー基礎風向方向曲げモーメント
(c) タワー基礎風向方向曲げモーメント
(d) ロータ軸力
(d) ロータ軸力
図 16　2 枚ブレードモデルの風応答解析結果
図 17　3 枚ブレードモデルの風応答解析結果
－17－
©CRIEPI
6.　試計算による風車の挙動評価例
本章には，開発コードを用いた数ケースの試
計算により，風車挙動を評価した例を示す．こ
こでは，F A S T では考慮不可能なタワーの抗力係
数の影響に着目するとともに，風向の変化とロー
タの回転数との関係を示す．
(a)2 枚ブレードモデル
6.1　タワーの抗力係数の影響
前章の計算条件のうち，0 に設定したタワー
の抗力係数を，0.6，0.8，1.0 に変化させた解
析を実施した．その際，タワー断面の代表長さ（外
径）を 2 枚および 3 枚ブレードモデルとも，高
さ方向に一律に 2.0m とした．これら新たに設定
した抗力係数と代表長さの値は，一般的な風車
の規模およびタワーの断面を円形と想定した場
(b)3 枚ブレードモデル
合，概ね妥当と考えられる．
図 18　タワー基礎風向方向曲げモーメント
解析結果として，タワー基礎の風向方向曲げ
モーメントの時刻歴を図 18 に，その定常状態に
表 8　タワー基礎風向方向曲げモーメント定常値
おける値すなわち時間 30 秒以降の平均値を表 8
曲げモーメント定常値（kNm）
2 枚ブレードモデル 3 枚ブレードモデル
1.0
1014
86
0.8
983
78
0.6
952
70
無視（前章の結果）
858
45
タワー抗力係数
に，それぞれ前章の結果と併せて示す．
図 18，表 8 より，2 枚，3 枚ブレードモデル
の定常状態におけるタワー基礎の風向方向曲げ
モーメントは，抗力係数を無視した場合のそれ
ぞれ 1.18 倍，1.91 倍となっている．
Y
本研究では，F A S T に添付された例題に基づい
て計算条件を設定したが，その値が概ね現実的
であること，また，強風時には風速条件がより
0度
厳しくなり，タワーの抗力係数の違いによる結
果の差異が増大すると考えられる．
-20 度
-10 度
風向
X
10 度
20 度
ロータの回転は
右ねじの方向
なお，ブレードの回転数はタワーの抗力係数
によって変化しない結果となった．
図 19　風向の設定
6.2　風向とロータ回転数との関係
図 19 に示すように，風向を -20，-10，10 およ
び 20 度に変化させた風応答解析を実施した．図
前章の解析条件では，風向を図 13 の X 軸方向，
には，ロータの回転方向を併せて示した．計算
すなわちロータ軸方向に一致させた．ここでは，
条件は，風向以外は前節の抗力係数 1.0 のケー
－18－
©CRIEPI
7.　まとめと今後の展開
本研究で得られた成果を以下に示す．
（1）既往の解析プログラムの多くが，計算の効
率化のためモード重畳法を採用し，自由度数を
制限している現状を鑑み，直接積分法に基づく
(a)2 枚ブレードモデル
有限要素解析コードを開発した．これにより，
自由度数を制限することなく，タワー，ブレー
ド各部の応答変位および断面力の算出を可能と
した．また，この開発コードは，ロータの回転
を含む風車の挙動を精度良く模擬するため，時々
刻々変化しながら作用する風荷重を準定常的に
近似する機能を有し，はり要素の有限回転を考
慮し得るものとした．
(b)3 枚ブレードモデル
（2）開発コードの有限回転および幾何学的非線
図 20　ロータ回転数
形性に関する解析機能について，既往の文献に
みられる，軸圧縮を受ける片持ち梁の座屈，す
表 9　ロータ回転数定常値
ロータ回転数定常値（rpm）
2 枚ブレードモデル 3 枚ブレードモデル
10
89.0
82.7
-10
89.0
82.7
20
87.0
80.6
-20
87.0
80.6
0（前節の結果）
90.0
83.3
なわちエラスティカ問題および固定放物線アー
風向（度）
チの面外座屈問題を解析し，その妥当性を示し
た．エラスティカ問題については，開発コード
による解を解析解と比較し，両者の一致を確認
した．また，固定放物線アーチの面外座屈問題
については，開発コードによる解が文献にみら
スと同様とした．
れる解と概ね一致することを確認した．両者の
解析結果として，図 20 にロータ回転数の時刻
間の差異は，開発コードについて計算条件を仮
歴を，表 9 に定常状態におけるロータ回転数と
定せざるを得なかったこと，および幾何剛性の
風向との関係を示す．
定式化が異なることによると考えられる．
図 20，表 9 より，2 枚，3 枚ブレードモデル
ともに，風向の絶対値の増加に伴い，ロータ回
（3）開発コードの風応答解析機能を検証する
転数は低下する傾向にあることがわかる．また，
ため，一様流が作用する 2 枚および 3 枚ブレー
本解析の条件下では，風向の正負とは無関係に，
ド風車の挙動を解析し，既往の解析プログラム
絶対値が等しい場合，ロータの回転数は一致す
FAST との間で解を比較した．その結果，開発コー
ることがわかる．
ドによって算出されたブレード先端の挙動には，
ロータの回転を示す定常的な変位がみられ，定
常状態におけるロータの回転数は，F A S T によ
る解と一致した．さらに，両コードの離散化手
－19－
©CRIEPI
法などの相違により過渡応答波形が異なるもの
参考文献
の，タワー基礎の風向方向曲げモーメントおよ
びロータの軸力は開発コードと FAST との間で概
1）独立行政法人新エネルギー・産業技術総合開
ね一致する結果となった．以上により，開発コー
発機構：日本型風力発電ガイドライン，http://
ドにより，ロータの回転を含む風車の風応答挙
www.nedo.go.jp/library/furyokuhoukoku/
動の解析が可能であることが示された．
index.html，2008.
2）土木学会：風力発電設備支持物構造設計指針・
今後，開発コードの特質，すなわち，ねじり
同解説 [2007 年版 ]，2007.
を含め自由度数を制限することなく，風車各部
3）例えば，独立行政法人新エネルギー・産業技
の変位および断面力を算出し得る機能を，設計
術総合開発機構：平成 18 年度風力発電利用率向
時の詳細な動的照査あるいは事故時の挙動分析
上調査委員会および故障・事故等調査委員会報
などに適用し，コードの信頼性向上を図る．
告書，2007．
その一方で，現状の開発コードには，
4）http://www.garradhassan.com/products/
・ブレードの死荷重，減衰を考慮した場合，解
ghbladed/index.php
の収束性，精度が低下する．
5）日本風力エネルギー協会：平成 17 年度第一
・ロータ軸とタワー頂部との二重節点を連結す
回風力発電システム技術講習会資料　講習 2　る回転ばねの自由度を全体座標系で指定するた
風車の機械要素設計，2005.
め，ロータ軸の傾きに伴う自由度方向の変化は
6）MOLENAAR, D.,“Cost-effective design and
厳密には考慮されない．
operation of variable speed wind turbines,”
といった課題がある．したがって，これらを改
PhD. Thesis, Delft Univ. Tech., Inst. for
善するため，さらなる開発が必要と考えられる．
Wind Energy, 2003.
また，空気力のモデル化において，翼素運動
7）h t t p : / / w i n d . n r e l . g o v / d e s i g n c o d e s /
量理論を反映するためのパラメータの導入も，
simulators/fast
併せて必要と考えられる．
8）牛山泉：風車工学入門，森北出版株式会社，
2002.
9）鷲津久一郎，宮本博，山田嘉昭，山本善之，
謝辞
川井忠彦共編：有限要素法ハンドブック II　応
本研究を実施するにあたり，多大なご助力を
用編，培風館，pp.128-149，1990.
下さいました，株式会社アーク情報システム先
10）前田幸雄，林正：立体骨組構造物の有限変
端技術センター数理解析部の，佐藤順一部長，
位解析，土木学会論文報告集，第 253 号，p p .
坂口剛氏，電力中央研究所地球工学研究所流体
13-27，1976.
科学領域の松宮央登氏に深甚なる感謝の意を表
11）Timoshenko, S. P. and J. M. Gere,
します．
Theory of Elastic Stability, McGraw-Hill,
pp.76-82, 1961.
－20－
©CRIEPI
電力中央研究所報告
〔不許複製〕
編集・発行人
e-mail
財団法人
電力中央研究所
地球工学研究所
千葉県我孫子市我孫子 1646
電 話 04 (7182) 1181 (代 )
[email protected]
発 行 所
財団法人 電力中央研究所
東京都千代田区大手町 1-6-1
電 話 03 (3201) 6601 (代 )
印 刷 所
株式会社 ユウワビジネス
東京都千代田区神田須田町 1-1
電話
03（ 3258） 9380
ISBN978-4-86216-888-7