基于Fluent的管道內(nèi)水合物漿液流動特性數(shù)值模擬

鐘一華，卿亞麗，楊海川，李璐伶，劉 婷，蔣 鵬

（1.天津大學(xué) 管理與經(jīng)濟學(xué)部，天津300072；2.深圳市燃氣集團股份有限公司，廣東 深圳518049；3.廣東大鵬LNG公司，廣東 深圳518040）

天然氣水合物是由甲烷、乙烷和二氧化碳等氣體分子與水形成的復(fù)雜籠型晶體[1,2]。在油氣輸送系統(tǒng)中，特別是在高壓低溫條件下，水合物顆粒極易生成，在流動過程中會不斷聚集，最終部分堵塞或完全堵塞管道，造成嚴重的設(shè)備與安全問題[3]。故有必要研究水合物漿液在管道中的流動特性和顆粒行為。

對管道中水合物漿液流動特性的研究主要基于實驗和理論計算。Carmargo等[4]在高壓（9 MPa）和低溫（258~338 K）環(huán)路中，對水合物在含瀝青質(zhì)的原油中的形成過程和流變特性進行了實驗。之后，Sinquin等[5]利用同樣的高壓環(huán)路，對天然氣水合物在油基體系中的流變特性進行了測試，并建立了考慮顆粒聚集和微絮凝破壞機理的水合物漿液流變行為數(shù)學(xué)模型，該模型已運用在CSMHyK軟件中。Gainville等[6]使用IFP-Lyre高壓環(huán)路研究了水合物漿液在層流和湍流中的壓降。Joshi等[7]基于?？松梨诘母邏涵h(huán)路實驗裝置，研究了水基天然氣水合物形成機制。雖然高壓環(huán)路實驗提供了天然氣水合物在管道中形成、成核和團聚等信息，但復(fù)雜的實驗設(shè)備和嚴格的高壓條件限制了對水合物漿液在管道流動堵塞方面的研究。利用制冷介質(zhì)（如HCFC、THF、TBAB和R11）替代碳氫水合物的漿液流動低壓實驗，逐漸成為研究水合物漿液流動特性和顆粒行為的有效方法[8]。但在實驗測量方面也有一定的局限性，如在測量顆粒的體積分數(shù)分布和速度矢量分布等方面存在技術(shù)困難。

隨著固液兩相理論的發(fā)展，利用計算流體力學(xué)（CFD）來研究水合物漿液在管道中的流動特性逐漸成為研究熱點。Jassim等[9]運用Fluent對天然氣輸送中的水合物沉積過程進行了模擬，但該研究僅考慮了天然氣對水合物的單向影響。Balakin等[10,11]利用STAR-CD對R11（一氟三氯甲烷）水合物漿液在管道中的流動情況進行了計算。Sule等[12]運用ANSYS CFX對水合物的形成過程進行了考察，并研究了流動速率、管道直徑和水相比例對管道溫度分布的影響。此外，王武昌等[13,14]對HCFC-141b和THF水合物漿液的流動行為進行了研究，并提出凝聚概率來判斷管道的流動安全性。韋雪蕾等[15]基于歐拉雙流體模型，模擬了不同粒徑的水合物漿液在水平管道中的流動情況。

目前，考慮了水合物顆粒流動過程中的粒徑變化和粘度特性的漿液流動數(shù)值模型較少，以及考慮了水合物顆粒密度和顆粒粒徑變化對管道壓降和顆粒體積分數(shù)分布影響的研究較少。本文以歐拉雙流體模型和RSM湍流模型為基準(zhǔn)，通過編制UDF對水合物顆粒粒徑和水合物漿液粘度進行自定義，基于ANSYS Fluent 14.0建立了R11水合物漿液在管道中的流動模型，并對水合物顆粒在較低流速下的沉積特性進行了分析。

1 數(shù)學(xué)模型與模擬

1.1 基本控制方程

采用歐拉雙流體模型來描述水合物漿液在管道內(nèi)的固液流動。假設(shè)流動過程只含有水合物顆粒及水兩相，水合物為光滑球形顆粒，忽略相間質(zhì)量和熱量傳遞，則液相和顆粒相的連續(xù)性方程可分別表示為：

式中，α為體積分數(shù)，%；ρ為密度，kg/m3； 為拉普拉斯算子； 為速度矢量，m/s；下標(biāo)s和l分別表示水合物顆粒相和液相。

液相和水合物相的動量守恒方程可用Navier-Stokes方程表示為：

