中子輸運計算集成軟件平臺的研究與設計

李 磊, 王世慶, 李 偉, 柳 建

(1.核工業(yè)西南物理研究院, 四川 成都 610225;2.中廣核研究院有限公司, 廣東 深圳 518000)

在中子輸運分析中，主要的方法有確定論方法以及非確定論方法[1]，其中非確定論方法即隨機模擬方法,具有適用復雜幾何機構、能進行精細模擬等優(yōu)勢，隨著現(xiàn)代計算機計算能力的大幅度提升，在核工程領域得到越來越多的應用[2]。由美國阿拉莫斯實驗室(LANL)開發(fā)的蒙特卡羅程序MCNP是隨機模擬程序中的佼佼者，已廣泛應用于輻射防護、反應堆設計、臨界裝置實驗等領域[3]。但是MCNP輸入卡使用文本結構，其內(nèi)容全部需要人工輸入，效率低并且容易出現(xiàn)輸入卡編寫錯誤，對工程人員的經(jīng)驗要求較高，在復雜幾何模型中，這些問題更加凸顯。

為解決MCNP輸入卡復雜、效率低、易出錯的問題，一些學者對MCNP的建模自動化進行了研究，并且開發(fā)了特定功能的應用工具。從最早的Visual Editor、Sabrina和Moritz[4]等模型可視化工具到MCAM[5]，極大地提高了程序使用的便利性。這些工具使用的方法也有多種，例如羅月童等人提出的BREP到半空間轉(zhuǎn)化算法[6]，張建生等人提出的UG模型空間樹轉(zhuǎn)化算法[7]，王寒冰提出的基于特征的BREP到CSG模型轉(zhuǎn)換算法[8]以及吳烔等人運用卷積神經(jīng)網(wǎng)絡(CNN)對CAonSTEP算法進行改進的ICAonSTEPS算法[9]等。

上面提到MCNP輔助工具都是基于大型三維建模軟件的二次開發(fā)，屬于重量級工具。隨著“第四代核能系統(tǒng)”研究不斷深入，需要對不同堆芯模型進行設計驗證，為滿足核反堆數(shù)字化設計需求[10]，將更多工作交由計算機完成。因此，亟需開發(fā)一套集成MCNP輸入卡自動生成、自動調(diào)用計算核心程序及計算結果后處理等功能于一體的輕量級軟件平臺，以降低MCNP程序使用門檻和時間成本，提高工作效率。本文工作依托于國家重點研發(fā)計劃“核安全與先進核能技術”重點專項，基于堆芯模型快速搭建、MCNP靈活輸入以及計算結果提取分析等功能，實現(xiàn)反應堆堆芯物理特性的快速分析[11]，提高項目研發(fā)效率。

1 集成軟件平臺的功能及關鍵模塊

集成軟件平臺主要由三部分組成，包括輸入卡自動生成模塊、計算核心調(diào)用模塊以及結果處理模塊。在整個中子輸運計算過程中，前后處理占了60%～80%的時間[12]，包括幾何模型構建、計算參數(shù)設置以及結果可視化分析等，相比于傳統(tǒng)的單獨處理方法，將其功能模塊集成后能節(jié)省大量時間。輸入卡生成模塊將實體幾何模型轉(zhuǎn)化為MCNP輸入卡中的幾何模型信息，再補全剩余卡片信息生成完整輸入卡文件。計算核心調(diào)用模塊中將上一模塊生成或外部空間導入的輸入卡文件使用MCNP程序進行計算，得到計算結果文件。結果處理模塊中進行計算結果的數(shù)據(jù)處理和分析，其處理的結果文件可以是計算模塊中產(chǎn)生的，也可以從外部空間中導入。整個集成軟件平臺的功能流程圖如圖1所示。下面對其中的幾個關鍵技術進行闡述。

圖1 集成軟件平臺的功能流程框圖Fig.1 The function flow chart of the integrated software platform

1.1 BREP文件到MCNP幾何輸入卡的自動生成

經(jīng)過對多種實體幾何模型文件類型的比較，最終選擇了可讀性高的邊界表示法(BREP)文件作為實體幾何模型源文件。BREP文件使用邊界來表示實體幾何模型，通過基本幾何元素(點、線、面、體等)來存儲幾何信息，同時依據(jù)已知的拓撲關系(體→面→環(huán)→邊→點)來構建各基本幾何元素間的連接關系，進而實現(xiàn)對實體模型的表示。BREP格式文件以文本形式存儲，便于讀寫。

自動生成輸入卡模塊程序?qū)REP文件轉(zhuǎn)化為MCNP輸入幾何卡，然后添加材料卡，數(shù)據(jù)卡等信息，生成完整的MCNP輸入卡文件。對如圖2所示的幾何模型，其BREP文件部分內(nèi)容如圖3所示。

