基于群體平衡理論的豎直管內水合物漿液流動特性模擬

姚淑鵬，李玉星，王武昌，宋光春，姜 凱，施政灼

（中國石油大學（華東）山東省油氣儲運安全省級重點實驗室，山東 青島 266580）

作為一種清潔能源，天然氣水合物（簡稱水合物）礦藏受到了越來越多國家的關注，水合物礦藏的開發和開采成為人們研究的熱點問題［1］。水合物礦藏主要分布在極地和深水區域，其中，深水區域約占總儲量的95%［2-3］。目前，世界上試開采的多為成巖的水合物礦藏，而對于總量大卻埋深淺、膠結性差的非成巖水合物礦藏，我國周守為院士提出了固態流化綠色開采方法［4-5］。

研究固態流開采過程中管內的多相流動，對了解豎直舉升管內水合物流動特性有重要意義，但實驗難度大、成本高，不宜推廣。近年來，隨著計算流體力學（CFD）的發展，數值模擬成為替代實驗進行相關研究的有效工具。Fatnes［6］采用Ansys CFX軟件對水平管內的水合物流動特性進行了模擬。Balakin等［7-8］通過CFD模擬，對CCl3F（R11）水合物在紊流流態純水體系中的流動沉積現象進行了研究，并通過求解該模型模擬了管內水合物顆粒的粒徑分布。江國業等［9］采用正交實驗，發現速度對壓降的影響最大。宋光春等［10-11］通過構建群體相平衡模型對R11水合物的流動特性進行了模擬研究，發現流速、水合物體積分數及水合物顆粒在管道內分布情況都會對水合物漿液的流動壓降產生影響。周守為等［12］通過有限差分迭代法對豎直井筒內的水合物的分解臨界點等流動特性進行了研究，對深水淺層水合物礦藏的開采提出了指導性意見。劉艷軍等［13］對水合物漿體垂直管輸的帶相變流動特性進行研究，發現水合物分解對漿體管道運輸具有減阻作用。上述研究沒考慮水合物顆粒在豎直管內流動輸送過程中的聚集與破碎，與實際過程存在一定差距，因此值得進一步完善。

本工作引入基于水合物顆粒聚集動力學的群體平衡模型對水合物漿液在豎直管內的流動特性進行了模擬，通過群體平衡模型，模擬了水合物顆粒在管道內的聚集與破碎，較好地模擬了水合物漿液在豎直管內的流動狀態，進而對不同平均流速、體積分數下的水合物漿液流動特性進行研究，為深水淺層水合物礦藏開采的流動安全保障和水合物漿液輸送技術提供理論支持。

1 模型建立

1.1 幾何模型

幾何模型的構建基于Balakin等開展的R11水合物漿液流動特性實驗所使用的管道，對它進行了三維模型的構建，管道長2.0 m，直徑4.52 cm。對構建好的幾何模型進行網格劃分，考慮到流體在邊界區域的物理參數變化相對比較劇烈，因此，對進口壁面進行了邊界層效應的處理，共設6層邊界層。對其他網格均以1 mm劃分。本次建模共有170 180個六面體網格，經過網格獨立性檢驗，網格質量為0.917。模型進口截面示意圖見圖1。

圖1 模型進口截面示意圖Fig.1 Schematic diagram of the inlet section of the model.

1.2 模型基本假設

建模過程中的基本假設有:水合物漿液在豎直管內流動過程等溫；忽略水合物漿液內相間質量傳遞，即不考慮水合物的生成及分解；流體是不可壓縮介質；水合物顆粒均為連續性介質；不考慮水合物在管壁上的黏附。

1.3 多相流模型及湍流模型

多相流模型選用歐拉-歐拉雙流體模型，由連續性方程（見式（1））、動量方程（見式（2））以及用來封閉方程組的本構方程組成。

本次模擬的水合物漿液的流動為固液兩相流，因此，需要重點考慮液固的耦合。在Fluent 模擬過程中，液固耦合由相間動量交換實現。本次模擬在計算相間動量交換時，主要考慮相間拖曳力，它由摩擦阻力和形狀阻力組成。相間拖曳力采用Gidaspow 模型［13］計算。

當αs≤20%時，采用Wen-Yu 公式（見式（3））計算。

當αs＞20%時，采用Ergun 公式（見式（4））計算。

則相間拖曳力可以由式（5）計算。

本文模擬中，還需要對水合物顆粒的黏度進行求解，見式（6）［14］。

R11水合物漿液與天然氣水合物漿液的表面黏度分別由Roscoe-Brinkman方程（見式（7））［15-16］和 Thomas方程（見式（8））［17］計算。

根據式（6）～式（8）可以分別編制水合物顆粒黏度的用戶自定義函數（UDF），用以計算R11水合物和天然氣水合物顆粒相的黏度。

