Experimental and Numerical Validation of the Improved Vortex Method Applied to CP745 Marine Propeller Model

Abstract The article presents a numerical analysis of the CP745 marine propeller model by means of the improved vortex method and CFD simulations. Both numerical approaches are validated experimentally by comparing with open water characteristics of the propeller. The introduced modification of the vortex method couples the lifting surface approach for the propeller blades and the boundary element method for the hub. What is more, a simple algorithm for determination of the propeller induced advance angles is established. The proposed modifications provide better results than the original version of the vortex method. The accuracy of the improved method becomes comparable to CFD predictions, being at the same time a few hundred times faster than CFD.


INTRODUCTION
The marine propeller design process is always iterative [4,7,8,3,11]. There are a lot of efficient algorithms for preparing correct designs. Older methods have relied on the lifting line theory, supplemented with pre-calculated correction factors for three dimensional effects [7,8], while more recent and more sophisticated algorithms involve lifting surface calculations to determine the blade pitch and camber directly [4,11]. Despite the fact that these methods are capable of giving good results, they still need validation of the designed geometry through determination of its hydrodynamic characteristics. This can be achieved either numerically or experimentally. Even most classical vortex methods [4,7] are able to provide reliable results for the design point of the propeller. However, it is sometimes important to know also the propeller performance in off-design conditions. This can be achieved experimentally [6,7,8] or numerically [4,7,9,13,11]. Model scale experiments are most reliable, but at the same time extremely expensive. Moreover, they need manufacturing of the propeller model. This is the reason why numerical methods are widely used -at least at the design stage.
Vast experience has been gained over the years on the application of ideal fluid models to the marine propeller design [4,9,13,3,11,19]. They are computationally cheap and efficient. However, even relatively complex vortex models with iterative relaxation of vortex wake [4,9,12,13] do not include true viscous effects, such as circulation reduction and frictional drag inclusion. As a result, these effects have to be taken into account either by artificial formulas, or by boundary layer calculations. However, this approach is unclear and may lead to significant errors -especially at high angles of attack, when flow separation may occur [4]. In such cases, more advanced methods, based on viscous fluid models, are capable of providing better solutions [15,16].
This paper presents results of experimental (Section 2) and numerical analyses of the CP745 propeller model. The propeller was designed and tested at the Ship Design and Research Centre [10] for cooperation with the upstream pre-swirl stator. Numerical analyses were performed by means of the improved vortex model and CFD approach. Vortex models are frequently used for determining hydrodynamic characteristics of the marine propeller. What is important, the presence of the hub is neglected in most applications of the lifting surface method for marine propellers [4,9,11]. The improved vortex method introduced in this paper (Section 3) couples the lifting surface approach for the propeller blades and the boundary element method for the hub. This approach is believed to be reasonable, as it requires much less computational effort than the full boundary element method applied both to the blades and hub, and at the same time still provides accurate results. Furthermore, a simple algorithm to determine the propeller induced advance angles with the use of the vortex method is introduced.
The latter numerical approach is performed by means of CFD calculations (Section 4). These calculations are much more time consuming, but they give more detailed results. In particular, they provide an opportunity to calculate points for lowest advance coefficients for which the vortex model calculations are unable to converge. Finally, both numerical methods and their results are described and compared with experimental data (Section 5).

OPEN WATER TEST GEOMETRY
The controllable pitch propeller model CP745 shown in Fig. 1 was subject to the tests conducted in the Ship Design and Research Centre towing tank.
The propeller was designed to cooperate with the upstream stator [10]. Its diameter (model scale) is D = 226.0 mm, the number of blades is Z = 5, the expanded area ratio is EAR = 0.757, the mean pitch ratio is P mean /D = 0.750 and the total skew angle is SKA = 15.85 ˚ . The propeller model was manufactured in accordance with the ITTC Recommended Procedures and Guidelines [5] at the Ship Hydrodynamics Division of the Ship Design and Research Centre.

COEFFICIENTS
The so-called open water characteristics of the propeller are typically expressed by means of the thrust coefficient K T , torque coefficient K Q , and efficiency η as functions of the advance coefficient J. The thrust coefficient is given by , (1) where T represents the propeller thrust, D stands for the propeller diameter, and n expresses the propeller rotational speed in terms of revolutions per second. The torque coefficient is defined as , (2) where Q represents the propeller torque. The efficiency of the propeller is given by , where the advance coefficient J is expressed by means of the propeller advance velocity V A , namely .

