粉末烧结多孔介质具有比重小、渗透性高、比表面积大以及对亲水性液体吸附能力强等优点^[1]，在多个领域具有广泛应用，如燃料电池^[2]、热管散热^[3]和石油工业^[4]等，研究多孔介质内液体输运具有重要意义。多孔介质内部结构复杂，通常将其简化为等截面毛细管、变截面毛细管、粗糙内表面毛细管或分形毛细管^[5‒7]，但与实际相差较大。球体随机堆积模型可以十分有效地模拟球形粉末烧结多孔介质^[8]，随机几何的出现与发展，使利用球体随机堆积模型研究多孔介质结构特征成为可能^[9]。本文针对球形粉末烧结多孔介质，采用理论分析和计算机模拟，对传统模型进行修正，推导一种准确性更高、适用范围更广的多孔介质结构特征理论模型，同时，将伯努利方程与所建模型耦合，构建多孔介质中液体爬升高度的控制方程，并以多孔介质液体爬升实验为实证对象，计算实验值与控制方程预测值的相关系数，评价控制方程和理论模型的有效性和适用性。

1.   多孔介质理论模型构建

根据多孔介质的定义，将多孔介质分为骨架和孔隙两部分^[10]。球形粉末烧结多孔介质的骨架由许多尺寸相近的球形粉末颗粒烧结而成，根据球体随机堆积理论，可将骨架简化为等直径球体随机堆积模型。基于随机几何理论，等直径球体随机堆积模型的比表面积定义为式（1）^[11]。

式中：S为模型比表面积；R为模型所用球体的半径；ρ_s为单位体积内球体数量，根据粉末烧结多孔介质的定义，${\rho _{\text{s}}} = 6(1 - \phi )/({\text{π }}d_{\text{s}}^{\text{3}})$，其中，ϕ为孔隙率，d_s为平均颗粒直径；rⁿ={r₁, r₂, ∙∙∙, r_n}表示n个球体的球心坐标；h(x‒r)为一个判断函数，当点x为球体内某一点时，h(x‒r)=1；g_n(rⁿ)为n点相关函数，与球体堆积方式有关，常用的球体堆积模型有硬球模型、完全穿透模型和部分穿透模型。

多孔介质内部结构复杂，难以获得准确的n点相关函数表达式。在传统模型中，通常令g_n(rⁿ)=0或g_n(rⁿ)=1，分别对应硬球模型和完全穿透模型^[11]，但适用范围较小，误差较大。为了获得更准确的n点相关函数的表达式，Blum和Stell^[12]建立了部分穿透模型的比表面积，引入不可变形系数（ε），如式（2）所示。

式（2）成立的条件是球体直径分布范围为（0，∞），球体随机堆积模型的配位数为n‒1，即任意两球体均互相接触。本文以粉末烧结多孔介质为研究对象，颗粒直径分布范围较窄，可视为等直径球体，而等直径球体随机堆积模型的最大概率配位数为6^[13]，因此，式（2）并不适用于球形粉末烧结多孔介质。于是，对n点相关函数进行修正，如式（3）所示。

另外，式（1）中存在无限项求和，对n点相关函数进行修正，不仅更适用于球形粉末烧结多孔介质，还可利用指数函数的泰勒级数展开式，将无限项求和化简为有限项求和，降低计算误差。将式（3）带入式（1）并整理，得到修正后的部分穿透模型比表面积，如式（4）所示。

进一步分析可知，当ε=1时，部分穿透模型变成硬球模型；当ε=0时，部分穿透模型变成完全穿透模型，因此，式（4）的适用范围更广。

2.   多孔介质模拟

利用球体随机堆积模型模拟球形粉末烧结多孔介质，并计算多孔介质的孔隙率和比表面积。具体步骤为：（1）利用离散建模软件EDEM在边长为1 mm的正方体空间内模拟球体随机堆积模型，球体直径分别为60、70、80、90和100 μm；（2）对空间连续压缩五次，模拟不同不可变形系数的部分穿透模型；（3）将模拟结果输出为球心在空间直角坐标系中的坐标；（4）利用蒙特卡罗法^[14]计算比表面积，在球体i表面随机生成N个点并计算各点到其他球体球心的距离，距离小于半径的点的数量记为M，则该球体的表面积为S_i=4πR²(N‒M)/N，对所有球体的表面积求和即可得到比表面积；（5）利用蒙特卡罗法计算孔隙率，在正方体空间中随机生成N个点并计算各点到球体球心的距离，距离小于半径的点的数量记为M，则孔隙率为ϕ=(N‒M)/N。不同直径球体堆积模型的模拟结果如图1所示，其中横坐标为孔隙率，纵坐标为比表面积。

图  1  孔隙率和颗粒直径对比表面积的影响

Figure  1.  Effect of porosity and particle diameter on the specific surface area

下载: 全尺寸图片 幻灯片

