Actuator disc methods for open propellers: assessments of numerical methods

The paper describes the assessment of two different actuator disc models as applied to the flow around open propellers. The first method is based on a semi-analytical approach returning the solution for the nonlinear differential equation governing the axisymmetric, steady, inviscid and incompressible flow around an actuator disc. Despite its low computational cost, the method does not require simplifying assumptions regarding the shape of the slipstream, e.g. the wake contrac- tion is not disregarded or prescribed in advance. Moreover, the presence of a tangential velocity in the wake as well as the spanwise variation of the load are taken into account. The second one is a commonly used procedure based on CFD techniques in which the effects of the propeller are syn-thetically described through a set of body forces distributed over the disc surface. Both methods avoid the difficulties and the computational costs associated with the resolution of the propeller blades geometrical details. The comparison is based on an in-depth error analysis of the two pro- cedures which results in a set of reference data with controlled accuracy. An excellent agreement has been documented between the two methods while the computational complexity is obviously very different. Among other things the comparison is also aimed at verifying the accuracy of the semi-analytical approach at each point of the computational domain and at quantifying the effect of the errors embodied in the two methods on the quality of the solution, both in terms of global and localperformanceparameters.Furthermore,thepaperprovidesasetofreferencesolutionswithcon-trolledaccuracythatcouldbeusedfortheverificationofnewandexistingcomputationalmethods.Finally,thecomputationalcostofthesemi-analyticalmodelisquantified,thusprovidingavaluableinformationtodesignerswhoneedtoselectacosteffectiveandreliableanalysistool.


