Performance Simulation of Marine Cycloidal Propellers: A Both Theoretical and Heuristic Approach

Abstract

The importance of mathematical and numerical simulation in marine engineering is growing together with the complexity of the designed systems. In general, simulation a makes it possible to improve the engineering design, reducing working time and costs of production as well. In this respect, the implementation of a simulation model for cycloidal propellers is presented. Cycloidal thrusters are being increasingly used in marine applications. Their best performance concerns low-speed applications, due to their ability to steer thrust in any direction. The proposed simulator is able to assess the performance of cycloidal propellers in terms of the generated thrust and torque, without resorting to consuming and demanding computational tools, such as CFD methods. This feature makes the presented model particularly suitable for the simulation in the time domain of the maneuverability of surface units, equipped with cycloidal propellers. In this regard, after embodying the implemented model in an already existing simulation platform for maneuverability, we show the most significant outputs concerning some simulated maneuvers, performed at cruise speed.

1. Introduction

Until the past century, the only way that naval architects had to predict the behavior of the system (intended as the ship or a part of it, such as the propulsion plant, or the auxiliary systems) they were working on, was to build a prototype and test it. Nowadays, thanks to the development and the progress of computer science, it is possible to shape not prototypes but simulators, based on mathematical laws, that can predict in advance the behavior not only of a single subsystem, but also of the whole ship. Indeed, one of the main advantages of mathematical and numerical simulation is the possibility to compare different design choices, so improving the engineering design and reducing working time and costs as well.

For example, making use of mathematical models, the hull performance can be analyzed under any weather conditions [1,2,3], it is possible to assess whether the designed machinery can guarantee the needed power [4], the magnitude and the direction of propulsion and steering forces can be predicted [5,6], or any kind of failure conditions can be analyzed [7,8,9].

Among the several thruster types, cycloidal propellers (CPs) are widely used on water tractors, ferries and some naval vessels. CPs are made of a set of vertical blade protruding from the hull and performing two main rotations: one around the rotor axis and one around the axis passing through the blade pivoting centre. Depending on their eccentricity value e, namely the ratio between the distance of the steering centre from the propeller axis and the radius of the rotor, they can be classified into true cycloidal ($\mathit{e}=1$), epicycloidal ($\mathit{e}<1$) and trochoidal propellers ($\mathit{e}>1$). Being able to generate almost the same thrust in every direction and combine both thrust and steering in a singel unit, CPs are very suitable for low-speed applications, such as dynamic positioning operations for instance.

To date, some investigations have been performed in order to model the behavior of such devices. In [10,11,12,13], some performances of cycloidal rudders have been shown, while propeller–hull interactions have been presented in [14], mainly in connection with maneuverability aspects. In this paper, the main features of a time domain simulation model for CPs are presented. The model is developed on a Matlab©Simulink platform and is able to calculate the time histories of the provided thrust and torque. The implemented model relies on a mixture of theoretical and empirical considerations. In particular, the propeller thrust and torque evaluation is based on the kinematics of each single blade, taking into account suitable correction factors in order to properly consider “dissipative” phenomena (such as interference between blades, the shielding induced by the half of the rotor which receives the oncoming flow, and the slight reduction of back thrust). The calibration of the simulator is carried out by comparing simulation outputs with real data found in open source. The final result is a simulation platform able to predict the performance of a CP in terms of generated thrust and torque, avoiding consuming and demanding computations, such as CFD methods. This is one of the significant aspects of the developed simulator that makes its use really effective when integrated into a platform for the simulation of ship motions.

In this regard, a first application of the developed simulator has been the validation of different thrust allocation logics of a DP system for a surface vessel equipped with a bow thruster and two cycloidal propellers at stern [15,16]. After that, in order to assess the reliability of the obtained simulation model, some maneuvers at cruise speed have been performed, embodying the CP simulator into a more complex dynamic model for maneuverability. In the following, the simulation outputs for one of these maneuvers are presented.

2. The Simulation Platform

As mentioned above, the present work is based on the simulation model which has been presented in [17,18] and concerns a 80 m long patrol vessel from the Italian Coast Guard. Such a simulation platform has been developed modular in order to be able to dealt with and possibly replace separately different sub-systems.

In modeling vessel dynamics, a crucial aspect is developing suitable mathematical models that reproduce the forces acting on the ship, as accurately as possible. In such a process, a complex system of coupled time-domain equations determines the evolution in time of each quantity which contributes to the ship dynamics. Usually, simulation platform are represented by means of flow–charts. In this work, the modular concept of the adopted platform is described by the mind map drawn in Figure 1, where the main part is represented by the motion equations of the ship (1) that mutually interact with the other sub-systems describing all the forces acting on the ship itself. Such forces mainly concern the interaction of the hull with the propulsion and steering systems as well as the environmental disturbs.

Whenever dealing with maneuvering problems, it is common to introduce two reference frames: the Earth-fixed reference frame $\{\mathrm{\Omega },{\underline{n}}_1,{\underline{n}}_2,{\underline{n}}_3\}$ and the body-fixed frame $\{O,{\underline{b}}_1,{\underline{b}}_2,{\underline{b}}_3\}$. Choosing the origin O as located on the mean water-free surface at midship, the main equations governing the ship motion are expressed as

$$\left\{\begin{array}{cc}\Delta (\dot{u}-x_G{r}^2-uv))\hfill & =X\hfill \\ \Delta (\dot{v}-x_G\dot{r}+ur))\hfill & =Y\hfill \\ I_z\dot{r}+m{x}_G(\dot{v}+ru)\hfill & =N\hfill \end{array}\right.$$

(1)

