0 引言

由于旋转效应的影响，直升机由阵风引起的气动弹性响应不同于固定翼飞机。对直升机进行阵风响应分析，既要考察孤立桨叶也考充分考虑机体的耦合效应，尤其是机体的惯性所带来的动载荷效应^[1]。当遭遇阵风时，考虑了机体耦合效应的桨叶可能会表现出不同的特性。为了提升直升机飞行品质，需要对结构的应力和桨盘拉力进行精确预测，这就需要一个能够将机体耦合效应体现在内的计算模型。对于稳定飞行状态，通常使用孤立桨叶系统来考察桨盘上气动力和结构载荷的变化。但是对于遭受阵风情况，机体的惯性也十分重要，此时孤立桨叶系统所计算的载荷显得保守。

目前，诸多阵风响应特性的模型多数以孤立桨叶为研究对象。Drees^[2]借助单桨叶模型考察了桨叶在不同阵风形式下的响应特性。他指出不同形式的阵风响应形式区别明显。Arcidiacono^[3]借助定常入流理论分析了有铰桨叶在正弦形式阵风下的响应。他指出MIL-S-8698中的阵风修正因子偏保守。Elliott^[4]使用基于傅里叶变换的有限元法和Theodorsen理论分析了桨叶的随机阵风响应。Xiang^[5]使用梁理论和准定常气动力模型考察了无铰的复材桨叶在前飞过程中的气动弹性响应。Yasue^[6]使用基于片条的准定常气动力模型计算气动力结合谐波平衡法考察了阵风响应特性并使用孤立桨叶风洞模型对计算结果进行了验证。Azuma^[7]运用叶素理论考察了考虑挥舞自由度的桨叶阵风响应并和孤立桨叶风洞实验进行了对比。由于桨叶工作的气动力条件较为复杂，尤其桨根和桨尖处由于攻角的变化而产生的气动力的变化更加频繁。因此，基于定常和准定常状态的气动力模型不能满足详细设计阶段对气动力精确度的要求，非定常气动力理论大显身手。Chopra^[8]借助中等变形梁理论建立了桨叶的刚弹耦合模型并结合非定常气动力考察了直升机悬停和前飞特性。董凌华^[9-10]借助中等变形梁理论和非定常气动力考察了倾转旋翼耦合系统的气动弹性特性。Yeo^[11]借助基于非定常气动力的CAMRAD II程序计算了多种旋翼的结构载荷。Park^[12]使用CFD-CSD松耦合策略发展了桨叶结构载荷的分析方法。王林鹏和戴玉婷^[13-14]借助中等变形梁和大变形梁理论结合非定常气动力考察了孤立桨叶的阵风响应特性。考察结果表明，中等变形梁理论结合非定常气动力的模型计算精度较好。Chopra^[15]建立了考虑机身耦合影响的弹性桨叶模型来考察其气动弹性响应，但是对于平稳飞行状态来说，机身惯性对桨叶载荷影响较小。但是在时域中，当直升机突然遭受阵风，由于机身惯性的存在，其耦合效应将对桨叶载荷产生明显影响。

为考察机身耦合效应对阵风响应的影响程度，本文借助基于Hamilton原理的中等变形梁理论，结合Pitt-Peters动态入流理论和Leishman Beddoes(L-B)动态失速模型^[16-18]，并添加机身俯仰和滚转自由度惯性效应，建立用于考察整机阵风响应特性的耦合模型。并借助建立的机身耦合模型，考察机身耦合模型在遭受阵风时的响应特性，以期为工程应用提供借鉴参考。

1 考虑机身耦合的气动弹性模型

图 1为考虑机体惯性效应的桨叶结构与气动力耦合策略图。阵风和动态入流作为速度矢量被添加到气动力模型中，该速度矢量将最终影响桨叶的时域中气动力攻角，从而影响桨叶气动力。机身俯仰和滚转惯性效应将被叠加进入桨盘动力学方程中，从而影响桨叶结构动力学方程的计算条件。结构模型根据计算条件计算桨叶挥舞和扭转自由度变形和速度，计算结果将作为反馈传递给气动力计算模型，进而用于更新气动力计算的初始条件，完成耦合循环。

