2. 中国科学院大学 北京 100049;
3. 中国科学院近代物理研究所 兰州 730000
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000, China
钍基熔盐堆(Thorium Molten Salt Reactor,TMSR)核能系统项目是中国科学院先导科技专项之一，其战略性目标是研发第四代裂变反应堆核能系统[1]。随着熔盐堆系统研究的不断深入，熔盐堆系统分析的需求愈发凸显，现有的反应堆系统分析程序大多是针对压水堆系统建立的两相流一维流动模型无法直接应用到熔盐堆的系统模拟当中。需要一个单相一维多工质的系统模拟程序用于熔盐堆的功率升降、流量变化、负载变化等工况模拟，验证系统设计参数。
Aspen HYSYS是面向油气生产、气体处理和炼油工业的模拟、设计和性能监测的流程模拟工具[2]，同时具有稳态和动态两种仿真模式。其技术广泛应用于石油开采、储运、天然气加工、石油化工、精细化工、制药、炼制等领域，在世界范围内石油化工模拟、仿真技术领域占主导地位[3]。该软件具有丰富的用户自定义和软件通讯接口，具有良好的软件的扩展、二次开发和与其它软件的链接功能[4]。HYSYS具有强大的物性计算包和物性预测能力[5]，可以用作单相和多相一维多工质的系统模拟。但是HYSYS中缺少反应堆中子动力学计算模块，该模块是反应堆与化工系统的关键区别，要实现对反应堆的仿真，需要将中子计算模块加入到HYSYS中。HYSYS允许用户按照自身需求建立物性包、增加新的反应方程以及控制参数的变化等，具有开发成为一款满足反应堆系统仿真需求的高效系统瞬态仿真软件的条件和潜力。
本文对HYSYS进行二次开发，将中子动力学计算模块与HYSYS进行耦合，并对耦合程序进行对比验证，旨在使HYSYS程序能够对反应堆系统各种稳态和瞬态工况进行可靠、准确模拟，为反应堆设计和安全分析提供帮助。
1 TMSR-SF1介绍TMSR-SF1是10 MWth固态燃料钍基熔盐堆的简称，是中国科学院TMSR专项的一部分。它具有发电效率高、安全性好、核燃料有效利用和防止核扩散等优点[6]。
1.1 TMSR-SF1系统和设备TMSR-SF1包括如下主要系统和设备：堆本体和堆内构件、反应堆冷却剂系统及设备、安全相关系统、核辅助系统、测量和控制系统等，分为带走堆芯热量的一回路、中间换热的二回路和将热量带走的空冷侧三个部分。
堆本体由堆芯燃料球、石墨反射层、堆芯支撑机构、石墨反射层压紧机构、堆芯围筒、主容器、控制棒及控制棒驱动机构、堆内测量机构等组成。采用直径为6 cm的TRISO (Tri-structural iso-tropic)包覆颗粒燃料球作为燃料。
一回路主要由主循环泵、中间换热器(Intermediate Heat Exchanger,IHX)、溢流罐、管路和辅助设备组成，采用2LiF-BeF2盐作为冷却剂。二回路主要由循环泵、IHX二回路侧、熔盐-空气换热器、溢流罐、管路和辅助设备组成，采用FLiNaK盐作为冷却剂。
正常运行时一回路中熔盐冷却剂以强迫循环的形式从堆芯带出核裂变能，并通过双熔盐热交换器把热能从一回路熔盐转移到二回路熔盐；然后二回路熔盐中的热量通过空气散热器转移到大气环境。
1.2 TMSR-SF1系统主要参数
|
表 1 MSR-SF1主要热工水力参数 Table 1 Main thermal-hydraulic parameters of TMSR-SF1. |
该耦合程序以Aspen HYSYS和中子动力学计算程序为基础，模拟反应堆系统在不同运行工况下的响应情况，重点观测系统温度等参数随时间变化。在中子动力学计算中，采用6群中子动力学方程模型，并考虑燃料温度变化对堆芯温度反馈的影响，具有较高的计算精度和效率，可以用于各种堆型的系统仿真。
2.1 HYSYS计算模型采用HYSYS程序进行热工水力计算，在HYSYS中，所有单元操作均采用集总参数模型。这里主要对仿真建模中涉及的数学计算模型进行简要的介绍，而不再对用到的各个元件和设备的具体的计算原理进行详细的描述。守恒关系是HYSYS中数学计算的基础。下面将给出质量、组分、能量等的动态守恒关系式，并简述软件的求解策略。
2.1.1 物料守恒质量守恒的基本关系式为：
| $ \frac{\text{d}({{\rho }_{\text{o}}}V)}{\text{d}t}={{F}_{\text{i}}}{{\rho }_{\text{i}}}-{{F}_{\text{o}}}{{\rho }_{\text{o}}} $ | (1) |
式中：Fi、Fo分别为进、出口流量；ρi、ρo分别为进、出口的流体密度；V为体积。
式(1)只是Aspen HYSYS中所使用方程的简化形式，在软件中考虑了反应、蒸发密度变化等更复杂的情况。
2.1.2 组分守恒组分守恒的基本关系为：系统中某一组分的变化等于流入系统的该组分减去流出该系统的组分再加上该组分发生反应的量。
对于多组分物流，组分j的平衡关系式如下：
| $ \frac{\text{d}({{C}_{j\text{o}}}V)}{\text{d}t}={{F}_{\text{i}}}{{C}_{j\text{i}}}-{{F}_{\text{o}}}{{\rho }_{j\text{o}}}+{{R}_{j}}V $ | (2) |
式中：Cji为入口物流中组分j的浓度；Cjo为出口物流中j的浓度；Rj为组分j的反应率。
在反应堆中，中子的影响对物质的改变非常小，而堆中物质相容性非常好，化学反应可以忽略不计，认为发生反应的量为零。
2.1.3 能量守恒能量流入或流出系统是通过对流或者传导，通过边界向系统中增加能量是通过热传导或辐射。
对于有散热的连续反应器，可以采用如下方程：
| $ \begin{align} & \frac{\text{d}}{\text{d}t}[(u+k+\phi )V]={{F}_{\text{i}}}{{\rho }_{\text{i}}}({{u}_{\text{i}}}+{{k}_{\text{i}}}+{{\phi }_{\text{i}}})- \\ & _{{}}^{{}}{{F}_{\text{o}}}{{\rho }_{\text{o}}}({{u}_{\text{o}}}+{{k}_{\text{o}}}+{{\phi }_{\text{o}}})+ \\ & _{{}}^{{}}Q+{{Q}_{r}}-(w+{{F}_{o}}{{P}_{o}}-{{F}_{i}}{{P}_{i}}) \end{align} $ | (3) |
式中：u、k、Φ分别为单位质量的内能、动能和势能；w为系统单位时间内所做的轴功；Po为容器的压力；Pi为进料的压力；Q为通过边界增加的热量；Qr为反应热。
2.1.4 求解策略HYSYS中采用欧拉隐式求解方法，对体积（压力-流量）、能量和组分的求解并不同步，以减少计算机求解工作量，他们分别对应的默认求解频率分别为1、2、10个时间步。
2.2 中子计算模型 2.2.1 点堆模型考虑到计算的速度、精确性以及计算模型的复杂程度，选择广泛采用的具有6组缓发中子的反应堆中子动力学方程模型：
| $ \frac{\text{d}P(t)}{\text{d}t}={{k}_{\text{eff}}}\left( \frac{{{k}_{\text{eff}}}-1}{{{k}_{\text{eff}}}}-\beta \right)\frac{P(t)}{l}+\sum\limits_{i=1}^{6}{{{\lambda }_{i}}{{c}_{i}}(t)} $ | (4) |
| $ \frac{\text{d}{{c}_{i}}(t)}{\text{d}t}={{k}_{\text{eff}}}\frac{{{\beta }_{i}}}{l}P(t)-{{\lambda }_{i}}{{c}_{i}}(t)\ \ \ \ \ i=1,2,\cdots ,6 $ | (5) |
总反应性表示为：
| $ \begin{align} & \rho (t)={{\rho }_{0}}+{{\rho }_{\exp }}(t)+({{\alpha }_{f}}+{{\alpha }_{m}})({{T}_{f}}-{{T}_{f0}})+ \\ & _{{}}^{{}}{{\alpha }_{\text{r}}}({{T}_{\text{r}}}-{{T}_{\text{r}0}})+{{\alpha }_{\text{c}}}({{T}_{\text{c}}}-{{T}_{\text{c}0}}) \end{align} $ | (6) |
式中：ρ0为初始反应性；ρexp为控制机构动作等原因所引起的反应性变化，通过该项可以引入反应性瞬变；Tf、Tr、Tc和T6分别表示燃料、反射层、下降管和上升管中冷却剂按质量取的平均温度；αf、αm、αr和αc分别表示燃料、慢化剂、反射层、冷却剂的温度系数。
程序中采用gear算法求解中子动力学方程，该算法可用于求解刚性微分方程组。程序模型中将反射层、燃料、冷却剂等划分成多个节点，反射层、燃料、冷却剂等每个节点相关的量，在程序中默认为均匀划分。
2.2.2 耦合方法选择Automation方式来实现相应的耦合功能。将中子动力学计算程序改编成动态链接库文件(Dynamic Link Library,DLL)，以供其他程序在需要的时候调用，该DLL可以返回修正的堆芯功率。进行HYSYS中建立反应堆系统模型时，用加热器代替堆芯对整个系统进行加热，堆芯部分轴向分成多段，加载不同的反应堆功率，热量通过换热器传递到二回路中。读取t时刻HYSYS 中流经堆芯的温度T(i,t)及添加的反应性ρ(t)后，传递到DLL中；DLL中进行燃料的传热计算和堆芯燃料、流体的温度反馈，返回修正后的堆芯功率到HYSYS中，经重新分配作为下一时刻t+1的堆芯功率P(i,t+1)，再利用HYSYS进行热工水力计算。
3 仿真验证为验证仿真计算结果的准确性，本文针对熔盐堆系统进行仿真，并将仿真计算结果与RELAP5的计算结果进行比较。RELAP5程序是由美国爱达荷国家工程实验室开发的热工水力系统分析程序，被广泛应用于事故分析和安全评审等方面[8]。在空气冷却器流量增加和流经空气冷却器的冷空气温度升高两种情况下，对流经堆芯的熔盐进出口温度和反应堆的反应性进行分析，模拟时间均为2000 s。
建立的系统如图 1所示。初始时刻反应堆均处于满功率运行状态，堆芯额定功率为10 MW，HYSYS软件模拟的反应堆满功率运行的进出口温度分别为875.6 K、904.5 K。
|
图 1 HYSYS模拟系统图 Fig. 1 System diagram of HYSYS simulation. |
具体工况描述如下：在0时刻，熔盐空气冷却器空气侧流量增加了3%，并维持在新流量水平上运行，其他条件保持不变。空气流量的增加通过增大空气侧进口压力使进出口压差增大的方法实现。
HYSYS仿真结果与RELAP5仿真结果对比如图 2所示。无耦合情况下模拟堆芯功率没有发生变化。从图 2中可以看出，有无耦合对HYSYS模拟结果影响较大，无耦合情况下HYSYS模拟所得到的堆芯进出口温度均呈现一直下降的状态，而耦合后HYSYS模拟的堆芯进出口温度则有温度的短时间小幅上升情况出现，且无耦合模拟结果温度变化速率和幅度明显快于耦合后模拟结果，这是因为耦合后加入了负温度反馈特性，对温度变化幅度和速率均产生了抑制变化的效果，使得温度的降低速率和幅度变慢，逐渐趋于新的平衡状态。耦合后HYSYS模拟结果与RELAP5模拟结果的总体趋势是一致的：初始时均有小幅的温度降低再升高过程，然后温度缓慢降低。出现该现象的原因是：熔盐空气冷却器的流量增大使得该换热器的空气侧换热能力增强，从二回路带走的热量增多，从而导致流经该换热器的二回路FLiNaK熔盐温度降低；二回路与一回路通过双熔盐换热器进行热量传递，二回路熔盐温度的降低引起一回路流经反应堆堆芯的FLiBe熔盐温度降低；由于反应堆的负温度反馈特性，温度的降低引入正的反应性，堆芯功率增加，使得短时间内对一回路熔盐加热功率大于冷却功率，熔盐温度升高；后由于温度的变化量减小，正反应性逐渐减小，冷却功率大于加热功率，温度逐渐降低。
|
图 2 空气流量增加时堆芯进出口温度随时间的变化 Fig. 2 Temperature change of the core inlet and outlet vs. time with air flow rate increasing. |
但是由于初始时刻，HYSYS进出口温度较RELAP5略低，在整个模拟时间内，HYSYS模拟得到的堆芯进出口温度均低于RELAP5的进出口温度。而且图 2中HYSYS温度变化曲线斜率略大于RELAP5，即温度变化略快于RELAP5，HYSYS耦合的负温度反应性系数大于RELAP5的反应性系数。这种情况的出现与热工计算及温度权重等因素有关。
3.2 空气冷却器空气温度升高具体工况描述如下：在0时刻，熔盐空气冷却器空气侧温度增加了10 K，空气流量保持不变。分析空气进口温度变化对反应堆运行的影响，温度增加可以直接在HYSYS中通过设定进口温度实现。
该工况下，HYSYS仿真结果与RELAP5仿真结果对比如图 3所示，无耦合情况下模拟堆芯功率没有发生变化。从图 3中可以看出，有无耦合对HYSYS模拟结果影响较大，无耦合情况下HYSYS模拟所得到的堆芯进出口温度均呈现不断升高的状态，而耦合后HYSYS模拟的堆芯进出口温度则有温度的短时间小幅下降情况出现，且无耦合模拟结果温度变化速率和幅度明显快于耦合后模拟结果，这是因为耦合后加入了负温度反馈特性，对温度变化幅度和速率均产生了抑制变化的效果，使得温度的升高速率和幅度变慢，逐渐趋于新的平衡状态。耦合后HYSYS模拟结果与RELAP5模拟结果的总体趋势是一致的：两者总体变化趋势相同，均出现初始时刻温度短暂升高后降低、再不断升高的情况，且温度升高速度不断降低，温度变化曲线趋于平缓。出现该现象的原因为：在相同的空气流量情况下，突然升高熔盐-空气换热器进口温度，使得换热器的冷却能力降低；二次侧FliNaK熔盐温度升高，进而导致一回路总FliBe熔盐温度升高；由于该反应堆具有负温度系数，温度的升高引入负反应性，堆芯加热功率减小，初始时功率降低量大于冷却能力降低量，即发热量小于冷却量，故一回路熔盐温度降低，进而导致一回路温度升高量降低，引入反应性减小，温度降低速度变慢；温度降至最小值时，冷却能力小于堆芯功率，一回路温度身高，由于负反馈的影响，温度升高逐渐变慢，但加热功率仍占主导地位。
|
图 3 入口空气温度升高时堆芯进出口温度随时间的变化 Fig. 3 Temperature change of the core inlet and outlet vs. time with inlet temperature of air increasing. |
但在温度从最低值开始升高时，RELAP5的升温速率大于HYSYS升温速率，可以得出，HYSYS耦合的负温度反应性系数大于RELAP5的反应性系数，这与热工计算及温度权重等因素有关。
3.3 阶跃引入反应性通常情况下反应堆测量系统和反应性控制系统存在一定的误差，导致实际运行时可能出现棒位略高于临界棒位而引入少量反应性。这里就针对反应堆引入正反应性的情况进行仿真。假设反应堆处于额定功率运行状态，在0时刻突然引入5 pcm的反应性，分析反应性小幅度扰动对反应堆运行的影响。
引入5 pcm反应性后HYSYS耦合仿真结果与RELAP5仿真结果的对比如图 4所示。
|
图 4 引入正反应性时堆芯进出口温度随时间的变化 Fig. 4 Temperature change of the core inlet and outlet vs. time with additional positive reactivity introduced. |
图 4中两种模拟所得堆芯进出口温度曲线的总体趋势是相同的，总体上均呈现升高的趋势。引入正反应性之后，反应堆堆芯功率增加，使得熔盐的温度升高，由于该反应堆负温度系数的作用，温度增加的速率受到抑制，随着时间的增加，温度变化速率逐渐减小，温度趋近于新的稳定值。但是HYSYS模拟得到的温度变化比RELAP5的温度变化快。
以上几个工况计算结果的对比表明，HYSYS计算的结果与RELAP5的结果趋势相同，计算结果相差不大，证明了通过HYSYS耦合进行系统仿真的可行性。
4 结语基于熔盐堆系统与普通化工系统的不同，本文对HYSYS软件进行深度扩展，在加入熔盐物性和修正换热系数计算公式的基础上，将点堆模块与HYSYS软件进行耦合，引入堆芯反馈特性。通过软件模拟和RELAP5仿真计算的对比结果来看，扩展后的HYSYS 软件成功地模拟了熔盐堆系统的瞬态过程。软件计算的回路温度随时间的变化趋势、温度峰值大小等均与RELAP5计算结果符合较好。
HYSYS 软件能够实现对熔盐堆系统的瞬态过程比较准确地模拟，证明软件适用于反应堆回路瞬态仿真，体现了其强大的可二次开发能力，显示了较强的可靠性。在缺乏熔盐堆系统专用仿真软件的现状下，对HYSYS软件功能进行深度挖掘，开发出具有较完善功能和模块的熔盐堆系统仿真软件是熔盐堆研究的一种有益探索，是对熔盐堆系统仿真软件的一个有效补充，对现阶段熔盐堆的热工水力和安全研究具有重要的意义。
HYSYS软件的耦合不仅利用了HYSYS软件物性计算能力强大和模拟计算准确等优势，而且弥补了通用热流体计算程序不能实现反应性反馈计算的不足，还能避免许多热工分析程序无法实现不同物性冷却剂混合、流体接触空气等情况的缺点，可以进行一些特殊现象的模拟。后续可以利用HYSYS软件进行更多典型工况和特殊物理现象的仿真，开发熔盐物性计算方法，丰富和完善其在熔盐堆中的应用。
致谢 该工作的完成得到了中国科学院上海应用物理研究所TMSR中心反应堆物理部计算组职工何龙的帮助，在此表示衷心感谢。| [1] |
江绵恒, 徐洪杰, 戴志敏. 未来先进核裂变能——TMSR核能系统[J].
中国科学院院刊, 2012, 27 (3) : 366 –374.
JIANG Mianheng, XU Hongjie, DAI Zhimin. Advanced fission energy program-TMSR nuclear energy system[J]. Bulletin of Chinese Academy of Sciences, 2012, 27 (3) : 366 –374. DOI: 10.3969/j.issn.1000-3045.2012.03.016 (  0)
|
| [2] |
蔡飞云, 刘江楠, 张强. 流程模拟在焦化吸收稳定系统中的应用[J].
化工自动化及仪表, 2014, 41 (3) : 311 –314.
CAI Feiyun, LIU Jiangnan, ZHANG Qiang. Application of process simulation in stable system for coking absorption[J]. Control and Instruments in Chemical Industry, 2014, 41 (3) : 311 –314. DOI: 10.3969/j.issn.1000-3932.2014.03.019 ( 0)
|
| [3] |
金大天, 曹义鸣, 王丽娜, 等. HYSYS的二次开发及其在多组分气体膜分离过程模拟中的应用[J].
膜科学与技术, 2012, 32 (5) : 87 –91.
KIM Daechun, CAO Yiming, WANG Lina, et al. Secondary development of HYSYS and its application in the simulation of multi-component gas separation by membrane[J]. Membrane Science and Technology, 2012, 32 (5) : 87 –91. DOI: 10.3969/j.issn.1007-8924.2012.05.016 ( 0)
|
| [4] |
唐延泽, 郑海伟, 陈玉爽, 等. HYSYS软件应用于熔盐冷却系统瞬态分析的可行性研究[J].
核技术, 2015, 38 (4) : 040601 .
TANG Yanze, ZHENG Haiwei, CHEN Yushuang, et al. Feasibility study of HYSYS code for molten salt cooled system in dynamics analysis[J]. Nuclear Techniques, 2015, 38 (4) : 040601 . DOI: 10.11889/j.0253-3219.2015.hjs.38.040601 ( 0)
|
| [5] |
郭广智. 石油化工动态模拟软件HYSYS[J].
石油化工设计, 1997, 14 (3) : 29 –33.
GUO Guangzhi. Petrochemical dynamic simulation software HYSYS[J]. Petrochemical Design, 1997, 14 (3) : 29 –33. ( 0)
|
| [6] |
秋穗正, 张大林, 苏光辉, 等. 新概念熔盐堆的固有安全性及相关关键问题研究[J].
原子能科学技术, 2009, 43 (S1) : 64 –75.
QIU Suizheng, ZHANG Dalin, SU Guanghui, et al. Research on inherent safety and relative key issues of a molten salt reactor[J]. Atomic Energy Science and Technology, 2009, 43 (S1) : 64 –75. ( 0)
|
| [7] |
TMSR卓越中心. 10 MWt固态燃料钍基熔盐实验堆概念设计报告[R]. 上海:中国科学院上海应用物理研究所, 2014
Center for Excellence in TMSR. The conceptual design report on 10 MWt solid thorium-based molten salt experiment reactor[R]. Shanghai:Shanghai Institute of Applied Physics, Chinese Academy of Sciences, 2014 ( 0)
|
| [8] |
刘余, 张虹. RELAP5程序耦合接口的开发[J].
核动力工程, 2009, 30 (6) : 38 –40.
LIU Yu, ZHANG Hong. Development of RELAP5 coupling interface[J]. Nuclear Power Engineering, 2009, 30 (6) : 38 –40. ( 0)
|