超全面解读BJH方法如何计算孔径分布过程

1. 用BJH方法计算介孔孔径分布

通常将孔径在2nm-50nm范围的孔称为介孔、将孔径大于50nm的孔称为大孔。液氮温度下氮气吸附最适合的孔径（孔宽）约在0.4nm-50nm范围内，随着温度控制和压力测量技术的进步，目前该方法也可用于测定更大的孔隙度，可用于孔径范围在2.0nm至100nm之间多孔固体的孔表征。

气体吸附法分析介孔和大孔的方法经多年完善现在已经形成为测试的标准[1](简称“标准”)，其包含的在液氮温度下氮气吸附法以及对得到的数据采用Barret，Joyner和Halenda方法[2]（简写为BJH方法）处理结果被人们普遍接受并被编制成软件，经适当设置或选项后就得到分析或测试 告。

下面按照“标准”的思路 介绍BJH方法计算孔径分布的要点和步骤。

1.1 BJH方法要点

77K时吸附在多孔固体表面上的氮气量是其压力的函数，测得的氮气吸附量可细分为膜厚变化量和毛细管凝聚或脱除量两部分。

吸附过程中，随着气体压力的上升，在多孔物质的表面及孔壁发生多层吸附并形成液膜、孔内发生毛细管凝聚并形成类似液体的弯月面。液膜厚度与压力、样品性质有关，并可用不同的算式描述；毛细管凝聚时的孔径与p / p 0 的关系可用Kelvin方程描述。据此，在某一假设下，从实测样品得到的一组吸附数据就可计算出其每两相邻数据（压力、吸附量）之间的膜厚体积变化量、发生毛细管凝聚或脱除的孔径体积变化量，进而得到某一假设下该样品的孔径、孔体积和孔表面积分布数据。

1.2 BJH方法计算孔径分布的三个步骤

“标准”介绍了BJH计算介孔孔径分布有如下3个步骤：

（1） 不论采用等温线的吸附支还是脱附支，数据点均按压力降低的顺序排列；

（2） 将压力降低时氮气吸附体积的变化归于两方面的贡献：

a． 由Kelvin方程对高、低两个压力计算出的尺寸范围内的孔隙中毛细管凝聚物的脱除；

b． 脱除了毛细管凝聚物后的孔壁上多层吸附膜的减薄；

（3）为准确计算孔径和孔体积，必须考虑毛细管凝聚物从孔隙中脱除时会残留多层吸附膜。

以上3个步骤在计算过程中是比较繁琐的，计算难点是确定“两方面贡献”的各自大小。

2. 用BJH方法推算孔径分布的具体过程

为保证计算顺利， 中引用了“标准”所附的原始数据，且按“标准”方法验算了全部推导数据。还对“标准”数据格式作了调整、 对部分行进行了省略，整理、验算后的数据见表 1、表2。

2.1 由氮气相对压力值计算对应的凝聚曲率、膜厚和孔径

表1第1列相对压力数值是测量得到的，按计算步骤第1点要求由高到低排列。对于77K氮气吸附，假定样品中的孔皆为圆柱状孔，用Kelvin方程和膜厚描述算式就能计算出各个相对压力时的孔径和膜厚，结果列于表1第2列至第10列，各列计算过程详述如下。

在液氮温度下以氮气为吸附气体时随着相对压力增高气体会在样品孔内凝聚，以圆柱状孔为例，气体相对压力与其在孔内凝聚曲率关系可用Kelvin方程表述，见算式1：

式中：rk是凝聚在孔隙中吸附气体的曲率半径（在上式中的数值单位已经换算为nm）；σl是液态凝聚物的表面张力（0.0088760N/m ）；Vml是液态凝聚物的摩尔体积（0.034752 l/mol ）；R是气体常数（8.314J·mol-1·K -1 ）；Tb分析温度（77.35K）；p是氮气的吸附平衡压力（式中取值单位不限，但要与p0单位保持一致）；p0是液氮温度下氮气的饱和蒸气压（压力单位不限）。

