一种加速在线生成核反应堆堆芯特征线的方法与流程

本发明涉及核反应堆堆芯数值计算领域，具体涉及一种加速在线生成核反应堆堆芯特征线的方法。

背景技术

中子输运方程又称玻尔兹曼输运方程，描述的是处于运动状态的中子的守恒关系，如式(1)所示：

式中：

ψ——中子角通量密度；

r——中子的空间位置向量；

ω——中子运动的角度向量；

e——中子的能量；

t——时间；

v——中子运动的速度；

——数学算符，其中i,j,k分别为x,y,z方向的方向向量；

∑t——反应总截面；

q——中子源项；

当讨论稳态问题时，认为中子通量密度不随时间变化，即：

且ψ(r,ω,e,t)＝ψ(r,ω,e)，此时方程变成：

在特征线上，式(2)可写成：

式中：

s——特征线坐标；

对能量进行离散可得特征线上的多群中子输运方程：

将求解区域划分为许多源区i，用离散后的角度m代替ω。在源区内式(4)可以写成：

式中：

k——方向为m穿过源区的特征线编号；

对源项采取平源近似，有：

qi,m(s)＝qi,m(6)

联立式(5)和式(6)，可以解得特征线上的中子角通量密度分布：

在实际特征线计算中，令s＝sk即可得到第k条特征线的出射中子角通量：

由式(8)可知，欲求解特征线方程首先应该知道特征线信息，包括特征线段长度sk和特征线段对应的截面信息。为获得这些信息，通常有两种方法可以选择，一种是生成模块化特征线并存储各模块的特征线参数；另一种是在线生成特征线并计算特征线参数。当采用第一种方法时，对于模块几何结构比较复杂、几何尺寸比较大、或者模块种类比较多的情况，都会大大增加计算机的存储负担。但是对于几何较为复杂的反应堆堆芯，我们采用在线生成特征线的方法，其核心思想是在特征线扫描的同时针对局部几何获取特征线参数。但是传统的在线生成特征线方法计算效率非常低，原因是在特征线扫描之前不存储任何关于特征线的参数，所有参数都必须在线计算，使用过后随机释放内存。这大大地限制了在线生成特征线方法在实际工程中的应用。

由前文可知，传统的模块化特征线方法需要在进行特征线计算之前存储模块内所有特征线的信息，这导致了它对于几何的适应性较差。在线生成特征线方法可以有效地拓展特征线方法的几何适应性，但是由于传统的在线生成特征线方法不存储特征线信息，每条长特征线上的计算参数都在扫描时在线计算生成，并在扫描结束后抛弃以释放计算机内存。在真实的核反应堆堆芯中，往往有成千上万个基本几何单元(如轻水堆中的栅元结构、球床堆中的燃料球结构等)排列在其中，每条长特征线往往会穿过几十至上百个这种基本几何单元。直接将在线生成特征线方法应用到这类几何会带来大量的重复计算，会大大降低程序的计算效率。

技术实现要素：

为了克服上述现有技术存在的问题，本发明的目的在于提供一种加速在线生成核反应堆堆芯特征线的方法，该方法适用于存在重复基本单元的核反应堆堆芯几何结构，考虑到现有反应堆堆型以及当前所规划的下一代反应堆堆型，或者包含燃料栅元，或者包含燃料球，而燃料栅元和燃料球在堆芯中都是大量重复存在的基本几何单元；因此本发明方法对于待模拟的核反应堆堆型没有要求；为了方便，以包含3个平源区的二维压水堆燃料栅元为例说明本发明的技术方案，但是本发明同样适合其他包含多个平源区的二维或三维几何。

为了达到上述目的，本发明采用如下技术方案：

一种加速在线生成核反应堆堆芯特征线的方法，包括如下步骤：

步骤1：典型压水堆燃料栅元由内到外分别为燃料区、气隙区和慢化剂区；描述一个典型压水堆燃料栅元的几何共需要8个面，分别为：上边界面、下边界面、左边界面、右边界面、前边界面、后边界面、气隙内表面和气隙外表面；其中气隙内表面和气隙外表面为圆柱面，其余面均为平面；若考虑二维压水堆燃料栅元，认为栅元几何在z方向无限，则上边界面和下边界面消失，同时，三维燃料栅元中的“面”在二维几何中变成了“线”；所以二维压水堆燃料栅元几何由6条线描述，分别为前边界线、后边界线、左边界线、右边界线以及气隙的内边界线、外边界线；其中气隙内边界线和外边界线为同心圆，其余边界线均为直线段。

