2.1 Multibody dynamics

The multibody system SIMPACK [Reference Lugner, Arnold and Vaculin13, 14] is used to set up and simulate the dynamic model of the rotor with its flexible rotor blades including large rigid body motions of the rotor hub and blade hinges as well as small deformations of the elastic structure. The development of the simulation package was originally initiated by DLR and later out-sourced for further development and commercial distribution.

SIMPACK is a general, non-linear three-dimensional multibody simulation tool and allows dynamic models for rotating and non-rotating systems consisting of rigid and elastic bodies to be built. The degrees-of-freedom are represented by a number of independent state variables [21] that define the motion of each body. In describing the kinematic behaviour, the motion or the position of the multibody system is studied with respect to the kinematic joints. Here, non-linear joint definitions provide the correct coupling of the bodies and introduce contributions from rotation w.r.t. stiffness and damping effects. For further modelling a huge element library for joints, force elements and other features is available. Also user force elements (UFEL) can be used to add dynamics by introducing ordinary differential equations (ODE) in terms of state space models, as done in this work for the indicial aerodynamics from section 2.2. Next to the simulation in time domain, the possibility of linearisation allows the system analysis by investigating the frequency and damping content. Here, the state space representation of the linear and time invariant system (LTI) is applied to derive the eigenbehaviour of the system.

If a dynamic multibody system is described by a minimum set of generalised coordinates q(t), the equations of motion form a set of non-linear second-order ODE [Reference Arnold, Burgermeister, Führer, Hippmann and Rill22] as given by

(1) \begin{align} M\!\left( q \right) \cdot \ddot q\!\left( t \right) = f\!\left( {q,\dot q\;,t} \right) \end{align}

with the mass matrix M(q) containing mass and inertia properties of all bodies as well as the vector of applied and gyroscopic forces f( $q,\dot q,\;t$ ) which depends on simulation time t. The equations of motion from Equation (1) can be also transferred into a generic state space form. If complex kinematics with constraints, connections in closed kinematic loops (see Fig. 2) or time-discrete state equations from discrete force and control elements need to be handled, the system motion is described by second-order differential-algebraic equations (DAE) with

(2) \begin{align} M\!\left( q \right) \cdot \ddot q\!\left( t \right) = f\!\left( {q,\dot q,\;t} \right) - {G^T}\!\left( {q,\;t} \right) \cdot \lambda \left( t \right) \end{align}

including the constraint forces -G ^T (q, t)·λ(t) that are determined by the constraint matrix G(q, t) and Lagrange multipliers λ(t).

Figure 2. 7AD rotor hub with articulated blade and closed kinematic loop to the swash plate.

Flexibility of the non-rotating blade is described in the multibody system SIMPACK either with a modal approach [Reference Schwertassek and Wallrapp23, Reference Wallrapp24] for available finite element models or on the basis of an intrinsic beam model. The following options allow the introduction of elastic rotor blades [14]:

1. (A) Implementation of complete elastic blade as one elastic body and additional geometric stiffness terms via standard FlexModal interface

2. (B) Implementation of elastic blade with connected elastic substructures via standard FlexModal interface

3. (C) Application of the intrinsic elastic beam model SIMBEAM

4. (D) Application of the rotor blade generator based on SIMBEAM

Option (A) enables the straight-forward use of finite element models of industrial model size. Currently, the FlexModal interface supports ABAQUS, ANSYS, MSC.NASTRAN, NX.NASTRAN, PERMAS and RADIOSS. The option (B) is advantageous to add further non-linear characteristics of the multibody joints that interconnect elastic substructures along the rotor blade. For beam models [Reference Surrey, Ortun, Truong and Wienke25], it can be easily combined with option (A). In addition, option (C) and (D) provide a solution, if a finite element code is not available or is not supported by FlexModal. In the present work, SIMPACK 2021 was used.

Option (A) which is based on the implementation of the complete elastic model together with additional geometric stiffness using the FlexModal interface is chosen for the 7AD rotor blades. Two FEM solutions from a preprocessing step with the finite element code MSC.NASTRAN are required in order to process the blade:

• • A modal solution provides the modal elastic model with natural frequencies and mode shapes for the non-rotating blade at Ω = 0 rad/s.

