基于黏性涡粒子/沙粒DEM的直升机沙盲建模

谭剑锋，何龙，于领军，周国臣

1. 南京工业大学 机械与动力工程学院，南京 211816 2. 中国空气动力研究与发展中心 低速空气动力研究所 旋翼空气动力学重点实验室，绵阳 621000 3. 陆军航空兵学院 航空机械工程系，北京 101123 4. 南京航空航天大学 航空学院， 南京 210016

直升机在沙漠、戈壁等环境下垂直起降和近地飞行，旋翼高速旋转产生的下洗流将掀起细沙颗粒和碎石，形成显著的“沙盲(Brownout)”现象，从而阻挡飞行员视线，导致直升机起降、贴地飞行困难，甚至因无法识别周围环境(如障碍物等)导致飞行事故。美国国防部航空安全技术报告表明：“沙盲”导致的弱视环境(Degraded Visual Environments, DVEs)引发了几乎一半的空军直升机坠机事故，同时也是美国陆军直升机损伤的主要因素。因此，研究直升机“沙盲”现象以及建立直升机“沙盲”现象分析方法显得尤为重要。

为研究直升机“沙盲”现象的形成过程， Whitehouse等基于模型旋翼开展“沙盲”试验以测试“沙盲”状态下的可视度，然而试验并未测量沙粒运动轨迹和流场特性。随后Tritschler等采用高速PIV技术研究旋翼流场干扰两相流状态下的沙粒流动，结果表明旋翼桨尖涡的相互融合显著增加了靠近地面的上洗速度，并影响沙粒聚集和扬起高度。为填补全尺寸直升机“沙盲”现象试验数据空白，美国陆军联合研究项目实验室的Wong等在美国陆军Yuma试验场开展了EH-60L的3D-LZ“沙盲”飞行试验，获得了全尺寸直升机“沙盲”形成过程数据。Leishman等采用单相和双相流测试技术研究旋翼“沙盲”状态下的沙粒起跳问题，结果表明受地面影响，旋翼桨尖涡发生拉伸、融合，并形成高度三维非定常特性，诱发地面涡的形成，从而导致沙粒起跳、悬浮。为进一步明确旋翼涡系的影响，2017年，Leishman等采用高速流体显示技术(Flow Visualization，FV)和高速两分量PIV研究旋翼地效下的旋翼涡系结构，研究表明受地面影响，旋翼桨叶内段涡系与桨尖涡存在速度差异，并在靠近地面融合，形成复杂的湍流结构，影响地面射流特性，最终影响“沙盲”形成。上述试验表明直升机“沙盲”现象受旋翼与地面非定常气动干扰影响显著。

为研究直升机“沙盲”现象的形成过程，众多研究者开展直升机“沙盲”现象的数值分析方法研究。Quackenbush等采用基于自由尾迹模型的CHARM软件研究旋翼靠近地面流场，并研究直升机“沙盲”现象，结果表明直升机贴地飞行状态下地面涡引发沙粒移动。随后，Whitehouse等基于自由尾迹Free-wake和面元法建立了旋翼/地面气动干扰模型，研究旋翼靠近地面的流场特性，并耦合沙粒运动模型和可视化模型模拟直升机“沙盲”现象。结果表明，此方法可实现直升机“沙盲”现象模拟，但并未考虑旋翼桨尖涡的黏性效应和地面黏性效应，也未考虑沙粒的摩擦碰撞特性。为此，Andrea采用基于面元-时间步进自由尾迹模型的ADPANEL软件研究旋翼靠近地面的非定常流动现象，并耦合粒子输运模型(Particle Transport Model，PTM)研究EH101和纵列式直升机“沙盲”特性，但方法中仍未考虑黏性效应、地面复杂流动以及沙粒起跳真实特性。此外，Keller等基于涡量输运模型(Vorticity Transport Model，VTM)和无黏地面模型研究旋翼靠近地面的复杂流动现象，并耦合基于沙粒起跳通量和起跳速度罚值的半经验模型构建沙粒输运模型PTM，从而模拟直升机“沙盲”现象，但仍未考虑沙粒间以及与地面的碰撞。为此，Leishman等采用时间精确自由尾迹模型和地面镜像模型模拟旋翼近地面流场，耦合基于沙粒输运方程的沙粒动力学模型，建立拉格朗日沙粒跟踪方法。方法采用Bagnold沙粒扬起速度罚值模型和沙粒摩擦速度罚值模型，并嵌入Shao沙粒夹带和起跳概率模型，模拟直升机“沙盲”现象。结果表明此方法可实现“沙盲”现象模拟，且旋翼与地面干扰的非定常尾迹结构对沙云特性影响显著。然而，方法仍未考虑旋翼尾迹和地面的无黏效应，且采用基于半经验的沙粒夹带模型、沙粒起跳模型、以及速度罚值的沙粒通量模型计算沙粒扬起过程，而未考虑沙粒碰撞存在的能量交换。为此，2020年，徐国华等采用CCM++软件和离散动力学EDEM软件开展了全尺寸直升机“沙盲”现象初步研究，但研究的沙粒数量较少(万级)，沙粒间距较大，因此在“沙盲”形成过程和计算方法方面仍需进一步研究。

