姚淑鵬 李玉星 王武昌 宋光春 施政灼 王小玉 劉帥

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

近年來，天然氣水合物以其巨大的儲(chǔ)量和綠色環(huán)保等優(yōu)勢(shì)受到了世界上越來越多的國(guó)家的關(guān)注[1-2]。天然氣水合物礦藏儲(chǔ)量非常大，據(jù)相關(guān)專家估計(jì)，世界上水合物礦藏的水合物儲(chǔ)量相當(dāng)于現(xiàn)在已經(jīng)探明的碳化合物的兩倍，其中絕大部分存在于深水區(qū)域，約占儲(chǔ)量的95%[3]。而深海天然氣水合物礦藏中，絕大部分是儲(chǔ)存于海底淺表層，以非成巖水合物礦藏的形式存在，有著埋深淺、非成巖、膠結(jié)性差的性質(zhì)。非成巖水合物的自身性質(zhì)決定了其開采有著很高的生態(tài)環(huán)境等風(fēng)險(xiǎn)及更高的開采難度[4]。為此，我國(guó)首次提出了固態(tài)流開采的非成巖水合物綠色開采方式。天然氣水合物固態(tài)流開采時(shí)，水合物在海底的環(huán)境條件下穩(wěn)定存在，通過機(jī)械采掘的方式將水合物碎化采集、二次破碎、預(yù)分離、加海水流化后進(jìn)入舉升管，以水合物漿液的形式舉升至海面的工程船上進(jìn)一步處理[5]。水合物漿液進(jìn)入舉升管后主要分為兩個(gè)階段，兩階段以水合物臨界分解點(diǎn)（即達(dá)到分解條件的位置）為界。水合物分解臨界點(diǎn)以下的舉升管段，水合物漿液未達(dá)到分解條件，以顆粒的形式與海水混合成固-液兩相的水合物漿液進(jìn)行輸送；當(dāng)水合物漿液輸送至水合物分解臨界點(diǎn)以上管段時(shí)，水合物顆粒會(huì)發(fā)生相變分解，產(chǎn)生氣態(tài)天然氣，此時(shí)漿液為固-液-氣三相流動(dòng)。

水合物漿液在豎直管道內(nèi)的研究對(duì)水合物固態(tài)流開采有著重要的意義，但其實(shí)驗(yàn)的成本太高，不適宜推廣使用。近年來，隨著計(jì)算流體力學(xué)（CFD）的發(fā)展，運(yùn)用CFD 數(shù)值模擬技術(shù)進(jìn)行水合物漿液管內(nèi)多相流動(dòng)研究有效地解決了實(shí)驗(yàn)成本較大、難以開展的問題。BALAKIN 等[6-7]對(duì)R11 水合物在湍流純水體系中的流動(dòng)沉積現(xiàn)象進(jìn)行了CFD模擬研究。FATNES[8]通過ANSYS CFX 軟件模擬了水平管內(nèi)水合物漿液的流動(dòng)特性。江國(guó)業(yè)等[9]通過正交模擬試驗(yàn)，對(duì)彎管內(nèi)水合物漿液流速、水合物顆粒粒徑及水合物體積分?jǐn)?shù)對(duì)水合物漿液流動(dòng)特性的影響進(jìn)行了研究。宋光春等[10-13]通過構(gòu)建群體相平衡模型對(duì)R11 水合物漿液在水平管中的聚集、破碎和沉積等流動(dòng)特性進(jìn)行了模擬研究。但是，當(dāng)前運(yùn)用CFD 技術(shù)對(duì)水合物漿液的研究主要集中在水平管，而與水合物開采相關(guān)的豎直管內(nèi)模擬研究還較少。周守為等[5]通過有限差分迭代的模擬對(duì)豎直井筒內(nèi)的水合物的固相和氣相含量的變化等水合物漿液流動(dòng)特性進(jìn)行了研究。劉艷軍等[14]通過CFD 模型中的Euler 多相流模型以及Finite-Rate/Eddy-Dissipation 模型對(duì)天然氣水合物漿體在垂直管輸中的分解現(xiàn)象對(duì)流動(dòng)特性的影響進(jìn)行了研究。上述學(xué)者對(duì)水合物漿液在豎直管內(nèi)的部分流動(dòng)特性進(jìn)行了分析研究，但其并沒有考慮水合物顆粒在流動(dòng)過程中的聚集和破碎，也不能較好地反映各變量對(duì)水合物漿液流動(dòng)特性影響的強(qiáng)弱，需要進(jìn)一步研究完善。

