Efficient flutter prediction based on Harmonic Balance and V-g methods

Liu Nan，Bai Junqiang，*，Liu Yan，Hua Jun

(1.School of Aeronautics，Northwestern Polytechnical University，Xi’an 710072，China;2.China Aeronautical Establishment，Beijing 100012，China)

Efficient flutter prediction based on Harmonic Balance and V-g methods

Liu Nan1，Bai Junqiang1，*，Liu Yan1，Hua Jun2

(1.School of Aeronautics，Northwestern Polytechnical University，Xi’an 710072，China;2.China Aeronautical Establishment，Beijing 100012，China)

An efficient flutter prediction method is proposed and applied in the analysis of transonic flutter boundary and influences of structural parameters.Based on frequency domain flutter analysis V-g method，an artificial damping term is applied to aeroelastic system to maintain harmonic motion.Thus，the frequency domain structural dynamic equations can be obtained.The aerodynamic describing function matrices can be acquired efficiently by harmonic balance method at a variety of harmonic frequencies.Combining aerodynamic describing function matrices and frequency domain structural equations，the aeroelastic stability problem can be transferred into a generalized eigenvalue problem.It is demonstrated that the flutter boundary of this method conforms with high-fidelity time-marching method and the analysis efficiency improves significantly.Moreover，when the structural parameters change，new flutter boundary can be easily obtained by generalized eigenvalue analysis in this method，other than time-marching method which needs numerous CFD/CSD coupled simulations.

flutter;frequency domain;harmonic balance;aerodynamic describing function;timemarching method

0 Introduction

Flutter is a typical dynamic aeroelastic phenomenon，which could cause catastrophic results.Transonic flutter prediction is greatly important for the design of transport aircrafts and fighters.Nowadays， Computational Fluid Dynamics(CFD)has been widely applied in transonic flutter predictions［1］.However，it is still very time-consuming to calculate instantaneous response of aeroelastic system by time-marching methods.However，it is still very time-consuming to calculate instantaneous response of aeroelastic system by time-marching methods.

To improve analysis efficiency，reduced order model has been widely applied in flutter analysis.Weiwei Zhang et al［2］applied ARX reduced order model to construct aeroelastic equation in state space form，which shows similar results with CFD/CSD simulation，except for transonic condition.Walter A.Silva［3］created multiple-input multiple-output reduced order aeroelastic system based on Volterra series，which established the relationship between modal inputs and generalized aerodynamic forces(GAFs).However，the precision of unsteady aerodynamic forces based on Volterra series is very sensitive to amplitude of input signals and time steps［4］.Weiwei Zhang et al［5］also predicted flutter boundary by nonlinear ROM based on RBF neural network model.

Except for the time domain method，frequency domain method is another important approach of flutter analysis［6］.Aerodynamic describing functions can be obtained by exciting aeroelastic system with harmonic disturbance of small amplitudes.Then，flutter velocity and frequency of different structural parameters can be obtained by V-g or P-K method，which shows great efficiency，comparing with time-marching methods［7］.However，it takes excessive computational costs to calculate harmonic responses of aeroelastic system，because it needs to carry on unsteady flow simulation for each structural mode every harmonic frequency.

In addition，harmonic balance(HB)method shows superior efficiency in the simulation of periodic unsteady flow［8-10］，which is widely applied in the simulation of turbomachines［11-12］，helicopter rotors［13-14］，dynamic derivations［15］，et al.Therefore，a frequency domain flutter analysis method based harmonic balance is proposed in this paper，which takes efficiency and precision into account.

1 Time-marching method

Reynolds-Averaged Navier-Stokes(RANS) equation is the highest-fidelity flow solver in current engineering applications，which can be written as:

where Q and R are the vectors of flow variables and residuals respectively.Third-order Roe scheme is applied in spatial discretization，while approximate factorization scheme is used in implicit time-marching method.And S-A turbulence model is adopted.

Two-degree-of-freedom structural dynamic equation can be written as:

where

where

Eq.(3)can be solved by second-order Euler prediction-correction method:

The flow chart of time-marching approach of aeroelastic system is shown in Fig.1.

Fig.1 Flow chart of time-marching m ethod of aeroelastic system圖1氣動彈性系統時間推進求解的流程

2 Harmonic balance method

Harmonic balance method is an efficient approach of periodic unsteady flow field.Because of the periodicity，the vector of flow variables and residuals can be deployed by Fourier series:

The time derivation of Q is:

Substituting Eq.(5)into Eq.(1)and rearranging，we can get a pair of algebraic equations:

Eq.(7)can be written into matrix form:

where

Eq.(8)is the governing equation of original harmonic balance method.In this method，the calculation of residual terms by Fourier transformation is very complex and tedious，which is as follows:

To resolve this problem，sub-time levels were introduced by Hall et al［11］to convert Q back to time domain through Fourier transformation:

where

{ti，i=1，2，...，2NH}are equally distributed sub-time levels in a cycle，and ti=2π(i－1)/2NH.

Substituting Eq.(10)into Eq.(8):

Both sides in Eq.(11)are multiplied by E－1:

To achieve an easy approach，the residual terms are simplified by:

where

In this step，nonlinear residual terms can produce terms whose frequencies are higher that NHω，which are included in low-frequency terms by Eq.(13).This rough simplification may produce non-physical solutions in strong nonlinear problems［16］.Therefore，the number of harmonics should be chose adequately to guarantee that high-frequency harmonic terms approach zero.Finally the governing equation of harmonic balance method can be obtained:

In order to solve Eq.(14)by time-marching method，pseudo time derivation is added:

Eq.(15)is the same as original RANS Eq.(1)，except for its middle term(which is also called harmonic balance source term).Therefore，harmonic balance solver can be obtained from original steady RANS solver.

Then the harmonic balance method is validated by NACA64A010 airfoil under forced harmonic motion［17］.The angle of attack changes as follows:

The flow conditions of two test cases are shown in Table 1.The comparison of unsteady aerodynamic forces between harmonic balance and time-marching methods are illustrated in Fig.2 and Fig.3.It is demonstrated that the unsteady forces calculated by harmonic balance and time-marching methods coincide with each other，which validates the precision of harmonic balance method.The time costs of time-marching and harmonic balance methods are shown in Table 2.

Table 1 Flow condition of NACA64A010 under forced motions表1 NACA64A010強迫運動算例的來流條件

Fig.2 Com parison of unsteady aerodynam ic forces between time-marching and harmonic balance methods in CT2 case圖2 CT2算例時間推進和諧波平衡方法計算結果對比

Fig.3 Com parison of unsteady aerodynam ic forces between time-marching and harmonic balance methods in CT5 case圖3 CT5算例時間推進和諧波平衡方法計算結果對比

Table 2 Time costs of time-marching and harmonic balance methods(unit:s)表2時間推進和諧波平衡方法計算時間對比(單位:s)

3 Frequency domain solver

The corresponding lift and pitch moment coefficients can be written as:

Therefore，the structural Eq.(2)can be transferred into frequency domain:

When the aeroelastic system is undergoing single harmonic motion，the plunge and pitch displacements can be written as:

Through aerodynamic describing function matrix F，the aerodynamic forces can be easily obtained by:

where

The aerodynamic describing functions can be obtained through exciting aeroelastic system by harmonic disturbance of small amplitude.Harmonic Balance method is applied to acquired aerodynamic describing functions efficiently.

However，the aeroelastic systemmay exhibit different behaviors(convergent oscillation，simple harmonic or divergent oscillation)in different freestream condi-tions.Therefore，an artificial damping coefficient j g K can be added into the aeroelastic system(19)to maintain harmonic motion:

The stability of aeroelastic system can be judged by the sign of artificial damping coefficient.Eq.(21)can be rearranged into a generalized complex eigenvalue problem:

where

At a fixed frequency ω，the matrix A will only change with freestream velocity.Different eigenvalue λ can be obtained in various velocity.When the real part of eigenvalue λReis equal to 1.0，the imaginary part of eigenvalue λImis the damping coefficient at a certain velocity.Thus we can get the relationship between damping coefficients and frequencies(or velocities).When damping coefficient is equal to 0，the frequency and velocity is flutter frequency and velocity respectively.

The procedure of frequency domain flutter analysis is as follows:

1)Calculate aerodynamic forces when the system are undergoing harmonic motionsat different frequencies;

2)Obtain aerodynamic transfer function matrix;

3)Construct generalized complex eigenvalue problem by aerodynamic transfer function matrix and structural dynamic equation;

4)Get the freestream velocity which makes the real part of eigenvalue λRe=1;

5)Obtain the relationship between damping coefficients and frequencies/velocities;

6)Acquire the flutter velocity and frequency.

4 Flutter analysis of Isogai wing

The Isogai Case A wing［18-19］is a standard configuration of flutter analysis.The structural parameters of this case are:xα=1.8，r2α=3.48，a=－2.0，ωh=ωα=100 rad/s，μ=60.The inviscid flutter velocity boundary shows a“S”shape in transonic Mach number，while the viscous flutter velocity boundary does not.

