我學的是核科學與技術專業，本科和碩士期間寫了不少程序，如計算堆芯中子動力學的，計算流體力學的。但多數僅針對具體問題，通用性有些不足。碩士修了一門仿真技術的課程，我的作業是求解點堆中子動力學程序，靠這個還拿了優秀。覺得這個程序還蠻通用的，對於同專業的搞個大作業，或用於畢設等或許有幫助。
詳細的推導、驗證過程、程序代碼見：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”的文件，可以在其中找到計算結果與解析解的比較和分析。計算程序還適用於快中子堆的求解。