方程引用的参数都是77K氮的,因此BJH方法只适用液氮温度下氮气吸附的孔径分布计算，不能用于其他温度下（如273K）的氮气吸附、也不能用于液氮温度下的其他气体（如氧气或氩气）吸附的孔径计算。将计算得到的曲率半径r K 值列入表1 第2 列，两相邻曲率半径r K 的平均值列入表1 第3 列。还要注意从上式计算得到的数值是凝聚在孔隙中的液态氮的曲率半径，并不是孔的半径。

多孔固体表面吸附气体后形成的液膜厚度与气体相对压力的关系是由样品性质决定的，不同样品其液膜厚度数学表达式可能不同，常见的有Haley方程等多种形式[3]，也可由试验或数学参比等温线获得。在仪器配套软件中也有多种选项。这里选用“标准”的：Haley方程t=0.354[-5/ln( p / p 0 )] 1/3 ( 用符号t表示液膜厚度，单位为nm)，式中0.354是单层液氮膜的厚度值，将计算得到的膜厚度值列入表1的第4列，将两相邻压力计算的膜厚变化量（写为Δ t ，单位为nm），列入表1的第5列。

对于圆柱状孔，其半径（写为rp表示，单位为nm）可由已经算得的膜厚和曲率半径加和得到，即由公式r p = r K + t 计算得到，将所得数值列入表1第6列，由两相邻压力计算得到的孔半径计算孔平均半径（写为，单位为nm），将所得数值列入表1的第7列。

圆柱状孔的直径（写为dp，单位为nm）由公式dp=2r p 得到，列入表1的第8列。由两相邻压力计算得到的孔直径计算平均直径（写为，单位为nm），列入表1的第9列。

为了后续便于从孔中脱附凝聚物体积推算孔体积，引入体积修正因子（写为Q）。Q值算式与假设的孔形状有关[4],[5]，且孔形状不同时差别很大，对于圆柱状孔Q值的几何意义是半径不同的圆柱体体积之比，也即等于两半径平方比，无单位。 见算式 2，其数值列入表1的第10列：

注意上述计算得到一共9列数值都只与相对压力有关而与吸附量无关。

2.2 由样品吸附气体值推算凝聚液体积

在液氮温度下氮气在多孔固体表面发生多层吸附并形成液膜、在孔内凝聚填充，所吸附的气体转变为液氮，并认为其物理性质同宏观液氮性质相同。

表2第11列吸附量（写为Δ V N ，标准状态下的氮气体积值，单位cm 3 ·g -1 STP ）是实测值。由第11列两相邻压力点计算得到吸附差值（写为Δ V ，单位cm 3 ·g -1 STP ），列入表2第12列。将Δ V 换算成液氮体积（写为Δ V L ，单位cm 3 ·g -1 ），用算式Δ V L = 0.0015468Δ V N 得到Δ V L 值 ，列入表2第13列，其中0.0015468是标准状态下氮气体积转换成77K液氮体积的因子。

实际上，上两小节中的计算都是简单的套用公式的数学推导或是数值单位变换过程。

2.3 由气体相对压力、气体吸附量值共同推算孔径、孔体积值和孔壁面积值

如前所述，测得的氮气吸附量可细分为膜厚变化量和毛细管凝聚或脱除量（写为Δ V k ，单位cm 3 ·g -1 ）两部分，用算式表示这两部分有Δ V L =Δ V t +Δ V k 关系，即要从换算成液氮体积吸附的每个差值中再求出每个膜厚变化值和孔隙凝聚液氮脱除值。笔者认为这是后续计算的重点和难点。

“标准”中有这样一段表述：Δ V t 为由未被毛细管凝聚物填充的孔壁面上的氮膜减薄而脱除的氮气体积。其第1行的数值为零，因为假定在最高压力时，所有介孔均被充满。在此特定实例中，由暴露孔壁的表面积来计算Δ V t 。

