Model Updating and Aeroelastic Correlation of a Scaled Wind Tunnel Model for Active Flutter Suppression Test

Abstract

This article presents a modal correlation and update carried out on an aeroelastic wind tunnel demonstrator representing a conventional passenger transport aircraft. The aim of this work is the setup of a corresponding numerical model that is able to capture the flutter characteristics of a scaled aeroelastic model designed to investigate and experimentally validate active flutter suppression technologies. The work described in this paper includes different finite element modeling strategies, the results of the ground vibration test, and finally the strategies adopted for modal updating. The result of the activities is a three-dimensional hybrid finite element model that is well representative of the actual aeroelastic behavior identified during the wind tunnel test campaign and that is capable of predicting the flutter boundary with an error of 1.2%. This model will be used to develop active flutter suppression controllers, as well as to perform the sensitivity analyses necessary to investigate their robustness.

1. Introduction

The use of scaled wind tunnel models to validate simulation results is well known in the field of aeroelastic predictions for both military and civil transport aircraft. The activity is not limited to the research investigations phase but can also support and mitigate the risk of the certification phase based on full-scale flight tests. This is particularly true in the case of unconventional configurations whose aeroelastic behavior is more difficult to predict. The scaling process, as well as the adopted strategies to realize aeroelastically scaled wind tunnel models, has been completely described by Bisplinghoff in [1]. The traditional approach is typically based on a condensation of the mass and stiffness distributions of the considered aircraft to those of beams representing the elements of the model with attached masses and separated aerodynamic sectors that do not appreciably contribute to the overall bending or torsional stiffness. Thus, the entire process is mainly reduced to the definition of the properties of the scaled beam-like structures to match the actual (scaled) stiffness distribution of the reference aircraft well. It is clear that, with this approach, the actual structural topology of the aircraft is not reproduced, and the final configuration of the scaled model is mainly driven by the manufacturing technologies, by the material properties, and by the need to install the sensors and actuators necessary to investigate the relevant aeroelastic phenomena.

The current tendency for modern aircraft to show enhanced flexibility, as well as to enable aeroelastic deflections to be used in a beneficial manner [2], makes the use of scaled aeroelastic wind tunnel models more important. On the other side, the adoption of more and more complex Active Control (AC) systems for aeroelastic control, together with the need to validate their performances during the wind tunnel test campaign, opens the door to new manufacturing technologies, such as 3D printing and alternative materials [3], but at the same time makes the level of aeroelastic correlation between the scaled model and the reference aircraft even more important. This is the case for the example of flutter [1], which is a dynamic instability that results from the coupling of aerodynamics and structural dynamics. It represents a complex phenomenon and, at the same time, a design constraint that is complicated due to its dependence on the coupling of close structural modes, and thus the mass and stiffness distributions within the structure. The level of representativity of the scaled model becomes a key factor, especially in the case where not only the phenomenon but also its active control is under investigation. The current advanced state of the art of AC sees Active Flutter Suppression systems as the main solution to this kind of instability [4]. Since flutter is an aggressive and abrupt kind of instability, it is urgent to figure out how precise and robust the AC law has to be. Its design and implementation rely strongly on the accuracy of the model representing the real behavior of the aircraft.

In this broad context, the activities described herein fit the development of a well-correlated mathematical model that is able to fully characterize the phenomenon and that can be used as the basis for the design of control laws. As noted in [4], the move from theoretical–numerical studies or studies based on highly simplified wind tunnels or flight models to the full implementation of the technology requires more experimental work. The investigation should comprehensively consider wind tunnels and in-flight tests, using scaled models that represent real full-size aircraft in terms of their configuration and complexity, as this would help identify technological deficiencies and development needs, especially from the perspective of safety and the associated uncertainty–reliability simulation and the test capabilities required. Many works can be found in the literature on the numerical/experimental correlation of wing structures and updates to airplane models. An interesting work on the aeroelastic scaling and optimization of a joined-wing aircraft is described in [5,6]. Some relevant work with a scope similar to that of this manuscript was carried out during the US-based Performance Adaptive Aeroelastic Wing (PAAW) project [7] and the European Flutter Free Flight Envelope Expansion for Economical Performance improvement (FlexOP) project [8,9,10], while a summary of both projects can be found in [11].

The Department of Aerospace Science and Technology (DAER) at Politecnico di Milano (PoliMi), in collaboration with the University of Washington, contributed to the so-called Active Flutter Suppression (AFS) project that aimed to validate active flutter suppression technologies. On the basis of its large amount of experience during the last 30 years in the design, manufacturing, and testing of wind tunnel aeroelastic models [12,13,14], DAER developed a dedicated wind tunnel model, called X–Dia, to be used for the AFS project and focused on these technologies, including uncertainties and robustness issues [15]. This paper reports the activity required to set up the numerical aeroelastic model and in particular the numerical vs. experimental aeroelastic correlation and updating phase. Traditionally, this activity is formulated as an optimization problem [16]. The state of the art of optimization software is extraordinarily advanced; many commercial optimizers allow very complex problems to be considered, finding the best minimum solution in the infinite set with a user-friendly and automated process. In this case, the design sensitivity and optimization of MSC/NASTRAN has been adopted [17].

