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

Authors: Y. P. Zhang; Y. M. Chen; J. K. Liu; G. Meng
Journal: Advances in Acoustics and Vibration
 (2011)

---

## Abstract

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.


---

## Body

## 1. 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 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–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–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.


## 2. Equations of Motions

The physical model shown in Figure1 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 ahb 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.Figure 1
Sketch of a two-dimensional airfoil.For cubic restoring forces with subsonic aerodynamics, the coupled equations for the airfoil in nondimensional form can be written as follows:(1)ξ..+xαα..+2ζξω̅U*ξ.+[ω̅U*]2ξ=-1πμCL(t)+P(t)bmU2,xαrα2ξ..+α..+2ζα1U*α.+[1U*]2(α+ηα3)=2πμrα2CM(t)+Q(t)mU2rα2,
where the superscript denotes the differentiation with respect to the nondimensional time t ,   defined as t=Ut1/b, and t1 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. CL(t) and CM(t) are the lift and pitching moment coefficients, respectively. For an incompressible flow, the expressions for CL(t) and CM(t) are given by
(2)CL(t)=π(ξ..-ahα..+α.)+2π[α(0)+ξ.(0)+(12-ah)α.(0)]ϕ(τ)+2π∫0tϕ(t-σ)×[α.(σ)+ξ..(σ)+(12-ah)α..(σ)]dσ,(3)CM(t)=π(12+ah)[α(0)+ξ.(0)+(12-ah)α.(0)]ϕ(τ)+π(12+ah)∫0tϕ(t-σ)×[α.(σ)+ξ..(σ)+(12-ah)α..(σ)]dσ+π2ah(ξ..-ahα..)-(12-ah)π2α.-π16α..,
where the Wagner function ϕ(t) is given by Jone’s approximation, ϕ(t)=1-ψ1e-ε1t-ψ2e-ε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–6] introduced the following four new variables(4)w1=∫0te-ε1(t-σ)α(σ)dσ,w2=∫0te-ε2(t-σ)α(σ)dσ,w3=∫0te-ε1(t-σ)ξ(σ)dσ,w4=∫0te-ε2(t-σ)ξ(σ)dσ.
System (1) can then be rewritten in a general form containing only differential operators as(5)c0ξ..+c1α..+c2ξ.+c3α.+c4ξ+c5α+c6w1+c7w2+c8w3+c9w4+c10G(ξ)=f(t),d0ξ..+d1α..+d2ξ.+d3α.+d4ξ+d5α+d6w1+d7w2+d8w3+d9w4+d10M(α)=g(t).
The coefficients c0,c1,…,c10;d0,d1,…,d10 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 vectorX=(x1,x2,…,x8)T, where the superscript “T” denotes the transpose of a matrix, with x1=α, x2=α., x3=ξ, x4=ξ., x5=w1, x6=w2, x7=w3, and x8=w4, the coupled equations given in (5) can be written as a set of eight first-order ordinary differential equations written in vector form(6)X.=Y(X,t).
This approach allows existing methods suitable for the study of ordinary differential equations to be used in the analysis. For more details of (5) and (6), please refer to [4–6].


