现在大多数分析仪器的测试都采用计算机控制、采集信息及用配套软件处理数据，这为分析和测试提供了很大方便，不过，这同时也带来一些风险，如若不考虑样品及测试的具体情况、贸然用软件默认的模板生成的测试报告其可信度就要特别注意。相对其他仪器，在使用静态体积法孔隙分析仪器的学生中发生此类风险会多一些，笔者认为其原因有二：一是不了解静态体积法气体吸附测试孔径分布是建立在人为假设上的；二是不了解孔径分布的相关参数设定对计算过程和对分析测试结果的影响。有吸附仪器应用专家直接指出：若简单地相信气体吸附仪预置的报告模板出的分析数据，对复杂微介孔材料的分析则可能得出错误或虚假的结论^[1]。

本文试着以中华人民共和国国家标准GB/T21650.2-2008《压汞法和气体吸附法测定固体材料孔径分布和孔隙度第2部分：气体吸附法分析介孔和大孔》^[2] (以下简称“标准”)思路介绍Barret、Joyner和Halenda方法^[3] (以下简写为BJH方法)计算介孔孔径分布的过程和步骤。帮助学生了解介孔吸附分析中参数的意义和数据处理的要点。

通常将孔径在2-50 nm范围的孔称为介孔、将孔径大于50 nm的孔称为大孔。液氮温度下氮气吸附最适合的孔径(孔宽)约在0.4-50 nm范围内，随着温度控制和压力测量技术的进步，目前该方法可用于孔径范围上限到约100 nm的多孔固体孔表征^[2]。

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

77 K时吸附在多孔固体表面上的氮气量是其压力的函数，测得的氮气吸附量可细分为膜厚变化量和毛细管凝聚或脱除量两部分。吸附过程中，随着气体压力的上升，在多孔物质的表面及孔壁发生单层、多层吸附并形成液膜，随后孔内发生毛细管凝聚并形成类似液体的弯月面。液膜厚度与气体压力、样品表面性质有关，并可用数学方程式描述^[4]；毛细管凝聚时的孔径与p/p₀ (p是氮气的吸附平衡压力、p₀是液氮温度下氮气的饱和蒸气压力，p/p₀是相对压力)的关系可用Kelvin方程描述。据此，在某一假设下，从实测样品得到的一组(相对压力、吸附量)吸附数据就可计算出其每两相邻数据之间的因膜厚变化引起的液膜体积变化量、因发生毛细管凝聚或脱除的相应孔径内凝聚体积的变化量，进而得到该假设下样品的孔径、孔体积和孔表面积分布数据。

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

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

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

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

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

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

上述3个步骤中的计算难点是确定上述“两方面贡献”的各自大小。

为保证计算准确和便于复核，笔者引用了“标准”附录A的原始数据，且按“标准”方法验算了所附推导出的数据，整理后的数据(本文分列成两个表，对部分行进行了省略)见表1、表2。

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

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

(1) $r_{\mathrm{k}}=-\frac{2 \sigma_{1} V_{\mathrm{ml}}}{R T_{\mathrm{b}} \ln \left(\frac{p}{p_{0}}\right)}=-\frac{0.953}{\ln \left(\frac{p}{p_{0}}\right)}=-\frac{0.414}{\lg \left(\frac{p}{p_{0}}\right)}$

式(1)中：r_k是凝聚在孔隙中吸附气体的曲率半径(在上式中的数值单位已经换算为nm)；σ_l是液态凝聚物液氮的表面张力(0.0088760 N·m^-1)；V_ml是液态凝聚物液氮的摩尔体积(0.034752 L·mol^-1)；R是气体常数(8.314 J·mol^-1·K^-1)；T_b是分析测试时的冷浴温度(77.35 K)。应该注意到式(1)中引用的参数都是77 K氮气的，因此BJH方法只适用液氮温度下氮气吸附的孔径分布计算，不能用于其他温度(如273 K)的氮气吸附，也不能用于77 K的其他气体(如氧气或氩气)吸附的孔径计算。

将得到的r_k值列入表1第2列，$ \overline {{r_{\rm{k}}}} $($ \overline {{r_{\rm{k}}}} $是两相邻曲率半径r_k的平均值，单位为nm)列入表1第3列。另外，还要注意计算得到的r_k是凝聚在孔隙中液氮的曲率半径，并不是孔的半径。

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

对于圆柱状孔，其半径可由已经算得的膜厚和曲率半径加和得到，即算式r_p= r_k + t，(式中r_p是圆柱状孔半径，单位为nm)将所得r_p值列入表1第6列，将$\overline {{r_{\rm{p}}}} $($\overline {{r_{\rm{p}}}} $是由两相邻压力计算得到的孔半径r_p的平均值，单位为nm)值列入表1的第7列。