The work is set out as follows: Section 2 presents the X–Dia model, in terms of its conceptual design, realization, and testing procedures to characterize its dynamic and aeroelastic behavior. The large experimental campaign has made measures available with which researchers can study the flutter mechanism and tune the mathematical model. Section 3 summarizes the results of the experimental modal analysis activities representing the targets for the updating process. Section 5 describes the process adopted to gain the final FEM model that optimally describes the real structure and its flutter onset, based on experimental data, proper modeling considerations, and the optimization algorithm. This last aspect is deeply considered in Section 6, outlining the general rules, the modeling constraints, and the critical points found in the process of the updating phase. Section 7 contains the comparison results of the modal and flutter analyses between the experimental data and the mathematical model predictions and, in light of the correlation evaluation methods, the necessary steps to achieve the optimum. In particular, through different sub-sections, we provide the evolution of the updating process of the FEM model based on an optimization algorithm to achieve an optimum correlation level and accurately determine the flutter boundary. The procedure stands out in the multi-stage development of the model due to the ever-increasing level of detail and structural interventions according to the aeroelastic response requirements. Finally, Section 8 presents the conclusions on the work that has been carried out and a brief reminder of the obtained results in terms of a comparison between the numerical and experimental flutter velocity.

2. The X–Dia Wind Tunnel Model

The X–DIA model was initially conceived by PoliMi to study different topics concerning flight mechanics, aeroelasticity, and active controls. Its configuration was subject to variations according to the intended purpose. The final model shows a more conventional layout, with a low positive swept wing and a T-tail without the presence of the canard, (the original configuration is reported in [12,13], while the final model is shown in Figure 1) to meet the aeroelastic demands of the activities. The model is representative of a generic aircraft, but it is not a scaled version of a specific aircraft. The design specifications were related to a generic configuration and wind tunnel test requirements.

The skeleton of the X–Dia model is composed of aluminum beams connected to each other. The fuselage beam part is in a tubular structure with a rectangular aluminum section. The wing spar is a beam with an omega section shape, while the vertical and horizontal ones are cross-shaped. Connections of structural elements are designed and milled from aluminum blocks. The wing spar is covered with discrete aerodynamic sectors, made by 3D printing technology using Windform composite SLS (Selective Laser Sintering) materials (developed and patented by CRP Technology) to guarantee high stiffness, so to avoid local modes. Both the wing spar and aerodynamic sectors are shown in Figure 2a. They are designed to allow internal space for the installation of the instrumentation. Two milled ribs support the aileron hinges and the servo-actuators, which are based on the Harmonic Drive electric motors (RSF-5B) coupled with an in-house developed servo-control law [18].

X–Dia is equipped with an on-board computer to manage the active control systems and data acquisitions. The embedded computer is equipped with an in-house developed Real Time Application Interface (RTAI), a real-time operating system based on Linux. The embedded computers, as well all the necessary acquisition cards, are based on the PC-104 form factor allowing for a small final size and weight of the complete embedded system, which is shown in Figure 2b. The main sensors are MEMS accelerometers by PCB, which were selected for their high accuracy at low frequencies and mainly because of being already pre-amplified. The corresponding channels could be read by a standard data acquisition system without the need of expensive and heavy IEPE modules, usually adopted in case of piezo-electric accelerometers. This choice allowed for a small embedded computer to be installed inside the model that is capable of managing the test in terms of acquisition and active control duties. By having the complete instrumentation on-board, the model is almost in a free–free condition inside of the wind tunnel, which is a beneficial circumstance for this kind of test since a wind tunnel is by definition a noisy environment.

Among the latest changes made to the new X–Dia model is that a safety device was installed at the tips of the wing, as shown in Figure 3, in view of the flutter wind tunnel tests. This system allows us to change the mass configuration so to eliminate the flutter instability. The device contains a small mass of 200 gr able to slide inside a pod from two on–off positions (respectively forward and backward) thanks to a linear pneumatic actuator. The slider length is about 15 cm. According to the position of the sliding mass, the aircraft is prone to flutter or not. The forward mass position has been designed to increase the flutter velocity above the maximum velocity obtainable in the wind tunnel. During the flutter wind tunnel tests, the mass is positioned in the backward position. Then, the maximum vertical acceleration of the tip is monitored. When a threshold value is exceeded the system is automatically switched on to move the mass in the forward position in both pods to immediately stop the flutter and to avoid damaging the model. The threshold value has been setup during the preliminary tests to 4g, a value found to be the best compromise to see a complete flutter development without reaching dangerous vibration levels.

The tip sensors data are processed separately, so that the left and right safety pods can work independently as they are activated only by local acceleration, practically developing two independent co-located anti-flutter devices. This strategy offers the possibility to face the onset of several kinds of flutter (symmetric, non-symmetric, body-freedom flutter, etc.).

During the wind tunnel test, carried out at PoliMi’s Large Wind Tunnel, the aircraft is free–free suspended in the wind tunnel by means of an elastic connection made by cables so that the structure can experience the rigid body modes with a frequency less than 1 Hz.

3. Ground Vibration Test (GVT)

The experimental campaign had a double purpose: identify the actual dynamic properties of the model as well identify the flutter mechanism. The first target has been achieved by a dedicated Ground Vibration Test (GVT) session [19,20,21]. With the evolution of the project, the mass and in some cases the stiffness distribution of the model changed due to the installation of different components. In particular, the following items have impacted the dynamic properties of the model: the installation of reinforcing rods to increase the stiffness out of plane of the original horizontal tail; the on-board data acquisition system, which heavily influences the model mass and CG position;and the flutter stopper safety devices that were installed on both the wing tips. Different GVT sessions have been carried out for tracking the linear dynamic properties and for periodically updating the intermediate numerical models adopted [19,20,21]. The two most relevant configurations in view of numerical vs. experimental correlation and model updating were selected, with the anti-flutter safety device installed on both the wing tips and the moving mass in rear position or in forward position. For these configurations, two different sensors setup have been adopted:

Setup 1:

all the lifting surfaces—i.e., wings, horizontal planes, and vertical tail—have been measured in the out-of-plane direction. A roving hammer impact test has been carried out, limiting the mass loading effect due to the mass of the accelerometers. Five fixed accelerometers have been installed upon the structure and used as references, while the other points have been used for excitation and computed with 10 averages, for a total of 46 measurements.

Setup 2:

the left wing has been instrumented with 28 fixed accelerometers that are equally distributed to measure in both in and out-of-plane direction. The excitation has been introduced by the hammer impact. This setup has been realized mainly to investigate deeply specific aspects, such as the dynamic behavior of the pod.

Two mass configurations have been considered, depending on the anti-flutter safety device:

1.

Safety On: the aircraft with the mass of both anti-flutter devices fully forward.

2.

Safety Off: the aircraft with the mass of both anti-flutter devices fully backward.

The aeroelastic demonstrator inside the wind tunnel chamber during the GVT is shown in Figure 4a, while the detailed experimental mesh adopted for the Setup 1 is shown in Figure 4b. A database of Frequency Response Functions (FRFs) has been acquired considering a sampling frequency of 512 Hz and a frequency resolution of 0.0625 Hz, while the anti-aliasing filters cut frequency is 80% of the Nyquist frequency, which is set to 200 Hz. The FRFs have been computed using the $H_v$ method to eliminate the influence of uncorrelated noise in the input and the output signals [22]. The wind tunnel configuration adopted during the GVTs is the same one used during the flutter tests, to reduce as much as possible the differences in the dynamic behavior of the model due to the suspension system.

Once the measurements were collected, a preliminary validation phase was planned to substantiate the quality of the measurements. Although many verification criteria are available in the literature [23], the Modal Assurance Criterion (MAC) was chosen [24]. The experimental validation was completed using Modal Phase Collinearity (MPC), Mean Phase Deviation (MPD), and Mode Overcomplexity Value (MOV) indexes.

Table 1. Comparison among frequencies with tip masses in forward and rear positions (BEN = bending; TOR = torsion; S = symmetric; AS = asymmetric).

ID	Mass forward (Hz)	Mass Rear (Hz)	ΔFrequency (%)	Shape
1	5.20	5.20	0.00	1st S-BEN
2	7.14	7.43	4.06	Left wing 1st TOR
3	7.78	7.88	1.29	Right wing 1st TOR
4	8.26	8.28	0.24	1st S-BEN htail + 1st S-BEN wing
5	9.00	9.10	1.11	1st AS-BEN htail + 1st S-BEN wing
6	10.43	10.38	−0.48	1st S-BEN htail + 1st S-BEN wing
7	12.55	12.49	−0.48	1st S-TOR htail torsion + 1st S-TOR wing
8	13.45	13.46	0.07	1st S-BEN htail
9	14.98	14.74	−1.60	2nd AS-BEN wing and 2nd AS-BEN htail
10	16.36	16.28	−0.49	2nd S-TOR wing + 2nd S-BEN htail
11	19.69	19.20	−2.49	2nd AS-TOR wing + 2nd AS-BEN htail +
				+ 1st AS-BEN vtail
12	24.12	23.35	−3.19	2nd AS-TOR wing + 2nd S-BEN htail

and

Table 2. Comparison among damping with tip masses in forward and rear positions.

ID	Mass forward (Hz)	Mass Rear (Hz)	ΔDamping (%)
1	0.13	0.17	30.77
2	1.22	1.14	−6.56
3	1.08	0.75	−30.56
4	1.03	1.30	26.21
5	1.65	3.07	86.06
6	0.48	0.39	−18.75
7	1.00	0.90	−10.00
8	2.65	2.73	3.02
9	1.60	1.23	−23.13
10	0.83	0.85	2.41
11	0.54	0.57	5.56
12	0.65	0.65	0.00

show the comparison among frequencies and damping measured with the tip masses in forward and rear positions, respectively. Figure 5 shows the comparison between measured mode shapes with tip masses in forward and rear positions represented in the form of a Cross Mac matrix. From the results it is possible to draw the following conclusions:

• As expected, the different position of the tip masses has a limited impact on the frequencies and damping.
• Looking at the Cross MAC matrix it is possible to see that the differences in the mode shapes are limited to the low frequencies mode shapes and mainly the wing torsional modes, due to the position of the moving masses positioned on the wing tip.
• Based on the GVT data, it appears that the presence of a light non-symmetry in the real model makes it difficult to measure separate global symmetric and anti-symmetric torsional modes. They appear as two different modes at about the same frequency, involving each wing half separately. This generated a poor correlation in the Cross-MAC matrices, as reported in Figure 6b,c.

This unusual behavior is the main source of difficulty in producing a final well-correlated model. However, it has been decided not to try to completely recover the symmetry of the model, nor to try emulating the same behavior in the numerical model, but rather to leave a non-perfectly matched torsional behavior to test the robustness of the active control laws for flutter suppression.

After the GVT, a preliminary wind tunnel test session was scheduled, which aimed at the identification of the flutter behavior, before proceeding with the detailed modal updating phase that is based on both GVT and experimental flutter results.

4. Preliminary Experimental Flutter Test

The wind tunnel flutter test campaign has been carried out in the PoliMi’s large Wind Tunnel, suspending the X–Dia model in a free–free configuration by means of elastic cables. The scope of this test campaign is to identify the flutter boundary keeping the overall aeroelastic trend under control, while the velocity increases.

