Parameterization of nonlinear aeroelastic reduced order models via direct interpolation of Taylor partial derivatives

The identiﬁcation of optimally sparse Taylor partial derivatives presents a new opportunity in eﬃcient nonlinear aerodynamic model reduction for complex aeroelastic systems. Unfortunately, for this class of reduced order model (ROM), the robustness that is observed in the linear regime to parameters including; dynamic pressure, control hinge linear stiﬀness, or even freeplay, is quickly forfeited given the amplitude (or velocity) dependency of the aerodynamic loading on the structure. In this paper, nonlinear sensitivity is addressed by interpolating a library of nonlinear unsteady aerodynamic ROMs across a compact subspace in dynamic pressure and freeplay magnitude.

The ROM, based on Lagrange interpolation of sparse higher-order Taylor partial derivatives, demonstrates excellent precision in modelling high amplitude

Introduction
One of the challenges in modeling aeroelastic systems with discrete structural nonlinearities such as, freeplay, in the transonic flow regime, is the potential for nonlinear aerodynamic loads that coexist with the structural nonlinearity [1], [2].While linearization of the nonlinear structural model is common practice in aeroelastic problems involving freeplay (employing techniques such as fictitious masses [3] or residual vectors [4]), the assumption of aerodynamic linearization may, for some transonic systems, be invalid with small parameter changes [5].In such cases, to avoid the exhaustive computational burden associated with using a computational fluid dynamics (CFD) code to resolve the forces at each time-step, a vast range of nonlinear unsteady aerodynamic reduced order models (ROMs) can be used [6], [7], [8], [9], [10], [11].
Considering that in the field of nonlinear aeroelasticity the nonlinear unsteady aerodynamic forces on the elastic structure are typically reduced to integrated quantities, the functional series approach (multi-variable Taylor series expansion or Volterra series) is well suited, describing the generalized nonlinear aerodynamic forces on the structure as a nonlinear dynamic function of the generalized structural displacements.The merits of this class of ROM are that it is relatively simplistic to implement with minimal modification to standard CFD software, it intrinsically captures nonlinearity in its functional form, and that the online computational savings are typically several order of magnitude compared to full-order aeroelastic simulations.Two primary drawbacks include the curse of dimensionality that is associated with identifying higher-order systems -traditionally limiting it from being applied to complex 3D aeroelastic systems, and a reduced ability to generalize for systems with strong nonlinearity.
In recent work by the authors [2], limitations related to the curse-ofdimensionality are addressed using sparsity promoting algorithms to signif-icantly reduce the amount of training data required to identify the higherorder polynomial terms.The study demonstrates that through Orthogonal Matching Pursuit, it is possible to efficiently identify the optimized s-sparse nonlinear unsteady aerodynamic ROMs, exemplified through application to a three-dimensional aeroelastic stabilator model experiencing high amplitude freeplay-induced limit cycles.The comparison shows excellent agreement between the ROM and the full-order model (FOM), demonstrating the feasibility of applying accurate, higher-order polynomial-based ROMs to complex nonlinear aeroelastic problems without incurring significant computational burdens.
In terms of the second limitation related to increased sensitivity and reduced ability to generalize, first one must consider the importance of generalization and what a reasonable expectation is for a given ROM to generalize.For example, consider the linear definition of functional series-based ROMs, i.e., the aerodynamic impulse response which, within the linear unsteady aerodynamic regime, is quite robust to certain parameter changes.These include; dynamic pressure (the generalized forces scale linearly for a fixed Mach number and angle-of-attack (AoA)), structural stiffness and damping parameters (e.g., related to control surfaces or pylon linkages), and even nonlinear structural parameters, such as, freeplay magnitude.As a result of this insensitivity, it is common practice to generate a linear ROM at a fixed Mach number and AoA (noting high sensitivity to these two parameters), and to then use the linearized ROM freely for the desired application.Unfortunately, in the nonlinear regime (e.g., high amplitude limit cycles and/or high-AoA maneuvers), there is no guarantee of such generalizability given the structural amplitude (or velocity) dependency of the nonlinear aerodynamic loading on the structure.Specifically, as the nonlinearity strength increases, the ROM can become sensitive to said parameters, hereby leading to a reduction in the ability of the ROM to generalize.With this in mind, it is desirable to generate a single nonlinear ROM that functions well within the aforementioned parametric sub-spaces, i.e., the objective of this research.This paper proposes that using Lagrange polynomials, it is possible to interpolate a library of linear and sparse nonlinear unsteady aerodynamic ROMs (in the form of Taylor partial derivatives), across a subspace defined in dynamic pressure and freeplay.The library contains unsteady aerodynamic ROMs that are tuned for optimal aeroelastic performance for various combinations of freeplay and dynamic pressure, referred to as sampled ROMs from now on.The interpolation scheme allows the sampled ROMs to be generated with disparity in polynomial-order, sparsity, and cardinality, i.e., the shape and sparsity of the derivative tensors to be interpolated are unimportant.The interpolated nonlinear unsteady aerodynamic ROM is applied to a three-dimensional aeroelastic stabilator model with freeplay undergoing high amplitude LCO at high transonic Mach numbers, demonstrating excellent precision with respect to the full-order aeroelastic model.
The remainder of the paper is constructed as follows; in Section 2 the ROM identification procedure and interpolation scheme are defined.The nonlinear aeroelastic system is described in Section 3 and practical details pertaining to implementation of the nonlinear unsteady aerodynamic ROM library and interpolation scheme in Section 4. The results are presented and discussed in Section 5, then final discussion and concluding remarks are given in Section 6.
2 Generation and Interpolation of the Unsteady Aerodynamic ROM Library

