宋光春，李玉星，王武昌，姚淑鵬，魏 丁，房祥鵬

（1.山東省油氣儲運(yùn)安全省級重點(diǎn)實(shí)驗(yàn)室 中國石油大學(xué)（華東），山東 青島 266580；2.泰能天然氣有限公司，山東 青島 266580；3.青島新奧燃?xì)庥邢薰?，山東 青島 266580）

研究管道內(nèi)水合物顆粒的聚集過程對深水油氣管道的流動安全保障具有重要意義［1-2］。目前，對管道內(nèi)水合物顆粒聚集過程的實(shí)驗(yàn)多局限在靜態(tài)條件下［3-4］，而研究流動條件下水合物顆粒的聚集過程則因?yàn)槿狈τ行z測手段而大大受限。因此，隨著近年來計算流體力學(xué)的發(fā)展，數(shù) 值模擬已成為研究水合物顆粒聚集過程的有效手段。

李瑩玉等［5］利用彈性碰撞理 論描述了水合物顆粒間的聚集，進(jìn)而模擬了流動過程中水合物顆粒的粒徑變化，然而該模擬只考慮了對心碰撞且忽略了由碰撞導(dǎo)致的顆粒破碎。引入水合物顆粒聚集粒徑計算模型［6-7］進(jìn)行水合物顆粒聚集模擬是國內(nèi)外學(xué)者普遍采用的一種方法［8-9］，但該方法取水合物顆粒間的聚集力為常數(shù)，且并未考慮水合物顆粒間的聚集效率。Balakin等［10-11］建立了水合物顆粒數(shù)量密度的群體平衡模型，并通過求解該模型模擬了管道內(nèi)水合物顆粒的粒徑變化。該模擬方法由于較為全面地考慮了水合物顆粒在流動過程中的聚集和破碎，因此模擬結(jié)果與工程實(shí)際更為接近。但他們在模型建立過程中進(jìn)行了簡化，假設(shè)了發(fā)生聚集的兩個水合物顆粒粒徑是相同的，因此該模型還有待進(jìn)一步完善。

針對以上數(shù)值模型模擬的不足和缺陷，本工作引入了基于水合物顆粒聚集動力學(xué)的群體平 衡模型，通過碰撞頻率和聚集效率描述水合物顆粒間的聚集，通過破碎頻率和破碎后子顆粒的粒徑分布函數(shù)描述水合物顆粒的破碎。再利用FLUENT 14.5 軟件，對固液兩相流模型和群體平衡模型進(jìn)行聯(lián)合求解，實(shí)現(xiàn)對管道內(nèi)水合物顆粒聚集過程的數(shù)值模擬。

1 數(shù)值模型

1.1 幾何模型

以 Balakin等［12］開展的管內(nèi) CCl3F（R11）水合物漿流動特性實(shí)驗(yàn)為模擬對象，將模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對比。根據(jù)實(shí)驗(yàn)所用環(huán)道建立了三維水平圓管模型（見圖1），管道長度為3 m，管道內(nèi)徑為45.2 mm。對幾何模型進(jìn)行六面體網(wǎng)格劃分，進(jìn)口壁面以1 mm為步長進(jìn)行劃分；在近壁面處，為處理邊界層效應(yīng)，采用8層網(wǎng)格加密，其余網(wǎng)格均以1 mm為步長進(jìn)行劃分。圖1幾何模型總計267 696個六面體網(wǎng)格，網(wǎng)格質(zhì)量0.913。對網(wǎng)格數(shù)量進(jìn)行獨(dú)立性檢驗(yàn)（見圖2）。從圖2可看出，當(dāng)網(wǎng)格數(shù)分別為527 616個和267 696個時，管流（2 m/s）充分發(fā)展段近壁面處速度梯度的相對偏差較?。?.3%），可見網(wǎng)格數(shù)量滿足獨(dú)立性要求。

1.2 多相流模型及湍流模型

