Effect of Aerodynamic Damping Approximations on Aeroelastic Eigensensitivities

Abstract

Aeroelastic sensitivities for the flutter solution are a crucial component of the multi-disciplinary optimization methods employed in modern aircraft design. This paper derives the aeroelastic sensitivities for different aerodynamic damping approximations—the p-k method, the g method and the generalized aeroelastic analysis method—highlighting the influence of the employed aerodynamic approximation on the required derivatives. The derived formulation for the determination of aeroelastic sensitivities by means of a direct method is verified for the case of a two-degree-of-freedom typical section model, where analytical aeroelastic sensitivities can be analytically obtained. For this and for an additional model, namely the AGARD 445.6 weakened wing, the significant effect of the aerodynamic damping approximation on the aeroelastic sensitivities is shown.

1. Introduction

The design of new aircraft entails a demand for high efficiency, which eventually is achieved by extensive multi-disciplinary design optimization (MDO) while considering constraints on stability and safety. Therefore, the MDO must ensure the aircraft to be free of flutter conditions [1] within the flight envelope and the aeroelastic stability margin as required by certification. Aeroelastic stability is typically expressed by the eigenvalues of the aeroelastic eigenproblem formulated in the Laplace domain. Hence, the eigensensitivities of the aeroelastic eigenproblem facilitate the design process for compliance with flutter constraints as they describe the change in stability by changing design parameters. This applies especially for gradient-based MDO, where sensitivities are required to fulfill imposed constraints and obtain the next iterated design [2,3].

This paper follows up on the work of Cardani and Mantegazza [4], as well as Bindolino and Mantegazza [5], who have thoroughly described the computation of eigenvalue and eigenvector derivatives for flutter eigenproblems. Murthy and Haftka [6] have given a detailed description on solving for the eigensensitivities using the direct and adjoint methods. Additionally, the commercial solver MSC Nastran provides the computation of eigenvalue sensitivities [7], based on [8], employing the doublet lattice method (DLM) [9].

Extending the achievements made so far, this work addresses the influence of aerodynamic damping approximations on the aeroelastic eigensensitivities. Unsteady aerodynamic methods providing the air loads for the aeroelastic eigenproblem by the aerodynamic transfer function matrix typically account for purely harmonic structural oscillations. Therefore, the present damping of the oscillations as described by the real part of the eigenvalues is not captured in the aerodynamic loads. This shortcoming affects the solution of the aeroelastic eigenproblem, leading to several aeroelastic solution methods. The well-known p-k method provides the flutter solution for purely harmonic air loads by iteratively matching the frequencies of the oscillations and the air loads [10]. The g method developed by Chen [11] extends the p-k method for small values of the real part of the eigenvalues by providing a first-order approximation of the analytic continuation of the aerodynamic transfer function matrix. Edwards et al. [12] showed that, for linear potential flow, the complete analytic continuation is valid for the aeroelastic analysis and termed it the generalized aeroelastic analysis method (GAAM).

Existing work regarding the computation of aeroelastic eigensensitivities in the context of gradient-based MDO has typically focused on the precision of the derivatives involved, obtained either analytically by automatic differentiation or numerically by finite differences or complex step formulas [13], but using the common p-k approximation when computing the aerodynamic damping term. This work focuses on the effect that different levels of physical approximation used for the determination of the aerodynamic damping term may have on the aeroelastic eigensensitivites, as this influence may overtake the approximation error introduced by the numerical method used for the computation of the derivatives. Thus, the present work focuses on differentiating the nonlinear aeroelastic eigenproblem, identifying the involved sensitivities of the structural dynamics and different aerodynamic damping approximations. To this aim, the aerodynamic sensitivities for the p-k method, the g method and GAAM are derived, which allows us to accurately obtain the aeroelastic eigensensitivities for each method and display the difference in the resulting eigenvalue sensitivities. This is demonstrated for a typical section model as well as for the AGARD 445.6 wing. Furthermore, the derived aeroelastic sensitivities are verified for the typical section model employing Theodorsen’s aerodynamic airfoil theory by obtaining the analytical derivatives.

In Section 2.1, the nonlinear aeroelastic eigenproblem is introduced, as well as the employed solution method in Section 2.2, and the investigated aerodynamic damping approximations are presented in Section 2.3. Section 3.1 and Section 3.2 describe the aeroelastic sensitivities and their solution by the direct method in physical and modal coordinates employing GAAM. The aerodynamic sensitivities of the p-k and g methods are derived in Section 3.3 and Section 3.4. Section 3.5 presents the sensitivity of the modal matrix for employing modal coordinates. The verification for the typical section model for the three approximations is shown in Section 4.1 and Section 4.2, showing the eigenvalue sensitivities for the AGARD 445.6 wing with varying sweep obtained with the p-k and g method. For the verification, the derivatives of the generalized Theodorsen aerodynamic airfoil theory are presented together with the derivatives of the generalized Theodorsen function in Appendix A and Appendix B.

2. Aeroelastic Stability

2.1. Aeroelastic Eigenvalue Problem

The aeroelastic governing equations couple the structural dynamics with the aerodynamic forces by second-order differential equations resulting from Newton’s second law. For obtaining linear stability, complex harmonic motion $\widehat{\mathit{x}}{\mathrm{e}}^{st}$ around a steady state is assumed. This transforms the second-order differential equations into the Laplace domain with the Laplace variable $s=\sigma +\mathrm{i}\omega$ as a system of nonlinear, algebraic equations

$$\begin{array}{cc}\hfill (s^2\mathit{M}+s\mathit{D}+\mathit{K}-\mathit{A}\left(s\right))\widehat{\mathit{x}}=\mathit{G}\left(s\right)\widehat{\mathit{x}}& \hfill =\mathbf{0},\hfill \end{array}$$

(1)

where the unsteady aerodynamic forces are described in a time-linearized manner by the aerodynamic transfer function matrix $\mathit{A}\left(s\right)$ multiplied by the amplitude of motion $\widehat{\mathit{x}}$. The complex-valued matrix $\mathit{A}\left(s\right)$ especially depends on steady-state flow parameters such as Mach number and dynamic pressure or Reynolds number according to the employed aerodynamic method. The mass matrix $\mathit{M}$, the damping matrix $\mathit{D}$ and the stiffness matrix $\mathit{K}$ are real symmetric matrices describing the structural properties. Equation (1) is a complex-valued, nonlinear eigenproblem with eigenvalues s and nonzero eigenvectors $\widehat{\mathit{x}}$. Consequently, the real part of the eigenvalues indicates growing, $\sigma >0$, or decaying, $\sigma <0$, oscillations. Hence, the instability onset is given for $\sigma =0$ denoted as the flutter onset.

The aeroelastic governing equations and, thus, Equation (1) are typically applied in modal coordinates in order to reduce the system size by modal truncation, where only a small number of structural mode shapes are considered. Therefore, the amplitudes of the physical displacements $\widehat{\mathit{x}}$ are transformed into the modal participation factors $\widehat{\mathit{q}}$ by the modal matrix

$$\begin{array}{cc}\hfill \mathsf{\Phi }& \hfill =\left[\begin{array}{ccc}\cdots & {\mathit{\phi }}_j& \cdots \end{array}\right],\hfill \end{array}$$

which consists of structural mode shapes $\mathit{\phi }$ obtained from linear modal analysis of the structural model; see Section 3.5. The transformed Equation (1) in modal coordinates results in

$$\begin{array}{cc}\hfill (s^2\tilde{\mathit{M}}+s\tilde{\mathit{D}}+\tilde{\mathit{K}}-\tilde{\mathit{A}}\left(s\right))\widehat{\mathit{q}}={\mathsf{\Phi }}^T\mathit{G}\left(s\right)\mathsf{\Phi }\widehat{\mathit{q}}=\tilde{\mathit{G}}\left(s\right)\widehat{\mathit{q}}& \hfill =\mathbf{0}\hfill \end{array}$$

(2)

