苗建新，曹良志,*，方 超，賀清明，曹啟祥，尹 苗

(1.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049；2.核工業(yè)西南物理研究院，四川 成都 610041)

采用球諧函數(shù)展開中子輸運(yùn)方程的角度變量，可避免離散縱標(biāo)方法中出現(xiàn)的射線效應(yīng)；采用有限元方法離散輸運(yùn)方程的空間變量，可減少傳統(tǒng)結(jié)構(gòu)網(wǎng)格對復(fù)雜幾何建模上的近似；將兩種方法結(jié)合起來求解中子輸運(yùn)方程可兼具二者的優(yōu)點(diǎn)。標(biāo)準(zhǔn)的Galerkin有限元方法適用于橢圓方程的求解，如二階中子輸運(yùn)方程，擴(kuò)散方程等；對于一階中子輸運(yùn)方程這種雙曲方程，直接采用Galerkin有限元方法會(huì)出現(xiàn)非物理的震蕩現(xiàn)象。國內(nèi)外針對球諧函數(shù)有限元方法的研究主要求解的是二階自共軛方程[1]和二階偶對稱方程[2]，但由于截面會(huì)出現(xiàn)在方程分母上，對于包含真空區(qū)域的問題無法進(jìn)行求解，限制了方法的應(yīng)用場景。Buchan等[3]針對求解一階中子輸運(yùn)方程的球諧函數(shù)有限元方法進(jìn)行了研究，解決了數(shù)值不穩(wěn)定現(xiàn)象，但目前還未應(yīng)用于三維大規(guī)模問題。

包層作為聚變堆中的重要結(jié)構(gòu)，承擔(dān)著氚增殖、能量轉(zhuǎn)換和輻射屏蔽等一系列功能。因此對包層的中子學(xué)分析顯得尤為重要。聚變堆包層由第一壁、上下蓋板、增殖單元、冷卻劑管道以及背板等結(jié)構(gòu)組成，存在較多的非規(guī)則幾何結(jié)構(gòu)，如采用基于結(jié)構(gòu)網(wǎng)格的確定論方法進(jìn)行計(jì)算要引入大量的近似，計(jì)算精度較低。蒙特卡羅方法雖能描述復(fù)雜的幾何結(jié)構(gòu)，但計(jì)算耗時(shí)較長。由于有限元方法對于復(fù)雜幾何的處理具備天然的優(yōu)勢，近年來，基于有限元方法的確定論程序開始應(yīng)用于包層中子學(xué)分析。如美國的商用軟件ATTILA[4]采用離散縱標(biāo)和間斷有限元方法，已應(yīng)用于國際熱核聚變實(shí)驗(yàn)堆(ITER)的中子學(xué)分析中。韓國的AETIUS[5]與ATTILA采用相同的方法，已應(yīng)用于韓國氦冷陶瓷試驗(yàn)包層的中子學(xué)分析中。目前國內(nèi)進(jìn)行聚變堆包層中子學(xué)計(jì)算大多采用蒙特卡羅程序，還未開發(fā)能應(yīng)用于包層中子學(xué)分析的有限元程序。

本文研究求解一階中子輸運(yùn)方程的球諧函數(shù)穩(wěn)定有限元方法，并基于此方法開發(fā)中子學(xué)分析程序NECP-FISH。最終將其應(yīng)用于氦冷陶瓷包層，并與蒙特卡羅程序NECP-MCX[6]的計(jì)算結(jié)果進(jìn)行對比分析。

1 球諧函數(shù)有限元方法

NECP-FISH程序采用球諧函數(shù)和多尺度穩(wěn)定有限元方法求解一階穩(wěn)態(tài)中子輸運(yùn)方程。

不考慮裂變材料，中子輸運(yùn)方程右端只有外源項(xiàng)和散射源項(xiàng)，多群輸運(yùn)方程不同能群通過散射源項(xiàng)進(jìn)行耦合，此時(shí)方程可寫為單能形式：

(1)

(2)

根據(jù)加權(quán)殘差法獲得一階輸運(yùn)方程的弱形式：

