乏燃料溶液系統(tǒng)臨界安全分析計(jì)算方法研究

銀華北 王永平 茍軍利 劉國明 祖鐵軍 尹 文 鄭友琦 杜夏楠

1（西安交通大學(xué) 西安 710049）

2（中國核電工程有限公司 北京 100840）

為實(shí)現(xiàn)我國閉式循環(huán)目標(biāo)，自主發(fā)展后處理技術(shù)、修建后處理廠具有重要意義。由于乏燃料后處理過程中核燃料會發(fā)生固-液相變、富集等復(fù)雜過程，涉及眾多設(shè)備和環(huán)節(jié)，存在發(fā)生系統(tǒng)達(dá)臨界事故的風(fēng)險(xiǎn)。因此，在乏燃料后處理廠中，對核臨界事故的預(yù)防和防護(hù)，是保證核安全的關(guān)鍵內(nèi)容。

在乏燃料溶液系統(tǒng)的臨界事故中，臨界瞬間功率快速上升，伴隨劇烈的輻解、相變傳質(zhì)過程，針對固體燃料的核反應(yīng)堆堆芯瞬態(tài)分析程序無法模擬上述過程。傳統(tǒng)的核臨界安全事故分析方法采用經(jīng)驗(yàn)公式估計(jì)裂變次數(shù)，精度及適用范圍有限。

針對溶液系統(tǒng)的上述特征，國內(nèi)外機(jī)構(gòu)開展了瞬態(tài)分析研究。英國原子能委員會首先發(fā)布了CRITEX程序，將溶液罐視作圓管，沿軸向分層進(jìn)行一維熱工水力-輻解氣體計(jì)算，與點(diǎn)堆動力學(xué)耦合進(jìn)行瞬態(tài)分析。法國與日本分別利用其溶液系統(tǒng)臨界實(shí)驗(yàn)裝置SILENE［1］與TRACY［2］進(jìn)行實(shí)驗(yàn)研究溶液系統(tǒng)臨界事故下的功率釋放與溶液變化，并開發(fā)了AGNES2［3］、TRACE［4］、FECTH［5］等 程 序。其 中，TRACE與AGNES2程序由日本原子能研究所先后開發(fā)，其中子學(xué)均采用點(diǎn)堆動力學(xué)計(jì)算，而熱工-水力采用R-Z二維軸對稱模型，將溶液對流近似等效為導(dǎo)熱，二者區(qū)別在輻解氣體模型，TRACE程序基于CRITEX模型，將氣泡速度近似為僅功率變化有關(guān)的函數(shù)，而AGNES2程序由氣泡在溶液中的受力平衡導(dǎo)出氣泡速度關(guān)系式；AGNES2程序在實(shí)驗(yàn)驗(yàn)證中精度優(yōu)于TRACE程序。FETCH程序由法國開發(fā)，采用球諧函數(shù)展開下的有限元方法進(jìn)行中子學(xué)計(jì)算，并直接使用兩相流計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）軟件進(jìn)行熱工水力-輻解氣體建模計(jì)算，且考慮了溶液流動對緩發(fā)中子先驅(qū)核分布的影響，其驗(yàn)證結(jié)果與實(shí)驗(yàn)值符合較好，然而由于采用了CFD模型，導(dǎo)致計(jì)算效率低，工程實(shí)用性差。

中國核電工程公司開發(fā)了CACCS程序［6］，動力學(xué)計(jì)算采用點(diǎn)堆模型，熱工水力與輻解氣體模型與AGNES2程序相似。清華大學(xué)開發(fā)了TCCHAR程序［7］，其中子學(xué)計(jì)算采用蒙特卡羅方法，熱工模型采用集總模型，輻解氣體模型在考慮產(chǎn)生、流入與流出的基礎(chǔ)上還考慮了氣泡大小在軸向的變化，其驗(yàn)證結(jié)果與實(shí)驗(yàn)值誤差在15%以內(nèi)。然而，在瞬態(tài)計(jì)算中，基于蒙特卡羅的中子學(xué)計(jì)算顯然存在計(jì)算效率低、瞬態(tài)中的小反應(yīng)性變化易被統(tǒng)計(jì)方差淹沒的缺點(diǎn)。

