Mathematical Modelling of Gimballed Tilt-Rotors for Real-Time Flight Simulation

Abstract

This paper introduces a novel gimballed rotor mathematical model for real-time flight simulation of tilt-rotor aircraft and other vertical take-off and landing (VTOL) concepts, which improves the previous version of a multi-purpose rotor mathematical model developed by ZHAW and Politecnico di Torino as part of a comprehensive flight simulation model of a tilt-rotor aircraft currently implemented in the Research and Didactics Simulator of ZHAW and used for research activities such as handling qualities studies and flight control systems development. In the novel model, a new formulation of the flapping dynamics is indroduced to account for the gimballed rotor and better suit current tilt-rotor designs (XV-15, V-22, AW-609). This paper describes the mathematical model and provides a generic formulation as well as a specific one for 3-blades proprotors. The method expresses the gimbal attitude but also considers the variation of each blade’s flapping due to the elasticity of the blades, so that the rotor coning angle can be represented. A validation of the mathematical model is performed against the available literature on the XV-15 Tilt-rotor aircraft and a comparison between the previous model is provided to show the improvements achieved. The results show a good correlation between the model and the reference data and the registered performance allow real-time flight simulation with pilot and hardware in the loop.

1. Introduction 1.1. Motivation for Research

Current developments in vertical take-off and landing (VTOL) technology such as the advancement in structures, automation and control and energy generation-storage-utilization are driving aircraft designers towards a wide range of innovative architectures that may better suit the transport of people and goods in high population density scenarios. The technological challenges the designers have to face are several, involving propulsion efficiency, noise, performance, safety and airworthiness and so forth. The industry and research institutes have so far responded to these challenges with multiple VTOL architectures, partially driven by a rather wide background of experimental solutions which have already been investigated in the past decades, as duly documented in References [1,2], yet enforced by renewed enthusiasm and business interest as well as improved design tools. It is fair to recognize that, in spite the wave of enthusiasm and commitment that invested the VTOL world between the 1960s and the late 1990s, most of the VTOL concepts were proven quite unsuccessful or did not make it to production. Among the many, however, thanks to the XV-15 Research Aircraft [3], tilt-rotors did make it through the research phases and they found application in the military world with the Bell Boeing V-22 Osprey and will soon reach civil certification with the Leonardo AW-609 (formerly BA-609). A tilt-rotor is an aircraft that can convert its rotors within a range exceeding 90 degrees, from helicopter mode (which allows vertical take-off) to airplane mode and therefore providing enhanced vehicle performance in terms of speed, endurance and fuel consumption. The tilt-rotor concept offers high potential for commercial aircraft, military platforms as well as unmanned aerial vehicle (UAV)/drone applications. The research interest in tilt-rotors is renewed today with the uprising of air taxi concepts combined with the need for more efficient and cost effective VTOL platforms and the innovation effort in the more electric and green propulsion field.

1.2. Tilt-Rotor Real-Time Simulation

Reliable flight simulation models are not easy to obtain, they require a considerable effort in terms of aerodynamic data processing and mathematical modelling and both activities involve a great deal of time, expertise and economical resources. Nevertheless, flight simulation models play a crucial role in flight control law design as well as handling qualities prediction of highly complex systems such as convertiplanes. Although many publications have been produced so far about tilt-rotors, especially related to the XV-15 and the V-22 aircraft, sources describing mathematical models for specific real-time flight simulation applications are limited. First substantial attempts of developing a generic tilt-rotor model for real-time flight simulation can be traced back to the XV-15 Tilt-Rotor Research Aircraft, within which the Generic Tilt Rotor Simulator (GTRS) was developed. In spite of some additional information regarding a simulation model of the experimental Boeing Vertol Model reported in 1973 in Reference [4], the Bell XV-15 (early Model 301) mathematical model is generally used as the baseline of current research and development activities. Its first development was described also in 1973 in Reference [5], while in Reference [6] Ferguson reports the last issue of the GTRS model available on public domain. The GTRS, however, is believed to have been further modified and extensively improved and it is currently being employed for flight simulation of the V-22 Osprey aircraft of both US Army and US Navy. If References [5,6] are compared, an important modification of the aerodynamic mathematical model can be noticed, mainly implemented because, as Ferguson reports, some model predictions did not match successfully the validation data and were found to be overly optimistic, especially in hover condition.

1.3. Gimballed Proprotors

Rotor Mathematical models of different levels of detail can be developed and the degree of complexity relates to the specific application. The specific case of tilt-rotor, however, requires some clarification: tilt-rotor aircraft generally makes use of gimballed proprotors (the name proprotor essentially comes from the combined propeller-rotor function this system provides in its operational range between helicopter and airplane mode). As reported in Reference [7], a gimballed rotor is basically a hingeless rotor attached to the rotor shaft through a gimbal, so the rotor head is able to tilt relative to the shaft with a mild elastic restraint and the hub tilt degree of freedom minimizes the blade flapping relative to the hub. As reported by Padfield in Reference [8] gimballed rotors generally use either a Universal joint or a constant velocity (homokinetic) joint. According to Reference [2], the XV-15 gimballed rotor exploits the Universal Joint principle (see Figure 1 and Figure 2).