〈w,φ〉?V+(L*w,φ)V=(w,s)V

(3)

本文所采用的多尺度有限元方法是子網(wǎng)格有限元方法中的一種。子網(wǎng)格有限元方法最早由Hughes[7]提出，Buchan在2008年第1次將其應(yīng)用到求解中子輸運(yùn)方程中[3]。此方法最重要的思想就是將未知量分解為兩部分：

φ=φc+φr

(4)

式中：第1部分φc為連續(xù)可解部分，定義在全局尺度上；第2部分φr為殘差部分，這部分在最終離散的方程中并不進(jìn)行求解。根據(jù)上述對未知量的分解，同時(shí)由于殘差部分定義在每個(gè)單元上，邊界條件只應(yīng)用于連續(xù)部分，所以式(3)的弱形式可寫為：

〈w,φc〉?V+(L*w,φc+φr)=(w,s)

(5)

同時(shí)殘差方程可表示為：

Lφr=s-Lφc

(6)

此殘差方程可近似求解，使用不同的近似方法產(chǎn)生了不同的子網(wǎng)格有限元方法。本文中的代數(shù)子網(wǎng)格方法(algebraic sub-grid scale method)采用代數(shù)方程對通量密度的殘差部分進(jìn)行近似：

φr=λ(s-Lφc)

(7)

式中，λ相當(dāng)于對L-1的近似。將式(7)代入到式(5)中，可獲得代數(shù)子網(wǎng)格方法的弱形式為：

〈φ,φc〉+(L*φ,φc)-(L*φ,λLφc)=

(φ,s)-(L*φ,λs)

(8)

上式中的未知量φc可由基函數(shù)展開φ(r,Ω)=φT(r,Ω)φ，此基函數(shù)可用分段多項(xiàng)式函數(shù)N(r)與球諧函數(shù)Y(Ω)的克羅內(nèi)克積的形式來表示：

φ(r,Ω)=N(r)?Y(Ω)

(9)

最終方程的離散形式可寫為：

(10)

在每個(gè)內(nèi)部單元和邊界單元上，式(10)中左端的每項(xiàng)均可計(jì)算得到并組裝到整體剛度矩陣中，最終通過求解線性方程組來獲得每個(gè)網(wǎng)格節(jié)點(diǎn)的通量密度。對此球諧函數(shù)有限元方法的詳細(xì)推導(dǎo)見參考文獻(xiàn)[8]。

2 程序開發(fā)及功能介紹

基于上述理論方法，并借助開源的前后處理平臺SALOME[9]，開發(fā)了中子學(xué)分析程序NECP-FISH，程序的整體框架如圖1所示，該程序具備從可視化建模、中子學(xué)分析到結(jié)果可視化展示的一系列功能。其中中子學(xué)分析涵蓋了中子通量密度、氚增殖比、中子釋熱的計(jì)算，主要由球諧函數(shù)有限元求解器實(shí)現(xiàn)。借助SALOME平臺開發(fā)程序的優(yōu)點(diǎn)主要是：首先，SALOME平臺集成了很多開源工具(如實(shí)現(xiàn)CAD建模的OpenCASCADE，網(wǎng)格剖分的Gmsh和NETGEN，可視化的ParaView)，借助其開發(fā)程序，可使程序的功能更完備；其次，SALOME平臺提供了用戶二次開發(fā)的接口，用戶通過開發(fā)接口程序，可方便地將自己的數(shù)值計(jì)算程序接入到平臺中形成完整的分析程序，可實(shí)現(xiàn)數(shù)據(jù)的內(nèi)部傳遞，方便用戶操作；最后，SALOME平臺的使用也較為簡單，它在Windows和Linux系統(tǒng)下均可使用。借助其開發(fā)的程序兼容性也較好。但另一方面，由于平臺集成的多為開源工具，相比于提供前后處理功能的商用軟件，在功能上還有一定欠缺。

圖1 NECP-FISH程序框架

2.1 可視化建模

