Further Study on One of the Numerical Methods for Pure Loss of Stability in Stern Quartering Waves

Abstract

The International Maritime Organization (IMO) finalized the second-generation intact stability criteria in 2022. However, an accurate and practical numerical method for stability loss has yet to be established. Therefore, a 6 DOF numerical model is further improved based on the previous study. Firstly, the rolling motion is simulated using a seakeeping model instead of the previous maneuvering mathematical model. Secondly, the roll-restoring variation is calculated directly considering the instantaneous wet hull instead of the previous pre-calculated method. Thirdly, transferring frequency to time is used to obtain heave and pitch motions, further considering yaw angle and sway velocity. Fourthly, the dynamic forces for sway, roll, and yaw motions are calculated, further considering the effect of the speed variation. Fifthly, the 6 DOF motions are used to determine the instantaneous wet hull, and the FK force and the hydrostatic force are calculated by the body’s exact method. Finally, a new conclusion is obtained that the sway and yaw motions’ effect on the ship speed loss, the relative longitudinal wave profile by the speed loss, the rudder angles, and the accompanying rudder forces in the rolling direction are significant, and much more than their centrifugal force or coupled force in the rolling direction.

1. Introduction

Dynamic stability failure modes could happen when the ship sails in the natural seas. The environment of the natural sea is complicated, and the existing intact stability code cannot guarantee the safety of the ship in the natural sea. For improving the safety level of the ship in the actual seas, the interim guideline on the second-generation intact stability criteria was approved by the International Maritime Organization (IMO) on 10 December 2020 [1], and her explanatory notes were finalized by the IMO in 2022 [2]. Five stability failure modes, including dead ship condition, pure loss of stability, parametric rolling, surf-riding/broaching, and excessive acceleration, three levels, including vulnerable Level 1, Level 2, direct assessment, and operational guidance, are considered in the stability criteria. An accurate numerical method is crucial for predicting stability failure modes directly. However, the physics of the failure modes had not yet been well understood, especially since there is still a dispute on the definition of pure loss of stability in stern quartering waves. An accurate and practical numerical method for stability loss has yet to be well established. Therefore, a mathematical model established previously by the authors [3] is further improved for calculating stability loss directly.

In the early stage, pure loss of stability is considered as the loss of a static roll-restoring arm in waves at the crest by some pioneer researchers, such as Paulling [4,5], Hamamoto and Nomoto [6], and Kuo and Vassalos [7], and their studies focused on how to calculate the roll-restoring arm in waves. Therefore, the vulnerable Level 1 and Level 2 criteria of pure loss of stability are drafted based on the roll-restoring arm in waves.

The natural sea is a random wave, and the stochastic process of stability loss in irregular waves is studied by Dunwoody [8], and Bulian and Francescutto [9]. The methods of failure probability are studied at this stage. The Grim effective wave is used by Umeda and Yamakoshi [10], a narrow-banded stochastic process is studied by Vermeer [11], the “critical wave” is used by Themelis and Spyrou [12], the up-crossing theory is used by Bulian et al. [13], and the split-time method is used by Belenky et al. [14].

With the development of hydrodynamics and the improvement of human cognition, stability loss is considered a capsizing mode of significant time-domain roll motion due to the lost roll-restoring arm at the crest and the staying time at the crest. Bassler pointed out that stability loss is not limited to the loss of the roll-restoring arm at the crest [15]. Some numerical methods for stability loss are reviewed by Neves [16]. A surge-roll coupled mathematical model for pure loss of stability in following seas was developed by Hashimoto [17], and the heave and pitch motions were further added by the authors, and a conclusion is confirmed that a small heeling angle should be set in the numerical simulations and the experiments for stability loss in following seas [18]. The time-domain simulation begins to be carried out in this stage.

With the improvement of human cognition, stability loss is considered not “pure” in astern waves. A new conclusion is pointed out by Umeda that the maneuvering sway and yaw motions could be the reason for the large rolling in stern quartering waves, and stability loss could not be pure [19,20]. Following that, their new mathematical model was submitted to IMO [21]. As a result, at least 4 DOF mathematical models are requested by IMO for predicting pure loss of stability in the time domain.

The significant rolling motion during stability loss is predicted by Liu et al. [22] using the CFD software/method, and their results are compared with the author’s experimental results in the reference [3]. However, the practical application of the CFD method in engineering is limited because a considerable time could be costly compared to the potential method. The vulnerable criteria produce conservative results, and their methods are much simpler than the direct numerical predicting method in the time domain. In addition, the vulnerable criteria of pure loss of stability are studied by some researchers [23,24,25].

The numerical methods without standard expressions could not be used in general. Therefore, a standard maneuvering prediction method is presented by Yasukawa and Yoshimura [26]. Some broaching prediction methods are given out by Hashimoto et al. [27] and Umeda et al. [28] based on the maneuvering prediction methods. Based on the maneuvering and broaching prediction methods, the first author tried to establish a surge-sway-roll-pitch (4 DOF) standard numerical method for predicting stability loss in astern waves [29,30]. As a result, a numerical method with surge-sway-heave-roll-pitch-yaw coupled is established. However, the roll-restoring arm in waves is pre-calculated, taking heave and pitch motions into account [3].

