非球形固體顆粒對彎管內壁的磨損機制

曹 騫,康 燦*,滕 爽,焦 儂,丁可金

(1.江蘇大學 能源與動力工程學院,江蘇 鎮(zhèn)江 212013;2.中國船舶集團第七〇四研究所,上海 200031)

作為采礦系統(tǒng)的關鍵部件之一，彎管起到連接各直管段的關鍵作用.在礦漿輸送過程中，固體顆粒并不完全跟隨液流運動，尤其在曲率較大的彎管中，顆粒與壁面將發(fā)生碰撞和摩擦，最終導致壁面材料去除，即磨損[1-3].單個顆粒單次作用于壁面而產生的磨損可以忽略，但長時間內大量顆粒作用于同一壁面位置，將造成壁面穿透甚至管路失效，不僅引起經濟損失，還會帶來嚴重的安全隱患[4-6].研究顆粒與管道內壁之間的相互作用具有重要的意義.在實際輸送過程中，礦石顆粒形狀各異，顆粒的形狀無法準確描述，也很難用統(tǒng)一的標準去評價顆粒形狀的不規(guī)則程度.目前，在固液兩相流動的研究中，針對顆粒形狀開展的研究較少，對不同形狀顆粒與管壁之間的摩擦與磨損的認識尚不充分[7].

針對管內固液兩相流動的研究主要從試驗與數值模擬兩個方面開展.Alajbegovic等[8]采用激光多普勒測速儀(LDA)分別測量了豎直管道內顆粒和流體的速度分布，進而分析了顆粒對管內湍流的抑制作用；Kesana等[9]借助探頭在水平放置的彎管內不同位置進行取樣，研究顆粒直徑和流體黏度對顆粒分布的影響.歐拉-歐拉和歐拉-拉格朗日方法是處理固液兩相流動時最常用的兩種數值模擬方法，前者將離散相處理為“擬流體”，在歐拉框架下對顆粒進行求解，該方法適用于顆粒濃度較高而直徑較小的情況；后者采用N-S方程(即Navier-Stokes方程，用以描述黏性不可壓縮流體的動量守恒)在歐拉框架下對連續(xù)相進行求解，同時運用牛頓第二定律在拉格朗日框架下對離散相進行求解[10].

早期關于顆粒與管壁之間摩擦及壁面磨損的研究主要依靠試驗，根據試驗結果建立經驗及半經驗的磨損預測模型，描述兩相流動狀態(tài)對壁面磨損的影響[11-12].Archard模型是最早建立、應用最為廣泛的磨損模型之一，該模型認為磨損量與顆粒和壁面接觸時的法向力和滑動距離成正比，與材料的硬度成反比[13].Finnie等[14-15]首次提出塑性材料的微切削理論，建立了考慮碰撞速度及碰撞角相關的磨損模型；Oka等[16-17]討論了材料類型、沖擊角、沖擊速度以及顆粒尺寸對磨損的影響，建立了Oka模型；在Zhang的模型中，額外考慮了顆粒形狀對壁面磨損的影響[18].隨著數值模擬技術的快速發(fā)展，采用磨損模型配合多相流模型來預測流體機械過流部件的磨損已成為重要方法之一[19-20].

本文中采用DEM與CFD耦合的方法，對一彎管內的固液兩相流動及固體顆粒對管壁的磨損進行數值模擬，借助商用離散元軟件EDEM和計算流體動力學軟件ANSYS Fluent，并充分考慮顆粒與顆粒、顆粒與流體以及顆粒與壁面之間的相互作用，模擬彎管內的流動特征、顆粒的運動規(guī)律及彎管壁面的定量磨損規(guī)律；通過與試驗結果進行對比，驗證模擬的有效性；進而改變顆粒形狀，研究顆粒形狀對顆粒運動及彎管內壁磨損的影響.

1 彎管磨損試驗方法

彎管磨損試驗在如圖1所示的磨損試驗臺上進行.該試驗臺主要由料箱、閘閥、驅動泵及試驗段構成.閘閥1和閘閥2分別用于控制清水的流量和顆粒的濃度，在試驗段的進口和出口各安裝1個壓力表以監(jiān)測管路壓力.固液兩相混合物在離心泵的驅動作用下進入彎管.試驗時所采用顆粒的形狀不規(guī)則，取一定數量的顆粒，通過測量顆粒的密度和質量，計算出顆粒的當量直徑約為4.5 mm，顆粒相的體積分數約為5%.試驗中采用內徑為40 mm，彎徑比(彎管中軸線曲率半徑與管徑之比)為1.5的彎管進行磨損試驗.