with the modal mass matrix $\tilde{\mathit{M}}$, modal damping matrix $\tilde{\mathit{D}}$, modal stiffness matrix $\tilde{\mathit{K}}$ and the generalized aerodynamic force matrix $\tilde{\mathit{A}}\left(s\right)$. It is noted that, without modal truncation, the solutions of Equation (2) are identical to those of Equation (1), with the eigenvectors related by $\widehat{\mathit{x}}=\mathsf{\Phi }\widehat{\mathit{q}}$. Furthermore, the modal approach eases the use of different spatial discretizations for the structural and aerodynamic disciplines by employing the common set of modal coordinates.

2.2. Solutions of the Aeroelastic Eigenproblem

A general method for obtaining the solutions of the nonlinear aeroelastic eigenproblem is by finding the roots of Equations (1) and (2) when considered as nonlinear functions of the vector of eigenpairs $\mathbf{s}={\left[s{\widehat{\mathit{x}}}^T\right]}^T$ and $\tilde{\mathbf{s}}={\left[s{\widehat{\mathit{q}}}^T\right]}^T$, respectively. In order to equal the number of unknowns, this approach requires one additional equation that normalizes the eigenvectors. The employed equations for the eigenvector normalization are ${\widehat{\mathit{x}}}^T\mathit{W}\widehat{\mathit{x}}-1=0$ and ${\widehat{\mathit{q}}}^T\tilde{\mathit{W}}\widehat{\mathit{q}}-1=0$, where the weighting matrices $\mathit{W}$ and $\tilde{\mathit{W}}$ are free to choose symmetric matrices including the identity matrix; see also Section 3.2. The evaluation of the aerodynamic transfer function matrix $\mathit{A}\left(s\right)$ is directly integrated in the function evaluation, ensuring the matched values of the aerodynamic transfer function matrix with the iterated eigenvalues s. In this work, the root finding is performed by a quasi-Newton method with a Krylov approximation for the inverse of the Jacobian [14]. The quasi-Newton approach requires good initial estimates in order to properly converge to the roots. Therefore, the structural eigenvectors and eigenfrequencies are used as starting solutions for the wind-off condition $\mathit{A}\left(s\right)=0$. Then, by stepwise increasing $\mathit{A}\left(s\right)$, e.g., by a velocity or density sweep, and solving the nonlinear eigenproblem at each step, the starting solutions are marched up to the final solution. Other starting points may be required to find static divergence [12] or the roots of fluid modes [15].

2.3. Aerodynamic Damping Approximations

The aerodynamic transfer function matrix $\mathit{A}\left(s\right)$ of Equations (1) and (2) serves as a general description of the linearized unsteady aerodynamic forces. Depending on the employed aerodynamic method, no explicit expression can be determined and it has to be determined numerically. Typically, aerodynamic methods make use of the nondimensional frequency, which is the reduced complex-valued frequency

$$\begin{array}{c}\hfill s^{*}=\frac{L}{V}s={\sigma }^{*}+\mathrm{i}{\omega }^{*},\end{array}$$

(3)

where L is the characteristic length, V the freestream velocity and ${\omega }^{*}$ the reduced frequency. In the general case, termed the generalized aeroelastic analysis method (GAAM) by Edwards and Wieseman [12], $\mathit{A}\left(s^{*}\right)$ is provided for the both the real and imaginary part of $s^{*}$ and, therefore, the aerodynamic damping is exactly incorporated. Recently, a modification to the extensively used doublet lattice method (DLM) [9] allowed for the application of GAAM to general nonplanar configurations in subsonic compressible flow [16]. However, if the evaluation of the aerodynamic transfer function matrix is limited to harmonic motion that is the evaluation at only the imaginary part of $s^{*}$

$$\begin{array}{c}\hfill {\mathit{A}}_k\left({\omega }^{*}\right)=\mathit{A}\left(\mathrm{i}{\omega }^{*}\right),\end{array}$$

(4)

this leads to the well-known p-k method [10]. In this case, the evaluation of the aerodynamic transfer function matrix disregards the damping coefficient ${\sigma }^{*}$ and, thus, it is a constant approximation for ${\sigma }^{*}\ne 0$. Another approximation is the g method by Chen [11], which extends the p-k method by a first-order approximation for ${\sigma }^{*}$ at ${\sigma }^{*}=0$, resulting in

$$\begin{array}{c}\hfill {\mathit{A}}_g\left(s^{*}\right)=\mathit{A}\left(\mathrm{i}{\omega }^{*}\right)+\frac{\partial \mathit{A}\left(\mathrm{i}{\omega }^{*}\right)}{\partial {\sigma }^{*}}{\sigma }^{*}=\mathit{A}\left(\mathrm{i}{\omega }^{*}\right)-\mathrm{i}\frac{\partial \mathit{A}\left(\mathrm{i}{\omega }^{*}\right)}{\partial {\omega }^{*}}{\sigma }^{*}={\mathit{A}}_k\left({\omega }^{*}\right)-\mathrm{i}\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\omega }^{*}}{\sigma }^{*}.\end{array}$$

(5)

Equation (5) shows that the g method uses the first-order derivative of ${\mathit{A}}_k\left({\omega }^{*}\right)$, which, in general, is available either analytically, by finite difference approximations or interpolation. This first-order approximation for ${\sigma }^{*}$ is found from the assumption that $\mathit{A}\left(s^{*}\right)$ is analytic and, thus, complex-differentiable at $s^{*}$ ([17] 1.9(ii)). Hence, the Cauchy–Riemann equations apply, relating the partial derivatives of $\mathit{A}\left(s^{*}\right)$ with respect to the real and imaginary part of $s^{*}$ by

$$\begin{array}{c}\hfill \frac{\partial \mathit{A}\left(s^{*}\right)}{\partial {\sigma }^{*}}=-\mathrm{i}\frac{\partial \mathit{A}\left(s^{*}\right)}{\partial {\omega }^{*}}.\end{array}$$

Employing this property for ${\sigma }^{*}=0$ results in the first-order term for ${\sigma }^{*}$ in Equation (5) and, subsequently, ${\mathit{A}}_g\left(s^{*}\right)$ is analytic at ${\sigma }^{*}=0$.

Recently, the p-L method, which recovers a true aerodynamic damping representation from the values at the imaginary axis, has been presented [18]. Even though the aeroelastic eigensensitivities produced by the p-L method could be embedded in the framework presented in this work, they can also be obtained by using its corresponding generalized state-space form. This formulation shall be presented in a separate work.

Overall, the p-k method, the g method and GAAM yield the same aerodynamic forces for ${\sigma }^{*}=0$, which implies that the three methods return the same flutter onset but their solutions are different for ${\sigma }^{*}\ne 0$. Furthermore, the methods differ in their analytic properties: the p-k method is not analytic, the g method is analytic at reduced frequencies with ${\sigma }^{*}=0$, and GAAM is considered to be analytic at typical values of $s^{*}$.

3. Aeroelastic Sensitivities

3.1. Derivatives of the Aeroelastic Eigenproblem

The aeroelastic sensitivities for the eigenvalues of the aeroelastic eigenproblem in physical coordinates, Equation (1), are obtained by differentiating Equation (1) with respect to a design parameter $\beta$. In terms of the dynamic aeroelastic matrix $\mathit{G}\left(s\right)$, the differentiated equation results in

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathit{G}\left(s\right)}{\mathrm{d}\beta }\widehat{\mathit{x}}+\mathit{G}\left(s\right)\frac{\mathrm{d}\widehat{\mathit{x}}}{\mathrm{d}\beta }=\frac{\partial \mathit{G}\left(s\right)}{\partial \beta }\widehat{\mathit{x}}+\frac{\partial \mathit{G}\left(s\right)}{\partial s}\widehat{\mathit{x}}\frac{\mathrm{d}s}{\mathrm{d}\beta }+\mathit{G}\left(s\right)\frac{\mathrm{d}\widehat{\mathit{x}}}{\mathrm{d}\beta }=0,\end{array}$$

(6)