1.1 气动力模型

气动模型部分由动态入流模型和动态失速模型组成，Pitt-Peters模型^[18]作为动态入流模型参与建模。Leishman-Beddoes^[19-22]模型是时域半经验动态失速模型。此模型中，流场被分为两个阶段——附着流阶段和分离流阶段。最终分离过程中的法向力系数计算公式的推导过程在之前发表论文中已有详细论述，不再赘述。最终法向力可写成分离状态下的法向力系数C_N^f(t)和脱落涡的附加法向力系数C_N^v(t)的和：

$ {C_N}\left( t \right) = C_N^f\left( t \right) + C_N^v\left( t \right) $

(1)

俯仰力矩系数为：

$ {C_m}\left( t \right) = C_m^f\left( t \right) + C_m^v\left( t \right) $

(2)

1.2 结构动力学模型

桨叶的结构动力学模型根据Hamilton方程进行推导^[23-26]。方程可写成：

$ \delta \mathit{\Pi } = \int\limits_{{t_1}}^{{t_2}} {\left( {\delta U - \delta T - \delta W} \right){\rm{d}}t} = 0 $

(3)

式中，U为应变能, T为动能, W为外力虚功。

桨叶的有限元形式可写成：

$ \delta \mathit{\Pi } = \int\limits_{{\psi _1}}^{{\psi _2}} {\sum\limits_{i = 1}^n {\left( {\delta {U_i} - \delta {T_i} - \delta {W_i}} \right){\rm{d}}\psi } } = \int\limits_{{\psi _1}}^{{\psi _2}} {\sum\limits_{i = 1}^n {\left( {{\Delta _i}} \right){\rm{d}}\psi } } = 0 $

(4)

式中n为每个桨叶的梁单元数量。

图 2为有限单元节点的自由度分布图。

在每个单元中，考虑了轴向拉伸自由度u, 挥舞弯曲、摆振弯曲自由度v_ss、w，以及扭转自由度$φ$。则有：

$ \left\{ {\begin{array}{*{20}{c}} {u\left( s \right)}\\ {{v_{ss}}\left( s \right)}\\ {w\left( s \right)}\\ {\varphi \left( s \right)} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{H_u}}&0&0&0\\ 0&H&0&0\\ 0&0&H&0\\ 0&0&0&{{H_\varphi }} \end{array}} \right\}{q_i} $

(5)

式中H为Hermite多项式, q_i定义如下:

$ q_i^{\rm{T}} = \left\{ {\begin{array}{*{20}{c}} {{u_1}}&{{u_2}}&{{u_3}}&{{u_4}}&{{v_{ss1}}}&{{{v'}_{ss1}}}&{{v_{ss2}}}&{{{v'}_{ss2}}}\\ {{w_1}}&{{{w'}_1}}&{{w_2}}&{{{w'}_2}}&{{\varphi _1}}&{{\varphi _2}}&{{\varphi _3}}&{{\varphi _3}} \end{array}} \right\} $

(6)

在离散形式下，每个单元Δ_i写成：

$ {\Delta _i} = \delta q_i^{\rm{T}}{\left( {M\ddot q + C\dot q + Kq - \mathit{\boldsymbol{F}}} \right)_i} $

(7)

式中，M为单元质量矩阵, C为单元阻尼矩阵, K为单元刚度矩阵, F为广义力向量。

将桨叶的每个单元的应变能和动能叠加在一起，可得到桨叶的整体应变能和动能方程。最终整个桨叶的积分形式可以写成：

$ \delta \mathit{\Pi } = \int\limits_{{\psi _1}}^{{\psi _2}} {\sum\limits_{i = 1}^n {\left[ {\delta q_i^{\rm{T}}{{\left( {M\ddot q + C\dot q + Kq - F} \right)}_i}} \right]{\rm{d}}\psi } } = 0 $

