みずほリサーチ＆テクノロジーズ技報 Vol.2 No.1

浮体式洋上風力発電に関わる数値シミュレーションの利用状況と MHRT の取り組み
坂本大樹i，吉村英人ii，眞鍋尚iii

Utilization of Numerical Simulation and Our Solution for Floating
Offshore Wind Turbine
Daiki SAKAMOTO, Hideto YOSHIMURA, Takashi MANABE
近年，脱炭素化社会の実現に向けて，浮体式洋上風力発電の開発が進められてきており，浮体に関するシミュレーション技術にも注目が高まっている．本稿では，浮体式洋上風車の設計に関わるシミュレーシ
ョンについての基礎理論をまとめたうえで，当該分野に関連する当部のソリューションである「船舶・海洋構造物動揺解析システム MIZUHO_MARIS-Ⅱ」を紹介する．
(キーワード): 浮体式洋上風力発電，動揺解析，三次元特異点分布法

1 はじめに

2 数値シミュレーションの適用事例
2.1 浮体式洋上風車のシミュレーションの流れ

世界的に脱炭素化社会の実現が叫ばれる現代において，我が国でも 2020 年に策定された「2050 年カー

浮体式洋上風車には風荷重や波・潮流の外力によ

ボンニュートラルに伴うグリーン成長戦略」では，電り浮体動揺が発生する．また浮体形状によっては，細力部門の脱炭素化の柱の一つとして，洋上風力発電

い部材が用いられることもあるため，波などの外力が取り上げられており，2030 年 10GW，2040 年 30～による弾性変形を考慮する必要も考えられる．そし45GW の案件形成を導入目標としている．て風車特有の問題として発電量を安定化させるため洋上風力発電の中でも，遠浅な海岸が少ない日本に行われるピッチ角制御もあり，これらの問題を一では，浮体式洋上風力のポテンシャルが大きいとさ体として解く連成解析が必要である．NEDO が発行れているが，着床式と比べて浮体式の開発は遅れてしている浮体式洋上風力ガイドブック 2)にもこれらおり，北九州市沖などで実証試験が行われている段に関する記載があり，解析フローを図 1 に示す．次階である．

節以降では各解析について簡単にまとめる．

今後はより事業化への流れが加速すると考えられ，施設の設計の際には安全基準を満たすためにシミュ

2.2 浮体動揺についてレーション技術が活用されるため，浮体に関するシ2.2.1

ミュレーション技術にもより注目が高まっていくと運動方程式

一般的な浮体動揺の解析では以下の運動方程式を考えられる．本稿では，そのような浮体式洋上風力発解くことになる．

電の設計に関連する数値シミュレーション技術を整𝑀𝑥̈ + 𝐶𝑥̇ + 𝐾𝑥 = 𝐹 + 𝐹 + 𝐹 + 𝐹

理したうえで，関連のある MHRT の取り組みについて紹介する．

i

サイエンスソリューション部

先進技術システムチーム

ii

サイエンスソリューション部

先進技術システムチーム 上席主任コンサルタント

iii

サイエンスソリューション部

先進技術システムチーム 次長 博士（工学）主任コンサルタント

1

博士（工学）(1)

みずほリサーチ＆テクノロジーズ技報 Vol.2 No.1

以下ではその流れを簡単にまとめる．
まず翼素理論では，各ブレードを複数の翼素に分解し，各翼素でのスラスト係数，接線力係数を評価す
る．各翼素でのスラスト力𝑑𝑇は式(4)のように表される．

𝑑𝑇 = 0.5𝜌𝑉

𝐵𝑐(𝐶 𝑐𝑜𝑠𝜙 + 𝐶 𝑠𝑖𝑛𝜙)𝑑𝑟(4)

ここで，𝐵：ブレード数，𝐶 ：揚力係数，𝐶 ：抗力係数，𝑟：翼半径，𝑐：翼弦長である．
上式を無次元化したスラスト係数𝐶 は式(5)のように表される．

