Time-domain non-linear aeroelastic analysis via a projection-based reduced-order model

ABSTRACT The aeroelastic phenomenon of limit-cycle oscillations (LCOs) is analysed using a projection-based reduced-order model (PROM) and Navier–Stokes computational fluid dynamics (CFD) in the time domain. The proposed approach employs incompressible Navier–Stokes CFD to construct the full-order model flow field. A proper orthogonal decomposition (POD) of the snapshot matrix is conducted to extract the POD modes and corresponding temporal coefficients. The POD modes are directly projected to the incompressible Navier–Stokes equation to reconstruct the flow field efficiently. The methodology is applied to a plunging cylinder and an aerofoil undergoing LCOs. This scheme decreases the computational time while preserving the capability to predict the flow field accurately. The ROM is capable of reducing the computational time by at least 70% while maintaining the discrepancy within 0.1%. The causes of LCOs are also investigated. The scheme can be used to analyse non-linear aeroelastic phenomena in the time domain with reduced computational time.

construct the full-order model flow field. A proper orthogonal decomposition (POD) of the snapshot matrix is conducted to extract the POD modes and corresponding temporal coefficients. The POD modes are directly projected to the incompressible Navier-Stokes equation to reconstruct the flow field efficiently. The methodology is applied to a plunging cylinder and an aerofoil undergoing LCOs. This scheme decreases the computational time while preserving the capability to predict the flow field accurately. The ROM is capable of reducing the computational time by at least 70% while maintaining the discrepancy within 0.1%. The causes of LCOs are also investigated. The scheme can be used to analyse non-linear aeroelastic phenomena in the time domain with reduced computational time.

INTRODUCTION
For the past five decades, traditional panel methods such as the doublet lattice method (DLM) and vortex lattice method (VLM) have been the industrial standards for aeroelastic analyses, ever since their introduction by Albano and Falkner (1,2) . These panel methods are inherently linear techniques and thus cannot provide accurate aerodynamic predictions in the presence of highly nonlinear flow characteristics. A well-known non-linear aeroelastic

Navier-Stokes CFD
Incompressible Navier-Stokes equation-based CFD is employed to analyse the aerodynamics. Moreover, as the proposed ROM scheme constructs a reduced-order basis and then multiplies the basis to a full-order equation, which is also known as the Galerkin projection, the fullorder model (FOM) should be pre-acquired before constructing the ROM. For the FOM used in this paper, the Navier-Stokes equation is written by considering the conservation of mass, momentum, and energy, expressed as Equation (1).
The total energy is expressed as Equation (3).

Reduced-order modelling of the incompressible Navier-Stokes equation
The Galerkin projection is employed to construct the ROM, based on Navier-Stokes CFD. For this projection, a basis vector is obtained via proper orthogonal decomposition (POD). Generally, two methods are used for the construction of ROMs: direct projection on the governing equation (22,23) and the projection of the basis vector on the discretised governing equation (24) . The continuous ROM (based on the direct projection method) is mathematically simpler than the discretised ROM. However, while the discretised ROM is computationally efficient, it is usually constructed by methods such as Gauss-Newton with approximation tensors (GNAT) and the discrete empirical interpolation method (DEIM), which in many cases require iterative projections (25) . In this paper, a continuous ROM is used for direct projection on the governing equation. The POD extracts the characteristics of larger-dimension systems (in this paper and many others, FOM results) by applying the least-squares principle. Thereafter, the extracted POD modes can be used as basis vectors for ROM. The pre-acquired results from the FOM are converted into a snapshot matrix with a size of N × S, as shown in Equation (4). N corresponds to the number of degrees of freedom (DoFs), while S corresponds to the number of time steps used.
: : : The fluid velocity and pressure are expressed by their average and perturbed values rather than their complete values. Generally, the flow field results possess a large number of DoFs; therefore, they are represented as a matrix. POD modes are extracted from the matrix via singular value decomposition (SVD), as shown in Equations (5)- (7).
The ROM of the incompressible Navier-Stokes equation is expressed as an inner product of the POD modes and the momentum equation, as shown in Equation (8).
where (a, b) = a bd . The velocity at a certain coordinate and time is expressed as Equation (9).