d_p值(d_p为圆柱状孔的直径，单位为nm)由算式d_p = 2r_p得到，列入表1的第8列。$\overline {{r_{\rm{p}}}} $值($\overline {{r_{\rm{p}}}} $是由两相邻压力计算得到的孔直径d_p的平均值，单位为nm)，列入表1的第9列。

为了后续便于从孔中脱附凝聚物体积推算孔体积，引入Q值(定义Q为体积修正因子，在3.2.3节稍作讨论)。对于圆柱状孔Q值的几何意义是两个等高但半径不同的圆柱体体积之比，也即等于两半径平方比，无单位。计算见式(2)，其数值列入表1的第10列：

(2) $Q=\left[\frac{\overline{r_{\mathrm{p}}}}{(\overline{r_{\mathrm{k}}}+\Delta t)}\right]^{2}$

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

77 K下氮气在多孔固体表面随着气体压力升高发生吸附，然后形成液膜，继而在孔内凝聚填充，此时所吸附的氮气转变为“液氮”，并认为这种“液氮”其物理性质同宏观液氮物理性质相同。

表2第11列V_N值(V_N是吸附量、标准状态下的氮气体积值，单位cm³·g^-1 STP)是实测后换算成标准状态值。由第11列两相邻压力点计算得到ΔV值(ΔV是两相邻压力吸附气体量V_N之差，单位cm³·g^-1 STP)，列入表2第12列。将ΔV换算成ΔV_L (ΔV_L是将ΔV换算成液氮体积值，单位cm³·g^-1)，用算式ΔV_L= 0.0015468ΔV得到ΔV_L值，列入表2第13列，算式中的系数0.0015468是标准状态下氮气体积转换成77 K液氮体积的换算因子。

相对压力p/p₀ (表1第1列、表2第1列)和吸附量V_N (表2第11列)都是原始测量数据。表1的第2列至第10列是由p/p₀套用算式逐列推导得到的；表2的12列、13列是由V_N套用算式逐列推导得到的。这些推导得来的数据都是只来自一列原始数据。需要特别注意：从14列至19列各列则须由V_N和p/p₀推导出的数值共同参与计算才可，且须计算完一行数值后才可继续下一行计算。

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

如前所述，ΔV_L可细分为ΔV_t (ΔV_t是两相邻压力因膜厚变化引起的液态吸附质体积变化量，单位cm³·g^-1)和ΔV_k (ΔV_k是两相邻压力因毛细管凝聚或脱除引起的液态吸附质体积变化量，单位cm³·g^-1)两部分，用算式表示这两部分有ΔV_L = ΔV_t + ΔV_k关系，即要从换算成液氮体积吸附的每个差值Δ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列及以后列的数值。

计算思路：先计算出孔体积，再利用孔体积与圆柱状孔壁面积关系计算裸露出的孔壁面积。

计算ΔV_p (ΔV_p是孔体积值，单位cm³·g^-1)：ΔV_p是从已经计算得到的ΔV_k值(表2第1行第15列)根据算式ΔV_p = ΔV_k × Q得到，数值记入第1行第16列。

计算Δa_p (Δa_p是ΔV_p含有的孔壁面积值，单位m²·g^-1)：圆柱体体积和圆柱体圆周表面积与圆柱体半径之间有关系式Δa_p = 2 × (ΔV_p/r_p) × 10³ (式中10³是单位变换因子)，数值列入第1行第18列。

表2第17列是ΣΔV_p (ΣΔV_p是累积脱除的孔体积值，单位cm³·g^-1)，其某行的值是对第16列的之前各行数值的累加；第19列是ΣΔa_p (ΣΔa_p是累积暴露的表面积值，单位m²·g^-1)，其某行的值是对第18列的之前各行数值的累加。

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

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

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

后续各行的第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方法计算孔径分布的基本数据(即“标准”的附录A)。

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

笔者对上述的“孔径值的中点”理解为是该体积增量对应的上、下孔径的算术平均值，即$\overline {{d_{\rm{p}}}} $；对“体积增量除以确定该增量的上、下孔径之差的商ΔV_p/Δd_p”理解为若测试点非常密集、变化量趋近无限小时即可用dV_p/dd_p微分形式表示；对“体积增量除以确定该增量上的上、下孔径对数值之差的商ΔV_p/Δlogd_p”理解为若测试点非常密集、变化量趋近无限小时即可用dV_p/dlogd_p微分形式表示。因此，作图时相应的轴标分别用dV_p/dd_p代替ΔV_p/Δd_p，dV_p/dlogd_p代替ΔV_p/Δlogd_p。

