带间隙非线性的机翼操纵面颤振特性研究

张易明， 何绪飞， 艾剑良

(复旦大学 航空航天系，上海 200433)

颤振是一种典型的气动弹性动稳定性问题，它从气流中持续不断的吸收能量，从而微小的扰动都会导致振荡幅值持续增大，使得某一个或多个运动模态克服结构和气动阻尼，最终不可抑制地达到破坏状态。颤振的临界速度范围通常比较狭窄，造成致命破坏的时间往往也很短。因此，对民用机的颤振特性开展分析具有重要意义。

操纵面间隙非线性是当今大型飞机颤振非线性研究的一个重点问题。飞机操纵面和主翼面之间制造公差、不恰当的安装及旋转部件之间机械连接会导致间隙非线性(freeplay nonlinearity)，即产生一个小范围的零刚度区域。当操纵面间隙增大到一定范围，飞机就可能在气动弹性稳定性包线内就发生颤振或极限环振荡(limit cycle oscillation, LCO)。现代飞机设计越来越多地采用非质量平衡的操纵面，依靠作动器的刚度和阻尼防止发生颤振和LCO。

然而，由于间隙导致气弹稳定性降低的现象仍有发生，进而导致颤振或LCO的情况。LCO的不利影响包括不舒适的飞行条件、操纵面振动和疲劳失效。

目前，国内外许多学者已开展了非线性间隙对操纵面颤振特性的影响方面的研究。Fichera等[1]探讨了操纵面间隙的影响，建立了带控制链间隙 T 型尾翼的气动弹性风洞模型。并建立了两个数值模型的气动弹性风洞模型：①通过驱动器反馈回路的集中质量单元来表征非线性，以构建状态空间系统来描述尾翼动力学；②用高阶谐波平衡法描述。Morino等[2]采用基于计算流体动力学的气动弹性分析方法模拟由间隙引起的操纵面LCO。Wang等[3]研究了不可压缩流动中具有间隙和三次俯仰刚度的典型二元翼型的非线性响应，并将微分变换法应用于气动弹性系统。Kholodar[4]提出了一种典型的喷气飞机气动弹性稳定包线的间隙分析方法。在副翼、升降舵、方向舵和外挂的连接处分别模拟了旋转间隙。Pereira等[5]采用分岔和HOS(Righer-order Statistic)法分析了带俯仰刚度硬化非线性和操纵面刚度间隙非线性的三自由度翼面的动力学响应。Ning等[6]采用自由界面模态综合法，通过连接关系的扩展，建立了适合于非线性气动弹性分析的控制方程。通过数值模拟和颤振风洞试验研究了有间隙的折叠操纵面(导弹翼面)的颤振特性。陈立勇[7]等分析了不可逆操纵面低速颤振的影响因素；杨智春[8]等对近年来带集中非线性环节的机翼非线性气动弹性问题进行了综述，总结和讨论了相关建模及分析方法；陈识等[9]针对含操纵系统间隙的飞机操纵面进行了振动及颤振特性研究，采用等效线性化方法进行了仿真并分析了其颤振特性。李宇飞等[10]研究了间隙非线性对二元翼面颤振特性的影响，利用描述函数法表征非线性间隙，计算了带间隙非线性的三自由度二元翼面极限环振荡现象。邰瑞雪等[11]针对飞机副翼间隙控制进行研究，并从制造工艺角度研究优化改进方案以减小副翼间隙。

现有研究多针对二元翼面开展间隙非线性对翼面颤振特性的影响，少有针对三维复杂机翼的研究，并且操纵面间隙对颤振极限环振荡触发速度和振荡频率的影响机理尚未完全研究清楚。因此，本文针对三维复杂机翼-副翼系统开展间隙非线性对其颤振特性影响的研究，并探索间隙非线性影响下LCO触发速度、振荡频率、和幅值间的关系。

1 考虑间隙非线性的三维操纵面颤振模型

为研究操纵面颤振特性，将机翼-副翼组合体拆分为机翼子结构及副翼子结构两部分，如图1所示。其中，机翼子结构包含发动机的固执结构；副翼子结构为自由舵面，具有6个刚体模态。各子结构节点位移可写为式(1)形式：

(a) 机翼子结构

(b) 副翼子结构图1 机翼及副翼子结构Fig.1 Substructures of wing and aileron

u=Φrqr+Φkqk+Φdqd

(1)

式中:u为节点位移；Φr，Φk，Φd分别为对质量节点矩阵归一化的刚体模态、低阶弹性模态、高阶弹性模态；qr，qk，qd分别为刚体模态坐标、低阶弹性模态坐标、高阶弹性模态坐标。

