{'fraction_non_alphanumeric': 0.032199775938708396, 'fraction_numerical': 0.01814173683639912, 'mean_word_length': 4.963793103448276, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A novel unsteady aerodynamic Reduced-Order Modeling method for transonic aeroelastic optimization

In aircraft design, structural optimization and uncertainty quantification concerning transonic aeroelastic issues are computationally impractical, because the iterative process requires great number of aeroelastic analysis. Emerging Reduced-Order Model (ROM) method is convenient for transonic aeroelastic analysis. However, current ROMs cannot be reused during iteration, thus time cost is still way too large. This study proposed an improved ROM suitable for Arbitrary Mode Shapes (ROM-AMS), which is reusable regardless the variation of design variables. By adopting Principal Component Analysis, ROM-AMS method can significantly reduce the number of basis mode shapes and improve the accuracy of flutter analysis. In an optimization case, the weight of cropped delta wing is reduced by 28.46%, and the efficiency is 900 times higher than traditional ROM method, which demonstrates the feasibility of this method in iterative design process.


Introduction
In aerospace field, design optimization and uncertainty qualification (UQ) play an important role. Design optimization is carried out to obtain a subtle design with minimized weight or other objective features, meanwhile subject to necessary constraints such as safety. UQ is applied to quantitatively estimate the uncertainty, for instance design variable uncertainty, and impact of such uncertainty on the quantities of interest (Beran et al., 2017). In essence, design optimization and UQ are all iterative process, during which values of design variables are constantly updated and large number of tests are conducted.
Of interest here is structural optimization/UQ concerning transonic flutter constrains. Flutter is essentially the dynamic instability of aeroelastic systems. To avoid such disastrous failure in structure, flight speed should be lower than the critical flutter velocity (so called flutter boundary). In flutter analysis, the acquisition of unsteady aerodynamic response occupies the longest computation time, thus a rapid and reliable aerodynamic calculation method is of most importance for iterative process. At present, two branches are applied to calculate aerodynamic force in UQ/optimization: Traditional AeroElascity (TAE) and Computational AeroElasticity (CAE) (Beran et al., 2017). TAE predicts aerodynamic force rapidly using lower-fidelity aerodynamic tools. For example, Doublet Lattice Method (DLM) (Jutte et al., 2014;Werter and De Breuker, 2016;Wan et al., 2003;De Leon et al., 2012), strip theory (Georgiou et al., 2014;Weisshaar, 1981) and ONERA stall model (Stanford and Beran, 2013a) were applied to aeroelastic design optimization in subsonic regime, while piston theory was frequently used in supersonic and hypersonic regime for optimization and UQ based optimization (such as Reliability-Based Design Optimization, RBDO) Beran, 2013b, 2012).
On the other hand, many aircrafts, such as large commercial transport planes and fighters, need to cruise or combat in the transonic regime. Due to the flow nonlinearity, aircraft in transonic regime would undergo a dangerous reduction of flutter boundary, namely so called "transonic dip". However, flow nonlinearity such as shock wave and separated boundary layers cannot be correctly captured by TAE. As a compromise, Prandtl-Glauert correction factor is incorporated in TAE method (Mallik et al., 2013). But this method is still invalid when Mach number is approaching 1. With the development of Computational Fluid Dynamics (CFD), CAE becomes an appealing method in transonic regime and has been applied to design optimization concerning static-state aeroelastic problems (Barcelos and Maute, 2008). To obtain flutter velocity, traditional CAE analysis should be casted at many trial freestream velocities until a critical velocity is found. Such implement is computational impractical for iterative process and mainly used for post verification at present (Schuster et al., 2003;Yurkovich, 2003). To ease the problem, the flutter analysis was treated as eigenvalue problem by using Schur method, because flutter is actually Hopf-Bifurcation (Badcock and Woodgate, 2008). This implement has been utilized for UQ  and UQ based optimization (Marques and Badcock, 2012). In addition, (Chen et al., 2004;Stanford et al., 2015) proposed a novel field-panel method which is conductive to transonic aeroelastic optimization. First, CFD solver is launched to obtain a background steady flow. Linearization about the steady background flow is then carried out for a range of reduced frequencies and interpolated onto a flat-plate wing mesh with a field panel scheme (Stanford et al., 2015). Above process produces Aerodynamic Influence Coefficients (AIC) governing the relationship between pressure and downwash at a set of reduced frequencies. The AIC is independent of the variation of structural parameters, therefore it is reusable during structural optimization.
However, aforementioned Schur methods requires significant modifications to existing codes (Yao and Marques, 2017). Instead, unsteady aerodynamic Reduced-Order Model (Dowell and Hall, 2001;Lucia et al., 2004;Raveh, 2005) is applied in (Yao and Marques, 2017) to model nonlinear aerodynamic response. Existing ROM can be summarized in two branches: the Proper Orthogonal Decomposition (POD)/Dynamic Mode Decomposition (DMD) (Rowley et al., 2009;Schmid, 2010) method and the system identification method (such as Volterra series (Milanese and Marzocca, 2009) and Auto Regressive with eXogenous input model, ARX). At present, system identification based ROMs have been applied to aeroelastic problems extensively, for instance, flutter analysis at high angle of attack (Zhang and Ye, 2007), revelation of complex fluid mechanisms such as frequency lock-in in transonic buffeting flow (Gao et al., 2017b) and low Reynolds number vortex-induced vibrations (Zhang et al., 2015b), transonic Limit Cycle Oscillation (LCO) prediction (Mannarino and Mantegazza, 2014;Zhang et al., 2016Zhang et al., , 2012, control low design for active flutter suppression (Chen et al., 2012), active control of transonic buffet flow (Gao et al., 2017a). In addition, a lot of novel nonlinear ROMs have been developed, for example: Kou proposed two layer ROM which can capture both linear and nonlinear aerodynamic responses , and Multi-kernel neural networks for nonlinear unsteady aerodynamic modeling (Kou and Zhang, 2017b). Maximilian Winter utilized neurofuzzy model to conduct unsteady aerodynamic computation in varying freestream conditions (Winter and Breitsamter, 2016).
Nevertheless, above method, including field-panel method and ROMs, are all fixed at Prescribed Mode Shapes (such ROM is referred as ROM-PMS) and not robust to mode shape variation in structure optimization. Existing ROMs can adapt to the change of modal mass or modal frequency (Song et al., 2011;Wang et al., 2008), but not adapt to the change of mode shapes. Marques demonstrated that neglecting the variation of mode shapes in structure-changeable cases will lead to misleading results . To solve the problem, Zhang et.al developed ROM suitable for arbitrary mode shapes (ROM-AMS) (Zhang et al., 2015a). First, a set of Radial Basis Functions (RBFs) are selected as basis mode shapes to linearly fit the physical mode shapes. Then, different from traditional ROM, ROM-AMS is constructed in basis mode coordinate. After being converted to real modal coordinate, ROM-AMS can replace CFD solver in aeroelastic analysis. For various structures with same aerodynamic shape and flow condition, ROM-AMS is reusable and robust to the variation of both modal frequency and mode shapes. Only one run of CFD process is required in structural optimization/UQ. Under the same framework, (Winter et al., 2017) proposed another ROM-AMS where Chebyshev polynomials are selected as basis mode shapes, which can capture the global features of physical mode shapes.
However, in (Winter et al., 2017;Zhang et al., 2015a), the number of required basis mode shapes is relatively large. Take ARX model as example, with same delay order, the number of parameters to be identified is square with the number of the basis modes. Consequently, excessive number of basis modes will significantly prolongs the time of training and modeling, especially in practical problems where the whole flight envelope in transonic Mach number interval is subject to the flutter boundary. However, less basis modes will enlarge the deviation of flutter prediction. So there is a tradeoff between the model simplicity and fitting accuracy. To make the tradeoff easier, PCA basis is adopted here as basis mode shapes. First, conduct parametric sampling and modal analysis on structures in design variable space. After this step, mode shapes belong to sample structures are selected as the input of PCA. Then, basis mode shapes are chosen from the PCA basis. Further analysis reveals that just small number of basis modes can reach desirable accuracy. Ultimately, ROM-AMS is applied in transonic aeroelastic design optimization for the first time. To verify the effectiveness of PCA basis, the variation of mode shapes is relatively large.

Basis mode shapes
In previous ROM-AMS, basis mode shapes are RBFs scattering in the wing plate. The shape of 2-dimensional RBF is shown in Fig1. After assigned correct coefficients, these RBFs can fit arbitrary mode shape in the wing plate. In (Zhang et al., 2015a), 18 RBFs were exploited, while the fitting accuracy is still less than perfect, which results in 9% deviation in flutter speed. If the number of basis modes continues to increase, the complexity and training time of the ROM will increase significantly. In light of this, the tradeoff seems difficult if RBFs are used as basis modes.

Fig1 Two dimension example and locations of RBF basis
To make the tradeoff easier, PCA is applied here to generate basis mode shapes. In some cases, PCA is also known as POD. First, parametric sampling is carried out in design variable space, and l sample structures are obtained in this step. Then, after modal analysis, n mode shapes  are acquired as the snapshots of PCA, where n l k  , and k is the number of selected modes of each structure. Vector  represents the displacements at each node in the finite element model. As the number and location of nodes may vary in different structures, to form the snapshot matrix, mode shapes of sample structures must be expressed under a standard mesh, namely the normalization of mode shapes, which can be achieved by interpolation. Hence, the snapshot matrix can be formed by normalized mode shapes (2) Hence N-dimensional mode shape  can be reduced to m-dimensional vector  . To determine a proper m, Modal Assurance Criteria is brought in and defined as below MAC=1 indicates a good agreement between simulated mode shape and real mode shapes. The minimum value of m which satisfy (3) is chosen as the number of basis modes, where  is defined by user. In some cases where N is far greater than n , eigenvalue decomposition for T PP will be difficult, whereas After the basis mode shapes being obtained, the j th mode shape of arbitrary structure in design variable space can be expressed as where ij  can be obtained by solving linear equations.
In summary, there are two reasons accounting for why PCA method can simulate mode shapes with fewer basis modes. 1) In contrast to RBF and Chebyshev polynomials, PCA basis are customized for specific problem. 2) Although sometimes the mode shapes vary greatly in iterative process, they still obey some basic patterns. For example, in section 3.1, the second order mode shapes of different structures perform a torsional commonness. Such basic patterns can be easily identified by PCA.

