Open-Water Thrust and Torque Predictions of a Ducted Propeller System with a Panel Method

This paper discusses several modelling aspects that are important for the performance predictions of a ducted propulsor with a low-order Panel Method. The aspects discussed are the alignment of the wake geometry, the influence of the duct boundary layer on the wake pitch, and the influence of a transpiration velocity through the gap. The analysis is carried out for propeller Ka4-70 operating without and inside a modified duct 19A, in which the rounded trailing edge is replaced by a sharp trailing edge. Experimental data for the thrust and torque are used to validate the numerical results. The pitch of the tip vortex is found to have a strong influence on the propeller and duct loads. A good agreement with the measurements is achieved when the wake alignment is corrected for the presence of the duct boundary layer.


Introduction
The ability to accurately predict the thrust and torque of a ducted propeller in open-water conditions is very important for a calculation method used in the design stage. RANSE methods have been progressively introduced for the calculation of ducted propeller systems, meeting considerable success in predicting open-water characteristics for the well-known Ka-series [1][2][3]. However, due to their relative complexity and time requirements, they are not yet routinely used in the design process, which is often still based on the use of inviscid flow methods.
Various numerical methods based on inviscid (potential) flow theory have been proposed for the analysis of ducted propellers. Examples are the combination of a Panel Method, also known as Boundary Element Method, to model the duct with a vortex lattice method for the propeller [4], and a Panel Method for the complete ducted propeller system operating in unsteady flow conditions including blade sheet cavitation [5]. Both methods applied a transpiration velocity model for the gap flow between propeller blade tip and duct inner surface and analysed a duct with a sharp trailing edge.
These references show that the application of inviscid flow models to ducted propellers, albeit of great usefulness, may meet some serious limitations related to the occurrence of flow regions where viscous effects cannot be ignored and have to be modelled in some way for the correct prediction of the ducted propeller thrust and torque. One of such region concerns the gap flow, which has a strong influence on the propeller and duct circulation distribution, and therefore, on the distribution of loading between propeller and duct, as studied in detail by Baltazar and Falcão de Campos [6]. In addition, there may be a considerable interaction between the vorticity shed from the propeller blade tips and the boundary layer developing on the duct inner side, as found in the works of Krasilnikov et al. [3] and Rijpkema and Vaz [7]. This effect has not been studied before with potential flow methods, and its importance is therefore unknown.
The purpose of this paper is to show the importance of an efficient and robust method for the vortex pitch alignment, including the effect of gap modelling and the effect of duct boundary layer. The computational results are compared with open-water data measured at MARIN for the ducted propeller Ka4-70 with P/D = 1.0 operating without and  inside a modified duct 19A [8]. The section geometry of this duct, denoted as duct 19Am was obtained by replacing the round trailing edge of the duct 19A [9] by a sharp trailing edge. This sharp trailing edge allows the application of a classical Kutta condition for the prediction of the duct circulation.
Details of the mathematical formulation of the Panel Method are shown in Section 2. The numerical models for the gap region and the interaction with the duct boundary layer are presented in Section 3. This section also discusses the two wake alignment methods that have been investigated, which are a rigid wake model with prescribed geometry and an iterative wake alignment model. The influence of these models on the results and the comparison with the experimental data is shown in Section 4. In Section 5, the main conclusions are drawn.