EXPERIMENTAL RESULTS
The open water test to determine propeller hydrodynamic characteristics was conducted by means of the Froude's method [8]. The propeller revolution rate n was constant during the test and equal to n = 21 rps. Changes of the advance coefficient J were achieved by setting a suitable towing carriage speed. Typically, the model scale ratios in propeller tests often exceed 10, as a result of which noticeable differences between propeller characteristics determined in the model scale and those for the real propeller might occur. In order to account for this effect, the test is conducted for at least two values of Reynolds number, namely for Re at which the self-propulsion test is conducted (but not less than 200 000) and for Re as high as possible [6]. The testing experience gained at the Ship Hydrodynamics Division of Ship Design and Research Centre says that it is sufficient to conduct open water tests at the Reynolds number value exceeding 500 000. However, preserving sufficiently high Reynolds number values is not enough to ensure the adequacy of results, and corrections due to the difference between the model and full-scale Reynolds numbers are to be applied. For typical propellers, the standard ITTC-78 method [6] may be utilised. During the tests, the distance between the water free surface and the propeller shaft axis should be equal to or greater than 1.5D.
The measurements were carried out in the Ship Design and Research Centre towing tank, the length × breadth × depth dimensions of which are 270 × 12 × 6 m, respectively. The maximal velocity of the test carriage is 12 m s −1 , and the basic measuring device used in the test is the dynamometer H29. Its measuring capacity is 400 N of thrust force and 15 Nm of torque. As for the maximal rate of propeller revolution, it is up to 25 rps.
The shell of the dynamometer (Fig. 2) is shaped in such a way as to preserve the lowest possible form drag. Prior to the test, control runs with a dummy hub without blades were conducted in order to determine the reference level for the test. Figure 3 presents the experimental open water characteristics of the CP745 propeller model. Points represent actual measurements, whereas solid lines correspond to fourth order polynomial interpolations. The experimental characteristics shown in Fig. 3 were used as the reference for comparing numerical predictions discussed further in Sections 3, 4 and 5.

LIFTING SURFACE CALCULATIONS LIFTING SURFACE CALCULATIONS
Propeller blades are discretised in a classic way, typical for lifting surface calculations [7,4,9,13], so the hydrodynamic singularities are not located on real blade surfaces, see Fig. 4. Instead, they are positioned on the propeller camber surface which would be the real surface provided that the propeller blades are infinitesimally thin. The finite thickness of the blades is modelled with proper distribution of sources, the intensities of which are determined using the stripwise applied thin profile theory [7,4].
The hydrodynamic load effects are modelled with a set of radial bound vortices supplemented with a proper system of free vortices [7,4,9]. The circulation distribution over the blade has to satisfy the kinematic boundary condition on the lifting surface. If the marine propeller is moving forward with constant velocity V A in the uniform ideal fluid, and rotates around its axis with constant angular velocity ω, the kinematic boundary condition is expressed as: , (5) where is the velocity induced at point i by all singularities in the system and can be expressed as , (6) In the above formula, A ij is the unit velocity induced by j-th singularity at i-th control point, and B j is the j-th singularity strength. The values of A ij depend on the singularity type.
In order to satisfy the boundary condition, a set of control points is defined on the lifting surface. Each control point provides one linear equation. Additionally, the Kutta condition at the trailing edge is applied to preserve the equality between the numbers of equations and unknowns (circulation of bound vortices).
Unlike most other lifting surface propeller calculations, this approach includes the presence of the hub, modelled by means of sources distributed over the hub surface in such a way as to satisfy the kinematic boundary condition on the hub [9]. Their intensities are regarded as additional unknowns when solving the linear equation system.
Once the strengths of the singularities are known, local total velocities can be determined and local pressures calculated. The ideal fluid forces are calculated through integration of the pressure distribution over the blade surface. Since the propeller blade is replaced with an infinitesimally thin vortex surface, there are no real back and face sides of the blade for pressure determination. The total velocity values for their pressure values can be determined in the following manner , where Γ i is the circulation of the nearest bound vortex and Δc is the distance between the neighbouring bound vortices.