從圖3中可以看到，BREP文件中包含了實體模型的幾何信息和拓撲信息。幾何信息中位置信息部分以標識“Locations”開始，其后的數(shù)字表示位置信息的數(shù)量。每個位置信息是一個 3×4的矩陣，描述三維空間的線性變換?！癝urfaces”標識符下包含幾何模型所有的曲面信息，圖3中的曲面信息為一個圓柱面的表示方式，三維正交坐標系中圓柱面的軸通過點(-6,-6,0)，方向為[0,0,1]，半徑為0.5，其參數(shù)方程為：

S(u,v)=P+r·(cos(u)·Dx+sin(u)·Dy)
+v·Dv,(u,v)∈[0,2π)×(-∞,∞)

(1)

其中，Dv，Dx，Dy一起組成三維正交坐標系，圓柱面的中心軸通過點P，方向為Dv，圓柱面的半徑為r?！癟Shapes”標識符下包含幾何模型中各基本幾何元素間的拓撲關系信息，因此可以獲取幾何實體與各曲面之間的關系?！癋a”表示面(face)，“So”表示實體(solid)，圖3中的拓撲信息表示一個由圓柱面和上下底面構成的圓柱體。

MCNP輸入卡中的幾何描述包括曲面卡和柵元卡兩部分。曲面卡包含幾何模型的所有曲面信息，相當于BREP文件中的“Surfaces”部分，因此直接將曲面信息轉(zhuǎn)化為MCNP輸入卡中的柵元卡信息，MCNP中曲面的表示規(guī)則見表1， MCNP中使用一般方程對曲面進行描述，而BREP中對曲面的描述則使用參數(shù)方程，因此生成MCNP的曲面卡信息需要獲得兩者之間的轉(zhuǎn)換關系。由于球面、圓柱面等曲面與坐標軸之間存在關系，其一般方程參數(shù)值可以根據(jù)BREP文件中的曲面參數(shù)可以直接確定，本文主要處理平面的參數(shù)方程與一般方程之間的關系，根據(jù)向量共面的條件可知：

(2)

其中，(x0,y0,z0)為平面上的點，(X1,Y1,Z1)和(X2,Y2,Z2)分別為平面上的向量Du和Dv的坐標。求解行列式再與平面一般方程比較，獲得兩者之間的轉(zhuǎn)換關系(式3)。

(3)

表1 MCNP曲面描述

柵元卡中的柵元由曲面定義的半空間通過正則運算組合而成，模型中的每個區(qū)域都必須定義，不能存在空隙。因此MCNP柵元構建的算法步驟為：

1)根據(jù)BREP文件中提取的體與面拓撲信息生成組成幾何模型的實體集，若某個實體存在位置變換信息，則將當前位置與位置變換矩陣Q相乘以確定實體的最終位置；

2)確定實體集中各個實體之間的關系，本文中只考慮一個實體完全包含另一個實體和兩個實體分開這兩種情況，最終生成實體樹表示各個實體之間的關系；

3)最后從葉子節(jié)點開始遍歷實體樹上的所有實體節(jié)點，生成“a ±f1±f2… #b1#b2…”格式的柵元信息，其中a為柵元號，f1、f2為當前組成實體曲面邊界的曲面號，±表示曲面方向，b1、b2為實體節(jié)點的所有子實體節(jié)點的柵元號。

使用模型轉(zhuǎn)換算法將圖2所示的實體模型的BREP格式表示轉(zhuǎn)換為MCNP輸入卡的幾何描述格式，最終結果如圖4所示，符合MCNP輸入卡的格式要求。

圖4 BREP→MCNP模型轉(zhuǎn)換結果Fig.4 The result of converting BREP file into MCNP geometry input card

1.2 對MCNP程序計算核心的調(diào)用

對自動生成的幾何卡補充材料、數(shù)據(jù)等信息后生成完整的輸入文件，在集成平臺內(nèi)部調(diào)用MCNP計算核心就能夠直接進行中子輸運計算，并捕獲程序輸出的計算過程信息，同步顯示到當前平臺的信息區(qū)，運行界面如圖5所示。

圖5 MCNP計算信息圖Fig.5 The MCNP calculation information

1.3 計算結果后處理

MCNP的計算結果是如圖6所示的文本格式，包含計算結果數(shù)據(jù)以及一些特殊字符串[13]。當前對MCNP計算結果的處理絕大多數(shù)仍使用傳統(tǒng)的手工數(shù)據(jù)分析方法，需要人員從結果文件中提取數(shù)據(jù)，再導入專業(yè)的數(shù)據(jù)分析軟件中進行處理，工作量大且效率低[14]。集成平臺中的后處理模塊較好地實現(xiàn)了數(shù)據(jù)提取和繪制圖表兩個功能。

圖6 MCNP結果文件片段Fig.6 The result file fragment of MCNP

