CC BY-NC-ND
0+ MODEL
Available online at www.sciencedirect.com
ScienceDirect
Publishing Services by Elsevier
International Journal of Naval Architecture and Ocean Engineering xx (2017) 1 — 15
http://www.journals.elsevier.com/international-journal-of-naval-architecture-and-ocean-engineering/
Development of a numerical simulation tool for efficient and robust
prediction of ship resistance
Geon-Hong Kim*, Sanghoon Park
Hyundai Maritime Research Institute, Hyundai Heavy Industries Co., Ltd, Republic of Korea
Received 7 June 2016; revised 5 December 2016; accepted 8 January 2017 Available online ■ ■ ■
Abstract
In this paper, a two-phase flow solver HiFoam® has been developed based on the OpenFOAM® to predict resistance of a ship in calm water. The VOF method of OpenFOAM® was reviewed and a simple flux limiter was introduced to enhance the robustness of the solver. The procedure for predicting ship motion was modified by introducing Quasi-Steady Fluid-Body Interaction (QS-FBI) with least square regression to improve the efficiency. Other minor factors were considered as well in terms of the efficiency and robustness. The HiFoam was applied to KCS and JBC simulations to validate its efficiency and accuracy by comparing the results to experimental data and STAR-CCM+. The HiFoam® was also applied to various ships and it showed good agreements to the experimental data.
Copyright © 2017 Society of Naval Architects of Korea. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Keywords: OpenFOAM; Resistance; CFD; VOF; QS-FBI; Flux limiter for VOF
1. Introduction
Computational Fluid Dynamics (CFD) has advanced rapidly in recent years and become as one of the most important technique in engineering fields. Especially the CFD technique plays an important role in shipbuilding industries by replacing experiments successfully at the early design stages. Potential flow codes have been applied to many marine hy-drodynamic problems such as resistance, seakeeping, propeller and wave load analysis since it provides solutions in a few minutes with moderate accuracy (Newman, 1992). As an interest in environment grows and the crude oil price rises, the demand on a ship with lower emission and fuel consumption has been increased as well. Consequently, recent ships have more complicated geometries with various appendages that can be hardly handled by the potential codes. Such a problem has been remedied by applying viscous solvers using free
* Corresponding author. Peer review under responsibility of Society of Naval Architects of Korea.
surface capturing techniques such as Volume of Fluid (VOF) or level set method. A number of general CFD tools like CFX, Fluent, and STAR-CCM+ have been introduced and applied to ship hydrodynamics problems successfully.
It is well known that such viscous solvers provide reasonable solutions to general CFD problems. Furthermore, as computing hardware and techniques on parallel computing developed very rapidly, CPU time for carrying out a simulation using the viscous solver has been significantly reduced. However, it still takes too long time to be practically used in industrial fields. For instance, a trim optimization analysis requires hundreds of simulations, which should be completed within limited time. If a simulation of a single case takes an hour, it will take more than 4 days to complete 100 simulations. A statistical technique such as Design of Experiment (DOE) might be considered to reduce the time but the computational efficiency should be considered prior to all other factors. Robustness is another important property of a CFD tool. If a CFD tool is too sensitive and easily diverged, it cannot be trusted and used in an industrial field though it provides highly accurate solutions.
http://dx.doi.org/10.1016/j.ijnaoe.2017.01.003
2092-6782/Copyright © 2017 Society of Naval Architects of Korea. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
Open source techniques have been attracted a lot of attention for the last decade since one can access its source codes to use or modify without paying for it. Many open source CFD codes like SU2 (Palacios et al., 2013, 2014) and PyFR (Witherden et al., 2014) have been introduced to CFD field. The OpenFOAM® is one of the most popular open source CFD codes since it was released to public in 2004. It provides libraries and applications with source codes written in C++ language for the solution of continuum mechanics problems. For ship hydrodynamics problems, the OpenFOAM® provides interFoam, a multi-phase solver using interface capturing method based on VOF scheme. It has been studied to predict the resistance of a ship on water including heave and pitch motions accurately by using the OpenFOAM® (Park et al., 2013; Lee, 2013), and they showed that the errors of resistance prediction lie within the range of ±2% relative to experimental data and other commercial CFD tools like STAR-CCM+ and FLUENT. Based on the feasibility studies on the OpenFOAM®, Hyundai Heavy Industries (HHI) began in earnest to use the OpenFOAM® for evaluating performance of a ship numerically since 2014 and developed HiFoam® (Lee, 2014), which is a program for carrying out CFD simulations automatically from pre to post processing based on the interFoam of OpenFOAM-2.1.x.
As we mentioned earlier, the accuracy of OpenFOAM® on predicting ship resistance has been validated a lot. However, it is not as robust and efficient as other commercial codes yet which are essential for being widely used in a shipbuilding company. In this paper, we focused on how to enhance robustness and efficiency of the OpenFOAM® solver, especially HiFoam, to make it comparable to commercial codes without losing accuracy.
2. Numerical methods
2.1. Governing equations
interFoam, a solver for two incompressible, isothermal immiscible fluids using a VOF phase fraction based on interface capturing approach, solves following continuity, transport of phase-fraction, and momentum equations.
P = aPi + (1 - a)Pg
m = am? + (1 — a)mg
where the subscripts I and g represent liquid and gas, respectively. T in Eq. (3) is the deviatoric stress tensor where the stress term can be expressed as following for an incompressible flow.
T = 2mS
S is the mean
rate of strain tensor defined as S = 0.5[Vu + (Vu)T]. For Newtonian and incompressible flows, the stress term in Eq. (3) can be further decomposed with the aid of identity of the divergence operation in vector calculus.
VT = v-(mVu)+Vu$Vm
fb in Eq. (3) is a body force term including gravity and surface tension, where the surface tension can be evaluated by using the Continuum Surface Force (CSF) model of Brackbill et al. (1992).
fs = skVa
where k is mean curvature of the free surface, obtained from following relation.
It is common to use modified pressure in VOF method for a convenience in applying boundary conditions for the pressure since it contains hydrostatic components.
pd = p — Pg-x
where g and x are gravity and position vectors, respectively. The momentum Eq. (3) may be re-written for a two-phase flow as following by using the above relations.
9(pu) 9t
+ V-(puu) - V-(mVu) - Vu-Vm -Vpd — g -xVp + skVa
V- u 0
Va „ . — + V- (ua) = 0
9 (pu) Vt
+ V-(puu) = —Vp + V-T + pfb
where a is phase (or volume) fraction that measures how much portion of a cell is filled with the heavier fluid, usually water. Value of a lies within the range of [0, 1] where 0 represents fully gaseous cell while 1 represents a cell fully immersed in the liquid. The properties of a fluid such as density and viscosity can be estimated by using the phase fraction as following.
2.2. Numerical solution procedure
The governing equations described in the preceding section are discretized by means of finite volume method on an unstructured grid. The van Leer scheme (van Leer, 1979) and linear interpolation scheme are used to discretize the convection and diffusion terms, respectively. The Local Time Stepping (LTS) scheme was selected for temporal discretization for an efficient quasi-steady simulation. The gradients of flow properties can be estimated by using either linear interpolation or least square method.
The interFoam uses Pressure-Implicit with Splitting of Operator (PISO) (Rhie and Chow, 1982) algorithm for the pressure—velocity coupling procedure, which is shown in
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
Fig. 1. Flow chart of the pressure—velocity coupling procedure of solver interFoam
Fig. 1. Once the phase fraction equation is solved by using face fluxes of the previous time step, pressure, velocity fields, and face fluxes are corrected based on the new phase fraction. However, the phase fraction is not updated after the face fluxes are corrected since it is excluded from the outer corrector, as shown in Fig. 1. As a result, it might cause an inconsistency in p fraction and pressure fields when a transient simulation is carried out. However, we did not modify the velocity-pressure coupling algorithm in the present work because a quasi-steady solution is of interest. It is worthwhile to note that Fig. 1 represents the PISO algorithm that is used in solvers of OpenFOAM®-2.1 .x where the step for updating phase fraction is included to the outer corrector loop in the later versions.
3. Improvement of efficiency and robustness
researchers have studied to tackle such problems. They have developed and applied additional source term (Rusche, 2002; Park and Rhee, 2014) and high-resolution schemes such as CICSAM (Ubink, 1997) and HRIC (Muzaferija et al., 1999) for suppressing numerical diffusion and predicting sharp interfaces. The interFoam uses following phase fraction equation instead of Eq. (2), which is re-arranged by Weller (2002).
9a / \ r — + V-(ua)+V-[ura(1 - a)] = 0
where ur = u — ug is the relative (or compression) velocity that can be modeled at cell faces as follows.
urf = nf min
where Ca is a model coefficient, f is volume flux and nf is face flux of the gradient of phase fraction defined as following.
The time step size, however, for updating the phase fraction field is strictly limited in terms of CFL number, where the value of phase fraction is expected to be bounded within the range of [0, 1]. A simple one-dimensional system was considered as shown in Fig. 2 to investigate the effect of time step size on phase fraction. The half of domain is initially filled with water, where the water is flowing into the domain from the left face and air flows out through the right face. The flow speed is assumed to have constant value of unity and the domain is subdivided into five cells with uniform width for simplicity. The free surface is located on cell 2, and we call this cell as an interface (or mixture) cell. The rest cells are fully filled with a single-phase fluid and we will refer those cells as fully wetted/dried cells depends on their content hereafter.
The flow of Fig. 2 is initially driven by convection only since no gradient of both pressure and velocity exist at the moment, thus the momentum equation can be simplified as following.
9(pu) 9t
-V-(puu)
3.1. Review on VOF method
In this section, we would like to discuss the improvement of robustness of a VOF simulation. In simulations of ship hydrodynamics, it is very important to capture an interface between two different mediums accurately and efficiently. The interface is captured by solving the phase fraction Eq. (2) which is mathematically hyperbolic. It is well known that such a conservation equation suffers from numerical diffusion when first order accurate scheme is applied while classical higher order schemes exhibit spurious oscillations near the discontinuous solutions (Sweby, 1984). Many earlier
freesurface
Cell 0 1 3 4
a 1 1 0.5 0 0
Fig. 2. Simple one-dimensional domain for investigating the behavior of phase fraction.
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
It is easy to infer from the above equation that the flow is accelerated if the density changes across a cell. The flow velocity at cell 3 in Fig. 2 might be approximated as follows.
,n+1 un
where Dpf 3 = pf 23 — Pf 34 is a difference in density across the cell 3 and the subscript f,23 represents the value on the face between cell 2 and 3. The face values can be estimated by using linear interpolation. For a simulation of water—air mixture, the density of water is almost 1000 times larger than air and it causes large difference in Dp" 3 thus increase in time step size Dt may cause significant change in speed of light fluid in the vicinity of interface. To eliminate such unexpected behavior, the flux limiter was applied on faces at which interface cell and single-phase cell meet each other. The basic idea of the limiter is that the main driver of flow acceleration is undesirable interpolation of density on faces of a single-phase cell where the density on its faces should be identical to the phase of the corresponding cell. Thus, an interpolated value on the face was replaced to a physically correct value. To accomplish this, interface cells were distinguished from fully wetted/dried cells once the phase fraction field is updated based on the given flow field using MULES (Multidimensional Universal Limited Explicit Solver). Then the phase fraction weighted flux fa is limited on the faces by using the following equation.
(fa) = ff -pos( af — 0.5)
where af is interpolated phase fraction on a face f. The pos function gives 1 if the value inside the parenthesis is positive and 0 if negative. The above limiter was tested by implementing to the interFoam and running a dam break simulation. Fig. 3 shows the resultant flow fields with and without flux limiter. The free surface is represented by white solid line and water is flowing from left to right of domain. We see that, in Fig. 3(b), the flow across the free surface further accelerates if flux limiter is not applied, as it is explained earlier, while the flow with flux limiter doesn't. Fig. 3(c) and (d) shows detailed flow fields of each case and high velocities in air region near free surface are observed. The present flux limiter is very simple but successfully distinguishes flow regions of light fluid from heavy fluid.
For the case of ship hydrodynamics, the flux limiter (17) helps to enhance the robustness at higher CFL number by suppressing the speed of light fluid near the free surface. Fig. 4(a) shows force distributions of KCS simulations with and without the flux limiter with Coa = 1.0 by using MULSE of OpenFOAM-2.1.x. The simulation without flux limiter diverged after a certain iterations while the other one with flux limiter converged. Fig. 4(b) compares the velocity contours near free surface and one can easily see very high speed at dry cells adjacent to the free surface when flux limiter is not applied. However, such jump of velocity near the free surface has been successfully suppressed by applying the flux limiter.
aoJ-J-1—.-J—
0 0.06 0.1 0.16 0.2 0.25 0.3
Distance, m
(a) Computational Domain (b) Velocity magnitude distribution on curve B
(c) Detail view of A, without flux limiter (d) Detail view of A, with flux limiter
Fig. 3. Comparison of results of dam break simulations with and without flux limiter.
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
(b) Velocity contours
Fig. 4. Comparison of the effect of flux limiter on ship hydrodynamic simulations.
3.2. Improved motion prediction procedure
A ship changes its attitude while it moves forward on water mainly due to change of pressure acting on its surface. Although the amount of change in attitude of a ship in dynamic state is relatively small, it causes change of pressure field, which affects to the resistance of a ship. Therefore, it is important to predict attitude accurately in dynamic state.
A model for mesh motion similar to Dynamic Fluid-Body Interaction (DFBI) of STAR-CCM+ (Muzaferija et al., 1999) has been applied but the motion is estimated not at every time step but after an interval of certain time steps. The motion is frozen during the period and the flow fields are updated only. The present model takes the interaction between fluid and body into account in a quasi-steady manner, thus we refer to such a model as Quasi-Steady Fluid-Body Interaction (QS-FBI). Lee (2013) showed that this method successfully collaborated with the LTS schemes and well predicted the motion of a ship, however, it takes relatively long time to predict the motion since it interactively finds the equilibrium state of the ship.
We tried to find an efficient way to predict the motion of a ship during a resistance simulation, and we found that the pitching moment My and the vertical force balance DFz, which are drivers of the ship motion, almost linearly varied as trim angle q and sinkage distance s are changed, as shown in Fig. 5(a) and (b). Thus, the least square method was employed to predict the final attitude of a ship in equilibrium state with less motion steps. The objective of the present method is to find q and s subject to following conditions.
My(q) = 0 (18)
Fz(s)-W = DFz(s)=0 (19)
where W is the weight of a ship. First we evaluate moments and forces at several different q and s conditions. Then, the new q and s are estimated to satisfy the above relations by means of a simple linear regression as follows.
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
Fig. 5. Histories of vertical force, pitching moment, and ship motions.
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
qnew _
q2 - q2 eMl-WM,
■My + q
sAFv- sDFz
--DFz + s
(20) (21)
where overbar denotes averaged value. To implement the above method to the present code, additional variables for storing the histories of 6 and s are required as well as corresponding moment and force histories. It might be cumbersome to modify the source code but it does not demand moment of inertia as an input when a simulation is carried out, since the trim angle is no longer evaluated from the motion of equation. Furthermore, the weight of a ship is estimated by integrating hydrostatic pressure on ship surfaces before the simulation begins. Thus, the present simulation does not ask any parameters related to the ship motions and it makes dynamic simulations much easier.
The present method was verified by applying to a simulation on KRISO Container Ship (KCS) at a speed of 24 knots in calm water with free heave and pitch motions (Van et al., 1998). Fig. 5(a) and (b) show AFZ with s and My with 6,
respectively. Those figures also compare results of the present method to those of the original QS-FBI method. Despite the motions were not determined by solving the motion of equations, the present method predicts the force balance and moment close to the QS-FBI results. Fig. 5(c) compares the motion histories with motion steps and we see that the present method saved about 38% of motion steps to predict the equilibrium state of attitude, where one motion step required 200 time steps. The original QS-FBI predicted smooth variation of ship motions while the present method showed more steep and sudden change of motions. The present method even showed overshooting behavior in prediction of trim angle, as shown in Fig. 5(c).
The efficiency of the present method in predicting ship motion is more evident when it is compared to the results of dynamic mesh simulations. To demonstrate this, a simulation using interDyMFoam of OpenFOAM-3.0.x with Co = 100 and Coa = 25 was conducted. Fig. 6(a) shows the force distribution with physical time and we saw that the dynamic mesh simulation should be carried out at least 80 s in physical time to get converged solution which took about 16,000 s of CPU time, while the present method took about 6000 s of CPU
(b) Trim angle distributions 6. Comparison of dynamic mesh and QS-FBI method for predicting ship attitudes.
S^ew _
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
time only. It is worthy to note that the present method is a pseudo-steady simulation and it is not plausible to compare the results to those of dynamic mesh simulation directly. Therefore we compared the predicted motion with respect to CPU time instead of physical time or iteration steps, as shown in
Fig. 6(b).
3.3. Treatment of corner flows
It is well known that sharp edges or corners are major sources of instability if the geometry is not properly modeled. We have also suffered from such problems when we carried out simulations using the OpenFOAM®. Fig. 7 shows the
(b) Stern aperture of twin LNGC Fig. 7. Geometrical sources of instability in ship hydrodynamics simulations.
corner and sharp edge where such unstable behavior is often observed in ship hydrodynamics analysis. It is obvious that a fluid is accelerated when it passes over the corners or sharp edges. If the grid in this region is distributed inappropriately to resolve such a sudden acceleration, pressure may be significantly dropped more than it is expected to be since the original OpenFOAM® solvers does not include any limiter to suppress such unphysical phenomena. Then computations easily diverged due to those regions as shown in Fig. 7.
To tackle such a problem, the corner flows were reviewed from the perspectives of geometry and numerical limiter. It is obvious that the grid size in flow direction should be small enough to prevent significant variations in pressure and velocity. It is also important how to divide the space near the corner into small finite volumes when the grid is generated. Fig. 8 shows two different types of grid arrangements.
Cell P in Fig. 8(a) meets cell N directly and the middle point between the two cells is not coincide with the face center. Even the middle point is not located on the face f and it can be often found in practice. Such arrangement of cells increases face skewness, which may harmfully affect to the
Fig. 8. Grid arrangements near corners.
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
solution. Fig. 8(b) shows an intermediate cell connecting two cells P and N adjacent to the corner. The intermediate cell may increase the number of cells but it plays an important role in preventing highly skewed faces, and therefore it can reduce error due to improper interpolations at corresponding cells. Fig. 9 compares the resultant pressure contour near the corner where a hull and transom meet to each other. We see that highly skewed face in Fig. 9(a) caused larger pressure drop in the corresponding region while intermediate cell seems to suppress such significant drop in this region, as shown in Fig. 9(b).
It is not plausible, however, to control grid arrangement near the corners or sharp edges in industrial practices. Another candidate to control the instability due to geometry is a numerical limiter. The basic idea for adopting a limiter is that the instability is mainly caused by overly accelerated flow near corners/sharp edges. Such flows generate excessive mass flow on corresponding cells and this leads unphysical pressure variation in pressure—velocity coupling procedure. To prevent such unexpected behavior, the overly accelerated flows might be clipped by using the clipping procedure first proposed by
Jasak et al. (2014). It was originally proposed to suppress acceleration of air in the vicinity of free surface but we found that the limiter is also effective to the flows over the corners and sharp edges. As they proposed in their original work, we also adopted 5 x u0 as an upper limit of velocity and tested that it works well without losing accuracy.
3.4. Remarks on local time stepping
There are two major control parameters for the LTS scheme: smoothing and damping coefficients, which have values between 0 and 1. A field of inverse of local time step, so called as reciprocal time step in the OpenFOAM®, is estimated as following.
(D)max'2rCoDV
where Co is Courant number. It is further limited based on the interface Courant number near the free surfaces. The reciprocal time step is smoothed spatially if a local time step size of one cell is larger than 1 + Csmooth times of its neighbor cell's value, where Csmooth is a smoothing coefficient. Once again, it is further limited if the damping coefficient Cdamp is less than unity. If a local time step size increases too much, then the reciprocal time step is damped to be
(¿) -(1 ~ ^(D)
Fig. 9. Comparison of computational results with different grid arrangements
near corner.
The effects of smoothing and damping coefficients were tested by carrying out KCS simulations for 24 knots. Both coefficients were varied from 0.2 to 0.8 and one coefficient was set to be unity while the other had been tested. Table 1 shows the results of test on the smoothing coefficient and the force distribution with iteration steps are shown in Fig. 10. It is worthwhile to note that the computations were diverged if damping coefficient was solely applied without smoothing coefficient, thus the results of test on damping coefficients are not presented here. The computational efficiency is seemed to be saturated if smoothing coefficient is larger than 0.4, and it becomes worse if the smoothing coefficient increases to be 0.8. Despite the error tends to grow as the smoothing coefficient increases, all results of given coefficients still lie within the tolerant error bound. Therefore, it would be eligible to select the smoothing coefficient value between 0.4 and 0.8. In the rest of this paper, we selected 0.6 as the smoothing coefficient of simulations while a damping coefficient remains the unity.
4. Computational results
4.1. KRISO Container Ship (KCS)
We first validated the present method by running simulations on KCS with various speeds and the results were compared to those of STAR-CCM+ and experimental data.
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
Table 1
Computational results with various smoothing coefficient values.
Smoothing coefficients
Iterations
Elapsed time, s
Value x 1000
Error, %
Attitudes
Trim, deg
Sinkage, x1e2m
Exp. 0.2 0.4 0.6 0.8
9800 8200 8200 9400
6014 4683 4540 5224
n.681 n.657 n.646
0.73 0.81 1.45 1.76
0.169 0.168 0.166 0.169 0.169
— 1.394
— 1.096
— 1.085
— 1.085 1.078
Iterations
Fig. 10. Effect of smoothing coefficients on resistance simulations.
The principal particulars of the objective ship and simulation conditions are shown in Tables 2 and 3, respectively.
Fig. 11 shows computational domain, which has a hex-ahedral shape with its dimension of —1.5 < x/LPP < 3, 0 < y/LPP < 1.5, and —1.5 < z/LPP < 1 in x, y, and z directions, respectively. The fluid flows along x-direction and positive y-axis points the starboard direction while z-axis points upward. According to the present coordinate system, a displacement in z-direction represents a sinkage distance and a trim angle is defined as a rotational angle about y-axis. As boundary conditions for velocity and turbulent properties, Dirichlet condition was specified on inlet patch while Neumann condition was applied on outlet patch. For the case of pressure field, the fixed value boundary condition was used for
Table 2
Principal particulars of KCS.
Particulars Full scale Model scale
Length between perpendiculars Lpp [m] 230.0 7.2786
Length of waterline lwl [m] 232.5 7.3577
Maximum beam of water line BWL [m] 32.2 1.0190
Depth D [m] 19.0 0.6013
Draft T [m] 10.8 0.3418
Displacement volume V [mn] 52,030 1.6490
Wetted surface area without rudder Swet [m2] 9424 9.4379
LCB (%LPP), fwd+ — 1.48 — 1.48
outlet patch while the zero gradient boundary condition was applied to the rest patches of the computational domain. To take the effect of motion on phase fraction field into account appropriately, the modified boundary conditions for the phase fraction developed by Lee (2013) was adopted on inlet and side patches. The mesh was generated by means of trimmed mesh method of STAR-CCM+, which resulted in about 1.15 million cells.
LTS scheme for temporal discretization was applied and second order schemes for spatial discretization. The initial flow field was estimated for the first 5000 iterations under fixed motion condition. If the flow field converged or iteration number exceeded 5000, then the motion of a ship came into play. In the present simulations, the motion was determined based on averaged force and moment for the last 100 time steps where the motion was fixed during the period to estimate the force and moment on a ship in quasi-steady state. The computation was finalized when the motion is converged. The number of sub-cycle for phase fraction was set to be three with two PISO corrector loops. As a turbulent model, density weighted k-omega SST model was used which was modified based on the original k-omega SST model in OpenFOAM® to reflect the effect of density change in turbulent field near free surface. The density weighted turbulence models have been included in the OpenFOAM-3.0.x or later (OpenFOAM).
We also carried out simulations using STAR-CCM+ v10.04 with the Equilibrium DFBI model for the body
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15 11
Table 3
Computational conditions for KCS simulations.
Flow properties
Density r [kg/m3] 999.5
Kinematic viscosity V [m2/s] 1.27 x 10—6
Speed conditions
Case No. 1 2 3 4 5 6
Speeds VM, [m/s] 0.915 1.281 1.647 1.922 2.196 2.379
Froude Number, Fn 0.108 0.152 0.195 0.227 0.260 0.282
Reynolds Number, Re 5.23 x 106 7.33 x 106 9.42 x 106 1.10 x 107 1.26 x 106 1.36 x 106
motion (STAR-CCM+ User Guide), which accelerates prediction in the equilibrium position of a ship motion. Both computations using HiFoam and STAR-CCM+ have been carried out using 64 cores of HPC cluster of HHI, which consists of Intel Xeon CPU's working at 2.60 GHz. Please note that we did not carry out the grid test since the efficiency and robustness of HiFoam are the main concern and the accuracy has been discussed a lot in the earlier studies in detail (Park et al., 2013; Lee, 2013, 2014).
Table 4 shows the estimated resistance and elapsed time of KCS simulations as well as the experimental data. Distributions of resistances, trim angles, and sinkage distances are shown in Fig. 12. As we can see, the present method predicts a resistance of a ship mostly in an hour on 64-core cluster while it maintains accuracy of about ±1% for KCS cases. Consequently, the simulations of six different Fn conditions ended within five and a half hours without losing accuracy or causing any diverged case. Fig. 9 also compares the various numerical results of Tokyo 2015 Workshop on CFD in Ship Hydrodynamics (Larsson et al., 2015) where a diamond symbol represents the averaged CT value of participants and error bars denote the range of the data. The present method predicts the trim angle very well compared to other numerical methods,
Fig. 11. Computational domain for KCS simulations.
but it tends to under-predict the sinkage distance in the highspeed regime while other numerical simulations provide good predictions.
4.2. Japan bulk carrier
Japan Bulk Carrier (JBC), a capsize bulk carrier designed by National Maritime Research Institute (NMRI), Yokohama National University and Ship Building Research Centre of Japan (SRC), was introduced at Tokyo 2015 A Workshop on CFD in Ship Hydrodynamics. It has a stern duct as an Energy Saving Device (ESD). Towing tank experiments with and without ESD were carried out and the resistance and self-propulsion data is available as well as PIV measurement near stern. We selected the JBC as another objective ship for validating the capability of HiFoam to apply to blunt ships with appendages. Table 5 shows the principal particulars of JBC where the LPP of model scale is 7.0 m.
Unlike KCS, only one speed condition is available for JBC where the Froude number is 0.142. The density and viscosity are p = 998.2 kg/m3 and n = 1.107 x 10~6 m2/s, respectively. The computational domain was set to be identical to the cases of KCS simulations. Fig. 13 shows grid arrangements on z = 0 plane and near stern with and without the ESD. Please refer to the earlier study (Kim and Jun, 2015) for more detailed information on grid test of JBC simulations using OpenFOAM®. The mesh was refined on free surface plane to capture Kelvin waves and near stern to resolve the wake structure appropriately. Consequently, resultant meshes for without and with ESD cases have 1.43 M and 1.52 M
Table 4
Comparison of computational results to experimental data.
Fn CT x 1000 Elapsed time [s]
EXP CCM+ HiFoam CCM+ HiFoam
Value Error Value Error
0.108 3.796 3.790 —0.16% 3.810 0.36% 2430 1507
0.152 3.641 3.623 —0.48% 3.653 0.32% 2451 1649
0.195 3.475 3.460 —0.44% 3.504 0.82% 2469 3107
0.227 3.467 3.491 0.69% 3.486 0.55% 2435 3779
0.260 3.711 3.659 — 1.39% 3.724 0.34% 2431 3562
0.282 4.501 4.510 0.19% 4.483 —0.41% 2438 2296
+ MODEL
12 G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
Froude Number Fn
(a) Total resistance coefficients
Froude Number Fn
( b ) Trim angle
(c) Sinkage distance
Fig. 12. Comparison of numerical simulations to experimental data.
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15 13
Table 5
Principal particulars of JBC.
Particulars Full scale
Length between perpendiculars lpp [m] 280.0
Length of waterline lwl [m] 285.0
Maximum beam of water line BWL [m] 45.0
Depth D [m] 25.0
Draft T [m] 16.5
Displacement volume V [m3] 178,369.9
Wetted surface area without ESD Swet [m2] 19,556.1
Wetted surface area with ESD 19,633.9
LCB (%LPP), fwd+ 2.5475
(c) Mesh near stem with ESD
Fig. 13. Grid arrangement of JBC simulations with and without ESD.
cells, respectively, which have about 24% more cells relative to KCS simulation cases. The density weighted k-omega SST model was applied to predict the turbulent properties, and Gauss linear scheme was adopted for the gradient operations.
Table 6 compares the computational results to experimental data. We see that results of present simulations agree very well to the experimental data, despite an additional appendage exists. Furthermore, the present method properly reflects the
effect of ESD, which helps to slightly reduce total resistance acting on the objective ship. In the present work, the simulations took about 2206 and 2657 s for without and with ESD cases, respectively. Especially the present simulation predicts trim angle very accurately to the measured data while sinkage distances have about 1.5% of errors relative to the experimental data. Therefore, we can conclude that the present method predicts the resistance of a blunt ship with additional appendage accurately and efficiently.
Fig. 14 compares the resultant axial velocity contours of without and with ESD cases to PIV measurement data at three different positions along x-direction near stern. The present method tends to under-predict the vertical flows near stern, which may due to the turbulence model. As it was discussed at the CFD Workshop in Tokyo 2015, Reynolds averaged turbulence models such as RSM and two-equation models tend to under-estimate the bilge vortex while LES overly estimates the vertical flow. However, the present method predicts important flow features well, especially disturbed flow field due to ESD at x/Lpp = 0.9843 in Fig. 14.
4.3. Validations on various ships
The main purpose of the present work is to develop an efficient and robust CFD tool for evaluating resistance of a ship. The present method was applied to four different types of ships with various speed conditions. We selected 4 container carriers, 4 VLCC's, 3 LNG/LPG carriers, and 2 PCTC's that have been tested and manufactured by HHI as validation cases. The Froude numbers of test cases lie within the range of 0.092 < Fn < 0.261 where total 82 cases have been tested. Fig. 15 compares the resultant total resistances of computational results to experimental data. Please note that the values are non-dimensionalized by dividing maximum resistance value of each type to encrypt data due to security regulations. The solid line indicates the experimental data and dashed lines represent the ±2% error boundaries. The results of present simulations are well aligned along the experimental data line and lie within the ±2% error lines while it occasionally predicts resistances with error of larger than this error boundaries. The meshes for the present numerical simulations consisted of 1.46 million cells in average where twin skeg LNG cases required up to 2.0 million cells. In terms of computational speed, it took about 5800 s in average for the computations without causing any diverged case.
Table 6
Computational results of JBC simulations without and with ESD.
Without ESD With ESD
Exp HiFoam Exp HiFoam
Value Error Value Error
CTM, X1000 Trim, %Lpp Sinkage, %LPP 4.289 -0.180 -0.086 4.311 -0.181 -0.085 0.52% 0.33% -1.55% 4.263 -0.180 -0.086 4.300 -0.181 -0.085 0.86% 0.39% -1.56%
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15
Fig. 14. Comparison of Axial velocity contours of PIV measurements and computational results.
Fig. 15. Estimated dimensionless total resistances of various types of ships.
5. Concluding remarks
HiFoam, a numerical analysis tool based on the Open-FOAM®, has been developed to evaluate performances of a ship numerically. The interface capturing scheme using VOF method was applied to reflect the effect of free surfaces and the phase fraction equation was solved explicitly. To extend its applicability to problems in industrial field, it is essential to ensure efficiency and robustness of the solver. A simple flux limiter was proposed in the vicinity of free surface to enhance the solver's robustness where a sudden change of free surface may cause an instability in VOF simulations. The proposed limiter was verified by applying it to a dam break problem and the limiter successfully distinguished one phase region from another. We also proposed QS-FBI mesh motion with least square method to accelerate prediction of ship motions and showed that it saved about 38% of computational time in predicting attitude of KCS. As minor factors, the local time stepping scheme was reviewed and corner flows as well to find the most appropriate set up for numerical simulations in term
+ MODEL
G.-H. Kim, S. Park / International Journal of Naval Architecture and Ocean Engineering xx (2017) 1—15 15
of efficiency and robustness. The present method was applied to KCS simulations to verify its feasibility and we could conclude that it was comparable to the commercial codes. The HiFoam was also applied to JBC simulations without and with ESD cases to see its applicability to a ship of low speed, blunt body with an appendage. As a result, the present method provided very close solution in both total resistance and attitudes to the experimental data. Finally, it has been applied to various types and speeds of ships including container carriers, VLCC's, LNG/LPGC's, and PCTC's. We see that the HiFoam is highly efficient and robust with good accuracy in predicting resistance of ships, which is competitive with other commercial CFD codes and fully capable to replace them.
References
Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for
modeling surface tension. J. Comput. Phys. 100, 335—354. Jasak, H., Vukcevic, V., Christ, D., 20—27 November 2014. Rapid free surface simulation for steady-state hull resistance with FVM using FOAM-extend. In: 30th Symposium on Naval Hydrodynamics, Hobart, Tasmania, Australia.
Kim, G.H., Jun, J.H., 2—4 Dec 2015. Numerical simulations for predicting resistance and self-propulsion performances of JBC using OpenFOAM and STAR-CCM+. In: Tokyo 2015 a Workshop on CFD in Ship Hydrodynamics. National Maritime Research Institute, Tokyo, Japan. Larsson, L., Stern, F., Visonneau, M., Nobuyuki, H., Takanori, H., Kim, J., 2—4 Dec 2015. In: Proceedings of Tokyo 2015 a Workshop on CFD in Ship Hydrodynamics. National Maritime Research Institute, Tokyo, Japan. Lee, S.B., 23 —26 June, 2014. Application of OpenFOAM to prediction of ship performances. In: 9th International OpenFOAM Workshop, Zagreb, Croatia.
Lee, S.B., 11—14 June, 2013. Application of OpenFOAM to prediction of hull resistance. In: 8th International OpenFOAM Workshop, Jeju, South Korea. Muzaferija, S., Peric, M., Sames, P., Schelin, T., 1999. A two-fluid Navier-Stokes solvers to simulate water entry. In: Twenty-Second Symposium on Naval Hydrodynamics. Newman, J.N., 14—18 Dec 1992. Panel methods in marine hydrodynamics. In: 11th Australasian Fluid Mechanics Conference, University of Tasmania, Hobart, Australia.
OpenFOAM, http://www.openfoam.org.
Palacios, F., Colonno, M.R., Aranake, A.C., Campos, A., Copeland, S.R., Economon, T.D., Lonkar, A.K., Lukaczyk, T.W., Taylor, T.W.R., Alonso, J.J., 2013. Stanford University Unstructured (SU2): An Open-Source Integrated Computational Environment for Multi-Physics Simulation and Design. AIAA Paper 2013-0287.
Palacios, F., Economon, T.D., Aranake, A.C., Copeland, S.R., Lonkar, A.K., Lukaczyk, T.W., Manosalvas, D.E., Naik, K.R., Padron, A.S., Tracey, B., Variya, A., Alonso, J.J., 2014. Stanford University Unstructured (SU2): Open-Source Analysis and Design Technology for Turbulent Flows. AIAA Paper 2014-0243.
Park, S.H., Rhee, S.H., 2014. Numerical diffusion decrease of free-surface flow analysis using source term in volume fraction transport equation. J. Comput. Fluids Eng. 19 (No.1), 15—20.
Park, S.H., Park, S.W., Rhee, S.H., Lee, S.B., Choi, J.E., Kang, S.H., 2013. Investigation on the wall function implementation for the prediction of ship resistance. Int. J. Nav. Archit. Ocean Eng. 5 (No.1), 33—46.
Rhie, C.M., Chow, W.L., 1982. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation. In: AIAA-82-0998, AIAA/ ASME 3rd Joint Thermophysics, Fluids, Plasma and Heat Transfer Conference, St. Louis, Missouri.
Rusche, H., 2002. Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions. Ph.D. thesis. Imperial College of Science, Technology & Medicine.
Sweby, P.K., 1984. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J. Numer. Analysis 21 (No.5), 995—1011.
Ubink, O., 1997. Numerical Prediction of Two Fluid Systems with Sharp Interface. Ph.D. thesis. University of London.
STAR-CCM+ User Guide, CD-adapco.
Van, S.H., Kim, W.J., Yim, G.T., Kim, D.H., Lee, C.J., 1998. Experimental investigation of the flow characteristics around practical hull forms. In: Proceedings 3rd Osaka Colloquium on Advanced CFD Applications to Ship Flow and Hull Form Design, Osaka Japan.
van Leer, B., 1979. Towards the ultimate conservative difference scheme. J. Comput. Phys. 32 (1), 101—136.
Weller, H.G., 2002. Derivation, Modeling, and Solution of the Conditionally Averaged Two-phase Flow Equations. Technical Report TR/HGW/02. Nabla Ltd.
Witherden, F.D., Farrington, A.M., Vincent, P.E., 2014. An open source framework for solving advection-diffusion type problems on streaming architectures using the flux reconstruction approach. Comput. Phys. Commun. 185, 3028—3040.