Fig.1 Test bench for wear圖1 磨損試驗臺

2 數值模型及模擬方案

2.1 流動控制方程

本研究中采用雙向耦合方法來處理顆粒與液相間的相互作用力，所有計算域的流體相均在絕對坐標系中進行求解.流體流動的控制方程如下：

連續(xù)性方程：

動量守恒方程：

其中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的比例構建.為保證彎管入口流動充分發(fā)展、顆粒分布均勻以及出口無回流，進口直管和出口的長度均設為管徑(D)的10倍.圖中θ為極角，定義彎管進口截面θ=0°，出口截面位于θ=90°位置.γ用于表征彎管壁面上點的位置，彎管外側脊線γ=0°，內側脊線γ=180° (-180°)，如圖3所示.

Fig.3 Computational model圖3 計算模型

采用六面體結構化網格對計算域進行離散，以保證計算的穩(wěn)定性以及計算結果的準確性.近壁面網格的y+分布在30～300之間.為驗證網格數對計算結果的影響，設計了網格數不同的5套計算域網格，在相同的邊界及工況條件下進行計算，比較不同網格數方案獲得的彎管進出口壓差，結果列于表1中.隨著網格數的增加，彎管進出口壓差的相對變化量減小.最終選擇了網格數約為53萬的網格方案(方案4).

表1 網格無關性驗證Table 1 Grid independence examination

2.4 計算方法及邊界條件設置

固液兩相之間的耦合計算在歐拉-拉格朗日框架下進行.利用ANSYS Fluent在歐拉坐標系中對液相進行求解，控制方程為基于雷諾平均的Navier-Stokes方程，采用Realizablek-ε湍流模型，通過設置入口湍動能(k)和湍動能耗散率(ε)使控制方程封閉，采用SIMPLEC算法耦合流場中的壓力和速度.計算域的進口采用速度進口邊界條件，出口設置為Outflow邊界.近壁面區(qū)采用標準壁面函數處理，壁面粗糙度設為0.046 mm.利用EDEM在拉格朗日坐標系中對固相進行求解.顆粒與顆粒和顆粒與壁面之間的接觸模型采用Hertz-Mindlin無滑移模型[20].顆粒與壁面的屬性列于表2中.通過耦合接口，實現EDEM和ANSYS Fluent之間的顆粒位置和動量等信息的傳遞[25].

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

2.5 磨損模型

目前大多數磨損模型基于經驗公式建立，無法將所有磨損因素考慮在內.本文中選用Zhang開發(fā)的用于計算沖蝕的經驗模型，該模型認為壁面磨損主要與顆粒形狀、碰撞速度及碰撞角度有關[18,26].該模型的表達式為式中，ER表示侵蝕率，為單位質量的顆粒碰撞壁面引起的材料損失；BH為壁面材料的布氏硬度；Fs為顆粒的形狀因子，對于尖銳顆粒，Fs=1.0，對于球形顆粒，Fs=0.2；Vp為顆粒碰撞壁面時的速度；α為用弧度表示的顆粒碰撞角度；經驗系數n=2.41，C=2.17×10-7.

本文中通過編程開發(fā)軟件Visual Studio完成模型編譯，并通過API (Application Programming Interface)嵌入EDEM，從而獲得顆粒與壁面的碰撞以及壁面磨損信息.

在EDEM中，通過單位面積的磨損深度(h)來表征壁面磨損情況：

其中m為與壁面碰撞的顆粒的質量，A為計算單元的面積，ρw為壁面材料的密度.