In both cases, gimballed rotors present a substantial difference with respect to conventional articulated blade rotor systems, since if the elastic property of the blades is neglected, in gimballed rotors blade centrifugal forces are carried through the central hub structure and are balanced by the other blades, yielding to the suppression of the single blade flapping motion. However, the gimballed rotor does present a certain stiffness: the inner stiffness of the gimballed hub does provide a spring restraint and allows the torque transfer between the shaft and the blades but as reported by Carlson in Reference [7], the nature of the resulting fundamental stress and vibration is supposed to be reduced and advantageous with respect to conventional hinge-less or bearing-less rotors. Furthermore, the blades of a gimballed rotor are in truth not infinitely rigid, so the blades will also be able to slightly flap independently despite being constrained to the gimballed hub. Consequently, if the responses of an articulated and a gimballed 3 blades rotor to a cyclic swashplate input are compared, the reader shall expect the response of the gimballed rotor to show reduced yet not absolutely absent 1/rev blade flapping motion due to the structural flexibility of the blade constrained to the hub. Furthermore, the hub loads are expected to present a 3/rev oscillation, being each blade’s airloads a function of the azimuth position and being the overall hub loads the result of a summation of the loads of the 3 blades. Thus, in gimballed rotors oscillatory loads are not completely cancelled and several are the studies focusing on the development of control systems to mitigate the effect of those oscillation within the overall aeroelastic behaviour of the aircraft, among which are References [9,10]. The mathematical modelling of gimballed rotors has been addressed in the past and several fruitful models were produced to properly address the aeromechanic and aeroelastic design of such rotors, the author can refer to comprehensive mathematical models such as Johnson’s CAMRAD [11] as well as specific tilt-rotor studies such as References [7,12], yet the complexity of these models has always been limiting their performance and its applicability to real-time flight simulation. As Heffley duly discussed in 1988 in Reference [13] for the conventional helicopter case, if a minimum required level of complexity is to be defined, mathematical models for piloted flight simulation shall include an adequate representation of the main rotor dynamics to properly reproduce the hover as well as the forward flight characteristic modes of the rotorcraft. Among all the complex phenomena involved in the rotor dynamics, Heffley reckons that the coupling between the rigid-body modes and the rotor flapping does affect the handling qualities of the rotorcraft and that the flapping is most important to a pilot-in-the-loop simulation due to the characteristic control lag following cyclic input. Feathering and lead-lag dynamics are instead not mentioned in Reference [13] among the minimum set of dynamics to be included in a mathematical model for piloted simulations. However, while a Lead-Lag and Feathering dynamic formulations might be worth implementing in mathematical models for piloted simulations of conventional articulated rotors (an extensive model was developed by Takahashi in Reference [14] and could be revised for real-time applications), these additional dynamics are usually neglected for the tilt-rotor case, being proprotors generally considered stiff-in plane. As reported in the previous section of this paper, the baseline reference for tilt-rotor mathematical modelling for flight simulation used a simplified and static model of the proprotors flapping dynamics that made it suitable for piloted simulations with the software and hardware available at that time: while Reference [5] implemented a linearised static model of the proprotor disk flapping attitudes (lateral and longitudinal cyclic attitudes and conicity), Reference [6] introduces a first order dynamic model of the overall proprotors Tip Path Plane (TPP) flapping motion:

_{C11C12C21C22a¯˙1b¯˙1}+_{A11A12A21A22a¯1b¯1}=_B1B2,

where_a¯1and_b¯1are the longitudinal and lateral attitudes of the TPP, and decouples the coning dynamics reducing it to a static relation:

_a0=0.75·T/_nb+_Kβ·_a¯0Ib·^Ω¯2+_Kβ,

where T is the rotor thrust perpendicular to the disk plane,_Ibis the blade inertia,_nbthe number of blades,_a¯0, the precone angle andΩ¯ the non-dimensional rotor angular velocity. In Equation (1), the direct terms_A11and_A22of the stiffness matrix are influenced by the blades collective pitch, the rotor advance ratio and spring constant_kG, which has a different value compared to_kβ in Equation (2), to account for the different spring characteristic of the gimbal and that of blade structure. The numerical values of_kGand_kβ are reported in Table 1.

However, according to what reported by Ferguson in Reference [6], since the required solution time would have exceeded that allowable for the real-time simulation using the Sigma 8 computer located at NASA ARC, the dynamic flapping model was not implemented in the simulation and the blade flapping due to cyclic inputs is assumed to occur instantaneously considering the rotor in an equilibrium condition, yielding to Equation (3). As reported by Ferguson in Reference [6], the requirement for that model was to “keep the computational loop time to less 70 ms in order to maintain real-time” (p. 7), and the “differential equations for blade flapping that would properly account for the rotor following time were determined to require a solution time in excess of that allowable” (p. 12). The authors of the present paper assume that the limits of the model in Reference [6] were related to the hardware and software set up available back then in 1980, the use of the same model with updated hardware would in principle allow the use of the dynamic flapping model shown in Equation (1). However, the authors of the present paper were not able to fully reproduce Ferguson’s model only based on Reference [6], due to the lack of a complete description of the mathematical derivation and the unclear definition of some variables, as reported in Reference [15].

_{A11A12A21A22a¯1b¯1}=_B1B2.

Relatively more recent works generally reduce the flapping dynamics to that of the TPP as already introduced in Reference [6], modelled as a first-order system, or, like McVicair in Reference [16], implement a detailed modelling of the blades’ flapping motion which allows to include the 3/rev oscillations of the hub loads but treats the proprotor as articulated rather than gimballed. On the other hand, Padfield suggests in Reference [8] the FLIGHTLAB gimbal formulation, which neglects the first harmonic flapping and the coning due to the blades’ flexibility and only accounts for the overall gimbal motion.

In previous developments of their tilt-rotor model for real-time simulation, the authors of the present paper used the GTRS described in Reference [6] as baseline for the overall aerodynamic modelling of the generic tilt-rotor but replaced the reference rotor model with a novel one which uses a multi-blade formulation of the rotor’s airloads and flapping motion and tested it successfully in the Research and Didactics Simulator (ReDSim) of ZHAW, a reconfigurable flight simulation platform for both helicopters and airplanes. This rotor model exploits a second order semi-analytical formulation of the flapping motion of each blade of the rotor and it is reported here in Equation (4) in short form for the j-th blade, being m the inertial term, c the damping, k the stiffness of the flapping hinge and_Fβthe external forcing input.

