大机动飞行是直升机，尤其是攻击和侦察直升机的重要飞行状态，开展直升机大机动飞行特性的研究，对于直升机飞行品质、飞行载荷以及气动噪声等的研究具有重要意义。大机动实时飞行仿真是研究大机动飞行特性的有效手段，它不仅具备较高的准确性，同时具备快速、实时的计算效率，满足日常工程设计的需求。

与稳定和有限机动飞行状态相比，直升机在大机动飞行过程中，操纵输入、运动速度、角速度以及姿态存在快速、大幅的变化，不仅导致旋翼气动力、桨叶运动存在非定常、非线性的特征，更重要的是，还会诱导旋翼的尾迹涡系产生额外的拉伸、压缩、倾斜、弯曲等动态畸变，严重影响旋翼桨盘入流分布和各部件气流环境，改变直升机大机动飞行中的动态响应特性。例如，机体俯仰滚转角运动诱导的尾迹畸变已被证明是导致俯仰/滚转交叉耦合响应“反号”问题的主要原因^[1-3]。

目前，在直升机机动飞行动力学建模中采用的旋翼尾迹模型或入流模型主要有尾迹畸变增广的动量理论^[3]、尾迹畸变增广的动态入流模型^[4-6]以及非定常自由尾迹模型^[7-13]。其中，尾迹畸变增广的Pitt-Peters动态入流模型由于其模型简单，计算效率高等优点而被广泛采用。然而，尾迹畸变增广的动量或动态入流模型中，机体非定常运动诱导的尾迹畸变采用预定的涡管尾迹模型描述，无法计入尾迹自诱导引起的各种畸变。另外，在该类模型中定义了尾迹弯曲系数K_Re，即尾迹畸变诱导的桨盘入流梯度与机体俯仰滚转角运动直接诱导的桨盘入流梯度之比，用于控制尾迹弯曲畸变与诱导入流梯度之间的正比关系。但不同研究者通过理论分析得到的K_Re不完全相同^{[6, 11, 14]}，另外文献^{[12, 14]}中分别通过理论和试验研究表明K_Re的大小随旋翼拉力系数、俯仰滚转角速率等状态参数变化。文献^[6]中采用尾迹畸变增广的动态入流模型进行飞行动力学建模，以UH-60直升机俯仰滚转机动飞行为算例，结果表明不同飞行状态需通过设置不同大小的K_Re使得异轴响应计算结果与飞行试验结果基本一致。因此尾迹畸变增广的动态入流模型的通用性受到限制。文献^{[10, 12]}中的研究表明，非定常自由尾迹模型能够较准确地捕捉机体非定常运动诱导的尾迹畸变和尾迹自诱导畸变，正确计算桨盘入流分布，但是在悬停和小速度飞行状态非定常自由尾迹模型的数值稳定性和收敛性的鲁棒性较差，并且其计算量相对较大，特别是难以实现实时机动飞行仿真。

针对上述现状，本文提出一种既能模拟大机动飞行中旋翼尾迹的拉伸、压缩、倾斜、弯曲等动态畸变行为，又能满足实时仿真快速、高效需求的旋翼动态尾迹模型。该模型采用具有5个刚体运动自由度和1个半径伸展收缩自由度的涡环单元来表示旋翼尾迹涡系，有效地简化了尾迹涡系运动描述的复杂程度，并局部采用涡环对空间点诱导速度的解析表达式来进行诱导速度场计算，极大地提高了计算效率。最后，以悬停状态旋翼拉力和俯仰角速度阶跃突增为算例，通过与尾迹畸变增广的Pitt-Peters动态入流模型的对比分析，验证了本文模型的正确性和实时性。

1 旋翼动态涡环尾迹建模 1.1 物理模型

本文在旋翼无限片桨叶假设的作用盘理论和经典固定尾迹模型的基础上，只考虑旋翼尾迹中圆柱涡系环量的水平分量，将圆柱涡系或涡管离散成若干个带状涡环，其长度分别为l₁,l₂,l₃，…，l_n，设涡管上涡强的线密度为γ，则第n个涡带的涡量为Γ_n=γl_n，然后将带状涡环简化为一个集中涡环，如图 1所示。当带状涡环长度l_n较小时，可较准确地模拟圆柱涡系在桨盘平面的诱导速度场。

1.2 动态涡环尾迹建模

为描述大机动飞行中尾迹的动态畸变，需建立尾迹涡环的生成、运动和涡核半径的变化模型。由于涉及涡环的运动，因此首先建立基于涡环中心的坐标系OX_OY_OZ_O和旋翼桨毂不旋转坐标系OXYZ。