为了准确描述孔径分布情况及绘制更直观的图表，特别是为了能表达最常见的4种孔径分布，从表1、表2中选取部分列生成表3，并在表3中增加和定义第20列(当前体积增量所对应孔径值的增量Δd_p)、21列(Δd_p的上、下孔径对数值之差Δlogd_p)、22列(即第9列$ \overline {{d_{\rm{p}}}} $的对数值log$ \overline {{d_{\rm{p}}}} $)、23列(体积增量除以确定该增量的上、下孔径之差的商ΔV_p/Δd_p)和第24列(体积增量除以确定该增量上的上、下孔径对数值之差的商ΔV_p/Δlogd_p)的计算方法和数值。

从第20列到第24列的各列计算过程数学定义清晰、运算简单，细节不再赘述。

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

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

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

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

有关c值的这些信息值得注意，或许是BJH方法在计算较小的介孔时误差比较大的原因之一。

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

(3) $\frac{\mathrm{d} V_{p}}{\mathrm{d} d_{p}}=\frac{1}{2.303 d_{p}} \times \frac{\mathrm{d} V_{p}}{\mathrm{d} \log d_{p}}$

人们把孔形状常分为典型的3类：平板型、圆筒型和球型。“标准”举的例子是将介孔全部看成圆柱状，在一些仪器厂家的分析软件中BJH方法没有孔形状选项，就此问题笔者请教过部分吸附仪器公司的应用专家，他们认为在介孔范围内其他形状孔与圆柱状孔的计算结果很相近，在计算时都选用了圆柱状孔模型。

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

BJH方法用于分析比较小的介孔(孔径小于10 nm)时会有低估孔径的现象，对更小的介孔(孔径小于4 nm)时会低估孔径达到20%^[7]；对于微孔(孔径小于2 nm)则不能使用BJH方法。这是由于微孔相邻孔壁间的相互作用很强，不能将所吸附的、液化的氮看成具有常规热力学性质的液体，不能再使用Kelvin方程^[2]描述孔内凝聚液与孔径关系。对于微孔内吸附质与孔径关系现在常用微孔填充、分子模拟方法或密度函数理论等定量描述^[8]。

介孔样品吸附测试通常得到Ⅳ型等温线，其典型特征是吸附分支和脱附分支不一致。人们对迟滞回线分类、回线形状与孔型的关系进行了很多研究和论述，详见相关文献^{[2, 5, 7, 9]}。

在用BJH方法处理有迟滞回线的数据时是选吸附分支还是选脱附分支？这是一个经常遇到的问题，面对选择应谨慎考虑。“标准”给出了如下选择建议：

(1)比较均匀的圆柱孔的相对简单的孔结构可能产生狭窄的H1型迟滞回线，此时往往采用脱附分支进行分析；

(2)如果观察到H2型回线，表明出现了连通、孔堵塞及相关的渗透现象，这意味着采用等温线的任何一部分都不完全稳妥，因为可能具有混合效应(即同时具有延迟凝聚和网络渗透)。如果采用一定的方法，考虑孔径对延迟凝聚现象，尤其是在孔隙流体的介稳态范围内的影响，则可以采用吸附支数据进行孔径分析；

(3)如果观察到所测等温线脱附支在p/p₀ = 0.42 (注意测试条件：77 K，N₂)附近时出现陡降(在3.5小节中进一步讨论)，则应用吸附支的数据计算孔径分布。

当要选用脱附支的数据计算孔径分布时要注意仔细观察吸附、脱附等温线图中的脱附支：在p/p₀= 0.42附近时是否有吸附量陡降的情况，若有陡降，在孔径分布图中将会出现不真实的峰(假峰)。

样品脱附时凝聚在孔中的吸附质随着p/p₀降低会发生解聚或蒸发，有复杂孔道的样品在脱附时可能出现吸附量陡降的情况。对于77 K吸附氮气条件下的样品，这种陡降通常在脱附到p/p₀ = 0.42附近时发生。这种陡降被认为是“张力强度效应”造成的，它与吸附温度和吸附质种类有关，同一样品在不同温度、用不同吸附质测试时发生陡降的相对压力值并不相同^{[2, 9]}。当用BJH方法计算有这种“陡降”脱附支的数据时，其所得数据若用在后续绘制孔径分布图时将会出现假峰。

前节采用了“标准”提供的数据进行计算，运算中所得的各个中间值和最终值没有出现负值或难于理解的情况，若今后用自己测试得到的数据推导时，有可能会在计算到孔径比较小时累积孔容积等出现负值的情况，计算到负值时终止^{[5, 10]}即可。出现负值^[4]有一些是人为原因。

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

用BJH模型分析数据时孔径下限受宏观模型限制不应低于2 nm，上限受测量条件限制目前比较可靠的是在100 nm左右，对低于2 nm或大于100 nm的孔径分布结果要特别慎重，这也是举例作孔径分布图时只选第4行至第39行(大致对应2 nm至100 nm孔径范围)共36对数据的原因。

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