基于贝叶斯正则化BP神经网络的 超临界二氧化碳流动阻力特性预测

伍 权，朱 萌，陈 磊，卿浩然，周琳刚，许 凯， 苏 胜，胡 松，刘 辉，向 军

（1.国家能源集团新能源技术研究院有限公司，北京 102209； 2.华中科技大学煤燃烧国家重点实验室，湖北 武汉 430074）

超临界二氧化碳（S-CO2）具有稳定、易得、无毒、压缩性好和流动性强等特点，相比于传统水循环，基于S-CO2工质的新型动力循环具有体积小、效率高、灵活性好[1]等优势。目前，基于S-CO2循环的研究在光热发电[2-3]、火力发电[4-5]、核能[6-7]以及余热利用[8-9]等方面发展迅速。S-CO2新型动力循环的工质阻力特性会对S-CO2循环效率以及各部件设计产生较大影响，若无法准确计算循环管路和主要部件的阻力，S-CO2系统将无法准确设计并高效可靠运行。因此，S-CO2的流动阻力特性可靠评估是S-CO2技术商业化推进中至关重要的一环。

由于湍流的脉动性和新型动力系统实验数据缺乏，现有研究尚未形成对S-CO2摩擦阻力特性的统一理论。Kirillov等人[10]认为管截面密度分布是影响非等温流体阻力特性的关键因素，提出了壁面密度（依靠内壁温计算的流体密度）和主流密度比修正项来对等温流体公式进行修正。Fang等人[11]利用现有文献中的S-CO2以及其他流体数据，采用大量函数类型进行回归测试，提出了壁面黏度与主流黏度比修正项以及膜密度与拟临界点密度比的指数修正项。张海松等[12]研究了S-CO2在拟临界点附近的阻力系数峰值现象，认为传热恶化会对流动阻力特性造成负面影响，并采用超临界沸腾数SBO和普朗特数比作为修正项计算S-CO2摩擦阻力。然而，上述各修正项可能在一定的参数范围内适用，而在其他参数范围将不再是决定性因素，这给S-CO2流动阻力计算的工程化应用造成困难，设计人员难以根据需要采取合理修正项。

基于此，人工神经网络算法被认为是实现S-CO2阻力预测的一种有效途径。该方法不考虑流动传热的基本理论，而是通过大量训练，学习数据之间的特征和规律，并基于该规律预测其他数据的特征。这种方式的泛化能力强，不再局限于某一特定公式，可以预测更广参数范围的阻力特征。同时，采用神经网络可以直接将测量值例如温度、压力等作为输入，阻力作为输出，不再求解不同状态下的特征参数例如雷诺数、普朗特数等，有利于工程化实现。Ma等人[13]收集了14篇文献中的超临界水传热数据，利用BP（误差反向传播）神经网络进行了数据训练以预测其对流换热系数，结果表明其均方差为4.13%，相关系数R为0.971 7。Balcilar等人[14]基于垂直光滑铜管中的R134a冷却数据，训练了多个神经网络模型以预测其传热和压降特性，预测误差可保证在5%以内。Shadloo等人[15]采用多层感知器神经网络对长管道中的气体/非牛顿流体两相流的阻力特性进行了研究，预测值与实验数据的偏差在大多数情况下不超过5%。现有研究中针对S-CO2阻力特性进行神经网络预测的相关文献相对较少，神经网络对于其流动阻力的预测能力有待验证。此外，现有非绝热S-CO2实验数据[12,16-18]主要集中在CO2临界点（7.38 MPa、31 ℃）[19-20]附近，且质量流速较低（<1 500 kg/(m2·s)）、流体温度较低（<200 ℃），无法满足光热、火电以及核能等领域的实际需求。

聚焦于解决S-CO2发电领域中S-CO2流动摩擦阻力特性的预测问题，本文对更宽参数范围（750~ 2 200 kg/(m2·s)、10~20 MPa、200~340 kW/m2、60~500 ℃）的S-CO2阻力特性进行实验测量，分析压力、质量流速以及热流密度对摩擦阻力的影响，并根据实验数据，建立贝叶斯正则化BP神经网络模型，以实现对S-CO2摩擦阻力特性的准确预测。

1 实验系统

1.1 S-CO2循环系统

S-CO2流动阻力特性研究实验系统如图1所示。该系统首先通过加压柱塞泵将储罐中的低压CO2（纯度99.99%）增压至S-CO2循环回路。S-CO2循环回路由循环泵克服流动阻力，S-CO2依次流过质量流量计、回热器冷端、预热器、实验段、回热器热端以及冷却器，最终再次进入循环泵形成完整循环。冷却器采用与常规火电站相同的水冷循环实现对S-CO2的冷却，以保证进入循环泵的S-CO2在50 ℃以内。

