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 [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 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.
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. [4–6] 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 [4–6].
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.
(a)
(b)
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.
(a)
(b)
(a)
(b)
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.
(a)
(b)
(c)
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.
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 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.
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.
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.
Appendix
We have the following Expressions of the Coefficients in (5):
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).