核反應(yīng)堆物理熱工耦合的瞬態(tài)分析方法研究

王昆鵬，許 超，李聰新，劉宇生，溫麗晶，王宏凱，劉 健

（生態(tài)環(huán)境部核與輻射安全中心，北京 100082）

核反應(yīng)堆瞬態(tài)分析是核反應(yīng)堆安全及事故分析的基礎(chǔ)，一座核反應(yīng)堆能否安全運行的關(guān)鍵是能否成功控制反應(yīng)堆功率在各種情況下的安全運行[1]。世界上3 次嚴(yán)重的核事故發(fā)生后，人們對核反應(yīng)堆的安全性提出了更高的要求，對核反應(yīng)堆安全分析的要求也更加嚴(yán)格[2]。而核事故的發(fā)生往往伴隨劇烈的功率變化，所以，從安全的角度出發(fā)，核反應(yīng)堆的設(shè)計和事故分析往往都希望得到最接近堆芯實際功率的功率值，從而準(zhǔn)確地模擬反應(yīng)堆功率的空間分布及瞬態(tài)變化[3]。另一方面，堆芯功率的變化又會影響堆芯燃料及冷卻劑等材料的溫度與密度參量，這些參量變化帶來的反饋作用又反過來影響中子通量密度，因此，考慮堆芯物理和熱工的耦合計算是反應(yīng)堆瞬態(tài)分析的重點[4]。

1 物理熱工耦合的基本方式

核反應(yīng)堆的瞬態(tài)分析需要進行堆芯物理和一回路系統(tǒng)的熱工耦合計算。物理熱工耦合的方法有很多種[5-7]，根據(jù)堆芯和一回路系統(tǒng)耦合點的不同，主要分為3類。

第一種耦合方式為內(nèi)部耦合，如圖1 所示。采用該耦合方式進行熱工水力計算時，計算包括堆芯在內(nèi)的所有一回路設(shè)備，而中子學(xué)計算只提供堆芯部分的功率分布。此種耦合方式的特點是中子程序和熱工程序互為邊界條件，計算結(jié)果為對方的初始條件，因此，該耦合方法便于利用已有中子學(xué)程序和系統(tǒng)程序開展物理熱工耦合計算。

第二種耦合方式為外部耦合，如圖2 所示。在該耦合方式中，中子程序計算包括控制棒和堆芯流體在內(nèi)的中子計算；系統(tǒng)程序僅包括除了壓力容器之外的熱工水力計算。因此，進行物理熱工耦合計算時，需要對中子程序進行改進，增加熱工反饋模塊。

圖1 內(nèi)部物理熱工耦合Fig.1 Internal physical thermal coupling

圖2 外部物理熱工耦合Fig.2 External physical thermal coupling

第三種耦合方式為混合耦合，即第一種方式和第二種方式的結(jié)合。在中子計算的同時進行了堆芯部分的熱工水力計算，而系統(tǒng)程序也計算了包括堆芯在內(nèi)的整個一回路的熱工水力，顯然利用此方法進行計算的結(jié)果更加精確，但計算時間較長。

物理熱工耦合的計算模型又有很多種，主要取決于中子動力學(xué)和熱工水力的計算方法。中子動力學(xué)的計算方法主要基于點堆動力學(xué)、中子擴散理論動力學(xué)和中子輸運理論動力學(xué)，熱工水力的計算方法包括點模型、單通道、子通道、系統(tǒng)程序等，國外主要核電發(fā)展國家關(guān)于物理熱工耦合的研究見表1，日本、墨西哥等國基于點堆模型開展了關(guān)于快堆、行波堆等堆型的研究工作[8]；德國、法國和美國都基于中子擴散動力學(xué)和子通道程序耦合的方式開展了氣冷快堆、VVER 反應(yīng)堆的瞬態(tài)研究[9]；美國和德國基于中子輸運動力學(xué)和計算流體動力學(xué)（CFD）程序開展了相關(guān)前瞻性的研究工作；我國也進行了類似的研究工作。

表1 國外物理熱工耦合的常見方法Table 1 Common methods of physical thermal coupling abroad