Unsteady Aerodynamic ROMs Based on Functional Series
The earliest works in nonlinear aeroelasticity consider the identification of Volterra kernels by applying unit-impulse/step functions (or variations of these) to excite the structural system within a full-order aerodynamic solver, and recording the "exact" linear / nonlinear aerodynamic impulse response [12], [13], [14], [15], [16], [17], [18], [19].More recently, it has been shown that identification of the nonlinear kernels using impulses can be overly rigid and comes with a range of limitations, in particular for CFD-based identification [20].An alternative approach is to randomly excite the full-order system and record the unsteady aerodynamic response.Then using any form of linear regression, the linear and nonlinear kernels can be derived.Examples include; time-delay neural networks [21], least squares (LS) [22], [20], or sparsity promoting LS algorithms [2], where the limitations described above can be largely overcome.To this effect, Balajewicz and Dowell [20] introduce sparse Volterra series for efficient generation of higher-order order transonic aeroelastic ROMs, which was recently extended by the authors [2] using sparsity promoting algorithms to drastically reduce the amount of training data needed to identify the higherorder system.Brown [23] presents a multi-input Volterra based approach for a 3-DOF aeroelastic system.A blended step is used as excitation and ℓ 1regularized least-squares is used to derive the kernels which are expressed in terms of Laguerre polynomials to reduce the number of coefficients to be identified.The performance is generally very good, however, it is noted that out-of-sample performance in terms of reduced velocity is challenging in the nonlinear regime, as to be expected.

Multi-Variable Taylor Series Expansion of the Unsteady Aerodynamic Forces
Assuming that the unsteady aerodynamic forces on the structure are a dynamic function of the structural displacement, i.e., there is no self-sustaining unsteadiness in the flow, the relationship can be written as where Q represents the aerodynamic forces, u represents the structural displacements and f () is some unknown function.This can be rewritten in discrete time, assuming that the relationship is time-invariant and asymptotically stable, as where u = {u n , u n−1 , u n−2 , . . .u n−k } T is a scalar valued function length k, with the subscript denoting the discrete time interval.Assuming that the system is mildly nonlinear and memory fading then k defines the number of time lags required for the nonlinear memory to fade and f (u) can be approximated using a multi-variable Taylor series expansion, given by which can be evaluated at u = a and written in multi-index form as where d p = (∂ p f )(a) is a p th -order tensor that contains the p th -order partial derivatives of f (a), i.e., d 1 is the gradient of f (a), d 2 is the Hessian matrix, and so on.Given that f (a) is not known a priori, the terms of the tensors d 0 , d 1 , . . ., d p are estimated from input-output training data, i.e., the coefficients are identified to minimize the error between the Taylor approximation and the true values by argmin||Q − T (u)||.