Introduction
The 'actuator disc method' (ADM) constitutes a widely employed design and/or analysis tool both in its analytical and CFD-based formulation. Besides the most famous and simple momentum theory (Glauert, 1935), a nonlinear variant of the actuator disc model has also been developed by Wu (1962) in his pioneering work. Since then Greenberg and Powers (1970) and Greenberg (1972) provided important improvements to the nonlinear ADM both in its theoretical and numerical aspects, and a detailed review of the most relevant analytical matters can be found in the book of Breslin and Andersen (1994).
More recently, Conway (1998) provided an exact and implicit solution to the flow around an actuator disc. The solution can be regarded as the flow induced by a distribution of ring vortices modelling the propeller wake. The solution is made explicit through an iterative and semianalytical procedure, and, for this reason, the method will be termed the 'semi-analytical actuator disc method' CONTACT R. Bontempo rodolfo.bontempo@unina.it (SA-ADM) from here on. The method deals with rotors characterized by heavy loads with non-uniform distributions and with slipstream rotation and contraction. Forgber, Lindner, & Schwarze, 2015;Hu, Li, Han, Cai, & Xu, 2016;Insinna, Griffini, Salvadori, & Martelli, 2014;Liu et al., 2015;Luo, Yu, Dai, Fang, and Fan, 2016;Manna, Benocci, & Simons, 2005;Manna, Vacca, & Deville, 2004;Montis et al., 2009;Shi, Xu, & Wei, 2016;Sun, Su, Wang, & Hu, 2016). CFD techniques are also widely employed to solve the Reynolds averaged Navier-Stokes (RANS) or Euler equations in a domain comprising a propeller that is represented through a body-force distribution over a disc area. For this reason we name this approach the 'CFD actuator disc method' (CFD-ADM). The main goal of this technique is to reduce the computational efforts related to the solution of the flow around the propeller blades. Moreover, since the screw is replaced by a simple disc, mesh generation efforts are greatly reduced. Obviously, such a method cannot predict the details of the flow around the rotor and the blades tip, but it can successfully address the global effects of the propeller on the hull flow, as witnessed by the short literature review reported hereafter.
At the end of the 1970s, Favin (1977, 1979) attempted to represent the propeller through an axisymmetric distribution of body forces weakly coupled to the incompressible RANS equations. Kawamura, Miyata, and Mashimo (1997) simulated the flow around five tanker hulls introducing the Nakatake (1989) simplified actuator disc propeller model into the WISDAM-V finite-volume method developed by Miyata, Zhu, and Watanabe (1992). In the above procedure, the interaction between the two flow models is handled through a time-marching procedure which should converge towards the steady self-propelling condition. Phillips, Turnock, and Furlong (2009) developed a coupled blade element-RANS procedure to determine the manoeuvring coefficients of a self-propelled ship hull travelling straight ahead. The analysis was carried out at a prescribed drift angle and for differing rudder angles. In their approach, the blade element method is used to model the effects of the propeller with a set of forces which are then inserted into the RANS domain. The latter may include the ship hull and the rudder. A circumferentially averaged nominal wake fraction is evaluated with a RANS code at different radial stations. The computed wake velocity is used to evaluate the propeller loads through the blade element method. Then, the calculated thrust and torque values are introduced in the CFD code as momentum sources distributed over a disc of finite thickness. This procedure is repeated until convergence. Choi, Min, Kim, Lee, and Seo (2010) examined the speed-power performance of various types of commercial ship hulls analysing their resistance and propulsion characteristics. In their study asymmetric body forces are used, while the effect of a finite number of blades is neglected. The adopted procedure starts by first evaluating the incoming flow velocity at the propeller plane by the commercial RANS solver ANSYS Fluent . Then, the obtained inflow is used as input data for a potential-flow solver. After that the thrust and torque distributions on the actuator disc plane are represented as known body forces in the RANS code by means of a 'user-defined function' (UDF). Obviously, the procedure is iterated until convergence. The UDF is a technique that allows any user to implement certain models into the general multi-purpose ANSYS Fluent code. In that context the UDF was used to insert the asymmetric bodyforce propeller model into the CFD code. Starke and Bosschers (2012) adopted an hybrid boundary element method-RANS approach to investigate the scale effects in ship powering performance. The viscous flow around the ship hull is coupled with the flow induced by the screw modelled with a boundary element method. Moreover, an extensive description of the coupling procedure details is also reported. In order to predict the free turning manoeuvre of a tanker-like ship, Broglia, Dubbioso, Durante, and Di Mascio (2013) coupled different propeller actuator disc models with an in-house developed URANS code. Two models were considered: a modified Hough and Ordway (1964) approach and a model based on blade element momentum theory.
The literature review presented above proves that the ADM is currently and successfully employed to study propeller wakes. This paper, which relies on the work of Conway (1998), quantifies the effect of the errors embodied in the SA-ADM and CFD-ADM on the quality of the solution, both in terms of global and local performance parameters. A code-to-code comparison has been carried out with the objective of establishing the computational cost required by the two methods to comply with a prescribed accuracy in the clean context offered by the selected test cases. In addition to what has just been described, further original contributions of the paper can be found. For example, the comparison is also aimed at verifying the accuracy of the SA-ADM at each point of the computational domain, not only in the far wake as previously reported in Conway (1998). Furthermore, the paper provides a set of reference solutions with controlled accuracy that could be used for the verification of new and existing computational methods. Finally, the computational cost associated with the SA-ADM method is quantified; a fact of the utmost importance for designers who have to select the optimal analysis tool to be integrated in a design system. The paper is structured as follows. Section 2 comprises a synthesis of the semi-analytical actuator disc of Conway (1998), aimed at pointing out the most important numerical issues discussed in the forthcoming error analysis. Then, the SA-ADM and CFD-ADM results are compared in Section 3. With this aim, the popular CFD suite ANSYS Fluent is adopted.

