Nonlinear Time-Domain Strip Theory Formulation for Low-Speed Manoeuvring and Station-Keeping

This paper presents a computer effective nonlinear time-domain strip theory formulation for dynamic positioning (DP) and low-speed manoeuvring. Strip theory or 2D potential theory, where the ship is divided in 20 to 30 cross sections, can be used to compute the potential coefficients (added mass and potential damping) and the exciting wave loads (Froude-Krylov and diffraction forces). Commercially available programs are ShipX (VERES) by Marintek (Fathi, 2004) and SEAWAY byAmarcon (Journée &Adegeest, 2003), for instance. The proposed method can easily be extended to utilize other strip theory formulations or 3-D potential programs like WAMIT (2004). The frequency dependent potential damping, which in classic theory results in a convolution integral not suited for real-time simulation, is compactly represented by using the state-space formulation of Kristiansen & Egeland (2003). The separation of the vessel model into a low-frequency model (represented by zerofrequency added mass and damping) and a wave-frequency model (represented by motion transfer functions or RAOs), which is commonly used for simulation, is hence made superfluous. Transformations of motions and coefficients between different coordinate systems and origins, i.e. data frame, hydrodynamic frame, body frame, inertial frame etc., are put into the rigid framework of Fossen (1994, 2002). The kinematic equations of motion are formulated in a compact nonlinear vector representation and the classical kinematic assumption that the Euler angles are small is removed. This is important for computation of accurate control forces at higher roll and pitch angles. The hydrodynamic forces in the steadily translating hydrodynamic reference frame (equilibrium axes) are, however, assumed to be linear. Recipes for computation of retardation functions are presented and frequency dependent viscous damping is included. Emphasis is placed on numerical computations and representation of the data from VERES and SEAWAY in Matlab/Simulink. For this purpose a Simulink add-in to the Marine Systems Simulator (MSS) at the Norwegian University of Science and Technology has been developed (Fossen et al., 2004).


Introduction
In 1949 Ursell published his famous paper on potential theory for determining the hydrodynamic coefficients of semicircular cross sections, oscillating in deep water in the frequency domain (Ursell, 1949).This was used as a rough estimation for zero speed ship applications.Motivated by this Grim (1953), Tsai (1959, 1960, 1961), Gerritsma (1960) and others used conformal mapping techniques like the Lewis conformal mapping to transform ship-like cross sections to semicircles such that more realistic hull forms could be calculated.Exciting wave loads were computed using undisturbed regular waves.Denis & Pierson (1953) published a superposition method to describe the irregular waves assuming that the sea could be described as a sum of many simple harmonic waves; each wave with its own frequency, amplitude, direction and random phase lag.The responses of the ship at zero speed were calculated for each of these individual harmonic waves and superposed.
The extension to forward speed was made available by Korvin-Kroukovsky & Jacobs (1957), and which was further improved in the sixties.Later Frank (1967) published a pulsating source theory to calculate the hydrodynamic coefficients of a cross section of a ship in deep water directly, without using conformal mapping.Keil (1974) published a theory for obtaining the potential coefficients in very shallow water using Lewis conformal mappings.The S-T-F strip theory (Salvesen et al., 1970), which accounts for forward speed as well as transom stern effects, is made available through the program ShipX (VERES) (Fahti, 2004).
Since strip theory determines the hydrodynamic coefficients from potential theory it is common to calculate the added resistance of a ship due to waves e.g. by using the integrated pressure method by Boese (1970) or the radiated energy method (Gerritsma & Beukelman, 1972).In roll it is common to use the viscous correction by Ikeda et al. (1978) based on semi-empirical methods.
For zero speed 3-D potential theory can be used to compute the hydrodynamic coefficients (WAMIT, 2004), while forward speed effects still are considered to be difficult to solve.Consequently, the 2-D approach (strip theory) is still very favorable for calculating the behavior of a ship at forward speed.For a more detailed discussion on advantages and disadvantages when comparing 2-D with 3-D theories; see Faltinsen & Svendsen (1990).