物理熱工耦合的計算方式包括，一步耦合（顯式耦合，如圖3 所示）和多步耦合（隱式耦合，如圖4 所示）。圖3 中，假定TN時刻的堆芯功率分布和熱工-水力狀態(tài)已經(jīng)求出，采用一次通過的過程得到TN+1時刻的結(jié)果。其優(yōu)點是整個過程只進行了一次中子學(xué)計算和一次熱工水力計算，沒有進行迭代，節(jié)省時間；缺點是時間步長不能太長，否則計算結(jié)果失真度很大。隱式耦合，在中子學(xué)計算和熱工-水力計算之間進行迭代計算，直到中子學(xué)計算和熱工水力計算收斂為止。這種耦合的方式計算結(jié)果較為準(zhǔn)確，然而計算時間較長。

圖3 一步耦合Fig.3 One-step coupling

圖4 多步耦合Fig.4 Multi-step coupling

2 基于點堆動力學(xué)的物理熱工耦合

在早期，研究人員主要采用簡單的“點堆”動力學(xué)近似模型[10，11]，對反應(yīng)堆瞬態(tài)動力學(xué)進行研究。該方法的特點是將整個反應(yīng)堆堆芯功率假定為一個平均值，該值隨時間進行變化。

點堆動力學(xué)的方程為：

和

式中，N(t)——中子密度；

Ci(t)——第i組緩發(fā)中子先驅(qū)核濃度；

λi——第i 組緩發(fā)中子先驅(qū)核的衰變常數(shù)；

βi——第i組緩發(fā)中子份額；

Λ——中子代時間；

β——緩發(fā)中子份額；

Q——源項；

ρ(t)——t時刻的反應(yīng)性。

ρ(t)的取值取決于反應(yīng)堆狀態(tài)，即燃料棒溫度、冷卻劑溫度、冷卻劑密度、空泡份額等，反應(yīng)性根據(jù)上述條件的變化而變化，反應(yīng)性的變化會影響功率的變化。因此，采用點堆動力學(xué)模型進行物理熱工耦合計算時，應(yīng)在每個時間步上對反應(yīng)性進行修正，其計算公式為：

式中， ρ0——初始時刻引入的反應(yīng)性；

αn——變量Tn的反應(yīng)性反饋系數(shù)；

ΔTn(t)——t時刻變量Tn的變化量。

針對不同的反應(yīng)堆，變量Tn都有對應(yīng)的計算公式，但是，計算公式中都包含功率。基于點堆動力學(xué)的物理熱工反饋計算模型如圖5所示。

基于點堆動力學(xué)的物理熱工耦合方法計算效率較高，一般能實現(xiàn)對反應(yīng)堆瞬態(tài)的快速分析，點堆動力學(xué)方程剛性較強,常采用全隱式龍格庫塔方法求解[12]。

3 基于時空中子動力學(xué)的物理熱工耦合

圖5 基于點堆動力學(xué)的物理熱工耦合計算流程Fig.5 Neutronic thermal-hydraulic coupling calculation flow based on point reactor model

點堆模型采用集總參數(shù)的方法，不能描述與空間有關(guān)的動力學(xué)效應(yīng)，僅能對小型緊耦合系統(tǒng)且偏離臨界狀態(tài)不遠的問題給出較滿意的計算結(jié)果；而像大型商用動力堆這類可能發(fā)生中子通量空間分布顯著畸變的情況，“點堆”模型已遠遠不夠。