ROM-AMS
Then the corresponding response, namely Generalized Aerodynamic Coefficients (GAC) are calculated as below Hence the excitation-response data set has been obtained. For the discrete-time multi-input multi-output system, the ARX model is selected to construct the unsteady aerodynamic ROM based on basis mode shapes where t denotes the time step, u stands for the input, namely the generalized displacements of basis modes, and f represents the output, namely the GAC of basis modes.

Ai and
Bi are the constant coefficient matrices to be estimated, na and nb are delay orders which is determined by user. Once the training data is obtained, Ai and Bi can be calculated by using the least square method, hence the unsteady aerodynamic ROM on basis mode shapes is constructed. As the scale of Ai and Bi grows quadratically with increasing m, the training time of ARX model should also quadratically increases. In light of this, the number of basis mode shapes should be reduced as much as possible.
However, it should be emphasized that current ROM is based on basis mode shapes, and cannot be used in real modal coordinate. Therefore generalized displacements and GAC should be transformed between real modal coordinate and basis modal coordinate Equation (9) transforms real generalized displacements to basis modal coordinate, and (10) transforms basis GAC back to real modal coordinate. Only by this way can ROM-AMS be coupled in real modal coordinate.

Fig2 Construction of ROM-AMS on basis modes
Fig3 Flowchart of ROM-AMS and flutter analysis Given a training signal in real mode coordinate for arbitrary structure, ROM-AMS can calculate responses without another run of CFD. Thus the excitation-response data set are obtained to train classical ROM in real modal coordinate, which can be written in continuous-time state-space form and the structural motion equation under modal coordinate is written as where structural characteristics M , G and K should be determined through modal analysis. Convert (12) to continuous-time state-space form Coupling (11) and (13), the aeroelastic equation can be written as Hence the stability analysis of aeroelastic system can be transformed to the eigenvalue problem of A. To analysis the flutter characteristics, vg and v plot are obtained by solving the eigenvalue of A under different dynamic pressure. The detailed derivation can be found in (Zhang et al., 2015b).