Identification of the Taylor Partial Derivatives
The procedure for identifying the unsteady unsteady aerodynamic ROMs is now described.This definition is for a single-input single-output system.The reader is referred to [20] for identification using a multi-input approach.

Band Limited Random Excitation
The first step in creating any unsteady aerodynamic ROM is to perturb the structural degrees-of-freedom (DOFs) within a full-order aerodynamic solver and to record the aerodynamic responses.In this approach, the structural model is excited using band limited random noise within a commercial finite volume CFD solver.Using the SISO identification procedure, each structural DOF is excited in isolation (i.e., nonlinear interactions between structural DOFs are neglected).The amplitude and frequency band of the excitation functions are chosen by estimating the LCO amplitudes and frequencies in the aeroelastic response.An important consideration is to observe a smooth transition from the undeformed structure and converged steady-state fluid forces.The reason being that any discontinuity (i.e., a step-like change in structural displacement) will cause spurious aerodynamic response information and inaccuracies in the ROM identification process.A hyperbolic tangent is applied to the the raw signal u j,r , ensuring a smooth transition to the modal excitation over the first 20 time intervals, given for n total input samples by where u j is the band limited random signal used to excite the j th structural DOF.

Input and Output Matrices
Identifying the coefficients of the derivatives is a linear problem [24], given according to the system inputs and outputs by where M is a matrix of inputs, Q is a vector of outputs, and D = {d 1 , d 2 , . . ., d p } contains the p th -order Taylor partial derivatives which are identified by solving the inverse problem.The construction of each of the inputs and outputs is now described.Considering a total of m structural DOFs, to construct the vector of outputs each j th structural DOF is excited individually using u j ∈ R n , j = 1, . . ., m (Eq.5) and the resultant full-order aerodynamic forces are projected onto each i th structural DOF To construct the matrix of inputs, first a lower left triangular circulant matrix is constructed from u j (truncated for k time lags) to give L j ∈ R n×k , then the p th -order multivariable Taylor series expansion of the rows gives M j ∈ R n×κ where κ = p 1 k+(p−1) p .

Identification of the Reduced Order Model using Pseudo-Inverses
The full set of partial derivatives (no sparsity) D ∈ R m×m can be identified using pseudo-inverses, evaluated according to where + is the Moore-Penrose Pseudo-Inverse.D ij ∈ R κ×1 contains the flattened derivative tensors corresponding to the generalised aerodynamic forces for Q i (u j ).The final full ROM is given by iterating through M j and Q ij , according to

Identification of the Reduced Order Model using Orthogonal Matching Pursuit
Initially, the OMP algorithm is explicitly defined.[25] Considering the linear problem for then the objective is to recover a sparse representation of x ∈ R nx , according to where ||x|| 0 is the ℓ 0 pseudo-norm, which is the number of non-zero elements in x.Assuming that x is s-sparse (s x ≥ ||x|| 0 ), it can be recovered exactly by OMP if A and x satisfy following inequality: where µ M is the mutual coherence of the columns of A and s x is the sparsity of x.From Eq. 11, x can be at most 1 2µ -sparse.The OMP algorithm is as follows To construct the Optimal Sparsity ROM (OS-ROM) Ds ∈ R m×m , M j and Q ij can be substituted into Eq.10, giving argmin where D ij s ∈ R κ×1 contains the optimized flattened sparse derivative tensors corresponding to the generalised aerodynamic forces for Q i (u j ).Eq. 12 can be solved directly using the OMP Algorithm 1 according to Finally, iterating through each structural DOF, the OS-ROM is given by Practically speaking, D s contains the optimal set of s-sparse coefficients of the p th -order unsteady aerodynamic ROM, satisfying the linear problem in Eq. 6 which can be used in a time-marching aerodynamic or aeroelastic simulation, to obtain the aerodynamic forces on the structure at each time interval.

