郎雪梅，姚柳眉，樊栓獅，李剛，王燕鴻，2

（1 華南理工大學化學與化工學院，廣東 廣州 510640； 2 廣東省先進絕緣涂料工程技術研究中心，廣東 珠海 519175；3 廣東省燃料電池技術重點實驗室，廣東 廣州 510640）

引 言

天然氣作為21 世紀最重要的能源之一被廣泛應用于交通運輸和生產(chǎn)生活[1]。隨著天然氣需求的不斷增加，安全、經(jīng)濟的天然氣儲運技術也更為重要。目前天然氣的儲運技術主要有管道運輸、液化儲運（LNG）、壓縮儲運（CNG）。管道儲運廣泛應用于工業(yè)和居民生活，適合長期持續(xù)的遠距離輸送，但工程投資高，靈活性差。天然氣液化儲運能夠提供約600 m3CH4(STP)/m3的高儲存容量[2]，目前主要用于遠洋輸送方面，但儲存溫度低以及連續(xù)不斷的氣體蒸發(fā)對生產(chǎn)設備和儲存裝置要求較高[3]。天然氣壓縮儲運是將天然氣加壓到20～25 MPa 進行儲運，適合于中小規(guī)模的儲運，靈活性高。但高壓儲運增加了運輸風險[4]。水合物儲運（NGH）是近幾年發(fā)展起來的天然氣儲運技術，因其固態(tài)儲存、儲氣密度高、經(jīng)濟性高、安全性高得到越來越多的關注[5-6]。

水合物儲運是利用水和天然氣在低溫高壓條件下生成固體水合物進行天然氣儲運的技術。在標準條件下1 m3水合物可以儲存180 m3的天然氣[7-9]，且可以在冰點以下和大氣壓條件下儲存數(shù)月[10-12]。根據(jù)日本三井公司的測算，在500 km 范圍內NGH 運輸成本比LNG 低18%～24%[13]，適合于短距離小規(guī)模的天然氣運輸。但天然氣水合物的生成速率慢，限制了NGH 技術的應用[14]。Fan 等[15-16]利用泡沫鋁、聚氨酯泡沫等多孔材料與動力學促進劑耦合，突破了水合物形成速度的障礙，實驗室內實現(xiàn)了20 min 內將水合儲甲烷密度提升至接近理論值。

從研究角度看，水合物儲運已具備實現(xiàn)工業(yè)化的潛力。但水合物的生長是傳質傳熱控制的反應，因此在放大實驗中存在諸多不確定因素[17]，利用Fluent 數(shù)值模擬可以將水合反應中的復雜生成傳熱過程數(shù)值化，清晰了解反應器內部熱量分布及反應程度，從而更好地指導水合儲氣反應進行放大。

Fluent 模擬可用于研究水合物分解[18]和水合物生成[19]。Neto 等[20]開發(fā)了一種預測海底管道水合物形成的Fluent 模型，該模型運用雙膜理論探究氣液相傳質過程，基于傳質動力學模型搭建水合物生成數(shù)學模型，但模擬結果未與實驗驗證。Song 等[21]采用了滲透率模型和反應動力學模型來研究砂質沉積物中水合物的非均勻生成，結果顯示初始流體分布不同對水合物生成過程有很大影響，模擬與實驗結果吻合良好。但關于水合物在反應器內的模擬研究未見報道。

本文針對水合物反應器中甲烷水合物生成傳熱過程，建立基于化學反應動力學模型和多孔材料內傳質傳熱模型的甲烷水合物生成傳熱數(shù)學模型。通過與Fan 等[15-16]的實驗結果比較驗證數(shù)學模型準確性。該模型有助于分析填充多孔材料的反應器內水合物的生成情況，指導水合物放大反應器的設計。

1 多孔材料內甲烷水合物生成傳熱數(shù)學模型

1.1 生成反應

水合物的生成需要低溫高壓的環(huán)境[9]，是氣體和液體生成固體的一個過程。假設本文溫壓條件下生成的甲烷水合物為sⅠ型水合物[22-23]，水合數(shù)n為5.75[9]。甲烷水合物的摩爾平衡如式(1)所示。

1.2 基本假設

多孔材料中甲烷水合物生成過程涉及溫壓變化、多相流以及化學反應。在數(shù)值模擬計算的過程中，需要應用水合物生成過程中的Arrhenius 方程、質量守恒方程和能量守恒方程。參照天然氣水合物有關的理論和現(xiàn)有的研究成果[24-25]，建立了數(shù)學模型，并做出下面的假設：

(1)忽略體系中的動力學促進劑，假設生成反應中只由液態(tài)水和甲烷氣生成水合物；

(2)模型建立只考慮三相，即水相、甲烷相和水合物相；

(3)忽略甲烷在水中的溶解效應和水蒸氣的存在[26]；

(4)假設多孔材料內的所有相在同一時刻都具有相同的壓力；

(5)忽略輻射換熱和節(jié)流效應，流體與固體之間保持局部熱平衡。

1.3 控制方程

本文選擇基于化學反應動力學模型和多孔材料內傳質傳熱模型[27]，建立甲烷水合物生成傳熱數(shù)學模型。模型中的水相、甲烷氣相和水合物相質量守恒方程如下：

水相

式中，K是多孔材料的絕對滲透率，μm2；krw是水相的相對滲透率，μm2；krg是氣相的相對滲透率，μm2；μw和μg分別為水和氣相的黏度，mPa·s；φ是孔隙率；ρw、ρg和ρh分別為水、氣和水合物相的密度，g/cm3；pg是氣相壓力，MPa；mw和mg分別為單位時間單位體積生成的甲烷水合物中消耗的水質量和氣質量，g/(cm3·s)；mh是單位時間單位體積生成的甲烷水合物質量，g/(cm3·s)；Sw、Sg和Sh分別為水相、甲烷氣相和水合物相的飽和度。

反應釜內部填充多孔材料，因此建立多孔介質模型。假設多孔材料是均勻分布，在模擬過程增加一個源項Si，該源項由兩部分組成，一個黏性損失項α和一個慣性損失項Cn。黏性損失項為水相、氣相和水合物相在各個方向上的滲透率，計算式如下[28]。

式中，Si為動量方程中i方向上的源項，kg/(m2·s2)；α為相對滲透率，m-2；ui為i方向流體速度，m/s；Cn為慣性阻力系數(shù)，m-1；u為流體速度，m/s；Dp為孔徑，m。

不同反應時間水相、甲烷氣相和水合物相飽和度的具體數(shù)值可從Fluent 中得到，它們之間關系如下。

式(4)中的質量源項mh采用Arrhenius 方程進行建模[29-30]，并結合水合物反應界面面積進行校正。模型搭建過程參考了Kim 等[31]的水合物分解動力學模型。當p≥pe時，可作為水合物生成數(shù)學模型。因此甲烷水合物質量源項mh可用式(9)計算。

式(13)有以下假設：水合物生長基于多孔材料孔隙表面；反應界面面積取決于基底面和氣相水相之間的接觸可能性。式中，Ap是多孔材料內的孔隙表面積以及水合反應過程中形成的水合物表面積，m2/m3，使用Kozeny-Carman 方程進行估算[34]，用式(14)計算。

式中，φe為有效孔隙率，定義為氣相和水相體積與全部孔隙體積的比值；Г是水合反應中活躍的孔隙表面積分數(shù)。絕對滲透率具有如下形式[35]，隨水合物飽和度的增加而減小。

多孔材料區(qū)域中氣體、水和水合物之間保持總能量平衡，流體與固體之間保持局部熱平衡。總能量是內能和動能之和，能量平衡也就是控制體內各相內能、相間流動能與反應熱的動態(tài)平衡，因此能量守恒方程經(jīng)過簡化，為

式中，C1和C2均為常數(shù)，C1=5.31×107J/kmol；C2=-15850 J/(kmol·K)。

針對上述數(shù)學模型構建三相質量源項和水合物能量源項的UDF 程序。UDF 程序利用Fluent 軟件中的預定義Define 宏命令實現(xiàn)。使用“#DEFINE”宏命令定義源項中全局變量，DEFINE_SOURCE 宏命令定義程序中三相質量源項和能量源項的變化，因此數(shù)學模型較為復雜，于是采用編譯UDF 的方式實現(xiàn)程序在Fluent中的裝載。

2 數(shù)值模擬條件及方法

2.1 物理模型與網(wǎng)格劃分

通過與實驗數(shù)據(jù)對比驗證建立的甲烷水合物生成傳熱數(shù)學模型準確性。根據(jù)文獻[15-16]報道的多孔材料水合物儲氣實驗進行建模。圖1為實驗中反應釜的模型，上方為氣體填充區(qū)域，下方為多孔材料區(qū)域；多孔材料區(qū)域的高度根據(jù)實驗過程填充多孔材料量和溶液量來確定，用該物理模型進行甲烷水合物生成傳熱模擬實驗。

圖1 物理模型Fig.1 Physical model

網(wǎng)格劃分工具采用Fluent Meshing，驗證模型的網(wǎng)格劃分采用先劃分表面網(wǎng)格再劃分體網(wǎng)格的方法，同時劃分邊界層和加密流體域連接之間的網(wǎng)格。最小網(wǎng)格尺寸為2 mm，體網(wǎng)格劃分采用了Poly-Hexcore網(wǎng)格，圖2為模型的網(wǎng)格劃分。

圖2 模型網(wǎng)格劃分Fig.2 Grid system of model

2.2 邊界條件及初始化

選擇聚氨酯泡沫（PU）和泡沫鋁（AF）作為多孔填料進行水合物儲甲烷實驗模擬驗證[15-16]。水合物設置選用冰的基本參數(shù)，再根據(jù)甲烷水合物的物理參數(shù)進行修改；在單元區(qū)域條件中根據(jù)多孔材料參數(shù)設置多孔區(qū)域條件，將多孔材料區(qū)域和氣體填充區(qū)域的交界面設置為interior 邊界條件。反應釜壁面無滑移且設置為溫度壁面條件。采用標準初始化，兩組模擬的運行參數(shù)見表1。

表1 模擬運行參數(shù)Table 1 Operating parameters of the simulation

多孔材料區(qū)域水合物生成時氣、水、水合物飽和度變化對釜內介質表觀熱導率λa的影響如式(20)。

式中，λw、λg、λh和λp分別為水、甲烷、甲烷水合物和多孔材料的熱導率，W/(m·K)。

2.3 數(shù)值模擬方法

本文采用商業(yè)計算流體力學軟件Ansys Fluent進行數(shù)值模擬，使用瞬態(tài)計算，時間步長根據(jù)實際需求調整。多相流模型選擇歐拉多相流，湍流模型選擇Realizableκ-ε湍流模型（反應釜內的計算單元存在因水合物生成而產(chǎn)生壓降的局部湍流，故使用湍流模型）。采用雙精度求解器、Phase Coupled SIMPLE 算法，除壓力項采用二階精度來提高計算準確性，湍流、動量和能量等方程中對流項均采用一階精度。松弛因子設定保持默認，以標準化殘差合理振蕩波動判斷收斂。

3 模擬結果及分析

3.1 模型驗證

圖3 為文獻[15-16]中甲烷水合物在PU 和AF 多孔材料中的生成動力學結果。用本文描述的數(shù)值模型方法進行模擬，模擬的時間根據(jù)文獻[15-16]中水合物快速生成的時間段進行選擇，分別模擬了33.5 min 和13.5 min。模擬初始填充氣體量、水量和多孔材料量按照實驗設定，模擬過程壁面溫度、反應釜內壓力和多孔材料孔隙度的設定均與兩組實驗過程保持一致。

圖3 甲烷水合物生成動力學[15-16]Fig.3 Kinetic of methane hydrate formation[15-16]

圖4(a)、(b)為甲烷水合物在聚氨酯泡沫PU 和泡沫鋁AF 中的生成動力學模擬數(shù)據(jù)與實驗數(shù)據(jù)的比較，絕對平均相對誤差AARD 分別為4.38% 和5.22%，可以看到實驗和模擬值生成趨勢接近，模擬儲氣量和實驗儲氣量基本一致，表明該數(shù)值模型建模中假設簡化合理，所建甲烷水合物生成傳熱模型適用于該體系。這兩組數(shù)據(jù)也可以明顯看到水合物在泡沫鋁多孔材料中比在聚氨酯泡沫多孔材料中生成更快。

圖4 模擬數(shù)據(jù)與實驗數(shù)據(jù)比較Fig.4 Comparison of simulation data and experimental data

3.2 多孔材料熱導率對水合物生成傳熱過程的影響

圖5 為聚氨酯泡沫內z=0 平面甲烷水合物體積分數(shù)隨時間的變化?？梢钥吹剿衔锷蓞^(qū)域在多孔介質中。從前15.0 min 云圖中看到反應初期水合物生成比較均勻，在整個多孔介質填充區(qū)域均有生成。15.0 min 以后，隨著生成量的增加，水合物呈現(xiàn)中部區(qū)域生成量低于靠近冷壁面區(qū)域，氣液相交界面生成量較多。圖6 為聚氨酯泡沫內z=0 平面甲烷水合物溫度隨時間的變化?？梢钥吹椒磻獌鹊臏囟确植寂c水合物生成量密切相關。10.0 min 之前，系統(tǒng)整體溫度較低，水合物生成較快，釋放出大量的熱。由于聚氨酯泡沫導熱差，隨著水合物生成量增大，多孔材料區(qū)域溫度快速升高，尤其是中心區(qū)域，超過276 K。到15.0 min時，除冷壁面附近，大部分區(qū)域溫度升至276 K 以上。受到溫度升高的熱抑制和釜內壓力降低影響，水合物生成速率開始逐漸減緩，冷量從壁面緩慢傳遞至反應釜中心，釜內溫度開始下降。到33.0 min 左右，水合物生成基本完成，溫度也恢復到環(huán)境溫度。

圖5 聚氨酯泡沫PU內z=0平面甲烷水合物體積分數(shù)隨時間的變化Fig.5 Methane hydrate volume fraction at z=0 plane with time variation in PU foam system

圖6 聚氨酯泡沫PU內z=0平面甲烷水合物溫度隨時間的變化Fig.6 Methane hydrate temperature at z=0 plane with time variation in PU foam system

圖7 為泡沫鋁內z=0 平面甲烷水合物體積分數(shù)隨時間的變化，可見在泡沫鋁內的水合物生成傳熱過程與聚氨酯泡沫中明顯不同。前10.0 min，泡沫鋁材料區(qū)域水合物的體積分數(shù)保持均勻分布。這是因為初期釜內溫度較低，泡沫鋁材料又能提供較多的成核位點[15]，水合物可以在整個區(qū)域內均勻生成。由于泡沫鋁熱導率為聚氨酯泡沫的6000倍，盡管水合物大量生成放熱，但并未出現(xiàn)PU 體系中出現(xiàn)的中心過熱現(xiàn)象，而是在整個AF 填充區(qū)域無較高溫升，見圖8 水合過程溫度云圖。AF 及時將中心區(qū)域的熱傳遞到壁面而被環(huán)境吸收，最高中心溫度為274 K。整個儲氣過程持續(xù)13.5 min 完成，較PU體系縮短20.0 min。

圖7 泡沫鋁AF內z=0平面甲烷水合物體積分數(shù)隨時間的變化Fig.7 Methane hydrate volume fraction at z=0 plane with time variation in AF foam system

圖8 泡沫鋁AF內z=0平面甲烷水合物溫度隨時間的變化Fig.8 Methane hydrate temperature at z=0 plane with time variation in AF foam system

圖9 為水合反應中心的水合物、氣、水飽和度和釜內介質表觀熱導率隨時間的變化，可從看出，隨著水合物飽和度增加和氣、水飽和度下降，兩個多孔材料體系中介質表觀熱導率都隨之增加。反應結束時PU體系和AF體系介質表觀熱導率分別為0.23 W/(m·K)和64.82 W/(m·K)。

圖9 水合反應中心的水合物、氣、水飽和度和介質表觀熱導率隨時間的變化Fig.9 Plot of Sh,Sg,Sw and λa of hydration reaction center with time

為了能夠更加直觀比較兩種多孔材料的導熱能力，對不同高度的水合物溫度數(shù)據(jù)進行匯總，圖10 為PU 和AF 體系內水合物溫度計算值取點位置，在反應釜的中線位置，z=0 mm 和x=0 mm 處，沿y軸均勻取點。均在多孔材料區(qū)域取兩個測溫點，一個測溫點在多孔材料中心區(qū)域，一個測溫點在多孔材料區(qū)域和氣體填充區(qū)域交界面，并按照兩點間距取氣體填充區(qū)域測溫點，因此PU 體系多孔材料區(qū)域每隔20 mm 取一個點，多孔材料區(qū)域2 個位置點，氣體填充區(qū)域4 個位置點；AF 體系每隔18 mm取一個點，多孔材料區(qū)域2個位置點，氣體填充區(qū)域5個位置點。圖11(a)、(b)分別為聚氨酯泡沫PU 體系和泡沫鋁AF 體系內不同高度甲烷水合物溫度隨時間的變化，可以明顯看出兩個體系最高溫度相差9.0 K。甲烷水合物生成時會放出大量的熱，造成周圍溫度上升。由于多孔材料填充區(qū)域是水合物生成區(qū)，溫度波動比較大，氣體填充區(qū)域溫度波動較為平緩。在PU 體系中，由于甲烷水合物初期大量生成，造成溫度上升，而PU 的熱導率較小，生成熱無法及時排出，從而出現(xiàn)了局部高溫。過高的溫度不僅減緩了水合物的生成，還造成周圍水合物的分解，分解吸收熱量又引起了溫度的下降。故而出現(xiàn)了溫度驟升驟降的現(xiàn)象。而在AF 體系中，AF 的熱導率高，熱量傳遞快，水合物生成熱及時通過泡沫鋁骨架導向壁面，故圖中的溫度波動更加平緩，溫差更小。PU 體系中，從5.0 min 到30.0 min，水合物生成中心區(qū)域一直保持較高溫度，最高溫度甚至升高至283 K，直到水合反應基本完成溫度才開始下降，而氣體填充區(qū)域溫度則基本保持在273.0～274.0 K 內；而AF 體系中，AF 填充的水合區(qū)域溫度基本保持在273.0～274.0 K 內，該區(qū)域溫度在快速反應階段略有升高但也很快降至環(huán)境溫度，反應釜整體最高溫升不超過1.5 K。

圖10 聚氨酯泡沫PU和泡沫鋁AF體系內水合物溫度計算值取點位置Fig.10 The point of hydrate temperature in PU and AF porous materials system

圖11 不同高度甲烷水合物溫度隨時間的變化Fig.11 Methane hydrate temperature variation with time at different heights

綜上，水合反應中，熱量傳遞是影響水合物生成速率的關鍵因素之一。填料熱導率越大，越有利于水合生成熱移除以及水合物生成。導熱不良時，容易在中心形成局部過熱，對水合物生成造成熱抑制。

4 結 論

本文建立了多孔材料內甲烷水合物生成傳熱數(shù)學模型，通過與實驗數(shù)據(jù)對比驗證了該模型的可靠性，并對不同熱導率的多孔材料內水合反應過程進行數(shù)值模擬，得到以下結論。

（1）模型驗證中AARD 小于6%，實驗和模擬儲氣量趨勢接近，表明該水合物生成傳熱模型基本準確，可靠性高。

（2）水合反應過程中，熱量傳遞是影響水合物生成速率的關鍵因素之一。

（3）填料的熱導率越大，越有利于水合生成熱的移除以及水合物的生成。導熱不良時，容易在中心部分形成局部過熱，對水合物生成造成熱抑制。

（4）該模型可用于水合反應器的設計和優(yōu)化。在進行甲烷水合反應器設計時，應重點考慮反應器內的傳熱控制。