Highly Accurate Solution of Limit Cycle Oscillation of an Airfoil in Subsonic Flow

The homotopy analysis method (HAM) is employed to propose a highly accurate technique for solving strongly nonlinear aeroelastic systems of airfoils in subsonic flow. The frequencies and amplitudes of limit cycle oscillations (LCOs) arising in the considered systems are expanded as series of an embedding parameter. A series of algebraic equations are then derived, which determine the coefficients of the series. Importantly, all these equations are linear except the first one. Using some routine procedures to deduce these equations, an obstacle would arise in expanding some fractional functions as series in the embedding parameter. To this end, an approach is proposed for the expansion of fractional function. This provides us with a simple yet efficient iteration scheme to seek very-high-order approximations. Numerical examples show that the HAM solutions are obtained very precisely. At the same time, the CPU time needed can be significantly reduced by using the presented approach rather than by the usual procedure in expanding fractional functions.


Introduction
Predicting amplitude and frequency of flutter oscillations of an airfoil via analytical and/or semianalytical techniques has been an active area of research for many years.The describing function technique [1], sometimes referred to as the harmonic balance (HB) or as linearization method, is a widely used method for obtaining an equivalent linear system such that traditional linear aeroelastic methods of analysis can then be employed [2,3].According to the number of considered harmonics, the HB method is called HB1 method when only the first harmonic is included, otherwise as the high-dimensional HB method.Lee et al. [4] studied the aeroelastic system by considering two dominant harmonics and by an improved HB1 method, respectively.Recently, the high-dimensional HB method was further improved to investigate the aeroelastic motions of an airfoil [5,6].Essentially, the incremental harmonic balance (IHB) method is a semianalytical method for nonlinear dynamic systems.It was used by Shahrzad and Mahzoon [7] and Cai et al. [8], respectively, to predict the amplitudes and frequencies of the LCOs of an airfoil in steady impressible flow.Recently, Chung et al. proposed a new incremental method and applied it to solve aeroelastic problems with freeplay [9] and hysteresis [10] structural nonlinearities, respectively.In addition, the center manifold theory, originally developed to qualitatively analyze nonlinear vibrations, was employed to obtain the approximations of airfoil LCOs [11,12].
The approximations obtained by HB1 method are relatively accurate for low wind speeds.However, the errors become larger and larger as the speed increases.In some nonlinear flutter cases, the HB1 method may cease to be valid.In principle, the high-dimensional HB method and the IHB method can give approximate solutions with any desired accuracy as long as enough harmonics are taken into account.Unfortunately, however, it becomes more and more difficult to implement either one of them when the number of considered harmonics increases.Likewise, using the center manifold theory can provide us with satisfactory approximations for the LCOs only in a small range of Advances in Acoustics and Vibration bifurcation values.When far away from the bifurcation points, results loose accuracy significantly or even become completely incorrect [13].Thus, it is necessary to develop new easier-to-use methods which can guarantee accuracy for high flow speeds and in more flutter cases, for example, weakly and strongly nonlinear systems.
Over the past decades, Liao developed the homotopy analysis method (HAM), which does not require small parameters and thus can be applied to solve nonlinear problems without small or large parameters [14][15][16].The main procedure is to construct a class of deformation equations in a quite general form by introducing an auxiliary parameter.Through these equations, nonlinear problems can be transformed into a series of linear subproblems, which can be solved much more easily step by step.Recently, the HAM has been used in various nonlinear problems [17][18][19][20][21].
In this study, the HAM is employed to propose an efficient and highly accurate approach for nonlinear aeroelastic motions of an airfoil.A major obstacle is met when deducing the high-order deformation equations, because the Taylor expansion of fractal functions is rather cumbersome.An approach is proposed to deal with this problem.This simple yet efficient method ensures an excellent efficiency of the HAM; hence, highly accurate solutions can be easily obtained for both weakly and strongly nonlinear aeroelastic systems.

