Table of Contents Author Guidelines Submit a Manuscript
Advances in Acoustics and Vibration
Volume 2011 (2011), Article ID 926271, 10 pages
http://dx.doi.org/10.1155/2011/926271
Research Article

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

1Nari Technology Development Limited Company, 20 High-Tech Road, Nanjing 210061, China
2Department of Mechanics, Sun Yat-sen University, 135 Xingang Road, Guangzhou 510275, China
3State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, China

Received 22 September 2010; Accepted 20 April 2011

Academic Editor: Kok Keong Choong

Copyright © 2011 Y. P. Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 [1416]. 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 [1721].

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 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 _, positive in the downward direction. The elastic axis is located at a distance 𝑎𝑏 from the midchord, while the mass center is located at a distance 𝑥𝛼𝑏 from the elastic axis. Both distances are positive when measured towards the trailing edge of the airfoil.

926271.fig.001
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:..𝜉+𝑥𝛼..𝛼+2𝜁𝜉𝜔𝑈̇𝜉+𝜔𝑈21𝜉=𝐶𝜋𝜇𝐿(𝑡)+𝑃(𝑡)𝑏𝑚𝑈2,𝑥𝛼𝑟2𝛼..𝜉+..𝛼+2𝜁𝛼1𝑈1̇𝛼+𝑈2𝛼+𝜂𝛼3=2𝜋𝜇𝑟2𝛼𝐶𝑀(𝑡)+𝑄(𝑡)𝑚𝑈2𝑟2𝛼,(1) where the superscript denotes the differentiation with respect to the nondimensional time 𝑡 ,   defined as 𝑡=𝑈𝑡1/𝑏, and 𝑡1 is the real time. 𝜉=_/𝑏 is the nondimensional plunge displacement; 𝜂 is the coefficient of cubic pitching stiffness; 𝑈is a nondimensional flow velocity defined as 𝑈=𝑈/(𝑏𝜔𝛼), and 𝜔 is given by 𝜔=𝜔𝜉/𝜔𝛼, where 𝜔𝜉 and 𝜔𝛼 are the natural frequencies of the uncoupled plunging and pitching modes, respectively; 𝜁𝜉 and 𝜁𝛼𝜓  are the damping ratios; 𝑟𝛼 is the radius of gyration about the elastic axis. 𝑃(𝑡) and 𝑄(𝑡) are the externally applied force and moment, 𝑚 is the airfoil mass per unit length and 𝜇𝜓  is the airfoil-air mass ratio. 𝐶𝐿(𝑡) and 𝐶𝑀(𝑡) are the lift and pitching moment coefficients, respectively. For an incompressible flow, the expressions for 𝐶𝐿(𝑡) and 𝐶𝑀(𝑡) are given by 𝐶𝐿(𝑡)=𝜋..𝜉𝑎..̇1𝛼+̇𝛼+2𝜋𝛼(0)+𝜉(0)+2𝑎̇𝛼(0)𝜙(𝜏)+2𝜋𝑡0×𝜙(𝑡𝜎)̇𝛼(𝜎)+..1𝜉(𝜎)+2𝑎..𝐶𝛼(𝜎)𝑑𝜎,(2)𝑀1(𝑡)=𝜋2+𝑎̇1𝛼(0)+𝜉(0)+2𝑎1̇𝛼(0)𝜙(𝜏)+𝜋2+𝑎𝑡0×𝜙(𝑡𝜎)̇𝛼(𝜎)+..1𝜉(𝜎)+2𝑎..+𝜋𝛼(𝜎)𝑑𝜎2𝑎..𝜉𝑎..𝛼12𝑎𝜋2𝜋̇𝛼16..𝛼,(3) where the Wagner function 𝜙(𝑡) is given by Jone’s approximation, 𝜙(𝑡)=1𝜓1𝑒𝜀1𝑡𝜓2𝑒𝜀2𝑡 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. [46] introduced the following four new variables𝑤1=𝑡0𝑒𝜀1(𝑡𝜎)𝛼(𝜎)𝑑𝜎,𝑤2=𝑡0𝑒𝜀2(𝑡𝜎)𝑤𝛼(𝜎)𝑑𝜎,3=𝑡0𝑒𝜀1(𝑡𝜎)𝜉(𝜎)𝑑𝜎,𝑤4=𝑡0𝑒𝜀2(𝑡𝜎)𝜉(𝜎)𝑑𝜎.(4) System (1) can then be rewritten in a general form containing only differential operators as𝑐0..𝜉+𝑐1..𝛼+𝑐2̇𝜉+𝑐3̇𝛼+𝑐4𝜉+𝑐5𝛼+𝑐6𝑤1+𝑐7𝑤2+𝑐8𝑤3+𝑐9𝑤4+𝑐10𝑑𝐺(𝜉)=𝑓(𝑡),0..𝜉+𝑑1..𝛼+𝑑2̇𝜉+𝑑3̇𝛼+𝑑4𝜉+𝑑5𝛼+𝑑6𝑤1+𝑑7𝑤2+𝑑8𝑤3+𝑑9𝑤4+𝑑10𝑀(𝛼)=𝑔(𝑡).(5) The coefficients 𝑐0,𝑐1,,𝑐10;𝑑0,𝑑1,,𝑑10 are given in the appendix, 𝑓(𝑡) and 𝑔(𝑡) are functions depending on initial conditions, Wagner’s function, and the forcing terms. The nonlinear restoring forces, 𝐺(𝜉) and 𝑀(𝛼), are expressed as 𝐺(𝜉)=𝛾𝜉3 and 𝑀(𝛼)=𝜂𝛼3, respectively, with 𝛾 and 𝜂 as coefficients.

