一種新的水合物有效導熱系數球形串并聯模型

孫始財 谷林霖 楊震東 李艷敏

山東省土木工程防災減災重點實驗室·山東科技大學土木工程與建筑學院

0 引言

天然氣水合物（以下簡稱水合物）主要存在于海洋沉積物以及陸地凍土地區，因其全球儲量巨大（約為2×1016m3）[1]，已成為國際認可的非常規戰略能源之一。自2002年，水合物開發從勘查階段進入勘查試采一體化階段，近年正逐步向商業開采階段邁進[2]。然而，由于水合物穩定域邊界很敏感，在開發過程中溫度發生變化時容易引起水合物分解或水合物二次生成，可能導致海底塌陷、滑坡和加劇溫室效應等不可逆的環境問題[3-5]或管道堵塞、井噴和爆炸等安全問題[6]，掌握水合物的熱物理學性質是安全、高效開發水合物資源的基礎。導熱系數是水合物重要的熱物理學性質參數之一，水合物分解前緣主要受到傳熱傳質的影響[7]，而傳熱的主導方式是導熱，因而有必要聚焦水合物導熱系數研究。

目前對水合物導熱系數的研究主要集中在實驗室測量和分子動力學模擬。例如，姚貴策等[8]實時測量在253～270 K、7 MPa時未經壓實的甲烷水合物導熱系數為0.44～0.46 W/(m·K)。黃犢子等[9-10]用十二烷基苯磺酸鈉（SDS）溶液合成了甲烷水合物樣品，并測得未經壓實和通過2 MPa軸向壓力壓實的甲烷水合物導熱系數分別約為0.33 W/(m·K)和0.57 W/(m·K)。上述研究[8-11]基本采用不壓實或較低壓力壓實的水合物樣品，其測得的甲烷水合物樣品導熱系數介于0.33～0.57 W/(m·K)。實際上，這些樣品并不是純甲烷水合物，其水合物固體顆粒之間存在孔隙，而孔隙中往往存在氣體或者液體（水），從而影響導熱系數測量精度。冰粉合成水合物（種冰法）的理論轉化率可達到近100%[12]，因此作為水合物制備方法而被廣泛使用，但其樣品中仍不可避免地存在或多或少的孔隙[13]。研究者們為了得到更純甲烷水合物的導熱系數，對冰粉進行預先壓實[14]或對已制備的樣品進行高壓壓實[15-17]，以減小水合物顆粒間隙。然而，即使經過這樣處理的樣品也無法保證沒有孔隙，即不能保證測量結果是純水合物的導熱系數。

為此，不少學者采用分子動力學模擬純甲烷水合物的導熱系數，其結果基本在0.56～0.75 W/(m·K)[18-21]，與溫度成正相關性[18-19,22]；且在較高壓力區間其與壓力有顯著的正相關[22]，而在相同溫壓條件下，晶穴占有率越大，導熱系數越大[18,22-23]。雖然分子動力學模擬從理論上排除水合物孔隙的影響，并考慮了晶穴占有率、甲烷分子與水籠子的碰撞和范德華作用強度等因素，但其計算過程復雜，過于理論化且使用不方便。

綜上，實驗測量結果與其測量樣品的品質關系甚大，而樣品的品質由實驗步驟、實驗方法等條件綜合決定。而分子動力學模型從理論上模擬水合物導熱系數，其過程不依賴實驗，但計算復雜，不利于工程實際應用。因此，筆者提出一種新的便于計算的水合物有效導熱系數球形串并聯模型，可以為水合物資源開發和衍生技術應用提供技術支撐。

1 球形串并聯模型

對自然界水合物導熱系數進行原位測量難度大，而實驗室非原位測量需要保壓保溫取樣保存技術[24]，其實施過程復雜，所以實驗室制備與原位測量一體化的方案更具有可行性。但實驗室制備的甲烷水合物樣品往往不是純甲烷水合物。因為甲烷水合物樣品的固有特性和實驗室樣品制備方法的有限性，導致甲烷水合物樣品中水合物顆粒間不可避免地存在孔隙。由于孔隙中殘存著氣體或水，其測得的導熱系數不是純甲烷水合物導熱系數[25]。例如，水的導熱系數約0.61 W/(m·K)[26]，與甲烷水合物導熱系數接近，而氣體的導熱系數0.03 W/(m·K)[26]比甲烷水合物導熱系數小一個數量級，故氣體的存在導致測得的導熱系數偏小。因此，筆者提出一種新的水合物導熱系數球形串并聯模型，以消除水合物顆粒間孔隙的影響。

1.1 模型假設

在樣品制備初期，水合物顆粒間的孔隙比較大；隨著水合物的生成，孔隙體積逐漸減小。當孔隙尺寸縮小到一定程度時，水合物逐漸停止生長[27-29]，此時水合物樣品內部的孔隙結構確定。本文認為此時水合物顆?？紫吨袩o流動的反應氣和反應水，系統處于穩定。為了簡化，將水合物顆粒間的孔隙假設為球形，且在樣品中隨機分布，球形孔隙外的區域皆為密實水合物，幾何模型如圖1-a所示。基于本文文獻[27-29]，水合物樣品中球形孔隙半徑最小值為9.2 nm。假設半徑為9.2 nm的球形孔隙是甲烷水合物不能繼續生成的最小孔徑，認為小于此孔徑的孔隙中不會再有水合物生成。熱阻網絡圖如圖1-b所示。由于球形孔隙的對稱性，幾何模型僅繪制了球形孔隙的1/8，并以邊長為L的立方體作為研究對象。通過此模型可以計算不同孔隙度下水合物樣品熱阻和導熱系數。球形孔隙的球心位于坐標軸原點，球形孔隙半徑為r1，熱流量Q方向沿z軸正方向，立方體的熱阻網絡圖如圖1-b所示。根據熱流路徑可將模型分為兩個區域分別計算。Ⅰ區熱流依次經過孔隙流體相（氣相或者水相）和水合物相，則區域Ⅰ的等效熱阻RI為孔隙流體相熱阻R1和水合物相熱阻Rh的串聯熱阻（圖1-b中左支路）；Ⅱ區僅含水合物相，即熱流僅途徑水合物相，故Ⅱ區等效熱阻RII為水合物相熱阻Rh（圖1-b中右支路）。Ⅰ區和Ⅱ區等效熱阻并聯獲得模型的總熱阻。應注意，熱阻與各區域（相）的導熱系數、熱流方向厚度及垂直熱流方向的截面積有關。因此，水合物相熱阻Rh和孔隙流體相熱阻R1會因其所屬區域（相）幾何特點不同而不相等。各區域（相）的幾何特點是由r1和L決定的。

圖1 球形串并聯模型示意圖

1.2 模型建立

多相、多組分材料的導熱能力通常用有效導熱系數描述。本文中有效導熱系數皆指包含水合物及孔隙流體（氣體或水）雙相的水合物樣品導熱系數。式（1）、（2）分別代表傅里葉定律和熱阻原理：

式中λe表示水合物樣品的有效導熱系數，W/(m·K)；δ表示沿熱流方向的單元體厚度，m；A表示垂直于熱流方向的單元體截面積，m2；Rmol表示模型的總熱阻，K/W；Rj表示j區域的等效熱阻，K/W。

首先根據圖1幾何模型和熱阻網絡模型積分得到區域Ⅰ、Ⅱ的等效熱阻方程為：

式中λh和λ1分別表示甲烷水合物和孔隙流體導熱系數，W/(m·K)；μ1,h表示僅與水合物及孔隙流體自身導熱系數相關的無因次中間變量；r1表示球形孔隙的孔徑，m；L表示單元體邊長，m；β表示孔隙度，指水合物樣品中孔隙體積V1占樣品總體積V0的百分數，即

Turner等[28]研究發現，當沉積物孔隙半徑大于60 nm時，甲烷水合物相平衡不會發生變化，即影響水合物形成條件的最大孔隙半徑為60 nm。因此當L為60 nm、r1為9.2 nm時，計算得到水合物樣品的最小孔隙度β為0.001 9。將式（2）～（5）代入式（1），得到求解水合物樣品有效導熱系數的球形串并聯模型方程：

1.3 參數取值

由于分子模擬的結果差異較大且過于理想化[18-23]；而實驗測量結果依賴于樣品品質和測量方法，這為參數選取制造了困難。根據式（6）可知，欲反算出純水合物導熱系數需要樣品有效導熱系數（λe）、樣品孔隙度（β）及孔隙流體導熱系數（λ1）這3個參數。Zhao等[14]采用預先壓實的冰粉在溫度273.35 K和壓力10 MPa環境下合成甲烷水合物樣品，并使系統溫度在263.15 K和278.15 K之間周期振蕩超過3次，最后對樣品老化3天以上，從而確保冰粉完全轉化為水合物。隨后，在確保整個測試過程水合物樣品不分解的前提下，采用雙螺旋結構的Hot Disk探頭測量不同溫度壓力下的甲烷水合物的導熱系數。其預先壓實的冰粉具有與冰相同的導熱系數，很好地消除了冰粉顆粒間隙的影響[14]；其后續溫度震蕩、老化等步驟使可能殘留的水（冰粉）得以完全轉化為水合物。這種樣品制備方法亦未采用SDS等表面活性劑，減少了測量干擾因素[9,15,17]。同時，Sun等[30-31]采取類似的方法制備水合物樣品進行非原位測量，認為該方法制備的樣品純度高，且甲烷水合物樣品的孔隙度介于19.8%～17.8%，故取平均孔隙度18.8%作為樣品孔隙度（β）。將Zhao等[14]在263.15～279.15 K測得的水合物樣品導熱系數0.484 6 ～ 0.494 6 W/(m·K)作為樣品有效導熱系數（λe）。由于此樣品在制備及測量過程中基本處于無液態水系統下，判斷水合物顆粒間孔隙流體為甲烷氣，故孔隙流體導熱系數（λ1）取值 0.03 W/(m·K)[26]。用式（6）計算得到純甲烷水合物導熱系數隨溫度變化關系（圖2），擬合這些計算值可以得到純甲烷水合物導熱系數的半經驗方程為：

圖2 純甲烷水合物導熱系數與溫度變化圖

式中T表示溫度，K。

式（7）在4.0～9.5 MPa壓力范圍內計算純甲烷水合物導熱系數的誤差為±0.007 W/(m·K)。Zhao等[14]認為在2～10 MPa范圍內，氣相壓力對導熱系數的影響很?。籖osenbaum等[20]亦發現在3.8～14.2 MPa范圍內，壓力與導熱系數相關性很小。故式（7）忽略一定范圍的壓力對導熱系數的影響。

將式（7）代入式（6），得到甲烷水合物樣品有效導熱系數球形串并聯模型為：

式中B=λh/λ1-1 ；C1=1.948 89 ；C2=1.208 99 ；C3=1.240 70。根據孔隙流體的不同（水或者氣體），λ1可取不同的參數，也可是關于溫度的函數。式（8）適用于計算水合物顆粒間孔隙流體為單一相的水合物樣品有效導熱系數。當水合物處于生成或分解過程時，孔隙流體是氣液兩相共存[32]，可能發生微對流傳熱，不再是簡單的導熱問題，亦不是本文討論的對象。

2 模型驗證

由于標準的純水合物樣品及其導熱系數尚未見公開報道，這給模型驗證帶來一定的困難。實驗測量結果受到測量方法、樣品制備方法和水合物顆粒間孔隙等多種因素影響，而純理論的分子動力學模擬忽略了一些實際因素而不能完全反映真實情況。鑒于此，把實驗測量結果和分子動力學模擬結果作為參考進行驗證。首先將純甲烷水合物導熱系數半經驗方程計算結果與實驗測量結果[17,20,33]和分子動力學模擬結果[18,20,22,34]進行對比分析。從圖3可以看出，Rosenbaum等[20]采用瞬態平面熱源法（TPS）測得的261.5～277.4 K、4.6～6.3 MPa條件下壓實的樣品導熱系數為0.68±0.01 W/(m·K)，比同溫度下半經驗方程[式（7）]計算結果低 0.08±0.01 W/(m·K)。除了測量誤差的原因，也可能是在有限的實驗溫度壓力范圍內，逐漸縮小的孔隙阻礙了水合物繼續形成，導致樣品中存在少量的低導熱系數氣體而使實驗結果偏小。Gupta等[33]和Li等[17]結合實驗測量結果計算出零孔隙度純甲烷水合物的導熱系數分別是 0.70±0.04 W/(m·K)和大于 0.70 W/(m·K)。由于他們均未明確說明溫度范圍，猜測可能是實驗溫度的范圍，即277.3～279.1 K和244.58～278.43 K。如圖3所示，半經驗方程計算結果比同溫度范圍內Li等[17]和Gupta等[33]的計算結果略高0.03～0.07 W/(m·K)。由于文獻的計算模型沒有考慮溫度對導熱系數的影響，本文模型考慮了此影響，而使模擬結果隨溫度增大而增大。但總體上來說，模擬結果比較接近文獻結果。

圖3 部分文獻與本文得到的純水合物導熱系數圖

除了文獻實驗測量結果，圖3還給出了部分文獻的分子動力學模擬結果[18,20,22,34]。從圖3可以看出，分子模擬結果與實驗測得的導熱系數變化規律基本一致。但是分子模擬結果之間差異比較大，僅在265～280 K溫度范圍較為接近，因此采用該溫度段分子模擬結果進行驗證。從圖3可以看出，半經驗方程計算結果（圖3中直線）比楊德偉等[18]的分子模擬結果約低0.01 W/(m·K)，比Rosenbaum等[20]分子模擬結果約高0.03 W/(m·K)，說明式（7）在263.15～279.15 K范圍的計算結果比較準確。

另外，采用球形串并聯模型［式（8）］計算本文文獻[14]條件下甲烷水合物有效導熱系數與文獻實驗測量結果的相對誤差為0.05%～0.24%（表1），說明本文模型計算精度比較高。

表1 球形串并聯模型模擬結果及誤差表

式中εerror表示相對誤差；λexper表示實驗測得的甲烷水合物導熱系數，W/(m·K)；λmodel表示球形串并聯模型計算結果，W/(m·K)。

3 模型應用

本文提出的球形串并聯模型既可以計算含不同孔隙流體及不同孔隙度的水合物樣品有效導熱系數，也可以根據實驗測得的有效導熱系數計算實驗樣品的孔隙度，籍以分析樣品制備的純度。

3.1 水合物樣品有效導熱系數

利用式（8）計算了孔隙僅含氣或僅含水的水合物有效導熱系數，λ1分別取值0.03 W/(m·K)[26]或0.57 W/(m·K)[35]，λh取值 0.765 4 W/(m·K)，結果如圖4所示。由于水的導熱系數與純水合物的導熱系數接近，所以孔隙流體僅為水時，水合物有效導熱系數與純水合物導熱系數接近。而氣體的導熱系數遠小于純水合物導熱系數，因而孔隙流體僅為氣體時，水合物有效導熱系數與純水合物導熱系數相比變化很大。相比于孔隙度參數，孔隙流體性質對水合物導熱系數的影響更大。所以，樣品制備過程如果水合物固體顆粒間隙比較大，必然影響測量精度。

圖4 有效導熱系數與孔隙度關系圖

3.2 水合物樣品的孔隙度分析

實驗室中水合物樣品制備方法直接決定樣品品質，進而影響導熱系數的測量精度。選取有代表性的幾組水合物導熱系數[8-9,14-17]，利用本文提出的球形串并聯模型分別計算其孔隙度，結果如圖5和表2所示，分析樣品制備純度。

表2 甲烷水合物樣品的平均孔隙度與浮動值的計算結果表

圖5 各文獻樣品孔隙度分析圖

從圖5及表2可知，姚貴策等[8]實驗采用的樣品未經壓實，其有效導熱系數最小，孔隙度介于20%～25%，平均孔隙度最大。而Waite等[16]持續用102 MPa壓實4天所得樣品的有效導熱系數最大，計算得到的平均孔隙度明顯小于其他樣品。Li等[17]在 4 MPa、12 MPa和 16 MPa 壓實壓力下樣品的平均孔隙度依次是17.23%、15.08%和14.05%，很明顯地隨壓實壓力的升高而減小。說明壓實方法在一定程度上能夠制備更純的樣品，但仍無法達到零空隙。溶解氣法[9,15,17]制備的水合物樣品孔隙度稍低（10%～18%），但是表面活性劑不參與水合物反應而殘留在水合物顆粒間孔隙中可能會影響測量結果，進而使計算得到的孔隙度偏小。在沒有很高的壓實壓力時，溶解氣法制備的樣品中殘余氣體較少。但考慮其樣品可能受到活性劑的污染，種冰法[14]樣品制備相對更好。從圖5和表2可以看出，預先壓實種冰法制備的樣品[14]平均孔隙度為18.73%，浮動值為±0.002 3，說明每次制備的樣品純度基本一致，這樣測量的結果更可靠。同樣是采用種冰法制備樣品，Waite等[16]實驗采用的樣品經過壓實后平均孔隙度僅為7.34%，但是浮動值明顯大于預先壓實冰粉的樣品[14]的浮動值，表明每次樣品的純度差異較大。種冰法可以使氣體通過冰粉顆粒間隙相對均勻擴散，形成的水合物樣品也相對均勻；而預先壓實冰粉后冰粉顆粒更加密實，水合物樣品中孔隙相對均勻和更小。在用SDS溶液制備的樣品[9,15,17]中，黃犢子[9]實驗采用的樣品其孔隙度的浮動值低，即樣品相對均一。Li等[17]壓實壓力和平均孔隙度都高于黃犢子[9]，但其浮動值接近，這可能是壓實過程中造成更多裂隙，甚至可能造成水合物晶體破損。李棟梁[15]實驗過程與Li等[17]類似，但浮動值差異較大；姚貴策等[8]實驗采用的樣品和Waite等[16]的制備方法差異很大，但浮動值接近。這些情況說明樣品的孔隙度受多種因素影響。

基于上述分析，對水合物樣品制備方法可以得到一些初步的認識。在每次制備的水合物樣品差異程度上，預先壓實冰粉優于一般壓實操作；在水合物樣品純度上，壓實樣品優于未壓實樣品，種冰法優于溶解氣法。

4 結論

1）本文提出的半經驗方程在263.15～279.15 K的溫度范圍內可以比較準確地計算純水合物導熱系數，在4.0～9.5 MPa壓力范圍內的誤差為±0.007 W/(m·K)。

2）利用球形串并聯模型計算預先壓實種冰法制備的水合物樣品導熱系數的值與實驗測量結果誤差為0.05%～0.24%，說明球形串并聯模型可以用于預測已知孔隙度、孔隙流體類型及溫度的樣品有效導熱系數。計算了含不同孔隙流體及不同孔隙度的水合物樣品有效導熱系數，發現孔隙流體類型對有效導熱系數影響更為明顯。

3）通過分析不同方法制備的水合物樣品的孔隙度和純度，初步獲得了高純度水合物樣品制備方法的一些認識，即未經壓實的水合物樣品孔隙度介于20%～25%，普通壓實后的樣品孔隙度介于7.34%～17.23%，孔隙度總體上隨壓實壓力的增大減小，從而可以提高樣品純度；樣品中殘余的表面活性劑使樣品孔隙度減小，但是影響樣品純度；預先壓實種冰法無需表面活性劑，可很好地制備高品質水合物樣品。