7 Nonlinear Equations of Motion

7.1 Limitations of Linear Equations of Motion

Although linearized equations of motion are useful in analyzing controls fixed stability and turning ability of straight line stable ships, they have several limitations. Linearized equations of motion:

cannot be used to analyze the turning ability of directionally unstable ships
cannot accurately predict the behavior of vessels in tight maneuvers
cannot capture speed loss in a turn that can be up to \(40\%\) for high speed destroyers
cannot include coupling effects between speed and steering equations

Nonlinear models are needed to overcome the above limitations. Typically this is done by retaining the higher order terms of the Taylor expansion of hydrodynamic forces that were previously neglected. Popular approaches include terms up to second and third order of the Taylor series expansion.

For most commercial ships, Taylor series beyond third order are generally not considered as they do not contribute significantly to the dynamics in practice. For defense vehicles, it is common to include upto fourth or fifth order terms to capture the intricate dynamics.

7.2 Nonlinear Equations of Motion

The nonlinear speed and steering equations are given by \(\eqref{eq-nonlinear-eom-1}\), \(\eqref{eq-nonlinear-eom-2}\) and \(\eqref{eq-nonlinear-eom-3}\)

\[\begin{align} (m - X_{\dot{u}}) \dot{u} = f_1(u, v, r, \delta) \label{eq-nonlinear-eom-1} \end{align}\]

\[\begin{align} (m - Y_{\dot{v}}) \dot{v} - (Y_{\dot{r}} - m x_G) \dot{r} = f_2(u, v, r, \delta) \label{eq-nonlinear-eom-2} \end{align}\]

\[\begin{align} (I_z - N_{\dot{r}}) \dot{r} - (N_{\dot{v}} - m x_G) \dot{v} = f_3(u, v, r, \delta) \label{eq-nonlinear-eom-3} \end{align}\]

where \(f_1(u, v, r, \delta)\), \(f_2(u, v, r, \delta)\) and \(f_3(u, v, r, \delta)\) are shown in \(\eqref{eq-f1}\), \(\eqref{eq-f2}\) and \(\eqref{eq-f3}\) and contain all the terms upto the third order of the Taylor series expansion.

\[\begin{align} f_1(u, v, r, \delta) =& T - R + X_u (u - U) + X_v v + X_r r + X_{\delta} \delta + X_{uu} (u-U)^2 + X_{uv} (u-U)v \nonumber \\ &+ X_{ur} (u-U)r + X_{u\delta} u \delta + X_{vv} v^2 + X_{vr} vr + X_{v\delta} v \delta + X_{rr} r^2 \nonumber\\ &+ X_{r\delta} r \delta + X_{\delta \delta} \delta^2 + X_{uuu} (u - U)^3 + X_{uuv} (u - U)^2v + X_{uur} (u - U)^2r \nonumber \\ &+ X_{uu\delta} (u - U)^2\delta + X_{uvv} (u - U)v^2 + + X_{uvr} (u - U)vr \nonumber \\ &+ X_{uv\delta} (u - U)v\delta + X_{urr} (u - U)r^2 + X_{ur\delta} (u - U)r\delta \nonumber \\ &+ X_{u\delta \delta} (u - U)\delta^2 + X_{vvv} v^3 + X_{vvr} v^2r + X_{vv\delta} v^2\delta + X_{vrr} vr^2 \nonumber \\ &+ X_{vr\delta} vr\delta + X_{v\delta \delta} v\delta^2 + X_{rrr} r^3 + X_{rr\delta} r^2 \delta + X_{r\delta \delta} r \delta^2 \nonumber \\ &+ X_{\delta \delta \delta} \delta^3 + mvr + mx_G r^2 \label{eq-f1} \end{align}\]

\[\begin{align} f_2(u, v, r, \delta) =& Y_u (u - U) + Y_v v + Y_r r + Y_{\delta} \delta + Y_{uu} (u-U)^2 + Y_{uv} (u-U)v \nonumber \\ &+ Y_{ur} (u-U)r + Y_{u\delta} u \delta + Y_{vv} v^2 + Y_{vr} vr + Y_{v\delta} v \delta + Y_{rr} r^2 \nonumber\\ &+ Y_{r\delta} r \delta + Y_{\delta \delta} \delta^2 + Y_{uuu} (u - U)^3 + Y_{uuv} (u - U)^2v + Y_{uur} (u - U)^2r \nonumber \\ &+ Y_{uu\delta} (u - U)^2\delta + Y_{uvv} (u - U)v^2 + + Y_{uvr} (u - U)vr \nonumber \\ &+ Y_{uv\delta} (u - U)v\delta + Y_{urr} (u - U)r^2 + Y_{ur\delta} (u - U)r\delta \nonumber \\ &+ Y_{u\delta \delta} (u - U)\delta^2 + Y_{vvv} v^3 + Y_{vvr} v^2r + Y_{vv\delta} v^2\delta + Y_{vrr} vr^2 \nonumber \\ &+ Y_{vr\delta} vr\delta + Y_{v\delta \delta} v\delta^2 + Y_{rrr} r^3 + Y_{rr\delta} r^2 \delta + Y_{r\delta \delta} r \delta^2 \nonumber \\ &+ Y_{\delta \delta \delta} \delta^3 - mur + Y_{o} + Y_{ou} (u-U) + Y_{ouu} (u - U)^2 \label{eq-f2} \end{align}\]