## 3. 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(7)Mx..+μx.+Kx+CW(x)+F(x)=0,
where x=[ξ,α]T, W(x)=[W1W2W3W4]T,(8)M=[c0c1d0d1],μ=[c2c3d2d3],K=[c4+c10c5d4d5+d10],C=[c6c7c8c9d6d7d8d9],
and F(x)=[0d10ηα3]T.Firstly, introduce a new time scale(9)τ=ωt,
where ω denotes the frequency of the LCO. Then, (7) becomes(10)ω2Mx′′+ωμx′+Kx+CW(x,ω)+F(x)=0,
where the superscript denotes the differentiation with respect to τ. Considering that LCOs are independent of initial conditions, one can adopt the following initial conditions:(11)x(0)=[ha]T,x′(0)=[β0]T.
The LCOs of system (10), (11) are periodic motions with frequency ω; thus, x can be expressed in a Fourier series(12)x=∑k=0∞(ckcos⁡kτ+sksinkτ),
where ck,sk are the coefficients in 2×1 vector form.Leta0,h0,ω0,β0, and x0(τ) 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(13)x0(τ)=[h0cos⁡τ+β0sinτa0cos⁡τ]T.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 x0(τ) to the exact solution, so do A(p),H(p),Ω(p),B(p) from the initial approximations a0,h0,ω0,β0 to a,h,ω,β, respectively.Based on (12), one may choose the linear auxiliary operator as(14)L[u(τ,p)]=ω02[∂2u(τ,p)∂τ2+u(τ,p)].
Thus(15)L[(cos⁡τsinτ)]=0.One may define the nonlinear operator according to (10),(16)N[u(τ,p),Ω(p)]=Ω2(p)M∂2u(τ,p)∂τ2+Ω(p)μ∂u(τ,p)∂τ+Ku(τ,p)+CW(u(τ,p),Ω(p))+F(u(τ,p)).
Using the two operators, a family of equations can then be constructed as(17)(1-p)L[u(τ,p)-x0(τ)]=λpN[u(τ,p),Ω(p)]
subject to the initial conditions(18)u(0,p)=[H(p)A(p)]T,∂u(τ,p)∂τ|τ=0=[B(p)0]T,
where the auxiliary parameter λ is a nonzero constant. Equations (17) and (18) are called the zeroth-order deformation equation.Whenp=0, (17) and (18) have the solution(19)u(τ,0)=x0(τ).
When p=1, they are exactly the same as (10) and (11) provided that(20)u(τ,1)=x(τ),A(1)=a,H(1)=h,Ω(1)=ω,B(1)=β.ExpandA(p),H(p),Ω(p),B(p), and u(τ,p) as the series(21)u(τ,p)=∑k=0∞uk(τ)pk,A(p)=∑k=0∞akpk,H(p)=∑k=0∞hkpk,Ω(p)=∑k=0∞ωkpk,B(p)=∑k=0∞βkpk.
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(22)x(τ)=∑k=0nuk(τ),a=∑k=0nak,h=∑k=0nhk,ω=∑k=0nωk,β=∑k=0nβk.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(23)L[uk+1(τ)-χk-1uk(τ)]=λRk(τ)
subject to the initial conditions(24)uk(0)=[hkak]T,uk′(0)=[βk0]T,
where(25)Rk(τ)=1(k-1)!∂k-1N[u(τ,p),Ω(p)]∂pk-1|p=0,(26)χk={0,k=1,1,k≥2.Due to the rule of solution expression and the linear operatorL, 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(27)Γkc(ak,hk,ωk,βk)=1π∫02πRk(τ)cos⁡τdτ=0,Γks(ak,hk,ωk,βk)=1π∫02πRk(τ)sinτdτ=0.
Solving (27), ak,hk,ωk, and βk are determined step by step as k increases.Note that whenk=0, Rk+1(τ) is essentially the right hand side of (10) with x=x0, 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 a0,h0,ω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 high-order approximations [22].


## 4. Expansion of Fractional Functions

A key procedure of implementing the HAM is to deduce the high-order deformation equation, that is, to obtainRk(τ) 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 α=cosτ into w1=∫0te-ε1(t-σ)α(σ)dσ yields a simple illustration(28)∫0te-ε1(t-σ)cos⁡(iωσ)dσ=ε1cos⁡(iωt)+iωsin(iωt)ε12+i2ω2-ε1ε12+i2ω2e-ε1t.
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:(29)W1(cos⁡(iωt))=W3(cos⁡(iωt))=ε1cos⁡(iωt)+iωsin(iωt)ε12+i2ω2,W2(cos⁡(iωt))=W4(cos⁡(iωt))=ε2cos⁡(iωt)+iωsin(iωt)ε22+i2ω2,W1(sin(iωt))=W3(sin(iωt))=-iωcos⁡(iωt)+ε1sin(iωt)ε12+i2ω2,W2(sin(iωt))=W4(sin(iωt))=-iωcos⁡(iωt)+ε2sin(iωt)ε22+i2ω2,
where Wj,j=1,2,3,4 correspond to the jth component of vector W, and i=1,2,…, Differentiating k times 1/[ε12+(∑i=0nωipi)2] with respect to p will result in a complicated term (i.e., [ε12+(∑i=0nωipi)2]k) in the denominator. Thereby, the expression of  ∂kW/∂pk  becomes more and more complex as k increases, which makes it quite tough to deduce high-order deformation equations.We take1/[ε12+(∑i=0nωipi)2] as an illustrative example to propose a means for expanding fractional functions. First of all, denote the denominator as
(30)ε12+(∑i=0nωipi)2=ε12+ω02+∑k=1n(∑i=0kωiωk-i)pk:=∑k=0nσkpk,
where σ0=ε12+ω02 and σk=∑i=0kωiωk-i,1≤k≤n. Taking the nth-order Taylor series of 1/(∑k=0nσkpk) as ∑k=0nθkpk, then one has(31)1∑k=0nσkpk=∑k=0nθkpk+O(pk+1).
Rewrite (31) as(32)(∑k=0nσkpk)[∑k=0nθkpk+O(pn+1)]=∑k=0n(∑i=0kσiθk-i)pk+O(pk+1)=1.
Equating the coefficients of pk results in(33)σ0θ0=1,∑i=0kσiθk-i=0,k=1,2,…,n.
Interestingly, (33) is always linear. That means it is rather easy to determine θk if σi are all known, i=0,1,2,…,k.


## 5. Numerical Examples

### 5.1. Main Results

The system parameters under consideration areμ=100, rα=0.5, ah=-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)=ξ(0)=ξ.(0)=0.Using analytical techniques developed for nonlinear dynamical systems, the linear flutter speed is found atU*=UL*=6.0385 [4, 5]. As U* increases beyond UL*, LCO arises, and thus UL* is a Hopf bifurcation point. Note that the flutter boundary UL* 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.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.
(a)(b)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.In the proposed method, the zeroth-order HAM approximation is essentially given by the HB1 method. The higher-order 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 whenU* is in [UL*,2UL*] 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 whenU* is larger than the secondary bifurcation value, about 2UL*, so is the presented approach, as shown in Figures 2 and 3.Figures4 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.The phase planes of LCOs of system (1) with U*=2UL*. Dots denote the 50th-order HAM solution with λ=-1, real lines the numerical result, and dashed lines the HB1 one.
(a)(b)The time history responses of system (1) with U*=1.5UL*. Dots denote the 50th-order HAM solutions with λ=-1 and real lines the numerical results.
(a)(b)More precisely, the 120th-order HAM solutions shown in Table1 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, an,hn,ω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.Table 1
Comparisons of the amplitudes and frequencies obtained by the HAM (λ=-1) with numerical solutions.
U*U*=1.5UL*U*=2UL*HAMFrequencyPitchPlungeFrequencyPitchPlungeN=400.0775636561280.137381457860.3568580.06586090.21848360.6945022N=800.0775636064760.137381511720.356858140.06578670.21856460.6964279N=1200.077563606470900.137381511730.356858150.06578330.21856850.6965209Numerical0.07756350.13738160.35685900.06578290.21856890.6965298Residues of (6) (U*=1.5UL*) with HAM solutions attained with λ=-1, (a): N=40; (b): N=80; (c): N=120, where the real and dash lines denote the residues of the first and the second equations, respectively.
(a)(b)(c)Very interestingly, it is found thatωi is independent of η, while ai,βi, and hi 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 η.