湍流模型則采用Fluent 15.0 軟件自帶的標準k-ε模型。

1.4 群體平衡模型

最初，群體平衡模型用來對人口數量的動態平衡進行描述。后來，Hulburt等［18］將群體平衡模型引入到化工領域中，通過跟蹤氣泡、液滴、顆粒等離散相的數量密度函數的變化，將離散相的成核、生長、聚集、破碎等微觀行為與它的體積、面積、粒徑等宏觀屬性聯系起來，這樣既保留了離散相的微觀屬性，又通過平均化的方式表征了連續相的運動，減少了計算量。

在群體平衡模型引入中，根據本文假設，可使用式（9）的群體平衡方程。

碰撞頻率方面，主要考慮由差速沉降與流動剪切造成的碰撞，并取兩者碰撞頻率之和作為水合物顆粒的實際碰撞頻率。其中，差速沉降碰撞頻率采用Camp等［19］提出的公式（見式（10））計算。

流動剪切碰撞頻率方面，當水合物顆粒小于Kolmogorov尺度時，其處于湍流耗散區，在這一區域內水合物的聚集行為主要受渦內局部剪切力的影響。此時，水合物顆粒的碰撞頻率可用Saffman等［20］提出的公式（見式（12））計算。

當水合物顆粒大于Kolmogorov尺度時，其處于湍流慣性區，在這一區域內水合物被主流場牽引運動。此時，水合物顆粒的碰撞頻率可用Abrahamson［21］提出的公式（見式（13））計算。

采用曲線模型計算水合物顆粒間的聚并效率，由于連續相為水，水合物顆粒間不存在液橋力，故聚并效率計算時主要考慮范德華力與流動剪切力之比。聚并效率按式（14）［22］計算。

對水合物顆粒破碎效應的模擬選用Fluent 15.0軟件自帶的lehr模型。

根據式（10）～（14）編制UDF，并以此計算群體平衡模型中的關鍵參數。

1.5 模型求解

通過Fluent 15.0 軟件，對模型進行求解。進口設定為速度進口，出口為壓力出口，出口壓力為0，壁面無滑移。考慮到水合物開采時，井筒內的實際工況重力方向與水合物漿液的流動方向相反，所以本次模擬中設定重力方向與流動方向相反；本次模擬選取二階迎風格式；當各因子的殘差收斂到10-5時，認為已經收斂，模擬結束。部分模擬參數見表1，模擬工況見表 2［23］。

表1 模型參數表Table 1 Model parameter

表2 模擬工況參數表Table 2 Simulation condition parameter

1.6 模型驗證

目前，由于豎直管內的水合物漿液的流動實驗較少，相關實驗數據尚不完備。為驗證本次模型的準確性，依據Balakin等［7-8］開展的水平管內R11水合物漿流動特性的實驗對模型進行驗證，主要以單位壓降作為驗證參數。模型驗證的工況見表3，單位壓降的驗證結果見表4。由表4可知，模型能夠較準確地模擬水合物漿液的流動狀態。

表3 模型驗證的工況參數Table 3 Condition parameters of model verification

表4 單位壓降的驗證結果Table 4 Verification results of unit pressure drop

2 流動特性研究

固態流開采過程中，采掘設備將水合物粉碎成細小顆粒后，將水合物顆粒與海水混合，然后通過密封管道輸送至海上平臺。水合物漿液在管道輸送的過程中，流經一段較長的豎直舉升管。因此，主要研究了水合物漿液的不同舉升速度，以及與海水不同比例混合的水合物濃度對豎直管內水合物漿液流動壓降、速度、濃度分布以及水合物粒徑分布等流動特性的影響。

2.1 豎直管內水合物漿液顆粒粒徑分布

水合物漿液在豎直管內的流動過程中，會存在水合物顆粒的聚集和破碎，從而導致水合物漿液在流動中產生的不同粒徑分布。圖2為不同入口流速、不同體積分數下出口截面x軸方向水合物顆粒的粒徑分布。

由圖2可知，9種工況下水合物顆粒粒徑在整個出口截面上均呈現對稱分布，同時，近壁面處水合物顆粒粒徑較大且粒徑變化梯度較大，管道中間部分水合物顆粒粒徑較小且粒徑較為均勻。根據管道內水合物顆粒的聚集動力學可知，水合物顆粒發生碰撞聚集主要是流動剪切造成的。水合物漿液在管道近壁面處流動剪切最強，使得水合物顆粒在近壁面處發生碰撞聚集的概率最高，導致水合物顆粒的粒徑增大。而在管道中間部分，水合物漿液的流動剪切作用較弱且水合物顆粒分布較為均勻，所以水合物顆粒發生碰撞聚集的概率較小，粒徑也較小。

