Aero-Propulsive Modeling for eVTOL Aircraft Using Wind Tunnel Testing with Multisine Inputs

I. Introduction ∗ any novel aircraft designs are currently under investigation to support a future Urban Air Mobility (UAM) M transportation system [ 1 , 2 ]. Vehicles supporting UAM operations require precise hover and efficient cruise capabilities as well as the ability to safely transition between flight regimes. Hybrid and electric vertical takeoff and landing (eVTOL) aircraft are a clear fit for this new transportation model. In general, eVTOL aircraft are a combination of traditional fixed-wing and rotary-wing aircraft leveraging certain attributes from each type of vehicle. Rotary-wing aircraft features provide the ability to takeoff and land vertically, hover, and precisely maneuver in confined areas. Longer endurance, better efficiency, and the ability to operate at high speeds is derived from fixed-wing aircraft attributes.

Distributed hybrid and electric propulsion technology used in many eVTOL aircraft concepts has further broadened the traditional aeronautical vehicle design space and resulted in many unique vehicle designs [3–9].

Although the operational utility of eVTOL aircraft has great potential, there are many research areas to be addressed prior to introduction into a UAM environment. One essential research area is accurate eVTOL vehicle aero-propulsive modeling enabling flight dynamics simulation development to support flight control system design and certification, research on practical flight operations, and many other areas. Efficient and accurate vehicle aero-propulsive model development, however, is challenged by several vehicle attributes, including: many control surfaces and propulsors, propulsion-airframe interactions, propulsor-propulsor interactions, high-incidence angle propulsor aerodynamics, vehicle instability, rapidly changing aerodynamics through transition, and large flight envelopes that need to be characterized by a global aero-propulsive model. Compared to many conventional aircraft, eVTOL aircraft designs exhibit greater aero-propulsive complexity and many interacting factors requiring development of new testing and model development strategies. For conventional fixed-wing aircraft, the propulsion aerodynamics and airframe aerodynamics can generally be studied separately because the interaction effects are relatively small. In contrast, integrated aero-propulsive modeling is performed for eVTOL aircraft because the propulsion aerodynamics, airframe aerodynamics, and propulsion-airframe interactional aerodynamics are generally highly coupled. Previous research has investigated methods for efficient eVTOL aircraft aero-propulsive modeling across their wide flight envelopes using computational fluid dynamics simulations [10, 11], wind tunnel testing [12–16], and simulated flight testing [17].

The present work builds on previous eVTOL aircraft modeling research to develop an effective aero-propulsive modeling strategy for eVTOL aircraft using novel, efficient wind tunnel testing. Consequently, the following discussion focuses on methods used to develop mathematical models from experimental data, termed aircraft system identification [ 18 – 20 ]. The distributed propulsion, tilt-wing aircraft modeling approach proposed in Ref. [ 16 ] is combined with a more efficient, hybrid wind tunnel experiment designed using modern response surface methods [ 21 , 22 ] and dynamic programmed test input (PTI) excitations [ 23 – 25 ]. This multilayered testing and modeling strategy is similar to previous fixed-wing aircraft modeling research using flight testing [ 26 ] and flight simulations [ 27 ]. The present work is also inspired by recent NASA Learn-to-Fly Project [ 28 , 29 ] wind tunnel testing [ 30 ] and flight testing [ 31 – 33 ] used for modeling fixed-wing aircraft and extends efficient model development strategies to eVTOL aircraft. Reference [ 30 ] presents a related wind tunnel testing and modeling approach leveraging a combination of response surface designs and PTI excitations used for fixed-wing aircraft aerodynamic modeling. The new contributions of this work include a new experimental design methodology and a novel, multistep model identification procedure tailored to the data collection methods, as well as application to a complex eVTOL aircraft with significant aero-propulsive coupling.

The paper is organized as follows: Section II introduces the experimental aircraft. Section III describes the experimental design methodology, followed by a description of the wind tunnel test and data processing methods in Sec. IV. The aero-propulsive modeling approach and the model identification strategy are described in Sec. V.

Section VI provides sample modeling results, followed by discussion of the results and modeling approach in Sec. VII.

Overall conclusions are summarized in Sec. VIII.

II. Aircraft The modeling approach developed in this paper is applied to the Langley Aerodrome No. 8 (LA-8) vehicle [ 9 ].

The LA-8, pictured in Fig. 1, is a subscale, tandem tilt-wing, distributed electric propulsion aircraft configuration built as a testbed for eVTOL aircraft technology. The LA-8 was developed at NASA Langley Research Center as one of several eVTOL research aircraft intended to explore their unique flight characteristics and resolve implementation challenges to help bring similar full-scale vehicles into mainstream operation. The LA-8 project has enabled research in rapid vehicle development [ 9 , 34 ], computational predictions [ 35 ], wind tunnel testing [ 15 , 36 , 37 ], aero-propulsive modeling [16, 17, 38], flight controls [39], and flight test strategies [17, 40].

∗ Information available online at https://evtol.news/aircraft [accessed 29 October 2021]

(a) LA-8 front view (b) LA-8 rear view Fig. 1 LA-8 mounted in the NASA Langley 12-Foot Low-Speed Tunnel.

The LA-8 is equipped with 20 control effectors, including two tilting wings, four elevons, four flaps, two ruddervators, and eight electric motors/propellers. Figure 2 is a photo of the LA-8 with annotations showing the vehicle propulsor and control surface definitions. Wing, elevon, flap, and ruddervator deflections are defined positive trailing edge downward.

Propellers 1, 3, 6, and 8 rotate clockwise and propellers 2, 4, 5, and 7 rotate counterclockwise, as viewed from the rear.

All propellers are 16-inch diameter, 8-inch pitch, fixed-pitch, 3-bladed propellers.

𝒏 𝒏 𝟐 𝟑 𝒏 𝟏 𝒏 front 𝟒 𝜹 𝒘 𝟏 𝜹 𝜹 𝒙 𝒆 𝜹 𝒃 𝒇 𝟏 𝜹 𝟏 𝒇 𝒆 𝟐 𝟐 𝒏 𝒏 𝟔 𝟕 𝒚 𝒏 𝒃 𝟓 𝒏 𝟖 𝜹 𝒘 𝟐 𝜹 𝒆 𝜹 𝜹 𝜹 𝟑 𝒆 𝒇 𝒇 𝟒 𝟒 𝟑 𝜹 𝜹 𝒓 𝒓 𝟏 𝟐 Fig. 2 LA-8 propulsor and control surface definitions.

III. Experiment Design The factors under test were divided into static factors and dynamic factors for the wind tunnel experiment design based on the frequency band of typical variations in flight, as well as the capabilities of the wind tunnel testing apparatus.

Static and dynamic test factor experiment designs were developed separately, but the designs were run simultaneously for test execution. The respective experiment design methods forming the hybrid testing strategy are discussed in the following subsections. Note that this study was exploratory in nature with an objective of refining the proposed testing approach to guide future wind tunnel tests for eVTOL vehicles. Consequently, the wind tunnel experiment was expected to yield excess data with the goal of determining the data collection requirements enabling development of satisfactory aero-propulsive models in an efficient manner.

A. Static Experiment Design The static test factors were angle of attack 𝛼 , angle of sideslip 𝛽 , front wing angle 𝛿 , and rear wing angle 𝛿 .

𝑤 𝑤 1 2 These factors were held at constant settings during data collection due to their slow relative movement and for operational convenience in wind tunnel testing. Moving the airflow orientation angles and wing angles dynamically changes the tare values and introduces additional dynamic modeling complexity. The factor settings for each of the four static test factors were independently commanded by the wind tunnel control system at each test point. Although previous research suggests parameterizing eVTOL aircraft models using body-axis velocity components [ 11 , 16 ], the test matrices were specified using 𝛼 and 𝛽 for ease of envelope definition and simplified integration into the wind tunnel test apparatus.

The freestream velocity was held at a constant setting for the testing described in this paper, but would need to be varied to develop an aero-propulsive model valid throughout the operational flight envelope.

A set of five sequential static test blocks was developed using design of experiments (DOE) [ 21 ] and response surface methodology (RSM) [ 22 ] theory to acquire the data necessary to identify and validate increasingly complex ® aero-propulsive models. Modeling block design was accomplished with the aid of Design-Expert , a commercially † available statistical software package. Previous model identification results for the LA-8 aircraft suggested that at least pure quadratic and two-factor interaction model terms are needed to characterize the slow moving test factors [ 16 ], but additional model complexity in the form of higher-order model terms was not investigated in that study. The hybrid testing strategy employed in this work allows identification of higher-order models, while still conducting efficient testing.

The static factor settings were chosen using sequential, completely randomized I-optimal response surface designs, which are designs that minimize the average prediction variance for a predefined model order over the range of factor settings [21, 22]. The four blocks used for modeling were I-optimal designs for: 1) a quadratic design model (up to pure quadratic and two-factor interaction model terms), 2) a cubic design model (up to pure cubic terms and all arrangements of cross terms up to a total of third order), 3) a quartic design model (up to pure quartic terms and all arrangements of cross terms up to a total of fourth order), and 4) a quintic design model (up to pure quintic terms and all arrangements of cross terms up to a total of fifth order).

The I-optimal block designs were composed of the minimum number of test points needed to fit a full model of the specified complexity. Three additional center points were also included in each block to allow estimation of pure error and to aid in stabilization of the prediction variance at the center of the experimental region [ 21 ]. Each sequential I-optimal design augmented the collective design from previous blocks to sequentially increase the identifiable model complexity, improve the model prediction, and avoid duplicating previously tested combinations of factor settings. All modeling blocks were executed to investigate the model complexity needed to characterize eVTOL aircraft, but not all blocks were expected to be needed to identify a satisfactory aero-propulsive model.

An additional fifth block consisting of 20 static test points selected using a random number generator was used as validation data withheld from model identification. The choice of randomized test points makes the validation test impartial to the experimental design. This validation block was found to provide a good estimate of prediction error while using a modest number of test points for the present application.

Figure 3 shows two-dimensional slices of the factor space for the four static test factors. Each block is plotted sequentially with points from the previous blocks to show how the higher complexity I-optimal designs fill the design space. The randomized validation test points are also shown on each plot. Figure 4 shows a four-dimensional representation for all static test factors and all test blocks. Although the factors are displayed in coded units in the figure, the variables would be converted into engineering units for test execution.

Assessment of the prediction variance of a response surface design for a given model structure provides insight into its precision of prediction and allows comparison of different response surface designs. The variance of the predicted response is   − 1 2 𝑇 𝑇 Var [ ˆ 𝑦 ( 𝒙 )] = 𝜎 𝒙 𝑿 𝑿 𝒙 (1) 0 0 where ˆ 𝑦 ( 𝒙 ) is the predicted response evaluated at the design space location 𝒙 expanded to the form of the model 0 0 structure, 𝑿 is a matrix composed of the designed test points in the form of the model structure, and 𝜎 is the measurement error variance [ 22 ]. From Eq. (1) , the prediction variance is a function of the experiment design, the model structure, the location in the design space, and the measurement facility error variance. The scaled prediction variance (SPV) and unscaled prediction variance (UPV) are used to compare experiment designs prior to experimentation, which removes the dependence on 𝜎 . SPV is defined as   − 1 𝑁 Var [ ˆ 𝑦 ( 𝒙 )] 𝑇 𝑇 SPV = = 𝑁 𝒙 𝑿 𝑿 𝒙 (2) 𝜎 where the number of test points 𝑁 penalizes a larger design size [ 22 ]. The SPV considers the prediction accuracy as † Information available online at https://www.statease.com/software/design-expert/ [accessed 29 October 2021] Fig. 3 Sequential two-dimensional slices of the coded factor space for the static test factors.

Fig. 4 Four-dimensional representation of the coded factor space for the static test factors.

well as the expense of test points when comparing designs. The UPV, defined as   − 1 Var [ ˆ 𝑦 ( 𝒙 )] 𝑇 𝑇 UPV = = 𝒙 𝑿 𝑿 𝒙 (3) 𝜎 provides an assessment of the prediction precision independent from the size of the experimental design.

