明平洲,李治剛,潘俊杰,安 萍,蘆 韡,劉 東,余紅星,孫玉發(fā)

(1.中國(guó)核動(dòng)力研究設(shè)計(jì)院核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室,四川 成都 610213;2.中國(guó)核動(dòng)力研究設(shè)計(jì)院,四川 成都 610213)

事故分析以及瞬態(tài)工況分析是反應(yīng)堆堆芯計(jì)算的重要組成部分。隨著計(jì)算機(jī)軟硬件的演變,多學(xué)科代碼的耦合能夠更真實(shí)地模擬反應(yīng)堆行為。中子物理和熱工水力是廣為研究的耦合計(jì)算方案,在堆芯尺度內(nèi),它能夠?qū)Χ研緩?fù)雜的瞬態(tài)工況進(jìn)行細(xì)致描述,為設(shè)計(jì)或系統(tǒng)安全分析提供多種參數(shù)。中子時(shí)空動(dòng)力學(xué)在擴(kuò)散計(jì)算近似條件下可以提供堆芯內(nèi)部瞬態(tài)過(guò)程中重要中子學(xué)參數(shù)隨時(shí)間的變化,通過(guò)耦合熱工反饋計(jì)算功能,提升計(jì)算準(zhǔn)確度,且瞬態(tài)工況下在數(shù)值求解時(shí)包含時(shí)間變量的離散和空間變量的離散。論文基于核動(dòng)力院現(xiàn)有的研究條件和計(jì)算方案[1],提取和集成自研軟件CORTH[2]內(nèi)的瞬態(tài)熱工反饋計(jì)算功能,為完整的堆芯耦合計(jì)算方案提供原型參考。除開(kāi)對(duì)瞬態(tài)熱工反饋計(jì)算功能的解讀,福清機(jī)組的插棒提棒問(wèn)題也將用于實(shí)際堆芯問(wèn)題的數(shù)值驗(yàn)證。

1 計(jì)算方案簡(jiǎn)介

中子時(shí)空動(dòng)力學(xué)模型如式 (1)所示。式中,下標(biāo)g代表能群,能群數(shù)為G;下標(biāo)i代表緩發(fā)先驅(qū)核組,組數(shù)為ND。?g(→r,t)為g群中子通量密度;ci(→r,t)為第i組先驅(qū)核濃度;vg為g群中子速度;Dg為g群擴(kuò)散系數(shù);Σr,g和Σf,g分別為g群移出截面和裂變截面,Σs,g′→g為g′群到g群的散射截面;ν為每次裂變釋放的中子數(shù);χg為瞬發(fā)中子裂變譜;χg,i為緩發(fā)中子譜份額;λi為先驅(qū)核衰變常數(shù);βi為緩發(fā)中子份額。

整個(gè)中子時(shí)空動(dòng)力學(xué)計(jì)算方案現(xiàn)階段以瞬態(tài)中子學(xué)計(jì)算為主干,獲得堆芯的三維功率分布信息。熱工反饋計(jì)算根據(jù)相應(yīng)信息計(jì)算得到堆芯內(nèi)部的溫場(chǎng)分布,從而引起堆芯各個(gè)位置處的截面數(shù)據(jù)變化,再影響瞬態(tài)中子學(xué)計(jì)算的結(jié)果。由于這些非線性關(guān)系,總的數(shù)值方法采用源迭代方法來(lái)消除非線性關(guān)系,獲得各個(gè)瞬態(tài)時(shí)刻的堆芯內(nèi)參數(shù)分布。相應(yīng)計(jì)算方案如圖1所示。

圖1 計(jì)算方案的流程圖Fig.1 Flow chart of the calculation scheme

