Aeronautical Engineering

We are IntechOpen, the world’s leading publisher of Open Access books Built by scientists, for scientists 4,800 122,000 135M Open access books available International authors and editors Downloads Our authors are among the 154 TOP 1% 12.2% Countries delivered to most cited scientists Contributors from top 500 universities Selection of our books indexed in the Book Citation Index in Web of Science™ Core Collection (BKCI) Interested in publishing with us? Contact [email protected] Numbers displayed above are based on latest data collected. For more information visit www.intechopen.com 18 Aeronautical Engineering Časlav Mitrović, Aleksandar Bengin, Nebojša Petrović and Jovan Janković University of Belgrade, Faculty of Mechanical Engineering, Serbia 1. Introduction Aerospace engineering is the primary branch of engineering concerned with the design, construction and science of flight vehicle. Consequently, they are usually the products of various technological and engineering disciplines including aerodynamics, propulsion, avionics, materials science, structural analysis and manufacturing. These technologies are collectively known as aerospace engineering. It is divided into two major and overlapping branches: aeronautical engineering and astronautical engineering. It is typically a large combination of many disciplines that makes up aeronautical engineering. The development and manufacturing of a modern flight vehicle is an extremely complex process and demands careful balance and compromise between abilities, design, available technology and costs. Aeronautical engineers design, test, and supervise the manufacture of aircraft. They also develop new technologies for use in aviation. Aeronautical Engineering is a chapter that encompasses challenging areas such as aircraft design, light-weight structures, stability and control of aeronautical vehicles, propulsion systems, and low and high speed aerodynamics. The field also covers their aerodynamic characteristics and behaviors, airfoil, control surfaces, lift, drag, and other properties. The chapter will include all our research and published papers. 2. Aerodynamics This chapter presents a selection of published scientific papers. The work was the subject of research in the following fields of aerodynamics: Subsonic, Transonic and Supersonic, High Angle of Attack, High Lift, Computational Fluid Dynamics, Wind Tunnel and Flight Testing and Helicopter Rotor Aerodynamics. 2.1 Unsteady motion of two dimensional airfoil in incompressible inviscid flow In this subchapter we shall be looking at many ways in which to solve the problem of unsteady incompressible flow over an aerofoil. The flow being incompressible is a great simplifier to the problem, this allows to take many of the results of steady flow as read. It is still however, not a trivial problem. www.intechopen.com 402 Mechanical Engineering 2.1.1 Introduction Numerical method for solution of unsteady flow over airfoil oscillating in incompressible, inviscid flow field is based on well known panel method in form developed for steady flow by Hess and Smith. Due to unsteadiness in addition to boundary conditions, additional conditions are necessary. Kelvin theorem and Kutta condition is to be applied in unsteady form. This results with a nonlinear problem. 2.1.2 Formulation of problem and boundary conditions Inviscid incompressible fluid flow is described by one partial differential equation which satisfy law of fluid continuity (Laplace equation):  2  2   0  x2  y 2 Equation is the same for steady and unsteady flow. Hence, methods for steady flow can be applied in a solution procedure for unsteady flow. When is the panel method applied, the problem is reduced on determination of the velocity field which must satisfy boundary condition of zero normal flow at the surface of airfoil. When the airfoil is moving relative to a steady undisturbed flow, this boundary condition must be expressed in form:    VF  t   VS  t    n  t   0      where VF is the fluid velocity, VS is the surface velocity and n is normal at the surface. When the airfoil is moving relative to a steady incompressible inviscid flow, circulation (t) around airfoil varying with time t. Total circulation around airfoil and vortex wake must be zero (Kelvin theorem), which is possible to achieve only if the circulation with intensity (/t)t is separated from trailing edge of airfoil in interval dt. This separate vorticity is carried out in downstream direction. There is direct relationship between this and Kutta condition for trailing edge. Kutta condition can be expressed as equality of airfoil upper and lower side pressures on the trailing edge, and on the basis of unsteady Bernoulli equation can be expressed as:   VU  VL 1    VU2  VL2   VU  VL 2 2 t     where U and L denote upper and lower side of trailing edge, respectively ( x l  1 ). Therefore, above mentioned consistence of Kutta condition and vortex shedding model can be clearly seen. Namely, during time above mentioned circulation (t) around airfoil is balanced by vortex shedding with intensity VU  VL , with mean velocity VU  VL 2 .     Therefore, unsteadiness of the problem is taken into consideration through unsteady form of kinematics boundary condition, Kelvin theorem, unsteady Bernoulli equation, and unsteady form of Kutta condition for trailing edge. www.intechopen.com 403 Aeronautical Engineering 2.1.3 Discretization and numerical solution procedure Solution of flow over airfoil (moving arbitrary), depending on time and starting from time t=0 is calculated in successive time intervals tk (to= 0, k=1,2,3...). Fig. 1 shows model for the time tk. Airfoil contour in time tk is replaced by N linear elements. (i)k and k are uniform source and circulation distributions on i-th element (i=1,2,...,N), where (i)k varies from one element to another, k is the same for each element on airfoil and k denotes time tk. Total circulation k is given as k (airfoil perimeter). An elementary vortex wake with length of k and pitch angle of k in regard to the x-axis (to the free stream direction) is attached to the trailing edge. Length k and angle k are arbitrary in first iteration. Their values will be determined as the part of the solution. Circulation in trailing edge vortex wake element is (w)k, where:  k   w  k   k   k 1 (1) Hence, circulation on the element is equal to the difference between circulations around airfoil in times tk-1 and tk, assuming that k-1 has been already determined. Vortex wake consist of concentrated vortices formed by vorticity shed at earlier times, which is assumed to be transformed into discrete vortices. Concentrated vortices is moving with resulting velocity calculated in the center of each vortex at each successive time interval. Therefore, strength and positions of discrete vortices are regarded as known at time tk, and there are N+3 unknowns (i)k (i=1,2,...,N), k, k and k at the time tk,. Fig. 1. Solution at time tk They are determined by satisfying following conditions:  N conditions of zero normal velocity component in external middle point at each segment of airfoil (Vnj )k  0 www.intechopen.com (2) 404 Mechanical Engineering where (Vnj )k is total normal velocity component at the external middle point of j-th. element in time tk.  condition of equal pressures in middle points of two elements on an airfoil on both sides of trailing edge (unsteady Kutta condition): (Vt1 )2k  (VtN )2k  2   k   k 1  /  tk  t k 1  (3) where (Vt1 )k is total tangent velocity component at middle point of first element and (VtN )k is total tangent velocity component at middle point of N-th element in time tk.  trailing edge vortex wake element length and direction (k and k ) are determined from condition that the element is tangent to the local resultant velocity and that its velocity is proportional to the local resulting velocity. If (Uw)k and (Ww)k are components of total velocity at middle point of trailing edge vortex wake element, excluding influence of the element on itself, then: 1 tan  k   Ww  k U w  k , 2 2  k  U w  k   Ww  k  2 t k  tk 1    (4) By taking into the consideration that the problem considers incompressible flow, formulas for velocities induced from source and vorticity distributions are the same as these for the steady case. Basic set of equations is nonlinear therefore following iterative procedure for its solution is accepted. Since k and k are assumed it leaves N linear equations from (2) and square equation from (3). N linear equations are solved to give (i)k (i=1,2,...,N) relatively to the k. Then k is determined from square equation (3). Once known (i)k (i=1,2,...,N) and k, (Uw)k and (Ww)k can be determined and replaced into (4) to obtain new values of k and k. Procedure is repeated until k and k are of requested accuracy. When intensities are determined, source and vortex distribution is known. Pressure coefficient follows from unsteady form of Bernoulli equation: cp  1  V2 2   2 2 U U  t (5) where V is total velocity on external side of airfoil, and  is velocity potential. Forces and moment coefficients are obtained by direct integration of pressure distribution. Since solution in time tk is determined, model is applied for the time tk+1 in the same way as this for time tk, with vortex vale computed in time tk. Distributed vorticity in trailing edge vortex wake element in time tk is now considered as concentrated in vortex with strength (w)kk at the position in time tk+1 determined by coordinates: 1 1 x   xTE  k   k cos k  U w  k  tk 1  t k  , z   zTE  k   k sin  k   Ww  k  t k 1  t k  2 2 Resultant velocity in the center of each next concentrated vortex in vortex wake is calculated by solution at time tk, and their position in time tk+1 directly follows. www.intechopen.com Aeronautical Engineering 405 2.1.4 Discussion of results Computer program was developed from formulated numerical method. Flow over NACA 0012 airfoil was modeled by using 20 linear elements. Airfoil oscillating at high frequency trough certain angle of attack is considered flow characteristics are computed in time interval from 0 to 30 sec. Fig. 2 shows static and dynamic lift curve resulting from the application of this numerical method. Dynamic lift curve is obtained by airfoil oscillating through angle of attack of 6, with amplitude 5 and reduced frequency k=0.5. Fig. 2. Static and dynamic lift curve Fig. 3. Influence of reduced frequency Dynamics lift curves for two reduced frequencies k=0.5 and k=0.15 based on obtained results are also shown (Fig. 3). In this case, airfoil oscillating through angle of attach of 2.5, with amplitude of 4. Dynamic lift curve shape depending on reduced frequency can be seen from the figure. Fig. 4. Computed vortex traces as a function of the reduced frequency k,   2.5  4 sin t  www.intechopen.com 406 Mechanical Engineering Vortex wake shape of airfoil oscillating through same angle of attack, with same amplitude, but with different reduced frequency obtained by this method is shown at Fig. 4. It is similar with one obtained by methods of visualization in tunnels. Numerical method presented in this paper leads to inviscid flow field calculation over airfoil moving arbitrary in time if it is supposed that flow stays attached and that it is separated art trailing edge. Method is completely general although only results of oscillating at high frequencies are presented. 2.2 An optimal main helicopter rotor projection model obtained by viscous effects and unsteady lift simulation This subchapter presents a model of airfoil which allows obtaining a method for optimal main helicopter rotor projection by viscous effect and unsteady lift simulation through algorithm and set of program entireties, applicable to ideological and main project of helicopter rotor. Numerical analysis considerations in this paper can be applied, on basis of a real rotors theoretical consideration, with sufficient accuracy in analysis and constructive realizations of helicopter rotor in real conditions. The method for unsteady viscous flow simulations by inviscid techniques is developed. This allows the to define such optimal conception model of aerodynamic rotor projection which is corresponding to rotor behavior in real condition and with sufficiently quality from the aspect of engineer use. 2.2.1 Introduction The basic aim is to determine the lift of helicopter rotor blades by using singularity method, that is 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 singularity method is based on the need for avoidance of very expensive experiments in the first stage of rotor planning by using contemporary aerodynamic analysis applied to available computer technique. It is necessary to modulate vortex wake, unsteady 2D flow field characteristics and blade dynamic characteristics for determination of aerodynamic forces acting upon helicopter rotor blade. In the first part airfoil is approximated by vortex and source panels, and boundary layer by layer of vortices changing their position during time. Vortex trail and trail of separation are modulated by free vortices. The separation is modulated by free vortices. The dependence of coefficient of lift on angle of attack is determined for 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 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 what it is required generation of new vortices that would all together satisfy the airfoil contour boundary conditions. In addition to a complex problem of analytical modeling these phenomena there is a problem of modeling interactions between effects of these phenomena. This approach to www.intechopen.com 407 Aeronautical Engineering helicopter rotor planning allows experimental investigations in tunnels and in flight to be final examinations. In that way, a very useful interaction between numerical calculation and experimental results is achieved. The finale result is obtaining of unsteady lift dependence of angle of attack. Model established in such a way is characteristic for the helicopter rotor blade airflow and it is based on the influence of the previous lifting surface’s wake influence on the next coming blade. 2.2.2 Foundations of the irrotational 2-D flow The planar potential flow of incompressible fluid can be treated in Cartesian coordinates x and y. Two dimensional potential incompressible flow is completely defined by the speed potential and stream function and is presented by Cauchy-Reimann equations:     0 x y i     0 x y where u and v are velocities in x and y directions respectively. The Laplace partial differential equations can also be introduced  2  0 and  2  0 . The fulfillment of Cauchy-Reimann conditions enables combining of the velocity potential and stream function w  z    ( x , y )  i ( x , y ) . 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  u  iv  V dz where u is its real part (velocity component in x-direction) and v is the imaginary part (velocity component in y-direction). Circulation and flow are equal to zero for any closed curve in complex plane. Complex potential has no singularities except at the stagnation point. 2.2.3 Simulation of the moving vortex If we assume that the moving vortex of intensity 0 is at the distance z0, than the perturbance potential of such a flow is: w( z)  V ze i  V  a2  i i a 2  i e  i 2a  V  sin     ln z  0 ln( z  z0 )  0 ln   z0  2 2  z z  and it’s complex velocity: V  V e i  V www.intechopen.com i0 ( a2 z2 ) i0 a2  i i 1   e a V sin       2    2 ( z  z0 ) z  a2  z2 2   z0   z  408 Mechanical Engineering In this case displacement of the stagnation point appears, which is now at the distance z  a  ei   . According to the request that complex velocity at the stagnation point must be zero: dw V  0 dz we can determine position z and according to that ' i.e. V ( aei  )  0 . The angle ' defines new position of stagnation point, while the difference of angles is  =  - ' . In order to keep the stagnation point at the steady position, a vortex of intensity 1 must be released.Intensity of the vortices that are released can be determined according to the Kelvin’s theorem d dt  0 . So every change in circulation must be compensated by vortex inside of cylinder/airfoil of the opposite sign, whose intensity is a function of the variation of circulation around the cylinder i    d dt  t and it is equal do the difference of the intensities before and after the free moving vortex is introduced 1  1  4 a  V sin   1   sin      . Now the total complex potential has the value: w( z)  w0 ( z)  w0 ( z)  w1 ( z) where: w0 ( z)  V ze i  V w0 ( z)  a2  i e  i 2aV sin     ln z z  a2  i 0 i ln( z  z0 )  0 ln   z0  , 2 2  z  w1 ( z )   a2  i1 i ln( z  z1 )  1 ln   z1  2 2  z  Complex velocity is: a2  i i e  2a  V  sin     z z2 2 i 0 ( a z 2 ) i1 1 i1( a2 z 2 ) i 1  0    2 ( z  z0 )  a2  2 ( z  z1 )  a2  2   z0  2   z1   z   z  V  V e i  V (6) After a period of time t 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 stagnation point and complex velocity at that point equal to zero, we determine a new position of the stagnation point by new angle ' . This new angle ' defines a new difference  =  - '. So the new position of the stagnation point is defined by distance z  a  ei   . This again causes www.intechopen.com 409 Aeronautical Engineering generation of a new vortex 2 which is equal do the difference of vortex intensity before and after displacing of the free moving vortex from time t1 to time t2.  2    4 aV sin       4 aV sin     Now the total complex potential has the value: w( z)  V ze  i  V  a2  i i a 2  i e  i 2 aV sin     ln z  0 ln( z  z0 )  0 ln   z0  2 2  z z   a2  i  a2  i i i  1 ln( z  z1 )  1 ln   z1   2 ln( z  z2 )  2 ln   z2  2 2  z 2  z  2  (7) Complex velocity is: V  V e i  V i0 ( a2 z2 ) i0 1 a2  i i   2 e a V sin         2 ( z  z0 ) z  a2  z2 2   z0   z  (8) i1( a z ) i 2 i 2 ( a z ) 1 i1 1    2 2 ( z  z1 ) a  2 ( z  z2 )  a2  2   z1  2   z2   z   z  2  2 2 2 Vortices travel down the flowfield by its velocity so that their position is determined by solving the system of equations: dxi   , dt  x i dyi   dt  y . i Position of the vortex at a new moment can be determined by dzm  Vm dt and its elementary displacement zm  zmn  zms  t  V . Now a new distance of each vortex as well as the trajectory of the vortex or any fluid particle is given by zmn  zms  t  V , or:  i    m  i0 i0 1 a2 zmn  zms  t  V e  i  V 2 e  i     2 z 2 ( z  z0 )  z z2   2  z  2 z0  a       n i m 1  i m    2 ( z  z )   z2    m 1  2  z  2 zm     a     www.intechopen.com (9) 410 Mechanical Engineering 2.2.4 Calculation According to the whole 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 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 t moving vortex changes the stagnation point position which must be compensated by introduction of new vortex which brings stagnation point back to its proper place new positions of all vortices are calculated by multiplying local velocities and adding these values to previous boundary condition of impermeability of the airfoil must be fulfilled vortex diffusion is simulated by variation of the vortex size and arbitrary step computer program is made so that flow parameters can be calculated for different angles of attack. 2.2.5 Analysis of the calculation Program was run first for characteristic cases of flow around the cylinder and than flow around airfoil. Due to computer limitations, a rather small number of streamlines is shown. Model is calculated for different angles of attack, different flowfield velocities, different intensities of the moving vortex and its starting position. The result of complete calculation for 2/3D simulation (include compare with photos from water tunel) is shown at figs. 5-12. Fig. 5. Fig. 6. Fig. 7. Fig. 8. www.intechopen.com 411 Aeronautical Engineering Fig. 9. Fig. 10. Fig. 11. Fig. 12. 2.2.6 Conclusion remarks Analysis of results achieved by unsteady lift modeling and viscous effects simulation method shows that they can be used with sufficient accuracy in rotor analysis and construction. 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 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 boundary layer simulation at large Reynolds’s number by local concentration of computational points. Achieved results imply the direction for the further development of this program    the model should be expanded for transonic flow in the aim of blade tip analysis a computational dispersion and gradual disappearing of vortex wake are possible by using viscous vortex shell, as more elegant solution than violent elimination of vortex wake used in this program an inclusion of curved vortex elements in aim of achieving better results from vortex wake self-induction aspect On the basis of analysis of presented model and program package for viscous effects and unsteady lift simulation it can be concluded that this subchapter presents an original scientific contribution, applicable in aerodynamic analyses of practical problems in helicopter rotor projecting. www.intechopen.com 412 Mechanical Engineering 2.3 Improved solution approach for aerodynamics loads of helicopter rotor blade in forward flight This subchapter presents the numerical model developed for rotor blade aerodynamics loads calculation. The model is unsteady and fully three-dimensional. Helicopter blade is assumed to be rigid, and its motion during rotation is modeled in the manner that rotor presents a model of rotor of helicopter Aerospatiale SA 341 “Gazelle”. 2.3.1 Introduction Helicopter rotor aerodynamic flow field is very complex, and it is characterized by remarkably unsteady behavior. The most significant unsteadiness appears during the forward flight. In that case, the progressive motion of helicopter coupled with rotary motion of rotor blades causes drastic variations of local velocity vectors over the blades, where the advancing or retreating blade position is of great significance. In first case, the local tip transonic flow generates, while in second, speed reversal appears. In addition, in forward flight, blades encounter wakes generated by forerunning blades and so encounter nonuniform inflow. The wake passing by the blade induces high velocities close to it causes changes in lifting force. Besides that, in horizontal flight blades constantly change pitch, i.e. angle of attack at different azimuths. Such angle of attack variations are very rapid, so that dynamic stall occurs, especially in case of retreating blades. 2.3.2 Dynamics In this paper the rotor of helicopter SA 341 “Gazelle” is modeled, at which the blades are attached to the hub by flap, pitch and pseudo lead-lag hinges. Blade motion in lead-lag plane is limited by dynamic damper, which permits very small maximum blade deflection. Fig. 13. Coordinate systems Due to such small angular freedom of motion, we can assume that there is no led-lag motion at all. Pitch hinge is placed between flap hinge and pseudo lead-lag hinge. According to that, the following frames have been selected for use: www.intechopen.com 413 Aeronautical Engineering     fixed frame F, with x-axis in the direction of flight, and z-axis oriented upwards; frame H, connected with the rotor and rotates together with it; it is obtained by rotating the F frame for a certain azimuth angle ; keeping the common z-axis; frame P, connected to the flapping hinge, so that the y-axis is oriented along the blade; its origin is displaced from the rotating axis for the value e , while it is rotated for the angle  with respect to the frame H; frame B, connected to the blade, displaced for the value e from the origin of P, and tilted for the value  (pitch angle) with respect to the P. In derivation of the equations of motion, the following was assumed:   the rotor does not vibrate, and its rotation velocity  is constant, the blade is considered absolutely rigid. With assumptions above mentioned, equation of blade flapping motion is:   B   2 B cos   mb e x g R sin   M A where B is moment of inertia about Px , mb mass of the blade, x g position of blade center of gravity in P frame, and M A is aerodynamic moment. 2.3.3 Aerodynamics The flow field is assumed to be potential (inviscid and irrotational) and incompressible. In that case, velocity potential satisfies the Laplace equation   0 . The equation is the same, both for steady and unsteady flows. Owing to that, methods for steady cases can be applied for the solution of unsteady flow problems, as well. Unsteadiness is introduced by:     unsteady boundary condition: V  VT  n  0   D 0 Dt  Kelvin theorem:  unsteady form of the Bernoulli equation: p  p   V2  V 2  2 0.5   t   where:  - is velocity potential, V - is absolute fluid velocity, VT - is lifting surface velocity,  n - is the normal of the lifting surface at a certain point, and  - is the bound circulation. In case of inviscid problems, it is necessary to satisfy Kutta condition at the trailing edge. Based on unsteady Bernoulli equation, the difference between upper and lower surface pressure coefficients, in case of the thin lifting surface, is: C P  C PU  C PL   www.intechopen.com VU2  VL2 V2  M  2   dl    2  V t  LE  (10) 414 Mechanical Engineering where integral should be calculated from leading edge to a certain point M at the surface, and  is the local bound vortex distribution. If we assume that spanwise velocities are small, the difference of velocity squares can be calculated as: VU2  VL2  2V (11) By substituting (11) in (10), the following equation can be obtained: C P   M  2   V    dl    2   t LE V   The Kutta condition can be expressed as the uniqueness of pressure coefficients at the trailing edge: TE  2   2    C p  2  V TE      dl V     V 2   TE t   0  t LE V    where  is the contour circulation which covers the lifting surface. Since it is impossible to be V   , the relation within the parentheses must be equal to zero and the expression for unsteady Kutta condition comes out directly as:   V TE t If the right hand-side part is substituted with (11) written for the trailing edge, we obtain: VU2TE  VL2TE VUTE  VLTE     VUTE  VLTE t 2 2   From this equation, it can be clearly seen that the variation of the lifting surface circulation in time can be compensated by releasing vortices of magnitude VUTE  VLTE at the velocity VUTE  VLTE 2 .     2.3.4 Discretization and numerical solution procedure The method for the solution of this problem is based on the coupling of the dynamic equations of blade motion with the equations of aerodynamics. Discretization in time is done by observing the flow around the blade in a series of positions that it takes at certain times t k , which are spaced by finite time intervals t at different azimuths. Also, discretization of the thin lifting surface is done by using the panel approach. By this method, the lifting surface is divided in a finite number of quadrilateral surfaces – panels. Vorticity distribution is discretized in a finite number of concentrated, closed quadrilateral linear vortices, in such a way that one side of the linear vortex is placed at the first quarter chord of the panel, and represents the bound vortex of the corresponding panel. www.intechopen.com Aeronautical Engineering 415 Fig. 14. The steady panel scheme The opposite side of the vortex is always placed at the trailing edge, while the other two sides are parallel to the flow. The wake is represented by quadrilateral vortex in the airflow behind the lifting surface. One side of it is connected to the trailing edge, while the opposite one is at the infinity. The other two sides (trailing vortices), which actually represent the wake, are placed parallel to the airflow. The vorticity of the quadrilateral vortex is equal to the sum of the vorticities of all bound vortices of the panels that correspond a certain lifting surface chord, but opposite in direction. Then the trailing edge vorticity is equal to zero. Model established in such a manner corresponds to the steady flow case. On the other hand, it can be very easily spread in order to include the unsteady effects. 2.3.4.1 Vortex releasing model The variation of the lifting surface position in time induces variation of circulation around the lifting surface as well. According to the Kelvin theorem, this variation in circulation must also induce the variation around the wake. According to the unsteady Kutta condition, this can be achieved by successive releasing of the vortices in the airflow. Fig. 15. Vortex releasing Suppose that the lifting surface has been at rest until the moment t, when it started with the relative motion with respect to the undisturbed airflow. The vortex releasing, as a way of circulation balancing, is done continually, and in such a way a vortex surface of intensity www.intechopen.com 416 Mechanical Engineering   t  is formed. At the next moment t  t , the flow model will look like in Figure 15. The circulation of the vortex element joined to the trailing edges is equal to the difference in circulations at moments t  t and t. Fig. 16. Unsteady panel scheme We will discretize the vortex “tail” by replacing it with the quadrilateral vortex loop, whose one side is at the trailing edge, and the opposite side is at the finite distance from the trailing edge (shed vortex). By this, we can obtain the final model for unsteady case. 2.3.4.2 Discretization of wake The established vortex-releasing model is appropriate for the wake modeling using the “free wake” approach. During the time, by continuous releasing of the quadrilateral vortex loops, the vortex lattice formed of linear trailed and shed vortices is created. The wake distortion is achieved by altering the positions of node points in time, by application of kinematics    relation ri  t  t   ri  t   Vi  t  t . The velocities of the collocation points are obtained as sums of the undisturbed flow velocity and velocities induced by other vortex elements of the flow field. Induced velocities are calculated using the Biot-Savart law. In order to avoid the problems of velocity singularities, line vortex elements are modeled with core. The core radius varies with the gradient of the bound circulation at the position where vortex is released. Fig. 17.Discretized wake www.intechopen.com 417 Aeronautical Engineering The wake influence at large distances from the blade is negligible, so it is possible to neglect the wake distortion at a sufficient distance from the rotor. The number of revolutions for which it is necessary to calculate the wake distortion can be determined as m  0.4  , where  is the advance ratio. After that, the wake shape is “frozen” in achieved state, and it moves by flow field velocity, keeping it for the rest of the time. The “frozen” part of the wake still influences the adjacent area in which it is still being distorted. After some distance, even the influence of the “frozen” part becomes negligible, and then it is eliminated from the model. 2.3.4.3 Definition of equation set The boundary condition of impermeability of the lifting surface should be satisfied at any moment of time, tk  k  0 ,1, 2 in a finite number of points of lifting surface   Vi  VT i   ni  0 ; i  1, 2, , n Points at which this condition must be satisfied are called the control points. One of them is placed on each panel, at the three-quarter chord panel positions. By this, at every moment of time, the number of lifting surface impermeability conditions is equal to the number of unknown values of circulations of bound vortices. The equations of motion  of the lifting surface are known, as well as the velocities VT i of all characteristic points,  and their normals ni as well. At each flow field velocity can be divided to the  point,   free stream velocity and perturbation velocity: Vi  V  wi . The perturbation velocity is induced by lifting surface and wake vortex elements. It is calculated by Biot-Savart law. At every moment, the wake shape and circulations of its vortex lines are known, and so the wake-induced velocity at every flow field point is known as well. On the other hand, the circulations of the bound vortices are unknowns (their positions are defined by the lifting surface shape). The boundary condition for the i-th control point can be written as: n  ai , j  ji  t   bi i 1 where ai , j are the coefficients depending of the blade geometry, and bi are the coefficients containing the influence of the wake and free-stream flows. This way, by writing equations for all control points, the equation set of the unknown bound circulations is obtained. Along with this equation set, the Kutta condition must be satisfied. In case of numerical solutions, it is customary to satisfy Kutta condition in vicinity of the trailing edge. According to that, we obtain discretized form of unsteady Kutta condition: V  n   t  t     t   0 ln t where  n ln is intensity of the distributed vorticity at the trailing edge panel  n , and ln is the trailing edge panel cord length. By adding the Kutta conditions to the equation set, an over-determined equation set is obtained. It can be reduced to the determined system by the www.intechopen.com 418 Mechanical Engineering method of least squares. After that, it can be solved by some of the usual approaches, by which the unknown values of circulations  i at the time t are obtained. 2.3.4.4 Determination of the aerodynamic force After unknown circulations  i are obtained, velocity at every point of the flow field is known, and we can use them for the determination of aerodynamic forces that act on the blade. The calculation aerodynamic force is necessary for the defining of the blade position at the next moment of time. The total aerodynamic force is calculated as the sum of forces acting on all panels. Aerodynamic force acting on a single panel can be defined by:    Fi  V  i  ef . BC  where BC is bound circulation vector and effective circulation can be defined by using: M i  1  i 1   i 1 dl  i  ef .   i        k   i  V t LE V t  k 1 4 where Mi denote point at quarter-chord position on the i-th panel. After determination of aerodynamic forces, moment of aerodynamic forces M A , necessary for blade flapping equation, can be calculated as sum of moments acting flapping hinge. 2.3.5 Results Presented method for blade air-loads calculation has been implemented in a computer code. Forward flight application is presented for SA 341 “Gazelle” helicopter at moderate forward speed condition (advance ratio  is 0.35). Fig. 18. Calculated free wake geometry;  = 0.35 (3D view) www.intechopen.com Aeronautical Engineering 419 By analyzing the drawing of the blade wakes, it can be concluded that model applied in this paper gives reasonable simulation of actual wake behavior, specially in the domain of wake boundaries, where wake roll-up occurs (although it is slightly underestimated compared with existing experimental data). In addition, larger wake distortion in the domains of the forerunning blades or their wakes is noticeable. Fig. 19. Circulation distribution on the rotor disk; =0.35 (different points of view) The program results (Figure. 19), show the difference in circulation distributions at different azimuths, as well as the disturbances caused when blades are passing the wakes of other blades, and the characteristic reversal flow domains. 2.3.6 Conclusion Obtained results can define suggestions for the future solution improvement. Firstly, by incorporating the transonic flow calculations, the advancing blade tip simulation would be more appropriate. Secondly, viscous interaction should be included as well, which would improve the wake roll-up simulation. The viscous vortex core simulation in time would improve the results concerning the wake vanishing effects far enough from the blade. Finally, introduction of the curvilinear vortex elements would give better results from the aspect of wake self-induction. With these enhancements, this solution should prove to be a useful and efficient tool to the rotorcraft performance evaluation. 3. Flight dynamics and control The flight dynamics mathematical model of an aircraft, which would be strictly determined, would comprise of a system of non-linear, non-stationary, partial differential equations. To simplify these equations we introduce a number of assumptions. Ignored are the elastic characteristics of the aircraft so the aircraft can be thought as a rigid body and, in that way, the dispersal of parameters is eliminated. Also, fuel consumption is disregarded and so the non-stationarity due to temporal change in helicopter mass is eliminated. 3.1 Mathematical models of helicopter flight dynamics The helicopter is specific in regards to other traffic-transportation means, not just by its structure but also by its motion possibilities. The helicopter can move vertically, float in the air, www.intechopen.com 420 Mechanical Engineering turn in place, move forward and lateral, and can perform these movements in combinations. Because of this, helicopter dynamics modeling and testing is a very complex problem. In the present, problems in helicopter flight dynamics are mostly solved in aid of modern computers. Though inevitable in many complex problems, computers do not make it possible to understand the physical nature of the problem. Fortunately, many problems considering helicopters can be analyzed without overly complex calculus and usually it is possible to obtain simple formulas. Though not suitable for calculus, these formulas, when designing the helicopter, enable a satisfactory interpretation of required aerodynamic and dynamic phenomena. The mathematical model described in this paper is related to three-dimensional (space) geometry and kinematics, and rigid body dynamics and fluid dynamics through which it moves. 3.1.1 Motion of supporting rotor blades To understand the flight dynamics of helicopters and determine dynamic moments and forces that act upon the helicopter, it is a necessity to preinvestigate the motion of supporting rotor blades. From a vast number of different types of helicopters, the single rotor helicopter was chosen that has its blades coupled with the main rotor by a hinge about which they can freely move. It should be noted that there are also rotors that have blade fixedly connected to the hub. 3.1.1.1 Equations of blade fluttering Rotor blades are regarded to as rigid body. The horizontal hinge is placed at length eR from rotation axis. The shaft rotates at angle speed =const, and the blade flutters at angle speed of d/dt. The axis that passes through the blade is parallels to the axis of inertia of the blade and passes through the hinge (Figure 20). Fig. 20. Blade fluttering scheme Paremeter R represents the length of the blade,  represents the flutter angle of the blade. Following a complex calculus, obtained are the equations for blade fluttering:   2 1      M Ay / J y   J x   cos   J x    sin   0 2 J y   sin   M z www.intechopen.com (12) 421 Aeronautical Engineering 3.1.1.2 Equations of blade lead-lag It is assumed =0 and that the blade is moving forward in relation to the vertical hinge by the lead-lag angle amount . The vertical hinge is placed at distance eR from the shaft axis. The coordinate system is placed as in the previous case. From this follows the equation for blade lead-lag:   2   2     M z / J z Fig. 21. Blade lead-lag scheme If the azimuth angle is described as =t, then follows: Mz d 2 d   2   2 d J z 2 d 3.1.1.3 Equation of blade climb It is assumed that flutter and lead-lag angles are equal zero. The blade step is the angle between the blade cross section chord and the plane of the hub, designated as k. Figure 22 shows the coordinate system attached to the blade. Equations of blade motion about longitudinal axis are: Fig. 22. Coordinate system at blade cross section k  2  k  Mx / J x  2 J z  k sin  k  M z J y  k cos k  J y k cos k  0 3.1.2 Rotor forces To project forces the following axis may be use: control axis, rotor disc axis which is normal to the rotor plane, that is, to the plane on which reside blade tips, and shaft axis. www.intechopen.com 422 Mechanical Engineering Once the axis is chosen, the remaining axis of the coordinate system will be normal to it and pointed lateral, that is, to the tail of the helicopter. Customarily, the force component along the chosen axis is said to be the tow force, the force component pointed towards the tail is said to be H force, and the force component pointed lateral is said to be Y force. If the force components are designated without subscripts, it is assumed they are determined relative to control axis, whereas subscripts “D” and “S” are used when relating to rotor axis, that is, shaft axis. Since flutter and mount angles are usually small (amounts greater than 10° are considered extreme), the relation between these components can be obtained: T  TD  TS H  H D  TD a1  H S  TS B1 , 3.1.2.1 Longitudinal equilibrium of forces Angle 1 is the longitudinal amplitude of a cyclic change in blade step; angle a1s is the angle between shaft and axis of rotor disc. After extensive calculus the expression for longitudinal amplitude of cyclic change in blade step is obtained: B1  M f  G  fR  H  hR  MS  a1 T  hR  MS Fig. 23. Drawing for determining longitudinal equilibrium of forces For e=0, it can be adopted that Ms=0 and Mf=0, and since T=G, follows: B1   f H  h G , p   Mf f D cos   G h G  hR Above equation has a simple physical interpretation: the amplitude of the longitudinal cyclic control must have such a value in order to position the direction of the resultant rotor force through the center of mass. 3.1.2.2 Lateral equilibrium of forces Angle A1 presents the amplitude of lateral cyclic change in blade step of the supporting helicopter rotor: www.intechopen.com 423 Aeronautical Engineering A1  - G  f1 R  MS  b1  Tt  ht R G  hR  M S Fig. 24. Drawing for determining lateral equilibrium of forces By replacing value A1 into corresponding equations, we obtain the value of angle , which determines the position of the fuselage.  Tt G  f1R  MS  b1  Tt  ht R  G G  hR  MS f1 , which means the rotor hub is h positioned vertically above the center of mass. All values of these determined angles are so called trimmed values. If MS=0 and ht=h, which can often be assumed, follows   3.1.3 Non-linear mathematical model of flight dynamics Mathematical modeling of helicopter motion is a very complex task and, therefore, it is necessary to introduce series of assumptions and approximations. Knowledge of motion of individual helicopter blades is not necessary for investigating dynamic characteristics of the helicopter, except in a special case, but rather for defining forces and moments in a disturbed flight it is sufficient enough to view the rotor as whole. Because of a great number of different helicopters, in this paper a single rotor helicopter was studied, that has its blades connected to the hub by hinges. As mentioned before, the helicopter can perform different motions and it would be very difficult to make a mathematical model that would combine all those motions. It is assumed the helicopter is airborne and in straightforward flight. It is required that the helicopter, at straight forward flight, has following velocity components: Wx, Wy, and Wz, at nominal values, and angle of turn , angle of roll , angle of climb , as long as the intensity of disturbance is in permitted limits. Figure 25 presents a schematic of the helicopter with a floating coordinate system tied to its center of mass, and Figure 26 presents a helicopter block diagram. www.intechopen.com 424 Mechanical Engineering Fig. 25. Schematic of helicopter Fig. 26. Helicopter dynamics block diagram After introducing a series of assumptions, such as:     helicopter mass is a constant value, helicopter is a rigid body, 0xz is a plane of symmetry, angle increments , ,  are too small, and so on; and we come to a non-linear mathematical model by deviations in the form: d  Wx  d  Wz  dt    1 f1 Wx , Wz ,  , u1 , u2   mg cos    m   1 f 2 Wx , Wz ,  , u1 , u2  WzN m   mg sin     dt m d  d    1  ,  , u , u ,    f 3 Wx , Wz , W z 1 2 dt dt Jy     d Wy   1 f d    dt www.intechopen.com   Wy ,  ,  , u3 , u4   WzN m  mg cos  mg sin  dt m d      dt 4    , , d    dt  d    dt   1  f 5 Wy ,  ,  , u3 , u4  J xz    Jx     1 f 6 Wy ,  ,  , u3 , u4  J xz    Jz  425 Aeronautical Engineering Where:     U1=B1 – amplitude of cyclic change in step in the longitudinal direction (regarding to longitudinal motion), U2=0 – change of collective step of the helicopter rotor blade (regarding to longitudinal motion), U3=A1 – amplitude of cyclic change in step in the lateral direction (regarding to lateral motion), and U4=t – change of collective step of the tail rotor (regarding to lateral motion). 3.1.4 Linearized mathematical model of flight dynamics In technical applications it has been shown that, with an acceptable accuracy, linearized mathematical models may be used under the condition that deviations of physical quantities from their nominal values are small. The outcome of adopted presumptions is that the output values, input values, and the vector of state for both longitudinal and lateral motion will be: X   X1  X9  T , X i   Xi1  Xi 6  , u   u 1  u4  T T (13) The vector equation of state for the linearized mathematical model with non-dimensional quantities, deviations, that is, quantities of state is shown in equation 13. Equation 14 present the equation at exit.  a11 a12 a  21 a22 0 0   a41 a42 0 X   0  0 0 0 0  0 0  0 0 a13 a23 0 a43 0 0 0 0 0 a14 a24 1 a44 0 0 0 0 0 1 0  0 Xi   0 0  0 0 0 0 0 a55 0 a75 0 a95 0 0 0 0 a56 0 0 0 0 0 0 0 0 0 1 a77 0 a97 0   b11 b12 b  0   21 b22 0 0  0   0  b41 b42 0 a59  X   0   0  0 0   a79 0 0  0 1  0   a99  0 0 0 0 0 0 a58 0 0 0 0 0 0  0 0  0 0   0 0  b53 b54  u  0 0  b73 b74   0 0   b93 b94  (14) 0 0  0  X 0 0 0 0 0 1 0 0 0  0 0 0 0 0 1 0 0  0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 Besides the way this is presented, in a form of common matrix, also should be noted that longitudinal and lateral motions are separated, because this was the condition for deriving this mathematical model. www.intechopen.com 426 Mechanical Engineering In equations above, equations for longitudinal motion are presented within the first for rows of the matrices, while the remaining five rows present the equation of state and equation of lateral motion. Designation used in above equation is: C*  a 13 1 / ix i z  mc cos 2 1  ixz a12  x w a14  xq a 21  zu a 23  mc sin  a 22  zw a 56  mc cos a 55  y v ˆ a W 59 a11  xu a 58  mc sin  a 41  mu  mw zu xN a 42  mw  mw zw ˆ z  mq  mw W xN q a 43  mw mc sin  a44 b11  xB1 ˆ z a24  W xN q b12  zB1 b 41  mw zB1  mB1   b 21  x0 b22  z0 b 42  mw z0  m0  b74   l  n ixz / ix  C * b94   n  l ixz / iz  C * a75   nv ixz / ix  lv  C * a77   lp  np ixz / ix  C * a95   nv  lv ixz / iz  C * a97   np  lpixz / iz  C * a99   nr  lr ixz / iz  C * b93   nA  lA ixz / iz  C * b 53  y A1 t  b54  yt b73  lA1  nA1 ixz / ix C * t t 1 a79   lr  nr ixz / ix  C t 1 * 3.1.5 Program results The program was tested on the example of a single rotor helicopter which main rotor blades are tied to the hub over hinges. The helicopter is described by the following input data: helicopter weight G=45042N, rotor abundance degree s=0.058, rotor radius R=8.1m, hub height coefficient h=0.25, drag coefficient =0.013, number of blades of the main rotor b=4, blade mass m=79.6kg, rotor operate mode coefficient =0.3, gradient of lift a=5.65, velocity of blade top R=208m/s, distance of blade mass center coefficient xg=0.45, distance of hinge from shaft eR=0,04R, and air density at flight altitude (100m) =1.215 kg/m3. For longitudinal motion the mathematical model in vector form is: X  A X  b u 0.066   0.0509 0.1323 0.0734 0.00263  0.1344  0.1216 1.2525   0 0.3  0.3578 0.9477     X X u  0  0 0 0 1  0      0.844  12.1 0  6.512  28.329 17.88  www.intechopen.com 427 Aeronautical Engineering  X1  X  2  X  X3     X4   u1  u  u2  Equation at exit is: 1 0 Xi   0  0 0 0 0 1 0 0  X 0 1 0  0 0 1  Xi1    Xi   Xi 2   Xi 3  Matrix A of the linearized model of helicopter for lateral motion is: 0  0.15125 0.0734  0 0 1  2.948 0 A   64.42  0 0 0   55.297 0 0.413  X  A X  B u X   X5 X6 X7 X8 0 0. 3  0 0  0 1.378   0 1  0 1.64  X9  T u  u3 u4  It can be seen that matrix A is singular. From the point of automation control this means the system has unlimited number of equilibrium states. 3.1.6 Conclusion A further analysis of the mathematical model can be made in order to investigate the dynamic and static properties, and to determine control that would guaranty the object execution of the required dynamic behavior. 3.2 Modeling and simulation of spin on the vuk-T sailplane The major assumptions that were a starting point in this subject are the adopted aerodynamic conception of the aircraft and the full system of equations for the aircraft motion in the case of spin. Evaluation of aerodynamic quotients and their derivatives, and thus, the equation system as a whole, is based on defined aircraft geometric characteristics and atmospheric conditions. The method developed for simulation of the real flight situation is based on the I and II non-linear equation systems. The aerodynimc coefficients used in systems are estimated from flight information. 3.2.1 Introduction A spin is an interesting manouvre, if only for the reason that at one time there stood to its discredit a large proportion of all airplane accidents that had ever occurred. It differs from other manouvres in the fact that wings are “stalled”. Wings are beyond the critical angle-offattack, and this accounts for the lack of control which the pilot experiences over the www.intechopen.com 428 Mechanical Engineering movements of the airplane while spinning. It is form of “auto-rotation”, which means that there is a natural tendency for the airplane to rotate of its own accord. In order to reach height, a sailplane pilot has to circle inside a thermal column in constant turns of a very small radius. If for any reason (severe turbulence, pilot's error, etc.) the speed of sailplane drops below the stalling speed, at such high angles of bank the sailplane will most probably fall into a spin. Spin is a very dangerous and unpleasant maneuver. Good spin recovery characteristics are the imperative for any modern sailplane. Fig. 27. VUK-T sailplane Vuk-T sailplane (Fig. 27) is a modern single seat, all composite sailplane for advanced pilot training and competitions, designed at the Belgrade Faculty of Mechanical Engineering. Beside the other complex analyses, spin characteristics have been analyzed thoroughly. Major assumptions that were starting points in this subject are the adopted aerodynamic concepts of the sailplane and the full system of equations for the aircraft motion in case of spin. Evaluation of aerodynamic quotients and their derivatives, and thus, the equation system as a whole, is based on defined aircraft geometric characteristics and atmospheric conditions. The method developed for simulation of the real flight situation is based on the first and second non-linear equation systems. 3.2.2 Non-linear mathematical models of flight dynamics Mathematical modelling of spin motion is a very complex task and, therefore, it is necessary to introduce series of assumptions and approximations. Mathematical model is based on the first and second non-linear equations systems, shown bellow. The first system (V=const and H=const) is:   a11 cos cos  p 1   a13q  a14  a15 cos  cos  cos    a21  a22r cos   a23 p sin   a24 cos sin   a25 p  a31rq  a32   a33 p  a34r  a35  a36 r  a41pq  a42   a43 p  a44r  a45q  a46  a47 q  a51pr  a52  www.intechopen.com  a53r  a54q  a55  a56  a57 429 Aeronautical Engineering   a61r sin   a62q cos    a71p  a72r cos  tan   a73q sin  tan  The coefficients are given by: q S V S V C z   a12  1 a13  1 a14  a15   C z nk    nk  t  V 2m 2m g S V S V a21  C y    a22  1 a23  1 a24  a25  C y yk    k  t  2m V 2m Iy  Iz SlV 2 a32  C l    a31  Ix 2I x a11   a33  SlV 2 SlV 2 S V 2 SlV 2 C lp   a34  C lr   a35  C l vk    k  t  a36  C l k    k  t  2I x 2I x 2I x 2I x a41  I pp SlV 2 SlV 2 SlV 2 Iz  Ix a42  C n   a43  C np   a44  C nr   a45   2I y 2I y 2I y Iy Iy a46  a51  a56  Ix  Iy Iz a52  SlV 2 SlV 2 C n vk    k  t  a47  C n k    k  t  2I y 2I y I pp Sb V 2 SbV 2 Sb V 2 a54  C mq   a55  C m   a53  C m   2I z 2I z 2I z Iz Sb V 2 Sb V 2 C m   a57  C m nk    nk  t  a61  1 a62  1 a71  1 a72  1 a73  1 2I z 2I z The second system (V and H are not constant) is:   b11 p cos cos  V tan  V V  b12  b13q  b14  b15  b16 cos  cos  cos  V cos  V   b21V   b22r cos   b23 p sin   b24 cos cos  V   b25  b26 V V V p  b31rq  b32 V 2   b33 V 2 p  b34 V 2r  b35 V 2  b36 V 2 r  b41pq  b42 V 2   b43 V 2 p  b44 V 2r  b45q  b46 V 2  b47 V 2 q  b51pr  b52 V 2   b53r  b54 V 2q  b55q  b56 V 2  b57 V 2   b61r sin   b62q cos  www.intechopen.com 430 Mechanical Engineering   b71p  b72r cos  tan   b73q sin  tan  sin  Vr  V 2 V 2 V  b81  b82V tan   b83Vq tan   b84  b85  b86 cos  cos  cos  cos H  b91V cos  sin   b92V sin  cos cos   b93V  cos sin  and the applied coefficients are: b11   S S C z   b12  1 b13  1 b14  g b15  1 b16   C z vk    hk  t  2m 2m b21  S S C y    b22  1 b23  1 b24  g b25  1 b26  C y vk    vk  t  2m 2m b31  Iy  Iz b32  Ix b35  Sl C l vk    vk  t  b36  Sl C l    vk  t  k 2I x 2I x b41  b44  b51  Ix  Iy Iz Sl Sl Sl C l   b33  C lp   b34  C l   2I x 2I x 2I x r Sl Iz  Ix Sl b42  C n   b43  C n   2I y Iy 2I y p I p p Sl Sl Sl C nr   b45   b46  C n vk    vk  t  b47  C n k    k  t  Iz 2I y 2I y 2I y b52  I pp Sb Sb Sb Sb C m   b53  b54  C mn   b55  C m   b56  C m   2I z 2I z 2I z 2I z Iz b57  b81   Sb C m nk    nk  t  b61  1 b62  1 b71  1 b72  1 b73  1 2I z S S C x   b82  1 b83  1 b84  1 b85  g b86  C x nk    hk  t  2m 2m b91  1 b92  1 b93  1 3.2.3 The general model Modeling was carried out using the software package MATLAB. Simulation was performed by using the SIMULINK module, having the special feature to simulate a dynamic system within a graphic mode, where the linear, non-linear, time-continuous or discrete multivariable systems having concentrated parameters can be analyzed. Simulation is achieved www.intechopen.com Aeronautical Engineering 431 by creating the SIMULINK model and implementing appropriate functions for the numeric solving differential equations. 3.2.4 Program results The program has been tested on the VUK-T sailplane project developed at the Institute of Aeronautics, at Faculty of Mechanical Engineering in Belgrade. The obtained numerical results fully agree, and in some cases are complementing with the flight test results obtained for this sailplane. Fig. 28. Basic diagram (left), SE greater 40% (right) www.intechopen.com 432 Mechanical Engineering Fig. 29. Basic diagram (left), SR smaller 40% (right) In the process of sailplane analysis, influence of variation of elevator and ruder areas has been tested very effectively through the previously described model. The results obtained for two different elevator (SE) and rudder (SR) areas are presented in Fig's. 28 & 29 respectively. The results present variations of height, speed, and the characteristic angular velocities and angles of the sailplane during the first 20 seconds of spin. Practical application of this software has shown that computational design in flight mechanics gives exceptional results, allowing the designer to follow the development of calculation with a large number parameters in detail, and continuously be able to fulfill the criteria of technical demands, performing all necessary optimizations. www.intechopen.com Aeronautical Engineering 433 4. Systems, subsystems and equipments 4.1 Modeling of transient state and pressure recovering in aircraft proportional hydraulic servo-actuators In the paper is presented computer analysis of hydraulic servo-actuator transient state between its engage and full pressure recovering. In different mathematical forms hydraulic actuator dynamic model is assumed without any geometrical and physical discontinuities and ambiguity of initial pressure conditions. Hydraulic actuator can be assumed with two serial connected compressible fluid flows controlled by supply and return variable fluid flow restrictors enclosed in control servo-valve and separated by actuator piston, expressed by equivalent mass, viscous damping and arbitrary external force. Presented mathematical model includes transient state of hydraulic actuator, which can be described and determined by ambiguity of the initial pressure conditions in actuator chambers as result of the final state of its previous operations and existing fluid leakage which can produce arbitrary value of initial pressure. Initial condition of pressure is primarily caused by external force, which is arbitrary value during actuator operations. Pressure surge in the moment of change direction of piston motion as result of geometric and flow asymmetry of actuator and its control servo-valve is included in the model also. Hydraulic actuator is usually assumed with compressible fluid flow including the effects of its viscosity. Fluid compressibility is assumed as quasi-static change of its density depending of static pressure. Each of the mentioned effects produces local pressure drop and surge and corresponding actuator operational time delay, which is the limiting factor of its cyclic velocity. Dominant influence on actuators time delay (less than 3% of unit step discrete control piston stroke time) is caused by fluid volumetric compressibility. In the paper following problems are treated: pressure discontinues changes in reverse of piston motion direction; nonlinear effects of actuator behavior; actuator dynamic model linearization. Any direction change of actuator motion produces pressure discontinuity in its source and return pipelines. This is caused by inversion of fluid flow which produce connection change between supply pipeline and actuator chambers. In the moment of fluid flow direction change each of actuator chambers inter-change connections with system pump and return pipeline and produce corresponding discrete change of pressure in actuator chambers. Possible pressure drop or surge is also caused by geometric asymmetry of servo valve. 4.1.1 Effects of system discontinuities Any change of direction of actuator motion produces pressure discontinuity in its both pipelines. This discontinuity is caused by inversion of fluid flow which produce the change of connection between supply pipeline and both actuator chambers. In the moment of change of fluid flow direction each of actuator chambers change connections with system pump and return pipeline, producing corresponding discrete changes of pressure in actuator chambers. Possible pressure drop or surge is also caused by geometric asymmetry of servo valve. These effects are explained on the following figures. On figure 30 is shown actuator motion asymmetry between direct and reverse modes. This asymmetry is result of pressure distribution along supply and return streamlines, shown on diagrams a) and c) on figure 32. Diagram a) corresponds to direct mode of actuator function www.intechopen.com 434 Mechanical Engineering represented by symmetrical pressure drops at supply and return branches of its servo-valve. The third step of pressure drop correspond to the actuator piston position and corresponds to the applied external force. On diagram c) is shown pressure distribution for reverse actuator mode. Main difference between these modes is in opposite directions of external force related to the streamline of fluid flow. For reverse mode external force support system pump like additional serial connected system sources. This fact appears on figure 30 as different curve gradient for direct and reverse modes. However, absolute value of gradient is greater for reverse mode. More data about gradient value will be shown corresponding to the figure 35. On figure 31 is presented actuator output on servo-valve control input assumed as transient step unit function. This approximation is very close to the real situation for digitally controlled actuators. On diagrams b) and d) of figure 32 are shown, respectively, equivalent pressure drops for direct and reverse actuator modes corresponding to the conventional mathematical modeling of hydraulic actuator with total pressure drop on servo-valve by neglecting effects at its supply and return parts as two separated fluid flows. This usual approximation cannot be accepted if in system are included effects of hydraulic pressure drop and surge caused by fluid compressibility. Fig. 30. Fig. 31. www.intechopen.com 435 Aeronautical Engineering Fig. 32. 4.1.2 Conventional system modeling Real state is described by pressure drops at supply and return parts of control servo valve and corresponding pressure difference caused by external force: pul  ps  pa piz  pr  p0 pa  pr  F Ak (15) sgn(Fxr ) Function sgn denotes both directions of external load action, represented as direct and reverse modes of actuator function. Equations (15) defines basic formulation of system dynamic model. Closed formulation of mentioned pressure drops cannot be determined without additional approximations. If we assume that fluid flow through servo-valve is a turbulent, approximate expressions of corresponding equivalent pressure differences in accordance to the figures b) and d) are defined for incompressible fluid flow by following relations: . Qs  Q0  Ak x k    br xr where pr'  2  . ( ps  p'a ) Qs  Q0   Ak xk   br xr 2  ( pr'  p0 ) (16) ka k x k  ps and p'a  a x k  p0 Ak Ak 4.1.3 Separate flow modeling Previous relations (15) can be expanded for approximately symmetric supply and return flow characteristics of servo-valve in the following form (with assumed value of p0=0) for direct and reverse modes: www.intechopen.com 436 Mechanical Engineering p  p0 1 F F 1 ) psr  s ) pr  p0  p  ( ps  pa  ps  p  ( ps  2 2 2 Ak Ak (17) Expressions (16) needs more attention on its meaningless. Corollary of expressions (16) is that nominal system pressure for zero external load is defined as (17).It means that hydraulic system pressure for zero load must be equal to the psr. In the cases of continuous actuator function previous relation holds. In addition, this expression are satisfied for all regimes in which effects of pressure surge can be neglected. In opposite cases, presented mathematical system formulation does not hold. Then subjected mathematical model is not compatible. Corresponding to the relations (16) pressure drop can be determined in the form (18), where pF   F / Ak represent pressure drop corresponding to external load. . Qs  Q0  Ak x k    br xr 2  . ( ps  p0  p ) Qs  Q0   Ak xk    br xr 2  ( ps  p0  p ) (18) Finally, static pressure in supply and return branches of actuator streamline can be expressed in expanded form:   F F 1 1 sgn xr  pr   ps  p0  sgn xr  pa   ps  p0  2 2 Ak Ak   (19) Relations (19) are similar to the relations (16). In previous discussion hydraulic system pump is assumed as strong one. It means that the system pump is able to takes up system pressure to the maximal nominal value. This assumption is valid except for existence of system model incompatibilities. This problem can be solved by assuming system pump as a weak one at the initial moment of actuator engage. It follows that any regime of small external load must be assumed as of weak pump. This statement arise from the fact that static pressure in hydraulic system at any moment of its function is caused by external load and pressure loses. As consequence of previous statements, corresponding boundary conditions at actuator pipeline inlet must be determined at initial moment as maximal pump flow. Caused value of system static pressure exists till the moment of pressure upgrading to its nominal system value. In that moment boundary conditions changes to the determined inlet pressure (equal maximal nominal value) and caused value of fluid flow ( second part of flow cross the relive valve), which corresponds to the strong system pump. If reduced valve is built in hydraulic system at actuator supply branch mentioned effects decreases. But in initial moment of actuator engage they can’t vanish completely. It means that each pump can’t be strong one for the whole possible regimes, spatially at its initial moment. 4.1.4 Boundary conditions Corresponding actuator block diagrams for the cases of weak and strong system pump are presented on figure (4) and relations (20) or (21). . Ak x k    br xr or in expanded form www.intechopen.com 2  . ( pp max  p ) Ak x k    br xr 2  ( pp max  p ) (20) 437 Aeronautical Engineering . Ak x k   sbrs xr 2  ( pp max  F     s s sgn xr )  br   br   br Ak (21) where index s denotes parameters of symmetric servo-valve. Relation (21) is final mathematical model form of control servo-valve pressure drop for the case of actuator symmetry. This formulation is well known. It must be noted that ppmax is correct term only for strong system pump. For the cases of weak pump ppmax becomes equal ps, with corresponding changes of boundary conditions formulation, which gives following model: Qp max  Ak x k   sbrs xr 2  ( pp max  F sgn xr ) Ak (22) Fig. 33. 4.1.5 Analysis of incompressible models Corresponding to discussion about strong and weak regimes of pump function, it is of interest to determine actuator abilities for suppression external load. Corresponding coefficient  is defined as ratio between maximal external acting load and maximal possible load which can be suppressed by the system pump. It is also shown correlation between nonlinear and formed linear actuator models with incompressible fluid flow. Effects of model nonlinearities and its corresponding linearisation for the case of incompressible flow are presented on figure 34. We can conclude from diagram that for usual reserve of actuator ability these effects are practically similar. This conclusion indicates that other effects are of higher influence than the effects of model linearisation. On figure 35 is compared nonlinear and linearised approximation on the whole domain of actuator piston stroke. Corresponding differences of actuator output between these models are presented on the figure 36. www.intechopen.com 438 Mechanical Engineering Fig. 34. Fig. 35. Fig. 36. Fig. 37. Fig. 38. Fig. 39. On figures 38 and 39 are shown diagrams of non dimensional coefficients of model linearisation corresponding to its nominal regimes. 4.1.6 Actuator transfer characteristics In this paper actuator simulation is presented for step and frequent harmonic input of actuator control servo-valve. Step input is used for actuator quasi stationary characteristics and parameters identification. Harmonic impute can be used for high cyclic actuator identification. www.intechopen.com Aeronautical Engineering 439 On diagram of figure 40 is presented computer simulation of nonlinear periodic actuator output piston stroke for harmonic control input, both in non dimensional form. If corresponding relative amplitude of servo-valve throttle is equal 0.50 depending of relative time, equal 1 for total piston stroke for its maximal possible velocity. It is easy to see that approximate harmonic output is recovered approximately after 5 cycles. System nonlinearity is shown on diagram of figure 41, which represents inverse simulation of system relative control input for harmonic relative output. Fig. 40. Fig. 41. Fig. 42. www.intechopen.com 440 Mechanical Engineering On diagram of figure 42 are presented simulation curves of relative actuator piston stroke position and its corresponding relative velocity as actuator system outputs for unit step control relative input if exists only inertial actuator load including equivalent total mass of aerodynamic control surfaces. Corresponding actuating relative time delay is less than 0,004. 4.1.7 Actuator modeling with assumed quasi-static fluid compressibility Pressure drop in hydraulic systems can be caused by small external load or local increasing of fluid flow to be greater than maximal possible pump source flow. To prevent this it is suggested to separate corresponding branch of actuator supply by corresponding reduction valve. In these cases possible pressure surge are not of high influence and is determined by equivalent actuator stiffness together with potential external load. Pressure increases proportionally with increasing of piston displacement corresponding to its velocity. This pressure increasing is too slower than for the cases of pressure surge caused by fluid compressibility. If actuator is assumed with quasi static compressible fluid flow, system model can be presented for symmetric supply and return branch of fluid flow in the following form:  ()br () (  xr )  ()br () (  xr ) 2  2  . . . . ( ps  pa )   Ak x k   Ak ( H cil  x k ) pa  c( pa  pr ) ( pr  p0 )   Ak x k   Ak ( H cil  x k ) pa  c( pa  pr ) (23) F  Ak ( pa  pr ) where are: flow coefficient, br equivalent geometric wide, xr position of control valve throttle, ps supply pressure of hydraulic system pump, pa static pressure in supply chamber of actuator cylinder, pr static pressure in return chamber of actuator cylinder, p0 static pressure in return pipeline, Ak area of actuator piston, coefficient of fluid compressibility, c coefficient of fluid leakage, Hcil piston stroke, coefficient of parasite volume of connected pipeline to actuator cylinder, xk position of piston and F applied external force (including inertial forces) to the actuator piston. Both signs in equations corresponds to the direct and reverse modes of actuator function. Fig. 43. www.intechopen.com Aeronautical Engineering 441 Fig. 44. On figure 43 are presented simulation of actuator piston relative stroke and relative static pressure in supply and return actuator chambers for usual mathematical form of actuator dynamic model. As it is explained in abstract, initial relative pressure values (equal 0.5) can exists in ideal model only. In real cases, initial values of relative pressure in actuator chambers is the result of actuator history and fluid leakage, which produces it’s ambiguity. Extreme case corresponds to the zero values of initial relative pressures. Corresponding simulation of supply and return relative pressures are presented on figures 45 and 46. Fig. 45. Fig. 46. www.intechopen.com 442 Mechanical Engineering 4.1.8 Conclusion All of the exposed diagrams are related to non dimensional ratio system coordinates, where are: k1 piston position, pa static pressure in actuator supply chamber, pr static pressure in actuator return chamber, r control servo valve throttle position, ratio of power reserve corresponding to applied external load. Presented system model enables its compatibility corresponding to the various initial conditions. On figure 14 are shown system simulation for incompressible fluid flow and corresponding model non compatibility of initial conditions. Possible pressure difference between supply and return actuator chambers which is not compensated by external force produces piston “shock” motion, which can not be described by incompressible flow system modeling. Pressure difference can be caused by various effects which produces fast changes of fluid static pressure and/or external load. Actuator locked position for longer time period also is the reason for described effects. Pressure drop or surge caused by fluid compressibility and initial condition discontinuity as result of closed control servo-valve throttle position are shown on the diagrams on figures 16 and 17. Piston position difference in relation with its incompressible model motion is defined on figure 44. Initial piston acceleration produces practical piston “shock” motion , expressed in the later as increasing static error of its position (less than 1% for usual types of fluids). Mentioned effects are of high interest for digitally driven actuators. Supply and return pressure surge is presented on figure 45. Supply and return pressure drop are presented on figure 46. 5. References Bengin A., Mitrović, Č., Cvetković, D., Bekrić, D., Pešić, S.,: Improved Solution Approach for Aerodynamics Loads of Helicopter Rotor in Forward Flight, Journal of Mechanical Engineering 54(2008)3, 170-178, UDK 533.661 , Cvetković, D., Mitrović, Č., Bengin A., Kostić, I.,: Mathematical Models of Helicopter Flight Dynamics, AIAA 2002-0529, (CA 94035-1000) 40th AIAA Aerospace Sciences Meeting & Exhibit 14–17 January 2002, Reno, Nevada, USA Green, W. L.,: Aircraft hydraulic system, John Wiley and sons, UK, 1985, Janković, J.,: Computer Analaysis and Simulation of Transient State and Pressure Recovering in Fast Cyclie Hydraulic Actuators, Proccedings of ICAS – 96, Sorento, Italy, 1996. Mitrović, Č., Cvetković, D., Bengin A., Bekrić, D.,: An Optimal Main Helicopter Rotor Projection Model Obtained by Viscous Effects and Unsteady Lift Simulation, Journal of Mechanical Engineering , (2010), vol. 56 / 6, 357-367, Mitrović, Č., Cvetković, D., Bengin A., Bekrić, D.,:Two Dimensional Airfoil in Incompresibile Inviscid Flow, 2nd Ankara International Aerospace Conference & Symposia AIA ‘98, Ankara, 1998. Mitrović, Č., Cvetković, D., Kostić, I., Subić ,A.,:Safety Aspects in Sport Flying: Analysis of Spin of a Training Sailplane, The 4th International Conference on the Engineering of Sport, Ref. No O_142, 3-6 September 2002, Kyoto, Japan Napomena: međunarodni, objavljen Rudinger, G.,: Wave Diagrams for Nonsteady Flow in Ducts, D. Van Nostrand company, Canada 1955. www.intechopen.com Mechanical Engineering Edited by Dr. Murat Gokcek ISBN 978-953-51-0505-3 Hard cover, 670 pages Publisher InTech Published online 11, April, 2012 Published in print edition April, 2012 The book substantially offers the latest progresses about the important topics of the "Mechanical Engineering" to readers. It includes twenty-eight excellent studies prepared using state-of-art methodologies by professional researchers from different countries. The sections in the book comprise of the following titles: power transmission system, manufacturing processes and system analysis, thermo-fluid systems, simulations and computer applications, and new approaches in mechanical engineering education and organization systems. How to reference In order to correctly reference this scholarly work, feel free to copy and paste the following: Časlav Mitrović, Aleksandar Bengin, Nebojša Petrović and Jovan Janković (2012). Aeronautical Engineering, Mechanical Engineering, Dr. Murat Gokcek (Ed.), ISBN: 978-953-51-0505-3, InTech, Available from: http://www.intechopen.com/books/mechanical-engineering/aeronautical-engineering InTech Europe InTech China University Campus STeP Ri Slavka Krautzeka 83/A 51000 Rijeka, Croatia Phone: +385 (51) 770 447 Fax: +385 (51) 686 166 www.intechopen.com Unit 405, Office Block, Hotel Equatorial Shanghai No.65, Yan An Road (West), Shanghai, 200040, China Phone: +86-21-62489820 Fax: +86-21-62489821 © 2012 The Author(s). Licensee IntechOpen. This is an open access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.