Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012 (2012), Article ID 563182, 21 pages
Research Article

On Nonlinear Simulation Methods and Tools for Evaluating the Performance of Ships and Offshore Structures in Waves

Ship Design Laboratory, National Technical University of Athens, 15773 Athens, Greece

Received 20 January 2012; Revised 1 June 2012; Accepted 8 June 2012

Academic Editor: Ioannis K. Chatjigeorgiou

Copyright © 2012 Shukui Liu and Apostolos Papanikolaou. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


This paper describes the development of alternative time domain numerical simulation methods for predicting large amplitude motions of ships and floating structures in response to incoming waves in the frame of potential theory. The developed alternative set of time domain methods simulate the hydrodynamic forces acting on ships advancing in waves with constant speed. For motions’ simulation, the diffraction forces and radiation forces are calculated up to the mean wetted surface, while the Froude-Krylov forces and hydrostatic restoring forces are calculated up to the undisturbed incident wave surface in case of large incident wave amplitude. This enables the study of the above waterline hull form effect. Characteristic case studies on simulating the hydrodynamic forces and motions of standard type of ships have been conducted for validation purpose. Good agreement with other numerical codes and experimental data has been observed. Furthermore, the added resistance of ships in waves can be calculated by the presented methods. This capability supports the increased demand of this type of tools for the proper selection of engine/propulsion systems accounting for ship’s performance in realistic sea conditions, or when optimizing ship’s sailing route for minimum fuel consumption and toxic gas emissions.

1. Introduction

The accurate prediction of the seakeeping behavior of ships and offshore structures in real seas is a demanding task for naval architects and of great practical interest to shipbuilders, owners/operators, as it affects both their design and operation.

Quasi 2D strip theory approaches to the seakeeping of ships were the first which delivered satisfactory results for practical applications to the prediction of wave-induced loads and motions, and they are widely used even today. With the rapid advance of computer technology in the 70ties, various frequency domain 3D approaches were also developed. But due to some inherent limitations, their successful applications are limited to certain extent. On the other hand, the time domain simulation methods become more and more popular as they enable the address of large amplitude ship motion problems which is very important for the design and the assessment of safe operation of modern ships and offshore structures operating in a variety of adverse environmental conditions.

Following the pioneering work of Finkelstein [1] and Cummins [2], many researchers investigated seakeeping problems by different time domain approaches and showed promising results for both linear and nonlinear problems of different complexity level. Lin and Yue [3] showed the applicability of a time domain Green function method to large amplitude ship motions. Following this formulation, Singh et al. [4] appeared to have obtained good results in some applications. However, this could not be confirmed in some other practical cases by other researchers; namely, when this method is applied to floating bodies with a flare at the waterline, which is common to modern ship designs, numerical problems may arise and computations fail. Duan and Dai [5] found that the commonly used panel method employing the transient Green function for a non-wall-sided floating body does not satisfy the mean-value theorem of definite integrals for the near-water surface panels and solved this problem by introducing an imaginary vertical surface, which encloses the hull surface in the fluid domain. This method works fine, unless the body has some bulb-like hull form, which exceeds the projection of the water plane. Zhang [6] used a similar scheme, but the matching surface is placed some distance away from the body and moving at the same speed as the ship, so that there should not be any problem with body’s shape. In this scheme, a part of the free surface is included in the inner domain, which needs updating at each time step. Yasukawa [7] and Kataoka and Iwashita [8] also used similar schemes to solve the problem. The difference between these methods basically lies in the different ways of treatment of the boundary condition in the far field, in simulating the free surface and the numerical schemes for solving the core equations. More recently, there is a trend of integrating CFD techniques into the hybrid method so as to study highly nonlinear phenomena. Iafrati and Campana [9] presented a hybrid method combining a CFD scheme using conventional grids and the BEM for potential-flow free-surface problems. Sueyoshi et al. [10] use particle methods in the inner domain and a boundary element method in the outer domain to study various wave-free surface problems. Lin et al. [11] presented recently a paper where they combined a viscous flow solver in the inner domain and potential flow solver in the outer domain. There is a slight overlap between the two introduced domains, which creates a matching domain.

The objective of this paper is to give an overview of recent developments of time domain numerical simulation methods developed at the Ship Design Laboratory of NTUA for predicting large amplitude motions of ships and floating structures in response to incoming waves in the frame of potential theory. Addressed shiplike bodies are assumed with zero or nonzero constant forward speed. Developed method(s) and related software tool(s) prove applicable to ship design and the assessment of the operation of ships and offshore structures in seaways.

In the course of the research, which formed the major part of the Ph.D. work of the first author [12], the authors developed a new time-domain transient Green function method and demonstrated its applicability by solving some fundamental hydrodynamic problems [13]. However, some inherent limitations of this method surfaced in the study of realistic shiplike forms with flared sections around the waterline. Thus, a hybrid method concept was further developed for the simulation of the hydrodynamic forces acting on realistic ship hull forms in waves [14, 15]. In this method, the fluid domain is decomposed into an inner and an outer part. The Rankine source method is applied in the inner domain, while the transient Green function method is used in the outer domain. This hybrid method works well with a relatively small number of panels compared to a pure Rankine source method, noting that the free-surface panelization is restricted between the body boundary and the control surface.

The calculation of the force components in the motion simulation is based on the following assumptions: if the incident wave amplitude is small (quasilinear case), then the Froude-Krylov and hydrostatic restoring forces are calculated by using ship’s geometric/hydrostatic data referring to ship’s mean floating position, whereas, if the incident wave amplitude is large, then the above force components are calculated up to the undisturbed incident wave surface by pressure integration over the instantaneous wetted surface of the moving body. The diffraction forces and radiation forces, however, are calculated up to the mean wetted surface by the developed hybrid method for all incident wave amplitude cases. Calculated force components are introduced into the equations of motions to predict the ship motions in the time domain.

