A fully nonlinear, dynamically consistent numerical model for solid-body ship motion. I. Ship motion with fixed heading

In this paper, we describe the details of our numerical model for simulating ship solid-body motion in a given environment. In this model, the fully nonlinear dynamical equations governing the time-varying solid-body ship motion under the forces arising from ship–wave interactions are solved with given initial conditions. The net force and moment (torque) on the ship body are directly calculated via integration of the hydrodynamic pressure over the wetted surface and the buoyancy effect from the underwater volume of the actual ship hull with a hybrid finite-difference/finite-element method. Neither empirical nor free parametrization is introduced in this model, i.e. no a priori experimental data are needed for modelling. This model is benchmarked with many experiments of various ship hulls for heave, roll and pitch motion. In addition to the benchmark cases, numerical experiments are also carried out for strongly nonlinear ship motion with a fixed heading. These new cases demonstrate clearly the importance of nonlinearities in ship motion modelling.


Introduction
For many decades, much effort has been devoted to modelling a ship's motion at sea in naval architecture (Lewis 1989). Accurate prediction of ship motions in real time is critical for improving operability and preventing large-amplitude ship motions, e.g. vessel capsize. This is particularly important to today's modern fast passenger ferries, high-powered naval vessels and cargo ships.
The earliest effort on modelling ship motion can be traced back more than half a century. For example, St Denis & Pierson (1953) first proposed a method to predict the statistics of ship responses to a realistic seaway. Using spectral methods developed in applied mathematics, they established a relationship between the spectral density of ship responses and the input ocean wave spectrum. Later, Peters & Stoker (1967) developed a different algorithm for thin ship body motion: a first-order theory using a systematic perturbation procedure with the ship's beam and unsteady motion assumed *Author for correspondence (ray.lin@navy.mil). comparable (of similar magnitudes). This approach was further improved by Newman (1961) with a set of small parameters and more appropriate body boundary conditions. Later, the strip theory was developed by Ogilvie & Tuck (1969) in which linear ship motion coefficients are introduced with the slender-body assumption. These coefficients include added mass and damping coefficients used in heave and pitch motions. The surface integration used in this theory complicates its computational implementation. The theory is also inconsistent: the formulation is applicable in the short-wavelength domains, while the slender-body approximation works in the long-wavelength domains. An interpolation theory (Maruo 1970) and a unified theory (Newman 1978) were then proposed to reduce this inconsistency.
Introduction of the more complex (but still linear) Neumann-Kelvin approach in numerical modelling greatly advanced the ship motion research (e.g. Dawson 1977;Magee 1994;Shin et al. 1997), because of its capability to model arbitrary ship surfaces. However, it has its own limitations, e.g. in solving a forward-speed Green function in finite water depth (Beck & Reed 2000).
When nonlinear terms are comparable to, or even stronger than the linear terms, new models are then needed. As an intermediate step to fully nonlinear models, several 'blending methods' were developed. For more details, we refer the reader to the ISSC report on extreme hull girder loading (ISSC 2000).
A different approach was introduced by Wilson et al. (1998) and Gentaz et al. (1999) to solve the Reynolds averaged Navier-Stokes (RANS) equations in the time domain in which an iterative method is used for steady solutions, and a timestepping method is used for unsteady solutions. The outcome of this approach is inconclusive, partly owing to insufficient numerical results. More seriously, it faces convergence problems for strongly nonlinear ship motion.
Several other ship motion models were developed based on potential flow, e.g. Lin et al. (1986). However, they only partially include nonlinearities, and also use many empirical and free parameters, thus limiting their applications to ship design.
In a more recent review by Beck & Reed (2000), approximately 80 per cent of ship motion models are based on the strip theory for its simple numerical implementation and its flexibility on ship hull forms. This implies also that most models are not suitable for strongly nonlinear ship motion.
To model strongly nonlinear ship motion, it is therefore necessary to develop a methodology completely different from the above traditional approaches. In particular, it should be dynamically and mathematically consistent, numerically efficient, and independent of a priori experimental data.
For this purpose, we have developed a new generation ship motion model. In the first phase of our effort (Lin & Kuang 2004, 2006Lin et al. 2005), we developed a 'steady ship motion model' in which only the nonlinear interactions of surface waves (including ship motion-generated waves) are simulated. The response of the ship hull to these interactions (e.g. dynamic pressure force and the buoyancy force from the displaced water) is not included in that model. In the rest of this paper we call it the phase 1 model.
To continue this effort, we expand the model to include an independent component for simulating six-degrees-of-freedom, solid-body motion under the forces calculated by the phase 1 model. A challenge in this effort is to accurately model these forces when the ship position (e.g. underwater volume) varies in time; that is, the position variation is determined by the six-degrees-of-freedom ship solid-body motion, but the motion changes are determined by these forces. To address these technical difficulties, we need to develop a suite of new numerical methodologies for this component. In particular, these methods must be numerically convergent and efficient to model dynamic dependences of the ship position and ship motion. The combination of this component and the steady ship motion model is called the 'Digital, Self-consistent, Ship Experimental Laboratory Ship Motion Model' (DiSSEL).
In this paper, we report for the first time the mathematical details and the numerical benchmark results of DiSSEL, focusing in particular on the solidbody motion component. The paper is organized as follows: the mathematical equations of ship motion are in §2. Section 3 provides the benchmark results with respect to experiments and independent numerical models. Discussions are given in §4.

