储氢合金氢气释放过程多物理场耦合建模分析

袁小野

（武汉理工大学汽车工程学院 湖北省现代汽车零部件技术重点实验室 武汉 430070）

金属氢化物储氢是指在一定的环境温度和环境压强的条件下，固态氢、气态氢气、储氢多孔介质合金进行质量交换、动量交换、并伴随着有反应热的能量交换的综合多物理场耦合过程［1］.2002以前，针对上述过程的研究主要以代数计算方式进行研究［2－4］.2002年以来，Delahaye，Hardy等人在 COMSOL Multiphysics 3.0－4.0的平台上进行了二维、二维轴对称的仿真结果，对于金属储氢反映过程的理论研究逐步向可视化、形象化的方向发展，同时也对金属储氢反映过程的偏微分方程的形式进行了进一步的完善［5－6］.然而，考虑到金属储氢吸收和释放反映过程中氢气密度，固态氢密度，容器温度以及压强等参数在容器内部的变化趋势主要受时间因素影响.即，对于某一时刻，物理量空间分布的波动性远小于时间波动性［7］.因此本研究旨在通过二维模型和一维模型的对比分析，验证一维有限元模型在氢气存储效率模型分析中的更广泛的参数分析范围.

1 基本原理

储氢容器的反应过程概括为气态氢气，固态氢气和能量3种物质及其相关物理属性的相互转换过程.其中气态氢气在视为多孔介质的金属储氢容器中满足气态渗流条件和动量守恒定律，其状态满足Darcy方程.

式中：ε为金属氢化物的孔隙率；vg为氢气的达西速度，m／s；ρg为气态氢的空间密度，kg／m3.氢气的渗透过程可以通过方程（2）来描述动量守恒定律，方程（2）中k（m2）为达西渗透率；p为达西压强，Pa；μ（Pa·s）为粘度；物理速度u＝εvg.同时气体在空隙中满足理想气体状态量方程（3）中Mg为气体摩尔质量，kg／mol；Rg为标准气体状态量常数，J／（mol·K）；T 为温度，K.

释放反应过程中，质量源m与金属氢化物的颗粒密度ρs（kg／m3）、相变激发能量Ed（J／mol）、转化率Cd（1／s）、平衡状态压强peq（Pa）、温度场T（K）有关，具体函数关系可以写为

固态氢与气态氢气的转换过程中的质量守恒定律可以通过方程（5）与方程（1）联立描述

反映过程中的能量传递过程符合多孔介质中气体热传递规律.

式中：ρ为混合状态等效密度，kg／m3；Cp为等效热容，J／（kg·K）；keff为等效热导率，W／（m·K）；ΔH 为相变潜热，J／kg.式（6）中左边第二项为对流传导，右边第一项为热传递.

从方程（1）～（6）的联立过程，可以总结出如图1所示的变量耦合关系图.

图1 释放反应过程变量耦合关系图

从图1中不难看出，氢气存储反应过程可以由独立变量描述.从系数项，质量源项对独立变量的函数关系中（式（2）～（4））也可以看出，描述储氢反应过程的偏微分方程组是具有高度非线性程度的.这也对储氢反应过程的高维仿真的收敛性带来了极大的挑战.

2 数值仿真模型

金属储氢释放反应过程根据上述基本原理，采用 COMSOL Multi－physics 5.0 中的 Darcy′s Law模块，Mass diffusion模块和heat transfer in porous media模块建立了二维仿真模型和一维仿真模型.对结果进行了对比分析.同时讨论了收敛性和模型的适用范围.

2.1 几何模型

仿真对空心圆柱形储氢元件进行了二维轴对称和一维建模，分别如图2，图3所示.其中，二维几何模型中，r和z分别表示圆柱坐标系的径向坐标和轴向坐标，L1＝0.75mm为氢气释放边界，L2＝0.57mm，L3，L4和D1分别为元件顶部，外壁，底部和储氢介质反应区域.忽略物理量在反应区域纵向空间分布的变化情况，可取反应区任意水平截线得出一维几何模型如图3所示.图3中p1，p2和L0分别对应储氢元件的氢气释放边界，储氢元件外壁和反应区域.

