[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Next Article in Journal
Autonomous UAV Chasing with Monocular Vision: A Learning-Based Approach
Previous Article in Journal
Measurement and Identification of Flame Describing Function (FDF) Based on Parallel Subsystem Model
Previous Article in Special Issue
Identification of High-Order Linear Time-Invariant Models from Periodic Nonlinear System Responses
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Modeling, Trim, and Linearization of a Side-by-Side Helicopter in Hovering Conditions

by
Francesco Mazzeo
1,
Marilena D. Pavel
2,*,
Daniele Fattizzo
3,
Emanuele L. de Angelis
3,† and
Fabrizio Giulietti
3,†
1
Departmentof Engineering “Enzo Ferrari”, University of Modena and Reggio Emilia, 41121 Modena, Italy
2
Department of Aerospace Engineering, Delft University of Technology, 2628 Delft, The Netherlands
3
Department of Industrial Engineering, University of Bologna, 47121 Forlí, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Aerospace 2024, 11(11), 927; https://doi.org/10.3390/aerospace11110927
Submission received: 1 October 2024 / Revised: 30 October 2024 / Accepted: 7 November 2024 / Published: 9 November 2024
(This article belongs to the Special Issue Vertical Lift: Rotary- and Flapping-Wing Flight)
Figure 1
<p>Side-by-side helicopter (technical data in <a href="#app1-aerospace-11-00927" class="html-app">Appendix A</a>). The rotorcraft has a take-off mass of 20.62 kg, a rotor radius of 0.5 m, and an overall size of 2.5 m.</p> ">
Figure 2
<p>Schematic representation of hub-body and blade frames of reference.</p> ">
Figure 3
<p>Simulink algorithm to partially decouple and solve rotor dynamics in a non-rotating frame of reference.</p> ">
Figure 4
<p>Blade configurations for a clockwise rotor with <math display="inline"><semantics> <mrow> <msub> <mi>N</mi> <mi>b</mi> </msub> <mo>=</mo> <mn>3</mn> </mrow> </semantics></math> and <math display="inline"><semantics> <mrow> <msub> <mi>N</mi> <mi>c</mi> </msub> <mo>=</mo> <mn>3</mn> </mrow> </semantics></math>.</p> ">
Figure 5
<p>Trim algorithm.</p> ">
Figure 6
<p>Pilot controls (<b>left</b>) and rotorcraft attitude (<b>right</b>) in trim condition at variable forward speeds.</p> ">
Figure 7
<p>Rigid body poles affected by lateral and longitudinal disc tilts, according to the numerical (black markers) and analytical (red markers) frameworks.</p> ">
Figure 8
<p>Frequency (<math display="inline"><semantics> <mi>ω</mi> </semantics></math>) of the rigid body poles affected by lateral and longitudinal disc tilts.</p> ">
Figure 9
<p>Rigid body poles with and without the effect of the lead-lag dynamics. The high-frequency poles are represented on the left side of the plot, with a different x-axis scaling.</p> ">
Figure 10
<p>Frequency (<math display="inline"><semantics> <mi>ω</mi> </semantics></math>) of the rigid body poles with the inclusion of a complete flap, dynamic lead–lag, and non-uniform inflow complexities.</p> ">
Figure 11
<p>Rigid body poles in the cases of uniform (black markers) and non-uniform dynamic inflow (green markers).</p> ">
Figure 12
<p>Mode participation with non-uniform inflow.</p> ">
Figure 13
<p>Main rotor poles in the complex plane.</p> ">
Figure 14
<p>Mode participation of the rotor-coupled poles.</p> ">
Figure 15
<p>Mode participation of the uniform inflow poles.</p> ">
Versions Notes

Abstract

:
In the present paper, a flight dynamics model is adopted to represent the trim and stability characteristics of a side-by-side helicopter in hovering conditions. This paper develops a numerical representation of the rotorcraft behavior and proposes a set of guidelines for trimming and linearizing the highly coupled rotor dynamics derived by the modeling approach. The trim algorithm presents two nested loops to compute a solution of the steady-state conditions averaged around one blade’s revolution. On the other hand, a 38-state-space linear representation of the helicopter and rotor dynamics is obtained to study the effects of flap, lead–lag, and inflow on the overall stability. The results are compared with an analytical framework developed to validate the rotorcraft stability and compare different modeling approaches. The analysis showed that non-uniform inflow modeling led to a coupled longitudinal inflow–phugoid mode which made the vehicle prone to dangerous instabilities. The flap and lead–lag dynamics introduced damping in the system and can be considered beneficial for rotor dynamics.

1. Introduction

Considerable research has been performed over the past 50 years to improve the fidelity of mathematical models for rotorcraft flight dynamics [1,2,3]. In this respect, different approaches have been developed, summarised in order of complexity: (1) the most classical one formulates rotor nonlinearity analytically by using blade element theory to calculate forces and moments generated by the rotor system [4]; (2) a second type of modeling technique considers numerical differential equations for both rigid body motion and blades’ flexibility, proposing an integral nonlinear formulation of the rotorcraft loads [5,6]; (3) a third type employs the principles of multi-body dynamics to describe rotorcraft behavior [7,8], an example of which can be seen in the known flight dynamics software FLIGHTLABTM (Release v3.6.0) [9]; and (4) a more refined approach consists of coupling computational fluid dynamics (CFD) with computational structural dynamics (CSD) solutions [10]. While each approach has proven its applicability for rotorcraft, most of the above-mentioned models use the common ’time scale separation’ assumption when calculating the trim and performing the linearization. This means that the rotor (i.e., inflow, flapping, and lead–lag) dynamics are typically stable and significantly faster than that of the rigid body, and a residualization method can be applied by employing singular perturbation theory [11]. This approach removes the need to estimate higher-order rotor dynamic states while retaining information on the influence of the residualized dynamics on the slower rigid body dynamics [12]. While the residualization procedure is commonly applied for the design of control systems, since it allows the state-space system to be simplified by separating slower and faster dynamics, this procedure should be carefully implemented in highly nonlinear systems such as rotorcraft where the time scale separation assumption is not always valid [13].
The main goal of this paper is to fundamentally understand the effects of blade flexibility on rigid body dynamics by implementing a novel method for performing trim and linearization in a highly coupled nonlinear system. Specifically, a side-by-side helicopter was studied and flight-tested by SAB Group s.r.l. and the University of Bologna to understand the potential of this configuration for future VTOL developments in Urban Air Mobility (UAM). A first modeling approach was developed for a preliminary analysis [14] involving an analytical representation of the vehicle flight dynamics. This framework, however, highlighted a disagreement between the predicted stability and the rotorcraft behavior. The 14 degrees of freedom (dof)analytical model predicted overall stable behavior, except for a slightly unstable, oscillatory phugoid mode. On the other hand, a qualitative analysis of the preliminary flight tests revealed a tendency of the rotorcraft to diverge from its hovering condition, behaving as a longitudinally unstable system with a real, positive pole. As expected, a typical unstable phugoid in rotorcraft flight dynamics appeared as a divergence of the pitch angle [15]. For this reason, a higher-complexity framework has been developed and described in this paper by including additional rotor effects and a numerical representation of the aerodynamic blade loads. It will be demonstrated that, while a numerical approach eliminates the need for complex algebraic expansions of the blade’s nonlinear structural dynamics, it also requires a methodology for in-place linearization and trim analysis. To this extent, three major issues were identified when dealing with numerical formulations:
  • The model describes the aerodynamic loads on each single blade by applying a blade element approach in a rotating frame of reference, and sums up their contributions in the two rotors. This single-blade representation leads to a loss of physical meaning when developing linear state-space representation;
  • Rotor dynamics obtained by a Lagrangian approach on the single-blade representation produce a high level of coupling between flap and lead–lag dynamics;
  • The presence of complex, nonlinear terms derived by numerical integration makes the linearization process impractical with an explicit version of the equations.
In order to achieve the main purpose of this paper, a set of guidelines has been provided to overcome the listed issues and provide a 38-state-space representation of the side-by-side helicopter flight dynamics. The methodology presented in this paper will provide motivation for a more in-depth analysis of mathematical modeling techniques to be used in future advanced rotorcraft configurations [16].
This paper is organized as follows: in the modeling section, the numerical model is presented for the case studied in this paper, with a focus on the rotor dynamics and its suitable representation. Trim and linearization methodologies are proposed for this type of mathematical representation and equilibrium results are summarized and compared with an analytical framework. The following section shows the results of the stability of rigid body dynamics, as affected by different levels of complexity: the effects of flap, lead–lag, and non-uniform inflow will be addressed in this section. Finally, the particular rotor dynamic characteristics of this configuration are reproduced by including the collective, advancing, and regressive flap–lead–lag modes and their coupling. This work terminates with some concluding remarks. Rotorcraft data and matrix calculations are reported in the Appendices A and B.

2. Materials and Methods

2.1. Modeling Overview

This case study is of a symmetric side-by-side helicopter, whose data (and their notation) are recalled in Appendix A. The rotorcraft has two counter-rotating, ducted, and synchronous rotors which were modeled as isolated ones. The rotors are located on each side of the fuselage and do not present any intermeshing section. For this reason and by considering the geometry of the case study (Figure 1), it is reasonable to assume that the aerodynamic interactions between the rotors can be neglected. The rotors rotate at a constant angular speed and they are actuated with a typical helicopter swashplate. The system is thus over-actuated, meaning that the single-rotor controls have to be mixed in order to reduce them to classical helicopter ones, i.e., the global collective pitch ( θ 0 ) for vertical motion, lateral ( A 1 s ) and longitudinal ( B 1 s ) cyclics for in-plane movement, and yaw control ( Δ B 1 s ) are obtained as follows:
θ 0 = θ 0 ( 1 ) = θ 0 ( 2 )
A 1 s = A 1 s ( 1 ) = A 1 s ( 2 )
B 1 s = B 1 s ( 1 ) + B 1 s ( 2 ) 2
Δ B 1 s = B 1 s ( 1 ) B 1 s ( 2 ) 2
where the superscript numbers (1) and (2) identify the rotor. The mathematical framework adopted in this paper is a 24 dof numerical model, which describes the flight dynamics of the rotorcraft under the hypothesis of a rigid body moving on a flat and non-rotating Earth. Second-order coupled flap–lag dynamics was integrated into the main rotor description and a non-uniform, first-order, dynamic inflow provided by Peters et al. [17] was included as well. The fuselage aerodynamic loads ( F f and M f ) are represented with a flat-plate approximation by using the software OpenVSP (Release 3.34.0) to calculate the parasite drag of the frontal, side, and top equivalent areas [14]. The effect of the shrouds is taken into account as well with the integration of semi-analytical results produced by Bourtsev et al. [18] on the Leishman representation of the thrust with ducted fans [19]). In their work, Bourtsev et al. provide a simple and analytical model for calculating the wake contraction factor ( a w ) of a fan-in-fin rotor, where the effect of the shroud depends on the characteristics of its geometry, such as lip radius, type, or depth. The total thrust of each rotor (T) is given by the rotor force component perpendicular to the tip-path plane ( T r ) plus the contribution of the shroud, such that
T = T r + T d where T d T = 1 1 2 a w
and a w is derived from Bourtsev’s results. For the sake of simplicity, the model does not take into account the aerodynamic interactions between rotor and shroud, and the contraction factor is kept constant for the duct geometry. The latter approximation makes the model suitable for hovering studies, while it lacks proper validation in the case of forward flight. The rotorcraft’s body frame of reference (centered in the center of gravity CG) and rotor directions are shown in Figure 1.
The numerical modeling approach described in this work allows for isolating the influence of rotor dynamic effects (flap, lead–lag, and inflow) on the overall rotorcraft stability. This single-blade representation can be employed for a deeper study of the blade design and vibration analysis, as well as to consider particular chord distributions and unbalanced rotors. The results are compared with an already developed and validated analytical model described in Ref. [14]. The model has 14 dof and is based on Talbot’s formulation for classical helicopter configuration [4]. The latter provides analytical expressions of the total forces and moments generated by the rotor (so not requiring numerical integration) and includes second-order tip-path plane dynamics and a uniform dynamic inflow. A comparison between the two modeling approaches is performed in this paper to validate the results obtained with the numerical approach at an equivalent level of complexity.

