用于确定储氢合金吸/放氢热力学参数的方法

摘要

本发明涉及一种合金材料关键性能参数的获取方法，特别是一种用于确定储氢合金吸/放氢热力学参数的方法，首先通过储氢合金吸放氢性能测试仪进行相关测试，获得不同温度下合金吸/放氢压力-浓度等温线数据，然后对这些实验数据进行一系列自动分析处理，最终输出储氢合金吸/放氢焓变、吸/放氢熵变和在设定温度下的吉布斯自由能变数据。本发明可以大大改善工程应用中对储氢合金成分-性能筛选的效率和分析结论的可靠性，解决了以往储氢合金吸/放氢热力学参数数据分散性大，出现错误的概率大，分析测试过程繁琐，数据输出不够直观等问题。

说明书

技术领域

本发明涉及一种合金材料关键性能参数的获取方法，特别是一种用于确定储氢合金吸/放氢热力学参数的方法。 

背景技术

上世纪六十年代人们发现一些合金能够快速大量地可逆吸放氢气，从而引发了世界范围内对储氢合金的研究热潮。目前储氢合金已在氢气的储运和压缩、同位素分离和储存、热泵与空调、镍氢电池、金刚石合成、氨合成、永磁体生产等众多领域得到广泛应用。在各种应用中，合金吸/放氢过程所伴随的焓变ΔH^o、熵变ΔS^o、吉布斯自由能变ΔG^o是必须弄清的重要参数，因为它们决定了不同应用领域对储氢合金的选择，以及某个温度下，反应的方向。比如，对于热泵的应用，一般要求储氢合金吸氢焓变的范围为-29.2~-37.6 kJ/molH₂；对于电池的应用，一般要求吸氢焓变的范围在-25~-50 kJ/molH₂；对于燃料电池汽车储氢的应用，一般要求放氢焓变范围在15~24 kJ/molH₂。由于氢化反应的熵变主要取决于分子氢向溶解氢的转变，因此不同合金的氢化反应其熵变基本相同，在-130 J/(molH₂K)左右。尽管如此，对于氢增压应用，获得比较精确的熵变值仍是非常必要的。 

对于合金吸/放氢过程焓变和熵变的传统获取方法是：首先通过Sieverts装置进行步进式充放氢实验，测定合金吸/放氢静态压力-浓度等温线，然后由实验者根据测得的等温线估计出不同温度下的反应平台压，依据范特霍夫关系拟合获得焓变和熵变。由于很多环节都是人为操作完成，特别是其中的反应平台压是人为估计得到的，因此采用传统方法无法实现计算机自动分析，一方面造成分析工作量大，另一方面导致相关过程随意性强，没有比较固定的标准，所获得的数据分散性较大，出现错误的概率大大增加，不利于得出可靠的研究结论。对于吉布斯自由能变，实验者一般也是根据已经得到的焓变和熵变手工计算得到，因此上述问题也同样存在。 

发明内容

本发明的目的是提供了一种用于确定储氢合金吸/放氢热力学参数的方法，快速、高效、准确，克服了以往储氢合金吸/放氢热力学参数数据误差大，可靠性低，分析过程工作量大，无法实现计算机自动分析的缺点。 

本发明的目的是通过以下技术方案实现的，一种用于确定储氢合金吸/放氢热力学参数的方法，包括以下步骤： 

（1）称取储氢合金样品放入储氢合金吸/放氢性能测试仪的样品室中； 

（2）由低温到高温，在不同温度下，对储氢合金样品进行储氢合金吸/放氢压力-浓度等温线的测试，获得储氢合金吸/放氢压力-浓度等温线测试数据；

（3）在MATLAB软件的Workspace中建立一个包含有多列测试数据的数组，将获得的储氢合金吸/放氢压力-浓度等温线测试数据中的样品室温度、氢浓度和氢压力数据按先后顺序分别填入对应的数组中得到实验数据数组；

