全隱式龍格庫塔法求解點堆動力學方程

王偉吉，葉金亮，方成躍

（1.海軍裝備部，北京 100071；2.中國艦船研究設計中心，湖北 武漢430064）

在反應堆動力學和反應堆安全分析中，必 須利用中子動力學方程對中子增殖進行計算。由于點堆中子動力學方程組是一個耦合的剛性一階微分方程組，求解困難；而且實際動力堆運行時，反應性溫度反饋系數不能被忽略，這給反應堆安全分析和運行現場的實時計算等帶來了很大的難度。因而，尋求一種高效、實時的數值方法來求解點堆動力學方程組就很有必要。圍繞著點堆動力學方程組的剛性問題，國內外學者發展了眾多的數值解法，主要有（1）線性多步法，比如 GEAR法、改進 GEAR 法［1］；（2）分段多項式逼近法，如拉格朗日插值型和埃爾米特型［2，3］；（3）基于積分方程的各種方法，如變量隱含積分法、Hansen方法［4］；（4）權重限制法；（5）外推和修勻技巧；（6）有限元法［5］；（7）泰勒級數展開法［6］；（8）冪級數法［7］；（9）指數基函數法［8］。

以上介紹的數值方法中，GEAR法和分段多項式逼近法計算精度較高，但GEAR方法較費時，而分段多項式法計算工作量較小，較省時間，但有時需把步長取得非常小，否則會帶來較大的誤差，這也導致計算時間拉長［9］。改進GEAR方法穩定域要比GEAR方法大，對步長要求不是那么苛刻，但是計算量較大，比較費時。

另外，冪級數法、指數基函數法計算速度較快、精度較高，適用性較強。

該文研究的基于高斯－勒讓特求積公式節點的全隱式龍格庫塔法，克服了點堆中子動力學方程組的剛性，保證了數值計算結果的穩定性和精確性，而且計算速度較快，計算過程簡捷，易于編程實現。

1 基于高斯－勒讓特求積公式節點的全隱式龍格庫塔法

不考慮外加中子源的S組緩發中子的點堆中子動力學方程組為

式中，i＝1，2，Λ，S；N（t）為中子密度，cm－3；t為時間，s；ρ（t）為反應性；β為緩發中子總份額；βi為第i組緩發中子份額；λi為第i組緩發中子先驅核衰變常數，s－1；Λ 為中子代時間，s；Ci（t）為第i組緩發中子先驅核平均濃度，cm－3；S為緩發中子先驅核組數。

將上述方程組表示為矩陣形式，有：

其中，Y（t）＝［N（t），C1（t），Λ，CS（t）］T

F（t）表示為如下（S＋1）×（S＋1）階矩陣：

由于F（t）為非常系數矩陣，且是病態的，因此式（3）是剛性方程。在熱堆線性反應性引入、快堆階躍反應性加入、快堆負線性反應性加入等輸入條件下，F（t）是強剛性的，求解式（3）非常困難。式（3）的全隱式龍格庫塔法求解形式為

式中，h為計算步長，單位s；ν是GLFIRK方法的節點數，δi為GLFIRK方法的GL節點，γi為權系數，B＝（βij）（i，j＝1，Λ，ν）為 GLFIRK方法的系數矩陣，f（t，Y）代表式（3）右端式子。令與式（7）聯立后可得

式（8）中，fj＝（tn＋δjh，pj）

這樣，對式（8）進行推理、化簡可得

其中，

在式（9）、式（10）、式（11）中，r為矩陣（B－1）T的特征值，Iν為ν階單位矩陣，l為迭代次數，e＝（1，Λ，1）T∈Rν，P＝（p1，Λ，pi，Λ，pv），D＝（p1，Λ，pi，Λ，pv）。

因此，由式（9）可迭代計算出Pl＋1，那么

將式（12）代入式（6）即可求得每步的數值解。

表1給出三級六階GLFIRK法系數，具體推導過程可參閱文獻［10］。另外，文獻［10］指出GLFIRK是B穩定的，并給出了證明過程，這里不再贅述。

該文以三級六階GLFIRK方法求解點堆動力學方程，并用MATLAB編制了計算程序，程序簡單實用，計算速度快。

表1 三級六階GLFIRK方法系數Table 1Three stages six orders GLFIRK coefficients

2 計算實例

