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

Abstract

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 𝑈₀ 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: Φ_𝑇(𝑥,𝑦,𝑧,𝑡)=Φ₀(𝑥,𝑦,𝑧,𝑡)+Φ(𝑥,𝑦,𝑧,𝑡),(2.2) where Φ₀ 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: 𝜕Φ∇²Φ(𝑝,𝑡)=0,(2.3)𝜕𝑛=𝑣^𝑛−𝜕Φ⁰𝜕𝑛on𝑆^𝑏𝜕(𝑡),𝑡>0,(2.4)₂Φ𝜕𝑡²+𝑔𝜕Φ𝜕𝑧=0on𝑆^𝑓Φ=𝜕Φ(𝑡),𝑡>0,(2.5)𝜕𝑡=0on𝑆^𝑓Φ,∇Φ,𝜕Φ(𝑡),𝑡=0,(2.6)𝜕𝑡→0on𝑆^∞,𝑡>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 𝐺(𝑝,𝑡;𝑞,𝜏)=2^∞0√𝑔𝑘𝑒^{𝑘(𝑧+𝜁)}𝐽₀(𝑘𝑅)sin√𝑔𝑘(𝑡−𝜏)𝑑𝑘,(2.9) where 𝑝(𝑥,𝑦,𝑧) and 𝑞(𝜉,𝜂,𝜁) are the field and source points; respectively, 𝐫_𝑝𝑞=(𝑥−𝜉)𝐢+(𝑦−𝜂)𝐣+(𝑧−𝜁)𝐤, and 𝐽₀ 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 −𝑈₀𝑥 [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+^𝛼𝜋/2₀−^{3𝜋/2𝜋/2}||𝐻𝑘¹,𝜃||²𝑘₁𝑘₁cos𝜃−𝑘cos𝜒+𝜌√1−4Ωcos𝜃𝑑𝜃8𝜋^{2𝜋−𝛼0}_𝛼0||𝐻𝑘₂,𝜃||²𝑘₂𝑘₂cos𝜃−𝑘cos𝜒√1−4Ωcos𝜃𝑑𝜃,(2.12) where 𝜒 is the incident wave heading, 𝜃 is the elementary wave angle, 𝜌 is the density of sea water, and 𝛼₀ is the critical angle (𝛼₀=arcos(1/(4Ω)) for Ω>1/4 and 𝛼₀=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: 𝑘_𝑗(𝜃)=𝐾⁰21−2Ωcos𝜃±√1−4Ωcos𝜃cos²𝜃⎛⎜⎜⎝+−∶∶𝑗𝑗==12⎞⎟⎟⎠.(2.15) When 𝑉=0, then Ω=0, 𝑘₁(𝜃)→∞ and 𝑘₂(𝜃)=𝑘. 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]: 𝐹¹⁽²⁾=𝜌𝑘²8𝜋^02𝜋𝐹||𝐻(𝜃)||²(cos𝜃−cos𝜒)𝑑𝜃;₂₍₂₎=𝜌𝑘²8𝜋^02𝜋||𝐻(𝜃)||²(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: Φ₀=𝑔𝜁𝜔𝑒^𝑘𝑧sin[𝑘(𝑥cos𝜒+𝑦sin𝜒)−𝜔𝑡],𝜁=𝜁𝑛𝑡𝑇𝜁^𝑎1−cos𝑛𝑡𝜋𝑇/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𝜋𝑈₀/𝑔.

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.

(a)

(b)

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.

(a)

(b)

(c)

(d)

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).

(a)

(b)

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.

(a)

(b)

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.

(a) Heave amplitude of S175 ship at 𝐹𝑛=0.275

(b) Pitch amplitude of S175 ship at 𝐹𝑛=0.275

(c) Added resistance of S-175 ship at 𝐹𝑛=0.275

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.

(a)

(b)

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.

(a)

(b)

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: 𝑥=𝑥𝑥₀=𝑥₀(1+𝑧)(𝑧≥0;modified)𝑦=⎧⎪⎪⎨⎪⎪⎩𝑏1−2𝑥₀(1+𝑧)(𝑥≤0,𝑧≥0;modifiedV2)𝐿²1−𝐷𝑧²1+2𝐿𝑥⁰²/51+2𝑥⁰𝑧𝑏1−2𝑥₀𝐿(𝐷₀−𝐷)²𝑧<0𝐿²1+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.

(a)

(b)

Figure 13 

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

(a)

(b)

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.

Acknowledgments

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.