CN103902784A - 用于超临界水堆瞬态核热耦合的安全分析计算装置 - Google Patents

Abstract

本发明公开了一种用于超临界水堆瞬态核热耦合的安全分析计算装置，属于核电厂设计与反应堆安全分析领域，适用于超临界水堆，具体为，利用双群时空动力学方程的时空离散求解，建立瞬态物理分析程序；通过通道之间的流量分配计算和各个通道内部流动换热求解，建立瞬态热工分析程序；在此基础上，将物理计算程序嵌入至瞬态系统分析程序，建立超临界瞬态耦合分析程序，分析系统的安全性。特别是利用双群中子时空动力学方程建立的时空动力学准静态解法模型与多通道分析方法相结合，能在保证精度的基础上，又较子通道模型有效的减少了计算时间，进而提高超临界水堆设计的经济性，符合安全分析的发展趋势。

Description

用于超临界水堆瞬态核热耦合的安全分析计算装置

技术领域

本发明涉及核电厂设计与反应堆安全分析领域，具体的涉及用于超临界水堆瞬态核热耦合的安全分析方法和装置。

背景技术

超临界水堆是在第四代核能系统国际研讨会上，被GIF（第四代国际核能论坛）选定的作为长远开发目标的六种堆型之一，也是唯一被选定的轻水堆型。循环工质在吸收核裂变热之后，直接进入汽轮机做功，从而实现核能向热能、热能再向动力能的转变过程。为确保核能、流体动能安全地转换成动力能，堆芯内部能量转换与流动换热非常关键，而对于采用直流循环设计的超临界水堆则要求更高。

超临界水堆目前依然面临诸多难题，核热耦合作用就是其中之一。水在超临界条件下的物态表现与在亚临界条件下的物态表现完全不同。虽然水在超临界条件下是单相流体，但它的物性变化非常剧烈，而且这么大的物性变化往往会集中在一个较小的几何尺度内。极端的物性变化会引发局部慢化剂密度的不均匀性大大超过了常规压水堆工况，局部水铀比易受扰动，从而不可避免地引起其堆芯物理特性发生显著改变，导致堆芯局部反应性及功率的扰动，从而进一步引起堆芯流动工质温度及物性的急剧变化，形成异于常规压水堆强烈的反馈特性，核热耦合作用可能诱发堆芯稳定性问题。

所以，确保核反应堆工作时核热耦合的安全性十分重要，对于安全性的控制和分析也是迫切需要解决的问题，而安全控制的重点在于超临界水堆耦合计算程序。

迄今为止，国内外围绕超临界水堆耦合计算程序建立而进行了一些工作。

国外研究学者多采用程序开发手段。截止目前，已经开展了一定的超临界水堆物理-热工耦合工作，并给出可行性验证和应用分析。例如，2005年，Yamaji A.等开发了适用于高温超临界水堆SCLWR-H的三维核热耦合计算程序，并将其用于三维堆芯计算。它包含一维单通道热工分析程序SPROD和三维中子物理计算程序COREBN。2006年，Yoo J.等开发了适用于超临界水冷快堆SWFR的三维核热耦合计算程序，它包含三维扩散程序SRAC与单通道热工水力计算程序。将燃料组件的所有冷却剂通道平均近似，一个代表平均燃料棒通道，一个代表热通道。随后，Reiss T.和Monti L.等开发了超临界水堆稳态耦合程序，并建立了ERANOS/TRACE/CFD多尺度耦合分析程序，用于满功率稳态工况下HPLWR全堆芯物理与热工耦合计算。

国内研究学者更多地开展了很多系统程序改进方面的工作，并逐渐深入至核热耦合的研究。例如，2009年，胡珀，杨燕华等对堆芯计算程序PARCS和热工水力程序RELAP5进行了适应性改造，应用于美国超临界水堆系统分析。2010年，单建强等借助MCNP中子物理计算程序和ATHAS子通道热工分析程序，开发了稳态耦合程序，应用于加拿大超临界水堆系统分析。它采用对计算机或者服务器性能要求很高的外耦合方式，具备复杂集合体建模、连续能量精确求解、多维热工计算的能力。2009年，刘晓晶、程旭等基于子通道热工分析程序COBRA-IV和物理计算程序SKETCH-N建立了中子物理与热工水力耦合分析模型，应用于双排燃料组件分析。耦合界面利用数据传递实现。由上海交通大学牵头、程旭担任首席科学家的973计划“超临界水堆关键科学问题的基础研究项目”研究团队进行了大量的超临界水堆核热耦合程序开发工作。其中，清华大学王侃、张鹏等采用蒙特卡洛方法，进行了超临界水堆三维物理计算以及物理-热工耦合研究；上海核工程研究设计院廖成奎、上海交通大学张少泓等结合超临界水热质传输特性研究，进行了超临界水堆核热耦合行为研究。2010年，中国核动力院安萍、姚栋等开发了核热耦合程序MCATHAS，应用于欧洲超临界水堆燃料组件分析；2009年，上海交通大学刘占权等提出了美国超临界水堆的一维稳态耦合模型；2011年，西安交通大学杨萍、曹良志等建立了超临界水堆三维中子物理热工耦合模型，应用于加拿大超临界坎杜堆的堆芯设计。华北电力大学周涛团队成员刘晓壮、李臻洋、孙灿辉、陈娟、罗峰、王晗丁也建立了超临界水堆热工分析程序TH-02、组件中子物理计算程序DRAD，并进行了稳态核热耦合分析。

上述各种方法中，蒙特卡罗计算模型具有强大的几何处理功能及连续能量中子截面库的应用，但影响其广泛应用的最主要原因是其计算耗时长，效率较低。而子通道模型虽然计算精确，但其同样存在计算耗时长，程序过于复杂的问题。单通道的计算模型虽说程序简单，计算速度较高，但是计算精度有所欠缺，使得用于超临界水堆瞬态核热耦合的安全分析方法和装置仍有所欠缺。

发明内容

为了克服上述问题，本发明人对现有的用于超临界水堆瞬态核热耦合的安全分析方法和计算装置进行了锐意研究，结果发现：通过设置一个安全控制装置，该装置包括瞬态物理模块、瞬态热工模块、瞬态耦合模块和瞬态分析模块，通过瞬态物理模块和瞬态热工模块计算各个时刻的相关参数，通过瞬态耦合模块将瞬态物理模块和瞬态热工模连接起来，将瞬态物理模块和瞬态热工模产生的数据交互传递，并且在传递过程中进行数据转换，便于瞬态物理模块和瞬态热工模的计算，更新计算的初始参数，同时，由瞬态分析模块接收瞬态物理模块和计算瞬态热工模计算出的数据，并用来系统的安全分析，以保证整个反应堆的安全，从而完成本发明。

本发明的目的在于提供以下方面：

（1）一种用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，该装置包括：瞬态物理模块、瞬态热工模块、瞬态耦合模块和瞬态分析模块；

瞬态物理模块，在启动后的第一个大时间步长内，所述瞬态物理模块接收从稳态物理装置传递来的信息，所述信息包括中子通量分布信息、中子密度信息和中子价值信息，根据接收到的信息计算得到轴向功率分布信息，并传出所述轴向功率分布信息和计算过程中得到的中间过程信息，所述中间过程信息包括中子通量分布信息、先驱核浓度信息和功率值信息；在第一个大时间步长结束后，在第二个或其后的大时间步长内，所述瞬态物理模块接收从瞬态耦合模块传递来的信息，所述信息包括中子通量分布信息、中子密度信息和中子价值信息，根据接收到的信息计算得到轴向功率分布信息，并传出所述轴向功率分布信息和计算过程中得到的中间过程信息，所述中间过程信息包括中子通量分布信息、先驱核浓度信息和功率值信息，进入下一个大时间步长时重复该过程；

瞬态热工模块，在启动后的第一个大时间步长内，所述瞬态热工模块接收稳态热工参数信息，所述稳态热工参数信息包括稳态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，根据接收到的稳态热工参数信息计算得到工质温度分布信息和工质密度分布信息，并传出计算得到工质温度分布信息和工质密度分布信息以及计算过程中得到的中间过程信息，所述中间过程信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度；在第一个大时间步长结束后，在第二个或其后的大时间步长内，所述瞬态热工模块接收从瞬态耦合模块传递来的瞬态热工参数信息，所述瞬态热工参数信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，根据接收到的瞬态热工参数信息计算得到工质温度分布信息和工质密度分布信息，并传出计算得到工质温度分布信息和工质密度分布信息以及计算过程中得到的中间过程信息，所述中间过程信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，进入下一个大时间步长时重复该过程；

瞬态耦合模块，其用于接收每个大时间步长内，瞬态物理模块和瞬态热工模块传出信息，对接收到的信息进行转换，并将转换所得的数据信息分别传给瞬态热工模块和瞬态物理模块；

