李震东 许 磊 包士毅
(浙江工业大学化工机械设计研究所 杭州 310023)
(过程装备及其再制造教育部工程研究中心 杭州 310023)
(浙江工业大学机械工程学院 杭州 310032)
管道是天然气集输管中重要的输运设备,其中的流体会携带固体颗粒撞击管道壁面,进而造成壁面质量损失,产生冲蚀磨损问题[1]。管道冲蚀磨损过程较为复杂,与管内介质、管道材料、管道所处环境等因素息息相关[2]。有研究指出,管道中弯头处的冲蚀速率是直管段的50 倍左右[3]。因此,较为准确地预测不同条件下管道弯头处的磨损速率,对于预防输气管泄漏、减少经济损失等问题有较大意义。
随着计算流体力学的发展,由Tsuji 等人[4]提出了计算流体动力学-离散单元法(computational fluid method and discrete element method,CFD-DEM)。该方法被应用于二维流化床研究,且在随后的多相流研究中得到广泛应用[5]。为了模拟更为真实的工况,一些颗粒形状模型被开发出来,例如:超椭球模型、组合颗粒模型以及多面体模型[6]被用来对非球形颗粒进行建模。对于磨损理论的研究早在二十世纪五六十年代就开始了。1960 年Finnie[7]提出了脆性材料冲蚀模型,并以一弹性球体垂直撞击脆性材料为例分析脆性材料的冲蚀机理,随后多种冲蚀模型被学者们提出,这些模型对后来学者们运用数值模拟方法来研究磨损问题发挥了重要作用。基于上述理论,一些学者对气固两相流中管内的非球形颗粒流动进行研究。Hilton 等人[8]采用超椭球模型研究了在水平气力输送管道中,颗粒形状对稀释流到段塞流这一转换过程中的影响。Kruggel-Emden 等人[9]除了基本的球形颗粒外还建立了4 种多面体模型,讨论了这5 种颗粒在气力带动过程中流动特性的差异。周甲伟[10]研究了颗粒形状和旋流强度对弯管磨损的影响,其颗粒形状采用组合颗粒法进行建模,此外,他们还研究了煤块颗粒破碎过程。Zeng 等人[11]则分析了弯管中V 型磨损带形成的原因,并且比较了不同多面体形式颗粒对弯管磨损的影响。
本研究采用CFD-DEM 方法耦合切向撞击能量模型[12](shear impact energy model,SIEM)来研究不同形状颗粒对90 °弯管冲蚀磨损特性的影响。首先将标准算例与实验数据进行比较验证上述方法的可靠性,其次,对比不同形状颗粒对弯管的磨损情况。此外,将气速这一重要因素考虑在内,对比不同气速条件下颗粒形状对弯管磨损的影响程度。
由于本论文算例均为稀气-固两相流,为了节省计算时间,故均采用单向耦合方法,流体控制方程如下所示:
连续相方程为
动量守恒方程为
式中,ρf为流体密度(假设流体不可压缩,所以为常数),u为流体流速,p为流体压力,μeff为流体有效粘度,g为重力加速度。
离散相运动可分为平移和旋转,可以用牛顿第二定律来计算,具体如下:
式中,m为颗粒质量,dv/dt为加速度,Fc为接触力,Fd为流体对固体的曳力,FlS为Saffman 升力,FlM为Magnus 升力,Fb为浮力,I为惯性矩,dω/dt为角加速度,Tc为接触扭矩,Tf为流体作用下的扭矩。采用Cundall 等[13]提出的线弹性-阻尼模型计算Fc:
式中,Fc,n为接触力的法向分量,Fc,t为接触力的切向分量,kn为法向弹性系数,kt为切向弹性系数,δn为法向位移,δt为切向位移,ηn为法向阻尼系数,ηt为切向阻尼系数,vn为法向速度,vt为切向速度。阻尼系数可根据Ting 等人[14]提出的方程得到,当接触力2 个分量满足式(8),则接触力切向分量按照式(9)求解。
式中,fs为滑动摩擦系数。
由于非球形颗粒在运动中存在各向异性,引入局部坐标系来获得非球形的转动惯量I,但接触扭矩Tc仍按全局坐标系计算求解,具体如式(10)所示。
式中,L为从颗粒中心到接触点的矢量。
除了常规的球形模型以外,超椭球模型、组合颗粒模型以及多面体模型[6]也被广泛应用于DEM 建模方法中。本文采用由Barr[15]提出的超椭球方程来对非球形颗粒进行建模,具体公式如下:
式中,s1和s2为颗粒形状系数,a、b、c为轴线方向的半轴长。
在单向耦合中,只考虑流体对固体所产生的力,其中包括曳力Fd、Saffman 升力FlS、Magnus 升力FlM和扭矩Tf。
曳力Fd采用Di Felice[16]所提出的曳力模型,具体如下:
式中,Fd0为颗粒无其他粒子作用时流体对其产生的阻力,αf为流体体积分数,CD为流体曳力系数,dp为颗粒直径,vf为气体速度,vp为颗粒速度,Rep,α为流体雷诺数,μf为流体粘度。
颗粒在流体作用下的线速度和角速度都较大,不可忽略,故Saffman 升力[17]必须考虑,其表达式如下:
式中,ωf为流体角速度,ClS为升力系数,可根据Mei[18]的研究成果得到,具体如式(18)所示,Rep为颗粒雷诺数。
Magnus 升力是由颗粒旋转而引起的,具体如下:
式中,Rep为Magnus 升力作用下的雷诺数;Rer为颗粒旋转的雷诺数;ωp为颗粒角速度;ClM为升力系数,可根据Oesterlé 等人[19]的研究成果得到,具体如式(23)所示。
颗粒浮力Fb如下:
式中,Vp为单个颗粒体积。
颗粒的转动可根据Rubinow 等人[20]的研究成果得到,具体如下:
式中,CR为旋转扭矩系数,可根据Rubinow 等人[20]以及Dennis 等人[21]的研究内容中获得,具体如式(27)所示。
采用SIEM 模型[12]来计算弯管磨损情况,磨损模型如下所示:
式中,W为材料被碰撞后所损失的体积,EShear为颗粒碰撞壁面后颗粒损失的切向动能,p为壁面塑性流动压力,t0为碰撞开始时间,t为接触持续时间,Fc,t为颗粒与壁面接触时对壁面的切向力,vt为切向速度。
本文所采用的几何模型以及实验数据均来自Chen 等人[22]的研究内容及成果,实验装置示意图如图1 所示。实验管材选用铝,质量流量为40 lb/d(2.08 ×10-4kg/s)的沙子从颗粒入口注入后在空气压缩机的带动下,冲击测试单元处的弯管,并从出口流出,气体速度为150 ft/s(45.72 m/s),弯管特定位置的磨损率通过轮廓测量仪来获取。
图1 Chen 等人[22]所采用的实验装置示意图
计算中取实验示意图中颗粒经过弯管竖直向上冲击弯管的直管以及后续弯管部分作为模拟模型,具体尺寸如图2 所示。
图2 管道几何参数示意图
连续相采用商用Fluent 软件进行模拟,流体运动基于压力求解器,用SIMPLEC 来求解气体控制方程,对流相采用QUICK 方法。离散相固体颗粒采用商用DEM 软件DEMSLab 来进行模拟,具体参数如表1 所示,弯管网格划分如图3 所示。
表1 模拟参数表
图3 模拟采用的网格
由于本文只针对弯管磨损进行分析,因此只对弯管处网格无关性进行分析。将弯管部分网格进行调整,分别为21 008、23 028、25 048、27 068、29 088个网格单元,在相同边界条件下得到弯管外拱处磨损率,并与实验数据进行比较,结果如图4 所示。从模拟结果上看,采用25 048 个网格的计算精度已经达到要求,因此选择该规格网格继续进行本文研究。
图4 不同网格数量对应弯管外拱磨损情况与Chen 等人[22]的实验数据对比结果
模拟中颗粒分布如图5 所示(颗粒被放大10倍),从图中可以看出,2 条黑色箭头为颗粒在冲击弯管、经过壁面反弹后颗粒流的主要流向,这时颗粒在气力的带动下,二次冲击弯管、进而造成V 型磨损带,这与Zeng 等人[11]的研究结论一致。
图5 验证算例的管道中颗粒分布
为了研究颗粒形状对弯管冲蚀磨损的影响,本文基于Chen 等人[22]的实验条件,在45.72 m/s 的气速条件下(管道流场如图6 所示),采用超椭球模型建立6种颗粒模型,具体形状及尺寸如图7(a)~ (f)所示。其中a、b、c半轴长单位为m,后文图中的(a)~(f)对应6 种不同形状的颗粒模型,并以上述6 种颗粒模型为研究对象探究不同颗粒形状对弯管磨损的影响程度,单个不同形状颗粒的质量均相同。
图6 验证算例的管道中速度场与压力场云图分布
图7 颗粒形状示意图及其具体尺寸
图8 上图为不同形状颗粒对弯管拱背的磨损率,下图则表示不同形状颗粒在气力的带动下多数颗粒在与壁面碰撞前其速度能达到的最大值以及颗粒对弯管内壁造成多次碰撞后的最低速度。总体来看,不同形状颗粒造成的拱背处磨损率变化趋势并无明显差异,磨损率会随着角度的增加,在40 °~50 °之间达峰值,随后会迅速下降。
图8 不同形状颗粒造成的弯管拱背处磨损率以及对应情况下颗粒的最大速度以及最小速度
对比图8(a)~(d)中4 种颗粒模型下弯管拱背的磨损率,可以看出,随着形状系数的增加(颗粒球形度为球形颗粒表面积与与颗粒相同体积的颗粒表面积之比,正方体约为0.806;颗粒越接近方形,球形度越低),弯管的最大磨损率会先小幅下降,再逐渐上升;图8 还显示出弯管的最大磨损位置会朝大角度方向偏移,最终会稳定在47 °左右。最大磨损率的变化趋势与颗粒最大速度的变化趋势相同,即随着形状系数的增大颗粒的最大速度先有所下降,然后逐渐增大,这也直接说明弯管最大磨损率与颗粒最大速度有关,而颗粒的最大速度与颗粒受到的曳力有关,曳力越大,冲击速度也越大。此外,图8中不同形状颗粒的速度极差也能反映颗粒经过弯管对其产生的切向撞击能的相对大小。总体看,随着形状系数的增加,颗粒的速度极差越大,管壁受到的切向撞击能也越大,这也正好符合图9 中(a)~(d)总磨损率逐渐上升这一规律。尽管(b)对弯管拱背最大磨损率略低于(a),但(b)的总磨损率却较大,这也恰好说明随着形状系数的增大,颗粒流对弯管的磨损范围也会相应增大。
图9 不同形状颗粒造成的弯管总磨损率
通过颗粒在管壁附近的姿态可判断其对管壁的冲蚀磨损方式。欧拉角便是描述颗粒姿态的方法,XYZ坐标系先后绕着Z轴、X′轴以及Z″轴旋转Ψ度、θ度以及φ度得到局部坐标系X‴Y‴Z″,如图10所示,Ψ、θ、φ分别为欧拉角3 个分量进动角、章动角和自转角,这样颗粒在某一位置的姿态便可通过欧拉角表示。通过统计多组数据中颗粒在弯管拱背附近的平均姿态以及欧拉角3 个分量的标准差来表示大多数颗粒在靠近管壁时的状态。将90 °弯管分成18 份,分别对每段位置的颗粒欧拉角进行统计,由于球形颗粒的各向同性,故不需考虑,统计结果如图11 所示。通过颗粒3 个分量的平均值,可以画出3 种不同形状系数颗粒在弯管拱背附近XY面的投影情况,如图12 所示。不同形状颗粒的姿态在5 °~20 °、65 °~85 °之间差异明显,当颗粒不为球形时,弯管的磨损是颗粒滑动摩擦和滚动摩擦共同作用的结果,随着颗粒逐渐趋于方形,颗粒发生滑动摩擦的造成的磨损量增多,弯管的磨损量也会升高。总之颗粒形状的改变直接影响着颗粒作用于管道壁面的撞击能大小、作用范围以及通过影响颗粒所受曳力间接影响弯管磨损。
图10 欧拉角示意图[12]
图11 不同形状系数颗粒欧拉角3 个分量的均值和标准差
图12 不同形状系数颗粒参数获取位置以及颗粒在XY 面的投影
进一步分析长椭球(e)和扁椭球(f)对弯管磨损的结果,由于长椭球颗粒与扁椭球颗粒绕Z轴旋转任意角度都不会改变姿态,因此仅需考虑颗粒欧拉角前2 个分量的平均值和标准差,结果如图13 所示。两种颗粒在弯管拱背附近XY面的投影情况,如图14 所示。从图13 可以看出,拱背不同位置的颗粒,尽管欧拉角2 个分量的标准差都较大,但由于颗粒关于Z轴与X轴中心对称,图14 中颗粒的投影仍可以代表大多数颗粒在XY面的投影。从图13可以看出,拱背附近颗粒,在气力带动下,长椭球颗粒更易发生滑动摩擦,扁椭球更易发生滚动摩擦,颗粒发生滑动摩擦损失的切向动能多于滚动摩擦,造成弯管的磨损也更严重。从图8中(e)和(f)的速度极差也可说明长椭球颗粒对管壁的切向撞击能要高于扁椭球,因此在38 °~55 °位置处,长椭球造成的磨损率高于扁椭球,而在其他位置略低。这可能与颗粒-管壁的接触次数和接触时间有关。图9(d)的总磨损率明显高于图9(e),仍可以由2种颗粒的速度极差来判断,即由于单个(d)颗粒损失的动能高于(e),经过多次累加,就会造成形状接近方形的颗粒对弯管的磨损多于长椭球颗粒。
图13 长椭球与扁椭球颗粒进动角与章动角的均值和标准差
图14 颗粒参数获取位置以及颗粒在XY 面的投影
图15 为不同形状颗粒对整个弯管磨损率影响的情况。当颗粒流为球形颗粒时,磨损严重区域呈椭圆形,管壁两侧的V 形磨损带还不明显;在形状系数增大过程中,椭圆形变大,随后磨损严重区域沿着管壁两侧扩展,在形状系数为8 时扩展区域最为明显。图15 结果也表明长椭球对弯管造成的磨损低于形状系数为8 的颗粒,但高于扁椭球。模拟云图也可以说明,随着形状系数的增加,弯管磨损严重区域逐步扩大,V 形磨损带也越发明显,长椭球对弯管产生的冲蚀磨损高于扁椭球,这也印证了之前的结论。图16 为不同颗粒在管内的位置分布,不同颗粒流竖直向上冲击弯管,随后沿着弯管外拱以及管壁两侧流出,造成了不同的冲蚀形貌,同样可以观察到上述现象。
图15 不同形状颗粒对弯管冲蚀磨损云图
以图16(a)、(d)、(e)和(f)4 种颗粒形状为研究对象,以v=45.72 m/s 和v=30.48 m/s 这2 种气速为入口边界条件来对比不同气速下颗粒形状对弯管的磨损程度。
图17 为不同颗粒在2 种气速条件下对弯管拱背的磨损率,左图为v=45.72 m/s 气速条件,右为v=30.48 m/s 气速条件。从图中可以看出,无论改变气速还是颗粒形状,弯管拱背的磨损趋势都并无太大差别,且气速变动并未引起最大磨损率位置发生改变。
图18 为不同形状颗粒对弯管总磨损率以及对应弯管的磨损云图。对比同种颗粒形状在不同速度下的弯管总磨损率的差值可以看出,气速的增大使得更接近方形颗粒的算例(d)中弯管磨损率的升高最为明显。综合来看,气速的增大会使弯管各部分磨损量均匀上升,弯管磨损分布并无较大改变,磨损分布的不同主要是颗粒形状不同所导致。
(1)相同质量流量以及气速条件下,同质量的非球形颗粒较球形颗粒造成的弯管内壁磨损率更大,且颗粒形状越接近方形,弯管磨损越严重,磨损范围越广。
(2)相比于扁椭球,长椭球颗粒造成的弯管磨损更大,这主要由于长椭球颗粒在管壁面更易发生滑动摩擦,扁椭球易发生滚动摩擦。
(3)形状近似为方形的颗粒进出弯管损失的动能要多于长椭球,这就导致长椭球颗粒对弯管冲击造成的弯管质量损失较少。
(4)气速的提升仅会均匀增大不同形状颗粒对弯管的磨损率,而磨损分布以及最大磨损率位置并无较大差异。