膏體料漿管道輸送中粗骨料顆粒運動規律分析

顏丙恒，李翠平，吳愛祥，王洪江，侯賀子

(1.北京科技大學 土木與資源工程學院，北京，100083；2.金屬礦山高效開采與安全教育部重點實驗室，北京，100083)

金屬礦山全尾砂膏體料漿作為一種高濃度復雜懸浮液其尾砂粒徑分布范圍廣，在添加粗骨料的條件下，管道輸送時顆粒-流體相互作用更為復雜，輸送過程中顆粒運動問題難以分析[1]。目前高濃度懸浮液中顆粒對膏體料漿的流動影響，主要集中在顆粒質量分數、顆粒級配方案對料漿流變參數的影響上[2-3]。在粗尾砂顆粒以及粗骨料顆粒的離析問題研究中，主要開展了靜置條件下的離析實驗研究[4]，而管道輸送中粗尾砂顆粒以及粗骨料顆粒受剪所引起的離析問題則沒有涉足。但已有研究發現，2 mm級別的粗尾砂顆粒，在屈服應力較小的高濃度懸浮尾砂料漿輸送中，因管道輸送過程對料漿具有一定的剪切變稀作用，會有一定的沉降現象[5-6]；但當料漿屈服應力超過一定值時，則對粗尾砂顆粒在管流輸送中的沉降起到明顯限制[7]。與直徑小于20 μm的尾砂顆粒構成的勻質流動模型相比，尾砂粒徑分布較廣的高濃度懸浮料漿，因粗尾砂顆粒受料漿屈服應力的限制在管道輸送過程中可不出現沉降，此狀態可視為偽勻質流動模型來研究[7-8]。現實生產中為了提高充填體的強度以及改善膏體的流動性可在膏體中添加粒徑更大的粗骨料顆粒，例如：水淬渣、碎石[1]這類較大直徑的固體顆粒。目前工業應用中粗骨料顆粒的上限直徑可以達到 15 mm的級別[9-10]。由于粗骨料顆粒不具有極細尾砂顆粒的流變活性，管道輸送中與全尾砂料漿構成的偽勻質體之間的力學相互作用復雜，粗骨料是否發生相對運動以及運動的規律尚未明晰。在已有相關實驗中，雖然采用了核磁共振成像技術(MRI)與電阻層析成像技術(ERT)[5-7]進行研究，但是其研究的料漿濃度低于全尾砂膏體充填料漿濃度，并且其研究的內容是尾砂料漿管道輸送時固體物料在管道截面上的濃度分布以及料漿流動模型。至于充填料漿中粗骨料顆粒直徑級別的顆粒濃度分布以及輸送過程中的沉降問題則沒有涉及。考慮到全尾砂膏體料漿的不透明性，在管道輸送中難以開展類似于研究黏塑性流體中顆粒運動的光學觀測實驗[11-12]，而采用MRI與ERT觀測管道雖然在技術原理上可行，但是實驗花費成本高，研究使用的實驗設備需要專門定制。為此，本文作者將全尾砂料漿構成的偽勻質體視為載流體，而粗骨料視為被承載固體顆粒。通過分析偽勻質體在管道內的流變特性，將含粗骨料膏體的管道輸送模型視為一種復合流動模型。以數值計算的研究方法來代替難以開展的物理實驗，通過構建模型進行數值模擬，分析粗骨料顆粒在管流過程中的相對運動規律，驗證含粗骨料顆粒的膏體充填料漿管道輸送模型并證明管道不同區域內粗骨料顆粒運動狀態的不同。

1 復合流動模型理論分析

1.1 全尾砂膏體料漿流變學模型