在建模過程中采用以下假設(shè)：流動過程等溫且忽略相間質(zhì)量傳遞，即不考慮水合物的生成及分解；水合物漿由水相和水合物顆粒相構(gòu)成，不考慮氣相和油相；水合物顆粒粒徑連續(xù)分布；不考慮水合物顆粒在管壁上的黏附過程。在上述假設(shè)下，采用的物理模型主要包括多相流模型、湍流模型和群體平衡模型。其中，多相流模型采用歐拉-歐拉雙流體模型，由連續(xù)性方程（見式（1））、動量方程（見式（2））控制方程和若干用于封閉方程組的本構(gòu)方程構(gòu)成。

式中，i為水相或水合物顆粒相；t為時間，s；ρ為密度，kg/m3；α為體積分?jǐn)?shù)，%；?為拉普拉斯算子；u為速度矢量，m/s；p為壓力，Pa；τi為應(yīng)力張量，Pa；Mi則為相間動量交換項(xiàng)，kg/（m·s）2。

水合物漿的液固耦合是多相流模型建立時考慮的重點(diǎn)，在模擬過程中，液固耦合通過相間動量交換實(shí)現(xiàn)。在計算Mi時，主要考慮相間 拖曳力（Mdi）和湍流擴(kuò)散力（Mti），并以此建立相間作用力模型以封閉多相流模型，見式（3）。

式中，ur為相間相對流速，m/s；μtm為湍流黏度，kg/（m·s）；σd為Planck擴(kuò)散系數(shù)；αs為水合物體積分?jǐn)?shù)，%；kls為動量傳遞系數(shù)，本工作采用Gidaspow模型［13］計算：當(dāng)αs≤20%時，由Wen-Yu公式（見式（4））計算；當(dāng)αs＞20%時，由Ergun公式（見式（5））計算。

式中，L為水合物顆粒粒徑，m；μl為水相動力黏度，kg/（m·s）；CD為曳力系數(shù)。

除了相間作用力模型，多相流模型的封閉還需確定水合物顆粒相的黏度μs，計算式見式（6）［14］：

式中，αl為水相體積分?jǐn)?shù)，%；μm為混合相黏度，kg/（m·s），可用 Roscoe-Brinkman 方程［15-16］（見式（7））計算得到。

根據(jù)式（6）～（7）即可編制用戶定義函數(shù)（UDF）并以此計算水合物顆粒相的黏度μs。湍流模型則采用標(biāo)準(zhǔn)k - ε模型，近壁面采用標(biāo)準(zhǔn)壁面函數(shù)進(jìn)行處理。

1.3 群體平衡模型

根據(jù)本工作假設(shè)，不考慮水合物顆粒的生成分解及水合物顆粒的管壁黏附且認(rèn)為水合物顆粒粒徑連續(xù)分布，因此可采用式（8）的群體平衡方程［17］：

式中，n（L，t）表示粒徑為L的水合物顆粒在t時刻的數(shù)量密度，m-3；β（L-L′，L′）表示粒徑分別為L-L′和L′兩水合物顆粒的碰撞頻率，m3/s；a為兩水合物顆粒發(fā)生碰撞后的聚并效率；S（L′）表示粒徑為L′的水合物顆粒的破碎頻率，s-1；b（L|L′）則表示粒徑為L′的水合物顆粒破碎后產(chǎn)生粒徑為L水合物顆粒的概率。

基于水合物顆粒聚集動力學(xué)，對式（8）中關(guān)鍵參數(shù)的計算公式進(jìn)行選取。

碰撞頻率方面，主要考慮由差速沉降和流動剪切造成的碰撞并取兩者碰撞頻率之和作為水合物顆粒的實(shí)際碰撞頻率。其中，差速沉降碰撞頻率（βDS）采用式（9）［18］計算：

式中，V為沉降速率，m/s，可由式（10）計算得到：

對于流動剪切碰撞頻率，當(dāng)水合物顆粒小于Kolmogorov尺度時，用式（11）［19］進(jìn)行計算；當(dāng)水合物顆粒大于Kolmogorov尺度時，用式（12）［20］進(jìn)行計算：

