巫英偉,賀亞男,章 靜,田文喜,蘇光輝,秋穗正

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院,陜西省先進(jìn)核能技術(shù)重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710049)

核能作為一種清潔、高效的能源形式,對滿足人類不斷增長的能源需求發(fā)揮著重要作用。核能技術(shù)的發(fā)展與安全運(yùn)行密不可分,這需要對核反應(yīng)堆進(jìn)行深入的研究和精確的分析。核反應(yīng)堆系統(tǒng)龐雜且運(yùn)行環(huán)境嚴(yán)苛,存在多物理場耦合的復(fù)雜現(xiàn)象,其涉及的物理過程涵蓋了中子輸運(yùn)、熱傳導(dǎo)、燃料性能演變、冷卻劑流動換熱等諸多方面,這導(dǎo)致反應(yīng)堆安全運(yùn)行面臨諸多挑戰(zhàn)性難題。采用多物理場耦合分析方法解決此類問題可提高核反應(yīng)堆運(yùn)行的安全性和經(jīng)濟(jì)性。

多物理場耦合分析方法通常采用高保真數(shù)值計(jì)算模型及先進(jìn)求解技術(shù),結(jié)合不同物理學(xué)模型,將各種物理過程和空間尺度進(jìn)行有效耦合,并依托超算平臺及大規(guī)模并行技術(shù),構(gòu)建多物理場耦合分析平臺,以全面、系統(tǒng)地模擬核反應(yīng)堆的行為。在過去十幾年里,多物理場耦合分析方法在核工程領(lǐng)域取得了顯著進(jìn)展,并被廣泛用于預(yù)測和模擬復(fù)雜系統(tǒng)的行為。美國開展了CASL[1]、NEAMS[2]和NRC-CRAB計(jì)劃[3],搭建了VERA[4]和CRAB等堆芯多物理場耦合計(jì)算的軟件集成平臺;歐洲以SALOME為基礎(chǔ)建立了堆芯多物理場耦合計(jì)算體系[5];韓國則建立了DeCART-MATRA軟件包以實(shí)現(xiàn)多物理場耦合的計(jì)算[6]。國內(nèi)方面,各科研院所與高校展開合作,開展了多物理場耦合高保真數(shù)值反應(yīng)堆研發(fā)工作。中國原子能科學(xué)研究院、中國科學(xué)院與北京科技大學(xué)建立了數(shù)值反應(yīng)堆原型系統(tǒng)CVR(China Virtual Reactor)[7]。北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所同上海核工程研究設(shè)計(jì)院、西安交通大學(xué)針對國和一號(CAP1400)建立了數(shù)值反應(yīng)堆[8]。中國核動力研究設(shè)計(jì)院同哈爾濱工程大學(xué)同樣也開展了相關(guān)工作[9]。

有限元方法作為廣泛應(yīng)用于各個工程技術(shù)領(lǐng)域的數(shù)值分析方法,在固體導(dǎo)熱、彈塑性力學(xué)、失穩(wěn)變形、熱工水力、裂紋擴(kuò)展等復(fù)雜現(xiàn)象的數(shù)值模擬分析具有一定的優(yōu)勢。因此,不少機(jī)構(gòu)基于有限元方法開發(fā)了燃料性能分析程序,如美國的BISON[10]和FALCON[11]、阿根廷的DIONISIO[12]、法國的ALCYONE[13]以及韓國的PRIME[14]等。除此之外,也有不少學(xué)者在商用有限元程序COMSOL和ABAQUS等的基礎(chǔ)上進(jìn)行二次開發(fā),實(shí)現(xiàn)對燃料性能的模擬[15-17]。有限元分析憑借其靈活的建模、準(zhǔn)確的數(shù)值求解等優(yōu)勢,將會是未來燃料性能分析的研究趨勢,因此其在各多物理場耦合框架中也作為主流的分析工具。