（4）将实验数据数组进行分割，得到对应不同温度且相互关联的储氢合金吸/放氢压力-浓度等温线数据数组；

（5）将每个储氢合金吸/放氢压力-浓度等温线数据数组中含零数据行去掉；

（6）将每个储氢合金吸/放氢压力-浓度等温线数据数组继续分割成吸氢和放氢两部分，分割点为每条储氢合金吸/放氢压力-浓度等温线的反应平衡压最大值位置；

（7）分别对每个储氢合金吸放氢压力-浓度等温线数据数组中吸氢和放氢数据进行保形插值处理，获得大量内插数据点；

（8）步骤（7）得到的内插数据点的横坐标表示氢浓度n，纵坐标表示反应平衡压p，对内插数据点进行 微分处理获得的极大值，其极大值所对应的反应平衡压作为吸/放氢压力-浓度等温线的反应平台压，相继获得不同温度下储氢合金吸放氢压力-浓度等温线反应平台压，将不同温度下储氢合金吸放氢压力-浓度等温线反应平台压中属于吸氢过程的构成一个数组，属于放氢过程的构成一个数组；

（9）将吸/放氢过程所有测试点对应的样品室温度取平均值，该平均值作为储氢合金吸/放氢压力-浓度等温线的反应温度，将储氢合金吸/放氢压力-浓度等温线的反应温度中属于吸氢过程的构成一个数组，属于放氢过程的构成一个数组，并保证各温度数组中数据顺序与上述相关反应平台压数组中的数据顺序一致； 

（10）依据范特霍夫公式，其中为反应平台压，为标准大气压，R表示气体常数，T表示反应温度，回归出吸/放氢反应的焓变ΔH^o、熵变ΔS^o，然后根据吉布斯自由能变计算公式ΔG^o=ΔH^o-T·ΔS^o得到吉布斯自由能变ΔG^o；

（11）依据Workspace建立的包含多列测试数据，以氢浓度为横坐标，以反应平衡压为纵坐标作图，其纵坐标采用自然对数坐标，得到不同温度的储氢合金吸/放氢压力-浓度等温线图，同时以1000/T为纵坐标，以对应的反应平台压的自然对数值为纵坐标作图，保证与前述等温线图的纵坐标压力值上下限相一致，同时对数据点进行范特霍夫公式拟合；

（12）输出储氢合金吸/放氢压力-浓度等温线图，范特霍夫曲线图，吸/放氢反应的焓变ΔH^o、熵变ΔS^o以及设定温度下吉布斯自由能变ΔG^o。

优选地，进行储氢合金吸/放氢压力-浓度等温线测试时，所使用样品量为1-5g。 

优选地，进行储氢合金吸放氢压力-浓度等温线测试时，样品室恒温温度选取3-6个不同温度，温度间距为10-30℃。 

优选地，MATLAB软件的Workspace中建立的测试数据数组包含有9-18列，并且按一列温度、一列氢浓度、一列氢压力的顺序排列。 

优选地，对吸氢和放氢数据分别进行保形插值处理，获得的内插数据点各为8000-10000个。 

优选地，上述步骤（4）至步骤（12）的数据处理过程均由MATLAB程序自动完成。 

优选地，上述步骤（4）至步骤（12）的数据处理过程根据数组运算法则进行，并且不存在任何循环结构。 

与现有技术相比，本发明具有以下有益效果：（1）可实现整个数据处理过程的计算机自动分析；（2）同时有数据图输出，保证分析过程的直观性；（3）这些数据图输出格式满足常规出版要求，可以直接用于研究报告和杂志论文中；（4）大大提高储氢合金吸/放氢热力学参数分析过程的效率；（5）大大提高储氢合金吸/放氢热力学参数的准确性和可靠性。 

附图说明

图1是本发明一种用于确定储氢合金吸/放氢热力学参数的方法的流程示意图； 

图2是本发明LaNi₅合金的吸/放氢压力-浓度等温线图；

图3是本发明LaNi₅合金的范特霍夫曲线图。

具体实施方式