Mathematical Formulation
2.1. Potential Flow Problem. Consider a propeller of radius R rotating with constant angular velocity Ω inside a duct both advancing with constant axial speed U along its axis in an incompressible ideal fluid of constant density ρ and kinematic viscosity ν otherwise 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 In a moving reference frame advancing and rotating with the propeller, we introduce a Cartesian coordinate system (x, y, z), 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 blade k = 1, and the zaxis completing the right-hand system. We use a cylindrical coordinate system (x, r, θ) related to the Cartesian system by the transformation: (1) Figure 1 shows the coordinate system used to describe the propeller geometry and the fluid domain around the ducted propeller.
In the reference frame moving with the propeller blades, the flow is steady, and assuming that the perturbation International Journal of Rotating Machinery 3 velocity due to the ducted propeller is irrotational then, the velocity field V (x, y, z) may be described by a perturbation potential φ(x, y, z) in the form: where is the undisturbed onset velocity in the moving reference frame.
The perturbation potential satisfies the Laplace equation: The boundary of the domain consists of the blade surfaces S B , the duct surface S D , and the hub surface S H . The kinematic boundary condition, is satisfied on the blade, duct, and hub surfaces, where ∂/∂n denotes differentiation along the normal and n is the unit vector normal to the surface directed outward from the body. At infinity, the flow disturbance due to the ducted propeller vanishes To allow for the existence of circulation around the propeller blades and the duct, vortex sheets are shed from the trailing edge of the blades and the duct. The boundary conditions on the vortex sheet surfaces S W are the tangency of the fluid velocity on each side of the sheet and the continuity of the pressure across the sheet. In steady flow, these conditions are where p is the pressure and the indices + and − denote the two sides of the vortex sheets, taken, respectively, on the side of the back (normally the suction side) and face (normally the pressure side) of the blade at the trailing edge, and on the inner side and outer side of the duct at the trailing edge. In the case of the propeller blade, the unit normal to the vortex sheet is defined pointing from the face (−) to the back (+) side of the blade. In the case of the duct, the unit normal to the vortex sheet is defined pointing from the outer (−) to the inner (+) side of the duct. In order to specify uniquely the circulation around the blades and duct, it is necessary to impose the Kutta condition at the blade trailing edge and at the duct trailing edge. The Kutta condition states that the velocity must remain finite at a sharp trailing edge. In the theoretical formulation, the Kutta condition, (8), is needed to guarantee that the vortex wakes are shed from the sharp trailing edge, and (7) should be applied at the vortex wake including the trailing edge.
The excess pressure p due to the fluid motion (with respect to pure hydrostatic conditions) may be determined from the Bernoulli's equation for the steady flow in the moving reference frame in the form where p ∞ is the pressure of the undisturbed inflow. Applying Green's second identity, assuming for the interior region to S B ∪S D ∪S H , φ = 0, we obtain the integral representation of the perturbation potential at a point p on the body surface, in the wake surfaces S W and the surface dipole strength at the blade and duct trailing edges. In potential flow theory, the dipole strength at a given point on the wake equals the circulation for a closed contour around the lifting surface (blade or duct) through that point.

Velocity, Pressure and Forces.
The velocity on the surface is obtained by differentiation of the surface potential distribution. From Bernoulli's equation, (9), the pressure coefficient C p can be determined from: where The inviscid thrust and torque on the ducted propeller are obtained by integration of the pressure distribution on the blade and duct surfaces. The propeller operation conditions are defined by a single nondimensional parameter: the advance coefficient J = U/nD, where n = Ω/2π is the rate of revolution and D = 2R the propeller diameter. The nondimensional thrust and torque of the ducted propeller system are given by the propeller thrust coefficient K TP , the duct thrust coefficient K TD , and the torque coefficient K Q : (12) where T P is the propeller thrust, T D the duct thrust, and Q the propeller torque. The ducted propeller efficiency is given by

Surface Discretisation.
For the numerical solution of the integral equation, (10), we discretise the blade surfaces S B , the duct surface S D , the hub surface S H , and the wake surfaces S W in bilinear quadrilateral panels. The blade is discretised in the spanwise radial direction by a number of strips, extending chordwise from the blade leading edge to the trailing edge. A cosine formula for the panel spacing in the radial and chordwise directions is used to increase the panel density near the blade extremities. The duct surface is divided into three regions: between, upstream, and downstream of the blades. The regions between and downstream of the blades are discretised along the axial direction with the same panel spacing of the blade tip section and blade wake tip section, respectively. The region upstream of the blades is discretised along the axial direction with a Vinokur stretching function [10]. Equidistant angular spacing is used in circumferential direction. For the discretisation of the hub surface, an elliptical grid generator is used [11]. The blade and duct wake surfaces are discretised in the spanwise and circumferential directions, extending downstream from the trailing edge the corresponding strips on the propeller blade and duct surfaces.

