# Flutter instability and active aeroelastic control with time delay for a two-dimensional airfoil

J.-J. Sinou

► To cite this version:

J.-J. Sinou. Flutter instability and active aeroelastic control with time delay for a two-dimensional airfoil. European Journal of Mechanics - A/Solids, 2022, 92, pp.104465. 10.1016/j.euromechsol.2021.104465. hal-03463910

HAL Id: hal-03463910

https://hal.science/hal-03463910

Submitted on 2 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Flutter instability and active aeroelastic control with time delay for a two-dimensional airfoil

J.-J. Sinou*

Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, École Centrale de Lyon, 69134 Écully, France
Institut Universitaire de France, 75005 Paris, France

## ARTICLE INFO

**Keywords:**

Stability analysis

Airfoil flutter

Active control

Time delay

## ABSTRACT

This paper studies the stability of a two-dimensional controlled airfoil system with time delay. First of all, a hybrid approach based on the frequency sweep test and the eigenvalue problem approximation using the Taylor series expansion of the delayed term is proposed to predict the critical time delay as well as the angular frequency of flutter instability. The efficiency of the proposed approach is illustrated by numerical examples for the prediction of self-sustaining vibrations of a phenomenological airfoil model with two degrees of freedom. One point of interest is to bring an understanding of the physical phenomena of the dynamic behavior involved in the appearance of flutter for the controlled airfoil system with time delay. Numerical results indicate that although one active controller is effective in suppressing airfoil flutter in a controlled system without time delay, control failure may happen depending on the time delay in control design. Secondly not only the effectiveness but also some drawbacks of two proposed controllers, namely the linear quadratic regulator state-feedback controller and the pole placement approach, are discussed through numerical simulations.

## 1. Introduction

Flutter is one of the classical problems of self-excited systems and one of the most important instability phenomenon in aeroelasticity. Due to interaction between the structural dynamic of the mechanical airfoil system and the surrounding airflow, instability due to flutter can occur, leading to a reduction in aircraft performance and, in the worst case, catastrophic failure of the structure.

Several solutions based on passive or active vibration control (Borglund and Kuttenkeuler, 2002; Lee et al., 2007; Tsushima and Su, 2017; Malher et al., 2017; Ouyang et al., 2021) have been recommended and studied for the suppression of flutter instability. The use of active feedback control has attracted the attention of many researchers (Ko et al., 1999; Platanitis and Strganac, 2004; Bhoi and Singh, 2005). A brief review on the available linear and nonlinear control techniques and capabilities for aeroelastic structural systems can also be found in Librescu and Marzocca (2005). Generally speaking the first step in engineering is to consider idealized controls for which no time delay is considered. In reality, time delays are often encountered in various engineering systems with active control. Some short time delays in control loop are unavoidable because of the dynamics involved in feedback controls to stabilize a system and the time it takes for effects to propagate through system components. It is well-known that the presence of this time delay may result in degradation in the control

efficiency and instability of the control system. Therefore, it is crucial to understand the effect of delays on the control system and the potential conditions for the unwanted occurrence of flutter instability in a controlled airfoil system. This topic has been the subject of particular attention by some researchers in recent years (Ramesh and Narayanan, 2001; Yuan et al., 2004; Librescu et al., 2005; Zhao, 2009, 2011; Luo et al., 2016).

From a theoretical point of view, the introduction of time delay in the dynamic equation of a mechanical system leads to a more complex analysis than in the case of the delay-free controlled system. Indeed, in the case of a linear time-invariant system, the characteristic equation becomes transcendental because of the exponential functions associated with the time-delays. To overcome this difficulty, a number of methodologies have been proposed to assess the stability of time-delay systems. In the present work, we use the classical substitute for the exponential time-delay terms by considering the approximation based on Taylor series expansion in the characteristic equation with time delay. One of the originality of the proposed work is to answer this question by using a hybrid strategy (Sinou and Chomette, 2021) based on the combination of the frequency sweeping test (Gu et al., 2003) with an expansion based approach of the delayed term for approximating the eigenvalues problem.

<footnote>* Correspondence to: Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, École Centrale de Lyon, 69134 Écully, France.
E-mail address: jean-jacques.sinou@ec-lyon.fr.</footnote>
The second objective of the present work is to predict the flut-
ter instability of a controlled airfoil system with time delay and to
propose a discussion on the impact of two active controllers (i.e. the
Linear Quadratic Regulator (LQR) state-feedback controller and the
pole placement technique) on the critical time delay.

This paper is organized as follows. First of all, the airfoil model
without control and the associated stability analysis are briefly dis-
cussed. Secondly the dynamic equation of controlled airfoil flutter is
established by considering two strategies based on the LQR state feed-
back design or the pole placement technique. Then the stability analysis
of time-delay systems and more specifically the determination of the
critical time delay based on a hybrid approach are undertaken. Finally
simulation results are discussed to demonstrate the relevance of the
proposed hybrid approach for predicting efficiency flutter instability
of delay-controlled airfoil system. Additionally some comments are
given on the impact of the choice of controllers input parameters on
the value of the critical time delay for the controlled airfoil system
and, consequently, on the occurrence of the flutter instability and its
dynamic behavior.

2. Airfoil model and stability analysis without control

The mechanical system under study corresponds to a two-degrees-
of-freedom airfoil model subjected to quasi-steady aerodynamic load-
ing (Zhao, 2009), as shown in Fig. 1. Considering the two degrees of
freedom which are the plunging deflection h (positive in the downward
direction) and the pitching angle α about the elastic axis (positive nose
up), the aeroelastic equations of motion can be expressed as

$$
\begin{bmatrix} m & m x_{\alpha} \\ m x_{\alpha} & I_{\alpha} \end{bmatrix} \begin{bmatrix} \ddot{h} \\ \ddot{\alpha} \end{bmatrix} + \begin{bmatrix} c_h & 0 \\ 0 & c_{\alpha} \end{bmatrix} \begin{bmatrix} \dot{h} \\ \dot{\alpha} \end{bmatrix} + \begin{bmatrix} k_h & 0 \\ 0 & k_{\alpha} \end{bmatrix} \begin{bmatrix} h \\ \alpha \end{bmatrix} = \begin{bmatrix} -L_{qs} \\ T_{qs} \end{bmatrix} \quad (1)
$$

where $L_{qs}$ and $T_{qs}$ are the quasi-steady aerodynamic lift and moment that can be expressed for a two-dimensional and incompressible flow such as

$$
L_{qs} = \pi \rho_a b s_p \left( b \left( \ddot{h} + V \dot{\alpha} - b \bar{a} \dot{\alpha} \right) + 2V \left( V \alpha + \dot{h} + b \left( \frac{1}{2} - \bar{a} \right) \dot{\alpha} \right) \right) \quad (2)
$$

$$
T_{qs} = \pi \rho_a b^2 s_p \left( b \bar{a} \dot{h} - V b \left( \frac{1}{2} - \bar{a} \right) \dot{\alpha} - b^2 \left( \frac{1}{8} + \bar{a}^2 \right) \ddot{\alpha} \right. \\
\left. \qquad + 2V \left( \frac{1}{2} + \bar{a} \right) \left( V \alpha + \dot{h} + b \left( \frac{1}{2} - \bar{a} \right) \dot{\alpha} \right) \right) \quad (3)
$$

The definitions and values of all the parameters are given in Table 1. V that defines the flow speed of the airfoil system is classically considered as the unfolding parameter for the determination of flutter instability in airfoil structures. To be noted that the expression of the aerodynamic model and the aeroelastic forces in Eqs. (1), (2) and (3) come from the original model of Theodorsen (1935). In the present case the simplified proposed model is suitable for incompressible flow and by considering the well known Theodorsen's circulation function C(k) equal to 1 (with the reduced frequency k given by k = ωb/V where ω is the circular frequency of the body motion, see for example Brunton and Rowley, 2013 for more details). It is admitted that the quasi-steady aerodynamics conditions corresponds to 0 ≤ k ≤ 0.05 which means more precisely that the unsteady effects are usually small.

The second-order differential equations describing the dynamic be-
havior of the system can be rewritten in the following matrix form

$$
M\ddot{x}(t) + C\dot{x}(t) + Kx(t) = 0 \quad (4)
$$

where $\mathbf{x} = [h \ \alpha]^T$. The mass, damping and stiffness matrices ($\mathbf{M}, \mathbf{C}, \mathbf{K}$) include the contribution of the structural airfoil and the quasi-steady aerodynamic lift and moment such as

$$
\mathbf{M} = \mathbf{M}_s + \mathbf{M}_{\text{aero}} = \begin{bmatrix} m & m x_{\alpha} \\ m x_{\alpha} & I_{\alpha} \end{bmatrix} + \begin{bmatrix} \pi \rho_a b^2 s_p & -\pi \rho_a b^3 s_p \bar{a} \\ -\pi \rho_a b^3 s_p \bar{a} & \pi \rho_a b^4 s_p \left(\frac{1}{8} + \bar{a}^2\right) \end{bmatrix} \quad (5)
$$

$$
\mathbf{C} = \mathbf{C}_s + \mathbf{C}_{\text{aero}} = \begin{bmatrix} c_h & 0 \\ 0 & c_{\alpha} \end{bmatrix} + \begin{bmatrix} 2\pi\rho_a b s_p V & \pi\rho_a b^2 s_p V \left(1 + 2\left(\frac{1}{2} - \bar{a}\right)\right) \\ -2\pi\rho_a b^2 s_p V \left(\frac{1}{2} + \bar{a}\right) & -2\pi\rho_a b^3 s_p V \bar{a} \left(\frac{1}{2} - \bar{a}\right) \end{bmatrix} \quad (6)
$$

$$
\mathbf{K} = \mathbf{K}_s + \mathbf{K}_{\text{aero}} = \begin{bmatrix} k_h & 0 \\ 0 & k_\alpha \end{bmatrix} + \begin{bmatrix} 0 & 2\pi\rho_a b s_p V^2 \\ 0 & -2\pi\rho_a b^2 s_p V^2 \left(\frac{1}{2} + \bar{a}\right) \end{bmatrix} \quad (7)
$$

where $\mathbf{M}_s$, $\mathbf{C}_s$ and $\mathbf{K}_s$ defined the structural mass, damping and stiffness matrices, respectively (we have $\mathbf{M}_s = \mathbf{M}_s^T$, $\mathbf{C}_s = \mathbf{C}_s^T$ and $\mathbf{K}_s = \mathbf{K}_s^T$) while the matrices $\mathbf{M}_{\text{aero}}$, $\mathbf{C}_{\text{aero}}$ and $\mathbf{K}_{\text{aero}}$ correspond to the additional contributions due to quasi-steady aerodynamic flow conditions (we have $\mathbf{M}_{\text{aero}} = \mathbf{M}_{\text{aero}}^T$, $\mathbf{C}_{\text{aero}} \neq \mathbf{C}_{\text{aero}}^T$ and $\mathbf{K}_{\text{aero}} \neq \mathbf{K}_{\text{aero}}^T$).

The stability of the airfoil system can be determined by considering
the eigenvalue problem corresponding to Eq. (1) that can be written as
follow

$$
(\lambda_k^2 \mathbf{M} + \lambda_k \mathbf{C} + \mathbf{K}) \Phi_k = \mathbf{0} \quad (8)
$$

where $(\lambda_k, \Phi_k)$ are the eigenvalues and eigenvectors of the problem. Stability is determined by considering the complex eigenvalues $\lambda_k$. Therefore, three situations are possible in relation to the concept of stability for the airfoil system :

* if all eigenvalues have negative real parts, the system is stable (i.e. $\Re(\lambda_k) < 0 \ \forall \lambda_k$).

* if at least one eigenvalue has positive real part with the associated imaginary part not equal to zero (i.e. $\exists \lambda_k \in \Re(\lambda_k) > 0$ and $\Im(\lambda_k) \neq 0$), the system is unstable and the airfoil encounters flutter. The imaginary part of the associated positive eigenvalue defines the angular frequency of the unstable mode. The value of the flow speed $V$ for which the airfoil system becomes unstable is denoted by the flutter speed $V_f$. The associated angular frequency (classically referred to as the flutter angular frequency) is denoted by $\omega_f$.

* if at least one eigenvalue has positive real part with the associated imaginary part equal to zero (i.e. $\exists \lambda_k \in \Re(\lambda_k) > 0$ and $\Im(\lambda_k) = 0$), the system is unstable and the airfoil encounters divergence. The value of the flow speed $V$ for which the airfoil system becomes unstable is denoted by the divergence speed $V_d$.