圖2 不同工況下的出口截面X軸方向水合物顆粒的粒徑分布Fig.2 Particle size distribution of hydrate particles in the outlet section X direction under different working conditions.

對比相同入口流速，不同水合物漿液濃度的水合物顆粒分布可知，在豎直管內，水合物漿液濃度越高，近壁面處的顆粒粒徑越大，且粒徑分布梯度也越大，不利于輸送的安全性。這是由于水合物漿液濃度越高，水合物顆粒的數目越多，發生碰撞的概率隨之增大，在流動剪切較強的近壁面處，水合物顆粒發生劇烈的碰撞聚集，使水合物顆粒粒徑不斷增大。而在流動剪切作用較弱的管道中間部分，水合物顆粒粒徑碰撞聚集的概率較小，分布較為均勻。對比相同水合物漿液濃度、不同入口流速的漿液粒徑分布可知，在豎直管內，水合物漿液流速越高，水合物顆粒的最大粒徑反而降低，這是因為雖然流速越高，流動剪切作用越強，水合物顆粒碰撞聚集的概率越大，但流動剪切越強，較大的水合物顆粒也越容易破碎，能維持的最大聚集粒徑也就越小，這對輸送的安全性有利。

2.2 豎直管內水合物漿液濃度分布

不同入口流速、不同體積分數下出口截面水合物顆粒濃度分布云圖見圖3。

由圖3可知，9種工況的模擬結果都為均勻懸浮流。這是因為在豎直管內，由于水合物漿液流動方向為逆重力方向，兩者可認為在一條直線上，所以沒有出現類似水平管中因重力影響而使濃度在整個截面上分布不對稱的情況［25］，而是對稱分布的。同時，9種工況的濃度分布都呈現出近壁面處高，管道中間部分低的分布狀況。這是由于近壁面處水合物漿液流速較低，而管道中間部分流速較高所致，流速較高會使水合物顆粒的分散系數增大，從而導致管道中間部分濃度較低且分布較均勻，而四周因為流速降低而造成不均勻分布，出現水合物顆粒的堆積。

當流速相同時，隨著水合物漿液中水合物顆粒濃度的增加，水合物顆粒濃度分布依然保持管道中間部分低，近壁面處高的分布狀態，管道中間部分濃度均勻且濃度值在初始設置濃度附近，管內濃度整體隨之增大。隨著漿液流速的增大，相同濃度下水合物漿液的濃度分布基本相同。這說明在豎直管道中，水合物漿液的流速對水合物漿液的濃度分布影響較小。

2.3 豎直管內水合物漿液流動壓降

水合物漿液在豎直管內的流動過程中，流動壓降是備受關注的重點問題之一，因此，對豎直管內9種 工況組合的單位壓降進行了研究，水合物漿液流動壓降見表5。

由表5可知，當水合物濃度相同時，隨著水合物漿液的流速增加，單位壓降也隨之增大。同時，隨著水合物漿液流速的增大，單位壓降的增大梯度也越來越大。以體積分數為10%的水合物漿液為例，當入口平均流速從1.0 m/s增加到1.5 m/s（增大50%）和2.0 m/s（增大100%）時，單位壓降分別增加了104.9%，237.2%。當水合物流速相同時，隨著水合物漿液濃度的增加，單位壓降也隨之增大。同時，隨著水合物漿液濃度的增加，單位壓降的增大梯度也越來越大。以入口平均流速1.5 m/s為例，當水合物體積分數由10%增長到20%（增大100%）和30%（增大200%）時，單位壓降分別增大了13.5%，34.2%。

此外，通過對9組實驗數據的觀察可知，單位壓降的增大主要是由流速增大導致水合物漿液在豎直管中流動摩阻增大造成的，較高的流速不利于輸送的經濟性。

圖3 不同工況下的出口截面濃度分布云圖Fig.3 Distribution cloud diagram of the outlet section concentrations under different working conditions.

表5 水合物漿液流動壓降Table 5 Flow pressure drop of hydrate slurry

3 結論

1）在豎直管道中，水合物顆粒粒徑在整個截面上是對稱分布的，且為近壁面處顆粒粒徑較大，管道中間部分顆粒粒徑較小。同時，水合物顆粒的最大粒徑隨水合物濃度的增加而增大，隨流速的增加而減小。

2）在豎直管道中，水合物漿液的濃度在整個截面上對稱分布，為均勻懸浮流。水合物濃度為近壁面處高，管道中間部分低。同時，隨著水合物濃度的增加，截面上水合物的濃度整體增加，且管道中間部分濃度均勻且濃度值在初始設置濃度附近，而流速對水合物濃度在橫截面上的分布影響不大。

符 號 說 明

a兩水合物顆粒發生碰撞后的聚并效率