In order to validate the developed theoretical methods and numerical schemes (computer codes), several characteristic case studies on the hydrodynamic forces and motions of standard type ships have been conducted. Results are compared with those of other authors and available experimental data.

During the conducted validation it was observed that the developed hybrid method, which is formulated in the earth-fixed coordinate system, is quite time consuming. Thus another scheme has also been developed as an alternative, namely, a hybrid method formulated in the body-fixed coordinate system. This method works much faster, despite the fact that an additionally arising line integral term in the formulation of the relevant boundary value problem complicates the calculations. This new approach has been applied to the study: the seakeeping performance of a modified Wigley hull with different above-water shapes at nonzero speed. Obtained results show the effect of variation of the above-waterplane ship hull form on motions at reasonable computational time; thus it proves that the developed method can be useful for ship design optimization and the assessment of ships’ performance in high seas.

Finally, an important part of the developed set of methods and software tools refers to the calculation of the added resistance and required powering of ships in waves [16]. This type of calculation proved recently very important in ship design and the assessment of ship’s operation in waves, namely, when selecting ship’s engine/propulsion system and considering ship’s performance in terms of sustainable service speed in realistic sea conditions. It also affects ship’s operation when optimizing the sailing route for minimum fuel consumption and toxic air emissions in view of green-sailing considerations [17].

2. Formulation

Consider a 3D floating body moving on the free water surface with speed 𝑈0 and undergoing 6 DOF motions in response to a regular, harmonic wave. An earth-fixed Cartesian coordinate system 𝑂-𝑋𝑌𝑍 is chosen with the 𝑋-𝑌, horizontal plane coincident with the undisturbed free surface, and the 𝑍-axis pointing upwards through the ship’s initial mass centre, as shown in Figure 1. The fluid is assumed to be ideal and the water depth infinite.

Figure 1: Coordinate system definition and decomposition of the fluid domain.

The motions of the ship may be determined through the motion of the body-fixed system 𝐺𝑥𝑦𝑧 relative to body-travelling system 𝑂𝑥𝑦𝑧. A total of six components are needed to define the motions, typically three translations, that is, surge, sway, and heave, and three rotations, that is, roll, pitch, and yaw. Following Newton’s 2nd Law, the six DOF motions of the travelling, rigid body in space are determined by solution of the following set of equations: 𝑑𝐏𝑑𝑡=𝑚𝑑𝐯𝑑𝑡=𝐅=𝐅𝐻𝑆+𝐅𝐈+𝐅𝑅+𝐅𝐷+𝐅𝑀𝑅𝑑𝐋𝑑𝑡=𝐋𝑑𝝎𝑑𝑡+𝝎×(𝐈𝝎)=𝐌=𝐌𝐻𝑆+𝐌𝐈+𝐌𝑅+𝐌𝐷,(2.1) where 𝐏=𝑚𝐯 and 𝐋=𝐈𝝎 are the linear momentum and angular momentum, respectively, 𝐅 and 𝐌 are the total force and the total moment about the mass center on the body, 𝑚 is the body’s mass, and 𝐯 is the absolute velocity vector of the mass center 𝐺 in the 𝑂𝑋𝑌𝑍 system; 𝐈 and 𝝎 are the inertia tensor and angular velocity about ship’s mass center 𝐺. The moments and products of inertia in 𝐈 are constants in the moving and rotating system 𝐺𝑥𝑦𝑧. The position of a rigid body in space is fully determined by the position of 𝐺 and the angular orientation of the 𝐺𝑥𝑦𝑧 system with respect to the body-travelling system 𝑂𝑥𝑦𝑧. In solving these equations, the position of the ship’s mass centre 𝐺 in the earth-fixed system is defined by the vector 𝑥𝐺(𝑡)=[𝑥𝐺,𝑦𝐺,𝑧𝐺]𝑇 and its velocity is expressed by the time derivative: 𝐯(𝑡)=𝑑𝐱𝐺(𝑡)/𝑑𝑡. To express a vector in the system 𝐺𝑥𝑦𝑧 with respect to 𝑂𝑋𝑌𝑍 system, as necessary for the large motion simulations, the use of relevant transformation matrices expressed by the Euler angles is needed [18].

On the right-hand side of (2.1), there are forces due to incident wave, diffraction, radiation, restoring, and other possible terms. Basically these forces/moments are calculated by integrating the pressure expressed by Bernoulli’s equation on the body surface. Other possible forces acting on the body, like mooring forces (in case of a zero speed, moored floating structure problem), can be added [19]. Since an exact, fully nonlinear model is quite time consuming and complicated for numerical computation, we restrict ourselves in simplified nonlinear simulations and the consideration of some of the most important and more tractable nonlinear effects. In particular, for simulating large amplitude motions, the incident wave forces (Froude-Krylov) and restoring forces are computed exactly over the undisturbed instantaneous wetted body surface and transferred into the motion equations, whereas the inherent hydrodynamic forces, namely, radiation and diffraction forces, are calculated up to mean wetted surface by using Bernoulli’s equation once the corresponding velocity potentials are solved.

2.1. Radiation and Diffraction

