2.1. Model reduction for aeroelastic problem

POD is a method to rebuild the behaviour of the overall system with a small number of degrees of freedom as in the case of unsteady flows for aeroelastic analysis [Reference Hall, Thomas and Dowell37, Reference Lieu, Farhat and Lesoinne38]. Here, POD is employed to construct a CFD-based ROM for Euler equations in active flutter suppression. Compared with the traditional aeroelastic linearised flow equations, the aerodynamic force is jointly determined by the structural deformation and the deflection of control surface. First, the unsteady Euler equations are linearised around a mean flow solution (equilibrium), with the flow solution satisfying the relation $\textbf{F}({\textbf{w}_0},{\rm{ }}{\textbf{u}_0},{\rm{ }}{\textbf{u}_{c0}})$ , where ${\textbf{w}_0}\,=\,{\boldsymbol{0}}$ , $\dot{\textbf{w}}_0\,=\,{\boldsymbol{0}}$ , $\textbf{u}_0\,=\,{\boldsymbol{0}}$ , $\dot{\textbf{u}}_0\,=\,{\boldsymbol{0}}$ , $\textbf{u}_{c0}\,=\,{\boldsymbol{0}}$ , $\dot{\textbf{u}}_{c0}\,=\,{\boldsymbol{0}}$ . Here, F indicates the nonlinear numerical flux resulting from the spatial discretisation, w is the vector of conservative flow variables, u is the structural displacement vector, $\textbf{u}_c$ denotes the control surface deflection. Supposing that $\partial \textbf{w},{\rm{ }}\partial \textbf{u},{\rm{ }}\partial \textbf{u}_c$ are small perturbations around the steady state variables. Then, one obtains the linearised flow equations about the steady-state solution $(\textbf{w}_0,{\rm{ }}\textbf{u}_0,{\rm{ }}\textbf{u}_{c0})$ through a first order Taylor’s expansion:

(1) \begin{align}\textbf{A}_{\boldsymbol{0}}\dot {\textbf w} + \textbf{Hw} + (\textbf{E} + \textbf{C})\dot {\textbf u} + \textbf{Gu} + \textbf{Lu}_c = \boldsymbol{0}\end{align}

where

\begin{align*}\textbf H = \frac{{\partial \textbf F}}{{\partial \textbf w}}({\textbf w_0},{\textbf u_0},{\textbf u_{c0}})\quad \textbf G = \frac{{\partial \textbf F}}{{\partial \textbf u}}({\textbf w_0},{\textbf u_0},{\textbf u_{c0}})\quad \textbf E = \frac{{\partial \textbf A}}{{\partial \textbf u}}{\textbf w_0}\quad \textbf C = \frac{{\partial \textbf F}}{{\partial \dot {\textbf u}}}({\textbf w_0},{\textbf u_0},{\textbf u_{c0}})\quad \textbf L = \frac{{\partial \textbf F}}{{\partial {\textbf u_c}}}({\textbf w_0},{\textbf u_0},{\textbf u_{c0}})\end{align*}

Here, A denotes the diagonal matrix of cell volumes. The matrices H, E, G, C and L are the first-order terms of a Taylor expansion of the non-dimensional numerical flux function around the equilibrium ( ${\textbf w_0},{\textbf u_0},{\textbf u_{c0}}$ ). The matrix H is the gradient of the numerical flux function with respect to the vector of fluid variables. The matrix E indicates the gradient of the cell volumes with respect to the generalised coordinates. Matrices G and C are the gradients of the flux function with respect to the generalised coordinates and their velocities, respectively. Finally, the matrix L indicates the gradient of the numerical flux function with respect to control surface deflection. Note that the matrices G, E and C need to be recomputed for modified structure.

In active flutter suppression system, the elastic influence of the control surface can be ignored, so the linearisation of the structural dynamic equation around an equilibrium state can be written as:

(2) \begin{align}\bar {\textbf M}\ddot {\textbf u} + {\bar {\textbf D}_0}\dot {\textbf u} + {\bar {\textbf K}_s}\textbf u = {\textbf P_0} \textbf w\end{align}

where

\begin{align*}{\bar {\textbf K}_0} = \frac{{\partial {\textbf f^{int}}}}{{\partial \textbf u}}({\textbf u_0},{\dot {\textbf u}_0})\quad {\bar {\textbf K}_s} = {\bar {\textbf K}_0} - \frac{{\partial {\textbf f^{ext}}}}{{\partial \textbf u}}({\textbf w_0},{\textbf u_0})\quad {\bar {\textbf D}_0} = \frac{{\partial {\textbf f^{int}}}}{{\partial \dot {\textbf u}}}({\textbf u_0},{\dot {\textbf u}_0})\quad {\textbf P_0} = \frac{{\partial {\textbf f^{ext}}}}{{\partial \textbf w}}({\textbf w_0},{\textbf u_0})\end{align*}

For the analysis of the system stability, the terms $\dfrac{{\partial {\textbf f^{ext}}}}{{\partial \textbf u}}$ and ${\bar {\textbf D}_0}$ can be neglected, which leads to the structural dynamic equation:

(3) \begin{align}\bar {\textbf M}\ddot {\textbf u} + {\bar {\textbf K}_0} \textbf u = {\textbf P_0} \textbf w\end{align}

