© J.M. Colomer et al., Published by EDP Sciences 2021

1 Introduction

The search for more efficient aircraft configurations leads aircraft designers to explore new aircraft concepts such as the blended wing body, the box wing, or the strut-braced wing. While the classical wing-body configuration is well known and understood, little is known about the in-flight behavior of these new aircraft concepts. In that context, the design, construction, and testing of unmanned, aeroelastically scaled models presents itself as a means of acquiring experimental knowledge on these new concepts. This mitigates the risk of manufacturing and testing a full scale aircraft both from the economical and operational point of view. Also, the use of aeroelastically scaled models allows seeing the effect of design modifications of existing aircraft — such as a new engine — on the flutter response. In this paper, we investigate the introductions of slight variations on the wing planform of flight demonstrators to minimize the discrepancies on the dynamic aeroelastic behavior despite the differences in compressibility conditions (i.e., Mach number).

Traditional aeroelastic scaling of flying demonstrators assumes that there exists flow similarity between the two scales, or that at least that the difference in the flow conditions (especially the compressibility effects) is negligible. However, it is often not possible to achieve the same Mach number (M) due to the operational limits of the unmanned flight demonstrators.

The fundamentals of classical aeroelastic scaling methods were described by Bisplinghoff et al. [1] in 1955. Scruton and Lambourne [2] extended the approach to include compressibility and heat transfer effects. French [3] proposed an optimization-based approach for matching the static response of finite element models for aeroelastic scaling purposes. To consider dynamic aeroelastic scaling, French and Eastep [4] introduced a two-step approach. First, the structure stiffness was matched through the minimization of the squared differences in static deflections. Then, the design of non-structural masses was determined via an optimization problem where the difference in mode shapes was minimized and the reduced modal frequencies were set as constraints.

In the field of modal optimization and stiffness matching, Pereira et al. [5] described a technique based on the optimization of the natural frequencies applied to the design of a joined-wing wind tunnel model. Bond et al. [6] showed that matching mode shapes is also a requirement, apart from natural frequencies. Richards et al. [7] compared a single-step direct modal matching to a decoupled approach. In the latter case, the stiffness is matched through static deflections and subsequently the modes and frequencies are matched by optimizing the nonstructural mass distribution.

Several researchers considered the nonlinearities in the static scaling, such as Bond et al. [6], who expanded the modern approach described by French and Eastep [4] to include geometrical nonlinearities in the aeroelastic scaling process. Ricciardi et al. [8] modified the two-step approach proposed by Richards et al. [7] to include the matching of nonlinear static deflections in the first optimization loop. Both of these efforts apply the proposed scaling methodology to the joined-wing SensorCraft model. Subsequently, Ricciardi et al. [9] investigated the use of a single-step method where linear and nonlinear static responses were matched, while satisfying the natural frequency constraints. Wan and Cesnik [10] presented a technique for the scaling of very flexible aircraft with geometrical nonlinearities. They extended the linear scaling factors and similarity rules to consider the case of aeroelastic scaling with geometrical nonlinearities.

More recently, Ricciardi et al. [11] presented a systematic approach for the design of aeroelastically scaled models. They used an optimization-based method to match vibration and buckling modes and eigenvalues, as well as to match a linear static response. They pointed out the problem of mode swapping and mixing during the optimization. Mas Colomer et al. [12, 13] described the static aeroelastic scaling when flow similarity cannot be achieved. In the context of the modern aeroelastic scaling (i.e., with assumed flow similarity), they addressed the problem of mode crossing by implementing a mode tracking strategy [14]. Cavallaro and Demasi [15] presented a detailed literature survey on aeroelastic scaling of the joined-wing SensorCraft models. Another review by Afonso et al. [16] focused on the state-of-the-art on nonlinear aeroelasticity of high aspect ratio wings. On the applications side, Spada et al. [17] applied the two-step method to the high aspect ratio wing of the reference aircraft [18] developed in the context of the NOVEMOR project under European Union’s 7th framework program. Pontillo et al. [19] presented a set of tools to design a cantilevered wing model that reproduces the aeroelastic static and dynamic behavior of a conceptual full size aircraft.

