A Fully Nonlinear, Dynamically Consistent Numerical Model for Ship Maneuvering in a Seaway

Abstract

This is the continuation of our research on development of a fully nonlinear, dynamically consistent, numerical ship motion model (DiSSEL). In this paper we report our results on modeling ship maneuvering in arbitrary seaway that is one of the most challenging and important problems in seakeeping. In our modeling, we developed an adaptive algorithm to maintain dynamical balances numerically as the encounter frequencies (the wave frequencies as measured on the ship) varying with the ship maneuvering state. The key of this new algorithm is to evaluate the encounter frequency variation differently in the physical domain and in the frequency domain, thus effectively eliminating possible numerical dynamical imbalances. We have tested this algorithm with several well-documented maneuvering experiments, and our results agree very well with experimental data. In particular, the numerical time series of roll and pitch motions and the numerical ship tracks (i.e., surge, sway, and yaw) are nearly identical to those of experiments.

1. Introduction

Ship motion with fixed heading is an overly simplification of general ship motion dynamics, for example, ship maneuvering in seaway is far more complex than the fixed heading. In our earlier work, Lin and Kuang [1] developed a fully nonlinear, dynamically consistent ship motion model (DiSSEL) that can accurately simulate ship motion with fixed heading. A continuing effort is therefore to expand DiSSEL for simulating ship motions when maneuvering in seaway.

Accurate prediction of ship motion when maneuvering in a seaway is one of the most challenging and important problems in ship motion dynamics. Much effort has been devoted in the past on modeling ship maneuvering in seaway. The earliest studies started from the simplest cases of simulation ship tracks (i.e., surge, sway, and yaw motions) in calm water. For example, Chaplin [2] and Folger et al. [3] predicted ship tracks in calm water based on empirical laws derived from experimental data. Söding [4] made a substantial step forward by including dynamics in ship motion modeling. He introduced a set of hydrodynamic maneuvering coefficients in simulation, though those coefficients are determined mainly from experimental data. Hooft and Pieffers [5] furthered this effort by including the effects of waves on ship motion (e.g., wave induced motions). In 1988, Lee [6] developed a potential flow method to calculate the rudder forces and the moments, thus making it possible to fully simulate ship maneuvering motion state. Nearly a decade late, Bailey et al. [7] proposed a hybrid approach in which experimental data and theoretical formulations (e.g., linearized or parameterized ship response function) are integrated into a single model, thus providing an enhanced description of fluid loading. More recently, several complicated computational models are developed for simulating complex cases, for example, the viscous flow models by Chau [8] and El Moctar [9].

These studies have greatly improved our understandings of ship tracks and ship maneuvering dynamics. The methodologies developed from these efforts are still widely used [10–15]. However, as pointed by Söding [10], there are also many limitations in these methods. For example, the potential flow methods ignore the dissipative effects from (tiny) fluid viscosity, small-scale turbulence, and flow separation. The viscous flow models are capable of evaluating these dissipative effects. But they are computationally difficult and expensive: numerical instabilities and computational capability limitations make direct simulation in the oceanic parameter domain impossible. In addition, inclusion of any a priori condition (e.g., those empirical extrapolations from experimental data) will affect generality of the model.

To utilize the mathematical simplicity of potential flow approaches and to include appropriate dissipative effects, Lin and Kuang [16] developed a physically more realistic and computationally more efficient model in which the “blocking theory” is used to numerically determine the rudder forces, the lift forces, the moments, and the Kutta conditions at the trailing edges. Based on this approach, Lin et al. [17] were able to model accurately ship track and heel motion when maneuvering in calm water. Their results are well benchmarked by experiments data.

The natural extension of the effort by Lin and Kuang [16] and by Lin et al. [17] is to develop a dynamically consistent model for predicting ship maneuvering in seaway. The model by Lin and Kuang [1] can be a very good candidate for this purpose for several reasons: there is no a priori condition nor parameterization in this model, and the model is fully nonlinear and dynamically consistent. In particular, prediction of ship motion in fixed heading from this model agrees well with experiments. For predicting ship maneuvering in seaway, several mechanisms need to be included (or reintroduced) into this model, including change of ship heading, motion acceleration, and propeller forces.