Stability loss is a complicated wave phenomenon, and nonlinear dynamics can be used to study its mechanics [31]. However, time-domain stability loss in astern seas is related to the ship’s maneuvering force, the rudder force, and the body’s exact roll-restoring force in waves. Therefore, it is still challenging to predict stability loss accurately in the time domain. Therefore, a 6 DOF numerical model is further established based on the previous research [3].

The improvement includes the following aspects: (1) The roll-restoring arm in waves is pre-calculated, taking heave and pitch motions into account in the previous research [3], and the mathematical model can be regarded as a 4 DOF model. In contrast, the roll-restoring arm in waves is calculated considering the instantaneous wet hull determined by 6 DOF ship motions, and the mathematical model can be regarded as a 6 DOF model. (2) The roll motion is simulated using a time-domain seakeeping model instead of the previous maneuvering mathematical mode. (3) The body’s exact FK and hydrostatic forces in the rolling direction are calculated directly, referring to the numerical method for parametric roll [32,33,34], and the body’s exact FK forces in the swaying and yawing direction are calculated directly. (4) The diffraction forces and radiation forces are calculated with a two-dimensional vectorizing method and an enhanced strip method at several constant speeds, and the speed-varied diffraction forces and radiation forces are considered by an interpolation method. (5) The speed loss due to the significant yaw motion is noted in this paper. (6) The roll-restoring arms’ comparison at both the starboard and port, and the different components of the restoring force variations GZ_W_FK + B, GZ_W_FK, and GZ_W_B are given out in this paper. (7) The numerical maximum roll angels at both the starboard and port are compared with the experimental results, and a good agreement is realized in this paper, while only the maximum roll angels at the port were given out with a large error compared with the experimental results in the previous research [3]. (8) The new phenomena of subharmonic rolling and yaw-roll coupling during pure loss of stability are defined in this paper. (9) A new conclusion about the sway and yaw motions’ effect on pure loss of stability is obtained.

2. Mathematical Model

2.1. Coordinate Systems

A horizontal body coordinate system was developed by Hamamoto and Kim [35] for studying stability in waves, and the coordinate system is shown in Figure 1. Their relations are shown in the authors’ previous reference [3]. The ship heading angle is ${\chi }_C$, and the heading angle of the incident wave $-{\chi }_C$ is added to the coordinate system. The instantaneous heading angle $\chi$ takes the yaw motion into account.

2.2. Mathematical Model

The ship motion is six degrees of freedom (DOF) in the real seas, and the mathematical model with different DOF is used for the engineering application. The mathematical model with 6 DOF is also used by many researchers. However, the dynamic stability failure models in waves are nonlinear phenomena. The nonlinear significant rolling motion during stability loss in astern waves is related to seakeeping, maneuvering, stability, thrust, and speed loss in waves, and few mathematical models can be utilized to predict stability loss perfectly. An improved 6 DOF mathematical model for predicting stability loss in astern waves is presented as follows. Equation (1) is the time-domain surge motion. The sway and yaw motions are maneuvering mathematical models for the course keeping, as shown in Equations (2) and (4). For using maneuvering hull forces as few as possible in the roll motion, Equation (3) is the time-domain roll motion with the seakeeping model instead of the previous maneuvering mathematical model in the author’s reference [3]. When a ship advances in astern waves with a high forward speed while the wavelength is near equal to the ship length, the encountered frequency could be very low, even zero. The radiation coefficients could be divergent at a low encountered frequency. Therefore, the heave and pitch motion with a low encountered frequency could be divergent. Therefore, the frequency-domain heave and pitch motions are pre-calculated by an enhanced strip method developed by Kashiwagi [36,37]. Transferring frequency to time is used to obtain the heave and pitch motions, as shown in Equations (5)–(8). The yaw motion is not considered in the frequency-domain heave and pitch motions. However, the instantaneous heading with yaw angle and sway velocity are considered in the time domain, as shown in Equation (9), which is an improvement compared with the authors’ previous study in the reference [3].

Equation (10) is the control equation for course keeping.

$$\left(m+A_{11}\right)\dot{u}-(m+A_{22}(u))vr=X_H+X_P+X_R(\delta )+F_1^{FK}({\xi }_G/\lambda ,\chi )$$

(1)

$$\left(m+A_{22}(u)\right)\dot{v}+\left(m+A_{11}\right)ur=Y_H+Y_R(\delta )+F_2^{FK}({\xi }_G/\lambda ,\zeta ,\phi ,\theta ,\chi )+F_2^{DF}({\xi }_G/\lambda ,{\chi }_C,u)$$

(2)

$$\begin{array}{l}(I_{xx}+A_{44})\stackrel{\cdot }{p}+D(p)+A_{42}(u)\stackrel{.}{v}+B_{42}(u)v+A_{46}(u)\stackrel{.}{r}+B_{46}r+C_{46}(\chi -{\chi }_C)=K_H+K_R(\delta )\\ +F_4^{FK+B}({\xi }_G/\lambda ,\zeta ,\phi ,\theta ,\chi )+F_4^{DF}({\xi }_G/\lambda ,{\chi }_C,u)\end{array}$$

(3)

$$\left(I_{zz}+A_{66}(u)\right)\dot{r}=N_H+N_R(\delta )+F_6^{FK}({\xi }_G/\lambda ,\zeta ,\phi ,\theta ,\chi )+F_6^{DF}({\xi }_G/\lambda ,{\chi }_C,u)$$

(4)