𝜎(1 − 𝑎)(𝐶 𝑐𝑜𝑠𝜙 + 𝐶 𝑠𝑖𝑛𝜙)
(5)
𝑠𝑖𝑛 𝜙
ここで，𝑀：質量マトリクス，𝐶：減衰マトリクス，また同様にして接線力𝑑𝑄 は式(6)のように表される．
𝐾：剛性マトリクスであり，𝐹 ，𝐹 ，𝐹 ，𝐹 はそれぞれ風荷重，流体力，係留力，浮力を表している．以下
𝑑𝑄 = 0.5𝜌𝑉 𝐵𝑐(𝐶 𝑠𝑖𝑛𝜙 − 𝐶 𝑐𝑜𝑠𝜙)𝑟𝑑𝑟(6)では𝐹 ，𝐹 ，𝐹 といった各外力解析方法について記
そして式(6)を無次元化した接線力係数𝐶 は式(7)の𝐶 =図 1 連成解析のフロー載する．
2.2.2ように表される．

𝜎(1 − 𝑎)(𝐶 𝑠𝑖𝑛𝜙 − 𝐶 𝑐𝑜𝑠𝜙)
(7)
𝑠𝑖𝑛 𝜙
なお𝜎はソリディティであり式(8)で表される．
𝐶 =

空力解析

空力解析の手法としては，翼素・運動量理論(BEM)
や P＆P モデル，一般化動的後流モデルなどが存在し，

𝜎 = 𝑐𝑁 2𝜋𝑟

洋上風力まわりの数値解析では BEM が用いられることが多い．

(8)

一方，運動量理論ではブレード前後の運動量の変BEM はブレード翼素の揚力係数と抗力係数からス化量からスラストと接線力を求める方法である．力

ラストと回転力を算出する翼素理論と風がロータ平𝑑𝑇はブレード前後の運動量の変化量と等しいため，面を通る前後での運動量の変化からスラストと回転式(9)のように表される．

力を推定する運動量理論を組み合わせた理論であり，

𝑑𝑇 = 𝜌𝐴𝑉(𝑉 − 𝑉 )簡易的な手法であるため風車の設計に広く用いられている．

(9)

ここで𝐴：ロータの断面積であり，まず空力解析の基本であるが，軸流速度の誘導係

𝐴 = 2𝜋𝑟𝑑𝑟

数𝑎を用いると，風車上流の一様流速𝑉 に対して，風(10)車のロータ位置における軸方向速度𝑉 ，ロータから十である．式(9)を用いつつ無次元化したスラスト係数

分に下流における軸方向速度𝑉 はそれぞれ式(2) ，(3)𝐶 を求めると式(11)のようになる．のように表すことができる．

𝐶 = 4𝑎(1 − 𝑎)

𝑉 = 𝑉 (1 − 𝑎)(2)

𝑉 = 𝑉 (1 − 2𝑎)

(3)

(11)

また環状部分に作用する接線力は，後流が得る回転モーメントに等しいため式(12)のように表される．

𝑑𝑄 = 𝜌𝑉𝐴 × 𝑟𝜔 × 𝑟

誘導係数𝑎は正の値をとるため，ロータの下流におけ(12)る風速は，上流の風速と比べて小さい値となる．このよって無次元化した接線力係数𝐶 は式(13)のようにような領域は風車ウエイクと呼ばれ，洋上風力のみ表される．

ならず風力発電全般にて，風車を設置する際には検𝐶 = 4𝑎′(1 − 𝑎)𝜆討が必要となる現象である．
BEM では各タイムステップにおいて，この誘導係(13)

なお𝑎′は角誘導係数であり式(14)のように表される．

数𝑎を求めたうえでスラストや回転力を算出する．
2

みずほリサーチ＆テクノロジーズ技報 Vol.2 No.1

𝑎 =𝜔2Ω

また軸方向流体力についてはモリソン式では考慮し(14)ないことが多い．

翼素理論・運動量理論にてそれぞれ求めたスラス

