Design Method for Contra-Rotating Propellers for High-Speed Crafts: Revising the Original Lerbs Theory in a Modern Perspective

The main theoretical and numerical aspects of a design method for optimum contrar-rotating (CR) propellers for fast marine crafts are presented. We propose a reformulated version of a well-known design theory for contra-rotating propellers, by taking advantage of a new fully numerical algorithm for the calculation of the mutually induced velocities and introducing new features such as numerical lifting surface corrections, use of an integrated modern cavitation/strength criteria, a modified method to consider different numbers of blades among the two propellers, and to allow for an unloading function in the search for the optimal circulation distribution. The paper first introduces the main theoretical principles of the new methods and then discusses the influence of the main design parameters on an emblematic example of application in the case of counter rotating propellers for a pod propulsor designed for fast planing crafts (35 knots and above).


Introduction
In this last decade, the general interest in counterrotating propellers (CRPs) has increased and different manufacturers are developing new propulsion systems based on double propellers, especially for their use in fast planing crafts and pleasure yachts (0.5-2 MW), but also in larger vessels, in the high-power range (2)(3)(4). Of course, there are other fields of applications where a CRP propulsion arrangement might be beneficial: for instance, electrical podded drives for large high-speed displacement navy ships or ro-ro ships for high-speed transportation.
Despite this broad range of application, only a limited number of theoretical design methods have appeared in literature about contrarotating propellers; in fact, apart from different practical semiempirical approaches, the early method proposed by Morgan in 1960 remains still a valid option to design their optimum wake-adapted geometries.
In this paper, we present the main theoretical elements of a modified CRP design method based on the general framework of Morgan but reformulated its most simplified passages, taking advantage of modern computational techniques, with the final intent to produce valid CRP designs for pod propulsors of high-speed planing crafts, as the one presented in Figure 1, designed to drive planing crafts at top speed of 35/45 knots with maximum input powers in excess of 1 MW.
The original design method initially developed by the authors [1] has been further developed and improved, over the years, in its theoretical aspects, resulting in the following main features.
(i) Morgan's theory combines the analytical computation of self-induced velocities at propeller disk given by Lerbs [2], with a numerical treatment of axial interference factors (which states the relationship International Journal of Rotating Machinery between disk and wake axial velocity components) given by Tachmindji [3] for an infinitely bladed propeller. The present method adopts a full numerical lifting line model (with slipstream contraction effect) for exact representation of wake velocities.
(ii) Cavitation and strength constraints are imposed by means of an iterative and automatic blade geometry optimization method, which verifies the margin on the maximum admissible local stress and cavitation inception. The adopted method was inspired by the lifting line/lifting surface propeller programs developed in the 80s by various propeller designers and researchers in Italy [4][5][6].
(iii) In Morgan's original method, lifting surface corrections are addressed following the well-known Ludwieg and Ginzel et al. [7] parametric formulations that permit to calculate the value of the camber just at mid chord, for a given standard chordwise circulation distribution. The present theory applies the numerical lifting surface corrections for camber and pitch based on a 3D vortex/source lattice model [8] that leaves great freedom in terms of assumed chordwise load distribution and treats the problem of the vortex wake alignment and contraction with consistent numerical schemes, validated in the case of a single propeller in [9].
(iv) Morgan's original theory can only consider optimum spanwise circulation distributions, while present theory offers the possibility to apply a nonoptimum circulation curve [10] in order to unload blade sections especially at tip and root where their working condition is often very close to cavitation limits, in particular for fast planing boats. Such kind of cavitation can lead to blade erosion and unpleasant phenomena like radiated noise and induced vibrations of stern structures. The function that describes the unloading factor at each blade radial position is applied to the optimum circulation distribution of the equivalent propeller; hence it cannot be customized for each of the two CR propellers, but equally applied to the equivalent one.
Having synthetized the main technical features of the present design method, we go on, in the following sections, discussing the key aspects of its theoretical formulation.
Finally some examples of practical applications of the design method on two different cases involving fast-planing boats with pod propulsors are presented and the results obtained are analysed in order to discuss the effect of a systematic perturbation of the main input design parameters upon the resulting blade geometry.

