第 27 回数値流体力学シンポジウム
A08-3

推進器影響を考慮した波浪中船体周り流れのシミュレーション
Numerical simulation of flows around a hull in regular waves with propeller effects
大橋 訓英, 海上技術安全研究所, 三鷹市新川 6-38-1, E-mail：k-ohashi@nmri.go.jp
日野 孝則, 横浜国立大学, 横浜市保土ヶ谷区常盤台 79-5, E-mail：hino@ynu.ac.jp
坂本 信晶, 海上技術安全研究所, 三鷹市新川 6-38-1, E-mail：sakamoto@nmri.go.jp

Kunihide Ohashi, National Maritime Research Institute, 6-38-1, Shinkawa, Mitaka, Tokyo
Takanori Hino, Yokohama National University, 79-5, Tokiwadai, Hodogaya, Yokohama, Kanagawa
Nobuaki Sakamoto, National Maritime Research Institute, 6-38-1, Shinkawa, Mitaka, Tokyo
Numerical simulation with regular head waves including propeller effects is carried out. Propeller effects arerepresented by body forces which are obtained by the simplified propeller model which is based on the infinite
blade theory. Ship heave and pitch motion equations are solved and ship motions are taken into account by agrid deformation. A practical tanker hull form is selected. Firstly, computed motion amplitudes and resistanceincrease coefficients are compared with measured data and resulted within a reasonable agreement. Then, flows
on the propeller plane with and without propeller effects are compared. Wake flows are changed with orbital
motions of the regular wave at a short wave length. On the other hand, in a long wave length, wake flows areeffected by ship motions and waves. Time averaged axial velocity which is also averaged in the propeller circle isincreased comparing with the steady condition without waves.

緒言
近 年, 数 値 流 体 力 学 (CFD) の 進 展 は め ざ ま し く,
Reynolds Averaged Navier-Stokes(RANS) によるシミュ
レーションは, 船舶性能設計の現場でも使用されている.
また, 船舶における地球温暖化ガスの削減スキームが国際条約として採択されたことにより, 実際の海域に近い状態における性能把握に向け, 波浪中での運動等, より複雑な流体現象を計算するニーズも高まりつつある.
著者らも, URANS による波浪中でのシミュレーション手法と数理モデル開発により, 船体運動等について計測結果との比較を通し, 妥当な結果を得てきた. (1)(2)(3)
ここでは, 波浪中での自由航走シミュレーションに向けて, 推進器影響を体積力で表すモデルに基づく, 推進器影響を考慮した波浪中シミュレーションにより, 推進器有無での流場や運動の変化について考察を行う.とができる. 自由表面モデルとして, 界面捕獲型の一つである一相の Level Set 法を適用する.
対流項 e, f , g, 拡散項 ev , f v , gv , 体積力項 H は下式の通りである.

βuβvβw 2
u + p vu 
 wu 
e= , f =  v 2 + p  , g =  wv  , (3)
uv

uwvw
w2 + p 0000




 
 τxx  v  τxy  v  τzx  fx  ev = τ , f = τ , g = τ , H = f  xy 
 yy 
 yz 
 y
τzx
τyz
τzz
fz

1.

2.

数値計算法

ここで β は擬似圧縮性パラメータである.


2.1 ナビエストークスソルバー SURF
本 計 算 で 使 用 し た CFD コ ー ド は 海 上 技 術 安 全 研究所で開発されている SURF(“Solution algorithm forUnstructured RaNS with FVM”)(4) (5) である. 支配方程式は 3 次元非圧縮レイノルズ平均ナビエストークス方程式であり, 擬似圧縮性を仮定し, 圧力と流場のカップリングを行う. また, 実時間と擬似時間を導入し, 非定常計算を行う.

τij = R1

上式で, 各変数は密度 ρ0 , 代表速度 U0 , 代表長さ L0 で無次元化され, (x, y, z) 方向の流速は (u, v, w) で表す. 実時間は t, 擬似時間は τ である. 自由表面流れでは, 圧力p を次式に変更する.

