旋翼桨叶结构载荷计算方法比较研究

吴 杰, 杨卫东, 虞志浩

(南京航空航天大学 直升机旋翼动力学重点实验室，南京 210016)

直升机旋翼动力学的挑战之一是精确地预测桨叶结构振动载荷。除了结构动力学建模与气动力建模精度对载荷有重要影响外，载荷计算方法也同样是非常重要的因素。Bielawa[1]最早采用力积分法与模态叠加法两种载荷计算方法对无铰式旋翼开始了这方面的研究。文章指出力积分法虽然能得到更准确的结果，并且能以更少的模态收敛，但其实现过程更为复杂。Thomas[2]利用CAMRAD对几种直升机旋翼桨叶(BO105，CH-34及SA349/2)进行了力积分法与曲率法的对比研究。他指出当剖面之间的桨叶结构特性相差不大时，曲率法能获得较好的载荷；而力积分法则依赖积分步长以及桨叶分段数的选择，并且弯矩预测比剪力预测需要更多的分段数和更小的积分步长。上世纪八十年代，反力法在著名的旋翼动力学综合分析平台2GCHAS中被引入。Lim等[3]利用2GCHAS以及CAMRAD/JA研究了UH-60A直升机结构载荷数据。文章认为曲率法的缺点在于高阶导数引起的数值精度损失；同时力积分法很难处理近桨根处的载荷计算，因为桨根处剖面特性变化通常较为剧烈；而反力法采用了有限元方程中组集之前的单元矩阵计算节点载荷，预测精度较高。

本文着重比较上述三种载荷计算方法对于桨叶结构载荷的影响。采用刚柔耦合模型[4]描述旋翼桨叶的动力学运动关系。该模型将铰接式旋翼铰与轴承处的转角抽象成三个独立自由度，作为刚体转角与桨叶弹性变形耦合，相较于采用小转角假设的传统有限元模型在瞬态响应计算及碰撞研究中具有明显的优势[5]。气动力由准定常气动模型计算得到，并计入所有刚体转角及弹性变形的影响。由于尾迹模型对于挥舞弯矩的预测至关重要[6]，本文采用自由尾迹入流模型计算桨盘入流。剖面气动力系数从翼型数据表中根据有效迎角与来流速度插值得到。离散后的气动力以广义力的形式表达，运动方程依据Hamilton原理推导建立。由于传统积分方法对于瞬态响应积分收敛较慢以及力积分法会引入数值累积误差，本文发展了一种新的基于Richardson外插技术的自适应变步长Legendre-Gauss数值积分算法，加速收敛速度，降低累积误差。

1 运动方程及求解方法

1.1 坐标系

(1)

其中，u,v,w为桨叶在变形前坐标系中度量的弹性变形，η,ζ为质点位于剖面坐标系中的坐标。矢量在不同坐标系中的坐标可经由转角余弦矩阵相互转换得到。比如变形前坐标系中的矢量坐标经过转换矩阵TDU转到变形后坐标系中，即：rD=TDUrU，其中TDU是关于Euler角的函数矩阵。

图1 铰与轴承的运动

为了充分考虑铰接式旋翼桨叶在铰及变距轴承处的惯性运动，引入三个刚体自由度，用以模拟铰与轴承处的桨叶转角，分别是刚体挥舞角(θ1+βp)、刚体摆振角(θ2)以及刚体变距角(θ3+θ0)，如图1所示。βp指桨叶预锥角，θ0为操纵变距。θ1、θ2与θ3分别是因惯性力与外力综合作用引起的挥舞铰、摆振铰以及变距轴承的刚体转角。为方便列写动能变分式，式(1)中的位置矢量通常需要经余弦转换矩阵转换到惯性系中。桨叶动能变分项最终由剖面动能变分沿桨叶展向积分得到：

(2)