Calibration of the ROM
The continuous ROM may be numerically unstable in a few cases. It is suspected that this instability is either caused by the use of the higher-order POD modes or due to its inherent instability (26) . Additionally, the difference between the optimal energy and the dynamic modes due to POD may also result in instability. Therefore, the calibration of temporal coefficients for ordinary differential equations proposed by Galletti is considered in this paper (27) . The temporal coefficient of the POD modes collected from the FOM is expressed as Equation (10).
The error in the ROM is expressed as Equation (11), whereã m (t) is the continuous spline function that interpolates the set of discretised temporal coefficients. This continuous spline function can be easily extracted from the snapshot matrix. Furthermore, a minimisation of the difference is also attempted.
The optimisation problem for the ROM is defined by Equations (12) and (13) where A k and B k are the calibrated ordinary differential equation (ODE) coefficients, and b k is the Lagrange multiplier. Moreover, the optimisation problem of the cost function J can be written as Equation (14) a Even though the above-mentioned traditional calibration methods for the temporal coefficients may be effective, they may not be accurately calibrated if the number of POD modes increases or if the FOM is initially irregular. It is suspected that their inaccuracy arises due to the use of an excessive number of variables during the optimisation. To overcome such drawbacks, a non-intrusive ROM (NIROM) can be considered. Popular NIROM techniques include the use of radial basis function interpolation and artificial neural networks. In many cases, NIROM techniques are based on POD modes. As POD-NIROM is a data-driven approach, it enables the calibration of successive temporal coefficients. Therefore, it will be examined in a future work, based on the proposed POD-ROM.

NUMERICAL RESULTS
The present ROM approach is verified using two examples, and the relevant procedures are presented in the following subsections. First, a flow field surrounding a plunging twodimensional cylinder is reconstructed. The temporal coefficients are calibrated, and their impact on the accuracy of the reconstructed flow field is also examined. Subsequently, a twodimensional aerofoil undergoing LCOs is considered. The accuracy and the reduction in the computational time required for the reconstruction of the complex non-linear flow field are examined. In this section, the POD-ROM results for pure aerodynamics are presented.

Description of the analysis
The reconstruction of the flow field surrounding a simple cylinder under the prescribed plunging motion is attempted by using the proposed ROM for its validation. A two-dimensional characteristic-based split transient incompressible Navier-Stokes solver is employed for the  aerodynamic analysis (28,29) . The plunging cylinder travels at 1m/s and the Reynolds number is 200, as illustrated in Fig. 1.
The FOM is constructed by the Navier-Stokes CFD solver, and the two-dimensional flow field is discretised using three-node triangular elements. A no-slip boundary condition is applied to the wall of the cylinder, and an arbitrary Lagrangian-Eulerian (ALE) scheme is used to deal with the plunging motion. Using the ALE scheme, the grid plunges in accordance with the prescribed motion in a rigid manner. The time step is fixed at 1 × 10 −4 s, and 126,178 degree of freedoms (DoFs) are used for the cylinder analysis, as presented in Fig. 2.