Kinematics
In this section the 6 DOF equations of motion are derived using vectorial mechanics with emphasis placed on keeping the kinematics nonlinear while linear theory is assumed for hydrodynamic forces and moments.In this context it is convenient to define the vectors without reference to a coordinate frame (coordinate free vector).A vector v / is defined by its magnitude and direction.The vector v / 0 in the point O decomposed in reference frame n is denoted as vn 0 , which is also referred to as a coordinate vector.
2.1.1.Coordinate Systems Three orthogonal coordinate systems are needed to described the motions in 6DOF; see Figures 1 and 2: is assumed fixed on the Earth surface with X n -axis pointing North, the Y n -axis pointing East, and the Z n -axis down of the Earth tangent plane.The n-frame position pnó[n, e, d] T and Euler angles ó[{, h, t] T are defined in terms of the vector (Fossen, 2002):  Ω Hydrodynamic frame (h-frame): The hydrodynamic forces and moments are defined in a steadily translating hydrodynamic coordinate system X h Y h Z h moving along the path of the ship with the constant speed U with respect to the nframe.For DP Uó0.The X h Y h -plane is parallel to the still water surface, and the ship carries out oscillations around the moving frame -axis is positive towards starboard, and the X h axis is positive forwards.This is also referred to as the equilibrium axis system (Bailey et al., 1998).The coordinate origin is denoted W while the h-frame generalized position vector is: ó[p, q, r] T with respect to the n-frame are denoted as: Hydrodynamic programs for strip theory define a local reference frame with axes different from the h-frame.Consequently, the data sets from these programs must be transformed to the h-frame by using rotation matrices.

Rotation Matrices
The notation Rb a é SO(3) implies that the rotation matrix Rb a between two frames a and b (from a to b) is an element in SO(3), that is the special orthogonal group of order 3: The group SO( 3) is a subset of all orthogonal matrices of order 3, i.e.SO(3):O(3) where O( 3) is defined as: Hence it follows that A principal rotation a about the i-axis is denoted as R i?

Kinematic Equations of Motion
The linear velocity v / w of W and the angular velocity u / hn of the h-frame with respect to the n-frame are given by: where u / bn is the angular velocity of the h-frame with respect to the n-frame, and r / w is the vector from O to W. Decomposing these vectors into the b-frame gives: Introducing the screw transformation: we see that: Transformation from the b-to the h-frame: The transformation from the b-frame to the h-frame is done in terms of the rotation matrix: where dt, dh, and d{ are oscillations of the b-frame with respect to the h-frame.This is based on the assumption that the yaw angle oscillations dt with respect to the mean heading angle g* 6 are small such that the rotation matrix in yaw is close to the identity matrix.For many applications, the roll and pitch oscillations will be small, such that: From equation ( 14) it follows that: . The velocity transformation between the h and b frames then becomes: where J*(d ) é R6"6 is a generalized velocity transformation matrix: The force b in the b-frame can be transformed to the force h in the h-frame in a similar manner by: Transformation from the b-to the n-frame: The velocity transformation between the b and n frames is: where the Euler angle rotation matrix (zyx-convention) between the n and b frame is defined as the product of the three principal rotations: The Euler rates satisfies: ˙óT ( ) b bn (24 where T ( ) é R3"3 is the Euler angle attitude transformation matrix: Consequently: where J( ) é R6"6 is the velocity transformation matrix: 2.1.4.Transformation of Hydrodynamic Data The hydrodynamic data is defined in the program-specific data-frame with origin W (common to the h-frame), and must be transformed to the h-frame axes by the following transformation matrix: where Rh data é SO(3) is the rotation matrix between the data-frame and the h-frame found by combining the principal rotations equation ( 6).Then, it follows that the generalized rigid body inertia matrix Mdata RB , added inertia matrix Adata(u), damping matrix Bdata(u), spring stiffness matrix Gdata, and forces data can be rotated to the h-frame axes (denoted by *) by: since TT T óI 6"6 .This implies that the signs of the data-components have to be corrected while all model matrices are unchanged.For n/2 rotations this will not be the case.
The rotation matrix Rh data between the data-frame and the h-frame coordinate systems are defined below for two of the strip theory programs in commercial use: ShipX (VERES) by Marintek: For VERES (Fathi, 2004) the Z-axis is positive upwards and the X-axis is positive towards the stern.The program is based on the S-T-F strip theory formulation Salvesen et al. (1970).Hence, the h-frame rotation matrix becomes: In VERES the longitudinal center of gravity LCG* is given relative AP (positive forwards) when inputted to the program while the outputted LCG on the data files is related to L pp /2 (positive backwards).The vertical center of gravity VCG is given relative the baseline (positive upwards).Let T denote the water line for the actual load condition relative to the baseline (positive upwards), see Figure 1.Hence: (1) If the vessel motion point W in VERES is chosen as G it follows that: (2) Alternatively, the vessel motion point W can be chosen at T implying that: The corresponding quantities Adata, Bdata, Gdata and force amplitudes/angles are also computed with respect to W. SEAWAY by Amarcon: For SEAWAY (Journe ´e & Adegeest, 2003) the Z-axis is positive upwards and the X-axis is positive towards the bow.Hence, the h-frame rotation matrix becomes: Similarly, as for VERES the vectors rb w and rb g must be computed from the SEAWAY data sets, see Journe ´e & Adegeest (2003) for details.

Kinetics
The generalized forces acting on the vessel are found by formulating Newton's 2nd law in b-frame coordinates.Since the hydrodynamic forces and moments are computed in h-frame coordinates these will be transformed to the b-frame.

b-Frame Representation
Consider the 6DOF rigid-body equations of motion (Fossen, 2002): where H é R6 is a vector of generalized hydrodynamic forces, é R6 is the generalized propulsion and control forces and M RB é R6"6 is the rigid-body system inertia matrix.The generalized forces and inertia matrix are computed with respect to the coordinate origin O: where m is the mass, I 3"3 is the identity matrix, S(rb g ) é SS(3) defined as: and The inertia matrix I c is computed with respect to the center of gravity G and it can easily be transformed to the origin O by using the parallel axes theorem:

h-Frame Representation
The b-frame equations of motion can be transformed to the h-frame by using equation ( 18).Consequently: where we have assumed that d is small.Hence equation ( 41) can be written: The hydrodynamic forces in the h-frame is denoted as * H óJ*(d )\ T H and *óJ*(d )\ T such that: where * denotes vectors and matrices in the h-frame and:

Hydrodynamic Forces
The generalized hydrodynamic forces are (Faltinsen, 1990): The radiation-induced forces and moments * R are functions of frequency and time.In the next section a transformation will be applied to remove the frequency dependent quantities.

Time-Domain Representation
In this section a state-space model for effective time-domain simulation of the 6 DOF equations of motion is derived.The time-domain representation is based on the assumption that the oscillations d of the b-frame with respect to the h-frame are small such that: Notice that this assumption is only made for the transformation of the linear generalized hydrodynamic forces from the h-frame to the b-frame, while the nonlinear kinematics between the n-frame and the b-frame equation ( 26) are preserved and hence valid for large Euler angles .This is important from a feedback control point of view, where the nonlinear kinematics is exploited in the design.

Cummins Equation
From equations ( 53)-( 56) and ( 57) it follows that the b-frame representation is: and Notice that we have assumed that g( )BG where G is the generalized stiffness matrix.Cummins (1962) has shown that the frequency dependent terms A(u) and B(u) can be removed by writing the equations of motion in the following form: From Ogilvie (1964) it follows that: K(q)ó 2 n 0 (B(u)ñB(ê)) cos(uq)du ( 64) where M A é R6"6 is a constant (frequency independent) generalized added inertia matrix evaluated at the infinity frequency, K(q) é R6"6 is a time-varying matrix of retardation functions which can be computed off-line using the A(u) or B(u) data sets and equations ( 64)-( 65), since K(q) is causal.Kristiansen & Egeland (2003) have developed a state-space formulation for the potential damping term in equation ( 62).Consider:

Linear Zero Speed State-Space Representation of Cummins Equation
where K(tñq) is the retardation function.For causal systems: If (q) is a unit impulse, then (t) given by equation ( 66) will be an impulse response function.Consequently, (t) can be represented by a linear state-space model: ˙óA r òB r , (0)ó0 ( 68) where (A r , B r , C r , D r ) are constant matrices of appropriate dimensions.Applying the Laplace tranformation to equations( 68)-( 69), the potential damping term can be written as: where D p (s) é R6"6 is a transfer function matrix.Notice that the filter: now contains the memory effect of the fluid.Cummins equation ( 62) is now written as a Linear Time-Invariant (LTI) model: M where is computed using the state-space model equations ( 68)-( 69).This gives the time-domain model: óC r òD r (76) Notice that the property: also holds for this model since the generalized added mass matrix M A is frequency independent and symmetric.

Marine Systems Simulator (Matlab/Simulink)
In the Marine Systems Simulator two Matlab m-files for postprocessing of the VERES data are provided (Fossen et al., 2004).These are: Veres2ABC.mComputes the model matrices, retardation function state-space models etc.Output file: ABC.mat Veres2force.mCreates a table of generalized diffraction/Froude-Krylov force coefficients.Output file: Forces_TF.dat The data flow is shown in Figure 3 where the S-175 container ship is used as case study.The data file Forces_TF.dat is used as input for the Simulink program shown in Figure 4 while the file ABC.mat must be manually loaded into work space.The numerical recipes used in the postprocessing of the data are described in Section 5.

Numerical Recipes
The numerical recipes used in Veres2ABC.mare described below.

Low and High Frequency Approximations
For zero speed Uó0 it follows that the frequency dependent potential damping satisfies: such that: At high frequencies the tail of the damping curve is given by (Journe ´e, 1993): where b ii is a constant and ) h typically is less than 5 rad/s for merchant ships, see Figure 5.
If strip theory data is used in the interval 0\u\ ) h , numerical computation of K ii (q) is significantly improved by using the high frequency approximation B ii (u) for u[) h .The retardation functions equation ( 79) and the corresponding 5th-order state-space model approximations equations ( 68)-( 69) are shown in Figure 6 for the S-175 container ship.