• • A static solution provides the geometric stiffness terms for the description in the relevant degrees-of-freedom of the rotating blade for Ω > 0 rad/s (centrifugal forces due to angular velocity).

FlexModal reads the model geometry, mass and stiffness matrices, natural frequencies, mode shapes, geometric stiffening terms and the used load case entries for the derivation of the geometric stiffening terms from the pre-processed finite element data. FlexModal then translates this input into a common flexibility description which is based on the standard input data (SID) format as described in Ref. [Reference Wallrapp24]. The file containing the flexible substructure is added in SIMPACK with marker locations which are the same as found at the node positions in the finite element model, changing a rigid into an elastic body. Regarding aspects of component modal synthesis and substructuring, any hybrid model consisting of rigid and flexible bodies can be built. Rigid bodies can be replaced by flexible ones by introducing a modal elastic model as described above or spring stiffness can be applied to a joint between two bodies from the MBS element library. SIMPACK adds the additional equations for the modal degrees-of-freedom to the set of DAE and solves the resulting equations.

2.2 Indicial aerodynamics

The unsteady aerodynamic model for rotor blade sections takes advantage of indicial response functions [Reference Leishman17], which capture the transient changes in aerodynamic variables within the Duhamel superposition integral. These functions model the unsteady aerodynamic response in time domain for variations in angle-of-attack and/or axial inflow. The model is based on Wagner’s function [Reference Leishman17] in combination with steady two-dimensional aerofoil data and applied to the multibody model. This function considers a two-dimensional thin aerofoil in impulsive motion from rest to uniform velocity U and is commonly used as a model for the circulatory lift development related to a change in angle-of-attack with fixed wing problems. According to Fung [Reference Fung26], the description of unsteady lift is extended from an impulsive to general motion including both pitch and heave motion by taking unsteady apparent mass terms for non-circulatory effects into account.

The fast aerodynamic model is based on radial independent strips and has been implemented using the approximation of R.T. Jones [Reference Arnold and Waitz18, Reference Arnold and Waitz19] for Wagner’s function Φ(τ) with

(3) \begin{align} \Phi(\tau) = 1 - 0.165{e^{ - 0.0455 \cdot \tau }} - 0.335{e^{ - 0.300 \cdot \tau }}\quad{\rm{where}}\,\tau = U \cdot t/b \end{align}

represents the non-dimensional reduced time. According to Equation (3), Wagner’s function describes the unsteady lift development towards the steady-state value for increasing distance τ traveled in semichords b during time t. The approximation is shown in Fig. 3 together with values from the state space formulation of Wagner’s function in the SIMPACK UFEL.

Figure 3. Wagner’s function for an incompressible fluid.

All indicial response functions affect the circulatory part of the unsteady lift, are implemented in terms of a state space model in the SIMPACK force element and extended for rotors under axial flow condition. Related contributions include blade aerofoil from the Wagner function [Reference Wagner27], axial dynamic inflow from a non-linear differential equation [Reference Leishman17] derived out of the results gained by Carpenter and Fridovich [Reference Carpenter and Fridovich28], and wake periodicity as modelled in Loewy’s function [Reference Loewy29] by a complex factor. The associated expressions for axial dynamic inflow and wake periodicity are listed as Equations (4) and (5), respectively. They introduce the rotor parameters of rotor angular velocity Ω, inflow λ _i , thrust coefficient C _T , pitch angular velocity ω and number of blades N _b . All circulatory contributions are accounted for by three dynamic states per radial strip, which can be configured in the model setup.

(4) \begin{align} \left( {\frac{{0.849}}{\Omega }} \right) \cdot {\dot \lambda _i} + 2 \cdot \lambda _i^2 = {C_T}\end{align}

(5) \begin{align} {e^{i \cdot 2\pi \cdot \omega /\left( {{N_b}\cdot\Omega } \right)}} = \cos \left( {\frac{{2\pi \cdot \omega }}{{{N_b}\cdot\Omega }}} \right) + i \cdot \sin \left( {\frac{{2\pi \cdot \omega }}{{{N_b}\cdot\Omega }}} \right)\end{align}