where $v_O=u{\underline{b}}_1+v{\underline{b}}_2$ denotes the linear velocity of O expressed in the body-fixed basis and $\omega =r{\underline{b}}_3$ is the angular velocity, $\Delta$ is the vessel displacement, $x_G$ is the longitudinal coordinate of gravity center w.r.t. $\{O,{\underline{b}}_1,{\underline{b}}_2,{\underline{b}}_3\}$, $\Delta$ is the vessel mass, $I_z$ is the moment of inertia about ${\underline{b}}_3$-axis, $\underline{R}=X{\underline{b}}_1+Y{\underline{b}}_2$ and $\underline{M}=N{\underline{b}}_3$ are the resultants of forces and the moments expressed in the $\underline{b}$-basis, respectively. In addition, $X=X_H+X_P+X_E$, $Y=Y_H+Y_P+Y_E$, $N=N_H+N_P+N_E$, where subscripts H, P, and E refer to hull, propellers, and environmental forces and moments respectively.

3. Cycloidal Propeller Model

Cycloidal propellers allows precise and stepless thrust generation since propulsion and steering forces can be generated and varied simultaneously. As a result of the rotation around its vertical axis, the same amount of thrust can be provided almost over ${360}^{\circ }$ by blades with hydrodynamically shaped profiles that assure a high degree of efficiency. In this section, a detailed description of the mathematical model developed for the computation of the thrust T and the torque $Q_P$ delivered by a CP is presented.

3.1. Blade Motion

The rapid and precise thrust variation of CPs is based on the kinematics of the blades (usually from 4 to 6 and equally spaced from each other) that move along a circular path, centered at the rotor center, and at the same time perform a superimposed pivoting motion around a suitable vertical axis. When the steering center overlaps the center of the rotor casing, the blades are not angled with respect to the tangent to the blade circular trajectory and no thrust is originated in this circumstance. Instead, if the steering center is moved away from the center of the rotor casing, the blades are set at a variable angle with respect to the tangent of their circular path, and thrust is generated. From top to bottom, Figure 2 shows an example of the installation of two CPs on a ship, with some details of the machinery. During the revolution motion, the maximum angle reached by the blades increases with the eccentricity e, defined as:

$$\mathit{e}=\frac{OC}{D/2}$$

(2)

where $OC$ is the distance between the center of the rotor O and the steering center C, and $D/2$ is the radius of the rotor. The motion of the pivot point (assumed as the center of mass) of the blade, relative to a stationary observer, results from the superimposition of the rotational movement of the rotor casing along a straight line representing the forward motion of the vessel. The pivot point follows the curve of a cycloid. The rolling radius of the cycloid is $D/2$ and the advance coefficient $\mathit{\lambda }$ is:

$$\mathit{\lambda }=\frac{V_A}{\pi nD}$$

(3)

where $V_A$ is the advance velocity and n is the rotor speed. During one revolution, the propeller travels a distance $\mathit{\lambda }D\pi$ in the direction of the vessel motion.

To generate thrust, the blades are angled with respect to the circular path described by their pivoting point. To achieve this, the steering center is moved from O to C. The resulting angle of attack leads to the generation of hydrodynamic lift and drag forces on each blade. The thrust provided by the propeller is the sum of such hydrodynamic forces, is always perpendicular to the line $OC$ and its intensity increases with the distance $OC$. By shifting the steering center C, it is possible to produce thrust in any direction and of different intensities. Therefore, the thrust provided by an epicycloidal propeller can be represented as a function of two plane polar coordinates:

• the geometric or driving pitch (between 0 and $0.8R$ for constructive limits): that is the distance (expressed as a percentage of the rotor radius $R=D/2$) between the steering center C and the center of the rotor O;
• the steering pitch (between $0^{\circ }$ and ${360}^{\circ }$): the angle between a fixed axis (with respect to the hull) and the line $OC$.

3.2. Thrust Generation

Changing the steering center, the blades acquire a certain attack angle, so generating corresponding lift and drag forces which give rise to the desired thrust. The hydrodynamic forces components acting transversely to the desired thrust direction cancel each other out. It is possible to produce thrust in any direction putting the steering center in the right position. The zero-thrust condition can be selected at any time, making the ship very safe to handle.

Each blade generates instantly a hydrodynamic force which is the sum of the lift (component of the hydrodynamic force, perpendicular to the oncoming flow) and the drag (parallel to the oncoming flow). The sum of all the hydrodynamic forces generated by all the blades gives rise to the corresponding total thrust. Also, each blade generates a corresponding torque which contributes to the total torque $\underline{M}$ acting on the propeller. For each choice of driving and steering pitch, there are corresponding curves of thrust and torque coefficients $K_S$ and $K_D$ as functions of the advance coefficient $\mathit{\lambda }$. The coefficients $K_S$ and $K_D$ are defined in analogy with the corresponding screw-propeller coefficients, respectively $K_T$ and $K_Q$, by the formulae reported in

Table 1. Cycloidal and screw propeller non-dimensional coefficients.

Coefficients	Cycloidal	Screw
advance coefficient	$\mathit{\lambda }=\frac{V_A}{\pi nD}$	$J=\frac{V_A}{nD}$
thrust coefficient	$K_S=\frac{T}{\frac{1}{2}\rho DL{\mathit{u}}^2}$	$K_T=\frac{T}{\rho n^2{D}^4}$
torque coefficient	$K_D=\frac{4M}{\rho D^{2}L{u}^2}$	$K_Q=\frac{Q_O}{\rho n^2{D}^5}$
efficiency	${\eta }_O=\frac{K_S}{K_D}\mathit{\lambda }$	${\eta }_O=\frac{K_T}{K_Q}\frac{J}{2\pi }$

where ${\mathit{V}}_{\mathit{A}}$ is the advance velocity, T is the propeller thrust, M and ${\mathit{Q}}_{\mathit{O}}$ are the CP and screw propeller torque respectively, L is the blade length, D is the rotor diameter, $\mathit{\rho }$ is the sea water density, L is the blade length, u is the tangential speed (u = nπD).

.

3.3. Kinematics of the Blade

