非球形固体颗粒对弯管内壁的磨损机制

曹 骞,康 灿*,滕 爽,焦 侬,丁可金

(1.江苏大学 能源与动力工程学院,江苏 镇江 212013;2.中国船舶集团第七〇四研究所,上海 200031)

作为采矿系统的关键部件之一，弯管起到连接各直管段的关键作用.在矿浆输送过程中，固体颗粒并不完全跟随液流运动，尤其在曲率较大的弯管中，颗粒与壁面将发生碰撞和摩擦，最终导致壁面材料去除，即磨损[1-3].单个颗粒单次作用于壁面而产生的磨损可以忽略，但长时间内大量颗粒作用于同一壁面位置，将造成壁面穿透甚至管路失效，不仅引起经济损失，还会带来严重的安全隐患[4-6].研究颗粒与管道内壁之间的相互作用具有重要的意义.在实际输送过程中，矿石颗粒形状各异，颗粒的形状无法准确描述，也很难用统一的标准去评价颗粒形状的不规则程度.目前，在固液两相流动的研究中，针对颗粒形状开展的研究较少，对不同形状颗粒与管壁之间的摩擦与磨损的认识尚不充分[7].

针对管内固液两相流动的研究主要从试验与数值模拟两个方面开展.Alajbegovic等[8]采用激光多普勒测速仪(LDA)分别测量了竖直管道内颗粒和流体的速度分布，进而分析了颗粒对管内湍流的抑制作用；Kesana等[9]借助探头在水平放置的弯管内不同位置进行取样，研究颗粒直径和流体黏度对颗粒分布的影响.欧拉-欧拉和欧拉-拉格朗日方法是处理固液两相流动时最常用的两种数值模拟方法，前者将离散相处理为“拟流体”，在欧拉框架下对颗粒进行求解，该方法适用于颗粒浓度较高而直径较小的情况；后者采用N-S方程(即Navier-Stokes方程，用以描述黏性不可压缩流体的动量守恒)在欧拉框架下对连续相进行求解，同时运用牛顿第二定律在拉格朗日框架下对离散相进行求解[10].

早期关于颗粒与管壁之间摩擦及壁面磨损的研究主要依靠试验，根据试验结果建立经验及半经验的磨损预测模型，描述两相流动状态对壁面磨损的影响[11-12].Archard模型是最早建立、应用最为广泛的磨损模型之一，该模型认为磨损量与颗粒和壁面接触时的法向力和滑动距离成正比，与材料的硬度成反比[13].Finnie等[14-15]首次提出塑性材料的微切削理论，建立了考虑碰撞速度及碰撞角相关的磨损模型；Oka等[16-17]讨论了材料类型、冲击角、冲击速度以及颗粒尺寸对磨损的影响，建立了Oka模型；在Zhang的模型中，额外考虑了颗粒形状对壁面磨损的影响[18].随着数值模拟技术的快速发展，采用磨损模型配合多相流模型来预测流体机械过流部件的磨损已成为重要方法之一[19-20].

本文中采用DEM与CFD耦合的方法，对一弯管内的固液两相流动及固体颗粒对管壁的磨损进行数值模拟，借助商用离散元软件EDEM和计算流体动力学软件ANSYS Fluent，并充分考虑颗粒与颗粒、颗粒与流体以及颗粒与壁面之间的相互作用，模拟弯管内的流动特征、颗粒的运动规律及弯管壁面的定量磨损规律；通过与试验结果进行对比，验证模拟的有效性；进而改变颗粒形状，研究颗粒形状对颗粒运动及弯管内壁磨损的影响.

1 弯管磨损试验方法

弯管磨损试验在如图1所示的磨损试验台上进行.该试验台主要由料箱、闸阀、驱动泵及试验段构成.闸阀1和闸阀2分别用于控制清水的流量和颗粒的浓度，在试验段的进口和出口各安装1个压力表以监测管路压力.固液两相混合物在离心泵的驱动作用下进入弯管.试验时所采用颗粒的形状不规则，取一定数量的颗粒，通过测量颗粒的密度和质量，计算出颗粒的当量直径约为4.5 mm，颗粒相的体积分数约为5%.试验中采用内径为40 mm，弯径比(弯管中轴线曲率半径与管径之比)为1.5的弯管进行磨损试验.