セミサブ型などの細い部材で構成される浮体につト係数・接線力係数が等しいものとすることで，誘導いては，波の変形解析の影響は小さく，粘性抵抗によ係数𝑎を求めたうえで，軸流速度や風車後流速度をる非線形抗力が支配的であるため，浮体の弾性変形求めることができる．また，翼端渦などによる抵抗をや共振現象を考慮することができるモリソン式を適考慮する場合は，それぞれの抵抗係数を求め，運動方用することが多い．

程式に組み込むことで表現する．以上が BEM の概要

上記のように浮体形状によって適する手法が異なである．

り，各手法の特徴をまとめると表 1 のようになる．

BEM の計算には計算対象とする風場の条件に合わ

なお浮体の弾性変形を考慮した動揺計算手法の検討

せた揚力係数や抗力係数を用いることが望ましいが， もなされており，最初から浮体を弾性体として評価
現実的には代表的なレイノルズ数における値のみをする弾性応答解析と各タイムステップで，最初は浮使用することが多く，計算精度に課題が残る．この問

体を剛体として各外力を求め，その後浮体を弾性体

題点を解決するために，より発展した計算方法であとし，求めた外力を与えて応力を求める 2 段階解析

る Generalized Dynamic Wake Model や CFD 計算を用の 2 つが存在する．詳しくは石原ら(2008)6)を参照さいた研究も存在する．また，風車の運転停止時のようれたい．

な場合は，風車が受ける外力𝐹 を抗力係数𝐶 と受風

表 1 浮体形状と解析手法

面積𝐴′を用いて，単純に式(15)のように評価することも考えられる．

𝐹 =
2.2.3

1
𝜌𝐶 𝑉𝐴′
2(15)流体解析

流体解析は大きく分けてポテンシャル理論に基づくものとモリソン式に基づくものに分けられる．そ2.2.4

れぞれの特徴を以下にまとめる．

2 種類に分けることができ，それぞれについて以下に① ポテンシャル理論設定した浮体形状まわりのポテンシャルを算出し，造波減衰と波強制力を求める手法である．部材によまとめる．
① 準静的解析る波の変形解析効果を考慮することができるため，準静的解析では設定した係留索に対して歪率と荷バージ型といった大型の浮体形状に対して用いられ重からなる係留索特性曲線を与え，各時刻においてることが多い．また浮体を剛体と仮定し任意形状に設定した係留索特性を参照し，係留ラインの上端点対応可能である．詳細は 3.2.3 に記載する．とアンカー点から係留索張力を求め運動方程式に反映させる手法のことである．

② モリソン式
モリソン式を用いると浮体に働く軸直交方向の流② 動的解析

体力は式(16)のように表される．

𝜋𝐷𝑓 = 𝑚𝑢 ̇ + 𝐶 , 𝜌(𝑢 ̇ − 𝑥 ̈ )4
1+ 𝐶 , 𝜌𝐷|𝑢 − 𝑥 ̇ |(𝑢2− 𝑥̇ )係留解析

係留解析は大きく分けて準静的解析と動的解析の一方で，動的解析は準静的解析に係留ラインに作用する慣性力・減衰力を考慮した解析方法であり，例
えば動的解析手法の一つであるランプドマス法では，(16)

係留索を複数のマスとバネ要素に分解し，要素ごとに運動方程式を解くものとなっている．係留索に作用する外力や係留索自体のダイナミクスを考慮でき

ここで，𝑚：部材の質量，𝜌：流体密度，𝑢 ：流体速度の部材直交成分，𝐶 , ：軸直交方向の付加質量係数， る一方で，計算が不安定化したり，計算時間が長期化
𝐶 , ：軸直交方向の抗力係数，𝐷：部材の直径である． したりする傾向が強い．また係留索をマスとバネで3

みずほリサーチ＆テクノロジーズ技報 Vol.2 No.1
モデル化しているため，各要素の回転方向の変位は細い部材で構成され各部材に働く応力を検討する際考慮できないという特徴がある．に適していると言える．部材を梁で表現した研究事国の指針では，「解析の対象に応じて準静的解析又例は存在するが，シェルやソリッドとして評価したは管海官庁が適当と認める動的解析を行い，浮体施計算事例はあまり見当たらない．

