一种超临界水堆堆芯三维瞬态性能分析方法

摘要

本发明公开了一种超临界水堆堆芯三维瞬态性能分析方法包括：步骤1：执行堆芯稳态计算，提供功率参数、热工参数，执行迭代耦合计算直至堆芯稳态功率参数和热工参数收敛；步骤2：提供堆芯初始状态和组件截面库；步骤3：执行瞬态计算，提供功率参数、热工参数，执行迭代耦合计算直至堆芯瞬态功率参数和热工参数收敛；步骤4：得到热组件栅元尺度的精细功率分布，进行热组件子通道分析，最终给出安全评价关键参数；步骤5：判断所述安全评价关键参数是否超过相应的瞬态或事故安全限值，实现了能够准确描述物理与热工-水力之间的反馈，提供较精确的堆芯三维功率分布，真实地模拟超临界水堆的瞬态过程和事故过程的技术效果。

说明书

技术领域

本发明涉及核反应堆设计研究领域，尤其涉及一种超临界水堆堆芯三维瞬态 性能分析方法。

背景技术

超临界水堆(SCWR)的运行工况在水的热力学临界点(374℃，22.1MPa)以 上，系统热效率高，经济性好，但同时，与常规水冷堆相比，超临界水堆运行参 数大幅提高，并且堆芯结构和冷却剂流程相对复杂，带来大量的堆芯物理、热工 -水力以及安全分析等技术难题，其中，强烈的核热耦合特性是超临界水堆堆芯 设计与性能分析的关键技术难题之一。

超临界水堆的反应性异常事件中，大多数与控制棒动作相关，控制棒动作异 常事件将引起超临界水堆堆芯功率分布畸变，伴随着强烈的物理与热工-水力耦 合效应，功率分布变化将直接影响超临界水堆燃料最大包壳温度，而包壳温度不 超限是超临界水堆非常重要的瞬态安全分析准则，对于功率分布异常事件，点堆 或一维的中子动力学模型无法正确描述功率的空间分布，无法准确计算功率随时 间的变化，为保证计算的包络性，就必须引入大量的保守性假设，此外，超临界 水堆由于本身就具备强核热耦合特性，其反应性及功率分布异常事件的瞬态过程 更加复杂，而保守方法不能提供其真实的响应过程，这使得设计人员难以对反应 性异常事件进行深入的研究和分析，只有耦合三维中子时空动力学与热工-水力 学进行瞬态过程分析，才能准确描述物理与热工-水力之间的反馈，提供较精确 的三维功率分布，真实地模拟超临界水堆的瞬态过程和事故过程。

现有技术中，对于超临界水堆功率分布异常事件，国际上采用的瞬态分析方 法，其中子动力学部分采用点堆或一维模型，并不具备堆芯核热耦合三维功率计 算能力，难以真实或精确地模拟例如反应性及功率分布异常事件的瞬态过程。

综上所述，本申请发明人在实现本申请实施例中发明技术方案的过程中，发 现上述技术至少存在如下技术问题：

在现有技术中，现有的超临界水堆功率分布瞬态分析方法存在不具备堆芯核 热耦合三维功率计算能力，难以真实或精确地模拟例如反应性及功率分布异常事 件的瞬态过程的技术问题。

发明内容

本发明提供了一种超临界水堆堆芯三维瞬态性能分析方法，解决了现有的超 临界水堆功率分布瞬态分析方法存在不具备堆芯核热耦合三维功率计算能力，难 以真实或精确地模拟例如反应性及功率分布异常事件的瞬态过程的技术问题，实 现了能够准确描述物理与热工-水力之间的反馈，提供较精确的堆芯三维功率分 布，真实地模拟超临界水堆的瞬态过程和事故过程的技术效果。

为解决上述技术问题，本申请实施例提供了一种超临界水堆堆芯三维瞬态性 能分析方法，所述方法包括：

步骤1：超临界水堆堆芯三维稳态物理-热工水力耦合模块执行堆芯稳态计 算，堆芯燃料管理模块向子通道热工-水力模块提供功率参数，子通道热工-水力 模块向堆芯燃料管理模块提供热工参数，执行迭代耦合计算直至堆芯稳态功率参 数和热工参数收敛；

步骤2：堆芯三维稳态物理-热工水力耦合模块向超临界水堆堆芯三维瞬态 物理-热工水力耦合模块提供堆芯初始状态和组件截面库；

步骤3：超临界水堆堆芯三维瞬态物理-热工水力耦合模块执行瞬态计算， 三维中子时空动力学模块向子通道热工-水力模块提供功率参数，子通道热工- 水力模块向三维中子时空动力学模块提供热工参数，执行迭代耦合计算直至堆芯 瞬态功率参数和热工参数收敛；

步骤4：在堆芯计算中对热组件进行功率重构得到热组件栅元尺度的精细功 率分布，再次利用子通道热工-水力模块进行热组件子通道分析，最终给出安全 评价关键参数；

步骤5：判断所述安全评价关键参数是否超过相应的瞬态或事故安全限值， 评价超临界水堆堆芯在瞬态过程下的安全性能。