Fig.1 Test bench for wear图1 磨损试验台

2 数值模型及模拟方案

2.1 流动控制方程

本研究中采用双向耦合方法来处理颗粒与液相间的相互作用力，所有计算域的流体相均在绝对坐标系中进行求解.流体流动的控制方程如下：

连续性方程：

动量守恒方程：

其中t为时间；ρf为流体密度；u为流体的速度矢量；∇为哈密尔顿算子；αf为流体相的体积分数；P为流体压力；μ为流体的动力黏度；g为重力加速度，取值为9.81 m/s2；F为颗粒对流体的作用力.

2.2 颗粒建模及运动方程

为考虑颗粒形状对壁面磨损的影响，此处引入球形度的概念，用以描述颗粒形状接近球体的程度.颗粒的球形度定义为与颗粒相同体积的球体的表面积与颗粒的表面积之比.颗粒的球形度越大，其形状越接近于球体.

颗粒球形度的表达式为

其中φ为颗粒的球形度；Vp为颗粒体积；Sp为颗粒的表面积.

在模拟中采用四面体(φ=0.67)、六面体(φ=0.81)、八面体(φ=0.84)、十二面体(φ=0.91)以及二十面体(φ=0.94)形状的颗粒.在EDEM中对非球形颗粒进行填充，即采用若干小颗粒填充获得目标形状颗粒，如图2所示.填充的颗粒数越多，填充后颗粒的形状越接近目标颗粒，但颗粒数的增多会加大计算机的运算负荷[21].兼顾计算的准确性和经济性，每个目标颗粒所填充的小颗粒数在250～300之间.

Fig.2 Models of non-spherical particles图2 非球形颗粒模型

颗粒的受力决定了颗粒的轨迹，进而影响弯管壁面的磨损特性.颗粒的运动通过牛顿第二定律来描述.颗粒主要受到重力、曳力FD、Saffman升力FS、Magnus升力FM、压力梯度力FP、虚拟质量力FV以及颗粒间的相互作用力FC.

式(4)中,mp表示颗粒的质量；dup/dt为颗粒的加速度.在式(5)中，I为颗粒的转动惯量；dω/dt为颗粒的角加速度，Tc为颗粒间接触产生的力矩，Tf为流体施加在颗粒上的力矩.

在颗粒受到的各力中，重力和曳力对颗粒运动的影响最大[22].阻力的表达式为

其中[23]：

上式中，CD为曳力系数，dp为颗粒直径，up为颗粒运动速度矢量，Rep为颗粒雷诺数.

固体颗粒的Stokes数(St)表征颗粒惯性力与曳力的相对大小，其表达式为

上式中，D为管道的当量直径.Stokes数越小，固体颗粒对流体的跟随性越好[24].

2.3 计算模型和网格划分

本文中所采用的弯管模型由进口段直管、弯管以及出口段直管构成，如图3所示.其中弯管部分由试验弯管按照1:1的比例构建.为保证弯管入口流动充分发展、颗粒分布均匀以及出口无回流，进口直管和出口的长度均设为管径(D)的10倍.图中θ为极角，定义弯管进口截面θ=0°，出口截面位于θ=90°位置.γ用于表征弯管壁面上点的位置，弯管外侧脊线γ=0°，内侧脊线γ=180° (-180°)，如图3所示.

Fig.3 Computational model图3 计算模型

采用六面体结构化网格对计算域进行离散，以保证计算的稳定性以及计算结果的准确性.近壁面网格的y+分布在30～300之间.为验证网格数对计算结果的影响，设计了网格数不同的5套计算域网格，在相同的边界及工况条件下进行计算，比较不同网格数方案获得的弯管进出口压差，结果列于表1中.随着网格数的增加，弯管进出口压差的相对变化量减小.最终选择了网格数约为53万的网格方案(方案4).