All the efforts mentioned above assume that there exists flow similarity between both scales, or that at least the differences in flow conditions (compressibility and viscous effects) are negligible. However, this may not be the case for flying demonstrators that have a limited flight envelope (both in altitude and speed). Therefore, a method is needed for the cases where a flying demonstrator is required to have the same flutter behavior (i.e., scaled flutter speeds, frequencies and modes) despite the differences in flow conditions,

In this paper, we first review the classical similarity criteria used for flight demonstrator scaling. Then, we present a new optimization-based approach that maximizes the similarity of the dynamic aeroelastic response when flow similarity cannot be achieved, focusing on differences in the Mach number. We focus in the Mach number only since the formulation of the aerodynamic model used — the doublet lattice method (DLM), included in Nastran — only considers corrections based on the Mach number.

To evaluate the proposed method, we apply it to the Common Research Model (CRM) wing [20], which is a well-known reference wing — corresponding to the in-flight shape — typically used to benchmark aerodynamic analysis methods. The goal is to optimize the wing planform to achieve the same flutter behavior in the scaled flight demonstrator despite the different flight Mach number.

With the application of the proposed method, we see that the error in the flutter speed, due to the Mach number difference, is effectively reduced. We also observe a significant improvement in the agreement in the evolution of the damping of the first flutter mode with airspeed through the V -g plot.

2 Methodology

The proposed method aims to match the aeroelastic response of two wings with different geometries. The method addresses both differences in external shape and internal structure geometry. To illustrate this, in Figure 1a we consider the geometry at rest of a reference wing (in blue) and another wing with a different chord (in red). Through the use of interpolation, we transfer the aeroelastic response of the reference wing to the new wing, as shown in Figure 1b. In this paper, we use optimization to match the two aeroelastic responses even if the geometries of the wing are not exactly coincident.

The aeroelastic equation of the reference aircraft (with uniform dimensions) is (1)

where M is the mass matrix, K is the structural stiffness matrix, A is the aerodynamic matrix, and x is the vectorof structural degrees of freedom (DOF). For simplicity, we assume that all matrices and vectors have uniform physical dimensions. The subscript r indicates the reference aircraft. We only use the aerodynamic matrix that multiplies the structural displacements, but the procedure would be equivalent ifother aerodynamic matrices (i.e., multiplying velocity and acceleration vectors) were considered. Also, if a harmonic response is considered, A represents a complex matrix that includes the effect of positions, velocities, and accelerations.

We write the same equation for the model (subscript m) whose aeroelastic response we want to match with the reference aircraft (2)

From that point, we derive the equations in a generic form, without distinguishing between the reference aircraft and the scaled model, as the procedure is the same for both cases. We consider that we can use the normal mode basis to express the structural displacements (3)

where η are the modal coordinates. Substituting equation (3) into equation (1) and left-multiplying by [Φ_r]^T, yields the generalized aeroelastic equation (4)

As described by Bisplinghoff [1] and further detailed by Ricciardi et al. [9], the aeroelastic response of two aircraft will be the same if their aeroelastic equations, once nondimensionalized, are coincident. Therefore, we nondimensionalize equation (4) and compare them to see which quantities have to be preserved to ensure the scaling of the aeroelastic response. We start by using the bi-orthogonality property of the modal problem to rewrite equation (4) as (5)

where represents the diagonal matrix of modal masses, and is the diagonal matrix of natural frequencies. By using the reference quantities for each scale (listed in Tab. 1), we write equation (5) in terms of these reference quantities and the nondimensional matrices as (6)

where the barred matrices are nondimensional and ^★ denotes the derivative with respect to nondimensional time. By reordering the multiplying factors comprising the reference quantities, we identify two nondimensional groups. These are the reduced frequency κ and the inertia ratio μ, which appear in the nondimensional equations, now instantiated for the reference aircraft r and the scaled model m, (7)

and (8)

Now, let us assume that the nondimensional modes of the model aircraft are such that they are equal to the modes of thereference structure when interpolated onto the model mesh, i.e., (9)

where [H] is the matrix that interpolates the structural displacements from the mesh of the reference aircraft to the mesh of the model aircraft. Substitutingequation (9) into equation (8), yields (10)

By comparing equation (10) to equation (7), we conclude that to achieve an equivalent response between the reference aircraft and the model, we have to ensure that (11) (12) (13) (14) (15)