式中，G為流場局部的絕對速度梯度，s-1；u為水合物顆粒的平均速度，m/s。

采用曲線模型計算水合物顆粒間的聚并效率，見式（13）［21］：

式中，H為表征范德華力大小的哈梅克常數(shù)，J；R為發(fā)生碰撞兩水合物顆粒的調(diào)和半徑，m。

計算破碎頻率時則主要考慮由流動剪切導(dǎo)致的破碎，計算式見式（14）：

式中，E和m為經(jīng)驗(yàn)常數(shù)。

采用二元分布作為水合物顆粒破碎后的粒徑分布函數(shù)［22-23］，表達(dá)式見式（15）：

根據(jù)式（9）～（15）編制UDF并以此計算群體平衡模型中的關(guān)鍵參數(shù)。

1.4 模型求解

利用FLUENT 14.5軟件，對固液兩相流模型和群體平衡模型聯(lián)合求解。管道入口邊界條件設(shè)定為速度入口，管道出口邊界條件設(shè)定為壓力出口，以二階迎風(fēng)差分格式離散動量方程，壓力-速度耦合使用SIMPLE算法，群體平衡方程使用離散法進(jìn)行求解，瞬態(tài)模擬時間步長設(shè)置為0.1 s，殘差設(shè)置為10-5。模型求解時所涉及的部分參數(shù)見表1。

表1 模型參數(shù)Table 1 Model parameters

2 結(jié)果與討論

2.1 模擬工況

本工作主要模擬不同條件下水合物顆粒的聚集過程，并分析管內(nèi)流 速和水合物體積分?jǐn)?shù)對水合物顆粒聚集過程的影響。模擬時的流速和水合物體積分?jǐn)?shù)取值均參考文獻(xiàn)［12］，具體見表2。

表2 模擬工況Table 2 Simulated working conditions

在使用群體平衡模型時，需給出水合物顆粒的初始粒徑分布。為保證可對比性，模擬的7組工況 的水合物顆粒初始粒徑均相同，即水合物初始粒徑均為5 μm左右［25-26］（不可破碎），初始粒徑分布見表3。

表3 水合物顆粒 初始粒徑分布Table 3 Initial particle size distribution of hydrate particles

2.2 實(shí)驗(yàn)驗(yàn)證

依據(jù)相關(guān)實(shí)驗(yàn)［12］，取流動壓降和水合物顆粒粒徑分布對數(shù)值模型的可行性和有效性進(jìn)行驗(yàn)證。表4為單位管長壓降實(shí)驗(yàn)值與模型模擬值的對比，由表4可知，工況5～7單位管長壓降實(shí)驗(yàn)值與模擬值的變化趨勢相同且兩者間的相對誤差小于18%。因此，建立的數(shù)值模型能較好地模擬水合物顆粒在管道內(nèi)流動時的壓降變化規(guī)律。

工況1～5水合物顆粒粒徑分布實(shí)驗(yàn)值與模擬值的對比見圖3。由圖3可知，由數(shù)值模型模擬得到的顆粒粒徑分布與通過實(shí)驗(yàn)獲得的顆粒粒徑分布規(guī)律相同，均近似呈對數(shù)正態(tài)分布。因此，建立的數(shù)值模型也可較好地模擬水合物顆粒在管內(nèi)流動時的粒徑分布情況。

表4 實(shí)驗(yàn)和模擬條件下單位管長壓降對比Table 4 Comparison of pressure gradient between experimental results and simulation results

圖3 水合物顆粒粒徑分布實(shí)驗(yàn)值與模擬值對比Fig.3 Comparison of experimental value and simulation value of hydrate particle size distribution.

2.3 低流速下水合物顆粒聚集過程分析