### 5.2. 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.5UL* as ϖ= 0.07756360647090. Denote the discrepancy between the nth-order HAM frequency solution and ϖ as e(n)=(∑k=0nω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.3 all 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.2 decreases 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.Figure 7
Theλ-curves with respect to ω, where U*=1.5UL*.


### 5.3. 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*=2UL* 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.Table 2
Comparisons of the amplitudes and frequencies given by the HAM Pade approximations (λ=-1) with the numerical solutions, where U*=1.5UL*.
[m, n] HAM PadeFrequencyRelative error (%)PitchRelative error (%)[5,5]0.077651009776630.110.137450510522940.05HAM 100.077703846285050.180.13726535168074−0.08[7,8]0.077567463573885e−30.13737312674052−6e−3150.077593885939524e−20.13735291378024−2e−2[10,10]0.077563945007644e−40.13738142730826−6e−5HAM 200.077571127404821e−20.13737399877222−5e−3HAMN= 1200.077563606470900.13738151173287Table 3
Comparisons of the amplitudes and frequencies given by the HAM Pade approximations (λ=-2) with the numerical solutions, where U*=2UL*.
[m, n] HAM PadeFrequencyPitch[5,5]0.0670111960.221689949HAM 100.07303670.156371061[8,8]0.065955600.21807163160.52075063.48085267[10,10]0.06587740.21854332HAM 207.63500−28.477895HAMN= 120,0.06578290.2185689


### 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 Section4. 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 Tn when using the routine procedure, and by Sn when employing the means presented in Section 4. Figure 8 shows the ratio between Tn and Sn versus varying n. When n>10, Tn is more than Sn 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.Table 4
The CPU time needed in seeking thenth-order HAM approximations, the parameter values are U*=1.5UL* and λ=-1.
Thenth HAM solutionThe CPU time needed (second)nThe routine procedureThe presented approach10142.815584.6201966.850/4680/239120/1393Figure 8
The ratio between the CPU timesTn and Sn, where λ=-1 and U*=1.5UL*.