为了叙述方便，给6条线分别编号如下：

其中左边界线为1号边界，前边界线为2号边界，右边界线为3号边界，后边界线为4号边界，气隙内边界线为5号边界，气隙外边界线为6号边界；

若以气隙内边界线和外边界线的圆心为坐标原点建立二维直角坐标系，根据实际压水堆燃料栅元的尺寸，写出各条边界线的方程如下：

1号边界：x＝-0.63

2号边界：y＝-0.63

3号边界：x＝0.63

4号边界：y＝0.63

5号边界：r＝0.41

6号边界：r＝0.47；

步骤2：在几何描述完成之后，对目标几何应用特征线方法进行计算之前，应该产生特征线；特征线是穿过目标几何的一簇平行的等距的直线，每条直线都被目标几何截成若干段，记作n段(n≥1)，这些特征线段是进行特征线计算的基本单位；特征线段的空间位置由角度和起始点坐标唯一确定，在二维几何中，角度被定义为特征线与x正方向的夹角，特征线的角度由用户指定，此处以说明；特征线段的起始点坐标为(x,y)的形式，对于同一条特征线，第n+1段特征线的起始点即为第n段特征线的终点；由于存在大量的特征线段，而每一条特征线的起始点各不相同，所以并不提前存储所有特征线段的起始点坐标，仅存储每条特征线第一段的起始点坐标。由经过数学换算，将穿过压水堆燃料栅元的平行特征线段统一表示为：

y＝0.5774x+b(9)

其中b特征线在边界上的截距，由特征线段的起始点确定；

步骤3：对于角度为30°的特征线，根据起始点的不同，从1号边界出发的特征线的终点有两种可能：到达4号边界或者到达6号边界，最终到达哪个边界取决于特征线在1号边界上的起始点坐标与1号边界上分界点之间的关系；从2号边界出发的特征线的终点有两种可能：到达6号边界或者到达3号边界，最终到达哪个边界取决于特征线在2号边界上的起始点坐标与2号边界上分界点之间的关系；从3号边界和4号边界出发的特征线将离开此基本几何单元，故不需考虑；从5号边界出发的特征线的终点有两种可能：到达5号边界或者到达6号边界，最终到达哪个边界取决于特征线在5号边界上的起始点坐标与5号边界上分界点之间的关系；从6号边界出发的特征线的终点有四种可能：到达6号边界、到达5号边界、到达4号边界以及到达3号边界，最终到达哪个边界取决于特征线在6号边界上的起始点坐标与6号边界上分界点之间的关系；具体对应关系见步骤4说明；

步骤4：求解各个边界上的分界点，所谓分界点，是指在特征线上从该点的两侧出发，将到达不同的边界；具体如下：

1号边界：联立式(10)和式(11)

(y-y0)＝0.5774(x+0.63)(10)

x²+y²＝0.47²(11)

考虑直线和圆相切的情况，解得：

y0＝0.1689或y0＝-0.9166(舍去)

当-0.63≤y<0.1689时，特征线将和6号边界相交；当0.1689≤y<0.63时，特征线将和4号边界相交；

2号边界：联立式(12)和式(13)

(y+0.63)＝0.5774(x-x0)(12)

x²+y²＝0.47²(13)

考虑直线和圆相切的情况，解得：

x0＝-0.1512或x0＝-2.0309(舍去)

当-0.63≤x<-0.1512时，特征线将和6号边界相交；当-0.1512≤y<0.63时，特征线将和3号边界相交；

5号边界：过原点作与特征线方向垂直的辅助线a，由于以特征线角度说明，所以辅助线a的角度为120°，辅助线a与5号边界有2个交点，分别在第二象限与第四象限，记第二象限的交点为a，第四象限的交点为b；将从a点逆时针到达b点所表示的圆弧记为将从b点逆时针到达a点所表示的圆弧记为从5号边界出发的特征线，根据其出发点的位置不同，将到达不同的目标边界；具体来说，从出发的特征线，将到达5号边界；从出发的特征线，将到达6号面边界；已知5号边界是半径为0.41的圆，所以点a和点b的坐标用极坐标表示为：a:b:

点a和点b确定之后，将各区间表示为：

对于从5号边界出发的特征线，设其起始点坐标为：(x,y)，则：当时，特征线将和5号边界相交；当或时，特征线将和6号边界相交；

6号边界：过原点作与特征线方向垂直的辅助线a，由于以特征线角度说明，所以辅助线a的角度为120°，辅助线a与6号边界有2个交点，分别在第二象限与第四象限，记第二象限的交点为a，第四象限的交点为d；将从a点逆时针到达d点所表示的圆弧记为将从d点逆时针到达a点所表示的圆弧记为直线a与5号边界也有2个交点，作经过这2个交点的5号边界的切线b和c，切线b和c与的交点记为b和c。过3号边界和4号边界的交点作平行于特征线方向的辅助线d，直线d与的交点记为e；这样，6号边界被分为等5段；已知6号边界是半径为0.47的圆，利用各条辅助线的方程，通过几何运算获得点a～点e的坐标，用极坐标表示为：a:b:(0.47,4.1762)，c:d:e:(0.47,1.0364)。

点a至点e确定之后，将各区间表示为：

对于从6号边界出发的特征线，设其起始点坐标为：(x,y)，则：当或时，特征线将和6号边界相交；当时，特征线将和5号边界相交；当或时，特征线将和3号边界相交；当时，特征线将和5号边界相交；

步骤5：存储上述各个边界的坐标区间与特征线目标边界的对应关系；

步骤6：在特征线扫描的过程中，根据所存储的对应关系，由特征线在某个边界上的起点坐标，迅速检索到其目标边界；即能快速地计算下一段特征线段的长度，将此参数提供给式(8)即能完成此段特征线计算；

下面针对典型压水堆栅元，假设特征线辐角为给定特征线的起始点坐标，描述此特征线的扫描以说明此过程：

假设特征线从1号边界的点a出发，给定a的坐标为：(-0.63,-0.31)，根据所存储的1号边界上的对应关系：-0.63≤y<0.1689，可知特征线的目标边界是6号边界；通过特征线与6号边界求交点得到点b的坐标(-0.4282,-0.1935)以及ab的长度0.233，与ab段的材料参数一起，即能够进行后续的特征线方程计算；

由点b的坐标，根据所存储的6号边界上的对应关系：可知特征线的目标边界是5号边界；通过特征线与5号边界求交点得到点c的坐标(-0.376,-0.1634)以及bc的长度0.0601，与bc段的材料参数一起，即能够进行后续的特征线方程计算；

由点c的坐标，根据所存储的5号边界上的对应关系：可知特征线的目标边界是5号边界；通过特征线与5号边界求交点得到点d的坐标(0.3296,0.244)以及cd的长度0.8148，与cd段的材料参数一起，即能够进行后续的特征线方程计算；

由点d的坐标，根据所存储的5号边界上的对应关系：或可知特征线的目标边界是6号边界；通过特征线与6号边界求交点能够得到点e的坐标(0.3818,0.2741)以及de的长度0.0601，与de段的材料参数一起，即能够进行后续的特征线方程计算；

由点e的坐标，根据所存储的6号边界上的对应关系：或可知特征线的目标边界是3号边界；通过特征线与3号边界求交点得到点f的坐标(1.26,0.4174)以及ef的长度0.8898，与ef段的材料参数一起，即能够进行后续的特征线方程计算；

通过以上过程，即完成此燃料栅元几何中的特征线扫描计算。

和现有技术相比较，本发明具备如下优点：

采取传统的在线生成特征线方法对全堆芯生成特征线时，一条长特征线被复杂的堆芯几何分割成了大量的特征线段，获取每条特征线段的长度信息的方式为：已知特征线起始点，全局检索特征线的目标边界(特征线的起始点坐标和方向求解与每个边界的距离，距离最短者为目标边界)。由于压水堆堆芯中存在成千上万个边界面，所以此过程将耗费海量的计算代价。

本发明考虑到在真实的压水堆堆芯中，充满了大量的燃料栅元结构，这些燃料栅元结构大量重复出现；将这些大量重复出现的几何结构视为典型几何单元，存储典型几何单元的特征信息，在获取目标边界的时候，不需要全局搜索目标边界，而是根据所存储的对应关系直接获取。这样使得生成特征线这一过程的计算代价节省了成千上万倍。