表1 网格无关性验证Table 1 Grid independence examination

2.4 计算方法及边界条件设置

固液两相之间的耦合计算在欧拉-拉格朗日框架下进行.利用ANSYS Fluent在欧拉坐标系中对液相进行求解，控制方程为基于雷诺平均的Navier-Stokes方程，采用Realizablek-ε湍流模型，通过设置入口湍动能(k)和湍动能耗散率(ε)使控制方程封闭，采用SIMPLEC算法耦合流场中的压力和速度.计算域的进口采用速度进口边界条件，出口设置为Outflow边界.近壁面区采用标准壁面函数处理，壁面粗糙度设为0.046 mm.利用EDEM在拉格朗日坐标系中对固相进行求解.颗粒与颗粒和颗粒与壁面之间的接触模型采用Hertz-Mindlin无滑移模型[20].颗粒与壁面的属性列于表2中.通过耦合接口，实现EDEM和ANSYS Fluent之间的颗粒位置和动量等信息的传递[25].

表2 材料属性及相互作用设置Table 2 Material properties and interaction parameters

2.5 磨损模型

目前大多数磨损模型基于经验公式建立，无法将所有磨损因素考虑在内.本文中选用Zhang开发的用于计算冲蚀的经验模型，该模型认为壁面磨损主要与颗粒形状、碰撞速度及碰撞角度有关[18,26].该模型的表达式为式中，ER表示侵蚀率，为单位质量的颗粒碰撞壁面引起的材料损失；BH为壁面材料的布氏硬度；Fs为颗粒的形状因子，对于尖锐颗粒，Fs=1.0，对于球形颗粒，Fs=0.2；Vp为颗粒碰撞壁面时的速度；α为用弧度表示的颗粒碰撞角度；经验系数n=2.41，C=2.17×10-7.

本文中通过编程开发软件Visual Studio完成模型编译，并通过API (Application Programming Interface)嵌入EDEM，从而获得颗粒与壁面的碰撞以及壁面磨损信息.

在EDEM中，通过单位面积的磨损深度(h)来表征壁面磨损情况：

其中m为与壁面碰撞的颗粒的质量，A为计算单元的面积，ρw为壁面材料的密度.

EDEM根据颗粒与壁面间的相对速度和力识别磨损区域的受力.切向累积能对应颗粒在壁面滑动以及小角度碰撞所产生的能量，法向累积能对应颗粒以大角度碰撞壁面产生的能量.切向累积接触能(Et)和法向累积接触能(En)的表达式为

式中，Ft为切向力；Vt为切向相对速度；δt为时间步长；Fn为法向力；Vn为法向相对速度.

3 DEM-CFD数值模拟方法的验证

为验证数值模拟方案的物理有效性，以Alajbegovic等在垂直管道内固液两相流动特性试验中获得的结果作为参照[8].Alajbegovic等利用激光多普勒测速仪测量了垂直管道内2 200 mm高度处颗粒和流体沿径向位置的轴向速度.此处首先构建了与文献[8]一致的垂直管道模型，设定的工况与文献[8]保持一致，验证模型如图4所示.借助本文构建的数值模拟模型对该问题进行了数值模拟，验证模型中各参数的设置列于表3中.

表3 验证模型参数Table 3 Parameters of validation model

Fig.4 Validation model with particles traveling in a vertical pipe图4 垂直管道验证模型

试验和模拟的结果对比如图5所示.其中，r为距中轴线的距离，R为管道半径.可以看出，在主流区内，流体速度明显高于颗粒速度，这是由于相间滑移所致，数值模拟结果与试验结果表现出的趋势一致.沿径向方向，颗粒和流体的速度均逐渐下降；受边界层的影响，靠近壁面的流体和颗粒的速度出现较大的变化梯度.试验结果中各个位置的颗粒速度略高于模拟结果中的对应值.对于液相速度，试验结果与模拟结果接近.综合考虑试验误差及数值模型中无法考虑的实际因素，本文中的数值模型有效.

