基于OpenFOAM的豎直管內(nèi)氣-液-水合物流動(dòng)特性研究

摘 要 針對(duì)天然氣水合物開采過程中涉及的復(fù)雜氣液固三相流動(dòng)問題，基于群體平衡模型，提出了新的水合物顆粒聚并破碎模型組合，同時(shí)考慮了氣泡的聚并破碎行為，運(yùn)用開源CFD軟件OpenFOAM對(duì)建立的模型進(jìn)行求解驗(yàn)證?？疾炝巳肟诹魉?、氣相和水合物顆粒體積分?jǐn)?shù)對(duì)豎直管內(nèi)氣-液-水合物流動(dòng)特性的影響。結(jié)果表明：所建立的模型準(zhǔn)確有效，各相的濃度分布、氣泡的平均粒徑徑向分布和水合物顆粒平均粒徑大小受水合物顆粒體積分?jǐn)?shù)的影響很大，壓力損失梯度隨著流速的增加先減小后增大，存在特定工況下的最優(yōu)流速，氣相體積分?jǐn)?shù)的增加會(huì)讓管內(nèi)流體的壓力損失梯度明顯減小。

關(guān)鍵詞 天然氣水合物 豎直管 多相流 群體平衡模型 OpenFOAM

中圖分類號(hào) TE311" "文獻(xiàn)標(biāo)識(shí)碼 A" "文章編號(hào) 0254?6094（2023）01?0067?09

天然氣水合物（俗稱“可燃冰”）是在高壓低溫條件下，由甲烷等氣體和水結(jié)合而成的類冰狀物。自然界中約98%的天然氣水合物資源聚集在海底沉積物中，約2%存在于永久凍土帶。相較于傳統(tǒng)化石能源，天然氣水合物清潔無污染、能量密度高、資源潛力巨大[1]。對(duì)天然氣水合物的勘探、開發(fā)和利用具有重要的環(huán)境、資源和戰(zhàn)略意義。目前，天然氣水合物的開采方法主要包括降壓法、加熱法、化學(xué)劑法、CO置換法及固體開采法等。自2011年起，我國(guó)先后在祁連山木里地區(qū)、南海神狐等地區(qū)進(jìn)行了5次天然氣水合物試驗(yàn)性開采，累計(jì)產(chǎn)氣量達(dá)1.17×106 m3[2]。無論采用哪種開采方式，要將所開采的天然氣或固體水合物顆粒輸送至海平面，管道中都會(huì)出現(xiàn)天然氣、水和水合物共存的氣液固三相流動(dòng)，使天然氣水合物輸送過程十分復(fù)雜[3]。因此，深入研究豎直管內(nèi)氣-液-水合物流動(dòng)特性對(duì)天然氣水合物的開發(fā)和利用是十分必要的。

近年來，隨著計(jì)算流體力學(xué)（CFD）的發(fā)展，數(shù)值模擬方法被廣泛應(yīng)用于多相流的研究中。眾多學(xué)者利用數(shù)值模擬方法，對(duì)管道內(nèi)水合物流動(dòng)特性進(jìn)行了研究。LI L等分別采用歐拉模型和CFD?PBM模型模擬較大氣相體積分?jǐn)?shù)下天然氣水合物在水力提升管道中的三相流，數(shù)值模擬和實(shí)驗(yàn)對(duì)比結(jié)果表明，CFD?PBM模型模擬的氣泡流動(dòng)狀態(tài)比歐拉模型更加均勻，流速和氣體濃度分布與實(shí)驗(yàn)值吻合較好[4]。劉艷軍等采用Euler模型對(duì)天然氣水合物漿體垂直管輸?shù)墓?液兩相流、氣-液-固三相流流動(dòng)特性進(jìn)行研究，發(fā)現(xiàn)天然氣水合物顆粒的速度分布、濃度分布隨高度的變化呈現(xiàn)出波動(dòng)-均勻-波動(dòng)的規(guī)律，水合物分解對(duì)漿體管道運(yùn)輸具有減阻作用[5]。YAO S P等采用種群平衡理論模擬了垂直管道中水合物漿體在達(dá)到臨界分解點(diǎn)之前的固-液兩相流動(dòng)，分析了不同工況下垂直管道中壓降、相分布和水合物粒度分布，發(fā)現(xiàn)水合物顆粒的粒徑和粒徑變化梯度在近壁面區(qū)域較大，水合物顆粒的最大粒徑隨水合物濃度的增加而增大，隨流速的增加而減小，壓力損失梯度隨流速和濃度的增加而增大，隨初始粒徑的增大而減小[6]。CHEN W等研究了水合物分解作用下垂直管道內(nèi)氣液固流動(dòng)行為，分析了管道直徑、入口速度和固相參數(shù)對(duì)天然氣水合物管道三相流動(dòng)的影響，結(jié)果表明，隨著管徑的增大，管道壓力損失梯度減小，管道入口流速對(duì)管道內(nèi)流場(chǎng)的流型結(jié)構(gòu)有很大的影響[7]。CAO X W等

