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

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

（中國石油大學(xué)（華東）山東省油氣儲運(yùn)安全省級重點(diǎn)實(shí)驗(yàn)室，山東 青島 266580）

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

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

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

1 模型建立

1.1 幾何模型

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

圖1 模型進(jìn)口截面示意圖Fig.1 Schematic diagram of the inlet section of the model.

1.2 模型基本假設(shè)

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

1.3 多相流模型及湍流模型

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

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

當(dāng)αs≤20%時，采用Wen-Yu 公式（見式（3））計(jì)算。

當(dāng)αs＞20%時，采用Ergun 公式（見式（4））計(jì)算。

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

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

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

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

湍流模型則采用Fluent 15.0 軟件自帶的標(biāo)準(zhǔn)k-ε模型。

1.4 群體平衡模型

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

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

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

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

當(dāng)水合物顆粒大于Kolmogorov尺度時，其處于湍流慣性區(qū)，在這一區(qū)域內(nèi)水合物被主流場牽引運(yùn)動。此時，水合物顆粒的碰撞頻率可用Abrahamson［21］提出的公式（見式（13））計(jì)算。

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

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

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

1.5 模型求解

通過Fluent 15.0 軟件，對模型進(jìn)行求解。進(jìn)口設(shè)定為速度進(jìn)口，出口為壓力出口，出口壓力為0，壁面無滑移?？紤]到水合物開采時，井筒內(nèi)的實(shí)際工況重力方向與水合物漿液的流動方向相反，所以本次模擬中設(shè)定重力方向與流動方向相反；本次模擬選取二階迎風(fēng)格式；當(dāng)各因子的殘差收斂到10-5時，認(rèn)為已經(jīng)收斂，模擬結(jié)束。部分模擬參數(shù)見表1，模擬工況見表 2［23］。

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

表2 模擬工況參數(shù)表Table 2 Simulation condition parameter

1.6 模型驗(yàn)證

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

表3 模型驗(yàn)證的工況參數(shù)Table 3 Condition parameters of model verification

表4 單位壓降的驗(yàn)證結(jié)果Table 4 Verification results of unit pressure drop

2 流動特性研究

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

2.1 豎直管內(nèi)水合物漿液顆粒粒徑分布

水合物漿液在豎直管內(nèi)的流動過程中，會存在水合物顆粒的聚集和破碎，從而導(dǎo)致水合物漿液在流動中產(chǎn)生的不同粒徑分布。圖2為不同入口流速、不同體積分?jǐn)?shù)下出口截面x軸方向水合物顆粒的粒徑分布。

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

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

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

2.2 豎直管內(nèi)水合物漿液濃度分布

不同入口流速、不同體積分?jǐn)?shù)下出口截面水合物顆粒濃度分布云圖見圖3。

由圖3可知，9種工況的模擬結(jié)果都為均勻懸浮流。這是因?yàn)樵谪Q直管內(nèi)，由于水合物漿液流動方向?yàn)槟嬷亓Ψ较颍瑑烧呖烧J(rèn)為在一條直線上，所以沒有出現(xiàn)類似水平管中因重力影響而使?jié)舛仍谡麄€截面上分布不對稱的情況［25］，而是對稱分布的。同時，9種工況的濃度分布都呈現(xiàn)出近壁面處高，管道中間部分低的分布狀況。這是由于近壁面處水合物漿液流速較低，而管道中間部分流速較高所致，流速較高會使水合物顆粒的分散系數(shù)增大，從而導(dǎo)致管道中間部分濃度較低且分布較均勻，而四周因?yàn)榱魉俳档投斐刹痪鶆蚍植?，出現(xiàn)水合物顆粒的堆積。

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

2.3 豎直管內(nèi)水合物漿液流動壓降

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

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

此外，通過對9組實(shí)驗(yàn)數(shù)據(jù)的觀察可知，單位壓降的增大主要是由流速增大導(dǎo)致水合物漿液在豎直管中流動摩阻增大造成的，較高的流速不利于輸送的經(jīng)濟(jì)性。

圖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 結(jié)論

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

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

符 號 說 明

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

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

CD曳力系數(shù)

dp固相顆粒粒徑，m

拖曳力，kg·m/s

G絕對速度梯度，s-1

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

g重力加速度，m/s2

H哈梅克常數(shù)

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

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

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

密度，m-3

p壓力，Pa

R發(fā)生碰撞兩水合物顆粒的調(diào)和半徑，m

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

t時間，s

u 速度矢量，m/s

u平均速度，m/s

V沉降速率，m/s

速度，m/s

相間的相對速度，m/s

z與流體性質(zhì)有關(guān)的常數(shù)

▽ 梯度

α體積分?jǐn)?shù)

β相間曳力常數(shù)

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

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

ρ密度，kg/m3

τ應(yīng)力張量，Pa

φ體積濃度

μ黏度，mPa·s

下角標(biāo)

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

l 液相

m 混合相

s 固相

TH 湍流耗散區(qū)

TG 湍流慣性區(qū)

參 考 文 獻(xiàn)

［1］ 宋廣喜，雷懷玉，王柏蒼，等. 國內(nèi)外天然氣水合物發(fā)展現(xiàn)狀與思考［J］. 國際石油經(jīng)濟(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］ 許紅，黃君權(quán)，夏斌，等. 最新國際天然氣水合物研究現(xiàn)狀及資源潛力評估（下）［J］. 天然氣工業(yè)，2005，25（6）:18-23.

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

［5］ 周守為，陳偉，李清平，等. 深水淺層非成巖天然氣水合物固態(tài)流化試采技術(shù)研究及進(jìn)展［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］ 江國業(yè)，王曉婭，孫鵬. 基于正交試驗(yàn)設(shè)計(jì)的水合物漿液流動特性數(shù)值模擬［J］. 科技導(dǎo)報(bào)，2014，32（13）:23-27.

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

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

［12］ 周守為，趙金洲，李清平，等. 全球首次海洋天然氣水合物固態(tài)流化試采工程參數(shù)優(yōu)化設(shè)計(jì)［J］. 天然氣工業(yè)，2017，37（9） :1-14.

［13］ 劉艷軍，唐孝蓉，胡坤. 天然氣水合物漿體分解對其在垂直管中流動特性影響的研究［J］. 化學(xué)通報(bào)，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］. 大連:大連理工大學(xué)，2013.

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

［17］ 趙鵬飛，王武昌，李玉星，等. 管道內(nèi)水合物漿流動的數(shù)值模型［J］. 油氣儲運(yùn)，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］ 陳鵬，劉福旺，李玉星，等. 水合物漿液流動特性數(shù)值模擬［J］. 油氣儲運(yùn)，2014，33（2）:160-164.