2.2. Frames of Reference

The mathematical framework adopts three main frames of reference (f.o.r.), as shown in Figure 2. The first one is the body f.o.r. ( [ x B y B z B ] ), a right-handed f.o.r. identified by the superscript B in the equations. It is centered in the center of gravity (CG) and coherent with the rotorcraft attitude: the lateral axis y B is directed towards the right-hand side rotor, while z B is directed towards the underneath of the vehicle. The second f.o.r. to be described is the local hub-body frame ( [ x H B y H B z H B ] ), identified by the superscript HB in the text. This right-handed f.o.r. is centered on each main rotor center of rotation, with the vertical axis z H B directed downwards and parallel to the rotor shaft. The horizontal axis x H B is instead perpendicular to the rotor shaft and pointed towards the nose of the rotorcraft. The third f.o.r. to be described is the local blade frame of reference (superscript bl): this is centered on the aerodynamic center of each blade section, coherent with the azimuthal position of the blade ( ψ ), and its flap ( β ), lead–lag ( ξ ), and twist ( θ b ) angles. The orientation of this f.o.r. is depicted in Figure 2. The Figure also provides the positive direction of the blade’s angles in a clockwise rotor. In general, a positive ψ starts from the rotorcraft’s tail and increases with the rotor’s sense of rotation, a positive β is defined when the blade flaps up, and a positive ξ is defined when the blade lags with respect to the sense of rotation. The general rotation matrix between two of these f.o.r. is represented by the letter R.

2.3. Forces and Moments

The main rotor forces and moments were modeled with a blade element approach by dividing each blade into a finite number of sections governed by 2D aerodynamics, and by integrating the infinitesimal contributions along the span. The main rotor forces are divided into three contributions: (1) aerodynamic ( F a ), (2) inertial ( F i ), and (3) centrifugal ( F c ). The main rotor moments are provided by their (1) aerodynamic ( M a ), (2) inertial ( M i ), and (3) flap hinge stiffness ( M s ) components. Considering the single-blade contribution of a Γ -rotating rotor [20] (where Γ = 1 for a counter-clockwise rotor and Γ = 1 for a clockwise rotor), the aerodynamic force in the hub-body frame of the i-th blade (subscript i) is computed as
F a i H B = r c R e 0 R H B b l ( d L b l + d D b l ) d r
where the infinitesimal force given by the sum of lift and drag ( d L and d D ), rotated from the blade to the hub-body frame of reference, is integrated along the blade span, from the root cutout to the blade tip. R H B b l = R H B b l ( ψ i , β i , ξ i , θ b i ) is the rotation matrix between the two mentioned frames; it is specific for each blade section and depends on its radial (r) and azimuthal ( ψ i ) locations. The infinitesimal lift and drag contributions can be summarized as
d L b l = 1 2 ρ | V a | 2 C L ( α i , R e i ) c Γ sin α i 0 cos α i ,
d D b l = 1 2 ρ | V a | 2 C D ( α i , R e i ) c Γ cos α i 0 sin α i
where C L and C D are the aerodynamic lift and drag coefficients of the specific blade section subjected to the Reynolds number R e = c | V a | ν and angle of attack α = arctan V a z Γ V a x . The airfoil’s aerodynamics were modeled by employing look-up tables of conventional NACA profiles and by extending them with semi-analytical formulas developed by Viterna [21]. The local twist angle is instead the parameter where the pilot controls enter into the computation. In particular,
θ b i = θ 0 + A 1 s cos ( ψ i + ψ θ ) + B 1 s sin ( ψ i + ψ θ ) + r R θ w K t ξ ξ i K t β β i
where θ 0 , A 1 s , and B 1 s are the collective pitch, lateral, and longitudinal cyclic input applied to the j-th rotor. The parameter ψ θ is used to consider the swashplate’s phase angle in the local twist definition. The aerodynamic velocity at the blade section V a = [ V a x V a y V a z ] is instead computed from the inertial velocity ( V I H B , see Ref. [6]) as
V a = R b l H B ( V I H B 0 0 v i V g )
where V g is the gust velocity in the hub-body frame. The local inflow velocity is
v i = Ω R λ 0 + λ s r R sin ψ i + λ c r R cos ψ i
and its coefficients λ = [ λ 0 λ s λ c ] are derived from Peters et al. [17]. Each blade is also subjected to inertial and centrifugal forces, which are, respectively, modeled as
F i i H B = m b η ξ ξ i ¨ Γ cos ψ i sin ψ i 0 , and F c i H B = 1 2 m b ( Ω R ) 2 Γ sin ψ i cos ψ i 0
where η ξ = R e P e L . The expressions of F i i and F c i were provided by Taamallah [6] and re-arranged for a Γ -rotating rotor. The aerodynamic moment in the hub-body frame of the i-th blade is computed as the integral of the infinitesimal contributions at each blade section, given by the cross product between the section position vector with respect to the hub center and the aerodynamic force. In mathematical terms,
M a i H B = r c R e 0 r × [ R H B b l ( d L b l + d D b l ) ] d r
Additional moments come from the flap hinge stiffness and inertial contributions and depend on the rotor disc tilt coordinates a 1 and b 1 (see Section 3.2). These are, respectively,
M s i H B = 1 1 e 0 R N b K s β 2 Γ b 1 a 1 0 , and M i i H B = N b 2 m b e 0 r G Ω 2 Γ b 1 a 1 0
Finally, the total external force and moment in the body f.o.r. applied by the j-th rotor is computed by the sum of the contributions of each blade plus the moment derived by the cross product between the rotor hub position and the external force. In particular,
F i B = R B H B i = 1 N b ( F a i + F c i + F i i )
M j B = r H B × F j B + R B H B i = 1 N b ( M a i + M s i + M i i )
The two rotors share the same design and operate at equal and constant angular speed; their only differences are their location with respect to the rotorcraft’s CG and the sense of rotation. The total external force and moment acting on the VTOL, expressed in body f.o.r., are
F = F ( 1 ) + F ( 2 ) + F f
M = M ( 1 ) + M ( 2 ) + M f