\[\begin{align} f_3(u, v, r, \delta) =& N_u (u - U) + N_v v + N_r r + N_{\delta} \delta + N_{uu} (u-U)^2 + N_{uv} (u-U)v \nonumber \\ &+ N_{ur} (u-U)r + N_{u\delta} u \delta + N_{vv} v^2 + N_{vr} vr + N_{v\delta} v \delta + N_{rr} r^2 \nonumber\\ &+ N_{r\delta} r \delta + N_{\delta \delta} \delta^2 + N_{uuu} (u - U)^3 + N_{uuv} (u - U)^2v + N_{uur} (u - U)^2r \nonumber \\ &+ N_{uu\delta} (u - U)^2\delta + N_{uvv} (u - U)v^2 + + N_{uvr} (u - U)vr \nonumber \\ &+ N_{uv\delta} (u - U)v\delta + N_{urr} (u - U)r^2 + N_{ur\delta} (u - U)r\delta \nonumber \\ &+ N_{u\delta \delta} (u - U)\delta^2 + N_{vvv} v^3 + N_{vvr} v^2r + N_{vv\delta} v^2\delta + N_{vrr} vr^2 \nonumber \\ &+ N_{vr\delta} vr\delta + N_{v\delta \delta} v\delta^2 + N_{rrr} r^3 + N_{rr\delta} r^2 \delta + N_{r\delta \delta} r \delta^2 \nonumber \\ &+ N_{\delta \delta \delta} \delta^3 - m x_G ur + N_{o} + N_{ou} (u-U) + N_{ouu} (u - U)^2 \label{eq-f3} \end{align}\]

Note that two subscript term refer to a second derivative term divided by \(2!\) and a three subscript term refers to a third order derivative term divided by \(3!\). The Taylor coefficient \(\frac{1}{n!}\) for the \(n^{th}\) order term is included into the definition of the hydrodynamic coefficient. Thus, \(Y_{uv}\) represents \(\frac{1}{2!}\frac{\partial^2Y}{\partial u \partial v}\) and \(Y_{uvr}\) represents \(\frac{1}{3!}\frac{\partial^3Y}{\partial u \partial v \partial r}\). Note that in the above formulation cross coupling of velocity and acceleration terms are not included. In \(\eqref{eq-f1}\), \(T\) represents the thrust from the propeller and \(R\) denotes the resistance of the hull. \(Y_{o}\), \(Y_{ou}\) and \(Y_{ouu}\) in \(\eqref{eq-f2}\) and \(N_{o}\), \(N_{ou}\) and \(N_{ouu}\) in \(\eqref{eq-f3}\) denote the influence of propeller on sway force and yaw moment respectively.

Exploiting the geometric symmetry of a ship between the port and starboard side, some of the coefficients can be neglected. For example, in the sway equation, it can be observed that \(Y_u = Y_{uu} = Y_{uuu} = 0\) as a change in surge speed alone should not cause any sway force. Similarly, all the coefficients of terms with even powers \(v\), \(r\) and \(\delta\) in sway and yaw equation must be zero. For example, \(Y_{\delta \delta} = 0\) since the sway force due to the rudder must be anti-symmetric about \(\delta = 0\). However, for the surge force all even powers of \(v\), \(r\) and \(\delta\) are retained and all odd powers are neglected. For example, \(X_{vvv} = 0\) as the surge force must develop in the same direction equally for port and starboard sway motions. Neglecting the corresponding hydrodynamic coefficients as discussed here yields the simplified expressions of \(f_1(u, v, r, \delta)\), \(f_2(u, v, r, \delta)\) and \(f_3(u, v, r, \delta)\) as shown in \(\eqref{eq-f1-new}\), \(\eqref{eq-f2-new}\) and \(\eqref{eq-f3-new}\)

\[\begin{align} f_1(u, v, r, \delta) =& T - R + X_u (u - U) + X_{uu} (u-U)^2 \nonumber \\ &+ X_{vv} v^2 + X_{vr} vr + X_{v\delta} v \delta + X_{rr} r^2 \nonumber\\ &+ X_{r\delta} r \delta + X_{\delta \delta} \delta^2 + X_{uuu} (u - U)^3 \nonumber \\ &+ X_{uvv} (u - U)v^2 + X_{uvr} (u - U)vr \nonumber \\ &+ X_{uv\delta} (u - U)v\delta + X_{urr} (u - U)r^2 + X_{ur\delta} (u - U)r\delta \nonumber \\ &+ X_{u\delta \delta} (u - U)\delta^2 + mvr + mx_G r^2 \label{eq-f1-new} \end{align}\]