Differential Evolution algorithm
Differential evolution algorithm (DE) (Storn and Price, 1997) is applied in present work to optimize a wing structure. Comparing with gradient algorithm, DE algorithm can capture a global optimal solution. First, the collection of all structural parameters is considered as an individual, namely x represents the ith parameter, namely the "gene" of individual.
D is the total gene number of an individual. Then, NP individuals are generated by randomly changing the structural parameters in the optimization space to construct the initial population XG0.
Undergoing 3 steps listed below, the next generation of population is obtained: (1) Mutation: In the gth generation XGg ( 0 XG for the first iteration), mutation is conducted should be created again until boundary constraint is satisfied. The collection of ( +1) i g v forms the first trial population XGnext1; (2) Crossover: compared with each other. The individual whose objective function is lower is selected into next generation. When the comparison is accomplished at all index, population of next generation XGg+1 is updated and ready for next iteration; Iteration of the 3 steps should be kept until convergence, as is shown in Fig4.

Numerical Example
The prototype of the model is cropped delta wing (Schairer and Hand, 1999) which is a typical experimental and numerical case for transonic LCO analysis (Attar and Gordnier, 2006;Chen et al., 2014;Gordnier, 2002;GORDNIER and MELVILLE, 2001;Peng and Han, 2011). The aerodynamic features are governed by leading-edge vortex and shock wave (GORDNIER and MELVILLE, 2001). The structural weight is about 0.24kg. As is shown in Fig5, with a beam and lumped mass installed, the original flat-plate wing is altered to parameter-changeable model. Because ROM-AMS is not robust to the variation of aerodynamic shape, the aerodynamic shape is assumed unchangeable during optimization. The framework of optimization is seeking a set of optimal structural parameters to define a lightest structure, meanwhile satisfying the flight envelop which is subject to flutter boundary. To verify the effectiveness of ROM-AMS, 7 design variables are deliberately chosen so that the modal frequency and mode shapes are sensitive to the variation of design variables. Among these variables: x1 and x2 define the position of both end of the beam; x3 defines the spanwise location of lumped mass which is installed at the leading edge, and x4 is the weight of lumped mass; x5 and x6 define the sectional shape of the beam, and x7 is the thickness of wing plate. The detailed definition and the variable space are shown in Fig 5(b) and Table1. To clearly illustrate the variables, Fig 5(b) magnifies the beam and wing plate thickness in an exaggerated way. Table2 compares experimental and computational modal frequencies, which verifies the accuracy of Nastran solver. Tip end position of the beam(ratio of chordtip) 15% 70% x3 Position of the lumped mass(ratio of span) 60% 85% x4 Weight of lumped mass(kg) 0 0.105 x5 Sectional height of the beam(mm) 0.5 2 x6 Sectional width of the beam(ratio of x5) 50% 200% x7 Thickness of plate(ratio of thickness in (Schairer and Hand, 1999)) 50% 100%