全尾砂膏體料漿濃度高、尾砂粒徑分布范圍廣，相比于低濃度全尾砂料漿在層流狀態輸送時粗尾砂顆粒沉降形成的分層流動模型[7]，高濃度的全尾砂料漿因濃度高、屈服應力大可使較大粒徑的尾砂顆粒在管道輸送過程中處于懸浮狀態而不會發生明顯的沉降現象。此時可以將高濃度全尾砂料漿視為一種偽勻質的懸浮液模型[13]。對于這種高濃度懸浮液，目前主要將其視為具有黏性、塑性、彈性、觸變性的復雜非牛頓流體來研究[1,3,14-15]。在管道輸送問題研究中，彈性現象的時間尺度較小[7]，依據Weissenberg數[16]判定其彈性效應可以不考慮。觸變性是一種發生在受剪流動過程中的動態現象，本文主要研究管道輸送過程中全尾砂膏體料漿流動問題，其流變屬性研究的前提是假設剪切流動達到了動態平衡的階段。本文研究粗骨料顆粒相對運動時，視為載流體的全尾砂膏體料漿其主要的非牛頓流變屬性為黏性與塑性。對此最常見的流變模型為Bingham流體和Herschel-Bulkley流體這2種模型。其流變模型分別如式(1)和式(2)所示：

式中：τ為剪切應力，Pa；τy為屈服應力，Pa；為剪切速率，s-1；ηc為塑性黏度，Pa·s；k為稠度系數，Pa·sn；n為冪律指數，n≤1。不失一般性，本文研究的全尾砂膏體料漿流變學模型采用2參數的Bingham流體模型，其在管道輸送時管內流速分布如圖1所示。圖1中：r為距離管道中心的半徑，m；R為管道的半徑，m；r/R表示相對距離管道中心的位置。流速與剪切變形速率曲線由屈服應力為80 Pa、塑性黏度為1.5 Pa·s計算所得。

從圖1可以看出：Bingham流體的管內輸送流動分為柱塞流動區域與剪切流動區域，柱塞流動區域流速基本相等，而剪切流動區域流速沿管徑方向逐漸降低。Bingham流體在不可壓縮流動的條件下，剪切速率與剪切變形張量的第二不變量相聯系。剪切速率的分布類比于剪切變形速率的分布，其值在中間柱塞流動區基本為0 s-1，沿管徑方向逐漸增大，在管壁處達到最大值，此時對應的壁面剪切應力為τw。

圖1 Bingham流體管道輸送流速分布圖Fig.1 Velocity profile of Bingham fluid in pipeline transportation

1.2 粗骨料膏體復合流動模型分析

依據Bingham流體在管內輸送時流速分布的情況可知：柱塞流動區域內的剪切速率非常低，而在靠近管壁的過程中逐漸增大。全尾砂顆粒與水構成的偽勻質懸浮液在管道輸送時，依據是否發生明顯的剪切流動，沿管道徑向可劃分為剪切流動區與非剪切流動區。

非剪切流動區域內顆粒周圍流場的流速基本不發生變化，流速梯度較小(見圖1中的柱塞流動區)。較小的流速梯度難以使粗骨料顆粒發生轉動而使其與偽勻質懸浮體發生相對運動。若顆粒在此區域滿足所受重力大于懸浮液在垂直方向上屈服應力對其產生的阻力，則顆粒會發生一定的沉降運動。

剪切流動區域內顆粒周圍的流速變化明顯，流速梯度較大(見圖1中的剪切流動區域)。顆粒兩側較大的流速梯度易使此區域內的粗骨料顆粒發生自旋運動，而使其與偽勻質懸浮體發生相對運動。若非剪切流動區范圍較大，則運動的粗骨料顆粒較多。由于靠近管壁粗骨料顆粒不受偽勻質懸浮體的限制，在輸送過程中的相對運動會使其與管壁存在碰撞與摩擦，加劇管壁的磨損。

對上述不同區域內粗骨料顆粒的相對運動進行定量研究時，需要分析粗骨料顆粒與偽勻質懸浮體的相互作用力關系。考慮到偽勻質懸浮體的Bingham流變屬性，粗骨料顆粒作為宏觀體積不可忽略的被承載固體顆粒，其與載流體之間的阻力、升力、虛質量力等相互作用力較為復雜[17]。在復合流動過程中難以定量化研究二者之間的相對運動，因此，使用數值計算的方式分析粗骨料顆粒在管輸過程中的運動規律是一種可行的方法。

1.3 數值計算方法