可見，國內(nèi)外現(xiàn)有程序中大部分采用點(diǎn)堆動力學(xué)進(jìn)行中子學(xué)計(jì)算，使用余弦函數(shù)與貝塞爾函數(shù)近似系統(tǒng)內(nèi)的功率分布與反應(yīng)性反饋權(quán)重，其優(yōu)勢在于計(jì)算效率高，然而由于大量乏燃料后處理的需求，后續(xù)后處理工業(yè)示范廠規(guī)模大幅增加，勢必加大溶液處理設(shè)備容量。而大容量溶液處理設(shè)備為保證系統(tǒng)安全大量布置控制毒物，并控制容器形狀保證幾何安全，上述假設(shè)不再適用，勢必引入誤差，對乏燃料存儲策略、臨界事故預(yù)防措施的制定帶來不確定性；此外，隨著計(jì)算機(jī)科學(xué)的發(fā)展，三維計(jì)算模型的計(jì)算時間逐漸下降至可接受范圍內(nèi)，基于時空動力學(xué)的溶液系統(tǒng)瞬態(tài)分析程序變得可能，從而進(jìn)一步提高安全分析結(jié)論的可靠性。因此，開發(fā)一套幾何適用范圍廣、模型精度高的臨界安全分析程序，對于提高核臨界事故的分析精度、保障乏燃料后處理廠的安全具有重要的意義。

為此，本研究開發(fā)了一套用于溶液系統(tǒng)的確定論并行三維瞬態(tài)分析程序Hydra-TD，主要包含截面制作模塊，中子動力學(xué)模塊及熱工-輻解氣體模塊，其中中子動力學(xué)模塊采用離散縱標(biāo)并行輸運(yùn)程序計(jì)算，在保持三維中子學(xué)計(jì)算幾何適用性廣、分辨率高的優(yōu)點(diǎn)的同時保證中子學(xué)計(jì)算效率；熱工水力-輻解氣體模塊在AGNES2模型基礎(chǔ)上進(jìn)行簡化以加速計(jì)算；熱工水力-輻解氣體-中子學(xué)耦合中考慮了溶液的膨脹、輻解氣體對熱散射的影響等因素，增加耦合精度。將程序與SILENE實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比驗(yàn)證的結(jié)果誤差較小，顯示程序具備臨界安全分析的能力。

本文對研究工作分為三個部分介紹，首先為程序的理論模型，介紹中子學(xué)截面制作、時空動力學(xué)模型和熱工-輻解氣體模型，以及中子學(xué)-熱工水力-輻解氣體模型的耦合方法；隨后利用法國SILENE實(shí)驗(yàn)裝置對本研究開發(fā)程序的模型進(jìn)行了檢驗(yàn)，其裂變次數(shù)、裂變率等結(jié)果精度較高；最后所得結(jié)論為本研究開發(fā)程序能夠正確模擬瞬發(fā)臨界第一裂變峰的功率變化，但無法還原功率下降階段的振蕩，輻解氣體模型可能需要進(jìn)一步改進(jìn)。

1 理論模型

1.1 截面制作與產(chǎn)生

本文采用多群近似進(jìn)行中子學(xué)計(jì)算，需要獲取溶液系統(tǒng)各材料各個狀態(tài)的有效多群截面。多群截面制作模塊流程如下：

首先，使用核數(shù)據(jù)處理程序NECP-Atlas［8］，從評價核數(shù)據(jù)庫出發(fā)，進(jìn)行共振重構(gòu)及線性化處理、多普勒展寬、不可分辨區(qū)有效自屏截面處理，熱中子散射計(jì)算和多群常數(shù)計(jì)算，制作問題無關(guān)的多群截面庫；然后，在非共振段使用溫度和背景截面插值獲得多群參數(shù)，在共振能量段使用超細(xì)群方法求解均勻問題的慢化方程，使用通量歸并獲得問題相關(guān)的多群參數(shù)。

式中：σi，k，g是第i區(qū)、k核素、多群g的有效自屏截面；σi，k，h是第i區(qū)、核素k、超細(xì)群h的截面；φi，h是第i區(qū)、超細(xì)群h的中子通量密度。