采用工況4～7分析低流速下水合物顆粒的聚集過程。圖4為較低流速（1.5 m/s）下，使用數(shù)值模型模擬和以濃度差為推動力的水合物顆粒生長模型分別計算得到的管道內(nèi)水合物顆粒平均粒徑的變化情況。由圖4可知，在較低流速下，由數(shù)值模型模擬得到的水合物顆粒平均粒徑變化過程可大致分以下4個階段：1）緩慢增長階段。該階段水合物顆粒粒徑較小，碰撞頻率較低，因此平均粒徑增長緩慢。在該階段，隨管道內(nèi)水合物體積分?jǐn)?shù)的升高，水合物顆粒平均粒徑的增長速度逐漸增大。這是因?yàn)楦唧w積分?jǐn)?shù)下水合物顆粒間更容易發(fā)生碰撞和聚集。2）快速增長階段。該階段水合物顆粒粒徑由小變大，碰撞頻率升高，水合物顆粒平均粒徑開始快速增長。在該階段，管道內(nèi)水合物的體積分?jǐn)?shù)越高，水合物顆粒平均粒徑開始快速增長的時間就越早，同時水合物顆粒平均粒徑的增長速度和所能達(dá)到的最大粒徑也越大。3）破碎主導(dǎo)階段。在階段1和階段2，顆粒聚集占據(jù)了主導(dǎo)地位，因此水合物顆粒平均粒徑呈增長趨勢。當(dāng)水合物顆粒粒徑增大到一定程度時，破碎過程將超過聚集過程占據(jù)主導(dǎo)地位，此后水合物顆粒進(jìn)入破碎主導(dǎo)階段，平均粒徑開始逐漸降低。值得注意的是，工況4由于水合物濃度較低、顆粒粒徑較小，故在平均粒徑變化過程中水合物顆粒的聚集過程和破碎過程能逐漸達(dá)到平衡，因此并未出現(xiàn)明顯的破碎主導(dǎo)階段。4）動態(tài)平衡階段。當(dāng)水合物顆粒粒徑降至一定程度時，聚集過程重新占據(jù)主導(dǎo)地位，水合物顆粒平均粒徑開始逐漸增大，直到管道內(nèi)水合物顆粒的破碎過程和聚集過程達(dá)到動態(tài)平衡。此后，管道內(nèi)水合物顆粒的平均粒徑逐漸維持穩(wěn)定，不再發(fā)生明顯變化。在該階段，管道內(nèi)水合物體積分?jǐn)?shù)越高，水合物顆粒所能維持的最終穩(wěn)定粒徑就越大。

圖4 低流速下由數(shù)值模型模擬和生長模型分別計算得到的管道內(nèi)水合物顆粒平均粒徑變化情況Fig.4 Variations of average particle sizes of hydrate particles in pipeline calculated by numerical model simulation and growth model at low flow rates.

圖4中使用的以濃度差為推動力的水合物顆粒生長模型見式（16）［11］：

式中，r為水合物顆粒粒徑的生長速率，m/s；kS為水合物顆粒的生長速率常數(shù)，m4/(s·mol)；cb為液相主體中R11的濃度，mol/m3；ceq為水合物顆粒與液相主體交界面處R11的濃度，mol/m3，當(dāng)水合物顆粒生長過程結(jié)束時ceq=cb。假設(shè)R11可全部溶解在水中且水合物生長過程中R11呈指數(shù)規(guī)律消

耗，則cb和ceq可分別通過式（17～18）進(jìn)行計算：

式中，cin為液相主體中R11的初始濃度，mol/m3；tmax為水合物顆粒生長過程結(jié)束所需的時間，s。

在使用上述生長模型時，本工作假設(shè)水合物顆粒的生長驅(qū)動力始終處于最大狀態(tài)，即始終有cb=cin。模型計算過程中，kS取1.3×10-9m4/(s·mol)［11］，cin則可根據(jù)文獻(xiàn)［12］報道的實(shí)驗(yàn)數(shù)據(jù)計算得到。從圖4可看出，即使在保持最大生長驅(qū)動力的情況下，各工況內(nèi)由水合物顆粒聚集導(dǎo)致的粒徑增長仍較由水合物顆粒生長導(dǎo)致的粒徑增長更為快速、劇烈。由水合物顆粒生長帶來的粒徑增長速率僅相當(dāng)于聚集過程中緩慢增長階段的粒徑增長速率。這說明，相對于顆粒生長，水合物顆粒間的聚集更容易導(dǎo)致顆粒粒徑增大、是流動過程中水合物顆粒粒徑增大的主要原因。