設の最大変位量及び最大ライン張力を計算すること」とあるように明確にどちらを用いるかの指定はない． 2.4 制御解析について
ただし一般的に準静的解析と比べると動的解析を行

風車は風況に適した発電や安全性の確保のためにった方が係留索張力は大きくなりやすいため，準静自動で制御を行っている．制御はトルク制御，ブレー的解析と動的解析にて用いる安全率を変えることで， ドピッチ制御，ヨー制御の 3 種類に分けることがでその差分を評価している．なお「ISO19901-7」では表き，それぞれについて簡単にまとめる．

2 に示すように推奨解析方法が記載されている他，ABS などの合成繊維ロープを用いた係留システムで① トルク制御は動的解析を求める船級協会も存在する．このようまず発電機回転数が定格以下の場合は，発電効率な現状を踏まえ，国の指針では準静的解析を基本とが最大となる周速比となるように発電機のトルク Qしつつ．対象となる浮体施設ごとに管海官庁が判断を式(17)に従って制御する．またこの場合ブレードピすることとしている．

7)ッチ角は 0 である．ここで周速比は風速と翼端速度

表 2 係留タイプによる推奨解析手法のまとめ
係留タイプ

限界状態 解析すべき状態

解析法

永久係留

終局

非損傷/破断動的

テンポラリー係留の比であり，風車の効率は周速比の関数で示される．

2)

𝑄=𝑘×Ω

(17)

過渡準静的又は動的

疲労

非損傷

動的

終局

非損傷/破断準静的又は動的

過渡準静的又は動的

なお，Ω は発電機回転数，𝜌は空気密度，𝑅はロータ

ーー半径，𝐶疲労

𝑘=𝜋𝜌𝑅 𝐶
2𝑟 𝜆𝜂

(18)

は最適収束比，𝜂 は発電機効率，𝑟は増速機 のギア比である．
2.3 構造解析について
構造解析の手法は大きく分けて以下の 2 つに分け② ブレードピッチ制御

ることができる．

発電機回転数が定格以上の場合，発電機出力が一定となるようにブレードピッチ角を調整する．これ① マルチボディダイナミクス(MBS)をブレードピッチ制御と呼ぶ．この場合は出力が定MBS は物体を複数の剛体要素に分割し，それぞれ格を維持するため，発電機のトルク𝑄′は式(19)で示すの物体について運動方程式をたてる．そして剛体要指示値とする．なお𝑃 は定格出力である．

素間の拘束条件を踏まえたうえで運動方程式を解く

𝑄′ =ことで，物体全体の運動を解析する方法のことであ𝑃Ω(19)

る．浮体式洋上風力の場合，風車構造をブレード，タまたブレードピッチ角は発電機の回転数が一定となワー，浮体基礎などそれぞれの要素ごとに剛体とみるように式(20)に示す PI 制御によって制御を行う．

なし，拘束は回転ジョイントや剛結で表現すること𝜃で風車全体の運動を解析することができる．この際= 𝜅(𝐾 𝑒(𝑡) + 𝐾 𝑢 (𝑡))(20)ジョイント部分には拘束力として反力が生じる．上

ここで𝑒(𝑡)，𝑢 (𝑡)は発電機回転数の偏差とその積分

述のとおりすべての部材を剛体とするため，バージ値，比例ゲイン𝐾 と積分ゲイン𝐾 は Hansen et al4)に型などの大型の浮体の解析に適している．

従い式(21) ，(22)により求める．また，𝜅はゲインス
ケジューリング関数であり，式(23)を用いる．

② 有限要素法(FEM)

𝐾 =構造物を複数の要素に分解し，各要素の変位や応力を求めるもの．支持構造物の任意断面における断面力の照査が可能となるため，セミサブ型といった
42𝐼 Ω 𝜉 𝜔
𝜕𝑃𝑟(− )
𝜕𝜃(21)

みずほリサーチ＆テクノロジーズ技報 Vol.2 No.12𝐼 Ω 𝜔𝜕𝑃𝑟(− )
𝜕𝜃1
𝜅=
𝜃1+𝜃

𝐾 =(22)(23)