EDEM根據顆粒與壁面間的相對速度和力識別磨損區(qū)域的受力.切向累積能對應顆粒在壁面滑動以及小角度碰撞所產生的能量，法向累積能對應顆粒以大角度碰撞壁面產生的能量.切向累積接觸能(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為管道半徑.可以看出，在主流區(qū)內，流體速度明顯高于顆粒速度，這是由于相間滑移所致，數值模擬結果與試驗結果表現出的趨勢一致.沿徑向方向，顆粒和流體的速度均逐漸下降；受邊界層的影響，靠近壁面的流體和顆粒的速度出現較大的變化梯度.試驗結果中各個位置的顆粒速度略高于模擬結果中的對應值.對于液相速度，試驗結果與模擬結果接近.綜合考慮試驗誤差及數值模型中無法考慮的實際因素，本文中的數值模型有效.

Fig.5 Comparison between numerical and test results圖5 模擬和試驗結果對比

4 結果分析及討論

4.1 彎管內流動及磨損分析

試驗中采用的固體顆粒形狀不一，平均當量直徑為4.5 mm.為了使模擬工況接近試驗工況，模擬中的固體顆粒由圖3中五種不同形狀系數的顆粒構成，每種顆粒所占的固相體積分數為1%.在輸送的固液兩相混合物中，固相體積分數為5%.

圖6所示為通過模擬得到管道內顆粒的質量(ms)隨啟動時間(t)的變化，即在管道內的固液兩相混合物開始流動的初始階段所獲得的結果.可以看出，自初始時刻開始，管道內的顆粒數目持續(xù)增加，在t=0.14 s時刻，管道內固體顆粒的總質量達到最大，而后管道內的顆粒質量達到動態(tài)平衡，流動狀態(tài)趨于穩(wěn)定.為了獲得穩(wěn)定的磨損模擬結果，對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)中可以明顯看出，流體在彎管內側出現局部加速現象，而在流經彎管后出現低速區(qū)，斷面流速分布不均勻.從圖7(b)的壓力分布來看，其基本上與圖7(a)所示的速度分布相呼應.

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

流體在經過90°彎管時，在離心力的作用下，出現強烈的二次流，這從圖7中的流動速度分布可以看出，該現象也與文獻[27]中獲得的結果相似.彎管內部的流體向曲率較大的外壁流動，由于流體流動具有連續(xù)性，外壁附近的流體被迫向彎管內部流動，形成復雜的空間流動結構[28].

沿彎管進口(0°截面)位置向彎管出口(90°截面)方向創(chuàng)建角度間隔為30°的截面.在彎管進口位置，流體受到離心力的作用較小，不存在二次流動，此處指向流線曲率中心的順壓梯度較大，流體的徑向速度指向內側壁面，如圖8所示.自30°截面位置開始，過流斷面上出現旋向相反的對稱渦，此時顆粒對流體的擾動較強，渦核向彎管內側壁面靠攏，低速運動的液體被旋渦卷吸，并隨著主流向下游發(fā)展[29].從旋渦形態(tài)的演變來看，渦核的位置不斷向內側壁面移動，同時渦的強度降低，而旋渦的形狀趨于平滑和規(guī)則，對稱性得到保持.同時，在流體從進口向出口流動的過程中，受到的離心力作用不斷累積，使得內側壁面周圍的流體速度不斷降低，同時外側壁面周圍流體速度不斷升高，如圖8所示.

Fig.8 Patterns of liquid flow on different cross sections圖8 不同截面的液流形態(tài)

顆粒的軌跡及運動矢量分布如圖9所示，顆粒與流體以相同的速度自進口均勻進入計算域內.由于邊界層的影響，靠近管壁的顆粒在上升過程中較主流區(qū)顆粒的速度低[30].顆粒進入彎管后，靠近外側壁面的顆粒速度逐漸下降；同時，位于主流區(qū)的顆粒在慣性作用下繼續(xù)向下游運動，受流場的影響速度下降，在與彎管外側壁面發(fā)生第一次碰撞后，速度顯著降低.顆粒速度矢量箭頭的長度表征顆粒所受曳力的大小，彎管內側壁面附近的顆粒自高速液流中獲得能量，與彎管外側壁面產生碰撞.在γ=50°～100°截面區(qū)間內，顆粒速度最低.顆粒離開彎管后，速度逐漸上升.

Fig.9 Particle trajectory and velocity vector distribution圖9 顆粒軌跡及運動矢量分布

部分顆粒與γ=±(45°～60°)區(qū)間壁面發(fā)生碰撞后，速度降低，與流體之間的速度差增大，所受到的曳力作用更強.顆粒的Stokes數表征顆粒慣性力與曳力的相對大小，顆粒的Stokes數越小，則顆粒對液流的跟隨性越強.在二次流的影響下，這些顆粒集中在γ=-30°～30°區(qū)間壁面附近，流出彎管后，二次流的作用減弱，顆粒逐漸呈分散狀態(tài).

4.1.2 彎管壁面的磨損特性

通過數值模擬獲得的彎管壁面磨損量分布如圖10所示.可以看出，整個管路的磨損集中在彎管壁面上，彎管外脊線及附近區(qū)域磨損較為嚴重.進口段直管無磨損產生，出口段直管與彎管交界處僅出現少量磨損.

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