近年來，研究人員已開始以多群、多維、節(jié)塊擴散模型作為首選進行反應(yīng)堆的安全事故分析和仿真模擬[13]。節(jié)塊擴散模型可以提供堆芯功率空間分布，提高了安全分析的準(zhǔn)確度，使對大型商用堆的分析達到了精度要求。但是，隨著對一些特殊問題研究的不斷深入，如具有體積小、高泄漏和含有強烈吸收體等特點的空間堆，體積小、冷卻劑中又含有大量空泡的溶液堆，以及其他一些由非對稱布置的控制棒移動誘導(dǎo)的反應(yīng)堆瞬態(tài)問題，對反應(yīng)堆安全分析提出了新的課題要求，擴散模型已不能給出精確的結(jié)果，從而需要有更高精度的模擬計算。所以，隨著反應(yīng)堆技術(shù)的發(fā)展，公眾對反應(yīng)堆的安全性和經(jīng)濟性有了更高的追求，使研究人員不再滿足基于擴散理論和利用低維的方法解決反應(yīng)堆動力學(xué)問題，而越來越趨向基于輸運理論的三維時空動力學(xué)方程的求解[14]，以期獲得對堆芯物理瞬態(tài)過程更精確的變化分析結(jié)果。并且核安全要求的是完整的堆芯瞬態(tài)分析，反應(yīng)堆的動力學(xué)受制于中子物理過程和堆芯多種狀況的耦合，事故安全反饋效應(yīng)是必須考慮的。

描述核系統(tǒng)動力學(xué)特性的方程組屬于剛性初值問題，求解時，要求對時間變量的離散問題具有良好的穩(wěn)定性和收斂性。對時空動力學(xué)方程組中的時間變量的處理有兩類方法，即間接近似格式和直接的差分格式。而對時間變量的離散處理通常采用無條件穩(wěn)定的全隱式向后差分格式。用全隱式向后差分格式處理以下常微分初值問題：

式中，n——時間步；

tn+1——第n+1步的時間；

tn——第n 步的時間；

Δtn+1——時間步長，Δtn+1=tn+1-tn；

t——時間。

關(guān)于時空動力學(xué)的研究，已經(jīng)有了很多研究成果，主要集中在對時間變量的處理，具體包括一些數(shù)值處理方法和相應(yīng)的瞬態(tài)計算程序。而各類時間變量離散方法中，因子分解法是將時空相關(guān)的中子通量密度直接分解為一個隨時間變化很緩慢的形狀函數(shù)和一個僅隨時間變化的幅函數(shù)，可以根據(jù)各自對時間的依賴性采取不同的時間步長，以大大減少計算時間，滿足三維時空中子動力學(xué)計算對速度的高要求。同時，因子分解在形式上是嚴(yán)格的，并不存在近似，所以，具有和直接離散方法相同的精度。

如式（6）對時空相關(guān)的中子通量直接進行了因式分解，代入三維時空動力學(xué)方程中可以得到有關(guān)形狀函數(shù)的方程：

式中，ψ( r,E,Ω,t )——描述中子通量密度的形狀函數(shù)；

T( t )——描述中子通量密度的幅函數(shù)；

Φ——中子通量密度；

r——中子在空間中的位置；

E——中子能量；

Ω——描述中子運動方向的向量。

絕熱近似方法完全忽略了式（6）中形狀函數(shù)和幅函數(shù)的時間導(dǎo)數(shù)項，并忽略了先驅(qū)核形狀分布的時間延遲，不區(qū)分緩發(fā)中子源和瞬發(fā)中子源項。

準(zhǔn)靜態(tài)方法[15，16]認(rèn)為形狀函數(shù)隨時間變化很緩慢，完全忽略式（6）中形狀函數(shù)的時間導(dǎo)數(shù)項，其他保持不變，時間變量t 在形狀函數(shù)方程中只呈現(xiàn)為一個參數(shù)，即：

改進準(zhǔn)靜態(tài)方法（Improved Quasistatic Method，IQS）則保留式（8）形狀函數(shù)的時間導(dǎo)數(shù)項，采用全隱式向后差分格式進行時間離散，即：

比較3 種方法可知，絕熱法作了太多近似，精度有待提高，該法已不再使用；準(zhǔn)靜態(tài)方法可以較好地描述漸近的瞬態(tài)變化，然而對于較快的瞬態(tài)變化，特別是階躍變化引起的瞬跳變化則不能給出精確解；而改進準(zhǔn)靜態(tài)方法考慮了形狀函數(shù)隨時間的變化，不論是瞬跳還是漸進的瞬態(tài)過程結(jié)果，都較準(zhǔn)確。

