一种应用于屏蔽计算曲面衔接中的S

一种应用于屏蔽计算曲面衔接中的s
n-mc耦合计算方法
技术领域
1.本发明涉及核反应堆屏蔽计算领域，是一种针对曲面几何（包含球面和椭球面）衔接下的离散纵标-蒙特卡罗耦合计算方法。


背景技术：

2.为了结合离散纵标方法（以下简称sn方法）求解速度快和蒙特卡罗方法（以下简称mc方法）几何建模精确的优点，上世纪70年代美国橡树岭国家实验室的emmet等提出了二维s
n-mc耦合方法，同时基于二维sn差分程序dot和三维mc程序morse开发了中间接口domino。该方法的核心思想是利用sn计算提供的角注量率转化为关于能量、空间和角度的累积分布函数作为mc续算的源项输入。在该方法的后续发展中，逐渐实现了三维sn程序和mc程序的在坐标系下平面、圆柱面的耦合计算。至2011年，美国田纳西大学的kulesza开发了disco软件包，可以实现二维三维下两类坐标系不同几何的耦合计算，代表该方法在常见几何下发展已比较成熟。由于反应堆几何结构通常以圆柱、棱柱等为主，常见几何的平面、柱面面源能够满足大部分情形下的耦合计算需求，但对某些特殊构件中的球形或椭球形部分，需要使用球面和椭球面作为续算面源。同时考虑到sn方法中坐标系网格无法准确的描述上述两种曲面几何，因此需要研制一种在坐标系下近似求解曲面面源的方法。


技术实现要素：