进一步的，所述安全评价关键参数包括：最大包壳壁面温度、芯块热焓。

进一步的，在堆芯三维中子学计算中使用少群截面，少群截面关于热工参数 及预设参数的变化采用分段插值的方式预先加工成截面库供堆芯三维中子学模 块使用，所述预设参数包括：燃料组件燃耗、有无控制棒。

进一步的，所述执行迭代耦合计算直至堆芯稳态功率参数和热工参数收敛具 体包括：

开始；

读取输入文件，初始化组件少群截面、堆芯中子通量；

进行子通道热工模块稳态计算；

进行更新截面；

进行三维中子学模块稳态计算；

进行子通道热工模块稳态计算；

判断堆芯功率分布是否收敛，若否则返回步骤重新进行更新截面；若是则稳 态计算结束，进入瞬态计算。

进一步的，所述超临界水堆堆芯三维瞬态物理-热工水力耦合模块执行瞬态 计算具体包括：

开始瞬态计算；

进行时间步更新，准备n时刻的初值；

进行子通道热工模块瞬态计算一个时间步；

进行更新截面；

进行三维中子时空动力学模块瞬态计算第一个时间步；

判断堆芯功率分布是否收敛，若否则返回执行时间步更新，准备n时刻的初 值；若是则n＝n+1；判断tn是否小于T，其中，tn为第n时步对应时间，T为预 先指定的瞬态计算至T时刻，若是则返回执行时间步更新，准备n时刻的初值； 若否则计算结束。

进一步的，所述超临界水堆堆芯三维稳态物理-热工水力耦合模块执行堆芯 稳态计算和所述超临界水堆堆芯三维瞬态物理-热工水力耦合模块执行瞬态计 算，均采用相同的空间求解方法：第二类边界条件节块格林函数方法。这样以保 证三维稳态-瞬态计算分析过程中中子学计算的自洽性。

进一步的，所述执行堆芯三维瞬态物理计算的三维中子时空动力学模块采用 的时间离散方法为向后欧拉方法。其中，向后欧拉方法具备绝对-稳定性，并且 便于程序的串行耦合开发。

本申请实施例中提供的一个或多个技术方案，至少具有如下技术效果或优 点：

由于采用了将超临界水堆堆芯三维瞬态性能分析方法设计为包括：步骤1： 超临界水堆堆芯三维稳态物理-热工水力耦合模块执行堆芯稳态计算，堆芯燃料 管理模块向子通道热工-水力模块提供功率参数，子通道热工-水力模块向堆芯燃 料管理模块提供热工参数，执行迭代耦合计算直至堆芯稳态功率参数和热工参数 收敛；步骤2：堆芯三维稳态物理-热工水力耦合模块向超临界水堆堆芯三维瞬 态物理-热工水力耦合模块提供堆芯初始状态和组件截面库；步骤3：超临界水 堆堆芯三维瞬态物理-热工水力耦合模块执行瞬态计算，三维中子时空动力学模 块向子通道热工-水力模块提供功率参数，子通道热工-水力模块向三维中子时空 动力学模块提供热工参数，执行迭代耦合计算直至堆芯瞬态功率参数和热工参数 收敛；步骤4：在堆芯计算中对热组件进行功率重构得到热组件栅元尺度的精细 功率分布，再次利用热工模块进行热组件子通道分析，最终给出安全评价关键参 数；步骤5：判断所述安全评价关键参数是否超过相应的瞬态或事故安全限值， 评价超临界水堆堆芯在瞬态过程下的安全性能的技术方案，即通过耦合三维中子 时空动力学与热工-水力学进行瞬态过程分析，中子动力学部分采用三维模型， 使其具备堆芯核热耦合三维功率计算能力，实现超临界水堆堆芯三维瞬态分析的 技术效果；所以，有效解决了现有的超临界水堆功率分布瞬态分析方法存在不具 备堆芯核热耦合三维功率计算能力，难以真实或精确地模拟例如反应性及功率分 布异常事件的瞬态过程的技术问题，实现了通过建立堆芯三维稳态-瞬态物理- 热工水力耦合计算流程，能够准确描述物理与热工-水力之间的反馈，提供较精 确的堆芯三维功率分布，真实地模拟超临界水堆的瞬态过程和事故过程的技术效 果。

进一步的，本方法可以真实地模拟超临界水堆堆芯瞬态过程，提高对超临界 水堆堆芯瞬态工况的认识，更加精确的估计安全裕量，改善堆芯设计并提高堆芯 性能。

附图说明

为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例和附图， 对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本 发明，并不作为对本发明的限定。

图1为超临界水堆堆芯三维稳态-瞬态物理-热工水力耦合计算流程示意图；

图2为超临界水堆堆芯三维瞬态性能分析方法稳态计算流程示意图；

图3为超临界水堆堆芯三维瞬态性能分析方法瞬态计算流程示意；

图4为向后欧拉方法瞬态中子学计算流程示意图。

具体实施方式