But, one might anticipate a mathematical difficulty of modeling encounter waves in ship maneuvering. The wave frequency and the wave number vector of the encounter waves vary with the ship speed and ship heading. The variation due to the latter is in particular significant in ship maneuvering. As we will find later in this paper, changes in the encounter frequencies do have profound impacts on ship-motion interactions and thus on the ship motion state. In this paper, we will discuss the details of this model expansion, the mathematical/numerical difficulties, and the algorithms resolving these difficulties. This paper is organized as follows: in Section 2 is an overview of DiSSEL. Section 3 describes the details of the relevant algorithm. Model results validation is given in Section 4, followed by the discussions.

2. Fundamentals of DiSSEL Nonlinear Ship Motion Model

DiSSEL is a hybrid-flow ship motion model in time domain: the surface waves are described by velocity potential, but the fluid viscous effect is included in the form of wave breaking and wave blocking (separation effects) [1, 16, 18]. Inclusion of the viscous effect makes DiSSEL satisfying the Kutta condition at a trailing edge, and velocity stagnation points (𝐕=0), thus allowing accurate evaluation of the lift force [17]. This hybrid-flow approach builds a solid foundation of our effort to expand the model for ship maneuvering in a seaway.

In DiSSEL, ship-wave interaction and ship solid body motion dynamics are modeled by two separate components. Interactions between the two dynamical processes are determined through the pressure force on the ship hull surface and the restoring force from the water displaced by underwater ship volume. For the details of the partial differential equations, the boundary conditions and the surface integrations, we refer the reader to Lin et al. [19], Lin and Kuang [1, 20].

To expand DiSSEL for a ship maneuvering in a seaway, it is necessary to include a propulsion mechanism. As the first step, this mechanism can be approximated from experimental data (empirical). This approximation does not affect the dynamical consistency of our model, since the propulsion is determined by the ship engine, independent of the ship-wave interactions. In the following, we focus on the differences in the equations and the physical variables that are different from those of the original DiSSEL model.

Another notice should be made on the definitions of the reference frames used in DiSSEL. In Lin and Kuang [1], two additional reference frames are introduced: the ship reference frame (SRF) and the model reference frame (MRF). The latter is defined to be the reference frame moving horizontally with the ship mass center, but with the origin on the mean free surface. When a ship maneuvers in seaway, its track direction can change. But MRF does not change its orientation. Therefore, the rotational effect from this direction change needs to be accounted for correctly. For details of the reference frames, we refer the reader to Lin and Kuang [1].

2.1. Ship-Wave Interaction Model

In DiSSEL, the surface waves are described by the following set of partial differential equations: ∇²𝜑≡∇^2ℎ𝜑+𝜕²𝜑𝜕𝜑𝜕²𝑧=0,(−𝐻≤𝑧≤𝜂),𝜕𝑡+12∇𝜑−𝐮^𝑠+𝐯_ℎ⋅∇𝜑+𝑝𝜌_𝑤−𝜕𝐮+𝑔𝑧_𝑠+𝐯_ℎ𝜕𝑡⋅𝐱−𝜈∇^2ℎ𝜕𝜂𝜑=0,(−𝐻≤𝑧≤𝜂),𝜕𝑡+∇^ℎ𝜂⋅∇_ℎ𝜑−𝐮_𝑠+𝐯_ℎ=𝜕𝜑𝜕𝑧,(𝑧=𝜂),(1) where 𝑝 is the pressure, 𝜌_𝑤 is the water density, 𝜈 is the turbulent viscosity (arising from the wave-breaking processes), 𝐱 is the position vector, 𝐻 is the water depth, 𝜑 is the velocity potential, 𝜂 is the free surface elevation, ∇(≡𝐱̂𝜕/𝜕𝑥+𝐲̂𝜕/𝜕𝑦+𝐳̂𝜕/𝜕𝑧) is the gradient, ∇_ℎ is the horizontal component part of ∇, and 𝐮_𝑠 is the ship maneuvering velocity (in SRF, it is the given forward speed). It should be noted that 𝜕𝐮_𝑠/𝜕𝑡 includes both magnitude and direction changes of the maneuvering velocity. In addition, 𝐯_ℎ is the horizontal component of translational velocity 𝐯_𝑐 for surge and sway (absent in the fixed heading scenario).

The impenetrable boundary conditions are𝜕𝜑𝜕𝑧=0(2) at the flat bottom 𝑧=−𝐻, and 𝐧̂_Σ⋅∇𝜑=𝐧̂_Σ⋅𝐮_𝑠+𝐯_𝑐+Ω×𝐱^Σ−𝐱_𝑐(3) on the ship surface 𝐱^Σ. In the above equation, 𝐧̂_Σ is the normal unit vector of the ship surface 𝐱^Σ, and 𝐱_𝑐 is the mass center of the ship hull.

