带控制截面机翼结构基于非线性能量阱的颤振抑制*

陈恒 王扬渝 金江明

(浙江工业大学机械工程学院, 杭州 310014)

引言

机翼颤振是飞行中常见可能带来灾难性后果的气动弹性现象,属于自激振动.机翼的颤振主要是由机翼的结构非线性和气动非线性相耦合引起.由于结构非线性和气动非线性的存在,振动幅值限制在一定范围内,其主要形式为极限环振动.极限环振动是气动弹性线性的重要问题,造成飞行器结构的危害,如F-16和F/A18战机在超音速飞行时均会产生极限环振动[1].对于极限环颤振,人们在理论与实验研究方面已经做了不少研究.Gilliat等[2]研究了结构非线性和失速状态下气动非线性的运动,得到极限环颤振.Jones和Roberts[3]识别了带双线性的非线性气动弹性结构的极限环振动.Tang和Dowell[4]研究了带挂装的机翼产生间隙非线性结构的振动响应.赵永辉等[5]进行了大展弦比夹芯翼大攻角颤振分析.崔鹏等研究了新型运输机的颤振特征,发现翼梢装置对机翼颤振的不利影响[6].周秋萍和邱志平研究了外挂连接具有初偏间隙非线性的机翼颤振问题,并对其颤振进行区间分析得到颤振边界曲线[7].

机翼结构往往带有间隙非线性、双线性和碰摩等强非线性特征,存在分岔和混沌等动力特征,减振难度显著增大.目前其减振方法包括主动控制和被动控制形式.主动控制需测量振动状态,通过控制策略反馈到控制器,实现对振动的抑制[8].对于非光滑振动系统,线性刚度吸振器不再适用[9],Vakakis等[10]提出非线性靶能量转移的概念,耦合非线性能量阱(nonlinear energy sink,以下简称NES),实现非线性被动控制.McFarland等[11]实验验证NES在强非线性系统中的减振效果.Lee等[12]使用NES成功抑制了范德波尔方程的极限环振动,并应用NES技术对机翼颤振进行抑制,研究了一个二自由度机翼结构通过NES的减振效果,得到三种不同的抑制机理[13].

另外,研究气动弹性系统在不同流速下的整体振动特征,常用的方法是计算其分岔图.对于预测和分析机翼结构分岔的方法主要有范式法,中心流方法[14],胞映射法[15],谐波平衡法[16]等.目前,数值法在气动弹性系统分析和分岔计算中等到广泛应用和发展,如G. Dimitriadis[17]用数值法计算了带间隙非线性的气动弹性结构的分岔,沈华勋[18]用数值法研究柔性机翼的稳定性问题.而对于气动弹性系统的强非线性,不同系统具有很强的独特性,难以用一种通用的方法实现分岔的分析.本文带控制界面的机翼存在采用双线这种非光滑的刚度,通过耦合NES实现其颤振抑制,通过基于弧长数值连续的方法构造分岔图,验证其减振效果.

1 二元机翼基于非线性能量阱的减振模型

考虑不可压缩流中的具有浮沉、俯仰和控制截面转动三个自由度的二元机翼颤振系统模型(如图1(a)所示,其中NES质量m=0).图l中e是空气动力焦点到弹性轴的距离;b是半弦长;a是从半弦点到弹性轴的距离;c是从半弦点到控制截面的距离.h,α和β分别为浮沉位移,俯仰位移和控制截面的转动位移.Kh和Kα表示h和α方向的弹性系数;c1和c2表示h和α方向的非线性刚度系数;Kβ为控制截面转动的刚度,用双线性表示(如图1(b)所示).U为气流速度.运动方程为:

(1)

其中M为机翼质量,Sα和Sβ分别为机翼主体和控制截面的质量不平衡,Iα和Iβ为机翼主体和控制截面的转动惯量.

图1 带控制截面二元机翼模型(a)机翼截面; (b)机翼控制截面双线性扭转刚度Fig. 1 Two-dimensional wing with control surface model(a) Schematic of the typical wing sections with control-surface; (b) Bilinear torsional stiffness of wing control-surface

(2)