In this section, the kinematical model of a blade is presented. In particular, a $2-$dimensional plane model has been adopted, propellers have been modeled as counter-rotating and two distinguished coordinate systems have been introduced: the first one is the hull-fixed frame, while the second one $\{O,{\underline{e}}_1,{\underline{e}}_2,{\underline{e}}_3\}$ is rotated clockwise, about the vertical axis passing through O and parallel to ${\underline{b}}_3={\underline{e}}_3$, by an angle $\beta \in [0,2\pi ]$ (the steering pitch) which determines (the perpendicular of) the steering force direction. The angle $\beta$ is related to the rudder pitch of the cycloidal propeller. The steering center C lies on the straight line passing through O and parallel to ${\underline{e}}_2$. The relationship between the bases $\left\{{\underline{b}}_i\right\}$ and $\left\{{\underline{e}}_i\right\}$, in accordance with Figure 3, is expressed as

$$\left\{\begin{array}{cc}{\underline{e}}_1\hfill & = cos\beta {\underline{b}}_1+sin\beta {\underline{b}}_2\hfill \\ {\underline{e}}_2\hfill & =-sin\beta {\underline{b}}_1+cos\beta {\underline{b}}_2\hfill \\ {\underline{e}}_3\hfill & ={\underline{b}}_3\hfill \end{array}\right.$$

(4)

During the rotation, the projection P of the blade shaft on the plane $\langle O,{\underline{b}}_1,{\underline{b}}_2\rangle$ describes a circumference having center O e radius R coinciding with the rotor radius. Such a circumference is parameterized by

$$P\left(\theta \right):\left\{\begin{array}{cc}x\hfill & =Rcos\theta \hfill \\ y\hfill & =Rsin\theta \hfill \\ z\hfill & =0\hfill \end{array}\right.$$

(5)

where $\theta$ denotes the angle (function of time) describing the revolution motion of the blade. The unit vector $\underline{t}$ tangent to the circular path of P has components in the vector basis $\left\{{\underline{b}}_i\right\}$ of the form

$$\underline{t}\left(\theta \right):\left\{\begin{array}{cc}{t}_1\hfill & =-sin\theta \hfill \\ t_2\hfill & = cos\theta \hfill \\ t_3\hfill & =0\hfill \end{array}\right.$$

(6)

Introducing the vector

$$(C-O)=s{\underline{e}}_2=-ssin\beta {\underline{b}}_1+scos\beta {\underline{b}}_2, s\in [0,0.8R]$$

(7)

the vector joining the steering center C with the point P can be expressed as

$$(P-C)=(Rcos\theta +ssin\beta ){\underline{b}}_1+(Rsin\theta -scos\beta ){\underline{b}}_2$$

(8)

The variable s is usually called driving pitch and controls the magnitude of the thrust. The unit vector orthogonal to $(P-C)$ and belonging to the plane $\langle O,{\underline{b}}_1,{\underline{b}}_2\rangle$ identifies with the unit vector of the blade chord and it is given by

$$\frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}=\frac{(-Rsin\theta +scos\beta ){\underline{b}}_1+(Rcos\theta +ssin\beta ){\underline{b}}_2}{\sqrt{{(-Rsin\theta +scos\beta )}^2+{(Rcos\theta +ssin\beta )}^2}}$$

(9)

The pivoting motion of the blade around its own vertical axis can be described by the angle $\alpha$ (function of time) between the unit vectors $\underline{t}$ and $\frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}$. Due to the relation

$$cos\alpha =\frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}\underline{t}=\frac{R+ssin(\beta -\theta )}{\sqrt{{(-Rsin\theta +scos\beta )}^2+{(Rcos\theta +ssin\beta )}^2}}$$

(10)

where the dot denotes the usual scalar product between vectors, choosing anticlockwise the positive direction of rotation around the blade shaft, the pivoting angle $\alpha$ can be defined as

$$\alpha =\pm arccos\left(\frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}\underline{t}\right) \mathrm{where} \begin{array}{cc}+& \mathrm{if}cos(\theta -\beta )\ge 0\\ -& \mathrm{if}cos(\theta -\beta )<0\end{array}$$

(11)

The above outlined kinematical model can be summarized by the following figure.

Supposing now that the vessel is moving, let ${\underline{v}}_O=\widehat{u}{\underline{b}}_1+\widehat{v}{\underline{b}}_2$ be the velocity of O (with respect to an Earth-fixed frame) expressed in the hull-fixed basis. Denoting by

$${\underline{v}}_p^{\prime }=-R\dot{\theta }sin\theta {\underline{b}}_1+R\dot{\theta }cos\theta {\underline{b}}_2$$

(12)

the velocity of the point P with respect to the body-fixed frame, the velocity of P with respect to the Earth-fixed frame is given by

$${\underline{v}}_p={\underline{v}}_p^{\prime }+{\underline{v}}_O+\underline{\omega }\wedge (P-O)=\left[\widehat{u}-R(\dot{\theta }+\dot{\psi })sin\theta \right]{\underline{b}}_1+\left[\widehat{v}+R(\dot{\theta }+\dot{\psi })cos\theta \right]{\underline{b}}_2$$

(13)

where $\underline{\omega }=\dot{\psi }{\underline{b}}_3$ is the angular velocity of the vessel. The velocity of the oncoming flow experienced at P by a blade-fixed observer is then $-{\underline{v}}_P$; its unit vector $\underline{\widehat{t}}$ is expressed as

$$\underline{\widehat{t}}=-\frac{{\underline{v}}_P}{|{\underline{v}}_P|}=-\frac{\left[\widehat{u}-R(\dot{\theta }+\dot{\psi })sin\theta \right]{\underline{b}}_1+\left[\widehat{v}+R(\dot{\theta }+\dot{\psi })cos\theta \right]{\underline{b}}_2}{\sqrt{{\left[\widehat{u}-R(\dot{\theta }+\dot{\psi })sin\theta \right]}^2+{\left[\widehat{v}+R(\dot{\theta }+\dot{\psi })cos\theta \right]}^2}}$$

(14)

Making use of the unit vector $\underline{\widehat{t}}$ it is possible to characterize the attack angle of the incident flow as

$$\widehat{\alpha }=\pi -arccos\left[\frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}·\underline{\widehat{t}}\right]$$

