本发明涉及核工程中燃耗数值模拟领域，具体涉及一种基于线积分有理近似方法求解微扰燃耗灵敏度系数的方法。
背景技术：
：反应堆中能量的释放来自于核素与中子发生反应和放射性核素自发衰变反应。而正是由于这些反应，使得反应堆的核素组成在不断的发生变化。燃耗计算的目标即在给定核素之间转换关系以及相应系数(包括中子反应截面，中子通量，衰变常数，分支比等)，和初始核素原子核密度向量的情况下，得到核素原子核密度向量随时间变化的特性。描述放射性核素自发衰变性质的参数直接来源于核评价数据库；核素中子反应截面参数一般是首先由核评价库加工得到点截面或者多群截面，再经中子输运计算后按照通量能谱归并为问题相关的单群截面；中子通量密度则是根据输运计算得到的各中子反应单群截面和中子反应的反应热，按照给定功率水平进行归一而得到。核评价数据库中存储的核素衰变性质参数，中子反应截面以及反应热来自于实验测量或者理论预测，均存在一定的不确定性，从而会导致燃耗计算结果的不确定性。定义响应为依赖于计算结果函数式，对该响应不确定性的评价可以通过计算不同输入参数的灵敏度系数的方法而实现。经过一系列理论推导，在仅考虑燃耗计算自身时，灵敏度系数的计算可归结为关于燃耗方程、扰动参数和响应的一个积分运算；在考虑燃耗-输运耦合计算时，该积分计算也是灵敏度系数表达式中的一部分。燃耗方程的求解算法有泰勒展开与截断、线性子链法(TTA)、切比雪夫有理近似(CRAM)和线积分有理近似(QRAM)等。泰勒展开与截断是较为经典的方法，具有效率较高的优势，然而该方法应用的前提是短寿期核素单独予以处理以降低燃耗矩阵的范数，导致其求解精度较低。切比雪夫有理近似和线积分有理近似适于带通量的燃耗方程求解，两种方法的有效性基于燃耗矩阵特征值分布于负实轴附近的特点，具有优异的求解精度和效率。线性子链法在衰变燃耗方程的求解上优于其他方法。灵敏度系数计算所涉及的积分运算，目前通常采用数值积分法，划分计算时间区间为若干小区间，每个小区间进行两次燃耗计算，并存储计算结果。为了达到计算精度需求，往往需要划分上百个小区间，因而产生了巨大的时间和存储开销。技术实现要素：为了克服上述现有技术存在的问题，本发明的目的在于提供一种基于线积分有理近似方法求解微扰燃耗灵敏度系数的方法，线积分有理近似算法属于燃耗方程求解的矩阵数值逼近类算法，具有较高的精度和效率；虽然其效率略低于切比雪夫有理近似方法，但是其求得的解精度随时间变化的稳定性更高；使用线积分有理近似方法可以展开原子核密度向量以及共轭原子核密度向量为求解时间区间上的向量函数式，进而代入积分项中求得微扰燃耗灵敏度系数，预期在达到同样计算精度的条件下该方法的效率和存储需求均优于数值积分方法。为了实现以上目的，本发明采取如下的技术方案：基于线积分有理近似方法求解微扰燃耗灵敏度系数的方法，包括如下步骤：步骤1：确定线积分有理近似的离散化参数：选取积分曲线Γ为公式(1)确定的抛物型曲线，根据燃耗矩阵特征值分布于负实轴附近的特性，确保所有特征值被该积分曲线围绕一次；z＝φ(u)＝μ(iu+1)2,u∈R公式(1)式中：z——复数变量u——实数变量φ(u)——映射实轴为复平面上抛物型积分曲线的函数，简称映射函数μ——控制曲线尺寸大小的因子选取离散积分点以及权重因子如下：zk=φ^(uk)=μ^(i·(k-n+12)h^+1)2k=1,...,n]]>式中：zk，ck——离散积分点以及权重因子——优化后的映射函数及其导函数——优化后的控制离散点选取的量i——虚数单位n——线积分有理近似阶数uk——离散积分点对应的实数μ^=π·n4t18Λ+1]]>h^=8Λ+1n]]>式中：t1——控制解稳定区间的量Λ——控制解稳定区间的量关于控制解稳定区间的量t1和Λ，其值的选取在后面的步骤3有进一步讨论；步骤2：选取有理近似阶数n为偶数，计算原子核密度向量函数以及共轭原子核密度向量函数，其中，原子核密度向量函数如下式：式中：N——原子核密度向量t——燃耗时间tf——燃耗步末时间I——单位矩阵A——燃耗方程系数矩阵，简称燃耗矩阵共轭原子核密度向量函数如下式：式中：N*——共轭原子核密度向量A*——燃耗矩阵的共轭矩阵其中共轭原子核密度向量N*(tf)依赖于响应R的定义。公式(4)可写做：式中：Pk——Re((zkI-tfA)-1N(0))，即式(4)第k个解分量的实部Qk——Im((zkI-tfA)-1N(0))，即式(4)第k个解分量的虚部xk，yk——积分点zk的实部与虚部ak，bk——权重因子ck的实部与虚部τ——t/tf，称为相对时间公式(5)可写做：式中：Pk*——Re((zkI-tfA*)-1N*(tf))，式(5)第k个解分量的实部Qk*——Im((zkI-tfA*)-1N*(tf))，式(5)第k个解分量的虚部步骤3：将原子核密度向量函数以及共轭原子核密度向量函数代入求解微扰燃耗灵敏度系数的积分公式中，进行解析积分，得到微扰燃耗灵敏度系数；微扰燃耗灵敏度系数由下面的积分公式计算得到：式中：R——响应——微扰燃耗灵敏度系数B——燃耗矩阵的扰动方向矩阵将公式(6)和公式(7)代入上式，使用解析的方式得到微扰燃耗灵敏度系数；当相对时间τ接近0时，原子核密度向量函数的精度较差；而相对时间τ接近1时，共轭原子核密度向量函数的精度较差；为了降低灵敏度系数的求解误差，采用相对时间分段的方法。以相对时间分三段为例，分别得到在相对时间区间[0.0,0.01]，[0.01,0.1]和[0.1,1.0]上的原子核密度向量函数，以及相对时间区间[0.99,1.0]，[0.9,0.99]和[0,0.9]上的共轭原子核密度向量函数。在三段相对时间区间上，控制解稳定区间的量Λ取为3，t1分别取为0.01，0.1以及1.0。公式(8)分五段积分，即[0.0,0.01]，[0.01,0.1]，[0.1,0.9]，[0.9,0.99]和[0.99,1.0]。上面的分段策略可以保证灵敏度系数计算的相对误差在千分位。对于不同的求解精度和效率要求，可以采取不同的分段数目以及相应的调整控制稳定区间的量的取值。和现有技术相比较，本发明具备如下优点：对于求解微扰燃耗灵敏度系数，本发明所提出的方法能够在保证精度需求的前提条件下，相比数值积分方法获得更高的效率和节省内存使用。以上述相对时间三段划分为例，本发明提出的方法的计算量相当于12次燃耗计算，存储量相当于300个原子核密度向量(按照有理近似阶数为24估计)。在达到相同精度的条件下，预计数值积分方法至少需要上千次燃耗计算的计算量以及相对应数目的原子核密度向量的存储量。因此本发明提出的方法可以作为微扰燃耗敏感性分析的有力工具。附图说明图1QRAM和CRAM对指数函数的逼近误差比较。具体实施方式下面结合附图对本发明做进一步详细说明：本发明使用线积分有理近似方法展开原子核密度向量函数与共轭原子核密度向量函数，以此为基础得到微扰燃耗灵敏度系数。如图1所示，线积分有理近似方法得到的解(QRAM)，其精度随时间变化的稳定性高于切比雪夫有理近似方法(CRAM)。具体包括如下步骤：步骤1：确定线积分有理近似的离散化参数。选取积分曲线Γ为公式(1)确定的抛物型曲线，根据燃耗矩阵特征值分布于负实轴附近的特性，可以确保所有特征值被该积分曲线围绕一次。因此可以应用矩阵指数基于复平面上线积分的定义式。z＝φ(u)＝μ(iu+1)2,u∈R公式(1)式中：z——复数变量u——实数变量φ(u)——映射实轴为复平面上抛物型积分曲线的函数，简称映射函数μ——控制曲线尺寸大小的因子线积分的计算借助于数值积分。根据平衡截断误差与数值积分误差的原则，选取离散积分点以及权重因子如下：zk=φ^(uk)=μ^(i·(k-n+12)h^+1)2k=1,...,n]]>式中：zk，ck——离散积分点以及权重因子——优化后的映射函数及其导函数——控制离散点选取的量，后面有表达式i——虚数单位n——线积分有理近似阶数uk——离散积分点对应的实数μ^=π·n4t18Λ+1]]>h^=8Λ+1n]]>式中：t1——控制解稳定区间的量Λ——控制解稳定性的量关于控制解稳定区间的量t1和Λ，其值的选取在后面的步骤3有进一步讨论。步骤2：计算原子核密度向量函数以及共轭原子核密度向量函数。由于采用了数值积分方法计算线积分，存在近似，所以使用约等号。其中，原子核密度向量函数如下式：式中：N——原子核密度向量t——燃耗时间tf——燃耗步末时间I——单位矩阵A——燃耗方程系数矩阵，简称燃耗矩阵有理近似阶数n一般取为偶数，此时离散积分点以及权重因子呈共轭复数对的形式出现，因此上式可以简化为：类似的，共轭原子核密度向量函数如下式：式中：N*——共轭原子核密度向量A*——共轭燃耗矩阵其中共轭原子核密度向量N*(tf)依赖于响应R的定义。通过求解一系列复数线性方程组，公式(4)可写做：式中：Pk——Re((zkI-tfA)-1N(0))，即式(4)第k个解分量的实部Qk——Im((zkI-tfA)-1N(0))，即式(4)第k个解分量的虚部xk，yk——积分点zk的实部与虚部ak，bk——权重因子ck的实部与虚部τ——t/tf，称为相对时间公式(5)可写做：式中：Pk*——Re((zkI-tfA*)-1N*(tf))，式(5)第k个解分量的实部Qk*——Im((zkI-tfA*)-1N*(tf))，式(5)第k个解分量的虚部步骤3：将原子核密度向量函数以及共轭原子核密度向量函数代入求解微扰燃耗灵敏度系数的积分公式中，进行解析积分，得到微扰燃耗灵敏度系数。微扰燃耗灵敏度系数可由下面的积分公式计算得到：式中：R——响应——微扰燃耗灵敏度系数B——燃耗矩阵的扰动矩阵将公式(6)和公式(7)代入上式，即可使用解析的方式得到微扰燃耗灵敏度系数。当相对时间τ接近0时，原子核密度向量函数的精度较差；而相对时间τ接近1时，共轭原子核密度向量函数的精度较差。为了降低灵敏度系数的求解误差，采用相对时间分段的方法。以相对时间分三段为例，分别得到在相对时间区间[0.0,0.01]，[0.01,0.1]和[0.1,1.0]上的原子核密度向量函数，以及相对时间区间[0.99,1.0]，[0.9,0.99]和[0,0.9]上的共轭原子核密度向量函数，在三段相对时间区间上，控制解稳定区间的量Λ取为3，t1分别取为0.01，0.1以及1.0。公式(8)分五段积分，即[0.0,0.01]，[0.01,0.1]，[0.1,0.9]，[0.9,0.99]和[0.99,1.0]。如下式：式中：na(τ)，nb(τ)，nc(τ)——t1分别取0.01，0.1和1.0时得到的以相对时间为因变量的原子核密度向量函数na*(τ)，nb*(τ)，nc*(τ)——t1分别取0.01，0.1和1.0时得到的以相对时间为因变量的共轭原子核密度向量函数对于式(9)的计算，以相对时间区间[0,0.01]上的积分为例：式中：B(i,j)——矩阵B的第i行第j列元素计算矩阵Ca：Ca(2i-1,2j-1)=∫00.01exc,i(1-τ)+xa,jτ(ac,icos(yc,i(1-τ))-bc,isin(yc,i(1-τ)))×(aa,jcos(ya,jτ)-ba,jsin(ya,jτ))dτ]]>Ca(2i-1,2j)=-∫00.01exc,i(1-τ)+xa,jτ(ac,icos(yc,i(1-τ))-bc,isin(yc,i(1-τ)))×(ba,jcos(ya,jτ)+aa,jsin(ya,jτ))dτ]]>Ca(2i,2j-1)=-∫00.01exc,i(1-τ)+xa,jτ(bc,icos(yc,i(1-τ))+ac,isin(yc,i(1-τ)))×(aa,jcos(ya,jτ)-ba,jsin(ya,jτ))dτ]]>式(10)可转化为：式中：M*c(i,:)——矩阵M*c的第i行Ma(j,:)——矩阵Ma的第j行上面的分段策略可以保证灵敏度系数计算的相对误差在千分位量级。对于不同的求解精度和效率要求，可以采取不同的分段数目以及相应的调整控制稳定区间的量Λ和t1的取值。当前第1页1 2 3