The sensor layout is shown in Figure 7b and includes 20 embedded MEMS monoaxial accelerometers in addition to one capacitive triaxial accelerometer located in the CG of the complete model. The test session was remotely managed from the wind tunnel control room by an internet cable connected to the on-board system and laid along the front suspension cable. Additionally, two cables for the computer and the servo-actuators power supply were used. Figure 7a shows the X–Dia model installed inside the wind tunnel chamber.

For the assigned wind-tunnel velocities reported in

Table 3. Wind tunnel test speeds for each configuration.

Speeds (m/s)
Backward mass	0	10	15	20	25	30	35	37.5	40
Forward mass	0	10	15	20	25	30	35		40

, both anti-flutter mass configurations have been considered, and the complete model has been self-excited by different kinds of input excitation produced by oscillating the ailerons: sine sweep, impulse, band-limited white noise, and square wave. The constant band-limited white noise allowed us to perform a preliminary real-time flutter identification to demonstrate that the model is flutter-free until the maximum planned velocity when the anti-flutter device is turned on.

The next wind tunnel test session was aimed at the fully experimental flutter identification, considering when the anti-flutter device was turned off. Two different methodologies were adopted. The first one is based on the use of PolyMAX method, implemented in the LMS Test-lab [20,25], which is the same approach adopted during the GVT. The results, in terms of frequency and damping values, were also computed using an in-house code, based on the ARX approach [26].

Once the maximum speed was reached, the aeroelastic system appeared to be close to a zero-damping condition, so the velocity incremental step was reduced to 0.5 m/s. The wind tunnel test continued looking for the flutter point that was achieved at 41.5 m/s.

5. Structural Models

The need to have numerical models that correctly reproduce the features of interest can be traced back to the possibility of simulating the real structure, predicting the results of the analyses, as well as developing, designing, and tuning its response. Even if a large and complete set of tools is available, creating a well correlated mathematical model is not a linear and straightforward process but an iterative operation and the result of a compromise. This is especially true for wind tunnel aeroelastic models that are typically realized by single beams to which many and different elements are connected, such as the aerodynamic sectors, the local stiff elements like ribs to support the electric motors driving the control surfaces, local masses represented by sensors, as well as the cables. The goal of the correlations is based on a comparison among measured and predicted quantities, able to fully capture the aeroelastic behavior, such as the modal frequencies and mode shapes, as well the classical $V-g$ plots, while the natural choice to numerically model this kind of aeroelastic wind tunnel models, especially where the load carrying member is a beam, by classical stick model, sometimes it happens that they are not able to fully capture the actual aeroelastic behavior. This is what happened in the case of the X–Dia model that was forced to adopt a hybrid model in place of a simple stick model.

Some hypotheses have been introduced in the construction of the model to simplify the subsequent analyses without excessively compromising the final results. At the beginning, a classical stick model based on beam elements (Timoshenko formulation, including warping DoF, for the wing spar modeling; Euler–Bernoulli for the fuselage), concentrated masses, rigid elements, and aerodynamic behavior of a lifting surface, described with the Doublet Lattice Method (DLM), has been built. These hypotheses are justified in the perspective of achieving a correlated model. The efforts have been headed to catch a sufficient detail level to allow the FEM model to describe the global trend. The properties of beam elements have been extracted from analytical section calculations, combined with experimental and CAD measures. In addition, only certain structural modes are involved in the flutter mechanism and for this reason, they are ones that are worth focusing on during the correlation and the updating process.

During the correlation phase, involving frequencies, mode shapes, and $V-g$ plots comparisons, it appeared that the need to increase the fidelity of the model to capture the effect of local structural variations cannot be captured by the simple stick model, as will be further described in Section 7. Therefore, a new model, called hybrid, was implemented.

One of the difficulties in modeling the actual aeroelastic behavior by means of a simple stick model is due to the use of an open, omega section for the wing spar. Both the poor torsional stiffness typical of the open sections and the location of elastic axis outside the wing plane, due to the omega-shape, generating a strong coupling between the in-plane motion and the torsional one, confer to the X-DIA wing dynamic a significant complexity. Despite the possibility in MSC/NASTRAN to use beam elements controlling the position of the elastic center of the spar section, an adequate tuning phase has been requested to identify the correct value for the shear stiffness factor so as to correctly model the bending-torsion coupling.

Furthermore, structural elements with marked 3D characteristics, such as the connection between the wing and the fuselage, have shown a dynamics that may not be grasped with the stick model. The previous considerations, based on the preliminary experimental correlation phase, have highlighted the need to modify the initial stick model by modeling the wing spar by plate elements and the wing–fuselage connection by solid elements. The result is a hybrid model, as represented in Figure 8, with different level of detail depending on the actual impact of these zones on the overall dynamic behavior.

Another local modification has been necessary to correctly model the details of the connection between the wing and the anti-flutter pods. It consists of a concentrated mass linked with fairly flexible connections to a rigid section. The difference in stiffness of the components, together with its large distance from the axis of the aircraft, makes the impact of this model on the modal properties significant. To manage this issue, a dedicated elastic element has been introduced to model the connection between the wing tip and the attachment of the pod. Tuning the stiffness value of the spring allows to take into account possible compliance effects.

It has to be noted that, despite all the modifications made to the model refer to structural properties, the resulting effect is also to improve the aeroelastic correlation, due to their impact on the mode shapes and as consequence on the unsteady aerodynamic forces. Interesting and beneficial studies could embrace the aerodynamic formulation of the anti-flutter device, modeled using MSC/NASTRAN slender body formulation, but it does not fall within the goals of this work [27].