m_βj¨+c_βj˙+k_βj=_Fβ,j.

The terms in Equation (4) can be expressed as

m=_Iβc=−_M1k=−_M2+_{mb rcg az2}+_KβFβ=_M3−_{mb rcg}(_az3−_gz),

where_M1,_M2,_M3represent the aerodynamic moment resulting from the integration along the blade span,_Iβis the blade moment of inertia,_mbis the blade mass,_rcgis the distance of the centre of gravity of the blade from the origin of the hub reference system,_gzis the acceleration of gravity and_az,2,_az,3 z-components of the acceleration due to the blade motion in space (see next section for more details). Equation (4) is computed for each blade and the result fed back into the inflow and airloads computation. Forces and moments are computed for each blade and summed up at the origin of the hub reference system. As reported in Reference [15], tests performed on the multi-purpose rotor mathematical model showed good correlation with the experimental data in terms of overall hub forces and moments for both the XV-15 tilt-rotor and the UH-60 helicopter cases. However, the model in Reference [15] was essentially formulated as a fully-articulated rotor in which a spring component_Kβ (see Equation (5)) in each blade’s flapping equation was introduced and tuned according to what reported in Reference [6] to better approximate the behaviour of a gimballed rotor. Therefore, this formulation in Reference [15] could not distinguish the stiffness of the gimbal, which influences the cyclic response of the proprotor, from the stiffness of each blade, which causes the rotor cone. As reported in Table 1, the stiffness of the gimbal is generally very small and tuned to provide adequate cyclic control of the rotor, therefore it is not strictly comparable to that of the metal blade-retention. The novel model presented in the following sections of this paper takes inspiration from the elegant and simple formulation of 2-bladed rotor flapping dynamics introduced by Chen in Reference [17] as well as the insightful descriptions of the gimbal rotors in Padfield’s work [8] and aims to modify the multi-articulated rotor formulation previously developed by the authors of Reference [15] to distinguish the contribution of the gimbal stiffness_KGand that of_Kβassociated to the elastic flapping of the blade and responsible for the coning of the rotor disk. While the extended formulation is described in the following sections of this paper, the resulting flapping equation can be written in rotating axis, for the j-th blade, as

m_βj¨+c_βj˙+_cstructΔ_βj˙+k_βj+_{kG βG,j}+_kβΔ_βj=_Fβ,j,

where m, c and_Fβ are unchanged with respect to Equation (5), while_cstructis an additional term accunting for the structural damping._βGrepresents the gimbal attitude in the rotating reference system and can be related to the overallβby the matrix_Bψ, a linear combination of the blades’ azimuthal distance

_βG=_Bψ·β.

βnow represents the overall angle between the blade and the plane perpendicular to the rotor mast (origin of the non-rotating reference system), while the actual blade flapping due to the blade elasticity and responsible of the disk coning is

Δβ=β−_βG.

This formulation is expected to better correlate with the dynamic response of the gimballed rotors without neglecting the elastic blade flapping motion and the resulting coning dynamics. Furthermore, the novel formulation allows the improvement of the response without overly increasing the complexity of the mathematical model and ensuring the real-time performance of the resulting Simulink code without making use of an expensive and less flexible hardware/software set up. 2. Mathematical Model 2.1. Reference Frames

The presence of a gimbal joint introduces an additional degree of freedom, so the set of reference systems used in Reference [15] and reported here in Figure 3 to describe the motion of an articulated rotor has been modified by adding a reference system which, according to the assumptions made later in Section 2.5, is aligned to the rotor tip-path-plane. The four reference frames are listed below.

Hub —Fixed, centered in the rotor hub, with:

• _xhin the hub plane and positive forward;

• _yhin the hub plane and positive outboard (starboard for the right rotor, port for the left rotor);

• _zhperpendicular to the hub plane and positive downwards.

Rotating

—Centered in the rotor hub, rotating with angular speed equal toΩ, and having:

• _xrcoincident with the projection of the blade axis on the hub plane and positive outwards;

• _zrperpendicular to the hub plane and positive upwards;

• _yrperpendicular to_xrand_zr(hence positive in the blade advancing direction).

Tip-path-plane (TPP)

—Centered in the rotor hub, rotating about_yrwith angular speed equal to−_β˙G; its axes are:

• _ytppcoincident with_yr(perpendicular to the blade and positive in its advancing direction);

• _ztppperpendicular to the rotor tip-path-plane and positive upwards.

• _xtppperpendicular to_ytppand_ztpp;

Blade

—Centered in the blade flapping hinge, rotating with angular speed equal to−Δβ˙; its axes are:

• _xbcoincident with the blade axis and positive from the flapping hinge to the blade tip;

• _ybcoincident with_ytppand_yr(perpendicular to the blade and positive in its advancing direction);

• _zbperpendicular to_xband_yb.

Transformations between different reference frames are obtained using the following rotation matrices:

• From Blade system to TPP system:

  _TΔβ=cosΔβ0−sinΔβ010sinΔβ0cosΔβ=cos(β−_βG)0−sin(β−_βG)010sin(β−_βG)0cos(β−_βG),

  where bothβand_βGare positive when the blade flaps upwards.

• From TPP system to Rotating system:

  _TβG=cos_βG0−sin_βG010sin_βG0cos_βG

