工业过程中的流体输运常常夹杂着不同种类、大小、形状和浓度的固体颗粒，颗粒随流体流动，冲击管道内壁、阀门和其他部件，可能导致设备故障甚至损坏[1-2]．管路中最易磨损的是弯头，颗粒在弯头内受二次流动的影响发生了流动方向的偏转，并造成局部严重的磨损问题；但是，对弯头内磨损机理的研究是一个复杂的问题，它涉及到诸多微观和复杂的两相流动及壁面侵蚀磨损机制的研究[3-4]．离散元(DEM)作为研究颗粒运动的计算机模拟方法在石油化工、农业装备等领域中广泛应用，能够准确计算颗粒与颗粒、颗粒与壁面的碰撞，并可以与流体运动进行耦合计算[5-10]．随着CFD-DEM方法的快速发展，文献[5-6]提出该数值方法来模拟氧化铝颗粒通过丙烯酸弯管的磨损，模拟与试验结果符合较好．文献[11]在考虑颗粒-颗粒和颗粒-壁面相互作用的基础上，提出了基于CFD-DEM的气固两相流磨损预测模型．CFD-DEM模拟方法在求解颗粒两相流流动及管壁磨损中具有较高准确性，但目前针对气固两相流弯管磨损的研究较多，针对液固两相流数值模型中各个机制系统的研究较少．本研究采用CFD-DEM耦合数值模拟方法，对承载液固两相流动的一个90º弯管内固体颗粒引起的弯头磨损进行了数值模拟研究．通过自定义功能将discrete random walk(离散随机游走)模型和颗粒对湍流调制作用定义到CFD-DEM模型中，建立了流体与颗粒全耦合的数值模型．采用考虑颗粒体积分数的欧拉-拉格朗日方法，同时考虑颗粒-颗粒、颗粒-壁面相互作用下单向与双向耦合方法、升力与其他次要力、重力施加方向对弯管磨损模拟的影响．最后，研究颗粒与湍流相互作用对粗颗粒造成弯头磨损预测的影响．1 数值方法和计算模型1.1　连续相计算方法连续相的运动主要由质量守恒方程和纳维-斯托克斯方程来计算，其中颗粒对流体的作用力Fs-w表示为Fs-w=-1V∑i=1nfs-w，式中：fs-w为作用在单个颗粒上流体的合力；V为流体单元的体积．研究对象为粗颗粒，且颗粒体积分数较大，须要计算流体体积分数．由于弯管内产生二次流，因此选用RNG k-ε湍流模型．由颗粒存在产生的湍动能和湍流耗散率源项，表示为Sk=1V∑k(vf+vf'-vp)Fk'; (1)Sε=ρfCε3ν2V∑ksk|vf+vf'-vp|2/dp, (2)式中：vf和vp分别为流体和颗粒的平均速度；vf'为流体的脉动速度；Fk'为流体作用在颗粒上的表面力，包括曳力、升力和虚拟质量力等；dp为颗粒粒径；ρf为流体密度；sk=CDRep/24为阻力因子，表示阻力系数与斯托克斯阻力的比值，其中颗粒雷诺数Rep=ρfdp|vf-vp|/μf，CD为阻力系数；Cε3=Cε2'CεpRep1.416，其中Cε2'≈Cε2=1.92，Cεp=0.058．1.2　离散相控制方程在拉格朗日坐标下，采用牛顿第二定律计算颗粒的平移和旋转运动，其控制方程为mpdvp/dt=Fdra+Fg+Fpr+Fvm+FSaf+FMag+Fc+Fbu;Ipdωp/dt=∑Tc+Tf，式中：Fdra为颗粒受到的曳力，对于多颗粒工况以及考虑颗粒体积分数，采用文献[12]提出的曳力模型；Fg为颗粒自身的重力；Fbu为颗粒受到的浮力；Fpr为压力梯度力；Fvm为虚拟质量力；FSaf为流体剪切力产生的Saffman升力；FMag为由于颗粒自旋产生的Magnus升力，对应计算公式通过自定义程序定义到模型中并计算；Fc为颗粒-颗粒、颗粒-壁面间接触力，采用Hertz-Mindlin(无滑移)接触模型；Ip为颗粒惯性力矩；Tc为颗粒接触力矩；Tf为流体作用产生的力矩；ωp，t，mp分别为颗粒角速度、时间和颗粒质量．1.3　离散随机游走模型湍流对颗粒运动的影响可以用离散随机行走(DRW)模型来模拟[13]．湍流速度波动服从高斯分布，通过以下关系得到，即u'=ξ1(u¯'2)1/2;v'=ξ2(v¯'2)1/2;w'=ξ3(w¯'2)1/2,式中：ξ1，ξ2，ξ3为服从高斯分布的随机数；u'，v'，w'为流体在x，y，z三个坐标下的脉动速度．由于RANS模型的简单性，假设湍流为各向同性，流体脉动速度的均方根为(u¯'2)1/2=(v'2)1/2=(w¯'2)1/2=(2k/3)1/2，式中k为湍流动能，用RNG k-ε模型得到．离散随机游走模型考虑了脉动速度对颗粒运动影响，与式(1)和(2)共同组成颗粒-流体湍流相互作用模型．1.4　磨损模型颗粒磨损是一个非常复杂的过程，还没有具体普遍适用性的模型，但基于不同条件下的应用，已经提出了许多预测磨损的经验公式．使用的磨损预测模型应考虑尽可能多的因素，以满足某些特定的使用要求．本研究考虑液固两相流动、弯头几何特性和材料特性，采用美国塔尔萨大学磨损腐蚀研究中心开发的磨损模型[4]，计算公式为ER=CBH-0.59Fsvp'dFθ；Fθ=∑i=15Aiθi，式中：ER为磨损率，定义为每千克撞击在壁面上的颗粒从壁面上除去的材料质量；BH为壁面材料的布氏硬度；Fs为颗粒的尖锐程度，对于圆角、半圆角和锋利的颗粒，取值为0.2，0.53和1；vp'为颗粒碰撞速度；θ为以弧度制表示的撞击角度；d=2.41为速度指数；C=2.17×10-7为经验常数；Ai (i=1，2，…，5)值分别为5.398 3，-10.106 8，10.932 7，-6.328 3，1.423 4．1.5　物理模型选用内径D为30 mm的标准弯管作为研究对象，为确保流体和颗粒可以充分发展，定义进口直管和出口管长度为20D，弯头弯径比R/D=1.5．颗粒工厂位置设置在距离进口2D的地方，弯头弯曲角度从0°~90°均分，几何模型如图1所示．该几何模型的网格划分采用O型网格法，考虑到须计算颗粒体积分数，边界层网格厚度太小会影响计算的准确性，边界层网格y+为100，保证壁面函数应用的合理性．弯管进口定义为5 m/s的速度入口，出口设置为压力出口，其参考静压为一个标准大气压，整个计算域壁面假设为无滑移壁面．在Fluent与EDEM耦合计算中采用相间耦合的SIMPL算法来耦合流场中的速度场和压力场，颗粒间的接触计算采用Hertz-Mindlin(无滑移)模型，颗粒与壁面间的接触采用添加了磨损模型的Hertz-Mindlin(无滑移)模型．颗粒密度为2 650 kg/m3；颗粒泊松比为0.17；颗粒入射速度为4 m/s；颗粒-颗粒恢复系数为0.95；颗粒-颗粒静摩擦系数为0.005；颗粒-颗粒滚动摩擦系数为0.5；壁面密度为7 800 kg/m3；泊松比为0.3；弹性模量为200 GPa；颗粒-壁面恢复系数为0.737；颗粒-壁面静摩擦系数为0.2；颗粒-壁面滚动摩擦系数0.3[14]．10.13245/j.hust.211008.F001图1几何模型及计算域网格本研究首先在EDEM软件中建立了颗粒模型，在Fluent中求解了稳态条件下弯道内的流场，然后将EDEM与Fluent进行双向耦合，颗粒的生成速率由模拟的颗粒体积分数确定．本研究设置EDEM计算时间步长为5×10-6 s，Fluent计算时间步长为5×10-5 s，计算总时长为3 s．本研究建立了网格数分别为7.441 8×104，1.383 16×105，2.570 08×105的三套模型，对比了三套网格在弯管进口处流体速度剖面，发现粗网格与参考网格误差为0.84%，细网格与参考网格误差为0.17%，最终采用参考网格进行模拟，网格如图1所示．1.6　计算模型验证文献[15]针对水平管水力输送进行了一系列试验，研究了不同粒径和体积分数的砂在水平管道中的输送．利用本研究所建立的全耦合CFD-DEM数值模型方法，对水平管道输送平均粒径为2.4 mm、体积分数为15%的颗粒进行了数值模拟．模拟水平管道的直径为D=53.2 mm，管长为10D．在Fluent和EDEM中都定义了周期边界条件．在Fluent中设置流体质量流量，采用标准k-ε湍流模型．Fluent和EDEM所用的时间步长分别为5×10-5 s和5×10-6 s．采用单向耦合方法模拟注入颗粒随流体的运动，时间为2.5 s，然后采用双向耦合法(CFD-DEM耦合)计算注入颗粒的瞬态状态5 s．图2给出了颗粒体积分数与y/D的关系图(y为沿y方向的直径坐标值)，并与试验数据[15]和之前的数值模拟数据[14]进行了比较，图中φ1为颗粒体积分数．在不同y方向获得颗粒体积分数的方案如下：沿管道创建一系列横截面，水平距离为100 mm，同样，在每个平面创建10条等距水平线；各个高度处的颗粒体积分数均由该参数沿直线对所有监测截面进行平均得到；最后对获得的数据进行时间平均，即在最后3 s内，以0.1 s的时间间隔对计算的颗粒体积分数进行时间平均．试验结果对应的条件为：φ1=0.15；vp=1.8 m/s．前人模拟结果对应的条件为：φ1=0.152；vp=1.8m/s．本研究模拟采用的颗粒体积分数和入口速度与试验条件保持一致．结果表明本文计算结果与试验结果符合较好．10.13245/j.hust.211008.F002图2水平管颗粒体积分数试验与模拟对比2 结果与讨论2.1　单向与双向耦合对弯头磨损预测的影响采用单向耦合方法模拟低颗粒体积分数条件下的液固两相流动需要较少的计算资源，并且加快了计算时间．然而，文献[16]认为：即使当固体体积分数低至1×10-3时，也应该考虑颗粒对流体的影响，并且弯管中流体与颗粒之间复杂的相互作用及离心效应使得颗粒局部聚集，对流场产生显著影响．本研究使用单向和双向耦合方法模拟了0.5%和2.5%的颗粒体积分数．弯头最外侧壁面处的磨损率如图3所示，图中α为角度．0.5%和2.5%的磨损率曲线非常相似，而磨损率本身增加了10倍左右，单向耦合和双向耦合预测的磨损率存在显著差异，当弯头弯曲角度大于40°时，单向耦合方法预测的磨损率明显高于双向耦合方法预测的磨损率．当采用单向耦合方法时，一部分颗粒直接冲击弯头45°~70°的区域．在双向耦合方法中，颗粒运动与流场有较好的跟随性．10.13245/j.hust.211008.F003图3弯头外侧壁面磨损率1—单向耦合，φ1=2.5%；2—双向耦合，φ1=2.5%；3—单向耦合，φ1=0.5%；4—双向耦合，φ1=0.5%．图4为弯管磨损率云图，由图可见：单向耦合方法下可以观察到两个磨损最严重的区域，而双向耦合在85°出现一个磨损峰．本研究结果与文献[17]进行了定性比较，磨损分布形貌试验符合较好．结果表明：单向耦合弯管磨损预测结果与试验具有较大差异，耦合方法会影响颗粒运动预测的精度和磨损结果，在弯头磨损模拟中，即使颗粒体积分数低至0.5%，也应采用双向耦合方法．10.13245/j.hust.211008.F004图4弯管磨损率云图2.2　颗粒受力对弯头磨损预测的影响弯头处的磨损主要由颗粒的分布决定，颗粒的受力严格控制着颗粒的运动轨迹．在本研究中，这些力分为两类，即升力(包括Saffmann力和Magnus力)和其他次要力(包括虚拟质量力和压力梯度力)．图5给出了不同组合力作用下沿弯头外侧壁面的磨损率．结果表明升力可以对预测的磨损率产生显著影响．在考虑升力的情况下，当弯曲角度大于30°时，磨损率相对较高，而次要力对弯头磨损率的影响很小．在考虑升力的情况下，颗粒-壁面碰撞频率增加，而次要力对碰撞频率的影响不大．10.13245/j.hust.211008.F0051—无升力+无次要力；2—升力+无次要力；3—无升力+次要力；4—升力+次要力．图5　弯头外侧壁面磨损率通过阐明弯头处的颗粒分布，可以进一步研究升力在弯头处磨损所起的作用．图6给出了两种情况下沿弯道不同截面的颗粒体积分数．在0°和30°断面内，颗粒扩散良好，升力对弯道内颗粒分布影响不大；在60°和90°横截面上的颗粒都聚集在弯头的外侧壁面区，但两种情况下粒颗粒集中区域的周向分布不同．在没有升力的情况下，颗粒在周向扩散得更广，集中区域的厚度更薄．在有升力的情况下，颗粒集中区域的体积分数比没有升力的情况下要高．结果表明：升力在周向能形成较大的颗粒集中区，因此沿弯头外侧壁面的磨损率较大．在考虑升力的情况下，这种差异会导致沿弯头背面的磨损更为严重．10.13245/j.hust.211008.F006图6颗粒体积分数分布2.3　重力方向对弯头磨损预测的影响在实际应用中，管道布置方式往往不尽相同，其主要有五种方式，如图7所示，其中H和V分别表示管路水平和竖直布置．管道布置方式变化会导致颗粒分布受重力方向影响，在不同重力施加方向对弯管造成的磨损会显著不同，但磨损变化规律基本相似．图7所示为五种弯头布置方式下，弯头最外侧弧线上磨损率随弯头角度的分布．弯头随着弯曲角度的增大，磨损率都逐渐增大，并且都在85°左右出现最大磨损峰值．V-H，重力朝上方式的弯管磨损率最小，为7.56 nm/s；V-H，重力朝下方式的磨损率最大，达到9.69 nm/s；其余三种弯管最大磨损量相近，基本都在8.2 nm/s左右．最小磨损率与最大磨损率相差21.9%．结果表明：重力在弯管磨损预测中起重要作用，若重力与颗粒的运动方向相同，则这些颗粒会被重力加速并且以更高的速度与弯头碰撞，造成较大的磨损；相反，颗粒会减速并以较低的速度与弯头发生碰撞，造成较小的磨损．本研究结果表明：弯头安装时应尽量避免采用V-H，重力朝下方式布置，从而可以减少弯头磨损20%左右．10.13245/j.hust.211008.F007图7弯头外侧壁面磨损率1—H-H；2—H-V，重力朝上；3—H-V，重力朝下；4—V-H，重力朝上；5—V-H，重力朝下．2.4　颗粒与湍流作用对弯头磨损预测的影响模拟中通过湍流调制模型来预测颗粒-湍流相互作用对弯头外侧壁面磨损的影响．由于颗粒与湍流的相互作用受颗粒体积分数的影响较大，因此对体积分数为2.5%和0.5%的颗粒进行了模拟比较．在无湍流调制模型、只有湍流对颗粒的影响模型(DRW模型)和DRW与湍流调制模型结合的情况下，对弯头磨损进行了预测．除颗粒体积分数外，其他参数与参考情况相同．对于1 mm颗粒，图8给出了三种情况预测的沿弯头外侧壁面弧线上的磨损速率．结果表明，不考虑湍流对颗粒的影响(无DRW模型)情况预测的弯头外侧壁面的磨损率较低，而其他两种情况的预测结果非常相似．对于体积分数为2.5%的1 mm颗粒，DRW模型在预测弯头内的磨损率方面效果显著，而在较高体积分数下，湍流调制模型在弯管磨损预测方面影响则比较微弱，对弯头磨损预测计算结果几乎没有影响．因此，在颗粒高体积分数工况下对弯管磨损的预测中可以忽略湍流调制模型的影响．10.13245/j.hust.211008.F008图8弯头外侧壁面处磨损率1—无湍流调制模型+DRW模型；2—湍流调制模型+DRW模型；3—无湍流调制模型+无DRW模型．通过沿弯头的一系列截面上的湍流动能(k)可以看出颗粒对湍流的影响，如图9所示．很明显，k由于颗粒的存在而受到抑制，特别是在颗粒集中的区域．降低k值可以减小流体的速度波动及由此产生的与颗粒速度的相互作用．然而，湍流动能的减小对最终磨损率的影响不大．这可能是由于颗粒尺寸相对较大，而颗粒的局部聚集在弯头处形成了高度集中的颗粒区域．决定颗粒碰撞速度的主要机制是由于颗粒在弯头中方向发生偏转而受到惯性效应，使得颗粒在弯头局部区域发生聚集增加了颗粒-颗粒间的碰撞．10.13245/j.hust.211008.F009图9湍动能分布在弯管内模拟了直径为1 mm、体积分数为0.5%的颗粒，图10给出了考虑湍流调制模型和没有考虑湍流调制模型预测的沿弯头外侧壁面的磨损率．虽然颗粒直径仍为1 mm，但考虑湍流调制模型预测的磨损率要高于未考虑湍流调制模型预测的磨损率，弯头弯曲角度在30°~45°间有较大差异．在考虑湍流调制模型作用下，颗粒碰撞频率较高．当弯头外侧附近的颗粒体积分数比为0.5%时，颗粒体积分数要低得多．在弯曲角度30°~60°之间，弯头附近的颗粒体积分数小于10%，颗粒-颗粒碰撞并没有在颗粒碰撞运动中起主导作用，在这种情况下，湍流调制模型对颗粒运动的影响不可以被忽略．当颗粒体积分数沿弯头弯曲角度增加时，颗粒-颗粒碰撞在决定颗粒碰撞运动和由此产生的磨损率方面超过了其他机制．10.13245/j.hust.211008.F010图10弯管外侧壁面磨损率3 结论a. 在承载固-液两相流动的弯管磨损研究中，无论是大体积分数工况还是小体积分数工况，双向耦合方法是必须要采用的计算方法．b. 在弯管最外侧壁面处升力对弯管磨损影响显著，而在弯管外侧轮廓线30°~90°处，没有升力时产生的磨损量更大，这主要是由于在升力作用下改变了颗粒分布及运动情况，使得颗粒更向弯管外侧周向聚集，局部体积分数较高．c. 不同重力施加方向对弯管造成的磨损规律基本相似，磨损随弯管弯曲角度的增大而增大，并且最大磨损在弯管85°左右，其中整个弯管最大和最小磨损量都在V-H型弯管中．在弯头安装过程中应尽量避免采用V-H，重力朝下方式布置，从而可以减少弯头磨损20%左右．d. 在大直径较高体积分数工况下，在弯管最外侧壁面上湍流调制模型对磨损影响不大，而DRW模型具有较为显著的影响，主要是颗粒碰撞占主导机制．在颗粒小体积分数工况下，颗粒湍流调制是必须要考虑的影响因素．