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

Phl_zovnf

3690人浏览 · 2024-11-04 23:15:28

Phl_zovnf · 2024-11-04 23:15:28 发布

参考文献：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$ ，令 $y=-ln(1-\varepsilon )^{\frac{1}{n}}$ ， $x=ln(\frac{P_{g}}{P_{eq}})\cdot t$ ，因此原方程 $-ln(1-\varepsilon )^{\frac{1}{n}}=Ca\cdot exp(\frac{-Ea}{RT})ln(\frac{P_{g}}{P_{eq}})t$ 可看为： $y^{\frac{1}{n}}=k\cdot x$ ，即 $y=(kx)^{n}$ ，假设 $n=1$ ，有 $y=kx$ ，因此拟合 $-ln(1-\varepsilon )$ 与 $ln(\frac{P_{g}}{P_{eq}})\cdot t$ 关系时，是一条过原点（0，0）的直线，但是实际用 $y=kx$ 拟合时，有可能直线不那么直，因此反应级数 $n\neq 1$ ，直接在Origin软件里面用 $y=(kx)^{n}$ 函数去拟合数据点，获得一系列级数 $n$ ，然后取平均值 $\bar{n}$ ，如文献Fig.5那样的操作。同时，我们也得到一系列斜率 $k$ 的值，由于 $k=C_{a}exp(\frac{-E_{a}}{RT})$ 是在不同温度下的值，该等式两边取对数，有 $lnk=lnC_{a}-\frac{E_{a}}{R}\cdot \frac{1}{T}$ ，也就是说可以做线性 $y=kx+b$ 拟合，由于 $E_{a}$ 一般单位是 $kJ/mol$ ， $kJ=1000J$ ，因此很多时候就把这 $1000$ 乘在 $\frac{1}{T}$ ，就变成了 $x$ 轴是 $\frac{1000}{T}$ ， $y$ 轴是 $lnk$ ，这样利用拟合得到的斜率和截距，就可以还原出吸氢激活能 $E_{a}$ 的值吸氢反应常数 $C_{a}$ 的值。

利用取值软件webplotdigitizer，获取文献中图5（a）数据，重新在Desmos中绘画，得到下以下散点图，并通过公式 $y=(kx)^{n}$ 去拟合数据点，得到下图带颜色的拟合曲线，可以看到拟合效果跟文献的基本是差不多的

利用取值软件webplotdigitizer，获取文献中图5（b）数据，重新在Desmos中绘画，得到下图中的红色散点图，通过线性拟合，得到拟合公式 $y=-0.892852x+0.716697$ ，跟文献的 $y =-0.89629*x + 0.73007$ ，非常接近。不过，通过上一步得到的一系列k的值，作出lnk与1000/T的散点图，即下图中的黑色点，然后再进行线性拟合，得到拟合公式 $y=-0.829709x+0.485381$ ，虽然看直线挺接近，单斜率和截取就跟文献差得比较多了。

以下是取值软件得到的数据以及根据公式得到的结果数据，与文献图片反推数据进行对比。

通过取值软件webplotdigitizer，提取以下文献图4吸氢动力学的点，

得到关于 $\varepsilon \sim t$ 的散点图，再作出 $-ln(1-\varepsilon )$ 与 $ln(\frac{P_{g}}{P_{eq}})\cdot t$ 的散点图，如下图十字“+”的那些数据点，并通过公式 $y=(kx)^{n}$ 去拟合数据点，得到下图不同颜色的虚线拟合曲线。这里，为了使得虚线接近文献的拟合曲线，每条虚线的 $P_{eq}$ 都进行了人为的调整。

选取的 $P_{eq}$ 的值在下表，得到新拟合的k和n的值，同时也做出lnk与1000/T的关系图，得到线性拟合公式 $y=-0.642934x-0.194409$ ，可以看到截距是负数，显然不太合理。另外，下表中“文献Peq公式（atm）”是根据文献图3的PCT曲线拟合公式 $lnP_{eq}=exp(-1889.6928/T+10.91058)$ 计算后得到不同温度下的 $P_{eq}$ 。

如果 $P_{eq}$ 的值就按文献Peq公式（atm）的值呢，得到以下拟合曲线

得到以下表格的值k(2nd)和n(2nd)，可见n的平均值是差不多的，

但是，作出lnk与1000/T的散点图，得到线性拟合公式 $y=-0.471908x-0.965063$ ，显然也不太好。

再利用取值软件webplotdigitizer，如下4个图获取文章动力学图的数值点，经过拟合，得到文献 $P_{eq}$ 的值得选取，分别为273.15K:6.18676MPa；263.15K:4.86797MPa；253.15K:3.54283MPa；243.15K:2.62342MPa。

再反推文献 $P_{eq}$ 的取值接近以下图中位置，约x=1.38wt%，与y轴平行的直线上。

然后导入到Desmos数学画图软件， $P_{eq}$ 取值改变（273.15K:6.18676MPa；263.15K:4.86797MPa；253.15K:3.54283MPa；243.15K:2.62342MPa。），得到下图虚线拟合曲线，得到n平均值n(3rd)=0.893。

进一步得到lnk拟合值和1000/T的线性拟合， $y=-0.749118x+0.242895$ ，合理些。

由此计算得到 $Ca=1.274934\, 1/s, \, Ea=6.228167kJ/mol$ ；而根据文献原式子 $y=-0.89629*x + 0.73007$ ，可推导出 $Ca=2.075225\, 1/s, \, Ea=7.452175kJ/mol$ 。可以看见，算出来的 $Ca$ 和 $Ea$ 跟文献还是有点差距的。

由此得到吸氢动力学方程： $-ln(1-\varepsilon )^{\frac{1}{0.894}}=1.275\cdot exp(\frac{-6228}{RT})ln(\frac{P_{g}}{P_{eq}})t$ ，仿真的公式就如下面所示。

$\frac{\partial c}{\partial x}=1.275\cdot exp(\frac{-6228}{RT})ln(\frac{P_{g}}{P_{eq}})n(1-\frac{c}{c_{max}})\left [ -ln(1-\frac{c}{c_{max}})^{\frac{n-1}{n}}\right ]$