(15)

3.4. Hydrodynamic Forces

Making use of some simplifying assumptions, a suitable model for evaluating the hydrodynamic forces generated by each blade is presented. It is supposed that the velocity of the incident flow is the same on the entire surface of the blade and coincides with $-{\underline{v}}_P$. Under such a condition, the lift and drag produced by each blade can be expressed as

$$\begin{array}{cc}\hfill \underline{L}& =C_L\frac{1}{2}{\rho }_{w}A{\left|{\underline{v}}_P\right|}^2\underline{\widehat{n}}\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill \underline{D}& =C_D\frac{1}{2}{\rho }_{w}A{\left|{\underline{v}}_P\right|}^2\underline{\widehat{t}}\hfill \end{array}$$

(16b)

where $C_L$ and $C_D$ are the lift and drag coefficients, respectively; ${\rho }_w$ is the sea water density; A is the blade lateral area; ${\underline{v}}_P$ is the oncoming flow velocity; $\underline{\widehat{t}}$ is the unit vector of the lift force (unit vector of the oncoming flow at P); and $\underline{\widehat{n}}$ is the unit vector of the drag force (perpendicular to $\underline{\widehat{t}}$).

The unit vector $\underline{\widehat{n}}$ can be determined by the following procedure, in which two main scenarios are distinguished:

• the attack angle $\widehat{\alpha }$ belongs to the interval $\left[0,\frac{\pi }{2}\right]$, namely the oncoming flow hits the blade from the front. In such a circumstance, the unit vector $\underline{\widehat{n}}$ is determined according to the requirements

  $$\underline{\widehat{n}}=\left\{\begin{array}{cc}{\underline{b}}_3\wedge \underline{\widehat{t}}\hfill & \mathrm{when} \underline{\widehat{t}}\wedge \frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}·{\underline{b}}_3>0\hfill \\ \hfill & \\ -{\underline{b}}_3\wedge \underline{\widehat{t}}\hfill & \mathrm{when} \underline{\widehat{t}}\wedge \frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}·{\underline{b}}_3<0\hfill \end{array}\right.$$

  (17)
• $\widehat{\alpha }\in \left[\frac{\pi }{2},\pi \right]$, the oncoming flow hits the blade from the back. In this case, $\underline{\widehat{n}}$ is singled out by the requests:

  $$\underline{\widehat{n}}=\left\{\begin{array}{cc}-{\underline{b}}_3\wedge \underline{\widehat{t}}\hfill & \mathrm{when} \underline{\widehat{t}}\wedge \frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}·{\underline{b}}_3>0\hfill \\ \hfill & \\ {\underline{b}}_3\wedge \underline{\widehat{t}}\hfill & \mathrm{when} \underline{\widehat{t}}\wedge \frac{{(P-C)}^{\perp }}{{|(P-C)}^{\perp }|}·{\underline{b}}_3<0\hfill \end{array}\right.$$

  (18)

As remaining particular cases, if $\widehat{\alpha }=0$ or $\widehat{\alpha }=\pi$ there is no lift while if $\widehat{\alpha }=\frac{\pi }{2}$ then $\underline{\widehat{n}}=\underline{\widehat{t}}$. The above described procedure allows us to determine the lift and drag provided by each single blade. The resultant hydrodynamic force generated by the epicycloidal propeller can be computed as the sum of all contributions given by each blade.

3.5. Torque Acting on the Rotor

Figure 4 shows the implemented model for the propulsion system. Input data of the model are the desired engine speed $RP{M}_{set}$, the desired propeller pitch s, and the desired thrust angle $\beta$. The output is the array ${\left[X_p,Y_p,N_p\right]}^{\top }$ composed by the longitudinal and lateral propeller forces and the resulting moment.

The dynamics of the shaft line is described by the equation:

$$\frac{\mathrm{d}n}{\mathrm{dt}}=\frac{1}{I_{tot}}\left(Q_{eng}-Q_p-Q_{fric}\right)$$

(19)

where n is the shaft speed; $I_{tot}$ is the total axial inertia taking into account: (i) engine, (ii) gears, (iii) shaft and (iv) propeller contributions; $Q_{eng}$ is the engine torque; $Q_{fric}$ represents frictions; and $Q_p$ is the propeller torque.

In order to compute the torque acting on the rotor, the Newton-Euler moments equations for each single blade and for the rotor are taken into account separately. Developed in the hull-fixed reference frame and with respect to the point O (center of the rotor), the Newton-Euler moments equation for each blade can be expressed as

$${\underline{M}}_O^H+{\underline{M}}_O^G+{\underline{M}}_O^R+{\underline{M}}_O^I=I_{G_B}^B\left(\underline{\dot{\omega }}\right)+\underline{\omega }\wedge I_{G_B}^B\left(\underline{\omega }\right)+m_B(G_B-O)\wedge {\underline{a}}_{G_B}$$

(20)

where ${\underline{M}}_O^H$, ${\underline{M}}_O^G$, ${\underline{M}}_O^R$, and ${\underline{M}}_O^I$ are the moments acting on the blade, respectively due to hydrodynamic, weight, reactive and inertial forces; $I_{G_B}^B$ is the inertia tensor of the blade w.r.t. its gravity center $G_B$ (which is assumed to coincide with the pivot point P); $\underline{\omega }=(\dot{\theta }-\dot{\alpha }){\underline{b}}_3$ is the blade angular velocity w.r.t. the hull-fixed frame; ${\underline{a}}_{G_B}$ is the acceleration of $G_B$ w.r.t. the hull-fixed frame; and $m_B$ is the blade mass. The moment of the hydrodynamic force ${\underline{M}}_O^H$ is given by:

$${\underline{M}}_O^H=(P-O)\wedge (\underline{L}+\underline{D})$$

(21)