2.2. Solid Body Motion Model

The ship solid body motion is decoupled into the translational motion of ship mass center 𝐱_𝑐 and the rotation about 𝐱_𝑐. The translational motion, 𝐯_𝑐, is described by the following equation in MRF: 𝑚_𝑠𝑑𝐯_𝑐𝑑𝑡+𝐷^𝜈𝐯_𝑐=𝐅,(4a) where 𝑚_𝑠 is the total ship mass and 𝐅 is the total force on the ship: 𝐅=𝐅^𝑝+𝐅^𝑔+𝐅^𝑟,(4b) where 𝐅^𝑝≡−_Σwet𝐧̂_Σ𝑝𝑑Σ(4c) is the (driving) pressure force on the ship hull,𝐅^𝑔=𝜌_𝑤𝑉_wet−𝑚_𝑠𝐠,(4d)

is the buoyancy force (restoring force). This force acts only in the vertical direction. But the rudder force 𝐅^𝑟 in (4b) is the new addition. We will discuss this in more detail in the next section. 𝑉_wet is under water volume.

The solid body rotation about the ship mass center is governed by the Liouville equation defined in SRF: 𝐈⋅𝑑𝛀𝑑𝑡+𝛀×(𝐈⋅𝛀)+𝐷^Ω(𝐈⋅𝛀)=𝚪,(5a) where 𝐈 is the moment of inertia tensor about the ship mass center 𝐱_𝑐 (and its components are denoted by 𝐼_(𝑖,𝑗)), 𝐷_𝜈 and 𝐷_Ω are integrated dissipative effects (e.g., drag) of the fluid on the ship hull, and Γ is the total moment (torque) on the ship:𝚪=𝚪^𝑝+𝚪^𝑔+𝚪^𝑟.(5b) It includes the (driving) pressure torque:𝚪^𝑝=−_Σwet𝑑Σ𝐱^Σ−𝐱_𝑐×𝐧̂_Σ𝑝,(5c) (where Σ_wet is total wetted ship surface), the restoring torque: 𝚪^𝑔=𝐱_wet−𝐱_𝑐×𝑚_𝑠−𝜌_𝑤𝑉_wet𝐠,(5d)

(where 𝐱_wet is the position vector in under water volum), and the torque Γ^𝑟 from the rudder (a new addition). It should be pointed out that the restoring torque is always horizontal, thus affecting only roll and pitch motions.

2.3. Rudder Force and Moment

