关于金属氢化物（储氢）PCT曲线拟合、ZBS有效导热系数模型、JMAK类型吸放氢动力学方程的笔记

参考文献：Experimental and numerical study of metal hydride beds with Ti0.92Zr0.10Cr1.0Mn0.6Fe0.4 alloy for hydrogen compressionhttps://www.sciencedirect.com/science/article/pii/S1385894723043851?via%3Dihub#s0010

一、PCT曲线拟合

根据以下文献内容，可知PCT曲线是用多项式来拟合的，同时，该拟合的PCT曲线方程是在某个参考的温度下实测的曲线，然后用多项式去拟合这个实测的PCT曲线。另外，如果根据实测PCT曲线中截取不同温度下平台中点 $P_{eq}$ ，制作 $lnP_{eq}\sim 1/T$ 曲线，然后用 $lnP_{eq}=\frac{-\Delta H}{R}\cdot \frac{1}{T}+\frac{\Delta S}{R}$ 或 $lnP_{eq}=\alpha +\beta \cdot \frac{1}{T}+\gamma \cdot \frac{1}{T^{2}}$ 去拟合。

根据文献，吸氢 $lnP_{eq}=\frac{-\Delta H}{R}\cdot \frac{1}{T}+\frac{\Delta S}{R}$ 直线拟合的比较好， $y=-1889.6928x+10.91058\Rightarrow \frac{-\Delta H}{R}=-1889.6928, \frac{\Delta S}{R}=10.91058$ ，解得 $\Delta H\approx 15.711 kJ/mol, \Delta S\approx 90.71 J/(mol\cdot K)$ 。放氢 $lnP_{eq}=\alpha +\beta \cdot \frac{1}{T}+\gamma \cdot \frac{1}{T^{2}}$ 拟合得更好，有 $y=7.91147+524.4724x-481225.91276x^2$ 。

将 $lnP_{eq}=\frac{-\Delta H}{R}\cdot \frac{1}{T}+\frac{\Delta S}{R}$ 和 $lnP_{eq}=\alpha +\beta \cdot \frac{1}{T}+\gamma \cdot \frac{1}{T^{2}}$ 两式进行对比， $\Delta H=-R\cdot (\beta -\frac{2\gamma }{T}), \Delta S=\alpha +\frac{\gamma }{T^2}$ ，代入 $\alpha ,\beta ,\gamma$ 的值，就得到以下关系：

根据文献中拟合的多项式数据，有如下曲线，蓝色是放氢PCT拟合，红色是吸氢PCT拟合。注意，公式是 $y=a_{1}\cdot x^1+a_{2}\cdot x^2+...+a_{12}\cdot x^{12}$ ，没有 $a_{0}$ 项，且注意y的单位应该是atm，不是MPa。

根据 $P_{eq}=P_{eq,ref}\cdot exp[\beta \cdot (\frac{1}{T}-\frac{1}{T_{ref}})+\gamma\cdot (\frac{1}{T^{2}-}\frac{1}{T_{ref}^{2}})]$ ，其中 $P_{eq,ref}$ 是上述拟合的曲线函数，可得到不同温度下的放氢PCT拟合曲线（如下橙色线和绿色线）

吸氢就按 $P_{eq}=P_{eq,ref}\cdot exp[\frac{-\Delta H}{R} \cdot (\frac{1}{T}-\frac{1}{T_{ref}})]$ 公式，得到不同温度下的吸氢PCT拟合曲线（如下紫色线和黑色线）。

二、ZBS有效导热系数模型

以下为ZBS模型公式，其中 $b=9.87$ ， $\lambda _{e}$ 是有效导热系数， $\lambda _{g}$ 是氢气导热系数， $\delta$ 是孔隙率，B是形状因子， $\lambda _{g,ref}$ 是氢气在0.1MPa下的参考导热系数， $l_{m}$ 是氢气的气体自由程， $r_{p}$ 是孔的平均大小，在这里笔者认为是颗粒与颗粒直接的空隙大小。

参考文献:Modified Zehner-Schlunder models for stagnant thermal conductivity of porous media(Int. J. Heat Mass Transfer. Vol. 37, No. 17, pp. 2751-2759, 1994)

按参考文献，考虑以上几何形式的热量通道，具体的推导过程不太会，这里提到原文里的结论：B是引入的形状因子，当 $B=0$ 时， $r=0$ ，代表没有固体颗粒域； $B=1$ ， $r^2+z^2=1$ ，说明底面是圆，绕z轴旋转， $A_{fs}$ 就是一个球面，固体颗粒域就是个球；当 $B=\infty$ 时， $r=1$ ，说明底面是矩形，绕z轴旋转， $A_{fs}$ 就是一个圆柱面。

原作者根据推导的公式（8），将B与孔隙率Φ用近似公式（9b）来替代公式（8），该文献作者则修正了近似公式的系数和指数，C=1.364，m=1.055

好了，到此，已知的量有孔隙率 $\phi$ ，推出形状因子B， $\lambda _{g}$ 还需要知道 $l_{m}$ 和 $r_{p}$ 的值，可惜第一篇参考文献没有给出这些参数具体的值，笔者在这里根据网上查的值和文献里面给的其他参数值进行推算，得到 $\lambda _{e}\approx 3.46 W/m/K$ ，感觉跟实际范围差不多。

三、JMAK型吸放氢动力学方程

关于JMAK方程怎么推导，可见：

储氢合金/金属氢化物CFD传质公式来源（仅针对JMAK方程）_jmak有效动力学参数-CSDN博客https://blog.csdn.net/qq_24800941/article/details/129191569为什么选JMAK方程，因为其拟合效果不错。

吸氢：

当温度 $T$ 恒定， $C_{a}exp(\frac{-E_{a}}{RT})$ 可视为常数 $k$ ，因此原 $JMAK$ 方程可看为： $y^{\frac{1}{n}}=k\cdot x,\frac{1}{n}\cdot lny=lnk+lnx, y=-ln(1-\varepsilon ),x=ln(\frac{P_{g}}{P_{eq}})\cdot t$ ，假设 $n=1$ ，有 $y=kx$ ，因此拟合 $-ln(1-\varepsilon )$ 与 $ln(\frac{P_{g}}{P_{eq}})\cdot t$ 关系时，是一条过原点（0，0）的直线，但是实际用 $y=kx$ 拟合时，有可能直线不过原点（0,0），因此反应级数 $n\neq 1$ ，即 $y^{\frac{1}{n}}=k\cdot x,\frac{1}{n}\cdot lny=lnk+lnx, y=-ln(1-\varepsilon ),x=ln(\frac{P_{g}}{P_{eq}})\cdot t$ ，需要对y轴 $-ln(1-\varepsilon )$ 和x轴 $ln(\frac{P_{g}}{P_{eq}})\cdot t$ 再进行对数，获得一系列级数 $n$ ，然后取平均值 $\bar{n}$ 。

然后固定时间 $t$ ，到达 $t$ 时刻时，不同温度 $T$ 下，有不同的反应分数 $\varepsilon$ ，将 $C_{a}ln(\frac{P_{g}}{P_{eq}})\cdot t$ 看成系数 $p$ ， $-ln(1-\varepsilon_{t} )^{\frac{1}{n}}$ 看成 $k$ ，所以有 $k=p\cdot exp(\frac{-E_{a}}{RT})$ ，两边取对数，得 $lnk=lnp-\frac{1}{RT}\cdot E_{a}$ ，所以作出 $lnk$ 与 $1000/T$ 的拟合直线，根据斜率得到吸氢激活能 $E_{a}$ 的值，根据截距得到吸氢反应常数 $C_{a}$ 。