By introducing a variable vector 𝐗=(𝑥1,𝑥2,,𝑥8)T, where the superscript “T” denotes the transpose of a matrix, with 𝑥1=𝛼, 𝑥2=̇𝛼, 𝑥3=𝜉, 𝑥4=̇𝜉, 𝑥5=𝑤1, 𝑥6=𝑤2, 𝑥7=𝑤3, and 𝑥8=𝑤4, the coupled equations given in (5) can be written as a set of eight first-order ordinary differential equations written in vector forṁ𝐗=𝐘(𝐗,𝑡).(6) 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 [46].

3. Homotopy Analysis Method

It is assumed that there is no external forces, that is, 𝑄(𝑡)=𝑃(𝑡)=0 in (1). For large values of 𝑡 when transients are damped out and steady solutions are obtained, we can let 𝑓(𝑡)=𝑔(𝑡)=0. Then, (5) can be rewritten in vector form as𝐌..𝐱̇+𝝁𝐱+𝐊𝐱+𝐂𝐖(𝐱)+𝐅(𝐱)=0,(7) where 𝐱=[𝜉,𝛼]T, 𝐖(𝐱)=[W1W2W3W4]T,𝑐𝐌=0𝑐1𝑑0𝑑1𝑐,𝝁=2𝑐3𝑑2𝑑3,𝑐𝐊=4+𝑐10𝑐5𝑑4𝑑5+𝑑10𝑐,𝐂=6𝑐7𝑐8𝑐9𝑑6𝑑7𝑑8𝑑9,(8) and 𝐅(𝐱)=[0𝑑10𝜂𝛼3]T.

Firstly, introduce a new time scale𝜏=𝜔𝑡,(9) where 𝜔 denotes the frequency of the LCO. Then, (7) becomes𝜔2𝐌𝐱+𝜔𝝁𝐱+𝐊𝐱+𝐂𝐖(𝐱,𝜔)+𝐅(𝐱)=0,(10) where the superscript denotes the differentiation with respect to 𝜏. Considering that LCOs are independent of initial conditions, one can adopt the following initial conditions:𝐱(0)=𝑎T,𝐱(0)=𝛽0T.(11) The LCOs of system (10), (11) are periodic motions with frequency 𝜔; thus, 𝐱 can be expressed in a Fourier series𝐱=𝑘=0𝐜𝑘cos𝑘𝜏+𝐬𝑘sin𝑘𝜏,(12) where 𝐜𝑘,𝐬𝑘 are the coefficients in 2×1 vector form.