式中，p為壓力，Pa；g為重力加速度，m/s2； 和 分別為顆粒相和液相的應(yīng)力張量，Pa；ps為顆粒壓力，Pa；M為相間動量交換系數(shù)，kg/(m·s)2；μl、μs和ζs分別為液體、固體和流體的剪切粘度，kg/(m·s)；I為單位向量。

顆粒壓力可表示為：

式中，Θs為顆粒溫度，K；ess為顆粒恢復(fù)系數(shù)，ess=0.9；g0,ss為顆粒徑向分布函數(shù)，其表達如式(9)所示；αmax為顆粒最大裝填系數(shù)。

假定混合相粘度（μm）與固相顆粒粘度（μs）和液相粘度（μl）成線性關(guān)系，其表達式為式(10)。對于R11水合物漿液，混合相粘度與固相分數(shù)有關(guān)，可用Brinkman方程[10]來擬合，其表達式為式(11)。結(jié)合式(10)和式(11)，則固體顆粒粘度（μs）可表示為式(12)。

相間動量傳遞量M主要考慮相間曳力作用，Mls可表示為：

式中，Kls為Gidaspow曳力模型的交換系數(shù)[16]；CD為曳力系數(shù)；Res為相對雷諾數(shù)；ds為固體顆粒粒徑，μm。

在流動中，水合物顆粒會經(jīng)歷形成、成長、聚并和聚集階段，導(dǎo)致粒徑不斷變化。水合物顆粒粒徑（ds）可用Camargo模型來進行計算[17]，其表達式為：

式中，d1為顆粒初始直徑，μm；γ為剪切速率，s-1；Fa為顆粒間粘附力，nN；fr為分形維數(shù)，由實驗值決定；Balakin等[11]利用群體平衡模型，結(jié)合實驗數(shù)據(jù)，經(jīng)回歸得到d1=7 μm，fr=1.83，F(xiàn)a=1.75 nN。

1.2 數(shù)值模擬

采用的彎曲管道三維模型如圖1所示，該模型依據(jù)Balakin等[11]的R11水合物漿液管內(nèi)流動實驗數(shù)據(jù)建模。其在實驗中測定了R11水合物在常壓環(huán)流管路中的阻力特性及流動形態(tài)，但對于水合物顆粒在管道中的分布仍未明晰。

圖1 管道三維模型

圖1 所示管道的內(nèi)徑為45.2 mm，大寫字母A-F表示管道截面編號。此外，利用ICEM對管道進行六面圖網(wǎng)格劃分，水平、彎曲及進出口端面的網(wǎng)格劃分情況也如圖1所示，對近壁處的網(wǎng)格進行加密處理，以提高計算精度。R11水合物漿液的基本性質(zhì)如表1所示。

表1 R11水合物漿液的兩相基本性質(zhì)

采用雷諾平均Navier-Stokes方程（RANS）來計算流體的湍流特性。本文應(yīng)用雷諾應(yīng)力模型（RSM）來求解水合物漿液在管道流動中的湍流現(xiàn)象。與標(biāo)準(zhǔn)k-ε模型相比，RSM不采用Boussinesq假設(shè)，而是直接計算雷諾應(yīng)力項，故求解時間更長。RSM模型表達式為：

利用Fluent 14.0對上述模型進行求解。入口邊界條件為水合物漿液的速率、湍流強度、水力直徑和顆粒相分數(shù)，此外設(shè)定參考運行壓力為101325 Pa，參考點為（0.35 m,0.814 m,0）。出口邊界條件設(shè)置為壓力出口，其表壓為0。在壁面，顆粒相和液相兩相無滑移現(xiàn)象。壓力-速度耦合使用PC-SIMPLE算法，并以一階迎風(fēng)算法對動量方程、體積分數(shù)、湍流動量和耗散方程進行離散求解，殘差設(shè)置為10-5。

1.3 網(wǎng)格獨立性驗證

對于水合物漿液在彎曲管道中的湍流流動模擬，為減少數(shù)值模擬網(wǎng)格誤差，以水為流體介質(zhì)，計算了網(wǎng)格總數(shù)從156128到1466172個范圍內(nèi)，流體壓降梯度的變化情況。模擬采用的流速為1.50 m/s，結(jié)果如圖2所示。由圖2可知，隨著網(wǎng)格數(shù)的增加，壓降梯度先急劇下降，當(dāng)網(wǎng)格數(shù)大于220000后，壓力梯度開始趨于平穩(wěn)，再增加網(wǎng)格數(shù)量，結(jié)果基本不變。因此，本文選取224721個網(wǎng)格作為后續(xù)計算網(wǎng)格數(shù)。