針對前面的復合流動模型分析，采用 Euler-Lagrangian方法研究粗骨料顆粒與偽勻質懸浮液之間的相對運動。在 Euler坐標系下處理連續相(Bingham流體)，在 Lagrangian坐標系下處理離散相(粗骨料顆粒)。粗骨料顆粒與Bingham流體復合流動時，使用傳統的 DPM 離散相追蹤方法難以確定所附加的作用力公式，故采用宏觀顆粒模型 MPM[18]確定所追蹤顆粒的位移與速度等物理量，基于流體動量變化率來求解粗骨料顆粒與Bingham流體之間的相互作用力[19]。

1.3.1 偽勻質懸浮液數值模擬方法

非牛頓流體數值計算的關鍵是確定剪切速率與表觀黏度之間的關系。本文基于ANSYS Fluent軟件，使用Herschel-Bulkley流變模型來描述Bingham塑性體，設定冪律指數n=1。針對Bingham塑性體具有屈服應力的特點，在求解表觀黏度時，為避免當剪切速率為 0 s-1時出現無窮大的情況，引入臨界剪切速率，其計算公式為：

式中：η為表觀黏度，Pa·s；為臨界剪切速率，s-1；剪切速率大于臨界剪切速率時，表觀黏度由式(3)計算，當剪切速率小于臨界剪切速率時，表觀黏度由式(4)計算。臨界剪切速率一般受非牛頓流體自身流變特性的影響，可通過靜態屈服應力測量實驗獲取[14-15]。

1.3.2 粗骨料顆粒-偽勻質懸浮液數值模擬方法

粗骨料顆粒的直徑往往大于多個流體計算域的網格尺寸。傳統的 DPM 方法在追蹤離散相顆粒時，視顆粒為質量點，不考慮顆粒之間的相互作用問題。對于顆粒與流體間的阻力與升力，需要指定特殊的模型來研究[19]。考慮到粗骨料顆粒直徑約為1 cm數量級，管道直徑為10 cm數量級，顆粒一般包含多個流體計算單元，難以采用 DPM顆粒追蹤方法研究。而宏觀顆粒模型MPM在研究顆粒時，允許顆粒包含多個計算單元，基于顆粒與流體間動量交換計算顆粒與流體間的阻力、升力以及扭矩[18]，同時所研究的流體介質可以為非牛頓流體[20]。計算時MPM模型首先通過顆粒與流體間的動量轉移以確定顆粒周邊的流場速度，之后求解顆粒所受到的一系列作用力。MPM 模型中單位質量下粗骨料顆粒受力的平衡方程[21]為

式中：vp為顆粒的速度，m/s；t為時間，s；Fbody為顆粒所受的體力項(例如：重力)，N；Fsurf為顆粒所受的面力項(例如：阻力、升力與扭矩等)，N；Fcoll為顆粒-顆粒，顆粒-管壁的碰撞力，N。在確定顆粒所受作用力后，便可在下一個計算時間步長內確定顆粒新的速度與位移。其中顆粒所受的面力Fsurf(阻力、升力以及扭矩)由顆粒周圍非牛頓懸浮液的速度場、壓力場以及剪切應力場計算出。對應于速度場、壓力場以及剪切應力場，流體作用于顆粒上的力與扭矩可通過3個分量：虛質量分量、壓力分量以及黏性分量來表示[22]，如下式所示：

式中：i為方向下標，代表在i方向上力的合成；Ri為i方向上粗顆粒所受流場的作用力與扭矩，N；Rm,i為i方向上的虛質量項分量，N；RP,i為i方向上的壓力項分量，N；Rυ,i為i方向上的黏性項分量，N。這3項分量分別由懸浮液流場中的速度場、壓力場和剪切應力場求得，具體計算公式分別如式(7)，(8)和(9)所示。虛質量項分量由顆粒所占體積內所有流體單元的動量改變速率積分求出，如下式所示：

式中：mf為計算單元的流體質量，kg；vf,i為流體在i方向上的速度分量，m/s；vp,i為顆粒在i方向上的速度分量，m/s；Δt為流場的計算時間步長，s。壓力項分量由顆粒表面附近的壓力場計算，如下式所示：