Let 𝑎0,0,𝜔0,𝛽0, and 𝐱0(𝜏) denote the initial approximations of 𝑎,,𝜔,𝛽, and 𝐱(𝜏), respectively. Due to solution expression (12) and initial conditions (11), the initial guess of solution can be chosen as𝐱0(𝜏)=0cos𝜏+𝛽0sin𝜏𝑎0cos𝜏T.(13)

The homotopy analysis method is based on such continuous variations, 𝐴(𝑝),𝐻(𝑝),Ω(𝑝),𝐵(𝑝), and 𝐮(𝜏,𝑝), that, as the embedding parameter 𝑝 increases from 0 to 1, 𝐮(𝜏,𝑝) varies from the initial guess 𝐱0(𝜏) to the exact solution, so do 𝐴(𝑝),𝐻(𝑝),Ω(𝑝),𝐵(𝑝) from the initial approximations 𝑎0,0,𝜔0,𝛽0 to 𝑎,,𝜔,𝛽, respectively.

Based on (12), one may choose the linear auxiliary operator as𝐿[]𝐮(𝜏,𝑝)=𝜔20𝜕2𝐮(𝜏,𝑝)𝜕𝜏2+𝐮(𝜏,𝑝).(14) Thus𝐿cos𝜏sin𝜏=0.(15)One may define the nonlinear operator according to (10),𝑁[]𝐮(𝜏,𝑝),Ω(𝑝)=Ω2𝜕(𝑝)𝐌2𝐮(𝜏,𝑝)𝜕𝜏2+Ω(𝑝)𝝁𝜕𝐮(𝜏,𝑝)𝜕𝜏+𝐊𝐮(𝜏,𝑝)+𝐂𝐖(𝐮(𝜏,𝑝),Ω(𝑝))+𝐅(𝐮(𝜏,𝑝)).(16) Using the two operators, a family of equations can then be constructed as𝐮(1𝑝)𝐿(𝜏,𝑝)𝐱0[𝐮](𝜏)=𝜆𝑝𝑁(𝜏,𝑝),Ω(𝑝)(17) subject to the initial conditions𝐮(0,𝑝)=𝐻(𝑝)𝐴(𝑝)𝑇,𝜕𝐮(𝜏,𝑝)|||𝜕𝜏𝜏=0=𝐵(𝑝)0𝑇,(18) where the auxiliary parameter 𝜆 is a nonzero constant. Equations (17) and (18) are called the zeroth-order deformation equation.

When 𝑝=0, (17) and (18) have the solution𝐮(𝜏,0)=𝐱0(𝜏).(19) When 𝑝=1, they are exactly the same as (10) and (11) provided thatΩ𝐮(𝜏,1)=𝐱(𝜏),𝐴(1)=𝑎,𝐻(1)=,(1)=𝜔,𝐵(1)=𝛽.(20)

Expand 𝐴(𝑝),𝐻(𝑝),Ω(𝑝),𝐵(𝑝), and 𝐮(𝜏,𝑝) as the series𝐮(𝜏,𝑝)=𝑘=0𝐮𝑘(𝜏)𝑝𝑘,𝐴(𝑝)=𝑘=0𝑎𝑘𝑝𝑘,𝐻(𝑝)=𝑘=0𝑘𝑝𝑘,Ω(𝑝)=𝑘=0𝜔𝑘𝑝𝑘,𝐵(𝑝)=𝑘=0𝛽𝑘𝑝𝑘.(21) As long as the parameter 𝜆 is properly chosen, all of these series are convergent at 𝑝=1. Then, the nth-order HAM solutions can be given as𝐱(𝜏)=𝑛𝑘=0𝐮𝑘(𝜏),𝑎=𝑛𝑘=0𝑎𝑘,=𝑛𝑘=0𝑘,𝜔=𝑛𝑘=0𝜔𝑘,𝛽=𝑛𝑘=0𝛽𝑘.(22)