下面结合附图对本发明做进一步详述： 

一种用于确定储氢合金吸/放氢热力学参数的方法，包括以下步骤：

（1）称取储氢合金样品放入储氢合金吸/放氢性能测试仪的样品室中； 

（2）由低温到高温，在不同温度下，对储氢合金样品进行储氢合金吸/放氢压力-浓度等温线的测试，获得储氢合金吸/放氢压力-浓度等温线测试数据；

（3）在MATLAB软件的Workspace中建立一个包含有多列测试数据的数组，将获得的储氢合金吸/放氢压力-浓度等温线测试数据中的样品室温度、氢浓度和氢压力数据按先后顺序分别填入对应的数组中得到实验数据数组；

（4）将实验数据数组进行分割，得到对应不同温度且相互关联的储氢合金吸/放氢压力-浓度等温线数据数组；

（5）将每个储氢合金吸/放氢压力-浓度等温线数据数组中含零数据行去掉；

（6）将每个储氢合金吸/放氢压力-浓度等温线数据数组继续分割成吸氢和放氢两部分，分割点为每条储氢合金吸/放氢压力-浓度等温线的反应平衡压最大值位置；

（7）分别对每个储氢合金吸放氢压力-浓度等温线数据数组中吸氢和放氢数据进行保形插值处理，获得大量内插数据点；

（8）步骤（7）得到的内插数据点的横坐标表示氢浓度n，纵坐标表示反应平衡压p，对内插数据点进行微分处理获得的极大值，其极大值所对应的反应平衡压作为吸/放氢压力-浓度等温线的反应平台压，相继获得不同温度下储氢合金吸放氢压力-浓度等温线反应平台压，将不同温度下储氢合金吸放氢压力-浓度等温线反应平台压中属于吸氢过程的构成一个数组，属于放氢过程的构成一个数组；

（9）将吸/放氢过程所有测试点对应的样品室温度取平均值，该平均值作为储氢合金吸/放氢压力-浓度等温线的反应温度，将储氢合金吸/放氢压力-浓度等温线的反应温度中属于吸氢过程的构成一个数组，属于放氢过程的构成一个数组，并保证各温度数组中数据顺序与上述相关反应平台压数组中的数据顺序一致； 

（10）依据范特霍夫公式，其中为反应平台压，为标准大气压，R表示气体常数，T表示反应温度，回归出吸/放氢反应的焓变ΔH^o、熵变ΔS^o，然后根据吉布斯自由能变计算公式ΔG^o=ΔH^o-T·ΔS^o得到吉布斯自由能变ΔG^o；

（11）依据Workspace建立的包含多列测试数据，以氢浓度为横坐标，以反应平衡压为纵坐标作图，其纵坐标采用自然对数坐标，得到不同温度的储氢合金吸/放氢压力-浓度等温线图，同时以1000/T为纵坐标，以对应的反应平台压的自然对数值为纵坐标作图，保证与前述等温线图的纵坐标压力值上下限相一致，同时对数据点进行范特霍夫公式拟合；

（12）输出储氢合金吸/放氢压力-浓度等温线图，范特霍夫曲线图，吸/放氢反应的焓变ΔH^o、熵变ΔS^o以及设定温度下吉布斯自由能变ΔG^o。

进行储氢合金吸/放氢压力-浓度等温线测试时，所使用样品量为1-5g。 

进行储氢合金吸放氢压力-浓度等温线测试时，样品室恒温温度选取3-6个不同温度，温度间距为10-30℃。 

MATLAB软件的Workspace中建立的测试数据数组包含有9-18列，并且按一列温度、一列氢浓度、一列氢压力的顺序排列。 

对吸氢和放氢数据分别进行保形插值处理，获得的内插数据点各为8000-10000个。 

上述步骤（4）至步骤（12）的数据处理过程均由MATLAB程序自动完成。 

上述步骤（4）至步骤（12）的数据处理过程根据数组运算法则进行，并且不存在任何循环结构。 

