基于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,