where the sought eigenvalue derivative $\mathrm{d}s/\mathrm{d}\beta$ appears by the chain rule with the partial derivative $\partial \mathit{G}\left(s\right)/\partial s$ evaluated at s. Moreover, Equation (6) contains the eigenvector derivative $\mathrm{d}\widehat{\mathit{x}}/\mathrm{d}\beta$ as well as the partial derivative $\partial \mathit{G}\left(s\right)/\partial \beta$ at s. Solving for $\mathrm{d}s/\mathrm{d}\beta$ in Equation (6) must factor in the also unknown eigenvector derivative $\mathrm{d}\widehat{\mathit{x}}/\mathrm{d}\beta$, leading to the direct method described in Section 3.2.

According to Equation (1), the partial derivative $\partial \mathit{G}\left(s\right)/\partial \beta$ aggregates the sensitivities of the aeroelastic model matrices to

$$\frac{\partial \mathit{G}\left(s\right)}{\partial \beta }=s^2\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\beta }+s\frac{\mathrm{d}\mathit{D}}{\mathrm{d}\beta }+\frac{\mathrm{d}\mathit{K}}{\mathrm{d}\beta }-\frac{\partial \mathit{A}\left(s\right)}{\partial \beta }$$

(7)

with the partial derivative of the aerodynamic transfer function matrix $\partial \mathit{A}\left(s\right)/\partial \beta$ evaluated at s. By categorizing the design parameters in structural and geometrical parameters, the partial derivative $\partial \mathit{A}\left(s\right)/\partial \beta$ only applies for geometrical design parameters affecting the aerodynamic shape. The partial derivative $\partial \mathit{G}\left(s\right)/\partial s$ comprises the derivative of the polynomial part for the structural dynamics and the partial derivative of the aerodynamic transfer function matrix, resulting in

$$\frac{\partial \mathit{G}\left(s\right)}{\partial s}=2s\mathit{M}+\mathit{D}-\frac{\partial \mathit{A}\left(s\right)}{\partial s}.$$

(8)

In the general case, GAAM, where the aerodynamic transfer function matrix is given in terms of the reduced complex-valued frequency $s^{*}$, the partial derivatives must account for the chained function of $s^{*}$; see Equation (3). Hence, the partial derivatives of the aerodynamic transfer function matrix in Equations (7) and (8) expand to

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}\left(s\right)}{\partial s}& \hfill =\frac{\partial \mathit{A}\left(s^{*}\right)}{\partial s^{*}}\frac{\partial s^{*}}{\partial s}=\frac{L}{V}\frac{\partial \mathit{A}\left(s^{*}\right)}{\partial s^{*}},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}\left(s\right)}{\partial \beta }& \hfill =\frac{\partial \mathit{A}\left(s^{*}\right)}{\partial \beta }+\frac{\partial \mathit{A}\left(s^{*}\right)}{\partial s^{*}}\frac{\partial s^{*}\left(s\right)}{\partial \beta }=\frac{\partial \mathit{A}\left(s^{*}\right)}{\partial \beta }+\frac{\partial \mathit{A}\left(s^{*}\right)}{\partial s^{*}}\frac{\partial (sL/V)}{\partial \beta },\hfill \end{array}$$

(10)

where the latter term $\partial (sL/V)/\partial \beta$ is zero unless $\beta$ refers to one of the characteristic values L, V. Importantly, it is noted that the partial derivatives with respect to s or $s^{*}$ demand the aerodynamic transfer function matrix to be analytic, i.e., complex-differentiable, at the evaluated complex-valued frequencies. Consequently, the derived equations must be reformulated for the application with the p-k and g methods; see Section 3.3 and Section 3.4.

Employing modal coordinates for the aeroelastic eigenproblem instead of physical coordinates leads to differentiating Equation (2), which can be expanded to

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\mathsf{\Phi }}^T}{\mathrm{d}\beta }\mathit{G}\left(s\right)\mathsf{\Phi }\widehat{\mathit{q}}+{\mathsf{\Phi }}^T\frac{\mathrm{d}\mathit{G}\left(s\right)}{\mathrm{d}\beta }\mathsf{\Phi }\widehat{\mathit{q}}+{\mathsf{\Phi }}^T\mathit{G}\left(s\right)\frac{\mathrm{d}\mathsf{\Phi }}{\mathrm{d}\beta }\widehat{\mathit{q}}+{\mathsf{\Phi }}^T\mathit{G}\left(s\right)\mathsf{\Phi }\frac{\mathrm{d}\widehat{\mathit{q}}}{\mathrm{d}\beta }& \hfill =\mathbf{0}\hfill \end{array}$$

(11)

by the product rule, revealing the derivative of the modal matrix $\mathrm{d}\mathsf{\Phi }/\mathrm{d}\beta$ and the modal eigenvector derivative $\mathrm{d}\widehat{\mathit{q}}/\mathrm{d}\beta$. The derivative $\mathrm{d}\mathit{G}\left(s\right)/\mathrm{d}\beta$ entails the eigenvalue derivative $\mathrm{d}s/\mathrm{d}\beta$ as well as the partial derivatives $\partial \mathit{G}\left(s\right)/\partial \beta$ and $\partial \mathit{G}\left(s\right)/\partial s$; cf. Equation (6). The modal matrix derivative comprising structural eigenvector derivatives is derived from the structural model as described in Section 3.5.

3.2. Direct Method with Eigenvector Normalization

The direct method solves Equations (6) and (11) for the eigenvalue and eigenvector derivatives by appending the eigenvector normalization in order to complete the set of equations [4,5,6]. In this work, the eigenvectors are normalized with a symmetric weighting matrix $\mathit{W}$ by

$${\widehat{\mathit{x}}}^T\mathit{W}\widehat{\mathit{x}}=1,$$

(12)

which is in line with the solution approach for the eigenproblem; see Section 2.2. It is acknowledged that this equation is a degenerate quadratic form for the complex-valued eigenvectors, which can yield a zero norm for nonzero vectors [6]. However, this condition did not occur in any of the investigations carried out.

Differentiating Equation (12) with respect to $\beta$ states the additional condition that the eigenvector normalization constraint remains fulfilled, resulting in the additional equation

$$\frac{\mathrm{d}{\widehat{\mathit{x}}}^T}{\mathrm{d}\beta }\mathit{W}\widehat{\mathit{x}}+{\widehat{\mathit{x}}}^T\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\beta }\widehat{\mathit{x}}+{\widehat{\mathit{x}}}^T\mathit{W}\frac{\mathrm{d}\widehat{\mathit{x}}}{\mathrm{d}\beta }=0.$$

(13)

Together with Equation (6), this renders the linear system of equations for the eigenvalue and eigenvector derivatives

$$\left[\begin{array}{cc}\frac{\partial \mathit{G}\left(s\right)}{\partial s}\widehat{\mathit{x}}& \mathit{G}\left(s\right)\\ 0& 2{\widehat{\mathit{x}}}^T\mathit{W}\end{array}\right]\left[\begin{array}{c}\frac{\mathrm{d}s}{\mathrm{d}\beta }\\ \frac{\mathrm{d}\widehat{\mathit{x}}}{\mathrm{d}\beta }\end{array}\right]=\left[\begin{array}{c}-\frac{\partial \mathit{G}\left(s\right)}{\partial \beta }\widehat{\mathit{x}}\\ -{\widehat{\mathit{x}}}^T\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\beta }\widehat{\mathit{x}}\end{array}\right],$$

(14)

where the bilinear forms in Equation (13) are summarized for a symmetric weighting matrix $\mathit{W}$. The right-hand side of this system of equations gathers the derivatives of the aeroelastic model matrices, which allows us to solve the system of equations for multiple right-hand sides and, thus, multiple design parameters. Alternatively, as discussed in [5], the adjoint system of equations can be solved for an indicator vector as the right-hand side selecting one component of the solution vector, such as the eigenvalue derivative. Subsequently, the eigenvalue derivatives result efficiently from the vector product of this solution with an arbitrary number of right-hand sides of Equation (14). For the eigenvalue derivative, this reads as