2.4 高流速下水合物顆粒聚集過程分析

利用工況1～3分析高流速下水合物顆粒的聚集過程。較高流速（2.4～3.9 m/s）下，通過數(shù)值模型模擬得到的管道內(nèi)水合物顆粒平均粒徑的變化情況見圖5。由圖5可知，不同于低流速（1.5 m/s）下的工況，高流速下管道內(nèi)水合物顆粒的平均粒徑變化過程只包含快速增長（階段Ⅰ）和動態(tài)平衡（階段Ⅱ）階段。

圖5 高流速下由數(shù)值模型模擬計算得到的管道內(nèi)水合物顆粒平均粒徑變化情況Fig.5 Variations of average particle size of hydrate particles in pipeline calculated by numerical model simulation at high flow rates.

根據(jù)水合物顆粒聚集動力學(xué)可知，流動剪切是造成水合物顆粒發(fā)生碰撞聚集的主要原因。管道內(nèi)流體流速越高，流動剪切就越強(qiáng)，因此高流速下水合物顆粒往往能更快、更劇烈地發(fā)生碰撞聚集。故高流速工況下水合物顆粒平均粒徑增長較快，粒徑變化過程中并未出現(xiàn)明顯的緩慢增長階段。另一方面，同工況4類似，工況1～3與工況5～7相比，水合物濃度較低、顆粒粒徑較小，故在粒徑變化過程中水合物顆粒的聚集過程和破碎過程能逐漸達(dá)到平衡，因此并未出現(xiàn)明顯的破碎主導(dǎo)階段。此外，從圖5還可看出，管道內(nèi)流體流速越高，快速增長階段水合物顆粒平均粒徑的增長速度和其所能達(dá)到的最大粒徑值就越大，這同樣是由高流速下水合物更容易發(fā)生碰撞聚集所致的。

3 結(jié)論

1）建立的數(shù)值模型能較好地模擬水合物顆粒在管道內(nèi)流動時的壓降變化規(guī)律以及粒徑分布情況。

2）低流速下管道內(nèi)水合物顆粒的聚集過程可大致分為緩慢增長、快速增長、破碎主導(dǎo)和動態(tài)平衡四個階段。高流速下管道內(nèi)水合物顆粒的聚集過程則只包括快速增長和動態(tài)平衡兩個階段。

3）在相同流速下，管道內(nèi)水合物體積分?jǐn)?shù)越高，水合物顆粒的平均粒徑就越大。在相同水合物體積分?jǐn)?shù)下，管道內(nèi)水合物顆粒的平均粒徑隨流速的升高而增大。

4）在管內(nèi)流動條件下，相對于顆粒生長，水合物顆粒間的聚集更容易導(dǎo)致顆粒粒徑增大、是流動過程中水合物顆粒粒徑增大的主要原因。

［1］ Sloan E D. Natural Gas Hydrates in Flow Assurance［M］. New York：Elsevier Science Ltd，2010：1-36.

［2］ Song Guangchun，Li Yuxing，Wang Wuchang，et al. Investigation of hydrate plugging in natural gas+diesel oil+water systems using a high-pressure fl ow loop［J］.Chem Eng Sci，2017，158：480-489.

［3］ Dieker L E，Aman Z M，George N C，et al. Micromechanical adhesion force measurements between hydrate particles in hydrocarbon oils and their modifications［J］.Energy Fuels，2009，23（12）：5966-5971.

［4］ Morrissy S A，Lim V W，May E F，et al. Micromechanical cohesive force measurements between precipitated asphaltene solids and cyclopentane hydrates［J］.Energy Fuels，2015，29（10）：6277-6285.