基于拉格朗日体系的黏性涡粒子方法(Viscous Vortex Particle Method, VVPM)能较好捕捉旋翼/地面干扰效应和复杂流场特性，本文将采用非定常面元/VVPM/黏性地面气动模型模拟旋翼时变尾迹结构和旋翼/地面非定常气动干扰效应，并基于离散动力学，建立沙粒DEM模型和沙粒-流场耦合模型，构建基于黏性涡粒子/沙粒DEM(Discrete Element Method)的直升机“沙盲”现象分析方法。通过与美国陆军Yuma试验场的EH-60L“沙盲”飞行试验测量结果和基于沙粒夹带模型、沙粒起跳模型、速度罚值沙粒通量模型的拉格朗日沙粒跟踪方法计算结果对比，验证本文方法计算“沙盲”的准确性。随后，研究直升机悬停和前飞状态下的“沙盲”形成过程和特性。

1 计算方法

1.1 旋翼/地面非定常气动干扰模型

直升机“沙盲”状态的旋翼/地面非定常干扰流场受旋翼尾迹和旋翼气动载荷影响显著。非定常面元法能较好模拟旋翼悬停和前飞状态的桨叶气动载荷，且相比于CFD具有较高的计算效率，已应用于单旋翼、双旋翼气动干扰分析。为此，本文采用非定常面元法模拟直升机旋翼桨叶气动特性，直升机旋翼流场势函数可表示为

(1)

式中：为桨叶表面汇；为桨叶表面偶极子；为桨叶表面；为尾迹涡面；为表面法矢量；为空间点位置矢量。

直升机旋翼桨叶边界条件表示为

(2)

(3)

根据Neumann边界条件和后缘Kutta条件将式(2)转为代数方程，计算得到旋翼桨叶汇和偶极子分布。根据非定常伯努利方程计算得到旋翼桨叶压力系数分布：

(4)

式中：为参考压强；为空气密度；为当地速度；为当地压强；为当地参考速度。

根据三维Navier-Stokes方程，采用黏性涡粒子法模拟旋翼尾迹涡量分布：

(5)

式中：为运动黏性系数；为旋翼尾迹涡矢量分布；为旋翼速度场。

为体现地面对旋翼尾迹结构和流场的干扰，采用基于第二类Fredholm方程的黏性地面气动模型：

(6)

式中：为地面切矢量；为直升机旋翼桨叶和尾迹诱导的滑移速度，即

(7)

采用四边形单元离散地面，根据式(6)黏性地面气动模型和离散非均匀Neumann问题的Green积分函数，耦合涡扩散方程计算得到地面涡量扩散到流体的涡量为

(8)

式中：为三维热扩散核函数；为积分面积。

耦合式(1)、式(3)、式(5)～式(7)实现旋翼桨叶气动模型、尾迹模型与黏性地面气动模型的相互诱导耦合，耦合式(5)和式(8)实现旋翼尾迹和黏性地面气动模型的涡耦合，形成旋翼/地面非定常气动干扰模型，实现对直升机“沙盲”状态下的旋翼/地面非定常干扰流场求解。

1.2 沙粒DEM模型

“沙盲”的形成过程是颗粒相为稀相两相流动过程，在拉格朗日体系下，追踪颗粒相的运动，从而得到沙云的形成特征。对于沙粒间的相互作用，离散单元法(Discrete Element Model)在离散度大、流动显著的颗粒流系统分析中得到广泛应用。基于离散动力学理论，沙粒运动方程由牛顿第二定律得出,并考虑沙粒受到的流体力、碰撞力、重力以及沙粒旋转：