共振計(jì)算中，針對超細(xì)群方法的特點(diǎn)，模塊采用了散射源加速及碰撞概率表插值以提高計(jì)算效率。散射源加速方法下，為避免計(jì)算源項(xiàng)時進(jìn)行大量超細(xì)群的累加，每次計(jì)算時從上次的散射源中減去高能群超過最高對數(shù)能降的散射源，加上低能群的散射源：

式中：Sjg為j區(qū)域內(nèi)其他能群散射至g能群的散射源；Δuf為對數(shù)能群寬度；Σsjk，g為j區(qū)域內(nèi)k核素在g能群宏觀散射截面；PNk，k為中子穿越Nk個能群后與核素k發(fā)生核反應(yīng)的概率；φjg是j區(qū)域g群通量。

為提高計(jì)算效率，多群截面制作模塊產(chǎn)生WIMS-69群各材料微觀截面后，使用NECP-Hydra程序［9］進(jìn)行一次溶液系統(tǒng)的穩(wěn)態(tài)計(jì)算，使用全問題能譜將各材料截面歸并為適用于熱譜堆芯的4群2階散射截面。在STACY燃料溶液罐問題中的穩(wěn)態(tài)計(jì)算結(jié)果表明這一歸并能夠保持對系統(tǒng)有效增殖因子keff及溶液內(nèi)通量分布的計(jì)算精度，keff偏差為十萬分之十四，如表1所示；將溶液區(qū)軸向分為三層，徑向分為4層，分群統(tǒng)計(jì)其通量分布（按蒙特卡羅程序輸出的一個源中子的統(tǒng)計(jì)結(jié)果進(jìn)行歸一），所得分群分布偏差最大為5%，如表2所示。

表1 STACY問題keff計(jì)算結(jié)果對比Table 1 Comparison on calculated keff of STACY problem

表2 STACY問題溶液區(qū)分群通量計(jì)算結(jié)果對比Table 2 Comparison on calculated neutron flux distribution of STACY problem

針對溶液系統(tǒng)的輻射分解現(xiàn)象，為考慮輻解氣體對熱散射的影響，上述計(jì)算時同時制作了游離氫的截面數(shù)據(jù)，在瞬態(tài)計(jì)算時將對應(yīng)的輻解氣體由水中氫轉(zhuǎn)化為游離氫。

1.2 中子時空動力學(xué)計(jì)算

本文在中子時空動力學(xué)計(jì)算中采用預(yù)估-校正準(zhǔn)靜態(tài)方法［10］，將中子通量密度與緩發(fā)中子先驅(qū)核密度分解為幅值與形狀函數(shù)分別求解，在大時間步上求解空間形狀分布，從而減少輸運(yùn)求解次數(shù)節(jié)省計(jì)算時間。

在多群、預(yù)估校正準(zhǔn)靜態(tài)的離散下，幅函數(shù)滿足精確點(diǎn)堆方程，其具體形式見文獻(xiàn)［10］，求解時，中子密度方程采用全隱式向后差分離散，緩發(fā)中子先驅(qū)核方程對其右端產(chǎn)生項(xiàng)βˉd(t)n(t)/Λ(t)采用亞當(dāng)斯方法離散，即作二階插值多項(xiàng)式在［tn，tn+1］內(nèi)積分，其余項(xiàng)采用解析積分；中子通量密度形狀函數(shù)的求解采用全隱式向后差分，需要求解的方程形式與有源次臨界問題類似。

在輸運(yùn)方程的求解上，為保證計(jì)算效率，本文采用圓柱幾何r-θ-z六面體網(wǎng)格下的KBA（Koch-Baker-Alcouffe）并 行離 散縱 標(biāo)（SN）方法，基 于NECP-Hydra程序［9］完成開發(fā)。在KBA并行算法下，輸運(yùn)問題被空間區(qū)域分解，根據(jù)各個SN角度上的空間依賴關(guān)系，依次在各個并行單元上逐層逐角度啟動掃描，如圖1所示。KBA并行算法能避免迭代格式退化，是結(jié)構(gòu)幾何下并行效率最高的算法之一。在500×500×500網(wǎng)格的輸運(yùn)問題中，NECP-Hydra在百核下能保持90%的并行效率。

