大展弦比柔性机翼非线性颤振研究

2019-07-05张忠源段静波

兵器装备工程学报 2019年6期

张忠源,段静波,路平,叶廷

(1.陆军工程大学石家庄校区无人机工程系，石家庄 050003； 2.石家庄铁道大学工程力学系，石家庄 050003；3.中国人民解放军71939部队，济南 250000)

当机翼攻角大于失速攻角时，机翼升力会迅速减小，对飞行器的安全造成极大威胁[1]。所以实际情况中，机翼所受气动力是非线性的，尤其对于大柔性机翼，受气动力作用，机翼变形更加明显[2]，极大提高了机翼产生失速的可能性。为解决这一问题，就要在计算大柔性机翼颤振时引入非线性气动力模型，本文将ONERA非线性气动力模型引入大展弦比柔性机翼颤振问题，同时基于半解析半数值的传递函数方法[3]，该方法求解过程简洁和统一，边界条件处理规范和方便，可以更加快捷、有效、准确地计算出机翼的颤振特性。

由于大柔性机翼扭转会产生非线性气动力，CC Xie等[4]针对这个问题进行了研究，用平面双点阵方法计算频域内的非定常气动力，忽略偏转翼的弯曲效应。在给定的载荷条件下，对系统进行气弹性稳定性分析。同样针对大展弦比大柔性机翼，密歇根大学C Cesnik[5]团队也进行了深入研究，并搭建了收集几何非线性气动弹性响应的数据实验平台，为飞机提供可在飞行中测量的特定气动弹性特征，例如，耦合的刚性、弹性体不稳定性，阵风期间的大的机翼偏转等。D Tang[6]将柔性机翼的气动弹性分析与风洞试验相结合，介绍了一种弹性载荷作用下柔性大展弦比翼型气动弹性模型的理论气动弹性模型。MJ Patil等[7,8]提出使用完整飞机模型的气动弹性特性来获得结果，由于机翼具有大柔性，飞机整体的飞行动态特性也会发生变化，并用严格的非线性气动弹性来分析。

本研究对大展弦比柔性机翼进行非线性颤振分析，求解思路如下：先将柔性机翼视为曲梁，结合非线性ONERA气动力模型得到机翼颤振微分方程，为传递函数法求解做好准备。然后，将颤振微分方程转换为状态空间方程形式，通过求解复特征值问题，获得了大展弦比柔性机翼在非线性气动力作用下的颤振速度和颤振频率。

1 柔性机翼非线性颤振分析模型

1.1 曲梁单元振动微分方程

将柔性机翼视为曲梁，根据文献[9]，考虑曲梁质心轴与弹性轴不重合的情形，可获得曲梁六自由度振动方程为：

(1)

如图1所示，曲梁曲率为常数R，u,w,v分别为曲梁单元沿坐标轴ξ,η,ζ三个方向的位移，ψζ,ψη,φξ分别为曲梁单元绕ζ,η,ξ轴的扭转角。ρ为曲梁的密度，A为曲梁单元截面积。E,G为拉伸弹性模量和剪切弹性模量，k为剪切刚度修正系数。ζa为弹性轴到质心轴的距离。Iζ,Iη,J分别为绕坐标轴ζ,η,ξ三个方向的惯性矩。t为时间。qξ(ξ,t)、qη(ξ,t)、qζ(ξ,t)为沿坐标轴ξ,η,ζ三个方向的分布力，mξ(ξ,t)、mη(ξ,t)、mζ(ξ,t)为绕坐标轴ξ,η,ζ三个方向的分布弯矩。

1.2 ONEAR非线性气动力模型

在忽略机翼重力影响条件下，机翼颤振时的外力为气动力产生的分布升力以及分布扭矩。本文采用片条理论进行非定常气动力计算。对于大展弦比大柔性机翼，在正常飞行期间局部迎角可能非常大，这将导致机翼失速，即强烈的非线性不稳定失速现象。在本章中，使用了ONEAR[10]动态失速模型，并且给出了作用在四分之一弦长处的升力和力矩的公式：

(2)

式中：L表示机翼升力；V是来流速度；α为瞬时攻角；φξ为有效攻角，相当于机翼实际扭转角；b为机翼半弦长；Cz表示空气动力系数；Cza表示线性气动力部分对应的气动力函数；Czb表示非线性气动力部分对应的气动力系数。参数下标z=L时，表示与升力有关的系数；当下标z=M时，表示与力矩有关的系数。由于ONEAR线性部分与经典Theodorson理论一致，所以其线性部分参数选取如表1所示。

表1 ONERA模型线性部分参数值

非线性气动力部分的有关参数有r1z，r2z，r3z。

为简化计算，通常用折线来近似代替静态空气动力曲线，如图2所示。

ΔCz可以表示为：

(3)

(4)

这里w1/4表示四分之一弦长处的位移，根据几何关系，w1/4可以表示为：

(5)

将式(5)代入式(4)并整理得：

(6)

将式(6)代入式(2)，可得：

z=LorM

(7)

单位展长上的升力和相应的俯仰力矩可以表示为：

(8)

1.3 机翼单元颤振微分方程

将式(8)代入式(1)，即可得到大变形柔性机翼的颤振微分方程。

(9)

2 传递函数法求解

对机翼的颤振微分方程(9)进行Fourier变换，并整理得：

(10)

其中，气动力系数CL和CM的Fourier变换可以表示为：

(11)

根据传递函数方法，定义状态向量：

将式(10)改写为状态空间形式方程：

(12)