• From Rotating system to Hub system:

  _Tψ=−cosψsinψ0sinψcosψ000−1,

  whereψ=Ωtis the blade azimuth angle, that is, the angle between_xrand the negative direction of_xh, and is positive according to the direction of the rotor angular speedΩ(counter-clockwise for the right rotor, clockwise for the left rotor).

2.2. Blade Position, Speed and Acceleration

Position, speed and acceleration of a generic blade section S are determined with respect to the frames of reference defined in Section 2.5. The following vectors are defined, representing, respectively, the aircraft linear speed, linear acceleration, rotational speed and rotational acceleration, with respect to the Hub reference frame:

_Vh=_uhvhwh;_ah=_u˙hv˙hw˙h;_ωh=_phqhrh;_ω˙h=_p˙hq˙hr˙h.

The position vector of the blade section in the Blade reference frame is

_rS,b=r00,

while the position of S in the TPP reference frame:

_rS,tpp=_TΔβ·_rS,b+e00,

where e is the flapping hinge offset. The speed of the blade section in the TPP reference frame is given by the following three contributions.

• Relative speed of S with respect to the TPP system, due to blade rotation around the flapping hinge:

  _VS,rel,tpp=_TΔβ·0−Δβ˙0×_rS,b.

• Linear speed of the origin of the TPP system (due to aircraft motion):

  _Vtpp=_TβG−1·_Tψ−1·_Vh.

• Tangential speed owing to reference frame rotation, which in turn is due to aircraft motion, gimbal motion and rotor spinning:

  _VS,tang,tpp=_ωtpp×_rS,tpp,

  where:

  _ωtpp=_TβG−1·_Tψ−1·_ωh+0−_β˙G0+00Ω.

Therefore, the total section speed in the TPP reference frame is expresse as

_VS,tpp=_VS,rel,tpp+_Vtpp+_VS,tang,tpp,

while the velocity of the blade section in the blade system results in

_VS,b=_TΔβ−1·_VS,tpp.

The blade section acceleration in the TPP reference frame is given by the sum of the following five contributions:

• relative acceleration of S (due to blade flapping) with respect to the TPP system, which in turn comprises a tangential and a centrifugal component:

  _aS,rel,tpp=_TΔβ·0−Δβ¨0×_rS,b+0−Δβ˙0×0−Δβ˙0×_rS,b,

• linear acceleration of the origin of the TPP system (due to aircraft motion):

  _atpp=_TβG−1·_Tψ−1·_ah.

• tangential acceleration owing to angular acceleration of the TPP reference frame:

  _aS,tang,tpp=_ω˙tpp×_rS,tpp,

  where:

  _ω˙tpp=_TβG−1·_Tψ−1·_ω˙h+0−_β¨G0+00Ω˙

• centrifugal acceleration, due to rotation of the TPP reference frame:

  _aS,centr,tpp=_ωtpp×_ωtpp×_rS,tpp.

• Coriolis acceleration:

  _aS,Cor,tpp=2_ωtpp×_VS,rel,tpp.

The resulting total blade section acceleration in the TPP system is:

_aS,tpp=_aS,rel,tpp+_atpp+_aS,tang,tpp+_aS,centr,tpp+_aS,Cor,tpp.

Finally, the total acceleration in the Blade reference frame can be obtained as:

_aS,b=_TΔβ−1·_aS,tpp.

2.3. Rotor Aerodynamics

The blades airloads are modelled according to the Blade Element Theory as suggested in Reference [18]. Each rotor blade is divided into a finite number of sections, the aerodynamic loads being determined for each section and then integrated along the blade span. The aerodynamic forces generated by a generic blade section are represented in Figure 4.

Only lift and drag are considered in the present model, while the pitching moment about the aerodynamic centre is neglected. In order to compute the aerodynamic forces, the (absolute) angle of attack of each blade section must be determined:

α=_ϑtwist+_ϑ0−_A1cosψ−_B1sinψ−_K1β−_α0+ϕ,

where:

• _ϑtwistis the twist angle, function of the position along_xb;

• _ϑ0,_A1and_B1are, respectively, the collective pitch, the lateral and the longitudinal cyclic;

• _K1=tan_δ3, being_δ3the pitch-flap coupling angle;

• _α0is the zero-lift angle of attack;

• ϕis the inflow angle.

The inflow angle identifies the direction of the airflow affecting the blade section, and can in turn be expressed as:

ϕ=arctan_{V⊥ V‖}.

_V⊥is the component of the airflow speed, in wind axes, which is perpendicular to the blade section; to determine it, the following effects are taken into account:

• blade motion_VS,b;

• rotor inflow speed

  _Vinflow=ΩR_λ0+rR_λ1ccosψ−rR_λ1ssinψ,

  in which the three inflow states_λ0,_λ1cand_λ1s can be computed according to Pitt-Peters dynamic Model as described in Reference [19].

_V‖is the component of the airflow speed, in wind axes, parallel to the blade section, for which only the blade motion_VS,b[2]is taken into account. Both_V⊥and_V‖are non-linear functions ofβ,_βG,β˙and_β˙G, thus making alsoϕandαnon-linear functions of such quantities. Assuming small angles, it is possible to write:

sinβ≈β;sin_βG≈_βG;andcosβ≈1;cos_βG≈1.

Consequently, neglecting higher-order infinitesimals and considering the hinge offset to be null (the full derivation is not reported here for the sake of brevity), the expressions of_V⊥and_V‖can be written as

_V⊥=_{B1 xuo}−r_λ1cΩ+_qhr−β_uhcosψ++−_{A1 xuo}−r_λ1sΩ+β_vh+_phrsinψ−ΩR_λ0−rβ˙+_whV‖=β_phr+_vhcosψ+−β_qhr+_uhsinψ−r−Ω+_rh.