瞬态分析模块，其用于接收瞬态物理模块在每个大时间步长输出的轴向功率分布信息和瞬态热工模块在每个大时间步长输出的工质温度分布信息和工质密度分布信息，并将接收到的信息与安全系数信息进行比较，若接收到的信息超出允许的安全范围，发出指令控制降低反应堆反应速率，并报警；若接收到的信息在允许的安全范围内，待接收到下一个大时间步长的信息后继续与安全系数信息进行比较。

（2）根据上述（1）所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，在瞬态物理模块中，轴向功率分布通过包括如下步骤的过程计算得到：

第一步，计算出小时间步长下不同时刻的堆芯裂变总功率，包括以下子步骤：

a.设定中子通量密度仅沿堆芯轴向变化，建立如下双群中子时空动力学方程(I)和(II)：

$\begin{array}{c}\frac{1}{V_g}\frac{\∂{\mathrm{\Φ}}_g(z,t)}{\∂t}=\\ \▿D_g{\▿\mathrm{\Φ}}_g(z,t)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Φ}}_g(z,t)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}{+\mathrm{\Σ}}_{g^{\′}\→g}\right]{\mathrm{\Φ}}_{g^{\′}}(z,t)+x_g\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)\end{array}---\left(I\right)$

$\frac{{\∂C}_I(z,t)}{\∂t}=\frac{1}{k_{\mathrm{eff}}}{\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Φ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(\mathrm{II}\right)$

其中g＝1,2，I＝1,2,...6；

b.把中子通量密度用与时间相关的幅函数以及与空间相关的形状函数的乘积表示，其表达式如下式(III)：

Φ_，(z，t)＝n(t)Ψ_g(z，)  (III)

进而得到如下双群中子时空动力学时空离散方程(IV)和(V)：

$\begin{array}{c}\frac{1}{V_g}\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}\mathrm{\Ψ}(z,t)+\frac{1}{V_g}\frac{d{\mathrm{\Ψ}}_g(z,t)}{\mathrm{dt}}n\left(t\right)=\\ n\left(t\right)\▿D_g\▿{\mathrm{\Ψ}}_g(z,t)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,t)n\left(t\right)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{\mathrm{\χ}}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\Ψ}}_{g^{\′}}(z,t)n\left(t\right)+x_g\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)\end{array}---\left(\mathrm{IV}\right)$

$\frac{d{C}_I(z,t)}{\mathrm{dt}}=\frac{1}{k_{\mathrm{eff}}}{n\left(t\right)\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(V\right)$

其中，g＝1,2I＝1,2,...6；

c.将双群中子时空动力学方程(I)和(II)两边同时乘以中子价值，并且沿轴向进行积分，合并化简，得出如下点堆中子动力学方程(VI)和(VII)：

$\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}=\frac{\mathrm{\ρ}\left(t\right)-\stackrel{\‾}{\mathrm{\β}}\left(t\right)}{\mathrm{\Λ}\left(t\right)}n\left(t\right)+\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I\left(t\right)---\left(\mathrm{VI}\right)$

$\frac{d{\stackrel{\‾}{C}}_I\left(t\right)}{\mathrm{dt}}=-{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I\left(t\right)+\frac{{\stackrel{\‾}{\mathrm{\β}}}_I\left(t\right)}{\mathrm{\Λ}\left(t\right)}n\left(t\right)---\left(\mathrm{VII}\right),$

其中，I=1,2,…6

式中n(t)表示中子密度(中子/cm³)，

表示总的缓发中有效份额，

表示总的先驱核数(原子/cm³)，

Λ表示瞬发中子每代时间(s)，

设定中子通量总价值保持不变，向后差分，对点堆中子动力学方程(VI)和(VII)进行时间离散，得到如下方程(VIII)和(IX)：

$\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}}=\frac{\mathrm{\ρ}(j+1)-\stackrel{\‾}{\mathrm{\β}}(j+1)}{\mathrm{\Λ}(j+1)}n(j+1)+\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I(j+1)---\left(\mathrm{VIII}\right)$

$\frac{C_I(j+1)-C_I\left(j\right)}{\mathrm{\Δt}}=-{\mathrm{\λ}}_I{C}_I(j+1)+\frac{{\mathrm{\β}}_I(j+1)}{\mathrm{\Λ}(j+1)}N(j+1)---\left(\mathrm{IX}\right)$

将方程(VIII)和(IX)化简，得到如下的功率计算式(X)和缓发中子先驱核浓度(XI)：

$n(j+1)=\frac{1}{1-\frac{\mathrm{\ρ}(j+1)-\mathrm{\β}}{\mathrm{\Λ}}\mathrm{\Δt}}n\left(t\right)+\frac{\mathrm{\Δt}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(j+1)}{1-\frac{\mathrm{\ρ}(j+1)-\mathrm{\β}}{\mathrm{\Λ}}\mathrm{\Δt}}---\left(X\right)$

$C_I(j+1)=\frac{C_I\left(j\right)}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}+\frac{{\mathrm{\β}}_I(j+1)\mathrm{\Δt}}{\mathrm{\Λ}(j+1)(1+{\mathrm{\λ}}_I\mathrm{\Δt})}N(j+1)---\left(\mathrm{XI}\right)$

进而计算出小时间步长下不同时刻的堆芯裂变总功率；

第二步，使用100%稳态工况时的径向功率分布因子进行稳态计算；

第三步，计算得出堆芯功率轴向分布，包括以下子步骤：

a.化简第一步中的公式(IV)和公式(V)，得到的形状函数时空离散方程(XII)和缓发中子先驱核浓度时空离散方程(XIII)，组成形状函数方程组：

$\frac{1}{V_g}\frac{d{\mathrm{\Ψ}}_g(z,t)}{\mathrm{dt}}=$

$\begin{array}{c}{\▿D}_g{\▿D}_g(z,t){\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,t)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\ψ}}_{g^{\′}}(z,t)+\frac{x_g}{n\left(t\right)}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)\\ -\frac{1}{V_{g}n\left(t\right)}\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}\mathrm{\Ψ}(z,t)\end{array}---\left(\mathrm{XII}\right)$

其中，g＝1,2

$\frac{{\∂C}_I(z,t)}{\∂t}=\frac{1}{k_{\mathrm{eff}}}{n\left(t\right)\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(\mathrm{XIII}\right)$

其中，I＝1,2,...6

b.向后差分，对形状函数方程组进行离散求解，得到如下形状函数时间离散方程(XIV)和(XV)：

$\begin{array}{c}\frac{1}{V_g}\frac{{\mathrm{\Ψ}}_g(z,j+1)-{\mathrm{\Ψ}}_g(z,j)}{\mathrm{\Δt}}=\\ \frac{d}{\mathrm{dz}}{D}_g\frac{d}{\mathrm{dz}}{\mathrm{\Ψ}}_g(z,j+1)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,j+1)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\Ψ}}_{g^{\′}}(z,j+1)---\left(\mathrm{XIV}\right)\\ +\frac{x_g}{n(j+1)}\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,j+1)-\frac{1}{V_{g}n(j+1)}\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}}\mathrm{\Ψ}(z,j+1)\end{array}$

$\frac{C_I(z,j+1)-C_I(z,j)}{\mathrm{\Δt}}=\frac{1}{k_{\mathrm{eff}}}n(j+1){\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left[v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,j+1)\right]-{\mathrm{\λ}}_I{C}_I(z,j+1)---\left(\mathrm{XV}\right)$

c.将形状函数时间离散方程组的第二个方程(XV)进行化简，得到缓发中子先驱核浓度的离散求解表达式，然后将其代入形状函数时间离散方程组的第一个方程(XIV)进行化简，得到如下各个轴向节点i处的形状函数求解表达式(XVI)：

${\mathrm{\Ψ}}^i(j+1)=\frac{[S_{\text{g}}^{\′\′i-1/2}(j+1)+S_g^{\′\′i+1/2}(j+1)]\·\frac{{{\mathrm{\Δ}}_z}^2}{2}+D^{i-1/2}{\mathrm{\Ψ}}^{i-1}(j+1)+D^{i+1/2}{\mathrm{\ψ}}^{i+1}(j+1)}{D^{i+1/2}+D^{i-1/2}+({\mathrm{\Σ}}_r^{i-1/2}+{\mathrm{\Σ}}_r^{i+1/2})\·\frac{{{\mathrm{\Δ}}_z}^2}{2}+(\frac{1}{V_{g}n(j+1)}\·\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}})\·{{\mathrm{\Δ}}_z}^2}---\left(\mathrm{XVI}\right)$

其中， $S_g^{\′\′}(j+1)=\left[\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left(v{\mathrm{\Σ}}_{{\mathrm{fg}}^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(j+1)\right)\right]\·[(1-\mathrm{\β})x_g+x_g\underset{I=1}{\overset{6}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_I{\mathrm{\β}}_I\mathrm{\Δt}}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}]\·\frac{1}{k_{\mathrm{eff}}}$