z
Fn2− u0i u0j であり, R はレイノルズ数 U0 L0 /ν, ν は動粘性係数, −u0i u0j はレイノルズ応力項である.
計算法については文献(4) (5) に詳しいので, 概略について述べる.
空間離散化は, 非構造格子を用い, セル中心に変数を配置した有限体積法によって評価する. セル形状は四面体, 三角柱, 四角錐, 六面体の組合せが可能である. 対流項は Flux-Difference-Splitting に基づく 2 次風上差分(MUSCL type), 拡散項は中心差分を用いる.
実時間微分は 2 次精度, 擬似時間微分は 1 次精度後退差分で近似し, 擬似時間については局所時間刻みによる収束加速を行う. ナビエストークス方程式を離散化及び線形化して得られる連立方程式はガウスザイデル法と多重格子法で解く.

∂q∗ ∂q ∂(e − ev ) ∂(f − f v ) ∂(g − gv )++++−H = 0∂t∂τ∂x∂y∂z(1)
TTq = [ p u v w ] , q∗ = [ 0 u v w ]p = p∗ +∂uj∂ui
∂xj + ∂xi

2.2 推進器影響のモデル化
推進器影響は, 文献(6) のポテンシャル理論に基づく, 翼数を無限大とした簡易プロペラモデルによる体積力で表す.

(2)pΓ(r, θ)Vθ 12 CD N c(r) 1 + (h(r)/r)2 Vox Voθ
fx (r, θ) =−r∆x
2πr∆x(4)√p∗ は場の圧力, Fn はフルード数 U0 / gL0 , z は高さ方向位置である. この変更により重力項を圧力に含めるこ1

Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A08-3

p
12 2
Γ(r, θ)Vx
2 CD N c(r) 1 + (h(r)/r) Voθft (r, θ) =+r∆x
2πr∆x(5)fy = ft cosθ, fz = ft sinθ(6)r, θ はモデルにおけるプロペラ面内の円筒座標系であり, Γ(r, θ) はプロペラ面を分割した要素上の流入速度に基づく方程式から得られる (r, θ) 位置での循環強さである.
Vx , Vθ はプロペラの誘導速度を含む流入速度の各方向成分, Vox , Voθ は流入速度を周方向に平均した成分である.
CD はプロペラ翼断面の粘性抗力の抵抗係数であり, Nは翼数, c(r) は半径方向各位置でのコード長さ, h(r) は半径方向各位置でのプロペラ後方の自由渦のピッチを表す.
流場と体積力のカップリングは次の手順による. まず,計算格子からプロペラ面上の流速 (u, v, w) を補間し, 体積力を計算する. 次に計算格子とプロペラ面が交差する領域を探索し, セルの交差面の面積と式 (4), (6) の体積力をかけあわせ, 推進器影響として考慮する. また, 擬似時間の繰返しステップの中でプロペラ流入速度を計算する
ことにより体積力も更新できる.

Fig. 1: Computational grids near hull
較できるが, サージの影響はほとんど見られない. 計算結果はいずれも概ね実験結果と一致していることが分かる.
Fig.4 に, 波長 λ/L = 0.6 でプロペラ影響なしの場合の船体固定座標でのプロペラ面位置伴流分布とベクトル図を示す. 入射波と船体の出会い周期を 4 分割した各時刻としている. また Fig.5 にプロペラ円内の平均流速とその時刻平均, さらにプロペラ位置 (x/L = 0.9825) での波高の無次元値と船体運動の無次元値の時刻歴を示す. 波
長が短い場合, 船体運動が小さくなるため, プロペラ面位置の伴流分布はほぼ入射波とともに変化し, 船尾縦渦も入射波に伴い伸縮する. プロペラ円内の平均流速も同様で, 時間平均値は入射波なしの場合より上昇する.

2.3 規則波モデルと運動モデル
規則波については, 文献(3) による吸収造波方法により計算領域内に入射させる. 吸収造波領域を計算格子外周から一定範囲に設定し, 波高に相当するレベルセット関数及び水粒子の速度と支配方程式を解くことで得られた値を, 計算領域外周からの距離に応じた sin の 2 次関数で混合させる.
運動モデルは文献(2) による手法であり, 重心周りの運動方程式を 1 次精度で離散化する. 擬似時間の繰返しス
テップの中で船体運動と流場をカップリングし, 運動の初期値は Adams-Bashforth 法, 修正値は Adams-Moulton法により与える. 船体運動は移動格子法に基づく船体近傍の格子変形により考慮する.