3.为了克服上述现有技术存在的问题，本发明的目的在于提供一种应用于屏蔽计算曲面衔接中的s
n-mc耦合计算方法，已球面为例，本发明在处理曲面的耦合计算问题时，在sn完成比面源区域略大部分的输运计算后，根据sn计算划分的z方向网格间隔将球面面源分为多层的圆形环带，计算圆环投影与sn划分xy平面网格的交点，根据交点位置确定球面穿过的sn网格，将这些网格中心点的角注量率近似视为穿过该网格的弧段上的角注量率，同时根据交点和z方向高度确定该弧段出射法向量与xyz坐标系的夹角余弦，从而得到不同位置出射的面源法向量。利用角注量率可以求得粒子类型(p)、面源(s)、能群(g)、环带分层(k)、网格(q)、方向（m）上的粒子出射流，再将出射流转化为累积分布函数（cdf），作为mc程序的源项概率分布输入，进行源抽样和面源之外区域的续算，相比于使用sn方法直接计算整个区域，耦合方法通过定义曲面面源，利用mc建模更精确地描述反应堆局部的曲面区域，计算结果更为精准；相比于mc方法直接计算整个区域，耦合方法将规则几何区域划分给sn计算，通过减少mc计算区域降低了mc方法的计算时间，提高了计算效率。
4.为了实现以上目的，本发明采取如下的技术方案予以实施：步骤1：离散纵标方法sn输运计算得到反应堆几何结构包含曲面面源区域在内的
网格角注量率；步骤2：根据sn的z方向分层将曲面划分为环带，计算环带投影与sn网格面的交点，根据交点数量将环带等分为相同数量的环面网格，进而获得环面上网格与sn网格的映射关系；具体映射关系用环面网格编号与sn网格编号来对应：环面网格编号格式，sn网格编号格式，其中与一一对应，号通过公式映射得到：
ꢀꢀꢀꢀ
（1）式中：——某一环带对应的分层号；——在某一环带上的环面网格编号；——sn输运计算划分在方向上划分的网格编号；——sn输运计算在x方向划分的下边界、上边界网格号；——sn输运计算在y方向划分的下边界、上边界网格号；步骤3：获得网格映射关系后，将对应sn网格中心点的角注量率视为从环面网格处出射的角注量率；同时由于仅考虑沿面源正法线方向的出射流，需要计算曲面上不同位置的法向量：根据步骤2中计算的环带投影与sn网格面的交点以及曲面的中心坐标获得交点的法向量：球面某点法向量：，椭球面某点法向量：，其中：——曲面与sn网格的交点坐标；——球面和椭球面的球心坐标；——椭球面的三个方向的轴长；将法向量和sn计算的离散方向m进行向量点积运算，结果为正判定为正向出射；步骤4：根据角注量率依次求得方向m、网格q、分层k、能群g、面s、粒子类型p上的出射流；曲面计算出射流的流程与平面、柱面类似，有以下三点不同点：（1）通过对同一分层上网格出射流的求和，获得某一分层出射粒子流；（2）对曲面而言，法向量随网格位置变动，需根据粒子飞行方向与不同位置的
法向量的点积判断是否为正向出射；（3）计算从某一网格出射的粒子流时，网格面积应为曲面上的网格面积，推导得到的公式如下：球面网格：
ꢀꢀꢀꢀ
（2）椭球面网格：（3）
ꢀꢀ
（4）（5）
ꢀꢀꢀꢀꢀꢀꢀꢀ
（6）
ꢀꢀ
（7）式中：——曲面某一分层上某网格的面积；——椭球壳表面积；——某一分层的高度；——椭球面三个方向的轴长；——某一分层上的曲面网格数目；——横向截面的椭圆周率，；——与横向截面垂直的纵向椭圆的椭圆周率；——椭球面在不同高度下的缩放系数；——极坐标系下某点所处位置的极角；
步骤5：根据粒子出射流，计算得到方向m、网格q、分层k、能群g、面s、粒子类型p的累积分布函数cdf；步骤6：利用累积分布函数cdf进行源抽样，进行mc续算；为确保抽样获得的粒子几何位置能均布在球或椭球的表面，先根据分层k的累积分布函数随机抽样得到高度z，根据高度z确定截面方程，对于球面而言，其在截面圆网格边界的分布与圆心角成正比，对于椭球而言，为简化抽样过程，近似认为其在截面椭圆网格边界的分布与圆心角成正比，具体几何位置的抽样公式如下：球面：
ꢀꢀꢀꢀꢀ
（8） （9） （10）（11）
ꢀꢀꢀꢀꢀꢀꢀ
（12）椭球面：
ꢀꢀꢀꢀꢀ
（13）（14）（15）
ꢀꢀꢀꢀ
（16）
ꢀꢀ
（17） （18）式中：——几何位置抽样时的粒子位置坐标；
——椭球面三个方向的轴长；——球面或椭球面某一分层z处截面的高度；、 ——球面或椭球面某一高度分层的下边界，上边界；、——某一截面椭圆的长轴、短轴；、——抽样时产生的随机数，在之间；——圆心角，；——曲面网格极坐标系下对应的极角下边界、上边界；——球面的半径；——某一截面圆的半径。与现有技术相比，本发明有如下优点：在本发明方法中，通过曲面网格与sn计算网格的映射关系，近似求解得到曲面面源的累积分布函数，在对某些特殊几何结构的屏蔽计算中，利用曲面耦合计算可以兼顾sn计算的高效率和mc对于复杂几何计算的高精度。
附图说明
5.图1为本发明方法流程。
6.图2为曲面网格分解示意图。
7.图3为球面耦合计算与sn直接计算结果比较。
具体实施方式
8.下面结合附图和具体实施方式对本发明做进一步详细说明：一种应用于屏蔽计算曲面衔接中的s
n-mc耦合计算方法，在反应堆屏蔽计算中，某些复杂几何区域的形状为球面或椭球面，一般sn方法的计算网格无法准确描述这类耦合面，需要采用相交网格映射的方式近似求解。以球面为例，在sn完成比面源区域略大部分的输运计算后，根据sn计算划分的z方向网格将球面分为多层的圆形环带，计算环带投影与sn网格面的交点，利用交点坐标和环带所在高度求得曲面某位置的法向量，同时将交点所在sn网格中心点的角注量率视为环带映射区域的出射角注量率。利用角注量率可以求得粒子类型、面源、能群、曲面分层、网格、方向上的出射流，再将出射流转化为累积分布函数，作为mc程序的源项输入，进行曲面之外区域的续算。
9.如图1所示，本发明方法包括如下步骤：步骤1：离散纵标方法sn输运计算得到反应堆几何结构包含曲面面源区域在内的网格角注量率；步骤2：如图2所示，将图2中四张图按箭头排序为
①
，
②
，
③
，
④
。参照图例
①
，
②
根据sn的z方向分层将曲面划分为圆形环带；将环带投影到sn计算的xy平面上，如图例
③
所示，虚线区域为环带投影穿过的sn网格，同时通过计算环带投影与sn网格面的交点，根据交点数
量将环带等分为相同数量的环面网格，根据交点与sn网格中心点的对应关系建立起获得环面网格与sn网格的映射关系。具体映射关系可以用环带网格编号与sn网格编号来对应：曲面网格编号格式，sn网格编号格式，其中与一一对应，号可通过公式映射得到：
ꢀꢀꢀꢀ
（1）式中：——某一环带对应的分层号；——在该环带上的曲面网格编号；——sn输运计算划分在方向上划分的网格编号。
10.——sn输运计算在x方向划分的下边界、上边界网格号；——sn输运计算在y方向划分的下边界、上边界网格号；步骤3：获得网格映射关系后，将对应sn网格中心点的角注量率视为从环面网格处出射的角注量率，如图2中图例
④
所示；同时由于仅考虑沿面源正法线方向的出射流，需要计算曲面上不同位置的法向量：根据步骤2中计算的环带投影与sn网格面的交点以及定义曲面的中心坐标获得交点的法向量：球面某点法向量：，椭球面某点法向量：，其中：——曲面与sn网格的交点坐标；——球面和椭球面的球心坐标；——椭球面的三个方向的轴长；将法向量和sn计算的离散方向m进行向量点积运算，结果为正判定为正向出射；步骤4：根据角注量率依次求得方向m、网格q、分层k、能群g、面s、粒子类型p上的出射流；曲面计算出射流的流程与平面、柱面类似，有以下三点不同点：（1）通过对同一分层上网格出射流的求和，获得某一分层出射粒子流；（2）对曲面而言，法向量随网格位置变动，需根据粒子飞行方向与不同位置的法向量的点积判断是否为正向出射；（3）计算从某一网格出射的粒子流时，网格面积应为曲面上的网格面积，推导得到
的公式如下：球面网格：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（2）椭球面网格： （3）（4）
ꢀꢀꢀꢀꢀꢀꢀ
（5） （6）
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（7）式中：——曲面某一分层上某网格的面积；——椭球壳表面积；——某一分层的高度；——椭球面三个方向的轴长；——某一分层上的曲面网格数目；——横向截面的椭圆周率，；——与横向截面垂直的纵向椭圆的椭圆周率；
——椭球面在不同高度下的缩放系数；——极坐标系下某点所处位置的极角；步骤5：根据粒子出射流，计算得到方向m、网格q、分层k、能群g、面s、粒子类型p的累积分布函数cdf；步骤6：利用累积分布函数cdf进行源抽样，进行mc续算；为确保抽样获得的粒子几何位置能均布在球或椭球的表面，先根据分层k的累积分布函数随机抽样得到高度z，根据高度z确定截面方程，对于球面而言，其在截面圆网格边界的分布与圆心角成正比，对于椭球而言，为简化抽样过程，近似认为其在截面椭圆网格边界的分布与圆心角成正比，具体几何位置的抽样公式如下：球面：
ꢀꢀꢀ
（8）
ꢀꢀꢀꢀꢀꢀ
（9） （10）（11）（12）椭球面：
ꢀꢀꢀ
（13）
ꢀꢀ
（14）（15）（16）（17）
ꢀꢀꢀꢀꢀ
（18）式中：——几何位置抽样时的粒子位置坐标；——椭球面三个方向的轴长；——球面或椭球面某一分层z处截面的高度；、 ——球面或椭球面某一高度分层的下边界，上边界；、——某一截面椭圆的长轴、短轴；——抽样时产生的随机数，在之间；——圆心角，；——曲面网格极坐标系下对应的极角下边界、上边界；——球面的半径；——某一截面圆的半径。
11.为验证本发明的有效性，图3给出了对于一个立方体不锈钢的点源屏蔽问题采用sn程序直接计算和通过球面面源耦合计算的结果，结果表明在面源附近区域，两种计算方式的中子注量率误差在1%左右，证明了采用上述求解过程的可行性。
12.利用本发明求解某些存在上述特殊曲面结构的反应堆屏蔽问题时，使用sn程序计算大范围结构简单的屏蔽区，经由曲面面源使用mc程序计算特殊结构区，可以减少屏蔽计算的时间，同时保证对曲面区域的计算精度。