1. Introduction
For over a century, the study of aeroelasticity has been pivotal in many real-life engineering quests (Livne Reference Livne2003; Shubov Reference Shubov2006; Beran, Stanford & Schrock Reference Beran, Stanford and Schrock2017). By investigating the combined effect of aerodynamic, elastic and inertial forces, countless detrimental phenomena have been identified and system designs rectified. Torsional divergence, flutter, limit-cycle oscillations, shock buffet and its structural response, called buffeting, are well-known examples witnessed on aircraft and all can (and have) lead to undesirable reduced overall performance, structural fatigue and worse. Questions on flow stability, specifically in the absence of structural dynamics, have similarly been studied for many decades (Theofilis Reference Theofilis2011). Owing to the invaluable insights gained over the years, the means of transportation that we use every day have been designed to be safer and more efficient with a reduced environmental footprint. Paramount to all these advancements are reliable yet fast and accurate methods to discover how and when these phenomena occur and, crucially, how to limit their impact or prevent them altogether. Numerical analysis, alongside experiments, has proven a powerful tool in this regard.
Shock buffet is an aerodynamic phenomenon characterised by self-excited flow oscillations due to shock-wave/boundary-layer interaction on aircraft wings at transonic speeds. Although it was first identified in the 1960s, a complete multidisciplinary explanation remains elusive, despite progress made continually by experiments and numerical investigations (Tijdeman Reference Tijdeman1977; Jacquin et al. Reference Jacquin, Molton, Deck, Maury and Soulevant2009; Dandois Reference Dandois2016). One of the first theoretical models by Lee (Reference Lee1990) proposed an acoustic feedback loop that sustained the shock oscillations, which was corroborated in experimental studies, although contradictory results and theories have also been put forward (Jacquin et al. Reference Jacquin, Molton, Deck, Maury and Soulevant2009; Feldhusen-Hoffmann et al. Reference Feldhusen-Hoffmann, Statnikov, Klaas and Schröder2018; Moise, Zauner & Sandham Reference Moise, Zauner and Sandham2022). The first global stability study on an aerofoil in turbulent flow linked shock buffet to an unstable eigenmode (Crouch et al. Reference Crouch, Garbaruk, Magidov and Travin2009), soon followed by similar findings on infinite straight and swept wings (Crouch, Garbaruk & Strelets Reference Crouch, Garbaruk and Strelets2019; Paladini et al. Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019a; He & Timme Reference He and Timme2021). Spanwise-localised pockets of shear-layer pulsation synchronised with an outboard-propagating shock oscillation, dubbed buffet cells, that are correlated with the sweep angle, were found to be a defining feature of three-dimensional shock buffet (Iovnovich & Raveh Reference Iovnovich and Raveh2015; Dandois Reference Dandois2016). The first global stability study on a finite-wing aircraft by Timme & Thormann (Reference Timme and Thormann2016) and Timme (Reference Timme2020) corroborated the insight gained from related experimental work (Dandois Reference Dandois2016; Sugioka et al. Reference Sugioka, Koike, Nakakita, Numata, Nonomura and Asai2018; Masini, Timme & Peace Reference Masini, Timme and Peace2020; Sugioka et al. Reference Sugioka, Nakakita, Koike, Nakajima, Nonomura and Asai2021).
The mutual interaction between an elastic wing structure and shock buffet has been investigated previously too. Experimental work on a finite wing found that weak shock-wave/boundary-layer interaction does not affect the structural response strongly, despite large flow field fluctuations (Steimle, Karhoff & Schröder Reference Steimle, Karhoff and Schröder2012). However, with increasing shock strength due to higher free-stream Mach number, strong fluid–structure coupling was observed, despite weaker fluctuations in the flow overall, and the aeroelastic system responded to the unsteady flow excitation. Importantly, it has been shown, for a pitch-plunge (and variants thereof) typical section aerofoil, that the introduction of an elastic structure has the ability to destabilise an otherwise stable flow (Nitzsche et al. Reference Nitzsche, Ringel, Kaiser and Hennings2019). This work also points out limitations of traditional flutter prediction methods, such as the so-called pk-type analysis (Hassig Reference Hassig1971; where p is the eigenvalue and k the reduced frequency), in tracing all relevant modes, prompting the development of more advanced methods. This observation is congruent with our discussion. Expanding global stability analysis by incorporating fluid–structure interaction is indeed an active research area going beyond transonic wings. Novel methods and physics, with unstable eigenmodes originating both in the fluid and structural degrees of freedom, have been explored for various configurations, such as a cylinder with splitter plate and spring-mounted plates and aerofoils (Pfister, Marquet & Carini Reference Pfister, Marquet and Carini2019; Negi, Hanifi & Henningson Reference Negi, Hanifi and Henningson2020; Pfister & Marquet Reference Pfister and Marquet2020; Moulin & Marquet Reference Moulin and Marquet2021). Another important earlier study to be pointed out is that in Lesoinne et al. (Reference Lesoinne, Sarkis, Hetmaniuk and Farhat2001) where the linearised physics were modelled in a similar approach to ours and applied to finite-wing flutter in inviscid flow.
This contribution builds upon previous work by Timme (Reference Timme2020), that studied the fluid-only shock-buffet instability on a finite-wing aircraft, and expands the investigation by Houtman & Timme (Reference Houtman and Timme2021), which introduced an elastic structure. We developed the fluid–structure coupled Jacobian operator, and its integration into the stability framework, to enable the search for dominant fluid modes while at the same time determining the impact of the elastic wing structure. Previous work used the so-called Schur complement method to trace aeroelastic modes describing the wing vibration (Bekemeyer & Timme Reference Bekemeyer and Timme2019; Badcock et al. Reference Badcock, Timme, Marques, Khodaparast, Prandina, Mottershead, Swift, Da Ronch and Woodgate2011). This formulation, which can be shown to be equivalent to legacy tools in conventional flutter analysis, is an efficient way of establishing the flutter boundary while making use of computational fluid dynamics functionality for the aerodynamics. However, it is not suitable for finding fluid instabilities due to the mathematical structure of modelling the coupled physics. The method favoured herein does not require any rearrangement or decomposition of the matrix eigenvalue problem. Choosing a coupled solution approach with an inner–outer iterative structure relying on the shift-invert spectral transformation necessitated the implementation of a novel preconditioner for the inner sparse iterative linear equation solver, while offering speed up compared with discipline-specific block-Jacobi-type preconditioning. In previous work, the aerostructural coupled system has been solved in the context of adjoint gradient computation for multidisciplinary optimisation. Discipline-specific preconditioners from the fluid and structural solvers were reused in Kenway, Kennedy & Martins (Reference Kenway, Kennedy and Martins2014) in a block-Jacobi fashion, noting that discarding the off-diagonal coupling terms allows for easier parallelisation. Block-Jacobi and Gauss–Seidel preconditioners for a three-field formulation were compared in Zhang & Zingg (Reference Zhang and Zingg2018). While this is an intricate discussion, clear performance gains were realised overall when including various coupling terms. A similar strategy of observing the discipline coupling is followed in our work.
This paper continues with a description of the physical models, linearised analyses and numerical details in § 2. The focus is on the coupled discrete Jacobian matrix operator and the adaptation of the inner–outer iterative eigenvalue solver, using the implicitly restarted Arnoldi method from the ARPACK linear algebra library and a bespoke sparse iterative linear solver, all made available in the industrial DLR-TAU code. The adoption of the Sherman–Morrison–Woodbury formula for the parallel inversion of block-arrowhead matrices in deriving a preconditioner for the coupled fluid–structure system and its benefits are discussed, too. Results for the NASA Common Research Model, introduced as a test case in § 3, are scrutinised in § 4, where the impact of the fluid–structure coupling on both the shock-buffet related direct and adjoint eigenmodes as well as the structural sensitivity of the instability to identify the core of the global dynamics, colloquially called wavemaker, are elucidated. A verification test case, together with an assessment of the numerical methods, is outlined in the appendices.
2. Theory and methods
2.1. Physical models of fluid and structure
The starting point of our work is the set of governing equations in semi-discrete form:
\begin{equation} \dot{\boldsymbol{w}} = \mathcal{R}(\boldsymbol{w}) , \end{equation}
where 
$\boldsymbol {w}$ is the state vector comprising two parts, 
$\boldsymbol {w}_{f}$ and 
$\boldsymbol {w}_{s}$, to represent the fluid and structural degrees of freedom, respectively, and 
$\mathcal {R}$ is the corresponding discretised residual operator. The expression 
$\dot {(\,{\cdot }\,)}$ denotes a temporal derivative. Equation (2.1) depends on a large number of parameters which are not explicitly stated herein.
The fluid system is assumed to be governed by the Reynolds-averaged Navier–Stokes (RANS) equations in three-dimensional space coupled with a suitable turbulence model. The fluid state vector 
$\boldsymbol {w}_{f}$ of conservative variables is given by 
$\boldsymbol {w}_{f} = [\rho, \rho \boldsymbol {u}_{f}, \rho E, \rho \tilde {\nu }]^{{\rm T}}$, where 
$\rho$ is the density, 
$\boldsymbol {u}_{f}$ the Cartesian velocity-field vector, 
$E$ the specific total energy and 
$\tilde {\nu }$ the working variable of the turbulence model, when using the negative Spalart–Allmaras model (Allmaras, Johnson & Spalart Reference Allmaras, Johnson and Spalart2012). The governing equations written in conservative integral form in the arbitrary Lagrangian–Eulerian formulation can be stated as
\begin{equation} \frac{\text{d}}{\text{d}t} \int_{\varOmega_{f}(t)} \boldsymbol{w}_{f} \,\text{d}V + \int_{\partial\varOmega_{f}(t)} (\boldsymbol{F} - \boldsymbol{F}_v - \boldsymbol{w}_{f} \dot{\boldsymbol{x}}_{f}) \boldsymbol{\cdot} \boldsymbol{n} \,\text{d}S = \int_{\varOmega_{f}(t)} \boldsymbol{S}\,\text{d}V, \end{equation}
where 
$\varOmega _{f}(t)$ is a time-dependent (due to structural motion) arbitrary control volume enclosed by 
$\partial \varOmega _{f}(t)$, the vectors 
$\boldsymbol {F}$ and 
$\boldsymbol {F}_v$ are the inviscid and viscous fluxes, respectively, 
$\boldsymbol {n}$ is the normal vector to 
$\partial \varOmega _{f}$ and 
$\boldsymbol {S}$ describes the source term of the turbulence model in our case. The term 
$\boldsymbol {w}_{f}( \dot {\boldsymbol {x}}_{f}\boldsymbol {{\cdot }} \boldsymbol {n})$, which accounts for flux balances in a moving domain as a result of the fluid–structure interaction, with 
$\dot {\boldsymbol {x}}_{f}$ as the Cartesian mesh velocity vector, has been exposed explicitly in (2.2). The RANS equations (plus turbulence model) provide the fluid part of (2.1) and 
$\mathcal {R}_f$ contains the spatial discretisation of the last two integrals and an additional term arising from the geometric conservation law (Thomas & Lombard Reference Thomas and Lombard1979), when dealing with the time-dependent discrete control volumes in the temporal derivative. The dimension of the fluid system is given by the number of mesh points times the number of conservative variables.
The structural system, defined in a suitable domain 
$\varOmega _s(t)$ and communicating with the fluid domain through the fluid/structure interface 
$\varGamma (t)$ (e.g. the wetted surface of the wing structure) (Pfister et al. Reference Pfister, Marquet and Carini2019; Negi et al. Reference Negi, Hanifi and Henningson2020; Lesoinne et al. Reference Lesoinne, Sarkis, Hetmaniuk and Farhat2001), is governed by the second-order ordinary differential equation
\begin{equation}
\boldsymbol{\mathsf{M}}\ddot{\boldsymbol{x}}_s + \boldsymbol{\mathsf{C}}\dot{\boldsymbol{x}}_s +
\boldsymbol{\mathsf{K}}\boldsymbol{x}_s = \boldsymbol{f},
\end{equation}
where matrices 
$\boldsymbol{\mathsf{M}}$, 
$\boldsymbol{\mathsf{C}}$ and 
$\boldsymbol{\mathsf{K}}$ represent the constant mass, damping and stiffness matrices, respectively. We discard structural damping throughout, only accounting for the aerodynamic damping. Vectors 
$\boldsymbol {f}$ and 
$\boldsymbol {x}_s$ are any present forces, specifically the aerodynamic pressure and friction forces acting on the interface 
$\varGamma$ in our case, and the structural coordinates, respectively. The deformation of the structure, 
$\boldsymbol {x}_s$, can be expressed as 
$\boldsymbol {x}_s=\boldsymbol{\varPhi} \boldsymbol {\eta }$, where 
$\boldsymbol{\varPhi}$ is a matrix consisting of the spatial orthogonal mode shapes of the structural system and 
$\boldsymbol {\eta }$ (modal coordinates) provides the time-dependent amplitude of each mode shape contributing to the deformation. The velocity of the structural coordinates is defined as 
${\boldsymbol {u}}_s\equiv \dot {\boldsymbol {x}}_s=\boldsymbol{\varPhi} \dot {\boldsymbol {\eta }}$ for ease of notation in the following. At the interface, 
$\varGamma$, the velocity of the structural points equals those of the fluid domain to enforce the no-slip and no-penetration condition 
$\boldsymbol {u}_s=\dot {\boldsymbol {x}}_{f}$ 
$(=\boldsymbol {u}_{f})$. Matrix 
$\boldsymbol{\varPhi}$, together with the wind-off structural frequencies, is obtained from an undamped free-vibration analysis of the underlying finite-element model. The structural equation (2.3) of second order is projected into the modal space (using a suitable set of 
$m$ dominant modes with lowest frequency) and rewritten as a first-order equation by standard linearisation (Tisseur & Meerbergen Reference Tisseur and Meerbergen2001), specifically 
$\dot {\boldsymbol {w}}_s = \mathcal {R}_s$, using 
$\boldsymbol {w}_s = [\boldsymbol {\eta }^T,\dot {\boldsymbol {\eta }}^T]^{{\rm T}}$ and defining the structural residual 
$\mathcal {R}_s$ as
\begin{equation} \mathcal{R}_s =
\boldsymbol{\mathsf{D}}\boldsymbol{w}_s + \vartheta \boldsymbol{\mathsf{E}} \varPhi^{T}
\boldsymbol{f}(\boldsymbol{w}_{f}, \boldsymbol{w}_s),
\end{equation}
where 
$\boldsymbol{\mathsf{D}} = [0, \boldsymbol{\mathsf{I}};-\varPhi ^{T} \boldsymbol{\mathsf{K}} \varPhi, 0]$ and 
$\boldsymbol{\mathsf{E}} = [0,\boldsymbol{\mathsf{I}}]^{\rm T}$, with 
$\boldsymbol{\mathsf{I}}$ being the 
$m\times m$ identity matrix. The mass ratio 
$\vartheta$ results from writing the equation in dimensionless form, consistent with the fluid equations, and is a function of reference fluid density and length in our formulation. If the mode shapes are mass normalised such that the generalised mass matrix becomes 
$\boldsymbol{\varPhi} ^{T} \boldsymbol{\mathsf{M}} \boldsymbol{\varPhi} = \boldsymbol{\mathsf{I}}$, the generalised stiffness matrix, 
$\boldsymbol{\varPhi} ^{T} \boldsymbol{\mathsf{K}} \boldsymbol{\varPhi}$, is a diagonal matrix containing the squares of the angular frequencies of each structural mode. Hence, the modal structural equations are uncoupled when no generalised aerodynamic forces, 
$\boldsymbol{\varPhi}^T\boldsymbol {f}$, are applied. Aerodynamic coupling leads to aeroelastic phenomena.
A time-invariant static aeroelastic solution of (2.1), 
$\bar {\boldsymbol {w}}=[\bar {\boldsymbol {w}}_{f}^T,\bar {\boldsymbol {\eta }}^T,0]^{{\rm T}}$, satisfying 
$\mathcal {R}(\bar {\boldsymbol {w}})=0$ and referred to as a coupled base state, is of special interest for global stability analysis and shall be detailed explicitly. The elastic restoring force, resulting from the static deformation, is in equilibrium with the generalised aerodynamic forces, and the modal structural equation simplifies to 
$\boldsymbol{\varPhi} ^T \boldsymbol{\mathsf{K}} \boldsymbol{\varPhi} \bar {\boldsymbol {\eta }} = \boldsymbol{\varPhi} ^T \boldsymbol {f}(\bar {\boldsymbol {w}}_{f}, \bar {\boldsymbol {\eta }})$ with the static deformation in physical coordinates in domain 
$\bar {\varOmega }_{s}$ given as 
$\bar {\boldsymbol {x}}_s=\boldsymbol{\varPhi} \bar {\boldsymbol {\eta }}$. Accordingly, the static deformation will define the equilibrium fluid domain 
$\bar {\varOmega }_{f}$ within which the inviscid and viscous fluxes (and source term of turbulence model) in (2.2) are evaluated. Velocities on the interface 
$\bar {\varGamma }$ are, of course, zero. Linear dynamic investigations effectively look at the perturbations around this coupled equilibrium state.
There are a few points to add to complement the preceding description. First, using RANS-level aerodynamics, the dimension of the fluid problem is much larger than that of the structural problem, 
$2m$, due to both the requirements on the resolution of the nonlinear fluid flow and the assumed linear nature of the structural model with constant (e.g. deformation-independent) mass and stiffness matrices. We note, however, that a nonlinear structure will also result in a large number of structural degrees of freedom according to the detail included in the finite-element model. Second, the innocuous appearing unsteady aerodynamic forces, 
$\boldsymbol {f}$, in the context of aircraft aeroelasticity and loads typically represent motion-induced (resulting from the wing motion itself linked to e.g. flutter instability) and gust-induced (with the external excitation originating away from the wing) contributions. As noted by Hall (Reference Hall and Dowell2022), a third unsteady aerodynamic force contribution results from the unstable shock-buffet flow field itself (with the excitation originating at the wing but not requiring its motion). This is a non-trivial discussion. With no attempt to be exhaustive, one aspect here is that the unsteady aerodynamic force and the excited wing motion will mutually interact and modulate one another, with e.g. amplitude-dependent (nonlinear) synchronisation taking place in certain scenarios. While a strict distinction between the motion- and buffet-induced contributions becomes questionable, particularly in an established shock-buffet/buffeting state, buffet loads are often initially derived from an assumed rigid wing. Herein, we will look at this challenge by extracting dominant global modes from a converged time-invariant solution at globally stable pre-onset (and mildly unstable) conditions. Finally, the interface between the fluid and structural systems plays an important role. Normally the structural coordinates do not conform with the fluid surface mesh coordinates, and hence a suitable mapping at the interface between the disciplines is needed, typically by splining the mode shapes to the surface mesh for the coupled analysis to allow the transfer of aerodynamic loads and structural deformations between the fluid and modal structural governing equations without the need for additional transformations.
2.2. Global stability analysis
With the governing equations of the fluid and structural systems defined, we can return to (2.1). Decomposing the unknown fluid and structural variables into static/steady base state and perturbation vectors via 
$\boldsymbol {w}(t) = \bar {\boldsymbol {w}} +\varepsilon \tilde {\boldsymbol {w}}(t)$ (with factor 
$\varepsilon \ll 1$) results, after some more manipulation, in the linearised equations of the form
\begin{equation} \dot{\tilde{\boldsymbol{w}}} = \boldsymbol{\mathsf{J}} \tilde{\boldsymbol{w}}, \end{equation}
where 
$\boldsymbol{\mathsf{J}} = \partial \mathcal {R}/\partial \boldsymbol {w}$ is the coupled Jacobian matrix, which is a large sparse matrix conveniently expressed as consisting of four blocks 
$\boldsymbol{\mathsf{J}}=[\boldsymbol{\mathsf{J}}_{ff}, \boldsymbol{\mathsf{J}}_{fs}; \boldsymbol{\mathsf{J}}_{sf}, \boldsymbol{\mathsf{J}}_{ss}]$, with 
$\boldsymbol{\mathsf{J}}_{sf} = \vartheta \boldsymbol{\mathsf{E}} \varPhi ^{T} \partial \boldsymbol {f}/\partial \boldsymbol {w}_f$ and 
$\boldsymbol{\mathsf{J}}_{ss} = \boldsymbol{\mathsf{D}} + \vartheta \boldsymbol{\mathsf{E}} \boldsymbol{\varPhi} ^{T} \partial \boldsymbol {f}/\partial \boldsymbol {w}_s$. In addition, the fluid-to-structure coupling matrix can be expressed as 
$\boldsymbol{\mathsf{J}}_{fs}=[\boldsymbol{\mathsf{J}}_{f\eta },\boldsymbol{\mathsf{J}}_{f\dot {\eta }}]$, where 
$\boldsymbol{\mathsf{J}}_{f\eta }=\partial \mathcal {R}_{f}/\partial {\boldsymbol {\eta }}$ and similarly 
$\boldsymbol{\mathsf{J}}_{f\dot {\eta }}=\partial \mathcal {R}_{f}/\partial {\dot {\boldsymbol {\eta }}}$. With these definitions, the fluid Jacobian matrix, 
$\boldsymbol{\mathsf{J}}_{ff}=\partial \mathcal {R}_{f}/\partial {\boldsymbol {w}_{f}}$, is obvious. When the standard exponential ansatz of the form 
$\tilde {\boldsymbol {w}}(t) = \hat {\boldsymbol {w}} {\rm e}^{\lambda t}$ is taken, (2.5) is recast as an eigenvalue problem,
\begin{equation} \lambda \hat{\boldsymbol{w}} = \boldsymbol{\mathsf{J}} \hat{\boldsymbol{w}}, \end{equation}
where 
$\lambda = \sigma + {\rm i}\omega$ denotes the eigenvalue, which is a complex number consisting of a decay/growth rate 
$\sigma$ and angular frequency 
$\omega$ (with 
$i$ as the complex unit 
$\sqrt {-1}$), and 
$\hat {\boldsymbol {w}}$ is the corresponding eigenvector containing the complex spatial amplitudes describing the coherent dynamics. When 
$\sigma$ is greater than zero, the system is said to be globally unstable and this is what must be prevented within the aircraft flight envelope. Hence, the solution of this eigenvalue problem is desired as a first step in the analysis (once a nonlinear equilibrium solution 
$\bar {\boldsymbol {w}}$ of (2.1) is available).
A complete study of global instability mechanisms requires the solution of the adjoint equations. These have now become ubiquitous in fluid mechanics, for application in error estimation, shape optimisation, flow control and many more (Luchini & Bottaro Reference Luchini and Bottaro2014). Hence, (2.6) is also solved in its (discrete) adjoint form such that
\begin{equation} \lambda^* \check{\boldsymbol{w}} = \boldsymbol{\mathsf{J}}^{\dagger} \check{\boldsymbol{w}}, \end{equation}
where the eigenvalue will be complex conjugate to that of the direct global problem with corresponding adjoint eigenvector 
$\check {\boldsymbol {w}}$. The adjoint Jacobian operator 
$\boldsymbol{\mathsf{J}}^{\dagger}$ is defined through a suitable inner product 
$\langle \boldsymbol {b}, \boldsymbol{\mathsf{J}}\boldsymbol {a}\rangle = \langle \boldsymbol{\mathsf{J}}^{{\dagger} }\boldsymbol {b}, \boldsymbol {a}\rangle$, where 
$\langle {\boldsymbol {a}},{\boldsymbol {b}}\rangle = {\boldsymbol {a}}^H \boldsymbol{\mathsf{Q}} {\boldsymbol {b}}$ describes the weighted inner product of two arbitrary vectors 
$\boldsymbol {a}$ and 
$\boldsymbol {b}$, with 
$\boldsymbol{\mathsf{Q}}$ expressing a positive–definite matrix. The adjoint operator can thus be given explicitly as 
$\boldsymbol{\mathsf{J}}^{\dagger} =\boldsymbol{\mathsf{Q}}^{-1}\boldsymbol{\mathsf{J}}^T\boldsymbol{\mathsf{Q}}$. We combine a norm that was chosen e.g. in rigid-wing shock-buffet investigations by Sartor, Mettot & Sipp (Reference Sartor, Mettot and Sipp2015) and Paladini et al. (Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019a) and in the study of a rotationally flexible cylinder with splitter plate in Basso et al. (Reference Basso, Hwang, Assi and Sherwin2021). Specifically, the choice of an inner product and its resulting norm physically represent a measure of energy in the coupled system
\begin{equation} \|\boldsymbol{w}\|^2=\langle{\boldsymbol{w}},{\boldsymbol{w}}\rangle = \int_{\varOmega_{f}} \boldsymbol{w}_f^H \boldsymbol{w}_f \,\text{d}V + \boldsymbol{\eta}^H \left(\boldsymbol{\varPhi}^T \boldsymbol{\mathsf{K}} \boldsymbol{\varPhi}\right) \boldsymbol{\eta} + \dot{\boldsymbol{\eta}}^H \left(\boldsymbol{\varPhi}^T \boldsymbol{\mathsf{M}} \boldsymbol{\varPhi} \right)\dot{\boldsymbol{\eta}}. \end{equation}
Hence, the weight matrix 
$\boldsymbol{\mathsf{Q}}$ contains the discrete cell volumes on its diagonal for fluid degrees of freedom as well as generalised stiffness and mass matrices for the structural degrees of freedom. Observe the equivalence between the modal and physical coordinates, specifically for the structural elastic strain energy, 
$\boldsymbol {\eta }^H (\boldsymbol{\varPhi} ^T \boldsymbol{\mathsf{K}} \boldsymbol{\varPhi}) \boldsymbol {\eta } = \boldsymbol {x}_s^H \boldsymbol{\mathsf{K}} \boldsymbol {x}_s$, and similarly for the structural kinetic energy, 
$\dot {\boldsymbol {\eta }}^H \dot {\boldsymbol {\eta }}=\dot {\boldsymbol {\eta }}^H (\boldsymbol{\varPhi} ^T \boldsymbol{\mathsf{M}} \boldsymbol{\varPhi} ) \dot {\boldsymbol {\eta }} = \boldsymbol {u}_s^H \boldsymbol{\mathsf{M}} \boldsymbol {u}_s$, with 
$\boldsymbol{\varPhi} ^T=\boldsymbol{\varPhi} ^H$ and 
$\boldsymbol{\varPhi}$ being mass normalised. Other inner products to define the energy norm, particularly for a compressible fluid, have been discussed in the literature and can be explored in the future (Bugeat et al. Reference Bugeat, Chassaing, Robinet and Sagaut2019; Yeh & Taira Reference Yeh and Taira2019). Due to the non-normality of the Jacobian operator, the set of direct and adjoint global modes form a bi-orthonormal basis where 
$\langle \check {\boldsymbol {w}}_i,\hat {\boldsymbol {w}}_j\rangle =\delta _{ij}$ (once appropriately normalised) with the standard definition of the Kronecker delta 
$\delta _{ij}$. Where applicable, this test of bi-orthonormality was done for all computed modes. Regions of high spatial amplitudes in the adjoint vector can be physically interpreted as those where a harmonic forcing (or initial condition) affects the direct eigensolution most, i.e. the global flow field is most receptive to such imposed perturbations. Related to this is resolvent analysis to identify the maximum response to harmonic forcing, optimised over all possible forcings (Houtman, Timme & Sharma Reference Houtman, Timme and Sharma2022, Reference Houtman, Timme and Sharma2023).
Another extension of the concepts comes from combining the direct and adjoint global modes to form what is often called the wavemaker, which shows the sensitivity of the dynamics to localised feedback and can reveal the origin (or core) of an instability (Giannetti & Luchini Reference Giannetti and Luchini2007). To be consistent and comparable to previous related work (Paladini et al. Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019a,Reference Paladini, Marquet, Sipp, Robinet and Dandoisb), the fluid wavemaker can be defined as the normalised pointwise product of the direct and adjoint eigenvectors, such that
\begin{equation} \theta_{f,i} = {\|\boldsymbol{\mathsf{R}}_i\check{\boldsymbol{w}}_{f}\|_2\,\|\boldsymbol{\mathsf{R}}_i\hat{\boldsymbol{w}}_{f}\|_2}, \end{equation}
where 
$\boldsymbol{\mathsf{R}}_i$ is a diagonal matrix to extract the relevant variables of the direct and adjoint solution vectors at a given location 
$i$ describing a discrete control volume. The vector 2-norm is indicated through 
$\|{\cdot }\|_2$, and, importantly, the direct and adjoint eigenvectors of the coupled system observe the earlier-stated bi-orthonormality to normalise the wavemaker. Similarly, following the work by Basso et al. (Reference Basso, Hwang, Assi and Sherwin2021), the wavemaker of the modal structural system can be stated as (defining 
$\boldsymbol {\chi }\equiv \dot {\boldsymbol {\eta }}$ for ease of notation)
\begin{equation} \theta_{s,j} = {|\boldsymbol{\mathsf{R}}_j\check{\boldsymbol{\eta}}|\,|\boldsymbol{\mathsf{R}}_j\hat{\boldsymbol{\eta}}|}+ {|\boldsymbol{\mathsf{R}}_j\check{\boldsymbol{\chi}}|\,|\boldsymbol{\mathsf{R}}_j\hat{\boldsymbol{\chi}}|} , \end{equation}
where diagonal matrix 
$\boldsymbol{\mathsf{R}}_j$ extracts the modal degrees of freedom denoted by subscript 
$j$. The interested reader is also referred to the work by Skene et al. (Reference Skene, Yeh, Schmid and Taira2022) that includes an insightful discussion on different definitions of the wavemaker in the literature.
2.3. Numerical approach
2.3.1. Static aeroelastic base state calculation
The RANS equations (plus turbulence model) are solved herein using the industrial DLR-TAU code which relies on a second-order, finite-volume, vertex-centred spatial discretisation (Schwamborn, Gerhold & Heinrich Reference Schwamborn, Gerhold and Heinrich2006). The inviscid fluxes are computed using a central scheme with matrix artificial dissipation. The Green–Gauss theorem is used to evaluate the gradients of flow variables required for viscous fluxes and source terms, where needed. Turbulence closure is provided by the negative Spalart–Allmaras model using the Boussinesq eddy-viscosity assumption. The far-field boundary is described as free-stream flow through a characteristic boundary condition consistent with the discretisation of the interior fluxes, while the no-slip, no-penetration condition on viscous walls is enforced strongly. Additionally, when a symmetry plane boundary condition is required (specifically, for our isolated wing test case to verify the implementation, see Appendix A), this is imposed by removing components of the momentum equations normal to the plane. The steady-state flow solution is converged to low residual levels approaching machine epsilon using the implicit backward Euler method with lower–upper symmetric Gauss–Seidel iterations and local time stepping for convergence acceleration. A geometric multi-grid method is also used to improve convergence rates.
The coupled aeroelastic problem requires solving the mass–spring modal structural equations with applied (generalised) aerodynamic forces iteratively in a staggered fashion updating the structural and fluid degrees of freedom in turn. The mode shapes were mapped one to one to the surface mesh for the fluid equations through interpolation using radial basis functions (see e.g. Michler Reference Michler2011). This allows the transfer of aerodynamic forces and structural deformations between the fluid and structural equations without the need for additional transformations. While the Newmark-beta scheme can be used to integrate the second-order structural equations in time, for the static aeroelastic coupling where the aerodynamic loads are balanced by the wing stiffness, specifically 
$\boldsymbol{\varPhi} ^{T} \boldsymbol{\mathsf{K}} \boldsymbol{\varPhi} \,\bar {\boldsymbol {\eta }} = \boldsymbol{\varPhi} ^{T}\boldsymbol {f}(\bar {\boldsymbol {w}}_{f},\bar {\boldsymbol {\eta }})$, the structural degrees of freedom are updated iteratively based on the latest loads estimate. The coupling between aerodynamics and structure for these nonlinear iterations (either for computing a base state or for unsteady time marching) is done using the tools provided through the FlowSimulator framework, see for instance Reimer, Heinrich & Ritter (Reference Reimer, Heinrich and Ritter2020). However, the fluid–structure coupling for linearised analyses in identifying dominant global modes has been integrated entirely in the TAU flow solver as part of this work with detail provided in the following.
2.3.2. Aeroelastic triglobal stability calculation
The Jacobian matrix blocks of the linearisation are evaluated on the statically deformed geometry for subsequent stability analysis. Matrices 
$\boldsymbol{\mathsf{J}}_{ff}$ and 
$\boldsymbol{\mathsf{J}}_{sf}$ are computed from a hand-derived, analytical formulation, while 
$\boldsymbol{\mathsf{J}}_{fs}$ is computed using a central finite-difference residual evaluation of the fluid equations on the deformed volume mesh, with mesh deformation applied through a radial basis function approach (see e.g. Michler Reference Michler2011). Specifically, a finite-difference step size of 
$10^{-5}$ is used throughout and the sensitivity of the computed global modes was found to be negligible. Also see previous work (e.g. Timme & Thormann Reference Timme and Thormann2016) where the frequency-domain aerodynamic response to structural forcing, applied through the columns of matrix 
$\boldsymbol{\mathsf{J}}_{fs}$, has been assessed with respect to time-marching pulse excitation. Forming 
$\boldsymbol{\mathsf{J}}_{ss}$ is trivial due to the linear modal nature of the structural system. Note that the viscous contribution to the aerodynamic forces 
$\boldsymbol {f}$ is discarded herein for stability analysis, leaving the pressure components only, and the impact of this simplification in the analysis of aircraft wing aeroelasticity was found to be negligible when compared with time-marching simulation that accounted for the full aerodynamic force vector including viscous terms (Belesiotis-Kataras & Timme Reference Belesiotis-Kataras and Timme2021). Details on the meaning of the different matrices and simplifications can be found in Badcock et al. (Reference Badcock, Timme, Marques, Khodaparast, Prandina, Mottershead, Swift, Da Ronch and Woodgate2011), while a thorough discussion of the underlying linear harmonic solver in DLR-TAU can be found in Thormann & Widhalm (Reference Thormann and Widhalm2013).
In previous work (Badcock et al. Reference Badcock, Timme, Marques, Khodaparast, Prandina, Mottershead, Swift, Da Ronch and Woodgate2011; Timme & Badcock Reference Timme and Badcock2011; Timme, Marques & Badcock Reference Timme, Marques and Badcock2011), the coupled fluid–structure system was solved using the Schur complement eigenvalue method. This method utilises the Schur complement matrix 
$\boldsymbol{\mathsf{S}}$ to solve for the structural part of the eigenvalue problem in (2.6). Rearranging gives
\begin{equation} \boldsymbol{\mathsf{S}}(\lambda)\hat{\boldsymbol{w}}_{s} =\lambda \hat{\boldsymbol{w}}_{s}, \end{equation}
where 
$\boldsymbol{\mathsf{S}}(\lambda )$ is expressed as
\begin{equation} \boldsymbol{\mathsf{S}}(\lambda) = \boldsymbol{\mathsf{J}}_{ss} - \boldsymbol{\mathsf{J}}_{sf}(\boldsymbol{\mathsf{J}}_{ff} - \lambda \boldsymbol{\mathsf{I}})^{{-}1} \boldsymbol{\mathsf{J}}_{fs}, \end{equation}
and 
$\lambda$ is an eigenvalue originating in the uncoupled structural system. This is an important condition, as will be seen below. The second part of the right-hand side in (2.12) is dubbed the interaction (or coupling) matrix 
$\boldsymbol{\mathsf{S}}^c=-\boldsymbol{\mathsf{J}}_{sf}(\boldsymbol{\mathsf{J}}_{ff} - \lambda \boldsymbol{\mathsf{I}})^{-1} \boldsymbol{\mathsf{J}}_{fs}$ and can be formed in both the frequency and time domain (Timme & Badcock Reference Timme and Badcock2011; Timme et al. Reference Timme, Marques and Badcock2011). It was noted by the authors that determining the values of the Schur complement matrix, specifically the aerodynamic coupling matrix 
$\boldsymbol{\mathsf{S}}^c$, requires considerable computational effort, if not done efficiently. As it depends on the steady-state solution and the eigenvalue, scanning over a large parameter space quickly becomes prohibitive, and thus approximations are needed. While multiple approaches to approximating matrix 
$\boldsymbol{\mathsf{S}}^c$ were offered, a linear spline interpolation is chosen as a surrogate herein while iterating to the eigensolution. Depending on using either 
$\lambda$ or 
$i\omega$ when forming matrix 
$\boldsymbol{\mathsf{S}}^c$, the Schur complement method can be shown to be equivalent to either a classical p or pk flutter analysis (Hassig Reference Hassig1971). It is also intuitive that the absence of aerodynamics leads to the uncoupled, linear eigenvalue problem for the structural dynamics, 
$\boldsymbol{\mathsf{D}}\hat {\boldsymbol {w}}_{s} =\lambda \hat {\boldsymbol {w}}_{s}$. Finally, the fluid part of the coupled system can be computed from the other equation of the original coupled matrix problem in (2.6):
\begin{equation} (\boldsymbol{\mathsf{J}}_{ff}-\lambda \boldsymbol{\mathsf{I}}) \hat{\boldsymbol{w}}_{f}={-}\boldsymbol{\mathsf{J}}_{fs}\hat{\boldsymbol{w}}_{s}, \end{equation}and these ideas, and equivalent relations for the adjoint eigenvalue problem, have been exercised in previous work (Bekemeyer & Timme Reference Bekemeyer and Timme2019).
There are two requirements for the Schur complement method to work; an appropriate shift based on the structural (wind-off) frequencies (or a previously converged solution) and the matrix decomposition with relevant eigenmodes from the matrix 
$\boldsymbol{\mathsf{J}}_{ss}$. With regards to applying the same method for buffet studies, whilst the former requirement can be obtained through engineering insight even for a fluid mode, it is not possible to obtain the latter when interest is in modes emerging from the fluid system. Specifically, we are not aware of a suitable matrix decomposition that would isolate the leading flow physics linked to the buffet phenomenon to similarly apply the efficient Schur complement solver. Hence, a coupled eigenvalue solver is needed, incorporating the full Jacobian matrix. This can be done by using the methods available through the ARPACK library.
The implicitly restarted Arnoldi method, as implemented in the ARPACK library, is a routine for finding a set of eigenvalues for large sparse matrices (Sorensen Reference Sorensen1992). Starting with a random search direction 
$\boldsymbol {v}$, the algorithm computes, in principle, the so-called Krylov subspace 
$\mathcal {K}_r = [\boldsymbol {v}, \boldsymbol{\mathsf{J}}\boldsymbol {v}, \boldsymbol{\mathsf{J}}^{2}\boldsymbol {v},\ldots, \boldsymbol{\mathsf{J}}^{r-1}\boldsymbol {v}]$ of dimension 
$r$. After each successive multiplication with the Jacobian operator 
$J$, the resulting vector is orthonormalised and added as column to a matrix 
$\boldsymbol{\mathsf{V}}_r$. Following projection, a few Ritz eigenvalues 
$\lambda _{r}$ and corresponding Ritz approximate eigenvectors 
$\hat {\boldsymbol {w}}=\boldsymbol{\mathsf{V}}_r\boldsymbol {y}_{r}$ (where eigenvector 
$\boldsymbol {y}_{r}$ is associated with 
$\lambda _{r}$) of the upper Hessenberg matrix 
$\boldsymbol{\mathsf{H}}_r = \boldsymbol{\mathsf{V}}^{H}_r \boldsymbol{\mathsf{J}} \boldsymbol{\mathsf{V}}_r$ are good approximations of those largest eigenvalues of 
$\boldsymbol{\mathsf{J}}$. A major advantage is that the algorithm needs only to compute the matrix-vector product 
$\boldsymbol{\mathsf{J}}\boldsymbol {v}$, so the full matrix 
$\boldsymbol{\mathsf{J}}$ does not need to be obtained and stored explicitly. The algorithm can also be run in shift-invert mode to find the largest eigenvalues closest to a user-defined complex-valued shift 
$\zeta$ by operating on vector 
$\boldsymbol {v}$ with 
$\boldsymbol{\mathsf{A}}^{-1}=(\boldsymbol{\mathsf{J}}-\zeta \boldsymbol{\mathsf{I}})^{-1}$ instead of 
$\boldsymbol{\mathsf{J}}$. To achieve this, a linear system of the size of the coupled problem needs to be solved and a practical way is using an iterative Krylov subspace solver, which again heavily relies on matrix-vector products now with the shifted Jacobian matrix 
$\boldsymbol{\mathsf{A}}$ and this can, in principle, be done matrix free for low memory footprints (Knoll & Keyes Reference Knoll and Keyes2004; Vevek, Houtman & Timme Reference Vevek, Houtman and Timme2022a). However, the matrix-forming approach is used herein. An inner–outer iterative eigenvalue method is thus established, whereby the outer implicitly restarted Arnoldi process is used to find the eigenvalues while an inner iterative solver enables the shift-invert operation. The preconditioned generalised conjugate residual method with inner orthogonalisation and deflated restarting (GCRO-DR) is chosen as the inner iterative solver (Xu, Timme & Badcock Reference Xu, Timme and Badcock2016).
2.3.3. Preconditioned sparse iterative solver
The GCRO-DR solver is an iterative method for solving systems of linear equations (Eiermann, Ernst & Schneider Reference Eiermann, Ernst and Schneider2000; Parks et al. Reference Parks, de Sturler, Mackey, Johnson and Maiti2006). It utilises the restarted Arnoldi method for solving the system 
$\boldsymbol{\mathsf{A}}\boldsymbol {x} = \boldsymbol {b}$ by seeking a solution vector 
$\boldsymbol {x}$ that minimises the norm of the residual, 
$\|\boldsymbol {b} - \boldsymbol{\mathsf{A}}\boldsymbol {x}\|$. Specifically, compared with the standard generalised minimal residual method, GCRO-DR recycles the Krylov subspace between restarts, with the potential to avoid convergence stall and indeed to accelerate convergence while reducing the required size of the subspace. Whichever Krylov subspace method is used, a preconditioned version is normally needed in practice to achieve convergence, with the preconditioner 
$\boldsymbol{\mathsf{P}}^{-1}$ being an approximation of 
$\boldsymbol{\mathsf{A}}^{-1}$. Ideally, the operator 
$\boldsymbol{\mathsf{P}}^{-1}$ should be as close to 
$\boldsymbol{\mathsf{A}}^{-1}$ as possible without incurring a large computational cost for both computing and applying the preconditioner. In doing so, using preconditioning reduces the condition number of the system and can offer significant speed up.
As in previous work, when running the linear harmonic solver in the chosen flow solver in a distributed-memory parallel computing environment, the implementation approximates the fluid-only matrix, 
$\boldsymbol{\mathsf{J}}_{ff}$, used for preconditioning by a block-Jacobi form. Specifically, when running 
$n$ parallel processes, the fluid Jacobian matrix (according to the partitioned mesh) is split into 
$n$ block rows (one block row per process including all entries needed for global matrix-vector products), and only the block-diagonal part of it, local to this process, is used for preconditioning. Note that it was found beneficial for reasons of rate of convergence and stability to factorise a matrix which combines Jacobian matrices arising from approximate first-order and exact second-order spatial schemes (McCracken et al. Reference McCracken, Da Ronch, Timme and Badcock2013; Vevek et al. Reference Vevek, Timme, Pattinson, Stickan and Büchner2022b). This matrix is then factorised by the incomplete lower–upper method with no fill in, denoted ILU(0). However, such approach poses a challenge for a fluid–structure coupled system in parallel, as the coupling matrices, 
$\boldsymbol{\mathsf{J}}_{fs}$ and 
$\boldsymbol{\mathsf{J}}_{sf}$, would be discarded, when using an equivalent discipline-level block-Jacobi formulation extended to the structural degrees of freedom. A new preconditioning approach, henceforth referred to as the arrowhead preconditioner, was developed in the parallel implementation to address this coupling challenge. Initial testing was done with a simple case, specifically the Goland wing, as the corresponding coupled Jacobian matrix was small enough to fit into memory of a single core and be factorised using ILU(0), without discarding the coupling blocks. The results from these sequential tests were then used as a benchmark to verify the implementation of the various code additions in parallel and assess the performance. Details of the methods are described in the appendices.
3. The NASA common research model
The NASA Common Research Model (CRM) resembles a modern passenger aeroplane and exists as both a physical model (for wind-tunnel testing) and a computational model. It was designed as a universal test case for researchers to compare new ideas and results (Vassberg et al. Reference Vassberg, Dehaan, Rivers and Wahls2008). The wing has an aspect ratio of 9, a taper ratio of 0.275 and a 
$35^\circ$ quarter-chord sweep angle. Herein, the scaled-down wind-tunnel wing/body/horizontal-tail version is discussed featuring a mean aerodynamic chord of 0.189 m with a full span of 1.586 m and reference area of 
$0.280\ {\rm m}^{2}$. The pylons and nacelles were discarded and the horizontal tail-setting angle was 
$0^\circ$. All design details including aerofoil data can be found in the cited reference.
An unstructured computational mesh has previously been generated for the half-span simulations in Timme (Reference Timme2020) with the SOLAR mesh generator (Martineau et al. Reference Martineau, Stokes, Munday, Jackson, Gribben and Verhoeven2006) and, upon mirroring with respect to the symmetry centre plane, the full-span case has approximately 
$12.3\times 10^6$ mesh points with 
$3.3\times 10^5$ points on solid walls. A viscous wall normal spacing of 
$y^+<1$ is ensured, eliminating the need for wall functions in the flow model. The spherical far-field boundary is located 100 semi-span lengths away.
The Reynolds number based on mean aerodynamic chord is 
$Re = 5.0\times 10^{6}$ and the reference free-stream Mach number is 
$M = 0.85$, chosen based on runs 153/182 of the European Transonic Windtunnel (ETW) test campaign (Lutz et al. Reference Lutz, Gansel, Waldmann, Zimmermann and Schulte am Hülse2016). In accordance with our previous work (Timme Reference Timme2020), and contrary to the experimental configuration where the transition location is fixed at 10 % local chord length, herein, a fully turbulent boundary layer is assumed. Run 182 measured the static deformation of the elastic wing at several angles of attack. For intermediate angles not measured, but required for our study around the shock-buffet onset condition, interpolation of the experimentally measured deformations was used (Keye & Gammon Reference Keye and Gammon2018). In contrast to the same previous work, herein the impact of both static and dynamic aeroelastic deformation is simulated (with validation results provided in the following) rather than pre-deforming the wing based on the given experimental data and running the fluid simulations on the then quasi-rigid wing. As such, additional information on the wind-tunnel configuration and conditions is required. The density ratio (that is, solid-to-fluid density) is approximately 
$5000$ and the loading factor (i.e. dynamic pressure over Young's modulus) is 
$0.33\times 10^{-6}$ for the case considered herein. The total temperature and total pressure during the runs were approximately 
$300$ K and 
$192$ kPa, respectively.
The finite-element model used for this study represents the NASA CRM wind-tunnel geometry (excluding the tunnel itself) and is publicly available. Following a modal structural analysis, the first 30 normal modes with lowest frequency are kept (Tinoco et al. Reference Tinoco2018) which is deemed sufficient considering that the model comes without pylon/nacelle, vertical tail plane or control surfaces. Importantly, the structural frequency range covers the shock-buffet frequency range previously identified in Timme (Reference Timme2020). Figure 1 shows some of these modes, conveniently scaled for visualisation purposes. As expected for such an aircraft model, the first two modes describe wing bending, as seen for the starboard wing in figure 1(a). Higher-frequency modes describe various combinations of wing twist and bending and more complex variations. Other modes primarily capture fuselage (modes 3 to 6) and tail (modes 9 to 11) motions, and these are then not expected to be dominant in this study.
Representative structural mode shapes for NASA CRM test case with corresponding wind-off frequencies. Surface colours indicate the modal deformation in 
${z}$-direction. Peculiar features on the aircraft surface are the various cut outs on the wind-tunnel model for experimental sensors and to house the instrumentation. (a) Mode 2, 40.94 Hz; (b) mode 20, 426.74 Hz; (c) mode 25, 526.93 Hz; (d) mode 26, 568.48 Hz; (e) mode 28, 641.54 Hz; ( f) mode 30, 679.62 Hz.
4. Results and discussion
We start with validating the fluid–structure coupled simulations with respect to available experimental data. This is followed by a flutter analysis and a recap of the aerodynamic stability analysis, although assessing the impact of model asymmetry herein not previously considered in Timme (Reference Timme2020). Last but not least, the majority of the discussion addresses the aeroelastic stability analysis. All results are stated in their non-dimensional form, based on the mean aerodynamic chord and reference free-stream values, unless explicitly specified otherwise. Following the insight gained in previous work, typical parameter settings for the present eigenvalue computations are summarised in table 1. The computation of the interaction matrix for the flutter analysis, which relies on the same linear harmonic solver, uses parameter settings in accordance with the inner iterations of the eigenvalue solver. Full-span calculations resulting in nearly 
$74\times 10^6$ complex-valued degrees of freedom are done on four compute nodes, each having twin Skylake 6138 processors, 40 hardware cores and 384 GB of memory.
Typical parameter settings used for inner–outer eigenvalue solver.
4.1. Validation of simulation set-up
The static aeroelastic solution of the aircraft model was computed at angles of attack 
$\alpha = 3.0^\circ$, 
$3.5^\circ$, 
$3.7^\circ$ and 
$3.75^\circ$ to match the loaded wind-tunnel shape (interpolated from run 182 in the ETW campaign) and surface pressure data. Figure 2 shows the wing bending and twist deformations at 
$\alpha = 3.75^\circ$ on the port and starboard wings (evaluated at 50 % chord) compared with those measured using stereo pattern tracking via markers distributed over the wing surface. These markers were affixed to the lower surface of the port wing (Lutz et al. Reference Lutz, Gansel, Waldmann, Zimmermann and Schulte am Hülse2016). Overall good agreement can be observed, on par with other fluid–structure coupled simulations (Tinoco et al. Reference Tinoco2018), with a maximum bending of approximately 18 mm and a washout twist of 
$-1.2^\circ$ at the wing tip. Note the variability in the experimental data themselves, specifically for the twist deformation in figure 2(b), previously reported in Keye & Rudnik (Reference Keye and Rudnik2015). The figure includes error bars describing the confidence interval encompassing 90 % of all recorded values for a single data point. The relevant deformation data for our study were available only for angles of attack 
$\alpha =3.0^\circ$ and 
$4.0^\circ$ during the pitch-pause polar. Considering that the difference in the fluctuation range for these two angles of attack was quite small, linear interpolation was used for the intermediate 
$\alpha =3.75^\circ$ (similar to obtaining the mean deformation itself). Importantly, an asymmetry between the port and starboard wing is also clearly visible in the numerical data. The cause of this asymmetry lies in the high fidelity of the finite-element model which takes into account the different cut outs on each wing to house experimental sensors and related instrumentation as well as the asymmetric inner structure of the model itself. These features propagate to the normal mode shapes and their frequencies. For instance, the first two modes are wing bending for port (at wind-off structural frequency 39.4 Hz) and starboard (40.9 Hz), respectively, with the latter seen in figure 1(a). The corresponding surface pressure coefficient of the steady base flow can be found in figure 3 with good agreement to the experimental measurements at most spanwise stations. Herein, the interest is below and around onset conditions of a global instability which means, in our modelling framework, the time-invariant, static aeroelastic solution is approximately equal to the time-averaged mean state, which is what is shown for the experimental data. A discussion on the seemingly missing experimental data points around the mid semi-span was given in Tinoco et al. (Reference Tinoco2018) and Timme (Reference Timme2020). Differences in the surface pressure between port and starboard wing are minor and not noticeable to the naked eye. Similar levels of agreement were found in the comparison for all other angles of attack used in this study. Overall, the results are reassuring that the fluid–structure coupling has been done correctly (such as defining the mass ratio 
$\vartheta$ and various scalings) and the models are of sufficient fidelity.
Wing deformation showing (
$a$) bending and (
$b$) twist over dimensionless span 
${\eta }$ at angle of attack 
${\alpha = 3.75^\circ }$ comparing experimental data from the ETW campaign and static aeroelastic results emphasising differences on port and starboard wings.
Comparison of experimental and numerical surface pressure coefficient data at angle of attack 
${\alpha = 3.75^\circ }$ and eight spanwise locations. Streamwise coordinate 
${x}$ is normalised by corresponding local chord length 
${c}$.
4.2. Flutter analysis
Initially, eleven samples of the interaction matrix 
$S^c$ were computed using the linearised frequency-domain solver, assuming simple harmonic structural motion at reduced frequencies in the range 
$\omega =0$ to 3 in increments of 0.3. The sampling frequency range corresponds to the wind-off structural frequencies, once made dimensionless based on the actual wind-tunnel flow conditions. Specifically, computing a matrix 
$\boldsymbol{\mathsf{S}}^c(\omega )$ involves first solving linear systems with the frequency-shifted fluid Jacobian matrix, 
$(\boldsymbol{\mathsf{J}}_{ff}-{\rm i}\omega \boldsymbol{\mathsf{I}})$, for each of the 
$m$ columns of 
$(\boldsymbol{\mathsf{J}}_{f\eta }+{\rm i}\omega \boldsymbol{\mathsf{J}}_{f\dot {\eta }})$ and, second, the projection with the matrix 
$\boldsymbol{\mathsf{J}}_{sf}$. Combining the 
$2m$ columns of the matrix 
$\boldsymbol{\mathsf{J}}_{fs}$ into 
$m$ columns has been discussed at length in other work and is not repeated herein (Badcock et al. Reference Badcock, Timme, Marques, Khodaparast, Prandina, Mottershead, Swift, Da Ronch and Woodgate2011; Timme, Badcock & Da Ronch Reference Timme, Badcock and Da Ronch2013). This was done at three angles of attack 
$\alpha =3.0^\circ$, 
$3.5^\circ$ and 
$3.75^\circ$. Upon inspection of the results at the highest angle of attack, additional samples were added in regions of significant activity, bringing the total to 23. Figure 4 shows both the absolute value and phase of some selected entries of matrix 
$\boldsymbol{\mathsf{S}}^c$, specifically of its submatrix 
$\boldsymbol{\mathsf{G}}=(\boldsymbol{\varPhi} ^{T} \partial \boldsymbol {f}/\partial \boldsymbol {w}_f)(\boldsymbol{\mathsf{J}}_{ff}-{\rm i}\omega \boldsymbol{\mathsf{I}})^{-1}(\boldsymbol{\mathsf{J}}_{f\eta }+{\rm i}\omega \boldsymbol{\mathsf{J}}_{f\dot {\eta }})$, at two angles of attack. Values at angle of attack 
$\alpha =3.0^\circ$ were found to be almost identical to those at 
$\alpha =3.5^\circ$ and are not shown here for clarity. At the lower angle of attack shown in the figure, the matrix elements describe a smooth trend with respect to the frequency. At the higher angle of attack, which corresponds to a shock-buffet condition in the fluid-only analysis, the entries have significant variation in absolute value and phase in the frequency range where the band of aerodynamic modes with increased decay rate is observed (Timme Reference Timme2020), concretely between approximately 
$\omega =2$ and 
$3$, indicating, first, a strong aerodynamic response to a structural forcing and, second, a strong coupling between the fluid and structure as discussed shortly. The cause lies both in the proximity of the sampling frequency to some eigenvalues of the fluid system and in the non-normality of the governing equations (He & Timme Reference He and Timme2020). For instance, the significant peaks at a forcing frequency of approximately 
$\omega =2.7$ coincide with one computed eigenvalue, labelled 
$c'$, of the band of shock-buffet modes shown in figures 5 and 7 below. It must be noted though that interpreting a linearised frequency-domain aerodynamic response to a structural forcing in the globally unstable flow at angle of attack 
$\alpha =3.75^\circ$ should be done carefully. Also observe the jumps in the phase going from 
$-{\rm \pi}$ to 
${\rm \pi}$ which is due to visualisation purposes, whereas the actual flutter calculations with the matrix 
$\boldsymbol{\mathsf{S}}^c$ are done with the equivalent Cartesian complex numbers instead of the modulus and phase shown in the figure.
Log magnitude (a–d) and phase (e–h) of complex entries of matrix 
$\boldsymbol{\mathsf{G}}$ at angles of attack 
${\alpha = 3.5^\circ }$ and 
${3.75^\circ }$ over sampling frequency 
${\omega }$. Indices above the plots refer to matrix entries 
${\mathsf{G}}_{ij}$. Symbols indicate sample locations.
Eigenvalues originating in structural system based on pk-type flutter analysis at 
${\alpha = 3.0^\circ }$, 
$3.5^\circ$, 
$3.7^\circ$ and 
${3.75^\circ }$ at the flow condition encountered in the wind-tunnel environment. Mode 28 at 
${\alpha = 3.75^\circ }$, denoted 
${\blacksquare }$, failed to converge and the proper value, computed with the Arnoldi method, is therefore also shown. A close-up view, indicated by the red box in (
$a$), is shown in (
$b$) for clarity. The fluid modes denoted by 
$b$ and 
$c'$ correspond to the leading buffet and destabilised modes, respectively, found with the Arnoldi method and are included to show their proximity to the structural modes.
Once the interaction matrices have been computed and assessed, they feed into the flutter stability calculation, specifically (2.11). Figure 5 shows the predicted eigenvalues from the flutter analysis at angles of attack 
$\alpha = 3.0^\circ$, 
$3.5^\circ$, 
$3.7^\circ$ and 
$3.75^\circ$ at the target flow condition encountered in the wind-tunnel environment. Since tracing the eigenmodes from the uncoupled system state is somewhat arbitrary and would depend on how the test point is reached in the wind tunnel, we chose to increase the density until the simulated flow condition matches the experiment, while keeping all other variables, specifically velocity, temperature and pressure, frozen at the target condition. The target corresponds to a density of approximately 
$1.53\ {\rm kg}\ {\rm m}^{-3}$. It should be emphasised that, while the chosen approach to tracing the modes is arbitrary and does not necessarily agree with how the target condition is reached in the wind tunnel, the mode tracing agrees with the direct calculation using the Arnoldi method, as shown in § 4.4. The figure itself appears rather busy, hence a moment is taken to explain it step by step. The eigenvalues are shown in the complex plane. While the focus is on the eigenmodes originating in the structural system for wind-off conditions, specifically those of the matrix 
$\boldsymbol{\mathsf{J}}_{ss}$, the figure also includes the shock-buffet modes, labelled 
$b$ and 
$c'$ hereafter consistent with previous work (Timme Reference Timme2020), in faint colour for angles of attack 
$\alpha =3.7^\circ$ (mode 
$b$) and 
$3.75^\circ$ (mode 
$c'$), to demonstrate the connection with the results presented in figure 7. Strictly speaking, denoting the aeroelastic modes as either structural or fluid modes would not be correct due to the coupled nature of the problem. For ease of writing, however, we will do so nevertheless, when the origin of the modes in the uncoupled case can be unambiguously traced to either the structure or the fluid, as is the case in our investigation.
In the figure, the eigenvalues at angles of attack 
$\alpha = 3.0^\circ$ and 
$3.5^\circ$ show no instabilities and the modes do not migrate significantly with increasing angle of attack. At angles of attack 
$\alpha = 3.7^\circ$ and 
$3.75^\circ$, it is evident that modes with frequencies below approximately 
$\omega =1.7$ have changed also very little compared with the lower angles of attack. At the same time, instabilities are found for different structural modes at these two higher angles of attack. Modes 29 and 30 are unstable at 
$\alpha = 3.7^\circ$, but decrease greatly in growth/decay rate at 
$\alpha = 3.75^\circ$. This strong decrease is also seen for modes 22, 23 and 25. In contrast, modes 19 to 21 are destabilised further to have an increased growth rate, bringing these modes into the unstable half-plane. However, mode 28 at angle of attack 
$\alpha =3.75^\circ$ (denoted by symbol 
$\blacksquare$), did not trace correctly despite all efforts and the additional samples in this frequency range to avoid heavy reliance on interpolation of the matrix 
$\boldsymbol{\mathsf{S}}^c$. Instead it jumped onto an odd state. This was attributed to the fact that the entries of the interaction matrix 
$\boldsymbol{\mathsf{S}}^c$, and therefore 
$\boldsymbol{\mathsf{S}}(\lambda )$ in (2.12), show extreme variation due to the shifted fluid Jacobian matrix, 
$(\boldsymbol{\mathsf{J}}_{ff}-\lambda \boldsymbol{\mathsf{I}})$, being close to singular due to the proximity to the mode labelled 
$c'$ in figure 7. Indeed, mode 
$c'$ crosses the imaginary axis when going from a fluid-only to the coupled system, passing in close distance to the troublesome mode 28 while increasing the value of fluid density. This provides another reason for the use of the Arnoldi method, besides the ability to compute fluid modes in the first place. It can distinguish between fluid and structural modes in highly contested regions. The correct mode 28 coming from the ARPACK calculation is therefore also included in the figure, denoted by the green square labelled ‘Mode 28’. It has been noted that the second term of matrix 
$\boldsymbol{\mathsf{J}}_{ss} = \boldsymbol{\mathsf{D}} + \vartheta \boldsymbol{\mathsf{E}} \boldsymbol{\varPhi} ^{T} \partial \boldsymbol {f}/\partial \boldsymbol {w}_s$ can often be discarded (Timme et al. Reference Timme, Marques and Badcock2011). This was confirmed herein by evaluating the sensitivity of the generalised forces, 
$\boldsymbol{\varPhi} ^{T}\boldsymbol {f}$, with respect to the modal amplitudes 
$\boldsymbol {\eta }$ (at fixed base-flow solution). The added term had negligible influence on the results. Viscous force contributions were also discarded for similar reasons, as outlined earlier.
To identify the composition of the structural entries of an eigenmode, we adopted and adapted the modal assurance criterion (MAC) (Allemang Reference Allemang2003), defined herein as
\begin{equation} \text{MAC}(i,j) = \frac{|\boldsymbol{e}_i^T\boldsymbol{\hat{\eta}}_j |^2}{(\boldsymbol{e}_i^T\boldsymbol{e}_i)\,\,(\boldsymbol{\hat{\eta}}_j^H\boldsymbol{\hat{\eta}}_j)}, \end{equation}
with 
$\boldsymbol {e}_i$ as linearly independent unit basis vectors to represent one free-vibration mode shape at a time and 
$\boldsymbol {\hat {\eta }}_j$ indicating the structural part of the coupled eigenmodes. In general terms, the MAC gives an indication of the consistency between modes (such as those measured experimentally and computed from a numerical model) and can take a value between 0 and 1, with 1 indicating fully consistent mode shapes. Figure 6 shows the MAC of the structural part, 
$\hat {\boldsymbol {\eta }}$, of the structural eigenmodes computed using the Schur complement method compared with the amplitudes of the wind-off normal modes for angles of attack of 
$\alpha =3.0^\circ$, 
$3.5^\circ$ and 
$3.75^\circ$. (Mode 28 at angle of attack 
$\alpha =3.75^\circ$ uses the solution from the Arnoldi method, as explained earlier.) Recall that the aircraft deformation results from 
$\boldsymbol {x}_s=\boldsymbol{\varPhi} \boldsymbol {\eta }$ with 
$\boldsymbol{\varPhi}$ as the column matrix of spatial orthogonal normal mode shapes. The structural part of the six fluid modes, specifically the pairs 
$a$, 
$b$ and 
$c$, computed using the Arnoldi method are included as well for the highest considered angle of attack as it is rather instructive. For the structural modes at all angles of attack, the MAC is dominated by the diagonal, indicating the largest contributions to the eigenmodes come from the corresponding normal modes. This also gives a means of identifying whether a mode originates in the fluid or structural system, as the MAC for a structural mode is likely to be dominated by a single entry. At angles of attack 
$\alpha = 3.0^\circ$ and 
$3.5^\circ$, no discernible difference is present and all off-diagonal entries are at least an order of magnitude lower, as no significant input from other modes is present. Noteworthy are the three pairs of modes, specifically modes 1/2, 14/15 and 19/20, which have a comparably high MAC between them. This can be explained by these pairs having similar mode shapes and/or close wind-off frequencies. For instance, modes 1 (39.4 Hz) and 2 (40.9 Hz) are both wing bending emphasising either port or starboard deformation. At angle of attack 
$\alpha = 3.75^\circ$, the MAC for modes 1 to 18 shows rather similar features compared with the two lower angles of attack. However, the structural parts of the aeroelastic eigenvectors 19 to 30 now come with a significant contribution from a wide range of normal modes. This observation agrees with the discussion around figure 5. Modes 20 to 25, the frequencies of which lie around the frequency of the leading fluid modes, seem to have a particularly salient interconnection, indicating that the aerodynamics related to the shock-buffet phenomenon can cause a strong coupling between structural degrees of freedom. In addition, a clear contribution gap of normal modes 3 to 6 and 9 to 11 to those active higher-frequency aeroelastic modes is noticeable. Upon inspection of the free-vibration mode shapes, it can be said that modes 3 to 6 are dominated by fuselage bending and modes 9 to 11 are dominant on the horizontal tail. Indeed, the shock-buffet unsteadiness on the wings (and its wake) does not heavily impact on the horizontal tail in this study. Finally, the structural parts of the shock-buffet modes (
$a$ to 
$c$), discussed in more detail below, have contributions from most normal modes (except those aforementioned fuselage and tail modes) and, in particular, from those in the same frequency range as the shock-buffet dynamics.
Visualisation of modal assurance criterion correlating the structural part of the aeroelastic (structural) modes with the amplitudes of the wind-off finite-element method (FEM) modes at different angles of attack. At angle of attack 
$\alpha =3.75^\circ$ the corresponding values of the leading aeroelastic (fluid) modes, labelled 
$a$, 
$b$ and 
$c$, are included (with mode 
$c'$ in final column).
4.3. Aerodynamic global stability analysis
The interest now turns to elucidating the impact of the asymmetrically deformed full-span wing geometry on the aerodynamic stability characteristics at angles of attack 
$\alpha =3.7^\circ$ and 
$3.75^\circ$. For the eigenspectra computed with the shift-invert Arnoldi method, multiple shifts were used to cover the relevant frequency range. Figure 7(a) shows the eigenspectra of the fluid-only system, as computed on the perfectly symmetric case (symmetric with respect to the fuselage centre plane) and the asymmetric case from our static aeroelastic simulation, cf. the deformation presented in figure 2. The full-span symmetric data are taken from Timme (Reference Timme2020) and the labelling follows accordingly. While the symmetric, pre-deformed geometry (corresponding to deformations measured in run 182 of the experimental campaign) gives two nearly identical unstable eigenvalues at approximately 
$\lambda = 0.16 + 2.37i$, the asymmetric geometry from the static coupling simulation results in two visibly distinct eigenvalues. In the symmetric case, the corresponding eigenvectors revealed symmetric/anti-symmetric coherent spatial features of equal amplitudes on both wings. An interpretation is that unsteadiness on the two wings is more or less independent (leaving some coupling effects from the global flow field aside). Interestingly, for the asymmetric, statically deformed geometry, on the other hand, the coherent spatial structures of the two unstable modes, while appearing similar to the coherent features on the symmetric geometry, now dominate one wing each, as presented in figure 8. The figure shows the magnitude of the unsteady surface pressure coefficient of the two unstable global shock-buffet modes, labelled b in figure 7, and a visualisation of coherent structures through the volumetric iso-surfaces of the real part of the 
$x$-momentum 
$\widehat {\rho u}$ at values of 
$\pm 0.75$. Note that the mode with the highest growth rate (in figure 8a) shows activity on the port wing. Recall from Timme (Reference Timme2020) that those modes are discrete realisations of the continuous band of medium-wavelength modes reported on infinite swept wings (Crouch et al. Reference Crouch, Garbaruk and Strelets2019; Paladini et al. Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019a; He & Timme Reference He and Timme2021). Hence, the fluid modes 
$a$ to 
$c$ (and additional modes visualised in Timme Reference Timme2020), linked to shock buffet on a finite swept wing, all come with similar coherent structures with a spanwise wavelength/wavenumber correlated with the frequency. Note that the distinct mode labelled 
$c'$ migrated significantly, approaching the unstable region, due to the increased physical realism of simulating an asymmetric static aeroelastic deformation and we will return to discussing this mode shortly.
Comparison of eigenspectra showing (
$a$) fluid-only results for symmetric (sym) vs asymmetric (asym) static deformation at angle of attack 
$\alpha = 3.75^\circ$ and (
$b$) asymmetric fluid-only vs coupled aeroelastic results, all computed with the Arnoldi method. Modes are labelled according to Timme (Reference Timme2020), along with mode 
$c'$, which migrates into the unstable half-plane in the coupled system. For reference in (
$b$), all faint-coloured fluid–structure interaction (FSI) modes are structural modes with the red-dashed box indicating the relevant region from figure 5.
Magnitude of unsteady surface pressure coefficient 
${|\hat {C}_p|}$ and volumetric iso-surfaces of real part of 
${x}$-momentum 
${\widehat {\rho u}}$ at values of 
${\pm 0.75}$ of 
$(a)$ leading and 
$(b)$ second unstable global modes from fluid-only stability analysis on statically deformed, asymmetric geometry at angle of attack 
$\alpha =3.75^\circ$. Underlying eigenvectors are scaled to unit length with respect to the inner product, specifically 
$\langle \hat {\boldsymbol {w}}_{f},\hat {\boldsymbol {w}}_{f}\rangle =1$. The base-flow zero-skin-friction line is also shown on the surface.
4.4. Coupled aeroelastic global stability analysis
The ramifications of including an elastic aircraft structure in the shock-buffet stability analysis are now addressed. Figure 7(b) gives the eigenspectra as computed by ARPACK for the fluid-only system and for the coupled system (and should be examined in unison with figure 5). Although this figure also shows results for angle of attack 
$\alpha = 3.7^\circ$ for completeness to demonstrate the impact a fluid–structure coupled approach can have in the vicinity of instability onset, the focus here is on angle of attack 
$\alpha =3.75^\circ$. Emphasising it again, all structural modes found with the earlier pk-type flutter method (bar the aforementioned difficulties with mode 28) discussed in § 4.2 were also identified by the Arnoldi method applied to the coupled system, further confirming that the implementation is correct. Multiple interesting features can be observed. The discussion focuses on three regions of the eigenspectrum, specifically the leading shock-buffet modes b, the unstable structural modes near the angular frequency 
$\omega = 1.8$, and the unstable modes 27, 28 and 
$c'$ at approximately 
$\omega = 2.7$.
First, the two unstable shock-buffet modes, denoted b, are also observed in the coupled system, albeit at slightly decreased frequencies and increased growth rates. This suggests that including an elastic structure could destabilise an otherwise stable fluid-only mode at a reduced angle of attack (Nitzsche et al. Reference Nitzsche, Ringel, Kaiser and Hennings2019). The results indicate that this change in onset angle of attack is small (particularly when compared with other typical factors having an impact on this type of simulation, such as turbulence modelling). However, it is important to point out that the migration of the leading buffet-related modes for the current test case is rather rapid (cf. figure 6 in Timme Reference Timme2020). Specifically, the leading modes 
$b$ only emerge from the dense cloud of spurious non-descript modes, to become identifiable from a global stability analysis, at an angle of attack of approximately 
$\alpha =3.6^\circ$. Also observe, as exemplified for the near-onset angle of attack 
$\alpha =3.7^\circ$ in our test case, that the leading shock-buffet modes 
$b$ are influenced in their migration by structural mode 26. These shock-buffet modes again emphasise a separate wing each, as pictured in figures 9(a) and 9(b), and the surface flow characteristics (not included in the plot for reasons of clarity) and the volumetric iso-surfaces are very similar to their fluid-only counterparts in figure 8. An interpretation could be that the shock-buffet unsteadiness remains the dominant physics even when including the elastic wing structure. Having said this, the structural part of the eigenvectors shows highest activity (i.e. deformation) on the same wing as the coherent spatial flow features, possibly revealing a coupling effect. To be clear, the unsteady wing deformation is visualised through the relation 
$\hat {\boldsymbol {x}}_s={\boldsymbol{\varPhi} }\,\hat {\boldsymbol {\eta }}$ and the real part of the dominant 
$z$-component (predominantly normal to the wing surface) is shown in the figure. To aid reproducibility of the results, when plotting the structural part of eigenvectors in physical space herein, the underlying vector 
$\hat {\boldsymbol {\eta }}$ was scaled so that the phase of its modal degree-of-freedom with the highest magnitude was forced to be zero. In addition, the plot of the MAC (cf. figure 6) describing the composition of the structural part of these fluid modes indicates a strong contribution from, and hence coupling with, most of the wind-off structural mode shapes (bar the aforementioned fuselage and empennage modes). The largest contributions come from modes 19 to 24, having wind-off frequencies close to the lower shock-buffet frequency range. Consequently, the structural motion described by these coupled fluid modes is a mix of several wind-off structural mode shapes that combine in a non-trivial manner. When visualised over a period of oscillation via 
$\tilde {\boldsymbol {x}}_s(t) = \boldsymbol{\varPhi} \,\hat {\boldsymbol {\eta }}\, {\rm e}^{\lambda t}+\text {c.c.}$ (with c.c. denoting the complex conjugate eigensolution, cf. figures 9(a) and 10a), the wing deformation resembles a stationary oscillation predominantly along the trailing edge towards the outboard region of the large scale coherent flow structures. The differences in deformation magnitude between port and starboard wings must be interpreted together with the static deformation in figure 2. Likewise, the adjoint modes (shown as a representative slice at constant span for the leading eigenmode in figure 9d) are again very similar to their fluid-only counterpart. As is often found with adjoint modes in external aerodynamics, coherent features reveal both a strong upstream support compared with their corresponding direct mode (in figure 9c) and little spatial overlap. A triangular structure is visible, resembling the observations for the span-periodic modes on infinite wings (Paladini et al. Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019a; He & Timme Reference He and Timme2020). Similarly, for the two-dimensional adjoint aerofoil mode (and accordingly the span-uniform mode on the infinite wing) it was argued that the oblique lines, impinging near the shock foot and therefore likely to be important in the global dynamics, coincide with so-called characteristic lines (Sartor et al. Reference Sartor, Mettot and Sipp2015).
Visualisation of (
$a$) leading and (
$b$) second shock-buffet modes of fluid–structure coupled system (labelled b in figure 7) at angle of attack 
$\alpha =3.75^\circ$. Surface contours show real part of deformation in 
${z}$-direction derived from structural part 
${\hat {\boldsymbol {\eta }}}$ of coupled eigenvector, while volumetric iso-surfaces illustrate real part of 
$x$-momentum 
${\widehat {\rho u}}$ at values of 
${\pm 0.75}$. Underlying direct eigenvectors are unit length with respect to the inner product, specifically 
$\langle \hat {\boldsymbol {w}},\hat {\boldsymbol {w}}\rangle =1$, whereas the adjoint eigenvector additionally satisfies bi-orthonormality, specifically 
$\langle \check {\boldsymbol {w}},\hat {\boldsymbol {w}}\rangle =1$. Slices of (
$c$) direct mode at dimensionless span 
$\eta = 0.66$, (
$d$) adjoint mode at 
$\eta = 0.55$ and (
$e$) momentum-only wavemaker at 
$\eta = 0.576$ are also shown. The inset of (
$e$) shows iso-surfaces of the wavemaker at values of 
$\theta _{f}=5\times 10^2$, 
$1\times 10^3$ and 
$1\times 10^4$. All slices include the base-flow sonic line (solid black). The base-flow zero-skin-friction line is also shown in (
$a$), (
$b$) and inset of (
$e$).
Visualisation of structural part of leading shock-buffet mode 
$b$ showing imaginary part of deformation in 
$z$-direction for (
$a$) direct (cf. the real part in figure 9a) and (
$b$) adjoint mode. Corresponding structural wavemaker is given in (
$c$). The base-flow zero-skin-friction line is also included for orientation.
Figure 9(e) visualises the fluid wavemaker (computed with only the momentum components of the direct and adjoint solution vectors) for the leading shock-buffet mode as a slice at constant span alongside three volumetric iso-surfaces over the wing surface. Discernible sensitivity to localised feedback can be pinpointed in (and above) the separation region behind the shock foot, having relatively little activity around the shock wave itself, with quite a narrow span extent overall. Close inspection of the inset plot shows the highest value of the wavemaker, 
$\theta _{f}=1\times 10^4$, right at the shock foot where shock-induced boundary-layer separation initiates. This is an important result of our study, showing the wavemaker related to shock buffet on a finite wing for the first time and extends the insight gained from previous aerofoil and infinite-wing studies (Paladini et al. Reference Paladini, Beneddine, Dandois, Sipp and Robinet2019a). This also reinforces recent arguments trying to establish the active localities involved in the instability. Paladini et al. (Reference Paladini, Marquet, Sipp, Robinet and Dandois2019b) identified the shock foot as the core of the instability and the separated boundary layer and the shock front as essential in the dynamics. Moise et al. (Reference Moise, Zauner and Sandham2022) bring arguments that feedback loops involving wave propagation mechanisms, as suggested by Lee's model (Lee Reference Lee1990), are not at the heart of aerofoil shock buffet and question the role of the shock wave as either passive or active in combination with the separated boundary layer. The structural wavemaker is shown in figure 10(c) with the direct and adjoint structural deformations given in figure 10(a,b). To be precise, for visualisation purposes, what is shown in figure 10(c) is a modification of the structural wavemaker, defined earlier for modal coordinates in (2.10), using physical Cartesian coordinates 
${\boldsymbol {x}}_s=(x_{s},y_{s},z_{s})$ instead, specifically,
\begin{equation} \theta_{s,i} = \|\boldsymbol{\mathsf{R}}_i\check{\boldsymbol{x}}_s\|_2 \, \|\boldsymbol{\mathsf{R}}_i\hat{\boldsymbol{x}}_s\|_2 + \|\boldsymbol{\mathsf{R}}_i\check{\boldsymbol{u}}_s\|_2 \, \|\boldsymbol{\mathsf{R}}_i\hat{\boldsymbol{u}}_s\|_2, \end{equation}
where the diagonal matrix 
$\boldsymbol{\mathsf{R}}_i$ extracts the direct and adjoint structural solution at a given spatial location 
$i$, and the vector 2-norm is indicated by 
$\|{\cdot }\|_2$. Considering the values of the respective wavemakers, this would suggest that the buffet phenomenon can be most effectively influenced by providing feedback through the fluid degrees of freedom. This would agree with the observation that the shock-buffet unsteadiness does not require structural vibration to initiate or self-sustain. Having said this, the spatial distribution of the structural wavemaker is interesting in that it emphasises the trailing edge region of the wing approximately where the large scale coherent flow perturbations are located. However, trying to relate the wavemaker to the stiffness and mass properties of the underlying wing structure is a rather intricate discussion. Nevertheless, our identification of the fluid and structural wakemakers could also be useful for realising effective control strategies for the buffet phenomenon, and the interested reader is referred to e.g. D'Aguanno, Schrijer & van Oudheusden (Reference D'Aguanno, Schrijer and van Oudheusden2019) and Sartor et al. (Reference Sartor, Minervino, Wild, Wallin, Maseland, Dandois, Soudakov and Vrchota2020).
Second, in figures 5 and 7(b), for the group of modes near angular frequency 
$\omega = 1.8$, the Arnoldi method does indeed identify three unstable structural modes, specifically modes 19, 20 and 21, also found by the earlier flutter analysis. The coherent flow features of these modes’ eigenfunctions (not explicitly shown herein for reasons of brevity) closely resemble those of the shock-buffet modes. Having said this, as noted earlier, a correlation between the frequencies of shock-buffet modes and their spatial amplitudes (specifically the spanwise wavelength of coherent flow features) was reported in Timme (Reference Timme2020), which is also seen in the flow features of the structural modes. For instance, the spanwise wavelengths are larger than those of the shock-buffet modes b, in accordance with their lower frequencies. In figure 6, the entries of the MAC for these modes are largest in their respective underlying wind-off structural modes, as is expected, while also revealing a strong coupling activity with neighbouring modes.
Third, in figure 7(b), for the group of eigenmodes near angular frequency 
$\omega = 2.7$, one must first be able to identify and distinguish the different modes unambiguously, considering that the pk-type Schur complement method was unable to trace mode 28 (cf. the discussion surrounding figure 5). The coherent flow features of eigenmodes 27 and 28 are similar to the shock-buffet modes b (and indeed modes 
$c$ and 
$c'$), specifically the one emphasising the port wing, but have smaller spanwise wavelengths, as is expected according to their respective higher frequencies (Timme Reference Timme2020). The entries of the MAC in figure 6 are dominated by the corresponding wind-off modes. This strong diagonal dominance in the MAC for the structural eigenmodes no matter the angle of attack has been observed throughout, e.g. it is also the case for the second group of interesting (unstable) modes near angular frequency 
$\omega = 1.8$. Hence, even though the Schur complement method failed for one mode, when traced from wind-off conditions, the success using the Arnoldi method clearly identifies it as the missing mode 28, originating in the structural system. Mode 28 is visualised in figure 11(a). Besides the now familiar coherent flow features, the structural contribution in the eigenvector mainly describes, in agreement with the discussion on the MAC, the wind-off structural mode seen in figure 1(e) that accounts for wing twist in the outer span stations. Mode 27 (not shown in figure 1) comes with a strong deformation on the horizontal tail plane. For complementary insight into these unstable structural modes such as a frequency syncing, the reader is referred to the unsteady time-marching fluid–structure coupled simulations on the same test case, while focusing both on the initial growth of disturbances and the nonlinear saturation, in Belesiotis-Kataras & Timme (Reference Belesiotis-Kataras and Timme2021). Turning the attention to the last unstable mode in the third group of interesting modes, the entries of the MAC for mode 
$c'$ (last column in figure 6) are largest for wind-off structural modes 28 and 29, suggesting it is indeed a fluid mode. However, its proximity to the structural modes gives it a notably different behaviour in the MAC than the other fluid modes, with a larger emphasis on the nearby structural modes. Compared with the other fluid modes, which all show a non-trivial combination of wind-off modes to isolate the structural activity where the flow unsteadiness sits (cf. figure 9a,b), mode 
$c'$ reveals dominant deformation from mode 28 on the starboard wing and a somewhat stronger deformation from mode 29 on the port side, as seen in figure 11(b). Note that wind-off structural mode 29 appears anti-symmetric to mode 30 shown in figure 1( f). The fluid mode 
$c'$, which has now migrated into the unstable half-plane due to strong fluid–structure interaction, has coherent flow features around the port wing and these are very similar to those of its fluid-only counterpart (not shown herein). Besides our earlier statement that the coupling of aerodynamics with an elastic wing structure could destabilise the global flow field earlier, a second finding in this regard is the ability to destabilise additional, otherwise stable, fluid modes.
Visualisation of (a) structural mode 28 and (b) unstable mode 
$c'$ of the coupled system. Variables and plotting styles are identical to those in figure 9.
To further quantify the strength of fluid–structure coupling in the different direct eigenmodes, we took inspiration from a visualisation idea presented in Negi, Hanifi & Henningson (Reference Negi, Hanifi and Henningson2021). In figure 12, the coupling ratio defined as 
${\langle \hat {\boldsymbol {w}}_s,\hat {\boldsymbol {w}}_s\rangle }/\langle \hat {\boldsymbol {w}},\hat {\boldsymbol {w}}\rangle$ (with 
${\langle \hat {\boldsymbol {w}}_s,\hat {\boldsymbol {w}}_s\rangle }= \hat {\boldsymbol {x}}_s^H \boldsymbol{\mathsf{K}}\hat {\boldsymbol {x}}_s + \hat {\boldsymbol {u}}_s^H \boldsymbol{\mathsf{M}} \hat {\boldsymbol {u}}_s$, cf. § 2.2) is plotted for a selection of eigenmodes in the shock-buffet frequency range. The ratio defines the weight of the structural and fluid components in the eigenvector and can also be used to identify whether a mode originates in the structural or fluid system. Generally, if its value is low, it is a fluid-dominated mode. Vice versa, if it is high, the structural components dominate suggesting an origin in the structural system. Importantly, considering the variety of possible vector norms that can be used, the coupling ratio is best interpreted for a chosen norm while looking at the trend for a changing significant parameter, such as angle of attack in our study. The figure reveals interesting features for the different angles of attack presented. Note, results for angle of attack 
$\alpha = 3.7^\circ$ are similar to those at 
$\alpha = 3.75^\circ$ and therefore not shown here. First, at angles of attack 
$\alpha = 3.0^\circ$ and 
$3.5^\circ$, both below buffet onset, all relevant modes in the spectrum have relatively high values in the coupling ratio, confirming that they are indeed structural modes and no strong coupling is present. Second, near buffet onset, all structural modes have approximately two to three orders of magnitude lower values in the coupling ratio due to more weight in the fluid entries. The physically relevant fluid modes 
$a$ to 
$c$ (only discernible for angle of attack 
$\alpha =3.75^\circ$) have the lowest coupling ratio, indicating they are indeed fluid modes with relatively lower weight in the structural entries. Since the projection of the unstable shock-buffet (fluid) eigenmodes onto the structural degrees of freedom is nevertheless non-trivial, one can expect an aeroelastic response. Third, fluid mode 
$c'$ at angle of attack 
$\alpha =3.75^\circ$ appears in the ratio similar to its neighbouring structural modes, distinguishing it from the other fluid modes. The strong coupling of fluid and structure, which eventually causes the destabilisation of this mode, further shows itself in this visualisation. Also, while the applied vector norm must be assumed, comparing with the incompressible pitching aerofoil results at transitional Reynolds numbers presented in Negi et al. (Reference Negi, Hanifi and Henningson2021) that indicate several orders of magnitude separation in the ratio between high-frequency fluid and low-frequency aeroelastic modes, our results would suggest stronger fluid–structure coupled interactions with fluid and structural modes in close frequency proximity.
Coupling ratio for coupled modes at various angles of attack. Fluid and structural modes are denoted with half- and fully filled markers, respectively.
5. Conclusion
The interaction between an elastic wing structure and the shock-buffet phenomenon that is a self-sustained flow unsteadiness even in the absence of wing vibration was investigated herein. For this purpose the necessary modifications to the linear harmonic incarnation of the industrial DLR-TAU solver were outlined that implemented a fluid–structure coupled formulation enabling its global mode computation using the shift-invert implicitly restarted Arnoldi method as outer iteration together with an inner sparse iterative linear equation solver. The critical component of the inner Krylov method is effective preconditioning that accounts for the discipline coupling and for that a novel approach based on the manipulation of a so-called block-arrowhead matrix was selected that offers speed up gains in parallel computation needed for a stiff aircraft case in turbulent transonic flow. These code extensions permitted the first computation of dominant fluid modes related to shock buffet alongside modes originating in the structural system, even in close frequency proximity, when a conventional pk-type flutter analysis can fail. The chosen academic test case was the well-known high-speed NASA CRM with available pressure and deformation data from a previous experimental wind-tunnel test campaign.
Validating the simulation set-up for a free-stream Mach number of 0.85, a chord Reynolds number of 
$5\times 10^6$ and several angles of attack below and in the vicinity of shock-buffet onset, static aeroelastic deformations together with the steady base-flow pressure distributions showed good agreement with measured wind-tunnel data, and those results were used as base state for all subsequent linearised analyses. Importantly, an asymmetric wing deformation (with respect to the fuselage centre plane), as a consequence of using a finite-element structural model of the actual wind-tunnel geometry, was observed, which effectively resulted in breaking up the symmetric/anti-symmetric pairs of full-span modes found in a previous study using almost perfect symmetric wing deformation. For routine flutter analysis, generalised aerodynamic influence coefficient matrices were computed with the linear harmonic solver with respect to modal structural excitation at distinct frequencies. It was shown how a conventional pk-type flutter method can fail tracing all structural modes unambiguously, owing to a strong interaction of the structural degrees of freedom with a passing fluid mode. The importance of having access to a coupled eigenmode solver, such as the one presented herein, was hence demonstrated. Consequently, the coupled approach succeeded in identifying all modes, both fluid and structure, no matter the flow condition. Near shock-buffet onset the overall dynamics becomes very active with several unstable structural modes appearing in the characteristic frequency range of the flow phenomenon. With the calculation of direct/adjoint pairs of eigenmodes, the core of the instability, colloquially called (fluid) wavemaker, was located right at the shock foot and its downstream separated boundary layer. This observation supported recent conclusions from the literature for buffeting aerofoils and infinite wings while at the same time extending those ideas to a finite wing. Similarly, the best region to introduce local feedback in the structural system points, somewhat intuitively, to the wing trailing edge at the span location of the buffet unsteadiness.
Overall it can be said that including the elastic structure adds physical realism to a high-fidelity aerodynamic analysis with the leading-order effects of the resulting aeroelastic study being asymmetric wing deformation, the potential for earlier shock-buffet onset and richer fluid–structure coupled dynamics. Otherwise pure fluid modes show a non-trivial projection onto the structural degrees of freedom in their coupled counterpart, hence giving an explanation for structural response, such as buffeting, resulting from unsteady flow. The work underpins the interest of the technical community towards more discipline coupling in this edge-of-the-flight-envelope regime.
Acknowledgements
An earlier version of this work was presented as paper AIAA 2021-0610. We thank both the German Aerospace Center and Airbus for access to the TAU flow solver and FlowSimulator tools (the latter used for static aeroelastic coupling). We also thank colleagues at the European Transonic Windtunnel (M. Schulz) for providing the data for the confidence intervals in figure 2(b). The simulation data that support the findings of this study are available from the authors upon reasonable request. The CRM geometry and finite-element model are openly available at https://commonresearchmodel.larc.nasa.gov/. Last but not least, we thank the University of Liverpool for computing time on the high-performance computing system.
Funding
The first author is grateful for the financial support by an Engineering and Physical Sciences Research Council (EPSRC) Industrial CASE scholarship (grant number 2271972 as part of EP/T517574/1) in partnership with Airbus. We also acknowledge support by the EPSRC through grant EP/R037027/1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Verification of coupled eigenmode solver
For debugging and verification purposes, the Goland wing test case was introduced. This wing is cantilevered with a constant 4 % thick, symmetric, parabolic-arc aerofoil and a (semi-)span of 20 ft and chord of 6 ft. The structure is modelled by a finite-element method, with details to be found in Beran et al. (Reference Beran, Khot, Eastep, Snyder and Zweber2004), and, through a free-vibration analysis, the four lowest-frequency wind-off modes are kept for the study. Those mode shapes mapped to the aerodynamic surface mesh and structural frequencies are provided in figure 13(a) for the sake of completeness. Two hexahedral fluid meshes are considered; one coarse mesh with approximately 20 000 points and one finer mesh with 400 000 points. Inviscid Euler flow at a reference free-stream Mach number of 
${M = 0.845}$ and zero degree angle of attack is assumed. The target altitude (to establish air density and velocity according to the standard atmosphere) for the aeroelastic analysis is 30 000 ft. A slip boundary condition is imposed at the solid wing surface, while that at the wing-root plane is symmetry and the far-field boundary is described as free-stream flow through a characteristic boundary condition.
Goland wing test case showing (a) free-vibration mode shapes and (b) part of the eigenspectrum. The eigenspectrum is computed by a pk-type method, showing the trace with respect to altitude and the eigenvalues at target altitude of 30 000 ft, compared with implicitly restarted Arnoldi method (IRAM) as implemented in DLR-TAU through coupling with ARPACK library. An exact p-type analysis, without approximations on the linearised aerodynamic response, is also shown.
Solely for debugging purposes, the relevant aeroelastic eigenmodes were first computed on the coarse mesh using three different methods, specifically the herein introduced coupled eigenmode solver (running both in sequential and parallel executions), the Schur complement approach and, after exporting the coupled matrix, MATLAB functions for sparse matrices (which also make use of the implicitly restarted Arnoldi method with shift-invert transformation, but with direct methods for matrix inversion), and the agreement was excellent throughout. The same modes were computed on the finer mesh, using the Schur complement approach as both a pk- and p-type method, with details discussed in Bekemeyer & Timme (Reference Bekemeyer and Timme2019), as well as our fluid–structure coupled inner–outer iterative method. Results are shown in figure 13(b). Note that direct matrix inversion is ruled out for the bigger mesh for obvious reasons. For the ARPACK calculation, 25 outer iterations and 1 restart were set with a convergence criterion of 
$10^{-6}$. The Krylov space for the inner iterations was 150 with 20 deflation vectors and a convergence criterion of 
$10^{-10}$. This was done twice at shifts 
$\zeta =0.05+0.15i$ and 
$0.05+0.45i$, corresponding to the vicinity of the wind-off structural frequencies. For the Schur complement approach, the p-type analysis was performed to check the pk-type results. To offer some explanation, a p-type analysis differs from a pk-type analysis in that it does not simplify the aerodynamic response to be simple harmonic but instead allows the aerodynamics to be included for the correct damping value (i.e. growth/decay rate). Essentially, the interaction matrix is evaluated as 
$\boldsymbol{\mathsf{S}}^c(\lambda )$ instead of 
$\boldsymbol{\mathsf{S}}^c(\omega )$. This ensures more accurate tracing as the method can properly represent the aerodynamics when departing from the imaginary axis. Indeed, the agreement between the different methods in the figure, specifically for modes 2 and 4 (having a larger distance to the imaginary axis) when comparing the new coupled solver and the p-type analysis, supports the notion that the coupled eigenmode solver was implemented correctly, opening up the tool to larger, more practical cases.
Appendix B. Details of arrowhead preconditioner
When using the fully coupled Jacobian matrix in parallel, block-Jacobi preconditioning discards the coupling matrices, 
$\boldsymbol{\mathsf{J}}_{sf}$ and 
$\boldsymbol{\mathsf{J}}_{fs}$. Therefore, incorporating these matrices in the formulation grants a better approximation of the inverse. This is possible by utilising an identity presented by Stanimirović, Katsikis & Kolundžija (Reference Stanimirović, Katsikis and Kolundžija2019), which gives the inverse of a block-arrowhead matrix based on the Sherman–Morrison–Woodbury formula.
Preprocessing stage of arrowhead preconditioner
Block-arrowhead matrices have the sparsity pattern of an arrow (hence the name)
\begin{equation} \boldsymbol{\mathsf{P}} = \begin{pmatrix}\boldsymbol{\mathsf{A}}_{1} & 0 & 0 & \dots & \boldsymbol{\mathsf{B}}_{1} \\ 0 & \boldsymbol{\mathsf{A}}_{2} & 0 & \dots &\boldsymbol{\mathsf{B}}_{2} \\ 0 & 0 & \ddots & & \vdots \\ \vdots & \vdots & &\boldsymbol{\mathsf{A}}_{n} & \boldsymbol{\mathsf{B}}_{n} \\ \boldsymbol{\mathsf{C}}_{1} & \boldsymbol{\mathsf{C}}_{2} & \dots & \boldsymbol{\mathsf{C}}_{n} &\boldsymbol{\mathsf{D}}, \end{pmatrix}\end{equation}
resembling our specific problem, where we define 
$\boldsymbol{\mathsf{A}}_{i} = (\boldsymbol{\mathsf{J}}_{ff, i}-\zeta \boldsymbol{\mathsf{I}})$ (i.e. the local diagonal blocks of the shifted fluid Jacobian matrix 
$\boldsymbol{\mathsf{J}}_{ff}$) and accordingly 
$\boldsymbol{\mathsf{B}}_{i} = \boldsymbol{\mathsf{J}}_{fs, i}$, 
$\boldsymbol{\mathsf{C}}_{i} = \boldsymbol{\mathsf{J}}_{sf,i}$ and 
$\boldsymbol{\mathsf{D}} = (\boldsymbol{\mathsf{J}}_{ss}-\zeta \boldsymbol{\mathsf{I}})$. Subscripts 
$i = 1,\dots, n$ denote the process number. Here, although all off-diagonal blocks of 
$(\boldsymbol{\mathsf{J}}_{ff} -\zeta \boldsymbol{\mathsf{I}})$ are discarded like before for the fluid-only problem, all elements of the other matrices are kept. The exact analytical inverse of 
$\boldsymbol{\mathsf{P}}$ can then be derived as
\begin{align} \boldsymbol{\mathsf{P}}^{{-}1} &=\begin{pmatrix} \boldsymbol{\mathsf{A}}_{1}^{{-}1} & 0 & 0 & \dots & 0 \\ 0 &\boldsymbol{\mathsf{A}}_{2}^{{-}1} & 0 & \dots & 0 \\ 0 & 0 & \ddots & & \vdots\\ \vdots & \vdots & & \boldsymbol{\mathsf{A}}_{n}^{{-}1} & 0 \\ 0 & 0 & \dots &0 & 0 \end{pmatrix}\nonumber\\ &\quad +\begin{pmatrix}\boldsymbol{\mathsf{A}}_{1}^{{-}1}\boldsymbol{\mathsf{B}}_{1}\\ \boldsymbol{\mathsf{A}}_{2}^{{-}1}\boldsymbol{\mathsf{B}}_{2}\\ \vdots \\\boldsymbol{\mathsf{A}}_{n}^{{-}1}\boldsymbol{\mathsf{B}}_{n}\\ -I \end{pmatrix} \boldsymbol{\cdot} \boldsymbol{\mathsf{F}} \boldsymbol{\cdot} (\boldsymbol{\mathsf{C}}_{1} \boldsymbol{\mathsf{A}}_{1}^{{-}1}\quad \boldsymbol{\mathsf{C}}_{2} \boldsymbol{\mathsf{A}}_{2}^{{-}1} \quad \dots\quad \boldsymbol{\mathsf{C}}_{n} \boldsymbol{\mathsf{A}}_{n}^{{-}1} \quad -I ),\end{align}
where matrix 
$\boldsymbol{\mathsf{F}}$ is defined as
\begin{equation} \boldsymbol{\mathsf{F}} = \left(\boldsymbol{\mathsf{D}} -
\sum_{i=1}^{n} \boldsymbol{\mathsf{C}}_{i}\, \boldsymbol{\mathsf{A}}_{i}^{{-}1}\,
\boldsymbol{\mathsf{B}}_{i}\right)^{{-}1},
\end{equation}
and 
$\boldsymbol{\mathsf{I}}$ is the identity matrix, both with the same (small) dimensions as matrix 
$\boldsymbol{\mathsf{D}}$. Computing the factors 
$\boldsymbol{\mathsf{A}}_{i}^{-1}\boldsymbol{\mathsf{B}}_{i}$ in (B2) essentially requires applying the ILU factorisation of the shifted fluid Jacobian matrix, 
$(\boldsymbol{\mathsf{J}}_{ff,i}-\zeta \boldsymbol{\mathsf{I}})$, to the 
$2m$ columns of matrix 
$\boldsymbol{\mathsf{J}}_{fs,i}$ (for each process locally). However, this needs to be done only once and for all, so computing 
$\boldsymbol{\mathsf{Y}}_i=\boldsymbol{\mathsf{A}}_{i}^{-1} \boldsymbol{\mathsf{B}}_{i}$ and 
$\boldsymbol{\mathsf{F}}$ can be done as a preprocessing step, described in Algorithm 1. The additional memory requirements are for matrix 
$\boldsymbol{\mathsf{Y}}$ of size 
$\boldsymbol{\mathsf{J}}_{fs}$ and matrix 
$\boldsymbol{\mathsf{F}}$ of size 
$\boldsymbol{\mathsf{J}}_{ss}$. Preconditioning an arbitrary vector 
$\boldsymbol {v}$ is explained in Algorithm 2. Compared with applying the block-Jacobi variant, besides operations on the negligible size of the structural system, only two additional matrix-vector products (one with matrix 
$\boldsymbol{\mathsf{J}}_{sf,i}$ and one with matrix 
$\boldsymbol{\mathsf{Y}}_i$), a vector–vector addition of the size of the fluid domain and one global sum of size of the structural problem are needed, which makes the computational overhead acceptable, with details provided below. Note, since the structure-to-fluid coupling matrix 
$\boldsymbol{\mathsf{J}}_{sf}$ is very sparse with the only non-zero entries coming from the surface points where the generalised forces are integrated, matrix-vector products with this matrix are relatively cheap. Similarly, the fluid-to-structure coupling matrix 
$\boldsymbol{\mathsf{J}}_{fs}$ is not a dense matrix as it results from residual evaluations on the deformed mesh using a radial basis function tool which linearly reduces the applied deformation to zero within a specified distance from the wall. Hence, disturbed volume mesh points (and consequently non-zero residuals) are confined to the vicinity of the wall.
Application stage of arrowhead preconditioner
The performance of different preconditioning approaches is now discussed. The speed up gained by using the block-arrowhead preconditioner compared with the block-Jacobi preconditioner (easily implemented by setting all elements of the coupling blocks, 
$\boldsymbol{\mathsf{J}}_{fs}$ and 
$\boldsymbol{\mathsf{J}}_{sf}$, in the preconditioner to zero) was insignificant for the inviscid Goland case without strong shock waves (see Appendix A). However, for the larger aircraft case with turbulent shock-wave/boundary-layer interaction, the inclusion of the arrowhead preconditioner gave significant speed up when compared with the block-Jacobi preconditioner. Figure 14 shows such a comparison of performance using the preconditioned GCRO-DR iterative solver algorithm for 20 successive inner linear solutions required by the shift-invert Arnoldi method. The convergence criterion was set to 
$10^{-7}$ with 
$120/20$ Krylov/deflation vectors. Two shifts, 
$\zeta = 0.05+2.63i$ and 
$2.63i$, were chosen with the latter deliberately defining a challenging problem due to its proximity to an eigenmode for angle of attack 
$\alpha =3.7^\circ$, but nonetheless represents a typical set-up encountered in many simulation scenarios and exemplifies the benefits of the arrowhead preconditioner.
Typical convergence behaviour of preconditioned GCRO-DR iterative solver (with 120 Krylov and 20 deflation vectors) showing 20 linear solution histories for NASA CRM eigenvalue problem at angle of attack 
$\alpha = 3.7^\circ$ with shifts (a) 
$\zeta = 0.05+2.63i$ and (b) 
$\zeta = 2.63i$, contrasting arrowhead and block-Jacobi preconditioners.
Figure 14(a) presents a benign scenario for shift 
$\zeta =0.05+2.63i$. Although the block-Jacobi method does converge, the arrowhead preconditioner shows consistent superior performance on average in terms of iteration count. This behaviour is also reflected when comparing computation times. A single application of the arrowhead preconditioner takes approximately 
$50\,\%$ longer when compared with block Jacobi due to the additional operations involved. However, a complete linear solution to the specified tolerance is roughly 
$25\,\%$ faster, due to preconditioning only accounting for a small part of the total computation time in each Krylov iteration with most of the cost coming from the orthogonalisation of the Arnoldi vectors, increasing with the size of the subspace. This speed up in computation time therefore matches the reduction in the average number of iterations needed to reach convergence. Figure 14(b) for shift 
$\zeta =2.63i$ not only demonstrates clear performance gains of the arrowhead preconditioner, but more importantly its robustness. Down to a convergence level of 
$3\times 10^{-6}$, the discipline-coupled preconditioner takes roughly half the number of iterations needed by block Jacobi, largely due to the superior rate of convergence in the initial iterations. Beyond that, the block-Jacobi preconditioner is entirely inadequate, as illustrated through the stalled iterations. Additionally, in our investigation, the block-Jacobi preconditioner never outperformed the arrowhead preconditioner, even when performance gains were low. While we show next that an inner convergence level in the range 
$10^{-7}$ to 
$10^{-5}$ is often sufficient to achieve reasonably converged results in the outer Arnoldi iteration, the robustness of the arrowhead method is preferred and was hence the default choice in our study.
Appendix C. Assessment of inner–outer iterative Krylov method
An assessment of convergence properties of the inner–outer iterative Krylov method follows next. Table 2 summarises the impact of the specified convergence level for the iterative Krylov linear solver on both the relative error of the eigenmode, 
$|1-\lambda /\lambda _{10^{-9}}|$, and the resulting outer residual norm, 
$\|\boldsymbol{\mathsf{J}}\hat {\boldsymbol {w}} - \lambda \hat {\boldsymbol {w}}\|$, to ensure the accuracy of the results, subject to the inherent iterative error of the methods exercised herein. The remaining settings are identical to those in § 4. The two chosen eigenmodes are the fluid mode labelled 
$c'$ and hard-to-converge structural mode 28, shown in figures 5 and 7. It is found that the residual norm is roughly an order of magnitude higher than the tolerance of the linear solver. Also, as usual, the frequency converges quicker than the growth rate. A tolerance of 
$10^{-3}$ is clearly insufficient, as the precision required on the eigenvalues, especially the growth rate, is greater than what can be reached. A tolerance of 
$10^{-7}$ to 
$10^{-5}$ is sufficiently precise for an engineering accuracy. The ultimate decision would depend on the scope of a specific study. For instance, choosing a tolerance of 
$10^{-9}$ could be useful to demonstrate consistency between different methods but would add little added insight into the physics that warrants the greatly increased computational cost. From these observations, a tolerance of 
$10^{-7}$ was selected as the best option.
Impact of inner solution tolerance on outer convergence at angle of attack 
$\alpha =3.75^\circ$ using shift 
$\zeta = 2.7i$. The eigenvalues correspond to fluid mode 
$c'$ and hard-to-converge structural mode 28 (cf. figures 5 and 7). The relative error is shown for the growth rate and frequency separately and calculated, e.g. for the real part, as 
$|1-\text {Re}(\lambda )/\text {Re}(\lambda _{10^{-9}})|$, with 
$\lambda _{10^{-9}}$ denoting the solution with tolerance 
$10^{-9}$ and bold decimal places indicating unconverged digits with respect to 
$\lambda_{10^{-9}}$.
Open access