本文通過基于水合物顆粒聚集動(dòng)力學(xué)的群體平衡模型和正交試驗(yàn)的方法對(duì)不同工況下的水合物漿液在豎直管內(nèi)的流動(dòng)特性進(jìn)行了模擬，主要研究了水合物固態(tài)流開采時(shí)，水合物漿液在臨界分解點(diǎn)以前的固液兩相流動(dòng)特性，并選取了水合物初始粒徑、水合物體積分?jǐn)?shù)以及水合物漿液流速三個(gè)因素進(jìn)行分析。研究結(jié)果可以為水合物固態(tài)流開采的輸送研究提供借鑒。

1 模型建立

1.1 幾何模型

本次模擬研究中，為更好地對(duì)所構(gòu)建模型進(jìn)行驗(yàn)證，基于Balakin 等開展的管內(nèi)CCl3F（R11）水合物漿流動(dòng)特性實(shí)驗(yàn)所使用管道的相關(guān)參數(shù)進(jìn)行了三維模型的構(gòu)建。最終，所建立的管道長(zhǎng)2.0 m，直徑4.52 cm。

模型構(gòu)建后進(jìn)行網(wǎng)格的劃分。由于流體在管道邊界區(qū)域處物理參數(shù)的變化相對(duì)比較劇烈，因此這里對(duì)進(jìn)口壁面進(jìn)行了邊界層效應(yīng)的加密處理，共設(shè)了6 個(gè)邊界層。除此之外的其他網(wǎng)格則均以1 mm劃分。最終，共建立了170 180 個(gè)六面體網(wǎng)格，網(wǎng)格的質(zhì)量為0.917。幾何模型的網(wǎng)格圖如圖1 所示。

圖1 幾何模型網(wǎng)格圖Fig.1 Mesh of geometric model

1.2 模型基本假設(shè)

在建模過程中有5 個(gè)基本假設(shè)，即：①認(rèn)為流體是不可壓縮介質(zhì)；②忽略水合物漿內(nèi)相間質(zhì)量傳遞，即不考慮水合物的生成及分解；③水合物漿液在立管內(nèi)流動(dòng)過程等溫；④不考慮水合物在管壁上的粘附；⑤認(rèn)為水合物顆粒均為連續(xù)性介質(zhì)。

1.3 多相流模型及湍流模型

本次模擬研究中，首先對(duì)多相流模型進(jìn)行了比選，最終選擇歐拉-歐拉雙流體模型，其包含了連續(xù)性方程、動(dòng)量方程和用來封閉方程組的本構(gòu)方程三部分。其中連續(xù)性方程和動(dòng)量方程分別為

式中：t為時(shí)間，s；i為不同相態(tài)；α為體積分?jǐn)?shù)；ρ為密度，kg/m3；?為拉普拉斯算子；u為速度矢量，m/s；τi為應(yīng)力張量，Pa；p為壓力，Pa；Mi為相間動(dòng)量交換項(xiàng)，kg/(m·s)2。

本次模擬研究為固液兩相流動(dòng)，流動(dòng)過程中不存在水合物的相變，因此模型中需要重點(diǎn)考慮液固的耦合。FLUENT 是通過相間動(dòng)量交換來實(shí)現(xiàn)液固耦合計(jì)算的，本次模擬在計(jì)算相間動(dòng)量交換的選取中主要考慮了相間拖曳力，其主要由形狀阻力和摩擦阻力構(gòu)成。模擬的相間拖曳力選取Gidaspow 模型[15]進(jìn)行計(jì)算。

其中，當(dāng)αs≤20%時(shí)，選用Wen-Yu 公式計(jì)算，其如式（3）所示

當(dāng)αs＞20%時(shí)，選用Ergun 公式計(jì)算，其如式（4）所示

則相間拖曳力可以由式（5）求得

式 中：β為相間曳力常數(shù)，kg/(m3·s)；分別是液相、固相的速度，m/s；φs為固相體積濃度；dp為固相顆粒粒徑，m；CD為曳力系數(shù)；為拖曳力，N；為相間的相對(duì)速度，m/s。

在模型的構(gòu)建中，水合物顆粒的黏度也需要進(jìn)行求解，求解公式為[16]

其中，R11 水合物漿液和天然氣水合物漿液的表面黏度μm分別由Roscoe-Brinkman 方程（7）[17]和Thomas 公式（8）[18]計(jì)算

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

模型中的湍流模型選用FLUENT 15.0 軟件中的標(biāo)準(zhǔn)k-ε模型，此處不再詳細(xì)介紹。

1.4 群體平衡模型

在群體平衡模型構(gòu)建時(shí)，根據(jù)本文的5 個(gè)假設(shè)，當(dāng)不存在水合物顆粒的相變以及水合物顆粒不在管壁粘附，并認(rèn)為水合物顆粒的粒徑為連續(xù)分布時(shí)，群體平衡方程可由式（9）表述。

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

本次模型中碰撞頻率主要考慮由差速沉降和流動(dòng)剪切造成的碰撞，且取兩者碰撞頻率的和作為水合物顆粒間的實(shí)際碰撞頻率。其中，選取CAMP等[19]提出的公式來計(jì)算差速沉降碰撞頻率βDS。

式中：V為沉降速率，模型中采用式（11）計(jì)算可得

對(duì)模型中流動(dòng)剪切碰撞頻率進(jìn)行計(jì)算，當(dāng)水合物顆粒粒徑小于Kolmogorov 尺度時(shí)，其處于湍流耗散區(qū)，其聚集主要受到渦內(nèi)局部剪切力影響。這里，選取SAFFMAN 等[20]提出的公式計(jì)算水合物顆粒的碰撞頻率

當(dāng)水合物顆粒大于Kolmogorov 尺度時(shí)，其處于湍流慣性區(qū)，受到主流場(chǎng)牽引而運(yùn)動(dòng)。這里，選取ABRAHAMSON 等[21]提出的公式計(jì)算水合物顆粒的碰撞頻率

以上兩式中：u為水合物顆粒的平均速度，m/s；G為流場(chǎng)局部的絕對(duì)速度梯度，s-1。

模型中選用曲線模型來計(jì)算水合物顆粒間的聚并效率。本次模擬中由于連續(xù)相是水，水合物顆粒間沒有液橋力作用，因此聚并效率計(jì)算時(shí)主要考慮范德華力和流動(dòng)剪切力之比。聚并效率可由式（14）計(jì)算[22]

式中：H為表征范德華力大小的哈梅克常數(shù)；k為表示兩顆粒發(fā)生碰撞后破碎可能性的參數(shù)；R為發(fā)生碰撞兩水合物顆粒的調(diào)和半徑。

模型構(gòu)建時(shí)，水合物顆粒破碎效應(yīng)選取FLUENT 15.0 軟件中的Lehr 模型，此處不再詳細(xì)介紹。

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

1.5 模型求解

本文模型利用FLUENT 15.0 軟件進(jìn)行求解。邊界條件設(shè)置：速度進(jìn)口、壓力出口（壓力為0）、壁面無滑移。模擬中具體的流速等參數(shù)下文會(huì)詳細(xì)列出。除此之外，模擬選取二階迎風(fēng)格式；且當(dāng)各因子的殘差收斂到10-5時(shí)，認(rèn)為已經(jīng)收斂，模擬過程結(jié)束。部分模擬參數(shù)如表1 所示，模擬工況如表2 所示。

表1 模型參數(shù)Tab.1 Parameters of model

表2 模擬工況參數(shù)Tab.2 Parameters of simulation condition

1.6 模型驗(yàn)證

目前，對(duì)水合物漿液在立管中流動(dòng)的實(shí)驗(yàn)數(shù)據(jù)仍較為缺乏，因此本文通過間接驗(yàn)證的方式對(duì)模型進(jìn)行了驗(yàn)證。首先采用BALAKIN 等[6-7]在水平環(huán)道內(nèi)R11 水合物漿液流動(dòng)研究中所得的實(shí)驗(yàn)數(shù)據(jù)對(duì)模型進(jìn)行了驗(yàn)證。驗(yàn)證中選取了環(huán)道內(nèi)流動(dòng)的單位壓降作為驗(yàn)證指標(biāo)。驗(yàn)證模擬的工況參數(shù)如表3 所示，驗(yàn)證結(jié)果如表4 所示。同時(shí)，為了驗(yàn)證模型中重力影響計(jì)算的準(zhǔn)確性，采用GILLIES 等[23]在豎直管內(nèi)進(jìn)行的砂和水的固液兩相流動(dòng)實(shí)驗(yàn)來進(jìn)行對(duì)比驗(yàn)證，并選取速度分布作為指標(biāo)。驗(yàn)證所用幾何模型為長(zhǎng)2.0 m、直徑0.103 m 的管道，創(chuàng)建方式和圖1 所示幾何模型相同。驗(yàn)證模擬的工況如表5 所示，驗(yàn)證結(jié)果如圖2 所示。經(jīng)過驗(yàn)證，本文所構(gòu)建的流動(dòng)模型能夠較為準(zhǔn)確地模擬水合物在豎直管內(nèi)的流動(dòng)特性。

表3 水合物流動(dòng)模型驗(yàn)證的工況參數(shù)Tab.3 Working condition parameters for hydrate flow model verification

表4 壓降驗(yàn)證結(jié)果Tab.4 Results of pressure drop verification

表5 重力作用驗(yàn)證的工況參數(shù)Tab.5 Working condition parameters of gravity action verification

圖2 速度分布實(shí)驗(yàn)和模擬結(jié)果對(duì)比Fig.2 Speed distribution experiment and comparsion of simulation results

2 流動(dòng)特性及正交試驗(yàn)研究

天然氣水合物漿液在豎直管內(nèi)流動(dòng)特性的研究對(duì)水合物固態(tài)流開采有著重要意義。水合物固態(tài)流開采過程中，水合物漿液在舉升管內(nèi)的流動(dòng)存在臨界分解點(diǎn)。水合物漿液在臨界分解點(diǎn)之前為固液兩相流動(dòng)。下文將通過FLUENT 對(duì)水合物漿液在豎直管道內(nèi)的流動(dòng)特性進(jìn)行模擬研究和正交試驗(yàn)分析，研究粒徑和濃度分布特性，以及初始粒徑、平均流速、水合物濃度對(duì)水合物漿液的流動(dòng)摩阻和最大聚集粒徑的影響程度。

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

在本文的研究中，雖然為固-液兩相流動(dòng)（不存在水合物顆粒的生成和分解），但由于流動(dòng)中存在水合物顆粒的聚集和破碎，水合物顆粒在流動(dòng)過程中會(huì)存在粒徑分布的變化，對(duì)水合物漿液的黏度等特性造成影響。

圖3 所示為工況3 出口截面水合物漿液分布云圖，由圖3 可以發(fā)現(xiàn)，水合物顆粒的粒徑在整個(gè)出口截面上呈現(xiàn)出均勻?qū)ΨQ的分布狀態(tài)，其他工況也有著類似的分布情況。這是由于在豎直管道中，水合物漿液的流動(dòng)方向與重力方向可認(rèn)為是在同一條直線上，因此重力對(duì)整個(gè)截面上水合物顆粒的作用是均勻分布的，進(jìn)而使得水合物顆粒粒徑分布在整個(gè)截面上是對(duì)稱的。除此之外，從圖3 和圖4 可以發(fā)現(xiàn)，各工況下的水合物粒徑均呈現(xiàn)出近壁面處水合物顆粒粒徑較大且粒徑變化梯度較大，而管中心部分水合物顆粒粒徑較小且粒徑較為均勻的分布規(guī)律。這可以通過水合物顆粒的聚集動(dòng)力學(xué)進(jìn)行分析。當(dāng)水合物漿液在豎直管道內(nèi)流動(dòng)時(shí)，流動(dòng)剪切是影響水合物碰撞聚集的主要因素。在管道近壁面區(qū)域水合物漿液流動(dòng)剪切作用最強(qiáng)，而管中心區(qū)域流動(dòng)剪切作用較弱，所以水合物顆粒在近壁面處發(fā)生碰撞聚集的概率要大于管中心處，因此水合物顆粒的粒徑分布呈現(xiàn)出中間低四周高的狀況。

圖3 出口截面水合物顆粒粒徑分布云圖（工況3）Fig.3 Cloud chart of hydrate particle size distribution in outlet section(working condition 3)

圖4 各工況下的出口截面粒徑分布Fig.4 Particle size distribution of outlet section under all kinds of working conditions

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

豎直管內(nèi)水合物漿液流動(dòng)安全研究中，水合物濃度分布為重要的研究?jī)?nèi)容之一。為此，結(jié)合模擬結(jié)果對(duì)水合物漿液濃度分布的特性進(jìn)行分析討論。

由圖5、圖6 可知，本次模擬各工況模擬結(jié)果均為均勻懸浮流，即水合物在整個(gè)管截面上均勻?qū)ΨQ分布。在豎直管道中，水合物漿液的流動(dòng)方向與重力方向處在一條直線上，因此不會(huì)出現(xiàn)水平管中水合物漿液分布上下不對(duì)稱的情況[24]。這是由于豎直管道管中心區(qū)域流速高而近壁面區(qū)域流速低的速度分布導(dǎo)致的。較高的流速會(huì)增大水合物顆粒的分散系數(shù)，從而導(dǎo)致中間濃度較低且分布較均勻，而近壁面處因?yàn)榱魉俳档蜁?huì)出現(xiàn)水合物顆粒的堆積，使?jié)舛仍龃蟆?/p>

圖5 出口截面水合物顆粒濃度分布云圖（工況3）Fig.5 Cloud chart of hydrate particle concentration distribution in outlet section(working condition 3)

2.3 豎直管內(nèi)水合物漿液流動(dòng)正交試驗(yàn)

在本次模擬研究中考慮的試驗(yàn)因子為水合物初始粒徑、體積分?jǐn)?shù)和平均流速，并對(duì)每個(gè)因子各選取了三個(gè)水平，如表6 所示。在分析研究中，選取對(duì)流動(dòng)的安全性和經(jīng)濟(jì)性有著重要影響的水合物漿液的流動(dòng)壓降和最大粒徑作為評(píng)價(jià)指標(biāo)。

圖6 各工況下的出口截面濃度分布Fig.6 Concentration distribution of outlet section under all kings of working conditions

表6 水合物流動(dòng)特性影響因子水平Tab.6 Impact factor level of hydrate flow behavior

根據(jù)模擬結(jié)果整理出表7、表8，并根據(jù)表7、表8 中所得的結(jié)果進(jìn)行計(jì)算得到表9、表10，表9、表10 中的極差代表了這個(gè)因素對(duì)評(píng)價(jià)指標(biāo)的影響強(qiáng)度，極差越大，其代表的因素對(duì)評(píng)價(jià)指標(biāo)的影響越大。因此分析可得，對(duì)豎直管內(nèi)水合物漿液流動(dòng)摩阻影響最大的因素為水合物漿液的平均流速，其在生產(chǎn)運(yùn)行時(shí)應(yīng)著重考慮，在保證輸送安全的基礎(chǔ)上選取較低的輸送流速能有效增加輸送的經(jīng)濟(jì)性，因素影響的強(qiáng)弱排列是：平均流速、水合物初始粒徑、體積分?jǐn)?shù)。同時(shí)，對(duì)于豎直管內(nèi)最大粒徑影響最大的為漿液中水合物的體積分?jǐn)?shù)，因此在保證漿液輸送的流動(dòng)安全時(shí)水合物濃度應(yīng)該著重考慮，因素影響的強(qiáng)弱排列是：體積分?jǐn)?shù)、水合物初始粒徑、平均流速。

表7 水合物流動(dòng)特性研究評(píng)價(jià)指標(biāo)（壓降）Tab.7 Hydrate flow behavior research and evaluation index(pressure drop)

表8 水合物流動(dòng)特性研究評(píng)價(jià)指標(biāo)（粒徑）Tab.8 Hydrate flow behavior research and evaluation index(particle size)

表9 水合物影響因素極差分析（壓降）Tab.9 Range analysis of hydrate influencing factors(pressure drop)Pa/m

表10 水合物影響因素極差分析（最大粒徑）Tab.10 Range analysis hydrate influencing factors maximum(particle size)mm

3 結(jié)論

水合物漿液在豎直管內(nèi)的流動(dòng)特性的研究對(duì)固態(tài)流開采的安全性和經(jīng)濟(jì)性有著重大意義。本文引入了基于水合物顆粒聚集動(dòng)力學(xué)的群體平衡模型和正交試驗(yàn)方法，通過CFD/FLUENT 模擬對(duì)水合物漿液在豎直管內(nèi)的流動(dòng)特性進(jìn)行了研究。通過對(duì)模擬結(jié)果的分析，主要得到以下結(jié)論：

（1）當(dāng)水合物漿液在豎直管內(nèi)流動(dòng)時(shí)，水合物顆粒的粒徑在整個(gè)截面上是對(duì)稱分布的，且存在近壁面區(qū)域水合物顆粒粒徑較大，管中心區(qū)域水合物顆粒粒徑較小的分布狀況。

（2）當(dāng)水合物漿液在豎直管內(nèi)流動(dòng)時(shí)，水合物的濃度在整個(gè)截面上對(duì)稱分布，且存在近壁面區(qū)域水合物濃度較高，而管中心區(qū)域水合物濃度較低、較為均勻的分布狀況。

（3）當(dāng)水合物漿液在豎直管內(nèi)流動(dòng)時(shí)，本次研究所考慮的三個(gè)因素中，對(duì)豎直管內(nèi)水合物漿液流動(dòng)摩阻影響最大的因素為水合物漿液的平均流速，對(duì)豎直管內(nèi)最大粒徑影響最大的因素為漿液中水合物的體積分?jǐn)?shù)。