Fig. 2 displays the evolutions of real parts and angular frequencies as a function of the flow speed V. For V = 0, the two modes, denoted by mode I and mode II, are stable with angular frequencies equal to ωI,0 = 10.94 rad/s and ωII,0 = 37.70 rad/s (and the associated real parts are equal to approximately ℜ(λI,0) = -0.318 and ℜ(λII,0) = -7.066, respectively). As the flow speed V increases, the real part ℜ(λI,V) first decreases and then increases slowly whereas the real part ℜ(λII,V) decreases continuously. The associated angular frequencies ωI,V and ωII,V tend to get closer. The critical Hopf bifurcation point where the system loses stability is observed at the flutter speed Vf = 23.46 m/s. As shown in Fig. 2, one complex eigenvalue crosses the complex plane imaginary axis reflecting the fact that the airfoil system moves from stable to unstable behavior. For clarity, stable modes are indicated by green dots while the unstable mode is indicated by red dots in Fig. 2. The flutter angular frequency ωf at the Hopf bifurcation point is equal to 24.32 rad/s and it comes from the evolution of mode I according to the flow speed. After the Hopf bifurcation point, the angular frequencies of the two modes remain close and cross at the flow speed V = 25.93 m/s (with ωI,V = ωII,V = 25.65 rad/s). Then the two angular frequencies tend to move away from each other for V > 25.93 m/s with ωI,V > ωII,V. For V > Vf the real parts are opposite: one mode is stable (i.e., the associated real part is negative and decreases versus the increase of the flow speed) whereas the other one is unstable (i.e., the associated real part is positive and increases versus the increase of the flow speed). In conclusion, the airfoil system is stable for the flow speed V ∈ [0; Vf] and encounters flutter instability for V > Vf.
<table><caption>**Table 1**

Geometrical and physical parameters of the airfoil system.</caption><thead><tr><th>Parameter</th><th>Notation</th><th>Value</th></tr></thead><tbody><tr><td>Span of the airfoil</td><td>s<sub>p</sub></td><td>1 m</td></tr><tr><td>Airfoil mass</td><td>m</td><td>2.049 kg</td></tr><tr><td>Dimensionless distance of the elastic axis from midchord</td><td>ā</td><td>-0.6847</td></tr><tr><td>Semichord of the airfoil section</td><td>b</td><td>0.135 m</td></tr><tr><td>Torsional stiffness</td><td>k<sub>a</sub></td><td>6.833 N/m/rad</td></tr><tr><td>Plunging stiffness</td><td>k<sub>h</sub></td><td>2844.4 N/m</td></tr><tr><td>Damping coefficient in pitching</td><td>c<sub>a</sub></td><td>0.036 N s/rad</td></tr><tr><td>Damping coefficient in plunging</td><td>c<sub>h</sub></td><td>27.43 N s/m</td></tr><tr><td>Air density</td><td>ρ<sub>a</sub></td><td>1.225 kg/m<sup>3</sup></td></tr><tr><td>Mass moment of inertia of the airfoil about the elastic axis</td><td>I<sub>a</sub></td><td>m<sub>a</sub>x<sub>a</sub><sup>2</sup> + 0.0517 kg m<sup>2</sup></td></tr><tr><td>Distance of the elastic axis from center of mass</td><td>x<sub>a</sub></td><td>(0.0873 - (1 + ā)b) m</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Schematic mechanical system of a two-dimensional airfoil.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Stability analysis of the two-dimensional airfoil without control (green = stable; red = unstable). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 3. Controlled airfoil system

For a given flight condition of $V > V_f$, the flutter instability may be suppressed by applying an active feedback control. However the effect of time delay which represents the lag between the time needed for control gains computation satisfying a specified control objective, is an important issue because stability switch may be observed with increase of time delay. In the following, the airfoil system regulated by state feedback control with time delay is presented. Stability analysis of the controlled system with a single time-delay and a methodology to determine the critical time delay for which the stability switch appears are also discussed.

#### 3.1. Airfoil system regulated by state feedback control with time delay

The dynamics of the airfoil system (4) regulated by state feedback control with time delay $\tau$ can be expressed by

$$M\ddot{x}(t) + C\dot{x}(t) + Kx(t) = b u(t - \tau) \quad (9)$$

$$u(t - \tau) = -f^T \dot{x}(t - \tau) - g^T x(t - \tau) \quad (10)$$

where $b$ is the control input distribution matrix. $u$ is the associated state feedback control vector with $f$ and $g$ being matrices that corresponds to proportional velocity and displacement feedback control gains. To be noted that the time delays for the velocity and displacement state feedback are assumed to be equal.

Substituting the expression of variables separation $x(t) = Ze^{st}$ (where $Z$ is a constant vector) in Eqs. (9) and (10) yields to the transcendental eigenvalue problem (Ram and Mottershead, 2007; Ram et al., 2009)

$$(s^2M + s(C + bf^T e^{-st}) + K + bg^T e^{-st})Z = 0 \quad (11)$$

The non-trivial solutions ($s_i, Z_i \neq 0$) for the closed-loop system with time delay are eigenpairs of Eq. (11). Note that we have to determine the proportional velocity and displacement control gains $f$ and $g$ in the context of airfoil flutter control. Different strategies exist. In the rest of the paper we propose to apply two well known strategies:

* The first one is based on the LQR state feedback design. Consider the continuous time linear dynamic system in the state space form (with $y = [x \ \dot{x}]^T$)

$$\dot{y} = Ay + Bw \quad (12)$$

where $y$ and $w$ represent respectively the state and control vectors of the system. The basic theory of the optimal linear quadratic regulator is to estimate the optimal gain matrix $G$ so that the state-feedback law $w = -Gy$ minimizes the quadratic cost function defined as

$$J_{LQR} = \int_{0}^{\infty} (y^T Q y + w^T R w) dt \quad (13)$$

where $y^T Q y$ and $w^T R w$ define the state cost with weight $Q$ and the control cost with weight $R$. In other words, LQR selects closed-loop poles that balance between state errors and control effort. As a result, the closed-loop poles correspond to the state-space eigenvalues of $A - BG$ computable by solving the Riccati equation (Kwakernaak and Sivan, 1972). Note that the gain matrix $G$ is directly linked to the matrices of proportional velocity and displacement control gains $f$ and $g$ previously defined in Eq. (10). Such a linear-quadratic regulator is simply an automated way to find an appropriate state feedback controller with respect to the initial system under consideration. One of the drawbacks is that there is no initial choice on the behavior of the controlled system as opposed to strategies based on a full state feedback (i.e. the LQR strategy is to place the pole locations so that the closed-loop system optimizes the cost function defined in Eq. (13)).
One of the critical design step in LQR methodology is choosing the real, symmetric, positive definite weights **Q** and **R** in order to tune the regulator. Classically, one can choose diagonal weights applying the Bryson's rule

$$ \mathbf{Q}_{ii} = \frac{1}{(y_i)_{\max}^2} \qquad (14) $$

$$ \mathbf{R}_{ii} = \frac{1}{(w_i)_{\max}^2} \qquad (15) $$

where $(y_i)_{\max}^2$ and $(w_i)_{\max}^2$ define the maximum desired values or the control input for component of the state/actuator signal. Although Bryson's rule generally works well, it is often, from a practical point of view, only the starting point of an iterative trial-and-error design procedure to obtain desirable properties for the closed-loop system.

* A second strategy classically applied is the use of full state feedback, also known as pole placement, which is sometimes preferred because there is a clearer relationship between the controller parameters and the behavior obtained for the controlled system. In the present work the classical pole placement algorithm (Kautsky et al., 1985; Valasek and Olgac, 1995) will be used to place the closed-loop poles at some specific locations by computing a state-feedback gain matrix. Consider the linear dynamic system given in Eq. (12), the pole placement approach computes a gain matrix **G** such that the state feedback **w** = −**G****y** forces the closed-loop poles to the desired locations. In other words, the state-space eigenvalues of **A** − **BG** match the predefined locations affecting the response as desired. In the case of multiple input systems the gain matrix **G** is not unique and choosing the best **G** values is not trivial. In the following, the closed-loop poles locations are chosen so that the complex conjugate eigenvalues associated with the stable mode remains unchanged (i.e. $\lambda_{3,4}^{\text{controlled}} = \lambda_{3,4}^{\text{uncontrolled}}$) while the two others complex conjugate eigenvalues associated with the unstable mode of the uncontrolled systems are placed such as

$$ \lambda_i^{\text{controlled}} = \gamma \Re(\lambda_i^{\text{uncontrolled}}) + \sqrt{-1} \Im(\lambda_i^{\text{uncontrolled}}) \quad \text{for } i = 1, 2 \quad (16) $$

where $\Re()$ and $\Im()$ define the real and the imaginary part of eigenvalues, respectively. The scalar coefficient $\gamma$ corresponds to the control gain applied on the real part of the uncontrolled poles. If $\gamma < 0$ the controlled system without time delay is stable, assuming that the uncontrolled system was initially unstable. Moreover the algorithm proposed by Kautsky et al. (1985) is applied to find a robust solution, such that the assigned closed-loop poles are insensitive to perturbations in the system data (i.e. the provided solution minimizes the sensitivity of the closed-loop poles to perturbations in the matrices **A** and **B**).

## 3.2. Stability analysis of time-delay systems and estimation of the critical time delay $\tau_c$

In this section we discuss a methodology to analyze the stability of the airfoil system regulated by state feedback control with time delay and to predict the critical time delay as well as the angular frequency of the associated unstable mode. The proposed approach is based on the determination of the specific point of the frontier stability curve provided by the frequency sweeping tests (Gu et al., 2003) that coincides with the evolution of the generalized eigenvalue problem for a Taylor series expansion up to the $p$th order of the closed-loop system with time-delay (Sinou and Chomette, 2021).

### 3.2.1. Preamble on the Taylor series expansion of the delayed term and the stability of the controlled system with time delay

First of all the eigenvalues of Eq. (11) can be approximated by converting the transcendental eigenvalue problem to an approximated eigenvalues problem obtained by Taylor series expansion. Substituting the Taylor series expansion of the delayed term $\mathbf{u}(t - \tau)$ up to the $p$th order

$$ \mathbf{u}(t - \tau) = \sum_{n=0}^{p} \frac{(-1)^n}{n!} \tau^n \mathbf{u}^{(n)}(t) \qquad (17) $$

where $\mathbf{u}^{(n)}$ denotes the $n$th derivative of $\mathbf{u}$ with respect to time, the closed-loop system defined by Eqs. (9) and (10) can be approximated by

$$ \mathbf{M}\ddot{\mathbf{x}}(t) + \mathbf{C}\dot{\mathbf{x}}(t) + \mathbf{K}\mathbf{x}(t) = -\mathbf{b}\mathbf{f}^T \sum_{n=0}^{p} \frac{(-1)^n}{n!} \tau^n \mathbf{x}^{(n+1)}(t) - \mathbf{b}\mathbf{g}^T \sum_{n=0}^{p} \frac{(-1)^n}{n!} \tau^n \mathbf{x}^{(n)}(t) \qquad (18) $$

The approximation of the generalized eigenvalues problem (11) by using the $p$th order Taylor expansion yields to the first-order realization defined by

$$ (\mathbf{P}_p - s\mathbf{Q}_p) \mathbf{V}_p = \mathbf{0} \qquad (19) $$

where $\mathbf{P}_p$ and $\mathbf{Q}_p$ are $(p+1) \times (p+1)$ matrices and $\mathbf{V}_p$ is a $(p+1)$ vector. Expression of $\mathbf{P}_p$ is given in Eq. (20) (see Box I) and expressions of $\mathbf{Q}_p$ and $\mathbf{V}_p$ in Eq. (21),

$$ \mathbf{Q}_p = \operatorname{diag} \left( \mathbf{I}, \dots, \mathbf{I}, (-1)^{(p-1)} \frac{\tau^p}{p!} \mathbf{b} \mathbf{f}^T \right); \quad \mathbf{V}_p = \left[ \mathbf{Z} \quad s\mathbf{Z} \quad \dots \quad s^p \mathbf{Z} \right]^T \quad (21) $$

