AWakeModel for the Prediction of Propeller Performance at Low Advance Ratios

A low order panel method is used to predict the performance of propellers. A wake alignment model based on a pseudounsteady scheme is proposed and implemented. The results from this full wake alignment (FWA) model are correlated with available experimental data, and results from RANS for some propellers at design and low advance ratios. Significant improvements have been found in the predicted integrated forces and pressure distributions.


Introduction
For the prediction of propeller performance, the panel method shows great advantage over other numerical tools.On one hand, the panel method recovers the real blade geometry as much as possible, instead of simplifying the propeller blades into lifting surfaces as the Vortex lattice method (VLM) does.On the other hand, the panel method is far more efficient than Reynolds-averaged Navier-Stokes (RANS) simulation, not only in a computational sense, but also on grid generation.Generally, the panel method strikes a good balance on accuracy and computational efficiency.Close to the design condition, the results from panel method are in good agreement with either experimental measurement or RANS simulation.
As the advance ratio becomes smaller, the pitch of the wake sheets becomes smaller.Consequently, the wake shed from one blade could be close to the consequent blade and, thus, strongly affect the pressure distribution on the consequent blade.Therefore, in order to correctly predict the propeller performance at low advance ratios, an accurate wake model is necessary.Intensive investigation has been taken in the past two decades on modeling the wake of a propeller.
Greeley and Kerwin [1] invented a fast wake alignment scheme (PSF-2 type alignment).The wake is aligned on a coarse grid first, and then interpolated onto a fine grid.Pyo and Kinnas [2] presented a fully aligned wake model with biquadratic dipole panels.The high order wake model, although gained success on predicting the rollup of the wake and also the propeller performance near design condition, was more computational intensive than those with low-order panels.Later, Lee and Kinnas [3] developed a 3D wake model based on low-order panel method for marine propellers, which emphasizes more on the modeling of a developed tip vortex cavity.The tip vortex cavity is modeled with a helical vortex tube, of which the complete geometry is determined from requiring the pressure on its surface to be equal to vapor pressure.For the fully wetted cases, the tip vortex core is still necessary.Lee and Kinnas [4] also applied the same wake model to fully unsteady nonaxisymmetric flows.He [5] proposed a fully aligned wake model based on a VLM code, MPUF3A.In this wake model, the tip vortex core is no longer necessary.
Independently, Politis [6] formulated the unsteady motion of propellers with an unsteady wake model in an inertial frame of reference.Propellers under different, unsteady operating conditions (impulsive start, inclined inflow, and heave motion) were simulated.Greco et al. [7] proposed a sophisticated steady wake model, including the fully aligned near wake, helical intermediate wake and far wake disk, for analyzing the steady performance of propellers.This wake model was thoroughly validated with experimental data for propeller near design conditions, but the off-design performance was not addressed.
In this paper, the same idea as that in He's [5] wake model is improved within the framework of a low-order panel method code, PROPCAV.The panel code without any aligned wake model is first applied to a 3D wing at low angle of attack (AOA) on the purpose of verification.
With the new wake model and also the PSF-2 type wake model, the panel code is then applied to a 5-bladed propeller DTNRC 4381, especially at low advance ratios, and a 3bladed propeller DTMB 4119.In the cases of a 3D wing at low AOA and propeller DTNRC 4381 at different advance ratios, the results from PROPCAV are correlated with RANS.In the cases of propeller, comparisons with available experimental measurements are also carried out.

Formulation
The velocity flow field is decomposed into two components: inflow velocity and perturbation velocity, written as follows: where q is the total velocity, U ∞ is the inflow velocity, and u is the perturbation velocity due to the presence of obstacles.

Governing Equations.
The methodology applied in this study is based on the inviscid potential flow theory: where φ is the perturbation potential.The perturbation velocity field, therefore, is governed by the Laplace equation: If the inflow is also irrotational, we can also define the total potential Φ as follows: Appling Green's second identity, the Laplace equation ( 3) can be written as the following boundary integrated form: where S H represents the surface of hydrofoil or propeller blades and S W represents the surface of trailing edge wake.G(p, q) is the Green's function, which is defined as 1/R(p; q) in 3D, and R is the distance between the two points p and q.A low-order panel method using constant dipole panels is adopted to solve (5).

