廖清云，史博會(huì)，宋尚飛，王君傲，姚海元，段 旭，宮 敬，聞 峰

（1. 中國(guó)石油大學(xué)（北京） 天然氣水合物國(guó)家重點(diǎn)實(shí)驗(yàn)室，油氣管道輸送安全國(guó)家工程實(shí)驗(yàn)室，石油工程教育部重點(diǎn)實(shí)驗(yàn)室，城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 102249；2. 中海油研究總院有限責(zé)任公司，北京 100028；3. 國(guó)家石油天然氣管網(wǎng)有限公司，北京 100007）

隨著深海油氣資源開(kāi)發(fā)的快速發(fā)展，氣體水合物的開(kāi)發(fā)利用及其在流動(dòng)保障中的防治研究逐漸深入，其中，氣體水合物生長(zhǎng)及分解過(guò)程的微觀機(jī)理成為了研究熱點(diǎn)?？紤]到常規(guī)實(shí)驗(yàn)手段在時(shí)間尺度和空間尺度上的限制，國(guó)內(nèi)外已有多位學(xué)者利用分子動(dòng)力學(xué)模擬進(jìn)行了微觀尺度的研究，包括水合物在純水體系中的生長(zhǎng)機(jī)理[1]、添加劑或雜質(zhì)（如無(wú)機(jī)鹽/水合物動(dòng)力學(xué)抑制劑/天然果膠/蠟分子等）的存在對(duì)水合物生長(zhǎng)機(jī)理的影響[2-5]、水合物在逐步升溫過(guò)程中的分解機(jī)理[6]、水合物在無(wú)機(jī)鹽中的分解機(jī)理[7]及不同溫度和壓力條件下抑制劑存在時(shí)水合物的分解機(jī)理[8]等。

在水合物生長(zhǎng)及分解過(guò)程的分子模擬研究中，均需在初始模型中添加根據(jù)實(shí)驗(yàn)數(shù)據(jù)[9]搭建的水合物晶胞結(jié)構(gòu)作為水合物生長(zhǎng)種子或水合物分解前體結(jié)構(gòu)。因水合物晶胞結(jié)構(gòu)中氫原子位置具有不確定性，在分子模擬中需根據(jù)Bernal-Fowler規(guī)則、體系構(gòu)型的勢(shì)能最小原則或氫鍵水分子朝向成比例原則等，進(jìn)一步確定氫原子坐標(biāo)。在構(gòu)建模擬模型后，需開(kāi)展能量最小化以降低體系能量并尋找體系勢(shì)能極小值以消除體系內(nèi)不合理接觸，隨后通過(guò)弛豫過(guò)程，即一段有限時(shí)長(zhǎng)的分子動(dòng)力學(xué)模擬，以進(jìn)一步調(diào)整能量最小化后的模型，使系統(tǒng)從初始配置逐漸演變達(dá)到平衡狀態(tài)，才能為分子模擬計(jì)算提供更為合理的初始構(gòu)型，該過(guò)程對(duì)獲得合理的分子模擬結(jié)果至關(guān)重要。

弛豫過(guò)程以分析體系的熱力學(xué)參數(shù)（如溫度、能量和密度等）及結(jié)構(gòu)參數(shù)（如序參數(shù)、均方位移和徑向分布函數(shù)等）為主[10]。然而，不同研究者選擇的弛豫系綜及弛豫時(shí)間各有不同：Tung等[1]以具有固定溫度的正則系綜（NVT）（T= 200 K）弛豫20 ps平衡氣-液-固體系；Li等[3]以NVT系綜（T= 276 K）弛豫400 ps平衡液-固體系；Xu等[4]以NVT系綜（T= 260 K）弛豫30 ps，及具有固定溫度及壓力的系綜（NPT）（T= 260 K，p= 15 MPa）弛豫30 ps平衡液-固體系；Choudhary等[6]以NPT系綜（p= 100 MPa）弛豫5 ns平衡水合物晶胞結(jié)構(gòu)；Xu等[7]以NVT系綜在模擬溫度弛豫50 ps平衡液-固體系；Kondori等[8]以NVT系綜（T= 273 K）弛豫50 ps平衡水合物晶胞結(jié)構(gòu)；張慶東等[11]首先以NVT系綜（T= 263.15 K、T= 273.15 K、T= 283.15 K）固定氧原子位置弛豫50 ps，隨后解除氧原子約束在相同條件下弛豫1 ns平衡水合物晶胞結(jié)構(gòu)?？傊?，NVT系綜在弛豫中被廣泛應(yīng)用，且弛豫時(shí)間通常設(shè)置在皮秒量級(jí)，但關(guān)于如何確定合理的弛豫方案的研究較少。

本文以甲烷（CH4）水合物晶胞籠型結(jié)構(gòu)的演變規(guī)律為分析參數(shù)，探討不同弛豫方案對(duì)CH4水合物晶胞構(gòu)型的影響，包括：弛豫系綜、弛豫時(shí)間、弛豫壓力與溫度等，目的是確定更為合理的弛豫方案以開(kāi)展CH4水合物的分子動(dòng)力學(xué)模擬，為水合物分子動(dòng)力學(xué)模擬研究提供理論基礎(chǔ)。

1 模擬模型和方法

1.1 CH4水合物晶胞模型

圖1 3 × 3 × 3 sI型CH4水合物晶胞示意Fig. 1 Schematic diagram of 3 × 3 × 3 sI methane hydrate unit cell

1.2 模擬設(shè)定

基于相關(guān)文獻(xiàn)研究[1,3,5,16]，對(duì)于分子間相互作用，確定選用TIP4P-Ew力場(chǎng)[17]描述H2O分子，OPLS-AA力場(chǎng)[18]描述CH4分子。各原子的電荷及L-J勢(shì)函數(shù)參數(shù)如表1所示，不同原子間的L-J勢(shì)函數(shù)參數(shù)值通過(guò)Lorentz-Berthelot混合規(guī)則計(jì)算。將H2O分子視為剛體，采用SHAKE算法[19]約束H2O分子的H-O鍵長(zhǎng)與H-O-H鍵角。范德華力采用L-J勢(shì)函數(shù)計(jì)算，截?cái)喟霃綖? nm。短程庫(kù)倫作用截?cái)喟霃綖?.9 nm，長(zhǎng)程庫(kù)倫相互作用通過(guò)pppm算法計(jì)算[20]，計(jì)算精度為10-4。模擬時(shí)間步長(zhǎng)為1 fs，每0.5 ps輸出一次熱力學(xué)參數(shù)與原子位置信息。利用Nose-Hoover恒溫器及恒壓器實(shí)現(xiàn)體系的溫度壓力控制，模擬體系在x、y和z方向上均采用了周期性邊界條件。

表1 L-J勢(shì)函數(shù)參數(shù)值與原子電荷值Table 1 L-J potential function parameter value and atomic charge value

1.3 CH4水合物晶胞構(gòu)型穩(wěn)定性判定方法

本文以一種改進(jìn)的籠型結(jié)構(gòu)識(shí)別方法判定CH4水合物晶胞構(gòu)型穩(wěn)定性，其具體流程如下：

（1）篩選成籠氧原子。當(dāng)相鄰氧原子間的間距小于0.35 nm，且相鄰的兩個(gè)氧原子及與其相連的一個(gè)氫原子間形成的角度小于30°，則認(rèn)為這兩個(gè)氧原子之間形成氫鍵[21]，當(dāng)一個(gè)氧原子上形成的氫鍵數(shù)大于等于3時(shí)則認(rèn)為其為成籠氧原子。

（2）基于氧原子篩選結(jié)果，進(jìn)行五元環(huán)的識(shí)別。將成籠氧原子視為無(wú)向圖的頂點(diǎn)，根據(jù)各氧原子之間的連接關(guān)系來(lái)判斷是否形成五元環(huán)。

（3）512籠結(jié)構(gòu)如圖2所示，綠色實(shí)線(xiàn)為512半籠的中心五元環(huán)，紅色實(shí)線(xiàn)為半籠中與中心五元環(huán)相連的五元環(huán)，藍(lán)色實(shí)線(xiàn)為完整的512半籠結(jié)構(gòu)?？紤]到其為對(duì)稱(chēng)結(jié)構(gòu)，可采取半籠拼接的方法進(jìn)行識(shí)別[22]。根據(jù)五元環(huán)的識(shí)別結(jié)果，判斷各五元環(huán)之間的連接關(guān)系，進(jìn)而確定以各五元環(huán)為中心的512半籠，再將得到的512半籠進(jìn)行拼接，最終得到512全籠。

圖2 512籠結(jié)構(gòu)示意Fig. 2 Schematic diagram of 512 cage structure

水合物晶體構(gòu)型內(nèi)氫原子的位置會(huì)影響氫鍵判定中的角度計(jì)算結(jié)果，并影響所識(shí)別的環(huán)籠結(jié)構(gòu)數(shù)目。因此，以該改進(jìn)的籠型識(shí)別方法分析弛豫過(guò)程中體系籠型結(jié)構(gòu)數(shù)目，若所識(shí)別的籠型結(jié)構(gòu)數(shù)目越接近理論值且波動(dòng)越小，則說(shuō)明該弛豫方案的設(shè)定越優(yōu)。對(duì)于3 × 3 × 3 sI型CH4水合物晶胞結(jié)構(gòu)而言，按理論分析應(yīng)有35個(gè)512籠，其連接關(guān)系如圖3所示，其中紅色小球?yàn)檠踉?，藍(lán)色實(shí)線(xiàn)為氫鍵。

本文的成本系數(shù)(Ck)是指該方案的生命周期成本(工程造價(jià))與所有方案成本總和之比。各方案的價(jià)值系數(shù)Vk = Fk / Ck，計(jì)算結(jié)果見(jiàn)表3?？梢?jiàn)在4種防腐方案中，D方案的價(jià)值系數(shù)最高，即采用功能性玻璃鋼復(fù)合材料作為火電廠(chǎng)煙囪防腐工程材料是最優(yōu)解決方案。

圖3 3 × 3 × 3 sⅠ型CH4水合物晶胞中512籠的正交示意(a)及透視示意 (b)Fig. 3 Orthogonal schematic diagram (a) and perspective schematic diagram (b) of 512 cages in 3 × 3 × 3 sⅠmethane hydrate unit cell

2 結(jié)果與討論

2.1 弛豫系綜及弛豫時(shí)間對(duì)CH4水合物晶胞構(gòu)型的影響

為研究弛豫系綜及弛豫時(shí)間對(duì)CH4水合物晶胞構(gòu)型的影響，設(shè)定弛豫方案1-1及弛豫方案1-2。弛豫時(shí)間均為500 ps；弛豫方案1-1設(shè)定為NVT系綜，弛豫方案1-2設(shè)定為NPT系綜；以水合物-液態(tài)水-甲烷氣三相平衡溫壓條件為基礎(chǔ)，10 MPa時(shí)的水合物-液態(tài)水-甲烷氣三相平衡溫度約為281 K[1]，故暫設(shè)定弛豫溫度為250 K，弛豫壓力為10 MPa。為避免瞬時(shí)結(jié)構(gòu)波動(dòng)對(duì)結(jié)果分析造成影響，將512籠識(shí)別數(shù)量在3 ps的時(shí)間窗口內(nèi)進(jìn)行平均后再進(jìn)行分析。

方案1-1（NVT，T= 250 K，下同）弛豫500 ps的512籠識(shí)別結(jié)果如圖4(a)所示。在500 ps弛豫過(guò)程中，512籠數(shù)量變化可分為3個(gè)階段，如表2所示。在58 ps內(nèi)快速達(dá)到近平衡狀態(tài)，在58 ps后隨著模擬時(shí)間的延長(zhǎng)逐漸趨于更穩(wěn)定的平衡狀態(tài)。此外，體系溫度在± 5%范圍內(nèi)波動(dòng)，如圖4(b)所示，表明體系溫度已達(dá)到平衡狀態(tài)。

表2 不同弛豫方案下512籠數(shù)量變化階段Table 2 Change stages of 512 cage number under different relaxation schemes

圖4 方案1-1弛豫500ps的512籠識(shí)別結(jié)果(a)及溫度變化(b)Fig. 4 512 cage recognition results (a) and temperature changes (b) of scheme 1-1 relaxation 500 ps

方案1-2（NPT，T= 250 K，p= 10 MPa，下同）弛豫500 ps的512籠識(shí)別結(jié)果如圖5(a)所示。在500 ps弛豫過(guò)程中，512籠數(shù)量在弛豫初期波動(dòng)較大，其變化可分為4個(gè)階段，如表2所示。弛豫在20 ps內(nèi)快速達(dá)到近平衡狀態(tài)，在20 ps后隨著模擬時(shí)間的延長(zhǎng)逐漸趨于512籠數(shù)量更高的近平衡狀態(tài)，在364 ps后達(dá)到更穩(wěn)定的平衡狀態(tài)。此外，體系溫度在± 5%范圍內(nèi)波動(dòng)，如圖5(b)所示，表明體系溫度已達(dá)到平衡狀態(tài)。

圖5 方案1-2弛豫500 ps的512籠識(shí)別結(jié)果(a)及溫度變化(b)Fig. 5 512 cage recognition results (a) and temperature changes (b) of scheme 1-2 relaxation 500 ps

由表2可知，500 ps弛豫結(jié)束后，兩種弛豫方案下體系中識(shí)別得到的512籠數(shù)量難以穩(wěn)定達(dá)到35個(gè)。這是由于H2O分子在模擬中不斷運(yùn)動(dòng)，其原子存在越過(guò)模擬盒邊界的情況，進(jìn)而導(dǎo)致邊界處H2O分子所形成的氫鍵數(shù)未達(dá)到成籠條件而被忽略。方案1-1和1-2弛豫后的體系結(jié)構(gòu)如圖6所示，綠色實(shí)線(xiàn)為氫鍵，模擬盒邊緣的各方向應(yīng)各有3個(gè)512籠，因本環(huán)籠識(shí)別代碼中的氫鍵篩選過(guò)程暫未考慮周期性邊界條件，致使方案1-1未識(shí)別出的512籠位于模擬盒的左下邊界及右中邊界，方案1-2未識(shí)別出的512籠位于模擬盒的右下邊界，故512籠數(shù)量最終呈現(xiàn)出在接近理論值附近波動(dòng)變化的現(xiàn)象。

圖6 方案1-1(a)及方案1-2(b) 弛豫后的體系結(jié)構(gòu)Fig. 6 512 cages of scheme 1-1 (a) and scheme 1-2 (b) after relaxation

對(duì)于弛豫系綜，對(duì)比方案1-1和方案1-2得到的512籠識(shí)別結(jié)果，如圖7所示。由圖7(a)可知，在方案1-2下，512籠數(shù)量能更快達(dá)到較高的數(shù)值，但其在364 ps內(nèi)的數(shù)值波動(dòng)較大。經(jīng)對(duì)比發(fā)現(xiàn)，在較短的模擬時(shí)長(zhǎng)內(nèi)，方案1-1得到的結(jié)果穩(wěn)定性更好。因此，可以認(rèn)為方案1-2（NPT系綜）下弛豫對(duì)籠型結(jié)構(gòu)的形成有一定的促進(jìn)作用，但在方案1-1（NVT系綜）下弛豫更利于結(jié)構(gòu)的穩(wěn)定。這一現(xiàn)象與Zhang等[23]的結(jié)論一致：在NPT系綜下水合物模擬的成核速率最快，而在NVT系綜下模擬得到的水合物結(jié)晶度更高。由圖7(b)可知，方案1-1和1-2弛豫過(guò)程中體系溫度變化相近，均在平衡值附近波動(dòng)。

圖7 方案1-1與方案1-2弛豫的512籠識(shí)別結(jié)果(a)及溫度變化(b)對(duì)比Fig. 7 Comparison of 512 cage recognition results (a) and temperature change (b) between scheme 1-1 and 1-2

對(duì)于弛豫時(shí)間，由圖4(a)與圖5(a)可知，弛豫時(shí)間越長(zhǎng)，越有利于獲得穩(wěn)定的CH4水合物晶胞構(gòu)型，但弛豫時(shí)間的增加會(huì)顯著延長(zhǎng)計(jì)算時(shí)間。適宜的弛豫時(shí)間，不僅可以避免因弛豫時(shí)間較短不能使體系達(dá)到穩(wěn)定平衡的缺點(diǎn)，還可以彌補(bǔ)弛豫時(shí)間較長(zhǎng)而導(dǎo)致計(jì)算時(shí)長(zhǎng)增加的缺陷。對(duì)于波動(dòng)較大的弛豫方案1-2，由于512籠數(shù)量在模擬初期變化明顯，需選取適當(dāng)較長(zhǎng)的弛豫時(shí)間以達(dá)到體系的平衡，考慮到其在100 ps后波動(dòng)幅度有所減小，且此時(shí)兩種弛豫方案內(nèi)識(shí)別得到的512籠數(shù)量均可達(dá)到32個(gè)，故選用100 ps作為本文后續(xù)分析研究的弛豫時(shí)間。

基于方案1-1及方案1-2，為進(jìn)一步研究弛豫溫度和壓力條件及系綜聯(lián)用對(duì)CH4水合物晶胞構(gòu)型的影響，設(shè)定多種弛豫方案，如表3所示。設(shè)定弛豫時(shí)間為100 ps，分別設(shè)置不同溫度和壓力以分析弛豫溫壓條件對(duì)弛豫過(guò)程的影響，并研究NVT及NPT系綜聯(lián)合應(yīng)用對(duì)弛豫過(guò)程的影響。

2.2 弛豫溫度對(duì)CH4水合物晶胞構(gòu)型的影響

表3中弛豫方案2-1至2-3，是在弛豫方案1-1基礎(chǔ)上，調(diào)整弛豫時(shí)間為100 ps，改變弛豫溫度設(shè)定的3個(gè)方案。方案2-1至2-3的弛豫溫度分別設(shè)置為200 K、250 K和300 K。3種不同弛豫溫度下的籠識(shí)別結(jié)果對(duì)比如圖8(a)所示。對(duì)于最低弛豫溫度200 K，在相同弛豫時(shí)間下，體系內(nèi)識(shí)別得到的512籠數(shù)量最高，且其波動(dòng)幅度最小。然而，對(duì)于最高弛豫溫度300 K而言，體系內(nèi)所識(shí)別的512籠數(shù)量波動(dòng)幅度最大。弛豫溫度為250 K時(shí)，512籠識(shí)別數(shù)量結(jié)果處于兩者之間。總體而言，在弛豫時(shí)間為75~100 ps時(shí)，弛豫方案2-1至2-3所識(shí)別的512籠數(shù)量逐漸接近。

圖8 方案2-1至2-3弛豫的512籠識(shí)別結(jié)果(a)及溫度變化(b)Fig. 8 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 2-1 to 2-3

表3 弛豫方案Table 3 Relaxation schemes

因水合物的生成條件是低溫高壓，在低溫條件下，水合物生成的驅(qū)動(dòng)力更強(qiáng)，故較低的溫度有利于水合物籠型結(jié)構(gòu)的穩(wěn)定，更利于籠型結(jié)構(gòu)識(shí)別分析。但隨著弛豫時(shí)間達(dá)到100 ps，不同弛豫溫度所達(dá)到的水合物晶胞構(gòu)型穩(wěn)定效果逐漸接近。此模擬結(jié)果與張慶東等[11]得到的結(jié)論一致，即弛豫溫度會(huì)影響CH4水合物的穩(wěn)定性，且弛豫溫度越低，其穩(wěn)定性越高；但對(duì)于已穩(wěn)定的水合物晶胞結(jié)構(gòu)，弛豫溫度的影響將不再明顯。因此，圖8(a)的模擬結(jié)果也再次佐證了本模擬體系內(nèi)CH4水合物晶胞構(gòu)型在弛豫100 ps時(shí)已基本達(dá)到穩(wěn)定。方案2-1至2-3弛豫過(guò)程中的溫度變化如圖8(b)所示，可知各體系的溫度基本在設(shè)定值附近波動(dòng)。

越低的弛豫溫度越有利于CH4水合物晶胞構(gòu)型的穩(wěn)定，而過(guò)高的弛豫溫度會(huì)減弱CH4水合物晶胞構(gòu)型的穩(wěn)定性。對(duì)于TIP4P-Ew水分子模型，冰點(diǎn)約為245 K[24]，若弛豫體系中含有液態(tài)水，需設(shè)置高于冰點(diǎn)的模擬溫度以避免冰的產(chǎn)生。對(duì)于本文所選參數(shù)，三相平衡溫度在10 MPa時(shí)為281 K，弛豫溫度需低于所選力場(chǎng)模型的相平衡溫度以避免水合物結(jié)構(gòu)失穩(wěn)。對(duì)于本文研究的體系，250 K的弛豫溫度較合理，為此，表3中的弛豫方案3-1、3-2、4-1、4-2和4-3均選取250 K作為弛豫溫度，以開(kāi)展后續(xù)研究。

2.3 弛豫壓力對(duì)CH4水合物晶胞構(gòu)型的影響

表3中弛豫方案3-1和3-2，是在弛豫方案1-2基礎(chǔ)上，調(diào)整弛豫時(shí)間為100 ps，改變弛豫壓力設(shè)定的兩個(gè)方案。方案3-1與3-2的弛豫壓力分別設(shè)置為10 MPa與30 MPa。方案3-1和方案3-2在100 ps內(nèi)的512籠識(shí)別結(jié)果對(duì)比如圖9(a)所示。數(shù)據(jù)表明，增大壓力后512籠數(shù)量的變化更加劇烈。在100 ps模擬末，30 MPa弛豫壓力下的512籠識(shí)別數(shù)量仍有顯著波動(dòng)。這是因?yàn)榉肿觿?dòng)力學(xué)模擬通過(guò)調(diào)整模擬盒尺寸并隨之縮放各原子位置以實(shí)現(xiàn)壓力的控制，在高壓下體系的體積調(diào)整較大，故其結(jié)構(gòu)的波動(dòng)變化也越劇烈，具體表現(xiàn)為在高壓下弛豫時(shí)，所識(shí)別的512籠數(shù)量穩(wěn)定性差。方案3-1和方案3-2下體系溫度的變化如圖9(b)所示，數(shù)據(jù)表明體系溫度基本在設(shè)定值附近波動(dòng)。因此，對(duì)于NPT系綜下弛豫壓力的確定，一方面，要避免高壓所引起的體系結(jié)構(gòu)變化；另一方面，應(yīng)滿(mǎn)足后續(xù)正式模擬壓力需求。為此，針對(duì)本文研究體系，為保證結(jié)構(gòu)穩(wěn)定，選取10 MPa作為NPT系綜弛豫壓力。

圖9 方案3-1和3-2弛豫的512籠識(shí)別結(jié)果(a)及溫度變化(b)Fig. 9 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 3-1 to 3-2

2.4 雙系綜弛豫方案對(duì)CH4水合物晶胞構(gòu)型的影響

經(jīng)前續(xù)模擬結(jié)果分析，在適宜的弛豫時(shí)間和弛豫溫度下，NVT系綜下的弛豫能得到穩(wěn)定的結(jié)構(gòu)，但水合物結(jié)構(gòu)達(dá)到一定數(shù)量所需的時(shí)間更長(zhǎng)；在適宜的弛豫溫度和壓力下，NPT系綜下可更快得到一定數(shù)量的水合物結(jié)構(gòu)，但水合物結(jié)構(gòu)的波動(dòng)變化較大；兩個(gè)系綜下的弛豫均可平衡體系溫度。因此，考慮將兩個(gè)系綜相結(jié)合，借助NVT系綜易于得到穩(wěn)定結(jié)構(gòu)的優(yōu)勢(shì)，加以NPT系綜對(duì)籠型結(jié)構(gòu)形成的促進(jìn)作用，得到表3中的弛豫方案4-1，總弛豫時(shí)長(zhǎng)為200 ps，首先在NVT下弛豫100 ps以穩(wěn)定初始結(jié)構(gòu)，再在NPT下弛豫100 ps以進(jìn)一步達(dá)到接近理論結(jié)構(gòu)的穩(wěn)定狀態(tài)。

方案4-1模擬過(guò)程中512籠識(shí)別數(shù)量的變化如圖10(a)所示。其亦可大致分為3個(gè)階段，如表4所示。弛豫在58 ps內(nèi)快速達(dá)到近平衡狀態(tài)，在58 ps后隨著模擬時(shí)間的延長(zhǎng)逐漸趨于更穩(wěn)定的平衡狀態(tài)。方案4-1弛豫過(guò)程中體系的溫度變化如圖10(b)所示，弛豫過(guò)程體系溫度基本在設(shè)定值附近波動(dòng)。因此，在弛豫方案4-1下，首先進(jìn)行NVT弛豫再進(jìn)行NPT弛豫可有效地平衡體系。

圖10 方案4-1弛豫的512籠識(shí)別結(jié)果(a)及溫度變化(b)Fig. 10 512 cage recognition results (a) and temperature changes (b) of relaxation scheme 4-1

表4 方案4-1 512籠數(shù)量變化階段Table 4 Change stages of 512 cage number of relaxation scheme 4-1

將方案4-1、4-2及4-3進(jìn)行對(duì)比，NVT弛豫后再進(jìn)行NPT弛豫，512籠數(shù)量在NPT弛豫初期會(huì)出現(xiàn)一定波動(dòng)，但在167 ps即達(dá)到33個(gè)512籠，且數(shù)值基本穩(wěn)定，這一階段快于如圖11(a)所示的單獨(dú)NVT弛豫，且其在100~200 ps弛豫時(shí)間內(nèi)波動(dòng)范圍明顯小于如圖11(b)所示的單獨(dú)NPT弛豫。在200 ps的弛豫時(shí)間下，單獨(dú)進(jìn)行NPT弛豫識(shí)別得到的512籠平均值更大，單獨(dú)進(jìn)行NVT弛豫識(shí)別得到的512籠平均數(shù)量最少。以上結(jié)果進(jìn)一步表明，NPT弛豫有利于籠型結(jié)構(gòu)的形成，但其結(jié)構(gòu)穩(wěn)定性較差，而NVT弛豫可優(yōu)化穩(wěn)定水合物結(jié)構(gòu)，兩者結(jié)合后有利于水合物弛豫平衡。

圖11 方案4-3 (a) 及方案4-2 (b) 與方案4-1弛豫的512籠識(shí)別結(jié)果對(duì)比Fig. 11 Comparison of 512 cage recognition results of relaxation between scheme 4-1 and scheme 4-3 (a) and scheme 4-2 (b)

此外，在NPT弛豫過(guò)程中，由于體系體積處于波動(dòng)變化，故以分段密度變化考察體系平衡狀態(tài)，對(duì)所選方案的弛豫結(jié)果進(jìn)行驗(yàn)證。對(duì)于NPT弛豫過(guò)程，從體系基本穩(wěn)定（512籠達(dá)到31個(gè)）后開(kāi)始每隔10 ps取一次密度的平均值，計(jì)算每段平均值的標(biāo)準(zhǔn)差，當(dāng)標(biāo)準(zhǔn)差小于0.001 g/cm3時(shí)[25]，認(rèn)為體系達(dá)到平衡狀態(tài)。對(duì)于單獨(dú)NPT弛豫，從模擬20 ps后（表2）開(kāi)始取密度平均值；對(duì)于NVT + NPT弛豫，從NVT弛豫結(jié)束后，即100 ps開(kāi)始取密度平均值。密度計(jì)算結(jié)果如表5所示，單獨(dú)NPT弛豫在200 ps內(nèi)分段密度標(biāo)準(zhǔn)差仍大于0.001 g/cm3，到500 ps時(shí)分段密度達(dá)到平衡狀態(tài)。而首先進(jìn)行NVT弛豫，再進(jìn)行NPT弛豫，可在200 ps時(shí)達(dá)到密度的平衡狀態(tài)。因此，與水合物籠型結(jié)構(gòu)數(shù)量的變化類(lèi)似，從分段密度的角度衡量，將NVT系綜及NPT系綜結(jié)合進(jìn)行弛豫，水合物結(jié)構(gòu)更易達(dá)到平衡狀態(tài)。

2.5 雙系綜弛豫后水合物分解溫度的模擬分析

經(jīng)不同弛豫方案的模擬比較，擬推薦采用雙系綜弛豫方案4-1，即首先進(jìn)行100 ps的NVT弛豫，再進(jìn)行100 ps的NPT弛豫，故對(duì)雙系綜弛豫后處于平衡態(tài)的水合物構(gòu)型進(jìn)行分解模擬以驗(yàn)證弛豫方案的適用性。采用空隙誘導(dǎo)法進(jìn)行水合物分解模擬，Choudhary等[6]的模擬表明，水合物分解溫度在空隙率為10%左右時(shí)處于平衡階段，此時(shí)的分解溫度可確定為該力場(chǎng)下的水合物分解溫度?；陔p系綜弛豫后的水合物構(gòu)型，隨機(jī)移除H2O分子及CH4分子共計(jì)115個(gè)，獲得空隙率約為8%的水合物結(jié)構(gòu)，如圖12所示，橙線(xiàn)為氫鍵，綠球?yàn)镃H4分子中的碳原子，藍(lán)線(xiàn)為隨機(jī)移除H2O分子后破壞的氫鍵，藍(lán)球?yàn)殡S機(jī)移除的CH4分子。

水合物分解模擬過(guò)程相關(guān)參數(shù)設(shè)置與弛豫過(guò)程一致，如1.2節(jié)所述。模擬壓力為10 MPa，模擬溫度由隨機(jī)移除分子后的平衡溫度185 K勻速升溫至350 K，升溫速率為0.165 K/ps，總模擬時(shí)長(zhǎng)為10 ns。體系的溫度及勢(shì)能變化如圖13所示，勢(shì)能突變位置即為相變點(diǎn)，故該模型在10 MPa下的分解溫度約為292 K。對(duì)比10 MPa下水合物實(shí)際分解溫度286.2 K，模擬分解溫度高約5.8 K。這一偏差是由于水合物分解溫度模擬研究依賴(lài)于模擬體系大小，3 × 3 × 3 sI型CH4水合物晶胞模擬系統(tǒng)較小，而較小體系模擬所得到的水合物分解溫度較高[6]。此外，Choudhary等[6]使用3 × 3 × 3 sI型CH4水合物晶胞進(jìn)行分解模擬時(shí)，所得分解溫度偏差范圍為5~35 K，故本模擬結(jié)果在偏差范圍內(nèi)且偏差較小。對(duì)于水合物分解過(guò)程，某一水合物籠坍塌后，水合物開(kāi)始分解，如圖14(a)所示（6.5 ns時(shí)），在隨后的1 ns內(nèi)完全分解，由于分解后CH4分子過(guò)飽和，形成氣液兩相體系，如圖14(b)所示（7.3 ns時(shí)），其中橙線(xiàn)為氫鍵及H2O分子，綠球?yàn)镃H4分子中的碳原子。以上模擬結(jié)果表明，雙系綜弛豫方案所得的穩(wěn)定水合物構(gòu)型可用于后續(xù)產(chǎn)生相的模擬分析，該弛豫方案可作為產(chǎn)生相模擬的基礎(chǔ)。

圖14 水合物分解模擬中6.5 ns (a)及7.3 ns (b)時(shí)體系構(gòu)型Fig. 14 System configuration at 6.5 ns (a) and 7.3 ns (b)during the hydrate decomposition simulation

3 結(jié)論

在進(jìn)行產(chǎn)生相分子模擬前，需要對(duì)體系進(jìn)行弛豫以合理化初始結(jié)構(gòu)。對(duì)于不同體系與模擬需求，弛豫過(guò)程的設(shè)置不盡相同。

（1）需選取適宜的弛豫時(shí)間。一定時(shí)間弛豫后，構(gòu)型變化并不明顯，過(guò)長(zhǎng)弛豫時(shí)間會(huì)顯著增加計(jì)算時(shí)長(zhǎng)。

（2）在NVT系綜下，弛豫有利于提高體系的結(jié)晶度，更低的弛豫溫度可更有效優(yōu)化并穩(wěn)定構(gòu)型，但溫度的設(shè)置需綜合考慮所選力場(chǎng)模型的相平衡條件及冰點(diǎn)。在體系基本穩(wěn)定后（合理的弛豫時(shí)間下），溫度影響不再顯著。故應(yīng)結(jié)合后續(xù)模擬需求考慮，選擇適中的弛豫溫度。

（3）NPT系綜下的弛豫有利于水合物結(jié)構(gòu)的形成，越高的弛豫壓力可在越短弛豫時(shí)間內(nèi)得到可識(shí)別的籠型結(jié)構(gòu)，但其穩(wěn)定性越差。弛豫壓力應(yīng)根據(jù)后續(xù)模擬壓力需求，結(jié)合相平衡條件及選定的弛豫溫度而確定。

（4）不同弛豫系綜所產(chǎn)生的弛豫效果不同。在穩(wěn)定溫度的基礎(chǔ)上，NVT弛豫更利于結(jié)構(gòu)的穩(wěn)定，而NPT系綜對(duì)籠型結(jié)構(gòu)的形成具有促進(jìn)作用。在體系弛豫時(shí)可優(yōu)先考慮NVT弛豫以平衡溫度與穩(wěn)定結(jié)構(gòu)，再結(jié)合NPT弛豫以促進(jìn)水合物結(jié)構(gòu)達(dá)到更接近理論數(shù)量的穩(wěn)定狀態(tài)。

因此，對(duì)于本文3 × 3 × 3 sI型CH4水合物晶胞，250 K的弛豫溫度及10 MPa的弛豫壓力可以得到較好的弛豫效果。在此溫度和壓力條件下，首先進(jìn)行100 ps NVT弛豫，再進(jìn)行100 ps NPT弛豫，即共200 ps的弛豫處理，所得到的平衡結(jié)構(gòu)可基本擬合水合物分解溫度及分解過(guò)程，具備開(kāi)展產(chǎn)生相分子模擬的適用性。

符號(hào)說(shuō)明

p―弛豫壓力，MPa；T―弛豫溫度，K；q―原子電荷值，e；σ―L-J勢(shì)函數(shù)參數(shù)，nm；ε―L-J勢(shì)函數(shù)參數(shù)，kJ/mol。