2.4. Rotor Flap–Lag Dynamics

The rotor dynamics is described by a system of second order ordinary differential equations (ODE), derived with a Lagrangian method by Tamallah (Ref. [22]) for a P-L-F (Pitch-Lag-Flap) small-scale helicopter rotor. The system has 6 degrees of freedom for each rotor and describes the coupled flap and lead-lag dynamics in a rotating frame of reference. In particular, for the i-th blade,
β i ¨ = 1 A 1 ( D i β ξ ˙ i + Q i β F i β ) ξ i ¨ = 1 A 2 ( D i ξ β ˙ i + Q i ξ F i ξ )
where A 1 = m b ( R e 0 ) 2 3 and A 2 = m b ( e F 2 + 2 e F r G + ( R e 0 ) 2 3 ) are constants which depend on the blade’s inertia. The terms D i β and D i ξ are nonlinear coupling terms whose expression can be found in Ref. [6] and is specific for the i-th blade. Similarly, Q i F i is the excitation of each ODE and represents the generalized force that creates the flapping and lead–lag motion [6]. For the sake of this paper, the flap–lag equations derived in this form present some criticalities: the equations are highly nonlinear and include numerical integration into the Q terms. In addition, the rotating frame representation is suitable for the rotorcraft simulation but at the same time, it loses physical meaning when studying the rotor stability and its regressive and advancing modes. In order to overcome these problems, Equation (19) was adapted to represent the tip-path plane and rotor center of gravity dynamics, which are indeed described by the flapping and lead-lag coefficients, in a partially decoupled system [23]. As a first step, Equation (19) was re-written in a non-rotating frame by describing the rotor flap and lead-lag angles with the Coleman representation for a 3-bladed rotor [24]. In particular, for the i-th blade,
β i = a 0 a 1 cos ψ i b 1 sin ψ i β i ˙ = a 0 ˙ a 1 ˙ cos ψ i b 1 ˙ sin ψ i + a 1 Ω sin ψ b 1 Ω cos ψ β i ¨ = a 0 ¨ a 1 ¨ cos ψ i b 1 ¨ sin ψ i + ( a 1 Ω 2 2 b 1 ˙ Ω b 1 Ω ˙ ) cos ψ i + ( b 1 Ω 2 + 2 a 1 ˙ Ω + a 1 Ω ˙ ) sin ψ i ξ i = ξ 0 ξ c cos ψ i ξ s sin ψ i ξ i ˙ = ξ 0 ˙ ξ c ˙ cos ψ i ξ s ˙ sin ψ i + ξ c Ω sin ψ ξ s Ω cos ψ ξ i ¨ = ξ 0 ¨ ξ c ¨ cos ψ i ξ s ¨ sin ψ i + ( ξ c Ω 2 2 ξ s ˙ Ω ξ s Ω ˙ ) cos ψ i + ( ξ s Ω 2 + 2 ξ c ˙ Ω + ξ c Ω ˙ ) sin ψ i
where the tip-path plane dynamics are described by the cone-shapedrotor coordinates, i.e., the coning angle a 0 and the lateral and longitudinal flap coordinates a 1 and b 1 . The motion of the rotor center of gravity is instead described by the lead–lag coordinates in the non-rotating frame, i.e., the collective ξ 0 , and the advancing/regressive lag, ξ c and ξ s . By substituting Equation (20) into (19), the rotor dynamics can be then expressed in the following form:
β ¨ = K 1 β β ˙ K 0 β β + C 1 β ξ ˙ + C 0 β ξ + E β ξ ¨ = K 1 ξ ξ ˙ K 0 ξ ξ + C 1 ξ β ˙ + C 0 ξ β + E ξ
where β = [ a 0 a 1 b 1 ] and ξ = [ ξ 0 ξ c ξ s ] . The expression of the K, C, and E terms can be found in Appendix B. In particular, K 0 and K 1 are constant matrices (at each time step) derived from the frame transformation, while C 0 and C 1 are representative of the level of coupling between flap and lead–lag dynamics. These matrices are highly nonlinear, and an explicit transformation from a rotating to a non-rotating frame was not feasible. This criticality will be addressed in Section 2.6 by approximating them with Taylor expansion. E β and E ξ are instead the excitation of the second-order ODE, represented by the flapping and lead–lag moments computed with numerical integration.
The system in Equation (21) represents the rotor dynamics in a non-rotating frame and can be solved with classical numerical methods [25]. However, the main goal of this paper is to fundamentally understand the effect of rotor dynamics on the overall stability of the rotorcraft and identify the major cause of longitudinal instabilities. For this reason, a highly coupled system is not suitable for isolating the influence of flap and lead–lag and a partially decoupled solution has to be found. To overcome this issue, the dynamic was solved with a Simulink scheme depicted in Figure 3. At each time step, t, the flap and lead–lag dynamics equations (Equation (21)) were solved separately by approximating the coupling and forcing terms with the values of the rotor flap and lead–lag coordinates at the previous time step ( t Δ t ). In the algorithm, ψ = [ ψ 1 ψ 2 ψ 3 ] is the vector of a single blade’s azimuthal coordinates used to compute the constant matrices ( K 0 and K 1 ). With the same azimuthal positions and the single blade flap and lead–lag angles, derived from the inverse Coleman transformation, the C 0 , C 1 , and E terms of the equations are computed and used to solve the dynamics. The algorithm allows for solving the rotor equations with a partial decoupling so one dynamic or the other can be excluded by simply cutting off its branch and isolating its effect on the overall system.

2.5. Trim