实验时，通过加压柱塞泵和泄压阀配合实现对实验循环压力的调节；通过控制循环泵频率实现对实验循环流量的调节；通过控制预热器和实验段的加热功率实现对实验段进口温度以及实验段热流密度的调节。系统设计参数为：最大工质压力30 MPa，最大加热功率150 kW，最大流量1 800 L/h。

1.2 实验段

实验段设计如图1所示。实验段管材为Inconel 625合金钢，管内径为12 mm，外径为20 mm，在2个电极板之间的加热段长度为1 870 mm。

除加热段外，两端各留有1 000 mm的过渡段，用于稳定加热段的流动。实验段采用电阻发热原理，将低电压大电流电源连接到实验段两端，即2个电极板处，钢管本身的电阻使其产生均匀的加热热流。同时，实验段管外壁包裹了保温层以降低散热损失。加热段沿管长方向均匀布置有28个热电偶，用于测量管外壁温度。实验段的入口和出口分别安装有测量工质温度的热电偶和测量工质压力及压差的取压孔。实验段的流向为垂直向上，与重力方向相反。

1.3 数据采集

实验中，工质温度、壁温测量分别采用Omega 工程公司生产的K型铠装热电偶、K型热电偶丝，测温精度为0.75%。工质压力和压差测量采用北京俄华通仪表技术有限公司生产的压力、压差变送器，测压精度为0.25%F.S。工质流量测量采用北京首科实华自动化设备有限公司生产的质量流量计，测量精度0.2%。通过模拟信号将上述数据传输至由研华科技（中国）有限公司生产的数据采集卡，并转换为数字信号传入电脑。该数据采集卡对模拟信号的识别精度为0.1%，数据采集间隔为2 s。当实验循环各关键参数稳定后，开始进行数据采集，选取稳定状态下的60 s数据取平均后作为该次测量的最终数据，以避免波动误差。

2 数据处理与神经网络模型

2.1 数据处理

2.1.1 阻力计算

由于实验管段竖直，因此实验系统中测量的压差包括摩擦阻力ΔPf、重力压降ΔPg和加速压降ΔPa。摩擦阻力ΔPf的测量值可由式(1)计算。

式中：ΔP为总压降，kPa。

摩擦阻力ΔPf还可通过理论计算得到：

式中：f为摩擦阻力系数；L为管长，m；D为管内径，m；G为流体的质量流速，kg/(m2·s)；ρ为流体密度，kg/m3。

对于重力压降ΔPg，由于S-CO2的密度变化在实验参数范围内是非线性的，因此直接通过简单平均处理存在较大误差，本文采用目前认可度较高的Ornatskiy公式[21]来计算：

式中：h为比焓，kJ/kg；g为重力加速度，m/s2。

对于加速压降ΔPa，采用式(4)计算：

流体比焓、密度等物性参数通过NIST REFPROP 9.1软件查询获取。实验段沿长度方向不同位置z的主流比焓为：

式中：hb,z为位置z的主流比焓，kJ/kg；z为沿长度方向距实验段入口的距离，m。

实验段的加热热流密度指内壁面传递给流体的热流密度：

式中：qwall为内壁面热流密度，kW/m2。

由于实验只能测量管外壁温，因此管内壁温需要通过间接计算得到。假定管周向和轴向的导热是均匀的，可将本问题简化为内热源一维稳态导热，最终内壁温计算式如下：

式中：t为温度，℃；λ为金属导热系数，kW/(m·℃)；d为管外径，m。

2.1.2 实验参数与不确定度

实验主要从质量流速、流体压力以及热流密度3个方面来探索摩擦阻力的差异。此外，流体温度和壁面温度也是测量参数。其中质量流速、热流密度和摩擦阻力等属于间接测量值，压力和温度等属于直接测量值。考虑到数据采集卡存在识别精度，即使是直接测量值仍需做不确定度计算。表1给出了各参数在实验时的数值范围以及不确定度。

表1 实验参数范围与不确定度 Tab.1 Range and uncertainty of experimental parameters

各参数不确定度采用式(8)[16]计算：

2.2 神经网络模型

BP神经网络具有信号的前向传播和误差的反向传播2个过程。当输入信号经过隐含层、输出层获得的输出信号与期望输出不符时，将输出误差反向从输出层逐层传回输入层，每层每个神经元所获得的误差数据将作为调整其自身权值的依据，最终使得误差逐渐下降至目标范围或不再下降，训练完成。

