Nomenclature
-
$\textbf{A}$
,
$\textbf{H}$
-
aeroelastic system matrices
-
${\textbf{A}_1}$
,
${\textbf{B}_1}$
,
${\textbf{C}_1}$
,
${\textbf{D}_1}$
-
system matrices of aerodynamic subsystem
-
${\textbf{A}_2}$
,
${\textbf{B}_2}$
,
${\textbf{C}_2}$
,
${\textbf{D}_2}$
-
system matrices of structural dynamic subsystem
-
${\bar{{\textbf{A}}}_i}$
-
influence matrices in rational function approximation (
$i = 0,1, \cdots, n + 2$
) -
${A_i}$
,
${B_i}$
-
approximation constants for Wagner’s function (
$i = 1,2$
) -
$a$
-
location of midchord in front of elastic axis, nondimensionalised with respect to semichord
$b$
-
$b$
-
semichord length of aerofoil, m
-
${\bar b_i}$
-
coefficients for inflow states (
$i = 1, \cdots, N$
) -
$c$
-
chord length of aerofoil, m
-
${\textbf{D}_a}$
-
aerodynamic sensitivity with respect to aerodynamic states
-
${\textbf{F}_i}$
-
coefficient matrices of aerodynamic governing equation (
$i = 1,2,3,4$
) -
${f_1}$
,
${f_2}$
-
actuator forces, N
-
$\mathcal{L}$
,
$\mathcal{M}$
-
aerodynamic lift and moment on aerofoil, N/m and N
-
${\textbf{M}_a}$
,
${\textbf{C}_a}$
,
${\textbf{K}_a}$
-
aerodynamic inertial, damping and stiffness matrices
-
${\textbf{M}_s}$
,
${\textbf{C}_s}$
,
${\textbf{K}_s}$
-
structural inertial, damping and stiffness matrices
-
$N$
-
number of aerodynamic states
-
$n$
-
number of aerodynamic lag terms
-
$p$
-
roots of aeroelastic characteristic equation
-
${\textbf{Q}_i}$
-
matrices in first-order equation of motion with delays (
$i = 1, \cdots, 7$
) -
$\textbf{q}$
-
vector of aerofoil motion
-
$\textbf{R}$
-
aerodynamic load vector
-
${\textbf{R}_0}$
-
constant load vector
-
$s$
-
Laplace variable
-
$t$
-
time, s
-
${U_\infty }$
-
freestream velocity, m/s
-
${\textbf{u}_1}$
,
${\textbf{x}_1}$
,
${\textbf{y}_1}$
-
control input, state and output of aerodynamic subsystem
-
${\textbf{u}_2}$
,
${\textbf{x}_2}$
,
${\textbf{y}_2}$
-
control input, state and output of structural dynamic subsystem
-
$\bar{\textbf{W}}$
,
$\bar{\textbf{c}}$
-
coefficient matrices of inflow differential equation
-
$\bar w$
-
downwash, m/s
-
$\textbf{x}$
-
state of aeroelastic system
-
$Z$
-
variable for Z-transform
Greek symbol
-
$\alpha $
-
rigid-body pitching of aerofoil, rad
-
${\beta _i}$
-
constants for rational function approximation (
$i = 1, \cdots, n$
) -
$\boldsymbol\lambda $
-
aerodynamic state vector
-
$\xi $
-
rigid-body plunging of aerofoil, m
-
$\rho $
-
air density,
${\rm{kg}}/{{\rm{m}}^3}$
-
${\tau _a}$
-
actuation delay, s
-
${\tau _s}$
-
measurement delay, s
-
${\phi _w}$
-
Wagner’s function
1.0 Introduction
Abundant studies in aeronautical engineering are focused on wing and aircraft configurations under various flow and flight conditions. Among the studies, numerical simulations have been the most convenient approach to predict wing and aircraft responses. For example, beam-, plate- and shell-based finite-element (FE) models have been developed to represent the aircraft structures, coupling with an appropriate aerodynamic model, such as the simple potential flow-based formulations (e.g. Theodorsen aerodynamics) or the sophisticated high-fidelity computational fluid dynamics (CFD) models (e.g. NASA’s FUN3D). The accuracy of these numerical aeroelastic formulations generally relies on the fidelity of structural and aerodynamic models and the coupling scheme between them. Even though CFD techniques continue to mature, there are still disadvantages relevant to their solutions, such as errors introduced by simplified boundary conditions [Reference Samani and Hixon1] and costs in updating grids for dynamic problems [Reference Neves, Lawson, Bennett, Khanal and Hoff2]. On the other hand, modern FE formulations can model the structural inertial and rigidity properties with satisfactory accuracy. However, accurately capturing the damping, especially the nonlinear damping of composite structures, remains challenging. Nonlinear structural damping is essential to aeroelastic tailoring and control [Reference Eugeni, Saltari and Mastroddi3] and may also impact the post-flutter limit-cycle oscillations of an aircraft [Reference Chen, Zhang, Zhou, Wang and Mignolet4]. Although modern FE schemes can usually achieve a great convergence with mesh refinement, they cannot eliminate the modeling error introduced by structural dynamic formulations due to their inherent assumptions from mechanics (e.g. kinematic assumption) and uncertainties from material properties and manufacturing process [Reference Glaessgen, Reeder, Sleight, Wang, Raju and Harris5, Reference Raju, Murthy and Patel6].
Direct experimental measurement of wing and aircraft responses from wind-tunnel testing can address several above issues inherent to numerical simulations and offer valuable insights in supporting aircraft design. For example, with the structural component being tested, the uncertainties embedded in material properties and possible manufacturing and assembling flaws can be revealed [Reference Glaessgen, Reeder, Sleight, Wang, Raju and Harris5, Reference Raju, Murthy and Patel6], which is one of the above challenges for numerical simulations. However, wind-tunnel experiments can become inconvenient in some cases. For example, not all wind tunnels can accommodate the full-scaled specimen of modern large transport aircraft (e.g. Boeing 747 and Airbus 380). Additionally, the energy and power constraints of wind tunnels have limited the experimental study of full-size supersonic and hypersonic vehicles. Even though experiments can be carried out with sub-scaled models in wind tunnels, the design and build of sub-scaled models need to satisfy complex similitude relationships, such as Mach, Reynolds, Strouhal, Froude and Euler numbers. Bushnell [Reference Bushnell7] discussed that among the scaling parameters, Reynolds number shortfalls and discrepancies, caused by the model-scale factor modified by various tunnel pressurisation and cryogenic mitigation approaches, constitute a major and long-standing flight-to-wind tunnel scaling issue. Moreover, the aeroelastic deformations between full- and sub-scaled models must be appropriately scaled for aeroelastic studies, where matching both Reynolds and Froude numbers is difficult [Reference Katz and Plotkin8]. Wan and Cesnik [Reference Wan and Cesnik9] studied the geometrically nonlinear scaling of flexible vehicles, where the matching of Froude and Reynolds numbers was achieved by properly scaling the air (gas) density, yet the gas viscosity for the two models was not considered. Although scaling of gas viscosity could be potentially achieved by compressing the air [Reference Katz and Plotkin8] or applying heavy gas in wind-tunnel testing on sub-scaled models [Reference Anders, Anderson and Murthy10], several properties of heavy gas, including transition, separation and shock-boundary layer interaction, differ from air [Reference Anders, Anderson and Murthy10], leading to a mismatch between wind-tunnel testing and full-scale flight testing. In general, it is almost impossible to satisfy all the similitude relationships between full- and sub-scaled models. Therefore, sub-scale models cannot reproduce the consistent behaviour as the full-scaled models due to the mismatch in either structure or flow conditions. Bushnell [Reference Bushnell7] also discussed other implementation issues regarding wind-tunnel testing with sub-scaled models, including wind-tunnel wall effect, flow disturbance in wind-tunnels, model mounting influences, etc. These effects must be corrected for the wind-tunnel testing results before they can be applied to study free-flight vehicles [Reference Hackett and Cooper11–Reference Bertram, Hoffmann, Goertz, Gebbink and Janssen13].
To address the above issues of traditional numerical and experimental methods, it is imperative to develop new dynamic experimental approaches. Functional and self-contained hardware and specimens can be designed to capture real-world challenges in computational simulations. A promising approach in this direction is hardware-in-the-loop (HIL). HIL has been successfully applied to aeronautical engineering, including morphing wing optimisation [Reference Boria, Stanford, Bowman and Ifju14], testing of aircraft powerplant fuel cell [Reference Bradley, Moffitt, Mavris, Fuller and Parekh15], altitude control and navigation [Reference Strub, Theodoulis, Gassmann, Dobre and Basset16–Reference Waszniowski, Hanzálek and Doubrava18], flight safety of unmanned aerial vehicles [Reference Elston, Argrow, Frew, Houston and Straka19, Reference Goppert, Shull, Sathyamoorthy, Liu, Hwang and Aldridge20], just to name a few. HIL can also be used to verify the developed software for aerospace systems [Reference Patel, Patel, Deodhare and Chetty21, Reference Fritz, Winter, Freund, Pflueger, Zeile, Eickhoff and Roeser22]. In the above research studies, the data communication between software and hardware is often one-way and not always in real-time, dwarfing HIL’s value in studying the behaviour of highly dynamic and coupled systems, such as the aeroelastic characteristics of flexible wings.
Real-time hybrid simulation (RTHS) can be considered a type of HIL when the hardware (or specimen) and software (or computational simulation) interact in real time. It is a novel, powerful and cost-effective experimental technique for examining the behaviour of complex, full-scale structural systems under realistic loading conditions [Reference Mahin and Shing23–Reference Shing, Bursi and Vannan25]. Studying aircraft system dynamics using an RTHS platform is a new concept in aeronautical engineering. It couples computational and experimental components together to understand the structural, aerodynamic, aeroelastic and control characteristics of an aircraft wing member or the complete vehicle. Su and Song [Reference Su and Song26] proposed and investigated a RTHS platform to study the aeroelastic characteristics of flexible wings and vehicles. In [Reference Su and Song26], the part that was less understood or difficult to model (e.g. the nonlinear structural dynamics of a slender wing) was considered the physical experiment component; the part that can be numerically simulated with reasonable accuracy while inconvenient and expensive to capture in experiments (e.g. the aerodynamic loads or rigid body motion of a complete aircraft) was considered as the numerical component. Data communications of kinematic and load information were established by using actuators and sensors to interface these two components in real time and capture the aeroelastic behaviour of the original system. It is worth noting that the division of numerical and experimental components can be different for individual aeroelastic systems. For example, the aerodynamics of a 2-D aerofoil with a small vibration in incompressible and inviscid flow can be well captured by the Theodorsen aerodynamics. However, the aerodynamic loads of a finite wing in turbulent flow may be computationally expensive to solve, even numerically. While more details need to be addressed, as a proof of concept, the aeroelasticity of a 2-D aerofoil was studied with the proposed RTHS platform. The parametric study showed that the proposed RTHS platform could accurately capture the transient response of the aerofoil as long as the levels of process and measurement noises and actuation delay were within a bounded range.
However, the actuation and sensor delays not only contribute to the errors in the measured transient response of the RTHS system, but also bring unwanted dynamics into the coupled system to change the RTHS system stability characteristics, for example, destabilising the RTHS system to cause experimental failure [Reference Maghareh, Dyke, Prakash and Rhoads27–Reference Hayati and Song29]. Although compensation algorithms have been developed to mitigate the influence of actuator lags during real-time structural testing [Reference Hayati and Song29, Reference Chen and Ricles30], the actuator and sensor delays cannot be completely eliminated. For the RTHS platform of aircraft systems, such inherent delays may alter the stability characteristics obtained through the coupled RTHS system, such as the measured flutter boundary. While the prior study [Reference Su and Song26] demonstrated the great potential of applying the proposed RTHS platform to capture the transient response of flexible wings, the impact of actuation and sensor delays on the stability characteristics measured through RTHS has not been examined. In order to address this gap, this study will advance the prior RTHS study on flexible wings [Reference Su and Song26] by quantifying these delays’ influence on the flutter behaviour obtained via RTHS. Additionally, three potential flow-based aerodynamic formulations, suitable for the cases studied here, are considered in this study as the numerical component to understand how the RTHS system behaves with different aerodynamic models via comparative analysis.
Figure 1. A 2-D rigid aerofoil section.
2.0 Theoretical formulation
2.1 Aeroelastic equation of motion
As shown in Fig. 1, the aeroelastic equation of motion of a thin aerofoil section with plunging (
$\xi $
, positive down, measured from the position where the spring is balanced by the aerofoil’s weight) and pitching (
$\alpha $
, positive nose up) degrees-of-freedom is given by
\begin{align}{\textbf{M}_s}\ddot{\textbf{q}}\!\left( t \right) + {\textbf{C}_s}\dot{\textbf{q}}\!\left( t \right) + {\textbf{K}_s}\textbf{q}\!\left( t \right) = \textbf{R}\!\left( t \right) + {\textbf{R}_0}\end{align}
\begin{align}\dot{\boldsymbol{\lambda}} \!\left( t \right) = {\textbf{F}_1}\ddot{\textbf{q}}\!\left( t \right) + {\textbf{F}_2}\dot{\textbf{q}}\!\left( t \right) + {\textbf{F}_3}\textbf{q}\!\left( t \right) + {\textbf{F}_4}\boldsymbol\lambda \!\left( t \right),\end{align}
where
\begin{align} \textbf{q}(t) &= \{\xi(t) \quad \alpha(t)\}^T \nonumber\\[5pt] \textbf{R}(t) &= \{- \mathcal{L}(t)\quad \mathcal{M}(t)\}^T. \end{align}
In Equation (2),
$\mathcal{L}$
and
$\mathcal{M}$
are the aerodynamic lift and moment applied on the aerofoil section, respectively.
${\textbf{M}_s}$
,
${\textbf{C}_s}$
,
${\textbf{K}_s}$
are the structural mass, damping, and stiffness matrices of the aerofoil-spring system, respectively.
${\textbf{R}_0}$
is a constant load considering the gravity force and non-zero pitching angle of the aerofoil when the torsional spring is not stretched, independent from the degrees of freedom. Equation (1) can be written in the state-space form, given by
\begin{align}\dot{\textbf{x}}\!\left( t \right) = \textbf{Ax}\!\left( t \right) + \textbf{H},\end{align}
where the state variable is
\begin{align} \textbf{x}(t) = \begin{Bmatrix} \textbf{q}^T(t)\;\;\;\; & \dot{\textbf{q}}^T(t)\;\;\;\; & \boldsymbol{\lambda}^T(t) \\ \end{Bmatrix}^T .\end{align}
The details of Equation (1)–(4) can be found in Su and Song [Reference Su and Song26].
The aerodynamic governing equation in Equation (1) is general and can represent several linear unsteady aerodynamic formulations. Several aerodynamic formulations, including the finite-state inflow theory [Reference Peters, Johnson, Friedmann and Chang31–Reference Peters and He33], Wagner’s function [Reference Jones34], and rational functional approximation [Reference Eversman and Tewari35, Reference Wright and Cooper36], can be cast into the same form of the aerodynamic governing equation [Reference Hang, Fei and Su37]. However, the aerodynamic state
$\boldsymbol\lambda $
and coefficient matrices (
${\textbf{F}_1}$
,
${\textbf{F}_2}$
,
${\textbf{F}_3}$
, and
${\textbf{F}_4}$
) need to be defined differently in individual formulations as shown in the Appendix. Also demonstrated in the Appendix, the aerodynamic loads of these linear unsteady aerodynamic formulations are calculated using Theodorsen-like equations, given by
\begin{align}\textbf{R}\!\left( t \right) = {\textbf{M}_a}\ddot{\textbf{q}}\!\left( t \right) + {\textbf{C}_a}\dot{\textbf{q}}\!\left( t \right) + {\textbf{K}_a}\textbf{q}\!\left( t \right) + {\textbf{D}_a}\boldsymbol\lambda \!\left( t \right),\end{align}
with definitions of
${\textbf{M}_a}$
,
${\textbf{C}_a}$
,
${\textbf{K}_a}$
, and
${\textbf{D}_a}$
can be found in the Appendix.
2.2 Real-Time hybrid simulation platform
The proposed real-time hybrid aeroelastic simulation platform [Reference Su and Song26] is illustrated in Fig. 2. The aerofoil kinematic data are measured at each time step from the vibration experiment. The kinematic data are fed to the unsteady aerodynamic simulation in the computer. The simulation calculates the aerodynamic loads in real-time, based upon the kinematic input, which are transformed into the actuator forces (
${f_1}$
and
${f_2}$
) to actuate the aerofoil motion. The aerofoil response is continuously measured and sent to the computer for the simulation, which closes the loop of RTHS for the aeroelastic system.
Figure 2. Sketch of the proposed RTHS platform.
The RTHS platform was simulated numerically by modeling aerodynamic and structural dynamic subsystems individually, where the data communication between the two subsystems was enabled to allow for the coupling (see Fig. 3). One can refer to Su and Song [Reference Su and Song26] for details of the subsystems.
Figure 3. RTHS with actuator and sensor function blocks.
Figure 4. RTHS with actuator and sensor measurement delays.
2.3 RTHS system stability with actuator and sensor measurement delays
In the prior study [Reference Su and Song26], it is shown that RTHS can capture the transient response of flexible wings, but the impact of actuation and sensor delays on the stability characteristics of the coupled RTHS system has not been examined. In this study, the two coupled state-space subsystems shown in Fig. 3 are considered to study the impact of actuator and sensor measurement delays (
${\tau _a}$
and
${\tau _s}$
, respectively) on the system stability. Although the actuator and sensor dynamics may also impact the overall RTHS stability, the current study assumes such an impact is less significant than the delays and does not consider these dynamics. One can refer to Fig. 4 for the RTHS platform with both delays. The equation of motion of the coupled aeroelastic system with the delays is given by
\begin{align} \begin{gathered} \textbf{M}_s \ddot{\textbf{q}}(t) + \textbf{C}_s \dot{\textbf{q}}(t) + \textbf{K}_s \textbf{q}(t) = \textbf{R}(t - \tau_a) + \textbf{R}_0 \\[5pt] \dot{\boldsymbol{\lambda}}(t - \tau_a) = \textbf{F}_1 \ddot{\textbf{q}}(t - \tau_a - \tau_s) + \textbf{F}_2 \dot{\textbf{q}}(t - \tau_a - \tau_s) + \textbf{F}_3 \textbf{q}(t - \tau_a - \tau_s) + \textbf{F}_4 \boldsymbol{\lambda}(t - \tau_a), \end{gathered} \end{align}
where the aerodynamic load is
\begin{align}\textbf{R}\!\left( {t - {\tau _a}} \right) = {\textbf{M}_a}\ddot{\textbf{q}}\!\left( {t - {\tau _a} - {\tau _s}} \right) + {\textbf{C}_a}\dot{\textbf{q}}\!\left( {t - {\tau _a} - {\tau _s}} \right) + {\textbf{K}_a}\textbf{q}\!\left( {t - {\tau _a} - {\tau _s}} \right) + {\textbf{D}_a}\boldsymbol\lambda \!\left( {t - {\tau _a}} \right).\end{align}
Equation (6) is then rewritten as
\begin{align}{\textbf{Q}_1}\dot{\textbf{x}}\!\left( t \right) + {\textbf{Q}_2}\dot{\textbf{x}}\!\left( {t - {\tau _a}} \right) + {\textbf{Q}_3}\dot{\textbf{x}}\!\left( {t - {\tau _a} - {\tau _s}} \right) = {\textbf{Q}_4}\textbf{x}\!\left( t \right) + {\textbf{Q}_5}\textbf{x}\!\left( {t - {\tau _a}} \right) + {\textbf{Q}_6}\textbf{x}\!\left( {t - {\tau _a} - {\tau _s}} \right) + {\textbf{Q}_7},\end{align}
where
\begin{equation} \begin{aligned} \textbf{Q}_1 & = \begin{bmatrix} \textbf{I}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & \textbf{M}_s\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{0} \end{bmatrix} \,, & \textbf{Q}_2 & = \begin{bmatrix} \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{I} \end{bmatrix} \,, \\[5pt] \textbf{Q}_3 & = \begin{bmatrix} \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & - \textbf{M}_a\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & - \textbf{F}_1\;\;\;\; & \textbf{0} \end{bmatrix} \,, & \textbf{Q}_4 & = \begin{bmatrix} \textbf{0}\;\;\;\; & \textbf{I}\;\;\;\; & \textbf{0} \\[5pt] - \textbf{K}_s\;\;\;\; & - \textbf{C}_s\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{0} \end{bmatrix} \,, \\[5pt] \textbf{Q}_5 & = \begin{bmatrix} \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{0} \\[5pt] \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{D}_a \\[5pt] \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{F}_4 \end{bmatrix} \,, & \textbf{Q}_6 & = \begin{bmatrix} \textbf{0}\;\;\;\; & \textbf{0}\;\;\;\; & \textbf{0} \\[5pt] \textbf{K}_a\;\;\;\; & \textbf{C}_a\;\;\;\; & \textbf{0} \\[5pt] \textbf{F}_3\;\;\;\; & \textbf{F}_2\;\;\;\; & \textbf{0} \end{bmatrix} \,, & \textbf{Q}_7 & = \begin{Bmatrix} \textbf{0} \\[5pt] \textbf{R}_0 \\[5pt] \textbf{0} \end{Bmatrix} . \end{aligned}\end{equation}
In fact, Equation (8) reduces to Equation (3) if
${\tau _a}$
and
${\tau _s}$
are both zero, i.e.
\begin{equation} \begin{aligned} \textbf{A} & = \!\left(\textbf{Q}_1 + \textbf{Q}_2 + \textbf{Q}_3\right)^{-1} \!\left(\textbf{Q}_4 + \textbf{Q}_5 + \textbf{Q}_6\right) \\[5pt] \textbf{H} & = \!\left(\textbf{Q}_1 + \textbf{Q}_2 + \textbf{Q}_3\right)^{-1} \textbf{Q}_7. \end{aligned}\end{equation}
Nonetheless, with
${\textbf{Q}_7}$
omitted for the stability study, the characteristic equation is
\begin{align}{\rm{det}}\!\left[ {\!\left( {p{\textbf{Q}_1} - {\textbf{Q}_4}} \right) + \!\left( {p{\textbf{Q}_2} - {\textbf{Q}_5}} \right){e^{ - p{\tau _a}}} + \!\left( {p{\textbf{Q}_3} - {\textbf{Q}_6}} \right){e^{ - p\!\left( {{\tau _a} + {\tau _s}} \right)}}} \right] = 0.\end{align}
which can be numerically solved to find its roots
$p$
(eigenvalues) and the associated eigenvectors. Finally, one needs to study the aeroelastic system’s stability with the change of freestream speed
${U_\infty }$
and delays (
${\tau _a}$
and
${\tau _s}$
). This is done by tracking the real part of the aeroelastic root
$p = p\!\left( {{U_\infty },{\tau _a},{\tau _s}} \right)$
. The occurrence when the real part of a root turns from negative to positive indicates the onset of instability of the hybrid aeroelastic simulation system.
3.0 Numerical studies
The numerical studies are carried out with a 2-D thin aerofoil (refer to Fig. 1), whose properties are listed in Su and Song [Reference Su and Song26]. For simplicity, no structural damping is considered in this work.
The finite-state inflow theory (FSI), Wagner’s function (Wagner), and Roger’s rational function approximation (RFA) are used to calculate the aerodynamic loads on the aerofoil in the real-time hybrid simulation. There are six aerodynamic states (inflow states) involved in the finite-state inflow theory. The Wagner’s function inherently has two aerodynamic states. Four lag terms are used for the RFA using Roger’s method, resulting in eight aerodynamic states. The reduced frequency range considered for the RFA is from 0 to 3.
3.1 Stability of nominal aeroelastic system
The aeroelastic stability characteristic of the aerofoil can be conveniently identified by solving the eigenvalue problem of the system matrix
$\textbf{A}$
in Equation (10). A sequence of
$\textbf{A}$
matrices can be formed in a range of freestream speed
${U_\infty }$
. The aerofoil’s flutter boundary is found when a root locus crosses the imaginary axis, as shown in Fig. 5. Without comparing the aerodynamic modes, the root locus of the aeroelastic modes obtained from the three models are similar to each other, and so is the flutter boundary, which is found to be 35.6 (FSI), 35.7 (Wagner) and 35.9 (RFA) m/s, respectively.
Figure 5. Nominal aeroelastic root locus without delay, with Wagner’s function,
${{\rm{U}}_\infty }$
from 0 (triangle) to 40 m/s (square).
Figure 6. Individual frequency and damping of RTHS system roots with equal delays, aerodynamic loads from Wagner’s function.
Figure 7. Flutter speed and frequency of RTHS with equal delays.
Note that the aerodynamic modes will be removed from the comparisons in the following studies.
3.2 Equal actuation and sensor measurement delays
Due to the complexity of the study that involves both the actuation and sensor measurement delays, one may assume the two delays are equal (i.e.
${\tau _a} = {\tau _s} = \tau $
). The RTHS system’s stability is studied with the change of freestream speed
${U_\infty }$
and delay
$\tau $
. Figure 6 plots the frequency and damping of the RTHS system’s roots vs. the freestream speed with selected delays (0, 0.5, 2, 5 and 15 ms, respectively), using the Wagner’s function. As can be seen, when
$\tau $
is small (e.g. 0.5 or 2 ms), the RTHS system does show a trend of frequency coalescence, just like the nominal aeroelastic system with no delays. However, this is not the case if
$\tau $
is larger than 2 ms (e.g. 5 or 15 ms), where the frequency coalescence effect is not apparent. While not all results are shown here, the RTHS systems with delays are reasonably consistent when modeled using the three aerodynamic models. The consistency can be observed from the flutter (instability) boundary predicted by the RTHS systems with delays shown in Fig. 7. In Fig. 7a, the straight lines define a
$ \pm 5{\rm{\% }}$
tolerance band about the flutter speed of the nominal aeroelastic system with no delays, where the aerodynamic loads are calculated using the finite-state inflow theory as a reference. With delays, the three aerodynamic formulations still agree well, with the Wagner’s function being the most conservative. If both the actuation and sensor measurement delays are limited to less than 4.0 ms, the instability boundary of the RTHS system has a relative error of less than
$5{\rm{\% }}$
, compared to the nominal aeroelastic system. This is within an acceptable range. For example, Tang and Dowell [Reference Tang and Dowell38] measured the flutter boundary of a flexible wing with a 3 to 5% relative error compared to the theoretical value. Obviously, the actuation and sensor measurement delays in RTHS can be tightened up if more measurement accuracy is needed. Figure 8 exemplifies how the root loci of the aeroelastic modes change with the delay, where the aerodynamic modes are no longer shown for the comparison. Essentially, the root loci shown in Fig. 8b with the zero delays are identical to those of the two aeroelastic modes shown in Fig. 5, since they are both the solution of the nominal aeroelastic system using the same aerodynamic formulation.
Figure 8. Root loci of RTHS system modes with finite-state inflow theory or Wagner’s function (colours are used to distinguish the root loci with different delays).
Figure 9. Flutter speed and frequency of RTHS with unsteady aerodynamics from finite-state inflow theory.
From Fig. 7b, one can see the jump of the flutter frequency of the RTHS system associated with the finite-state inflow theory, which indicates a discontinuity of mode shape with the change of delay. This phenomenon is further exemplified in Fig. 8a by the shift of locus between 0.5 and 0.6 ms. On the other hand, the model with Wagner’s function (as well as RFA, but not shown) does not show the locus shift (see Fig. 8b).
3.3 Individual actuation and sensor measurement delays
In this case, the actuation and sensor measurement delays are considered independent from each other, both varying from 0 to 15 ms. Figure 9 provides the surface plot of flutter boundary and flutter frequencies of the RTHS system, using the FSI only. A flat plateau of the flutter boundary surface can be observed when the actuation and sensor measurement delays are less than 4 ms. In contrast, the flutter boundary decreases linearly with the increase in delay. Figure 10 plots the contour of the flutter boundary surfaces from the three aerodynamic models, where the values are ratios of the flutter speeds of the three RTHS cases compared to the nominal aeroelastic system with the finite-state inflow theory. From the plot, one can see the sensitivity of the flutter boundary versus the individual delays. Again, when the delays are less than 4 ms, the flutter boundary is insensitive to the delays. When the delays are above it, the sensitivity to both delays becomes a constant value. One can also observe that the flutter boundary has a very similar sensitivity to the individual delays. Figure 10 clearly shows the required actuation and sensor measurement delays to meet the accuracy of the stability boundary of the RTHS system. Finally, Fig. 11 compares the flutter boundary predicted by the different aerodynamic models in the studied ranges of the actuation and sensor measurement delays where the maximum relative error among the different models is less than
$2.5{\rm{\% }}$
.
Figure 10. Contour of flutter speed surface plots (blue dash-dot: finite-state inflow theory; red solid: Wagner’s function; black dash: RFA).
Figure 11. Comparison of flutter speed of RTHS with different aerodynamic models.
3.4 Consistency verification
The mode shape of the RTHS aeroelastic system may change with the delays. It is important to understand how the delays may impact the unstable mode, especially around the flutter boundary. Let the flutter mode of the nominal aeroelastic system be
${\bar{\textbf{x}}_0}$
, the similarity between the flutter mode with delays and
${\bar{\textbf{x}}_0}$
can be evaluated by
\begin{align}MAC\!\left( {{\tau _a},{\tau _s}} \right) = \frac{{\langle \bar{\textbf{x}}\!\left( {{\tau _a},{\tau _s}} \right),{{\bar{\textbf{x}}}_0}\rangle \cdot \langle {{\bar{\textbf{x}}}_0},\bar{\textbf{x}}\!\left( {{\tau _a},{\tau _s}} \right)\rangle }}{{\langle \bar{\textbf{x}}\!\left( {{\tau _a},{\tau _s}} \right),\bar{\textbf{x}}\!\left( {{\tau _a},{\tau _s}} \right)\rangle \cdot \langle {{\bar{\textbf{x}}}_0},{{\bar{\textbf{x}}}_0}\rangle }},\end{align}
where
$\bar{\textbf{x}}\!\left( {{\tau _a},{\tau _s}} \right)$
are the flutter mode with actuation and sensor measurement delays. Operator
$\langle \cdot, \cdot \rangle $
is the inner product of two vectors, which is not commutative for the complex mode shapes. Figure 12 provides a surface plot and its contour of the MAC number of the flutter modes, using the Wagner’s function for aerodynamic modeling. Note that only the structural mode components (i.e. pitching, plunging and their rates) are involved in the comparison. One may keep the MAC number above 0.99 if both delays are less than 2.3 ms, as shown by the shaded region in Fig. 12b.
Figure 12. MAC number of flutter modes of RTHS system.
Figure 13. Simulation of RTHS system with actuation and sensor measurement delays.
As a verification of the stability analysis formulation and results of RTHS systems, a time-domain simulation of the RTHS with both actuation and sensor measurement delays is carried out. The RTHS system can be referred to in Fig. 4. The aerofoil is placed in a freestream of 34 m/s, with zero initial plunging and pitching displacements. For simplicity, the actuation and sensor measurement delays are considered equal. From Fig. 13, when the delay is 3.5 ms, the aeroelastic response of the aerofoil converges, while the response diverges if the delay is 4.5 ms. Therefore, the RTHS system turns to be unstable when the delay increases from 3.5 ms to 4.5 ms. This behaviour agrees well with the flutter boundary predicted and plotted in Fig. 7, where the flutter boundary is about 34 m/s with a delay of around 4 ms. Figure 13 also indicates that the RTHS system with a delay less than 4 ms can maintain a similar stability property as the nominal aeroelastic system. If one wants to measure more accurate transient responses using the RTHS system, the delays must be further reduced.
4.0 Conclusion
The stability of RTHS systems for the aeroelasticity of flexible wings was studied in this paper. In the proposed RTHS platform, the aerodynamic loads on the flexible wings (represented by a typical 2-D aerofoil) were predicted by three well-known 2-D unsteady aerodynamic formulations. Actuation and sensor measurement are essential in RTHS systems, and the delays developed in their operations inherently impact the measured stability characteristics of the aeroelastic system. In this study, actuation and sensor measurement delays were modeled in the RTHS platform, where the aeroelastic roots of the RTHS system were solved and evaluated for stability. By the parametric study performed in this paper, it has been found that the delays can alter the stability characteristics of the RTHS system, changing the measured aeroelastic stability boundary and unstable modes with different delay amounts. However, the RTHS platform can still accurately capture the nominal system’s stability characteristics when both delays were small (e.g. less than 4.0 ms). This finding will guide the implementation of aeroelastic RTHS platforms by providing important specifications on the future delay compensation design. Moreover, a consistent RTHS system’s stability boundary was obtained when different aerodynamic models were applied as the numerical component, providing cross-verification among these models. Finally, the RTHS system’s stability boundary with delays was also verified using time-domain simulations. In conclusion, the results obtained in this study have quantified the influence of actuation and sensor measurement delays on RTHS system stability and offer important guidance on the future implementation of the proposed RTHS platform.
Competing interests
The authors declare none.
Appendix
The aerodynamic subsystem can be represented by the unified Equation (1b) when the three aforementioned aerodynamic formulations are applied individually. However, the definition and physical meaning the aerodynamic states
$\boldsymbol\lambda $
are different in different aerodynamic formulations. In the finite-state inflow theory, the weighted summation of the aerodynamic states accounts for the additional induced velocity due to the wake, i.e.
\begin{align}{\lambda _0} = \frac{1}{2}\mathop \sum \limits_{i = 1}^N {\bar b_i}{\lambda _i} = \frac{1}{2}\bar{\textbf{b}}\boldsymbol\lambda, \end{align}
where
$N$
is the number of inflow (aerodynamic) states defined on the aerofoil, and
${\bar b_i}$
are coefficients that can be obtained by the least-square method [Reference Peters, Karunamoorthy and Cao32]. With the Wagner’s function, the two aerodynamic states of the aerofoil are
\begin{equation} \begin{aligned} \lambda_1(t) & = A_1 B_1 \frac{U_\infty}{b} \int_0^t \bar w_{\frac{3}{4}c}(\sigma) e^{- B_1 \frac{U_\infty}{b} \!\left(t - \sigma\right)} \textrm{d}\sigma \\[5pt] \lambda_2(t) & = A_2 B_2 \frac{U_\infty}{b} \int_0^t \bar w_{\frac{3}{4}c}(\sigma) e^{- B_2 \frac{U_\infty}{b} \!\left(t - \sigma\right)} \textrm{d}\sigma \,, \end{aligned}\end{equation}
where
${\bar w_{\frac{3}{4}c}}$
is the downwash at three quarters from the leading edge.
${A_1}$
,
${A_2}$
,
${B_1}$
and
${B_2}$
are constants to approximate the Wagner’s function. The summation of these two states is the integral of convolution, given by
\begin{align}{\lambda _1} + {\lambda _2} = - \mathop \int \nolimits_0^t {\bar w_{\frac{3}{4}c}}\!\left( \sigma \right)\frac{{d{\phi _w}\!\left( {t - \sigma } \right)}}{{d\sigma }}d\sigma, \end{align}
where
${\phi _w}$
is the Wagner’s function. With the RFA, the aerodynamic states are lagged terms from the aerofoil kinematics, given by
\begin{align}{\lambda _i}\!\left( s \right) = \dfrac{s}{{s + \!\left( {\frac{{{U_\infty }}}{b}} \right){\beta _i}}}\textbf{q}\!\left( s \right),{\rm{\;\;\;\;\;\;\;\;}}\!\left( {i = 1,2, \cdots, n} \right),\end{align}
where
$s$
is the Laplace variable.
$n$
is the number of lag terms.
${\beta _i}$
are pre-determined positive real numbers. If one follows the Roger’s method, the number of aerodynamic states is
$N = 2n$
.
Since the physical meanings of the aerodynamic states are different, the corresponding coefficient matrices also vary. For the finite-state inflow theory, the coefficient matrices in Equation (1b) are
\begin{equation} \begin{aligned} \textbf{F}_1 & = \bar{\textbf{W}}^{-1} \bar{\textbf{c}} \begin{Bmatrix} 1 & b \!\left(\dfrac{1}{2} - a\right) \end{Bmatrix} \\[5pt] \textbf{F}_2 & = \bar{\textbf{W}}^{-1} \bar{\textbf{c}} \begin{Bmatrix} 0 & U_\infty \end{Bmatrix} \\[5pt] \textbf{F}_3 & = \textbf{0}_{N \times 2} \\[5pt] \textbf{F}_4 & = \bar{\textbf{W}}^{-1} \!\left(- \frac{U_\infty}{b} \textbf{I}_{N \times N}\right), \end{aligned}\end{equation}
where the coefficients
$\bar{\textbf{W}}$
and
$\bar{\textbf{c}}$
are both defined in Peters et al. [Reference Peters, Karunamoorthy and Cao32]. For the Wagner’s function, the coefficient matrices are
\begin{equation} \begin{aligned} \textbf{F}_1 & = \begin{bmatrix} 0\;\;\;\; & 0 \\[5pt] 0\;\;\;\; & 0 \end{bmatrix} \\[5pt] \textbf{F}_2 & = \begin{bmatrix} A_1 B_1 \dfrac{U_\infty}{b}\;\;\;\; & A_1 B_1 U_\infty \!\left(\dfrac{1}{2} - a\right) \\[5pt] A_2 B_2 \dfrac{U_\infty}{b}\;\;\;\; & A_2 B_2 U_\infty \!\left(\dfrac{1}{2} - a\right) \end{bmatrix} \\[5pt] \textbf{F}_3 & = \begin{bmatrix} 0\;\;\;\; & A_1 B_1 \dfrac{U_\infty^2}{b} \\[5pt] 0\;\;\;\; & A_2 B_2 \dfrac{U_\infty^2}{b} \end{bmatrix} \\[5pt] \textbf{F}_4 & = \begin{bmatrix} - B_1 \dfrac{U_\infty}{b}\;\;\;\; & 0 \\[5pt] 0\;\;\;\; & - B_2 \dfrac{U_\infty}{b} \end{bmatrix} . \end{aligned}\end{equation}
With the RFA, the matrices take the form of
\begin{equation} \begin{aligned} \textbf{F}_1 & = \textbf{0}_{2n \times 2} \\[5pt] \textbf{F}_2 & = \begin{bmatrix} \textbf{I}_{2 \times 2}, \cdots, \textbf{I}_{2 \times 2} \end{bmatrix}^T \\[5pt] \textbf{F}_3 & = \textbf{0}_{2n \times 2} \\[5pt] \textbf{F}_4 & = - \frac{U_\infty}{b} \begin{bmatrix} \beta _1 \textbf{I}_{2 \times 2}\;\;\;\; & \;\;\;\; & \\[5pt] \;\;\;\; & \ddots\;\;\;\; & \\[5pt] \;\;\;\; &\;\;\;\; & \beta _n \textbf{I}_{2 \times 2} \end{bmatrix}. \end{aligned}\end{equation}
The aerodynamic load is the output of the aerodynamic subsystem, as shown in Equation (5). For the finite-state inflow theory, the aerodynamic output matrices are
\begin{equation} \begin{aligned} \textbf{M}_a & = \pi \rho b^2 \begin{bmatrix} - 1\;\;\;\; & ba \\[5pt] ba\;\;\;\; & - b^2 \!\left( \dfrac{1}{8} + a^2 \right) \end{bmatrix} \\[5pt] \textbf{C}_a & = 2 \pi \rho b U_\infty \begin{bmatrix} - 1\;\;\;\; & - b \!\left( 1 - a \right) \\[5pt] b \!\left( \dfrac{1}{2} + a \right)\;\;\;\; & b^2 \!\left( \dfrac{1}{2} - a \right) a \end{bmatrix} \\[5pt] \textbf{K}_a & = 2 \pi \rho b U_\infty^2 \begin{bmatrix} 0\;\;\;\; & - 1 \\[5pt] 0\;\;\;\; & b \!\left( \dfrac{1}{2} + a \right) \end{bmatrix} \\[5pt] \textbf{D}_a & = - \pi \rho b U_\infty \begin{Bmatrix} - 1 \\[5pt] b \!\left( \dfrac{1}{2} + a \right) \end{Bmatrix} \bar{\textbf{b}}. \end{aligned}\end{equation}
For the formulation with the Wagner’s function, the aerodynamic output matrices are
\begin{align} \textbf{M}_a &= \pi \rho b^2 \!\left[\begin{array}{c@{\quad}c} - 1 & ba \\[5pt] ba & - b^2 \!\left( \dfrac{1}{8} + a^2 \right)\end{array} \right]\nonumber\\[5pt] \textbf{C}_a & = \pi \rho b U_\infty \!\left[\begin{array}{c@{\quad}c} - 2 \phi_w(0) & - b - 2 b \!\left( \dfrac{1}{2} - a \right) \phi_w(0) \\[5pt] 2 b \!\left( \dfrac{1}{2} + a \right) \phi_w(0) & - b^2 \!\left( \dfrac{1}{2} - a \right) + 2 b^2 \!\left( \dfrac{1}{4} - a^2 \right) \phi_w(0) \end{array} \right]\nonumber\\[5pt] \textbf{K}_a &= 2 \pi \rho b U_\infty^2 \!\left[\begin{array}{c@{\quad}c} 0 & - \phi_w(0) \\[5pt] 0 & b \!\left( \dfrac{1}{2} + a \right) \phi_w(0) \end{array} \right] \nonumber\\[5pt] \textbf{D}_a & = 2 \pi \rho b U_\infty \!\left[\begin{array}{c@{\quad}c} - 1 & - 1 \\[5pt] b \!\left( \dfrac{1}{2} + a \right) & b \!\left( \frac{1}{2} + a \right) \end{array} \right]\end{align}
With the RFA, the aerodynamic output matrices are
\begin{equation} \begin{aligned} \textbf{M}_a & = \frac{1}{2} \rho U_\infty^2 \bar{\textbf{A}}_2 \!\left( \frac{b}{U_\infty} \right)^2 \\[5pt] \textbf{C}_a & = \frac{1}{2} \rho U_\infty^2 \bar{\textbf{A}}_1 \!\left( \frac{b}{U_\infty} \right) \\[5pt] \textbf{K}_a & = \frac{1}{2} \rho U_\infty^2 \bar{\textbf{A}}_0 \\[5pt] \textbf{D}_a & = \frac{1}{2} \rho U_\infty^2 \begin{bmatrix} \bar{\textbf{A}}_3 & \cdots & \bar{\textbf{A}}_{n+2} \end{bmatrix}, \end{aligned}\end{equation}
where
${\bar{\textbf{A}}_i}$
(
$i = 0,1, \cdots, n + 2$
) are the fitted aerodynamic influence matrices using the least-square method.
