侯朋朋，王春華，潘 振，商麗艷，韋雪蕾，劉志明

（遼寧石油化工大學(xué)石油天然氣工程學(xué)院，遼寧撫順113001）

在中國天然氣水合物市場深度分析及“十三五”發(fā)展戰(zhàn)略報告中重點指出，要加快天然氣水合物試開采進(jìn)度，推動產(chǎn)業(yè)化進(jìn)展。目前，國際上對天然氣水合物的勘探研究一直處于探索階段，近期我國已在南?！吧窈Ｓ颉痹囬_采成功，但現(xiàn)在對于天然氣水合物開發(fā)基礎(chǔ)和技術(shù)方法的研究仍需進(jìn)一步探索[1]。

在“地層流體抽取”試采法中，將天然氣水合物從海底輸送至海面的水力輸送垂直硬管是核心設(shè)備之一。近年來，國內(nèi)外許多學(xué)者對天然氣水合物漿液的多相流動進(jìn)行了數(shù)值模擬[2-6]和實驗研究[6-11]。宋光春等[12]利用水合物顆粒數(shù)量密度的群體平衡模型（PBM），分析了該模型下的水合物顆粒變化，但未考慮氣相的直徑變化。丁麟等[13]采用雙流體模型和歐拉方法對氣體水合物漿液系統(tǒng)的多相流流動進(jìn)行了模擬，結(jié)果表明，水合物的形成與壓降和流體的相變有關(guān)，然而并沒有對相變后的流動進(jìn)行分析。Xu H L等[14]以海底天然氣水合物的固態(tài)開采法系統(tǒng)中的開采管為研究對象，對管道內(nèi)水合物漿液的固液兩相流動進(jìn)行了CFD模擬，得到管徑、流量、顆粒粒徑等參數(shù)對開采系統(tǒng)的影響規(guī)律，由于是基于固液兩相不發(fā)生相變的情況下進(jìn)行的模擬，因此對水合物漿液的流動還需考慮氣液固三相流動。

目前，國內(nèi)外對海底天然氣水合物管道輸送的相關(guān)研究文獻(xiàn)較少。由于海底天然氣水合物開采輸送過程中不斷有氣體分解出，因此本文以水合物漿液的氣相為主導(dǎo)系統(tǒng)，流動模型中代入PBM方程，考慮了水合物分解后氣相的聚集與破碎對管道內(nèi)流動特性的影響，并對漿液的湍動能耗散率、氣相的初始濃度與速度影響氣泡分布情況進(jìn)行了分析，從而對深海水合物開采管道的漿液流變特性提供了一定的基礎(chǔ)數(shù)據(jù)和技術(shù)支持。

1 數(shù)值模擬模型的建立

1.1 幾何模型及網(wǎng)格劃分

計算模型為長度3 m，管徑300 mm的垂直管，利用ICEM對模型計算域進(jìn)行網(wǎng)格劃分，采用六面體網(wǎng)格，橫截面上網(wǎng)格劃分采用“古錢幣”形式將垂直管劃分為128 250個單元的三維結(jié)構(gòu)網(wǎng)格（見圖 1）。

圖1 模型三維結(jié)構(gòu)網(wǎng)格Fig.1 3D mesh of model

1.2 邊界條件

運用Fluent軟件建立PBM模型，定義進(jìn)口邊界為速度入口，速度為漿體輸送速度，因為湍動能和湍動耗散率難以估計，用計算獲得的湍流強度和水力直徑代替，湍流強度取4%，水力直徑為0.3 m。定義出口邊界條件為充分發(fā)展的邊界條件。固液兩相參數(shù)采取文獻(xiàn)[15]中工況參數(shù)，氣相用PBM模型對直徑為0.5～2.0 mm的氣泡分五組進(jìn)行模擬。在壁面處，流體相采用壁面函數(shù)法和無滑移邊界條件，顆粒相在壁面設(shè)置有滑移條件。顆粒間碰撞歸還系數(shù)為0.95，并加以重力加速度等環(huán)境因素。

1.3 數(shù)值計算模型

1.3.1 控制方程 將垂直管內(nèi)氣液固三相流動過程簡化為恒溫體系，即不考慮能量方程，同時采用Euler多相流模型模擬垂直管內(nèi)氣液固三相間的相互作用，水合物漿液的三相看作充滿計算區(qū)域的連續(xù)介質(zhì)進(jìn)行處理，通過連續(xù)方程和動量方程求解每一相，得到控制方程分別為連續(xù)性方程式(1)和動量方程式(2):

式中，ρ 為液相密度，kg/m3；t為時間，s；ui、uj為速度分量，m/s；xi、xj分別為 x和 y坐標(biāo)，m；p為壓力，Pa；μ為動力黏度，Pa·s；τij=-ρ為 Reynolds應(yīng)力，、為對應(yīng)的脈動量。式（1）、（2）實質(zhì)上是質(zhì)量守恒方程和動量守恒方程在流體力學(xué)中的具體體現(xiàn)。1.3.2 湍流模型 流體在垂直管道水利提升時速度變化較大，同時產(chǎn)生較大的湍動能以及湍動耗散率，本文采用realizableεk-ε湍流模型。運輸方程為式（3）：

式中，k為湍動能；ε為湍動耗散率；Gk為由于平均速度梯度引起的湍動能k的產(chǎn)生項；模型常數(shù)C1ε、C2ε、Cμ、σk和 σε的 取 值 為 1.44、1.92、0.09、1.00、1.30；μt為湍動黏度，可表示成 k和 ε的函數(shù)關(guān)系式（4）:

1.3.3 相群平衡模型 氣泡直徑大小的變化與聚并及破碎過程有關(guān)，需要添加一個平衡方程來描述氣泡在體系中的平衡，該方程稱為相群平衡模型[16]。研究中氣泡間會發(fā)生連續(xù)不斷的反應(yīng)，并對氣泡大小進(jìn)行分組，這樣可對漿體中氣泡流動進(jìn)行模擬研究，從而對水合物漿液氣液固多相流動機(jī)理進(jìn)行更直觀準(zhǔn)確的模擬研究，氣泡群體平衡方程可用式(5)表示：

式中，n(v,t)為氣泡分布函數(shù)；a(v,v’)為氣泡聚并函數(shù)；b(v)為氣泡破碎函數(shù)。

群體平衡模型的求解方法有離散法、標(biāo)準(zhǔn)動量法和動量積分法。本文采用離散法，此模型第k組氣泡的輸運方程為式（6）:

式中，ρg為氣相密度，kg/m3；Bk,c為氣泡聚并生成函數(shù)；Dk,c為氣泡聚并消亡函數(shù)；Bk,b為氣泡破碎生成函數(shù)；Dk,b為氣泡破碎消亡函數(shù)。

天然氣水合物分解出的氣泡破碎聚并的原因主要考慮為湍流碰撞，文獻(xiàn)[17]中Luo與Lehr的氣泡聚并破碎模型被廣泛應(yīng)用。本文采用Luo等基于破碎機(jī)制建立的氣泡破碎模型對水合物分解出的氣泡分布大小進(jìn)行模擬。

1.4 模擬方法可靠性

1.4.1 網(wǎng)格無關(guān)性驗證 為保證CFD數(shù)值模擬結(jié)果的準(zhǔn)確性并盡量減少計算量，首先進(jìn)行網(wǎng)格無關(guān)性驗證。選取內(nèi)徑為300 mm的垂直管，網(wǎng)格尺寸為 4.0、5.0、6.0、6.5 mm。圖 2為在入口氣泡平均直徑為0.45 mm、氣含率為0.3條件下所得氣泡平均直徑沿垂直管軸向的分布。4組模擬結(jié)果具有相同的氣泡大小分布趨勢，網(wǎng)格尺寸為6.0 mm與6.5 mm的軸向平均氣泡大小基本相同，此時所得到的氣泡大小與另外兩組存在一定誤差，因此為方便計算，采取6.0 mm的網(wǎng)格尺寸。

圖2 網(wǎng)格無關(guān)性檢驗Fig.2 Grid independence test

1.4.2 實驗驗證 依據(jù)相關(guān)實驗[18]，選取氣泡直徑分布和流動壓降對數(shù)值模型的可行性進(jìn)行驗證。表1為不同氣相體積分?jǐn)?shù)與不同流速的7種工況。表2為氣泡從直井2.500 mm遞減到0.200 mm的6組氣泡直徑分布。

表1 模擬工況Table 1 Simulate condition table

表2 氣泡初始直徑分布Table 2 Initial bubble diameter distribution

表3為單位管長壓降模擬值與實驗值的對比。由表3可知，工況5-7單位管長模擬值與實驗值的變化趨勢相同且兩者間的相對誤差均小于12%。因此該模型可較好地模擬水合物在管道中流動的壓降流動變化規(guī)律。

表3 實驗和模擬條件下單位管長壓降對比Table 3 Comparison of pressure drop of unit tube length under test and simulation conditions

圖3為工況1—4氣泡直徑分布模擬值與實驗值的對比。由圖3可知，模擬得到的氣泡直徑分布與實驗值氣泡直徑分布規(guī)律相同，均呈現(xiàn)正態(tài)分布。其中在0.5 mm以下的氣泡占比較大，大于1.5 mm的氣泡占比較少，因此該數(shù)值模型可較好的模擬水合物顆粒在分解出氣相后的氣液固三相流動。

圖3 氣泡直徑分布實驗值與模擬值對比Fig.3 Comparison of experimental and simulated values of bubble diameter distribution

2 結(jié)果與分析

2.1 湍動能耗散率對流動特性的影響

圖4為氣含率0.1時，不同湍動耗散率下軸向距離氣泡個數(shù)與氣泡體積分?jǐn)?shù)分布。由圖4(a)可知，氣泡個數(shù)隨軸向距離呈現(xiàn)先增后減的趨勢，這是由于氣泡在入口處發(fā)生明顯的破碎現(xiàn)象，由于漿液在入口處受到湍動能耗散率的影響較大，使氣泡碰撞破碎強度增強，氣泡個數(shù)隨之增多；在管道軸向距離2 m處氣泡的聚并與破碎作用趨于平衡，氣泡個數(shù)保持不變。當(dāng)湍動能耗散率較低，即漿液湍動較弱時，氣泡的破碎作用相對較弱，氣泡的數(shù)量增加速率較慢，隨著湍動能耗散率增大，氣泡的破碎作用增強。圖4(b)為氣含率0.1時改變漿液湍動耗散率得到的氣泡體積分?jǐn)?shù)分布。由圖4(b)可知，由于氣含率0.1時水合物顆粒剛開始分解，因此漿液的湍動能耗散率的變化只是使氣泡大小及其分布寬窄有所變化，氣泡的最高體積分?jǐn)?shù)在0.2～0.3，變化較為平穩(wěn)。

圖4 氣含率0.1、不同湍動耗散率下軸向距離氣泡個數(shù)與氣泡體積分?jǐn)?shù)分布Fig.4 Gas holdup 0.1,the number of bubbles and the number of bubbles at axial distance under different turbulent dissipation rates

圖5 為氣含率0.3時改變湍動能耗散率計算的氣泡大小和氣泡體積分?jǐn)?shù)分布。由圖5可知，在該氣含率條件下，如果不外加機(jī)械能量，即湍動能在0.5 m2/s3時氣泡體積分?jǐn)?shù)表示的氣泡大小分布呈現(xiàn)氣泡不均的二次聚集，氣泡分布很寬；當(dāng)增大湍動能耗散率時，氣泡破碎作用增強，氣泡大小分布變?yōu)檩^小的單峰分布。當(dāng)湍動能增大到1.5 m2/s3以上時，漿液流動由氣泡不均勻分布轉(zhuǎn)變?yōu)榫鶆蚍植嫉耐牧?。由此可見，氣含率相同時，由于耗散率的不同，可導(dǎo)致流動處于不同的流動區(qū)域。

圖5 氣含率0.3、不同湍動耗散率下軸向距離氣泡個數(shù)與氣泡體積分?jǐn)?shù)分布Fig.5 Gas holdup 0.3,the number of bubbles and the number of bubbles at axial distance under differ ent turbulent dissipation rates

2.2 氣含率對流動特性的影響

對氣泡尺寸以2.500～0.200 mm遞減的方式分成5組，設(shè)置Ratio Exponent值為3.25。將水合物漿中初始?xì)夂史謩e取0.10、0.25、0.30和0.35四種工況，5組氣泡沿軸面流態(tài)分布所占?xì)庀囿w積分?jǐn)?shù)分布如圖6所示。

圖6 不同氣含率不同尺寸的氣泡沿軸向的體積分?jǐn)?shù)分布Fig.6 Axial distribution of bubbles with different gas holdup at different sizes

在氣含率為0.10和0.25的條件下，在管道入口2 m前，氣泡處于不斷聚并的狀態(tài)，氣泡Sauter直徑不斷增大，但增加趨勢變緩。在氣含率為0.35的條件下，氣泡在管道入口1 m左右即達(dá)到聚并和破碎平衡，之后Sauter直徑不再發(fā)生變化；在氣含率為0.10和0.25時氣泡破碎速率較低，由于湍流渦體的動能很小，不足以克服小氣泡的表面張力，因此在氣含率較小時氣泡破碎率較低；氣含率大于0.25時，直徑在2.5 mm的氣泡聚并速率增加較快，氣泡聚并作用增強。

圖7為在不同氣含率的條件下，管道中軸面上氣泡平均直徑分布云圖。由圖7可知，在垂直入口處的初始?xì)馀葜睆?.55 mm左右，氣泡聚并作用不明顯，由于氣泡聚并在管道中逐漸占主導(dǎo)作用，因此氣泡會迅速增大，增至1.50 mm左右，在氣泡聚并和破碎達(dá)到平衡時，氣泡直徑趨于穩(wěn)定，在管道出口處氣泡直徑變化不大。在氣含率較低時，氣泡間距離較大，以小氣泡為主，氣泡聚并作用很弱，氣泡的湍動速度不足以使氣泡相互碰撞，因此由湍動渦體引起的氣泡聚并作用基本可以忽略；當(dāng)氣含率在0.25到0.30時，氣泡個數(shù)增多，氣泡碰撞頻率和聚并頻率增大，導(dǎo)致氣泡直徑變大且分布較寬；當(dāng)氣含率進(jìn)一步增大至0.35時，由于氣泡聚并作用形成大氣泡的體積分?jǐn)?shù)明顯增大，漿液的湍動耗散率隨之增大，因此氣泡直徑大小分布沿軸向變化更快，圖7(a)中，在管道入口1/4處，氣泡直徑接近2 mm，氣相結(jié)構(gòu)與低氣含率條件下相比發(fā)生明顯轉(zhuǎn)變。

圖7 不同氣含率的管道軸向位置氣泡平均直徑Fig.7 Average diameter of bubbles in axial positions of pipes with different gas holdup

2.3 進(jìn)氣速度對氣泡變化的影響

圖8 、9分別為改變表觀氣相、表觀液相和固相速度，其他兩相不變時氣泡和漿液上升速度徑向分布的模擬計算結(jié)果。由圖8（a）和圖9（a）可以看出，表觀氣速為0.25 m/s時，氣泡和漿液上升速度沿徑向分布較為均勻，管道中心氣泡三相之間的滑移速度與單氣泡上升速度接近；當(dāng)流動進(jìn)入氣體分解量較大后，漿液速度和氣泡速度上升徑向分布不均勻性增加，三相之間的滑移速度明顯大于單氣泡流動速度，與進(jìn)氣速度成線性關(guān)系，這主要是分解出的氣泡聚并成大氣泡的加速作用造成的。由圖8（b）、（c）和圖 9（b)、（c）可以看出，相同進(jìn)氣速度下，表觀液速較高時，漿液速度和氣泡上升速度徑向分布較均勻，三相間的滑移速度與單氣泡上升速度接近，表觀液速降到1.02 m/s時，漿液速度和氣泡上升速度沿徑向分布均勻下降，三相滑移速度增大。

圖8 氣液固三相對氣泡速度徑向分布影響Fig.8 Influence of gas-liquid solid three relative bubble velocity radial distribution

圖9 氣液固三相對漿液速度徑向分布影響Fig.9 Influence of r adial distr ibution of gas-liquid solid three relative slurry velocity

3 結(jié) 論

（1）采用PBM模型模擬水合物漿液的流動，充分考慮氣泡直徑分布變化。該模型中，氣含率和漿液湍動能耗散率是影響氣泡大小分布的重要參數(shù)。當(dāng)漿液的湍動耗散率增大時，氣泡破碎速率增大，使氣泡直徑減小。

（2）氣含率的改變對氣泡碰撞與聚團(tuán)作用比較明顯，隨著氣含率的增加，Sauter平均直徑達(dá)到平衡越快，氣泡體積分?jǐn)?shù)在0.30～0.35時，發(fā)生聚團(tuán)作用顯著，使天然氣水合物漿液提升速度明顯加快。

（3）氣相的速度變化對水合物漿液的速度提升影響較大，在氣液固三相體系中，氣泡行為對體系的流體力學(xué)行為有決定性影響，改變表觀氣速、表觀液速和固相速度都會導(dǎo)致氣含率發(fā)生改變，進(jìn)而導(dǎo)致氣泡大小分布和相間作用發(fā)生改變。增大表觀氣速，中心區(qū)域的漿液速度有較大幅度的增加，并沿徑向位置遞減。