$$\frac{\mathrm{d}s}{\mathrm{d}\beta }=\left[\begin{array}{cccc}1& 0& \cdots & 0\end{array}\right]\left[\begin{array}{c}\frac{\mathrm{d}s}{\mathrm{d}\beta }\\ \frac{\mathrm{d}\widehat{\mathit{x}}}{\mathrm{d}\beta }\end{array}\right]={\left[\begin{array}{c}-\frac{\partial \mathit{G}\left(s\right)}{\partial \beta }\widehat{\mathit{x}}\\ -{\widehat{\mathit{x}}}^T\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\beta }\widehat{\mathit{x}}\end{array}\right]}^T\left({\left({\left[\begin{array}{cc}\frac{\partial \mathit{G}\left(s\right)}{\partial s}\widehat{\mathit{x}}& \mathit{G}\left(s\right)\\ 0& 2{\widehat{\mathit{x}}}^T\mathit{W}\end{array}\right]}^T\right)}^{-1}\left[\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right]\right).$$

For the aeroelastic eigenproblem in modal coordinates, the linear system of equations is obtained according to Equation (14), where the terms are replaced by their modal counterparts and the modal eigenvectors are normalized by the symmetric weighting matrix $\tilde{\mathit{W}}$, resulting in

$$\left[\begin{array}{cc}\frac{\partial \tilde{\mathit{G}}\left(s\right)}{\partial s}\widehat{\mathit{q}}& \tilde{\mathit{G}}\left(s\right)\\ 0& 2{\widehat{\mathit{q}}}^T\tilde{\mathit{W}}\end{array}\right]\left[\begin{array}{c}\frac{\mathrm{d}s}{\mathrm{d}\beta }\\ \frac{\mathrm{d}\widehat{\mathit{q}}}{\mathrm{d}\beta }\end{array}\right]=\left[\begin{array}{c}-\frac{\partial \tilde{\mathit{G}}\left(s\right)}{\partial \beta }\widehat{\mathit{q}}\\ -{\widehat{\mathit{q}}}^T\frac{\mathrm{d}\tilde{\mathit{W}}}{\mathrm{d}\beta }\widehat{\mathit{q}}\end{array}\right].$$

(15)

The terms containing the modal matrix derivative—see Equation (11)—are included in the right-hand side by $\partial \tilde{\mathit{G}}\left(s\right)/\partial \beta$.

3.3. Aeroelastic Sensitivities for the p-k Method

As described in Section 2.3, the aerodynamic transfer function matrix of the p-k method is not analytic. Therefore, the partial derivative with respect to the Laplace variable s does not exist, which is a requirement for the aeroelastic derivatives determined in Section 3.1. Consequently, the aeroelastic sensitivities are derived for nonanalytic aerodynamic transfer function matrices by splitting the argument s into its real and imaginary part. Considering the dynamic aeroelastic matrix $\mathit{G}\left(s\right)$ as a function of $\omega$ and $\sigma$ renders the differentiation of Equation (1) to

$$\begin{array}{c}\hfill \frac{\partial \mathit{G}\left(s\right)}{\partial \beta }\widehat{\mathit{x}}+\frac{\partial \mathit{G}\left(s\right)}{\partial \sigma }\widehat{\mathit{x}}\frac{\mathrm{d}\sigma }{\mathrm{d}\beta }+\frac{\partial \mathit{G}\left(s\right)}{\partial \omega }\widehat{\mathit{x}}\frac{\mathrm{d}\omega }{\mathrm{d}\beta }+\mathit{G}\left(s\right)\frac{\mathrm{d}\widehat{\mathit{x}}}{\mathrm{d}\beta }=0,\end{array}$$

(16)

where s still serves as the argument summarizing $\sigma$ and $\omega$. In contrast to Equation (6), the partial derivatives of $\mathit{G}\left(s\right)$ to $\sigma$ and $\omega$ appear together with the derivatives $\mathrm{d}\sigma /\mathrm{d}\beta$ and $\mathrm{d}\omega /\mathrm{d}\beta$, which correspond to the real and imaginary part of the eigenvalue derivative. The partial derivatives of $\mathit{G}\left(s\right)$ are obtained as

$$\begin{array}{c}\hfill \frac{\partial \mathit{G}\left(s\right)}{\partial \sigma }=2s\mathit{M}+\mathit{D}-\frac{\partial \mathit{A}\left(s\right)}{\partial \sigma },\end{array}$$

(17)

$$\begin{array}{c}\hfill \frac{\partial \mathit{G}\left(s\right)}{\partial \omega }=\mathrm{i}2s\mathit{M}+\mathrm{i}\mathit{D}-\frac{\partial \mathit{A}\left(s\right)}{\partial \omega },\end{array}$$

(18)

where the partial derivatives of $\mathit{A}\left(s\right)$ with respect to $\sigma$ and $\omega$ are found. This allows us to evaluate the nonanalytic formulations of the aerodynamic transfer function matrix, where the partial derivatives $\partial \mathit{A}\left(s\right)/\partial \sigma$ and $\partial \mathit{A}\left(s\right)/\partial \omega$ are differently provided.

For the p-k method, the partial derivative $\partial \mathit{A}\left(s\right)/\partial \sigma$ is zero because of the constant approximation for ${\sigma }^{*}\ne 0$, and the partial derivative $\partial \mathit{A}\left(s\right)/\partial \omega$ is the partial derivative of ${\mathit{A}}_k\left({\omega }^{*}\right)$, including the partial derivative of the reduced frequency ${\omega }^{*}$. The same applies to the partial derivative with respect to $\beta$, which is obtained analogously to Equation (10).

In summary, the partial derivatives in the case of the p-k method result in

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}\left(s\right)}{\partial \sigma }& \hfill =\frac{L}{V}\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\sigma }^{*}}=\mathbf{0},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}\left(s\right)}{\partial \omega }& \hfill =\frac{L}{V}\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\omega }^{*}},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}\left(s\right)}{\partial \beta }& \hfill =\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial \beta }+\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\omega }^{*}}\frac{\partial (\omega L/V)}{\partial \beta }.\hfill \end{array}$$

(21)

In order to solve Equation (16) for the eigensensitivities, the eigenvector normalization from Equation (12) is required as a second equation. Subsequently, both equations are split into their real and imaginary parts, resulting in four equations for the unknown eigensensitivities $\mathrm{d}\sigma /\mathrm{d}\beta$ and $\mathrm{d}\omega /\mathrm{d}\beta$, as well as the real and imaginary counterparts of $\mathrm{d}\widehat{\mathit{x}}/\mathrm{d}\beta$, namely $\Re \left(\mathrm{d}\widehat{\mathit{x}}/\mathrm{d}\beta \right)$ and $\Im \left(\mathrm{d}\widehat{\mathit{x}}/\mathrm{d}\beta \right)$. By expanding the complex-valued multiplications with the eigenvector derivative into real and imaginary components, the four equations establish the direct method for the real and imaginary parts of the eigensensitivities:

$$\left[\begin{array}{cccc}\Re \left(\frac{\partial \mathit{G}\left(s\right)}{\partial \sigma }\widehat{\mathit{x}}\right)& \Re \left(\frac{\partial \mathit{G}\left(s\right)}{\partial \omega }\widehat{\mathit{x}}\right)& \Re \left(\mathit{G}\left(s\right)\right)& -\Im \left(\mathit{G}\left(s\right)\right)\\ \Im \left(\frac{\partial \mathit{G}\left(s\right)}{\partial \sigma }\widehat{\mathit{x}}\right)& \Im \left(\frac{\partial \mathit{G}\left(s\right)}{\partial \omega }\widehat{\mathit{x}}\right)& \Im \left(\mathit{G}\left(s\right)\right)& \Re \left(\mathit{G}\left(s\right)\right)\\ 0& 0& 2\Re \left({\widehat{\mathit{x}}}^T\mathit{W}\right)& -2\Im \left({\widehat{\mathit{x}}}^T\mathit{W}\right)\\ 0& 0& 2\Im \left({\widehat{\mathit{x}}}^T\mathit{W}\right)& 2\Re \left({\widehat{\mathit{x}}}^T\mathit{W}\right)\end{array}\right]\left[\begin{array}{c}\frac{\mathrm{d}\sigma }{\mathrm{d}\beta }\\ \frac{\mathrm{d}\omega }{\mathrm{d}\beta }\\ \Re \left(\frac{\mathrm{d}\widehat{\mathit{x}}}{\mathrm{d}\beta }\right)\\ \Im \left(\frac{\mathrm{d}\widehat{\mathit{x}}}{\mathrm{d}\beta }\right)\end{array}\right]=\left[\begin{array}{c}-\Re \left(\frac{\partial \mathit{G}\left(s\right)}{\partial \beta }\widehat{\mathit{x}}\right)\\ -\Im \left(\frac{\partial \mathit{G}\left(s\right)}{\partial \beta }\widehat{\mathit{x}}\right)\\ -\Re \left({\widehat{\mathit{x}}}^T\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\beta }\widehat{\mathit{x}}\right)\\ -\Im \left({\widehat{\mathit{x}}}^T\frac{\mathrm{d}\mathit{W}}{\mathrm{d}\beta }\widehat{\mathit{x}}\right)\end{array}\right].$$

