Prediction of the Open-Water Performance of Ducted Propellers with a Panel Method

João Manuel Baltazar 1,* ID , Douwe Rijpkema 2, José Falcão de Campos 1 and Johan Bosschers 2 1 Department of Mechanical Engineering, Instituto Superior Técnico, Av. Rovisco Pais 1, Universidade de Lisboa, 1049-001 Lisboa, Portugal; falcao.campos@tecnico.ulisboa.pt 2 Maritime Research Institute Netherlands, 2 Haagsteeg, 6708 PM Wageningen, The Netherlands; d.r.rijpkema@marin.nl (D.R.); j.bosschers@marin.nl (J.B.) * Correspondence: joao.baltazar@tecnico.ulisboa.pt; Tel.: +351-218-419-289


INTRODUCTION
The calculation of the flow around a ducted propeller system in open-water with a panel code has been the subject of investigation by IST and MARIN, see Baltazar and Falcão de Campos (2009), Baltazar et al. (2011Baltazar et al. ( , 2013Baltazar et al. ( , 2015)).In these studies a low-order panel method has been used to predict the open-water diagram of a ducted propeller system, where several modelling aspects have been analysed for the improvement of the inviscid (potential) flow solution.The investigation comprehended the influence of the Kutta condition, gap flow model and wake model on the performance prediction.A similar method has also been implemented by MARIN in a in-house panel code (Bosschers et al. 2015).
Since the blade trailing-edge has a sharp geometry and the duct has a thick round trailing-edge, an alternative Kutta condition was proposed.In this new Kutta condition, the chordwise location for pressure equality on both sides of the duct trailing-edge has a strong influence on the propeller force prediction.From the comparison between the closed gap model and the gap model with transpiration ve-locity based on the work presented by Hughes (1997), similar results were found.However, the loading predictions of the ducted propeller system were found to be critically dependent on the blade wake pitch, especially at the tip.In this way, a simple model for the interaction between the blade wake and the boundary layer on the duct inner side was implemented in combination with the wake alignment model (Baltazar et al. 2011(Baltazar et al. , 2015)).With this model, a reasonable to good agreement was obtained between the inviscid predictions and the experimental data from openwater tests (Baltazar et al. 2015).However, significant differences were still seen in the open-water predictions at low advance ratios.
It is known that the prediction of the propeller performance at bollard pull conditions is important in the design stage.In the present work, a comparison between the results obtained by the panel method with RANS calculations is made to obtain a better insight on the viscous effects of the ducted propeller and on the limitations of the inviscid flow model.The comparison focus mainly at the low advance ratios.The paper is organised as follows: a description of the numerical methods is given in Section 2; the comparison of the inviscid predictions with the RANS calculations and the experimental open-water data is presented in Section 3; in Section 4 the main conclusions are drawn.

Choice of Coordinate System
Consider a propeller of radius R rotating with constant angular velocity Ω inside a duct and advancing with constant axial speed U along its axis in an incompressible ideal fluid of constant density ρ at rest in a domain extending to infinity in all directions.The propeller is made of K blades symmetrically distributed around an axisymmetric hub.The duct is also considered to be axisymmetric of inner radius at the propeller plane R d ≥ R which defines a gap height h = R d − R. The flow field is steady in a reference frame rotating with the propeller blades around its axis.Figure 1 shows the coordinate system used to describe the propeller geometry and the fluid domain around the ducted propeller.We introduce a Cartesian coordinate system (x, y, z) rotating with the propeller blades, with the positive x-axis direction opposite to the propeller axial motion, the y-axis coincident with the propeller reference line, passing through the reference point of the root section of the reference blade, and the z-axis completing the right-hand-system.The unit vectors in this system are denoted by ( e x , e y , e z ).We use a cylindrical coordinate system (x, r, θ) with unit vectors ( e x , e r , e θ ) related to the Cartesian system by the transformation y = r cos θ, z = r sin θ. (1) The undisturbed onset velocity in the rotating frame is (2)

Panel Code PROPAN
PROPAN is a IST in-house panel code which implements a low-order potential-based panel method for the calculation of the incompressible potential flow around marine propellers.The code has been widely used in the calculation of the three-dimensional potential flow around ducted propellers, (Baltazar and Falcão de Campos 2009, Baltazar et al. 2011, 2013, 2015).
Applying Green's second identity, assuming for the interior region to the blade surface S B , duct surface S D and hub surface S H that the perturbation potential is equal to zero, we obtain the integral representation of the perturbation potential φ at a point p on the body surface, where G (p, q) = −1/R (p, q), R (p, q) is the distance between the field point p and the point q on the boundary S B ∪ S D ∪ S H ∪ S W .With the ∂φ/∂n q on the surfaces S B , S D and S H known from the Neumann boundary condition on the body surface, the Equation ( 3) is a Fredholm integral equation of the second kind in the dipole distribution µ(q) = −φ(q) on the surfaces S B , S D and S H .The Kutta condition yields the additional relationship to determine the dipole strength ∆φ(q) in the wake surfaces S W .
For the numerical solution of Equation (3), we discretise the blade surface S B , the duct surface S D , the hub surface S H and the wake surfaces S W in bi-linear quadrilateral elements which are defined by four points in the surface.The integral equation, Equation (3), is solved by the collocation method with the element centre point as collocation point.We assume a constant strength of the dipole and source distributions on each element.The influence coefficients are determined analytically using the formulations of Morino and Kuo (1974).To reduce the computational time, far field formulas are also used in the calculation of the influence coefficients.
To determine the dipole strength of the wake an iterative pressure Kutta condition at the blade and duct trailing edges is applied.However, different forms of the Kutta condition are considered, since the blade trailing-edge has a sharp geometry, whereas the duct trailing-edge presents a blunt round geometry.For the blade trailing-edge, equal pressure on the collocation points of the panels adjacent to the trailing-edge on the upper and lower sides is considered.
For the duct trailing-edge, the chordwise location of pressure equality on both sides is specified, which controls the strength of the shed vorticity.Due to the possible occurrence of flow separation, a constant pressure distribution downstream of the Kutta points is assumed in this model.Note that the potential flow solution at the duct trailing edge satisfies the integral equation, Equation ( 3), but the corresponding pressure distribution is disregarded aft the Kutta points.A detailed description of the iterative Kutta condition model can be found in Baltazar et al. (2015).
For the wake geometry, two wake models are considered: a rigid wake model (RWM) and a wake alignment model (WAM).In the rigid wake model, the geometry of the wake surfaces is specified empirically.For the blade wake, the pitch of the vortex lines is assumed constant along the axial direction and equal to the blade pitch.For the duct, the wake leaves the trailing edge at the bisector.In the wake alignment model, the blade wake pitch is aligned with the local fluid velocity, while the radial position of the vortex lines is fully prescribed.For the modelling of the interaction between the blade wake and the boundary layer on the duct inner side, a correction to the blade wake pitch near the tip is introduced.The gap between the blade tip and the duct inner surface is modelled as a rigid surface, named in this study as the gap strip, and where the transpiration velocity is neglected.In this closed gap model, the boundary condition on the gap panels sets the source strength to cancel the normal component of the inflow velocity lead-ing to the Neumann boundary condition, Equation ( 4).The gap strip extends from the blade tip to the duct inner side.
The description of the wake and gap models can be found in Baltazar et al. (2015).

RANS Code ReFRESCO
ReFRESCO (www.For open-water (steady) calculations the equations are solved in the body-fixed reference frame which is rotating with velocity Ω.For the boundary conditions a uniform flow at the inlet and constant pressure at the outer boundary is applied.At the outlet an outflow condition of zero downstream gradient is used.For the propeller blades, duct and hub, a non-slip boundary condition is set.A rotational velocity is prescribed to the blades and shaft, while the duct does not rotate.For all RANS calculations presented in this paper, the k − ω SST 2-equation eddy-viscosity model proposed by Menter Menter (1994) is used.A second-order convection scheme (QUICK) is used for the momentum equations and a first-order upwind scheme is used for the turbulence model equations.A fine boundary resolution is applied to obtain y + ∼ 1.

General
Results are presented for the propeller Ka4-70 inside the duct 19A.The propeller Ka4-70 is a four-bladed propeller of the Kaplan type.A pitch-diameter ratio of P/D = 1.0 is considered.The duct 19A has a length-diameter ratio of 0.5.The gap between the duct inner side and the blade tip is uniform and equal to 0.8% of the propeller radius.The particulars of the Ka-series and the duct section geometry can be found in Kuiper (1992).
Figure 2 shows the grids used for the calculations with panel code PROPAN and RANS code ReFRESCO.The discretisations of the PROPAN grid are: 50×26 on each blade, 190×160 on the duct and 75×80 on the hub.For the wake alignment model considered in the inviscid computations, the wake geometry is obtained after 5 iterations using an angular step of 4 degrees.In combination with the wake alignment model, a duct boundary layer correction is applied (Baltazar et al. 2011).An iterative pressure Kutta condition is applied which, in general, converged after 3 iterations to a precision of |∆C p | ≤ 10 −3 .
The RANS results presented in this report referred to a grid with 26.8 million cells.The ducted propeller is tested for a range of advance coefficients between 0.0 and 0.8 corresponding to Reynolds numbers from 3.5×10 5 to 3.8×10 5 , where the Reynolds number is defined based on the propeller blade chord length at 0.7R and the resulting onset velocity at that radius.Other used quantities are the vorticity ω = ∇ × V and the pressure coefficient: where V is the fluid velocity vector, p is the static pressure, p ∞ the undisturbed static pressure, ρ the fluid density and The propeller operating conditions are defined by the advance coefficient J = U/(nD), where n = Ω/2π is the rate of revolution and D = 2R the propeller diameter.The non-dimensional thrust and torque of the ducted propeller system are given by the propeller thrust coefficient K T P , the duct thrust coefficient K T D and the torque coefficient K Q : where T P is the propeller thrust, T D the duct thrust and Q the propeller torque.The ducted propeller efficiency is given by η = U (T P + T D )/(2πnQ).

Comparison Between PROPAN and ReFRESCO
This section presents the comparison of the calculations with the panel code PROPAN and with the RANS code ReFRESCO for the ducted propeller.The comparison is presented for three advance coefficients: J = 0.1, 0.2 and 0.5.The advance coefficients J = 0.1 and 0.2 refer to highly loaded conditions from where significant differences are still obtained with the present model, see Baltazar et al. (2011Baltazar et al. ( , 2015)).The inviscid thrust and torque is compared with the experimental data in Table 1.Results are presented for the rigid wake model (RWM), the wake alignment model (WAM) and the wake alignment model with a reduced pitch for the gap strip.
Significant differences are seen for the propeller thrust and torque at low advance ratios for both the rigid and aligned wake models.The present wake alignment model includes a correction in the axial velocity due to the interaction between the blade wake and the duct boundary layer.In the present calculations a power law velocity profile with exponent equal to 1/7 is considered and a thickness equal to δ/R = 4% is assumed.In this study estimation of the boundary layer thickness from the ReFRESCO computations has not been done due to the complex interaction between the blade tip vortex and duct boundary layer.However, the differences in the propeller forces suggest that higher corrections to the blade wake pitch may be needed and a new case is analysed, where the gap strip is rotated from the leading edge in order to reduce its pitch.For the present calculations the pitch of the gap strip is assumed to be equal to P/D = 0.9, whereas the blade pitch is constant and equal to P/D = 1.0.In Figure 3 a detail of the gap strip and the blade wake sheet is presented.We note that the gap strip is modelled as a rigid surface and is disconnected from the wake alignment model.The blade wake pitch near the tip may be controlled by the duct boundary layer correction.However, for low advance ratios large corrections are needed and this has led to divergence of the Kutta condition and non-smooth surface grids.Therefore, the reduction of the gap strip pitch has proven to be a robust technique and can be applied at low advance ratios.As expected, a higher reduction in the propeller thrust and torque is obtained with the WAM using a reduced pitch for the gap strip and approaches the results to the experimental data (Table 1).
For a detailed analysis between the inviscid and RANS calculations, the wake geometry is compared with the vorticity field, and the blade and duct pressure distributions are presented for the three advance ratios in Figures 4 to 10.A comparison between the viscous vorticity field and the three inviscid wake models at the planes x/R = 0.2 and x/R = 0.4 downstream from the propeller is presented.A reduction in the pitch of the tip vortex is seen when changing from the rigid wake model to the wake alignment model.However, the assessment of the correct location of the tip vortex core from the ReFRESCO calculations is difficult to make due to the interaction between the tip vortex and the duct boundary layer, creating a region a viscous flow region at the duct inner side.
The comparison of the blade and duct pressure distributions are presented along the chordwise direction s/c at the radial section r/R = 0.95 and circumferential position θ = 0 degrees, respectively.For the advance ratios J = 0.1 and 0.2, an improvement in the agreement between the inviscid and viscous pressure distributions is obtained when us- ing the wake alignment model with reduced pitch for the gap strip.This comparison shows the influence of the tip vortex pitch on the prediction of the pressure distribution, especially on the duct inner side downstream of the propeller.A reduction in the suction peak at the blade leading edge is also visible, since the blade wake sheet strongly affects the local flow direction.For the pressure distribution at J = 0.5, a good agreement is achieved with the wake alignment at the blade suction peak and duct pressure distribution downstream from the propeller.In this case, no significant improvements are obtained when combining the wake alignment model with the reduced gap pitch.For this advance coefficient (J = 0.5), the assumption of a duct thickness equal to δ/R = 4% is sufficient for the correct prediction of the propeller and duct load.However, two exceptions are observed in the comparison of the pressure distributions: at the blade leading edge near the tip and in the duct inner side.For the blade pressure, a larger suction peak is obtained with the inviscid model.This suction peak decreases with the reduction of the blade wake pitch at the tip.For the duct pressure at the inner side, local pressure minima are observed in the viscous computations, which are related to the passage of the tip vortices from the different blades.This effect is not captured by the inviscid calculations due to the gap model, where the blade wake is attached to the duct inner side.The correlation between the position of the tip vortex core and the peaks of low pressure is illustrated in Figure 10 for the plane z = 0, which cor- responds to the circumferential position θ = 0 degrees.In this figure the inviscid wake geometries are compared with the viscous total vorticity field along the longitudinal direction.A good agreement is obtained with the aligned wakes, except near the blade wake tip, where some differences are still observed.

Comparison With Experimental Data
In this section the predicted thrust and torque coefficients are compared with experimental data available from openwater tests (Kuiper 1992).The inviscid calculations are presented for the wake alignment model assuming a constant duct boundary layer (δ/R = 4%) for all advance ratios.In addition, the wake alignment model is also applied considering a reduced pitch for the gap strip equal to P/D = 0.9.The ReFRESCO calculations are also in- cluded in the comparison.Figure 11 illustrates the comparison of the thrust and torque coefficients with the experiments.A section viscous drag coefficient of 0.007 and suppression of the chordwise component of the blade section lift are considered for all inviscid computations.This suppression models the effect of flow separation which eliminates the non-physical suction peaks at the leading edge in the potential flow theory.No viscous drag correction to the duct thrust has been applied.As expected, the propeller thrust and torque are well predicted for the advance ratios higher than 0.3 when using the wake alignment model without gap pitch correction.For the advance ratios lower than 0.3 a significant improvement in the comparison with the experiments is obtained with the wake alignment model using a reduced pitch for the gap strip.For example, at J = 0.1 the differences between the measured and the predicted propeller thrust reduce from 20% to 7% by applying the gap pitch reduction.A smaller influence of the gap pitch is observed for the higher advance coefficients, which is due to the decrease of the tip vortex strength.The duct thrust coefficient agrees well with the measurements for low advance coefficients.For high advance coefficients, an over-prediction of the duct thrust is seen, which is probably due to the occurrence of flow separation on the outer side of the duct and it is not modelled in the inviscid method.
A better agreement of the propeller forces with the experimental data is obtained with the viscous calculations using code ReFRESCO for advance ratios up to 0.7.In this range the agreement is in the order of 1%.

CONCLUSIONS
In this paper a comparison of panel method with RANS for open-water performance predictions of a ducted propeller system is made.The investigation focused on the improvement of the inviscid performance predictions near bollard pull conditions, which are important in the design of such systems.Therefore, the comparison with RANS computations is made to obtain a better insight on the viscous effects and on the limitations of the inviscid flow model.The comparison of the blade and duct pressure distributions shows that the obtained blade wake pitch near the tip using the wake alignment model with duct boundary layer correction is not sufficient to correctly predict the propeller loads.Since the loading predictions of the duct propeller system are critically dependent on the blade wake pitch especially near the tip, an additional correction is proposed for the gap strip pitch.In this way, a reduction in the blade wake pitch is imposed and the agreement of the propeller forces with the experiments is improved significantly.This technique has proven to be robust, allowing a better control of the quality of the blade wake and duct grids, which is crucial for convergence of the wake alignment model and Kutta condition.
refresco.org) is a community-based open-usage CFD code for the maritime world.It solves the multiphase (unsteady) incompressible RANS equations, complemented with turbulence models, cavitation models and volume fraction transport equations for different phases,(Vaz et al. 2009).The equations are discretised using a finite-volume approach with cell-centred collocation variables.A strong-conservation form and a pressurecorrection equation based on the SIMPLE algorithm is used to ensure mass conservation.The implementation is facebased, which allows grids with elements consisting of an arbitrary number of faces and hanging nodes.The code is parallelised using MPI and sub-domain decomposition, and runs on Linux workstations and HPC clusters.For turbulence modelling, RANS/URANS, SAS and DES approaches can be used.ReFRESCO is currently being developed within a cooperation lead by MARIN.

Figure 2 :
Figure 2: Overview of the surface grids used for the inviscid calculations with PROPAN (top) and RANS calculations with ReFRESCO (bottom).

Figure 3 :
Figure 3: Detail of the gap strip and blade wake sheet.WAM with reduced gap pitch equal to P/D = 0.9.

Figure 11 :
Figure 11: Open-water diagram.Comparison between numerical and experimental results.

Table 1 :
Inviscid thrust and torque coefficients for J = 0.1, 0.2 and 0.5 and comparison with experimental data.Influence of the wake model.