Abstract
A modified incremental harmonic balance method is presented to analyze the aeroelastic responses of a 2-DOF airfoil aeroelastic system with a nonsmooth structural nonlinearity. The current method, which combines the traditional incremental harmonic balance method and a fast Fourier transform, can be used to obtain the higher-order approximate solution for the aeroelastic responses of a 2-DOF airfoil aeroelastic system with a nonsmooth structural nonlinearity using significantly fewer linearized algebraic equations than the traditional method, and the dominant frequency components of the response can be obtained by a fast Fourier transform of the numerical solution. Thus, periodic solutions can be obtained, and the calculation process can be simplified. Furthermore, the nonsmooth nonlinearity was expanded into a Fourier series. The procedures of the modified incremental harmonic balance method were demonstrated using systems with hysteresis and free play nonlinearities. The modified incremental harmonic balance method was validated by comparing with the numerical solutions. The effect of the number of harmonics on the solution precision as well as the effect of the free-play and stiffness ratio on the response amplitude is discussed.
1. Introduction
Nonlinear phenomena can be frequently encountered in the aerospace industry. It is important to predict the nonlinear aeroelastic characteristic of airfoils [1], which is the reason that several operational aircraft continue to experience limit cycle oscillations (LCOs) within their flight envelopes, often resulting in performance degradation [2]. Therefore, it is of great significance to perform nonlinear flutter analysis.
The harmonic balance (HB) method was first used to analyze the nonlinear aeroelastic responses of structures by Woolston et al. [3] and Shen [4]. This method has been widely used for nonlinear systems, as it transforms the problem into a set of nonlinear algebraic equations using a truncated Fourier series, and is generally applicable for weak and strong nonlinear systems. The classical harmonic balance method, in which the assumed solution contains only one major harmonic in the analysis, can also be called the describing function method and cannot predict the high-order harmonic responses [2, 5]. However, for the most engineering problems, such as aeroelastic problems, the solutions obtained by the classical harmonic balance method are not accurate. Some researchers have proposed improved harmonic balance methods. For the nonlinear aeroelastic analysis of a 2-DOF airfoil with a cubic nonlinearity, a larger number of harmonics can be used to describe the quadratic bifurcation. However, with the increase in the number of harmonics, the derivation of nonlinear terms becomes too complex. The computational cost of the harmonic balance method increases with the square of their order. The higher the order of the problem, the larger the required number of HB coefficients to describe its solution [6]. Vio et al. [7] introduced a search procedure using genetic algorithms to evaluate the coefficients of a harmonic balance solution, showing that the genetic algorithms could provide high-quality initial guesses for the harmonic coefficients. Dimitriadis et al. [8] proposed a harmonic shifting technique to allow a higher-order harmonic balance to accurately estimate both branches of the limit cycles occurring after the second bifurcation. Furthermore, the number of unknowns is larger than the number of nonlinear algebraic equations obtained using the HB method, which cannot be solved directly. For this issue, Lim et al. [9] used Newton’s method to reduce the complexity of the HB method. Chen et al. [10] proposed a new method that translated the problem to a minimization problem and then analyzed the problem by solving the gradient equations. Niu et al. [11] developed a modified harmonic balance method to analyze the nonlinear aeroelasticity of two-degrees-of-freedom airfoils with cubic restoring forces, for which the particle swarm optimization method, with a high calculation efficiency, was used to solve the minimization problem. Evidently, modifications of the HB method are required to achieve more accurate solutions, often resulting in more complex procedures. Berci et al. [12] proposed a novel combined approach whereby the transient response was obtained from the multiple timescales method, in which the asymptotic periodic behavior was corrected using the HB method.
The incremental harmonic balance (IHB) method [13], combining the incremental method and the harmonic balance method in the numerical calculation, was presented by Lau. It has been widely applied to analyze several kinds of nonlinear equations. Lau et al. [14] solved the piecewise-nonlinear equations with linear rigidity and obtained the amplitude-frequency curves as well as the relationship between the harmonic constant and the external excitation frequency. Lau analyzed the dynamic characteristics of a piecewise linear stiffness system and extended it to an asymmetric piecewise linear stiffness system [15, 16]. Shen et al. [17] extended the IHB to analyze the nonlinear dynamics of a spur gear pair, in which the time-varying stiffness and static transmission error were represented by a multiorder harmonic series through a Fourier expansion. They found the general forms of the periodic solutions, which were useful for obtaining solutions with arbitrary precision. Liu et al. [18] investigated the aeroelastic system of an airfoil with nonlinear hysteresis in the pitching degree of freedom using the IHB. They extended the method of undetermined coefficients by expanding the hysteresis nonlinearity into a series to deduce the IHB iterative scheme. Liu et al. [19] studied a sort of elastic mass with asymmetric hysteresis characteristics using IHB. They obtained the analytical linearized algebraic equations of the system, and the coefficients of the algebraic expression were determined by an incremental procedure and an iterative process of the regulated variable based on the algebraic equation. Huai et al. [20] improved the IHB to determine periodic solutions of bilinear hysteretic systems, in which an improved continuation method was used, called the two-point tracing algorithm, where the examination of the stability at the turning point made the calculation more efficient for tracing the amplitude-frequency response. Niu et al. [21] combined the harmonic balance method with the incremental harmonic balance approach to obtain higher-order approximate steady-state solutions for strongly nonlinear systems, which could simplify the calculation process for high-order nonlinear terms.
Most of the studies mentioned above involved systems with smooth nonlinearities, such as quadratic, cubic, or even higher-order nonlinearities, which could be solved using the HB or IHB methods. It is well known that it is easier to handle smooth nonlinearities than nonsmooth ones. Furthermore, as the number of degrees of freedom (DOFs) of the system changes, the difficulty of the analysis and calculations dramatically increases, which is inconvenient for obtaining solutions. Moreover, nonsmooth nonlinearities exist widely in engineering fields, such as free-play and hysteresis nonlinearities, which can be expressed in the form of nonsmooth functions. Not surprisingly, systems with free-play or hysteresis nonlinearities can rarely be solved by the HB or IHB methods due to the difficulty of handling nonsmooth functions [18]. The main reason is that it is very complex to expand nonsmooth functions into the Fourier series. Thus, the development of an effective and general method for the Fourier series expansion of multi-DOF systems with nonsmooth nonlinearities is indispensable.
In this paper, a modified incremental harmonic balance method based on the fast Fourier transform and the traditional incremental harmonic balance method is presented for two-degrees-of-freedom aeroelastic airfoil motion with nonsmooth structural nonlinearities, in which the solution is based on multiple harmonics. After the alternative incremental and iterative processes, the analytic algebraic equations are obtained by applying the proposed IHB, in which the nonsmooth nonlinearity is expanded into a Fourier series. Finally, the procedures of the modified incremental harmonic balance method were demonstrated using systems with hysteresis and free-play nonlinearities. The modified incremental harmonic balance method was validated by comparing with the numerical solutions obtained by the Runge–Kutta method. The effect of the number of harmonics on the solution precision and the effect of the free-play and stiffness ratio on the response amplitude are analyzed.
2. Equations of Motions
A two-dimensional airfoil is shown in Figure 1, which oscillates in pitch and plunge. The model is free to rotate on the X-Z plane and translate in the vertical direction. The plunge displacement is denoted by , and represents the pitch angle about the elastic axis. is the length of the semichord. is the distance between the center of mass of the airfoil section and the elastic axis. is the distance between the elastic axis and the semichord point. and are the plunge stiffness and pitch stiffness, respectively.
Consequently, the aeroelastic equations of the airfoil can be recast in a nondimensional form as follows [22]:where is the nondimensional displacement, is the nondimensional time, is the physical time, , where and are the natural frequencies of the uncoupled plunging and pitching modes, respectively, is the radius of gyration about the elastic axis, is the nonlinear moment, and and are the unsteady aerodynamic lift and moment, respectively, which can be modeled based on the unsteady vortex lattice model. Figure 2 shows the discrete model of the vortex lattice for a two-dimensional airfoil section. To establish the aeroelastic equations, the unsteady lift and the moment are recast into the following form:where is the air density, is the dimensionless flow velocity, is the dimensionless location of the jth vortex, is the dimensionless strength of the bound vortices, and is the total number of bound vortices. Further details can be found in Ref. [22].
Generally, the aerodynamic force based on the unsteady vortex lattice model is referred to as the reduced-order force [23–25]. Using of equations (1) and (2), one can obtain the aeroelastic equation. For brevity, the aeroelastic equation can be rewritten as follows:where , , , , , , and , which is the generalized coordinate vector associated with the aerodynamic force (further information can be found in Ref. [22]). , where and are the number of real and complex conjugate eigenvalues of the aerodynamic system, respectively, and R is the number of reserved aerodynamic modes. The relevant elements in the matrix can be found in Ref. [22]. Furthermore, .
If the nonlinearity is due to backlash in loose or worn control surface hinges, the nonlinearity exhibits free-play. Furthermore, if the friction and backlash must be considered, the nonlinearity is usually in the form of hysteresis [18], which is a nonsmooth function. Figure 3 shows the illustration of nonsmooth function, which can be expressed as follows:where and .
3. Modified Incremental Harmonic Balance Method Based on Fast Fourier Transform
3.1. Traditional Incremental Harmonic Balance Method
The first step in the IHB method is the Newton–Raphson procedure. Equation (3) is transformed by introducing , resulting in the following:where denotes the vibrational frequency. The double () and single () primes represent the second and first derivative with respect to , respectively. denotes the current state of the vibration. The neighboring state can be expressed by adding the corresponding increments as follows:where and are the small terms. Substituting equations (6) into (5), all higher-order terms are neglected, such as and , which is described in Ref. [26]. This procedure is different from the HB. In the HB method, the state can be approximated by a Fourier series, i.e.,.A higher-order harmonic balance refers to cases where . Expressing by the first-order Taylor series approximation yields the following linearized incremental equation:where , which is the error vector. tends to zero as approaches the exact solutions.
It is assumed thatwhere and , where .
Equations (8) and (9) are written in the matrix form as follows:where , .
Substituting equations (10) into (7), and implementing the Galerkin procedure, a set of linearized algebraic equations are obtained as follows:where , , , , , , and and are the first-order Taylor series approximations of the nonlinear forces caused by and , respectively. In this paper, the nonlinear force is only related to the pitch angle. Thus, the elements related to the pitch angle in and are determined below, and the remaining elements are zero. The explicit forms of and can be worked out as follows:where and represents the number of solutions of the equations and in the interval . Here , , and the solutions of and are assumed to be , which are sorted in the ascending order.
In equation (12), and denote the relationship between , , and in the interval , , , , respectively. They can be expressed as follows:where in equation (13) is referred to as equation (8). The above integral can be obtained directly by MATLAB. The matrices and vectors depend upon the initial solutions. Equation (11) contains equations and unknowns, which cannot be determined directly by solving these equations. Generally, one can fix one of the Fourier coefficients to be zero [18], and then equation (11) describes a set of equations in terms of the increments and , which can be solved iteratively.
3.2. Modified Incremental Harmonic Balance Method
The assumed solution given in the previous section (equation (8)), using the traditional incremental harmonic balance method, depends on the number of harmonics and frequencies. On the one hand, if the number of harmonic balance coefficients is large, the dimensions of matrices in equation (11) become very large, and the computational cost increases with the number of harmonics. On the other hand, if the initial frequency is far from the true solution, it will result in an incorrect solution. Thus, the determination of the number of harmonics and frequencies is the key to the precision of the incremental harmonic balance method.
The basic idea of the modified incremental harmonic balance method is that a fast Fourier transform is carried out to extract the dominant frequency components in the response of the system using a numerical method, which determines the harmonic number and the dominant frequency of the approximate solution. The dominant frequency components are then applied to the incremental harmonic balance method to solve the linearized equations, which effectively determines the approximate solution of the nonlinear system and improves the solution precision.
To reduce the dimensions of matrices in equation (11), the assumed solution is written as follows:where is determined by a fast Fourier transform. , , , and denote the coefficients of the harmonics.
According to the steps discussed above, the solutions for an aeroelastic system with a nonsmooth structural nonlinearities can be obtained. Figure 4 shows the flowchart of the proposed method.
3.3. Comparison with Traditional IHB Method
There are two key differences between the two IHM methods:(1)If the solutions of the multi-DOF system contain high-order harmonics, the number of linearized algebraic equations in terms of and that must be solved using the traditional IHB method will be very large, and the solution process will be time-consuming. However, the modified IHB method can obtain the higher-order approximate solution with significantly fewer linearized algebraic equations for multi-DOF systems because it uses equations (14) and (15), in place of equations (8) and (9). The parameters in equations (14) and (15) can be determined by a fast Fourier transform. The number of the linearized algebraic equations, i.e., equation (11), that must be solved using the current method is smaller than that by the traditional IHB.(2)When the nonlinear term is nonsmooth, its Fourier series expansion in the modified IHB method is much simpler, as an effective and general method for the Fourier series expansion of nonsmooth nonlinearities was developed. The explicit expressions of the first-order Taylor series approximations of the nonsmooth nonlinear term were derived, i.e., equations (12) and (13), allowing the nonsmooth nonlinearity to be expanded as a Fourier series. This approach can provide guidance for representing other nonlinearities.
4. Validation and Result Analysis
A description of a two-dimensional airfoil and the corresponding parameters can be found elsewhere [28]. In structural dynamics, the dynamic model of a complex structure is often reduced to a simple one with a few degrees of freedom by modal reduction. A similar procedure has been applied to simplify the models of unsteady aerodynamics [27]. It may be possible to predict the unsteady aerodynamic response of a complex system over a wide frequency range based on the reduced-order model. The modal reduction, which can extract the most important aerodynamic modes, can be applied to reduce the order of the aerodynamic equations [22]. In this study, two aerodynamic modes are reserved, and are the generalized coordinate vectors associated with the aerodynamic force. The nondimensional linear flutter velocity for the linear aeroelastic system (i.e., ) was . The linear flutter velocity was nearly the same as that reported previously [28]. In Ref. [28], the responses are obtained by using Runge–Kutta and are stable periodic solutions. In this section, the modified incremental harmonic balance method is applied to directly simulate the responses of the 2-DOF airfoil aeroelastic systems with nonsmooth structural nonlinearities.
4.1. Validation
The incremental harmonic balance method was mainly used to determine the periodic motion of nonlinear systems. To demonstrate that the modified incremental harmonic method can be used to determine the response of the aeroelastic system with a nonsmooth nonlinearity, the plunge and pitch responses of the two-dimensional airfoil were obtained through two examples: one with a hysteresis nonlinearity, and one with a free-play nonlinearity.
For the example with a hysteresis nonlinearity, the following parameters were employed: , , (i.e., ). All the elements were zero in the initial condition, except that (approximately ). The numerical solutions were obtained by the Runge–Kutta method, after which a fast Fourier transform was performed to extract the dominant frequency components in the responses, which determined the number of harmonics and the frequency of the approximate solution. The numerical simulation and the frequency spectrum are shown in Figure 5.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
As shown in Figure 5, 0.4524 was the key frequency component of the plunge response. The other responses consisted mainly of two frequency components, i.e., 0.4524 and 1.357, where 1.357 is three times larger than 0.4524. First, supposing that the assumed solution of the aeroelastic equation with a nonsmooth structural nonlinearity contained 1, , terms, the solution was obtained according to the method described in Section 3. Next, because the second frequency response in Figure 5 was three times larger than the first frequency response, a solution containing 1, , , , and terms was assumed. With this approach, the effect of harmonics on the accuracy of the solutions could be determined. The two different assumed solutions yielded the following approximations. When the assumed solution contained 1, , and terms, the following approximation was obtained:
When the assumed solution contained 1, , , , and , the following approximations were obtained:
The numerical solutions were obtained by the Runge–Kutta method. Figures 6 and 7 show comparisons of the periodic solution obtained using the current method and that obtained using the numerical method. Differences between the approximation obtained using the assumed solution containing 1, , and terms and the numerical method were evident. However, the approximation obtained using the assumed solution containing 1, , , , and terms showed excellent agreement with the numerical solution obtained using the Runge–Kutta method. For the plunge response, the approximation obtained using both assumed solution forms agreed well with the numerical solution. Because 0.4524 was the key frequency component for the plunge response in Figure 5, the approximation containing 1, , and terms could cover all the dominant frequencies of the plunge response, which could accurately reflect the system response. The other responses consisted mainly of frequency components of 0.4524 and 1.357. Thus, the approximation containing 1, , , , and terms was needed to obtain an accurate result.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
For the free-play nonlinearity, the following parameters were employed: . , and (i.e., ). All the elements were zero initially, except for (approximately ). The numerical solutions obtained by the Runge–Kutta method and the fast Fourier transform of the responses are shown in Figure 8.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
As shown in Figure 8, the responses consisted mainly of 0.4398 and 1.307 frequency components, where 1.307 was approximately three times greater than 0.4398. Approximations containing 1, , and terms and 1, , , , and terms were used to represent the first and second frequency components, respectively. Using these approximations, the effects of the harmonics on the accuracy of the solutions were determined. When the responses contained higher-order harmonics, the number of terms in the linearized algebraic equations obtained by the current method was significantly less than that obtained by the traditional IHB. For instance, the expansion in the traditional IHB was replaced with in the current method. For the current system with nonsmooth structural nonlinearities, the number of the linearized algebraic equations, i.e., equation (11), by the traditional IHB was 28. However, only 20 equations were required for the proposed method. Thus, solutions with higher accuracy and fewer linearized algebraic equations could be obtained by the proposed method.
For the solution containing 1, , and terms, the approximate solutions were as follows:
For the solution containing 1, , , , and terms, the approximate solution were as follows:
The numerical solutions were obtained using the Runge–Kutta method. Figures 9 and 10 show the comparisons of the periodic solutions of the current method and the numerical method. The approximation containing 1, , and terms and the numerical solutions have a marginal deviation. However, the approximation containing 1, , , , and terms showed excellent agreement with the numerical solutions.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4.2. Parametric Analysis
To analyze the effects of and on the response of the nonlinear aeroelastic system, the response curves for different values of and are discussed. The approximation containing 1, , , , and terms was adopted with the proposed method in the following simulations.
4.2.1. Effect of and on the Responses of the Aeroelastic System with Hysteresis
The effects of and on the response amplitude curves of the aeroelastic system with hysteresis were investigated. First, a fixed value of was set, and was varied to examine its effect on the response. Next, a fixed value of was set, and was varied. The comparisons of amplitude are shown in Figures 11–14.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
As shown in Figures 11 and 12, the approximation containing 1, , , , and terms was in excellent agreement with the numerical solutions. Furthermore, the nondimensional velocity range over which the solution was periodic decreased as increased. As increased, the nondimensional velocity over which the solution was periodic increased. For example, when was 0.2, the nondimensional velocity at which the solution was periodic was 1.5. When was 0.5, the nondimensional velocity at which the solution was periodic was 1.7.
As shown in Figures 13 and 14, the solution obtained by the proposed method was in good agreement with the numerical solution. As the nondimensional velocity varied from 1.6 to 1.9, there were periodic solutions for different values of . The response amplitudes increased with for the same nondimensional velocity. For instance, when varied from to , the plunge and pitch amplitude varied from 0.007623 to 0.1526 and from 0.1526 to 0.2175, respectively.
4.2.2. Effect of on the Responses of the Aeroelastic System with Free Play
The effect of in the aeroelastic system with free play on the response amplitude curve was investigated, and the response amplitude curves are shown in Figures 15 and 16.
(a)
(b)
(a)
(b)
When was set to zero, the aeroelastic system with hysteresis becomes a system with free play. As shown in Figures 15 and 16, the nondimensional velocity range over which the solution was periodic increased with . With equals , the solution was periodic when the nondimensional velocity approached 1.6, which was close to the linear flutter velocity. Thus, a system with less free play can be considered a system with a weak nonlinearity or even a linear system. Similarly, the response amplitudes increased with .
As shown in Figures 14 and 16, for the system with hysteresis, the solution was periodic when the nondimensional velocity approached 1.6. For the free play system, the solution was periodic when the nondimensional velocity approached 1.1. Thus, the hysteresis nonlinearity could enhance the flutter velocity. The response amplitude of the system with a free play nonlinearity was slightly larger than that of the system with hysteresis. The reason was that the effective stiffness of the free play nonlinearity was less than that of the hysteresis nonlinearity.
5. Conclusions
In this paper, a modified incremental harmonic balance method is presented for an aeroealstic system with a nonsmooth structural nonlinearity. The procedure for the modified incremental harmonic balance method was demonstrated for systems with hysteresis and free play nonlinearities. The validity of the modified incremental harmonic balance method was demonstrated by comparing with the numerical solutions. In addition, the influence of the parameters on the nonlinear aeroelasticity was also studied. The following conclusions were drawn:(1)A modified incremental harmonic balance method was presented to analyze the responses of a 2-DOF airfoil aeroelastic system with a nonsmooth structural nonlinearity. Not only could the periodic solutions be obtained, but also the calculation process was simplified.(2)The proposed approach can be applied in other nonsmooth cases, especially those arising in aeroelastics. The application of the incremental harmonic balance method was extended.(3)For a given , the response amplitudes increased with . For a given nondimensional velocity, the hysteresis nonlinearity could lead to a lower response amplitude than the free play nonlinearity did.
6. Discussion
Since the response of the aeroelastic model in this paper is stable, the method in this paper is limited to the stable periodic solution. Going forward, more studies are required, including the unstable periodic solution. The homology method may obtain the unstable periodic solution. Therefore, it may be suitable for the response of aeroelastic system with nonsmooth structural nonlinearities by integration of the incremental harmonic balance method and the homotopy method. The application of the incremental harmonic balance method may be expanded.
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 are no conflicts of interest regarding the publication of this paper.