不同相體系下H2S水合物生長過程分子動力學模擬

羅澤東，李愛蓉，劉發(fā)平，陸 原

（1.西南石油大學 化學化工學院，四川 成都 610500；2.中海油（天津）油田化工有限公司，天津 300450）

目前我國已探明氣田中有一百多個高含硫氣田（H2S體積分數(shù)在2%～10%之間），其累計探明儲量已超過1012m3，其中普光氣田天然氣中H2S含量高達13%～18%[1]。高含硫天然氣比普通天然氣在較高溫度和較低壓強下更容易形成氣體水合物[2]。其原因是在天然氣中H2S氣體的存在會促進水合物的形成，且一旦形成水合物就很難去除，在生產(chǎn)、運輸和加工過程中造成管道堵塞、采氣管道腐蝕和安全隱患等問題。因此，研究H2S水合物的形成過程尤為重要。但是，水合物的成核與生長過程是一個微觀現(xiàn)象，很難通過實驗進行觀察。同時，H2S具有腐蝕性和劇毒性使得從實驗層面研究其水合物的形成過程更具有挑戰(zhàn)性。

目前，分子動力學模擬（MD模擬）可以從分子層面考察氣體水合物成核與生長的微觀規(guī)律[3-8]。近年來已有文獻報道借助MD模擬探究了純H2S水合物的成核與生長行為。Liang等[9]模擬了在230～250 K的溫度范圍內(nèi)H2S水合物的微觀成核過程，結(jié)果表明，相對于CH4水合物，H2S水合物的成核速率更大，其原因在于H2S氣體在水中具有更大的溶解度。Liang等[10]采用局部過熱方法，熔化一半H2S水合物晶體，進而隨機產(chǎn)生H2S氣相、H2S溶液相和H2S水合物固相的三相體系與H2S溶液相、H2S水合物固相的兩相體系；再采用局部過熱與局部過冷的溫度調(diào)控方式，分別研究兩個體系下H2S水合物的生長過程。模擬結(jié)果表明，三相體系和兩相體系下H2S水合物的生長速率遠遠大于CH4水合物，這主要歸因于H2S分子在水中的高溶解度和遷移率。但是，從初始模型與研究手段來看，這種局部過熱耦合過冷的方式明顯加速了H2S水合物的生長速率。同時，水合物晶體種子也會促進水合物籠對溶液中H2S的吸附。這可能掩蓋了不同體系下H2S水合物更為真實的生長過程，所揭示的H2S水合物生長的微觀機理并不完善。

為了進一步探究溫度對不同體系下H2S水合物生長過程的影響，本文通過MD模擬，分別在H2S氣相、溶液相和水合物相的三相體系，以及H2S溶液相和水合物相的兩相體系中，研究了H2S水合物在較寬溫度范圍內(nèi)的晶體生長過程。就生長厚度、擴散系數(shù)、局部密度等因素對H2S水合物生長速率的影響進行了深入分析。

1 模型建立與模擬方法

1.1 模型建立

1.1.1 氣液固三相模型

H2S水合物屬于Ⅰ型水合物，晶格常數(shù)為1.1877 nm。本文H2S水合物單晶胞的氧原子坐標取自X射線衍射實驗結(jié)果[11]，然后對單晶胞的3個方向擴胞形成2×2×10單元的超晶體，并將其放置于2.4 nm×2.4 nm×12.0 nm的正交模擬盒子（z方向長度為12.0 nm）。本文在Liang等[7]關(guān)于H2S水合物生長模型的基礎(chǔ)上，只保留少部分晶體種子，以此建立三相與兩相體系下H2S水合物的生長模型。其中，三相建模細節(jié)為：（1）固定z方向左邊（2.4 nm × 2.4 nm × 2.4 nm）三維區(qū)域中的水合物晶體，剩余部分在350 K和10 MPa，等溫等壓系綜（NPT系綜）中，熔化2.5 ns后，產(chǎn)生H2S氣相、H2S溶液相和H2S水合物固相三相體系；（2）繼續(xù)固定左邊的晶體，熔化區(qū)從350 K分別降溫至285 K、275 K、270 K、265 K、260 K和250 K，降溫速率為0.75 K/ps，使熔化區(qū)逐漸平衡至各生長溫度；（3）隨后固定熔化區(qū)，左邊晶體分別在各個溫度下恒溫100 ps，以釋放晶體區(qū)與熔化區(qū)之間的應(yīng)力。

形成的H2S氣相、溶液相、水合物固相的三相初始模型如圖1所示。本文所有的MD模擬程序采用LAMMPS[12]軟件進行。

圖1 H2S氣液固三相初始MD模型

1.1.2 液固兩相模型

為了使初始構(gòu)型中不包含H2S氣相，故隨機移除H2S水合物單晶胞中一半的H2S分子，然后對單晶胞的3個方向擴胞形成2 × 2 × 10單元的H2S水合物超晶體，此時超晶體的晶穴占有率為50%。并將其放置于2.4 nm × 2.4 nm × 12.0 nm的正交模擬盒子（z方向長度為12.0 nm）。熔化和平衡過程的模擬細節(jié)與三相模型相同，H2S溶液相與H2S水合物固相的兩相初始模型如圖2所示。

圖2 H2S液固兩相初始MD模型

1.2 模擬方法

水分子間的相互作用勢選用TIP4P/2005勢函數(shù)[13]，而H2S分子之間的相互作用采用Kamath[14]三點模型，如表1所示，這些參數(shù)能夠準確反映H2S水合物的相平衡，并且在以往的研究中得到了廣泛的應(yīng)用[4,5]。H2O分子與H2S分子的Lennard-Jones（LJ）交叉作用參數(shù)由Lorentz-Berthelot混合規(guī)則計算，長程庫倫作用采用pppm算法[15]求解，LJ交叉作用與長程庫倫作用的截斷半徑均選取1.0 nm。在整個模擬過程中應(yīng)用三維周期邊界條件，時間步長為2 fs。對于產(chǎn)出階段設(shè)定NPT系綜，在6個初始溫度和10 MPa下正式模擬100 ns，分別研究不同體系下H2S水合物的生長情況。在每個溫度下均進行3次不同初速度隨機數(shù)種子的模擬，以求取數(shù)據(jù)平均值。

表1 LJ勢函數(shù)參數(shù)及原子電荷分配

采用3種不同的措施來跟蹤H2S水合物的生長過程，具體為：（1）氧原子的徑向分布函數(shù)（RDF）與H2S分子的配位數(shù)；（2）H2S和H2O分子的均方位移曲線（MSD）；（3）H2S分子密度分布。

RDF是指距離任意參考粒子α中心r處出現(xiàn)β粒子的幾率，β粒子可以與α粒子相同也可以不同，可以定義為式(1)。對g(r)在[0,r0]進行積分，可計算得到參考粒子在r0范圍內(nèi)的臨近配位數(shù)。

式中，dN表示r→r+dr間的粒子數(shù)目；ρ為系統(tǒng)的密度，g/cm3。

均方位移如式(2)：

式中，MSD(t)為分子的均方位移，10-20m2；尖括號表示對總原子數(shù)的平均；ri(t)和ri(0)分別表示i粒子在t時刻的位置和在0時刻的位置，其中時間t的單位為ns。

對于固體分子來說，運動受到限制，其均方位移在0附近波動；液體分子的運動相對自由，其均方位移隨模擬時間增加而增大。根據(jù)愛因斯坦的擴散定律，有式(3)：

式中，D為分子擴散系數(shù)，10-11m2/s。

1.3 模型驗證

圖3(a)是溫度為260 K、壓力為10 MPa下氧原子的RDF。第一波峰在0.1～100 ns內(nèi)均位于0.273 nm，表明水合物晶胞中五邊形的邊長為0.273 nm，與文獻數(shù)據(jù)吻合[3,4]。圖3(b)為260 K、10 MPa下與H2S分子相關(guān)的臨近配位數(shù)。在距離中心H2S分子的0.572 nm范圍內(nèi)所有的H2O分子被定義為水臨近配位數(shù)，同時在距離中心H2S分子的0.912 nm范圍內(nèi)所有的H2S分子被定義為硫臨近配位數(shù)[6]。從圖3(b)可以看出，H2S分子的水臨近配位數(shù)達到20個，表明構(gòu)成512小籠中水分子個數(shù)為20個，與水合物晶體的幾何結(jié)構(gòu)保持一致。硫臨近配位數(shù)逐漸減小至13，H2S分子的第一臨近硫配位數(shù)與文獻[6]保持一致，表明離中心H2S分子的臨近H2S分子達到13個。由此可見，在高濃度H2S氣相區(qū)，H2S氣體被消耗用于水合物晶體的生長。

圖3 260 K、10 MPa下氧原子徑向分布函數(shù)(a)與H2S分子的配位數(shù)(b)演變

2 結(jié)果與討論

2.1 氣液固三相體系下H2S水合物的生長過程

2.1.1 構(gòu)型分析與生長速率

通過NPT模擬，在285 K、275 K、270 K、265 K、260 K和250 K的初始溫度下探究三相體系下H2S水合物晶體生長的微觀過程。圖4為該體系下H2S水合物100 ns時的生長快照。在100 ns的時間跨度內(nèi)，6個初始溫度在相平衡穩(wěn)定區(qū)以內(nèi)均發(fā)生晶體生長，并且從圖4發(fā)現(xiàn)在3個較低溫度下，最終構(gòu)型在界面處出現(xiàn)無定型的水籠結(jié)構(gòu)和帶有局部特征的H2S氣泡；而在圖4的3個較高溫度下，高溫加劇了H2S氣泡的聚集程度，最終體系具有明顯相分離現(xiàn)象，新晶相呈現(xiàn)長程有序分布。

圖4 H2S水合物在三相體系不同初溫下100 ns時的生長快照

為了更深入理解三相體系下水合物晶體生長過程中的水分子運動情況，本文在250 K下，隨機追蹤一個H2S分子在20～80 ns內(nèi)周圍水分子的微觀運動現(xiàn)象，如圖5所示。

圖5 某個H2S分子的周圍水結(jié)構(gòu)在20～80 ns的演變過程

由圖5可知，在第27 ns時，18個水分子通過氫鍵作用形成半籠型包圍這個H2S分子，36 ns后，周圍水分子形成比較完整的512籠型。在48 ns后，512籠繼續(xù)吸引周圍水分子，以此籠為基礎(chǔ)連接一個六邊形。對于水合物籠的拓撲結(jié)構(gòu)來說，水分子構(gòu)成的六邊形是形成大籠的初始條件。然而，在63 ns時六邊形消失，多個水分子形成更復(fù)雜的無定型團簇。即使這里的512小籠在10 ns間隔內(nèi)并未消失，但在遠離界面處的液相中，長達數(shù)十納秒并未出現(xiàn)完整的單晶胞水合物。這是因為雖然低溫可以提供較大的驅(qū)動力促使晶體成核，但三相體系下水分子運動受到限制，擴散較慢，水分子需要頻繁更新自身坐標才能到達成核位點，因此需要更長時間生長出單晶胞。該溫度下生成的無定型團簇是一個亞穩(wěn)態(tài)中間體，無定型團簇的形成將延遲晶體的進一步生長。這與Jacobson等[16,17]認為高驅(qū)動力下，非晶核（即無定型團簇）在動力學上可能比晶核更有利的結(jié)果保持一致，因為非晶核形成的自由能勢壘較低。

在三相體系下，溫度對H2S水合物生長速率的影響表現(xiàn)為其生長速率先增大后減小的趨勢，其結(jié)果如表2所示。250 K時水合物的生長厚度為3.6 nm，為最小值；260 K時水合物的生長厚度為6.4 nm，為最大值，此時水合物的最大生長速率為0.120 nm/ns，且在260～270 K下，H2S水合物的生長速率相差不大。文獻報道的MD模擬與實驗結(jié)果均表明水合物新晶體的最大生長速率會出現(xiàn)在某個中等溫度區(qū)間內(nèi)[18,19]，這說明在三相體系下溫度對氣體水合物生長過程的影響具有相似性。本文計算得到的H2S水合物生長速率（260 K，0.120 nm/ns）大于同體系下CH4水合物的生長速率（260 K，0.054 nm/ns）[20]。

表2 三相體系下H2S水合物于不同溫度的生長數(shù)據(jù)

2.1.2 擴散系數(shù)

為了考察H2S分子和H2O分子在三相體系下運動的快慢，分別計算了H2O和H2S分子在不同溫度下的均方位移，結(jié)果如圖6所示。由圖6可知，隨生長溫度升高，H2O與H2S分子的MSD值逐漸增大。與H2O分子相比，20 ns后H2S分子的MSD曲線出現(xiàn)較大波動，表明高溫加速H2S分子運動，使溶液相中的H2S分子朝向H2S氣相區(qū)運動，相分離狀態(tài)更加明顯，從而延緩晶體的進一步生長。

根據(jù)式(3)，通過擬合圖6中H2O與H2S分子的MSD曲線，得出各溫度下H2O與H2S分子的擴散系數(shù)，如表3所示。H2S分子的擴散系數(shù)約為H2O分子擴散系數(shù)的8～10倍，且H2O分子的擴散系數(shù)（275 K，2.472 ×10-10m2/s）小于文獻中純水的自擴散系數(shù)（275 K，11.49 × 10-10m2/s）[21]，表明晶體生長過程中水的運動受到限制，H2S水合物的生長速率主要受H2S分子擴散的影響。而文獻報道關(guān)于CH4水合物的三相生長過程中，發(fā)現(xiàn)CH4的擴散系數(shù)小于H2O的擴散系數(shù)[20,22]?？梢姡琀2S分子在水相中的擴散能力大于CH4分子，使得三相體系下H2S水合物的生長速率大于CH4水合物。因此氣體水合物的生長速率與氣體分子在水相中的擴散能力較強相關(guān)。

圖6 三相體系下H2O分子(a)與H2S分子(b)均方位移隨時間的變化

表3 三相體系下不同溫度H2O分子與H2S分子的擴散系數(shù)

2.1.3 H2S的密度分布

通過計算H2S分子在z方向的質(zhì)量密度分布可以觀察出新晶體在空間中排布的對稱性，如圖7所示。

圖7 三相體系下z方向H2S分子的密度分布

由圖7可知，在285 K下，新生長出的晶體中H2S密度分布呈現(xiàn)上下波動、周期性復(fù)現(xiàn)趨勢，與晶體種子區(qū)的密度變化趨勢保持一致。而250 K下，H2S分子的密度出現(xiàn)局部最大值，表明低溫下較低的擴散系數(shù)使得密度分布不均。初始溫度增加可以強化H2O分子與H2S分子的擴散，使得新晶相呈現(xiàn)長程有序分布，與圖4不同溫度下的水合物生長快照保持一致。

2.2 液固兩相體系下H2S水合物的生長過程

2.2.1 構(gòu)型分析與生長速率

從圖8可以看出，250 K下的最終構(gòu)型呈現(xiàn)無定型的水籠結(jié)構(gòu)（黃色方框區(qū)域）；其余5個溫度的最終構(gòu)型表現(xiàn)為新晶相長程有序，僅有少數(shù)H2S分子溶解在水中。與三相體系不同的是，在較高溫下生長構(gòu)型并未出現(xiàn)相分離狀態(tài)。

圖8 H2S水合物在兩相體系不同初溫下50 ns時的生長快照

表4為不同初始溫度下H2S水合物的兩相生長數(shù)據(jù)。250 K下水合物生長厚度為3.0 nm，為最小值；260 K下水合物生長厚度為6.9 nm，達到最大值。除250 K外，其余5個生長溫度對應(yīng)的生長速率相差不大，約為0.220 nm/ns。與三相體系相比，兩相體系下溫度對H2S水合物生長過程的影響不是很敏感，在初始溫度為260 K條件下，兩相體系下生長速率0.228 nm/ns大于三相水體系的生長速率0.120 nm/ns。

表4 兩相體系下H2S水合物于不同溫度的生長數(shù)據(jù)

為進一步觀察兩相體系的H2S水合物晶體生長，跟蹤了260 K下固液界面處晶體厚度的演變過程。圖9為兩相體系下0.5～12 ns內(nèi)H2S水合物晶體生長過程，第7 ns時，新生長出的厚度為1.6 nm，瞬時生長速率可達0.228 nm/ns，該體系下表現(xiàn)為界面處的逐層生長。表明兩相體系下的界面遷移行為可促進H2S水合物的生長，而三相體系下帶有局部特征的H2S氣泡卻延緩了晶體生長。

圖9 兩相體系下H2S水合物在固液界面處晶體厚度的演變

2.2.2 擴散系數(shù)

圖10為兩相體系下H2O與H2S分子均方位移曲線MSD隨時間的變化。

圖10 兩相體系下H2O分子(a)與H2S分子(b)均方位移隨時間的變化

圖10(a)中H2O分子的MSD隨著時間的變化迅速增加；而從圖10(b)可以看出，10 ns后，H2S分子的均方位移曲線變得較為平緩，這是因為兩相體系下H2S水合物的生長速率較大，溶液中的H2S分子在較短時間內(nèi)被水合物籠捕捉完成，最終體系變?yōu)镠2S水合物和純水兩相。

通過擬合圖10中H2O與H2S分子的MSD曲線，得出兩相體系下H2O和H2S分子的擴散系數(shù)，如表5所示。

表5 兩相體系下不同溫度H2O分子與H2S分子的擴散系數(shù)

兩相體系下，H2O分子和H2S分子的擴散系數(shù)十分接近。表明兩相體系中主體分子與客體分子的擴散能達到相互協(xié)同的效果，進而促進H2S水合物晶體生長，增加了H2S水合物的生長速率。

2.2.3 H2S的密度分布

圖11為兩相體系下100 ns時，z方向H2S分子的密度分布。與250 K相比，285 K下z方向H2S分子密度峰值較大且分布有序。這表明，通過增加初始溫度，可使得新生長出的水合物晶體變得長程有序，與三相體系下H2S分子的密度分布相似。

圖11 兩相體系下z方向H2S分子的密度分布

3 結(jié)論

本文采用MD模擬，研究了三相和兩相體系下溫度對H2S水合物生長過程的影響，并借助生長構(gòu)型、均方位移曲線和密度分布等手段對比討論了不同溫度條件下H2S水合物生長過程的差異，得到主要結(jié)論如下：

（1）三相體系下，溫度對H2S水合物生長過程的影響較大。無定型結(jié)構(gòu)出現(xiàn)在低溫區(qū)域，而具有長程有序的新晶體出現(xiàn)在高溫區(qū)域。250～285 K溫度范圍內(nèi)，H2S水合物的最大生長速率會出現(xiàn)在某個中等溫度區(qū)間內(nèi)。

（2）兩相體系下，溫度對H2S水合物生長過程的影響不敏感。除250 K，其余溫度下生長出的新晶體均表現(xiàn)為長程有序，且生長速率相差不大，約為0.220 nm/ns。

（3）H2S水合物在兩相體系下的生長速率大于三相體系。三相體系下H2S分子的擴散系數(shù)約為H2O分子擴散系數(shù)的8～10倍，較大的擴散系數(shù)使得體系產(chǎn)生相分離狀態(tài)，不利于晶體生長；而在兩相體系下H2O分子與H2S分子的擴散能達到相互協(xié)同的效果，促進晶體生長。

本研究結(jié)論可為含硫氣藏天然氣開采中水合物的防治提供理論參考。尤其是H2S分子溶解在水相的體系，可以考慮添加能夠減弱H2O與H2S分子間相互擴散能力的抑制劑，在水合物生長初期抑制其形成。