In the frame of potential theory, the unsteady flow field can be described by a velocity potential: Φ𝑇(𝑥,𝑦,𝑧,𝑡)=Φ0(𝑥,𝑦,𝑧,𝑡)+Φ(𝑥,𝑦,𝑧,𝑡),(2.2) where Φ0 is the incident wave potential, Φ is the disturbed flow potential, 𝑡 is time, and 𝑝(𝑥,𝑦,𝑧) is a point in the flow field. In the fluid domain Ω(𝑡), Φ(𝑝,𝑡) satisfies Laplace’s equation and a set of boundary conditions: 2Φ(𝑝,𝑡)=0,(2.3)𝜕Φ𝜕𝑛=𝑣𝑛𝜕Φ0𝜕𝑛on𝑆𝑏,𝜕(𝑡),𝑡>0(2.4)2Φ𝜕𝑡2+𝑔𝜕Φ𝜕𝑧=0on𝑆𝑓,(𝑡),𝑡>0(2.5)Φ=𝜕Φ𝜕𝑡=0on𝑆𝑓,(𝑡),𝑡=0(2.6)Φ,Φ,𝜕Φ𝜕𝑡0on𝑆,,𝑡>0(2.7) where 𝑛 is the unit normal vector pointing out of the fluid domain Ω(𝑡) and 𝑣𝑛 is the normal component of the body motion velocity.

We transfer the solution of the above boundary value problem (BVP) ((2.3)–(2.7)) by use of Green’s 3rd theorem into a set of integral equations for the sought velocity potential. For the solution of the deduced integral equations, we employ the time domain transient Green function method (TDGF) or alternatively a derivative thereof, namely, the HYBRID method. Furthermore, the HYBRID method can be formulated either in the earth-fixed coordinate system or in the ship-travelling coordinate system; the different versions will be herein named as HYBRID I and HYBRID II.

First, the flow field is decomposed into two parts, that is, the inner domain denoted as I, enclosed by the body boundary, the control surface, and a part of the free surface, and the outer domain denoted as II, enclosed by the control surface, the remaining free surface, and the boundary surface at infinity (Figure 1).

The integral equation for solving the velocity potential by TDGF is as follows: 2𝜋Φ(𝑝,𝑡)+𝑆𝑏(𝑡)𝜕Φ(𝑞,𝑡)𝜕𝑛𝑞1𝑟𝑝𝑞1𝑟𝑝𝑞1𝑟𝑝𝑞1𝑟𝑝𝑞𝜕Φ𝜕𝑛𝑞𝑑𝑠𝑞=𝑡0𝑑𝜏𝑆𝑏(𝜏)𝐺𝜕Φ𝜕𝑛𝑞𝜕𝐺Φ𝜕𝑛𝑞𝑑𝑠𝑞+1𝑔𝑤𝑙(𝜏)𝐺𝜕Φ𝜕𝐺𝜕𝜏Φ𝑉𝜕𝜏𝑁𝑑𝑙𝑞,(2.8) where the time-dependent Green function is defined as 𝐺(𝑝,𝑡;𝑞,𝜏)=20𝑔𝑘𝑒𝑘(𝑧+𝜁)𝐽0(𝑘𝑅)sin𝑔𝑘(𝑡𝜏)𝑑𝑘,(2.9) where 𝑝(𝑥,𝑦,𝑧) and 𝑞(𝜉,𝜂,𝜁) are the field and source points; respectively, 𝐫𝑝𝑞=(𝑥𝜉)𝐢+(𝑦𝜂)𝐣+(𝑧𝜁)𝐤, and 𝐽0 is the zero-order Bessel function.

When applying the HYBRID method, a matching surface is introduced to split the whole fluid domain into an inner domain and an outer domain. In the inner domain I, the Rankine source method is used to solve the flow field. The integral equation takes the form of 𝑆𝑐+𝑆𝑏+𝑆𝑓ΦI𝜕(𝑞,𝑡)𝜕𝑛𝑞1𝑟𝑝𝑞1𝑟𝑝𝑞𝜕𝜕𝑛𝑞ΦI(𝑞,𝑡)𝑑𝑠𝑞=2𝜋ΦI(𝑝,𝑡)𝑝𝑆I,(2.10) where 𝑛 is the unit normal vector pointing outward of the inner domain. q(ξ,η,ζ,τ) is the source point, p(x,y,z,t) is the field point; the denotation 𝑆𝑏, 𝑆𝑐, and 𝑆𝑓 represent, respectively, the body surface, the control surface, and the free surface. In the outer domain II, the transient time-domain Green function is employed to solve the disturbed potential on 𝑆𝑐. In case of earth-fixed coordinate system, the integral equation is expressed as 2𝜋ΦII(𝑝,𝑡)+𝑆𝑐ΦII𝜕(𝑞,𝑡)𝜕𝑛𝑞1𝑟𝑝𝑞1𝑟𝑝𝑞1𝑟𝑝𝑞1𝑟𝑝𝑞𝜕ΦII𝜕𝑛𝑞𝑑𝑠𝑞=𝑡0𝑑𝜏𝑆𝑐𝐺𝜕ΦII𝜕𝑛𝑞ΦII𝜕𝐺𝜕𝑛𝑞𝑑𝑠𝑞𝑝𝑆𝑐.(2.11)

In order to decrease the computational time, in present work we use the linearized body boundary condition for surface-piercing bodies. In the forward speed case, the 𝑚𝑗 term is introduced to account for the forward speed effect on the boundary condition. The 𝑚𝑗 term can be either calculated directly or approximated by assuming that the steady flow may be represented by the undisturbed stream 𝑈0𝑥 [20].

2.2. Incident Wave and Restoring Forces

In classical linear seakeeping theory, the incident wave pressure is integrated over the mean wetted surface of the vessel. However, it has been observed in experiments and numerical studies that the nonlinear effects of the incident wave and of ship’s hydrostatic restoring forces are significant compared with those of other hydrodynamic force components, such as diffraction and radiation forces; thus it is important to include these effects in the motion simulation accurately. In this paper, when simulating the response to large amplitude incident waves, the pressure of the incident wave is taken into account up to the instantaneous undisturbed wave surface.

In the implementation of this concept, a very fine local panelization is implemented for the calculation of these two force terms. At each time instant, every panel’s geometric information is updated according to the motion results and the wetting of the panels under the instantaneous wave profile is determined by checking 𝑧𝜁(𝑥,𝑦).