$+\frac{1}{V_g\mathrm{\Δt}}{\mathrm{\Ψ}}_g\left(j\right)+\frac{1}{n(j+1)}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_I{x}_g{C}_I(z,j)}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{g^{\′}\→}{\mathrm{\Ψ}}_{g^{\′}}(j+1)---\left(\mathrm{XVII}\right)$

上述公式XVII为中子产生项，可以用来计算产生中子的数量，

对得到形状函数进行归一化判断，判断准则为中子总价值的计算值与基准值之间的相对误差小于10^-5；每一个中等时间步长进行一次归一化判断，若形状函数符合归一化标准，则用该形状方程的解去求解中子通量和缓发中子先驱核浓度，然后转入下一个时刻的计算；若形状方程的解不满足归一化条件，则进行“形状函数迭代更新”，并重复求解和判定过程；

d.根据由点堆中子动力学方程求解得到的幅值函数，以及形状函数，得到如下用于计算不同时间步长下轴向功率分布的公式(XVIII)：

q_l ⁱ(j+1)＝n(j+1)·Ψⁱ(j+1)·ξ₂  (XVIII)

其中，q_l ⁱ(j+1)表示在j+1时刻轴向节点i的燃料棒线功率密度(W/m)；

ξ₂表示时空离散的功率修正因子。

（3）根据上述（1）所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，在瞬态热工模块中，工质温度分布和工质密度分布通过包括如下步骤的过程计算得到：

第一步，建立瞬态流量分配计算模型，进行流量分配；

第二部，计算各个通道的总压降，以面积为权重求解所有通道的平均压降，其计算公式如下述公式(XIX)：

$\mathrm{\Δ}{p}_{\mathrm{av}}=\frac{\underset{i}{\mathrm{\Σ}}\mathrm{\Δ}{p}_{i,t}{A}_i}{\underset{i}{\mathrm{\Σ}}{A}_i}---\left(\mathrm{XIX}\right)$

第三步，判断平均压降的收敛性，若满足收敛条件，即

时，完成流量分配，其中，收敛准则ε取值为10^-5，Δp_i是指各个通道的总压降；

第四步，计算不同时刻燃料芯块温度及包壳温度，其计算公式如下述公式(XX)：

$\frac{\∂T}{\∂t}=\frac{\mathrm{\ΔQ}}{\mathrm{\ρ}\·c_p}---\left(\mathrm{XX}\right)$

其中，ρ表示密度(kg/m³)，

c_p表示定压比热容(J/(kg·K))，

ΔQ表示燃料芯块内热源释热量与外表面热通量的差值，或者包壳内、外表面的热通量差值(W)；

第五步，计算冷却剂与慢化剂的温度分布，其计算公式如下述公式(XXI)：

Q＝(T₁-T₂)/ΣR  (XXI)

其中，Q表示燃料棒外表面或水棒内、外壁面的释热功率(W)；

T₁表示传热过程起始点的温度(℃)；

T₂表示传热过程终止点的温度(℃)；

ΣR表示换热热阻(W/℃)，

第六步，求解工质温度分布，

冷却剂通道内的工质流动满足如下三大守恒方程：

质量守恒方程(XXII)：

$\frac{\∂w_t}{\∂t}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\∂w_i}{\∂t}---\left(\mathrm{XXII}\right)$

动量守恒方程(XXIII)：

$\frac{1}{A_i}\frac{\∂w_i}{\∂t}+\frac{\∂}{\∂z}\left(\frac{{w_i}^2}{{{\mathrm{\ρA}}_i}^2}\right)=-\frac{\∂P}{\∂z}\±\mathrm{\ρ}\stackrel{\‾}{g}-\frac{f{w}_i\left|w_i\right|}{2d_e{{\mathrm{\ρA}}_i}^2}---\left(\mathrm{XXIII}\right)$

能量守恒方程(XXIV)：

$\mathrm{\ρ}\frac{\∂h}{\∂t}+\frac{\∂}{\∂z}\left(\frac{w_{i}h}{\mathrm{\ρ}{A}_i}\right)=q_v^{\mathrm{channel}}\left(t\right)---\left(\mathrm{XXIV}\right)$

其中，z表示轴向节点高度(m)；

t表示时间(s)；

ρ表示冷却剂密度(kg/m³)；

f表示摩擦阻力系数；

h表示冷却剂焓值(J/kg)；

w_t表示冷却剂总流量(kg/s)；

w_i表示冷却剂通道i的流量(kg/s)；

A_i表示冷却剂通道i的流通面积(m²)；

d_e表示水力当量直径(m)；

表示t时刻单位体积工质所吸收的热量(W/m³)。

（4）根据上述（3）所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，建立瞬态流量分配计算模型的过程包括如下步骤：

将动量守恒方程(XXIII)沿通道长度进行积分，得到各冷却剂通道内部的通道流量瞬态变化方程组(XXVI)：

$\frac{L}{A_1}\frac{{\∂w}_1}{\∂t}=\mathrm{\Δ}{p}_{t,1}\mathrm{\Δ}{p}_{a,1}-\mathrm{\Δ}{p}_{f,1}-\mathrm{\Δ}{p}_{e1,1}$

$\frac{L}{A_2}\frac{{\∂w}_2}{\∂t}={\mathrm{\Δp}}_{t,2}-\mathrm{\Δ}{p}_{a,2}-{\mathrm{\Δp}}_{f,2}-{\mathrm{\Δp}}_{e1,2}$

……

$\frac{L}{A_{14}}\frac{\∂w_{14}}{\∂t}={\mathrm{\Δp}}_{t,14}-{\mathrm{\Δp}}_{a,14}-{\mathrm{\Δp}}_{f,14}-{\mathrm{\Δp}}_{\mathrm{el},14}---\left(\mathrm{XXVI}\right)$

其中，L表示通道长度(m)，

Δp_t表示并联通道的总压降(Pa)，

Δp_a表示加速压降(Pa)，

Δp_f表示摩擦压降(Pa)，

Δp_el表示提升压降(Pa)，

将质量守恒方程(XXII)和通道流量瞬态变化方程组(XXVI)联立，得到瞬态流量分配计算模型。

（5）根据上述（1）所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，瞬态核热耦合模块包括中子截面数据库，

中子截面数据库接收瞬态物理模块计算得到的中子通量分布信息，根据中子通量分布信息计算出轴向功率分布信息，并将轴向功率分布信息输入至瞬态热工模块；

中子截面数据库接收瞬态热工模块计算得到的轴向冷却剂温度分布信息、慢化剂温度分布信息，根据轴向冷却剂温度分布信息、慢化剂温度分布信息计算出中子截面分布信息，并将中子截面分布信息输入至瞬态物理模块；

中子截面数据库接收瞬态物理模块计算得到的幅值函数，根据幅值函数计算出总功率，并将总功率输入至瞬态热工模块。

瞬态核热耦合模块包括中子截面数据库，中子截面数据库接收瞬态物理模块传递的截面数值，根据截面数值计算出轴向功率分布，并向瞬态热工模块传递轴向功率分布；

中子截面数据库接收瞬态热工模块传递的轴向工质温度信息和工质密度信息，并将工质温度信息和工质密度信息传递到瞬态物理模块。

（6）根据上述（1）所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，安全系数信息是指燃料包壳允许达到的最大温度、燃料芯块中心最高温度和最高系统压力，其中，燃料包壳允许达到的最大温度是1260℃，燃料芯块中心允许达到的最高温度是2590℃，系统压力允许达到的最大压力值是30.3MPa。

（7）一种用于控制超临界水堆瞬态核热耦合安全的方法，其特征在于，该方法通过用于超临界水堆瞬态核热耦合的安全分析计算装置实现。

根据本发明提供的用于超临界水堆瞬态核热耦合的安全分析计算装置，中的瞬态物理模块加入反应随空间的变化量，使该装置计算结果更能反映实际的参数分布，与此同时设置瞬态耦合模块，使得该装置能在保证精度的基础上，又较子通道模型有效的减少了程序的运行时间，使该装置具有更好的可靠性和经济性，该装置具体具有如下有益效果：

第一、该装置的瞬态物理模块，采用双群中子时空动力学方程离散求解，在计算中子通量密度时引入了空间变量，使其成为时间和空间的函数，使得计算过程更符合实际情况，同时也提高了反应堆的经济性；

第二、该装置设置中子截面数据库，增加了该装置计算的可靠性；

第三、该装置在瞬态热工模块中加入流量分配程序，使得该装置能在保证精度的基础上，又较子通道模型有效的减少了计算时间；

附图说明

图1示出根据本发明一种优选实施方式的耦合流程图；

图2示出根据本发明一种优选实施方式的瞬态物理模块计算流程图；

