Journal of Marine Science and Technology (2019) 24:884–901

https://doi.org/10.1007/s00773-018-0594-7

ORIGINAL ARTICLE

Development of a structured overset Navier–Stokes solver

with a moving grid and full multigrid method
Kunihide Ohashi1 · Takanori Hino2 · Hiroshi Kobayashi1 · Naoyuki Onodera3 · Nobuaki Sakamoto1

Received: 15 February 2018 / Accepted: 29 August 2018 / Published online: 8 September 2018
© JASNAOE 2018

Abstract
An unsteady Reynolds averaged Navier–Stokes solver with a structured overset grid method has been developed. Velocity–
pressure coupling is achieved using an artificial compressibility approach, spatial discretization is based on a finite-volume
method, and convective fluxes are evaluated by a third-order upwind scheme based on flux-difference splitting. Body motions
are considered using the grid deformation technique and grid velocities in the convective term. Viscous fluxes are evaluated
via second-order centered differencing. The full multigrid (FMG) method is applied to obtain fast convergence. The cell flag
on a coarse grid level is determined using the cell flag on a fine grid level. In the coarse grid-level calculations at the FMG
stage, the data are interpolated until the finest grid level is achieved at an overset update interval. Then, the data are updated
based on the overset relations at the finest grid level and then transferred to a coarser grid level. For free surface treatments,
a single-phase localized level-set method is employed. The computations for flows around a hull form, including an unsteady
simulation with regular waves, are demonstrated.

Keywords Overset grid · Full multigrid method · Moving grid · Localized level-set method

1 Introduction are mainly compared with experimental data, and the

errors between the calculated results and the measured
In the field of hydrodynamics, numerical simulations are values are summarized. The 2010 workshop was held in
widely utilized to develop ship hull shapes and energy- Gothenburg of Sweden [3], and new challenging test cases
saving devices, especially at the design stages. In addition, were set. Several cases in head sea waves with freedom
numerical simulations play an important role when deter- of ship motions were tested, and the local flow quantities
mining experimental conditions, including the measurement (which include the turbulent kinetic energy and Reyn-
of flow positions and quantities (e.g., [1]). olds stresses) were compared. Additionally, the flow field
The progress of numerical methods in ship hydrody- around a ship hull with the roll decay motions were exam-
namics can be observed in series of workshops in this ined. The assessment from the 2010 workshop indicated
field. In the 2005 workshop [2], several hull forms under that the standard deviation of the resistance coefficient
steady towing conditions were selected as the standard from the experimental towing condition data was 2.1%,
test cases, whereas oblique towing and wave diffraction which is slightly reduced from that of the 2005 workshop.
without ship motion conditions were selected as the spe- The two equation turbulent models, including the alge-
cial test cases. The resistance coefficients and wave fields braic Reynolds stress model, have become the standard
models. Reynolds stress components are consistent with
* Kunihide Ohashi experimental data when using a reasonably fine grid. The
k‑[email protected] resistance coefficient and motions in waves have a cer-
tain degree of accuracy compared with the measured data,
1
National Maritime Research Institute, 6‑38‑1 Shinkawa, however, the errors are larger than the values associated
Mitaka, Tokyo, Japan
with calm water conditions. The latest workshop for ship
2
Yokohama National University, 79‑5 Tokiwadai, Hodogaya, hydrodynamics was held in 2015 [4], and further complex
Yokohama, Kanagawa, Japan
cases were presented. A ship with motions under regular
3
Japan Atomic Energy Agency, 148‑4, 178‑4, Wakashiba, waves with a certain wave angle was selected, and a bulk
Kashiwa, Chiba, Japan

13
Vol:.(1234567890)
Journal of Marine Science and Technology (2019) 24:884–901 885

carrier hull form with an energy-saving device under a et al. verified and validated the numerical method at the
self-propulsion condition was tested, which means that conditions with and without the energy-saving device [17].
the marine propeller was working behind the main hull. An overset grid method with a structured grid has sev-
In a more complex case, a twin screw vessel, which has a eral advantages. Once the baseline grids are generated
complex geometry in a free running condition with regular and assembled, the modification of shapes of appendages,
waves, was also set. The assessment of the most recent such as energy-saving devices, is completed by replacing
workshop is now being performed. the appendage grid and updating the overset information.
The unstructured grid and overset grid methods are Many variations of configurations can be treated in a shorter
the two major methods in the ship hydrodynamics fields. time. Additionally, a structured grid method can be easily
An unstructured grid that is fit to the body shapes is the extended to higher order discretization methods.
first choice for calculating complex geometries. Recently, In this work, a new structured solver based on the artifi-
a second-order discretization in space and time with the cial compressibility approach is developed. The solver can
finite volume method and the volume of fluid method for treat overset grids with relation information, which is gener-
free surface treatments has become the general method ated by the in-house system [18], and it is also suitable for
applied for ship hydrodynamics (e.g., [5, 6]). The open the design process of energy-saving devices. The discretiza-
source code for the Computer Aided Engineering (CAE), tion of convective and viscous fluxes with a moving grid
openFOAM solver is also utilized to determine the steady technique is described, and the linearization process for the
towing condition. In addition, a wave decomposition system equation is introduced. The full multigrid method
model based on the phase field equation is applied to (FMG) [19, 20] is applied to obtain fast convergence. A
openFOAM solver [7–9]. Regarding the unstructured grid special treatment must be used for the FMG to apply the
method in the field of ship hydrodynamics, a local sub- overset grid method. The cell flag to be solved or unsolved
division of hexahedral cells is mainly employed because is inherited from finer grid levels to coarser grid levels on
of its affinity for the features of flows around a hull. The the multigrids. In the coarse-grid level calculations of the
local subdivision hexahedral method can easily adapt FMG stage, the data transfers are performed until reaching
the local grid refinement to resolve the boundary layers the finest grid level at an overset update interval, and then
or the air–water interface. The disadvantage associated data interpolation is performed to reflect the overset rela-
with using the unstructured grid method is that once the tions. For a case in which the body surface is overlapped
body shape changes, the computational grid must be re- between the overset grids, an overlapped region is eliminated
generated from the beginning. Modern ships often have based on the priority of the overset relation to calculate the
energy-saving devices in the aft part of the body surface, correct forces and moments. Free surface treatment which
and the design process for energy-saving devices is usually is based on the single-phase localized level-set method, the
based on trial and error. Thus, the re-generating process physical models that include turbulence and propulsion are
may not be efficient. also described. To achieve simulations that include ship
Recently, an overset grid method was re-examined and motions in waves, the wave model is implemented and body
applied to the ship hydrodynamics field. Carrica et al. motions are obtained by solving the motion equations. The
applied the single-phase level-set method [10] with the body motions are accounted for by a moving grid technique
dynamic overset grid method, including 6DOF ship motions with a grid-deforming methodology. To avoid updates of
to the flows around ships in waves [11]. The Suggar code the overset information with the grid-deforming process,
[12] is used for the overset grid method and coupled with the the grid-overlapped regions around the body move with
structured solver. Sadat-Hosseini et al. extended the method the body motion without altering the relative locations. For
to the flows around the tanker hull form and conducted a regions away from the body, the amount of deformation
detailed analysis [13]. Broglia et al. applied the in-house gradually decreases with distance from the body. Thus, the
structured solver and the overset grid method to maneuver- grids that are far from the body are fixed. This method is
ing condition [14]. Propeller models which represent pro- adapted to avoid the computational load using the dynamic
peller effect with body forces are examined on simulations overset grid method.
of turning motions of a ship. The complex configuration First, the present method with the FMG is applied to steady
which includes the twin propeller shafts and the brackets, flows around a hull with several appendages to examine the
the bilge keels and the rudder is composed using the grids influence of an overset update interval for the solution inter-
of the appendages and the overset grid method. Di Mascio polation and transfer of the multigrid method. Additionally,
et al. applied the single-phase level-set method to the flows the verification based on the generalized Richardson extrapo-
around the standard test case ship (Series 60) and a surface lation method is performed. Second, the propeller model is
combatant [15]. Broglia et al. computed the flow around a applied to the flows around a hull with and without the energy-
high speed craft [16] using the same method [15]. Korkmaz saving devices, and the applicability of the propeller model

13
886 Journal of Marine Science and Technology (2019) 24:884–901

( )
is examined in Sect. 3.1.2. Next, computations with regular 𝜕uj
is defined as 𝜏ij = − u�i u�j , R is the Reynolds
1 𝜕ui
+
waves on several wavelengths including the ship motions, are R 𝜕xj 𝜕xi

