氢的安全高效储存是氢能应用的关键。金属氢化物储氢具有储氢密度高、安全性好、氢气纯度高等优点^[1-3]，是一种具有很大发展潜力的车载储氢方式。然而，金属氢化物粉末导热性差，使金属氢化物反应器内部热量传递缓慢，影响了吸放氢反应速率^[4-6]。基于研究成本、安全性与方便施加不同操作条件等方面的考虑，数值模拟是研究储氢反应器性能非常有效的方法。国内外研究者针对金属氢化物储氢反应器的热质传递特性开展了较多的模拟研究工作。

鲍泽威等^[7]建立了金属氢化物反应器的2维轴对称模型，发现降低氢化物床层与换热管壁面之间的接触热阻，提高氢化物床的有效导热系数，可以提高传热效率，加快吸氢反应。Bhouri等^[8]基于数值模拟方法对储氢反应器的吸氢过程进行了优化研究，所研究反应器内混合装填了LaNi_4.3Al_0.4Mn_0.3合金和Li–Mg–N–H储氢材料。袁晟毅等^[9]建立了内置螺旋换热管的金属氢化物反应器的3维数学模型，分析了不同操作参数下反应器的传热性能及吸氢反应特性。Chippar等^[10]提出了一种基于分室化方法的层叠式金属氢化物反应器，通过数值仿真发现反应器的分室化使得反应床层内温度分布更均匀并可以提高吸氢反应速率。尧兢等^[11]建立了相变材料蓄热的金属氢化物储氢反应器模型，模拟研究了操作参数及相变材料物性参数对反应器工作过程的影响。Chibani等^[12]建立了一个2维数学模型来研究金属氢化物反应器内的吸氢反应过程，发现反应速率常数、活化能、床层导热系数和供氢压力对床层温度和吸氢反应的影响较大。Lewis等^[13]提出在金属氢化物反应器内部集成压花板式换热器（embossed plate heat exchanger，EPHX），并通过数值模拟方法研究了此反应器的传热性能及吸氢反应特性。Chandra等^[14]设计了一台内置锥形翅片和换热管的储氢反应器，其内部可以装填5 kg的LaNi₅合金，通过数值仿真发现锥形翅片具有更大的表面积，而且锥形翅片具有漏斗效应可以更有效地装填储氢合金粉末，从而具有更好的强化换热效果。

已有研究者针对金属氢化物反应器开展的一些模拟研究工作，主要针对吸氢过程，对放氢过程则关注较少。金属氢化物储氢反应器的放氢特性对氢燃料电池的运行会有较大的影响。一些研究者基于数值模拟方法分析了换热流体温度、出口压力、对流换热系数等参数对金属氢化物反应器放氢过程的影响^[15-17]。Chung等^[18]模拟了镁基金属氢化物储氢罐的放氢过程。然而，以上研究以LaNi₅和Mg₂Ni合金为储氢材料，两种合金分别存在储氢量低、成本高和放氢温度高等问题，限制了其商业应用。而钛系储氢合金在车载储氢领域具有较大优势，同时也属于中国的优势特色资源。但关于钛系金属氢化物储氢反应器放氢过程中床层热质传递特性的模拟研究还未见报道。

作者建立了装填Ti_0.95Zr_0.05Mn_1.55V_0.45Fe_0.09储氢合金的金属氢化物反应器放氢过程的多物理场耦合模型，分析了不同换热流体温度、流速和氢气排出压力对氢化物床层热质传递特性的影响，并比较了优化操作参数和基准操作参数两种操作条件下反应器的放氢性能。

1 反应器数学模型

研究的圆柱形反应器内部填充有Ti_0.95Zr_0.05Mn_1.55V_0.45Fe_0.09合金与膨胀石墨组成的复合压块，如图1（a）所示。其中，复合压块中膨胀石墨的质量含量为5%。此储氢合金具有活化容易、吸放氢反应速率快、平衡压力适中等优点，且能在常温下实现放氢^[19-20]。此合金与膨胀石墨制作成的复合压块在经过1000次吸放氢反应循环后依然可保证高的储氢量和有效导热系数。所建立的反应器模型为二维轴对称模型，求解区域包含复合压块床层、换热管和气体自由空间，如图1所示。为建立描述放氢过程模型，做出如下假设：1）金属氢化物和氢气的热物性恒定；2）复合压块床层为多孔介质，局部热平衡是有效的；3）忽略辐射传热^[18]；4）反应器与外界绝热，仅与换热流体发生热交换。