VISCOUS EFFECTS
The no-slip boundary condition, which holds in real fluids, is not covered by hydrodynamic singularities. Hence the viscous effects have to be additionally modelled in an artificial way. In the present approach the blade profile viscous drag C D is calculated as [9]: , (8) where F D is the drag force acting on the whole blade section, ρ is the density of water, c is the blade profile chord length at considered radius, t is the maximum thickness of profile section at considered radius, and R e is the Reynolds number for the considered blade section.
The boundary layer growth over the propeller blade leads to the reduction in effective pitches and cambers of blade section profiles. A simple formula for pitch and camber reduction is applied [4]: , (9) where Δφ is the pitch angle reduction, expressed in radians, f is the maximum geometrical camber of blade section profile at considered radius. Finally, the maximum camber is reduced in the following way: , Importantly, as the pitch and camber corrections contain only geometrical quantities, they can be applied at the geometry discretisation stage.

FREE VORTEX SYSTEM
Lifting surface calculations should include the propeller vortex wake as part of the system of singularities. Most classical applications assume a true helical surface of constant pitch [7,4]. However, both the experiment [4], and viscous fluid calculations [13] show that free vortices deform and do not hold the true helical surface shape. Despite this discrepancy, the classical approach gives quite good results near the design point. However, in order to determine correct propeller characteristics for other conditions, a special procedure for free vortex sheet shape determination is to be adopted. In most solutions to this problem [4,9,12] the propeller free vortex sheets are discretised as sets of short straightline segments. Total velocities are calculated at their endpoints and then they are convected by dr = U dt to their new locations. The time step dt is determined upon the arbitrary maximum distance dm, prescribed by the user before calculations as dt = dm/U. Here, the calculations were conducted with dm = 0.02 D. The time step may vary between subsequent iterations.
If the velocity vector is expressed in the cartesian coordinate system, this leads to artificial radial expansion of free vortices. In order to avoid this problem, a cylindrical coordinate system is used instead. This simple approach allows to achieve reasonable results for high and moderate advance coefficients. The propeller thrust and torque values converge after 3-5 iterations.

DETERMINATION OF INDUCED ADVANCE ANGLE
The radial distribution of the induced advance angles is one of the essential quantities, usually determined during the propeller blade design. Its correctness is merely achieved indirectly by confirmation of the designed propeller's load values at the design point. In this paper a simple method to determine the induced advance angle is proposed.
The induced velocities are computed over the propeller blade during the lifting surface calculations. The mean integral value is defined as , (11) where U a is the mean axial induced velocity and u a (c) is the local value of the axial induced velocity at a particular chord station. A similar formula can be formulated for the tangential induced velocities. It is noteworthy that the mean value of the induced velocity refers to the velocity induced on the propeller disc at the corresponding radius which appears in the propeller lifting line model. If so, the induced advance angle can be determined in the following way , The induced advance angles can also be validated by the dynamic criterion. This means that the thrust and torque values determined by integrating pressure forces should be equal to those calculated upon the lifting line theory, namely where Γ is the accumulated circulation of bound vortices at the corresponding radius, and V w is the local inflow velocity, defined as , (14) The results obtained by means of the above method are encouraging. What is important, this method can be used for constructing more sophisticated methods for propeller design, especially those cooperating with upstream stators. coefficients (J < 0.3). This is because vortex model calculations are unable to converge. Further comparisons are discussed towards the end of this article (Section 5).

CFD TURBULENCE MODELLING
The RAS (Reynolds-Averaged Simulation) approach is utilised to model the averaged turbulent flow. A closed system of equations can be formulated for the incompressible fluid [20]. It consists of the averaged continuity equation , (16) and the steady-state vector Reynolds equation, provided that the Boussinesq assumption holds , (17) In the above equations, p e denotes the effective pressure p e = ρ −1 + 2/3 k, and ν e is the effective viscosity ν e = ν t + ν. The number and form of additional equations depend on the adopted turbulence model. For the two-equation k-ω SST model [14] we have two more equations. One of them is the steady-state equation for transport of kinetic energy of velocity fluctuations k, in the form The other is the steady-state equation for transport of turbulence frequency ω, namely , The first blending function F 1 is defined to be close to 1 near the wall and 0 far from it

