RANS prediction of the KVLCC2 tanker in head waves

: The present study is devoted to the computation of the KVLCC2 tanker in head wave with free heave and pitch motion. A RANS solver using finite-volume discretization and free-surface capturing approach is employed for the computation. Free ship motion is captured with a mesh deformation approach. Three different wave lengths (0.6Lpp, 1.1Lpp and 1.6Lpp) are computed. We focus on numerical uncertainty estimation in this paper. For each test case, three different meshes and at least three different time steps have been used to access both time and spatial discretization error. Additional computations with different setups aimed at identifying different numerical discretization errors will also be performed. It is demonstrated that special attention needs to be paid to time discretization. To keep the same time accuracy, time step needs to be reduced on fine mesh for such kind of unsteady free-surface computation involving important pitch or roll motion.


INTRODUCTION
With the progress made both in simulation software and computer hardware, CFD computation for steady ship hydrodynamic application can now be performed routinely with reasonable time and sufficient accuracy even at full scale. However, accurate prediction for applications involving strong interaction between ship motion and free-surface such as the stability study of high speed vehicle remains a challenging task for CFD computation. Another example for such application is the simulation of ship motion in waves, namely the so-called seakeeping simulation. Potential flow solver in frequency or time domain is commonly employed for this kind of problem. A review of different numerical approaches using potential flow model can be found in the report of ITTC seakeeping committee [1] . Although potential flow solver has been continuously improved for seakeeping simulation, extreme wave condition such as green water and slamming impact loading, deck load, whipping, breaking waves, dynamic structural loading, etc. can not be handled properly with potential flow assumption. In such situation, RANS simulation can provide a more reliable prediction. However, RANS simulation for seakeeping applications requires an accurate prediction both in time and in space. Due to excessive requirement in computational resources, validation of RANS computation for seakeeping applications remains hurdle. Weymouth et al. [2] have performed a validation work for a modified Wigley model and have shown a good agreement between RANS computation and measurement over a wide range of wave frequency and wave height. Another ship model frequently employed for RANS seakeeping validation is the DTMB 5415 frigate model. It has been chosen as one of the test case for the Tokyo 2005 CFD Workshop [3] . 4 groups including that of the authors have performed computation for the DTMB 5415 model with fixed position in waves. Although standard deviation for 0 th amplitude for total resistance and heave force are 4.5% and 6.7% respectively, that for the pitching moment is as high as 70%. Simulation for this test case is also reported by Carrica et al [4] where an additional test case with deeper waves is also presented. Unsteady breaking bow wave is observed in their computation for the case with higher wave steepness. The authors have performed further investigation for this model recently [5] and [6] with detailed comparison with available experiments on the forces coefficients, the ship motion history, and the phase-averaged free-surface elevations. Good agreement between numerical simulation and measurement data has been observed at both moderate (Fr=0.28) and high (Fr=0.41) Froude numbers. We have also performed a numerical simulation for a model container ship in head waves and in oblique waves [7] for different wave lengths. It has been found that at high wave frequency, added resistance in waves is under-estimated compared with the measurement data. We continue to investigate this issue in the present study. The test case chosen is the KVLCC2 tanker in head waves. This is one of the test cases chosen for the Gothenburg 2010 workshop. Measurement data are not yet available. We will focus mainly on numerical uncertainty investigation in the present paper. Three different wave lengths (0.6Lpp, 1.1Lpp and 1.6Lpp) have been simulated. Different mesh size and time step have been employed to access numerical discretization error both in time and in space.