傳統(tǒng)系統(tǒng)分析程序與子通道程序框架固定、缺乏靈活的程序架構(gòu)以及先進(jìn)的數(shù)值算法和物理模型,而先進(jìn)的多物理場耦合平臺可以實(shí)現(xiàn)各物理模型的靈活植入,充分發(fā)揮多物理場耦合的優(yōu)勢,可彌補(bǔ)傳統(tǒng)分析程序在軟件結(jié)構(gòu)及算法上的不足。因此,搭建多物理場耦合平臺,針對耦合問題中的關(guān)鍵技術(shù)開展研究,對加快我國自主化多物理場耦合平臺的開發(fā)進(jìn)程具有重要意義。

西安交通大學(xué)核反應(yīng)堆熱工水力研究室(Nuclear THermal-hydraulic Laboratory, XJTU-NuTHeL)結(jié)合多物理場耦合平臺與有限元方法的優(yōu)勢,開展了高保真模擬工具的開發(fā),包括系統(tǒng)分析程序NUSAC(Nuclear System Analysis Code)和子通道分析程序FLARE(Fully-coupled transient code for Liquid-metal-cooled Advanced REactor)[18-19]、針對多種燃料的燃料性能分析程序BEEs[20-22],并構(gòu)建了基于有限元方法的多物理場耦合模擬體系[23],實(shí)現(xiàn)了堆內(nèi)關(guān)鍵現(xiàn)象的精細(xì)分析。

本文旨在介紹XJTU-NuTHeL為建立多物理場耦合平臺所進(jìn)行的熱工流體計(jì)算模型的開發(fā)、燃料性能分析技術(shù)的研究以及多物理場耦合框架的建立等工作,為核反應(yīng)堆多維度多物理場耦合高保真數(shù)值模擬分析提供有益參考。

1 熱工流體計(jì)算模型開發(fā)

核反應(yīng)堆熱工水力分析是核反應(yīng)堆設(shè)計(jì)、安審、運(yùn)行的重要組成部分,因此熱工流體的精確計(jì)算對于核能的發(fā)展具有重要意義。單相及兩相流動是核反應(yīng)堆正常運(yùn)行和瞬態(tài)工況下非常重要的物理現(xiàn)象,準(zhǔn)確模擬流動現(xiàn)象是核反應(yīng)堆系統(tǒng)分析程序的關(guān)鍵。針對不同的尺度,堆芯熱工水力分析可分為系統(tǒng)安全分析和子通道分析兩種。隨著數(shù)值計(jì)算方法和多物理場耦合平臺的發(fā)展,開發(fā)模塊化、更高精度的熱工水力成為國際熱點(diǎn),如RELAP7[24]、SAM[25]、MOLTRES[26]。相較而言,當(dāng)前國內(nèi)在高階、高精度、高保真的熱工水力分析程序開發(fā)方面研究較少。XJTU-NuTHeL采用高精度間斷有限元方法實(shí)現(xiàn)了單相及兩相物理模型計(jì)算,開發(fā)了核反應(yīng)堆系統(tǒng)安全分析程序NUSAC和子通道分析程序FLARE。

1.1 核反應(yīng)堆系統(tǒng)兩相流分析模型開發(fā)

在核反應(yīng)堆系統(tǒng)流動傳熱模擬中,XJTU-NuTHeL基于先進(jìn)多物理場耦合框架針對壓水堆和先進(jìn)反應(yīng)堆開展了一系列研究,開發(fā)了核反應(yīng)堆系統(tǒng)安全分析程序NUSAC。在單相流模型中,NUSAC采用高階連續(xù)伽遼金有限元離散和SUPG(Streamline-Upwind-Petrov-Galerkin)與PSPG(Pressure Stabilizing-Petrov Galerkin)穩(wěn)定性算法處理了不可壓縮流體流動數(shù)值不穩(wěn)定性問題[27]。其中單相流守恒方程為:

(1)

式中:ρ、t、u、p、T和z分別為密度、時間、流速、壓力、溫度和z軸位置;ψ為形函數(shù);圓括號(F,ψ)代表函數(shù)F在求解域內(nèi)的積分離散;g、f和De分別為重力加速度、阻力系數(shù)和水力直徑;τPSPG和τSUPG分別為與擾動相關(guān)的穩(wěn)定系數(shù);H、cp和q?分別為流體的焓、比定壓熱容和體積釋熱率。