本发明提供了一种超临界水堆堆芯三维瞬态性能分析方法，解决了现有的超 临界水堆功率分布瞬态分析方法存在不具备堆芯核热耦合三维功率计算能力，难 以真实或精确地模拟例如反应性及功率分布异常事件的瞬态过程的技术问题，实 现了能够准确描述物理与热工-水力之间的反馈，提供较精确的堆芯三维功率分 布，真实地模拟超临界水堆的瞬态过程和事故过程的技术效果。

本申请实施中的技术方案为解决上述技术问题。总体思路如下：

采用了将超临界水堆堆芯三维瞬态性能分析方法设计为包括：步骤1：超临 界水堆堆芯三维稳态物理-热工水力耦合模块执行堆芯稳态计算，堆芯燃料管理 模块向子通道热工-水力模块提供功率参数，子通道热工-水力模块向堆芯燃料管 理模块提供热工参数，执行迭代耦合计算直至堆芯稳态功率参数和热工参数收 敛；步骤2：堆芯三维稳态物理-热工水力耦合模块向超临界水堆堆芯三维瞬态 物理-热工水力耦合模块提供堆芯初始状态和组件截面库；步骤3：超临界水堆 堆芯三维瞬态物理-热工水力耦合模块执行瞬态计算，三维中子时空动力学模块 向子通道热工-水力模块提供功率参数，子通道热工-水力模块向三维中子时空动 力学模块提供热工参数，执行迭代耦合计算直至堆芯瞬态功率参数和热工参数收 敛；步骤4：在堆芯计算中对热组件进行功率重构得到热组件栅元尺度的精细功 率分布，再次利用热工模块进行热组件子通道分析，最终给出安全评价关键参数； 步骤5：判断所述安全评价关键参数是否超过相应的瞬态或事故安全限值，评价 超临界水堆堆芯在瞬态过程下的安全性能的技术方案，即通过耦合三维中子时空 动力学与热工-水力学进行瞬态过程分析，中子动力学部分采用三维模型，使其 具备堆芯核热耦合三维功率计算能力，实现超临界水堆堆芯三维瞬态分析的技术 效果；所以，有效解决了现有的超临界水堆功率分布瞬态分析方法存在不具备堆 芯核热耦合三维功率计算能力，难以真实或精确地模拟例如反应性及功率分布异 常事件的瞬态过程的技术问题，实现了通过建立堆芯三维稳态-瞬态物理-热工水 力耦合计算流程，能够准确描述物理与热工-水力之间的反馈，提供较精确的堆 芯三维功率分布，真实地模拟超临界水堆的瞬态过程和事故过程的技术效果。

为了更好的理解上述技术方案，下面将结合说明书附图以及具体的实施方式 对上述技术方案进行详细的说明。

实施例一：

在实施例一中，提供了一种超临界水堆堆芯三维瞬态性能分析方法，请参考 图1-图4，所述方法包括：

步骤1：超临界水堆堆芯三维稳态物理-热工水力耦合模块执行堆芯稳态计 算，堆芯燃料管理模块向子通道热工-水力模块提供功率参数，子通道热工-水力 模块向堆芯燃料管理模块提供热工参数，执行迭代耦合计算直至堆芯稳态功率参 数和热工参数收敛；

步骤2：堆芯三维稳态物理-热工水力耦合模块向超临界水堆堆芯三维瞬态 物理-热工水力耦合模块提供堆芯初始状态和组件截面库；

步骤3：超临界水堆堆芯三维瞬态物理-热工水力耦合模块执行瞬态计算， 三维中子时空动力学模块向子通道热工-水力模块提供功率参数，子通道热工- 水力模块向三维中子时空动力学模块提供热工参数，执行迭代耦合计算直至堆芯 瞬态功率参数和热工参数收敛；

步骤4：在堆芯计算中对热组件进行功率重构得到热组件栅元尺度的精细功 率分布，再次利用热工模块进行热组件子通道分析，最终给出安全评价关键参数；

步骤5：判断所述安全评价关键参数是否超过相应的瞬态或事故安全限值， 评价超临界水堆堆芯在瞬态过程下的安全性能。

其中，在本申请实施例中，所述安全评价关键参数包括但不限于：最大包壳 壁面温度、芯块热焓。

其中，在实际应用中，上述安全评价关键参数对应的安全限值分别为：最大 包壳壁面温度瞬态工况下安全限值为850℃，事故工况下安全限值为1260℃；芯 块热焓瞬态工况下安全限值为170cal/g，事故工况下安全限值为230cal/g。(事故 工况下，安全准则要求不出现过度的堆芯损坏；瞬态工况下，安全准则要求不出 现系统性燃料棒损坏、燃料芯块损坏。)

其中，在本申请实施例中，在堆芯计算中使用少群截面，少群截面关于热工 参数及预设参数的变化采用分段插值的方式预先加工成截面库供堆芯模块使用， 所述预设参数包括但不限于：燃耗、有无控制棒。

其中，在本申请实施例中，所述执行迭代耦合计算直至堆芯稳态功率参数和 热工参数收敛具体包括：

开始；

读取输入文件，初始化组件少群截面、堆芯中子通量；