定义fJ为子结构对接自由度上的对接内力，Ψd为子结构剩余模态。

为处理高阶模态，引入双协调动态子结构法[12]。该方法基于瑞莱-李兹原理，用部件关于其对接边界反力的剩余模态集来近似替代截尾主模态[13]。由双协调子结构法可得式(2)近似关系(详细力学原理及论述过程可参见文献[13] ，本文中式(2)为文献中式(15)的变体)

Φdqd=ΨdfJ

(2)

式中：Ψd由子结构有限元模型及约束求得。

对每一个子结构，气动弹性动力学方程为

(3)

式中:Mg,Kg分别为质量矩阵及刚度矩阵；faero为作用于节点上的气动力。将式(1)代入式(3)可得动力学方程如式(4)

(4)

(5)

(6)

对于机翼及舵面子结构的对接节点，其对应点位移相等、作用力相反的双协调条件可表达为

(7)

当需考虑对接间隙时，引入对接点之间的位移差δ，带间隙的位移协调条件改写为

uJg,wing=uJg,aileron+δ

(8)

考虑对接点位移、内力协调条件，则可由式(4)得到消去节点内力的装配方程

(9)

式中:T为装配矩阵；fNL(δ)为间隙处节点的非线性内力。特别地，机翼-副翼系统仅有一个间隙环节，可定义标量形式的为副翼偏转角度δ及舵面偏转恢复力矩fNL(δ)，其关系由式表示

(10)

式中:d为间隙大小；kδ为舵偏刚度，对确定的机翼-副翼系统，可先进行估算，而后通过开展地面振动试验令计算模态频率与测试模态频率总误差最小的方式对其值进行修正。

2 带间隙非线性的操纵面颤振频域分析方法

2.1 非线性气动弹性方程

对具有单个舵偏间隙非线性环节的机翼-副翼系统，式可简化为如下形式

(11)

引入位移的简谐振荡假设

(12)

将式(12)代入式(11)，可得气动弹性方程的频域形式

(13)

式中:ω为系统颤振频率；Uq为模态坐标振幅；Uδ为舵偏振幅。式(13)右端非线性项使得V-g方法难以应用于该问题的颤振计算，因此引入间隙刚度描述函数对其进行线性化处理。

2.2 舵偏间隙刚度的描述函数

对fNL(δ)作傅里叶展开，可得到舵偏角频率ωδ倍数的简谐振荡叠加形式

(14)

各展开分量中，频率为ωδ的简谐振荡为fNL(δ)的主谐波响应。保留主谐波响应，可获得与时间无关的频域表达式

(15)

舵偏间隙刚度的描述函数:

(16)

(17)

(18)

则描述函数式(16)可简写为

(19)

引入描述函数后，式(13)可写为

(20)

式(20)具有线性气动弹性方程的形式，可利用V-g方法进行计算。

分别由机翼及副翼子结构的局部模态计算得到模态气动力影响系数矩阵Qq,wing和Qq,aileron，通过子结构装配矩阵T得到系统总体模态气动力影响系数矩阵

(21)

(22)

(23)

(24)

(25)

3 带间隙非线性的操纵面颤振时域分析法

式(13)中，气动力影响系数矩阵Q是在给定模态后，通过频域亚声速偶极子格网法(Doublet-Lattice Method, DLM)给出的，因此若要在时域内求解(13)式，首先要将Q矩阵拟合到时域中。时域和拉氏域有如下关系

(26)

式中,L为Laplace变换，s为Laplace算子。定义减缩频率k

k=ωL/V

(27)

式中:L为参考弦长;V为风速;ω为简谐振荡频率。气动力影响系数矩阵Q是k的函数，为使其先拟合至拉氏域再转换至时域，假设拟合后的气动力矩阵具有如下形式

(28)

式中:A0,A1,A2,D,E为系数矩阵，R为滞后根，参数p满足如下关系

p=ik=sL/V

(29)

选取一系列减缩频率k，使得

(30)

由最小状态(Minimum-State, MS)法进行时域气动力拟合，得到A0,A1,A2,D,E；随后将式(29)代入式(28)可得

(31)

设总模态位移为

(32)

则系统的时域气动弹性方程为

(33)

为消除上式中Laplace算子，引入增广状态

(34)

则有时域表达式

(35)

将式(34)、式(35)代入式(33)，可得机翼-副翼系统在时域下的状态空间形式气动弹性方程

(36)

式中,各系数矩阵表达式如下

(37)

式(36)、式(37)、式(10)即为带间隙的机翼-副翼系统时域仿真所需的方程。

4 带间隙机翼-副翼系统颤振的数值仿真

4.1 机翼-副翼系统结构模型及气动模型

