Energy-Saving Devices in Ship Propulsion: Effects of Nozzles Placed in Front of Propellers

Abstract

The hydrodynamic effect exerted by a nozzle placed in front of a KP505 propeller on the propulsive performances is studied by using extensive numerical simulations. The influence of a NACA 0015 nozzle with a chord length of 0.3 of the propeller diameter, D, mounted at 0.2 D in front of the propeller plane is studied for a various range of relevant nozzle diameters and different angles of attack. A detached eddy simulation (DES)-based hybrid technique implemented on the ISIS-CFD finite volume solver of the Numeca’s Fine^TM/Marine environment is proposed to fit the purpose. Systematically conducted simulations have proven that the net thrust reflecting the overall drag, which includes the nozzle, depends on the duct size. The duct presence determines two regions of the inflow into the propeller. One is the inner region of the nozzle where the high-speed flow exists because of the contraction of the duct. The other is the outer region of the nozzle where the flow decelerates due to the duct wake. Lower- and higher-pressure coefficients on the suction and pressure sides, cover a significantly wider area than those of the case without the nozzle, leading therefore to greater thrust and torque. The existence of a critical attack angle for which the magnitude of the relative axial force becomes maximum for the smallest nozzle diameter has been noticed.

1. Introduction

Greenhouse gas (GHG hereafter) emissions trap and hold heat in the atmosphere and contribute to climate change. Ships emitted approximately 1100 million tons of CO₂ in 2007, responsible for 3.5% of the global emissions. Based on the observed trend, the emissions in 2050 will be doubled compared to 2010. This scenario is in sharp contrast to the aim of reducing global GHG amount by at least 50% by 2050 and further to net zero in the second half of this century to mitigate the effects of climate change. To reduce the environmental impact of the maritime transport, a major improvement is required to reduce the energy consumption and increase the propulsive efficiency. This can be achieved by improving the ship design. An Energy Efficiency Design Index (EEDI) was imposed based on sea trials as the most important policy measure to reduce GHG emissions from shipping. Depending on the type and size, ships are not allowed to exceed a threshold for emitted CO₂ per ton-mile.

Sea trial procedures for EEDI mostly include actions to adjust the values in calm water conditions. The calm water focus is the traditional approach in the field and the research of hull shapes and propeller design had focused on optimization based only on calm water conditions. This procedure does not necessarily lead to a ship with the lowest energy consumption over its lifetime as it underestimates not only the non-linear effects due to waves, but also the local flow features that influence the propeller work. Global warming effects will contribute to more severe sea states in the future which will lead to higher risks of ship capsizing, loss of cargo and oil spill harming the marine ecosystem. Ship size increases constantly, so the added resistance in waves encounters new extremes as the ship length now approaches the length of large rogue waves. Under such circumstances, the propulsive performances need to be enhanced regardless of the hull geometry and the navigation conditions.

Obviously, aside from the main engine, the propeller represents the other key component that determines the efficiency of the fuel consumption for a ship. Since it also represents the main source of vibrations at board, several experimental [1,2] and theoretical studies [3,4,5] were carried out for establishing a reliable and stable method for its optimal design. In terms of the theoretical studies, both potential and viscous-based numerical approaches have been proposed over time to solve the problem for the calm water and in real sea navigation cases [6,7,8,9,10,11].

Both wetted and cavitating propellers were considered to fit the purpose and some encouraging results have been reported so far [12,13,14]. Moreover, various experimental and numerical studies were continuously attempting at defining atypical propellers such as tip rake propellers [15,16,17] or contracted load tip propellers [18,19] meant to contribute to the improvement of the overall propulsive performances.

Lowering the vibration level can be achieved either by improving the propeller design or by optimizing the wake field behind the hull. Since the first method requires validation in the cavitation tunnel which is rather expensive and time-consuming, the attention has focused towards improving the ship wake. A uniform wake proved to efficiently contribute to a reduction of the risk of cavitation and to lowering the level of propeller-transmitted excitations. Besides, a uniform wake may determine a significant increase of the propeller efficiency.

Intensive studies have proven that a uniform wake can be obtained by placing several flow control devices in front of the propeller. Although they were initially devised for improving the propulsive efficiency, finally it has been proven that the cavitation behavior, vibration levels and maneuverability performances could also be improved based on an appropriate design process. Notwithstanding that control devices such as wake equalizing nozzles, Grothues spoilers, stator fins, or different combinations of them were studied, no clear conclusions on their efficiency could be withdrawn so far. Wake equalizing nozzles previously studied in References [20,21,22,23] may accelerate the inflow into the upper region of the propeller where the axial velocity is small relative to the corresponding lower region, and in doing so, it improves the uniformity of the wake over the entire propeller disk. Moreover, a well-designed equalizer proved to contribute to a significant decrease of the amount of flow separation at the stern. In certain working conditions, the solution may generate additional thrust, reduce propeller-induced vibrations and improve the overall steering characteristics because of the increased incoming axial speed.

Wake equalizers were firstly proposed by Schneekluth [20], where it was described how the wake equalizing duct accelerates the inflow of the propeller in the area of powerful wake, thus reducing the required propulsive power. Schneekluth detailed the mechanism of accomplishing this task by improving propeller efficiency, generating a lift with a forward force component on the foil section, thus reducing the flow separation at the aft body. The author also proved that the propeller-excited vibrations can by reduced by decreasing the propeller tip loading. A study of the energy-saving in the powering characteristics was carried out in Reference [21]. Two different hull forms were generated from the original hull form of the vessel to optimize the stern flow of the vessel. A partial wake equalizing duct was investigated by measuring the resistance and self-propulsion characteristics and by visualization the flow features. The author concluded that the partial wake equalizing duct concept with an appropriate stern design affects not only the flow characteristics at the aft end, but also the propulsion characteristics.

Çelik [22] computed the effect on the propulsive performance of a chemical tanker by a wake equalizer duct (WED) and concluded that propeller characteristics and the ship resistance are slightly affected by the WED, but an additional thrust could be obtained. The author reported that the maximum gain obtained by using an appropriate WED design is about 10%. A comprehensive work on the subject was performed by Go et al. [23], who numerically investigated the effect of a nozzle placed in front of a propeller on the propulsion performances. Numerical simulations were performed for a wide range of diameters and angles of attack and it was eventually shown that a duct with a high angle of attack provides faster inflow into the propeller than in the case without a duct, in which it may improve the thrust of the propeller, but simultaneously burdens the torque.