demonstrated in Sect. 3.2. Finally, we conclude the present number, 𝜈 is the kinematic viscosity coefficient, and −u�i u�j is
study.
the Reynolds stress component, which is determined by the
turbulence model.
2 Numerical method The integrational form of Eq. 1 is expressed as Eq. 5
using the finite volume method and Gauss integral theorem.
2.1 Governing equations

𝜕t ∫ ∫ ∫ 𝜕𝜏 ∫ ∫ ∫ ∫∫∫
𝜕 𝜕
The governing equation is the incompressible three-dimen- qdV + q∗ dV − HdV
sional Reynolds-averaged Navier–Stokes equation using the Vi,j,k Vi,j,k Vi,j,k

artificial compressibility approach to couple pressure and ( )

∮
+ (e − ev )nx + (f − f v )ny + (g − gv )nz dS = 0.
velocities. The dual time stepping method is applied to an
unsteady computation. 𝜕Vi,j,k
(4)
v
𝜕q 𝜕q∗ v
𝜕(e − e ) 𝜕(f − f ) 𝜕(g − g ) v
+ + + + − H = 0. Finally, we derive the discretized form of the following
𝜕t 𝜕𝜏 𝜕x 𝜕y 𝜕z
[ ]T [ ]T equation for a structured grid with a cell-centered layout:
q = 0 u v w , q ∗= p u v w .
(1) 𝜕Vi,j,k qi,j,k 𝜕Vi,j,k q∗i,j,k
+ − Vi,j,k Hi,j,k
Equation 1 is non-dimensionalized by the reference density 𝜕t 𝜕𝜏
[ ]i+1∕2 [ ]j+1∕2 [ ]k+1∕2
𝜌0 , reference velocity U0 , and reference length L0 . (u, v, w) + E − Ev i−1∕2 + F − Fv j−1∕2 + G − Gv k−1∕2 = 0,
represent the velocity in the (x, y, z) directions, respec- (5)
tively. Temporal time is expressed by t, and artificial time is
where ±1∕2 means directions of each cell face, and E , F , G
expressed by 𝜏 . The pressure p is modified in the following
are the convective fluxes, Ev, Fv, Gv are the viscous fluxes.
equation for free surface flows:

p = p∗ +
z
. 2.2 Convective flux with the moving grid technique
Fn2 (2)
The convective flux E is expressed as
where p∗ is the pressure in a computational domain, Fn is
the Froude number, and z is the vertical coordinate from the ⎡ 𝛽U ⎤
static water plane ( z = 0). The gravity term is included with ⎢u(U − Ug ) + pSx⎥
E=⎢
v(U − Ug ) + pSy⎥
, (6)
the pressure in this modification. ⎢ ⎥
The convective terms e, f , g, viscous terms ev, f v, gv and ⎣w(U − Ug ) + pSz⎦
additional body force term H are defined as follows:
where U = uSx + vSy + wSz , Ug = ug Sx + vg Sy + wg Sz . The
⎡ 𝛽u ⎤ ⎡ 𝛽v ⎤ convective flux Ei+1∕2 is evaluated using the upwinding
⎢(u − ug )u + p⎥ ⎢ (v − vg )u ⎥
e =⎢ scheme based on the flux difference splitting scheme of Roe
(u − ug )v ⎥
, f =⎢
(v − vg )v + p⎥
,
⎢ ⎥ ⎢ ⎥ [21],
⎣ (u − ug )w ⎦ ⎣ (v − vg )w ⎦
1[ ]
⎡ 𝛽w ⎤ Ei+1∕2 = E(qL ) + E(qR ) − |A|(qL − qR ) i+1∕2 , (7)
⎢ (w − wg )u ⎥ 2
g=⎢
(w − wg )v ⎥
, (3)
⎢ ⎥ where qL and qR are the left and right interface values of q at
⎣(w − wg )w + p⎦ the i + 1∕2 face, respectively. E(qL ) and E(qR ) are evaluated
⎡0⎤ ⎡0⎤ ⎡0⎤ ⎡0⎤ by qL , qR and i + 1∕2 face area vector (Sx , Sy , Sz ) with grid
⎢𝜏 ⎥ ⎢𝜏 ⎥ ⎢𝜏 ⎥ ⎢f⎥ velocity Ug . |A| is obtained in the following way. First, A is
ev =⎢ xx⎥, f v = ⎢ xy⎥, gv = ⎢ zx⎥, H = ⎢ x⎥, defined as follows:
𝜏 𝜏 𝜏 f
⎢ xy⎥ ⎢ yy⎥ ⎢ yz⎥ ⎢ y⎥
⎣𝜏zx⎦ ⎣𝜏yz⎦ ⎣𝜏zz⎦ ⎣fz⎦
⎡0 𝛽Sx 𝛽Sy 𝛽Sz ⎤
𝜕E ⎢ Sx U − Ug + uSx uSy uSz ⎥
A= ∗ =⎢ ⎥.
where (ug , vg , wg ) are the grid velocities caused by the mov- 𝜕q ⎢ S y vSx U − Ug + vSy vS z ⎥
⎣ Sz wSx wSy U − Ug + wSz ⎦
ing grid, 𝛽 is a parameter of the artificial compressibility
approach and 𝛽 = 1.0 is given in the present computation. 𝜏ij (8)

13
Journal of Marine Science and Technology (2019) 24:884–901 887

Ug Ug Ug
The eigenvalues of A above are U − Ug, U − Ug, U − 2
+ c, where 𝛬 = diag(U − Ug , U − Ug , U − 2
+ c, U − 2
− c) ,
Ug
and U − 2 − c, where c is the pseudo-speed of sound, which |A| is also expressed as follows:
is defined as follows: |A| = R|𝛬|L, (15)
√ Ug
Ug 2
(9) where |𝛬| = diag(|U − Ug |, |U − Ug |, |U − + c|,
c = (U − ) + 𝛽(Sx2 + Sy2 + Sz2 ). 2
2 |U −
Ug
− c|).
2
The right and left eigenvectors of A are derived as follows: The values qL and qR on the i + 1∕2 face are evalu-
� � � � ated using the MUSCL[22] approach with the third-order
⎡0 0 U
𝛽 c + 2g
U
− 𝛽 c − 2g ⎤
⎢ � � � �⎥ upwinding scheme:
⎢ lx mx 𝛽Sx + u U + c − Ug
𝛽Sx + u U − c − 2 ⎥
Ug