实施例

制备储氢合金LaNi₅，其原料金属为La 99 wt%，Ni 99.9 wt%，采用氩气保护下的真空感应熔炼，熔炼后将合金装入真空度为1 Pa的石英管中进行退火，并在1000 ℃下保温6 h后迅速淬入水中。退火样品经化学分析，以验证成分的准确性。最后，称取大约1g的样品放入储氢合金吸放氢性能测试仪的样品室中，在恒温水浴的环境中，进行40℃、60℃和80℃的吸/放氢压力-浓度等温线的测试。测试获得如下表1，表2和表3的数据中。然后打开MATLAB软件在Workspace窗口中建立名为PCT的数组，将这些等温线实测数据拷入PCT数组中，同时确定吉布斯自由能变的计算温度为T=298 K。在MATLAB的命令窗（Command window）中执行命令HSG=cal_HSG(PCT,T)后计算得到吸氢反应的焓变为ΔH^o=-29.6071 kJ/molH₂，熵变为ΔS^o=-0.1076 kJ/(molH2K)，吉布斯自由能变为ΔG^o_298 K=2.4511 kJ/molH₂；放氢反应的焓变为ΔH^o=30.1288 kJ/molH₂，熵变为ΔS^o=0.1076 kJ/(molH₂K)，吉布斯自由能变为ΔG^o_298 K= -1.9413 kJ/molH₂，输出吸/放氢压力-浓度等温线图（图2）和范特霍夫曲线图（图3）。 

表1 40 ℃下LaNi5合金吸/放氢压力-浓度等温线测试数据 

。

表2 60 ℃下LaNi5合金吸/放氢压力-浓度等温线测试数据 

。

表3 80 ℃下LaNi5合金吸/放氢压力-浓度等温线测试数据 

。

 其中命令HSG=cal_HSG(PCT,T)完成如下过程： 

（1）  将PCT数组进行分割，得到对应不同恒温温度的相互关联的3个储氢合金吸/放氢压力-浓度等温线数据数组；

（2）  将每个储氢合金吸/放氢压力-浓度等温线数据数组中含零数据行去掉；

（3）  将每个储氢合金吸/放氢压力-浓度等温线数据数组继续分割成吸氢和放氢两部分，分割点为每条储氢合金吸/放氢压力-浓度等温线的反应平衡压最大值位置；

（4）  分别对每个储氢合金吸放氢压力-浓度等温线数据数组中吸氢和放氢数据进行保形插值处理，分别获得10000个内插数据点；

（5）  对所有这些内插数据点进行微分处理，并寻找其极大值所对应的平衡压，该反应平衡压作为吸/放氢压力-浓度等温线的反应平台压，通过这种方法相继获得不同温度的储氢合金吸/放氢压力-浓度等温线反应平台压，将这些反应平台压属于吸氢过程的构成一个数组，属于放氢过程的构成一个数组；

（6）  将每个吸/放氢过程所有测试点对应的样品室温度取平均值，该平均值作为储氢合金吸/放氢压力-浓度等温线的反应温度，将这些反应温度属于吸氢过程的构成一个数组，属于放氢过程的构成一个数组，并保证各温度数组中数据顺序与上述相关反应平台压数组中的数据顺序一致；

（7）  然后依据范特霍夫公式和数组运算法则，回归出吸/放氢反应的焓变ΔH^o、熵变ΔS^o，然后根据吉布斯自由能变计算公式ΔG^o=ΔH^o-T·ΔS^o得到设定温度下的吉布斯自由能变ΔG^o；

（8）  以氢浓度为横坐标，以反应平衡压为纵坐标作图，且纵坐标采用自然对数坐标，得到不同温度的储氢合金吸/放氢压力-浓度等温线图，同时以1000/T为纵坐标，以对应的反应平台压的自然对数值为纵坐标作图，保证与吸/放氢压力-浓度等温线图的纵坐标压力值上下限相一致，同时对数据点进行范特霍夫公式拟合；

