全隐式龙格库塔法求解点堆动力学方程
王伟吉，叶金亮²，方成跃2
（1.海军装备部，北京100071；2.中国报船研究设计中心，医北武汉430064）
要：强剧性问题时数值求解点滩中子动力学方程组的难点之一.该文用基于高斯-勒让特求积公式节点的全隐式龙格库塔法（简称GLFIRK）求解点堆动力学方程组，该方法是B稳定的，面且计算精度高，对于E级GLFIRK，其计算精度为2E阶，该文在阶跃、线性和正弦等不同反应性加人条件下对点堆动 力学方程组进行了计算，计算结果表明，该方法计算精度高、计算速度较快、适应能力较好，可满足一定的工程应用要求.
关键调：高斯-勒让特求积公式节点：全隐式龙格库塔法：点堆动力学方程
中国分类号：TL11 文章标志码：A
文章编号：0258-0918（2014)03-0289-07
TheFull ImplicitRunge-KuttaMethodofSolving Point Reactor Kinetics Equations
WANG Wei-ji' YE Jin-liang* FANG Cheng-yue*
(1. Armament Department of Navy Beijing L0o071 China;2. China Ship Study and Design Centre Wuhan of Habei Prov 430064 China)
Abstraet; The strong stiff problem is one dificult point of numerical solving of PointReactor Kinetics Equstions. The paper solves Point Reactor Kinetics Equations usingFull Implicit Runge-Kutta method based on Gauss-Legendre Quadrature nodes (GLFIRK). The method is B stability and its puting precision is high the S stagesGLFIRK method’s puting precision is 2E. The paper calculates the Point ReactorKinetics Equations in step ramp. and sinusoidal reactivity insertion. And the resultshows its puting accuracy is high puting speed is good and applicability isstrong. which meet some demand of practical engineering.
Key words: Gauss-Legendre Quadrature Nodes (GL); Full Implicit Runge-kutta; PointReactor Kinetics
在反应堆动力学和反应堆安全分析中，必须利用中子动力学方程对中子增殖进行计算.
由于点堆中子动力学方程组是一个耦合的刚性一阶微分方程组，求解困难；而且实际动力堆运行时，反应性温度系数不能被忽略，这给反 应堆安全分析和运行现场的实时计算等带来了很大的难度.因而，寻求一种高效、实时的数值方法来求解点堆动力学方程组就很有必要，围绕着点堆动力学方程组的刚性问题，国内外学者发展了众多的数值解法，主要有（1）线性多步 法，比如GEAR法、改进GEAR法；（2）分段多项式逼近法，如拉格朗日插值型和埃尔米特型²-；（3）基于积分方程的各种方法，如变量隐含积分法、Hansen方法：（4）权重限制法； （5）外推和修匀技巧：（6）有限元法；（7）泰勒级数展开法；（8）幂级数法”；（9）指数基函数法
C.（t）为第i组缓发中子先驱核平均浓度，cm；S为缓发中子先驱核组数.
将上述方程组表示为矩阵形式，有：
(3)
(4)
其中，Y（t）=[N（t)，C（t)，A，Cs（t)]
F（t）表示为如下（S1）×（S1）阶矩阵：
(5)
以上介绍的数值方法中，GEAR法和分段多项式逼近法计算精度较高，但GEAR方法较 费时，面分段多项式法计算工作量较小，较省时间，但有时需把步长取得非常小，否则会带来较大的误差，这也导致计算时间拉长口.改进GEAR方法稳定城要比GEAR方法大，对步长要求不是那么苛刻，但是计算量较大，比较费时.
由于F（r)为非常系数矩阵，且是病态的，因此 式（3）是刚性方程.在热堆线性反应性引人、快维阶联反应性加人、快堆负线性反应性加人等输人条件下，F（r）是强刚性的，求解式（3）非常困难，式（3）的全隐式龙格库塔法求 解形式为
另外，级数法、指数基函数法计算速度较快、精度较高，适用性较强.
(6)
该文研究的基于高斯-勒让特求积公式节学方程组的刚性，保证了数值计算结果的稳定 点的全隐式龙格库塔法，克服了点堆中子动力性和精确性，而且计算速度较快，计算过程简捷，易于编程实现.
式中，A为计算步长，单位sv是GLFIRK方法的节点数，为GLFIRK方法的GL节点，Y 为权系数，B=（β）（i，j=1，A，v）为GLFIRK方法的系数矩阵，f（t，Y）代表式（3）右端式子.
1基于高斯-勒让特求积公式节点 的全隐式龙格库塔法
令p，=Y hB K，，与式（7）联立后可得
不考虑外加中子源的S组缓发中子的点堆中子动力学方程组为
(8)
1式（8）中.f =（r.8h，p;）
这样，对式（8）进行推理、化简可得
(2)
(6)
式中，i=1.2，A.S；N（t）为中子密度，cm；t 为时间，s：p（t）为反应性：3为缓发中子总份额；3.为第i组缓发中子份额；入：为第：组缓发中子先驱核衰变常数，s；A为中子代时间，s；
其中，
(10)
CPU所需总时间（Intel Core2Duo CPU E7500@2.93GHz）.由表3可以看出，步长h=0.005s时数值解与精确解吻合得很好，基本保持七位有效数字一致，面h=0.1s时较差.CPU所需总时间在快速减少，步长取得越 另外，随着时间步长的增加，计算得到收敛解小，求解得到的数值解精度越高，越接近解析解，与此同时，所需计算时间就越长，因此，在满足一定精度要求的条件下，可适当增加 时间步长，这样可以缩短计算时间，以满足实时性要求.
(11)
在式（9）、式（10）、式（11）中，r为矩阵（B-1）的特征值，1.为v阶单位矩阵，为选代次数，e=(1.，1)∈R P=（p，A，pA，p) D=(p A，pA，p).
因此，由式（9）可选代计算出P1，那么
K，=f(t 8;h，p)=（Fp;)｜=a（12)将式（12）代人式（6）即可求得每步的数值解.
表2阶联正反应性输入热堆中子密度N（r）随时间变化Table 2 Neutron densityN(r) v. s. Time with positivestep reactivity insertin in thermal reactor
表1给出三级六阶GLFIRK法系数，具体推导过程可参阅文献[10].另外，文献[10]指 出GLFIRK是B稳定的，并给出了证明过程，这里不再赞述.
t/s 0.1 1. 799 528 h=0.005.A=0.01 sA=0.1 s 1.799 528 2.026 1500.2 1. 851 268 1.851 268 1.851 268 1. 782 2660.3 1.900 405 1 960 405 1 900 405 1 921 4150.4 1.947 593 1 947 593 1. 947 593 1.941 1960.5 1. 993 313 1 993 313 1. 993 313 1. 995 2610.6 2.037 922 2.037 922 2.037 922 2.037 3290.7 2.081 692 2.081 692 2.081 692 2. 081 8730.8 2. 124 832 2.124 832 2. 124 832 2. 124 7770.9 2.167 504 2.167 504 2.167 504 2.167 5211.0 2. 209 841 2.209 840 2.209 840 2.209 835ssd - 0.047 0.031 0 016
该文以三级六阶GLFIRK方法求解点堆动力学方程，井用MATLAB编制了计算程序，程序简单实用，计算速度快.
Table 1 Three stages six orders GLFIRK coefficients 表1三级大阶GLFIRK方法系数
2 10 15 30 24 36 24÷ 2 10 36 30 1518
例2计算大正反应性输人后，在瞬发临界的情况下，热难内中子密度随时间的变化.输人的正反应性p-β=0.007，N（0）- 1.0cm，计算1s内中子密度随时间的变化，热堆参数同例1.计算结果如表3、表4所示.表3第二列为解析解，其余三列为不同时间步长时间步长下数值解与精确解之间的相对误差. 下由GLFIRK方法求解的数值解.表4为不同相对误差|E，| 精确解一数值解|精确解 ×100%，由表4可知，随着计算步长的增大，相对误差在增 加，而计算时间在缩短.另外，对于大正反应性输人，点堆动力学方程刚性已不明显，但GLFIRK方
2计算实例
例1计算阶跃正反应性输人后热堆内中子密度N（t）随时间的变化.热堆参数为B：2.66×10;1. 491× 10;1.316× 10;2. 849 × 10; 8. 96 × 10; 1. 82 × 10; 1代时间A=2.0×10-缓发中子总份额β= 0.012 7;0.0317;0.115;0.311;1.4;3.87;中子0.007.输人阶联正反应性p=0.003N（0）=1.0cm"；计算1s内中子密度随时间的变化.计算结果如表2所示.表2第二列为解析解， 其余三列为不同时间步长下由GLFIRK方法求解的数值解，最后一行表示计算得到收做解
法仍能在大步长条件下（h=0.1s）将计算精度维持在10，由此可见GLFIRK方法求解非刚性微分方程时计算步长的选取范围较之解刚性微分方程要宽些.
表3大正反应性输入热堆中子密度N[r）随时同变化 Table 3 Neutren densityN() v. s. Time with greatpositive reactivity insertion in thermal rector
t/s 精确新 A=0.005 s A =0.01 s A=0.1 s0.1 4. 522 303 E1 4. 522 304 E1 4.522 304 E1 4.522 483 E10.2 1. 597 257 1.597 258 1. 597 258 1. 597 3585 188 961 E2 5. 188 965 E2 5. 188 965 F2 5. 189 439 E20.3 E2 1. 667 288 E2 E2 E20.4 1. 667 286 E1 E3 1 667 288 E) 1. 667 490 E30.5 5 345 879 E) 5.345 888 5 345 888 5. 346 699 E1. 713 190 1. 713 193 1. 713 193 1. 713 5050.6 5.489 493 E4 5. 489 506 E 5. 489 507 E4 5. 490 673 EI0.7 Et Et E E40.8 1. 758 905 ES 1. 758 910 ES 1. 758 910 ES 1. 759 337 ES0.9 5. 635 701 5.635 715 5.635 716 5.637 2551.805 726 ES 1. 805 732 ES 1. 805 732 ES 1. 806 280 ES1.0 E6 E5 E6 E6Cpu.s 0.109 0.062 0.016
表4不同时间步长下的相对误差
Table 4
Relative Error in different Time step
0.1 t/s A0.005s 3 095E-7 A=0.01 s 3. 206E-7 3. 979E-6 =0. 1 s0.2 4. 238E-7 4.351E-7 6. 302E-50.3 6 893E-7 7 001E-7 9. 218E-50.4 0.5 1. 001E-6 1. 596E-6 1. 012E-6 1 607E-6 1.240E-40.6 1. 982E-6 1.993E-6 1. 840E-40.8 0.7 2. 457E-6 2 691E-6 2.468E-6 2. 702E-6 2. 148E-4 2. 454E-40.9 2.568E-6 2.579E-6 2. 757E-41.0 3. 106E-6 3.117E-6 3.065E-4
例3计算阶跃负反应性输人后，堆内中子密度随时间的变化.输人的负反应性p=
0.007，N（0）=1.0cm²，计算1s内中子密度随时间的变化，热堆参数同例1.计算结果见表5.表5第二列为精确解，其余三列为在不同计算步长下由GLFIRK计算得到的数值解， 最后一行为不同计算步长下GLFIRK计算所花时间.由表5可知，在步长较小时，数值解与精确解吻合的相当好，但在大步长下数值解较之精确解出现很大偏差，有些数值甚至是错误 的，可见，在阶跃负反应性输入时，点堆方程刚性较强，这时GLFIRK计算步长宜取较小值，以满足精度要求，否则会出现较大偏差甚至错误数值.
表5阶联负反应性输入热堆 中子密度N（）随时间变化
Table 5 Neutron density N(t) v. s. Time withnegative step reactivity insertion in thermal reactor
t/s 精确好 =0. 005 s =0.01 s A=0. 1 s0.1 4. 900 708 E-1 4 900 696 E-1 4. 900 696 E-1 1. 353 332 E-10.2 4. 809 743 E-1 4. 809 732 E-1 4. 809 732 E-1 7.328 059 E-10.3 4. 727 747 4. 727 736 4. 727 734 2.939 9274. 652 903 E-1 E-1 E-1 5. 922 079 E-10.4 E-1 4. 652 893 E-1 4. 652 893 E-1 E-10.5 4. 583 873 E-1 4. 583 863 E-1 4. 583 863 E-1 3. 682 849 E-10 6 4 519 650 4 519 639 4. 519 640 5. 159 2834. 459 467 E-1 4. 459 457 E-1 4. 459 457 E-1 4. 005 364 E-10.7 E-1 E-1 E-1 E-10.8 4.402 732 E-1 4. 402 722 E-1 4 . 402 722 E-1 4 725 090 g-10.9 4. 348 978 4. 348 968 4.348 968 4.120 1154. 297 830 E-1 4. 297 820 E-1 4.297 820 E-1 4 460 286 E-11. 0 E-1 E-1 E-1 E-1Cpu8 0. 062 0. 032 0 016
例4计算线性反应性输人后，热堆中子密度随时间变化.输人的反应性p=0.0007t，N（0）=1.0cm，计算10s内中子密度随时 间的变化.热堆参数同例1.计算结果见表6、图1.表6第二列为精确解，其余三列为不同时间步长下由GLFIRK方法求解的数值解，最
续表
后一行为不同计算步长下GLFIRK所花时间，图1中相对误差的对数值是取10为底的对数值.
t/s A-1.0X b- A= =10 s 0.001 s 0.01 * 0.1 s3 1. 668 750 1.668 792 1 669 170 1. 672 957 1.711 684 E0 E) E0 E0 E04 2. 228 442 2.228 516 2.229 181 2.235 849 2.304 4513. 277 2513.277 398 3.278 720 3. 291 986 3. 429 614 E) 0 ED E0 E)5 E) B) E0 E0 036 5 582 063 5 582 405 5.585 586 5 617 544 5 953 368 EO EO E0 E01. 216 798 1.216 915 1.217 968 1.228 569 1.342 455 E0E1 E1 E1 E1 E18 4. 278 625 4.279 336 4.285 701 4 350 074 5.073 565 E1 E1 E1 E1 E19 4.875 186 4.876 958 4 892 815 5. 055 040 7.120 6664. 511 585 4. 517 211 4. 567 723 3 E2 E2 15 108 965 16. 986 083 E2 E210 ES E5 ES ES E5Cpu. 26. 328 3.375 0.39 0. 047s
非常接近精确解，面当计算步长逐步增大时，相 由表6、图1可知，在步长较小时，数值解对误差快速增加，严重偏离解析解，此时计算发散.因此，在线性反应性输人时，点堆方程刚性保证计算精度，但这样就增加了计算时间，实时 很强，此时，GLFIRK宜取较小的时间步长以性较差.
表6线性反应性验入热维中子密度N（r）随时间变化Table 6 Neutron density N(r) v.s. Time with positiveramp reactivity insertion in thermal reactor
/s确解 A=1.0X 10* 0. 001 A= 0.01 h= 0.1 s A=1. 130 833 1. 130 849 1.130 997 1. 132 479 1. 147 5061 E0 E0 E0 E0 E92 1. 338 200 1. 338 226 1. 338 458 1. 340 782 1. 364 445 EO E0 Eo E0 E0
图1不网时间步长下的相对误差曲线
Fig. 1Relative Error Curve in different time step
的100倍时仍能保证5~7位有效数字相同.图2为p=0.23，T=5s时的正弦反应性输人后中子密度随时间的变化规律.从图2可以看 出，中子密度随时间的增长不停地振荡，振幅逐渐加大，该结果与文献[13]得出的结论一致.
例5计算正弦反应性输人后，热堆中子密度随时间的变化.p=psin（πt/T）.T为输人正弦反应性的周期，sp为反应性幅度，热堆参数同例1.计算结果见表7、图2.由表7可 以看出，GLFIRK方法能准确地计算出峰值及其对应的时间，而且当其步长为Hermite方法