The side-by-side helicopter dynamics can be described by a set of 24 equations, involving 6 for the rigid body motion, 12 for the rotor, and 6 for the dynamic inflow of the two main rotors. Computing the trim condition means solving the system in a steady-state configuration, thus setting all the first- and second-order derivatives to equal zero. The full set of equations to be solved is
m t o U ˙   = m t o ω × U + F g + F I ω ˙   = M ω × I ω β ¨ ( 1 / 2 ) = [ K 1 β β ˙ K 0 β β + C 1 β ξ ˙ + C 0 β ξ + E β ] ( 1 / 2 ) ξ ¨ ( 1 / 2 ) = [ K 1 ξ ξ ˙ K 0 ξ ξ + C 1 ξ β ˙ + C 0 ξ β + E ξ ] ( 1 / 2 ) λ ˙ ( 1 / 2 ) = Ω ( M λ L 1 L 2 ) 1 λ + Ω M λ 1 C aero where U ˙ = ω ˙ = 0 ; β ¨ = β ˙ = 0 ; ξ ¨ = ξ ˙ = 0 ; λ ˙ = 0 ;
and F g is the weight vector force. The inflow dynamics are provided by Peters [17], and the expression of the matrices M λ , L 1 , and L 2 can be found in this reference, as well as C a e r o . The latter is a vector of the aerodynamic coefficients of thrust, rolling, and pitching moments. The calculation of the equilibrium conditions with a numerical, single-blade model requires a nonconventional approach. The presence of numerical integrations in the mathematical framework makes it very impractical to solve the system with a symbolic routine. In addition, with the rotor loads being time-dependent, the trim solution is not unique for each flight condition. Having a single-blade representation of the rotor disc leads to a dependency between the rotor dynamics and the current blades’ configuration ( ψ ), and it has been observed that the tip-path plane and rotor CG coordinates depend on the blade location at which the numerical algorithm has started. Flap wobbling and rotor vibrations are observed such that β and ξ change during one revolution. For these reasons, a general steady-state solution to the trim problem which is able to keep the rotorcraft in equilibrium cannot be found if considering a unique blade setup. In order to find a condition that minimizes the average rotor acceleration in each blade ψ location, the trim has to be computed based on an average rotor force and moment in a N c number of blade symmetric configurations. For example, by setting the parameter N c = 3 , the trim is solved as an average between the equispaced blade configurations depicted in Figure 4.
In general, considering ψ 0 = [ ψ 1 ψ 2 ψ 3 ] as the vector containing the initial position of each blade, the k-th configuration can be computed as
ψ k = ψ 0 + k 2 π N b N c where k = 0 , 1 , 2 , . . . N c 1 ;
The trim algorithm is described in Figure 5. The solution to the problem is found by iteratively solving the steady-state flap–lag equations, the steady-state 6 dof equations of motion (EoM), and the non-uniform steady-state inflow. In order to help the convergence of the algorithm, it was observed that the iterations should be split into an internal loop, which solves the EoM system and the steady inflow for fixed rotor coordinates, and an external loop which solves and updates at each step the flap and lead-lag values. Each block of nonlinear equations was solved with a Newton–Raphson algorithm. The algorithm works as follows:
(a)
The algorithm starts with the initial condition of the 26 unknowns (6 flap coefficients, 6 lead–lag coefficients, 6 inflow coefficients, 2 collective pitches, 4 cyclic controls, roll ( ϕ ), and pitch ( θ ) attitudes).
(b)
The rotor dynamic equations are solved for an N c finite number of equispaced blades’ configurations and the disc tilt and rotor CG coordinates are derived for each of them.
(c)
The average values of β k and ξ k are computed.
(d)
Knowing the average rotor coordinates, the main rotor forces and moments are computed with Equations (15) and (16). Again, this computation is repeated for N c blade configurations.
(e)
The average main rotor loads are computed.
(f)
With the average rotor disc tilts, CG coordinates and main rotor loads, the system of steady-state equations of motion is solved. Two additional Equations, (1) and (2), are included to take the control mix into account and close the system.
(g)
With the new state obtained on the averaged EoM, the inflow dynamics is solved. The equations of the 3-states inflow model are summarised in Peters et al. [17];
(h)
The internal loop that solves EoM and inflow is repeated until the stopping criterion (maximum relative error between the rotorcraft state, inflow, and controls at each internal iteration) is satisfied.
(i)
Once the internal loop has converged, the rotorcraft state is updated and a new external iteration is started. The external loop is repeated until the stopping criterion (maximum relative error between the TPP and rotor CG coordinates at each external iteration) is satisfied.
The results at a variable forward speed are reported in Figure 6, together with a comparison with the 14 dof analytical model previously developed by the authors [14]. The curves are typical helicopter trim curves, with a collective pitch control of approximately 10 at hover which has a decreasing trend as the forward speed increases (up to the minimum power condition). A positive longitudinal cyclic control (push effort for the pilot) is adopted at a positive forward speed together with a nose-down, negative-pitch attitude. The two models have a good agreement in both attitude and control results, with a small discrepancy in the pitch angle.

2.6. Linearization

In order to study the flight dynamic properties and open-loop stability of a rotorcraft, the numerical model has to be linearized around a specific trim condition. While the development of state-space linearized models for the analytical frameworks has already been described in several works [26,27,28], the process of linearizing a numerical, single-blade model, such as the one presented in this paper, still needs to be properly stated. The higher degrees of freedom, the numerical integration of forces and moments along each single blade, and the high level of coupling between flap and lead–lag dynamics create challenges in the derivation of an approximated linear model that describes the rotorcraft behavior around specific flight conditions. In this section, a methodology to overcome this problem and to develop a 38-state linear model is presented. The 38 ordinary differential equations describing the rigid body motion, rotor, and inflow dynamics are
m t o U ˙ = m t o ω × U + F g + F I ω ˙ = M ω × I ω ϕ ˙ = p + q sin ϕ tan θ + r cos ϕ tan θ θ ˙ = q cos ϕ r sin ϕ β ˙ ( 1 / 2 ) = d β ( 1 / 2 ) d β ˙ ( 1 / 2 ) = K 1 β d β K 0 β β + C 1 β d ξ + C 0 β ξ + E β ( 1 / 2 ) ξ ˙ ( 1 / 2 ) = d ξ ( 1 / 2 ) d ξ ˙ ( 1 / 2 ) = K 1 ξ d ξ K 0 ξ ξ + C 1 ξ d β + C 0 ξ β + E ξ ( 1 / 2 ) λ ˙ ( 1 / 2 ) = Ω ( M λ L 1 L 2 ) 1 λ + Ω M λ 1 C a e r o
where the first 8 equations represent the rigid body motion plus the Euler angle definition. The second-order flap and lead–lag dynamics described in a non-rotating frame of reference in Equation (21), with a set of 6 equations per each rotor, are reduced to 24 first-order ODEs by applying a separation of variables. In particular, d β = [ a ˙ 0 a ˙ 1 b ˙ 1 ] and d ξ = [ ξ ˙ 0 ξ ˙ c ξ ˙ s ] are the variables added to reduce the second- to first-order ODEs. The system of equations is linearized and reduced to the form
x ˙ = A x + B τ
where A, B, x , and τ are the state and control matrices and the state and control vectors. In particular,
x = u w q θ v p ϕ r β ( 1 ) β ( 2 ) d β ( 1 ) d β ( 2 ) ξ ( 1 ) ξ ( 2 ) d ξ ( 1 ) d ξ ( 2 ) λ ( 1 ) λ ( 2 )
τ = θ 0 A 1 s B 1 s Δ B 1 s
In order to reduce Equation (24) to (25), the small perturbation theory is applied. Each general state x is described as a perturbation from the trim condition, thus x = x 0 + δ x , where δ x is a small perturbation and the superscript 0 identifies the trim values. The second-order terms are neglected and the equations are written in terms of small perturbations, to describe the linear dynamics of the system around a specific trim condition. The external forces and moments ( F = [ X Y Z ] and M = [ L M N ] ) and the excitation of the rotor flap and lead–lag dynamics ( E β and E ξ ) are the forcing terms of the respective ODEs. These terms are represented by analytic functions of the disturbed motion variables and their derivatives using Taylor’s theorem. They can be expressed as their trim value plus a sum of contributions from each state and control variable, such that
F = F 0 + F u u + F w w + . . . + F Δ B 1 s Δ B 1 s
where the first-order derivative about a generic perturbation δ x of the x state is computed numerically with a central finite difference formula as
F x = F x i | x 0 = F ( x 0 + δ x ) F ( x 0 δ x ) 2 δ x
The same approach is also applied to the inflow and the coupling terms in the rotor dynamics equations ( C 0 and C 1 ) which are nonlinear terms representing the coupling between the flap and lead–lag dynamics of each rotor. In general, the forcing and coupling terms in the rotor equations (C and E ) are nonlinear functions of the rotor states themselves and involve numerical integration for their computation. This makes it impossible to rearrange these terms in an explicit form with the rotor states, and Taylor’s approximation becomes a viable solution for their linearization. The state matrix A will indeed be a function of the trim condition and the stability derivatives which describe the linear response of the perturbed states. As a matter of fact, with the presented methodology, the stability of the system is described not only by the classical rigid body derivatives but also by the coupling body–rotor–inflow derivatives and the ones related to the flap and lead–lag forcing and coupling terms. The general structure of the A matrix is as follows:
A = A B ( 8 × 38 ) A R ( 24 × 38 ) A λ ( 6 × 38 ) = A B B ( 8 × 8 ) A B R ( 8 × 24 ) A B λ ( 8 × 6 ) A R B ( 24 × 8 ) A R R ( 24 × 24 ) A R λ ( 24 × 6 ) A λ B ( 6 × 8 ) A λ R ( 6 × 24 ) A λ λ ( 6 × 6 )
where A B , A R , and A λ are, respectively, the linearized versions of the rigid body, rotor, and inflow dynamics, and they are composed of a set of 3 submatrices for each of them. Considering, for example, the rigid body dynamics, A B B represents the 8-state linear model of the rigid body equations of motion (excluding the heading contribution), while A B R and A B λ represent the coupling effect of the flap, lead–lag, and inflow dynamics on the body motion. On the other hand, A R B for example, represents the coupling effect of the rigid body motion on the rotor dynamics (flap and lead–lag). While the expression of A B B and the normalization of the stability derivatives can be found in Padfield [15], A B R (and similarly A B λ ) is made by the stability derivatives derived by the total forces and moments, perturbed with a flap or lead–lag perturbation, i.e.,
A B R = X a 0 ( 1 ) X a 1 ( 1 ) X b 1 ( 1 ) X ξ 0 ( 1 ) Z a 0 ( 1 ) Z ξ 0 ( 1 ) M a 0 ( 1 ) M ξ 0 ( 1 ) 0 0 L a 0 ( 1 ) L ξ 0 ( 1 ) 0 0 N a 0 ( 1 ) N ξ 0 ( 1 )
The flap and lead–lag dynamics, instead, are represented by the linearized version of the rotor dynamics equations, i.e., A R . Similarly to the body dynamics, A R is made of pure rotor–rotor dynamics ( A R R ) plus the coupling effects of body and inflow on the rotor itself. The structure of the matrix A R is derived from the linearization of Equation (21), and is structured as follows:
A R = Z 24 × 8 Z 6 × 6 I 6 Z 6 × 6 Z 6 × 6 K 0 β K 1 β C 0 β , 0 C 1 β , 0 Z 6 × 6 Z 6 × 6 Z 6 × 6 I 6 K 0 ξ K 1 ξ C 0 ξ , 0 C 1 ξ , 0 Z 24 × 6 + Z 6 × 1 C 0 x β ξ 0 + E x β ( 1 ) C 0 x β ξ 0 + E x β ( 2 ) Z 6 × 1 C 0 x ξ β 0 + E x ξ ( 1 ) C 0 x ξ β 0 + E x ξ ( 2 )
where C 0 β , 0 and C 0 ξ , 0 are the coupling coefficients of the flap and lead–lag equations evaluated in the trim condition, while C 0 x β , C 0 x ξ , E x ξ , and E x ξ are the stability derivatives of the coupling and forcing terms, computed with respect to the generic perturbed state x as mentioned by Equation (29). Symbols Z i × j and I i represent, respectively, a zero matrix of dimensions i × j and an identity matrix of dimensions i × i . All of these terms are derived for both the main rotors 1 and 2, and thus they are structured as, for example,
K 0 β = K 0 β ( 1 ) Z 3 × 3 Z 3 × 3 K 0 β ( 2 )
It can be observed that A R is split, for the sake of clarity, into two matrices: the first one depends exclusively on the trim condition, while the second depends on the stability derivatives. It is indeed a 24 × 38 matrix and the general notation applied for both terms is with the x-perturbed state in the last subscript. In addition, β 0 and ξ 0 are the flap and lead–lag coordinates at the trim, considering the average solution found in the trim section. In general, in the i-th column of this matrix, the derivatives are computed for the i-th perturbed state and for the rotor indicated by the superscript. A similar methodology was applied to the inflow dynamics to derive the linearized version of the inflow equations represented by the matrix A λ .
In order to study the open-loop stability of the rotorcraft, an analysis of the eigenvalues and eigenvectors of the matrix A is performed. In general, 38 different poles were identified and related to the rigid body, flap, lead–lag, and inflow dynamics, and the effects of different degrees of complexity were studied.