Fig.5 Comparison between numerical and test results图5 模拟和试验结果对比

4 结果分析及讨论

4.1 弯管内流动及磨损分析

试验中采用的固体颗粒形状不一，平均当量直径为4.5 mm.为了使模拟工况接近试验工况，模拟中的固体颗粒由图3中五种不同形状系数的颗粒构成，每种颗粒所占的固相体积分数为1%.在输送的固液两相混合物中，固相体积分数为5%.

图6所示为通过模拟得到管道内颗粒的质量(ms)随启动时间(t)的变化，即在管道内的固液两相混合物开始流动的初始阶段所获得的结果.可以看出，自初始时刻开始，管道内的颗粒数目持续增加，在t=0.14 s时刻，管道内固体颗粒的总质量达到最大，而后管道内的颗粒质量达到动态平衡，流动状态趋于稳定.为了获得稳定的磨损模拟结果，对t=0.50 s后的结果进一步进行分析.

Fig.6 Variation of particle mass in the pipe in initial stage of transport图6 初始阶段管道内的颗粒质量

4.1.1 弯管内的流动特征

上述工况下的液相速度和压力分布分别如图7和图8所示.液体的流动与颗粒的运动密切相关，从而影响固体颗粒对壁面的磨损.在图7(a)中可以明显看出，流体在弯管内侧出现局部加速现象，而在流经弯管后出现低速区，断面流速分布不均匀.从图7(b)的压力分布来看，其基本上与图7(a)所示的速度分布相呼应.

Fig.7 Velocity and pressure distributions of liquid图7 液相的速度分布和压力分布

流体在经过90°弯管时，在离心力的作用下，出现强烈的二次流，这从图7中的流动速度分布可以看出，该现象也与文献[27]中获得的结果相似.弯管内部的流体向曲率较大的外壁流动，由于流体流动具有连续性，外壁附近的流体被迫向弯管内部流动，形成复杂的空间流动结构[28].

沿弯管进口(0°截面)位置向弯管出口(90°截面)方向创建角度间隔为30°的截面.在弯管进口位置，流体受到离心力的作用较小，不存在二次流动，此处指向流线曲率中心的顺压梯度较大，流体的径向速度指向内侧壁面，如图8所示.自30°截面位置开始，过流断面上出现旋向相反的对称涡，此时颗粒对流体的扰动较强，涡核向弯管内侧壁面靠拢，低速运动的液体被旋涡卷吸，并随着主流向下游发展[29].从旋涡形态的演变来看，涡核的位置不断向内侧壁面移动，同时涡的强度降低，而旋涡的形状趋于平滑和规则，对称性得到保持.同时，在流体从进口向出口流动的过程中，受到的离心力作用不断累积，使得内侧壁面周围的流体速度不断降低，同时外侧壁面周围流体速度不断升高，如图8所示.

Fig.8 Patterns of liquid flow on different cross sections图8 不同截面的液流形态

颗粒的轨迹及运动矢量分布如图9所示，颗粒与流体以相同的速度自进口均匀进入计算域内.由于边界层的影响，靠近管壁的颗粒在上升过程中较主流区颗粒的速度低[30].颗粒进入弯管后，靠近外侧壁面的颗粒速度逐渐下降；同时，位于主流区的颗粒在惯性作用下继续向下游运动，受流场的影响速度下降，在与弯管外侧壁面发生第一次碰撞后，速度显著降低.颗粒速度矢量箭头的长度表征颗粒所受曳力的大小，弯管内侧壁面附近的颗粒自高速液流中获得能量，与弯管外侧壁面产生碰撞.在γ=50°～100°截面区间内，颗粒速度最低.颗粒离开弯管后，速度逐渐上升.

Fig.9 Particle trajectory and velocity vector distribution图9 颗粒轨迹及运动矢量分布

