中子時空動力學內瞬態熱工反饋計算的分析

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

(1.中國核動力研究設計院核反應堆系統設計技術重點實驗室,四川 成都 610213;2.中國核動力研究設計院,四川 成都 610213)

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

1 計算方案簡介

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

整個中子時空動力學計算方案現階段以瞬態中子學計算為主干,獲得堆芯的三維功率分布信息。熱工反饋計算根據相應信息計算得到堆芯內部的溫場分布,從而引起堆芯各個位置處的截面數據變化,再影響瞬態中子學計算的結果。由于這些非線性關系,總的數值方法采用源迭代方法來消除非線性關系,獲得各個瞬態時刻的堆芯內參數分布。相應計算方案如圖1所示。

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

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

中子學計算方面主要分為擴散計算和截面更新計算兩部分。其中擴散計算在時空動力學中主要指代瞬態擴散計算,它可以通過固定源問題下的格林函數粗網節塊法進行求解[3]。截面計算則采用讀表插值的方法,利用組件程序提供的截面庫,根據各個瞬態時刻的狀態參數來更新堆芯內部的截面數據。現階段計算方案中根據計算問題的性質和側重點,選定冷卻劑密度和燃料有效溫度兩個狀態參數。瞬態計算在進入之前通過穩態計算或讀取狀態數據庫來提供初始時刻的通量、先驅核的分布,對應的熱工反饋計算從圖1可以明確,位于每個瞬態時刻的擴散計算與截面更新計算之間,該內容將在論文中進行詳細解讀和分析。

2 熱工反饋計算的分析

2.1 概述

熱工反饋計算主要根據三維全堆芯的功率分布或通量分布,計算得到堆芯內部的溫場分布。與中子學部分耦合之后,相應的熱工參數將作為狀態點,用于截面計算,反過來影響堆芯的功率分布。通過多次數值迭代計算達到平衡。這里熱工反饋計算從CORTH程序中進行提取,這里將對相應的計算進行解讀和分析。

2.2 熱工反饋計算

考慮到工程計算的易用性,選擇單通道模型來進行熱工反饋計算,此時熱工反饋計算部分的軸向分段應與中子學部分保持一致,但僅僅考慮活性區;徑向上也與中子學的粗網格保持相同大小,例如中子學瞬態擴散計算按照組件進行粗網格的劃分。由于按照組件進行了通道劃分,因此整個熱工反饋計算由兩層循環進行數值計算,如圖2所示。

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

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

2.3 冷卻劑密度和燃料有效溫度的計算

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

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

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

芯塊控制體:

氣隙控制體:

包殼控制體:

慢化劑控制體:

以上能量守恒方程中S為內熱源,離散化之后可得到相應的三對角形式方程:

相關系數A、B、C、α均使用已知量計算,涉及到控制體相關的幾何參數、熱物性和換熱系數等內容。對于棒狀燃料元件,瞬態當前時刻的冷卻劑溫度為前一時刻與當前時刻冷卻劑溫度的平均值。冷卻劑密度可以通過物性函數的轉換,將冷卻劑溫度轉化為密度信息。

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

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

3 數值實驗

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

固定步長運行方式主要指代龍格庫塔時間離散的計算步長取為固定,實驗中取為0.1 s;而自適應步長運行方式表明采用嵌入低階格式以估計當前時步的計算誤差,在計算過程中可能放大或縮小計算步長,減少計算次數。整個計算過程中輸出模擬瞬態時間500 s內的相對功率變化趨勢,反映該物理問題的瞬變過程。

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

圖4 瞬態過程的相對功率變化圖Fig.4 Relative power trend of transient process

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

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

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

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

4 結論

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