当转角位小角度时力矩Mα≈eLα,Mβ≈eLβ.

当振动结构耦合本征非线性刚度的结构会形成能量单向不可逆转移,达到减振效果[10],非线性能量阱结构本质特征是要求耦合结构的刚度为本征非线性,如三次方非线性刚度的NES结构可以形成非线性能量阱结构.在二元机翼模型上增加非线性能量阱结构(单自由度的三次方非线性刚度的NES结构),实现对机翼颤振的抑制(如图1(a)所示,其中NES质量m≠0).在原来三自由度的基础上增加一个自由度,形成四自由度系统.其运动方程如下:

(3)

方程(3)的无量纲形式如下:

ελ(y′-δα′-v′)+C(y-δα-v)3=0

δελ(δα′+v′-y′)+C(y-δα-v)3=0

εv″+ελ(v′+δα′-y′)+C(v+δα-y)3=0

(4)

其中δ=0.3.

2 减振前后数值响应的整体特征

为了描述机翼结构减振前后的整体响应特征,在不同风速下分别数值模拟系统响应,记录响应稳定后的极值形成峰值-峰值图.构造峰值-峰值图的具体步骤为,给定初始条件,风速从[0.8 1.0]范围内200等分,在每个不同风速下数值模拟系统的响应,在位移风速坐标下记录τ∈[8000 8500]内响应的极值点(通过数值响应可认定若系统为稳定,则当时间τ>8000时响应已趋于稳定).考虑初始大扰动和初始小扰动两种初始值情况构造峰值-峰值图,小扰动初值为在浮沉自由度y=0.001,y′=0,其他自由度α=β=v=0,α′=β′=v′=0,其数值响应的峰值-峰值图如图2所示(黑色点为减振前响应极值,红色点表示减振后响应极值).大扰动初值为在浮沉自由度y=0.1,y′=0,其他自由度α=β=v=0,α′=β′=v′=0,其数值响应的峰值-峰值图如图3所示(黑色点为减振前响应极值,蓝色点表示减振后响应极值).

减振前系统,对于大初值和小初值始扰动,其峰值-峰值图保持基本相同(图2和图3黑色点).初步反应了系统在耦合NES之前的分岔结构,可以发现当u*=0.88时系统发生Hopf分岔.由峰值-峰值图可以发现当风速小于0.903时,控制截面振动幅值小于0.3系统响应为光滑振动,如果风速大于0.903系统导致非光滑振动.当风速介于0.903和0.918之间时响应为周期运动形成极限环.当风速介于0.918与0.922之间以及0.923与0.941之间将出现混沌现象.当风速大于0.903除上述区域外系统都将发生准周期运动.

图2 小初始值下减振前后峰值-峰值图对比Fig. 2 Peak-peak diagram of the sytem with/without NESat small initial conditions

图3 大初始值下减振前后峰值-峰值图对比Fig. 3 Peak-peak diagram of the sytem with/without NES atlarge initial conditions

图2表示小扰动初值情况下减振前后系统的峰值-峰值图对比.从图2可以看出当初值为小扰动时,系统减振效果非常显著,减振后Hopf分岔点为0.892明显大于减振前的0.88,这表明如果风速小于0.892,增加NES后系统可以完全抑制颤振.如果风速介于0.892和0.908之间,则减振后系统为极限环运动,相对于减振前其极限环的幅值大大减少.当风速值大于0.908则系统减振后为准周期振动,但是其响应的幅值远小于减振前.

图3表示大扰动初值情况下减振前后系统的峰值-峰值图对比.从图3可以看出当初值为大扰动时,Hopf分岔点同样为0.892大于减振前的0.88,这表明如果风速小于0.892,增加NES后系统可以完成抑制颤振.风速介于0.892和0.919之间,则减振后系统为极限环运动,相对于减振前其极限环的幅值降低.当风速值大于0.919则系统减振前后振幅无明显变化,减振效果不佳.

3 减振前后系统分岔图的计算和对比

数值连续法用于求解带一个或多个参数的非线性问题[19],在求解气动弹性力学问题中已有相应的软件包如AUTO[20]和MatCont[21],Roberts等[22]利用其求解双线性问题的解,Dimitriadis求解了间隙非线性机翼振动的分岔问题[18].数值连续法主要用于求解如下形式的非线性代数方程的解:

f(x,λ)=0

(5)

其中f表示非线性函数向量,x表示未知数向量,λ表示变量.给定方程在λ0时已知解x0,利用数值连续法寻找在变量λ0+Δλ时的解x0+Δx,其中Δλ为微小变化.

在本文中需要求解的是常微分方程,其形式为:

(6)

其中x=[x1(t),…,xn(t)]T,f=[f1,…,fn],λ为变量参数.使用数值连续法需要把常微分方程转化为代数方程.基于系统的周期解分析分岔问题,把方程(5)转化为一系列的非线性代数方程:

F(x(t),λ)=0

(7)

其中F表示非线性函数.在特定的参数λ下,给定一个假设的初始值x0,基于牛顿迭代可以得到其精确解x.

(8)

x1=x0+Δx

(9)

其中Δx表示对x0的改进,经过多次迭代可以得到精确值.

把常微分方程转化为代数方程常用的方法为打靶法,这实际上是一种边界值方法.若方程的解x(t)是周期解,那么数值连续法就转化为求解如下的边界值问题.

x(0)=x(T)

(10)

常微分方程转化为:

F(x,T,λ)=x(0)-x(T)=0

(11)

给定初始值x(0),通过数值积分得到x(T),计算函数F的值.对于求解微分方程的响应可以使用龙格库塔等数值方法.假设系统参数为λ0时,响应为x0(t)周期为T0,那么当参数为λ=λ0+Δλ时,

=-F(x0(0),T0,λ)

(12)

(13)

方程(11)表示n个方程有n+1个未知数.在这个方程中,使用相位固定法,最简单的就是假设初值中一个元素始终为0,比如x1(0)=0,这样未知数就降为n个.

(14)

其中δx是微小量,δxi表示为n×n矩阵δ×I的第i列.当x(0)的第一个元素均为0时,计算雅可比矩阵从i=2开始,结果矩阵为n×(n-1).

上述方法为基本的数值连续法,称为自然参数连续.对于非线性问题,基于自然参数的连续法求解中其解的分支可能出现折返,交叉甚至消失.出现这类情况会产生奇异点,自然参数的连续法将失效.更加有效的方法是采用弧长连续(arclength continuation).弧长数值连续选用参数s不再采用系统参数λ,对应的步长采用Δs.当数值连续进行到分岔附近则需要减少步长,以免丢失真实的解.当远离分岔时若采用很小步长,则影响方法的计算效率.故需采用变步长的策略[15].弧长为x(0),T和λ对应s的改变速率.

(15)

把Δλ看成另外一个未知数,可以得到:

(16)

(17)

此时系统为n+1各方程,可以求解n+1个未知数.

(18)

(19)

对于周期解,另外一个重要性质是其稳定性.对其初始条件施加一个微小扰动,运动偏离原来周期运动则为不稳定,反之为稳定.对于一个周期为T的运动轨迹x(t,x(0)),其在t=T时的雅可比矩阵为:

(20)

对初始条件施加微扰动Δx(0),对扰动后的解x(T,x(0)+Δx(0))泰勒展开:

Δx(T)=ΦT(x(0))Δx(0)+o(‖Δx(0)‖2)

(21)

其中Δx(T)=x(T,x(0)+Δx(0))-x(T,x(0)),表示经过一个周期后运动的偏离量,如果经过m个周期.

Δx(mT)=[ΦT(x(0))]mΔx(0)+o(‖Δx(0)‖2)

(22)

周期运动的稳定性取决于矩阵ΦT的2n个特征值,即Floquet乘子.若Floquet乘子至少有一个值大于1,则周期运动为不稳定,反之则系统稳定.