Viscous Effects
In general it is difficult to compute accurate estimates of the viscous damping unless CFD is used.A frequently used approximation is to assume a viscous ramp for the diagonal terms (Bailey et al., 1998): where b i [0 can be chosen for instance as 5-20% of the maximum value max S (B ii (u)).Hence: Alternatively, the viscous damper could be chosen as an exponential function:

Numerical Results for the S-175 Container Ship
This section presents frequency dependent added mass and potential damping including viscous effects (Figures 8-9), time series of hydrodynamic excitation forces (Figure 10) and vessel motions in heave, roll and pitch (Figure 11).The retardation functions are given in Figure 6.The results are based on data from ShipX (VERES).The main particulars of the S-175 container ship are given in Table 1, while the body plan is shown in Figure 7.The numerical results are computed for beam seas with the JONSWAP wave spectrum using significant wave height H 13 ó5 m and peak frequency u p ó0.56 rad/s.The wave spreading factor was set to 4.

Conclusion
In this paper a computer effective nonlinear time-domain strip theory formulation for low-speed applications has been presented.The proposed model can be used to simulate dynamic positioning (DP) and low-speed manoeuvres for ships and floating rigs in real time.The model is also well suited for feedback control design since it incorporates the effect of varying sea states.
The nonlinear kinematic equations of motion are represented in terms of Euler angles while the kinetics (rigid-body dynamics and hydrodynamic generalized forces) is assumed to be linear.This is a common assumption when using strip theory code to solve the radiation and diffraction problems.In this context, it has been focused on representing the frequency dependent potential damping terms by a state-space representation such that on-line computations of the retardation functions can be avoided.The transformations for motions and potential coefficients, between different coordinate systems and origins, have been put into a rigid framework using vectorial mechanics.The result is a 6 DOF vectorial representation of the equations of motion.
The potential coefficients (added mass and potential damping) as well as the exciting wave loads (Froude-Krylov and diffraction forces) have been computed using ShipX (VERES) by Marintek for the S-175 container ship.The numerical results for this ship with discussions are presented at the end of the paper.

Longitudinal Model
The longitudinal terms can be represented in matrix form according to:

Lateral Model
Where the added mass terms are given by:

Figure 1 .
Figure 1.Definitions of coordinate origins: W (water line), G (centre of gravity), and O (equations of motion).The h-frame is located in W and the b-frame is located in O.The variables LCG, VCG, and T are outputs from VERES while LCO and VCO are user inputs.

Figure 2 .
Figure 2. Definitions of ship motions in the b-frame.

2 )Ω
Body-fixed frame (b-frame): The b-frame X b Y b Z b is fixed to the hull, see Figure 2. The coordinate origin is denoted O and is located on the center line a distance LCO relative to L pp /2 (positive backwards) and a distance VCO relative to the baseline (positive upwards).The center of gravity G with respect to O is located at rb g ó[x g , y g , z g ] T while the h-frame origin W with respect to O is located at rb w ó[x w , y w , z w ] T .The X b -axis is positive toward the bow and the X h -axis is parallel to the mean X b -axis, the Y b -axis is positive towards starboard, and the Z b -axis is positive downward as shown in Figure 2. Consequently, the body-fixed b-frame carries out oscillations around the steadily translating hframe.The b-frame linear velocities vb o ó[u, v, w] T in O and angular velocities b bn product î is defined in terms of the matrix S(rb o ) é SS(3) (skewsymmetric matrix of order 3) such that: since T\ T óT.This follows directly from the results presented in Section 2.1.3.Notice that the principal rotations equation (6) for an angle n gives a diagonal Rh data and D strip theory program can be used to compute frequency dependent added mass A*(u), potential damping B*(u), and the generalized Froude-Krylov * FK and diffraction forces * D .These terms are all computed in the hydrodynamic reference frame.The generalized radiation-induced forces then becomes: * R óñA*(u) ¨*ñB*(u) ˙*ñg*( *) where g*( *) is the restoring forces and moments.This gives the h-frame representation: óJ*(d ) T g*( *)