Scheme for Direct Interpolation of the ROM Library
This inability of the nonlinear ROM identified at a single set of operating conditions to generalize to new operating conditions is a limitation that is addressed by interpolating a library of ROMs (Taylor partial derivatives).The objective is to minimize the number of sampled ROMs that need to identified (given the associated computational cost), while ensuring that the twodimensional subspace can be resolved with sufficient accuracy.
For an N -dimensional parameter space, containing the parameters χ 1 , χ 2 , . . ., χ N , for N R > 1 sampled ROMs ( Di T ) that are constructed at operating points X i = {χ 1i , χ 2i , . . ., χ N i }, then the set of operating points of these ROMs Λ is given by Then the objective is to interpolate the bases of ROMs to construct a new ROM at the operating point X NR+1 / ∈ X.Before direct interpolation is possible, disparity in the tensor dimensions must be accounted for by projecting each ROM onto a base vector DTI0 that contains the indices of all the coefficients contained in the ROM library.Given that each p th -order tensor in Di T can be generated with a different number of lag terms, the base vector is created by identifying the ROM with the maximum number of terms for each derivative order separately, given by where p , k mp is the number of lag terms used to generate the largest p th -order tensor in the basis, and therefore, DTI0 contains the index of every coefficient in the bases.Lagrange polynomials are used to interpolate the ROM basis which, considering that each parameter can be defined by an interpolation scheme of different order p L1 , p L2 , . . ., p LN , are given at the new location j by where where χ N k N and χ N l N are the operating conditions at the sampled ROM locations and χ N j is the operating conditions of the interpolated ROM.

Modified AGARD 445.6 Wing
The AGARD 445.6 wing is a well known transonic benchmark case, with experiments conducted in the NASA transonic dynamic wind tunnel.The model consists of a tapered swept wing (see Figure 1 In this paper, the wing is modified to represent an all-movable control surface (first presented by Carrese et al. [29]), with a torsional spring added to a centered node at the root, which is free to rotate about the pitch axis.The torsional spring contains a zero-stiffness dead-zone and a nominal stiffness of k δ = 500 Nm/rad otherwise, as depicted in Fig. 1(b).

Nonlinear Aeroelastic Equations-of-Motion
The equations of motion for an aeroelastic system with concentrated structural nonlinearity in discrete (nodal) coordinates for a given freestream Mach number M ∞ are given as where M v and K v are the structural mass and stiffness matrices, and v = {v 1 , v 2 , . . ., v N } T is the displacement vector of N degrees-of-freedom, and v is the time derivative of v.The external force F v = {F v1 , F v2 , . . ., F vN } T is the aerodynamic force vector in nodal coordinates.The freeplay loads are given by the term F c (δ) which takes the form where δ is the rotational displacement of the root about the freeplay hinge axis and 2δ s is the total rotational freeplay magnitude.
The system described by Eq. ( 19) can be reduced by considering modal coordinates, such that the structural motion is approximated as the linear superposition of a subset of m normal modes Φ v due to generalised displacement ξ.Given the freeplay nonlinearity, the mode shapes in Φ v cannot properly account for localized displacements in the region of the nonlinear hinge.In this work, the fictitious masses (FM) method proposed by Karpel and Newman [3] is used to improve the representation of these local deformations in the set of low frequency modes.A large fictitious mass is added to the DOF of the mass matrix where the discrepancy in localized displacements occurs, then the normal mode shapes are obtained from free vibration analysis and used in the aeroelastic simulation.The baseline FM modes Φ B are derived using ANSYS MAPDL, yielding the generalised system in baseline fictitious mass coordinates where , and Q = Φ T B F v is the generalized aerodynamic force vector.The first eight eigenvalues of the baseline FM modes and baseline normal modes are given in Table 1 where excellent agreement with the results of Carrese et al. [29] can be observed.For comprehensive numerical validation and derivation of the FM model, see Hale et al. [30].
At this point it is important to note that in the general definition of the Taylor series expansion of the unsteady aerodynamic forces, the structural displacements defined by u (Eq. 1) are equivalent to ξ and will be referred to as such from now on.

Computational Fluid Dynamics Model
For the FOM, the generalized aerodynamic force vector Q is obtained using the commercial finite-volume Navier-Stokes solver ANSYS Fluent 2023 R1.The Euler equations for transient flowfields are solved via a coupled pressure-based solver with implicit second-order spatial and first-order temporal discretization of the flowfields with Rhie-Chow: distance-based flux interpolation.The convergence criteria are set to 1 × 10 −4 for the scaled residuals at each timestep.The investigation is conducted on a structured grid of 70 × 10 3 elements, with a minimum orthogonal quality of 0.032.It is important to note that this numerical mesh is validated against experimental campaign [26] via linear stability analysis [28], [27] for the unmodified AGARD wing.Grid deformation is facilitated using a diffusion-based approach.The Modal Projection and force Reconstruction (MPR) method [31] is used to project the structural mode shapes onto the fluid grid.MPR includes a robust interpolation scheme that accounts for disparity in the grid topologies, and conserves forces and moments.
3.4 Unsteady Aerodynamic ROM Library and Interpolation Scheme

Generalized Aerodynamic Forces
The ROMs are generate with knowledge of the modal LCO response.Three separate sets of generalized aerodynamic forces are generated; using generalized displacements that are based on the natural frequencies and modal LCO amplitudes for i) the lowest dynamic pressure q ∞ = 3768 [Pa] ( ξA ), ii) the highest dynamic pressure of q ∞ = 4012 [Pa] with δ s = 0.5 • ( ξB ), and iii) the highest dynamic pressure with δ s = 1 • ( ξC ).Table 2 summarizes the frequencies and maximum amplitudes of the random excitation of the generalized displacements for each mode.Four modes are used in the identification procedure which is the minimum number required for the basis to provide a good representation of the full-order structural model [30].Figure 3(a) presents an example of the band limited random excitation of mode 2 ξ 2,C for a total of n = 400 samples.Figure 3(b) presents the total generalized aerodynamic forces in mode 1 and mode 2, comparing Q F OM , the linear Q ROM and the third-order Q 42 OS−ROM .Although the forces in mode 1 are well captured by the linear model, the third-order ROM demonstrates a more than 2× reduction in error.In mode 2 it is clear that the third-order OS-ROM provides superior performance with a 3 -4× reduction in error compared to the linear ROM, nearly perfectly overlaying the FOM aerodynamic response.

