本发明涉及核反应堆堆芯设计和反应堆物理计算领域，具体涉及一种改进反应堆堆芯三维中子通量密度精细分布的方法。
背景技术：
：反应堆堆芯设计和运行需要准确快速地计算出反应堆及相关的设备内三维中子通量密度分布的情况。目前广泛采用的传统反应堆物理分析计算方法，随着反应堆堆芯设计越来越复杂，安全要求越来越高，逐渐不能满足工程计算要求。所谓传统反应堆物理分析计算方法又称为“两步法”，第一步是在全反射边界条件下对各种非均匀组件进行多群中子输运计算，得到组件的少群等效均匀化群常数及不连续因子等物理量；第二步，是根据前一步生成的均匀化参数，采用粗网节块方法对堆芯进行少群中子扩散或输运计算，附加功率重构等计算，便可实现对三维堆芯功率分布的模拟。“两步法”由于只考虑有限孤立组件，导致真实堆芯布置中组件之间的相互影响无法充分考虑，以及组件的燃耗历史效应等重要因素，其理论已经限制了对计算精度的追求。“一步”全堆芯非均匀中子输运计算可以从根本上解决上述问题。目前的计算条件下，直接对大型核反应堆堆芯进行三维中子输运计算还不现实。目前的二维/一维耦合方法被提出，核心思想是耦合二维、一维计算来替代三维计算。但是目前该方法中关于泄漏项的处理有些简单，需要做一些改进。技术实现要素：为了改进反应堆三维中子通量密度精细分布获取精度，本发明的目的在于提供一种改进反应堆堆芯三维中子通量密度精细分布的方法，本发明方法将二维/一维耦合计算中径向泄漏项的平近似改为正交多项式展开，轴向泄漏项用通量计算出来的分布来代替平近似，能够很好的改善反应堆堆芯三维中子通量密度的精细分布的计算精度，给反应堆堆芯设计和安全提供可靠的信息。为了实现上述目的，本发明采用如下技术方案：一种改进反应堆堆芯三维中子通量密度精细分布的方法，步骤如下：步骤1：对于所涉及的反应堆堆芯进行几何建模，划分计算区域，离散角度空间，生成特征线信息，指定各计算区域的材料，获取材料的宏观截面参数，为计算区域的中子通量密度通量、反应堆边界条件、特征值设置初始值；a)先将反应堆堆芯按照轴向划分成若干层，每一层划分成若干计算区域；b)根据计算需求对每一层计算区域的三维角度空间进行离散；c)在每一层内生成特征线信息；d)根据所涉及反应堆堆芯材料信息，读取各层中每一个计算细区的宏观截面参数；e)对于二维、一维计算区域的中子通量密度，反应堆边界条件，特征值赋初值；步骤2：更新径向泄漏项，根据步骤1中所得几何尺寸、材料信息计算一维离散纵标差分方程，求得每个栅元的每一群中子通量密度，求得轴向泄漏项；在二维/一维耦合算法中，径向泄漏项是由二维求解器计算提供的，因此在求解径向一维方程时，认为该项是已知源项；一般计算中，每一层的宽度为10cm-20cm，二维提供的径向泄漏项在每一层内是常数，引入这样一个平泄漏近似对于计算精度是有影响的；所以给一维计算提供的径向泄漏项是对二维提供的平泄漏做了勒让德多项式多项式展开，得到径向泄漏沿轴向的分布，对于最后的计算结果有明显改善；步骤3：更新轴向泄漏项，根据步骤1中所得几何尺寸、材料信息、特征线信息对轴向的每一层进行二维输运计算，迭代求解每一层的二维中子通量密度，计算栅元均匀化截面，计算径向泄漏项；在计算过程中，利用平源区通量计算得到一个栅元内轴向泄漏项的空间分布；在二维/一维耦合算法中，轴向泄漏项是由一维轴向求解计算提供的，因此在求解径向二维方程时，认为该项是已知源项；步骤2中一维计算的区域选择的是一个栅元，所以上述轴向泄漏项在一个栅元中是一个平分布；真实情况是轴向泄露在一个栅元中也存在分布；利用上一个迭代中二维计算得到的通量分布计算得到栅元内轴向泄漏项的空间分布，对于整个二维一维耦合计算有精度上的改善；用平源区上下表面的标通量之差计算得到轴向泄漏项在栅元内的分布；平源区的上下表面的标通量则是通过将通量沿轴向展开得到的；同样采用二阶的勒让德多项式，计算得到上下表面的标通量；最终，栅元内每个平源区的差值组成的分布乘以栅元的轴向泄漏项得到栅元内轴向泄漏项的差值；从平分布近似到有分布的精确描述，提高了计算的精度；步骤4：判断特征值和三维中子通量密度是否收敛，如果不收敛，则转到步骤2继续迭代，并且计算中采用上一次迭代中求得的中子通量密度和特征值，直至中子通量密度和特征值收敛，就能够得到三维中子通量密度精细分布。与现有技术相比，本发明有如下突出优点：1.一维计算中对径向泄漏项进行多项式展开，使用带分布的径向泄漏项，比平泄露近似精度高。2.二维计算中求得轴向泄漏项的二维分布，使得计算更接近真实情况。附图说明图1为二维/一维耦合计算流程图。图2为二维/一维堆芯划分示意图。图3为反应堆堆芯几何布置示意图。图4为计算区域角度空间离散示意图。图5为计算区域特定方向特征线示意图。图6为径向泄漏项多项式展开示意图。图7为轴向泄漏项多项式展开示意图。图8为三维非均匀压水堆堆芯相对棒功率分布图。具体实施方式本发明将二维/一维耦合计算中径向泄漏项的平近似改为正交多项式展开，轴向泄漏项用通量计算出来的分布来代替平近似，它能够很好的改善反应堆堆芯三维中子通量密度的精细分布的计算精度，给反应堆堆芯设计和安全提供可靠的信息，具体实施方式如下所示。图1所示为二维/一维耦合方法的总体流程图。步骤1：对于所涉及的反应堆堆芯进行几何建模，划分计算区域，离散角度空间，生成特征线信息，指定各计算区域的材料，获取材料的宏观截面参数，为计算区域的中子通量密度通量、反应堆边界条件、特征值设置初始值，图2为二维/一维堆芯划分示意图；f)先将反应堆堆芯按照轴向划分成若干层，每一层划分成若干计算区域，图3为压水堆堆芯几何布置示意图；g)根据计算需求对每一层计算区域的三维角度空间进行离散，如图4所示的角度空间离散，图中用了9个离散方向等效八分之一角度空间，θ代表极角，代表辐角；h)在每一层内生成特征线信息，如图5所示的某区域中某一个方向的所有特征线，图中实线为特征线，虚线为两条相邻特征线的间距中心线，是出射中子角通量，是入射中子角通量密度，si,k是特征线的长度，δAk是特征线线宽；i)根据所涉及反应堆堆芯材料信息，读取各层中每一个计算细区的宏观截面参数；j)对于二维、一维计算区域的中子通量密度，反应堆边界条件，特征值赋初值；步骤2：更新径向泄漏项，根据步骤1中所得几何尺寸、材料信息计算一维离散纵标差分方程，求得每个栅元的每一群中子通量密度，求得轴向泄漏项。在二维/一维耦合算法中，径向泄漏项是由二维求解器计算提供的，因此在求解径向一维方程时，认为该项是已知源项。一般计算中，每一层的宽度在10cm-20cm左右，二维提供的径向泄漏项在每一层内是常数，引入这样一个平泄漏近似对于计算精度是有影响的。所以本发明中在给一维计算提供的径向泄漏项是对二维提供的平泄漏做了勒让德多项式多项式展开，得到径向泄漏沿轴向的分布，对于最后的计算结果有明显改善；步骤3：更新轴向泄漏项，根据步骤1中所得几何尺寸、材料信息、特征线信息对轴向的每一层进行二维输运计算，迭代求解每一层的二维中子通量密度，计算栅元均匀化截面，计算径向泄漏项；在计算过程中，利用平源区通量计算得到一个栅元内轴向泄漏项的空间分布；在二维/一维耦合算法中，轴向泄漏项是由一维轴向求解计算提供的，因此在求解径向二维方程时，认为该项是已知源项。步骤2中一维计算的区域选择的是一个栅元，所以上述轴向泄漏项在一个栅元中是一个平分布；真实情况是轴向泄露在一个栅元中也存在分布。本发明利用上一个迭代中二维计算得到的通量分布计算得到栅元内轴向泄漏项的空间分布，对于整个二维一维耦合计算有精度上的改善。本发明用平源区上下表面的标通量之差计算得到轴向泄漏项在栅元内的分布。平源区的上下表面的标通量则是通过将通量沿轴向展开得到的。同样采用二阶的勒让德多项式，计算得到上下表面的标通量。最终，栅元内每个平源区的差值组成的分布乘以栅元的轴向泄漏项得到栅元内轴向泄漏项的差值。从平分布近似到有分布的精确描述，提高了计算的精度；步骤4：判断特征值和三维中子通量密度是否收敛，如果不收敛，则转到步骤2继续迭代，并且计算中采用上一次迭代中求得的中子通量密度和特征值，直至中子通量密度和特征值收敛，就可以得到三维中子通量密度精细分布。步骤1中二维/一维耦合算法中的几何描述、计算区域划分包括轴向分层、径向划分子区、子区内划分计算区域、所有计算区域的角度空间离散、子区内特征线生成均需要根据不同的反应堆堆芯和不同的计算条件选择合适的方案。一般压水堆堆芯轴向分为10到15层，径向每一层分为9到25个子区，每个子区中几百个计算区域不等，通常使用48到80个角度方向去离散角度空间，角度空间的离散方案确定了特征线的方向，特征线线宽一般采用0.001cm到0.05cm。步骤2中求解离散纵标差分方程，径向泄漏项的展开，推导过程如下所示。三维直角坐标系中，稳态多群中子输运方程如公式(1)所示，Ω·▽ψg(r,Ω)+Σt,g(r)ψg(r,Ω)＝Qg(r,Ω),g＝1,...,G公式(1)式中：ψg——第g能群中子角通量密度g——能群标识；G——能群总数；Q——中子输运方程源项。以特征值问题为例，第g能群源项的具体形式为：式中：φg——第g能群中子通量密度——第g能群的裂变源项；——第g能群的散射源项。经过角度离散后，在第g群、第m离散方向上的形式如公式(2)所示。根据反应堆堆芯进行的几何划分，对每个栅元进行积分，并除以栅元面积，最终整理得：ξmΔx[ψg,m,x+p(z)-ψg,m,x-p(z)]+ηmΔy[ψg,m,y+p(z)-ψg,m,y-p(z)]+μmdψg,mp(z)dz+Σt,g,p(z)ψg,mp(z)=Qgp(z)]]>在该方程中，径向泄漏项定义为：TLg,m,pRadial(z)=ξmΔx[ψg,m,x+p(z)-ψg,m,x-p(z)]+ηmΔy[ψg,m,y+p(z)-ψg,m,y-p(z)]]]>其物理意义是第k层沿轴向在离散方向m上的泄漏。在二维/一维耦合算法中，径向泄漏项是由二维求解器计算提供的，因此在求解径向一维方程时，认为该项是已知源项。一般计算中，每一层的宽度在10cm-20cm左右，二维提供的径向泄漏项在每一层内是常数，引入这样一个平泄漏近似对于计算精度是有影响的。所以本发明中在给一维计算提供的径向泄漏项是对二维提供的平泄露做了多项式展开，得到径向泄漏沿轴向的分布，对于最后的计算结果有明显改善。由于径向泄漏项在每一层是平分布，为了得到径向泄漏项精细的轴向分布，本发明将径向泄漏项沿轴向展开，采用勒让德多项式，表达式如下所示。Ln(ξ)=Σi=02linPi(ξ)]]>whereP0＝1,P1＝ξ,P2＝3ξ2-1/4l0k=TLg,kRadial]]>l1k=Δzkdzk-1(TLg,k+1Radial-TLg,kRadial)+zk+1(TLg,kRadial-TLg,k-1Radial)]]>l2k=(Δzk)2d(Δzk-1+Δzk)(TLg,k+1Radial-TLg,kRadial)-(Δzk+1+Δzk)(TLg,kRadial-TLg,k-1Radial)]]>wherezi＝(Δzk+Δzi)(Δzk+2Δzi)d＝(Δzk-1+Δzk)(Δzk+Δzk+1)(Δzk-1+Δzk+Δzk+1)其中，L表示某个方向的泄漏项，经验表明采取勒让德二阶近似就可以取得很好的结果；0阶、1阶、2阶的系数分别如上述所示。得到这三个系数后，就可以轻松计算得到每个细网上的径向泄漏值。图6为径向泄漏项展开前后示意图，展开之前在每层内径向泄漏项是平分布，展开之后可以计算得到每个平源区网格上的径向泄漏项。一维Sn差分方程求解步骤：1)用勒让德多项式沿轴向展开径向泄漏项，得到每个差分网格的径向泄漏项；2)根据上次计算得到的Keff和通量更新右端源项；3)根据差分格式一次计算每个网格的、每一个方向、每一个能群的角通量；4)判断步骤3)计算得到的通量是否收敛，不收敛则继续进行2)-3)的迭代计算，若通量收敛则结束计算。收敛准则如下：||φ‾gk-φ‾gk-1φ‾gk||≤ϵ1]]>ε1为很小的正数，也就是收敛条件。步骤3中根据步骤1中反应堆堆芯进行的几何划分，对轴向的每一层进行积分，并处以网格宽度，最终整理得：在该方程中，定义轴向泄漏项为：其物理意义是第k层沿轴向在离散方向m上的泄漏。在二维/一维耦合算法中，轴向泄漏项是由一维轴向求解计算提供的，因此在求解径向二维方程时，认为该项是已知源项。步骤2中一维计算的区域选择的是一个栅元，所以上述轴向泄漏项在一个栅元中是一个平分布。真实情况是轴向泄露在一个栅元中也存在分布。本发明利用上一个迭代中二维计算得到的通量分布计算得到栅元内轴向泄漏项的空间分布，对于整个二维一维耦合计算有精度上的改善。TLg,m,kAxial(x,y)=μmΔzk[ψg,m,k+1/2(x,y)-ψg,m,k-1/2(x,y)]]]>TL_shapekAxial(x,y)=μmΔzk[φk+1/2(x,y)-φk-1/2(x,y)]]]>φ(x,y,ξ)=Σi=02linPi(ξ)]]>其中,P0＝1,P1＝ξ,P2＝3ξ2-1/4φ0k=φk]]>φ1k=Δzkdzk-1(φk+1-φk)+zk+1(φk-φk-1)]]>φ2k=(Δzk)2d(Δzk-1+Δzk)(φk+1-φk)-(Δzk+1+Δzk)(φk-φk-1)]]>其中,zi＝(Δzk+Δzi)(Δzk+2Δzi)d＝(Δzk-1+Δzk)(Δzk+Δzk+1)(Δzk-1+Δzk+Δzk+1)φk+1/2(x,y)＝φ(x,y,1/2)φk-1/2(x,y)＝φ(x,y,-1/2)上述表达式中，本发明用平源区上下表面的标通量之差计算得到轴向泄漏项在栅元内的分布。平源区的上下表面的标通量则是通过将通量沿轴向展开得到的。同样采用二阶的勒让德多项式，计算得到上下表面的标通量。最终，栅元内每个平源区的差值组成的分布乘以栅元的轴向泄漏项得到栅元内轴向泄漏项的差值。从平分布近似到有分布的精确描述，提高了计算的精度。图7为轴向泄漏项展开前后示意图，重构之前在每个栅元内轴向泄漏项是平分布，展开之后可以计算得到每个平源区网格上的轴向泄漏项。MOC方法是将公式(2)转换成沿着某条特征线的微分方程，如下所示：认为Σt,g,Qg在某个区域内是常数，则公式(6)是一个常微分方程，可以有解析解。步骤1中二维的计算区域在合适的范围内，MOC方法可以得到非常精确的输运解。沿着某一段特征线的平均中子角通量密度表达式如下所示：对同一方向所有特征线的平均中子角通量密度进行体积加权平均得到该计算区域这个方向的平均中子角通量密度，对该计算区域所有方向的中子角通量密度进行加权平均就可以得到该计算区域的平均中子通量密度。由于公式(7)的右端含有源项，该源项需要通过中子通量密度进行求解，所以一般采用迭代计算求解MOC方程。二维MOC方程求解流程如下：1)更新轴向泄漏项，该泄漏项由一维计算提供，然后根据上述重构方式，得到轴向泄漏项在一个栅元内的空间分布；2)根据上一步计算得到的Keff、入射角通量和标通量更新右端源项；3)根据MOC方程解的表达形式迭代计算每个平源区的标通量；4)判断平源区标通量是否收敛，若不收敛，继续迭代求解，若收敛更新径向泄漏项和栅元均匀化截面；判断所有能群是否都已经计算，对各个能群进行迭代，直到所有区域通量收敛。步骤4判断特征值和三维通量是否收敛，如果不收敛，则转到步骤3继续迭代，并且计算中采用上一次迭代中步骤2-步骤3求得的中子通量密度和特征值，反复迭代直至中子通量密度和特征值收敛，就可以得到三维中子通量精细分布。所述收敛条件为：||φgi,(l)-φgi,(l-1)φgi,(l)||≤ϵ1,||keff(l)-keff(l-1)keff(l)||≤ϵ2]]>其中：为第l次迭代第i区第g群通量为第l次迭代特征值ε1，ε2为两个很小的正数，也就是收敛条件。大量数值验证结果显示，本发明具有可靠的精度、很好的效率、对于平泄漏近似的改善非常有效，适应工程实际中的计算要求。图8是典型的压水堆堆芯三维功率分布，其中外围黑色部分是反射层没有功率，内部的控制棒区域也没有功率。当前第1页1 2 3