(22)

Taken together, Equation (22) is equivalent to Equation (14) if the aerodynamic transfer function matrix is analytic. However, Equation (22) enables the separate evaluation of the partial derivatives with respect to the real and imaginary part of s, as required in the nonanalytic case. In addition, the linear system of equations for the eigensensitivities employing modal coordinates is obtained by substituting the terms in Equation (22) with their modal counterparts analogous to Equation (15).

3.4. Aeroelastic Sensitivities for the g Method

For the g method, similar to the p-k method, the aeroelastic derivatives of Section 3.1 are not applicable, since the g method is only analytic at $\sigma =0$. Therefore, it is necessary to separately provide the partial derivatives of the aerodynamic transfer function matrix and to employ the direct method stated in Equation (22) in order to solve for the eigensensitivities.

This requires the partial derivatives of the aerodynamic transfer function matrix appearing in Equations (7), (17) and (18). They are obtained by differentiating the aerodynamic transfer function matrix of the g method, ${\mathit{A}}_g\left(s^{*}\right)$, defined in Equation (5), and by accounting for the chained partial derivative of the reduced quantities ${\sigma }^{*}$, ${\omega }^{*}$ to the dimensional quantities $\sigma$, $\omega$, respectively.

The partial derivative of ${\mathit{A}}_g\left(s^{*}\right)$ to ${\sigma }^{*}$ results in the linear factor $-\mathrm{i}\partial {\mathit{A}}_k\left({\omega }^{*}\right)/\partial {\omega }^{*}$ for ${\sigma }^{*}$ in ${\mathit{A}}_g\left(s^{*}\right)$. For the partial derivative with respect to ${\omega }^{*}$, this factor yields a second-order partial derivative of ${\mathit{A}}_k\left({\omega }^{*}\right)$ to ${\omega }^{*}$, which, together with the first-order partial derivative, renders the partial derivative of ${\mathit{A}}_g\left(s^{*}\right)$ to ${\omega }^{*}$. Similarly, the partial derivative with respect to $\beta$ includes the mixed partial derivative of ${\mathit{A}}_k\left({\omega }^{*}\right)$ to ${\omega }^{*}$ and $\beta$. Moreover, by the chain rule, it includes the derivatives $\partial {\mathit{A}}_g\left(s^{*}\right)/\partial {\sigma }^{*}$ and $\partial {\mathit{A}}_g\left(s^{*}\right)/\partial {\omega }^{*}$, together with the partial derivatives of ${\sigma }^{*}$ and ${\omega }^{*}$ to $\beta$, which only apply if $\beta$ refers to L or V.

In summary, the partial derivatives in the case of the g method result in

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}\left(s\right)}{\partial \sigma }& \hfill =\frac{L}{V}\frac{\partial {\mathit{A}}_g\left(s^{*}\right)}{\partial {\sigma }^{*}}=-\mathrm{i}\frac{L}{V}\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\omega }^{*}},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}\left(s\right)}{\partial \omega }& \hfill =\frac{L}{V}\frac{\partial {\mathit{A}}_g\left(s^{*}\right)}{\partial {\omega }^{*}}=\frac{L}{V}\left(\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\omega }^{*}}-\mathrm{i}\frac{{\partial }^2{\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {{\omega }^{*}}^2}{\sigma }^{*}\right),\hfill \end{array}$$

(24)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial \mathit{A}\left(s\right)}{\partial \beta }& =\frac{\partial {\mathit{A}}_g\left(s^{*}\right)}{\partial \beta }=\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial \beta }-\mathrm{i}\frac{{\partial }^2{\mathit{A}}_k\left({\omega }^{*}\right)}{\partial \beta \partial {\omega }^{*}}{\sigma }^{*}\hfill \\ \hfill & +\left(\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\omega }^{*}}-\mathrm{i}\frac{{\partial }^2{\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {{\omega }^{*}}^2}{\sigma }^{*}\right)\frac{\partial (\omega L/V)}{\partial \beta }-\mathrm{i}\frac{{\partial }^2{\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {{\omega }^{*}}^2}\frac{\partial (\sigma L/V)}{\partial \beta }\hfill \end{array}\end{array}$$

(25)

3.5. Sensitivities of the Modal Matrix

The sensitivities of the modal matrix require determination of the derivatives of the structural mode shapes in wind-off condition and neglecting any structural damping terms with respect to a design parameter $\beta$. This is equivalent to finding the sensitivities of the structural eigenproblem providing the mode shapes.

The structural mode shapes $\mathit{\phi }$ are the eigenvectors of the generalized eigenproblem

$$\begin{array}{cc}\hfill \mathit{K}\mathit{\phi }& \hfill =\lambda \mathit{M}\mathit{\phi }\hfill \end{array}$$

(26)

resulting from the structural dynamics equations for harmonic motion. The structural eigenfrequencies are given by $\sqrt{\lambda }$ and the structural eigenvectors are mass-normalized by

$$\begin{array}{c}\hfill {\mathit{\phi }}^T\mathit{M}\mathit{\phi }=1,\end{array}$$

(27)

which renders the modal mass matrix $\tilde{\mathit{M}}={\mathsf{\Phi }}^T\mathit{M}\mathsf{\Phi }$ as the identity matrix.

Differentiating Equations (26) and (27) with respect to $\beta$ allow us to solve for the structural eigenvalue and eigenvector sensitivities with the direct method:

$$\begin{array}{cc}\hfill (\frac{\mathrm{d}\mathit{K}}{\mathrm{d}\beta }-\lambda \frac{\mathrm{d}\mathit{M}}{\mathrm{d}\beta })\mathit{\phi }-\mathit{M}\mathit{\phi }\frac{\mathrm{d}\lambda }{\mathrm{d}\beta }+(\mathit{K}-\lambda \mathit{M})\frac{\mathrm{d}\mathit{\phi }}{\mathrm{d}\beta }& \hfill =\mathbf{0},\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill 2{\mathit{\phi }}^T\mathit{M}\frac{\mathrm{d}\mathit{\phi }}{\mathrm{d}\beta }+{\mathit{\phi }}^T\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\beta }\mathit{\phi }& \hfill =0.\hfill \end{array}$$

(29)

The latter equation requires the mass matrix to be symmetric in order to combine both bilinear forms containing the eigenvector derivative. Moreover, it enforces the condition imposed by Equation (27) with a changing parameter $\beta$. This constitutes a linear system of equations for the eigenvalue and eigenvector sensitivities

$$\begin{array}{cc}\hfill \left[\begin{array}{cc}-\mathit{M}\mathit{\phi }& \mathit{K}-\lambda \mathit{M}\\ 0& 2{\mathit{\phi }}^T\mathit{M}\end{array}\right]\left[\begin{array}{c}\frac{\mathrm{d}\lambda }{\mathrm{d}\beta }\\ \frac{\mathrm{d}\mathit{\phi }}{\mathrm{d}\beta }\end{array}\right]& =\left[\begin{array}{c}-(\frac{\mathrm{d}\mathit{K}}{\mathrm{d}\beta }-\lambda \frac{\mathrm{d}\mathit{M}}{\mathrm{d}\beta })\mathit{\phi }\\ -{\mathit{\phi }}^T\frac{\mathrm{d}\mathit{M}}{\mathrm{d}\beta }\mathit{\phi }\end{array}\right],\hfill \end{array}$$