where the hydrodynamic force is described in terms of lift and drag. Expressing all vectors in the basis $\left\{{\underline{b}}_i\right\}$ as $(\underline{L}+\underline{D})=f_1{\underline{b}}_1+f_2{\underline{b}}_2$ and $(P-O)=R{\underline{b}}_{1}cos\theta +R{\underline{b}}_{2}sin\theta$, one has:

$${\underline{M}}_O^H=(R{f}_{2}cos\theta -R{f}_{1}sin\theta ){\underline{b}}_3$$

(22)

The weight force moment is given by:

$${\underline{M}}_O^G=(P-O)\wedge m\underline{g}=m_{B}g(Rsin\theta {\underline{b}}_1-Rcos\theta {\underline{b}}_2)$$

(23)

where $\underline{g}=g{\underline{b}}_3$ is the gravity acceleration. In order to evaluate the inertial force moment ${\underline{M}}_O^I$, it is necessary to assess the dragging and Coriolis forces acting on the blade and their associated moments ${\underline{M}}_O^S$ and ${\underline{M}}_O^C$. After that, the total inertial forces moment can be expressed as the sum:

$${\underline{M}}_O^I={\underline{M}}_O^C+{\underline{M}}_O^S$$

(24)

The moments ${\underline{M}}_O^S$ and ${\underline{M}}_O^C$ are calculated by means of the well-known formulae:

$$\begin{array}{cc}\hfill {\underline{M}}_O^C& ={\underline{M}}_P^C+{\underline{F}}^C\wedge (O-P)\hfill \end{array}$$

(25a)

$$\begin{array}{cc}\hfill {\underline{M}}_O^S& ={\underline{M}}_P^S+{\underline{F}}^S\wedge (O-P)\hfill \end{array}$$

(25b)

where ${\underline{M}}_P^C$ and ${\underline{M}}_P^S$ are the moments, w.r.t. the pivot point P, of the Coriolis and the dragging forces acting on the blade, while ${\underline{F}}^C$ and ${\underline{F}}^S$ denote the resultants of the Coriolis and dragging forces. The Coriolis force acting on any single blade is given by:

$${\underline{F}}^C=-2{\int }_{\mathcal{B}}k{\underline{\omega }}_{\psi }\wedge {\underline{v}}_Q\mathrm{d}\tau$$

(26)

where ${\underline{\omega }}_{\psi }=\dot{\psi }{\underline{b}}_3$ is the angular velocity of the vessel, ${\underline{v}}_Q={\underline{v}}_P+\underline{\omega }\wedge (Q-P)$ (with $\underline{\omega }=(\dot{\theta }-\dot{\alpha }){\underline{b}}_3$) is the velocity of a generic point Q of the blade w.r.t. the hull-fixed frame and k is the mass density of the blade. After implementing calculations, one gets the expression:

$$\begin{array}{cc}\hfill {\underline{F}}^C& =-2m_B{\underline{\omega }}_{\psi }\wedge {\underline{v}}_P-2m_B{\underline{\omega }}_{\psi }\wedge \left[\underline{\omega }\wedge (G_B-P)\right]=-2m_B{\underline{\omega }}_{\psi }\wedge {\underline{v}}_{G_B}\hfill \\ \hfill & =-2m_B(\dot{\psi }{v}_{G_{B}1}{\underline{b}}_2-\dot{\psi }{v}_{G_{B}2}{\underline{b}}_1)\hfill \end{array}$$

(27)

where ${\underline{v}}_{G_B}=v_{G_{B}1}{\underline{b}}_1+v_{G_{B}2}{\underline{b}}_2$ denotes the velocity of the gravity center of the blade w.r.t. the hull-fixed frame. By definition, the moment ${\underline{M}}_P^C$ of the Coriolis force with respect to the pivot point P is given by:

$$\begin{array}{cc}\hfill {\underline{M}}_P^C& =-2{\int }_{\mathcal{B}}k(Q-P)\wedge [{\underline{\omega }}_{\psi }\wedge \left[{\underline{v}}_P+\underline{\omega }\wedge (Q-P)\right]]\mathrm{d}\tau \hfill \\ \hfill & =-2{\int }_{\mathcal{B}}k(Q-P)\wedge [{\underline{\omega }}_{\psi }\wedge {\underline{v}}_P]\mathrm{d}\tau -2{\int }_{\mathcal{B}}k(Q-P)\wedge [{\underline{\omega }}_{\psi }\wedge \left[\underline{\omega }\wedge (Q-P)\right]]\mathrm{d}\tau \hfill \end{array}$$

(28)

Making use of the results shown in [19], the moment ${\underline{M}}_P^C$ can be expressed as:

$${\underline{M}}_P^C=-2m_B(G_B-P)\wedge [{\underline{\omega }}_{\psi }\wedge {\underline{v}}_P]-\underline{\omega }\wedge I_P^B\left({\underline{\omega }}_{\psi }\right)-{\underline{\omega }}_{\psi }\wedge I_P^B\left(\underline{\omega }\right)+I_P^B(\underline{\omega }\wedge {\underline{\omega }}_{\psi })$$

(29)

where now $I_P^B$ denotes the inertia tensor of the blade w.r.t. the pivot point P. Expression (29) holds in general. In our case, since the considered mathematical model is two-dimensional and in view of the assumption $G_B\equiv P$, it is easily seen that all terms appearing in Equation (29) vanish, so having ${\underline{M}}_P^C=0$.

For what concerns the dragging force, by definition it reads as:

$$\begin{array}{cc}\hfill {\underline{F}}^s& =-{\int }_{\mathcal{B}}k[{\underline{a}}_O+{\underline{\omega }}_{\psi }\wedge [{\underline{\omega }}_{\psi }\wedge (Q-O)]+{\underline{\dot{\omega }}}_{\psi }\wedge (Q-O)]\mathrm{d}\tau \hfill \\ \hfill & =-m_B{\underline{a}}_O-m_B{\underline{\omega }}_{\psi }\wedge [{\underline{\omega }}_{\psi }\wedge (G_B-O)]-m_B{\underline{\dot{\omega }}}_{\psi }\wedge (G_B-O)\hfill \\ \hfill & =-m_B\left(a_{O1}{\underline{b}}_1+a_{O2}{\underline{b}}_2\right)+m_B{\underline{\omega }}_{\psi }\wedge \left(R\dot{\psi }sin\theta {\underline{b}}_1-R\dot{\psi }cos\theta {\underline{b}}_2\right)+m_B\left(R\ddot{\psi }sin\theta {\underline{b}}_1-R\ddot{\psi }cos\theta {\underline{b}}_2\right)\hfill \\ \hfill & =-m_B\left(a_{O1}{\underline{b}}_1+a_{O2}{\underline{b}}_2\right)+m_B{\dot{\psi }}^2\left(Rcos\theta {\underline{b}}_1+Rsin\theta {\underline{b}}_2\right)+m_{B}R\ddot{\psi }\left(sin\theta {\underline{b}}_1-cos\theta {\underline{b}}_2\right)\hfill \end{array}$$

(30)

where ${\underline{a}}_O=a_{O1}{\underline{b}}_1+a_{O2}{\underline{b}}_2$ is the acceleration of the rotor center O w.r.t. the Earth–fixed frame. The moment ${\underline{M}}_P^S$ of the dragging force w.r.t. the pivot point P is given by:

$$\begin{array}{cc}\hfill {\underline{M}}_P^S& ={\int }_{\mathcal{B}}k(Q-P)\wedge [{\underline{a}}_O+{\underline{\omega }}_{\psi }\wedge [{\underline{\omega }}_{\psi }\wedge (Q-O)]+{\underline{\dot{\omega }}}_{\psi }\wedge (Q-O)]\mathrm{d}\tau \hfill \\ \hfill & =-m_B\left(G_B-P\right)\wedge {\underline{a}}_O-{\int }_{\mathcal{B}}k(Q-P)\wedge [{\underline{\omega }}_{\psi }\wedge [{\underline{\omega }}_{\psi }\wedge \left[\left(Q-P\right)+\left(P-O\right)\right]\left]\right]\mathrm{d}\tau +\hfill \\ \hfill & -{\int }_{\mathcal{B}}k\left(Q-P\right)\wedge [{\underline{\dot{\omega }}}_{\psi }\wedge \left[\left(Q-P\right)+\left(P-O\right)\right]]\mathrm{d}\tau \hfill \\ \hfill & =-m_B\left(G_B-P\right)\wedge {\underline{a}}_O-{\underline{\omega }}_{\psi }\wedge I_B\left({\underline{\omega }}_{\psi }\right)-m_B\left(G_B-P\right)\wedge {\underline{\omega }}_{\psi }\wedge [{\underline{\omega }}_{\psi }\wedge \left(P-O\right)]+\hfill \\ \hfill & -I_B^P\left({\underline{\dot{\omega }}}_{\psi }\right)-m_B\left(G_B-P\right)\wedge [{\underline{\dot{\omega }}}_{\psi }\wedge (P-O)]\hfill \\ \hfill & =-m_B\left(G_B-P\right)\wedge [{\underline{a}}_O+{\underline{\omega }}_{\psi }\wedge [{\underline{\omega }}_{\psi }\wedge \left(P-O\right)]+{\underline{\dot{\omega }}}_{\psi }\wedge (P-O)]-{\underline{\omega }}_{\psi }\wedge I_P^B\left({\underline{\omega }}_{\psi }\right)-I_P^B\left({\underline{\dot{\omega }}}_{\psi }\right)\hfill \end{array}$$

(31)

Again, since the model is two-dimensional and $G_B\equiv P$, only the last term does not vanish, so yielding:

$${\underline{M}}_P^S=-I_P^B\left({\underline{\dot{\omega }}}_{\psi }\right)=-I_{33}\ddot{\psi }{\underline{b}}_3$$

(32)

where $I_{33}$ is the moment of inertia of the blade w.r.t. the vertical axis passing for P. In order to implement Equation (25), the terms ${\underline{F}}^S\wedge (O-P)$ and ${\underline{F}}^C\wedge (O-P)$ need to be calculated:

$$\begin{array}{cc}\hfill {\underline{F}}^S\wedge (O-P)& =mR\left(a_{O1}sin\theta -a_{O2}cos\theta \right){\underline{b}}_3-m{R}^2\ddot{\psi }{\underline{b}}_3\hfill \end{array}$$

(33a)

$$\begin{array}{cc}\hfill {\underline{F}}^C\wedge (O-P)& =-2mR\left(\dot{\psi }{v}_{G_{B}2}sin\theta +\dot{\psi }{v}_{G_{B}1}cos\theta \right){\underline{b}}_3\hfill \end{array}$$

(33b)

Now, inserting all the obtained results into Equation (24), we end up with the final expression of the moment w.r.t. O of the inertial forces:

$${\underline{M}}_O^I=-2mR\left(\dot{\psi }{v}_{G_{B}2}sin\theta +\dot{\psi }{v}_{G_{B}1}cos\theta \right){\underline{b}}_3-I_{33}\dot{\psi }{\underline{b}}_3+mR\left(a_{O1}sin\theta -a_{O2}cos\theta \right){\underline{b}}_3-m{R}^2\ddot{\psi }{\underline{b}}_3$$

(34)

For our two-dimensional model with $G_B=P$, the first two terms on the right side of (20) are:

$$\begin{array}{cc}\hfill I_{G_B}^B\left(\underline{\dot{\omega }}\right)& =I_{33}\left(\ddot{\theta }-\ddot{\alpha }\right){\underline{b}}_3\hfill \end{array}$$

(35a)

$$\begin{array}{cc}\hfill \underline{\omega }\wedge I_{G_B}^B\left(\underline{\omega }\right)& =0\hfill \end{array}$$

(35b)

$$\begin{array}{cc}\hfill m_B\left(G_B-O\right)\wedge a_{G_B}& =m_B\left(R{\underline{a}}_{G2}cos\theta -R{\underline{a}}_{G1}sin\theta \right){\underline{b}}_3\hfill \end{array}$$