(8)

1.3 阵风模型

首先假设阵风在空间内不变，即阵风在桨叶各个位置的速度一致。阵风强度随时间变化。考虑两种典型的阵风形式——脉冲阵风和斜坡阵风。

(1) 脉冲阵风。图 3(a)为脉冲形式阵风。阵风方向为垂向，阵风强度为W，作用时间为t₁-t₀。

(2) 斜坡阵风。图 3(b)为斜坡形式阵风。其最大强度为W。

1.4 桨叶气动弹性响应

桨叶的气动弹性效应是气动和结构的耦合计算。

首先，动态入流谐波分量速度可以通过求解动态入流模型的状态空间方程计算^[18]。通常求解的谐波分量均为无量纲形式, 计算所得诱导速度v也是无量纲结果, 总的诱导速度需要无量纲结果进行还原：v=vΩR。

图 4为叶素截面各参数的几何关系图。其中入流角β_in可通过下式计算：

$ {\beta _{{\rm{in}}}} = \arctan \left[ {\frac{{\nu - {v_{\rm{h}}}}}{{\mathit{\Omega r + }{\mathit{V}_f}\cos \left( {\mathit{\Omega t} - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)}}} \right] $

(9)

式中，v_h为桨叶挥舞速度，可通过结构模型得到；N为桨叶的数量; j为每个桨叶的编号。

由于挥舞变形和扭转变形是桨叶旋转的主要变形，而摆振变形数值较小，所以此次研究中不考虑摆振变形。

气动攻角α的表达式可以写成：

$ \alpha = {\theta _0} - {\beta _{{\rm{in}}}} + \varphi $

(10)

式中, $φ$为桨叶的扭转变形，可通过结构模型计算；θ₀为桨叶操纵变距，其表达式为：

$ {\theta _0} = {\theta _{75}} + {\theta _{tw}}\left( {\frac{x}{R} - 0.75} \right) + {\theta _{1c}}\cos \psi + {\theta _{1s}}\sin \psi $

(11)

其中，θ₇₅为总距, θ_tw为预扭角, θ_1c为横向周期变距, θ_1s纵向周期变距。

1.5 计算过程

在计算诱导速度之前，所有计算参数都需要先进行初始化。根据入流条件，更新气动力模型中的攻角。将攻角带入气动力模型计算法向力系数C_N和力矩系数C_m。

叶素法向力F_N、弦向力F_C和力矩M可通过下式进行计算：

$ \begin{array}{l} {F_N} = 0.5\rho V_l^2{C_N}{c_l}\\ {F_C} = 0.5\rho V_l^2{C_C}{c_l}\\ M = 0.5\rho V_l^2{C_m}c{l^2} \end{array} $

(12)

式中，ρ为空气密度, V_l为当地速度, l为叶素长度。

叶素升力F_L和阻力F_D以及桨叶垂直方向拉力F_Z可表示为：

$ \begin{array}{l} {F_L} = {F_N}\cos \alpha - {F_C}\sin \alpha \\ {F_D} = {F_N}\sin \alpha + {F_C}\cos \alpha \\ {F_Z} = {F_L}\cos {\beta _{{\rm{in}}}} - {F_D}\sin {\beta _{{\rm{in}}}} \end{array} $

(13)

通过对叶素上力和力矩进行累加，可求得单个桨叶上载荷，将所有桨叶载荷叠加将求得桨盘拉力F_{TT_all}以及横向力分量F_{LL_all}和纵向力分量F_{MM_all}，各分量表达式可表示为：

$ {F_{TT{\rm{\_all}}}} = \sum\limits_1^N {\sum\limits_1^n {{F_Z}} } $

(14)

$ {F_{LL{\rm{\_all}}}} = \sum\limits_1^N {\sum\limits_1^n {\left[ {{F_Z}r\cos \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)} \right]} } $

