Abstract
The lack of stability is a problem encountered when applying the classical POD-Galerkin method to problems of unsteady compressible flows around a moving structure. To solve this problem, a hybrid reduced-order model named POD-ARX is constructed in this paper. The construction of this model involves two steps, including first extracting the fluid modes with the POD technique and then identifying the modal coefficients with the ARX model. The POD modes with the block of all modified primitive variables are extracted from the system response to the training signal. Once the POD modes are obtained, the snapshots are projected on these modes to determine the time history of modal coefficients and the resulting modal coefficients are used to identify the parameters of ARX model. Then, the ARX model is used to predict the modal coefficients of the system response to the validation signal. Sample two-dimensional aerodynamic force calculations are conducted to demonstrate this method. Results show that this method can produce a stable and accurate prediction to the aerodynamic response with significant improvement of computational efficiency for linear and even some nonlinear aerodynamic problems. In addition, this method also shows good wide-band characteristics by using the “3211” multistep signal as the training signal.
1. Introduction
In recent years, with the rapid development of the computer technology, computational fluid dynamics (CFD) has become a more and more common tool in the analysis of flow physics. The ability of performing high-fidelity unsteady flow simulation makes the CFD solver capture complex flow phenomenon accurately, such as separated flows and shock waves. However, due to the huge cost of time and computational resources, nowadays, the CFD technique is still not very suitable for analyzing fluid-structure interaction problems or other unsteady problem that needs repeated calculations. Therefore, it’s very important to develop a surrogate model with high levels of both computational efficiency and accuracy.
In 1990s, a series of unsteady aerodynamic reduced-order models (ROM) were proposed. Compared to the direct CFD numerical simulation, these models can not only improve the computational efficiency significantly but also have the satisfying computational accuracy. Current researches on ROM can be generally divided into two groups. The first branch employs the system identification (SI) methods to build the relationship between the input and output data, including the autoregressive with exogenous input (ARX) model [1, 2], the Volterra series [3–5], and the neural networks [6–9]. The second branch is based on the eigenmodes of the flow field, including the proper orthogonal decomposition (POD), which projects the governing equations onto these eigenmodes to obtain a low-order dynamic model [10–12]. Due to the good performance in the accuracy and efficiency, the POD technique has become an active area of the ROM research and has been widely used in optimal design [13, 14], control design [15, 16], aeroelastic analysis [17–19], etc. However, due to the influence of inner product, simplification of the original governing equations, and the lack of dissipation in numerical schemes, the conventional POD ROMs are usually unstable and additional stabilization strategies have to be adopted to solve these problems [20–24]. In addition, some combinations between POD and system identification methods have also been proposed in recent years [25–27].
Lucia and Beran [26] proposed a hybrid approach by combining POD and Volterra theory. In this method, the POD technique is used to extract the fluid modes and the Volterra theory is employed to identify the modal coefficients of the flow field. This method has been demonstrated on a two-dimensional subsonic inviscid flow over a bump with forcing and been successfully used to predict the limit-cycle oscillation behavior over an elastic panel in supersonic flow [26, 27]. Although this method has shown its efficiency and accuracy in aerodynamic response prediction, it also presents some drawbacks. First, this method extracted the POD modes from the startup response of full-order model to validation signal itself and it showed in [27] that the POD basis derived from the impulse response data was not adequate for modeling the aeroelastic problems accurately, which means that these modes are difficult to represent the dominant feature of the flow generated from other input signals with different frequencies. Second, Lucia and Beran [26] and Lucia et al. [27] conducted the orthogonal decomposition for each variable separately. Although this method can provide the maximum POD power structure for each variable, it also makes the system relatively complicated and heavy. In addition, Placzek [28] points out that since the variables are correlated with each other, a part of physical problems may be ignored with this method. Therefore, it is preferable to form a set of eigenmodes from blocks containing all the fluid variables. But the traditional conservative variables are not appropriate because the Navier-Stokes (N-S)/Euler equations are not quadratic and the Galerkin projection will yield an inadequate implicit form of the modal coefficients.
Fortunately, some researchers, trying to apply the POD method in the compressible flow, provide some excellent ideas. Placzek [28] and Placzek et al. [29] introduced the modified primitive variables into the analysis of the nonlinear compressible flows around a rigid body with motion by POD. With the modified primitive variables, the N-S/Euler equations can be transformed to an explicit format about the modal coefficients by the Galerkin projection. However, the resulting reduced-order model lacks some dissipation and has to adopt the correction method to produce a stable response, which brings some new difficulties, such as the choice of the correction method for different flow conditions and the complex computation of the parameters.
In order to overcome the disadvantages of the above methods, a new POD-ARX model is constructed in this paper. The model is constructed by extracting the modes by POD technique and identifying the modal coefficients by the ARX model. As a difference model, the ARX model is very easy to implement with existing unsteady CFD codes and extremely efficient computationally. It is also well suitable for the multi-input/multioutput (MIMO) identification procedure and has been extensively used in aeroelastic problems [1, 2, 30–33]. Besides, recent researches have proved that the ARX model can be used in the nonlinear problems of an airfoil in pitching motion under large angle of attack in the near-instable flows [34–36]. The modified primitive variables are applied as a block to conduct the proper orthogonal decomposition, and POD modes are extracted from the snapshots obtained by full-order system response to the training signal with a wide range of frequencies. Finally, the performance of the method was validated by a subsonic inviscid flow sample around the moving NACA0012 airfoil.
2. Construction of the Reduced-Order Model
2.1. Overview of Proper Orthogonal Decomposition
The aim of POD is to find a set of optimal orthogonal basis to provide a best approximation to the behavior of the full-order system dynamics. Since the snapshots are usually centered, the problem can be transformed into finding the best basis to approximate the fluctuations of the snapshots around a mean state. Consequently, the problem in the discrete domain can be expressed as follows [26]: where represents the full-order base solution, which can be a steady CFD solution or an average solution of the snapshots used to extract the POD modes; represents a linear transformation between the full-order solution and the reduced-order solution ; and the modal matrix is made up of fluid modes , which is also called the POD modes. The modal coefficients make up the column of the matrix .The optimal basis functions can be yielded by solving the following eigenvalue problem [26]: where is an () matrix, is the number of snapshots and is the number of data points in each snapshot. Snapshots, which are also called samples, represent the solutions of full-order system dynamics at different time. These solutions are generally collected to provide a good variety of flow field behaviors. is the matrix of eigenvectors of , and the vectors make up the column of . is the corresponding nonnegative diagonal matrix arranged in descending order that represents the eigenvalues of the system. Thus, once the snapshot matrix is created, the eigenvalues and eigenvectors of the eigen-equation (2) can be solved. However, the matrix is so huge that it requires massive memory storage and is very time consuming to solve. In practice, the snapshot matrix is usually redundant and can be eliminated by resizing the eigenvectors in and eigenvalues in . Finally, the modal matrix can be solved using (3). The eigenvalue can be interpreted as the weight of contribution to each mode in the POD reconstruction and the “energy” captured by the retained modes relative to the whole set of eigenvalues of the correlation matrix can be defined as follows: where represents the “energy” captured by the retained POD modes. Since the eigenvalue sequence descends rapidly, generally, the first few POD modes can capture more than ninety percent of the total “energy,” which is usually enough to describe the system physics. More details about the method can be found in reference [26].
2.2. Model Reduction
Once the POD basis functions have been identified using the method of snapshots, the Euler equations must be recast to solve the modal coefficients in lieu of the full-system variables. For incompressible flows, this is generally accomplished using the Galerkin method. But for compressible flows, the reduced-order model of the Euler equations is generally difficult to solve due to the resulting implicit formulation. To solve the problem, Placzek et al. [29] introduces modified primitive variables into the unsteady N-S equations, which are then expressed as a polynomial form and adequate for the Galerkin projection.
The Euler equations in a moving frame with the modified primitive variables can be written as follows [29]: where and is the density, is the velocity vector of the fluid, is pressure, is the velocity of the mesh and can be described by , where and represent the translation velocity and the angular velocity of the mesh, respectively. is the position vector relating to the moving frame, is the heat capacity ratio.
Equation (5) can also be rewritten in the quadratic form: where and are defined by
Introduce into (6) and project the two sides of the system equations (6) onto each POD mode ; (6) becomes a quadratic ordinary differential equation: where and the parameters , , , , and are constants related to the basis functions, which can be computed from the analytical expressions and the details can be found in reference [28, 29].
2.3. Overview of ARX
The ARX model is a linear difference model, which describes the response of a system as a sum of scaled previous outputs and scaled values of inputs to the system. For a multi-input/multioutput system, the model can be written as follows: where is the vector of system output at the th step, is the vector of system input at the th step, matrices and are the constant coefficients to be identified, and are called the ARX model orders. According to the equation, the system response at given time can be expressed as an algebraic series of multiplications and additions.
The accuracy of the model depends highly on the inputs used to obtain the training data. There must be as much information about the system’s dynamics as possible to be packed into the training set of data for the identification procedure. By comparing the “3211” multistep input signal with other different input signals, Cowan [33] points out that the 3211 input signal is easy to implement in experiments and can excite the best frequency response. Once the system inputs and outputs are given, the matrices and can be calculated by the least-squared method and the construction of the ARX model is completed.
2.4. Construction of the POD-ARX Reduced-Order Model
It should be noted that the reduced-order model constructed in Section 2.2 is instable because the model lacks artificial dissipation in the process of construction. To reproduce the correct behavior of the original full-order model, Placzek et al. [29] contrasts several correction methods and recommends that the Tikhonov regularization method be a more robust method for the nonautonomous system. Using this method, the accuracy and stability of the response have been significantly improved for both the short- and the long-term behavior of the full-order model. Although this method is an identification method, it still needs the calculation of the parameters , , , , and , which is time consuming and complex. In addition, the choice of the parameters of the regularization method is also difficult.
The objective of this paper is to determine the modal coefficients of (8) using the system identification model with the motion signal being the system input and the modal coefficients being the system output. Furthermore, when , , and are small relative to , , , and , their projections and are also small. That is to say, (8) can be treated as a linear or weak nonlinear system. Therefore, the ARX model is chosen here as the identification model. Then a hybrid ROM can be developed in this paper by combining the POD technique and the ARX model. In this method, the POD technique is used to extract the fluid modes from the snapshots and the ARX model is used to identify the modal coefficients of the flow field. The details of this method are described as follows: (1)Observe and record the snapshots of the system for a predetermined input signal (training signal)(2)Extract the POD basis functions from the snapshots gathered in step (1) by POD technique(3)Project the snapshots onto the basis functions to obtain the time history of the modal coefficients(4)Construct the ARX model using the SI technique, with the input signal in step (1) and the time history of the modal coefficients in step (3) being taken as the input and output data, respectively(5)Predict the modal coefficients according to a new input signal by the ARX model obtained in step (4). Then construct the resulting flow field using the modal coefficients and the corresponding POD modes
3. CFD Code Validation
The full-order solver uses the cell-centered finite volume method to solve the Euler equation. The AUSM+ scheme is used to discretize the computation domain [37], while the implicit LU-SGS scheme is used for temporal integration [38]. For the unsteady calculation, the dual-time stepping method is adopted, with the fourth-order Runge-Kutta scheme for the subiteration.
The solver was validated with the pitching motion of NACA0012 in transonic flow. The airfoil was forced to pitch about its quarter chord and the motion is described as where is the mean angle of attack, is the amplitude of the motion, and the reduced frequency is defined by where is the free-stream velocity of the flow and is the chord of the airfoil. In this case, deg, deg, , and . Figure 1 is the comparison of the aerodynamic coefficients obtained from CFD codes and experiment [39], which shows that the CFD codes are reliable for the calculation of unsteady aerodynamic forces.
(a) Lift coefficients
(b) Moment coefficients
4. Results
4.1. Wide-Band Characteristics of the ROM
Using the POD-ARX hybrid reduced-order model outlined in Section 2.4, the aerodynamic characteristics of a two-dimensional NACA0012 airfoil with pitching motion is investigated here. A C-shape structure mesh is used and the mesh consists of nodes, which means the DOFs of the full-order model are . The full-order Euler solver is employed to provide snapshots for the POD modes, as well as the base solution for comparison with the ROMs results.
The model will be investigated in a subsonic flow with and the motion of the airfoil will be limited to the pitching motion described as (11) at deg.
Cowan [33] compared several different training signals and pointed out that the “3211” multistep input signal has a wide range of frequency. Therefore, the “3211” input signal, including a total of 310 time steps, is adopted here and shown in Figure 2. A set of dimensionless flow field is extracted at each time from startup to provide snapshots for proper orthogonal decomposition. The frequency characteristic of the “3211” multistep signal is shown in Figure 3. It is found that the training signal has a good coverage to the frequency band from nearly 40 Hz to 90 Hz, with the dominant frequency being 65 Hz.
Equation (4) is adopted here to evaluate the “energy” captured by the POD modes. It is found that the first 8 POD modes contain more than 99% “energy” of the flow field. Consequently, the first 8 POD modes are extracted here to conduct the model reduction.
The snapshots are projected onto the basis functions to determine the time history of the modal coefficients. Then, the resulting data is used as the training data to identify the parameters and of the ARX model. A variety of ARX model orders ( and ) are tried until the best fit for the training data is found. The model orders and are ultimately chosen as the best fit for the training data. Then, the “3211” multistep signal is input to the newly constructed ARX model to test if it could accurately predict the modal coefficients of multistep response. Figure 4 shows the comparison of the multistep time history of modal coefficients obtained from the Euler solution and those from the ARX model. Note that the former were obtained by directly projecting the Euler solution onto the basis functions. We can see that the reduced-order model fits the training data extremely well. For a quantitative analysis, a relative error is defined to evaluate the performance of the ROM: where and represent the ROM output and the output of the Euler solver, respectively. According to the above equation, the relative errors of different modes are both less than 0.1%.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
First, the validation signal described by (10) is set to be Hz. The comparison of modal coefficients obtained from Euler and ROM solution to the oscillation signal is shown in Figure 5. Obviously, the results by the model match well with those from the Euler solution, with the relative error of each coefficient being 0.3%, 0.59%, 4.47%, 1.13%, 4.98%, 2.2%, 3.8%, and 2.9%, respectively. Besides, it can be found that the first 2 modes achieve better agreement than the other modes. The reason for this phenomenon may lie in the fact that the first two modes can capture more linear characteristic of the system (the dominant characteristic of this linear system) while the last several modes maintain more nonlinear characteristic.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
For a more intuitive evaluation to the result of the ROM, the flow field is reconstructed according to (1). The pressure distributions on the airfoil predicted by the CFD and ROM at four instants in time T/4, 2T/4, 3T/4, and T during the 5th period are shown in Figure 6. Apparently, the pressure distribution on the airfoil obtained from the ROM shows excellent agreement with the result from CFD.
(a)
(b)
(c)
(d)
In addition, the time histories of lift coefficients and moment coefficients are also compared between the ROM and CFD by the integral of the surface pressure on the airfoil, shown in Figure 7. The relative error of the lift coefficients and the moment coefficients is 1.1% and 4.2%, respectively. It is evident that the model has accurately predicted the aerodynamic response.
(a) Lift coefficients
(b) Moment coefficients
In order to further exhibit the characteristic of this model, two different validation signals with Hz and Hz are input to the ROM to predict the response. The aerodynamic coefficients and the pressure distributions on the airfoil at four instants in time T/4, 2T/4, 3T/4, and T during the 4th period predicted by CFD and ROM are shown in Figures 8–11. The results show that the pressure distributions and the aerodynamic coefficients from the ROM both agree well with those from the Euler solver at Hz (the relative error of the lift coefficient and moment coefficient is 1.1% and 4.9%). By contrast, the relative error of the lift coefficients and moment coefficients is 1.5% and 6.8% at Hz. Obviously, due to the frequency being far away from the dominant frequency of the “3211” signal, the errors are slightly larger at Hz, but the accuracy is still acceptable in unsteady aerodynamic calculation. Therefore, we can see that the ROM developed in this paper performs well in a wide range of frequencies, mainly resulting from the good combination of the model with the “3211” multistep signal.
(a) T/4
(b) 2T/4
(c) 3T/4
(d) T
(a) Lift coefficients
(b) Moment coefficients
(a)
(b)
(c)
(d)
(a) Lift coefficients
(b) Moment coefficients
4.2. Nonlinear Performance of the ROM
To evaluate the performance of the POD-ARX ROM in nonlinear problems, in this section, the ROM will be investigated under a larger mean angle of attack. Figure 12 is the lift coefficient curve of the airfoil when . It is obvious that the system begin to exhibit the nonlinear aerodynamic characteristics when deg.
Figure 13 is the pressure distributions on the airfoil predicted by CFD and ROM at four instants in time T/4, 2T/4, 3T/4, and T during the 4th period when deg and Hz. The corresponding aerodynamic results are shown in Figure 14. The results show that the pressure distributions and the aerodynamic coefficients from the ROM both agree well with those from the Euler solver (the relative error of the lift coefficient and moment coefficient is 2.1% and 5.9%). However, it should be noted that due to the influence of the nonlinearity of the system, the “3211” multistep signal cannot fully excite the physical characteristics of the flow field which is very important for the extraction of the POD modes. Therefore, in this case, the first 500 steps (nearly one and a half period) of the validation signal are chosen as the training signal. In addition, more POD modes () and more ARX model orders ( and ) are chosen to conduct the model reduction. In fact, in the case of large mean angle of attack, it is not only difficult to fully excite all the physical characteristics of the flow field by the training signal but also difficult to obtain a stable steady-state base flow. When we further increase the mean angle of attack to 8 deg, it is found that the steady-state flow becomes unstable and the ROM suffers from a failure.
(a)
(b)
(c)
(d)
(a) Lift coefficients
(b) Moment coefficients
5. Computational Efficiency of the ROM
To evaluate the efficiency of the POD-ARX method in a more intuitive way, in this section, we will make a comparison about the computational time of the full-order solver and ROM. Results from the case in Section 4.1 are summarized in Table 1, where all the computations are run on a personal computer (CPU: 3.2 GHZ, Memory: 8.0 GB) and the calculation time is measured by the CPU time covering five periods of motion.
Apparently, the cost of identifying the ARX model is almost negligible. When the time for training the input signal by the full-order CFD solver (3056 seconds) and extracting the POD modes (5497 seconds) is taken into account, the whole model reduction process will cost 8553 seconds. Therefore, the computational efficiency can be improved by almost fifty percent. Once the ROM is built up, the calculation time is almost negligible when applying the model to the corresponding problems.
6. Conclusions
A new hybrid POD-ARX reduced-order model is developed in this paper by combining the POD technique and the ARX model. The method involves extracting the fluid basis functions with the POD technique and using the ARX model to identify the modal coefficients. The method is tested on a NACA0012 airfoil with pitching motion in subsonic inviscid flow. First, the motion at small angle of attack ( deg) is conducted. The “3211” multistep signal is used as the training signal to obtain the fluid snapshots and finally eight basis functions are extracted by the POD technique. Then, an eight-output ARX model is constructed. Three validation signals with different frequencies are input to the ARX model, respectively. The resulting time histories of the modal coefficients and the aerodynamic coefficients of the flow field both achieve good agreement with those obtained from the full-order CFD solver. When ignoring the time for obtaining the snapshots and extracting the POD modes, the computational time of the ROM can be reduced by nearly five orders of magnitude. Second, to evaluate the nonlinear performance of the method, a pitching motion at a larger angle of attack ( deg) is conducted. Since the nonlinear characteristics of the flow field, the “3211” multistep signal fails to excite the full characteristics of the flow field and hence the first 500 steps of the validation signal is used as the training signal. A total of 16 basis functions are extracted to construct the ROM. The results obtained from the ROM have a good agreement with those from the Euler solver.
The results have proven that the POD-ARX model can provide stable and accurate predictions for the unsteady aerodynamic response efficiently. Furthermore, compared with the POD-Volterra method, the POD-ARX model has good wide-band characteristics due to the easy combination with the “3211” multistep signals in the linear range, which is especially suitable to problems with uncertain frequencies or a wide range of frequencies. In addition, using all the primitive variables as a block also benefits to reduce the numbers of the identification terms. And compared with the correction method proposed in reference [29], the ARX model is very efficient and easy to implement without any additional calculations of the system parameters. However, it should be noted that for nonlinear problems, an appropriate training signal and a stable steady-state base flow is very important to the ROM. The initial investigation about this method is made just for two-dimensional aerodynamic problems in this paper. Future work will concern on the three-dimensional problems and the aeroelastic analysis.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Strategic Priority Research Program of Chinese Academy of Sciences (Grant no. XDA17030100).