$$(m+A_{33}(u))\stackrel{\cdot \hspace{0.17em}\cdot }{\zeta }+B_{33}(u)\stackrel{\cdot }{\zeta }+C_{33}\zeta +A_{35}(u)\stackrel{\cdot \hspace{0.17em}\cdot }{\theta }+B_{35}(u)\stackrel{\cdot }{\theta }+C_{35}\theta =F_3^{FK}(u)+F_3^{DF}(u)$$

(5)

$$(I_{yy}+A_{55}(u))\stackrel{\cdot \hspace{0.17em}\cdot }{\theta }+B_{55}(u)\stackrel{\cdot }{\theta }+C_{55}\theta +A_{53}(u)\stackrel{.\hspace{0.17em}.}{\zeta }+B_{53}(u)\stackrel{\cdot }{\zeta }+C_{53}\zeta \hspace{0.17em}=F_5^{FK}(u)+F_5^{DF}(u)$$

(6)

$${\zeta }_G(t)={\zeta }_{Ga}(u)\cos [2\pi \cdot ({\xi }_G/\lambda )-{\delta }_3(u)]$$

(7)

$$\theta (t)={\theta }_a(u)\cos [2\pi \cdot ({\xi }_G/\lambda )-{\delta }_5(u)]$$

(8)

$${\xi }_G/\lambda =(u(t)\times \cos (\chi (t))-v(t)\sin (\chi (t))-C_W)\cdot t/\lambda$$

(9)

$$\dot{\delta }=\left\{-\delta -K_P(\chi -{\chi }_C)-K_P{T}_{D}r\right\}/T_E$$

(10)

The subscripts H, R, and P refer to the hull, rudder, and propeller, respectively. Note: the subscripts P in Equations (10) and (35) have different meanings.

The hydrodynamic force acting on a ship hull, the resistance in calm water, the propeller thrust, and the steering rudder forces are shown in the authors’ previous reference [3].

2.3. Excited Forces by the Wave

The body’s exact FK and hydrostatic forces in the roll direction should be taken into account for predicting parametric roll. The authors pre-calculate roll restoring variation for predicting stability loss in the previous study [3]. Here, the body’s exact FK and hydrostatic forces in the rolling direction are calculated directly. The body’s exact FK forces in the swaying and yawing direction are also calculated directly. In addition, the diffraction forces and radiation forces are calculated with a two-dimensional vectorizing method at several constant speeds, and the speed-varied diffraction forces and radiation forces are considered by an interpolation method. These are the improvements compared with the authors’ previous study in the reference [3].

A ship is advancing with a constant speed $u$ and oscillating with a circular frequency ${\omega }_e$ in deep water. With the assumption of linearized potential flow, the velocity potential and the excited wave forces $F_j^{FK}$ and $F_j^{DF}$ are given out by Kashiwagi et al. in the reference [36,37] and their seakeeping book according to the slender ship theory. The excited wave forces $F_j^{FK}$ and $F_j^{DF}$ are re-obtained detailed as follows, according to the reference [36,37]. $n_j$ denotes the $j-th$ component of the unit normal vector, and $S_H(x)$ is the sectional contour at station $x$. The method of subsection integration is used from step 2 to step 3 in Equation (13).

$$\begin{array}{l}{F}_j^{FK}({\xi }_G/\lambda ,u,\zeta ,\theta ,\chi )=\frac{\rho g{\zeta }_w}{i\omega }{\displaystyle {\int }_L\left\{{\displaystyle {\int }_{S_H}{n}_j}(i{\omega }_e-u\frac{\partial }{\partial x}){\phi }_{0}d\ell \right\}dx}\\ =\frac{\rho g{\zeta }_w}{i\omega }{\displaystyle {\int }_L\left\{{\displaystyle {\int }_{S_H}\left[(i{\omega }_e+iku\cos \beta )n_j{e}^{-kz-ik[x\cos \chi +y\sin (-\chi )]}\right]}d\ell \right\}dx}\\ =\rho g{\zeta }_w{\displaystyle {\int }_L{e}^{-ikx\cos \chi }\left\{{\displaystyle {\int }_{S_H}{e}^{{}^{-kz-ik[x\cos \chi +y\sin (-\chi )]}}{n}_j}d\ell \right\}dx}\hspace{1em}\hspace{1em}(j=2,3,4,5,6)\end{array}$$

(11)

$$\begin{array}{l}{F}_j^{FK}(t)=\left|F_j^{FK}({\xi }_G/\lambda ,u,\zeta ,\theta ,\phi ,\chi )\right|\cdot \cos \left\{2\pi \cdot ({\xi }_G/\lambda )-{\delta }_j^{FK}\right\}\hspace{1em}\hspace{1em}(j=2,3,4,5,6)\\ F_1^{FK}(t)=-\rho g{\zeta }_{w}k\cos \chi {\displaystyle {\int }_{AE}^{FE}{C}_1(x)S(x)e^{-kd(x)/2}\sin k({\xi }_G+x\cos \chi )dx}\\ C_1=\frac{\sin (k\sin \chi \cdot B(x)/2)}{k\sin \chi \cdot B(x)/2}\end{array}$$

(12)