As previously explained in Sinou and Chomette (2021) and Ram et al. (2009), the controlled system without time delay (i.e. $\tau = 0$) has only $2q$ eigenvalues (with $q$ the number of degrees of freedom). On the contrary, if we consider the same system with a time delay $\tau \neq 0$ there is an infinity of eigenvalues but only $2q$ of them affect the dynamics of the system. These last ones are called the primary eigenvalues while the others are secondary eigenvalues. One of the issues is therefore to be able to distinguish between primary and secondary eigenvalues. A solution can be to follow the evolution of the primary resonances by continuously increasing the time delay $\tau$ from the reference controlled system without time delay. When the time delay is increased, if a real part of the eigenvalues related to the primary resonances becomes greater than zero, the system becomes unstable and the corresponding time delay defines the critical time delay $\tau_c$. Moreover, the imaginary part associated with the eigenvalue with positive real part obviously corresponds to the angular frequency of the unstable mode. However considering high order Taylor series expansion may lead to ill-conditioned matrix $\mathbf{Q}_p$ for $\tau \to 0$ making it impossible in some cases to track primary resonances. An alternative is then to combine the determination of the eigenvalues of the system (19) with the determination of the stability boundary defined by the frequency sweeping tests (Gu et al., 2003) that is briefly described in the following section.

### 3.2.2. Preamble on the frequency sweeping tests

Before providing a methodology for rapidly detecting the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ for which the controlled system becomes unstable, stability criteria based on the frequency sweeping tests is briefly discussed. We consider the reformulation of the problem (9) and (10) in the state space form $\dot{\mathbf{y}}(t) = \mathbf{A}_0\mathbf{y}(t) + \mathbf{A}_1\mathbf{y}(t - \tau)$ where $\mathbf{y} = [\mathbf{x} \quad \mathbf{x}]^T$.

Supposing that the controlled airfoil system is stable at $\tau = 0$. Let $\operatorname{rank}(\mathbf{A}_1) = r$ and define

$$ \tau_{lim}^{i} = \begin{cases} \min_{1 \le k \le n} \frac{\theta_k^i}{\omega_k^i} & \text{if } \lambda_i \left( (j\omega_k^i \mathbf{I} - \mathbf{A}_0) , \mathbf{A}_1 \right) = e^{-j\theta_k^i} \\ \infty & \text{if } \rho \left( (j\omega \mathbf{I} - \mathbf{A}_0) , \mathbf{A}_1 \right) > 1 \quad \forall \omega > 0 \end{cases} \qquad (22) $$
$$ \mathbf{P}_{\rho} = \begin{bmatrix} \mathbf{0} & \mathbf{I} & \mathbf{0} & \dots & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{I} & \dots & \mathbf{0} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \dots & \mathbf{I} \\ \mathbf{K} + \mathbf{b}\mathbf{g}^T & \mathbf{C} + \mathbf{b}\mathbf{f}^T - \tau\mathbf{b}\mathbf{g}^T & \mathbf{M} - \tau\mathbf{b}\mathbf{f}^T + \frac{\tau^2}{2}\mathbf{b}\mathbf{g}^T & \dots & (-1)^{(p-1)} \frac{\tau^{p-1}}{(p-1)!} \mathbf{b}\mathbf{f}^T + (-1)^p \frac{\tau^p}{p!} \mathbf{b}\mathbf{g}^T \end{bmatrix} \quad (20) $$

Box I.

and

$$ \tau_{lim} = \min_{1 \le i \le r} \tau_{lim}^{i} \quad (23) $$

$\rho(\mathbf{A})$ corresponds to the spectral radius of the matrix $\mathbf{A}$. $\tau_{lim}$ corresponds to the delay margin of the controlled system: the controlled system with time delay defined in Eqs. (9) and (10) is stable for all $\tau \in [0; \tau_{lim}]$ but becomes unstable at $\tau = \tau_{lim}$.

### 3.2.3. Addressing the problem of determining the critical time delay $\tau_c$ from the delay margin $\tau(\omega)$

The aim of this section is to discuss an efficient strategy for addressing the problem of obtaining the delay margin for a predefined angular frequency interval of interest [$\omega_{min}; \omega_{max}$] and to predict the critical time delay $\tau_c$ (for which the controlled system becomes unstable) and the associated angular frequency $\omega_c$ that corresponds to the angular frequency of the unstable mode (i.e. the airfoil controlled system with time delay $\tau_c$ will oscillate at $\omega_c$).

Based on the previous Sections 3.2.1 and 3.2.2, the following statements can be made:

* The generalized eigenvalues problem given in Eq. (11) can be rewritten as the pseudo-polynomial

$$ \det(s\mathbf{I} - \mathbf{A}_0 - \mathbf{A}_1 e^{-\tau s}) = 0 \quad (24) $$

and Eq. (19) corresponds to its approximation by using the $p$th order Taylor expansion.

* Considering a specific angular frequency interval of interest [$\omega_{min}; \omega_{max}$], if

$$ \rho((j\omega\mathbf{I} - \mathbf{A}_0), \mathbf{A}_1) > 1 \quad \text{for all } \omega \in [\omega_{min}; \omega_{max}] \quad (25) $$

it can be concluded that the controlled system is stable for all $\tau \in [0; \infty)$. Otherwise the computation of the problem (22) generates a nonzero estimate of the delay margin over the finite interval of interest $\omega = [\omega_{min}; \omega_{max}]$. This delay margin function (denoted by $\tau(\omega)$ for simplicity in the rest of the paper) can be determined by sampling the angular frequency interval of interest [$\omega_{min}; \omega_{max}$], and at each sample, computing the generalized eigenvalues such as $\tau(\omega) = \min_{1 \le i \le n} \frac{\theta_i}{\omega}$ with $\theta_i = -j^{-1}\log(\lambda_i)$ and $\theta_i \in [0; 2\pi]$ for any eigenvalue $\lambda_i$ of the matrix $(j\omega\mathbf{I} - \mathbf{A}_0)^{-1}\mathbf{A}_1$.

Considering the evolutions of eigenvalues versus $\tau$, defined by $\lambda_i(\tau)$, the crossing of the characteristic roots with the delay margin $\tau(\omega)$ represents the transition from the stable region to the unstable one. However multiple crossings may exist and since only the primary eigenvalues are of interest, there are only a finite number of crossing with $\tau(\omega)$ to be considered (to be more precise there are $q$ crossings of interest). Moreover the pair of interest $(\tau, \omega)$ that gives the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ of the self-excited controlled system that is subjected to flutter instability for $\tau > \tau_c$ corresponds to the coincidence between the delay margin $\tau(\omega)$ and the eigenvalues $\lambda_i(\tau)$ such as

$$ \Re(\lambda_i(\tau_-(\omega))) < 0 \quad \text{and} \quad \Re(\lambda_i(\tau_+(\omega))) > 0 \quad (26) $$

and

$$ \frac{|\Re(\lambda_i(\tau(\omega))) - \omega|}{\omega} < \epsilon \quad (27) $$

In fact $(\tau_c, \omega_c)$ corresponds to the crossing of the delay margin $\tau(\omega)$ with the primary eigenvalues of the controlled system and Eq. (26) reflects the fact that one real part of eigenvalues becomes greater than zero at this crossing. Eq. (27) verifies that the eigenvalue $\lambda_i$ is in the neighborhood of the crossing between the delay margin $\tau(\omega)$ and the eigenvalues $\lambda_i(\tau)$. This second criterion defined in Eq. (27) can be seen as an estimator of the quality of the prediction of $(\tau_c, \omega_c)$. The more this criterion tends towards zero, the better the prediction will be. $\epsilon$ corresponds to the residual from which the quality of the prediction is assumed to be sufficient.

Note that it is not necessary to compute the eigenvalues and to follow their evolution on the full range of interest $\tau = [0; \tau_{max}]$ with $\tau_{max}$ the maximum value fixed for the time delay (or more specifically $\tau = [\tau_{lim}; \tau_{max}]$ due to the fact that the controlled system is stable for all $\tau \in [0; \tau_{lim}]$, as previously discussed in Section 3.2.2). Indeed the computation of the eigenvalues is only done on the curve corresponding to the delay margin $\tau(\omega)$ (available through the frequency sweeping test method, see Section 3.2.2) for the prediction of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$. If the $p$th order Taylor expansion is not sufficient and provides a not precise enough approximation of the generalized eigenvalues of Eq. (19), it is possible that no pair $(\tau, \omega)$ checks Eqs. (26) and (27). In this case, it is then necessary to consider solving the problem (19) with a Taylor expansion of a higher order. The corollary is that if a pair $(\tau, \omega)$ verifies the relations (26) and (27) for a given $p$th order Taylor expansion, then this indicates that this order approximation is sufficient to provide reliable information about the stability of the controlled system with time delay and the prediction of the critical time delay $\tau_c$ for which the controlled system becomes unstable. Of course, this is also directly related to the choice of the value of the parameter $\epsilon$.

We can also point out that, in the case of a controlled system with only one unstable mode, only one pair $(\tau_c, \omega_c)$ normally verifies the relations defined by Eqs. (26) and (27) (i.e. the primary eigenvalue associated with the unstable mode of the controlled system). If by accident the solution of Eqs. (26) and (27) is not unique this means that a secondary eigenvalue also intersects the delay margin $\tau(\omega)$ with a change of sign of its real part from negative to positive. In practice, this case is very unlikely. However, if this happens, it is still possible to select the physical case of interest by distinguishing the primary and secondary eigenvalues via the evolution of these eigenvalues by continuously decreasing the time delay $\tau$ from the critical time delay identified for each eigenvalue. When $\tau \to 0$ the eigenvalue of the physical case of interest will be consistent with the eigenvalue of the reference controlled system without time delay.

## 4. Application and numerical examples

In this section the relevance and effectiveness of the strategy proposed in Section 3.2.3 will be first discussed. We remind that the airfoil system controlled without time delay is stable and that the main objective is to estimate the critical time delay $\tau_c$ for which the system
become unstable, as well as the angular frequency $\omega_c$ of unstable oscillations. Moreover, for both LQR control and pole placement method, a more detailed study will be carried out in order to highlight some recommendations to be taken into account due to the impact of the control gain on the critical delay $\tau_c$.

### 4.1. Validation of the proposed methodology for predicting the critical time delay $\tau_c$ from the delay margin $\tau(\omega)$

We will only consider the case of the controlled airfoil system via LQR state feedback. The validation of the proposed methodology for predicting the critical time delay $\tau_c$ from the delay margin $\tau(\omega)$ would be the same by using the full state feedback controller via the pole placement technique.

So to demonstrate the efficiency of the proposed strategy in Section 3.2.3 for predicting the critical time delay $\tau_c$ for the controlled airfoil system, we consider the case of the unstable system without control for $V = 30$ m/s. As a reminder the uncontrolled airfoil system is stable for $V \in [0; V_f]$ with $V_f = 23.46$ m/s and it becomes unstable for $V > V_f$. The values of control gains f and g for the case under study are provided in Table 2. These LQR controllers are obtained with $R = I_{2 \times 2}$ and $Q = \text{diag}(\frac{1}{0.01^2}, \frac{1}{0.01^2}, 0, 0)$. In this case the goal is to reach the position $(h, \alpha) = (0, 0)$ with the maximal value $h_{max} = 0.01$ m and $a_{max} = 0.01$ rad. Note that we do not constrain the value of these two variables with such a condition, but we encourage the system to keep the variables as close to zero as possible. The associated poles for the uncontrolled and controlled systems without time delay are also provided in Table 2 (note that the state-space eigenvalues $\lambda_{1,2}$ will always be associated with the unstable mode of the uncontrolled system in the rest of the paper for clarity). Table 3 gives estimation of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ for which the controlled airfoil system becomes unstable. These results are based on the proposed strategy developed in Section 3.2.3 by considering various orders for the Taylor expansion in Eqs. (19), (20) and (21). Table 3 provides also the values of $\mathfrak{R}(\lambda_i(\tau_-(\omega)))$ and $\mathfrak{R}(\lambda_i(\tau_+(\omega)))$ which are used as indicators to detect the transition from stable to unstable behavior of the controlled airfoil system (i.e. these two quantities must be of opposite signs, see Eq. (26)), and the value $|\mathfrak{R}(\lambda_i(\tau(\omega))-\omega)|$ which is the second criterion to check for ensuring that one is in the vicinity of the intersection between the delay margin $\tau(\omega)$ and the primary eigenvalue $\lambda_i$ whose real part changes sign. It is very clear that by increasing the order of expansion of Taylor the value $|\mathfrak{R}(\lambda_i(\tau(\omega))-\omega)|$ tends to zero and thus we improve the detection of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$.

In order to verify the validity of the estimate $(\tau_c, \omega_c)$ based on the proposed methodology, two kinds of verification are conducted:

* at first the evolution of the delay margin $\tau(\omega)$ based on the frequency sweeping approach as well as the evolution of the primary and secondary eigenvalues based on the Taylor expansion of $\rho$th order (i.e. imaginary parts of eigenvalues versus the time delay $\tau$) are shown in Figs. 3. The symbols o correspond to the angular frequency of modes I and II (with $\omega_{I,V} > \omega_{II,V}$ at $V = 30$ m/s) for the controlled airfoil system without time delay. It appears without any ambiguity that the prediction of $(\tau_c, \omega_c)$ (defined by different colored crosses according to the Taylor expansion orders) is indeed in the vicinity of the crossing of the delay margin with the primary eigenvalues of the controlled system. For orders 5 and 6, the estimation of the crossover is excellent which is in total agreement with the value of the proposed criterion defined in Eq. (27) for the estimation of the quality of the prediction (see also the values $|\mathfrak{R}(\lambda_i(\tau(\omega))-\omega)|$ in Table 3). To be noted that the proposed prediction of $(\tau_c, \omega_c)$ for order 3 is a little bit worse (less close to the crossover). This fact is also in agreement with the results presented in Table 3 (indeed the quality of the prediction $(\tau_c, \omega_c)$ was estimated to be worse for this order). For

<table><caption>**Table 2**
Poles for the uncontrolled and controlled systems without time delay and estimation of f and g for the controlled airfoil system with LQR at $V = 30$ m/s.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>λ<sub>1,2</sub><sup>uncontrolled</sup></td><td>4.40 ± 27.14i</td></tr><tr><td>λ<sub>3,4</sub><sup>uncontrolled</sup></td><td>-20.71 ± 24.50i</td></tr><tr><td>λ<sub>1,2</sub><sup>controlled</sup></td><td>-6.46 ± 27.95i</td></tr><tr><td>λ<sub>3,4</sub><sup>controlled</sup></td><td>-20.77 ± 24.09i</td></tr><tr><td>f</td><td>[16.91 7.94 ; 7.94 4.92]</td></tr><tr><td>g</td><td>[375.53 255.13 ; 48.00 85.03]</td></tr></tbody></table>