As a result of the linearisation all terms in_βGand_β˙G in Equation (32) disappear and two linear functions of onlyβandβ˙ are obtained. Moreover, these functions are identical to those derived in Reference [15] for an articulated rotor having hinge offsete=0. Sinceϕis itself a non-linear function of_V⊥and_V‖, a linearisation is performed by means of a first-order Taylor expansion around the point(β=0,β˙=0):

ϕ(β,β˙)≈ϕ(β=0,β˙=0)+∂ϕ∂β˙_|β=0β˙=0β˙+∂ϕ∂β_|β=0β˙=0β

Thus, the inflow angle becomes a linear function having the following form:

ϕ=_ϕ1β˙+_ϕ2β+_ϕ3

where:

_ϕ1=_{n1 d3n32}+_d32ϕ2=_{n2 d3}−_{n3 d2n32}+_d32ϕ3=arctan_{n3 d3}

and

_n1=−r_n2=−_uhcosψ+_vhsinψ_n3=(_{B1 xuo}−r_λ1c)Ω+r_qhcosψ+(−_{A1 xuo}−r_λ1s)Ω+r_phsinψ−ΩR_λ0+_whd2=rcosψ_ph−sinψ_qhd3=cosψ_vh+sinψ_uh−r(_rh−Ω).

Equation (33) substituted into Equation (28), yields:

α=_α1β˙+_α2β+_α3

where:

_α1=_ϕ1α2=−_K1+_ϕ2α3=_ϑtwist+_ϑ0−_A1cosψ−_B1sinψ−_α0+_ϕ3.

2.4. Aerodynamic Forces and Moments

The blade section angle of attack being expressed according to Equation (28), if the airfoil lift and drag curves are known, lift and drag can be determined for each blade section. Since both_Cland_Cddepend onαand the latter is a function ofβandβ˙, their expressions are linearised using a first-order Taylor expansion in the neighbourhood of(β=0,β˙=0):

_Cl≈_{Cl |β=0β˙=0}+d_Cldα∂α∂β˙_|β=0β˙=0β˙+d_Cldα∂α∂β_|β=0β˙=0β=_{Cl |β=0β˙=0}+_{Clα |β=0β˙=0α1}β˙+_{Clα |β=0β˙=0α2}β

_Cd≈_{Cd |β=0β˙=0}+d_Cddα∂α∂β˙_|β=0β˙=0β˙+d_Cddα∂α∂β_|β=0β˙=0β=_{Cd |β=0β˙=0}+_{Cdα |β=0β˙=0α1}β˙+_{Cdα |β=0β˙=0α2}β.

Onceϕ,αand the aerodynamic coefficients are known, the aerodynamic forces can finally be computed. The lift and drag generated by each blade section are, respectively

dL=12ρ^V2c_Cldr;dD=12ρ^V2c_Cddr,

whereV=_V⊥2+_V‖2is the airflow speed. By projecting the lift and drag into the Blade reference frame, the following forces and moments are obtained:

• aerodynamic force along_zb:

  d_Fa=dDsinϕ+dLcosϕ=12ρ^V2c(_Cdsinϕ+_Clcosϕ)dr

• aerodynamic force along−_yb:

  d_Qa=dDcosϕ−dLsinϕ=12ρ^V2c(_Cdcosϕ−_Clsinϕ)dr,

• aerodynamic moment around−_yb:

  d_Ma=d_Far=12ρ^V2c(_Cdsinϕ+_Clcosϕ)rdr,

• aerodynamic moment−_zb:

  d_Na=d_Qar=12ρ^V2c(_Cdcosϕ−_Clsinϕ)rdr.

For implementation in the flapping dynamics equation, the terms inβandβ˙in the expression of_maare collected as

d_Ma=d_M1β˙+d_M2β+d_M3,

with:

d_M1=12ρ^V2c(_Cdα(_α3)sinϕ+_Clα(_α3)cosϕ)_α1rdrd_M2=12ρ^V2c(_Cdα(_α3)sinϕ+_Clα(_α3)cosϕ)_α2rdrd_M3=12ρ^V2c(_Cd(_α3)sinϕ+_Cl(_α3)cosϕ)rdr.

The overall aerodynamic loads generated by each blade (as shown in Figure 5) are subsequently obtained by integrating the above forces and moments along the whole length of the blade:

_Fa=_∫0Rd_Fa;_Qa=_∫0Rd_Qa;_Ma=_∫0Rd_Ma;_Na=_∫0Rd_Na

_M1=_∫0Rd_M1;_M2=_∫0Rd_M2;_M3=_∫0Rd_M3;

Integration along the blade is operated symbolically, using the method of the trapezoids. LetdFbe a generic aerodynamic load generated by a blade section, and n the number of intervals in which the blade has been divided; the length of each interval is designated asΔx. The overall aerodynamic load generated by the blade is then:

F=_∫0RdF≈Δx2dF(_x1)+2dF(_x2)+2dF(_x3)+⋯+2dF(_xn)+dF(_xn+1).

Airfoil Aerodynamic Coefficients

As the blades of a tilt-rotor aircraft proprotors operate in a wide range of angles of attack, the complete airfoil polar curves (i.e., for angles of attack spanning from−180to+180degrees) need to be included in the model inizialization file. Available data for the NACA 64-208 airfoil covered only a limited range of angles of attack; therefore, the missing parts of the polar curves have been reconstructed using the empirical approach described hereafter.

• The airfoil lift and drag coefficients are expressed by tables that are function of the angle of attack while for ranges of angle of attack outside available data the approximation proposed by Hoerner in Reference [20] is used:

  _Cl=_kClsinαcosα_Cd=_kCdsin^α2,

  where the coefficients_kCland_kCd are tuned to match the reference data in Reference [21].

