Study on JFNK Global Solution Method of of Full-core Thermal Sub-channel Model
-
摘要: 反应堆热工子通道模型详细考虑了轴向流动、横向交混，湍流交混等多种耦合因素，是堆芯热工水力分析的关键模型，但是这些因素为子通道数值模拟带来了困难和挑战。为了提高热工子通道模型的计算效率和收敛性，本文基于Jacobian-Free Newton-Krylov（JFNK）全局求解方法（以下简称JFNK方法），开发了全堆芯热工子通道模型的全局求解框架，并基于现有程序的模型和框架建立了基于物理预处理的残差系统，增强JFNK方法的收敛速率。结果表明，JFNK方法的计算效率是固定点迭代方法的5倍，且JFNK方法的效率优势随着收敛精度的提高会更加明显。因此，对于复杂热工子通道模型，JFNK方法有着不错的潜力和效率优势。Abstract: The thermal sub-channel model of reactor takes into account many coupling factors such as axial flow, lateral mixing and turbulent mixing in detail, and it is the key model of thermal hydraulic analysis of the core. However, these factors bring difficulties and challenges to the numerical simulation of sub-channel. In order to improve the computational efficiency and convergence of the thermal sub-channel model, this paper develops a global solution framework for the thermal sub-channel model of the whole core based on the Jacobian-Free Newton-Krylov (JFNK) global solution method (hereinafter referred to as JFNK method), and establishes a residual system based on physical pretreatment according to the existing code model and framework to enhance the convergence rate of JFNK method. Numerical solutions show that the speedup of JFNK methods is about five compared with the fixed-point iterative method used in the original sub-channel models and there are higher efficiency for JFNK methods as convergence precision improves, which indicates the potential and efficiency advantages of JFNK methods for complicated thermal sub-channel model.
-
Key words:
- Full-core thermal hydraulics /
- Sub-channel analysis /
- JFNK /
- Physical precondition
-
0. 引 言
热工子通道模型是包含轴向流动、横向交混、湍流交混等因素的三维复杂全堆芯流动传热耦合模型，是堆芯热工水力分析的关键模型。目前热工子通道分析程序普遍采用固定点迭代方法求解子通道模型，无论是采用“横向流”求解方案的COBRA-III[1]、COBRA-IIIC[2]、COBRA-V[3]等子通道程序，还是采用压力梯度求解方案的COBRA-EN[4]、FLICA-3F[5]以及中国核动力研究设计院自主开发的CORTH[6]等子通道程序普遍是各个守恒方程之间反复迭代并交换耦合数据，最终迭代至子通道模型收敛。因此模型越复杂，耦合数据越多，迭代次数也会越多，从而引起程序的计算效率和收敛性降低。为了解决计算效率和收敛性的问题，很多学者针对子通道程序的加速方法做了一些研究。例如COBRA-EN[4]程序在固定点迭代的框架下采用Newton-Raphson方法求解质量方程，相比于压力梯度求解方案程序的收敛性更强，计算更稳定。CTF[7]程序在求解压力修正量时，首先用Newton-Raphson方法将质量和能量方程局部线性化得到压力修正方程，然后根据压力矩阵规模的大小选择采用高斯消元法或者Krylov子空间方法求解得到压力修正量，进一步提高计算效率和收敛性。此外，FLICA-4[8]也采用广义极小残量法（GMRES）和稳定的双共轭梯度法 （Bi-CGSTAB）求解压力。然而这些研究普遍是在固定点迭代框架下进行的，依旧是各个方程分别求解的，各个方程之间仍采取反复迭代。
由于Jacobian-Free Newton-Krylov（JFNK）方法是一种全局求解方法，而且有着不错的收敛性和计算效率，近年来已开始尝试用于求解复杂热工子通道耦合问题。H.Esmaili[9]基于半隐式物理预处理建立了环形燃料组件热工水力子通道模型的残差，采用JFNK方法进行求解，结果与之前相比，总压降误差小于2%，其他关键热工参数吻合良好，但是对于棒状燃料组件并不适用。桑迪亚国家实验室基于CTF子通道程序开发残差版本CTF-R[10-11]，采用JFNK方法进行求解，与CTF程序相比，其有着更快的收敛速度和更强的收敛性，但是因为其只建立了液相流动模型的残差，所以只能模拟单相流动传热。目前JFNK方法求解全堆芯热工子通道模型的研究比较少，具有很强的研究价值。因此本文建立了基于全堆芯热工子通道模型的残差系统，实现了全堆芯热工子通道模型的JFNK全局求解。
1. 热工子通道耦合模型
本文采用了液相热不平衡和汽液两相速度不等的两相流模型。假设流体不可压缩且两相压力相等，流体的过冷物性参数和饱和物性参数分别是焓和压力的函数。基于以上假设，模型的守恒方程[12]如下：
（1）混合连续性方程：
$$ \frac{{\partial \rho }}{{\partial t}} = \frac{{\partial W}}{{\partial Z}} + \frac{1}{A}\mathop \sum \limits_{k \ne i} {s_{ik}} \cdot {w_{ik}} $$ (1) 式中，ρ为密度；t为时间；W为轴向质量流速；Z为轴向高度；A为子通道截面积；wik为横向质量流速；sik为横向交混间隙；下标i和k分别表示子通道i和k的变量。由于采用等质量的湍流交混模型[13]，混合连续性方程只包括质量流量的时间变化项、轴向和横向变化项。
（2）混合能量方程：
$$ \begin{split} \frac{\partial }{{\partial t}}&\left( {\rho {H_i}} \right) + \frac{1}{A}\mathop \sum \limits_{k \ne i} \left[ {\frac{{{s_{ik}}}}{{{L_{ik}}}}{k_{\text{t}}}\left( {{H_i} - {H_k}} \right)} \right] + \frac{1}{A}\frac{\partial }{{\partial Z}}\left( {W \cdot {H_i} \cdot A} \right) + \\ &\frac{1}{A}\mathop \sum \limits_{k \ne i} \left( {{s_{ik}} \cdot {w_{ik}} \cdot {H_{ik}}} \right) = \frac{1}{A} \cdot q \\[-20pt] \end{split}$$ (2) 式中，H为混合焓；Lik为相邻通道间质心距离；q为功率密度；kt为湍流热扩散系数。混合能量方程包括混合焓的时间变化项、湍流交混项、轴向和横向变化项和源项。
（3）液相能量方程：
$$\begin{split} \rho \left( {1 - \alpha } \right) \cdot& \frac{{\partial {H_{{\text{l}}i}}}}{{\partial t}} + \left( {1 - x} \right) \cdot \left[W\frac{{\partial {H_{{\text{l}}i}}}}{{\partial Z}} + \frac{1}{A}\mathop \sum \limits_{k \ne i} {s_{ik}} \cdot {w_{ik}} \cdot {H_{{\text{l}}ik}} +\right.\\ &\left.\frac{1}{A}\mathop \sum \limits_{k \ne i} \frac{{{s_{ik}}}}{{{L_{ik}}}}{k_{\rm{t}}}\left( {{H_{{\text{l}}i}} - {H_{{\text{l}}k}}} \right)\right] = \varGamma + \frac{1}{A}\chi \cdot q \\[-20pt] \end{split}$$ (3) 式中，Hl为液相焓；α、x分别为流体的空泡份额和含气率；
$ \text{ }\chi $ 为加热液相热流系数；$\varGamma $ 为气相冷凝成液相导致液相焓变化。液相能量包括液相焓的时间变化项、轴向和横向变化项、湍流交混项、冷凝项和源项。（4）轴向动量方程：
$$ \begin{aligned} - \frac{{\partial P}}{{\partial Z}} =& \frac{{\partial {W_i}}}{{\partial t}} + \rho g\cos \theta + \frac{1}{A} \cdot \frac{\partial }{{\partial Z}}(Av_i'{W_i}^2) + \\ &\frac{1}{A}\mathop \sum \limits_{k \ne i} {s_{ik}}v_{ik}'{W_{ik}} \cdot {w_{ik}} + {\tau _{\text{w}}} + {\tau _{\text{z}}} +\\ &\frac{1}{A}\mathop \sum \limits_{k \ne i} \frac{{{s_{ik}} \cdot {\mu _{\text{t}}}}}{{{L_{ik}}}}(v_i'{W_i} - v_k'{W_k}) \end{aligned} $$ (4) 式中，μt为湍流动力粘度；τw为壁面摩擦损失；τz为定位格架处局部阻力损失；
$v_i'$ 为有效动量迁移比容，$v_i' = \dfrac{{{x^2}{v_{\rm{g}}}}}{\alpha } + \dfrac{{{{\left( {1 - x} \right)}^2}{v_{\text{l}}}}}{{1 - \alpha }}$ ，vg、vl分别为气相和液相比体积；g为重力加速度；θ为子通道与垂直方向夹角。轴向动量方程包括压力项、轴向动量的时间变化项、重力项、轴向和横向变化项、壁面摩擦项、格架阻力项和湍流交混项。（5）横向动量方程：
$$ - \frac{{\partial P}}{{\partial {y_{ik}}}} = \frac{{\partial {w_{ik}}}}{{\partial t}} + \tau _{\text{w}}' + \frac{\partial }{{\partial Z}}\left( {v_{ik}' \cdot {W_{ik}} \cdot {w_{ik}}} \right) + \frac{\partial }{{\partial {y_{ik}}}}\left( {v_{ik}' \cdot w_{ik}^2} \right) $$ (5) 式中，
$\tau _{\text{w}}'$ 为子通道间横向摩擦损失；yik为横向距离。横向动量方程包括压力项、横向动量的时间变化项、横向摩擦力项、轴向和横向变化项。2. JFNK全局求解方法
2.1 全堆芯热工子通道残差系统
如图1所示，本文基于交错网格的有限体积法对稳态守恒方程（1）~方程（5）进行离散，空间离散采用一阶迎风格式。子通道耦合界面上的变量通过相连子通道间的流通截面积权重F来计算，即
${F_{ik}} = {m_1} \cdot {F_i} + \left( {1 - {m_1}} \right) \cdot {F_k}$ ，其中，${m_1} = \dfrac{{{A_k}}}{{{A_i} + {A_k}}}$ 。横向交混量N由流向标识因子确定，即${N_{ik}} = {m_2} \cdot {N_i} + (1 - {m_2}) \cdot {N_k}$ ，${m_2} = \left\{ {\begin{array}{*{20}{c}} {1 \;\;{\text{ }}{w_{ik}} \geqslant 0} \\ {0 \;\;{\text{ }}{w_{ik}} < 0} \end{array}} \right.$ 。对于通道i、轴向控制体j，具体离散模型见开放科学（OSID）标识码中文档。对于上述控制方程组成的耦合系统，传统固定点迭代方法的求解流程如图2所示。图2中显示：先求解混合焓和液相焓方程更新流体的温度场，再联立求解轴向动量方程和横向动量方程得到流体的速度场和压力场。可清晰地看到该求解方法中存在多层迭代，而随着问题复杂度的增加，每层迭代过程的收敛速率会降低，再加上多层迭代之间变量反复更新，最终可能导致整个热工耦合问题的收敛速率下降，严重影响程序的计算效率。为了提高计算效率，本文基于JFNK方法，实现全堆芯热工子通道模型的全局求解。本文选取混合焓、轴向质量流速、横向质量流速、轴向压降、横向压差∆Pik作为JFNK全局求解的残差变量，基于离散方程建立了全堆芯热工子通道模型的残差系统，如下所示：
$$ {\boldsymbol{F}}\left( {\boldsymbol{x}} \right) = {\left[ {{F_1},{F_2},{F_3},{F_4},{F_5}} \right]^{\rm{T}}} $$ (6) 式中，
${\boldsymbol{x}} = {\left[ {H,W,{w_{ik}},\Delta P,\Delta {P_{ik}}} \right]^{\rm{T}}}$ ；F1、F2、F3、F4和F5分别为离散方程中等号左端项。2.2 基于原有程序框架的物理预处理
在JFNK方法[14]中，首先采用Newton法将非线性方程组（6）局部线性化，如式（7）所示。
$$ {\boldsymbol{J}}\left({\boldsymbol{x}}\right)\cdot \delta {\boldsymbol{x}} = -{\boldsymbol{F}}\left({\boldsymbol{x}}\right) $$ (7) $$ {\text{ }}{{\boldsymbol{x}}^{k + 1}} = {{\boldsymbol{x}}^k} + \delta {\boldsymbol{x}} $$ 式中，
${\boldsymbol{J}}\left( {\boldsymbol{x}} \right)$ 是Jacobian矩阵；上标$k$ 为Newton迭代次数。然后采用Krylov子空间方法求解线性化的线性方程组，由于Krylov子空间方法仅需要Jacobian矩阵与向量乘积，而并没有明确需要Jacobian矩阵的具体形式，因此采用一阶差分代替微分从而构造Jacobian矩阵与向量乘积，即Jacobian-free技术，如式（8）所示。$$ {{\boldsymbol{P}}^{ - 1}} \cdot {\boldsymbol{J}}\left( {\boldsymbol{x}} \right) \cdot {\boldsymbol{v}} \approx \frac{{{{\boldsymbol{P}}^{ - 1}}{\boldsymbol{F}}\left( {{\boldsymbol{x}} + \varepsilon \cdot {\boldsymbol{v}}} \right) - {{\boldsymbol{P}}^{ - 1}}F\left( {\boldsymbol{x}} \right)}}{\varepsilon } $$ (8) 式中，
${\boldsymbol{v}}$ 为Krylov子空间的基向量；ε为微小扰动量。JFNK方法的计算量主要来自Krylov方法求解线性方程组，为了加快Krylov方法的收敛速率，本文基于原有固定点迭代框架进行了物理预处理。固定点迭代思路是将矩阵A可分解为2个矩阵P和Q，即A=P−Q，经过化简可以得到迭代矩阵：$$ \begin{aligned} &{\boldsymbol{F}}\left( {\boldsymbol{x}} \right) = {\boldsymbol{Ax}} - {\boldsymbol{b}} = \left( {{\boldsymbol{P}} - {\boldsymbol{Q}}} \right){\boldsymbol{x}} - {\boldsymbol{b}} = 0 \\ & {{\boldsymbol{x}}^{{\text{out}}}} = {{\boldsymbol{P}}^{ - 1}}\left( {{\boldsymbol{Qx}} + {\boldsymbol{b}}} \right) \\[-14pt] \end{aligned} $$ (9) 将矩阵P作为预处理矩阵，预处理后的残差
${{\boldsymbol{P}}^{ - 1}}{\boldsymbol{F}}\left( {\boldsymbol{x}} \right)$ 、${{\boldsymbol{P}}^{ - 1}}{\boldsymbol{F}}\left( {{\boldsymbol{x}} + \varepsilon \cdot {\boldsymbol{v}}} \right)$ 可以得到简化：$$ \begin{aligned} & {{\boldsymbol{P}}^{ - 1}}{\boldsymbol{F}}\left( {\boldsymbol{x}} \right) = {{\boldsymbol{P}}^{ - 1}}\left[ {\left( {{\boldsymbol{P}} - {\boldsymbol{Q}}} \right){\boldsymbol{x}} - {\boldsymbol{b}}} \right] \\ & ={\boldsymbol{x}} - {{\boldsymbol{P}}^{ - 1}}\left( {{\boldsymbol{Qx}} + {\boldsymbol{b}}} \right) = {\boldsymbol{x}} - {{\boldsymbol{x}}^{{\text{out}}}} \end{aligned} $$ (10) $$ {{\boldsymbol{P}}^{ - 1}}{\boldsymbol{F}}\left( {{\boldsymbol{x}} + \varepsilon \cdot {\boldsymbol{v}}} \right) = {\boldsymbol{x}} + \varepsilon \cdot {\boldsymbol{v}} - {{\boldsymbol{x}}^{{\text{out}},\varepsilon }} $$ (11) 式中，
${{\boldsymbol{x}}^{{\text{out}}}}$ 、${{\boldsymbol{x}}^{{\text{out}},\varepsilon }}$ 分别为${\boldsymbol{x}}$ 、${\boldsymbol{x}} + \varepsilon \cdot {\boldsymbol{v}}$ 作为输入时一次固定点迭代的输出结果。同时式（8）可以化简为：$$ {{\boldsymbol{P}}^{ - 1}} \cdot {\boldsymbol{J}}({\boldsymbol{x}}) \cdot {\boldsymbol{v}} \approx \frac{{{{\boldsymbol{x}}^{{\text{out}}}} - {{\boldsymbol{x}}^{{\text{out}},\varepsilon }}}}{\varepsilon } + {\boldsymbol{v}} $$ (12) 此外模型中的压力场与横向质量流速场在数量级上存在较大差异，为了避免这种差异，本文基于模型中每一个物理场变量的范数对向量
$ {\boldsymbol{x}} $ 进行Scaling处理[15]。经过以上处理构建了基于原有模型框架的物理预处理残差系统。2.3 求解流程和迭代策略
JFNK求解流程如图3所示[16]。JFNK求解的最终收敛准则如下所示：
$$ \frac{{\Vert {{\boldsymbol{F}}}^{\prime }\left({{\boldsymbol{x}}}^{k}\right)\Vert }_{2}}{{\Vert {{\boldsymbol{F}}}^{\prime }\left({{\boldsymbol{x}}}^{0}\right)\Vert }_{2}}\leqslant {\rm{rtol}} $$ (13) 式中，|| ||2为L2范数；
${\boldsymbol{F'}}\left( {\boldsymbol{x}} \right)$ 为预处理后的残差系统，即$ {\boldsymbol{F'}}\left( {\boldsymbol{x}} \right){\boldsymbol{ = }}{{\boldsymbol{P}}^{ - 1}}{\boldsymbol{F}}\left( {\boldsymbol{x}} \right) $ ；${\left\| {{\boldsymbol{F'}}\left( {{{\boldsymbol{x}}^0}} \right)} \right\|_2}$ 为初始范数；rtol为残差收敛准则。为了保证Krylov迭代的效率和收敛性，每个Newton迭代步中Krylov迭代的收敛准则可基于Newton迭代的残差范数收敛情况自动获得，如式（13）所示，其中，ηk为强制项系数[17]，可由式（14）计算。
$$ {\left\| {{\boldsymbol{J'}}\left( {{{\boldsymbol{x}}^k}} \right) \cdot \delta {{\boldsymbol{x}}^k} + {\boldsymbol{F'}}\left( {{{\boldsymbol{x}}^k}} \right)} \right\|_2} \leqslant {\eta _k}{\left\| {{\boldsymbol{F'}}\left( {{{\boldsymbol{x}}^k}} \right)} \right\|_2} $$ (14) $$ {\eta _k} = \frac{{\left| {{{\left\| {{\boldsymbol{F'}}\left( {{{\boldsymbol{x}}^k}} \right)} \right\|}_2} - {{\left\| {{\boldsymbol{J'}}\left( {{{\boldsymbol{x}}^{k - 1}}} \right) \cdot \delta {{\boldsymbol{x}}^{k - 1}} + {\boldsymbol{F'}}\left( {{{\boldsymbol{x}}^{k - 1}}} \right)} \right\|}_2}} \right|}}{{{{\left\| {{\boldsymbol{F'}}\left( {{{\boldsymbol{x}}^{k - 1}}} \right)} \right\|}_2}}} $$ (15) 式（8）中微扰量ε可由式（15）计算[18]：
$$ \varepsilon =\frac{{b}_{\text{mach}}\sqrt{\left(1+||{{\boldsymbol{x}}}^{k}||_{2}\right)}}{{||{\boldsymbol{v}}}||_{2}} $$ (16) 式中，bmach为机器精度的平方根；J′(x)为预处理后的Jacobian矩阵。
3. 结果对比与分析
为了验证JFNK方法的数值特性，本文用百万千瓦级压水堆全堆芯算例开展稳态热工子通道模型求解方法对比研究。全堆芯组件编号和径向功率分布如图4所示。轴向功率采用对称余弦功率分布，其他参数如表1所示。
图 4 全堆芯组件编号及径向功率分布Figure 4. Assembl Number and Radial Power Distribution of Full Core参数 量值 燃料组件数目 157 燃料组件边长/cm 21.504 燃料棒数目（每个燃料组件） 264 导向管和控制棒数目（每个燃料组件） 25 燃料棒外直径/cm 0.95 导向管外径/cm 1.245 棒中心间距/cm 1.26 参考热流密度/(W·cm−2) 55.6534 堆芯入口流量/(m3·h−1) 64066.2 冷却剂入口温度/℃ 292.4 堆芯出口压力/105Pa 155 轴向功率因子峰值 1.55 轴向加热区起始点高度/cm 15.25 轴向加热区长度/cm 366.0 本文在径向上将每个组件划分为2×2网格（图5），在轴向方向上划分为27个网格，分别用固定点迭代方法和JFNK方法求解，计算结果如图6~图9所示。图6是堆芯出口的混合焓和轴向质量流速分布图，根据堆芯径向功率分布特性以及出口焓和质量流速的分布特性，本文选取了2个组件进行结果对比分析，编号分别为79和52。由图可知：JFNK方法与固定点迭代方法的计算结果吻合得很好。组件79轴向质量流速先增大后减小，而组件52是一直增大的；混合焓2个组件都是增加的，只是组件79增加得更快，混合焓也更高；压力变化两者基本一致，从出口到入口都呈现下降趋势，在有定位格架位置压力变化更加剧烈。
图 7 组件79和52对应通道轴向质量流速比较Figure 7. Comparison of Axial Mass Fluxes in Assemblies 79 and 52
图 8 组件79和52对应通道轴向混合焓比较Figure 8. Comparison of Axial Mixing Enthalpies in Assemblies 79 and 52为了对比分析不同耦合方法的计算效率，图10给出了固定点迭代方法和JFNK方法在范数精度1×10−6下的收敛曲线，表2给出了这2种耦合方法在不同范数精度下的计算时间。可以看出，JFNK方法与固定点迭代方法相比，其需要更少的非线性迭代次数。固定点迭代方法收敛曲线会缓慢下降，并且在收敛到一定精度后会产生振荡，而JFNK方法可快速光滑收敛。在范数精度10−6下，JFNK方法的计算效率是固定点迭代方法的5倍左右，并且JFNK方法的效率优势会随着收敛精度的提高而更明显。
表 2 不同收敛精度对应的2种方法的收敛性比较Table 2. Comparison of Convergence of Two Methods with Different Convergence Accuracies收敛精度 固定点迭代方法 JFNK全局求解方法 加速比 时间/s 迭代次数 时间/s Newton迭代次数 1×10−3 45.0 95 16.97 4 2.65 1×10−4 89.8 213 24.92 4 3.60 1×10−5 153.4 336 34.91 5 4.39 1×10−6 194.5 454 39.10 5 4.97 此外为了提高固定点迭代方法的收敛速率和改善收敛性，通常需要匹配相应的松弛因子。在本文模型中，相邻通道的横向压降使用了松弛因子φ，即
$\Delta P_{ik}^m = \varphi \cdot \Delta P_{ik}^{m - 1} + \left( {1 - \varphi } \right) \cdot \Delta P_{ik}^m$ ，其中m为当前迭代步。在固定点迭代方法中，不同松弛因子会给迭代收敛速率带来较大影响。图11给出了不同松弛因子下固定点迭代收敛曲线，表3给出了不同松弛因子下固定点迭代和JFNK收敛特性和计算效率。由图表可见：松弛因子取0.6~0.8，固定点迭代能够收敛，当持续增大松弛因子，收敛性会越来越差；而当松弛因子取0.5时，迭代出现了发散现象。因此固定点迭代中松弛因子选取不合适会使收敛变慢甚至发散。但是对于JFNK方法，改变松弛因子对JFNK求解的计算时间没有太大的影响，此外松弛因子为0.5时固定点迭代方法会发散，而在JFNK方法中却可以快速收敛。由此可见，JFNK方法有更好的鲁棒性和高效性，一定程度上避免了固定点迭代中针对不同问题要选取有效松弛因子的困难。
图 11 不同松弛因子下固定点迭代收敛曲线Figure 11. Convergence Curves of Fixed Point Iteration with Different Relaxation Factors表 3 不同松弛因子下的2种方法的迭代次数和计算时间Table 3. Iterative Number and Computation Time of Two Methods with Different Relaxation Factorsφ 固定点迭代方法 JFNK全局求解方法 固定点迭代步 时间/s Newton迭代步 时间/s 0.50 无法计算 5 35.9 0.60 454 194.5 5 39.7 0.70 994 422.4 5 39.2 0.80 2376 1172.8 11 54.3 4. 结 论
本文建立了基于全堆芯热工子通道模型的残差系统，实现了三维全堆芯热工子通道模型的JFNK全局求解。相比于传统的固定点迭代方法，JFNK方法拥有更好的计算效率和收敛优势，且随着收敛精度的提高或计算规模的扩大，收敛优势会更加明显。对于本文模型，传统的固定点迭代收敛过程会出现数值振荡，且松弛因子的选取将严重影响收敛速率，甚至可能发散，而JFNK方法由于各个方程之间同步更新变量，其收敛快速且收敛曲线光滑无振荡，具有更好的鲁棒性和高效性。本文建立了组件级别的三维全堆芯热工子通道模型，计划未来将JFNK方法扩展用于三维全堆芯热工子通道pin-by-pin级别模型的计算研究。
-
图 4 全堆芯组件编号及径向功率分布
Figure 4. Assembl Number and Radial Power Distribution of Full Core
图 7 组件79和52对应通道轴向质量流速比较
Figure 7. Comparison of Axial Mass Fluxes in Assemblies 79 and 52
图 8 组件79和52对应通道轴向混合焓比较
Figure 8. Comparison of Axial Mixing Enthalpies in Assemblies 79 and 52
图 11 不同松弛因子下固定点迭代收敛曲线
Figure 11. Convergence Curves of Fixed Point Iteration with Different Relaxation Factors
参数 量值 燃料组件数目 157 燃料组件边长/cm 21.504 燃料棒数目（每个燃料组件） 264 导向管和控制棒数目（每个燃料组件） 25 燃料棒外直径/cm 0.95 导向管外径/cm 1.245 棒中心间距/cm 1.26 参考热流密度/(W·cm−2) 55.6534 堆芯入口流量/(m3·h−1) 64066.2 冷却剂入口温度/℃ 292.4 堆芯出口压力/105Pa 155 轴向功率因子峰值 1.55 轴向加热区起始点高度/cm 15.25 轴向加热区长度/cm 366.0 
下载: 导出CSV
表 2 不同收敛精度对应的2种方法的收敛性比较
Table 2. Comparison of Convergence of Two Methods with Different Convergence Accuracies
收敛精度 固定点迭代方法 JFNK全局求解方法 加速比 时间/s 迭代次数 时间/s Newton迭代次数 1×10−3 45.0 95 16.97 4 2.65 1×10−4 89.8 213 24.92 4 3.60 1×10−5 153.4 336 34.91 5 4.39 1×10−6 194.5 454 39.10 5 4.97
下载: 导出CSV
表 3 不同松弛因子下的2种方法的迭代次数和计算时间
Table 3. Iterative Number and Computation Time of Two Methods with Different Relaxation Factors
φ 固定点迭代方法 JFNK全局求解方法 固定点迭代步 时间/s Newton迭代步 时间/s 0.50 无法计算 5 35.9 0.60 454 194.5 5 39.7 0.70 994 422.4 5 39.2 0.80 2376 1172.8 11 54.3
下载: 导出CSV
-
[1] ROWE D S. COBRA-III: a digital computer program for steady state and transient thermal--hydraulic analysis of rod bundle nuclear fuel elements: BNWL-B-82[R]. Richland: Battelle Pacific Northwest Labs. , 1971. [2] ROWE D S. COBRA IIIC: digital computer program for steady state and transient thermal-hydraulic analysis of rod bundle nuclear fuel elements: BNWL-1695[R]. Richland: Battelle Pacific Northwest Labs. , 1973. [3] STEWART C W, WHEELER C L, CENA R J, et al. COBRA-IV: the model and the method: BNWL-2214[R]. Richland: Pacific Northwest Laboratories, 1977. [4] BASILE D, BEHGI M, CHIERICI R, et al. COBRA-EN: Code System for Thermal-Hydraulic Transient Analysis of Light Water Reactor Fuel Assemblies and Cores: RSICC Code Package PSR-507[R]. U.S.: Oak-Ridge, 2001. [5] PLAS R. FLICA III M-reactors or test loops thermohydraulic computer code: CEA-N--2418[R]. Gif-sur-Yvette: CEA Centre d'Etudes Nucleaires de Saclay, 1984. [6] DENG J, LU Q, LIU Y, et al. Review of sub-channel code development for pressurized water reactor and introduction of CORTH-V2.0 sub-channel code[J]. Progress in Nuclear Energy, 2020, 125: 103373. doi: 10.1016/j.pnucene.2020.103373 [7] HU J W, SALKO JR R K, WYSOCKI A J. CTF 4.0 theory manual[R]. Oak Ridge: Oak Ridge National Lab. , 2019. [8] 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 Engineering and Design, 2000, 200(1-2): 139-155. doi: 10.1016/S0029-5493(99)00332-5 [9] ESMAILI H, KAZEMINEJAD H, KHALAFI H, et al. Subchannel analysis of annular fuel assembly using the preconditioned Jacobian-free Newton Krylov methods[J]. Annals of Nuclear Energy, 2020, 146: 107616. doi: 10.1016/j.anucene.2020.107616 [10] PORTER N W. Development of a novel residual formulation of CTF and application of parameter estimation techniques[D]. Raleigh: North Carolina State University, 2018. [11] PORTER N W, MOUSSEAU V A, AVRAMOVA M N. CTF-R: a novel residual-based thermal hydraulic solver[J]. Nuclear Engineering and Design, 2019, 348: 37-45. doi: 10.1016/j.nucengdes.2019.04.006 [12] PLAS R. FLICA III M-reactors or test loops thermohydraulic computer code[R]. Gif-sur-Yvette: CEA Centre d'Etudes Nucleaires de Saclay, 1984. [13] MOORTHI A, SHARMA A K, VELUSAMY K. A review of sub-channel thermal hydraulic codes for nuclear reactor core and future directions[J]. Nuclear Engineering and Design, 2018, 332: 329-344. doi: 10.1016/j.nucengdes.2018.03.012 [14] KNOLL D A, KEYES D E. Jacobian-free Newton–Krylov methods: a survey of approaches and applications[J]. Journal of Computational Physics, 2004, 193(2): 357-397. doi: 10.1016/j.jcp.2003.08.010 [15] CHAN T F, VAN DER VORST H A. Approximate and incomplete factorizations[M]//KEYES D E, SAMEH A, VENKATAKRISHNAN V. Parallel Numerical Algorithms. Dordrecht: Springer, 1997: 167-202. [16] ZHOU X F, ZHONG C M, ZHANG Y Y. Jacobian-free Newton Krylov two-node coarse mesh finite difference based on nodal expansion method for multiphysics coupled models[J]. Annals of Nuclear Energy, 2022, 168: 108915. doi: 10.1016/j.anucene.2021.108915 [17] EISENSTAT S C, WALKER H F. Choosing the forcing terms in an inexact Newton method[J]. SIAM Journal on Scientific Computing, 1996, 17(1): 16-32. doi: 10.1137/0917003 [18] PERNICE M, WALKER H F. NITSOL: a Newton iterative solver for nonlinear systems[J]. SIAM Journal on Scientific Computing, 1998, 19(1): 302-318. doi: 10.1137/S1064827596303843 [19] 广东核电培训中心. 900MW压水堆核电站系统与设备[M]. 北京: 原子能出版社, 2005: 611. [20] 梁志滔. 压水堆核电站堆芯子通道分析[D]. 广州: 华南理工大学, 2011. -
下载:
计量
- 文章访问数: 1810
- HTML全文浏览量: 278
- PDF下载量: 392
- 被引次数: 0