(30)

which is then solved for each eigenvector and corresponding eigenvalue included in the modal matrix. Subsequently, the derivative of the modal matrix is given by

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\mathsf{\Phi }}{\mathrm{d}\beta }& \hfill =\left[\begin{array}{ccc}\cdots & \frac{\mathrm{d}{\mathit{\phi }}_j}{\mathrm{d}\beta }& \cdots \end{array}\right].\hfill \end{array}$$

(31)

4. Results

4.1. Verification of Aeroelastic Sensitivities

The verification of the analytically derived aeroelastic sensitivities is conducted for the two-degree-of-freedom (DOF) typical section model by analyzing the convergence error to finite difference approximations. The investigated typical section model employs Theodorsen’s aerodynamic airfoil theory [19] with the generalized Theodorsen function developed by Edwards [20,21]. This allows us to derive the analytical partial derivatives of the aerodynamic transfer function matrix presented in Section 3.1 and, moreover, in Section 3.3 and Section 3.4 for the p-k and g methods. The analytical derivatives of the aerodynamic transfer function matrix and the generalized Theodorsen function are derived in Appendix A and Appendix B. Additionally, two degrees of freedom facilitate the employment of the modal approach without modal truncation in order to verify the derivatives developed for the aeroelastic eigenproblem in modal coordinates—see Section 3.1 and Section 3.5—by comparison to the solutions obtained in physical coordinates.

The 2-DOF typical section model constitutes a two-dimensional airfoil section with a plunge and pitch degrees of freedom, denoted with subscripts h and $\alpha$, in incompressible flow; see Figure 1. The air loads are described by the lift force and the aerodynamic moment around the elastic axis, where the section model is elastically restrained with a vertical stiffness $k_h$ and a torsional stiffness $k_{\alpha }$. Therefore, the nonlinear aeroelastic eigenproblem for the typical section model results in

$$\begin{array}{cc}\hfill \left(s^2\left[\begin{array}{cc}m& S_{\alpha }\\ S_{\alpha }& I_{\alpha }\end{array}\right]+\left[\begin{array}{cc}{k}_h& 0\\ 0& k_{\alpha }\end{array}\right]-\mathit{A}\left(s^{*}\right)\right)\left[\begin{array}{c}{\widehat{x}}_h\\ {\widehat{x}}_{\alpha }\end{array}\right]& \hfill =\left[\begin{array}{c}0\\ 0\end{array}\right],\hfill \end{array}$$

(32)

where the mass matrix includes the total mass m, the first moment of mass $S_{\alpha }$ around the elastic axis and the moment of inertia $I_{\alpha }$. The aerodynamic transfer function matrix $\mathit{A}\left(s^{*}\right)$ provides the negative lift force and the aerodynamic moment in the Laplace domain from the amplitudes of motion ${\widehat{x}}_h$, ${\widehat{x}}_{\alpha }$. Based on Theodorsen’s aerodynamic airfoil theory [19], commonly described in the aeroelastic literature [22,23], the unsteady aerodynamic forces are extended for growing and decaying oscillations by employing the generalized Theodorsen function $C\left(s^{*}\right)$ [20]. Expressed in dimensional form with reference to the unit span of the section, the aerodynamic transfer function matrix is given by

$$\begin{array}{c}\mathit{A}\left(s^{*}\right)=\rho V^2\pi ({s^{*}}^2\left[\begin{array}{cc}-1& eb\\ eb& -(\frac{1}{8}+e^2)b^2\end{array}\right]\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+s^{*}\left[\begin{array}{cc}-2C\left(s^{*}\right)& \left(-1-2C\left(s^{*}\right)(\frac{1}{2}-e)\right)b\\ 2C\left(s^{*}\right)(\frac{1}{2}+e)b& (\frac{1}{2}-e)\left(2C\left(s^{*}\right)(\frac{1}{2}+e)-1\right)b^2\end{array}\right]\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left[\begin{array}{cc}0& -2C\left(s^{*}\right)b\\ 0& 2C\left(s^{*}\right)(\frac{1}{2}+e)b^2\end{array}\right]),\end{array}$$

(33)

which includes the dynamic pressure $\frac{1}{2}\rho V^2$ and the area per unit span $2b$. The reference length L for the reduced frequencies is the half chord length b. The elastic axis is positioned relative to the mid-point by the nondimensional parameter e. The generalized Theodorsen function $C\left(s^{*}\right)$ is provided in Equation (A5) in Appendix B.

The aerodynamic transfer function matrix for the p-k method ${\mathit{A}}_k\left({\omega }^{*}\right)$ is found by evaluation of Equation (33) with ${\sigma }^{*}=0$, neglecting the real part of the reduced complex-valued frequency $s^{*}$; see Equation (4). For the g method, in Equation (5), the additionally required partial derivative $\partial {\mathit{A}}_k\left({\omega }^{*}\right)/\partial {\omega }^{*}$ is obtained by applying the analytic property of $\mathit{A}\left(s^{*}\right)$ at ${\sigma }^{*}=0$, relating

$$\frac{\partial {\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\omega }^{*}}=\frac{\partial \mathit{A}\left(\mathrm{i}{\omega }^{*}\right)}{\partial {\omega }^{*}}=\mathrm{i}\frac{\partial \mathit{A}\left(\mathrm{i}{\omega }^{*}\right)}{\partial s^{*}},$$

where $\partial \mathit{A}\left(s^{*}\right)/\partial s^{*}$, given in Equation (A1) in Appendix A, requires the derivative of the generalized Theodorsen function derived in Appendix B.

The employed parameter values for the section model are listed in

Table 1. Parameter values for the typical section model. The structural parameters m, $S_{\alpha }$, $I_{\alpha }$, $k_h$ and $k_{\alpha }$ are given per unit span, $\rho$ is the air density, e the relative position of the elastic axis and b the half chord length.

m (kg/m)	${\mathit{S}}_{\mathit{\alpha }}$ (kg)	${\mathit{I}}_{\mathit{\alpha }}$ (kg·m)	${\mathit{k}}_{\mathit{h}}$ (N/m2)	${\mathit{k}}_{\mathit{\alpha }}$ (N)	b (m)	e	$\mathit{\rho }$ (kg/m3)
292.4823	73.1206	113.482	9.1396 $\times {10}^5$	4.1965 $\times {10}^5$	1	−0.15	1.225

.

In Figure 2, the eigenvalue solutions for the section model obtained by the p-k method, the g method and GAAM are shown over a velocity sweep from 0 m/s to 300 m/s. Overall, the three methods yield similar results. From Figure 2b,c, the flutter onset is found for the second eigenvalue $s_2$, where the real part ${\sigma }_2$ is zero at the freestream velocity $V=212.2$ m/s. As expected, the flutter onset is identically identified by all three methods. However, for other velocities where $\sigma \ne 0$, the eigenvalue solutions exhibit slight variations. For the p-k method, the differences are more pronounced, while the g method deviates little from GAAM. In fact, for high velocities where $\sigma$ sufficiently increases from zero, small deviations for the g method are observed. These results confirm the main distinction between the aerodynamic damping approximations in which the real part of the Laplace variable is included to different degrees.

At each velocity, the derivatives of the eigenvalue solutions with respect to a design parameter are obtained by employing the direct method; see Section 3.2. The design parameter for the verification of the derived equations is selected to be the half chord length b. Hence, the quantities of interest are how the eigenvalues change with changing b. Varying the half chord length implies a change in shape, requiring the partial derivative of the aerodynamic transfer function matrix $\partial \mathit{A}\left(s\right)/\partial \beta$. Moreover, since the half chord length is the reference length, the additional terms regarding the partial derivative of the reduced frequency with respect to L are included. Therefore, the design parameter b is well suited to verify the derived eigenvalue derivatives. However, the structural derivatives with respect to b are neglected in the following verification and the structural parameters are considered to be independent of b.

