NOMENCLATURE
- $a_i$
temporal coefficient
- $\tilde{a_i}$
continuous spline function
- e
energy
- $N_M$
number of POD modes
- p
pressure
- q
dynamic pressure
- $\bar{q}$
average dynamic pressure
- Re
Reynolds number
- u
flow velocity
- $\bar{u}$
flow average velocity
- $\underline{\underline{W}}$
snapshot matrix
Greek Symbols
- $\gamma$
eigenvalue
- $\rho$
density
- $\Phi$
POD mode
- $\partial$
partial operator
- $\bigtriangledown$
Delta operator
- $\tau$
shear force
1.0 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(Reference Albano and Rodden1,Reference Falkner2) . 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 phenomenon known as limit cycle oscillations (LCOs) is caused by the highly non-linear characteristics of either or both the structure and aerodynamics. From a structural point of view, non-linearities such as the dry friction between pylons(Reference Denegri3), free-play of actuators(Reference Tang, Bartels, Chen and Liu4) and geometric non-linearity may induce LCOs. Regarding aerodynamics, LCOs can be induced by transonic shock movements(Reference Knipfer and Schewe5) and flow separation at high angles of attack.
Recent investigations by Poirel, Metivier, Rudmin, Yuan et al(Reference Poirel, Harris and Benaissa6,Reference Metivier, Dumas and Poirel7,Reference Poirel and Yuan8,Reference Poirel, Metivier and Dumas9,Reference Poirel and Mendes10,Reference Rudmin, Benaissa and Poirel11,Reference Yuan, Poirel and Wang12,Reference Poirel and Mendes13) addressed the LCOs of a clean NACA 0012 aerofoil. These studies employed wind-tunnel experiments and computational fluid dynamics (CFD) analyses of an aerofoil in the transitional Reynolds number regime. The researchers observed LCOs with small to large amplitude depending on different Reynolds numbers, and it was suspected that these LCOs were induced by laminar separation. Translational Reynolds numbers of $ 10^4<\mathrm{Re}<10^6 $ are known to be highly non-linear; their characteristics may result in complex viscous phenomena such as laminar boundary-layer separation, the transition of the laminar shear layer, and flow reattachment(Reference Mueller14,Reference Gad-el-Hak15) .
Due to the above-mentioned highly non-linear aerodynamic characteristics of LCOs, it is recommended that analyses be conducted using non-linear CFD solvers such as Navier–Stokes CFD. However, the fluid–structure interaction (FSI) analysis coupled with the time-domain Navier–Stokes CFD requires a significant amount of time and is computationally extensive, thereby limiting its use in aeroelastic analyses. Even though computing power has grown exponentially over the last few decades, the computational time required for CFD in the aeroelastic field remains excessive. Many researchers have attempted to reduce the computational time by employing reduced-order models (ROMs) and parallel computing(Reference Kheirandish, Beppu and Nakamichi16,Reference Hall, Thomas and Dowell17) . However, parallel computing is not capable of reducing the total accumulative computational resources; rather it reduces the computational time via using simultaneous calculations across multiple central processing units (CPUs). Conversely, ROM reduces the order of the system, thereby diminishing the computational time. Consequently, ROM has been actively investigated since 1990 for use in aeroelastic as well as aerodynamic, pollution, medical and electronic fields. Hall et al endeavoured to apply ROM to aeroelastic systems; the frequency-domain LCO analysis of a two-dimensional aerofoil via Euler CFD and the proper orthogonal decomposition–reduced-order model (POD-ROM) technique is an example of their endeavours(Reference Hall, Thomas and Dowell17). Additionally, Lucia developed the Volterra theory and POD-ROM hybrid methodology for aeroelastic systems(Reference Lucia, Beran and Silva18). Other studies have been performed to analyse incompressible flows using POD-ROM. Park et al and Burkadt et al implemented the POD-ROM technique for incompressible, viscous cavity flow(Reference Park and Lee19,Reference Burkardt, Gunzburger and Lee20) . Galletti et al successfully implemented POD-ROM using temporal coefficient calibration for an incompressible flow past a square cylinder(Reference Galleti, Bruneau, Zannetti and Iollo21). It was concluded that ROM can reduce the computational time when analysing highly non-linear aeroelastic systems, using Navier–Stokes CFD.
This paper presents the development and application of POD-ROM for incompressible Navier–Stokes CFD. POD-ROM is preliminarily validated using a reconstruction of the flow field around a plunging cylinder. Thereafter, the reconstructed flow field is evaluated in terms of its accuracy and the reduction in the computational time. After the validation of the proposed POD-ROM, the flow field of a NACA 0012 aerofoil under LCOs at subsonic speed is analysed. Subsequently, this reconstructed flow field is also evaluated in terms of its accuracy and the reduction in the computational time. To the best of the authors’ knowledge, this paper is one of the first attempts reported in literature to employ POD-ROM based on the Galerkin projection to analyse highly non-linear LCOs in incompressible regimes in the time domain.
2.0 FORMULATION
2.1 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 full-order 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).
where
The total energy is expressed as Equation (3).
2.2 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(Reference Sirovich22,Reference Noack, Afanasiev and Morzynski23) and the projection of the basis vector on the discretised governing equation(Reference Chaturantabut and Sorensen24). 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(Reference Carlberg and Farhat25). 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\times 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)=\int_{\Omega}a . b d\Omega}$ . The velocity at a certain coordinate and time is expressed as Equation (9).
2.3 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(Reference Rowley26). 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(Reference Gallletti, Bottaro, Bruneau and Iollo27). 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 $\tilde{a}_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^{\prime}_k$ and $B^{\prime}_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)
where
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.
3.0 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 two-dimensional cylinder is reconstructed. The temporal coefficients are calibrated, and their impact on the accuracy of the reconstructed flow field is also examined. Subsequently, a two-dimensional 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.
3.1 Flow field reconstruction of a plunging cylinder
3.1.1 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(Reference Zienkiewicz and Cordina28,Reference Nithiarasu, Cordina and Zienkiewicz29) . 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\times10^{-4}$ s, and 126,178 degree of freedoms (DoFs) are used for the cylinder analysis, as presented in Fig. 2.
3.1.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.
a Accumulated energy.
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 is used. The computational time and average discrepancy are presented as functions of the number of snapshots in 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.
a Due to divergence in the optimisation, the temporal coefficients are not calibrated.
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® i7-7700K processor at 4.37GHz. Therefore, the computational time is reduced by 83% in this example.
3.2 Flow field reconstruction for an aerofoil
3.2.1 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(Reference Poirel, Harris and Benaissa6,Reference Metivier, Dumas and Poirel7,Reference Poirel and Yuan8,Reference Poirel, Metivier and Dumas9,Reference Poirel and Mendes10,Reference Rudmin, Benaissa and Poirel11,Reference Yuan, Poirel and Wang12,Reference Poirel and Mendes13) . 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
a Due to divergence in the optimisation, the temporal coefficients are not calibrated.b Total elapsed time without temporal coefficient calibration.
The aerofoil is analysed using incompressible Navier–Stokes CFD with the OpenFOAM PIMPLE dynamic grid solver. The typical time step ranges from 3 $\times10^{-5}$ s to 7 $\times10^{-5}$ s. In the computation, “nOuterCorrectors” is set to 500, “nCorrectors” is set to 3, “nNonOrthogonalCorr” 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.
3.2.2 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 et al(Reference Poirel, Harris and Benaissa6,Reference Metivier, Dumas and Poirel7,Reference Poirel and Yuan8,Reference Poirel, Metivier and Dumas9,Reference Poirel and Mendes10,Reference Rudmin, Benaissa and Poirel11,Reference Yuan, Poirel and Wang12,Reference Poirel and Mendes13) , who reported that the laminar separation and separation bubble were the causes of the LCOs.
3.2.3 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.
a Accumulated energy.
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® 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.
4.0 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.
ACKNOWLEDGEMENTS
This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF), funded by the Ministry of Science, ICT and Future Planning (2017R1A2B4004105).
SUPPLEMENTARY MATERIAL
To view supplementary material for this article, please visit https://doi.org/10.1017/aer.2020.59