\[\begin{align} f_2(u, v, r, \delta) =& Y_v v + Y_r r + Y_{\delta} \delta + Y_{uv} (u-U)v \nonumber \\ &+ Y_{ur} (u-U)r + Y_{u\delta} u \delta \nonumber\\ &+ Y_{uuv} (u - U)^2v + Y_{uur} (u - U)^2r \nonumber \\ &+ Y_{uu\delta} (u - U)^2\delta + Y_{vvv} v^3 + Y_{vvr} v^2r + Y_{vv\delta} v^2\delta + Y_{vrr} vr^2 \nonumber \\ &+ Y_{vr\delta} vr\delta + Y_{v\delta \delta} v\delta^2 + Y_{rrr} r^3 + Y_{rr\delta} r^2 \delta + Y_{r\delta \delta} r \delta^2 \nonumber \\ &+ Y_{\delta \delta \delta} \delta^3 - mur + Y_{o} + Y_{ou} (u-U) + Y_{ouu} (u - U)^2 \label{eq-f2-new} \end{align}\]

\[\begin{align} f_3(u, v, r, \delta) =& N_v v + N_r r + N_{\delta} \delta + N_{uv} (u-U)v \nonumber \\ &+ N_{ur} (u-U)r + N_{u\delta} u \delta \nonumber\\ &+ N_{uuv} (u - U)^2v + N_{uur} (u - U)^2r \nonumber \\ &+ N_{uu\delta} (u - U)^2\delta + N_{vvv} v^3 + N_{vvr} v^2r + N_{vv\delta} v^2\delta + N_{vrr} vr^2 \nonumber \\ &+ N_{vr\delta} vr\delta + N_{v\delta \delta} v\delta^2 + N_{rrr} r^3 + N_{rr\delta} r^2 \delta + N_{r\delta \delta} r \delta^2 \nonumber \\ &+ N_{\delta \delta \delta} \delta^3 - m x_G ur + N_{o} + N_{ou} (u-U) + N_{ouu} (u - U)^2 \label{eq-f3-new} \end{align}\]

Notice that the coefficients of the higher order terms in the above equations are known as hydrodynamic coefficients instead of hydrodynamic derivatives (as was done for their linear counterparts). This is because the higher order coefficients are obtained by curve fitting of the data rather than by evaluating slopes. For example, \(Y_{vvv}\) is evaluated by curve fitting the \(Y\) vs \(v\) data rather than evaluating the third order slope of the \(Y\) vs \(v\) curve at \(v=0\).

The nonlinear speed and steering equations along with the kinematic equations can be rearranged to express in the state space form (discussed in detail in later chapters) as shown in \(\eqref{eq-state-space-neom}\).

\[\begin{align} \dot{u} &= \frac{f_1(u, v, r, \delta)}{(m - X_{\dot{u}})} \nonumber \\ \dot{v} &= \frac{(I_z - N_{\dot{r}})f_2(u, v, r, \delta) + (Y_{\dot{r}} - m x_G)f_3(u, v, r, \delta)}{(m - Y_{\dot{v}})(I_z - N_{\dot{r}}) - (N_{\dot{v}} - m x_G)(Y_{\dot{v}} - m x_G)} \nonumber \\ \dot{r} &= \frac{(N_{\dot{v}} - m x_G)f_2(u, v, r, \delta) + (m - Y_{\dot{v}})f_3(u, v, r, \delta)}{(m - Y_{\dot{v}})(I_z - N_{\dot{r}}) - (N_{\dot{v}} - m x_G)(Y_{\dot{v}} - m x_G)} \label{eq-state-space-neom}\\ \dot{x}_{0O} &= u \cos(\psi) - v \sin(\psi) \nonumber \\ \dot{y}_{0O} &= u \sin(\psi) + v \cos(\psi) \nonumber \\ \dot{\psi} &= r \end{align}\]

7.3 Solving Nonlinear Equations of Motion

Together these coupled six equations can be solved to obtain the complete solution of the maneuvering motions in calm water for specified control input \(\delta\). As analytical solutions may not exist in most cases, these equations are usually solved numerically.

In general these equations can be solved using implicit or explicit schemes. Explicit methods calculate the state of a system at a later time from the state of the system at the current time, while implicit methods find a solution by solving an equation involving both the current state of the system and the later one. Explicit schemes are fast as the next states are evaluated purely from the knowledge of the current states. Implicit schemes are more computationally expensive as a matrix inversion is usually required to solve for the next states. For schemes achieving the same order of accuracy, implicit schemes are numerically more stable than explicit schemes. Thus, the time step \(h\) can be larger in implicit schemes as compared to explicit schemes.

