基于歐拉-歐拉方法的細長管道內磨粒流加工數值模擬研究

劉 宵， 尹 洪 超， 穆 林

(大連理工大學 能源與動力學院 遼寧省復雜能源轉化和高效利用重點實驗室，遼寧 大連 116024)

0 引 言

固液兩相流體的流動和相互作用在工業生產過程中非常普遍.如石油、泥漿、雨水在管道中的輸運過程，流體中夾雜的固體顆粒會在管道中沉積、磨削壁面[1]；在傳熱傳質領域，在流體中加入固體顆粒，會提高流體的對流換熱系數，使傳熱過程得到強化[2].20世紀60年代，美國宇航局為去除航天用液壓閥體等裝置的內部毛刺，提出一種利用兩相流拋光零部件表面的加工工藝，即磨粒流加工工藝(abrasive flow machining，AFM)[3].

近年來，AFM工藝不斷完善，表現出拋光效率高、適應性強、對工件損傷低、表面壓應力提高等優勢，并得到了廣泛關注和應用[3].國外研究方面，Williams等[4]以AFM的加工方式對密封板的鋁-鋼接觸部位進行拋光加工，極大提升了金屬板密封、防泄漏的效果.Davies等[5]研究了不同黏度的流體與不同粒徑的顆粒對AFM加工的影響.Petri等[6]考慮了影響磨粒流加工的諸多參數，設計出適合各種流體和被加工材料的AFM預測系統.國內研究方面，周迪鋒等[7]以CFD-DEM耦合模型研究了磨粒撞擊次數和撞擊速度對于拋光效果的影響.李琛[8]以可實現k-ε湍流模型和軟性磨料的方法研究了磨粒流加工過程中速度、壓力、湍流動能、湍流耗散率等諸多流體參數的分布特征.高航等[9]指出，磨粒粒度特性以及不同黏度的流體流動特性對AFM加工過程具有顯著影響，仍需進行深入研究.

另一方面，基于加工機理的仿真建模和數值分析過程可以分為流體流動和顆粒運動兩個方面.其中，以流體流動為研究重點的主要目的是為了得到整個流場內的壓力、速度、湍流強度等分布情況，這是影響磨粒流加工特性的宏觀因素，可以對加工后的結果進行定性預測.此類研究大多選擇了基于拉格朗日方法的Mixture模型[10]或VOF模型[11]，適用于復雜曲面加工、磨粒表面的追蹤以及湍流脈動的研究[12].而以顆粒運動為研究重點，則可以更好地得到局部磨削程度的規律，以及顆粒對壁面磨削量的預測.此類研究大多采用了DPM(discrete phase model)或DEM(discrete element method)，分別適合模擬顆粒軌跡追蹤和顆粒之間的相互作用規律.本文的研究對象為某船用高壓氧輸送管道，該管道是一種狹長的異形管道，適合從流體的角度來研究AFM加工效果.相關領域的學者們大多采用歐拉-拉格朗日的方法進行研究[13]，先求解流體場，再將部分流體參數代入求解固體顆粒的運動，這樣可能會導致固體顆粒的運動參數不準確，因此本文以歐拉-歐拉方法為原理構建仿真計算模型，該方法采用兩套方程分別對流體和固體進行求解，并以相對速度、作用力等進行耦合，避免流、固參數的混淆.同時，考慮AFM加工過程中可能的傳熱及傳質特性對加工過程的影響，重點研究顆粒粒徑、磨料入口溫度、磨料流體物性參數及流動特性等對壁面磨削特性的影響規律.

1 數學模型

數學計算模型以CFD理論為基礎，包括兩相流方程、湍流方程、能量方程.在固液兩相流動的耦合計算方面，采用基于歐拉-歐拉原理的求解方法，將固體相(solid phase)視為流體，與液相(liquid phase)相互滲透，考慮固體顆粒相之間以及固液兩相之間的相互作用.對固體、流體分別采用獨立的方程進行求解，其連續方程和動量方程如下：