圖1 KBA并行算法示意圖Fig.1 Diagram of KBA parallel sweep algorithm

除NECP-Hydra中已經(jīng)采用的切比雪夫外推加速、擴(kuò)散綜合加速方法［11］外，針對上述時間相關(guān)固定源問題，本文還采用了如下加速方法［12］：

1）對固定源項(xiàng)歸一：

式中：S0為外中子源的歸一化形狀分布；·為在問題相空間區(qū)域內(nèi)積分值。

2）進(jìn)行類似于冪法的迭代過程直至迭代收斂。

式中：M為包含泄漏、消失與散射貢獻(xiàn)的輸運(yùn)算子；F為裂變源算子。在瞬態(tài)計(jì)算過程中，每次輸運(yùn)計(jì)算的通量矩分布、k特征值都被保留，作為下次計(jì)算的迭代初值，在通量形狀變化不大時可大幅減少迭代次數(shù)。

對于精確點(diǎn)堆方程系數(shù)及形狀函數(shù)方程時間相關(guān)項(xiàng)中的角度相關(guān)項(xiàng)，本文采用與離散縱標(biāo)方法內(nèi)部儲存一致的PN通量矩展開表示，即：

為使形式簡潔，方程（6）省去了空間變量，方程（7）、（8）省去了時間和空間變量。其中：φk，lg為g能群k，l階通量矩；Ylk(Ω)為對應(yīng)的球諧函數(shù)；ψ*g為求解系統(tǒng)穩(wěn)態(tài)共軛輸運(yùn)方程得到的共軛角通量；ψ?g為歸一化后的中子角通量形狀函數(shù)。

由于瞬態(tài)過程中系統(tǒng)各處的溶液組分隨時變化，為便于考慮溶液的各種物理變化，本研究使用微觀截面進(jìn)行中子動力學(xué)計(jì)算，在每個時間步按當(dāng)?shù)販囟炔逯但@得各核素的微觀截面，隨后結(jié)合當(dāng)前區(qū)域的核子密度計(jì)算得到宏觀截面，代入中子輸運(yùn)模塊求解。

1.3 熱工-輻解氣體計(jì)算

在溶液系統(tǒng)臨界過程中，由于功率升高，溶液溫度上升，體積、濃度也發(fā)生變化。同時，裂變反應(yīng)產(chǎn)生的高能碎片與燃料溶液中的分子發(fā)生碰撞，產(chǎn)生輻照裂解氣體（H2、O2、N2等），并形成微泡；此外，中子和伽馬射線一并與分子的相互作用對輻解氣體的產(chǎn)生亦有貢獻(xiàn)［1］。上述因素在瞬態(tài)過程中影響熱散射、多普勒效應(yīng)及核子密度分布，其計(jì)算精度關(guān)系到中子輸運(yùn)計(jì)算精度。

對于溶液中輻解氣體計(jì)算，經(jīng)調(diào)研國內(nèi)外相關(guān)研究，本文參考AGNES2程序［3］的模型，建立如下氣泡輸運(yùn)方程：

式中：Fi，j為（i，j）控制體內(nèi)的空泡體積份額；νt為能量轉(zhuǎn)換系數(shù)，m6·J-1·mol-1，Pi，j為控制體內(nèi)功率密度，W·m-3；Ci，j為控制體內(nèi)輻解氣體濃度，mol·L-1；C0為輻解氣體臨界釋放濃度；θ為單位階躍函數(shù)；vi，j為氣泡軸向遷移速度，m·s-1；G為輻解氣體產(chǎn)生率，mol·J-1；τ為氣體溶解時間。