Equations of Motions
The physical model shown in Figure 1 is a two-dimensional airfoil, oscillating in pitch and plunge, which has been employed by many authors.The pitch angle about the elastic axis is denoted by α, positive with the nose up; the plunge deflection is denoted by h, positive in the downward direction.The elastic axis is located at a distance a h b from the midchord, while the mass center is located at a distance x α b from the elastic axis.Both distances are positive when measured towards the trailing edge of the airfoil.
For cubic restoring forces with subsonic aerodynamics, the coupled equations for the airfoil in nondimensional form can be written as follows: where the superscript denotes the differentiation with respect to the nondimensional time t , defined as t = Ut 1 /b, and t 1 is the real time.ξ =h /b is the nondimensional plunge displacement; η is the coefficient of cubic pitching stiffness; U * ←− is a nondimensional flow velocity defined as U * = U/(bω α ), and ω is given by ω = ω ξ /ω α , where ω ξ and ω α are the natural frequencies of the uncoupled plunging and pitching modes, respectively; ζ ξ and ζ α ψ are the damping ratios; r α is the radius of gyration about the elastic axis.P(t) and Q(t) are the externally applied force and moment, m is the airfoil mass per unit length and μ ψ is the airfoil-air mass ratio.C L (t) and C M (t) are the lift and pitching moment coefficients, respectively.For an incompressible flow, the expressions for C L (t) and C M (t) are given by where the Wagner function φ(t) is given by Jone's approximation, φ(t) = 1 − ψ 1 e −ε1t − ψ 2 e −ε2t with the constants as ψ 1 = 0.165, ψ 2 = 0.335, ε 1 = 0.0455, and ε 2 = 0.3.Due to the existence of the integral terms in (3), (1) is a system of integrodifferential equations.In practice, the integral and the nonlinear terms make it difficult to analytically study the dynamic behavior of the system.In order to eliminate the integral terms, Lee et al. [4][5][6] introduced the following four new variables ( The coefficients c 0 , c 1 , . . ., c 10 ; d 0 , d 1 , . . ., d 10 are given in the appendix, f (t) and g(t) are functions depending on initial conditions, Wagner's function, and the forcing terms.The nonlinear restoring forces, G(ξ) and M(α), are expressed as G(ξ) = γξ 3 and M(α) = ηα 3 , respectively, with γ and η as coefficients.By introducing a variable vector X = (x 1 , x 2 , . . ., x 8 ) T , where the superscript " T " denotes the transpose of a matrix, with x 1 = α, x 2 = α, x 3 = ξ, x 4 = ξ, x 5 = w 1 , x 6 = w 2 , x 7 = w 3 , and x 8 = w 4 , the coupled equations given in (5) can be written as a set of eight first-order ordinary differential equations written in vector form This approach allows existing methods suitable for the study of ordinary differential equations to be used in the analysis.

Homotopy Analysis Method
It is assumed that there is no external forces, that is, Q(t) = P(t) = 0 in (1).For large values of t when transients are damped out and steady solutions are obtained, we can let f (t) = g(t) = 0.Then, (5) can be rewritten in vector form as where and F(x) = 0 d 10 ηα 3 ] T .Firstly, introduce a new time scale where ω denotes the frequency of the LCO.Then, (7) becomes where the superscript denotes the differentiation with respect to τ. Considering that LCOs are independent of initial conditions, one can adopt the following initial conditions: The LCOs of system ( 10), ( 11) are periodic motions with frequency ω; thus, x can be expressed in a Fourier series where c k , s k are the coefficients in 2 × 1 vector form.Let a 0 , h 0 , ω 0 , β 0 , and x 0 (τ) denote the initial approximations of a, h, ω, β, and x(τ), respectively.Due to solution expression (12) and initial conditions (11), the initial guess of solution can be chosen as The homotopy analysis method is based on such continuous variations, A(p), H(p), Ω(p), B(p), and u(τ, p), that, as the embedding parameter p increases from 0 to 1, u(τ, p) varies from the initial guess x 0 (τ) to the exact solution, so do A(p), H(p), Ω(p), B(p) from the initial approximations a 0 , h 0 , ω 0 , β 0 to a, h, ω, β, respectively.
Based on (12), one may choose the linear auxiliary operator as Thus One may define the nonlinear operator according to (10), Using the two operators, a family of equations can then be constructed as subject to the initial conditions where the auxiliary parameter λ is a nonzero constant.Equations (17) and ( 18) are called the zeroth-order deformation equation.