Kinematic Boundary Condition.
Applying the kinematic boundary condition to the solid surfaces, we obtain the following: 2.2.2.Pressure Kutta Condition.At trailing edge, Δφ w is unknown, which is subjected to the pressure Kutta condition: where the + sign denotes the suction side panel, and thesign denotes the pressure side panel.The physical meaning of ( 7) is that free vortex sheets cannot bear pressure jump.
For steady state, (7) leads to the modified Morino-Kutta condition [10]: where r TE is the vector from the control point of (−) panel to the (+) panel at the trailing edge.
For unsteady problem in 2D, or even steady problem in 3D, (7) has to be enforced through the unsteady Bernoulli equation: An iterative scheme has to be applied in order to enforce (9), proposed by Kinnas and Hsin [11], namely the Iterative pressure Kutta condition (IPK).

Evaluation of Local Velocity.
The local perturbation velocities in the computational domain are not required for evaluating the pressure on the blade surface, but evoked by the wake model.Theoretically, the local perturbation velocity u at any point x p can be evaluated through the following singular integration: where G(p, q) is the same as Green's function that is defined in (5).Practically, the numerical evaluation of (10) casts two questions: where and how to evaluate the velocity.If the field point p is located right on the vortex sheet, (10) contains singular integration.For low-order panel method, evaluating the velocity at the control point of a panel can avoid the singularity.Lee and Kinnas [3] implemented this idea in their unsteady wake model.In order to align the wake, the velocities at the nodal points of a panel are interpolated from the velocities at the control points.Repaneling was carried out at each time step.
However, the velocities are finite at the control points only if the panels are regularly and smoothly distributed in space.At the late stage of the unsteady free wake evolution, or in certain intermediate steps to determine the steady wake, the geometry of the vortex sheet could be entangled.The velocity evaluated at the control point of a panel could lead to a huge deviation from the exact value.Numerical tests show that evaluating the velocity at nodal point with certain desingularized vortex kernel can be more robust than the schemes mentioned above.
In 2D, the dipole-induced velocity in ( 10) is equivalent to the induced velocity due to a vorticity distribution with strength φ/∂s: where R ⊥ is the vector which is perpendicular to (11), which led to the vortex blob method in the study of the evolution of 2D free vortex sheet.This method, although simple, was successful in many 2D applications.The vortex-induced velocity is limited to a finite value, and the inherent Kelvin-Helmholtz instability is suppressed.
In 3D, the vortex sheet is discretized into constant dipole elements, which are equivalent to vortex loops with constant vorticity strength.Therefore, induced velocity due to the vortex loop can be evaluated using the Biot-Savart law: Replacing the Biot-Savart kernel R/R 3 , with the Rosenhead-Moore kernel R/(R 2 + δ 2 ) 3/2 , leads to the 3D extension of the vortex blob method.Studies by Ramsey [13], He [5], and Lindsay and Krasny [14] show that this method can give out reasonable results for the 3D evolution of vortex sheets.