## 5.1. Main Results

The system parameters under consideration areμ=100, rα=0.5, ah=-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)=ξ(0)=ξ.(0)=0.Using analytical techniques developed for nonlinear dynamical systems, the linear flutter speed is found atU*=UL*=6.0385 [4, 5]. As U* increases beyond UL*, LCO arises, and thus UL* is a Hopf bifurcation point. Note that the flutter boundary UL* 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.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.
(a)(b)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.In the proposed method, the zeroth-order HAM approximation is essentially given by the HB1 method. The higher-order 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 whenU* is in [UL*,2UL*] 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 whenU* is larger than the secondary bifurcation value, about 2UL*, so is the presented approach, as shown in Figures 2 and 3.Figures4 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.The phase planes of LCOs of system (1) with U*=2UL*. Dots denote the 50th-order HAM solution with λ=-1, real lines the numerical result, and dashed lines the HB1 one.
(a)(b)The time history responses of system (1) with U*=1.5UL*. Dots denote the 50th-order HAM solutions with λ=-1 and real lines the numerical results.
(a)(b)More precisely, the 120th-order HAM solutions shown in Table1 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, an,hn,ω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.Table 1
Comparisons of the amplitudes and frequencies obtained by the HAM (λ=-1) with numerical solutions.
U*U*=1.5UL*U*=2UL*HAMFrequencyPitchPlungeFrequencyPitchPlungeN=400.0775636561280.137381457860.3568580.06586090.21848360.6945022N=800.0775636064760.137381511720.356858140.06578670.21856460.6964279N=1200.077563606470900.137381511730.356858150.06578330.21856850.6965209Numerical0.07756350.13738160.35685900.06578290.21856890.6965298Residues of (6) (U*=1.5UL*) with HAM solutions attained with λ=-1, (a): N=40; (b): N=80; (c): N=120, where the real and dash lines denote the residues of the first and the second equations, respectively.
(a)(b)(c)Very interestingly, it is found thatωi is independent of η, while ai,βi, and hi 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 η.


## 5.2. 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.5UL* as ϖ= 0.07756360647090. Denote the discrepancy between the nth-order HAM frequency solution and ϖ as e(n)=(∑k=0nω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.3 all 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.2 decreases 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.Figure 7
Theλ-curves with respect to ω, where U*=1.5UL*.


## 5.3. 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*=2UL* 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.Table 2
Comparisons of the amplitudes and frequencies given by the HAM Pade approximations (λ=-1) with the numerical solutions, where U*=1.5UL*.
[m, n] HAM PadeFrequencyRelative error (%)PitchRelative error (%)[5,5]0.077651009776630.110.137450510522940.05HAM 100.077703846285050.180.13726535168074−0.08[7,8]0.077567463573885e−30.13737312674052−6e−3150.077593885939524e−20.13735291378024−2e−2[10,10]0.077563945007644e−40.13738142730826−6e−5HAM 200.077571127404821e−20.13737399877222−5e−3HAMN= 1200.077563606470900.13738151173287Table 3
Comparisons of the amplitudes and frequencies given by the HAM Pade approximations (λ=-2) with the numerical solutions, where U*=2UL*.
[m, n] HAM PadeFrequencyPitch[5,5]0.0670111960.221689949HAM 100.07303670.156371061[8,8]0.065955600.21807163160.52075063.48085267[10,10]0.06587740.21854332HAM 207.63500−28.477895HAMN= 120,0.06578290.2185689


## 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 Section4. 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 Tn when using the routine procedure, and by Sn when employing the means presented in Section 4. Figure 8 shows the ratio between Tn and Sn versus varying n. When n>10, Tn is more than Sn 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.Table 4
The CPU time needed in seeking thenth-order HAM approximations, the parameter values are U*=1.5UL* and λ=-1.
Thenth HAM solutionThe CPU time needed (second)nThe routine procedureThe presented approach10142.815584.6201966.850/4680/239120/1393Figure 8
The ratio between the CPU timesTn and Sn, where λ=-1 and U*=1.5UL*.


## 6. 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, whenU* increases beyond the secondary point, that is, U*=2.1UL*. 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 Figure8 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.



---

*Source: 926271-2011-06-23.xml*