(9)

(10)

式中：为沙粒颗粒质量；为沙粒速度矢量；为沙粒惯量；为沙粒角速度矢量；为沙粒力矩矢量。通过采用显式积分技术，计算沙粒的运动轨迹，得到沙粒在空间的分布形态。

(11)

(12)

(13)

式中：为滚动摩擦系数；为接触点到沙粒质量中心的距离；为沙粒角速度单位矢量。

图1 沙粒接触碰撞模型Fig.1 Contact model of sand particle

1.3 沙粒-流场耦合模型

基于沙粒-流场耦合模型，将沙粒DEM动力学模型嵌入直升机旋翼/地面非定常干扰流场，实现直升机“沙盲”现象分析。由于直升机“沙盲”的沙粒直径为100 μm量级，因此忽略沙粒对旋翼/地面非定常干扰流场的影响，但需体现旋翼/地面非定常干扰流场对沙粒动力学特性的影响。为此，将流场作用于沙粒的气动力耦合至沙粒动力学方程，从而实现沙粒-流场耦合。

沙粒与流场的气-固两相间作用主要通过相间的作用力来耦合，沙粒受到流体力为

(14)

式中：为空气密度；为流场速度；为沙粒移动速度；为沙粒等效面积。

根据低速流体的Stokes定律，计算得到沙粒在空气中的阻力系数为

(15)

式中：沙粒雷诺数为

(16)

其中：为沙粒等效直径。

旋翼/地面干扰流场与沙粒DEM耦合计算流程如图2所示。图中：为旋翼/地面干扰流场计算时间；Δ为计算步长；为计算终止时间；Δ为沙粒DEM计算时间步长：

基于上述沙粒-流场耦合模型，将直升机旋翼/地面非定常气动干扰分析方法与沙粒DEM模型耦合，建立直升机“沙盲”现象数值分析方法。

图2 直升机“沙盲”计算流程图Fig.2 Computer flow diagram of helicopter brownout

2 直升机“沙盲”特性

2.1 旋翼/地面干扰流场验证

为了验证本文直升机“沙盲”现象数值分析方法计算旋翼/地面干扰流场的准确性，采用美国马里兰大学旋翼/地面干扰试验，模型旋翼的桨叶数目、半径、扭转角、弦长、桨尖马赫数、翼型分别为2、0.086 m、0°、0.019 6 m、0.08、圆弧翼型。旋翼总距和离地高度分别为12°、1.5。

图3给出了旋翼不同径向位置处的径向速度随高度变化，以及试验结果和CFD的OVERTUNS软件计算结果(其中为径向速度，为垂向速度)。各径向位置处，计算得到的速度随高度变化趋势与试验测量结果和CFD计算结果基本吻合。随着高度增加，径向速度先快速增加后减小，形成典型的地面射流现象。此外，随着径向位置增加，径向速度峰值对应的高度逐渐减小，速度峰值先增加后略有减小，且速度峰值和其对应的高度与试验结果均吻合较好。因此，本文直升机“沙盲”现象数值分析方法较好模拟旋翼/地面干扰下的流场特性。

图3 旋翼/地面干扰下的径向速度变化特性Fig.3 Varations of radial velocities during rotor-ground interaction

2.2 前飞状态EH-60L“沙盲”验证

为验证本文直升机“沙盲”现象数值分析方法的准确性，采用Wong和Tanner在美国陆军 Yuma试验场YPG的Oasis着陆区EH-60L直升机“沙盲”飞行试验。飞行测试试验区如图4所示，着陆区包括7条宽度不一的跑道，其中飞行试验在跑道1、3、5、7中进行，其他跑道设置了包括电线杆、电缆、小型建筑、汽车等障碍物，以便测量EH-60L机载传感器灵敏度。试验通过由多台相机构成的立体拍摄系统记录EH-60L飞行时“沙盲”照片，并获得“沙云”位置和形态。

图4 “沙盲”试验的Oasis着陆区示意图Fig.4 Schematic of Oasis landing zone used for brownout experiments