Wake Alignment Scheme.
The free vortex sheets shed from the trailing edge from a propeller blade are material surfaces, which have to be aligned with the local flow velocity.Alignment schemes play very important roles in determining the spatial geometry or distribution of the free vortex sheets.For the steady and unsteady cases, the alignment schemes are numerically different.
Engineering community is more interested in the steady propeller performances.As an analogy to the Euler-explicit scheme for the unsteady cases, many researchers align the vortex sheet using the velocity at the upstream nodal points: Δx, Numerical tests show that this scheme is stable and converges fast.However, this scheme has the same deficiency as the Euler-explicit scheme for the unsteady alignment: it continuously enlarges the rotating radius around the center of rollingup, and thus cannot predict the correct location and size of the rolling-up region.Therefore, in order to get better prediction of the wake geometry, alignment scheme in a manner of trapezoidal rule is preferred.However, direct application of the trapezoidal rule turned out to be unstable.
Tian and Kinnas [15] proposed a pseudounsteady alignment approach when simulating the leading edge vortex (LEV) for delta wings with sharp leading edge.In the LEV model, the accuracy of the size of the rollup region is more important than that in the free wake model.In the study of this paper, although no LEV model is applied for predicting the propeller performance, the same idea is adopted for aligning the propeller wake.
As shown in Figure 1, consider a material line defined as follows: We have which leads to Discretizing ( 16) with finite difference scheme for the slopes and with the Euler-Explicit scheme for the unsteady terms, we have where q y,i = 1/2(q n y,i−1 + q n y,i ), q x,i = 1/2(q n x,i−1 + q n x,i ), and i , and β = Δt/Δt * i , we have Notice that ( 18) is the alignment scheme for vortex sheet under horizontal uniform inflow.In other words, the inflow is parallel to the x-axis.The idea behind this approach is that the fixed coordinate has to be in the same direction of the dominant inflow, therefore, (18) behaves in an upwind manner numerically.For example, the slope of the curve at current point is evaluated using current point and the upstream point.However, in the case of the flow around a propeller, the dominant inflow would not be in the axial direction, but along helices with constant pitch.The rotational component of the inflow is important.The flow around the propeller is solved relative to a coordinate system which rotates with the propeller.In other words, the U ∞ in (1) should be replaced by U ∞ − ω × r, where ω is the angular velocity of the propeller, and r denotes the radial vector from the rotation axis to any points on the blades of the propeller.The essence of the alignment scheme described above is to decompose the total velocity vector into two components: (1) component parallel to the inflow and (2) component normal to the inflow; then keep the length of the vortex segment in the inflow direction unchanged.The advance ratio J S is defined as J S = U ∞ /nD, where n and D are the rotational velocity in revolutions per second, and diameter of the propeller respectively.
Consider a point under cylindrical coordinates x i−1 , r i−1 , θ i−1 , with Δθ being the grid size parameter in the streamwise direction.If the wake is only aligned with the inflow, the consequent point will be at x i = x i−1 + RJ s Δθ/π; r i = r i−1 and θ i = θ i−1 + Δθ.Therefore, the effective inflow direction in the Cartesian system on the segment connecting the (i − 1)th and the ith points is where in Cartesian system Δs i is defined as Now the averaged perturbation velocity u i = 1/2(u i + u i+1 ) can be decomposed into two components: where u i,s is the scalar projection of u i in Cartesian system onto s i , and u i,n is the vector component of u i normal to s i .
Defining Δt * i = |Δs i |/( u i,s Δθ/ω prop + |Δs i |), and β = Δt/Δt * i we obtaine the alignment scheme for the wake of propeller: where x n i denotes the coordinates of the ith point at nth time step.It is easy to show that if the projection of the vector 22) is implemented in the full wake alignment (FWA) model.Although (22) can be used for unsteady problems, in the study of this paper, only steady alignment is addressed.Also, the symmetry or periodicity of the flow field is exploited.
When applying (22) for the steady alignment, the wake shape is recalculated through (22) repeatedly without resolving for the potentials until reaching the maximum number of iterations, or reaching a given convergence criterion.Then, the influence coefficients, due to the wake are reevaluated.With the updated influence coefficients due to the wake, the potentials on the blades and also on the wake are resolved for.Repeating this procedure will lead to converged wake shape and converged forces on the blades.

Application to 3D Wing at Low Angle of Attack with Straight Wake
Before applying the FWA model, the basic framework of current panel code is first tested with a 3D wing at low angle of attack (AOA).A tip-swept wing, which is expected to resemble a propeller blade with a large skew angle, is investigated.As shown in Figure 2(a), sweep is superimposed onto an elliptic planform, which has an aspect ratio AR = 1.57.All the sections of the wing have NACA 0010 thickness form without camber.A straight wake without alignment is attached at the trailing edge of the wing.In this case, no LEV model is applied.The results at low AOA (AOA=2DEG) from current BEM model is validated with RANS, which is obtained through commercial software STAR-CCM+.
k − ω SST turbulence model is adopted for the RANS simulation, and the wall y+ is mainly controlled in the range between 30 and 80.For the convenience of grid generation, the trailing edge of the wing is made to be round, which has a radius of 0.5% of the chord length at each section.The RANS case took around 3 hours with 16 Intel Xeon 2.54 GHz CPUs to converge.The BEM code took 3 minutes on a single CPU of the same type to get the results.
Three sections are taken to compare the pressure distribution.As shown in Figures 2(b)-2(d), the pressure distributions on all the sections are in good agreement with both methods.