Fundamental Assumptions and "Equivalent Propeller"
Concept. As already mentioned before, the present design procedure is a revisitation of the method developed by Morgan [11], which makes use of the analytical lifting line theory of Lerbs for a single propeller [2]. The design theory of two contrarotating propellers is more complicated than that of a single one, due to the mutual hydrodynamic interference between the two propellers; for this reason some engineering simplifying assumptions were introduced by Morgan in order to reduce the problem to that of a single propeller and conveniently exploit the analytical solution of Lerbs. This section gives a brief description of the Morgan's theoretical framework whose fundamental hypotheses can be summarized as follows.
(i) Both the propellers are to rotate at the same rpm.
(ii) Each propeller delivers 50% of the total thrust and absorbs 50% of the total torque.
(iii) Both propellers are moderately loaded, and the wake does not change appreciably. Along the axis.
(iv) The self and mutual interference of the vortex sheets shed by each propeller blade is negligible.
(v) Blades are represented by lifting lines, and lifting surface corrections are introduced as a correction a posteriori.
(vi) Any unsteady effect is neglected; for example, mutually induced velocities (fore and aft blades) are supposed to be time averaged.
(vii) The bound circulation on the blade can be represented by a sine Fourier series.
The key point of this idealized approach lies in the definition of the so-called equivalent propeller, that is, an optimum "virtual" propeller which produces 50% of the total required thrust and absorbs 50% of the total torque, having the hydrodynamic pitch angle equal to the average of the hydrodynamic pitch angles of the two actual propellers. In this way, to find the optimum circulation, Morgan initially reduces the problem of the two CR propellers to a single equivalent one, by ideally bringing the longitudinal distance between the two propellers to zero (the two propeller discs are located at the same longitudinal coordinate), so that the mutually induced velocities become independent from the longitudinal distance of the propellers and their value can be found as function of the bound circulation by the induction factors given by the classical Lerbs theory for single wake adapted propellers.
The velocity components diagram, considered for the fore and aft propellers, is reported in Figure 1.
Using momentum theory, considering symmetry conditions and averaging the value of the induced velocities over a complete blade turn, the following relationships are found International Journal of Rotating Machinery 3 Figure 1: ZF 4000 Large Pod, using contrarotating propellers designed with the presented method. between self-induced and mutually induced (interference) velocities: where the subscripts have the following meanings: 1 = at forward propeller, 2 = at rear propeller, a = axial component, t = tangential component, s = self-induced velocity, i = interference induced velocity.
Moreover, since unsteady effects are not taken into account, the factors f a and f t for obtaining the circumferentially velocity components averages are needed, while factors g a and g t explicitly express the influence of the actual axial distance between the two propellers and of the slipstream contraction, respectively, in a global sense.
Using Stokes' law, Lerbs derived the following expression for factor f : where G * is the total nondimensional circulation, V s is the ship speed, and x represents the nondimensional radial position of the considered blade section. As mentioned in the introduction, the factors, g a and g t , are calculated by Morgan through the simplified model of Tachmindji, which considers an infinitely bladed propeller without hub (zero hub radius). Moreover, Morgan's theory simply considers (g a ) 1 = (g a ) 2 .
As described later in this section, the present work makes use of a fully numerical method for calculating the induction factors which overcome the limitations deriving from the crude assumptions made by Tachmindji to derive an analytical expression for them.
In principle, the different velocity components acting on the forward propeller are And similarly for the rear one, where V a stands for the (averaged) speed of advance.
Under the leading assumption of the equivalent propeller, that is, longitudinal distance between propellers zero (g a = 0, g t = 0, and V a1 = V a2 = V a ), equal torque of both propellers (u s1 = u s2 ), (3) and (4) can be rewritten as follows, to find the axial and tangential velocity components and the hydrodynamic pitch angle on the fore propeller: and on the aft one, Hence, the hydrodynamic pitch for the equivalent propeller, calculated as the averaged value on the front/rear ones, can be derived by neglecting second-order terms as follows: Following Lerbs' [2] work, it is possible to approximate the integral equation for induced velocities by means of the following expression: where Z = propeller number of blades, x h = nondimensional hub radius, h a m , h t m = Lerbs' improper integrals function of (β i ) eq and radial position.
In order to reduce the analytical problem to the numerical solution for a discrete number of unknowns, that is, the circulation amplitude coefficients appearing in (8), the nondimensional (continuous) circulation over the lifting line 4 International Journal of Rotating Machinery of the equivalent propeller has been approximated through a sum of n odd sinusoidal functions: where ϕ is related to the nondimensional radial coordinate x through the following equation: When (2) for circumferential factor f is used into (7) for tan (β i ) eq (kinematics condition), in combination with (8), and (9) for explicating induced velocities and total circulation, an integral-differential equation is obtained, whose unknowns are the Fourier coefficient G * m : By writing the previous equation in n radial points, it is possible to solve an n-sized system in the n unknowns G * m . Since the term u a /u t is a function of the total circulation G * which, in turn, depends on the hydrodynamic pitch angle β eq , (11) represents a nonlinear system of equations and must be solved by an iterative process. As a first guess for the value of u a /u t , one can use the following simplified relation: where (tan β i ) eq can be calculated as in which tan (β i ) eq = V a /ωr, w 0 = effective wake fraction, w x = local wake fraction, η I = ideal efficiency of a single propeller. This ideal efficiency η I can be found solving Kramer's equations when the following design parameters are known: T = propeller thrust, Z = number of propeller blades, D = maximum propeller diameter, n = propeller revolutions (equal on both props), A E /A O = expanded area ratio (initial estimation).
With this initial value for u a /u t , system (11) is solved (after having numerically approximated Lerbs' improper integrals h a m , h t m ) and the total circulation G * is determined. Then a new value for induced velocities ratio u a /u t is derived through (8) and compared with the original one; if the absolute difference between the value obtained in the current iteration and that of the previous iteration is outside a certain tolerance, system (11) is solved again with the updated value. This iterative scheme goes on until convergence is achieved.
Once the system (11) is solved and total circulation G * for the equivalent propeller is found, the ideal inviscid thrust coefficient C * Tsi is calculated from the following equation, which is derived integrating along the radius the contribution to thrust of the lift force produced by each blade section expressed by direct application of the Jukoski theorem: The ideal thrust calculated in this manner is then compared to the requested value given as input T opportunely corrected to compensate viscous losses. If those values do not agree conveniently, a new value for (tan β i ) eq is used for the solution of (11) at each radial position. The new tentative value of this parameter is estimated by multiplying the old value by a scalar factor directly proportional to the difference between the requested and the calculated thrust. Again this higher-level iterative procedure is repeated until the convergence is found on the given thrust. Usually three/four iterations are sufficient to achieve a satisfactory convergence limit.
If fore and aft propellers have a different number of blades, the aforementioned procedure is computed twice since two equivalent propellers with two different blade numbers need to be calculated.

Nonoptimum Load Circulations.
As already mentioned, through the limit concept of EqP, it is possible to reduce the problem of the simultaneous design of two coupled CR propellers into a more conventional problem of design of a single optimum propeller. The well-known Lerbs [2] lifting line design method for single propeller, valid also in the case of nonoptimum circulation distribution, has been adapted in this work also to the equivalent propeller. The application is not straightforward and requires some modifications, in relation to the peculiar hydrodynamic operation of a CR propellers set.
The optimum circulation distribution calculated with the method outlined in the previous paragraph is modified by an unloading function λ G , variable across the radius. Special care must be taken to ensure a very smooth behavior of this unloading curve that is reflected on a regular distribution of a modified final circulation, and so to avoid unrealistic peaks predicted by the induction factor theory. Applying the variable unloading factor λ G to the optimum circulation G * , the modified circulation distribution F * is found. Also this new curve, similarly to its parent function, can be expressed in terms of a sine series, in which the n terms F * m are to be found by the numerical procedure. In formulae, International Journal of Rotating Machinery 5 As for the optimum propeller, the value of the design ideal thrust should be met by scaling this new generated load distribution through a constant factor k. The final unloaded circulation distribution so obtained will be referred to as L * and defined as through (15)- (16), it is possible to express the formula for the final ideal thrust coefficient from (14), valid for the unloaded counterrotating propellers, keeping the usual Morgan assumption of the equivalence between the two induction factors f a e f t : introducing the k scale factor and differentiating the thrust coefficient with respect to radial coordinate, the following expression is obtained: which, as before, integrated along the radius becomes This (19) must be solved with respect to the unknown scale factor k and hence can be rearranged as While the circulation F * (and its sinusoidal component amplitudes F * m ) is known and kept constant during the searching of the solution, the unknowns to be found by an iterative converging algorithm are the scale factor k and the functions h t m , h a m , which in turn depend on the hydrodynamic pitch angle of the trailed vortical wake β i , at each radial position x considered in the calculation.
With reference to the velocity diagrams of Figure 1, yet following the Morgan EqP model, the relations between the self-and mutually induced velocities of (7) can be rewritten as a function of L * , as follows: in which The solution method developed by the authors and implemented in the relative design code is detailed in the following. The value of k factor is calculated using an iterative procedure in which at the first step the system (20) is solved assuming the induction factors h a , h t , i a i t , the axial and tangential velocity components V a , V t , and the pitch distribution β i , previously found in Block A and valid for the optimum propeller. Contemporary, as proposed by Lerbs [2], the induced velocities are neglected.
Once (20) is solved for the scaling factor k, the corresponding hydrodynamic characteristics and the velocity field generated by the unloaded circulation L * (16) can be found through the following expressions: The knowledge of the new induced velocities (23) will permit the calculation of the hydrodynamic angle β i and hence the new induction functions as defined by Lerbs [2]. This procedure is repeated until the current value of the ideal thrust coefficient C * Tsi (19) meets the design value given as input. Usually three or four iterations are needed to converge within a sufficient tolerance.

Load Distribution for Actual
Propellers. The next phase of the design process deals with the calculation of the actual propellers, starting from the equivalent propeller just solved.
Lerbs derived the following relationship between the selfinduced velocities: International Journal of Rotating Machinery where δ and ξ are the average contraction ratio of slipstream δ at the aft propeller and the average circulation factor, used to ensure equal torque on each propeller, respectively.
Substituting these values into (1) for the axial and tangential relative velocities, it is possible to derive the following expressions for the hydrodynamic pitch angles tan (β i ) 1 and tan (β i ) 2 : where The aforementioned values δ and ξ are obtained from the following formulae: An iterative approach is needed to solve the two equations for δ and ξ, because of the coupling terms. The procedure, then, starts by assigning zero to both of these parameters (hence their averaged values δ and ξ are also zero); then it calculates their actual values by (27) and computes the new averaged value δ and ξ until convergence on such parameters is achieved. Factors δ and ξ are then applied to (25) and (26) to calculate the hydrodynamic pitch angles tan (β i ) 1 and tan (β i ) 2 for the Actual Propellers.
With these angles, the total circulation G * (or the modified F * ), and the induced velocities u a1 , u t1 , u a2 , u t2 , it is possible to derive the load distribution along the blade with the following equations: in which C L = section lift coefficient; c = section chord.

Calculation of Blade Outline and Section
Shape. The next stage deals with the determination of the blade geometry in terms of chord length, thicknesses, pitch, and camber which ensure the requested section lift coefficient while ensuring cavitation and strength constraints. There are several approaches when solving the foregoing problem, but the authors follow the works by Connolly [12] for the calculation of blade stresses and the method proposed by Grossi [6] for cavitation issues; the last procedure is based upon an earlier work by Castagneto and Maioli [4] where minimum pressure coefficients on a given standard blade section are computed.
The semiempirical simplified blade model proposed by Connolly considers radial stresses σ r due to bending moment and radial stresses σ c due to centrifugal forces: the former are determined by developed thrust and absorbed torque, while the latter are determined by propeller revolutions. These quantities are related to propeller geometry and dynamic parameters by means of the following equations: where ct 2 represent thickness ratio product (equal to blade sectional modulus), T and Q are thrust and torque, respectively, tabulated coefficients A 1 , A 2 are function of nondimensional radial position, P is section blade pitch (φ is pitch angle), and i stands for tangent of the rake angle (positive if aft rake). Total radial stress is given by the sum of σ r and σ c , and its value must be lower than the maximum admissible stress value σ adm , that is: where σ Y = yield stress (for nickel aluminium bronze alloy equal to 550 MPa); K S = safety factor to take into account uncertainties of an unsteady load on the propeller and to allow for fatigue effects; K rob = tuning factor to calibrate the relative weight of the strength criteria over cavitation criteria. Substituting (31) into (29) leads to the following strength criteria equation: Having regard to cavitation criteria, the present theory makes use of a simplified cavitation inception rule, where minimum pressure coefficient C min P on each blade section has to be higher than the cavitation number σ I 0 times a safety factor K P used to ensure a certain design margin (usually K P <1): For minimum pressure coefficient estimation, the authors follow Castagneto and Maioli's work where thin profile theory for standard propeller foils and mean lines is corrected by the results of experimental tests done by them at the cavitation tunnel: where h 1 , h 2 , h 3 are given by Castagneto and Maioli in a tabular form being only function of thickness and camber distribution, and C L f , C Lα are lift coefficient developed by blade section by camber and angle of attack, respectively. The foregoing equation can be rewritten in order to make use of known quantities as C L c (given by (28)) and to consider for (33): where and p = C L f /(C L f /C Lα ) express the proportion of lift coefficient developed by camber on the total lift and are selected by the propeller designer.
Since the cavitation-strength problem (illustrated with the foregoing equations ((29)-(35)) depends on interrelated unknowns, this is solved by making use of an integrated approach described in the following.
At a first step, centrifugal radial stress σ c (which depends on the unknown blade thickness) is set to zero and equation (32) is solved for the blade section modulus ct 2 ; then the third-order (35) is solved to determine the unknown quantity xso obtaining chord length cand, in turn, blade thickness t.
Once this first guess for thickness distribution is derived, it possible to evaluate the actual value for σ c (30) and the foregoing procedure is repeated until satisfactory convergence is achieved on the parameters involved.
(u as ) 2 (u as ) 2 Some additional constraints are needed to ensure a minimum blade thickness at tip sections (usually 3% of propeller radius) and a maximum thickness over chord ratio at blade root (normally 15% or less).
The designer is allowed to define different cavitation and strength margins for fore and aft propellers by setting different values for K P and K rob .

Lifting Surface Corrections.
A numerical procedure for determining the lifting surface corrections is presented following the work by Greeley and Kerwin [13]. The aim is to find the correction to a first guess camber and pitch of the blade sections calculated after the lifting line design program, ensuring a given load distribution over the chords and the same total design thrust. A Cartesian coordinate system is fixed on the propeller; the x axis is the axis of the propeller shaft with its positive pointing upstream, the y axis is an axis on the first blade arbitrarily selected so as to pass through the mid chord at r = 0.7 R with its positive outward, and the z axis is the third axis and has its positive in accordance with the righthand rule. It is also possible to define a cylindrical coordinate system (shown in Figure 2) in which a point in space can be defined by (r, θ, x) where r is the radial coordinate, θ is the angular coordinate, and x is the same as defined above.
Fundamental assumption for the development of the theory can be summarized as follows. (ii) Sources/sinks are sued to allow for thickness effects.
(iii) The lattice is located on the mean surface of the propeller.
(iv) The effect of slipstream contraction and wake alignment in terms of pitch and radial distribution of the free vortices is considered following Greeley-Kerwin [13] in the more generalized version of Grassi and Brizzolara [9].
(v) The fluid is inviscid and incompressible.
(vi) The flow is steady and axisymmetric. The strength of the bound circulation (found with the lifting line model) is given by θ r (r) and θ l (r) express the angular position of the blade section leading edge and trailing edge in a cylindrical reference system of Figure 3; G r is the nondimensional circulation distributed on the blade mean surface, on a structured mesh of discrete (concentrated) vortex segments and G is the total nondimensional spanwise circulation defined as The strength of the helical free vortices can be easily calculated as stated in the following: By combining (37) and (39), one obtains the relation between the strength of the free vortices and the strength of the spanwise vortices: It is now possible to calculate the strength of the chordwise vortices on the blade surface: As previously stated, the continuous distribution of circulation and sources described above is replaced by a discrete lattice of vortex/sources elements through a sub-division in the radial and chordwise direction. In order to align the mean surface of the blade with the local flow velocity, the calculation of the velocity V i induced by this vortex/sources mesh is needed: by the application of Biot-Savart, law one obtains where S and dl are, respectively, a vector from the vortex segment to the point P and the elementary vector tangent to the vortex line.   Figure 4: Chordwise coordinate transformation Ψ (as from [13]).
By further defining A 3 as the area between the trailing edge and a generating line along which a lifting line would be placed, it is possible to rewrite (42) as follows: Since this gives the induced velocity by a single blade, it is necessary to sum all the contributions due to the Z propeller blades.
The bound circulation distribution Γ r in the previous equations can be expressed by the following equation: where γ(r) = Γ(r)/(θ τ − θ λ ); C 0 , A 0 , and A n are constant coefficients and, in general, are functions of r; Ψ is the new chordwise coordinate and is defined as per Figure 4: The first term of (44) represents a constant load distribution due to angle of attack; the third, some arbitrary distribution in the form of a sine series with amplitude functions A n . The first three coefficients are coupled by the following relation: so that (37) is satisfied. The effect of thickness is then addressed by introducing a source-sink system distributed over the blade as a lattice whose strength σ(r, θ) is calculated by means of the classical linear approximations from airfoil theory: where H s is the velocity induced at the point P by a unit source located at point r, θ.
Once the total-induced velocity V i tot = V i + V i thickness is computed, it is possible to prescribe the boundary condition to be satisfied on the blade surface: where α(r) is the angle of attack, f p is the camber along the chord x c , V r is the resultant inflow velocity to the blade section, U is the resultant-induced velocity from lifting line theory, and (V i tot ) n is the induced velocity normal to blade chord. To determine α, the following integral is computed starting from the leading edge to the trailing edge, being f p equal to zero: Then a new nose-tail line can be drawn based on the new pitch and a new camber distribution along the chord can be computed.
From the numerical point of view, the foregoing boundary condition is implemented through the computation of the total-induced velocity at a defined number of control points; the actual lattice geometry is then displaced accordingly along cylindrical sections. The vortex/source lattice is then recomputed and induced velocities updated until a satisfactory convergence is achieved on the geometry itself.
An iterative approach is needed performing two computing loops: an inner loop to align the singularities located on the blade and wake with the resulting flow and an outer loop which takes care of adjusting the circulation distribution in order to deliver the prescribed thrust.

Calculation of Induced Velocities with a Fully Numerical
Lifting-Line Method. In order to calculate the entire velocity field in the wake of the first propeller, a fully numerical lifting line model with slipstream contraction effect has been included in the design methodology. The expression for induced velocities obtained by Lerbs, in fact, holds only at the propeller disk and makes it impossible to calculate exactly the induced velocities in the wake of the forward propeller. Morgan made use of the Tachmindji infinitely bladed propeller model to calculate the axial velocity component induced by each propeller on the other. In the present work, a more accurate model has been introduced, consistent with the lifting line method used to design wake-adapted propellers, where induced velocities are computed by making use of the Biot-Savart law to the discrete bound vortex elements and to the free trailing vortices. The velocity field induced on a point P (defined by vector r = P − P 0 ) by a three-dimensional vortex d l = (P − P 0 ) of strength Γ can be derived through a modified version [14] of the well-known Biot-Savart law: where According to the foregoing expression, it is possible to derive the velocity vector induced at a point P by a lifting line (discretized by N parts of length Δ = R/N) and by a trailing helical vortex line shed in the wake by the lifting line, once defined a cylindrical coordinate system (r, Ψ, z) fixed on the propeller, as represented in Figure 5. This vortex line has the following geometrical parameters: r 0 = radius, β i0 = pitch angle, r 0 M = longitudinal extension.
Which are related through the following equation: This helical vortex line can be divided into N regular intervals, so obtaining N + 1 singularity elements located at points: (n − 1) n = 1, 2, . . . , N + 1. (53) In order to achieve a good numerical convergence on the computed velocities, the parameters M and N have to be related to the pitch angle of the trailed vortex line β i0 , the number of blades g, and the ratio r/r 0 , which defines the relative position between the trailing vortex line and the radial position of the point where the velocities are evaluated. Once defined the basic vortical elements, it is possible to model the velocity field induced by a propeller by superimposing several horseshoe lifting lines composed by a bound vortex lying on the blade and two trailing vortices directed downstream. Lifting line circulation is related to the free vortex circulation by dividing the lifting line itself into S equal intervals (so obtaining S + 1 trailing vortex) whose intensity can be written as follows: Given the bound circulation distribution Γ s (r), the hydrodynamic pitch angle distribution β i0 (r) (both computed for the actual propeller as described in Chapter 2.3), the number of blades g, the propeller, and hub radius R and r h , the numerical procedure calculates the velocity induced on the point P through the described pure numerical model.
A comparison can be made between numerical calculations and analytical formulation (Lerbs' theory) which is valid rigorously only on lifting lines and on bisector lines of the propeller disk between them. Figure 6 summarizes this comparison in terms of tangential and axial induced velocities (u t , u a ) on the lifting line, showing a very good agreement between analytical formulation and numerical approach.
The same analysis can be applied to calculate the average factors ( f t , f a ) derived by Lerbs through (2) and the interference factors (g a ) computed following Tachmindji [3]. The average factors can be defined as follows: where average velocity components (u a , u t ) are defined as follows: Integrals in (56) can be computed numerically by evaluating N + 1 velocity components at equally spaced angular positions from Ψ = δ to Ψ = 2π/Z − δ: computation at Ψ = 0 and Ψ = 2π/Z should be avoided in order to eliminate the numerical instability near the lifting line (except for the propeller disk).  Comparison between numerical and analytical formulation for f t and f a is shown in Figure 7 for one of the application cases described in the following (see Section 4) where a very good agreement between distributions confirms the validity of the numerical model developed. Moreover the  initial assumption about axial and tangential factors ( f t = f a ) is confirmed.
Interference factor g a (defined as the ratio between the circumferential mean axial component induced by a given propeller at a certain radial position on its disc and the induced components at a different axial position z-see (58) and (63)) can be computed only through a numerical  approach; Lerbs analytical formulation, in fact, holds only at propeller disk: As described in previous Section 2.1, Morgan follows a simplified approach by Tachmindji which holds for optimum infinitely bladed propeller model without any hub effect; in this formulation interference, factors are related to x, z/R, and λ i , where It is interesting to compare interference factor computed through numerical formulation with that calculated by Tachmindji; results are shown in Figure 8. The difference between the two formulations near the hub can be easily explained considering that Tachmindji's infinite bladed propeller analytical approach distributes the circulation on the whole propeller disk, so the circulation distribution effectively goes to zero only at x = 0. The difference above x = 0.6, instead, is to be attributed to the infinite number of blades assumption inherent to the analytical formulation used by Tachmindji.
The fully numerical approach followed by the authors for evaluating interference factor allows a better estimation of blade pitch especially in the outer radii blade region.
The specific effect of the hub on the circulation and local blade geometry has been widely investigated and commented in another work [15] in which a fully numerical lifting line method with mirror images has been compared with the present method for the design of contrarotating propellers. As a global result, it was verified by a different method for the analyses (vortex lattice, panel method, and RANSE) that the hub effect is important for the pitch and camber distribution at the blade root but does not particularly influence the global hydrodynamic performance of the propeller.

Practical Importance of Lifting Surface Corrections.
To highlight the relative importance of the numerical or parametric lifting surface correction methods, the design code has been run two times on the same design case with the pitch and camber correction methods. The reference design case corresponds to the propellers built and tested in the cavitation tunnel described in the next section for validation purposes. Figures 9 and 10 present the pitch and camber distributions, respectively, as obtained by using the uncorrected present lifting line method; the approximated lifting surface corrections of Morgan, Silovich, and Denny, as parametrically defined by Oossanen Van [16]; finally by applying the fully numerical lifting surface correction method described in Section 2.5.
As evident from the Figures 9 and 10, large corrections are applied to the pitch angles of the fore propeller sections and to the camber of both propellers sections. The trend of corrections calculation with the numerical lifting surface method is in general to increase the pitch and camber at the root of the blades for both propellers. The (unrealistically) steep increase of pitch and camber of outer sections (x > 0.9) closest to the tip is a known drawback of the lifting surface method, caused by the rapid chord length reduction at this location, and it is normally truncated on smaller values at tip.

Experimental Validation
Two propellers designed with the proposed method have been built and tested at the cavitation tunnel of the University of Genoa with an arrangement, specially built for the purpose by ZF Marine [17].
The main characteristics of the tested propellers, designed with an early version of the present method, are given in Table 1.
Results of the experimental measurements are synthesized in Figures 11 and 12 that show the variation of the total thrust and torque coefficients versus the advance coefficient (calculated with the diameter of the fore propeller). All the coefficients have been normalized, dividing the measured or calculated thrust and torque by the given value used in the design.   The reported numerical predictions were obtained using a fully numerical lifting surface analyses method specifically developed for CR propellers [18]. The predictions made with the lifting surface method are very much consistent with the values given as input to the design code, giving practically the same values. The comparison with the measurements done at design point (J/J des = 1) highlights an underprediction of the total thrust by almost 10% with a better agreement on the prediction of total torque, lower than about 5% of the corresponding measured value. Higher thrust and torque continue to be measured, with respect to the numerical predictions, in general at advance ratio lower than the design one. The accuracy level of the experiments, though, is not easy to assess and some source of uncertainties remain in the calculation of the effective advance and the speed with two   CR propellers and a relatively large pod inside the cavitation tunnel.
Recent feedback from full scale applications of propeller sets designed with the present method on several high-speed crafts has been very satisfactory, confirming the reliability (within practical accuracy) of the developed numerical methods.

Parametric Design Cases
The design tool described in the foregoing paragraphs has been applied to six design cases in order to underline the influence of the input parameters/settings on the resulting propeller geometry. The target applications are two planing crafts equipped with pod propulsion and running at different top speeds: design cases are grouped into two classes (A and B) depending on target speed. Design input data common to each group are described in Table 2.

Design Group A.
Design group A shows the effect of strength K rob and cavitation K P margins, according to (31) and (33), on blade chord length and thickness distributions. The foregoing parameters have no major influence on pitch and camber distributions; therefore, these curves are not shown in the following graphs but they will be given in Section 4.2 where the effect of design speed, number of blades, and unloading factor on blade geometry will be discussed.
Input data for propellers of design group A are listed in Table 3.
Obtained blade geometry is shown in the graphs of Figures 13, 14, 15, and 16.    As already mentioned, a given maximum thickness over chord ratio (common to all design cases) is imposed at blade root in order to avoid the prescription of small chord lengths by the design method, due to the zero-circulation value at the hub: therefore, the computed thickness and chord-length values at this radial position essentially depend only on the prescribed t/c and resulting blade sectional modulus. This last parameter is not influenced by the cavitation margin as it is demonstrated in cases A-1 and A-3 that show the same thickness and chord-length at the hub: at the outer radial positions, case A1 (lowest cavitation margin) shows higher chord length values and lower thicknesses when compared with case A-3. On the other hand, case A-2 (higher strength margin) presents the largest chord lengths and thicknesses distributions at all radii.
Finally, the best efficiency, as expected, is achieved by case A-3 which was obtained by imposing the minimum cavitation and strength margins.

Design Group B.
As described in the foregoing paragraph, design group B has been generated with the aim of discussing the influence of design speed, number of blades, and the type of loading curve on resulting blade geometry and propeller performance. Input data for propellers of group B are listed in Table 4.
Resulting blade geometry in terms of blade outline and thickness distribution is shown in Figures 17, 18, 19, and 20.
As expected, the five bladed set results in lower chord length/thicknesses compared to the four bladed cases.     Unloaded and optimum circulation cases (B-1 and B-3) consistently show differences on designed camber and pitch in the area where the optimum circulation is unloaded, that is, at the tip and at the root, especially on pitch, as per Figures 21 and 22.
The resulting blade geometry in terms of pitch/diameter and maximum camber/chord ratio is shown in Figures 21,  22, 23 and 24: in order to discuss the effect of design speed on the mentioned distributions, also values obtained for Case A-1 (the 19 most representative of group A) are included.
The highest design speed (= minimum thrust) required in group A results in lower camber/chord length distribution and higher pitch/diameter when compared to group B.
Finally, the effect of a nonoptimum design circulation obtained by unloading the curve at the root and tip regions as shown in Figures 25 and 26 is to decrease pitch and camber ratios at inner and outer radii and to increase it around midspan (B1 versus B3).
The loss of propeller efficiency due to the unloading procedure is quite small compared to the original load distribution, but the effective reduction on local cavitation and strength of hub and tip vortexes is important.

Conclusions
An integrated method based on lifting line/lifting surface theory for the design of modern high-speed counterrotating propeller has been presented in its principal theoretical aspects. New theoretical elements have been introduced with respect to the classical theory of Morgan/Lerbs, namely, the exact lifting surface correction, the exact calculation of mutual velocity induction coefficients, a new way of imposing strength and cavitation criteria, a new method to consider different blades number on the two propellers, and an unloaded circulation distribution.
The method has proven to be very effective in the design of cavitation-free counterrotating propeller optimized for efficiency, as presented in the validation study by means of model tests in Section 3.
The effect of main design parameters introduced in the new design method has been discussed through application examples, in particular the effect of the strength and cavitation margins and the effect of an unloading function applied to the optimum circulation distribution. Additionally, the effect of the blade number and design speed is also considered and analysed through dedicated design cases.
In principle (according (32) to (35)), strength and cavitation margins are both influencing blade sectional modulus: however, from the numerical point of view, the influence of the strength margin is much more significant as demonstrated by design cases A1 and A3. Accordingly, once the same strength margin level is set, the effect of a change of the cavitation margin is to increase chord lengths and, at the same time, to reduce thicknesses (since blade sectional modulus has to maintain its value regulated by the strength margin).
As expected (referring to Section 4.2), the effect of an increased number of blades is to reduce blade chord/thickness ratio and the pitch/diameter camber/chord ratios, due to a requested lower load on each blade.
The introduction of an unloading curve (which results in reduced pitch/diameter and camber/chord ratio at root and tip radii) has proven to be effective to avoid cavitation at both blade extremities, without significantly affecting open water efficiency.
Experience achieved so far with the new design method have been very satisfactory, and its use is encouraged by the authors, as recently verified also at full scale in a recent installation on a planing pleasure yachts with two pod propulsors. This first application confirmed the predicted performance of the propeller designed by the authors, improving its top speed by nearly two knots with respect to twin screw conventional shaft line arrangement; the increase in top speed is expected to be even more significant (up to 4 knots) if pod propulsors were fitted in the original Vshape stern hull form, which was on the contrary modified with semitunnels. The difference in bare hull resistance between the two stern arrangements has been predicted by viscous flow numerical solvers, with a model similar to those validated by Brizzolara and Villa [19].
Future developments foresee the enhancement of the method with new numerical procedure to determine the optimum load distribution on the two propellers based on a fully numerical lifting line method which might overcome