程序的可視化建模通過SALOME平臺集成的CAD和網(wǎng)格剖分工具實(shí)現(xiàn)。程序可直接導(dǎo)入設(shè)計(jì)好的幾何文件進(jìn)行建模，也可對分開導(dǎo)入的幾何模型進(jìn)行裝配從而建立最終的計(jì)算模型。對于結(jié)構(gòu)簡單的基準(zhǔn)題，用戶還可直接構(gòu)建幾何體再通過布爾運(yùn)算進(jìn)行建模。對于建立好的幾何模型，NECP-FISH程序可采用不同的網(wǎng)格剖分方法，目前程序有限元求解器支持的非結(jié)構(gòu)網(wǎng)格類型為二維三角形和三維四面體網(wǎng)格。利用程序建立的CFETR氦冷陶瓷包層的中子學(xué)分析模型如圖2所示。

圖2 NECP-FISH程序建立的模型

上述的建模和剖分網(wǎng)格操作均可通過程序提供的用戶圖形界面可視化完成，同時(shí)，為更方便完成中子學(xué)計(jì)算，NECP-FISH程序在SALOME平臺中采用Python結(jié)合Qt庫的方式開發(fā)了接口程序，通過接口程序的可視化界面，將多群截面、求解參數(shù)、邊界條件、中子源以及剖分出的網(wǎng)格存儲(chǔ)在HDF5格式的輸入文件中。相比于傳統(tǒng)程序的文本輸入，HDF5文件的特點(diǎn)是可分級分類存儲(chǔ)數(shù)據(jù)，其存儲(chǔ)相同大小的數(shù)據(jù)量所占的內(nèi)存空間更小。最終生成的輸入卡片可直接提供給有限元求解器進(jìn)行計(jì)算。

2.2 中子學(xué)分析

程序的中子學(xué)分析功能是由球諧函數(shù)有限元求解器(PN-FEM)實(shí)現(xiàn)的，有限元求解器可計(jì)算出各網(wǎng)格節(jié)點(diǎn)上的中子通量密度，再根據(jù)氚增殖比以及釋熱率的計(jì)算公式可給出相應(yīng)物理量的計(jì)算結(jié)果。

1) 有限元求解器

NECP-FISH程序的核心是有限元求解器。求解器是基于函數(shù)庫進(jìn)行開發(fā)的，圖3所示為求解器的整體框架。Fortran通用函數(shù)庫和輸入輸出模塊將輸入文件的內(nèi)容讀入，再通過球諧函數(shù)模塊、有限元模塊、材料截面和源模塊構(gòu)造出要求解的問題。問題的求解主要是由開源數(shù)學(xué)庫psblas[10]組成的線性代數(shù)求解模塊實(shí)現(xiàn)的。通過建立不同模塊的函數(shù)庫，可使程序結(jié)構(gòu)清晰易于理解，同時(shí)基于函數(shù)庫開發(fā)程序也更快捷。

圖3 求解器的整體框架

求解器的中子學(xué)計(jì)算流程如圖4所示，整體可分為預(yù)處理、求解和后處理。在預(yù)處理過程中讀取輸入并根據(jù)計(jì)算核數(shù)決定是否進(jìn)行空間區(qū)域分解。求解過程主要是進(jìn)行多群迭代，多群迭代是能量由高到低依次對每個(gè)能群進(jìn)行求解。由于計(jì)算的問題是屏蔽問題，所以不存在裂變源迭代。同時(shí)由于群內(nèi)散射移到了方程左端，所以求解過程中不需要散射源迭代。程序利用前面能群計(jì)算出的通量密度矩更新散射源后即可進(jìn)入單群方程求解，之后在每個(gè)單元上計(jì)算單元矩陣和向量再組裝到整體剛度矩陣和向量中。組裝好的線性代數(shù)方程組采用廣義極小殘余算法(GMRES)進(jìn)行迭代求解。對于沒有上散射的問題，只需1次能群掃描，對于包含上散射的問題，需從存在上散射的最高能群進(jìn)行多次能群掃描，直至多群通量密度收斂。最后在后處理模塊中實(shí)現(xiàn)氚增殖比、中光子釋熱等物理量的計(jì)算，并生成可視化文件。