The fundamental solution of Equation (4) is found from integral tables for rational functions [Reference Bronstein and Semendjajew30] and leads to a kernel expression λ _i (t)=f(C _T , Φ _λ ) for dynamic inflow [Reference Arnold31] that contains the axial inflow function Φ _λ (τ _λ ) with

(6) \begin{align} \Phi _\lambda ( {{\tau _\lambda }}) = 1 - {e^{ - 3.3315 \cdot {\tau _\lambda }}}\quad{\rm{where}}\,{\tau _\lambda } = \sqrt {{C_T}} \cdot \Omega \cdot t\end{align}

represents the non-dimensional reduced time depending on rotor parameters. This function is a measure for the inflow change related to a pitch input and introduces a transition effect into the kernel expression with a similar curve shape as depicted in Fig. 3 for Wagner’s function, but ranging from zero to one. Since the mathematical setup of the axial inflow function resembles Equation (3), it can be iteratively solved together with the approximated Wagner function to adapt the motion induced lift from Wagner in terms of the thrust input for the inflow induced lift from the axial inflow function.

The influence of added rotor inflow on the development of unsteady lift is obtained by applying pitch step excitations in the SIMPACK force element and evaluating the temporal development for the angle-of-attack which is a measure for lift and can be compared to the pitch input. Here, the indicial response for the Wagner function and axial dynamic inflow interacts at each radial station r depending on two separate scales for reduced time, namely τ and τ _λ . The contribution of Wagner’s function with τ (rotor in-plane dynamics related to U = Ω $\cdot$ r) reaches the steady state in physical time earlier than the one of axial inflow with τ _λ (rotor out-of-plane dynamics related to $\surd$ C _T = $\surd$ 2·λ _i ). As illustrated in Fig. 4 for a two-bladed rotor, dynamic axial inflow reduces the effective angle-of-attack for the rotor section, whilst a wing section without inflow follows the basic Wagner function and reaches the commanded pitch angle. Hence, the onset of rotor inflow decreases the angle-of-attack towards a steady value which represents the so-called lift deficiency [Reference Miller32] of a rotor under axial inflow. To complete the extension of Wagner’s function to axial rotor flow, the periodicity term from Loewy’s function [Reference Loewy29] is implemented in the state space model, which was carefully tested for harmonic pitch oscillation to gain aerodynamic transfer functions for comparison of the basic Wagner function with the extended models. Figure 5 shows the lift amplitudes and phase angles versus rotor harmonics for the radial section r of a four-bladed rotor (Ω = 62.83 rad/s, r = 1.0m, b = 0.05m, C _T,0 = 0.008263) oscillating around mid chord. The influence of inflow and periodicity related to 4/rev can be seen in both lift amplitude and phase angle. As predicted by Loewy, reduced damping is found due to a smaller amplitude [Reference Ham1, Reference Johnson2] and a less negative angle for multiples of the blade passing frequency, i.e. 4/rev and 8/rev. For the used parameters, 1/rev equals a reduced frequency of k = 0.05 which is defined as:

(7) \begin{align} k = \omega \cdot \,b/U\quad {\rm{with}}\quad U = \Omega \cdot r\end{align}

Figure 4. Temporal development of angle-of-attack (AoA) for pitch step input due to axial inflow dynamics.

Figure 5. Unsteady lift amplitudes and phase angles of a four-bladed rotor obtained with Wagner function and extensions for rotor inflow ( $\texttt{+}$ ) and wake periodicity ( $\texttt{++}$ ).

2.3 Aeroelastic assessment

Standard features in MBS SIMPACK allow the aeroelastic assessment in a tightly coupled manner on equation level of the rotor blades for axial inflow condition in time domain and by using linearised results. In the latter approach which is based on linearisation, consecutive analysis of the eigenbehaviour and stability rating from damping values can be followed here, since the system matrices for the symmetric rotor with axial inflow remain time invariant and the unsteady aerodynamic forces are considered in the dynamic state set. The analysis is performed with the intrinsic eigensolver and delivers typical flutter results in terms of frequency and damping plots. In general, these results allow the verification of frequency placement and stability for the aeroelastic rotor system.