2.3. Drift Force and Added Resistance

For the calculation of the quasi-second-order drift forces (zero speed case) and added resistance in waves, we follow [21] far-field approach. In the more general, non-zero forward speed case, the added resistance is expressed as 𝑅𝐴𝑊=𝜌8𝜋𝛼0𝜋/2+𝛼𝜋/203𝜋/2𝜋/2||𝐻𝑘1||,𝜃2𝑘1𝑘1cos𝜃𝑘cos𝜒+𝜌14Ωcos𝜃𝑑𝜃8𝜋2𝜋𝛼0𝛼0||𝐻𝑘2||,𝜃2𝑘2𝑘2cos𝜃𝑘cos𝜒14Ωcos𝜃𝑑𝜃,(2.12) where 𝜒 is the incident wave heading, 𝜃 is the elementary wave angle, 𝜌 is the density of sea water, and 𝛼0 is the critical angle (𝛼0=arcos(1/(4Ω)) for Ω>1/4 and 𝛼0=0 for Ω1/4).

The complex function 𝐻(𝑘𝑗,𝜃), known as the Kochin function, describing the elementary waves radiated from the ships is given by 𝐻𝑘𝑗=,𝜃𝑆𝜙𝜕𝜕𝑛𝜕𝜙𝐺𝜕𝑛𝑗(𝜃)𝑑𝑠,(2.13) where 𝐺𝑗𝑘(𝜃)=exp𝑗(𝜃)𝑧+𝑖𝑘𝑗(𝜃)(𝑥cos𝜃+𝑦sin𝜃)(2.14) and 𝑘𝑗(𝜃), 𝑗=1,2, are the unsteady wave numbers: 𝑘𝑗𝐾(𝜃)=0212Ωcos𝜃±14Ωcos𝜃cos2𝜃+𝑗=1𝑗=2.(2.15) When 𝑉=0, then Ω=0, 𝑘1(𝜃) and 𝑘2(𝜃)=𝑘. Thus the wave systems are reduced to the ring wave only, then the drift force in the horizontal plane may be expressed as follows [22]: 𝐹1(2)=𝜌𝑘28𝜋02𝜋||||𝐻(𝜃)2𝐹(cos𝜃cos𝜒)𝑑𝜃;2(2)=𝜌𝑘28𝜋02𝜋||𝐻||(𝜃)2(sin𝜃sin𝜒)𝑑𝜃,(2.16) where 𝐻(𝜃)=𝑆[]𝑛exp𝑘𝑧+𝑖𝑘(𝑥cos𝜃+𝑦sin𝜃)𝜙𝑘𝑧+𝑖cos𝜃𝑛𝑥+𝑖sin𝜃𝑛𝑦𝜕𝜙𝜕𝑛𝑑𝑠.(2.17)

3. Numerical Implementations

For the numerical implementation of the theoretical formulations, the boundary element method (BEM) is adopted. The boundary of the inner domain is discretized by a number of quadrilateral or triangular panels. On each panel the potential value or source strength is assumed to be constant. In the implementation of the developed TDGF method, the Gauss elimination method is used to find the inverse of the influence matrix, while for the employed HYBRID method the generalized minimum residual (GMRES) method is adopted to obtain the numerical solution of the governing integral equations. After the potentials and their spatial partial derivatives on the panels are obtained, the pressure is calculated by Bernoulli’s equation and the hydrodynamic forces are obtained by integrating the pressure over the wetted hull surface. In the following, the 6 DOF motions of the ship are simulated by using an iterative prediction-correction scheme until a convergence is arrived, upon which the simulation will march to the next step. In order to obtain stable and accurate predictions, the Chimera grid concept is introduced [23], in which two panel systems are set up in the beginning of the simulation and information exchange takes place between the two panel systems at every time step.

The flowchart of the numerical procedure and implemented code for HYBRID is shown in Figure 2. Details are given in papers published before (e.g., [15, 20, 23]). In this section some specific considerations are highlighted.

Figure 2: Flowchart of the HYBRID method program, with forward speed.
3.1. Modeling of Initial Transients

In order to smoothly introduce the incident wave disturbance into the numerical scheme and mitigate the effect of initial transients on the steady response to an incident regular wave, the velocity potential of the incident wave is modified as follows: Φ0=𝑔𝜁𝜔𝑒𝑘𝑧[]𝑡sin𝑘(𝑥cos𝜒+𝑦sin𝜒)𝜔𝑡,𝜁=𝜁𝑛𝑇𝑎1cos𝑡𝜋𝜁𝑛𝑇/2𝑡<𝑛𝑇𝑎𝑡𝑛𝑇.(3.1) This is particularly important when simulating horizontal plane motions, especially in short waves.

For steady problem, similar scheme can be applied to the acceleration of the speed and the period should be set as 𝑇=8𝜋𝑈0/𝑔.

3.2. Free-Surface Condition

The free surface can be simulated either by using expression (2.5) or updating the dynamic and kinematic condition of the free surface; respectively, 𝜕Φ𝜕𝑡=𝑔𝜁,𝜕𝜁=𝜕𝑡𝜕Φ.𝜕𝑧(3.2) Though essentially the same, numerically the former one uses all historical data thus increase the memory burden for storing while the latter needs a prediction-correction scheme thus increase computational time. If a higher-order free-surface condition is employed, the partial spatial derivatives on the free-surface panels are needed; they can either be calculated by solving another set of integral equations or by using a finite-difference scheme.

In case of forward speed present, the Chimera grid system is used so as to reduce the computational burden and give accurate description of the near-field free surface. Figure 3 shows an example panelization used in the computation.