在兩相流模型中,NUSAC采用了一維兩流體六方程模型和經(jīng)過廣泛驗(yàn)證的與RELAP5相同的本構(gòu)模型,其中兩相流守恒方程為:

(2)

式中:下標(biāo)l、g、i、w分別代表液相、汽相、汽-液交界面、壁面;α和Г分別為空泡份額和單位體積相變流量;F、e、Q分別為阻力、比內(nèi)能和熱源;h*和h′分別為流體與壁面?zhèn)髻|(zhì)的相焓和流體與相間傳質(zhì)的相焓。

同時,在兩相流數(shù)值方法中,NUSAC引入間斷重構(gòu)伽遼金數(shù)值離散方法和Roe-type對流項(xiàng)數(shù)值算法實(shí)現(xiàn)了高階空間離散格式,顯著減小了兩相流的數(shù)值擴(kuò)散。在單相和兩相流數(shù)值求解過程中,NUSAC通過自動微分方法形成雅可比矩陣,縮短了程序開發(fā)周期,避免了人為錯誤構(gòu)建雅可比矩陣引入的誤差[18]。此外,NUSAC采用了PJFNK(Preconditioned Jacobian Free Newton-Krylov)全耦合求解方法,在保證精度的同時提高了JFNK方法的收斂性[19]。

XJTU-NuTHeL通過一系列國際通用的兩相流測試基準(zhǔn)題對NUSAC進(jìn)行了廣泛而全面的驗(yàn)證。通過水龍頭問題驗(yàn)證了兩相流方程的正確性,如圖1所示。在此基礎(chǔ)上,NUSAC模擬了Bartolomei實(shí)驗(yàn)的若干組工況[28],并與實(shí)驗(yàn)值和RELAP5計(jì)算值進(jìn)行了對比,如圖2所示,驗(yàn)證了NUSAC處理兩相過冷沸騰工況的能力。NUSAC通過模擬FRIGG實(shí)驗(yàn),驗(yàn)證了NUSAC處理飽和沸騰工況的能力,如圖3所示。過冷沸騰工況和飽和沸騰工況含有大量本構(gòu)關(guān)系式,如汽-液界面阻力、汽-液界面?zhèn)鳠岬?通過這兩個基準(zhǔn)實(shí)驗(yàn)也驗(yàn)證了NUSAC本構(gòu)模型植入的正確性。

圖1 水龍頭問題驗(yàn)證[19]

圖3 FRIGG飽和沸騰實(shí)驗(yàn)驗(yàn)證[19]

NUSAC也可用于反應(yīng)堆系統(tǒng)穩(wěn)態(tài)和瞬態(tài)計(jì)算,根據(jù)美國核能署主蒸汽管道斷裂事故中的堆芯設(shè)計(jì)參數(shù)[29],模擬了穩(wěn)態(tài)及降功率瞬態(tài)的響應(yīng),計(jì)算結(jié)果如圖4所示。圖4中:T_l為一回路流體溫度,K;T_second為蒸汽發(fā)生器二次側(cè)流體溫度,K。可以看出,系統(tǒng)回路穩(wěn)態(tài)溫度分布符合預(yù)期,驗(yàn)證了NUSAC系統(tǒng)模擬的能力。同時,在穩(wěn)態(tài)工況下引入降功率瞬態(tài)變化,在堆芯通道內(nèi)冷卻劑進(jìn)出口溫度和包殼最高溫度分布響應(yīng)如圖5所示,NUSAC計(jì)算值與RELAP5計(jì)算值符合良好,驗(yàn)證了NUSAC的正確性。

a——計(jì)算對象示意圖;b——冷卻劑溫度分布云圖

a——冷卻劑進(jìn)出口溫度;b——包殼最高溫度

1.2 液態(tài)金屬快堆子通道分析模型開發(fā)