In order to obtain the eigenvalue derivatives, GAAM requires the partial derivatives $\partial \mathit{A}\left(s\right)/\partial s^{*}$ and $\partial \mathit{A}\left(s^{*}\right)/\partial b$, given by Equations (A1) and (A2) in Appendix A. The partial derivative $\partial {\mathit{A}}_k\left({\omega }^{*}\right)/\partial b$, for the p-k and g methods, is the derivative $\partial \mathit{A}\left(\mathrm{i}{\omega }^{*}\right)/\partial b$, where the real part ${\sigma }^{*}$ is neglected. The g method additionally requires the second-order partial derivatives ${\partial }^2{\mathit{A}}_k\left({\omega }^{*}\right)/{\partial {\omega }^{*}}^2$ and ${\partial }^2{\mathit{A}}_k\left({\omega }^{*}\right)/\partial {\omega }^{*}\partial b$, as discussed in Section 3.4. These are again obtained from the analytic property of $\mathit{A}\left(s^{*}\right)$, resulting in

$$\begin{array}{cc}\hfill \frac{{\partial }^2{\mathit{A}}_k\left({\omega }^{*}\right)}{{\partial {\omega }^{*}}^2}& \hfill =\frac{{\partial }^2\mathit{A}\left(\mathrm{i}{\omega }^{*}\right)}{{\partial {\omega }^{*}}^2}=-\frac{{\partial }^2\mathit{A}\left(\mathrm{i}{\omega }^{*}\right)}{{\partial s^{*}}^2},\hfill \\ \hfill \frac{{\partial }^2{\mathit{A}}_k\left({\omega }^{*}\right)}{\partial {\omega }^{*}\partial b}& \hfill =\frac{{\partial }^2\mathit{A}\left(\mathrm{i}{\omega }^{*}\right)}{\partial {\omega }^{*}\partial b}=\mathrm{i}\frac{{\partial }^2\mathit{A}\left(\mathrm{i}{\omega }^{*}\right)}{\partial s^{*}\partial b},\hfill \end{array}$$

where the second-order partial derivatives of $\mathit{A}\left(s^{*}\right)$ are evaluated at ${\sigma }^{*}=0$.

For the verification, the half chord length of the section model is changed by step sizes $\Delta b$ ranging from 1 × 10${}^{-6}$ m to $\Delta b_{max}=0.05$ m. Subsequently, for a given velocity, the eigenvalue solutions for each modified section model are obtained by solving the nonlinear eigenproblem. Figure 3 displays the change in the eigenvalues for changing half chord length by $\pm \Delta b$ at a freestream velocity $V=209.6$ m/s. The velocity is close to the flutter onset and, thus, the eigenvalue $s_2$ is close to the imaginary axis where $\sigma =0$. To illustrate the results, at this velocity, the second eigenvalue would become less stable with increasing half chord length and eventually become unstable. Additionally, the predicted linear changes in the eigenvalues by the eigenvalue derivatives are shown as lines attached to the eigenvalues of the unchanged section model. The depicted lines display the scaled eigenvalue derivatives by multiplication with $\pm \Delta b_{max}$. From visual judgement, the resulting tangent lines show very good agreement with the slopes of the nonlinear curves at the solution of the unchanged model. In addition to the visual comparison, the numerical values of the obtained eigenvalue derivatives are given in

Table 2. Eigenvalue derivatives with respect to the half chord length b for the typical section model employing the p-k method, the g method and GAAM at $V=209.6$ m/s.

	$\mathbf{d}{\mathit{s}}_1/\mathbf{d}\mathit{b}$ (rad/(m·s))	$\mathbf{d}{\mathit{s}}_2/\mathbf{d}\mathit{b}$ (rad/(m·s))
p-k method	$-44.180995-9.676179\mathrm{i}$	$31.725084-13.803641\mathrm{i}$
g method	$-54.545970-0.113813\mathrm{i}$	$45.695638-15.883591\mathrm{i}$
GAAM	$-54.064094+0.513874\mathrm{i}$	$45.905266-16.045078\mathrm{i}$

.

By comparing the result between the p-k method, g method and GAAM in Figure 3, the difference in the eigenvalue solution for the unchanged section model is noted, as is also shown in Figure 2. This difference is more pronounced at $s_1$ with greater distance from the imaginary axis. Hence, the nonlinearly changing eigenvalues exhibit differences, which are the greatest for the p-k method, since the unchanged eigenvalue solution already shows the greatest offset. The g method shows very little difference to GAAM for the nonlinearly and linearly obtained second eigenvalues $s_2$. For the first eigenvalue, a small offset in the unchanged eigenvalue is observed, with little differences in the nonlinearly obtained eigenvalues for greater values of $+\Delta b$. Although the unchanged eigenvalues for the g method and GAAM are very close, a small but visible difference in the visualized eigenvalue derivatives at the first eigenvalue is observed, as reported in Table 2, reflecting the difference in the underlying aerodynamic damping approximation. This is one of the main results of the current work, pointing out the effect of the aerodynamic damping approximation on the aeroelastic eigensensitivities, even though the predicted eigensolution to the flutter equation is very similar.

Finally, the derived analytical eigenvalue derivatives are verified by the analysis of the relative error to forward finite difference approximations. Considering the eigenvalues as a function of $\beta$ renders the finite difference approximations using nonlinearly obtained eigenvalues evaluated at $\beta +\Delta \beta$. Subsequently, the relative error is given by

$${\eta }_{\mathrm{d}s/\mathrm{d}\beta }=|\frac{\mathrm{d}s/\mathrm{d}\beta -\left(s\right(\beta +\Delta \beta )-s(\beta \left)\right)/\Delta \beta }{\mathrm{d}s/\mathrm{d}\beta }| = |\frac{O\left(\Delta \beta \right)}{\mathrm{d}s/\mathrm{d}\beta }|,$$

(34)

where the latter equality results from Taylor’s theorem, leading to a first-order approximation error if the evaluated derivatives are the analytical derivatives. Consequently, the relative error converges linearly with decreasing $\Delta b$ for the analytical derivatives.

Figure 4 displays the relative errors of Equation (34) for the eigenvalue derivatives at two freestream velocities $V=209.6$ m/s and $V=241.2$ m/s. For the three investigated aerodynamic damping approximations, the relative errors separated in real and imaginary parts exhibit an exponential decline for the decreasing step size $\Delta b$. Examining the slopes shows a first-order reduction in the relative error, which is in line with the stated condition for analytical derivatives in Equation (34). For very small step sizes $\Delta b$, the relative error deviates from the first-order reduction because of the increasing numerical error in the finite difference approximations due to numerical cancellation. From these results, it is concluded that the obtained eigenvalue derivatives are in agreement with the finite difference approximations, verifying the derived equations for the eigenvalue sensitivities.

To this end, it remains to verify the eigenvalue sensitivities employing modal coordinates. For this, it is sufficient to compare the eigenvalue derivatives obtained in modal coordinates without modal truncation to the verified solutions obtained in physical coordinates. Figure 5 shows this comparison for the design parameter b employing GAAM. In Figure 5a, the eigenvalue solutions are displayed for a velocity sweep from 0 m/s to 300 m/s. Figure 5b,c show the real and imaginary parts of the eigenvalue derivatives at each velocity. The eigenvalue derivatives are obtained by the direct method in physical coordinates—see Equation (14)—and in modal coordinates—see Equation (15). No differences are observed in Figure 5, verifying the equality of the solutions of the direct method in physical and modal coordinates without modal truncation.

4.2. Aeroelastic Sensitivities of the AGARD 445.6 Wing

In this section, the derived aeroelastic sensitivities are applied to the AGARD 445.6 weakened wing model for varying wing sweep. The AGARD 445.6 weakened wing model is a widely known aeroelastic test case [24,25]. It is a swept wing clamped at the wing root with a quarter chord sweep of $\Lambda =45$${}^{\circ }$, a root chord length of $c=0.5578$ m, a span of 0.762 m and a taper ratio of 0.66. The model comprises a finite element model (FEM) with 200 shell elements and a doublet lattice model (DLM) with 200 boxes, which both are obtained from [26] for a Mach number of 0.9 and air density $\rho$ = 0.31 kg/m³.