建立了水平和上斜管道內(nèi)氣-液-水合物多相紊流系統(tǒng)的數(shù)學(xué)模型，主要考慮了水合物顆粒的碰撞、團(tuán)聚和破碎機(jī)制，分析了氣體流速、料漿流速、水合物體積分?jǐn)?shù)和管道傾角對(duì)流型和粒徑分布的影響，發(fā)現(xiàn)水合物顆粒傾向于團(tuán)聚，并分布在管道底層，水合物體積分?jǐn)?shù)和水合物漿體流速的增加會(huì)引起流型轉(zhuǎn)變和顆粒團(tuán)聚強(qiáng)化，水合物的粒度分布在向上傾斜的部分更加對(duì)稱，最終達(dá)到團(tuán)聚破碎的平衡狀態(tài)[8]。上述研究雖然考察了兩相或三相流情況下，入口流速、各相體積分?jǐn)?shù)及管道直徑等因素對(duì)流動(dòng)的影響，但沒有考慮到氣泡和水合物顆粒粒徑共同變化的情況。

筆者基于群個(gè)體平衡模型提出了新的聚并破碎模型組合，在開源CFD軟件OpenFOAM中植入并進(jìn)行求解驗(yàn)證，在考慮氣泡和水合物顆粒聚并破碎行為同時(shí)發(fā)生的情況下，模擬研究豎直管內(nèi)氣-液-水合物三相流的流動(dòng)特性。重點(diǎn)分析了入口流速、氣相體積分?jǐn)?shù)和水合物顆粒體積分?jǐn)?shù)對(duì)豎直管內(nèi)氣泡和水合物顆粒的平均粒徑分布、濃度分布和壓力損失梯度的影響，為深海天然氣水合物的安全高效開發(fā)和利用提供幫助。

1 數(shù)值計(jì)算模型

1.1 雙流體模型

雙流體模型（Two?Fluid Model，TFM）也被稱為歐拉-歐拉模型，模型中將連續(xù)相和分散相均視為連續(xù)介質(zhì)，各相遵循各自的連續(xù)性方程和動(dòng)量方程。雙流體模型是描述多相系統(tǒng)的宏觀模型，由于計(jì)算量相對(duì)較小，計(jì)算效率較高，目前在工業(yè)模擬中得到廣泛應(yīng)用。為了便于研究氣-液-水合物多相流動(dòng)特性與流動(dòng)參數(shù)的關(guān)系，在建模過程中不考慮水合物的生成和分解、氣相在液相中的溶解和析出、環(huán)境溫度變化的影響，因此氣-液-水合物之間無傳熱傳質(zhì)過程。其瞬時(shí)控制方程如下：

連續(xù)性方程

+▽·（αρ[u][→]）=0 （1）

動(dòng)量方程