液態(tài)金屬快堆由于具有良好的固有安全性以及核燃料增殖能力,被認(rèn)為是最有發(fā)展前景的第4代反應(yīng)堆[30]。為了保證快堆的安全運(yùn)行,有必要對其堆芯進(jìn)行熱工水力安全分析。子通道分析方法是核反應(yīng)堆堆芯熱工水力的重要分析手段,具有比系統(tǒng)分析方法更高的計(jì)算精度,而且相較于CFD分析方法,其計(jì)算資源消耗更小。因此,XJTU-NuTHeL建立了適用于液態(tài)金屬快堆的子通道分析模型,開發(fā)了全耦合子通道瞬態(tài)分析程序FLARE。

FLARE的控制方程如下。

質(zhì)量方程:

(3)

軸向動量方程:

(4)

橫向動量方程:

(5)

能量方程:

(6)

通過上述四方程模型,FLARE在計(jì)算中考慮了通道間的橫向交混、湍流交混以及橫向?qū)?。除此之?FLARE具有鈉、鉛以及鉛鉍的物性狀態(tài)函數(shù),可以對上述3種冷卻劑進(jìn)行計(jì)算。同時,由于液體金屬快堆具有特殊的繞絲結(jié)構(gòu),FLARE本構(gòu)模型中還考慮了繞絲對通道幾何、摩擦阻力以及湍流交混的影響,如圖6所示。程序基于有限元與間斷有限元方法相結(jié)合的混合方法對控制方程進(jìn)行離散,并通過PJFNK算法對離散方程組進(jìn)行全耦合求解。

圖6 帶繞絲燃料組件典型示意圖

本文采用19棒束實(shí)驗(yàn)[31]對FLARE進(jìn)行驗(yàn)證,圖7為通道編號劃分以及測量通道示意圖。圖8展示了FLARE對組件出口歸一化溫度的計(jì)算結(jié)果以及與SUBAC程序[32]和Pronghorn-SC程序[33]的計(jì)算結(jié)果對比。由圖8可見,在高流量下軸向?qū)嵴贾鲗?dǎo),并由于幾何結(jié)構(gòu)不同,角、邊和中心通道的出口溫度存在著明顯的差異。隨著入口流量的減小,橫向?qū)嶂饾u發(fā)揮作用,各通道的出口溫度逐漸趨于一致。從圖8可看到,FLARE的計(jì)算結(jié)果與實(shí)驗(yàn)值符合良好,相對誤差在10%以內(nèi),并且與SUBAC以及Pronghorn-SC計(jì)算結(jié)果接近,說明FLARE能夠?qū)σ簯B(tài)金屬快堆堆芯內(nèi)的熱工水力參數(shù)進(jìn)行準(zhǔn)確計(jì)算。

圖7 棒束實(shí)驗(yàn)通道編號示意圖

a——高流量工況;b——中流量工況;c——低流量工況

2 燃料性能分析技術(shù)研究

核反應(yīng)堆燃料作為反應(yīng)堆內(nèi)裂變場所以及阻擋裂變產(chǎn)物釋放的第1道屏障,開展其服役性能分析與評估至關(guān)重要。傳統(tǒng)的燃料性能分析程序大多采用1.5維的有限差分方法開發(fā),受限于方法的適用性,此類程序的研究對象僅限于棒狀燃料?？紤]到有限元方法在固體力學(xué)方面的廣泛應(yīng)用及其建模對象的靈活性,XJTU-NuTHeL基于有限元開發(fā)了燃料性能分析程序BEEs。該程序不僅具備針對傳統(tǒng)棒狀燃料在穩(wěn)、瞬態(tài)下的分析功能[34-35],還能夠針對事故容錯燃料(Accident Tolerant Fuels, ATF)[20]以及其他幾何形狀(如環(huán)形[22]、板狀[21,36]等)的燃料開展多物理場耦合分析。

2.1 包覆顆粒彌散燃料性能分析