（9）  最后输出储氢合金吸/放氢压力-浓度等温线图，范特霍夫曲线图，吸/放氢反应的焓变ΔH^o、熵变ΔS^o以及设定温度下吉布斯自由能变ΔG^o。

以下是上述数据处理过程的MATLAB程序： 

function HSG=cal_HSG(PCT,T)

PCT1=PCT(:,1:3);

        [i1,~]=find(PCT1==0);

        if isempty(i1)

            PCT1=PCT1;

        else

            PCT1=PCT1(1:i1-1,:);

        end

        t1=PCT1(:,1);

        c1=PCT1(:,2);

        p1=PCT1(:,3);

        [~,j1]=max(p1);

        abt1=mean(t1(1:j1));

        det1=mean(t1(j1:end));

        abHM1=c1(1:j1);

        abP1=p1(1:j1);

        abx1=linspace(abHM1(1),abHM1(end),10000);

        abfitResults1=pchip(abHM1,abP1);

        aby1 = ppval(abfitResults1, abx1);

        [~,abk1]=max(diff(abx1)./diff(log(aby1)));

        abPe1=aby1(abk1);

        deHM1=c1(j1:end);

        deP1=p1(j1:end);

        dex1=linspace(deHM1(end),deHM1(1),10000);

        defitResults1=pchip(deHM1,deP1);

        dey1 = ppval(defitResults1, dex1);

        [~,dek1]=max(diff(dex1)./diff(log(dey1)));

        dePe1=dey1(dek1);

        %PCT2

        PCT2=PCT(:,4:6);

        [i2,~]=find(PCT2==0);

        if isempty(i2)

            PCT2=PCT2;

        else

            PCT2=PCT2(1:i2-1,:);

        end

        t2=PCT2(:,1);

        c2=PCT2(:,2);

        p2=PCT2(:,3);

        [~,j2]=max(p2);

        abt2=mean(t2(1:j2));

        det2=mean(t2(j2:end));

        abHM2=c2(1:j2);

        abP2=p2(1:j2);

        abx2=linspace(abHM2(1),abHM2(end),10000);

        abfitResults2=pchip(abHM2,abP2);

        aby2 = ppval(abfitResults2, abx2);

        [~,abk2]=max(diff(abx2)./diff(log(aby2)));

        abPe2=aby2(abk2);

        deHM2=c2(j2:end);

        deP2=p2(j2:end);

        dex2=linspace(deHM2(end),deHM2(1),10000);

        defitResults2=pchip(deHM2,deP2);

        dey2 = ppval(defitResults2, dex2);

        [~,dek2]=max(diff(dex2)./diff(log(dey2)));

        dePe2=dey2(dek2);

        %PCT3

        PCT3=PCT(:,7:9);

        [i3,~]=find(PCT3==0);

        if isempty(i3)

            PCT3=PCT3;

        else

            PCT3=PCT3(1:i3-1,:);

        end

        t3=PCT3(:,1);

        c3=PCT3(:,2);

        p3=PCT3(:,3);

        [~,j3]=max(p3);

        abt3=mean(t3(1:j3));

        det3=mean(t3(j3:end));

        abHM3=c3(1:j3);

        abP3=p3(1:j3);

        abx3=linspace(abHM3(1),abHM3(end),10000);

        abfitResults3=pchip(abHM3,abP3);

        aby3 = ppval(abfitResults3, abx3);

        [~,abk3]=max(diff(abx3)./diff(log(aby3)));

        abPe3=aby3(abk3);

        deHM3=c3(j3:end);

        deP3=p3(j3:end);

        dex3=linspace(deHM3(end),deHM3(1),10000);

        defitResults3=pchip(deHM3,deP3);

        dey3 = ppval(defitResults3, dex3);

        [~,dek3]=max(diff(dex3)./diff(log(dey3)));

        dePe3=dey3(dek3);

        abt=[abt1;abt2;abt3];

        det=[det1;det2;det3];

        abPe=[abPe1;abPe2;abPe3];

        dePe=[dePe1;dePe2;dePe3];

        subplot(1,2,1);

        hold on

        plot(c1,p1,'-rs','DisplayName','40 ℃');

        plot(c2,p2,'-bo','DisplayName','60 ℃');

        plot(c3,p3,'-k>','DisplayName','80 ℃');

        set(gca,'YScale','log','FontSize',18,'FontName','times new roman');

        % Create xlabel

        xlabel('Hydrogen Content (wt.%)','FontSize',20,'FontName','Times New Roman');

        % Create ylabel

        ylabel('Pressure (MPa)','FontSize',20,'FontName','Times New Roman');

        set(gca,'YLim',[0.1 10])

        box(gca,'on');

        % Create legend

        legend2 = legend(gca,'show');

        set(legend2,'YColor',[1 1 1],'XColor',[1 1 1],...

        'Position',[0.283040364583333 0.719823340900192 0.109375 0.144770736733642]);

