2.1. Coordinate Systems

The coordinate systems shown in Figs. 1 and 2 are defined. First, we define the hub coordinate system \(F_h\), which is fixed to the fuselage and centered at the rotor hub. Each axis of \(F_h\) is parallel to the fuselage axes, with the positive \(x_h\)-axis pointing toward the nose, the positive \(y_h\)-axis pointing to starboard, and the positive \(z_h\)-axis pointing downward relative to the fuselage.

Subsequently, we define the blade coordinate system \(F_b\), which rotates with the blades. For the sake of simplifying the analysis, the flapping hinge is assumed to be at the same location as the rotor hub, with no hinge offset. The \(x_b\)-axis of \(F_b\) is in the spanwise direction of the blade (positive toward the tip), the \(y_b\)-axis is in the chordwise direction (positive in the direction of rotation), and the \(z_b\)-axis is chosen to form a right-handed system.

The main rotor is assumed to rotate at a constant angular velocity \(\Omega\) clockwise as viewed from above, and the azimuth angle \(\psi = \Omega t\) is defined with the nose direction as the reference. Additionally, the flap angle \(\beta\) is defined as the angle formed by the \(x_b\)-axis of the rotor blade with the \(x_h\)-\(y_h\) plane, with the upward (coning) direction taken as positive.

2.2. Equation of Motion for Flap Angle

The equation of motion for the flap angle of the rotor blade is derived from the balance of moments about the flapping hinge. Considering a minute blade element \(dr\) located at a distance \(r\) from the hub, we define the inertial force acting on it as \(dT_{in}\) and the aerodynamic force be \(dT_{ae}\) as shown in Fig. 2. The moment balance is expressed by the following equation:

Next, the aerodynamic force \(dT_{ae}\) acting on the element is determined. For the velocity components of the element shown in Fig. 3, the rotor rotation axis parallel component \(\mu_P\) and the perpendicular component \(\mu_T\), non-dimensionalized by the tip speed \(R\Omega\), are as follows:

Based on blade element theory, using the lift curve slope \(a\) and the chord length \(c\), \(dT_{ae}\) is approximated by the following equation:

2.3. Inflow Ratio Model

Under the high-speed flight conditions targeted in this study, it is essential to appropriately evaluate the air inflow condition (inflow ratio) to the rotor disk plane. Assuming steady level forward flight (\(\mu_y = 0, \mu_z = 0\)), Fig. 4 shows the relationship between the inclination angle \(i\) of the tip path plane (TPP) and the airspeed. Under this assumption, the non-dimensionalized airspeed \(\mu_\infty\) of the helicopter is simply equal to the forward speed component:

The advance ratio \(\mu\) and the inflow ratio \(\lambda\) are defined by the following equations as a superposition of the inflow component due to the airspeed of the helicopter and the induced velocity \(\nu\) associated with the rotation of the rotor.

In contrast, regarding the high-speed flight that is the main focus of this study, we assume that \(\mu_y\) is small, while \(\mu_\infty \sin i\) becomes non-negligible relative to \(\nu\). As shown in Figs. 5(a) and (b), with the increase in flight speed, the influence of the induced velocity \(\nu\) relatively decreases, and instead, the influence of the inflow component \(\mu_\infty \sin i\) due to the forward tilt of the fuselage becomes dominant. Therefore, as a model for the inflow ratio in the high-speed regime, a format that accounts for all terms is used, as shown in the following equation:

3.1. Targeted Helicopter and Rotor Parameters

This study targets the uncrewed helicopter “Fazer” with the flybar removed ⁴.

Its primary parameters are listed in Table 1. In this paper, all physical parameters and specifications are standardized to SI units, and angles are measured in radians.

3.2. Steady Forward Flight Conditions and Equation of Motion

This section examines how flight speed influences the characteristics of the flap angle \(\beta\). Note that the general governing equations prior to any simplification are provided in Appendix A. To simplify the discussion and derive the specific equation of motion for steady forward flight (Eq. (13)), we make the following assumptions:

1. ・

   The flight speed is constant, and the helicopter is flying forward. That is, \(\mu_x \neq 0\) and \(\mu_y = \mu_z = 0\).

2. ・

   The attitude of the helicopter is kept constant. That is, \(p = q = r = 0\).

3. ・

   \(\theta_0\), \(\theta_{1c}\), and \(\theta_{1s}\) are constant.

4. ・

   The blade parameters \(m_r, \beta, \theta, a, c\) are constant regardless of the spanwise position.

Under these assumptions, the equation of motion for the rotor flap angle is as follows:

The coefficients of \(\beta\) and \(\dot{\beta}\) on the left-hand side of Eq. (13) are functions of \(\psi\). This indicates that the equation of motion for the rotor flap angle constitutes a periodic system with time-varying coefficients. According to Floquet’s theory, the response of a periodic system subjected to periodic external forces (such as cyclic pitch control) contains not only the periodic components of the input but also higher-order harmonic components that are integer multiples of the input period. Especially under high-speed flight conditions where the contribution of higher-order components becomes significant, an approximation that limits the Fourier series expansion of \(\beta\) to only the zeroth and first orders, as in the following equation, is inadequate:

On the other hand, consider the case where flight speed is low, and the higher-order terms of \(\mu_x\), as well as the higher-order terms of the products of \(\mu_x\) with \(\theta\), \(\beta\), and their time derivatives, are sufficiently small. In this case, the equation of motion for the flap angle can be approximated as a constant-coefficient differential equation as follows:

The benchmark trim parameters and control inputs shown in Table 2 are uniquely determined from the steady-state flight equilibrium of the entire helicopter. Specifically, the trim state was calculated by constraining the vehicle velocities to the command values and the angular velocities to zero. The remaining unknown variables, including control inputs and attitude angles, were obtained by simultaneously and numerically solving the equations of motion (for the flap angle, vehicle translation, and vehicle rotation) based on the Fourier series expansion equations shown in Appendix A. This numerical solution was obtained using the “fsolve” command in Maple, which performs iterative calculations. To ensure a fair comparison of how the expansion order \(n\) affects flapping accuracy and aerodynamic calculations, this trim routine is executed using the conventional first-order flapping model (\(n = 1\)), incorporating fuselage aerodynamic drag. Once these benchmark control inputs (\(\theta_0, \theta_{1c}, \theta_{1s}\)) and attitude parameters are determined, they are kept completely constant across all simulations. For each target speed (\(\mu_x = 0.04\) and \(\mu_x = 0.25\)), we then evaluate the multi-harmonic flapping dynamics (Eq. (17)) for each order \(n\) under the exact same control benchmark. This isolated evaluation protocol allows us to evaluate the pure sensitivity of the aerodynamic force/moment estimations to the model order \(n\) under a consistent physical flight state.

3.3. Comparison of Numerical Calculation Results

Numerical solutions to the equation of motion (Eq. (13)) for the flap angle \(\beta\) under steady forward flight conditions are presented herein. We evaluate the steady-state response of \(\beta\) at both low-speed (\(\mu_x = 0.04\)) and high-speed (\(\mu_x = 0.25\)) conditions. Fig. 6 illustrates the variations in \(\beta\) as a function of the azimuth angle \(\psi\). Fig. 6 indicates that in low-speed flight, the periodic variation is dominated by the first-order component, with negligible higher-order harmonic content. These results validate the first-order expansion in Eq. (15) for low flight speeds. In contrast, in high-speed flight, \(\beta\) possesses harmonic components with respect to \(\psi\).

Figure 7 presents the power spectral density (PSD) of the flapping angle \(\beta\). Given the rotor angular velocity \(\Omega = 89\) rad/s, the fundamental rotor rotation frequency (1 /rev) corresponds to 14.2 Hz. In the figure, vertical dashed lines indicate the 1 /rev, 2 /rev, and 3 /rev frequencies for reference purposes. We find that in the low-speed regime (\(\mu_x = 0.04\)), the PSD exhibits spectral peaks only at the DC and 1 /rev frequencies. This result validates the treatment of Eq. (16) as a constant-coefficient differential equation and confirms the adequacy of the first-order expansion in Eq. (15) for low-speed conditions. In contrast, for the high-speed regime (\(\mu_x = 0.25\)), significant peaks emerge at the 2 /rev and 3 /rev frequencies in addition to the DC and 1 /rev components. These spectral characteristics demonstrate that as the advance ratio \(\mu_x\) increases, the influence of the azimuth-dependent periodic coefficients in Eq. (13) becomes non-negligible. This quantitatively proves that conventional first-order approximation models are insufficient for capturing rotor dynamics in high-speed flight regimes.

4.1. Higher-Order Harmonic Series Expansion of Flap Angle

To accurately capture the flap dynamics described by Eq. (13), we expand the flap angle \(\beta\) into a higher-order harmonic series. Letting \(n\) denote the highest order of the harmonic terms considered, the expansion is given by

We verify the modeling fidelity of the \(n\)-th order model by comparing its solution \(\beta_n(\psi)\) with the true value \(\beta(\psi)\) obtained through direct numerical integration of Eq. (13). Fig. 8 compares the flap angle responses for both low-speed and high-speed flight conditions. Additionally, we employ the mean absolute error defined by Eq. (18) as a quantitative index:

In contrast, under the high-speed flight condition (\(\mu_x = 0.25\)), expanding from a first-order to a higher-order model yields a substantial improvement in fidelity. Although accuracy further increases from \(n=2\) to \(n=3\), no significant gain is observed beyond the third-order. These results demonstrate that the conventional first-order approximation is inadequate for capturing the response characteristics in high-speed regimes. We conclude that achieving sufficient fidelity requires at least a second-order model, with a third-order or higher expansion being ideal for high-speed flight analysis.

4.2. Computational Implementation and Engineering Applications

The multi-harmonic constant-coefficient formulation presented in Appendix A offers two major computational advantages depending on the intended simulation or control scenario.

Steady-State Algebraic Solution:

For steady forward flight (where transient dynamics have decayed), the time-derivatives of the Fourier coefficients can be set to zero (\(\ddot{(\cdot)} = 0\), \(\dot{(\cdot)} = 0\)). This simplifies the differential equations in Appendix A to a set of static linear algebraic equations in the form of a matrix equation:

Transient State-Space Simulation:

When transient responses are needed, the constant-coefficient system can be easily converted into a standard linear time-invariant (LTI) state-space model. Unlike Eq. (13), which requires an extremely fine integration step size to resolve the high-frequency azimuth angle \(\psi = \Omega t\), this LTI model can be simulated using a significantly larger time step (e.g., \(10\) ms), matching the update rate of standard flight control computers.

Expected practical applications of this reduced-cost, high-fidelity model include:

1. ・

   Application of modern flight control (e.g., inverse dynamics transformation): Advanced control frameworks are based on state-space representations. Our constant-coefficient model provides this format directly, enabling high-speed control that accounts for higher-order flapping coupling and lateral-longitudinal cross-coupling without the use of computationally intensive time-varying models.

2. ・

   Real-time hardware-in-the-loop (HIL) simulation: In HIL systems running on embedded flight controllers, rotor flap dynamics can be executed at standard control update rates (100–400 Hz) rather than at kilohertz update rates.

3. ・

   Online rotor load/vibration monitoring: Estimation of dynamic blade loads and high-order fatigue stress (2 /rev and 3 /rev components) on board during flight, which can be linked to active vibration control, such as higher harmonic control.

4.3. Impact on Evaluation of Rotor Aerodynamic Characteristics

We express the aerodynamic force components acting on the rotor in the hub coordinate system (\(x_h, y_h, z_h\))—namely thrust \(T\), lateral force \(Y\), and drag \(H\)—as follows:

To evaluate the pure sensitivity of the aerodynamic force and moment estimations to the model order \(n\), we conduct the calculations using the identical benchmark trim parameters (vehicle attitude and control inputs) fixed in Section 3.2. We take the aerodynamic characteristics obtained by directly integrating the time-varying Eq. (13) as the reference values (true values) for comparison with those obtained from the constant-coefficient \(n\)-th order models. Based on the findings in the previous section, we limit the comparison to models of order up to the third order, as no significant improvement is expected for \(n > 3\). Fig. 10 presents the normalized estimation accuracy for aerodynamic forces and moments under the high-speed flight condition of \(\mu_x = 0.25\).

As shown in Fig. 10(a), the estimation accuracy of the aerodynamic forces significantly deteriorates with the first-order model, with substantial errors particularly in the lateral force \(Y\) and drag \(H\). This occurs because the asymmetry of the lift distribution between the advancing and retreating sides intensifies at higher forward speeds. Thus, the resulting higher-order flap angle components contribute significantly to the aerodynamic force balance. Since the first-order model neglects these components, it fails to represent this physical phenomenon. In contrast, the second- and third-order models show good agreement with the true values. It should be noted that while the normalized value for the lateral force \(Y\) at \(n=3\) (\(99.64\%\)) appears slightly less accurate than at \(n=2\) (\(100.12\%\)), this difference of less than \(0.5\%\) is not a physical degradation due to higher-order expansion. This minute discrepancy is a numerical artifact caused by the truncation of higher-order coupling terms and the limits of numerical integration (e.g., discretization and round-off errors) in Eq. (20). From an engineering perspective, both models are considered to have fully converged to the true value. We therefore conclude that a second-order model provides sufficient accuracy for describing the helicopter’s translational motion at high speeds.

Similarly, Fig. 10(b) confirms that the first-order model is inadequate for evaluating rotor aerodynamic moments. All components (\(M_{xh}\), \(M_{yh}\), \(M_{zh}\)) exhibit large errors because the first-order model fails to capture the azimuthal variation in the flap angle distribution. Conversely, adopting a second-order or higher model significantly improves estimation accuracy, with the second- and third-order models yielding nearly identical results.

These findings demonstrate that the conventional first-order model is insufficient for evaluating rotor aerodynamic characteristics in high-speed flight, underscoring the need for a second-order or higher-harmonic expansion. Furthermore, considering the trade-off between computational cost and accuracy, we conclude that a second-order model is the most practical and reasonable choice for high-speed flight analysis.