包覆顆粒彌散燃料是眾多ATF燃料選型方案之一,結(jié)構(gòu)包含三向同性TRISO(TRI-structural ISOtropic fuel)球形包覆顆粒、NITE-SiC(Nano-Infiltration and Transient Eutectic-phase Silicon Carbide)基體,如圖9所示。TRISO顆粒包含UO2芯核和4層包覆層,各包覆層都有特殊的功能。第1層為疏松熱解碳層(Buffer),能夠?yàn)闅鈶B(tài)裂變產(chǎn)物提供儲存空間并緩沖芯核的熱膨脹與輻照腫脹;第2層為內(nèi)致密熱解碳層(Inner Pyrolytic Carbon, IPyC),主要是防止裂變產(chǎn)物對SiC層的侵蝕,并承受部分內(nèi)壓;第3層為化學(xué)氣相沉積法制成的碳化硅層(Chemical Vapor Deposited Silicon Carbide, CVD-SiC),其阻擋固態(tài)裂變產(chǎn)物的能力較強(qiáng),能夠大幅度改善燃料滯留裂變產(chǎn)物的能力。此外,SiC的高導(dǎo)熱性和比熱容也能顯著降低燃料的峰值溫度。第4層為外致密熱解碳層(Outer Pyrolytic Carbon, OPyC),與IPyC材質(zhì)相同,主要作用是保護(hù)承壓的SiC層,并且在SiC失效時阻擋裂變產(chǎn)物的釋放。包覆顆粒彌散燃料具有良好的導(dǎo)熱性能和耐輻照性能,并能有效阻止裂變產(chǎn)物的擴(kuò)散。

圖9 TRISO包覆顆粒彌散燃料形式

包覆顆粒彌散燃料結(jié)構(gòu)極其復(fù)雜,高燃料體積份額的燃料芯塊可包含上千個TRISO顆粒,而每個TRISO顆粒具有多層結(jié)構(gòu),給開展其燃料模擬帶來了較大挑戰(zhàn)性。此外,對包覆顆粒彌散燃料進(jìn)行行為模擬需要處理復(fù)雜的非線性模型,同時考慮到輻照、熱、力多場耦合作用,增加了分析的復(fù)雜性。XJTU-NuTHeL基于有限元方法開發(fā)了包覆顆粒彌散燃料性能分析程序,建立了燃料芯塊-代表性體積元-TRISO顆粒的多尺度模擬方法(圖10)?？紤]到TRISO顆粒對NITE-SiC基體的變形影響可忽略,模擬過程中首先對移除顆粒的芯塊總體進(jìn)行計(jì)算并基于溫度分布等結(jié)果劃分子區(qū)域,再從子區(qū)域中選取代表性體積單元進(jìn)行計(jì)算分析,對其中TRISO顆粒性能進(jìn)行精細(xì)評估。

圖10 燃料性能多尺度模擬方法示意圖

XJTU-NuTHeL選用國際原子能機(jī)構(gòu)(IAEA)針對高溫氣冷堆公開的基準(zhǔn)題(IAEA Coordinated Research Program, CRP-6)[37]對程序進(jìn)行了系列驗(yàn)證。其中基準(zhǔn)題案例CASE-8是輻照過程中溫度和內(nèi)壓周期變化條件下的TRISO燃料顆粒,因涉及輻照-熱-力耦合而被國際上眾多程序選為驗(yàn)證算例。圖11展示了IPyC和CVD-SiC層內(nèi)側(cè)的切向應(yīng)力結(jié)果,由圖11可知,本程序的計(jì)算結(jié)果與國際上其他程序的計(jì)算結(jié)果[38-40]符合良好,驗(yàn)證了程序進(jìn)行燃料模擬的能力。

a——IPyC層內(nèi)側(cè)的切向應(yīng)力;b——CVD-SiC層內(nèi)側(cè)的切向應(yīng)力

圖12示出周向非均勻換熱條件下包覆顆粒彌散燃料性能分析的典型結(jié)果。程序?qū)崿F(xiàn)了燃料芯塊-代表性體積元-TRISO顆粒的多尺度熱力學(xué)三維精細(xì)模擬(圖12a～c),獲取了燃料各部分的溫度和應(yīng)力分布情況。燃料芯塊總體溫度受外表面周向非均勻換熱影響而呈現(xiàn)“偏心”分布。由于熱解碳材料隨中子輻照而產(chǎn)生收縮,CVD-SiC層在IPyC和OPyC的輻照變形下受壓,同時隨著燃耗的加深在Buffer層和IPyC層之間會出現(xiàn)間隙。此外,程序?qū)崿F(xiàn)了針對包覆顆粒彌散燃料的失效評估,得到燃料失效概率的變化情況(圖12d)。