势能变分项由中等变形梁理论[7]推导得到。桨叶抽象为细长柔性Bernoulli梁，变形过程中剖面与桨叶参考轴线始终保持垂直，不计横向剪切。为了更精确地对弹性扭转建模，本文在势能变分项推导中未使用阶次准则。文献[8]指出，保留弹性扭转的高阶项(三阶或以上)对于保持力积分法与模态叠加法结果的一致性非常重要。

与传统有限元建模方法不同，由于模型中刚柔耦合的存在，方程的建立与桨叶剖面振动载荷的计算需要作出相应的修改。气动力虚功需要以对应于刚体自由度与弹性自由度的广义力形式给出，因为它对刚体自由度同样做功。外力变分项可写成：

(3)

其中，R为质点在旋转桨毂坐标系中的坐标；TRD是矢量从变形后坐标系到旋转桨毂坐标系中的转换矩阵，因此也是刚体转角以及弹性变形的函数；φ为弹性扭转角；Fa与Ma分别指变形后坐标系中的气动力与力矩。实际计算中，它们不仅包含采用准定常气动模型得到气动环量力与力矩，还考虑了基于Greenberg薄板理论的非环量力与力矩分量。非环量部分是迎角变化率及弹性扭转的函数，提供气动阻尼，对气弹响应计算中的数值稳定性至关重要。

1.2 动力学方程

在得到动能变分、势能变分以及外力虚功后，桨叶非线性动力学微分方程组依据Hamilton原理建立起来。桨叶运动由形函数离散成广义节点自由度q的形式(由刚体自由度与Chopra 15自由度[9]组成)。动力学方程的最终形式可表示为：

(4)

(5)

为避免Jacobi阵的求逆操作，通过求解关于方程单步迭代差值的线性方程组推进牛顿迭代算法[11]，减少了程序计算量。同时，修改的牛顿法只在程序开始时计算Jacobi阵并在随后的迭代中当作常量处理；但在本文中为保证收敛方向，当收敛速度变慢时程序会适时更新Jacobi阵。

1.3 载荷计算方法

旋翼动力学中一般使用三种计算方法预测桨叶结构载荷。力积分法[1,8,12]将剖面结构载荷沿外段桨叶从计算参数点(即待求载荷的径向位置)到桨尖积分得到参考点结构载荷。因而这种方法会累积所有可能的由数值算法或剖面特性差异引起的误差。根据Hamilton原理与牛顿定律的等价性，桨叶惯性载荷可直接由动能变分对应的自由度变分分量得到[8]。即惯性载荷可表示为：

(6)

基于这种相对简单的对应关系，求解动能变分中的非线性项TNL子函数可在力积分法计算载荷时被再次调用，极大地降低了载荷程序实现的复杂度。由于质量阵也来自于动能变分，一种新的计算惯性载荷的方法可以推导出来：

(7)

其中，Q是由所有桨叶运动自由度组成的列向量，在本文中写成(θ1,θ2,θ3,u,v,w,φ,v′,w′)T的形式，二阶导数为其对应加速度项；M为桨叶运动自由度对应的离散之前的质量阵；g为体现式(6)中关系的抽象函数。

力积分法计算中，外段剖面对计算参考点x0的弯矩包括两部分：外段剖面力对计算参考点的力矩以及外段剖面内力对于本地轴线形成的弯矩；其中剖面力又由惯性载荷与气动外力组成。最终弯矩沿桨叶展向积分得到。以变形前坐标系中的挥舞弯矩为例：

(8)

其中，Δw=w-w0, Δx=(x+u)-(x0+u0)，下标0表示该变量是对应于计算参考点处的运动变形。因为传统Legendre-Gauss数值积分方法在有瞬态气动力作用时收敛速度较慢，本文将Richardson外插技术与之结合，发展了一种新的数值积分算法，以加快积分的收敛速度并尽可能地减少累积误差。

曲率法[2,4]的思想可以简洁地解释为弯矩只与梁的弯曲程度有关。该方法需要对桨叶弹性变形作二次偏导数，因而基于传统模态法求解响应时须保留足够多的模态才能达到合理的精度。本文采用隐式方法在位形空间中求解方程，因此不存在这种因模态截断引起的问题。仍以挥舞弯矩为例：