(35c)

Inserting all the above calculated contributions into Equation (20), we obtain the explicit expression of the reactive moment:

$${\underline{M}}_O^R=-{\underline{M}}_O^H-{\underline{M}}_O^G-{\underline{M}}_O^I+I_{G_B}^B\left(\underline{\dot{\omega }}\right)+\underline{\omega }\wedge I_{G_B}^B\left(\underline{\omega }\right)+m_B\left(G_B-O\right)\wedge {\underline{a}}_{G_B}$$

(36)

Inserting the reactive moments (36) acting on each single blade into the moments equation for the rotor, the engine torque can be calculated as:

$$M_O^E=\sum _{i=1}^n{\left({\underline{M}}_O^R\right)}_i·{\underline{b}}_3-{\underline{M}}_O^I·{\underline{b}}_3+I_O^{rot}\left({\underline{\dot{\omega }}}_r\right)·{\underline{b}}_3$$

(37)

where ${\underline{M}}_O^I$ is the inertial forces moment acting on the rotor, $I_O^{rot}$ is the inertia tensor of the rotor w.r.t. its center O, ${\underline{\omega }}_r=\dot{\theta }{\underline{b}}_3$ is the angular velocity (w.r.t. the hull–fixed frame) of the rotor and n is the number of blades. The moment ${\underline{M}}_O^I$ can be calculated as already made for the blades, namely as the sum ${\underline{M}}_O^I={\underline{M}}_O^C+{\underline{M}}_O^S$. In this case, since the center of the rotor O is fixed w.r.t. the hull, the same arguments as in [20] can be applied so obtaining the general explicit expressions:

$${\underline{M}}_O^C=-{\underline{\omega }}_r\wedge I_O^{rot}\left({\underline{\omega }}_{\psi }\right)-{\underline{\omega }}_{\psi }\wedge I_O^{rot}\left({\underline{\omega }}_r\right)+I_O^{rot}({\underline{\omega }}_r\wedge {\underline{\omega }}_{\psi })$$

(38)

$${\underline{M}}_O^S=-m_{rot}(G_{rot}-O)\wedge {\underline{a}}_O-{\underline{\omega }}_{\psi }\wedge I_O^{rot}\left({\underline{\omega }}_{\psi }\right)-I_O^{rot}\left({\underline{\dot{\omega }}}_{\psi }\right)$$

(39)

with $I_O^{rot}$ the inertia tensor of the rotor w.r.t. its center O (which is assumed to coincide with the gravity center $G_{rot}$) and $m_{rot}$ is the mass of the rotor. Again, since the model is two-dimensional and $G_{rot}=O$, we have ${\underline{M}}_O^C=0$ and ${\underline{M}}_O^S=-I_O^{rot}\left({\underline{\dot{\omega }}}_{\psi }\right)$.

Inserting now all the obtained results into Equation (37), we obtain the final expression for the engine torque:

$$M_O^E=\sum _{i=1}^n{\left({\underline{M}}_O^R\right)}_i·{\underline{b}}_3+I_O^{rot}\left(\ddot{\psi }{\underline{b}}_3\right)·{\underline{b}}_3+I_O^{rot}\left(\ddot{\theta }{\underline{b}}_3\right)·{\underline{b}}_3$$

(40)

In conclusion, a detailed model for the evaluation of forces and torques acting on the CP has been expound. The relevance of the presented approach concerns the possibility to easily change the propeller characteristics (number, length and shape of blades as well as rotor diameter) and evaluate the corresponding performance variations.

4. Numerical Modeling and Validation: Free Running Test

The mathematical model illustrated in Section 3 has been used to develop a Matlab©Simulink simulator for cycloidal propellers. In this section, the main features and the validation of such simulator are presented.

In order to simplify the simulation platform, some hypotheses have been assumed: $\left(i\right)$ the propeller has been considered in free-running conditions; $\left(ii\right)$ the problem has been assumed stationary; $\left(iii\right)$ the model has been implemented two dimensional; $\left(iv\right)$ linear superposition of the contributions of each blade in terms of generated forces and moments has been adopted.

The propeller model input data are:

• the propeller geometry (length, chord and orbit diameter of the blade—see

  Table 2. Geometric parameters of the propeller.

  Parameter	Value
  Number of blades	5
  $\mathrm{Rotor}\mathrm{diameter} \left(\mathrm{m}\right)$	$3.2$
  $\mathrm{Blade}\mathrm{length} \left(\mathrm{m}\right)$	$2.65$
  $\mathrm{Blade}\mathrm{chord} \left(\mathrm{m}\right)$	$0.7744$
  $\mathrm{Maximum}\mathrm{tip}\mathrm{thickness} \left(\mathrm{m}\right)$	$0.242$

  );
• the sea water characteristics (viscosity and density);
• the lift $C_L$ and drag $C_D$ coefficients of the blade;
• the rotor speed and maximum available pitch;
• the steering pitch angle $\beta$ ($0^{\circ }$ in forward direction, ${180}^{\circ }$ in astern condition) and the driving pitch s (expressed as a percentage of the rotor radius).

In the present case study, data are:

The whole simulation model consists of a set of identical subsystems, each of them representing the behavior of a single blade. Making use of (16) and (40), the components of total thrust and torque are calculated in the basis $\left\{{\underline{b}}_i\right\}$.

Free-running characteristics have been evaluated through a simulation campaign where $K_S$ and $K_D$ coefficients have been evaluated, in accordance with Table 1, for several $\mathit{\lambda }\in (0,0.6)$ and $s\in (30\%,90\%)$ values. In particular, the evaluation of coefficients $K_S$ and $K_D$ has been performed in the pitch range from $30\%$ to $90\%$, with steps of $20\%$. Results are reported in Figure 5 that shows the comparison between propeller manufacturer data (dashed lines) and the coefficients $K_S$ and $K_D$ obtained through simulation (continuous lines). The available data concern an existing cycloidal propeller, with the same geometry of the simulated one.