Based on the above literature survey, it may be concluded that the performance exerted by a nozzle placed in front of a propeller as an energy-saving device should naturally be evaluated by taking into consideration the ship hull forms. However, given the multitude of the parameters that usually contribute to the performances of a tandem energy-saving device (ESD) and a propeller (nozzle diameter, the angle of attack, the relative distance between the two and so on), a series of systematic studies are still required to grasp a better understanding on how this ensemble may work efficiently from the propulsive point of view.

The present paper describes a set of numerical simulations of the three-dimensional (3D) viscous flow around a NACA0015 nozzle placed in front of the five-bladed KP505 propeller. Four different diameters and five attack angles were considered for studying the cumulated influence of the duct on the hydrodynamic performances of the propeller. For comparison purposes, a set of computations is initially performed for the propeller working alone in open water at various advance coefficients ranging from 0.1 to 1. The finite-volume ISIS-CFD solver of the Fine^TM/Marine suite is employed in all the computations reported here. Closure to turbulence is achieved by using a hybrid DES-shear stress transport (SST) model. The rotating frame technique is used for the single propeller computations, whereas sliding grids are used when the nozzle is also considered. Numerical uncertainty analyses were preceding the solutions in both cases to reveal the level of the simulation errors. In the present study, a uniform inflow condition is imposed without considering a certain ship to investigate the net effect of a duct on the propulsive performance. Notwithstanding that the inflow condition induced by a hull is crucial to choose the ESD geometry, the present work aims at providing a methodology for investigation of the performances of a nozzle-propeller ensemble regardless of the ship type. Based on the mutual interaction between the incoming flow particulars determined by the particulars of the nozzle and propeller, a better physical insight into the flow is expected to help in designing an optimal ESD.

2. Computational Milestones

2.1. Propeller and Nozzle Geometries, Computational Domain and Mesh Particularities

The right-handed five-blade KP505 propeller model, whose main particulars are tabulated in

Table 1. Main particulars of the KP505 propeller.

Main particulars	Value
Diameter, $D$(m)	0.250
Pitch/Diameter mean, $P/D$	0.950
Expanded area ratio, $Ae/Ao$	0.800
Hub Ratio	0.180
Rake at tip (m)	0
Effective skew (degrees)	32
Direction of rotation	Right-hand
Section	NACA66

, is considered in the present study. The propeller shown in Figure 1a, designed by Korea Research Institute of Ships and Ocean Engineering (KRISO) for the KCS container ship, is working in a pulling condition, the same as in the open water experimental tests. The flow around the 250 mm diameter is simulated for a set of 19 incoming flow velocities so that the resulted values of the advance coefficient will be between 0.1 and 1.

The computational domain consists of a cylindrical volume whose diameter is 5 times the propeller diameter, and length is 12 times the propeller diameter, out of which 4 diameters are for the upstream of the propeller. The open water single propeller was computed by using the rotating frame technique, while when the nozzle was considered, sliding grids were used. The nozzle depicted in Figure 1b has a NACA0015 profile. The chord length of the duct is $0.3D$ and its trailing edge is located $0.2D$ upstream from the center of the propeller. A set of five angles of attack of the duct,$\theta$, ranging from $0°$ to $20°$, increased in equal steps of $5°$, is considered in the numerical investigation. The nozzle diameter denoted by $D_{N }$ in Figure 1b was successively varied from $0.7D$ to $D$ in steps of $0.1D$.

Several grids were generated to host the discretization of the differential equations of the flow. The unstructured hexahedral HEXPRESS grid generator of the Numeca’s Fine^TM/Marine package was employed to generate the computational meshes. The generator has a CAD direct import capability which allows the modification or the decomposition of the geometry. This facility makes possible the split of the solid bodies in subsequent sub-surfaces regardless of their geometric complexity. The fully automatic hexahedral grid generation was accomplished as the buffer cell and boundary layer insertion for high-quality cells inside. The dimensionless wall distance $y^{+}$ on the propeller blades, hub and cap, as well as on the nozzle walls, was imposed at a subunit value, such that the viscous sublayer be properly resolved. Areas of clustering cells were placed around the blade and nozzle edges, where the flow gradients are expected to be high, as shown in Figure 1c,d. Since the development of physical parameters in the propeller wake is important for a comprehensive understanding of the whole flow mechanism, annular volumes of cell clustering were placed in the wake of the propeller tips and hub, as depicted in Figure 1e.

Figure 1f shows the two sliding grids used for the flow computation around the nozzle-propeller configuration. A first grid is surrounding the shaft, whereas the second is for the propeller. The sliding grid method works in parallel and the different sub-domains can be distributed arbitrarily over the processors. No explicit interpolation was used to find the states on the sliding faces; instead, the coupling algorithm identified real cells that were used as neighbors for the cells in the other sub-domain. The fluxes were computed across the sliding faces independently of the appertaining sub-domain and therefore without special concern related to conservation of mass and momentum. Moreover, the cells had to have the same level of refinement on both sizes of the sliding domains, such that the transmission of information between the components be carried out with the least possible errors. To achieve that purpose, the minimum cell size has been chosen at the same value on the interface of the sliding grids, as Figure 1e shows.

2.2. Numerical Highlights

The 3D ISIS-CFD unsteady viscous flow solver was employed in the present study. The simulation was accomplished in a global approach in which the momentum and mass conservation equations written with respect to a Cartesian system of coordinates are solved. The solver is based on the finite-volume method to build the spatial discretization of the transport equations on unstructured grids. The face-based method is generalized to unstructured meshes composed of arbitrary volume shapes. Fluxes are built using the AVLSMART scheme, which is based on the third-order QUICK scheme. Velocity field is obtained from the momentum conservation equation and the pressure field is extracted from the mass conservation constraint transformed into a pressure equation. Picard’s procedure is used for the linearization of the equations. The whole discretization is fully implicit in space and time and is formally second-order accurate. Closure to turbulence was achieved through a hybrid model, which used the shear stress transport (SST)-based DES. The model provides the accuracy of large eddy simulation (LES) for highly separated flow regions and computational efficiency of the Reynolds Averaged Navier-Stokes (RANS) equations in the near-wall region.

