WEC-Sim (Ver.2) Theory
訳: 今井康貴（imaiy@cc.saga-u.ac.jp）

理論
本節は WEC-SIM コードの基本的な理論を概説する。

はじめに
WEC のモデル化は、入射波、装置運動、PTO 機構および係留の間の相互作用を含む。WEC-SIM
は、動力性能と設計を最適するためラディエーションおよびディフラクション法[2][3]を使用する。
一般に、ラディエーションとディフラクション法は、周波数領域における境界要素法(BEM)ソルバを用いて流体力を求め、線形重ね合わせを使用して、時間領域における系の動特性を予測する。
座標系
図は 3 次元の浮体式ポイントアブソーバに波が入射する様子を示す。さらに WEC-SIM における座標系および 6 つの運動自由度（DOF）を示す。WEC-SIM 座標系では、角度ゼロの場合、X 軸に波が伝達する方向と定義する。Z 軸は垂直の上向きとする。また、Y 軸方向は右手法則で定義する。ベクトルと行列では、方向 1,2,3 がそれぞれサージ(x), スウェイ(y)、ヒーブ(z)と一致する。
方向 4,5,6 がそれぞれロール(Rx), ピッチ(Ry), ヨー(Rz)に一致する。

境界要素法
流体力を求める一般的なアプローチは、線形波仮定の使用である。つまり波浪場は入射波成分、ラディエーション成分、ディフラクション成分の和と仮定する。WEC-SIM では、これらの成分
は、周波数領域のポテンシャル流れ BEM ソルバ（例えば WAMIT[4]、AQWA-FER[5]およびNemoh[6])から得られる線形係数を用いてモデル化される。BEM による解はラプラス方程式を速度ポテンシャルを用いて解くことにより得られる。それは、流れが非粘性、非圧縮、渦なしと仮定する。周波数領域 BEM の理論の詳細は[4]にある。WAMIT、AQWA、Nemoh の BEM ソルバ
による解は BEMIO によって*.h5 データ構造に変換され、WEC-SIM は無次元流体力係数をその*.h5 データからインポートする。あるいは、*.h5 データ構造は、ユーザーによって手動で作成される。
次に、WEC-SIM は流体力係数を以下のようにスケーリングする。ρは水の密度、ωは周波数（rad/sec）、g は重力である:
Fx (ω ) =Fx (ω )ρg,A (ω ) =A (ω )ρ, B (ω ) =B (ω )ρω, K=K

ρg(1)

時間領域計算系の動的応答は、WEC の運動方程式を解くことにより求められる[3][7]。浮体の運動方程式は、その重心に関して次式で与られる:1mXɺɺ = Fext + Frad + FPTO + Fv + FB + Fm

(2)

ここで、 Xɺɺ は装置の（並進および回転）加速度ベクトル、 m は質量ベクトルである。右辺は以下の通り
Fext ：波強制力ベクトル、
Frad ：ラディエーション力ベクトル
FPTO ：PTO により発生する力のベクトル
Fv ：減衰力ベクトル
FME ：モリソン力ベクトル
FB ：復原力ベクトル
Fm ：係留力ベクトル
Fext と Frad は周波数領域の BEM ソルバで計算される。ラディエーション項は、付加質量と波減衰から構成され、それぞれ浮体の加速度および速度に関係する。波強制項は、物体が存在しない場合に入射波によって生成されるフルード=クリロフ力成分と、浮体存在に起因するディフラクション成分から構成される。
WEC-SIM は規則波および不規則波シミュレーションに使用可能であるが、 Fext と Frad は正弦的な定常応答シナリオとランダム海象仮定で求められていることに注意する。正弦定常応答は、単純な WEC に規則波が入射する設計にしばしば使用される。しかし、ランダム海象シミュレーシ
ョンあるいは係の流体メモリ効果が不可欠なシミュレーションでは、浮体の流体メモリ遅延力を表わすために重畳積分法が推奨される。
正弦定常応答シナリオ
このアプローチは、系の応答が正弦定常的と仮定し、規則波シミュレーションにのみ有効である。
ラディエーション項は、与えられた波周波数における付加質量および波ラディエーション減衰項をから計算される。
Frad = − A (ω ) Xɺɺ − B (ω ) Xɺ(3)