$$\begin{array}{l}{F}_j^{DF}(u,\chi )=\frac{\rho g{\zeta }_w}{i\omega }{\displaystyle {\int }_L\left\{{\displaystyle {\int }_{S_H}{n}_j}(i{\omega }_e-u\frac{\partial }{\partial x}){\phi }_{7}d\ell \right\}dx}\hspace{1em}(j=2,3,4,5,6)\\ =\frac{\rho g{\zeta }_w}{i\omega }{\displaystyle {\int }_L\left\{{\displaystyle {\int }_{S_H}(}i{\omega }_e{n}_j{\phi }_7-u\cdot n_j\frac{\partial {\phi }_7}{\partial x})d\ell \right\}dx}\\ =\frac{\rho g{\zeta }_w}{i\omega }{\displaystyle {\int }_L\left\{{\displaystyle {\int }_{S_H}(}i{\omega }_e{n}_j{\phi }_7+u{\phi }_7\frac{\partial \hspace{0.17em}{n}_j}{\partial x})d\ell \right\}dx}\\ =\frac{\rho g{\zeta }_w}{i\omega }{\displaystyle {\int }_L\left\{{\displaystyle {\int }_{S_H}\left[i{\omega }_e{n}_j{k}_0{e}^{-k_{0}z-i{k}_{0}x\cos \chi }(i\sin (-\chi ){\phi }_2+{\phi }_3)\right]d\ell }\right\}dx}\\ +\frac{\rho g{\zeta }_w}{i\omega }{\displaystyle {\int }_L\left\{{\displaystyle {\int }_{S_H}\left[u{k}_0{e}^{-k_{0}z-i{k}_{0}x\cos \chi }(i\sin (-\chi ){\phi }_2+{\phi }_3)\frac{\partial \hspace{0.17em}{n}_j}{\partial x}\right]}d\ell \right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2+{\phi }_3)(n_j+\frac{u}{i{\omega }_e}\frac{\partial \hspace{0.17em}{n}_j}{\partial x}\right]d\ell }\right\}dx}\hspace{1em}\end{array}$$

(13)

$$\begin{array}{l}{F}_2^{DF}(u,\chi )=-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2+{\phi }_3)(n_2+\frac{u}{i{\omega }_e}\frac{\partial \hspace{0.17em}{n}_2}{\partial x}\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2{n}_2+{\phi }_3{n}_2)\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{i\sin (-\chi )(-\rho {\displaystyle {\int }_{S_H}{\phi }_2{n}_{2}d\ell )+(-\rho {\displaystyle {\int }_{S_H}{\phi }_3{n}_{2}d\ell })}\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{i\sin (-\chi )\left[A_{22}(x)+\frac{1}{i{\omega }_e}{B}_{22}(x)\right]+\left[A_{23}(x)+\frac{1}{i{\omega }_e}{B}_{23}(x)\right]\right\}dx}\end{array}$$

(14)

$$\begin{array}{l}{F}_3^{DF}(u,\chi )=-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2+{\phi }_3)(n_3+\frac{u}{i{\omega }_e}\frac{\partial \hspace{0.17em}{n}_3}{\partial x}\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2{n}_3+{\phi }_3{n}_3)\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}{\phi }_3{n}_{3}d\ell +i\sin (-\chi )(-\rho {\displaystyle {\int }_{S_H}{\phi }_2{n}_{3}d\ell )}}\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{\left[A_{33}(x)+\frac{1}{i{\omega }_e}{B}_{33}(x)\right]+i\sin (-\chi )\left[A_{23}(x)+\frac{1}{i{\omega }_e}{B}_{23}(x)\right]\right\}dx}\end{array}$$

(15)

$$\begin{array}{l}{F}_4^{DF}(u,\chi )=-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2+{\phi }_3)(n_4+\frac{u}{i{\omega }_e}\frac{\partial \hspace{0.17em}{n}_4}{\partial x}\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2{n}_4+{\phi }_3{n}_4)\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{i\sin (-\chi )(-\rho {\displaystyle {\int }_{S_H}{\phi }_2{n}_{4}d\ell )+(-\rho {\displaystyle {\int }_{S_H}{\phi }_3{n}_{4}d\ell })}\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{i\sin (-\chi )\left[A_{24}(x)+\frac{1}{i{\omega }_e}{B}_{24}(x)\right]+\left[A_{34}(x)+\frac{1}{i{\omega }_e}{B}_{34}(x)\right]\right\}dx}\end{array}$$

(16)

$$\begin{array}{l}{F}_5^{DF}(u,\chi )=-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2+{\phi }_3)(n_5+\frac{u}{i{\omega }_e}\frac{\partial \hspace{0.17em}{n}_5}{\partial x}\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2+{\phi }_3)(-x{n}_3-\frac{u}{i{\omega }_e}{n}_3)\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(-1)(x+\frac{u}{i{\omega }_e})(i\sin (-\chi ){\phi }_2+{\phi }_3)n_3\right]d\ell }\right\}dx}\\ ={\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }(x+\frac{u}{i{\omega }_e})\left\{-\rho {\displaystyle {\int }_{S_H}{\phi }_3{n}_{3}d\ell +i\sin (-\chi )(-\rho {\displaystyle {\int }_{S_H}{\phi }_2{n}_{3}d\ell )}}\right\}dx}\\ ={\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }(x+\frac{u}{i{\omega }_e})\left\{\left[A_{33}(x)+\frac{1}{i{\omega }_e}{B}_{33}(x)\right]+i\sin (-\chi )\left[A_{23}(x)+\frac{1}{i{\omega }_e}{B}_{23}(x)\right]\right\}dx}\end{array}$$

(17)