从图1中可以看出，随压缩次数的增加，孔隙率降低，比表面积先增大后减小。以直径为90 μm的球体堆积模型为例，分别计算每次压缩后，球体堆积模型中任意两球体球心间的距离，若距离小于90 μm，则说明两球体的接触方式为部分穿透。图2为部分穿透球体数量，从图中可以看出，硬球模型中任意两球体球心间的距离均大于90 μm。第一次压缩后，部分穿透球体数量为63对，说明大部分球体间仍为点接触，球体堆积更加密集，模型体积减小使比表面积增大；随压缩次数增加，部分穿透球体数量迅速增加，比表面积减小。

图  2  压缩次数对部分穿透球体数量的影响

Figure  2.  Effect of compression times on the number of the partial penetration spheres

下载: 全尺寸图片 幻灯片

3.   模型模拟验证

利用部分穿透模型预测多孔介质比表面积，需确定多孔介质的不可变形系数。通过上述模拟数据分析发现，不可变形系数与孔隙率有关。图3以孔隙率为横坐标，以不可变形系数为纵坐标，绘制了不可变形系数的模拟结果，虚线为线性拟合方程。从图中可以看出，不可变形系数与孔隙率相关系数在0.98以上，具有较强的线性相关性，因此不可变形系数的预测方程如式（5）所示。

图  3  不可变形系数线性拟合

Figure  3.  Linear fitting of the impenetrability parameter

下载: 全尺寸图片 幻灯片

为了评价多孔介质比表面积理论模型的准确性，将拟合结果带入式（4），利用传统模型和本文所建立的模型预测球体堆积模型的比表面积，并与模拟结果行对比。图4以直径为60 μm的球体堆积模型为例，绘制了模拟结果和各模型预测结果的比表面积。从图中可以看出，随孔隙率减小，硬球模型的误差增大，完全穿透模型的误差降低。当孔隙率约为0.45时，硬球模型的误差约为0.2%，完全穿透模型的误差约为50%；当孔隙率约为0.10时，硬球模型的误差约为100%，完全穿透模型的误差约为20%。部分穿透模型始终具有较高的准确性，最大误差约为6%。因此，与传统模型相比，式（4）的准确性更高，适用范围更广。

图  4  模拟结果和模型预测结果比较

Figure  4.  Comparison of the simulation and model predictions results

下载: 全尺寸图片 幻灯片

4.   模型实验验证

多孔介质液体爬升过程与比表面积有关。当润湿表面的表面能低于干燥表面时，液体具有沿干燥表面扩散的趋势，表面润湿所释放的表面能是多孔介质液体爬升的主要动力来源，多孔介质比表面积越大，单位体积多孔介质润湿所释放的表面能越多，多孔介质液体爬升速度越快。同时，如果将多孔介质孔隙视为一系列微元管道，管内定常流动分析指出，流量与管径的四次方成正比，多孔介质比表面积越大，孔隙尺寸越小，通过多孔介质的流量越小，多孔介质内液体爬升速度越慢。因此，以实际多孔介质液体爬升实验为实证对象，将伯努利方程与所建多孔介质比表面积理论模型耦合，对模型在实际应用中的可靠性进行验证。

首先对多孔介质内液体爬升过程进行假设：（1）液体在多孔介质内的流动是沿长度方向的一维流动；（2）液体为不可压缩的牛顿流体；（3）初始时刻的动能忽略不计，重力势能记为0。根据不可压缩粘性流体的伯努利方程可知，多孔介质内液体爬升高度的控制方程如式（6）所示。

式中：ΔP为液面两侧压差，ρ为液体密度，v和h分别为液体爬升速度和高度，g为重力加速度，h_w为液体爬升过程中的能量损失。

Eisfeld和Schnitzlein^[15]将24种能量损失预测方程与实验数据进行对比，评价出准确性最高、适用范围最广的是Reichelt提出的预测方程，如式（7）所示。

式中：$D$为多孔介质的水力直径，μ为动力粘度。Reichelt指出，采用尺寸相近的球形粉末制备多孔介质，最外一层粉末的孔隙率偏高，流动阻力较小，因此需要考虑多孔介质边界效应对能量损失的影响。

为了获得多孔介质液面两侧压差的表达式，White^[16]以液体爬升平衡位置为研究对象，从能量平衡角度推导出多孔介质的等效管径为r_eff=2ϕ/S，并指出多孔介质液面两侧压差与管径为等效管径的毛细管束模型液面两侧压差相等。因此，由拉普拉斯方程可知液面两侧压差可表示为式（8）所示。

式中：r₁和r₂为两相互垂直截面的曲率半径，γ_lv为气液表面张力，θ₀为平衡接触角。需要注意的是，在平衡位置，液体爬升速度为0，伯努利方程可简化为ΔP=ρgh，液体仅受竖直向上的表面张力和竖直向下的重力作用，受接触角滞后性的影响，平衡接触角小于水平表面的表观接触角，液面两侧压差增大，液体爬升平衡高度增加。