Graphical presentation of the distribution of prediction variance throughout the design space is an effective way to assess experimental designs. Each test block design used for this work was assessed using fraction of design space (FDS) plots [ 41 ]. FDS graphs depict the prediction variance distribution over the design space in a concise manner, where the prediction variance metrics are plotted against the FDS encompassing a prediction variance less than or equal to a particular value. It is also useful to consider the FDS including a particular model precision, quantified by the confidence interval half-width 𝛿 [ 42 – 44 ]. The model precision 𝛿 normalized by the response standard deviation 𝜎 plotted against FDS provides further insight into the prediction capability of the model developed from a particular experiment design, prior to conducting the experiment. For this study, a design was deemed to be adequate for fitting a particular model complexity if 𝛿 / 𝜎 was less than two for greater than 95% of the design space. The prediction variance ∗ threshold PV used to determine the FDS within a given model precision level is   𝛿 / 𝜎 ∗ PV = (4) 𝑡 𝛼 / 2 , 𝑁 − 𝑝 where 𝑁 is the number of test points, 𝑝 is the number of parameters in the model, and 𝛼 is the significance level (chosen as 𝛼 = 0 . 05 for this study).

Table 1 lists the number of test points for each block, as well as the FDS with 𝛿 / 𝜎 ≤ 2 for quadratic, cubic, and quartic evaluation models. An adequate FDS ( FDS ≥ 0 . 95 ) for a normalized model precision 𝛿 / 𝜎 ≤ 2 is obtained with a design order one power larger than the evaluation model order. Figures 5, 6, and 7 show the UPV, SPV, and 𝛿 / 𝜎 threshold values against FDS for a quadratic, cubic, and quartic model evaluation model order, respectively. Each block is analyzed sequentially, meaning that the analysis for each block includes the design points in the current block and all previous blocks. The UPV and 𝛿 / 𝜎 threshold curve decreases in value and becomes more uniform (flat) as the design order increases. The SPV curves for each block are similar for the quadratic evaluation model. For the cubic and quartic evaluation models, a lower and more uniform SPV is obtained with a design complexity greater than the evaluation model complexity.

Table 1 Cumulative experimental design properties for each test block Design Block Cumulative Quadratic Model Cubic Model Quartic Model Block Order Points Points FDS with 𝛿 / 𝜎 ≤ 2 FDS with 𝛿 / 𝜎 ≤ 2 FDS with 𝛿 / 𝜎 ≤ 2 1 Quadratic 18 18 0.448 0.000 0.000 2 Cubic 23 41 0.999 0.215 0.000 3 Quartic 38 79 1.000 0.998 0.151 4 Quintic 59 138 1.000 1.000 0.994 Fig. 5 FDS plots for a quadratic evaluation model.

Fig. 6 FDS plots for a cubic evaluation model.

Fig. 7 FDS plots for a quartic evaluation model.

B. Dynamic Experiment Design The ten control surfaces (four elevons, four flaps, and two ruddervators) and eight propulsors present on the LA-8 aircraft were treated as dynamically changing factors at each static test point. Orthogonal phase-optimized multisine inputs [ 18 , 23 – 25 ] were the excitation input type, or PTI, used for these 18 individual control effectors. Multisine inputs have been used successfully in previous aircraft system identification work to simultaneously characterize the effects of many individual control effectors [45–48].

A multisine input is defined as a sum of multiple sinusoidal functions with different amplitudes, frequencies, and phase angles, where the frequencies are chosen to encompass the frequency range corresponding to the system dynamics of interest. To make all inputs orthogonal in both the time domain and frequency domain, the multisine signal for each 𝑗 th control effector is assigned sinusoids with a unique subset of discrete harmonic frequency indices 𝐾 selected from 𝑗 the complete set of available frequency indices, 𝐾 . The available frequencies are 𝑓 = 𝑘 / 𝑇, 𝑘 = 1 , 2 , ..., 𝐾 , where 𝑇 is 𝑘 the fundamental period. For 𝑚 total control effectors, each 𝑗 th input signal 𝒖 is defined as 𝑗   ∑︁ √︁ 2 𝜋𝑘 𝒕 𝒖 = 𝐴 𝑃 sin + 𝜙 𝑗 = 1 , 2 , ..., 𝑚 (5) 𝑗 𝑘 𝑘 𝑇 𝑘 ∈ 𝐾 𝑗 where 𝐴 is the signal amplitude, 𝑃 is the 𝑘 th power fraction, 𝜙 is the 𝑘 th phase angle defined on the interval (− 𝜋, + 𝜋 ] , 𝑘 𝑘 and 𝒕 is the time vector containing 𝑁 discrete points. The relative peak factor (RPF), defined as 

 max 𝒖 − min 𝒖 / 2
1 𝑗 𝑗 RPF 𝒖 = (6) √ √︃ 𝑗 𝑇 𝒖 𝒖 / 𝑁 𝑗 𝑗 is the range of input amplitude divided by the root-mean-square of the signal, referenced to the peak factor for a single-frequency sinusoid. Minimum RPF is preferred for system identification flight testing to prevent perturbing the aircraft far from the reference flight condition. For wind tunnel testing, minimum RPF also approximately minimizes the peak-to-peak amplitude of the multisine signals, which keeps the rate of change low. This is important for repeated dynamic use of the control surface actuators and electric motors during wind tunnel testing. The relative peak factor for a multisine signal is minimized by optimizing the phase angles using the simplex algorithm because the optimization problem is not convex [ 23 ]. Minimizing the RPF by phase-shifting the sinusoidal components of a multisine input does not affect the excitation input energy or the mutual orthogonality of the inputs.

For this study, individual multisine signals were generated for each of the LA-8 control surfaces and propulsors, for a total of 18 different excitation signals. The multisine design process followed similar steps to those shown in Ref. [ 17 ] for simulated LA-8 flight experiments. The difference in the design process for this current work is a consideration placed on identification of pure quadratic and interaction model terms among control effectors, in addition to the linear model terms considered in Ref. [ 17 ]. Several harmonic components were assigned to each control surface and propulsor multisine signal, where the overall frequency range was set to between 0.05 Hz and 1.8 Hz in accordance with frequencies that would be used in flight. The propulsion harmonic components were focused into lower frequencies below 1.2 Hz to adhere to the lower bandwidth of the propulsors [ 38 ]. Focusing the propulsor excitation frequencies below this limit avoids over-stressing the motors and prevents the propulsor excitation power from being attenuated by the motor dynamics. Six different sets of multisine input frequency components were considered for the experiment by varying the fundamental period 𝑇 from 30 seconds to 180 seconds in 30-second increments. A larger fundamental period results in a finer frequency resolution, Δ 𝑓 = 1 / 𝑇 Hz, which allows assigning more frequency components to each individual multisine signal. The design with the shortest fundamental period ( 𝑇 = 30 seconds) had three harmonic components assigned to each propulsor and control surface signal; the design with the longest fundamental period ( 𝑇 = 180 seconds) had 16 harmonic components assigned to each propulsor signal and 18 harmonic components assigned to each control surface signal. Because the starting phase angles for each harmonic component in the non-convex RPF optimization are generally chosen randomly in (− 𝜋, + 𝜋 ] , a different set of phase angles optimized for minimum RPF is generally obtained each time a multisine signal is designed. Multisine optimization with randomly chosen starting phase angles was performed 30 times for each different set of frequency components and the design with the shortest time to decrease the maximum absolute pairwise correlation among up to quadratic and two-factor interaction control effector model terms was selected to compare to signals developed with different fundamental periods.

Multisine signals are orthogonal in the time-domain at integer multiples of 𝑇 , which might be interpreted to suggest using a multisine design with a small fundamental period. However, obtaining high quality modeling results requires low correlation rather than zero correlation [ 25 ], meaning that good modeling results can be obtained by using a data collection time shorter than the fundamental period of the multisine signal. Also, using a larger number of frequency components provides more diverse dynamic information, which has been shown to improve modeling results [ 45 ].

Following an approach similar to Refs. [ 17 , 25 ], correlation metrics were used as criteria to assess the quality of each multisine design as data collection time progresses.

Correlation between two signals can be assessed using the pairwise correlation coefficient. The correlation coefficient 𝑟 between two signals, 𝒙 and 𝒙 , is defined as 𝑖 𝑗 𝑖 𝑗 𝑇 ( 𝒙 − ¯ 𝑥 ) ( 𝒙 − ¯ 𝑥 ) 𝑖 𝑖 𝑗 𝑗 𝑟 = (7) √︁ 𝑖 𝑗 √︁ 𝑇 𝑇 ( 𝒙 − ¯ 𝑥 ) ( 𝒙 − ¯ 𝑥 ) ( 𝒙 − ¯ 𝑥 ) ( 𝒙 − ¯ 𝑥 ) 𝑖 𝑖 𝑖 𝑖 𝑗 𝑗 𝑗 𝑗 where ¯ 𝑥 and ¯ 𝑥 are the respective mean signal values. A correlation coefficient value of zero means the signals are 𝑖 𝑗 uncorrelated, or orthogonal, and an absolute correlation coefficient of one indicates that the signals are completely correlated, or linearly dependent. Lower absolute correlation generally results in better modeling results [ 25 ]. A correlation coefficient between signals used for modeling with magnitude greater than 0.9 indicates that data collinearity, or correlation between signals high enough to cause corrupted model identification, may be encountered [ 18 , 20 ].

Another metric that can be used to assess signal correlation is the variance inflation factor (VIF). For the signal 𝒙 , the 𝑗 respective VIF is VIF = (8) 𝑗 1 − 𝑅 𝑗 where 𝑅 is the coefficient of determination obtained through creating a regression model of 𝒙 as a function of all other 𝑗 𝑗 regressor signals. A VIF value greater than 10 suggests that data collinearity may be present [ 18 , 42 , 49 ]. The 𝑟 , VIF , 𝑖 𝑗 𝑗 and 𝑅 metrics only quantify correlation between pairs of signals and, thus, cannot diagnose collinearity among more 𝑗 than two signals [18, 50].

An alternative method that can be used to assess multiple correlation between more than two inputs is analysis 𝑇 of the eigenvalues of 𝑿 𝑿 , where 𝑿 is a matrix composed of column vectors of the signals in the regression model.

𝑇 The inverse of the 𝑿 𝑿 matrix is required to compute the ordinary least-squares regression solution. The ratio of the 𝑇 maximum eigenvalue and minimum eigenvalue, 𝜅 = 𝜆 / 𝜆 , is the condition number of the 𝑿 𝑿 matrix. A value max min of 𝜅 close to one indicates low multiple correlation whereas a large value of 𝜅 indicates an ill-conditioned estimation problem due to data collinearity. Values of 𝜅 indicating adverse effects from data collinearity range anywhere from 100 to 100,000 depending on the particular data set [18, 20, 49, 50].

The evolution of correlation metrics over time for four multisine designs with a different fundamental period is shown in Fig. 8. A linear and full quadratic model (every possible linear, quadratic, and two-factor interaction model (a) Maximum absolute 𝒓 (linear model) (b) Maximum absolute 𝒓 (quadratic model) 𝒊 𝒋 𝒊 𝒋 (c) Maximum VIF (linear model) (d) Maximum VIF (quadratic model) 𝒋 𝒋 2 2 (e) Maximum 𝑹 (linear model) (f) Maximum 𝑹 (quadratic model) 𝒋 𝒋 𝑻 𝑻 (g) Condition number of 𝑿 𝑿 (linear model) (h) Condition number of 𝑿 𝑿 (quadratic model) Fig. 8 Input signal correlation metrics against time.

term) for all the control surfaces and propulsors are considered for this analysis. This represents a conservative analysis performed prior to the experiment that assumes all candidate regressors are included in the model. Many candidate model terms considered here, particularly for the full quadratic model, were expected to be excluded through model structure determination after data collection. The PTI excitations were run in a continuous loop for testing, meaning that different set points contained different portions of the full PTI signals. Consequently, the analysis here shows the highest value of each metric for time segments evaluated throughout the full duration of the periodic signal, as opposed to just from the start of the signal at 𝑡 = 0 .