+▽·（α ρ[u][→][u][→]）=-α▽p+▽·（ατ）+αρ[g][→]+[F][→]（2）

其中，α表示體積分?jǐn)?shù);ρ表示密度;φ表示不同的相，φ=d表示分散相，φ=c表示連續(xù)相;[u][→]表示相應(yīng)相瞬時(shí)運(yùn)動(dòng)速度;p為局部壓力;[g][→]為重力加速度;τ為粘性應(yīng)力;[F][→]為相間作用力。由于管內(nèi)氣體和水合物含率較低，因此氣泡和水合物顆粒為分散相，海水為連續(xù)相。

在雙流體模型中氣-液-水合物相的耦合是通過相間作用力實(shí)現(xiàn)的，因此相間作用力是建模的重點(diǎn)內(nèi)容。為了使模擬更加接近實(shí)際情況，本課題考慮的相間作用力包括曳力、升力、壁面潤(rùn)滑力、湍流分散力和虛擬質(zhì)量力。

曳力是各相間最重要的作用力，單位體積的分散相所受曳力F為：

F=K[u][→]（3）

其中，K表示動(dòng)量傳遞系數(shù);[u][→]表示相間相對(duì)速度。由于氣液固三相間的物性差別，動(dòng)量傳遞系數(shù)模型的選擇也有所不同，氣-液相間采用Tomiyama模型[9]，即：

K=α α ρ|[u][→]|（4）

式中 C——曳力系數(shù);

d——分散相粒徑大小;

K——Tomiyama模型動(dòng)量傳遞系數(shù)。

液-水合物和氣-水合物相間選用Gidaspow模型[10]，即當(dāng)分散相體積分?jǐn)?shù)小于20%時(shí)，由

Wenamp;Yu模型計(jì)算;分散相體積分?jǐn)?shù)大于20%時(shí)，由Ergun模型計(jì)算。具體如下：

K=αα ρ|[u][→]|（1-α） （5）

KDE=150+1.75|[u][→]| （6）

式中 d——連續(xù)相粒徑大小;

K、K——Wenamp;Yu模型和Ergun模型動(dòng)量傳遞系數(shù);

μ——連續(xù)相粘度。

除了曳力外，分散相在流體中運(yùn)動(dòng)還會(huì)受到與其運(yùn)動(dòng)方向相垂直的升力F：

F=Cαρ[u][→]×（▽×[u][→]）（7）

式中 C——升力系數(shù)。

許多研究人員針對(duì)不同工況提出了多種升力系數(shù)模型，文中選取WANG T F等[11]基于CFD回歸法得到的模型，因OpenFOAM中無此模型，故需要單獨(dú)植入。

靠近壁面的粒子周圍流場(chǎng)是非對(duì)稱的，會(huì)產(chǎn)生一種將粒子推離壁面的力，這種力稱為壁面潤(rùn)滑力F：

F=Cα ρ|[u][→]-（[u][→]·[η][→]）|[η][→]（8）

式中 C——壁面潤(rùn)滑力系數(shù)，采用Frank模型[12];

[η][→]——壁面的單位法向量。

湍流分散力F表示連續(xù)相湍動(dòng)對(duì)分散相的影響，其作用效果是使分散相沿管道徑向分布更加均勻，表達(dá)式如下：

F=K▽?duì)粒?）

其中，K與所選擇的湍流分散力模型有關(guān)，氣-液相間采用Lopez de Bertodano模型[13]，液-固相間使用Burns模型[14]。

虛擬質(zhì)量力F也稱為附加質(zhì)量力：

F=Cα ρ（10）

其中，C為虛擬質(zhì)量力系數(shù)。對(duì)于球形粒子，氣相和水合物顆粒的虛擬質(zhì)量力系數(shù)分別取0.25和0.5。

1.2 湍流模型和顆粒動(dòng)力學(xué)理論