R =⎢ � 2 � � � ⎥,
qLi+1∕2 =
2 5
qi+1 + qi −
1
q ,
⎢ l m 𝛽S + v U + c − Ug 𝛽Sy + v U − c − Ug ⎥ 6 6 6 i−1
⎢ y y y
� 2 � � 2 �⎥ (16)
⎢ l m 𝛽S + w U + c − Ug 𝛽Sz + w U − c − Ug ⎥ 2 5 1
⎣ z z ⎦ qRi+1∕2 = qi + qi+1 − q .
z 2 2 6 6 6 i+2
(10)
The fluxes of F and G are similarly obtained by replacing
1 the suffix i with j, k.
L=
2𝛽c{𝛽(Sx2 + Sy2 + s2z ) + U(U − Ug )}
2.3 Viscous flux
⎡ 2𝛽c{mx (vSz − wSy ) + my (wSx − uSz ) + mz (uSy − vSx )}
⎢ 2𝛽c{−l (vS − wS ) − l (wS − uS ) − l (uS − vS ))}
⎢ x z y y x �z z � y
Ug
x The viscous flux Ev is expressed as follows:
⎢ 2 2 2 2
U + 𝛽(Sx + Sy + Sz ) − c + 2 U
⎢ � � ⎡ 0 ⎤
⎢ U
−{U 2 + 𝛽(Sx2 + Sy2 + Sz2 ) − −c + 2g U} ⎢𝜏xx Sx + 𝜏xy Sy + 𝜏zx Sz⎥
⎣
Ev = ⎢
𝜏 S + 𝜏yy Sy + 𝜏yz Sz⎥
. (17)
2𝛽c{my (U − Ug )w + my 𝛽Sz − mz (U − Ug )v − mz 𝛽Sy } ⎢ xy x ⎥
2𝛽c{−ly (U − Ug )w −�ly 𝛽Sz +�lz (U − Ug )v + lz 𝛽Sy } ⎣𝜏zx Sx + 𝜏yz Sy + 𝜏zz Sz⎦
U
−𝛽 2g − c Sx
�U � The velocity gradient on a face is calculated using the diver-
𝛽 2g + c Sx
gence theorem for the control volume constructed around the
2𝛽c{mz (U − Ug )u + mz 𝛽Sx − mx (U − Ug )w − 𝛽mx Sz } face i + 1∕2 as shown in Fig. 1. The face centers i and i + 1
2𝛽c{−lz (U − Ug )u − �lz 𝛽Sx + �
lx (U − Ug )w + lx 𝛽Sz } are located at the cell center, thus the cell center values are
U
−𝛽 2g − c Sy utilized (red circles in Fig. 1). Other variables are obtained
�U � by averaging the values that surround the point.
𝛽 2g + c Sy
The viscous flux of Fv , and Gv is similarly obtained by
2𝛽c{mx (U − Ug )v + mx 𝛽Sy − my (U − Ug )u − my 𝛽Sx } ⎤ replacing the suffix i with j, k.
2𝛽c{−lx (U − Ug )v −�lx 𝛽Sy +�ly (U − Ug )u + ly 𝛽Sx ⎥
⎥
U
−𝛽 2g − c Sz ⎥, 2.4 Time discretization
�U � ⎥
𝛽 2g + c Sz ⎥
⎦ Time discretization is separated into temporal t and pseudo
(11) time 𝜏 . The temporal time step (n) is discretized using three-
level backward differencing, and the first-order backward
where (lx , ly , lz ) and (mx , my , mz ) are the vectors that satisfy Euler scheme is applied to the pseudo-time step (m).
the following relations:
ly mz − lz my = Sx , lz mx − lx mz = Sy , lx my − ly mx = Sz
(12)
and
lx Sx + ly Sy + lz Sz = 0, mx Sx + my Sy + mz Sz = 0. (13)

Using the eigenvectors R and L , A is expressed as follows:

A = R𝛬L, (14)
Fig. 1  Control volume for the velocity gradient at i + 1/2

13
888 Journal of Marine Science and Technology (2019) 24:884–901

(n+1) (n+1)
3Vi,j,k (n) (n)
qi,j,k − 4Vi,j,k (n−1) (n−1)
qi,j,k + Vi,j,k qi,j,k ⎛ 3V (n+1)(m) V (n+1)(m)
𝜕Vi qi ⎜ i,j,k +
i,j,k
≈ , ⎜ 2𝛥t 𝛥𝜏i,j,k
𝜕t 2𝛥t ⎝
(n+1) (18)
𝜕Vi q∗i Vi,j,k 𝛥q∗(m+1)
≈
i,j,k
, 𝜕E(m)
i+1∕2
𝜕E(m)
i−1∕2
𝜕E(m)
vi+1∕2
𝜕E(m)
vi−1∕2
𝜕𝜏 𝛥𝜏 + − − +
𝜕q∗(m)
i
𝜕q∗(m)
i
𝜕q∗(m)
i
𝜕q∗(m)
i

where 𝛥q∗(m+1) = q∗(m+1) − q∗(m) . 𝜕F(m)

j+1∕2
𝜕F(m)
j−1∕2
𝜕F(m)
vj+1∕2
𝜕F(m)
vj−1∕2
i,j,k i,j,k i,j,k + − − +
The pseudo-time step process is iterated at each temporal 𝜕q∗(m)
j
𝜕q∗(m)
j
𝜕q∗(m)
j
𝜕q∗(m)
j
time step to satisfy the divergence free condition, and the ⎞
𝜕G(m) 𝜕G(m) 𝜕G(m) 𝜕G(m)
local time stepping approach is adopted for the pseudo-time +
k+1∕2
−
k−1∕2
−
vk+1∕2
+
vk−1∕2
⎟𝛥q∗(n+1)(m+1)
∗(m) ⎟ i,j,k
step 𝛥𝜏 to obtain fast convergence, 𝜕q∗(m)
k
𝜕q∗(m)
k
𝜕q∗(m)
k
𝜕qk ⎠
� �i+1∕2 � �j+1∕2
Vi,j,k = E(m) − Ev (m) i−1∕2 + F(m) − Fv (m) j−1∕2
𝛥𝜏i,j,k = CFL ∑ Ug
. (19) � �k+1∕2
faces 0.5(�U − 2
� + c) (n+1)
+ G(m) − Gv (m) k−1∕2 + Vi,j,k Hi,j,k
⎛ 𝜕E(m) 𝜕E(m) ⎞
− ⎜ ∗(m) − ⎟𝛥q∗(n+1)(m+1)
i+1∕2 vi+1∕2
2.5 Determination of grid velocities
⎜ 𝜕q ∗(m) ⎟
𝜕qi+1 ⎠
i+1,j,k
⎝ i+1
The grid velocity Ug i+1∕2 used to satisfy the geometrical con- ⎛ 𝜕E(m) ⎞
𝜕E(m)
servation law is derived from the volume where the face + ⎜ ∗(m) −
i−1∕2 vi−1∕2
⎟𝛥q∗(n+1)(m+1)
i + 1∕2 sweeps. Ug i+1∕2 is evaluated using the same discre- ⎜ 𝜕q ∗(m) ⎟
𝜕qi−1 ⎠
i−1,j,k
⎝ i−1
tization as the temporal time step. The swept volume at face
⎛ 𝜕F(m) 𝜕F(m) ⎞
i + 1∕2 and time step n is expressed as 𝛥Vi+1∕2
n
, and the time
− ⎜ ∗(m) − ⎟𝛥q∗(n+1)(m+1)
j+1∕2 vj+1∕2
⎜ 𝜕q ∗(m) ⎟ i,j+1,k
step n + 1 is 𝛥Vi+1∕2
n+1
. Then, the grid velocity Ug i+1∕2 is ⎝ j+1 𝜕qj+1 ⎠

obtained as follows: ⎛ 𝜕F(m) 𝜕F(m) ⎞

+ ⎜ ∗(m) − ⎟𝛥q∗(n+1)(m+1)
j−1∕2 vj−1∕2
⎜ 𝜕q ∗(m) ⎟
𝜕qj−1 ⎠
i,j−1,k
(n+1)
3𝛥Vi+1∕2 (n)
− 𝛥Vi+1∕2 ⎝ j−1
Ug i+1∕2 = . (20) ⎛ 𝜕G(m) ⎞
2𝛥t 𝜕G(m)
−⎜ ⎟𝛥q∗(n+1)(m+1)
k+1∕2 vk+1∕2
−
⎜ 𝜕q∗(m) ∗(m) ⎟
𝜕qk+1 ⎠
i,j,k+1
The grid velocities in the j, k directions are similarly ⎝ k+1
obtained by replacing the suffix i with j, k. ⎛ 𝜕G(m) 𝜕G(m) ⎞
+⎜ ⎟𝛥q∗(n+1)(m+1)
k−1∕2 vk−1∕2
−
⎜ 𝜕q∗(m) ∗(m) ⎟ i,j,k−1
2.6 Linearization of equations ⎝ k−1
𝜕qk−1 ⎠

(3Vi(n+1)(m) q(n+1)(m) − 4Vi(n) q(n) + Vi(n−1) q(n−1) )

The following linearization is performed for the convective −
i,j,k i,j,k i,j,k
.
and viscous fluxes: 2𝛥t
(22)
𝜕E(m)
E(m+1) = E(m) + ∗(m) 𝛥q∗(m+1) ,
𝜕q The symmetric Gauss–Seidel (SGS) method is applied to
(21) solve Eq. 22.
𝜕E(m)
v
E(m+1)
v
= E(m)
v
+ 𝛥q∗(m+1) .
𝜕q∗(m)
2.7 Applying the overset grid method
The convective term in Eq. 21 is evaluated using the first- The weight values for the overset interpolation are deter-
order upwinding scheme. mined by an in-house system [18]. The details of the system
Finally, we have the linearized and implicit forms of Eq. 5 can be found in [18], and a summary is provided below:
as follows:
1. The priority of the computational grid is set.
2. The cells of a lower priority grid and those inside a body
are identified (called in-wall cells herein).

13
Journal of Marine Science and Technology (2019) 24:884–901 889

eliminated. The ratio where the face is completely covered

by the stern fin becomes zero, and the face that is partly cov-
ered by the stern fin has a specific value (right-hand side of
Fig. 3). The figures show that the present method produces
accurate results

2.8 Full multigrid method

Fig. 2  Schematic diagram for the elimination of lapped body surfaces
The multigrid method can obtain fast convergence by reducing
the residual, which has a large wavelength on a coarse grid.
3. Receptor cells for which the flow variables must be The agglomeration type is applied to the present method. A
interpolated from donor cells are defined. Two layers coarse grid is defined for two cells in each (i, j, k) direction,
of cells on a higher priority grid and facing the outer thus, one cell on the coarse grid is constructed by eight cells
boundary are set as receptor cells to satisfy the third- on the fine grid.
order discretization of the NS solver. Additionally, two The solution ql on the fine grid level (l) is transferred to the
cells that neighbor the in-wall cells, cells of a lower pri- coarse grid level (l + 1) as follows:
ority grid and cells inside the domain of a higher priority
grid are also set as the receptor cells. l
ql+1 = Tl+1 ql , (23)
4. The weight values for the overset interpolation are deter-
mined by solving the inverse problem based on the Fer- where Tl+1
l
is defined as follows:
guson spline interpolation.
(24)
l
Tl+1 ql = (𝛴ql Vl )∕Vl+1 .
Flow variables of the receptor cell are updated when the
The equation to be solved on the coarse grid level is
boundary condition is set. The forces and moments are inte-
expressed as follows:
grated on the higher priority grid to eliminate the lapped
region on body surfaces. First, the cell face of the lower 𝜕ql+1 Vl+1
= −Rl+1 (ql+1 ) − Pl+1 , (25)
priority grid is divided into small faces. Second, the small 𝜕𝜏
face is projected to the cell face of the higher priority grid
where Rl+1 is the residual on the coarse grid and Pl+1 is the
using the normal vector of the higher priority face. Then, the
forcing term, which is given as follows:
2D solid angle is calculated and the small face is determined
to be in or out of the higher priority face (Fig. 2). Once the
small face is in the higher priority face, the area ratio is set to
Pl+1 = Ql+1
l
Rl (ql ) − Rl+1 (ql+1 ), (26)
zero. Finally, the area ratio is integrated on the higher prior- where Ql+1 is the transfer operator for the residual and
ity face, and then used to integrate the forces and moments
l
defined by the sum of the eight cells in the fine grid,
on the lower priority face.
Figure 3 shows one result of the eliminated area ratio. Ql+1 Rl (ql ) = 𝛴(Rl ). (27)
The main hull, refined hull grid and stern fin grid are
l

lapped. First, the lapped region between the main hull and The solution on the coarse grid is expressed as q+l+1, then the
refined hull grid is eliminated. The area ratio of the main correction to the fine grid level is defined as follows:
hull becomes zero (left-hand side of Fig. 3). Then, the
lapped region between the refined hull and stern fin grid is q+l = ql + Il+1
l
(q+l+1 − ql+1 ), (28)

where Il+1
l
is an interpolation operator for a correction.
The computation starts from the coarsest grid on the full
multigrid method (FMG) and the solution is interpolated to
the sequentially finer grids after a certain number of steps
at each stage (Fig. 4). By constructing the better initial solu-
tion for the finest grid using the grid sequence, the FMG is
expected to achieve fast convergence [19, 20]. When this
FMG technique is combined to the overset grid method, the
overset interpolation relations have to be defined for each
grid level. Since this is sometimes very difficult, the pre-
Fig. 3  Eliminated area ratio on ship hull surfaces sent solver adopts a novel method for the solution exchange

13
890 Journal of Marine Science and Technology (2019) 24:884–901

The extrapolation approach is applied to avoid a singularity

in which the solid wall condition prevents free surface move-
ment and causes large deformation near the wall. The level-set
function 𝜙 near a solid wall region is corrected by solving the
extrapolation equation as follows:
𝜕𝜙 𝜕𝜙 𝜕𝜙 𝜕𝜙
− Cx − Cy − Cz = 0, (29)
𝜕𝜏 𝜕x 𝜕y 𝜕z
where (Cx , Cy , Cz ) are the vectors of the gradient of the dis-
Fig. 4  FMG sequence and the reflection method for the overset rela- tance from a solid wall.
tion
The dynamic free surface condition can be approximated
using the following two conditions in the single-phase method.
between overset girds during the computations on a coarser First, the velocity gradients normal to the free surface are zero.
grid stage of the FMG sequence (Fig. 4). In this approach, To satisfy the condition, the velocity components are extrapo-
the overset interpolation relations are defined in the finest lated to the air region where 𝜙 < 0 by solving the following
grid level only and the overset updates in a coarser grid equation:
stage are performed in the designated interval of steps (over-
𝜕q 𝜕q 𝜕q 𝜕q
set update interval). In the overset update process, first the + Cx + Cy + Cz = 0, (30)
solution on the current grid level is interpolated to the fin- 𝜕𝜏 𝜕x 𝜕y 𝜕z
est grid using the interpolation operator Il+1 l
defined above. where (Cx , Cy , Cz ) are the tangential vectors near the wall
The interpolated solutions at the finest grid are exchanged region and the vector of gradient 𝜙 away from the wall.
based on the overset relation on the finest level. Then the Second, the pressure on the free surface is given as
updated solutions are transferred back to the current gird atmospheric pressure:
level using the solution transfer operator Tl+1 l
. The influence ZC
of the overset update interval is also investigated in the pre- p= on the free surface, (31)
Fn2
sent research.
The cell flag to be solved or unsolved on a coarse grid where atmospheric pressure is assumed to be zero and ZC
level is determined using the cell flag on a fine grid level. An is the z-coordinate of the interface. As shown in Fig. 6, the
unsolved cell flag is included in the cell on a coarse grid, and closest point on the interface from the cell (i) is denoted as
the cell is flagged as an unsolved cell (Fig. 5). Furthermore, a A and the closest point from the cell (i + 1) is denoted as
cell flag on the coarser grid level is similarly determined. Once B. In addition, the intersection of the interface and the line
the receptor cell value on the finest grid is included in coarser connecting the cell centers is denoted as C. From the rela-
grid level, the correction to the finer grid from a coarse grid tion, ZC is given as follows:
is treated as zero.
|𝜙i |Zi+1 + |𝜙i+1 |Zi
ZC = . (32)
2.9 Free surface model |𝜙i | + |𝜙i+1 |

In the remaining air region, pressure is extrapolated using

The single-phase localized level-set method [23, 24] is applied Eq. 30.
to capture free surfaces around bodies except employing a
third-order upwinding scheme for the convective flux of the
level-set function and the extrapolation approach near the solid 2.10 Propulsion modeling
wall and to the air region.
The propeller model based on the potential theory [25] is
applied to achieve the propulsion condition. The propeller

Fig. 5  Flagging of cells with multigrid levels Fig. 6  Pressure condition on a free surface

13
Journal of Marine Science and Technology (2019) 24:884–901 891

effect is accounted for by the body forces that are calcu-

lated using the following equations: m(ȧ � − b� r� + c� q� ) = Fx�
√ m(ḃ � − c� p� + a� r� ) = Fy�
1
𝛤 (r, 𝜃)V𝜃 C Nc(r) 1 + (h(r)∕r)2 Vox Vo𝜃
2 D m(ċ � − a� q� + b� p� ) = Fz�
fx (r, 𝜃) = − , (36)
r 2𝜋r Ix ṗ � + (Iz − Iy )q� r� = Mx�
(33)
Iy q̇ � + (Ix − Iz )r� p� = My�
1
√
2
𝛤 (r, 𝜃)Vx C Nc(r) 1 + (h(r)∕r)2 Vo𝜃 Iz ṙ � + (Iy − Ix )p� q� = Mz� ,
ft (r, 𝜃) = + 2 D
, (34)
r 2𝜋r

fy = ft cos 𝜃, fz = ft sin 𝜃, (35) where m is the system mass, and Ix , Iy , and Iz are the prod-
ucts of the system mass and the moments of inertia in the
where r and 𝜃 are the cylindrical coordinates at the propeller
principal directions, respectively. (a� , b� , c� ) is the linear
plane. The propeller circle is divided into fan-shaped seg-
velocity of a motion, and (p� , q� , r� ) is the angular velocity
ments at (r, 𝜃). 𝛤 (r, 𝜃) is the vortex strength, and 𝛤 (r, 𝜃) is
of a motion on the body-fixed coordinate.
determined using the equation based on the boundary condi-
The coordinate transformation can be performed using
tion at the segment (r, 𝜃). Vx , and V𝜃 are the inflow veloci-
Euler angles [11]. Except for the angular velocities from the
ties to the segment, and Vox , and Vo𝜃 are the circumferential
earth coordinate to the body-fixed coordinate, the transfer
averaged velocities. CD is a drag coefficient that is given by
matrix is presented as follows:
the empirical formula, N is the blade number of a propeller,

⎡ cos 𝜓 cos 𝜃 sin 𝜓 cos 𝜃 − sin 𝜃 ⎤ (37)

Jb = ⎢ − sin 𝜓 cos 𝜙 + sin 𝜙 sin 𝜃 cos 𝜙 cos 𝜓 cos 𝜙 + sin 𝜙 sin 𝜃 sin 𝜓 sin 𝜙 cos 𝜃 ⎥.
⎢ ⎥
⎣ sin 𝜃 sin 𝜓 + cos 𝜙 sin 𝜃 cos 𝜓 − sin 𝜙 cos 𝜓 + cos 𝜙 sin 𝜃 sin 𝜓 cos 𝜃 cos 𝜙 ⎦

c(r) is a chord length at each radial direction, and h(r) is the For the angular velocities, the transfer matrix is defined as
pitch of a free stream vortex. follows:
The coupling between the flow field and the propeller
model is performed by the velocities Vx and V𝜃 . First, the ⎡1 0 − sin 𝜃 ⎤
velocity at the segment center is interpolated from the com- Jb = ⎢ 0 cos 𝜙 cos 𝜃 sin 𝜙 ⎥. (38)
⎢ ⎥
putational grid. Then, the 𝛤 (r, 𝜃) is determined, and the body ⎣ 0 − sin 𝜙 cos 𝜃 cos 𝜙 ⎦
forces are calculated. Finally, the cross-section between the
computational grid and the propeller segment is searched,
As an example, hydrodynamic forces (Fx , Fy , Fz ) on the
and the forces at the cell are derived by multiplying the body
earth-fixed coordinate are transformed to the body-fixed
forces and the cross-sectional area. The coordinate transfer
coordinate as follows:
from the body-fixed coordinate to the earth-fixed coordinate
or its inverse is performed using Euler angles as described in
the following section.
(39)

⎡ Fx� ⎤ ⎡ cos 𝜓 cos 𝜃 sin 𝜓 cos 𝜃 − sin 𝜃 ⎤⎡ Fx ⎤

⎢ Fy� ⎥ = ⎢ − sin 𝜓 cos 𝜙 + sin 𝜙 sin 𝜃 cos 𝜙 cos 𝜓 cos 𝜙 + sin 𝜙 sin 𝜃 sin 𝜓 sin 𝜙 cos 𝜃 ⎥⎢ Fy ⎥.
⎢ ⎥ ⎢ ⎥⎢ ⎥
⎣ Fz� ⎦ ⎣ sin 𝜃 sin 𝜓 + cos 𝜙 sin 𝜃 cos 𝜓 − sin 𝜙 cos 𝜓 + cos 𝜙 sin 𝜃 sin 𝜓 cos 𝜃 cos 𝜙 ⎦⎣ Fz ⎦

2.11 Arbitrary body motions Hydrodynamic moments (Mx , My , Mz ) on the earth-fixed

coordinate are similarly transformed.
Body motions are defined as prescribed or obtained by solving Except for the angular velocities from the body-fixed
the motion equation. The motion equation for the body-fixed coordinate to the earth-fixed coordinate , the transfer matrix
coordinate is as follows: is presented as follows:

13
892 Journal of Marine Science and Technology (2019) 24:884–901

⎡ cos 𝜃 cos 𝜓 sin 𝜙 sin 𝜃 cos 𝜓 − cos 𝜙 sin 𝜓 cos 𝜙 sin 𝜃 cos 𝜓 + sin 𝜙 sin 𝜓 ⎤
Je = ⎢ cos 𝜃 sin 𝜓 sin 𝜙 sin 𝜃 sin 𝜓 + cos 𝜙 cos 𝜓 cos 𝜙 sin 𝜃 sin 𝜓 − sin 𝜙 cos 𝜓 ⎥. (40)
⎢ ⎥
⎣ − sin 𝜃 sin 𝜙 cos 𝜃 cos 𝜙 cos 𝜃 ⎦

For the angular velocities, the transfer matrix is defined as three-level backward differencing, and the body accelera-
follows: tion is added to the pressure.

⎡ 1 sin 𝜙 tan 𝜃 cos 𝜙 tan 𝜃 ⎤ 2.12 Wave model

Je = ⎢ 0 cos 𝜙 − sin 𝜙 ⎥. (41)
⎢ ⎥
⎣ 0 sin 𝜙∕ cos 𝜃 cos 𝜙∕ cos 𝜃 ⎦ The incoming waves are generated in the region that is des-
ignated near the boundaries of a computational domain. The
The accelerations of motions are obtained by solving the flow variables are blended between the solution qsolved and
Eq. 36, and the velocity ḣ , which is representative of all analytical value qwave through the blending function 𝛼 , as
components, is derived from the following equation: follows:

ḣ n+1 = ḣ n + 𝛥tḧ n+1 . (42) q = (1 − 𝛼)qsolved + 𝛼qwave . (45)

The velocity ḣ of the body-fixed coordinate is transformed The level-set function is also blended in the region as
into the earth-fixed coordinate by matrices (40) and (41). follows:
Then, the displacement of a predictor step is calculated using
the third-order Adams–Bashforth scheme [11]:
𝜙 = (1 − 𝛼)𝜙solved + 𝛼𝜙wave . (46)
qwave and 𝜙wave are obtained from the equation which is based
𝛥t ( ̇ n+1 ) on the wave theory.
hn+1 = hn + 23h − 16ḣ n + 5ḣ n−1 . (43)
12 𝛼 is the weight function which is based on the distance
In a strong coupling which means that the displacement is from the outer boundary as follows:

) xm ≤ x ≤ xo ,
updated in pseudo-time steps, the displacement h is calcu- {
1.0 (
sin 2 x −x xb ≤ x ≤ xm
lated using the third-order Adams–Moulton scheme [11]: 𝛼x = (47)
n 𝜋 x−xb
𝛥t ( ̇ n+1 ) m b
hn+1 = hn + 5h + 8ḣ n − ḣ n−1 . (44)

) ym ≤ y ≤ yo .
12 {
1.0 (
sin 2 y −y yb ≤ y ≤ ym
The body motions are accounted for by the grid defor- 𝛼y = n 𝜋 y−yb (48)
mation, which depends on the weight function and the m b

distance from solid surfaces. The computational grids at

a certain amount of the distance from the solid surface Similar to the wave damping method, xb, and yb are the coor-
are deformed by the displacement. Then, the deforma- dinates from which the blending region starts, and xo, and yo
tion becomes zero with the function of the power of sin. are the locations of the boundaries. xm , and ym are arbitrary
The boundary condition on the body surface becomes the locations from which the constant value (1.0) starts.
body velocity, which is given by the grid velocity with

Table 1  Division numbers of the computational grids

Grid G3 G2 G1
IM × JM × KM IM × JM × KM IM × JM × KM

Bilge keel × 2 25 × 45 × 29 33 × 65 × 37 45 × 93 × 53
Aft fin × 2 37 × 45 × 31 49 × 65 × 41 65 × 93 × 53
Rect. fin1,2 37 × 25 × 45 53 × 33 × 65 73 × 49 × 89
Duct 45 × 45 × 25 61 × 65 × 37 85 × 89 × 57
Rudder 49 × 69 × 35 69 × 97 × 49 97 × 137 × 65
Rect. 73 × 57 × 51 101 × 81 × 73 145 × 117 × 105
Hull 81 × 85 × 57 113 × 121 × 77 161 × 169 × 113
Fig. 7  Computational grids

13
Journal of Marine Science and Technology (2019) 24:884–901 893

Table 2  Number of solved cells Grid Fine Medium Coarse

Total Solved Total Solved Total Solved

Grid G2
Bilge keel × 2 73,728 69,632 9216 8704 1152 896
Aft fin × 2 122,880 116,736 15,360 14,592 1920 1536
Rect. fin1 106,496 65,901 13,312 7878 1664 1
Rect. fin2 106,496 66,044 13,312 7889 1664 0
Duct 138,240 124,936 17,280 15,436 2160 1440
Rudder 313,344 297,152 39,168 36,996 4896 3894
Rect. 576,000 269,380 72,000 31,412 9000 242
Hull 1,021,440 908,720 127,680 111,726 15,960 11,526
Grid G1
Bilge keel × 2 210,496 193,600 26,312 24,200 3289 2662
Aft fin × 2 306,176 281,600 38,272 35,200 4784 3872
Rect. fin1 304,128 191,450 38,016 23,034 4752 278
Rect. fin2 304,128 191,819 38,016 23,102 4752 193
Duct 413,952 360,636 51,744 44,314 6468 4177
Rudder 835,584 792,246 104,448 98,598 13,056 10,706
Rect. 1,737,216 770,403 217,152 91,078 27,144 2097
Hull 3,010,560 2,644,082 376,320 326,184 47,040 35,488

3 Calculated results be updated. The FMG level is started with the medium grid
at the grid G2 . The solved cell number is remained at the
3.1 Steady double model flow coarse grid level with the grid G1, therefore, the FMG level
is started with the coarse grid at the grid G1.
3.1.1 Ship hull with several appendages with the FMG The overset update interval for the FMG method with
the overset grid method that the solutions are interpolated
The FMG and overset grid methods are applied to the flows and transferred based on the multigrid method is exam-
around a tanker hull shape [26] with several appendages. The ined at the grid G2 . Figure 8 shows the time history of the
Reynolds number is 3.0 × 106, and EASM [27] is applied to L2-norm of the pressure with changed in the overset update
calculate Reynolds stresses. The minimum spacing on the interval from 5 to 40. Figure 9 shows a comparison of the
wall satisfies y+ ≤ 1 in this case. Table 1 shows the division total resistance coefficients. The total resistance coeffi-
number of computational grids in each direction, and Fig- cients is non-dimensionalized by 12 𝜌U02 L02 . Spikes of the
ure 7 shows the aft part of the ship hull of the Grid G2 . The
grids are arranged with the priority of overset interpolation interval 5
in Table 1, and three grid sets G1, G2 and G3 are composed interval 10
interval 20
by the computational grids increasing
√ the division numbers interval 40
with the ratio approximately 2 considering the multigrid 10
-2

10-2
agglomeration of cells. As appendages, two bilge keels, two
aft hull fins, one duct situated just before a propeller, and a
L2 norm of p

-3
10
rudder are arranged with refined rectangular grids, which -3
10
resolve the flows around the aft fins and the rudder. The -4
10
computational grid consists of 10 grids with approximately 100 200 300

1.0 million cells at G3, 2.7 million cells at G2 and 7.6 mil-
lion cells at G1. 10-4
Three levels of the multigrid method are applied. Table 2
shows the total and solved cell number with the multigrid
level at the grids G2 and G1. The solved cell number of the 0 500 1000 1500
refined rectangular grid of the grid G2 becomes zero at the Step number
coarse grid level with the cell flagging rule of the present
method, thus the data of the refined rectangular grid cannot Fig. 8  Time history of the L2 norm (𝛥p)

13
894 Journal of Marine Science and Technology (2019) 24:884–901

L2-norm are observed when the interpolation and the cor- larger after 301 time steps, in which the computational grid
rection of the multigrid are performed with overset update changes to the fine grid using intervals of 20 and 40. The
intervals. The difference of the L2-norm between inter- total resistance coefficient shows a small difference with
vals of 5 and 10 is small, whereas the L2-norm becomes the overset update interval change until 300 time steps,
whereas a larger fluctuation with intervals of 20 and 40
appears after 301 time steps. From the above results, the
0.005 interval 5 overset update interval is 10.
interval 10
interval 20 Figures 10, 11 and 12 present comparisons of the results
interval 40 with and without the FMG using the grids G2 and G1. Addi-
0.004
0.002 tionally, the time histories of the total and frictional resist-
Total resistance coef.

ance coefficients without the FMG and MG methods at the

0.003 0.0015
grid G1 are shown in Fig. 11 for the comparison. The time
history of the L2-norm of the FMG present a spike at 301
0.001
time steps on the grids G2 and G1, and 601 time steps on the
0.002 100 200 300 400
grid G1 when the computational grid shifts from coarser grid
level to finer grid level (Fig. 10). Limited time is required
0.001 with the FMG method on the computation at coarser grid
level, and then the slope of the FMG becomes equivalent
to the result without using the FMG at the fine grid level
0
0 500 1000 1500 (Fig. 12). The steps from 500 to 800, where the total resist-
Step number
ance coefficient nearly converges on the grid G2 , and the
steps from 500 to 800 of the MG method and the steps from
Fig. 9  Convergence history of the total resistance coefficient 800 to 1100 of the FMG method on the grid G1 are magnified

Fig. 10  Comparison of the L2 Full-Multigrid

norm of 𝛥p (left: grid G2, right: Multigrid
Full-Multigrid
grid G1) Multigrid
-2 -2
10 10
L2 norm of p

L2 norm of p

-3 -3
10 10

10-4 10-4

Grid G2 Grid G1

0 500 1000 1500 0 500 1000 1500

Step number Step number

Fig. 11  Convergence history of Total Resi. Coef.(Full-Multigrid)

Total Resi. Coef. Multigrid)
the resistance coefficient (left: Total Resi. Coef.(wo Multigrid)
0.005 Total Resi. Coef.(Full-Multigrid) 0.006 Frictional Resi. Coef. (Full-Multigrid)
grid G2, right: grid G1) Total Resi Coef. (Multigrid) Frictional Resi. Coef. (Multigrid)
Frictional Resi. Coef. (wo Multigrid)
Frictional Resi. Coef. (Full-Multigrid)
Frictional Resi. Coef. (Multigrid) 0.005 0.00148
0.004
0.00146
0.00166

0.00144
0.004
Resistance coef.

Resistance coef.

0.00142
0.003 0.00164 800 900 1000 1100

0.00148
0.003
0.00146
0.00162
0.002 500 600 700 800 0.00144

0.002 0.00142
500 600 700 800

0.001
0.001

Grid G2 Grid G1
0 0
0 500 1000 1500 0 500 1000 1500
Step number Step number

13
Journal of Marine Science and Technology (2019) 24:884–901 895

Fig. 12  Comparison of the 10000 25000

elapsed time (left: grid G2, Full-Multigrid Full-Multigrid

right: grid G1) Multigrid Multigrid
8000 20000

Elapsed Time[sec.]

Elapsed Time[sec.]
6000 15000

4000 10000

2000 5000

Grid G2 Grid G1
0 0
0 500 1000 1500 0 500 1000 1500
Step number Step number

in Fig. 11. The total and frictional resistance coefficients frictional resistance coefficient from the total resistance
using the FMG and MG methods converge rapidly, thus, coefficient. The verification for the frictional and pres-
the present MG and FMG methods show the effectiveness. sure resistance coefficients is also performed. The fric-
The elapsed time at 650 steps on the grid G2, where the total tional resistance coefficient of each grid is 1.560 × 10−3 ,
resistance coefficient converges within ± 0.5% from the aver- 1.189 × 10−3 and 1.078 × 10−3 from the grid G3 to the grid
aged value are compared. The FMG method takes 3271 s, G1. The order of accuracy p is estimated as 3.482 and again
whereas the multigrid method (FMG vs. MG) takes 5118 s, the theoretical accuracy pth is assumed 2.0. The uncer-
thus the FMG succeeds in reducing the computational time tainty USN is estimated USN %S1 = 16.2 by the uncorrected
by approximately 37%. Additionally, the elapsed time at 650 approach and USN %S1 = 5.90 with the correction factor.
steps on the grid G1 with the MG method and at 1050 steps The uncertainty USN %S1 is estimated 60.4 using the factor
with the three level FMG method where the total resistance of safety approach [30]. Similarly, the pressure resistance
coefficient converges within ±0.5% from the averaged value coefficient of each grid is 4.210 × 10−4 , 4.410 × 10−4 and
are also compared. The FMG method takes 7628 seconds, 3.660 × 10−4 . Since the the pressure resistance coefficient
whereas the multigrid method takes 9840 seconds on the exhibits the oscillatory divergence, the grid uncertainty
grid G1. Consequently, the FMG method reduces the elapsed cannot be obtained.
time by approximately 23% even in the case of using finer
grids. 3.1.2 Propeller working condition
The verification based on the generalized Richardson
extrapolation methods [28–30] is performed using the The present method with the propeller model is applied to
three grids G3, G2 and G1. Since the resistance coefficients the flows around a hull with an appendage [4]. A stern duct
converge to the fourth digit, the iterative error can be neg- is equipped to the hull as an energy-saving device, and
ligible in comparison with the grid uncertainty UG . Con- the effect of the duct is also examined. Table 3 shows the
sequently, the uncertainty USN is estimated using the grid division number of computational grids in each direction,
uncertainty UG . The total resistance coefficient of each and Fig. 13 shows the grids near the aft hull. The three
grid is 1.981 × 10−3 , 1.630 × 10−3 and 1.444 × 10−3 from grids, which include the duct, stern tube and hull, are over-
the grid G3 to the grid G1. The total resistance coefficient lapped with the overset method. The experimental results
shows the monotonic convergence, and the estimated order and computational conditions are available on the work-
of accuracy p becomes 1.832. When the theoretical accu- shop websites. The Reynolds number is 7.248 × 106 , and
racy pth is assumed 2.0, the uncertainty USN is estimated EASM [27] is applied. The minimum spacing on the wall
USN %S1 = 17.75 based on the uncorrected approach [28] satisfies y+ ≤ 1 which is derived from the restriction of
and USN %S1 = 1.90 with the correction factor [28]. Mean- the low-Reynolds number turbulence model. The propeller
while, USN %S1 is 24.3 using the factor of safety method diameter is non-dimensionalized by the ship length, and
[30] based on the parameter P = p∕pth . Further discus-
sions may be required on the scatter of the uncertainties
among various verification methods. The uncertainty USN Table 3  Division number of the Grid IM × JM × KM
computational grid
values vary between USN %S1 = 1.90 and USN %S1 = 24.3 in
Duct 225 × 129 × 81
the present case. The time history of the frictional resist-
Stern tube 133 × 145 × 81
ance coefficient can be found in Fig. 11, and the pressure
Hull 101 × 241 × 169
resistance coefficient can be evaluated by subtracting the

13
896 Journal of Marine Science and Technology (2019) 24:884–901

Figures 14 and 15 show the axial velocity contours and

cross-flow vectors in the propeller working condition with-
out the stern duct. From Figs. 14, 15, 16 and 17, the figure
on the left-hand side is the measured result, and the figure
on the right-hand side is the calculated result. The body
forces, which express the propeller effect, accelerate the
axial velocity and induce the rotational vectors. Figures 16
and 17 show the same variables with and without the stern
duct. The computational results can reproduce the aspect in
which the acceleration of the axial velocity with the stern
duct becomes lower than the case without the stern duct
and the rotational components of the right-hand side and
lower part becomes smaller than the rotational components
without the stern duct.
Table 4 presents a comparison of the propeller rotational
speed based on the dimensionalized value, thrust and torque
Fig. 13  Computational grid coefficients for the condition without the stern duct. Table 5
presents a comparison of the condition with the stern duct.
non-dimensionalized diameter is 0.029, and the propeller The computation can capture the condition in which the pro-
rotational speed is balanced to equal the resistance at the peller rotational speed with the stern duct becomes slower
ship point which is defined at the 2015 workshop [4]. than the speed without the stern duct and the thrust and

Fig. 14  Comparison of the axial

velocity for propeller working
conditions without the stern
duct

Fig. 15  Comparison of the

cross-flow vectors for propeller
working conditions without the
stern duct

13
Journal of Marine Science and Technology (2019) 24:884–901 897

Fig. 16  Comparison of the

axial velocity for the propeller
working conditions with the
stern duct

Fig. 17  Comparison of the

cross-flow vectors for propeller
working conditions with the
stern duct

Table 4  Comparison of the propeller characteristics without the stern Table 6  Wavelengths and heights of regular head waves
duct
Case Wavelength Wave height
Variable Meas. Comp.
C1 0.6506 0.0102
Prop. rot. speed 7.8 7.71 C2 0.8507 0.0128
Thrust coef. 0.217 0.217 C3 1.1497 0.0203
Torque coef. 0.0279 0.0282 C4 1.3708 0.0245
C4 1.9505 0.0323

Table 5  Comparison of the propeller characteristics with the stern

duct 3.2 Unsteady flow with regular head waves
Variable Meas. Comp.
Finally, the case of an unsteady flows is examined. The
Prop. rot. speed 7.5 7.59 unsteady computation with a main hull and a rudder in regu-
Thrust coef. 0.233 0.226 lar head waves is selected. The benchmark measured data are
Torque coef. 0.0295 0.0289 available from the workshop site [4]. The Reynolds number
is 1.074 × 107, the Froude number is 0.261 and the k − 𝜔SST
torque coefficients become larger. The present results are model is applied. The wavelengths and heights of the regular
consistent with the measured results, and the effectiveness head waves are shown in Table 6. All data are non-dimen-
of the flows is consistent with the energy-saving device. sionalized by the ship length. The wave height increases

13
898 Journal of Marine Science and Technology (2019) 24:884–901

Table 7  Division number of the computational grids with the wavelength. The computational grids are refined
Grid IM × JM × KM rectangular grids designed to the flows around the rudder,
the rudder and hull grids, and the background rectangular
Refined rect. 49 × 49 × 49 grid. The three rectangular grids are chosen in accordance
Rudder 49 × 89 × 33 with the wavelength. Table 7 shows the division number of
Hull 257 × 161 × 161 computational grids in each direction. The grids are tabled
Rect. (C1) 185 × 89 × 97 in the order of the priority of the overset grid method. The
Rect. (C2,C3) 153 × 89 × 97 rectangular grid for Case C1 has the largest division number
Rect. (C4,C5) 117 × 89 × 97 IM to resolve the shorter wavelength, and at least 40 points

Fig. 18  Instantaneous view of

the free surface around a hull

13
Journal of Marine Science and Technology (2019) 24:884–901 899

lie in each wavelength that the numbers of points are larger the background rectangular grid to the hull grid, and the
than the number of Ref. [10]. Additionally, more than ten waves generated by the interaction between the incoming
grid points rest in the wave height. All of the grids except the waves and the hull are also observed in the background rec-
background rectangular grids are maintained, and the overset tangular grid. Thus, the present overset grid method shows
information is generated by replacing the three background the effectiveness of the unsteady computation with the
rectangular grids. A wall function is applied to the body incoming waves and ship motions.
surfaces in all cases for the stable computation, even in the Figure 19 shows a comparison with an added resistance
case of higher wave heights C4 and C5. coefficient. The added resistance coefficient is introduced by
Heave and pitch motions are free in all cases, and the non- subtracting the resistance associated with calm water con-
dimensionalized temporal time steps Δt is set to 0.005. The ditions from the resistance associated with waves. In other
regions where the overset relations are composed deform words, the added resistance coefficient indicates an increase
with the body motions to maintain the overset information, in resistance from the calm water condition. The present
and the amount of deformation gradually decreases with computational results show reasonable consistency with the
distance from the body surface. Such method is adapted to measured results and clearly show the tendency of the peak
avoid the computational load using the dynamic overset grid value at 𝜆∕L = 1.1497.
method. The first amplitudes of the Fourier analysis for the ship
Figure 18 shows the free surface near the hull. The motions are compared in Fig. 20. The current results are con-
incoming wave height becomes larger as the wavelength sistent with the measured results. The amplitude of the heave
increases. The wave height contours smoothly connect from motion has a peak value at 𝜆∕L = 1.1497, where the added
resistance coefficient also has a peak value. The amplitude of
3 the pitch motion becomes larger as the the wavelength ratio
Meas. increases. Figure 21 shows the comparison of the flow field
+ Comp.
of the condition 𝜆∕L = 1.15 at the position x∕L = 0.9825
2.5 +
which is the propeller plane in the earth fixed coordinate
2
according to the quarter of the encounter period. The cal-
culated results reproduce the time-dependent change of the
RAW/4 g 2B2/L

+
axial flow including the effect of the incoming waves and
1.5
ship motions, and the cross flow vectors of the calculated
+ results trace the aspect of the SPIV measurement[31].
1

+
0.5 + 4 Conclusions
0
0.5 1 1.5 2 A new structured grid solver that can manage the over-
/L set grid method is developed. The governing equation is
based on the artificial compressibility approach, and the
Fig. 19  Comparison of an added resistance coefficient dual time-stepping method is applied to calculate unsteady

Fig. 20  Comparison of the first 1

+
1.5

amplitude of motion +
Meas.
+ + Comp.

+
1

+
0.5
/k

+
z/

0.5

+ Meas.
+ + Comp. +
0
0.5 1 1.5 2
0
0.5
+ 1 1.5 2
/L /L

Heave motion Pitch motion

13
900 Journal of Marine Science and Technology (2019) 24:884–901

flows. The grid velocity is included in the convective

term with the moving grid approach for body motions.
The FMG technique is extended to the overset grid sys-
tem. Cell flagging for the overset procedure is transferred
from the fine grid to the coarse grid. The prolongation
and restriction procedures of the multigrid method are
performed in the coarse grid level computations of the
FMG stage with an overset update interval. The method
to eliminate the body lapped region caused by the overset
grid method is introduced and shows the effectiveness for
a typical sample.
The overset update interval for the FMG with the overset
grid method is examined for the case of a ship hull with sev-
eral appendages. The comparisons of the L2-norm and the
convergence history of the resistance coefficient show that
the overset update interval of 10, which is the interval of the
solution transfer and interpolation, shows good performance.
The FMG is examined on the two grid sets which have the
different grid sizes, and the present FMG method succeeds
in reducing the computational time by approximately 37%
on the medium grid size and 23% even in the case of using
finer grid compared with the results without using the FMG.
The fundamental verification study is performed using the
three grid sets increasing
√ the division numbers with the ratio
approximately 2 . The total and frictional resistance coef-
ficients of this case take the monotonic convergence, and the
uncertainties USN are obtained using the three different ways
based on the generalized Richardson extrapolation method.
The proposed solver with the propeller and explicit
algebraic stress turbulence models is applied to flows with
the energy-saving device. The aspects of the experimental
data that correspond to the propulsion condition are repro-
duced in the present computations. Thus, the coupling of
the propeller model with the solver via body forces pro-
duces good results, and has the ability to estimate ship per-
formance, including the effects of energy-saving devices.
Finally, the unsteady flows around a container ship with
ship motions in regular head waves are analyzed. The cell
size of the background rectangular grids, which generates
regular head waves, is adjusted to the wavelength and over-
lapped based on the wavelength. The computational grid is
modified by replacing the rectangular grids. Regular waves
are generated in a region designated near the boundaries
of the computational domain. The incoming wave height
increases as the wavelength increases as intended. The
amplitude of the ship motion, which is obtained by solv-
ing the motion equation, is compared to the experimental
data, and the results are consistent. The regions in which
the overset relations are composed deform with the body
motions to maintain the overset information. The present
method works well and can avoid the dynamic overset
Fig. 21  Comparisons of flow fields according to the encounter period method and the computational load. The added resistance
(left: measured data, right: calculated data) coefficient, which is introduced by subtracting the resistance

13
Journal of Marine Science and Technology (2019) 24:884–901 901

with a calm water condition from the resistance with waves, motions of KVLCC2 with fixed and free surge in short and long
is also consistent with the measured results. Additionally, the head waves. Ocean Eng 59:240–273
14. Broglia R, Dubbioso G, Durante D, Di Mascio A (2012) Simu-
comparison of the flow field is performed. The calculated lation of turning circle by CFD: analysis of different propeller
results reproduce the time-dependent change of the axial models and their effect on manoeuvring prediction. Appl Ocean
flow and cross flow vectors comparing with the results of Res 39:1–10
the SPIV measurement. 15. Di Mascio A, Broglia R, Muscari R (2006) On the application of
the single-phase level set method to naval hydrodynamic flows.
Comput Fluids 36:868–886
Acknowledgements The authors would like to express the deepest 16. Broglia R, Durante D (2017) Accurate prediction of complex
appreciation to Prof. Toda of Osaka Univ. for providing the results of free surface flow around a high speed craft using a single-phase
the measurement. This work was supported by JSPS KAKENHI under level set method. Comput Mech. https​://doi.org/10.1007/s0046​
Grant Number JP16K06919. 6-017-1505-1
17. Korkmaz KB, Orych M, Larsson L (2015) CFD predictions
including verification and validation of resistance, propulsion
and local flow for the Japan bulk carrier (JBC) with and without
an energy saving device. In: Proc. Tokyo 2015 workshop on CFD
References in ship hydrodynamics
18. Kobayashi H, Kodama Y (2016) Developing spline based over-
1. Hino T, Hirata N, Ohashi K, Toda Y, Zhu T, Makino K, Takai set grid assembling approach and application to unsteady flow
M, Nishigaki M, Kimura K, Anda M, Shingo S (2016) Hull form around a moving body. J Math Syst Sci 6:339–347. https​://doi.
design and flow measurements of a bulk carrier with an energy- org/10.17265​/2159-5291/2016.09.001
saving device for CFD validations. Proc. PRADS 19. Ferziger JH, Perić M (2012) Computational methods for fluid
2. Hino T (2005) Proc. of CFD workshop Tokyo dynamics, 3rd edn. Springer, Berlin
3. Larsson L, Stern F, Visonneau M (2014) Numerical ship hydro- 20. Farmer J, Martinelli L, Cowles G, Jameson A (1995) Fully non-
dynamics: an assessment of the Gothenburg 2010 workshop. linear CFD techniques for ship performance analysis and design.
Springer, Berlin. https​://doi.org/10.1007/978-94-007-7189-5 AIAA-95-1690-CP
4. Larsson L, Stern F, Visonneau M, Hirata N, Hino T, Kim J (2015) 21. Roe PL (1986) Characteristic-based schemes for the Euler equa-
Proc. Tokyo 2015 workshop on CFD in ship hydrodynamics tions. Ann Rev Fluid Mech 1986:18
5. Queutey P, Visonneau M (2007) An interface capturing method 22. Anderson WK, Thomas JL, van Leer B (1986) Comparison of
for free-surface hydrodynamic flows. Comput Fluids 36:1481– finite volume flux vector splitting for Euler equations. AIAA J
1510. https​://doi.org/10.1016/j.compf​l uid.2006.11.007 24(9):1453–1460
6. ReFRESCO. http://www.refre​sco.org 23. Hino T (1999) An interface capturing method for free surface
7. Kim SE, Rhee B, Shan H, Gorski J, Paterson E, Maki K (2010) flow computations on unstructured grids. J Soc Nav Archit Jpn
Prediction of turbulent free-surface flows around surface ships 186:173–183
using a scalable unstructured finite-volume based RANS solver. 24. Peng T, Merriman B, Osher S, Zhao H, Kang M (1999) A PDE-
In: Proc. Gothenburg 2010 workshop on CFD in hydrodynamics based fast local level set method. J Comput Phys 155:410–438
8. Jasak H, Vukčević V, Christ D (2014) Rapid free surface simula- 25. Ohashi K (2018) Development of numerical method to simulate
tion for steady-state hull resistance with FVM using OpenFOAM. flows around a ship at propulsion conditions in regular waves
In: Proc. of 30th symposium on naval hydrodynamics, pp 548–554 coupling with the ship propulsion plant model. Appl Ocean Res
9. Vukčević V, Jasak H, Malenica S (2016) Decomposition model for 73:141–148. https​://doi.org/10.1016/j.apor.2018.02.006
naval hydrodynamic applications. Part II: Verification and valida- 26. Proc. of CFD workshop Tokyo 1994
tion. Ocean Eng 121:76–88 27. Rumsey CL, Gatski TB (2001) Recent turbulence advances applied
10. Carrica PM, Wilson RV, Stern F (2007a) An unsteady single- to multielement airfoil computations. J Aircr 38(5):904–910
phase level set method for viscous free surface flows. Int J Numer 28. ITTC-Recommended Procedures and Guidelines 7.5-03-01-01
Methods Fluids 53:229–256 (2017) Uncertainty analysis in CFD verification and validation
11. Carrica PM, Wilson RV, Noack RW, Stern F (2007b) Ship motions methodology and procedures
using single-phase level set with dynamic overset grids. Comput 29. Wilson R, Shao J, Stern F (2004) Discussion: criticisms of the
Fluids 36:1415–1433 correction factor verification method. J Fluid Eng 126:704–706
12. Noack RW (2005) SUGGAR: a general capability for moving 30. Xing T, Stern F (2010) Factors of safety for Richardson extrapola-
body overset grid assembly, AIAA paper 2005-5117. In: 17th tion. J Fluid Eng 132(6):061403
AIAA computational fluid dynamics conference 31. Wu PC, Hossain MA, Kyaw HA, Matsumoto A, Toda Y (2018)
13. Sadat-Hosseini H, Wu PC, Carrica PM, Kim H, Toda Y, Stern F Forces, ship motions and nominal wake for KCS ship model in
(2013) CFD verification and validation of added resistance and regular head waves. In: Proc. of JASNAOE Spring Meeting 25

13