Figures 8a and 8b show the maximum absolute pairwise correlation value, max (| 𝑟 |) , for each multisine design as a 𝑖 𝑗 function of time for a linear and quadratic evaluation model, respectively. The maximum correlation for each different multisine design shows a similar rate to decrease below 0.5, which is achieved in roughly 10 seconds for the linear model and roughly 25 seconds for the quadratic model. Figure 8a shows that the pairwise correlation values are zero at integer multiples of the fundamental period of the respective multisine design for a linear model. For the quadratic model, the correlation pairwise correlation values decrease until passing the fundamental period for the multisine design where the pairwise correlation value remains relatively constant thereafter. Multisine designs with a larger fundamental period ultimately achieve lower correlation among quadratic model regressors. Similar character is observed for the maximum VIF and 𝑅 variation with time shown in Figures 8c-8f. For a linear model, all multisine designs achieve a maximum 𝑗 𝑗 VIF of less than 10 in roughly 7 seconds. For a quadratic model, the multisine designs with a longer fundamental period achieve a maximum VIF of less than 10 in roughly 40 seconds.

𝑇 Figures 8g-8h show the condition number 𝜅 of 𝑿 𝑿 for each multisine design as a function of time. For the linear model, all multisine designs achieve a condition number less than 100 in approximately 7 seconds at nearly the same rate and, thereafter, the condition numbers continue to decrease at a slower rate. For the quadratic model, the condition number for each multisine design is seen to decrease until reaching the fundamental period, and a nearly constant condition number is observed thereafter. The multisine designs with longer fundamental periods achieve a condition number less than 1,000 in roughly 40 seconds. Again, this analysis assumes that all candidate regressors are included in the model. If a subset of the model terms is selected for the model, then the condition number is typically reduced.

The overall takeaway from Fig. 8 is that multisine designs with different fundamental periods reduce correlation metric values at a similar rate up to the respective fundamental period of each design. Also, multisine designs with longer fundamental periods obtain lower correlation metrics over time for a quadratic model. Informed by these time-dependent correlation analysis results and previous research showing the benefits of increased frequency resolution [ 45 ], the multisine design with the largest fundamental period ( 𝑇 = 180 seconds) was selected to be used for the wind tunnel experiment.

The input spectra for the final set of orthogonal phase-optimized multisine signals with a fundamental period of 𝑇 = 180 seconds is shown in Fig. 9. There are 308 total harmonic components, with 16 frequencies assigned to each propulsor and 18 frequencies assigned to each control surface in an alternating manner. The overall frequency range is between 𝑓 = 0 . 05 Hz and 𝑓 = 1 . 756 Hz with a frequency resolution of Δ 𝑓 = 1 / 𝑇 = 0 . 00556 Hz. The input min max spectra plot shows that the propulsor harmonic components are in a lower frequency range, reflecting that the input excitations were designed to be within the bandwidth of the propulsion system. Figure 10 shows the first 20 seconds of the input excitation signals normalized to have a maximum absolute value of one. The RPF values for the propulsor inputs were below 1.32 and the RPF values for the control surface inputs were below 1.60. A gain is applied to scale each input signal to a sufficient amplitude to obtain a good signal-to-noise ratio for model identification.

Fig. 9 Multisine input spectra for the LA-8 control effectors.

Fig. 10 Normalized multisine inputs for the LA-8 control effectors.

IV. Data Collection and Signal Processing ‡ The experiment described in the previous section was executed at the NASA Langley 12-Foot Low-Speed Tunnel.

The static test points, consisting of 𝛼 , 𝛽 , 𝛿 , and 𝛿 settings, were run using an automated data collection procedure 𝑤 𝑤 1 2 developed for previous testing using DOE/RSM techniques [ 12 , 15 ]. The PTI excitations for control surface and propulsor commands were run in a continuous loop while traversing through the static test matrix, but data were only collected at each static test point. Experiments used for this paper were executed at a dynamic pressure of 3.5 psf (freestream airspeed of 54.3 ft/s at standard sea-level conditions), with the test factor ranges shown in Table 2. This condition represents a high-speed transition phase of flight for the LA-8 aircraft. Note that multiple dynamic pressure settings need to be tested to develop a full-envelope aero-propulsive model.

Data collected for model identification included sting orientation, control surface deflection angles, propulsor rotational speeds, and applied forces and moments. The six force and moment components were measured using a strain gage balance, and a propeller rotational speed measurement was provided by the electronic speed controller. Direct control surface position measurements were not available, so a calibration curve was developed for each control surface to convert the desired deflection angles to commanded pulse width modulation (PWM) signal values for testing. Data were collected at a sample rate of 50 Hz for 60 seconds at each static test point, with less overall time expected to be needed for model identification.

Table 2 Test factor ranges at ¯ 𝒒 = 3 . 5 psf ( 𝑽 = 54 . 3 ft/s) Factor(s) Units Minimum Maximum 𝛼 deg − 6 + 6 𝛽 deg − 5 + 5 𝛿 , 𝛿 deg 0 + 25 𝑤 𝑤 1 2 𝛿 , 𝛿 , 𝛿 , 𝛿 deg − 15 + 15 𝑒 𝑒 𝑒 𝑒 1 2 3 4 𝛿 , 𝛿 , 𝛿 , 𝛿 deg 0 + 20 𝑓 𝑓 𝑓 𝑓 1 2 3 4 𝛿 , 𝛿 deg − 15 + 15 𝑟 𝑟 1 2 𝜂 , 𝜂 , ..., 𝜂 𝜇 s 1400 1500 1 2 8 The dynamic nature of the data used for this study required multiple signal processing steps to condition the data for model identification. First, the propulsor speed and control surface deflection angle signals were found to have a time lag relative to the force and moment measurements. The time lag was determined by finding the peak of the cross-correlation function between the control surface and propulsor signals and the dominant force or moment component where the ‡ Information available online at https://researchdirectorate.larc.nasa.gov/12-foot-low-speed-tunnel-12-ft-lst/ [accessed 29 October 2021] control effects manifested. The cross-correlation was computed using data collection runs executed with a reduced number of control effectors being dynamically excited to clarify time lag estimation. After determining the time lag, time skew corrections were made for the signals used for model identification.

The force and moment signals were also found to contain significant measurement noise and suspected structural modes. An analog sixth-order low-pass Butterworth anti-alias filter with a cutoff frequency of 10 Hz was applied to the balance measurements prior to sampling. Nonetheless, residual frequency content outside of the excitation frequencies was observed in the data. To mitigate these effects, the model identification from the dynamic data was performed using only information content contained in the dynamic excitation frequency range, as will be discussed further in Sec. V.

§ For time-domain model validation, data were smoothed using a zero-phase shift digital filtering technique, with a digital sixth-order low-pass Butterworth filter with a cutoff frequency of 3.1 Hz applied both forward and backward in time [ 51 ]. The cutoff frequency was selected to preserve lower frequency information associated with the control effector excitation frequencies, but reject most higher frequency noise. The frequency response of the zero-phase shift digital smoother is shown in Fig. 11 (phase is not shown because it is zero across the frequency range). The linear magnitude at 1.756 Hz (the highest PTI excitation frequency) is greater than 0.999 and − 20 dB attenuation is achieved at 3.7 Hz.

Fig. 11 Frequency response for the zero-phase shift digital low-pass Butterworth smoother.

Sample measured and smoothed force and moment data in the time and frequency domain are shown in Figs. 12-13.

¶ The power spectra shown in Fig. 13 were computed using Welch’s method. The power spectra of the original force and moment measurements show evidence of structural modes and measurement noise above the PTI excitation frequencies, which are avoided by only performing model identification using the frequency band corresponding to the control effector excitation frequencies, indicated by the shaded areas on the plot. The figures also show that the smoother suppresses the structural modes and measurement noise above the PTI excitation frequencies while preserving the lower-frequency information attributed to the dynamic excitation inputs.

Furthermore, the reaction torque effects in the rolling and yawing moment signals must be accounted for because the propulsors are being dynamically excited [ 17 , 38 ]. The combination of aero-propulsive moments and moments due to the angular momentum rate of the propulsors is measured by the balance. Therefore, to model the aero-propulsive moments, the angular momentum rate effects must be removed. This is accomplished using the expressions 𝑛 𝑛 𝑝 𝑝 ∑︁ ∑︁ ¤ ¤ 𝐿 = 𝐿 + ¤ ℎ = 𝐿 + ℎ cos 𝛿 = 𝐿 + 𝐼 Ω cos 𝛿 (9) 𝐸 𝑥 𝐸 𝑝 𝑤 𝐸 𝑝 𝑝 𝑤 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘 = 1 𝑘 = 1 𝑛 𝑛 𝑝 𝑝 ∑︁ ∑︁ ¤ ¤ 𝑁 = 𝑁 + ¤ ℎ = 𝑁 + ℎ (− sin 𝛿 ) = 𝑁 + 𝐼 Ω (− sin 𝛿 ) (10) 𝐸 𝑧 𝐸 𝑝 𝑤 𝐸 𝑝 𝑝 𝑤 𝑘 𝑘 𝑘 𝑘 𝑘 𝑘 = 1 𝑘 = 1 where 𝐿 and 𝑁 are the uncorrected balance measurements. The angular momentum of the 𝑘 th propulsor about 𝐸 𝐸 its axis of rotation is ℎ = 𝐼 Ω , where 𝐼 is the moment of inertia of the rotating portion of the propulsor and 𝑝 𝑝 𝑝 𝑝 𝑘 𝑘 𝑘 𝑘 Ω = 2 𝜋𝑛 is the rotation rate in radians per second, with clockwise rotation when viewed from behind the rotating 𝑝 𝑘 𝑘 propulsor being positive. For use in the above equations, the angular momentum rate for each propulsor is rotated into the aircraft body axes through the corresponding wing angle 𝛿 , and then summed to compute the components of net 𝑤 𝑘 ¤ ¤ angular momentum rate for all 𝑛 propulsors ( ℎ , ℎ ). Note that if the vehicle orientation was dynamically changing 𝑝 𝑥 𝑧 during data collection, propulsor gyroscopic effects would also need to be taken into account.

§ Information available online at https://www.mathworks.com/help/signal/ref/filtfilt.html [retrieved 5 April 2022] ¶ Information available online at https://www.mathworks.com/help/signal/ref/pwelch.html [retrieved 5 April 2022]

Fig. 12 Sample balance force and moment measurement histories ( 𝒗 = 𝒘 = 0 ft/s, 𝜹 = 𝜹 = 12 . 5 deg).

𝒘 𝒘 1 2 Fig. 13 Sample balance force and moment measurement power spectra ( 𝒗 = 𝒘 = 0 ft/s, 𝜹 = 𝜹 = 12 . 5 deg).

𝒘 𝒘 1 2 V. Aero-Propulsive Modeling Approach Aero-propulsive modeling for this effort focuses on developing a polynomial representation of the aero-propulsive forces and moments as a function of vehicle state and control variables. Factors under test, or close variants, are defined as explanatory variables, and the model is identified from the wind tunnel data collected using the techniques described in Secs. III-IV. Aero-propulsive modeling for tilt-wing, distributed electric propulsion aircraft requires a different approach compared to conventional fixed-wing and rotary-wing aircraft modeling approaches. eVTOL vehicles can be considered a fixed-wing/rotary-wing hybrid suggesting that a combination of modeling approaches will facilitate suitable model development. Accordingly, the modeling variable formulation described in Ref. [ 16 ], and applied here, merges appropriate fixed-wing and rotary-wing modeling attributes with new strategies to model the LA-8 and other similar vehicles.

The modeled responses are the dimensional body-axis aero-propulsive forces 𝑋, 𝑌 , 𝑍 in lbf and moments 𝐿, 𝑀, 𝑁 in ft · lbf. The explanatory variables are defined as the body-axis velocity components 𝑣, 𝑤 in ft/s, propeller speeds 𝑛 , 𝑛 , ..., 𝑛 in revolutions per second, wing angles 𝛿 , 𝛿 in radians, elevon deflections 𝛿 , 𝛿 , 𝛿 , 𝛿 in radians, 1 2 8 𝑤 𝑤 𝑒 𝑒 𝑒 𝑒 1 2 1 2 3 4 flap deflections 𝛿 , 𝛿 , 𝛿 , 𝛿 in radians, and ruddervator deflections 𝛿 , 𝛿 in radians. Note that it is important to 𝑓 𝑓 𝑓 𝑓 𝑟 𝑟 1 2 3 4 1 2 perform modeling with explanatory variables centered on a reference value to maintain low correlation among candidate regressors in time or spacial domain analysis and to avoid spectral leakage in frequency domain analysis [ 52 ]. For this work, the centering reference values were the explanatory variable mean values in the data used for modeling.