and (16)

The information for the first three conditions can be obtained through a simple normal mode analysis, while the information necessary toevaluate the last condition can be obtained from an aeroelastic analysis module such as Nastran’s solver option SOL 145.

The aerodynamic shape and the flow similarity (Mach and Reynolds numbers) are required to guarantee the equality of the nondimensional aerodynamic matrices. Traditionally, as explained by Ricciardi et al. [9] and Pires [21], compressibility and viscous effects are neglected in the design of scaled flight demonstrators. For the cases where the internal structure architecture of the scaled demonstrator is different from that of the reference aircraft, the nondimensional modal masses, frequencies, and shapes are usually obtained by optimizing the scaled model parameters, namely structural thicknesses and nonstructural masses [9].

As previously mentioned, the nondimensional aerodynamic matrices depend on the aerodynamic shape and the flow conditions. In the case of airliners, the reference aircraft flies at high transonic Mach numbers (M ≈ 0.8), whereas the flight demonstrator flies at near-incompressible conditions (M ≈ 0.3).

Fig. 1 From the geometry at rest (a) of the reference wing (blue) and the wing under design (red) with a different geometry, a given aeroelastic response can be interpolated from the reference wing to the new one (b).

2.1 Wing planform optimization

We formulate an optimization problem that matches the dynamic aeroelastic response as closely as possible through the minimization of the squared error between the aerodynamic matrices of both models. The aeroelastic equations of motion are typically used to find the aircraft flutter modes, speeds, and frequencies. The flutter points correspond to the airspeed values for which the aeroelastic oscillations are undamped. Below the flutter airspeed, these oscillations are damped, whereas they are divergent above the flutter airspeed. In those cases, the equation of linear aeroelasticity is solved for harmonic oscillatory solutions. Speeds are expressed as {ẋ } = iω{x} and accelerations as {ẍ} = −ω²{x}. This allows writing the harmonic solution as a complex matrix that when multiplied by the displacements gives the aerodynamic forces due to displacements, speeds, and accelerations. This complex aerodynamic matrix is computed for each frequency ω, and depends on the Mach number, M. Taking this into consideration, and considering that we have a reference aircraft (subscript r) and the model we want to optimize to have the same aeroelastic behavior (subscript m), we write the equation of aeroelasticity for each aircraft as (17)

and (18)

where [Ā_h](X_a, κ, M) is the generalized aerodynamic matrix. This is a complex matrix that yields the aerodynamic forces due to displacements, speeds, and accelerations for the harmonic solution case when multiplied by the displacements. This complex matrix depends on the aerodynamic surface geometry (X_a), the reduced frequency (κ), and the Mach number (M).

As seen on the previous section, aeroelastic similarity between two aircraft requires each term in equation (18) to be equal to its counterpart in equation (17). The left-hand side is matched through modal optimization [14]. The right-hand side would be equal if flow similarity existed and the aerodynamic shape were preserved. If flow similarity cannot be achieved, then we attempt to find the model design parameters affecting the right-hand side of equation (18) that maximize the similitude with the corresponding term in equation (17). Throughout the rest of this paper, the term “flow similarity” applies to equality in Mach number only, since the aerodynamic model we use (DLM, included in Nastran) considers corrections based on the Mach number uniquely.

To maximize the flow similarity, we solve an optimization problem that determines the design variables for the model wing planform that minimize the difference between the two terms mentioned above. The objective function for this problem is (19)

This function quantifies the error between the two aerodynamic models through the sum of the squared L² norms of the difference between the aerodynamic matrices for a set of reduced frequencies (κ_i, ∀i ∈ 1, …, N).

2.2 Equivalent optimization problem formulation

Typically, the nondimensional generalized matrices [Ā_hr] and [Ā_hm] would be used to evaluate the objective function in equation (19). However, it may be difficult and error-prone from the user point of view to establish the dimensional transformation matrices for large Nastran models. Therefore, we use the generalized matrices, [Ă_h], which are direct outputs from MSC Nastran (in the Nastran documentation these are denoted as [Q_hh ]; see the Appendix for more details on how these matrices are computed.) Using these matrices, we define an optimization problem that is equivalent to the one expressed by the minimization of the objective function in equation (19).

