直升机悬停状态下的振动计算

姚 颂,芮筱亭,王景弘

(南京理工大学 能源与动力工程学院，南京 210094)

振动是影响直升机发展的主要问题之一.对于大多数现代直升机而言，机动能力和最大前进速度受到过度振动的限制.发动机、变速箱、尾桨和旋翼质量不平衡都是振动激励源.然而直升机振动的主要来源是主旋翼.振动对直升机操作的几乎所有方面都是有害的，包括飞行员疲劳和结构完整性；降低电子设备的可靠性和操作准确度；影响相机和测量设备等机器的精度.旋翼桨叶可以模拟为旋转梁[1]，由于弹性变形和刚体运动的耦合，旋转梁的振动特性与非旋转梁的振动特性有很大差异[2].沿着梁长度变化的离心力改变了梁的整体刚度(应力刚化效应)，使固有频率和模态发生变化[3].众多研究人员针对不同的课题，采用不同的方法研究了旋转梁的动态特性.例如Southwell原理、积分矩阵法、Rayleigh Ritz法[4]、Galerkin法[5]和扰动法，以及动态刚度法，这涉及到幂级数解及使用Frobenius方法计算相应微分方程的解[6].传统的有限元法(Conventional finite element method,CFEM)[7]已被用于许多分析旋翼固有频率的问题.一种基于低自由度模型的谱有限元法(Spectral finite element method，SFEM)[8]被用于分析旋转锥形梁的运动特性.传递矩阵法用于非旋转锥形梁[9]和带裂纹的旋转梁[10]的自由振动计算.

有限元法和多体系统动力学是武器系统动力学的重要基础.然而，使用这些方法对多刚柔体耦合的复杂武器系统的振动特性计算可能带来很大的计算量和由计算不良条件引起的计算奇点.多体系统传递矩阵方法(MSTMM)[11-13]可以避免复杂线性多系统特征值问题的计算错误，因为系统总矩阵阶次较低，显著提高了振动特性的计算速度[14].本文基于MSTMM，建立了耦合直升机四叶片旋翼/机身的多体系统动力学模型.推导出整体传递方程，并模拟了旋翼振动特性.为了证明所提出方案的准确性，与商业软件ANSYS Workbench仿真结果进行比较.

1 建模与传递矩阵推导

1.1 直升机拓扑模型

图1中所示的典型直升机可以建模为多刚性-柔性耦合系统，直升机主要分为旋翼系统、机身和机尾3个部分.

图1 典型直升机模型

拓扑结构如图2所示，四桨叶/机身耦合系统的动力学模型在空间内振动.全局惯性系oxyz用于描述系统运动.本文将直升机模型简化为8个元件，分别是4个旋转Euler-Bernoulli梁组成的无铰接柔性桨叶、桨毂、传动轴、机身和尾梁8个部分.元件编号1至4号是桨叶，考虑其在挥舞方向、摆振方向及扭转方向的振动，桨叶通过柔性连接到5号元件桨毂，桨毂处理为一个空间刚体，桨毂连接至元件6传动轴，传动轴处理为一根围绕中心轴旋转的梁，传动轴与元件7机身相连，机身处理为一个刚体，机身与元件8尾梁为刚性连接.旋翼与传动轴以同一转速旋转，旋转角速度为Ω.

图2 直升机拓扑模型

1.2 系统总传递矩阵的推导

递矩阵法使用状态矢量描述力学系统中任一点的力学状态，坐标系方向及符号约定如图3所示，状态矢量是由两部分元素构成的一个列阵，一是该点状态变量的位移及角位移(x,y,z,θx,θy,θz)，二是内力及内力矩(qx,qy,qz,mx,my,mz).多自由度无阻尼振动系统可由模态振动叠加而成，某阶模态下固有频率ωk对应的状态矢量可以表示为z=Zke-iωkt,Zk为模态坐标下的状态矢量，其状态变量为固有频率ωk对应的位移及内力的模态坐标，为了书写简洁将上标k省略，即Zi,j=[X,Y,Z,Θx,Θy,Θz,Mx,My,Mz,Qx,Qy,Qz]T，其中i为体元件编号，j为铰元件编号.

图3 坐标系及符号约定