b(L|L') 粒徑為L'的水合物顆粒破碎后產生粒徑為L水合物顆粒的概率

CD曳力系數

dp固相顆粒粒徑，m

拖曳力，kg·m/s

G絕對速度梯度，s-1

國民在畜禽產品消費成本方面普遍不高，這得益于國內廉價的飼料與勞動力，但在推進畜牧業現代化發展中，這種觀念可能成為一種主觀障礙。為有效解決問題，需進一步加快物聯網技術在畜牧業發展中的應用。另外，在將RFID物聯網技術等應用服務于現代畜牧業發展中，不可避免地會出現一些行業類似的問題，如標準制定問題、成本控制問題以及信息安全問題等。但隨著現代技術的進步、政府的有效引導以及市場的不斷完善，此類問題必將迎刃而解。

g重力加速度，m/s2

H哈梅克常數

L，L'水合物顆粒粒徑，m

M相間動量交換項，kg/(m·s)2

n(L，t) 粒徑為L的水合物顆粒在t時刻的數量

密度，m-3

p壓力，Pa

R發生碰撞兩水合物顆粒的調和半徑，m

S(L') 粒徑為L'水合物顆粒的破碎頻率，s-1

t時間，s

u 速度矢量，m/s

u平均速度，m/s

V沉降速率，m/s

速度，m/s

相間的相對速度，m/s

z與流體性質有關的常數

▽ 梯度

α體積分數

β相間曳力常數

β(L-L'，L') 粒徑分別為L-L'和L'的水合物顆粒的碰撞頻率，m3/s

βDS差速沉降碰撞頻率，s-1

ρ密度，kg/m3

τ應力張量，Pa

φ體積濃度

μ黏度，mPa·s

下角標

i，j水相或水合物顆粒相

l 液相

m 混合相

s 固相

TH 湍流耗散區

TG 湍流慣性區

參 考 文 獻

［1］ 宋廣喜，雷懷玉，王柏蒼，等. 國內外天然氣水合物發展現狀與思考［J］. 國際石油經濟，2013，21（11）:69-76.

［2］ Boswell R，Collett T S. Current perspectives on gas hydrate resources［J］. Energy Environ Sci，2011，4（4）:1206-1215.

［3］ 許紅，黃君權，夏斌，等. 最新國際天然氣水合物研究現狀及資源潛力評估（下）［J］. 天然氣工業，2005，25（6）:18-23.

［4］ 周守為，陳偉，李清平. 深水淺層天然氣水合物固態流化綠色開采技術［J］. 中國海上油氣，2014，26（5）:1-7.

［5］ 周守為，陳偉，李清平，等. 深水淺層非成巖天然氣水合物固態流化試采技術研究及進展［J］. 中國海上油氣，2017，29（4）:1-8.

［6］ Fatnes E D. Numerical simulations of the flow and plugging behaviour of hydrate particles［D］. Bergen:University of Bergen，2010.

［7］ 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.

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

［9］ 江國業，王曉婭，孫鵬. 基于正交試驗設計的水合物漿液流動特性數值模擬［J］. 科技導報，2014，32（13）:23-27.

［10］ 宋光春，李玉星，王武昌，等. 基于群體平衡理論的管內水合物漿流動特性數值模擬［J］. 化工進展，2018，37（2）:561-568.

［11］ 宋光春，李玉星，王武昌，等. 基于群體平衡理論的水合物聚集動力學模型［J］. 化工進展，2018，37（1）:80-87.

［12］ 周守為，趙金洲，李清平，等. 全球首次海洋天然氣水合物固態流化試采工程參數優化設計［J］. 天然氣工業，2017，37（9） :1-14.

［13］ 劉艷軍，唐孝蓉，胡坤. 天然氣水合物漿體分解對其在垂直管中流動特性影響的研究［J］. 化學通報，2018，81（3）:267-273.

［14］ Ding J，Gidaspow D. A bubbling fluidization model using kinetic theory of granular flow［J］. AIChE J，1990，36（4）:523-538.

［15］ 王繼紅. 冰漿的管道輸送熱流動特性［D］. 大連:大連理工大學，2013.

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

［17］ 趙鵬飛，王武昌，李玉星，等. 管道內水合物漿流動的數值模型［J］. 油氣儲運，2016，35（3）:272-277.

［18］ Hulburt H M，Katz S. Some problems in particle technology［J］. Chem Eng Sci，1964，19（8）:555-574.

［19］ Camp T R，Stein P C. Velocity gradients and internal work influid motion［J］. J Bsn Soc Civ Eng，1943，30（4）:219-237.

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

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

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

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

［24］ And X L，Logan B E. Collision frequencies between fractal aggregates and small particles in a turbulently sheared fluid［J］.Environ Sci Technol，1997，31（4）:1237-1242.

［25］ 陳鵬，劉福旺，李玉星，等. 水合物漿液流動特性數值模擬［J］. 油氣儲運，2014，33（2）:160-164.