Substituting (21) into (17) and (18), differentiating (17) and (18) k times, dividing the differentiations by k! and then letting 𝑝=0, one can obtain the kth-order (𝑘1) deformation equation𝐿𝐮𝑘+1(𝜏)𝜒𝑘1𝐮𝑘(𝜏)=𝜆𝐑𝑘(𝜏)(23) subject to the initial conditions𝐮𝑘(0)=𝑘𝑎𝑘𝑇,𝐮𝑘𝛽(0)=𝑘0𝑇,(24) where𝐑𝑘1(𝜏)=𝜕(𝑘1)!𝑘1𝑁[]𝐮(𝜏,𝑝),Ω(𝑝)𝜕𝑝𝑘1||||𝑝=0,(25)𝜒𝑘=0,𝑘=1,1,𝑘2.(26)

Due to the rule of solution expression and the linear operator 𝐿, 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𝚪𝑐𝑘𝑎𝑘,𝑘,𝜔𝑘,𝛽𝑘=1𝜋02𝜋𝐑𝑘𝚪(𝜏)cos𝜏𝑑𝜏=0,𝑠𝑘𝑎𝑘,𝑘,𝜔𝑘,𝛽𝑘=1𝜋02𝜋𝐑𝑘(𝜏)sin𝜏𝑑𝜏=0.(27) Solving (27), 𝑎𝑘,𝑘,𝜔𝑘, and 𝛽𝑘 are determined step by step as 𝑘 increases.

Note that when 𝑘=0, 𝐑𝑘+1(𝜏) is essentially the right hand side of (10) with 𝐱=𝐱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 𝑎0,0,𝜔0, and 𝛽0 can be obtained by using the Newton-Raphson method. Importantly, (27) is always linear as 𝑘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 obtain 𝑅𝑘(𝜏) 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) 𝑘 times, dividing them by 𝑘!, and then setting 𝑝=0. This kind of approach is based on the classical theories of the Taylor series. In our study, however, using this method to expand 𝐂𝐖(𝐮(𝜏,𝑝),Ω(𝑝)) will cost a large amount of computational resources. For example, substitution of 𝛼=cos𝜏 into 𝑤1=𝑡0𝑒𝜀1(𝑡𝜎)𝛼(𝜎)𝑑𝜎 yields a simple illustration𝑡0𝑒𝜀1(𝑡𝜎)𝜀cos(𝑖𝜔𝜎)𝑑𝜎=1cos(𝑖𝜔𝑡)+𝑖𝜔sin(𝑖𝜔𝑡)𝜀21+𝑖2𝜔2𝜀1𝜀21+𝑖2𝜔2𝑒𝜀1𝑡.(28) For large values of 𝑡, the second term in (28) approaches to zero and is neglected since only steady solutions (LCOs) are taken into account. Therefore, the integrations 𝐖(cos(𝑖𝜔𝑡)) and 𝐖(sin(𝑖𝜔𝑡)) can be expressed as follows, respectively:𝐖1(cos(𝑖𝜔𝑡))=𝐖3=𝜀(cos(𝑖𝜔𝑡))1cos(𝑖𝜔𝑡)+𝑖𝜔sin(𝑖𝜔𝑡)𝜀21+𝑖2𝜔2,𝐖2(cos(𝑖𝜔𝑡))=𝐖4=𝜀(cos(𝑖𝜔𝑡))2cos(𝑖𝜔𝑡)+𝑖𝜔sin(𝑖𝜔𝑡)𝜀22+𝑖2𝜔2,𝐖1(sin(𝑖𝜔𝑡))=𝐖3=(sin(𝑖𝜔𝑡))𝑖𝜔cos(𝑖𝜔𝑡)+𝜀1sin(𝑖𝜔𝑡)𝜀21+𝑖2𝜔2,𝐖2(sin(𝑖𝜔𝑡))=𝐖4=(sin(𝑖𝜔𝑡))𝑖𝜔cos(𝑖𝜔𝑡)+𝜀2sin(𝑖𝜔𝑡)𝜀22+𝑖2𝜔2,(29) where 𝐖𝑗,𝑗=1,2,3,4 correspond to the jth component of vector 𝐖, and 𝑖=1,2,, Differentiating 𝑘 times 1/[𝜀21+(𝑛𝑖=0𝜔𝑖𝑝𝑖)2] with respect to 𝑝 will result in a complicated term (i.e., [𝜀21+(𝑛𝑖=0𝜔𝑖𝑝𝑖)2]𝑘) in the denominator. Thereby, the expression of  𝜕𝑘𝐖/𝜕𝑝𝑘  becomes more and more complex as 𝑘 increases, which makes it quite tough to deduce high-order deformation equations.