式中：p為壓力，N；d2為顆粒表面附近流體單元的近似面積，m2；r為從流體單元中心指向顆粒中心的徑矢；ri為徑矢r在i方向的直角坐標分量，m。黏性項分量由顆粒表面附近的剪切應力場計算，如下式所示：

式中：τji為作用在與j方向垂直的平面上并沿i正方向的剪切應力；rj為徑矢r在j方向的直角坐標分量，m。求得上述3項分量后，MPM模型將與上述計算值相等的源項添加至流體的相應控制方程中。同時 MPM模型采用硬球碰撞模型求解粗顆粒之間的相互碰撞以及粗顆粒與管壁之間的碰撞。

2 數值模型建立與無關性驗證

2.1 數值模型建立

2.1.1 管道計算域尺寸參數設置

本次數值實驗所構建的計算域為最基本的直管道模型。采用MPM模型時，顆粒直徑與網格尺寸的比值，對計算結果影響較大。考慮到數值計算的可行性，選取了3種不同的顆粒直徑與網格尺寸的比值，分別為3:1，5:1和8:1，構建了3種管道計算域模型。3種管道計算域均由六面體結構網格生成，相應的尺寸參數見表1。

表1 管道尺寸參數表Table 1 Parameters of pipe size

2.1.2 設定模擬參數

數值模擬設定的管道平均流速為1 m/s，Bingham流體對應的流變參數如下：屈服應力為80 Pa，塑性黏度為1.5 Pa·s，臨界剪切速率0.01 s-1，偽勻質懸浮液的混合密度為1 800 kg/m3。對視為Bingham流體的偽勻質懸浮液，在管輸過程中沿徑向方向上剪切應力分布不同，導致相應的表觀黏度隨管徑位置發生變化。判別牛頓流體在管道內流動狀態的雷諾數公式不適應于非牛頓流體的管道流態判別[23]。對Bingham流體而言，可使用有效黏度求解其在管道流動時的雷諾數。其定義為管壁處的剪切應力與管道內平均剪切速率之比，如下式所示：

式中：ηe為有效黏度，Pa·s；U為管道內平均流速，m/s；D為管道內徑，m。

管道壁面處的剪切應力可由 Buckingham方程求出，如下式所示：式中：右邊括號內第3項為高階項，當τy/τw=0.5時忽略掉此項后方程求解相對誤差為 5.9%，當τy/τw=0.4時忽略掉此項后方程求解相對誤差為1.8%，因此，在略去第3項后，仍可以得到足夠精確的數值解[24]。由式(10)與式(11)可以求得，Bingham流體的有效黏度如式(12)所示，相應的雷諾數如式(13)所示。

式中：Re為Bingham流體的雷諾數；ρf為密度，kg/m3。將Bingham流體的各項參數以及管道尺寸參數代入式(12)和(13)，求得在屈服應力為80 Pa，塑性黏度為1.5 Pa·s的條件下，管道內的雷諾數為63.53遠遠小于下臨界雷諾數 2 320[25]。管道內的偽勻質懸浮液在剪切流動過程中，其流動狀態可以視為層流。

2.2 無關性驗證

2.2.1 網格無關性驗證

對不同的管道計算域模型進行了網格無關性驗證。數值實驗時為避免管道入口處流速場的變化影響顆粒的運行軌跡，所加載的顆粒位置距離管道入口邊界0.2 m，并位于管道的非剪切流動區中心。粗骨料顆粒直徑為15 mm，密度為2 700 kg/m3。計算時間步長為10-5s，計算時長為0.5 s。求取3種不同顆粒與網格尺寸比下，計算時長為0.5 s時，X，Y和Z3個坐標軸方向的位移分量，相關數據如表2所示。

表2 網格無關性驗證Table 2 Independence verification of mesh size

由表2可知：顆粒與網格尺寸比為5:1時，顆粒在X，Y和Z這3個坐標軸方向上的位移分量與8:1時的結果相差不大。使用MPM模型時顆粒與網格尺寸比越大，顆粒追蹤結果越精確，但計算成本越高。受制于所研究問題的實際情況以及計算成本問題，本文所選擇的顆粒與網格尺寸比為5:1。