进行子通道热工模块稳态计算，根据初始组件少群截面、堆芯中子通量，计 算堆芯内慢化剂密度和冷却剂密度分布；

进行更新截面，根据慢化剂密度和冷却剂密度分布，插值计算组件少群截面；

进行三维中子学模块稳态计算，根据新的组件少群截面，计算获得组件平均 功率分布，更新堆芯中子通量；

进行子通道热工模块稳态计算，根据新的组件少群截面、堆芯中子通量，计 算堆芯内慢化剂密度和冷却剂密度分布；

判断功率分布是否收敛，若否则返回步骤重新进行更新截面；若是则稳态计 算结束，进入瞬态计算。

其中，在本申请实施例中，所述超临界水堆堆芯三维瞬态物理-热工水力耦 合模块执行瞬态计算具体包括：

开始瞬态计算；

进行时间步更新，准备n时刻的初值；

进行子通道热工模块瞬态计算一个时间步，根据初始组件少群截面、堆芯功 率分布，计算堆芯内慢化剂密度和冷却剂密度分布；

进行更新截面，根据慢化剂密度和冷却剂密度分布，插值计算组件少群截面；

进行三维中子时空动力学模块瞬态计算第一个时间步，根据新的组件少群截 面，计算获得堆芯功率分布；

判断功率分布是否收敛，若否则返回执行时间步更新，准备n时刻的初值； 若是则n＝n+1，判断第n时步对应时间tn是否小于预先指定的瞬态计算时间T， 若是则返回执行时间步更新，准备n时刻的初值；若否则计算结束。

进一步的，所述超临界水堆堆芯三维稳态物理-热工水力耦合模块执行堆芯 稳态计算和所述超临界水堆堆芯三维瞬态物理-热工水力耦合模块执行瞬态计 算，均采用相同的空间求解方法：第二类边界条件节块格林函数方法。这样以保 证三维稳态-瞬态计算分析过程中中子学计算的自洽性。

第二类边界条件节块格林函数方法求解思路如下：

待求解的三维中子时空动力学方程形式为：

$\left(\begin{array}{c}\frac{1}{v_g}\frac{\partial {\phi }_g(\overrightarrow{r},t)}{\partial t}-▿·D_g▿{\phi }_g(\overrightarrow{r},t)=\underset{g^{\prime }\ne g}{\Sigma }{\Sigma }_{s,g^{\prime }\to g}{\phi }_{g^{\prime }}(\overrightarrow{r},t)+{\chi }_g(1-\beta )\underset{g^{\prime }=1}{\overset{G}{\Sigma }}v{\Sigma }_{f,g^{\prime }}{\phi }_{g^{\prime }}(\overrightarrow{r},t)\\ +\underset{i=1}{\overset{ND}{\Sigma }}{x}_{g,j}{\lambda }_i{C}_i(\overrightarrow{r},t),g=1,\mathrm{...},G\\ \frac{\partial {\phi }_g(\overrightarrow{r},t)}{\partial t}={\beta }_i\underset{g^{\prime }=1}{\overset{G}{\Sigma }}{\mathrm{v\Sigma }}_{f,g^{\prime }}{\phi }_g(\overrightarrow{r},t)-{\lambda }_i{C}_i(\overrightarrow{r},t),i=1,\mathrm{...},ND\end{array}\right)---\left(1\right)$

其中，g表征能群，能群总数为G；v_g表征第g群中子速度，单位cm·s^-1； 表征第g群中子通量，单位cm^-2·s^-1；D_g表征第g群扩散系数，单位cm； Σ_r,g和Σ_f,g分别表征第g群移出截面和裂变截面，Σ_s,g′→g表征第g′群到第g群的散 射截面，单位cm^-1；χ_g表征瞬发中子裂变谱；ν表征每次裂变释放中子数；i表征 缓发中子先驱核分组，先驱核组总数为ND；χ_g,i表征第i组缓发中子谱分额；λ_i表征第i组缓发中子先驱核衰变常数，单位s^-1；表征第i组缓发中子先驱 核浓度，单位cm^-3；β_i表征第i组缓发中子份额。

为求解瞬态方程组(1)，可记Δt＝t-t₀，y₁＝y(t)，y₀＝y(t₀)，其中t表征当前待求解时刻，t₀表征上一时刻。 f(t,y)表征方程组(1)中非时间偏导数项的总和，则方程组(1)可简写为：

$\frac{\partial y}{\partial t}=f(t,y)---\left(2\right)$

方程的时间离散采用A-稳定向后欧拉格式，可得：

y₁＝y₀+Δt×f(t,y₁)(3)

应用向后欧拉格式，先驱核浓度方程离散为：

类似地，中子通量方程离散为：

将式(4)、式(5)代入中子扩散方程，得到式(6)，问题转化为一个固定源问题：

$-▿·D_g▿{\phi }_g(\overrightarrow{r},t)+{\Sigma }_{r,g}{\phi }_g(\overrightarrow{r},t)=-\frac{1}{v_g}\frac{{\phi }_g(\overrightarrow{r},t)}{\Delta t}+{\hat{Q}}_g(\overrightarrow{r},t)+S_g(\overrightarrow{r},t_0)---\left(6\right)$