本研究中氣液水合物三相的湍流模型均選用目前廣泛使用的標(biāo)準(zhǔn)k?ε模型，水合物顆粒的粘度則采用顆粒動(dòng)力學(xué)理論（Kinetic Theory of Granular Flow，KTGF）來計(jì)算，這里不再詳細(xì)介紹。

1.3 群體平衡模型

群體平衡模型（Population Balance Model，PBM）是用來描述多相體系中分散相粒徑隨時(shí)間和空間的分布情況，目前得到廣泛應(yīng)用，RAMKRISHNA D[15]給出了較為通用的形式：

其中，n（v）是數(shù)密度函數(shù)，其具體數(shù)值與變量v有關(guān);[u][→]是分散相速度;a（v，v′）表示聚并速率;b（v）表示破碎速率;β（v，v′）是子粒子分布函數(shù)。式（11）右邊涉及的聚并速率、破碎速率和破碎后子粒子分布函數(shù)是影響求解數(shù)密度函數(shù)的重要因素。聚并和破碎過程主要受分散相周圍流場(chǎng)的影響，又與分散相物性息息相關(guān)，涉及到非常復(fù)雜的多層次多尺度物理機(jī)制。

筆者采用概率模型來得到水合物顆粒的聚并速率a，水合物顆粒之間的聚并速率等于其碰撞頻率h（d，d）和聚并效率λ（d，d）的乘積，即：

a=h（d，d）·λ（d，d）" （12）

水合物顆粒之間的碰撞是多種機(jī)制引起的，總的碰撞頻率為所有機(jī)制引起的碰撞頻率之和，由于實(shí)際輸送的水合物顆粒粒徑較大，筆者主要考慮湍流隨機(jī)脈動(dòng)和差速沉降兩種機(jī)制。在雷諾數(shù)較大的情況下湍流隨機(jī)脈動(dòng)是引起碰撞的主要機(jī)制，當(dāng)水合物顆粒尺寸大于Kolmogorov尺度時(shí)，水合物顆粒間的碰撞頻率可用ABRAHAMSON J[16]提出的模型進(jìn)行計(jì)算：

h=2π（R+R）（13）

式中 R——兩個(gè)相互碰撞水合物顆粒的半徑;

U——兩個(gè)相互碰撞水合物顆粒的平均速度的平方。

當(dāng)水合物顆粒尺寸小于Kolmogorov尺度時(shí)，水合物顆粒間的碰撞頻率采用SAFFMAN P G和TURNER J S[17]提出的模型進(jìn)行計(jì)算：

h=

（14）

式中 d——水合物顆粒直徑;

ε——湍動(dòng)能耗散率;

ν——運(yùn)動(dòng)粘度。

差速沉降引起的碰撞頻率h采用CAMP T R和STEIN P C[18]提出的模型進(jìn)行計(jì)算：

h=（d+d）|V-V|（15）

V=R（16）

式中 R——粒子的半徑;

V——水合物顆粒瞬時(shí)沉降速度;

μ——連續(xù)相動(dòng)力粘度。

確定碰撞頻率模型后還需考慮聚并效率，文中聚并效率選取VEN T和MASON S G[19]提出的曲線模型此，模型在水合物漿液流動(dòng)的研究中已得到廣泛的使用[6，8，20]，如下：

λ（d，d）=k

（17）

R=（18）

式中 H——哈梅克常數(shù);

k——與流體性質(zhì)有關(guān)的參數(shù);

R——當(dāng)量半徑。

式（17）中相關(guān)參數(shù)取值參考SONG G C等[20]的工作。

水合物顆粒的破碎速率取決于固體顆粒的大小和最小的湍流渦旋，本研究選取MARCHISIO D L等[21]提出的破碎速率模型：

b（v）=cvεd" （19）

其中，c、c、c、c為常數(shù)，分別取值0.000 6、-1.25、0.75、1.0。

子粒子分布函數(shù)的形式取決于顆粒的強(qiáng)度、形態(tài)等特性及破碎機(jī)理等多種因素，文中水合物顆粒破碎后的子粒子分布函數(shù)選用均勻分布。