对上一表述笔者这样理解：确定第1个液膜厚度所占体积值和第1个孔解脱附体积值是计算的关键，假定在最高压力时，所有孔均被充满，设想此时样品没有纯粹的非孔表面，最高点和次高点脱除的液氮体积（Δ V L ）全部是由符合某些尺寸的孔脱除的，这部分孔的孔壁表面在相对压力次高点时才裸露出来，因此最高点与次高点之间的膜厚变化Δ V t =0 （记入表2 第1 行第14 列），相应的孔内凝聚液氮脱除体积数值Δ V k =Δ V L -Δ V t =Δ V L （记入表2 第1 行第15 列）。两数值在表2 中用加粗字标出。

后续行中两相邻压力点膜厚变化值并不为0 ，计算后续行第14、15列数值要 知道裸露孔壁面积值，故从第14列开始需要计算好上一行全部数值后才可开始计算当前行第14列及以后列的数值。

2.3.1 确定最初脱附的裸露孔壁表面积，完成数据第1行的计算

首先计算出解凝聚体积，再利用孔体积与圆柱状孔壁面积关系计算裸露出的孔壁面积。

计算解凝聚体积（写为Δ V p ，单位cm 3 g -1 ）：Δ V p 是从已经计算得到的Δ V k 值（表2 第1 行第15 列）根据算式Δ V p =Δ V k * Q 得到，数值记入第1 行第16列。