Figure 3: Example of free-surface panelization when using Chimera grid system.

4. Results and Discussion

4.1. Diffraction Problem of a Submerged Spheroid with Forward Speed

In the work of Iwashita and Ohkusu [24], an investigation has been carried out on a submerged spheroid with length-to-breadth ratio 𝐿/𝐵=5 and depth-to-breadth ratio 𝑑/𝐵=0.75 with 𝑑 measured from the free surface to the body center. As the added resistance is herein of interest, a fine mesh, especially near the free surface, with 400 panels, is arranged to represent the simple body. Both the TDGF method and the NEWDRIFT code (Papanikolaou [25], a 3D numerical panel method which is based on slender body theory assumptions for the forward speed effects (see [26])) are applied to this case study. First-order numerical results, that is, amplitudes of the wave exciting forces in longitudinal and vertical direction in head seas, are shown in Figure 4. For the 𝑥-direction exciting force, NEWDRIFT and Iwashita’s results appear to agree well, while the present TDGF method result deviates from them in the long-wave range, though the experimental data locate in between them. For the 𝑧-direction force component, a good agreement is observed despite NEWDRIFT’s peak value shift a little from others. This might be due to the fact that while the other two numerical methods take into account forward speed effects exactly, NEWDRIFT takes them into account approximately. Figures 5 and 6 show the numerical results of the added resistance due to the diffraction potential only at different speed. The results of present TDGF method agree well with Iwashita’s results.

Figure 4: Amplitudes of the wave exciting forces in longitudinal and vertical direction in head seas: 𝐹𝑛=0.2.
Figure 5: Diffraction-added resistance: 𝐹𝑛=0.2.
Figure 6: Diffraction-added resistance: 𝐹𝑛=0.3.
4.2. Motion of a Semi-Submersible

This study was conducted for a standard-type semi-submersible body, which has been investigated extensively in previous ITTC studies [27, 28]. In a time domain simulation method the excitation does not necessarily need to be steady, but will be in general an irregular wave excitation; thus, it is important to validate the time histories of the motions. Figure 7 shows the calculated motion histories of the studied ITTC semi-submersible in regular wave 𝑇=3 s, 𝜒=135. The definition of the incident wave is shown in (3.1). Three wave periods prove sufficient to reach a steady-state response after the decay of initial transient effects. Whereas the roll and yaw motions show a steady, periodic response, it is not so with the simulated surge and sway motion due to some slight “drift away” behavior, which may become more significant in short waves. This is actually what we should expect in the simulation of horizontal plane motions in a numerical (or physical) wave basin. Numerically this “drift-away” effect can be removed by applying some light artificial mooring/restoring.

Figure 7: Calculated time histories of surge, sway, roll, and yaw motions; 𝑇=3 s, 𝜒=135.

In this study, about 700 panels are used to represent 1/4 of the whole body. We implement the present hybrid method to examine the quality of the obtained velocity potential with respect to quasi-second-order effects represented by the drift force calculation. Calculated drift forces are significant in short waves, and as in short waves, the drift forces are mainly due to the diffraction effects we present in the following results for the diffraction drift forces in the horizontal plane. As shown in Figure 8, a good agreement between the results from NEWDRIFT and present HYBRID time domain method is observed, despite some deviation in very short waves. This may be due to the differences between the potential solvers (NEWDRIFT in the frequency domain versus HYBRID in the time domains) or drift force formulations (near-field method in NEWDRIFT versus far-field method in HYBRID).

Figure 8: Diffraction drift force in longitudinal and transverse directions; 𝜒=135.
4.3. Motions of S175 Ship at 𝐹𝑛=0.275 and Different Wave Conditions: HYBRID I

The S-175 container ship has been investigated in various benchmark studies by ITTC members since the late 70ties [29, 30]; thus it is still a very good model for validation purposes due to the richness of relevant data. A panelization consisting of nearly 800 panels was prepared to represent half of the ship in simulating the radiation and diffraction problem. The diffraction problem of ITTC S175 hull is studied by the HYBRID I method. Figure 9 shows the wave exciting force/moment results compared to the results from the 3D frequency domain panel code NEWDRIFT. The agreement is overall very good, though the results from the present method are slightly lower than those from NEWDRIFT. It should be noted that the NEWDRIFT code is based on the zero-speed Green function method and forward speed effects are taken into account in an approximate way via slender-body theory assumptions.

Figure 9: Wave exciting force of S175 ship at 𝐹𝑛=0.275.

The HYBRID I method has been also applied to simulate ship motions, under both small-amplitude and large-amplitude assumptions, in which the Froude-Krylov and hydrostatic restoring forces/moments are calculated dependent on the motion model. Figure 10 shows the comparison of results from the present hybrid method against available experimental data and results from NEWDRIFT. For small amplitude incident wave simulation case, the present method gives some results lower than NEWDRIFT and is actually closer to experimental data. For large-amplitude simulations, the incident wave steepness varies systematically; namely, 𝐴/𝜆=0.01,0.02, and 0.04 (noted as CS1, CS2, and CS3, resp., in the graph) where 𝐴 is the wave amplitude. The motion amplitudes decrease gradually as the wave steepness increases. This is physically meaningful, considering the quickly increased damping and restoring due to the above water flared hull form of S175 ship. On the other side, for 𝐴/𝜆=0.04 the resulting peak values of heave and pitch motions are much lower than experimental results, and the steepness of the RAO becomes smaller. For 𝐴/𝜆=0.01 case, in long-wave range the heave motion is very close to experimental data while the pitch motion amplitude is higher. Considering that this wave is quite flat, the deviation of RAOs from results based on small amplitude motion assumption clearly shows the importance of using different models. It should be noted that the wave steepness of the experimental data, which was done for validation of the linear numerical methods then, is not known.

