作为第四代先进核能系统的一种，铅铋堆主回路采用池式设计^[1]。在文献中提及，池式反应堆设计会导致一系列与液态金属冷却剂有关的复杂热工水力学和安全问题，如腔室以及冷热池的热分层现象。现已开发出多个程序：David Pialla等人针对池式钠冷快堆，开展了系统分析程序CATHARE与CFD程序TRIO_U的耦合以及系统分析程序ATHLET与CFD程序OPENFOAM的耦合，以进行精确的瞬态分析^[2]；Zhi-Xing Gu等人基于自主开发的二维CFD程序，并耦合自主开发的二维燃料棒传热程序及点堆动力学模型，开发了二维热工水力计算程序MPC_LBE^[3]。上述为典型的多维程序，虽然能够对池式反应堆的热分层现象进行精确的模拟，但是其所需的计算资源较为庞大，特别是对反应堆进行长时间模拟时，这个问题显得尤为突出。

与多维程序相比，系统分析程序能够大幅减小计算资源消耗，且对堆芯的模拟也较为精确。为此，国际上已经开发或二次开发了几个大型商用系统安全分析程序，如SAS4A/SASSYS^[4]、FAST^[5]以及CATHARE^[6]。其中，SAS4A与FAST分别是美国与瑞典开发出的针对快堆的专用系统分析程序，可对铅铋堆、钠冷快堆等进行稳态及瞬态模拟；CATHARE是法国开发的反应堆系统分析程序，其早期主要是针对水堆的系统安全分析而开发的大型商用程序，后来为了适用于铅基反应堆，进行了二次开发，推出了衍生版本。但是这些商用程序的源代码不易获得，相较于自主开发的小型微机系统分析程序，在针对某些特殊反应堆时不能进行灵活的修改。

在本研究中，将基于二维轴对称柱坐标的燃料棒传热程序以及零维点堆动力学模型集成到自主开发的一维单相CFD程序中，实现“中子-热工”耦合，得到了一款用于铅铋冷却反应堆瞬态安全分析的多物理场耦合程序。并通过加速器驱动的次临界系统（XADS反应堆）进行稳态验证以及XADS的失束（BT）事故进行瞬态验证。

1 程序模型

1.1 一维CFD

1.1.1 控制方程及其离散

由于堆芯的多个通道以及反应堆的多个环路在本质上是一个多维热工水力问题，因而本文采用集**数法，使用一维CFD程序进行模拟。对铅铋堆中冷却剂流动及传热的数值计算，由一维单相不可压缩牛顿型流体的连续性方程、动量方程及能量方程模拟给出，即

$ \dfrac{{\partial \rho }}{{\partial t}} + \dfrac{{\partial \rho v}}{{\partial {\textit{z}}}} = 0 $ (1)

$ \dfrac{{\partial \rho v}}{{\partial t}} + \dfrac{{\partial \rho vv}}{{\partial t}} = - \dfrac{{\partial p}}{{\partial {\textit{z}}}} + \dfrac{\partial }{{\partial {\textit{z}}}}\left( {\mu \dfrac{{\partial v}}{{\partial {\textit{z}}}}} \right) + \rho g{\text{cos}}\theta - \dfrac{f}{D}\left( {\dfrac{1}{2}\rho {v^2}} \right) $ (2)

$ \dfrac{{\partial \rho {c_p}T}}{{\partial \tau }} + \dfrac{{\partial \rho {c_p}T_{\rm{v}}}}{{\partial {\textit{z}}}} = \dfrac{\partial }{{\partial {\textit{z}}}}\left( {{\lambda _{{\rm{LBE}}}}\dfrac{{\partial T}}{{\partial {\textit{z}}}}} \right) + {q_{{{\rm{v}}}}} $ (3)

式中 $ ：\rho $为密度； $ v $为速度； $ \mu $为动力粘度； $ p $为压力； $ g $为重力加速度； $ \theta $为流体流动方向与重力方向之间的夹角； $ f $为阻力系数，其值由经验公式给出^[7-8]； $ D $为流道水力直径； $ T $为温度； $ {c}_{p} $为比热容； $ \lambda $为导热系数； $ {q}_{{{\rm{v}}}} $为体积热流密度。

对于能量方程（3）中的源项 $ {q}_{\rm{v}} $，其在堆芯流道中，为燃料棒包壳与冷却剂之间的换热体积功率；在换热器流道中，为一二回路之间的换热体积功率；在腔室及冷、热池及管道中，忽略管壁与外界的热交换，因而其值为零。

