基于DEM的直升机沙盲加速计算方法.pdf

http:/DOI:10.13700/j.bh.1001-5965.2021.0450基于 DEM 的直升机沙盲加速计算方法谭剑锋1，*，韩水1，王畅2，于领军3(1.南京工业大学机械与动力工程学院，南京211816；2.中国空气动力研究与发展中心旋翼空气动力学重点实验室，绵阳621000；3.陆军航空兵学院航空机械工程系，北京101123)摘要：直升机沙盲数值模拟是研究沙盲演化特性的重要手段，而沙盲由众多动力学特性复杂的沙粒构成，这导致沙盲数值模拟复杂且计算量庞大。基于离散单元法(DEM)和沙粒动力学方程，将沙粒映射至背景网格实现加速计算，并将背景网格分裂为多子区再次加速计算，构建背景网格映射-分裂加速计算模型，且耦合沙粒接触碰撞模型、沙粒-流场耦合模型、旋翼/地面气动干扰模型，提出基于 DEM 的直升机沙盲加速计算方法。通过与美国陆军 EH-60L 着陆-起飞沙盲测试结果对比表明：所提方法能准确捕捉着陆-起飞状态的直升机沙盲，且相比于沙盲直接模拟方法，所提方法计算量显著减小。直接模拟方法的计算量随沙粒数量抛物线增加，而所提方法计算量随沙粒数量线性增加。当沙粒数量大于 1107时，相比于仅背景网格映射模型加速方法，所提方法计算量减小70.29%。关键词：沙盲；离散单元法；映射-分裂模型；加速方法；直升机中图分类号：V221.52；TB553文献标志码：A文章编号：1001-5965（2023）06-1352-10受地面干扰，起飞或着陆状态下的直升机旋翼产生较大下洗流和外洗流，推动沙粒移动和扬起，形成直升机沙盲现象。直升机沙盲不仅引发弱视环境1，威胁飞行安全，且扬起的沙粒易进入发动机，与桨叶碰撞摩擦，导致直升机核心部件加速损坏，因此，直升机沙盲已成为影响飞行安全的重要因素。由于直升机沙盲涉及复杂的旋翼流场和沙粒运动，因此，开展直升机沙盲试验存在较大难度，目前国内外主要通过基于沙盲数值模拟方法研究直升机沙盲形成特性。文献 2-4 采用基于雷诺平均Navier-Stokes（ReynoldsaverageNavier-Stokes,RANS）方程的计算流体力学(computationalfluiddynamics,CFD)方法计算旋翼/地面干扰流场，并耦合基于连续流体假设的沙粒运动模型研究沙盲现象。然而，此方法并未考虑单个沙粒运动特性。随后，文献5-6 将时间精确 RANS 方法和拉格朗日沙粒跟踪方法耦合，研究沙盲形成过程，但方法中采用 Bagnold沙粒扬起速度罚值模型、Shao 沙粒夹带和起跳概率模型等经验模型，且此方法计算量较大。为此，文献 7 将自由涡流法-雷诺平均 Navier-Stokes 方程(freevortexmethod-ReynoldsaverageNavier-Stokes,FVM-RANS)混合流场计算方法与拉格朗日沙粒跟踪 方 法 耦 合，并 嵌 入 基 于 图 形 处 理 器（graphicsprocessingunit,GPU）的并行算法，以提高直升机沙盲数值模拟效率，但此方法整体计算量仍然较大。文献 8 耦合 CCM+流体软件和 EDEM 软件分析直升机沙盲特性，然而庞大的计算量限制了沙粒颗粒数量。由于直升机沙盲现象需计算旋翼/地面非定常收稿日期：2021-08-09；录用日期：2021-08-27；网络出版时间：2021-09-1611：28网络出版地址： J.北京航空航天大学学报，2023，49（6）：1352-1361.TAN J F，HAN S，WANG C，et al.Accelerated computational method of helicopter brownout based on DEMJ.Journal of BeijingUniversity of Aeronautics and Astronautics，2023，49（6）：1352-1361（in Chinese）.2023年6月北京航空航天大学学报June2023第49卷第6期JournalofBeijingUniversityofAeronauticsandAstronauticsVol.49No.6干扰流场和靠近地面的旋翼尾迹结构，而 CFD 方法存在计算量大、尾迹数值耗散偏大的问题。为此，国外众多研究者采用如自由尾迹等涡方法计算旋翼/地面非定常干扰流场，以提高沙盲数值模拟效率。文献 9-11 将自由尾迹模型与沙粒输运模型耦合，研究 EH101 直升机沙盲现象。文献 12 将涡量输运模型和沙粒输运模型耦合，研究纵列式直升机沙盲特性。此外，基于 Bagnold 和 Shao 等半经验沙粒假设模型，文献 13-14 将时间精确自由尾迹模型、地面镜像模型、拉格朗日沙粒跟踪方法耦合研究沙盲两相流特性。虽然相比于基于 RANS 的CFD 方法，涡方法计算旋翼/地面干扰流场效率更高，然而庞大的沙粒动力学计算量是影响直升机沙盲数值模拟效率的重要因素。为此，文献 15 采用包括 Gaussian 方法、k 均值方法、Osiptsovs 方法等聚类方法提高直升机沙盲模拟效率。此类方法将局部沙尘颗粒集聚，并当量为一个等效沙粒簇，从而减小沙粒动力学方程求解数量，降低计算量。文献 16 将聚类方法与快速多极子方法(fastmultipolemethod,FMM)和 GPU 并行算法耦合，以提高沙盲数值计算效率，但此类方法存在由沙粒集聚当量化导致的沙粒轨迹累积误差，且计算效率与计算精度存在反比关系。文献 17 采用基于半隐式欧拉方程的沙粒运动模型和预先计算的流场耦合，并嵌入消 息 传 递 接 口(messagepassinginterface,MPI)和GPU 并行技术，以提高沙盲数值模拟效率。然而，为方便实现并行和提高计算效率，此方法忽略沙粒接触碰撞。由于沙粒-沙粒、沙粒-地面接触碰撞将改变沙粒动力学特性和沙云特性，增加聚类方法的累积误差和沙盲数值模拟计算量，因而给聚类方法和基于 MPI 和 GPU的并行方法带来了较大挑战。为此，本文基于离散单元法(discreteelementmethod，DEM)和沙粒动力学方程，将沙粒映射至背景网格以实现加速计算，并将背景网格分裂为多子区以再次加速计算，构建背景网格映射-分裂加速计算模型。随后耦合沙粒接触碰撞模型、沙粒-流体耦合模型、旋翼/地面气动干扰模型，提出直升机沙盲加速计算方法。通过与美国陆军 EH-60L 着陆-起飞沙盲试验对比，验证本文方法的准确性，并通过与 DEM 直接模拟方法、背景网格加速模型对比，验证本文方法的计算效率。1计算方法1.1沙粒离散动力学模型直升机复杂的沙云形态由众多沙粒构成，因此，通过拉格朗日体系下的大规模沙粒运动轨迹分析研究直升机沙盲形成特性。为此，基于离散动力学理论18-19和 DEM，并根据牛顿第二定律，可得到沙粒动力学方程为mdVpdt=Fi+Ff+Fg（1）Idpdt=T（2）m VpIpTFfFiFg式中：、和分别为沙粒质量、速度、惯量和角速度；、和分别为沙粒力矩矢量、流场诱导的气动力、沙粒-沙粒和沙粒-地面间的碰撞力及重力。对式(1)和式(2)进行时间推进求解，得到各时间点下的单个沙粒运动特性，获得沙粒空间分布，从而可分析直升机沙盲形成过程。1.2沙粒接触碰撞模型Fni沙粒-沙粒接触碰撞和沙粒-地面接触碰撞显著改变沙粒的运动轨迹和运动特性，因此，在分析直升机沙盲现象和沙云形成过程中，需考虑沙粒接触碰撞特性。由于沙场环境下沙粒数量庞大，单颗沙粒将与其他沙粒碰撞，因此，针对每个沙粒，需计算与周围沙粒的接触碰撞力。基于Hertz-Mindlin 模型20-21，沙粒接触碰撞模型由弹簧-阻尼器-滑块系统构成，如图 1 所示。沙粒与沙粒的接触碰撞法向力为Fni=43Y3/2nRn25/6Snmvreln（3）YRmnnvrelnSn式中：、和 分别为沙粒的当量弹性模量、半径、质量、法线方向的变形量、单位矢量、相对法向速度、法向刚度和等效刚度。Fti沙粒与沙粒的接触碰撞切向力为Fti=Sttt25/6Snmvrelt（4）法向碰撞力Fni碰撞力Fi碰撞力Fi切向碰撞力FtiziyixizjyjxjRnt摩擦系数沙粒j沙粒i切向阻尼切向刚度法向刚度法向阻尼图1沙粒碰撞的弹簧-阻尼器-滑块系统模型Fig.1Contactmodelofsandparticlewithspring-dashpotsystem第6期谭剑锋，等：基于 DEM 的直升机沙盲加速计算方法1353ttvreltSt式中：、和分别为沙粒 i 的切向变形量、单位切矢量、相对切向速度和切向刚度。r沙粒的滚动摩擦力为r=r?Fni?Rii（5）rRii式中：、和分别为摩擦系数、沙粒 i 中心与接触点的间距和角速度矢量。沙粒-地面接触碰撞与沙粒-沙粒接触碰撞相似，将地面等效为刚性，且考虑沙粒的弹性模量、剪切模量、质量等参数，采用上述接触碰撞模型计算沙粒-地面接触碰撞。1.3沙粒-流场耦合模型直升机沙盲是由沙粒和旋翼流场构成的两相流动，因此，直升机沙盲数值模拟需考虑沙粒-流场耦合作用。由于沙床沙粒直径通常为 100m 量级，微小尺度的沙粒对大尺度的直升机旋翼流场影响较小，可忽略沙粒对直升机旋翼流场的影响，但大尺度旋翼流场对微小尺度沙粒运动特性影响较大。为此，需考虑直升机旋翼/地面干扰流场对沙粒的耦合作用。基于单向耦合思想，将旋翼流场诱导的沙粒气动力嵌入沙粒动力学方程式(1)，实现沙粒-流场耦合模拟。Ff直升机旋翼/地面干扰作用下的沙粒气动力可由沙粒阻力特性表述为Ff=12air?vpvf?(vpvf)CDA（6）airvfvpA式中：和分别为流场密度和速度；和 分别为沙粒速度和面积。CD根据 RANS 方程定律，沙粒气动阻力系数为CD=24Re(1+0.15Re0.687)（7）Re式中：沙粒雷诺数为Re=?vpvf?de（8）de其中：为沙粒等效直径；为空气黏性系数。vf式(6)和式(8)的旋翼流场速度可由直升机旋翼/地面气动干扰模型计算得到：vf=wV13K(xsyi)idV+wSKG(xs z)(z)dS（9）wV13K(xsyi)idVrSKG(xs z)(z)dSyiKGzS式中：为旋翼尾流诱导项，xs为沙粒的位置矢量；为旋翼诱导项；为地面诱导项；为涡粒子对称光滑参数；K 为 Biot-Savart 核函数；和 i分别为涡粒子位置和涡矢量；为速度势函数；为 Green 函数的梯度；和 分别为地面涡面位置和涡矢量；为地面面积。根据直升机旋翼非定常面元气动模型22，式(9)中的速度势函数 可表示为(x,y,z,t)=14wSBn(1rt)dS 14wSB(1rt)dS+14wSwn(1rt)dS（10）SBSWnrt式中：和 分别为直升机旋翼桨叶的汇和偶极子分布；、和 分别为旋翼桨叶面积、涡面面积、法矢量和距离。直升机流场可由基于黏性 RANS 方程的黏性涡粒子方法23计算得到，因此，式(9)中的涡粒子位置和涡矢量为dydt=u(y,t)（11）ddt=u+2（12）u式中：为旋翼流场速度矢量。地面对旋翼流场的影响由黏性地面气动模型模拟24计算得到，因此，式(9)中的地面涡矢量为n2wSKG(z z)(z)dS=wV13K(zyi)idV+（13）1.4沙盲加速计算方法O(N2)采用直接模拟方法可计算式(1)、式(3)和式(4)的沙粒接触碰撞力和动力学方程，且易实现并行计算，但由于沙粒相互间均需计算接触碰撞，因此，整体沙粒动力学计算量为(N 为沙粒数量)。沙粒数量较小时，沙粒动力学计算量较小，然而由于直升机沙盲形成范围较大，沙粒数量庞大，因此，采用直接模拟方法将很难实现沙粒数量庞大的沙盲模拟。为此，针对 DEM 和直升机沙盲特性，首先建立背景网格映射模型，如图 2 所示。由于沙粒 i 仅与周围沙粒产生接触碰撞，与相对较远的沙粒不存在碰撞。为此，将计算域划分为若干的背景网格，并将沙粒 i 映射至背景网格(ig,jg,kg)。由此，沙粒i 与周围沙粒接触碰撞转化为背景网格(ig,jg,kg)的内部沙粒碰撞计算和背景网格(ig,jg,kg)与相邻背景网格(ig1ig+1,jg1jg+1,kg1,kg+1)的沙粒碰撞计算，其中相邻背景网格数量最多为 26。由此，基于背景网格计算式(1)、式(3)、式(4)等沙粒接触碰撞力，沙粒接触碰撞数量计算次数大大减小，直升机沙盲数值仿真计算效率得到改善。由于沙场沙粒直径通常为 100m 量级，背景网格(ig,jg,kg)与相邻背景网格的沙粒接触碰撞仅发生在较小的相邻边界区域，为此，建立背景网格分裂模型(见图 3)，对背景网格映射模型进行再次加1354北 京 航 空 航 天 大 学 学 报2023年速计算。采用二维背景网格以更清晰阐述背景网格分裂模型。将背景网格(ig,jg)分为 5 个区域，分别为 A、B、C、D、E 区，其中 A、B、C、D这 4 个区在背景网格(ig,jg)的 4 个边界，宽度为沙粒最大直径，而 E 区为背景网格(ig,jg)中央区域。因此，E 区的沙粒不与背景网格(ig,jg)的相邻网格沙粒接触碰撞，因而无需计算其接触力，仅需计算自身内部接触碰撞力。A、B、C、D 这 4 个区的沙粒将与相邻背景网格相应区域的沙粒碰撞，因而仅需计算局部区域的沙粒接触碰撞力。如背景网格(ig,jg)的 A 区沙粒将与背景网格(ig1,jg)的 C 区沙粒、背景网格(ig1,jg+1)的 D 区沙粒、背景网格(ig1,jg1)的 B 区沙粒、背景网格(ig,jg+1)的 D 区沙粒、背景网格(ig,jg1)的 B 区沙粒产生接触碰撞，而与其他区域无接触碰撞。因此，背景网格(ig,jg)的沙粒接触碰撞计算次数明显减小，从而减小沙粒接触碰撞计算量。耦合背景网格映射模型和背景网格分裂模型提出背景网格映射-分裂加速计算模型，实现沙盲的 DEM 加速计算。将上述背景网格映射-分裂加速计算模型、沙粒-流场耦合模型、沙粒接触碰撞模型、沙粒动力学模型、以及旋翼/地面气动干扰模型耦合，构建基于DEM 的直升机沙盲加速计算方法，实现直升机沙盲的数值模拟。2直升机沙盲加速计算特性分析2.1本文方法验证为验证直升机沙盲加速计算方法的准确性，计算 EH-60L 直升机着陆-起飞状态下的沙云形成过程，并与美国陆军 Yuma 沙盲试验25对比。试验着陆区如图 4 所示(为摄像机的可视角)，试验中采用立体拍摄系统测量直升机着陆过程各时间点的沙云形态。184.5 m相机1观测区域161.5 m110 m130.5 m240 m相机2跑道123456跑道773.2 mxy38.1 m30.5 m12.2 m45.7 m27.4 m61.0 mSN图4Yuma 试验着陆区Fig.4LandingzoneinYumaflighttest依据直升机参数和 Yuma 试验场地参数，建立直升机和地面数值模型，其中直升机质量为 7394kg，旋翼桨叶数目为 4，前倾角为 4，转速为 258r/min，桨叶半径为 8.18m，弦长为 0.5274m，后掠角为 20，翼型为 SC1095 和 SC1095R8，地面模型为 160m160m。根据实测的跑道沙粒分布区间，沙床由直径(1100m)和位置随机分布的球形沙粒构成。EH-60L 直升机着陆过程各时间点的位置和飞行轨迹如图 5 和图 6 所示，其中 x 为直升机纵向方向，z 为直升机高度方向，R 为旋翼的半径。计算得到时间为 2，6，17，25s 的沙云位置、试验结果25如图 7 和图 8 所示，图中同时给出拉格朗日沙粒跟踪方法计算结果13。图中红色球体表示沙粒背景网格沙粒碰撞沙粒碰撞沙粒碰撞(i)(ig,jg)图2背景网格映射模型Fig.2Backgroundgridmappingmodel碰撞(ig1,jg+1)(ig1,jg1)(ig1,jg)(ig+1,jg)(ig,jg)CE碰撞碰撞碰撞碰撞碰撞(ig,jg1)碰撞碰撞(ig,jg+1)(ig+1,jg1)背景网格(ig+1,jg+1)CD DA AB B图3背景网格分裂模型Fig.3Backgroundgridsplittingmodel10010203002461001020z/Rx/Rt/s飞行测试(z)本文(z)飞行测试(x)本文(x)图5着陆过程中 EH-60L 位置Fig.5EH-60Lspositionduringapproaching第6期谭剑锋，等：基于 DEM 的直升机沙盲加速计算方法1355飞行测试结果，黑色方块表示拉格朗日沙粒跟踪方法计算的结果，黄色为本文方法计算的结果。计算得到 EH-60L 直升机着陆-起飞状态下各时间点沙云俯视图和侧视图与飞行测试结果基本一致，且计算的沙云轮廓与飞行测量轮廓比较吻合，因此，本文方法能较好捕捉直升机着陆-起飞状态的沙云轮廓变化过程。t=2s 时，拉格朗日沙粒跟踪方法计算的沙云前端超前于飞行测试结果，而本文方法计算的沙云前端与试验测试结果基本吻合(见图 7(a)。此外，本文方法计算的沙云高度与拉格朗日沙粒跟踪方法计算结果基本相同(见图 8(a)，均略高于飞行测试结果。t=6s 时，本文方法计算得到沙云形态与拉格朗日沙粒跟踪方法计算结果和飞行测试结果基本相似，但 2 种方法计算得到的沙云前端位置均超前于飞行测试结果(见图 7(b)。此外，拉格朗日沙粒跟踪方法计算得到沙云高度明显低于飞行测试结果，而本文方法计算的沙云高度与飞行测试结果更加吻合(见图 8(b)。主要原因在于拉格朗日沙粒跟踪方法采用基于半经验模型的 Bagnold 沙粒扬起速度罚值模型、沙粒摩擦速度罚值模型、Shao 沙粒夹带和起跳概率模型模拟沙粒接触碰撞，而未考虑沙粒-沙粒和沙粒-地面的直接接触碰撞，从而影响沙粒碰撞、起跳、扬起和悬浮过程。t=17s 时，直升机飞行高度降至最低，接近着陆状态。相比于前述几个时间点，此时计算得到的沙云形态与飞行测试结果吻合的更好(见图 7(c)，沙云前端位置与飞行测试结果更接近。与前述时间相似，由于拉格朗日沙粒跟踪方法采用基于半经验的沙粒起跳等模型，计算得到的沙云前端超前于飞行测试结果，高度明显低于飞行测试结果，而本文方法直接模拟沙粒的接触碰撞，计算的沙云高度与飞行测试结果更吻合。xEH-60L160 m160 myz本文飞行轨迹飞行测试飞行方向着陆区域图6EH-60L 着陆轨迹Fig.6TrajectoryofEH-60Lduringapproaching(a)t=2 s(b)t=6 s(c)t=17 s(d)t=25 s飞行测试飞行测试飞行测试飞行测试拉格朗日沙粒跟踪方法拉格朗日沙粒跟踪方法拉格朗日沙粒跟踪方法拉格朗日沙粒跟踪方法本文方法本文方法本文方法本文方法直升机直升机xyzxyzxyzxyz图7不同时间的沙云俯视图Fig.7Dustcloudatdifferenttimefromtopview(a)t=2 s(b)t=6 s(c)t=17 s(d)t=25 s飞行测试飞行测试飞行测试飞行测试拉格朗日沙粒跟踪方法拉格朗日沙粒跟踪方法拉格朗日沙粒跟踪方法拉格朗日沙粒跟踪方法本文方法本文方法本文方法本文方法xyzxyzxyzxyz图8不同时间的沙云侧视图Fig.8Dustcloudatdifferenttimefromsideview1356北 京 航 空 航 天 大 学 学 报2023年t=25s 时，直升机起飞后快速爬升，拉格朗日沙粒跟踪方法计算得到的沙云前端位置明显滞后于飞行测试结果，而本文方法计算的沙云前端位置与飞行测试结果吻合较好(见图 7(d)。此外，拉格朗日沙粒跟踪方法计算的沙云高度低于飞行测试结果，而本文方法计算的沙云高度与飞行测试结果更加吻合。同时，随着沙云向后方移动，拉格朗日沙粒跟踪方法计算的沙云高度降低，且明显低于飞行测试结果，而本文方法计算的沙云高度保持较好，与飞行测试结果更吻合。主要原因为沙粒接触碰撞不仅发生在沙粒与地面，也发生在空间运动的沙粒，因而沙粒接触碰撞影响沙粒的移动、跳跃、轰击、扬起和悬浮等特性。为更准确分析本文方法的计算精度，根据相对误差 RMS 分析计算结果与飞行测试结果误差。RMS=vuut1Ntext1(Ntexti=1(yiyi,text)2)（14）式中：Ntext、yi和yi,text分别为飞行测试数据总数、仿真值和飞行测试值。根据飞行测试的沙云位置、拉格朗日沙粒跟踪方法计算的沙云位置及本文方法计算的沙云位置，结合式(14)，得到相对误差如表 1 所示。相比于拉格朗日沙粒跟踪方法，本文方法的 RMS 均更小。主要原因在于拉格朗日跟踪沙粒方法计算的沙云前端和高度均略偏离飞行测试结果，而本文方法可较好捕捉沙云轮廓。此外，相比于俯视图的沙云轮廓，本文方法计算的沙云高度误差减小的更明显。由此表明，本文方法能较好捕捉着陆-起飞状态的沙云形状。表1沙云计算相对误差Table1RMSsofcomputationaldustclouds时间/s方法RMS相对误差/%俯视侧视俯视侧视2拉格朗日沙粒跟踪方法0.34380.203本文方法0.29590.207210213.91289.8196拉格朗日沙粒跟踪方法0.23750.5621本文方法0.155780.22734.43359.53917拉格朗日沙粒跟踪方法0.60810.944本文方法0.13490.810577.79814.18625拉格朗日沙粒跟踪方法0.67221.1513本文方法0.47170.926429.82619.5382.2本文方法效率基于上述飞行状态，采用 OpenMP 并行和 Inteli7-9700CPU 测试本文方法效率。沙云计算时间步长为 1.0105s，计算时长为 0.001s。分别采用直接模拟方法、基于背景网格映射模型、基于背景网格映射-分裂模型计算沙云。由于本文方法的背景网格映射-分裂模型未涉及假设，因而计算精度未受影响。直接模拟方法、基于背景网格映射模型和基于背景网格映射-分裂模型的模拟计算时间如表 2 所示。相比于直接模拟方法，基于背景网格映射模型和基于背景网格映射-分裂模型的计算时间均显著减小。原因为 2 种方法仅需计算相邻背景网格沙粒接触碰撞，而无需计算所有沙粒接触碰撞，因而计算量明显减小。当沙粒数量小于 1106时，基于背景网格映射模型和基于背景网格映射-分裂模型的计算时间基本相似，原因为背景网格边界 4 个分区的沙粒接触碰撞计算所占比例较小，对沙盲模拟计算量影响较小。但当沙粒数量超过 1106时，基于背景网格映射-分裂模型的计算时间明显小于基于背景网格映射模型，且沙粒数量超过 1107时，计算时间减小 70.29%。图 9 为直接模拟方法、基于背景网格映射模型、基于背景网格映射-分裂模型的沙盲模拟计算时间随沙粒数量变化曲线。直接模拟方法的计算时间随沙粒数量呈现抛物线增加，主要原因为直接表2加速方法计算时间Table2Computationaltimeofacceleratedmethods沙粒数量计算时间/s直接模拟法基于背景网格映射模型基于背景网格映射-分裂模型1104246.7432.2234.2251044055.0146.6042.97110514304.82147.69139.225105381144.40398.08333.9111061493581.00795.13701.25510675405520.0010243.893857.501107117761400.0015521.874610.95第6期谭剑锋，等：基于 DEM 的直升机沙盲加速计算方法1357O(N2)模拟计算量为。然而基于背景网格映射模型和基于背景网格映射-分裂模型的计算时间随沙粒数量呈线性增加，主要原因为 2 种方法仅需计算相邻背景网格的沙粒接触碰撞。此外，相比于基于背景网格映射模型，基于背景网格映射-分裂模型的沙盲数值模拟计算时间随沙粒数量变化斜率明显更小，因而基于背景网格映射-分裂模型的沙盲数值模拟加速计算方法可有效减小直升机沙盲模拟时间。直接模拟法基于背景网格映射模型基于背景网格映射-分裂模型基于背景网格映射模型基于背景网格映射-分裂模型12840时间/107 s沙粒数量/104(a)3种加速计算方法15101001 00050050沙粒数量/10415101001 00050050(b)2种加速计算方法181260时间/103 s图9不同方法的计算时间随沙粒数量变化Fig.9Variationsofcomputationaltimeofdifferentmethodswithnumberofparticles2.3EH-60L 直升机着陆-起飞的沙盲演化特性直升机着陆-起飞过程中不同时间点的沙云和尾迹结构如图 10 和图 11 所示。着陆状态下，t=2s 时，旋翼桨尖涡向后移动并靠近地面(见图 10 中蓝色所示)，受到地面的干扰，旋翼桨尖涡向两侧扩散，并形成两侧卷起桨尖涡，从而诱发两侧外洗流和上洗流，推动直升机两侧沙粒移动、碰撞、跳跃、扬起，并在上洗流作用下悬浮，以及在卷起桨尖涡诱导下形成回流和悬浮，最终形成两侧沙云。直升机旋翼前端桨尖涡向下移动，受地面干扰作用，在机头下方形成弧形地面涡，诱发地面射流(见图 12，r 为距离旋翼中心径向距离，R 为旋翼半径)，从而推动沙粒向前移动、碰撞、跳跃、扬起，并在旋翼下洗流作用下形成回流和悬浮，最终构成直升机前端沙云。受飞行高度和飞行速度影响，直升机前端沙云覆盖旋翼中后部区域。t=6s 时，飞行高度逐渐降低，旋翼桨尖涡与地面干扰更加明显。旋翼前端桨尖涡与地面干扰形成更为明显的弧形地面涡，且相比于 2s 状态，位置向前移动，因而形成位置更靠前的地面射流，从而(a)t=2 s(b)t=6 s(c)t=17 s(d)t=25 s沙云沙云沙云沙云回流回流回流旋翼桨尖涡旋翼桨尖涡旋翼桨尖涡旋翼桨尖涡xyzxyzxyzxyz图10着陆状态 EH-60L 沙云俯视图Fig.10DustcloudforEH-60Lduringapproachingfromtopview(a)t=2 s(b)t=6 s(c)t=17 s(d)t=25 s沙云沙云沙云沙云回流回流扬起地面地面地面地面xyzxyzxyzxyz直升机直升机直升机直升机图11着陆状态 EH-60L 沙云侧视图Fig.11DustcloudforEH-60Lduringapproachingfromsideview1358北 京 航 空 航 天 大 学 学 报2023年推动沙粒形成位置更靠前的沙云，覆盖直升机旋翼前部，也覆盖直升机驾驶舱区域，形成弱视环境。t=17s 时，直升机飞行高度降至最低。相比上述几个时间点，旋翼桨尖涡与地面干扰更加明显，形成的弧形地面涡向两侧扩散区域更大，向机头前方移动距离更远，从而在直升机机头前方形成更为显著的地面射流，且射流速度显著增加(见图 12)。因此，沙粒移动速度更大、碰撞更为突出、扬起和回流更为严重，沙粒悬浮数量更多，形成浓度更大的沙云。沙云位置向机头前方移动，完全覆盖直升机机头和座舱，从而构成更为严重的弱视环境。随着时间继续增加，直升机开始爬升，旋翼桨尖涡与地面干扰形成的机头前方弧形地面涡减弱，且两侧扩散区域也同时减小。当 t=23s 时，旋翼与地面干扰区域向机身后方移动，从而显著减小机头前方地面径向速度，地面射流逐步消逝，机头前方沙云逐步消除。部分扬起的沙粒向机身后方移动，但未覆盖直升机前半部分和座舱，因而弱视环境逐渐消除。3结论1）计算得到 EH-60L 直升机着陆-起飞状态不同时间点沙云前端位置和高度均与飞行测试结果吻合较好，本文方法较好捕捉着陆-起飞状态的沙云形态。2）基于直接模拟的直升机沙盲数值模拟计算量随沙粒数量呈现抛物线增加，而基于背景网格映射-分裂模型的直升机沙盲数值模拟加速方法计算量随沙粒数量呈现线性增加，计算量显著减小。3）相比基于背景网格映射模型，基于背景网格映射-分裂模型的直升机沙盲模拟计算量随沙粒数量增加的斜率明显减小，且采用 OpenMP 并行和时间步长 1105s，沙粒数量超过 1107的着陆状态，计算量减小 70.29%。4）旋翼桨尖涡与地面干扰在机头前方形成弧形地面涡，诱发地面射流。随着直升机高度降低，机头前方的地面射流更加显著，诱发沙粒移动、碰撞、跳跃、扬起，并在上洗流和下洗流作用下回流和悬浮，形成机头前方弧形沙云，覆盖直升机机头和座舱，形成弱视环境。5）随着直升机爬升，飞行高度增加，旋翼桨尖涡与地面干扰减弱，机头前方地面射流削弱，沙云逐步减弱，并向机身后方移动，弱视环境消除。参考文献（References）SZOBOSZLAYZ,DAVISB,FUJIZAWAB.Degradedvisualen-vironmentmitigation(DVE-M)program,Yuma2016flighttrialsinbrownoutC/PresentedattheAHSInternational73rdAnnualFor-um&TechnologyDisplay.Alexandria:AHS,2017.1RYERSONC,HAEHNELR,KOENIGG,etal.Visibilityenhance-ment in rotorwash cloudsC/43rd AIAA Aerospace SciencesMeetingandExhibit.Reston:AIAA,2005.2WENRENYH,WALTERJ,FANM,etal.Vorticityconfinementand advanced rendering to compute and visualize complexflowsC/44thAIAAAerospaceSciencesMeetingandExhibit.Re-ston:AIAA,2006.3HAEHNELRB,MOULTONMA,WENRENW,etal.Amodeltosimulate rotorcraft-induced brownoutC/Proceedings of the 64thAnnual Forum of the American Helicopter Society.Alexandria:AHS,2008:589-601.4THOMASS,LAKSHMINARAYANVK,KALRATS,etal.Eu-lerian-LagrangiananalysisofcloudevolutionusingCFDcoupledwith a sediment tracking algorithmC/The American HelicopterSociety67thAnnualForum.Alexandria:AHS,2011:298-315.5KUTZBM,GUNTHERT,RUMPFA,etal.Numericalexamina-tionofamodelrotorinbrownoutconditionsC/Presentedatthe70thAnnualForumoftheAmericanHelicopterSociety.Alexan-dria:AHS,2014:2450-2461.6THOMASS,AMIRAUSM,BAEDERJ.GPU-acceleratedFVM-RANShybridsolverforsimulatingtwo-phaseflowbeneathahov-eringrotorC/The69thAnnualForumoftheAmericanHelicopterSociety.Alexandria:AHS,2013:2462-2484.7胡健平,徐国华,史勇杰,等.基于CFD-DEM耦合数值模拟的全尺寸直升机沙盲形成机理J.航空学报,2020,41(3):123363.HUJP,XUGH,SHIYJ,etal.Formationmechanismofbrown-outinfull-scalehelicopterbasedonCFD-DEMcouplingsnumeric-al simulationJ.Acta Aeronautica et Astronautica Sinica,2020,41(3):123363（inChinese）.8KELLERJD,WHITEHOUSEGR,WACHSPRESSDA,etal.Aphysics-basedmodelofrotorcraftbrownoutforflightsimulationap-plicationsC/TheAmericanHelicopterSociety62ndAnnualFor-um.Alexandria:AHS,2006:1097-1107.9WACHSPRESSDA,WHITEHOURSEGR,KELLERJD,etal.Ahighfidelitybrownoutmodelforreal-timeflightsimulationsandtrainersC/TheAmericanHelicopterSociety65thAnnualForum.Alexandria:AHS,2009:1281-1304.101062261014051015t=2 st=6 st=17 st=23 sz/mVx/(ms1)图12EH-60L 旋翼前端流场径向速度（r/R=1.0）Fig.12RadialvelocityatforwardpartofEH-60L(r/R=1.0)第6期谭剑锋，等：基于 DEM 的直升机沙盲加速计算方法1359ANDREAAD,SCORCELLETTIF.Enhancednumericalsimula-tionsofhelicopterlandingmaneuversinbrownoutconditionsC/TheAmericanHelicopterSociety66thAnnualForum.Alexandria:AHS,2010:1084-1101.11PHILLIPSC,KIMHW,BROWNRE.Theflowphysicsofheli-copter brownoutC/Present at the American Helicopter Society66thAnnualForum.Alexandria:AHS,2010:1273-1291.12SYALM,LEISHMANJG.Predictionsofbrownoutdustcloudscompared to photogrammetry measurementsJ.Journal of theAmericanHelicopterSociety,2013,58(2):1-18.13GOVINDARAJANBM,LEISHMANJG.Predictionsofrotorandrotor/airframeconfigurationaleffectsonbrownoutdustcloudsJ.JournalofAircraft,2016,53(2):545-560.14GOVINDARAJANB,LEISHMANJG,GUMEROVNA.Evalu-ationofparticleclusteringalgorithmsinthepredictionofbrownoutdustcloudsC/Presentedatthe67thAnnualForumoftheAmeric-anHelicopterSociety.Alexandria:AHS,2011:325-347.15HU Q,GUMEROV N A,SYAL M,et al.Toward improvedaeromechanicssimulationusingrecentadvancementsinscientificcomputingC/Presentedatthe67thAnnualForumoftheAmeric-anHelicopterSociety.Alexandria:AHS,2011:2105-2119.16PORCR,MIGLIOE,PAROLININ,etal.HPCsimulationsofbrownout:Anoninte