3. Results and Discussion

3.1. Rigid Body Dynamics

The rigid body stability of the side-by-side helicopter in hovering conditions is studied by analyzing the effects of single degrees of complexity. In particular, the investigation starts by comparing the numerical and analytical mathematical frameworks in the same, standard and simplest configuration, to validate and highlight the differences between the two formulations. The analytical representation was described in Ref. [14] and validated by comparing the implementation of the main rotor forces and moments with a classical helicopter configuration. The comparison showed a good agreement and a reasonable physical meaning of the results for a low-complexity framework. The first comparison is made at the lowest level of modeling complexity, hence with:
  • A uniform and dynamic inflow ( λ = [ λ 0 0 0 ] in both rotors);
  • Uniform flap dynamics ( β = [ a 0 0 0 ] in both rotors);
  • Lead–lag neglected ( ξ = [ 0 0 0 ] in both rotors).
Starting from this condition, different levels of complexity are included, and their effects on the rotorcraft stability (system poles and eigenvectors) are pointed out. At first, complete flapping dynamics are considered in both models, then the lead–lag effect, and finally, the influence of a dynamic, non-uniform inflow.
Figure 7 represents the rigid body poles on the complex plane. The two different modeling approaches are distinguished with red and black markers, while the plain/empty symbolsidentify the level of complexity. The first comparison is made with uniform flap dynamics, meaning that lateral and longitudinal disc tilts are neglected and β = [ a 0 0 0 ] . On the other hand, the plain symbols introduce the effect of a 1 and b 1 in the rigid body stability and, in this case, β = [ a 0 a 1 b 1 ] . The mode identification is made through eigenvector analysis and by recalling classical aircraft nomenclature for low- and high-frequency modes. The frequency of the poles is represented in Figure 8, highlighting a clear difference between six low-frequency modes and two higher-frequency ones. Two couples of oscillatory dynamics identify slightly unstable, low-frequency phugoid and dutch-roll modes. The low-frequency behavior is also characterized by stable spiral and heave subsidence. From the high-frequency point of view, a couple of stable short-period and roll dynamics are reported. The plot shows a good agreement between the analytical and numerical representations, contributing to the validation of the second. The typical predicted modes are similar and the same stability properties arise. The complete flapping dynamics represented by the solid markers in the plot introduce a first level of complexity, identified as second-order flapping dynamics with lateral and longitudinal disc tilts ( a 1 and b 1 ). The two models behave similarly, with the tip-path plane dynamics affecting the oscillatory poles by increasing their frequency and adding lateral and longitudinal damping to the high-frequency system (short period and roll). This has an overall beneficial effect on the rotorcraft stability. The short period and roll frequency are directly linked to the lateral and longitudinal damping derivatives M q and L p [15] which are reported in Table 1 for the two cases, together with the rotor disc coordinates with a perturbed pitch rate q. The rotorcraft reacts to a disturbance in angular rate with a damping moment around the perturbed axis. The different modeling approach to the flap dynamics leads to different disc tilts which are directly linked to the damping derivative. Indeed, the flap hinge stiffness and inertial moments described in Equation (14) are proportional to a 1 and b 1 which are, for a small perturbation of pitch rate (the same happens with a perturbation of roll rate p), different in the two models (see Table 1). The discrepancy is, in absolute values, quite small, and the disc tilts maintain the same order of magnitude, as well as the damping derivatives.
A second level of complexity is the lead–lag, modeled in the numerical framework. The methodology presented in this work to linearize numerical models allows us to partially decouple the flap and lead–lag dynamics which are represented by a nonlinear system of second-order ODEs. By implementing the scheme in Figure 3, one can easily neglect the lead–lag effect both on the rigid body and the rotor dynamics itself. Figure 9 reports the rigid body poles in the cases with and without blades’ lead–lag motion. While a small stabilizing effect is observed on the oscillatory modes (phugoid and dutch-roll), very small changes are observed in the frequency of these poles (Figure 10). The roll and short period are instead more affected by the presence of the lead–lag, and the frequency of these two stable and real dynamics moves from 5.8 to 10.7 rad/s and from 7.2 to 15 rad/s for, respectively, the roll and short period modes. The reason why this higher complexity affects the high-frequency poles is related to the coupling between the tip-path plane and the rotor center of gravity dynamics. The lead–lag angle is not directly responsible for an additional force/moment acting on the blades but affects the value of the flap. Table 2 reports the values of rotor disc tilt and CG coordinates with a pitch perturbation Δ q = 0.01 rad/s and the damping derivatives with and without lead–lag. It is observed that, while b 1 is similar in the two cases, a 1 significantly changes in the case of coupled lead–lag dynamics. The latter affects the absolute value of hinge stiffness and inertial moments computed by Equation (14) and increases the longitudinal damping M q when the lead–lag is “on”. The same process occurs with a roll rate perturbation p, and L p . In addition, the small variation in oscillatory poles can be linked to the variation in the coning angle, which is, in the case of a lagging rotor, larger.
The last and highest level of complexity can be achieved by adding a non-uniform inflow model, i.e., by considering, for both rotors, λ = [ λ 0 λ s λ c ] . Figure 11 represents the rigid body poles in the cases of uniform and non-uniform dynamic inflow. Three different x-axis scalings represent the low-, high-, and very-high-frequency poles in the complex plane. The very high-frequency plot shows the inflow poles derived from the 38-state linearization method. These six poles (three per rotor) are real and stable dynamics representing the uniform inflow at 33.5 rad/s and the lateral/longitudinal inflow modes at 126 rad/s ( 0.5 / r e v ). It is observed that the presence of a non-uniform inflow model affects mostly the short period, dutch-roll, and phugoid dynamics.
Considering the first one, an increase in the longitudinal damping is observed, linked to the mode coupling between the short period and longitudinal inflow λ c (see Figure 12 for the mode participation). A similar coupling effect which is highly detrimental to rotorcraft stability is the one between the longitudinal inflow and the phugoid mode. Even with a smaller mode participation, the longitudinal inflow contribution has the effect of moving the phugoid poles from a pure oscillatory, slightly unstable, situation, to a couple of real stable/unstable poles of similar frequency, as shown in Figure 10. On the other hand, the dutch-roll stability remains unchanged although its frequency decreases by about 50%. Coupling with both lateral and longitudinal inflows is observed in the dutch-roll and roll poles as well, although without significant changes in the lateral stability. It can be assumed that being the configuration that is symmetrical with respect to the longitudinal plane, the effects of the lateral inflows of the two rotors balance themselves, while the longitudinal contribution is detrimental to the low-frequency longitudinal stability.