本文BP神经网络模型采用MATLAB编程实现，设计了2层前馈神经网络结构（图2），可以实现对任意“一个有限空间到另一个有限空间的连续映射函数”的拟合[22]。更高的层数虽然理论上拟合能力更强，但存在过拟合以及训练难度增加的问题。隐含层神经元采用tansig非线性激活函数，输出层神经元采用purelin线性激活函数：

在常规BP神经网络中，一般采用神经网络响应的均方差作为目标函数，通过最优化算法实现均方差的逐渐减小，当均方差小于所设定的范围或不再继续下降时，认为神经网络已完成训练。然而，该方式易出现过拟合现象，使得神经网络泛化能力较差，无法达到较好的预测效果。与常规BP神经网络不同，贝叶斯正则化算法将目标函数修改为：

式中：α和β为正则化系数；Ew为神经网络中所有权值的均方差；Ed为常规神经网络的目标函数。

设计α和β的意义在于平衡神经网络拟合效果与泛化能力之间的矛盾。当α过大时，训练具有更好的泛化能力，但拟合精度不高，可能忽略了较多的数据特征；当β过大时，对训练集的响应均方差最小，但容易导致过拟合而对非训练集无法实现较好的预测效果。在贝叶斯理论中，假定训练集与权重集的先验概率服从高斯分布，最优的权值应具有最大后验概率，而最大后验概率在上述条件下可以等价于最小化F，求得F最小时的α和β为：

式中：γ为有效网络参数的数量；n为数据个数。该算法更适合训练具有噪声和数据量少的数据集[23]。由于阻力实验各参数均存在不确定度，且数据集规模对于可训练亿级数据的神经网络来说相对较小，因此采用贝叶斯正则化BP神经网络作为阻力预测模型，以提高模型泛化能力。

在模型中，将流体压力、质量流速、热流密度和主流温度4个参数作为输入参数，摩擦阻力系数作为输出参数。为提高模型预测精度，对输入参数和输出参数进行了归一化处理，归一化后各参数范围均在0~1。归一化方法如下：

实验总共包含194组数据集，模型训练集设定为165组，测试集为29组（随机选取）。隐含层神经元个数的选取采用经验公式和重复训练测试联合完成，常用公式为[24]：

式中：Nhidden为隐含层神经元个数；Nin、Nout分别为输入、输出参数个数；X为经验系数，一般取1~10。本模型中，输入参数有4个，输出参数为1个，因此推荐的隐含层神经元个数为3~12。为选择最佳隐含层神经元个数，进行重复训练测试，每种神经元个数方案重复训练100次（每次的测试集是随机的）并对100次训练的性能指标取平均作为最终值来对比不同神经元个数的效果。性能评价指标为Pearson相关系数R：

式中：x、y分别为模型计算值和实验测量值；n为数据个数。相关系数R的取值范围为0~1，其值越高，表明模型预测效果越好。

图3显示了不同隐含层神经元个数对相关系数R的影响。由图3可以发现，初始时随着隐含层神经元个数的增加，R逐渐增加，但当隐含层神经元个数超过11时，R几乎不再增加。因此，选取隐含层神经元个数为11，以确保拟合精度并避免过拟合。在该隐含层神经元个数下的100次重复训练中，选取R最高的训练神经网络作为后续计算分析使用。

3 结果与讨论

3.1 贝叶斯正则化BP神经网络与常规BP神经网络泛化能力对比

为验证贝叶斯正则化BP神经网络在预测S-CO2流动阻力时的泛化优势，采用新测试数据集（不在原165组训练集和29组测试集之内）进行性能测试。新测试数据集的参数范围见表2。

表2 新测试数据集的参数范围 Tab.2 Parameter range of new test data set

新测试数据集对于模型属于未知数据集，总共35组，其中21组为文献数据，剩余14组为本实验系统提供。此外，有27组（21组文献+6组本实验系统）为超过原始194组数据集参数范围限制的数据（即超限点）。

图4给出了贝叶斯正则化BP神经网络与常规BP神经网络对新测试数据集的预测结果对比。常规BP神经网络除未使用贝叶斯正则化算法外，其余设定以及训练均与贝叶斯正则化BP神经网络相同。图4中E±10%表示与实验值误差在10%以内的预测值个数占总个数的百分比，E±30%表示与实验值误差在30%以内的预测值个数占总个数的百分比。

