天然氣水合物分解動(dòng)力學(xué)模型研究

雷昊 米洪剛 馬勇新 戚志林 嚴(yán)文德

(1.中海石油(中國(guó))有限公司湛江分公司，廣東 湛江 524057;

(2.重慶科技學(xué)院復(fù)雜油氣田勘探開發(fā)重慶市重點(diǎn)實(shí)驗(yàn)室，重慶 401331)

天然氣水合物是21世紀(jì)的新能源，主要分布在深海沉積物或陸域永久凍土中，由天然氣客體分子與水分子籠形晶格在高壓低溫條件下形成的類冰狀的結(jié)晶物質(zhì)。分子籠形晶格使用比各自籠子直徑更小的弱極性客體分子來填充［1-3］。“本征動(dòng)力學(xué)”主要是對(duì)天然氣水合物分解動(dòng)力學(xué)系統(tǒng)的“反應(yīng)”進(jìn)行系統(tǒng)研究［4-5］，且參數(shù)變化由分解期間的攪拌速率、溫度和壓力共同影響。在超出一定攪拌速率的情況下對(duì)2種不同的動(dòng)力學(xué)機(jī)制之間的變化進(jìn)行觀察，即天然氣水合物的分解動(dòng)力學(xué)。近些年，研究者對(duì)實(shí)驗(yàn)過程提出了改進(jìn)，采用改進(jìn)的模型對(duì)數(shù)據(jù)進(jìn)行驗(yàn)證［6-10］。針對(duì)以上問題，研究者進(jìn)行了分子動(dòng)力學(xué)(MD)的甲烷水合物溶解過程模擬。結(jié)果發(fā)現(xiàn)，水合物的分解是一個(gè)非平衡的過程，故無法使用平衡熱力學(xué)進(jìn)行描述。

1 天然氣水合物分解動(dòng)力學(xué)模型

1.1 氣體水合物模型

模型中涉及變量包括:氣體水合物平衡壓力p*，分解吸收的焓熱量Δhβ-LOW，水和客體分子化學(xué)勢(shì)μ［11-12］。在液態(tài)水存在的狀態(tài)下，平衡條件如下:

文中，各變量的上標(biāo)H表示水合物，上標(biāo)L表示液態(tài)水，下標(biāo)G表示客體分子，下標(biāo)W表示水分子，下標(biāo)O表示純物質(zhì)，上標(biāo)β表示空水合物晶格，空水合物晶格的水化學(xué)位之差為ΔμβOW-L，客體分子晶體的水化學(xué)位之差為

通過以下方程組計(jì)算水在不同結(jié)構(gòu)中變相的熱質(zhì)特性:

式中:h— 焓熱量，J/mol;

cp— 等壓熱容，MPa-1;

δ—與溫度相關(guān)的線性導(dǎo)數(shù);

v— 容度，m3;

k—空水合物晶格的等溫可壓縮性，s-1。

T— 溫度，K。

關(guān)于水合物相的化學(xué)勢(shì)客體分子的貢獻(xiàn)ΔμH-βW定義為:

式中:vi—不同孔穴類型的容度;θi—分子i的籠形空間占有率。

使用朗格繆爾系數(shù)Ci以及在平衡態(tài)下相應(yīng)的氣相逸度fG計(jì)算籠形空間占有率:

相對(duì)于徑向坐標(biāo)r，有式(8):

客體參數(shù)εi和籠格參數(shù)r0i的物理意義是中央井電位的飽和深度以及零交叉點(diǎn)的徑向位置。根據(jù)依賴?yán)矢窨姞栂禂?shù)理想的范特霍夫行為，確定了一個(gè)精確溫度系數(shù)Ci:

式中Ci為水合物孔穴的朗格繆爾系數(shù)，MPa-1。

對(duì)溫度的朗格繆爾系數(shù)的對(duì)數(shù)求導(dǎo)來計(jì)算客體分子的解吸焓熱值，利用CWP，獲得式(11):

1.2 水合物相的幾何特性

應(yīng)用統(tǒng)計(jì)速率理論(statistial rate theory of interfacial transport，簡(jiǎn)稱SRT)時(shí)需要用到客體分子的平均接觸面積σG，其估計(jì)如下:

σi為第i類型孔穴內(nèi)部的客體分子平均接觸面積，根據(jù)分子i的籠形空間占有率θi及格子的化學(xué)計(jì)算系數(shù)vi，得到式(14):

利用氣水合物單元晶體內(nèi)部水分子的數(shù)值~NUCW以及單元分子和格子體積vUCW(單元尺寸aUC孔穴直徑 ri)，得到 σW的表達(dá)式［14］。

1.3 狀態(tài)方程式

根據(jù)使用馬蒂亞斯和科普曼修正的能量參數(shù)的Redlich-Kwong-Soave狀態(tài)方程，對(duì)yW是在假定理想的Raoult-Dalton行為下進(jìn)行直接估計(jì):

1.4 溶解度

假設(shè)vLG，∞為無限稀釋的偏摩爾體積，水分子的客體溶解度xG的表達(dá)式如下:

要求使用狀態(tài)方程計(jì)算氣相逸度fvG，亨利系數(shù)HLG及參數(shù) A，B，C:

由于連續(xù)解析和溶解模型(consecutive desorption and melting model，簡(jiǎn)稱CDM)在當(dāng)前公式中是二元系統(tǒng)(水+1個(gè)客體分子)，以下的表述是有效的:

1.5 擴(kuò)散系數(shù)

為了得到依賴溫度的擴(kuò)散系數(shù)DL的正確描述，使用了來自于實(shí)驗(yàn)數(shù)據(jù)參數(shù)的阿列紐斯相關(guān)性。

2 改進(jìn)的天然氣水合物分解動(dòng)力學(xué)模型

2.1 模型基本思想和邊界條件

CDM的基本思想是假設(shè)天然氣水合物的分解局部連續(xù)，且是重復(fù)“反應(yīng)”的分解過程，即:(1)一部分水籠客體分子在固體—流體界面解析;(2)空水籠是局部溶解。在固體—流體界面暴露出一個(gè)新的客體分子，并且過程將持續(xù)直到不穩(wěn)定因素不存在或水合物完全分解。

模型理論推導(dǎo)的邊界條件定義如下:(1)單個(gè)客體水合物(SI和SII);(2)沒有外部熱傳質(zhì)的限制;(3)等溫分解過程。

CDM模型內(nèi)采用上述分解過程的基本假設(shè)來估計(jì)分解動(dòng)力學(xué)，同時(shí)使用碰撞頻率模型描述和溶解動(dòng)力學(xué)解吸來估計(jì)2個(gè)過程的動(dòng)力學(xué)常數(shù)，對(duì)界面?zhèn)鬟f解吸過程的統(tǒng)計(jì)速率理論(SRT)［15］進(jìn)行描述。

2.2 改進(jìn)的模型

公式(22)闡述了固體水合物變?yōu)橐簯B(tài)水和水分子客體的氣水合物本征動(dòng)力學(xué)，公式(23)闡述了固體水合物變?yōu)橐簯B(tài)水以及在欠飽和環(huán)境溶解中的客體分子氣水合物本征動(dòng)力學(xué)。

分解(DC)方程

溶解(DL)方程為:

式中:σW—水的平均接觸面積，m2;

σG— 客體分子平均接觸面積，m2;

n˙ — 摩爾通量，mol/(s·m2);

DLW來估計(jì);

DLG—水相內(nèi)部客體分子的自擴(kuò)散系數(shù)，m2/s;

ρW— 液體密度，mol/m3;

R— 通用氣體常數(shù)，8.314 J/(mol·K);

x—液相質(zhì)量分?jǐn)?shù)，%，下標(biāo)G、W分別代表客體分子和水分子;

y—?dú)庀噘|(zhì)量分?jǐn)?shù)，%，下標(biāo)G、W分別代表客體分子和水分子;

ΔhβOW-L— 空晶體格焓熱量變化，J/mol;

θG—?dú)怏w分子的籠形空間占有率，%;

fG— 客體分子的液相逸度，Pa。

上式中，上標(biāo)“*”表示對(duì)應(yīng)的等溫三相(水合物-液體-氣體)平衡條件。

在條件限制下，˙NG為本征動(dòng)力學(xué)的客體分子分解的總速率，可以利用整體幾何水合物-流體界面AH-F乘法計(jì)算［16］:

式(22)、(23)、(24)是關(guān)于動(dòng)力學(xué)理論問題的一般描述。

3 CDM模型討論

3.1 CDM模型的邊界條件

在開始討論CDM結(jié)果之前，先定義幾個(gè)邊界條件:(1)假設(shè)流體的壓力和相應(yīng)的組合物在分解過程中連續(xù);(2)液體中具有游離流體;(3)水合物-流體邊界層無傳熱和傳質(zhì)限制的理想攪拌釜式反應(yīng)器特性。由于氣水合物分解包含整個(gè)宏觀驅(qū)動(dòng)在氣水合物-液態(tài)界面局部穿過分解邊界層，進(jìn)而產(chǎn)生了一個(gè)等溫分解過程。

3.2 時(shí)間常數(shù)

通過特征時(shí)間常數(shù)給出了客體分子分解的單個(gè)速率常數(shù)，其值取自式(22)和(23)的反應(yīng)常數(shù)和的倒數(shù)值:

式中，τ是水分子從水合物到液態(tài)的時(shí)間改變測(cè)量值，通過觀察得出利用解吸過程進(jìn)入液態(tài)比各自解吸進(jìn)入液態(tài)慢3個(gè)數(shù)量級(jí);同樣，空晶體格的溶解發(fā)生比解吸過程快幾個(gè)數(shù)量級(jí)。

根據(jù)以上觀點(diǎn)，認(rèn)為水合物溶解過程比分解的過程更慢，即客體分子的解吸過程具有速率限制。

3.3 分子驅(qū)動(dòng)勢(shì)能ΔΘ

基于依賴溫度的本征動(dòng)力學(xué)分解，根據(jù)客體分子集聚逸度和平衡逸度fG、f*G來定義客體分子驅(qū)動(dòng)勢(shì)能ΔΘ理論。

不同溫度的驅(qū)動(dòng)勢(shì)函數(shù)定義依賴于具體的氣體釋放速率。首先，隨著溫度的升高分解率出現(xiàn)升高趨勢(shì);其次，在溶解的過程中，線性趨勢(shì)很明顯，且隨著驅(qū)動(dòng)的增加，氣體釋放速率將會(huì)變慢。

3.4 分子反應(yīng)速率常數(shù)kΘD

進(jìn)一步討論，分子反應(yīng)速率常數(shù)kΘD用式(27)來定義:

使用客體分子在分解過程中的最大可能分解速率極限的測(cè)量值來解釋速率常數(shù)，進(jìn)一步討論依賴于溫度的分解反應(yīng)。

式中，EA表示分解過程中活化能，ΔhH-βG是客體分子從滿晶格籠到液相的分解焓熱量值。

根據(jù)熱力學(xué)觀點(diǎn)，活化能是關(guān)于系統(tǒng)跨越化學(xué)勢(shì)阻礙ΔμA-S＞0的能量需求。自發(fā)反應(yīng)進(jìn)行的條件要求為ΔμF-A＜0。對(duì)于ΔμA-S＞0的熱平衡叫做活化能，這里的熱平衡是指反應(yīng)熱。

3.5 勢(shì)函數(shù)

在實(shí)際應(yīng)用中，通過直接引用驅(qū)動(dòng)函數(shù)的手段來闡述實(shí)驗(yàn)性氣水合物分解動(dòng)力學(xué)，公式(31)闡述了2個(gè)隨機(jī)的勢(shì)函數(shù)ψ1和ψ2之間的關(guān)系:

式中，微分項(xiàng)是關(guān)于ψ2— ψ1平面的 HLV平衡曲線的斜率。為了證明在動(dòng)力學(xué)相變應(yīng)用過程中，本征動(dòng)力學(xué)是否起了重要作用，在實(shí)驗(yàn)過程中，它是通過全程觀察記錄固體被液體層潤(rùn)濕而成為游離氣體水合物的過程。因此，物理過程在水合物-流體界面發(fā)生更類似于客體分子直接水合物-液態(tài)解吸的情況。

為了獲得更接近本征動(dòng)力學(xué)的實(shí)驗(yàn)數(shù)據(jù)，長(zhǎng)度d采用液態(tài)水相DWG附近的內(nèi)部分解速率常數(shù)KCG，以及客體分子擴(kuò)散系數(shù)來定義:

在水合物-液態(tài)界面下，這個(gè)值等于本征動(dòng)力學(xué)在限速情況下擴(kuò)散邊界層的厚度，d相當(dāng)于本征動(dòng)力學(xué)在限速情況下晶體直徑。

4 CDM預(yù)測(cè)與MD模擬對(duì)比實(shí)例

對(duì)CDM預(yù)測(cè)與MD模擬結(jié)果進(jìn)行對(duì)比，在290 K時(shí)，應(yīng)用公式(23)與公式(24)進(jìn)行CDM的預(yù)測(cè)。

在圖1右邊區(qū)域，當(dāng)溶解客體的xG在接近飽和值290 K和15.96 MPa，即HLV的等溫平衡點(diǎn)時(shí)，CDM模型的分解速率急劇降低。從圖中可以看出，這是一種客體吸附控制的轉(zhuǎn)變方法:(1)超過了轉(zhuǎn)變范圍;(2)水合物晶格熔化方法;(3)客體分子溶化隨摩爾分?jǐn)?shù)遞減。

從結(jié)果看，在客體溶解的低摩爾分?jǐn)?shù)下，CDM與MD模型的分解速率有很好的一致性，且2個(gè)模型的邊界條件相似。

圖1 CDM預(yù)測(cè)與MD模擬結(jié)果對(duì)比

5 結(jié)語

本次研究應(yīng)用本征氣體水合物分解動(dòng)力學(xué)改進(jìn)模型對(duì)天然氣水合物成分的分解過程進(jìn)行預(yù)測(cè)，結(jié)論如下:

(1)在相同驅(qū)動(dòng)梯度條件下，天然氣水合物的分解速度比各成分單獨(dú)分解的速度相對(duì)更慢。

(2)用客體解吸驅(qū)動(dòng)梯度和氣水合物分解速率的上極限表示的水合物晶格熔解過程具有獨(dú)立性，應(yīng)用分子分解速率流標(biāo)準(zhǔn)化公式能得到不同水合物形式的分解速率常數(shù)。

(3)對(duì)客體分子的低濃度分解速率進(jìn)行模型預(yù)測(cè)與MD模擬，結(jié)果顯示出2個(gè)模型良好的一致性。

［1］Windmeier C，Oellrich L R.Theoretical Study of Gas Hydrate Decomposition Kinetics-Model Development［J］.J.Phys.Chem.，2013，117:10151-10161.

［2］VysniauskusA，Bishnoi P R.A Kinetic Study of Methane Hydrate Formation［J］.Chem.Eng.Sci.，1983，38(7):1061-1072.

［3］GuptaA.Methane Hydrate Dissociation Measurements and Modelling:the Role of Heat Transfer and Reaction Kinetics［G］.Colorado School of Mines，2007:101-108.

［4］Kim H C.A Kinetic Study of Gas Hydrate Decomposition［D］.Canada:University of Calgary，1985:1-25.

［5］Kim H C，Bishnoi P R，Heidemann R A，et al.Kinetics of Methane Hydrate Decomposition［J］.Chem.Eng.Sci.，1987，42(7):1645-1653.

［6］Clarke M.Determination of the Intrinsic Kinetics of Gas Hydrate Decomposition by Particle-Size Analysis［D］.Canada:University of Calgary，2000:1-50.

［7］Clarke M，BishnoiP R.Determination of the Intrinsic Rate of Ethane Gas Hydrate Decomposition［J］.Chem.Eng.Sci.，2000，55(21):4869-4883.

［8］Clarke M，Bishnoi P R.Determination of the Activation Energy and Intrinsic Rate Constant of Methane Gas Hydrate Decomposition［J］.Can.J.Chem.Eng.，2001，79(1):143-147.

［9］Clarke M，Bishnoi P R.Measuring and Modeling the Rate of Decomposition of Gas Hydrates Formed from Mixtures of Methane and Ethane［J］.Chem.Eng.Sci.，2001，56(16):4715-4724.

［10］Clarke M，Bishnoi P R.Determination of the Intrinsic Rate Constant and Activation Energy of CO2Gas Hydrate Decomposition Using In-Situ Particle Size Analysis［J］.Chem.Eng.Sci.，2004，59(14):2983-2993.

［11］Sloan E D，KohC A.Clathrate Hydrates of Natural Gases［M］.3rd.CRC Press:Boca Raton，F(xiàn)L，2007:95-105.

［12］Makogon Y F，Dunlap W.Hydrates of Hydrocarbons［M］.Tulsa:Pennwell Publishers，1997:1-30.

［13］Holder G D，Zetts S P，Pradhan N.Phase Behaviour in Systems Containing Clathrate Hydrates［J］.Rev.Chem.Eng.，1988，5(1):1-70.

［14］BazantM Z，Trout B L.A Method to Extract Potentials from the Temperature Dependence of Langmuir Constants for Clathrate-Hydrates［J］.Physica A，2001，300(1/2):139-173.

［15］Ward C A，F(xiàn)indlay R D，Rizk M.Statistical Rate Theory of Interfacial Transport.I.Theoretical Development［J］.J.Chem.Phys.，1982，76(11):5599-5605.

［16］周錫堂，陶鮮花，龐重軍.天然氣水合物分解動(dòng)力學(xué)研究進(jìn)展［J］.天然氣化工，2006，31(5):70-73.