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

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

引 言

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

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

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

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

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

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

1.1 生成反應(yīng)

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

1.2 基本假設(shè)

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

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

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

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

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

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

1.3 控制方程

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

水相

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

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

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

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

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

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

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

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

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

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

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

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

通過與實驗數(shù)據(jù)對比驗證建立的甲烷水合物生成傳熱數(shù)學(xué)模型準(zhǔn)確性。根據(jù)文獻[15-16]報道的多孔材料水合物儲氣實驗進行建模。圖1為實驗中反應(yīng)釜的模型，上方為氣體填充區(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）作為多孔填料進行水合物儲甲烷實驗?zāi)M驗證[15-16]。水合物設(shè)置選用冰的基本參數(shù)，再根據(jù)甲烷水合物的物理參數(shù)進行修改；在單元區(qū)域條件中根據(jù)多孔材料參數(shù)設(shè)置多孔區(qū)域條件，將多孔材料區(qū)域和氣體填充區(qū)域的交界面設(shè)置為interior 邊界條件。反應(yīng)釜壁面無滑移且設(shè)置為溫度壁面條件。采用標(biāo)準(zhǔn)初始化，兩組模擬的運行參數(shù)見表1。

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

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

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

2.3 數(shù)值模擬方法

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

3 模擬結(jié)果及分析

3.1 模型驗證

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

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

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

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

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

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

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

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

圖7 為泡沫鋁內(nèi)z=0 平面甲烷水合物體積分?jǐn)?shù)隨時間的變化，可見在泡沫鋁內(nèi)的水合物生成傳熱過程與聚氨酯泡沫中明顯不同。前10.0 min，泡沫鋁材料區(qū)域水合物的體積分?jǐn)?shù)保持均勻分布。這是因為初期釜內(nèi)溫度較低，泡沫鋁材料又能提供較多的成核位點[15]，水合物可以在整個區(qū)域內(nèi)均勻生成。由于泡沫鋁熱導(dǎo)率為聚氨酯泡沫的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內(nèi)z=0平面甲烷水合物體積分?jǐn)?shù)隨時間的變化Fig.7 Methane hydrate volume fraction at z=0 plane with time variation in AF foam system

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

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

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

為了能夠更加直觀比較兩種多孔材料的導(dǎo)熱能力，對不同高度的水合物溫度數(shù)據(jù)進行匯總，圖10 為PU 和AF 體系內(nèi)水合物溫度計算值取點位置，在反應(yīng)釜的中線位置，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 體系內(nèi)不同高度甲烷水合物溫度隨時間的變化，可以明顯看出兩個體系最高溫度相差9.0 K。甲烷水合物生成時會放出大量的熱，造成周圍溫度上升。由于多孔材料填充區(qū)域是水合物生成區(qū)，溫度波動比較大，氣體填充區(qū)域溫度波動較為平緩。在PU 體系中，由于甲烷水合物初期大量生成，造成溫度上升，而PU 的熱導(dǎo)率較小，生成熱無法及時排出，從而出現(xiàn)了局部高溫。過高的溫度不僅減緩了水合物的生成，還造成周圍水合物的分解，分解吸收熱量又引起了溫度的下降。故而出現(xiàn)了溫度驟升驟降的現(xiàn)象。而在AF 體系中，AF 的熱導(dǎo)率高，熱量傳遞快，水合物生成熱及時通過泡沫鋁骨架導(dǎo)向壁面，故圖中的溫度波動更加平緩，溫差更小。PU 體系中，從5.0 min 到30.0 min，水合物生成中心區(qū)域一直保持較高溫度，最高溫度甚至升高至283 K，直到水合反應(yīng)基本完成溫度才開始下降，而氣體填充區(qū)域溫度則基本保持在273.0～274.0 K 內(nèi)；而AF 體系中，AF 填充的水合區(qū)域溫度基本保持在273.0～274.0 K 內(nèi)，該區(qū)域溫度在快速反應(yīng)階段略有升高但也很快降至環(huán)境溫度，反應(yīng)釜整體最高溫升不超過1.5 K。

圖10 聚氨酯泡沫PU和泡沫鋁AF體系內(nèi)水合物溫度計算值取點位置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

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

4 結(jié) 論

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

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

（2）水合反應(yīng)過程中，熱量傳遞是影響水合物生成速率的關(guān)鍵因素之一。

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

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