Abstract
In order to calculate the limit cycle oscillations and bifurcations of nonlinear aeroelastic system, the problem of finding periodic solutions with maximum vibration amplitude is transformed into a nonlinear optimization problem. An algebraic system of equations obtained by the harmonic balance method and the stability condition derived from the Floquet theory are used to construct the general nonlinear equality and inequality constraints. The resulting constrained maximization problem is then solved by using the MultiStart algorithm. Finally, the proposed approach is validated, and the effects of structural parameter uncertainty on the limit cycle oscillations and bifurcations of an airfoil with multiple nonlinearities are studied. Numerical examples show that the coexistence of multiple nonlinearities may lead to low amplitude limit cycle oscillation.
1. Introduction
A great attention has been paid to limit cycle oscillations (LCOs) of nonlinear aeroelastic systems, and there are many numerical methods to predict the amplitude and frequency of LCOs. For example, the classical flutter problems were studied by using the harmonic balance method and the incremental harmonic balance method (IHBM) [1]. A frequency relation that depends on airfoil parameters was derived in [2]. In [3], the airfoil flutter with multiple strong nonlinearities was investigated via the incremental harmonic balance method. Besides the researches of finding LCOs, the bifurcation analysis of LCOs is another important part of nonlinear aeroelastic investigations [4–6]. The bifurcation of an airfoil with both structural and aerodynamic nonlinearities in a supersonic flow field was investigated in [7] by using the center manifold theory and complex normal form method. Chen and Liu [8] studied the supercritical as well as subcritical Hopf bifurcation of an airfoil with cubic nonlinearity.
In structural dynamics, some physical parameters may vary in an uncertain way, and uncertainty quantification is increasingly becoming an important part of computational science. Many approaches have been developed to solve uncertainty quantification problems. The most popular methods in the literature for quantification of aeroelastic uncertainty are sampling-based schemes. For instance, the Monte Carlo method and a response surface method were adopted in [9] to map the relationship between the input variables and the output parameters. Nevertheless, traditional Monte Carlo-based simulation is typically too expensive. To overcome this limitation, some other methods [10] have been developed. For example, the effects of the linear and cubic stiffness coefficients on LCOs and bifurcations of an airfoil with multiple nonlinearities were analyzed in [11] through the equivalent linearization method. Recently, in [12] the polynomial chaos expansion (PCE) was selected to quantify the effects of parameter uncertainty on response variability in aeroelastic system, and the limit of the PCE was also revealed. More recently, the LCOs and the Hopf bifurcations of an airfoil with uncertainties from the system parameter and the initial condition were investigated in [13] via the stochastic collocation method. Most recently, since the maximum vibration amplitude of LCOs is useful to analyze the characters of nonlinear systems, the bounds and peak responses of the LCOs were determined in [14] by applying the continuation method along with the incremental harmonic balance method. Moreover, in [15] this method was extended to detect bifurcations of an airfoil flutter system with hysteresis nonlinearity.
The maximum vibration amplitude of nonlinear system with parameter uncertainty is of particular interest to researchers. The application of the harmonic balance method or its variants to obtain the LCOs results in a system of nonlinear algebraic equations which may be resolved by a nonlinear solver such as Newton-Raphson algorithm. However, the nonlinear solver must be used recursively in order to perform parametric studies. Furthermore, the harmonic balance method in [14] that sweeps the uncertainty parameter space with a fixed step would become time-consuming subject to increasing the number of uncertainty parameters. In addition, the combined effects of parameter variations on the nonlinear dynamics cannot be investigated via the continuation method when all the parameters studied vary simultaneously. Obviously, the computation of the maximum vibration amplitude itself is just an optimization problem taking the maximum vibration amplitude as the object function. Moreover, global optimization techniques have evolved a lot in recent years especially for the MultiStart algorithm. Therefore, the objective of the present work is to develop a methodology to determine the LCOs with maximum vibration amplitude of nonlinear aeroelastic system.
Remainder of this paper is organized as follows: the general formulation of the proposed method for determining the LCOs and bifurcations is presented in Section 2. Validations of the method are then conducted in Section 3, and the LCOs and bifurcations of an airfoil with multiple nonlinearities and uncertainties are also analyzed. Finally, concluding remarks are presented and discussed in Section 4.
2. The Proposed Method
A detail method (named constrained optimization harmonic balance method) which combines the harmonic balance method and the Floquet theory [16–18] along with a GlobalSearch algorithm [19] (called the OQNLP MultiStart algorithm) is proposed in this section. The proposed method is formulated as a constrained, nonlinear optimization problem. As described previously, this paper is devoted to study the maximum vibration amplitude in nonlinear structure. Therefore, the maximum vibration amplitude of nonlinear structure can be set as the optimization objective of the nonlinear optimization problem.
Within the framework of the nonlinear constraints optimization theory, the nonlinear algebraic equations derived from the harmonic balance method are treated as the generally nonlinear equality constraints, and the nonlinear inequality constraint to determine the stability of limit cycle is obtained by using the Floquet theory. With these nonlinear equality and inequality constraints, the MultiStart algorithm [19] is selected to find the maximum vibration amplitude of nonlinear systems. It should be mentioned that the proposed method is the first to employ a global search optimization framework with nonlinear equality and inequality constraints.
2.1. The Nonlinear Equality Constraints Derived from the HBM
To obtain the LCOs of nonlinear system, the harmonic balance method is adopted. To present the method, the following equation of motion is considered: where , , and are, respectively, the generalized mass, damping, and stiffness matrices, and are, respectively, the displacement and force vectors, and is the vibration frequency. It should be mentioned that contains both internal and external nonlinear forces.
In the HBM, can be written as a Fourier series up to the th term as follows:
Substituting (2) into (1) and applying the Galerkin procedure [20, 21] yield the following nonlinear function: where corresponds to the harmonic coefficients of the nonlinear forcing term; and are, respectively, defined by
The difficulty with solving (3) is in finding a relationship between and since the harmonic coefficients of the nonlinear forcing term are implicit functions of the harmonic coefficients of the displacement. To overcome this difficulty, the alternating frequency time technique is employed.
The key idea of the harmonic balance method is to find the unknown harmonic coefficients in (3). As the conventional HBM is used, (3) is a set of nonlinear equations being directly solved by a Newton-Raphson-type procedure. However, when the response of the nonlinear structure is desired over a range of frequencies, the solution of (3) must be repeated for each frequency. Furthermore, the application of the Newton-Raphson-type method requires that the number of unknown variables is equal to the number of nonlinear equations; that is, (3) is a well-defined nonlinear system. If the nonlinear system of algebraic equations is underdefined, in which the number of the unknown variables is greater than the number of nonlinear equations in (3), the root finding algorithm cannot be used to determine the unknown harmonic coefficients. Fortunately, periodic solution can be obtained in another way once the nonlinear function of (3) is satisfied. Therefore, unlike the traditional implementations of the harmonic balance method, the nonlinear function of (3) is used to construct the nonlinear equality constraints of optimization problem.
2.2. The Nonlinear Inequality Constraint with Stability Analysis
After the limit cycle is constructed, the stability analysis of such limit cycle should be conducted. The stability analysis can be performed by studying the complex eigenvalues (Floquet multipliers) of the monodromy matrix [18] associated with the limit cycle of (1). In order to compute the monodromy matrix, (1) can be rewritten in the state space form as follows: where denotes the state space vector.
Dividing both sides of (5) by which is the initial condition of , some rearrangements yield with , where represents the identity matrix.
The monodromy matrix is obtained by time integration over one period of the differential system shown in (6). In the following, the fourth-order Runge-Kutta method is selected to solve the differential equation of (6). Using as the initial integration point, the state space vector at time is computed as follows: where is the integration step, denotes the state space vector evaluated at time , and , , , and are, respectively, defined by
Differentiating (7) with respect to leads to where
Finally, the eigenvalues (Floquet exponents) of the matrix evaluated at time are obtained [18]: , . The Floquet multipliers corresponding to the Floquet exponents can be used to determine the stability of LCOs. With the application of Floquet theory, the necessary and sufficient conditions for solution stability require (or ) for all . Therefore, the following stability criteria can be inferred: where the vectors and are defined as , and , respectively.
Using the Floquet theory, the stability of limit cycle is examined on a basis of eigenvalues resulting from the initial value problem (6) with the identity matrix as initial condition. Finally, the stability study results in the nonlinear inequality constraint.
2.3. Optimization with GlobalSearch
As our objective of optimization is to find the optimal solution aiming to maximize the vibration amplitude of the nonlinear system of (1), thus the nonlinear equality constraints of (3) and the nonlinear inequality constraint of (11) have to be combined and solved simultaneously. Therefore, the following nonlinear optimization problem can be formulated: where is a set of design parameters and/or uncertainty parameters and is the unknown frequency to be determined.
The solution of this nonlinear optimization problem with respect to a vector of unknown gives a resonance frequency and a vector of harmonic coefficients at a set of parameter values . Their accurate and effective calculation is a very important problem, which is discussed in the following section.
The choice of optimization algorithm is very important, because the final results are usually dependent on the specific algorithm in terms of accuracy and local minima sensitivity. Evolutionary algorithms are less sensitive to local minima; however, they are time-consuming, and constraints have to be included as a penalty term to the objective function. On the other hand, gradient-based algorithms can lack in global optimality but allow multiple constraints and are more robust, especially for problems in which a large number of constraints are prescribed. In this investigation, the advanced OptQuest nonlinear programs (OQNLP) MultiStart gradient-based algorithm from [19] is implemented, and the gradients are approximated by finite differences. For the purpose of completeness sake, short descriptions of optimization algorithms used in this study are briefly outlined next.
2.3.1. The OQNLP MultiStart Algorithm for Global Optimization [19]
The authors in [19] described the OQNLP MultiStart algorithm for global optimization which enable them to find feasible solutions to a system of nonlinear constraints more efficiently. The OQNLP MultiStart algorithm uses the OptQuest Scatter search algorithm to generate candidate starting points for a local NLP solver. This paper adopts the OQNLP MultiStart algorithm to seek a feasible solution to a system of nonlinear constraints.
The optimization problem of (12) can be rewritten as the general nonlinear optimization formulation as follows: where is the vector of length optimization variables, is the objective function, which returns a scalar value, and the vector function returns a vector of length containing the values of the equality and inequality constraints evaluated at .
For constrained optimization problem, most methods were based on penalty function methods that transform into an unconstrained function , consisting of a sum of the objective and the constraints weighted by penalties. By using penalty function methods, the objective is inclined to guide the search toward the feasible solution. The exact penalty function [22–24], which is the most well-known exact penalty function, is used as a merit function for evaluating candidate starting points. For the problem of (13), this function is where are positive penalty weights and the function equals the absolute violation of the th constraint of problem (13) at the point .
A flowchart of the OQNLP MultiStart algorithm is shown in Algorithm 1 in which stands for the starting point generator and represents the candidate starting point produced. The OQNLP MultiStart algorithm uses the scatter search algorithm [25] to generate a set of potential start points. For a description of the scatter search algorithm, see [25]. Starting from the point , the local NLP solver produces the final point . The function Update Locals is then used to process and store solver output and the updated penalty weights, and is also calculated.
| ||||||||||||||||||||||||||||||||||||||||||||||
There are two stages of the algorithm. The algorithm performs , and iterations for the Stages 1 and 2, respectively. At each iteration, the starting point generator produces the candidate starting point , and this point is also used to calculate the exact penalty value . After finishing iterations of Stage 1, the local solver is called at the best point that has the smallest value of in Stage 1. In Stage 2, the MultiStart algorithm runs the local solver starting at the points that pass the distance and merit filters. For the distance filter, the distance factor is used to determine whether a point is in an existing basin of attraction. For the merit filter, the threshold factor and wait cycle are used to update threshold. The basin radius of the distance filter and the threshold of the merit filter are updated during the optimization process. A basin radius decreases after wait cycle consecutive start points are within the basin. These two filters help the MultiStart algorithm to run with a few points and with high success rate. If the local solver runs starting from the point that passes the distance and merit filters, it can yield a positive exit flag, which indicates convergence. A detailed analysis of the algorithm can be found in [19].
2.3.2. Sequential Quadratic Programming
The Sequential Quadratic Programming (SQP) methods [22–24] are known to be powerful when solving problems with significant nonlinearities, so the SQP method will be used in the MultiStart procedure. In the following, the principle of the SQP method is briefly introduced.
The main idea of the SQP method is to formulate the following Quadratic Programming (QP) subproblem to be solved: where is defined as the search direction and denotes the Hessian of the Lagrangian function which is updated using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) updated formula.
The method used to solve the QP subproblem is the active set strategy. The solution of (15) is used to form a search direction for a line search procedure. In other words, the solution is used to form the next iteration:
The step length parameter is determined by an appropriate line search procedure so that a sufficient decrease in a merit function is obtained. The algorithm above only outlines the principles of the SQP method. Details of the SQP procedure can be found in the literatures [22–24] and references therein.
This section proposes a method that returns the maximum vibration amplitude of a structure. The result is obtained by a concatenation of 3 methods: the harmonic balance method to turn the dynamical problem into an algebraic one, the Floquet theory to determine the stability of the solution, and finally a MultiStart algorithm to find the one point that satisfies the HBM equalities plus the stability constraint inequality and which maximizes the vibration amplitude.
It should be emphasized that the use of the optimization algorithm described in Section 2.3 is particularly important. According to our optimization experiments, other optimization methods such as Genetic Algorithm are very difficult to find an optimal solution that is better than the optimal solution obtained by the MultiStart algorithm.
3. Application to the Nonlinear Aeroelastic System
In this section, efficiency of the proposed methodology based on the combination of the Harmonic Balance Method with the Floquet theory will be firstly demonstrated through numerical simulations of an aeroelastic system with the flow speed fixed. Then, the LCOs and bifurcations of the aeroelastic system against uncertain parameters from nonlinear stiffness properties and the flow speed will be analyzed.
3.1. Flutter Model
A classical airfoil with cubic nonlinear which has been used extensively in the literatures [1, 7, 8, 11, 14, 15, 26] is used. The equation of motion for the airfoil is given by where , ; , represent, respectively, the plunge and pitch degree of freedom (DOF); are the nonlinear plunge and pitch cubic stiffness coefficient, respectively, is the flow speed. The flutter model described in (17) adopts steady aerodynamics theory.
To compare with the previous results obtained in [14], the system parameters under consideration are chosen as
3.2. Calculating the LCOs at a Given Flow Speed
In order to validate the proposed method, only the plunge nonlinearity is considered. Parameter values used in this simulation are , , and . These values were selected to obtain a configuration that yields a subcritical instability [12]. For subcritical instability, two LCOs can coexist for a given flow speed, one of them being unstable. In the following, numerical simulations are presented to capture the multiple LCOs at the given flow speed.
While looking for multiple solutions with the fixed flow speed, there are two categories of optimization variables, namely, the harmonic coefficients and the vibration frequency. It should be noted that the objective function in (12) is the maximization of the vibration amplitude for searching the two LCOs, while the unstable LCO is obtained by reversing the sign of the inequality constraint of (12).
Based on the OQNLP MultiStart algorithm described in Section 2.3.1, optimizations are then performed to find these LCOs and the parameters of the optimization algorithm are listed as follows: the number of trial points was chosen to be 1000, and the usual value of 200 for the number of Stage 1 points has been taken. The maximum number of generations allowed for the SQP algorithm was 600, while the function tolerance and the nonlinear constraints tolerance were both set to 10−6. 600 numerical integration points were chosen to perform time integration by virtue of the fourth-order Runge-Kutta method. The number of harmonics retained is 10. A frequency range from 0.2 to 5 Hz was considered, and the vibration displacements covered by this frequency range are to be maximized.
The numerical simulations were performed on a Lenovo notebook computer with Intel Core processors i3-330M of 2.13 GHz and 2 GB DDR3 RAM. For calculating the stable LCO, the proposed method requires 112 seconds to converge, while the CPU time to calculate the unstable LCO is 132 seconds. The optimization results are given in Figure 1. Note that only the individual amplitudes of the odd harmonics have been plotted, since the even harmonics are equal to zero. In addition, it is clear that the first and the third harmonic components are dominant and others are very small.
In order to make the comparisons and to demonstrate the existence of the unstable LCO, direct numerical integration has been performed for the unstable LCO at . Two cases of initial condition were considered, and all simulations were carried out over 1000 periods. For the first case, the initial condition was provided by the HBM solution of the proposed method, and Figure 2(a) displays the temporal evolution of numerical integration solution. As demonstrated in Figure 2(a), the unstable LCO jumps to the stable LCO with high-amplitude after a finite period of time. For the second case, starting from the initial displacement which is less than the amplitude of unstable LCO, the evolution process is given in Figure 2(b). It can be observed in Figure 2(b) that the unstable LCO will go to the stable equilibrium solution, thus confirming the existence of the unstable LCO.
(a) The first case
(b) The second case
Table 1 summarizes the nonlinear equality and inequality constraints that were evaluated at these optimal solutions. In Table 1, denotes the maximum absolute error of the nonlinear algebraic equations of (3), and indicates the Floquet multipliers obtained from the monodromy matrix. Conforming with the Floquet theorem, there is always one Floquet multiplier 1 in Table 1. It can be clearly observed that the nonlinear algebraic equations in (3) and the stability condition in (11) are satisfied, and the largest real part of the Floquet multipliers confirms one of the attained LCOs to be unstable.
Figure 3 depicts the time series and phase portraits of two simultaneously LCOs at , indicating that the main difference between them is the level of the amplitudes, with the stable LCO about four times bigger than the unstable LCO. Overall, the proposed method can be viewed as a root finding algorithm for calculating the LCOs at a given flow speed.
(a)
(b)
The absolute errors in the time domain over one period of vibration between the results of proposed method and that of the time integration method are shown in Figure 4. It is observed from Figure 4 that the discrepancies of the results between the proposed method and the time integration method are quite small.
(a) The unstable LCO
(b) The stable LCO
3.3. Nonlinear Dynamic Response of the Aeroelastic System with Uncertainties
In the following, three cases of uncertainty are investigated by utilizing the proposed method. In each of the considered cases, the flow speed is varied from 2 to 12, and the values of and are included in the range , . The three cases are modeled as follows:
Case 1. Uncertainties come from the cubic plunging stiffness and the flow speed, , , .
Case 2. The cubic pitching stiffness and the flow speed are uncertain, , , .
Case 3. Both the flow speed and the values of cubic plunging and pitching stiffness vary, , , .
With the help of the increment harmonic balance method along with the continuation method, Cases 1 and 2 were studied in a previous paper by Chen et al. [14]. However, for Case 3 all the parameters studied vary simultaneously, and the number of unknown variables is by four () more than the number of the nonlinear harmonic balance equations. The nonlinear harmonic balance equations are the so called underdefined system which cannot be solved using the nonlinear root finding algorithm. Therefore, the combined effects of parameter variations on the peak vibration amplitude cannot be investigated via the increment harmonic balance method in [14].
By using the proposed method, the maximum vibration amplitudes and bifurcation points are found as follows.
(a) Searching the Maximum Vibration Amplitudes. In order to obtain the vibration peaks, the vibration amplitudes of the LCOs are maximized. As explained previously, the unknown parameters that have to be determined are the harmonic coefficients , the vibration frequency , the flow speed , and the cubic plunging stiffness and/or the cubic pitching stiffness of the stable LCOs.
After the MultiStart algorithm stops successfully, numerical results are obtained. The optimization results are reported in Table 2 and Figure 5. Observe in Table 2 that for Cases 1 and 2 the vibration amplitudes reach its maximum at the lower bound of the cubic plunging (pitching) stiffness. However, when and vary simultaneously, a maximum vibration amplitude for Case 3 is found at the upper bound of and at the lower bound of . As for the influences of the flow speed, the maximum vibration amplitudes occur at the highest flow speed considered for both Cases 2 and 3, while for Case 1 the vibration peak locates at the flow speed of . In addition, in Figure 5 analysis of each harmonic contribution notices the emergence of third and fifth harmonics for Case 2.
The nonlinear equality and inequality constraints that were evaluated at these optimal solutions are listed in Table 3. As observed from Table 3, the maximum absolute error of the nonlinear algebraic equations for all these LCOs considered is . Moreover, all the moduli of the Floquet multipliers are less than 1, so the LCOs related to the three cases obtained by the proposed method are stable.
The time series and phase portraits related to these optimal solutions are illustrated in Figure 6. For Case 1, as can be seen from Figure 3 in [14], the upper bound at is at the middle of 0.5 and 0.55, which is the same as the maximum vibration amplitude 0.52619 shown in Figure 6(a). However, for Case 2 the mean value and standard deviation of the amplitude of the LCO at given in Figure 6 and Figure 10 of [14] are 0.6 and 0.04, respectively. The sum of the mean value and standard deviation is 0.64 which is less than the maximum amplitude 0.68357 found by the proposed method. Hence, the advantage of the proposed method is demonstrated.
(a)
(b)
From Figure 6 the combined effect of the cubic plunge and pitch stiffness coefficients and on the maximum vibration amplitude is clearly revealed. As can be seen in Figure 6, the smallest value of the maximum vibration amplitude for the three cases corresponds to Case 3. This indicates that the coexistence of the plunge and pitch nonlinearity may lead to low-amplitude LCO compared to the cases in which structural nonlinearity is applied only in the plunge or pitch DOF.
(b) Finding the Bifurcation Points. For computing the bifurcation points, the stability condition shown in (11) is set as . The vibration displacement is to be maximized. Numerical optimization results are shown in Table 4 and Figure 7. For Case 1, a bifurcation point is detected at . This results is similar to that in Figure 1 from [14], where a turning point between occurs. Therefore, the LCO for Case 1 is not only the bifurcation point but also the turning point. As observed in Figures 2 and 7 of [14], a supercritical bifurcation for Case 2 exists at . The optimization results shown in Table 4 are identical with that reported in [14]. Furthermore, for Case 3 the optimization solution identifies the flow speed of as the bifurcation point.
The nonlinear equality and inequality constraints related to these optimal solutions are given in Table 5. For the three cases, there are always two multipliers equal to 1, which the optimization algorithm intends to find. In addition, it can be seen that the nonlinear equality and inequality constraints are satisfied. Therefore, the proposed method provides a way to detect bifurcation.
The time series and phase planes associated with these optimization solutions are presented in Figure 8. Observe in Figure 8 that a significant magnitude of vibration response on the plunge DOF is observed for Case 1 as the flow speed becomes close to 2.3277. However, for Case 2 a maximum vibration amplitude of 0.00010081 m was experienced at the flow speed of , and the aeroelastic system oscillates with negligible amplitude.
(a)
(b)
With respect to the computational cost of this example, the proposed method spent 208.4375, 195.6094, and 221.6094 seconds of CPU time for the three cases studied here. It should be noted that the conventional continuation method cannot deal with the problem which the uncertainty variables vary simultaneously. However, using the proposed method the uncertainty variables are treated as the optimization variables to investigate the effects of uncertainty and nonlinearity.
4. Conclusions
An efficient method is presented for finding the LCOs and bifurcations in nonlinear aeroelastic system. The proposed method which combines three methods is formulated as a constrained, nonlinear optimization problem. The nonlinear algebraic equations derived from the harmonic balance method are applied to form the nonlinear equality constraints, while the Floquet theory is used to construct the stability constraint inequality and help determine the bifurcation points. Then, the MultiStart algorithm is used to optimize the response amplitude within the specified range of physical parameters.
The efficiency and ability of the proposed method are demonstrated, and the combined effects of multiple nonlinearities and uncertainties on the maximum vibration amplitudes of LCOs and bifurcations of an airfoil are analyzed. It is illustrated that the proposed approach can be used to find the vibration peaks, the bifurcation points, and the multiple solutions. In addition, the presence of multiple nonlinearities may results in low amplitude LCO compared to that structural nonlinearity is applied only in the plunge or pitch DOF. Future investigations will be devoted to extend the presented method to chaotic motion of nonlinear aeroelastic system.
Acknowledgments
The author is grateful to the anonymous referees for their valuable comments. This study has been financially supported by Natural Science Foundation of China (Project no. 10904178).