中子輸運(yùn)計(jì)算集成軟件平臺(tái)的研究與設(shè)計(jì)

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

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

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

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

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

1 集成軟件平臺(tái)的功能及關(guān)鍵模塊

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

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

1.1 BREP文件到MCNP幾何輸入卡的自動(dòng)生成

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

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

從圖3中可以看到，BREP文件中包含了實(shí)體模型的幾何信息和拓?fù)湫畔ⅰ缀涡畔⒅形恢眯畔⒉糠忠詷?biāo)識(shí)“Locations”開(kāi)始，其后的數(shù)字表示位置信息的數(shù)量。每個(gè)位置信息是一個(gè) 3×4的矩陣，描述三維空間的線性變換。“Surfaces”標(biāo)識(shí)符下包含幾何模型所有的曲面信息，圖3中的曲面信息為一個(gè)圓柱面的表示方式，三維正交坐標(biāo)系中圓柱面的軸通過(guò)點(diǎn)(-6,-6,0)，方向?yàn)閇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一起組成三維正交坐標(biāo)系，圓柱面的中心軸通過(guò)點(diǎn)P，方向?yàn)镈v，圓柱面的半徑為r。“TShapes”標(biāo)識(shí)符下包含幾何模型中各基本幾何元素間的拓?fù)潢P(guān)系信息，因此可以獲取幾何實(shí)體與各曲面之間的關(guān)系。“Fa”表示面(face)，“So”表示實(shí)體(solid)，圖3中的拓?fù)湫畔⒈硎疽粋€(gè)由圓柱面和上下底面構(gòu)成的圓柱體。

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

(2)

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

(3)

表1 MCNP曲面描述

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

1)根據(jù)BREP文件中提取的體與面拓?fù)湫畔⑸山M成幾何模型的實(shí)體集，若某個(gè)實(shí)體存在位置變換信息，則將當(dāng)前位置與位置變換矩陣Q相乘以確定實(shí)體的最終位置；

2)確定實(shí)體集中各個(gè)實(shí)體之間的關(guān)系，本文中只考慮一個(gè)實(shí)體完全包含另一個(gè)實(shí)體和兩個(gè)實(shí)體分開(kāi)這兩種情況，最終生成實(shí)體樹(shù)表示各個(gè)實(shí)體之間的關(guān)系；

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

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

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

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

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

圖5 MCNP計(jì)算信息圖Fig.5 The MCNP calculation information

1.3 計(jì)算結(jié)果后處理

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

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

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

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

2 集成軟件平臺(tái)測(cè)試結(jié)果

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

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

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

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

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

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

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

圖10 堆芯臨界計(jì)算結(jié)果Fig.10 Calculation results of core criticality

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

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

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

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

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

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

3 結(jié)束語(yǔ)

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

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