ここで A (ω ) と B (ω ) はそれぞれ付加質量およびラディエーション減衰行列である。 ω は波周波数 rad/sec であり、 Xɺ は浮体の速度ベクトルである。
自由表面形状は、与えられた波高、波周波数、水深における線形波の仮定に基づいて与えられる。
規則波の波励起力は次式で与えられる。
 Hi ωt Fext = ℜ  R f
FX (ω ) e ( ) 2

ℜ は式の実部を表す。 R f はランプ関数、 H は波高、 FX(4)は波強制力ベクトルであり、力の振幅と位相で表される。
重畳積分系の流体メモリ効果を取り入れるため、カミンズ方程式[8]に基づく重畳積分計算が使用される。
ラディエーション項は次式で計算される。
2

tFrad = − A∞ Xɺɺ − ∫ K ( t − τ ) Xɺ (τ ) dτ

(5)0

ここで A∞ は周波数無限大の負荷質量行列、 K はラディエーションインパルス応答関数である。
規則波の場合、波強制力は最終節で記述される式で計算される。自由表面は、不規則波の場合、多数の規則波成分の線形重合で構築される。それは、波スペクトルで特徴づけられる。波スペク
トルは波周波数ごとの波エネルギー分布を記述し、有義波高およびピーク波周期で特徴づけられる。不規則波の波強制力は、すべての周波数にわたる積分の実部であり、次式で表される。
 ∞i (ω + φ )Fext = ℜ  R f ∫ FX (ωγ ) e γ
2S (ωγ ) dωγ 
 0


(6)

ここで S は入射波スペクトルであり、 φ は任意の位相角である。
Krの状態空間表現

ラディエーション重畳積分を状態空間(SS)形式で表わすことが望ましい[9]。これは劇的に計算速度を増加させ、かつ従来の制御方法(線形状態空間モデルに依存する)の適用を可能にする[10]。
1 セットの偏微分方程式を解いて K r を求め、1 セットの常微分方程式で線形の状態空間を構築して K r を近似する。一般に、それは以下の近似で表される。
Xɺ r ( t ) = Ar X r ( t ) + Br ζɺ ( t ) ,X r ( 0) = 0(7)t∫ K ( t − τ ) dτ ≈ C X ( t ) + D ζɺ ( t )rrrr(8)−∞Ar Br Cr Dr は時間非依存、入力、出力であり、 ζは系への入力である。

状態空間行列からの K r 計算
単一入力の状態空間モデルのインパルス応答は次式で表される。
xɺ = Ar x + Br u(9)y = Cr x(10)

これは強制力の無い応答（ u = 0 ）と同じである。
初期条件 Br をもつ場合は次式になる。
xɺ = Ar x, x ( 0 ) = Br(11)y = Cr x(12)非ゼロの D 行列をもつ連続系のインパルス応答は、t=0 で無限である。したがって、t=0 で CrBrは連続値をもつ。線形時不変系（LTI 系）の一般解は次式で与えられる。
tx ( t ) = e Ar t x ( 0 ) + ∫ e r (A t −τ )Br u (τ ) dτ(13)0

ここで e A は行列指数関数である。 K r は次式で求まる。
rK r ( t ) = Cr e Ar t Br3(14)

ラプラス変換と伝達関数ラプラス変換は一般的な積分変換である。時間領域の関数 f ( t ) を複素数 s の関数 F(s)に変換する線形演算子であり、次式で計算される。
∞

F ( s ) = ∫ f ( t ) e − st dt(15)

0f ( t ) の微分のラプラス変換は次式になる。

f ( t ) − ste dt0 dt∞sF ( s ) = ∫(16)ラプラス変換はいくつかの便利な関係がある。最初の関係は初期値定理および最終値定理である。
f ( 0+ ) = lim sF ( s )(17)f ( ∞ ) = lim sF ( s )

(18)

dmy
d m −1 yd nu
d nu+ a1 m −1 + ⋯ + am y = b0 n + b1 n + ⋯ + bnu
mdt
dt
dt
dt

(19)

s →∞s →0

次式で表される線形入出力系を考える。

u が入力、 y が出力を示す。これをラプラス変換すると、微分方程式は 2 つの多項式になる。
A ( s ) = s m + a1s m −1 + ⋯ + am −1s + am

(20)

B ( s ) = b0 s + b1s(21)

nn −1+ ⋯ + bn −1s + bnA ( s ) は系の特性方程式である。伝達関数 G ( s ) は次式で定義される。

G(s) =Y (s)U (s)=s m + a1s m −1 + ⋯ + am −1s + amb0 s n + b1s n −1 + ⋯ + bn −1s + bn