圖4 求解器的中子學(xué)計(jì)算流程圖

2) 程序空間并行

隨著網(wǎng)格數(shù)的增多，計(jì)算規(guī)模迅速增大，因此需利用高性能計(jì)算平臺并行計(jì)算以提高計(jì)算效率。由于球諧函數(shù)很難實(shí)現(xiàn)各角度矩的解耦，所以程序只采用空間并行。并行策略是基于消息傳遞接口(MPI)的分布式并行。并行求解示意圖如圖5所示。首先在預(yù)處理過程中進(jìn)行空間區(qū)域分解，非結(jié)構(gòu)網(wǎng)格的空間區(qū)域分解算法十分復(fù)雜，通用的做法是使用專門的程序包。NECP-FISH程序采用的是METIS程序包，將非結(jié)構(gòu)網(wǎng)格各節(jié)點(diǎn)按特定編號排列好，再調(diào)用相應(yīng)函數(shù)即完成區(qū)域分解。這樣就將整個(gè)計(jì)算問題劃分為幾個(gè)小的區(qū)域分配到各核上，并明確了核與核之間要通信的網(wǎng)格節(jié)點(diǎn)編號。區(qū)域分解后相當(dāng)于在每個(gè)計(jì)算核心上組裝并求解整體線性方程組的一部分，迭代求解過程中的矩陣和向量相乘以及向量點(diǎn)積計(jì)算要進(jìn)行核間的通信。最終統(tǒng)一區(qū)域交界處各點(diǎn)的計(jì)算結(jié)果即完成并行計(jì)算。

圖5 程序并行求解示意圖

2.3 結(jié)果可視化

結(jié)果的可視化需生成存儲(chǔ)計(jì)算結(jié)果的特定格式文件，該文件包括網(wǎng)格信息、離散節(jié)點(diǎn)或單元上的計(jì)算結(jié)果。NECP-FISH程序可視化結(jié)果包括中子通量密度分布以及釋熱分布，采用XDMF格式。該格式采用XML格式存儲(chǔ)數(shù)據(jù)路徑，再讀取HDF5文件中的數(shù)據(jù)來完成結(jié)果的可視化。相比于傳統(tǒng)的VTK格式文件，XDMF格式的文件可并行輸出，同時(shí)數(shù)據(jù)存儲(chǔ)所占的內(nèi)存更小。文件可導(dǎo)入到SALOME平臺中的ParaVis模塊中查看，也可采用其他可視化軟件進(jìn)行查看。

3 氦冷陶瓷包層計(jì)算結(jié)果

將NECP-FISH程序應(yīng)用于聚變示范堆氦冷陶瓷包層的中子學(xué)分析中，此包層由核工業(yè)西南物理研究院設(shè)計(jì)[11]。由于包層內(nèi)部為增殖單元對稱排布，且增殖單元涵蓋了整個(gè)包層的全部材料，因此分析單個(gè)增殖單元即可反映整個(gè)包層的中子學(xué)特性。本文選取具有代表性的包層增殖單元進(jìn)行計(jì)算。圖6為增殖單元的材料布置和邊界條件，氚增殖劑為Li4SiO4，中子倍增劑為Be，兩者做成U型球床的結(jié)構(gòu)，中間用冷卻板間隔開。背板及冷卻板的材料均采用低活化的鐵素體鋼(CLF-1)。在第一壁前增加了2 mm的鎢涂層，同時(shí)在鎢涂層前定義了各向同性的均勻體源。

根據(jù)圖6的幾何布置和尺寸，利用NECP-FISH程序建立計(jì)算的三維模型，并剖分出非結(jié)構(gòu)的四面體單元，如圖7所示。

圖6 增殖單元的材料布置和邊界條件

a——建立的三維增殖單元模型；b——程序剖分的網(wǎng)格