PROPELLER DISCRETISATION AND GRID CONVERGENCE
The propeller blade is discretised both, in radial and chordwise directions. The geometry approximations of the radii are determined with sine spacing, namely: , (15) where i is the actual index, r 0 is the dimensionless propeller core radius, and A is the number of radial discretisation stations. Further, the propeller geometry parameters, such as chord length, maximum camber, etc., are interpolated among these radii by means of second order polynomials. As for chordwise discretisation, uniform spacing is adopted. This is a reasonable approach, as no cavitation calculations are performed. The face and back side ordinates of blade sections are interpolated along the radius in a similar manner as the geometrical characteristics. Figure 5 presents the results calculated for the propeller design point (advance coefficient J = 0.442) with the prescribed vortex wake geometry (true helical surface of mean pitch ratio H/D = 0.5581). Based on these results, it was decided that the 20×20 singularity grid preserves sufficient accuracy and optimal time of calculations. Figure 6 shows the results of lifting surface calculations. The label 'Standard' refers to the calculations performed for propeller blades alone, whereas the label 'Improved' corresponds to calculations for the blade-hub configuration. Significant influence of the hub presence on the results is clearly visible. The thrust and torque coefficients for the case with hub are closer to experimental results and so is the efficiency. This is possible because of coupling of the lifting surface approach for the propeller blades and the boundary element method for the hub. Moreover, the improved vortex method takes advantage of the proposed simple algorithm to determine the propeller induced advance angles. However, it is not possible to calculate points for lowest advance where ,

RESULTS
The distance from the wall is represented by the parameter y. The SST turbulence model redefines the eddy viscosity ν t in order to avoid over-prediction of shear stresses near the wall, namely , Finally, the second blending function F 2 , being similar to F 1 , is defined as .

EQUATION DISCRETISATION
All the considered equations can be regarded as the generic scalar ϕ transport equation. The steady-state version of this equation is given by , where Γ stands for diffusion coefficient and S ϕ represents sources or sinks of the transported variable ϕ. The left-hand side of Equation (24) represents convection, while the first term on the right-hand side describes diffusion.
The corresponding discretised version of the Finite Volume Method [2] for the transport equation (24) is: , The control volume V around the centroid P can be of any shape, provided that it is convex and consists of f planar faces. The face area vector S f points outward from the control volume V. Most importantly, the cell face variables ϕ f need to be interpolated from the centroid values of the control volumes.
The divergence schemes include both convection terms and other diffusive terms. All of them take advantage of Gauss integration. The variable ϕ is interpolated to the cell face variables ϕ f by the second-order limited scheme. This make it possible to limit towards first-order upwind in regions of rapidly changing gradients. Additionally, specialised versions of these schemes are utilised for the velocity to calculate a single limiter applied to all velocity components. In order to maintain boundedness of the steady-state solution, bounded variants of convective schemes are considered.
The surface normal gradient term ( ϕ) f · S f , which is present in the diffusive terms, is evaluated at the cell face that connects two cells. In order to maintain the second-order accuracy for non-orthogonal meshes, apart from the orthogonal scheme, a non-orthogonal correction is considered. Moreover, the surface normal gradient is needed to calculate Laplacian terms using Gaussian integration.
The gradient terms, such as p, are discretised by means of Gaussian integration. The interpolation of values from cell centres to face centres is achieved by central differencing. A cell limiting idea is introduced here in order to improve the boundedness and stability. This idea makes it possible to limit the face gradient within the bounds of gradients calculated for adjacent cell centres.
Finally, linear interpolation of values from cell centres to face centres is used, which is a standard approach.

Mesh
The flow domain shown in Fig. 7 is divided into two parts, namely the rotating propeller and the steady ambient volume. Both are discretised separately and merged by means of arbitrary mesh interfaces. Computational grids consist of mostly hexahedral elements. The mesh of the propeller can be inspected in Fig. 8.  Table 1 shows the detailed mesh statistics for the considered flow geometry. In order to ensure that the flow near the walls is properly resolved, thin layers around the blades, hub and shaft are generated. The quality of the mesh near the blades is inspected in terms of the average and maximal values of y + distribution. The former value being , whereas the latter is defined as the maximum value of y + on the blade , In the above formulas, |S| stands for the area of blade surface S. The average ȳ + value on the blades is 0.16 and the maximal value is 3 for all advance coefficients J. Figure 9 demonstrates mesh convergence by showing the influence of the number of nodes on the efficiency η. It is obvious that increasing the number of mesh nodes above 10 × 10 6 has negligible effects on the efficiency η. This is because the results are nearly constant above 10 × 10 6 nodes. The mesh convergence test is shown for the advance coefficient J = 0.59.

