`Advances in Acoustics and VibrationVolume 2011, Article ID 926271, 10 pageshttp://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

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.

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: where the superscript denotes the differentiation with respect to the nondimensional time  ,   defined as , and 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 where the Wagner function is given by Jone’s approximation, with the constants as , , , and .

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 System (1) can then be rewritten in a general form containing only differential operators as The coefficients 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 and , respectively, with and as coefficients.

By introducing a variable vector , where the superscript “T” denotes the transpose of a matrix, with , , , , , , , and , 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. 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, in (1). For large values of when transients are damped out and steady solutions are obtained, we can let . Then, (5) can be rewritten in vector form as where , and .

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, can be expressed in a Fourier series where are the coefficients in vector form.

Let , and 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

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 to the exact solution, so do from the initial approximations to , respectively.

Based on (12), one may choose the linear auxiliary operator as ThusOne 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.

When , (17) and (18) have the solution When , they are exactly the same as (10) and (11) provided that

Expand , and as the series As long as the parameter is properly chosen, all of these series are convergent at . 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 , one can obtain the kth-order () deformation equation subject to the initial conditions where

Due to the rule of solution expression and the linear operator , the right hand side of (23) should not contain the first harmonics and , because they can result in the so-called secular terms as and , respectively. To this end, let Solving (27), , and are determined step by step as increases.

Note that when , is essentially the right hand side of (10) with , 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 , and can be obtained by using the Newton-Raphson method. Importantly, (27) is always linear as , 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 . 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 into yields a simple illustration 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 and can be expressed as follows, respectively: where correspond to the jth component of vector , and , Differentiating times with respect to will result in a complicated term (i.e., ) 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 as an illustrative example to propose a means for expanding fractional functions. First of all, denote the denominator as where and . Taking the nth-order Taylor series of as , then one has Rewrite (31) as Equating the coefficients of results in Interestingly, (33) is always linear. That means it is rather easy to determine if are all known, .

#### 5. Numerical Examples

##### 5.1. Main Results

The system parameters under consideration are , , , , , , , and .

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 and .

Using analytical techniques developed for nonlinear dynamical systems, the linear flutter speed is found at [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.

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 , dashed lines represent HB1 results, and real lines denote numerical solutions.
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 , 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 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 , 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.

Figure 4: The phase planes of LCOs of system (1) with . Dots denote the 50th-order HAM solution with , real lines the numerical result, and dashed lines the HB1 one.
Figure 5: The time history responses of system (1) with . Dots denote the 50th-order HAM solutions with 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.

Table 1: Comparisons of the amplitudes and frequencies obtained by the HAM () with numerical solutions.
Figure 6: Residues of (6) () with HAM solutions attained with , (a): ; (b): ; (c): , 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 as = 0.07756360647090. Denote the discrepancy between the nth-order HAM frequency solution and as . 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 , , , and all approach to , while the one with does not. As decreases from −0.5 to −1 and further to −1.2, the convergent rate of the HAM solution increases. However, as decreases even a little, the convergent rate decreases (). It can even lead to the misconvergence of the HAM series (). 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, is a good choice.

Figure 7: The -curves with respect to , where .

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 and , the HAM series are disconvergent at . 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 () with the numerical solutions, where .
Table 3: Comparisons of the amplitudes and frequencies given by the HAM Pade approximations () with the numerical solutions, where .
##### 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 . 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.

Table 4: The CPU time needed in seeking the nth-order HAM approximations, the parameter values are and .
Figure 8: The ratio between the CPU times Tn and Sn, where and .

#### 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, . 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.

#### 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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
15. S. J. Liao, “On the homotopy analysis method for nonlinear problems,” Applied Mathematics and Computation, vol. 147, no. 2, pp. 499–513, 2004.
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.
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.
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.
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.
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.
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.
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.
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.
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.