Identified models typically use the test factors as the explanatory variables; however, due to unique characteristics of eVTOL aircraft, test facility integration limitations, and operational convenience, analysis was instead performed by defining certain alternative explanatory variables for modeling. For example, as mentioned in Secs. III-IV, wind tunnel testing was performed by varying 𝛼 and 𝛽 directly for operational convenience. Although testing was performed with experimental factors of 𝛼 and 𝛽 , modeling was performed using body-axis velocity components 𝑣 and 𝑤 to make the modeling strategy valid from hover through forward flight. Since body-axis velocity components are closely related to airflow angles, this variable change does not affect the ability to identify a model from the data because the regressors retain their low correlation when testing is conducted at a constant forward speed [ 16 ]. Similarly, testing was performed by varying motor PWM commands ( 𝜂 , 𝜂 , ... , 𝜂 ), but modeling was performed using measured propeller rotational 1 2 8 speeds ( 𝑛 , 𝑛 , ... , 𝑛 ) because propeller speed is more relevant to describe propeller aerodynamics.

1 2 8 The following subsections provide an overview of the methods used for model structure development, parameter estimation, and model validation. The model structure identification and parameter estimation methods used for this † work were adapted from the System IDentification Programs for AirCraft (SIDPAC) software toolbox. The model identification approach occurs in two stages. First, the continuous time series data at each static test point, or set point, are analyzed using frequency domain methods to determine the model structure and parameter estimates associated with the dynamically changing explanatory variables ( 𝑛 , 𝑛 , ..., 𝑛 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 ) using the 1 2 8 𝑒 𝑒 𝑒 𝑒 𝑓 𝑓 𝑓 𝑓 𝑟 𝑟 1 2 3 4 1 2 3 4 1 2 aero-propulsive force and moments ( 𝑋 , 𝑌 , 𝑍 , 𝐿 , 𝑀 , 𝑁 ) as the response variables. Then, in the second step, the parameter estimates for the models identified from the dynamically changing explanatory variables become the responses to be modeled as a function of the static explanatory variables ( 𝑣 , 𝑤 , 𝛿 , 𝛿 ). Model structure determination and parameter 𝑤 𝑤 1 2 estimation are performed in a weighted least-squares formulation using the parameter uncertainty estimates from the first model identification step. This framework ultimately yields a two-layer model to predict the aero-propulsive force and moments as a function of the static and dynamic explanatory variables, as will be explained further throughout the remainder of the paper.

To the best knowledge of the authors, this paper presents the first development and application of a compound modeling strategy leveraging frequency-domain and weighted least-squares regression methods for model structure determination to combine the static and dynamic data information content for aggregate model identification. Previous work presented in Refs. [ 26 , 27 ] implemented a similar two-layered aerodynamic modeling approach using time-domain ordinary least-squares regression methods. Reference [ 33 ] developed a global aerodynamic model using a tabulated representation of local parameters estimated using the same frequency-domain equation-error method used in this work.

A. Part I: Frequency-Domain Modeling at Each Set Point The first step in the modeling process involves analyzing the continuous time series data collected during each individual set point. As mentioned in Sec. IV, the sampled data contained measurement noise and structural dynamics outside of the range of the PTI excitations. Application of frequency domain estimation techniques allows model identification to be performed only in the frequency band associated with the input excitations, which effectively smooths the modeling data [ 18 , 52 , 53 ]. Additional benefits of model identification in the frequency domain include accurate parameter uncertainty estimation, increased computational speed, estimation of nearly unbiased parameter estimates when the regressors contain noise, and least-squares weighting based on frequency components as opposed to individual data points.

Model identification was performed using the equation-error method formulated in the frequency domain [ 18 , 52 , 54 ].

The regressor and response data are detrended and transformed into the frequency domain using a Fourier transform technique leveraging time-domain cubic interpolation and the chirp- 𝑧 transform to produce a high accuracy transform with an arbitrary frequency range and resolution [ 18 , 55 ]. For this study, the transform frequency range was selected to match the excitation input design with a fundamental period of 𝑇 = 180 s (see Sec. III.B). The Fourier transform † Information available online at https://software.nasa.gov/software/LAR-16100-1 [accessed 29 October 2021] frequencies were set between 𝑓 = 0 . 05 Hz and 𝑓 = 1 . 756 Hz with a frequency resolution of Δ 𝑓 = 0 . 00556 Hz, min max resulting in 𝑀 = 308 transform frequencies. A parameter sensitivity study indicated that a coarser frequency resolution resulted in parameter estimation errors, as explained in [ 52 ], whereas a finer frequency resolution provided minimal additional parameter estimation accuracy.

1. Model Structure Determination Multivariate orthogonal function (MOF) modeling, described in Refs. [ 18 , 56 ], was selected as the model structure identification technique for this work. Although MOF modeling was originally developed for model structure determination using real-valued data in the time domain (or spacial domain), the method can also be formulated in the frequency domain, while still allowing nonlinear model terms to be considered for the model structure [ 52 , 57 , 58 ].

The key modifications, developed in Refs. [ 57 , 58 ], are applied to the candidate regressor and response variable data used in the MOF modeling algorithm. The nonlinear candidate regressors are assembled in the time domain using detrended explanatory variables and are subsequently detrended again before transforming the candidate regressor data into the frequency domain. After the Fourier transform is applied to the candidate regressor data, the real and imaginary ˜ components of the complex candidate regressor matrix 𝑿 are concatenated to form a real vector: " # ˜ 𝑅𝑒 ( 𝑿 ) ′ 𝑿 = (11) ˜ 𝐼𝑚 ( 𝑿 ) Similarly, the detrended response variable vector is transformed to the frequency domain and then the real and imaginary components of the complex response variable vector ˜ 𝒛 are concatenated to form a real vector: " # 𝑅𝑒 ( ˜ 𝒛 ) ′ 𝒛 = (12) 𝐼𝑚 ( ˜ 𝒛 ) ′ ′ The real candidate regressor matrix 𝑿 and response variable vector 𝒛 assembled from complex data are then used in the standard MOF modeling algorithm, which is now briefly summarized.

The multivariate orthogonal function (MOF) modeling algorithm [ 18 , 56 ] starts by orthogonalizing a predefined set of candidate regressors using an algorithm such as Gram-Schmidt orthogonalization or QR decomposition. Orthogonal regressors are convenient for model structure development because of the ability to independently assess the potential of the orthogonalized candidate regressors to model the response variable. This facilitates selecting only the model terms that significantly contribute to model effectiveness. The orthogonal regressors are then ranked from greatest to least decrease in the mean squared fit error (MSFE): 𝑇 MSFE = ( 𝒛 − ˆ 𝒚 ) ( 𝒛 − ˆ 𝒚 ) (13) 𝑁 In other words, this ranks the orthogonalized regressors from highest to lowest ability to improve the model. Candidate orthogonal regressors are brought into the model structure in this order.

Deciding which terms to include in the final model can then be done using one or more statistical metrics. A common threshold for MOF modeling is to minimize the predicted squared error (PSE) [ 18 , 59 ]. The PSE is the sum of the MSFE and a model complexity penalty related to the number of terms included in the model, 𝑝 PSE = MSFE + 𝜎 (14) max 𝑁 where 𝑝 is the number of terms in the current model structure and 𝜎 is an estimate of the upper-bound of mean max squared error for the model prediction of data not used to develop the model. The quantity 𝜎 is estimated here using max the variance between the measured response 𝒛 and mean measured response ¯ 𝒛 : 𝑁 ∑︁ 2 2 𝜎 = [ 𝒛 ( 𝑖 ) − ¯ 𝒛 ] (15) max 𝑁 − 1 𝑖 = 1 After the orthogonal regressors were ranked by their ability to reduce the MSFE, the cutoff for model term addition was chosen to be the candidate model term where the PSE was minimized. When the orthogonalized regressors are ranked as stated above, the PSE metric is guaranteed to have a single global minimum [18].

For this work, MOF modeling was applied separately to the time series data collected at each individual set point.

Although the model structure determined for each response at each set point was similar, there were some differences due to the experimental nature of the study and having a rough estimate of 𝜎 . The model structure for each individual max response used to identify the final model parameter estimates at each set point was constrained to be the same. The model terms selected for the final model for each response were the terms that the MOF algorithm included in the model for a majority of the set points.

2. Parameter Estimation After determining the model structure for each response, the final parameter values were estimated using ordinary least-squares regression with the complex regressor and response data. Ordinary least-squares regression is used to ˜ estimate a vector of 𝑝 unknown model parameters in a vector 𝜽 for a given model ˜ 𝒚 = 𝑿𝜽 , where ˜ 𝒚 is the length 𝑀 ˜ complex model response vector and 𝑿 is a 𝑀 × 𝑝 matrix consisting of column vectors of the complex regressors assumed to be error-free [ 18 ]. The regression equation, including the complex measured response variable ˜ 𝒛 corrupted by constant variance, zero-mean, and uncorrelated complex error ˜ 𝝂 , is given as: ˜ ˜ 𝒛 = 𝑿𝜽 + ˜ 𝝂 (16) For complex least-squares parameter estimation, the optimal estimate of the unknown parameters 𝜽 is determined by minimizing the cost function:

† ˜ ˜ 𝐽 ( 𝜽 ) = ˜ 𝒛 − 𝑿𝜽 ˜ 𝒛 − 𝑿𝜽 (17) It follows that the solution to compute an optimal estimate of the unknown real-valued parameters is h  i   − 1 † † ˆ ˜ ˜ ˜ 𝜽 = 𝑅𝑒 𝑿 𝑿 𝑅𝑒 𝑿 ˜ 𝒛 (18) ˆ where 𝜽 is a vector of 𝑝 estimated parameters. The modeled response variable vector is: ˆ ˜ ˆ ˜ 𝒚 = 𝑿 𝜽 (19) ˆ ˆ A length 𝑝 vector of standard errors 𝒔 ( 𝜽 ) corresponding to the estimated parameters 𝜽 is given as: √︄           h  i † − 1 † ˆ ˆ ˆ ˜ ˜ 𝒔 ( 𝜽 ) = 𝑅𝑒 ˜ 𝒛 − ˜ 𝒚 ˜ 𝒛 − ˜ 𝒚 diag 𝑅𝑒 𝑿 𝑿 (20) 2 𝑇 ( 𝑓 − 𝑓 ) max min ˆ This form of 𝒔 ( 𝜽 ) accounts for the fact that a frequency range narrower than zero to the Nyquist frequency is used for analysis [52].

After complex least-squares parameter estimation is completed, an additional step is needed to identify the bias term in a model equation because the detrended data used to estimate the model parameters contain only dynamic ˆ ˆ information [ 18 , 52 ]. The bias parameter estimate 𝜃 is found as the mean value of ( 𝒛 − 𝑿 𝜽 ) , where 𝒛 is the measured 𝑜 response variable in the time domain, 𝑿 is a matrix consisting of column vectors of the regressors in the time domain, ˆ and 𝜽 is the model parameter vector estimated previously using complex least-squares regression. The bias parameter standard errors were estimated accounting for colored residuals using the method described in Refs. [ 18 , 60 ] to compute a more representative estimate of parameter uncertainty.