由图4可以看出：贝叶斯正则化BP神经网络的相关系数、E±10%和E±30%等指标均优于常规神经网络；对于超限点，贝叶斯正则化BP神经网络同样具有较好的预测性能，而常规BP网络的误差则相对较大。综上，贝叶斯正则化BP神经网络相比常规BP神经网络具有更好的泛化能力，可满足工业应用要求。

3.2 贝叶斯正则化BP神经网络模型与典型摩擦 系数公式对比

在以往的研究中，非等温流体（管壁存在加热或冷却）经验公式主要是在等温流体公式的基础上进行修正，考虑管截面方向的物性差异。Filonenko公式[25]是目前应用最广的等温流体公式，在其后的很多研究者基于该公式进一步修正，获得非等温条件下的经验公式。Kirillov等人[10]给出了基于密度比修正的非等温流体公式，张海松等[12]给出了基于超临界沸腾数SBO和普朗特数比的非等温流体公式。Churchill[26]和Haaland[27]给出了另外2种常见的等温流体公式。为验证BP神经网络模型的实际预测效果，将其计算结果与表3中的5个典型摩擦系数公式计算结果进行对比，如图5所示。

表3 典型摩擦系数公式 Tab.3 Typical friction coefficient formula

由图5可知：当摩擦阻力较小时，各公式的预测效果均较好，但当摩擦阻力逐步增加后，预测值和实验值差距开始变大；在本实验数据范围内，3个等温流体公式的预测值基本一致，3个评价指标值也基本相同；5个典型公式计算值与摩擦阻力实验值的相关系数均在0.98~0.99，属于极高相关；在5个典型公式中，Kirillov等人的非等温流体公式预测效果较好，E±10%=28.35%，E±30%=86.08%；与5个典型公式相比，BP神经网络模型的预测效果最好，相关系数可达0.998 6，E±10%=98.97%，E±30%=100%。

3.3 不同边界条件对S-CO2摩擦阻力的影响

图6给出了不同质量流速G（765.56、1 445.56、2160.71 kg/(m2·s)）下S-CO2摩擦阻力实验值以及预测值。由图6可知，随着质量流速G的增加，S-CO2摩擦阻力增加，且摩擦阻力-主流温度曲线的斜率也逐渐增大。这是由于摩擦阻力与质量流速的平方G2呈正比关系。温度变化主要影响摩擦阻力系数和密度，将主流温度作为自变量时，其斜率主要受G2（正相关）、密度（负相关）以及摩擦阻力系数（正相关）的影响。因此，G增加，其斜率也逐渐增大。

图7给出了不同热流密度（208.56、279.14、335.10 kW/m2）下S-CO2摩擦阻力实验值以及预测值。由图7可知，随热流密度的增加，摩擦阻力下降。热流密度主要影响了管截面上的流体物性分布，即对于同一主流温度，热流密度主要影响摩擦阻力系数。随着热流密度的增加，根据实验现象，其摩擦阻力系数减小。因此，在BP神经网络模型中，考虑热流密度的影响是必要的。

图8给出了不同流体压力（10.47、15.53、19.52 MPa）下S-CO2摩擦阻力实验值以及预测值。由图8可以看出，随着压力的增加，摩擦阻力逐渐减小，且该现象在较低压力下更明显。根据摩擦阻力公式，当管长L、管径D和质量流速G不变时，摩擦阻力主要受摩擦阻力系数f（正相关）和密度ρ（负相关）影响。随压力的增加，S-CO2密度增加，黏度增加。密度增加会使S-CO2体积流量下降，流速下降，摩擦阻力下降；但黏度增加会导致流体间切应力增加，流体间能量损耗增大，摩擦阻力系数增加。综合来看，密度影响更关键（可能的原因是黏度对摩擦阻力系数的影响并非线性正比例关系且密度同样影响摩擦阻力系数），因此摩擦阻力下降。

4 结 论

1）S-CO2摩擦阻力随热流密度的增加而减小，即其摩擦阻力系数随热流密度的增加而减小，这有利于S-CO2循环发电的大规模工程应用。S-CO2摩擦阻力随流体压力的增加而减小，主要原因在于压力增加导致流体密度增加，而黏度增加导致的阻力增加不占主导地位。S-CO2摩擦阻力随质量流速的变化情况符合传统规律。

2）完成数据训练的贝叶斯正则化BP神经网络模型能够有效预测S-CO2的摩擦阻力，该方法泛化能力强、拟合度高且可以避免求解不同状态下的SCO2特征参数，简化了工程应用，为S-CO2阻力预测提供了新的思路。