Figure 3 .
Figure 3. Flow chart showing the numerical computations.

Figure 4 .
Figure 4. Marine Systems Simulator (MSS): 6 DOF equations of motion including wave excitation forces represented in Simulink for real-time simulation.

Figure 5 .
Figure 5. Experimental data showing the high-frequency approximation and viscous ramp for B 33 (u).

Figure 6 .
Figure 6.The retardation functions K ii are computed from equation (79) using trapezoidal integrations and plotted as a function of time t.The corresponding state-space models (5thorder) are plotted on top.This clearly shows that the state-space model accurately describe the retardation functions.

Figure 7 .
Figure 7. Body plan of the S-175 container ship.

Figure 8 .
Figure 8. Frequency dependent added mass A 11 , A 22 , A 33 and potential damping B 11 , B 22 , B 33 for the S-175 container ship.Circles indicates VERES data points while the solid line is due to interpolation.For the B ii -data the high-frequency approximation b ii /u3 is applied.In addition a viscous ramp function is added to the B ii -plots.For surge, A 11 and B 11 , added mass is chosen as 10% of the mass while B 11 simply is a viscous ramp.Added resistance data can further be used to improve damping in surge.

Figure 9 .
Figure 9. Frequency dependent added mass A 44 , A 55 , A 66 and potential damping B 44 , B 55 , B 66 for the S-175 container ship.Circles indicates VERES data points while the solid line is due to interpolation.For the B ii -data the high-frequency approximation b ii /u3 is applied.In addition a viscous ramp function is added to the B ii -plots.For roll, the viscous effect to due bilge keels (and possible anti-roll tanks) is included in the B 44 -plot.

Figure 10 .
Figure 10.Diffraction and Froude-Krylov Forces and moments in 6 DOF versus time simulated in Simulink using VERES data tables.

Table 2 .
S-T-F Strip Theory Coefficients.