數(shù)據(jù)提取功能通過檢索結果文件中關鍵字定位數(shù)據(jù)所在位置來實現(xiàn)，如圖6中的“k(coll)”字符串，其后為具體的數(shù)據(jù)值。但從圖6中可以看到，數(shù)據(jù)值之間可能存在“|”等特殊字符，因此每獲取一行數(shù)據(jù)，需要使用正則表達式過濾數(shù)據(jù)中的特殊字符并將行數(shù)據(jù)分隔為一維數(shù)組，最終處理完所有數(shù)據(jù)后獲得一個二維數(shù)組，將二維數(shù)組保存到Excel格式的中間文件中[15]，完成數(shù)據(jù)提取工作。

繪制圖表功能利用Java繪圖工具包JFreeChart[16]將結果文件中提取的數(shù)據(jù)繪制為折線圖。折線圖繪制過程中，先獲取數(shù)據(jù)集對象DataSet，將繪圖數(shù)據(jù)添加到數(shù)據(jù)集中，然后調(diào)用JFreeChart的API生成并顯示折線圖表。

2 集成軟件平臺測試結果

完整的集成軟件平臺程序的界面如圖7所示，左側為程序中各模塊的操作區(qū)，各個模塊既可以一起使用，也可以單獨使用。右側為信息區(qū)，主要顯示各模塊中生成的各種文件信息、過程信息等。

圖7 程序界面Fig.7 The program interface

為測試軟件中各模塊的功能以及效率，對某小型堆堆芯的臨界系數(shù)Keff進行計算，以獲取最佳燃料富集度。所選取堆芯得具體參數(shù)如表2所示[17]，該堆型采用7×7的堆芯布局，包括37個燃料組件，每個燃料組件為標準的壓水堆燃料組件[18]，如圖8所示，其中B為可燃毒物棒，G為導向管，I為儀表管，余下的為UO2燃料棒。計算中使用富集度為2.4%的UO2燃料棒。

表2 小型堆堆芯幾何參數(shù)

圖8 燃料組件布局Fig.8 The arrangement of the fuel assembly

本計算中子源位于堆芯中心，源強為20 000，模擬300代中子循環(huán)下堆芯臨界系數(shù)Keff的變化，計算完成后提取結果文件中各中子代循環(huán)下的Keff值，并保存到Excel文件中，提取數(shù)據(jù)如圖9所示。然后使用繪圖功能生成圖10所示的結果曲線圖。當程序計算到100代循環(huán)時，臨界系數(shù)Keff已趨于穩(wěn)定，因此為節(jié)省時間，后續(xù)計算取100代中子循環(huán)的結果為有效結果進行分析。

圖9 臨界計算數(shù)據(jù)Fig.9 Results of criticality calculation

圖10 堆芯臨界計算結果Fig.10 Calculation results of core criticality

可以看出，本軟件處理數(shù)據(jù)較為方便，并且后續(xù)容易根據(jù)用戶需求開發(fā)更豐富的后處理功能。

從圖10可以看出，當前計算模型下，堆芯臨界系數(shù)在1.07上下波動，此時堆芯處于超臨界狀態(tài)。為了獲得臨界狀態(tài)，在相同幾何布局下，調(diào)整燃料富集度分別為1.0%、1.5%、2.0%，進行了多輪計算，以獲得堆芯達到臨界狀態(tài)時的燃料富集度。對比計算結果如圖11所示。

圖11 不同富集度下的臨界系數(shù)Keff比較Fig.11 Comparison of the critical coefficient Keff under different enrichment degrees

由圖11可知，隨著燃料富集度的提升，堆芯臨界系數(shù)Keff逐漸增大，該堆型的堆芯達到臨界狀態(tài)時的燃料富集度在1.5%～2.0%。在此富集度區(qū)間再進行多輪計算，最終得到該堆型堆芯臨界的最佳燃料富集度為1.85%，結果如圖12所示。

圖12 富集度1.85%的Keff變化曲線Fig.12 Variation of the Keff in every generation with the enrichment of 1.85%

整個測試過程顯示，采用集成軟件平臺，完成完整一輪堆型驗證比原有方式節(jié)約2人·天，總體效率提高約60%。

3 結束語

通過對BREP文件格式和MCNP程序輸入輸出文件的深入研究，采用Python和Java聯(lián)合編程技術，將幾何模型轉(zhuǎn)輸入卡的處理、中子輸運計算核心的調(diào)用、計算結果的圖形化后處理實現(xiàn)了一體化集成，搭建了中子輸運仿真集成軟件平臺，平臺中各個模塊既可以關聯(lián)使用，又能單獨使用，具有較高的靈活性。

通過堆芯臨界分析項目對集成軟件平臺的有效性進行了驗證。結果顯示，該軟件平臺在MCNP輸入文件準備中，解決了人工建模易出錯、耗費時間長的問題，在后處理中解決了數(shù)據(jù)處理繁瑣的問題，整個仿真流程的效率提高60%，并且降低了對工程分析人員的技術要求。這為滿足日益增加的反應堆物理隨機分析場景，提供了一個便捷有效的工具。