Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 44, 4, pp. 881-906, Warsaw 2006 TILTROTOR MODELLING FOR SIMULATION IN VARIOUS FLIGHT CONDITIONS Marek Miller Institute of Aviation, Warsaw e-mail: m.miller@ilot.edu.pl Janusz Narkiewicz Institute of Aeronautics and Applied Mechanics, Warsaw University of Technology e-mail: jnark@meil.pw.edu.pl Acomputermodel of a tilt-rotor has been developed for calculating perfor- mance, simulating flight and investigating stability and control. Themodel is composed of a fuselage, wings, an empennage, engine nacelles and ro- tors. Tiltrotor equations of motion have been obtained by summing up inertia, gravity and aerodynamic loads acting on each part of the aircraft. Aerodynamic loads at wings, empennage and rotor blades have been cal- culated using a quasisteadymodel. For rotor induced velocity, the Glauert model has been used. The influence of the rotor inflow wing and empen- nage aerodynamic loads has been found using the actual value of induced velocity. The computer program of tilt-rotor model has been developed in theMatLab environment.The sub-programs for load calculation have been supplemented by modules for calculation of trim states and stability and control matrices. During the first stage of model investigation, steady fli- ght conditionswere calculated,which gave insight into rotorcraftbehaviour andmodel quality. Key words: tiltrotor, rotorcraft, flight control, simulation Notations Indicies a – aerodynamic b – inertia c – related to three-dimensional body (fuselage, nacelle) f – fuselage 882 M. Miller, J. Narkiewicz g,G – gravity h – horizontal stabilizer i – inertia I =1 – element (wing, nacelle, rotor) on right side of aircraft I =2 – element (wing, nacelle, rotor) on left side of aircraft J =1,2,3 – rotor blade index of i-th rotor m – aerodynamic moment, moving element, mass n – nacelle p – aircraft, fuselage r – rotor Matrices and vectors A – rotation matrix of moving element, A=A(φ,ϕ,γ) AG – aircraft rotation matrix Ar – transformation matrix of coordinate system in airfoil AV – velocity matrix C – control vector, C= [τ,δw,δh,δv,Θ] ⊤ Ix – inertia matrix of x-th element k() – coefficients of inflow due rotors and wings Qxy – y-th loads of x-th element, Qxy = [Fxy,Mxy] ⊤ V – vector of linear velocity in Opxpypzp, V = [U,V,W ] ⊤ V a – section velocity V c – vector of linear velocity of three-dimensional body U – induced velocity in Opxpypzp W – wind velocity in Opxpypzp ∆W() – inflow due rotors and wings in Opxpypzp X – vector of aircraft motion, X = [V ,Ω,xg,Φ] ⊤ = [U,V,W,P,Q,R,xg,yg,zg,Φ,Θ,Ψ] ⊤ Y x – state vector of x-th element g – vector of gravity acceleration gp – vector of gravity acceleration in Opxpypzp rCG – position of C.G. in Opxpypzp rnCG – position of C.G. in moving element coordinate system rn – position of moving element coordinate system in Opxpypzp xg – translation of aircraft relative to ground, xg = [xg,yg,zg] ⊤ Ω – vector of angular velocity in Opxpypzp, Ω= [P,Q,R] ⊤ Ωp – velocity matrix of plane (non-moving element) Ωm – velocity matrix of moving element Φ – Euler’s angles described as vector, Φ= [Φ,Θ,Ψ]⊤ Θ – control of rotor swash-plates, Θ= [Θ01,Θ02,Θ11,Θ12,Θ21,Θ22] ⊤ Tiltrotor modelling for simulation... 883 δw – inclination angles of wing flaps (ailerons), δw = [δw11,δw12,δw21,δw22] ⊤ δh – inclination angles of elevators, δh = [δh1,δh2,δh3] ⊤ δv – inclination angles of rudders, δv = [δv1,δv2] ⊤ ω – angular velocity of rotors, ω= [ω1,ω2] ⊤ τ – angle of nacelle tilt Coefficients Ar – characteristic area Cx(α,δ) – drag coefficient (in 2D flow) Cz(α,δ) – lift coefficient (in 2D flow) Cmy(α,δ) – aerodynamic pitch moment coefficient (in 2D flow) Cxc(α,δ) – drag coefficient (in 3D flow) Cyc(α,δ) – side force coefficient (in 3D flow) Czc(α,δ) – lift coefficient (in 3D flow) Cmxc(α,δ) – aerodynamic bankmoment coefficient (in 3D flow) Cmyc(α,δ) – aerodynamic pitch moment coefficient (in 3D flow) Cmzc(α,δ) – aerodynamic yaw moment coefficient (in 3D flow) Rr – characteristic length Va – modulus of motion velocity of section Va = |V a| c(y) – chord (characteristic dimension) g – gravity acceleration mx – mass of x-th element α – angle of attack β – slip angle δ – angle of inclination of steering element ρ – air density 1. Introduction In the recent years, increasing interest in tilt rotor technology is observed among rotorcraft community. The development of V-22 aircraft spurred a gre- at research and development effort in tilt rotor aerodynamics (McVeigh et al., 2004), acoustics (Polak and George, 1998), flight control (Weakley et al., 2003), aeroelasticity (Nixon et al., 2003), etc. The next tilt rotor project, Bell- AgustaBA 609, initiated a tilt-rotor design effort in Europe. Recently, several European research projects concerning tilt rotor technology are going on or have already been completed (Cicale, 2003), and some new initiatives appe- ared. Despite the expanding interest in development and applications of the 884 M. Miller, J. Narkiewicz tilt-rotor, there are notmanypapers published in generally accessed literature presenting a comprehensive approach to the tiltrotormodelling and simulation (Polak and George, 1998; McVicar and Bradley, 1992; Srinivas and Chopra, 1996). It was a background for developing atWarsaw University of Technolo- gy, a computermodel of the tiltrotor for flight simulation and analysis of trim, stability and control in various flight conditions. 2. General model of mobile objects Atiltrotor computermodelwas constructed using the genericmodel ofmobile objects, developed at Warsaw University of Technology, to simplify model derivation and simulation of motion for various air, water and land vehicles. This software in the MatLab environment is composed of a main program and sub-programmes performing computations frequently done in dynamics of airplanes, helicopters, sea vessels, etc. The possibility of the general approach stems from the fact that inertia, gravity and aerodynamic/hydrodynamic loads that act on the vehicles con- sidered are described by the same formulae in reasonably chosen coordinate systems. As a consequence, using this software, computer models of six de- grees of freedom motion of various vehicles may be prepared in a fast and efficient way by a proper selection of systems of coordinates and application of supporting subroutines. Within this software, a vehicle ismodelled as the base part (”fuselage”) to which other parts/elements of it are joined. These parts may be fixed to the fuselage or may rotate and/or translate relative to it. Equations of motion are derived in the main system of coordinates fixed to the fuselage by summing loads acting on all elements of the vehicle. These loads are calculated in local coordinate systems fixed to the elements and then transformed to the system of coordinates of the fuselage. It is easy to change the number of elements of a model by adding new or deleting the existing elements and enriching the methods of calculating loads by application of various methods. The computer program for modelling vehicle dynamics is constructed in the samemodular way, by applying subroutines for calculating each type of loads and their transformation to the main ”fuselage” system of coordinates. Various subprograms were developed for frequently performed operations like calculation of inertia matrices, flow velocities in different coordinate sys- tems, angles of attack, slip angles, table-look procedures for aerodynamic co- efficients, etc. The generic model was used for the tilt rotor and sea vessel analysis. Tiltrotor modelling for simulation... 885 3. Tiltrotor model The tiltrotor model investigated in this study was based on V-22 Osprey air- craft. Due to the modular structure of the software, it is possible to adjust themodel to the concepts of tilt-wing, a partly tiltingwing or other rotorcraft with different empennage, wings, nacelles, etc. The tiltrotor (Fig.1) model is composed of a fuselage, two wings with two trailing edge flaps at each wing, nacelles mounted at the wings tips, rotors mounted in front of the nacelles and a horizontal stabilizer with three flaps and two vertical fins with one rudder mounted on each fin. All components of the aircraft are rigid. The rotors have three blades fixed to the shaft by a pitch bearing. The pitch angles of rotor blades are controlled by swash-plates for constant (collective) and harmonic (cyclic) components. Fig. 1. Configuration of the tiltrotor The tiltrotor equations of motion are derived in the aircraft coordinate system Opxpypzp fixed to the fuselage (Fig.2). The centre Op of the aircraft system is placed in the point, where rotor shafts intersect with the fuselage plane of symmetry. The Opxp axis lies in the longitudinal symmetry plane of the fuselage, parallel to the horizontal reference plane of the aircraft and is directed ”to the cockpit of the aircraft”. The Opzp axis lies in the longitudinal plane of symmetry of the fuselage and is directed down to the undercarriage. The Opyp axis is pointing right, while looking along the Opxp axis. The two other systems of coordinates important for simulation of aircraft motion (Fig.2) are the ground coordinate system Ogxgygzg fixed to the ground and the gravity coordinate system OGxGyGzG related to the gravity accelera- tion. The centre of the ground coordinate system Og may be placed in any selected point in space (for instance: a selected point on the airfield). The aircraft motion is simulated from this point. The Ogzg axis is vertical, it co- incidences with the direction of the Earth gravity acceleration and is directed 886 M. Miller, J. Narkiewicz Fig. 2. Systems of coordinates of an aircraft along its positive value. The axes Ogxg and Ogyg lie in the horizontal plane, i.e. the plane perpendicular to the direction of gravity acceleration. The Ogxg lies along the local geographical meridian pointing north, and the Ogyg axis points east. The centre OG of the gravity system of coordinates coincides with the centre of fuselage Op. The axis of the gravity system is parallel to the axis of the ground coordinate system. The gravity coordinate system is used for describing the position and attitude of the aircraft. The gravity system is translated fromthe inertia systembythevector xG(t)= [xG(t),yG(t),zG(t)] ⊤, which is a function of time. The vector xG(t) describes the translation of the tiltrotor in space. It should be noted here that the point OG may not be the centre of gravity of the aircraft. The relation between the gravity and the aircraft system of coordinates is described by Euler’s angles of rotation. The relation between the coordinates in both systems is given by an equation xp =AG(Ψ,Θ,Φ)xG (3.1) where the rotation matrix AG has the form AG = (3.2) =    cosΨ cosΘ sinΨ cosΘ −sinΘ cosΨ sinΘsinΦ− sinΨ cosΦ sinΨ sinΘsinΦ+cosΨ cosΦ cosΘ sinΦ cosΨ sinΘcosΦ+sinΨ sinΦ sinΨ sinΘcosΦ− cosΨ sinΦ cosΘcosΦ    Aircraft motion is described by the vector X = [V ,Ω,xg,Φ] ⊤ = [U,V,W,P,Q,R,xg,yg,zg,Φ,Θ,Ψ] ⊤ (3.3) as a composition of four vectors: Tiltrotor modelling for simulation... 887 – translation velocity V = [U,V,W ]⊤ – rotation (rates) Ω= [P,Q,R]⊤ – Euler’s angles written in a vector form Φ= [Φ,Θ,Ψ]⊤ – translation of the aircraft relative to the ground xg = [xg,yg,zg] ⊤. 3.1. Equation of motion Equations of motion of the aircraft are derived using d’Alambert princi- ple, summing up at the point Op all loads (forces and moments) acting on the fuselage, wings, control surfaces, nacelles, and rotors. The system of six equations ofmotion is obtained, whichmay be grouped as two subsystems for forces and moments acting on the fuselage and wings (index p), two nacelles (indices n1, n2) and two rotors (indices r1, r2) Fp+Fn1+Fn2+Fr1+Fr2 =0 (3.4) Mp+Mn1+Mn2+Mr1+Mr2 =0 Each element of the above equations consists of inertia, aerodynamic and gravity parts Q() = [F (),M()] ⊤ =Q()i+Q()a+Q()g (3.5) 3.2. Inertia loads The expression for inertia loads is obtained from the conservation of mo- mentum, which after performing some mathematical manipulations, may be written in the matrix form Qpb = IpẎ p+ΩpIpY p (3.6) where the tilt-rotor state vector Y p = [U,V,W,P,Q,R] ⊤ is composed of com- ponents of aircraft translation velocity and rates. The inertia matrix has the form Ip =          mp 0 0 0 Szp −Syp 0 mp 0 −Szp 0 Sxp 0 0 mp Syp −Sxp 0 0 −Szp Syp Ixp −Ixyp −Ixzp Szp 0 −Sxp −Ixyp Iyp −Iyzp −Syp Sxp 0 −Ixzp −Iyzp Izp          (3.7) 888 M. Miller, J. Narkiewicz and the velocity matrix is Ωp =          0 −R Q 0 0 0 R 0 P 0 0 0 −Q P 0 0 0 0 0 −W V 0 −R Q W 0 −U R 0 P −V U 0 −Q P 0          (3.8) Expression (3.6) describes inertia loads acting on the fuselage and on parts of the rotorcraft fixed to it, i.e. wings and control surfaces. For parts of the aircraft rotating relative to the fuselage, i.e. nacelles (Fig.3) and rotors, the inertia loads are calculated as Qmb = Im(Ẏ p+ Ẏm)+(Ωp+Ωm)Im(Y p+Ym) (3.9) In (3.9), the inertia matrices of rotating elements are calculated in the fuselage system of coordinates. Fig. 3. Nacelle coordinate systems Additionalvelocitymatrices Ωm areadded, containing inageneral case the rates and velocities of these elements relative to the fuselage. In the tiltrotor case, they have the form Tiltrotor modelling for simulation... 889 Im =          mm 0 0 0 Szm −Sym 0 mm 0 −Szm 0 Sxm 0 0 mm Sym −Sxm 0 0 −Szm Sym Ixm −Ixym −Ixzm Szm 0 −Sxm −Ixym Iym −Iyzm −Sym Sxm 0 −Ixzm −Iyzm Izm          (3.10) Ωm =          0 −Rm Qm 0 0 0 Rm 0 −Pm 0 0 0 −Qm Pm 0 0 0 0 0 0 0 0 −Rm Qm 0 0 0 Rm 0 −Pm 0 0 0 −Qm Pm 0          The nonlinear parts of equation (3.9) contain all accelerations acting on the rotating elements, including gyroscopic effects. 3.3. Gravity loads Thegravity forces andmoments are calculated first in the centres of gravity of the fuselage and aircraft elements. Next, they are transformed to the cen- tre Op of the fuselage systemof coordinates.The vector of gravity acceleration in the gravity system of coordinates has components g= [0,0,g]⊤ (3.11) It is rotatedwith respect to the fuselage systemof coordinates using the trans- formation gp =AG(Ψ,Θ,Φ)g (3.12) Themasses of fuselage, wings and empennage are accounted for together and the gravity loads of these parts are calculated as fuselage gravity loads Fpg = mpgp = mpAGg (3.13) Mpg = rCG×Fpg = rCG× (mpAGg) where the vector rCG describes the position of C.G. of the fusela- ge/wings/empennage relative to Op. The positions of C.G. of other parts of the airplane in the local systems of coordinates are calculated as rCG = rn+A(φ,θ,ψ)rnCG (3.14) where: rn is the vector of C.G. of the given element relative to the fuselage centre, A(φ,θ,ψ) – general description of the matrix of rotation of the local 890 M. Miller, J. Narkiewicz system of coordinates (fixed to the element) relative to the fuselage system of coordinates (it should be defined separately for each rotating element of the tilt-rotor), rnCG – vector of the position of CG of the element in the local system of coordinates. The gravity loads acting on other aircraft elements are transferred to the fuselage system of coordinates by using formulae Fng = mngp = mnAGg (3.15) Mng = [rn+A(φ,θ,ψ)rnCG]×Fng = [rn+A(φ,θ,ψ)rnCG]× [mnAGg] 3.4. Aerodynamics loads In the generic model, two methods of calculation of quasi-steady aerody- namic loads are used: a 2Dmodel for elongated elements composed of airfoils and a 3Dmodel for solid bodies. Fig. 4. Two-dimensional flowmodel The two dimensional flow (Fig.4) is used for wings, rotor blades and hori- zontal andvertical stabilizers. For better description, some coordinate systems are introduced: – coordinate system of the mounting point of the element ONxNyNzN – local airfoil coordinate systems Oaxayaza – coordinate system fixed to the airfoil aerodynamic centre Oacxacyaczac – coordinate system connected with the geometrically twisted airfoil Oagxagyagzag – coordinate system defined by the vector of velocity of the local airfoil Oavxavyavzav – coordinate system of local inflow of the airfoil OVxVyVzV . Tiltrotor modelling for simulation... 891 The transformation of thementioned above coordinate systems can be descri- bed as xN = ry+rac+AagAavAaVxV = ry +rac+ArxV where Ar =AagAavAaV . The elongated element is divided along the span for cross sections inwhich the two dimensional flow is assumed. In each cross section, the instant total flow velocity (Fig.4) is calculated as V ()i =A −1 r() A −1(V +Ω×r()−AGW +U()+∆W()) (3.16) taking into account the velocities of: a) motion of the fuselage and their aircraft elements in the inertia coordi- nate system b) motion of the air W due to wind, gusts, etc. c) inflow due rotors and wings. Ar() matrix is a substitute matrix describing rotations of the coordinate sys- tems in airfoil (Fig.4) from the systemof themounting point of the element to the system of local velocity inflow (Ar =AagAavAaV ). The inflow due rotors and wings ∆W () is calculated as ∆W () = k()(V +Ω×r()−AGW +U()) (3.17) where k() is a coefficient assumed for a given element. On rotor blades, the rotor induced velocity U() obtained from Glauert’s model is also taken into account while calculating flow velocity. Aerodynamic loads acting in the ”aerodynamic centre” (AC) of the cross sections are cal- culated using aerodynamic coefficients for the airfoil section obtained for the actual airfoil angle of attack α and deflection δ of flaps (if they exist) by a table look procedure. In the cross sections, the vectors of aerodynamic forces andmoments are calculated in the flow coordinate system as dP = [dD,0,dL]⊤ dM = [0,dM,0]⊤ (3.18) where: — drag dD = 1 2 ρc(y)V 2a Cx(α,δ)dy (3.19) — lift dL = 1 2 ρc(y)V 2a Cz(α,δ)dy (3.20) —moment dM = 1 2 ρc2(y)V 2a Cmy(α,δ)dy (3.21) 892 M. Miller, J. Narkiewicz Fig. 5. Three-dimensional airflow The loads are transferred to theelement local systemof coordinates, integrated along the span and transferred to the fuselage system of coordinates. The fuselage and nacelles (Fig.5) are treated as three-dimensional bodies. Aerodynamic loads, as in the 2D case, are calculated in the centre of local, flow coordinate system P()a = 1 2 ρArV 2C()c(α,β) (3.22) M()a = 1 2 ρArRrV 2Cm()c(α,β) using the local instant velocity V c = [Uc,Vc,Wc] ⊤ in the body aerodynamic centre. The angle of incidence α and angle of slip β are calculated as α = asin Wc √ U2c +V 2 c +W 2 c β = asin Vc √ U2c +V 2 c (3.23) and the table lookprocedure is used for obtaining theaerodynamic coefficients. The loads from the flow system of coordinates are transformed to the element system of coordinates using the rotation matrix AV =    −cosβcosα −sinβ cosβ sinα −sinβcosα cosβ sinβ sinα −sinα 0 −cosα    (3.24) Finally, as in the 2D case, the loads are transformed from the element local system of coordinates to the fuselage system of coordinates. Tiltrotor modelling for simulation... 893 4. Details of modelling tiltrotor parts 4.1. Fuselage The inertia loads acting on the fuselage are calculated from formula (3.6), gravity loads from (3.13) and the aerodynamic loads from (3.22). The iner- tia matrix covers the fuselage, wings and empennage – elements fixed to the fuselage. 4.2. Wings The tiltrotor wings have prescribed planeform, twist and airfoil distribu- tion along the span. At each wing, there are two flaps (ailerons) controlled individually. Flap deflection angles δw = [δw11,δw12,δw21,δw22] ⊤ form a part of control variables in the simulation of aircraft motion. The aerodynamic lo- ads are calculated using the 2D model. The induced velocities of rotors are included into velocity of wing airflow calculated in sections along the span. 4.3. Horizontal stabilizer The horizontal stabilizer forms a part of the empennage and is mounted at the fuselage tail. There are three flaps (elevators) at the trailing edge of the stabiliser, controlled individually. Thehorizontal stabilizermay have arbitrary airfoil shape and twist angle distribution along the span. The aerodynamic loads are calculated using the 2D model. Influence of the induced velocity of rotors is taken into account in the stabilizer sections along the stabilizer span with the time delay due to the distance travelled from the rotor to empennage. The elevator deflection angles, written as the vector δh = [δh1,δh2,δh3] ⊤ are a part of the vector of control variables. 4.4. Vertical stabilizers The vertical stabilizers aremounted at the tip of the horizontal stabilizer. There is one flap (rudder) at the trailing edge of each of the two stabilizers, controlled individually. The vertical stabilizers may have arbitrary airfoil sha- pe and twist angle distributions along the span. The aerodynamic loads are calculated using the 2Dmodel. Influence of the induced velocity of the rotors is taken into account in the proper sections along the stabilizer span with the time delay due to the distance travelled from the rotor to empennage. The rudder deflection angles δv = [δv1,δv2] ⊤ are a part of the vector of control variables. 894 M. Miller, J. Narkiewicz 4.5. Engine nacelles The engine nacelles are placed at the tip of each wing. They may rotate about the axis perpendicular to the fuselage plane of symmetry. The nacelles inertia loads are calculated using expression (3.10) and gravity loads using expression (3.15). Aerodynamic loads are calculated using the 3Dmodel. The velocity of airflow around the nacelles contains inflow from the rotors and a component due to rotation of nacelles about the fuselage axis. The angle of nacelle rotation τ is included into the set of control variables. 4.6. Rotors The rotors rotation axes are perpendicular to the axis of nacelle rotation relative to the fuselage.When the rotor axes are in the horizontal position, the right rotor rotates clockwise, and the left counter-clockwise looking from the rear of the fuselage. For calculations of inertia loads, the rotors are treated as rotating discs and in the hub system of coordinates, the inertia matrices have the form Iri =          mri 0 0 0 0 0 0 mri 0 0 0 0 0 0 mri 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Izri          (4.1) The final form of inertia matrices of rotors (3.10)1 results from nacelle and rotors rotation about the fuselage axis. Each rotor has three blades mounted to the shaft by pitch bearings. The pitch of the rotor blades is controlled by the swash-plate, resulting in collective and periodic control in the form θij = θ0i+θ1icos ( Ωit+ 2π 3 j ) +θ2i sin ( Ωit+ 2π 3 j ) (4.2) The blades may have arbitrary planeform, twist and airfoil distribution along the span. The aerodynamic loads are calculated using strip theory with qu- asisteady aerodynamic loads using the table-look procedure for calculations of aerodynamic coefficient of airfoils (2D model described above). The induced flow is calculated usingGlauert’smodel.The control of rotor swash-platesmay be written in the vector form Θ = [Θ01,Θ02,Θ11,Θ12,Θ21,Θ22] ⊤ containing the sub-set of tiltrotor control variables. The control vector of the tiltrotor C= [τ,δw,δh,δv,Θ] ⊤ (4.3) Tiltrotor modelling for simulation... 895 consists of (respectively) tilt angle of nacelles, angles of deflection of ailerons, elevators and rudders, pitch control of rotors, which is 16 control states. The angular velocity of the rotor is assumed constant. 5. Calculation of trim in steady flight Including inertia loads into the general form of equations of motion, the tilt- rotor equations of motion have the form IpẎ p+ΩpIpY p+ 2 ∑ i=1 (IniẎ ni+ΩniIniY ni)+ 2 ∑ i=1 (IriẎ ri+ΩriIriY ri)= (5.1) =Qpg+ 2 ∑ i=1 Qnig + 2 ∑ i=1 Qrig +Qpa+ 2 ∑ i=1 Qnia+ 2 ∑ i=1 Qria In a steady flight, the accelerations are zero and the equations of motion are reduced to a system of 6 algebraic equations Q=−ΩpIpY p− 2 ∑ i=1 (IniẎ ni+ΩniIniY ni)− 2 ∑ i=1 (IriẎ ri+ΩriIriY ri)+ (5.2) +Qpg+ 2 ∑ i=1 Qnig+ 2 ∑ i=1 Qrig +Qpa+ 2 ∑ i=1 Qnia+ 2 ∑ i=1 Qria =0 For trim conditions the equations of motion have the general form Q(V p,Ωp,Φ,τ,δw,δh,δv,Θ)=0 (5.3) The values of trim controls in a steady flight are calculated using the Levenberg-Marquardt method to minimise total loads (5.3) acting on the til- trotor. This approach (Miller, 2004) allows one to obtain the trim states for the cases when the number of calculated trim parameters is arbitrary, i.e. less, equal or greater than the number of equations of motion. Minimizing thenonlinear functions numerically leads to compution of local minima, which may not be the real solution of the trim from the physical point of view. To avoid such cases, the total loads acting on the tiltrotor in a steady flight were monitored. Simulations of the tiltrotor flight were done using controls calculated in the trim procedure to prove that the parameters obtained by the minimisation method were correct. 896 M. Miller, J. Narkiewicz 5.1. Data for simulation The aim of simulation performed in this studywas to check the validity of themodel andtoget insight into tiltrotor behaviour in the trimmedstate.Data of V-22 were used in the simulation. A part of the design datawas taken from accessible literature (e.g., Miller and Narkiewicz, 2003; [11]). For parameters with no available data, the values were assumed as for corresponding parts of a similar aircraft (Miller, 2004). Base dimensions of the simulated tiltrotor are shown in Fig.6 (Miller and Narkiewicz, 2003) and given in Table 1. Fig. 6. Dimensions of V-22 Table 1.Tiltrotor data Rotor System Number of blades 3 Blades tip speed m/s (fps) 201.75 (661.90) Diameter m (ft) 11.58 (38.00) Disc area m2 (ft2) 210.70 (2268.00) Weights Take off kg (lbs) 15032 (33140) Dimensions Length, fuselage m (ft) 1748 (57.33) Width, rotors turning m (ft) 25.55 (83.33) Width, horizontal stabilizer m (ft) 5.61 (18.42) Height, nacelles fully vertical m (ft) 6.63 (21.76) Height, vertical stabilizer m (ft) 5.38 (17.65) Tiltrotor modelling for simulation... 897 5.2. Results of simulations Before simulating motion of the complete tiltrotor, separate modules and complete codes were debugged. Next, for the tiltrotor data, the results of sim- plified cases were checked for consistency with the proper reactions on the input. For instance, analysis of the influence of each control surfaces (Miller, 2004; Miller and Narkiewicz, 2001) on the tiltrotor flight was done for the selected flight phases (Swertfager and Martin, 1992). The tiltrotor flight si- mulations were done for three tiltrotor flight modes: helicopter, airplane and conversion. The longitudinal, symmetrical cases of the flight are presented in this paper. The side velocity V , angular velocities P, Q, R roll Φ and yaw Ψ angles as well as the rudder deflection angles δv = δv1 = δv2 were assu- med zero. The same values were assumed for the deflection of wing flaps δw = δw11 = δw12 = δw21 = δw22, angles of elevators δh = δh1 = δh2 = δh3 and collective control of the rotor pitch Θ0 = Θ01 = Θ02 for flaps on the wings and elevator. The longitudinal cyclic control of the rotor was assumed symme- trical and Θ1 = Θ11 = Θ12, whereas the lateral cyclic control unsymmetrical Θ2 = Θ21 =−Θ. During the trim calculations for the assumed values of forward U and vertical W flight velocities, the tiltrotor pitch angle Θ, deflection of flaps on thewing δw and elevators δh, collective Θ0 and cyclic Θ1 pitch of rotor blades and the nacelle tilt angle τ were computed. The airplane, conversion and helicopter modes of the tiltrotor flight were considered. The results of a steady flight in the airplanemode (horizontal with vertical climb) are given in Fig.7-Fig.10. The forward flight velocity was changed within the range U =10-180m/s, and thevertical velocitywas changedwithin the range W =−20-20m/s. The tiltrotor pitch angle obtained from simulations was equal 0 for the assumed flight conditions, and it was not presented in graphs. The calculated deflection angles of wing flaps (Fig.7) and elevators (Fig.8) did not exceed the values of available control surface deflections of V-22. Fig. 7. Deflection of flaps 898 M. Miller, J. Narkiewicz Fig. 8. Deflection of elevators Thecollective control of the rotor (Fig.9) is almost proportional to forward speed.When the forward speed is low (below 60m/s), theminimal tilt angle of nacelles (Fig.10) assuring proper values of the lift is about 70◦. Thedeflection of wing flaps is maximum for low forward speed, to balance the inclination of lift from the rotors. The deflection of elevators ismaximum,when the aircraft is in the conversionmode because of the necessity to provide a proper tiltrotor pitch moment. It has the maximum value when the conversion stops and the rotor axis is in the horizontal position. Fig. 9. Collective pitch Fig. 10. Tilt angle of nacelles 5.3. Steady forward flight with vertical climb in conversion mode In the conversionmode, a steady forwardflightwith vertical climb velocity with possibile deflection angle of the nacelle was simulated. The forward flight velocity was changed within the range U =20-180m/s. Tiltrotor modelling for simulation... 899 As in the previous case, the values of calculated deflection angles of wing flaps (Fig.11) and elevators (Fig.12) do not exceed the values of available control for V-22. The collective pitch of rotor is approximately proportional to forward speed. The deflection of wing flaps (Fig.11) is maximum for low forward speed. When the forward speed is low, the minimum tilt angle of nacelles (Fig.14) is about 55◦ assuring theproper aerodynamic lift. It becomes smaller (about 43◦) when the forward speed increases. Fig. 11. Deflection of flaps Fig. 12. Deflection of elevators Fig. 13. Collective pitch Fig. 14. Tilt angle of nacelles 900 M. Miller, J. Narkiewicz 5.4. Steady forward flight with vertical climb in helicopter mode In the helicoptermode, a steady forward flightwith vertical climbwith the tilt angle of the nacelles 90◦ was simulated. The forward flight velocity was changed within the range U =−70-70m/s. The vertical velocity was changed within the range W =−20-20m/s. Fig. 15. Collective pitch Fig. 16. Cyclic control (longitudinal) Fig. 17. Deflection of flaps For the assumed steadyflight conditions, the negligible small value of pitch angle of the tiltrotor is obtained. The collective control of rotor swash-plates (Fig.15) increases with the increase of vertical speed. The cyclic control of rotor pitch (Fig.16) varies in the opposite way: when the vertical speed in- creases, the pitch angle of tiltrotor is also stabilized. The inclination angles Tiltrotor modelling for simulation... 901 Fig. 18. Deflection of elevators of wing flaps (Fig.17) and elevators (Fig.18) obtained in calculations do not exceed available values for V-22. In the range of negative forward speed, the influence of flap deflection on tiltrotor motion is not substantial, but it beco- mes noticeable at positive forward speed greater than 10-20m/s. The sign of deflection of control surfaces depends on the direction of vertical speed. 5.5. Side flight in helicopter mode In the side flight in the helicopter mode, the forward flight velocity was changed within the range U = −70-70m/s, and the side flight velocity was changed within the range V =−20-20m/s. The tiltrotor vertical velocity the W =0. The lateral cyclic pitch control was symmetrical Θ1 = Θ11 = Θ12 and longitudinal cyclic control – unsymmetrical Θ2 = Θ21 = −Θ22. The aircraft pitch and yaw angles, pitch and yaw rates and both rudder deflections were assumed zero. Fromequilibriumconditions, the roll angle of tiltrotor Φ, inclination angles of wing flaps δw = δw11 = δw12 = δw21 = δw22, inclination angles of elevators δh = δh1 = δh2 = δh3 and collective control of the rotor swash-plate Θ0 = Θ01 = Θ02 were calculated. The results are shown in Fig.19-Fig.24. In these flight conditions, the tiltrotor roll angle depends on side velocity. The deflections of flaps and elevators are small. The value of the cyclic pitch control depends on the value of side speed and its sign. Fig. 19. Deflection of flaps 902 M. Miller, J. Narkiewicz Fig. 20. Deflection of elevators Fig. 21. Collective pitch Fig. 22. Roll angle of tiltrotor Fig. 23. Cyclic control (lateral) Fig. 24. Cyclic control (longitudinal) Tiltrotor modelling for simulation... 903 5.6. Stability and controllability For stability analysis, of the tiltrotor the equations of motion are nume- rically linearized with respect to state and control variables for given trim values. This allows analysis of the control matrix and calculation of eigenva- lues and eigenvectors for examination stability. This option of model analysis allows one to investigate the stability and controllability of the aircraft. For the assumed tiltrotor design data, both stability and controllability analysis was made for the whole forward speed range of the tiltrotor. In the range of velocity considered, complex eigenvalueswithnegative real partswere obtained only for low velocity of forward flight (from 60 to 70m/s), when the conversion mode occurs. In other forward speed, negative real parts of eigenvalues were obtained. These results show that in the steady flight, the aircraft is stable in the whole range of flight velocities. Analysing the control matrix, it may be stated that the control variables influence the related loads with minor cross-coupling effects. For the airplane mode, this is summarized in Table 2. Table 2.Controllability analysis in the airplane mode Loads Control elements airplane mode transition helicopter mode Fx rotor collective pitch nacelles tilt angle, rotor collective pitch rotor longitudinal cyclic pitch Fy rudders (no lateral cyclic pitch) rudders (no lateral cyclic pitch) rotor lateral cyclic pitch Fz flaps, elevators (not nacelles tilt nor longitudinal cyclic pitch used) flaps, nacelles tilt, rotor collective pitch, elevators (smaller influence) rotor collective pitch Mx flaps, elevators, rotor collective pitch deflection of flaps, rotor collective control asymmetric rotor collective pitch My elevators (rotor longitudinal cyclic pitch not used) deflection of elevators, rotor collective control rotor longitudinal cyclic pitch (symmetric) Mz rotor collective pitch, rudders rotor collective pitch, rudders rotor collective pitch, rotor longitudinal cyclic pitch 904 M. Miller, J. Narkiewicz 5.7. Results of numerical simulations For the calculated steady flight parameters, simulations of flight were car- ried out to check their accuracy. The steady flight in the airplane mode was simulated for: U = 180m/s, V = W = P = Q = R = 0, Φ = Θ = Ψ = 0 and the trim parameters obtained from calculations were: deflection of wing flaps δw = −6.75 ◦, elevators δh = 7.23 ◦ and collective pitch of the rotor Θ0 =45.91 ◦. The tiltrotor displacements are presented in Fig.25 (horizontal) and in Fig.26 (vertical). Thevariations of attitude angles are very small (below 0.1◦). Fig. 25. Horizontal displacement in the airplanemode Fig. 26. Vertical displacement in the airplanemode It can be seen that during 1800mdistance, the flight altitude decreased by about 0.1m and side translation was about 0.1m to the left. These values are attributed to numerical errors in calculation of the trim values. Similar results were obtained in the transition and helicopter modes. 6. Conclusions The tiltrotor computer model was developed for flight simulation, trim stabi- lity and control analysis. The model is composed of rigid elements: fuselage, wings, empennage, rotors, but due to the modularity of the code these as- sumptions may be easily released. The design parameters of V-22 tiltrotor were used for simulations. Some data of the tiltrotor had to be assumed, and there was no possibility to compare the results of numerical simulations with the flight data. On the grounds of results of calculations performed, it may be concluded that the developedmodel of the tiltrotor works properly. Tiltrotor modelling for simulation... 905 A tiltrotor is a complex rotorcraft, and several simplifying assumptions had to be applied. Theymight be released by adjusting themodel for specific needs of a particular helicopter. Acknowledgement The researchwas carried out as a part of a project supported by the State Com- mittee for ScientificResearch”The research on the control of tiltrotor in selected flight stages”, grant No. 5T12C06724. References 1. CicaleM., 2003,ERICA: theEuropeantiltrotordesignandcritical technology projects, 29th European Rotorcraft Forum, Friedrichshafen, Germany 2. McVeigh M.A., Nagib H., Wood T., Kiedaisch J., Stalker A., Wy- gnanski I., 2004,Model and full scale tiltrotor download reduction tests using active flow control, 60th Annual Forum, 1, Baltimore,MD 3. McVicar J.S.G., Bradley R., 1992, A generic tilt-rotor simulation model with parallel implementation and partial periodic trim algorithm, 18th Euro- pean Rotorcraft Forum, Avignon, France 4. Miller M., 2004, The Control of Tiltrotor in Chosen Flight Stages, Ph.D. Thesis,Warsaw, Poland 5. Miller M., Narkiewicz J., 2001, Control of tiltrotor aircraft, 4th National Rotorcraft Forum,Warsaw, Poland 6. Miller M., Narkiewicz J., 2003, Simulation of tiltrotormotion, 29th Euro- pean Rotorcraft Forum, Friedrichshafen, Germany 7. Nixon M.W. et al., 2003, Aeroelastic stability of a four-bladed semi- articulated soft-in-plane tiltrotor model, 59 AHS Forum, 110, Phoenix, AZ, USA 8. PolakD.R.,GeorgeA.R., 1998,Flowfield and acousticmeasurements from amodel tiltrotor in hover, Journal of Aircraft, 35, 6, 921-929 9. SrinivasV., Chopra I., 1996,Validation of a comprehensive aeroelastic ana- lysis for tiltrotor aircraft,AHS 52nd Annual Forum,Washington, D.C., USA 10. Swertfager T.A., Martin S. Jr., 1992, The V-22 for SOF, 48th American Helicopter Society Annual Forum, Washington, D.C. 11. V-22 Osprey Technical Specifications, http://www.boeing.com/rotorcraft/ military/v22/v22spec.htm 12. Weakley J.M., Kleinhesselink K.M., Mason D.H., Mitchell D.G., 2003, Simulation evaluation of V-22 degraded-mode flying qualities, 59 AHS Forum, 135, Phoenix, AZ, USA 906 M. Miller, J. Narkiewicz Modelowanie tiltrotora dla symulacji w różnych stanach lotu Streszczenie Opracowano symulacyjnymodel statku powietrznego typu tiltrotor przeznaczony do symulacji lotu oraz analizy osiągów, stabilności i sterowania.Model wiropłata zło- żony jest z kadłuba, skrzydeł, usterzenia ogonowego, gondoli silnikowych i wirników. Równania ruchu zostały uzyskaneprzez sumowanie obciążeń od sił bezwładności, gra- witacyjnych i aerodynamicznychdziałających na każdy element statku powietrznego. Obciążenia aerodynamiczne skrzydeł, stateczników i łopat wirników zostały obliczo- ne z zastosowaniem quasistacjonarnego modelu opływu. Do wyznaczania prędkości indukowanej wirników zastosowano model Glauerta. Wpływ strumienia zawirniko- wego na skrzydła i stateczniki jest obliczany z wykorzystaniem aktualnej wartości prędkości indukowanej wirników. Program do modelowania wiropłata został opraco- wany w środowiskuMatLab. Program zbudowany jest z modułów obliczeń obciążeń poszczególnych elementówwiropłata, które wykorzystywane są również dowyznacza- nia warunków lotu ustalonego, stateczności i sterowności. Podczas pierwszego etapu badań wyznaczono warunki ustalonego lotu tiltrotora w różnych konfiguracjach, co pozwoliło zbadać zachowanie i potwierdzić poprawnośćmodelu. Manuscript received February 9, 2006; accepted for print July 5, 2006