時(shí)間離散方面,龍格庫(kù)塔方法是精度較高的常微分方程數(shù)值解法。相關(guān)研究在隱式歐拉格式基礎(chǔ)上采用精度從1階到4階的對(duì)角線隱式龍格-庫(kù)塔方法 (DIRK),并利用Richardson外推和嵌入低階方法實(shí)現(xiàn)了時(shí)間步長(zhǎng)的自適應(yīng),該方法具有良好的穩(wěn)定性和較高精度[3]。根據(jù)式 (1),相應(yīng)的數(shù)學(xué)模型可以簡(jiǎn)化為式 (2)。式 (3)則是利用式 (2)推導(dǎo)得到的一階龍格庫(kù)塔數(shù)值形式,即向后歐拉格式。實(shí)際計(jì)算方案可以采用更精確的二階DIRK方法進(jìn)行時(shí)間離散。每一級(jí)的求解與向后歐拉格式時(shí)間離散后的求解形式是相似的,區(qū)別只在于替換每一級(jí)的時(shí)步和初值。

中子學(xué)計(jì)算方面主要分為擴(kuò)散計(jì)算和截面更新計(jì)算兩部分。其中擴(kuò)散計(jì)算在時(shí)空動(dòng)力學(xué)中主要指代瞬態(tài)擴(kuò)散計(jì)算,它可以通過(guò)固定源問(wèn)題下的格林函數(shù)粗網(wǎng)節(jié)塊法進(jìn)行求解[3]。截面計(jì)算則采用讀表插值的方法,利用組件程序提供的截面庫(kù),根據(jù)各個(gè)瞬態(tài)時(shí)刻的狀態(tài)參數(shù)來(lái)更新堆芯內(nèi)部的截面數(shù)據(jù)?，F(xiàn)階段計(jì)算方案中根據(jù)計(jì)算問(wèn)題的性質(zhì)和側(cè)重點(diǎn),選定冷卻劑密度和燃料有效溫度兩個(gè)狀態(tài)參數(shù)。瞬態(tài)計(jì)算在進(jìn)入之前通過(guò)穩(wěn)態(tài)計(jì)算或讀取狀態(tài)數(shù)據(jù)庫(kù)來(lái)提供初始時(shí)刻的通量、先驅(qū)核的分布,對(duì)應(yīng)的熱工反饋計(jì)算從圖1可以明確,位于每個(gè)瞬態(tài)時(shí)刻的擴(kuò)散計(jì)算與截面更新計(jì)算之間,該內(nèi)容將在論文中進(jìn)行詳細(xì)解讀和分析。

2 熱工反饋計(jì)算的分析

2.1 概述

熱工反饋計(jì)算主要根據(jù)三維全堆芯的功率分布或通量分布,計(jì)算得到堆芯內(nèi)部的溫場(chǎng)分布。與中子學(xué)部分耦合之后,相應(yīng)的熱工參數(shù)將作為狀態(tài)點(diǎn),用于截面計(jì)算,反過(guò)來(lái)影響堆芯的功率分布。通過(guò)多次數(shù)值迭代計(jì)算達(dá)到平衡。這里熱工反饋計(jì)算從CORTH程序中進(jìn)行提取,這里將對(duì)相應(yīng)的計(jì)算進(jìn)行解讀和分析。

2.2 熱工反饋計(jì)算

考慮到工程計(jì)算的易用性,選擇單通道模型來(lái)進(jìn)行熱工反饋計(jì)算,此時(shí)熱工反饋計(jì)算部分的軸向分段應(yīng)與中子學(xué)部分保持一致,但僅僅考慮活性區(qū);徑向上也與中子學(xué)的粗網(wǎng)格保持相同大小,例如中子學(xué)瞬態(tài)擴(kuò)散計(jì)算按照組件進(jìn)行粗網(wǎng)格的劃分。由于按照組件進(jìn)行了通道劃分,因此整個(gè)熱工反饋計(jì)算由兩層循環(huán)進(jìn)行數(shù)值計(jì)算,如圖2所示。

圖2 熱工反饋計(jì)算的偽碼描述Fig.2 Pseudo code description for the thermal feedback calculation