FLOW SOLVER
Computation has been performed with the ISIS-CFD flow solver developed by EMN (Equipe Modélisation Numérique, i.e. CFD Department of the Fluid Mechanics Laboratory). Turbulent flow is simulated by solving the incompressible Reynolds-averaged Navier-Stokes equations (RANS). The flow solver is based on finite volume method to build the spatial discretization of the transport equations. The velocity field is obtained from the momentum conservation equations and the pressure field is extracted from the mass conservation constraint, or continuity equation, transformed into a pressure-equation. In the case of turbulent flows, additional transport equations for modeled variables are discretized and solved using the same principles. The gradients are computed with an approach based on Gauss's theorem. Non-orthogonal correction is applied to ensure a formal first order accuracy. Second order accurate result can be obtained on a nearly symmetric stencil. Inviscid flux is computed with a piecewise linear reconstruction associated with an upwinding stabilizing procedure which ensures a second order formal accuracy when flux limiter is not applied. Viscous flux are computed with a central difference scheme which guarantee a first order formal accuracy. We have to rely on mesh quality to obtain a second order discretization for the viscous term. Free-surface flow is simulated with a multi-phase flow approach. Incompressible and non-miscible flow phases are modeled through the use of conservation equations for each volume fraction of phase/fluid. Implicit scheme is applied for time discretization. Second order three-level time scheme is employed for time-accurate unsteady computation. Velocity-pressure coupling is handled with a SIMPLE like approach. Ship free motion can be simulated with a 6 DOF module. Some degree of freedom can be fixed as well. An analytical weighting mesh deformation approach is employed when free-body motion is simulated. Several turbulence models ranging from one-equation model to Reynolds stress transport model are implemented in ISIS-CFD. Most of the classical linear eddy-viscosity based closures like the Spalart-Allmaras one-equation model, the two-equation k-ω SST model by Menter [5] , for instance are implemented. Two more sophisticated turbulence closures are also implemented in the ISIS-CFD solver, an explicit algebraic stress model (EASM) [6] and a Reynolds stress transport model. The EASM model is employd in the present study. Wall function is implemented for two-equation turbulence model.

Description of test case
The test case chosen for the present study is the test case 1.4a proposed for the Gotherburg 2010 Workshop in Ship Hydrodynamics. The geometry is the well known KVLCC2 tanker. Measurement under similar condition with the same scale factor 100 are being performed at INSEAN and at Osaka University. The two measurements differ only by the fact that surge motion is free in the measurement performed at Osaka University. With a Froude number of 0.142 at service speed, ship speed at model scale is 0.797 m/s, giving a Reynolds number of 2.55 10 6 . In the present study, we do not perform the case with free surge motion. Only pitch and heave motions are free. Three different wave lengths (0.6Lpp, 1.1Lpp and 1.6Lpp) investigated in the measurement have been simulated. Wave height H=0.06m (0.288T) is the same for all cases. The depth of the original KVLCC2 model with a scale factor of 100 is 0.3m. With a draft T=0.208m, green water on deck is observed in the measurement for the case with long wave length. To avoid such undesirable effect, the deck of the ship model is extended to 0.35m in the measurement performed in INSEAN. The same modification on ship geometry is applied in numerical simulation. With this modification, monitoring of pressure loading on the deck does not show any green water on deck during the simulation for all test cases.
Since the sway, yaw and roll motion is fixed, only half boat is simulated by using a mirror plane. The computational domain is extended to about 1Lpp before the ship where wave boundary condition is applied. Up and down boundaries are located at about 0.5LPP and 1.5 LPP respectively where pressure boundary condition is applied. Side boundary is also treated as a mirror plane and is located at about 1.8Lpp. Downstream boundary is located at about the same distance after the boat where classical exit boundary condition is applied. Unstructured hexahedral mesh is generated with a commercial software Hexpress. Three different meshes are generated for each case to access spatial discretization error. Based on the parameter used for mesh generation, the refinement ratio from the coarse grid to the medium grid and from the coarse grid to the fine grid are about 1.33 and 1.67 respectively. The grid is refined in the longitudinal direction near the free surface such that there is about 60 nodes per wave length and 10 nodes in wave height on the coarse grid. Based on our experiences, with such spatial resolution, incident waves can be captured with sufficient accuracy.

Resistance prediction in calm water
In order to be able to evaluate added resistance in waves, it is necessary to compute resistance in calm water. Such computation can also provide an indication to spatial discretization error. Table 1 gives the predicted resistance in Newton obtained on three different grids with half boat. Number of cells for each grid is also given in the table. The last line gives the extrapolated value and the observed order of convergence based on the grid refinement ratio mentioned above. Based on the extrapolated value, numerical uncertainty is about 1%, 2% and 3% respectively for the fine grid, the medium grid and the coarse grid without taking into account any safety factor. This is a typical range uncertainty level for resistance prediction using the ISIS-CFD code with a recommended setup. The medium grid is the grid density commonly used for routine application.