根据YPG试验场Oasis着陆区沙粒样本分析获得沙粒主要分布在0～1 000 μm。由于大直径沙粒重量大，移动和空气中悬浮困难，而占主要成分的小直径沙粒重量轻，移动和悬浮容易，因此，本文分析中建立长度为160 m、宽度为160 m的地面模型，采用球形沙粒，直径为1～100 μm，且沙粒直径和位置随机分布。

EH-60L直升机旋翼由4片后掠20°桨叶构成，桨叶半径、弦长、翼型、后掠角分别为8.18 m、0.527 4 m、SC1095和SC1095R8、20°，桨叶采用非线性负扭,如表1所示。EH-60L直升机平均前飞速度为37 km/h，飞行重量为7 394 kg，飞行轨迹如图5和图6所示。

表1 EH-60L直升机参数Table 1 Parameters of EH-60L helicopter

图5 各时间点直升机空间位置Fig.5 Positions of helicopter at different time points

图6 直升机飞行轨迹和Oasis着陆区模型Fig.6 Trajectory of helicopter and model of Oasis landing zone

计算得到此飞行状态下各时间点直升机沙云位置、飞行测试结果、以及基于拉格朗日沙粒跟踪方法的计算结果如图7～图9所示。

图7 沙云形态t=0 sFig.7 Shape of dust clouds at t=0 s

本文方法计算得到各时间“沙云”俯视图和侧视图位置均与飞行测试的“沙云”结果比较吻合，且计算获得的“沙云”轮廓与飞行测试轮廓一致。由此表明本文方法能较好的模拟直升机近地飞行的“沙云”轮廓和形成过程。

在=0,2 s，本文方法和拉格朗日沙粒跟踪方法计算得到“沙云”前端位置基本相同，均略超前于飞行测试结果(见图8(a))。=4 s时，拉格朗日沙粒跟踪方法计算得到“沙云”前端位置明显超前于飞行测试结果，而本文方法计算得到的“沙云”前端位置与飞行测试结果一致(见图9(a))。此外，拉格朗日沙粒跟踪方法计算得到的“沙云”向后收缩，与飞行测试结果的“沙云”扩张形式存在明显差异(见图9(a))，而本文方法计算得到的“沙云”向两侧扩散，与飞行测试结果一致。因此，相比于拉格朗日沙粒跟踪方法，本文方法能更好捕捉“沙云”前端位置和形状变化特性。

图8 沙云形态t=2 sFig.8 Shape of dust clouds at t=2 s

图9 沙云形态t=4 sFig.9 Shape of dust clouds at t=4 s

在=0,2,4 s，拉格朗日沙粒跟踪方法计算的“沙云”高度低于飞行测试结果 (见图7(b)、图8(b)、图9(b))，尤其在=10～16，而本文计算方法计算的“沙云”高度与飞行测试结果更接近，“沙云”轮廓与飞行测试结果基本一致。此外，随着“沙云”向后移动，拉格朗日沙粒跟踪方法计算得到的“沙云”高度先增加，后快速降低，由此表明直升机下方的沙粒扬起，后快速下沉，而本文方法计算得到的沙粒先扬起，而后悬浮，并与飞行测试的沙云形态特性一致，且“沙云”高度变化特性与飞行测试结果更吻合。

因此，相比基于自由尾迹模型和沙粒起跳罚值等模型的拉格朗日沙粒跟踪方法，本文直升机“沙盲”数值分析方法计算得到的直升机“沙云”现象精度更高，与美国陆军YPG试验场Oasis着陆区EH-60L“沙盲”飞行试验结果更吻合。

2.3 EH-60L直升机悬停“沙盲”形成特性

基于上述EH-60L直升机参数和地面参数，计算模拟EH-60L悬停状态(=1.5)“沙盲”形成过程，结果如图10所示。

=0.5 s时旋翼/地面干扰形成的地面射流对沙粒移动影响时间较短，暂未形成比较明显的“沙云”现象。受到旋翼桨尖涡与地面干扰的非定常流场作用(见图11，其中为垂直方向速度)，沙粒向外侧移动，在旋翼桨尖涡处形成明显的空白区，且在旋翼桨尖涡之间形成堆积(见图10(a))，主要原因在于旋翼桨尖涡正下方速度较大，沙粒移动速度较快，而桨尖涡之间的沙粒受到流场作用起跳，并受沙粒间碰撞作用而停止移动形成堆积。

图10 悬停状态EH-60L“沙盲”形成过程Fig.10 Process of brownout for EH-60L in hover

