陈 俊，罗 彬，吴照国，孔晓宁，林 萌,杨燕华,

(1.国家电投集团科学技术研究院有限公司， 国家能源核电软件重点实验室， 北京102209；2.上海交通大学核能科学与工程学院，上海200240)

COSINE子通道均相流分析程序的研发

陈 俊1，罗 彬1，吴照国1，孔晓宁1，林 萌2,杨燕华1,2

(1.国家电投集团科学技术研究院有限公司， 国家能源核电软件重点实验室， 北京102209；2.上海交通大学核能科学与工程学院，上海200240)

COSINE一体化软件包的子通道安全分析程序cosSubc基于子通道控制体三维网格模型，采用轴向及横向的热工水力控制方程，包括两流体和均相流两种求解算法。本文介绍了子通道均相流程序的物理模型和数值算法，并用cosSubc均相流程序和参考程序COBRA-TF分别对典型1000MW核电厂稳态算例进行计算分析，结果表明：cosSubc均相流程序与COBRA-TF吻合较好，具备堆芯子通道的热工水力计算能力。

COSINE；cosSubc；自主研发； 热工水力；均相流

核电设计软件用于核电厂设计、安全分析及运行管理，是支撑我国先进核电技术持续创新发展的关键重大共性技术，与CAP1400技术研发项目、CAP1700预研项目及压水堆重大共性技术项目目标的实现均有直接关系。核电设计软件的自主化是实现大型先进压水堆核电站重大专项总体目标的重要保证和必要条件。只有实现了核电设计软件的自主化，才能完全实现先进核电技术的自主化。基于此背景，国家核电技术有限公司北京软件技术中心启动了我国第一批自主开发的核电厂设计与安全分析软件。

COSINE(Core and System Integrated Engine for design and analysis，堆芯设计和安全分析一体化软件包)程序是我国第一批自主开发的核电厂设计与安全分析软件之一，具有堆芯核设计、堆芯热工设计、系统安全分析软件等功能。其中子通道安全分析程序命名为cosSubc，包括两流体模型程序和均相流模型程序，本文所涉及的为cosSubc均相流程序。

目前国内外已有子通道程序，部分基于均相流模型[1-3]，部分基于漂移流模型[4-7]，还有一些程序基于两流体模型[8,9]，另外部分学者对现有程序进行改进应用到不同堆型的程序开发与应用[10,11],还有学者应用现有的程序技术及改进分析[12,13]。其中大多数程序在严重瞬态或强烈扰动引起较大横流是预测结果较差，且不能给出控制容积内详细的两相流动参数分布[6]。为给出计算区域内每个位置的各状态参数的详细分布，需要发展一种适用性更广的计算方法。本文首先介绍cosSubc均相流程序的模型及算法；其次利用cosSubc均相流程序对1000MW核电厂稳态算例进行计算分析。

1 cosSubc均相流程序的模型与算法

cosSubc均相流程序基础算法基于两流体程序Newton-Raphson算法[14]，并在此基础上改进应用到均相流模型；该程序能模拟计算区域内任何一点的各状态参数。

本程序假设气液间达到热平衡和力平衡，气液两相压力相等；并假定存在两相时气液两相处于饱和状态，过冷过热状态通过相应公式修正。

1.1 守恒方程

1.1.1 质量守恒方程

cosSubc均相流程序的质量守恒方程主要包括混合相连续方程，不凝气体总连续方程，不凝气体组分连续方程组成：

第i种不凝气体的质量守恒方程:

(1)

总不凝气体质量守恒方程:

(2)

混合相的质量守恒方程可以展开如下：

(3)

式中：αg——气相体积份额；

ρnc,k——不凝气体k组分密度；

ρm——混合密度；

1.1.2 动量守恒方程

cosSubc均相流程序的动量守恒方程基于准三维网格，其动量守恒方程包括轴向流动方向上的动量守恒方程以及垂直于轴向的横向上的动量守恒方程。

混合相轴向动量守恒方程：

(4)

混合相正交横向动量守恒方程：

(5)

式中：下标i——速度流动的方向；

vm,j——j方向的横向速度；

τk,ij——k相j方向在i方向的作用应力；

1.1.3 能量守恒方程

cosSubc均相流程序的各相之间假设热平衡，则其各相之间的温度相等，气相和液相之间处于相同压力场下。为便于程序算法设计及数据处理，将能量方程的内能表达式转换为焓的表达式，相应的转换公式为：

u=h-pv

(6)

则其焓的能量方程可表示为混合相能量守恒方程：

(7)

式中：hm——混合相的焓；

p——总压；

Φm——其他源项。

能量方程中的dp/dt项表征了动能输运、压力作功对能量输运的影响，其中总能(e+uiui/2)方程中的动能项通过动量方程消去，考虑压力作功-，形成焓形式的能量方程中的dp/dt项。