部分颗粒与γ=±(45°～60°)区间壁面发生碰撞后，速度降低，与流体之间的速度差增大，所受到的曳力作用更强.颗粒的Stokes数表征颗粒惯性力与曳力的相对大小，颗粒的Stokes数越小，则颗粒对液流的跟随性越强.在二次流的影响下，这些颗粒集中在γ=-30°～30°区间壁面附近，流出弯管后，二次流的作用减弱，颗粒逐渐呈分散状态.

4.1.2 弯管壁面的磨损特性

通过数值模拟获得的弯管壁面磨损量分布如图10所示.可以看出，整个管路的磨损集中在弯管壁面上，弯管外脊线及附近区域磨损较为严重.进口段直管无磨损产生，出口段直管与弯管交界处仅出现少量磨损.

Fig.10 Distribution of the amount of wear over the inner wall of the pipe图10 壁面磨损量分布

在弯管外脊线即外侧壁面的中心线附近取3条线，分别对应γ=0°、-5°和5°位置，研究外侧面磨损随断面位置θ的变化规律，结果如图11所示.可以看出，侧面上各位置的磨损量随着θ的增大先增大后减小，脊线上(γ=0°)在θ=50°～80°范围内磨损较为严重，在θ=60°处，磨损量达到最大值.γ=-5°和5°上各位置的磨损量较为接近，严重磨损区域的磨损量较γ=0°时减小30%左右，磨损量最大位置分别位于θ=50°和70°处.

Fig.11 Distribution of the amount of wear at different positions图11 各位置磨损量分布

图12所示为试验前后弯管出口处的图片.为便于观察壁面磨损，试验前将壁面喷涂蓝漆.经过40 h的试验，将试验段弯管取下，可以看到弯管内壁的蓝漆都已被磨蚀，呈现出明显的粗糙形态.在数值模拟结果中，脊线位置的磨损量最大.运用高压水射流切割技术，将弯管沿着脊线切开，结果如图13所示.可以看出，弯管内侧壁面已被破坏，表面粗糙，说明此处存在一定程度的磨损.弯管外侧面较为光滑，但弯管外侧壁面的厚度较直管部分明显减小；同时，在弯管出口与直管的交界处出现明显沟槽，这与自弯管出口到直管进口的壁面磨损量陡降有关.

Fig.12 Images of the outlet of the elbow pipe before and after the test图12 试验前后弯管出口壁面的形态

Fig.13 Images of the divided elbow pipe after the test图13 试验后弯管剖开图

为便于将模拟与试验结果进行对比，引入无量纲相对磨损率WR，用以衡量各磨损量的相对大小.相对磨损率WR表示为

其中，W为各位置的磨损深度，WMax为壁面上磨损的最大深度.试验和模拟中弯管脊线上各位置的相对磨损率如图14所示，根据图3中所示的断面划分方法，自θ=0°断面开始，随着θ不断增大，相对磨损率均先增大后减小，试验结果中，在θ=50°时相对磨损率最大；模拟结果显示在θ=60°时磨损最为严重.由于试验中弯管进口接近驱动泵的出口，故弯管进口的流速沿径向分布不均匀；而数值模拟中采用了速度进口边界条件，流速设为均匀分布.流速分布是引起试验与模拟结果差异的主要原因.在试验中，流速较高的液体携带颗粒通过弯管时，由于颗粒所受的离心力较模拟结果大，所以颗粒更趋向于碰撞靠近弯管进口的壁面，从而造成θ=50°位置磨损严重.从弯管出口(θ=90°)到出水管的进口处(θ=100°)，相对磨损率迅速下降，试验和模拟结果较一致.

Fig.14 Distributions of relative wear rate obtained numerically and experimentally图14 试验和模拟获得的相对磨损率分布