圖2 網(wǎng)格獨立性驗證

1.4 模型驗證

為驗證本文所建數(shù)值模型的可靠性，利用Fluent 14.0，在水合物體積分數(shù)分別為0、5%、14%和38%的條件下，對流速為0.80~3.90 m/s的管道阻力特性進行了模擬，并將水合物漿液流動壓降梯度的預(yù)測值與Balakin等[18]測得的實驗壓降梯度進行了對比，如圖3所示。由圖3可知，在相同的體積分數(shù)下，壓降梯度的模擬值與實驗值在不同流速下均吻合較好，最大誤差在10%范圍內(nèi)，顯示出所建模型的正確性，故可用該模型來研究水合物漿液在管道流動的壓降及濃度場分布特性。另一方面，隨著流速增加，壓力梯度也呈上升趨勢，該現(xiàn)象可由Darcy-Weisbach[19]方程來描述，即管道內(nèi)壓降梯度與流體速率存在正相關(guān)關(guān)系。當(dāng)水合物體積分數(shù)在0~14%時，在相同流速下的壓降基本相同，但當(dāng)水合物體積分數(shù)達到38%，壓降明顯增大，此時水合物粘度系數(shù)增大，表現(xiàn)出非牛頓流體特性。

圖3 不同流速和水合物顆粒體積分數(shù)下壓降梯度實驗值和模擬值對比

2 結(jié)果及討論

2.1 水合物顆粒密度的影響

在水合物顆粒體積分數(shù)為38%、流速為1.50 m/s的條件下，采用歐拉雙流體模型和RSM湍流模型，考察了水合物顆粒密度分別為800 kg/m3、1138 kg/m3和1300 kg/m3時對水合物漿液壓降梯度和顆粒平均粒徑的影響，如圖4所示。

隨著水合物顆粒密度從800kg/m3增加到1300kg/m3，壓降梯度從548 Pa/m上升到1160 Pa/m，這主要是由于流固交換系數(shù)導(dǎo)致阻力增大和水合物顆粒沉積所致。同時可知，增大顆粒密度，水合物平均直徑從約98.5 μm略微上升至101 μm，其原因為顆粒沉積引起固體顆粒團聚，使得顆粒粒徑增加。

水合物顆粒的體積分數(shù)分布反映了顆粒的聚集特性。圖5顯示了在水合物顆粒密度分別為800 kg/m3和1300 kg/m3情況下，管道不同位置截面處（截面位置見圖1中A~F）水合物顆粒體積分數(shù)分布圖。此外，水合物顆粒垂直于流動方向的速度矢量圖也顯示在了圖中。模擬初始條件為水合物顆粒體積分數(shù)和流速分別為38%和1.50 m/s。當(dāng)顆粒密度為800 kg/m3時，截面A~F靠近中心處的水合物濃度高于壁面處的濃度，這主要是由于二次流的影響，將水合物顆粒向中心集聚。當(dāng)水合物密度為1300 kg/m3，由于離心作用，水合物顆粒向近壁面累積，使截面A的底部和截面B的左側(cè)水合物顆粒體積分數(shù)較高，其最大值為38.6%。另一方面，對于截面C到截面F，顆粒分散相對均勻，水合物顆粒體積分數(shù)分布差異較小。

圖5 水合物顆粒密度對其體積分數(shù)分布的影響

2.2 水合物顆粒粒徑的影響

水合物顆粒的粒徑也會影響其流動阻力特性和體積分數(shù)分布。根據(jù)Camargo顆粒粒徑計算模型（式(18)），水合物顆粒在管道中流動受剪切速率等影響，會發(fā)生聚并，引起粒徑變化。由圖4(b)所示，在水合物顆?；鶞?zhǔn)密度下，其平均粒徑為100 μm。此外，還考察了顆粒平均粒徑在70~120 μm變化時，水合物顆粒的壓降梯度及體積分數(shù)分布的變化，結(jié)果如圖6所示。模擬初始條件為：水合物漿液的入口流速為1.50 m/s，水合物顆粒體積分數(shù)為38%。當(dāng)粒徑為70 μm時，壓降梯度值約為1010 Pa/m，當(dāng)粒徑為120 μm時，壓降梯度降低至860 Pa/m。隨著水合物顆粒平均直徑的增加，水合物顆粒壓降梯度減小。這主要是因為大粒徑的水合物向管壁沉積，使得管道內(nèi)有效流速增加，雷諾數(shù)增加，故壓降出現(xiàn)降低趨勢。