1.1 复合压块床层

Ti_0.95Zr_0.05Mn_1.55V_0.45Fe_0.09合金与氢气的反应动力学方程为^[21]：

$ \frac{{\partial {\omega _{\rm{H}}}}}{{\partial t}} = - {K_0}\exp \left( { - \frac{{{E_{\rm{A}}}}}{{{k_{\rm{B}}}T}}} \right)\left( {\frac{{p - {p_{\rm{{eq,d}}}}}}{{{p_{\rm{{eq,d}}}}}}} \right)\left( {{\omega _{{\rm{H}},\max }} - {\omega _{\rm{H}}}} \right) $

(1)

式中，ω_H为储氢量，t为时间，K₀为反应常数，E_A为活化能，k_B为波尔兹曼常数，T为床层温度，p为氢气压力，p_eq,d为脱氢反应平衡压力，ω_H,max为最大储氢量。

可通过下式计算p_eq,d：

$\frac{{{p_{\rm{{eq,d}}}}}}{{{p_0}}} = \exp \left( { - \frac{{\Delta {H_{\rm{d}}}}}{{RT}} + \frac{{\Delta {S_{\rm{d}}}}}{R}} \right)$

(2)

式中，p₀为参考压力，ΔH_d为放氢反应焓，R为通用气体常数，ΔS_d为放氢反应熵。

床层能量方程为：

$\begin{aligned}[b] &\left[ {\left( {1 - \varepsilon } \right){c_{{\rm{p}},{\rm{s}}}}{\rho _{\rm{s}}} + \varepsilon {c_{\rm{{p,g}}}}{\rho _{\rm{g}}}} \right]\frac{{\partial T}}{{\partial t}} + \nabla \cdot \left( { - {\lambda _{\rm{s}}}\nabla T} \right)= \\ & - \left( {1 - \varepsilon } \right){\rho _{\rm{s}}}{\alpha _{\rm{M}}}\frac{{\partial {\omega _{\rm{H}}}}}{{\partial t}}\frac{{\Delta {H_{\rm{d}}}}}{{{M_{\rm{g}}}}} \end{aligned} $

(3)

式中，c_p,s为床层的比热容，ε为孔隙率，ρ_s为床层的密度，c_p,g为氢气比热容，ρ_g为氢气的密度，λ_s为床层的有效导热系数，α_m为复合压块中储氢合金的质量分数。

1.2 换热管

换热管壁的能量方程为：

${\rho _{\rm{w}}}{c_{{\rm{p}},{\rm{w}}}}\frac{{\partial {T_{\rm{w}}}}}{{\partial t}} = \nabla \cdot \left( {{\lambda _{\rm{w}}}\nabla {T_{\rm{w}}}} \right)$

(4)

式中，ρ_w为换热管壁的密度，c_p,w为换热管壁的定压比热容，T_w为换热管温度，λ_w为换热管壁的导热系数。

1.3 气体自由空间

气体自由空间中氢气流速较低，所占体积很小，忽略对流换热作用，其能量方程为：

${\rho _{\rm{g}}}{c_{{\rm{p,g}}}}\frac{{\partial {T_{\rm{g}}}}}{{\partial t}} = \nabla \cdot \left( {{\lambda _{\rm{g}}}\nabla {T_{\rm{g}}}} \right)$

(5)

式中，ρ_g为氢气密度，c_p,g为氢气的定压比热容，T_g为氢气温度，λ_g为氢气的导热系数。

1.4 初始与边界条件 1.4.1 初始条件

假定氢化物初始反应分数和初始温度分布均匀。初始反应分数设为1.0，初始温度等于换热流体温度（T_f）。

1.4.2 边界条件

绝热壁面：

$\frac{{\partial T}}{{\partial {{n}}}} = \frac{{\partial {T_{\rm{w}}}}}{{\partial {{n}}}} = \frac{{\partial {T_{\rm{g}}}}}{{\partial {{n}}}} = 0$

(6)

式中，n为外法线矢量。

氢气出口边界：

$\frac{{\partial {T_{\rm{g}}}}}{{\partial {{n}}}} = 0$

(7)

换热管与换热流体之间换热壁面：

${\lambda _{\rm{w}}}\frac{{\partial {T_{\rm{w}}}}}{{\partial {{n}}}} = h\left( {{T_{\rm{w}}} - {T_{\rm{f}}}} \right)$

(8)

式中，h为对流换热系数。

1.5 换热系数

