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

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

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

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

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

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

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

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

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

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

1 理論模型

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

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

首先，使用核數(shù)據(jù)處理程序NECP-Atlas［8］，從評(píng)價(jià)核數(shù)據(jù)庫(kù)出發(fā)，進(jìn)行共振重構(gòu)及線(xiàn)性化處理、多普勒展寬、不可分辨區(qū)有效自屏截面處理，熱中子散射計(jì)算和多群常數(shù)計(jì)算，制作問(wèn)題無(wú)關(guān)的多群截面庫(kù)；然后，在非共振段使用溫度和背景截面插值獲得多群參數(shù)，在共振能量段使用超細(xì)群方法求解均勻問(wèn)題的慢化方程，使用通量歸并獲得問(wèn)題相關(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ì)算中，針對(duì)超細(xì)群方法的特點(diǎn)，模塊采用了散射源加速及碰撞概率表插值以提高計(jì)算效率。散射源加速方法下，為避免計(jì)算源項(xiàng)時(shí)進(jìn)行大量超細(xì)群的累加，每次計(jì)算時(shí)從上次的散射源中減去高能群超過(guò)最高對(duì)數(shù)能降的散射源，加上低能群的散射源：

式中：Sjg為j區(qū)域內(nèi)其他能群散射至g能群的散射源；Δuf為對(duì)數(shù)能群寬度；Σsjk，g為j區(qū)域內(nèi)k核素在g能群宏觀散射截面；PNk，k為中子穿越Nk個(gè)能群后與核素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ì)算，使用全問(wèn)題能譜將各材料截面歸并為適用于熱譜堆芯的4群2階散射截面。在STACY燃料溶液罐問(wèn)題中的穩(wěn)態(tài)計(jì)算結(jié)果表明這一歸并能夠保持對(duì)系統(tǒng)有效增殖因子keff及溶液內(nèi)通量分布的計(jì)算精度，keff偏差為十萬(wàn)分之十四，如表1所示；將溶液區(qū)軸向分為三層，徑向分為4層，分群統(tǒng)計(jì)其通量分布（按蒙特卡羅程序輸出的一個(gè)源中子的統(tǒng)計(jì)結(jié)果進(jìn)行歸一），所得分群分布偏差最大為5%，如表2所示。

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

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

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

1.2 中子時(shí)空動(dòng)力學(xué)計(jì)算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)于溶液中輻解氣體計(jì)算，經(jīng)調(diào)研國(guó)內(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；τ為氣體溶解時(shí)間。

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

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

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

1.4 耦合方法

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

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

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

本研究采用的計(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）隨后在大時(shí)間步內(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ì)算的時(shí)間步進(jìn)流程如圖3所示。

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

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

按照§1的理論模型，本文編制了圓柱幾何乏燃料溶液系統(tǒng)臨界安全時(shí)空瞬態(tài)分析程序。為了驗(yàn)證該程序的計(jì)算精度，對(duì)SILENE實(shí)驗(yàn)裝置的部分實(shí)驗(yàn)進(jìn)行了建模計(jì)算。SILENE實(shí)驗(yàn)裝置是一個(gè)硝酸鈾酰溶液的均勻?qū)嶒?yàn)堆（235U富集度達(dá)92.7%），位于法國(guó)Burgundy的Valduc中心，主要用于臨界事故研究，但還可以用于核加熱、劑量測(cè)量、臨界訓(xùn)練。堆芯在一個(gè)很大的混凝土房間中央，是一個(gè)小的環(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，長(zhǎng)度為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個(gè)CPU核心進(jìn)行計(jì)算。中子學(xué)計(jì)算中，劃分的網(wǎng)格為r-θ-z方向24×1×210，角度采用層對(duì)稱(chēng)S4求積組。熱工輻解氣體計(jì)算中，對(duì)容器壁的網(wǎng)格進(jìn)行細(xì)化，對(duì)容器上部空氣部分網(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，通過(guò)以2 cm·s-1的速度提升反應(yīng)性引入棒，共引入約0.97 $的反應(yīng)性。對(duì)該問(wèn)題進(jìn)行時(shí)間步長(zhǎng)敏感性分析，分別在脈沖階段取大時(shí)間步為62.5 ms、31.25 ms、20 ms下計(jì)算，其裂變率曲線(xiàn)對(duì)比如圖5所示，時(shí)間步為31.25 ms與20 ms的裂變率曲線(xiàn)已重合，可見(jiàn)在脈沖階段取大時(shí)間步為31.25 ms能夠使結(jié)果收斂。計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量結(jié)果的對(duì)比如圖6所示，溶液高度、溫度及空泡份額的變化和反應(yīng)性變化繪于圖7。隨著反應(yīng)性引入，系統(tǒng)功率開(kāi)始指數(shù)上升，使溶液溫度上升，密度下降，溶液高度升高；同時(shí)輻解氣體開(kāi)始積累，到達(dá)臨界濃度后析出，使反應(yīng)性快速下降。

圖5 S2-300實(shí)驗(yàn)各時(shí)間步長(zhǎng)模擬計(jì)算裂變率結(jié)果對(duì)比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é)果（彩圖見(jiàn)網(wǎng)絡(luò)版）(a) 脈沖期間裂變率與溫度變化，(b) 全過(guò)程裂變率與溫度變化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ì)算獲得的裂變率曲線(xiàn)趨勢(shì)在第一功率峰處與實(shí)驗(yàn)值曲線(xiàn)一致，平均溫度曲線(xiàn)位于實(shí)驗(yàn)值波動(dòng)區(qū)間內(nèi)。計(jì)算值與實(shí)驗(yàn)值的詳細(xì)對(duì)比列于表2，峰值裂變率誤差為22%，總裂變次數(shù)、峰值前裂變次數(shù)、功率倍增時(shí)間與峰值時(shí)間的誤差分別為10%、2%、4%和3%，均在10%以?xún)?nèi)，以上結(jié)果表明，程序正確模擬了瞬態(tài)過(guò)程中溶液系統(tǒng)的復(fù)雜的多物理過(guò)程。

然而，由圖6（b）中與實(shí)驗(yàn)值的比較結(jié)果可知，現(xiàn)階段程序還無(wú)法精確模擬第一裂變峰后的功率振蕩，因此未能模擬第二裂變峰，計(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é)語(yǔ)

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

致謝本文研究?jī)?nèi)容受?chē)?guó)防科工局乏燃料后處理科研專(zhuān)項(xiàng)中“后處理廠(chǎng)核臨界安全事故研究”項(xiàng)目資助，特此致謝。

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