其中，ge(ξ,ω)=0，转移矩阵Fe(ω,V)为12×12的方阵，其非零元素为：

Fe(8,10)=-1,Fe(9,10)=1

(13)

边界条件为：

Mbeηe(ξ=0,ω)+Nbeηe(ξ=1,ω)=γe(ω)

(14)

其中，γe(ω)为由位移和力组成的列向量，Mbe、Nbe分别为单元边界条件选择矩阵，γe(ω)、Mbe、Nbe表达式为：

式(12)的传递函数解为：

ηe(ξ,ω)=He(ξ,ω,V)γe(ω)

(15)

其中，

He(ξ,ω,V)=eFe(ω,V)ξ[Mb+NbeFe(ω,V)]-1

(16)

由于大变形后的机翼划分为若干常曲率曲梁单元来描述整个机翼，其求解方程可借鉴有限元方法的思想进行组集。曲梁单元截面上的内力为：

(17)

其中：Mη、Mζ分别为曲梁单元绕η、ζ轴弯矩；Tξ分别为曲梁单元绕ξ轴的扭矩；Nξ为沿ξ轴的轴力；Qη和Qζ分别为沿η、ζ轴的剪力。

将曲梁单元内力写成矩阵形式：

Qe(ξ)=Qeη(ξ)ηe(ξ,ω)

(18)

(19)

将式(15)代入式(18)，可的达到曲梁单元端点处的内力：

(20)

按照有限元组集方法对各节点统一编号，并进行组集拼接，可得到机翼整体平衡方程为：

K(ω,V)γ(ω)=f

(21)

其中，

(22)

其中：K(ω,V)可视为整体刚度矩阵；γ(ω)可视为整体节点位移向量；f为各单元节点内力拼装成的向量。由于本文中将机翼的气动力与机翼作为一个完整的系统来考虑，除此之外，机翼没有受到其他外力作用，因而根据单元节点内力与外载荷平衡，可得出

f=0

(23)

当机翼颤振时，γ(ω)应有非零解，此时须满足条件

det[K(ω,V)]=0

(24)

由于K(ω,V)为复矩阵，其行列式值等于零的必要条件为矩阵行列式值的实部与虚部均为零，即

(25)

上式两个方程包含两个未知变量V和ω，因而求解式(25)可得到V和ω的解。其中，V即为机翼颤振速度，ω即为机翼颤振频率。通常满足式(25)条件的解可能有多个，其中V最小的一组解即为机翼的颤振特性。

3 算例分析与讨论

前两节推导出了非线性气动力影响下的机翼颤振求解方法，为研究非线性气动力与线性气动力对机翼颤振造成的差异，选取某大柔性机翼模型进行非线性气动力和线性气动力的颤振计算，机翼模型参数如表2所示，通过改变机翼线密度，展弦比，抗弯刚度和抗扭刚度等参数，对比线性与非线性气动力模型下的颤振结果的差异。

表2 大柔性机翼模型参数

这里取3种典型变形进行分析，记变形一、变形二和变形三，对应的翼尖位移取半展长的3.125%、6.250%、12.50%，分别对3种变形使用线性和非线性颤振分析。计算得到3种变形下的线性和非线性的机翼颤振速度和颤振频率，如图3和图4所示。

使用表2中的参数进行非线性气动力条件下的颤振计算，通过比较线性气动力和非线性气动力求解的颤振计算结果，可以发现，运用非线性气动力求解的颤振速度要低于线性气动力求解的颤振速度，是因为随着弹性变形的增加，机翼的气动力系数减小，其相应的气动特性也降低，随之颤振速度也降低。颤振频率通过对比发现，随着变形程度的加大，变化趋势相同，但是非线性气动力计算得到的颤振频率变化更加剧烈。

将气动弹性系统中的弯曲刚度、扭转刚度、弹性轴位置、半展现比和机翼线密度等参数作为变量，以表2中的参数为基准，以变形条件6.25%为例，对比线性气动力和非线性气动力求解的颤振计算结果的变化规律。定义垂直弯曲刚度与表2中垂直弯曲刚度值的比为弯曲刚度比：ηb=EIb/EIζ，扭转刚度与表2中扭转刚度值的比：ηt=GJt/GJ，改变垂直弯曲刚度和扭转刚度，得到的线性与非线性颤振速度变化趋势如图5、图6所示。随着垂直弯曲刚度和扭转刚度的增加，颤振速度也随之提高，增加垂直弯曲刚度和扭转刚度可以有效改善其气动弹性稳定性，可见在进行大柔性无人机设计时需要着重注意弯曲和扭转刚度的设计。

密度和展弦比是大展弦比机翼的重要设计参数(见图7、图8)，以表2中的参数为基准，以变形条件6.25%为例，机翼颤振速度随机翼半展现比的变化比较简单，非线性和线性的结果都是单调递减的，也就是说展现比越大，柔性机翼的气动弹性稳定性越差。颤振速度随机翼密度的增加而减小，但是线性和非线性的计算误差随着密度增加而增大，这是由于随着机翼的质量增加，所受气动力增加，导致两种气动力的计算结果误差变大。

最后来研究弹性轴在机翼中的相对位置对大柔性无人机气动弹性稳定的影响，大柔性无人机的弹性轴在弦向位置是非常重要的气动弹性设计量，如图9所示，横轴表示弹性轴在弦长相对位置的百分比，纵轴为颤振速度，随着弹性轴后移，颤振速度降低较快，即弹性轴后移严重削弱了大柔性无人机的气动弹性稳定性，非线性颤振速度相比于线性计算结果降低的比较慢。