(9)

其中，EIy为挥舞 刚度；EA为拉伸刚度；θt指剖面预扭角；eζ指剖面重心在剖面坐标系中摆振方向的偏移。值得注意的是，与式(8)不同，式(9)中的弯矩是在变形后坐标系中度量的；因此在与力积分法比较时，需要将二者的载荷矢量统一到同一坐标系中。由于飞行实测载荷数据是在变形后坐标系中得到的，实际计算时式(8)的结果同样需要经过TDU转换到变形后坐标系中，而式(9)则不必。

反力法[12]基于有限元原理，因此结果取决于动力学有限元建模的精度。在有限元组集过程中，节点内力在相邻单元之间是相互抵消的；但在最终方程求解之后，由组集之前的相应单元矩阵形成的梁单元方程是非平衡系统，其余量即为相邻单元作用于节点的结构载荷。反力法不能有效预测有限元非节点位置处的结构载荷。尽管如此，反力法仍是一种计算代价较小的载荷计算方法。节点反力可写为：

(10)

其中，下标EL指该变量是原始单元矩阵或向量，并非直接取自方程组集之后的总体矩阵或向量。挥舞弯矩是Freaction中对应于w′自由度的分量，摆振弯矩即是对应于v′自由度的分量。类似的，反力法的载荷结果是在变形前坐标系中度量的，在与别的方法以及实验结果比较时需要进行转换。

2 算例与讨论

采用BO105模型桨叶[13]作为第一个算例对比在无气动外力作用下三种载荷计算方法的结果。BO105直桨叶安装在无铰式桨毂上。给定桨尖挥舞方向初速度10 m/s，径向位置21.6%R处的挥舞弯矩、摆振弯矩以及扭转弯矩分别如图2/3/4 所示。可以看出，即使在这种剧烈振动状态下，三种方法的结果仍能吻合得很好。这说明了三种方法在对于无铰式旋翼桨叶纯结构振动载荷具有相同的精度。

第二个算例采用SA349/2[15]直升机飞行实测载荷数据验证三种载荷计算方法。SA349/2是装有非线性减摆器的铰接式旋翼，其桨叶具有先进几何外形(预扭及非对称翼型等)及非线性剖面特性(剖面重心偏置等)。选取试验中的飞行状态3 (μ=0.198,CT/σ=0.062) ，比较研究三种载荷计算方法。

图2 挥舞弯矩，剖面位置21.6%R, BO105

图5 挥舞弯矩，剖面位置80%R, SA349/2

图6 摆振弯矩，剖面位置46%R, SA349/2

径向位置80%R处挥舞弯矩及46%R处摆振弯矩分别如图5/6所示。三种方法都能预测到弯矩的大致趋势。对于挥舞弯矩，力积分法在整个周期上都预测偏高；并且弯矩随时间变化的曲线较平滑，相位细节有损失。这种预测偏高的现象同样在摆振弯矩出现，如图6所示。对于桨叶中段的摆振弯矩，曲率法与反力法都比力积分法获得了更为合理的结果。在桨叶前进边，三种方法都预测过高。除了线性当量阻尼模型对于非线性减摆器建模的不准确以外，准定常气动模型可能不足以预测这种飞行状态下桨叶前行边气动载荷的计算。因为中等前飞状态下，桨涡干扰的存在降低了准定常气动模型的气动力预测精度。这种情况需要采用非定常气动模型[4]或计算流体动力学模型对气动力进行建模[8,15]。

3 结 论

本文建立了考虑刚体与弹性运动耦合的旋翼刚柔耦合动力学模型。依据Hamilton原理与牛顿定律的等价性，发展了一种易于实现的力积分方法。为加快积分收敛速度及减少积分误差，开发了一种新的基于Richardson外推法的Legendre-Gauss数值积分算法。着重比较直升机常用的三种载荷计算方法，以无铰式BO105模型桨叶以及铰接式SA349/2直升机为研究对象，分析对比了三种载荷方法的预测精度与适应范围。基于上述研究，本文总结如下：