$$\begin{array}{l}{F}_6^{DF}(u,\chi )=-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2+{\phi }_3)(n_6+\frac{u}{i{\omega }_e}\frac{\partial \hspace{0.17em}{n}_6}{\partial x}\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(i\sin (-\chi ){\phi }_2+{\phi }_3)(x{n}_2+\frac{u}{i{\omega }_e}{n}_2)\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }\left\{-\rho {\displaystyle {\int }_{S_H}\left[(x+\frac{u}{i{\omega }_e})(i\sin (-\chi ){\phi }_2+{\phi }_3)n_2\right]d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }(x+\frac{u}{i{\omega }_e})\left\{i\sin (-\chi )(-\rho {\displaystyle {\int }_{S_H}{\phi }_2{n}_{2}d\ell )}-\rho {\displaystyle {\int }_{S_H}{\phi }_3{n}_{2}d\ell }\right\}dx}\\ =-{\zeta }_w\omega {\omega }_e{\displaystyle {\int }_L{e}^{-k_{0}z-i{k}_{0}x\cos \chi }(x+\frac{u}{i{\omega }_e})\left\{i\sin (-\chi )\left[A_{22}(x)+\frac{1}{i{\omega }_e}{B}_{22}(x)\right]+\left[A_{32}(x)+\frac{1}{i{\omega }_e}{B}_{32}(x)\right]\right\}dx}\end{array}$$

(18)

The lift force of the after and forward sections is also considered based on the strip theory.

2.4. Roll Restoring Force Variation

The roll-restoring arm GZ_W_FK + B, GZ_W_FK, and GZ_W_B based on the Froude–Krylov assumption can be defined as follows.