颗粒碰撞弯管壁面的接触能量分布如图15所示.可以看出，弯管外侧面受到的切向累积能量远高于法向累积能量.所以弯管壁面磨损产生的主要原因是颗粒对壁面的小角度碰撞和切削[31].提取弯管外侧脊线区域在0.1 s时长内颗粒与壁面的碰撞数据.图16所示为在弯管不同位置颗粒与壁面的碰撞速度和碰撞角度.在θ=0°～20°区域内，颗粒对该区域的碰撞速度较高，但碰撞频率较低，因此该区域无明显磨损.在θ=30°～60°区域内，受外侧壁面低速流体的影响，颗粒的碰撞速度不断下降，而该区域的碰撞角度较大，大多数颗粒与壁面的首次碰撞发生在该区域，在θ=30°位置，颗粒碰撞壁面的角度达到最大值，这与图15中法向能量最大的区域基本吻合.在θ=60°～90°区域内，颗粒在运动过程中与壁面产生碰撞后，还会与其他颗粒产生碰撞，因此沿弯管出口，碰撞角度逐渐减小，颗粒做接近直线的波浪式运动.

Fig.15 Distribution of contact energy over the surface of the elbow pipe图15 弯管表面接触能量分布

Fig.16 Collision velocity and collision angle as particles impact the pipe wall图16 颗粒与壁面的碰撞速度及角度

4.2 颗粒形状对壁面磨损的影响

为探究颗粒形状对弯管壁面磨损的影响，在模拟中设置五种不同工况，每种工况中仅取一种颗粒形状，其他条件与标准工况保持一致.图17所示为不同颗粒形状工况下外脊线上的磨损量分布.弯管脊线上各位置的磨损量总体上随着θ的增大先增大后减小.五种工况中，四面体由于球形度最小，颗粒尖锐，这是四面体颗粒与壁面产生碰撞时会造成更严重磨损的主要原因.二十面体颗粒引起的磨损略高于十二面体.在不同颗粒形状工况下，弯管最大磨损量位置不同，当管内颗粒为四面体时，所对应的最大磨损量位置位于θ=80°的脊线上，当管内颗粒为二十面体时，最大磨损位置位于θ=50°处，随着颗粒球形度的增大，出现最大磨损量的位置更靠近弯管进口.

Fig.17 Effect of particle shape on wear of the elbow pipe图17 颗粒形状对弯管磨损量的影响

随着颗粒球形度的增大，颗粒与流体的接触面积减小，所受到的曳力作用减弱，颗粒的随流性越差.取模拟工况中球形度相差最大的两种颗粒进行分析，图18所示为四面体颗粒和二十面体颗粒在流场中的角速度和速度分布.颗粒进入弯管后，受二次流和其他颗粒的影响，开始做低速旋转运动，颗粒与壁面发生碰撞后反弹，引起受力失衡，从而颗粒以较高的角速度进行旋转.在外侧壁面θ=50°的位置，二十面体的角速度较大，这是由于二十面体颗粒的随流性最差，主流区的颗粒在惯性作用下，与θ=50°的位置发生碰撞，造成50°位置附近的壁面出现较为严重的磨损.相反，四面体颗粒的随流性最好，主流区的颗粒脱离流线的时刻也最晚，因此外侧壁面θ=80°位置的磨损最为严重.颗粒跟随液流的能力不仅影响磨损的位置，也影响磨损的程度[32].由图18可知，四面体颗粒的运动速度大于二十面体颗粒，因此碰撞壁面时的速度较高，造成更严重的磨损.

Fig.18 Angular velocity and velocity distributions of two group of particles with different shapes图18 两种颗粒的角速度及速度分布

5 结论

a.采用DEM和CFD耦合的模拟方法可以准确预测固液两相流体在弯管内的运动特征，所采用的磨损模型可以较为合理的预测弯管壁面磨损的位置与程度.

b.弯管内存在强烈的二次流空间流动结构，颗粒在二次流的影响下，集中于弯管外壁中心线附近，所以弯管外壁面中心区域的磨损较为严重.

c.在多颗粒混合输送工况下，沿着弯管出口方向，壁面磨损量先增大后减小，在θ=60°弯管截面处的磨损量最大，与试验结果较一致.

d.单一颗粒形状条件下，随着颗粒球形度的增大，颗粒的随流性越差，磨损严重区域向弯管进口方向迁移，壁面的磨损程度先减小后增大；在球形度为0.91时，壁面的磨损量最小.