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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 基于點(diǎn)堆動(dòng)力學(xué)的物理熱工耦合

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

點(diǎn)堆動(dòng)力學(xué)的方程為：

和

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

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

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

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

Λ——中子代時(shí)間；

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

Q——源項(xiàng)；

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

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

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

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

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

針對(duì)不同的反應(yīng)堆，變量Tn都有對(duì)應(yīng)的計(jì)算公式，但是，計(jì)算公式中都包含功率?；邳c(diǎn)堆動(dòng)力學(xué)的物理熱工反饋計(jì)算模型如圖5所示。

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

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

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

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

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

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

式中，n——時(shí)間步；

tn+1——第n+1步的時(shí)間；

tn——第n 步的時(shí)間；

Δtn+1——時(shí)間步長，Δtn+1=tn+1-tn；

t——時(shí)間。

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

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

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

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

Φ——中子通量密度；

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

E——中子能量；

Ω——描述中子運(yùn)動(dòng)方向的向量。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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