上述方程中的未知量，除空泡體積份額、輻解氣體量外，還包括氣泡遷移速度。然而，溶液中輻解氣體氣泡的遷移過程受溶液對流、相間曳力、浮升力、氣泡的生長、合并、溶解、破裂等多種因素影響，精確模擬其運(yùn)動需進(jìn)行精細(xì)建模的兩相流數(shù)值計(jì)算，將極大影響計(jì)算效率。為此，本文中的輻解氣體遷移速度采用近似關(guān)系式，通過一系列預(yù)先選取的兩相流工況數(shù)值計(jì)算結(jié)果對當(dāng)?shù)毓β拭芏?、?dāng)?shù)厝芤簞恿W(xué)黏度與當(dāng)?shù)貕毫Φ纫蜃兞繑M合得到。

對于系統(tǒng)熱工計(jì)算，本文在溶液區(qū)及容器區(qū)分別計(jì)算導(dǎo)熱，并在區(qū)域交界面上將對流、輻射換熱等效為導(dǎo)熱系數(shù)進(jìn)行計(jì)算：

式中：ρ為材料密度；Cp為定壓比熱；T為溫度；λeff為導(dǎo)熱系數(shù)或等效導(dǎo)熱系數(shù)；Φ?為體積熱源即功率密度。

1.4 耦合方法

在耦合計(jì)算中，中子動力學(xué)模塊求解獲得各個時刻的功率分布，熱工-輻解氣體模塊求解獲得各個時刻的溫度、質(zhì)量密度及輻解氣體濃度分布，兩者通過耦合接口按體積或質(zhì)量權(quán)重轉(zhuǎn)化為所需網(wǎng)格上的分布；質(zhì)量密度與輻解氣體濃度分布通過如下計(jì)算轉(zhuǎn)化為中子學(xué)計(jì)算所需的核子密度。

式中：Ninuc為核素密度，l表示溶液初始包含核素的集合；εl為網(wǎng)格內(nèi)溶液體積份額；G為輻解氣體產(chǎn)生率，mol·J-1；Φ?為功率密度；NA為阿伏伽德羅常數(shù)；δ為Kronecker符號，僅當(dāng)下標(biāo)相同時為1，其余時刻為0；l to gas表示溶液中輻解為氣體的核素集合；nigas∈gas為輻解氣體平均每分子內(nèi)各核素的原子數(shù)；cgas為輻解氣體的摩爾濃度分布。上式的物理含義為溶液在升溫膨脹的同時被輻解使一部分核素析出并轉(zhuǎn)移，其核素密度按物質(zhì)的量混合至當(dāng)?shù)鼐W(wǎng)格內(nèi)，如圖2所示。

此外，隨溶液膨脹與輻解氣體析出，溶液液位將發(fā)生變化。耦合計(jì)算中，熱工-輻解氣體模塊根據(jù)溶液受熱膨脹后密度與輻解氣體所占體積得出總體積，從而得出溶液液位，耦合接口中結(jié)合液位在上溢網(wǎng)格內(nèi)獲得當(dāng)?shù)厝芤号c空氣的體積份額，隨后利用當(dāng)?shù)厝芤嘿|(zhì)量密度獲得當(dāng)?shù)睾怂孛芏?，以考慮該效應(yīng)：

式中：Voverflow為發(fā)生溶液上溢的網(wǎng)格；air表示氣體包含核素的集合；εair為網(wǎng)格內(nèi)空氣體積份額；ρa(bǔ)ir為當(dāng)?shù)乜諝饷芏?。式?4）和（15）的物理含義為：在溶液由于膨脹上溢處將溶液與空氣按各自的體積份額混合至當(dāng)?shù)鼐W(wǎng)格，如圖2所示。

圖2 中子學(xué)與熱工水力-輻解氣體耦合中的體積混合示意圖Fig.2 Volume mixing operation when coupling neutronics with thermal-hydraulics / radiolysis gas

本文在耦合中采用算子分裂方法，即將溫度、輻解氣體場與中子場分離解耦，先后求解后更新至下一時間步。算子分裂方法是對多物理非線性耦合問題的線性化，其誤差隨時間步減小而減小。

本研究采用的計(jì)算流程如下：首先進(jìn)行一次穩(wěn)態(tài)前向計(jì)算，以獲取初始邊界條件；隨后進(jìn)行共軛計(jì)算，以共軛通量作為準(zhǔn)靜態(tài)下分離形狀函數(shù)與幅函數(shù)的權(quán)重函數(shù)；進(jìn)入瞬態(tài)計(jì)算后：

