铅铋自然循环热工水力分析程序开发

2022-10-10程道喜杜开文齐晓光翟伟明马晓瑶张学锋

原子能科学技术 2022年9期

程道喜，杜开文，于松，齐晓光，翟伟明，周平，马晓瑶，张学锋

(中国原子能科学研究院，北京 102413)

铅基快堆(LFR)作为第4代核能系统选定的最具有发展前景的反应堆堆型之一，具有良好的增殖核燃料和嬗变核废料的潜力，并拥有突出的经济性和安全性。铅或铅铋作为反应堆冷却剂，拥有良好的自然循环特性，对于铅基反应堆事故工况下的余热导出，可以完全依靠一回路系统冷却剂的自然循环来实现；甚至小型铅基反应堆可以设计成纯自然循环反应堆。因此铅或铅铋的自然循环特性对于铅基反应堆的正常运行和事故工况余热导出均具有非常重要的意义。

目前针对铅铋自然循环的分析，国内外开展了基于已有程序的适用性评价以及程序开发。如瑞典皇家理工学院Ma等[1-2]开展了基于铅冷热工水力回路TALL的自然循环特性分析实验，模拟了瞬态事故，包括失热阱事故，一、二次失流事故，超功率以及启动、停止的瞬态过程等，并利用TRAC/AAA程序等进行了相应的数值模拟分析，验证了程序对于铅铋系统的模拟能力；意大利ENEA的Agostini等[3]在CHEOPE装置上开展了自然循环实验并利用RELAP5 MOD3.2开展了相应的分析，取得了较好的结果。Brogohain等[4]开发了铅铋自然循环程序，并利用矩形回路自然循环实验结果进行了初步验证。Zhao等[5]开发了纯自然循环铅冷反应堆的稳态热工水力分析程序，并与MARS程序开展了初步对比验证；Fiori和Zhou[6]验证了RELAP5/SCDAP对于铅铋热工水力行为的预测能力。陆道纲等[7]基于MATLAB Simulink编制了含有散热损失模型的铅铋自然循环回路计算程序，并利用实验结果进行了验证，结果表明铅铋回路散热会对分析结果产生较大的影响。

目前，对于铅或铅铋作为冷却剂的自然循环特性尚未得到充分认识，上述针对铅铋系统的分析程序最初均是针对水等介质开发，对于铅铋系统的适用性仍需进行深入研究及验证。为更好支持铅铋反应堆自然循环热工水力分析工作，本文在研究国内外快堆系统分析程序模型的基础上，建立堆芯、热交换器等一回路关键部件的热工水力模型，开发集成铅铋工质物性参数及流动换热关系式的自然循环瞬态分析程序，并基于铅铋回路内自然循环实验数据对程序进行初步验证。

1 数学物理模型

本程序主要开展一回路系统的自然循环热工水力分析，不包括中子模块，并且将二回路视为边界条件，给定换热器二次侧入口流量及温度。

1.1 一回路水力模型

闭合回路内自然循环的水力计算，需在整个冷却剂流动的闭合回路上对动量方程进行积分。每个控制单元动量方程为：

(1)

其中：Δz为控制单元长度；a为流动面积；W为质量流量；t为时间；ρ为流体密度；p为压力；g为重力加速度；θ为方向角度；f为摩擦阻力系数；De为当量直径；K为局部阻力系数；i为控制单元号。在整个冷却剂闭合环路上积分可得：

(2)

其中：j为环路上设备或部件号；ji代表第j个部件的第i个控制单元；Δppump为泵压头。

程序中泵模型采用离心泵。泵压头可通过额定转速下的泵特性曲线H0=f(W0)，结合泵转速与扬程、流量的关系计算得到，如式(4)所示。

(3)

(4)

其中：H为泵扬程；n为转速；下标0代表额定参数。

1.2 堆芯热工水力模型

堆芯采用平行多通道模型，每个通道代表1类组件。燃料元件内热量通过导热方式进行传递，包壳与组件内冷却剂之间通过对流换热进行热量传递。柱坐标下，燃料元件导热方程为：

(5)

堆芯组件内冷却剂能量方程为：

(6)

组件盒壁的能量方程为：

(7)

其中：下标b、f、s分别代表燃料元件、组件盒内冷却剂以及组件盒壁；T为温度；cb为燃料或包壳的比热；cs为组件盒壁的比热；r为径向坐标；z为轴向坐标；λb为燃料或包壳导热系数；λf为冷却剂导热系数；qbv为体积释热率；h为冷却剂比焓；hfb、afb分别为元件棒与冷却剂的换热系数和换热面积；hsf、asf分别为组件盒内壁与冷却剂的换热系数和换热面积；Φs为组件盒壁与外部环境的热交换量，即：

