Strojniski vestnik - Journal of Mechanical Engineering 56(2010)6, 357-367 UDC 629.735.45 Paper received: 23.03.2007 Paper accepted: 02.03.2010 An Optimal Main Helicopter Rotor Projection Model Obtained by Viscous Effects and Unsteady Lift Simulation Caslav Mitrovic1* - Aleksandar Bengin1 - Dragan Cvetkovic2 - Dragoljub Bekric1 1 University of Belgrade, Faculty of Mechanical Engineering, Serbia 2 University "Singidunum" Belgrade, Faculty of Informatics and Management, Serbia This paper presents an airfoil model which provides a method for optimal main helicopter rotor projection by viscous effect and unsteady lift simulation through algorithm and set of program entireties, applicable to the ideological and main project of helicopter rotor. Based on real rotors theoretical consideration and the, numerical analysis considerations in this paper can be applied with sufficient accuracy in the analysis and constructive realizations of helicopter rotor in real conditions. The method for unsteady viscous flow simulations by inviscid techniques is developed. The aim of this paper is to determine helicopter rotor blade lift with the highest possible accuracy by using a singularity method and to define an optimal conception model of aerodynamic rotor projection corresponding to rotor behavior in real conditions and with sufficient quality from the aspect of engineer use. © 2010 Journal of Mechanical Engineering. All rights reserved. Keywords: unsteady lift, viscous effect simulation, potential flow, lifting surface, panel method, free-wake model, main helicopter rotor aerodynamics 0 INTRODUCTION The basic aim of this paper is to determine the lift of helicopter rotor blades by using the singularity method, and to develop an optimal model of aerodynamic rotor planning by simulating unsteady flow and viscous effects. The model is considered corresponding to rotor behavior in real conditions and with high quality from the aspect of engineer use. The idea of unsteady lift modeling and viscous effects simulation by using the singularity method is based on the need to avoid expensive experiments in the first stage of rotor planning by using contemporary aerodynamic analysis applied to the available computer technique. It is necessary to modulate vortex wake, unsteady 2D flow field characteristics and blade dynamic characteristics in order to determine the aerodynamic forces acting upon the helicopter rotor blade. First, the airfoil is approximated by vortex and source panels, and boundary layer by layer of vortices changing their position during a certain period of time. Vortex trail and the trail of separation are modulated by free vortices. The separation is modulated by free vortices. The dependence of lift coefficient on the angle of attack is determined for the known motion of an airfoil. The behavior of one separated vortex in an article time model is spread into a series of positions at a certain number at different vortices. The separated flow velocity profile is approximated by the superposition of displacement thickness of these vortices coupled with the potential solution model. After every time step the position of free vortices is changed, for which a generation of new vortices that would together satisfy the airfoil contour boundary conditions is required. Secondly, the model for 3D for complete flow field is developed and is applied to the calculation of flow over the helicopter blade. Helicopter rotor blade is divided into a certain number of segments. Each segment defines an airfoil whose flow conditions are similar to those described around an airfoil. The flow over the entire helicopter rotor blade can be defined by using certain intersection characteristics. This can form the basis for unsteady force acting upon the blade, the rotor and the helicopter study. In addition to the complex problem of analytical modeling of these phenomena there is also an issue of modeling interactions between the effects of these phenomena. This approach to helicopter rotor planning allows experimental investigations in tunnels and in flight to be final examinations. In that way, a very useful * Corr. Author's Address: Faculty of Mechanical Engineering, Kraljice Marije St. 16, 11000 Belgrade, Serbia, cmitrovic@mas.bg.ac.rs interaction between a numerical calculation and experimental results is achieved. In addition to the standard decompression and synthesis method for the purposes of this paper the singularity method, the panel method, the vortex surface method and the PIC (particle in cell) method will be used. The final result is the achievement of unsteady lift dependence of the angle of attack. A model established in such a way is characteristic of the helicopter rotor blade airflow and is based on the influence of the previous lifting surface's wake influence on the next coming blade. According to such an analysis and a chosen model, a computer program, used for the helicopter blade airflow analysis is developed. 1 FOUNDATIONS OF THE IRROTATIONAL 2-D FLOW The planar potential flow of incompressible fluid can be treated in Cartesian coordinates x and y. If the physical plane is mapped to the complex plane by z = x + iy where i = V-T. The symmetrical point with respect to the x-axis is z = x - iy . A two-dimensional potential incompressible flow is completely defined by the speed potential and stream function and is presented by Cauchy-Reimann Eqs. [1] and [2]: 9 6 9w „ j 9w 96 „ -!- + —!- = 0 and —--- = 0 , (1) 9 x 9 y 9 x 9 y where u and v are velocity in x and y direction respectively. The Laplace partial differential equations can also be introduced: V26 = 0 and V> = 0 , (2) where V is the Laplace's operator: ^ 92 92 V2 =-r +-r . (3) 9x 9y (3) The fulfillment of Cauchy-Reimann conditions enables combining the velocity potential and stream function: w(z ) = 6( x y) +i y( x, y) (4) of complex variable z = x + iy . This complex function entirely defines the planar potential flow of incompressible fluid as a function of a complex coordinate. The complex analytical function w (z), called the complex flow potential, always has a unique value for the first derivative. This derivative of the complex potential is equal to the complex velocity at that point, i.e.: dw dz = u - iv = V (5) where u is its real part (velocity component in x-direction) and v is the imaginary part (velocity component in ^-direction). Circulation and flow equal zero for any closed curve in a complex plane. Complex potential has no singularities except at the stagnation point. 2 SIMULATION OF THE MOVING VORTEX If we assume that the moving vortex of intensity I0 is at the distance z0, then the perturbance potential of such a flow is: w(z) = Vœze-i'a + Vœ — e+'a + i2aV • sin(oc-P)lnz + z ir„, f a2 _ +i^°ln(z - z0)ln I —-z( 2^ 2ïï I z (6) and its complex velocity: - a2 i V = V e-io-V — e+io +— 2a•V • sin(a-B) + œ œ 2 œ V ' ' z z 'T0(a 7 z2) (7) iTq 1 _ + (z-z0) + zQ a _ In this case a displacement of the stagnation point appears, which is now at the distance: z = a • eiP . (8) According to the request that complex velocity at the stagnation point must be zero: = V = 0, dz (9) the position z can be determined and according to that $ i.e.: V (aeip) = 0. (10) The angle $ defines a new position of the stagnation point, while the difference of angles is: A$= $- (11) In order to keep the stagnation point at a steady position, a vortex of intensity Ij must be released. Intensity of the vortices that are released can be determined according to the Kelvin's theorem: 2 dr/ dt = 0. Therefore, every change in circulation must be compensated by the vortex inside the cylinder/airfoil of the opposite sign, whose intensity is a function of the variation of circulation around the cylinder: r = ~(drj dt) At (12) and it is equal to the difference of the intensities before and after the free moving vortex is introduced: r =AT1 = a■V^\_sin(a~P[)-sin(a-p). (13) Now the total complex potential has the value: w( z) = Wo( z) + Wr0( z) + wr( z), (14) where: w0 (z) = V„ze-ia + V„ — e+ia + i2aV„sin(a - ß)\m ir, in, I a wr0(z) = —-Mn(z-zo)-—ln 2k o> ~ "'I " zo 2k I z Wr (z) = ir~ ln(z - zi) - '-f1 ln I a— zi 1 2k 2k I z or: w(z) = V ze-ia +Vm— e+ia + i2a -Vm sin (a-ß)ln z + z ir0 w x ir0, f a2 _ +_L ]n( z - z0) r-^ln--zc 2k 2k I z iri w s ir, f a2 _ +_Lin(z - zi)^ln--z 2k 2k ^ z Complex velocity is: V„ e-ia - V„ e+'a + - 2a-V„ z z ir0 1 + irp(a7z 2) + (15) 2k (z-z0) ir1 1 irl(a2/z2) 2k\--z0 2k (z - z1) 2k i — _ z After a period of time At induced velocity at the trailing edge is a consequence of the disposition of both moving and released vortex. Setting again the condition that a point at the circle is the stagnation point and complex velocity at that point equals zero, a new position of the stagnation point by new angle $ is determined. (18) This new angle p defines a new difference AP = P - P. Thus, the new position of the stagnation point is defined by distance: z = a ■ eip. (17) This again causes a generation of a new vortex r2 which is equal to the difference of vortex intensity before and after displacing the free moving vortex from time t1 to time t2 [1]: r2 = Ar = 4k aVa sin (a- P) - -4k aVixsin (a - p). Now the total complex potential has the value: w(z) = Wo (z) + Wro (z) + wrj (z) + wFi (z) , (19) or: 2 w( z) = V, ze-ia + Vm —e+ia + i2aVm sin (a-p)ln z + z ir0 w x ir0, f a2 _ +_K in( z - zo)-Kln I--zo 2k 2k ^ z ir., , , ir. f a2 _ +_L ln( z - z1) - —l- ln|--z, 2k 2k ^ z ir2 w x ir2, f a2 _ +z - z2) - —Knl--z2 2k 2k ^ z Complex velocity is: V = V„e-ia - V„ ar e+ia + - 2a ■ V„ sin (a -p) + z z ir0 1 + r y z2) i---- 2k (z-z0) 2k|IT -z0) ir, 1 r y z2) ß)+ 2k (z - z1) 2k\- z1 ) (16) ir, 1 +—--- . ir2(a2/z2) (20) 2k (z-z2) 2kI a!-z General equations of the following can now be given: • the stagnation point position: z = a ■ (21) • intensity of the released vortex: rm = 4ka ■ V„ [sin (a-$m )- sin (a - $)], (22) • total circulation: rE=r+r0 +fjrm (23) 2 • circulation inside the cylinder/airfoil: - = r + — (24) • complex potential: w(z) = Wo (z) + w- (z) + w- (z) , or: i(r"rm ) (25) w(z ) = V» ze-ia + V» — e+'a + z 2u ln z + +^in(z-zo)-Ufa--zo | + 2n 2u V z -X ^ ln( z - zm) - ^ lnf a- - zm 2u 2u V z complex velocity: „2 ^ = V = V e-ia-V —e dt m » » z2 __ 2u -zo) ir0 2u| z--- z0 a i (-r-rm ) 2u z + (26) —m 1 2U (z -zm} 2uf z-z 2m a Vortices travel down the flowfield by its velocity so that their position is determined by solving the system of Eqs. [1] and [2]: dxi d dt ck dyL = _C dt Cy (27) The position of the vortex at a new moment can be determined by: dZm = Vmdt , (28) and its elementary displacement: Azm = zm - zm = At-V. m m m (29) Fig. 1. Thin vortex surface Now a new distance of each vortex as well as the trajectory of the vortex or any fluid particle is given by [1]: zm = zm+At-V z" = zm +A• r(& Fig. 7. Knot dorming avoidance Fig. 8 Change of distribution and directions of closed vortex nooses 1 £<*> 2 new vortex new boundary layer front new vortex /7777*^ boundary layer Fig. 9. Vortices forming at solid boundary surfaces • a circulation within the boundary layer which passes through a given point at the surface in a unit of time is streamwise and given by equation [7] : dr rV -1 m ( )d \u(rc — = -Jo \n x co (jl)\ut (ll)dl = ^~^~ • 5 HELICOPTER BLADE VORTEX WAKE MODELING Helicopter blade vortex wake modeling includes each previous consideration. That way the changed discrete vortex wake model is consistent with: • a panel model, • a viscous effects model, • a boundary layer simulation model, • a streamline separation model, • a vortex wake model. Therefore, this modeling concept is reasonably adequate for engineering practice. Vortex wake composed of a net of vortex line segments with in general, different vortex intensities is approximated in the following manner: • a change of circulation around the carrying surface is balanced by releasing vortex fibers from carrying surface into a vortex wake; • a new line of vortex line segments occurs in radial direction in each time; • a change of circulation in radial direction around the carrying surface induces vortex fibers in chord direction; • these vortex fibers can also be approximated by vortex line segments; • knot points of vortex net are taken for characteristic points of vortex wake; • vortex wake deformation is accomplished by displacing characteristic points of vortex wake during time; • characteristic point velocity is equal to the sum of undisturbed flow field velocity and the velocity induced by other vortex elements of flow field. Therefore, vortex wake shape corresponding to reality is accomplished. Viscosity is also taken into consideration by the process of additional modeling of vortex wake rolling up, taking place on the edges of vortex wake in real conditions. Vortex line segments change its positions during the deformation process. In this case there is a possibility for the knot points to be found very close to the other vortex line segments, especially in areas of highly deformed vortex wake such as the roll up area and reverse flow area. If a vortex line segment is modeled as potential vortex fiber, it will induce extremely high velocities in very close point. Therefore, it is necessary to use vortex line segments with vortex shell for vortex wake modeling. This vortex shell is used for eliminating singularities and has no physical significiance. Then, velocities in points of flow field close to the vortex line segment have finite values. The edges of the vortex wake are modeled by vortex line segments with smaller vortex shell radius, while internal vortex elements are modeled by greater radius elements, not as a presentation of real physical effects, but as a good approximation of a vortex surface with continually distributed velocities. A value for vortex shell radius of elements at the ends of vortex wake is taken as rc = 0.00275 R, where R is rotor radius. Vortex shells avoid big differences between velocities at the rotor blade in areas close to the vortex wake. Vortex wake effects at a certain distance from the blade can be neglected. It is understandable that the size of the deformation area depends on the regime of helicopter flight. Therefore, as a competent parameter, a coefficient of rotor work regime can be adopted. The area of calculation of vortex wake deformation can be defined through the number of revolutions made by rotor as /u= 0.4 m"1, where m is the number of revolutions for which it is necessary to calculate vortex wake deformation. Afterwards, the vortex wake is 'frozen' in the achieved state and continues to flow downstream, which means that their knot points have the same velocity as undisturbed flow. Therefore, induced velocities in these points are not necessary to compute. The "frozen" part of the vortex wake, however, still effects the vortex wake deformation area that it borders on until the distance does not exceed certain value. Furthermore, the "frozen" part can also be neglected by simple elimination from the model. In this way, the achieved discrete model is consistent with the panel model and the releasing vortices model. 6 VORTICITY IN INCOMPRESSIBLE VISCOUS FLOW In many cases, viscosity of the fluid can be neglected. If it has to be taken into account, some approximations and simplifications can be made, without s significant affect on the final results. Incompressible fluid flow is governed by Navier-Stokes equation: Au = 0, (38) du dt + (u -V) u = - Vp + vV2u. (39) The vorticity of the flow is defined as: o = Vx u . (40) By taking the curl of Eq. (39) the following is obtained: ^ + (u- V)o = o-Vu + vV2o. (41) For two-dimensional flows o- Vu = 0 and Eq. (40) reduces to: do , _2 + (u - V)o = vV o. (42) If o is known than u can be computed using Biot-Savart law. Thus: (x)= x i? 2n i(x^dx. (43) (45) The essence of the inviscid-vortex method is to replace by: o(x) = YriS(* " X), (44) where S are functions approximating the Dirac S -function and r,- is circulation of the i-the vortex. It has to satisfy the kinematics condition: dr _ (u'x-up dt 2 ' where: • drdt - is the change of circulation in the shear layer moving over the stagnation point, • ux', ux" - are velocities on the upper and lower side of the shear layer. In order to satisfy the equation of motion of inviscid flow: D0 = (o-V)-V + vV2o _ 0. (46) The velocity of every vortex must be determined by the value of flow velocity at its instant position: dx, dt L = u(x, t ). (47) Let us assume that a rather small number of vortices per unit length represent the free shear layer. The known assumption of the no-slip condition at the wall (airfoil) leads to the model of wake simulation. If the wake is properly simulated, it will leave airfoil at the trailing edge. In order to obtain proper simulation of viscous effects, it is necessary to assume that the size of the vortex is proportional to the thickens of the wake ct « /8. If the wake length is l then the number of vortices can be defined as: 8-1 l n «—- «-2 ct2 8-f2 f2 where Re is Reynolds number. 6.1 Mapping of a Circular Cylinder to Airfoil Mapping of the circular cylinder to airfoil is done by Joukowsky transformation: where: Ç = z +—, z z = se's + pe' (e-a) (49) (50) is the parametric equation of the function of mapping of 8 and mapping derivative: dÇ i a2 dz z2 (51) 7 MODEL AND CALCULATION ANALYSIS According to the mentioned analysis a numerical model was established. This model is based on the following: • vortices move along the flowfield by the flowfield velocity, • trajectory of every vortex is defined, • for every point in the flowfield it is necessary to determine the value of complex potential and complex velocity, • at a certain moment in a defined initial point a free moving vortex is simulated, • after every time interval At moving vortex changes the stagnation point position which must be compensated by an introduction of a new vortex which brings stagnation point back to its proper place, • new positions of all the vortices are calculated by multiplying local velocities and adding these values to previous ones, • viscous effects should be simulated by generating vortices on the airfoil to satisfy no-slip condition; these introduced vortices must be moved by velocity which is defined by an inviscid part of the equation of motion, • boundary condition of impermeability of the airfoil must be fulfilled, • vortex diffusion is simulated by a variation of the vortex size and arbitrary step, • computer program is made so that flow parameters can be calculated for different angles of attack. The program was first run for characteristic cases of flow around the cylinder and than for flow around airfoil. Due to computer limitations, a rather small number of streamlines is shown. A model has been calculated for different angles of attack, different flowfield velocities, different intensities of the moving vortex and its starting position. The result of a complete calculation for 2D/3D simulation (include comparison with photos from water tunel) is shown in Figs. 10 to 21. Fig. 10. Simulation from water tunel at smaller angle of attack Fig. 11. Moving vortex acting airfoil at leading edge (smaller angle of attack) Fig. 12. Simulation from water tunel at greater angle of attack 8 CONCLUSIONS lift Analysis of results achieved by unsteady modeling and viscous effects simulation 2 method shows that they can be used with sufficient accuracy in rotor analysis and construction. edge (greater angle of attack) Fig. 14. Moving vortex passing airfoil near below (zero angle of attack) Fig. 15. Moving vortex passing airfoil far above (smaller angle of attack) The aim of this particular simulation is to use advantages of vortex methods. For example, vortex methods use the description of flow field of the smallest range; aerodynamic forces can be obtained with a small number of vortices. On the other hand, singular vortex distribution can be accurately determined by using data obtained in small time range. Vortex methods also permit a boundary layer simulation at a large Reynolds's number by a local concentration of computational points. Fig. 16. Moving vortex passing airfoil near above (negative angle of attack) Fig. 17. Moving vortex acting airfoil at leading edge (negative angle of attack) Fig. 18. Moving vortex passing airfoil far above (greater angle of attack) The achieved results imply the direction for further development of this program: • the model should be expanded for transonic flow in the aim of blade tip analysis, Fig. 19. Moving vortex passing airfoil far above (negative angle of attack) mm Fig. 20. 3D free wake geometry (side view) Fig. 21. 3D free wake geometry (isometric view) • a computational dispersion and gradual disappearing of vortex wake are possible by using viscous vortex shell, as a more elegant solution than violent elimination of vortex wake used in this program, • an inclusion of curved vortex elements in the aim of achieving better results from vortex wake self-induction aspect, • from the elastic point of view the blade should be modulate as deformable by using a finite element method. 9 REFERENCES [1] Mitrovic, Č. (2002). Modeling of unsteady helicopter rotor lift. University of Belgrade press, Belgrade. [2] Bengin, A., Mitrovic, Č., Cvetkovic, D., Bekric, D., Pešic, S. (2008). Improved solition approach for aerodynamics loads of helicopter rotor blade in forwerd flight. Strojniski vestnik - Journal of Mechanical Engineering, vol. 54, no. 3, p.170-178. [3] Colin, O. (1973). Unsteady thin-airfoil theory for subsonic flow. AIAA Journal, vol. 11, no. 2, p. 205-209. [4] Sears, W.R. (1976). Unsteady motion of airfoils with boundary layer separation. AIAA Journal, vol. 14, no. 2, p. 216-220. [5] Katz, J., Plotkin, K. (1991). Low-speed aerodynamics - from wing theory to panel methods, McGraw-Hill, New York. [6] Geissler, W. (1978). Nonlinear unsteady potential flow calculations for three-dimensional oscillating wings. AIAA Journal, vol. 16, no. 11, p. 1168-1174. [7] Fishelov, D. (1990). A New Vortex Schemefor Viscous Flows, Journal of Computational Physics, vol. 27, p. 211-222.