数理解析研究所講究録
第 1989 巻 2016 年 8-158船舶の作る波の CFD 解析
横浜国立大学大学院工学研究院

$*$

日野孝則

Takanori Hino
Faculty of Engineering
Yokohama National University\S 1. はじめに水面上を航行する船舶の流体力学的性能を評価するためには自由表面影響を考慮する必要がある。 船体に働く流体抵抗の中で造波抵抗の占める割合は船の種類や船速にもよるが、 概ね 10% から 40% 程度である。 また、 実海域を航行する船の安全性においても海洋波の影響は極めて重要である。
したがって、 船舶流体力学の分野では、 古くから自由表面流れを扱ってきた。 しかし、後述するように自由表面条件は罪線形であり、 これが理論解析を困難にしてきた。 直進する船体に関しては、 微小波高を仮定して自由表面条件を線形化してポテンシャル理論をベースにした造波抵抗理論が展開された。 物体表面条件についても船の幅が微小と仮定して線形化近似を適用した、Michell 理論 [1] を初めとして多くの手法が開発された。 後に電子計算機の出現に伴って、 物体形状をパネルで近似し、 線形自由表面条件を満足するように吹き出し分布を求める Neumann-Kelvin 問題 [2] などが定式化された。
これらの造波抵抗理論により、 進行する船の後方に生じる八の字に広がる波紋 (Kelvinパターンと呼ばれる) が計算され、 船の作る波の特性が解明されてきた。 一方で造波抵抗の算定精度は実験結果との間に乖離があり、実用的な抵抗評価手法としては限界があった。
近年、 電子計算機の急速な発達によって、 計算流体力学 (CFD, Computational FuidDynamics) 手法が開発された。 これはナビエストークス方程式の数値解法であり、 計算領域を格子で分割して離散化し、 支配方程式と境界条件を満たすように格子点における〒 240-8501 横浜市保土ケ谷区常盤台 79-5

9圧力や速度の値を数値的に計算する手法である。 流体力学における実験、 理論に続く第三の解析アプローチとして期待されており、 流体物理や流体工学の種々の分野に適用されている。

船舶流体力学分野においても CFD の適用が進んでれており、 その中で自由表面の扱いも研究されている。 数値モデルとしては、 非線形自由表画条件を精度良く扱うことが重要となる。
ここでは、 船の作る波の CFD 解析における自由表薗の扱いを概観するとともに、 いくつかの計算例を紹介する。

\S 2. CFD 解析における自由表面モデル
船舶流体力学における CFD の支配方程式は 3 次元非圧縮ナビエストークス方程式であり、 空気と水の二相流れとして定式化すると、 次のように表すことができる。
$u_{i}$

を速度、$x_{i}$を空間座標として、 連続の条件は、(1)$\frac{\partial u_{j}}{\partial x_{j}}=0$であり、 運動方程式は$\frac{\partial(p_{l}u_{i})}{\partial t}+\frac{\partial}{\partial x_{J}\prime}[\rho_{i}u_{i}u_{j}+p\delta_{ij}-\tau_{ij}]-F_{i}=0, i=1, 2, 3$

ここで、

$t$

は時聞、

$p$

は圧力で、 粘性応力$\tau_{ij}$は(3)$\tau_{ij}=2\mu_{l}\mathcal{S}_{ij}^{*}\equiv\mu_{l}\{(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}})-\frac{2}{3}\frac{\partial u_{k}}{\partial x_{k}}\delta_{ij}\}$である。 また、 外力として、と粘性係数$\mu$$x_{3}$方向の重力$g$(2)のみを考えて、$F_{l}$’

lg i3、とし、 密度$\delta$$=\rho$$\rho$