The flow was accelerated from rest up to a rotational speed of 600 revolutions per minute (rpm) over a certain period of time based on a semi-sinusoidal ramp in all the cases reported here. The axial velocity was kept unchanged at a value that corresponded to the advance coefficient the computation was performed at. The time step was chosen depending on the minimum grid size, such that the resulting Courant number be less than unit in all the computations, as it will be described in the following. All the computations were carried out in two successive runs. In the first run, closure to turbulence model was achieved by the classic $K-\omega$ SST model of Menter [24]. Computations were performed for about 20 seconds, until the thrust stabilized completely. The output of the first run was then used as input for the second run, which employed the hybrid DES-SST turbulence model. The time step of the second run was decreased four to five times, so that the Courant number remained below 0.2.

All the solutions discussed in the present paper were computed in parallel on 240 processors in double precision, by splitting the computational domain into multiple connected sub-domains having approximately the same number of unknowns. The METIS partitioning algorithm was used to accomplish the division, and the communication of faces’ data between sub-domains was performed according to the Message Passing Interface (MPI) standard. All the computations employed the DES turbulence model, and all the reported solutions were computed for 25 seconds physical time when the propeller had already performed 250 complete rotations. For the finest grid case, a number of about 251,500 cells were computed on each processor. All the computations were performed on a high-performance computer with 624 Intel E5 2680 v3 processors with 12 cores, with a frequency of 2.5 and 3.3 GHz for the turbo boost regime. In the finest mesh computational case, an average of 0.171 milliseconds per mesh point was necessary.

2.3. Boundary Conditions

The boundary conditions for the simulations reported here are either Dirichlet or Neumann type. That is, a Dirichlet condition is imposed for the incoming velocity at the upstream, at the downstream and on the exterior boundaries, whereas a Neumann condition for pressure $\left(\partial p/\partial n=0\right)$ is imposed at the downstream. The no-slip condition is used on all the solid boundaries. In ISIS-CFD, the choice of default free-stream values is:${\omega }_{\infty }=\lambda U_{\infty }/D, {\mu }_{t\infty }={10}^{-3}{\mu }_t$ and $k_{\infty }={\mu }_{t\infty }{\omega }_{\infty }/\rho$, where $U_{\infty }$ is the oncoming axial velocity at the far inlet. Here, $\lambda$ is a factor of proportionality, which is recommended to be between 1 and 10. The limits are comparable with values recommended for boundary layer flows, therefore $\lambda =10$ was chosen [24]. It should be noticed that free-shear layers are more sensitive to small free-stream values of ${\omega }_{\infty }$ and larger values of $\omega$ are recommended in the free-stream: at least $\lambda =40$ for mixing layers, increasing up to $\lambda =80$ for round jets. Wall boundary conditions are taken as $k=0$ and $\omega =60\mu /\beta \rho {\left(\mathsf{\Delta }y\right)}^2$, where $\Delta y$ is the distance of the first point away from the wall and $y^{+}<1$, thus allowing to directly resolve the viscous sublayer. This condition led to a maximum cell dimension of 10^–6 meters measured normal to the solid walls. The growth ratio of the cells through the inflation layers is 1.15 for all solid surfaces of the propeller in all computational cases.

3. Results and Discussions

3.1. Propeller Working in Open Water (POW)

For the given computational domain size, four different grids were generated for the POW grid convergence test. Let them be denoted by G1–G4. G1 consists of 7.12 × 10⁶ cells, G2 of 15.01 × 10⁶ cells, whereas G3 and G4 consist of 31.32 × 10⁶ and 60.3 × 10⁶ cells, respectively. The test is performed for a set of three advance coefficients, namely J = 0.2, J = 0.7 and J = 1.0. The grid convergence test was carried out in terms of the computed thrust and torque coefficients compared with the experimental data reported in References [25,26], as

Table 2. Grid sensitivity test. Comparisons between the thrust and torque coefficients measured (EFD) and computed (CFD) on the four meshes (G1–G4).

		${\mathit{K}}_{\mathit{T}}$			$10{\mathit{K}}_{\mathit{Q}}$
J		0.2	0.7	1	0.2	0.7	1
EFD		0.435	0.185	0.022	0.622	0.311	0.096
CFD	G1	0.4489	0.1909	0.0227	0.6410	0.3189	0.0981
	$\left|\epsilon \right| \left(\%\right)$	3.096	3.091	3.084	2.964	2.477	2.141
	G2	0.4471	0.1901	0.0226	0.6360	0.3174	0.0979
	$\left|\epsilon \right| \left(\%\right)$	2.706	2.683	2.655	2.201	2.016	1.941
	G3	0.4245	0.1894	0.0225	0.6344	0.3171	0.0975
	$\left|\epsilon \right| \left(\%\right)$	2.473	2.323	2.222	1.955	1.924	1.538
	G4	0.4292	0.1868	0.0216	0.6298	0.3064	0.0948
	$\left|\epsilon \right| \left(\%\right)$	1.351	0.964	1.852	1.238	1.501	1.266

shows. The figures tabulated in Table 2 show a monotonic decrease of the computational error for both thrust and torque coefficients with the increase of the resolution of the discretization. The same conclusion may be drawn from Figure 2, which depicts the evolution of the relative error with the d-th root order of the number of cells, Nc, of the four grids considered. Here, d stands for the dimensionality of the flow problem.

The verification and validation (V&V) were performed only for the propeller working in open water since there is no experimental data available for the nozzle-propeller tandem. Any consistent study should use a minimum of three solutions to evaluate the convergence with respect to input parameter, so four unstructured grids with a similar topology were generated, such that the dimensionless distance to the wall of the nearest cell center remained $y_C^{+}<1$. Grid dependency tests were conducted using the four grids described above, which enabled two separate grid studies to be performed and compared. The first study gives estimates for grid errors and uncertainties on G4 using the three finest grids, G2–G4, whereas the second gives estimates for grid errors and uncertainties on G3 using the three coarsest grids, G1–G3.

The grid refinement ratio $r_{Gi}={\left(N_C\left(G_4\right)/N_C\left(G_i\right)\right)}^{1/3}$, its number of cells $N_C$, and ${\left(y_C^{+}\right)}_{\max }$ computed according to References [27,28] are tabulated in

