随着堆芯计算的高精度需求以及堆芯燃料的高燃耗发展，严格追踪燃耗过程中的核素成份变化十分重要。有关燃耗问题的一个重要部分便是求解燃耗方程：
| $\begin{align} & \frac{\operatorname{d}{{N}_{i}}}{\operatorname{d}t}=\sum\limits_{i\ne j}{{{\lambda }_{ij}}{{N}_{j}}}-{{\lambda }_{i}}^{\rm eff}{{N}_{i}} \\ & {{N}_{i}}(0)={{N}_{0}}\, \ i=1, \cdot \cdot \cdot, n \\ \end{align}$ | (1) |
式中：Ni为第i种核素的核密度；N0为初始核素密度；n为核素的总数；λij为中子诱发反应率和自发衰变率之和；λijeff为第i种核素的消失率。
目前燃耗方程的求解主要有两类算法：一是基于单条燃耗链求解的线性子链法(Transmutation Trajectory Analysis, TTA)；二是基于燃耗矩阵求解的矩阵指数法[1]。前者是将燃耗过程中复杂的非线性链简化为一系列的线性燃耗链进行解析求解，其优势是避免了不同核素之间由于半衰期、反应截面等参数差距过大而造成的刚性问题，但燃耗过程中核素反应关系复杂，拆解建立燃耗链的工作相当复杂繁重，计算效率相对较低；后者则针对燃耗过程中每一种核素，考虑该核素的每种反应率，建立燃耗矩阵，基于燃耗矩阵求解，常用的求解方法有Chebyshev有理近似方法(Chebyshev Rational Approximation Method, CRAM)、Padé近似方法、Krylov Subspace方法[2-4]等，这类求解方法充分考虑到多种核素之间的相互作用，没有物理上的近似，但缺点是由于燃耗矩阵的刚性不可避免地会对短寿命核素的计算存在误差，同时时间步长的选取也可能存在困难。
本文基于这两类算法，分别采用TTA、CRAM、Padé近似以及Krylov Subspace方法，对三种规模的燃耗系统进行燃耗计算，并通过对比分析，选择高精度下精细步长的16阶系数CRAM的结果作为参考，比较各计算方法的性能优劣。
1 燃耗计算方法本文研究过程中，不考虑燃耗过程中中子通量及反应截面的更新，这对于本文的研究目的并没有影响，此时式(1)是一阶常微分方程组。
1.1 TTATTA的核心是将复杂的核素反应关系网分解为一系列的线性链，基于深度优先查找，从一种初始核密度非零的核素开始，考虑所有可能的反应方式，建立线性链，从而计算链中每种核素的核密度[5]。
假设一条线性链其有效衰变常数为λijeff, i=1, …, n，只有第一种核素有非零的初始核密度x1(0)，则经Δt时间后，第n种核素的核密度[3]为：
| $ {x_n}(\Delta t) = {x_1}(0){B_n}\sum\limits_{i = 1}^n {\alpha _i^n{{\rm{e}}^{-\lambda _i^{{\rm{eff}}}\Delta t}}} $ | (2) |
其中：
| ${B_n} = \prod\limits_{j = 2}^n {b_{j - 1, j}^{{\rm{eff}}}} \;, \;b_{i, j}^{{\rm{eff}}} = \frac{{{b_{i, j}}{\lambda _i} + {\sigma _{i, j}}\phi }}{{\lambda _i^{{\rm{eff}}}}}$ | (3) |
| $\alpha _{i}^{n}={\prod\limits_{j=1}^{n-1}{\lambda _{j}^{\text{eff}}}}/{\prod\limits_{\begin{smallmatrix} j=1 \\ j\ne i \end{smallmatrix}}^{n}{(\lambda _{j}^{\text{eff}}-\lambda _{i}^{\text{eff}})}}\;$ | (4) |
由式(4)可观察到，当反应中存在环链时，该方法无法应用，同时如果两核素的有效衰变常数非常接近时也容易引起数值不稳定，因此本文解决方法是在拆解建立线性链的过程中，遇到环链则断开该环链，当两核素的有效衰变常数相对比值小于1.0×10-6时，将其一核素的有效衰变常数变为1.00001倍。截断误差设定为：
| $ b_{n, j}^{{\rm{eff}}}{P_n} < \frac{{{x_1}(0)}}{{{x_{tot}}(0)}}\rm cutoff $ | (5) |
| $ {P_n} = \frac{{{x_{n + 1}}(t)}}{{{x_1}(0)}}\left| {_{\lambda _{n + 1}^{{\rm{eff}}}{\rm{ = }}0}} \right. $ | (6) |
式中：cutoff为设定的截断份额，本文设定cutoff=10-15。
当bn, jeffPn满足式(5)时，该种核素被忽略，意味着该条燃耗链至此也将断开。
1.2 矩阵指数法式(1)可以写成矩阵形式：
| $\frac{{{\mathop{\rm d}\nolimits} {\pmb{N}}}}{{{\mathop{\rm d}\nolimits} t}} = {\pmb{AN}}\;, \;{\pmb{N}}(0) = {{\pmb{N}}_0}$ | (7) |
式中：N为核素密度矢量；A为燃耗矩阵。
式(7)的通解为：
| $ {\pmb{N}}(t) = {{\mathop{\rm e}\nolimits} ^{{\pmb{A}}t}}{{\pmb{N}}_0}\;\;, \;\;{\pmb{N}}(0) = {{\pmb{N}}_0} $ | (8) |
显然燃耗方程中燃耗矩阵指数eAt的求解便是燃耗计算的关键，其求解方法决定了整个燃耗计算的精度和效率。
1.2.1 CRAMCRAM源于1969年，Cody、Meinardus和Varga[6]将e-x在[0, ∞)进行近似处理。近似函数rk, l的表达式为：
| ${{r}_{k,l}}={{{p}_{k}}(x)}/{{{p}_{l}}(x)}\;$ | (9) |
通过Remes algorithm或Carathéodory-Fejér method[7]可以确定k和l[8]。根据留数定理，将式(9)写成如下表达式：
| ${r_{k, k}}(x) = {\alpha _0} + \sum\limits_{i = 1}^k {\frac{{{\alpha _i}}}{{x - {\theta _i}}}} $ | (10) |
式中：α0为x取无穷时取值；αi为留数；θi为极点。由于αi、θi共轭出现，核密度的取值为实数，故式(10)的计算可以缩减一半：
| ${r_{k, k}}(x) = {\alpha _0} + 2Re(\sum\limits_{i = 1}^{k/2} {\frac{{{\alpha _i}}}{{x - {\theta _i}}})} $ | (11) |
该方法的计算精度和效率与k的取值紧密相关。本文直接引用14阶和16阶系数[9-10]，并通过Carathéodory-Fejér方法获得12阶系数。图 1给出高精度下在负实轴上16、14和12阶系数的CRAM结果r16, 16(x)、r14, 14(x)和r12, 12(x)与ex的误差图。从图 1可以观察到，系数阶数越高，精确程度也越高。
|
图 1 负实轴上不同阶数CRAM的计算结果rk, k(x)与ex的误差图 (a) (ex-r16, 16(x))，(b) (ex-r14, 14(x))，(c) (ex-r12, 12(x)) Figure 1 Plots of errors of ex and rk, k(x) based on CRAM with different orders on the negative real axis (a) (ex-r16, 16(x)), (b) (ex-r14, 14(x)), (c) (ex-r12, 12(x)) |
因此式(8)可写为：
| $\mathit{\boldsymbol{N}}(t) = {\alpha _0}{\mathit{\boldsymbol{N}}_0} + 2Re(\sum\limits_{i = 1}^{k/2} {\frac{{{\alpha _i}}}{{\mathit{\boldsymbol{A}}t - {\theta _i}}})} {\mathit{\boldsymbol{N}}_0}$ | (12) |
对eAt的(p, q)阶Padé近似定义式[11]为：
| $ {{\mathop{\rm e}\nolimits} ^{\mathit{\boldsymbol{A}}t}} \approx \frac{{{N_{pq}}(\mathit{\boldsymbol{A}}t)}}{{{D_{pq}}(\mathit{\boldsymbol{A}}t)}} $ | (13) |
其中：
| ${N_{pq}}(\mathit{\boldsymbol{A}}t) = \sum\limits_{j = 0}^p {\frac{{(p + q - j)!p!}}{{(p + q)!j!(p - j)!}}{{(\mathit{\boldsymbol{A}}t)}^j}} $ |
| ${D_{pq}}(\mathit{\boldsymbol{A}}t) = \sum\limits_{j = 0}^q {\frac{{(p + q - j)!q!}}{{(p + q)!j!(q - j)!}}{{( - \mathit{\boldsymbol{A}}t)}^j}} $ |
由于Padé近似方法中存在一系列的展开式，当||At||较大时，可能会出现收敛速度慢和圆整误差等问题，因此在Padé近似方法中需应用缩放乘方技术[12]提高其准确率，即选定参数m，使得式(14)中||A/m||充分小，且为便于计算设定m满足m=2k。
| ${{\rm{e}}^{\mathit{\boldsymbol{A}}t}} = {({{\rm{e}}^{\mathit{\boldsymbol{A}}t/m}})^m}$ | (14) |
本文中采用变阶数的对角Padé近似方法（即p=q），阶数可选值包括[3, 5, 7, 9, 13]，根据一定计算原则[12]确定其阶数。
1.2.3 Krylov Subspace方法Krylov Subspace方法的基本思想是将大型稀疏的燃耗矩阵投影到低阶的Krylov子空间，然后再通过Padé近似方法或CRAM等经典矩阵指数法来求解降阶子空间矩阵。
该方法的基本原理是首先通过Arnoldi过程获得Krylov Subspace方法的标准正交基[
| $ \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_m} = {\mathit{\boldsymbol{V}}_m}{\mathit{\boldsymbol{H}}_m} + {h_{m + 1}}_{, m}{\vec v_{m + 1}}{\vec e_m}^{\rm{T}} $ | (15) |
式中：Vm=[
Hm∈Rm×m为Arnoldi过程产生的Hessenberg上三角矩阵。
通过式(15) Hessenberg矩阵可以改写为：
| $ {\mathit{\boldsymbol{H}}_m} = {\mathit{\boldsymbol{V}}_m}^{\rm{T}}\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_m} $ | (16) |
设式(17)来近似N(t+Δt)：
| $ \mathit{\boldsymbol{N}}(t + \Delta t) \approx {\mathit{\boldsymbol{N}}_q}(t + \Delta t) = {\mathit{\boldsymbol{V}}_m}y $ | (17) |
式中：y是m维系数组合。
考虑到
| $ \begin{array}{l} y = {({\mathit{\boldsymbol{V}}_m}^T{\mathit{\boldsymbol{V}}_m})^{-1}}{\mathit{\boldsymbol{V}}_m}^{\rm{T}}{\mathit{\boldsymbol{N}}_q}(t + \Delta t)\\ \;\; = {({\mathit{\boldsymbol{V}}_m}^T{\mathit{\boldsymbol{V}}_m})^{-1}}{\mathit{\boldsymbol{V}}_m}^{\rm{T}}{{\rm{e}}^{\mathit{\boldsymbol{A}}\Delta }}^t\mathit{\boldsymbol{N}}(t)\\ \;\; = {\mathit{\boldsymbol{V}}_m}^{\rm{T}}{{\rm{e}}^{\mathit{\boldsymbol{A}}\Delta }}^t\mathit{\boldsymbol{N}}(t) \end{array} $ | (18) |
又
| $ \mathit{\boldsymbol{N}}(t) = \beta {\vec v_1} = \beta {{\rm{V}}_m}{\vec e_1} $ | (19) |
其中：
| $ \beta {\rm{ = }}{\left\| {\mathit{\boldsymbol{N}}\left( t \right)} \right\|_2};\;{\overrightarrow v _1} = \mathit{\boldsymbol{N}}\left( t \right)/\beta ;\;\;{\overrightarrow e _1} = {\left[{1, 0, \cdots, 0} \right]_m} $ |
因此将式(16)、(18)、(19)带入式(17)中可得：
| $ {\mathit{\boldsymbol{N}}_q}(t + \Delta t) = \beta {\mathit{\boldsymbol{V}}_m}{{\rm{e}}^{{\mathit{\boldsymbol{H}}_m}\Delta t}}{\vec e_1} $ | (20) |
至此，便可将燃耗矩阵降至m维，同时文献[14]中指出，式(21)的修正式能够进一步提高精度：
| $ {\mathit{\boldsymbol{N}}_q}(t + \Delta t) = \beta {V_{m + 1}}{{\rm{e}}^{{{\mathit{\boldsymbol{\bar H}}}_{m + 1}}\Delta t}}{\vec e_1} $ | (21) |
其中：
|
(22) |
式中：hm+1, m为Arnoldi过程获得的Hessenberg矩阵元素。
2 数值结果及分析将前文所述的4种求解方法应用到3种不同规模的燃耗算例中。3种燃耗算例均在通量为3.0×1014cm-2·s-1下燃耗，总燃耗时长为300 d（共19个燃耗步，燃耗步长分别设定[1 d×2, 2 d×4, 10d×3, 20 d×8, 50 d×2]）。算例1只有34种核素，包括22种锕系核素和12种裂变产物；算例2和算例3分别有385种和1306种核素。算例1的初始核素组成为234U:235U:238U:239Pu=0.03:3.10:96.49: 0.38；算例2包括5种初始核素，其核素核密度为16O:234U:235U:238U:239Pu=66.71:0.01:1.00:0.01:32.27；算例3中初始核素采用燃耗过的燃料，包含153种核素。在对比分析中发现，采用精细步长，高精度下CRAM和Padé近似计算结果非常一致，CRAM计算又十分迅速，因而选取高精度下16阶系数的CRAM (CRAM16)精细步长(0.25 d)的计算结果作为参考，同时核密度小于1.0×10-30 b-1·cm-2的核素将被忽略。
图 2为3种算例在不同求解方法下的核密度相对误差，其中Krylov Subspace方法在算例1和2中采用的子空间维数分别为5和50，同时由于该方法中内部燃耗步长是基于局部截断误差估计[13]自动计算的，而在大规模燃耗系统中步长会出现趋于无穷的现象，因此该方法不适用于算例3。由图 2(a)和(b)可观察到，在中小规模燃耗系统下，各个方法的结果都很好，CRAM和Padé近似效果最佳。而对于图 2(c)，需要注意的是，如同Pusa在文献[2]中指出的，Padé近似方法在大规模燃耗系统中会出现不合理的核素密度结果，在计算中发现如果采用双精度确实会产生该现象，但在高精度下Padé近似的结果如图 2(c)所示，与CRAM的结果相当一致。3种算例下，TTA和高精度精细步长的CRAM16误差较大，事实上由于本文中TTA采取遇环链便断开此链的处理方式，所得到的核素密度相应偏低。
|
图 2 核素密度相对误差 (a)算例1，(b)算例2，(c)算例3 Figure 2 Relative errors of nuclide density (a) Case 1, (b) Case 2, (c) Case 3 |
由于CRAM和TTA适用于各燃耗规模，进一步考察其重要核素相对误差的均根值(Root Mean Square, RMS)随燃耗深度增加的变化情况，如图 3所示。其中算例1直接考察所有核素的相对误差的RMS，算例2和算例3分别考虑50和119种重要核素。由图 3(a)、(b)可观察到，CRAM和TTA方法下中小燃耗系统中重要核素整体相对误差的RMS随燃耗加深整体呈减小的趋势，受燃耗深度的影响较小。TTA方法后期有所上升，猜测是由于燃耗步长的增加引入的误差增大。图 3(c)表明大规模燃耗系统下，重要核素的相对误差变化受燃耗深度的影响较小。同时，3种算例下重要核素相对误差的RMS结果也表明，CRAM整体误差稳定性高于TTA，而且CRAM的系数越高，其误差稳定性越好，也越适用于大规模燃耗系统。
|
图 3 CRAM与TTA重要核素相对误差的RMS (a)算例1，(b)算例2，(c)算例3 Figure 3 RMS of the relative errors of important nuclides with CRAM and TTA (a) Case 1, (b) Case 2, (c) Case 3 |
注意到高精度下Padé近似的准确度较高，双精度下CRAM和TTA的应用性较强，进一步考察其步长稳定性。通常情况下，在燃耗方程求解中，燃耗步长取值越小，计算结果就越精细，但过于精细的燃耗步长对于计算效率是不利的。因此为了保证计算结果的精确度和计算效率，需要兼顾燃耗步长的选取以及其对计算结果的影响。
仍然选用前文3种算例模型，设定在通量3.0×1014cm-2·s-1下等步长燃耗至100 d，各等步长分别设为[0.5 d, 1 d, 2 d, 5 d, 10 d, 20 d, 25 d, 50 d, 100d]，参考值为各方法下更为精细的步长(0.25 d)的计算结果，分别考察第100 d得到的各方法下核素密度受步长变化的影响情况。如图 4和5给出算例2模型分别在高精度和双精度下Padé近似方法获得的4种重要锕系核素和3种重要裂变产物的步长稳定性变化情况。从图 4中可以观察到，Padé近似方法在高精度下核密度结果几乎不受步长变化的影响，步长稳定性非常高，而在双精度下受步长变化的影响也在10-9以内，因此通常对于中小规模的燃耗计算，采用Padé近似方法时可采取较大步长以提高效率。图 6为双精度下算例3模型下3种不同阶数的CRAM和TTA的步长稳定性情况，锕系核素以238U为例，裂变产物以135Xe为例。由此可知，CRAM的步长稳定性远优于TTA的步长稳定性，CRAM不论锕系核素还是裂变产物其核密度受燃耗步长变化的影响较小，误差维持在10-9量级以下，而且系数阶数越高，其稳定性越好。其他算例模型的结果与之类似，不作赘述。
|
图 4 算例2模型中Padé近似在高精度下重要核素的步长稳定性 (a)重要锕系核素，(b)重要裂变产物 Figure 4 Step stability of important nuclides with Padé approximation with high precision in case 2 (a) Important actinides, (b) Important fission products |
|
图 5 算例2模型中Padé近似在双精度下重要核素的步长稳定性 (a)重要锕系核素，(b)重要裂变产物 Figure 5 Step stability of important nuclides with Padé approximation with double precision in case 2 (a) Important actinides, (b) Important fission products |
|
图 6 算例3模型中CRAM和TTA下重要核素的步长稳定性 (a) 238U，(b) 135Xe Figure 6 Step stability of important nuclides with CRAM and TTA in case 3 (a) 238U, (b) 135Xe |
最后考察各方法在效率方面的表现，燃耗模块基于FORTRAN平台进行编程计算，所有算例均在2.10 GHz的计算机下完成，分别统计前文3种算例模型在双精度下各求解方法的单步燃耗(10 d)耗时情况，如表 1所示。其中CRAM对于矩阵求逆部分，由于采取了矩阵稀疏化[15-16]处理策略，并通过高斯迭代求解，即使在大规模燃耗系统中计算也非常迅速，系数阶数的增长几乎没有增加时间成本。TTA由于需要不断的拆解建立燃耗链耗时稍长，Padé近似中时间主要消耗在n阶矩阵多项式相除的高斯求解中，但如果采用等燃耗步长，则可省去后续的高斯求逆成本，整体计算效率可提升。Krylov Subspace方法内部子步长是基于局部误差估计自动划分，同时降阶矩阵指数用Padé近似计算，耗时相对较长。
|
表 1 各方法下单步燃耗耗时 Table 1 Computing time of single burnup for all methods |
分别采用CRAM、Padé近似方法、Krylov Subspace方法以及TTA对三种不同规模的燃耗系统进行燃耗分析，并考察核素密度相对误差、重要核素相对误差的RMS随燃耗的变化、步长稳定性和效率等相关参数，结果表明：
1) CRAM适用各种规模下的燃耗系统，计算精确度高，而且系数阶数越高，精确度越高，步长稳定性越好，同时计算效率最优。
2) Padé近似方法在中小规模的燃耗系统中表现优良，大规模的燃耗系统中需要采用高精度计算，而且在高精度下，其精确度非常高，步长稳定性十分优良。
3) Krylov Subspace方法在双精度下比较适合于中小规模的燃耗系统，同时求解中还需借用CRAM或Padé近似来进行与Hessenberg矩阵相关的矩阵计算，计算精确度较高，但不适用于大规模燃耗矩阵。
4) TTA方法同样适用于各规模的燃耗矩阵，但目前由于线性链的简化限制，相对参考值精确度较低，各规模下的精度程度相当，步长稳定性也不如Padé和CRAM。
本研究对于不同燃耗计算系统的开发及应用具有较好的借鉴意义。所研究的矩阵指数法将与蒙特卡罗输运相结合，采用预估修正等燃耗策略，实现蒙特卡罗燃耗耦合计算功能。
| [1] |
Cetenar J. General solution of bateman equations for nuclear transmutations[J]. Annuals of Nuclear Energy, 2006, 33: 640-645. DOI:10.1016/j.anucene.2006.02.004 |
| [2] |
Pusa M, Leppänen J. Computing the matrix exponential in burnup calculations[J]. Nuclear Science and Engineering, 2010, 164(Pt23): 140-150. DOI:10.13182/NSE09-14 |
| [3] |
Isotalo A E, Aarnio P A. Comparison of depletion algorithms for large systems of nuclides[J]. Annals of Nuclear Energy, 2011, 38(2): 261-268. DOI:10.1016/j.anucene.2010.10.019 |
| [4] |
范文玎, 孙光耀, 张彬航, 等. 基于切比雪夫有理逼近方法的蒙特卡罗燃耗计算研究与验证[J]. 核技术, 2016, 39(4): 040501. FAN Wending, SUN Guangyao, ZHANG Binhang, et al. Research and verification of Monte Carlo burnup calculations based on Chebyshev rational approximation method[J]. Nuclear Techniques, 2016, 39(4): 040501. DOI:10.11889/j.0253-3219.2016.hjs.39.040501 |
| [5] |
吴明宇, 王事喜, 杨勇, 等. 回溯算法在燃耗计算中的应用[J]. 原子能科学技术, 2013, 47(7): 1127-1132. WU Mingyu, WANG Shixi, YANG Yong, et al. Application of backtracking algorithm to depletion calculations[J]. Atomic Energy Science and Technology, 2013, 47(7): 1127-1132. DOI:10.7538/yzk.2013.47.07.1127 |
| [6] |
Cody W J, Meinardus G, Varga R S. Chebyshev rational approximations to exp (-x) in[0, +∞) and applications to heat-conduction problems[J]. Journal of Approximation Theory, 1969, 2(1): 50-65. DOI:10.1016/0021-9045(69)90030-6 |
| [7] |
Schmelzer T. Carathéodory-Fejér approximation, MATLAB central[CP/OL]. 2008. http://www.mathworks.com/matlabcentral.
|
| [8] |
Carpenter A J, Ruttan A, Varga R S. Extended numerical computations on the"1/9"conjecture in rational approximation theory[J]. Lecture Notes in Mathematics, 1984, 1105(4): 383-411. DOI:10.1007/BFb0072427 |
| [9] |
Pusa M. Correction to partial fraction decomposition coefficients for Chebyshev rational approximation on the negative real axis[J/OL]. arXiv: 1206. 2880[math. NA], 2012. http://www.oalib.com/paper/4122717
|
| [10] |
张竞宇, 马亚栋, 陈义学, 等. CRAM在放射性核素存量计算中的应用[J]. 核技术, 2017, 40(8): 080502. ZHANG Yadong, MA Yadong, CHEN Yixue, et al. Research and verification of Monte Carlo burnup calculations based on Chebyshev rational approximation method[J]. Nuclear Techniques, 2017, 40(8): 080502. DOI:10.11889/j.0253-3219.2017.hjs.40.080502 |
| [11] |
Moler C, Loan C F. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[J]. SIAM Review, 2003, 45(1): 3-49. DOI:10.1137/S00361445024180 |
| [12] |
Higham N J. The scaling and squaring method for the matrix exponential revisited[J]. SIAM Review, 2009, 51(4): 747-764. DOI:10.1137/090768539 |
| [13] |
Sidje R B. Expokit:a software package for computing matrix exponentials[J]. ACM Transactions on Mathematical Software, 1998, 24(1): 130-156. DOI:10.1145/285861.285868 |
| [14] |
Yamamoto A, Tatsumi M, Sugimura N. Numerical solution of stiff burnup equation with short half lived nuclides by the Krylov Subspace method[J]. Journal of Nuclear Science and Technology, 2007, 44(2): 147-154. DOI:10.1080/18811248.2007.9711268 |
| [15] |
James R B, Donald J R. Sparse matrix computations[M]. New York: Academic Press, 1976, 3-22. DOI:10.1016/B978-0-12-141050-6.50006-4
|
| [16] |
Pusa M, Leppänen J. Solving linear systems with sparse gaussian elimination in the Chebyshev rational approximation method[J]. Nuclear Science and Engineering, 2013, 175(3): 250-258. DOI:10.13182/nse12-52 |