2.2.2 時間步長無關性驗證

顆粒追蹤的結果除了受顆粒與網格尺寸比的影響外，還受到流場計算時間步長的影響。在確定合適的時間步長時需要綜合考慮顆粒的弛豫時間、庫朗數以及計算成本這3個方面。MPM模型建議時間步長應小于其顆粒弛豫時間的1%左右[19]。綜上所述選取了3種不同的時間步長來進行驗證，計算時長0.5 s，X，Y和Z3個坐標軸方向的位移分量如表3所示。

表3 時間步長無關性驗證Table 3 Independence verification of time step

因加載顆粒位于管道軸心，Z軸方向上顆粒兩側所受作用力基本為對稱分布。在0.5 s內，Z軸方向上的位移非常小基本可以不考慮。由表3可知：時間步長取為1×10-5s時，顆粒在X，Y和Z3個軸方向上的位移分量與 5×10-6s時的結果相比雖有一定的出入，但是相差不是很大。考慮到本文意在分析粗骨料顆粒的運動規律，時間步長取為1×10-5s的計算精度滿足本文需要，并且具有可行性。

3 粗骨料顆粒運動規律分析

3.1 不同流態區內粗骨料顆粒位移與線速度分析

類似于上文中無關性驗證，本節在距入口邊界0.2 m處沿Z軸方向與Y軸方向間隔0.01 m加載一個粗骨料顆粒，具有相同編號的粗骨料顆粒在同一個計算文件中，顆粒加載如圖2所示。

圖2 顆粒加載位置及所屬計算案例文件示意圖Fig.2 Diagram of location of particles for calculation cases

3.1.1 粗骨料顆粒沿X軸方向位移與線速度分析

粗骨料顆粒沿Z軸方向分布，計算時長為 0.5 s時，顆粒沿X軸方向(管道輸送軸向)的位移以及相應X軸線速度分量如圖3所示。

圖3 粗骨料顆粒沿X軸方向相對位移與線速度分量Fig.3 Relative displacement and linear velocity component of coarse aggregate particles alongX-axis

由圖3可知：沿Z軸方向，0，0.01和0.02 m處所加載的粗骨料顆粒其沿X軸相對位移基本一致，并且X軸線速度分量均為1.44 m/s左右，這與圖1中柱塞流區的流速基本一致。而0.03 m與0.04 m處加載的粗骨料顆粒其沿X軸相對位移以及線速度分量與前3個位置相比差異較明顯，并且越靠近管壁差異更為明顯。可見粗骨料顆粒在剪切流動區內由于兩側流速差的作用在X軸方向(管道軸向)上與偽勻質體之間存在一定的相對運動。

3.1.2 粗骨料顆粒沿Z軸方向位移與線速度分析

粗骨料顆粒沿Z軸方向分布，計算時長0.5 s時，顆粒沿Z軸方向(管徑方向)的位移以及相應Z軸線速度分量如圖4所示。

圖4 粗骨料顆粒沿Z軸方向相對位移與線速度分量Fig.4 Relative displacement and linear velocity component of coarse aggregate particles alongZ-axis

由圖4可知：沿Z軸方向，0，0.01和0.02 m處所加載的粗骨料顆粒，沿Z軸相對位移與線速度值基本一致。粗骨料顆粒的Z軸相對位移以及Z軸線速度分量均非常小可以忽略不計，說明非剪切流動區域內粗骨料顆粒在輸送過程中沿管道徑向基本不存在相對運動，被偽勻質懸浮液很好地限制。而在剪切流動區域內粗骨料顆粒的Z軸相對位移與Z軸線速度分量較明顯，因此，在輸送過程中與偽勻質懸浮液之間存在相對運動。0.04 m處的顆粒由于靠近管壁，受到管壁限制作用其值小于0.03 m處的顆粒徑向位移。

3.1.3 粗骨料顆粒沿Y軸方向位移與線速度分析