According to Ref. [Reference Arnold, Burgermeister, Führer, Hippmann and Rill22], the linearisation of the equations of motion from Equation (1) is the key to the linear stability and eigenmode analysis close to an equilibrium and to other methods of linear system analysis. Writing the right-hand side for the vector of applied and gyroscopic forces as f( $q,\dot q,\;u$ (t)) with system inputs u(t) the linearisation about an equilibrium given by $q\!\left( t \right) = {q^*},\;\;\dot q = 0,\;\;u = 0$ yields

(8) \begin{align} M\!\left( {{q^*}} \right) \cdot \ddot q = - {\rm{D}} \cdot \dot q - {\rm{K}} \cdot q + B \cdot u\left( t \right) \quad{\rm{with}}\quad D & = - \frac{{\partial f}}{{\partial \dot q}}\left( {{q^*},\;0,\;0} \right),\,K = - \frac{{\partial f}}{{\partial q}}\left( {{q^*},\;0,\;0} \right)\\ {\rm{and}} \quad B & = \frac{{\partial f}}{{\partial u}}\left( {{q^*},\;0,\;0} \right) \nonumber \end{align}

where the Jacobians $\partial f/\partial q$ , $\partial f/\partial \dot q$ and $\partial f/\partial u$ are approximated by classical finite differences.

3.1 Dynamic and aerodynamic models

The mechanical setup of the 7AD rotor system is depicted in Fig. 6 and consists of the flexible rotor blades and a rigid model for the DLR-ONERA rotor hub that is improved with measured values [Reference Hoffmann and Opitz33] for pitch link and lagging spring stiffness as well as lagging damping. As shown in Fig. 2, the rotor hub is modelled with SIMPACK features to close the kinematic loop composed of lead-lag hinge, flap hinge and pitch bearing via the pitch link to the swash plate. The flexibility of the non-rotating blade is described in the multibody system with a modal approach from a real modes solution and additional geometric stiffness contributions for the rotating blade are considered by static load cases for the relevant degrees of freedom, both computed in a preprocessing step with the finite element software MSC.NASTRAN. Furthermore, the centrifugal pitching moment is considered at each radial marker via a SIMPACK Expression. The four bladed rotor has a diameter of 4.2m, mean blade chord of 0.14m and a solidity of 0.085. Both the flap and lead-lag hinges are located at a radius of 0.075m and the blade attachment is found at the radial station of 0.275m. The blade mass is approximately 3.18kg. Figure 7 shows the MSC.NASTRAN beam model of the 7AD rotor blade along the straight blade axis and curved blade tip with additional rigidly connected nodes to model trailing and leading edge positions.

Figure 6. Multibody model of articulated 7AD rotor.

As described in a previous study [Reference Arnold4], measured vibration data is available from a ground vibration test (GVT) with rowing hammer technique for the clamped blade comprising modal properties in terms of natural frequencies and mode shapes for the non-rotating 7AD rotor blade. The comparison of simulation results for the clamped beam model with the experimental modal data takes the measured frequencies into account and global stiffness parameters of the finite element models are adjusted with respect to modulus of elasticity and shear, if necessary. The natural frequencies for the elastic blade modes obtained with the finite element and multibody model are given in Table 2 in terms of a relative frequency deviation with respect to the GVT results as reference values.

Table 2. Comparison of calculated frequencies with GVT data for clamped rotor blade (normalised to f_GVT)

Figure 7. Finite element model with rigidly connected leading and trailing edge.

Each rotor blade is equipped with 34 SIMPACK User Force Elements along the straight blade axis and the parabolic tip section. They provide the unsteady aerodynamic forces from circulation based on the extended Wagner function and related non-circulatory enhancements for the general motion of the aerofoil section taking heave and pitch into account. The unsteady aerodynamic contributions based on circulation include blade aerofoil, axial rotor inflow and wake periodicity from rotor blade interaction. Here, indicial functions allow to separate these effects into three cases according to Table 3 and to study their individual impact on rotor blade flutter. Similar to the previous aeroelastic assessment of the 7AD rotor [Reference Arnold4], a pitch setting according to a blade thrust of 827N at the nominal rotor speed of 100% is chosen for the test cases.

Table 3. Test cases for unsteady aerodynamic model

3.2 Dynamic simulation results