The rudder force 𝐅^𝑟 in (4b) and the moment Γ^𝑟 in (5b) are the new physical variables included in this expansion. Since the rudder position is often defined relative to the ship hull, it is therefore very convenient to evaluate them in SRF. Therefore, this leads to the following formulations: 𝚪𝐅^𝑟=_𝐴𝐧̂^𝑟𝑝^𝑟𝑑𝑠=𝐹^𝑟𝑥𝐱̂_𝑠+𝐹^𝑟𝑦𝐲̂_𝑠+𝐹^𝑟𝑧𝐳̂_𝑠,(6)_𝑟=_𝐴𝐑×𝐧̂^𝑟𝑝^𝑟𝑑𝑠+𝐑^𝑡×𝐅^Ω=Γ^𝑟𝑥𝐱̂_𝑠+Γ^𝑟𝑦𝐲̂_𝑠+Γ^𝑟𝑧𝐳̂_𝑠+𝐱−𝐱_𝐶×𝐅^Ω,(7) where 𝑝^𝑟 is the potential flow pressure; but 𝐴=𝐴^Σ−𝐴^block (𝐴^Σ is the wetted rudder surface, and 𝐴^block is the blocking effected area) [16], 𝐧̂^𝑟 is the unit normal vector of 𝐴, 𝐱_𝐶 is the position vector of the center of the circle tangent to the point 𝐱 on the ship trajectory, 𝐑^𝑡 is the radius vector of ship trajectory to the point 𝐱, and𝐑≡𝑅_𝑥𝐱̂_𝑠+𝑅_𝑦𝐲̂_𝑠+𝑅_𝑧𝐳̂_𝑠(8) is the rotating arm vector of the rudder. The force 𝐅^Ω in (7) describes the effect of maneuvering directional changes in SRF and is, therefore, dependant on the changes in the rudder angle 𝜹 (defined relative to the ship hull):𝐅^Ω=_𝐴𝜌2𝛀^∗×𝐯_𝑐+𝛀^∗×𝛀^∗×𝐱+𝑑𝛀^∗𝛀𝑑𝑡×𝐱𝑑𝑠,_∗=𝑑𝜷𝑑𝑡+Ω^𝑧𝐳̂_𝑠,(9) where 𝜷 is the drifting angle vector (the angle between the ship’s horizontal velocity vector and the ship’s longitudinal axis). If the rudder can be approximated as a plate rotating about the vertical direction of the ship hull with a rudder deflection angle 𝛿, the horizontal components of 𝐅^𝑟 can be substantially simplified:𝐹^𝑟𝑥=−𝐶_𝑙12𝜌𝑢^2𝑟𝐹𝐴sin²𝛿+_𝐴𝑛^𝑟𝑥̂𝑝^𝑟𝑤𝑑𝑠,_𝑟𝑦=𝐶_𝑙14𝜌𝑢^2𝑟𝐴sin2𝛿+_𝐴𝑛^𝑟𝑦̂𝑝^𝑟𝑤𝑑𝑠,(10) where 𝑢_𝑟 is the effective inflow speed; 𝑛^𝑟𝑥̂ and 𝑛^𝑟𝑦̂ are the 𝑥 and 𝑦 components of the normal vector of the rudder surface, respectively; 𝑝^𝑟𝑤 is pressure on rudder surface due to the incident waves; 𝐶_𝑙 is the lift curve slope:𝐶_𝑙=2𝜋1+(2/(𝐿/𝑑)),(11) where 𝐿 and 𝑑 are the length and the width of the rudder, respectively [21]. This theory has been tested by many experiments and is broadly used in ship dynamics. It should be pointed out that the separation effects are included by using 𝐴 (not 𝐴^Σ) in (6)–(10).

By (10), the vertical moment is Γ^𝑟𝑧=𝑅_𝑥𝐹^𝑟𝑦−𝑅_𝑦𝐹^𝑟𝑥≈𝑅_𝑥𝐹^𝑟𝑦=𝐶_l14𝜌𝑢^2𝑟𝐴sin2𝛿+_𝐴𝑛^𝑟𝑦̂𝑝^𝑟𝑤𝑑𝑠𝑥_𝑟−𝑥_𝑐,(12) since 𝑅_𝑦 is in general much smaller than 𝑅_𝑥. In (12), (𝑥_𝑟−𝑥_𝑐) is the longitudinal distance between the ship mass center 𝑥_𝑐 and the rudder mass center 𝑥_𝑟.

The effective inflow into the rudder is derived from the ship forward speed and the accelerated flow in the propeller slipstream. The latter depends on the propeller’s mechanical design. In the cases discussed in this paper, we choose an empirical estimate 𝐮_𝑟=1.6𝐮_𝑠,(13) from experimental data for naval vessels with Froude numbers ranging from 0.2 to 0.4 [22].

2.4. Ship Track Kinematics

Kinematics of ship track variation can be best described in an inertial reference frame, for example, the earth-fixed reference frame (ERF) (𝐱̂_𝑒,𝐲̂_𝑒,𝐳̂_𝑒). The ship moving direction angle, 𝜃^∗, is the sum of the yaw angle 𝜃_𝑧 (the rotational angle in 𝑧 direction, and 𝑑𝜃_𝑧/𝑑𝑡=Ω_𝑧 obtained from (5a)–(5d)) and the drift angle 𝛽:𝜃𝛽^∗=≡t𝜃a_𝑧n+₋₁𝛽,||𝐮_𝑠||𝜈+^𝑦𝜈_𝑥,(14) where 𝜈_𝑥 and 𝜈_𝑦 are ship translation motion in the 𝑥 and 𝑦 directions in ship coordinate, respectively. The sketch map can be find in Figure 5 page 198 [23, Figure 5, page 198]. The track 𝐱_track of the ship in ERF is defined as 𝐱_track=𝑥_track𝐱̂_𝑒+𝑦_track𝐲̂_𝑒.(15) In particular,𝑑𝐱_track𝑑𝑡=𝐀^𝑚⋅𝐮_𝑠+𝐯_ℎ≡𝐮,𝐀_𝑚=cos𝜃^𝑧−sin𝜃_zsin𝜃_𝑧cos𝜃_z.(16)

