Acta Polytechnica doi:10.14311/AP.2017.57.0355 Acta Polytechnica 57(5):355–366, 2017 © Czech Technical University in Prague, 2017 available online at http://ojs.cvut.cz/ojs/index.php/ap CONTRA-ROTATING PROPELLER AERODYNAMICS SOLVED BY A 3D PANEL METHOD WITH COUPLED BOUNDARY LAYER Vít Štorch∗, Jiří Nožička Department of Fluid Dynamics and Thermodynamics, CTU in Prague, Technická 4, Prague, Czech Republic ∗ corresponding author: vit.storch@fs.cvut.cz Abstract. The aerodynamics of contra-rotating propellers is a complex three-dimensional problem of an unsteady flow, which is often approached by assuming numerous simplifications. Presented computational model combines a 3D panel method with a force-free vortex wake and a two-dimensional two-equation boundary layer model, in an attempt to capture all the main contributing elements of the flow physics. An emphasis is placed on the interaction of the viscous boundary layer region with the inviscid region and the development of a portable method of their coupling. The kinematics of a force-free vortex wake is supplemented with a vortex aging due to a diffusion. Extra attention is paid to the process of the blade passing through the wake of another blade. To demonstrate the ability of the numerical model, several test cases are presented illustrating the reaction of the system to various operational conditions. Keywords: 3D panel method, boundary layer, force-free wake, contra-rotating propeller, strong interaction. 1. Introduction Numerical tools based on inviscid flow methods de- veloped for specific types of open rotor aerodynamics are among the top choices for design and analysis purposes. The numerical model described in this pa- per has been developed for analysis of contra-rotating propellers but it is general enough to be applicable to a wide range of engineering problems of the external flow past lifting bodies, such as wing configurations, vertical and horizontal axis wind turbines and single propellers. Contra-rotating propellers (CRP) were subject of intensive research in the past but due to the mechan- ical and design complexity of the system, its usage is limited (e.g. coaxial helicopters and ship propeller systems). Nowadays, contra-rotating propellers are gaining new attention due to emerging fields of appli- cation in various unmanned aerial vehicles and also, for example, as a means of reducing an airliners tur- boprop engine fuel consumption. Existing numerical models such as that of Leishman and Ananthan [1] or those described in the review by Coleman [2], often approach the problem using the blade element momentum theory with Prandtl tip correction, or a lifting surface model with a vortex wake, neglecting the boundary layer and propeller thickness. The presented method aims to include as many aspects of the real CRP flow as possible. The main component of the model is an unsteady 3D panel method emitting a force-free wake. Viscous effects are assumed to be contained within a thin boundary layer region which is simulated by a two equation boundary layer model in an integral form. The 2D boundary layer model chosen for the imple- mentation and modification is the one used in the airfoil analysis tool XFOIL, well documented by Drela and Giles [3] and Drela [4]. It was chosen for its proven and reliable results even in the presence of a strong viscous-inviscid interaction, compared to one-equation models such as that of Thwaites. A strong interaction between the inviscid flow region and boundary layer occurs in the regions of transition, separation and laminar bubbles especially in low Reynolds number flows. Aim of the presented work is to develop and verify a CRP design and analysis tool which can provide reliable performance prediction for a wide range of operating conditions. To demonstrate the capabilities of the numerical model, it is being tested and verified on cases of single and contra-rotating propellers. 2. 3D panel method 3D panel method serves as a tool for solving Laplace’s elliptic partial differential equation: ∇2Φ = ∆Φ = ∂2Φ ∂x2 + ∂2Φ ∂y2 + ∂2Φ ∂z2 = 0. (1) Several elementary solutions, which can conve- niently be used to construct a new solution using superposition satisfying arbitrary boundary condition, exist. Suitable elementary solutions in three dimen- sions include a vortex filament, a point source and a point doublet and their various distributions over surfaces. Chosen elementary solutions used in the presented method are: flat constant source panels and quadrilateral vortex rings. Since any combina- tion of elementary solutions is a valid solution of the Laplace’s equation, the problem is reduced to find- ing the right combination of elementary solutions to satisfy the boundary condition. 355 http://dx.doi.org/10.14311/AP.2017.57.0355 http://ojs.cvut.cz/ojs/index.php/ap Vít Štorch, Jiří Nožička Acta Polytechnica The surface of each blade is discretized by a struc- tured grid of quadrilateral panels containing both constant source and vortex ring elementary solutions. These panels can generally be twisted, which will re- sult in small gaps between neighbouring flat constant source panels described by Hess [5]. The vortex ring panels take the twisted panel geometry into account. Velocity formulation of the 3D panel method is used. The boundary condition on the surface of the body takes form of a zero normal flow: ~n ·~c = 0. (2) Te boundary condition is enforced at collocation points, which are placed in centre point of each panel. For each panel i, a single boundary condition is written as ~ni · [ ~ci,rel + ~ci,ind + N∑ j ~bijΓj + N∑ j ~aijσj ] = 0, (3) where ~ni is the i-th panel normal, ~ci,rel is the relative velocity of the i-th panel collocation point to the free stream flow, ~ci,ind is the induced velocity at i- th collocation point by all vortex wakes present in the domain, ~aij is the induced velocity of j-th source panel with a unit source strength on i-th collocation point, ~bij is the induced velocity of j-th vortex ring panel with unit circulation on i-th collocation point, and finally Γj and σj are the unknown j-th panel circulation and source strength. To reduce the number of unknowns, the source strength σj of each panel is computed from the relative velocity according to Ashby [6]: σj = −~nj ·~cj,rel. (4) The Kutta condition, which requires a zero net cir- culation at the trailing edge, is met by attaching a row of wake panels to the trailing edge with the following circulation, calculated based on the circulations of upper and lower vortex ring panels near the trailing edge: Γwake = −Γupper + Γlower. (5) During each time step, the solution to the system of N linear equations (eq. 3) has to be found, where N is the number of collocation points in the domain and also the number of unknown circulation values. 3. Unsteady force-free vortex wake The wake is formed by vortex filaments, which are emitted from the trailing edge of each blade. During each time step, the wake is convected downstream under an influence of the induced and free stream velocity. Such model of wake produces a wake sheet aligned with the local velocity, which results in no forces acting on the imaginary wake surface (force- free wake). Each wake panel, represented solely by a vortex ring, receives a circulation value based on the Kutta condition upon its creation behind the trailing edge. The circulation is retained by the wake panel as it is convected downstream until its extinction at a predefined age. As a result, the wake contains the history of the blade bound circulation. To prevent unphysical induced velocities close to the wake surface and to enable a simulating diffusion, Lamb-Oseen 2D vortex core model [7] modified for use with the 3D vortex filament, is used in the following form: ~c = Γ 4π ~r1 ×~r2 |~r1 ×~r2|2 ~r0 · ( ~r1 |~r1| − ~r2 |~r2| ) · ( 1 − exp ( −1.2526 |~r1 ×~r2|2 R2C|~r1|2 )) , (6) where ~r1 and ~r2 are position vectors of a point in space relative to the starting point (1) and ending point (2) of the vortex filament. Vector ~r0 is a vector from point (1) to point (2). Diffusion process is simulated by a wake growth model described by Saffman [8], which simulates the effect of a dissipation over time and is introduced by the core growth model: RC = √ 4νt. (7) The procedure of forming new wake row is shown in Fig. 1. During the transition between two time steps, each blade rotates by the appropriate angle increment. Wakes stay at the same position which results in their detachment from the trailing edges. Then, the wake nodes are convected in the direction of a local velocity calculated as the sum of the free stream velocity and body induced and wake induced velocities. The last task before the transition to a new time step is completed, consists of filling the gap between wakes and blades with a new row of wake panels which will obtain circulation from the Kutta condition. 4. 2D boundary layer A two-dimensional boundary layer model using two governing equations in integral form as described by Drela and Giles [3], and Drela [4], was chosen as the basis of the boundary layer model used in the present computational method. Details of the boundary layer model and an example of an alternative way of cou- pling to a 2D panel method is discussed in a previous work of the authors [9]. The pair of governing non-linear first order ordinary differential equations of the 2D boundary layer model is formed by the integral momentum equation and kinetic energy shape parameter equation: dθ dξ + (2 + H) θ ue due dξ = Cf 2 ; (8) θ dH∗ dξ + H∗(1 −H) θ ue due dξ = 2CD −H∗ Cf 8 . (9) 356 vol. 57 no. 5/2017 Contra-Rotating Propeller Aerodynamics Figure 1. Kinematics of a force-free wake. H and θ were chosen as primary variables and the remaining variables are defined as functions of the primary variables H∗ = H∗(H,θ), Cf = Cf (H,θ) and CD = CD(H,θ). These functions are defined by closure equations that were implemented in the form found in [3]. Edge velocity ue is treated as a function of displacement thickness δ∗ and inviscid edge velocity ue inv. Primary attention is paid to this relation in the upcoming sections. An auxiliary equation is solved together with gov- erning equations. In the laminar region, auxiliary equation is the formula for detecting the transition using perturbation amplification ratio eñ based on the Orr-Sommerfeld equation. It describes the growth of Tollmien Schlichting waves based on the state of the boundary layer: dñ dξ = dñ dReθ m + 1 2 l 1 θ , (10) where: dñ dReθ = 1 100 (( 2.4H − 3.7 + 2.5 tanh(1.5H − 4.65) )2 + 0.25)1/2; (11) l = 6.54H − 14.07 H2 ; (12) m = ( 0.058 (H − 4)2 H − 1 − 0.068 )1 l . (13) When ñ reaches a predefined critical value ñcrit the solver switches to turbulent closure equations. The eñ transition model is often used for a transition prediction in low Reynolds number airfoil flow [3, 4]. It is assumed that both propellers are subjected to a low turbulence external flow and the downstream one passes the wake layer for only small fraction of the revolution. By changing the value of ñ, the free stream turbulence intensity can be accounted for. In the turbulent region, the auxiliary equation for the amplification ratio is replaced by the lag- entrainment equation for the shear stress coefficient. Green et al. [10] proposed a lag-entrainment method which estimates the shear stress coefficient based on its equilibrium value: θ(3.15 + H + 1.72 H−1 ) Cτ dCτ dξ = 4.2(C0.5τEQ −C 0.5 τ ), (14) where: CτEQ = H∗ 0.015(H − 1)3 (1 −Us)H3 ; (15) Us = H∗ 2 ( 1 − 4 3 H − 1 H ) . (16) To solve the boundary layer equations, a simple downstream marching algorithm is used. The gov- erning equations are discretized using a two-point central differencing scheme and solved by the Newton iteration method, station after station, during the downstream pass of the boundary layer. The single auxiliary and two governing equations are rewritten to produce zero on the right hand side. Left hand sides of the respective equations are labelled as functions F1,F2 and F3. Jacobian matrix J is shown for the case of the turbulent closure, where the square root of the shear stress coefficient is solved for together with H and θ: J =   ∂F1 ∂Hi ∂F1 ∂θi ∂F1 ∂C0.5 τ,i ∂F2 ∂Hi ∂F2 ∂θi ∂F2 ∂C0.5 τ,i ∂F3 ∂Hi ∂F3 ∂θi ∂F3 ∂C0.5 τ,i   . (17) After initializing boundary layer variables to sensi- ble values borrowed from the upstream station, the Newton iteration at i-th station proceeds as follows: J   Hi,new −Hi,oldθi,new −θi,old C0.5τ,i,new −C 0.5 τ,i,old   =  −F1−F2 −F3   . (18) 4.1. Viscous drag determination The viscous drag can be determined according to the wake layer properties far downstream [11]: Dvisc = %θ∞c2∞. (19) To prevent the need of tracking streamlines on the surface of the three-dimensional wake and computing the wake shear layer, an approximation was used in form of the Squire-Young formula, which extrapolates the behaviour of the wake in the far field based on the state of the boundary layer just behind the trailing edge. In the current implementation, the viscous drag is calculated as [11]: Dvisc = %θT.E. ( cT.E. c∞ )0.5(5+HT.E.) c2∞, (20) where cT.E. is the inviscid velocity just behind the trailing edge and HT.E. is the shape factor at the trailing edge. 357 Vít Štorch, Jiří Nožička Acta Polytechnica 4.2. 2D Boundary layer viscous-inviscid coupling The described two equation boundary layer model fails to produce a solution in most circumstances if a prescribed edge velocity distribution is used. This is due to a singular behaviour of the governing equations described by Goldstein for conditions near a separa- tion [12]. To obtain a converged solution, interaction between the boundary layer and inviscid solution must be provided in form of a relation between the edge velocity and displacement thickness. A two dimensional airfoil analysis tool XFOIL by M. Drela [4] uses simultaneous solutions of boundary layer equations and 2D panel method equations in one large system. Extending this approach to a three- dimensional body is associated with many difficulties. The boundary layer coupling implemented in the presented model uses a different, quasi-simultaneous, solution approach. This topic is covered by Veld- man [13], who discussed the criteria for convergence and existence of the solution and described a simple interaction law [14]: ue i,NEW = ue i,OLD + 4c∞ π∆ξ (δ∗i,NEW − δ ∗ i,OLD). (21) The first step that was performed in the search of a suitable viscous - inviscid interaction model was an exact determination of linear coefficients dij describing the change of edge velocity at i-th node due to the change of the displacement thickness at j-th node: dij = ∂ue i ∂δ∗j . (22) These linear coefficients form a response matrix, which can be used for a precise determination of the new velocity distribution, when the shape of the air- foil changes (with assumption of small displacements of the airfoil surface). An example of such response matrix is presented in Fig. 2. It was calculated using a numerical differentiation for NACA 0012 airfoil at 5◦ angle of attack, discretized by 17 surface nodes. Im- portant feature of the matrix is a dominating positive main diagonal with still significant negative subdiago- nal and superdiagonal. Other members of the matrix have less significant values. Each coefficient on the main diagonal dii shows the response of the edge ve- locity at i-th surface node to its own displacement in the normal direction. Therefore the coefficients from the main diagonal can be used in an interaction law: ue i,NEW = ue i,OLD + dii(δ∗i,NEW − δ ∗ i,OLD). (23) The interaction law shown above, which is used in the present computational model, is a result of an effort to produce an accurate local linear interaction coefficient which increases the convergence rate of the solver. One method that was initially tested for calculating the linearized interaction coefficient dii was using a 8.9 −3.1 −0.2 −0.1 0 0 0 0 0 0 0 0 0 −0.1 −0.1 −3 9 −5.2 12 −5.2 0.2 −0.3 −0.1 −0.1 0 0 0 0 0 0 −0.2 −0.4 0.1 −0.6 −1.3 −3.9 12.1 −5.8 0.3 −0.4 −0.1 −0.1 0 0 0 0 −0.1 −0.3 0.6 −1.2 0.7 −1.4 1 −5.7 14.2 −7.2 0.5 −0.5 −0.1 −0.1 0 0 −0.1 −0.2 0.3 −0.2 −1 0.7 −1.3 0.8 0 −6.9 17.6 −9.3 0.8 −0.8 −0.1 0 −0.1 −0.1 0.3 −0.2 −0.1 −1 0.9 −1.5 1.2 −0.4 0 −9.3 23.5 −13.5 1.6 −0.8 0 −0.1 0.3 −0.1 −0.2 0 −1.3 1.1 −1.9 1.7 −0.2 −0.5 −0.3 −13.9 35 −21.9 2.1 −0.2 0.4 0.1 −0.1 −0.1 0 −1.7 1.6 −3.1 3.1 −0.2 −0.2 −1 −1.1 −26 53.1 −23.2 0.4 0.5 0.1 0.1 −0.1 0 −2.9 2.8 −6.3 6.7 −0.1 0.1 −0.8 −2.1 −3.1 −49.4 16.5 5.6 1.8 1.4 0.6 −0.1 0.1 −6.2 6.3 2.8 −3.2 0 −0.2 0.1 0.1 0.6 −1.6 12.6 14 −13.9 −0.7 −0.8 −0.1 −0.2 2.8 −3.1 1.6 −1.9 0 −0.1 −0.1 0 0.9 1.1 −0.8 −5.8 19.2 −9.5 −0.2 −0.4 −0.2 1.6 −1.9 1.1 −1.4 0 −0.2 −0.1 0.5 0 −0.1 0.5 0.4 −7.5 16.2 −7.2 0 −0.4 1.1 −1.5 0.9 −1.1 −0.1 −0.2 0.4 −0.1 −0.1 0.2 0.1 −0.2 0.4 −6.5 13.6 −5.8 0 0.7 −1.3 0.7 −1.1 −0.2 0.4 −0.2 −0.2 0.1 0 0.1 0 −0.3 0.3 −5.6 11.8 −5 0.9 −1.4 0.7 −1.3 0.6 −0.3 −0.2 0 0 0 0 0 −0.1 −0.3 0.2 −4.9 10.6 −3.6 −1.3 −0.7 0.2 −0.4 −0.3 0 0 0 0 0 0 0 −0.1 −0.2 0.1 −4.6 11.1 −5.2 8.9 −3.1 −0.2 −0.1 0 0 0 0 0 0 0 0 0 −0.1 −0.1 −3 9 i=1 2 3 ... ... ... ... ... ... ... ... ... ... ... ... ... N d ij j=1 2 3 ... ... ... ... ... ... ... ... ... ... ... ... ... N Figure 2. Typical matrix of velocity response to displacement thickness numerical differentiation. This method is computa- tionally expensive, since the interaction coefficient is determined based on the edge velocity of the origi- nal geometry and of a modified geometry, where the i-th point has been moved outwards in the normal direction by a small displacement value ∆δ∗i : dii = ∂ue i ∂δ∗i ≈ umode i −ue i ∆δ∗i . (24) A fast way of approximating the interaction coeffi- cient dii based on only one inviscid solution was found to be as follows: dii = 2ue i,inv (ξi − ξi−1) . (25) 4.3. Portable boundary layer model The local linear interaction coefficient dii only esti- mates the local response of the edge velocity to a local change of the displacement thickness. Further bound- ary layer passes, with inviscid flow solution updates in-between, are necessary to arrive at a converged solution. Even with an accurate linear interaction coefficient, between 10 and 20 passes are required for a converged solution. This is a known problem of quasi-simultaneous viscous-inviscid interaction meth- ods [14]. Repeatedly calculating an inviscid solution between passes can be easily performed with a two-dimensional panel method, but this approach would be too time consuming for a three dimensional panel method with a sufficient panel density. A completely different solution approach was devel- oped for the boundary layer coupling to a 3D panel method. To make the boundary layer model easily portable, a replacement inviscid model, which allows to approximate the edge velocity based on initial sur- face velocity distribution and displacement thickness alone, is proposed. 4.4. Replacement inviscid model Each column of the response matrix dij contains in- formation about the response of the edge velocity to 358 vol. 57 no. 5/2017 Contra-Rotating Propeller Aerodynamics 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.1 −0.05 0 0.05 0.1 Detail ∆ δ j * Detail 2:1 c ∞ x [m] y [m ] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Detail Detail 2:1 ∆ u e i −1∆ u e i ∆ u e i +1 Velocity distribution of NACA 0012, α = 5°, c ∞ = 1 m/s with and without surface jump x [m] u e [ m /s ] original airfoil airfoil with surface jump Figure 3. Response of velocity to surface jump. Note: Surface displacement has been exaggerated for illus- tration purposes a single node displacement. One such response is demonstrated in Fig. 3. As can be seen, neighbouring nodes of the node j with surface jump ∆δ∗j are subjected to a drop in the edge velocity of about half the magnitude of the ve- locity growth at j-th node. This can also be observed by comparing the values of the subdiagonal and su- perdiagonal with the main diagonal of the response matrix. Based on these observations, a new replacement inviscid model, which can be used for boundary edge velocity updates between boundary layer passes, has been formulated: ue i = ue i,orig + dii(δ∗i − 0.5δ ∗ i−1 − 0.5δ ∗ i+1) (26) Although the replacement inviscid model does not account for global effects of boundary layer thickening, such as shifting of the stagnation point, its estima- tion of the edge velocity behaviour is sufficient to produce a boundary layer solution surprisingly sim- ilar to the XFOIL results (see Fig. 4). When the replacement inviscid model is used together with the boundary layer model, the only input parameters needed are the velocity distribution with surface co- ordinates. This makes the resulting boundary layer truly portable, with a small impact on the result ac- curacy. 5. Time-stepping calculation procedure The process begins with a sudden start of the rotation of blades without any wake. The wakes gain one addi- tional row of panels during each time step. Each blade can have assigned a different rotational frequency f 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 2 3 4 5 6 7 8 9 x 10 −3 x/b [1] δ* /b [ 1 ] Suction side B.L. displacement thickness, NACA 0012, α=5°, Re=5e5 XFOIL Quasi−simultaneous B.L. with repeated 2D panel method updates Portable boundary layer with replacement inviscid model Figure 4. Displacement thickness distributions on the suction side of NACA 0012 airfoil calculated using different methods and even different axis of rotation. During the i-th iteration step, the following subsequent operations are performed (1.) All blades are rotated to their new positions de- fined by the position angle φk i = φk i−12πfk∆t. At this new position, the relative velocity due to the rotation is calculated at each collocation point on the blade surface. (2.) Wakes become detached from the blades due to the rotation in the previous step. Each wake is con- vected downstream and attached to the appropriate blade via a new panel row as described in section 3. (3.) The velocity induced by each wake on itself, other wakes and all blades is calculated in this step using an appropriate vortex core size given by the growth model. (4.) Next, the system of equations given by a 3D panel method (section 2) is solved. (5.) Once the distribution of the vortex ring circula- tion and source panel’s strengths are known, the velocities induced by the blades on the wake nodes are calculated (6.) When the portable boundary layer model is en- abled, several passes of the boundary layer solver along streamlines are performed (see Fig. 5). The resulting displacement thickness is used to modify the original shape of the blade. Second calculation of the 3D panel method is carried out using newly modified geometry. (7.) After evaluating performance parameters of each blade, time step (i + 1) can begin. 6. Practical considerations Aerodynamics drives the design of a blade only from a certain radius towards the blade tips. The root section of a propeller serves for the attachement of the blade to the shaft and is shaped with regards to 359 Vít Štorch, Jiří Nožička Acta Polytechnica Figure 5. Surface streamlines of a heavily loaded blade. stiffness and strength. To avoid numerical difficulties in the region of a blunt trailing edge, a replacement root section is modelled for each blade. A procedure that can be described as a wake blow- off is used for low advance ratios, where the free stream velocity is not sufficient to convect the newly created wake downstream away from the propeller. When the wake stays in the plane of the propeller disc, it loses stability and collapses into a vortex cloud (Fig. 8a). To help producing a stable solution and speed up the time for establishing a converged wake shape, the wake is blown off at the beginning of the time stepping (Fig. 8b) using some predefined advance ratio, which is gradually lowered towards the required low value during the time stepping. On the one hand, during the passage of the blade through a wake sheet, the influence of wake vortex filaments upon blade in terms of an induced velocity is without a presence of unrealistic velocities thanks to the growing vortex core immediately behind the trailing edge. On the other hand, vortex wake nodes in a close proximity of blade surface panels can ex- perience unphysical velocities of an unbound magni- tude induced by the surface panels. These unreal- istic velocities occur due to the discrete modelling of surface vorticity and can destroy the wake shape completely. To solve this problem, blade surface vor- tex filaments feature a finite viscous core model as well. The core does not grow and its core radius is determined according to the local panel size, small enough to prevent influencing the 3D panel method inviscid solution. To decrease the computational time, the wake length is limited by a maximum age of each wake filament. A typical calculation takes 200 time steps, while the wake age is limited to 100 time steps. Be- cause the boundary layer calculation takes significant time, only last 40 time steps, used for performance evaluation are calculated with boundary layer enabled. Previous time steps, during which the wake shape develops, are calculated as inviscid. These two mea- sures alone can decrease computational time more than five times while having a negligible effect on results. a) Wake blow off disabled, λ = 0 b) Wake blow off enabled, λ = 0 Figure 8. Wake behind single propeller. 7. Performance parameters of single and contra-rotating propellers For performance comparisons, non-dimensional values of the advance ratio λ, thrust coefficient cT , and power coefficient cP of a single propeller operating under uniform free stream velocity c∞ can be defined as follows [15]: λ = c∞ πfD ; (27) cT = T ρf2D4 ; (28) cP = P ρf3D5 ; (29) η = Tc∞ P = λπ cT cP . (30) Contra rotating propellers used for the experimen- tal verification were measured [16] only in static thrust conditions, where efficiency is always zero. In a propulsion system, the efficiency of producing a thrust in static conditions is expressed by the Figure of Merit “FoM” which is defined, for a single propeller, 360 vol. 57 no. 5/2017 Contra-Rotating Propeller Aerodynamics Fine mesh size 38x29 Medium mesh size 22x19 Coarse mesh size 14x13 Figure 9. Three meshes used in mesh sensitivity study 0 0.1 0.2 0.3 0.4 0 0.05 0.1 c T [ 1 ] λ [1] Thrust sensitivity to mesh, propeller 22x18 0 0.1 0.2 0.3 0.4 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 c P [ 1 ] λ [1] Power sensitivity to mesh, propeller 22x18 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 η [ 1 ] λ [1] Efficiency sensitivity to mesh, propeller 22x18 Coarse mesh Medium mesh Fine mesh Coarse mesh Medium mesh Fine mesh Coarse mesh Medium mesh Fine mesh Figure 10. Sensitivity of propeller performance curves to mesh by Eq. 31: FoM = T 3/2 P √ 2ρA (31) The Figure of Merit definition for contra-rotating propeller systems is discussed in detail by Leishman and Anathan in [1]. The authors conclude that “fun- damentally, any definition of FoM can be adopted, as long as the definition is used to compare like-with- like” and presents a very complex FoM definition that takes into account the relative thrust sharing of the two rotors and unequal disk loadings. For the purpose of this work, the following simpler definition of the FoM is used (index 1 marks the upstream propeller, index 2 downstream propeller): FoM = (T1 + T2)3/2 (P1 + P2) √ 2ρ max(A1,A2) (32) Equations 33,34 and 35 describe dimensionless pa- rameters that are used for evaluating the overall per- formance of the contra-rotating propeller system. cT = T1 + T2 ρ0.25(f21 + f22 )(D41 + D42 ) ; (33) cP = P1 + P2 ρ0.25(f31 + f32 )(D51 + D52 ) ; (34) η = (T1 + T2)c∞ P1 + P2 . (35) Blade mesh Blade + wake Average panels time/step 14×13 2444 2.1 s 22×19 3876 4.7 s 38×29 6844 15.12 s Table 1. Computational cost vs. mesh size. 8. Mesh sensitivity A mesh sensitivity study was performed for the case of a single propeller because the absence of the second propeller allowed computing a finer mesh without an excessive computational time. Each blade of the single propeller was discretized by 14 × 13 (coarse) 22 × 19 (medium) and 38 × 29 (fine) mesh panels (see Fig. 9). Upstream propeller from the contra-rotating propeller set was used for this study. Based on the results (Fig. 10 and Tab. 1), a medium mesh was selected for a further blade geometry representation. 9. Experimantal verification results The experimental verification of the numerical method was performed using data from a previous experimen- tal research of propellers [16, 17]. Contra rotating propellers were measured only in static thrust condi- 361 Vít Štorch, Jiří Nožička Acta Polytechnica 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.02 0.04 0.06 0.08 0.1 0.12 c T [ 1 ] λ [1] Single propeller APC 20x13 thrust validation 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.01 0.02 0.03 0.04 0.05 0.06 c P [ 1 ] λ [1] Single propeller APC 20x13 power validation 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 η [ 1 ] λ [1] Single propeller APC 20x13 efficiency validation inviscid 3D panel method viscous (3D p.m. + B.L.) experiment inviscid 3D panel method viscous (3D p.m. + B.L.) experiment inviscid 3D panel method viscous (3D p.m. + B.L.) experiment Figure 11. Comparison of numerical and experimental data in case of a single propeller tions. For a verification of the numerical model under forward flight conditions, results of a single propeller in a wind tunnel were used. 9.1. Verification using single propeller wind tunnel results Based on the results obtained during recent series of measurements of model propellers in a wind tun- nel [17], data from the measurement of an APC 20×13 propeller were selected for the verification of the nu- merical model in a single propeller test scenario. Due to vibrations, the measured data from the wind tunnel contains some scattering of points, as can be seen from the error bars (Fig. 11). Geometry of the propeller blades was scanned using a manual contact scanner. The coefficient of thrust is predicted well, except a region near the static thrust (λ = 0), where the numer- ical results are too optimistic. This is most probably connected to a degraded accuracy of the simplified coupling between the panel method and boundary layer where a significant trailing edge separation is already present. Overpredicted thrust also results in an overprediction of the power coefficient near the static thrust. It appears that the inviscid calcula- tion power prediction is more accurate near the static thrust. This is however an accidental match, since the power due to the thrust overprediction is compen- sated by missing viscous power losses. The inviscid calculation was performed using the 3D panel method with the boundary layer disabled. In the region of the peak efficiency (λ = 0.1−0.2), the panel method with the boundary layer a provides best match with the experiment. 9.2. Verification using measured contra-rotating propellers in static thrust regime To test the performance of the computational method, the data from the static thrust measurement of a pair of 22-inch carbon fibre propellers were used for 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure of Merit, propeller distance 7% rotation frequency ratio f 1 /(f 1 +f 2 ) [1] F o M [ 1 ] FoM − experiment FoM − inviscid calculation FoM − 3D p.m. + B.L. Figure 12. Figure of Merit vs. ratio of frequency of rotation a comparison. The upstream propeller has a des- ignation 22 × 18 and is marked by an index 1 and the downstream propeller’s size is 22×20 and is de- noted by an index 2. The measurement described in [16] provides the thrust, torque and rpm of each propeller for various rotational speed ratios. Al- though the pair of propellers was designed to give a peak performance for 1:1 speed ratio at the de- sign flight speed, at a static thrust regime, the mea- sured Figure of Merit showed a peak performance when the upstream propeller was rotating consider- ably faster. Analysis of numerical results showed some degree of a trailing edge separation on both propellers even at the peak of Figure of Merit. At extreme speed ratios, where one of the propellers rotates more than twice faster than the other, the faster pro- peller tends to suffer a large separation on most of the blade span. This is supported by the ex- perimental results, where the Figure of Merit falls rapidly towards the extreme speed ratios. Consid- 362 vol. 57 no. 5/2017 Contra-Rotating Propeller Aerodynamics 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 rotation frequency ratio f 1 /(f 1 +f 2 ) [1] c T , c P [ 1 ] Total thrust and Power coefficients, propeller distance 7% c T − experiment c P − experiment c T − inviscid calculation c P − inviscid calculation c T − 3D p.m. + B.L. c P − 3D p.m. + B.L. Figure 13. Total thrust and power coefficients of the system vs. ratio of frequency of rotation 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0 0.05 0.1 0.15 0.2 0.25 Thrust coefficient, propeller distance 7% rotation frequency ratio f 1 /(f 1 +f 2 ) [1] c T [ 1 ] c T1 − experiment c T2 − experiment c T1 − inviscid calculation c T2 − inviscid calculation c T1 − 3D p.m. + B.L. c T2 − 3D p.m. + B.L. Figure 14. Thrust coefficient of each propeller vs. ratio of frequency of rotation ered that the separation and stall are still among very problematic phenomena for the simulation, the experimental and numerical results of the Figure of Merit and total thrust and torque (Figures 12 and 13) match quite nicely. When the perfor- mance of each propeller is analysed separately (Fig- ures 14 and 15), the match between the experiment and numerical results is much weaker. This can be attributed to the interference between the mea- sured pair of thrusts and torques due to a friction between coaxial shafts, because the global perfor- mance of the contra-rotating system is not influ- enced. 10. Contra-rotating propeller numerical results 10.1. Effects of propeller distance The simple momentum theory predicts a subtle in- crease of performance when the propeller distance is increased. This is due to additional air being sucked 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 Power coefficient, propeller distance 7% rotation frequency ratio f 1 /(f 1 +f 2 ) [1] c P [ 1 ] c P1 − experiment c P2 − experiment c P1 − inviscid calculation c P2 − inviscid calculation c P1 − 3D p.m. + B.L. c P2 − 3D p.m. + B.L. Figure 15. Power coefficient of each propeller vs. ratio of frequency of rotation 0 10 20 30 40 50 60 0.6 0.65 0.7 0.75 0.8 0.85 Figure of Merit vs. propeller distance relative propeller distance d/D [%] F o M [ 1 ] FoM − 3D p.m. + B.L. FoM − experiment Figure 16. Figure of Merit vs. propeller distance in between the two propellers, thus decreasing the disc loading required for a given thrust. Experimental results show a little change to the total thrust and power with increasing propeller distance, but most sources agree (e.g. [2, 18]), that the propeller distance influences the thrust and torque ratio between the two propellers noticeably. With the increasing propeller distance, the upstream propeller becomes less and less influenced by the downstream propeller induced velocity, which results in a higher loading of the up- stream rotor. The higher loaded upstream propeller produces slightly stronger slipstream with a higher velocity, which, in turn, unloads the downstream pro- peller. Numerical study was performed with the same set of propellers as used in the previous experimental verification. Both propellers were rotating at the same speed under static thrust conditions (i.e. hover) and their distance d varied between 7% and 60% of the propellers diameter. There are only two data points for two propeller distances available from the experi- mental investigation and these were included in the graphs. Figure 16 shows a small rise of the Figure of Merit with a propeller distance as supported by the momentum theory. Notable thrust redistribution 363 Vít Štorch, Jiří Nožička Acta Polytechnica 0 10 20 30 40 50 60 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 Thrust coefficient of each propeller vs. propeller distance relative propeller distance d/D [%] c T [ 1 ] c T1 − 3D p.m. + B.L. c T2 − 3D p.m. + B.L. c T1 − experiment c T2 − experiment Figure 17. Thrust coefficient of each propeller vs. propeller distance showing thrust sharing Figure 18. Wake behind CRP with axial inflow between propellers is visible in Fig. 17 for propeller distances below 20%D. For higher propeller distances, the thrust of the upstream propeller remains constant, while the downstream propeller thrust slowly rises, probably due to dissipation effects of the upstream propeller wake. 10.2. Transition between horizontal flight and hover The transition between the horizontal flight and hover is a routine manoeuvre performed by many unmann ed aerial vehicles (UAVs). The three-dimensional un- steady character of the presented computational model allows simulating these complex conditions, where models with axi-symmetric flow assumption are not usable. The presented simulation starts with a contra- rotating propeller set, used in previous cases, rotating at the same speed under an axial inflow of the free stream velocity producing a low advance coefficient λ = 0.06. The free stream flow velocity component in the direction perpendicular to the axis of the rotation is increased from zero to twice the axial component. The angle of inflow β therefore changes from 0° to 63°. Figure 19. Wake behind CRP with 60° inflow 0 10 20 30 40 50 60 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 Figure of Merit and efficiency vs. inflow angle Free stream inflow angle β [°] F o M , η [ 1 ] FoM − 3D p.m. + B.L. η − 3D p.m. + B.L. Figure 20. Figure of Merit and efficiency vs. inflow angle The shape of the wake forming behind the propellers is shown in Figures 18 and 19. As can be seen from the results, the direction of the free stream velocity has a low impact on thrust coefficients. All parameters seem to be fairly insensi- tive to the inflow angle. Due to a very low advance ratio, both the Figure of Merit and efficiency were used for evaluation. Thrust coefficients, Figure of Meri and efficiency remain almost constant for inflow angles below 30° and gradually begin rising for higher inflow angles (see Figures 20 and 21) It is generally known that helicopters and other similar rotorcraft of- ten consume significantly less power for forward flight than hovering, which is in accordance with the results above. 11. Conclusions Presented numerical model is capable of a fast analy- sis of contra-rotating propellers respecting the mutual influence of the individual blades. Apart from pro- viding average thrusts and torques of each blade, the computational method gives an insight into the char- acter of unsteady loads acting on the blades, calculates complex wake shapes and provides information about 364 vol. 57 no. 5/2017 Contra-Rotating Propeller Aerodynamics 0 10 20 30 40 50 60 0.1 0.105 0.11 0.115 0.12 0.125 0.13 0.135 0.14 0.145 Thrust coefficients vs. inflow angle Free stream inflow angle β [°] c T [ 1 ] c T1 − 3D p.m. + B.L. c T2 − 3D p.m. + B.L. Figure 21. Thrust coefficients vs. inflow angle viscous forces, transition position and possible sep- aration regions. The force-free wake component of the model performs well also in static thrust condi- tions, where many simpler models brake down. The method copes well (although with a reduced accu- racy) with flows past a maximum lift angle of attack, where a separation occurs. This is visible in the ex- perimental verification of both single propeller and contra-rotating setups. Due to its computational time requirements, the method fits between simple blade element momentum methods and finite volume solvers. It is ideal for fine tuning geometry of contra-rotating blades and analysis of their performance. Further im- provements can be made in the treatment of separated boundary layer regions to extend the accuracy of the calculation near the limits of the operational envelope. Thanks to the versatility of the numerical method, it can be applied with small changes to any case of wings and blades under an arbitrary motion. List of symbols A propeller disk area [m2] ~aij source panel influence vector [m] ~bij vortex ring panel influence vector [m−1] ~c flow velocity vector [m s−1] c∞ free stream velocity magnitude [m s−1] CD dissipation coefficient [1] cT coefficient of thrust [1] cP coefficient of power [1] Cf skin friction coefficient [1] Cτ shear stress coefficient [1] CτEQ equilibrium shear stress coefficient [1] D propeller diameter [m] Dvisc viscous drag force [N] d distance between propellers [m] dii interaction coefficient [s−1] FoM Figure of Merit [1] f frequency of rotation [s−1] H shape parameter [1] H∗ kinetic energy shape p [1] L lift force [N] M number of spanwise panels [1] N number of streamwise panels [1] ~n normal vector [1] ñ TS wave amplification exponent [1] P shaft power [W] p0 Total pressure [Pa] Re Reynolds number [1] Reθ momentum thickness Reynolds number [1] RC radius of viscous core [m] ~r position vector [m] T thrust [N] t time [s] ue B.L. edge velocity [m s−1] x,y,z space coordinates [m] α angle of attack [°] β inflow angle [°] Γ circulation of a vortex filament [m2 s−1] δ∗ displacement thickness [m] η propeller efficiency [1] θ momentum thickness [m] λ advance ratio [1] ξ B.L. surface coordinate [m] ν kinematic viscosity [m2 s−1] % density [kg m−3] σ source strength density [s−1] Φ Potential [m2 s−1] Acknowledgements This work was supported by the Grant Agency of the Czech Technical University in Prague, grant SGS16/068/OHK2/1T/12 and the Centre of 3D volumet- ric anemometry — COLA, reg. no. CZ.2.16/3.1.00/21569 — Operating Program Prague competitiveness. References [1] J. G. Leishman, S. Ananthan. Aerodynamic optimization of a coaxial proprotor. In 62nd American Helicopter Society International Annual Forum. 2006. [2] C. P. Coleman. A survey of theoretical and experimental coaxial rotor aerodynamic research. Tech. Rep. 3675, NASA, 1997. [3] M. Drela, M. Giles. Viscous-inviscid analysis of transonic and low Reynolds number airfoils. AIAA Journal 25(10):1347–1355, 1987. doi:10.2514/3.9789. [4] M. Drela. XFOIL – An analysis and design system for low Reynolds number airfoils. Low Reynolds number aerodynamics 1989. doi:10.1007/978-3-642-84010-4_1. [5] A. Smith, J. Hess. Calculation of non-lifting potential flow about arbitrary three-dimensional bodies. Tech. Rep. E.S. 40622, Douglas Aircraft Company, Inc., Long Beach, CA, 1962. [6] D. Ashby. Potential flow theory and operation guide for the panel code PMARC_14. NASA, Ames Research Center, Moffet Field, CA, 1999. TM-1999-209582. [7] H. Lamb. Hydrodynamics, Fourth edition. Cambridge University Press, London, 1916. doi:10.5962/bhl.title.18729. 365 http://dx.doi.org/10.2514/3.9789 http://dx.doi.org/10.1007/978-3-642-84010-4_1 http://dx.doi.org/10.5962/bhl.title.18729 Vít Štorch, Jiří Nožička Acta Polytechnica [8] P. G. Saffman. Vortex Dynamics. Cambridge Monographs on Mechanics. Cambridge University Press, 1992. doi:10.1017/cbo9780511624063. [9] V. Štorch, J. Nožička. On viscous-inviscid interaction for boundary layer calculation using two-equation integral method. In Studentská tvůrčí činnost 2015. CTU in Prague, Faculty of Mechanical Engineering, 2015. [10] J. E. Green, D. Weeks, J. Brooman. Prediction of turbulent boundary layers and wakes in compressible flow by a lag-entrainment method. H.M.S.O, London, 1977. [11] R. Eppler. About classical problems of airfoil drag. Aerospace Science and Technology 7(4):289–297, 2003. doi:10.1016/s1270-9638(03)00029-4. [12] S. Goldstein. On laminar boundary-layer flow near a position of separation. The Quarterly Journal of Mechanics and Applied Mathematics 1(1):43–69, 1948. doi:10.1093/qjmam/1.1.43. [13] A. E. P. Veldman. New, quasi-simultaneous method to calculate interacting boundary layers. AIAA Journal 19(1):79–85, 1981. doi:10.2514/3.7748. [14] A. E. Veldman. A simple interaction law for viscous– inviscid interaction. Journal of Engineering Mathematics 65(4):367–383, 2009. doi:10.1007/s10665-009-9320-0. [15] Q. R. Wald. The aerodynamics of propellers. Progress in Aerospace Sciences 42:85–128, 2006. doi:10.1016/j.paerosci.2006.04.001. [16] V. Štorch, J. Nožička, M. Brada. Experimental setup for measurement of contra-rotating propellers. In Topical Problems of Fluid Mechanics 2017, pp. 285–294. edited by Šimurda D., Bodnár T., 2017. ISSN 2336-5781, doi:10.14311/tpfm.2017.036. [17] J. Filipský, V. Štorch. Comparison of propeller analysis methods and experimental data. In Engineering Mechanics 2014, pp. 172–175. Brno University of Technology, 2014. ISSN 1805-8248. [18] M. Ramasamy. Measurements comparing hover performance of single, coaxial, tandem and tilt-rotor configurations. In AHS 69th Annual Forum. 2013. 366 http://dx.doi.org/10.1017/cbo9780511624063 http://dx.doi.org/10.1016/s1270-9638(03)00029-4 http://dx.doi.org/10.1093/qjmam/1.1.43 http://dx.doi.org/10.2514/3.7748 http://dx.doi.org/10.1007/s10665-009-9320-0 http://dx.doi.org/10.1016/j.paerosci.2006.04.001 http://dx.doi.org/10.14311/tpfm.2017.036 Acta Polytechnica 57(5):355–366, 2017 1 Introduction 2 3D panel method 3 Unsteady force-free vortex wake 4 2D boundary layer 4.1 Viscous drag determination 4.2 2D Boundary layer viscous-inviscid coupling 4.3 Portable boundary layer model 4.4 Replacement inviscid model 5 Time-stepping calculation procedure 6 Practical considerations 7 Performance parameters of single and contra-rotating propellers 8 Mesh sensitivity 9 Experimantal verification results 9.1 Verification using single propeller wind tunnel results 9.2 Verification using measured contra-rotating propellers in static thrust regime 10 Contra-rotating propeller numerical results 10.1 Effects of propeller distance 10.2 Transition between horizontal flight and hover 11 Conclusions List of symbols Acknowledgements References