例1［11］計算階躍正反應性輸入后熱堆內中子密度N（t）隨時間的變化。熱堆參數為βi：2.66×10－4；1.491×10－3；1.316×10－3；2.849×10－4；8.96×10－4；1.82×10－4；λi：0.012 7；0.031 7；0.115；0.311；1.4；3.87；中子代時間Λ＝2.0×10－5緩發中子總份額β＝0.007，輸入階躍正反應性ρ＝0.003；N（0）＝1.0cm－3；計算1s內中子密度隨時間的變化。計算結果如表2所示。表2第二列為解析解，其余三列為不同時間步長下由GLFIRK方法求解的數值解，最后一行表示計算得到收斂解CPU所需總時間（Intel Core2Duo CPU E7500＠2.93GHz）。由表3可以看出，步長h＝0.005s時數值解與精確解吻合得很好，基本保持七位有效數字一致，而h＝0.1s時較差。另外，隨著時間步長的增加，計算得到收斂解CPU所需總時間在快速減少。步長取得越小，求解得到的數值解精度越高，越接近解析解，與此同時，所需計算時間就越長。因此，在滿足一定精度要求的條件下，可適當增加時間步長，這樣可以縮短計算時間，以滿足實時性要求。

表2 階躍正反應性輸入熱堆中子密度N（t）隨時間變化Table 2 Neutron densityN（t）v.s.Time with positive step reactivity insertion in thermal reactor

例2計算大正反應性輸入后，在瞬發臨界的情況下，熱堆內中子密度隨時間的變化。輸入的 正 反 應 性ρ＝β＝0.007，N（0）＝1.0cm－3，計算1s內中子密度隨時間的變化，熱堆參數同例1。計算結果如表3、表4所示。表3第二列為解析解，其余三列為不同時間步長下由GLFIRK方法求解的數值解。表4為不同時間步長下數值解與精確解之間的相對誤差。相對誤差100%，由表4可知，隨著計算步長的增大，相對誤差在增加，而計算時間在縮短。另外，對于大正反應性輸入，點堆動力學方程剛性已不明顯，但GLFIRK方法仍能在大步長條件下（h＝0.1s）將計算精度維持在10－4，由此可見GLFIRK方法求解非剛性微分方程時計算步長的選取范圍較之解剛性微分方程要寬些。

表3 大正反應性輸入熱堆中子密度N（t）隨時間變化Table 3 Neutron densityN（t）v.s.Time with great positive reactivity insertion in thermal reactor

表4 不同時間步長下的相對誤差Table 4 Relative Error in different Time step

例3計算階躍負反應性輸入后，堆內中子密度隨時間的變化。輸入的負反應性ρ＝0.007，N（0）＝1.0cm－3，計算1s內中子密度隨時間的變化，熱堆參數同例1。計算結果見表5。表5第二列為精確解，其余三列為在不同計算步長下由GLFIRK計算得到的數值解，最后一行為不同計算步長下GLFIRK計算所花時間。由表5可知，在步長較小時，數值解與精確解吻合的相當好，但在大步長下數值解較之精確解出現很大偏差，有些數值甚至是錯誤的。可見，在階躍負反應性輸入時，點堆方程剛性較強，這時GLFIRK計算步長宜取較小值，以滿足精度要求，否則會出現較大偏差甚至錯誤數值。

表5 階躍負反應性輸入熱堆中子密度N（t）隨時間變化Table 5 Neutron density N（t）v.s.Time with negative step reactivity insertion in thermal reactor

例4計算線性反應性輸入后，熱堆中子密度隨時間變化。輸入的反應性ρ＝0.000 7t，N（0）＝1.0cm－3，計算10s內中子密度隨時間的變化。熱堆參數同例1。計算結果見表6、圖1。表6第二列為精確解，其余三列為不同時間步長下由GLFIRK方法求解的數值解，最后一行為不同計算步長下GLFIRK所花時間。圖1中相對誤差的對數值是取10為底的對數值。

由表6、圖1可知，在步長較小時，數值解非常接近精確解，而當計算步長逐步增大時，相對誤差快速增加，嚴重偏離解析解，此時計算發散。因此，在線性反應性輸入時，點堆方程剛性很強，此時，GLFIRK宜取較小的時間步長以保證計算精度，但這樣就增加了計算時間，實時性較差。

表6 線性反應性輸入熱堆中子密度N（t）隨時間變化Table 6 Neutron density N（t）v.s.Time with positive ramp reactivity insertion in thermal reactor

續表