a——芯塊總體溫度;b——代表性體積元應(yīng)力;c——TRISO應(yīng)力;d——燃料失效概率

2.2 板狀燃料元件性能分析

板狀燃料元件具有結(jié)構(gòu)緊湊、功率密度高、換熱能力強(qiáng)等特點(diǎn),在研究堆、船用動力堆等堆型具有廣泛應(yīng)用。同時,板狀燃料元件包殼與芯體直接接觸,對溫度、應(yīng)力應(yīng)變演變及裂變氣體釋放等均有影響。根據(jù)燃料元件芯體的結(jié)構(gòu),可將板狀燃料分類為單片式燃料元件和彌散型燃料元件。其中,單片式燃料元件采用整片的陶瓷燃料板作為燃料芯體,其芯體厚度一般較薄;彌散型燃料元件則采用陶瓷燃料顆粒彌散在金屬基體中的復(fù)合材料作為燃料芯體。

針對板狀燃料的上述特點(diǎn),XJTU-NuTHeL基于有限元平臺建立了燃料元件的材料物性和行為模型,采用全耦合方法實(shí)現(xiàn)了燃料元件性能分析[41]。為了準(zhǔn)確模擬彌散型芯體復(fù)合結(jié)構(gòu),建立了一系列等效物性模型,包括等效熱導(dǎo)率、等效熱膨脹系數(shù)、等效楊氏模量、等效泊松比、等效本征應(yīng)變等。圖13示出壽期初和壽期末的燃料元件中心截面處的燃耗、溫度和Von Mises應(yīng)力分布。在壽期初,燃料元件溫度分布是力學(xué)變形的主導(dǎo)因素,在功率分布和冷卻劑自下而上的對流換熱的綜合影響下,元件上部溫度較高,熱膨脹劇烈,導(dǎo)致壽期初中上部應(yīng)力較大。而在燃耗較深的情況下,長期輻照造成的本征應(yīng)變在力學(xué)方面逐漸占據(jù)了主導(dǎo)作用,元件中部燃耗最深的部位產(chǎn)生了很大的輻照腫脹,導(dǎo)致壽期末燃料元件中部的應(yīng)力顯著上升。

a——壽期初;b——壽期末

在燃料元件及輻照-力-熱耦合分析的基礎(chǔ)上,開發(fā)鋯合金的氧化和吸氫模型,實(shí)現(xiàn)了核反應(yīng)堆惡劣環(huán)境下燃料包殼的氧化和吸氫模擬。針對1/4燃料元件模擬獲得的代表性計(jì)算結(jié)果如圖14所示。氧化層厚度分布主要受溫度影響,在元件中部最厚。吸氫首先在氧化部位發(fā)生,然后沿濃度梯度和溫度梯度擴(kuò)散。最終燃料包殼邊緣處氫濃度達(dá)到飽和固溶氫,并析出形成氫化物。

圖14 燃料包殼氧化和吸氫模擬計(jì)算結(jié)果

3 多物理場耦合框架網(wǎng)格映射技術(shù)方案及評估

隨著高精度數(shù)值算法、大規(guī)模并行計(jì)算及計(jì)算機(jī)軟件技術(shù)的發(fā)展,實(shí)現(xiàn)反應(yīng)堆不同系統(tǒng)部件精細(xì)建模及多物理場分析,獲取高保真計(jì)算結(jié)果成為可能。當(dāng)前,國內(nèi)外已開發(fā)多個多物理場耦合平臺用于反應(yīng)堆多場耦合計(jì)算,如MOOSE[42-43]、SALOME[44-45]、MpCCI[46]、COMSOL[47]、LIME[48]和CVR[7,49]等。傳統(tǒng)的直接耦合方式存在著數(shù)據(jù)交換效率偏低、耦合程序需采用相同網(wǎng)格等不足,需要建立多物理場耦合框架,采用準(zhǔn)確高效的網(wǎng)格映射算法,實(shí)現(xiàn)反應(yīng)堆多物理場耦合問題在統(tǒng)一框架下的求解。為此,XJTU-NuTHeL基于MOOSE/LIBMESH平臺搭建了堆芯多物理場耦合框架,用于實(shí)現(xiàn)反應(yīng)堆的高保真多物理場耦合計(jì)算。