The dynamic simulation results for hover are gained from the application of linearisation features in MBS SIMPACK and show the development of the rotor frequencies versus nominal rotor speed for the dynamic rotor setup. As preparatory steps, an equilibrium computation for the derivation of the dynamic model states followed by a preload computation are performed to adjust the rotor model to the chosen rotor speed. The normalised frequencies of the 7AD rotor are found in Fig. 8 and correlate well with other results [Reference Surrey, Wendisch and Wienke34] available inside DLR. Here, the plotted modes are numbered in ascending frequency order and the colour code describes the major motion contribution in terms of blade lag (red), flap (blue) and torsion (green) at nominal rotor speed of 100%. Hence, the mode assigned as first elastic torsion T1 represents the flap bending mode F4 for rotor speeds below the nominal operation speed, and vice versa. The dashed lines show the rotor harmonics up to 8/rev and allow the verification of the frequency placement for the 7AD dynamic rotor setup.

Figure 8. Rotor frequencies for the dynamic model from linearisation (aerodynamics switched off).

3.3 Aeroelastic simulation results

The aeroelastic simulation results for rotor hover are shown in Figs 9–11 and comprise the development of the 7AD rotor frequencies and natural damping versus nominal rotor speed for the three studied aeroelastic setups. Here, the used definition of natural damping is the ratio of the negative real part and the absolute value of the eigenvalues. Again, the setups are processed using available SIMPACK linearisation features together with two preceding simulation steps to derive an equilibrium for all dynamic states including indicial aerodynamics and then to preload the rotor with the appropriate aerodynamic and centrifugal loads for the chosen rotor speed. As listed in Table 3, the flutter diagrams are extracted for three test cases which differ in the unsteady aerodynamic model for circulation comprising the basic Wagner function and two extensions for axial rotor inflow and wake periodicity. Here, the unsteady motion parts from the SIMPACK User Force Element are considered by aerodynamic lag states which are found as real parts in the linearisation result and, thus, influence the blade damping. As mentioned before, the dynamic rotor setup includes a measured value for lagging damping which introduces the large damping found with both lead-lag modes L1 and L2. Similar to an earlier assessment of the 7AD rotor [Reference Arnold4] utilising Wagner’s function and apparent mass contributions for unsteady aerodynamics, a pitch setting of around 3.70° at nominal rotor speed is chosen for the three cases. Due to the refined aerodynamic modelling which is especially found in the curved region of the blade tip, the new results show some differences. These comprise a first flap bending frequency F1 which is now slightly larger than 1/rev, inverted damping values for the lead-lag modes L1 and L2 as well as a coupling of L2 with the flap bending mode F3 for large rotor speeds.

Figure 9. Rotor frequencies and damping for the aeroelastic model from linearisation (Case 1/Wagner).

Figure 10. Rotor frequencies and damping for the aeroelastic model from linearisation (Case 2/Wagner $\texttt{+}$ ).

Figure 11. Rotor frequencies and damping for the aeroelastic model from linearisation (Case 3/Wagner $\texttt{++}$ ).

The simulation results for Case 1 based on Wagner’s function and non-circulatory apparent mass consider the unsteady aerodynamic contributions for the blade aerofoil and are depicted in Fig. 9. The flutter diagrams for the aeroelastic model comprise rotor frequencies and damping. In general, both curves allow the verification of frequency placement and stability for the aeroelastic rotor system. In comparison to the rotor frequencies from the dynamic simulation as presented in Fig. 8, all bending modes show larger frequencies except for mode F4. Here, the frequency slightly decreases at 100% of nominal rotor speed and indicates an aeroelastic coupling of bending mode F4 with torsion mode T1 for which the frequency increases. Above nominal rotor speed, the frequency of T1 is decreased what indicates now a coupling with bending mode F3 towards larger rotor speeds. In contrast, the frequency curves of both lead-lag modes remain unaltered. As already implied by the frequency behaviour of the aeroelastic model, the damping curves reveal a classical bending-torsion coupling of the third flap bending mode F3 and the first torsion mode T1 as known for articulated rotor blades. This mode pair becomes unstable at 151% rotor speed. The flutter coupling can be verified by switching off F3 and then, the negative damping of mode T1 disappears. Similar to the frequency curves of F4 and T1, the damping values for flap bending mode F4 represent those of torsion T1 below the nominal operation speed, and vice versa. In addition, a slight coupling is found between flap mode F3 and lead-lag mode L2 starting around the nominal rotor speed of 120%, which remains stable. Finally, the first flap mode F1 shows very large values of natural damping of around 30% and hence, is not depicted in the plot.

