Development and Application of Core pin-by-pin Calculation of Sub-channel Analysis Code LINDEN
-
摘要: 为实现子通道分析软件LINDEN的pin-by-pin计算功能，提高LINDEN软件开展堆芯热工水力分析的计算精度。本研究开发了全堆芯pin-by-pin计算自动建模功能，解决了大型稀疏矩阵的存储和求解问题，开展了157组件反应堆全堆芯pin-by-pin计算及分析。计算结果表明，LINDEN软件可以开展全堆芯pin-by-pin计算，获得堆芯内冷却剂温度分布、偏离泡核沸腾比（DNBR）、质量流速分布等关键热工水力参数的详细分布结果。相邻组件的功率差异会影响组件内的热工水力参数分布。
-
关键词:
- 子通道分析 /
- 自动建模 /
- 大型稀疏矩阵 /
- pin-by-pin
Abstract: To realize the pin by pin calculation function of the sub-channel analysis code LINDEN and improve the accuracy of core thermal hydraulic analysis by LINDEN code, an automatic modeling function has been developed for the full core pin-by-pin calculation, and the storage and solution problems of large sparse matrices have been solved. The full core pin-by-pin calculation and analysis of 157 assemblies were carried out. The calculation results indicate that LINDEN can perform full core pin-by-pin calculations, and the detailed distribution results of key thermal hydraulic parameters such as coolant temperature distribution, DNBR, mass flow rate distribution can be obtained. The power difference between adjacent assemblies affects the distribution of thermal hydraulic parameters within the assembly.-
Key words:
- Sub-channel analysis /
- Automatic modeling /
- Large sparse matrix /
- pin-by-pin
-
0. 引 言
在反应堆堆芯热工水力分析中，通常采用子通道分析软件来进行堆芯关键热工水力参数的计算[1]。受限于计算机速度，传统的子通道分析软件在开展堆芯热工水力分析时，会对堆芯进行简化建模，无法直接获得堆芯内详细真实的热工水力参数分布。开展全堆芯pin-by-pin热工水力计算，能够真实精确地完成堆芯几何建模，提高反应堆堆芯热工水力计算精度。
得益于计算机技术的发展，开展全堆芯pin-by-pin尺度的分析计算成为可能。反应堆堆芯内存在复杂的多物理场耦合特性，开展高保真热工水力计算是其中的重要组成部分。宾夕法尼亚州立大学开发的压水堆堆芯热工水力分析软件COBRA-TF[2]具备全堆芯pin-by-pin计算功能。基于该软件，Avramova[3]等在NURESIM平台上完成了子通道程序CTF及中子物理程序DYN3D的集成和耦合，并开展了弹棒事故分析。Magedanz[4]等利用耦合程序TORT-TD/CTF开展了沸水堆燃料组件的高保真物理热工耦合计算。Bennett[5]等开展了蒙特卡罗中子物理程序MCNP6和子通道程序CTF的耦合研究工作，并对该耦合程序进行了测试和应用。GUO[6]等实现了蒙特卡罗程序RMC与子通道程序CTF的耦合计算，开展了压水堆全堆芯物理热工耦合计算，分析了中子数及硼浓度等参数的影响。刘宙宇[7]等人开展了高保真中子学和物理-热工耦合计算方法研究，研制了物理-热工耦合计算程序NECP-X/CTF，分析了燃料棒导热模型及间隙导热系数等对计算结果的影响。史敦福[8]等人开展了蒙卡中子输运程序JMCT和子通道程序COBRA-EN的耦合计算研究，由于COBRA-EN不具备开展pin-by-pin尺度的计算功能，耦合程序开展了对反应堆组件的模拟计算，论证了耦合计算在反应堆分析中的重要性。
子通道程序的pin-by-pin计算功能对于提升堆芯关键热工水力参数的精确性具有重要意义，也是开展高保真物理热工耦合计算的基础。本文基于子通道分析软件LINDEN，开发了全堆芯pin-by-pin自动建模功能，解决了大型稀疏矩阵的存储和求解问题，实现了LINDEN软件的pin-by-pin计算功能，并对157组件反应堆堆芯进行了pin-by-pin计算和分析。
1. LINDEN软件介绍
LINDEN软件采用C++语言编制[9]，可用于计算正常和事故工况下堆芯内冷却剂的流量、压力、比焓值和含汽率等热工水力参数。LINDEN软件基于子通道模型开发，可以考虑相邻通道间的质量、动量和能量的交换与迁移。软件采用四方程的漂移流模型，包括混合质量守恒方程、混合能量守恒方程、混合动量守恒方程以及液相能量守恒方程。其中，液相能量守恒方程用于模拟过冷沸腾过程中液相的热不平衡态。软件通过对计算对象划分子通道控制体，使用空间均匀化的平均值来定义每个控制体内的参数，从而建立起守恒方程组，并在特定的边界条件下对方程组进行求解，得到各个控制体的热工水力参数。最后，根据计算得到了局部参数，利用临界热流密度（CHF）关系式计算堆芯内的偏离泡核沸腾比（DNBR）分布。
1.1 质量守恒方程
质量守恒方程需要考虑相邻通道间的横向流动。LINDEN软件采用的质量守恒方程为：
$$ \frac{{\partial {G_i}}}{{\partial Z}} + \frac{1}{{{S_i}}}\sum\limits_{k \leftrightarrow i} {{l_{ik}}{g_{ik}}{\delta _{ik}}} = 0 $$ (1) 式中，$ G $为轴向冷却剂质量流速；$ l $为相邻通道间隙；$ g $为横向冷却剂质量流速；Z为轴向高度；δ为横向流量方向的标识；S为间隙宽度；下标i为当前通道编号；下标k为相邻通道编号。
1.2 混合动量守恒方程
混合动量守恒方程需要考虑轴向和横向2个方向的动量守恒，包括轴向动量守恒方程和横向动量守恒方程。
轴向动量守恒方程为：
$$\begin{split} &\frac{{\partial \left( {V'G_i^2} \right)}}{{\partial Z}} + \frac{1}{{{S_i}}}\sum\limits_{k \leftrightarrow i} {{l_{ik}}\left[ {V{{_{ik,{\mathrm{a}}}^{^{{\prime}}}}}G_{ik,{\mathrm{a}}}{g_{ik}}{\delta _{ik}} - } \right.}\\ &\left. {{{\left( {\mu + {\mu _{\text{t}}}} \right)}_{ik}}\frac{{\left( {V_k^{^{{\prime}}}{G_k} - V_i^{^{{\prime}}}{G_i}} \right)}}{{{L_{ik}}}}} \right] \\ &=- \frac{{\partial {p_i}}}{{\partial Z}} - \left( {\frac{{{f_i}\overline {{G_i}} {G_i}}}{{2{\rho _i}{D_{{\text{h}}i}}}} + \frac{{{K_i}\left| {{G_i}} \right|{G_i}}}{{2{\rho _i}}}} \right) - {\rho _i}{g_{\mathrm{a}}} \end{split} $$ (2) 式中，$ V' $为冷却剂比体积；下标$ {\mathrm{a}} $表示轴向方向；$ L $为相邻通道中心距；$ \mu $为层流扩散因子；$ {\mu _{\text{t}}} $为湍流扩散因子；下标t表示湍流；$ p $为冷却剂压力；$ \rho $为冷却剂密度；$ {g_a} $为重力加速度；$ f $为阻力系数；$ K $为格架阻力系数；$ {D_{\text{h}}} $为水力直径；Z为轴向高度。
横向动量守恒方程为：
$$ - \frac{{\partial P}}{{\partial {y_{ik}}}} = \frac{{{f_i}\overline {{g_{ik}}} {g_{ik}}}}{{2\rho {D_{\text{h}}}}} + \frac{{\partial \left( {V_{ik}'{G_{ik}}{g_{ik}}} \right)}}{{\partial Z}} + \frac{{\partial \left( {V_{ik}'{g_{ik}}{g_{ik}}} \right)}}{{\partial {y_{ik}}}} $$ (3) 式中，$ y $为横向距离。
1.3 能量守恒方程
能量守恒方程为：
$$ \frac{{\partial \left( {{G_i}{H_i}{S_i}} \right)}}{{\partial Z}} = \frac{{\partial {G_i}}}{{\partial Z}}{S_i}{H_i} + {G_i}{S_i}\frac{{\partial {H_i}}}{{\partial Z}} $$ (4) 式中，$ {H_i} $为冷却剂焓值。
1.4 液相能量守恒方程
能量守恒方程为：
$$\begin{split} &(1-x)\left\{\left({G}_{i}\frac{\partial {H}_{\text{l}i}}{\partial Z}\right)+\frac{1}{{S}_{i}}\left[{\displaystyle \sum _{k\leftrightarrow i}{l}_{ik}({H}_{\text{l}ik}-{H}_{\text{l}i}){g}_{ik}{\delta }_{ik}}+\right.\right. \\ &\left. {\left. {\sum\limits_{k \leftrightarrow i} {\frac{{{l_{ik}}}}{{{L_{ik}}}}\left( {k - {k_t}} \right)\left( {{H_{li}} - {H_{lk}}} \right)} } \right]} \right\} ={q_{{\text{cond}},i}} + \frac{{{P_{{\text{h}}i}} \cdot {q_{{\text{wl}}}}}}{{{S_i}}} + \frac{{{q_i}}}{{{S_i}}} \end{split} $$ (5) 式中，$ k $为热扩散系数；$ {q_{{\text{cond}}}} $为冷凝引起的热量交换；$ {P_{{\text{h}}i}} $为通道加热周长；$ {q_{{\text{wl}}}} $为燃料棒表面线功率；$ {P_{\text{r}}} $为通道加热长度；$ {q_i} $为辐照引起的热量交换；下标l表示液相。
2. 自动建模功能开发
开展全堆芯pin-by-pin计算时，子通道数量和通道间耦合关系数量将达到数万至数十万的级别，采用人工建模的方法会耗费大量的时间，也会增加建模过程中出现人为错误的概率。结合子通道程序的使用方法和建模方式，开发适用于子通道程序的自动建模功能，实现全堆芯pin-by-pin计算的建模功能，是开展全堆芯pin-by-pin计算的前提条件。
在进行子通道程序建模时，需要考虑通道之间的横向流动。为此，在进行子通道编号之后，需要进行通道耦合关系的定义。本文以157组件为例，介绍自动建模功能的开发过程。
图1给出了157组件反应堆堆芯燃料组件编号情况，图2给出了每个燃料组件内的子通道编号情况。对于157组件，共有50868个子通道。从图1和图2可以看到，堆芯内的燃料组件和燃料组件内的子通道均采用了连续编号。在计算子通道几何参数（流通面积A、润湿周长$ {P_{\text{W}}} $、加热长度$ {L_{\text{H}}} $）时，采用了以燃料棒/导向管为基准的面积份额方法。如图3所示，组件内考虑了燃料棒和导向管2种结构，对燃料棒/导向管采用连续编号，每根燃料棒/导向管与4个子通道相关联，通过确定燃料棒/导向管与子通道的份额关系，计算子通道的几何参数。
$$ A=\displaystyle\sum_n^{ }\alpha_i\left(P^2-\frac{\text{π}d_i^2}{4}\right)\quad\ \ \ \ \ i=1,2,...,n $$ (6) $$ P\mathrm{_W}=\sum\limits_n^{ }\alpha_i\text{π}d_i\text{ }\;\; \ \ \ \ \ \qquad\quad i=1,2,...,n $$ (7) $$ L\mathrm{_H}=\sum\limits_n^{ }\alpha_i\text{π}d_i\text{ }\;\ \ \ \ \ \qquad\, \; \quad i=1,2,...,n $$ (8) 式中，$ P $为燃料棒中心距；$ d_i $为燃料棒/导向管直径；$ {\alpha _i} $为子通道内燃料棒/导向管的份额；$ n $为与子通道相邻的燃料棒/导向管总数。
子通道间的相邻关系包括通道编号、通道中心距、通道间隙等，除了考虑组件内子通道的相邻关系，还考虑了组件之间子通道的相邻关系。
3. 数据存储与求解算法
在进行动量守恒方程和能量守恒方程的求解过程中，需要建立通道数×通道数的系数矩阵。当通道数量达到数万级别时，无法直接创建对应通道数量的二维数组，这是小规模计算时不会出现的问题。为了解决系数矩阵的存储问题，考虑到动量方程/能量方程系数矩阵是大型稀疏矩阵的特点，采用三元组法对系数矩阵进行计算和存储。
对于n×n阶的大型稀疏矩阵，假设其中含有m个非零项。用3个一维数组来存放该系数矩阵，即：LINE(I)——存入非零项I的行数；COL(I)——存入非零项I的列数；Y(I)——存入非零项的数值。
通过采用该方法进行稀疏矩阵的存储，可以将二维数组转换为3个一维数组，避免了数组太大造成的内存超限问题。图4中给出了3×3阶矩阵的存储示例，A表示系数矩阵；$ {a_{i,j}} $为系数矩阵系数。
LINDEN软件的求解流程如图5所示。首先，通过读取输入文件确定计算对象；然后，完成初始化；依次求解能量守恒方程、横向动量守恒方程、轴向动量守恒方程以及质量守恒方程，完成堆芯冷却剂流量分配和温度分布的计算；通过多次迭代计算，确保计算结果满足收敛条件；在获得堆芯内的热工参数分布之后，结合CHF关系式，完成堆芯DNBR分布的计算。
在进行能量方程和动量方程求解时，需要求解关于压力和焓值的大型稀疏矩阵，为了减小计算量，只针对系数矩阵中的非零项进行存储和计算，并采用高斯-赛德尔迭代法[10]进行求解，求解过程如下：
$$ \begin{split} & \left( \begin{array}{*{20}{c}}a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,n-2} & a_{1,n-1} & a_{1,n} \\ a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,n-2} & a_{2,n-1} & a_{2,n} \\ & \vdots & & & & \vdots & \\ a_{n-1,1} & a_{n-1,2} & a_{n-1,3} & \cdots & a_{n-1,n-2} & a_{n-1,n-1} & a_{n-1,n} \\ a_{n,1} & a_{n,2} & a_{n,3} & \cdots & a_{n,n-2} & a_{n,n-1} & a_{n,n}\end{array} \right)\times \\ & \left(\begin{array}{l}x_1 \\ x_2 \\ \vdots \\ x_{n-1} \\ x_n\end{array}\right)=\left(\begin{array}{l}b_1 \\ b_2 \\ \vdots \\ b_{n-1} \\ b_n\end{array}\right)x_i^{(k+1)} \\= & \frac{1}{a_{i,i}}\left[b_i-\sum\limits_{j=1}^{i-1}a_{i,j}x_j^{(k+1)}-\sum\limits_{j=1}^{i-1}a_{i,j}x_j^{(k)}\right]\\[-1pt] \end{split} $$ (9) 式中，ai,i为主对角元素；$ x $为方程的解向量；上标k、k+1为迭代步数。在计算过程中，只针对非零项进行计算。
4. 计算结果与分析
本研究实现了157组件反应堆堆芯的自动建模功能开发，并进行了计算。堆芯内燃料组件布置和编号如图1所示。每个燃料组件内共划分324个子通道，子通道编号如图2所示。表1给出了157组件反应堆堆芯几何参数和边界条件。
表 1 堆芯布置与关键参数Table 1. Geometirc Parameters of Test Section参数 数值 燃料组件数量 157 燃料棒直径/mm 9.5 燃料棒加热长度/m 3.65 燃料棒中心距/mm 12.6 堆芯功率/MW 3113 堆芯流量/(kg·s−1) 12494.1 冷却剂入口温度/℃ 288.3 堆芯出口压力/MPa 15.5 图6给出了157组件反应堆堆芯内燃料组件的径向功率因子分布情况，可以看到，中心区域组件功率高于边缘组件。图7给出了子通道径向功率因子分布情况，子通道径向功率分布为中心对称分布，中心区功率高于边缘区。
图8~图10分别给出了1号组件、50号组件及79号组件出口冷却剂的温度分布情况。由于堆芯内各组件的功率分布为中心区高于边缘区，1号组件的径向功率因子低于50号组件和79号组件，其冷却剂出口温度也低于另外2个组件。此外，从图8、图9中可以看到，1号组件和50号组件的出口冷却剂温度分布呈现出了上下对称而左右不对称的现象，结合组件径向功率因子分布，可以发现，位于这2个组件左右两侧的组件，其径向功率因子存在差异，与之相邻的右侧组件的功率要高于左侧组件，因此，组件出口冷却剂温度出现了右侧温度高于左侧的现象。对于79号组件，该组件位于堆芯中心，与之相邻的上下左右侧的组件径向功率因子是对称的，且上下侧组件功率因子高于左右侧组件。因此，该组件的出口冷却剂温度是对称的，且上下区域的冷却剂温度高于左右区域。
图 8 1号组件出口冷却剂温度分布Figure 8. Temperature Distribution of Coolant at the Outlet of Assembly 1
图 9 50号组件出口冷却剂温度分布Figure 9. Temperature Distribution of Coolant at the Outlet of Assembly 50
图 10 79号组件出口冷却剂温度分布Figure 10. Temperature Distribution of Coolant at the Outlet of Assembly 79图11给出了部分子通道（子通道170、16046、25442）的冷却剂焓值轴向分布结果的对比情况。这3个通道分别位于1号组件、50号组件及79号组件。可以看到，冷却剂焓值沿轴向高度不断增大。靠近出口侧，通道轴向功率为0，由于横向流动的影响，子通道温度出现了小幅度的下降。
图 11 子通道冷却剂焓值分布对比Figure 11. Comparison of Coolant Enthalpy Distribution of Sub-channel图12给出了部分子通道（子通道170、16046、25442）的冷却剂空泡份额轴向分布结果的对比情况，3个子通道分别位于1号组件、50号组件及79号组件的同一位置。可以看到，子通道170的最大空泡份额远小于子通道16046和子通道25442。虽然3个子通道的径向功率因子相同，但是1号组件的径向功率因子小于其他2个组件，子通道170的热量小于其他2个子通道。靠近出口侧，通道轴向功率为0，冷却剂空泡份额不断减小。
图 12 子通道冷却剂空泡份额分布对比Figure 12. Comparison of Coolant Void Fraction Distribution of Sub-channel图13给出了部分子通道（子通道170、16046、25442）的冷却剂DNBR分布结果的对比情况。由于各通道的轴向功率分布具有中间高、两侧低的特点，DNBR分布也呈现出了中心区域低于进出口区域的特点。组件径向功率因子越高，子通道DNBR越小。
5. 结 论
本研究开发了自动建模功能，解决了全堆芯pin-by-pin计算中大型稀疏矩阵的存储和求解问题，基于子通道分析程序LINDEN实现了压水堆堆芯pin-by-pin计算功能，开展了157组件反应堆堆芯pin-by-pin计算，并对计算结果进行了分析。主要结论如下：
（1）解决了LINDEN软件自动建模和大型系数矩阵的存储求解问题，实现了LINDEN软件全堆芯pin-by-pin计算功能。
（2）基于LINDEN软件，开展了全堆芯pin-by-pin计算，计算结果合理可信。
（3）由于组件间横向流动的影响，相邻组件的功率分布差异会影响组件的热工水力参数分布结果，全堆芯pin-by-pin计算能够获得更加精确的热工水力参数。
-
图 8 1号组件出口冷却剂温度分布
Figure 8. Temperature Distribution of Coolant at the Outlet of Assembly 1
图 9 50号组件出口冷却剂温度分布
Figure 9. Temperature Distribution of Coolant at the Outlet of Assembly 50
图 10 79号组件出口冷却剂温度分布
Figure 10. Temperature Distribution of Coolant at the Outlet of Assembly 79
图 11 子通道冷却剂焓值分布对比
Figure 11. Comparison of Coolant Enthalpy Distribution of Sub-channel
图 12 子通道冷却剂空泡份额分布对比
Figure 12. Comparison of Coolant Void Fraction Distribution of Sub-channel
表 1 堆芯布置与关键参数
Table 1. Geometirc Parameters of Test Section
参数 数值 燃料组件数量 157 燃料棒直径/mm 9.5 燃料棒加热长度/m 3.65 燃料棒中心距/mm 12.6 堆芯功率/MW 3113 堆芯流量/(kg·s−1) 12494.1 冷却剂入口温度/℃ 288.3 堆芯出口压力/MPa 15.5 
下载: 导出CSV
-
[1] 俞冀阳,贾宝山. 反应堆热工水力学[M]. 第二版. 北京:清华大学出版社,2011:200-210. [2] AVRAMOVA M N, SALKO R K. CTF theory manual: ORNL/TM-2016/430[R]. Oak Ridge: Oak Ridge National Lab. , 2016. [3] PÉRIN Y, VELKOV K. CTF/DYN3D multi-scale coupled simulation of a rod ejection transient on the NURESIM platform[J]. Nuclear Engineering and Technology, 2017, 49(6): 1339-1345. doi: 10.1016/j.net.2017.07.010 [4] MAGEDANZ J, PERIN Y, AVRAMOVA M, et al. High-fidelity multiphysics simulation of BWR assembly with coupled TORT-TD/CTF[C]//Proceedings of the Conference on Advances in Reactor Physics - Linking Research, Industry and Education. Knoxville: American Nuclear Society, Inc, 2012. [5] BENNETT A, AVRAMOVA M, IVANOV K. Coupled MCNP6/CTF code: development, testing, and application[J]. Annals of Nuclear Energy, 2016, 96: 1-11. doi: 10.1016/j.anucene.2016.05.008 [6] GUO J J, LIU S C, SHANG X T, et al. Coupled neutronics/thermal-hydraulics analysis of a full PWR core using RMC and CTF[J]. Annals of Nuclear Energy, 2017, 109: 327-336. doi: 10.1016/j.anucene.2017.05.041 [7] 刘宙宇,曹良志. 高保真物理-热工耦合计算方法研究及应用[J]. 原子核物理评论,2020,37(3):797-803. [8] 史敦福,李康,秦桂明,等. 蒙卡中子输运程序JMCT和子通道热工水力程序COBRA-EN耦合计算[J]. 强激光与粒子束,2017,29(3):036007. [9] 白宁,朱元兵,任志豪,等. 子通道分析程序LINDEN的开发与初步验证[J]. 原子能科学技术,2013,47(S1):299-301. [10] 郑成德. 数值计算方法[M]. 北京:清华大学出版社,2010:60-72. -
下载:
计量
- 文章访问数: 261
- HTML全文浏览量: 79
- PDF下载量: 36
- 被引次数: 0