B. Part II: Weighed Least-Squares Aggregate Modeling The results from Part I of the model identification yield parameter estimates and uncertainty estimates associated with the dynamic explanatory variables at each set point. Part II of the model identification procedure developed for this work uses the Part I modeling results to develop a response surface model for each parameter identified in Part I as a function of the static explanatory variables held constant at each set point. Contrary to Part I, the explanatory and response variable data for this step are real and have no meaningful time dependence. Also, each response data point has an associated uncertainty estimate accurately determined using the methods implemented in Part I. This suggests using a weighted least-squares framework for identification of the model structure and parameter estimates. The weights are stored in a diagonal matrix 𝑾 with the diagonal elements being the inverse of the variance associated with each response observation (i.e., the uncertainty estimates obtained using the methods described in Sec. V.A.2):     𝑤 0 . . . 0 1 / 𝑠 0 . . . 0         0 𝑤 . . . 0 0 1 / 𝑠 . . . 0         𝑾 = = (21) . . . . . .

 .   .  . . . . . . . .

 .   .  . . . . . .

        0 0 . . . 𝑤 0 0 . . . 1 / 𝑠 𝑁    𝑁  1. Model Structure Determination The MOF modeling algorithm [ 18 , 56 ], summarized in Sec. V.A.1, is again used for model structure determination in Part II of model identification; however, certain modifications are needed to accommodate the fact that the variance estimates associated with each response measurement are known and not constant, which leads to a weighted least-squares estimation problem. First, after generation of the candidate regressors, but before orthogonalization, the candidate regressors and response data for each data point are multiplied by the square root of the diagonal elements of the diagonal weight matrix 𝑾 , or the inverse of the estimated standard error at each data point, as: √ √ √     𝑤 𝑥 𝑤 𝑥 . . . 𝑤 𝑥 𝑥 / 𝑠 𝑥 / 𝑠 . . . 𝑥 / 𝑠 1 11 1 12 1 1 𝑝 11 1 12 1 1 𝑝 1     √ √ √     𝑤 𝑥 𝑤 𝑥 . . . 𝑤 𝑥 𝑥 / 𝑠 𝑥 / 𝑠 . . . 𝑥 / 𝑠  2 21 2 22 2 2 𝑝   21 2 22 2 2 𝑝 2  ′′     𝑿 = = (22) . . . . . .

. .

    . . . . . . . .

. .

 . . .   . . .      √ √ √     𝑤 𝑥 𝑤 𝑥 . . . 𝑤 𝑥 𝑥 / 𝑠 𝑥 / 𝑠 . . . 𝑥 / 𝑠 𝑁 𝑁 1 𝑁 𝑁 2 𝑁 𝑁 𝑝 𝑁 1 𝑁 𝑁 2 𝑁 𝑁 𝑝 𝑁     √     𝑤 𝑧 𝑧 / 𝑠 1 1 1 1     √     𝑤 𝑧 𝑧 / 𝑠     2 2 2 2 ′′     𝒛 = = (23) . .

    . .

 .   .      √     𝑤 𝑧 𝑧 / 𝑠 𝑁 𝑁 𝑁 𝑁     Recall that the response variable data in this step are the model parameter estimates from the frequency-domain modeling performed at each set point in Part I. Additionally, the PSE expression is reformulated as: 1 𝑝 𝑇 PSE = ( 𝒛 − ˆ 𝒚 ) 𝑾 ( 𝒛 − ˆ 𝒚 ) + 𝑐 (24) 𝑁 𝑁 Here, the weight matrix 𝑾 serves as a model-independent error variance estimate. Due to the scaling by 𝑾 , the response 2 2 measurement error variance ˆ 𝜎 is equal to one. Thus, inclusion of the scale factor 𝑐 = 𝑐 ˆ 𝜎 = 𝜎 is a similar concept to max the scale factor used in Refs. [ 61 , 62 ] when a model-independent measurement error variance estimate was available from wind tunnel testing. A value of 𝑐 = 10 was found to be a good value to develop models with minimum prediction error.

2. Parameter Estimation After determining the model structure for each response, weighted least-squares regression is used to compute the final parameter estimates in ordinary regressor space. Weighted least-squares regression estimates a vector 𝜽 of 𝑝 unknown model parameters for a given model 𝒚 = 𝑿𝜽 , where 𝒚 is the length 𝑁 model response vector and 𝑿 is a 𝑁 × 𝑝 matrix consisting of column vectors of regressors assumed to be measured without error [ 18 , 20 , 49 ]. The regression equation, including a measurement of the response variable 𝒛 , corrupted by zero-mean, uncorrelated error 𝝂 with non-constant variance, is: 𝒛 = 𝑿𝜽 + 𝝂 (25) For weighted least-squares parameter estimation, the optimal estimate of the unknown parameters 𝜽 is determined by minimizing the cost function 𝑇 𝐽 ( 𝜽 ) = ( 𝒛 − 𝑿𝜽 ) 𝑾 ( 𝒛 − 𝑿𝜽 ) (26)

where 𝑾 is the diagonal matrix given in Eq. (21) and each diagonal element is the inverse of the variance associated with each response data point which, here, is a parameter estimate from Part I. Consequently, response data with lower uncertainty have an increased influence for computing the least-squares solution. It follows that the solution to compute an optimal estimate of the unknown parameters is   − 1 𝑇 𝑇 ˆ 𝜽 = 𝑿 𝑾 𝑿 𝑿 𝑾 𝒛 (27) ˆ ˆ where 𝜽 is a vector of 𝑝 estimated parameters. The length 𝑝 vector of standard errors 𝒔 ( 𝜽 ) corresponding to the ˆ estimated parameters 𝜽 is given as: √︂ h i
− 1 𝑇 ˆ 𝒔 ( 𝜽 ) = diag 𝑿 𝑾 𝑿 (28) An alternative way to compute the weighted least-squares solution is to multiply each normal equation by the square root of the respective diagonal element in 𝑾 (or the inverse of the measurement standard error) and then use ordinary ′′ least-squares regression to compute the parameter estimates. In other words, the same solution is obtained using 𝑿 ′′ and 𝒛 from Eqs. (22)-(23) in an ordinary least-squares regression estimator.

C. Final Model Validation Model fit metrics and modeling residuals alone do not provide information about the model predictive capability.

Assessment of model performance using validation data not used for modeling provides a more reliable estimate of model prediction accuracy. Final model validation is performed by comparing the measured response for validation data to the response predicted by the model for the same explanatory variable inputs. Further assessment is performed by analyzing the prediction residuals between the measured and predicted response, 𝒆 = 𝒛 − ˆ 𝒚 . Residuals and their statistical properties can be given further interpretability by normalization. The error normalization metric used in this work is the range of response variable measurements used to develop the local model, range ( 𝒛 ) = 𝒛 − 𝒛 .

max min Range normalization provides a fair comparison between prediction error metrics for different response variables used for aircraft modeling because longitudinal responses are generally biased above or below zero and lateral-directional responses are generally centered about zero. The normalized residual vector is defined as: 𝒛 − ˆ 𝒚 ∗ 𝒆 = (29) range ( 𝒛 ) Similarly, the normalized root-mean-square modeling error (NRMSE) is defined as: √︄ 𝑇 1 ( 𝒛 − ˆ 𝒚 ) ( 𝒛 − ˆ 𝒚 ) NRMSE = (30) range ( 𝒛 ) 𝑁 For this study, the goal was to develop models minimizing prediction error, where a value of approximately 5% or less for NRMSE was considered to be adequate based on analyst judgment and previous aerodynamic modeling studies conducted in the experimental facility used for wind tunnel testing.

VI. Results This section presents sample results for the aero-propulsive models identified for the LA-8 aircraft. The results presented here only consider models identified at ¯ 𝑞 = 3 . 5 psf, which corresponds to an airspeed of 54.3 ft/s at standard sea-level conditions. For eVTOL vehicle simulations, models valid at multiple dynamic pressure (or airspeed) settings throughout the flight envelope are needed [15, 16].

A. Part I Modeling Results Following the aero-propulsive modeling approach described in Sec. V, for Part I, a model structure for each force and moment component was developed using the time series data collected at each set point as a function of the dynamically changing explanatory variables. The candidate regressors included in the model structure for a majority of the individual set point models were included in a uniform, final aero-propulsive model structure used to identify the local parameter values for each set point. In other words, only regressors modeling a significant portion of the variation in the response variables throughout the majority of the static test variable space were included, in accordance with the model structure determination strategy discussed in Sec. V.A.1. The individual set point model structures for each aero-propulsive force and moment component were determined to be: 𝑋 = 𝑋 𝛿 + 𝑋 𝛿 + 𝑋 𝛿 + 𝑋 𝛿 + 𝑋 𝛿 + 𝑋 𝛿 + 𝑋 𝛿 + 𝑋 𝛿 + 𝑋 𝑛 + ...

𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝑛 1 𝑒 1 𝑒 2 𝑒 3 𝑒 4 𝑓 1 𝑓 2 𝑓 3 𝑓 4 1 1 2 3 4 1 2 3 4 (31) 𝑋 𝑛 + 𝑋 𝑛 + 𝑋 𝑛 + 𝑋 𝑛 + 𝑋 𝑛 + 𝑋 𝑛 + 𝑋 𝑛 + 𝑋 𝑛 2 𝑛 3 𝑛 4 𝑛 5 𝑛 6 𝑛 7 𝑛 8 𝑜 2 3 4 5 6 7 8 𝑌 = 𝑌 𝛿 + 𝑌 𝛿 + 𝑌 𝛿 + 𝑌 𝛿 + 𝑌 𝛿 + 𝑌 𝛿 + 𝑌 𝑛 + 𝑌 𝑛 + 𝑌 𝑛 + ...

𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑟 𝛿 𝑟 𝑛 1 𝑛 2 𝑛 3 𝑓 1 𝑓 2 𝑓 3 𝑓 4 𝑟 1 𝑟 2 1 2 3 1 2 3 4 1 2 (32) 𝑌 𝑛 + 𝑌 𝑛 + 𝑌 𝑛 + 𝑌 𝑛 + 𝑌 𝑛 + 𝑌 𝑛 4 𝑛 5 𝑛 6 𝑛 7 𝑛 8 𝑜 4 5 6 7 8 𝑍 = 𝑍 𝛿 + 𝑍 𝛿 + 𝑍 𝛿 + 𝑍 𝛿 + 𝑍 𝛿 + 𝑍 𝛿 + 𝑍 𝛿 + 𝑍 𝛿 + 𝑍 𝛿 + ...

𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑟 𝑒 1 𝑒 2 𝑒 3 𝑒 4 𝑓 1 𝑓 2 𝑓 3 𝑓 4 𝑟 1 1 2 3 4 1 2 3 4 1 (33) 𝑍 𝛿 + 𝑍 𝑛 + 𝑍 𝑛 + 𝑍 𝑛 + 𝑍 𝑛 + 𝑍 𝑛 + 𝑍 𝑛 + 𝑍 𝑛 + 𝑍 𝑛 + 𝑍 𝛿 𝑟 𝑛 1 𝑛 2 𝑛 3 𝑛 4 𝑛 5 𝑛 6 𝑛 7 𝑛 8 𝑜 𝑟 2 1 2 3 4 5 6 7 8 𝐿 = 𝐿 𝛿 + 𝐿 𝛿 + 𝐿 𝛿 + 𝐿 𝛿 + 𝐿 𝛿 + 𝐿 𝛿 + 𝐿 𝛿 + 𝐿 𝛿 + 𝐿 𝑛 + ...

𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝑛 1 𝑒 1 𝑒 2 𝑒 3 𝑒 4 𝑓 1 𝑓 2 𝑓 3 𝑓 4 1 1 2 3 4 1 2 3 4 (34) 𝐿 𝑛 + 𝐿 𝑛 + 𝐿 𝑛 + 𝐿 𝑛 + 𝐿 𝑛 + 𝐿 𝑛 + 𝐿 𝑛 + 𝐿 𝑛 2 𝑛 3 𝑛 4 𝑛 5 𝑛 6 𝑛 7 𝑛 8 𝑜 2 3 4 5 6 7 8 𝑀 = 𝑀 𝛿 + 𝑀 𝛿 + 𝑀 𝛿 + 𝑀 𝛿 + 𝑀 𝛿 + 𝑀 𝛿 + 𝑀 𝛿 + 𝑀 𝛿 + 𝑀 𝛿 + ...

𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑟 𝑒 1 𝑒 2 𝑒 3 𝑒 4 𝑓 1 𝑓 2 𝑓 3 𝑓 4 𝑟 1 1 2 3 4 1 2 3 4 1 (35) 𝑀 𝛿 + 𝑀 𝑛 + 𝑀 𝑛 + 𝑀 𝑛 + 𝑀 𝑛 + 𝑀 𝑛 + 𝑀 𝑛 + 𝑀 𝑛 + 𝑀 𝑛 + 𝑀 𝛿 𝑟 𝑛 1 𝑛 2 𝑛 3 𝑛 4 𝑛 5 𝑛 6 𝑛 7 𝑛 8 𝑜 𝑟 2 1 2 3 4 5 6 7 8 𝑁 = 𝑁 𝛿 + 𝑁 𝛿 + 𝑁 𝛿 + 𝑁 𝛿 + 𝑁 𝛿 + 𝑁 𝛿 + 𝑁 𝛿 + 𝑁 𝛿 + 𝑁 𝛿 + ...

𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑒 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑓 𝛿 𝑟 𝑒 1 𝑒 2 𝑒 3 𝑒 4 𝑓 1 𝑓 2 𝑓 3 𝑓 4 𝑟 1 1 2 3 4 1 2 3 4 1 (36) 𝑁 𝛿 + 𝑁 𝑛 + 𝑁 𝑛 + 𝑁 𝑛 + 𝑁 𝑛 + 𝑁 𝑛 + 𝑁 𝑛 + 𝑁 𝑛 + 𝑁 𝑛 + 𝑁 𝛿 𝑟 𝑛 1 𝑛 2 𝑛 3 𝑛 4 𝑛 5 𝑛 6 𝑛 7 𝑛 8 𝑜 𝑟 2 1 2 3 4 5 6 7 8 Although pure quadratic and two-factor interaction candidate regressors were considered for Part I model structure determination, the final model describing the variation of the force and moment components at each test point is linear in the dynamic explanatory variables. This makes sense because the variables were excited over a relatively small range of values. The method, however, is also capable of modeling nonlinear effects, such as propulsor-control surface interactions, which can be considered for future studies with a greater excitation range. As was found in Ref. [ 16 ], tilt-wing vehicles exhibit significant nonlinearity with body-axis velocity and wing angle variation, which is captured in the Part II model identification strategy.

ˆ Figure 14 shows the 𝑍 parameter estimates and error bars of ± 2 𝑠 ( 𝜽 ) computed using the parameter estimation techniques described in V.A.2 against the amount of data collection time used for modeling at a sample set point.

The parameters appear to converge to reasonably consistent values by the time 40 seconds has elapsed, which was also determined to be a satisfactory data collection time in a simulated flight test system identification study for the LA-8 aircraft [ 17 ]. A similar parameter convergence rate was observed for the parameters corresponding to the other aero-propulsive forces and moments, and at different set points. In view of the parameter convergence analysis shown in Fig. 14, the findings in Ref. [ 17 ], and with the knowledge that 40 seconds is the amount of time needed to complete two full periods of the lowest frequency sinusoidal component of the multisine signals ( 𝑓 = 0 . 05 Hz), 40 seconds was min selected as the amount of data collection time to use for modeling at each set point. Although the actual data collection time at each set point was greater than 40 seconds, providing additional information for modeling, a goal of this study was to inform future wind tunnel testing efforts with practical, efficient data collection strategies. Consequently, modeling was performed with the recommended data collection time of 40 seconds based on the present reasoning, as opposed to using all available data. Note that the parameter estimates shown in Fig. 14 contain asymmetries that are not apparent from the LA-8 vehicle configuration, such as differences in the reflected flap parameter estimates. This is a result of manufacturing differences between the clockwise and counterclockwise propellers, which resulted in a significant difference in thrust production between the propeller variants [ 37 ]. Since the propulsion-only and propulsion-airframe interaction effects are significant, this propulsion asymmetry is manifested in many of the model terms [16].

Figure 15a shows the model fit in the frequency domain using a modeling time of 40 seconds at a sample set point.

The corresponding control surface deflection and propulsor rotational speed signals are shown in Fig. 16. A good model fit is observed for each response. The 𝑅 metric shown on the subplot for each response is approximately 95% or higher, indicating that most of the variation of the response variable about its mean value is characterized by the model.

Figure 15b shows the corresponding time-domain model fit compared to the smoothed, measured aero-propulsive forces and moments. The modeled responses are close to the measured responses, indicating that the model is able to describe a large amount of the variation in each response. The modeling performance at other set points was similar.

Fig. 14 Variation of 𝒁 parameter estimates with data collection time used for modeling at a sample set point ( 𝒗 = 𝒘 = 0 ft/s, 𝜹 = 𝜹 = 12 . 5 deg).

𝒘 𝒘 1 2 (a) Model fit in the frequency domain (b) Model fit in the time domain Fig. 15 Comparison of response data and model fit at a sample set point ( 𝒗 = 𝒘 = 0 ft/s, 𝜹 = 𝜹 = 12 . 5 deg).

𝒘 𝒘 1 2 Fig. 16 Control effector signals at a sample set point ( 𝒗 = 𝒘 = 0 ft/s, 𝜹 = 𝜹 = 12 . 5 deg).

𝒘 𝒘 1 2 B. Part II Modeling Results After identification of the model parameters associated with the dynamic explanatory variables at each set point, models were developed to characterize the variation in those model parameters across the static variable space. As discussed in Sec. V.B, the parameter estimates associated with the dynamic explanatory variables were treated as the response variables for Part II of the modeling approach. In Part II, the static test variables associated with each set point ( 𝑣 , 𝑤 , 𝛿 , 𝛿 ) were the explanatory variables. The models for each Part II response variable were developed using 𝑤 𝑤 1 2 weighted least-squares regression techniques because uncertainty estimates for each response variable data point were available from the Part I parameter estimation results.

As an example of a polynomial model structure developed in this modeling stage, the final model structure identified for 𝑍 using MOF modeling in a weighted least-squares formulation was: 𝛿 𝑓 𝑍 = 𝑍 𝑤 + 𝑍 𝛿 + 𝑍 𝑤𝛿 + 𝑍 𝛿 + 𝑍 (37) 𝛿 𝛿 𝛿 𝑤 𝛿 𝑤 𝛿 𝛿 𝑓 𝑓 𝑓 1 𝑓 1 𝑓 2 𝑤 𝑓 1 1 1 1 1 1 1 ( 𝑤 ) ( 𝛿𝑤 ) ( 𝑤 𝛿𝑤 ) ( 𝛿 ) ( 𝑜 ) 𝑤 1 1 The bias parameters (e.g., 𝑍 ) had the largest number of model terms because they characterize the direct effects of the 𝑜 static explanatory variables on the aero-propulsive forces and moments. For example, a model structure identified for the 𝑍 model was: 𝑜 2 2 𝑍 = 𝑍 𝑤 + 𝑍 𝑣 + 𝑍 𝛿 + 𝑍 𝛿 + 𝑍 𝑤 + 𝑍 𝑤𝑣 + 𝑍 𝑣 + ...

𝑜 𝑜 𝑜 𝑜 𝑤 𝑜 𝑤 𝑜 𝑜 𝑜 ( 𝑤 ) ( 𝑣 ) ( 𝛿𝑤 ) 1 ( 𝛿𝑤 ) 2 2 ( 𝑤𝑣 ) 2 ( 𝑤 ) ( 𝑣 ) 1 2 2 2 𝑍 𝑤𝛿 + 𝑍 𝛿 + 𝑍 𝑤𝛿 + 𝑍 𝛿 𝛿 + 𝑍 𝛿 + 𝑍 (38) 𝑜 𝑤 𝑜 𝑜 𝑤 𝑜 𝑤 𝑤 𝑜 𝑜 1 2 𝑤 2 1 2 2 𝑤 ( 𝑤 𝛿𝑤 ) 1 ( 𝑤 𝛿𝑤 ) ( 𝛿𝑤 𝛿𝑤 ) 2 ( 𝑜 ) ( 𝛿 ) ( 𝛿 ) 1 2 1 2 𝑤 𝑤 1 2 For use of the final model equations to predict the aero-propulsive forces and moments as a function of all static and dynamic explanatory variables, each of the models identified in Part II are used to compute the values of the parameters in Eqs. (31) - (36) as a function of the centered static explanatory variables ( 𝑣 , 𝑤 , 𝛿 , 𝛿 ). Then, the total 𝑤 𝑤 1 2 aero-propulsive forces and moments are computed using the Eqs. (31) - (36) polynomial expressions as a function of the centered dynamic explanatory variables ( 𝑛 , 𝑛 , ..., 𝑛 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 , 𝛿 ).

1 2 8 𝑒 𝑒 𝑒 𝑒 𝑓 𝑓 𝑓 𝑓 𝑟 𝑟 1 2 3 4 1 2 3 4 1 2 Several models were developed to explore the utility of using a different number of modeling blocks and different model complexities. The same general modeling approach was used for each model, with the only difference being the data volume associated with the sequential blocks of static test points described in Sec. III.A. Figure 17 shows the modeling NRMSE ( NRMSE ), the validation NRMSE ( NRMSE ), and the cumulative execution time for each block 𝑚 𝑣 (including the data collection time from all preceding blocks) with a different model complexity determined based on the analysis presented in Sec. III.A and with the knowledge that at least quadratic model, or up to pure quadratic and two-factor interaction model terms, are needed to model the LA-8 aircraft over the static explanatory variable range tested for this study [ 16 ]. Consequently, up to a quadratic model is used for Block 1, even though the model evaluation metrics do not meet the FDS criteria explained in Sec. III.A. Blocks 2-4 are assigned the highest model complexity meeting the FDS criteria for this comparison (see Table 1 and the accompanying discussion). The test execution time shown in Fig. 17b indicates the amount of time needed to complete the data collection with a sample time of 40 seconds per set point, including the amount of time to conduct a static tare run. Block 1 generally has the highest validation NRMSE value for each response; there is also a significant increase in the values of the validation NRMSE compared to the modeling NRMSE for the Block 1 model, suggesting that the model is deficient due to the sparsity of data points.

For the 𝑋 , 𝑌 , and 𝑍 responses, the Block 4 model has the lowest prediction error. For 𝐿 and 𝑀 , the Block 2 and Block 3 models, respectively, have the lowest prediction error. The prediction error for 𝑁 is similar for each number of modeling blocks. For the Block 2 to 4 models, the respective modeling NRMSE and validation NRMSE for each response have similar values, which provides confidence that model development was successful.

(a) Modeling and validation NRMSE (b) Cumulative test execution time Fig. 17 Comparison of NRMSE metrics and test time for each modeling block with different model complexity.

Balancing prediction performance, test execution time, and model parsimony, using data up to Block 2 (a cubic design order) with a quadratic model complexity was selected as the final modeling strategy. Given these results, it is recommended to use a cubic I-optimal experiment design with up to a quadratic polynomial structure for modeling similar aircraft over a similar range of test variables. For the Block 2 quadratic model, the prediction error, quantified by the NRMSE metric, is less than 6% for every response, and 𝑋 , 𝑌 , 𝑍 , 𝑀 , and 𝑁 have NRMSE values of approximately 𝑣 𝑣 3% or less, indicating that high-quality models have been developed.

Figure 18a shows the model fit and model prediction compared to the smoothed, measured aero-propulsive force and moment histories. Figure 18b shows the same plot zoomed in on a single validation set point, which better conveys the dynamic prediction capability. The model fits and model predictions are close to the measured responses, indicating that the model is able to describe a large amount of the variation in each response. Figure 19 shows a history of normalized modeling residuals and normalized validation residuals. The modeling and prediction residuals have similar magnitudes, supporting the claim that good predictive models have been identified.

(a) Modeling and validation data (two modeling blocks) (b) Sample validation set point Fig. 18 Model fit and model prediction compared to smoothed, measured response data.

Fig. 19 Normalized modeling and prediction residuals against time.