Figure 10: Motion amplitude of S175 ship at 𝐹𝑛=0.275.

Interestingly, taking reference to another source [31], the experimental data of a similar container ship under different wave conditions are revealed (shown in the following as Figures 11(a) and 11(b)). When studying these data, it is observed that as the wave amplitude increases, especially when the amplitude is very large, there is, in general, a trend of RAO’s shift (“bending”) to the longer-wave side and the amplitude decreases, which is similar to the observations in the present study.

Figure 11: (a) The effect of wave height on the response amplitude operator of heave (𝜒=180 deg., 𝐹𝑛=0.239), (b) The effect of wave height on the response ampilitude operator of pitch (𝜒=180 deg., 𝐹𝑛=0.239).

In this case study, the 𝑚𝑗 terms were based on the undisturbed basic flow assumption to account for the steady potential effect on the oscillatory motions.

The added resistance of the ITTC S175 ship in head seas was also calculated based on using the potential and motion data obtained from present hybrid method (noted as HYBRID). Results calculated by the panel code NEWDRIFT (noted as NDfar) and short-wave range corrections based on Faltinsen et al.’s formula [32] and Kuroda et al. [33] and Tsujimoto et al.’s formula [34], respectively, are also shown in Figure 10(c) for comparison (noted as SW1 and SW2, resp.). The calculated results agree well with the experimental data [35] and other numerical results.

4.4. Motions of S175 Ship at 𝐹𝑛=0.275 and Different Wave Conditions: HYBRID II

The motion of the S175 container ship is also studied by using the HYBRID II method. For small-amplitude motion case, we assume the gravity (mass) centre is on the calm waterplane at the midship. For large-amplitude motion case, we estimate the radiation and diffraction forces up to the mean wetted surface but the hydrostatic and Froude-Krylov part exactly up to the instantaneous wavy surface (independent finer mesh is used for this part) about the actual gravity center. The 𝑚𝑗 term is either calculated through solving the relevant integral equations or is calculated in an approximate way, namely, based on the undisturbed basis flow assumption for the steady flow.

The numerical results for the heave and pitch motions are shown in Figure 12. During the numerical simulation, the wave amplitude is set constant 𝐴/𝐿=0.01. When using a small amplitude model, with either way of 𝑚𝑗 term computation, the heave motion is overestimated in long-wave range. However, when using the large amplitude model, which introduces the exact calculation of the Froude-Krylov and restoring forces, the heave motion in the long-wave range is improved and gets closer to experimental data. Some identified problem of the present method appears to be that when compared to experimental data, there is a shift of the peak value of the heave RAO curve for 𝜆 comparable to 𝐿 and this shift cannot be improved by applying the large amplitude model. For the pitch motion, when the small amplitude model is used, we observed also some shift when compared to experimental data. But by applying the large amplitude motion model, the pitch motion results were obviously improved, the shift becomes weak, and the results in all the studied range closely match the experimental results. This shift may be attributed to the 𝑚𝑗 term calculation, with this effect being mainly around the peak range, which is also observed in others’ computations [36]. It is herein confirmed that the 𝑚𝑗 term estimation based on simplified undisturbed basis flow assumption is quite reasonable.

Figure 12: Heave and pitch motion amplitude calculation of S175 ship; 𝐹𝑛=0.275.

It is understood that around the resonance range, the hydrodynamic forces are comparable to the hydrostatic ones and their balance is quite complicated. When we neglect the line-integral term (similar to the last term in integral equation (2.8), but defined in a body-translating system), some error is introduced, which we cannot estimate numerically for the time being. As shown in the calculation of the wave resistance of the Wigley hull [12], despite the highly oscillatory performance which is common to TDGF-based calculation methods, there is a definite deviation from experimental data, which seems to be responsible for the present shift. As pointed out by Lin and Kuang [37], the dissipating far-field boundary condition is no longer suitable for strongly nonlinear cases. Since the ship-generated wave energy is proportional to the Froude number, computational instabilities will occur if the far-field boundaries do not accurately estimate the radiated ship wave energy for high Froude numbers.

An important aspect of HYBRID II method refers to the improvement of the computational efficiency of the HYBRID I method. For forward speed problem, the hull will be moving; thus at every time step the free surface near the ship will change, so that the influence matrices need to be updated at every time step. It takes about 30 seconds on a PC computer for preparing the matrices, depending on how complicated the problem is. After the Chimera grid concept is introduced, though the panels that are far away from the ship are fixed, thus do not need to be updated at every time step, there is still a considerable amount of data processing at each time step. For one simulation, it takes more than 5 hours for the forced motion problem or even more than 10 hours for simulating motion problems on a regular PC hardware with Intel Core 2 QUAD CPU(Q8200 2.33 GHz). As the method is essentially a potential flow solver, this is not satisfactory.

Solving the problem in the body-travelling coordinate system, which results in a panelization system that is quite similar to the zero speed problems, things greatly improve. For a typical motion simulation in head seas condition, it is possible to finish a simulation run within 2 hours. But now the occurring problems are sifted to two other points:(1)the free surface condition update is complicated,(2)the introduced complication with the calculation of the waterline integral term which appears in (2.8).

The second point also affects the first point internally through the solved couple of integral equations.

4.4.1. The Effect of the above Waterplane Hull Form

