A Modified Incremental Harmonic Balance Method for 2-DOF Airfoil Aeroelastic Systems with Nonsmooth Structural Nonlinearities

A modiﬁed 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 signiﬁcantly 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 simpliﬁed. Furthermore, the nonsmooth nonlinearity was expanded into a Fourier series. The procedures of the modiﬁed incremental harmonic balance method were demonstrated using systems with hysteresis and free play nonlinearities. The modiﬁed incremental harmonic balance method was validated by comparing with the numerical solutions. The eﬀect of the number of harmonics on the solution precision as well as the eﬀect of the free-play and stiﬀness ratio on the response amplitude is discussed.


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]. erefore, it is of great significance to perform nonlinear flutter analysis. e harmonic balance (HB) method was first used to analyze the nonlinear aeroelastic responses of structures by Woolston et al. [3] and Shen [4]. is 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. e 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 highorder 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. e computational cost of the harmonic balance method increases with the square of their order. e 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 highquality 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. e 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 piecewisenonlinear 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. ey 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. ey 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. ey 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]. e main reason is that it is very complex to expand nonsmooth functions into the Fourier series. us, 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. e modified incremental harmonic balance method was validated by comparing with the numerical solutions obtained by the Runge-Kutta method. e 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.

Equations of Motions
A two-dimensional airfoil is shown in Figure 1, which oscillates in pitch and plunge. e model is free to rotate on the X-Z plane and translate in the vertical direction. e plunge displacement is denoted by h, and α represents the pitch angle about the elastic axis. b is the length of the semichord. x a b is the distance between the center of mass of the airfoil section and the elastic axis. ab is the distance between the elastic axis and the semichord point. k h and k α 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]: 2 Mathematical Problems in Engineering where ω h and ω a are the natural frequencies of the uncoupled plunging and pitching modes, respectively, r α is the radius of gyration about the elastic axis, M n (α(t)) is the nonlinear moment, and L and M ea 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 twodimensional airfoil section. To establish the aeroelastic equations, the unsteady lift L and the moment M ea are recast into the following form: where ρ a is the air density, V non is the dimensionless flow velocity, κ j is the dimensionless location of the jth vortex, Γ is the dimensionless strength of the bound vortices, and M 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][24][25]. Using of equations (1) and (2), one can obtain the aeroelastic equation. For brevity, the aeroelastic equation can be rewritten as follows: where q � q 1 · · · q m q m+1 q m+2 · · · q m+n T , which is the generalized coordinate vector associated with the aerodynamic force (further information can be found in Ref. [22]).
, where m and n are the number of real and complex conjugate eigenvalues of the aerodynamic system, respectively, and R is the number of reserved aerodynamic modes. e 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:

Traditional Incremental Harmonic Balance
Method. e first step in the IHB method is the Newton-Raphson procedure. Equation (3) is transformed by introducing τ � ωt, resulting in the following: Mathematical Problems in Engineering where ω denotes the vibrational frequency. e double ( ″ ) and single ( ′ ) primes represent the second and first derivative with respect to τ, respectively. (X 0 , ω 0 ) denotes the current state of the vibration. e neighboring state can be expressed by adding the corresponding increments as follows: where ΔX and Δω are the small terms. Substituting equations (6) into (5), all higher-order terms are neglected, such as Δω 2 M € X 0 and 2ω 0 ΔωMΔ € X, which is described in Ref. [26]. is procedure is different from the HB. In the HB method, the state can be approximated by a Fourier series, i.e.,x(t) � x 0 + N n�1 [x 2n− 1 cos(nωt) + x 2n sin(nωt)]. A higher-order harmonic balance refers to cases where N > 1. Expressing F(X) by the first-order Taylor series approximation yields the following linearized incremental equation: where R � − (ω 2 0 MX 2 0 + ω 0 DX 0 + KX 0 ) − F(X 0 ), which is the error vector. R tends to zero as (X 0 , ω 0 ) approaches the exact solutions.
It is assumed that where C s � 1cos τ cos 2 τ· · ·cos(N c − 1)τ sin τ sin 2 τ· · · sin N s τ] and Equations (8) and (9) are written in the matrix form as follows: Substituting equations (10) into (7), and implementing the Galerkin procedure, a set of linearized algebraic equations are obtained as follows: Airfoil Weak Figure 2: Discrete model of vortex lattice for a two-dimensional airfoil section [22].
0 S T KSdτ, and K NL and R NL are the first-order Taylor series approximations of the nonlinear forces caused by F ′ (X 0 ) and F(X 0 ), respectively. In this paper, the nonlinear force is only related to the pitch angle. us, the elements related to the pitch angle in K NL and R NL are determined below, and the remaining elements are zero. e explicit forms of K NL and R NL can be worked out as follows: where C stmp � 1 cos θ cos 2 θ · · · cos(N c − 1)θ sin θ sin 2 θ · · ·sin N s θ] and P represents the number of solutions of the equations |α(τ)| � δ 0 and |α(τ)| � δ in the interval (0, 2π). Here θ 0 � 0, θ P+1 � 2π, and the solutions of |α(τ)| � δ 0 and |α(τ)| � δ are assumed to be θ i (i � 1, . . . , P), which are sorted in the ascending order.
In equation (12), sgn a and sgn b denote the relationship between α(τ), δ 0 , and δ in the interval θ 0 θ 1 , θ 1 θ 2 , · · ·, θ P θ P+1 , respectively. ey can be expressed as follows: where α(τ) in equation (13) is referred to as equation (8). e above integral can be obtained directly by MATLAB. e matrices and vectors depend upon the initial solutions. Equation (11) contains (2 + R)(2N + 1) equations and (2 + R)(2N + 1) + 1 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 ΔA, which can be solved iteratively.

Modified Incremental Harmonic Balance Method.
e 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. us, the determination of the number of harmonics and frequencies is the key to the precision of the incremental harmonic balance method. e 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.
e 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 c k is determined by a fast Fourier transform. a jk , b jk , Δa jk , and Δb jk 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.

Comparison with Traditional IHB Method.
ere 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 ΔA and Δω that must be solved using the traditional IHB method will be very large, and the solution process will be timeconsuming. 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). e parameters in equations (14) and (15) can be determined by a fast Fourier transform. e 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. e explicit expressions of the first-order Taylor series approximations of the Mathematical Problems in Engineering nonsmooth nonlinear term were derived, i.e., equations (12) and (13), allowing the nonsmooth nonlinearity to be expanded as a Fourier series. is approach can provide guidance for representing other nonlinearities.

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. e 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 q 1 , q 2 are the generalized coordinate vectors associated with the aerodynamic force. e nondimensional linear flutter velocity for the linear aeroelastic system (i.e., M n (α(t)) � 0) was V non � 2.0 l. e 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.