换热管与换热流体之间的换热系数h由下式计算：

$h{\rm{ }} = \frac{{Nu \cdot \lambda }}{d}$

(9)

其中，努赛尔数Nu为：

$Nu = \frac{{\left( {f/8} \right) \cdot Re \cdot Pr}}{{1.07 + 12.7{{\left( {f/8} \right)}^{\frac{1}{2}}}\left( {P{r^{\frac{2}{3}}} - 1} \right)}} \cdot {\left( {\frac{{{\mu _{\rm{b}}}}}{{{\mu _{\rm{w}}}}}} \right)^n}$

(10)

式中，Re为雷诺数，Pr为普朗特数，f为摩擦系数， $\;{\mu _{\rm{b}}} $ 为管中心处的流体黏度， $\;{\mu _{\rm{w}}} $ 为管壁处的流体黏度。其中，f由下式求得：

$ f=(1.82\mathrm{lg}\;Re-1.64{)}^{-2}$

(11)

1.6 模型求解与验证

金属氢化物反应器模型采用COMSOL Multiphysics 5.4软件求解。模型选用三角形非结构性网格，经过网格独立性验证后确定的网格单元数为71882。求解器选择直接求解器MUMPS，相对容差限定为0.0001，所用参数如表1所示^[21]。

为了对所建立的模型进行验证，将计算结果与Bürger等^[19]与Capurso等^[20]的实验结果进行对比。在Bürger的实验中反应器内填充的材料为Ti_0.95Zr_0.05Mn_1.55V_0.45Fe_0.09合金与膨胀石墨，其中膨胀石墨的质量含量为5%。选择对比的参数为床层的温差变化和放氢质量，对比结果如图2所示。在Capurso的实验中反应器内填充的材料为Ti_0.95Zr_0.05Mn_1.55V_0.45Fe_0.09合金颗粒。选择对比的参数为储氢罐平均温度，对比结果如图2所示。从图2中可以发现数值结果与实验结果存在微小差异，但仍能准确地预测放氢反应完成所需的时间，模型的有效性得到了验证。

2 结果与讨论 2.1 床层温度与反应分数分布

本文中定义的基准操作参数为：流体温度为323.15 K、流体平均流速为3 m/s，氢气排出压力为0.5 MPa。在研究中，除非特别注明，否则按以上参数设置。在基准操作参数条件下，反应床层的温度和平均反应分数的分布如图3所示。图3（a）为床层温度分布随时间的变化情况。由于放氢过程为吸热反应，反应床层的温度在初始阶段迅速降低，最低温度值达到288 K以下。此外，可以发现换热管附近床层的温度高于其余区域。在时间为9000 s以内，温度分布随时间变化并不明显，这说明在该时间段反应器内还在持续进行着放氢反应，致使反应器温度没有明显升高。而在12000 s时反应器的温度快速提高，说明此时吸热反应几乎停止，使得床层在换热流体的加热下温度快速升高。

图3（b）为床层反应分数分布随时间的变化情况。反应分数在换热管壁附近较低，因为此处的换热效果更好，导致反应速率更快。在200 s时，换热管壁附近的放氢反应基本全部完成，但其他区域的反应分数仍然很低。从图3（b）中可以观察到反应锋面，随着反应的进行，反应锋面从换热管向反应器壁面逐渐移动。

2.2 操作参数的影响

为了保证车载燃料电池系统正常运行，根据美国能源部设定的目标，储氢系统的最低输出压力为0.3 MPa。图4给出了在氢气排出压力分别为0.3、0.5、0.7和0.9 MPa时床层平均反应分数的瞬态变化情况。如图4所示，排出压力越高，所需要的反应时间越长。当氢气压力从0.3 MPa增加到0.9 MPa时，放氢反应时间从7740 s增加到20200 s。这是由于放氢反应速率与氢气压力成反比关系，排出压力的增加抑制了放氢过程。

燃料电池的废热可以用于加热金属氢化物储氢罐。为了研究换热流体温度（T_f）对平均反应分数的影响，根据燃料电池运行温度设定T_f分别为313.15、323.15、333.15、343.15和353.15 K。图5为不同流体温度下平均反应分数随时间变化图，可看到随着T_f的升高，放氢反应速率越来越大。当T_f从313.15 K提高到353.15 K，反应时间缩短了约61%。这是因为较高的T_f会导致较高的平衡压力，从而加速放氢反应；又因为T_f越高，换热管温度与床层温差越大，传热更加迅速，使反应速率提高。