If the rudder deflection angle 𝛿 is given, then the rudder force (10) and moment (12) can be determined. If the ship velocity is calculated, for example, via (4a)–(4d) and (5a)–(5d), then the track can be determined via time integration of (16). On the other hand, given the navigation history, that is, the ship track 𝐱_track(𝑡), the total ship horizontal velocity, 𝐮, and horizontal acceleration can be calculated by simple time differentiation, and the direction angle 𝜃^∗ can be also determined via𝜃^∗=tan⁻¹𝐮𝐮⋅⋅𝐲𝐱̂̂^𝑒_𝑒.(17) The the rudder angle 𝛿 can be obtained from horizontal component of (4a)–(4d), vertical component of (5a)–(5d), (10), and (16) using an iterative method. Finally we want to point out that, in DiSSEL ship motion model, hydrodynamic interactions between the hull and the rudder are properly implemented, in addition to the implementation of the interactions of the six-degree freedom of the ship motions. These nonlinear interactions are obtained by solving a set of the nonlinear equations (1)–(7).

2.5. Centrifugal Effect

Ship track direction changes generate a finite centrifugal force. In a very simple mathematical description, the force can be described as a vertical rotation about the arc center 𝐶 of the track curve: 𝐮_𝑠=𝜔𝐳̂_𝑒×𝐱_track−𝐱_𝐶,(18) where 𝜔 is the rotation speed about 𝐶 and 𝐱_𝐶 is the position vector of the point 𝐶 in ERF. Therefore,𝜕𝐮_𝑠𝜕𝑡=𝜔²𝐳̂_𝑒×𝐳̂_𝑒×𝐱_track−𝐱_𝐶.(19) Thus,𝑝𝜌+𝜕𝜕𝐮𝑡^𝑠⋅𝐱^track=1𝜌𝑝+12𝜌𝜔²||𝐳̂_𝑒×𝐱_track−𝐱_𝐶||²≡𝜌𝑝̃.(20) In other words, the changes in the direction of the ship motion can be described by the effective pressure 𝑝̃ in (20). Given the ship track 𝐱_track(𝑡), the radius |𝐱_track−𝐱_𝐶| can then be calculated with the standard curvature formulation: ||𝐱_track−𝐱_𝐶||=1𝜅(𝑡),𝜅(𝑡)=||𝑑𝑥^track/𝑑𝑡𝑑²𝑦_track/𝑑𝑡²−𝑑𝑦_track/𝑑𝑡𝑑²𝑥_track/𝑑𝑡²||𝑑𝑥^track/𝑑𝑡²+𝑑𝑦_track/𝑑𝑡²^3/2.(21) Its influence on the solid body rotation of the ship can be computed via integration of the effective pressure over the ship surface.

3. Stable and Conservative Algorithm for Ship Maneuvering

When a ship maneuvers in seaway, its direction and speed change over time. This leads to time-varying encounter waves, thus potentially to numerical difficulties in conserving the momentum of the entire system. This will be a source of numerical instabilities. Therefore, we will discuss first the nature of the numerical difficulties (based on current ship motion models) and then introduce our new stable algorithm that conserves the momentum of the system.

3.1. Dynamical Imbalance from Rapidly Changing Encounter Frequencies

As pointed out by Lin and Kuang [24], the numerical balance between the driving forces and the restoring forces is the key to model correctly ship motions. Lin and Kuang [1] further showed that DiSSEL can well maintain such dynamic balances for fixed heading of various ship hulls in a wide range of environments.

To expand DiSSEL for modeling ship maneuvering, a natural consideration is whether the algorithm used in the model is able to maintain accurately the dynamic balance between the driving forces and the restoring forces on the ship hull when ship heading changes in time. The consideration arises from the very fact that such changes also alter the environment observed on the ship. For example, in SRF, a surface wave can be generally described as 𝜂(𝑥,𝑦,𝑡)=^𝑁_𝑀𝑖=1_𝑗=1𝐴_𝑖𝑗sin𝐤_𝑖𝑗⋅𝐱−𝐱_𝑐+𝜔^en_𝑖𝑗𝑡+𝜙_𝑖𝑗,(22) where 𝐴_𝑖𝑗 is the amplitude of the mode (𝑖,𝑗), 𝐤_𝑖𝑗 is the wave number vector, 𝜔^en_𝑖𝑗 is the encounter frequency, and 𝜙_𝑖𝑗 is the phase. The integer pair (𝑀,𝑁) defines the total number of the wave modes used in numerical simulation. The encounter frequency 𝜔^en_𝑖𝑗 is defined as 𝜔^en_𝑖𝑗=𝜔^sea_𝑖𝑗1−||||𝐜𝐮_𝑖𝑗||||cos𝛼_sea𝑖𝑗,(23) where 𝐜_𝑖𝑗 is the phase velocity of the wave mode, 𝛼^sea_𝑖𝑗 is the angle between 𝐜_𝑖𝑗 and 𝐮, and 𝜔^sea_𝑖𝑗 is the incident wave frequency in ERF. For any given incident wave, 𝛼^sea_𝑖𝑗 changes with the ship heading direction.