6. Model Updating Strategy

The wide and complete database, represented by the GVT results as well as by the preliminary wind tunnel flutter tests, described in Section 3 and Section 4, is the basis for the correlation process and later the updating phase. It must be taken into account that, especially for complex aeroelastic wind tunnel models, there are many components or unknown parameter values that prevent the best correlation level from being obtained. Unmodeled dynamics could be the cause for not reaching an acceptable correlation level. A typical solution to cover these discrepancies is fitting the experimental data changing model parameters in a least-square sense. This is at the core of the updating method herein employed. The modeling errors causing low correlation can be filled by the implementation of an optimization problem aiming at finding the model parameters values to minimize the numerical vs. experimental errors [28,29,30].

The so-called model updating is not straightforward and in most cases requires significant user experience. Various successful model updating methods have been developed to work on the matrices representing the dynamics of the system [29,31]. In the case of the X–DIA model, the updating problem has been formulated like a classical structural optimization process and solved thanks to the design sensitivity and optimization capabilities of SOL200 of MSC/NASTRAN [17]. This choice is mainly due to the possibility to use the same code and the same FEM models for both the updating process and the aeroelastic assessment.

6.1. Optimization Problem Setup

The optimization problem setting is based on some hypotheses:

• Many optimizers give the possibility to base the error to minimize many kinds of response variables and flutter velocity as well. Even if the scope of these activities is to capture the aeroelastic behavior with special attention on the flutter point, a purely structural-based formulation has been preferred. A restricted number of variables coming into play can guarantee a better convergence and stronger robustness of the results. The aeroelastic effects have been recovered in a second phase. Thus, the objective function is defined to minimize the differences in the natural frequencies and the differences in the mode shapes between the numerical model and the GVT results.
• The model configuration subjected to optimization has been limited only to the anti-flutter mass in backward position, which is the one prone to flutter. It would have been possible to carry out a double optimization on both configurations, but the actual difference between the two is due to the position of a single mass. In this way, the second mass configuration could be used as a check for the correctness of the optimal solution.

The optimization problem requires us to define an error function to be minimized and consists of the differences between numerical model outputs and the desired values, namely the experimental data. The function passed to the optimizer, thanks to DEQATN and DRESP2 MSC/NASTRAN cards [17], has the form:

$$J\left(f_1,f_2,f_3\right)={\left(\frac{f_{GV{T}_1}-f_1}{f_{GV{T}_1}}\right)}^2+{\left(\frac{f_{GV{T}_2}-f_2}{f_{GV{T}_2}}\right)}^2+{\left(\frac{f_{GV{T}_3}-f_3}{f_{GV{T}_3}}\right)}^2$$

(1)

No weight-sum squared function has been used, since the actual cost function contains the dimensionless relative errors between ${\mathbf{f}}_{{\mathbf{GVT}}_{\mathit{i}}}$, the measured natural frequencies, and ${\mathbf{f}}_{\mathit{i}}$, the numerical frequencies, and it has been assumed that they all have the same impact on the flutter mechanism and the same accuracy. Equation (1) is composed of the first three structural modes frequencies since they are the ones taking part in the aeroelastic instability. Indeed, the X–DIA presents a typical coupled bending–torsional flutter instability: the 1st is a bending mode, the 2nd and 3rd are two coincident torsional modes, anti-symmetric and symmetric, respectively.

One of the key aspects during model updating is represented by the choice of the design variables, i.e., the properties of the structure that are subject to modifications during the process. They have to be representative of the uncertainties resulting from the modeling process, but at the same time, they have to guarantee an updated model representative of the physical one, while few design variables allow to recover systematic error in global stiffness and mass representation, a higher number of design variables allows us to correct local modeling errors.

The design variables considered here are the thickness of the plates of the wing spar, chosen for the possibility to recover all the unmodeled details that in some way can impact on the stiffness distribution. However, to avoid unreasonable thickness changes, it is necessary to impose some user-defined constraint. The mass of the model is experimentally known and fixed, so accordingly to the thickness the density must also vary in such a way that the total mass does not change. The DLINK card of MSC/NASTRAN has provided an optimization constraint to specify the density as a function of independent variables according to an adaptation law stated by (Equation (2)):

$$\rho (t_1,t_2,t_3)=\frac{c_1}{t_1}+\frac{c_2}{t_2}+\frac{c_3}{t_3}$$

(2)

where $t_1,t_2$, and $t_3$ are thicknesses of spar parts (central and vertical webs and side flanges, see Figure 9) and $c_1,c_2$, and $c_3$ are constant coefficients, introduced by means of DTABLE MSC/NASTRAN card, linking the changes in thicknesses to the total mass of the structure, in order to have it constant. Each part of the spar can have an independent change (fixed spanwise), but the updated value of the density has to consider all of these, each with a different weight based on their portion of the overall volume.

$$c_{\mathit{i}}=\frac{V_{no{m}_{\mathit{i}}}}{V_{tot}}{t}_{no{m}_{\mathit{i}}}{\rho }_{nom}$$

where:

• $V_{no{m}_1}$ is the volume relating to part 1 of the omega-shaped spar with nominal values;
• $V_{tot}=V_{no{m}_1}+V_{no{m}_2}+V_{no{m}_3}$ is the total volume in the nominal configuration;
• $t_{no{m}_1}$ is the nominal thickness relating to part 1 of the omega-shaped spar;
• ${\rho }_{nom}$ is the density in nominal configuration of the spar material.