ここで，𝐼 はロータ，増速機，発電機の回転子の慣性モーメントの合計であり，𝜉 は制御系の減衰，𝜔は制御系の周期である．
なお，一般に実風車の制御アルゴリズムは，風車の図 2性能に直結するため，各風車メーカーは公開していMARIS における浮体動揺解析フローステムであり，計算条件の設定～解析～解析結果のないことが多いが，ここでは山口ら(2018)8)がまとめ可視化を全て同一のシステム上で行うことができる．

た標準的な風車モデルについて紹介している．また制御の与え方によっては，浮体の動揺量が増幅する

3.2.2

ネガティブダンピングと呼ばれる現象が発生するた浮体形状の入力

まずは解析を行う浮体形状を設定する．浮体形状め注意する必要がある．

は簡易的な箱型や船型形状が設定可能である．また
節点と要素を設定することで任意形状の浮体を設定③ ヨー制御することも可能である．

実際の風力発電では時々刻々と風向は変動する．
風車が無駄なく風を受けるために，風車の向きを風3.2.3

向に追従させる制御のことをヨー制御という．

流体力波力の算出設定された浮体形状に対して，流体力波力の算定を行う．MARIS ではポテンシャル理論に基づく 3 次

3 MHRT の取り組み

元特異点分布法を用いて流体力波力の算定を行っており，以下に簡単な流れを記載する．

本分野に関連する MHRT の取り組みとして，「船舶・海洋構造物動揺解析システム MIZUHO_MARIS-① 基礎方程式と境界条件Ⅱ（以下 MARIS と記載）」を用いた浮体動揺解析が挙一般にディフラクションポテンシャル関数𝜙 と浮げられ，以下に紹介する．

体の動揺を𝑖モードとした際のラディエーションポテンシャル関数𝜙 とすると，流体領域Ωにおける基礎

3.1 MARIS の概要

方程式と，自由表面𝑆 ，水底𝑆 ，浮体表面𝑆 ，原点かMARIS は 3 次元特異点分布法を用いて流体力，波らの水平距離を 𝑅 としたときの無限遠方での各境界力を求め，浮体係留構造物の動揺を解析するシステ

条件は以下のように表される．

ムであり，MHRT では MARIS の開発・販売や同システムを利用した受託解析を行っている．これまで沈・ディフラクション

埋トンネル・ケーソンの設置に関する安全性評価，係

𝜕 𝜙
𝜕 𝜙
𝜕 𝜙++= 0 (in Ω)𝜕𝑥𝜕𝑦𝜕𝑧
𝜕𝜙𝜙 𝜔−= 0 (𝑜𝑛 𝑆 )𝜕𝑧𝑔
𝜕𝜙= 0 (𝑜𝑛 𝑆 )𝜕𝑧
𝜕𝜙= 𝑣 (𝑜𝑛 𝑆 )𝜕𝑛𝜕𝜙lim √𝑅− 𝑖𝑘 𝜙 = 0→
𝜕𝑅(when 𝑅 → ∞)留船の安定性評価など，海洋構造物の安定性，安全性評価など幅広く活用されてきている．
3.2 解析機能
3.2.1

解析の流れ

解析の流れとしては，図 2 に示すとおり，まずは浮体の形状を入力し，計算に必要なパラメータを設置する．その後流体力・波力を算出し，算出された流体力波力を用いて時系列シミュレーションを行う形となっている．各項目について詳細を以下に記載する．なお MARIS は GUI・ソルバーが一体となったシ5(24)みずほリサーチ＆テクノロジーズ技報 Vol.2 No.1・ラディエーション

ついても同様にして求めることができる．
なおここで用いた Green 関数𝐺(𝑃, 𝑄)は，𝑃, 𝑄の座標の𝜕 𝜙𝜕 𝜙𝜕 𝜙++= 0 (in Ω)𝜕𝑥𝜕𝑦𝜕𝑧

他に水深ℎ，波数𝑘，波周波数に依存する関数であり，級数形と積分形の 2 種類が存在する．それぞれのGreen 関数を以下に示す．

𝜕𝜙𝜙 𝜔−= 0 (𝑜𝑛 𝑆 )𝜕𝑧𝑔