Table 3. Properties of the grid sets used in the simulations.

Grid	${\mathit{N}}_{\mathit{C}}$	${\mathit{r}}_{\mathit{G}\mathit{i}}$	$1/{\mathit{r}}_{\mathit{G}\mathit{i}}$	${\left({\mathit{y}}_{\mathit{C}}^{+}\right)}_{\mathbf{max}}$
G1	7,121,306	0.491	1.00	0.491
G2	15,011,242	0.629	1.282	0.364
G3	31,320,461	0.804	1.638	0.213
G4	60,305,214	1.00	2.038	0.182

. Modeling errors were defined based on the general rule which states that for an arbitrary flow parameter φ, the interval which contains the modeling error $\delta$ for 95 out of 100 cases is given by $E_{\phi } - U_{val} \le {\delta }_{\phi } \le E_{\phi } + U_{val}$, where $U_{val}$ is the uncertainty value and $E_{\phi }$ is the comparison error defined as $E_{\phi }=S_{\phi }-D_{\phi }$. Here, $S_{\phi }$ and $D_{\phi }$ are the values of the numerical solution and the experimental measured data, respectively. Since neither $S_{\phi }$ nor $D_{\phi }$ are free of errors, numerical $U_S$ and experimental $U_D$ uncertainties could therefore be established.

Effectively, the grid uncertainty was evaluated by using the methodology described in Reference [29] for monotonic convergence. In the following, the verification and validation will be described for the thrust coefficient only. The three advance ratio values tabulated in Table 2 will be considered here as well. The finer grids G4 and G3, which were built by halving the initial cell sizes of G2 and G1 on the boundaries, will be considered first. The grid ratio ($r_G$), the associated relative error between the $K_T$ computed on the finest mesh G4 and the second finest mesh G3, ${\epsilon }_{43}\%C_{T4}$, the ratio between the estimated order of convergence and the theoretical order of convergence, $p_G/p_{G, th},$ the grid uncertainty, $U_G\%G_4$, the experimental uncertainty, $U_D\%D,$ and the validation uncertainty, $U_V\% ,$ are tabulated in

Table 4. Verification for the thrust coefficient for $J=0.2, 0.7, 1$, grids G4, G3.

J	${\mathit{K}}_{\mathit{T}, \mathit{G}4}$	${\mathit{K}}_{\mathit{T}, \mathit{G}3}$	${\mathit{r}}_{\mathit{G}}$	${\mathit{\epsilon }}_{43}\mathbf{\%}{\mathit{K}}_{\mathit{T}}$	${\mathit{p}}_{\mathit{G}}/{\mathit{p}}_{\mathit{G}, \mathit{t}\mathit{h}}$	${\mathit{U}}_{\mathit{G}}\mathbf{\%}{\mathit{G}}_4$	${\mathit{U}}_{\mathit{D}}\mathbf{\%}\mathit{D}$	${\mathit{U}}_{\mathit{V}}\mathbf{\%}$
0.2	0.4292	0.4245	1.925	1.122	1.347	1.788	1.0	2.049
0.7	0.1868	0.1894	1.925	1.359	1.211	1.635	1.0	1.917
1.0	0.0216	0.0225	1.925	0.370	1.081	1.551	1.0	1.845

. The $p_{G, th}$ in Table 4 is the theoretical order of accuracy, which is given by the order of convection scheme, whereas the validation uncertainty is expressed as $U_V=\sqrt{{\left(U_G\%G_4\right)}^2+{\left(U_D\%D\right)}^2}$. The relative error between the finest mesh and the experimental data for the coefficient of total resistance on the finest grid G4 is smaller than the Richardson-based numerical uncertainty (see Table 2 and Table 4), and therefore the $K_T$ prediction may be considered as being validated.

Similarly, for the second verification, the results are tabulated in

Table 5. Verification for the thrust coefficient for $J=0.2, 0.7, 1$, grids G2, G1.

$\mathit{J}$	${\mathit{K}}_{\mathit{T}, \mathit{G}2}$	${\mathit{K}}_{\mathit{T}, \mathit{G}1}$	${\mathit{r}}_{\mathit{G}}$	${\mathit{\epsilon }}_{21}\mathbf{\%}{\mathit{K}}_{\mathit{T}}$	${\mathit{p}}_{\mathit{G}}/{\mathit{p}}_{\mathit{G}, \mathit{t}\mathit{h}}$	${\mathit{U}}_{\mathit{G}}\mathbf{\%}{\mathit{G}}_4$	${\mathit{U}}_{\mathit{D}}\mathbf{\%}\mathit{D}$	${\mathit{U}}_{\mathit{V}}\mathbf{\%}$
0.2	0.4471	0.4489	2.108	0.390	1.453	1.917	1.0	2.162
0.7	0.1901	0.1909	2.108	0.408	1.302	1.781	1.0	2.043
1.0	0.0226	0.0227	2.108	0.429	1.147	1.633	1.0	1.915

.

Based on the V&V and grid sensitivity test findings above, grid G4 is further chosen for the following POW simulation. The G4 is shown in Figure 1c. Special effort has been made to ensure the fulfillment of the requirements concerning the orthogonality, smoothness and clustering around the blade edges, the tip and the blade intersection with the hub. The resulting number of cells for each blade face is approximately $7.5\times {10}^5$. A viscous layer on 48 cells, which is corresponding to covering the boundary layer thickness, is introduced to keep the y⁺ below unit.

POW computations were performed based on a rotating frame technique for a series of 19 advance coefficients, i.e.,$0.1\le J\le 1$. As mentioned above, all the simulations were performed at the same rotational velocity of 10 rps, so the different advance coefficients were obtained by modifying the incoming axial velocity. A comparison between the experimental data provided in Reference [26] and the corresponding numerical solution is tabulated in

Table 6. Comparison between the computed and measured open water coefficients.

$\mathit{J}$	${\mathit{K}}_{\mathit{T}}$			$10{\mathit{K}}_{\mathit{Q}}$			${\mathit{\eta }}_0$
	CFD	EFD	$\left|{\epsilon }_{K_T}\right| \left(\%\right)$	CFD	EFD	$\left|{\epsilon }_{K_Q}\right| \left(\%\right)$	CFD	EFD	$\left|{\epsilon }_{{\eta }_0}\right| \left(\%\right)$
