一种适用于核反应堆的输运燃耗耦合计算的方法，包括如下步骤：1、由燃耗步初的原子核密度进行输运计算，得到粗网参数、微观反应率；2、由微观反应率和原子核密度进行燃耗计算，得到燃耗步末预估的原子核密度；3、由预估的原子核密度进行输运计算，得到粗网参数、微观反应率；4、燃耗步内划分燃耗CMFD子步，线性插值保存的粗网参数；5、由步骤4的粗网参数和步骤3的细网中子通量，更新燃耗CMFD子步的微观反应率；6、用燃耗CMFD子步上的微观反应率进行燃耗计算，得到燃耗步末准确的原子核密度；7、判定燃耗步数目是否与输入值一致以判定计算是否结束；本发明方法在保证极高精度的前提下，极大地扩大燃耗计算的步长，减小核反应堆整个寿期内计算的时间。
1.一种适用于核反应堆的输运燃耗耦合计算的方法，其特征在于：包括如下步骤:步骤1：根据燃耗步初的原子核密度和输入卡片提供的几何信息进行输运计算，输运计算时引入粗网加速方法即CMFD加速方法，保存燃耗步初粗网加速计算时得到的粗网参数，得到燃耗步初的各个燃耗区的微观反应率即各个燃耗区的核素的微观截面与相应燃耗区中子标通量密度的乘积，其中输运计算能够得到细网中子角通量密度，粗网加速方法能够更新粗网的中子标通量密度，来加快细网中子角通量密度的迭代收敛；具体包括如下步骤：
1)从数据库中读取各个核素的原始多群微观截面信息；
2)从输入卡片中读取材料中包含的核素原子核密度、温度和几何信息；
3)基于1)和2)所得到的各个核素的原始多群微观截面、原子核密度和温度信息，采用子群方法进行共振自屏计算，得到各个核素的多群有效自屏截面，具体的计算公式如下：式中：
σx,g,iso--核素iso、能群g、反应道χ的截面；
σx,iso(E)--核素iso、能量E处的反应道χ的截面；
4)根据3)得到的信息采用MOC特征线方法进行中子输运计算，得到各个区域的中子角通量密度，具体的计算公式如下：式中：
ψg(r,Ω)——第g能群区域r、角度Ω的中子角通量密度；
5)对公式(2)中计算得到的中子角通量密度进行全角度空间上积分，得到细网上的平均中子标通量密度 具体的计算公式如下：式中：
6)在进行MOC的迭代计算时，求解粗网平衡方程对应的粗网差分方程，得到粗网通量，计算公式如公式(4)所示；通过粗网差分计算前后粗网格平均中子标通量密度之比同步更新粗网格内所有细网格上的中子标通量密度，用来更新下一次MOC计算时方程右端的源项，能够加速MOC的迭代计算，计算公式如公式(5)所示；在MOC计算结束之后，保存此时粗网差分方程的参数，具体的计算公式如下：式中：
——在粗网网格n上，第g能群，方向为u的正方向上的净中子流；
——在粗网网格n上，第g能群，方向为u的负方向上的净中子流；
——在粗网网格n上，第g能群的粗网上的中子通量密度；
——在粗网网格n上，第g′能群的粗网上的中子通量密度；
——在粗网网格n上，从第g′能群散射到第g能群的散射截面；
——更新后的平源细区i第g群的平均标通量，用以更新第l+1次矩阵MOC计算的源项；
——第l次矩阵MOC计算后的CMFD计算得到的粗网n的第g群平均标通量；
——第l次矩阵MOC计算得到的粗网n的第g群均匀化标通量；
——第l次矩阵MOC计算所得到的平源细区i第g群的平均标通量；
7)由3)、4)、5)分别得到需要的各个核素的微观截面与相应区域中子标通量密度的乘积，即各个区域的各个核素的微观反应率；
步骤2：使用步骤1计算得到的燃耗步初的微观反应率和原子核密度进行燃耗计算，燃耗计算时，在两个燃耗步内划分燃耗子步，在每个燃耗子步内只进行燃耗计算，不进行输运计算，进行微观反应率的更新，得到燃耗步末的原子核密度，该原子核密度称为预估步的原子核密度；具体包括如下步骤：
1)从燃耗数据库中读取各个核素的衰变常数、反应道分支比、裂变核素的裂变产额；
2)读取在输运计算中保存下来的各个区域的各个核素的原子核密度，以及输运计算得到的各个区域各个核素的微观反应率；
3)根据1)、2)得到的信息，采用切比雪夫有理近似方法求解燃耗方程，具体的计算公式如下：式中：
步骤3：使用步骤2中得到的预估步的原子核密度进行输运计算，计算公式如公式(1)、公式(2)和公式(3)所示，输运计算时引入粗网加速方法即CMFD加速方法，计算公式如公式(4)和公式(5)所示，保存燃耗步末粗网计算时得到的粗网参数，得到燃耗步末的微观反应率即核素的微观截面与中子标通量密度的乘积，其中中子通量密度是细网上的中子标通量密度；
步骤4：在两个燃耗步中间划分燃耗CMFD子步，在每个燃耗CMFD子步上，对燃耗步初和燃耗步末的粗网参数进行线性插值，得到每个燃耗CMFD子步上的粗网参数；
步骤5：在每个燃耗CMFD子步上，使用步骤4中插值得到的粗网参数，对公式(4)进行求解，进行粗网计算，更新粗网的通量，并使用步骤3得到的燃耗步末的细网上的中子标通量密度，依照公式(5)更新每个燃耗CMFD子步上的细网上的中子标通量密度，进而更新每个燃耗CMFD子步上的微观反应率；
步骤6：使用步骤5中每个燃耗CMFD子步上更新的微观反应率即核素微观截面与更新后的中子标通量密度的乘积，通过公式(6)进行燃耗计算，并且相邻燃耗CMFD子步间，同样划分燃耗子步进行计算，最终得到燃耗步末的原子核密度，该原子核密度即需要的燃耗步末准确的原子核密度；
步骤7：判定燃耗步数目是否与输入值一致，如果燃耗步与输入值一致，则计算结束；如果不一致，返回至步骤1，进入到下一个燃耗步的计算。
一种适用于核反应堆的输运燃耗耦合计算的方法
技术领域
[0001] 本发明涉及核反应堆物理计算技术领域，具体涉及一种适用于核反应堆的输运燃耗耦合计算的方法。
背景技术
[0002] 在设计条件和计算精度要求不断提高的今天，核反应堆设计领域希望通过精确建模，并最大程度的减少上游核反应堆计算过程中引入的近似与误差，实现高保真的反应堆物理计算，从而对反应堆进行精确分析模拟。
[0003] 高保真的反应堆物理计算主要包括中子输运计算和燃耗计算两部分。中子输运计算根据原子核密度场、温度场等得到中子通量密度场，进而提供堆芯功率分布；燃耗计算根据某一时刻输运计算得到的微观中子反应率，得到该时刻的原子核密度场，如此进行不断迭代，进而得到核反应堆的功率分布等随时间的变化规律。
[0004] 在高保真的反应堆物理计算中，中子输运计算极为耗时，因此输运燃耗耦合策略至关重要。优秀的输运燃耗耦合策略，能够在保证极高精度的前提下，有效地减少整个燃耗寿期内，中子输运计算的次数，进而极大地减少整个反应堆物理计算的时间。
[0005] 燃耗计算过程中，最为常见的输运燃耗耦合策略是起点近似方法。所谓的起点近似方法就是认为两种燃耗步间的微观反应率(核素的微观截面和中子标通量密度的乘积)保持不变，然后从寿期初一直计算到寿期末。由于在整个寿期中，微观反应率在不断变化，显然起点近似只有划分极细的燃耗步，才能获得一个准确的结果。在起点近似的基础上，诞生了传统意义上的带预估校正的中子输运燃耗耦合策略。带预估校正的输运燃耗耦合策略能较大地扩大燃耗步长，并且保持计算的精度，并且几种预估校正的输运耦合策略的效果极为接近。
[0006] 在传统的带预估校正的输运耦合策略中，最常用的一种预估校正输运燃耗耦合策略是假设当前燃耗步长内的微观反应率(核素的微观截面和中子标通量密度的乘积)与原子核密度是线性变化关系：首先当前燃耗步初始时刻的原子核密度进行输运计算得到初始时刻的微观反应率，然后用该微观反应率求解当前燃耗步长内的燃耗方程得到燃耗步长末的原子核密度，称为预估步的原子核密度，然后用预估步的原子核密度再进行输运计算，得到燃耗步长末的微观反应率，用该微观反应率再求解一次当前燃耗步的燃耗方程，得到燃耗步长末的原子核密度，称为校正步的原子核密度。将预估步的原子核密度和校正步的原子核密度的平均值作为下一个燃耗步的原子核密度。
[0007] 传统的输运燃耗耦合策略，在组件或者堆芯布满毒物时，必须将燃耗的步长控制在很小的值，才能获得一个理想的结果，却也因此需要进行很多次中子输运计算，这在实际的高保真的反应堆物理计算中是极其耗时的，大大地延长了计算时间。
[0008] 预估校正方法在常规的不带可燃毒物棒的燃料组件计算中，精度是较高的，但是带可燃毒物棒的燃料组件乃至堆芯计算时，为了保证预估校正方法的准确性，必须将划分较细的燃耗步，因此在整个寿期的燃耗计算中，会进行多次的输运计算，大大增加了计算时间。发明内容：
[0009] 为了克服上述现有技术存在的问题，本发明的目的在于提供一种新的适用于核反应堆的输运燃耗耦合计算的方法，在保证极高精度的前提下，极大地扩大燃耗计算的步长，减少高保真的核反应堆物理输运计算的时间，从而减小核反应堆整个寿期内计算的时间。
[0010] 为了实现上述目的，本发明将采取以下技术方案予以实施：
[0011] 一种适用于核反应堆的输运燃耗耦合计算的方法，包括如下步骤:
[0012] 步骤1：根据燃耗步初的原子核密度和输入卡片提供的几何信息进行输运计算，输运计算时引入粗网加速方法即CMFD加速方法，保存燃耗步初粗网加速计算时得到的粗网参数，得到燃耗步初的各个燃耗区的微观反应率即各个燃耗区的核素的微观截面与相应燃耗区中子标通量密度的乘积，其中输运计算能够得到细网中子角通量密度，粗网加速方法能够更新粗网的中子标通量密度，来加快细网中子角通量密度的迭代收敛；具体包括如下步骤：
[0013] 1)从数据库中读取各个核素的原始多群微观截面信息；
[0014] 2)从输入卡片中读取材料中包含的核素原子核密度、温度和几何信息；
[0015] 3)基于1)和2)所得到的各个核素的原始多群微观截面、原子核密度和温度信息，采用子群方法进行共振自屏计算，得到各个核素的多群有效自屏截面，具体的计算公式如下：
[0016]
[0017] 式中：
[0018] σx,g,iso--核素iso、能群g、反应道χ的截面；
[0019] ΔEg--第g群的能量宽度；
[0020] σx,iso(E)--核素iso、能量E处的反应道χ的截面；
[0021] φ(E)--能量为E的中子通量密度；
[0022] 4)根据3)得到的信息采用MOC特征线方法进行中子输运计算，得到各个区域的中子角通量密度，具体的计算公式如下：
[0023]
[0024] 式中：
[0025] Ω——角度方向；
[0026] ——梯度算子；
[0027] ψg(r,Ω)——第g能群区域r、角度Ω的中子角通量密度；
[0028] g——第g能群；
[0029] G——能群总数；
[0030] Σt,g(r)——第g群的宏观总截面；
[0031] Qg(r,Ω)——中子输运方程源项；
[0032] 由此得到各个计算区域的中子角通量密度；
[0033] 5)对公式(2)中计算得到的中子角通量密度进行全角度空间上积分，得到细网上的平均中子标通量密度 具体的计算公式如下：
[0034]
[0035] 式中：
[0036] M——离散方向总数；
[0037] ωm——求积组
[0038] 且有：
[0039]
[0040] 由此得到各个计算区域的中子标通量密度；6)在进行MOC的迭代计算时，求解粗网平衡方程对应的粗网差分方程，得到粗网通量，计算公式如公式(4)所示；通过粗网差分计算前后粗网格平均中子标通量密度之比同步更新粗网格内所有细网格上的中子标通量密度，用来更新下一次MOC计算时方程右端的源项，能够加速MOC的迭代计算，计算公式如公式(5)所示；在MOC计算结束之后，保存此时粗网差分方程的参数，具体的计算公式如下：
[0041]
[0042] 式中：
[0043] n——差分粗网网格
[0044] u(u＝x,y,z)——方向标识；
[0045] ——差分粗网n在u(u＝x,y,z)方向上的网格宽度；
[0046] ——在粗网网格n上，第g能群，方向为u的正方向上的净中子流；
[0047] ——在粗网网格n上，第g能群，方向为u的负方向上的净中子流；
[0048] ——在粗网网格n上，第g能群，第r种均匀化截面；
[0049] ——在粗网网格n上，第g能群的粗网上的中子通量密度；
[0050] ——在粗网网格n上，第g′能群的产生截面；
[0051] ——在粗网网格n上，第g′能群的粗网上的中子通量密度；
[0052] ——在粗网网格n上，从第g′能群散射到第g能群的散射截面；
[0053] χg——第g能群的裂变谱；
[0054] keff——中子平衡方程的特征值；
[0055]
[0056] 式中：
[0057] ——更新后的平源细区i第g群的平均标通量，用以更新第l+1次矩阵MOC计算的源项；
[0058] ——第l次矩阵MOC计算后的CMFD计算得到的粗网n的第g群平均标通量；
[0059] ——第l次矩阵MOC计算得到的粗网n的第g群均匀化标通量；
[0060] ——第l次矩阵MOC计算所得到的平源细区i第g群的平均标通量；
[0061] 7)由3)、4)、5)分别得到需要的各个核素的微观截面与相应区域中子标通量密度的乘积，即各个区域的各个核素的微观反应率；
[0062] 步骤2：使用步骤1计算得到的燃耗步初的微观反应率和原子核密度进行燃耗计算，燃耗计算时，在两个燃耗步内划分燃耗子步，在每个燃耗子步内只进行燃耗计算，不进行输运计算，进行微观反应率的更新，得到燃耗步末的原子核密度，该原子核密度称为预估步的原子核密度；具体包括如下步骤：
[0063] 1)从燃耗数据库中读取各个核素的衰变常数、反应道分支比、裂变核素的裂变产额等燃耗数据；
[0064] 2)读取在输运计算中保存下来的各个区域的各个核素的原子核密度，以及输运计算得到的各个区域各个核素的微观反应率；
[0065] 3)根据1)、2)得到的信息，采用切比雪夫有理近似方法求解燃耗方程，具体的计算公式如下：
[0066]
[0067] 式中：
[0068] N(t)--原子核密度关于时间的表达式；
[0069] N(0)--初始时刻的原子核密度
[0070] A--燃耗矩阵；
[0071] Γ--曲线积分，逆时针绕矩阵tA所有特征值一圈；
[0072] I--单位矩阵；
[0073] i--虚数；
[0074] 步骤3：使用步骤2中得到的预估步的原子核密度进行输运计算，计算公式如公式(1)、公式(2)和公式(3)所示，输运计算时引入粗网加速方法即CMFD加速方法，计算公式如公式(4)和公式(5)所示，保存燃耗步末粗网计算时得到的粗网参数，得到燃耗步末的微观反应率即核素的微观截面与中子标通量密度的乘积，其中中子通量密度是细网上的中子标通量密度；
[0075] 步骤4：在两个燃耗步中间划分燃耗CMFD子步，在每个燃耗CMFD子步上，对燃耗步初和燃耗步末的粗网参数进行线性插值，得到每个燃耗CMFD子步上的粗网参数；
[0076] 步骤5：在每个燃耗CMFD子步上，使用步骤4中插值得到的粗网参数，对公式(4)进行求解，进行粗网计算，更新粗网的通量，并使用步骤3得到的燃耗步末的细网上的中子标通量密度，依照公式(5)更新每个燃耗CMFD子步上的细网上的中子标通量密度，进而更新每个燃耗CMFD子步上的微观反应率；
[0077] 步骤6：使用步骤5中每个燃耗CMFD子步上更新的微观反应率即核素微观截面与更新后的中子标通量密度的乘积，通过公式(6)进行燃耗计算，并且相邻燃耗CMFD子步间，同样划分燃耗子步进行计算，最终得到燃耗步末的原子核密度，该原子核密度即需要的燃耗步末准确的原子核密度；
[0078] 步骤7：判定燃耗步数目是否与输入值一致，如果燃耗步与输入值一致，则计算结束；如果不一致，返回至步骤1，进入到下一个燃耗步的计算。
[0079] 与常规的方法相比，本发明有如下突出优点：
[0080] 在燃耗计算过程中，不再假设两个燃耗步间的微观反应率为常值。本发明方法在校正步计算时，通过划分燃耗CMFD子步，在每个燃耗CMFD子步上通过线性插值粗网参数参数，进行粗网计算，更新了每个燃耗CMFD子步上的中子通量密度，捕捉到了微观反应率随燃耗的变化，因此可以扩大燃耗步长的划分，减小整个寿期燃耗计算的时间。
附图说明
[0081] 图1为本发明方法流程图。具体实施方式：
[0082] 下面结合附图和具体实施方式对本发明做进一步详细说明：
[0083] 如图1所示，本发明一种适用于核反应堆的输运燃耗耦合计算的方法，包括如下步骤:
[0084] 步骤1：根据燃耗步初的原子核密度和输入卡片提供的几何信息进行输运计算，输运计算时引入粗网加速方法即CMFD加速方法，保存燃耗步初粗网加速计算时得到的粗网参数，得到燃耗步初的各个燃耗区的微观反应率即各个燃耗区的核素的微观截面与相应燃耗区中子标通量密度的乘积，其中输运计算能够得到细网中子角通量密度，粗网加速方法能够更新粗网的中子标通量密度，来加快细网中子角通量密度的迭代收敛；具体包括如下步骤：
[0085] 1)从数据库中读取各个核素的原始多群微观截面信息；
[0086] 2)从输入卡片中读取材料中包含的核素原子核密度、温度和几何信息，输入卡片中已经对求解问题的几何进行详细的区域剖分，包括燃料棒栅元的径向区域划分和周向区域划分，一般情况下针对普通燃料棒栅元，燃料区域径向划分为3区，周向划分为8区；针对带有可燃毒物Gd棒的燃料栅元，由于其具有较强的空间自屏效应，燃料区域径向划分为10区，周向划分为8区，并且标记有裂变材料的区域为燃耗区；
[0087] 3)基于1)和2)所得到的各个核素的原始的多群微观截面、原子核密度和温度信息，采用子群方法进行共振自屏计算，得到各个核素的多群有效自屏截面，具体的计算公式如下：
[0088]
[0089] 式中：
[0090] σx,g,iso--核素iso、能群g、反应道χ的截面；
[0091] ΔEg--第g群的能量宽度；
[0092] σx,iso(E)--核素iso、能量E处的反应道χ的截面；
[0093] φ(E)--能量为E的中子通量密度；
[0094] 4)根据3)得到的信息采用MOC特征线方法进行中子输运计算，得到各个区域的中子角通量密度，具体的计算公式如下：
[0095]
[0096] 式中：
[0097] Ω——角度方向；
[0098] ——梯度算子；
[0099] ψg(r,Ω)——第g能群区域r、角度Ω的中子角通量密度；
[0100] g——第g能群；
[0101] G——能群总数；
[0102] Σt,g(r)——第g群的宏观总截面；
[0103] Qg(r,Ω)——中子输运方程源项；
[0104] 由此得到各个计算区域的中子角通量密度；
[0105] 5)对公式(2)中计算得到的中子角通量密度进行全角度空间上积分，得到细网上的平均中子标通量密度 具体的计算公式如下：
[0106]
[0107] 式中：
[0108] M——离散方向总数；
[0109] ωm——求积组
[0110] 且有：
[0111]
[0112] 由此得到各个计算区域的中子标通量密度；
[0113] 6)在进行MOC的迭代计算时，求解粗网平衡方程对应的粗网差分方程，得到粗网通量，计算公式如公式(4)所示。通过粗网差分计算前后粗网格平均中子标通量密度之比同步更新粗网格内所有细网格上的中子标通量密度，用来更新下一次MOC计算时方程右端的源项，可以加速MOC的迭代计算，计算公式如公式(5)所示。在MOC计算结束之后，保存此时粗网差分方程的参数，具体的计算公式如下：
[0114]
[0115] 式中：
[0116] n——差分粗网网格
[0117] u(u＝x,y,z)——方向标识；
[0118] ——差分粗网n在u(u＝x,y,z)方向上的网格宽度；
[0119] ——在粗网网格n上，第g能群，方向为u的正方向上的净中子流；
[0120] ——在粗网网格n上，第g能群，方向为u的负方向上的净中子流；
[0121] ——在粗网网格n上，第g能群，第r种均匀化截面；
[0122] ——在粗网网格n上，第g能群的粗网上的中子通量密度；
[0123] ——在粗网网格n上，第g′能群的产生截面；
[0124] ——在粗网网格n上，第g′能群的粗网上的中子通量密度；
[0125] ——在粗网网格n上，从第g′能群散射到第g能群的散射截面；
[0126] χg——第g能群的裂变谱；
[0127] keff——中子平衡方程的特征值；
[0128]
[0129] 式中：
[0130] ——更新后的平源细区i第g群的平均标通量，用以更新第l+1次矩阵MOC计算的源项；
[0131] ——第l次矩阵MOC计算后的CMFD计算得到的粗网n的第g群平均标通量；
[0132] ——第l次矩阵MOC计算得到的粗网n的第g群均匀化标通量；
[0133] ——第l次矩阵MOC计算所得到的平源细区i第g群的平均标通量。
[0134] 7)由3)、4)、5)分别得到需要的各个核素的微观截面与相应区域中子标通量密度的乘积，即各个区域的各个核素的微观反应率；
[0135] 步骤2：使用步骤1计算得到的燃耗步初的微观反应率和原子核密度进行燃耗计算，燃耗计算时，在两个燃耗步内划分燃耗子步，在每个燃耗子步内只进行燃耗计算，不进行输运计算，进行微观反应率的更新，得到燃耗步末的原子核密度，该原子核密度称为预估步的原子核密度；具体包括如下步骤：
[0136] 1)从燃耗数据库中读取各个核素的衰变常数、反应道分支比、裂变核素的裂变产额等燃耗数据；
[0137] 2)读取在输运计算中保存下来的各个区域的各个核素的原子核密度，以及输运计算得到的各个区域各个核素的微观反应率；
[0138] 3)根据1)、2)得到的信息，采用切比雪夫有理近似方法求解燃耗方程，具体的计算公式如下：
[0139]
[0140] 式中：
[0141] N(t)--原子核密度关于时间的表达式；
[0142] N(0)--初始时刻的原子核密度
[0143] A--燃耗矩阵；
[0144] Γ--曲线积分，逆时针绕矩阵tA所有特征值一圈；
[0145] I--单位矩阵；
[0146] i--虚数；
[0147] 步骤3：使用步骤2中得到的预估步的原子核密度进行输运计算，计算公式如公式(1)、公式(2)和公式(3)所示，输运计算时引入粗网加速方法即CMFD加速方法，计算公式如公式(4)和公式(5)所示，保存燃耗步末粗网计算时得到的粗网参数，得到燃耗步末的微观反应率即核素的微观截面与中子标通量密度的乘积，其中中子通量密度是细网上的中子标通量密度；
[0148] 步骤4：在两个燃耗步中间划分燃耗CMFD子步，在每个燃耗CMFD子步上，对燃耗步初和燃耗步末的粗网参数进行线性插值，得到每个燃耗CMFD子步上的粗网参数；
[0149] 步骤5：在每个燃耗CMFD子步上，使用步骤4中插值得到的粗网参数，对公式(4)进行求解，进行粗网计算，更新粗网的通量，并使用步骤3得到的燃耗步末的细网上的中子标通量密度，依照公式(5)更新每个燃耗CMFD子步上的细网上的中子标通量密度，进而更新每个燃耗CMFD子步上的微观反应率；
[0150] 步骤6：使用步骤5中每个燃耗CMFD子步上更新的微观反应率即核素微观截面与更新后的中子标通量密度的乘积，通过公式(6)进行燃耗计算，并且相邻燃耗CMFD子步间，同样划分燃耗子步进行计算，最终得到燃耗步末的原子核密度，该原子核密度即需要的燃耗步末准确的原子核密度；
[0151] 步骤7：判定燃耗步数目是否与输入值一致，如果燃耗步与输入值一致，则计算结束；如果不一致，返回至步骤1，进入到下一个燃耗步的计算。