<table><caption>**Table 3**
Estimation of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ for the controlled airfoil system with LQR at $V = 30$ m/s.</caption><thead><tr><th>Order</th><th>τ<sub>c</sub> (s)</th><th>ω<sub>c</sub> (rad/s)</th><th>|∫(λ<sub>i</sub>(τ(ω))-ω)|/ω</th><th>∫(λ<sub>i</sub>(τ<sub>-(</sub>(ω)))</th><th>∫(λ<sub>i</sub>(τ<sub>+(</sub>(ω)))</th></tr></thead><tbody><tr><td>3</td><td>0.0270</td><td>38.29</td><td>0.32620</td><td>-0.00044</td><td>0.00080</td></tr><tr><td>4</td><td>0.0281</td><td>38.61</td><td>0.07774</td><td>-0.00153</td><td>0.00226</td></tr><tr><td>5</td><td>0.0280</td><td>38.78</td><td>0.01410</td><td>-0.00077</td><td>0.00012</td></tr><tr><td>6</td><td>0.0280</td><td>38.77</td><td>0.00144</td><td>-0.00131</td><td>0.00172</td></tr></tbody></table>

the interested reader this can be simply explained by the fact that the time delay value for which the controlled system becomes unstable via the Taylor expansion is not well predicted for a Taylor expansion of order 3.

* on the other hand a temporal integration of the controlled airfoil system given in Eqs. (9) and (10) has been carried out for two values of time delay $\tau$ just before and just after the critical time delay $\tau_c$ (i.e. $\tau = \tau_c(1 \pm \chi)$ with $\chi = 0.01$). Results are shown in Figs. 4. Not surprisingly growing (decreasing, respectively) of the initial oscillations are observed for $\tau = \tau_c(1 + \chi)$ ($\tau = \tau_c(1 - \chi)$, respectively) indicating that the controlled airfoil system with time delay $\tau$ is unstable (stable, respectively). Moreover the angular frequency $\omega$ of the temporal vibrational response of the controlled airfoil system is estimated to be approximately equal to 38.8 rad/s which is in perfect agreement with the prediction $\omega_c$ found by the proposed method (see Table 3). Note that to obtain the vibrational oscillations of the controlled system a small disturbance on the displacement vector around zero is added on the airfoil system.

All these results unambiguously demonstrate the relevance of the proposed strategy to determine not only the critical time delay $\tau_c$ for which the controlled system becomes unstable but also to predict the associated angular frequency $\omega_c$ at which the airfoil system oscillates. It is also illustrated the interest to take into account the impact of the time delay $\tau$ in the control feedback and its impact on the stability of the controlled airfoil system because of the potential unstable behavior that can emerge due to a non-zero time delay $\tau$.

As a complement Fig. 5 shows the evolution of the two primary eigenvalues versus the time delay $\tau$ (by applying a classical stability analysis of the controlled system with time delay via the Taylor series expansion, as indicated in Section 3.2.1). The symbols □ and o define the starting point (i.e. $\tau = 0$ s) and the end point (i.e. $\tau = 0.04$ s) for the evolution of eigenvalues associated to modes I and II, respectively. It clearly appears that the controlled airfoil system without time delay is stable and becomes unstable for $\tau > \tau_c$. The unstable mode corresponds to mode I. The eigenvalue associated with mode II remains practically unchanged for the considered time delay interval.

Note that the simulation results proposed in Figs. 3–5 are not necessary to determine the value of critical time delay $\tau_c$ via the hybrid proposed methodology. They serve only to provide evidence that the hybrid proposed methodology is operating correctly.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution of $(\tau, \omega)$ based on the Taylor expansion of $p$th order for $V = 30$ m/s (--- 3rd order; --- 4th order; --- 5th order; --- 6th order; + prediction of $(\tau_c, \omega_c)$; $- \tau(\omega)$ based on the frequency sweeping approach; ◦ angular frequencies for the controlled airfoil system without time delay) (a) global view (b) zoom for the detection of $(\tau_c, \omega_c)$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Time responses of the controlled system via LQR for $V = 30$ m/s and $\chi = 0.01$ (a) before the critical time delay with $\tau_c (1 - \chi)$ and (b) after the critical time delay with $\tau_c (1 + \chi)$.</figcaption></figure>

## 4.2. Further study and discussion on the use of LQR state-feedback controller

Simulation results on the use of LQR controller by considering different speeds $V$ from 25 m/s to 60 m/s, with an iteration step of 1 m/s are undertaken. The gain matrices $f$ and $g$ are obtained with $R = I_{2 \times 2}$ and $Q = \text{diag}(\frac{1}{0.01^2}, \frac{1}{0.01^2}, 0, 0)$ as for the previous study presented in Section 4.1. Table 4 provides the values of $f$ and $g$ for each value of the speed $V$. It can be noted a continuous and significant evolution of each component of the two matrices according to the speed $V$ at which the LQR control is carried out. The associated poles for the controlled airfoil system are also provided in Table 5. Table 6 gives the prediction of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ for the controlled airfoil system with LQR, as well as the values of both criteria defined in Eqs. (26) and (27) that are used for the estimation of $(\tau_c, \omega_c)$. Note that the residual $\epsilon$ that defines the quality of the prediction is fixed at $10^{-2}$ for the rest of the study. It can already be noticed that the hybrid methodology based on the sweeping frequency test combined with Taylor series expansion of orders 5 or 6 allows the prediction of $(\tau_c, \omega_c)$ with the chosen quality criterion $\epsilon$. Note that the results provided in Table 6 correspond to the lowest order that verifies

$$ \frac{|\Im(\lambda_i(\tau(\omega)) - \omega)|}{\omega} < \epsilon. $$ It is highlighted that the critical time delay $\tau_c$ can evolve significantly depending on the speed $V$ (in fact it is more exactly directly dependent on the two gain matrices $f$ and $g$). In this particular case studied, increasing the speed $V$ decreases the critical time delay $\tau_c$ with the smallest value equal to 0.0150 s (for $V = 60$ m/s) which corresponds to a decrease of more than 2 times compared to the value obtained for $V = 25$ m/s. At the same time a significant increase of the angular frequency $\omega_c$ of the unstable mode is observed with a transition from 33.94 rad/s for $V = 25$ m/s to 69.76 rad/s for $V = 60$ m/s.

As previously mentioned in Section 3.1, one of the crucial steps in the LQR methodology is the choice of $Q$ and $R$. In order to examine the influence of these settings on the critical time delay $\tau_c$, four cases are investigated as shown in Table 7. The first case corresponds to the previous case already studied. The second and third cases are based on 10 times larger and smaller values for $h_{max}$ and $\alpha_{max}$ for the LQR controller, keeping as an objective to reach the position $(h, \alpha) = (0, 0)$. For the last case we choose to limit not only $h$ and $\alpha$ but also the associated speeds $\dot{h}$ and $\dot{\alpha}$. Note that we choose diagonal weights applying the Bryson's rule with $R = I_{2 \times 2}$ for the four cases. These four conditions are not intended to give a complete overview of the influence of the choice of $Q$ but only to illustrate its impact on the
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Evolutions of the two modes versus the time delay $\tau$ (green = stable; red = unstable). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

stability of the controlled airfoil system with time delay. Figs. 6(a) shows the poles placement by considering different speeds $V$ from 25 m/s to 60 m/s for all these cases. The symbol letters “I” and “F” define the initial starting point ($V = 25$ m/s) and the final point ($V = 60$ m/s) for each quantity of interest. The placement of the poles for the uncontrolled airfoil system are also provided by the red symbol o. Figs. 6(b) gives the prediction of both the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ for the four cases. It appears without any ambiguity that the choice of $Q$ can drastically change the value of the time delay from which the controlled airfoil system becomes unstable. Moreover, the angular frequency of unstable oscillations is also directly impacted by this choice which can be a determining factor when designing the airfoil system with respect to its vibratory behavior. For the interested reader, it has also been verified that there is no sensitivity on the generated poles of the controlled system to small perturbations in $Q$ (i.e. $(1 \pm \epsilon)Q$ with $\epsilon = 10^{-8}$).

In conclusion one of the most critical design steps in the LQR methodology is to appropriately tune the regulator for the design space of interest and to select the quantity $Q$ (and $R$ even if this case has not be illustrated for the sake of brevity) to achieve desirable properties for the closed-loop system, including the critical time delay $\tau_c$, which is an essential factor for robust design of system control in the field of mechanical engineering.

### 4.3. Discussion on the use of the pole placement technique

The main objective of this section is to discuss some results on the sensitivity of the critical delay $\tau_c$ for the controlled airfoil system via the pole placement technique and thus some precautions to take into consideration. In the following, the closed-loop poles locations are chosen so that $\lambda_{3,4}^{\text{controlled}} = \lambda_{3,4}^{\text{uncontrolled}}$ (i.e. the complex conjugate eigenvalues associated with the stable mode remains unchanged) while the two others complex conjugate eigenvalues associated with the unstable stable mode of the uncontrolled systems are placed as indicated in Eq. (16) with $\gamma = -1$.

At first we propose to examine the sensibility of the critical time delay $\tau_c$ by performing 100 simulations with very small perturbations on the closed-loop poles locations $\lambda_{1,2}^{\text{controlled}}$ such as $(1 + \epsilon) \lambda_{1,2}^{\text{controlled}}$ with $\epsilon \in [-1; 1] \times 10^{-8}$. We focus more specifically on the controlled airfoil system for $V = 30$ m/s. Fig. 7 shows the associated predictions of $\tau_c$ and $\Delta\omega_c$ defined by $\Delta\omega_c = \frac{\omega_c - \omega_0}{\omega_0}$ where $\omega_c$ is the angular frequency at $\tau_c$ and $\omega_0$ corresponds to the angular frequency of the controlled system without time delay (i.e. $\omega_0 = 27.14$ rad/s at $V = 30$ m/s for mode

<table><caption>**Table 4**
Estimation of $f$ and $g$ for the controlled airfoil system with LQR and by considering different speeds $V$ from 25 m/s to 60 m/s.</caption><thead><tr><th>V (m/s)</th><th>f</th><th>g</th></tr></thead><tbody><tr><td>25</td><td>[12.24 5.48 ; 5.48 3.03]</td><td>[225.94 142.31 ; 30.14 42.74]</td></tr><tr><td>26</td><td>[13.05 5.87 ; 5.87 3.30]</td><td>[250.64 159.92 ; 32.26 48.43]</td></tr><tr><td>27</td><td>[13.94 6.31 ; 6.31 3.62]</td><td>[278.40 180.18 ; 35.05 55.40]</td></tr><tr><td>28</td><td>[14.90 6.81 ; 6.81 4.00]</td><td>[308.84 202.97 ; 38.59 63.78]</td></tr><tr><td>29</td><td>[15.89 7.36 ; 7.36 4.44]</td><td>[341.41 228.06 ; 42.90 73.63]</td></tr><tr><td>30</td><td>[16.91 7.94 ; 7.94 4.93]</td><td>[375.53 255.13 ; 48.00 85.03]</td></tr><tr><td>31</td><td>[17.92 8.56 ; 8.56 5.47]</td><td>[410.63 283.87 ; 53.90 98.00]</td></tr><tr><td>32</td><td>[18.92 9.19 ; 9.19 6.05]</td><td>[446.20 313.93 ; 60.57 112.56]</td></tr><tr><td>33</td><td>[19.89 9.84 ; 9.84 6.68]</td><td>[481.82 345.02 ; 68.00 128.72]</td></tr><tr><td>34</td><td>[20.82 10.50 ; 10.50 7.34]</td><td>[517.13 376.88 ; 76.16 146.48]</td></tr><tr><td>35</td><td>[21.72 11.17 ; 11.17 8.04]</td><td>[551.86 409.25 ; 85.03 165.85]</td></tr><tr><td>36</td><td>[22.57 11.83 ; 11.83 8.77]</td><td>[585.82 441.93 ; 94.56 186.82]</td></tr><tr><td>37</td><td>[23.37 12.49 ; 12.49 9.53]</td><td>[618.83 474.76 ; 104.74 209.38]</td></tr><tr><td>38</td><td>[24.14 13.15 ; 13.15 10.32]</td><td>[650.82 507.57 ; 115.53 233.52]</td></tr><tr><td>39</td><td>[24.86 13.80 ; 13.80 11.13]</td><td>[681.68 540.23 ; 126.90 259.25]</td></tr><tr><td>40</td><td>[25.54 14.45 ; 14.45 11.96]</td><td>[711.38 572.64 ; 138.83 286.53]</td></tr><tr><td>41</td><td>[26.18 15.08 ; 15.08 12.82]</td><td>[739.90 604.70 ; 151.29 315.38]</td></tr><tr><td>42</td><td>[26.78 15.71 ; 15.71 13.68]</td><td>[767.23 636.32 ; 164.25 345.75]</td></tr><tr><td>43</td><td>[27.35 16.32 ; 16.32 14.57]</td><td>[793.36 667.44 ; 177.69 377.66]</td></tr><tr><td>44</td><td>[27.88 16.93 ; 16.93 15.47]</td><td>[818.33 698.02 ; 191.59 411.06]</td></tr><tr><td>45</td><td>[28.38 17.52 ; 17.52 16.38]</td><td>[842.15 723.00 ; 205.92 445.96]</td></tr><tr><td>46</td><td>[28.85 18.11 ; 18.11 17.31]</td><td>[864.85 757.33 ; 220.67 482.32]</td></tr><tr><td>47</td><td>[29.29 18.68 ; 18.68 18.24]</td><td>[886.47 786.00 ; 235.82 520.13]</td></tr><tr><td>48</td><td>[29.70 19.25 ; 19.25 19.18]</td><td>[907.05 813.98 ; 251.36 559.38]</td></tr><tr><td>49</td><td>[30.09 19.80 ; 19.80 20.13]</td><td>[926.61 841.25 ; 267.27 600.03]</td></tr><tr><td>50</td><td>[30.46 20.34 ; 20.34 21.09]</td><td>[945.21 867.80 ; 283.54 642.07]</td></tr><tr><td>51</td><td>[30.80 20.87 ; 20.87 22.05]</td><td>[962.89 893.62 ; 300.16 685.49]</td></tr><tr><td>52</td><td>[31.13 21.39 ; 21.39 23.02]</td><td>[979.67 918.71 ; 317.11 730.25]</td></tr><tr><td>53</td><td>[31.44 21.90 ; 21.90 23.99]</td><td>[995.62 943.06 ; 334.39 776.34]</td></tr><tr><td>54</td><td>[31.72 22.40 ; 22.40 24.96]</td><td>[1010.75 966.69 ; 351.98 823.73]</td></tr><tr><td>55</td><td>[32.00 22.90 ; 22.90 25.94]</td><td>[1025.12 989.59 ; 369.89 872.42]</td></tr><tr><td>56</td><td>[32.26 23.38 ; 23.38 26.91]</td><td>[1038.76 1011.77 ; 388.10 922.37]</td></tr><tr><td>57</td><td>[32.50 23.86 ; 23.86 27.90]</td><td>[1051.71 1033.25 ; 406.62 973.57]</td></tr><tr><td>58</td><td>[32.73 24.32 ; 24.32 28.88]</td><td>[1063.99 1054.02 ; 425.42 1026.00]</td></tr><tr><td>59</td><td>[32.95 24.78 ; 24.78 29.86]</td><td>[1075.66 1074.11 ; 444.52 1079.64]</td></tr><tr><td>60</td><td>[33.17 25.23 ; 25.23 30.84]</td><td>[1086.73 1093.52 ; 463.90 1134.47]</td></tr></tbody></table>

I). Two categories of results are clearly visible. The first one gives an estimate of the critical time delay around 0.0075 s while the second one indicates critical time delay around 0.0175 s, a value more than twice as high as the previous one. Also it appears very clearly that for the first category there is an increase of the angular frequency of the unstable mode $\omega_c$ compared to the value for the controlled system without time delay (i.e. $\Delta\omega_c > 0$) while for the second category a decrease of $\omega_c$ is obtained (i.e. $\Delta\omega_c < 0$). These results illustrate a potential difficulty of using the pole placement technique as they suggest a strong sensitivity of the critical time delay $\tau_c$ according to pole locations. In order to better understand the phenomena involved and the reasons for such different behaviors we suggest comparing the gain matrix $G$ for each simulation. Indeed the pole placement approach computes a gain matrix $G$ in order to force the poles of the closed loop to the desired locations and then this gain matrix $G$ has a direct influence on the estimation of the critical time delay $\tau_c$.

In order to define a similarity scale between two $2q \times q$ gain matrices taken from the 100 previous simulations the correlation coefficient is used. Considering two gain matrices $G^a$ and $G^b$ for the $a$th and $b$th simulations, respectively, the associated scalar correlation coefficient $r_{a,b}$ is defined by

$$r_{a,b} = \frac{\sum_{i=1}^{2q} \sum_{j=1}^{q} (G_{ij}^a - \bar{G}^a)(G_{ij}^b - \bar{G}^b)}{\sqrt{\left(\sum_{i=1}^{2q} \sum_{j=1}^{q} (G_{ij}^a - \bar{G}^a)^2\right) \left(\sum_{i=1}^{2q} \sum_{j=1}^{q} (G_{ij}^b - \bar{G}^b)^2\right)}} \quad (28)$$

where $\bar{G}^a$ and $\bar{G}^b$ corresponds to the mean of all values in array $G^a$ and $G^b$, respectively. $r_{a,b}$ can range from -1 to 1. If $r$ is close to 0,
<table><caption>**Table 5**
Poles for the uncontrolled airfoil system and the controlled airfoil system with LQR and by considering different speeds *V* from 25 m/s to 60 m/s.</caption><thead><tr><th>V (m/s)</th><th>λ<sup>uncontrolled</sup><sub>1,2</sub></th><th>λ<sup>uncontrolled</sup><sub>3,4</sub></th><th>λ<sup>controlled</sup><sub>1,2</sub></th><th>λ<sup>controlled</sup><sub>3,4</sub></th></tr></thead><tbody><tr><td>25</td><td>1.20 ± 25.21i</td><td>-16.02 ± 26.00i</td><td>-6.94 ± 26.86i</td><td>-15.52 ± 24.92i</td></tr><tr><td>26</td><td>1.92 ± 25.68i</td><td>-17.04 ± 25.63i</td><td>-6.51 ± 27.08i</td><td>-16.78 ± 24.76i</td></tr><tr><td>27</td><td>2.59 ± 26.10i</td><td>-18.01 ± 25.31i</td><td>-6.30 ± 27.30i</td><td>-17.89 ± 24.60i</td></tr><tr><td>28</td><td>3.23 ± 26.47i</td><td>-18.94 ± 25.02i</td><td>-6.25 ± 27.52i</td><td>-18.91 ± 24.44i</td></tr><tr><td>29</td><td>3.83 ± 26.82i</td><td>-19.84 ± 24.75i</td><td>-6.31 ± 27.74i</td><td>-19.87 ± 24.27i</td></tr><tr><td>30</td><td>4.40 ± 27.14i</td><td>-20.71 ± 24.50i</td><td>-6.46 ± 27.95i</td><td>-20.77 ± 24.09i</td></tr><tr><td>31</td><td>4.95 ± 27.45i</td><td>-21.55 ± 24.26i</td><td>-6.66 ± 28.17i</td><td>-21.64 ± 23.91i</td></tr><tr><td>32</td><td>5.49 ± 27.74i</td><td>-22.37 ± 24.02i</td><td>-6.92 ± 28.38i</td><td>-22.47 ± 23.73i</td></tr><tr><td>33</td><td>5.97 ± 28.01i</td><td>-23.17 ± 23.78i</td><td>-7.20 ± 28.59i</td><td>-23.28 ± 23.54i</td></tr><tr><td>34</td><td>6.46 ± 28.28i</td><td>-23.96 ± 23.54i</td><td>-7.52 ± 28.80i</td><td>-24.07 ± 23.34i</td></tr><tr><td>35</td><td>6.93 ± 28.54i</td><td>-24.72 ± 23.30i</td><td>-7.84 ± 29.01i</td><td>-24.83 ± 23.13i</td></tr><tr><td>36</td><td>7.38 ± 28.79i</td><td>-25.47 ± 23.06i</td><td>-8.18 ± 29.22i</td><td>-25.59 ± 22.91i</td></tr><tr><td>37</td><td>7.82 ± 29.03i</td><td>-26.21 ± 22.81i</td><td>-8.52 ± 29.42i</td><td>-26.32 ± 22.69i</td></tr><tr><td>38</td><td>8.25 ± 29.27i</td><td>-26.94 ± 22.56i</td><td>-8.87 ± 29.63i</td><td>-27.05 ± 22.46i</td></tr><tr><td>39</td><td>8.67 ± 29.50i</td><td>-27.66 ± 22.30i</td><td>-9.22 ± 29.83i</td><td>-27.76 ± 22.21i</td></tr><tr><td>40</td><td>9.08 ± 29.73i</td><td>-28.36 ± 22.03i</td><td>-9.57 ± 30.03i</td><td>-28.46 ± 21.96i</td></tr><tr><td>41</td><td>9.48 ± 29.95i</td><td>-29.06 ± 21.76i</td><td>-9.92 ± 30.22i</td><td>-29.16 ± 21.70i</td></tr><tr><td>42</td><td>9.87 ± 30.16i</td><td>-29.75 ± 21.47i</td><td>-10.27 ± 30.42i</td><td>-29.84 ± 21.42i</td></tr><tr><td>43</td><td>10.25 ± 30.38i</td><td>-30.43 ± 21.18i</td><td>-10.62 ± 30.61i</td><td>-30.52 ± 21.14i</td></tr><tr><td>44</td><td>10.63 ± 30.58i</td><td>-31.10 ± 20.87i</td><td>-10.96 ± 30.81i</td><td>-31.19 ± 20.84i</td></tr><tr><td>45</td><td>11.00 ± 30.79i</td><td>-31.77 ± 20.55i</td><td>-11.30 ± 31.00i</td><td>-31.85 ± 20.53i</td></tr><tr><td>46</td><td>11.36 ± 30.99i</td><td>-32.43 ± 20.22i</td><td>-11.64 ± 31.18i</td><td>-32.51 ± 20.21i</td></tr><tr><td>47</td><td>11.72 ± 31.19i</td><td>-33.08 ± 19.88i</td><td>-11.97 ± 31.37i</td><td>-33.16 ± 19.87i</td></tr><tr><td>48</td><td>12.07 ± 31.38i</td><td>-33.73 ± 19.52i</td><td>-12.30 ± 31.55i</td><td>-33.80 ± 19.52i</td></tr><tr><td>49</td><td>12.41 ± 31.57i</td><td>-34.37 ± 19.15i</td><td>-12.63 ± 31.74i</td><td>-34.44 ± 19.15i</td></tr><tr><td>50</td><td>12.75 ± 31.76i</td><td>-35.01 ± 18.76i</td><td>-12.96 ± 31.92i</td><td>-35.07 ± 18.76i</td></tr><tr><td>51</td><td>13.08 ± 31.95i</td><td>-35.64 ± 18.35i</td><td>-13.28 ± 32.09i</td><td>-35.71 ± 18.36i</td></tr><tr><td>52</td><td>13.41 ± 32.13i</td><td>-36.27 ± 17.92i</td><td>-13.57 ± 32.27i</td><td>-36.33 ± 17.93i</td></tr><tr><td>53</td><td>13.74 ± 32.31i</td><td>-36.89 ± 17.48i</td><td>-13.91 ± 32.44i</td><td>-36.95 ± 17.49i</td></tr><tr><td>54</td><td>14.06 ± 32.49i</td><td>-37.51 ± 17.01i</td><td>-14.22 ± 32.62i</td><td>-37.57 ± 17.02i</td></tr><tr><td>55</td><td>14.38 ± 32.66i</td><td>-38.13 ± 16.51i</td><td>-14.53 ± 32.79i</td><td>-38.18 ± 16.53i</td></tr><tr><td>56</td><td>14.69 ± 32.84i</td><td>-38.74 ± 15.99i</td><td>-14.84 ± 32.95i</td><td>-38.79 ± 16.01i</td></tr><tr><td>57</td><td>15.00 ± 33.01i</td><td>-39.35 ± 15.44i</td><td>-15.14 ± 33.12i</td><td>-39.40 ± 15.47i</td></tr><tr><td>58</td><td>15.31 ± 33.18i</td><td>-39.95 ± 14.85i</td><td>-15.44 ± 33.28i</td><td>-40.00 ± 14.88i</td></tr><tr><td>59</td><td>15.62 ± 33.34i</td><td>-40.55 ± 14.23i</td><td>-15.74 ± 33.45i</td><td>-40.60 ± 14.26i</td></tr><tr><td>60</td><td>15.92 ± 33.51i</td><td>-41.15 ± 13.56i</td><td>-16.03 ± 33.61i</td><td>-41.20 ± 13.60i</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Illustration of the choice of Q for the controlled airfoil system with LQR and by considering different speeds *V* from 25 m/s to 60 m/s (a) poles (b) estimation of τ_c and ω_c (+ case 1; + case 2; + case 3; + Case 4; o uncontrolled system). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

there is no relationship between the values for the two selected samples (i.e. no similarity between the gain matrices $G^a$ and $G^b$ of the two selected simulations). If $r_{a,b}$ is close to -1 or 1, it indicates the strongest possible agreement. Note that a correlation coefficient $r_{a,b}$ greater than zero means a positive relationship while a value less than zero indicates a negative relationship. Fig. 8 shows the scalar correlation coefficients $r_{a,b}$ for $a, b = 1, \dots, 100$. It appears that the correlation coefficient is very close to 1 whatever simulations are chosen and compared. It means there is a strong similarity between the gain matrices $G^a$ and

$G^b$ that come from two different simulations. However, two groups appear very clearly (see the dark brown and white squares for each $r_{a,b}$). This illustrates the fact that the two categories for the prediction of the critical time delay $\tau_c$ and the angular frequency $\omega_c$ previously seen in Fig. 7 are therefore directly links to the controllers g and f which belong one of the two groups estimate of $G$.

In order to bring elements of explanation on the physical phenomena involved to explain the prediction of two groups of $\tau_c$ so different, and thus to better understand the indirect effect of the choices of g and
<table><caption>**Table 6**
Estimation of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ for the controlled airfoil system with LQR and by considering different speeds $V$ from 25 m/s to 60 m/s.</caption><thead><tr><th>V (m/s)</th><th>Order</th><th>τ<sub>c</sub> (s)</th><th>ω<sub>c</sub> (rad/s)</th><th>[S(λ<sub>i</sub>(τ(ω))−ω)/ω]</th><th>Re(λ<sub>i</sub>(τ<sub>−</sub>(ω)))</th><th>Re(λ<sub>i</sub>(τ<sub>+</sub>(ω)))</th></tr></thead><tbody><tr><td>25</td><td>5</td><td>0.0360</td><td>33.94</td><td>0.00740</td><td>−0.00066</td><td>0.00143</td></tr><tr><td>26</td><td>6</td><td>0.0337</td><td>34.78</td><td>0.00130</td><td>−0.00014</td><td>0.00073</td></tr><tr><td>27</td><td>6</td><td>0.0319</td><td>35.69</td><td>0.00861</td><td>−0.00001</td><td>0.00008</td></tr><tr><td>28</td><td>6</td><td>0.0304</td><td>36.66</td><td>0.00355</td><td>−0.00017</td><td>0.00012</td></tr><tr><td>29</td><td>6</td><td>0.0291</td><td>37.69</td><td>0.00175</td><td>−0.00042</td><td>0.00016</td></tr><tr><td>30</td><td>6</td><td>0.0280</td><td>38.77</td><td>0.00143</td><td>−0.00054</td><td>0.00018</td></tr><tr><td>31</td><td>6</td><td>0.0270</td><td>39.87</td><td>0.00146</td><td>−0.00065</td><td>0.00008</td></tr><tr><td>32</td><td>6</td><td>0.0260</td><td>41.00</td><td>0.00153</td><td>−0.00042</td><td>0.00026</td></tr><tr><td>33</td><td>6</td><td>0.0252</td><td>42.15</td><td>0.00191</td><td>−0.00038</td><td>0.00019</td></tr><tr><td>34</td><td>6</td><td>0.0245</td><td>43.30</td><td>0.00272</td><td>−0.00039</td><td>0.00004</td></tr><tr><td>35</td><td>6</td><td>0.0238</td><td>44.45</td><td>0.00426</td><td>−0.00015</td><td>0.00014</td></tr><tr><td>36</td><td>6</td><td>0.0231</td><td>45.59</td><td>0.00983</td><td>−0.00012</td><td>0.00004</td></tr><tr><td>37</td><td>6</td><td>0.0223</td><td>46.73</td><td>0.00999</td><td>−0.00007</td><td>0.00009</td></tr><tr><td>38</td><td>6</td><td>0.0220</td><td>47.86</td><td>0.00460</td><td>−0.00008</td><td>0.00021</td></tr><tr><td>39</td><td>6</td><td>0.0215</td><td>48.97</td><td>0.00294</td><td>−0.00012</td><td>0.00033</td></tr><tr><td>40</td><td>6</td><td>0.0210</td><td>50.07</td><td>0.00215</td><td>−0.00011</td><td>0.00049</td></tr><tr><td>41</td><td>6</td><td>0.0205</td><td>51.16</td><td>0.00180</td><td>−0.00035</td><td>0.00041</td></tr><tr><td>42</td><td>6</td><td>0.0201</td><td>52.23</td><td>0.00144</td><td>−0.00015</td><td>0.00077</td></tr><tr><td>43</td><td>6</td><td>0.0197</td><td>53.29</td><td>0.00122</td><td>−0.00003</td><td>0.00103</td></tr><tr><td>44</td><td>6</td><td>0.0193</td><td>54.34</td><td>0.00111</td><td>−0.00019</td><td>0.00101</td></tr><tr><td>45</td><td>5</td><td>0.0190</td><td>55.41</td><td>0.00980</td><td>−0.00012</td><td>0.00116</td></tr><tr><td>46</td><td>5</td><td>0.0186</td><td>56.42</td><td>0.00902</td><td>−0.00007</td><td>0.00134</td></tr><tr><td>47</td><td>5</td><td>0.0183</td><td>57.44</td><td>0.00831</td><td>−0.00058</td><td>0.00096</td></tr><tr><td>48</td><td>5</td><td>0.0180</td><td>58.44</td><td>0.00779</td><td>−0.00064</td><td>0.00103</td></tr><tr><td>49</td><td>5</td><td>0.0177</td><td>59.43</td><td>0.00735</td><td>−0.00078</td><td>0.00100</td></tr><tr><td>50</td><td>5</td><td>0.0174</td><td>60.41</td><td>0.00700</td><td>−0.00083</td><td>0.00106</td></tr><tr><td>51</td><td>5</td><td>0.0171</td><td>61.38</td><td>0.00672</td><td>−0.00058</td><td>0.00141</td></tr><tr><td>52</td><td>5</td><td>0.0168</td><td>62.34</td><td>0.00642</td><td>−0.00092</td><td>0.00118</td></tr><tr><td>53</td><td>5</td><td>0.0166</td><td>63.30</td><td>0.00622</td><td>−0.00067</td><td>0.00152</td></tr><tr><td>54</td><td>5</td><td>0.0163</td><td>64.24</td><td>0.00600</td><td>−0.00080</td><td>0.00148</td></tr><tr><td>55</td><td>5</td><td>0.0161</td><td>65.18</td><td>0.00589</td><td>−0.00005</td><td>0.00233</td></tr><tr><td>56</td><td>5</td><td>0.0159</td><td>66.11</td><td>0.00567</td><td>−0.00068</td><td>0.00178</td></tr><tr><td>57</td><td>5</td><td>0.0157</td><td>67.03</td><td>0.00557</td><td>−0.00019</td><td>0.00235</td></tr><tr><td>58</td><td>5</td><td>0.0155</td><td>67.94</td><td>0.00537</td><td>−0.00102</td><td>0.00161</td></tr><tr><td>59</td><td>5</td><td>0.0152</td><td>68.85</td><td>0.00529</td><td>−0.00048</td><td>0.00221</td></tr><tr><td>60</td><td>5</td><td>0.0150</td><td>69.76</td><td>0.00513</td><td>−0.00117</td><td>0.00160</td></tr></tbody></table>

<table><caption>**Table 7**
Values of Q and R for the four cases under study.</caption><thead><tr><th>Case</th><th>R</th><th>Q</th></tr></thead><tbody><tr><td>1</td><td>I<sub>2x2</sub></td><td>diag(1/0.012, 1/0.012, 0, 0)</td></tr><tr><td>2</td><td>I<sub>2x2</sub></td><td>diag(1/0.12, 1/0.12, 0, 0)</td></tr><tr><td>3</td><td>I<sub>2x2</sub></td><td>diag(1/0.0012, 1/0.0012, 0, 0)</td></tr><tr><td>4</td><td>I<sub>2x2</sub></td><td>diag(1/0.012, 1/0.012, 1/0.12, 1/0.12)</td></tr></tbody></table>

f, we choose to study one case for each group in more detail. The first one denoted by 30+ corresponds to the case for which $\Delta\omega_c > 0$ and the second one denoted by 30− for the case $\Delta\omega_c < 0$. We recall that in each case we consider the controlled airfoil system for $V = 30$ m/s. Table 8 provides the matrices of proportional velocity and displacement control gains f and g for one specific simulation result of both cases 30+ and 30−. Figs. 9 give the evolution of the imaginary parts of the primary and secondary eigenvalues versus $\tau$ based on the Taylor expansion of $p$th order (for $p = 3, ..., 6$) as well as the delay margin $\tau(\omega)$ based on the frequency sweeping approach for both cases 30+ and 30−. Additionally Table 8 gives the prediction of the critical time delay $\tau_c$ based on the proposed methodology previously described in Section 3.2.3 for different orders of the Taylor expansion. Moreover the estimation of the quality of the prediction (i.e. $\frac{\int(\lambda_i(\tau(\omega)) - \omega)}{\omega}$) and the values of $Re(\lambda_i(\tau_-(\omega)))$ and $Re(\lambda_i(\tau_+(\omega)))$ are also given. Showing Figs. 9 it can be noted that the evolutions of the imaginary part of the unstable mode are opposite for cases 30+ and 30−. For case 30+ an increase of the imaginary part with the time delay $\tau$ is observed, while a decrease for case 30− is shown with a crossing between modes I and II at $\tau \approx 0.002$ s. For more clarity evolution of both mode I (i.e. the mode of the controlled airfoil system that becomes unstable for $\tau > \tau_c$) and

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Prediction of the critical time delay $\tau_c$ and $\Delta\omega_c$ for the pole placement technique and by considering the variability $(1+\epsilon)\lambda_{1,2}^{\text{controlled}}$ with $\epsilon \in [-1; 1] \times 10^{-8}$.</figcaption></figure>

mode II versus the time delay is shown in Figs. 10 for cases 30+ and 30−. Without any doubt it is demonstrated here that the two values of the predicted critical time delay $\tau_c$ for cases 30− and 30+ are indeed related to two different dynamic behaviors of the controlled airfoil system and more particularly to a different evolution of mode I versus the time delay $\tau$, this being directly related to the choice of the control gains f and g.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Correlation coefficient $r_{a,b}$ for the gain matrices G based on 100 simulations for the controlled airfoil system with the pole placement technique at $V = 30$ m/s.</figcaption></figure>

<table><caption>**Table 8**
Estimation of f and g for the two cases $30_+$ and $30_-$ and by considering the controlled airfoil system with the pole placement technique.</caption><thead><tr><th>Case</th><th>f</th><th>g</th></tr></thead><tbody><tr><td>30<sub>+</sub></td><td>[-5.35 29.68 ; -5.83 22.95]</td><td>[-689.93 2058.48 ; -62.97 1023.68]</td></tr><tr><td>30<sub>-</sub></td><td>[-3.32 39.13 ; -6.56 20.92]</td><td>[-696.55 2978.33 ; -913.60 1003.73]</td></tr></tbody></table>

In addition, it is worth noting that all the results presented in Table 9 and in Figs. 9 (via the crossing between the evolution of imaginary parts and the delay margin $\tau(\omega)$) demonstrate the effectiveness and relevance of the proposed approach to predict the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$. As a complement Figs. 11 give the temporal vibration responses just before and just after each critical time delay $\tau_c$ of both cases $30_-$ and $30_+$. The angular frequency of the unstable mode is estimated to be approximately equal to 36.7 rad/s and 23.4 rad/s for case $30_+$ and $30_-$, respectively. These results are in perfect agreement with the previous analysis. A full discussion of the validity of the proposed method for predicting $\tau_c$ and $\omega_c$ from the delay margin $\tau(\omega)$ as well as a detailed investigation of the results provided in Table 9 and in Figs. 9 and 11 are not proposed in this section since the argument would be quite similar and redundant with the explanation already given in Section 4.1. Note that these results are not in disagreement with the fact that the algorithm used (Kautsky et al., 1985) finds a robust solution, such that the assigned closed-loop poles are insensitive to perturbations in the system data (i.e. the provided solution for the controlled system without time delay minimizes the sensitivity of the closed-loop poles to perturbations in the input matrices of the airfoil system). In this present case the algorithm is not robust against the occurrence of a flutter instability for the controlled system with time delay since small perturbations in the assigned poles can induce different values of the critical time delay and, consequently, an evolution of the flutter instability and the dynamic behavior of controlled airfoil system with time delay.

Then the dynamic behavior of the controlled airfoil system in relation to the evolution of the speed $V$ from 25 m/s to 35 m/s is briefly

<table><caption>**Table 9**
Estimation of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ for the two cases $30_-$ and $30_+$ and by considering the controlled airfoil system with the pole placement technique.</caption><thead><tr><th>Case</th><th>Order</th><th>τ<sub>c</sub> (s)</th><th>ω<sub>c</sub> (rad/s)</th><th>|R(λ<sub>i</sub>(τ(ω)) - ω|<br>ω</th><th>R(λ<sub>i</sub>(τ<sub>-</sub>(ω)))</th><th>R(λ<sub>i</sub>(τ<sub>+</sub>(ω)))</th></tr></thead><tbody><tr><td>30<sub>+</sub></td><td>3</td><td>0.0173</td><td>36.58</td><td>0.00488</td><td>-0.00421</td><td>0.00221</td></tr><tr><td>30<sub>+</sub></td><td>4</td><td>0.0174</td><td>36.67</td><td>0.00109</td><td>-0.00428</td><td>0.00180</td></tr><tr><td>30<sub>+</sub></td><td>5</td><td>0.0174</td><td>36.68</td><td>0.00006</td><td>-0.00466</td><td>0.00148</td></tr><tr><td>30<sub>+</sub></td><td>6</td><td>0.0174</td><td>36.68</td><td>0.00001</td><td>-0.00380</td><td>0.00001</td></tr><tr><td>30<sub>-</sub></td><td>3</td><td>0.0075</td><td>23.45</td><td>0.00001</td><td>-0.00493</td><td>0.00612</td></tr><tr><td>30<sub>-</sub></td><td>4</td><td>0.0075</td><td>23.45</td><td>0.00005</td><td>-0.00427</td><td>0.00678</td></tr><tr><td>30<sub>-</sub></td><td>5</td><td>0.0075</td><td>23.45</td><td>0.00005</td><td>-0.00424</td><td>0.00681</td></tr><tr><td>30<sub>-</sub></td><td>6</td><td>0.0075</td><td>23.45</td><td>0.00005</td><td>-0.00424</td><td>0.00681</td></tr></tbody></table>

undertaken to confirm the analysis and results previously proposed on the use of the pole placement technique. Fig. 12(a) gives the evolution of the time delay $\tau_c$ and the associated angular frequency $\omega_c$ in relation to the evolution of the speed $V$ from 25 m/s to 35 m/s with a step of 1 m/s. For a given speed two categories, denoted by the symbols o and □ in Fig. 12(a), are observed for the prediction of the critical time delay $\tau_c$ and the angular frequency $\omega_c$. This result is perfectly consistent with the previous observations made for the case $V = 30$ m/s. Note that the figure summarizes all the results obtained for each speed between [25;35] m/s: for each speed, several simulations have been performed and only one result per category is shown. Table 11 gives the prediction of the critical time delay $\tau_c$ and the angular frequency of the unstable mode $\omega_c$ based on the proposed methodology and by using a 5th order Taylor expansion, as well as the estimation of the quality of the prediction (i.e. $|R(\lambda_i(\tau(\omega)) - \omega|/\omega$) and the values of $R(\lambda_i(\tau_-(\omega)))$ and $R(\lambda_i(\tau_+(\omega)))$. Table 10 provides the associated controllers g and f. For the reader comprehension, the two categories are denoted by $V_-$ and $V_+$ (with $V$ evolving between 25 m/s and 35 m/s with a step of 1 m/s) respecting the denomination previously chosen in the study for the specification of the indices $-$ and $+$. As previously explained for the case $V = 30$ m/s, the two dynamic behaviors of the controlled airfoil system
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Evolution of $(\tau, \omega)$, for the controlled airfoil system with the pole placement technique, based on the Taylor expansion of $p$th order for (a) case 30$_+$ and (b) case 30$_-$ (- - - 3rd order; --- 4th order; - - - 5th order; --- 6th order; + prediction of $(\tau_c, \omega_c)$; — $\tau(\omega)$ based on the frequency sweeping approach; ◦ angular frequencies for the controlled airfoil system without time delay). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Evolution of (a) mode I and (b) mode II versus the time delay $\tau$ (green = stable; red = unstable) for the controlled airfoil system with the pole placement technique (case 30$_+$ denoted by ○ and case 30$_-$ denoted by □). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

versus the critical time delay are linked to the choice of the control gains f and g. Fig. 12(b) illustrates this fact by showing the coefficient of similarity between gain matrices G from different simulations in relation to the evolution of the speed V from 25 m/s to 35 m/s with a step of 1 m/s (i.e. $V_-$ and $V_+$ with $V = [25 : 1 : 35]$ m/s). We see very clearly the two categories $V_-$ and $V_+$ with a perfect correlation on all the speeds for each category. These results illustrate therefore all the precautions to be taken when using a controller based on the pole placement technique. Indeed, different dynamic behaviors of the controlled airfoil system can appear potentially leading to a more or less early appearance of flutter instability in relation to the time delay. Moreover this may induce a strong variation of the flutter frequency $\omega_c$ which can be prejudicial in the context of designing a airfoil system with active control.

Finally Fig. 13 illustrates the prediction of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ of the unstable mode for the controlled airfoil system with the pole placement technique and by considering three values for the control gain applied on the real part of the uncontrolled poles (i.e. $\gamma = -0.5$, $\gamma = -1$ and $\gamma = -2$; see Eq. (16) for more details). To be noted that the value of $\gamma$ previously used in the study for the pole placement technique is equal to -1.

This parametric study is performed by considering 15000 simulations for different speeds V from 25 m/s to 60 m/s. Simulation results unsurprisingly illustrate the complexity associated with the potential evolution of $(\tau_c, \omega_c)$ and demonstrate all the precautions to be taken for the suppression of flutter instability by using active feedback control. Firstly it can be noted that the control gain $\gamma$ has an expected effect with an overall shift of $(\tau_c, \omega_c)$ to the right by increasing the absolute value of $\gamma$. However, it is important to keep in mind the results presented earlier, which showed that the phenomena involved are complex, so an increase in the absolute value of $\gamma$ (for a given value of the speed V) will not necessarily result in an increase in the critical time delay $\tau_c$. Secondly, the two categories previously denoted by $V_-$ and $V_+$ can also be guessed on the set of results associated with each control gain $\gamma$. For results that belong to group $V_+$, a significant variation in the angular frequency $\omega_c$ is also to be reported. This implies that the controller has a significant impact not only on the value of the critical time delay but also that the flutter frequency can drastically be modified by the controller due to unavoidable time delays (in comparison to the flutter frequency of the uncontrolled or controlled system without time delay). This can be problematic for engineers in the design process of such an airfoil system with regard to the problem of flutter instability and
<figure><img src="image_12.png" /><figcaption>Fig. 11. Time responses of the controlled airfoil system with the pole placement technique for (a,b) case 30+ and (c,d) case 30− (a,c) before the critical time delay with $τ_c (1 - χ)$ and (b,d) after the critical time delay with $τ_c (1 + χ)$ for $χ = 0.01$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Prediction of (a) the critical time delay $τ_c$ and the associated angular frequency $ω_c$ and (b) the correlation coefficient for the controlled airfoil system with the pole placement technique and by considering different speeds $V$ from 25 m/s to 35 m/s.</figcaption></figure>

the associated self-excited vibrations, and must therefore be taken into account and looked at carefully when designing a controller to try suppressing flutter instability. In some cases, the value of the critical

time delay $τ_c$ is very low for the three values of $γ$ (see all the couples $(τ_c, ω_c)$ on the left side of Fig. 13). Due to the fact that such short time delays in control loop are unavoidable in practice, particular
<table><caption>**Table 10**
Estimation of f and g for the controlled airfoil system with the pole placement technique and by considering different speeds V the from 25 m/s to 35 m/s.</caption><thead><tr><th>Case</th><th>f</th><th>g</th></tr></thead><tbody><tr><td>25</td><td>[-8.81 34.24; -1.95 13.61]</td><td>[-777.44 2902.02; -687.91 803.53]</td></tr><tr><td>25+</td><td>[-9.14 34.90; -4.42 13.94]</td><td>[-777.25 2153.15; 51.46 803.87]</td></tr><tr><td>26</td><td>[-7.58 35.50; -2.90 15.26]</td><td>[-761.51 2917.93; -729.30 841.87]</td></tr><tr><td>26+</td><td>[-8.66 32.47; -5.75 16.34]</td><td>[-759.21 2124.62; 15.46 848.11]</td></tr><tr><td>27</td><td>[-5.33 30.94; -9.23 15.71]</td><td>[-729.12 2855.42; -792.23 900.23]</td></tr><tr><td>27+</td><td>[-7.82 31.39; -5.96 18.19]</td><td>[-741.55 2106.99; -6.70 891.52]</td></tr><tr><td>28</td><td>[-5.29 37.37; -4.89 18.20]</td><td>[-729.11 2946.98; -819.26 921.17]</td></tr><tr><td>28+</td><td>[-5.77 33.82; -2.65 18.68]</td><td>[-729.62 2104.92; 17.57 922.07]</td></tr><tr><td>29</td><td>[-4.39 38.74; -5.28 19.70]</td><td>[-713.14 2966.82; -862.30 961.79]</td></tr><tr><td>29+</td><td>[-4.74 33.78; -2.00 20.05]</td><td>[-713.61 2089.43; 6.68 962.74]</td></tr><tr><td>30</td><td>[-3.52 39.95; -5.76 21.12]</td><td>[-697.15 2985.73; -907.58 1003.42]</td></tr><tr><td>30+</td><td>[-3.61 34.23; -0.93 21.21]</td><td>[-697.33 2077.09; -0.93 1003.60]</td></tr><tr><td>31</td><td>[-1.85 37.06; -10.05 21.63]</td><td>[-672.62 2960.27; -974.03 1053.41]</td></tr><tr><td>31+</td><td>[-2.75 34.10; -0.49 22.53]</td><td>[-681.31 2059.68; -14.85 1046.04]</td></tr><tr><td>32</td><td>[-1.73 41.48; -7.45 23.61]</td><td>[-664.84 3017.40; -1008.67 1089.00]</td></tr><tr><td>32+</td><td>[-2.26 33.16; -0.93 24.14]</td><td>[-665.96 2035.75; -37.29 1090.98]</td></tr><tr><td>33</td><td>[-2.37 42.36; -8.14 26.25]</td><td>[-658.83 3031.81; -1063.90 1142.53]</td></tr><tr><td>33+</td><td>[-1.27 36.92; 3.31 25.15]</td><td>[-631.13 2054.74; -22.64 1121.58]</td></tr><tr><td>34</td><td>[-1.62 43.37; -8.68 27.44]</td><td>[-643.50 3050.54; -1116.07 1187.61]</td></tr><tr><td>34+</td><td>[-1.50 34.78; 1.69 27.32]</td><td>[-640.44 2016.93; -49.10 1184.66]</td></tr><tr><td>35</td><td>[-0.90 44.39; -9.19 28.61]</td><td>[-628.51 3069.75; -1169.47 1233.64]</td></tr><tr><td>35+</td><td>[0.18 36.40; 3.73 27.52]</td><td>[-600.10 2013.58; -56.30 1210.22]</td></tr></tbody></table>

<table><caption>**Table 11**
Estimation of the critical time delay $\tau_c$ and the associated angular frequency $\omega_c$ for the controlled airfoil system with the pole placement technique and by considering different speeds V the from 25 m/s to 35 m/s.</caption><thead><tr><th>Case</th><th>τ<sub>c</sub> (s)</th><th>ω<sub>c</sub> (rad/s)</th><th>[∫(λ<sub>i</sub>(τ<sub>c</sub>(ω))·ω)]/ω</th><th>∫(λ<sub>i</sub>(τ<sub>c</sub>(ω)))</th><th>∫(λ<sub>i</sub>(τ<sub>+</sub>(ω)))</th></tr></thead><tbody><tr><td>25</td><td>0.0030</td><td>23.37</td><td>0.00001</td><td>-0.00517</td><td>0.00512</td></tr><tr><td>25+</td><td>0.0119</td><td>31.50</td><td>0.00357</td><td>-0.59784</td><td>0.13228</td></tr><tr><td>26</td><td>0.0040</td><td>23.33</td><td>0.00004</td><td>-0.00410</td><td>0.00659</td></tr><tr><td>26+</td><td>0.0130</td><td>32.16</td><td>0.00005</td><td>-0.00581</td><td>0.00365</td></tr><tr><td>27</td><td>0.0041</td><td>24.20</td><td>0.00005</td><td>-0.00779</td><td>0.00069</td></tr><tr><td>27+</td><td>0.0143</td><td>33.22</td><td>0.00001</td><td>-0.00542</td><td>0.00555</td></tr><tr><td>28</td><td>0.0061</td><td>23.27</td><td>0.00004</td><td>-0.00717</td><td>0.00014</td></tr><tr><td>28+</td><td>0.0163</td><td>35.78</td><td>0.00010</td><td>-0.00286</td><td>0.00810</td></tr><tr><td>29</td><td>0.0070</td><td>23.23</td><td>0.00011</td><td>-0.00557</td><td>0.00093</td></tr><tr><td>29+</td><td>0.0171</td><td>37.19</td><td>0.00001</td><td>-0.00014</td><td>0.01085</td></tr><tr><td>30</td><td>0.0078</td><td>23.25</td><td>0.00005</td><td>-0.00125</td><td>0.00976</td></tr><tr><td>30+</td><td>0.0178</td><td>38.78</td><td>0.00010</td><td>-0.00359</td><td>0.00209</td></tr><tr><td>31</td><td>0.0074</td><td>24.19</td><td>0.00022</td><td>-0.00300</td><td>0.00201</td></tr><tr><td>31+</td><td>0.0182</td><td>40.09</td><td>0.00010</td><td>-0.00332</td><td>0.00773</td></tr><tr><td>32</td><td>0.0088</td><td>23.60</td><td>0.00006</td><td>-0.00404</td><td>0.00703</td></tr><tr><td>32+</td><td>0.0185</td><td>41.02</td><td>0.00030</td><td>-0.00305</td><td>0.00147</td></tr><tr><td>33</td><td>0.0095</td><td>23.53</td><td>0.00005</td><td>-0.00134</td><td>0.00973</td></tr><tr><td>33+</td><td>0.0191</td><td>44.55</td><td>0.00157</td><td>-0.00192</td><td>0.00066</td></tr><tr><td>34</td><td>0.0099</td><td>23.71</td><td>0.00005</td><td>-0.00131</td><td>0.00977</td></tr><tr><td>34+</td><td>0.0188</td><td>45.03</td><td>0.00153</td><td>-0.00174</td><td>0.00086</td></tr><tr><td>35</td><td>0.0102</td><td>23.90</td><td>0.00367</td><td>-0.00123</td><td>0.00034</td></tr><tr><td>35+</td><td>0.0190</td><td>46.77</td><td>0.00010</td><td>-0.00294</td><td>0.00814</td></tr></tbody></table>

attention should be paid to the choice of the gain matrix. An alternative advanced strategy for active aeroelastic control may be used in these specific cases. To avoid the sensitivity of the controllers to time-delay a different strategy for active aeroelastic control with time delay might be to define an effective controller for a given constant time delay that represents the known offset between the time required to compute gains/forces for a desired control objective based on sensor signals and physical actuators to satisfy the control requirements. The development of such a control strategy is beyond the scope of this work, the reader can refer to Singh's article (Singh, 2015) for a complete study on the subject and the problem addressed in this article.

## 5. Conclusion

In this work, the efficiency of a hybrid methodology, based on the frequency sweep test and the eigenvalue problem approximation using the Taylor series expansion of the delayed term, for flutter instability of

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Prediction of $(\tau_c, \omega_c)$ for various control gains (black: $\gamma = -0.5$; red: $\gamma = -1$; blue: $\gamma = -2$) with the pole placement technique and by considering different speeds V from V = 25 m/s to 60 m/s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

airfoil controlled systems is illustrated. This computational technique is more specifically tested to predict the critical time delay of a two-dimensional airfoil controlled system with time delay as well as the associated flutter angular frequency. Moreover the use of the hybrid methodology allows not only an efficient and fast prediction of the critical time delay but also an estimation of the different dynamic behaviors that can occur during the appearance of a flutter instability for a controlled system with time delay. It is illustrated without any ambiguity that the presence of time delay in controllers may induce detrimental degradation of the controlled airfoil system with the appearance of flutter instability. Finally, a discussion is conducted on the dependence of the critical time delay versus some specific parameters of the LQR state-feedback controller and the pole placement technique.

Some non-exhaustive interesting outlooks for future research are as follows:

* One of the main challenges for engineers working on active control and flutter is to provide practical applications of active control to suppress or mitigate flutter instability phenomena in the real world. The study proposed in this paper focuses only on a simplified theoretical model (i.e. a two-dimensional airfoil system), in order to illustrate clearly the use of the proposed hybrid approach for predicting the critical time delay and to illustrate the potential limitations or recommendations when applying well-known linear controllers such as the LQR approach or the pole placement. This study does not address the question of how control forces can be applied to airfoil systems in the real world. Classically, the control forces are most often generated by the deflection of the trailing edge control surface. Therefore, it would be necessary to develop in future research more complex airfoil systems such as the well-known three-dimensional airfoil model with trailing edge control surface to address this issue. Moreover the issue of the development and evaluation of control laws in practical cases is an important task. Indeed control laws can be described in different ways such as ignoring actuator dynamics, introducing actuator dynamics with a structural filter or using only actuator dynamics. It is well-known that actuator dynamics can cause phase lag in the closed-loop.

* A directly related point concerns the control robustness that is achieved in practical cases at a cost of increased control surface activity for the airfoil system. Robust control theory relates more specifically to the view of controller design that explicitly deals
with uncertainty. The question is then to be able to design robust controls to operate correctly even if uncertain parameters or disturbances are considered. To be noted that defects in the geometry of the mechanical system under study or dispersion in the material properties are inherent in the manufacturing process, which therefore brings possible variability with respect to a deterministic target reference. In addition, a system is often subjected to changes due to degradation, thermal effects during service, or other factors that are difficult to identify precisely. These variabilities can have an impact on the key characteristics of controlled airfoil systems, causing deviations in performance from an optimal initial design, and thus unexpected flutter and vibration deterioration. Therefore the introduction of uncertainties and the propagation and analysis of the associated results of the controlled airfoil system with time delay are subjects of primary importance to adequately address the real problems faced by engineers and their need to accurately and robustly predict flutter instability.

* The proposed study focuses on flutter instability and active aeroelastic control with time delay for a airfoil system subjected to quasi-steady aerodynamic loading. The extension to the case of an unsteady aerodynamic model should be interesting for active flutter suppression from a more general point of view.

* The proposed study investigates only the notion of the onset of the flutter instability (i.e. the study is limited to the linear region). The extension to the study for estimating the nonlinear responses of airfoil systems and the effect of the time delay in control design on the vibrational magnitudes would be an interesting topic. This would naturally require the consideration of non-linearities in structural stiffness or aerodynamics.

## Declaration of competing interest

The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgment

J.-J. Sinou acknowledges the support of the Institut Universitaire de France. The author approved the version of the manuscript to be published.

## References

Bhoi, N., Singh, S.N., 2005. Control of unsteady aeroelastic system via state dependent riccati equation method. J. Guid. Control Dyn. 28, 78-84.

Borglund, D., Kuttenkeuler, J., 2002. Active wing flutter suppression using a trailing edge flap. J. Fluids Struct. 16 (3), 271-294.

Brunton, S.L., Rowley, C.W., 2013. Empirical state-space representations for Theodorsen's lift model. J. Fluids Struct. 38, 174-186.

Gu, K., Kharitonov, V.L., Chen, J., 2003. Stability of Time-Delay Systems. Birhauser Publication, Boston.

Kautsky, J., Nichols, N.K., Vandooren, P., 1985. Robust pole assignment in linear state feedback. Internat. J. Control 41 (5), 1129-1155.

Ko, J., Strganac, T.W., Kurdila, A.J., 1999. Adaptive feedback linearization for the control of a typical wing section with structural nonlinearity. Nonlinear Dynam. 18, 289-301.

Kwakernaak, H., Sivan, R., 1972. Linear Optimal Control Systems. Wiley-Interscience, John Wiley & Sons, pp. 1-924.

Lee, Y.S., Vakakis, A., Bergman, L., McFarland, D.M., Kerschen, G., 2007. Suppressing aeroelastic instability using broadband passive targeted energy transfers, part 1: theory. AIAA J. 45, 693-711.

Librescu, L., Marzocca, P., 2005. Advances in the linear/nonlinear control of aeroelastic structural systems. Acta Mech. 178, 147-186.

Librescu, L., Marzocca, P., Silva, W.A., 2005. Aeroelasticity of 2-D lifting surfaces with time-delayed feedback control. J. Fluids Struct. 20, 197-215.

Luo, M., Gao, M., Cai, G., 2016. Delayed full-state feedback control of airfoil flutter using sliding mode control method. J. Fluids Struct. 61, 262-273.

Malher, A., Touze, C., Doare, O., Habib, G., Kerschen, G., 2017. Flutter control of a two-degrees-of-freedom airfoil using a nonlinear tuned vibration absorber. J. Comput. Nonlinear Dyn. Am. Soc. Mech. Eng. 12, 051016.

Ouyang, Y., Gu, Y., Kou, X., Yang, Z., 2021. Active flutter suppression of wing with morphing flap. Aerosp. Sci. Technol. 110, 106457.

Platanitis, G., Strganac, T.W., 2004. Control of a nonlinear wing section using leading and trailing edge surfaces. J. Guid. Control Dyn. 27, 52-58.

Ram, Y.M., Mottershead, J.E., 2007. The receptance method in active vibration control. AIAA J. 45, 562-567.

Ram, Y.M., Singh, A.N., Mottershead, J.E., 2009. State feedback control with time delay. Mech. Syst. Signal Process. 23, 1940-1945.

Ramesh, M., Narayanan, S., 2001. Controlling chaotic motions in a two-dimensional airfoil using time-delayed feedback. J. Guid. Control Dyn. 239, 1037-1049.

Singh, K.V., 2015. Active aeroelastic control with time delay for targeted flutter modes. Aerosp. Sci. Technol. 43, 281-288.

Sinou, J-J., Chomette, B., 2021. Active vibration control and stability analysis of a time-delay system subjected to friction-induced vibration. J. Sound Vib. 500, 116013.

Theodorsen, T., 1935. General Theory of Aerodynamic Instability and the Mechanism of Flutter. NACA TR 496, pp. 291-311.

Tsushima, N., Su, W., 2017. Flutter suppression for highly flexible wings using passive and active piezoelectric effects. Aerosp. Sci. Technol. 65, 78-89.

Valasek, M., Olgac, N., 1995. Efficient eigenvalue assignments for general linear MIMO systems. Automatica 31, 1605-1617.

Yuan, Y., Yu, P., Librescu, L., 2004. Aeroelasticity of time-delayed feedback control of two-dimensional supersonic lifting surfaces. J. Guid. Control Dyn. 27, 795-803.

Zhao, Y.H., 2009. Stability of a two-dimensional airfoil with time-delayed feedback control. J. Fluids Struct. 25, 1-25.

Zhao, Y.H., 2011. Stability of a time-delayed aeroelastic system with a control surface. Aerosp. Sci. Technol. 15, 72-77.