0.10	0.4765	0.482	1.154	0.6706	0.677	0.954	0.1131	0.113	0.088
0.15	0.4451	0.458	2.898	0.6385	0.646	1.175	0.1665	0.169	1.501
0.20	0.4292	0.435	1.351	0.6298	0.622	1.238	0.2170	0.223	2.765
0.25	0.4152	0.412	0.771	0.6091	0.589	3.300	0.2713	0.278	2.470
0.30	0.3788	0.387	2.165	0.5416	0.557	2.843	0.3342	0.332	0.658
0.35	0.3710	0.361	2.695	0.5214	0.531	1.841	0.3965	0.379	4.414
0.40	0.3480	0.336	3.448	0.5003	0.497	0.660	0.4431	0.431	2.731
0.45	0.3009	0.310	3.024	0.4734	0.466	1.563	0.4554	0.477	4.743
0.50	0.2834	0.285	0.565	0.4464	0.437	2.106	0.5055	0.519	2.671
0.55	0.2581	0.259	0.349	0.4122	0.405	1.747	0.5484	0.561	2.298
0.60	0.2379	0.235	1.219	0.3822	0.376	1.622	0.5946	0.597	0.404
0.65	0.2109	0.209	0.901	0.3399	0.343	0.912	0.6423	0.631	1.759
0.70	0.1868	0.185	0.964	0.3074	0.311	1.171	0.6774	0.665	1.831
0.75	0.1573	0.161	2.352	0.2733	0.278	1.720	0.6872	0.691	0.553
0.80	0.1332	0.137	2.853	0.2418	0.247	2.151	0.7019	0.705	0.442
0.85	0.1145	0.111	3.057	0.2106	0.215	2.089	0.7163	0.697	2.694
0.90	0.0860	0.083	3.488	0.1783	0.181	1.514	0.6812	0.654	3.993
0.95	0.0522	0.053	1.533	0.1380	0.140	1.449	0.5719	0.575	0.542
1.00	0.0216	0.022	1.852	0.0948	0.096	1.266	0.3622	0.373	2.982

and depicted in Figure 3, where the experimental data are drawn with lines and the computed solutions by symbols. The errors are shown by bars whose lengths are drawn proportionally to $\left|\epsilon \right|$, which is computed in Table 6.

Whenever the open propeller simulation features are analyzed, the local flow characteristics in various cross-planes around the acting propeller should be considered. This is also needed as a basis for comparisons with the computational cases with nozzle when the performances of the ensemble will be discussed. In this respect, Figure 4 proposes a comparison of the instantaneous non-dimensional axial velocity fields computed in three different cross-planes placed at $x=-0.5R$, $x=0$ and $x=0.5R$ respectively, for the three advance ratios of $J=0.5$, $J=0.7$ and $J=1.0$. It is worth mentioning that all numerical solutions are reported to a Cartesian system of coordinates whose origin is placed in the center of the propeller. Obviously, the comparison reflects the increase of the streamwise velocity that corresponds to the higher advance ratio coefficients compared to the lower J values owing to the constant rotational speed. Nonetheless, discrete regions of lower axial velocities, which correspond to the vortices’ cores, are seen around the tips of the blades in the propeller plane and behind it. Figure 4 also shows that higher gradients of the axial velocity are detected in the heavier loading conditions, as expected.

Obviously, the same conclusion may be drawn from the instantaneous non-dimensional streamwise velocity contours drawn in the longitudinal vertical plane of the propeller depicted in Figure 5. Tip vortex locations are clearly marked by lower axial velocity, whose values are lower in the near-wake of the blades, as expected. The gradients of the axial velocity in the cores of vortices get smaller downstream in the wake because of the viscous dissipation. Apart from that, a rope-like swirl vortex is released by the hub and propagates in the downstream of the propeller. Its intensity is decreasing with the increase of the advance ratio, thus confirming the theory. The periodicity of the tip vortices is almost constant, a fact which suggests that the numerical dissipation induced by the solver is small enough.

Directly related to the instantaneous velocity is the instantaneous vorticity shown in Figure 6, since the latter is defined as twice that of the rotation vector, i.e., the curl of the velocity vector. Strong tip-released structures are shed in the wake, correlating with local maxima of turbulent kinetic energy. The good agreement between the velocity and vorticity is explainable since they derive one from the other. This can be seen in Figure 6, which depicts the instantaneous vorticity contours computed for the same advance coefficients as in the non-dimensional axial velocity representation. The figure shows that regardless of the J value, the near-field is dominated by well-defined coherent tip and root vortices. They manifest as spirals of a regular shape released by each blade. The author believes that this is the merit of the proper mesh clustering of the wake in those regions, as shown in Figure 1e.

On the contrary, owing to a rather poor meshing of the area in between the two, the less intensive blade trailing edge wake looks rather underestimated. The tip and root vortical structures, well-defined inside the near-wake, become unstable and eventually vanish to form the far wake. The mechanism of their formation was previously described by the author in [14] for the tip vortices and in [30] for those released by the intersection between the blades and hub, therefore no more assertions will be provided here. The vortex cores’ positions correspond strictly to the restricted regions of lowest pressure values in the wake domain, as well as to the velocities and turbulent kinetic energy fields. Tip vortices breakdown at a certain distance from the propeller. The location of transition to instability strongly depends on spiral-to-spiral distance as the number of blades affects the onset. Although independent from each other, the perturbation resulting from the destabilization of tip vortices can cause the hub vortex to become unstable as well.

In the following, a comparative analysis of the topology of the vortices released by the open water KP505 propeller is proposed in Figure 7, where six different computed advance ratios varying from 0.1 to 1.0 are considered. The geometries of the vortices were defined based on the $Q$ criterion.

Vortical structures shown in Figure 7 are drawn as iso-surfaces of $Q=15,000$ and colored by helicity. The propeller wake consists of five pairs of helical tip and root vortices (one originating from each blade). In all the depicted cases, the tip vortical structures are observed until the end of the refinement grid, where the mesh resolution becomes too coarse and the vortices are lost due to diffusion, regardless of the advance ratio. Obviously, this is a merit of the hybrid DES-SST turbulence model employed in the present study. Following a fairly small contraction of the radial location, caused by the acceleration of the flow behind the propeller, the helices formed by the tip vortices remain located on a cylinder of almost constant radius up to the end of the grid refinement block. Vortices that form at the root of the blades that eventually merge into the hub vortex are also visible for all advance ratio computations.