4.1 Aerodynamic describing function matrix

The aerodynamic describing function matrix is calculated byharmonic balance method to improve efficiency.Take h－/b=0.001 and ω=0.10 as an example.Tables 3 and 4 are the lift describing functionsand pitch moment describing functionsat two Mach number (Ma=0.60 and Ma=0.87)respectively.

Table 3 Aerodynam ic describing functions at Ma=0.60表3 Ma=0.60時的氣動描述函數

Table 4 Aerodynam ic describing functions at Ma=0.87表4 Ma=0.87時的氣動描述函數

It is illustrated that first-order harmonic is enough for subsonic case and more harmonics are necessary for transonic case to ensure the precision of aerodynamic describing functions.

4.2 Flutter boundary

Fig.4 Com parison of inviscid flutter boundaries between time-marching and frequency domain methods圖4時間推進方法和頻域法得到的無粘顫振邊界對比

Fig.5 Comparison of viscous flutter boundaries between time-marching and frequency domain methods圖5時間推進方法和頻域法得到的粘性顫振邊界對比

Fig.6 Inviscid damping coefficient versus speed index at Ma=0.88圖6 Ma=0.88時無粘阻尼系數隨來流速度的變化趨勢

The inviscid and viscous flutter boundaries are shown in Figs.4 and 5.It is demonstrated that the difference between time-marching and frequency domain methods is insignificant，even in trmansonic range. The inviscid damping coefficient at Ma=0.88 is exhibited in Fig.6，from which we can obtain the reason of“S”shape flutter boundary.The damping coefficient of plunge mode exceeds zero at a low flow velocity.But as the velocity increases，the damping coefficient returns back and falls below zero at a moderate velocity.Then the damping coefficient of pitch mode exceeds zero at a high velocity.

4.3 Influence of structural parameters

If the structural parameters undergo some change，flutter boundary can be easily obtained by generalized eigenvalue analysis in frequency domain method，other than time-marching method，which must carry on numerous CFD/CSD simulations.Therefore，the influences of mass ratio and frequency ratio on viscous flutter boundary are investigated below，based on the frequency domain method.

The influences of mass ratio μ and frequency ratio ωh/ωαare shown in Fig.7 and 8 respectively.It is illustrated that higher mass ratio is preferable in subsonic region，while lower mass ratio is better in transonic region of Isogai A case.The frequency ratio also has a great impact on flutter boundary.When the frequency ratio is close to 1.0，the flutter velocity decreases significantly.

Fig.7 Influence of mass ratio on flutter velocity圖7質量比對顫振速度的影響

Fig.8 Influence of frequency ratio on flutter velocity圖8頻率比對顫振速度的影響

5 Conclusion

1)Comparing with time-marching method，harmonic balance method is an efficient approach of periodic unsteady flow;

2)A frequency domain flutter analysis method based on harmonic balance is proposed，which shows high efficiency and precision;

3)When the structural parameters undergo some change，new flutter boundary can be obtained quickly by generalized eigenvalue analysis in the proposed frequency domain method，other than time-marching method，which must carry on a series of CFD/CSD coupled simulations.

［1］Schuster D M，Liu D D，Huttsell L J.Computational aeroelastic: success，progress，challenge［J］.Journal of Aircraft，2003，40 (5):843-856.

［2］Zhang Weiwei，Ye Zhengyin.Effect of control surface on airfoil flutter in transonic flow［J］.Acta Astronautica，2010，66:999-1007.

［3］Silva W A.Simultaneous excitation of multiple-input multipleoutput CFD-based unsteady aerodynamic systems［C］//48th AIAA/ ASME/AS-CE/AHS/ASC Structures，Structural Dynamics，and Materials Conference.Honolulu，Hawaii.AIAA 2007-1988.

［4］Raveh D E.Reduced-order models for nonlinear unsteady aerodynamics［J］.AIAA Journal，2001，39(8):1417-1429.

［5］Zhang Weiwei，Wang Bobin，Ye Zhengyin，et al.Efficient method for limit cycle flutter analysis by nonlinear aerodynamic reducedorder models［J］.AIAA Journal，2012，50(5):1019-1028.

［6］Ueda T，Dowell E H.Flutter analysis using nonlinear aerodynamic forces［C］//Proceedings of the 23rd AIAA/ASME/ASCE/AHS Structures，Structural，Dynamics，and Materials Conference，New Orleans，LA.AIAA 82-0728-CP.