附图说明

图1为5号边界上的分界点。

图2为6号边界上的分界点。

图3为特征线扫描过程。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步详细说明：

本实施例一种加速在线生成核反应堆堆芯特征线的方法，包括如下步骤：

步骤1：典型压水堆燃料栅元由内到外分别为燃料区、气隙区和慢化剂区。描述一个典型压水堆燃料栅元的几何共需要8个面，分别为：上边界面、下边界面、左边界面、右边界面、前边界面、后边界面、气隙内表面和气隙外表面。其中气隙内、外表面为圆柱面，其余面均为平面。若考虑二维压水堆燃料栅元，认为栅元几何在z方向无限，则上、下边界面消失。同时，三维燃料栅元中的“面”在二维几何中变成了“线”。所以可以说：二维压水堆燃料栅元几何可由6条线描述。分别为前、后、左、右边界线以及气隙的内、外边界线。其中气隙内、外边界线为同心圆，其余边界线均为直线段。

为了叙述方便，不妨给6条线分别编号。

其中左边界线为1号边界，前边界线为2号边界，右边界线为3号边界，后边界线为4号边界，气隙内边界线为5号边界，气隙外边界线为6号边界。

若以气隙内、外边界线的圆心为坐标原点建立二维直角坐标系，根据实际压水堆燃料栅元的尺寸，可以写出各条边界线的方程如下。

1号边界：x＝-0.63

2号边界：y＝-0.63

3号边界：x＝0.63

4号边界：y＝0.63

5号边界：r＝0.41

6号边界：r＝0.47

步骤2：在几何描述完成之后，对目标几何应用特征线方法进行计算之前，应该产生特征线。特征线是穿过目标几何的一簇平行的等距的直线，每条直线都被目标几何截成n段(n≥1)。这些特征线段是进行特征线计算的基本单位。特征线段的空间位置由角度和起始点坐标唯一确定。在二维几何中，角度被定义为特征线与x正方向的夹角。特征线的角度由用户指定，此处以为例。特征线段的起始点坐标为(x,y)的形式，对于同一条特征线，第n+1段特征线的起始点即为第n段特征线的终点。由于存在大量的特征线段，而每一条特征线的起始点各不相同，所以并不提前存储所有特征线段的起始点坐标。仅存储每条特征线第一段的起始点坐标。由经过简单的数学换算，可以将穿过压水堆燃料栅元的平行特征线段统一表示为：

y＝0.5774x+b(9)

其中b由特征线段的起始点确定。

步骤3：对于角度为30°的特征线，根据起始点的不同，从1号边界出发的特征线的终点有两种可能：到达4号边界或者到达6号边界，最终到达哪个边界取决于特征线在1号边界上的起始点坐标与1号边界上分界点之间的关系；从2号边界出发的特征线的终点有两种可能：到达6号边界或者到达3号边界，最终到达哪个边界取决于特征线在2号边界上的起始点坐标与2号边界上分界点之间的关系；从3号边界和4号边界出发的特征线将离开此基本几何单元，故不需考虑；从5号边界出发的特征线的终点有两种可能：到达5号边界或者到达6号边界，最终到达哪个边界取决于特征线在5号边界上的起始点坐标与5号边界上分界点之间的关系；从6号边界出发的特征线的终点有四种可能：到达6号边界、到达5号边界、到达4号边界以及到达3号边界，最终到达哪个边界取决于特征线在6号边界上的起始点坐标与6号边界上分界点之间的关系。具体对应关系见步骤4说明；

步骤4：求解各个边界上的分界点，所谓分界点，是指在特征线从该点的两侧出发，将到达不同的边界；具体如下：

1号边界：联立式(10)和式(11)

(y-y0)＝0.5774(x+0.63)(10)

x²+y²＝0.47²(11)

考虑直线和圆相切的情况，解得：

y0＝0.1689或y0＝-0.9166(舍去)

当-0.63≤y<0.1689时，特征线将和6号边界相交；当0.1689≤y<0.63时，特征线将和4号边界相交。

2号边界：联立式(12)和式(13)

(y+0.63)＝0.5774(x-x0)(12)

x²+y²＝0.47²(13)

考虑直线和圆相切的情况，可以解得：

x0＝-0.1512或x0＝-2.0309(舍去)