In previous studies, it was shown that the developed hybrid method is capable of accounting for the effect of above waterplane hull shape changes; thus, it can be practically applied to the optimization of the above water hull shape of ships. For demonstration purpose, we apply here the developed numerical method and computer code to a Wigley hull and two modifications of the basic hull with respect to the above waterplane hull forms, which are defined by the following expressions: 𝑥=𝑥0(1+𝑧)(𝑧0;modied)𝑥=𝑥0𝑏(1+𝑧)(𝑥0,𝑧0;modiedV2)𝑦=12𝑥0𝐿2𝑧1𝐷21+2𝑥0𝐿2/51+2𝑥0𝑧𝐿(𝐷0𝐷)2𝑏𝑧<012𝑥0𝐿21+2𝑥0𝐿2/5𝑧0.(4.1) The modified Wigley hull V2 has a slight flare introduced in the bow/stern regions and an overhang at the stern. The hybrid method is applied to predict the heave and pitch motion of this hull, and results are shown in Figures 13 and 14. Figure 13 shows heave and pitch motions at a small wave steepness; 𝐴/𝐿=0.01. Under this condition, results of applying the large amplitude model are almost identical to the results of the small-amplitude model, as the wave amplitude is small and the flare near the waterplane is also very limited so that the exact calculation of incident wave forces and restoring forces does not affect much the motion.

Figure 13: Heave and pitch amplitude of three different Wigley hulls; 𝐹𝑛=0.2, Amp/𝐿=0.01.
Figure 14: Heave and pitch amplitude of three different Wigley hulls; 𝐹𝑛=0.2, Amp/𝐿=0.02.

Under a steeper wave condition, 𝐴/𝐿=0.02, more deviations show up for both heave and pitch motions, as shown in Figure 14. For the heave motion, the amplitude of the modified hull has decreased compared to the other two hulls due to its large projected area on the vertical direction. On the other hand, the pitch motion of the modified hull V2 is obviously smaller than for the other two hulls in the long-wave range. These calculations offer some valuable information regarding the applicability of these methods in the preliminary design stage. When considering the actual sea conditions encountered in ship’s service route, we may be able to determine the optimal route or even the optimal hull shape for specific routes.

We should keep in mind that all the three Wigley hulls are actually very narrow compared to actual shiplike forms, with 𝐿/𝐵=10, and the introduced flares are also quite conservative compared to a real ship. Due to this reason, the 𝑚𝑗 term based on undisturbed basis flow assumption is employed.

5. Summary and Conclusions

In this paper, the formulation and validation of a newly developed time-domain transient Green function method and a hybrid method are presented. The developed methods are applied to different types of bodies (mathematical bodies, offshore platform, ship hull forms) under different operational conditions. Numerical results cover the radiation/diffraction problem, 6 DOF motions, either of small amplitude or large amplitude, either with zero speed or nonzero speed forward speed, and drift force/added resistance in waves. Through these (and similar others shown in the listed references) validation studies, the developed methods/codes proved to be valuable design tools for the ship hull-form assessment and optimization, as they are capable of studying the seakeeping characteristics of hulls with different above-water shapes and for a variety of forward speed and incident wave conditions.

The way ahead of the development work of our system is to account for nonlinearities of the radiation/diffraction disturbance on the free surface, to explore alternative numerical concepts for the matching of the free surface conditions and to further improve the efficiency of the numerical implementation of the methods by parallel programming in the new computer—cluster environment of NTUA-SDL.