图3示出根据本发明一种优选实施方式的时间步长示意图。

附图标号说明：

1-大时间步长

2-中等时间步长

3-小时间步长

具体实施方式

下面通过结合附图对本发明进行详细说明，本发明的特点和优点将随着这些说明而变得更为清楚、明确。

在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。

在本文中出现的：“功率分布”是指反应堆燃料棒内部不同节点的热功率的分布情况；

“轴向功率因子”是指的是沿反应堆燃料棒轴向方向的功率分布情况；

本文中出现的各个符号的意义和单位参见下表

在本发明一种优选的实施方式中，提供用于超临界水堆瞬态核热耦合的安全分析计算装置，该装置包括：瞬态物理模块、瞬态热工模块、瞬态耦合模块和瞬态分析模块；

在一个优选的实施方式中，瞬态物理模块在启动后的第一个大时间步长内，接收从稳态物理装置传递来的信息，所述信息包括中子通量分布信息、中子密度信息和中子价值信息，根据接收到的信息计算得到轴向功率分布信息，并传出所述轴向功率分布信息和计算过程中得到的中间过程信息，所述中间过程信息包括中子通量分布信息、先驱核浓度信息和功率值信息；在第一个大时间步长结束后，在第二个或其后的大时间步长内，所述瞬态物理模块接收从瞬态耦合模块传递来的信息，所述信息包括中子通量分布信息、中子密度信息和中子价值信息，根据接收到的信息计算得到轴向功率分布信息，并传出所述轴向功率分布信息和计算过程中得到的中间过程信息，所述中间过程信息包括中子通量分布信息、先驱核浓度信息和功率值信息，进入下一个大时间步长时重复该过程；其中，瞬态物理模块中的计算过程为：采用因式分解方法，对时空动力学方程进行时空离散求解，将时空动力学方程求解转换为幅值函数与形状函数的求解，如图2中所示，具体步骤如下：

（1）利用稳态物理计算，获得初始时刻的参数分布，包括：初始中子通量、初始中子密度和中子价值等，并用来计算形状函数分布；

（2）利用接收到的当前时刻中子通量分布，进行截面计算，更新形状函数分布，每一个中等时间步长更新一次；

（3）利用当前时刻中子截面数值，进行各个小时间步长的点堆中子动力学方程求解；

（4）利用形状函数的归一化判断，若满足归一化条件，得到当前大时间步长的形状函数分布，若不满足，从新进行形状函数迭代更新；

（5）更新中子通量分布、先驱核浓度及功率，进入下一时刻循环，直至完成整个大时间步长内的计算。

在更进一步优先选的实施方式中，瞬态物理模块计算轴向功率分布是通过包括如下步骤的过程计算得到：

第一步，计算出小时间步长下不同时刻的堆芯裂变总功率，包括以下子步骤：

a.设定中子通量密度仅沿堆芯轴向变化，建立如下双群中子时空动力学方程(I)和(II)：

$\begin{array}{c}\frac{1}{V_g}\frac{\∂{\mathrm{\Φ}}_g(z,t)}{\∂t}=\\ \▿D_g{\▿\mathrm{\Φ}}_g(z,t)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Φ}}_g(z,t)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}{+\mathrm{\Σ}}_{g^{\′}\→g}\right]{\mathrm{\Φ}}_{g^{\′}}(z,t)+x_g\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)\end{array}---\left(I\right)$

$\frac{{\∂C}_I(z,t)}{\∂t}=\frac{1}{k_{\mathrm{eff}}}{\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Φ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(\mathrm{II}\right)$

其中g＝1,2，I＝1,2,...6；

b.把中子通量密度用与时间相关的幅函数以及与空间相关的形状函数的乘积表示，其表达式如下式(III)：

Φ_g(z,t)＝n(t)Ψ_g(z,t)  (III)

进而得到如下双群中子时空动力学时空离散方程(IV)和(V)：

$\begin{array}{c}\frac{1}{V_g}\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}\mathrm{\Ψ}(z,t)+\frac{1}{V_g}\frac{d{\mathrm{\Ψ}}_g(z,t)}{\mathrm{dt}}n\left(t\right)=\\ n\left(t\right)\▿D_g\▿{\mathrm{\Ψ}}_g(z,t)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,t)n\left(t\right)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{\mathrm{\χ}}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\Ψ}}_{g^{\′}}(z,t)n\left(t\right)+x_g\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)\end{array}---\left(\mathrm{IV}\right)$

$\frac{d{C}_I(z,t)}{\mathrm{dt}}=\frac{1}{k_{\mathrm{eff}}}{n\left(t\right)\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(V\right),$

其中，9，1，2f，1，2，...6；

c.中子价值是表征反应堆内中子群体特性的物理量，中子在反应堆内位置不同，对反应堆的链式反应影响也是不同的，中子价值正是表征这种差异的物理量，根据初始时刻的中子价值定义，将双群中子时空动力学方程(I)和(II)两边同时乘以中子价值，并且沿轴向进行积分，合并化简，得出如下点堆中子动力学方程(VI)和(VII)：

$\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}=\frac{\mathrm{\ρ}\left(t\right)-\stackrel{\‾}{\mathrm{\β}}\left(t\right)}{\mathrm{\Λ}\left(t\right)}n\left(t\right)+\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I\left(t\right)---\left(\mathrm{VI}\right)$

$\frac{d{\stackrel{\‾}{C}}_I\left(t\right)}{\mathrm{dt}}=-{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I\left(t\right)+\frac{{\stackrel{\‾}{\mathrm{\β}}}_I\left(t\right)}{\mathrm{\Λ}\left(t\right)}n\left(t\right)---\left(\mathrm{VII}\right),$

其中，I=1,2,…6

式中n(t)表示中子密度(中子/cm³)，

表示总的缓发中有效份额，

表示总的先驱核数(原子/cm³)，

Λ表示瞬发中子每代时间(s)，

设定中子通量总价值保持不变，向后差分，对点堆中子动力学方程(VI)和(VII)进行时间离散，得到如下方程(VIII)和(IX)：

$\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}}=\frac{\mathrm{\ρ}(j+1)-\stackrel{\‾}{\mathrm{\β}}(j+1)}{\mathrm{\Λ}(j+1)}n(j+1)+\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I(j+1)---\left(\mathrm{VIII}\right)$

$\frac{C_I(j+1)-C_I\left(j\right)}{\mathrm{\Δt}}=-{\mathrm{\λ}}_I{C}_I(j+1)+\frac{{\mathrm{\β}}_I(j+1)}{\mathrm{\Λ}(j+1)}N(j+1)---\left(\mathrm{IX}\right)$

将方程(VIII)和(IX)化简，得到如下的功率计算式(X)和缓发中子先驱核浓度(XI)：

$n(j+1)=\frac{1}{1-\frac{\mathrm{\ρ}(j+1)-\mathrm{\β}}{\mathrm{\Λ}}\mathrm{\Δt}}n\left(t\right)+\frac{\mathrm{\Δt}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(j+1)}{1-\frac{\mathrm{\ρ}(j+1)-\mathrm{\β}}{\mathrm{\Λ}}\mathrm{\Δt}}---\left(X\right)$

$C_I(j+1)=\frac{C_I\left(j\right)}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}+\frac{{\mathrm{\β}}_I(j+1)\mathrm{\Δt}}{\mathrm{\Λ}(j+1)(1+{\mathrm{\λ}}_I\mathrm{\Δt})}N(j+1)---\left(\mathrm{XI}\right)$

进而计算出小时间步长下不同时刻的堆芯裂变总功率；

第二步，使用100%稳态工况时的径向功率分布因子进行稳态计算；

第三步，计算得出堆芯功率轴向分布，包括以下子步骤：

a.化简第一步中的公式(IV)和公式(V)，得到的形状函数时空离散方程(XII)和缓发中子先驱核浓度时空离散方程(XIII)，组成形状函数方程组：

$\begin{array}{c}\frac{1}{V_g}\frac{{\mathrm{d\Ψ}}_g(z,t)}{\mathrm{dt}}=\\ \▿D_g{\▿\mathrm{\Ψ}}_g(z,t)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,t)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\Ψ}}_{g^{\′}}(z,t)---\\ +\frac{x_g}{n\left(t\right)}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)-\frac{1}{V_{g}n\left(t\right)}\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}\mathrm{\ψ}(z,t)\end{array}\left(\mathrm{XII}\right)$

其中，g＝1,2