・級数形 Green 関数

𝜕𝜙= 0 (𝑜𝑛 𝑆 )𝜕𝑧𝐺(𝑃, 𝑄) =(25)∙ [𝑌 (𝑘𝑅) − 𝑖𝐽 (𝑘𝑅)] + 4

𝜕𝜙= 𝑣 (𝑜𝑛 𝑆 )𝜕𝑛lim √𝑅→𝜕𝜙− 𝑖𝑘 𝜙
𝜕𝑅(𝑘+𝑣 )𝑘 ℎ−𝑣 ℎ+𝑣

(30)

cos[𝑘 (𝑧 + ℎ)] cos[𝑘 (𝑧 + ℎ)] ∙ 𝐾 (𝑘 𝑅)

・積分形 Green 関数=0𝐺(𝑃, 𝑄) =(when 𝑅 → ∞)(𝑗 = 1,2, ⋯ ,6)(𝜇 + 𝑣)𝑒

ここで，𝑛 は浮体の j 方向に対する単位法線ベクトルであり，𝑣 は浮体表面の条件として与えられる複素流+2𝜋速を表す．また，𝑣 は入射ポテンシャルの法線方向微𝜕𝜙𝜕𝑛

1 1+ + 2𝑃. 𝑉.
𝑅 𝑅′

cosh[𝜇(𝜉 + ℎ)] cosh[𝜇(𝑧 + ℎ)]∙ [𝐽 (𝜇𝑅)]𝑑𝜇(𝜇𝑠𝑖𝑛ℎ𝜇𝑘 − 𝑣𝑐𝑜𝑠ℎ𝜇𝑘)(31)(𝑘 − 𝑣 ) cosh[𝑘(𝜉 + ℎ)] cosh[𝑘(𝑧 + ℎ)]∙ [𝐽 (𝑘)]𝑑𝜇𝑘 ℎ−𝑣 ℎ+𝑣

ここで𝑃. 𝑉.は主値積分を表し，𝑘 ：n 次の波数，𝐽 ：分として以下のように表される．

𝑣 =2𝜋(𝑣 − ℎ )𝑐𝑜𝑠ℎ[𝑘(𝜉 + ℎ)]𝑘 ℎ−𝑣 ℎ+𝑣0 次の J-Bessel 関数，𝑌 ：0 次の y-Bessel 関数，𝐾 ：(26)0 次の変形 Bessel 関数である．
級数系 Green 関数の場合，収束が速いが，P と Q に一致する場合には発散する可能性があるという特徴② 3 次元特異点分布法がある．MARIS では，P と Q の距離によって 2 種類

3 次元特異点分布法では，式(24)，式(25)を満足すの Green 関数を使い分けることとしている．

る Green 関数を定義し，支配方程式を境界条件の下で解くことでディフラクションポテンシャル，ラデ③ 流体力・波力の算出

ィエーションポテンシャルを求める．Green 関数を𝐺(𝑃, 𝑄)とすると任意の点におけるディフラクション

流体中に存在する構造物には，波の作用による波ポテンシャル𝜙 ，ラディエーションポテンシャル𝜙

圧・波力と，構造物が動揺した際に流体から受ける動は以下のように表される．

水圧が働く．波圧は構造物が動かない仮定のもとで14𝜋1𝜙 (𝑃) =4𝜋𝜙 (𝑃) =𝜎 (𝑄)𝐺(𝑃, 𝑄)𝑑𝑆

(27)

𝜎 (𝑄)𝐺(𝑃, 𝑄)𝑑𝑆

(28)

入射波の散乱波すなわちディフラクションポテンシ
ャルと既知である入射波のポテンシャルの和として求められ，波力を構造物表面で積分することで波力を求めることができる．一方で動水圧は付加質量とここで，𝑃 = (𝑥, 𝑦, 𝑧)：流体中の任意点，𝑄 = (𝜉, 𝜂, 𝜍)：造波減衰係数として表現されることが多く，MARIS