Validation.
e 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: V non � 1.5, 2δ � 0.0174 rad, δ 0 � 1.5δ (i.e., k 0 � 0.2k α ). All the elements were zero in the initial condition, except that α(0) � 0.052 rad (approximately 3 ∘ ). e 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. e numerical simulation and the frequency spectrum are shown in Figure 5.
As shown in Figure 5  Mathematical Problems in Engineering equation with a nonsmooth structural nonlinearity contained 1, cos(τ), sin(τ) 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, cos(τ), cos(3τ), sin(τ), and sin(3τ) terms was assumed. With this approach, the effect of harmonics on the accuracy of the solutions could be determined. e two different assumed solutions yielded the following approximations. When the assumed solution contained 1, cos(τ), and sin(τ) terms, the following approximation was obtained: When the assumed solution contained 1, cos(τ), cos(3τ), sin(τ), and sin(3τ), the following approximations were obtained: e 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, cos(τ), and sin(τ) terms and the numerical method were evident. However, the approximation obtained using the assumed solution  Mathematical Problems in Engineering containing 1, cos(τ), cos(3τ), sin(τ), and sin(3τ) 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, cos(τ), and sin(τ) terms could cover all the dominant frequencies of the plunge response, which could accurately reflect the system response. e other responses consisted mainly of frequency components of 0.4524 and 1.357. us, the approximation containing 1, cos(τ), cos(3τ), sin(τ), and sin(3τ) terms was needed to obtain an accurate result.
For the free-play nonlinearity, the following parameters were employed: V non � 1.1. δ � 0.0174 rad, and δ 0 � δ (i.e., k 0 � 0). All the elements were zero initially, except for α(0) � 0.052 rad (approximately 3 ∘ ). e numerical solutions obtained by the Runge-Kutta method and the fast Fourier transform of the responses are shown in Figure 8.  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, cos(τ), and sin(τ) terms and 1, cos(τ), cos(3τ), sin(τ), and sin(3τ) 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 x j0 � 1 + a j1 cos(τ) + a j2 cos(2τ) + a j3 cos(3τ) + b j1 sin(τ) + b j2 sin(2τ) + b j3 sin(3τ) in the traditional IHB was replaced with x j0 � 1 + a j1 cos(τ) + a j3 cos(3τ)+ b j1 sin(τ) + b j3 sin(3τ) 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.
us, solutions For the solution containing 1, cos(τ), cos(3τ), sin(τ), and sin(3τ) terms, the approximate solution were as follows: e 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. e approximation containing 1, cos(τ), and sin(τ) terms and the numerical solutions have a marginal deviation. However, the approximation containing 1, cos(τ), cos(3τ), sin(τ), and sin(3τ) terms showed excellent agreement with the numerical solutions.

Parametric Analysis.
To analyze the effects of k 0 and δ on the response of the nonlinear aeroelastic system, the response curves for different values of k 0 and δ are discussed.

Effect of k 0 and δ on the Responses of the Aeroelastic System with Hysteresis.
e effects of k 0 and δ on the response amplitude curves of the aeroelastic system with hysteresis were investigated. First, a fixed value of δ � 0.0174/2rad was set, and k 0 was varied to examine its effect on the response. Next, a fixed value of k 0 � 0.2 was set, and δ was varied. e comparisons of amplitude are shown in Figures 11-14.
As shown in Figures 11 and 12, the approximation containing 1, cos(τ), cos(3τ), sin(τ), and sin(3τ) terms was in excellent agreement with the numerical solutions. Furthermore, the nondimensional velocity range over which the solution was periodic decreased as k 0 increased. As k 0 increased, the nondimensional velocity over which the solution was periodic increased. For example, when k 0 was 0.2, the nondimensional velocity at which the solution was periodic was 1.5. When k 0 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 δ. e response amplitudes increased with δ for the same nondimensional velocity. For instance, when δ varied from 0.00174/2 rad to 0.0174 rad, the plunge and pitch amplitude varied from 0.007623 to 0.1526 and from 0.1526 to 0.2175, respectively.

Effect of δ on the Responses of the Aeroelastic System with Free Play.
e 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.
When k 0 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 0.00174/2 rad, the solution was periodic when the nondimensional velocity approached 1.6, which was close to the linear flutter velocity. us, 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.
us, the hysteresis nonlinearity could enhance the flutter velocity. e response amplitude of the system with a free play nonlinearity was slightly larger than that of the system with hysteresis. e reason was that the effective stiffness of the free play nonlinearity was less than that of the hysteresis nonlinearity.

Conclusions
In this paper, a modified incremental harmonic balance method is presented for an aeroealstic system with a nonsmooth structural nonlinearity. e procedure for the modified incremental harmonic balance method was demonstrated for systems with hysteresis and free play nonlinearities. e 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. e 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) e proposed approach can be applied in other nonsmooth cases, especially those arising in aeroelastics. e application of the incremental harmonic balance method was extended. (3) For a given k 0 , 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.

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. e homology method may obtain the unstable periodic solution. erefore, 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. e application of the incremental harmonic balance method may be expanded.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.