Advanced aeroelastic robust stability analysis with structural uncertainties

Robust flutter analysis described in this paper is based on the robust control theory framework. Therefore, a time-domain linear fractional transformation representation of the perturbed aeroelastic system is modeled. Then, the robust stability is analyzed by means of the structured singular value μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu$$\end{document}, which is defined as an alternative measure of robustness. Robust flutter analysis deals with aeroelastic (or aeroservoelastic) stability analysis taking structural dynamics, aerodynamics and/or unmodeled system dynamics uncertainties into account. Flutter is a well-known dynamic aeroelastic instability phenomenon caused by an interaction between structural vibrations and unsteady aerodynamic forces, whereby the level of vibration may trigger large amplitudes, eventually leading to catastrophic failure of the structure. The primary motivation of the robust flutter analysis is that this method allows the computation of the worst-case flutter velocity which can support, for example, the flight test program by a valuable robust flutter boundary. This paper addresses the issue of an approach for aeroelastic robust stability analysis with structural uncertainties with respect to physical symmetric and asymmetric stiffness perturbations on the wing structure by means of tuning beams.


Introduction
The primary aim of this paper is to investigate the impact of the stiffness uncertainties of the wings in spanwise direction and handle each wing separately in case of symmetric and especially asymmetric stiffness distribution. The asymmetric stiffness perturbation is mentioned in [3] within the robust flutter analysis for aeroelastic systems as an important issue for future investigations. It is quite possible that poor levels of precision with respect to manufacturing capabilities a small difference of bending and/or torsional stiffness may occur between the two wing structures which -in worst-case scenario-can be significant enough to cause an unpredictable coupling between symmetric and asymmetric modes. Therefore, it is essential to model such an asymmetry by means of an adequate uncertainty in physical stiffness model for each wing separately. This paper starts with a brief mathematical introduction on the LFT representation of uncertain systems and provides an understanding of the definition and interpretation of within the robust stability analysis.

Linear fractional transformation and -analysis
The linear fractional transformation is a common framework for robust stability analysis of complex systems based on the small gain theorem [1]. An LFT is an interconnection of operators arranged in a feedback manner. Using a linear complex partitioned operator the LFT, ( , Δ) is defined as the closed-loop transfer matrix from system input u to system output y as the upperloop LFT of system closed by the norm-bounded structured matrix of perturbations Δ ( Fig. 1):

3
An equivalent expression ( , Δ) is defined as the closedloop transfer matrix from system input u to system output y as the lower-loop LFT of system closed by Δ as illustrated in Fig. 2: It can be seen from the Eq. (2) that denotes the transfer function (matrix) from input u to output y of the nominal plant while , and content the information how the perturbation is embedded in the nominal plant. After the uncertainty block Δ is extracted from the nominal plant, the system can be described as a nominal plant with an artificial feedback-loop. From the Eq. (2) it becomes apparent that the LFT is wellposed if and only if the inverse of ( − Δ) exist [4] (Fig. 3).

Structured singular value
The structured singular value is an exact indicator of robust stability for systems with structured uncertainties (perturbations). It is a function of the complex transfer function matrix and ̄(⋅) denotes the maximum singular value of the argument [2] (2) In this context is robustly stable with respect to which is nor m-bounded by scalar ∈ ℝ such that ‖Δ‖ ∞ ≤ , ∀Δ ∈ if and only if ( ) ≤ 1 . In the -framework the model is usually weighted to normalize the normbounded uncertainty set Δ to unity For ≤ 1 there is no perturbation within exists that will destabilize the system. This state depicts that the true system dynamics are stable, assuming the nominal model dynamics with its set of uncertainty operators (modeling errors) are able to capture the dynamic behaviour of the true system.

Aeroelastic model
In this paper a condensed FE model of the FLEXOP demonstrator aircraft [5] for the numerical demonstration of the robust flutter analysis is used. The full FE model ( > 600000 nodes) comprises the wing, fuselage and empennage. The wing is modeled by a high-fidelity FE model comprising beam, surface and solid elements. The condensation of the FE model has been performed by the Guyan-reduction, also known as static condensation. The reduced model consists of 303 structural nodes and 1818 degrees of freedom (DoF). The aerodynamic model is based on vortex lattice method (VLM) [11] for steady aerodynamics and doublet lattice method (DLM) [10] for unsteady aerodynamics. A detailed overview of both the structural and aerodynamic model of the aircraft (Fig. 4) is described in [5].