(1)

(2)

(3)

(4)

式中：α為體積分數，ρ為密度，t為時間，v為速度，m為傳遞的質量，S為源項，p為壓力，τ為應力應變張量，g為重力加速度，F為外部體積力，下標s、l分別表示固相和液相，下標sl和ls分別表示從固相到液相和從液相到固相傳遞(或作用)關系，lift表示上升，vm表示虛質量.兩相間的相互作用力Rsl和Rls分別見式(5)和(6)：

Rsl=Ksl(vs-vl)

(5)

Rls=Kls(vl-vs)

(6)

式中：Ksl和Kls為兩相之間的動量交換系數，與摩擦、黏滯等效應有關.

式(3)和(4)中，Flift為升力，Fvm為虛質量力，表示固相相對于液相相對運動時，固體顆粒受到的慣性作用，計算公式分別為

Flift=-0.5αsρl(vl-vs)×(×vq)

(7)

(8)

液相和固相方程之間的平衡關系通過Rsl、Rls、Flift、Fvm耦合進行計算.

湍流方程采用基于雷諾平均的Navier-Stokes方程(RANS)，湍流動能方程和耗散率方程如下：

(9)

(10)

式中：k為湍流動能；vi為速度；μ為流體的動力黏度；μt為湍動黏度系數；Gk為由速度梯度產生的湍流動能；ε為湍流耗散率；σk為湍流動能普朗特數，取1；σε為湍流耗散普朗特數，取1.3；C1、C2為經驗常數，分別取1.44和1.92；ν為流體運動黏度.由于本文研究的流體介質具有不同黏度特性，為了更精確求解湍流黏度，式(9)和式(10)中的μt不宜取為常數，而采用下式計算：

(11)

(12)

固相和液相流動過程的能量傳遞按下式計算：

式中：h為比焓，u為熱力學能，q為熱通量，Q為熱交換量.

2 物理模型

2.1 幾何模型

本文以某船用高壓氧輸送管道為研究對象，構建了直管、三通管和U形管3種管型結構的物理模型，其幾何參數見圖1，管道外徑和內徑均取R=18 mm和r=9 mm，L為從流道入口到出口的長度，取L=6 m.管子材料為銅鎳合金，物性參數見表1.

圖1 物理模型

表1 管子物性參數

2.2 模擬設置

以ANASYS18.2系列軟件為模擬工具，首先使用ICEM-CFD進行幾何模型的建立和網格劃分，區分流體計算域及固體計算域，由于模型尺寸較長，為提升計算速度，在長直段的網格劃分稀疏，在接頭部及彎部的網格劃分較為密集.3種管型的網格數量在220×104～260×104.

選擇Fluent軟件為仿真計算工具，開啟Realizablek-ε湍流模型、能量方程、多相流模型.根據流速、顆粒濃度的相關研究可知[8，13]，速度會對AFM加工效果有較大影響，而顆粒濃度的影響較小.為實現較好的磨削效果同時避免顆粒沉積現象的出現，通常選取較高的入口速度以及較高的磨粒濃度[14].因此入口速度設置為30 m·s-1，液相設置為水.固相設置為SiC顆粒，體積分數為10%，為研究不同顆粒粒徑對AFM加工過程的影響，選取了3種顆粒直徑d=10，50，100 μm分別進行模擬.磨粒流的初始溫度分別設置為293、323、363 K.出口邊界條件為流動出口.內壁面的表面粗糙度取50 μm，摩擦因數0.5；外壁面設置為對流換熱表面，對流換熱系數為20 W·m-2·K-1，環境溫度T0為293 K.

2.3 模型驗證

為確保仿真計算結果合理，分別進行網格無關性驗證和數值模型可靠性驗證.對直管幾何模型進行網格加密處理后，得到的壓強分布與未加密時基本吻合，見圖2，表明本文劃分的網格數量能夠滿足計算精度要求.此外，將李俊燁等[15]對噴嘴小孔的磨粒流加工實驗結果與本文方法所建立的噴嘴小孔模型仿真計算結果進行對比，見圖3.結果表明，計算得到的速度分布與實驗結果之間的最大偏差為9.8%，可以認為本文選擇的模型是準確的.