One could clearly observe the complications and its numerical difficulties from (23): changes in 𝐯_ℎ and 𝜃_𝑧 (and therefore 𝛼^sea_𝑖𝑗) are dynamically determined via the driving forces and the restoring forces on the ship hull which depend on the encounter frequency 𝜔^en_𝑖𝑗. But changes in 𝜔^en_𝑖𝑗 depend also on 𝐯_ℎ and 𝜃_𝑧. Therefore, errors in the force balances could potentially lead to numerical instabilities and false numerical solutions in ship motion simulation.

To help understand this, we consider a very mathematically simple case: a ship maneuvering in a single environment wave (𝑀=𝑁=1) (and is called Case 1 in the rest of the paper). The ship heading direction changes at a steady rate (with the rudder angle 𝛿=15^∘). The value of the ship speed is also fixed and slow, with the Froude number 𝐹_𝑟=0.266. We first model this simple system with the existing version of DiSSEL [1]. The numerical driving forces (moments) and the restoring forces (moments) for heave, pitch, and roll are shown in Figure 1. From the figure one can observe that, as soon as the ship direction starts to change, the balance between the forces (moments) is lost numerically.

(a)

(b)

(c)

Figure 1 

The numerical driving and restoring forces (moments) for heave (a), pitch (b), and roll (c) from the original DiSSEL. The ship hull is EHF, with a rudder angle 𝛿=15∘ and the Froude number 𝐹𝑟=0.266. The incident wave height is 𝑤ℎ=1.496 m, and its wavelength normalized by ship length is 𝜆=1.0 (Case 1). The superscripts “𝑝” and “𝑞” represent the driving forces (moments) and restore forces (moments) which associated with (4c), (4d), and (5c), and (5d).

What is the cause of this imbalance? One possibility is the time variation of the encounter frequency 𝜔^en. To verify this, we consider Case 2: an artificial system in which 𝜔^en is fixed, but the ship heading direction change remains the same as in the previous case. Numerical results for this case are shown in Figure 2. Clearly, the force balance is well maintained for this artificial system. Of course, this is not a mathematical proof, but it is sufficient to demonstrate the key problem for retaining numerical force balances in ship motion modeling.

(a)

(b)

(c)

Figure 2 

Same as Figure 1, but with the encounter frequency fixed to 𝜔en=𝜔sea (Case 2).

3.2. Adaptive Method for the Time-Varying Encounter Frequencies

The above numerical experiments can guide our effort in developing a new algorithm to maintain dynamic balance for time-varying encounter frequencies. In particular, we intend to exploit the fact that the force balances can be well achieved for a steady encounter frequency.

A simple and effective approach is to apply the same numerical algorithm in a “fictitious” reference frame in which the encounter frequency is time invariant. To explain this approach, we will start from a single-wave scenario.

In this case, we introduce a “fictitious” time 𝑡^∗ and its corresponding time step Δ𝑡^∗, such that ||||𝜔^en(𝑡)Δ𝑡𝐤||||=||||𝜔^seaΔ𝑡^∗𝐤||||.(24) From (23), Δ𝑡^∗=|||1−|||𝐮𝐜|||cos𝛼^sea|||Δ𝑡.(25) But (25) could lead to Δ𝑡^∗=0 when 1−|||𝐮𝐜|||cos𝛼^sea=0,(26) that is, |𝐜|=||||||𝐮𝐜|||cos𝛼^sea|||.(27)