$\frac{{\∂C}_I(z,t)}{\∂t}=\frac{1}{k_{\mathrm{eff}}}{n\left(t\right)\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(\mathrm{XIII}\right)$

其中，I＝1,2,...6

b.向后差分，对形状函数方程组进行离散求解，得到如下形状函数时间离散方程(XIV)和(XV)：

$\begin{array}{c}\frac{1}{V_g}\frac{{\mathrm{\Ψ}}_g(z,j+1)-{\mathrm{\Ψ}}_g(z,j)}{\mathrm{\Δt}}=\\ \frac{d}{\mathrm{dz}}{D}_g\frac{d}{\mathrm{dz}}{\mathrm{\Ψ}}_g(z,j+1)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,j+1)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\Ψ}}_{g^{\′}}(z,j+1)---\left(\mathrm{XIV}\right)\\ +\frac{x_g}{n(j+1)}\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,j+1)-\frac{1}{V_{g}n(j+1)}\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}}\mathrm{\Ψ}(z,j+1)\end{array}$

$\frac{C_I(z,j+1)-C_I(z,j)}{\mathrm{\Δt}}=\frac{1}{k_{\mathrm{eff}}}n(j+1){\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left[v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,j+1)\right]-{\mathrm{\λ}}_I{C}_I(z,j+1)---\left(\mathrm{XV}\right)$

c.将形状函数时间离散方程组的第二个方程(XV)进行化简，得到缓发中子先驱核浓度的离散求解表达式，然后将其代入形状函数时间离散方程组的第一个方程(XIV)进行化简，得到如下各个轴向节点i处的形状函数求解表达式(XVI)：

${\mathrm{\Ψ}}^i(j+1)=\frac{[S_{\text{g}}^{\′\′i-1/2}(j+1)+S_g^{\′\′i+1/2}(j+1)]\·\frac{{{\mathrm{\Δ}}_z}^2}{2}+D^{i-1/2}{\mathrm{\Ψ}}^{i-1}(j+1)+D^{i+1/2}{\mathrm{\ψ}}^{i+1}(j+1)}{D^{i+1/2}+D^{i-1/2}+({\mathrm{\Σ}}_r^{i-1/2}+{\mathrm{\Σ}}_r^{i+1/2})\·\frac{{{\mathrm{\Δ}}_z}^2}{2}+(\frac{1}{V_{g}n(j+1)}\·\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}})\·{{\mathrm{\Δ}}_z}^2}---\left(\mathrm{XVI}\right)$

其中， $S_g^{\′\′}(j+1)=\left[\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left(v{\mathrm{\Σ}}_{{\mathrm{fg}}^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(j+1)\right)\right]\·[(1-\mathrm{\β})x_g+x_g\underset{I=1}{\overset{6}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_I{\mathrm{\β}}_I\mathrm{\Δt}}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}]\·\frac{1}{k_{\mathrm{eff}}}$

$+\frac{1}{V_g\mathrm{\Δt}}{\mathrm{\Ψ}}_g\left(j\right)+\frac{1}{n(j+1)}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_I{x}_g{C}_I(z,j)}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{g^{\′}\→}{\mathrm{\Ψ}}_{g^{\′}}(j+1)---\left(\mathrm{XVII}\right)$

上述公式XVII为中子产生项，可以用来计算产生中子的数量，

对得到形状函数进行归一化判断，判断准则为中子总价值的计算值与基准值之间的相对误差小于10^-5；每一个中等时间步长进行一次归一化判断，若形状函数符合归一化标准，则用该形状方程的解去求解中子通量和缓发中子先驱核浓度，然后转入下一个时刻的计算；若形状方程的解不满足归一化条件，则进行“形状函数迭代更新”，并重复求解和判定过程；

d.根据由点堆中子动力学方程求解得到的幅值函数，以及形状函数，得到如下用于计算不同时间步长下轴向功率分布的公式(XVIII)：

q_l ⁱ(j+1)＝n(j+1)·Ψⁱ(j+1)·ξ₂  (XVIII)

其中，q_l ⁱ(j+1)表示在j+1时刻轴向节点i的燃料棒线功率密度(W/m)；

ξ₂表示时空离散的功率修正因子，

进而计算出轴向功率分布。

在一个优选的实施方式中，瞬态热工模块在启动后的第一个大时间步长内接收稳态热工参数信息，所述稳态热工参数信息包括稳态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，根据接收到的稳态热工参数信息计算得到工质温度分布信息和工质密度分布信息，并传出计算得到工质温度分布信息和工质密度分布信息以及计算过程中得到的中间过程信息，所述中间过程信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度；在第一个大时间步长结束后，在第二个或其后的大时间步长内，所述瞬态热工模块接收从瞬态耦合模块传递来的瞬态热工参数信息，所述瞬态热工参数信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，根据接收到的瞬态热工参数信息计算得到工质温度分布信息和工质密度分布信息，并传出计算得到工质温度分布信息和工质密度分布信息以及计算过程中得到的中间过程信息，所述中间过程信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，进入下一个大时间步长时重复该过程；瞬态物理模块和瞬态热工模块同时开始，在一个大时间步长结束时同时完成计算和传输过程。

在一个优选的实施方式中，借鉴稳态热工水力程序，编制瞬态热工水力计算程序，其具体步骤包括：

（1）利用稳态工况热工参数分布信息和轴向功率分布信息，编制输入文件，将输入文件输入至瞬态热工模块；

（2）定义瞬态及其时间步长，确定相关变量随时间的变化情况，其中，相关变量为热工计算需要用到的所有变量，包括：燃料芯块温度、包壳表面温度、冷却剂和慢化剂的温度和密度等；

（3）利用瞬态流量分配模型，编制通道之间的流量分配计算模块；具体为通过瞬态流量分配模型传递的各个通道入口流量数据生成流量分配计算模块；

（4）通过嵌入多通道热工水力计算模块，完成当前时刻的热工平衡调整；

（5）更新堆芯内部的热工参数分布，进入下一时刻的热工水力计算；每经过一个大时间步长进行一次堆芯内部的热工参数分布更新，更新该参数是通过瞬态耦合模块传递的数据；

（6）重复上述步骤（2）-（5），直至完成所有时间步长下的热工水力计。

在进一步的优选实施方式中，瞬态热工模块的工质温度分布和工质密度分布是通过包括如下步骤的过程计算得到：

第一步，建立瞬态流量分配计算模型，进行流量分配；

第二部，计算各个通道的总压降，以面积为权重求解所有通道的平均压降，其计算公式如下述公式(XIX)：

$\mathrm{\Δ}{p}_{\mathrm{av}}=\frac{\underset{i}{\mathrm{\Σ}}\mathrm{\Δ}{p}_{i,t}{A}_i}{\underset{i}{\mathrm{\Σ}}{A}_i}---\left(\mathrm{XIX}\right)$

第三步，判断平均压降的收敛性，若满足收敛条件，即

时，完成流量分配，其中，收敛准则ε取值为10^-5，Δp_i是指各个通道的总压降；

第四步，计算不同时刻燃料芯块温度及包壳温度，其计算公式如下述公式(XX)：

$\frac{\∂T}{\∂t}=\frac{\mathrm{\ΔQ}}{\mathrm{\ρ}\·c_p}---\left(\mathrm{XX}\right)$

其中，ρ表示密度(kg/m³)，

c_p表示定压比热容(J/(kg·K))，

ΔQ表示燃料芯块内热源释热量与外表面热通量的差值，或者包壳内、外表面的热通量差值(W)；

第五步，计算冷却剂与慢化剂的温度分布，其计算公式如下述公式(XXI)：

Q＝(T₁-T₂)/ΣR  (XXI)

其中，Q表示燃料棒外表面或水棒内、外壁面的释热功率(W)；

T₁表示传热过程起始点的温度(℃)；

T₂表示传热过程终止点的温度(℃)；

ΣR表示换热热阻(W/℃)，

第六步，求解工质温度分布（请补充求解过程），

冷却剂通道内的工质流动满足如下三大守恒方程：

质量守恒方程(XXII)：

$\frac{\∂w_t}{\∂t}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\∂w_i}{\∂t}---\left(\mathrm{XXII}\right)$

动量守恒方程(XXIII)：

$\frac{1}{A_i}\frac{\∂w_i}{\∂t}+\frac{\∂}{\∂z}\left(\frac{{w_i}^2}{{{\mathrm{\ρA}}_i}^2}\right)=-\frac{\∂P}{\∂z}\±\mathrm{\ρ}\stackrel{\‾}{g}-\frac{f{w}_i\left|w_i\right|}{2d_e{{\mathrm{\ρA}}_i}^2}---\left(\mathrm{XXIII}\right)$

能量守恒方程(XXIV)：

$\mathrm{\ρ}\frac{\∂h}{\∂t}+\frac{\∂}{\∂z}\left(\frac{w_{i}h}{\mathrm{\ρ}{A}_i}\right)=q_v^{\mathrm{channel}}\left(t\right)---\left(\mathrm{XXIV}\right)$

其中，z表示轴向节点高度(m)；

t表示时间(s)；

ρ表示冷却剂密度(kg/m³)；

f表示摩擦阻力系数；

h表示冷却剂焓值(J/kg)；