subplot(1,2,2);

abplot=plot(1000./(abt+273.15),log(abPe./0.101325),...

    'MarkerFaceColor',[1 0 0],...

    'MarkerEdgeColor',[1 0 0],...

    'Marker','square',...

    'LineStyle','none',...

    'DisplayName','absorption');

hold on

set(gca,'FontSize',18,'FontName','times new roman');

box(gca,'on');

% Get xdata from plot

abxdata = get(abplot, 'xdata');

% Get ydata from plot

abydata = get(abplot, 'ydata');

% Make sure data are column vectors

abxdata = abxdata(:);

abydata = abydata(:);

% Find x values for plotting the fit based on xlim

abxplot = linspace(min(abxdata)-0.1,max(abxdata)+0.1);

% Find coefficients for polynomial (order = 1)

abfit = polyfit(abxdata, abydata, 1);

% Evaluate polynomial

abyplot = polyval(abfit, abxplot);

% Plot the fit

abLine = plot(abxplot,abyplot,'DisplayName','fitting line','Color',[0 0 0]);

HSG(1,1)=-abs(abfit(1)*8.314472);

HSG(1,2)=-abs(abfit(2)*8.314472)/1000;

deplot=plot(1000./(det+273.15),log(dePe/0.101325),'MarkerFaceColor',[1 1 1],...

    'MarkerEdgeColor',[1 0 0],...

    'Marker','square',...

    'LineStyle','none',...

    'DisplayName','desorption');

% Get xdata from plot

dexdata = get(deplot, 'xdata');

% Get ydata from plot

deydata = get(deplot, 'ydata');

% Make sure data are column vectors

dexdata = dexdata(:);

deydata = deydata(:);

% Find x values for plotting the fit based on xlim

dexplot = linspace(min(dexdata)-0.1,max(dexdata)+0.1);

% Find coefficients for polynomial (order = 1)

defit = polyfit(dexdata, deydata, 1);

% Evaluate polynomial

deyplot = polyval(defit, dexplot);

% Plot the fit

deLine = plot(dexplot,deyplot,'LineStyle','--','DisplayName','fitting line',...

    'Color',[0 0 0]);

set(gca,'YLim',log([0.1 10]./0.101325),'XTickLabel',{'2.7','3','3.3'},...

    'XTick',[2.7 3 3.3]);

% Create xlabel

xlabel('1000 / \itT\rm (K^-^1)','FontSize',20,'FontName','Times New Roman');

% Create ylabel

ylabel('ln(\itp / \rmp_0)','FontSize',20,'FontName','Times New Roman');

% Create legend

legend1 = legend(gca,'show');

set(legend1,'YColor',[1 1 1],'XColor',[1 1 1],...

'Position',[0.709309895833333 0.653714054642302 0.1552734375 0.190623390005152]);

HSG(2,1)=abs(defit(1)*8.314472);

HSG(2,2)=abs(defit(2)*8.314472)/1000; 

HSG(:,3)=HSG(:,1)-T*HSG(:,2);