圖2 直管壓力曲線

圖3 不同模型速度曲線的對比

3 結果分析

3.1 粒徑對拋光過程的影響

Preston方程[16]是使用較為廣泛的定量評價材料磨削加工的經驗公式.式中的壓力、速度、時間3個主要變量也是磨粒流加工過程中的重要變量，通過流體參數即可直接描述管壁的磨削效果.其表達式為

(15)

式中：Δz為材料去除量，kP為Preston系數，定義λ=vp，并通過λ的大小反映磨削程度的多少.

不同粒徑的顆粒與壁面發生碰撞時，會有不同的碰撞形態.當顆粒粒徑大于壁面的表面粗糙度時，其碰撞形式為彈性碰撞；顆粒的粒徑小于壁面的表面粗糙度時，顆粒可以進入相鄰凸起之間的縫隙中，除彈性碰撞外，還會產生犁、削的效果[17].除此之外，大顆粒具有較大的動能，單個顆粒與壁面碰撞時的磨削量更多，λ更大.因此，大粒徑的顆粒具有更好的磨削效果，但是在磨削至一定程度后，彈性碰撞的作用效果不再明顯，需更換小粒徑的顆粒進行犁、削，才會使壁面粗糙度進一步降低.3種管型λ分布見圖4.粒徑d=10 μm時，U形管的λ曲線波動較為明顯，這主要是由于U形管的幾何模型特點，流道分成內側(曲率較小一側)和外側(曲率較大一側)，見圖1(c).流體在轉彎處的內、外側會產生較大的速度、壓力差，并會對接下來的流動過程產生擾動.此外，小粒徑顆粒對壁面施加的彈性碰撞和犁、削綜合作用，導致局部壓力、速度波動較大，進而導致局部的λ波動較大.而直管與三通管中的流體均沒有受到因流道轉彎產生的內、外側擾動，因此局部λ變化較為穩定.出口壓力設置為0，入口壓力根據過程的壓降進行耦合計算，因此總體上λ呈下降趨勢，但在入口處，流體進入管內后，其流動會向充分發展段過渡，使得近壁區的速度增加，進而導致λ在入口位置先呈現出增加的趨勢.

圖4 不同管型λ曲線

相關研究[8，10]側重于流體和固體顆粒對壁面的磨削綜合作用效果，因此通常將流體的速度和壓力假設為固體顆粒的速度和壓力.本文采用歐拉-歐拉模型進行模擬，對流體和固體顆粒分別分析了λ的變化趨勢.總體上看，流體與固體顆粒的差異并不明顯.但是在U形管的局部，以及三通管的匯合處，λ差異較為明顯.這是因為相比于流體，固體顆粒具有更大的慣性，同時流體是連續的，固體顆粒是不連續的，因此固體顆粒運動狀態的變化具有一定的滯后性.在穩定的流段中，流體與固體顆粒可以保持相同的運動狀態，但在不穩定的流段中，固體顆粒與流體在經歷轉向、碰撞等過程后，會產生相對的滑移速度，使得局部的運動狀態產生明顯的差異.此外，粒徑的大小也會影響流體和固體顆粒對壁面的作用差異.粒徑越小，其慣性越小，宏觀運動的形式越接近流體，因此流體和固體顆粒的λ差異越小.根據3種管型管內流動的特性可以推測，在更復雜模型的AFM加工過程中，流體與磨粒的流動差異更加明顯，不能以流體的運動參數代替固體的運動參數.

3.2 入口溫度對拋光過程的影響

U形管的流道既有長直部分，又有局部突變部分，而對溫度分布的研究與管型關系不大，因此選擇U形管為研究對象.為避免流體介質出現凝固、沸騰等相變現象，同時又要突出溫度的差異性影響，液相流體水分別選取3個不同的初始溫度Tw=293，323，363 K，表2為3種溫度下水的物理性質.