(15)

$ {F_{MM{\rm{\_all}}}} = \sum\limits_1^N {\sum\limits_1^n {\left[ {{F_Z}r\sin \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)} \right]} } $

(16)

对各力分量进行无量纲处理，以便重新输入入流模型, 开始新一轮计算。其中，桨盘拉力系数C_T, 滚转力矩系数C_l和俯仰力矩系数C_m可由下式计算：

$ \begin{array}{l} {C_T} = \frac{{{F_{TT{\rm{\_all}}}}}}{{\rho {\rm{ \mathsf{ π} }}{R^2}{{\left( {\mathit{\Omega }R} \right)}^2}}}\\ {C_l} = \frac{{{F_{LL{\rm{\_all}}}}}}{{\rho {\rm{ \mathsf{ π} }}{R^3}{{\left( {\mathit{\Omega }R} \right)}^2}}}\\ {C_m} = \frac{{{F_{MM{\rm{\_all}}}}}}{{\rho {\rm{ \mathsf{ π} }}{R^3}{{\left( {\mathit{\Omega }R} \right)}^2}}} \end{array} $

(17)

通过求解状态空间方程，可计算v₀、v_c、v_s，则入流模型的诱导速度可以通过下式计算：

$ \begin{array}{*{20}{c}} {{v_{ji}} = \left[ {{{\bar \nu }_0} + {{\bar \nu }_{\rm{c}}}\frac{{{r_i}}}{R}\cos \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right) + } \right.}\\ {\left. {{{\bar \nu }_s}\frac{{{r_i}}}{R}\sin \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)} \right]\mathit{\Omega }R} \end{array} $

(18)

当考虑阵风时，将阵风看做模型中的速度分量，添加进模型，则此时入流角可以写成：

$ {\beta _{{\rm{in}}}} = \arctan \left[ {\frac{{\nu - {v_h} - {v_{{\rm{gust}}}}}}{{\mathit{\Omega }r + {V_f}\cos \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)}}} \right] $

(19)

式中，v_gust为阵风速度。

1.6 机体桨叶耦合

本次主要考虑机体的惯性对桨叶变形和载荷变化的耦合效应，在耦合模型中考虑了机体俯仰自由度和滚转自由度的惯性。定义桨盘在俯仰方向的倾角为θ_z0、滚转方向的倾角为θ_g0。当直升机处于前飞状态，将前飞速度分解成水平分量V_x和垂直分量V_z, 则桨盘倾角俯仰自由度的大小为：

$ {\alpha _S} = 0.5{\rm{ \mathsf{ π} }} - \arctan \left( {\frac{{{V_x}\cos {\theta _{z0}} - {V_z}\sin {\theta _{z0}}}}{{{V_x}\sin {\theta _{z0}} + {V_z}\cos {\theta _{z0}}}}} \right) $

(20)

V_f可由下式计算：

$ {V_f} = {V_x}\cos {\theta _{z0}} + {V_z}\cos {\theta _{z0}} $

(21)

机身的各自由度加速度计算式为：

$ {A_Z} = \frac{{{F_{TT{\rm{\_all}}}}\cos {\theta _{z0}} - {M_{ji}}g}}{{{M_{ji}}}} $

(22)

$ {A_X} = \frac{{{F_{TT{\rm{\_all}}}}\sin {\theta _{z0}} - {F_{zu}}}}{{{M_{ji}}}} $

(23)

$ {A_{{\rm{alaph}}}} = \frac{{{F_{MM{\rm{\_all}}}}}}{{{I_{fy}}}} $

(24)

$ {A_{{\rm{gunzhuan}}}} = \frac{{{F_{LL{\rm{\_all}}}}}}{{{I_{gz}}}} $

(25)

其中，A_Z为垂向加速度，A_X为水平方向加速度，A_alaph为俯仰角加速度，A_gunzhuan为滚转角加速度，M_ji为直升机机身质量，F_zu为机身阻力，I_fy为俯仰惯性矩，I_gz为滚转惯性矩。