Equations of motion for the aeroelastic system
The equations of motion for the nominal aeroelastic system in time domain can expressed in a matrix form as which describes a system of N linear ordinary differential equations (ODEs) with N degrees of freedom (DoF) in the FE model where x(t) is the displacement vector with translational and rotational DoFs of the nodes, , and if and only if ( ) ≤ 1.  The AIC matrix can be computed by several aerodynamic theories, such as DLM. In this paper, the subsonic unsteady aerodynamic forces have been modeled by means of DLM. Based on small disturbance hypothesis, DLM solves the linearized potential flow equation and obtains the aerodynamic forces under the assumption that aerodynamic surfaces oscillate harmonically. The nondimensional Laplace variable s is denoted s = g + ik where g is the damping and k is the reduced frequency. On the assumption of harmonic aerodynamic loads the nondimensional Laplace variable s becomes: where is the frequency of vibration and c ref the reference chord. Note that the dependence on the Mach number of the AIC matrix will be omitted from now on for conciseness. Using mode displacement method the physical displacement vector x(t) can be represented as a linear combination of m linearly independent vectors (mode shapes) leading to the approximation where m is the number of eigenvectors. Combination of the Eqs. (6), (7) and (8) results in the following reduced-order dynamics with where m is the generalized mass matrix, m is the generalized viscous damping matrix, m is the generalized stiffness matrix and m (ik) the generalized AIC matrix.

Rational function approximation
The generalized AIC matrix m (ik) ∈ ℂ m x m is a set of matrices which are calculeted for a set of suitable values of reduced frequency k. Thus, in order to compute AIC for any desired reduced frequency and perform time domain analysis (statespace representation), the AIC matrix in the frequency-domain has to be transformed into the Laplace domain and consequently into the time domain. One possible way is to fit the frequency dependent AIC matrix with rational functions in a least-squares sense. This method is called Rational Function Approximation (RFA). In this paper the Roger's formulation [6] is used to approximate the AIC matrix m (ik): The RFA Eq. (14) can be interpreted as a general two-part approach for aerodynamic loads based on quasi-steady and lag contributions. 0 m , 1 m and 2 m are ℝ m x m real coefficient matrices representing the contribution of acceleration, velocity and displecement of the rigid and flexible degrees of freedom on the aerodynamic loads denoting the quasi-steady part of the approximation. The L i m ∈ ℝ m x m matrices with the predefined poles i , i = 1, 2, ..., n p , are responsible for the lagging behavior of the unsteady flow. This is referred to time lag effect.
For time domain representation the equation system in (14) can be rearranged as follows [7] where (10)  For the lag states x L ∈ ℝ n lag x 1 the following ODE with ̇ as input can be derived [7]: The resulting generalized aerodynamic forces are then An alternative new approach for the generation of a generalized state-space aeroservoelastic model based on tangential interpolation is presented in [12]. Compared to RFA approach, this method enables a minimal order realization with interpolation of the unsteady aerodynamics by avoiding any selection of poles i [12].

LFT model of the perturbed system
Consider again the generalized equations of motion for the aeroelastic response of the aircraft with unstaedy aerodynamics by combining the Eqs. (9) and (20): The state-space respresentation of the Eq. (21) is formulated with the generalized states , ̇ and the unsteady aerodynamic states x L :