1.2

1.2
Meas. Surge Fixed
Meas. Surge FreeComp. Surge Fixed1

計算条件
計算対象は, 文献(3) と同条件下における, 向波規則波中のタンカー船型 (KVLCC2) である. 模型船スケールで船長は L = 3.2m, レイノルズ数は 2.55 × 106 , フルード数は0.142 である. 計算領域は片舷の場合, −1.5 ≤ x/L ≤ 3.5,−2.0 ≤ y/L ≤ 0, −2.0 ≤ z/L ≤ 0.1 であり, 原点は船首と静止水面が重なる位置とした. Fig.1 に船体近傍の計算格子を示す. 計算格子は局所に細分化された 6 面体より構成され, 船体近傍, 水面近傍と船尾付近は解像度を上げるため, より小さな 6 面体としている. 船体表面に垂直方向の最小格子間隔は 9.0 × 10−6 とした. プロペラ影
響を考慮する場合, y 対称面で計算格子の鏡像をとることで両舷化した. 片舷の格子総数は約 410 万点, 両舷で約820 万点である. プロペラ直径は D/L = 0.0306 であり,プロペラ回転数は平水中の模型自航点での無次元回転数n = 58.3 を与えた. また, プロペラによる体積力は 1 回転に相当する時間ステップで更新している.
乱流モデルには非線形の代数モデルである EASM(7) を適用し, 入射波として, 向波で波高 h/L = 0.01875 と 3 種の波長 λ/L = 0.6, 1.1, 1.6 を与え, 計算結果考察の際の時刻ゼロは, 波の頂点が船首と重なる時刻としている.
船体運動は z 方向の上下運動 (ヒーブ) と y 軸周りの回転運動 (ピッチ) の 2 自由度であり, ヒーブは z 方向と同じく正負をとり, ピッチは重心回りで船尾が沈む回転を正とした. 無次元時間刻みは ∆t = 0.003 − 0.005 とした.

3.

0.8/kz/

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

Meas. Surge Fixed
Meas. Surge FreeComp. Surge Fixed

100.5

11.50200.5

11.5/L/LHeavePitch2Fig. 2: Comparison of heave and pitch motion amplitudes
10Meas. Fixed Surge
Meas. Free Surge
Comp. Fixed Surge

98Raw/ g 2(B2/L)

7654321000.5

11.52/LFig. 3: Comparison of resistance increase coefficients

計算結果
Fig.2 にプロペラ非作動時における運動の振幅の無次元値について, 実験値の解析結果(8) との比較を, Fig.3 に同じく規則波なしの状態からの抵抗増加の無次元値の比較を示す. 実験結果では前後方向の運動 (サージ) 有無も比4.

Fig.6 と Fig.7 に波長 λ/L = 1.1 におけるプロペラ影響なしでの結果を示す. λ/L = 1.1 は入射波と船体運動の同調点付近の条件であり, λ/L = 0.6 の場合と比較し, 船体運動は大きくなり, 入射波と船体運動の干渉により, プロ
2Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A08-3

0.8w wave w/o Propeller
Average w Wave w/o Propeller
w/o Wave and Propeller

0.7

Uave

0.6

0.5

0.4t/T e = 0t/T e = 1/4
t/Te=1/4

0.3

00.1

0.2

t/Te=2/4

0.3

t/Te=3/4

0.4

0.5

0.6

0.7

Time(sec.)

0.6

1
Heave
Pitch
Wave height

0.5

0.01

0.3

-0.5

-0.3

-1

-0.6

z//k00h/L

0.005
0-0.005t/T e = 1/2t/T e = 3/4t/Te=1/4-1.5

Fig. 4: Axial flow contours and cross flow vectors(λ/L = 0.6, ∆U = 0.1)

00.1

0.2

0.3

t/Te=1/2

t/Te=3/4

0.4

0.6

0.5

0.7

-0.9

-0.01

TIme(sec.)

Fig. 5: Time history of averaged flow in propeller circle,motions and wave height (λ/L = 0.6)