• The missing portions of the polar curves are determined using Equation (46), and introducing corrections in order to best fit the available aerodynamic data.

• The derivatives of the aerodynamic coefficients with respect to the angle of attack (_Clαand_Cdα) are computed using the centred difference method:

  ^f′(_xj)≈f(_xj+1)−f(_xj−1)_xj+1−_xj−1.

2.5. Flapping Dynamics Model for Gimballed Rotors Instead of explicitly modelling the gimbal dynamics and increasing the complexity of the existing model, the kinematic constraint on the blades is taken into account by introducing the following two assumptions:

• The flapping angle of each blade is computed by adding to the gimbal attitude,_βG, another contributionΔβ, representing the variations of flapping angle in each of the N blades due to the blade flexibility or to a coning hinge:

  _βj=_βGj+_Δβjj=1,…,N.

• The gimbal plane is assumed to be parallel to the rotor tip-path-plane, hence:

  _βGj=_β1ssin_ψj+_β1ccos_ψjj=1,…,N,

• the angles_β1sand_β1c in Equation (49) are the lateral and longitudinal attitudes of the TPP and are obtained using the multi-blade coordinate transformation (MBC) truncated to the first order, therefore

  _β1s=2N∑j=1N_βjsin_ψj;_β1c=2N∑j=1N_βjcos_ψj.

Strictly speaking, the second hypothesis would be verified only ifΔβwas identically null. However, the approximation is reasonably acceptable when considering that the coning stiffness of the XV-15 rotor hub is in the order of¹⁰⁵ Nm/deg (Reference [6]). Evaluating_βGfor each blade, and substituting in its expression the definitions of the multi-blade coordinates, the following matrix is obtained:

_βG1βG2⋮_βGN=2N1cos(_ψ2−_ψ1)…cos(_ψN−_ψ1)cos(_ψ1−_ψ2)1…cos(_ψN−_ψ2)⋮⋮⋱⋮cos(_ψ1−_ψN)cos(_ψ2−_ψN)…1_β1β2⋮_βN,

which can be written in compact form as

_βG=_Bψ·β

and where

_Bψi,j=2Ncos(_ψj−_ψi)i,j=1,…,N.

Since the azimuth angle for each blade is given by

_ψj=Ωt+(j−1)2πNj=1,…,N

the following can be written:

_Bψi,j=2Ncos(j−i)2πNi,j=1,…,N.

Consequently,_βGis here related to the overall blade flapping angleβthrough a matrix whose coefficients are constant in time and depend only on the number of blades. In an analogous way,

_β˙G=_Bψ·β˙.

The flapping dynamics equation for each blade is written in the following form:

mβ¨+cβ˙+_cstructΔβ˙+kβ+_{kG βG}+_kβΔβ=_Fβ

where:

• m, c and_Fβare

  m=_Iβc=−_M1Fβ=_M3−_{mb rcg}(_az3−_gz).

• k contains the contributions of aerodynamic and inertial loads to system stiffness:

  k=−_M2+_{mb rcg az2}.

  In a generic flight condition, k, c and_Fβare in general different for each blade.

• _cstructtakes into account the system structural damping. As evaluating the structural damping is extremely difficult, it is a common practice to express this term as a function of other system parameters. In the present model,_cstruct has been expressed as a function of the viscous equivalent critical damping ratio, following the formulation proposed in Reference [22]:

  _cstruct=0.3·2m_kβm.

• _kGcontains the contributions to system stiffness dependent on_βG:

  _kG=_kG+_{mb rcg azG},

  where_kGis the torsional stiffness of the gimbal joint and_azG is z-component of acceleration seen by the blade due to the gimbal motion, function of the angular speed of the rotorΩand the yaw rate of the rigid body, in hub axes,_rh

  _azG=Ω−Ω+_rhrcg.

• _kβcan represent either the stiffness of the hub coning hinge (if present), or the structural stiffness.

By substitutingΔβ=β−_βG and Equation (52) in the flapping dynamics equation, and rewriting everything in matrix form, an equation in the following form is obtained:

Mβ¨+Cβ˙+Kβ=_Fβ.

In which the following terms appear:

• mass matrix:

  M=_IβI

  whereIis the identity matrix.

• damping matrix:

  C=_Ca+_cstruct(I−_Bψ),

  with:

  _Ca=_Ca1 0…00_Ca2 …0⋮⋮⋱⋮00…_{CaN Caj} =−_M1j .

• Stiffness matrix:

  K=_Ka+_kβI−(_kβ−_kG)_Bψ,

  with:

  _Ka=_Ka1 0…00_Ka2 …0⋮⋮⋱⋮00…_{KaN Kaj} =−_M2j +_{mb rcg az2j}.

• Non-homogeneous term:

  _Fβ=_{Fβ1 Fβ2} ⋮_{FβN Fβj} =_M3j −_{mb rcg}(_az3j−_gzj ),

where_M1j ,_M2j ,_M3j ,_az2j,_az3j are computed for the j-th blade using Equations (44) and (27), while_gzj is the gravity acceleration along_zbfor the j-th blade. For a three-bladed rotor such as the XV-15 one, in particular,_Bψreduces to:

_Bψ=23−13−13−1323−13−13−1323

and Equation (63) becomes

_Iβ000_Iβ000_{Iββ¨1β¨2β3}¨+133_c1+_{cstructcstructcstructcstruct}3_c2+_{cstructcstructcstructcstruct}3_c3+_{cstructβ˙1β˙2β3}˙++132_kG+3_k1+_kβkβ−_kGkβ−_kGkβ−_kG2_kG+3_k2+_kβkβ−_kGkβ−_kGkβ−_kG2_kG+3_k3+_kββ1β2β3=_{Fβ1 Fβ2 Fβ3} .