旋翼桨毂坐标系原点在桨毂中心，XY在桨毂平面内，X轴在直升机纵向对称剖面内指向机头反方向，Z轴垂直XY平面向上，Y轴方向由右手法则确定。涡环坐标系原点位于涡环中心，X_OY_O平面在涡环平面内。Z_O轴垂直X_OY_O平面向上。定义涡环相对于桨毂坐标系的滚转和俯仰姿态角分别为θ_x和θ_y，则涡环坐标系到旋翼桨毂坐标系的转换关系为

$\left[ \matrix{ X \hfill \cr Y \hfill \cr Z \hfill \cr} \right] = \left[ {\matrix{ {cos{\theta _y}} & {sin{\theta _x}sin{\theta _y}} & {sin{\theta _y}cos{\theta _x}} \cr 0 & {cos{\theta _x}} & { - sin{\theta _x}} \cr { - sin{\theta _y}} & {sin{\theta _x}cos{\theta _y}} & {cos{\theta _x}cos{\theta _y}} \cr } } \right]\left[ \matrix{ {X_O} \hfill \cr {Y_O} \hfill \cr {Z_O} \hfill \cr} \right]$

(1)

(1) 涡环运动模型

涡环具有5个刚体运动自由度和1个半径伸展收缩自由度，其中刚体运动包括3个线运动和俯仰滚转2个角运动，涡环线运动由涡环中心点在桨毂坐标系内的位置坐标(x,y,z)和线运动速度(x′,y′,z′)描述，角运动由涡环相对于桨毂坐标系的俯仰角θ_y和滚转角θ_x以及俯仰滚转角速率θ_y′,θ_x′描述，涡环半径伸展收缩自由度由涡环半径r及其伸展收缩速率r′描述。第i个涡环运动状态量向量可以表示为R_i=[x_i,y_i,z_i,r_i,θ_xi,θ_yi]^T，整个旋翼尾迹运动状态量向量可以写为R=[R₁ ^T,…,R _i ^T,…,R _N ^T] ^T，其中N为涡环个数。第i个涡环运动方程可以表示为

$\frac{d{{R}_{i}}}{dt}={{u}_{i}}\left( R \right)+{{v}_{i}}({{u}_{h}},{{R}_{i}})$

(2)

式中：u_i=[u_xⁱ,u_yⁱ,u_zⁱ,u_rⁱ,u_θxⁱ,u_θyⁱ] ^T 为尾迹自诱导的涡环运动速度；v_i (u_h,R_i)为旋翼线运动和角运动诱导的相对运动速度；旋翼运动速度和角速度列向量为u_h=[u_h,v_h,w_h,p_h,q_h,r_h] ^T。

(2) 旋翼涡环的生成

由前文描述可知，涡环本质上是由离散的带状涡环简化得到，而旋翼单位时间拖出的带状涡环长度与环量可以写成

$l=\Delta t({{v}_{0}}+{{w}_{h}})$

(3)

$\Gamma =\Delta t\gamma ({{v}_{0}}+{{w}_{h}})$

(4)

式中：v₀为当前桨盘平面的平均诱导速度；Δt为时间步长；γ为旋翼圆柱涡系涡强水平分量。将旋翼单位时间拖出的带状涡带集中为新的涡环，就得到新涡环(i=1) 的状态量R₁。

(3) 涡环尾迹诱导的运动速度确定方法

首先将第i个涡环按等方位角确定M个节点，然后计算所有涡环对第i个涡环上所有节点处的诱导速度。对于涡环的线运动速度u_xⁱ,u_yⁱ,u_zⁱ由M个节点的速度平均值确定。为了确定u_rⁱ，u_θxⁱ和u_θyⁱ，将M个节点的诱导速度转换到第i个涡环的当地坐标系，u_rⁱ即为M个节点上涡环径向速度平均值，而对于u_θxⁱ及u_θyⁱ，则先利用刚体角运动速度分布的特点，由M个节点上的垂向运动速度分布，通过最小二乘法确定其涡环坐标系下的俯仰滚转角速率，再转换为涡环俯仰滚转姿态角变化率u_θxⁱ,u_θyⁱ。

(4) 涡环涡核模型

本文采用文献^[15]提供的涡核半径扩散公式，第i个涡环涡核半径随时间的变化为

${{r}_{c}}^{i}=\sqrt{{{r}_{0}}^{2}+4\alpha \upsilon \delta i\Delta t}$

(5)

式中：r₀为初始涡核半径；α=1.256 4为Oseen常数；υ为空气运动黏性系数；相对等效黏性系数δ=1+a₁Re_v，其中a₁为湍流黏性扩散系数，Re_v为涡雷诺数。