Parametrization over flight speed
The -framework determines the stability over a range of airspeed to specify the onset of flutter. The generalized equations of motion for the aeroelastic response (22) are a function of the flight speed V, such that perturbations of this parameter can be integrated into the system as a linear fractional transformation.
Consider an additive perturbation V on the nominal velocity V Separate the terms in the system dynamics (22) that involve V: The additional input and outputs signals are introduced into the nominal aeroelastic system given in (22) to include the perturbation in velocity to the nominal dynamics in a feedback manner (Fig. 5): The nominal flutter problem modeled by means of perturbation to flight speed V in Eq. (28) can be solved as a -computation. The tranfer function matrix (s) which relates the input signal w V to the output signal z V can be derived from the Eq. (28): To make the perturbation on flight speed V which has the same unit as V , to unity norm-bound constraint ‖ V‖ ∞ ≤ 1 , the transfer function matrix (s) has to be scaled by a weighting W V , where The scaled plant transfer function matrix (s) is given by: and in frequency domain with s = i where (i ) is also called the flutter loop transfer function matrix of the uncertain system (Fig. 6).
Now the nominal flutter speed V nom, flutter can be determined via -computation by means of the following algorithm.
For < 1 there is no perturbation within exists that will destabilize the closed-loop system. Beacuse of well-known difficulties of exact real/mixed computation, which poses an NP-hard problem, efficiently computable upper bound of is used within the analysis.  The algoritm set out above should therefore give an introduction into the -framework. The nominal flutter speed V nom flutter can also be calculated by means of standard flutter solution techniques such as p−k-method or g -method. For example, the p−k-method [8] is widely used in flight flutter test campaigns. The g-method proposed by Chen [9] promises a reliable damping prediction by including a first order of damping term in the flutter equation. For the robust aeroelastic stability analysis additional uncertainties such as in structural dynamics, aerodynamics and/or unmodelled system dynamics have to be defined and subsequently included in the linear system to introduce modeling errors between the numerical model and the physical aircraft (Fig. 7).

Tuning beams
In this paper, account is taken of the uncertainties in the structural dynamics or, more explicitly, the structural stiffness model which has a significant impact within the lower frequency range and flutter analysis respectively. For the realization of the stiffness parameter variations tuning beams have been generated with respect to the condensed FE model. This approach is suitable for varying of stiffness parameters by only adjusting material properties of the tuning beams like bending and torisonal stiffness while avoiding a modification of the full FE model. The tuning beams are 3D beams and may also be called "space beam" elements [13].
The beam element resists force in any direction and moment about any axis. The element stiffness matrix el ∈ ℝ 12x12 of the 3D beam in local coordinate system is given by [13]: where the submatrices 11 el , 12 el , 21 el and 22 el are defined as follows: where E is elastic modulus, G is shear modulus, A is crosssectional area, I yy and I zz are principal moments of inertia of A, I p is the torsional constant and L is the length of the 3D beam element. The assembly of global stiffness matrix of the tuning structure requires that each of the element stiffness matrix has to be transformed into the global coordinate system and subsequently integrated into the global stiffness matrix of the nominal FE model in the following way: where tuned is the physical stiffness matrix of the tuned FE model, is the nominal physical stiffness matrix of the FE model and beam i is the physical stiffness matrix of the i th tuning beam. This type of representation of the stiffness matrix in Eq. 51 by superposition of nominal model and the Fig. 7 Arbitrarily oriented 3D beam element with 12 nodal degrees of freedom (DoF) in local and global coordinate system [13] 1 3 tuning structure is used for the uncertainty definition in the next section. The global stiffness matrix with its substructure components (nominal model and tuning structure) is visualized in Fig. 8.

Uncertainty modeling
In this study, the wing structure is affected by the uncertainty of the physical stiffness matrix by means of tuning beams' parameters. Two kinds of variations of physical stiffness characteristics have been defined. The first approach is characterized by symmetric variation of physical stiffness with respect to torsional rigidity ( GI p ) and/or flexural rigidity (bending stiffness) ( EI yy ), i.e. the tuning beams on the right and left wing structure have the same material properties to include symmetric uncertainties into the nominal model. The second approach denotes asymmetric variation of stiffness parameters of the tuning beams on the right and left-wing structure to represent asymmetric stiffness distribution in the nominal model.
Considering the perturbation of the physical stiffness matrix, the parametric additive uncertainties can be described as follows: It should be noted that is the modal (eigenvector) matrix of the nominal system. This assumption is reasonable for small perturbations and can be validated by modal correlation analysis between nominal and perturbed system.
The extended LFT model of the new perturbed system, which includes feedback signals relating the perturbation to flight speed and uncertainties described in the Eqs. (52) and (53), may be defined analogous to the Eq. (28): where the new submatricest K ∈ ℝ (2m+n lag ) x (m⋅n beam ) and K ∈ ℝ (m⋅n beam ) x (2m+n lag ) are given as follows T h e s u b m a t r i c e s V,K ∈ ℝ (3m+n lag )x(3m+n lag ) , K,V ∈ ℝ (m⋅n beam )x(3m+n lag ) and K ∈ ℝ (m⋅n beam )x(m⋅n beam ) are zero matrices.
The robust flutter problem modeled by means of perturbation to flight speed V and stiffness model uncertainty in Eq. (54) can be solved by a -computation. The tranfer function matrix K (s) which relates the input signals w V and w K to the output signals z V and z K can be derived from the Eq. (54): Once the transfer function matrix K (s) has been determined, the robust flutter speed V rob flutter can be determined within the -framework by means of the following algoritm.