Solution of the Integral Equation.
The integral equation, (10), is solved by the collocation method with the element centre point as collocation point. We assume a constant 3.3. Wake Models. Two wake models are considered: a rigid wake model and a wake alignment model for the blade wake where the pitch of the vortex lines are aligned with the local fluid velocity.
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. However, contraction of the blade wake is prescribed by an empirical formulation, following Hoshino [13] for the open propeller case. The radial variation of the helicoidal lines is occurring in a transition wake region, with a length of 2R, and is kept constant further downstream. For the duct, the wake leaves the trailing edge at the bisector. The value of the dipole strength at the blade and duct trailing edges is determined by the application of an iterative pressure Kutta condition, requiring that the pressure is equal on both sides of the blade and duct panels adjacent to the trailing edge.
In the wake alignment model, the corner points of the blade wake grid panels are displaced with the mean fluid velocity. At the (n+1)th iteration, the geometry in cylindrical coordinates of the wake strip i+1 can be determined by using an Euler scheme: where V x , V r , and V θ are the components of the mean vortex sheet velocity along the axial, radial, and circumferential directions, respectively, and Δt is the time step for the Euler vortex convection scheme. The velocity components are calculated from the integral equation of the velocity,   obtained by taking the gradient of (10). The velocities are evaluated at a staggered grid of control points, which are located equidistantly from consecutive dipole singularities of the wake. The first wake strip (i = 1) corresponds to the blade trailing edge. To control the wake alignment stability, the radial coordinates of the blade wake grid are kept constant during the iterative process. Hence, where The new axial and circumferential coordinates x (n+1) i+1 and θ (n+1) i+1 are calculated using an Euler scheme: The nondimensional time step Δθ = ΩΔt is introduced, which can also be expressed in terms of the number of time steps per propeller revolution N θ = 2π/Δθ. Note that in the wake alignment model, the downstream part of the duct is remeshed at each iteration of the propeller and duct vortex wakes alignment.

Gap Flow Models.
Two different models for the potential flow in the gap region are considered: a closed gap with zero gap width and a gap flow model with transpiration velocity.
In the closed gap model, the blade and duct grids match in the tip region. In this case, the flow is not allowed to pass between the blade tip and the duct inner surface. The blade wake connects to the duct inner surface down to the duct trailing edge.
In the gap model with transpiration velocity, a partial flow between the blade tip and the duct inner surface is allowed to pass in the gap region. The procedure proposed by Hughes [14] is implemented in the current method, where the gap flow is treated as a two-dimensional orifice. The velocity in the gap region can be related to the difference in pressure across the blade tip using Bernoulli's equation. The reduction in the flow from losses in the orifice is expressed in terms of an empirical discharge coefficient from the total volume flow rate through the gap Q G and the pressure difference across the gap Δp as A typical value as used by Hughes [14] is C Q = 0.84. From (18), the mean velocity through the gap at a given chordwise location can be expressed in terms of the difference in pressure between the pressure and suction sides of the blade at that location: where ω R is the mean relative flow velocity through the gap and ΔC p is the pressure difference at the blade tip in terms of the pressure coefficient defined by (11).
To incorporate this expression for the mean velocity into the Panel Method, an additional strip of panels along the blade tip is introduced which closes the gap. This wake strip continues past the duct trailing edge, connecting the duct wake with the blade wake. The dipole strength on the last wake strip is also determined from the Kutta condition. In the case of a closed gap, the boundary condition on the gap panels sets the source strength to cancel the normal component of the inflow velocity, (5). To allow for the existence of a transpiration velocity through the gap, the boundary condition on the gap strip becomes where n c is the unit normal vector to the mean camber line at the gap strip on the same chordwise position of the panel.
The value of ΔC p is obtained by an iterative procedure. First, the potential flow problem is solved as if the gap were completely closed and (5)  kinematic boundary condition on the gap strip. The pressure distribution is then computed by differentiating the potential using a finite difference scheme. The potential flow problem is solved again with the kinematic boundary condition on the gap strip specified by (20). The pressures on the gap strip are then recomputed and used to update the boundary condition on the gap panels. The solution process will be repeated until a specified convergence criterion is met.