(8)

其中：hsa、asa分别为组件盒壁与堆内冷却剂的换热系数和换热面积；R为组件盒壁外部与环境大气的热阻，根据外壁空气的自然对流散热来计算。

1.3 热交换器模型

图1 热交换器模型节点划分Fig.1 Nodal diagram for heat exchanger

热交换器为管壳式换热器，可采用单管模型来进行模拟。热交换器径向上可划分为4个区域：一次侧流体、换热管、二次侧流体以及换热器壳体。热交换器模型节点划分如图1所示。考虑对流换热以及流体轴向导热，忽略换热管以及壳体轴向导热，分别列出4个区域的控制方程：

一次侧流体：

(9)

二次侧流体：

(10)

(11)

壳体：

(12)

其中：下标p、se、t、sh分别代表一次侧流体、二次侧流体、换热管以及壳体；ct为换热管壁固体比热；csh为壳体固体比热；mp为一次侧流体质量流量；mse为二次侧流体质量流量；a对于流体为流通面积，对于换热管或壳体为横截面积；hpt、apt分别为一次侧流体与换热管换热系数及面积；hst、ast分别为二次侧流体与换热管之间的换热系数及面积；hss、ass分别为二次侧流体与壳体之间的换热系数及面积；λ为导热系数；Φ为壳体与外部的热交换量，对于“浸泡”在冷却剂内部的热交换器，如池式堆中间热交换器(IHX)等，Φ为壳体与堆池内冷却剂的热交换量，对于外部热交换器，Φ为壳体(添加保温层)与环境的热交换器量，表达式与式(8)相同。对于圆形的热交换器壳体，壳体外部与环境大气的热阻R的计算方式为：

(13)

其中：dsh1、dsh2、λsh分别为壳体内径、外径和导热系数；dh1、dh2、λsh分别为保温层(若有)内径、外径和导热系数；ha为保温层外表面与大气的自然对流换热系数。

1.4 管道模型

管道模型中考虑流体轴向导热以及管道外壁与周围介质的热交换，忽略管壁轴向导热。管道外壁与周围介质的热交换同样采用式(8)中的方法。

1.5 流动及换热模型

本程序中采用的铅铋换热关系式列于表1，圆管及环隙流动阻力关系式以及典型结构局部阻力采用Idelchik阻力手册关系式[8]，棒束阻力采用Cheng-Todreas关系式[9]。

表1 铅铋换热关系式Table 1 Lead-bismuth heat transfer correlation

2 求解方法及程序结构

2.1 求解方法

为保证稳定性，动量方程及各部件的能量方程采用半隐格式进行离散。考虑到动量方程的求解与整个系统温度分布有关，能量方程的求解与质量流量有关，需采用合适的方法来对动量方程和能量方程解耦。在自然循环建立到稳定的过程中，回路内冷却剂的热工水力参数变化较为缓慢，同时在研究的温度范围内，各物性参数变化率不大，因此对于此类缓慢变化的物理过程，可通过先求解动量方程，再利用求解动量方程的结果来求解能量方程，并循环迭代直至满足收敛准则的方法来获得最终解。本程序中求解步骤如下：

1) 利用前一时间步的温度分布、流量求解动量方程，得到环路内冷却剂的流量；

2) 根据求解动量方程获得的质量流量，从堆芯开始，依次求解各部件的能量方程，获得温度分布；

3) 根据新的温度分布重新求解动量方程，得到冷却剂流量，并再次依次求解能量方程，循环迭代多次直至满足收敛标准(温度、流量相对变化小于限值)，最终获得的当前时间步的各部件温度分布及流量等；

目前地方高校科研成果专利保护和转化不足的主要原因之一就是资金不足。由于专利转化的成本高，特别是我国高等院校资金主要源于国家拨款，国家对科研的支持是当前高校科研的主要动力。缺乏良好的资本市场和金融环境，严重阻碍了科技成果转化的速度与水平。相比之下，日本不仅在全国各类大学中建立知识产权管理本部，并给予资金上的大力支持，同时成立了区域知识产权战略总部，在产学研合作开发项目给予资金上的支持。这些资金扶持的举措积极推进了大学的创新能力。缺乏风险投资、足量投入，科技成果转化市场的有序良性发展难以得到保障和维系。

4) 对于各时刻重复1～3过程直至结束时刻。瞬态计算求解流程如图2所示。

图2 瞬态计算求解流程Fig.2 Flowchart of transient solving

2.2 程序构架