Numerical results
For the numerical demonstration of the proposed uncertainty description 11 tuning beams for each wing ( n beam = 22 ) have been defined. Each wing consists of 60 structural nodes. The beams are placed within the nondimensional span range = [0.25, 0.75] . For the modal model the first 15 modes ( m = 15 ) have been considered within the frequency range of 2.9 ≤ f ≤ 35.6 Hz. In this context, the considered frequency grid for -computation must be sufficiently dense in order not to miss a critical frequency, which can lead to a prohibitive computational cost. For the following case studies, a frequency grid of 50 logarithmically spaced points between 3.0 and 12.0 Hz has been considered with regard to the nominal flutter results. The computed values are then interpolated within the same frequency range using linear interpolation with a step size of Δf = 0.05 Hz.
Symmetric stiffness perturbations I Additive uncertainty of 10 2 N m 2 with respect to bending stiffness EI yy for each tuning beam on the left and right wing structure

Analysis results
The numerical example demonstrates that even small stiffness parameter variations on the wing structures have a major impact on the onset of the flutter. There are considerable differences between flutter speeds with respect to bending and torsional stiffness variations. The impact of bending stiffness variations on the flutter margin is much stronger than the torsional stiffness uncertainties. This was to be expected due to the flutter results of the nominal model as shown in the graphic below (Fig. 16) where normalized modal participation factors of the flutter mode are plotted. The flutter mode involves a strong coupling between the first wing bending and the first wing torsion mode (critical modes) whereas the first wing bending mode at f = 2.93 Hz dominates the motion at the flutter condition. Consequently, bending stiffness variations have a higher effect on the onset of the flutter condition. The defined uncertainty in Case I reduces the flutter speed by roughly 8% compared to the nominal flutter speed given in Table 1 whereas the symmetric uncertainty of the torsional stiffness defined in Case II leads to a reduction of the flutter speed by roughly 6.2% . The Case III depicts a combination of the previous two cases and leads to a further reduction of the flutter speed by roughly 3.8% compared to the Case I. A further important point relating to key characteristics of the results refers to the asymmetry of the stiffness distribution as an uncertainty integrated into the nominal model. Comparison of the flutter speeds between Cases I and IV leads to the conclusion that an asymmetric bending stiffness uncertainty has an greater influence on the flutter margin then an assumption of symmetric bending stiffness uncertainty (reduction of flutter speed by 10.2% ), whereas the flutter speeds in Cases II and IV related to the symmetric and asymmetric torsional stiffness uncertainties remain relatively unchanged. Comparison of the results between Cases III and VI also reflects that a combination of asymmetric stiffness uncertainties related to bending and torsion parameters has a stronger effect on the flutter margin in comparison to the assumption of symmetric stiffness uncertainties in the model. Fig. 9 Structural nodes of condensed FE model with tuning beams

Conclusions
In this work a robust aeroelastic stability analysis within the -framework has been carried out. Therefore, an LFT model of the perturbed aeroelastic system in state-space form parametrized over flight speed has been developed.
Here, the account is taken of the uncertainties in the structural stiffness model which has a significant impact within the lower frequency range and flutter analysis respectively. For the realization of the stiffness parameter variations tuning beams have been generated with respect to the condensed FE model. This approach is suitable for varying of stiffness parameters by only adjusting material properties of the tuning beams avoiding to intervene in the full FE model. The study is limited to the wings of the aircraft and focussed on the investigation of physical stiffness uncertainties of the wings in spanwise direction and handle f rob