A Simple Model for the Interaction of the Blade Wake with the Duct Boundary Layer.
A simple model for the interaction between the blade wake and the boundary layer on the duct inner side is implemented as well. Due to the duct boundary layer, a reduction in the axial velocity may be taken into account in the convection of blade vorticity in the wake alignment model. Considering δ as the duct boundary layer thickness and assuming a power law distribution for the velocity profile, we have To correct the potential flow velocities in the blade wake, a linear variation of the axial velocity is considered in the gap region, R ≤ r ≤ R d . The corrected axial velocity at the duct surface, V x (0), is calculated by extrapolation from its value at the gap, V x (h), and at the edge of the duct boundary layer, V x (δ). The velocity at the edge of the duct boundary layer is identical to the velocity on the duct surface in the potential flow method. A correction factor for the velocity at the duct surface can be derived from (21) and is given by  found in Kuiper [15]. The section geometry of the duct 19Am is given in Bosschers and van der Veeken [8], which also contains the results of the experiments performed at MARIN in the Netherlands. The gap between the duct inner surface and the blade tip is uniform and equal to 0.83% of the propeller radius. The profiles of the duct 19A and 19Am are compared in Figure 2. The comparison of the experimental open-water characteristics of the ducted propeller system with duct 19A [15] and with duct 19Am [8] is presented in Figure 3. Figure 4 shows a typical panel arrangement of propeller Ka4-70 inside the duct 19Am with an open-water hub. Calculations were carried out with the rigid wake model and wake alignment model. The wake geometries were obtained after 5 iterations of the wake alignment model using N θ = 90 time steps per revolution, which leads to an angular step of 4 degrees. For the ducted propeller case, a transpiration velocity gap model with a discharge coefficient equal to C Q = 0.84 [14] and a closed gap model were considered. A specified tolerance of |δ(ΔC p ) GAP | < 10 −2 on the gap strip control points was applied as convergence criterion, corresponding to approximately 20 iterations for the transpiration velocity.
In general, the iterative pressure Kutta condition converged after 5 iterations to a precision of |ΔC p | ≤ 10 −3 . In the iterative procedure to deal with nonlinearities of the numerical method, the outer loop is the wake alignment model, the middle loop is the transpiration velocity gap model, and the inner loop is the iterative pressure Kutta condition.

Ka4-70 without Duct.
First, the numerical results are presented for the Ka4-70 propeller without duct in uniform flow. In this case, the rigid wake model and the same wake vortex pitch alignment model as used for the ducted propeller is applied in combination with, or without, an empirical contraction. Figure 5 shows the comparison of the thrust and torque coefficients with experimental data from the open-water tests. A section viscous drag coefficient of 0.007 and a lower limiting value for the chordwise component of the blade section force between sections r/R = 0.7 and r/R = 1.0 were used for all computations. The limiting value models the effect of flow separation which eliminates the unphysical suction peaks at the leading edge in the potential flow theory. The limiting value has been set Table 1: Variation with the duct boundary layer thickness of the relative differences between the predicted and experimental (exp) thrust and torque coefficients. Ka4-70 inside duct 19Am at J = 0.5. Transpiration velocity gap model. to zero. For a flat plate, this would imply the suppression of the leading edge suction force. By applying a leading edge suction force correction, a significant increase of the torque coefficient at low advance ratio and a decrease in the thrust coefficient are obtained. Small differences between the wake models are seen with the exception of the low advance ratios. In this case, a reduction of the force coefficients is observed with wake alignment models for low advance ratios. These differences are related to the difficulty in obtaining converged and smooth aligned wake surfaces when large induced velocities are present. In general, a fair-to-good agreement is seen between the numerical predictions and the experimental data. The reasons for the underprediction of the force coefficients at low advance ratios are probably the limitations of the wake alignment model and the flow separation from the propeller blades, which is modelled by a simple suppression of the suction force. Figure 6 presents the panel arrangement including the aligned wake for the advance ratio J = 0.5. To control the convergence and smoothness of the wake grid, the tip vortex roll-up phenomenon is not modelled with the present alignment method. Figure 7 shows the blade circulation obtained from the wake dipole strength and the chordwise pressure distribution at the radial section r/R = 0.95. Similar circulation and pressure distributions are obtained between all wake models. Please note that large suction peaks near the blade leading edge are obtained in the inviscid flow solution, which are physically not realistic.

