1、第58 卷第3期2024年3月原子能科学技术Atomic Energy Science and TechnologyVol.58,No.3Mar.2024三维球床几何稀疏条数长特征线加速方法郭建，郭炯，李富2*，严睿，邹杨！（1.中国科学院上海应用物理研究所，上海2 0 18 0 0；2.清华大学核能与新能源技术研究院先进核能技术协同创新中心先进反应堆工程与安全教育部重点实验室，北京10 0 0 8 4)摘要：特征线法具有非常强的几何适应性，可用于三维球床高温气冷堆全堆芯求解，但存在选代次数多、计算速度慢的缺点。本文将长特征线加速方法应用于三维球床高温气冷堆以解决非规则几何数值加速的问题，基于
2、三维模型特征线布置较二维模型稠密的分析，提出了稀疏条数长特征线加速方法，极大地减少了加速方程的计算量，在不降低角度离散精度的前提下，获得了非常好的加速效果。通过基准题对加速参数的选取方式进行了研究，条数稀疏度取35、长特征线长度取2.0 cm左右、加速选代步取206 0 步可获得良好的加速效果。小型轻水堆三维基准题和球床堆芯简化模型的计算结果表明，采用稀疏条数长特征线加速可获得7 倍左右的时间加速比，此时对应的迭代步加速比为2 0 倍左右。关键词：球床高温气冷堆；特征线法；长特征线加速方法；中子输运中图分类号：TL325doi:10.7538/yzk.2023.youxian.0358Spar
3、se Macro-track Transport Acceleration MethodGUO Jian,GUO Jiong”,LI Fu.*,YAN Rui,ZOU Yang(1.Shanghai Institute of Applied Physics,Chinese Academy of Sciences,Shanghai 201800,China;2.Institute of Nuclear and New Energy Technology,Collaborative Innovation Centerof Advanced Nuclear Energy Technology,Key
4、 Laboratory of Advanced Reactor Engineering andSafety of Ministry of Education,Tsinghua University,Beijing 100084,China)Abstract:The method of characteristics has very strong geometric adaptability and canbe used to solve the 3-D whole core problem of high temperature gas-cooled pebble-bedreactor.Ho
5、wever,it suffers from the disadvantages of thousands of iterations and slowcalculation speed.At present,there have been some attempts to reduce calculation timefor method of characteristics,including coarse mesh finite difference(CMFD),generalcoarse mesh rebalance(GCMR)and general coarse mesh finite
6、 difference(GCMFD).The CMFD method is applied to realize the acceleration of regular region like squares.Furthermore,to break through the limitation of macroscopic regular geometry for thenumerical acceleration method,researchers successively propose the GCMR method and文献标志码：Afor 3-D Pebble-bed Geom
7、etry文章编号：10 0 0-6 9 31(2 0 2 4)0 3-0 6 54-0 8收稿日期：2 0 2 3-0 5-17；修回日期：2 0 2 3-0 8-16基金项目：中国科学院引才计划青年项目（SINAP-YCJH-202303）；中国科学院战略性先导科技专项（XDA02010000）；中国科学院前沿科学重点研究（QYZDY-SSW-JSC016）；国家自然科学基金（12 0 0 52 90）*通信作者：李富第3期郭建等：三维球床几何稀疏条数长特征线加速方法the GCMFD method,but these methods are very complex to implemen
8、t in irregulargeometries and some approximation factors must be introduced.Based on the analysis ofmore dense tracks are arranged for the 3-D model,the sparse macro-track transportacceleration method was proposed,which greatly reduces the calculation amount of theacceleration equation and obtains a
9、very good acceleration efficiency without reducing theaccuracy.The 2-D method of characteristics divides the angular space into M azimuthalangles and P pole angles.The 3-D method of characteristics generally uses a level sym-metric quadrature set.In 2-D geometry,the angle is relatively more dense at
10、 the poles.As a result,the 2-D macro-track transport acceleration can select a small number ofazimuth angles with reflection symmetry and one pole angle as the solution directions ofthe acceleration equation.The sparse macro-track transport acceleration method imple-mented in this paper used the sam
11、e angles as the method of characteristics.The numberof tracks for 3-D problems is an order of magnitude higher than for 2-D problems.Thatis why reducing the number of 3-D tracks is a promising method for reducing the amountof accelerated equation calculation.This paper applied the macro-track transp
12、ort acceler-ation method to 3-D whole core problem of high temperature gas-cooled pebble-bed reac-tor to solve the problem of irregular geometric numerical acceleration.The sparsity,macro-track length limit,and the maximum iteration step limit of the macro-track meth-od are the main factors that aff
13、ect the acceleration performance.The selection method ofacceleration parameters was studied through two benchmark problems,Excellent speed-up ratio can be obtained by setting the sparsity to 3-5,the macro-track length to around2.0 cm,the acceleration iteration steps limit to 20-60.The results of 3-D
14、 small lightwater reactor benchmark and simplified pebble-bed reactor core model show that about7 times the time speedup ratio and about corresponding 20 times the iterative step speed-up ratio can be achieved.Key words:high temperature gas-cooled pebble-bed reactor;method of characteristics;macro-t
15、rack transport acceleration method;neutron transport特征线法是实现球床高温气冷堆复杂几何全堆芯中子输运最有前景的方法。在数学上，特征线法是一种广泛应用的偏微分方程解析解法，2 0 世纪7 0 年代，AskewL2将偏微分方程特征线法和中子输运方程的离散纵标法、碰撞概率法结合起来，发展出了求解中子输运方程的特征线法。中子输运特征线法在求解区域按照离散纵标法产生的角度方向稠密布置大量的特征线，角通量在特征线上逐段解析求解，并通过积分得到网格标通量。中子输运方程的求解在特征线上完成，因而特征线法对求解区域几何和材料分布没有任何要求，非常适合于具有随
16、机分布燃料球的高温气冷堆中子输运求解。迭代次数多和计算速度慢是限制特征线法应用的主要因素。因此，特征线法通常作为组655件计算中的二维输运计算方法或者全堆芯计算2-D/1-D方法 3 的径向二维输运计算方法。WIMS4、C A SM O E5 和HELIOSL6 等组件计算软件均采用特征线法作为非均匀几何中子输运计算方法。为了解决特征线法计算速度慢的问题，CASMO、O p e n M O C 7 和DeCARTE81等程序采用粗网有限差分（CMFD)方法降低选代次数，可实现具有宏观方形规则区域的加速，显著减少了计算时间。基于2-D/1-D方法的MPACTC9和NECP-XL101等程序除了采
17、用CMFD方法进行数值加速外，还实现了在集群计算机上的大规模并行，能够在分钟量级给出单个算例的结果。为了突破数值加速方法中宏观规则几何的局限，Yamamotol提出了广义粗网再平衡(GCMR)方法，从理论上指出GCMR方法可应656用于非规则几何的粗网格，然而没有给出具体可行的方法。柴晓明等 12 提出了广义粗网有限差分(GCMFD)方法并引人了用于计算非规则几何耦合系数的宽度因子，但宽度因子无法从理论上直接给出而只能在迭代过程中逐渐逼近，此外，非规则几何区域网格边界面积和体积通常只能近似计算，这都减弱了GCMFD方法的加速效果。Smith131 提出了长特征线加速方法(macro-track
18、 transport acceleration method),不需要计算网格界面的中子流，具有和特征线法同样的几何处理能力，可真正实现任意几何特征线法加速。本文在长特征线加速方法中引入稀疏条数长特征线方法，可在不降低加速方程精度的前提下极大地缩短计算时间，解决三维几何加速方法计算量大的问题。1传统长特征线加速方法长特征线加速方法将连续几个线段合并为一个长特征线段，如图1所示，实线段表示长特征线。根据特征线法中的角通量，见式（1），可以得到长特征线上的源，见式（2）：Yig.m.k(so.out)=Vi.g.m.(so.in)etig.n.k+2.i.(sp.in)qg.m.k.1-eg.mk
19、.式中：q为长特征线上的源；pq为从起点p到终原子能科学技术第58 卷点q；亚为角通量；Q为网格源强；Z为中子总截面;t为光学厚度；i为网格编号；g为能群；m为离散角度方向；k为特征线编号；sin为线段的起点；Sout为线段的终点。光学厚度的计算公式为：Tg.m.k=(3)s=p.p+1,.q式中：s为特征线段编号；1为特征线段长度；s=p，p 十1，q为长特线从起点p到终点q穿过的所有特征线段。长特征线上的源可在特征线扫描过程中计算得到，然后可构造低阶的基于特征线法的加速方程。长特征线上的角通量为：(4)Tg.m.k.网格标通量用穿过网格的所有长特征线平均角通量修正为：d(m)i.g=0(0
20、)m.kEICM.K)igm,kEI(M.K)式中：为网格标通量，上角括号为迭代次数，0为启动长特征线加速的初始值，n为跳出加速时的迭代结果;I(M,K)为穿过网格i的所有角度方向的特征线段。(1)长特征线源用其穿过的所有网格源来修正为：ipg,m.k（n）(2)qm.g.k.风(5)=qm.g.k.风（0)Sp,P+1,s=.p+1,q(6)abj+1g,m2,k,平g,mj+1中i.j.8平g.mQi.j1.8QmQi.j.gQi+1,i.8,j-1Fig.1 Calculation method of grid scalar flux and macro-track source for
21、 macro-track transport acceleration methodj-1i-11a网格标通量计算方法示意图；b长特征线源计算方法示意图图1长特征线加速方法网格标通量和长特征线源计算方法+1i-11+1第3期郭建等：三维球床几何稀疏条数长特征线加速方法式中，Q为网格内的源。长特征线加速方法的计算流程如图2 所示。Smith13基于轻水堆1/4堆芯二维问题研究了长特征线方法的加速效果，在迭代步上得到了2 0 50 倍的加速效果。长特征线加速方法是一种非线性加速方法，本质上是对输运方程特征法的低阶近似，具有良好的稳定性，因此可以用于扩散非线性加速方法发散时问题的加速。长特征线加速方
22、法的缺点是计算量要明显高于其他加速方法，加速方程求解过于耗时导致其在时间上的加速效果弱于迭代步上的加速效果。开始初始化几何特征线布置和追踪长特征线初始化特征线扫描计算网格标通量、长特征线源式(2)计算kef和残差否收敛？是结束图2 长特征线加速方法计算流程Fig.2Flow chart of macro-tracktransport acceleration method2稀疏条数长特征线加速方法为了解决长特征线加速方法计算量大的问题，二维方法中通常采用稀疏角度的方式。二维特征线法将角度方向划分为M个方位角和P个极角，总的方向数为MP,方位角按照等角度间隔的方式划分，极角有多种划分方式，包括：
23、等角度、等权重和TY求积组 141等。三维特征线法一般采用旋转对称的求积组，常用的求积组有Carlson 求积组和 Lee求机组等 15。图3对比了二维和三维特征线法求积组的不同，可以看出二维方法中求积方向在两极相对稠密，因此二维特征线法可以选择具有反射对657称性的少量几个方位角和1个极角作为加速方程的求解方向，此时的网格标通量计算公式为：(7)m.kEI(Mrer.Pi.K)式中,I(Mrer,Pi,K)为选取的旋转对称方位角Mrer和极角Pi穿过网格i的特征线段。aa一二维特征线法求积组，16 个方位角，6 个极角；长特征线加速b一三维特征线法旋转对称求积组，S8图3二维和三维特征线法的
24、求积组对比计算长特征线源Fig.3Comparison of quadrature set式(6)for 2-D and 3-D method of characteristics长特征线扫描式(4)计算网格标通量式(5)计算ket和残差否收敛？是b三维问题相比二维问题特征线数量高一个数量级，本文提出了稀疏条数长特征线加速方法来降低三维加速方法的计算量，并引入了条数稀疏度概念，如图4所示。三维几何中，按照每t个相邻特征线取一的方式选择用于求解加速方程的特征线，这里的t即为条数稀疏度，此时网格标通量的计算公式为：(8)m.kEI(M,K,)采用稀疏条数长特征线加速方法相比稀疏角度的优点是不会降低
25、加速方程在角度离散上的精度。一维离散纵标法可以证明，S阶方程的精度相当于Pn-i阶，对于二维和三维问题虽然不能得到相同的结论，但求积方向越多方程的精度越高仍然成立。稀疏条数长特征线加速方法不降低特征线方向数，只在空间上减少特征线的数量，由于三维几何特征线布置相对二维几何是非常稠密的，减少特征线数量不会明显影响式（8)特征线数量积分的准确性。3加速方法参数选取条数稀疏度、长特征线长度和长特征线迭代步是影响加速效果的主要因素。本文基于小658原子能科学技术第58 卷a特征线方法，M=12;b型轻水堆基准题对长特征线加速方法的加速效果进行了分析。小型轻水堆基准题 16 是一个边长为50 cm的正方体
26、，计算模型划分为50 X50X50个边长为1cm的正方体小网格，三维特征线间距为0.2 cm，长特征线最大长度限制为2 cm，每次启动加速方法迭代50 步后退出。在进行加速参数研究前对特征线间距进行敏感性分析以避免特征线布置的疏密影响加速效果，计算结果如图5所示，有效增殖因数kefr偏差为相对特征线间距为0.0 5cm时的ker计算结果，可以看出keff偏差随着特征线间距降低而趋于稳定，本文特征线间距是根据kef偏差小于5 pcm的原则选取的。50r小型轻水堆面心立方球床4030上20F1000.0图5ker偏差随特征线间距变化Fig.5kef difference variation wit
27、h track spacing条数稀疏度对加速效果的影响如表1所列，条数稀疏度为1时的总扫描时间（包括特征线扫描和长特征线扫描)增加了接近1倍，简单计算可以得到长特征线每次扫描所需时间与特征线扫描所需时间相当。条数稀疏度不大于5时，随着条数稀疏度增加总扫描时间快速降低，稀疏角度长特征线加速方法，M=4；c 一利图4长特征线布置方法Fig.4Macro-track arrangement method迭代步和加速迭代步均不变，条数稀疏度进一步变大时加速效果开始出现明显减弱，迭代步和加速迭代步也大幅度增加。说明条数稀疏度大于5时，由于长特征线数量减少，式（8)的计算准确度会严重偏离真实情况。长特征
28、线长度对加速效果的影响如表2 所列，计算过程过中条数稀疏度为3，长特征线长度越长单次扫描所需时间越短，长特征线长度越长则低阶加速方程相比特征线方程近似越多，需要的迭代步越多，总的加速效果由以上两个方面因素共同决定，表现为总扫描时间先缩短后增加。长特征线迭代步用来控制每次启动加速方法后的迭代次数。长特征线迭代步对加速效果的影响如表3所列，计算过程中条数稀疏度为3,长特征线长度为2.0 cm。长特征线迭代步不大于6 0 步时，迭代步随长特征线迭代步增加而降低，之后即使进一步增加长特征线送代步，迭代步也保持不变。总扫描时间的最小值未出现在长特征线迭代步为6 0 步时，而是出现在2 0 步时，原因是总
29、的扫描时间中长特征线扫0.10.2特征线间距/cm稀疏条数长特征线加速方法，M=12，t=30.30.40.5描占很大一部分，长特征线迭代步的增多抵消了迭代步的降低。综合以上分析，条数稀疏度取35，长特征线长度取2.0 cm左右，长特征线迭代步取2 0 60步时可以获得比较好的加速效果。使用最佳参数时，三维长特征线加速方法在计算时间上可以获得7 倍左右的加速效果，此时在迭代步上可以获得2 0 倍左右的加速效果；只考虑迭代步最多时可以得到2 5倍的加速效果，与二维方法代步加速效果基本一致。第3期郭建等：三维球床几何稀疏条数长特征线加速方法表1长特征线加速方法加速效果随条数稀疏度的变化Table
30、1Macro-track transport acceleration performance with respect to sparsity加速总扫描条数稀疏度送代步277111211311411511612732828.9659特征线长特征线内存/MB送代步255355050315503274550295055028365502785600275516002.73714002726kelf扫描时间/s扫描时间/s13.6597.813.6167.113.684.713.655.613.641.614.836.739.583.2134.663.7时间s342.34611.4180.7398
31、.3369.2455.1951.5422.7398.30.9757490.9757630.9757630.9757630.9757630.9757630.975 7630.975 7630.975 763不收敛表2 长特征线加速方法加速效果随长特征线长度的变化Table2Macro-track transport acceleration performance with respect to length of macro-track长特征线送代步长度/cm2771.0121.5112.0112.5143.0153.5194.0204.5225.026加速送代步600550550700750
32、950100011001300特征线内存/MB扫描时间/s2553320330082950284828212.799279427772.766长特征线扫描时间/s14.8137.313.696.113.684.717.388.118.587.823.5104.224.7108.627.2115.532.1131.0总扫描时间/s342.34152.13109.6698.33105.39106.34127.69133.28142.67163.17keff0.9757490.9757630.9757630.975 7630.9757640.9757640.9757640.9757640.97576
33、40.975 764表3长特征线加速方法加速效果随长特征线选代步的变化Table 3Macro-track transport acceleration performance with respect to macro-track iteration step长特征线送代步送代步277103420203015401350116010701080109010加速送代步340400450520550600700800900特征线内存/MB扫描时间/s2553295029502950295029502950295029502950长特征线扫描时间/s42.056.724.763.618.570.7
34、16.180.813.684.712.492.812.4107.512.4122.312.4138.5总扫描时间/s342.3498.7588.3389.2196.8798.33105.11119.88134.61150.87ker0.9757490.9757620.9757630.975 7630.9757640.975 7640.9757640.9757640.9757640.9757646604球床几何加速效果使用球床几何模型对稀疏条数长特征线加速方法的加速效果进行了验证，受限于服务器的计算能力，构造了1个由17 2 个面心立方堆积的燃料球的简化模型，如图6 所示，模型边长为6 0 cm
35、，a反射层球床区域a一正方形球床几何模型；b一一球床区域燃料球堆积模型图6 三维面心立方规则堆积球床计算模型Fig.6Calculation model of 3-D pebble-bed geometrywith face-centered regularly packing pebbleTable 4 Macro-track transport acceleration performance in 3-D pebble-bed geometry条数长特征线稀疏度长度/cm一一41.051.042.052.043.053.044.054.045.055.0原子能科学技术第58 卷中心有一个
36、边长为31.46 cm的正方形空腔区域，燃料球半径为3cm,球床区域填充率为0.6 2 5，三维特征线间距为0.1cm。燃料球使用小型轻水堆基准题堆芯截面 16，反射层使用小型轻水堆基准题反射层截面，燃料球间隙为空腔，蒙特卡罗程序计算得到的kef=0.84310（0.0 0 0 16）。计算结果如表4所列，计算过程中长特征线迭代步固定为50 步，重点研究了条数稀疏度和长特征线长度对加速效果的影响。从计算结果可以看出，使用稀疏条数长特征线加速后，迭代步可以获得最多30 多倍的加速效果，特征线追踪时间可以获得最多7 倍左右的加速效果，此时迭代步加速效果为2 0 倍左右。三维面心立方堆积球床的计算证
37、明稀疏条数长特征线加速方法在包含空腔的非规则几何中同样可以获得比较好的加速效果。表4长特征线加速方法对三维球床几何的加速效果特征线长特征线送代步内存/MB4471665131948241879141871201830191 816191794221806221 788281 790291778总扫描扫描时间/s扫描时间/s一一10.2656.3718.9447.0911.0545.1615.7937.2715.0045.5515.0036.6517.3647.2917.3636.6122.1052.4022.8943.20kef时间/s352.8066.6366.0356.2153.0660.
38、5551.6564.6553.9774.5066.090.843340.843.330.843.330.843330.843.330.843330.843.330.843330.843330.843 330.843335结论为解决球床高温气冷堆三维特征线全堆芯输运问题的几何复杂性和迭代速度慢的问题，本文采用长特征线加速方法进行了加速研究，主要结论如下。1）将长特征线加速方法应用于三维非规则球床几何并提出了稀疏条数长特征线加速方法，极大地降低了加速方法的计算时间，显著改善了加速效果。2）稀疏条数长特征线加速方法在三维非规则球床几何扫描时间上可以获得7 倍左右的加速，此时迭代步上可以获得2 0 倍
39、左右的加速，仅考虑迭代步时最多可以得到30 多倍的加速效果，总体来说迭代步加速效果与二维特征线法相近。3）通过基准题研究，得到了所测试球床问题加速参数的选取范围。条数稀疏度取35、长特征线长度取2.0 cm左右、长特征线迭代步取2 0 6 0 步时可以获得比较好的加速效果。参考文献：1郭建，郭炯，王黎东，等针对球床的三维特征线第3期郭建等：三维球床几何稀疏条数长特征线加速方法方法几何建模研究 J原子能科学技术，2 0 18，9KOCHUNAS B,DOVVNAR T J,LIU Z.Par-52(3):412-419.allel 3-D method of characteristics in
40、 MPACTGUO Jian,GUO Jiong,WANG Lidong,et al.R.Chicago:American Nuclear Society,2013.Geometry modeling of 3D method of characteris-10马党伟,刘宙宇，赵晨，等NECP-X中角度与特tics for pebble-bedJ.Atomic Energy Science征线并行的研究与实现 原子能科学技术，and Technology,2018,52(3)：412-419(i n C h i-2018,52(3):405-411.nese).MA Dangwei,LIU Z
41、houyu,ZHAO Chen,et al.2ASKEW J.A characteristics formulation of the neu-tron transport equation in complicated geometriesR.England:United Kingdom Atomic EnergyAuthority,1972.3张广春，刘杰，YANGWONSIK.PROTEUS-MOC在TREAT试验堆稳态中子学计算中的应用 核技术，2 0 2 0，43（4）：0 40 0 0 8.ZHANG Guangchun,LIU Jie,YANG WONSIK.Applicatio
42、n of PROTEUS-MOC for thesteady state simulation of TREATJJ.NuclearTechniques,2020,43(4):040008(in Chinese).4ALAM S B,ALMUTAIRI B,KUMAR D,etal.Convergence studies using method of charac-teristics solver for the reduced-moderation waterreactor modelC/ANS Pacific Basin NuclearConference.USA:Dinesh Kuma
43、r,2018.5KNOTT D.Validation of the CASMO-4 trans-port solutionsCJ/Proc.Int.Meeting Mathe-matics Methods and Supercomputing in NuclearApplications.S.1.:s.n.J,19 9 3.6WEMPLE C A,GHEORGHIU H N M,STAMMLER R J J,et al.Recent advances inthe HELIOS-2 lattice physics codeCJ/Proceed-ings of the International
44、Conference on ReactorPhysics,Nuclear Power:A Sustainable Resource.Hungary:s.n.,2008.7BOYD W,SHANER S,LI L,et al.The Open-MOC method of characteristics neutral particletransport codeJJ.Annals of Nuclear Energy,2014,68:43-52.8JOO HG,CHO J Y,KIM K S,et al.Methodsand performance of a three-dimensional w
45、hole-core transport code DeCARTCJ/Physor 2004.United States:Ls.n.J,2004.661Development and application of angle and rayparallel in NECP-XJ.Atomic Energy Scienceand Technology,2018,52(3):405-411(in Chi-nese).11 YAMAMOTO A.Generalized coarse-meshrebalance method for acceleration of neutrontransport ca
46、lculations J.Nuclear Science andEngineering,2005,151(3)：2 7 4-2 8 2.12柴晓明，姚栋，王侃，等使用广义粗网有限差分方法加速中子输运特征线方法。计算物理，2010,27(4):541-547.CHAI Xiaoming,YAO Dong,WANG Kan,etal.Generalized coarse-mesh finite differenceacceleration for method of characteristics inunstructured meshesJ.Chinese Journal ofComputati
47、onal Physics，2 0 10,2 7（4):541-547(in Chinese).13 SMITH K S.Full-core,2-D,LWR core calcula-tion with CASMO-4ECJ/International Confer-ence on the New Frontiers of Nuclear Technolo-gy:Reactor Physics,Safety and High-perform-ance Computing.Korea:Ls.n.J,2002.14 YAMAMOTO A,TABUCHI M,SUGIMURAN,et al.Deriv
48、ation of optimum polar anglequadrature set for the method of characteristicsbased on approximation error for the BickleyfunctionJ.Journal of Nuclear Science andTechnology，2 0 0 7,44(2):12 9-136.15杜书华输运问题的计算机模拟M长沙：湖南科学技术出版社，19 8 9.16 TAKEDA T,IKEDA H.3-D neutron transportbenchmarksJ.Journal of Nuclear Science andTechnology，19 9 1,2 8(7):6 56-6 6 9.