BOUNDARY CONDITIONS
The main boundary conditions are the following: • Inlet. The specified constant velocity Ū according to || Ū || = V A = JnD is directed perpendicularly to the inlet surface, together with the zero normal pressure gradient ∂p/∂n = 0. Because of the two-equation k − ω SST turbulence model, two additional variables are necessary, namely k and ω. The low turbulence intensity case is assumed, meaning that the turbulence intensity τ t = 1% and the viscosity ratio ν t / v = 1 are considered. • Outlet. The constant pressure distribution is assumed here, along with the zero normal velocity gradient. This is justified simply because the outlet surface is located far from the propeller. • Wall. The impermeability and no-adhesion conditions are specified, i.e. the slip condition. This is true for external walls in order to minimise their influence on the propeller. As for the propeller blades, hub and shaft, the no-slip condition is applied, i.e. the impermeability and adhesion requirements are forced. This, however, is true in the rotating frame of reference. The flow in the near wall region is modelled by means of the scalable wall function to provide near wall boundary conditions for the mean flow. • Interfaces. An arbitrary cyclic mesh interface is considered to allow for coupling between the stationary outer domain and the rotating inner domain of the propeller. In order to save the computational time, compared to the transient rotorstator interaction, a steady state approach is used. This is commonly referred to as multiple reference frame simulation.

SOLUTION METHOD
The governing equations (16) -(23) are solved by means of the open source CFD software, namely OpenFOAM [17,18]. More specifically, the SIMPLE algorithm (semiimplicit method for pressure-linked equations) is used to solve pressure-velocity coupling. The pressure equation is solved by means of the GAMG solver (generalised geometricalgebraic multi-grid) with DIC (Diagonal incomplete-Cholesky) smoother. For the velocity fields, as well as for turbulent quantities, such as k and ω, a standard solver using a GS (Gauss-Seidel) smoother is utilised.
The under-relaxation factors are used to improve the stability of the solution. This is important when solving steady-state flows. The variable ϕ n+1 in the next P iteration is limited by means of the However, the lower the under-relaxation factors, the slower the convergence. The here assumed under-relaxation factors are listed in Table 2. Figure 10 presents the efficiency η according to Equation (3), thrust coefficient K T (Equation (1)), and torque coefficient K Q (Equation (2)) as functions of the advance coefficient J (Equation (4)). Dots correspond to the numerical prediction and solid lines to the experimental data. Very good agreement is visible for the whole range of the considered advance coefficients J. Further comparisons with the vortex method are discussed towards the end of this article (Section 5).

CONCLUSIONS
Several numerical analyses of the CP745 marine propeller model have been conducted and compared with measurements carried out in the Ship Design and Research Centre towing tank. The experimental characteristics ('Experiment') shown in Fig. 11 were used as the reference for comparing the standard vortex method ('Standard'), the improved vortex method ('Improved') and the turbulent computational fluid dynamics approach ('CFD') by means of the Reynolds-Averaged simulation.
It is clearly visible that the standard vortex method for blades alone gives the worst results in comparison with the experimental data. The situation is, however, significantly better for the improved vortex method, owing to the coupling of the lifting surface approach for the propeller blades and the boundary element method for the hub. This is true for the torque and thrust coefficients as well as for the efficiency. Furthermore, the improved vortex method utilises the established simple algorithm to determine the propeller induced advance angles. The proposed modification provides better agreement with the experiment, compared to the original version of the vortex method. It is noteworthy that the accuracy of the improved method becomes comparable to CFD predictions. The CFD calculations are the most accurate. What is more, the CFD approach allows to calculate points even for lowest advance coefficients J < 0.3, for which the vortex model calculations are unable to converge. However, the CFD method is a few hundred times slower than the vortex methods. All CFD calculations were performed on the i7 930 2.66 GHz processor (3 of 4 cores engaged). The calculation for a single advance coefficient takes 18-19 hours. At the same time the vortex method requires only a few minutes on i5 4590 3.30 GHz (1 of 4 cores engaged). This makes the vortex methods still attractive, keeping in mind that there is still room for improvement.