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

Authors: Ray-Qing Lin; Weijia Kuang
Journal: Modelling and Simulation in Engineering
 (2011)

---

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


---

## Body

## 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 Section2 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 (V=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:(1)∇2φ≡∇h2φ+∂2φ∂2z=0,(-H≤z≤η),∂φ∂t+[12∇φ-(us+vh)]⋅∇φ+pρw+gz-∂(us+vh)∂t⋅x-ν∇h2φ=0,(-H≤z≤η),∂η∂t+(∇hη)⋅[∇hφ-(us+vh)]=∂φ∂z,(z=η),
where p is the pressure, ρw is the water density, ν is the turbulent viscosity (arising from the wave-breaking processes), x is the position vector, H is the water depth, φ is the velocity potential, η is the free surface elevation, ∇(≡x̂∂/∂x+ŷ∂/∂y+ẑ∂/∂z) is the gradient, ∇h is the horizontal component part of ∇, and us is the ship maneuvering velocity (in SRF, it is the given forward speed). It should be noted that ∂us/∂t includes both magnitude and direction changes of the maneuvering velocity. In addition, vh is the horizontal component of translational velocity vc for surge and sway (absent in the fixed heading scenario).The impenetrable boundary conditions are(2)∂φ∂z=0
at the flat bottom z=-H, and (3)n̂Σ⋅∇φ=n̂Σ⋅[us+vc+Ω×(xΣ-xc)]
on the ship surface xΣ. In the above equation, n̂Σ is the normal unit vector of the ship surface xΣ, and xc 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 centerxc and the rotation about xc. The translational motion, vc, is described by the following equation in MRF: (4a)msdvcdt+Dνvc=F,
where ms is the total ship mass and F is the total force on the ship:    (4b)F=Fp+Fg+Fr,
where (4c)Fp≡-∬Σwetn̂ΣpdΣ
is the (driving) pressure force on the ship hull,(4d)Fg=(ρwVwet-ms)g,is the buoyancy force (restoring force). This force acts only in the vertical direction. But the rudder forceFr in (4b) is the new addition. We will discuss this in more detail in the next section. Vwet is under water volume.The solid body rotation about the ship mass center is governed by the Liouville equation defined in SRF:(5a)I⋅dΩdt+Ω×(I⋅Ω)+DΩ(I⋅Ω)=Γ,
where I is the moment of inertia tensor about the ship mass center xc (and its components are denoted by I(i,j)), Dν and DΩ are integrated dissipative effects (e.g., drag) of the fluid on the ship hull, and Γ is the total moment (torque) on the ship:(5b)Γ=Γp+Γg+Γr.
It includes the (driving) pressure torque:(5c)Γp=-∬ΣwetdΣ(xΣ-xc)×n̂Σp,
(where Σwet is total wetted ship surface), the restoring torque: (5d)Γg=(xwet-xc)×[(ms-ρwVwet)g],(wherexwet is the position vector in under water volum), and the torque Γr 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 forceFr in (4b) and the moment Γr 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)Fr=∬An̂rprds=Fxrx̂s+Fyrŷs+Fzrẑs,(7)Γr=∬AR×n̂rprds+Rt×FΩ=Γxrx̂s+Γyrŷs+Γzrẑs+(x-xC)×FΩ,
where pr is the potential flow pressure; but A=AΣ-Ablock (AΣ is the wetted rudder surface, and Ablock is the blocking effected area) [16], n̂r is the unit normal vector of A, xC  is the position vector of the center of the circle tangent to the point x on the ship trajectory, Rt is the radius vector of ship trajectory to the point x, and(8)R≡Rxx̂s+Ryŷs+Rzẑs
is the rotating arm vector of the rudder. The force FΩ 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):(9)FΩ=∬Aρ[2Ω*×vc+Ω*×(Ω*×x)+dΩ*dt×x]ds,Ω*=dβdt+Ωzẑs,
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 Fr can be substantially simplified:(10)Fxr=-Cl12ρur2Asin2δ+∬An̂xrpwrds,Fyr=Cl14ρur2Asin⁡2δ+∬An̂yrpwrds,
where ur is the effective inflow speed; n̂xr and n̂yr are the x and y components of the normal vector of the rudder surface, respectively; pwr is pressure on rudder surface due to the incident waves; Cl is the lift curve slope:(11)Cl=2π1+(2/(L/d)),
where L and d 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 A (not AΣ) in (6)–(10).By (10), the vertical moment is   (12)Γzr=RxFyr-RyFxr≈RxFyr=(Cl14ρur2Asin⁡2δ+∬An̂yrpwrds)(xr-xc),
since Ry is in general much smaller than Rx. In (12), (xr-xc) is the longitudinal distance between the ship mass center xc and the rudder mass center xr.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(13)ur=1.6us,
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) (x̂e,ŷe,ẑe).  The ship moving direction angle, θ*, is the sum of the yaw angle θz (the rotational angle in z direction, and dθz/dt=Ωz obtained from (5a)–(5d)) and the drift angle β:(14)θ*=θz+β,β≡tan-1(νy|us|+νx),
where νx and νy are ship translation motion in the x and y directions in ship coordinate, respectively. The sketch map can be find in Figure 5 page 198 [23, Figure 5, page 198]. The track xtrack of the ship in ERF is defined as (15)xtrack=xtrackx̂e+ytrackŷe.
In particular,(16)dxtrackdt=Am⋅(us+vh)≡u,Am=[cos⁡θz-sin⁡θzsin⁡θzcos⁡θz].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 xtrack(t), the total ship horizontal velocity, u, and horizontal acceleration can be calculated by simple time differentiation, and the direction angle θ* can be also determined via(17)θ*=tan-1(u⋅ŷeu⋅x̂e).
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 centerC of the track curve:
(18)us=ωẑe×(xtrack-xC),
where ω is the rotation speed about C and xC is the position vector of the point C in ERF. Therefore,(19)∂us∂t=ω2ẑe×[ẑe×(xtrack-xC)].
Thus,(20)pρ+∂us∂t⋅xtrack=1ρ[p+12ρω2|ẑe×(xtrack-xC)|2]≡p̃ρ.
In other words, the changes in the direction of the ship motion can be described by the effective pressure p̃ in (20). Given the ship track xtrack(t), the radius |xtrack-xC| can then be calculated with the standard curvature formulation: (21)|xtrack-xC|=1κ(t),κ(t)=|(dxtrack/dt)(d2ytrack/dt2)-(dytrack/dt)(d2xtrack/dt2)|[(dxtrack/dt)2+(dytrack/dt)2]3/2. Its influence on the solid body rotation of the ship can be computed via integration of the effective pressure over the ship surface.