Application to Propeller DTNSR 4381
The panel code is applied to a five-bladed propeller DTNSR 4381, which is investigated by Boswell [8] through experimental measurement.The design advance ratio J s of this propeller is 0.889.
For the simulations with the panel code, two wake models are adopted: PSF-2 wake model, which is a fast wake alignment scheme developed by Greeley and Kerwin, [1] and a full wake alignment (FWA) scheme, which is implemented in this study.For all the advance ratios, wakes of one revolution are modeled with the FWA scheme.For the PSF-2 model, wakes are trimmed at 2.5 R downstream with an ultimate wake disk.The hub of this propeller is also modeled through BEM panels, but the hub vortex core is omitted in this study.
RANS simulations of the same propeller are carried out with commercial code FLUENT, for J s = 0.883, 0.6, 0.5 and 0.4.K-ω SST turbulence model is adopted.QUICK scheme is used for spatial discretization, and SIMPLEC scheme is applied for pressure correction.1.5 million hexahedral cells are used to simulate one sector of the domain with periodic boundary condition.It took 12 to 14 hours on 32 Intel Xeon 2.54 GHz CPUs for the residuals to converge to 1E-6.
In the experimental study by Boswell [8], only integrated forces were measured.Figure 3 compares the integrated forces predicted by PROPCAV and RANS, and the experimental measurements.We can see that that the predicted KT from the panel code with either wake model agree the well with both experimental data and RANS near design J s .However, as the advance ratio becomes smaller, the KT predicted by the PSF-2 alignment starts deviating from the International Journal of Rotating Machinery  Figure 4: KT and 10KQ for propeller DTNSR 4381 from PROPCAV with different number of spanwise panels by chordwise panels on the blade and different streamwise grid size on the wake, where 6DEG stands for 60 streamwise wake panels per turn, and 10DEG stands for 36 streamwise wake panels per turn.
experimental data and the RANS results.At the meantime, the KT from fully aligned wake model is still in good agreement with the experimental data.The panel code with either wake model tends to slightly overestimate the KQ near the design condition, but the results from PSF-2 alignment deteriorate at low advance ratios.Figure 4 shows the results of convergence analysis of the BEM/FWA model: different spatial discretization on the blade and on the wake leads to almost the same KT and KQ. Figure 5 shows the aligned wake geometry at two different advance ratios (J s = 0.889 and J s = 0.5).Clearly at low advance ratio, the wake sheets are strongly rolledup.Figure 6 plots the wake geometry predicted by PSF-2 alignment at the same advance ratios as in Figure 5. Obviously, PSF-2 alignment did not predict the contraction of the transition wake radius.The contraction rate of the transition wake is a user input parameter in PSF-2 alignment.We simply set it as zero because this parameter differs from case to case.On the other hand, Figure 3 indicates that near the design condition, whether contraction of the wake radius has been predicted (FWA) or not (PSF-2), the predicted total performance of the propeller will not be affected significantly.Figure 7 shows sections of the aligned wake cut at the z = 0 plane, from different J s with different spatial discretization on the blade, where propeller frame of reference is defined as: x is along the axial direction, y is along the radial direction passing the middle of the root of the key blade, and z is normal to x and y directions.The rollup of the tip vortex can be observed.Also, the wake geometry with different spatial discretizations converged to similar results.Notice that in the current wake model, no tip vortex core model is adopted.More detailed comparisons of the results between BEM and RANS are carried out in two different blade sections at two different advance ratios, as shown in Figure 8.The correlations among different numerical methods are remarkably good.Results from panel method with fully aligned wake model are closer to that of RANS, comparing with the PSF-2 wake model.This is consistent with the results for the total forces.Also for the FWA model, different spatial grids yield almost coincident results.The pressure distribution from 30 × 60 and 30 × 80 panels is not shown in Figure 8 so that the difference of the pressure from RANS and BEM/FWA can be clarified.
In the case of 30 × 80 panels on the blade, the computational time for PROPCAV with PSF-2 wake model was only 1 minute on a single Intel Xeon 2.54 GHz CPU.The FWA model for the same discretization on the blade with 60 streamwise wake panels took 30 minutes to converge on the same type of CPU.Comparing with RANS, it is a big advantage to be able to achieve the results with the same order of accuracy in such a short time.It is also easy to parallelize the FWA model.For the wake model itself, linear parallel speedup could be achieved.
Figure 9 lists the convergence history of the FWA model for the Propeller 4381 at J s = 0.5.The averaged moving distance of the points at the last streamwise station on the wake is monitored.Notice that once the wake geometry is updated, the boundary integral equation has to be resolved for the corresponding dipole strength on the wake.However, resolving for the unknowns is not carried out at each step, as mentioned before, inversion of the coefficient matrix is only evoked at selected steps.
It should be pointed out that when the advance ratio is low enough, the leading edge vortex (LEV) could be developed, and affect the pressure distribution and total integrated forces.Greeley [16] proposed a semiempirical LEV model for vortex lattice method.Tian and Kinnas [15] started to predict the LEV effects in a more rational way.At the same time, the boundary layer effects on the pressure distribution on the blade would be noticeable at International Journal of Rotating Machinery  high loadings.In order to obtain the correct KT and KQ, the boundary layer effects have to be taken into account correctly as well.Sun and Kinnas [17] developed a viscous/inviscid interaction approach by coupling the panel code with a boundary layer solver, XFOIL.
The panel code in this study neither included the LEV model nor coupled with the boundary layer solver.The viscous effects are considered by using an empirical viscous pitch correction [18] and a skin friction coefficient, which corresponds to that of a flat plate at the same Reynolds number based on the chord length.It is expected that incorporating the LEV model and the viscous/inviscid interactive scheme into current code will make current results even better.