分子动力学理论指出，多孔介质润湿过程是一个动态变化过程，固液表面接触角与单个分子的运动状态和接触线的分布有关^[17]。由平衡位置伯努利方程可知，液体最终爬升高度与平衡接触角的余弦值成正比，因此以平衡接触角为基准。Li等^[17]制备了玻璃珠多孔介质并进行液体爬升实验，通过分析实验数据，给出与液体爬升速度有关的动态接触角（θ）的表达式：cosθ=cosθ₀‒ζv/γ_lv，式中：ζ为材料摩擦系数，与多孔介质的材料和液体性质有关。当液体爬升到平衡位置时，动态接触角与平衡接触角相等，如式（9）所示。

将式（7）和式（9）带入式（6），整理得到关于爬升高度的一阶微分方程，如式（10）所示，其中初始条件如式（11）所示。

式中：h为液体爬升高度，h′为液体爬升高度对时间（t）的导数。将多孔介质比表面积的理论模型代入上式，利用MATLAB编程，使用有限差分法求解不同时刻多孔介质内液体爬升高度。与传统方程相比，该式综合考虑了多孔介质的内部结构、颗粒平均配位数、不可变形系数、边界效应以及动态接触角的影响，较为准确地给出球形粉末烧结多孔介质中液体爬升高度的预测值。

为方便计算，令初始条件为h=10^‒9，t=0，计算过程中其他参数与对应实验相同。以文献[18]为例，文献中的实验参数为d_s=73.8 μm，ϕ=0.3773，ε=0.8332，S=45.41 mm^‒1，γ_lv=75.6 mN·m^‒1，ρ=1.04 g·cm^‒3，μ=1.16 mPa·s，ζ=20 Pa·s，θ=77.3°，D=1.43 mm，g=9.8 m·s^‒2，其中，文献中给出的接触角为水平表面的表观接触角，如果令接触角θ=77.3°，最终液体爬升平衡高度偏小，与理论分析结论相同。因此考虑到接触角滞后性的影响，取θ=58.0°。不同液体爬升实验的实验参数如表1所示。

表  1  不同文献中液体爬升实验参数

Table  1.  Experimental parameters of the actual capillary rise test in the different references

文献	颗粒平均直径， ds / μm	不可变形系数， ε	比表面积， S / mm‒1	孔隙率， ϕ	表面张力， γlv / (mN·m‒1)	动力粘度， μ / (mPa·s)	平衡接触角， θ0 / (°)
[18]	73.8000	0.8332	45.41	0.3773	75.58	1.16	58
[19]	63.5000	0.6266	49.71	0.2996	72.00	0.89	48
[20]	0.0254	0.6277	0.12	0.3000	72.75	0.95	60

下载: 导出CSV 

| 显示表格

为了评价预测结果的准确性和多孔介质比表面积理论模型的适用性，将预测结果与实验数据进行对比，利用Pearson相关系数（r_p）表征实验值和预测值之间的线性相关性，如式（12）所示。

式中：n为数据量，h_ds(i)和h_dm(i)分别是第i个预测值和实验值，${\bar h_{{\text{ds}}}}$和${\bar h_{{\text{dm}}}}$分别为预测值和实验值的平均值，r_p越接近于1，说明预测值的准确性越高。图5为液体爬升高度预测值与相应实验值之间的数据对比和误差分析，其中横坐标为实验测量结果，纵坐标为方程预测结果。由于三个实验最终平衡高度相差较大，为了更清晰地显示不同实验的实验数据，采用双对数坐标轴。图中虚线所示为±20%的误差范围，可以看出，大部分数据误差均小于20%，且相关系数均大于0.99，说明液体爬升高度控制方程和多孔介质比表面积理论模型的准确性和适用性较好。

图  5  预测结果和实验结果比较

Figure  5.  Comparison of the predicted and experimental results

下载: 全尺寸图片 幻灯片

5.   结论

（1）针对球形粉末烧结多孔介质，对n点相关函数g_n(rⁿ)进行修正，将无限项求和化简为有限项求和，推导出了针对球形粉末烧结多孔介质的比表面积理论模型。

（2）利用EDEM模拟不同粒径、不同孔隙率的多孔介质，拟合出不可变形系数（ε）的表达式。通过分析模拟数据，发现与传统模型相比，所建模型的准确性更高，适用范围更广。当孔隙率为0.10时，硬球模型的误差在100%以上；当孔隙率为0.45时，完全穿透模型的误差在50%以上，所建模型的误差始终低于6%。

（3）将伯努利方程与所建多孔介质比表面积理论模型耦合，在综合考虑多孔介质的内部结构、颗粒平均配位数、不可变形系数、边界效应以及动态接触角等因素影响的情况下，建立了多孔介质内液体爬升高度的控制方程。实际液体爬升测试显示，控制方程预测值与实验值的相关系数在0.99以上，说明比表面积模型具有较好的准确性和应用价值。