在彎管外脊線即外側壁面的中心線附近取3條線，分別對應γ=0°、-5°和5°位置，研究外側面磨損隨斷面位置θ的變化規(guī)律，結果如圖11所示.可以看出，側面上各位置的磨損量隨著θ的增大先增大后減小，脊線上(γ=0°)在θ=50°～80°范圍內磨損較為嚴重，在θ=60°處，磨損量達到最大值.γ=-5°和5°上各位置的磨損量較為接近，嚴重磨損區(qū)域的磨損量較γ=0°時減小30%左右，磨損量最大位置分別位于θ=50°和70°處.

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

圖12所示為試驗前后彎管出口處的圖片.為便于觀察壁面磨損，試驗前將壁面噴涂藍漆.經過40 h的試驗，將試驗段彎管取下，可以看到彎管內壁的藍漆都已被磨蝕，呈現出明顯的粗糙形態(tài).在數值模擬結果中，脊線位置的磨損量最大.運用高壓水射流切割技術，將彎管沿著脊線切開，結果如圖13所示.可以看出，彎管內側壁面已被破壞，表面粗糙，說明此處存在一定程度的磨損.彎管外側面較為光滑，但彎管外側壁面的厚度較直管部分明顯減?。煌瑫r，在彎管出口與直管的交界處出現明顯溝槽，這與自彎管出口到直管進口的壁面磨損量陡降有關.

Fig.12 Images of the outlet of the elbow pipe before and after the test圖12 試驗前后彎管出口壁面的形態(tài)

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].提取彎管外側脊線區(qū)域在0.1 s時長內顆粒與壁面的碰撞數據.圖16所示為在彎管不同位置顆粒與壁面的碰撞速度和碰撞角度.在θ=0°～20°區(qū)域內，顆粒對該區(qū)域的碰撞速度較高，但碰撞頻率較低，因此該區(qū)域無明顯磨損.在θ=30°～60°區(qū)域內，受外側壁面低速流體的影響，顆粒的碰撞速度不斷下降，而該區(qū)域的碰撞角度較大，大多數顆粒與壁面的首次碰撞發(fā)生在該區(qū)域，在θ=30°位置，顆粒碰撞壁面的角度達到最大值，這與圖15中法向能量最大的區(qū)域基本吻合.在θ=60°～90°區(qū)域內，顆粒在運動過程中與壁面產生碰撞后，還會與其他顆粒產生碰撞，因此沿彎管出口，碰撞角度逐漸減小，顆粒做接近直線的波浪式運動.

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所示為四面體顆粒和二十面體顆粒在流場中的角速度和速度分布.顆粒進入彎管后，受二次流和其他顆粒的影響，開始做低速旋轉運動，顆粒與壁面發(fā)生碰撞后反彈，引起受力失衡，從而顆粒以較高的角速度進行旋轉.在外側壁面θ=50°的位置，二十面體的角速度較大，這是由于二十面體顆粒的隨流性最差，主流區(qū)的顆粒在慣性作用下，與θ=50°的位置發(fā)生碰撞，造成50°位置附近的壁面出現較為嚴重的磨損.相反，四面體顆粒的隨流性最好，主流區(qū)的顆粒脫離流線的時刻也最晚，因此外側壁面θ=80°位置的磨損最為嚴重.顆粒跟隨液流的能力不僅影響磨損的位置，也影響磨損的程度[32].由圖18可知，四面體顆粒的運動速度大于二十面體顆粒，因此碰撞壁面時的速度較高，造成更嚴重的磨損.

Fig.18 Angular velocity and velocity distributions of two group of particles with different shapes圖18 兩種顆粒的角速度及速度分布

5 結論

a.采用DEM和CFD耦合的模擬方法可以準確預測固液兩相流體在彎管內的運動特征，所采用的磨損模型可以較為合理的預測彎管壁面磨損的位置與程度.

b.彎管內存在強烈的二次流空間流動結構，顆粒在二次流的影響下，集中于彎管外壁中心線附近，所以彎管外壁面中心區(qū)域的磨損較為嚴重.

c.在多顆?；旌陷斔凸r下，沿著彎管出口方向，壁面磨損量先增大后減小，在θ=60°彎管截面處的磨損量最大，與試驗結果較一致.

d.單一顆粒形狀條件下，隨著顆粒球形度的增大，顆粒的隨流性越差，磨損嚴重區(qū)域向彎管進口方向遷移，壁面的磨損程度先減小后增大；在球形度為0.91時，壁面的磨損量最小.