表2 不同溫度下水的物性參數

由圖5可看出，溫度Tw=293 K時的加工效果最好，這主要是兩方面原因導致的.從液相來看，溫度較低時，流體的密度較高，對管壁可以產生更大的擠壓力；同時，該溫度下流體表面張力較大，對于固體顆粒有更強的包裹性，在近壁區會對壁面產生更大的彈性力.從固體顆粒來看，金屬材料以及SiC顆粒在低溫時，密度也會增大，從而單個顆粒具有的動能更大；此外，低溫時材料的韌性較低而脆性較大，SiC顆粒碰撞壁面的粗糙部位時，產生的沖擊力更易使粗糙部位被磨削掉.

圖5 不同溫度下的λ曲線

AFM加工過程中，流體、固體顆粒與管壁摩擦，會有大量熱量產生.本文模擬加工時間1 h，在外界環境溫度T0=293 K的情況下，Tw與T0的溫差ΔT越小，則磨削過程中的熱積聚量越多.在管子不同L處的溫度分布基本均勻，只在進出口部位有所差異，因此選取了L=0和L=3 m兩個位置的截面進行分析.圖6中，T表示截面各處的溫度.隨著加工時間的延長，管壁溫度會逐漸升高，如前文所述，溫度的升高會對加工產生不利影響.因此在加工過程中，對管壁適當進行冷卻，或采取間斷加工的方式，有利于實現更好的加工效果.

圖6 截面溫度云圖

3.3 流體介質對拋光過程的影響

本節同樣選取U形管為研究對象，對不同流體介質對拋光過程影響進行研究，見圖7，流體介質分別為水和酒精時，λ具有一定的差異，為進一步研究這種差異產生的原因，本文選擇了流體的密度和動力黏度作為研究對象.靜壓和動壓計算公式中都包含密度變量，因此密度變化會對壓強分布造成影響.動力黏度表示流體中發生相對運動所需要的力的大小，從而會對流速、壓力產生影響.為此，補充兩種假想流體進行模擬.流體A的動力黏度與酒精相同，其他物理性質與水相同；流體B的密度與酒精相同，其他物理性質與水相同，4種流體的主要物性參數見表3.

表3 不同流體物性

圖7 不同流體的λ曲線圖

結果表明，流體的密度越大，對管壁的擠壓程度就越大，從而壁面受到的壓力越大.因此密度主要是通過影響壓強大小進而影響λ的大小.流體的動力黏度反映了流體與管壁產生相對運動時的摩擦力大小.流體的動力黏度大，宏觀表現為黏滯性高、流動所受的阻力大，損失的動能用于磨削壁面，即壁面受到了更大的切削力作用，因此當流體的動力黏度升高時，磨削效果會提升.綜合密度與動力黏度影響因素，流體A的動力黏度雖然小于流體B，但密度較比于流體B的更大，因為在本文模擬工況中密度起主導作用，流體A的λ更高.根據結果可以預測到，選擇密度大、動力黏度高的流體作為介質進行AFM加工，可以獲得較好的加工效果.

4 結 論

(1)磨粒粒徑對AFM加工效果有顯著影響.磨粒粒徑大于壁面表面粗糙度時，單位時間內的材料磨削量較大；磨粒粒徑小于壁面表面粗糙度時，U形管道流動的局部λ波動較大.磨粒粒徑增大使得流體相和固體相在流道的突變部位對管壁的作用差異明顯增強.

(2)分析了入口溫度對AFM加工的影響規律，得到了不同入口溫度下的AFM加工過程熱量積聚分布云圖，結果表明降低加工溫度、維持管壁恒溫更有利于加工過程.

(3)研究了不同流體介質對AFM加工效果的影響.以流體的密度和動力黏度為變量，對水、酒精以及兩種假想流體進行模擬，得到了以高密度、高動力黏度的流體作為磨粒流介質更有利于AFM的結論.