The scope of this optimization setting is to grasp the correct inertia allocation of the wing spar. The density constraint has been included to fix the total value of the mass and allow the algorithm to distribute the inertia independently (mass properties are shown in

Table 4. Mass properties of the FEM models.

Model	Mass (kg)	${\mathsf{x}}_{\mathbf{CG}}$ (m)	${\mathsf{y}}_{\mathbf{CG}}$(m)	${\mathsf{z}}_{\mathbf{CG}}$ (m)
Stick	59.42	−0.784	0	0.077
Hybrid	59.49	−0.777	0	0.077
Opt ver.1	59.49	−0.777	0	0.077
Opt ver.2	59.46	−0.779	0	0.077
${\mathbf{I}}_{\mathbf{xx}}(\mathbf{kg}·{\mathbf{m}}^{\mathbf{2}})$	${\mathbf{I}}_{\mathbf{yy}}(\mathbf{kg}·{\mathbf{m}}^{\mathbf{2}})$	${\mathbf{I}}_{\mathbf{zz}}(\mathbf{kg}·{\mathbf{m}}^{\mathbf{2}})$	${\mathbf{I}}_{\mathbf{xy}}(\mathbf{kg}·{\mathbf{m}}^{\mathbf{2}})$	${\mathbf{I}}_{\mathbf{xz}}(\mathbf{kg}·{\mathbf{m}}^{\mathbf{2}})$	${\mathbf{I}}_{\mathbf{yz}}(\mathbf{kg}·{\mathbf{m}}^{\mathbf{2}})$
14.52	25.65	37.17	0	2.94	−0.05
14.52	25.73	37.23	0	2.92	0
14.55	25.74	37.28	0	2.92	0
14.53	25.74	37.26	0	2.92	0

).

The differences in the mode shapes between the numerical model and the GVT data could be minimized and included in the optimization loop as well, using for example the related MAC values. The cost function weighted MAC term has the form of Equation (3).

$$J_{MAC}=W_{shap{e}_i}(1-MA{C}_{red})$$

(3)

where the $MA{C}_{red}$ is an entry of the ${\mathbf{MAC}}_{\mathbf{red}}$ matrix corresponding to a reduced set of key nodes $W_{shap{e}_i}$ is the weight matrix containing possibly different weights for the entries of the ${\mathbf{MAC}}_{\mathbf{red}}$.

Despite modal shapes having a great impact upon the flutter mechanism, this last term has not been included in the overall cost function but just used as a check after every optimization loop to track the mode shapes correlation obtained. Indeed, including the mode shapes into the objective function significantly increases the computational cost of the optimization process due to the need to face the non-linearities introduced.

6.2. Updating Phase

As already mentioned, the goal is to correctly capture the flutter behavior, but the optimization carried out involves only modal parameters (Section 6.1), i.e., modes at zero speed. When the aerodynamic forces come into play, the matrices of the system change according to the Mach number and the reduced frequency. To also include the flutter contribution within the procedure, the optimization runs have been performed while constantly monitoring the aeroelastic response. Optimization runs have been alternated with flutter analyses to check the current state and to tune the process in order to fit the experimental data. This allowed us to gradually study the flutter prediction compared to the experimental one. Thus, the union of these two techniques leads to a complete updating analysis of the aeroelastic behavior of the wind tunnel model.

The results of an updating phase strictly depends on the choice of the design variables and their allowable variation. Thus, starting from the hybrid finite element model, different configurations of design variables have been explored.

In a first step, a limited number of design variables, i.e., the thickness of the horizontal and vertical walls of the omega, kept constant spanwise, was adopted. The algorithm converged in four iterations, with a reduction of the cost function of one order of magnitude from $4.37\times {10}^{-2}$ to $6.48\times {10}^{-3}$. In a second step, a more aggressive optimization was applied: a larger number of design variables was introduced, which also impacted the stiffness of the wing–fuselage connection block, the shell’s transverse flexibility, and aerodynamic effect of the anti-flutter device. The algorithm required three iterations and reduced the cost function of another order of magnitude up to $2.06\times {10}^{-4}$. In particular, the algorithm showed that the shell transverse behavior is a key factor in representing the real model dynamics. Initially, the shell shear factor was listed as a design variable and was included in the optimization loop, but the algorithm tendency to increase it drastically suggested to the authors to not involve it in the optimization loop, but rather to consider infinite transverse rigidity shell elements as the starting point for the second optimization phase.

The model evolution is presented in Figure 10 with the variables trend in

Table 5. Design variables trend during optimization.

Design Variables	Initial Value	Final Value
1st Optimization loop
$t_1$	3.4 mm	4 mm
$t_2$	3.4 mm	4 mm
$t_3$	3 mm	2.59 mm
2nd Optimization loop
$t_1$	4 mm	4.06 mm
$t_2$	4 mm	4.06 mm
$t_3$	2.59 mm	2.61 mm
E of wing-fuselage connection block	7 GPa	4.5 GPa
E of the wingspar shell elements	7.2 GPa	7 GPa

and the obtained results are discussed in next Section 7.

7. Correlation and Updating Results

During the activities herein presented, four FEM models have been developed: a stick model (Section 5), a hybrid model (Section 5), and optimized models version 1 and version 2 (Section 6.2).

The correlation starts reporting the comparison of the low frequency modes in terms of frequencies and mode shapes.

Table 6. Frequencies comparison between the experimental modes and corresponding numerical ones of the models developed before the optimization process.