(22)

状態、入力、出力、フィードスルー行列が既知の場合、系の伝達関数は次式で計算できる。
G ( s ) = C ( sI − A) B + D
−1

(23)

系の周波数応答は周波数範囲において、s に jσ を代入することにより得られる。一般的に、G ( jσ )の大きさおよび位相はボード線図で表現される。
Realization Theory （実現理論）
流体のラディエーション係数の状態空間実現は時間領域(TD)で追求できる。
これは、系のインパルス応答関数サンプルの行列アセンブリから最小オーダーおよび離散時間状態行列(Ad、Bd、Cd、Dd)を見つけることである。この問題は離散時間系のために扱いが連続時間より簡単である。離散時間系のインパルス応答関数である理由は、系のマルコフ・パラメタから示される。
Kɶ r ( tk ) = Cd Adk Bd(24)tk = k ∆t 、 k = 0, 1, 2,… ∆t はサンプル周期である。t=0 で無限の値に帰着し、系の因果律を維持するため、上記の方程式はフィードスルー行列を含まない。
実現を得る最も一般的なアルゴリズムは、Kung [11]によって提案されたように、インパルス応答4関数をハンケル行列上で特異値分解(SVD)を行なうことである。系のオーダーおよび状態空間パラメタは、SVD の固有値数およびの係数から決定される。SVD は次式で実行される。
 K r ( 2 ) K r ( 3)K ( 3) K r ( 4 )H = r ⋮⋮0
 K r ( n )⋯ K r ( n )⋯0 ⋱⋮ ⋯0 

H = U ∑V *(25)(26)

ここで H はハンケル行列、 ∑ はハンケル特異値（降順）の対角行列である。
ハンケル特異値の検討は、特異状態はわずかであり、大きな精度損害なしモデルを縮小できることを明らかにした [10][12]。ここでは状態空間パラメタ計算における SVD 法のさらなる詳細は議論されない。[10][12]を参照すること。

波スペクトル
規則波を生成する能力は、特定条件におけるモデルの応答を観察する機会を提供する。
一定の波高および周期の海象は水槽試験外ではほとんど発生しない。
通常の海象はランダム波の時系列によってより正確に表わされる。ランダム時系列は異なる振幅および周期の様々な形状の重ね合わせをモデル化する。
波の重ね合せはスペクトルで特徴づけられる。
統計分析によって、スペクトルは、有義波高、ピーク周期、風速、フェッチ長のようなパラメタで特徴づけられる。
海洋産業で使用される一般的なスペクトルは次節で議論される。
WEC-SIM で使用される海洋スペクトルの一般形態は、次式で与えられる。
S ( f ) = Af −5 exp  − Bf −4 (27)

ここで f は波周波数 [Hz]、exp は指数関数を示す。

Pierson-Moskowitz Spectrum
最も単純なスペクトルの 1 つは[13]で提案された。風が大きなエリア一帯に定常的に長期間吹いた後、波が風と平衡すると仮定した。これは完全に荒れた海の概念である。「長期間」はおよそ10,000 周期であり、「大きなエリア」はおよそ 5,000 波長である。
スペクトルは次式で計算される。
S( f ) = 5  fp  
α PM g 2 −5f exp  −   4
 4  f  
( 2π )

A=4α PM g 25, B = f p444
( 2π )(28)(29)

α PM = 0.0081 、 g は重力加速度、 f p はスペクトルのピーク周波数である。しかし、このスペクトル表現は、ユーザーが有義波高を定義できない。WEC-SIM では、パワーマトリクスを生成する
ため、αPM 係数は目的海象の有義波高に合致するように計算される。αPM は次式で計算され5る。
α PM =∞

H m2 016 ∫ S * ( f ) df

(30)

0

 5  f p 4 f −5 exp  −   
 4  f  

(31)

mk = ∫ f k S ( f ) df , k = 0,1, 2,…

(32)

g

S* ( f ) =2( 2π )4スペクトルモーメントは次式で定義される。
∞

0m0 は自由表面の変位であり、スペクトルの有義波高 H m 0 に使用される。

(33)

H m 0 = 4 m0