Accuracy test of Reduced Order Model
The accuracy of ROM-PMS will be examined first in this section. Original flat plate wing is applied here for the comparison to experimental result. To construct ROM, signals in Fig 6 are utilized as excitation of the first 4 structural modals of flat plate wing. The freestream Mach number is 0.87. (GORDNIER and MELVILLE, 2001) indicates that the viscous influence is small for lower amplitude of structural deflection (just like the training signals), hence Mach number is adequate to define an identical flow as the experiment, and atmospheric parameter at 10km altitude is given as the free stream condition.
To verify the grid resolution, two sets of grids are generated: baseline and refined grid. Both grids are based on the surface mesh shown in Fig 7, which contains 10976 cells in one side. The resolution of surface mesh is slightly higher than its counterpart in (Peng and Han, 2011). Different from baseline grid, refined grid contains 30 layer of prisms in the vicinity of the surface mesh. Baseline grid contains about 0.4 million cells, while refined grid contains 1.2 million. Results of ROM-PMS method as well as experiment are listed in Table 3, which indicate that ROM-PMS method using baseline grid reaches enough accuracy and can be regarded as standard solution. The accuracy of ROM-PMS was also proven in (Zhang et al., 2015a). In addition, the flutter dynamic pressure of linearized aerodynamic method is listed in Table3 marked by p-k method (Schairer and Hand, 1999). Because of the transonic feature, linearized method cannot obtain accurate result. To keep consistent with experiment, psi (1.0 psi=6895.0 Pa) is applied in Table 3 as the unit of dynamic pressure.