We take 1/[𝜀21+(𝑛𝑖=0𝜔𝑖𝑝𝑖)2] as an illustrative example to propose a means for expanding fractional functions. First of all, denote the denominator as 𝜀21+𝑛𝑖=0𝜔𝑖𝑝𝑖2=𝜀21+𝜔20+𝑛𝑘=1𝑘𝑖=0𝜔𝑖𝜔𝑘𝑖𝑝𝑘=𝑛𝑘=0𝜎𝑘𝑝𝑘,(30) where 𝜎0=𝜀21+𝜔20 and 𝜎𝑘=𝑘𝑖=0𝜔𝑖𝜔𝑘𝑖,1𝑘𝑛. Taking the nth-order Taylor series of 1/(𝑛𝑘=0𝜎𝑘𝑝𝑘) as 𝑛𝑘=0𝜃𝑘𝑝𝑘, then one has1𝑛𝑘=0𝜎𝑘𝑝𝑘=𝑛𝑘=0𝜃𝑘𝑝𝑘𝑝+𝑂𝑘+1.(31) Rewrite (31) as𝑛𝑘=0𝜎𝑘𝑝𝑘𝑛𝑘=0𝜃𝑘𝑝𝑘𝑝+𝑂𝑛+1=𝑛𝑘=0𝑘𝑖=0𝜎𝑖𝜃𝑘𝑖𝑝𝑘𝑝+𝑂𝑘+1=1.(32) Equating the coefficients of 𝑝𝑘 results in𝜎0𝜃0=1,𝑘𝑖=0𝜎𝑖𝜃𝑘𝑖=0,𝑘=1,2,,𝑛.(33) Interestingly, (33) is always linear. That means it is rather easy to determine 𝜃𝑘 if 𝜎𝑖 are all known, 𝑖=0,1,2,,𝑘.

5. Numerical Examples

5.1. Main Results

The system parameters under consideration are 𝜇=100, 𝑟𝛼=0.5, 𝑎=0.5, 𝜁𝛼=𝜁𝜉=0, 𝜔=0.25, 𝑥𝛼=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 at 𝑈=𝑈𝐿=6.0385 [4, 5]. As 𝑈 increases beyond 𝑈𝐿, LCO arises, and thus 𝑈𝐿 is a Hopf bifurcation point. Note that the flutter boundary 𝑈𝐿 is independent of the nonlinear coefficient 𝜂. Lee et al. [4] found a secondary Hopf bifurcation as 𝑈 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.

fig2
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.
926271.fig.003
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 when 𝑈 is in [𝑈𝐿,2𝑈𝐿] 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 𝑈.

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 𝑈 is larger than the secondary bifurcation value, about 2𝑈𝐿, 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.