Parametric Subspace
The parameter space is discretized according to nine separate sampled ROM locations Dij T (Fig. 4).The ROMs are specifically targeted at freeplay values 0.5 • ≤ δ s ≤ 1 • , and dynamic pressures 3768[Pa] ≤ q ∞ ≤ 4102[Pa], capturing the nonlinear instability region up to 96% of the flutter boundary.

Hyperparameter Optimization
Hyperparameter optimization is conducted for each sampled ROM, i.e., tuning the ROM to optimize aeroelastic performance at the discrete location of the subspace.For the lowest dynamic pressure linear models are generated using pseudo-inverses where ξA is used to generate the generalized forces.For the highest dynamic pressure, OS-ROMs are generated using ξB for δ s = 0.5 • and, ξC for δ s = 0.75 • and δ s = 1 • .To reduce the computational burden associated with ROM generation, the ROMs for the mid-point dynamic pressures are simply defined as some ratio between the ROM generated for the maximum and minimum dynamic pressures, according to The aeroelastic optimization problem uses the objective argmin||δ F OM (t)− δ ROM (t)|| which is quantified according to where δ is the rotational aeroelastic response at the root hinge node.The hyperparameters of the aeroelastic optimization are summarized in Table 3.

Direct Interpolation Scheme
The discretization and interpolation schemes for the subspace are as follows; i) 2 × 2 discretization of the subspace with bilinear interpolation (ROM 1 4 ), and ii) 3 × 3 discretization of the subspace with biquadratic interpolation(ROM 2  9 ).Given that both linear and third-order ROMs are used, and that the second-order and third-order partial derivatives are identified using the same number of lag terms, then the base vector (Eq.27) is created using The interpolated ROM is then derived for the new freeplay δ s I and dynamic pressure q ∞ I from Eq. 17, for interpolation scheme of order p L by where q ∞ I − q ∞k q ∞l − q ∞k (27)