1）先更新截面進(jìn)行預(yù)估求解，獲取中子通量及先驅(qū)核濃度的形狀函數(shù)；

2）隨后在大時間步內(nèi)進(jìn)行形狀函數(shù)插值更新點(diǎn)堆參數(shù)；

3）進(jìn)行點(diǎn)堆求解獲取校正通量及功率分布；

4）從而進(jìn)行熱工-輻解氣體計(jì)算；

5）根據(jù)溫度分布進(jìn)行截面插值，根據(jù)質(zhì)量密度及輻解氣體含量更新核素密度分布從而歸并得到宏觀截面。

瞬態(tài)計(jì)算的時間步進(jìn)流程如圖3所示。

圖3 瞬態(tài)時間步進(jìn)流程Fig.3 Transient time stepping process of Hydra-TD

2 數(shù)值驗(yàn)證

按照§1的理論模型，本文編制了圓柱幾何乏燃料溶液系統(tǒng)臨界安全時空瞬態(tài)分析程序。為了驗(yàn)證該程序的計(jì)算精度，對SILENE實(shí)驗(yàn)裝置的部分實(shí)驗(yàn)進(jìn)行了建模計(jì)算。SILENE實(shí)驗(yàn)裝置是一個硝酸鈾酰溶液的均勻?qū)嶒?yàn)堆（235U富集度達(dá)92.7%），位于法國Burgundy的Valduc中心，主要用于臨界事故研究，但還可以用于核加熱、劑量測量、臨界訓(xùn)練。堆芯在一個很大的混凝土房間中央，是一個小的環(huán)形水槽，如圖4所示。不銹鋼容器高度為1 m，外徑360/368 mm，外壁厚度4 mm，內(nèi)徑70/76 mm，內(nèi)壁厚度3 mm，底部厚度36 mm，頂部厚度30 mm。反應(yīng)性引入棒為環(huán)形鎘棒，厚度為1 mm，內(nèi)外均有1 mm厚的不銹鋼包殼，直徑60/66 mm，長度為1 130 mm［1］。

圖4 SILENE實(shí)驗(yàn)裝置縱截面示意圖Fig.4 Diagram of axial section of SILENE facility

計(jì)算設(shè)備采用Linux 18.04 LTS 64位系統(tǒng)，處理器是Intel Xeon Gold 6230 2.10GHz，內(nèi)存空間為8 G，采用7個CPU核心進(jìn)行計(jì)算。中子學(xué)計(jì)算中，劃分的網(wǎng)格為r-θ-z方向24×1×210，角度采用層對稱S4求積組。熱工輻解氣體計(jì)算中，對容器壁的網(wǎng)格進(jìn)行細(xì)化，對容器上部空氣部分網(wǎng)格進(jìn)行粗化，其網(wǎng)格為r-z方向54×134。設(shè)置臨界濃度為10 mol·m-3，輻解氣體產(chǎn)生率為1.5×10-7mol·J-1，能量轉(zhuǎn)換系數(shù)為1×10-7mol·J-1［3］。

本文所計(jì)算實(shí)驗(yàn)為S2-300實(shí)驗(yàn)，其初始液位高度為38.91 cm，通過以2 cm·s-1的速度提升反應(yīng)性引入棒，共引入約0.97 $的反應(yīng)性。對該問題進(jìn)行時間步長敏感性分析，分別在脈沖階段取大時間步為62.5 ms、31.25 ms、20 ms下計(jì)算，其裂變率曲線對比如圖5所示，時間步為31.25 ms與20 ms的裂變率曲線已重合，可見在脈沖階段取大時間步為31.25 ms能夠使結(jié)果收斂。計(jì)算結(jié)果與實(shí)驗(yàn)測量結(jié)果的對比如圖6所示，溶液高度、溫度及空泡份額的變化和反應(yīng)性變化繪于圖7。隨著反應(yīng)性引入，系統(tǒng)功率開始指數(shù)上升，使溶液溫度上升，密度下降，溶液高度升高；同時輻解氣體開始積累，到達(dá)臨界濃度后析出，使反應(yīng)性快速下降。