は、 水と空気それぞれの値$(\rho_{i}, \mu_{l})=\{\begin{array}{ll}(p_{w}, \mu_{w}) 水中 (\rho_{a}, \mu_{a}) 空中\end{array}$(4)を切り替える。

自虫表面境界条件は、 運動学的条件、力学的条件の 2 つから構成される。 運動学的条件
は、 $h(x_{1}, x_{2}, x_{3}\rangle t)=0$ を自由表面形状とすれば、$\frac{\partial h}{\partial t}+u_{i}\frac{\partial h}{\partial x_{i}}=0$(5)10である。 これは、 自由表面上の流体粒子は自由表面上に留まり続ける、 と解釈される。
一方の力学的条件は、 以下のように定式化される。 応カテンソル $T_{ij}=-p\delta_{ij}+\tau_{ij}$ を、空中で$T_{a}$、水中で$T_{w}$として、 自由表画 $h(x_{1}, x_{2}, x_{3};t)=0$ の法線ベクトルを図 1 のように n(水面から水中へ向かうものを$n_{w\backslash }$空中へ向かうものを na)、接線ベクトルを$t$と定義すれば，法線方向の応カバランスは$n\cdot T_{w}n_{w}+n\cdot T_{a}n_{a}=\sigma(\nabla\cdot n)$である。 ここで、$\sigma$(6)は表面張力である。 $n=n_{w}=-n_{a}$ なので、 (6) 式は、$n\cdot T_{w}n-n\cdot T_{a}n=\sigma(\nabla\cdot n)$(7)となる。 また、 接線方向の応カバランスは以下の式で表される。
$t\cdot T_{w}n_{w}+t\cdot T_{a}n_{a}=0$(8)$n_{a}=-n$$n_{w}=n$$T_{w}$

図 1 自由表面の法線ベクトルと応力

CFD による船舶まわりの自由表面流れの解析においては、水と空気の二相を対象とする場合と水の部分のみを扱う単相を対象とする場合がある。 二相流アプローチでは、 流体を密度が異なる二つの流体が混合された一つの流体として解析することが多く、 その場合には運動学的条件により自由表面形状を決定し、 力学的条件は密度が変化するとして連続的に扱うことで自動的に満たされる。

単相アプローチでは、 水面における力学的条件を設定する必要がある。 空気の流れを無視すると、大気圧 (一定) を$p_{a}$として、 空気の応力$T_{a}$

は、$T_{a}=p_{a}\delta_{ij}$(9)であり、 水の応力も、 レイノルズ数 $Rn$ が大きいとして粘性項を無視すれば、$T_{w}=p\delta_{ij}$(10)$p=p_{a}-\sigma(\nabla\cdot n)$(11)となり、 結局法線方向の応カバランスは11で近似される。 表面張力を無視すれば(12)$p=p_{a}$であり、 水面上の圧力は大気圧に等しいことになる。

同様に、 接線方向の応カバランスは$t\cdot T_{w}n=0$(13)$\frac{\partial u_{i}}{\partial n}=0$(14)となり、 通常以下のように近似される。

これらの力学的境界条件は、 水面位置で満足する必要があるが、 水面位置はもう一つの境界条件である運動学的条件で規定されるので、 非線形境界条件となる。 この非線形性が自由表面を理論的に扱うことを困難にしており、 CFD などの数値的アプローチが期待される所以でもある。
自由表面を扱う数値モデルとしては、 大きく分けて、Euler 法と Lagrange 法がある。Euler 法
は水面形状を Euer 的に扱い、 計算格子上で水面位置が規定される。 水面位置と格子位置の関係によって、界面適合法と界面捕獲法に分けられる。
界面適合法は図 2(a) に示すように、 水面形状に沿うように計算格子を適合させる。 したがって、非定常流れでは計算格子が時間依存となるので、 移動格子を用いる必要がある。 また、 船舶のように水面を貫通する物体がある場合は格子は物体形状にも同時に適合する必要があり、 格子生成/格子変形は複雑な手続きとなる。 水面境界と計算格子が一致するので、 境界条件の計算精度が高いという利点がある。 反面、 前述の格子生成の手間に加え、 自由表面の大きな変形に対しては格子生成が著しく困難になるという短所がある。

界面捕獲法では、 計算格子は空間に固定され、 水面形状は何らかのマーカー関数を用いて格子上に定義される。 マーカー関数としては波高関数、 VOF(Volume of Fluid) 関数 (各セルの中の水の体積率) [3]、レベルセット関数 (水面からの符号付き距離) [4] などが用いられる。 運動学的条件を満たすために、 それぞれの関数に応じた方程式が用いられる。

(b)(a)
図 2 界面適合法 (a) と界面捕獲法 (b)Lagrange 法の代表的なアプローチが粒子法であり、 ナビエストークス方程式を Lagrange 的に解くために、 空間に粒子を分布させその動きを追跡する。 粒子分布から水面形状を求めることができる。 最大の特徴は、 計算格子を用いない点であり、 原理的に水面の大変形にも対応が可能であ12る。 一方で、 計算格子がないことから物体表面での流体力の積分における精度維持には注意が必要である。 粒子法では、 SPH(Smooth Particle Hydrodynamics) 法 [5] や MPS(Moving ParticleSemi-Implicit) 法 [6] などが代表的な手法である。

\S 4. CFD による計算例
CFD 手法による船の造る波の計算例をいくつか紹介する。
最初の例は、 直進するコンテナ運搬船のまわりの波形の計算例 [7] である。 速度は $230m$ の長さの実船で 24 ノット相当、 速度 $U$ と船長 $L$ および重力加速度$g$

による Froude 数 $(U/\sqrt{9^{L}})$ で0.26 である。 計算は長さが約 7. $3m$ の模型船を想定しており、 動粘性係数を$(UL/\nu)$$\nu$とした Reynolds 数は $1.4\cross 10^{7}$ である。 計算格子の点数は約 500 万点である。 計算の自由表面モデルには界面捕獲法を用いている。 図 3 に模型実験で計測された波紋図の計算結果との比較を示す。 船体から八の字状に広がる波の Kelvin パターンが再現されており、船首、船尾および水線面の変曲点から生じる複雑な波の成分も実験とよく対応している。

図 3 コンテナ船まわりの波紋の比較 (上 :実験、 下:計算)
次の例は、 肥大船モデルのまわりの波の計算例 [8] である。 船型は数式表示された肥大船モデルであり、 Froude 数を $0.3843$、Reynolds 数を $3.8\cross 10^{6}$ としている。 これは、 速度 1. $697m/s$ で進行する長さ 2. $0m$ の模型船まわりの流れに相当する。 図 4 は、 実験時の波形を計算値と比較したものである。 肥大船を高い Froude 数で走らせているので、 実験で見られる船首波や船体側方の波は13大変形し波崩れ現象を伴っている。界面捕獲法による自由表面モデルの適用により、 計算でも造波特性をある程度までシミュレートしている。

図 4 肥大船モデルまわりの波形 (左:実験、 右:計算))近年の高速船では、 船首の水面下に船首バルブと呼ばれる突起を装備して、 造波の低減を図ることが一般化している。 船首バルブの効果は、 主船体の船首で発生する波に対して、 逆位相の波を船首バルブで発生させて、 2 つの波の干渉効果により波高を低減させる、 と説明されている。 ここでは、 CFD によって、 船首バルブの効果を解析した。 対象はコンテナ船型で、 船速などのパラメータは上記と同じである。 図 5 は船首バルブの有無による計算波形の差を比較したものである。 船首バルブが無い場合 (上) では、 船首の波高が高く、 その後ろの波の山谷も船首バルブ付き (下) よりも大きいことが分かる。 船型形状の変更による流れ場や流体力の相違を解析することは、 CFD による船型設計において SBD(Simulation Based Design. シミュレーションをベースにした設計) を可能にするための重要要素である。