粗骨料顆粒沿Y軸方向分布，計算時長0.5 s，顆粒沿Y軸方向(豎直方向)的位移以及相應Y軸線速度分量如圖5所示。

圖5 粗骨料顆粒沿Y軸方向相對位移與線速度分量Fig.5 Relative displacement and linear velocity component of coarse aggregate particles alongY-axis

由圖5可知：不同Y軸坐標的粗骨料顆粒，其沿Y軸方向的相對位移與線速度在剪切流動區與非剪切流動區內差異較為明顯。比較-Y軸-0.03 m和-0.04 m處，+Y軸0.03 m和0.04 m處，Z軸0.03 m和0.04 m處的相對位移與線速度可以發現，-Y軸分布的顆粒其值最大，而+Y軸分布的顆粒其值最小。其原因是Y軸方向(豎直方向)的重力與-Y軸方向上促使顆粒運動的流速梯度同向，加速了顆粒向下運動。

3.2 剪切流動區內粗骨料顆粒運動原因分析

粗骨料顆粒在X軸方向上(管道軸向)不同的位移與線速度差異是由于Bingham流體管道輸送流速分布特性造成的。剪切流動區域內沿管徑方向上相對位移產生的原因是由于流速沿管徑方向上的梯度差造成粗骨料顆粒發生自旋運動而引起的。沿Z軸分布以及沿Y軸負方向分布的粗骨料顆粒其X軸線速度分量的流速梯度分別為dvx/dZ和dvx/dY，受沿管徑方向流速梯度的影響，處于剪切流動區域的粗骨料顆粒在管道輸送時存在由流場給予的扭矩值。其將使粗骨料顆粒存在一定的自旋運動，產生相應的角速度。在剪切流動區域內依據角速度判定的右手定則，由X軸速度分量引起的角速度方向分別沿Y軸負方向以及Z軸負方向。流速梯度引起的粗骨料顆粒自旋示意圖如圖6所示。

圖6 流速梯度引起的顆粒自旋示意圖Fig.6 Diagram of spin of particles caused by flow rate gradient

Z軸方向分布與-Y軸方向分布的粗骨料顆粒，計算時長0.5 s，顆粒沿Y軸與Z軸的旋轉角速度分量如圖7所示。

圖7Z軸分布與Y軸分布時顆粒角速度分量Fig.7 Angular velocity component of particles distributed alongZ-axis andY-axis

由圖7可知：粗骨料顆粒在剪切流動區域內的旋轉角速度明顯大于非剪切流動區域內的旋轉角速度。流速梯度引起的顆粒自旋是造成顆粒在剪切流動區內與偽勻質懸浮液發生相對運動的主要原因。在剪切流動區域內顆粒在Z軸分布時，其對應的Y軸角速度分量為負值，在-Y軸分布時，其對應的Z軸角速度分量為負值，與圖6中推斷的顆粒自旋角速度方向相對應，定性證明了MPM模型計算結果的可信性。

4 結論

1) 將高濃度、寬粒級全尾砂膏體料漿視為偽勻質懸浮液模型，分析其流變學屬性，采用 2參數的Bingham模型研究其管道輸送時的流速分布特性。

2) 針對粗骨料膏體中添加的粗骨料顆粒，依據Bingham流體的管道輸送特性，通過分析粗骨料顆粒在剪切流動區與非剪切流動區內顆粒的主要運動形式，構建了粗骨料膏體管道輸送的兩分區復合流動模型。

3) 采用MPM數值模型，研究粗骨料顆粒與偽勻質懸浮體復合流動時的復雜相互作用。通過分析顆粒在X軸、Y軸、Z軸的相對位移與相對線速度。得出粗骨料顆粒在非剪切流動區域內與偽勻質懸浮體之間基本不發生相對運動，而在剪切流動區域內存在較明顯的徑向偏移與軸向差速運動。

4) 解釋了粗骨料顆粒在剪切流動區域內發生相對運動的主要原因，即沿管徑方向的流速梯度引發粗骨料顆粒自旋運動而造成的。數值計算了顆粒的自旋角速度分量，并與理論分析做了定性的對比，所得結果具有一定的可信性。