The discrepancies between real and simulated data appearing in Figure 5 have been supposed to be due to the stated simplifying assumptions about the interactions among the blades. In order to overcome this issue, some correction factors have been introduced, based on physical considerations: a shielding effect that can affect the oncoming flow for some blades and the interference of each blade with the others.

4.1. Shielding Correction

This correction concerns the shielding effect experienced by the blades that are in the half circumference not directly exposed to the oncoming water flow. Figure 6 gives a qualitative idea of how the flow is deviated by the blades. A corresponding correction factor, consisting in a matrix of corrective coefficients $w_s\left(s\right)<1$, has been introduced in order to reduce the velocity of the oncoming flow in the part of the rotor not directly invested by the flow itself. The corrective factors depend only on the driving pitch values and not on the advance coefficient $\mathit{\lambda }$ and are implemented as it follows

$$\widehat{u}=\left\{\begin{array}{cc}\widehat{u}{w}_s\left(s\right)\hfill & \mathrm{if}cos\theta <1\hfill \\ \widehat{u}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(41)

where $\widehat{u}$ is the velocity component of the rotor center along ${\underline{b}}_1$.

4.2. Interference Correction

As sketched in Figure 7, each blade influences the flow of the adjacent blade, so modifying the angle of attack of the incident flow itself. This interference among the blades has been modeled by reducing the attack angle of the incoming flow by a suitable quantity $w_i$. This correction depends on the advance coefficient $\mathit{\lambda }$ and the pitch values

$${\widehat{\alpha }}_{corr}=\widehat{\alpha }{w}_i(s,\mathit{\lambda })$$

(42)

where ${\widehat{\alpha }}_{corr}$ is the angle of attack defined in (15).

Figure 8 shows the values of free-running propeller coefficients, obtained by applying the above mentioned corrections. Although the proposed corrections are purely empirical, the graph in Figure 8 exhibits a good agreement between simulated and real data.

5. Simulation Results

In this section results concerning a simulated maneuver are presented. In this contest, two counter rotating CPs have been fitted on a vessel, whose dynamics were known. In the proposed simulation, initially the vessel is moving forward, the propeller pitch is set at $s=50\%$ and the shaft speed is constant. At a certain instant the steering pitch $\beta$ is required to move from 0 to ${20}^{\circ }$. Figure 9 shows the required and delivered shaft speed. Generally, this kind of propellers work at constant shaft speed and the thrust is controlled through the pitch value s, as for controllable pitch screw propellers. Figure 10 shows the ship motion, made dimensionless with respect to the vessel length. From top to bottom, advance motion, side drift and heading are respectively shown. Figure 11 shows the vessel speed, the components of the linear and angular vessel velocities in the body-fixed basis, and the vessel drift angle. As expected, at instant $100s$ when a twenty-degree change in steeering is required and kept constant for the rest of the maneuver, the ship begins to rotate and drift, it slows down, while drift and rotation velocities increase until they stabilize to constant values. The delivered forces and moment are reported in Figure 12, where they are expressed in the $\underline{b}-$basis. It is worth noticing that the model is able to evaluate the lateral forces (16) generated by each single propeller. Although two counter-rotating thrusters ensure the compensation of the lateral forces in the case of the forward navigation, this is no longer true during the maneuverings where the evaluation of such forces is an important aspect. Moreover, the implemented kinematic approach allows us to observe the asymmetry of $X_p$, that is the force component along ${\underline{b}}_1$, during the evolution. This is due to the fact that the two drifting thrusters are actually undergoing two different inlet flows. Such effect is also reflected in the torque. Indeed, time histories of required torque by each propeller at the shaft are reported in Figure 13. In this case, some differences can be evaluated during the vessel rotation where load on the portside shaft is different from the starboard one. The integration of the propeller models together with the propulsion system, Figure 14, allows us to assess the matching between available power, represented by the underlying area of the black curve, and the power required by the propellers at every time–step. Portside and starboard required powers (red and blue lines respectively) are clustered in a small area of the motor load diagram. Finally, the required fuel flow rate time history is reported in Figure 15, in terms of percentage of its maximum value, allowing us to compute fuel consumption.

6. Conclusions

This work presents a simulation approach for both low and high speed manoeuvring of ships equipped with cycloidal propellers. The real strength comes from not having to calculate the propeller fluid dynamics, avoiding demanding computations that would make very difficult and complex an effective simulation of the whole ship propulsion plant behaviour (in the proposed approach, CFD method is just used for the evaluation of the lift and drag coefficients of the single blade). However, reliable results, concerning both the steady state and transient performance of the cycloidal propeller, are achieved. This is essentially due to the rigorous description of the motion of each blade and by introducing specific empirical correction factors that can be used for a preliminary performance estimation of several cycloidal propulsion units, characterized by different lengths and number of blades. In this sense, the propulsion simulator can be regarded as a parametric one. Indeed, open water diagrams can obtained for a wide range of cycloidal propellers only by changing the rotor diameter, number and length of blades. Through appropriate insights concerning the correction coefficients and therefore a dedicated hydrodynamic analysis, the simulator can reproduce the open water performance maps of other types of cycloidal thrusters, for which the literature and the industry provide very little information. Right because the lack of data, the proposed approach could be useful also to predict the behaviour of the ship during the design phase in terms of general performances, forces generation, response times and evaluation of energy/fuel consumptions. Moreover, a training platform for personnel could be implemented on this basis. Next developments should include the integration of the hydrodynamic interaction between the hull and the cycloidal propeller and vice versa, as well as the calculation of the hydrodynamic resistance of the cycloidal propellers intended as hull appendages, mainly by comparing empirical corrections with CFD results. Further improvements could include the study of different sizes of the main engine, represented by a thermodynamic model coupled to the cycloidal thruster model, in order to better analyse the engine-propeller dynamics (especially during transient conditions). Finally, the present research should be completed through a proper validation of the simulation approach, by means of experimental data of a cycloidal propulsion system installed on board a real ship.