3.2. Rotor Dynamics and Flap/Lead–Lag Coupling

The linearization methodology presented in this paper allows a 38-state space representation of the rotorcraft flight dynamics to be developed, including the rigid body, rotor, and inflow modes. The out-of-diagonal submatrices of the state matrix A (Equation (30)) describe the coupling between the three dynamics, while A B B , A R R , and A λ λ are the pure dynamic effects. In this section, the rotor–rotor modes are studied and represented for one single rotor. It is observed that no differences occur in hovering conditions between MR1 and MR2, and no particular coupling is present. Consider the complete numerical framework: complete, coupled, second-order flap and lead–lag dynamics and first-order non-uniform inflow description. A set of 12 oscillatory and stable poles is computed for each rotor, 6 related to flap and 6 related to lead–lag dynamics. Figure 13 shows the coupled collective, advancing, and regressive modes on the complex plane. The frequency and damping of the same eigenvalues are reported in Table 3. In general, the system is characterized by stable collective flap dynamics located at the angular frequency of the rotor, and a couple of advancing/regressive modes that have the same real part and a ± Ω on the imaginary axes with respect to the collective mode. A similar trend is observed on the lead–lag, which is a stable high-frequency dynamic, more dampened with respect to the flap. The high damping related to the lead–lag is connected to the spring damping coefficient ( K D ξ ) adopted in this model (see Table A2).
The identification of the rotor poles was achieved through the eigenvector analysis. In particular, the collective modes were characterized by the main contribution of the collective flap and lead–lag coefficients ( a 0 and ξ 0 ). On the other hand, advancing and regressive dynamics are mainly related to the rotor disc tilt and non-uniform lead–lag coefficients, i.e., a 1 , b 1 , ξ c , and ξ s . It can be observed from the plots of the mode participation (Figure 14) that a strong coupling was present between rotor dynamics, and the identification of the flap and lead–lag poles was not straightforward. Indeed, while the lead–lag always had a stronger contribution, the flap poles were identified by the frequency of the collective mode (close to Ω ) and its respective advancing/regressive poles. Coupling with rigid body modes (mainly p and q) is also present in the regressive flap.
Finally, the methodology for linearizing highly coupled numerical models highlighted the coupling between the flap, lead–lag, and inflow as well. In particular, Figure 15 shows the mode participation uniform inflow poles (also plotted in Figure 13). The latter are strongly coupled with the collective flap and lead–lag, as well as the rigid body dynamics. Indeed, a contribution from the rolling rate p is present in one of the two poles, while the other couples with the heave velocity w. This observation is coherent with the side-by-side configuration, where an inflow variation can create both a heave and roll movement.

4. Conclusions

The study presented in this paper highlighted the effects of a numerical modeling approach on side-by-side helicopter flight dynamics by introducing increased levels of complexity on relevant rotorcraft dynamic features (coning, flap, lead-l-ag, and inflow). To this aim, a 24 dof numerical model was obtained by integrating second-order flap and lead–lag dynamics, with a non-uniform, first-order, dynamic inflow and integral computation of single-blade forces and moments. The model was partially validated with a 14 dof analytical model on the trim and low-complexity features. The most outstanding results are summarized as follows:
  • A suitable trim methodology was developed. The algorithm is made of two nested loops and computes a trim solution averaged around one revolution. Trim curves at variable forward speed show a good agreement with the reference analytical results and provide typical helicopter behavior.
  • A methodology for linearizing highly coupled numerical models has been presented and validated by comparing the results with the analytical framework. The core issues that had to be addressed were (1) the loss of physical meaning in linearizing single-blade representations of rotor dynamics, (2) the high level of coupling between the flap and lead–lag, and (3) the presence of nonlinear complex terms derived from numerical integration. The guidelines to overcome these issues include (1) a non-rotating frame representation of rotor dynamics, (2) the partial decoupling of flap and lead–lag by approximating coupling and forcing terms with a previous time step solution, and (3) the introduction of additional stability derivatives linked to the coupling and excitation of rotor modes. Following these guidelines, a 38-state-space linear representation of rotorcraft dynamics has been developed. Rigid body, rotor, and inflow coupling effects are included in the representation.
  • As a result of the linearization method, the effects of single levels of complexity and rotor dynamics have been isolated and studied separately. It was observed that a dynamic disc tilt significantly increases the short period and roll damping and the lead–lag has a beneficial influence on the high-frequency modes. It was also observed that, by increasing the number of degrees of freedom in the flight dynamics framework, longitudinal instabilities arise, as observed from preliminary flight tests. Indeed, a coupled longitudinal inflow–phugoid mode led to dangerous instabilities.
  • Rotor dynamics have been addressed as well by showing the presence of coupled flap–lead–lag poles which behave as collective, advancing, and regressive stable modes. The collective flap is an oscillatory mode characterized by the angular frequency of the rotor, while the lead–lag is a higher-frequency dynamic denoted by a higher damping coefficient. The inflow has also been studied, showing coupling with the rotor collective dynamics and the rigid body heave and rolling motion.
In conclusion, the presented paper highlighted the capabilities of high-complexity numerical models in comparison with a simpler analytical representation of a side-by-side helicopter. A complete trim and linearization methodology was presented and flight dynamic characteristics were addressed.

Author Contributions

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

Funding

This study was carried out within the MOST-Sustainable Mobility National Research Center and received funding from the European Union Next-GenerationEU (PIANO NAZIONALE DI RIPRESA E RESILIENZA (PNRR)-MISSIONE 4 COMPONENTE 2, INVESTIMENTO 1.4-D.D. 1033 17/06/2022, CN00000023). This manuscript reflects only the authors’ views and opinions; neither the European Union nor the European Commission can be considered responsible for them. The publication of this study was funded by the Delft Technical University, through the Institutional Open Access Program (IOAP) with MDPI.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Acknowledgments

The authors acknowledge and thank SAB Group S.r.l. for providing the investigated rotorcraft prototype and sharing relevant data.

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A

Table A1. Side-by-side helicopter data.
Table A1. Side-by-side helicopter data.
DescriptionSymbolValue
Take-Off Mass [kg] m t o 20.62
Size [m]-2.5
Inertia moment with regard to x B [ kgm 2 ] I x x 3.532
Inertia moment with regard to y B [ kgm 2 ] I y y 2.222
Inertia moment with regard to z B [ kgm 2 ] I z z 5.342
Inertia product with regard to x B [ kgm 2 ] I y z 0
Inertia product with regard to y B [ kgm 2 ] I x z -0.052
Inertia product with regard to z B [ kgm 2 ] I x y -0.001
Table A2. Side-by-side helicopter’s rotor data (same design, counter-rotating).
Table A2. Side-by-side helicopter’s rotor data (same design, counter-rotating).
DescriptionSymbolValue
Sense of rotation Γ ± 1
Number of blades N b 3
Radius [m]R0.505
Mean chord [m]c0.051
Solidity ratio [-] σ 0.0964
Angular velocity [rpm] Ω 2400
Total hinge offset [m] e 0 0.075
Flap hinge offset [m] e F 0.0075
Root cutout from the flap hinge [m] r c 0.01
Blade mass [kg] m b 0.1613
Blade center of gravity with regard to the hub [m] r G 0.224
Spring restraint coefficient due to flap [Nm/rad] K s β 162
Spring restraint coefficient due to lag [Nm/rad] K s ξ 0
Spring damping coefficient due to flap [Nms/rad] K D β 0
Spring damping coefficient due to lag [Nms/rad] K D ξ 5
Pitch–Lag coupling ratio [-] K t ξ 0
Pitch–Flap coupling ratio [-] K t β 0
Blade twist coefficient [rad] θ w 0
Longitudinal incidence angle [rad] i s 0
Lateral incidence angle [rad] i c 0
Hub position of the MR1 with regard to body axes [m] r H B ( 1 ) [0 −0.645 0.066]
Hub position of the MR2 with regard to body axes [m] r H B ( 2 ) [0 0.645 0.066]