w_t表示冷却剂总流量(kg/s)；

w_i表示冷却剂通道i的流量(kg/s)；

A_i表示冷却剂通道i的流通面积(m²)；

d_e表示水力当量直径(m)；

表示t时刻单位体积工质所吸收的热量(W/m³)。

在一个优选的实施方式中，建立瞬态流量分配计算模型的过程包括如下步骤：

将动量守恒方程(XXIII)沿通道长度进行积分，得到各冷却剂通道内部的通道流量瞬态变化方程组(XXVI)：

$\frac{L}{A_1}\frac{{\∂w}_1}{\∂t}=\mathrm{\Δ}{p}_{t,1}\mathrm{\Δ}{p}_{a,1}-\mathrm{\Δ}{p}_{f,1}-\mathrm{\Δ}{p}_{e1,1}$

$\frac{L}{A_2}\frac{{\∂w}_2}{\∂t}={\mathrm{\Δp}}_{t,2}-\mathrm{\Δ}{p}_{a,2}-{\mathrm{\Δp}}_{f,2}-{\mathrm{\Δp}}_{e1,2}$

……

$\frac{L}{A_{14}}\frac{\∂w_{14}}{\∂t}={\mathrm{\Δp}}_{t,14}-{\mathrm{\Δp}}_{a,14}-{\mathrm{\Δp}}_{f,14}-{\mathrm{\Δp}}_{\mathrm{el},14}---\left(\mathrm{XXVI}\right)$

其中，L表示通道长度(m)，

Δp_t表示并联通道的总压降(Pa)，

Δp_a表示加速压降(Pa)，

Δp_f表示摩擦压降(Pa)，

Δp_el表示提升压降(Pa)，

将质量守恒方程(XXII)和通道流量瞬态变化方程组(XXVI)联立，得到瞬态流量分配计算模型。

瞬态耦合模块，其用于接收每个大时间步长内，瞬态物理模块和瞬态热工模块传出信息，对接收到的信息进行转换，并将转换所得的数据信息分别传给瞬态热工模块和瞬态物理模块；

在一个优选的实施方式中，瞬态核热耦合模块包括中子截面数据库，

中子截面数据库接收瞬态物理模块计算得到的中子通量分布信息，根据中子通量分布信息计算出轴向功率分布信息，并将轴向功率分布信息输入至瞬态热工模块；

中子截面数据库接收瞬态热工模块计算得到的轴向冷却剂温度分布信息、慢化剂温度分布信息，根据轴向冷却剂温度分布信息、慢化剂温度分布信息计算出中子截面分布信息，并将中子截面分布信息输入至瞬态物理模块；

中子截面数据库接收瞬态物理模块计算得到的幅值函数，根据幅值函数计算出总功率，并将总功率输入至瞬态热工模块。

瞬态核热耦合模块包括中子截面数据库，中子截面数据库接收瞬态物理模块传递的截面数值，根据截面数值计算出轴向功率分布，并向瞬态热工模块传递轴向功率分布；

中子截面数据库接收瞬态热工模块传递的轴向工质温度信息和工质密度信息，并将工质温度信息和工质密度信息传递到瞬态物理模块。

在一个优选的实施方式中，瞬态分析模块用于接收瞬态物理模块在每个大时间步长输出的轴向功率分布信息和瞬态热工模块在每个大时间步长输出的工质温度分布信息和工质密度分布信息，并将接收到的信息与安全系数信息进行比较，若接收到的信息超出允许的安全范围，发出指令控制降低反应堆反应速率，并报警；若接收到的信息在允许的安全范围内，待接收到下一个大时间步长的信息后继续与安全系数信息进行比较。

在一个优选的实施方式中，安全系数信息是指燃料包壳允许达到的最大温度、燃料芯块中心最高温度和最高系统压力，其中，燃料包壳允许达到的最大温度是1260℃，燃料芯块中心允许达到的最高温度是2590℃，系统压力允许达到的最大压力值是30.3MPa，即通过本发明提供的方法，控制系统的各项参数在允许的范围内。

如图3中所示，本发明中，瞬态热工模块中的大时间步长和瞬态物理模块中的大时间步长时间长度相等，大时间步长1选取的时间长度为0.1s；中等时间步长2选取的时间长度为0.01s；小时间步长3选取的时间长度为0.00002s。相同种类的时间步长之间没有间隔，各个大时间步长顺序相连，各个中等时间步长顺序相连，各个小时间步长顺序相连，且一个大时间步长的启始时刻也是一个中等时间步长的启始时刻，还是一个小时间步长的启始时刻。

在一个优选的实施方式中，由于计算某一瞬态发生时堆芯各种参数的变化，所以必须先进行相应的稳态计算，以给出瞬态发生之前的参数作为瞬态计算的边界条件，稳态计算是瞬态计算的前提和基础，只有完成了稳态计算才能进行相应的瞬态计算。

在一个优选的实施方式中，如图1所示，多通道瞬态耦合的耦合过程包括，瞬态物理模块的计算，其用于时空动力学方程求解，其将计算出的轴向功率分布信息传输给瞬态分析模块，同时将轴向功率分布信息和中间过程信息输出，由瞬态耦合模块处理，进行参数的更新；根据瞬态定义更新参数后，将参数信息输出给瞬态热工模块，首先通过建立动态流量分配模型，完成不同时刻各个通道之间的入口流量分配；通道内的瞬态热工模型，用于完成不同时刻各个通道内部的换热计算以及流动计算，输出工质温度分布信息和工质密度分布信息给瞬态分析模块，同时将工质温度分布信息和工质密度分布信息和中间过程信息输出给瞬态耦合模块，由瞬态耦合模块计算处理，并将处理后的数据输出给瞬态物理模块；瞬态耦合模块，用于完成瞬态物理模块计算的数据和瞬态热工模块计算的数据之间的链接，从而完成瞬态耦合计算过程，每个大时间步长内完成一次耦合模块的数据传递；其中，根据瞬态定义更新参数部分，属于瞬态耦合模块部分。

本发明的目的在于提供一种用于控制超临界水堆瞬态核热耦合安全的方法，该方法通过用于超临界水堆瞬态核热耦合的安全分析计算装置实现。

本发明的目的在于提供一种用于超临界水堆瞬态核热耦合的程序计算装置，该装置包括上述瞬态物理模块、瞬态热工模块和瞬态耦合模块。

本发明的目的还可以在于提供一种用于超临界水堆瞬态核热耦合的程序计算方法，该方法通过用于超临界水堆瞬态核热耦合的程序计算装置实现。

以上接合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。本发明的保护范围以所附权利要求为准。

Claims (7)

1.一种用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，该装置包括：瞬态物理模块、瞬态热工模块、瞬态耦合模块和瞬态分析模块；

瞬态物理模块，在启动后的第一个大时间步长内，所述瞬态物理模块接收从稳态物理装置传递来的信息，所述信息包括中子通量分布信息、中子密度信息和中子价值信息，根据接收到的信息计算得到轴向功率分布信息，并传出所述轴向功率分布信息和计算过程中得到的中间过程信息，所述中间过程信息包括中子通量分布信息、先驱核浓度信息和功率值信息；在第一个大时间步长结束后，在第二个或其后的大时间步长内，所述瞬态物理模块接收从瞬态耦合模块传递来的信息，所述信息包括中子通量分布信息、中子密度信息和中子价值信息，根据接收到的信息计算得到轴向功率分布信息，并传出所述轴向功率分布信息和计算过程中得到的中间过程信息，所述中间过程信息包括中子通量分布信息、先驱核浓度信息和功率值信息，进入下一个大时间步长时重复该过程；

瞬态热工模块，在启动后的第一个大时间步长内，所述瞬态热工模块接收稳态热工参数信息，所述稳态热工参数信息包括稳态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，根据接收到的稳态热工参数信息计算得到工质温度分布信息和工质密度分布信息，并传出计算得到工质温度分布信息和工质密度分布信息以及计算过程中得到的中间过程信息，所述中间过程信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度；在第一个大时间步长结束后，在第二个或其后的大时间步长内，所述瞬态热工模块接收从瞬态耦合模块传递来的瞬态热工参数信息，所述瞬态热工参数信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，根据接收到的瞬态热工参数信息计算得到工质温度分布信息和工质密度分布信息，并传出计算得到工质温度分布信息和工质密度分布信息以及计算过程中得到的中间过程信息，所述中间过程信息包括瞬态工况下的燃料芯块温度、包壳表面温度、冷却剂温度、慢化剂温度、冷却剂密度和慢化剂密度，进入下一个大时间步长时重复该过程；

瞬态耦合模块，其用于接收每个大时间步长内，瞬态物理模块和瞬态热工模块传出信息，对接收到的信息进行转换，并将转换所得的数据信息分别传给瞬态热工模块和瞬态物理模块；