The linearised CFD-based aeroelastic equations of Eq. (1) and Eq. (3) are too large for aeroelastic system. Therefore, the POD is used to reduce the full order models (FOMs). Denote $\{ {\textbf x^k}\} $ , $k = 1,2,3 \ldots m$ , a set of data, with ${\textbf x^k}$ is the n-dimensional space, and m is the number of snapshots. The POD method searches an m-dimensional proper orthogonal subspace, ${\boldsymbol \Psi} \in {\textbf R^{n \times m}}$ , to minimise the mapping errors from ${\boldsymbol \Psi}$

(4) \begin{align} \textbf G = \mathop {\min }\limits_{\boldsymbol \Phi} \sum\limits_{k = 1}^m {\left\| {{\textbf x^k} - {\boldsymbol \Omega} {\boldsymbol \Omega} {^T}{\textbf x^k}} \right\|} = \sum\limits_{k = 1}^m {\left\| {{\textbf x^k} - {\boldsymbol \Psi} {{\boldsymbol \Psi} ^T}{\textbf x^k}} \right\|} ,\quad {{\boldsymbol \Omega} ^H}{\boldsymbol \Omega} = {\textbf I} \end{align}

The minimisation problem is equivalent to [Reference Lieu39]:

(5) \begin{align} \textbf H = \mathop {\max }\limits_{\boldsymbol \Phi} \sum\limits_{k = 1}^m {\frac{{\left\langle {{{\left( {{\textbf x^k},{\boldsymbol \Omega} } \right)}^2}} \right\rangle }}{{{{\left\| {\boldsymbol \Omega} \right\|}^2}}}} = \sum\limits_{k = 1}^m {\frac{{\left\langle {{{\left( {{\textbf x^k},{\boldsymbol \Psi} } \right)}^2}} \right\rangle }}{{{{\left\| {\boldsymbol \Psi} \right\|}^2}}}} ,\quad {{\boldsymbol \Omega} ^H}{\boldsymbol \Omega} = \textbf I\end{align}

The constraint optimisation problem in Eq. (5) is transformed into the following Lagrange equation

(6) \begin{align} J\left( {\boldsymbol \Omega} \right) = \sum\limits_{k = 1}^m {{{\left( {{\textbf x^k},{\boldsymbol \Omega} } \right)}^2}} - \lambda \left( {{{\left\| {\boldsymbol \Omega} \right\|}^2} - 1} \right)\end{align}

Solving the partial derivative of the objective function $J(\boldsymbol \Omega )$ with respect to $\boldsymbol \Omega $ gives

(7) \begin{align} \frac{d}{{d{\boldsymbol \Omega} }}J\left( {\boldsymbol\Omega} \right) = 2{\textbf X}{\textbf X^H}{\boldsymbol \Omega} - 2\lambda {\boldsymbol \Omega} \end{align}

where $\textbf X = \{ {\textbf x_1},{\textbf x_2}, \cdots {\textbf x_m}\} \in {\textbf R^{n \times m}}$ is a matrix containing m snapshots as columns. By setting Eq. (7) to zero; thus, the following equation is obtained:

(8) \begin{align} \left( {\textbf X{\textbf X^T} - \xi \textbf I} \right){\boldsymbol \Psi} = 0\end{align}

The problem is transformed into solving the eigenvalue problem of the POD kernel, $\textbf X{\textbf X^H}$ . The eigenvalue problem has very large size, as $\textbf X{\textbf X^H} \in {\textbf R^{n \times n}}$ . Because $\textbf X{\textbf X^H}$ and $\textbf X{\textbf X^H}$ have the same eigenvalues, so we can obtain ${\boldsymbol \Psi} $ as:

(9) \begin{align} \left\{ \begin{array}{l}\textbf X{\textbf X^T}\textbf V = \textbf V{\boldsymbol\Lambda} \\[5pt] {\boldsymbol \Psi} = {\textbf {XV}}{\boldsymbol \Lambda ^{ - 1/2}}\end{array} \right.\end{align}

where ${\boldsymbol \Psi} = [{{\boldsymbol \psi} _1},{{\boldsymbol \psi} _2}, \ldots ,{{\boldsymbol \psi} _m}]$ , ${\boldsymbol \Lambda} = diag({\xi _1},{\xi _2} \cdots {\xi _m})$ , ${\xi _1} \ge {\xi _2} \ge ,\ldots, \ge {\xi _m}$ . The value of ${\xi _i}$ represents the contribution of the i-th snapshot to the original system. To build an aeroelastic ROM, it is possible to retain the first r-order POD modes ${{\boldsymbol \Psi} _r} = [{{\boldsymbol \psi} _1},{{\boldsymbol \psi} _2},\ldots,{{\boldsymbol \psi} _r}]$ while retaining most of the energy of the original system. By projecting the full-order series ${{\boldsymbol x}^{n \times 1}}$ on the r-order POD modes ${{\boldsymbol \Psi} _r} = [{{\boldsymbol \psi} _1},{{\boldsymbol \psi} _2},\ldots,{{\boldsymbol \psi} _r}]$ , we can reduce the full order system to a reduced r-order system

(10) \begin{align} \left\{ \begin{array}{l}{{\dot {\textbf x}}_r} = {\boldsymbol \Psi} _r^T{\textbf A_a}{{\boldsymbol \Psi} _r}{\textbf x_r} + {\boldsymbol \Psi} _r^T{\textbf B_a}\textbf y{\rm{ + }}{\boldsymbol \Psi} _r^T{\textbf L_a}{\textbf u_c}\\[5pt] {\textbf f^{ext}} = {\textbf C_a}{\textbf x_r}\end{array} \right.\end{align}

where ${\textbf A_a} = -\textbf A_0^{ - 1}\textbf H$ , ${\textbf B_a} = - \textbf A_0^{ - 1}\left[ {\textbf E + {\textbf {C G}}} \right]$ , ${\textbf L_a} = - \textbf A_0^{ - 1}\textbf L$ , ${\textbf C_a} = \textbf C{{\boldsymbol \Psi} _r}$ , $\textbf y = {\left[ {\dot {\textbf u}}\;{\textbf u} \right]^T}$ .

The above steps outline the process of generating an unsteady aerodynamics ROM in active flutter suppression. The aerodynamic force is jointly determined by the structural deformation and the deflection of control surface. The resulting aeroelastic model is obtained by coupling the structural dynamic equations, Eq. (3), with Eq. (10) for the ROM of the fluid.

2.2. LQR control design

A Linear Quadratic Regulator (LQR) regulator is used in current study for the active flutter control. The LQR assumes the availability of the entire state vector. In this work practice, some elements of the state vector may not be measurable, so sub-optimal control is generally used. The measured inputs are the displacement and velocity at wing tip where the sensors are located. First, the LQR state feedback design method states that

(11) \begin{align} {\textbf u_c} = - {\textbf K^*}{\textbf x_{ase}}\end{align}

where ${\textbf x_{ase}} = {\left[ {{\textbf x^r},{\rm{ }}\textbf y} \right]^T}$ . The suboptimal control design method assumes that only a few linear combinations of system states can be directly measured from the sensors. Let $\tilde {\textbf K} = {\textbf K^*}\textbf C_{ase}^T/({\textbf C_{ase}}\textbf C_{ase}^T)$ , ${\textbf C_{ase}} = \left[ {\textbf 0,{\rm{ }}{\textbf C_a}} \right]$ . The output vector feedback control law is formulated as

(12) \begin{align} {\textbf u_c} = - \tilde {\textbf K}{\textbf C_{ase}}{\textbf x_{ase}} = - \tilde {\textbf {Ky}}\end{align}

For a conventional LQR method, aiming to minimise the cost function

(13) \begin{align} J = \sum\limits_{k = 0}^\infty {\left\{ {{\textbf x^T}\textbf {Qx} + \textbf u_c^T\textbf R{\textbf u_c}} \right\}} \end{align}

where Q is a semi-positive definite weighting symmetric matrix, and R is a positive definite weighting symmetric matrix. The optimal problem of minimising the objective function J can be converted to the solution of:

(14) \begin{align} \begin{array}{l}\tilde {\textbf K} = {\textbf R^{ - 1}}\textbf {BP}\\[5pt] \textbf{PA} + {\textbf A^T}\textbf P - \textbf {PB}{\textbf R^{ - 1}}{\textbf B^T}\textbf P + \textbf Q = 0\end{array}\end{align}

where P is the solution of the Riccati equation, which can be calculated by many iterative algorithms. Like the conventional LQR design method, the optimal state feedback gains were obtained in this work using MATLAB /SIMULINK.

In the standard CFD-based POD/ROM construction procedure for active flutter suppression, when the structure displacement is obtained, the control law which designed by LQR is used to control the deflection of control surface at each time step, and add the influence of the control surface deflection to the aerodynamic force into the aeroelastic system to suppress flutter.

3.1. Structural dynamic reanalysis method

It is an important aspect to recalculate the eigenvalues of structure in some structural optimisation problems. For large-scale eigenvalue problem, one main obstacle is the high computational cost, and they require large computer memories and computing time. To overcome the obstacle, a number of methods including direct methods and approximate methods, have been proposed to ease the eigenvalue reanalysis [Reference Song, Chen and Sun40]. Direct methods are commonly based on the Sherman-Morrison-Woodbury (SMW) formula [Reference Akgün, Garcelon and Haftka41, Reference Huang, Wang and Li42]. But its application is limited by the rank of the modification, and also suffer from high computational costs. Different from direct methods, approximate methods aim at obtaining the modal characteristics of modified structures [Reference Beliveau, Cogan, Lallement and Ayer43–Reference Kirsch, Bogomolni and Sheinman45], in which the mode shapes are approximated as a linear combination of basic mode shapes. The advantage in doing so is a significant reduction of computational costs, and the applicability is extended for high-rank modifications of structures. The extended Kirsch combined method [Reference Chen and Yang44, Reference Chen, Yang and Lian46] is an efficient method for the case of large modifications of the structural parameters. More details relevant to this work are given next.

Consider an initial structure with stiffness matrix ${\textbf K_0}$ and mass matrix ${\textbf M_0}$ . The corresponding mode shapes ${\boldsymbol{\phi}} _0^i$ and modal frequencies $\lambda _0^i$ are calculated by solving the equations:

(15) \begin{align} {\textbf K_0}{\boldsymbol{\phi}} _0^i = \lambda _0^i{\textbf M_0}{\boldsymbol \phi} _0^i\end{align}

${\textbf K_0}$ and ${\textbf M_0}$ are perturbed into the form ${\textbf K_0} + \Delta \textbf K$ and ${\textbf M_0} + \Delta \textbf M$ , $\textbf K = {\textbf K_0} + \Delta \textbf K$ , $\textbf M = {\textbf M_0} + \Delta \textbf M$ , in which $\Delta \textbf K$ and $\Delta \textbf M$ are the perturbations in the stiffness matrix and mass matrix, respectively. The eigenvalue problem for the modified structure becomes:

(16) \begin{align} {\textbf K}{{\boldsymbol{\phi}} ^i} = {\lambda ^i}{\textbf M}{{\boldsymbol{\phi}} ^i}\end{align}

where ${\lambda ^i}$ and ${{\boldsymbol{\phi}} ^i}$ are the i-th eigenvalue and eigenvector of the modified structure, respectively.

Extended Kirsch combined method use the second-order eigenvector terms [Reference Chen, Wu and Yang47] as the basis vectors in the following modal shape reduced basis:

(17) \begin{align} {{\boldsymbol{\phi}} ^i} = {\boldsymbol{\phi}} _B^i{\textbf z^i},\quad {\textbf z^i} = {(z_0^i,{\rm{ }}z_1^i,{\rm{ }}z_2^i)^T} \in {{\textbf R}^{3 \times 1}}\end{align}

where

(18) \begin{align} {\boldsymbol{\phi}} _B^i = \left[ {{\boldsymbol{\phi}} _0^i,{\rm{ }}{\boldsymbol{\phi}} _1^i,{\rm{ }}{\boldsymbol{\phi}} _2^i} \right]\end{align}

(19) \begin{align} \lambda _1^i = {({\boldsymbol{\phi}} _0^i)^T}\Delta {{\textbf K}{\boldsymbol{\phi}}} _0^i - \lambda _0^i{({\boldsymbol{\phi}} _0^i)^T}\Delta {\textbf M}{\boldsymbol{\phi}} _0^i{\rm{ }}\end{align}

(20) \begin{align} {\boldsymbol{\phi}} _1^i = \sum\limits_{s = 1,s \ne i}^n {\frac{1}{{\lambda _0^i - \lambda _0^s}}} \left[{({\boldsymbol{\phi}} _0^s)^T}(\Delta {\textbf K} - \lambda _0^i\Delta {\textbf M}){\boldsymbol{\phi}} _0^i\right]{\boldsymbol{\phi}} _0^s - \frac{1}{2}\left[{({\boldsymbol{\phi}} _0^i)^T}\Delta {\textbf M}{\boldsymbol{\phi}} _0^i\right]{\boldsymbol{\phi}} _0^i = {{\boldsymbol{\phi}} _0}{\textbf Z}_1^i\end{align}

(21) \begin{align} \lambda _2^i & = {({\boldsymbol{\phi}} _0^i)^ T}\Delta {\textbf K}{\boldsymbol{\phi}} _1^i - \lambda _0^i{({\boldsymbol{\phi}} _0^i)^T}\Delta {\textbf M}{\boldsymbol{\phi}} _1^i - \lambda _1^i{({\boldsymbol{\phi}} _0^i)^T}{{\textbf M}_0}{\boldsymbol{\phi}} _1^i - \lambda _1^i{({\boldsymbol{\phi}} _0^i)^T}\Delta {\textbf M}{\boldsymbol{\phi}} _0^i{\rm{ }} \nonumber \\ i,s & = 1,2,3 \ldots ..,{\rm{ }}\;number{\rm{ }}\;of{\rm{ }}\;modes\end{align}

(22) \begin{align} {\boldsymbol{\phi}} _2^i & = \sum\limits_{s = 1,s \ne i}^n {\frac{1}{{\lambda _0^i - \lambda _0^s}}} \left[{({\boldsymbol{\phi}} _0^s)^T}(\Delta {\textbf K} - \lambda _0^i\Delta {\textbf M}){\boldsymbol{\phi}} _1^i - \lambda _1^i{({\boldsymbol{\phi}} _0^s)^T}({{\textbf M}_0}{\boldsymbol{\phi}} _1^i + \Delta {\textbf M}{\boldsymbol{\phi}} _0^i)\right]{\boldsymbol{\phi}} _0^s \nonumber \\ & \quad - \frac{1}{2}\left[{({\boldsymbol{\phi}} _0^i)^T}\Delta {\textbf M}{\boldsymbol{\phi}} _1^i + {({\boldsymbol{\phi}} _1^i)^T}({\textbf M_0}{\boldsymbol{\phi}} _1^i + \Delta {\textbf M}{\boldsymbol{\phi}} _0^i)\right]{\boldsymbol{\phi}} _0^i = {{\boldsymbol{\phi}} _0}{\textbf Z}_2^i\end{align}

where $\lambda _1^i$ and $\lambda _2^i$ are the i-th first-order and second-order eigenvalue of the modified structure, $\boldsymbol{\phi} _1^i$ and $\boldsymbol{\phi} _2^i$ are the i-th first-order and second-order eigenvector of the modified structure, respectively. The coefficient vector, ${\textbf z^i}$ , contains three unknowns (for a second-order perturbation). Substituting Eq. (20) and Eq. (22) into Eq. (17), ${\boldsymbol \Phi} $ can be written as

(23) \begin{align} {\boldsymbol \Phi} = \left[ {\begin{array}{*{20}{c}}{{\boldsymbol{\phi}} _0^{1}\ldots{\boldsymbol{\phi}} _0^i}\\{\boldsymbol{0}}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{{\boldsymbol{\phi}} _0^1\ldots{\boldsymbol{\phi}} _0^i}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{\boldsymbol{0}}\\{{\boldsymbol{\phi}} _0^1\ldots{\boldsymbol{\phi}} _0^i}\end{array}} \right]{\left[ {\begin{array}{{c}}{{\rm {\textbf I}}_0^1\ldots{\rm {\textbf I}}_0^i}\\{\boldsymbol{0}}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{{\textbf Z}_1^1 \ldots {\textbf Z}_1^i}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{\boldsymbol{0}}\\{{\textbf Z}_2^1 \ldots {\textbf Z}_2^i}\end{array}} \right]^T}\left[ {\begin{array}{{c}}{{{\textbf z}^1}}\\{\boldsymbol{0}}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{{{\textbf z}^2}}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{\boldsymbol{0}}\\{{{\textbf z}^3}}\end{array}} \right] = {\Phi _0}{\textbf Z}\end{align}

where

(24) \begin{align} {\textbf Z} = {\left[ {\begin{array}{*{20}{c}}{{\textbf I}_0^1 \ldots {\textbf I}_0^i}\\{\boldsymbol{0}}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{{\textbf Z}_1^1 \ldots {\textbf Z}_1^i}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{\boldsymbol{0}}\\{{\textbf Z}_2^1 \ldots {\textbf Z}_2^i}\end{array}} \right]^T}\left[ {\begin{array}{{c}}{{{\textbf z}^1}}\\{\boldsymbol{0}}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{{{\textbf z}^2}}\\{\boldsymbol{0}}\end{array}\begin{array}{{c}}{\boldsymbol{0}}\\{\boldsymbol{0}}\\{{{\textbf z}^3}}\end{array}} \right]\end{align}

Substituting Eq. (17) into the modified analysis equations Eq. (16), and premultiplying by ${({\boldsymbol{\phi}} _B^i)^T}$ , one obtains:

(25) \begin{align} {({\boldsymbol{\phi}} _B^i)^T}({\textbf K_0} + \Delta {\textbf K}){\boldsymbol{\phi}} _B^i{{\textbf z}^i} = {\lambda ^i}{({\boldsymbol{\phi}} _B^i)^T}({\textbf M_0} + \Delta \textbf M){\boldsymbol{\phi}} _B^i{z^i}\end{align}

Introducing the notation

(26) \begin{align} {\textbf K}_R^i = {({\boldsymbol{\phi}} _B^i)^T}({{\textbf K}_0} + \Delta {\textbf K}{\boldsymbol{\phi}})\end{align}

(27) \begin{align} {\textbf M}_R^i = {({\boldsymbol{\phi}} _B^i)^T}({{\textbf M}_0} + \Delta {\textbf M}){\boldsymbol{\phi}} _B^i\end{align}

and substituting Eq. (26) and (27) into Eq. (25), we can obtain a set of ( $3 \times 3$ ) matrix equation

(28) \begin{align} {\textbf K}_R^i{{\textbf z}^i} = {\lambda ^i}{\textbf M}_R^i{{\textbf z}^i}\end{align}

Thus, the coefficient vector ${{\textbf z}^i}$ is evaluated from Eq. (28). The i-th eigenvector of the modified structure is obtained by substituting ${{\textbf z}^i}$ into Eq. (17) and Z is obtained by substituting ${{\textbf z}^i}$ into Eq. (24). It needs to be noted that the first-order eigenvector of Eq. (28) is only used [Reference Chen and Yang44].

Finally, the i-th eigenvalue of the modified structure, $\lambda _K^i$ , is computed using the Rayleigh quotient:

(29) \begin{align} \lambda _K^i = \frac{{{{({{\boldsymbol{\phi}} ^i})}^T}({{\textbf K}_0} + \Delta {\textbf K}){{\boldsymbol{\phi}} ^i}}}{{{{({{\boldsymbol{\phi}} ^i})}^T}{{\textbf M}_0}{{\boldsymbol{\phi}} ^i}}}\end{align}

3.2. Efficient reduced-order framework

In section 2.1, the unsteady aerodynamics ROM should be regenerated for each structural modification. Therefore, it cannot be used for active/passive hybrid flutter suppression. In the proposed efficient reduced-order framework, the initial structural mode shapes ${{\boldsymbol \Phi} _0}$ is taken as the basic mode shapes for basic POD/ROM construction. For the modified structure, the physical displacement ${\boldsymbol \Theta} $ of the wing can be written as

(30) \begin{align} {\boldsymbol \Theta} = {\boldsymbol \Phi} {\textbf u}\end{align}

Substituting Eq. (23) into Eq. (30), the physical displacement of the wing can also be written as

(31) \begin{align} {\boldsymbol \Theta} = {{\boldsymbol \Phi} _0}(\textbf {Zu})\end{align}

and ${{\textbf u}_b} = {\textbf {Zu}}$ , ${{\textbf u}_b}$ is artificially defined as basic generalised displacements.

For the modified structure in the efficient reduced-order framework, the matrices G, E and C of the linearised flow solver, Eq. (1), substituting the relation ${{\textbf u}_b} = {\textbf {Zu}}$ , the matrices may be rewritten in terms of the vector of basic generalised displacements:

\begin{align*}{\textbf G} = \frac{{\partial {\textbf F}}}{{\partial {\textbf u}}}({{\textbf w}_0},{{\textbf u}_0},{\dot {\textbf u}_0}) = \frac{{\partial {\textbf F}}}{{{{\textbf Z}^{ - 1}}\partial {{\textbf u}_b}}}({{\textbf w}_0},{{\textbf u}_0},{\dot {\textbf u}_0}) = {\textbf Z}{{\textbf G}_b}\end{align*}

(32) \begin{align} {\textbf E} = \frac{{\partial {\textbf A}}}{{\partial {\textbf u}}}{{\textbf w}_0} = \frac{{\partial {\textbf A}}}{{{{\textbf Z}^{ - 1}}\partial {{\textbf u}_b}}}{{\textbf w}_0} = {\textbf Z}{{\textbf E}_b}\end{align}

\begin{align*}{\textbf C} = \frac{{\partial {\textbf F}}}{{\partial \dot {\textbf u}}}({{\textbf w}_0},{{\textbf u}_0},{\dot {\textbf u}_0}) = \frac{{\partial {\textbf F}}}{{{{\textbf Z}^{ - 1}}\partial {{\dot {\textbf u}}_b}}}({{\textbf w}_0},{{\textbf u}_0},{\dot {\textbf u}_0}) = {\textbf Z}{{\textbf C}_b}\end{align*}

where $\textbf{G}_{b}$ , $\textbf{E}_{b}$ , $\textbf{C}_{b}$ are the first order terms in a Taylor series expansion of the basis reduced r-order aeroelastic model. In flutter suppression, since the elastic influence of the control surface can be ignored, the gradient of the numerical flux function with respect to control surface deflection $\textbf{L}_{a}$ is not changed. Now, the ROM in the proposed efficient reduced-order framework for active/passive hybrid flutter suppression is written as

(33) \begin{align} \left\{ \begin{array}{l}{{\dot {\textbf x}}_r} = {\boldsymbol \Psi} _r^T{\textbf A}{{\boldsymbol \Psi} _r}{{\textbf x}_r} + {\boldsymbol \Psi} _r^T{\textbf Z}{{\textbf B}_b}{\textbf y} + {\boldsymbol \Psi} _r^T{{\textbf L}_a}{{\textbf u}_c}\\[5pt] {{\textbf f}^{ext}} = {\textbf Z}{{\textbf C}_b}{{\boldsymbol \Psi} _r}{{\textbf x}_r}\end{array} \right.\end{align}

where ${\textbf A} = - {\textbf A}_0^{ - 1}{\textbf H}$ , ${{\textbf B}_b} = - {\textbf A}_0^{ - 1}\left[ {{{\textbf E}_b} + {{\textbf C}_b}\; {{\textbf G}_b}} \right]$ , ${\textbf y} = {\left[ {\dot {\textbf u}}\;{\textbf u} \right]^T}$ , ${{\textbf C}_b} = {{\textbf P}_b}$ .

The structural dynamic equations of modified structure in the efficient framework are:

(34) \begin{align} \bar {\textbf M}\ddot {\textbf u} + \bar {\textbf {Ku}} = {{\textbf f}^{ext}}\end{align}

here, $\bar {\textbf M} = {{\textbf Z}^T}{\boldsymbol \Phi} _0^T({{\textbf M}_0} + \Delta {\textbf M}){\Phi _0}{\textbf Z}$ , $\bar {\textbf K} = {{\textbf Z}^T}{\boldsymbol \Phi} _0^T({{\textbf K}_0} + \Delta {\textbf K}){{\boldsymbol \Phi} _0}{\textbf Z}$ . In addition, when $\Delta {\textbf K} = {\boldsymbol{0}}$ and $\Delta {\textbf M} = {\boldsymbol{0}}$ , Z is the unit matrix, Eq. (33) and Eq. (34) are equivalent to Eq. (10) and Eq. (3), respectively. For active/passive hybrid flutter suppression, adding a mass balance and changing the structural stiffness were made by varying the structural parameters, both the structural model and the time-consuming CFD-based ROM have to be reconstructed to ensure accuracy. When the aeroelastic response is obtained, the control surface rotation is updated at each time step by control law which is designed by LQR and then add the influence of deflection of control surface to the aerodynamic force into the aeroelastic system to achieve the purpose of flutter suppression. The aeroelastic response of the modified structure is then obtained by the proposed efficient reduced-order framework, without the expensive time-consuming reconstruction procedure of the new POD basis. Therefore, the proposed framework can significantly reduce the computational cost. The proposed framework provides a potential powerful tool for active/passive hybrid flutter suppression that would be convenient to design more effective controllers to improve aircraft fatigue life and enhance aircraft performance.

4.1. POD/ROM solver validation

The accuracy of POD/ROM solver is validated by the AGARD 445.6 aeroelastic wing model [Reference Yates48, Reference Silva, Chwalowski and Perry49]. The aerofoil section is a NACA 65A004. The density of wing material is 381.98 kg/m³. The elastic modulus in the spanwise direction and the chordwise direction is 3.151Gpa and 0.416GPa, respectively. The shear modulus is 0.4392GPa, and the Poisson’s ratio is 0.31 [Reference Zhong and Xu50]. The structural model based on shell elements, shown in Fig. 1(a), consists of 231 nodes and 200 elements, and the thickness distribution of the wing was governed by aereofoil shape. For the fluid model, a multi-block structured mesh was employed, as shown in Fig. 1(b). The spatial convergence of the CFD mesh was analysed in Zhou et al. [Reference Qiang, Chen, Li and da Ronch51], which reported a good agreement of the results from both the medium and fine grids. In this paper, the total number of grid points on the medium grid, herein used, is 223,146 ( $99 \times 49 \times 46$ ).

Figure 1. AGARD 445.6 wing: (a) structural model and trailing-edge control surface, and (b) surface CFD mesh.

Figure 2 shows the flutter boundary computed by the full-order CFD/CSD aeroelastic model and POD/ROM solver, experimental data [Reference Yates48] are also included. The agreement between the POD/ROM, CFD/CSD aeroelastic model and experimental data is good for all Mach numbers considered (0.499 to 1.141), including the well-known transonic dip of the flutter speed. The accuracy of full-order CFD/CSD aeroelastic model and POD/ROM solver have been evaluated over the years in a number of aeroelastic studies [Reference Zhou, Chen and Li16, Reference Chen, Li, Zhou, da Ronch and Li34, Reference Li, da Ronch, Chen and Li36, Reference Zhou, Li, da Ronch, Chen and Li52, Reference Chen, Sun and Li53].The good agreement indicates that the POD/ROM solver has near the same accuracy as the nonlinear CFD model at small or moderate perturbation.

Figure 2. AGARD 445.6 wing flutter boundary.

4.2. Accuracy evaluation of the efficient aeroelastic CFD-based POD/ROM

The accuracy and efficiency of the efficient CFD-based POD/ROM has been successfully verified in [Reference Chen, Li, Zhou, da Ronch and Li34]. The AGARD 445.6 wing structural model is divided into four sections along spanwise direction, as shown in Fig. 1(a). Each section of the identical material properties, but material properties at different sections may be different. The material properties are assumed to vary as

1. Section 1, ${E_1} = (1 + 3\varepsilon ){E_0}$ , ${\rho _1} = (1 + 3\varepsilon ){\rho _0}$

2. Section 2, ${E_2} = (1 + 2\varepsilon ){E_0}$ , ${\rho _2} = (1 + 2\varepsilon ){\rho _0}$

3. Section 3, ${E_3} = (1 + \varepsilon ){E_0}$ , ${\rho _3} = (1 + \varepsilon ){\rho _0}$

4. Section 4, ${E_4} = {E_0}$ , ${\rho _4} = {\rho _0}$

where ${E_0}$ and ${\rho _0}$ are the Young’s modulus and density of the original wing structure, respectively. It is worth observing that the choice above is used for demonstration purposes and shall not be taken here as a limiting case of the present method.

Three cases for the improved AGARD 445.6 wing are considered: $\varepsilon$ = 1/12, 1/6, 1/3. For conciseness, the flutter speed at Ma = 0.960 (air density is 0.06341 kg/m³) and angle-of-attack AOA = 0deg (the flutter speed is 294m/s for standard AGARD 445.6 wing model) using the efficient aeroelastic CFD-based POD/ROM is presented in Table 1. Reference data is taken from the exact model in which both the structural model and the CFD-based POD/ROM have to be reconstructed for each modified structure. Although the flutter speed error increase for increasing level of the structural modifications, it is worth observing that the maximum difference is limited to about 2%. It indicates that the efficient CFD-based POD/ROM has good accuracy for flutter speeds prediction, even for very large range ( $\varepsilon$ = 1/3) of structural parameter variation.

Table 1. Flutter speed obtained by exact model and the efficient aeroelastic CFD-based POD/ROM

4.3. Accuracy evaluation of the efficient reduced-order framework

After evaluating the accuracy of POD/ROM solver and the efficient aeroelastic CFD-based POD/ROM, the whole proposed efficient reduced-order framework for active/passive hybrid flutter suppression will be evaluated using an improved AGARD 445.6 wing model in this section.

The structural model of the improved AGARD 445.6 wing were modified to accommodate a trailing-edge control surface at the tip of the wing, with dimensions equal to 20% of the wing span and 30% of the wing chord, is shown in Fig. 1(a). Two sensors were added at the wing tip to measure the displacements and velocities of the wing. The motion of the control surface is treated as an additional mode shape. The mode shape for the control surface deflection on CFD surface, is shown in Fig. 3. It can be seen that the control surface deflection on the CFD surface keeps the form of control surface deflection of structural model. It needs to be noted that rotations of the control surface are small, so modeling gaps can be ignored [Reference Zhou, Li, da Ronch, Chen and Li52].

Figure 3. The mode shape for the control surface deflection on the CFD surface (in red).

For simplifying the description, only the aeroelastic responses at Ma = 0.960 (air density is 0.06341Kg/m³) and AOA = 0deg are presented in this paper. For the three cases ( $\varepsilon$ = 1/12, 1/6, 1/3), the freestream speed is set to 340m/s, 360m/s and 395m/s, respectively, which are about 9% higher than the flutter speed. The weight matrices, Q ₁ = 0.1I ₆₀ × ₆₀ and R = 10, were used for the LQR synthesis. Figures 4–6 show the aeroelastic responses obtain by proposed efficient reduced-order framework. Reference data are representative of the exact system, where the aeroelastic CFD-based POD/ROM has to be reconstructed for each modified structure. For all cases, although the discrepancy between the two methods increases for increasing level of the structural modifications, the good agreement between the exact system and the proposed framework are founded. The good agreement indicate that the proposed efficient reduced-order framework has good predictive capability for aeroelastic response. As can be seen from Figs. 4–6, at same freestream speed condition, the aeroelastic response of passive control alone continuously diverges that the system is unstable, and the aeroelastic responses of the hybrid control is convergent that the system is stable. It indicates that the proposed framework is ideally suitable for active/passive hybrid flutter suppression to design more effective controllers with good accuracy. In addition, Figs. 4-6 also indicate that the efficient aeroelastic CFD-based POD/ROM has good accuracy for aeroelastic responses prediction.

Figure 4. Comparison of the aeroelastic response for $\varepsilon$ = 1/12 (V _∞=340m/s).

Figure 5. Comparison of the aeroelastic response for $\varepsilon$ = 1/6 (V _∞=365m/s).

Figure 6. Comparison of the aeroelastic response for $\varepsilon$ = 1/3 (V _∞=395m/s).

To further demonstrate the predictive capability of the proposed efficient reduced-order framework, the fifth modal coordinate related to the rotation of the trailing-edge control surface is presented in Fig. 7. The reference data are taken from the exact system. It can be seen that in all cases, the rotation of the trailing-edge control surface is consistent. It indicates again that the proposed framework has good accuracy for predict aeroelastic response, without reconstructing a new set of CFD-based POD/ROM for modified structure. All the above comparisons indicate that the proposed efficient reduced-order framework can capture the aeroelastic response with good accuracy, and ideally suitable for active/passive hybrid flutter suppression to design more effective controllers with high accuracy. It should be noted that the proposed framework would work well for small to moderate deformations but not for large variations.

Figure 7. Trailing-edge control surface rotation.

Finally, a quantification of flutter speed errors for the proposed efficient reduced-order framework and exact system is presented in Table 2. It is not unexpected that the accuracy of the proposed framework degrades with the increasing structural modifications. However, it is worth observing that the maximum difference is limited to about 2.5%. For all cases, compared with passive control alone, the flutter speed increased by 36.3%, 36.0% and 37.4%, respectively. Compared with active control alone [Reference Zhou, Li, da Ronch, Chen and Li52], the flutter speed increased by 7.3%, 13.9% and 27.3%, respectively. It is indicate that compared with passive flutter suppression alone [Reference Chen, Li, Zhou, da Ronch and Li34] and active flutter suppression alone [Reference Zhou, Li, da Ronch, Chen and Li52], the proposed framework can greatly expand the flutter boundary. All the above comparisons indicate that the proposed efficient reduced-order framework is very suitable for active/passive hybrid flutter suppression to design more effective controllers with good accuracy.

Table 2. Flutter speed obtained by exact system and the efficient reduced-order framework

4.4. Efficiency evaluation of the efficient reduced-order framework

The computational efficiency is one of the most important criteria of the proposed efficient reduced-order framework. All analyses were performed on a Windows 10 system PC with Intel^® Core^(TM) i7-9700K CPU (3.60 GHz, 8 cores, but only one core used) and 32 GB RAM.

For the exact system, each modified structure need regenerate a set of POD modes, requiring about 16h per configuration. For the three cases, required about 48h. In contract, the set of POD modes is generated only once for the proposed efficient reduced-order framework. For the proposed efficient reduced-order framework, it only requires about 16h. It is reasonable to consider about 100 cases and 20 values of the freestream dynamic pressure to assess the aeroelastic stability. With information reported in Table 3, the exact system would require over 1,600 CPU hours (this consists of 16h × 100 and 2.69s × 20 × 100, totaling 1,609.49h), whereas the proposed framework just require over 17 CPU hours (16h × 1 and 2.92s × 20 × 100, totaling 17.62h). It is obvious that the expected speed-up is proportional to the number of modified structural models, the computational advantage of the efficient reduced order framework becomes more significant the more the number of modified structural model is considered. All the demonstration cases indicated that the proposed efficient reduced-order framework has high computational efficiency and can be applied for active/passive hybrid flutter suppression to design more effective controllers.

Table 3. Computational cost of the exact system and the efficient reduced order framework