網(wǎng)格映射技術(shù)主要包括兩方面內(nèi)容,一方面是外部程序與多物理場耦合框架之間的網(wǎng)格映射,另一方面是框架內(nèi)程序之間的網(wǎng)格映射。在外部程序數(shù)據(jù)映射方面,本文采用了逐一映射的方式,即在框架內(nèi)構(gòu)建與外部程序采用相同網(wǎng)格劃分形式的數(shù)據(jù)存儲網(wǎng)格,之后將外部程序的計(jì)算結(jié)果根據(jù)點(diǎn)/單元逐個映射到網(wǎng)格中。網(wǎng)格映射方案如圖15所示。

a——求解域;b——外部求解離散方式及值分布;c——有限元網(wǎng)格及值分布

框架內(nèi)程序之間的網(wǎng)格映射采用網(wǎng)格映射算法實(shí)現(xiàn)。為了獲得堆芯多物理場耦合計(jì)算中高效映射算法,本文采用單塊板狀燃料作為對象,研究了不同算法的傳遞誤差和在千萬級-百萬級網(wǎng)格之間的網(wǎng)格映射效率。網(wǎng)格映射算法的效率評估結(jié)果列于表1。網(wǎng)格映射算法效率從高到低依次為:反距離權(quán)重插值[50]、網(wǎng)格函數(shù)傳遞[51]、L2映射[52]、最近鄰傳遞。需要說明的是,為了避免計(jì)算資源對效率評估結(jié)果的影響,本研究中對每種算法所采用的計(jì)算資源均為該方法在超算節(jié)點(diǎn)上測試所需的最少計(jì)算資源,超算單個節(jié)點(diǎn)的核數(shù)為64。

表1 網(wǎng)格映射算法效率評估結(jié)果

在傳遞精度方面,同樣選取板狀燃料作為研究對象,假設(shè)傳遞的變量分別為余弦形式的函數(shù),源網(wǎng)格為六面體網(wǎng)格(HEX8),目標(biāo)網(wǎng)格為四面體網(wǎng)格(TET4),源網(wǎng)格與目標(biāo)網(wǎng)格節(jié)點(diǎn)數(shù)分別為6 080和65 567。不同映射算法的傳遞誤差列于表2。其中,L2誤差為誤差平方之和,H1誤差為變量梯度誤差平方之和。由表2可見,網(wǎng)格函數(shù)傳遞算法在該類網(wǎng)格上具有最小的傳遞誤差,L2映射次之,最近鄰傳遞的誤差最大。綜合考慮,網(wǎng)格函數(shù)傳遞算法能夠兼顧精度和效率,在結(jié)構(gòu)化網(wǎng)格中映射變量時更具優(yōu)勢。

表2 不同映射算法的傳遞誤差

4 堆芯多物理場耦合框架搭建及典型應(yīng)用

基于上述映射算法,XJTU-NuTHeL開發(fā)了堆芯多物理場耦合框架,并實(shí)現(xiàn)了板狀燃料性能分析程序BEEs-Plates[21,36]、蒙特卡羅中子物理程序OpenMC[53]和一維流體域分析程序NUSAC[18,27]在框架內(nèi)的集成,初步實(shí)現(xiàn)了開源CFD軟件OpenFOAM[54]的耦合。除此之外,在該框架中還實(shí)現(xiàn)了由中國核動力研究設(shè)計(jì)院開發(fā)的熱工水力程序CORTH[55]和中子學(xué)程序CORCA[56]的耦合計(jì)算。框架的結(jié)構(gòu)如圖16所示。