這里遍歷各個(gè)組件進(jìn)行熱工參數(shù)計(jì)算,從計(jì)算機(jī)編程角度進(jìn)行解釋:首先根據(jù)中子學(xué)代碼提供的功率信息對(duì)每個(gè)組件的各個(gè)軸向分段求解表面平均熱流密度dens,然后調(diào)用Cal FuelSig函數(shù)求解各個(gè)網(wǎng)格處 (由組件與軸向分段進(jìn)行劃分)冷卻劑的密度和燃料有效溫度。由于區(qū)分了熱工反饋計(jì)算過(guò)程內(nèi)的數(shù)值變量和中子學(xué)單學(xué)科代碼所依賴的數(shù)值變量,所以最后在反饋計(jì)算完成之后需要進(jìn)行耦合數(shù)據(jù)的賦值操作。

2.3 冷卻劑密度和燃料有效溫度的計(jì)算

瞬態(tài)單通道模型不考慮橫向流動(dòng),因此按照燃料組件遍歷時(shí)各個(gè)通道內(nèi)的質(zhì)量流速和壓降比較明確,不需要求解動(dòng)量守恒方程便可完成熱工反饋所需參數(shù)的計(jì)算。瞬態(tài)單通道模型從燃料元件中心開(kāi)始將通道在徑向上分為6個(gè)控制體,所以某個(gè)通道的各個(gè)分段冷卻劑溫度信息需要保存6個(gè),與商業(yè)軟件PRINA類似。其中,燃料芯塊劃分為4個(gè)控制體,包殼和外部慢化劑分別作為1個(gè)控制體。

圖3中燃料棒徑向分為6個(gè)控制體,其中1,2,3為燃料芯塊,4為燃料芯塊和氣隙 (如果存在),5為包殼,6為冷卻劑。按照各個(gè)控制體生成三對(duì)角矩陣,求解得到冷卻劑溫度等信息。此時(shí)瞬態(tài)單通道模型中采用的假設(shè)為:不考慮軸向?qū)?不考慮流體勢(shì)能,動(dòng)能為常數(shù)且均勻,且流體流速方向在同一高度共線。因此瞬態(tài)單通道模型中只需要求解能量守恒方程便可完成反饋計(jì)算。

圖3 燃料元件的控制體劃分Fig.3 Mesh partition of fuel element

芯塊控制體:

氣隙控制體:

包殼控制體:

慢化劑控制體:

以上能量守恒方程中S為內(nèi)熱源,離散化之后可得到相應(yīng)的三對(duì)角形式方程:

相關(guān)系數(shù)A、B、C、α均使用已知量計(jì)算,涉及到控制體相關(guān)的幾何參數(shù)、熱物性和換熱系數(shù)等內(nèi)容。對(duì)于棒狀燃料元件,瞬態(tài)當(dāng)前時(shí)刻的冷卻劑溫度為前一時(shí)刻與當(dāng)前時(shí)刻冷卻劑溫度的平均值。冷卻劑密度可以通過(guò)物性函數(shù)的轉(zhuǎn)換,將冷卻劑溫度轉(zhuǎn)化為密度信息。

瞬態(tài)情況下需要考慮上一時(shí)刻的燃料溫度信息來(lái)計(jì)算導(dǎo)熱系數(shù)以及當(dāng)前時(shí)刻燃料有效溫度的初始狀態(tài)。計(jì)算燃料芯塊溫度常常包含有芯塊外表面溫度、中心溫度以及有效溫度等內(nèi)容,同時(shí)求解燃料芯塊溫度時(shí)需要區(qū)分板狀燃料元件 (軍用)和棒狀燃料元件,然后按照軸向分段對(duì)各個(gè)分段進(jìn)行遍歷求解。計(jì)算方案中集中在棒狀燃料元件的求解,在徑向上將燃料元件進(jìn)行分區(qū),如圖2所示分為三個(gè)分區(qū),整個(gè)芯塊溫度的計(jì)算便圍繞這些分區(qū)展開(kāi),徑向上利用一維熱傳導(dǎo)方程進(jìn)行求解。

