核反應堆系統多維度多物理場耦合有限元分析研究

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

(西安交通大學 核科學與技術學院,陜西省先進核能技術重點實驗室,陜西 西安 710049)

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

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

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

傳統系統分析程序與子通道程序框架固定、缺乏靈活的程序架構以及先進的數值算法和物理模型,而先進的多物理場耦合平臺可以實現各物理模型的靈活植入,充分發揮多物理場耦合的優勢,可彌補傳統分析程序在軟件結構及算法上的不足。因此,搭建多物理場耦合平臺,針對耦合問題中的關鍵技術開展研究,對加快我國自主化多物理場耦合平臺的開發進程具有重要意義。

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

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

1 熱工流體計算模型開發

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

1.1 核反應堆系統兩相流分析模型開發

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

(1)

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

在兩相流模型中,NUSAC采用了一維兩流體六方程模型和經過廣泛驗證的與RELAP5相同的本構模型,其中兩相流守恒方程為:

(2)

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

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

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

圖1 水龍頭問題驗證[19]

圖3 FRIGG飽和沸騰實驗驗證[19]

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

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

a——冷卻劑進出口溫度;b——包殼最高溫度

1.2 液態金屬快堆子通道分析模型開發

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

FLARE的控制方程如下。

質量方程:

(3)

軸向動量方程:

(4)

橫向動量方程:

(5)

能量方程:

(6)

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

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

本文采用19棒束實驗[31]對FLARE進行驗證,圖7為通道編號劃分以及測量通道示意圖。圖8展示了FLARE對組件出口歸一化溫度的計算結果以及與SUBAC程序[32]和Pronghorn-SC程序[33]的計算結果對比。由圖8可見,在高流量下軸向導熱占主導,并由于幾何結構不同,角、邊和中心通道的出口溫度存在著明顯的差異。隨著入口流量的減小,橫向導熱逐漸發揮作用,各通道的出口溫度逐漸趨于一致。從圖8可看到,FLARE的計算結果與實驗值符合良好,相對誤差在10%以內,并且與SUBAC以及Pronghorn-SC計算結果接近,說明FLARE能夠對液態金屬快堆堆芯內的熱工水力參數進行準確計算。

圖7 棒束實驗通道編號示意圖

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

2 燃料性能分析技術研究

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

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

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

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

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

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

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

a——IPyC層內側的切向應力;b——CVD-SiC層內側的切向應力

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

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

2.2 板狀燃料元件性能分析

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

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

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

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

圖14 燃料包殼氧化和吸氫模擬計算結果

3 多物理場耦合框架網格映射技術方案及評估

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

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

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

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

表1 網格映射算法效率評估結果

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

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

4 堆芯多物理場耦合框架搭建及典型應用

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

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

圖17 BEEs-Plates/OpenMC/NUSAC耦合程序結構及數據傳遞的方向

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

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

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

5 結論

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

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