湧き出し点，𝜎(𝑄)：湧き出し強さ，𝑆：物体表面である．上式から𝜎 (𝑃)に関する以下の積分方程式が得らでも同様である．また付加質量力と造波減衰力を合わせてラディエーション流体力と呼んでいる．
まずラディエーション流体力であるが，波周波数れる．

𝜕𝜙 (𝑃)𝜎 (𝑃)=−𝜕𝑛2
1+4𝜋に依存し，付加質量力は構造物の動揺する加速度，造波減衰力は速度にそれぞれ比例するものとなってお(29)り，流体力𝐹 は付加質量係数𝑀 ，造波減衰係数𝑁 を𝜎 (𝑄)𝐺(𝑃, 𝑄)𝑑𝑆用いて以下のように表される．

これを解いて𝜎 を求め，式(27)に代入することで流体𝐹 = −𝑀 𝑋̈ − 𝑁 𝑋̇

内の任意の点におけるディフラクションポテンシャ
ルが求まる．またラディエーションポテンシャルに(32)

一方で流体力𝐹 はラディエーションポテンシャル𝜙6

みずほリサーチ＆テクノロジーズ技報 Vol.2 No.1
を用いて以下のように表される．

可能であり，それぞれの運動方程式を以下に記載する．

𝐹 = 𝑅𝑒 𝜔 𝜌𝑋𝜙𝑟𝑗 𝑛𝑖 𝑑𝑠 𝑒𝑥𝑝



周波数一定型運動方程式(33)

𝑀 (𝜎) + 𝑚 𝑥̈ + 𝐶 (𝜎)𝑥̇ + 𝐾 𝑥= 𝜌 ∫ 𝑅𝑒(𝜙𝑟𝑗 )𝑛𝑖 𝑋̈ + 𝜔 𝐼𝑚(𝜙𝑟𝑗 )𝑛𝑖 𝑋̇ 𝑑𝑠(38)= 𝐹 (𝑡)

よって，付加質量係数𝑀 ，造波減衰係数𝑁 は以下の遅延関数型運動方程式ように表される．

𝑀 =𝜌𝑅𝑒(𝜙𝑟𝑗 )𝑛𝑖 𝑑𝑠

𝑁 = 𝜌𝜔𝐼𝑚(𝜙𝑟𝑗 )𝑛𝑖 𝑑𝑠

𝑀 (∞) + 𝑚 𝑥̈ +(34)+ 𝐶 (𝜎)𝑥̇ + 𝐾 𝑥 = 𝐹 (𝑡)1
𝑀 (∞) = 𝑀 (𝜎) +𝑁 (𝑡)𝑠𝑖𝑛𝜎𝑑𝑡
𝜎2
𝐿 (𝑡) =𝑁 (𝜎)𝑐𝑜𝑠𝜏𝑑𝑡𝜋

次に波力であるが，波圧𝑃を物体表面で積分して以下のように表すことができる．

𝐹 (𝑡) = −(39)

(40)
(41)

周波数一定型運動方程式が特定の周波数𝜎について𝑃(𝑡) 𝑛𝑗 𝑑𝑠の運動方程式となっているのに対して，遅延関数型運動方程式は不変付加質量𝑀 (∞)や遅延関数𝐿 (𝑡)𝑆(35)= 𝑖𝜌𝜔 cos(−𝜔𝑡 + 𝛽)

𝐿 (𝑡 − 𝜏)𝑥̇ 𝑑𝜏

𝜙 𝑛𝑗 𝑑𝑠を用いることで，全周波数に対する運動方程式とな𝑆っている．

ここでディフラクションポテンシャルは各要素において得られるため，上式の積分は次のようになる．
𝑚𝐹 (𝑡) = 𝑖𝜌𝜔 cos(−𝜔𝑡 + 𝛽)𝑖=13.3 解析事例

(36)
𝛥𝑆

ここでは MARIS における解析事例を紹介する．

𝜙𝑖 𝑛𝑖𝑗 𝛥𝑑𝑆𝑖

3.3.1

ここで，𝑛 は𝑖要素に関する𝑗方向の法線ベクトルでMARIS では任意の 3 次元形状の浮体に対して解析あり，以下のように表される．