fig4
Figure 4: The phase planes of LCOs of system (1) with 𝑈=2𝑈𝐿. Dots denote the 50th-order HAM solution with 𝜆=1, real lines the numerical result, and dashed lines the HB1 one.
fig5
Figure 5: The time history responses of system (1) with 𝑈=1.5𝑈𝐿. Dots denote the 50th-order HAM solutions with 𝜆=1 and real lines the numerical results.

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 𝑛 >120, 𝑎𝑛,𝑛,𝜔𝑛, and 𝛽𝑛 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.

tab1
Table 1: Comparisons of the amplitudes and frequencies obtained by the HAM (𝜆=1) with numerical solutions.
fig6
Figure 6: Residues of (6) (𝑈=1.5𝑈𝐿) with HAM solutions attained with 𝜆=1, (a): 𝑁=40; (b): 𝑁=80; (c): 𝑁=120, where the real and dash lines denote the residues of the first and the second equations, respectively.

Very interestingly, it is found that 𝜔𝑖 is independent of 𝜂, while 𝑎𝑖,𝛽𝑖, and 𝑖 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 𝑈=1.5𝑈𝐿 as 𝜛= 0.07756360647090. Denote the discrepancy between the nth-order HAM frequency solution and 𝜛 as 𝑒(𝑛)=(𝑛𝑘=0𝜔𝑘)𝜛. Figure 7 shows the 𝜆-curves with respect to 𝜔. Considering that the longitudinal coordinate refers to the logarithm of |𝑒(𝑛)|, 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.

926271.fig.007
Figure 7: The 𝜆-curves with respect to 𝜔, where 𝑈=1.5𝑈𝐿.
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 𝑈=2𝑈𝐿 and 𝜆=2, the HAM series are disconvergent at 𝑝=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.

tab2
Table 2: Comparisons of the amplitudes and frequencies given by the HAM Pade approximations (𝜆=1) with the numerical solutions, where 𝑈=1.5𝑈𝐿.
tab3
Table 3: Comparisons of the amplitudes and frequencies given by the HAM Pade approximations (𝜆=2) with the numerical solutions, where 𝑈=2𝑈𝐿.
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 𝑘 times, dividing them by 𝑘!, and then setting 𝑝=0. Denote the CPU time needed in seeking the nth-order HAM approximations by 𝑇𝑛 when using the routine procedure, and by 𝑆𝑛 when employing the means presented in Section 4. Figure 8 shows the ratio between 𝑇𝑛 and 𝑆𝑛 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.

tab4
Table 4: The CPU time needed in seeking the nth-order HAM approximations, the parameter values are 𝑈=1.5𝑈𝐿 and 𝜆=1.
926271.fig.008
Figure 8: The ratio between the CPU times Tn and Sn, where 𝜆=1 and 𝑈=1.5𝑈𝐿.

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, when 𝑈 increases beyond the secondary point, that is, 𝑈=2.1𝑈𝐿. 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 𝜃𝑖’s could be illconditioned or singular, which could result in additional numerical error.

Appendix

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

𝑐01=1+𝜇,𝑐1=𝑥𝛼𝑎𝜇,𝑐2=2𝜇1𝜓1𝜓2+2𝜁𝜉𝜔𝑈,𝑐3=1𝜇1+12𝑎1𝜓1𝜓2,𝑐4=2𝜇𝜀1𝜓1+𝜀2𝜓2,𝑐5=2𝜇1𝜓1𝜓2+12𝑎𝜀1𝜓1+𝜀2𝜓2,𝑐6=2𝜇𝜀1𝜓11𝜀112𝑎,𝑐7=2𝜇𝜀2𝜓21𝜀212𝑎,𝑐82=𝜇𝜀21𝜓1,𝑐92=𝜇𝜀22𝜓2,𝑐10=𝜔𝑈2,𝑑0=𝑥𝛼𝑟2𝛼𝑎𝜇𝑟2𝛼,𝑑1=1+1+8𝑎28𝜇𝑟2𝛼,𝑑2=1+2𝑎𝜇𝑟2𝛼𝜀1𝜓1+𝜀2𝜓2,𝑑3=12𝑎2𝜇𝑟2𝛼14𝑎21𝜓1𝜓22𝜇𝑟2𝛼+2𝜁𝛼𝑈,𝑑4=1+2𝑎𝜇𝑟2𝛼𝜀1𝜓1+𝜀2𝜓2,𝑑5=1+2𝑎𝜇𝑟2𝛼1𝜓1𝜓21+2𝑎12𝑎𝜓1𝜀1𝜓2𝜀22𝜇𝑟2𝛼,𝑑6=1+2𝑎𝜓1𝜀1𝜇𝑟2𝛼1𝜀112𝑎,𝑑7=1+2𝑎𝜓2𝜀2𝜇𝑟2𝛼1𝜀212𝑎,𝑑8=1+2𝑎𝜓1𝜀21𝜇𝑟2𝛼,𝑑9=1+2𝑎𝜓2𝜀22𝜇𝑟2𝛼,𝑑10=1𝑈2.(A.1)