Reconstruction of the flow field
From the FOM, which is constructed for 0-100s, the fully converged latter half (50-100s) is used to construct the snapshot matrix. The snapshot matrix is constructed by collecting the FOM results at 0.01-s interval; a total of 5,000 snapshots are used for the ROM. A total of eight POD modes are constructed via SVD. These eight modes contain 99.8% of the total  energy of the flow field. These eight POD modes and the average flow field are depicted in Fig. 3, and the energy ratio in terms of the POD modes used is presented in Table 1. Thereafter, the temporal coefficients are calibrated to improve the accuracy. The calibrated and non-calibrated temporal coefficients are shown in Fig. 4, where the calibrated temporal coefficients exhibit better similarity to the FOM. Even though a few temporal coefficients yield reasonable accuracy without calibration, several modes still require this calibration to achieve sufficient accuracy. Figure 5 depicts the phase plot among the temporal coefficients. It is observed that the temporal coefficients are successfully calibrated, and the resulting flow field is significantly similar to that predicted by the FOM. The POD modes and corresponding temporal coefficients are acquired, and the flow field is reconstructed. The flow field reconstructed using the proposed ROM and its original counterpart are presented in Fig. 6. The average discrepancy between the reconstructed fluid velocity field and that of the original is found to be 0.1%. The discrepancy in terms of the POD modes is presented in Table 2. As a greater number of POD modes are used, the flow field discrepancy decreases. However, if 10 or more modes are used, the average discrepancy increases. This phenomenon occurs as the use of too many POD modes leads to divergence of the temporal coefficients. The computational time required for the POD-ROM is presented as a function of the number of POD modes used in Table 3. In this case, the number of POD modes used does not affect the computational time greatly, except for the temporal coefficient calibration. Since the current temporal coefficient calibration scheme relies on optimisation, the use of a larger number of POD modes leads to a rapid increase in the computational time. On the other hand, the rest of the computations are less affected by the number of POD modes used, because the snapshot matrix is small. However, it is found that the computational time required for most of the processes increases slightly when a greater number of POD modes   Table 4. The number of snapshots used affects the computational time greatly via the POD mode and coefficient construction, and the computational load of the temporal coefficient calibration increases rapidly. However, the average discrepancy is not affected by the number of snapshots used as long as the time length of the snapshot results collected is sufficient. The number of snapshots (the time interval and total time length) should be sufficient to capture the energy; That is, it should be capable of capturing the oscillations of the higher modes (determined by the time interval) and the lowest one (determined by the total time length). Also, the number of snapshots collected should be sufficiently larger than the number of POD modes used. Otherwise, the energy ratio obtained for the POD modes may not be accurate, resulting in a large discrepancy from the FOM. In this paper, the smallest number of snapshots was used, with 1,000 snapshots being found to be sufficient.
The results indicate that the proposed ROM requires a computational time of 4h for a 50-s flow field sample, including the construction of the POD modes and the temporal coefficients,  and their calibration. However, note that the ROM can be constructed only after the FOM is acquired. On the contrary, the Navier-Stokes CFD FOM requires 24h to analyse the 50s of flow field, using the same single-core Intel R i7-7700K processor at 4.37GHz. Therefore, the computational time is reduced by 83% in this example.

Description of the analysis
This section presents the analysis of a NACA 0012 aerofoil undergoing LCOs. The geometry and characteristics of the aerofoil are modified from the wind-tunnel test and the CFD analyses conducted by Poirel, Metivier, Rudmin and Yuan et al (6,7,8,9,10,11,12,13) . The aerofoil has a chord length of 0.156m and is constrained in all directions, except for the pitch spring including damping. The aerofoil travels at a constant speed of 10m/s, and the Reynolds number is approximately 100,000. The current setup is illustrated in Fig. 7 The aerofoil is analysed using incompressible Navier-Stokes CFD with the OpenFOAM PIMPLE dynamic grid solver. The typical time step ranges from 3×10 −5 s to 7×10 −5 s. In the computation, "nOuterCorrectors" is set to 500, "nCorrectors" is set to 3, "nNonOrthogonal-Corr" is set to 2 and the maximum Courant number is set to 1. The flow field is discretised using two-dimensional three-node triangular elements. The grid consists of 27,816 DoFs and  is illustrated in Fig. 8. The structural movement is modelled by 6-DoF rigid-body motion, and for convenience when constructing the ROM, an overset scheme is used to interpolate the velocity and pressure data on the overset grid. The dynamic grid capability provided by OpenFOAM is used to build the FOM. OpenFOAM creates dynamic grids in which the far-field boundaries are attached at the initial position and grids as the aerofoil rotates. The dynamic grids rotate the dicretisation proportionally at the specified region along with the aerofoil. However, for convenience, when constructing the ROM, an overset scheme is integrated in the offline stage. The velocity and pressure components of FOM (OpenFOAM dynamic grid results) are interpolated on the stationary background grids. For the ROM construction, only the background grids (overset, stationary and non-deforming grids) are used. By adopting the overset scheme, the challenging task of implementing the dynamic grid capability in the ROM is avoided and grid deformation is eliminated. The background grid consists of 291,428-DoF three-node triangular elements. A schematic of the background discretisation is shown in Fig. 9.

Limit-cycle oscillations of the aerofoil
The pitch response of the LCOs is shown in Figs 10 and 11. As indicated by the pitch response, the aerofoil exhibits a periodic and limited amplitude response. The cause of the LCOs is evident, as the flow separation at the trailing edge creates a lower-pressure zone in the relevant area. This region of flow separation is shown in Fig. 12. The lower-pressure zone created by the flow separation is highly non-linear and facilitates pitch movement. These observations corresponds to those of the research conducted by Poirel, Rudmin and Yuan   (6,7,8,9,10,11,12,13) , who reported that the laminar separation and separation bubble were the causes of the LCOs.