$ \begin{array}{*{20}{c}} {{F_{{\rm{mid}}}} = {F_z}\cos \left[ {{\theta _{z0}}\sin \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)} \right] \cdot }\\ {\cos \left[ {{\theta _{g0}}\cos \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)} \right]} \end{array} $

(26)

方程(14)-(16)可写成：

$ {F_{TT{\rm{\_all}}}} = \sum\limits_1^N {\sum\limits_1^n {{F_{{\rm{mid}}}}} } $

(27)

$ {F_{LL{\rm{\_all}}}} = \sum\limits_1^N {\sum\limits_1^n {\left[ {{F_{{\rm{mid}}}}r\cos \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)} \right]} } $

(28)

$ {F_{MM{\rm{\_all}}}} = \sum\limits_1^N {\sum\limits_1^n {\left[ {{F_{{\rm{mid}}}}r\sin \left( {\mathit{\Omega }t - \frac{{2j{\rm{ \mathsf{ π} }}}}{N}} \right)} \right]} } $

(29)

桨盘状态可以根据下式进行更新：

$ {V_Z} = {V_Z} + {A_Z}\Delta t $

(30)

$ {V_X} = {V_X} + {A_X}\Delta t $

(31)

$ {\theta _{z0}} = {\theta _{z0}} + 0.5{A_{{\rm{alpha}}}}{\left( {\Delta t} \right)^2} $

(32)

$ {\theta _{g0}} = {\theta _{g0}} + 0.5{A_{{\rm{gunzhuan}}}}{\left( {\Delta t} \right)^2} $

(33)

图 5为整个计算流程图。

2 模型验证 2.1 弧立桨叶模型验证

为保证模型计算精度，对所建模型进行算例验证。以无铰式旋翼bo105为例, 旋翼桨叶数为4，旋翼半径4.9377 m, 转速383 r/min(40.12346 rad/s), 桨叶翼型截面为NACA0015，弦长0.395 m，实度σ=0.1，C_T/σ=0.07。线密度m₀=6.4683 kg/m, 预扭角β_p=0°。具体的剖面性能参数如表 1。

首先对结构部分做了验证，表 2为本文建立的中等变形梁模型的计算结果与文献计算结果的对比。从表 2可以看出，本文的计算结果与文献中的计算结果相比最大误差最大为3.14%，说明本文的中等变形梁模型有足够的精度，可以用来研究大展弦比机翼的气动弹性的响应特性。

对气动力模型的验证以及操纵角配平的验证见文献[14]，文中对计算气动力计算的验证有详细的描述。验证耦合模型时，计算了前进比为0.2的飞行状态，来考察孤立桨叶模型的建模的合理性。桨叶的旋转速度为383 r/min，操纵角为：θ₀=7.1°+1.2°cosψ－3.1°sinψ。图 6为本文计算结果和文献[26]中同样计算条件的计算结果对比图。从图中可以看出本文计算结果和文献[26]计算结果曲线重合度较好，可以说明孤立桨叶模型建模正确。

2.2 机身桨叶耦合模型验证

当直升机在定直飞行过程中，由于自身处于稳定状态，机身由惯性产生的动载荷对桨叶的载荷和变形影响较小，此时可近似认为考虑机身耦合的模型退化成孤立桨叶模型，利用这一点，借助孤立桨叶模型来验证机身耦合效应模型的可信度。计算条件下旋翼转速为383 r/min。耦合模型机身质量为3000 kg, 俯仰自由度惯性矩为18000 kg ·m ·s², 滚转自由度惯性矩为8200 kg ·m ·s²。计算了前进比为0.2和0.35两种前飞状态下的耦合模型的翼尖变形和翼根剪切力。图 7和图 8分别为前进比为0.2和0.35情况下，孤立桨叶模型和机身耦合模型的计算结果对比图，可以看出考虑耦合效应的模型的结果曲线和孤立桨叶模型曲线吻合度较高，可以说明机身耦合模型建模合理，可以用于考察直升机遭受阵风情况下桨叶气动弹性响应问题。

