Predictions of ship maneuverability based on virtual captive model tests

Maneuverability is an important hydrodynamic performance of a ship, and should be taken into account during the ship design stage. The present study of Computational Fluid Dynamic (CFD) calculations aims to offer a numerical tool for maneuvering prediction with high accuracy. The virtual captive model tests for a model scale KCS container ship are conducted using unsteady Reynolds-averaged Navier–Stokes (RANS) computation to obtain the full set of linear and nonlinear hydrodynamic derivatives in the 3rd-order Abkowitz model. The numerical uncertainty analysis is carried out for the pure sway and yaw–drift tests to verify the numerical accuracy. It is concluded that the lower order Fourier coefficients are preferred in the computation of the hydrodynamic derivatives. Moreover, part of the computed hydrodynamic forces and moments are compared with the available captive model test data, and good agreement is obtained. By substituting the computed hydrodynamic derivatives into the mathematical model, the standard turning and zigzag maneuvers are predicted. By comparing the predicted maneuvering results with the available experimental data and the prediction results by others, it is demonstrated that acceptable prediction accuracy can be achieved with the present method, which shows the effectiveness of the present method in predicting ship maneuverability.


Introduction
Maneuverability is an important hydrodynamic performance of a ship. The International Maritime Organization (IMO) has promulgated the Standards for Ship Manoeuvrability (IMO, 2002) to ensure the safety of navigation. To assess ship maneuverability at the design stage, reliable methods for predicting ship maneuverability are required. The Manoeuvring Committee of the International Towing Tank Conference (ITTC, 2008) and the Workshop on Verification and Validation of Ship Manoeuvring Simulation Methods (SIMMAN, 2008(SIMMAN, , 2014 collected and sorted all methods for prediction of ship maneuverability in practical applications. In general, maneuverability predictions using free running model tests (FRMT) such as zigzag and turning circle tests are believed to be most close to the reality. At the SIMMAN workshops, many organizations such as Maritime Research Institute Netherlands (MARIN) provided FRMT data for the standard ship models KVLCC1 & 2, KCS and DTMB 5415 (Quadvlieg, 2017). Besides, the system-based prediction method, which uses the mathematical model of maneuvering motion with the hydrodynamic coefficients obtained from empirical formulae (EMP) (Kijima, Katsuno, Nakiri, & Furukawa, 1990), CONTACT Lu Zou luzou@sjtu.edu.cn captive model tests, numerical computations, or system identifications (Sutulo & Guedes Soares, 2014;Tran Khanh, Ouahsine, Naceur, & El Wassifi, 2013;Zhang & Zou, 2011), is widely applied. Among the methods for determining the hydrodynamic coefficients, the captive model test method is regarded as the most reliable one, which includes an oblique towing test, rudder force test, planar motion mechanism (PMM) test, and circular motion test (CMT), etc. The advantage of the captive model test is that it can obtain all the linear and nonlinear hydrodynamic coefficients in the mathematical model. However, it also has some drawbacks, such as the fact that specific test facilities are required, and it is not convenient for the evaluation and optimization of ship maneuverability at the design stage. Nowadays, with the rapid development of computer techniques and the Computational Fluid Dynamics (CFD) method, CFD-based numerical methods have been successfully applied in the naval architecture (Zhang, Zhang, Tezdogan, Xu, & Lai, 2018), especially in the investigation of ship maneuverability (Stern, Yang, Wang, Sadat-Hosseini, & Mousaviraad, 2013;Sun, Su, Wang, & Hu, 2016). There are two types of CFD application in maneuvering predictions. One is the direct numerical simulation of standard maneuvers in the full time domain with steering rudder(s) and rotating propeller(s) (Dubbioso, Durante, Di Mascio, & Broglia, 2016;Wang, Zou, & Wan, 2017). However, the requirement of huge computational resources and the heavy computational burden limit its practical application. Another CFD application is more practical. It first determines the maneuvering hydrodynamic forces by conducting virtual captive model tests, i.e. simulating the captive model tests with a Reynolds-averaged Navier-Stokes (RANS)based solver, and then predicts the standard maneuvers using the mathematical models with the obtained hydrodynamic coefficients. Two classic types of mathematical models -the Abkowitz model (Abkowitz, 1964) and the Maneuvering Modeling Group (MMG) model (Yasukawa & Yoshimura, 2015) -for predicting ship maneuverability are generally adopted. More and more researchers are trying to use CFD methods for predicting the hydrodynamic derivatives. Otzen and Simonsen (2014) and Simonsen, Otzen, Klimt, Larsen, and Stern (2012), identified a reduced test matrix to be applied for PMM simulations, and standard maneuvers were conducted for the appended KCS ship model by using the Abkowitz model. Sakamoto, Carrica, and Stern (2012) conducted unsteady RANS (URANS) simulations of static and dynamic PMM tests for an un-appended surface combatant model 5415. Sung, Park, and Jun (2015) developed an easy procedure of virtual captive model tests to obtain the linear and nonlinear hydrodynamic derivatives for un-appended KCS and KVLCC1 & 2. He et al. (2016) computed the linear hydrodynamic derivatives of the KVLCC2 ship model by numerically simulating pure sway and pure yaw tests and performed the standard free running maneuvers using MMG mathematical model. Jin, Duffy, Chai, Chin, and Bose (2016) investigated the influence of scale effects on maneuvering coefficients for a KVLCC2 hull by using a URANS solver. Veedu and Krishnankutty (2016) investigated the effect of vessel speed on the hydrodynamic derivatives of the container ship S175. Shenoi, Krishnankutty, and Selvam (2016) numerically studied the maneuvering qualities of S175 using the derivatives obtained by a RANS-based solver, with the nonlinear and roll-coupled effects being taken into account. Liu, Zou, and Zou (2017) analyzed the effect of dynamic sinkage and trim on the hydrodynamic forces and hydrodynamic derivatives by numerically simulating static captive model tests. Guo and Zou (2017) conducted a system-based investigation on four degrees of freedom (4-DOF) ship maneuvering in calm water for the ONR tumblehome model by RANS simulation of captive model tests. A similar study for a destroyer hull was presented in Alexandersson, Korkmaz, and Mazza (2017). The virtual captive model tests have become a robust tool for acquiring the hydrodynamic derivatives/coefficients in the mathematical model of maneuvering motion at ship design phase. But few get the full set of the hydrodynamic derivatives/coefficients by this method. Most researchers used a modular model such as the MMG model in which the hull-rudder-propeller interaction coefficients were obtained by empirical formulae (He et al., 2016;Shenoi et al., 2016) or model tests (Sung et al., 2015). Others using the Abkowitz model simplified the mathematical model (Otzen & Simonsen, 2014;Simonsen et al., 2012). So far, no full set of hydrodynamic derivatives in the Abkowitz model has been presented for the KCS container ship. Most recently, Quadvlieg (2017) pointed out that the comparisons among numerical results from different institutes for the KCS ship under the same maneuvering conditions were not satisfactory. There was too much scatter in the predicted maneuvering parameters. To verify the reliability of CFD-based tools in maneuvering predictions better, more numerical simulations are needed to expand the database. Under this background, the present study aims to predict KCS ship's maneuverability using the Abkowitz model with the full set of linear and nonlinear hydrodynamic derivatives all determined by virtual captive model tests. The Abkowitz model is preferred since the effects of the hull-rudder-propeller interactions are automatically included and there is no need to determine the hull-rudder-propeller interaction coefficients. The paper is organized as follows. Firstly, grid and timestep dependency studies are carried out for a pure sway case and a yaw-drift case to estimate the numerical error and uncertainty due to grid discretization and time step. Then systematic computations are performed, and the predicted hydrodynamic forces and moments are presented and discussed. Some of the predicted forces and moments are compared with the available experimental data for validation purposes. Finally, by using the hydrodynamic derivatives determined from the virtual captive model tests, simulations of standard maneuvers are carried out and the results are compared with FRMT and CFD-based prediction results in the literature to validate the present numerical method and the established mathematical model.

Coordinate systems
Two right-handed coordinate systems are adopted, as shown in Figure 1. The body-fixed coordinate system o-xyz is used for computing the hydrodynamic forces. Its x-axis is pointing to the bow, y-axis pointing to starboard. The origin locates at the intersection of the water-line plane and the center-line plane at the mid-ship section. It is also the captive point in the captive model tests. The earth-fixed coordinate system o 0 -x 0 y 0 z 0 is set up to define the ship motion. The o 0 -x 0 y 0 plane is fixed on the undisturbed free surface and the z 0 -axis is pointing vertically downwards. u, v, and r are the surge velocity, sway velocity, and yaw rate in the body-fixed coordinate system; ψ is the heading angle (ψ = r). X and Y are the external force components in the longitudinal and transverse directions, respectively; N is the external moment about the ship's vertical axis. β is the drift angle defined by β = tan -1 (-v/u); δ is the rudder angle defined as positive when it turns to the port. U is the ship speed defined as U = √ u 2 + v 2 .

The Abkowitz model
The Abkowitz model for maneuvering motion in 3-DOF, i.e. surge, sway, and yaw is: where m is the ship mass, x G is the longitudinal position of ship's gravity center in the body-fixed coordinate system. I z is the moment of inertia about the z-axis.u,v anḋ r are the corresponding surge acceleration, sway acceleration, and yaw acceleration. Strøm-Tejsen and Chislett (1966) provides a full set of the expansion equations for the hydrodynamic forces and moment with detailed hydrodynamic derivations: X = X * + Xuu + X u u + X uu u 2 + X uuu u 3 + X vv v 2 + X vvu v 2 u + X rr r 2 + X rru r 2 u + X vr vr + X vru vr u + X δδ δ 2 + X δδu δ 2 u + X vδ vδ where u = u − u 0 is the disturbance in surge velocity. u 0 is the surge velocity in the initial state of forward motion. X u , X uu , X uuu , etc. are called hydrodynamic derivatives.

Captive model tests
There are two modes of captive model tests, 'static' and 'dynamic'. In the static mode, the model is constrained to move with a constant captive speed. Typical static captive model tests include a static drift test with a given drift angle β, a static rudder test with a given rudder angle δ, and a rudder-drift combined test with a given drift angle β and a given rudder angle δ. A static drift test is used to determine the derivatives relative to Similarly, the hydrodynamic derivatives relative to δ (X δδ ,Y δ ,Y δδδ ,N δ ,N δδδ ) and the cross-coupled derivatives relative to v and δ (X vδ ,Y vδδ ,Y vvδ ,N vδδ ,N vvδ ) can be determined from the static rudder test and the rudder-drift combined test, respectively. In these static tests, the captive speed is U c = U.
In the dynamic mode, typically the PMM tests, the model is forced to oscillate in a sinusoidal way so that various combinations of sway and yaw motions are achieved, including pure sway test, pure yaw test, yaw-drift, and yaw-rudder combined tests. In the present study, the pure sway test is used to determine the acceleration derivatives (Yv,Nv) with v oscillating harmonically (v = −v max cos ωt,v =v max sin ωt where the subscripts 'max' denote the amplitude, ω is the angular frequency), while the pure yaw test is used to determine the derivatives relative to the yaw rate, r (X rr ,Y r ,Y rrr ,Yṙ,N r ,N rrr ,Nṙ) with harmonic oscillating yaw motion (ψ = −ψ max cos ωt, r = r max sin ωt,ṙ =ṙ max cos ωt). The yaw-drift combined test is used to determine the cross-coupled derivatives (X vr ,Y vrr ,Y vvr ,N vrr ,N vvr ) with r oscillating harmonically and a prescribed drift angle β superposed (ψ = −ψ max cos ωt + β, v = −U c sin β), while the yaw-rudder combined test is used to determine the cross-coupled derivatives (X δr ,Y δrr ,Y δδr ,N δrr ,N δδr ) with r oscillating harmonically (ψ = −ψ max cos ωt, r = r max sin ωt,ṙ =ṙ max cos ωt) and a prescribed rudder angle δ superposed. In the dynamic tests, U ≈ U c for ψ max = y max ω/U c < < 1.

Computation of the hydrodynamic derivatives
For the static captive model tests, Equations (2)-(4) can be simplified as follows: Static drift: Static rudder: Rudder-drift: The hydrodynamic derivatives in Equations (5)-(7) can be obtained as polynomial coefficients by the leastsquare curve fitting method when the hydrodynamic forces and moments X, Y, N on the ship model have been acquired.
A pure sway test is mainly conducted to determine the linear acceleration derivatives (Yv,Nv) by imposing dynamically varying v andv. For the pure sway test, Equations (2)-(4) can be simplified as: For other dynamic captive model tests, the simplified mathematical models similar to Equation (8) can be obtained.
For the dynamic captive model tests, the Fourier series (FS) method (Sakamoto et al., 2012) is adopted for deriving the hydrodynamic derivatives from the hydrodynamic forces and moment X, Y, N. Since the sway and yaw motions are prescribed by sine and cosine functions, the forces and moment can be reconstructed as a Fourier series with the angular frequency ω of sway/yaw motion: where f (here X, Y, N) represents the time series of forces and moment; f cn and f sn are the nth-order Fourier sine and cosine coefficients, respectively.
On the other hand, for a pure sway test, the harmonic forms determined by substituting the model motion equation (v = −v max cos ωt,v =v max sin ωt) into the simplified mathematical model Equation (8) are: By using the same method, expressions similar to Equation (10) for other dynamic captive model tests can also be obtained. The Fourier sine and cosine coefficients in these expressions are given in Table 1. The functions of the N component are not listed in the table since they are the same as those of the Y component, and can be obtained by replacing 'Y' with 'N'. For the same reason, the functions in the yaw-rudder test are not listed in the table since they can be obtained by replacing 'v' in the functions of the yaw-drift test with 'δ'.
The 'Multiple-Run (MR)' approach (Yoon et al., 2015) is used to evaluate the hydrodynamic derivatives from these FS coefficients. This method uses curve-fittings of the mathematical models for sine and cosine coefficients similarly with the method for evaluating the derivatives of static captive tests mentioned above. The dynamic captive model tests are repeated over a range of test parameters of interest (e.g. v max for pure sway or r max for pure yaw test) to construct polynomial functions that  Table 1) for pure yaw motion is curve-fitted to the 1st-order Fourier sine coefficient Y s1 obtained from a series of pure yaw (varied r max ) tests to determine Y r and Y rrr . While all derivatives are available from the 0th-or the 1st-order coefficients, a few derivatives are also available from the higher order (the 2nd-or the 3rd-order) coefficients. The former is named as the MRL method (marked with an asterisk in Table 1) and the latter as the MRH method (Yoon et al., 2015).
The surge derivatives such as X u ,X uu , and X uuu in Equations (2)-(4) are determined by repeating the straight ahead case for a range of the captive speed U c with the constant propeller revolutions (Strøm-Tejsen & Chislett, 1966). The longitudinal force X, which represents the difference between the propeller thrust and the ship resistance, will vary as a function of speed u. From Equation (2) it follows that: The hydrodynamic derivatives X u , X uu , X uuu can be obtained directly as the coefficients of the 3rd-order polynomial. Derivatives such as X vv ,X rr ,  (1) is numerically solved using Euler's difference method (Obreja, Nabergoj, Crudu, & Păcuraru-Popoiu, 2010) to simulate the standard maneuvers. It should be noted that all governing equations are made non-dimensional by the speed U c , ship length L pp , and water density ρ. This provides the following nondimensional forms: Since all the variables in Equation (12) are nondimensional, the resultant hydrodynamic derivatives are non-dimensional as well.

Numerical method
To numerically simulate the captive model tests, the CFD software STAR-CCM+ (CD-Adapco, 2014) is applied. The code solves an incompressible unsteady RANS equation, closed by modeling the Reynolds stress tensor using the Realizable k-ε turbulence model (Shih, Liou, Shabbir, Yang, & Zhu, 1995). The reason to choose the Realizable k-ε model is that this model has been extensively used for predicting the hydrodynamic forces on a maneuvering ship. As noted by Quérard, Temarel, and Turnock (2008), the k-ε model for prediction of the hydrodynamic coefficients of ships is quite economical in terms of CPU time, compared to, for example, the SST turbulence model, which increases the required CPU time by nearly 25%.
The VOF method (Hirt & Nichols, 1981) is applied to capture the position of the free surface. A SIM-PLE algorithm (Patankar & Spalding, 1972) is employed for pressure-velocity coupling. A finite volume method (FVM) is used to discretize the flow domain into a finite number of control volumes (CVs) corresponding to computational grid cells. The temporal discretization is based on a 1st-order Euler difference, and the spatial discretization is performed with a 2nd-order upwind scheme for the convection term and secondary gradient contribution for the diffusion term. Mean flow quantities near the wall are simulated according to an all-y+ wall treatment where blended wall function is adopted (CD-Adapco, 2014). The all-y+ wall treatment is a hybrid treatment that attempts to emulate the high-y+ wall treatment for coarse meshes, and the low-y+ wall treatment for fine meshes. It is also formulated with the desirable characteristic of producing reasonable solutions when the wall-cell centroid falls within the buffer region of the boundary layer. Particularly, ship speed and drift angle are specified with the Dynamic Fluid Body Interaction (DFBI) Translation and Rotation module. It involves actual displacement of mesh vertices, which can model the 6-DOF motion of rigid body within the fluid system. Then the resultant force and moment acting on the body due to all influences are calculated, and the flow field is updated to find the new position and orientation of the ship. More details about the DFBI formulation can be found in Ohmori (1998).
The studies by Yao (2015) and Zipfel, Wagner, and Abdel-Maksoud (2014) showed that comparing the full propeller model by RANS simulation with a propeller body force model based on the actuator disk concept (Stern, Kim, Patel, & Chen, 1988), a fair overall agreement for the predicted total force can be achieved. Therefore, in the present study the rotating propeller is modeled as the body force radially distributed in the propeller disk with axial and tangential components. The body force is specified in a non-iterative manner in which the ship wake on the body force is neglected. The lateral force due to the propeller is ignored in the present study under the assumption that the propeller lateral force is a very small part of the total transverse force.

Ship model
The ship form studied in this paper is the KCS container ship model with bulbous bow and transom stern, appended with a suspended rudder and a propeller. The designed approaching speed is 24 knots in full scale. It is one of the benchmark ships used at SIMMAN for CFD validation analysis. The 1/52.667 scaled ship model with a length between perpendiculars of L pp = 4.3671 m and draught d = 0.20506 m is chosen. The particular dimensions are listed in Table 2. More details of the geometrical properties can be found in the SIMMAN (2014) website.

Computation cases
In order to obtain the full set of linear and nonlinear hydrodynamic derivatives in the 3rd-order Abkowitz model (Equations (2)

Computational domain and boundary conditions
A general view of the computational domain and boundary condition is depicted in Figure 2. The domain is set as 2L pp in front of the ship, 2.0L pp behind the ship, 1.2L pp above the undisturbed free surface, 2.5L pp in the transverse direction and below the undisturbed free surface, respectively (ITTC, 2011). The inlet, top, bottom, and side boundaries are treated as velocity inlets. The hydrostatic pressure is specified at the outlet plane. On the hull and rudder surfaces, the non-slip condition is specified. The pressure resistance fluctuation is usually found to be caused by the wave reflection at the non-physical side boundaries. Therefore, a numerical damping at the side boundaries is applied to remove the fluctuation.

Mesh generation
The grid structure around the ship is shown in Figure 3. The unstructured predominant hexahedral mesh is applied to discretize the computational domain. Near the ship surface, orthogonal prismatic cells are generated to improve the accuracy of the flow resolution. Since the blended wall function is applied, prism layers are used to achieve y+ values around 30. The free surface, the bow and stern of the hull, and the rudder parts are refined by using volume refinement boxes in grid generation in order to resolve the flow. The virtual disk of the propeller is refined to ensure the proper cell size in the areas where it is needed.

Numerical uncertainty analysis
Before proceeding to systematic computation, it is necessary to perform a verification study. In the present study, the verification study is conducted following the Grid Convergence Index (GCI) method developed by Roache (1998) and described in Celik, Ghia, Roache, and Freitas (2008). This method is currently used and recommended by the American Institute of Aeronautics and Astronautics (AIAA). The discretization errors caused by grid size and time step are evaluated for the two dynamic cases, one pure sway and one yaw-drift, as indicated in Table 4, assuming that the numerical uncertainties for the other cases are of the same order.
All mesh quantities are given as percentages in terms of a base size, in order to refine the grid as systematically as possible. The grid refinement is achieved by applying a refinement factor r G = √ 2 to the base size. Three grid sets (coarse, medium, fine) are adopted in the study. The fine, the medium, and the coarse grids consist of approximately 2.81M, 1.57M and 0.6M cells, respectively.  The refinement ratio of the time step is the same as refined mesh, that is r t = √ 2. Three sets of time steps are T/566, T/400, and T/283, respectively, where T = 1/f is the period of the oscillating motion. The mesh convergence analysis is performed with a fine time step, and the time-step convergence study is executed with a fine grid. The total computing time for coarse, medium, and fine grid density or time step is about 12, 24, and 36 hours, respectively.
The changes in solutions between two successive parameters are defined as: where φ k1 , φ k2 , and φ k3 correspond to the solutions with fine, medium, and coarse input parameters. The subscript k refers to the input parameters (i.e. 'G' for grid size and 't' for time step). The apparent order p k of the method is expressed by: The extrapolated values φ 21 ext can be calculated by: The approximate relative error between medium and fine solutions e 21 a and the extrapolated relative error e 21 ext can be computed as follows: The fine-grid or fine-time convergence index GCI 21 fine which represents the numerical uncertainty is calculated by: Traditionally, the verification of the static simulations is focused on constant forces or moments, while the dynamic simulations cover verification of time series for forces and moments on the time level such as in Simonsen and Stern (2008). However, from the ITTC (2014) Guidelines, it can be seen that instead of making the verification on the time series directly, it is recommended to approximate the time series of the forces and moments with FS coefficients, as is done in Sakamoto et al. (2012). Table 5 summarizes the results of the grid-spacing and time-step convergence study for only Y c1 , Y s1 , Y c3 , N c1 , N s1 , and N c3 in the pure sway case, since the corresponding hydrodynamic derivatives are obtained from those FS coefficients. As can be seen from Table 5, the lower order FS coefficients (Y c1 , Y s1 , N c1 , and N s1 ) have much smaller approximate relative errors, extrapolated relative errors, and uncertainty values than those of higher order FS coefficients (Y c3 and N c3 ) for grid resolution. The numerical uncertainties in the fine-grid solution GCI 21 fine for those lower order FS coefficients are less than 5%, while for higher order ones they are more than 25%. For the time-step convergence, all of those six coefficients are fairly insensitive to the size of the time step, since the maximum GCI 21 fine reaches only 3.25%. Table 6 summarizes the results of the grid-spacing and time-step convergence study for some FS coefficients in the yaw-drift case. Like the pure sway case above, the lower order FS coefficients have much smaller uncertainty values than those of higher order ones for both grid and time resolutions. The maximum numerical uncertainty among X s1 , Y 0 , Y s1 , N 0 , and N s1 is less than 11% in the mesh convergence study, while it is less than 2% in the time-step convergence study. From the numerical uncertainty analysis, it can be concluded that the lower order FS coefficients have much less numerical uncertainties. Therefore, MRL methods with the 0th-or the 1st-order FS coefficients are preferred to compute the hydrodynamic derivatives in the present study.
From the grid and time dependency studies, it can be seen that the approximate relative errors between medium and fine solutions in MRL methods are small, while the computational time is about doubled. Thus a medium grid density is chosen for the subsequent simulations to maintain an affordable computation cost, which contains about 1.6M cells. Besides, the medium time step (T/400) is chosen during the numerical simulation of dynamic captive model tests.

Results and discussion
After the verification study, systematic simulations of the captive model tests as listed in Tables 3 and 4 are carried out, and the total computational time required is about 20 days.

Hydrodynamic forces and moments
The hydrodynamic forces and moments on the ship are the main concern for evaluating the hydrodynamic derivatives. The computation results are validated by comparing with the Experimental Fluid Dynamics (EFD) data from the FORCE Technology's standard PMM tests (Simonsen, 2014). Figures 4-6 show the computation results of hydrodynamic forces and moments of the static tests in comparison with the EFD data. In general, the overall trend shows that the computation results of the nondimensional hydrodynamic forces are in good agreement with the experimental data. From Figure 4(a), it can be seen that although the computed trend of longitudinal forces X' at the lowest speed (Fr = 0.156) for positive drift angle are different from the experimental data, it is almost the same for negative drift angle. The average relative error E (E = [(EFD-CFD)/EFD]×100%) for longitudinal forces X' at Fr = 0.156 is about 8.9%; the standard deviation about the relative error is about 4.9%, which is acceptable. In the static rudder case at larger positive rudder angle and lowest captive speed, some discrepancies between CFD and EFD can be seen for Y' and N' (Figure 5(b) and (c)). The discrepancies might be due to the body force model of propeller. The propeller inflow is strongly dependent on the wake field. The wake flow changes significantly at lower speeds and large drift angles due to complex vorticity, while the simplified propeller model does not really reflect the difference in propeller load distribution over the disk. Thus, the predicted distribution of the propeller slipstream velocity that hits the rudder might not be accurate, which will affect the rudder performance and further influence X', Y', and N'. The experimental uncertainty could also be one of the causes.

Static mode
For the static drift case, X' increases when the captive speed decreases. This is because the propeller RPMs are kept constant at service speed self-propulsion point, and the propeller delivers too much thrust at the lower speeds. The three curves of Y' are almost coincident with each other, the same as those of N' (see Figure 4(b) and (c)), indicating that the non-dimensional transverse force and yaw moment are not sensitive to captive speed in static drift case. The influence of drift angle on global loads is well captured even at larger drift angles outside the linear region. For the static rudder case, approximately, the computed Y' and N' increase linearly with the increase of rudder angle. Y' and N' become nonlinear at the rudder angle larger than ± 20° (Figure 5(b) and (c)). For the rudder-drift case, the computed X', Y', and N' at the lower captive speeds (Fr = 0.156, 0.201) are validated by comparing with the experimental data. It can be seen that the slopes of the Y' and N' curves are both increasing as the captive speed decreases. Figure 7 shows the pure sway results at different maximum sway velocity and captive speeds. It can be seen that at the same captive speed, although the trend of the forces is almost the same for different transverse velocities, the amplitude of hydrodynamic forces increases as the sway velocity increases. Overall agreement between CFD and EFD is satisfactory, except the peak values in yaw moments. The time histories of forces and moment of pure yaw cases at different captive speeds and yaw rates are depicted in Figure 8. The overall trend shows that the computation results agree well with the experimental data. The amplitude of Y' and N' increases as the yaw rate increases at each captive speed. Figure 9 shows the results of the yaw-drift case. Due to the drift angle, asymmetric transverse force and yaw moment can be observed. The changing trend of the hydrodynamic forces is almost the same for different drift angles, but the amplitude of the forces increases as drift angle increases. Some discrepancy can be seen for the longitudinal force coefficients X' compared with EFD (Figure 9(a) and (d)), which might be due to the body force model for the propeller. The noise in the experiment and the relative small values may also lead to error. Figure 10 shows the results of the yaw-rudder case. Like the yaw-drift case, asymmetric transverse force and yaw moment can be observed.

Simulation of the standard maneuvers
The hydrodynamic derivatives with respect to the rudder angle, sway velocity, and yaw rate are obtained based on the virtual captive model tests. The results are listed in Table 7. It should be noted that for the surge inertia term X u an empirical value (Simonsen et al., 2012) is chosen, which is about 5% of the non-dimensional mass of the ship. The non-dimensional longitudinal gravity center x G , ship mass m', and moment of inertia I z are calculated (see Table 2) and also given in Table 7. Figures 11 and 12 show the simulation results (CFD-Present) of 35°turning circle and 10°/10°and 20°/20°z igzag maneuvers by using the hydrodynamic derivatives listed in Table 7 (Otzen & Simonsen, 2014;Simonsen et al., 2012) are also presented for comparison. In PMM-FORCE, the full set of PMM tests (3-DOF) was conducted to evaluate the maneuvers using the standard approach, while in CFD-FORCE, based on the PMM-FORCE data, a reduced test matrix was applied for the PMM computations with the same RANS code as the present study. Both PMM-FORCE and CFD-FORCE were conducted to determine the hydrodynamic derivatives in an integrated mathematical model, which is different from the mathematical model adopted in the present study. According to the SIMMAN Instructions In Table 8, the maneuvering parameters advance (A D ), tactical diameter (D T ), and overshoot angles (OSA) predicted by the system-based methods (CFD-Present, CFD-FORCE, and PMM-FORCE) are compared with those from FRMT (FRMT-MARIN). Table 8 also gives the relative comparison error E (%) of maneuvering parameters from the three system-based methods compared with those from FRMT, defined as (1 − Pred./FRMT)×100%, where 'Pred.' denotes the parameter predicted by the system-based methods.

CFD-FORCE)
From Figure 11 it can be seen that all the system-based prediction methods over-predict obviously the advance  Figure 11. Trajectories of ± 35°turning circle maneuvers.
and tactical diameter compared with FRMT. The large difference in turning circle maneuver between FRMT and those obtained from virtual or physical captive model tests may be ascribed to many reasons, such as the differences in the selected mathematical models, effect of 'solidity of captive model tests', the roll motion, the differences of the propeller models, and model scale effects, etc. (Quadvlieg, 2017). As can be seen in Table 8, all the relative comparison errors of both advance and tactical diameters by the CFD-Present method are smaller than those by CFD-FORCE, but larger than those by PMM-FORCE, except for the tactical diameter in the −35°t urning circle maneuver. The maximum and minimum comparison errors of the tactical diameter compared with FRMT occur in the simulations of CFD-FORCE within about −41% for the PS and −24% for the SB, respectively. The comparison error of the tactical diameter in the simulations of CFD-Present for the SB is about −25%, just slightly larger than that of CFD-FORCE. As for the advance, based on an overall quantitative comparison for the two CFD methods, better agreement can be achieved by the CFD-Present method within about −13% minimum relative error for the PS and −12% for the SB. Moreover, all the relative comparison errors of advance (maximum value of E equal to −17%) are smaller than that of tactical diameter (maximum value of E equal to −41%), indicating that the advance is better predicted.
In Table 8, the comparison errors of OSAs by the three system-based methods are positive, demonstrating that all the system-based methods under-predict the 1st and 2nd OSAs compared with FRMT. From Figure 12(a) (10°/10°) and Figure 12(b) (−10°/10°), it can be seen that compared with CFD-FORCE and PMM-FORCE, the 1st and 2nd OSAs predicted by the CFD-Present method show the best agreement with the minimum relative comparison errors of 2.3% and 15.2% for 10°/10°, 9.41% and 26.15% for −10°/10°zigzag maneuvers, respectively. The 1st and 2nd overshoot times predicted by the CFD-Present method are much closer to FRMT. From Figure 12(c) (20°/20°) and Figure 12(d) (−20°/20°), it can be seen that the comparison errors of the three systembased methods show little difference for both 1st and 2nd OSAs. Although the CFD-Present method gives a worse prediction of 1st and 2nd OSAs than PMM-FORCE, it predicts better than FORCE-CFD.
The cause of the differences in maneuvering motion predictions between CFD-FORCE and CFD-Present might be the test matrix reduction in CFD-FORCE and the distinction of the adopted mathematical models between the FORCE methods and present method. Besides, different numerical errors by different CFD methods could decrease the accuracy of maneuvering motion predictions to different extents.

Conclusions
By using an unsteady RANS solver, virtual captive model tests including static rudder, static drift, rudder-drift, pure sway, pure yaw, yaw-rudder, and yaw-drift tests at different captive speeds are carried out for the KCS ship model. Firstly, numerical uncertainty analysis is conducted. It is concluded that the values of numerical uncertainty are relatively small for the lower order Fourier coefficients, hence the lower order Fourier coefficients are preferred during the computation of the hydrodynamic derivatives. The computation results of nondimensional longitudinal force, transverse force, and yaw moment are validated by comparing with the available experimental data. It is demonstrated that the hydrodynamic forces and moments acting on a ship model in the captive model tests can be accurately predicted by using the proposed numerical model except for the longitudinal force in some cases mainly due to the simplified body force model for the propeller.
A full set of 65 linear and nonlinear hydrodynamic derivatives in the 3rd-order Abkowitz model are determined from the virtual captive model test. By using the hydrodynamic derivatives, 35°turning circle and 10°/10°, 20°/20°zigzag maneuvers are simulated, and the maneuvering characteristics of the ship are predicted. For validation of the prediction procedure developed in the present study, the numerical results are compared with the data of free running model tests from MARIN and those of CFD computations and PMM tests from FORCE Technology. From the comparisons it is confirmed that although the present method is less accurate than the PMM-based method overall, it gives a better prediction than the CFD method of FORCE Technology. These results show that the present method is effective in determining the hydrodynamic derivatives in the mathematical model of ship maneuvering motion, and can be used as an alternative to the captive model tests for predicting ship maneuverability at the ship design stage. However, it should be pointed out that the body force model seems to underestimate the changes of the propeller thrust and torque due to oblique propeller inflow in large drift angle and low speed simulations. Further analysis about the body force model and full propeller model in oblique flow is needed to get more accurate results. Moreover, the present work focuses on 3-DOF maneuvering motions and neglects the roll effect; while in fact ship roll is non-negligible during ship maneuvering, especially for a container ship. Therefore, a 4-DOF model with roll effect should be adopted in a future study. In addition, the present research is carried out at model scale and needs to be extended to full scale. In a future study, the virtual captive model tests can be conducted for models with different scales to investigate possible scale effects on the maneuvering derivatives. In this way, full-scale prediction can be achieved by extrapolation.