图2 储氢元件二维轴对称几何模型

图3 储氢元件一维几何模型

2.2 模型参数

图1 中所示释放反应过程中的变量耦合关系的系数见表1.

表1 释放反应过程模拟参数列表

2.3 边界条件

容器上部和侧部边界设置为零通量，同时在求解域中设置质量源m.设置出口边界为V0用于模仿恒定抽气速度的工作状态.在扩散对流PDE模块中，使用初始条件ρs＝ρs0模拟容器满载固态氢时的初始状态.同时，由于固体在容器中只会通过相变方式转换为气体，因此设置所有边界通量为零且对全部求解域设置为负数质量源.能量传递方程中，初始温度为T0，假设储氢设备保温良好，因此设置全部边界热绝缘，且在全部求解域中设置热源项模拟反映放热和加热过程.

3 结果与讨论

根据第三部分的参数设置和边界条件设置，采用自由三角形对二维模型进行网格剖分得到模型空间自由度为29 703的网络模型.对一维模型采用单元数量为500的固定分布网格剖分，得到自由度为2 503的网络模型.采用基于BDF算法的时变求解器求解释放反应从初始状态到稳定状态0～30s的稳定状态的物理场分布，得到如图4～图6所示的结果.

3.1 达西压强对比

二维轴对称和一维仿真模型的达西压强空间分布规律如图4所示.可以看到，储氢元件的最大压强位置在出口处，最小压强在储氢元件的外壁处.这种分布规律是由于设定恒定抽气速度的边界条件所引起.同时，不难看出，达西压强的空间分布波动情况小于10－6.因此，可以近似认为，储氢元件中压强在空间中呈均匀分布.

图4 氢气释放过程达西压强空间分布

由二维轴对称模型和一维模型对比分析可知，二者的最大压强在同一数量级，微小的差别可以归因于一维模型缺失了一个空间维度所致.

进一步地，对氢气释放过程中压强的时间变化规律进行了研究.图5的计算结果展现了释放反应从初始状态到稳定状态0～30s的稳定状态的达西压强时变规律.从图5中可以看出，反应过程开始后，达西压强从反应初始状态的5MPa初始压强，在极短的时间内骤降至稳定状态（约为114kPa）.达到稳定反应后，随着时间的推移，储氢元件中的固态氢不断转化为气态氢.从图5的对数之间坐标图中可以看出，达西压强与随着时间增加呈线性增长趋势.

3.2 收敛性讨论

从图4的对比分析可以看出，二维模型和一维模型的计算结论具有极高的吻合程度.由于物理量的空间分布波动性远远小于时间波动性且在金属储氢元件释放反应过程中工程上更为关注不同时刻的储氢元件的整体性能，因此一维模型具有更高的参数适应度和收敛情况.图6中对比了相同参数设置情况下，二维模型和一维模型在0～50min时间内的收敛情况.从图6中不难看出，对于相同的参数设置，二维的模型的时间步长倒数波动极大，当模拟反应过程时间超过800s后，模型的时间步长倒数快速上升，最终会导致模型容差溢出无法解算.相反，一维模型的收敛图在55s后，时间步长逐渐收敛于恒定值，解算效率是二维模型的数百倍.一维模型的参数适应性更高，更广泛地适用于金属储氢释放反应过程中储氢元件物理场分析.

图5 达西压强－时间曲线

图6 二维模型与一维模型的收敛性对比

4 结束语

金属储氢的存储反应过程是一个伴随吸热过程的气态氢气到固态氢的物理相变过程.反应过程的效率与初始条件以及热源条件等因素密切相关.物理状态的剧烈变化发生在初始阶段，后续过程反应相对平和.需要在容器设计过程中采取特殊措施确保反应过程初始时的安全性.一维模型和二维模型计算结果的时间分布具有一致趋势.以为模型维的收敛性和求解效率远远由于二维模型，适用于更宽松的参数范围.同时由于物理量空间分布的波动性较小，对于储氢反应过程关注物理量时变特性的研究，一维模型具有更高的价值和适用性.

［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.