The SA-ADM
Consider a cylindrical coordinate system (ζ , σ , θ). The actuator disc centre is located at (0,0,0), its radius is σ ad and the velocity of the free stream is U ∞ (see Figure 1). Suppose the flow to be axisymmetric, steady, incompressible and inviscid. Since the flow is axisymmetric, the Stokes stream function , defined as can be conveniently introduced. In the above equation, the velocity vector is u = (u, v, w). Moreover, for the class of axisymmetric problems the vorticity vector reads Furthermore, by means of Equation (1), the θ-component of the vorticity can be written as Note that the minus sign in the first-order derivative term makes the above equation different from the well-known axisymmetric formulation of the Poisson equation in cylindrical coordinates. Using the momentum and the angular momentum equations (Wu, 1962), the above Equation (3) returns This is the partial and nonlinear differential equation that governs the through-flow across the device. In particular, is the angular velocity, while the load distribution H is defined as In the above equation, H = p/ρ + (u 2 + v 2 + w 2 )/2 is the Bernoulli constant and σ s (ζ ) is the slipstream location defining the radial extension of the wake as a function of the streamwise coordinate ζ . Thus, the wake is mathematically defined as the space region The conditions at infinity that have to be associated with Equation (4) are (Wu, 1962): The Stokes stream function definition (1) can be easily employed to clarify the meaning of the previous conditions. From Equation (4) it is easily understood that, to formulate the differential problem (4)-(5) completely, two key physical quantities, i.e. the angular velocity (or the advance coefficient J ) and the distribution of the load, have to be known in advance. Indeed the function σ s (ζ ), defining the wake shape, is not known beforehand and it has to be computed through the iterative method detailed in Section 2.2. Eventually, as soon as the (u, v) components of the velocity field have been obtained, the tangential velocity component w can be evaluated by the well-known angular momentum equation

Exact solution of the flow around an actuator disc
As thoroughly described in Conway (1998) and the references therein, the flow around an actuator disc is modelled with the help of ring vortices (see Figure 2). In more detail, the stream function and the velocity field (u , v ) associated with a single ring vortex of radius r, strength κ and located at ζ = z and σ = 0, are (Basset, 1888;Lamb, 1932) where for ζ − z ≥ 0 the positive sign holds and vice versa. The quantity J, appearing in the above equation, represents a Bessel function of the first kind. The solution for can be written as the superposition of the flow induced by a distribution of ring vortices modelling the propeller wake. The density strength of this distribution is named γ ad (ζ , σ ), where the subscript 'ad' stands for 'actuator disc', meaning that the ring vortices are used to describe the flow induced by the disc. The exact flow solution is then obtained by integrating Equation (7) over the wake, i.e. (Conway, 1998) Specifically, the first term in Equation (10) represents the flow promoted by the actuator disc, whilst the second term is the flow induced by the free stream.
Consider now the requirements that must be met by γ ad (ζ , σ ) so that Equation (10) can be regarded as the solution of the differential problem (4)-(5). As previously stated, the distribution of the load and the advance coefficient have to be preliminarily prescribed. Moreover, these two physical quantities are directly related to the tangential vorticity distribution in the wake through Equation (4). This means that the actuator disc ring vortex distribution, i.e. the density strength γ ad (ζ , σ ) appearing in Equation (10), has to reproduce the prescribed wake vorticity field ω θ (ζ , σ ) induced by the disc H( ) distribution (see Equation 4). Note that, for the through-flow field (u, v, 0), ω θ is the only vorticity component that is different from zero, and the vorticity strength vector of a ring vortex is also directed in the θ -direction (see Figure 2). From these considerations it can be easily inferred that the following equation must hold: γ ad (ζ , σ ) = ω θ (ζ , σ ) (see Conway, 1998). The conditions at infinity (5a) and (5b) are obviously satisfied by the overall flow solution (10). Finally, integrating Equations (8) and (9) over the wake space region, the disc induced velocities u ad and v ad read