1.3 快速涡环诱导速度计算方法

采用2种方式实现涡环对任意点诱导速度的计算。第一种涡环按等方位角离散为若干个直涡段，通过直涡段对空间一点诱导速度叠加来计算涡环的诱导速度。该方法具有解析表达式，可方便地引入涡核修正，但计算量相对较大。第二种方法利用涡环对空间任一点的解析表达式计算空间诱导速度。对于当地涡环坐标系OX_OY_OZ_O，任取一点M极坐标为(ρ_M,θ_M,Z_M)，涡环半径为ρ_A，则涡环对点M的垂向与径向诱导速度分别为

$\begin{align} & {{V}_{z}}=\frac{\Gamma }{4\pi }(-\frac{4}{\sqrt{{{C}_{1}}+{{C}_{2}}}}{{E}_{1}}(\frac{\pi }{2},k)+ \\ & \frac{\text{ }2({{l}_{z}}^{2}+{{\rho }_{M}}^{2}-{{\rho }_{A}}^{2})}{{{({{C}_{1}}+{{C}_{2}})}^{3/2}}}{{E}_{3}}(\frac{\pi }{2},k)) \\ \end{align}$

(6)

$\begin{align} & {{V}_{r}}=-\frac{\Gamma {{l}_{Z}}}{8\pi {{\rho }_{M}}}(-\frac{4}{\sqrt{{{C}_{1}}+{{C}_{2}}}}{{E}_{1}}(\frac{\pi }{2},k)+ \\ & \frac{4{{C}_{1}}}{{{({{C}_{1}}+{{C}_{2}})}^{3/2}}}{{E}_{3}}(\frac{\pi }{2},k)) \\ \end{align}$

(7)

式中：C₁=ρ_A²+l_z²+ρ_M²；C₂=2ρ_Aρ_M；l_z=Z_M；k²=2 C₂/( C₁+ C₂)；E₁为第一类完全椭圆积分；E₃为第三类完全椭圆积分。由式(6，7) 可知，只需建立模数k的椭圆积分表，即可通过查表计算空间任一点诱导速度，提高诱导速度计算效率。由于推导出的积分公式很难引入涡核半径，因此综合第一和第二种方式计算空间任一点诱导速度，以保证精度和效率。

2 算例及初步验证

为了在算例分析中与尾迹畸变增广的Pitt-Peters动态入流模型进行初步对比验证，首先简单介绍该模型。Pitt-Peters动态入流模型假设旋翼桨盘诱导入流分布为

$\lambda \left( \bar{r},\psi \right)={{\lambda }_{0}}+{{\lambda }_{1c}}\bar{r}cos\psi +{{\lambda }_{1s}}\bar{r}sin\psi $

(8)

式中：λ₀为平均诱导入流系数；λ_1s为横向诱导入流梯度系数；λ_1c为纵向诱导入流梯度系数。其增广后的入流方程可以表示为^[6]

$M\left[ {\begin{array}{*{20}{c}} {{{\mathop \lambda \limits^\circ }_0}} \\ {{\lambda _{1s}}} \\ {{\lambda _{1c}}} \end{array}} \right] + V{L^{ - 1}}\left[ {\begin{array}{*{20}{l}} {{\lambda _0}} \\ {{\lambda _{1s}}} \\ {{\lambda _{1c}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{C_T}} \\ { - {C_L}} \\ { - {C_M}} \end{array}} \right]$

(9)

式中：M为显在质量矩阵;V为质量流量参数矩阵；C_T为气动拉力系数；C_L为气动滚转力矩系数；C_M为气动俯仰力矩系数；L=L＇+K_ReΔL为尾迹弯曲增广后的入流增益矩阵，其中，L＇为Pitt-Peters动态入流模型的原始增益矩阵，ΔL为入流-尾迹畸变耦合矩阵，ΔL与尾迹倾斜参数X，尾迹横向和纵向弯曲曲率κ_s,κ_c，尾涡间距参数S相关。为了描述尾迹动态畸变过程，该模型引入了如下一阶惯性环节

$\tau \left[ {\begin{array}{*{20}{c}} {\mathop X\limits^\circ } \\ S \\ {{\kappa _c}} \\ {{\kappa _s}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} X \\ S \\ {{\kappa _c}} \\ {{\kappa _s}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} X \\ S \\ {{\kappa _c}} \\ {{\kappa _s}} \end{array}} \right]_{qs}}$

(10)