𝑛, =𝑛,𝑛, =𝑛,𝑛, =𝑛,𝑛 , = (𝑦 − 𝑦 )𝑛 , − (𝑧 − 𝑧 )𝑛 ,𝑛 , = (𝑧 − 𝑧 )𝑛 , − (𝑥 − 𝑥 )𝑛 ,𝑛 , = (𝑥 − 𝑥 )𝑛 , − (𝑦 − 𝑦 )𝑛 ,形状作成並びにその他の条件設定を行うことができる．図 3 にバージ型，セミサブ型の浮体形状についてメッシュを作成した例を示す．
なおメッシュは各浮体の喫水面以下について作成し(37)

た．
3.3.2

流体力波力並びに時系列計算

ここで，(𝑛 , , 𝑛 , , 𝑛 , )は𝑖要素の外向法線ベクトル，設定した浮体形状，計算条件に対して解析を行う．

(𝑥 , 𝑦 , 𝑧 )は𝑖要素の中心の座標，(𝑥 , 𝑦 , 𝑧 )は構造物

計算開始ボタンをクリックすると実行ウインドウがの重心座標である．

立ち上がり，計算が実行される．計算条件設定画面と式(36)によって特定の周波数𝜔の規則波における計算実行時の実行ウインドウの例を図 4 に示す．

波力が求まることになる．実際には，浮体に作用する波は不規則波であるため，周波数スペクトルや方向3.3.3

計算結果の可視化

分布関数を使用して不規則波を生成し，生成した波計算が正常終了すると，計算結果を GUI 上から確に対して波力を求めている．詳しくは眞鍋(2003)9)を認することができる．計算結果としては，浮体動揺量参照されたい．

や設定した係留索の張力，防舷材の反力などについて時刻歴解析結果や周波数応答図が確認できる．計

3.2.4

時系列シミュレーション算結果の表示例を図 5 に示す．

上で求めた流体力波力を用いて浮体動揺の時系列シミュレーションを行う．解く運動方程式は周波数一定型運動方程式と遅延関数型運動方程式から選択

7

みずほリサーチ＆テクノロジーズ技報 Vol.2 No.1

図 5MARIS 計算結果表示例

引 用 文 献1) 国土交通省海事局：浮体式洋上風力発電施設技術
基準安全ガイドライン，(2020)

図 3 浮体式洋上風車メッシュ作成図2) 新エネルギー・産業技術総合開発機構：浮体式洋（上：バージ型 下：セミサブ型）

上風力発電技術ガイドブック，(2018)3) P.J.Moriarty，A.C.Hansen：AeroDyn Theory Manual，
National Renewable Energy Laboratory，(2005)
4) 喜屋武拓麻ら：3 次元 CFD と BEM による風車解析の比較，第 28 回数値流体力学シンポジウム，(2014)
5) 松隈秀和，宇都宮智昭：マルチボディダイナミク
スによる洋上風力発電用浮体基礎の動揺解析，土木学会，(2010)6) 石原孟ら：浮体の弾性変形を考慮した動揺予測モ
デルの開発，第 30 回風力エネルギー利用シンポ
ジウム，(2008) ，221-224

図 4 計算条件設定画面と実行ウインドウ

7) 洋上風力発電施設検討委員会：洋上風力発電設備に関する技術基準の統一的解説，(2020)

4 おわりに

8) 山口ら：風車制御が発電出力と支持物に作用する荷重に与える影響に関する研究，第 25 回風工学

本稿では，脱炭素化社会の実現に向けて近年注目

シンポジウム，(2018)が高まってきている浮体式洋上風力発電における数値シミュレーションの適用事例についてまとめた．

9) 眞鍋尚：遅延関数を用いた浮体の動揺シミュレー

また当該分野に関連する MHRT の取り組みとして，ション，富士総研技報，第 8 巻第 1 号(2003) ，116-

MARIS の紹介を行った．

131

MHRT では継続的に MARIS の改良を行っており，今後はより浮体式洋上風車の解析に適用するための改良を行い，当該分野の技術開発に貢献していきた
いと考えている．

8