分別采用NECP-FISH和NECP-MCX對此模型進(jìn)行計(jì)算。NECP-FISH的多群數(shù)據(jù)庫采用FENDL-2.1的MATXS格式數(shù)據(jù)庫[12]，并通過TRANSX處理獲得多群宏觀截面，能群結(jié)構(gòu)為Vitamin-J 175群的能群結(jié)構(gòu)。散射截面展開階數(shù)為P4，球諧函數(shù)展開階數(shù)為P15，剖分的四面體網(wǎng)格數(shù)為292 593。NECP-MCX采用FENDL-2.1的ACE格式數(shù)據(jù)庫，投入了100組，每組1 000 000個(gè)粒子。圖8為NECP-FISH計(jì)算出的中子通量密度的空間分布，隨著與中子源距離的增大，中子通量密度下降了2個(gè)數(shù)量級。

圖8 NECP-FISH計(jì)算獲得的中子通量密度空間分布

各區(qū)域的平均積分通量密度與NECP-MCX計(jì)算結(jié)果的對比如圖9所示，通量密度的偏差在距源較近的區(qū)域稍大，所有區(qū)域的相對偏差均在6%以內(nèi)。同時(shí)，圖中還給出了蒙特卡羅程序MCNP[13]的計(jì)算結(jié)果，兩蒙特卡羅程序的相對偏差均在0.4%以下。

圖9 各區(qū)域平均積分通量密度

對此問題進(jìn)行氚增殖比的計(jì)算，氚增殖比的定義為1個(gè)歸一的聚變中子在包層中產(chǎn)生的氚原子個(gè)數(shù)，產(chǎn)氚主要是通過中子與包層中的含鋰材料發(fā)生反應(yīng)實(shí)現(xiàn)的，因此NECP-FISH程序計(jì)算氚增殖比的公式為：

(11)

式中：Σt,Li4SiO4(E)為硅酸鋰材料的宏觀產(chǎn)氚截面；Sn為中子源強(qiáng)。NECP-FISH程序計(jì)算的此增殖單元的整體氚增殖比為1.339，NECP-MCX的計(jì)算結(jié)果為1.332，相對統(tǒng)計(jì)方差為0.006 7%。MCNP的計(jì)算結(jié)果為1.334。氚增殖比的分布及NECP-FISH與NECP-MCX程序計(jì)算結(jié)果的相對偏差列于表1。

表1 氚增殖比的計(jì)算結(jié)果

各區(qū)域中子釋熱的結(jié)果列于表2，中子釋熱是由通量密度與截面相乘再進(jìn)行積分獲得的：

表2 增殖單元各區(qū)域中子釋熱對比結(jié)果

Kn=?Σh(E)φ(r,E)drdE

(12)

式中，Σh(E)為宏觀釋熱截面。由上式可知，中子釋熱偏差與通量密度的偏差直接相關(guān)，NECP-FISH與NECP-MCX計(jì)算結(jié)果的最大偏差出現(xiàn)在靠近中子源的區(qū)域，所有區(qū)域的偏差都在6%以內(nèi)。在包層核熱計(jì)算方面，前文所提到的確定論程序ATTILA與蒙特卡羅程序最大相對偏差達(dá)到了28%[14]。

4 結(jié)論

本文研究了球諧函數(shù)有限元方法，并基于此方法開發(fā)了中子學(xué)分析程序NECP-FISH。該程序具備可視化建模、中子學(xué)計(jì)算和結(jié)果可視化功能。利用NECP-FISH對聚變堆氦冷陶瓷包層進(jìn)行了建模，并計(jì)算了中子通量密度、氚增殖比以及中子釋熱3個(gè)重要的中子學(xué)參數(shù)。結(jié)果表明，NECP-FISH與蒙特卡羅程序NECP-MCX計(jì)算的TBR的相對偏差為0.56%，各區(qū)域平均積分通量密度最大相對偏差為鎢涂層的-5.27%，其余區(qū)域均在5%以內(nèi)，中子釋熱最大相對偏差為-5.79%。結(jié)果證明了NECP-FISH程序能應(yīng)用于聚變堆包層中子學(xué)分析。某些區(qū)域偏差較大的原因可能是使用的多群數(shù)據(jù)庫與連續(xù)能量數(shù)據(jù)庫的差別，后續(xù)將進(jìn)行進(jìn)一步的研究。