随着换热流体平均速度的增加，换热流体在换热管内的对流换热系数h增大，导致传热流体与氢化物床层之间的换热速率增大。根据换热流体的常用流速范围，确定换热流体平均流速为1～7 m/s。当流体平均流速为1、3、5、7 m/s时，其相应h分别为6770.11、16547.58、25255.70、33434.34 W/（m²·K）。因此，增大流体平均流速提高了对流传热系数，加速了氢气的释放过程，这种现象在图6中可以观察到。当氢气流速从1 m/s增加到3 m/s时，反应时间从11200 s降低到了10740 s。但当流体流速超过3 m/s时，由于氢化物床与换热管壁之间的接触热阻成为整个传热过程的主要热阻，增加流速的效果已不明显。而且当流体平均流速高于5 m/s时，会导致流动阻力过大。因此，流体平均流速取为3 m/s就可使所研究反应器具有较高的放氢速率。

2.3 操作参数优化

由以上分析可知，优化的操作参数为：氢气排出压力为0.3 MPa、换热流体温度为353.15 K和换热流体平均流速为3 m/s。在此优化的操作条件下，计算得到h=19081.11 W/（m²·K）。在该条件下进行模拟仿真并与基准操作参数下的模拟结果进行对比，结果如图7所示。从图7发现，当放氢反应进行到3000 s时，优化操作参数下的床层平均反应分数为26.37%，而在基准操作参数下的床层平均反应分数为57.81%。

美国能源部设定的目标为：到2020年，车载储氢装置的重量储氢量达到4.5%、体积储氢量达到30 g /L、平均放氢速率达到0.08 g/s。本文所设计反应器的储氢量约为263 g，重量储氢量和体积储氢量分别为0.93%和50.37 g /L。在优化操作参数下，反应器释放出90%的储氢量所耗时为4060 s，平均放氢速率为0.06 g/s。重量储氢量未达到要求的原因是该材料本身的重量储氢量较低，仅仅只有1.5%。然而，目前实用的储氢合金的储氢量都较低，远远满足不了目标要求。此外，想要提高储氢装置的平均放氢速率可以采用多个反应器单元并联使用的方案。

在优化操作参数下，在2维模型中选取高度为300 mm，半径分别为15、25、35和45 mm等4个位置，读取了其反应分数与床层温度，结果如图8所示。

图8（a）给出了不同位置处的床层反应分数随时间变化趋势。从图8（a）可以发现，越接近换热管处放氢反应进行得越快。当半径分别为15、25、35、45 mm时，反应结束的时间分别为700、1620、3340、4780 s。从图8（b）可以发现，反应开始后，各位置温度迅速下降到275.7 K。此温度所对应的平衡压力与排气压力是近似相等的。之后被换热流体逐渐加热，慢慢升高至353.15 K（换热流体温度）。

图9为点半径15 cm，高度300 cm优化操作的结果。从图9可以发现，在反应初始阶段（约0～4 s），由于平衡压力与排气压力之间的差值很大，导致放氢反应速率非常高，相应的床层温度变化率也很大。在此阶段，换热流体所携带的热量还来不及被床层吸收，换热流体的温度对放氢反应过程没有太大影响。但随着床层温度逐渐降低到275.7 K时，放氢反应会停止。从图8（b）可以发现，接近换热管处床层首先被加热，之后向反应器壁面处逐渐被加热。相应地，放氢反应也从换热管处向壁面处逐渐完成，图3（a）中观察到的反应锋面在优化操作条件下也会出现。在此阶段，传热过程对放氢反应有着较大的影响。

3 结　论

建立了金属氢化物储氢反应器的数学模型，通过与实验结果对比，验证了模型的有效性，得到了在不同工况下反应器的热质传递特性。主要结论如下：

1）当氢气压力从0.3 MPa增加到0.9 MPa时，放氢反应时间从7740 s增加到20200 s；当换热流体温度从313.15 K升高到353.15 K时，放氢时间缩短了约61%。说明降低氢气排出压力，提高传热流体入口温度均能够提高反应速率，缩短反应器放氢反应的时间。

2）基准参数下的放氢时间为10740 s，优化参数下的放氢时间为4720 s，反应时间缩短了约56%，对操作参数的优化能够大大提高反应器的放氢速率。

3）放氢反应的初始阶段主要受氢气压力影响，但仅仅在反应进行的前4 s属于这个阶段。而随后的反应过程主要是受传热过程控制。