図 5 船首バルブの効果のシミュレーション (上:バルブなし、 下:バルブつき))14
最後の例は操縦運動時の船体まわりの流れの計算である。 船舶の安全性を考える際には、 その操縦特性が重要であり、 船舶の操縦性に関しては安全基準の遵守が求められている。 操縦性評価に必要となる操縦運動時に船体に働く流体力の評価にも CFD の適用が期待されている。 図 6 は操縦運動の基本モードである、 斜航時 (左) および旋回時 (右) の貨物船まわりの流れの計算値である。 斜航時には流れは船首左舷から入って船尾右舷へ向かうので、 波高も左舷側の方が大きい。 一方、旋回時の流れば、 船首左舷から入った流れは船底で方向を変えて船尾左舷から出て行くので、 波高分
布も複雑になる。

斜航旋回

図 6 斜航時 (左) および旋回時 (右) の船体まわりの流れ\S 5. おわりに船舶まわりの流れに対する計算流体力学 (CFD) 手法の適用における自由表面モデルの扱いを概観するとともに、 いくつかの計算例を紹介した。

船舶流体力学において自由表面は極めて重要な要素であり、 CFD による自由表面流れ解析の適用も多岐にわたっている。 海洋波のシミュレーションや波浪中で運動する船体まわりの非定常造波など、 本稿で触れなかった適用例については、 当該分野の関係文献などを参考にされたい。

参考文献[1] Michell, J. H., 1898. “The Wave Resistance of a Ship”. Phil. Mag., 45(5), pp. 106123.
[2] Brard, R., 1972, “The representation of a given ship form by singularity distriburions
when the boundary condition on the free surface is linearized”. J. Ship
, 16, pp.
${\rm Res}.$79-92.

15[3] Hirt, C.W. and Nichols, B.D., 1981, ‘Volume of Fluid (VOF) Method for the Dynamics
of Free Boundaries J. Comput. Phys., Vol.39, pp.201-225.
[4] Sussman, M., Smereka, P. and Osher, S., 1994, “ A Level Set Approach for ComputingSolutions to Incompressible Two-Phase Flow J. Comput. Phys., Vol.114, pp.146-159.
[5] Paredes, R., Imas, L., 2011, “Application of SPH in Fluid Structure Interaction Problems Involving &eesurface Hydrodynamics 11th Int. Conf. Fast Sea bansp., Honolulu,HI, 2011.
[6] Khayyer, A., Gotoh,
2011, (Enhancement of stability and accuracy of the moving
particle semiimplicit method J. Comput. Phys., Vol. 230, pp. 3093-3118.
[7] Hino, T. et al., 2010 “Flow Simulations Using Navier-Stokes Solver SURF Proc.
Gothenburg 2010 Workshop.
[8] Hino, $T$ , 2004,“Numerical Simulations of Breaking Waves around an Advancing Shipby an Unstructured NS Solver” , Proc. 25th Symp. Naval Hydrodynamics.
$H_{\rangle}$