Advances in Acoustics and Vibration
When p = 0, ( 17) and ( 18) have the solution When p = 1, they are exactly the same as (10) and (11) provided that Expand A(p), H(p), Ω(p), B(p), and u(τ, p) as the series As long as the parameter λ is properly chosen, all of these series are convergent at p = 1.Then, the nth-order HAM solutions can be given as Substituting ( 21) into ( 17) and (18), differentiating (17) and ( 18) k times, dividing the differentiations by k! and then letting p = 0, one can obtain the kth-order (k ≥ 1) deformation equation subject to the initial conditions where Due to the rule of solution expression and the linear operator L, the right hand side of (23) should not contain the first harmonics sin τ and cos τ, because they can result in the so-called secular terms as τ cos τ and τ sin τ, respectively.To this end, let Solving ( 27), a k , h k , ω k , and β k are determined step by step as k increases.Note that when k = 0, R k+1 (τ) is essentially the right hand side of (10) with x = x 0 , and the integrations in essence correspond to a harmonic balancing procedure.Therefore, ( 27) is actually the algebraic equation deduced by the HB1 method.It is nonlinear and independent upon λ.The solutions of a 0 , h 0 , ω 0 , and β 0 can be obtained by using the Newton-Raphson method.Importantly, (27) is always linear as k ≥ 1, which implies it is rather easy to obtain highorder approximations [22].

Expansion of Fractional Functions
A key procedure of implementing the HAM is to deduce the high-order deformation equation, that is, to obtain R k (τ) in our study.In most literature about the HAM, authors suggest differentiating the zeroth-order deformation equations (i.e., (17) and (18) in this paper) k times, dividing them by k!, and then setting p = 0.This kind of approach is based on the classical theories of the Taylor series.In our study, however, using this method to expand CW(u(τ, p), Ω(p)) will cost a large amount of computational resources.For example, substitution of For large values of t, the second term in (28) approaches to zero and is neglected since only steady solutions (LCOs) are taken into account.Therefore, the integrations W(cos(iωt)) and W(sin(iωt)) can be expressed as follows, respectively: where  the expression of ∂ k W/∂p k becomes more and more complex as k increases, which makes it quite tough to deduce high-order deformation equations.
We take 1/[ε 2 1 +( n i=0 ω i p i ) 2 ] as an illustrative example to propose a means for expanding fractional functions.First of all, denote the denominator as where Taking the nth-order Taylor series of 1/( n k=0 σ k p k ) as n k=0 θ k p k , then one has Rewrite (31) as Equating the coefficients of p k results in Interestingly, (33) is always linear.That means it is rather easy to determine θ k if σ i are all known, i = 0, 1, 2, . . ., k.