Figure 8 shows the non-dimensional dynamic pressure distribution on both faces of the propeller computed for $J=0.9$. The neglect of the static pressure explains the relative similarity of the distribution on the blades shown in the figure. High dynamic pressure areas that are located around the tips and behind the leading edges on the pressure side are found. Strong negative gradients are seen around each relative radius. On the contrary, areas of significant negative pressure are detected on the suction side, which are originating at about 0.8R and extend downwards to the blade root, almost following the reference line. The lowest pressure values correspond to the root areas, regardless of the blade position.

3.2. Propeller-Duct Ensemble Working in Open Water

In the following, the open water flow is studied by considering the NACA0015 nozzle placed in front of the propeller, as depicted in Figure 1b. The chord length of the duct is 0.3D and its trailing edge is located $0.2D$ upstream from the center of the propeller. A set of five angles of attack of the duct, $\theta$, ranging from 0° to 20°, increased in equal steps of 5°, is considered in the numerical investigation. The different attack angles were obtained by rotating the hydrodynamic profile around its trailing edge. Four nozzles’ equally increased diameters were also considered. They ranged from $D_N=0.7D$ to $D_N=1.0D$. The simulation was performed for an advance coefficient of 0.7. For consistency reasons, the flow was advancing with an axial velocity of 1.75 m/s, exactly as in the self-propeller case. The duct existence induces significant changes not only in the velocity field in the propeller plane, but also in the pressure fields on the blades, as it will be shown later when the local flow will be discussed. Those changes consequently determine modifications of the propeller thrust, torque and efficiency, as depicted in Figure 9.

The figure shows the variation of the three parameters with the diameter and the inclination angle of the duct. For comparison reasons, the computed and measured values of the force coefficients that correspond to the no-duct propeller case are also drawn. Regardless of the nozzle diameter,$D_N$, the values of $K_T$, $K_Q$ and ${\eta }_0$ show the same tendency of increasing patterns with respect to the attack angle $\theta$, except for the efficiency that corresponds to the $\theta =0°$ and $\theta =5°$ cases, in which losses of less than 2% were observed. Force and moment coefficients corresponding to the propeller-nozzle case are larger than the corresponding values for the propeller without the duct, except the efficiency computed for the above-mentioned two attack angle cases. The present findings resemble those of Go et al. [23], who performed a similar investigation, but for a higher axial velocity.

The most significant increase of the thrust and torque coefficients corresponds to the smaller duct diameters, except for the efficiency computed for the two lower values of the attack angle, as mentioned above. In the considered range of $D_N$, the differences in magnitudes of corresponding coefficients are rather small for a fixed value of $\theta$. However, as $D_N$ increases, the magnitude of the coefficient increase becomes smaller, as Figure 10 shows.

Figure 11 depicts a comparison of the pressure distributions on the nozzle and the propeller suction and pressure sides, respectively. The flow solution is computed in this case for nozzles with $D_N = 0.9D$ inclined at $\theta = 0°$, $10°$ and $20°$, respectively. Comparing the solutions drawn in Figure 11 to that shown in Figure 8 for the single propeller case, one may notice the wider area of lower pressure on the suction side than in the case without the duct. Besides, the non-dimensional positive pressure covers a much wider area on the pressure side, but with a rather different distribution on the blade surface. Obviously, the lower and higher pressures on the suction and pressure sides are ultimately responsible for the larger thrust and torque than in the case without a nozzle. Another important issue is related to the tip areas of the pressure sides of the blades, where the pressure does not reach the maximum value, which could be clearly seen in Figure 8. That is seemingly due to the nozzle, which is exerting a significant influence on the flow in the propeller plane. A sharp gradient of pressure is found around the leading edge between the relative radius of 0.4 and 0.9 as a result of the sharp velocity gradient related to the interface of the axial velocity layers.

The suction side pressure distribution is somehow consistent with the no-duct propeller computation. Nevertheless, in the $D_N = 0.9D$ at $\theta = 10°$ and $\theta = 20°$ cases, wider areas involving lower- and higher-pressure coefficients on corresponding surfaces compared to the case without the duct are present, as shown in Figure 11. Numerical solutions not reported here have shown that these augmented areas for the case of $D_N = 0.9D$ are smaller than in the case $D_N=0.7D$, which somehow confirms the difference in force coefficients between $D_N = 0.9D$ and $D_N = 0.7D$ at higher attack angles.

Figure 12 proposes a comparison of the instantaneous non-dimensional axial velocity fields computed for $D_N/D=0.7, J=0.7$ at 0.5R upstream of the propeller (left), in the propeller plane (middle) and at 0.5R downstream of the propeller (right) for $\theta =0°$, $\theta =10°$ and $\theta =20°$, respectively. The suction effect induced by the propeller affects the upstream flow, thus regions of axial velocity gradients distribute almost uniformly along the circumferential direction, rotated at angles that correspond to the positions of the blades. The acceleration induced by the duct contributes to an increase of the axial velocity in the slipstream compared to the case without the duct shown in Figure 4b. As a consequence, the theoretical border between the slipstream and freestream becomes more perceptible in the wake of the propeller, regardless of the attack angle of the nozzle. Velocity contours formed by the root-side vortices behind the propeller, which were rather prominent in the case without the duct, are slightly swept by the duct-induced faster flow and are loose in intensity. As $\theta$ increases, significant gradients are found inside both inner and outer corresponding areas of the nozzle, as comparatively shown in Figure 12b,c.

Next, a flow-field analysis is performed by means of the velocity and pressure distributions in the wake of the propeller. Figure 13 shows typical instantaneous velocity fields in the vertical plane of $y=0$ for the single propeller as well as for the duct-propeller arrangement drawn for different angles of attack in the computational case of $D_N=0.9D$. The x-axis and z-axis are normalized with respect to the propeller radius, R, and lie parallel to the main flow and vertically upward, respectively. In spite of the rather poor mesh resolution in the mid-area of the propeller, both freestream and slipstream are clearly distinguished in the single-propeller case shown in Figure 13a. Previous studies [23,31], in which the turbulence was modeled by using the $k-\omega$ SST model, reported failures not only in capturing the tip vortices, but also in establishing a clear difference between the freestream and slipstream in the wake region.