固定源方程(6)采用基于第二类边界条件的节块格林函数方法求解，求解思 路如下：

1、引入第二类边界条件格林函数得到横向积分偏中子通量方程的积分解；

2、利用节块界面非均匀通量和流的连续性条件，得到界面平均净中子流的 相应矩阵；

3、偏中子通量、先驱核浓度、横向泄漏项均采用2阶勒让德多项式展开， 采用剩余权重法，求解偏中子通量展开系数；

4、求解节块中子平衡方程得到中子平均通量，进而修正偏中子通量展开系 数。

节块格林函数求解过程中，堆芯被划分成若干个核截面相同的基本计算单 元，称之为节块。将节块作为空间离散的最小单元，k表征节块编号，x，y，z 方向坐标的定义域为：

$x\in [-a_x^k,a_x^k],y\in [-a_y^k,a_y^k],z\in [-a_z^k,a_z^k]---\left(7\right)$

以x方向为例说明。沿y,z两个方向对式(6)在节块k内积分，即可得到x方 向偏中子通量方程式：

$-D_g\frac{{\partial }^2}{\partial x^2}{\phi }_{gx}^k(x,t)+{\Sigma }_{r,g}^k{\phi }_{gx}^k(x,t)=-\frac{1}{v_g}\frac{{\phi }_{gx}^k(x,t)}{\Delta t}+{\hat{Q}}_{gx}^k(x,t)-L_{gx}^k(x,t)+S_{gx}^k(x,t_0)---\left(8\right)$

式中，分别表征横向积分通量、 源项、固定源项、横向泄漏项。

引入格林函数，则得到偏中子通量方程的积分解：

$\left(\begin{array}{c}{\phi }_{gx}^k(x,t)+\underset{-a_x^k}{\overset{a_x^k}{\int }}\frac{1}{v_g}\frac{{\phi }_{gx}^k(x_0,t)}{\Delta t}{G}_{gx}^k(x,x_0){\mathrm{dx}}_0=\\ \underset{-a_x^k}{\overset{a_x^k}{\int }}[{\hat{Q}}_{gx}^k(x_0,t)-L_{gx}^k(x_0,t)+S_{gx}(x_0,t_0)]G_{gx}^k(x,x_0){\mathrm{dx}}_0\\ -G_{gx}^k(x,a_x^k)J_{gx}^k(a_x^k,t)+G_{gx}^k(x,-a_x^k)J_{gx}^k(-a_x^k,t)\end{array}\right)---\left(9\right)$

式中，表征基于第二类边界条件的格林函数。

非均匀通量和净中子流在相邻节块界面上连续。非均匀通量连续的表达式 为：

$f_{gx+}^k{\phi }_{gx}^k(a_x^k,t)=f_{gx-}^{k+1}{\phi }_{gx}^{k+1}(-a_x^{k+1},t)---\left(10\right)$

其中，表征k节块x坐标方向正端点的g群不连续因子。表征x坐 标方向k+1节块的负端点不连续因子。

将式(9)代入式(10)，同时利用净中子流连续条件，可得：

$\begin{array}{c}{f}_{gx+}^k{T}_{gx}^k{J}_{gx}^{k-1}(a_x^{k-1},t)-[f_{gx+}^k{R}_{gx}^k+f_{gx-}^{k+1}{R}_{gx}^{k+1}]J_{gx}^k(a_x^k,t)\\ +f_{gx-}^k{T}_{gx}^k{J}_{gx}^{k-1}(a_x^{k-1},t)=f_{gx-}^{k+1}{[GQ-]}_{gx}^{k+1}\left(t\right)-f_{gx+}^k{[GQ+]}_{gx}^k\left(t\right)\end{array}---\left(11\right)$

式(11)给出了在k节块x坐标方向内部界面上的净中子流耦合关系。在外部 界面处，中子流条件可由堆芯边界条件获得，进而得到该列节块的净中子流响应 矩阵。以堆芯x方向左边界为例：

1、由节块边界反射对称条件，可得：

$J_{gx}^1(-a_x^1,t)=0---\left(12\right)$

2、由入射流为0条件，可得：

$J_{gx}^{in,1}(-a_x^1,t)=\frac{1}{4}\left({\phi }_{gx}^1\right(-a_x^1,t)f_{gx-}^1+2·J_{gx}^1(-a_x^1,t))=0---\left(13\right)$

将式(9)代入式(13)，可得：

$(R_{gx}^1{f}_{gx-}^1+2)·J_{gx}^1(-a_x^1,t)-T_{gx}^1·J_{gx}^1(a_x^1,t)f_{gx-}^1=-2·{[GQ-]}_{gx}^1{f}_{gx-}^1---\left(14\right)$

3、由边界通量为0条件，可得：

${\phi }_{gx}^1(-a_x^1,t)f_{gx-}^1=0---\left(15\right)$

将式(9)代入式(15)，可得：