Main Results.
The system parameters under consideration are μ = 100, r α = 0.5, a h = −0.5, ζ α = ζ ξ = 0, ω = 0.25, x α = 0.25, γ = 0, and η = 80.Numerical solutions of ( 6) can be obtained by the fourth-order Runge-Kutta method.Without special statement, the numerical solutions in this paper are obtained  subject to the initial conditions as α(0) = 1 • and α(0 Using analytical techniques developed for nonlinear dynamical systems, the linear flutter speed is found at U * = U * L = 6.0385 [4,5].As U * increases beyond U * L , LCO arises, and thus U * L is a Hopf bifurcation point.Note that the flutter boundary U * L is independent of the nonlinear coefficient η. Lee et al. [4] found a secondary Hopf bifurcation as U * increases further, where a jump of the amplitudes is detected (see Figures 2 and 3).Liu et al. [6] used the high-dimensional HB method to study the secondary Hopf bifurcation and found that to capture the secondary bifurcation, as many as 9 (or 5 dominant) harmonics have to be considered.
In the proposed method, the zeroth-order HAM approximation is essentially given by the HB1 method.The higherorder approximations only contribute a higher precision.Thus, it is not capable of detecting the second bifurcation at the present state.Even so, validity and high efficiency of the proposed method can be observed when U * is in [U * L , 2U * L ] or so.Figures 2 and 3 show the comparisons of the 50th-order HAM solutions with the HB1 and the numerical results.The HAM solutions are almost the same as the numerical ones, while the differences of the HB1 results increase rapidly with increasing U * .
The HAM approximation is based on the first HB method, because the first HAM approximation is the HB1 solution.Since the HB1 method is incapable of tracking the LCOs when U * is larger than the secondary bifurcation value, about 2U * L , so is the presented approach, as shown in Figures 2 and 3.
Figures 4 and 5 show the phase planes of LCOs and the time history responses of the nonlinear aeroelastic system, respectively.Again, the accuracy of the HAM solution can be demonstrated.Even though the phase plane is very complex, for example, the pitch LCO, the HAM is still capable of tracking it.Note that the numerical solution plotted in Figure 5 is obtained using the fourth-order Runge-Kutta method with initial values given by the HAM solution.
More precisely, the 120th-order HAM solutions shown in Table 1 are compared with the numerical ones.Excellent agreement can also be observed.The higher the order the HAM approximations are obtained to, the more accurate the solution is.For example, the relative difference between the 120th-order HAM solution and the numerical one is less than 0.001%.As shown in Figure 6, the residues of ( 6) with HAM solutions converge rapidly to 0. The absolute errors of residues with respect to the 40th-order, 80th-order, and 120th-order HAM solutions are at the order of 10 −8 , 10 −12 , and 10 −16 , respectively.Furthermore, as n >120, a n , h n , ω n , and β n are all small quantities compared with 10 −16 .Roughly speaking, the 120th-order HAM solution can be considered to be correct to 15 decimal places.Note that it is tough to obtain such a highly accurate solution using some numerical techniques, including the RK method.
Very interestingly, it is found that ω i is independent of η, while a i , β i, and h i are in inverse proportion to √ η.Thus, the convergence of series ( 21) is independent of η, and the proposed method can work for both weakly and strongly nonlinear problems.Furthermore, it is proved that the frequency of the LCOs of aeroelastic system (5) is independent of η, while the amplitudes are in inverse proportion to √ η.

Choosing the Auxiliary Parameter.
The HAM series are dependent upon the auxiliary parameter λ.For the choice of the value of λ, one should think about two aspects, that is, whether the series converge and the convergent rate.Liao [16] suggested a technique via plotting the curves of the attained HAM solutions versus different values of λ, namely, the λ-curves.From Table 1, one can assume the angular frequency of system (1) with U * = 1.5U * L as = 0.07756360647090.Denote the discrepancy between the nthorder HAM frequency solution and as e(n) = ( n k=0 ω k ) − .Figure 7 shows the λ-curves with respect to ω. Considering that the longitudinal coordinate refers to the logarithm of |e(n)|, one knows the HAM solutions attained with λ = −0.5, λ = −1, λ = −1.2, and λ = −1.3all approach to , while the one with λ = −1.5 does not.As λ decreases from −0.5 to −1 and further to −1.2, the convergent rate of the HAM solution increases.However, as λ = −1.2decreases even a little, the convergent rate decreases (λ = −1.3).It can even lead to the misconvergence of the HAM series (λ = −1.5).Therefore, on one hand one would expect to choose λ small enough to accelerate the convergence of the HAM series.On the other hand, it is prone to choose an improper one.In this study, λ = −1 is a good choice.

Homotopy-Padé Technique.
In order to achieve faster convergence of HAM series, currently, researchers introduced some optimal approaches and developed the optimal approaches [23,24].Also, the homotopy-Padé technique was proposed to accelerate the convergence of HAM series, [25].
In order to obtain the [m, n] Pade approximation of the HAM series, one should first compute all (m+n)th-order HAM approximations.Therefore, the [m, n] Pade approximations for the frequency and the pitch amplitudes are compared with their corresponding (m+n)th-order HAM solutions, respectively, as shown in Tables 2 and 3. Table 2 shows that the Pade approximations are more accurate than the corresponding HAM solutions, especially when m and n are relatively large.That implies the homotopy-Padé technique can really accelerate the convergence of the HAM series.As Table 3 shows, when U * = 2U * L and λ = −2, the HAM series are disconvergent at p = 1.While the HAM-Pade approximations can still approach to the highly accurate solution.In such a case, the convergent region of the HAM series is enlarged by the homotopy-Padé technique.5.4.About the CPU Time.Next, we will discuss why it is necessary and worthwhile to employ the approach for series expansion of fractional function, as shown in Section 4. The usually adopted procedure for deducing the higher-order deformation equations is differentiating the zeroth-order deformation equations k times, dividing them by k!, and then setting p = 0. Denote the CPU time needed in seeking the nth-order HAM approximations by T n when using the routine procedure, and by S n when employing the means presented in Section 4. Figure 8 shows the ratio between T n and S n versus varying n.When n >10, T n is more than S n by one order of magnitude.Moreover, it increases more and more rapidly as n increases.The presented technique can indeed save a large amount of computational effort.Table 4 shows the comparison of the respective CPU time needed in obtaining the nth-order HAM solution, even seeking the 120th-order solution.

Conclusions
Based on the HAM, we have proposed an approach for obtaining highly accurate approximations for LCOs of strongly nonlinear aeroelastic systems.An easy-to-use approach is proposed to tackle the difficulty in expanding fractional functions into the Taylor series.With the help of this approach, the HAM approximations can be obtained to a very high order and hence can provide solutions to any desired accuracy.The attained HAM solutions are almost the same as the numerical results.Since it is tough to achieve solutions to such high precision, even via the numerical solutions, thus our approaches can be used to validate other solution methods.Also, numerical examples demonstrate that the presented approaches are valid for both weakly and strongly nonlinear aeroelastic systems.These imply that the presented approaches could be applicable in more nonlinear problems, especially those with fractional functions.As mentioned above, the first HAM approximation is in essence the HB1 solution.Note also that the HB1 method is incapable of obtaining the LCO solutions, when U * increases beyond the secondary point, that is, U * = 2.1U * L .Therefore, the presented HAM fails in seeking the solution after the secondary point.In order to do so, one could give the initial solution guess (i.e., (11)) with the third harmonics, so that the initial solution can be determined.In addition, both Figure 8 and Table 4 show that the presented technique can indeed save a large amount of computational effort.The HAM approximations can be obtained to as high as 120th-order within less than half an hour at a microcomputer.As long as the auxiliary parameter is properly chosen, the 100th-order HAM solutions are precise to more than 14 decimals, as implied by Figures 6  and 7, respectively.It is fair to say the presented approach is capable of providing solution to very high precision.
As for the proposed approach for expanding fractional functions, the problem about the robustness of the calculation should be paid special attention in practicing.For example, the coefficient matrix of θ i 's could be illconditioned or singular, which could result in additional numerical error.

Appendix
We have the following Expressions of the Coefficients in (5):

Figure 2 :
Figure 2: Comparisons of the 50th-order HAM solutions for LCO amplitudes with HB1 results and numerical ones.Dots denote the HAM solutions obtained with λ = −1, dashed lines represent HB1 results, and real lines denote numerical solutions.

Figure 3 :
Figure 3: Comparisons of the 50th-order HAM solutions for LCO frequencies with HB1 results and numerical ones.Dots denote the HAM solutions obtained with λ = −1, dashed lines represent HB1 results, and real lines denote numerical solutions.

Figure 4 :Figure 5 :
Figure 4: The phase planes of LCOs of system (1) with U * = 2U * L .Dots denote the 50th-order HAM solution with λ = −1, real lines the numerical result, and dashed lines the HB1 one.

Figure 8 :
Figure 8: The ratio between the CPU times T n and S n , where λ = −1 and U * = 1.5U * L .

Table 1 :
Comparisons of the amplitudes and frequencies obtained by the HAM (λ = −1) with numerical solutions.

Table 2 :
Comparisons of the amplitudes and frequencies given by the HAM Pade approximations (λ = −1) with the numerical solutions, where U * = 1.5U * L .

Table 3 :
Comparisons of the amplitudes and frequencies given by the HAM Pade approximations (λ = −2) with the numerical solutions, where U * = 2U * L .

Table 4 :
The CPU time needed in seeking the nth-order HAM approximations, the parameter values are U * = 1.5U * L and λ = −1.