Application to Propeller DTMB 4119
Jessup [9] carried out experimental measurement of the wake flow of a three-bladed propeller, DTMB 4119.The BEM code with the FWA model is applied to this propeller, in order to see the effects of the desingularized parameter δ in the Rosenhead-Moore kernel R/(R 2 + δ 2 ) 3/2 when evaluating the dipole induced velocities.Three different values of δ are selected (0.03; 0.05; 0.075).
Figure 10 compares the predicted radial location of tip vortex for propeller DTMB 4119 at design advance ratio with the experimental measurement by Jessup [9].All the results are close to the experimental data.It seems that smaller values of δ lead to better prediction of the location a 3D wing at low angle of attack.The predicted pressure distribution agrees well with the RANS results.In case of the propeller working at design condition, both PSF-2 and FWA wake models predict reasonable KT and KQ.The results from FWA model are still in good agreement with the experimental data even as the advance ratio becomes smaller.Strong rollup of the wake at low J s is predicted by the FWA model.Detailed pressure comparisons with RANS show that FWA significantly improves the simulated pressure distribution.
The FWA wake model is also partially validated with the experimental measurement of the wake flow of propeller DTMB 4119 by Jessup [9].Correct radial location of tip vortex is predicted.Sensitivity study indicates that different values of the desingularized parameter δ only have limited effects on the predicted total performance of propellers.
On the evaluation of KT and KQ, the viscous effects are simplified with an empirical viscous pitch correction [18] and a skin friction coefficient.In the near future, a viscous/inviscid interaction model will be included into the current model in order to treat the viscous effects in a more rational way.A LEV model will also be added into the current wake model, so that more accurate pressure distribution towards the tip of the propeller blade can be obtained.
Correlation of the FWA model with the latest experimental data can also be very helpful to improve the methods described in this paper.Typically the PIV measurement of the wake flow of the INSEAN E779A propeller by Di Felice et al. [19] reveals many detailed vortical structures.An investigation will be carried out if the geometry of this propeller is made available.

Figure 1 :
Figure 1: Schematic plot for the pseudounsteady approach of alignment.

Figure 2 :
Figure 2: Geometry of the wing and comparison of pressure distributions with RANS.

Figure 7 :
Figure 7: Sections on the aligned wake cut at z = 0 plane at different J s (0.889, 0.7, 0.5 and 0.4) with different discretization on the blade.

Figure 10 :δFigure 11 :
Figure 10: Radial location of tip vortex for propeller DTMB 4119 at design J s (0.833), from different values of δ, the triangles represent the experimental measurement by Jessup [9].