1.2 数值算法

cosSubc均相流程序采用半隐式算法，方程计算过程如下：首先利用轴向和横向的混合动量守恒方程得到压力和流量之间关系式；其次利用Newton-Raphson算法[14]联立质量和能量守恒方程得到求解各主变量的变化率线性方程组；第三，利用各网格之间压力关系求解全场压力；最后更新其他变量及分析收敛。以此反复。

1.2.1 求解算法

求解算法的第一步为流量与压力的关系式。利用半隐式方法数值离散，横向和轴向的动量守恒方程可以表示为：

轴向动量守恒方程：

caxiWm,axi=aaxi+baxiΔp

(8)

横向动量守恒方程：

clatWm,lat=alat+blatΔp

(9)

式中：aaxi、alat——动量方程显式项；

baxi、blat——压力项系数；

caxi、clat——流量项系数。

其次，将动量方程带入混合质量守恒方程和混合能量守恒方程，利用改进的Newton-Raphson算法求解相应主变量。如果不存在不凝气体，线性方程组可以表示为：

Jwwδyw=-εw

(10)

其中:Jww为系数矩阵，其相应的表达式为：

(11)

(12)

如果存在不凝气体，线性方程组表示为：

(13)

式中：Jgg——不凝气体矩阵；

Jgw、Jwg——伴随矩阵；

δyg、δyw——解向量；

εw、εg——系统残差向量。

利用各控制体之间压降关系，结合质量守恒方程、能量守恒方程和动量守恒方程得到全场压力矩阵。单个控制体的压力方程可以表示为：

(14)

Newton-Raphson算法基于守恒方程的残差，其算法设计可适用于稳态和非稳态工况；在稳态计算时，其相应的稳定工况由守恒方程的迭代满足残差收敛为准。与普通的稳态算法相比，此稳态工况计算没有忽略时间项，其表达式严格满足守恒方程等式；但由于时间变化对工况稳定性有一定影响，其计算效率要稍低于普通的稳态算法。

1.2.2 稳定性分析

线性方程组计算时假设的初始值与最终解可能存在差异，在其计算过程中会存在收敛性和稳定性问题。堆芯子通道均相流程序中采用步长调整、多重循环迭代等方法提高系统的稳定性和收敛问题。

2 程序算例验证

通过对典型1000MW核电厂的堆芯稳态工况进行模拟分析，并与国际上知名的子通道程序COBRA-TF进行对比，验证cosSubc均相流程序。

2.1 初始条件及几何参数

根据对称性，将堆芯简化为八分之一模型；该模型有14条通道(channel)构成，每条通道与其他若干通道间由间隙(gap)相连，同时也与若干燃料棒相接触。模型网格节点如图1所示。

图1 模型网格示意图Fig.1 Grid of model(a) 通道与燃料棒;(b) 通道流动与间隙

图中CH1代表通道1，CHi代表第i个通道，BCin代表入口边界，BCout代表出口边界。

典型1000MW核电厂堆芯稳态算例的入口流量边界和出口压力边界条件如表1所示。

表1 堆芯计算边界条件及几何参数

计算考虑轴向壁面摩擦，定位格架局部阻力；通道间考虑横向局部阻力，湍流交混。主要模型如表2所示，COBRA-TF采用的模型参数与cosSubc相同。轴向不同位置有不同定位格架类型，不同定位格架的阻力系数不同。壁面传热定热流量指相应的燃料棒给定热功率。物性采用国际通用工业用水和水蒸气热力性质计算公式IAPWS-95。

表2 模型选择及参数

2.2 结果与分析

利用cosSubc对典型1000MW核电厂的八分之一堆芯稳态算例运行工况进行热工水力子通道分析，并将计算结果与COBRA-TF程序计算结果进行对比。

由于通道数目较多，选取靠近堆芯中心的通道1、通道2和通道3作为研究对象，并分别从沿程流量、压降以及温度分布进行分析。

图2给出了三个通道流量在轴向流动方向的变化趋势。其中红线为参考程序COBRA-TF计算结果，黑线为cosSubc计算结果。

对计算结果分析，发现沿流动方向，流量呈下降趋势；分析认为此三个通道均靠近堆芯中心，其流量下降因其为热通道缘故。

图2 通道1、通道2、通道3的沿程流量Fig.2 on-way flow of channel 1, channel 2 and channel 3(a) 通道1沿程流量；(b) 通道2沿程流量；(c) 通道3沿程流量

另外对比cosSubc和COBRA-TF计算结果，沿流动方向流量均呈锯齿形态。对比算例设置，发生脉冲上升点都处于定位格架所在位置。由于冷热通道的存在，流量会不断从热通道流向冷通道，而定位格架使各通道的流量通过局部阻力的形式再分配；在定位格架处，由于冷通道流量相对于热通道较大，根据动量守恒方程，流量与压降呈一定的线性关系，受到定位格架限制较明显，使在此处冷通道的流量通过压降的形式分流到热通道，使热通道压力脉冲式上升。由于流通面积因素，通道2的流量基本上等于通道1和通道3的两倍，详细见图1。