［5］ 李瑩玉，江國業(yè)，陳世一. 基于CFD的集輸管道內(nèi)水合物聚集行為仿真分析［J］.計算機(jī)應(yīng)用與軟件，2013，30（7）：101-103.

［6］ Muhle K. Flock Stability in Laminar and Turbulent Flow. In：Coagulation and Flocculation：Theory and Applications［M］.New York：Marcel Dekker，1996.

［7］ Camargo R，Palermo T. Rheological properties of hydrate suspensions in an asphaltenic crude oil［C］// International Conference on Gas Hydrates，Yokohama：ICGH4，2002.

［8］ 陳鵬，劉福旺，李玉星，等. 水合物漿液流動特性數(shù)值模擬［J］.油氣儲運(yùn)，2014，33（2）：160-164.

［9］ Fatnes E D. Numerical simulations of the flow and plugging behavior of hydrate particles［D］.Bergen：University of Bergen，2010.

［10］ Balakin B V，Hoffmann A C，Kosinski P. Experimental study and computational fluid dynamics modeling of deposition of hydrate particles in a pipeline with turbulent water flow［J］.Chem Eng Sci，2011，66（4）：755-765.

［11］ Balakin B V，Hoffmann A C，Kosinski P. Population balance model for nucleation，growth，aggregation，and breakage of hydrate particles in turbulent fl ow［J］.AIChE J，2010，56（8）：2052-2062.

［12］ Balakin B V，Pedersen H，Kilinc Z，et al. Turbulent fl ow of freon R11 hydrate slurry［J］.J Pet Sci Eng，2010，70（3/4）：177-182.

［13］ Ding Jianmin，Gidaspow D. A bubbling fluidization model using kinetic theory of granular fl ow［J］.AIChE J，1990，36（4）：523-538.

［14］ 王繼紅. 冰漿的管道輸送熱流動特性［D］.大連：大連理工大學(xué)，2013.

［15］ Pabst W. Fundamental considerations on suspension rheology［J］.P R Soc A，2004，48（1）：6-13.

［16］ 趙鵬飛，王武昌，李玉星，等. 管道內(nèi)水合物漿流動的數(shù)值模型［J］.油氣儲運(yùn)，2016，35（3）：272-277.

［17］ Hulburt H M，Katz S. Some problems in particle technology：A statistical mechanical formulation［J］.Chem Eng Sci，1964，19（8）：555-574.

［18］ Camp T R，Stein P C. Velocity gradients and internal work in fl uid motion［J］.J Bsn Soc Civ Eng，1943，30（4）：219-237.

［19］ Saffman P G，Turner J S. On the collision of drops in turbulent clouds［J］.J Fluid Mech，1956，1：16-30.

［20］ Abrahamson J. Collision rates of small particles in a vigorously turbulent fluid［J］.Chem Eng Sci，1975，30（11）：1371-1379.

［21］ Ven T G M V D，Mason S G. The microrheology of colloidal dispersions Ⅶ . Orthokinetic doublet formation of spheres［J］.Colloid Polym Sci，1977，255（5）：468-479.

［22］ Zhang Jianjun，Li Xiaoyin. Modeling particle-size distribution dynamics in a flocculation system［J］.AIChE J，2003，49（7）：1870-1882.

［23］ 李振亮. 基于群體平衡的活性污泥絮凝動力學(xué)［D］.重慶：重慶大學(xué)，2014.

［24］ Li Xiaoyin，Logan B E. Collision frequencies between fractal aggregates and small particles in a turbulently sheared fl uid［J］.Environ Sci Technol，1997，31（4）：1237-1242.

［25］ 呂曉方，吳海浩，史博會，等. 流動體系中二氧化碳水合物堵管時間實(shí)驗(yàn)研究［J］.實(shí)驗(yàn)室研究與探索，2013，32（11）：197-202.

［26］ Turner D J，Miller K T，Sloan E D. Direct conversion of water droplets to methane hydrate in crude oil［J］.Chem Eng Sci，2009，64（23）：5066-5072.