## 2.1. Ship-Wave Interaction Model

In DiSSEL, the surface waves are described by the following set of partial differential equations:(1)∇2φ≡∇h2φ+∂2φ∂2z=0,(-H≤z≤η),∂φ∂t+[12∇φ-(us+vh)]⋅∇φ+pρw+gz-∂(us+vh)∂t⋅x-ν∇h2φ=0,(-H≤z≤η),∂η∂t+(∇hη)⋅[∇hφ-(us+vh)]=∂φ∂z,(z=η),
where p is the pressure, ρw is the water density, ν is the turbulent viscosity (arising from the wave-breaking processes), x is the position vector, H is the water depth, φ is the velocity potential, η is the free surface elevation, ∇(≡x̂∂/∂x+ŷ∂/∂y+ẑ∂/∂z) is the gradient, ∇h is the horizontal component part of ∇, and us is the ship maneuvering velocity (in SRF, it is the given forward speed). It should be noted that ∂us/∂t includes both magnitude and direction changes of the maneuvering velocity. In addition, vh is the horizontal component of translational velocity vc for surge and sway (absent in the fixed heading scenario).The impenetrable boundary conditions are(2)∂φ∂z=0
at the flat bottom z=-H, and (3)n̂Σ⋅∇φ=n̂Σ⋅[us+vc+Ω×(xΣ-xc)]
on the ship surface xΣ. In the above equation, n̂Σ is the normal unit vector of the ship surface xΣ, and xc 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 centerxc and the rotation about xc. The translational motion, vc, is described by the following equation in MRF: (4a)msdvcdt+Dνvc=F,
where ms is the total ship mass and F is the total force on the ship:    (4b)F=Fp+Fg+Fr,
where (4c)Fp≡-∬Σwetn̂ΣpdΣ
is the (driving) pressure force on the ship hull,(4d)Fg=(ρwVwet-ms)g,is the buoyancy force (restoring force). This force acts only in the vertical direction. But the rudder forceFr in (4b) is the new addition. We will discuss this in more detail in the next section. Vwet is under water volume.The solid body rotation about the ship mass center is governed by the Liouville equation defined in SRF:(5a)I⋅dΩdt+Ω×(I⋅Ω)+DΩ(I⋅Ω)=Γ,
where I is the moment of inertia tensor about the ship mass center xc (and its components are denoted by I(i,j)), Dν and DΩ are integrated dissipative effects (e.g., drag) of the fluid on the ship hull, and Γ is the total moment (torque) on the ship:(5b)Γ=Γp+Γg+Γr.
It includes the (driving) pressure torque:(5c)Γp=-∬ΣwetdΣ(xΣ-xc)×n̂Σp,
(where Σwet is total wetted ship surface), the restoring torque: (5d)Γg=(xwet-xc)×[(ms-ρwVwet)g],(wherexwet is the position vector in under water volum), and the torque Γr 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 forceFr in (4b) and the moment Γr 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)Fr=∬An̂rprds=Fxrx̂s+Fyrŷs+Fzrẑs,(7)Γr=∬AR×n̂rprds+Rt×FΩ=Γxrx̂s+Γyrŷs+Γzrẑs+(x-xC)×FΩ,
where pr is the potential flow pressure; but A=AΣ-Ablock (AΣ is the wetted rudder surface, and Ablock is the blocking effected area) [16], n̂r is the unit normal vector of A, xC  is the position vector of the center of the circle tangent to the point x on the ship trajectory, Rt is the radius vector of ship trajectory to the point x, and(8)R≡Rxx̂s+Ryŷs+Rzẑs
is the rotating arm vector of the rudder. The force FΩ 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):(9)FΩ=∬Aρ[2Ω*×vc+Ω*×(Ω*×x)+dΩ*dt×x]ds,Ω*=dβdt+Ωzẑs,
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 Fr can be substantially simplified:(10)Fxr=-Cl12ρur2Asin2δ+∬An̂xrpwrds,Fyr=Cl14ρur2Asin⁡2δ+∬An̂yrpwrds,
where ur is the effective inflow speed; n̂xr and n̂yr are the x and y components of the normal vector of the rudder surface, respectively; pwr is pressure on rudder surface due to the incident waves; Cl is the lift curve slope:(11)Cl=2π1+(2/(L/d)),
where L and d 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 A (not AΣ) in (6)–(10).By (10), the vertical moment is   (12)Γzr=RxFyr-RyFxr≈RxFyr=(Cl14ρur2Asin⁡2δ+∬An̂yrpwrds)(xr-xc),
since Ry is in general much smaller than Rx. In (12), (xr-xc) is the longitudinal distance between the ship mass center xc and the rudder mass center xr.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(13)ur=1.6us,
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) (x̂e,ŷe,ẑe).  The ship moving direction angle, θ*, is the sum of the yaw angle θz (the rotational angle in z direction, and dθz/dt=Ωz obtained from (5a)–(5d)) and the drift angle β:(14)θ*=θz+β,β≡tan-1(νy|us|+νx),
where νx and νy are ship translation motion in the x and y directions in ship coordinate, respectively. The sketch map can be find in Figure 5 page 198 [23, Figure 5, page 198]. The track xtrack of the ship in ERF is defined as (15)xtrack=xtrackx̂e+ytrackŷe.
In particular,(16)dxtrackdt=Am⋅(us+vh)≡u,Am=[cos⁡θz-sin⁡θzsin⁡θzcos⁡θz].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 xtrack(t), the total ship horizontal velocity, u, and horizontal acceleration can be calculated by simple time differentiation, and the direction angle θ* can be also determined via(17)θ*=tan-1(u⋅ŷeu⋅x̂e).
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 centerC of the track curve:
(18)us=ωẑe×(xtrack-xC),
where ω is the rotation speed about C and xC is the position vector of the point C in ERF. Therefore,(19)∂us∂t=ω2ẑe×[ẑe×(xtrack-xC)].
Thus,(20)pρ+∂us∂t⋅xtrack=1ρ[p+12ρω2|ẑe×(xtrack-xC)|2]≡p̃ρ.
In other words, the changes in the direction of the ship motion can be described by the effective pressure p̃ in (20). Given the ship track xtrack(t), the radius |xtrack-xC| can then be calculated with the standard curvature formulation: (21)|xtrack-xC|=1κ(t),κ(t)=|(dxtrack/dt)(d2ytrack/dt2)-(dytrack/dt)(d2xtrack/dt2)|[(dxtrack/dt)2+(dytrack/dt)2]3/2. 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(22)η(x,y,t)=∑i=1N∑j=1MAijsin⁡{kij⋅(x-xc)+ωijent+ϕij},
where Aij is the amplitude of the mode (i,j), kij is the wave number vector, ωijen is the encounter frequency, and ϕij is the phase. The integer pair (M,N) defines the total number of the wave modes used in numerical simulation. The encounter frequency ωijen is defined as (23)ωijen=ωijsea[1-|ucij|cos⁡αijsea],
where cij is the phase velocity of the wave mode, αijsea is the angle between cij and u, and ωijsea is the incident wave frequency in ERF. For any given incident wave, αijsea changes with the ship heading direction.One could clearly observe the complications and its numerical difficulties from (23): changes in vh and θz (and therefore αijsea) are dynamically determined via the driving forces and the restoring forces on the ship hull which depend on the encounter frequency ωijen. But changes in ωijen depend also on vh and θz. 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 (M=N=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 Fr=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.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 Fr=0.266. The incident wave height is wh=1.496 m, and its wavelength normalized by ship length is λ=1.0 (Case 1). The superscripts “p” and “q” represent the driving forces (moments) and restore forces (moments) which associated with  (4c), (4d), and (5c), and (5d).
(a)(b)(c)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.Same as Figure1, but with the encounter frequency fixed to ωen=ωsea (Case 2).
(a)(b)(c)


### 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” timet* and its corresponding time step Δt*, such that (24)|ωen(t)Δtk|=|ωseaΔt*k|.
From (23), (25)Δt*=|1-|uc|cos⁡αsea|Δt.
But (25) could lead to Δt*=0 when (26)[1-|uc|cos⁡αsea]=0,
that is,
(27)|c|=||uc|cos⁡αsea|.To avoid this numerical singularity, we modify (25) by adding a truncation limit:(28)Δt*={Δτ*≡ε,if|1-|uc|cos⁡αsea|Δt≤ε|1-|uc|cos⁡αsea|Δt,Δτ*=0,otherwise.
Substituting (28) into (22) and setting M=N=1 (single-model wave), we have(29a)η(x,y,t)=Asin⁡{k⋅(x-xc)+μωsea(t*-τ*)},
where
(29b)μ=sgn⁡{1-|uc|cos⁡αsea},andτ*=∑Δτ*, t*=∑Δt*.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 Figure3. By comparing this with Figure 1, we can observe clearly the improvement by the new algorithm.Same as Figure1 (Case 1), but the forces are calculated by the new algorithm defined in (28) and (29a)-(29b).
(a)(b)(c)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 Fr=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°, Fr=0.2146 in the regular incident wave wh=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, x is the initial direction, and y is perpendicular to x direction. (c) Same as (a), but for the driving and restoring forces (moments).Figure 5
(a) Same as Figure4(a), but for the case Fr=0.368, δ=14.4°, wh=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:(30)η(x,y,t)=∑i=1N∑j=1MAijsin⁡{kij⋅(x-xc)+ωijent+ϕij}=∑i=1N∑j=1MAijsin⁡{kij⋅(x-xc)+ϕij+μωijsea(tij*-τij*)},
where Δtij* and Δτij* are defined similarly as in (28) for each wave mode cij and αijsea.To maintain computational efficiency, we introduce an effective incident wave frequency and wave amplitude:(31)ωsea≡∑ijM,NAijωijseaAΣ,AΣ≡∑ijM,NAij.
With these definitions, we can then define a single effective “fictitious” time step(32)Δt*≡1ωseaAΣ∑i,jM,NAijωijseaΔtij*.
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 Δt*changes in time because of time-varying Δtij*. Therefore, (30) can be rewritten as(33)η(x,y,t)=∑i=1N∑j=1MAijsin⁡{kij⋅(x-xc)+ϕij+μωijsea(cij*t*-τij*)},
where(34)cij*≡tij*t*.
Again, Δt and Δtij* (and therefore Δt* 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Δtij* 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 Figure6. 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ωp=0.586 and the significant wave height, SWH=6 m.The numerical driving and restoring forces (moments) from the improved DiSSEL for Case 4, but with the irregular incident waves specified in Figure6.
(a)(b)(c)


## 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(22)η(x,y,t)=∑i=1N∑j=1MAijsin⁡{kij⋅(x-xc)+ωijent+ϕij},
where Aij is the amplitude of the mode (i,j), kij is the wave number vector, ωijen is the encounter frequency, and ϕij is the phase. The integer pair (M,N) defines the total number of the wave modes used in numerical simulation. The encounter frequency ωijen is defined as (23)ωijen=ωijsea[1-|ucij|cos⁡αijsea],
where cij is the phase velocity of the wave mode, αijsea is the angle between cij and u, and ωijsea is the incident wave frequency in ERF. For any given incident wave, αijsea changes with the ship heading direction.One could clearly observe the complications and its numerical difficulties from (23): changes in vh and θz (and therefore αijsea) are dynamically determined via the driving forces and the restoring forces on the ship hull which depend on the encounter frequency ωijen. But changes in ωijen depend also on vh and θz. 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 (M=N=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 Fr=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.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 Fr=0.266. The incident wave height is wh=1.496 m, and its wavelength normalized by ship length is λ=1.0 (Case 1). The superscripts “p” and “q” represent the driving forces (moments) and restore forces (moments) which associated with  (4c), (4d), and (5c), and (5d).
(a)(b)(c)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.Same as Figure1, but with the encounter frequency fixed to ωen=ωsea (Case 2).
(a)(b)(c)


## 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” timet* and its corresponding time step Δt*, such that (24)|ωen(t)Δtk|=|ωseaΔt*k|.
From (23), (25)Δt*=|1-|uc|cos⁡αsea|Δt.
But (25) could lead to Δt*=0 when (26)[1-|uc|cos⁡αsea]=0,
that is,
(27)|c|=||uc|cos⁡αsea|.To avoid this numerical singularity, we modify (25) by adding a truncation limit:(28)Δt*={Δτ*≡ε,if|1-|uc|cos⁡αsea|Δt≤ε|1-|uc|cos⁡αsea|Δt,Δτ*=0,otherwise.
Substituting (28) into (22) and setting M=N=1 (single-model wave), we have(29a)η(x,y,t)=Asin⁡{k⋅(x-xc)+μωsea(t*-τ*)},
where
(29b)μ=sgn⁡{1-|uc|cos⁡αsea},andτ*=∑Δτ*, t*=∑Δt*.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 Figure3. By comparing this with Figure 1, we can observe clearly the improvement by the new algorithm.Same as Figure1 (Case 1), but the forces are calculated by the new algorithm defined in (28) and (29a)-(29b).
(a)(b)(c)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 Fr=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°, Fr=0.2146 in the regular incident wave wh=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, x is the initial direction, and y is perpendicular to x direction. (c) Same as (a), but for the driving and restoring forces (moments).Figure 5
(a) Same as Figure4(a), but for the case Fr=0.368, δ=14.4°, wh=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:(30)η(x,y,t)=∑i=1N∑j=1MAijsin⁡{kij⋅(x-xc)+ωijent+ϕij}=∑i=1N∑j=1MAijsin⁡{kij⋅(x-xc)+ϕij+μωijsea(tij*-τij*)},
where Δtij* and Δτij* are defined similarly as in (28) for each wave mode cij and αijsea.To maintain computational efficiency, we introduce an effective incident wave frequency and wave amplitude:(31)ωsea≡∑ijM,NAijωijseaAΣ,AΣ≡∑ijM,NAij.
With these definitions, we can then define a single effective “fictitious” time step(32)Δt*≡1ωseaAΣ∑i,jM,NAijωijseaΔtij*.
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 Δt*changes in time because of time-varying Δtij*. Therefore, (30) can be rewritten as(33)η(x,y,t)=∑i=1N∑j=1MAijsin⁡{kij⋅(x-xc)+ϕij+μωijsea(cij*t*-τij*)},
where(34)cij*≡tij*t*.
Again, Δt and Δtij* (and therefore Δt* 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Δtij* 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 Figure6. 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ωp=0.586 and the significant wave height, SWH=6 m.The numerical driving and restoring forces (moments) from the improved DiSSEL for Case 4, but with the irregular incident waves specified in Figure6.
(a)(b)(c)


## 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 Figures1 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”) timet* 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 t*. 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 Figure4 (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.



---

*Source: 356741-2011-10-09.xml*