Bretschneider Spectrum ブレトシュナイダー・スペクトル有義波高およびピーク波周波数の 2 パラメタ・スペクトルである。与えられた有義波高について、成長あるいは減衰も含め、一連の条件をカバーするためにピーク周波数が変更できる。
一般に、パラメタは風速(最も重要)、風向、フェッチ（吹送距離）、暴風雨前線の位置に依存する。
The spectrum is given as S(f)=Hm024(1.057fp)4f?5exp[?54(fpf)4]
スペクトルは次式で与えられる。
S( f ) = 5  f p 4 
4H m021.057 f p ) f −5 exp  −   (4 4  f  A=

4H m0251.057 f p ) ≈ H m 0 2 f p4(416B = (1.057 f p ) ≈45 4fp
4(34)
(35)
(36)

H m 0 は有義波高であり、一般に大きい方から 1/3 の平均値として定義される。

JONSWAP (Joint North Sea Wave Project) スペクトル※これが WEC-SIM のデフォルトのスペクトルである。
スペクトルは Hasselmann らにより提案された[14]。オリジナルは次式である。
 5  f p 4 −5−    y Γfexp4
 4  f  ( 2π )

(37)

4
  f
 −1  f 
Γ = exp  −  p 2σ      (38)S( f ) =α jg2

0.07, f ≤ f p
0.09, f > f p

σ =

A=

α jg2

( 2π )
64(39)
(40)

B=

5 4fp
4(41)

α j は無次元変数であり、風速とフェッチ長の関数である。

ほとんどの計測海象のスペクトル形状を平均表記するため、経験的な近似を行った。α j を希望の有義波高に近似するため、次式を計算する必要がある。
αj =H m02(42)∞16 ∫ S * ( f ) df
0

S* ( f ) = 5  f p 4 −5−    yΓ
fexp
4
 4  f  ( 2π )g2(43)

Spectrum purposed at ITTC （ITTC スペクトル）
JONSWAP スペクトルの別形状が第 17 回 International Towing Tank Conference (ITTC)で提案された。それは次式で提案された。
S* ( f ) = 5  f p 4 
 5  f p 4 310−5Γ2 4 −5−    yΓfexp−y≈Hffexp m0p444
 4  f  
 4  f  ( 2π )( 2π ) ( 0.834Tp )155

H m024  f −1  f 
Γ = exp  −  p 2σ      (45)0.07, f ≤ f p
0.09, f > f p(46)310

H m 0 2 f p4(47)5 4fp4(48)σ =
A=(44)( 2π )
B=4

以下に、 α j を近似した JONSWAP スペクトルと ITTC スペクトルを示す。2 つスペクトルは非常によく合致する。

Ramp Function（ランプ関数）
ランプ関数 (Ramp function: Rf) はシミュレーションの初期段階で強い移流を避けるために必要であり、波強制力の計算に用いられる。ランプ関数は次式で与えられる。
1 πt  1 + cos  π +   ,2tr  Rf =  1

t はシミュレーション時間、 tr はランプ時間である。

7t<1trt≥1tr(49)

Power Take-off (PTO) Forces （PTO 力）
Linear Spring-damper PTO 線形バネ・減衰 PTO
PTO 機構は線形のバネ・減衰系で表わされる。反力は次式で与えられる。
FPTO = − K PTO X rel − CPTO Xɺ rel(50)K PTO は PTO の剛性、 CPTO は PTO の減衰、 X rel と Xɺ rel はそれぞれ 2 体の相対変位と相対速度である。
PTO で消費されたパワーは次式で与えられる:
PPTO = − FPTO Xɺ rel = ( K PTO X rel Xɺ rel + CPTO Xɺ rel 2 )(51)

しかし、2 体の相対運動および速度はπ/2 によって位相が異なるため、変位と速度の積の時間平均はゼロになる。そのため、吸収されたパワーは次式になる。
PPTO = C PTO Xɺ rel 2(52)

油圧 PTOPTO 機構は油圧系としてモデル化され、反力は次式で与えられる。
FPTO = ∆ piston Apiston(53)∆ piston は油圧ピストンの差圧であり、 Apiston はピストン断面積である。PTO の absorbed power は次式になる。
PPTO = − FPTO Xɺ rel(54)Xɺ rel は 2 体の相対速度を表す。

Mechanical PTO （機械 PTO）
PTO 機構をダイレクトドライブのリニア発電機としてモデル化する[15]。反力は次式で与えられる。
 π FPTO = λ i τ pm  fd sq(55)τ pm は磁極ピッチ（磁石の距離。磁極中心間の距離）、λ fd は d 軸巻線のフラックスリンクである。

線の磁束鎖交である。また、 isq ステータ q 方向の電流である。
PTO の吸収パワーは次式で表される。
PPTO = − FPTO Xɺ rel(56)Mooring Forces（係留力）
係留荷重は線形の準静的な係留ばねを使用して表わされるか、Moordyn [1] (lumped-mass 法による運動モデル)から計算された、係留力の使用により表される。
線形の準静的な係留剛性が使用される場合、係留荷重は次式で表される。
Fm = − K m X − Cm Xɺ

K m は係留剛性、 Cm は減衰である。 X(57)および Xɺ はそれぞれ物体の変位および速度を表す。

MoorDyn で係留する場合、係留系の係船索は、それぞれ節点で接続される等長の線分要素に分割される。線分要素の質量は、重力と浮力、流体力荷重、海底との接触反力とともに節点節点にお8ける集中荷重とする。流力抗力および付加質量はモリソン式で計算される。係船索の軸剛性は各線分要素に線形剛性を適用することによりモデル化する。ここでは張力のみ考慮する。
集中質量の離散化によって発生する共振を減衰させるため、減衰項も各要素に適用される。曲げおよびねじり剛性を無視する。節点が海底以下に変位する場合、海底接触は鉛直剛性と減衰力によって表わした。

Additional Added-mass & Damping Forces 補足付加質量&減衰力抗力の寄与を追加し、実験値と比較して数値モデルを校正するため、付加質量と減衰力をモリソン式を用いて一次／二次の項として追加する。
Linear & Quadratic Damping Forces 一次・二次の減衰力
系の減衰影響は、運動方程式に追加する一次・二次の減衰項で指定される。
1Fv = −Cld Xɺ rel − Cd ρ AD Xɺ Xɺ
2(58)

Cld は一次減衰係数、 Cd は（二次）粘性減衰係数、 ρ は流体の密度、 AD は断面積である。

一般に、WEC 動力学では粘性影響を考慮する必要がある。特に線形モデルを使用する場合、この影響を無視すると、系の発電を過大に評価する。
粘性減衰をモデル化する一般的な方法は、運動方程式に二次の減衰項を加えることである。
装置用の粘性抵抗係数は注意深く選択しなければならない[2][3]。
しかし、それは、装置形状、スケール、船体と周囲流れの相対速度に依存する。
Re および KC 数が小さい場合、抗力係数は大きくなる。
抗力係数の経験的なデータは様々な文献および規準に記載されている。
しかし、利用可能なデータは既存の単純な幾何形状に制限されることがある。
現実的なポイントアブソーバ形状では、流体力は水槽試験や規定運動の数値シミュレーションにより評価されなければならない。
Morison Elements モリソン要素
モリソン式は、振動流中で柱に似た構造物に作用する流体力を仮定する。ポテンシャル流れによる影響と粘性影響を考える。柱は、直径(D)が波長 λw に比べて小さいことを意味する。一般に D/λ<0.1~0.2 である。この条件が満足されない場合、ディフラクション効果を考慮する必要がある。
細長形状を仮定すると、結果として生じる力は修正モリソン式で接近できる。[17] 物体上の各要素の公式は次式で与えられる。
1FME = ρ∀vɺ + ρ∀Ca ( vɺ − Xɺ ) + Cd ρ AD ( v − Xɺ ) v − Xɺ
2v は流体粒子速度、 ∀ は排水容積である。

現在、WEC-Sim はモリソン要素を計算する際、浮力効果を考えていないことに注意する。
9(59)Nonlinear Hydrodynamic Forces 非線形流体力
線形モデルは船体運動と波の振幅が波長に比べて小さいことを仮定する。瞬間の水面上昇と船体位置によって生じる非線形流体力を説明するため、弱非線形アプローチが適用される。
BEM による線形の流体力係数の計算ではなく、各時間ステップにおいて没水面に沿って各パネルの静圧と動圧を積分することにより、非線形の浮力およびフルード・クリロフ力成分を求める。
速度場と圧力場の決定に線形波理論が使用されるため、平均水位より上にある没水パネルの値は非現実的に大きくなる。
これを修正するため、Wheeler ストレッチ法が使用される。[18] これは流速と圧力を計算する際、水深に等しい高さを持つことを(瞬間の水面に基づいた)水柱に強いる。
z* =D(D + z)
D +η

−D

(60)

D は平均水深、 η は瞬時の水面 z である。

非線形 WEC-Sim は、スラミングや砕波等の高度に非線形な流体イベントを対象としていない。
References 参考文献
（略）

10