GVT #	Freq (Hz)	Stick #	Freq (Hz)	Error (%)	Hybrid #	Freq (Hz)	Error (%)
1	5.20	1	5.19	0.06	1	5.59	7.49
2	7.43	2	7.33	1.29	2	6.85	7.78
3	7.88	3	7.47	5.17	3	6.87	12.84
5	9.10	6	8.94	8.21	5	8.93	1.79
6	10.38	5	8.79	3.37	7	11.10	6.91
7	12.49	9	12.25	1.94	8	12.19	2.40
8	13.46	10	14.43	7.23	11	14.42	7.13

compares the first seven GVT frequencies with the ones obtained using the stick and the hybrid models, respectively. It is possible to see that the hybrid model better captures the dynamic properties than the stick one, even if an updating phase is still necessary due to the large frequency error present.

Figure 11 shows the Cross MAC matrix comparing the mode correlation among the free-suspended GVT mode shapes and the free-suspended numerical ones based on the stick model and the hybrid model, respectively. The choice to build the MAC matrices with the free-suspended mode shapes can be traced back to their capability to mostly represent the wind tunnel model behavior since both experimental and numerical modal analyses have been conducted with the same model configuration for the aeroelastic tests.

Starting from this situation, the updating activities aim at the modification of the structural model to improve the correlation in the rear mass configuration since they are the most important factor from the flutter behavior point of view.

Modal Correlation

Starting from the stick model results, it can be seen that after the modifications are introduced into the model, the correlation with the bending modes is strongly improved, in particular with the first mode that is the one involved in flutter instability. The experimental–numerical correlation level has been measured and evaluated by means of model validation techniques (Section 3).

Natural Frequencies Correlation

Attention has been focused on the lower frequency modes because they are the most significant for the dynamic behavior of the structure, especially for the flutter phenomenon. The wind tunnel tests have shown that the main flutter instability involves the 1st and 3rd modes, respectively, out-of-plane symmetric bending and symmetric torsional of the wing.

Table 7. Natural frequencies comparison. The presented modes involve mainly wing shapes, with anti-flutter mass in backward position.

Mode Type	GVT (Hz)	Stick (Hz)	Hybrid (Hz)	Opt 1 (Hz)	Opt 2 (Hz)
1st symm. bending	5.20	5.19	5.59	5.55	5.27
1st anti-s torsional	7.43	7.33	6.85	7.28	7.42
1st symm. torsional	7.88	7.47	6.87	7.31	7.47
Error (%)
1st symm. bending	-	0.06	7.49	6.72	1.34
1st anti-s torsional	-	1.29	7.78	2.02	0.13
1st symm. torsional	-	5.17	12.84	7.23	5.20

compares the frequencies of the first three modes of the mass-backward configuration since they are the ones used for the optimization process. The flutter results have been checked using different sets of modes, reaching a total of 40 modes to compute the response and looking at the possible variations. The result remained unaffected.

Starting from the stick model results, a very limited error can be noted (wheb dealing with modal frequencies, a threshold error of 1–2% is usually considered). Nevertheless, its prediction for flutter behavior is very bad, as described below in this section, demonstrating that a good correlation level on modal frequencies is not sufficient to also have the same level for aeroelastic ones. The transition from the stick model to the hybrid one has not produced better results concerning the the frequency errors but has paved the way for a better description in terms of mode shapes. For this reason the updating process has been applied to the hybrid model, thus setting up an optimization problem aiming at the minimization of the frequency errors of the first three low frequency modes, i.e., the first bending and two torsional ones working mainly on the stiffness of the wing spar and the wing–fuselage connection.

During the first step, updating only three model parameters has updated, i.e., the three wing spar thicknesses describing the Omega section, keeping them fixed spanwise, so to correct the wing stiffness distribution in a global fashion. In the objective function the frequencies to be reached have been set equal to $f_{GV{T}_1}=5.2$ Hz, $f_{GV{T}_2}=7.45$ Hz, and $f_{GV{T}_3}=7.45$ Hz.

As already anticipated, the reason why the frequencies $f_{GV{T}_2}$ and $f_{GV{T}_3}$ do not coincide with GVT’s results is due to the fact that the numerical model is perfectly symmetric showing two torsional modes (i.e., symmetric and anti-symmetric) with the same modal participation of left and right half wings. On the contrary, the actual physical model shows a slightly asymmetric behavior able to split any torsional mode into two separate modes with different participation on left and right side half wings (see Figure 12). This partially justifies the high error on the second torsional mode. In terms of the final design variables, the solution obtained shows a thickening of the vertical webs and side flanges of the spar, which intends to increase the torsional inertia and so the related frequencies, and a lightening of the central web that tries to decrease the bending frequency.

During the second step, a more extended model tuning was applied to improve the correlation level. The set of design variables were defined looking at the sensitivity of modal results showing that:

• Decreasing the Young modulus of the wing root connection elements allows us to decrease the first bending frequency.
• Decreasing the Young modulus or the thickness of the top flange of the spar mainly decreases the bending frequency.
• Changing the webs and flanges thicknesses of the spar impact on the torsional frequencies.

At the same time, further modification to the finite element model was applied, such as:

• Transverse shear flexibility set in order to have an infinitely rigid in transverse direction plate and trying to increase the torsional stiffness of the wing.
• Classical and not burdensome slender body aerodynamics formulation of the pod was introduced during the flutter analysis, to better reproduce the flutter behavior of the updated model.

The second updating step allowed us to improve the correlation level that will be now discussed in detail. Looking at the modal frequencies reported in Table 7 and the Cross MAC matrices ([24]) of Figure 13, it is possible to draw the following remarks:

• With the last optimization run, the error committed on the modal frequencies falls below an adequate threshold, thus achieving an excellent level.
• Starting from Hybrid model, the updating process has managed to reduce the error at every step, as can be seen in Table 7. The algorithm has impacted the target frequencies of 0.72% for the 1st symmetric bending, 6.28% for the 1st anti-symmetric torsional, and 6.4% for the 1st symmetric torsional, during the first optimization loop; the values were 5.05%, 1.92%, and 2.19% during the second loop.
• The 1st first bending mode has been captured very well, both in terms of frequency and shape. This is very important because the flutter mechanism involves the exactly 1st bending and the 2nd torsional (mode # 3).
• Due to the presence of two split torsional modes, as identified during the GVT campaign and reported Figure 12b,c, due to the fact that the actual wind tunnel model is not perfectly symmetric, the correlation results are still poor.

Despite this low correlation related to the torsion modes, the frequency error reduction as well the model improvement allowed us to capture the global aeroelastic behavior and it will be now discussed.

Concerning the flutter instability, commonly it is studied in terms of $V-g$ and $V-f$ plot typically obtained using the well known p-k method [32], where the aeroelastic damping g is the quantity tracked to identify flutter onset; the instability is associated with negative damping, reflected in a positive value of g.

In such a case, these values have been extracted through post-processing activities. No dedicated Operational Modal Analyses were carried out.

Figure 14 and Figure 15 depict the trend of the three modes involved in the flutter mechanism and compare the experimental results with the numerical ones. The preliminary wind tunnel tests had the scope to identify the flutter; no active flutter suppression controllers were available at that time. Thus, the experimental velocities were limited by the flutter point (Table 3). The error in predicting the flutter velocity is summarized in

Table 8. Flutter velocity and flutter frequency committed error comparison.

	GVT	Stick	Hybrid	Opt 1	Opt 2
${\mathbf{V}}_{\mathbf{f}}(\mathbf{m}/\mathbf{s})$	41.5	45.5	32.5	42	40.5
$\epsilon (\%)$	-	9.6	21.7	1.2	2.4
${\mathbf{f}}_{\mathbf{flutter}}\left(\mathbf{Hz}\right)$	6.2	5.73	6.47	6.41	6.41
$\epsilon (\%)$	-	7.58	4.35	3.38	3.38

. During these tests, the presence of unsteady aerodynamic forces practically made the left/right torsion difference no more significant and almost disappeared during the identification.

Looking at Table 7 and Table 8 it is possible to say that both step 1 and step 2 updated models have produced significant results in terms of error reduction with respect to the experimental flutter test results. In the first instance, the hybrid model alone did not improve the error in flutter velocity prediction but did prepare the ground for the optimization process. In fact, after the second updating loop, the final model was the best compromise between the errors in the flutter velocity and frequency prediction and the aeroelastic behavior represented in the $V-g$ plot.

Of special interest to us is the correlation obtained after the updating process in the case of the model with the mass of the safety flutter device in the forward position, which has not been directly involved in the updating. Indeed, since in reality the only difference between the two configurations is represented by the position of two concentrated masses on the tips, it is reasonable to expect some progress in the correlation level, demonstrating the global meaning of the model changes made.

Concerning the modal parameters, the modes participating in flutter refined from time to time their correlation, in particular the 1st bending mode (

Table 9. Natural frequencies comparison. The anti-flutter mass is in forward position.

Mode type	GVT (Hz)	Stick (Hz)	Hybrid (Hz)	Opt 1 (Hz)	Opt 2 (Hz)
1st symm. bending	5.20	5.22	5.43	5.45	5.22
1st anti-s. torsional	7.14	7.11	6.29	6.69	6.85
1st symm. torsional	7.78	7.14	6.57	6.9	6.98
Error (%)
1st symm. bending	-	0.38	4.4	4.8	0.38
1st anti-s. torsional	-	0.42	11.9	6.3	4.06
1st symm. torsional	-	8.23	15.55	11.31	10.28

and Figure 16). For the two torsional modes, the improvements were limited by the non-symmetry of the real model, as explained previously.

Looking at the flutter behavior reported in Figure 17 and Figure 18, the $V-g$ and $V-f$ plots have confirmed what was expected: the updating process applied to the FEM model with the safety mass in the backward positions produced an updated model able to capture the global aeroelastic behavior also in the case of the safety masses in the forward position.

8. Conclusions

This paper summarized the activity carried out to finalize the aeroelastic model of the X–DIA wind tunnel model. It proposed an overview for model updating and its role in aircraft development.

The first part of the work focused on FEM model developments. Starting from a typical stick model usually adopted for aeroelastic models, the model was improved into a hybrid model based on a plate model of the main spar that is more suitable to capture the structural behavior of the wing. This choice was mainly due to the particular omega shape of the wing cross-section that naturally couples the in-plane bending to the torsion modes. Then, solid elements were adopted to model the blocks connecting the wing spar to the fuselage. Special attention was devoted to the modeling of the anti-flutter pods due to their strong impact on the wing dynamic response.

After the correlation process with the GVT results, the updating process was started, conceived, and implemented as a typical structural optimization problem. Two different updating loops were performed, with different numbers of design variables. The correlation and updating process was complicated by the slight non symmetry of the physical model, making it impossible to reach a perfect correlation among the numerical and experimental torsion modes. Despite these difficulties, the updating process was successfully completed, thus obtaining a numerical dynamic model able to perfectly reproduce the global aeroelastic behavior measured during the wind tunnel testing, with a final error on the flutter velocity around 1.2%. The updated dynamic model will be used as the starting point to generate the state-space aeroelastic model to be used for the development of dedicated flutter suppression control strategies.