The dynamic problem is thus described by a system of second-order differential equations in which the only states are the flapping angles of each blade, and their derivatives. 2.6. Hub Loads The loads transmitted by each rotor to the airframe constitute the main output of the rotor model. The following components are considered:

• aerodynamic loads;
• blade inertia;
• blade weight;
• forces and moments generated by the flapping hinge.

Once all forces and moments acting on each single blade are determined, the constraint reactions in the flapping hinge can be computed as

R→=−_F→ext−_{mb a→cg}−_F→int,

where_F→extdenotes the external loads (i.e., those related to aerodynamics and gravity),_F→intindicates the internal loads produced by the flapping hinge, and_a→cgis the acceleration of the blade center of gravity. Due to Newton’s third law of motion, the loads transmitted by the blade to the rest of the aircraft, in the Blade reference frame, are simply:

F→=−R→.

Subsequently, the vectorF→obtained for each blade is transported to the center of the rotor hub and referred to the Hub frame of reference. Eventually, the overall hub loads are computed by summing up contributions from each blade (namely_Xh,_Yh,_Zh,_Lh,_Mh,_Nh).

3. Software Implementation 3.1. Tilt-Rotor Simulation Model

The novel model of the flapping dynamics of the gimballed rotor is meant to be part of the comprehensive tilt-rotor mathematical model for real-time flight simulation developed by ZHAW and Politecnico di Torino and implemented in the Research and Didactics Simulator (ReDSim). The comprehensive model has already been introduced in Reference [15] and for the sake of brevity will not be described in detail again in this paper. The tilt-rotor model is developed in Simulink environment and structured in functional blocks accounting for the main components of the aircraft. The novel gimballed rotor model replaces the previous one within the overall model and no major modification to the code are needed. This structure allows to easily introduce new features in the Simulink code and run the new version on the simulator right away with no need for compilation and upload of any embedded software features on dedicated hardware. According to the tilt-rotor configuration, the rotor model is implemented twice in the model, once for the right rotor, once for the left one, being all variable of interest duly transformed from body to rotor reference systems accounting for the nacelle conversion angle.

3.2. Numerical Integration

The mathematical model of the rotor is meant to be implemented in Simulink in discrete form. The integration of the flapping dynamics Equation (63) is performed using a fourth-order Runge-Kutta solver:

_yn+1=_yn+16_k1+2_k2+2_k3+_k4,

with:

_k1=hf(_tn,_yn)_k2=hf_tn+h2,_yn+_k12_k3=hf_tn+h2,_yn+_k22_k4=hf(_tn+h,_yn+_k3)

where h corresponds to the simulation sampling time_Ts. The current implementation uses a fixed sampling rate of 400 Hz for the numerical solution of flapping dynamics equations. As a matter of fact, the numerical stability of the solver is related to the high natural frequency of the blade flapping motion, which is affected by the great stiffness value_Kβ (see Table 1). The fourth-order Runge-Kutta was found to ensure the correct solution of the equations using a sampling rate for which the hardware/software set-up used, can ensure soft real-time. On the other hand, the equations of motion of the aircraft as well as the rotors’ inflow dynamics can be integrated using a discrete Euler backwards method with a constant sampling rate of 400 Hz as well.

3.3. Stand-Alone Model

In order to evaluate the behaviour of the gimballed rotor model and compare it to that of the articulated one, a Simulink model was created to simulate an isolated rotor, mounted on a test rig and subject to the desired external condition (in terms of air density and airstream velocity components) and the desired inputs (in terms of collective, lateral cyclic and longitudinal cyclic blade pitch inputs). A flow chart of the stand-alone rotor model architecture is provided in Figure 6.

4. Results and Discussion

The standalone model was used to compare the response of the previous rotor model (articulated version) and the proposed one, the novel gimballed rotor model. The reference rotor geometry is the XV-15 rotor for both the model options and all data was derived from References [2,6]. Results are presented below in this section.

4.1. Rotor Performance

The proposed gimballed rotor model was compared with the available data in terms of thrust coefficient_CT, power coefficient_CPand Figure of MeritFM . The reference data is derived from the experimental test conducted by NASA on the isolated XV-15 Proprotor and reported in Reference [21]. As shown in Figure 7, the model matches the experimental data for the thrust coefficient very closely, while Figure 8 shows that the power coefficient is slightly over-predicted at high thrust values. As a result the predicted figure of merit (see Figure 9) is slightly lower than what experimentally registered, as well. However, as shown in Figure 10, the results of the proposed rotor model in terms of static thrust and power coefficient of the uninstalled rotor, match what predicted by the reference model and validated against experimental wind tunnel data in Reference [6] The results, however, show a good correlation with the experimental data and further improvements can be made with a finer tuning of the airfoil’s aerodynamic parameters.

4.2. Blades Flapping Motion and Disk Coning Here below the time histories of the flapping motion of the rotor blades are reported, for three different conditions:

• Hover
• Low-speed flight in helicopter mode
• High-speed flight in airplane mode

For the sake of brevity, only the overall flapping angleβof the j-th blade of the rotor and the resulting coning angle_β0 (relative to the sole flapping motion and not accounting for the rotor pre-cone angle) are reported, to prove the thesis of this study. The input parameters are reported in Table 2 and are derived from the trim values of the overall tilt-rotor model in those three specific conditions.