ペラ面位置の伴流分布も大きく変動する. 時刻 t/T e = 0
では, 船体は上昇し, 船首が沈み込むため船体に引きずられる形で境界層は薄くなる. 時刻 t/T e = 1/2 では船体は沈下し, 船尾が沈み込むため, 境界層は伸長し, 船尾縦渦は明瞭でなくなる. プロペラ円内の平均流速にも影響が現れ, 入射波による変動に加え, 文献(9) で指摘されているように, ピッチ運動の影響が大きいように見られる. 時間平均値も入射波なしの場合より, 大幅に上昇する.
次に, Fig.8 と Fig.9 に波長 λ/L = 0.6 でプロペラ影響を考慮した場合の結果を示す. プロペラ面流速には, 入射波による加減速が見られ, ほぼ, プロペラ影響なしの場合(Fig.4) にプロペラの加速が重ねられた分布となることが分かる. プロペラ面内の平均流速とプロペラ推力は, プロペラ影響なしの場合と比較し, 入射波と若干の位相差をもち変動する. 短波長での, プロペラ作動時におけるプロペラ面内の平均流速の時間平均値は, 入射波なしの場合とほぼ同等となることから, 例えば自航要素等について,入射波なしの定常状態結果を用いて, 時間平均による諸量を解析することも可能なことが想定される.
最後に, Fig.10 と Fig.11 に波長 λ/L = 1.1 でプロペラ影響を考慮した場合の結果を示す.
プロペラ影響なしの場合 (Fig.6,Fig.7) と比較し船体運動と入射波の位相差が変わり, プロペラ面位置の伴流分布についても運動による影響が変化することにより, プロペラ影響なしの場合より時間差が見られる. 文献(10) の結果と比較し, 船体運動はプロペラ影響有無で差が大きいが, 文献(10) には時刻歴の結果がなく, 船型差も含めて,検討の余地がある. プロペラ推力とプロペラ円内平均流速についても入射波との位相差をもって変動し, 時間平均値は, プロペラ影響なしの場合より小さくなるものの,入射波なしの値と比較し, 高くなる傾向は残る. そのため, 長波長において, 入射波なしの定常状態結果を用いて時間平均による諸量を解析する場合, 入射波と船体の干渉を考慮していないことに留意する必要があると考えられる.の流場や運動の変化について考察した結果, 下記のことが言える.

• 推進器なしでの, 向波中の船体運動の振幅の無次元値と入射波なしの状態からの抵抗増加量の無次元値を実験結果と比較し, 概ね一致する結果が得られた.
• 推進器有無での船尾伴流を可視化した結果, 短波長
では, 伴流分布は入射波により変化し, プロペラ円内平均流速と推力は入射波と若干の位相差をもち変動する. プロペラ円内平均流速の一周期の時間平均値は入射波なしの場合とほぼ同等になる.
• 長波長では, 船体運動が大きくなり, 船体運動と入射波が干渉することで境界層は伸縮し, プロペラ円内平均流速と推力は大きく変動する。プロペラ円内平均流速の一周期の時間平均値は入射波なしの場合より高くなる.
参考文献(1) Larsson, L., Stern, F., Visonneau, M., “A Workshop on Numerical Ship Hydrodynamics Gothenburg 2010”(2) Sakamoto, N., Ohashi, K., “Unsteady wake analysis for Series 60 (Cb=0.6) pitching and heaving in
regular head seas by URANS,” 第 26 回 数値流力シンポジウム, (2012)(3) Ohashi, K., Sakamoto, N., Hino, T., “Numerical
Simulation of Flows around KVLCC2 Hull Formwith Ship Motions in Regular Waves,” Proc. ofECCOMAS MARINE 2013
(4) Hino, T., “Navier-Stokes Computations of Ship
Flows on Unstructured Grids,” Proc. of the 22nd
Symp. on Naval Hydro., (1998)(5) Hino, T., “An Interface Capturing Method forFree Surface Flow Computations on UnstructuredGrids,” J. of the Soc. Naval Archit. Japan , Vol.186(2003), pp. 173-183.

結言
波浪中での自由航走シミュレーションに向けて, 推進器影響を体積力で計算する手法を開発し, 推進器有無で

5.

3Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A08-3

0.8w wave w/o Propeller
Average w Wave w/o Propeller
w/o Wave and Propeller