If we were to compute [Ā_h], then [Ă_h ] would be nondimensionalized using a factor that is a function of the physical reference quantities of each aircraft (r or m). This process is the same for both the reference aircraft and the flight demonstrator, and therefore it is only described once for a generic matrix [Ă_h ]. Then, matrix [Ă_h] is nondimensionalized as (20)

Instead of computing [Ā_hr] and [Ā_hm], we use (21)

and (22)

both computed at the same physical scale — thus the same — to define a new optimization problem. Subsequently, we prove that the solution of this newly defined problem is the same as the solution of the problem defined at the end of Section 2.1. The advantage of this approach is that we use matrices that are a direct Nastran output, without having to uniformize (i.e., make them independent of the translation and rotation DOFs units, through the use of a reference length) and nondimensionalize them.

For the purpose of solving a different optimization problem that does not require direct computation of uniform and nondimensional matrices, and whose solution is the same as the minimum of the objective function on equation (19), we start by using the property that (23)

We start by writing a new objective function f_eq as (24)

Now, by substituting equation (20) into equation (24), we write f_eq as a function of the generalized matrices and , (25)

Substituting equation (19) into equation (25) yields (26)

Then, from equations (23) and (26), we get (27)

and therefore (28)

Now we can solve the same optimization problem by using the generalized matrices [Ă_hr] and [Ă_hm] built at the same scale and using the same normal modes [Φ]. The only difference between the two models is the planform shape and the Mach number.

2.3 Numerical tools

Now we describe the tools used to implement the strategy introduced above. For the aeroelastic analysis, we use MSC Nastran [22]. By using the SOL 145 solution sequence, we compute the flutter response of a wing given its structural finite element model, the definition of the aerodynamic surface planform, and a subset of the finite element model nodes to interpolate the displacements from the structural model to the aerodynamic model grid.

To generate the outer mold line (OML) wing geometry according to the design variables we use the GMSH [23] tool. Once the OML is generated for a particular design, the structural mesh is adapted accordingly using the mesh interpolation method proposed by Rendall and Allen [24].

The optimizer we use to solve this problem is COBYLA [25], which is a gradient-free, trust-region optimization method.

3 Example problem

In this section, we apply the wing planform optimization strategy we developed to an existing flutter model of the CRM wing [20]. In this application, we find the wing planform design that minimizes the error between the aerodynamic model of that wing at a given Mach number and the baseline CRM wing model at a different Mach number.

3.1 Description of the model

The model that we use as an example is based on a MSC Nastran model of the CRM wing¹. This model represents the wing of an airliner similar in dimensions and properties to the one of a Boeing 777 [26]. The wingbox is composed ofan upper and lower skin panels, leading and trailing edge spars, and ribs. Spar and rib caps are also included in the model. The masses of the wing parts on the leading and trailing edges that are not explicitly represented in the finite element model are represented by lumped masses attached to the spars using multipoint constraints (MPCs). Since the fuel mass inthe wing contributes significantly to the flutter behavior, the density of the structural material is artificially increased to account for the maximum fuel weight of a Boeing 777 — similar in dimensions to the CRM wing.

The wing aerodynamic model consists of MSC Nastran CAERO1 panels. This model computes the aerodynamic forces using the doublet-lattice method (DLM). Both structural and aerodynamic models are coupled within Nastran through surface splines that interpolate the forces and displacements [22].

Figure 2 depicts the wingbox structural model, which extends only from the leading edge spar to the trailing edge spar.

Figure 3 shows the dimensions of the baseline CRM wing. We use this baseline design as the reference wing flying at thereference Mach number for the optimization problem in the next section.

Fig. 2 The finite element model of the CRM wingbox consists of 12475 shell elements, 14134 bar elements, and 10804 nodes.

Fig. 3 Wing planform design variables: wing root, break, and tip chords (cr, cb, and ct respectively)and sweep angle Λ.

3.2 Optimization problem

In this section, we apply the optimization strategy described in Sections 2.1 and 2.2 to the model presented in Section 3.1. In this case, the baseline CRM wing configuration, shown in Figure 3, is taken as the reference aircraft, and its generalized aerodynamic matrix is computed at M = 0.8. The design variables in the optimization are the chord length and the leading edge position of both the root section and the wing tip, for a total of four design variables, as shown in Figure 3. The Mach number for this wing is M = 0.3. Thus, the wing OML changes between the two aircraft but the internal structure remains the same.