以某型飞机机翼-副翼系统开展数值仿真，机翼系统由主翼面、发动机、副翼操纵面构成，按前文所述，为研究操纵面颤振特性，将系统分为机翼(包含发动机)及副翼两个子结构(见图1)。舵面偏转角度与舵偏恢复力矩间的间隙非线性关系如图2所示。计算气动力所用网格如图3所示。

图2 舵偏转恢复力矩的非线性函数Fig.2 Non-linear function of flap-restoring moment

图3 气动力计算网格Fig.3 CFD grid

副翼结构的刚度特性由舵偏片、舵弯片决定，如图4所示。舵偏片及舵弯片截面参数列如表1所示。

图4 副翼结构Fig.4 Structure of aileron

表1 刚度原件的截面参数Tab.1 Cross section parameters of stiffness units

为计算系统气动力影响系数矩阵Q，首先构造系统的子结构局部模态：构造机翼子结构局部模态，如图5所示。选取机翼前20 阶根部固支弹性模态和全部附着模态，在每个模态中加入位移强制为零的舵面节点。图5为机翼第4 阶一阶弯曲模态加入位移强制为零的舵面子结构节点。同理，选取舵面子结构的前10 阶自由模态和全部附着模态，并在每个模态中加入位移强制为零的机翼节点；舵面第7 阶一弯模态加入位移强制为零的机翼子结构节点，如图6所示。

图5 机翼子结构局部模态Fig.5 Local mode of wing substructure

图6 舵面子结构局部模态Fig.6 Local modes of rudder substructure

4.2 地面振动试验与结构动力学仿真结果

前面提到，舵偏刚度kδ一般按先估算，后修正的方式处理。根据振动试验，选取kδ使仿真计算的主要结构模态频率与地面试验测得的频率尽可能相近。数值仿真试验中，舵偏刚度kδ的初始估算值为1.419 N·m/rad；经反复迭代，舵偏刚度的修正值为1.547 N·m/rad。“4.1”节所定义的机翼-副翼系统在无间隙状态下的主要模态振型及频率，如表2所示。表中第2列、第3列分别为地面振动试验结果及根据修正后的kδ所得出的数值计算结果。结果显示，在修正的kδ值下，数值计算所得各阶振型及频率均与地面试验较好地吻合，证明了颤振模型在不考虑刚度间隙情况下的准确性。

表2 地面振动试验及数值仿真结果Tab.2 Results of round vibration test and numerical simulation

(38)

(a) 结构频率曲线

(b) 部分模态频率受间隙影响范围图7 间隙舵偏刚度下结构频率计算结果Fig.7 Calculation results of structural frequency with freeplay in deflection stiffness

计算结果显示，机翼一扭(第9阶)、舵面偏转模态(第11阶)和舵面弯曲模态(第15阶)以及它们的混合模态(第10阶、第17阶)在舵偏振幅略微大于间隙时出现了显著变化。同时，系统舵面偏转模态(第11阶)频率随间隙增大最低下降至线性系统的一扭模态(第9阶)后不再继续下降，但舵面弯曲模态(第15阶)频率会下降直至低于原第14阶模态(发动机模态)，从而使得这两阶频率顺序交换。

4.3 非线性颤振的频域仿真

图8 非线性颤振分析图Fig.8 Analysis diagram for Nonlinear flutter

图9 极限环风速-对数相对幅值频域计算图Fig.9 Wind speed-logarithm relative amplitude of limit cycle by frequency domain calculation

计算结果显示，随着风速增大，机翼-副翼系统经历了四个阶段。

第一阶段：Vf<21.18，系统无颤振解；

4.4 非线性颤振的时域仿真

取大气密度ρ=1.225 kg/m3，参考弦长L=0.25 m，滞后根R=[0.2,0.5,1.0,1.5]，间隙大小δ=0.01 rad初条件为模态坐标q|t=0=0，舵偏刚度kδ=1.547 N·m/rad，对机翼-副翼系统进行时域仿真。机翼-副翼系统在风速Vf从0逐渐增大的过程中经历了如下阶段变化，如图10所示。

图10 舵面偏转的时域响应及极限环响应Fig.10 Time domain response and limit cycle response of rudder deflection modal coordinates

第一阶段：Vf<20.75，系统无颤振；

第四阶段：Vf∈[32.00,32.50)，系统颤振幅值、频率快速变化；

提取出现稳定极限环时舵偏运动的振幅、频率和风速，得到LCO曲线，如图11所示。除舵面偏转模态外，系统其余各阶模态均有类似的结果。时域仿真发现，系统存在多个频率的稳定极限环，不同LCO之间具有较为固定且范围狭窄的转换风速。当系统风速处于后一个LCO风速范围内时，系统在时域响应中会先经历前面一个(或几个)LCO状态却不能保持稳定，然后过渡到最终的稳定LCO。