基于BEEs-Plates/OpenMC/NUSAC耦合程序[23]開展了JRR-3研究堆組件級核-熱-力耦合計(jì)算。為了考慮燃耗對燃料性能的影響,在每一計(jì)算時刻將預(yù)先完成的燃耗計(jì)算中的結(jié)果映射到網(wǎng)格中,之后傳遞給BEEs-Plates進(jìn)行燃料性能計(jì)算。程序采用Picard迭代方式進(jìn)行求解,程序結(jié)構(gòu)如圖17所示,計(jì)算總時長為2.0×107s。對25%富集度的U3Si2-Al彌散燃料組件的中子物理計(jì)算結(jié)果如圖18所示?？梢钥吹?由于計(jì)算中在組件長度(z軸方向)和寬度方向(x軸方向)設(shè)置了反射邊界條件,組件邊緣區(qū)域的功率密度顯著高于中心區(qū)域。由于235U消耗速度較快,功率越高的區(qū)域鈾原子密度越小。

圖17 BEEs-Plates/OpenMC/NUSAC耦合程序結(jié)構(gòu)及數(shù)據(jù)傳遞的方向

圖18 組件功率密度(a)和235U原子密度分布(b)

部分燃料性能參數(shù)的計(jì)算結(jié)果如圖19所示。由圖19a可看出,在組件包殼與連接部件接觸區(qū)域容易發(fā)生應(yīng)力集中現(xiàn)象。當(dāng)組件的平均燃耗達(dá)到64.77 GW·d/tU時,組件最大應(yīng)力達(dá)到329.96 MPa。從圖19b中體積應(yīng)變的分布結(jié)果可知,組件的最大體積應(yīng)變?yōu)?0%,約為最小體積應(yīng)變的10倍。此外,與截面上功率的分布趨勢一樣,組件的體積應(yīng)變分布也呈現(xiàn)出中心低、邊緣高的趨勢,說明溫度對燃料的應(yīng)變影響較大。

圖19 組件應(yīng)力(a)和體積應(yīng)變(b)分布

5 結(jié)論

針對核反應(yīng)堆高精度數(shù)值模擬及多場耦合分析需求,基于開源有限元計(jì)算框架MOOSE,XJTU-NuTHeL團(tuán)隊(duì)完成了相關(guān)熱工流體計(jì)算模型開發(fā),構(gòu)建了系統(tǒng)分析程序NUSAC,其中一維單相有限元離散計(jì)算依靠相關(guān)穩(wěn)定性算法,一維兩相計(jì)算采用高精度間斷有限元離散,實(shí)現(xiàn)了核反應(yīng)堆系統(tǒng)相關(guān)穩(wěn)態(tài)及瞬態(tài)計(jì)算;構(gòu)建了子通道分析程序FLARE,可高精度準(zhǔn)確預(yù)測液態(tài)金屬快堆堆芯內(nèi)的熱工水力參數(shù)。開展了相關(guān)燃料性能分析技術(shù)研究,構(gòu)建了燃料性能分析程序BEEs,實(shí)現(xiàn)了涉及包覆顆粒彌散、板狀等多種類型、多類工況的燃料高維度精細(xì)化模擬。搭建了多物理場耦合框架,測試評估框架內(nèi)不同物理場求解程序的網(wǎng)格映射與數(shù)據(jù)傳遞算法,實(shí)現(xiàn)了板狀燃料組件級核-熱-力耦合計(jì)算。

綜上所述,針對核反應(yīng)堆多維度多物理場耦合分析的有限元技術(shù)可以重構(gòu)傳統(tǒng)核反應(yīng)堆熱工水力分析程序,實(shí)現(xiàn)在數(shù)值算法、離散精度等方面的優(yōu)化,為核反應(yīng)堆多維度耦合問題奠定了基礎(chǔ);可以優(yōu)化傳統(tǒng)燃料性能分析技術(shù),突破傳統(tǒng)的棒狀燃料元件的研究范圍,充分應(yīng)用到各類先進(jìn)燃料概念的開發(fā)驗(yàn)證中,高效解決了核反應(yīng)堆工程中的各類精細(xì)化模擬問題;可以為多物理場耦合問題提供合理統(tǒng)一的計(jì)算框架,通過在不同物理場程序間提供的網(wǎng)格映射與數(shù)據(jù)傳遞接口,充分考慮不同物理場之間的影響,有助于更好地理解和研究核反應(yīng)堆復(fù)雜多場環(huán)境下的相關(guān)現(xiàn)象。