The commercial solver MSC Nastran [7] is employed in order to extract the mass matrix $\mathit{M}$ and stiffness matrix $\mathit{K}$ as well as the modal matrix $\mathsf{\Phi }$ by performing a linear modal analysis. In this model, the first four structural mode shapes are considered, which are the first wing bending and first wing torsion, followed by the second wing bending and second wing torsion. The structural eigenfrequencies are listed in

Table 3. Parameter values and structural eigenfrequencies for the AGARD 445.6 weakened model. The reference length is the half chord length at the wing root, $L=c/2$.

Mach Number	$\mathit{\rho }$ (kg/m3)	$\Lambda$ (${}^{\circ }$)	c (m)	L (m)	Eigenfrequencies (rad/s)
0.9	0.31	45	0.5578	0.2789	59.4, 294.4, 310.7, 597.5

together with the parameter values. Furthermore, the aerodynamic transfer function matrices ${\mathit{A}}_k\left({\omega }^{*}\right)$ employing DLM are extracted from MSC Nastran, sampled at 17 reduced frequencies ${\omega }^{*}$ ranging from 0.001 to 5, clustered densely in the region up to 1. An interpolation method based on piecewise cubic functions provides the evaluation of ${\mathit{A}}_k\left({\omega }^{*}\right)$ at additional frequencies, which are required in obtaining the eigenvalue solutions. Moreover, the interpolation method enables the efficient evaluation of the derivatives of ${\mathit{A}}_k\left({\omega }^{*}\right)$ with respect to ${\omega }^{*}$, which are, therefore, continuous up to the second order. The aerodynamic transfer function matrices include the spline matrix, which is required for mapping the structural displacement onto the aerodynamic grid as well as the aerodynamic forces onto the structural grid using infinite plate spline interpolation [27]. Taken together, the extracted matrices allow us to solve the nonlinear aeroelastic eigenproblem in modal coordinates; cf. Equation (2).

In Figure 6, the eigenvalue solutions for the AGARD 445.6 weakened model obtained by the p-k method and the g method are shown over a velocity sweep from 0 m/s to 350 m/s. The flutter onset is found by both methods for the first eigenvalue $s_1$ at the freestream velocity $V=193.9$ m/s, as shown in Figure 6b,c. For smaller freestream velocities, the p-k method and the g method yield very similar solutions. The solutions start to differ for higher freestream velocities, where the real part of the eigenvalues $\sigma$ increases. This is most pronounced for the first and second eigenvalues, $s_1$ and $s_2$, which exhibit the most damping and, thus, significant deviations between both methods are observed for freestream velocities greater than the flutter onset velocity.

The sweep angle of the AGARD 445.6 weakened model is varied by applying a shear transformation for a step size $\Delta \Lambda$ without affecting the other parameters, such as root chord length, span and taper ratio. The sweep angle is changed for step sizes ranging from 4.67 × 10${}^{-8}$${}^{\circ }$ to $\Delta {\Lambda }_{max}=4.67{}^{\circ }$. For each step size, MSC Nastran is employed for acquiring $\mathit{M}$, $\mathit{K}$, $\mathsf{\Phi }$ and ${\mathit{A}}_k\left({\omega }^{*}\right)$, as described previously. Subsequently, the nonlinear aeroelastic eigenproblem is solved for each sweep angle in order to obtain the nonlinearly changing eigenvalues for varying sweep. Furthermore, the extracted mass matrix $\mathit{M}$ and stiffness matrix $\mathit{K}$ are used for forward finite difference approximations of the model sensitivities $\mathrm{d}\mathit{M}/\mathrm{d}\beta$ and $\mathrm{d}\mathit{K}/\mathrm{d}\beta$. From these sensitivities, the sensitivity of the modal matrix $\mathrm{d}\mathsf{\Phi }/\mathrm{d}\beta$ is obtained following Section 3.5. Additionally, the sensitivity of the aerodynamic transfer function matrix $\mathrm{d}{\mathit{A}}_k\left({\omega }^{*}\right)/\mathrm{d}\beta$ is determined by forward finite difference approximations for each sampled frequency. With these defined model sensitivities, the eigenvalue sensitivities are solved by employing the direct method for separated real and imaginary parts, as stated in Equation (22) in model coordinates. Hence, the linear system of equations includes the additional terms regarding the derivative of the modal matrix; see Equation (11).

Figure 7 shows the nonlinearly changing eigenvalues as well as the predicted linear change from the eigenvalue sensitivities for varying sweep at the flutter onset velocity. Similar to Section 4.1, the lines representing the eigenvalue sensitivities are defined by the multiplication of the eigenvalue sensitivities with $\pm \Delta {\Lambda }_{max}$. The eigenvalue sensitivities are shown for the model sensitivities obtained for the step size $\Delta \Lambda$ = 5.9 × 10${}^{-5}$${}^{\circ }$, resulting from a convergence study, shown in Figure 8. For the flutter onset velocity, the p-k and g methods yield the same eigenvalue $s_1$ for the unchanged model, as is confirmed in Figure 7a. However, the nonlinearly obtained eigenvalues for changing sweep differ for both methods, which is correctly reflected by the observed difference in the linearly predicted eigenvalues. Hence, the eigenvalue sensitivities are different for the p-k method and the g method, although the eigenvalue solutions are the same and, moreover, have a real part of zero.

For the second eigenvalue $s_2$—see Figure 7b—the difference between both methods results primarily from the offset between the eigenvalue solutions for the unchanged model. As a result, different eigenvalue sensitivities are obtained—see also Figure 8, where numerical values are compared. The difference in the eigenvalue solutions of the third and fourth eigenvalue shown in Figure 7c,d is less pronounced. This is in line with the previous observation that the first and second eigenvalues are more affected by the investigated aerodynamic damping approximations.

In Figure 9, the eigenvalue solutions for changing sweep are shown for the freestream velocity $V=230.7$ m/s, above the flutter onset. In this case, the difference between both methods results from the offset for the first and second eigenvalue, $s_1$ and $s_2$, of the unchanged model. Moreover, the eigenvalue sensitivities as well as the nonlinearly obtained eigenvalues reveal the difference in the change in the aeroelastic stability resulting from the different aerodynamic damping approximations.

5. Conclusions

In this work, the aeroelastic eigensensitivities are derived from the nonlinear aeroelastic eigenproblem accounting for different aerodynamic damping approximations, which are the p-k method, the g method and the generalized aeroelastic analysis method (GAAM). Moreover, the solution method for the aeroelastic eigensensitivities employing the direct method is presented in the context of different aerodynamic damping approximations. The nonlinear aeroelastic eigenproblem is considered in both physical and modal coordinates, leading to additional terms regarding the modal matrix sensitivities in the latter case.

The derived aeroelastic eigensensitivities are verified for the typical section model employing the three presented aerodynamic damping approximations applied to the generalized Theodorsen aerodynamic airfoil theory. The verification is performed by analytical derivatives of the involved sensitivities regarding a shape design parameter. This includes the first and second analytical derivatives of the generalized Theodorsen function. Additionally, the verification includes the approach in modal coordinates.

The influence of the aerodynamic damping approximations on the resulting eigenvalue sensitivities is discussed by means of the typical section model as well as the AGARD 445.6 wing model with varying sweep. The presented results display differences in the obtained eigenvalue sensitivities for the three investigated aerodynamic damping approximations, with the p-k method showing the most pronounced deviations. The differences in the eigenvalue sensitivities can be attributed in part to the offset in the underlying eigensolutions, resulting from increasing absolute values of the damping coefficients. However, the results demonstrate that, for very close or equal eigensolutions, the aerodynamic damping approximations also affect the eigenvalue sensitivities. Moreover, this applies to eigensolutions with zero damping coefficients.

For the application of flutter constraints in gradient-based MDO, the presented results indicate the preference of GAAM and the g method over the p-k method. Although the three methods identify the identical flutter onset where flutter constraints become active, GAAM and the g method are able to provide more accurate sensitivities because of the improved aerodynamic damping. Future work will address the application of the derived approach for the recently presented analytic continuation of DLM [16] as well as the application for high-fidelity CFD methods.