Mathematical model
This section provides the details of the mathematical equations, numerical algorithms and formulation. The description is focused on the solid-body motion component. There is only a brief summary of the phase 1 model; its details are in Lin & Kuang (2004, 2006 and Lin et al. (2005).

(a) Model reference frames
To integrate the solid-body motion component with the phase 1 model, one must examine the dynamical equations in two different reference frames. One is moving horizontally with the mass centre of the ship, with the origin set at the mean free surface (hereafter called the model reference frame). The other is attached to the ship body, with the origin at the ship mass centre (hereafter called the ship reference frame). The first is used in the phase 1 model. The latter is convenient for calculating the ship solid-body motion. The moment of inertia of the ship and the numerical grid of the hull are invariant in this reference frame.
Ship movement can be decoupled into a three-dimensional translational motion (u c + v c ) of the mass centre x c of the ship, and a three-dimensional rotation motion U about x c . In this description, u c is the ship velocity vector moving in the calm water, and v c is the translational motion vector responding to surface waves, where v x , v y , and v z are the surge, the sway, and the heave motions, respectively. Similarly, the angular velocity vector U is In naval engineering, U x , U y and U z are called roll, pitch and yaw, respectively.
In the model reference frame the x-axis is the ship heading direction (from the stern to the bow), the y-axis points to the port and the z-axis is upward.
The ship mass centre x c = (0, 0, z c ) moves vertically with the velocity i.e. the heave motion of the ship. In addition, the reference frame moves horizontally with the velocity The ship reference frame is often convenient to solve for the ship solid body motions (e.g. Goldstein 1980). In this reference frame the position vector is denoted by where the subscript s implies the quantity defined in the ship reference frame, and the superscript T means the transpose of the vector. Therefore, the position vector x in the model reference frame and x s in the ship reference frame can be transformed from each other via where A is the transformation matrix and can be expressed with the Euler angles q (2.3) Other Euler angle definitions can also be used as long as no two consecutive angles are defined relative to the same axis (Goldstein 1980). In particular, A −1 = A T , and its time variation is d dt For a rigid ship its boundary x S s does not vary in time. However, in the model reference frame, its time variation is given by In DiSSEL, equations (2.1), (2.4) and (2.5) are used to update the ship boundary x S .

(b) Hydrodynamic equations of the surface waves
The surface waves are solved in the model reference frame. The fluid is assumed inviscid and incompressible. Therefore, the fluid flow can be described by a single velocity potential 4, and can be solved by the following nonlinear hydrodynamic equations: in −H ≤ z ≤ h (h is the free surface elevation and H is water depth). In addition, there are three boundary conditions In equation (2.7), g is the gravitational acceleration vector, p the total pressure, V the gradient operator and V h the horizontal gradient. In equation (2.10),n S is the normal unit vector of the ship surface x S . The acceleration term in equation ( is from the acceleration of the model reference frame. The inviscid fluid approximation works well for large-scale motions. However, dissipation is important for small-scale flow (also called the sub-grid processes if their spatial scales are smaller than numerical resolutions), which arises from nonlinear interactions of large-scale waves, and/or from other turbulent processes, such as wave-breaking processes. To accommodate this, a small dissipative coefficient n is introduced in equation (2.7). It also acts to ensure numerical convergence with moderate resolution while retaining correct large-scale flow structures of interest. For more discussion on this dissipative effect, we refer the reader to . It should be pointed out that this approach has been commonly used in computational fluid dynamics for many decades (e.g. Spiegel & Veronis 1960;Yanai 1983). It was also successful in our earlier simulation (Lin et al. 2005;Lin & Kuang 2006).
In addition to the boundary conditions (2.8), (2.9) and (2.10), we also need conditions on the numerical domain boundaries. In our model, the far-field radiation conditions (as x → ∞) are used on the borders of the numerical domain that is confined by the forward boundary S f and the side and aft boundaries S a . Thus, the radiation boundary condition on S f is where 4 e and h e are the velocity potential and the free surface elevation of the environment (which vanish in calm water). The condition (2.11) implies that the waves generated by the ship do not radiate ahead of the ship (provided that the ship speed is faster than the wave propagation speed). Open boundary conditions are applied on S a that permit waves to travel freely across the boundaries where the over line represents the spatial average along the boundaries. Equation (2.12) also ensures conservation of mass for incompressible flow.

(c) Dynamic equations for the ship solid-body motion
A rigid ship motion can be decomposed into the translational motion v c of its mass centre x c and the rotation motion U about x c (Goldstein 1980). The mass centre of the ship is defined by where m s is the total mass of the ship, V s the total ship volume and r s the density of the ship. In particular dx c dt = u c + v c . (2.14) In the model reference frame, x c = (0, 0, z c ), and the variation of z c is given by equation (2.1).
The kinematics of the rotation in the model reference frame is also very simple. Given any point x S on the ship boundary, we have Equations (2.14) and (2.15) describe the six-degrees-of-freedom, solid-body motion of the ship. The dynamics of the ship motion is governed by the classical mechanics equations: where I is the tensor of the moment of inertia of the ship; D v and D U are the integrated dissipative effects (e.g. drag) of the fluid on the ship hull, and F and G are the net force and the net torque (moment) on the ship hull, respectively. The net force (moment) includes contributions from the dynamic pressure on the ship hull, and from the displaced water In equations (2.18) and (2.19), r w is the water density, x wet the geometric centre of V wet , S wet and V wet are the wetted surface and the underwater volume of the ship, respectively. The moment of inertia is (d is the Kronecker delta function).
It should be pointed out that the two dissipative coefficients, D v and D U , are from various sources, such as friction (ignored in the large-scale potential flow approximation) owing to small-scale flow generated by, for example, bilge keels and wave-breaking. In many cases, the dissipation owing to the bilge keels is very strong, sometimes reaching approximately 20-25% of the net moment on the ship. The small-scale flow, and thus the dissipation, depends on the incident waves and the ship profiles as well (Lin & Kuang 2008). Friction from wavebreaking is in general weak and contributes only approximately 1 per cent of the ship acceleration (Lin & Kuang 2007).
In numerical modelling of the ship motion accurate computation of the forces and torques (moments) is critical. This depends on several important factors. By equations (2.18) and (2.19), we can observe that accurate determination of the forces and the torques (moments) depends on two issues: an accurate wetted surface S wet (and therefore the underwater volume V wet ), and an appropriate numerical algorithm for surface integration. In addition, the dissipative effect must also be modelled consistently.
The wetted surface S wet varies in time, and can be determined if x c and q are known. They can be determined if F and G are given, but F and G depend on S wet via equations (2.18) and (2.19). Therefore, the relationships among the system variables are very nonlinear. In other words, all degrees of ship motion couple to each other.

(d) Numerical methods
In DiSSEL, the numerical state vector S is A hybrid scheme is used in the model to solve S. In this method, the velocity potential 4 and the free surface elevation h are approximated by finite Fourier expansions, e.g.
The ship surface x S s is described by either a finite-element mesh (unstructured grid) or a finite-difference mesh (structure body fix grid), and is fixed in the ship reference frame.
The numerical workflow is as follows: First given the initial state vector S(t 0 ), the nonlinear wave-wave interactions are solved via fast Fourier transform (FFT) transform between the spectral space and the physical space. The wetted surface S wet at t 0 is determined from S(t 0 ). Then the forces F and the torques G are evaluated via equations (2.18) and (2.19). After those calculations, a third-order Runge-Kutta scheme is used to update the state vector S(t 0 + Dt).
The advantage of the hybrid method is obvious. In the ship reference frame, the mesh is fixed and invariant, and is therefore determined before the time integration. Updating of 4 and h is carried out by the pseudo-spectral algorithm. Therefore, the central processing unit (CPU) time (flops) at each time step is determined by the FFT transforms, which is of the order (N log N ) 2 . Compared with that of N 4 for the regular finite-difference scheme, this is computationally much more efficient. Another advantage of the hybrid method is its numerical accuracy. The ship hull is much smaller than that of the entire numerical domain. The grid sizes for surface wave calculations are not appropriate for the numerical integration of equations (2.18) and (2.19). However, the finite element/difference mesh x S s can provide sufficient resolution for the surface integration. Update of x S is determined by equation (2.5), i.e. by z c , U (and therefore A). Together with the updated free surface elevation h, we can update the wetted surface S wet . Update of the normal vectorn S can be calculated either from the updated x S , or more efficiently, directly from Convergence of the numerical solutions is discussed in detail in Lin et al. (2005) and in Lin & Kuang (2006). Details on the integrated dissipative effects in equations (2.16) and (2.17) are described in Lin & Kuang (2007, 2008.

Model benchmark results
To validate the DiSSEL, we have carried out simulations for many cases with either experimental data or results from independent numerical ship motion models, e.g. LAMP (Large Amplitude Ship Motion Programme; Lin et al. 1986). In these cases, the ship motion generated wave patterns, nonlinear wave profiles around the ship, pressure distributions and unsteady ship motions (roll, pitch and heave) are compared and cross-examined. Reported here are selected cases, which are the extreme scenarios of the benchmarked examples. A fixed heading is assumed in all benchmark cases. Reported cases include ship motions in Sea State 3 (SWH = 0.9144 m, T p = 10.35 s), Sea State 4 (SWH = 1.9 m, T p = 10.75 s) and Sea State 5 (SWH = 3.25 m, T p = 11.85 s; note: SWH, significant wave height; T p, wave period).

(a) X-craft catamaran
The X-craft catamaran is an example of a high-speed ship. Its Froude number reached Fr = 0.77 in experiments, and thus is appropriate for testing the model's capabilities for very fast ship motion (one of the challenges for ship motion modelling). In this benchmark case, we simulated ship motion with three different Froude numbers: Fr = 0.23, 0.38, 0.77. Our results are then compared with experimental data and simulation results from an independent ship motion model LAMP. Figure 1 shows the results of the pitch motion in head seas for Sea State 4 and Fr = 0.77 from DiSSEL (solid line), experiments (dashed line) and LAMP (circles with dashed line). Compared with the LAMP results, our simulations are much closer to the experimental data. In particular, the amplitudes of the two sets of results are comparable. Better agreement is found for smaller Froude numbers as does the agreement between LAMP and experiment improve in these cases. Similar conclusions apply to the roll motion.
We also plot in figure 2, the ship motion-generated wave distributions for the three Froude numbers. Though there is no experimental measurement for the cases, there are results from other ship motion models (only for small Froude numbers), which our results agree well with (Lin et al. 2005;Lin & Kuang 2006). The wave distributions in figure 2 demonstrate clearly the robustness of our ship motion model for high-speed ships.

(b) HSS trimaran
The HSS Trimaran is another type of high-speed craft. We simulate motion of this craft in bow seas (120 • between ship heading and incident waves) with Sea State 5 and Fr = 0.51. In this case, there are both roll and pitch motions and it is thus ideal for understanding (nonlinear) interactions between different components of ship solid-body motion.
The results for the pitch motion and the roll motion are shown in figures 3 and 4, respectively. The line types for different results are defined similarly to those in figure 1.
The simulation results (from DiSSEL and LAMP) of the pitch motion are similar and agree well with the experimental data, as shown in figure 3, but the roll motion results are different. The time series from DiSSEL agree well with that of the experiment (in both amplitude and frequency). However, the agreement between that of the LAMP and the experiment is poor and becomes satisfactory only for smaller Froude numbers (Fr < 0.4). The benchmark results indicate that accurate simulation of roll motion is much more difficult than that of pitch motion. This is because the roll motion depends strongly on nonlinear shipwave interactions. These interactions become stronger as the forward speed of the ship increases. The good agreement between DiSSEL and experimental results for both pitch and roll motions demonstrates that these nonlinear interactions are properly implemented into DiSSEL.
(c) Planning craft propulsion system demonstrator Ship motions of the propulsion system demonstrator (PSD) in the highspeed planning regime are an example of a strongly nonlinear ship motion problem. In particular, the craft's underwater volume is reduced substantially from the displacement regime (low speed) to the planning regime (high speed). Consequently, maintaining dynamical balances accurately at different stages is a numerical challenge. Since the balances determine the trim angle and the vertical displacement, numerical errors in these balances will lead to erroneous and often divergent solutions. Because of the difficulty, there has not been any ship motion model capable of modelling this craft in the planning regime. We select this benchmarking case to demonstrate the capabilities of DiSSEL to simulate strongly nonlinear ship motion. The mean values of the vertical displacement and trim angle related to the ship forward speed are approximately equal to those in calm water (e.g. , and the heave and pitch motions in the following figures only include those arising from the ship-wave and wave-wave interactions. Figures 5 and 6 show the heave and pitch motions at Fr = 0.526 (the planning regime) and Sea State 3 in head seas, respectively. The numerical time series (the solid line) agrees nearly perfectly with the experimental data (the dashed line) for both pitch and heave motions. In addition, we have carried out simulations with even higher designed planning speeds (Fr ≤ 1.56). The results for Fr = 1.56 and Sea State 3 are shown in figures 7 and 8. There are no experimental measurements for benchmarking at these conditions. However, the numerical results are reasonable because the heave motion follows the incident waves with a small phase lag (as shown in figure 7), consistent with the empirical conclusions from experiments. Therefore, the results here can be used as 'prediction' for future experiments.
In figures 9 and 10, we plot the force/moment (from the dynamic pressure) and the restoring force/moment (from the buoyancy) from the simulations. From these figures, we can observe that the two kinds of forces (moments) are comparable in magnitude, but nearly opposite in direction, indicating the establishment of the dynamical balance in the numerical solutions. Planning craft sometimes exhibit low freeboard and bow submergence, often called 'Plow-In', as shown in figure 11 by the DiSSEL simulation. Plow-In occurs when the frequency of the ship motion generated waves is similar to the ship natural frequency. Numerical simulations show that Plow-In occurs when the forward speed is between 2.83 m s −1 (Fr = 0.27) and 3.24 m s −1 (Fr = 0.31). This condition is very similar to that observed in experiments . We want to emphasize here that this is the first successful numerical simulation of Plow-In.

(d) ONR hull form
Ship motion is also sensitive to the hull geometry, under and above the water (defined in calm water). The motions of two ship hulls with identical underwater geometry (in calm water) can be very different if the above-water geometries are not the same. This difference is significantly amplified in the resonant state in which the incident wave frequency is the same as the ship natural frequency, as shown in experimental results (figure 12). Simulation of the ship motion therefore requires correct evaluation of the forces and the moments on motion-dependent (thus time-varying) wetted ship surface. To demonstrate DiSSEL capabilities of resolving forces and moments with time-varying wetted surfaces, we simulate the roll motions of the ONR series Tumblehome hull and Flare hull. Both have the same underwater geometry (in calm water), but the above water geometries are very different; one is slanted inward and the other is outward. Figure 12 shows the results of the normalized roll motions (scaled by the wave steepness) of the two hulls. From the figure, we can observe that the numerical solutions agree very well with the experimental data. We can also observe that the difference in the normalized roll motion of the two hulls reaches a maximum at resonance. Obviously, one could not obtain the correct results if the underwater geometry is assumed invariant (as often used in the linear ship motion model).

Discussion
In this paper, we provide an overview of the DiSSEL ship motion model, and the detailed description of the ship solid-body motion component. Included in this model are all wave-wave and ship-wave interactions. Consequently, all six components of the ship solid-body motion are dependent on each other and any change in one component will affect the rest via nonlinear interactions of the system. In addition, dissipative processes from the blocking effect (Lin & Kuang 2008) and from sub-grid effects (e.g. turbulence, wave-breaking, subscale processes) on the ship motion are included in the dynamical equations (2.16) and (2.17). In particular, the dissipative effect is formulated based on fundamental physics and mathematical properties. There is no other parametrization (e.g. added mass) or linearization applied in our model. With our approach, a single velocity potential is needed to model all nonlinear interactions.
The numerical algorithm used in our model is a hybrid of pseudo-spectral and finite-difference/finite-element algorithms, capable of resolving very complex ship hull geometries and broad surface wave spectra. The algorithm is computationally efficient and robust (Lin & Kuang 2006).
Many benchmarking tests have been performed with the DiSSEL model. Reported in this paper is only a small fraction of them, which are chosen for the following reasons: (i) well-documented experimental and/or independent numerical simulation results; (ii) strongly nonlinear effects in ship motion; and (iii) numerically challenging problems. We also provide new simulation results (case c with Fr = 1.56) for future benchmark and prediction purposes. In particular, we would like to point out that the DiSSEL model can be used to simulate ship motion with any arbitrary forward speed, i.e. with the Froude number values far larger than those reported in this paper.
In all benchmark cases our simulation results agree well with experimental data, while the agreement of independent models varies. In general, our simulation results outperform those of independent models.
Our work is a step further towards providing a realistic ship manoeuvring model. By integrating the ship solid-body motion component and the surface wave model, we are able to determine ship motion, the motion generated waves and the responses to ship-environment interactions concurrently. However, the current version of DiSSEL is not complete. Future development includes simulation of ship motion for a ship manoeuvring in a seaway.