摘 要:建立了求解弹性机翼阵风响应的计算方法。在计算中,通过采用数值方法求解三维非定常Euler方程来获得气动特性;采用模态叠加的方法考虑弹性影响,实现了流体力学和弹性力学的耦合计算。通过对刚性机翼在攻角突然增大的阵风作用下的响应历程计算和二维NLR7301翼型的极限环振荡计算,对计算方法进行了验证。此后在“1-cos”阵风响应的计算中考虑弹性效应影响,先是只考虑了结构变形的前三个基本模态,弹性机翼气动力响应的计算结果与刚性机翼的响应计算结果有比较大的区别,弹性机翼阵风响应的升力峰值低于刚性机翼,这与文献中的结果是一致的。最后在计算中考虑了高阶弹性模态,计算结果表明:考虑高阶模态后,机翼气动力计算结果的总体变化趋势与只考虑前三个模态时基本一致,但结果中出现了高频的波动,波动的频率与高阶模态本身的频率有关。
1 引言
动力响应分析在现代高速飞机的设计过程中占有非常重要的地位,所谓动力响应,是指飞机在飞行过程中由于外力作用随时间任意变化引起的动力效应,如大气紊流和阵风引起的交变力,由于机翼和发动机尾流等引起的结构激振作用等[1-3]。在飞机设计计算中必须综合考虑气动弹性效应,以防止出现不利的气动弹性现象。
目前,在工程上的动力响应计算中,大多采用的是数值分析的方法,文献[4]采用以格林函数法为基础的边界元方法,结合时间历程法,求解位势方程,得到了超声速机翼-尾翼对突风作用的动态响应;文献[5]采用数值求解N-S方程的方法得到了某弹性翼身组合体在离散阵风形式下的动力响应。本研究以截面均为NACA0012翼型的展弦比Λ=5的平直弹性机翼为研究对象,采用数值求解非定常Euler方程的方法来分析其在“1-cos”形式阵风作用下的动力响应问题。
2 计算模型
2·1 流动控制方程
对于绕机翼的三维非定常湍流流动,流动的主控方程采用三维非定常欧拉(Euler)方程,在计算坐标系下可写为
2·2 弹性计算模型
为推导弹性计算模型,如图1所示,设作用在飞机弹性面上的载荷为p,是空间和时间的函数,扰动将产生附加气动力,设为Δp。若记弹性面上单位面积的质量为m(x,y),弹性面的总面积为S,弹性面上各点偏离其平衡位置的z轴向位移为w(x,y,t),由此可得z轴向力、对x,y轴力矩的扰动平衡方程如下[1]:
式中τ称为虚拟时间。在虚拟时间步中,计算本文采用隐式近似AF分解的方法[8]计算推进;在方程空间导数的离散中,对流项采用了Roe的通量差分来离散[6]。
对于阵风响应问题,计算时需要通过引入“网格速度”来计入阵风的影响[3]。根据相对运动的思想,如果网格的速度为u,相当于计算域在网格不动的情况下整体都受到-u的来流作用。因此,图2中阵风的作用与整个计算域网格以速度wg向下运动是一致的,模拟图2中的阵风可以通过对整个计算网格附加向下运动的网格速度wg来实现,具体反映到方程中,即计算中需考虑坐标对时间的导
3·2 计算网格及网格变形处理方法
本文计算中的机翼为截面均为NACA0012、翼型的展弦比Λ=5的平直弹性机翼,计算时生成C-H型网格,网格维数为:159(流向)×40(法向)×33(展向)。由于在计算中机翼会产生弹性变形,所以本文采用了动网格技术,即机翼变形后,网格也随之产生变形,网格变形采用的是文献[9]中介绍的类似弹簧变形的算法,网格的变形速度应计入式(13)中。如图4所示给出了机翼在俯仰方向变形前后展向某一截面内的网格变形情况。
3·3 验证算例
为验证算法,先计算刚性机翼在攻角突然增大的阵风下的动态响应,即初始状态下机翼以速度V∞水平飞行,突然受到一上升气流作用,相当与翼型攻角突然增加Δα(Δα=arctg(wg/V∞))。这一问题是非定常空气动力学研究中的一个基本问题,因为机翼攻角的突然变化相当于攻角的一个阶跃信号,此时机翼的气动力变化情况通常被视为系统在阶跃信号下的响应,称为指示函数(IndicialFunction)。计算时机翼取截面均为NACA0012、翼型的展弦比Λ=5的平直机翼,并取wg/V∞=0·08,相当于机翼攻角突然增加约α=4°(0·08弧度)。对两个不同马赫数(M∞=0·3和M∞=0·8)下的阵风响应进行数值模拟。图5所示给出了这两个马赫数下的升力响应计算结果随无量纲时间s =2tV∞/c的变化历程,图中记为“Calc.(Present)”,并与文献[3]中的计算结果(图中记为“Calc.(Ref.[3])”)进行了比较,可以看出本文计算结果与文献[3]中的计算结果吻合,从而验证了计算方法中的非定常流动计算方法、阵风处理方法和机翼气动特性计算方法的可靠性。
为进一步验证弹性力学与流体力学耦合计算的可靠性,采用上述算法来计算二维翼型的极限环振荡LCO(Limiting Cycle Oscillation)这一典型的流固耦合问题。以文献[10]中的NLR7301翼型LCO计算为例,图6所示给出了计算模型示意图,翼型的弹性运动被简化成了如下的两个自由度的模型,一个自由度h对应于翼型的沉浮运动,另一个自由度α对应于翼型的俯仰运动:
这一简化的弹性计算模型相当于式(7)中广义坐标q1= h, q2=α的情况,式中的mh,Sα,Dh,Dα,Kh和Kα等为结构参数,具体取值参见文献[10]。
本文对M∞=0·753,α0=0·6°的算例进行了计算,图7所示给出了计算出的h随无量纲时间的变化历程,从中可以清楚地看出极限环振荡的发展过程。该变化趋势与文献[10]中的计算结果一致,由进一步计算振荡的周期值知计算出LCO的频率为35·23 Hz,与文献[10]中的实测结果32·85 Hz接近。该结果进一步验证了本文算法的有效性。
4 弹性机翼在“1-COS”型阵风作用下的响应计算结果
针对图3所示的“1-cos”型阵风,机翼仍取截面均为NACA0012、翼型的展弦比Λ=5的平直机翼,来流马赫数M∞=0·3,W0=0·05,先进行刚性机翼计算,图8中实线(记为“Rigid”)为此时的机翼升力系数响应计算结果,在阵风作用下,机翼升力先是升高,然后随着阵风速度的减小而逐渐降程低。接下来考虑弹性计算,先只取式(8)给出的三个模态,此时的计算结果如图8中的虚线(记为“Elastic-i(i=1~3)modes considered”)所示,可以看到,考虑弹性模态后,机翼升力值同样会在阵风作用下升高,但升力峰值小于刚性状态,这是由于机翼在z方向的弹性变形减弱了阵风的作用造成的。随着阵风速度的减小,此时机翼的升力同样会降低,并出现负值,但后期逐渐恢复,与文献[5]中的计算结果吻合。图9所示给出了俯仰力矩的计算结果比较。
接下来考虑式(9)中的高阶模态对阵风响应计算结果的影响,计算中模态无量纲频率ω4=1·57,阻尼系数ζ4=0。此时的升力和俯仰力矩响应计算结果如图8和图9中的点划线示(记为“Elastic-i(i =1~4)modes considered”)。从图中可以看到,考虑高阶模态后,机翼升力和俯仰力矩的计算结果的总体变化趋势与只考虑前三个模态时一致,但结果中出现了高频的波动,这显然与高阶模态的频率ω4有关。
从结果中可以看到,此时的广义坐标q4值有明显振荡,无量纲周期值大致为T 3·9,实际上,由式(7)中常微分方程的特性知,该周期值与ω4应有关系:T =2π/ω4。这一分析结果与文献[5]中的结果是定性一致的。
5 结 论
本文建立了求解弹性机翼阵风响应的计算方法。在计算中,通过采用数值方法求解三维非定常Euler方程来获得机翼的气动特性;采用模态叠加的方法来考虑弹性影响,实现了流体力学方程与弹性力学方程的耦合计算方法。
在此需要指出,本文中为简便起见,仅采用了抛物振型来作为机翼的高阶模态进行了处理,在实际工程计算中,应根据振动试验或专业软件的计算结果来给定机翼的高阶弹性模态振型和频率。