To avoid this numerical singularity, we modify (25) by adding a truncation limit:Δ𝑡^∗=⎧⎪⎨⎪⎩Δ𝜏^∗≡𝜀,if|||1−|||𝐮𝐜|||cos𝛼^sea|||1−|||𝐮|||Δ𝑡≤𝜀𝐜|||cos𝛼^sea|||Δ𝑡,Δ𝜏_∗=0,otherwise.(28) Substituting (28) into (22) and setting 𝑀=𝑁=1 (single-model wave), we have𝜂(𝑥,𝑦,𝑡)=𝐴sin𝐤⋅𝐱−𝐱_𝑐+𝜇𝜔^sea𝑡^∗−𝜏^∗,(29a) where 𝜇=sgn1−|||𝐮𝐜|||cos𝛼^sea,(29b)

and 𝜏^∗=∑Δ𝜏^∗, 𝑡^∗=∑Δ𝑡^∗.

To test the applicability of this approach, we reexamine Case 1. The force balance is now well retained with the new algorithm as shown in Figure 3. By comparing this with Figure 1, we can observe clearly the improvement by the new algorithm.

(a)

(b)

(c)

Figure 3 

Same as Figure 1 (Case 1), but the forces are calculated by the new algorithm defined in (28) and (29a)-(29b).

We test further the new algorithm with more complicated cases that have experimental data. Due to space limitations, we show here only two of them: the EHF maneuvering in regular waves with the largest rudder angle 𝛿=30^∘ (Case 3) and with the fastest forward speed 𝐹_𝑟=0.368 (Case 4). The results for the two cases are shown in Figures 4 and 5. Because in the experiments, the ship motions, wave elevations, and trajectory motions were collected by separate systems, all the data sets could not be aligned in time accurately enough to know the initial conditions of the system. Therefore, the benchmarking comparison focuses only on the overall characteristics of the motion and not a point-by-point prediction comparison. As shown in these figures, our numerical results agree very well with the experimental data. The good agreement is because the dynamical balances are now correctly modeled in the simulation, as shown in Figures 4(c) and 5(c).

Figure 4 

(a) Roll and pitch motions of EHF with 𝛿=30∘, 𝐹𝑟=0.2146 in the regular incident wave 𝑤ℎ=5.984 m, 𝜆=2.0 (Case 3). The numerical results from the improved DiSSEL are on the top, and the experimental data are in the bottom. (b) Same as (a), but for the ship track in Earth coordinate, 𝑥 is the initial direction, and 𝑦 is perpendicular to 𝑥 direction. (c) Same as (a), but for the driving and restoring forces (moments).

Figure 5 

(a) Same as Figure 4(a), but for the case 𝐹𝑟=0.368, 𝛿=14.4∘, 𝑤ℎ=2.992 m, 𝜆=1.0 (Case 4). (b) Same as Figure 4(b), but for Case 4. (c) Same as Figure 4(c), but for Case 4.

This new algorithm can be easily extended to arbitrary two-dimensional irregular waves. The mathematical details are shown as follows: expression (29a) is now modified to irregular waves:𝜂(𝑥,𝑦,𝑡)=^𝑁_𝑀𝑖=1_𝑗=1=_𝑁𝐴_𝑖𝑗sin𝐤_𝑖𝑗⋅𝐱−𝐱_𝑐+𝜔^en_𝑖𝑗𝑡+𝜙_𝑖𝑗_𝑀𝑖=1_𝑗=1𝐴_𝑖𝑗sin𝐤_𝑖𝑗⋅𝐱−𝐱_𝑐+𝜙_𝑖𝑗+𝜇𝜔^sea_𝑖𝑗𝑡^∗𝑖𝑗−𝜏^∗𝑖𝑗,(30) where Δ𝑡^∗𝑖𝑗 and Δ𝜏^∗𝑖𝑗 are defined similarly as in (28) for each wave mode 𝐜_𝑖𝑗 and 𝛼^sea_𝑖𝑗.