［7］He S，Yang Z，Gu Y.Transonic limit cycle oscillation analysis using aerodynamic describing functions and superposition principle［J］.AIAA Journal，2014，52(7):1393-1403.

［8］Woodgate M A，Badcock K J.Implicit harmonic balance solver for transonic flow with forced motions［J］.AIAA Journal，2009，47 (4):893-901.

［9］Im D，Kwon J，Park S.Periodic unsteady flow analysis using a diagonally implicit harmonic balance method［J］.AIAA Journal，2012，50(3):741-745.

［10］Blanc F，Roux F X，Jouhaud J C.Harmonic-balance-based codecoupling algorithm for aeroelastic systems subjected to forced excitation［J］.AIAA Journal，2010，48(11):2472-2481.

［11］Hall K C，Thomas J P，Clark W S.Computation of unsteady nonlinear flows in cascades using a harmonic balance technique［J］.AIAA Journal，2002，40(5):879-886.

［12］Weiss J M，Subramanian V，Hall K C.Simulation of unsteady turbomachinery flows using an implicitly coupled nonlinear harmonic balance method［C］//Proceedings of ASME Turbo Expo 2011.Vancouver，British Columbia.GT 2011-46367.

［13］Ekici K，Hall K C，Dowell E H.Computationally fast harmonic balance methods for unsteady aerodynamic predictions of helicopter rotors［C］//46th AIAA Aerospace Sciences Meeting and Exhibit.Reno，Nevada.AIAA 2008-1439.

［14］Xu Jianhua，Song Wenping，Wang Long.Application of harmonic balance method in forward flight simulation for helicopter rotors［J］.Acta Aerodynamica Sinica，2013，31(5):546-553.(In Chinese)許建華，宋文萍，王龍.諧波平衡法在旋翼前飛繞流數值模擬中的應用研究［J］.空氣動力學學報，2013，31(5):546-553.

［15］Ronch A D，McCracken A J，Badcock K J.Linear frequency domain and harmonic balance predictions of dynamic derivatives［J］.Journal of Aircraft，2013，50(3):694-707.

［16］Liu L，Thomas J P，Dowell E H，et al.A comparison of classical and high dimensional harmonic balance approaches for a Duffing oscillator［J］.Journal of Computational Physics，2006，215:298-320.

［17］Davis S S.NACA64A010(NASA AMES model)oscillatory pitching［R］.AGARD Report No.702.

［18］Isogai K.On the transonic-dip mechanism of flutter of a sweptback wing［J］.AIAA Journal，1979，17(7):793-795.

［19］Alonso J J，Jameson A.Fully-implicit time-marching aeroelastic solutions［C］//32nd Aerospace Sciences Meeting＆Exhibit.Reno，NV.AIAA-94-0056.

V211.3

A

0258-1825(2016)05-0631-07

基于諧波平衡法和V-g法的高效顫振預測分析

劉南1，白俊強1，*，劉艷1，華俊2
(1.西北工業大學航空學院，陜西西安710072;2.中國航空研究院，北京100012)

提出了一種高效的顫振預測分析方法，將之應用于跨聲速顫振邊界分析及結構參數影響研究中。本方法基于頻域顫振分析V-g方法，為氣彈系統提供一定的人工阻尼使之保持簡諧運動狀態，從而將結構動力學方程轉換到頻域內。然后通過高效的諧波平衡法得到一系列簡諧運動頻率下的氣動力描述函數矩陣，結合頻域結構方程，將氣彈系統的穩定性問題轉化為廣義特征值求解。結果表明:本方法計算得到的顫振邊界與高精度的時間推進方法非常吻合，分析效率有了明顯的提升，而且當結構參數發生變化后，只需進行若干次廣義特征值求解即可得到新的顫振邊界，無需像時間推進方法一樣開展大量的氣動/結構耦合數值模擬。

顫振;頻域;諧波平衡;氣動力描述函數;時間推進方法

2015-09-13;

2015-11-10

劉南(1989-)，男，碩士，主要研究方向:飛行器氣動設計，氣動彈性.E-mail:revolution890926@163.com

白俊強*(1971-)，教授，博士生導師.E-mail:junqiang@nwpu.edu.cn

劉南，白俊強，劉艷，等.基于諧波平衡法和V-g法的高效顫振預測分析(英文)［J］.空氣動力學學報，2016，34(5):631-637.

10.7638/kqdlxxb-2015.0178 Liu N，Bai J Q，Liu Y，et al.Efficient flutter prediction based on Harmonic Balance and V-g methods［J］.Acta Aerodynamica Sinica，2016，31(5):631-637.

10.7638/kqdlxxb-2015.0178