计算裸露的孔壁面积（写为Δ a p ，单位m 2 g -1 ）：圆柱体体积、柱体表面积与圆柱体半径之间有关系式Δ a p =2*/*10 3 （式中10 3 是单位变换因子），数值列入第1 行第18列。

表2 第17 列是 累积脱除的孔体积（写为ΣΔ V p ，单位m 2 g -1 ），其某行的值是对第16 列之前各行数值的累加；第19 列是 累积暴露的表面积（写为ΣΔ a p ，单位m 2 g -1 ），其某行的值是对第18 列之前各行数值的累加。

表2 第1 行中第16 列、18 列是第1 个得到的数值，是初值，其累加值就等于初值，故第1 行的第16 列和第17 列、第18 列和第19 列数值相同。这样就依次得到第1行的第15列至第19列的数值。

2.3.2 确定第2行膜厚变化体积和孔的凝聚体积，完成数据第2行及以后各行的计算

得到表2第1行第14列至第19列数值后，从第2行乃至最后一行的计算遵循标准规定的步骤（1.2小节中的（2）、（3）点）要求即可得到全部需要的数值。

注意：表2 第2 行第14 列数值（用加粗斜体字下划线标出）是本行后续列的计算起点，其值是根据算式Δ V t =0.85*10 -3 Δ t *ΣΔ a p 得到的。注意此处ΣΔ a p 数值（用斜体加粗字标出）是取自上一行19 列的值，不是当前行的。Δ V t 算式中第1个系数取0.85的原因将在3.2.2小节中试着说明和讨论；第2个系数10 -3 是算式中ΣΔ a p 单位取m 2 g -1 、Δt 单位取nm和计算结果Δ V t 单位取cm 3 g -1 时的单位变换因子。

后续各行的14 列值也都是从上一行第19 列值按算式Δ V t =0.85*10 -3 Δ t *ΣΔ a p 计算得到的。

第2行及以后各行第15至19列的计算过程：因有Δ V p =Q*Δ V k 关系，从第2行15列数据推算出Δ V p 值，记入第16列；将当前Δ V p 值与以前所有Δ V p 值加和得到ΣΔ V p ，记入第17列；因有Δ a p =2Δ V p / r p 关系，推算出新裸露出的孔壁面积，记入第18列；将当前Δ a p 值与以前所有的Δ a p 加和得到ΣΔ a p ，记入第19列。这样就依次算出了第2行和以后各行的第15列至第19列数值。

至此，完成了表2 的计算，得到了BJH 方法计算孔径分布的基础数据（即“标准”的附表）。

3. 将BJH方法推算结果列表和作图

为了得到明了的数据和绘制更直观的分布图，还需要对“标准”的附表数据进行再整合（见表3 ）。

“标准”关于作图有如下表述：计算出的孔径分布可以采用多种方式来表达，最常见的有4种：小于（或大于）孔径累积孔体积、增量孔体积对孔径、微分孔体积对孔径和对数微分孔体积对孔径。累积分布是指在特定的孔径范围内将大于或小于当前孔径的孔隙总体积作图或列表；孔体积增量分布是将计算出的两个连续孔径之间的绝对孔体积与用于计算当前增量的孔径值的中点作图或列表；微分是将体积增量除以确定该增量的上、下孔径之差的商与该增量的孔径值的中点作图或列表；对数微分分布是将体积增量除以确定该增量上的上、下孔径对数值之差的商与对增量的孔径值的中点 作图或列表。现在，孔的表面积也常用累积分布来表达并与孔体积作在同一张图上。

笔者对上述的“孔径值的中点”理解为是该体积增量对应的上、下孔径的平均值，即;对“体

积增量除以确定该增量的上、下孔径之差的商Δ V p /Δ d p ”理解为若测试点非常密集、变化量趋近无限小时即可用d V p /d d p 表示; 对“体积增量除以确定该增量上的上、下孔径对数值之差的商Δ V p /Δlog d p ”理解为若测试点非常密集、变化量趋近无限小时即可用d V p /dlog d p 表示。因此，作图时相应的轴标分别

用d V p /d d p 代替Δ V p /Δ d p 、d V p /dlog d p 代替Δ V p /Δlog d p 。

3.1 利用得到数据列表和作图

为了能清晰用列表方法表达最常见的4 种孔径分布，从表1 、表2 中选取个别列生成表3 ，并在表3 中增加并定义第20 列（当前体积增量所对应孔径值的增量Δ d p ）、21 列（该增量上的上、下孔径对数值之差Δlog d p ）、22 列 （即第9列的对数值log）、23 列（体积增量除以确定该增量的上、下孔径之差的商Δ V p /Δ d p ） 和第24列（体积增量除以确定该增量上的上、下孔径对数值之差的商Δ V p /Δlog d p ） 的计算方法和数值。从第20列到第24列的数学定义清晰、运算简单，细节不再赘述。

由列表转为作图时各列和图轴的关系见表4。根据坐标轴与数列对应关系，根据需要选取指定列作图。现将常见的两幅图附后：图1是等温吸附线图（选两列的全部数据）；图2是双纵轴的孔径-累积孔体积/孔表面积叠加图(选三列的第4行至第39行、对应孔径2～100nm范围的数据)。

双纵轴方法作图能更清晰和简明和地表示孔径分布，也更便于对比和判断，故经常会用到。孔体积、孔面积数据加和时有从大于孔径或小于孔径两个方向，图2选用的是大于孔径累积。

3.2 对数据列数值计算过程和作图的一些思考

在计算孔径分布过程和将数值再整合过程中感觉有些细节要给予特别注意，笔者试着探讨如下。

3.2.1 计算Δ V t 时算式中系数0.85的说明

Δ V t 算式中系数0.85来自 献[2]中的变量c， 献对c值的范围和变化规律做了较详细的讨论，简要概括如下：c 是一个与孔径有关的变量，并不是一个常数。考虑到孔径与p / p 0 的对应关系，得出了c 值对p / p 0 变化的趋势，当p / p 0 越来越小，直至趋近于0 时，c 则会越来越大，直至趋近于1 。还指出当c 值的变化在±0.05 范围内，对计算的孔隙体积分布的影响几乎可以忽略不计。为了简化计算，人为规定c 值是常数，并给出了最佳的c 值结果：c= 0.85 。其它 献对c 值也有介绍，有的仪器公司将其定义为弯曲液面校正值，也等于0.85 [3]。

有关c 值的这些内容要注意，或许这是这个方法在计算较小的介孔时误差会比较大的原因之一。

3.2.2 孔体积孔径微分和对数微分分布曲线峰值不等的原因

微分孔体积分布曲线和对数微分分布曲线的峰值出现的位置和数值通常是不相等的，这可从表3中的相关数据列看出：微分孔体积分布最大值出现在第23列第31行；对数微分孔体积分布最大值出现在第24列第12行（均用斜体加粗字下划线标出）。有的同学对这两种分布以为只是图中坐标横轴的变化，不应引起纵轴对应的峰值变化，但从图中看到的两个峰值明显不同，对此常有同学提出疑问，不等的原因是因为微分及对数微分有如算式3的数学关系，并不是简单的横轴标度变化。

3.2.3 圆柱状孔与其他形状孔的差别

人们把孔形状常分为典型的3 类：平板形、圆筒形和球形。“标准”将介孔的形状全部看成圆柱状，仪器厂商的分析软件中对BJH方法也都是用圆柱状孔，没有其他如狭缝（平板）、球形等孔模型选项，就此问题笔者请教过部分仪器公司技术人员，他们认为在介孔范围内其他形状孔与圆柱状孔的计算结果很相近。现在，不论是进口仪器还是国产仪器软件分析介孔时都没有孔型这一选项。

平板形、圆筒形和球形孔发生凝聚时的Kelvin方程[4]、体积修正因子Q[4],[5]的算式会因孔径不同其数学式差别很大，另也有人采用采用“无模型法”[6]计算的。今后若进一步研读 献、整理和分析 献所附数据，或许能看到一组吸附数据用不同孔模型处理时其孔径分布相近程度究竟有多大。

4. 小结

回顾前节的推导和计算过程、观察得到的数据、图表和几点探论，就会发现推导过程虽很繁琐但都还顺利、结果都还合理、没有出现负值或不可理解的地方，这完全是因为我们选用了“标准”所用的数据且按其步骤逐步推导的原因。今后若用自己测试得到的数据，很可能会在计算到孔径比较小时出现不合理情况如负值，出现负值时终止[4],[7]即可，再尝试查找原因（ 如孔径相距较大 [3]等） 。

我们还应注意到数据处理量的变化：原始测试数据是43组、计算后得到数据是42组、作等温吸附线图时选用了全部43组原始测试数据，作孔径分布图时则只选取了中间36组数据。

用BJH模型分析数据时孔径下限受宏观模型限制不应低于2nm、上限受测量条件限制目前比较可靠的是在100nm左右，对低于2nm或大于100nm的孔径分布要特别慎重，这也是举例作孔径分布图时只选第4行至第39行的原因。早期 献介绍Kelvin方程适用的孔径上限只有50nm左右[5]。

至此，我们完整了解了计算BJH孔径分布的推导过程，知道了引用了那些物理参数及条件、对孔隙做了哪些几何假设和忽略、为简化计算做了哪些近似以及受到的仪器测量条件的哪些限制，明白这些之后，对制作和研读分析测试 告、评价分析测试结果会有一定帮助，还可根据对已有测试数据的分析和判断，为后续调整测试参数、改进测试方法提供指导。

参 考 献

[1] 朱庆山, 黄 来, 周素红, 邹涛, 李凤霞, 王勇, 谢朝晖. 压汞法和气体吸附法测定固体材料孔径分布和孔隙度第2 部分：气体吸附法分析介孔和大孔：中华人民共和国国家标准,GB/T21650.2-2008. 2008-04-16

[2] Barrett, E. P.; Joyner, L. G.; Halenda, P. P. J. Am. Chem. Soc 1951,73,373-380.

[3] 金彦任, 黄振兴. 吸附与孔径分布. 北京: 国防工业出版社, 2015: 57?58,128.

[4] 严继民张启元高敬琮吸附与凝聚- 固体的表面与孔（第二版）. 北京: 科学出版社, 1986: 168?171,200.

[5] S.J. 格雷格, K.S.W.辛. 吸附、比表面与孔隙率. 高敬琮, 刘西尧,. 北京: 化学工业出版社, 1989: 121, 141-151.

[6] 严继民, 张启元. 分子科学与化学研究，1983，(2), 75.

[7] 严继民, 张启元. 化学学 ，1977，35 (1,2), 13.

点击阅读原 ，提交计算需求！