圖5 S2-300實(shí)驗(yàn)各時間步長模擬計(jì)算裂變率結(jié)果對比Fig.5 Comparison of calculated results of fission rate for S2-300 experiment with different time step length

圖6 S2-300實(shí)驗(yàn)裂變率與溫度變化計(jì)算結(jié)果（彩圖見網(wǎng)絡(luò)版）(a) 脈沖期間裂變率與溫度變化，(b) 全過程裂變率與溫度變化Fig.6 Calculated results of fission rate and temperature for S2-300 experiment (color online)(a) The change of fission rate and temperature in the pulse duration, (b) The change of fission rate and temperature in the whole process

計(jì)算獲得的裂變率曲線趨勢在第一功率峰處與實(shí)驗(yàn)值曲線一致，平均溫度曲線位于實(shí)驗(yàn)值波動區(qū)間內(nèi)。計(jì)算值與實(shí)驗(yàn)值的詳細(xì)對比列于表2，峰值裂變率誤差為22%，總裂變次數(shù)、峰值前裂變次數(shù)、功率倍增時間與峰值時間的誤差分別為10%、2%、4%和3%，均在10%以內(nèi)，以上結(jié)果表明，程序正確模擬了瞬態(tài)過程中溶液系統(tǒng)的復(fù)雜的多物理過程。

然而，由圖6（b）中與實(shí)驗(yàn)值的比較結(jié)果可知，現(xiàn)階段程序還無法精確模擬第一裂變峰后的功率振蕩，因此未能模擬第二裂變峰，計(jì)算所得溶液溫度在第二裂變峰后偏低，功率下降偏慢。由圖7中空泡反應(yīng)性反饋與空泡份額的變化可知，可能的原因是程序高估了空泡的反饋，或?qū)馀葸w移速度存在低估，使負(fù)反饋偏強(qiáng)，這將在后續(xù)重點(diǎn)研究和改進(jìn)。

圖7 S2-300實(shí)驗(yàn)反應(yīng)性與溶液變化計(jì)算結(jié)果Fig.7 Calculated results of reactivity and variation of solution system for S2-300 experiment

表3 S2-300實(shí)驗(yàn)功率變化特性計(jì)算結(jié)果Table 3 Calculated power feature of S2-300 experiment

3 結(jié)語

本文基于超細(xì)群共振方法、預(yù)估-校正準(zhǔn)靜態(tài)方法及導(dǎo)熱微分方程與輻解氣體輸運(yùn)方程，開發(fā)了一套用于圓柱幾何溶液系統(tǒng)的臨界安全分析的并行三維瞬態(tài)分析程序，并對SILENE實(shí)驗(yàn)裝置的S2-300實(shí)驗(yàn)進(jìn)行了建模計(jì)算，結(jié)果表明，程序能夠正確模擬瞬發(fā)臨界第一裂變峰的功率變化；程序現(xiàn)階段無法還原S2-300的在第一裂變峰后的功率振蕩現(xiàn)象，由溶液空泡份額下降較慢得知可能是輻解氣體相關(guān)模型的影響，這部分內(nèi)容將在后續(xù)工作中重點(diǎn)研究和改進(jìn)。

致謝本文研究內(nèi)容受國防科工局乏燃料后處理科研專項(xiàng)中“后處理廠核臨界安全事故研究”項(xiàng)目資助，特此致謝。

作者貢獻(xiàn)聲明銀華北：程序開發(fā)、驗(yàn)證計(jì)算、結(jié)果分析、起草文章；王永平：程序開發(fā)及理論指導(dǎo)，文章整體設(shè)計(jì)和修改；茍軍利：熱工水力-輻解氣體模塊開發(fā)；劉國明：提供實(shí)驗(yàn)數(shù)據(jù)，結(jié)果分析；祖鐵軍：截面模塊程序開發(fā)及理論指導(dǎo)；尹文：截面制作程序模塊開發(fā)；鄭友琦：三維時空動力學(xué)程序開發(fā)理論指導(dǎo)；杜夏楠：程序截面模塊驗(yàn)證指導(dǎo)。