Time discretization error estimation
Uncertainty estimation for resistance in waves prediction is much difficult. In addition to usual time and spatial discretization error of the Navier-Stokes equation, additional sources of error exist such as time discretization error for the resolution of free body motion, additional discretization induced by the treatment of deforming grid, sampling uncertainty etc. We first examine the time discretization error by performing computation with the same grid using different time steps. Data sampling is performed over the same time interval to reduce sampling uncertainty. Table 2 presents the 0 th and the 1 st harmonic amplitudes for the resistance in waves for the half boat simulation. The first column is the number of time step per period. The last line is the extrapolated value using a least square fit approach assuming a second order accuracy. The observed order of convergence is also given in the same line. The reason why the assumed order of convergence is used for the uncertainty estimation is that for some case, monotonic convergence behavior is not observed. Hence, estimation based on observed order of convergence is not possible. Alternative estimation procedure must be used, which makes it difficult to compare different cases. By using an assumed second order accuracy for uncertainty estimation, result obtained will not be always reliable for small time steps. However, we believe that it still conveys a trustful information for large time step most of the cases. By comparing the numerical prediction with the extrapolated value, one can obtain the numerical uncertainty for each case. This estimation is also indicated in percentage on each column. Again safety factor is not used for this estimation to make the comparison for different cases more straightforward.   Tables 3 and 4 show similar results obtained with the medium grid and with the fine grid respectively. It is observed that for a given time step, time discretization error depends strongly on the grid density. While 125 time steps per period gives quite acceptable prediction with about 3% uncertainty on the coarse grid, the same time step induces an error of about 10% on the medium grid and on the fine grid. To reduce the time discretization error to 3%, which can be considered as acceptable for engineering application, one needs to reduce the time step by twice. On the other hand, opposite trend is observed for the 1 st harmonic amplitude. With the same time step, time discretization error decreases when the grid is refined. 200 time steps per period can be considered small enough for 1 st harmonic amplitude. But for 0 th harmonic amplitude, at least 250 time steps per period is necessary.

Spatial discretization error estimation
Reliable uncertainty estimation for spatial discretization error for the three test cases investigated in the present study is difficult. For the case with wave length=1.6Lpp, the 0 th harmonic predicted with three different grids using the finest time step differ by about 2.9%, while the time discretization error is as high as 4.6% on the fine grid. It is not surprising that an oscillating convergence behavior is observed. Smaller time discretization error is observed for the 1 st harmonic amplitude. A monotonic convergence behavior is observed. Table 11 shows the estimated uncertainty based on the prediction obtained with each grid using the smallest time step. Although safety factor is not taken into account, the prediction obtained with the medium grid using the smallest time step can be considered as acceptable. Similar behavior is observed for the case with wave length=1.1Lpp (table 12). Due to high time discretization error (about 3%), oscillating convergence behavior is observed for the 0 th harmonic amplitude. However, data difference is relatively small (about 2.4%), indicating a good spatial resolution. Higher level of uncertainty is observed for the 1 st harmonic amplitude. But the level is still acceptable for all grids. For the case with short wave length=0.6Lpp shown in table 13, spatial discretization error is smaller than sampling error. Results are nearly grid independent both for the 0 th and for the 1 st harmonic amplitude. Numerical prediction can be considered as accurate enough, even with the coarse grid.