瞬态分析模块，其用于接收瞬态物理模块在每个大时间步长输出的轴向功率分布信息和瞬态热工模块在每个大时间步长输出的工质温度分布信息和工质密度分布信息，并将接收到的信息与安全系数信息进行比较，若接收到的信息超出允许的安全范围，发出指令控制降低反应堆反应速率，并报警；若接收到的信息在允许的安全范围内，待接收到下一个大时间步长的信息后继续与安全系数信息进行比较。

2.根据权利要求1所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，在瞬态物理模块中，轴向功率分布通过包括如下步骤的过程计算得到：

第一步，计算出小时间步长下不同时刻的堆芯裂变总功率，包括以下子步骤：

a.设定中子通量密度仅沿堆芯轴向变化，建立如下双群中子时空动力学方程(I)和(II)：

$\begin{array}{c}\frac{1}{V_g}\frac{\∂{\mathrm{\Φ}}_g(z,t)}{\∂t}=\\ \▿D_g{\▿\mathrm{\Φ}}_g(z,t)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Φ}}_g(z,t)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}{+\mathrm{\Σ}}_{g^{\′}\→g}\right]{\mathrm{\Φ}}_{g^{\′}}(z,t)+x_g\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)\end{array}---\left(I\right)$

$\frac{{\∂C}_I(z,t)}{\∂t}=\frac{1}{k_{\mathrm{eff}}}{\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Φ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(\mathrm{II}\right)$

其中g＝1,2I＝1,2,...6；

b.把中子通量密度用与时间相关的幅函数以及与空间相关的形状函数的乘积表示，其表达式如下式(III)：

Φ_g(z,t)＝n(t)Ψ_g(z,t)  (III)

进而得到如下双群中子时空动力学时空离散方程(IV)和(V)：

$\begin{array}{c}\frac{1}{V_g}\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}\mathrm{\Ψ}(z,t)+\frac{1}{V_g}\frac{d{\mathrm{\Ψ}}_g(z,t)}{\mathrm{dt}}n\left(t\right)=\\ n\left(t\right)\▿D_g\▿{\mathrm{\Ψ}}_g(z,t)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,t)n\left(t\right)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{\mathrm{\χ}}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\Ψ}}_{g^{\′}}(z,t)n\left(t\right)+x_g\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)\end{array}---\left(\mathrm{IV}\right)$

$\frac{d{C}_I(z,t)}{\mathrm{dt}}=\frac{1}{k_{\mathrm{eff}}}{n\left(t\right)\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(V\right)$

其中，g＝1,2I＝1,2,...6；

c.将双群中子时空动力学方程(I)和(II)两边同时乘以中子价值，并且沿轴向进行积分，合并化简，得出如下点堆中子动力学方程(VI)和(VII)：

$\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}=\frac{\mathrm{\ρ}\left(t\right)-\stackrel{\‾}{\mathrm{\β}}\left(t\right)}{\mathrm{\Λ}\left(t\right)}n\left(t\right)+\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I\left(t\right)---\left(\mathrm{VI}\right)$

$\frac{d{\stackrel{\‾}{C}}_I\left(t\right)}{\mathrm{dt}}=-{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I\left(t\right)+\frac{{\stackrel{\‾}{\mathrm{\β}}}_I\left(t\right)}{\mathrm{\Λ}\left(t\right)}n\left(t\right)---\left(\mathrm{VII}\right)$

其中，I=1,2,…6

式中n(t)表示中子密度(中子/cm³)，

表示总的缓发中有效份额，

表示总的先驱核数(原子/cm³)，

Λ表示瞬发中子每代时间(s)，

设定中子通量总价值保持不变，向后差分，对点堆中子动力学方程(VI)和(VII)进行时间离散，得到如下方程(VIII)和(IX)：

$\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}}=\frac{\mathrm{\ρ}(j+1)-\stackrel{\‾}{\mathrm{\β}}(j+1)}{\mathrm{\Λ}(j+1)}n(j+1)+\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{\stackrel{\‾}{C}}_I(j+1)---\left(\mathrm{VIII}\right)$

$\frac{C_I(j+1)-C_I\left(j\right)}{\mathrm{\Δt}}=-{\mathrm{\λ}}_I{C}_I(j+1)+\frac{{\mathrm{\β}}_I(j+1)}{\mathrm{\Λ}(j+1)}N(j+1)---\left(\mathrm{IX}\right)$

将方程(VIII)和(IX)化简，得到如下的功率计算式(X)和缓发中子先驱核浓度(XI)：

$n(j+1)=\frac{1}{1-\frac{\mathrm{\ρ}(j+1)-\mathrm{\β}}{\mathrm{\Λ}}\mathrm{\Δt}}n\left(t\right)+\frac{\mathrm{\Δt}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(j+1)}{1-\frac{\mathrm{\ρ}(j+1)-\mathrm{\β}}{\mathrm{\Λ}}\mathrm{\Δt}}---\left(X\right)$

$C_I(j+1)=\frac{C_I\left(j\right)}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}+\frac{{\mathrm{\β}}_I(j+1)\mathrm{\Δt}}{\mathrm{\Λ}(j+1)(1+{\mathrm{\λ}}_I\mathrm{\Δt})}N(j+1)---\left(\mathrm{XI}\right)$

进而计算出小时间步长下不同时刻的堆芯裂变总功率；

第二步，使用100%稳态工况时的径向功率分布因子进行稳态计算；

第三步，计算得出堆芯功率轴向分布，包括以下子步骤：

a.化简第一步中的公式(IV)和公式(V)，得到的形状函数时空离散方程(XII)和缓发中子先驱核浓度时空离散方程(XIII)，组成形状函数方程组

$\begin{array}{c}\frac{1}{V_g}\frac{d{\mathrm{\Ψ}}_g(z,t)}{\mathrm{dt}}=\\ \▿D_g\▿{\mathrm{\Ψ}}_g(z,t)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,t)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\Ψ}}_{g^{\′}}(z,t)+\frac{x_g}{n\left(t\right)}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,t)\\ -\frac{1}{V_{g}n\left(t\right)}\frac{\mathrm{dn}\left(t\right)}{\mathrm{dt}}\mathrm{\Ψ}(z,t)\end{array}---\left(\mathrm{XII}\right)$

其中，g＝1,2

$\frac{{\∂C}_I(z,t)}{\∂t}=\frac{1}{k_{\mathrm{eff}}}{n\left(t\right)\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,t)-{\mathrm{\λ}}_I{C}_I(z,t)---\left(\mathrm{XIII}\right)$

其中，I＝1,2,...6

b向后差分，对形状函数方程组进行离散求解，得到如下形状函数时间离散方程(XIV)和(XV):

$\begin{array}{c}\frac{1}{V_g}\frac{{\mathrm{\Ψ}}_g(z,j+1)-{\mathrm{\Ψ}}_g(z,j)}{\mathrm{\Δt}}=\\ \frac{d}{\mathrm{dz}}{D}_g\frac{d}{\mathrm{dz}}{\mathrm{\Ψ}}_g(z,j+1)-{\mathrm{\Σ}}_{r,g}{\mathrm{\Ψ}}_g(z,j+1)+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}[\frac{1}{k_{\mathrm{eff}}}{x}_g(1-\mathrm{\β})v_g{\mathrm{\Σ}}_{\mathrm{fg}}+{\mathrm{\Σ}}_{g^{\′}\→g}]{\mathrm{\Ψ}}_{g^{\′}}(z,j+1)---\left(\mathrm{XIV}\right)\\ +\frac{x_g}{n(j+1)}\underset{I}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_I{C}_I(z,j+1)-\frac{1}{V_{g}n(j+1)}\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}}\mathrm{\Ψ}(z,j+1)\end{array}$

$\frac{C_I(z,j+1)-C_I(z,j)}{\mathrm{\Δt}}=\frac{1}{k_{\mathrm{eff}}}n(j+1){\mathrm{\β}}_I\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left[v{\mathrm{\Σ}}_{f,g^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(z,j+1)\right]-{\mathrm{\λ}}_I{C}_I(z,j+1)---\left(\mathrm{XV}\right)$

c.将形状函数时间离散方程组的第二个方程(XV)进行化简，得到缓发中子先驱核浓度的离散求解表达式，然后将其代入形状函数时间离散方程组的第一个方程(XIV)进行化简，得到如下各个轴向节点i处的形状函数求解表达式(XVI)：

${\mathrm{\Ψ}}^i(j+1)=\frac{[S_{\text{g}}^{\′\′i-1/2}(j+1)+S_g^{\′\′i+1/2}(j+1)]\·\frac{{{\mathrm{\Δ}}_z}^2}{2}+D^{i-1/2}{\mathrm{\Ψ}}^{i-1}(j+1)+D^{i+1/2}{\mathrm{\ψ}}^{i+1}(j+1)}{D^{i+1/2}+D^{i-1/2}+({\mathrm{\Σ}}_r^{i-1/2}+{\mathrm{\Σ}}_r^{i+1/2})\·\frac{{{\mathrm{\Δ}}_z}^2}{2}+(\frac{1}{V_{g}n(j+1)}\·\frac{n(j+1)-n\left(j\right)}{\mathrm{\Δt}})\·{{\mathrm{\Δ}}_z}^2}---\left(\mathrm{XVI}\right)$

其中， $S_g^{\′\′}(j+1)=\left[\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}\left(v{\mathrm{\Σ}}_{{\mathrm{fg}}^{\′}}{\mathrm{\Ψ}}_{g^{\′}}(j+1)\right)\right]\·[(1-\mathrm{\β})x_g+x_g\underset{I=1}{\overset{6}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_I{\mathrm{\β}}_I\mathrm{\Δt}}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}]\·\frac{1}{k_{\mathrm{eff}}}$