Aeroelastic Time Integration
The aeroelastic system is solved using the RMIT in-house Fluid-Structure Interaction code PyFSI.Full-order aeroelastic solutions are achieved by marching Eq. 22 forward in time, where the wing transient structural motion is solved using Newmark-β time-integration.Newton-Raphson iterations are used to converge the state-dependent freeplay load within each time-step by minimizing error in the stiffness matrix.The nonlinear fluid loads Q are resolved at every time-step using the CFD model described above.
For nonlinear unsteady aerodynamic ROM solutions, substituting Eq, 6 into Eq.22, the aeroelastic equation-of-motion becomes where M B contains the aeroelastic response in FM coordinates, updated at every time-step by marching the system of equations forward in time using the approach described above.The linear scaling of the generalized aerodynamic forces with dynamic pressure is entirely valid for the interpolated ROM.A numerical time-step of ∆t = 0.001 s is used in all simulations.

Results and Discussion
In this section the results of the aeroelastic ROM are presented and discussed.All simulations are conducted at M ∞ = 0.96 and α 0 = 0 • with varying dynamic pressure and freeplay values.

Aeroelastic Hyperparameter Optimization
Table 4 summarizes the results of the aeroelastic hyperparameter grid search.
For the linear ROMs the optimal number of lag terms is relatively consistent, as to be expected.For the OS-ROMs, the optimal number of lag terms increases as the freeplay magnitude increases, i.e., as the strength of the nonlinearity increases the memory fading nature of the system implies that more lags (or components in terms of Volterra kernels) are required.It is quite remarkable that high-precision ROMs can be identified with less than 15 total coefficients, where the full model is defined by hundreds.Figure 5 presents the phase portraits of the aeroelastic response at the freeplay hinge at the locations of the sampled ROMs.Excellent precision at these discrete locations can be observed, as to be expected.It should be noted that at δ s = 0.5 [ • ], q ∞ = 3768 [Pa], the system is actually marginally stable -captured by the FOM and ROM D11 T .

Out-of-Sample Performance
To demonstrate the ROM out-of-sample performance without interpolation the linear ROM D31 T and third-order OS-ROM D33 T are considered and the dynamic pressure and freeplay magnitude are varied.Figure 6(a) [32] presents the LCO amplitude at the hinge rotational axis as a function of dynamic pressure with a freeplay of δ s = 1 • .It can be seen that, as expected, the third-order OS-ROM does not scale particularly well with dynamic pressure, certainly not as well as what would be expected for an aeroelastic problem that is characterized by linear unsteady aerodynamics.That being said, considering the complexity and nonlinearity in the system, the result is reasonable and in line with the recent findings of Brown et al. [23].The linear ROM that has been tuned for aeroelastic performance at the lowest dynamic pressure of interest performs very poorly when scaled with dynamic pressure.
Figure 6(b) presents the LCO amplitude as a function of freeplay magnitude, dynamic pressure remains constant at the exact values that the linear ROM D31 T and third-order OS-ROM D33 T were calibrated to.The third-order OS-ROM generalizes quite well, aside from an over prediction of the amplitude of the LCO for the lowest freeplay value of δ s = 0.5 • .The linear ROM also performs very well, however, is unable to capture the stable response at δ s = 0.5 • .Despite generally good performance in terms of new freeplay values, the over prediction at the lowest freeplay value means that in order to derive a nonlinear ROM that is highly accurate across the entire parameter space, interpolation is necessary.