最終燃料芯塊的中心溫度計(jì)算得到后,利用加權(quán)方法,聯(lián)合燃料芯塊的外表面溫度獲得燃料芯塊的有效溫度值,計(jì)算公式為:

3 數(shù)值實(shí)驗(yàn)

使用福清5&6號(hào)機(jī)組第1循環(huán)壽期初的瞬態(tài)例題對(duì)整個(gè)計(jì)算方案進(jìn)行數(shù)值驗(yàn)證。該例題為熱態(tài)零功率的初始工況 (零功率條件下熱工反饋較小,但可以觀察數(shù)值穩(wěn)定性和敏感性等內(nèi)容),其中N1棒組以1.2步每秒的速度插入,183.3 s時(shí)再以1.2步每秒的速度提出,366.6 s離開(kāi)真實(shí)堆芯的活性區(qū),統(tǒng)計(jì)堆芯500 s的瞬變過(guò)程。實(shí)驗(yàn)中真實(shí)堆芯的軸向分為20段,并利用商業(yè)SMART軟件進(jìn)行對(duì)比驗(yàn)算。

固定步長(zhǎng)運(yùn)行方式主要指代龍格庫(kù)塔時(shí)間離散的計(jì)算步長(zhǎng)取為固定,實(shí)驗(yàn)中取為0.1 s;而自適應(yīng)步長(zhǎng)運(yùn)行方式表明采用嵌入低階格式以估計(jì)當(dāng)前時(shí)步的計(jì)算誤差,在計(jì)算過(guò)程中可能放大或縮小計(jì)算步長(zhǎng),減少計(jì)算次數(shù)。整個(gè)計(jì)算過(guò)程中輸出模擬瞬態(tài)時(shí)間500 s內(nèi)的相對(duì)功率變化趨勢(shì),反映該物理問(wèn)題的瞬變過(guò)程。

分析一:由于相對(duì)功率值偏小,縱坐標(biāo)的數(shù)值取對(duì)數(shù),便于對(duì)比數(shù)值的變化趨勢(shì)。從圖4中可以看到,在瞬態(tài)時(shí)間200多秒時(shí)SMART計(jì)算結(jié)果出現(xiàn)了相對(duì)功率先降低后升高的波谷情況,子通道模型下固定步長(zhǎng)或自適應(yīng)步長(zhǎng)的計(jì)算結(jié)果與SMART參考解的趨勢(shì)一致,相對(duì)功率均經(jīng)歷了先降后升的過(guò)程。利用數(shù)據(jù)處理工具比較單通道模型下固定步長(zhǎng)與自適應(yīng)步長(zhǎng)的相對(duì)功率值差異,兩者的相對(duì)偏差最大為0.98%,說(shuō)明兩者吻合較好;比較單通道模型下固定步長(zhǎng)與SMART計(jì)算結(jié)果的差異,相對(duì)功率值的相對(duì)偏差最大為17.5%,最小為10.3%,即計(jì)算結(jié)果存在一定偏差。值得注意的是相對(duì)功率的波谷處SMART計(jì)算結(jié)果與論文內(nèi)子通道模型的計(jì)算結(jié)果偏差約為1%左右,說(shuō)明瞬態(tài)模擬的關(guān)鍵事件位置不存在較大誤差,這里考慮計(jì)算結(jié)果在其他瞬態(tài)時(shí)刻的偏差可能來(lái)源于商業(yè)軟件SMART內(nèi)置的部分參數(shù)與數(shù)值實(shí)驗(yàn)中使用的控制參數(shù)不一致導(dǎo)致。

圖4 瞬態(tài)過(guò)程的相對(duì)功率變化圖Fig.4 Relative power trend of transient process