2.400
Next, to verify the accuracy of ROM-AMS and its capability for large variation of mode shapes, two typical structures are chosen from design space. Their corresponding parameters are var1 and var2. The mode shapes of two structures are compared in Fig8. Clearly, the modal frequency and mode shapes vary significantly with the change of structural parameter, and higher order modal gives more complex shape. To obtain basis mode shapes, Latin hypercube sampling is executed and 25 samples are created.  Fig10. Fig11 indicates that more basis modes give higher fitting accuracy. Take var1 as an example, when the basis number is lower than 7, the fitting accuracy is not adequate, especially for higher order modals which are difficult to fit. As the number of selected basis increasing, MACs of all modes are approaching to 1. 8 basis modes are used in this article, thus the fitting accuracy is adequate.

Optimization case
In transonic regime, it is flutter boundary that constrains flight envelop. For example, to design a flight envelop for original flat-plate wing, critical dynamic pressure under 5 Mach numbers is calculated at first. The 5 Ma numbers are {0.79, 0.83, 0.87, 0.91, 0.95}. The maximum safe dynamic pressure is assumed 75.6% of critical dynamic pressure according to the typical criteria that the critical flutter speed should be at least 15% higher than the maximum allowable speed. Then, the corresponding safe altitude, namely the flight envelop, of original model under 5 Mach numbers is plotted in Fig 14, which divides airspace into safe zoom and dangerous zoom.
In the following optimization, flight envelop of original flat-plate wing is given as the benchmark for constraint. Only those structures whose flight envelop is lower than benchmark can be identified as feasible structure. For example, the envelop marked with "safe envelop" is lower than benchmark envelop, which means this envelop has a broader safe zoom, and structure having such envelop is safer than original wing. On contrary, the flight envelop across the benchmark boundary should be marked with "dangerous" because the safe zoom is narrower when Ma<0.9. In some other cases, flight envelop of a structure is much lower than the benchmark, which implies structural redundancy and the weight can be further reduced. In summary, the structure should be deliberately designed so that satisfying the constraint without redundant weight.