IQS具體到瞬態(tài)時間的離散，形狀函數(shù)和幅函數(shù)對時間變量不同程度的依賴性提供了利用不同時間間隔的機會。形狀函數(shù)只有在需要時才進行更新計算，計算時間之間的間隔選取比幅函數(shù)的時間間隔大得多，在不降低解的精度下，使計算時間大幅度減小?？梢园凑詹煌臅r間間隔劃分分別求解形狀函數(shù)和幅函數(shù)，從而大大減少計算時間，提高了計算效率，解決了直接用離散方法計算速度過慢的問題。同時，規(guī)劃時間步長和各時間步長等級之間的變化成為改進準(zhǔn)靜態(tài)方法的重要內(nèi)容之一。

基于時空中子動力學(xué)的物理熱工耦合，反應(yīng)性截面是物理熱工反饋的關(guān)鍵。因瞬態(tài)計算過程分析的時間足夠短，通常以秒為單位，在這樣的時間尺度下可以認(rèn)為燃料的燃耗是恒定不變的，反應(yīng)堆中的毒物也不隨時間發(fā)生變化，各動力學(xué)參數(shù)變化非常緩慢，這樣各計算區(qū)域的截面僅取決于燃料及冷卻劑的溫度和密度。另外，對于控制棒區(qū)，還要考慮控制棒是否插入。因此，完整的截面群常數(shù)計算公式可以表示為：

式中，f1(Bu)——有關(guān)燃耗的多項式；

f2(v)——有關(guān)溫度的多項式；

f3(p)——有關(guān)相對功率密度的多項式；

f4(Cr)——控制棒的影響變化量。

各參數(shù)的計算公式和具體堆型有關(guān)，由具體堆型確定。

基于三維時空中子動力學(xué)模型的物理熱工耦合計算流程如圖6所示，而反應(yīng)堆的狀態(tài)需要熱工程序進行計算，從而獲得反應(yīng)堆燃料的溫度分布、冷卻劑的溫度和密度分布參數(shù)，可以利用式（9）計算出反應(yīng)性截面，從而進行中子動力學(xué)的計算。

圖6 基于時空中子動力學(xué)的物理熱工耦合計算流程Fig.6 Neutronic thermal-hydraulic coupling calculation flow based on space-time kinetic model

4 物理熱工耦合的發(fā)展趨勢

鑒于計算能力的迅速提高以及新的反應(yīng)堆廣泛采用MOX 燃料，設(shè)計燃耗深度也在不斷提高，燃料的各向異性非常顯著。因此，采用更精確的模型，即中子學(xué)計算采用蒙特卡羅方法求解中子輸運方程[6,17]；熱工水力計算采用子通道模型，甚至采用計算流體動力學(xué)（CFD）方法進行物理熱工耦合計算，都是未來的研究趨勢。另一方面，提高原有輸運程序的計算速度、改進熱工模型、開發(fā)人機接口界面、開發(fā)通用堆型的物理熱工耦合程序也是未來的研究熱點方向。

5 結(jié)論

本文總結(jié)了反應(yīng)堆瞬態(tài)分析中物理熱工耦合的主要方法，分析了各種耦合方法的優(yōu)缺點，結(jié)論如下：

（1）反應(yīng)堆堆芯和一回路系統(tǒng)的耦合方式可以分為內(nèi)部耦合、外部耦合和混合耦合3種方式，其中內(nèi)部耦合方式便于利用已有程序；外部耦合方式便于系統(tǒng)程序建立計算模型；混合耦合方式精度最高，但計算效率較低；

（2）物理熱工耦合的方式根據(jù)同一個時間步長內(nèi)是否迭代可以分為一步法和多步法，其中，一步法計算速度較快，但計算步長不能太大；多步法計算精度較高，但效率較低；

（3）根據(jù)物理熱工模型維數(shù)的不同，可以劃分為基于點堆動力學(xué)的耦合方法和基于時空中子動力學(xué)的耦合方法。早期受制于計算條件，一般采用基于點堆動力學(xué)的耦合方法。而時空中子動力學(xué)耦合模型，尤其是基于蒙特卡羅方法的時空中子動力學(xué)模型因廣泛的堆芯適應(yīng)性及高計算精度，并且可獲得各時間步長下堆芯三維功率、溫度等參數(shù)的分布，是未來物理熱工耦合模型的發(fā)展方向。