给定初始的一个周期解,基于上述的弧长数值连续,可以得到系统在不同参数下周期解,及其稳定性分析,构造系统减振前后的分岔图,如图4所示.图6中红线和绿线表示减振前系统的分岔图,其中红线为稳定解,绿线为不稳定解.蓝线及灰色虚线为减振后的分岔图,蓝线为稳定解,虚线表示不稳定解.图5为减振后分岔图和峰值-峰值图的对比,图5左边一栏为下分支与数值响应的峰值对比,图5右侧一栏为上分支和数值响应的对比.由图5可得当分岔图为稳定分支则其幅值和数值响应的峰值非常吻合,其不稳定分支可出现两种情形即准周期运动(小扰动初值风速大于0.908)和混沌现象(大扰动初值风速大于0.957).由图4可以得到增加NES后系统Hopf分岔点出现在μ*=0.892,这和峰值-峰值图的结论一致.在Hopf分岔点出现以前表示系统可以完全抑制颤振.当0.892<μ*<0.908时,分岔图中存在一个稳定下分支,当系统响应为小扰动时运动会趋于下分支.另外0.894<μ*<0.90时存在一个稳定的上分支,当系统响应为小扰动时响应会趋于上分支.此时减振前后均为极限环运动,减振措施可以使得极限环的幅值更小.由于系统双线性存在,在不稳定分支对于初值为小扰动才生准周期运动,对于大扰动初值则出现混沌现象.

图4 减振前后分岔图对比Fig. 4 Comparison of Bifurcation diagram for the systems with/without NES

通过峰分岔图计算对比,结合峰值-峰值图,发现添加NES作为减振措施当系统为小扰动时可以实现很好的颤振抑制效果.对于大扰动的情景μ*<0.892时可以完全抑制颤振,0.892<μ*<0.918时可以降低极限环振动幅值.而当μ*>0.918则无明显减振效果抑制失败,从分岔图及峰值-峰值图可以发现响应幅值与减振前基本保持相同,不会导致减振失败的负面效应.机翼颤振抑制大致可分为三种情况:(a)风速小于Hopf分岔点出现之前颤振完全抑制,图6所示为抑制前后的时间序列对比;(b)0.892<μ*<0.908时,减振前后分岔图均为稳定的极限环运动,振动幅值大幅降低,其时间序列对比如图7所示.(c)μ*>0.908初始小扰动,减振前后分岔图均不稳定分支.如图8所示,减振前响应为非光滑的准周期振动,而减振后为光滑小幅值准周期运动.

图5 分岔图与峰值-峰值图对比Fig. 5 Comparison between Peak-Peak diagram and bifurcation diagram

图6 μ*=0.888,初始小扰动减振前后时间序列对比Fig. 6 Comparison of time series for the systems with/without NES when μ*=0.888 with small initial conditions

图7 μ*=0.895,初始大扰动减振前后时间序列对比Fig. 7 Comparison of time series for the systems with/without NES when μ*=0.895 with large initial conditions

图8 μ*=0.935,初始小扰动减振前后时间序列对比Fig. 8 Comparison of time series for the systems with/without NES when μ*=0.935 with small initial conditions

4 结论

带控制截面的机翼结构简化为三自由度二元刚性气动弹性模型,机翼结构耦合非线性能量阱,实现对机翼颤振的抑制.建立机翼减振前后的运动方程,通过数值模拟构造减振前后峰值-峰值图,初步反应其整体的减振效果.减振后Hopf分岔点为0.892明显大于减振前的0.88,这表明如果风速小于0.892,增加NES后系统可以完成抑制颤振.小扰动情况或者大扰动初值在风速低于0.919时,系统减振效果显著.

通过弧长数值连续法和结合Floquet算子构造减振前后的分岔图,并研究了其稳定性.通过Hopf分岔点的推迟及对比减振前后稳定和不稳定分支的存在,验证了三类不同形式的颤振抑制:颤振完全抑制(风速小于Hopf分岔点出现之前),极限环部分抑制(减振前后分极限环运动幅值大幅降低),准周期运动抑制(减振前响应为非光滑的准周期振动,而减振后为光滑的小幅值准周期运动).

1Bunton R, Denegri Jr C. Limit cycle oscillation characteristics of fighter aircraft.JournalofAircraft,EngineeringNote, 2000,37:916～918

2Jones D P, Roberts I, Gaitonde A L. Identification of limit cycles for piecewise nonlinear aeroelastic systems.JournalofFluidsandStructures, 2007,23:1012～1028

