點堆中子動力學方程求解程序

我學的是核科學與技術專業，本科和碩士期間寫了很多程序，如計算堆芯中子動力學的，計算流體力學的。但多數僅針對具體問題，通用性有些不足。碩士修了一門仿真技術的課程，個人做業是求解點堆中子動力學程序，靠這個還拿了優秀。以爲這個程序還蠻通用的，對於同專業的搞個大做業，或用於畢設等或許有幫助。

詳細的推導、驗證過程、程序代碼見：https://github.com/ikheu/point_reactor

1 概述

計算中子通量密度的瞬變特性，對反應堆動力系統的運行安全分析與仿真而言十分重要。點堆中子動力學模型是反應堆動態學中最經常使用的方法。

點堆中子動力學模型描述了中子密度和緩發中子先驅核濃度隨時間變化的規律，基本方程如式（1）所示。因爲瞬發中子壽命與緩發中子壽命間存在着數量級的差異，耦合的點堆動力學微分方程組存在着剛性，這給數值求解帶來了必定困難。常規的顯式方法，如歐拉法、龍格庫塔法在計算該問題時會存在穩定性問題，這要求時間步長需取得很小，會帶來計算耗時以及大的累計偏差。針對點堆動力學方程組的剛性問題，發展了許多種處理方法，包括：（1）線性多步法，如GEAR算法及其改進形式；（2）基於積分方程的方法，如Hansen方法；（3）分段多項式逼近法，如埃爾米特插值型；（4）權重限制法；（5）指數基函數法等。 

$\left\{ \begin{gathered}
\frac{{dn(t)}}{{dt}} = \frac{{\rho (t) - \beta }}{\Lambda }n(t) + \sum\limits_{i = 1}^I {{\lambda _i}{C_i}(t)} + q \hfill \\
\frac{{d{C_i}(t)}}{{dt}} = \frac{{{\beta _i}}}{\Lambda }n(t) - {\lambda _i}{C_i}(t){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} i = 1,2,...,I \hfill \\ 
\end{gathered} \right.　　　　　　　　　　　　　　　　　(1)$

式中，${n(t)}$爲中子密度； 

         ${{C_i}(t)}$爲第i組緩發中子先驅核的濃度；

         ${\rho (t)}$爲反應性；

         $\Lambda$爲中子代時間；

         ${{\lambda _i}}$爲第i組緩發中子先驅核衰變常數；

         ${{\beta _i}}$爲第i組緩發中子份額；

         $\beta$爲緩發中子的總份額；

         $q$爲外加中子源項。

本文使用一階泰勒多項式積分方法求解點堆動力學方程，該方法能很好地解決方程的剛性問題，具備計算精度高、適用性好的特色。

2 一階泰勒多項式積分方法

詳細的推導過程有些繁瑣，這裏只說明其中兩個關鍵的步驟。

第一個是積分處理。合併式(1)，並在$\left[ {{t_n},{t_{n + 1}}} \right]$區間內積分，可得到式（2）：

$n({t_{n + 1}}) - n({t_n}) = \int_{{t_n}}^{{t_{n + 1}}} {\frac{{\rho (t)}}{\Lambda }n(t)}  - \sum\limits_{i = 1}^6 {[{C_i}({t_{n + 1}}) - {C_i}({t_n})]}  + q \cdot h　　　　　　　　　　　　　　(2)$

第二個是對中子密度做一階泰勒展開：

$n(\tau ) = n({t_{n + 1}}) + {n^,}({t_{n + 1}})(\tau  - {t_{n + 1}})　　　　　　　　　　　　　　　　　　　　(3)$

通過一系列的推導，可得到中子密度的表達式：

$n({t_{n + 1}}) = \frac{{n({t_n}) + q \cdot h + \sum\limits_{i = 1}^6 {{C_i}({t_n})}  - \sum\limits_{i = 1}^6 {\exp ( - {\lambda _i}h){C_i}({t_n})}  + \frac{{({F_2} - \sum\limits_{i = 1}^6 {\frac{{{\beta _i}}}{\Lambda }{G_{2,i}})} (\sum\limits_{i = 1}^6 {{\lambda _i}\exp ( - {\lambda _i}h){C_i}({t_n})}  + q)}}{{(1 - \sum\limits_{i = 1}^6 {\frac{{{\lambda _i}{\beta _i}}}{\Lambda }{G_{2,i}}} )}}}}{{1 - {F_1} + \sum\limits_{i = 1}^6 {\frac{{{\beta _i}}}{\Lambda }{G_{1,i}}}  - ({F_2} - \sum\limits_{i = 1}^6 {\frac{{{\beta _i}}}{\Lambda }{G_{2,i}})} (\frac{{\rho ({t_{n + 1}}) - \beta }}{\Lambda } + \sum\limits_{i = 1}^6 {\frac{{{\lambda _i}{\beta _i}}}{\Lambda }{G_{1,i}}} )/(1 - \sum\limits_{i = 1}^6 {\frac{{{\lambda _i}{\beta _i}}}{\Lambda }{G_{2,i}}} )}}　　　　　　　　　　　　　(4)$

截斷偏差是該方法的主要偏差來源，可考慮更高階數的泰勒展開，固然這會使模型更加複雜（上式已經很複雜了），求解過程也不夠簡便。

下圖爲本方法的計算流程圖。

首先給出初始條件，即${t_n} = 0$時的${n(t)}$、 ${{C_i}(t)}$，肯定已知項${\rho (t)}$、 ${{\lambda _i}}$、${{\beta _i}}$、$q$等；而後給定計算條件：計算時長與時間步長；逐時刻依次計算下一時刻${t_{n + 1}} = {t_n} + \Delta t$ 時的$n({t_{n + 1}})$、${n^,}({t_{n + 1}})$與${C_i}({t_{n + 1}})$，直至達到設置的計算時間。

認爲反應堆初始狀態是穩定的，給定初始的中子密度$n(0)$，根據式（1）可得：

$\frac{{d{C_i}(t)}}{{dt}} = 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}  \Rightarrow {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {C_i}(0) = \frac{{{\beta _i}}}{{{\lambda _i}\Lambda }}n(0)　　　　　　　　　　　　　　　　　　　　(5)$

使用老掉牙的Fortran90程序編制的計算程序，博客裏甚至沒有這個語言的着色，IT行業估計也沒人用這個，因此代碼就不貼了。能夠在文章開頭的github連接裏找到。 

3 計算結果

這裏給出了四種反應性引入，包括正反應性階躍、負反應性階躍、反應性線性變化、反應性正弦變化情形下，典型熱中子堆的中子密度響應曲線：

 github連接中有個「點堆模型與求解.pdf」的文件，能夠在其中找到計算結果與解析解的比較和分析。計算程序還適用於快中子堆的求解。