=1.0 s时，旋翼地面射流对沙粒的作用时间增长，靠近旋翼的沙粒堆积仍然存在，而受到最外侧旋翼桨尖涡形成的地面射流作用，沙粒开始起跳、扬起。时间增加到2.0 s，沙粒起跳高度逐渐增加，扬起现象更加明显。时间=4.0 s时，受地面射流和桨尖涡作用，细小颗粒扬起高度进一步增加，并沿着旋翼速度场移动。受重力和流体力的作用，沙粒向旋翼上方和内侧移动，并逐步形成沙云。时间增加到7.0 s后，沙粒起跳更加明显，扬起的沙粒形成非常明显的沙云，并逐步包围EH-60L直升机，从而形成典型的“沙盲”现象。

图11 悬停状态EH-60L/地面干扰流场Fig.11 Flowfield of EH-60L/ground interaction in hover

2.4 EH-60L直升机前飞状态“沙盲”形成特性

基于上述EH-60L直升机参数和地面参数，计算模拟EH-60L前飞状态(=1.0)沙盲形成过程，结果如图12所示。

直升机旋翼桨尖涡与地面干扰产生的两个卷起桨尖涡，推动直升机两侧沙粒向外移动，而中间的沙粒向后移动，前部沙粒向前移动(见图12(a)和图12(b))。时间=0.5 s，直升机两侧和前方的流场作用时间较短，沙粒暂未扬起。时间增加到1.0 s，受到直升机两侧地面卷起桨尖涡的作用，沙粒向外移动速度更明显，且受卷起桨尖涡向上诱导作用(见图13(a)，其中为平面合速度)，沙粒出现起跳、碰撞、轻微扬起。时间增加到2.0 s时，直升机两侧地面卷起桨尖涡所产生的沙粒起跳、扬起现象更为明显。同时，受直升机前方卷起桨尖涡诱导的向前地面射流作用(见图13(b)，其中为平面合速度)，机头前部的沙粒向前移动和起跳。当时间增加到5.0 s，受两侧卷起桨尖涡的显著影响，沙粒扬起现象显著，形成两侧沙云现象，直升机前部的沙粒扬起现象也更为明显。当时间增加7.0～10.0 s，直升机两侧的沙粒扬起速度更快，在重力和流体力作用下，沙粒悬浮并逐步回落，形成典型的回流现象。受直升机前方卷起桨尖涡作用，直升机机头前部的沙粒扬起较为严重，出现悬浮的沙云现象，最终形成典型前飞“沙盲”现象。

图13 前飞状态EH-60L/地面干扰流场Fig.13 Flowfield of EH-60L/ground interaction in forward flight

3 结 论

耦合旋翼黏性涡粒子和黏性地面气动模型，建立基于离散动力学的沙粒DEM模型和沙粒-流场耦合模型，构建直升机“沙盲”现象分析方法，并应用于美国陆军Yuma试验场EH-60L“沙盲”计算，得到如下结论：

1) 建立的直升机“沙盲”现象分析方法较准确模拟EH-60L全尺寸直升机近地飞行状态下的“沙云”轮廓。

2) 本文方法计算得到各时间的“沙云”形状、扬起位置、高度与飞行测试结果基本一致，并能较准确体现各飞行状态下的“沙云”形成过程。

3) 相比基于自由尾迹模型和沙粒起跳罚值等模型的拉格朗日沙粒跟踪方法，本文建立的直升机“沙盲”现象分析方法计算得到的“沙云”轮廓精度更高，与飞行测试结果更吻合。

4) 悬停状态下的旋翼桨尖涡与地面干扰形成地面射流，推动沙粒移动堆积，并导致沙粒碰撞起跳，在旋翼/地面干扰流场作用下扬起，形成沙云，覆盖直升机前方视野，形成沙盲现象。

5) 前飞状态下的旋翼桨尖涡与地面干扰形成两侧扩散卷起桨尖涡和机头前方卷起涡，推动沙粒移动、碰撞起跳，并在桨尖涡诱导作用下扬起，覆盖直升机前方视野，形成“沙盲”现象。

后续研究可将旋翼CFD方法和“沙盲”求解方法耦合，更好体现现代旋翼桨叶先进翼型对“沙盲”现象的影响，并可进一步研究旋翼构型与气动设计对“沙盲”现象的影响规律。