0.7

Uave

0.6

t/T e = 0t/T e = 1/4

0.5

0.4

t/Te=1/4

0.3

0

0.1

0.2

0.3

0.4

t/Te=2/4

0.5

0.6

0.7

t/Te=3/4

0.8

0.9

1

1.1

Time(sec.)
HeavePitch
Wave height1

0.6

0.5

0.3

00-0.5

-0.3

-1

-0.6

0.01

Fig. 6: Axial flow contours and cross flow vectors(λ/L = 1.1, ∆U = 0.1)

0h/Lt/T e = 3/4z/t/T e = 1/2/k0.005-0.005t/Te=1/4-1.5

0t/Te=1/2t/Te=3/4

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

11.1-0.9-0.01

TIme(sec.)(6) Ohashi, K., Hirata, N., Hino, T., “A comparativestudy of body force models representing effects ofcontrarotating propellers,” Proc. 105th Meeting of
the West-Japan Society of Naval Architects, (2003),
pp. 55-64,

Fig. 7: Time history of averaged flow in propeller circle,
motions and wave height (λ/L = 1.1)

(7) Sadat-Hosseini, H., Wu, P.C., Carrica, P.M., Kim,H., Toda, Y., Stern, F., “CFD verification and validation of added resistance and motions of KVLCC2with fixed and free surge in short and long headwaves,” Ocean Engineering, Vol.59 (2013), pp.240273.
(8) , “回転楕円体周り流場の DES,” 日本流体力学会年会 2011
(9) 塚田, 日夏, 長谷川, “波浪中における非定常船尾伴流の計測,” 関西造船協会誌, 第 228 号, (1997) , pp.1520(10) Nakamura, S., Naito, S., “Propulsive Performanceof a Container Ship in Waves,” Naval Arch. andOcean Engineering, 15, (1977) , pp.24-48
t/T e = 0t/T e = 1/4

t/T e = 1/2

t/T e = 3/4

Fig. 8: Axial flow contours and cross flow vectors(λ/L = 0.6, ∆U = 0.1)

4

Copyright c 2013 by JSFM

第 27 回数値流体力学シンポジウム
A08-3

1.3

1.40E-03
Total Velocity
Average Total Velocity
w/o Wave
Propeller Thrust

1.30E-03

1.2

1.20E-03
1.10E-03

1.1

9.00E-04
1ThrustUave

1.00E-03

8.00E-04
7.00E-04

0.9

6.00E-04
5.00E-04
t/Te=1/4

0.8

00.1

0.2

t/Te=2/4

0.3

0.4

t/Te=3/4

0.50.60.7

4.00E-04

Time(sec.)

0.6

1Heave
Pitch
Wave height

0.5

0.01

0.3

-0.3

z/0

1.3

1.30E-03

1.2-0.005

1.20E-03

-0.6

-1

0

0.1

0.2

t/Te=1/2

0.3

0.4

1.10E-03t/Te=3/4

0.5

0.6

0.7

-0.9

1.1

-0.01

1.00E-03

Uave

t/Te=1/4-1.5

1.40E-03
Total Velocity
Average Total Velocity
w/o Wave
Propeller Thrust

TIme(sec.)

9.00E-04

Fig. 9: Time history of averaged flow in propeller circle,
motions and wave height (λ/L = 0.6)

1Thrust-0.5/k0h/L

0.005
08.00E-04
7.00E-04

0.9

6.00E-04
5.00E-04
t/Te=1/4

0.80t/Te=2/4

t/Te=3/4

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.91

1.1

4.00E-04

Time(sec.)

0.6

1
HeavePitch
Wave height

0.5

0.01

0.3

-0.5

-0.3

-1-0.6

z//k

00h/L

0.005
0-0.005t/Te=1/4-1.5

0t/Te=1/2t/Te=3/4

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

1.1-0.9-0.01

TIme(sec.)

Fig. 11: Time history of averaged flow in propeller
circle, motions and wave height (λ/L = 1.1)
t/T e = 0t/T e = 1/4

t/T e = 1/2

t/T e = 3/4

Fig. 10: Axial flow contours and cross flow vectors(λ/L = 1.1, ∆U = 0.1)

5Copyright c 2013 by JSFM