Ma
where f is structural mass,  is the penalty factor which is usually a large number. C is penalty function: when the constraints are satisfied, C=0，otherwise C=1. In this work, constraints are: 1) flight envelop is lower than benchmark; 2) no single degree of freedom flutter (SDOFF) arises; 3) MAC values of all the modals should not be lower than 0.98. Fig15 is the history of objective function. As C(x) maintains 0 during iteration, objective function is actually the total mass of structure. In contrast to original flat-plate model with mass of 0.24kg, weight of optimal structure is 0.1717kg, which is reduced by 28.46%. Fig16 is the history of design variables (nondimensionalized by optimal values), and the optimal design variables are listed in Table 6.  The horizontal axis represents the sample number, and vertical axis means the actual value of the genes (namely design variables). It can be found that in the early stages of evolution, genetic values scatter in design space. However, with the evolution pushing forward, the population gradually converges, and finally huddle in a narrow interval (except x2 and x6), which indicates the algorithm is nearing termination.
However, not all variables can eventually converge. For example, the populations of x2 and x6 still scatter when G=150 as is shown in Fig17. The reason is: these variables don't have substantial impact on objectives and constrains, which means they can be arbitrary value. These variables cannot evolve since there is no different between 'good gene' and 'bad gene'. It is difficult to know in advance which variables will have an important impact in preliminary design stage, they must be found by observing their variance. In subsequent design process, these untypical variables should be excluded so that the optimization will be clearly-oriented. in initial generation, and G50 , G150 stand for the 50 th and the 150 th (final) generation respectively. It can be found that the flight envelop gradually boosts and attaches to the benchmark boundary, and the space for further optimization is decreasing. In G150, the flight envelop of optimal structure is almost equal to the benchmark envelop, which indicates the optimization is ready to terminate. The flutter results of optimal structure is verified as is shown in Table7 and Table8. 4530 different structures were analyzed during optimization, and the time cost on PC with 8 cores of 3.6GHz and 8G RAM is 57h33min. Table9 compares the time cost for one run of structural analysis through ROM-PMS and ROM-AMS. For a single structure, calculation of aerodynamic responses by CFD takes up the majority of cost time. It is the repetition of flow solving that obstructs ROM-PMS from practical optimization, because the computational cost will be 5.91 years for 4530 structures in above example, in contrast to 57h33min for ROM-AMS.  Lastly, the conclusion in  is verified as below, which revealed the fact that the ignorance of variation of mode shape in optimization will lead to erroneous results. Assuming in above optimization, the mode shapes are fixed and just the modal frequency is altered during optimization. The fixed mode shapes are taken from the original flat plate wing. The flutter analysis of optimal structure with fixed mode shape at Ma=0.87 is conducted as below, whereas modal frequency keeps real value. Real mode shapes and assumed fixed mode shapes are compared in Fig19. The outcome reveals that the flutter speed will be 207.4m/s once the mode shapes are fixed, which is well below the standard solution 297.02m/s (in table8). The reason is apparent: the variation of nodal line in Mode 2 is not considered without ROM-AMS. Fig20 demonstrates the significant deviation of flutter speed in vgω plot. Therefore, to conduct authentic aeroelastic optimization, variation of mode shapes must be considered.
(a) real frequency and mode shape (b) real frequency and fixed mode shape

Conclusion
This article proposed an unsteady aerodynamic reduced-order modeling method which is suitable for arbitrary structure (arbitrary mode shape) in design variable space, and applies this method to structural design optimization subject to transonic flutter constraint. The ROM is based on basis mode shapes by which the real mode shapes of parameter-changeable structure can be fitted, therefore it can be reused even if values of design variables are altered. In this work, Principal Component Analysis of structure samples in design variable space is first conducted, then PCA basis are selected as basis mode shapes. The accuracy tests of two representative structures demonstrate that just small number of PCA basis are required for accurate modeling, in comparison of other types of basis modes. The flutter predictions of those structures agree well with the result of traditional ROM with prescribed mode shapes. In optimization case where 4530 structures are analyzed, structural weight is reduced by 28.46%, while the efficiency is 900 times faster than traditional modeling method. The accuracy test of optimum structure manifests same order precision with traditional modeling method in transonic regime. Lastly, Marques's conclusion is verified, which demonstrates that the variation of mode shapes cannot be neglected in iterative process. At present, ROM-AMS is constructed on the same aerodynamic shape and flow parameters (such as Mach number and static angle of attack). For example, to conduct aeroelastic analysis in a flight envelop, several ROMs should be built at representative Mach numbers. In this case, PCA based approach has greater advantage because less basis number requires less training time.

Acknowledgments
This work was supported by the National Science Fund for Excellent Young Scholars

Appendix:
In this appendix, the accuracy of ROM-AMS is verified by Euler solver. The Critical Flutter Velocity (CFV) intervals calculated by CFD solver under each Mach number are listed as below, where the percentage in Table A