$R_{gx}^1·J_{gx}^1(-a_x^1,t)f_{gx-}^1-T_{gx}^1·J_{gx}^1(a_x^1,t)f_{gx-}^1=-2·{[GQ-]}_{gx}^1{f}_{gx-}^1+{\Phi }_{gx-}^1---\left(16\right)$

4、由90°旋转对称条件，可得：

$J_{gx}^1(-a_x^1,t)=-J_{gy}^1(-a_y^1,t)---\left(17\right)$

上述边界条件适用于堆芯径向，可将堆芯x方向和y方向节块统一求解。

中子平衡方程可通过在节块k内对式(6)体积分获得：

$\begin{array}{c}({\Sigma }_{r,g}^k+\frac{1}{v_g{\mathrm{\Delta t}}_j}){\bar{\phi }}_g^k\left(t\right)={\bar{S}}_g^k\left(t_0\right)+{\bar{\hat{Q}}}_g^k\left(t\right)\\ -\underset{u=x,y,z}{\Sigma }\frac{1}{2a_u^k}[J_{gu}^k(a_u^k,t)-J_{gu}^k(-a_u^k,t)]\end{array}---\left(18\right)$

式中，u＝x,y,z，分别表征节块体积平均通量、源项、 固定源项。

基于格林函数的偏中子通量方程积分解式(9)、净中子流耦合关系式(11)、基 于节块方法的中子平衡方程式(18)组成了偏积分通量分布、界面净流、节块平均 通量的方程组。以坐标u表征x，y，z坐标，偏通量、源、横向泄漏等一维分布 量通过2阶勒让德多项式展开，可得：

${\phi }_{gu}^k(u,t)=\underset{n=1}{\overset{3}{\Sigma }}{\phi }_{gun}^k{p}_{un}^k\left(u\right)$

${\hat{Q}}_{gu}^k(u,t)=\underset{n=1}{\overset{3}{\Sigma }}{\hat{Q}}_{gun}^k{p}_{un}^k\left(u\right)---\left(19\right)$

$L_{gu}^k(u,t)=\underset{n=1}{\overset{3}{\Sigma }}{L}_{gun}^k{p}_{un}^k\left(u\right)$

式中：

$p_{u1}^k\left(u\right)=1$

$\begin{array}{cc}{p}_{u2}^k\left(u\right)=u/a_u^k& u\in [-a_u^k,a_u^k]---\end{array}\left(20\right)$

$p_{u3}^k\left(u\right)=\frac{3}{2}{\left(\frac{u}{a_u^k}\right)}^2-\frac{1}{2}$

将(19)代入式(11)至式(17)，可得净流响应矩阵。采用剩余权重法，求解 偏中子通量展开系数，权重函数如(20)式所示。偏中子通量求解的矩阵形式为：

$\begin{array}{c}\left\{\right[A_u]+\frac{\left[G_{g,u}^{uu,k}\right]}{v_g{\mathrm{\Delta t}}_j}\}{\underset{‾}{\phi }}_{{}^{gu}}^k\left(t\right)=\left[G_{gu}^{uu,k}\right]\{\widehat{\underset{‾}{Q}}\left(t\right)-{\underset{‾}{L}}_{gu}^k\left(t\right)+{\underset{‾}{S}}_{gu}^k\left(t_0\right)\}\\ -{\underset{‾}{G}}_{gu}^{+u,k}{J}_{gu}^k(a_u^k,t)+{\underset{‾}{G}}_{gu}^{-u,k}{J}_{gu}^k(-a_u^k,t)\end{array}---\left(21\right)$

式中，下划线代表3维向量。表征展开系数向量。与同理可得。和分别表征和的 勒让德展开系数；[A_u]和为3×3方阵，其系数满足：

${\left({\underset{‾}{G}}_{gu}^{\pm u,k}\right)}_n={\int }_{-a_u^k}^{a_u^k}{\mathrm{duG}}_{gu}^k(u,\pm a_u^k)p_{un}^k\left(u\right)---\left(22\right)$

${\left[A_u\right]}_{mn}=\underset{-a_u^k}{\overset{a_u^k}{\int }}{p}_{u,m}^k\left(u\right)p_{u,n}^k\left(u\right)du---\left(23\right)$

${\left[G_{gu}^{uu,k}\right]}_{mn}=\underset{-a_u^k}{\overset{a_u^k}{\int }}{p}_{u,m}^k\left(u\right)\underset{-a_u^k}{\overset{a_u^k}{\int }}{p}_{u,n}^k\left(u_0\right)G_{gu}^k(u,u_0){\mathrm{du}}_{0}du---\left(24\right)$

采用2阶近似方法求解横向泄漏项的展开系数。假设在相邻3个节块内，横 向泄漏项关系式均成立，则可由3个节块平均横向泄漏获得3个方程，以确定展 开系数。系数表达式为：

$L_{gu1}^k={\bar{L}}_{gu}^k$

