本发明针对核反应堆堆芯中子学计算领域，提出了一种针对反应堆中子扩散方程的非均匀几何变分节块方法。
背景技术：
核反应堆中子学计算研究以核反应堆堆芯为应用对象，其堆芯由许多不同种类的组件构成。根据堆型的不同，组件内部的几何结构和材料布置复杂多变。因此，实际的反应堆中子学问题是一个三维非均匀几何的中子学问题。对核反应堆进行快速、精确的中子学计算，是反应堆设计和校核的基本保障。目前在堆芯物理设计过程中，主要采用均匀几何的扩散方程求解方法。
变分节块法最早由美国西北大学的e.e.lewis教授提出，是中子学计算方法的杰出代表，具备扎实的工程应用背景。其主要应用作有美国阿贡国家实验室开发的variant程序和nodal程序，法国原子能委员会的eranos程序，以及美国爱达荷国家实验室的instant程序等等。
变分节块法以二阶偶宇称形式的中子扩散方程为出发点，方程呈现椭圆方程的形式，有利于garlerkin方法的应用，且更适合有限元方法的空间离散。变分节块法的计算思想是：首先通过变分方法在均匀几何求解区域建立包含二阶中子输运方程和自然边界条件的泛函；然后采用标准正交多项式进行ritz离散，同时利用球谐函数实现角度展开，并构造响应矩阵；最后分别在三维堆芯的各个节块内分别求解响应矩阵方程；节块之间以流和其高阶矩耦合，最终得到问题区域的中子通量密度分布。变分节块法消除了横向积分，离散对象直接针对三维中子通量密度分布。因此，变分节块法不需要精细功率重构技术，只需将最终求得的中子通量密度矩代入展开式就可以得到中子通量密度分布。然而，传统的变分节块法程序仅具备均匀节块的处理能力，不足以精细描述节块内部的非均匀几何，无法避免节块均匀化带来的误差。
随着科学技术的发展和计算机水平的提高，人们已经越来越开始重视减少近似和假设，追求更高精度的中子学计算方法，消除均匀化过程是中子学发展的必然趋势。因此，针对反应堆中子扩散方程的非均匀几何变分节块方法对于中子学计算具有十分重要的意义。
技术实现要素：
为了克服上述现有技术存在的问题，本发明的目的是提出一种针对反应堆中子扩散方程的非均匀几何变分节块方法，它能够精细描述核反应堆的非均匀栅元结构；该方法将基于变分节块法，采用等参有限元来处理节块内部的精细几何结构，实现非均匀几何的中子扩散方程求解。
为了实现上述目的，本发明采取了以下技术方案予以实施：
一种针对反应堆中子扩散方程的非均匀几何变分节块方法，步骤如下：
步骤1：首先根据公式(1)中二阶偶宇称扩散方程建立包含公式(3)中子通量密度φ和中子流密度j的泛函，泛函中包含节块内部的中子守恒关系以及节块表面的流连续性条件：
针对某一特定能群，在扩散近似下，二阶偶宇称扩散方程为：
式中：
φ—节块内部中子通量密度；
ω—方位角向量；
σt—中子宏观总截面；
σa—中子宏观吸收截面；
q—中子源项；
根据变分原理，在由若干节块组成的整个非均匀求解区域上，对应扩散方程的泛函写作各个节块内部及其表面上泛函的叠加贡献：
式中：
f[φ,j]—整个非均匀几何求解区域内的泛函；
fv[φ,j]—单个节块内部的泛函；
v—节块的编号；
而扩散近似下的各节块泛函
其中γ是外部边界；
步骤2：有限元形状函数g(x,y)能够用来描述曲边几何结构，因此利用x‐y方向的有限元形状函数g(x,y)，分片常量多项式h(x,y)，z方向的正交多项式fz(z)和f′z(z)对节块内部中子通量密度φ、节块表面中子流j分别展开，实现非均匀节块的几何和材料描述功能：
式中：
t—转置符号；
fz(z)—节块内部轴向标准正交多项式向量；
g(x,y)—x‐y方向有限元形状函数向量；
φ—节块内部中子通量密度的展开矩向量；
f′z(z)—节块x‐y表面的轴向标准正交多项式向量；
fγ(γ′)—节块x‐y表面的径向标准正交多项式向量；
j±γ(γ′,z)—节块x‐y表面中子流密度展开矩向量，其中展开矩代表了展开系数的值；j±γ(γ′,z)是关于径向方向上的自变量γ′＝x,y和轴向方向的自变量z的函数：当γ＝x时γ′＝y，j±x(y,z)代表节块左侧和右侧的表面中子流密度展开矩向量；当γ＝y时γ′＝x，j±y(x,z)代表节块下侧和上侧的表面中子流密度展开矩向量；
j±z(x,y)—节块z表面中子流密度展开矩向量，它是关于径向方向上的自变量x,y的函数；
δz—节块z方向上的高度；
h(x,y)—节块z表面分片常量，对应于节块内部第e个面积为ae的有限元网格，它满足
步骤3：将公式(4)至公式(6)中的离散表达式代入各个节块内部的泛函公式(3)，得到表征中子通量密度展开系数φ和中子流密度展开系数j之间关系的响应矩阵方程公式(8)、公式(10)：
将公式(4)至公式(6)代入公式(3)，得节块内部泛函的离散形式：
根据变分原理，对公式(7)取φ的一阶变分为0，得扩散近似形式的矩阵方程：
其中：
其中：
-1—矩阵的求逆；
φ—节块内部中子通量密度展开矩向量；
j—节块表面净中子流密度展开矩向量；
q—中子源项展开矩向量；
对公式(7)取j的一阶变分为0，得节块表面偶宇称中子通量密度展开矩
联立公式(8)和公式(9)：
式中：
j—净中子流密度展开矩向量；
u—中子流源项展开矩向量；
其中：
为了将响应矩阵表达成通用形式，利用变量替换关系式
将公式(10)写为响应矩阵形式：
式中：
j+—出射中子流展开矩向量；
j-—入射中子流展开矩向量；
式中：
步骤4：对公式(14)、公式(8)所代表的响应矩阵方程，利用红-黑迭代的方法进行迭代求解，最终得到整个非均匀几何求解区域的中子通量密度分布φ(r)和中子流密度分布j±γ(γ′,z)、j±z(x,y)，从而完成针对反应堆中子扩散方程的非均匀几何变分节块方法。
与现有技术相比，本发明具有如下突出优点：
1.传统的节块方法绝大多数仅具备均匀节块的描述能力，本发明通过采用等参有限元对各节块内部的中子通量密度分布进行空间离散，可以有效地描述各个节块的内部的非均匀曲边几何结构，实现非均匀扩散计算的能力。
2.本发明通过在节块的轴向表面采用分片常量进行空间离散，能够实际描述非均匀节块表面的轴向泄漏分布，消除传统节块方法中轴向表面的均匀化过程，更加切合物理实际，同时提高计算精度。
附图说明
图1为等参有限元描述下的压水堆栅元非均匀几何。
具体实施方式
下面结合具体实施方式对本发明作进一步详细说明。该方法采用标准的源迭代的方法进行外迭代。针对群内迭代，具体计算流程包括以下几步：
步骤1：首先根据公式(1)中二阶偶宇称扩散方程建立包含公式(3)中子通量密度φ和中子流密度j的泛函，泛函中包含节块内部的中子守恒关系以及节块表面的流连续性条件：
针对某一特定能群，在扩散近似下，二阶偶宇称扩散方程为：
式中：
φ—节块内部中子通量密度；
ω—方位角向量；
σt—中子宏观总截面；
σa—中子宏观吸收截面；
q—中子源项；
根据变分原理，在由若干节块组成的整个非均匀求解区域上，对应扩散方程的泛函可以写作各个节块内部及其表面上泛函的叠加贡献：
式中：
f[φ,j]—整个非均匀几何求解区域内的泛函，其中φ代表非均匀几何求解区域的中子角通量密度，j代表非均匀几何求解区域的表面中子流密度；
fv[φ,j]—单个节块内部的泛函，其中φ代表节块内部的中子角通量密度，j代表节块表面的中子流密度；
v—节块的编号；
而扩散近似下的各节块泛函
其中γ是外部边界；
步骤2：有限元形状函数g(x,y)能够用来描述如图1所示的曲边几何结构，因此利用x‐y方向的有限元形状函数g(x,y)，分片常量多项式h(x,y)，z方向的正交多项式fz(z)和f′z(z)对节块内部中子通量密度φ、节块表面中子流j分别展开，实现非均匀节块的几何和材料描述功能：
式中：
t—转置符号；
fz(z)—节块内部轴向标准正交多项式向量；
g(x,y)—x‐y方向有限元形状函数向量；
φ—节块内部中子通量密度的展开矩向量；
f′z(z)—节块x‐y表面的轴向标准正交多项式向量；
fγ(γ′)—节块x‐y表面的径向标准正交多项式向量；
j±γ(γ′,z)—节块x‐y表面的中子流密度展开矩向量，其中展开矩代表了展开系数的值；它是关于径向方向上的自变量γ′＝x,y和轴向方向的自变量z的函数：当γ＝x时γ′＝y，j±x(y,z)代表节块左侧和右侧的表面中子流密度展开矩向量；当γ＝y时γ′＝x，j±y(x,z)代表节块下侧和上侧的表面中子流密度展开矩向量；
j±z(x,y)—节块z表面的中子流密度展开矩向量，它是关于径向方向上的自变量x,y的函数，脚标±z分别对应了节块的顶侧和低侧表面；
δz—节块z方向上的高度；
h(x,y)—节块z表面分片常量，对应于节块内部第e个面积为ae的有限元网格，它满足
步骤3：将公式(4)至公式(6)中的离散表达式代入各个节块内部的泛函公式(3)，得到表征中子通量密度展开系数φ和中子流密度展开系数j之间关系的响应矩阵方程公式(8)、公式(10)：
将公式(4)至公式(6)代入公式(3)，得节块内部泛函的离散形式：
根据变分原理，对公式(7)取φ的一阶变分为0，得扩散近似形式的矩阵方程：
其中：
其中：
-1—矩阵的求逆；
φ—节块内部中子通量密度展开矩向量；
j—节块表面净中子流密度展开矩向量；
q—中子源项展开矩向量；
对公式(7)取j的一阶变分为0，得节块表面偶宇称中子通量密度展开矩
联立公式(8)和公式(9)：
式中：
j—净中子流密度展开矩向量；
u—中子流源项展开矩向量；
其中：
为了将响应矩阵表达成通用形式，利用变量替换关系式
将公式(10)写为响应矩阵形式：
式中：
j+—出射中子流展开矩向量；
j-—入射中子流展开矩向量；
式中：
为使表达简便，省略展开基函数中的x,y,z自变量，公式(8)至公式(16)中的各系数矩阵可写作：
其中：
g|±x—g(x,y)函数在节块+x,-x表面上的取值；
g|±y—g(x,y)函数在节块+y,-y表面上的取值；
fz|±z—fz(z)函数在节块+z,-z表面上的取值；
其他符号的代表意义与前文相同。
步骤4：对于求解区域内各不同种类的节块，根据公式(11)，公式(15)至公式(28)分别计算各节块的响应矩阵；
步骤5：针对特定能群，利用红‐黑扫描的方式对响应矩阵方程公式(14)进行迭代求解，得到出、入射偏中子流密度矩j+、j-；在求得出、入射偏中子流密度矩j+、j-后，代入公式(8)，即可解得节块内部中子通量密度矩φ，从而可由公式(4)得知节块内部的中子通量密度分布φ(r)；
步骤7：进行下一能群的计算，能够求解整个非均匀求解区域内各群的中子通量密度分布和有效增值因子，完成针对反应堆中子扩散方程的非均匀几何变分节块方法。