通过以上分析对比，cosSubc计算结果与参考程序COBRA-TF在热通道沿程流量计算方面吻合较好。

图3给出了三个通道压降在轴向流动方向的变化趋势。红线为COBRA-TF计算结果，黑线为cosSubc计算结果。

图3 通道1、通道2、通道3的沿程压降Fig.3 on-way pressure drop of channel 1, channel 2 and channel 3(a) 通道1沿程压降；(b) 通道2沿程压降；(c) 通道3沿程压降

从图中可以看出，cosSubc和COBRA-TF计算的三个通道压降均呈阶梯型下降趋势，吻合较好。对比定位格架和压降相应的位置，发现存在定位格架处，压降下降数值要大，而无定位格架点，压降下降较缓慢。这是由于无定位格架处，压降最主要是由于重力和沿程阻力引起的压降；当相应位置存在定位格架时会额外引入较大的局部阻力，导致压降下降较快；而通过定位格架点时定位格架引起的局部阻力消失，压降速率降低。

图4给出了三个通道温度在轴向流动方向的变化趋势。红线为COBRA-TF计算结果，黑线为cosSubc计算结果。

从图中可以看出，cosSubc计算的三个通道温度分布与参考程序COBRA-TF计算结果吻合较好。另外对所有通道分析，发现通道1、通道2和通道3的温升要高于其他通道；其中，通道3最高温度为336.04℃，为所有14个通道的最高值，通道2为335.83℃，通道1为335.69℃。从温升上看，最小DNBR值应该在通道3附近的燃料棒周围。

图5给出了所有燃料棒的最小DNBR。从不同燃料棒最小DNBR分布可以看出，越靠近堆芯中心，最小DNBR值越小；越偏离中心，最小DNBR越大。另外，cosSubc和参考程序COBRA-TF计算结果均表明最小DNBR出现在第四个燃料棒，cosSubc计算结果为2.58848，COBRA-TF计算结果为2.83212。对比图1，燃料棒4处于通道3周围，这与图4的温升曲线结果一致。

对比参考程序COBRA-TF和cosSubc，发现cosSubc计算值要整体小于COBRA-TF对应的参数值。根据W-3公式[15]，DNBR是临界热流密度与燃料棒上最大热流密度比值；由于燃料棒定功率曲线，COBRA-TF与cosSubc的最大热流密度相同，而临界热流密度与相应位置点的压力，流量，平衡含气率，等效水力直径，饱和焓以及通道入口焓同时决定。

COBRA-TF和cosSubc采用相同的初始条件和几何参数，其等效水力直径、通道入口焓两个程序相同；而影响临界热流密度主要参数有压力，流量，平衡含气率和饱和焓。根据W-3公式[15]，临界热流密度随压力增加而减小，随流量和焓的增加而增加，与平衡含气率关系不明确。

图4 通道1、通道2、通道3的沿程温度分布Fig.4 on-way temperature of channel 1, channel 2 and channel 3(a) 通道1沿程温度分布；(b) 通道2沿程温度分布；(c) 通道3沿程温度分布

图5 所有燃料棒的最小DNBRFig.5 minimum DNBR of all rods

对比图3，发现cosSubc的沿程压降要稍大于COBRA-TF；由于初始条件相同，两程序具有相同压力边界，压降越大其相应的压力越大，因此在相同位置cosSubc压力大于COBRA-TF，而压力越大其计算的临界热流密度会偏小。对比图2，cosSubc的通道流量要小于COBRA-TF；而临界热流密度与流量成正比，流量越大临界热流密度越大。

DNBR与临界热流密度成正比，临界热流密度越小，DNBR越小；而临界热流密度随压力增加而减小，随流量和焓的增加而增加。综上分析，cosSubc计算值要整体小于COBRA-TF对应的参数值。

cosSubc和参考程序COBRA-TF计算结果存在偏差主要是因为采取不同求解算法的缘故；cosSubc基于Newton-Raphson算法，采用混合焓，压力、流量作为主变量进行求解，并利用平衡含气率与主变量，结合相应的状态方程求解辅助变量，这与COBRA-TF算法存在较大差异。

3 结论

本文介绍了自主研发的具备自主知识产权的COSINE热工水力子通道程序cosSubc均相流模型的物理模型和数值算法。通过对典型1000MW核电厂稳态工况进行八分之一堆芯建模分析，并将计算结果与COBRA-TF程序结果进行对比，结果表明：cosSubc能够进行大型压水堆堆芯热工水力分析计算；相应的流量、压降和温度分布等与参考程序COBRA-TF吻合较好；