$$\begin{array}{l}G{Z}_W\_FK+B=-F_4^{FK+B}({\xi }_G/\lambda ,\zeta ,\phi ,\theta ,\chi )/W\\ G{Z}_W\_FK\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=-F_4^{FK}({\xi }_G/\lambda ,\zeta ,\phi ,\theta ,\chi )/W\\ G{Z}_W\_B\hspace{1em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=-F_4^B({\xi }_G/\lambda ,\zeta ,\phi ,\theta ,\chi )/W\end{array}$$

(19)

The similar symbols of GZ used in this paper are defined as follows.

$$\begin{array}{l}GZ\_Sway=(A_{42}(u)\stackrel{.}{v}+B_{42}(u)v)/W\\ GZ\_Yaw=(A_{46}(u)\stackrel{.}{r}+B_{46}r+C_{46}(\chi -{\chi }_C))/W\\ GZ\_KR=(-K_R(\delta ))/W\end{array}$$

(20)

2.5. Roll Damping Moment

The roll damping moment is calculated by Equation (21) with the linear and cubic roll damping coefficients denoted by $\alpha$ and $\gamma$.

$$D(p)=(I_{xx}+A_{44})(2\alpha \cdot p+\gamma \cdot p^3)$$

(21)

The roll damping coefficients are obtained by carrying out the rolling decay test.

$$\alpha =\frac{2A}{T_{\phi }}=\frac{A}{\pi }\sqrt{\frac{m\cdot g\cdot GM}{I_{xx}+A_{44}}}$$

(22)

$$\gamma =C\cdot \frac{4T_{\phi }}{3{\pi }^2}{(\frac{180}{\pi })}^2=C\frac{8}{3\pi }\sqrt{\frac{I_{xx}+A_{44}}{m\cdot g\cdot GM}}{(\frac{180}{\pi })}^2$$

(23)

Takahashi’s method [38] is used to take the effect of speed as follows.

$$A=2.51995\times (1+0.8(1-e^{-10Fn}))$$

(24)

$$C=0.0005096\times (1+0.8(1-e^{-10Fn}))$$

(25)

2.6. Hull Forces in Still Water

Although the expressions of hull forces in still water are generally used for predicting ship maneuvering and broaching, for the system of the mathematical model in this paper, the expressions of the hull forces in still water $X_H,Y_H,K_H$, and $\hspace{0.17em}{N}_H$ are rewritten as follows by referring [26,28]. The unified expressions of the nondimensional maneuvering coefficients are rewritten by the authors in the reference [3], and the symbols of maneuvering coefficients are also shown in the reference [3]. The effect of the high-order coefficients on the surge motion is investigated by the authors, and its effect is generally small [3]. Therefore, only resistance is remained in $X_H$. The seakeeping model of the roll motion is used in this paper, and only the heel-induced hydrodynamic force is considered in $K_H$.

$$X_H=-R(u)$$

(26)

$$\begin{array}{l}{Y}_H=\frac{1}{2}\rho L_{pp}d{U}^2(Y_v^{\prime }\cdot v^{\prime }+Y_r^{\prime }\cdot r^{\prime }+Y_{\phi }^{\prime }\cdot \phi \\ \hspace{1em}\hspace{1em}+Y_{vvv}^{\prime }\cdot v^{\prime }{}^3+Y_{vvr}^{\prime }\cdot v^{\prime }{}^2{r}^{\prime }+Y_{vrr}^{\prime }\cdot v^{\prime }{r}^{\prime }{}^2+Y_{rrr}^{\prime }\cdot r^{\prime }{}^3)\end{array}$$

(27)

$$K_H=\frac{1}{2}\rho L_{pp}{d}^2{U}^2(K_{\phi }^{\prime }+Y_{\phi }^{\prime }\frac{OG}{d})\cdot \phi$$

(28)

$$\begin{array}{l}{N}_H=\frac{1}{2}\rho L_{pp}^{2}d{U}^2(N_v^{\prime }\cdot v^{\prime }+N_r^{\prime }\cdot r^{\prime }+N_{\phi }^{\prime }\cdot \phi \\ \hspace{1em}\hspace{1em}+N_{vvv}^{\prime }\cdot v^{\prime }{}^3+N_{vvr}^{\prime }\cdot v^{\prime }{}^2{r}^{\prime }+N_{vrr}^{\prime }\cdot v^{\prime }{r}^{\prime }{}^2+N_{rrr}^{\prime }\cdot r^{\prime }{}^3)\end{array}$$

(29)

Nondimensional sway velocity $v\prime$, and yaw rate $r\prime$ are expressed as follows:

$$v^{\prime }=\frac{v}{U}\hspace{0.17em},\hspace{1em}{r}^{\prime }=\frac{r{L}_{pp}}{U}\hspace{0.17em}$$

(30)

2.7. Rudder Force

Rudder forces are essential parts of the mathematical model for pure loss of stability in stern quartering waves, and the rudder force in the rolling direction is focused on in this paper. Although the expressions of rudder forces are generally used for predicting ship maneuvering and broaching, for the system of the mathematical model in this paper, the expressions of rudder forces $X_R,Y_R,\hspace{0.17em}{N}_R$, and $\hspace{0.17em}{K}_R$ with twin rudders are rewritten as follows by referring [39]. The average values of $\stackrel{\_\_}{{\gamma }_R}$, $\stackrel{\_\_}{{\ell }_R^{\prime }}$ are used in this paper. (S: starboard; P: port)

$$X_R=-(1-t_R)F_N\sin \delta$$

(31)

$$Y_R=-(1+a_H)F_N\cos \delta$$

(32)

$$N_R=-(x_R+a_H{x}_{HR})F_N\cos \delta$$

(33)

$$K_R=(z_R+a_H{z}_{HR})F_N\cos \delta$$

(34)

$$F_N=\frac{1}{2}\rho (A_{R(S)}+A_{R(P)})(u_R^2+v_R^2)f_{\alpha }\sin {\alpha }_R$$

(35)

$$\begin{array}{l}{u}_R=\epsilon \left(1-w_p\right)u\\ \hspace{1em}\hspace{0.17em}\sqrt{\eta {\left\{1+\kappa \left(\sqrt{1+\frac{8K_T(J_P)}{\pi J_P^2}}-1\right)\right\}}^2+1-\eta }\end{array}$$

(36)

$$v_R=\stackrel{\_\_}{{\gamma }_R}U(\beta -\stackrel{\_\_}{{\ell }_R^{\prime }}{r}^{\prime })$$

(37)

$${\alpha }_R=\delta -\stackrel{\_\_}{{\gamma }_R}\frac{U}{u_R}(\beta -\stackrel{\_\_}{{\ell }_R^{\prime }}{r}^{\prime })$$

(38)

$$\stackrel{\_\_}{{\gamma }_R}=\frac{1}{2}({\gamma }_{R(s)}+{\gamma }_{R(P)})$$

(39)

$$\stackrel{\_\_}{{\ell }_R^{\prime }}=\frac{1}{2}({\ell }_R^{\prime }{}_{(s)}+{\ell }_R^{\prime }{}_{(P)})$$

(40)

$$f_{\alpha }=\frac{6.13\mathsf{\Lambda }}{2.25+\mathsf{\Lambda }}\hspace{1em}$$

(41)

$$\epsilon =\frac{1-w_R}{1-w_P}$$

(42)

$$\eta =\frac{D_P}{H_R}\hspace{0.17em}$$

(43)

$$\beta =\arctan (\frac{-v}{u})\hspace{0.17em}$$

(44)

$$U=\sqrt{u^2+v^2}$$

(45)

2.8. Propeller Thrust and the Resistance in Still Water

The propeller thrust and the resistance in still water are used for time-domain surge motion, and the propeller thrust $X_P\hspace{0.17em}$ with twin propellers, and the hull resistance in still water $R\hspace{0.17em}$ is also given out as follows for the system of the mathematical model in this paper.

$$X_P=2\times (1-t_P)T$$

(46)

$$T=\rho n_P{}^2{D}_P^4{K}_T(J_P)$$

(47)

$$J_P=\frac{(1-w_P)\hspace{0.17em}u}{n_P{D}_P}$$

(48)

$$R=\frac{1}{2}\rho S_F{u}^2{C}_T(\frac{u}{\sqrt{g{L}_{PP}}})$$

(49)

3. Experiments and Subject Ship

The ONR (The Office of Naval Research in USA) tumblehome hull used in the authors’ previous reference [3] is still used in this paper. The principal particulars and the system parameters of the ONR tumblehome are given out in

Table 1. Principal particulars of the ONR tumblehome.

Items	Ship	Model
Length: LPP	154.0 m	3.800 m
Draft: d	5.494 m	0.136 m
Breadth: B	18.8 m	0.463 m
Depth: D	14.5 m	0.358 m
Displ.: W	8507 ton	127.8 kg
CB	0.535	0.535
GM	1.48 m	0.037 m
OG	−2.729 m	−0.067 m
LCB	−2.569 m	−0.063 m
Tφ	14.0 s	2.199 s
κyy/LPP	0.25	0.25
Κzz/LPP	0.25	0.25
$2\times A_R$	2 × 23.74 m2	2 × 0.0145 m2
${\delta }_{\max }$	35 degrees	35 degrees

and

Table 2. The parameters used in this mathematical model.

Items	Values	Items	Values
$m_x/m$	0.015	$t_R$	0.300
$K_P$	3.000	$a_H$	0.312
$T_D$	0.000	$x_R+a_H{x}_{HR}$	−87.4
$T_E$	0.100	$z_R+a_H{z}_{HR}$	−8.740
$t_P$	0.250	$\kappa$	0.350
$w_P$	0.150	$\epsilon$	1.170
$z_H/d$	0.852	$\stackrel{\_\_}{{\gamma }_R}$	0.494
$\mathsf{\Lambda }$	1.180	$\stackrel{\_\_}{{\ell }_R^{\prime }}$	−0.378
		$\eta$	0.982

.

The free-running experiment was conducted at China Ship Scientific Research Center, one member of the International Towing Tank Conference (ITTC) organization. The experiment of pure loss of stability was carried out in the seakeeping basin (length: 69 m, breadth: 46 m, depth: 4 m) according to the test guidelines in the ITTC.

The ship model was free running in astern waves with twin propellers and twin rudders. The 6 DOF optical fiber gyroscope was placed on the ship model to measure the roll, pitch, and yaw angles during the free running. An onboard system was connected to an onshore control computer by wireless to record the roll, pitch, yaw, rudder angles, and propeller rotation speed. A servo-needle wave height sensor was attached to a steel bridge which is 78 m in length and spans over the basin to measure the wave elevation at the middle of the basin. A total station system is used to record the position of the ship model with a specified propeller rate in calm water, and then one specified propeller rate corresponding to one nominal speed in calm water is obtained. The same specified propeller rate in calm water is used to achieve the same nominal speed in the free-running experiments in astern waves.

First, the model is kept with a zero forward speed and an initial heading angle near the wave maker manually by two workmen. Next, the waves are generated by the wave-making system. Then, the propeller rate is increased up to the specified value after receiving the order from the onshore control computer. Finally, the ship model is released free near one wave crest with its initial heading when the wave train propagates far enough, and then the model automatically runs in astern waves with its specified propeller rate and autopilot course.

A PD control system is used for course keeping which reacts according to the bias between the instantaneous heading angle measured by the gyroscope and the autopilot course, and the yaw velocity measured by the gyroscope. The rudder gain is set by the experience according to the reaction of course keeping in calm water.

The detail of the procedure of the experiment can be referred to the reference [3,18,40].

4. Simulations and Discussions

The roll-restoring variation becomes significant when the wavelength is nearly equal to the ship length, and the surging amplitude becomes significant when the wavelength is larger than the ship length in astern waves. One wave condition is selected for study stability loss in astern waves, which is the same as that in the reference [3]. The wave condition is H/Lpp = 0.05 and λ/Lpp = 1.25. The autopilot heading is set with χ_C = 30 degrees. The following sections will not mention the wave condition repeatedly. Some other results are also submitted to IMO by the first author as a member of IMO [41,42,43,44,45]. However, the opinions expressed in this paper only reflect the academic research and do not reflect the views of China or UK.

4.1. Roll-Restoring Arm

The roll-restoring arms based on the Froude–Krylov assumption are calculated considering the instantaneous wet hull. The roll-restoring arms in the port are further studied compared with the authors’ previous study [3]. GZ_still is the righting arm in still water, GZ_W_FK + B_static is the roll-restoring arms in waves with a static balance method, and GZ_W_FK + B_Fn is the roll-restoring arms in waves with an enhanced strip method at different constant speeds. Here, the heeling angles of 10 degrees (starboard) and −10 degrees (port)are set, respectively. The GZ loss in the starboard is significant when the midship section is located on the crest, as shown in Figure 2a. In contrast, the GZ loss in the port is significant if the midship section is on the upslope, as shown in Figure 2b. This is one of the reasons why stability loss occurs on the upslope of the waves. The roll-restoring arms in the time domain will be further discussed in the next part.

The ship velocity without/with the yaw motion and wave celerity are shown in Figure 3a,b. The wave celerity is much faster than the ship forward speed. The actual forward speed is periodically varied with a significant amplitude around a mean speed when the surge motion is considered. The mean speed is nearly equal to the nominal speed without the yaw motion or with a small yaw motion, as shown in Figure 3a. In contrast, the average speed is significantly smaller than the speed in calm water with a large yaw motion at this condition, as given out in Figure 3b, because the effective thrusts in the heading direction become small due to the significant yaw motion. The discussion on speed loss due to the significant yaw motion is an improvement compared with the authors’ previous study [3].

GZ_still and GZ_W_FK + B with and without the surge are shown in Figure 4a,b. The different components of the roll-restoring arms GZ_W_FK and GZ_W_B are shown in Figure 5a,b. Here, the heeling angles of 10 and −10 degrees are set, respectively. The retention time at the crest could be extended by the surge motion. Therefore, the retention time of stability loss could be extended at the crest at the starboard. Figure 5b shows that the stability loss at the port is significant at the upslope of the waves. The hydrostatic force GZ_W_B is near the same at the starboard and port in astern seas with a heading of 30 degrees, as given out in Figure 5a, while the wave exciting FK force GZ_W_FK is significantly different at the starboard and port in astern seas with a heading of 30 degrees, as given out in Figure 5b. As we all know, the roll-restoring arms, including GZ_W_B, GZ_W_FK, and GZ_W_FK + B, should be the same at the starboard and port in following seas because the ship hull is symmetric to the center line, and the roll-restoring arm in waves becomes small at the crest and becomes large at the trough when the wavelength is near equal to the ship length in following seas because the roll-restoring arm is related with the area moment of waterplane and the waterplane becomes small at the crest and becomes large at the trough. However, the wave exciting FK force GZ_W_FK could be different in oblique waves. The main reason why the roll variation GZ_W_FK + B at the starboard in Figure 4a is significantly different from that at the port in Figure 4b is that the wave exciting FK force GZ_W_FK is significantly different at the starboard and port in astern seas with a heading of 30 degrees, as given out in Figure 5b. The comparison of the roll-restoring arms at both the starboard and port, and the different components of the restoring force variations GZ_W_FK + B, GZ_W_FK, and GZ_W_B are given out in this section which is an improvement compared with the authors’ previous study [3].

4.2. Contributions from Sway, Yaw, and Rudder

The contributions of some key elements in the rolling mathematical model to the roll moment are shown in Figure 6, Figure 7 and Figure 8 by comparing with the roll-restoring arms in waves GZ_W_FK + B. The results are simulated with constant heeling angles of 20 and −20 degrees, respectively. As shown in Figure 6, we can find that the rudder produces a significant negative contribution in the rolling direction. One of the reasons is that the value of $z_R+a_H{z}_{HR}$ is large due to a high gravity center with a small GM, as shown in Table 2. The second reason is that the rudder’s normal force $F_N$ is large with a high inflow velocity at a high forward speed. The third reason is that the yaw motion becomes large, and the rudder angle becomes essential to keep the course. The coupled sway force has a slightly negative contribution in the rolling direction, as given out in Figure 7, but the coupled yawing forces are not tiny, as given out in Figure 8, but their negative and positive contribution is nearly the same.

The comparisons of maximum roll angles during stability loss between mathematical models with and without the sway and yaw motions are conducted. The 6 DOF numerical method can predict the roll angle, as given out in Figure 9a. In comparison, the 4 DOF numerical method underestimates the roll angle without the sway and yaw motions. This is because the 4 DOF method, ignoring the sway and yaw motions, means the rudder force is also ignored with zero rudder angle at a constant heading angle. However, the rudder force is a crucial factor for stability loss at high ship speeds, as given out in Figure 6. The effect of the coupled sway and yaw force on the maximum roll angles during stability loss is generally small, as shown in Figure 10a,b. Because the coupled sway and yaw force in the rolling direction is insignificant, as shown in Figure 7 and Figure 8. That is to say, the sway and yaw motions’ effect on stability loss comes from the rudder force because the yaw motion and the rudder angle are significant. The rudder angle is larger than 15 degrees to keep the course, as shown in Figure 11. The significant yaw motion also affects the surge force and the thrust in the heading direction, and the speed lost is significant, as given out in Figure 3b. Then the relative wave profiles are affected.

Therefore, the mathematical model should take the sway and yaw motions into account. Their effects on the ship speed loss, the relative longitudinal wave profile due to the speed loss, the rudder angle, and accompanying rudder forces in the rolling direction are significant and much more than their coupled or centrifugal force in the rolling direction.

4.3. New Phenomenon during Pure Loss of Stability

The time-domain rolling, pitching, yawing motions, and rudder angle in the tests are analyzed. Figure 11a shows a periodic stable rolling motion, while Figure 11b shows an unstable rolling motion. Figure 11c shows that the fourth and the seventh roll periods are two times the encounter period, which means subharmonic rolling occurs during stability loss. This phenomenon was also observed by Kan et al. [46] and named period bifurcation. The significant roll developments toward the port occur with the large starboard yaw angle in Figure 11d. The ship heading angle is 43 degrees, and the rudder angle is the maximum of 35 degrees to keep course during the ship capsizing, as shown in Figure 11d, and this could be a typical pattern of capsizing due to yaw-roll coupling with loss of stability.

The 6 DOF time-domain ship motions are calculated with the numerical method mentioned in this paper. The 6 DOF numerical method can predict the maximum rolling angles, as given out in Figure 12. However, the new phenomenon during stability loss, such as unstable rolling motion and subharmonic rolling, still cannot be captured by the numerical method.

The above discussion considers the rolling motion in the starboard and port with the improved 6 DOF numerical method, which is also an improvement, while the previous study [3] only considers the port.

5. Conclusions

The mathematical model of stability loss is updated based on the authors’ previous mathematical model, and the following remarks can be made:

(1)

The actual forward speed should be considered by the surge motion.

(2)

The yaw motion is significant and results in a ship speed loss which affects the relative longitudinal wave profile in the time domain. The rudder angles are substantial due to the considerable yaw motion, and then the rudder forces in the rolling direction are significant. The above effects from the coupled sway and yaw motion are much more than their centrifugal or coupled force in the rolling direction.

(3)

This paper’s improved 6 DOF mathematical model can be utilized for stability loss in astern waves at this stage.

As unstable rolling motions, subharmonic rolling and yaw-roll coupling could exist during stability loss, and some research could be carried out on these complicated phenomena.