分析二:統(tǒng)計(jì)在工作站服務(wù)器HP Z800a上的計(jì)算時(shí)間,固定步長(zhǎng)的計(jì)算時(shí)間為5 353.78 s,自適應(yīng)步長(zhǎng)完成計(jì)算則耗時(shí)為314.38 s。實(shí)驗(yàn)中記錄得到固定步長(zhǎng)情況下熱工反饋計(jì)算次數(shù)為50 000,自適應(yīng)步長(zhǎng)情況下熱工反饋計(jì)算次數(shù)為236。這表明自適應(yīng)步長(zhǎng)在本例題中能夠較多減少離散計(jì)算的時(shí)間點(diǎn)個(gè)數(shù),所以串行情況下本瞬態(tài)例題自適應(yīng)步長(zhǎng)的計(jì)算效率優(yōu)于固定步長(zhǎng)的計(jì)算效率,且正確性上兩者幾乎保持一致。

分析三:瞬態(tài)時(shí)間250 s左右相對(duì)功率值達(dá)到波谷,但實(shí)際上控制棒組在180 s左右插入到堆芯底部再提出。統(tǒng)計(jì)瞬態(tài)時(shí)間范圍內(nèi)熱工參數(shù)的數(shù)值大小,整個(gè)瞬態(tài)時(shí)間區(qū)間內(nèi)冷卻劑密度的變化和燃料有效溫度的數(shù)值變化在本例題中并不明顯,但變化趨勢(shì)與插棒提棒動(dòng)作比較吻合。這里以軸向第十層的組件為例,其坐標(biāo)為 (6,6)和 (8,8),(8,8)位置對(duì)應(yīng)堆芯徑向正中心組件。繪制對(duì)應(yīng)節(jié)塊在計(jì)算過(guò)程中冷卻劑密度和燃料有效溫度的變化趨勢(shì)圖,圖5表明冷卻劑密度變化比燃料有效溫度更為劇烈,180 s左右降到最低,且徑向上控制棒組附近的熱工參數(shù)變化幅度比堆芯正中心組件更大;從參數(shù)敏感性上來(lái)看,燃料有效溫度受控制棒組在軸向上的移動(dòng)較小,其中一段瞬態(tài)時(shí)間的數(shù)值未明顯變化,符合本例題代表的物理問(wèn)題特點(diǎn)。

綜上所述,本文闡述的計(jì)算方案具備在中子時(shí)空動(dòng)力學(xué)計(jì)算過(guò)程中提供瞬態(tài)熱工反饋計(jì)算功能,且通過(guò)真實(shí)堆芯例題可以初步歸納,現(xiàn)階段在計(jì)算結(jié)果上符合該例題對(duì)應(yīng)真實(shí)瞬態(tài)過(guò)程的變化趨勢(shì)。盡管熱工參數(shù)的影響遵循堆芯內(nèi)部的物理過(guò)程,但計(jì)算精度上與SMART商業(yè)軟件存在著少量偏差,這些偏差可能來(lái)源于多個(gè)因素,并不影響原型方案的有效性。

圖5 熱工反饋參數(shù)的變化圖Fig.5 The variation of thermal feedback parameters

4 結(jié)論

本文在中子時(shí)空動(dòng)力學(xué)的計(jì)算過(guò)程中引入了瞬態(tài)熱工反饋計(jì)算功能,通過(guò)對(duì)熱工反饋計(jì)算內(nèi)容的解讀和分析,說(shuō)明了整個(gè)計(jì)算方案在堆芯問(wèn)題域內(nèi)的求解流程,為中國(guó)核動(dòng)力研究設(shè)計(jì)院相關(guān)專業(yè)軟件的研制提供原型參考。與此同時(shí),真實(shí)堆芯的數(shù)值實(shí)驗(yàn)初步論證了整個(gè)原型計(jì)算方案的正確性和有效性。隨著研究的進(jìn)一步深入和拓展,子通道模型以及CFD方法將對(duì)熱工反饋計(jì)算進(jìn)行替換和對(duì)比驗(yàn)證,同時(shí)先進(jìn)的并行計(jì)算等內(nèi)容也將引入到計(jì)算方案中進(jìn)行計(jì)算效率的增強(qiáng)。