圖6 水合物顆粒平均粒徑對壓降梯度的影響

圖7 顯示了水合物顆粒平均粒徑對不同截面處的水合物顆粒體積分數(shù)分布的影響。隨著顆粒平均直徑的增大，同一截面上的顆粒體積分數(shù)分布存在顯著差異。在平均粒徑為70 μm時，顆粒在不同截面處的分布較為均勻，但當(dāng)平均粒徑增大到120 μm時，在所示的六個截面處，近壁面的水合物顆粒體積分數(shù)高于管道中心處，這主要是大顆粒水合物的聚集，并在離心力的作用下，又在管壁處聚集，引起了管道流通面積的減少，甚至可能導(dǎo)致管道堵塞。

圖7 水合物顆粒平均粒徑對其體積分數(shù)分布的影響

2.3 水合物的沉積特性

油氣管道中水合物顆粒的聚集與沉積是對安全運行的巨大挑戰(zhàn)。水合物沉積會增大壓降梯度，抑制流體流動，導(dǎo)致管道堵塞。在較高的流速下，流體流動呈現(xiàn)均勻狀態(tài)，但當(dāng)平均流動速度小于臨界值時，水合物顆粒便開始沉積在管道中。為量化沉積厚度，以水合物體積分數(shù)為14%計算，設(shè)定水合物漿液流速為0.10～0.40 m/s，水合物顆粒的密度和粒徑如表1所示，對水合物在管道中的體積分數(shù)分布進行了計算。

圖8 顯示了在不同流速（0.10 m/s、0.25 m/s、0.40 m/s）下截面F的水合物體積分數(shù)分布云圖。在截面的上部區(qū)域，水合物體積分數(shù)非常低，接近于0。然而，在截面的中間區(qū)域，水合物濃度開始變得均勻分散。水合物體積分數(shù)呈向截面底部逐漸增加的趨勢。當(dāng)流速分別為0.10 m/s和0.25 m/s時，最大水合物體積分數(shù)為0.55，水合物顆粒已在管道底部形成了固定床層，當(dāng)流速增大到0.40 m/s時，固定床層消失。如圖9所示，流速從0.10 m/s增加到0.40 m/s時，沉積厚度由14 mm降低至0 mm。這主要是由于以下原因?qū)е碌模篟11水合物顆粒的密度大于水，水合物顆粒沉積到以重力為主的底部區(qū)域；在橫截面的中部區(qū)域，顆粒浮力和阻力的組合與重力競爭，使顆粒分布均勻；由于流速的增加，湍流強度及剪切力增強，水合物顆粒的分散能力增強，沉積厚度降低甚至消失。

圖8 不同流速下的水合物顆粒體積分數(shù)分布（F截面）

圖9 水合物顆粒沉積厚度與流速的關(guān)系

3 結(jié)論

基于歐拉雙流體模型和RSM湍流模型，并引入了Brinkman水合物粘度模型和Camargo粒徑模型，利用Fluent 14.0對管道內(nèi)水合物漿液流動和沉積特性進行了三維數(shù)值模擬。網(wǎng)格獨立性分析發(fā)現(xiàn)，當(dāng)網(wǎng)格數(shù)大于220000后，壓降梯度開始趨于平穩(wěn)。對比實驗測定的管路壓降梯度，發(fā)現(xiàn)在不同水合物體積分數(shù)和流速下的壓降梯度的模擬值與實驗值誤差均小于10%，驗證了所建模型的準(zhǔn)確性。

水合物漿液的壓降梯度和顆粒平均尺寸均隨著水合物密度的增加而增大，但水合物顆粒粒徑的增加使得壓降梯度降低。此外，增大水合物顆粒密度和平均粒徑均促進了水合物向管壁聚集，出現(xiàn)了較大的體積分數(shù)分布梯度。水合物顆粒沉積特性表明，對于體積分數(shù)為14%的水合物漿液，當(dāng)流體速率小于0.40 m/s時，水合物顆粒會沉積在管壁上，引起管道堵塞。當(dāng)流速為0.10 m/s時，床層厚度可達14 mm，隨著流速的增加，水合物顆粒沉積床層高度減小。