The optimization problem consists in minimizing the objective function defined by equation (24) without constraints. In that case, we consider the objective function for one reduced frequency only, κ₁ = 0.183, which is theone corresponding to the first flutter mode of the baseline CRM wing at the reference conditions. Table 2 summarizes the objective function, design variables, and the problem parameters.

4 Results

In this section, we first present and discuss the optimization results. Then, we analyze the effects of the change in the Mach number on the baseline wing flutter behavior and analyze how the optimized wing for the new Mach number reduces the error in flutter behavior.

The objective function history is plotted against the number of iterations in Figure 4a for the solution of the problem detailed in Section 3.2. The objective function value of the baseline wing is 2.7 ×10⁻⁵ m⁶, and the optimal value after COBYLA performs 50 iterations is 1.1 × 10⁻⁵ m⁶, which represents a substantial reduction.

The optimal wing shows an increase in both the root and tip chords, but the increment is larger at the root, as shown in Figure 4b. The leading edge of both sections moves forward. The increase of both chords results in a wing area increment, which is consistent with the fact that the flight demonstrator Mach number M_m is lower, since a decrease in Mach number (in the subsonic regime) decreases the aerodynamic forces. Therefore, the optimizer increases the area to produce equivalent aerodynamic forces at a lower Mach number. In Figure 5a we see how the geometry of the wingbox is updated according to the optimized planform, while keeping the same layout as the reference wingbox in Figure 5b.

Since our goal is to reduce the error in the flutter response of a wing when the Mach number changes, wefirst compare the baseline wing aeroelastic response at the reference conditions (M_r) to the baseline wing at the model Mach number (M_m). Then, we do the same comparison between the reference case and the optimized wing at M_m by computing the V -g plot for each case (Fig. 6), which allows us to determine the flutter speed for each case.

From the flutter responses shown in Figure 6, we see an improvement with the optimized wing with respect to the baseline wing on the flutter speed. From the flutter response, we quantify the error in flutter speed both in absolute and relative terms for the baseline and optimized wing with respect to the same reference case (Tab. 3).

In Table 3, we see that the error in the flutter speed when using the baseline wing is 7.79%, while it decreases to 2.13% when using the optimized wing for the model Mach number. Thus, we see a significant improvement on the error in the flutter speed when considering the optimized design with the compressibility conditions of the scaled model.