3 阵风响应分析

本文主要考察机身惯性耦合效应在阵风作用下的影响，所以借助孤立桨叶模型的计算结果与机身耦合模型的计算结果进行对比。分别计算了悬停状态、前进比为0.2和0.35前飞状态下机体耦合模型的阵风响应情况。在前进比为0.2时, 操纵角为θ₀=7.1°+1.2°cosψ－3.1°sinψ。在前进比为0.35时, 操纵角为θ₀=10°+2°cosψ－8°sinψ。

3.1 悬停状态的阵风响应

在悬停状态下，用孤立桨叶模型和机身耦合模型分别计算遭受脉冲阵风和斜坡阵风时桨叶的响应情况。

图 9为阵风模型数据图。图 10(a)和图 10(b)分别为两种模型在遭受两种阵风情况下桨盘的拉力系数曲线。考虑机体影响的桨盘拉力系数和孤立桨叶的计算结果相差较小，但是在遭受阵风时，这一差值会有所增加，尤其考虑机体耦合效应的模型，桨盘拉力系数下降明显。图 11(a)和图 11(b)为两种模型的桨根剪切力的对比图。剪切力有与桨盘拉力系数相似的趋势，在遭受阵风时，机体耦合模型中桨根受剪切力较小。从上述对比情况来看，机身的耦合效应，会使阵风情况下的桨盘拉力下降，同时桨叶受载情况缓解。此处应重点关注拉力的下降，这将影响直升机的飞行稳定性。

3.2 前飞状态的阵风响应

在前飞状态下同样考察了脉冲和斜坡阵风情况下桨叶气动弹性系统响应特性。在图 12中给出前飞状态阵风模型图。图 13为由阵风引起的翼尖挥舞位移的时域曲线。对比结果可以看出，孤立桨叶模型计算的翼尖位移结果较耦合模型保守。图 14为两种阵风情况下的翼根剪力曲线，从图中可以看出机身耦合模型的计算结果小于孤立桨叶模型的结果，机身耦合效应对计算结果影响较为明显。图 15为两种阵风模型下翼尖的挥舞位移的频率响应曲线。从图 15(a)和图 15(b)中可以看出，在遭受脉冲阵风时，转速的四倍频会有明显响应，这一情况在孤立桨叶模型中并不明显。图 16为前进比为0.2时，在两种阵风情况下，两种计算模型中的桨盘拉力系数的变化曲线。从图中可以看出，在遭受阵风时，桨盘拉力系数均会变小，考虑机身耦合效应的模型的计算结果更加明显。考虑机身耦合效应模型的计算结果较孤立桨叶模型计算结果小17%。

图 17和图 18为高前进比0.35时，直升机的变形响应和桨盘拉力系数变化。从图中可以看出在高前进比情况下桨尖变形较低前进比大。而高前进比情况下，桨盘在遭受阵风后，拉力系数损失较为严重。孤立桨叶模型所反映的计算结果过于保守，没有将损失情况完全刻画出来，将不利于直升机稳定性判定。以上结果可以看出，在前飞状态下，尤其飞行速度较大情况下，在遭受阵风时，桨盘拉力会明显损失，应充分考虑机体耦合效应，以获得精确的计算结果。

4 结论

从孤立桨叶模型和机身耦合模型计算结果对比情况来看，有以下结论：

(1) 在悬停状态遭受阵风时，桨盘的上的拉力系数和剪切力均减小。考虑机身耦合效应的计算结果小于孤立桨叶模型。

(2) 在遭受阵风时，桨叶的桨尖变形增加。孤立桨叶模型计算的拉力系数结果比机身耦合模型计算结果大。在高前进比中，耦合效应对计算结果的影响较为明显，不可忽略。