VII. Discussion The preceding results show that good models have been identified for the LA-8 eVTOL aircraft using the hybrid testing and compound modeling approach developed in this paper. Furthermore, the final aero-propulsive model at ¯ 𝑞 = 3 . 5 psf was developed from a total of 48 minutes of test execution time. If the test technique and modeling approach were to be applied to model the LA-8 transition envelope tested in Refs. [ 15 , 16 ] by executing testing and developing models at several dynamic pressure settings, then a global transition model could be developed in a single day of wind tunnel testing. This approach is significantly faster than previous testing conducted to characterize the LA-8 aircraft using purely static DOE/RSM testing, which was already vastly more efficient compared to using one-factor-at-a-time testing [ 15 ]. The amount of test time needed for modeling data collection and tare runs for a purely static DOE/RSM wind tunnel experiment to develop aero-propulsive models at ¯ 𝑞 = 3 . 5 psf described in Ref. [ 16 ] took approximately 267 minutes (4.45 hours). The hybrid DOE/RSM and PTI wind tunnel test technique described in this paper is 5.6 times more efficient in terms of test execution time for the data used to develop the model. Note that the models developed in Ref. [ 16 ] included additional quadratic and interaction model terms for the control effectors; however, that is largely because of the reduced range of excitation for those variables in this study. The control effector amplitude used for the results presented in this paper could be increased without requiring additional test time.

In addition to reduced test time, the models developed using the approach presented in this paper allow characterization of additional important aero-propulsive phenomena for tilt-wing aircraft. The tilting wings add significant modeling challenges not seen in conventional aircraft designs because the propellers, wings, and wing-fixed control surfaces all change orientation with respect to the modeling frame of reference in the body axes. Furthermore, there will also be interactions with the magnitude and orientation of oncoming airflow, suggesting that up to four-factor interactions may be present among the explanatory variables. One way of handling this complexity is to develop a different aero-propulsive model at each combination of wing angle settings, treating each different combination of wing angle orientation as a vehicle configuration change. This method would be ideal when a transition wing angle schedule has been defined; however, this is impractical when the identified aero-propulsive model is used to develop the flight controller and informs the transition wing angle schedule. In previous LA-8 modeling work [ 16 ], the propeller, control surface, wing, and airflow interactions were modeled using only two-factor interaction effects and were limited by the purely static DOE/RSM testing strategy. Alternatively, the experimental design and compound modeling strategy developed in this paper supports modeling additional interactions among the explanatory variables. This is possible because of the way the model parameters associated with the dynamic test variables are modeled as a function of the static test variables.

As an example of how the compound modeling approach supports characterization of additional complexity, consider Eq. (37) , where 𝑍 is modeled as a function of 𝑤 and 𝛿 including variation with a 𝑤𝛿 interaction and a 𝛿 𝛿 𝑤 𝑤 𝑓 1 1 𝑤 quadratic term. Note that these are effectively cubic model terms because of the multiplication of the 𝑍 parameter by 𝛿 𝑓 𝛿 . For a static DOE/RSM experiment supporting up to a quadratic model, the only model terms involving 𝑤 and 𝑓 𝛿 that could model an interaction with 𝛿 are 𝛿 𝑤 and 𝛿 𝛿 . Suppose also that a larger range of the dynamic 𝑤 𝑓 𝑓 𝑓 𝑤 1 1 1 1 1 test variables were tested and the interaction term 𝑛 𝛿 was included in the Part I model. Then, by the same process, 2 𝑓 the 𝑛 𝛿 interaction could be modeled as a function of 𝑤 and 𝛿 including nonlinear model terms. This additional 2 𝑓 𝑤 1 1 parameterization flexibility enabled by the testing and modeling approach described herein provides significant benefits to modeling complex eVTOL aircraft. Furthermore, the range of each static explanatory variable, determined as described in Ref. [ 15 ], could be increased to model a larger range of aero-propulsive phenomena. This might require additional design complexity beyond the Block 2 cubic experiment design order, but this increase in static experiment design complexity could be accommodated with only modest increases in the amount of test time (cf. Fig. 17b) while allowing characterization of a broader range of flight conditions and aero-propulsive complexity using the same model identification framework.

When using the modeling approach developed in this work, a few limitations and drawbacks should be considered.

First, the models are based on a quasi-steady assumption, where aero-propulsive effects at the current point in time are only dependent on the current states and controls. Second, the data were collected in the wind tunnel with zero vehicle angular velocity, and consequently, no damping terms could be identified. Identification of models dependent on vehicle angular rates and the history of the explanatory variables will be needed to improve model predictive capability in dynamic maneuvering. Also, integrating the PTI excitation technique into an existing wind tunnel test facility typically requires significant time and engineering effort. After the PTI excitation capability is integrated, however, the test efficiency gains are significant. Finally, the compound model identification approach is more complex than ordinary least-squares regression analysis, which is generally used for model development from conventional static DOE/RSM experiments and is readily available in commercial statistical software packages. However, the model identification † approach described in this paper is still tractable to implement using algorithms available in SIDPAC to develop

models leveraging the much more efficient hybrid test technique. Ultimately, the ability to develop a higher complexity aero-propulsive model and significantly reduce the amount of required test time can easily justify the steps required to implement the new approach.

VIII. Conclusions

A novel wind tunnel testing and model identification strategy was developed and applied for eVTOL aircraft aero-propulsive modeling to support flight dynamics simulation development. The LA-8 tandem tilt-wing eVTOL aircraft that was used for this exploratory study exhibits aero-propulsive complexity representative of many current and future eVTOL aircraft. A wind tunnel experiment was developed in two parts. A four-factor static experiment was designed in a set of test blocks allowing identification of response surface models of increasing complexity for slow-moving airflow and wing angle test variables. A dynamic experiment composed of 18 simultaneous orthogonal phase-optimized multisine signals was designed for the propulsors and control surfaces, and executed at each set point.

A cubic I-optimal response surface design for the static test factors executed simultaneously with dynamic excitations running for 40 seconds at each static test point was found to be a good testing strategy for model development, balancing test time and model prediction capability. Using this approach, all data collection at a reference airspeed condition can be completed in under 48 minutes, which is significantly faster than previous eVTOL aircraft modeling efforts using purely static test techniques.

The model identification approach consisted of two parts in accordance with the data collection strategy. First, a local model was identified at each set point as a function of the dynamically changing explanatory variables using multivariate orthogonal function modeling and ordinary least-squares regression in the frequency domain. The model structure for each force and moment component was chosen to be identical for each set point based on the model terms deemed significant for a majority of the set points. Second, response surface models were developed for each model parameter associated with the dynamic test variables identified in the first stage of modeling as a function of the static test variables, using multivariate orthogonal function modeling and weighted least-squares regression leveraging available uncertainty estimates. The identified models were shown to have good predictive capability and small normalized model fit error. In addition, the model form supports parameterization of nonlinear aero-propulsive effects that cannot be captured in alternative modeling approaches.

Accurate and efficient modeling of eVTOL aircraft with inherently complex aero-propulsive characteristics is essential for effective flight control system design, handling qualities evaluation, and realistic flight simulation. These activities are required to enable safe and effective eVTOL aircraft flight operations. Modern eVTOL aircraft designs have many control effectors and complex aero-propulsive coupling, which require both efficient testing and general nonlinear modeling capability. The wind tunnel testing and modeling approach described in this work addresses these requirements.

Acknowledgments

This research was funded by the NASA Aeronautics Research Mission Directorate (ARMD) Transformational Tools and Technologies (TTT) project. Photography support was provided by Lee Pollard. Wind tunnel testing and LA-8 vehicle support was provided by Gregory Howland, David North, and Steven Geuther.

References

[1] Antcliff, K. R., Whiteside, S. K. S., Kohlman, L. W., and Silva, C., “Baseline Assumptions and Future Research Areas for Urban Air Mobility Vehicles,” AIAA SciTech 2019 Forum , AIAA Paper 2019-0528, Jan. 2019. https://doi.org/10.2514/6.2019-0528.

[2] Patterson, M. D., Antcliff, K. R., and Kohlman, L. W., “A Proposed Approach to Studying Urban Air Mobility Missions Including an Initial Exploration of Mission Requirements,” AHS International 74th Annual Forum & Technology Display , May 2018.

[3] Johnson, W., Silva, C., and Solis, E., “Concept Vehicles for VTOL Air Taxi Operations,” AHS Technical Conference on Aeromechanics Design for Transformative Vertical Flight , Jan. 2018.

[4] Silva, C., Johnson, W., Antcliff, K. R., and Patterson, M. D., “VTOL Urban Air Mobility Concept Vehicles for Technology Development,” 2018 Aviation Technology, Integration, and Operations Conference , AIAA Paper 2018-3847, Jun. 2018.

https://doi.org/10.2514/6.2018-3847.

[5] Johnson, W., and Silva, C., “NASA Concept Vehicles and the Engineering of Advanced Air Mobility Aircraft,” The Aeronautical Journal , Vol. 126, No. 1295, 2022, pp. 59–91. https://doi.org/10.1017/aer.2021.92.

[6] Rothhaar, P. M., Murphy, P. C., Bacon, B. J., Gregory, I. M., Grauer, J. A., Busan, R. C., and Croom, M. A., “NASA Langley Distributed Propulsion VTOL Tilt-Wing Aircraft Testing, Modeling, Simulation, Control, and Flight Test Development,” 14th AIAA Aviation Technology, Integration, and Operations Conference , AIAA Paper 2014-2999, Jun. 2014. https: //doi.org/10.2514/6.2014-2999.

[7] Saeed, A. S., Younes, A. B., Cai, C., and Cai, G., “A survey of hybrid Unmanned Aerial Vehicles,” Progress in Aerospace Sciences , Vol. 98, 2018, pp. 91–105. https://doi.org/10.1016/j.paerosci.2018.03.007.

[8] Kim, H. D., Perry, A. T., and Ansell, P. J., “A Review of Distributed Electric Propulsion Concepts for Air Vehicle Technology,” 2018 AIAA/IEEE Electric Aircraft Technologies Symposium , AIAA Paper 2018-4998, Jul. 2018. https://doi.org/10.2514/6.2018- 4998.

[9] North, D. D., Busan, R. C., and Howland, G., “Design and Fabrication of the Langley Aerodrome No. 8 Distributed Electric Propulsion VTOL Testbed,” AIAA SciTech 2021 Forum , AIAA Paper 2021-1188, Jan. 2021. https://doi.org/10.2514/6.2021-1188.

[10] Murphy, P. C., Buning, P. G., and Simmons, B. M., “Rapid Aero Modeling for Urban Air Mobility Aircraft in Computational Experiments,” AIAA SciTech 2021 Forum , AIAA Paper 2021-1002, Jan. 2021. https://doi.org/10.2514/6.2021-1002.

[11] Simmons, B. M., Buning, P. G., and Murphy, P. C., “Full-Envelope Aero-Propulsive Model Identification for Lift+Cruise Aircraft Using Computational Experiments,” AIAA AVIATION 2021 Forum , AIAA Paper 2021-3170, Aug. 2021. https: //doi.org/10.2514/6.2021-3170.

[12] Busan, R. C., Rothhaar, P. M., Croom, M. A., Murphy, P. C., Grafton, S. B., and O’Neal, A. W., “Enabling Advanced Wind-Tunnel Research Methods Using the NASA Langley 12-Foot Low Speed Tunnel,” 14th AIAA Aviation Technology, Integration, and Operations Conference , AIAA Paper 2014-3000, Jun. 2014. https://doi.org/10.2514/6.2014-3000.

[13] Murphy, P. C., and Landman, D., “Experiment Design for Complex VTOL Aircraft with Distributed Propulsion and Tilt Wing,” AIAA Atmospheric Flight Mechanics Conference , AIAA Paper 2015-0017, Jan. 2015. https://doi.org/10.2514/6.2015-0017.

[14] Murphy, P. C., Simmons, B. M., Hatke, D. B., and Busan, R. C., “Rapid Aero Modeling for Urban Air Mobility Aircraft in Wind-Tunnel Tests,” AIAA SciTech 2021 Forum , AIAA Paper 2021-1644, Jan. 2021. https://doi.org/10.2514/6.2021-1644.

[15] Busan, R. C., Murphy, P. C., Hatke, D. B., and Simmons, B. M., “Wind Tunnel Testing Techniques for a Tandem Tilt- Wing, Distributed Electric Propulsion VTOL Aircraft,” AIAA SciTech 2021 Forum , AIAA Paper 2021-1189, Jan. 2021.