Fig. 4 Evolution of the objective function (24) with the number of optimization iterations (a). Optimized planform (red) and baseline (blue) (b.

Fig. 5 The geometry of the optimized wingbox is updated according to the wing planform (a) while keeping the same layout as the reference wing (b).

Fig. 6 Baseline CRM wing V -g plots at M = 0.8 and for the optimized planform at M = 0.3; the optimized planform at model conditions brings the flutter speed closer to the one of the reference wing at reference conditions.

5 Validation of the interpolated modes

Since the interpolation matrix H for the transfer of modal shapes depends on the coordinates of both the structural mesh of the reference configuration and the ones of the structural mesh of the optimized model, we perform a qualitative validation of the interpolated shapes on the final design.We do this by representing both the reference modes on the baseline mesh and their interpolated counterparts at once for each modal shape to verify qualitatively that the nature of each mode is well preserved on the optimized design. In Figure 7 we see that the bending and torsion displacements of the first five reference modes considered (in blue) are well reproduced on the mesh of the optimized planform once interpolated (in red).

Fig. 7 First five natural vibration modes of the structure of the reference structure (blue) and their interpolation onto the structure mesh of the optimized planform (red).

6 Conclusions

Flight testing aeroelastically scaled models is a way of obtaining experimental knowledge on the in-flight behavior of innovative aircraft concepts. Traditional aeroelastic scaling of flying models considers that flow similarity exists or, at least, that the flow differences are negligible. However, that is not the case of airliners that fly in the transonic regime.

In this paper, we presented a method that maximizes the similarity in the dynamic aeroelastic response between two wings — despite different Mach numbers — by optimizing the planform of one of them.

Using the CRM wing as a test case, we first evaluated the error when the same wing is used at the flight demonstrator conditions (M_m = 0.3) to reproduce the reference wing aeroelastic behavior at the reference aircraft Mach number (M_r = 0.8). Then, we applied the proposed optimization method to that case, with the baseline design as a starting point. The proposed approach reduced the error in flutter speedfrom 7.79% to 2.13%.

This method allows the design of aeroelastically scaled models when the differences in the compressibility conditions cannot be neglected. The greater the change in these conditions, the greater the changes in the planform, so this method should only be considered when the resultant planform is still representative of the reference aircraft concept. Given the formulation of the DLM used here, this method only holds for the dynamic aeroelastic response expressed by the flutter model.

The approach proposed in this paper focuses on the similarity of the aerodynamic matrices used for flutter analyses only. It assumes that similarity of the normal modes and scaled natural frequencies is achieved by a separated optimization of the structural and mass properties to try to satisfy equations (11)–(13). This structural optimization should be performed after the optimal planform is determined by the method proposed in the current approach.

Appendix

A.1 Aerodynamic matrices in Nastran

In this section, we describe how the aerodynamic matrices [ă_h], required to evaluate the objective function in equation (24), are computed by Nastran using the SOL 145 solution sequence. We first describe all the necessary elements to its construction and then we assemble them to build [ă_h].

According to the Nastran Aeroelastic Analysis User’s Guide [22], the dimensionless vertical (i.e., normal to the panel) velocity is (A.1)

where [A_jj] is the aerodynamic influence coefficient matrix, p_j is the lifting pressure on element j, and q is the flight dynamic pressure. In this section, subindices j and k refer to the aerodynamic panels and grid points respectively. For the dynamic case, we write the dimensionless normal velocity (or downwash) as (A.2)

where and are the downwash matrices, κ is the reduced frequency, and {u_k} is the vector of displacements at aerodynamic grid points. By integrating the pressure at the panels, {p_j }, we get the forces on the aerodynamic grid points, (A.3)

where [S_kj] is the integration matrix, which depends on the surface of the panels.

By combining equations (A.1)–(A.3), we write the forces on the aerodynamic grid points as a function of the aerodynamic grid deflection as (A.4)

Now we express the displacements of the aerodynamic grid points as a function of the structural displacements {x_a } as (A.5)

where [G_ka] is the displacement interpolation matrix. By using the same interpolation matrix we get the forces on the structure points as (A.6)

By using equations (A.4)–(A.6), we write the aerodynamic forces on the structure points as a function of the displacements on the structure points as (A.7)

Nastran then computes the generalized aerodynamic matrix, (A.8)

If we were to compute [ă_h], we would uniformize and generalize (in the sense of the normal modes) the aerodynamic matrix [Q_aa] as (A.9)

From equations (A.8) and (A.9), knowing that and , where is the identity matrix, we conclude that (A.10)

Indeed, [ă_h] can be obtained as the matrix [Q_hh] computed by Nastran, using the SOL 145 solution sequence along with the CAERO1 cards for the definition of the aerodynamic panels.

References

All Tables

All Figures

	Fig. 1 From the geometry at rest (a) of the reference wing (blue) and the wing under design (red) with a different geometry, a given aeroelastic response can be interpolated from the reference wing to the new one (b).
In the text

	Fig. 2 The finite element model of the CRM wingbox consists of 12475 shell elements, 14134 bar elements, and 10804 nodes.
In the text

	Fig. 3 Wing planform design variables: wing root, break, and tip chords (cr, cb, and ct respectively)and sweep angle Λ.
In the text

	Fig. 4 Evolution of the objective function (24) with the number of optimization iterations (a). Optimized planform (red) and baseline (blue) (b.
In the text

	Fig. 5 The geometry of the optimized wingbox is updated according to the wing planform (a) while keeping the same layout as the reference wing (b).
In the text

	Fig. 6 Baseline CRM wing V -g plots at M = 0.8 and for the optimized planform at M = 0.3; the optimized planform at model conditions brings the flutter speed closer to the one of the reference wing at reference conditions.
In the text

	Fig. 7 First five natural vibration modes of the structure of the reference structure (blue) and their interpolation onto the structure mesh of the optimized planform (red).
In the text