在求解上述方程时，我们通常采用交错网格对计算空间进行离散，如图1所示：将温度、压力等标量信息储存于主网格节点上，将速度信息储存于主网格节点的界面上，故在控制体 $ (j-{1}/{2},j+{1}/{2}) $上对连续性方程及能量方程离散；在控制体 $ (j-1,j) $上对动量方程进行离散。

图 1. 交错网格

Fig. 1. Staggered grid

下载图片 查看所有图片

首先对连续性方程进行离散。我们认为液态铅铋合金为不可压缩流体，其密度仅由温度所决定；其导热系数及动力粘度等热物性也仅由温度所决定；其比热容为一恒定值。本文所使用的液态铅铋合金的各项物性参数均取自于基准报告^[9]。

在连续性方程中，使用二阶迎风格式对对流项进行离散，得到离散后的连续性方程如式（4）所示，其中 $ A $为流道截面积。

$ \rho _{j + 1/2}^nv_{j + 1/2}^{n + 1}{A_{j + 1/2}} - \rho _{j - 1/2}^nv_{j - 1/2}^{n + 1}{A_{j - 1/2}} = 0 $ (4)

在连续性方程离散的基础上，对能量方程进行离散。离散后的能量方程如式（5）所示，式中 $ \Delta z $表示控制体长度。其扩散项使用差分格式离散，对流项使用二阶迎风格式离散，并将离散后的扩散项和对流项代入式（5），可整理得到一个线性方程组，如式（6）所示。

$ \begin{split} \Delta {{\textit{z}}_j}{A_j}\rho _j^n{c_p}_j^n\dfrac{{(T_j^{n + 1} - T_j^n)}}{{\Delta t}} - &\Delta {{\textit{z}}_j}{A_j}{(q_v^{})_j} + \left[ {A_{j + 1/2}}\rho _{j + 1/2}^n({c_p})_{j + 1/2}^nT_{j + 1/2}^{n + 1}v_{j + 1/2}^n - {A_{j - 1/2}}\rho _{j - 1/2}^n({c_p})_{j - 1/2}^nT_{j - 1/2}^{n + 1}v_{j - 1/2}^n \right]= \\ & \left[ {{A_{j + 1/2}}\lambda _j^n\dfrac{{T_{j + 1/2}^{n + 1} - T_j^{n + 1}}}{{\Delta {{\textit{z}}_j}/2}} - {A_{j - 1/2}}\lambda _j^n\dfrac{{T_j^{n + 1} - T_{j - 1/2}^{n + 1}}}{{\Delta {{\textit{z}}_j}/2}}} \right] \end{split} $ (5)

$ {a_j}T_j^{n + 1} = {a_{j + 2}}T_{j + 2}^{n + 1} + {a_{j + 1}}T_{j + 1}^{n + 1} + {a_{j - 1}}T_{j - 1}^{n + 1} + T_{j - 2}^{n + 1} + {b_j} $ (6)

对动量方程进行离散时，同样的，扩散项使用差分格式离散，对流项使用二阶迎风格式离散，并对源项中的速度平方做线性化处理，得到离散后的动量方程，如式（7）所示。最后，经整理可得到动量方程的代数形式，其为一个线性方程组，如式（8）所示。

$ \begin{split} & \dfrac{{\Delta {{\textit{z}}_{j - 1/2}}\rho _{j - 1/2}^{n + 1}{A_{j - 1/2}}}}{{\Delta t}}(v_{j - 1/2}^{n + 1} - v_{j - 1/2}^n) + [\rho _j^{n + 1}v_j^n{A_j}v_j^{n + 1} - \rho _{j - 1}^nv_{j - 1}^n{A_{j - 1}}v_{j - 1}^{n + 1}] = \left( {{A_j}\mu _j^{n + 1}\dfrac{{v_{j + 1/2}^{n + 1} - v_{j - 1/2}^{n + 1}}}{{\Delta {{\textit{z}}_j}}} - {A_j}\mu _{j - 1}^{n + 1}\dfrac{{v_{j - 1/2}^{n + 1} - v_{j - 3/2}^{n + 1}}}{{\Delta {{\textit{z}}_{j - 1}}}}} \right) -\\ & (p_j^{n + 1} - p_{j - 1}^{n + 1}) + {A_{j - 1/2}}\rho _{j - 1/2}^{n + 1}\Delta {{\textit{z}}_{j - 1/2}}g\cos \theta + \dfrac{f}{{2D}}{A_{j - 1/2}}\rho _{j - 1/2}^{n + 1}\Delta {{\textit{z}}_{j - 1/2}}{(v_{j - 1/2}^n)^2} -\dfrac{f}{D} {A_{j - 1/2}}\rho _{j - 1/2}^{n + 1}\Delta {{\textit{z}}_{j - 1/2}}v_{j - 1/2}^nv_{j - 1/2}^{n + 1} \end{split} $ (7)