(a) 风速-幅值曲线

(b) 风速-频率曲线图11 极限环时域仿真计算结果Fig.11 Time domain simulation results of limit cycle

4.5 频域与时域计算结果对比

将LCO曲线的时域仿真结果和频域计算结果对比,如图12所示。

(a) 相对风速-对数相对舵偏幅值

(b) 风速-频率图12 频域与时域计算对比Fig.12 Comparison of frequency domain and time domain calculation

通过时域方法获得的三种频率的极限环分别与通过频域方法获得的三个颤振分枝对应，且其频率吻合度较高(见图12(a))。计算获得的极限环产生时的气流速度(时域法20.75 m/s，频域法21.18 m/s)、系统发散时的气流速度(时域法37.00 m/s，频域法38.01 m/s)亦吻合度较高。

两种算法计算得到的极限环转变模态间有一定差异，且差异随风速增大而愈发显著(见图12(b))。尤其在大风速下，时域方法计算得到的转变速度及幅值都明显小于频域方法的结果。

分析其原因，在使用频域法时引入的描述函数式(16)忽略了非线性刚度环节，因此频率响应的高频部分有所缺失，从而造成响应幅值偏低。在Vf=35 m/s，系统响应为剧烈振荡过渡阶段，被忽略的高频谐波成分占据能量比例更大，可能造成更大的计算误差。

5 间隙对机翼-副翼系统振荡特性影响

为研究间隙对机翼-副翼系统振荡特性的影响，进行了低速风洞试验。试验对象为三维矩形机翼-副翼系统，操纵面转动自由度带有间隙，试验安装如图13所示。舵面偏转间隙产生装置(联轴器)，通过改变联轴器间隙槽的大小模拟不同的间隙角度，如图14所示。试验中分别于副翼后缘、机翼梢部后缘、机翼梢部主梁布置加速度传感器，于机翼内部布置角位移计。试验风速范围为10～45 m/s。试验中采用4种间隙大小(包含无间隙情况)，不同间隙工况下颤振或极限环发散速度，如表3所示。舵面加速度响应及舵偏角位移响应幅值随风速变化的规律，如图15～图16所示。

图13 矩形机翼风洞试验Fig.13 Wind tunnel test of rectangular wing

图14 间隙产生装置Fig.14 Freeplay generator

图15 舵面加速度响应幅值Fig.15 Amplitude of aileron acceleration

图16 舵偏角响应幅值Fig.16 Amplitude of aileronrotating angle

表3 不同间隙条件下的LCO发生风速及发散风速Tab.3 Wind velocity of LCO occurring and diverge under different freeplay

试验结果显示，线性系统(间隙为0°)颤振速度大于45.5 m/s，颤振频率约13.5 Hz。间隙的存在使系统接近发散的风速显著降低，且间隙越大，颤振临界风速越低。

在试验间隙工况下，LCO出现速度均在33.0～33.5 m/s范围内，间隙大小的变化对LCO出现速度影响不显著；LCO幅值远大于非LCO幅值，幅值随风速增大而增大，且间隙越大，LCO幅值随风速变化的斜率越大。

6 结 论

本文针对带间隙非线性的机翼-副翼系统的颤振特性开展了理论与试验研究。利用双协调动态子结构法对间隙影响的缩比机翼-副翼结构分别进行建模，建立了三维机翼带间隙操纵面的非线性气动弹性方程。利用地面试验验证了不考虑间隙非线性时计算方法的准确性；进一步分别在频域内用谐波平衡法、在时域内利用有理函数拟合和时域推进方法进行了仿真计算，得到间隙非线性极限环解。两种方法在LCO生成气流速度及系统发散气流速度的计算结果有一定吻合度，在LCO模态转变的规律上显示出一致性，互相验证了算法的有效性。同时，两种方法的计算结果也存在差异，从差异产生的机理上看，主要源于对非线性的舵面偏转恢复力局的处理上：在频域法中为使系统方程具备线性特征而引入了描述函数，将非线性恢复力进行傅里叶展开并略去了非线性刚度环节频域响应的高频成分，这一处理使计算得到的LCO转换模态响应幅值偏低、在非LCO的剧烈振荡过渡阶段引入较大计算误差。最后，通过风洞试验研究了三维矩形机翼操纵面偏转间隙对系统振荡特性的影响，结果显示间隙使系统出现明显的LCO，在LCO转换阶段振荡幅值出现剧烈变化；并且间隙的存在使系统振荡发散风速显著降低，该结果说明间隙的存在增加了系统破坏的风险。