Limit Cycle Amplitude
The out-of-sample performance of the two interpolated ROMs is first investigated in terms of the ability to capture the transonic limit cycle amplitudes.
In Figure 7 the contours demonstrate that the biquadratic ROM 2 9 is able to capture the entire parameter space with excellent precision while the bilinear ROM 1  4 is less accurate.The discrepancies are larger in regions of the space that are furthest from the sampled ROM 4  1 locations.The direction of the contours indicates that the LCO amplitude is more sensitive to dynamic pressure than to freeplay amplitude.Figure 8 presents the LCO amplitude predictions for the entire parameter space.As mentioned previously, at the lowest dynamic pressure and freeplay the system is found to be stable and the sampled ROM is tuned to capture this.Although this was not by design, it is an excellent finding that the inclusion of this ROM in the library allows robust identification of the nonlinear flutter point.Figure 10(b) demonstrates that at δ s = 0.5 • the bilinear ROM 1  4 captures the LCO amplitudes at the root rotational axis with very good precision and reasonable precision for δ s = 0.625 • .However, for freeplay values above this, discrepancies can be observed through under prediction of the LCO amplitude as the dynamic pressure increases beyond q ∞ = 3846 [Pa].Figure 8(b) presents the LCO amplitudes modeled using the biquadratic ROM 2  9 .Excellent agreement can be observed between the FOM and ROM solutions -demonstrating that the nonlinear parameter space can indeed be captured with excellent precision using a biquadratic interpolation scheme.Although not presented here, almost identical trends are observed for the tip response using ROM 1  4 and ROM 2  9 .
(a) (b) Figure 9 presents phase portraits at the centroids of biquadratic parameter space, i.e., at the furthest distances from the sampled ROMs.At the lower dynamic pressure region q ∞ = 3846 [Pa], the bilinear ROM 1 4 and biquadratic ROM 2  9 both capture the LCO with excellent precision for δ s = 0.625 • .With a freeplay of δ s = 0.875 • the biquadratic ROM 2  9 is generally in very good agreement with the FOM, however, slightly over predicts the amplitude, making it more conservative.The bilinear ROM 1  4 slightly under predicts the amplitude and over predicts the depth of the inflections at the rotational velocity maxima.
For the higher dynamic pressure of q ∞ = 4012 [Pa] the biquadratic ROM 2   9   captures the limit cycle with excellent precision for freeplay of δ s = 0.625 • , while the bilinear ROM 1  4 under predicts the amplitude and over predicts the inflection points (at the maximum velocity regions of the cycle).With freeplay of δ s = 0.875 • , the bilinear ROM 1  4 severely under predicts the amplitude and also does not capture the sharpness of the inflection points.The biquadratic ROM 2  9 performs well, however, still slightly under predicts the velocity magnitude.

Limit Cycles in Generalized Coordinates
To probe the ROM performance further, the generalized displacements are now investigated comparing the FOM to the biquadratic ROM 2 9 solutions.Figure 10 presents bifurcation diagrams for each mode with freeplay δ s = 0.875 • , i.e., plotting the maxima and minima of the generalized displacements at the five dynamic pressure validation locations.It can be seen that mode 1 is captured with excellent precision.Modes 2 and 3 are also captured well, however, some asymmetry can be observed at the highest dynamic pressure in the ROM which leads to a slight over prediction and there is a general under prediction at the lowest dynamic pressure.The dynamics of mode 4 present more interesting quasi-periodic behavior.It is quite impressive that the biquadratic ROM 2  9 is able to capture these dynamics quite well, including bifurcations that appear to occur between q ∞ = 3846 [Pa] and q ∞ = 4012 [Pa], i.e., the transition from a five period to three-period response.The lowest dynamic pressures present more error for mode 4, which is likely due to the low amplitude at this speed making it difficult to resolve.These discrepancies to not impact the global response in nodal coordinates.
Phase portraits of the generalized displacements are presented in Fig. 11 for the dynamic pressure q ∞ = 4012 [Pa].For the lower freeplay δ s = 0.625 • the biquadratic ROM 2 9 performance is very good for all modes.For the larger freeplay δ s = 0.875 • mode 1 is captured with excellent precision while mode 2 and mode 3 demonstrate small discrepancies -not surprising given that this is at the centroid the most nonlinear quadrant of the parameter space.For mode 4 the five-period response is captured, however, asymmetry and associated discrepancies in amplitude can be observed.From these findings it is clear that the discrepancies that are observed for this case in Fig. 9 are the result of error in modes 2-4.The accuracy of the ROM could be potentially be improved using a multi-input identification framework, or certainly by optimizing the ROM in generalized coordinates, rather than for a single location in nodal coordinates.