Appendix B

K 1 β = 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3 1 0 2 Ω sin ψ 1 2 Ω cos ψ 1 0 2 Ω sin ψ 2 2 Ω cos ψ 2 0 2 Ω sin ψ 3 2 Ω cos ψ 3
K 0 β = 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3 1 0 Ω 2 cos ψ 1 + Ω ˙ sin ψ 1 Ω 2 sin ψ 1 Ω ˙ cos ψ 1 0 Ω 2 cos ψ 2 + Ω ˙ sin ψ 2 Ω 2 sin ψ 2 Ω ˙ cos ψ 2 0 Ω 2 cos ψ 3 + Ω ˙ sin ψ 3 Ω 2 sin ψ 3 Ω ˙ cos ψ 3
For the sake of this paper, the same transformation is applied for both the flap and lead–lag, meaning that K 1 ξ = K 1 β and K 0 ξ = K 0 β .
C 1 β = 1 A 1 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3 1 D 1 β 0 0 0 D 2 β 0 0 0 D 3 β 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3
C 0 β = 1 A 1 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3 1 D 1 β 0 0 0 D 2 β 0 0 0 D 3 β 0 Ω sin ψ 1 Ω cos ψ 1 0 Ω sin ψ 2 Ω cos ψ 2 0 Ω sin ψ 3 Ω cos ψ 3
C 1 ξ = 1 A 2 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3 1 D 1 ξ 0 0 0 D 2 ξ 0 0 0 D 3 ξ 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3
C 0 ξ = 1 A 2 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3 1 D 1 ξ 0 0 0 D 2 ξ 0 0 0 D 3 ξ 0 Ω sin ψ 1 Ω cos ψ 1 0 Ω sin ψ 2 Ω cos ψ 2 0 Ω sin ψ 3 Ω cos ψ 3
E β = 1 A 1 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3 1 Q 1 β F 1 β Q 2 β F 2 β Q 3 β F 3 β
E ξ = 1 A 2 1 cos ψ 1 sin ψ 1 1 cos ψ 2 sin ψ 2 1 cos ψ 3 sin ψ 3 1 Q 1 ξ F 1 ξ Q 2 ξ F 2 ξ Q 3 ξ F 3 ξ

References

  1. Johnson, W. Helicopter Theory; Princeton University Press: Princeton, NJ, USA, 1980. [Google Scholar]
  2. Kunz, D.L. Comprehensive Rotorcraft Analysis: Past, Present, and Future. In Proceedings of the 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials Conference, AIAA 2005-2244, Austin, TX, USA, 18–21 April 2005. [Google Scholar]
  3. Chen, R.; Yuan, Y.; Thompson, D. A review of mathematical modelling techniques for advanced rotorcraft configurations. Prog. Aerosp. Sci. 2021, 120, 100681. [Google Scholar] [CrossRef]
  4. Talbot, P.D.; Tinling, B.E.; Decker, W.A.; Chen, R.T.N. A Mathematical Model of a Single Main Rotor Helicopter for Piloted Simulation. NTRS-NASA Technical Reports Server, 1982, NASA-TM-84281, Document ID: 19830001781. Available online: https://ntrs.nasa.gov/citations/19830001781 (accessed on 8 November 2024).
  5. Howlett, J.J. UH-60A Black Hawk Engineering Simulation Program. Volume 1: Mathematical Model. NTRS-NASA Technical Reports Server, 1981, NASA-CR-166309, Document ID: 19840020737. Available online: https://ntrs.nasa.gov/citations/19840020737 (accessed on 8 November 2024).
  6. Taamallah, S. A flight dynamics helicopter UAV model for a single pitch lag flap main rotor. In Proceedings of the 36th European Rotorcraft Forum, Paris, France, 7–9 September 2010. [Google Scholar]
  7. Bauchau, O.A.; Kang, N.K. A Multibody Formulation for Helicopter Structural Dynamic Analysis. J. Am. Helicopter Soc. 1993, 38, 3–14. [Google Scholar] [CrossRef]
  8. Savino, A.; Cocco, A.; Zanotti, A.; Tugnoli, M.; Masarati, P.; Muscarello, V. Coupling mid-fidelity aerodynamics and multibody dynamics for the aeroelastic analysis of rotary-wing vehicles. Energies 2021, 14, 6979. [Google Scholar] [CrossRef]
  9. Du, Val, R.; He, C. FLIGHTLABTM modeling for real-time simulation applications. Int. J. Model. Simul. Sci. Comput. 2017, 8, 1743003. [Google Scholar]
  10. Quon, E.; Smith, M.J.; Whitehouse, G.R.; Wachspress, D.A. Hierarchical Variable Fidelity Methods for Rotorcraft Aerodynamic Design and Analysis. In Proceedings of the American Helicopter Society 67th Annual Forum, Virginia Beach, VA, USA, 3–5 May 2011. [Google Scholar]
  11. Skogestad, S.; Postlethwaite, I. Multivariable Feedback Control: Analysis and Design; John Wiley & Sons: Hoboken, NJ, USA, 2005. [Google Scholar]
  12. Jun, D.; Cocco, A.; Saetti, U.; Juhasz, O. Flight Dynamics of a Coaxial Compound Helicopter with Rotor-On-Rotor Interactional Aerodynamics. In Proceedings of the 50th European Rotorcraft Forum, Marseille, France, 10–12 September 2024. [Google Scholar]
  13. Pavel, M.D.; den Ouden, W.; Bertolani, G.; Giulietti, F. Understanding the effects of rotor dynamics on helicopter incremental non-linear controllers. In Proceedings of the 48th European Rotorcraft Forum, Winterthur, Switzerland, 6–8 September 2022. [Google Scholar]
  14. Mazzeo, F.; Pavel, M.D.; Fattizzo, D.; Bertolani, G.; de Angelis, E.L.; Giulietti, F. Flight dynamic modeling and stability of a small-scale side-by-side helicopter for Urban Air Mobility. Aerosp. Sci. Technol. 2024, 148, 109117. [Google Scholar] [CrossRef]
  15. Padfield, G.D. Helicopter Flight Dynamics: The Theory and Application of Flying Qualities and Simulation Modelling; John Wiley & Sons: Hoboken, NJ, USA, 2008. [Google Scholar]
  16. de Angelis, E.L.; Giulietti, F.; Rossetti, G.; Turci, M.; Albertazzi, C. Toward Smart Air Mobility: Control System Design and Experimental Validation for an Unmanned Light Helicopter. Drones 2023, 7, 288. [Google Scholar] [CrossRef]
  17. Peters, D.A.; Haquang, N. Dynamic inflow for practical applications. J. Am. Helicopter Soc. 1988, 33, 64–68. [Google Scholar] [CrossRef]
  18. Bourtsev, B.N.; Selemenev, S.V. Fan–in–Fin Performance at Hover Computational Method. In Proceedings of the 26th European Rotorcraft Forum, The Hague, The Netherlands, 26–29 September 2000. [Google Scholar]
  19. Leishman, J.G. Principles of Helicopter Aerodynamics; Cambridge University Press: New York, NY, USA, 2000. [Google Scholar]
  20. Choi, H.S.; Kim, E.T.; You, D.I.; Shim, H. Improvements in small-scale helicopter rotor modeling for the real-time simulation of hovering flight. Trans. Jpn. Soc. Aeronaut. Space Sci. 2011, 54, 229–237. [Google Scholar] [CrossRef]
  21. Viterna, L.A.; Janetzke, D.C. Theoretical and Experimental Power from Large Horizontal-Axis Wind Turbines; NASA-TM-82944; NASA Lewis Research Center: Cleveland, OH, USA, 1982.
  22. Taamallah, S. Small-Scale Helicopter Blade Flap-Lag Equations of Motion for a Flybarless Pitch-Lag-Flap Main Rotor. In Proceedings of the AIAA Modeling and Simulation Technologies Conference, Portland, OR, USA, 8–11 August 2011. [Google Scholar]
  23. Chen, R.T.N. Effects of Primary Rotor Parameters on Flapping Dynamics. NTRS-NASA Technical Reports Server, 1980, NASA-TP-1431, Document ID: 19800006879. Available online: https://ntrs.nasa.gov/citations/19800006879 (accessed on 8 November 2024).
  24. Coleman, R.P.; Feingold, A.M. Theory of Self-Excited Mechanical Oscillations of Helicopter Rotors with Hinged Blades. NASA Langley Aeronautical Laboratory, 1957, NACA TN 3844, Document ID: 19930084623. Available online: https://ntrs.nasa.gov/api/citations/19930084623/downloads/19930084623.pdf (accessed on 8 November 2024).
  25. Butcher, J.C. Numerical Methods for Ordinary Differential Equations; John Wiley & Sons: Hoboken, NJ, USA, 2016. [Google Scholar]
  26. Diftler, M. UH-60A helicopter stability augmentation study. In Proceedings of the 14th European Rotorcraft Forum, Milan, Italy, 2–11 September 1988. [Google Scholar]
  27. Blackwell, J.; Feik, R.; Perrin, R. Identification of rotor dynamic effects in flight data. In Proceedings of the 15th European Rotorcraft Forum, Amsterdam, The Netherlands, 12–15 September 1989. [Google Scholar]
  28. Aponso, B.L.; Johnston, D.E.; Johnson, W.A.; Magdaleno, R.E. Identification of higher-order helicopter dynamics using linear modeling methods. J. Am. Helicopter Soc. 1994, 39, 3–11. [Google Scholar] [CrossRef]