3Tang D, Dowell E H. Experimental and theoretical study of gust response for a wing-store model with freeplay.JournalofSoundandVibration, 2006,295:659～684

4赵永辉,胡海岩. 大展弦比夹芯翼大攻角颤振分析. 振动工程学报, 2004,17(1):25～30 (Zhao Y H, Hu H Y. Flutter Analysis of a high-aspect-ratio sandwich wing under large angle of attack.JournalofVibrationEngineering, 2004,17(1):25～30 (in Chinese))

5崔鹏,韩景龙. 新型运输机机翼的颤振特性分析. 振动工程学报, 2011,24(2):192～199 (Cui P, Han J L. Flutter analysis of new transport—type wings.JournalofVibrationEngineering, 2011,24(2):192～199 (in Chinese))

6周秋萍,邱志平. 机翼带外挂系统极限环颤振的区间分析. 航空学报, 2010,31(3):514～518 (Zhou Q P, Qiu Z P. Interval analysis for limit cycle flutter of a wing with an external store.ActaAeronauticaetAstronauticaSinica, 2010,31(3):514～518 (in Chinese))

7Tanaka N. Impact vibration control using a semi-active damper.JournalofSoundandVibration, 1992,158(2):277～292

8Tanaka N. Experiment of the dynamic damper with a preview action.JSMEInternationalJournal, 1987,53(487):650～655

9Lee Y S, Kerschen G, Vakakis A F, et al. Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment.PhysicaD, 2005,204(1):41～69

10 Kerschen G, Vakakis A F, Panagopoulos P N, et al. Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment.PhysicaD, 2005,204(1):41～69

11 Lee Y S, Vakakis A F, Bergman L A, et al. Suppression of limit cycle oscillations in the van der pol oscillator by means of passive nonlinear energy sinks (NESs).StructuralControlandHealthMonitoring, 2006,13:41～75

12 Lee Y S, Vakakis A F, Bergman L A, et al. Suppressing aeroelastic instability using broadband passive targeted energy transfers, Part 1: Theory.AIAAJournal, 2007,45(3):693～711

13 Vio G A, Cooper J E. Limit cycle oscillation prediction for aeroelastic systems with discrete bilinear stiffness.InternationalJournalofAppliedMathematicsandMechanics, 2005,1:100～119

14 Liu L, Wong Y S, Lee B H K. Application of the centre manifold theory in non-linear aeroelasticity.JournalofSoundandVibration, 2000,234(4):641～659

15 Levitas J, Weller T, Singer J. Poincare-like simple cell mapping for non-linear dynamical systems.JournalofSoundandVibration, 1994,176:641～662

16 Raghothama A, Narayanan S. Non-linear dynamics of a two-dimensional airfoil by incremental harmonicbalance method.JournalofSoundandVibration, 1999,226:493～517

17 Dimitriadis G. Bifurcation analysis of aircraft with structural nonlinearity and freeplay using numerical continuation.JournalofAircraft, 2008,45(3):893～905

18 沈华勋. 飞翼式柔性飞机纵向动力学建模与稳定性分析. 动力学与控制学报, 2016,14(3):241～246 (Sheng H X. Longitudinal danamic modeling and stability analysis of very flexible flying wings.JounalofDynamicsandControl, 2016,14(3):241～246 (in Chinese))

19 Kubicek M. Dependence of solution of nonlinear systems on a parameter.ACMTransactionsonMathematicalSoftware2, 1976:98～107

20 Doedel E J, Champneys A R, Fairgrieve T F, et al. AUTO 97-AUTO 2000: Continuation and bifurcation software for ordinary differential equations (with HomCont): User′s Guide, Technology Report, Concordia University, Montreal, Canada, 2000

21 Govaerts W, Dhooge A, Kuznetsov Y A. MATCONT: A matlab package for numerical bifurcation of ODEs.ACMTransactionsonMathematicalSoftware, 2003,29(2):141～16

22 Roberts I, Jones D P, Lieven N A J, et al. Analysis of piecewise linear aeroe-lastic systems using numerical continuation.ProceedingsoftheIMechEPartGJournalofAerospaceEngineering, 2002,216(1):1～11