Additional investigations
In the present study, we determine the time step based on the encounter frequency. Results obtained indicate that the encounter frequency is not the only criterion that needs to be taken into account for seakeeping simulation. While 94 time steps per period can give a time step independent solution for the case with wave length=0.6Lpp, the time discretization error can be as high as 30% for the case with wave length=1.1Lpp on the fine grid. Additional computations have been performed in order to find out the sources of numerical error that lead to such high level of numerical uncertainty. Computations will mainly be performed with the coarse grid with parameters corresponding to the case with wave length=1.6Lpp.
We first perform a computation without free ship motion. This computation gives a uncertainty level of 3.3% and 1.3% respectively for the 0 th and 1 st harmonic amplitudes for the resistance with the coarse time step (94 time steps per period). The case with wave length=1.1Lpp gives a still lower level of uncertainty (0.1% and 1.7%). These observations suggest that incident waves as well as its interaction with a fixed boat can be simulated with sufficient accuracy using about 100 time steps per period.
The next computations are intended to determine numerical discretization error due to the use of deforming mesh. We first perform a mono-fluid computation with an imposed sinusoidal pitch motion with an amplitude (0.04 rad) and period (1.4127 second) corresponding to the case with wave length=1.6Lpp. This simulation gives respectively 1.6% and 0.6% of uncertainty for the 0 th and 1 st harmonic amplitudes for the resistance with the largest time step. Such low level of uncertainty suggests that the discretization error of the Navier-Stokes equation on deforming mesh is not responsible for the high level of uncertainty observed in the seakeeping computation under similar condition. Then we perform a similar computation with free surface. Incident wave is not taken into account to avoid possible phasing error. This computation gives 0.8% and 9.3% of uncertainty for the 0 th and 1 st harmonic amplitudes respectively for the resistance with the largest time step. This result reveals that high amplitude of pitch motion with free surface can induce important time discretization error due to the use of deforming mesh. Another computation with an amplitude of 0.032 rad and period of 1.1206 second corresponding to the case with wave length=1.1Lpp gives 0.3% and 3.9% of uncertainty for the 0 th and 1 st harmonic amplitudes respectively. This level of uncertainty is much lower than the case for seakeeping prediction (13.0% and 3.6%).
Another source of error that needs to be quantified is the time discretization error due to the resolution of free ship motion. To quantify this error, we perform a simulation with free surge motion in waves. Other degree of freedoms are frozen. With this kind of simulation, we can use a rigid mesh motion, which eliminates the numerical discretization error due to deforming mesh. For this computation, a constant thrust is applied to the ship. Hence, 0 th harmonic amplitude of the resistance can not be used to access numerical uncertainty. The uncertainty observed on the 1 st harmonic amplitude of the resistance is very small (1.3%). 0 th and 1 st harmonic amplitudes for the velocity is also very small (1.2% and 1.1% respectively). However, for the pitch moment, 5.2% of uncertainty is observed for largest time step with an observed order of convergence of 1.2 for the 0 th harmonic amplitude. Similar order of convergence is also observed for the 0 th harmonic amplitude for the velocity, while a second order time discretization is used. Such an unexpected behavior requires a further investigation. For the 1 st harmonic amplitude for the pitch moment, the uncertainty level for the largest step is small (1.7%).

CONCLUSIONS
The KVLCC2 tanker in head waves has been simulated for three different wave lengths. We focus on numerical uncertainty estimation in this paper. For each case, computation has been performed with three different meshes and at least three different time steps. Time discretization error can be clearly quantified except for the case with short wave length for which time discretization error becomes smaller than sampling error. For the case with high wave length for which ship motion amplitude is high, time discretization error becomes very important. In addition, it increases when fine grid is used. Two critical issues for numerical discretization have been identified. The first one is the high time discretization error induced by the association of strong pitch/roll motion, deforming mesh, and free surface. Such kind of time discretization error increases as the grid is refined. A new criterion for time step is required for such kind of simulation. This issue has not been investigated in the present study. The second issue that has an important impact on numerical accuracy is the time discretization error for the resolution of body motion. Although a second order time discretization scheme is employed, observed order of convergence is closer to first order than to second order for most cases. Associated time discretization error seems to be quite high. Further investigations on this issue are required in future studies. Due to relatively high time discretization error and sampling error, estimation for spatial discretization error for seakeeping prediction is difficult. Nevertheless, we can consider that a grid that can provide an accurate enough prediction for a steady flow can also provide an accurate enough spatial solution for unsteady flow computation. With the grid density commonly used for routine application, time discretization error is higher than spatial discretization error when less than 200 time steps per period are employed when ship motion is important. An unnecessary fine grid is not recommended for such kind of simulation to avoid high time discretization error.