Figure 1. Side-by-side helicopter (technical data in Appendix A). The rotorcraft has a take-off mass of 20.62 kg, a rotor radius of 0.5 m, and an overall size of 2.5 m.
Figure 1. Side-by-side helicopter (technical data in Appendix A). The rotorcraft has a take-off mass of 20.62 kg, a rotor radius of 0.5 m, and an overall size of 2.5 m.
Aerospace 11 00927 g001
Figure 2. Schematic representation of hub-body and blade frames of reference.
Figure 2. Schematic representation of hub-body and blade frames of reference.
Aerospace 11 00927 g002
Figure 3. Simulink algorithm to partially decouple and solve rotor dynamics in a non-rotating frame of reference.
Figure 3. Simulink algorithm to partially decouple and solve rotor dynamics in a non-rotating frame of reference.
Aerospace 11 00927 g003
Figure 4. Blade configurations for a clockwise rotor with N b = 3 and N c = 3 .
Figure 4. Blade configurations for a clockwise rotor with N b = 3 and N c = 3 .
Aerospace 11 00927 g004
Figure 5. Trim algorithm.
Figure 5. Trim algorithm.
Aerospace 11 00927 g005
Figure 6. Pilot controls (left) and rotorcraft attitude (right) in trim condition at variable forward speeds.
Figure 6. Pilot controls (left) and rotorcraft attitude (right) in trim condition at variable forward speeds.
Aerospace 11 00927 g006
Figure 7. Rigid body poles affected by lateral and longitudinal disc tilts, according to the numerical (black markers) and analytical (red markers) frameworks.
Figure 7. Rigid body poles affected by lateral and longitudinal disc tilts, according to the numerical (black markers) and analytical (red markers) frameworks.
Aerospace 11 00927 g007
Figure 8. Frequency ( ω ) of the rigid body poles affected by lateral and longitudinal disc tilts.
Figure 8. Frequency ( ω ) of the rigid body poles affected by lateral and longitudinal disc tilts.
Aerospace 11 00927 g008
Figure 9. Rigid body poles with and without the effect of the lead-lag dynamics. The high-frequency poles are represented on the left side of the plot, with a different x-axis scaling.
Figure 9. Rigid body poles with and without the effect of the lead-lag dynamics. The high-frequency poles are represented on the left side of the plot, with a different x-axis scaling.
Aerospace 11 00927 g009
Figure 10. Frequency ( ω ) of the rigid body poles with the inclusion of a complete flap, dynamic lead–lag, and non-uniform inflow complexities.
Figure 10. Frequency ( ω ) of the rigid body poles with the inclusion of a complete flap, dynamic lead–lag, and non-uniform inflow complexities.
Aerospace 11 00927 g010
Figure 11. Rigid body poles in the cases of uniform (black markers) and non-uniform dynamic inflow (green markers).
Figure 11. Rigid body poles in the cases of uniform (black markers) and non-uniform dynamic inflow (green markers).
Aerospace 11 00927 g011
Figure 12. Mode participation with non-uniform inflow.
Figure 12. Mode participation with non-uniform inflow.
Aerospace 11 00927 g012
Figure 13. Main rotor poles in the complex plane.
Figure 13. Main rotor poles in the complex plane.
Aerospace 11 00927 g013
Figure 14. Mode participation of the rotor-coupled poles.
Figure 14. Mode participation of the rotor-coupled poles.
Aerospace 11 00927 g014
Figure 15. Mode participation of the uniform inflow poles.
Figure 15. Mode participation of the uniform inflow poles.
Aerospace 11 00927 g015
Table 1. Rotor disc tilt with a perturbed pitch rate Δ q = 0.01 rad/s and lateral/longitudinal damping derivatives in the analytical and numerical models.
Table 1. Rotor disc tilt with a perturbed pitch rate Δ q = 0.01 rad/s and lateral/longitudinal damping derivatives in the analytical and numerical models.
AnalyticalNumerical
a 0 [deg] 0.55 0.53
a 1 [deg] 0.0041 0.0047
b 1 [deg] 0.0103 0.0144
M q [-]−3.60−7.22
L p [-]−3.57−5.83
Table 2. Effect of lead–lag dynamics on rotor flap and lead–lag coordinates with a perturbed pitch rate Δ q = 0.01 rad/s and damping derivatives.
Table 2. Effect of lead–lag dynamics on rotor flap and lead–lag coordinates with a perturbed pitch rate Δ q = 0.01 rad/s and damping derivatives.
Lead–Lag OFFLead–Lag ON
ξ 0 [deg]0 1.56
ξ c [deg]0 0.065
ξ s [deg]0 0.009
a 0 [deg] 0.53 0.69
a 1 [deg] 0.0047 0.023
b 1 [deg] 0.0144 0.0125
M q [-]−7.2−15
L p [-]−5.8−10.7
Table 3. Main rotor poles, frequency, and damping.
Table 3. Main rotor poles, frequency, and damping.
           PoleFrequency [rad/s]Damping [s]
Uniform Inflow (w coupling)−55.755.7-
Uniform Inflow (p coupling)−57.657.6-
Collective Flap 34 ± 265 i 2670.12
Regressive Flap 34 ± 15 i 370.9
Advancing Flap 34 ± 518 i 5190.07
Collective Lead–lag 223 ± 607 i 6470.34
Regressive Lead–lag 223 ± 356 i 4200.53
Advancing Lead–lag 223 ± 859 i 8870.25
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Mazzeo, F.; Pavel, M.D.; Fattizzo, D.; de Angelis, E.L.; Giulietti, F. Numerical Modeling, Trim, and Linearization of a Side-by-Side Helicopter in Hovering Conditions. Aerospace 2024, 11, 927. https://doi.org/10.3390/aerospace11110927

AMA Style

Mazzeo F, Pavel MD, Fattizzo D, de Angelis EL, Giulietti F. Numerical Modeling, Trim, and Linearization of a Side-by-Side Helicopter in Hovering Conditions. Aerospace. 2024; 11(11):927. https://doi.org/10.3390/aerospace11110927

Chicago/Turabian Style

Mazzeo, Francesco, Marilena D. Pavel, Daniele Fattizzo, Emanuele L. de Angelis, and Fabrizio Giulietti. 2024. "Numerical Modeling, Trim, and Linearization of a Side-by-Side Helicopter in Hovering Conditions" Aerospace 11, no. 11: 927. https://doi.org/10.3390/aerospace11110927

APA Style

Mazzeo, F., Pavel, M. D., Fattizzo, D., de Angelis, E. L., & Giulietti, F. (2024). Numerical Modeling, Trim, and Linearization of a Side-by-Side Helicopter in Hovering Conditions. Aerospace, 11(11), 927. https://doi.org/10.3390/aerospace11110927

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop