儲氫合金氫氣釋放過程多物理場耦合建模分析

袁小野

（武漢理工大學(xué)汽車工程學(xué)院 湖北省現(xiàn)代汽車零部件技術(shù)重點實驗室 武漢 430070）

金屬氫化物儲氫是指在一定的環(huán)境溫度和環(huán)境壓強的條件下，固態(tài)氫、氣態(tài)氫氣、儲氫多孔介質(zhì)合金進行質(zhì)量交換、動量交換、并伴隨著有反應(yīng)熱的能量交換的綜合多物理場耦合過程［1］.2002以前，針對上述過程的研究主要以代數(shù)計算方式進行研究［2－4］.2002年以來，Delahaye，Hardy等人在 COMSOL Multiphysics 3.0－4.0的平臺上進行了二維、二維軸對稱的仿真結(jié)果，對于金屬儲氫反映過程的理論研究逐步向可視化、形象化的方向發(fā)展，同時也對金屬儲氫反映過程的偏微分方程的形式進行了進一步的完善［5－6］.然而，考慮到金屬儲氫吸收和釋放反映過程中氫氣密度，固態(tài)氫密度，容器溫度以及壓強等參數(shù)在容器內(nèi)部的變化趨勢主要受時間因素影響.即，對于某一時刻，物理量空間分布的波動性遠小于時間波動性［7］.因此本研究旨在通過二維模型和一維模型的對比分析，驗證一維有限元模型在氫氣存儲效率模型分析中的更廣泛的參數(shù)分析范圍.

1 基本原理

儲氫容器的反應(yīng)過程概括為氣態(tài)氫氣，固態(tài)氫氣和能量3種物質(zhì)及其相關(guān)物理屬性的相互轉(zhuǎn)換過程.其中氣態(tài)氫氣在視為多孔介質(zhì)的金屬儲氫容器中滿足氣態(tài)滲流條件和動量守恒定律，其狀態(tài)滿足Darcy方程.

式中：ε為金屬氫化物的孔隙率；vg為氫氣的達西速度，m／s；ρg為氣態(tài)氫的空間密度，kg／m3.氫氣的滲透過程可以通過方程（2）來描述動量守恒定律，方程（2）中k（m2）為達西滲透率；p為達西壓強，Pa；μ（Pa·s）為粘度；物理速度u＝εvg.同時氣體在空隙中滿足理想氣體狀態(tài)量方程（3）中Mg為氣體摩爾質(zhì)量，kg／mol；Rg為標準氣體狀態(tài)量常數(shù)，J／（mol·K）；T 為溫度，K.

釋放反應(yīng)過程中，質(zhì)量源m與金屬氫化物的顆粒密度ρs（kg／m3）、相變激發(fā)能量Ed（J／mol）、轉(zhuǎn)化率Cd（1／s）、平衡狀態(tài)壓強peq（Pa）、溫度場T（K）有關(guān)，具體函數(shù)關(guān)系可以寫為

固態(tài)氫與氣態(tài)氫氣的轉(zhuǎn)換過程中的質(zhì)量守恒定律可以通過方程（5）與方程（1）聯(lián)立描述

反映過程中的能量傳遞過程符合多孔介質(zhì)中氣體熱傳遞規(guī)律.

式中：ρ為混合狀態(tài)等效密度，kg／m3；Cp為等效熱容，J／（kg·K）；keff為等效熱導(dǎo)率，W／（m·K）；ΔH 為相變潛熱，J／kg.式（6）中左邊第二項為對流傳導(dǎo)，右邊第一項為熱傳遞.

從方程（1）～（6）的聯(lián)立過程，可以總結(jié)出如圖1所示的變量耦合關(guān)系圖.

圖1 釋放反應(yīng)過程變量耦合關(guān)系圖

從圖1中不難看出，氫氣存儲反應(yīng)過程可以由獨立變量描述.從系數(shù)項，質(zhì)量源項對獨立變量的函數(shù)關(guān)系中（式（2）～（4））也可以看出，描述儲氫反應(yīng)過程的偏微分方程組是具有高度非線性程度的.這也對儲氫反應(yīng)過程的高維仿真的收斂性帶來了極大的挑戰(zhàn).

2 數(shù)值仿真模型

金屬儲氫釋放反應(yīng)過程根據(jù)上述基本原理，采用 COMSOL Multi－physics 5.0 中的 Darcy′s Law模塊，Mass diffusion模塊和heat transfer in porous media模塊建立了二維仿真模型和一維仿真模型.對結(jié)果進行了對比分析.同時討論了收斂性和模型的適用范圍.

2.1 幾何模型

仿真對空心圓柱形儲氫元件進行了二維軸對稱和一維建模，分別如圖2，圖3所示.其中，二維幾何模型中，r和z分別表示圓柱坐標系的徑向坐標和軸向坐標，L1＝0.75mm為氫氣釋放邊界，L2＝0.57mm，L3，L4和D1分別為元件頂部，外壁，底部和儲氫介質(zhì)反應(yīng)區(qū)域.忽略物理量在反應(yīng)區(qū)域縱向空間分布的變化情況，可取反應(yīng)區(qū)任意水平截線得出一維幾何模型如圖3所示.圖3中p1，p2和L0分別對應(yīng)儲氫元件的氫氣釋放邊界，儲氫元件外壁和反應(yīng)區(qū)域.

圖2 儲氫元件二維軸對稱幾何模型

圖3 儲氫元件一維幾何模型

2.2 模型參數(shù)

圖1 中所示釋放反應(yīng)過程中的變量耦合關(guān)系的系數(shù)見表1.

表1 釋放反應(yīng)過程模擬參數(shù)列表

2.3 邊界條件

容器上部和側(cè)部邊界設(shè)置為零通量，同時在求解域中設(shè)置質(zhì)量源m.設(shè)置出口邊界為V0用于模仿恒定抽氣速度的工作狀態(tài).在擴散對流PDE模塊中，使用初始條件ρs＝ρs0模擬容器滿載固態(tài)氫時的初始狀態(tài).同時，由于固體在容器中只會通過相變方式轉(zhuǎn)換為氣體，因此設(shè)置所有邊界通量為零且對全部求解域設(shè)置為負數(shù)質(zhì)量源.能量傳遞方程中，初始溫度為T0，假設(shè)儲氫設(shè)備保溫良好，因此設(shè)置全部邊界熱絕緣，且在全部求解域中設(shè)置熱源項模擬反映放熱和加熱過程.

3 結(jié)果與討論

根據(jù)第三部分的參數(shù)設(shè)置和邊界條件設(shè)置，采用自由三角形對二維模型進行網(wǎng)格剖分得到模型空間自由度為29 703的網(wǎng)絡(luò)模型.對一維模型采用單元數(shù)量為500的固定分布網(wǎng)格剖分，得到自由度為2 503的網(wǎng)絡(luò)模型.采用基于BDF算法的時變求解器求解釋放反應(yīng)從初始狀態(tài)到穩(wěn)定狀態(tài)0～30s的穩(wěn)定狀態(tài)的物理場分布，得到如圖4～圖6所示的結(jié)果.

3.1 達西壓強對比

二維軸對稱和一維仿真模型的達西壓強空間分布規(guī)律如圖4所示.可以看到，儲氫元件的最大壓強位置在出口處，最小壓強在儲氫元件的外壁處.這種分布規(guī)律是由于設(shè)定恒定抽氣速度的邊界條件所引起.同時，不難看出，達西壓強的空間分布波動情況小于10－6.因此，可以近似認為，儲氫元件中壓強在空間中呈均勻分布.

圖4 氫氣釋放過程達西壓強空間分布

由二維軸對稱模型和一維模型對比分析可知，二者的最大壓強在同一數(shù)量級，微小的差別可以歸因于一維模型缺失了一個空間維度所致.

進一步地，對氫氣釋放過程中壓強的時間變化規(guī)律進行了研究.圖5的計算結(jié)果展現(xiàn)了釋放反應(yīng)從初始狀態(tài)到穩(wěn)定狀態(tài)0～30s的穩(wěn)定狀態(tài)的達西壓強時變規(guī)律.從圖5中可以看出，反應(yīng)過程開始后，達西壓強從反應(yīng)初始狀態(tài)的5MPa初始壓強，在極短的時間內(nèi)驟降至穩(wěn)定狀態(tài)（約為114kPa）.達到穩(wěn)定反應(yīng)后，隨著時間的推移，儲氫元件中的固態(tài)氫不斷轉(zhuǎn)化為氣態(tài)氫.從圖5的對數(shù)之間坐標圖中可以看出，達西壓強與隨著時間增加呈線性增長趨勢.

3.2 收斂性討論

從圖4的對比分析可以看出，二維模型和一維模型的計算結(jié)論具有極高的吻合程度.由于物理量的空間分布波動性遠遠小于時間波動性且在金屬儲氫元件釋放反應(yīng)過程中工程上更為關(guān)注不同時刻的儲氫元件的整體性能，因此一維模型具有更高的參數(shù)適應(yīng)度和收斂情況.圖6中對比了相同參數(shù)設(shè)置情況下，二維模型和一維模型在0～50min時間內(nèi)的收斂情況.從圖6中不難看出，對于相同的參數(shù)設(shè)置，二維的模型的時間步長倒數(shù)波動極大，當模擬反應(yīng)過程時間超過800s后，模型的時間步長倒數(shù)快速上升，最終會導(dǎo)致模型容差溢出無法解算.相反，一維模型的收斂圖在55s后，時間步長逐漸收斂于恒定值，解算效率是二維模型的數(shù)百倍.一維模型的參數(shù)適應(yīng)性更高，更廣泛地適用于金屬儲氫釋放反應(yīng)過程中儲氫元件物理場分析.

圖5 達西壓強－時間曲線

圖6 二維模型與一維模型的收斂性對比

4 結(jié)束語

金屬儲氫的存儲反應(yīng)過程是一個伴隨吸熱過程的氣態(tài)氫氣到固態(tài)氫的物理相變過程.反應(yīng)過程的效率與初始條件以及熱源條件等因素密切相關(guān).物理狀態(tài)的劇烈變化發(fā)生在初始階段，后續(xù)過程反應(yīng)相對平和.需要在容器設(shè)計過程中采取特殊措施確保反應(yīng)過程初始時的安全性.一維模型和二維模型計算結(jié)果的時間分布具有一致趨勢.以為模型維的收斂性和求解效率遠遠由于二維模型，適用于更寬松的參數(shù)范圍.同時由于物理量空間分布的波動性較小，對于儲氫反應(yīng)過程關(guān)注物理量時變特性的研究，一維模型具有更高的價值和適用性.

［1］ROSI N L，ECKERT J，EDDAOUDI M，et al.Hydrogen storage in microporous metal－organic frameworks［J］.Science，2003，300：1127－1129.

［2］DARKRIM F L，MALBRUNOT P，TARTAGLIA G P.Review of hydrogen storage by adsorption in carbon nanotubes.International［J］Journal of Hydrogen Energy，2002，27（2），193－202.

［3］MAT M D，KAPLAN Y.Numerical study of hydrogen absorption in an Lm－Ni 5hydride reactor［J］.International Journal of Hydrogen Energy，2001，26（9）：957－963.

［4］LEVESQUE D，GICQUEL A，DARKRIM F L，et al.Monte carlo simulations of hydrogen storage in carbon nanotubes［J］.Journal of Physics：Condensed Matter，2002，14（40）：9285.

［5］HARDY B J，＆ ANTON D L.Hierarchical methodology for modeling hydrogen storage systems［J］.Part II（2009）：Detailed models，International Journal of Hydrogen Energy，2009，34（7）：2992－3004.

［6］DELAHAYE A，AOUFI A，GICQUEL A，PENTCHEV I.Improvement of hydrogen storage by adsorption using 2－D modeling of heat effects［J］.AIChE Journal，2002，48（9）：2061－2073.

［7］MUTHUKUMAR P，RAMANA S V.Numerical simulation of coupled heat and mass transfer in metal hydride－based hydrogen storage reactor［J］.Journal of Alloys and Compounds，2009，472（1）：466－472.