$+\frac{1}{V_g\mathrm{\Δt}}{\mathrm{\Ψ}}_g\left(j\right)+\frac{1}{n(j+1)}\underset{I=1}{\overset{6}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_I{x}_g{C}_I(z,j)}{1+{\mathrm{\λ}}_I\mathrm{\Δt}}+\underset{g^{\′}=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\Σ}}_{g^{\′}\→}{\mathrm{\Ψ}}_{g^{\′}}(j+1)---\left(\mathrm{XVII}\right)$

上述公式XVII为中子产生项，可以用来计算产生中子的数量，

对得到形状函数进行归一化判断，判断准则为中子总价值的计算值与基准值之间的相对误差小于10^-5；每一个中等时间步长进行一次归一化判断，若形状函数符合归一化标准，则用该形状方程的解去求解中子通量和缓发中子先驱核浓度，然后转入下一个时刻的计算；若形状方程的解不满足归一化条件，则进行“形状函数迭代更新”，并重复求解和判定过程；

d.根据由点堆中子动力学方程求解得到的幅值函数，以及形状函数，得到如下用于计算不同时间步长下轴向功率分布的公式(XVIII)：

q_l ⁱ(j+1)＝n(j+1)·Ψⁱ(j+1)·ξ₂  (XVIII)

其中，q_l ⁱ(j+1)表示在j+1时刻轴向节点i的燃料棒线功率密度(W/m)；

ξ₂表示时空离散的功率修正因子。

3.根据权利要求1所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，在瞬态热工模块中，工质温度分布和工质密度分布通过包括如下步骤的过程计算得到：

第一步，建立瞬态流量分配计算模型，进行流量分配；

第二部，计算各个通道的总压降，以面积为权重求解所有通道的平均压降，其计算公式如下述公式(XIX)：

$\mathrm{\Δ}{p}_{\mathrm{av}}=\frac{\underset{i}{\mathrm{\Σ}}\mathrm{\Δ}{p}_{i,t}{A}_i}{\underset{i}{\mathrm{\Σ}}{A}_i}---\left(\mathrm{XIX}\right)$

第三步，判断平均压降的收敛性，若满足收敛条件，即

时，完成流量分配，其中，收敛准则ε取值为10^-5，Δp_i是指各个通道的总压降；

第四步，计算不同时刻燃料芯块温度及包壳温度，其计算公式如下述公式(XX)：

$\frac{\∂T}{\∂t}=\frac{\mathrm{\ΔQ}}{\mathrm{\ρ}\·c_p}---\left(\mathrm{XX}\right)$

其中，ρ表示密度(kg/m³)，

c_p表示定压比热容(J/(kg·K))，

ΔQ表示燃料芯块内热源释热量与外表面热通量的差值，或者包壳内、外表面的热通量差值(W)；

第五步，计算冷却剂与慢化剂的温度分布，其计算公式如下述公式(XXI)：

Q＝(T₁-T₂)/ΣR  (XXI)

其中，Q表示燃料棒外表面或水棒内、外壁面的释热功率(W)；

T₁表示传热过程起始点的温度(℃)；

T₂表示传热过程终止点的温度(℃)；

ΣR表示换热热阻(W/℃)，

第六步，求解工质温度分布，

冷却剂通道内的工质流动满足如下三大守恒方程：

质量守恒方程(XXII)：

$\frac{\∂w_t}{\∂t}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{\∂w_i}{\∂t}---\left(\mathrm{XXII}\right)$

动量守恒方程(XXIII)：

$\frac{1}{A_i}\frac{\∂w_i}{\∂t}+\frac{\∂}{\∂z}\left(\frac{{w_i}^2}{{{\mathrm{\ρA}}_i}^2}\right)=-\frac{\∂P}{\∂z}\±\mathrm{\ρ}\stackrel{\‾}{g}-\frac{f{w}_i\left|w_i\right|}{2d_e{{\mathrm{\ρA}}_i}^2}---\left(\mathrm{XXIII}\right)$

能量守恒方程(XXIV)：

$\mathrm{\ρ}\frac{\∂h}{\∂t}+\frac{\∂}{\∂z}\left(\frac{w_{i}h}{\mathrm{\ρ}{A}_i}\right)=q_v^{\mathrm{channel}}\left(t\right)---\left(\mathrm{XXIV}\right)$

其中，z表示轴向节点高度(m)；

t表示时间(s)；

ρ表示冷却剂密度(kg/m³)；

f表示摩擦阻力系数；

h表示冷却剂焓值(J/kg)；

w_t表示冷却剂总流量(kg/s)；

w_i表示冷却剂通道i的流量(kg/s)；

A_i表示冷却剂通道i的流通面积(m²)；

d_e表示水力当量直径(m)；

表示t时刻单位体积工质所吸收的热量(W/m³)。

4.根据权利要求3所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，建立瞬态流量分配计算模型的过程包括如下步骤：

将动量守恒方程(XXIII)沿通道长度进行积分，得到各冷却剂通道内部的通道流量瞬态变化方程组(XXVI)：

$\frac{L}{A_1}\frac{{\∂w}_1}{\∂t}=\mathrm{\Δ}{p}_{t,1}\mathrm{\Δ}{p}_{a,1}-\mathrm{\Δ}{p}_{f,1}-\mathrm{\Δ}{p}_{e1,1}$

$\frac{L}{A_2}\frac{{\∂w}_2}{\∂t}={\mathrm{\Δp}}_{t,2}-\mathrm{\Δ}{p}_{a,2}-{\mathrm{\Δp}}_{f,2}-{\mathrm{\Δp}}_{e1,2}$

……

$\frac{L}{A_{14}}\frac{\∂w_{14}}{\∂t}={\mathrm{\Δp}}_{t,14}-{\mathrm{\Δp}}_{a,14}-{\mathrm{\Δp}}_{f,14}-{\mathrm{\Δp}}_{\mathrm{el},14}---\left(\mathrm{XXVI}\right)$

其中，L表示通道长度(m)，

Δp_t表示并联通道的总压降(Pa)，

Δp_a表示加速压降(Pa)，

Δp_f表示摩擦压降(Pa)，

Δp_el表示提升压降(Pa)，

将质量守恒方程(XXII)和通道流量瞬态变化方程组(XXVI)联立，得到瞬态流量分配计算模型。

5.根据权利要求1所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，瞬态核热耦合模块包括中子截面数据库，

中子截面数据库接收瞬态物理模块计算得到的中子通量分布信息，根据中子通量分布信息计算出轴向功率分布信息，并将轴向功率分布信息输入至瞬态热工模块；

中子截面数据库接收瞬态热工模块计算得到的轴向冷却剂温度分布信息、慢化剂温度分布信息，根据轴向冷却剂温度分布信息、慢化剂温度分布信息计算出中子截面分布信息，并将中子截面分布信息输入至瞬态物理模块；

中子截面数据库接收瞬态物理模块计算得到的幅值函数，根据幅值函数计算出总功率，并将总功率输入至瞬态热工模块。

瞬态核热耦合模块包括中子截面数据库，中子截面数据库接收瞬态物理模块传递的截面数值，根据截面数值计算出轴向功率分布，并向瞬态热工模块传递轴向功率分布；

中子截面数据库接收瞬态热工模块传递的轴向工质温度信息和工质密度信息，并将工质温度信息和工质密度信息传递到瞬态物理模块。

6.根据权利要求1所述的用于超临界水堆瞬态核热耦合的安全分析计算装置，其特征在于，安全系数信息是指燃料包壳允许达到的最大温度、燃料芯块中心最高温度和最高系统压力，其中，燃料包壳允许达到的最大温度是1260℃，燃料芯块中心允许达到的最高温度是2590℃，系统压力允许达到的最大压力值是30.3MPa。

7.一种用于控制超临界水堆瞬态核热耦合安全的方法，其特征在于，该方法通过用于超临界水堆瞬态核热耦合的安全分析计算装置实现。

Images