系统总传递方程由总传递方程和几何方程构成，总传递方程的状态矢量由所有边界点的状态矢量构成.主传递方程中的每一个系数矩阵为从相应边界点到系统根的传递路径上所有元件传递矩阵的连乘积.几何方程用来描述同一元件不同输入点之间的几何关系，其系数矩阵可以在推导主传递方程的过程中依据同一元件不同输入点之间的几何关系自动导出.对于每一个分支，都可以根据同一个体元件不同输入点的几何关系推导出一组几何方程，几何方程的个数等于系统树梢个数减1.几何方程的系数矩阵为由从边界端到分叉体元件的输入端Ik(k=1,2,…，L)路径上各个元件的传递矩阵的连乘积，再左乘(-Hj,1)(若k=1)，或左乘(Hj,k)(若k=2,…，L).

在线性系统中ZO=UjZI，其中：ZI为输入点的状态矢量，ZO为输出点状态矢量，Uj为元件j的传递矩阵，由牛顿力学定律推导得到，总传递矩阵Uall为沿传递路径各元件传递矩阵依次连乘之积.最终整理可以得到UallZall=0，其中Zall为各边界点的状态矢量和，Uall即为总传递矩阵.

如图2拓扑图所示，该模型共有4个边界点，分别为Z1,0,Z2,0,Z3,0,Z4,0,Z8,0,以Z8,0为输出点，系统的传递矩阵推导如下：

Z8,0=U8Z7,8=U8U7Z7,6=U8U7U6Z5,6=

U8U7U6(U5,I1Z5,1+U5,I2Z5,2+

U5,I3Z5,3+U5,I4Z5,4)=U8U7U6(U5,I1U1Z1,0+

U5,I2U2Z2,0+U5,I3U3Z3,0+U5,I4U4Z4,0)=

U8U7U6U5,I1U1Z1,0+U8U7U6U5,I2U2Z2,0+

U8U7U6U5,I3U3Z3,0+U8U7U6U5,I4U4Z4,0,

(1)

或

T8-1Z1,0+T8-2Z2,0+T8-3Z3,0+T8-4Z4,0-Z8,0=0,

(2)

其中：

T8-1=U8U7U6U5,I1U1,T8-2=U8U7U6U5,I2U2,

T8-3=U8U7U6U5,I3U3,T8-4=U8U7U6U5,I4U4.

模型中有一个分叉点，4个分支，所以共有3个几何方程.第1个几何方程如下：

-H5,1U1Z1,0+H5,2U2Z2,0=0,

(3)

或

G5-1Z1,0+G5-2Z2,0=0,

(4)

其中：

G5-1=-H5,1U1Z1,0，G5-2=H5,2U2Z2,0.

第2个几何方程如下：

-H5,1U1Z1,0+H5,3U3Z3,0=0,

(5)

或

G5-1Z1,0+G5-3Z3,0=0,

(6)

其中：

G5-1=-H5,1U1Z1,0,G5-3=H5,3U3Z3,0.

第3个几何方程如下：

-H5,1U1Z1,0+H5,4U4Z4,0=0,

(7)

或

G5-1Z1,0+G5-4Z4,0=0,

(8)

其中：

G5-1=-H5,1U1Z1,0,G5-4=H5,4U4Z4,0.

联立式(2)、(4)、(6)、(8)，可得

Uall|30×60Zall|60×1=0,

(9)

其中：

2 主要元件传递矩阵的推导

2.1 旋转梁的传递矩阵的推导

图4 旋转梁挥舞方向振动的符号约定和坐标系

(10)

式中rH为梁的输入端距离旋转轴的距离，其中离心力F(x)可以表示为

(11)

(12)

梁的动能T由下式给出：

(13)

Hamilton原理可以表示如下:

(14)

将式(10)、(13)代入式(14)中，可以得到挥舞方向振动的运动微分方程：

(15)

将式(15)中位移改写为谱坐标形式，即z(x,t)=Z(x)eiωt，并代入式(11)、(12)，化简方程可得

ZIV(x)-(0.5C1+rHC1x+C2)Z″(x)+
C1(rH+x)Z′(x)+C3Z(x)=0,

(16)

其中:

方程(16)的根是使用幂级数中的Frobenius方法计算的，并且可以假设该微分方程的解如下：

(17)

将式(17)代入式(16)中，可以得到:

(18)

上式也可以写为如下形式：

C2(k+i+1)(k+i+2)ai+3+C1rH(k+i+1)2ai+2+

[C1(k+i)+0.5C1(k+i-1)(k+i)+C3]ai+1}xk+i.

(19)

代入i=-4到式(19)中，可以发现只有a1得系数非零，其他的项都等于0.因此，将i=-4代入式(19)中，发现a1是最低阶非消失系数，得到k(k-1)(k-2)(k-3)a1=0,其中a1≠0, 那么k只有0，1，2和3这4个根满足方程k(k-1)(k-2)(k-3)=0，并且a1是任意数，这里假设a1=1, 式(19)中的未知系数的一般形式可以用以下表达:

(20)

分别将i=-3,i=-2,i=-1及i=0代入式(20)中，依次得到:

(21)

从式(19)、(20)可以确定所有未知系数，对于k=0,1,2,3,微分方程式(19)的解可以用4个独立线性解的总和乘以任意常数表示如下:

Z(x)=A1f(x,0)+A2f(x,1)+A3f(x,2)+A4f(x,3),

(22)

式中A=[A1,A2,A3,A4]T是任意常数，方程(22)中k=0,1,2,3的4个独立解f(x,k)可以用函数的形式表示如下:

f(x,0)=a1(0)x0+a3(0)x2+a4(0)x3+a5(0)x4+…,

f(x,1)=a1(1)x1+a3(1)x3+a4(1)x4+a5(1)x5+…,

f(x,2)=a1(2)x2+a3(2)x3+a4(2)x4+a5(2)x5+…,

f(x,3)=a1(3)x3+a3(3)x4+a4(3)x5+a5(3)x6+….

(23)

由材料力学可得旋转梁在此平面内的弯矩和剪切力，分别为:

(24)

2rHx-x2)Z′(x)-F0Z′(x).

(25)

在MSTMM中，桨叶挥舞方向z-x平面内的边界条件为[Z,Θy,My,Qz]T，根据线性关系可以得到Θy=Z′，联立式(24)、(25)，可以写成如下形式Z(x)=D(x)A.

(26)

其中:

H1=EI[f‴(x,0)+γf′(x,0)],

H2=EI[f‴(x,1)+γf′(x,1)],

H3=EI[f‴(x,2)+γf′(x,2)],

H4=EI[f‴(x,3)+γf′(x,3)],

γ=(0.5C1x2+rHC1x+C2).

由式(26)可以得到在x=l1作为输入端的状态适量为ZI=[D(l1)]A，从而这组未知的系数向量可以表示为A=[D(l1)]-1ZI，那么位于x=l2的输出端的状态矢量为

ZO=[D(l2)]A=[D(l2)][D(l1)]-1ZI=URBZI,

(27)

式中URB=[D(l2)][D(l1)]-1为旋转梁在挥舞平面内横向振动的传递矩阵.

同理，在摆振运动平面内也以此方法获得其传递矩阵，但在离心力部分略有不同，由于离心力的作用方向分为摆振方向和展向分量，因此运动微分方程为

(28)

非圆截面梁的扭转通过材料力学可以得知微段的扭转角φ与扭矩T的关系为

(29)

式中:G为剪切模量；h为截面长边；b为截面短边；β为矩形截面杆扭转时的因数，当h/b远大于10时，β近似等于0.333,It=βhb3.在扭转方向上由于螺桨效应的影响，即由离心力Ω引起回复力矩，扭转振动的运动微分方程为

(30)

其中Iα为微段绕扭转轴的转动惯量.

通过求解式(30)这个微分方程，最终可以得到长为x的旋转梁扭转振动的传递矩阵如下：

(31)

2.2 传动轴的传递矩阵的推导

图5 旋转轴模型

(32)

代入模态坐标u(z,t)=U(z)eiωt和v(z,t)=V(z)eiωt，并引入量纲一的长度ξ=z/L与微分算子D=d/dξ，整理后式(33)变为

(33)

对式(33)求解后分别得到U(ξ)和V(ξ)的通解:

(34)

A1-A8与B1-B8是两组独立的8个常系数，其之间的关系如下：

(A1,A2,A3,A4)=kα(B1,B2,B3,B4),
(A5,A6,A7,A8)=kβ(B5,B6,B7,B8),

(35)

其中：

由式(32)可以得出扭转角，弯矩与剪切力的表达式如下：

(36)

(37)

(38)

将式(34)、(35)分别代入到式(36)～式(38)中，整理后可以得到x=l1处的状态矢量为

(39)

式中ZI=[U,V,θx,θy,Mx,My,Sx,Sy]T，以x=l2为输出端，可以得到：

(40)

(41)

式中US为旋转轴的传递矩阵.

3 计算与仿真对比验证

3.1 旋转梁的固有频率验证

为验证MSTMM计算旋转梁频率的准确性，分别对转速为0，20，40 rad/s时悬臂梁的频率计算，使用MATLAB编程计算与Ansys Workbench仿真结果的前五阶频率对比.梁的参数分别为，长度L=6 m，横截面长b=0.5 m，高h=0.04 m，密度ρ=334 kg/m3，弹性模量E=7.1 MPa，剪切模量G=2.730 8 MPa，泊松比μ=0.3. ANSYS Workbench的模型有88 413个节点，网格单元数为15 600.

计算结果与对比见表1，3种转速下前五阶频率结果误差较小.取3次相同计算条件下的平均时间，使用MATLAB编程的计算时间为6.3 s，ANSYS Workbench的计算时间是55.9 s.计算速度提高了8.8倍.

表1 悬臂梁在3种转速下MSTMM计算与ANSYS仿真的固有频率对比

续表1

3.2 旋转轴的频率计算结果验证

(42)

表2 圆截面悬臂旋转轴的固有频率(EIxx=EIyy)

3.3 悬停直升机整机频率计算

为了验证MSTMM方法计算的准确性，对图2直升机模型进行验证计算.直升机各部件的参数如下，桨叶密度ρ=334 kg/m3，长度L=6 m，横截面长b=0.5 m，高h=0.04 m，弹性模量E=7.1 MPa，剪切模量G=2.730 8 MPa；桨毂密度ρ=7 850 kg/m3，尺寸为长和宽为0.5 m，高0.3 m；传动轴密度ρ=7 850 kg/m3，尺寸为直径0.16 m,，高为1 m，弹性模量E=2.5 MPa，剪切模量G=9.615 4 MPa；机身密度ρ=170 kg/m3，长3 m，宽1.8 m，高1.6 m，质心与机身传动轴连接点在X轴上的投影距离0.4 m；尾梁密度ρ=140 kg/m3，截面高和宽为0.8 m，长为4 m，弹性模量E=2 MPa，剪切模量G=7.692 3 MPa.尾梁末端采用固定约束，旋翼系统的转速为36.651 9 rad/s. ANSYS Workbench模型的节点数88 260，16 056个单元. 前12阶频率见表3.

表3 使用MSTMM方法计算旋转状态直升机固有频率与ANSYS仿真结果对比

计算结果与仿真结果基本一致.3次相同计算条件下的平均时间，使用MATLAB编程的计算时间为11.5 s，ANSYS Workbench的计算时间是81.2 s.计算速度提高了7.1倍.由此可以使用MSTMM方法计算无约束条件即悬停状态下直升机模型的固有频率，前8阶的计算结果见表4.

表4 悬停状态下直升机模型固有频率

4 结 论

1)本文推导了旋转梁在挥舞、摆振和扭转3个自由度的运动方程，得到其传递矩阵，算例通过与ANSYS仿真结果对比，误差在1.5%以内，验证了结果的准确性.

2)推导了旋转轴的传递矩阵，算例与参考文献结果对比，误差不超过0.15%，验证了准确性.

3)MSTMM可以快速准确地计算动力学问题.并且相较于传统动力学软件Workbench ANSYS，结果相近，计算时间大大缩短，本文算例计算时间缩短了7倍以上.

4)对直升机旋翼/机身耦合系统验证计算，与仿真结果基本一致，并得到了悬停状态下的前8阶固有频率，对直升机设计具有一定参考价值.