Ka4-70 inside Duct 19Am.
Secondly, the effect of the gap and the wake alignment models are investigated for propeller Ka4-70 inside duct 19Am. Calculations are presented for the rigid wake model and for the wake alignment model without and with correction of the axial velocity due to the duct boundary layer. The thickness of the duct boundary layer δ/R is roughly assumed to be 4%, corresponding to h/δ = 0.21. A power law velocity profile with n = 7 is also assumed. The boundary layer thickness is a critical parameter and, since there is yet no information on suitable values to be used in this model, a sensitivity study is presented at the end of this section. Figure 8 shows the comparison of the predicted thrust and torque coefficients with experimental data from the open-water tests. A section viscous drag coefficient of 0.007 was used for all computations for the propeller. Lacking any adequate method to estimate the duct viscous drag under the influence of a propeller, no viscous drag correction to the duct thrust has been applied in the present comparison. The differences between the closed and transpiration gap models are very small in this case. The comparison with the experiments shows that the numerical results obtained with the rigid wake model and wake alignment model without duct boundary layer correction over-predict the propeller thrust and torque coefficients. By applying the axial velocity correction at the tip in the wake alignment model, a significative reduction of the propeller thrust and torque coefficients is observed. A good agreement with the measurements may be achieved. The duct thrust coefficient agrees well with the measurements for low advance coefficients. The correction for the duct boundary layer results in an increase of the duct thrust. However, for high advance ratios, all models over-predict the duct thrust, which is probably due to the occurrence of flow separation on the outer side of the duct, and it is not modelled in the present method. The differences between the measured and the predicted ducted propeller efficiencies are around 6% for the advance ratios lower than J = 0.6. Figure 9 presents the panel arrangement including the aligned wake with the duct boundary layer correction. A strong reduction of the blade wake pitch is seen at the tip. The pitch angle of the blade wake at the axial positions x/R = 0.5 and x/R = 2.0 is illustrated in Figure 10. The correction in the axial velocity due to the duct boundary layer introduces a large reduction of the vortex pitch at the blade wake tip, which is responsible for the significant decrease of propeller thrust and torque and increase of duct thrust. Figure 11 presents the circulation distributions on the propeller blade and duct. The duct circulation is shown as function of the angular position between two consecutive blade wakes at the duct trailing edge. The pressure distributions are shown in Figures 12 and 13 for the gap model with transpiration velocity and closed gap model, respectively. The pressure distributions are shown as function of the chordwise position s/c, where s is defined along the cylindrical section chord for the propeller blade and as the axial distance along the meridional duct section, from the leading edge to the trailing edge. In the closed gap model, the blade circulation at the tip is equal to the duct circulation discontinuity. Note that a finite circulation at the blade tip is also obtained in the gap model with transpiration velocity. Similar results are obtained with the two gap models. Small differences are seen in the blade tip  and duct circulations. In the closed gap model, the flow is not allowed to pass through the gap and increases the blade tip and duct circulations. For the pressure distributions and propeller forces, similar results are obtained between the two gap models. From the comparison between the wake models, a decrease in the blade circulation is obtained with the wake alignment model, which has the opposite effect on the duct circulation. Although the correction due to the duct boundary layer affects only the blade wake geometry near the tip, a very strong influence on the propeller and duct circulation distributions is obtained. For the blade pressure distribution, no suction peak is seen in the wake alignment model with duct boundary layer correction, indicating a different incidence angle to the blade section. The differences in the duct pressure distributions on the duct inner side in front of the propeller explain the differences in duct thrust. A sensitivity study for the influence of the duct boundary layer thickness on the predicted thrust and torque coefficients has been performed using the transpiration gap model. The relative differences between the numerical results and the experimental open-water thrust and torque values are presented in Table 1 for various values of the duct boundary layer thickness. It is observed that the influence is very significant and that an increase in boundary layer thickness always leads to a reduction of propeller thrust and torque.

Conclusions
A low-order Panel Method has been used to predict the thrust and torque of propeller Ka4-70 without and inside a modified duct 19A in open-water conditions. Calculations were carried out with a rigid wake model and an iterative wake alignment model for the blade wake pitch. For the open propeller, a fair-to-good agreement is seen between the numerical predictions and the experiments for all wake models. For the ducted propeller, the loading predictions of the duct propeller system were found to be critically dependent on the blade wake pitch, especially at the tip. A model to incorporate the influence of the duct boundary layer on the pitch of tip vorticity was found necessary to correctly predict the open-water characteristics of the ducted propeller. This mechanism was not studied before using panel methods but seems plausible from the physical point of view. The quantitative aspects of this mechanism will be investigated in more detail by comparing the results of the present panel method with RANS simulations.