本程序采用模块化结构，包括物性模块、输入输出模块、部件模块、关系式模块、求解模块等，程序结构如图3所示。程序采用Fortran95语言编制。

图3 程序结构Fig.3 Structure of code

3 程序初步验证

图4 铅铋自然循环实验装置Fig.4 Experimental device for natural circulation of lead-bismuth

本程序采用铅铋回路内自然循环实验数据进行初步验证。中国原子能科学研究院开展了铅铋自然循环实验研究工作，获得了不同功率启动自然循环以及功率台阶变化对于铅铋自然循环影响的实验数据。该实验装置主要由矩形闭合铅铋回路以及导热油冷却回路等组成，如图4所示。模拟组件采用电加热元件模拟堆芯燃料元件发热，热交换器为列管式换热器，一次侧(管侧)为铅铋，二次侧(壳侧)为导热油S750(氢化三联苯)，最终热阱为空冷器。

回路控制单元划分如图5所示。加热器径向控制体划分时，电热元件内部划分为3层(最内侧为发热丝，其余两层为MgO层)，包壳1层。加热器轴向加热区和非加热区通过输入参数中给定加热丝轴向体积释热率分布来实现模拟；热交换器径向分为一次侧、二次侧、换热管以及外壁。对于所有回路设备及管道外壁面散热考虑壳体、管道与外部空气通过自然对流方式来模拟，自然对流的传热关系式采用文献[14]中的关系式。热交换器二次侧导热油物性见生产商物性参数表[15]。

3.1 自然循环启动瞬态分析

在铅铋自然循环启动实验中，铅铋回路内工质从初始静止等温状态开始，瞬时启动模拟组件电加热和热交换器二次侧导热油循环泵，维持电加热功率、二次侧流量以及二次侧入口温度不变。图6为自然循环启动过程瞬态温度和流量的变化。

从图6可知，在电加热功率投入后，回路内铅铋自然循环流量迅速上升，100 s左右达到峰值流量，同时，电加热组件出口温度也同样迅速上升达到峰值，温度峰值时间稍早于流量峰值。

图5 回路控制单元划分Fig.5 Discretization of loop division of control units of loop

图6 自然循环启动过程瞬态温度和流量的变化Fig.6 Transient temperature and flow rate in natural circulation startup

在相同的二次侧边界条件下，电加热功率越大，组件峰值温度越高，最后形成的稳定自然循环流量也越大，冷热段温差也越大。由于自然循环流动较弱，瞬态过程中，电加热组件出口流体到达热交换器入口需要一定时间，因此热交换器入口温度的响应较组件出口要相对滞后；同样，热交换器出口流体到达电加热组件入口也需要一定时间，电加热组件入口温度下降较热交换器出口滞后。在散热较小的情况下，当稳定的自然循环建立后，冷却剂温度分布也逐步达到稳态，最终组件出口与换热器进口，组件入口与换热器出口之间的温度差异变得很小。

从图6可知，本程序计算结果能较为正确地反映铅铋自然循环回路内关键位置冷却剂温度的瞬态变化规律及相对关系，温度计算值与实验值偏差小于5%，流量计算值与实验值的偏差在自然循环建立的最初阶段稍大，其余时刻均小于10%。在实验中，由于不可避免的热损失以及热交换器向二次侧的导热，回路内初始存在温度差，从而使得自然循环建立速度变化。这一初始温度差在计算中不能准确反映，导致计算结果在自然循环建立阶段的偏差稍大。

3.2 功率台阶变化影响分析

功率台阶变化影响实验中，铅铋工质在一定电加热功率下建立稳定的自然循环后，不断改变电加热功率台阶，铅铋回路内工质的自然循环流动将随着功率台阶的不同而发生变化，组件进出口以及热交换器的温度也将响应功率的变化。如图7所示为功率台阶变化瞬态实验中，组件电加热功率随时间的变化，功率台阶分别为5、8和12 kW。

图7 组件功率台阶变化Fig.7 Step change of heater power

图8为功率台阶变化时瞬态温度和自然循环流量响应。随着功率以台阶形式变化，温度与流量有相应的响应。当组件电加热功率逐步上升时，自然循环流量增加，组件出口温度先于换热器出口温度上升；当组件电加热功率逐步下降时，自然循环流量降低，组件出口温度又先于换热器出口温度降低。对比结果可知，程序计算得到的结果与实验结果符合较好，能正确预测功率台阶变化时，组件进出、换热器进出口等关键位置的温度响应情况以及回路内自然循环流量的瞬态变化规律，其中温度相对偏差小于5%，自然循环流量相对偏差小于10%，进一步验证了本程序能较好模拟铅铋自然循环瞬态过程。