SA-ADM: solution strategy
As stated in the previous subsection, Equation (10) is the exact solution of the through-flow around an actuator disc. However, the implicitness of this expression prevents the direct evaluation of via (10), the unknowns σ s (ζ ) and γ ad (ζ , σ ) depending upon the Stokes stream function. The computation of is then carried out through an iterative and semi-analytical procedure thoroughly detailed in Conway (1998) and summarized hereinafter. Some improvements in the numerical approach are also introduced.
As detailed in Conway (1998), the method can deal with an H distribution of the polynomial type where the value at (ζ = 0, σ = σ ad ) is termed σ ad . If H( ) is supposed to be C 1 continuous then, from Equations (13) and (4), ω θ can be written as In order to integrate Equations (10)-(12) exactly, the tangential vorticity, appearing on the left-hand side of Equation (14), can conveniently be expressed through the following polynomial expansion (Conway, 1998): In fact, with the help of this alternative ω θ form, the radial integral in Equations (10)-(12) can be exactly evaluated, thus obtaining The quantities which appear in the above three equations, can be evaluated through recursive relations for each value of the integers ξ , μ and ν (see . In order to integrate Equations (16)-(18), a further difficulty has to be tackled. In fact the function σ s (ζ ) describing the shape of the wake is still unknown and must be determined as part of the whole solution. One way to proceed is as follows. Firstly, on account of the fact that the Stokes stream function is constant on the wake edge, one can exploit Equation (16) at σ s (ζ ). By doing so, an integral equation for σ s (ζ ) is obtained: At this point the solution procedure can be outlined as follows.
(1) Assume that provisional values of the discrete counterparts of the σ s (ζ ) and {A n (ζ )} n=0,N distributions are known for n z values of the axial coordinate.
(2) Then replace the continuous and ω θ distributions with their discrete representations on a grid spanning the wake with n zs × n rs points in the axial and radial directions, respectively. Specifically, integrate Equation (16) in the axial direction with an adaptive quadrature scheme to obtain the discrete distribution at all mesh points. Moreover, with the help of Equation (14) the azimuthal component of the vorticity can also be evaluated at all mesh points. (3) A least-squares minimization procedure is then applied to parametrize the ω θ distributions, computed at the previous step, leading to a new set of {A n (ζ )} n=0,N . (4) Update the slipstream σ s (ζ ) from Equation (20) evaluated at (ζ , σ s (ζ )). (5) Go to item (1) until convergence of the slipstream σ s (ζ ) is reached.
The unknowns σ s (ζ ) and are initialized neglecting the presence of the disc. Note that, in comparison to previously proposed implementations of the method (Conway, 1998;Greenberg & Powers, 1970), in the actual version of the semi-analytical procedure there is no need to approximate the updated slipstream shape σ s (ζ ) (see point (4)) through a least-squares approach. By so doing, the accuracy of the solution is improved further.
In the following, the SA-ADM is employed to simulate the flow around a propeller with J = 0.5 and characterized by a parabolic versus load distribution, i.e. M = 2 in Equation (13). From the same Equation (13), it can be readily understood that H(0) = a 0 . Moreover, in order to enforce a vanishing load at the disc edge and at the centre, the following two conditions have to be satisfied: a 0 = 0 and the sum of the a m parameters for m = 1, . . . , M has also to be zero, i.e. M m=1 a m = 0. Thus, for a parabolic distribution of the load, a 0 = 0, a 1 = −a 2 and whereb = σ 2 ad a 1 /(U ∞ σ ad ).

Comparison between the SA-ADM and the CFD-ADM
As discussed in the introduction, the CFD-ADM is a commonly employed analysis method that provides the solution of the Navier-Stokes or Euler equations in a domain comprising a propeller represented through a body-force distribution acting over a disc area. From a computational point of view, these distributions are obtained by pursuing several approaches, such as iterative and interactive coupling with lifting-surface, boundary element or blade element methods. However, many researchers have often applied a prescribed body-force distribution over the disc area (see for instance Dai, Gorski, & Haussling, 1991;Hoekstra, 2006;Piquet, Queutey, & Visonneau, 1987;Yang, Hartwich, & Sundaram, 1991). In this paper a CFD-ADM method is implemented in the popular CFD suite Fluent by introducing source terms in the axial and tangential momentum equations. In the following section, the results obtained by modelling the rotor through the CFD-ADM are compared with those of the semi-analytical procedure described in Section 2, assuming identical body-force distributions. To this end, an error analysis is conducted for both methods in Sections 3.1 and 3.2. Finally, in Section 3.3 a detailed comparison between the results of the two procedures is carried out.

Analysis of the error for the SA-ADM
As reported in Conway (1998), the evaluation of the Stokes stream function via the exact solution of the governing equation (10) is prevented due to the implicitness of this equation. In particular, the slipstream location σ s and the density strength γ ad all depend upon . The solution can be made explicit through the semianalytical procedure developed by Conway (1998) and further improved by Manna (2013, 2014). As detailed in Section 2.2, this procedure introduces three numerical parameters, i.e. n z , n zs , n rs , whose effects on the solution are investigated in the following error analysis. In addition to these parameters, the effects of some other numerical quantities are also investigated. One of these is the error related to the numerical integration of Equations (16)-(18). This task is accomplished with the help of an adaptive quadrature scheme (Lyness, 1970), which allows a required tolerance value (quad_err in Table 1) to be prescribed in advance. Another numerical parameter is the degree N of the polynomial appearing on the right-hand side of Equation (15). Finally, the last numerical parameter involved in the semi-analytical process is the value of the residue res which has to be achieved to exit the iterative procedure. The residue is defined in terms of two consecutive slipstream radii evaluated at ζ = 15σ ad : For the sake of conciseness, all numerical parameters are simultaneously made to increase in a six-level range as shown in Table 1. Assuming case SA-6 of Table 1 as reference, the relative errors on the performance coefficients and the discrete L 2 errors for u are reported in Table 2. The following classical definitions are adopted for the performance coefficients: where T is the thrust experienced by the rotor, P is the power transferred by the propeller to the fluid and η is the propulsive efficiency. On inspecting Table 2, it is immediately apparent that the errors are monotone decreasing functions of the refinement parameters which are simultaneously increased from SA-1 to SA-6. Furthermore, from a practical point of view, it should be observed that the error magnitude, even for the coarsest combination of   Case no.

Grid convergence analysis for the CFD-ADM
The CFD suite Fluent is a general purpose and widespread analysis software that can handle a variety of complex problems. In particular, it is possible to activate, in a disc shaped region, the source terms in the axial and tangential momentum equations to simulate the presence of the rotor. With this aim, a 2D computational domain is generated (see Figure 3) and discretized through a structured mesh. The domain is bounded by an inlet and an outlet boundary placed at ∓25σ ad from the plane of the disc, and by an upper wall located at 50σ ad from the ζ -axis. At the inlet a uniform axial velocity U ∞ is imposed, whereas the outlet is treated as a fully-developed flow boundary. An inviscid wall treatment is used for the upper surface. Finally, the axial and the tangential body-force radial profiles, obtained for the reference case SA-6, are prescribed at the actuator disc plane as source terms in the momentum equations. In particular, the computational domain is subdivided into an inner and an outer domain. The first one is the 10σ ad × 3σ ad region around the propeller (see Figure 3) and it is characterized by a uniform density of the mesh. To reduce the computational cost, in the outer domain the density of the mesh is gradually decreased going away from the inner domain. The effects  of the extension of the inner and outer domains are also investigated through a separate study whose results are omitted for the sake of briefness. The Euler equations are solved employing the well-known SIMPLE algorithm and a second-order discretization scheme. A grid convergence analysis is carried out employing four different structured meshes with an increasing number of nodes (see Table 3). In more detail, the uniform discrete spacing h of the inner domain mesh is doubled moving from the coarser case (CFD-1) to the finer one (CFD-4). Consequently, the number of cells quadruples at each grid level. The relative errors in the performance coefficients and the discrete L 2 errors in u, v and w are shown in Figure 4, where the CFD-4 solution is assumed as reference. The overall trend of the grid convergence study is more than satisfactory. In fact, the rate of error decay is very close to the theoretical value of two. With only the exception of the thrust coefficient, the mesh density appears adequate for the truncation errors to lie in the so-called asymptotic range.

Comparison of the results
In this section, the results of the SA-ADM and the CFD-ADM will be compared with each other both in terms of local and global quantities. Obviously, the most refined solutions, i.e. SA-6 and CFD-4, will be employed. Also note that, until now, the SA-ADM has been verified only at downstream infinity (see Conway, 1998) through an asymptotic approach. Hence, the comparison reported hereinafter also aims at verifying the SA-ADM in the whole computational domain. The analysis begins by looking at Figure 5, which reports the radial profiles of the radial, axial and tangential perturbation velocities evaluated through the SA-ADM and the CFD-ADM. As can easily be understood from this figure, the differences between the results of the two methods are extremely small. In particular, no significant differences seem to be present in any velocities components for ζ /σ ad = −1 and ζ /σ ad = 1.
In order to extend the assessment of the two proposed numerical methods to field data, the axial and radial velocity contours are compared in Figure 6. The top half part of the figure shows the contours obtained with the SA-ADM, while the bottom half is related to  the CFD-ADM. Figure 7 reports the traces in the (ζ , σ )plane of the streamlines. Moreover, Figure 7 also highlights the possibility of fully taking into account the wake contraction through the SA-ADM. The last two Figures 6  and 7 show that the agreement between the two methods is extremely good in the whole domain. However, some small differences also appear, especially in the proximity of the wake edge. This discrepancy can easily be detected by looking at the shape of the axial velocity iso-line u/U ∞ = 0.9 in the wake edge region.
To highlight these aspects further, Figure 8 reports a close-up view of the axial velocity radial distributions for different axial stations. As clearly shown, the discontinuity in the radial derivative of the axial velocity, which takes place at the wake edge, is smoothed in the CFD-ADM; a phenomenon whose impact on the solution could be relieved by locally refining the mesh.
Finally, in Table 4, the differences between the two methods are presented in a quantitative fashion both in terms of performance coefficients and velocity components.
From the wide comparison carried out between the SA-ADM and the CFD-ADM it can be concluded that the agreement between the two methods is excellent in the whole computational domain. This means that they could be used equivalently in the analysis of propellers. However, the methods are obviously characterized by very different computational costs; a fact that makes the   Table 4. Comparison between SA-ADM (SA-6) and CFD-ADM (CFD-4).

C T C P η
Relative errors (%) 9.183 × 10 −2 1.282 × 10 −1 3.637 × 10 −2 ζ /σ ad = −1 ζ /σ ad = 1 ζ /σ ad = 2 L 2 errors u/U ∞ 2.655 × 10 −4 3.981 × 10 −3 3.932 × 10 −3 v/U ∞ 2.589 × 10 −4 1.634 × 10 −4 3.533 × 10 −5 w/U ∞ -3.897 × 10 −3 4.254 × 10 −3 SA-ADM extremely attractive for design purposes. This is particularly true in the first stages of a design procedure when a very fast analysis tool is preferable to sweep the design space quickly. The computational cost related to the SA-ADM method is quantified in Figure 9 which, for all cases (see Tables 1 and 3), reports the values of C P , η and of the CPU time. All simulations have been performed on an Intel Xeon CPU E5-1620 v2 3.70 GHz. Figure 9 also shows a shaded area representing a 1% (±0.5%) deviation band from a reference value evaluated as the average of the SA-6 and CFD-4 results. Figure 9 obviously verifies that the computational cost of the CFD-ADM is bigger than that of the SA-ADM. Moreover, on inspecting these figures it can easily be inferred that, in order to comply with a prescribed accuracy of ±0.5%, the solution SA-2 should be selected for the SA-ADM. The run time of this SA-2 case is seven seconds, a value that can be considered small enough for a preliminary design procedure.

Conclusions
This paper has presented the application of the SA-ADM to the analysis of marine propellers. Although the method relies on the simplified assumptions of axisymmetric and inviscid flow regimes, it still provides significant improvements when compared to the widely employed momentum theory or to less advanced actuator disc models available in the literature. Most noticeably, the method allows the convergence of the wake to be accounted for, an important feature that is disregarded in all linear approaches. Additional properties concern the ability of the method to deal with heavy loads of arbitrary radial distribution.
The method has been applied to propellers characterized by a parabolic load distribution whose performance, both in terms of global and local parameters, is compared with those obtained through a CFD actuator disc model. The latter method, which is frequently used to analyse the flow around marine propellers, represents the effects of the impeller through a set of axisymmetric body forces. The comparison relies on a thorough analysis of the errors of the two methods providing a set of reference solutions characterized by a controlled accuracy. An excellent agreement has been documented between the results of the two methods while the computational complexity is obviously very different. For these reasons, the semi-analytical method is the ideal candidate tool to be adopted in a preliminary design procedure based on the repeated analysis scheme. Furthermore, the collected data may be of interest for the in-depth verification of other numerical approaches.
Ongoing research deals with the assessment of the performance of the method when an arbitrarily shaped hub is introduced. The coupling of the SA procedure with a blade resolving modulus providing the propeller load is also envisaged.

Disclosure statement
No potential conflict of interest was reported by the authors.