$ a_{j - 1/2}^{\rm{v}}v_{j - {1}/{2}}^{n + 1} = a_{j - 5/2}^{\rm{v}}v_{j - {5}/{2}}^{n + 1} + a_{j - 3/2}^{\rm{v}}v_{j - {3}/{2}}^{n + 1} + a_{j + 1/2}^{\rm{v}}v_{j + {1}/{2}}^{n + 1} + a_{j + 3/2}^{\rm{v}}v_{j + {3}/{2}}^{n + 1} + c_{j - 1/2}^{\rm{v}} $ (8)

1.1.2 程序算法及流程

由于在SIMPLE算法中，我们既假设了压力场又假定了速度场，这两者之间不一定协调，所以会对迭代收敛的速度产生影响。因此，在这里采用SIMPLER算法^[10]，其主要步骤如图2所示：首先通过上一时刻的速度场作为本次迭代的初始速度场，求解出假拟速度场 $ \hat {v} $；通过这个假拟的速度场求解出待修正的压力场 ${p}^{*}$；接下来通过待修正的压力场求解动量方程，得到待修正的速度场 $ {v}^{*} $；通过待修正的速度场求解压力修正方程，获得修正的压力场 $p$；通过修正的压力场对速度场进行修正，得到修正后的速度场 $ v $；最后通过连续性方程检验修正后的速度场是否收敛。

图 2. 算法流程

Fig. 2. Algorithm process

下载图片 查看所有图片

1.1.3 修正方程的构建

在SIMPLER算法中需要构建压力修正方程以对压力场进行修正，进而对速度场进行修正。我们用 $ {p}^{*}、{v}^{*} $表示待修正的压力及速度， $ {p}'、{v}' $表示压力场的修正量及速度场的修正量，则可将修正后的压力及速度表示为