Acknowledgments

This work is supported by the National Natural Science Foundation of China (10772202, 10102023, 10972241, 11032005), Doctoral Program Foundation of Ministry of Education of China (20090171110042), and Educational Commission of Guangdong Province of China (34310018).

References

  1. B. H. K. Lee, S. J. Price, and Y. S. Wong, “Nonlinear aeroelastic analysis of airfoils: bifurcation and chaos,” Progress in Aerospace Sciences, vol. 35, no. 3, pp. 205–334, 1999. View at Publisher · View at Google Scholar · View at Scopus
  2. J. K. Liu and L. C. Zhao, “Bifurcation analysis of airfoils in incompressible flow,” Journal of Sound and Vibration, vol. 154, no. 1, pp. 117–124, 1991. View at Google Scholar · View at Scopus
  3. Y. M. Chen and J. K. Liu, “On the limit cycles of aeroelastic systems with quadratic nonlinearities,” Structural Engineering and Mechanics, vol. 30, no. 1, pp. 67–76, 2008. View at Google Scholar · View at Scopus
  4. B. H. K. Lee, L. Liu, and K. W. Chung, “Airfoil motion in subsonic flow with strong cubic nonlinear restoring forces,” Journal of Sound and Vibration, vol. 281, no. 3–5, pp. 699–717, 2005. View at Publisher · View at Google Scholar · View at Scopus
  5. L. P. Liu and E. H. Dowell, “The secondary bifurcation of an aeroelastic airfoil motion: effect of high harmonics,” Nonlinear Dynamics, vol. 37, no. 1, pp. 31–49, 2004. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  6. L. Liu, E. H. Dowell, and J. P. Thomas, “A high dimensional harmonic balance approach for an aeroelastic airfoil with cubic restoring forces,” Journal of Fluids and Structures, vol. 23, no. 3, pp. 351–363, 2005. View at Publisher · View at Google Scholar · View at Scopus
  7. P. Shahrzad and M. Mahzoon, “Limit cycle flutter of airfoils in steady and unsteady flows,” Journal of Sound and Vibration, vol. 256, no. 2, pp. 213–225, 2002. View at Publisher · View at Google Scholar · View at Scopus
  8. M. Cai, J. K. Liu, and J. Li, “Incremental harmonic balance method for airfoil flutter with multiple strong nonlinearities,” Applied Mathematics and Mechanics (English Edition), vol. 27, no. 7, pp. 953–958, 2006. View at Publisher · View at Google Scholar · View at Scopus
  9. K. W. Chung, C. L. Chan, and B. H. K. Lee, “Bifurcation analysis of a two-degree-of-freedom aeroelastic system with freeplay structural nonlinearity by a perturbation-incremental method,” Journal of Sound and Vibration, vol. 299, no. 3, pp. 520–539, 2007. View at Publisher · View at Google Scholar · View at Scopus
  10. K. W. Chung, Y. B. He, and B. H. K. Lee, “Bifurcation analysis of a two-degree-of-freedom aeroelastic system with hysteresis structural nonlinearity by a perturbation-incremental method,” Journal of Sound and Vibration, vol. 320, no. 1-2, pp. 163–183, 2009. View at Publisher · View at Google Scholar · View at Scopus
  11. L. Liu, Y. S. Wong, and B. H. K. Lee, “Application of the centre manifold theory in non-linear aeroelasticity,” Journal of Sound and Vibration, vol. 234, no. 4, pp. 641–659, 2000. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  12. Q. Ding and D. L. Wang, “The flutter of an airfoil with cubic structural and aerodynamic non-linearities,” Aerospace Science and Technology, vol. 10, no. 5, pp. 427–434, 2006. View at Publisher · View at Google Scholar · View at Scopus
  13. J. Grzedziński, “Limitation of application of the center manifold reduction in aeroelasticity,” Journal of Fluids and Structures, vol. 21, no. 2, pp. 187–209, 2005. View at Publisher · View at Google Scholar · View at Scopus
  14. S. J. Liao, “A kind of approximate solution technique which does not depend upon small parameters-II: an application in fluid mechanics,” International Journal of Non-Linear Mechanics, vol. 32, no. 5, pp. 815–822, 1997. View at Google Scholar · View at Scopus
  15. S. J. Liao, “On the homotopy analysis method for nonlinear problems,” Applied Mathematics and Computation, vol. 147, no. 2, pp. 499–513, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  16. S. J. Liao, Beyond Perturbation: Introduction to Homotopy Analysis Method, Chapman & Hall/CRC Press, Boca Raton, Fla, USA, 2003.
  17. T. Pirbodaghi, M. T. Ahmadian, and M. Fesanghary, “On the homotopy analysis method for non-linear vibration of beams,” Mechanics Research Communications, vol. 36, no. 2, pp. 143–148, 2009. View at Publisher · View at Google Scholar · View at Scopus
  18. Y. M. Chen and J. K. Liu, “A study of homotopy analysis method for limit cycle of van der Pol equation,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 5, pp. 1816–1821, 2009. View at Publisher · View at Google Scholar · View at Scopus
  19. B. T. Kennedy, D. C. Weggel, D. M. Boyajian, and R. E. Smelser, “Closed-form solution for a cantilevered sectorial plate subjected to a twisting tip moment,” Mechanics Research Communications, vol. 35, no. 8, pp. 491–496, 2008. View at Publisher · View at Google Scholar · View at Scopus
  20. Y. M. Chen and J. K. Liu, “Uniformly valid solution of limit cycle of the Duffing-van der Pol equation,” Mechanics Research Communications, vol. 36, no. 7, pp. 845–850, 2009. View at Publisher · View at Google Scholar · View at Scopus
  21. Y. M. Chen and J. K. Liu, “Homotopy analysis method for limit cycle flutter of airfoils,” Applied Mathematics and Computation, vol. 203, no. 2, pp. 854–863, 2008. View at Publisher · View at Google Scholar · View at Scopus
  22. Y. M. Chen and J. K. Liu, “Homotopy analysis method for limit cycle oscillations of an airfoil with cubic nonlinearities,” JVC/Journal of Vibration and Control, vol. 16, no. 2, pp. 163–179, 2010. View at Publisher · View at Google Scholar · View at Scopus
  23. S. J. Liao, “An optimal homotopy-analysis approach for strongly nonlinear differential equations,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 8, pp. 2003–2016, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  24. Z. Niu and C. Wang, “A one-step optimal homotopy analysis method for nonlinear differential equations,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 8, pp. 2026–2036, 2010. View at Publisher · View at Google Scholar · View at Scopus
  25. Y. H. Qian and S. M. Chen, “Accurate approximate analytical solutions for multi-degree-of-freedom coupled van der Pol-Duffing oscillators by homotopy analysis method,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 10, pp. 3113–3130, 2010. View at Publisher · View at Google Scholar · View at Scopus