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

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

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

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

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

1 球諧函數有限元方法

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

不考慮裂變材料，中子輸運方程右端只有外源項和散射源項，多群輸運方程不同能群通過散射源項進行耦合，此時方程可寫為單能形式：

(1)

(2)

根據加權殘差法獲得一階輸運方程的弱形式：

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

(3)

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

φ=φc+φr

(4)

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

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

(5)

同時殘差方程可表示為：

Lφr=s-Lφc

(6)

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

φr=λ(s-Lφc)

(7)

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

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

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

(8)

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

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

(9)

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

(10)

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

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

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

圖1 NECP-FISH程序框架

2.1 可視化建模

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

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

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

2.2 中子學分析

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

1) 有限元求解器

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

圖3 求解器的整體框架

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

圖4 求解器的中子學計算流程圖

2) 程序空間并行

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

圖5 程序并行求解示意圖

2.3 結果可視化

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

3 氦冷陶瓷包層計算結果

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

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

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

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

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

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

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

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

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

(11)

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

表1 氚增殖比的計算結果

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

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

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

(12)

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

4 結論

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