当-0.63≤x<-0.1512时，特征线将和6号边界相交；当-0.1512≤y<0.63时，特征线将和3号边界相交。

5号边界：如图1所示，过原点作与特征线方向垂直的辅助线a，由于以特征线角度为例，所以辅助线a的角度为120°，辅助线a与5号边界有2个交点，分别在第二象限与第四象限，记第二象限的交点为a，第四象限的交点为b。将从a点逆时针到达b点所表示的圆弧记为将从b点逆时针到达a点所表示的圆弧记为从5号边界出发的特征线，根据其出发点的位置不同，将到达不同的目标边界。具体来说，从出发的特征线，将到达5号边界；从出发的特征线，将到达6号面边界。已知5号边界是半径为0.41的圆，所以点a和点b的坐标可以用极坐标表示为：a:b:

点a和点b确定之后，可以将各区间表示为：

对于从5号边界出发的特征线，设其起始点坐标为：(x,y)，则：当时，特征线将和5号边界相交；当或时，特征线将和6号边界相交。

6号边界：如图2所示，过原点作与特征线方向垂直的辅助线a，由于以特征线角度为例，所以辅助线a的角度为120°，辅助线a与6号边界有2个交点，分别在第二象限与第四象限，记第二象限的交点为a，第四象限的交点为d。将从a点逆时针到达d点所表示的圆弧记为将从d点逆时针到达a点所表示的圆弧记为直线a与5号边界也有2个交点，作经过这2个交点的5号边界的切线b和c，切线b和c与的交点记为b和c。过3号边界和4号边界的交点作平行于特征线方向的辅助线d，直线d与的交点记为e。这样，6号边界被分为等5段。已知6号边界是半径为0.47的圆，利用各条辅助线的方程，通过简单的几何运算可以获得点a～点e的坐标，用极坐标可以表示为：a:b:(0.47,4.1762)，c:d:e:(0.47,1.0364)。

点a至点e确定之后，可以将各区间表示为：

对于从6号边界出发的特征线，设其起始点坐标为：(x,y)，则：当或时，特征线将和6号边界相交；当时，特征线将和5号边界相交；当或时，特征线将和3号边界相交；当时，特征线将和5号边界相交。

步骤5：存储上述各个边界的坐标区间与特征线目标边界的对应关系；

步骤6：在特征线扫描的过程中，根据所存储的对应关系，由特征线在某个边界上的起点坐标，迅速检索到其目标边界。即可快速地计算下一段特征线段的长度，将此参数将提供给式(8)即可完成此段特征线计算。下面针对典型压水堆栅元，假设特征线辐角为给定特征线的起始点坐标，描述此特征线的扫描以说明此过程。

如图3所示，假设特征线从1号边界的点a出发，给定a的坐标为：(-0.63,-0.31)，根据所存储的1号边界上的对应关系：-0.63≤y<0.1689，可知特征线的目标边界是6号边界。通过特征线与6号边界求交点可以得到点b的坐标(-0.4282,-0.1935)以及ab的长度0.233，与ab段的材料参数一起，即可进行ab段的特征线方程计算。

由点b的坐标，根据所存储的6号边界上的对应关系：可知特征线的目标边界是5号边界。通过特征线与5号边界求交点可以得到点c的坐标(-0.376,-0.1634)以及bc的长度0.0601，与bc段的材料参数一起，即可进行bc段的特征线方程计算。

由点c的坐标，根据所存储的5号边界上的对应关系：可知特征线的目标边界是5号边界。通过特征线与5号边界求交点可以得到点d的坐标(0.3296,0.244)以及cd的长度0.8148，与cd段的材料参数一起，即可进行cd段的特征线方程计算。

由点d的坐标，根据所存储的5号边界上的对应关系：或可知特征线的目标边界是6号边界。通过特征线与6号边界求交点可以得到点e的坐标(0.3818,0.2741)以及de的长度0.0601，与de段的材料参数一起，即可进行de段的特征线方程计算。

由点e的坐标，根据所存储的6号边界上的对应关系：或可知特征线的目标边界是3号边界。通过特征线与3号边界求交点可以得到点f的坐标(1.26,0.4174)以及ef的长度0.8898，与ef段的材料参数一起，即可进行ef段的特征线方程计算。

通过以上过程，即可完成此燃料栅元几何中的特征线扫描计算。