储氢合金/金属氢化物床层有效导热系数的数学模型

最近看到一篇有关“储氢合金/金属氢化物床层有效导热系数的数学模型”的论文，文章DOI：10.1016/j.energy.2023.127085，文章提到的数学物理模型还算好理解一些，文章作者是西安交通大学杨福胜教授团队的学生，笔者也跟该文作者进行了联系，加深了公式的理解，在此特意分享给各位感兴趣的大佬。

 一、物理模型简图和假设

         文章里，对储氢合金床层作了以下几个假设：

1. 粒子是均匀的球体，粒子床是简单的立方体堆叠；

2. 颗粒处于点接触状态；

3. 忽略了通过颗粒接触的热传导，因为颗粒床的有效导热系数在真空中极低；

4. 颗粒床中的传热处于稳定状态。

    下面是设想的特征单元体，也就是本次文章储氢合金床层的数学模型，其特点就是一个边长为颗粒直径D的立方体中，有两个半径为R的半球形，D=L=2R。通过截面图分析，导热路径经过固体颗粒的路径长度为Ls，经过气体路径长为Lg，Lg+Ls=D=2R=L。

  二、一些热流量等物理量的数学推导

        2.1 几何数学关系

    根据截面图，通过几何数学关系，得到以下关系：   

    由于模型相当于是等效气膜厚度+等效颗粒床层厚度构建的物理模型，根据传热学中热阻串联原理，根据几何关系等效气膜厚度Lg，等效颗粒床层厚度为Ls，根据热流密度定义式，得，可解得:。

进一步根据数学关系，可得:

。

        2.2 导热量的计算分析：

      再根据傅立叶定律，得到以下等式：

     通过对dQ从角度0到σ（0≤θ≤π/2）进行积分，积分过程如下，得到文章里面的等式 

​​​​

 

    特别地，当σ=π/2时，Q的积分结果如下所示：

    分析纯气体的导热路径，导热面积S上相当于正方形面积减去圆的面积，易得以下等式 ：

     特征单元体中总的热流量是以下两者的和Qall，易得下列公式：

    并定义η1、η2，η1的物理意义是“颗粒-气膜-颗粒”导热路径中导热量与整个特征单元的总导热量的占比，η2的物理意义为角度0到σ积分的“颗粒-气膜-颗粒”导热路径中导热量与整个特征单元的总导热量的占比，η1、η2公式如下所示，这两个公式为超越方程：

        定义η3为有效气膜比，公式如下：

        注意η3中，是指η2=0.9时的有效气体的圆柱体体积/特征单元体总气体体积。

         2.3 η2公式的简化，以及η3如何计算的问题

         根据文献，η2的公式简化、推导过程和最终表达式如下：

    η2=0.9时，求出σ的值，得到对应Vgas,eff的值，再用η3的定义式求出η3的值。对于Vgas,eff的几何体积求法如下，是圆柱体减去球缺的体积，Vgas,all是正方体减去球型的体积。  

         2.4 η1、η3的计算结果以及它们的物理意义 

    通过代人具体的λs和λg的值，可解出η1的值，在Excel中输入公式，便可获得与文献表格中的η1数据。同时可以看到，η1均大于0.95，由于λs一般都大于1W/(m·K)，因此，可以说是等效传热路径中，纯气路传热只占<5%，而95%以上均是“颗粒—气膜——颗粒”路径。 

~文献结果：

~Excel计算结果：

    由于本人不太会解lnx-x=C的方程，因此借助数学工具desmos，当η2=0.9,、λs=1，λg=0.0258时，解得1-cosθ=0.5563，(sinθ)^2=1-(1-(1-cosθ))^2=0.80313031。因此得η3=0.321，与文献里表格1的η3=0.306接近，但不相近。我们在算一下 λs=5，λg=0.0258时，η3=2.00，文献η3=1.98；  λs=10，λg=0.0258时，η3=0.166，文献η3=1.65；  λs=20，λg=0.0258时，η3=0.139，文献η3=1.41。

三、等效导热系数的计算 

     根据上面这个图的分析，有三种导热路径，情况(b)是通过颗粒之间接触的纯颗粒之间导热，情况(a)是“颗粒-气膜-颗粒”路径，根据之前分析，这种情况占据主要地位(>95%)，最后的情况(c)是纯气体路径。

    热流密度q=λΔT/Δx，假设床层孔隙率为ε，颗粒之间接触比为δ，α为床单元中纯气体占总气体的比例。当接触比为1时，等效导热系数为λs，由于接触比为δ，因此路径体积比为λs*δ，由于孔隙率为ε，即有ε的比例非纯颗粒之间传递热量，又由于ε的比例的体积中，有α比例为纯气体，因此折算有αε为真正纯气体的导热比例，颗粒之间为(1-αε)的比例，所以总的路径体积比为δ*(1-αε)，所以情况(b)：q_b=λs*δ*(1-αε)*ΔT/L=λs(1-αε)δΔT/L。同理，情况(c)：q_c=λg*αε*ΔT/L=λgαεΔT/L。

    情况(a)中，根据之前分析，，因为非接触比为(1-δ)，颗粒之间为(1-αε)的比例，所以该情况总的路径体积比为(1-αε)*(1-δ)，所以情况(a)：q_a=λe*(1-δ)*(1-αε)*ΔT/L=(1-αε)(1-δ)/(Ls/L/λs+Lg/L/λg)ΔT/L。

     总热流密度q=q_a+q_b+q_c，得到以下关系：

    根据之前分析，q_a占95%以上，纯气体路径还是要归并到“颗粒-气膜-颗粒”中，床层孔隙率为ε，颗粒之间接触比为δ，(1-ε)δ表示的是非孔隙部分“颗粒-颗粒”路径的占比，由于纯气体路径归并，因此是整体1减去(1-ε)δ的比例=1-(1-ε)δ，这些都是“颗粒-气膜-颗粒”导热路径的比例，所以路径(a)的总路径体积比为1-(1-ε)δ，真空中，纯靠颗粒之间导热量很低，推断δ≈0，因此还是得到以下公式：

,。

     同时，定义系数N=1/η3，η3的物理意义是有效气膜比，因此等效的有效气膜厚度Lg'=Lg*η3=Lg/N，同时孔隙率ε应该理解为特征单元体中，总气膜厚度Lg占整个特征单元体长度L的比值，也就是说ε=Lg/L，Lg'=L*(ε/N)，又因为Ls+Lg'=L，所以有1-ε=Ls/L，Ls=L*(1-ε)。根据上面分析，整个特征单元长度L=L'=Lg'(等效的有效气膜厚度)+Ls，因此根据公式，显然有：

        N=1/η3,文献作者通过拟合数据（类似表格1里面的数据，x轴为k=λs/λg，y轴为N=1/η3）得到以下关系式。

因此只需知道孔隙率ε，颗粒导热系数λs，气体导热系数λg，即可得到等效导热系数λe。

有兴趣的读者可以自行翻看原文，文章连接如下：

Effective thermal conductivity of metal hydride particle bed: Theoretical model and experimental validationhttps://www.sciencedirect.com/science/article/pii/S0360544223004796?via%3Dihub