It is worth mentioning that in both referenced studies, wall functions were used as boundary conditions on the solid surfaces. On the contrary, in the present study, the DES-SST model based on an implicit splitting of the computational domain into two separate zones was employed. In the first region near the solid walls, the conventional RANS equations were solved. Inside the second region, the filtered Navier-Stokes equations of the LES approach were solved. Since the hybrid DES-SST model only works with no-slip conditions for velocity on the solid boundaries, therefore, without employing any wall functions, it appears that the merit in properly resolving the vortical structures resides in the more careful treatment of the boundaries.

When the propeller-nozzle arrangement is considered, several differences in the flow fields may be emphasized, as Figure 13b–e shows. With the increase of the angle of attack, significant modifications of the freestream become visible. Starting with $\theta =10°$, the tip vortices appear to move slightly over the outer surface of the duct, as depicted in Figure 13d,e. For further increases of $\theta$ from 10° to 20°, the vortices become stronger and rather unstable, seemingly due to the interference with the stronger vortex released by the trailing edge of the duct.

This behavior is even more prominent when the duct diameter is smaller, as shown in Figure 14, which shows the instantaneous velocity vector fields in the propeller wake computed at an advance coefficient of $J=0.7$ for $D_N/D=0.7$. As described above, the decrease of the nozzle diameter induces modifications of the streamwise velocity distribution not only in the propeller plane, but also in its wake.

Generally speaking, a duct with a smaller diameter and a high angle of attack provides faster inflow into the propeller than the case without the duct, which improves the thrust of the propeller, as shown in Figure 9. As mentioned above, there are two regions of the inflow into the propeller plane, which are determined by the nozzle. The first one corresponds to the inner region of the duct where a high-speed flow rate is detected because of the contraction of the duct. Within the inner region, the root vortices become unstable and eventually pair. The other region is placed outside of the nozzle, where the flow decelerates due to the nozzle wake. Obviously, inside the outer region, lower speed flow induced by the duct wake leads to a lower advance ratio, which ultimately increases the thrust. The combined effects of the inner and outer regions are more prominent with the increase of the angle of attack.

The above statement is sustained by the contours of pressure normalized with the atmospheric pressure computed at $J=0.7$ for the single propeller and for the propeller-nozzle arrangement in the $y=0$ plane at $D_N/D=0.7,$ and depicted in Figure 15. The comparison reveals that the extension of the higher pressure around the propeller tip increases from the single propeller case to the propeller-duct flow case as the nozzle angle of attack increases. Lower-pressure well-defined areas that correspond to the core of the vortices can be easily distinguished. The pressure gradients inside those areas get stronger as the angle of attack of the nozzle increases. Other negative pressure areas, directly dependent on the increase of the angle of attack, are seen in the wake of the hub, a fact that may explain the favorable condition for vortex pairing.

A comparison of the vorticity fields drawn in the longitudinal plane of $y=0$ is further proposed in Figure 16. The vorticity contours computed for the $D_N/D=0.7$ case are drawn for all five attack angles considered in the present study. Obviously, the nozzle placed upstream of the propeller also determines significant changes in the vorticity developed in the wake, as may be seen if the plots are compared with the single propeller case shown in Figure 6b. The near-wake of the duct leads to a modification of the tip vortices’ distribution and intensity. As the angle of attack increases, the freestream change becomes considerable because of the stronger influence of the vorticity released by the nozzle. Both root and tip vortices become unstable, so that for larger attack angles, the vortex pairing becomes imminent, as Figure 16c–e proves. It is also worth mentioning that the attack angle increase determines a stronger swirl vortex.

To get more insight into the impact of the duct on the propeller behind, a local flow analysis is proposed in the following by checking the nozzle influence on the wake structure, based on a close-up analysis of the interaction between the two. As explained above, when flow around a single propeller was considered, one could find that regardless of the advance ratio, the trajectory of tip vortices in the longitudinal plane inclined slightly towards the hub on the whole. This fact suggested a contraction of the trajectories of the tip vortices, which became more significant when the advance ratio decreased. Since things change when the duct of a certain diameter and angle of attack is placed in front of the propeller, two more investigations are therefore further proposed. The first one aims at clarifying the influence of the duct diameter, whereas the second is focused on the influence of its attack angle.

Figure 17 shows the vorticity computed in the $D_N/D=0.7$ and $D_N/D=1.0$ cases, respectively. In both cases considered here, the attack angle is 10°. To easily emphasize the differences between the two duct diameters, the vorticity contours are superposed. The $D_N/D=0.7$ solution is drawn in cyan color, whereas the $D_N/D=1.0$ solution is in yellow. The helical pitch of the vortices, which was constant in the single propeller case, shows a fairly progressive augmentation in the wake as the nozzle diameter decreases. Obviously, this is due to the axial flow acceleration that is determined by the smaller diameter duct. Since the angle of attack was kept unchanged, no spin was observed for the cores, as shown in the close-up view depicted in Figure 17b, where the horizontal axes of cores are denoted by (I-I) to (IV-IV). For clarity reasons, only the first four tip vertices are considered. The contraction toward the center line of the radial positions of cores found when the single propeller was analyzed is present here as well (curve A-A).

Figure 18 shows the vorticity computed for the $D_N/D=0.7$ duct inclined at three angles of attack, i.e., 0°, 10° and 20°, respectively. To emphasize the differences between the numerical solutions, the vorticity contours are again overdrawn. The 0° solution is drawn in white color, the 10° in cyan and the 20° solution in yellow. The modification of the attack angle generates instabilities both in time and space, mainly in the slip stream region. The increase of the attack angle leads to violent separations, which determines a significant change of the vorticity field at the downstream of the propeller, as shown in Figure 18a. The modification affects not only the tip vortices, but also the root ones. The most prominent consequence of the modification of the attack regards the helical pitch of the vortices, which loses its invariability. When the attack angle increases, the contraction line of the cores moves from the A-A position to B-B, as shown in Figure 18b. Apart from that, a spin rotation of the cores of the vortices can be observed. That is, the axes that connect the cores denoted by (I-I) and (II-II) in Figure 18b rotate back and forth, with further consequences on their further development in space.