cosSubc均相流程序具备堆芯子通道的热工水力计算能力。cosSubc属于COSINE热工水力部分的子通道程序，是热工水力方向的核心程序之一，本文对cosSubc研发的思路和方法，为具备自主知识产权的COSINE热工水力程序的开发和应用奠定了基础。

致谢

感谢项目国家科技重大专项“核电关键设计软件评估”资助(2016ZX06004-002)。

[1] Rowe D S. COBRA IIIC: A Digital Computer Program for Steady-State and Transient Thermal-hydraulic Analysis of Rod Bundle Nuclear Fuel Elements[R]. BNWL-1695,1973.

[2] Stewart C W, Wheeler C L, et al. COBRA-IV: The Model and the Method[R]. BNWL-2214, Battelle Pacific Northwest Laboratiories, 1973.

[3] 白宁, 朱元兵, 任志豪, 等. 子通道分析程序LINDEN的开发与初步验证[J].原子能科学技术, 2013,47增刊：209-301.

[4] Shan J Q, Zhang B, Li C Y, et al. SCWR subchannel code ATHAS development and CANDU-SCWR analysis[J]. Nuclear Enginnering and Design, 2009, 239(10): 1979-1987.

[5] Toumi I, Bergeron A, Gallo D, et al. FLICA-4: a three-dimensional two-phase flow computer code with advanced numerical methods for nuclear applications[J]. Nuclear Enginnering and Design, 2000:139-155.

[6] 刘伟, 白宁, 朱元兵，等. 反应堆热工水力子通道分析程序ATHAS的研发[J]. 核科学与工程, 2014,34(1):59-66.

[7] 刘伟, 朱元兵, 白宁，等. 热工水力子通道分析程序ATHAS的稳态验证[J]. 核科学与工程, 2014,34(2):187-192.

[8] Joseph M K, Charles W S, Judith M C. VIPRE-02: A Two Fluid Thermal-hydraulics Code for Reactor Core and Vessel Analysis: Mathematical Modeling and Solution Methods[J]. Nuclear Technology, 1992,100.

[9] Thurgood M J, Kelly J M. COBRA TF: Model Description[R]. NRC advanced Code Review Group Meeting, Silver Spring, Maryland, 1980.

[10] 傅晟威, 许志红, 杨燕华. 基于SCWR堆芯结构的子通道程序开发与应用[J]. 原子能科学技术, 2011,45(3):345-350.

[11] 臧金光, 闫晓, 黄善仿, 曾小康, 黄彦平. 计算流体力学应用于子通道分析的初步方法研究[J]. 原子能科学技术, 2014,48(6):991-996.

[12] 毛文龙, 游松波, 俞冀阳. TACR燃料通道的子通道分析[J]. 核科学与工程, 2006,26(3):204-208.

[13] 史国宝, 唐家欢, 王杨定. 流道阻塞后的流场模型的开发及应用[J]. 核科学与工程, 2000,20(2):117-122.

[14] Cuervo D, Avramova M, Ivanov, K, Miro R. Evaluation and enhancement of COBRA-TF efficiency for LWR calculations[J]. Annals of Nuclear Energy, 2006, 33(9):837-847.

[15] 林诚格, 郁祖盛,欧阳予. 非能动安全先进压水堆核电技术[M]. 北京: 原子能出版社, 2010.

DevelopmentofHomogeneousSubchannelCodeofCOSINE

CHENJun1,LUOBin1,WUZhao-guo1,KONGXiao-ning1,LINMeng2,YANGYan-hua1,2

(1. State Power Investment Corporation Research Institute, National Energy Key Laboratory of Nuclear Power Software Beijing, 102209,China；2. State Nuclear Power Software Development Center, Shanghai Jiao Tong University,Shanghai, 200240, China)

COSINE is a self-reliant development code package of core and system integrated engine for design and analysis .Based on 3-dimension grids, subchannel code of COSINE, named cosSubc, is solved by axial and transverse thermal hydraulic governing equations,and the calculation is accurate,reliable and stable. The cosSubc code contains two solver algorithm, two fluid model and homogeneous model, and the homogeneous model’s algorithm in cosSubc is introduced in this paper. Then the COBRA-TF and cosSubc homogeneous code were used to calculate the 1000MW nuclear power plant at the steady-state conditions. The comparison results show that the two codes agree well with each other. The cosSubc homogeneous code is suitable for thermal hydraulics analysis of large pressurized water reactor.

COSINE; cosSubc; Self-reliant development; Thermal hydraulic; Homogeneous

2017-01-29

核电关键设计软件评估重大专项 (No.2016ZX06004-002)

陈 俊(1982—)，男，浙江义乌人，工程师，博士，从事核电热工水力方向研究

TL33

：A

：0258-0918(2017)04-0684-07