式中：τ为时间常数矩阵。其具体细节可以参见文献^{[4, 6]}。为了将本文模型计算结果与上述模型进行对比，利用式(8) 将本文模型计算的桨盘诱导入流分布通过最小二乘法识别出λ₀，λ_1c，λ_1s。本文以BO-105直升机旋翼为例，其主要的参数为：半径4.94 m，旋翼转速40 rad/s，桨叶弦长0.27 m，桨叶片数为3片。

2.1 悬停状态算例

首先计算稳定悬停状态以验证本文模型的收敛性和桨盘平均诱导入流系数λ₀计算准确性。图 2给出了悬停状态尾迹几何形状及其收敛过程，拉力系数C_T为0.006。图 2(a)为尾迹的几何形状。图 2(b)中均方根值(Root mean square，RMS)是旋翼尾迹涡环状态量在前后两个旋转周期内的残差。可见，涡环尾迹呈向下收缩状态，与真实情况相符； RMS最后收敛趋近于0，表明尾迹几何形状已经收敛于周期稳态解。

表 1给出了悬停状态不同拉力系数涡环尾迹模型中桨盘平均诱导入流系数λ₀。计算结果与动态入流模型结果的对比。可见，该拉力变化范围内，涡环尾迹模型计算结果略大于动态入流模型。

2.2 悬停状态拉力系数突增算例

为验证本文模型对于平均诱导入流系数动态响应特性计算的正确性，以悬停状态拉力系数阶跃突增为算例(拉力系数在第1 s由0.006突增至0.008) ，进行对比验证。图 3(a)给出了λ₀动态响应时间历程计算结果的对比。由图可知，稳态及瞬态响应过程中两者基本一致。图 3(b~e)给出了该过程本文模型计算的尾迹几何形状动态响应过程。

2.3 悬停状态旋翼俯仰角速度突增算例

为验证本文模型对于机动飞行状态桨盘诱导入流梯度动态响应特性计算的正确性，以悬停状态旋翼俯仰角速度阶跃突增为算例，进行对比验证。图 4(a~c)给出了拉力系数C_T=0.006，第1 s突增抬头角速率q分别为10，20，30 °/s，至第3 s时，纵向诱导入流梯度系数λ_1c动态响应计算结果的对比。在计算中，旋翼的拉力系数保持不变，气动俯仰力矩保持为0。

由图 4(a~c)可知，本文模型的计算结果与取不同K_Re大小的动态入流模型的稳态响应幅值基本一致。时间常数也基本一致；纵向诱导入流系数变化相对旋翼俯仰角运动存在延迟特性，动态入流理论认为该特性是由空气附加惯性造成，并引入一阶惯性环节来描述。对于本算例的3个俯仰角速率，俯仰角速率越大，所对应的K_Re值越大，在1~1.5之间，与相关理论推导结果较为接近。图 4(d)给出了第1 s突增(抬头角速度30 °/s)至第3 s，拉力系数C_T分别等于0.004,0.006,0.008时的纵向诱导入流梯度系数λ_1c的动态正响应历程。由图 4(d)可知，在相同俯仰角速率下，K_Re值随拉力系数增加而增加，拉力系数越大，涡环尾迹模型计算的λ_1c动态响应变化规律越平滑，其变化曲线越接近于动态入流模型，当拉力系数较小时，动态入流模型用一阶惯性环节来表示尾迹畸变的动态过程，并不一定准确。图 5给出了悬停状态，C_T=0.006，旋翼在第1 s突增俯仰角速度(抬头速率30 °/s)至第3 s时的尾迹几何形状动态响应。由图可以看出，尾迹存在明显的拉伸、弯曲等畸变。

2.4 计算效率分析

以悬停状态拉力系数突增为例，在装有Intel Core(TM) i7-2600,3.4 GB处理器、4 GB内存的个人计算机上，对涡环尾迹模型进行时间步进计算，其中步长约为0.052 36 s，而CPU计算每一时间步长平均耗时约为0.048 s。可见，涡环尾迹模型在机动飞行状态状态计算效率较高，基本可以满足实时仿真的要求。

3 结论

(1) 本文提出了一种新的旋翼动态尾迹模型，该模型采用6自由度涡环单元来表示旋翼尾迹涡系，可描述机动飞行中旋翼尾迹动态畸变行为，具有应用于实时飞行仿真的潜力。通过与其他理论模型的对比分析，验证了本模型的正确性。

(2) 通过与尾迹畸变增广的动态入流模型对比分析表明，旋翼尾迹在俯仰角速度突增时存在明显的拉伸、弯曲等畸变，拉力系数与俯仰突增角速度大小均会对尾迹畸变参数K_Re产生影响，对于动态入流模型不同状态应取不同 K_Re值，而本文不需要人为设置该参数，具有更好的通用性。