氣泡的聚并速率和破碎速率均采用LEHR F等[22]提出的模型，其中氣泡的破碎形式為二元破碎。

1.4 模型求解及參數(shù)設(shè)置

本研究使用開源CFD軟件OpenFOAM中的multiphaseEulerFoam求解器，將式（7）升力模型和水合物顆粒的聚并速率、破碎速率模型植入到其中，最終完成TFM?PBM耦合模型求解。文中建立了長(zhǎng)10 m、直徑300 mm的三維豎直圓管的幾何模型，對(duì)幾何模型進(jìn)行六面體網(wǎng)格劃分，近壁面處進(jìn)行加密處理，共劃分為481 600個(gè)網(wǎng)格，如圖1所示。

為了得到比較合理的結(jié)果，筆者將氣泡和水合物顆粒尺寸從小到大劃分為15個(gè)區(qū)間（表1、2）。模擬過程中主要參數(shù)和模擬工況分別見表3、4。

2 結(jié)果與討論

2.1 模型驗(yàn)證

筆者根據(jù)BALAKIN B V等[23]的實(shí)驗(yàn)結(jié)果，對(duì)液固兩相流中水合物顆粒粒徑分布和壓力損失梯度的數(shù)值模擬結(jié)果的有效性進(jìn)行了驗(yàn)證;根據(jù)WANG T F等[24]的實(shí)驗(yàn)結(jié)果，對(duì)氣液固三相流中氣泡粒徑分布數(shù)值模擬結(jié)果的有效性進(jìn)行了驗(yàn)證。由圖2可知，水合物顆粒粒徑分布、壓力損失梯度和氣泡粒徑分布的模擬結(jié)果與實(shí)驗(yàn)結(jié)果分布規(guī)律基本相同。因此，所建立的各類模型有效，應(yīng)用于后續(xù)模擬分析是可行的。

2.2 氣泡和水合物顆粒粒徑分布

圖3為不同工況下管道出口截面處氣泡和水合物顆粒粒徑徑向分布規(guī)律，粒徑大小用Sauter平均直徑表示。工況1、2、6為氣相和水合物顆粒體積分?jǐn)?shù)都為10%時(shí)，不同入口流速下（1.0、1.5、2.0 m/s）出口截面處氣泡和水合物顆粒大小徑向分布規(guī)律;工況3、4、5為入口流速2.0 m/s時(shí)，不同水合物顆粒體積分?jǐn)?shù)下（10%、20%、30%）出口截面處氣泡和水合物顆粒大小徑向分布規(guī)律;工況3、6、7為入口流速2.0 m/s時(shí)，不同氣相體積分?jǐn)?shù)下（5%、10%、15%）出口截面處氣泡和水合物顆粒大小徑向分布規(guī)律。

由圖3a可以發(fā)現(xiàn)，在近壁面處不同工況下氣泡粒徑較小并且變化不大，而管道中心氣泡粒徑在不同工況下變化較大。隨著入口速度增加，氣泡粒徑呈現(xiàn)出減小的趨勢(shì)，這與許多研究者的結(jié)論相同。隨著水合物顆粒體積分?jǐn)?shù)的增大，氣泡的粒徑隨之明顯增大，并且大氣泡位于管道中部區(qū)域，小氣泡位于液膜區(qū)域，這反映出水合物顆粒的存在對(duì)氣泡的聚并、破碎行為影響很大。隨著氣相體積分?jǐn)?shù)的增加，氣泡粒徑的變化不大，說明氣相濃度的變化對(duì)氣泡大小的影響較小。

從圖3b可以看出，各工況下水合物顆粒在近壁面處粒徑較小，在整個(gè)出口截面上呈現(xiàn)雙峰趨勢(shì)。隨著入口流速和氣相體積分?jǐn)?shù)的增加，水合物顆粒粒徑變小，但變化范圍不大。隨著入口水合物顆粒體積分?jǐn)?shù)的增加，水合物顆粒平均粒徑明顯增大，可見水合物顆粒體積分?jǐn)?shù)對(duì)水合物顆粒粒徑的影響較大。這反映了水合物顆粒的聚并、破碎速率與流動(dòng)剪切密切相關(guān)，在流動(dòng)剪切作用較強(qiáng)的區(qū)域和情況下，水合物顆粒的破碎占主導(dǎo)地位。

2.3 氣相和水合物顆粒濃度分布

圖4為不同工況下管道軸向截面氣相和水合物顆粒濃度分布云圖，截面為靠近出口部分?？梢园l(fā)現(xiàn)，所有不同工況下水合物顆粒和氣相都不均勻地分布在截面上。上述結(jié)果與YAO S P等[6]模擬水合物顆粒濃度分布的結(jié)果有所不同，主要存在以下兩個(gè)方面的原因：一是未考慮除曳力外其他相間作用力對(duì)分散相的影響，特別是升力和湍流分散力;二是氣相的存在增加了擾動(dòng)，使水合物顆粒濃度分布大不相同。

通過對(duì)各工況下濃度分布云圖的分析可知，隨著入口流速的增加，除管道中心區(qū)域氣相更加集中外，其他區(qū)域的氣相分布更加均勻，水合物顆粒靠近壁面。隨著水合物顆粒體積分?jǐn)?shù)的增加，氣相更加向著中心區(qū)域聚集，并與靠近壁面區(qū)域氣相濃度差別增大，分界線也更加清晰。隨著氣相體積分?jǐn)?shù)的提高，氣相濃度分布不均勻性增強(qiáng)。同時(shí)，可以發(fā)現(xiàn)，水合物顆粒濃度分布與氣相濃度分布區(qū)域和形態(tài)正好相反，這反映了水合物顆粒與液相的跟隨性較好。

2.4 壓力損失梯度

圖5分別展示了入口流速、水合物顆粒體積分?jǐn)?shù)、氣相體積分?jǐn)?shù)對(duì)壓力損失梯度的影響。由圖5a可知，管道內(nèi)壓力損失梯度隨入口流速的增加先降低后升高，這與徐海良等[25]所得結(jié)果是一致的，存在一個(gè)壓力損失梯度最低值的最優(yōu)流速區(qū)間。由圖5b可知，管內(nèi)壓力損失梯度隨水合物顆粒體積分?jǐn)?shù)增加而變大。這反映了水合物顆粒體積分?jǐn)?shù)增加導(dǎo)致混合密度變大，垂直方向上靜壓損失增大;同時(shí)，隨水合物顆粒體積分?jǐn)?shù)增加，水合物顆粒間碰撞概率、水合物顆粒與管壁間碰撞概率變大，能量損失增加。由圖5c可知，隨著氣相體積分?jǐn)?shù)的增加，壓力損失梯度降低。這是由于氣相體積分?jǐn)?shù)的增加，垂直管中靜壓損失和混合粘度降低，因此，整個(gè)流動(dòng)過程的能量損失就減小了。

3 結(jié)論

3.1 在忽略水合物生長(zhǎng)和分解、氣相溶解和析出的情況下，所建立模型能夠準(zhǔn)確可靠地進(jìn)行氣液水合物三相流動(dòng)的模擬研究。

3.2 氣泡平均直徑沿管道徑向呈單峰或者雙峰分布，近壁面處氣泡粒徑一直較小，氣泡平均直徑隨入口混合流速的增加而減小，隨著水合物顆粒體積分?jǐn)?shù)的增大而明顯增加，隨著氣相體積分?jǐn)?shù)的增加，氣泡粒徑的變化不明顯。水合物顆粒平均直徑在整個(gè)截面上呈現(xiàn)雙峰分布，隨著入口流速和氣相體積分?jǐn)?shù)的增加，水合物顆粒粒徑變小，隨著入口水合物顆粒體積分?jǐn)?shù)的增加，水合物顆粒粒徑明顯增大。

3.3 氣相和水合物顆粒不均勻地分布在軸向截面上，氣相隨著入口流速和水合物顆粒濃度增加向著管道中心區(qū)域聚集，隨著氣相體積分?jǐn)?shù)的提高，分布不均勻性增強(qiáng)，水合物顆粒濃度分布與氣相濃度分布的區(qū)域和形態(tài)相反。其中水合物顆粒體積分?jǐn)?shù)對(duì)濃度分布的影響最大。

3.4 豎直管道內(nèi)混合流體的壓力損失梯度隨入口流速的增加先降低后升高，存在特定工況下的最優(yōu)流速，隨著水合物顆粒體積分?jǐn)?shù)增加壓力損失梯度增大，隨著氣相體積分?jǐn)?shù)的增加壓力損失梯度明顯降低，這反映出了氣泡的減阻效應(yīng)。

參 考 文 獻(xiàn)

[1]" "ZHANG Z G，WANG Y，GAO L F，et al.Marine gas hydrates： future energy or environmental killer？[J].Energy Procedia，2012，16：933-938.

[2]" "祝有海，龐守吉，王平康，等.中國(guó)天然氣水合物資源潛力及試開采進(jìn)展[J].沉積與特提斯地質(zhì)，2021，41（4）：524-535.

[3]" "李長(zhǎng)俊，黃婷，賈文龍.深水天然氣水合物及其管道輸送技術(shù)[J].科學(xué)通報(bào)，2016，61（22）：2449-2462.

[4]" "LI L，XU H L，YANG F Q.Three?phase flow of submarine gas hydrate pipe transport[J].Journal of Central South University，2015，22：3650-3656.

[5]" "劉艷軍，唐孝蓉，胡坤.天然氣水合物漿體分解對(duì)其在垂直管中流動(dòng)特性影響的研究[J].化學(xué)通報(bào)，2018，81（3）：267-273.

[6]" "YAO S P，LI Y X，WANG W C，et al.Numerical simulation of hydrate slurry flow characteristics in vertical pipes based on population balance theory[J].International Journal of Oil Gas and Coal Technology，2020，25（3）：319-339.

[7]" "CHEN W，XU H L，KONG W Y，et al.Study on three?phase flow characteristics of natural gas hydrate pipeline transmission[J].Ocean Engineering，2020，214：107727.1-107727.17.

[8]" CAO X W，YANG K R，LI T T，et al.Numerical simulation of hydrate particle behaviors in gas?liquid flow for horizontal and inclined pipeline[J].Case Studies in Thermal Engineering，2021，27：101294.1-101294.13.

[9]" "TOMIYAMA A，CELATA G P，HOSOKAWA S，et al.Terminal velocity of single bubbles in surface tension force dominant regime[J].International Journal of Multiphase Flow，2002，28（9）：1497-1519.

[10]" GIDASPOW D.Multiphase Flow and Fluidization：Continuum and Kinetic Theory Descriptions[M].New York：Academic Press，1994：35-37.

[11]" "WANG T F，WANG J F，JIN Y.Experimental Study and CFD Simulation of Hydrodynamic Behaviours in an External Loop Airlift Slurry Reactor[J].The Canadian Journal of Chemical Engineering，2004，82（6）：1183-1190.

[12]" "FRANK T.Advances in computational fluid dynamics（CFD）of 3?dimensional gas?liquid multiphase flows[C]//Simulation of Complex Flows（CFD）?Applica?tions and Trands.Wiesbaden：NAFEMS Seminar，2005：1-18.

[13]" "LOPEZ D B M.Turbulent bubbly two?phase flow in a triangular duct[D].Troy：Rensselaer Polytechnic Inst?itution，1992.

[14]" "BURNS A D，F(xiàn)RANK T，HAMILL I，et al.The Favre averaged drag model for turbulent dispersion in Eulerian multi?phase flows[C]//5th International Conference on Multiphase Flow.Janpan：ICMF，2004：1-17.

[15] RAMKRISHNA D.Population Balances：Theory and Applications to Particulate Systems in Engineering[M].San Diego：Academic Press，2000.

[16]" "ABRAHAMSON J.Collision rates of small particles in a vigorously turbulent fluid[J].Chemical Engineering Science，1975，30（11）：1371-1379.

[17]" "SAFFMAN P G，TURNER J S.On the collision of drops in turbulent clouds[J].Journal of Fluid Mechanics，1956，1（1）：16-30.

[18]" "CAMP T R，STEIN P C.Velocity gradients and inte?

rnal work in fluid motion[J].Journal of the Boston Society of Civil Engineers，1943，30：219-230.

[19]" "VEN T，MASON S G.The microrheology of colloidal dispersions VII. Orthokinetic doublet formation of spheres[J].Colloid and Polymer Science，1977，255（5）：468-479.

[20]" "SONG G C，LI Y X，WANG W C，et al.Numerical simulation of hydrate slurry flow behavior in oil?water systems based on hydrate agglomeration modelling[J].Journal of Petroleum Science and Engineering，2018，169：393-404.

[21] MARCHISIO D L，VIGIL R D，F(xiàn)OX R O.Implementation of the quadrature method of moments in CFD codes for aggregation?breakage problems[J].Chemical Engineering Science，2003，58（15）：3337-3351.

[22] LEHR F，MILLIES M，MEWES D.Bubble?Size distributions and flow fields in bubble columns[J].AIChE Journal，2002，48（11）：2426-2443.

[23]" BALAKIN B V，PEDERSEN H，KILINC Z，et al.Turbulent flow of freon R11 hydrate slurry[J].Journal of Petroleum Science and Engineering，2010，70（3?4）：177-182.

[24] WANG T F，WANG J F，YANG W G，et al.Experimental study on bubble behavior in gas?liquid?solid three?phase circulating fluidized beds[J].Powder Technology，2003，137（1?2）：83-90.

[25]" "徐海良，孫思聰，楊放瓊.天然氣水合物三相流段管道壓力損失分析[J].斷塊油氣田，2021，28（5）：661-666.

（收稿日期：2022-04-14，修回日期：2023-01-13）

Investigation of Gas?Liquid?Hydrate Flow Characteristics

in Vertical Pipe Based on OpenFOAM

ZHOU Hong?li， WANG Chun?sheng

（School of Petroleum Engineering， Northeast Petroleum University）

Abstract" "Aiming at the complicated gas?liquid?solid flow involved in the exploitation of natural gas hydrate， a new combination of hydrate particle coalescation?breakage model was proposed based on the population balance model and the coalescation?breakage behavior of bubbles was considered. The model was solved and verified with OpenFOAM in the CFD software， including the investigation of the effects of inlet velocity， gas phase and hydrate particle volume fraction on gas?liquid?hydrate flow characteristics in vertical tubes. The results show， the established model is accurate and effective， and the concentration distribution of each phase， the radial distribution of the average particle size of bubbles and the average particle size of hydrate particles are most affected by the volume fraction of hydrate particles. The pressure loss gradient decreases first and then increases with the increase of the flow rate， and there is an optimal flow rate under certain conditions. With the increase of gas volume fraction， the pressure loss gradient decreases obviously.

Key words" "natural gas hydrate， vertical pipe， multiphase flow， population balance model， OpenFOAM

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目（21076043）。

作者簡(jiǎn)介：周洪利（1988-），碩士研究生，從事復(fù)雜流體流動(dòng)與數(shù)值模擬研究，zhou_hli@163.com。

引用本文：周洪利，王春生.基于OpenFOAM的豎直管內(nèi)氣-液-水合物流動(dòng)特性研究[J].化工機(jī)械，2023，50（1）：67-75.