Figure 19 proposes a comparison of the vortical structures developed behind the propeller in the computational case of $D_N/D=0.7$ for all the attack angles of the duct considered in the present study. For reasons related to the clarity of graphical representation, only the solution corresponding to the sliding sub-domain whose mesh was drawn in red in Figure 1f is considered. Since the vortices’ topology is defined as iso-surfaces of $Q=15,000,$ as in the single propeller computation, a simple comparison between the solutions reveals that the duct presence determines a significant modification of the vortex propagation in space. As discussed previously, the increase of the attack angle of the nozzle led to a noticeable modification of the velocity and vorticity fields, so the inter-influence between the duct and propeller ultimately determines an attenuation of the vortices’ propagation in space.

The net thrust realized by the propeller has to be analyzed in conjunction with the drag forces exerted by the duct, shaft, hub and blades since they are parts of the flow equation. Based on numerical solutions not reported here, the shaft and hub drags are not significantly dependent on the duct size and its attack angle. Consequently, only the thrust exerted by the blades and the drag of the duct are drawn in Figure 20a,b, respectively. Figure 20a provides the propeller thrust computed for the five attack angles of the nozzle. For reasons related to an accurate quantitative estimation of the duct contribution, the thrust computed for the single propeller denoted by $T_0$ was found as well. As expected, the drag of the nozzle, denoted by $R$ in Figure 20b, increases proportionally with the attack angle regardless of the duct size. Obviously, this is due to the violent flow separation that has been shown in Figure 18a. On the other hand, for a certain angle of attack, the magnitude of the duct drag increases as the nozzle diameter increases, as also shown in Figure 20b.

The net axial force of the corresponding nozzle-propeller system can then be simply obtained by subtracting the drag of the duct from the thrust realized by the blades of the propeller, as depicted in Figure 20c. In contrast to the monotonic increase of the blade thrust and duct drag with the angle of attack, the magnitude of the total axial force has a depletion point, which may be considered to be the critical angle of attack at which the magnitude of the total axial force is maximum. It is worth mentioning that the critical angle of attack provides the largest gain in all of the cases considered here, regardless of the duct diameter. The most efficient duct proved to be the one with $D_N/D=0.7$ and $\theta =10°$, for which the efficiency reported to the thrust of the single propeller working in open water, expressed as $\eta =\left[\left(T-R\right)/T_0\right]\times 100$, amounts to about 10.15%.

4. Concluding Remarks

The present paper described a set of numerical simulations of the 3D viscous flow around a NACA0015 nozzle placed in front of the five-bladed KP505 propeller. Four different diameters and five attack angles were considered for studying the cumulated influence of the duct on the propeller performances. A set of computations was initially carried out at first for the single propeller working in open water. Nineteen advance coefficients ranging from 0.1 to 1 were considered. The finite-volume ISIS-CFD solver of the Fine^TM/Marine package was employed in the reported computations. Closure to turbulence was achieved by using a hybrid DES-SST model.

The rotating frame technique was used for the single propeller computations, whereas the sliding grids method was used in the nozzle-propeller simulations. Grid dependency and numerical uncertainty analyses were performed to check the level of the numerical errors, thus proving the robustness of the solver. For that purpose, a series of more than fifty computations was carried out to clarify not only the particulars of the local flow around the KP505 propeller working in open water either alone or in tandem with the duct, but also to estimate the balance of forces.

The role of an ESD in providing additional thrust, usually observed in the self-propulsion tests performed in the towing tank, was confirmed through a series of hybrid RANSE-LES-based calculations. Numerical solutions using the fully resolved propeller are particularly important for the final verification of the optimal configurations since the influence of the ESD on the propeller performance, which could be significant, is usually neglected in the current design practice. The present study has proven that an accurate analysis of the propeller functioning is finally needed to prove the reliability of the adopted solution. Last but not least, it seems to be necessary for the naval architect to verify that the improvements in terms of propulsive efficiency are not nullified by a bad functioning point of the given propeller, or to assess, consequently, the need of a redesign of the propeller itself to better address the new functioning point and inflow. The influence of the different ESD geometries on propeller performance and self-propulsion coefficients for an advance velocity which corresponds to the design ship speed should necessarily be verified also in off-design conditions since the ESD and the propeller are equally contributing to the foreseen gain in efficiency. With this turnaround at hand, based on the numerical solutions and the associated discussions provided above, the following concluding remarks can be put forward:

• The present research aimed at providing a methodology needed by a naval architect when an efficient ESD is sought for an optimal design not only of the propulsive system, but also for the active surface of the maneuvering system.
• Numerical simulations were conducted in several scenarios in which the nozzle diameter and attack angle were systematically modified within certain limits initially imposed.
• The computation conditions considered here proved that the nozzle with $D_N/D=0.7$ and $\theta =10°$ may lead to an increase of the resulting axial force, ultimately of the propulsion efficiency, of a little more than 10%, compared to the POW case.
• For the considered inflow axial velocity, the most effective nozzles proved to be those with the smaller diameters compared to the propeller diameter.
• Angles of attack larger than 15° did not assure a clear improvement of performances, seemingly because of the significant flow that takes place.
• The above remarks are only valid for the streamwise velocity chosen in this study. Therefore, the extrapolation of the above conclusions to other velocities may be hazardous. This is due to the differences induced by the different transposition criteria used in the dimensional analysis. As far as the Strouhal number is the most relevant parameter for the propeller, for the nozzle, the Reynolds criterion is the key parameter that describes the viscous flow.
• The present study did not take into account other relevant geometric parameters, such as the relative distance between the nozzle and propeller and the chord length of the nozzle. Since those parameters may influence the overall performance of the flow, more studies are required to completely clarify the subject.
• The reported research was focused only on the nozzle-propeller system, without taking into account the ship hull. However, a complete investigation should also consider the hull, knowing its strong influence on the flow that enters in the propulsor.
• Although the research reported here is an unsteady one, the solutions were discussed for the converged solution only. Nevertheless, the time history of the velocity and pressure fields provided by such investigation may be further used for an extended analysis not only on the propulsive performances, but also in the effective rudder design process. Moreover, a detailed description of the flow field behind the propeller may offer the possibility of avoiding the cavitation phenomenon on the active surface of the governing system.