圖1 不同時間步長下的相對誤差曲線Fig.1 Relative Error Curve in different time step

例5計算正弦反應性輸入后，熱堆中子密度隨時間的變化。ρ＝ρ0sin（πt／T），T為輸入正弦反應性的周期，s，ρ0為反應性幅度，熱堆參數同例1。計算結果見表7、圖2。由表7可以看出，GLFIRK方法能準確地計算出峰值及其對應的時間，而且當其步長為Hermite方法的100倍時仍能保證5～7位有效數字相同。圖2為ρ0＝0.2β，T＝5s時的正弦反應性輸入后中子密度隨時間的變化規律。從圖2可以看出，中子密度隨時間的增長不停地振蕩，振幅逐漸加大，該結果與文獻［13］得出的結論一致。

表7 正弦反應性輸入熱堆中子密度N（t）峰值及其對應時間Table 7 Peak of Neutron DensityN（t）and Corresponding Time with Sinusoidal Reactivity Insertion

圖2 正弦反應性輸入后中子密度隨時間的變化Fig.2 Neutron Density v.s.Time with Sinusoid Reactivity Insertion

3 結論

該文用GLFIRK方法在熱堆、快堆的典型反應性輸入條件下對點堆動力學方程進行求解，求解結果表明該方法的計算精度很高、計算速度較快、適應能力較強，滿足一定的工程應用要求。主要有以下幾點結論。

（1）對于剛性不顯著或者非剛性情況，如熱堆階躍小正反應性引入、階躍大正反應性引入，GLFIRK方法能滿足計算精度和實時性要求，從例1、例2的計算結果可以看出。

（2）對于剛性顯著的情況，如熱堆階躍負反應性引入，GLFIRK方法仍能滿足計算精度和實時性要求，從算例3的計算結果可以看出。此時，計算步長不宜取大值，以h＝0.001s為宜，這時可兼顧計算精度和實時性兩方面要求。

（3）對于剛性很強的情況，如熱堆線性反應性引入，GLFIRK方法不能同時兼顧計算精度和實時性要求，從例4的計算結果可以看出。此時，取步長h＝0.001s計算可以得到滿足工程應用計算精度的數值，但實時性要差些。

（4）對于復雜反應性引入情況，如正弦反應性引入，GLFIRK方法可以滿足計算精度和實時性要求。從例5計算結果可以看出。

綜上所述，GLFIRK方法是計算穩定、數值精度較高、適應力較好，可滿足一定的工程應用要求。

［1］ 肖紅，王侃.求解點堆動態方程的IGEAR方法［J］.核科學與工程，2008，38（2）.

［2］ J.P.Hennart Piecewise Polynomial Approximations for Nuclear Reactor Point and Space Kinetics［J］.Nuclear Science and Engineering，1977，64：875－901.

［3］ Kueng Yeh Polynomial Approach to Reactor Kinetics Equations［J］.Nuclear Science and Engineering，1978，66：235－242.

［4］ Hansen K.F.，Koen B.V..Little W.W .Stable Numerical Solutions of the Reactor Kinectics Equations［J］.Nuclear Science and Engineering，1965，22：51－59.

［5］ 趙志遠.用有限元法解點反應堆動力學方程［J］.原子能科學與技術，1981，15（6）：656－663.

［6］ 陳昌友.一個新的求解點堆中子動力學方程組的數值方法［J］.核科學與工程，1998，18（4）：364－370.

［7］ Basken J，Lewins J.D.Power series Solutions of the Reactor Kinetics Equations［J］.Nuclear Science and Engineering，，1996，122：407－416.

［8］ 黎浩峰，陳文振，朱倩，羅磊 .點堆中子動力學方程的指數基函數法求解［J］.核動力工程，2009，30（4）.

［9］ 經滎清，李君利，羅經宇 .點動態方程幾種解法的比較［J］.清華學報，1995，35（3）：7－12.

［10］ 李慶楊 .常微分方程數值解法［M］.北京：清華大學出版社，1990：87－120.

［11］ 胡大璞 .克服點堆中子動力學方程剛性的新方法—端點浮動法［J］.核動力工程，1993，14（2）.

［12］ 濮繼龍 .點堆動態方程的半解析解［J］.核動力工程，1983，4（1）.

［13］ Peinetti F.，Nicolino C..Ravetto P.Kinetics of a point reactor in the presence of reactivity Oscillations［J］.Ann Nucl Energy，2006，33：1189－1195.