$\begin{array}{c}{L}_{gu2}^k=\frac{a_u^k}{2d_u}[(2a_u^{k-1}+a_u^k)·(a_u^{k-1}+a_u^k)({\bar{L}}_{gu}^{k+1}-{\bar{L}}_{gu}^k)\\ +(a_u^k+a_u^{k+1})·(a_u^k+2a_u^{k+1})({\bar{L}}_{gu}^k-{\bar{L}}_{gu}^{k-1})]\end{array}---\left(25\right)$

$L_{gu3}^k=\frac{1}{2d_u}{\left(a_u^k\right)}^2·[(a_u^{k-1}+a_u^k)({\bar{L}}_{gu}^{k+1}-{\bar{L}}_{gu}^k)+(a_u^{k+1}+a_u^k)({\bar{L}}_{gu}^{k-1}-{\bar{L}}_{gu}^k)]$

节块的平均泄漏可由净中子流求得：

$\begin{array}{c}{\bar{L}}_{gu}^k=\frac{1}{2a_u^k}{\int }_{-a_u^k}^{a_u^k}{L}_{gu}^k(u,t)du\\ =\frac{1}{2a_v^k}[J_{gv}^k(a_v^k,t)-J_{gv}^k(-a_v^k,t)]+\frac{1}{2a_w^k}[J_{gw}^k(a_w^k,t)-J_{gw}^k(-a_w^k,t)]\end{array}---\left(26\right)$

式中，表征平均横向泄漏，$d_u=(a_u^{k-1}+a_u^k)(a_u^{k+1}+a_u^k)(a_u^{k-1}+a_u^k+a_u^{k+1}),$上 划线代表节块体平均值。

基于第二类边界条件的节块格林函数瞬态求解的基本方程由式(11)至式(17) 组成的净流响应矩阵、中子平衡方程式(18)和偏通量展开系数方程式(21)构成， 可采用源迭代方法求解。

进一步的，所述执行堆芯三维瞬态物理计算的三维中子时空动力学模块采用 的时间离散方法为向后欧拉方法。其中，向后欧拉方法具备绝对-稳定性，并且 便于程序的串行耦合开发。

其中，采用向后欧拉方法的瞬态中子学计算流程，请参考图4，具体为：首 先转变上一时刻的相关量；然后准备截面；然后更新源项；然后依次求解各方向 界面平均净流；然后求解节块平均通量和偏中子通量系数；然后判断所有群能是 否满足要求，若不满足则返回执行更新源项，若满足，则判断裂变源是否收敛， 若不收敛，则返回执行更新源项，若收敛则更新先驱核浓度；然后进行下一时刻 计算。

下面结合附图和实施例对本发明作进一步详细的说明。

本发明提供了一种超临界水堆堆芯三维瞬态性能分析方法，应用本发明所述 分析方法的堆芯三维稳态-瞬态物理-热工水力耦合计算流程如图1所示，应用本 发明所述分析方法的稳态计算流程如图2所示，应用本发明所述分析方法的瞬态 计算流程如图3所示，图1、图2、图3仅仅给出了本发明提出概念的一种实施 例，下面结合图1、图2、图3和具体实施方式对本发明作进一步说明：

如图1所示，超临界水堆堆芯三维稳态-瞬态物理-热工水力耦合计算流程包 括(1)SCWR堆芯三维稳态物理-热工水力耦合计算、(2)稳态模块向瞬态模块 数据传递、(3)SCWR堆芯三维瞬态物理-热工水力耦合计算、(4)热组件子通道 热工-水力计算、(5)关键参数安全评价五个步骤，其中，SCWR堆芯三维稳态物 理-热工水力耦合计算(1)包括堆芯燃料管理模块计算稳态堆芯和子通道热工- 水力模块计算稳态堆芯两个部分，稳态模块向瞬态模块数据传递(2)传递的数 据包括初始堆芯状态和组件截面库，SCWR堆芯三维瞬态物理-热工水力耦合计算 (3)包括三维中子时空动力学模块计算瞬态堆芯和子通道热工-水力模块计算瞬 态堆芯两个部分，热组件子通道热工-水力计算(4)利用热工模块，对功率重构 后的热组件进行子通道计算分析，给出最大包壳壁面温度(MCST)，关键参数安 全评价(5)通过判断MCST是否超过相应的瞬态或事故安全限值，评价超临界水 堆堆芯在瞬态过程下的安全性能。

如图2所示，超临界水堆堆芯三维瞬态性能分析方法稳态计算流程中，三维 中子学模块和子通道热工模块分别进行堆芯稳态计算，并执行迭代耦合计算直至 功率收敛。

如图3所示，超临界水堆堆芯三维瞬态性能分析方法瞬态计算流程中，子通 道热工模块和三维中子时空动力学模块分别进行堆芯瞬态计算，并执行迭代耦合 计算直至功率收敛。

上述超临界水堆堆芯三维瞬态性能分析方法，建立起了堆芯三维稳态-瞬态 物理-热工水力耦合计算流程，能够准确描述物理与热工-水力之间的反馈，提供 精确的堆芯三维功率分布，实现真实地模拟超临界水堆的瞬态过程和事故过程的 目的。

下面举个具体的例子进行介绍：