Flow field reconstruction of the aerofoil by the ROM
The flow field of an aerofoil is reconstructed using the proposed POD-ROM, and from the FOM constructed for 0-100s, the fully converged latter half (50s) of which is used to construct the snapshot matrix. The snapshot matrix is constructed by collecting the FOM in 0.01-s intervals, and a total of 5,000 snapshots are used. Subsequently, the FOM result is interpolated on the ROM overset grid, where the number of DoFs increases significantly. A total of 200   POD modes representing 99.95% of the energy are used to reconstruct the flow field. The first eight POD modes and the average flow field are depicted in Fig. 13. The energy ratio is presented as a function of the number of POD modes in Table 5. Thereafter, the flow field is reconstructed using the POD modes and corresponding temporal coefficients. The original and reconstructed pitch angle responses are presented in Figs 14 and 15, respectively. The reconstructed and original flow fields around the aerofoil are shown in Fig. 16. From Figs 14-16, it is evident that the flow field is reconstructed accurately because the pitch, pitch phase and fluid velocity contours exhibit only minor discrepancies. The average RMS discrepancy of the fluid velocities at an arbitrary time step (t = 55s) is found to be    0.061%, which is satisfactory. Table 6 presents the average discrepancy in terms of the number of POD modes used. The discrepancies decrease as the number of POD modes is increased. However, when a greater number of POD modes is used, this decreasing trend moderates.
Also, Fig. 17 shows the average discrepancy in terms of time, revealing an arbitrarily distributed trend but where the average discrepancy increases slightly as the time increases. It is suspected that this is because the POD coefficients are not calibrated and those for a few higher POD modes may diverge in time.
Regarding the computational power required, the 50s of FOM construction requires 39h using a single-core Intel R i7-7700K processor at 4.37GHz. Contrarily, ROM analysis requires 11h when using the same CPU and core. The offline interpolation process requires 7h, the snapshot matrix construction takes 1.7h, calculating the POD modes and temporal coefficients requires 4.7h and the flow field reconstruction requires 4.6h. Thus, the computational time is reduced by 72%. The computational time required when using different numbers of POD modes is presented in Table 7. Unlike the cylinder example, the LCO simulation contains a larger number of degrees of freedom and the computational time increases as the number of POD modes is increased. However, the snapshot matrix construction and ROM reconstruction are only slightly affected by the number of POD modes, whereas the POD mode and coefficient construction is greatly affected.
The proposed POD-ROM is found to be accurate for the reconstruction of the flow field generated by CFD. Throughout the relevant procedures, the computational time is significantly reduced, thereby making use of CFD for analysing aeroelasticity is computationally feasible.

CONCLUSION
A POD-ROM is developed for accurate and computationally feasible aeroelastic analyses. A precise analysis of the LCO caused by flow separation is performed via incompressible Navier-Stokes CFD. The proposed POD-ROM is also examined by reconstructing the flow field of a cylinder and an aerofoil undergoing LCO. By applying the relevant procedures, the flow fields are reconstructed and found to be similar to those obtained by FOM. Additionally, the computing resources required are significantly reduced. In the case of the cylinder, FOM requires a computational time of 24h for the analysis of a 50-s flow field, whereas ROM requires 4h, which is a reduction of 83%. For reconstructing the flow field of an aerofoil undergoing LCO, the FOM requires 39h, whereas the ROM requires 11h, which is a reduction of 72%. The authors are planning to extend the present ROM to fluid-structure interactions (FSIs) in the future. The final version of the ROM for FSI will include POD coefficient calibration/interpolation/extrapolation using an artificial neural network (ANN). Such an extended ROM for FSI will be capable of efficiently analysing transient phenomena such as LCOs or flutter of aircraft. Once the ROM framework including FSI is constructed, the ROM will be used for efficient and accurate aeroelastic analyses using CFD. Aircraft manoeuvres and structural deflection will be obtained. The main purposes of constructing the ROM for aeroelastic analysis may be parametric investigation. The resulting ROM is expected to produce reliable results across parameters such as time or airspeed with significantly decreased computational time.