The simulation results for Case 2 utilising Wagner’s function, rotor inflow and non-circulatory apparent mass consider unsteady aerodynamic contributions for blade aerofoil and axial rotor inflow. The resulting flutter diagrams are found in Fig. 10, which are compared to the results of Case 1 that neglected the rotor inflow. Here, the rotor frequencies for the flap bending modes F2 and F4 show larger frequencies for Case 2, whilst the frequencies for F3, F5 and the torsion mode T1 decrease. As seen in the frequency increase of flap bending mode F4 and decrease of torsion mode T1 at 100% of nominal rotor speed, their aeroelastic coupling is less pronounced with axial inflow. Above the operation speed of 100%, the frequency of T1 is decreased which indicates again the coupling with bending mode F3 towards larger rotor speeds. All other modes comprising F1, L1 and L2 remain more or less unchanged in their frequency behaviour. As before, the damping curves reveal the bending-torsion coupling of the third flap bending mode F3 and the first torsion mode T1 with flutter onset at 158% of nominal rotor speed. As for the frequency curves of F4 and T1, the damping values for flap bending mode F4 also represent those of torsion T1 below the nominal operation speed, and vice versa. In general, the torsion mode T1 is more damped for all rotor speeds and elastic flap bending modes (F2–F5) are less damped than before, whilst damping remains rather the same for the rigid flap mode F1 as well as for both lead-lag modes L1 and L2. Again, a minor coupling of flap mode F3 and lead-lag mode L2 develops at the nominal rotor speed of 120%, which is stable within the studied speed range. Also, the flap mode F1 is not found in the plot, since it reaches large damping values around 30%. To summarise the comparison to Case 1 in terms of rotor blade flutter in hover, it can be stated for the assessment of the 7AD rotor that the introduction of axial rotor inflow with Case 2 decreases the damping of the elastic flap bending modes, but increases torsion damping. As a result, the flutter onset is shifted to a larger rotor speed.

The simulation results for Case 3 using Wagner’s function, rotor inflow, wake periodicity and non-circulatory apparent mass consider the unsteady aerodynamic contributions for blade aerofoil, axial rotor inflow and wake periodicity. The resulting flutter diagrams from Fig. 11 are compared to those of Case 2 that do not consider wake periodicity. The obtained frequency results are the same for all modes as found before with Case 2 except for the first torsion mode T1 which shows a very small increase between the nominal rotor speeds of 120% and 140%. Also, the computed damping results are the same as known from Case 2 except for a slight decrease for flap mode F2 at 60% of nominal rotor speed. As before, the rigid flap mode F1 shows large natural damping values around 30% and is not plotted. As a summary of the comparison to Case 2, it can be stated for the flutter assessment of the 7AD rotor in hover that the added wake periodicity within Case 3 does neither change the frequency nor the damping behaviour of the aeroelastic rotor system. Hence, the flutter onset at the nominal rotor speed of 158% is the same as found for Case 2.

Finally, flutter onset gained for the three test cases of the 7AD rotor in hover is compared in Table 4. The lowest flutter onset is found for unsteady aerodynamics limited to the blade aerofoil at a nominal rotor speed of 151% (Case 1/Wagner), whilst the two other cases with added rotor inflow and wake periodicity both show flutter at 158% (Case 2/Wagner $\texttt{+}$ , Case 3/Wagner $\texttt{++}$ ). Here, the influence of rotor inflow plays the major role, since it decreases the damping of the flap bending modes. As a consequence, the unstable torsion mode within the critical flutter coupling is increased in damping due to a weaker coupling to the participating flap mode and flutter onset is shifted to an approximately 5% larger rotor speed. For the studied 7AD rotor configuration in hover, added wake periodicity within Case 3 does neither change frequency nor damping and hence, does not affect the aeroelastic coupling in terms of a change in flutter speed.

Table 4. Flutter onset for unsteady aerodynamic models