https://doi.org/10.2514/6.2021-1189.

[16] Simmons, B. M., and Murphy, P. C., “Aero-Propulsive Modeling for Tilt-Wing, Distributed Propulsion Aircraft Using Wind Tunnel Data,” Journal of Aircraft , published online 2 March 2022. https://doi.org/10.2514/1.C036351.

[17] Simmons, B. M., “System Identification for eVTOL Aircraft Using Simulated Flight Data,” AIAA SciTech 2022 Forum , AIAA Paper 2022-2409, Jan. 2022. https://doi.org/10.2514/6.2022-2409.

nd [18] Morelli, E. A., and Klein, V., Aircraft System Identification: Theory and Practice , 2 ed., Sunflyte Enterprises, Williamsburg, VA, 2016.

[19] Tischler, M. B., and Remple, R. K., Aircraft and Rotorcraft System Identification: Engineering Methods with Flight Test nd Examples , 2 ed., American Institute of Aeronautics and Astronautics, Reston, VA, 2012. https://doi.org/10.2514/4.868207.

nd [20] Jategaonkar, R. V., Flight Vehicle System Identification: A Time-Domain Methodology , 2 ed., American Institute of Aeronautics and Astronautics, Reston, VA, 2015. https://doi.org/10.2514/4.102790.

th [21] Montgomery, D. C., Design And Analysis of Experiments , 8 ed., John Wiley & Sons, Inc., Hoboken, NJ, 2013.

[22] Myers, R. H., Montgomery, D. C., and Anderson-Cook, C. M., Response Surface Methodology: Process and Product th Optimization Using Designed Experiments , 4 ed., John Wiley & Sons, Hoboken, NJ, 2016.

[23] Morelli, E. A., “Multiple Input Design for Real-Time Parameter Estimation in the Frequency Domain,” 13th IFAC Conference on System Identification , Aug. 2003. https://doi.org/10.1016/S1474-6670(17)34833-4.

[24] Morelli, E. A., “Flight-Test Experiment Design for Characterizing Stability and Control of Hypersonic Vehicles,” Journal of Guidance, Control, and Dynamics , Vol. 32, No. 3, 2009, pp. 949–959. https://doi.org/10.2514/1.37092.

[25] Morelli, E. A., “Practical Aspects of Multiple-Input Design for Aircraft System Identification Flight Tests,” AIAA AVIATION 2021 Forum , AIAA Paper 2021-2795, Aug. 2021. https://doi.org/10.2514/6.2021-2795.

[26] Favaregh, N., and Landman, D., “Global Modeling of Pitch Damping from Flight Data,” AIAA Atmospheric Flight Mechanics Conference and Exhibit , AIAA Paper 2006-6145, Aug. 2006. https://doi.org/10.2514/6.2006-6145.

[27] Omran, A., Newman, B., and Landman, D., “Global Aircraft Aero-Propulsive Linear Parameter-Varying Model Using Design of Experiments,” Aerospace Science and Technology , Vol. 22, No. 1, 2012, pp. 31–44. https://doi.org/10.1016/j.ast.2011.05.008.

[28] Heim, E. H., Viken, E., Brandon, J. M., and Croom, M. A., “NASA’s Learn-to-Fly Project Overview,” 2018 Atmospheric Flight Mechanics Conference , AIAA Paper 2018-3307, Jun. 2018. https://doi.org/10.2514/6.2018-3307.

[29] Riddick, S. E., “An Overview of NASA’s Learn-to-Fly Technology Development,” AIAA SciTech 2020 Forum , AIAA Paper 2020-0760, Jan. 2020. https://doi.org/10.2514/6.2020-0760.

[30] Murphy, P. C., and Brandon, J., “Efficient Testing Combining Design of Experiment and Learn-to-Fly Strategies,” AIAA Atmospheric Flight Mechanics Conference , AIAA Paper 2017-0696, Jan. 2017. https://doi.org/10.2514/6.2017-0696.

[31] Morelli, E. A., “Real-Time Global Nonlinear Aerodynamic Modeling for Learn-To-Fly,” AIAA Atmospheric Flight Mechanics Conference , AIAA Paper 2016-2010, Jan. 2016. https://doi.org/10.2514/6.2016-2010.

[32] Morelli, E. A., “Practical Aspects of Real-Time Modeling for the Learn-To-Fly Concept,” 2018 Atmospheric Flight Mechanics Conference , AIAA Paper 2018-3309, Jun. 2018. https://doi.org/10.2514/6.2018-3309.

[33] Morelli, E. A., “Autonomous Real-Time Global Aerodynamic Modeling in the Frequency Domain,” AIAA SciTech 2020 Forum , AIAA Paper 2020-0761, Jan. 2020. https://doi.org/10.2514/6.2020-0761.

[34] McSwain, R. G., Geuther, S. C., Howland, G., Patterson, M. D., Whiteside, S. K., and North, D. D., “An Experimental Approach to a Rapid Propulsion and Aeronautics Concepts Testbed,” NASA TM–2020-220437, Jan. 2020.

[35] Geuther, S. C., and Fei, X., “LA-8 Computational Analysis and Validation Studies Using FlightStream,” AIAA SciTech 2021 Forum , AIAA Paper 2021-1191, Jan. 2021. https://doi.org/10.2514/6.2021-1191.

[36] Geuther, S. C., North, D. D., and Busan, R. C., “Investigation of a Tandem Tilt-wing VTOL Aircraft in the NASA Langley 12-Foot Low-Speed Tunnel,” NASA TM–2020–5003178, Jun. 2020.

[37] Simmons, B. M., and Hatke, D. B., “Investigation of High Incidence Angle Propeller Aerodynamics for Subscale eVTOL Aircraft,” NASA TM–20210014010, May 2021.

[38] Simmons, B. M., “System Identification for Propellers at High Incidence Angles,” Journal of Aircraft , Vol. 58, No. 6, 2021, pp.

1336–1350. https://doi.org/10.2514/1.C036329.

[39] Cook, J., and Gregory, I., “A Robust Uniform Control Approach for VTOL Aircraft,” VFS Autonomous VTOL Technical Meeting and Electric VTOL Symposium , Jan. 2021.

[40] North, D. D., “Flight Testing of a Scale Urban Air Mobility Technology Testbed,” AIAA SciTech 2021 Forum , AIAA Presentation, Jan. 2021. URL https://ntrs.nasa.gov/citations/20205010998, accessed April 2022.

[41] Zahran, A., Anderson-Cook, C. M., and Myers, R. H., “Fraction of Design Space to Assess Prediction Capability of Response Surface Designs,” Journal of Quality Technology , Vol. 35, No. 4, 2003, pp. 377–386. https://doi.org/10.1080/00224065.2003.

11980235.

[42] Anderson, M. J., and Whitcomb, P. J., RSM Simplified: Optimizing Processes Using Response Surface Methods for Design of nd Experiments , 2 ed., CRC Press, Boca Raton, FL, 2017. https://doi.org/10.1201/9781315382326.

[43] Anderson, M. J., Adams, W. F., and Whitcomb, P. J., “How to Properly Size Response Surface Method Experiment (RSM) Designs for System Optimization,” 2016. URL https://cdn.statease.com/media/public/documents/how_to_properly_size_ experiments_aimed_at_system_optimization.pdf, [retrieved April 2022].

[44] Whitcomb, P., “FDS—A Power Tool for Designers of Optimization Experiments,” Stat-Teaser Newsletter , Stat-Ease, Inc., Sep.

2008, pp. 1–3. URL https://cdn.statease.com/media/public/documents/statteaser-0908.pdf, [retrieved April 2022].

[45] Morelli, E. A., and Smith, M. S., “Real-Time Dynamic Modeling: Data Information Requirements and Flight-Test Results,” Journal of Aircraft , Vol. 46, No. 6, 2009, pp. 1894–1905. https://doi.org/10.2514/1.40764.

[46] Morelli, E. A., “Flight Test Maneuvers for Efficient Aerodynamic Modeling,” Journal of Aircraft , Vol. 49, No. 6, 2012, pp.

1857–1867. https://doi.org/10.2514/1.C031699.

[47] Grauer, J. A., and Boucher, M. J., “Identification of Aeroelastic Models for the X-56A Longitudinal Dynamics Using Multisine Inputs and Output Error in the Frequency Domain,” Aerospace , Vol. 6, No. 2, 2019, pp. 1–25. https://doi.org/10.3390/ aerospace6020024.

[48] Grauer, J. A., and Boucher, M. J., “System Identification of Flexible Aircraft: Lessons Learned from the X-56A Phase 1 Flight Tests,” AIAA SciTech 2020 Forum , AIAA Paper 2020-1017, Jan. 2020. https://doi.org/10.2514/6.2020-1017.

th [49] Montgomery, D. C., Peck, E. A., and Vining, G. G., Introduction to Linear Regression Analysis , 5 ed., John Wiley & Sons, Hoboken, New Jersey, 2012.

[50] Belsley, D. A., Kuh, E., and Welsch, R. E., Regression Diagnostics: Identifying Influential Data and Sources of Collinearity , John Wiley & Sons, Hoboken, NJ, 2004.

[51] Gustafsson, F., “Determining the Initial States in Forward-Backward Filtering,” IEEE Transactions on Signal Processing , Vol. 44, No. 4, 1996, pp. 988–992. https://doi.org/10.1109/78.492552.

[52] Morelli, E. A., and Grauer, J. A., “Practical Aspects of Frequency-Domain Approaches for Aircraft System Identification,” Journal of Aircraft , Vol. 57, No. 2, 2020, pp. 268–291. https://doi.org/10.2514/1.C035599.

[53] Morelli, E. A., “Practical Aspects of the Equation-Error Method for Aircraft Parameter Estimation,” AIAA Atmospheric Flight Mechanics Conference and Exhibit , AIAA Paper 2006-6144, Aug. 2006. https://doi.org/10.2514/6.2006-6144.

[54] Klein, V., “Aircraft Parameter Estimation in Frequency Domain,” 4th Atmospheric Flight Mechanics Conference , AIAA Paper 1978-1344, Aug. 1978. https://doi.org/10.2514/6.1978-1344.

[55] Morelli, E. A., “High Accuracy Evaluation of the Finite Fourier Transform Using Sampled Data,” NASA TM-110340, Jun.

1997.

[56] Morelli, E. A., “Global Nonlinear Aerodynamic Modeling Using Multivariate Orthogonal Functions,” Journal of Aircraft , Vol. 32, No. 2, 1995, pp. 270–277. https://doi.org/10.2514/3.46712.

[57] Morelli, E. A., “Transfer Function Identification using Orthogonal Fourier Transform Modeling Functions,” AIAA Atmospheric Flight Mechanics (AFM) Conference , AIAA Paper 2013-4749, Aug. 2013. https://doi.org/10.2514/6.2013-4749.

[58] Morelli, E. A., and Cooper, J., “Frequency-Domain Method for Automated Simulation Updates Based on Flight Data,” Journal of Aircraft , Vol. 52, No. 6, 2015, pp. 1995–2008. https://doi.org/10.2514/1.C033121.

[59] Barron, A. R., “Predicted Squared Error: A Criterion for Automatic Model Selection,” Self-Organizing Methods in Modeling , Farlow, S. J., Ed., Marcel Dekker, Inc., New York, NY, 1984, pp. 87–104.

[60] Morelli, E. A., and Klein, V., “Accuracy of Aerodynamic Model Parameters Estimated from Flight Test Data,” Journal of Guidance, Control, and Dynamics , Vol. 20, No. 1, 1997, pp. 74–80. https://doi.org/10.2514/2.3997.

[61] Morelli, E. A., and DeLoach, R., “Response Surface Modeling using Multivariate Orthogonal Functions,” 39th AIAA Aerospace Sciences Meeting and Exhibit , AIAA Paper 2001-0168, Jan. 2001. https://doi.org/10.2514/6.2001-168.

[62] Morelli, E. A., and DeLoach, R., “Wind Tunnel Database Development Using Modern Experiment Design and Multivariate Orthogonal Functions,” 41st Aerospace Sciences Meeting and Exhibit , AIAA Paper 2003-653, Jan. 2003. https://doi.org/10.

2514/6.2003-653.