The results are compared with those of the previous version of the model (articulated rotor version). As clearly shown in Figure 11, Figure 12 and Figure 13, the novel gimbal model behaves as expected, therefore a considerable reduction of the coning angle is experienced, if compared to the behaviour of the previous rotor model (articulated version), in virtue of the effect of_Kβ. The coning angle_β0 is obtained averaging the flapping of the N blades of the rotor around 1 revolution. Precisely, for the sake of clarity Figure 14, Figure 15 and Figure 16 show the flapping motion of the single j-th blade, but the other two blades would shown the same motion, with a different phase angle, due to the different azimuthal position. Accordingly, the Simulink stand-alone model is run for each condition, the results are collected in the Matlab workspace and processed with a Matlab script to produce the plots. What reported here below are the time histories of the responses of interest for 1 s of simulation, in which the initial second is omitted not to account for the transient and only show the resulting stabilised response.

4.3. Code Performance in Pilot-In-The-Loop Simulations

The proposed rotor model was implemented in the comprehensive tilt-rotor model and tested in the ReDSim at ZHAW. The overall Flight Simulation Platform and has been described in Reference [15] and no major changes were introduced, except the replacement of the previous rotor model with the new proposed one. To allow pilot-in-the-loop simulations, the Simulink^® model must be able to maintain real-time. This is done with a proprietary C++ S-Function that, for each iteration step, compares the time of the simulation with the time of the computer. For this particular test, the S-Function was called at 400 Hz (every 2.5 ms) and the difference between simulation time and real-time was registered. The data is reported here in Figure 17. The histogram of the time difference in Figure 18 shows that about 99.87% of the samples have less than 2.5 ms of time difference. The mean value of the distribution will also depend on the resolution of the function used to suspend the execution of the current thread until the time-out elapses. This is deemed to be acceptable for pilot-in-the-loop simulation, where the main requirement is represented by the time delay between the pilot input and aircraft response (as per “Allowable airplane response delay” in MIL-F-8785C). The time history in Figure 17 shows a few peaks of values greater than 2.5 ms, one reaching 18 ms. These peaks are assumed to be caused by other services that are running on the Windows computer and will be investigated in the future. Several mitigations are possible such as: removal of windows services, stop antivirus at run time, increase task priority, and so forth. It is obvious from Figure 18 that the execution of one iteration of the model can be achieved in less than 2.5 ms, since the simulation is always able to recover from a delay caused by the operating system.

5. Conclusions and Future Developments

Based on the baseline offered by Ferguson in Reference [6] about the GTRS, a novel, more detailed and efficient model was developed. Compared to the reference, the novel model introduces a higher order flapping dynamic formulation, it accounts for the gimbal constraint to the blades flapping motion, it exploits a dynamic modelling of the rotor’s inflow and offers a more detailed multi-blade description of the rotor’s air-loads. As shown in the previous section of this paper, the novel mathematical formulation allows to increase the accuracy of the results (compared to the previous articulated version of the model) without overly increasing the complexity of the model or compromising its real-time capability. Thus, the novel gimballed rotor mathematical model improves the previous rotor model and the comprehensive tilt-rotor flight simulator of ZHAW and Politecnico di Torino, which can now be used for extensive tilt-rotor piloted simulations, from helicopter mode, through conversion mode to airplane mode. The novel mathematical formulation is believed to offer additional insights in the tilt-rotor real-time simulation topic as well as a new baseline for the development of other models for different aircraft architectures. However, the authors acknowledge that due to the lack of data on the actual coning angles of the XV-15 rotors, a comparison against experimental data could not be performed at this time. Future works will involve an extensive review of the aerodynamic database of the tilt-rotor simulation model aiming for a more complete validation against available data of the XV-15, as well as handling qualities studies performed in the ReDSim and flight control systems research and development activities. Precisely, an extensive handling qualities assessment will be performed according to ADS-33 as well as specific tilt-rotor-based MTEs, as suggested by Padfield in Reference [8]. Eventually, further studies involving tilt-rotor inceptor architectures are envisaged such as combined Collective-And-Throttle-Lever solutions.

Element	Simbol	Unit	Value	Reference
blade	kβ	[Nm/deg ]	244,047.23	[6]
gimbal	kG	[Nm/deg ]	305.05	[6]

Element	Symbol	Unit	Cond. 1	Cond. 2	Cond. 3
True Air Speed	TAS	[kts]	0	60	180
Nacelle Angle	iN	[deg]	90	90	0
Rotor Angular speed	Ω	[rpm]	589	589	517
Hub Inlet X-Velocity	uh	[m/s]	0	30.359	−2.998
Hub Inlet Z-Velocity	uh	[m/s]	0	−5.359	−92.551
Collective pitch input (after mixing unit)	θ0	[deg]	6.98	4.82	28.26
Longitudinal cyclic pitch input (after mixing unit)	B1	[deg]	−1.783	0.262	1.500

Author Contributions

Conceptualization, F.B., P.C. and G.G.; methodology, F.B. and P.C.; software, A.A., F.B. and P.C.; validation, A.A. and F.B.; formal analysis, A.A., F.B.; investigation, A.A. and F.B.; resources, F.B., G.G. and P.C.; data curation, A.A.; writing-original draft preparation, F.B.; writing-review and editing, A.A., P.C. and G.G.; visualization, A.A. and F.B.; supervision, G.G. and P.C.; project administration, G.G. and P.C.; funding acquisition, G.G., P.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

This activity is part of a Doctoral Research Program (Ph.D. in Aerospace Engineering at Politecnico di Torino) and a Master Degree Thesis (M.Sc. in Aerospace Engineering at Politecnico di Torino) supported by Zurich University of Applied Sciences, Centre for Aviation.

Conflicts of Interest

The authors declare no conflict of interest.