To maintain computational efficiency, we introduce an effective incident wave frequency and wave amplitude:𝜔^sea≡∑^{𝑀,𝑁𝑖𝑗}𝐴_𝑖𝑗𝜔^sea_𝑖𝑗𝐴^Σ𝐴_Σ,≡^𝑀,𝑁_𝑖𝑗𝐴_𝑖𝑗.(31) With these definitions, we can then define a single effective “fictitious” time stepΔ𝑡^∗≡1𝜔^sea𝐴^Σ𝑀,𝑁_𝑖,𝑗𝐴_𝑖𝑗𝜔^sea_𝑖𝑗Δ𝑡^∗𝑖𝑗.(32) The notions of the effective incident wave frequency and the “fictitious” time are the same as those in the one-dimensional regular wave formulations. But the actual physical meanings are different. In particular, 𝜔^sea is given and invariant throughout the simulation. But Δ𝑡^∗changes in time because of time-varying Δ𝑡^∗𝑖𝑗. Therefore, (30) can be rewritten as𝜂(𝑥=,𝑦,𝑡)_𝑁_𝑀𝑖=1_𝑗=1𝐴_𝑖𝑗sin𝐤_𝑖𝑗⋅𝐱−𝐱_𝑐+𝜙_𝑖𝑗+𝜇𝜔^sea_𝑖𝑗𝑐^∗𝑖𝑗𝑡^∗−𝜏^∗𝑖𝑗,(33) where𝑐^∗𝑖𝑗≡𝑡^∗𝑖𝑗𝑡^∗.(34) Again, Δ𝑡 and Δ𝑡^∗𝑖𝑗 (and therefore Δ𝑡^∗ by definition) must satisfy the CFL condition in the simulation.

It should be pointed that this method is more efficient than iterative methods not only because it calculates the mean Δ𝑡^∗𝑖𝑗 once in each time step, but also because it does not impose any convergence requirement and thus ensures a stable solution.

Since there is no experimental data on ship maneuvering in irregular seas with which to compare, we will only examine the dynamical balances in the numerical solutions, with no attempt made for benchmarking purposes. For this, we reconsider Case 4, but replace the regular wave with the irregular wave specified by the spectral shape in Figure 6. The numerical results are summarized in Figure 7. From the figures we can observe that the dynamical balances are well established in the numerical solutions for roll, pitch, and heave motions, that is, the driving forces (moments) are comparable in magnitude and opposite in sign to the restoring forces (moments).

Figure 6 

Irregular incident wave spectrum with the peak frequency 𝜔𝑝=0.586 and the significant wave height, SWH=6 m.

(a)

(b)

(c)

Figure 7 

The numerical driving and restoring forces (moments) from the improved DiSSEL for Case 4, but with the irregular incident waves specified in Figure 6.

4. Conclusion

When a ship maneuvers in a seaway, arbitrary incident waves measured on the ship are time varying (i.e., time-varying encounter frequencies). This time variation depends on the ship motion state, which, in turn, changes in time under the influence of the incident waves. Therefore, it is critical to maintain the dynamical balances between the driving forces (moments) and the restoring forces (moments) in a numerical simulation for the accurate modeling of ship motion. Such dynamical balance depends on the appropriate evaluation of the incident waves and the corresponding forces on the ship hull and thus requires well-designed numerical approaches to model the nonlinear dependencies between the environment and the ship motion.

We have demonstrated numerically that direct time integration of the ship-wave interactions cannot retain properly the balance when the ship heading changes in time. This heading change results in time-varying encounter wave frequencies. We anticipated that this time-varying encounter frequency could lead to phase misfits between the driving forces and the restoring forces, thus destroying the dynamical balances in the simulation and therefore the numerical stability. This is supported by our numerical simulations of Cases 1 and 2 (shown in Figures 1 and 2).

To eliminate the time variation of the encounter frequency (mainly due to ship heading direction change), we developed a new adaptive algorithm with a mathematical (but “fictitious”) time 𝑡^∗ defined in (32). This is equivalent to rescaling the numerical time (and thus improving numerical temporal resolution) dynamically with the ship motion states. In this new system, the encounter frequencies are steady as measured in 𝑡^∗. We have tested this new algorithm for many cases. All numerical results show that the dynamical balances between the driving forces (moments) and the restoring forces (moments) are well established.

This new algorithm has been implemented in the DiSSEL nonlinear ship motion model. With this improved model, we are able to reproduce satisfactorily the results from all known laboratory experiments. Results for two of such experiments are shown in Figure 4 (with the largest rudder angle) and Figure 5 (with the largest forward speed) in regular waves. We also tested the improved model for ship maneuvering in irregular waves specified in Figure 6, and the dynamical balances are retained as well, as shown in Figure 7.

Acknowledgments

The authors thank Terry Applebee, Y. T. Shen, and J. Gorski for their comments and suggestions to this study. They also thank the Seakeeping Division at NSWCCD for providing the experimental data to benchmark and validate their model. This work is supported by Naval Surface Warfare Center Independent Laboratory In-House Research (ILIR) Program.