7.3.1 Euler Scheme

Considering the state space \(\{x\} = [u \quad v \quad r \quad x_{0O} \quad y_{0O} \quad \psi]^T\), the equations of motion depicted in \(\eqref{eq-state-space-neom}\) can be expressed as

\[\begin{align} \{\dot{x}\} = f(t, \{x\}) \end{align}\]

Here, \(\{x\}\) is the unknown function of time. Euler scheme is the simplest explicit scheme where the next state is obtained from the current state by extrapolation along the slope calculated at the current state. Discretizing time \(t_{n+1} = t_{n} + h\) for \(n = 0, 1, 2, ...\), the solution \(\{x\}_{n+1}\) from Euler scheme is given by \(\eqref{eq-euler-scheme}\).

\[\begin{align} \{x\}_{n+1} = \{x\}_{n} + h f\left(t_n, \{x\}_n\right) \label{eq-euler-scheme} \end{align}\]

Given the initial condition \(\{x(t_0)\} = \{x_0\}\), the state space at the next time steps can be evaluated using \(\eqref{eq-euler-scheme}\).

7.3.2 Runge-Kutta \(4^{th}\) order Scheme

One of the popular explicit schemes is the Runge-Kutta \(4^{th}\) order scheme. Discretizing time \(t_{n+1} = t_{n} + h\) for \(n = 0, 1, 2, ...\), the solution \(\{x\}_{n+1}\) is given by

\[\begin{align} \{x\}_{n+1} = \{x\}_{n} + \frac{h}{6} \left(\{k_1\} + 2\{k_2\} + 2\{k_3\} + \{k_4\}\right) \label{eq-rk4} \end{align}\]

where

\[\begin{align} \{k_1\} &= f\left(t_n, \{x\}_n\right) \\ \{k_2\} &= f\left(t_n + \frac{h}{2}, \{x\}_n + \frac{h}{2} \{k_1\}\right) \\ \{k_3\} &= f\left(t_n + \frac{h}{2}, \{x\}_n + \frac{h}{2} \{k_2\}\right) \\ \{k_4\} &= f\left(t_n + h, \{x\}_n + h \{k_3\}\right) \end{align}\]

Thus, given the initial condition \(\{x(t_0)\} = \{x_0\}\), the state space at the next time steps can be evaluated using \(\eqref{eq-rk4}\). This scheme too can be thought of as the extrapolation of current state along a slope to reach the next state. However, unlike the Euler scheme that estimates the slope only at the current state, Runge-Kutta scheme estimates the slope as a weighted average of slopes evaluated at the beginning, midpoint and end of the interval.

There are variants of this scheme that adopt an adaptive time step too to increase the speed (read more here). This adaptive time version is available in most standard programming languages such as MATLAB or Python and can be readily used. In MATLAB, the Runge-Kutta scheme is provided by the function ode45. In Python, the corresponding function is solve_ivp provided through the scipy library.

In the explicit methods as discussed above, the required level of accuracy can be achieved using a sufficiently small time step \(h\). The Euler scheme with a well chosen time step is good enough to simulate the typical maneuvers of a ship. This is due to the fact that the derivatives \(\dot{u}\), \(\dot{v}\) and \(\dot{r}\) typically vary slowly with time due to the large inertia of the ships.

7.4 Memory Effects

Notice that in the Taylor expansion considered for the hydrodynamic forces, a dependence on time was not assumed. This means that the hydrodynamic forces and moments at an instant depend only on the current values of the motion, velocity and acceleration of the vessel and does not retain a memory of the past. Such a process is referred to as a Markov process. Note that in general hydrodynamic forces do not conform to the Markov assumption as waves generated (at the air water interface) at a current time do have effect on its dynamics at a later time. If the motions of the vehicle are sufficiently slow, the memory effects can be disregarded. Since maneuvering motion time scales are much longer than seakeeping motion time scales (occurring due to the presence of waves), the Markov assumption can be made for the maneuvering motions.

Notice that this is also the reason that the frequency dependence is disregarded for the added mass terms (\(X_{\dot{u}}\), \(Y_{\dot{v}}\), \(Y_{\dot{r}}\), \(N_{\dot{v}}\) and \(N_{\dot{r}}\)) and the damping terms (\(X_{u}\), \(Y_{v}\), \(Y_{r}\), \(N_{v}\) and \(N_{r}\)). The frequency dependent added mass terms can be calculated at the zero frequency using frequency domain panel method codes like WAMIT, ANSYS AQWA, Nemoh (open source), HydRA (free to use) etc. Due to this slow motion assumption, the hydrodynamic derivatives are also sometimes referred to as slow derivatives.