$ \left\{ {\begin{array}{*{20}{c}} {v = {v^{\text{*}}} + v'} \\ {p = {p^{\text{*}}} + p'} \end{array}} \right. $ (9)

由于 $ {v}^{*} $是求解动量方程得来的，因而其满足动量方程，并且假设修正后的 $ v $也满足动量方程，将 $ v{、v}^{*} $分别代入动量方程，并使二者相减，忽略相邻控制体的影响，可得到修正后的速度，如式（10）所示，其中 $ {a}_{j-1/2}^{v} $表示动量方程中以 $ {v}_{j-{1}/{2}}^{n+1} $为研究对象时 $ {v}_{j-{1}/{2}}^{n+1} $项的系数。

$ {v_{j - 1/2}} = v_{j - 1/2}^{\text{*}} + \dfrac{{{A_{j - 1/2}}}}{{a_{j - 1/2}^{\rm{v}}}}\left( {p_{j - 1}' - p_j'} \right) $ (10)

将修正后的速度代入到离散后的连续性方程中，可以得到压力修正方程

$ b_j^{\rm{p}}p_j' = b_{j - 1}^{\rm{p}}p_{j - 1}' + b_{j + 1}^{\rm{p}}p_{j + 1}' + d_j^{\rm{p}} $ (11)

$ \left\{ \begin{gathered} b_{j - 1}^{\rm{p}} = A_{j - 1/2}^2\dfrac{{\rho _{j - 1/2}^n}}{{a_{j - 1/2}^{\rm{v}}}} \\ b_{j + 1}^{\rm{p}} = A_{j + 1/2}^2\dfrac{{\rho _{j + 1/2}^n}}{{a_{j + 1/2}^{\rm{v}}}} \\ b_j^{\rm{p}} = b_{j - 1}^{\rm{p}} + b_{j + 1}^{\rm{p}} \\ d_j^{\rm{p}} = \rho _{j - 1/2}^nv_{j - 1/2}^*{A_{j - 1/2}} - \rho _{j + 1/2}^nv_{j + 1/2}^*{A_{j + 1/2}} \\ \end{gathered} \right. $ (12)

除了建立速度及压力的修正方程外，在SIMPLER算法中，还需要求解假拟速度，并通过假拟速度求解待修正的压力场。假拟速度通过上一迭代步的速度场 $ {v}^{0} $计算得到，如式（13）所示，其中 $ \hat {v} $为假拟速度, $ {\left({c}_{j-1/2}^{v}\right)}^{p} $为动量离散方程中 $ {c}_{j-1/2}^{v} $项里除去压力项后剩余的部分。

$ {\hat v_{j - 1/2}} = \dfrac{{a_{j - 5/2}^{\rm{v}}v_{j - {5}/{2}}^0 + a_{j - 3/2}^{\rm{v}}v_{j - {3}/{2}}^0}}{{a_{j - 1/2}^{\rm{v}}}} + \dfrac{{a_{j + 1/2}^{\rm{v}}v_{j + {1}/{2}}^0 + a_{j + 3/2}^{\rm{v}}v_{j + {3}/{2}}^0 + {{\left( {c_{j - 1/2}^{\rm{v}}} \right)}^p}}}{{a_{j - 1/2}^{\rm{v}}}} $ (13)

将上一迭代步的速度场 $ {v}^{0} $代入到动量离散方程，并将其与连续性方程联立求解，可以得到用于计算待修正压力场的压力计算方程

$ \left\{ \begin{gathered} b_j^{\rm{p}}p_j' = b_{j{{ - 1}}}^{\rm{p}}p_{j - 1}' + b_{j + 1}^{\rm{p}}p_{j + 1}' + d_j^{\rm{p}} \\ b_{j - 1}^{\rm{p}} = A_{j - 1/2}^2\dfrac{{\rho _{j - 1/2}^n}}{{a_{j - 1/2}^{\rm{v}}}} \\ b_{j + 1}^{\rm{p}} = A_{j + 1/2}^2\dfrac{{\rho _{j + 1/2}^n}}{{a_{j + 1/2}^{\rm{v}}}} \\ b_j^{\rm{p}} = b_{j - 1}^{\rm{p}} + b_{j + 1}^{\rm{p}} \\ d_j^{\rm{p}} = \rho _{j - 1/2}^n{{\hat v}_{j - 1/2}}{A_{j - 1/2}} - \rho _{j + 1/2}^n{{\hat v}_{j + 1/2}}{A_{j + 1/2}} \\ \end{gathered} \right. $ (14)

1.2 燃料棒传热

在整个反应堆中，燃料棒是温差最大、温度变化最为剧烈的地方，为了精确模拟其内部温度分布及变化，在主回路CFD网格的基础上，采用沿径向及轴向的二维网格对燃料棒空间进行划分。在燃料棒中，其芯块、气隙及包壳内部的传热为热传导，包壳与冷却剂之间的传热为热对流^[9]，其能量方程为

$ \dfrac{{\partial \rho {c_p}T}}{{\partial \tau }} = \dfrac{1}{r}\dfrac{\partial }{{\partial {\textit{z}}}}\left( {{\lambda _{\rm{m}}}r\dfrac{{\partial T}}{{\partial r}}} \right) + \dfrac{\partial }{{\partial {\textit{z}}}}\left( {{\lambda _{\rm{m}}}\dfrac{{\partial T}}{{\partial {\textit{z}}}}} \right) + {q_{{v}}}\left( {r,{\textit{z}}} \right) $ (15)

式中： $\lambda _{\rm{m}} $为燃料棒的导热系数。

同样的，我们需要对燃料棒的空间区域进行网格划分，并在此网格的基础上对能量方程进行离散，如式（16）所示，式中 $ \Delta {\textit{z}} $为控制体轴向高度， $ \Delta r $为控制体径向长度。当然，在燃料棒的轴向以及径向上会出现多个边界，在这里，我们认为包壳与冷却剂接触的边界为对流传热边界，其与冷却剂之间的换热由式（17）表示，式中 $ h $为包壳与冷却剂之间的对流换热系数； ${T}_{{\rm{LBE}}}$为在轴向上与包壳控制体所对应的冷却剂控制体的温度； $ {T}_{{{N}}+1/2} $为包壳与冷却剂接触的边界即包壳外表面温度。除此之外，其他边界与外部区域的传热几乎可以忽略，因此认为它们都是绝热边界。

$ \begin{split} &\Delta {{\textit{z}}_{i,j}}\rho _{i,j}^n({c_p})_{i,j}^n(T_{i,j}^{n + 1} - T_{i,j}^n)\dfrac{{r_{i + 1/2,j}^2 - r_{i - 1/2,j}^2}}{2} = \Delta {{\textit{z}}_{i,j}}\Bigg[\lambda _{i,j}^n{r_{i + 1/2}}\dfrac{{T_{i + 1/2,j}^{n + 1} - T_{i,j}^{n + 1}}}{{{r_{i,j}}/2}} - \lambda _{i,j}^n{r_{i - 1/2}}\dfrac{{T_{i,j}^{n + 1} - T_{i - 1/2,j}^{n + 1}}}{{{r_{i,j}}/2}}\Bigg]\Delta t + \\ &\quad \quad {r_{i,j}}\Delta {r_{i,j}}\left[ {\lambda _{i,j}^n\dfrac{{T_{i,j + 1/2}^{n + 1} - T_{i,j}^{n + 1}}}{{\Delta {{\textit{z}}_{i,j}}/2}} - \lambda _{i,j}^n\dfrac{{T_{i,j}^{n + 1} - T_{i,j - 1/2}^{n + 1}}}{{\Delta {{\textit{z}}_{i,j}}/2}}} \right]\Delta t + \Delta {{\textit{z}}_{i,j}}({q_{{v}}})_{i,j}^n\dfrac{{r_{i + 1/2,j}^2 - r_{i - 1/2,j}^2}}{2} \end{split} $ (16)

$ \lambda _{N,j}^n\left( {\dfrac{{\partial T}}{{\partial r}}} \right)_{N + 1/2,j}^{n + 1} = {h_j}\left[ {T_{{\rm{LBE}},j}^{n + 1} - T_{N + 1/2,j}^{n + 1}} \right] $ (17)

1.3 换热器传热

在换热器中同时流动着一回路及二回路的冷却剂，存在两个回路之间的换热。在1.1中，我们所开发的CFD模型是针对铅铋堆一回路的流动传热，因此需要对二回路的流动传热进行单独计算。在这里，由于铅铋堆还处于研究阶段，换热器大多采用的是单相管壳式换热器，其内部各回路冷却剂流动方式如图3所示。对换热器中的二回路冷却剂进行计算时，认为其是一维单相流动的，且计算区域为开式（从换热器二回路进口到出口），入口处流量保持不变。因此只需要计算传热，流速由固定不变的质量流量给出。因此，其控制方程及离散形式与式（3）、（5）一致。二者之间的换热关系如式（18）所示，式中 $ ：Q $为二者之间传递的热量； $ K $为二者之间的换热系数； $ S $为二者之间的传热面积； $ {d}_{{\rm{i}}} $、 $ {d}_{{\rm{o}}} $分别为内、外管内径； $ {h}_{1} $、 $ {h}_{2} $分别为一、二回路对流换热系数； $ \lambda $为换热管壁热导率^[11]。

图 3. 管壳式换热器

Fig. 3. Shell and tube heat exchanger

下载图片 查看所有图片

$ \left\{ \begin{gathered} Q = KS\left( {{T_{{\rm{LBE}}}} - {T_{\rm{W}}}} \right) \\ K = \dfrac{1}{{\dfrac{1}{{{h_1}}}\dfrac{{{d_{\rm{o}}}}}{{{d_{\rm{i}}}}} + \dfrac{{{d_{\rm{o}}}}}{{2\lambda }}{\text{ln}}\left( {\dfrac{{{d_{\rm{o}}}}}{{{d_{\rm{i}}}}}} \right) + \dfrac{1}{{{h_2}}}}} \\ \end{gathered} \right. $ (18)

式中： $T_{\rm{W}} $代表二回路工质水的温度。

1.4 点堆动力学

为了进行瞬态模拟，需要对中子动力学进行计算，由于铅铋堆属于快堆，因而采用忽略中子空间效应、只考虑时间效应的点堆动力学模型，其方程如式（19）所示。对于上式的求解，在这里采用一阶泰勒多项式对 $ n\left(t\right) $进行展开的半解析积分算法进行求解^[12]，如式（20）~（26）所示。

$ \left\{\begin{array}{l} \dfrac{{\rm{d}} n(t)}{{\rm{d}} t}=\dfrac{[\rho(t)-\beta]}{{\mathit{\Lambda}}} n(t)+\displaystyle\sum_{i=1}^{6} \lambda_{i} C_{i}(t)+q(t) \\ \dfrac{{\rm{d}} C_{i}(t)}{{\rm{d}} t}=\dfrac{\beta_{i}}{{\mathit{\Lambda}}} n(t)-\lambda_{i} C_{i}(t),\;\;i=1,2, \cdots, 6 \end{array}\right. $ (19)

$ n({t_{n + 1}}) = \dfrac{{\left\{ n({t_n}) + q({t_{n + 1}})h + \displaystyle\sum\limits_{i = 1}^6 {{C_i}({t_n})} - \displaystyle\sum\limits_{i = 1}^6 {{\lambda _i}\exp ( - {\lambda _i}h){C_i}({t_n})} \right\} + \dfrac{{\left\{ {{F_2} - \displaystyle\sum\limits_{i = 1}^6 {\dfrac{{{\beta _i}}}{{{\mathit{\Lambda}} }}{G_{2,i}}} } \right\}\left\{ {\displaystyle\sum\limits_{i = 1}^6 {{\lambda _i}\exp ( - {\lambda _i}h){C_i}({t_n})} + q({t_{n + 1}})} \right\}}}{{1 - \displaystyle\sum\limits_{i = 1}^6 {\dfrac{{{\lambda _i}{\beta _i}}}{{{\mathit{\Lambda}} }}{G_{2,i}}} }}}}{{\left\{ {1 - {F_1} + \displaystyle\sum\limits_{i = 1}^6 {\dfrac{{{\beta _i}}}{{{\mathit{\Lambda}} }}{G_{1,i}}} } \right\} - \left\{ {{F_2} - \displaystyle\sum\limits_{i = 1}^6 {\dfrac{{{\beta _i}}}{{{\mathit{\Lambda}} }}{G_{1,i}}} } \right\}\left\{ {\dfrac{{\rho ({t_{n + 1}}) - \beta }}{{{\mathit{\Lambda}} }} + \displaystyle\sum\limits_{i = 1}^6 {\dfrac{{{\lambda _i}{\beta _i}}}{{{\mathit{\Lambda}} }}} {G_{1,i}}} \right\}\left( {1 - \displaystyle\sum\limits_{i = 1}^6 {\dfrac{{{\lambda _i}{\beta _i}}}{{{\mathit{\Lambda}} }}} {G_{2,i}}} \right)}} $ (20)

$ n^{\prime}\left(t_{n+1}\right)=\dfrac{\left(\dfrac{\rho\left(t_{n+1}-\beta\right)}{{\mathit{\Lambda}}}+\displaystyle\sum_{i=1}^{6} \dfrac{\lambda_{i} \beta_{i}}{{\mathit{\Lambda}}} G_{1, i}\right) n\left(t_{n+1}\right)+\left(\displaystyle\sum_{i=1}^{6} \lambda_{i} \exp\left(-\lambda_{i} h\right) C_{i}\left(t_{n}\right)+q\left(t_{n+1}\right)\right)}{1-\displaystyle\sum_{i=1}^{6} \dfrac{\lambda_{i} \beta_{i}}{{\mathit{\Lambda}}} G_{2, i}} $ (21)

$ {C_i}\left( {{t_{n + 1}}} \right) = {\text{exp}}\left( { - {\lambda _i}h} \right){C_i}\left( {{t_n}} \right) + \dfrac{{{\beta _i}}}{{{\mathit{\Lambda}} }}\left[ {{G_{1,i}}n\left( {{t_{n + 1}}} \right) + {G_{2,i}}n'\left( {{t_{n + 1}}} \right)} \right] $ (22)

$ {F_1} = \displaystyle\int_{{t_n}}^{{t_{n + 1}}} {\dfrac{{\rho \left( \tau \right)}}{{\mathit{\Lambda}} }} {\rm{d}}\tau = \dfrac{{h\left[ {\rho \left( {{t_{n + 1}}} \right) + \rho \left( {{t_n}} \right)} \right]}}{{2\Delta }} $ (23)

$ {F_2} = \displaystyle\int_{{t_n}}^{{t_{n + 1}}} {\dfrac{{\rho \left( \tau \right)}}{{\mathit{\Lambda}} }} \left( {\tau - {t_{n + 1}}} \right){\rm{d}}\tau = - \dfrac{{{F_1}h}}{2} $ (24)

$ {G_{1,i}} = \exp \left( { - {\lambda _i}h} \right)\displaystyle\int_{{t_n}}^{{t_{n + 1}}} {\exp } \left[ {{\lambda _i}\left( {\tau - {t_n}} \right)} \right]{\rm{d}}\tau $ (25)

$ {G_{2,i}} = \exp \left( { - {\lambda _i}h} \right)\displaystyle\int_{{t_n}}^{{t_{n + 1}}} {\left( {\tau - {t_{n + 1}}} \right)\exp } \left[ {{\lambda _i}\left( {\tau - {t_n}} \right)} \right]{\rm{d}}\tau $ (26)

1.5 闭式并联多通道及多环路

闭式并联多通道及多环路本质上是属于多维的热工水力问题，使用集**数法的一维CFD并不能得到每个通道冷却剂的温度。但在对燃料棒传热进行计算时，又需要我们知道每个通道冷却剂的温度。为了解决上述问题，我们认为每个通道任意时刻的冷却剂流量占堆芯总流量的份额保持不变。当主回路的一维CFD模块在模拟堆芯时，将堆芯所有通道用集**数法等效为一个热力学平均管，其与燃料棒间的对流换热等于各通道冷却剂与对应燃料棒间对流换热的总和。对这个等效的管道进行一维的流动及传热模拟，计算得到当前时刻的温度及流速，并用于下一时刻堆芯热力学平均管的流动及传热模拟。除此之外，堆芯每一个通道的冷却剂流速由相同时刻的热力学平均管的冷却剂流速计算得出。通过堆芯热力学平均管当前时刻各节点的流速及每个通道占总流量的份额，可以得到每一个通道各网格节点的流速。然后在CFD网格的基础上，通过每一个通道各个节点的流速，单独求解能量方程，即式（3），得到当前时刻各个通道燃料及冷却剂的温度，并将这个温度用于下一时刻各个通道的传热模拟。即堆芯热力学平均管只用于求解堆芯冷却剂流速，而燃料棒及冷却剂温度则由单独求解能量方程所给出。并且，在单独计算各通道能量方程时，与每一个通道相对应的燃料棒联立计算，即通过式（17）将二者进行隐式耦合计算，同时得到通道冷却剂及燃料棒的温度信息。

除堆芯多通道外，如果在主回路中有多台换热器，就可能会出现每台换热器换热能力不同的情况，这时采用集**数法同样也不够准确，因而对于每一个环路上的换热器采用与前面所述的闭式并联多通道一致的模拟方法，这里不再另述。特别的，在单独计算每台换热器的一、二回路的能量方程时（二者的能量方程均与式（3）一致），将二者联立，通过式（18）进行隐式耦合计算。

2 程序耦合

根据上述物理、数学模型，我们将程序最为核心的部分分为四个模块：计算燃料棒传热及通道冷却剂传热的堆芯模块；基于集**数法及热力学平均管，计算主回路冷却剂流动及传热的CFD模块；计算换热器一、二回路传热的换热器模块；计算中子动力学的点堆动力学模块。它们之间的耦合过程如图4所示。

图 4. 程序主要模块耦合

Fig. 4. Main coupling module of the code

下载图片 查看所有图片

首先，在堆芯模块与CFD模块之间的耦合中，前者为后者提供了各通道的包壳表面温度及冷却剂温度用以计算堆芯热力学平均管与燃料棒之间的换热，而后者为前者提供了基于堆芯热力学平均管的流动情况，并通过流量分配得出各个通道的流速，用于计算传热；其次，在堆芯模块与点堆动力学模块的耦合中，堆芯模块为后者提供了各通道燃料棒的温度及冷却剂通道温度，用于燃料及冷却剂反馈的计算，而后者又为堆芯模块提供了中子密度用于计算热源；最后，在换热器模块与CFD模块的耦合中，前者为后者提供了各环路上换热器一、二回路温度用以计算换热器一回路的热力学平均管与二回路之间的换热，而后者为前者提供了每一台换热器一回路的冷却剂流速，以计算各个换热器一、二回路的传热。另外，在程序中，将泵简化为单个节点，泵压头为泵节点提供外部压头，若压头为零，则为自然循环。

3 程序验证

3.1 计算模型

为了验证本研究开发的系统安全分析程序的可靠性与准确性，我们采用OECD/NEA发布的ADS失束事故国际基准例题及 D’Angelo等人所发布的基准测试报告^[9,13]，对本程序进行验证。为了验证程序整体耦合，我们建立了一个完整的池式铅铋冷却反应堆主回路系统，其结构如图5（a）所示，包括了堆芯、热池、冷池、下腔室、泵及泵管道、换热器。在这个系统中，我们在堆芯处设置了3个通道，其中，我们关注的是1、2通道，每个通道的功率及流量分配如表1所示，并且我们认为各个通道的燃料棒功率在轴向上的轮廓相同。需要说明的是，设置的两台换热器的一、二回路流量都分别相同。

图 5. XADS堆芯结构及计算模型

Fig. 5. XADS core structure and calculation model

下载图片 查看所有图片

在计算模拟中，我们将主回路系统等效为图5（b）所示的计算模型，将冷池、热池及下腔室等效为一个圆管道进行计算。为了与基准保持一致，燃料棒线功率密度、点堆动力学参数、燃料棒与冷却剂之间的对流换热系数、燃料及冷却剂的反馈关系以及LBE冷却剂、芯块、包壳、气隙的物性均采用基准报告^[9,13]里的数据。

3.2 结果分析

3.2.1 稳态计算结果

首先使用稳态模拟来验证堆芯传热与CFD之间的耦合，在这里，我们使用Origin工具将本文开发程序计算所得数据添加到基准报告^[13]的数据中，如图6所示。其中，图6（a）为通道1（平均通道），图6（b）为通道2（热通道）。从对比图来看，燃料芯块温度呈近余弦分布，且本文开发程序所得结果无论是在数值还是趋势上，均与基准报告中10款通过验证的系统分析程序计算结果符合较好。且从计算结果来看，通道1、2计算数据均在基准报告^[13]中的10款程序的最大误差（1.8%、3%）以内。图7（a）、（b）分别为通道1及通道2燃料棒的二维剖面稳态温度分布。

图 6. 堆芯稳态轴向温度分布

Fig. 6. Steady-state axial temperature distribution of core

下载图片 查看所有图片

图 7. 燃料棒剖面稳态温度分布

Fig. 7. Steady-state temperature distribution of fuel profile

下载图片 查看所有图片

3.2.2 瞬态计算结果

对于瞬态模拟的验证，我们使用XADS反应堆的一个典型的瞬态事故：失束事故。此瞬态事故的条件为：在反应堆满功率运行下，将加速器束流切断，失束时间分别1、3、6、12 s，然后恢复加速器束流，观察堆芯功率及各关键节点温度的变化。使用本开发程序对上述瞬态工况进行模拟，并将计算数据与基准报告对比。

首先是通道1，在失束时间分别为1、3、6、12 s时，归一化功率与基准报告[9]中的数据对比如图8所示；燃料中平面温度及通道活性区冷却剂出口温度与基准报告^[9]中的数据对比如图9所示，其中图9（a）为燃料芯块中平面中心温度，图9（b）为燃料芯块中平面外表面温度，图9（c）为活性区出口温度。从对比图可以看出，在束流切断后，堆芯功率及温度迅速下降，且功率快速降至初始时的10%左右，当束流恢复后又逐渐上升至初始值。从计算数据来看，失束1、3、6、12 s时，燃料中平面中心下降温度与基准报告^[9]中结果最大偏差分别为2.9%、0.5%、4.6%、0.8%，活性区冷却剂出口的最大偏差分别为1.5%、2.8%、1.6%、1.2%，如表2所示。

图 8. 归一化功率

Fig. 8. Normalized power

下载图片 查看所有图片

图 9. 平均通道燃料中平面温度及活性区冷却剂出口温度

Fig. 9. Average pin fuel temperature at fuel zone midplane and outlet coolant temperature

下载图片 查看所有图片

由于基准报告^[13]中只记录了热通道在失束时间为1 s、6 s时的温度演变数据，因此对于通道2，选取时间为1 s、6 s的数据与基准报告进行对比，如图10所示，其中图10（a）为燃料芯块中平面中心温度，图10（b）为活性区出口温度。从对比图来看，热通道的计算结果与基准报告中结果的误差较平均通道稍大，但也保持在较小且可接受的范围内。从计算数据来看，失束1 s、6 s时，燃料中平面中心下降温度与基准报告^[13]中的结果最大偏差分别为4.3%、2.9%，活性区冷却剂出口的最大偏差分别为5%、5%，如表3所示。

图 10. 热通道燃料中心中平面温度及活性区冷却剂出口温度

Fig. 10. Hottest pin fuel centerline temperature at fuel zone midplane and outlet coolant temperature

下载图片 查看所有图片

4 结　论

由于多维计算程序在对铅铋堆主回路进行长时间的精确模拟时，消耗的计算资源较为庞大，且大型商用计算程序存在不能灵活修改的问题。因而本文选择了对铅铋堆系统分析程序的自主开发进行了初步研究。在本研究中，自主开发了一维CFD程序，并将二维燃料棒传热模型及点堆动力学模型耦合在其中，初步形成了一款适用于铅铋堆的系统分析程序。为了确保开发模型的准确性，使用了OECD/NEA发布的ADS失束事故国际基准例题进行程序验证。验证结果显示，各项关键参数在数值及趋势上均与国际上的同类程序计算结果吻合较好，且计算结果最大误差均保持在5%以内，表明本文开发模型具有较高的可信度。与多维程序以及商用程序相比，该程序计算资源的消耗小且建模灵活，可对铅铋堆堆芯热工水力及安全问题进行初步模拟分析。