应用本发明所述超临界水堆堆芯三维瞬态分析方法，分析了超临界水堆热态 零功率(HZP)和热态满功率(HFP)状态下各一个瞬态问题，以检验三维瞬态分 析方法的适用性。问题堆芯为超临界水堆CSR1000堆芯，HZP和HFP状态下堆芯 初始功率分别为2300W和2300MW。瞬态问题为弹棒问题，弹棒位置为第II流程 组件位置E11。HZP弹棒问题中，弹出棒初始位置距活性区底部为7.7cm，在0.1s 内弹至终止位置420cm。HFP弹棒问题中，弹出棒初始位置距活性区底部为210cm， 在0.1s内弹至终止位置420cm。HZP和HFP问题的弹出棒价值分别为683pcm和 293pcm。利用所述三维瞬态分析方法，首先执行CSR1000堆芯三维稳态物理-热 工水力耦合模块执行堆芯稳态计算，堆芯燃料管理模块向子通道热工-水力模块 提供功率参数，子通道热工-水力模块向堆芯燃料管理模块提供热工参数，执行 迭代耦合计算直至堆芯稳态功率参数和热工参数收敛；然后利用堆芯三维稳态物 理-热工水力耦合模块向堆芯三维瞬态物理-热工水力耦合模块提供堆芯初始状 态和组件截面库；然后执行CSR1000堆芯三维瞬态物理-热工水力耦合模块执行 瞬态计算，三维中子时空动力学模块向子通道热工-水力模块提供功率参数，子 通道热工-水力模块向三维中子时空动力学模块提供热工参数，执行迭代耦合计 算直至堆芯瞬态功率参数和热工参数收敛；然后在堆芯计算中对热组件进行功率 重构得到热组件栅元尺度的精细功率分布，再次利用热工模块进行热组件子通道 分析，最终给出安全评价关键参数-最大包壳壁面温度；最后判断最大包壳壁面 温度是否超过事故安全限值，评价CSR1000堆芯在弹棒瞬态过程下的安全性能。 三维瞬态分析结果表明，HZP弹棒问题的堆芯峰值功率为0.68NP(额定功率 2300MW)，最大包壳壁面温度为440℃，低于事故工况下安全限值1260℃；HFP 弹棒问题的堆芯峰值功率为1.74NP，最大包壳壁面温度为1075℃，低于事故工 况下安全限值1260℃。对CSR1000堆芯弹棒问题的三维瞬态分析表明，无论是 HFP弹棒还是HZP弹棒，事故过程中最大包壳壁面温度峰值均低于事故安全限值 1260℃，满足安全要求。

上述本申请实施例中的技术方案，至少具有如下的技术效果或优点：

由于采用了将超临界水堆堆芯三维瞬态性能分析方法设计为包括：步骤1： 超临界水堆堆芯三维稳态物理-热工水力耦合模块执行堆芯稳态计算，堆芯燃料 管理模块向子通道热工-水力模块提供功率参数，子通道热工-水力模块向堆芯燃 料管理模块提供热工参数，执行迭代耦合计算直至堆芯稳态功率参数和热工参数 收敛；步骤2：堆芯三维稳态物理-热工水力耦合模块向超临界水堆堆芯三维瞬 态物理-热工水力耦合模块提供堆芯初始状态和组件截面库；步骤3：超临界水 堆堆芯三维瞬态物理-热工水力耦合模块执行瞬态计算，三维中子时空动力学模 块向子通道热工-水力模块提供功率参数，子通道热工-水力模块向三维中子时空 动力学模块提供热工参数，执行迭代耦合计算直至堆芯瞬态功率参数和热工参数 收敛；步骤4：在堆芯计算中对热组件进行功率重构得到热组件栅元尺度的精细 功率分布，再次利用热工模块进行热组件子通道分析，最终给出安全评价关键参 数；步骤5：判断所述安全评价关键参数是否超过相应的瞬态或事故安全限值， 评价超临界水堆堆芯在瞬态过程下的安全性能的技术方案，即通过耦合三维中子 时空动力学与热工-水力学进行瞬态过程分析，中子动力学部分采用三维模型， 使其具备堆芯核热耦合三维功率计算能力，实现超临界水堆堆芯三维瞬态分析的 技术效果；所以，有效解决了现有的超临界水堆功率分布瞬态分析方法存在不具 备堆芯核热耦合三维功率计算能力，难以真实或精确地模拟例如反应性及功率分 布异常事件的瞬态过程的技术问题，实现了通过建立堆芯三维稳态-瞬态物理- 热工水力耦合计算流程，能够准确描述物理与热工-水力之间的反馈，提供较精 确的堆芯三维功率分布，真实地模拟超临界水堆的瞬态过程和事故过程的技术效 果。

进一步的，本方法可以真实地模拟超临界水堆堆芯瞬态过程，提高对超临界 水堆堆芯瞬态工况的认识，更加精确的估计安全裕量，改善堆芯设计并提高堆芯 性能。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本 创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意 欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明 的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等 同技术的范围之内，则本发明也意图包含这些改动和变型在内。