(1) 在纯结构系统给定初始仿真运动条件下，由于不考虑气动力，三种载荷计算方法都能胜任结构振动载荷的预测，包括桨叶剖面挥舞弯矩、摆振弯矩及扭转力矩。

(2) 曲率法与反力法易于实现。在有气动外力作用条件下，二者较力积分法能获得频率更为丰富幅值更为合理的载荷结果。反力法无法预测有限元非节点处的结构载荷。依靠相邻节点之间的形函数插值得到的载荷是不准确的，桨叶建模时需要人为添加有限元节点。

(3) 气动模型对于旋翼振动载荷至关重要，气动力的预测能力需要进一步提高。从对比结果看，三种载荷计算方法都未能准确预测桨盘前行边的振动载荷，这与前行边的桨叶剖面气动力预测精确不够有关。

参 考 文 献

[1]Bialawa Richard L. Blade stress calculation-mode deflection vs. force integration[J]. Journal of the American Helicopter Society, 1978, 23(3),10-16.

[2]Maier T H. An examination of helicopter rotor load calculations[C]. The American Helicopter Society National Specialists’ Meeting on Rotorcraft Dynamics, Arlington, Texas, November 1989.

[3]Lim J W, Analytical investigation of UH-60A flight blade airloads and loads data[C]. The 51th Forum of the American Helicopter Society, Fort Worth, Texas, May 1995.

[4]WANG Hao-wen, GAO Zheng. Rotor vibratory load prediction based on generalized forces[J]. Chinese Journal of Aeronautics, 2004,17(1):28-33.

[5]韩东,高正，王浩文,等. 直升机桨叶刚柔耦合特性及计算方法分析[J]. 航空动力学报, 2006,21(1):36-40.

HAN Dong, GAO Zheng, WANG Hao-wen. Analysis of the computational method and the helicopter blade characteristics with coupled rigid and elastic Motions [J]. Journal of Aerospace Power, 2006,21(1):36-40.

[6]Bousman W G, Maier T H, An investigation of helicopter rotor blade flap vibratory loads[C]. The 48th Annual Forum of the American Helicopter Society, Washington, D.C., June 1992.

[7]Hodges D H, Dowell E H. Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades[R]. NASA TN D-7818, December,1974.

[8]Datta A, Fundamental Understanding, Prediction and validation of rotor vibratory loads in steady level flight[D]. PhD thesis in Aerospace Engineering, 2004.

[9]Chopra I, Sivaneri N T, Aeroelastic stability of rotor blades using finite element analysis[R]. NASA-CR-166389, August 1982.

[10]虞志浩, 杨卫东, 张呈林. 基于Broyden法的旋翼多体系统气动弹性分析[J]. 航空学报, 2012,33(12):2171-2182.

YU Zhi-hao, YANG Wei-dong, ZHANG Cheng-lin. Aeroelasticity analysis of rotor multibody system based on Broyden method[J]. Chinese Journal of Aeronautics, 2012,33(12):2171-2182.

[11]Ascher U M,Petzold L R. Computer methods for ordinary differential equations and differential-algebraic equations[M]. SIAM press, 1998.

[12]Floros M W, Elastically tailored composite rotor blades for stall alleviation and vibration reduction[D]. PhD thesis in Aerospace Engineering, December 2000.

[13]Bir G, Chopra I. University of maryland advanced rotorcraft code (UMARC) theory manual[M].Vol3.July 20, 1994.

[14]Heffernan R, Gaubert M, Structural and aerodynamic loads and performance measurements of an SA349/2 Helicopter with an Advanced Geometry Rotor[R]. NASA TM-88370, 1986.

[15]Abhishek A, Analysis, Validation, Prediction and fundamental understanding of rotor blade loads in an unsteady maneuver[D]. PhD thesis in Aerospace Engineering, 2010.