Computational Savings
To ensure a stable limit cycle is achieved, 3000 numerical time-steps are required.Any CFD-based aerodynamic solutions are run on an Intel Xeon Gold 6152 CPU using 16 cores, while any ROM solutions on a single core (i.e., insignificant resources).Table 5 presents the computational cost associated with the various solution strategies presented in this paper.In terms of simulation time, the ROM solutions are three orders of magnitude faster.However, the time taken to generate the ROM must also be taken into account.The ROM 1  4 takes approximately 350 CPU-h to generate, while the ROM 2 9 takes approximately 675 CPU-h to generate.

Summary and Conclusion
Non-parametric nonlinear model reduction methods, such as, those based on polynomial functionals, can suffer from reduced ability to generalize.For complex nonlinear problems this means that convenient relations that are observed in the linear regime, such as the linear relationship between generalized force and dynamic pressure, or insensitivity to linear/nonlinear control stiffness properties, can become redundant.In this paper an approach for parameterizing nonlinear non-parametric unsteady aerodynamic ROMs is presented, for application to complex nonlinear aeroelastic systems.The approach considers direct interpolation (using Lagrange polynomials) of a an aeroelastic-tuned polynomial-based unsteady aerodynamic ROM library, consisting of a combination of linear and optimally sparse nonlinear tensors of Taylor partial derivatives, equivalent to Volterra kernels.The primary features of the framework and interpolation scheme are as follows: 1. Through OMP it is possible to efficiently generate the nonlinear ROMs in the library.2. Each ROM in the library can be of different order, sparsity (or lack-of), sparsity pattern, and cardinality.3. The ROM library can be relatively sparse and the interpolation scheme is simple.4. A completely new ROM is generated for each new set of operating conditions and thus the interpolation is only done once, before running the simulation.
The example given considers a 3D aeroelastic stabilator model with freeplay undergoing high amplitude LCO, where the nonlinear unsteady aerodynamic ROM is parameterized in dynamic pressure and freeplay magnitude.A bilinear interpolation of the subspace performs very well for low freeplays and dynamic pressures, while to model the entire subspace with high accuracy a biquadratic scheme is necessary.The online computational savings are more than three orders of magnitude.The time taken to generate the ROM is also of significance.It is estimated that it takes 8-16 times longer to generate the nonlinear ROM library, than generating the linearized aerodynamic impulse responses for this same system.
In conclusion, this work demonstrates that it is possible to generate a highly accurate nonlinear unsteady aerodynamic ROM that is as robust to parameter changes as a linearized ROM within the linear regime, albeit with a larger but reasonable offline computational cost.Generating the same ROM without sparsity promotion would be computationally prohibitive.

Fig. 1 :
Fig. 1: a) Modified AGARD 445.6 wing geometry specifications and b) hinge stiffness as a function of rotation

Fig. 3 :
Fig. 3: Cross validation data; a) band limited random excitation of mode 2, and b) comparison between the FOM and ROM generalized aerodynamic forces for Q4(ξ2)

Fig. 4 :
Fig. 4: Two-dimensional subspace and locations of the sampled ROMs

Fig. 5 :
Fig. 5: Comparison between FOM and ROM limit cycles at the sampled ROM locations

Fig. 9 :
Fig. 9: Phase portraits of the root rotational axis LCO at the centroids of the biquadratic parameter space

Fig. 10 :
Fig. 10: Bifurcation diagrams in generalized coordinates comparing FOM to ROM 2 9 solutions displaying a) maxima of the LCO, and b) minima of the LCO (logarithmic scale)

Table 1 :
Natural frequencies for the root stiff and free cases calculated directly with fictitious masses k δ = 500 [Nm/rad] k δ = 0 [Nm/rad]

Table 2 :
Parameters of the band limited random excitation

Table 3 :
Hyperparameter space for the aeroelastic optimization

Table 5 :
Computational cost for the FOM and ROM solutions