技術領域：
本發(fā)明涉及核反應堆設計和反應堆物理計算領域，具體涉及一種利用耦合蒙特卡羅方法模擬ads系統(tǒng)瞬態(tài)問題的方法。
背景技術：
：
加速器驅(qū)動次臨界反應堆系統(tǒng)(ads)是一種新型核能反應堆系統(tǒng)，通過加速器產(chǎn)生高能質(zhì)子流，轟擊散裂靶產(chǎn)生散裂中子，作為外中子源驅(qū)動次臨界堆芯，嬗變其中的核廢料。該系統(tǒng)有中子通量密度高和能譜較硬的特點，其u/pu核素裝載份額的減少會導致堆芯的有效緩發(fā)中子份額減少，也會導致燃料的多普勒溫度反饋效應急劇減弱，降低了超功率過程中的負反饋機制。該型反應堆的瞬態(tài)特性特殊，需要對其進行精確的計算以保證堆芯設計的安全。
現(xiàn)有的反應堆瞬態(tài)分析方法中，中子動力學的計算為確定論方法，將時空中子動力學方程經(jīng)全隱式向后差分后離散為求解各時間步上的固定源形式的中子學方程，得到每個時間步上的通量。現(xiàn)有方法中，多采用求解擴散方程的形式得到通量，精確度上還有待提高。
蒙特卡羅方法利用模擬粒子在堆內(nèi)的行為進行計算，最后統(tǒng)計得到通量參數(shù)。由于其屬于統(tǒng)計學方法，根據(jù)數(shù)學上的大數(shù)定律和中心極限定律可確定，當樣本數(shù)量足夠大時，能夠收斂于真解。將蒙特卡羅方法用于求解中子輸運方程可得到更為精確的通量值，用于動力學可得到更為精確的動力學計算結(jié)果。
針對加速器驅(qū)動次臨界問題，在實際將蒙特卡羅用于動力學各時間步下的中子學求解時，由于需要在計算初始給出共軛通量作為歸一化預估通量的權(quán)重，需要對蒙特卡羅計算進行處理，另外由于問題處于低次臨界度下，模擬時中子有接近于1的增殖比，難以截斷，造成模擬時間長且計算不穩(wěn)定，容易發(fā)散。
從以上分析可看出，要針對加速器驅(qū)動次臨界系統(tǒng)進行精確的瞬態(tài)分析計算，需要利用耦合蒙特卡羅方法進行分析，而耦合蒙特卡羅方法需要解決共軛通量計算和低次臨界度下模擬不穩(wěn)定的問題。
技術實現(xiàn)要素：
：
為了克服上述現(xiàn)有技術存在的問題，本發(fā)明的目的在于提供一種利用耦合蒙特卡羅方法模擬ads系統(tǒng)瞬態(tài)問題的方法，采用轉(zhuǎn)置矩陣和構(gòu)造迭代加速的操作，解決共軛通量計算和低次臨界度下模擬不穩(wěn)定的問題；與全確定論的瞬態(tài)計算方法相比，具有計算精度高的優(yōu)點，同時解決了在耦合蒙卡方法中的問題；該方法將動力學部分與中子輸運部分進行外耦合，方便使用各種蒙卡程序進行替換，具有通用性。
為達到上述目的，本發(fā)明采取了以下技術方案予以實施：
一種利用耦合蒙特卡羅方法模擬ads系統(tǒng)瞬態(tài)問題的方法，步驟如下：
步驟1：采用預估校正改進準靜態(tài)方法處理帶緩發(fā)中子先驅(qū)核的中子時空動力學方程，得到完全的預估通量方程，如公式(1)：
式中：
n--時間步序號；
ω--空間角度；
r--空間位置；
e--中子能量；
σ't,n+1--n+1時刻的形式總截面，包括n+1時刻的總截面和時間相關總截面；
σs,n+1(r,e',ω'→e,ω)--n+1時刻，位置r處，中子從e'能量ω'角度散射到e能量ω角度的散射截面；
f't,n+1--n+1時刻的形式裂變產(chǎn)生截面，包括n+1時刻的瞬發(fā)中子產(chǎn)生截面和n+1時刻的時間相關的緩發(fā)中子產(chǎn)生截面；
sn--n時刻的形式外源，包括實際外源，緩發(fā)中子先驅(qū)核相關源和時間相關源；
同時在預估校正改進準靜態(tài)方法中，需要對預估通量進行歸一化，其權(quán)重因子為共軛通量；
歸一化公式為：
式中：
v--中子速度；
ψ--形狀函數(shù)分布；
步驟2：根據(jù)步驟1公式推導要求，按照公式(4)計算0時刻的共軛通量
式中：
ψ*--共軛通量；
σ--總截面；
σs(r,e,ω→e',ω')--在位置r處，中子從e能量ω角度散射到e'能量ω'角度的散射截面；
νσf(r,e)--在位置r處，能量為e的中子的裂變產(chǎn)生截面；
χ(e')--裂變產(chǎn)生中子能量為e'的占比；
ψ*(r,e',ω')--在位置r處，能量為e'，角度為ω'的共軛通量；
通過轉(zhuǎn)置裂變矩陣和散射矩陣的方式，保持原邊界條件不變，進行前向計算，所得標通量即為對應的共軛通量，其分布能夠用于對動力學計算中預估通量的計算；
步驟3：根據(jù)步驟1推導得出的公式(1)，對公式(1)進行蒙特卡羅輸運計算，其中外源由n時刻的輸運計算得到，n＝0時即為初始時刻實際外源；截面以蒙特卡羅指定材料的形式給出；目標為計算n+1時刻的預估通量
式中：
m--迭代次數(shù)序號；
km--迭代第m步時的裂變源系數(shù)；
km-1--迭代第m-1步時的裂變源系數(shù)；
s0--歸一化后的外源；
φm--迭代第m步時通量
φm-1--迭代第m-1步時的通量
l--泄漏、吸收、散射項算子；
b--裂變產(chǎn)生項算子；
<>--關于全相空間的積分算子；
在蒙特卡羅中子模擬中，將裂變反應視為吸收，在第m-1次計算時模擬計算相應通量，并通過式(6)計算迭代系數(shù)km-1，然后通過式(5)構(gòu)造新的源項，不改變模擬問題的材料組成，進行第m次計算，反復迭代，直到兩次迭代中的km-1和km的差小于10-6時，視為收斂，此時的φm即為n-1時刻的預估通量
步驟4：根據(jù)步驟3求得的n+1時刻的預估通量
式中：
ψn+1--第n+1時刻的形狀函數(shù)；
由預估校正改進準靜態(tài)方法，在形狀函數(shù)已知的情況下，能夠求出n+1時刻的點堆參數(shù)，再由點堆動力學方程能夠得到n+1時刻的通量幅度tn+1；最終可得到n+1時刻的通量：
φn+1＝tn+1ψn+1公式(8)
重復步驟3步驟4繼續(xù)計算n+2時刻的通量，依次類推，直到計算到目標時刻。進而得到通量隨時間和空間的分布φ(r,e,t)，完成反應堆瞬態(tài)計算。
本發(fā)明與現(xiàn)有技術相比具有如下優(yōu)點：
1、采用耦合蒙特卡羅方法，且為外部耦合，可以提高瞬態(tài)計算中各時間步上中子輸運計算的準確度，另外還能靈活選取蒙特卡羅求解程序。
2、將裂變源迭代加速技巧用于蒙特卡羅求解，解決了低次臨界度下固定源問題模擬緩慢且不穩(wěn)定的問題，使蒙特卡羅方法更好地適用于加速器驅(qū)動次臨界系統(tǒng)的瞬態(tài)分析計算。
附圖說明
圖1為利用耦合蒙特卡羅方法模擬ads系統(tǒng)瞬態(tài)問題的方法的計算流程圖。
具體實施方式
下面結(jié)合附圖和具體實施例對本發(fā)明作進一步詳細說明：
本發(fā)明在原有的確定論的預估校正準靜態(tài)計算方法中，將瞬態(tài)分析中分時間步的中子輸運計算用蒙特卡羅方法進行耦合，通過轉(zhuǎn)置裂變截面和散射截面處理共軛通量計算，通過引入裂變源迭代因子，加速并穩(wěn)定蒙特卡羅中子模擬，達到預期效果。
如圖1所示，本發(fā)明的具體實施步驟如下：
步驟1：采用預估校正改進準靜態(tài)方法處理帶緩發(fā)中子先驅(qū)核的中子時空動力學方程，得到完全的預估通量方程：
式中：
n--時間步序號；
ω--空間角度；
r--空間位置；
e--中子能量；
σ't,n+1--n+1時刻的形式總截面，包括n+1時刻的總截面和時間相關總截面；
σs,n+1(r,e',ω'→e,ω)--n+1時刻，位置r處，中子從e'能量ω'角度散射到e能量ω角度的散射截面；
f't,n+1--n+1時刻的形式裂變產(chǎn)生截面，包括n+1時刻的瞬發(fā)中子產(chǎn)生截面和時間相關的緩發(fā)中子產(chǎn)生截面；
sn--n時刻的形式外源，包括實際外源，緩發(fā)中子先驅(qū)核相關源和時間相關源；
同時在預估校正改進準靜態(tài)方法中，需要對預估通量進行歸一化，其權(quán)重因子為共軛通量；
歸一化公式為：
式中：
v--中子速度；
ψ--形狀函數(shù)分布；
步驟2：根據(jù)第1步公式推導要求，計算共軛通量
中子輸運穩(wěn)態(tài)前向方程表達式如下：
中子輸運穩(wěn)態(tài)共軛通量表達式如下：
式中：
ψ*--共軛通量；
σ--總截面；
σs(r,e,ω→e',ω')--在位置r處，中子從e能量ω角度散射到e'能量ω'角度的散射截面；
νσf(r,e)--在位置r處，能量為e的中子的裂變產(chǎn)生截面；
χ(e')--裂變產(chǎn)生中子能量為e'的占比；
ψ*(r,e',ω')--在位置r處，能量為e'，角度為ω'的共軛通量；
由公式(3)和公式(4)可知，在0時刻，通過轉(zhuǎn)置散射截面，和改變裂變中子產(chǎn)生截面與裂變中子能譜的積分順序，再通過前向計算求解方法求得所得標通量，即為對應的共軛通量
步驟3：
根據(jù)步驟1推導得出的公式，對式公式(1)進行蒙特卡羅輸運計算，其中外源sn由n時刻的輸運計算的通量φn得到，n＝0時即為初始時刻實際外源；截面以蒙特卡羅指定材料的形式給出；通過目標為計算n+1時刻的預估通量
式中：
m--迭代次數(shù)序號；
km--迭代第m步時的裂變源系數(shù)；
km-1--迭代第m-1步時的裂變源系數(shù)；
s0--歸一化后的外源；
φm--迭代第m步時通量
φm-1--迭代第m-1步時的通量
l--泄漏、吸收、散射項算子；
b--裂變產(chǎn)生項算子；
<>--關于全相空間的積分算子；
在蒙特卡羅中子模擬中，將裂變反應視為吸收，在第m-1次計算時模擬計算相應通量φm-1，并通過式(6)計算迭代系數(shù)km-1，然后通過式(5)構(gòu)造新的源項，不改變模擬問題的材料組成，將源項改為以上計算的源進行第m次計算，反復迭代，直到兩次迭代中的km-1和km的相對誤差小于10-6時，視為收斂，此時的即為n+1時刻的預估通量
步驟4：
根據(jù)步驟3求得的n+1時刻的預估通量
ψn+1--第n+1時刻的形狀函數(shù)；
由預估校正改進準靜態(tài)方法，在形狀函數(shù)已知的情況下，可以求出n+1時刻的點堆參數(shù)，由此可以得到n+1時刻的通量幅度tn+1；最終可得到n+1時刻的通量：
φn+1＝tn+1ψn+1公式(8)
重復步驟3可繼續(xù)計算n+2時刻的通量，依次類推，直到計算到目標時刻；進而可得到通量隨時間和空間的分布φ(r,e,t)，完成反應堆瞬態(tài)計算。