The present work covers the first author’s Ph.D. and part of his postdoc research; the financial support of NTUA-SDL and continuous guidance by Professors A. Papanikolaou and W. Duan (Harbin Engineering University) is acknowledged. The authors also thank Dr. Ogawa for the permission to cite his research results in the current work.


  1. A. B. Finkelstein, “The initial value problem for transient water waves,” Communications on Pure and Applied Mathematics, vol. 10, pp. 511–522, 1957. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  2. W. E. Cummins, “The impulsive response function and ship motions,” Schiffstechnik, vol. 9, pp. 124–135, 1962. View at Google Scholar
  3. W. M. Lin and D. Yue, “Numerical solutions for large-amplitude ship motions in the time-domain,” in Proceedings of the 18th Symposium on Naval Hydrodynamics, pp. 41–66, 1990.
  4. S. P. Singh, D. Sen, and D. G. Sarangdhar, “3D seakeeping computation for different hull forms for design evaluations,” in Proceedings of the 15th International Conference on Hydrodynamics in Ship Design, Safety and Operation, pp. 245–256, Gdnask, Poland, 2003.
  5. W. Duan and Y. Dai, “Time-domain calculation of hydrodynamic forces on ships with large flare,” International Shipbuilding Progress, vol. 46, no. 446, pp. 223–232, 1999. View at Google Scholar · View at Scopus
  6. S. Zhang, “A hybrid boundary-element method for non-wall-sided bodies with or without forward speed,” in Proceedings of the 13th International Workshop on Water Waves and Floating Bodies (IWWWFB'98), pp. 179–182, 1998.
  7. H. Yasukawa, “Application of a 3-D time domain panel method to ship seakeeping problems,” in Proceedings of the 24th Symposium on Naval Hydrodynamics, pp. 376–392, 2003.
  8. S. Kataoka and H. Iwashita, “Estimations of hydrodynamic forces acting on ships advancing in the calm water and waves by a time-domain hybrid method,” Journal of the Society of Naval Architects of Japan, vol. 196, pp. 123–138, 2004. View at Google Scholar
  9. A. Iafrati and E. F. Campana, “A domain decomposition approach to compute wave breaking (wave-breaking flows),” International Journal for Numerical Methods in Fluids, vol. 41, no. 4, pp. 419–445, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  10. M. Sueyoshi, H. Kihara, and M. Kashiwagi, “A hybrid technique using particle and boundary-element methods for wave-body interaction problems,” in Proceedings of the 9th International Conference on Numerical Ship Hydrodynamics, vol. 1, pp. 241–252, 2007.
  11. W. M. Lin, H. Chen, and S. Zhang, “A hybrid numerical method for wet deck slamming on a high speed catamaran,” in Proceedings of the 10th International Conference on Fast Sea Transportation, 2009.
  12. S. K. Liu, Numerical simulation of large amplitude ship motions and applications to ship design and safe operation [Ph.D. thesis], Ship Design Lab., National Technical University of Athens, 2011.
  13. S. Liu, A. Papanikolaou, and W. Duan, “A time domain numerical simulation method for nonlinear ship motions,” Journal of Harbin Engineering University, vol. 27, no. 2, pp. 177–185, 2006. View at Google Scholar · View at Scopus
  14. S. K. Liu and A. Papanikolaou, “A time-domain hybrid method for calculating hydrodynamic forces on ships in waves,” in Proceedings of the 13th International Maritime Association of the Mediterranean (IMAM'09), Istanbul, Turkey, 2009.
  15. S. K. Liu and A. Papanikolaou, “Time-domain hybrid method for simulating large amplitude motions of ships advancing in waves,” Journal of the Korean Society of Naval Architects and Ocean Engineers, vol. 3, no. 1, pp. 72–79, 2011, Proceeding of the ITTC Benchmark Study and Workshop on Seakeeping, 2010. View at Publisher · View at Google Scholar
  16. S. Liu, A. Papanikolaou, and G. Zaraphonitis, “Prediction of added resistance of ships in waves,” Ocean Engineering, vol. 38, pp. 641–650, 2011. View at Google Scholar
  17. G. Papatzanakis, A. Papanikolaou, and S. K. Liu, “Optimization of routing with uncertainties,” in Proceedings of the 14th International Maritime Association of the Mediterranean (IMAM'11), Genoa, Italy, 2011.
  18. J. O. Kat de, “The numerical modeling of ship motions and capsizing in severe seas,” Journal of Ship Research, vol. 34, no. 4, pp. 289–301, 1990. View at Google Scholar
  19. S. K. Liu and A. Papanikolaou, “Application of a nonlinear time domain hybrid method to the study of a semi-submersible in waves,” in Proceedings of the 22nd International Offshore (Ocean) and Polar Engineering Conference, Greece, 2012.
  20. S. K. Liu and A. Papanikolaou, “A time-domain hybrid method for calculating hydrodynamic forces on ships in waves,” in Proceedings of the 14th International Maritime Association of the Mediterranean (IMAM'11), Genoa, Italy, 2011.
  21. H. Maruo, “Resistance in waves,” in 60th Anniversary Series, vol. 8, pp. 67–102, The Society of Naval Architects of Japan, 1963. View at Google Scholar
  22. H. Maruo, “The drift of a body floating on waves,” Journal of Ship Research, vol. 4, no. 3, pp. 1–10, 1960. View at Google Scholar
  23. S. K. Liu and A. Papanikolaou, “Application of Chimera grid concept to simulation of the free-surface boundary condition,” in Proceedings of the 26th International Workshop on Water Waves and Floating Bodies (IWWWFB'11), 2011.
  24. H. Iwashita and M. Ohkusu, “The Green function method for ship motions at forward speed,” Journal of Ship Technology Research, vol. 39, no. 2, 1992. View at Google Scholar
  25. A. Papanikolaou, “On integral-equation-methods for the evaluation of motions and loads of arbitrary bodies in waves,” Ingenieur-Archiv, vol. 55, no. 1, pp. 17–29, 1985. View at Publisher · View at Google Scholar · View at Scopus
  26. N. Salvesen, E. O. Tuck, and O. Faltinsen, “Ship motions and sea loads,” Transactions of the SNAME, vol. 78, pp. 250–287, 1970. View at Google Scholar
  27. ITTC Ocean Engineering Committee, Report of Ocean Engineering Committee, 1984.
  28. ITTC Ocean Engineering Committee, Report of Ocean Engineering Committee, 1987.
  29. ITTC Seakeeping Committee, “Comparison of results obtained with compute programs to predict ship motions in six-degrees-of-freedom and associated responses,” in Proceedings of the 15th International Towing Tank Conference (ITTC'78), pp. 79–92, 1978.
  30. Y. Kim, “Comparative study on linear and nonlinear ship motion and loads,” in Proceedings of the ITTC Workshop on Seakeeping, pp. 283–257, 2010.
  31. Y. Ogawa, “A study on nonlinear wave loads of a large container carrier in rough seas,” in Proceedings of the 10th International Symposium on Practical Design of Ships and Other Floating Structures, 2007.
  32. O. M. Faltinsen, K. J. Minsaas, N. Liapis, and S. O. Skjørdal, “Prediction of resistance and propulsion of a ship in a seaway,” in Proceedings of the 13th Symposium on Naval Hydrodynamics, pp. 505–529, 1980.
  33. M. Kuroda, M. Tsujimoto, and T. Fujiwara, “Investigation on components of added resistance in short waves,” Journal of the Japan Society of Naval Architects and Ocean Engineers, vol. 8, pp. 171–176, 2008. View at Google Scholar
  34. M. Tsujimoto, K. Shibata, and M. Kuroda, “A practical correction method for added resistance in waves,” Journal of the Japan Society of Naval Architects and Ocean Engineers, vol. 8, pp. 147–154, 2008. View at Google Scholar
  35. T. Takahashi, “A practical prediction method for added resistance of a ship in waves and the direction of its application to hull form design,” Transactions of the West-Japan Society of Naval Architects, vol. 75, pp. 75–95, 1988. View at Google Scholar
  36. H. B. Bingham and H. D. Maniar, “Computing the double-body m-terms using a B-spline based panel method,” in Proceedings of the 11th International Workshop on Water Waves and Floating Bodies (IWWWFB'96), Institut Für Schiffbau, Hamburg, Germany, 1996.
  37. R.-Q. Lin and W. J. Kuang, “A nonlinear method for predicting motions of fast ships,” in Proceedings of the 10th International Conference on Fast Sea Transportation, Athens, Greece, October 2009.