基于CAD技術(shù)的特征線中子輸運(yùn)計(jì)算程序開發(fā)

陳珍平，王電喜，何 桃，王國(guó)忠，鄭華慶，F(xiàn)DS團(tuán)隊(duì)

（1.中國(guó)科學(xué)技術(shù)大學(xué)，安徽 合肥 230027；2.中國(guó)科學(xué)院核能安全技術(shù)研究所，安徽 合肥230031）

中子輸運(yùn)問(wèn)題一直是核反應(yīng)堆物理分析與輻射屏蔽設(shè)計(jì)研究的重點(diǎn)問(wèn)題之一。隨著反應(yīng)堆技術(shù)的快速發(fā)展，在未來(lái)進(jìn)行先進(jìn)反應(yīng)堆設(shè)計(jì)與分析時(shí)對(duì)中子學(xué)模擬計(jì)算程序的要求也越來(lái)越高，主要體現(xiàn)在計(jì)算精度、計(jì)算效率以及處理幾何復(fù)雜度等方面，因此開發(fā)具有自主知識(shí)產(chǎn)權(quán)的、快速精確的中子學(xué)計(jì)算程序?qū)⒕哂蟹浅Ｖ匾囊饬x。特征線方法[1］（Method of Characteristics，MOC）是近年來(lái)在中子輸運(yùn)計(jì)算與數(shù)值模擬研究中熱點(diǎn)研究的方法之一。該方法是從中子輸運(yùn)方程的全微分形式出發(fā)，沿著中子飛行軌跡（稱為特征線）進(jìn)行積分求解的一種確定論方法。由于該方法同時(shí)具有碰撞概率法和離散縱標(biāo)法的優(yōu)點(diǎn)，且理論上不受計(jì)算模型幾何形狀的限制，近十多年來(lái)已逐漸成為國(guó)際上求解中子輸運(yùn)方程的主要確定論方法之一。

長(zhǎng)期以來(lái)，阻礙特征線方法能否廣泛應(yīng)用的主要問(wèn)題在于能否將特征線理論與有效的幾何處理方法結(jié)合起來(lái)。本文在FDS團(tuán)隊(duì)自主研發(fā)的大型集成多功能中子學(xué)計(jì)算與分析系統(tǒng)VisualBUS[2－3］框架下，利用核與輻射輸運(yùn)計(jì)算自動(dòng)建模軟件 MCAM[4－6］的幾何處理能力進(jìn)行特征線方法幾何預(yù)處理，可高效、快速地完成幾何建模、網(wǎng)格剖分和特征線生成等幾何處理功能。最后，利用基準(zhǔn)例題對(duì)特征線數(shù)值計(jì)算方法和程序進(jìn)行了充分地檢驗(yàn)，其結(jié)果驗(yàn)證了本研究幾何處理方法和計(jì)算程序的正確性與可靠性。

1 特征線理論

特征線方法是將中子輸運(yùn)方程沿著中子飛行軌跡和特定的離散方向進(jìn)行輸運(yùn)求解的一種方法。在本文中，對(duì)方位角離散采用Filippone[7］提出的離散方法；對(duì)極角離散耦合了LO離散和TY離散[8］方法；對(duì)幾何模型區(qū)域進(jìn)行非結(jié)構(gòu)網(wǎng)格剖分，劃分為許多平源近似區(qū)，在每個(gè)平源區(qū)里的中子源項(xiàng)被認(rèn)為是常數(shù)。

圖1 特征線示意圖Fig.1 The characteristic lines of MOC

根據(jù)前面的基本假設(shè)，單群特征線中子輸運(yùn)方程可表示為：

其中：

s是特征線長(zhǎng)度，Σt，i是區(qū)域i的宏觀總截面，Ψi，k（s，Ωm）是區(qū)域i沿著特征線方向長(zhǎng)度為s處的中子角通量密度，Qf（Ωm）是區(qū)域i的中子源項(xiàng)。

方程（1）為一階線性微分方程，可以解析求解得到：

式（2）中Ψini，k（Ωm）是特征線k入射到區(qū)域i里入射點(diǎn)處的入射角通量密度。

根據(jù)方程（2）可以得到，特征線k穿出區(qū)域i的出射角通量密度為：

式（3）中si，k是特征線k穿過(guò)區(qū)域i的總長(zhǎng)度；將特征線k上所有點(diǎn)的角通量密度進(jìn)行積分平均，則可以獲得特征線k上任一點(diǎn)的平均角通量密度：

同理，將穿過(guò)區(qū)域i的所有特征線上的平均通量密度在i區(qū)進(jìn)行體積加權(quán)平均，可獲得區(qū)域i的平均角通量密度：

式中δAk為特征線間距。

將區(qū)域i各個(gè)離散方向的角通量密度進(jìn)行加權(quán)求和，得到區(qū)域i的平均中子標(biāo)通量密度：

式中ωm是方向Ωm對(duì)應(yīng)的離散權(quán)重，M為空間離散方向總數(shù)目。

2 幾何預(yù)處理方法

2.1 幾何建模

FDS團(tuán)隊(duì)自主研發(fā)的核與輻射輸運(yùn)計(jì)算自動(dòng)建模軟件MCAM實(shí)現(xiàn)了多種商用工程CAD軟件（CATIA、AutoCAD、UG等）與多種輻射輸運(yùn)計(jì)算程序（MCNP、TRIPOLI[9］等）之間的接口。一方面，可以直接由計(jì)算機(jī)輔助設(shè)計(jì)（CAD）模型生成完整的輻射輸運(yùn)計(jì)算程序的輸入模型，包括空腔模型、材料、源和計(jì)數(shù)信息等；另一方面，可以解析輻射輸運(yùn)計(jì)算程序的輸入模型，生成CAD模型并可視化，供分析和修正。MCAM具有非常強(qiáng)大的幾何建模能力，能夠方便地建立各種復(fù)雜幾何的CAD模型。目前MCAM已經(jīng)在一些復(fù)雜核裝置（如FDS系列聚變反應(yīng)堆[10］、ITER 模型[11－13］等）的中子學(xué)建模和計(jì)算分析中得到了廣泛的應(yīng)用。

本文基于MCAM提供的造型功能和內(nèi)嵌的布爾運(yùn)算功能，在MCAM平臺(tái)上進(jìn)行二次開發(fā)了MOC幾何建模模塊，通過(guò)調(diào)用該模塊的接口，可以快速、精確地進(jìn)行幾何建模，且能方便地構(gòu)造出如反應(yīng)堆燃料棒、燃料組件和反應(yīng)堆堆芯等各種常見的幾何模型。同時(shí)，利用MOC建模模塊所構(gòu)造的幾何模型可以通過(guò)調(diào)用MCAM底層的渲染接口在MCAM的GUI中渲染出來(lái)（如圖2），可非常清晰、直觀地表達(dá)出模型的幾何結(jié)構(gòu)，便于用戶進(jìn)行模型正確性檢查。同時(shí)，還為幾何模型添加了材料屬性編輯功能，即用戶可以對(duì)圖2中的幾何模型進(jìn)行材料屬性添加，則MOC程序在進(jìn)行輸運(yùn)計(jì)算時(shí)，可以直接根據(jù)幾何模型的材料屬性來(lái)獲取材料的各種截面信息。最終構(gòu)造好的幾何模型（已添加材料屬性）可保存為標(biāo)準(zhǔn)的CAD文件格式（sat，step，igs等），MOC 程序做中子輸運(yùn)計(jì)算時(shí)可以直接將上述文件導(dǎo)入進(jìn)行輸運(yùn)計(jì)算，即本文開發(fā)的MOC程序支持CAD文件直接導(dǎo)入進(jìn)行輸運(yùn)計(jì)算的功能。

圖2 MCAM幾何建模Fig.2 Geometry modeling of MCAM

2.2 特征線跟蹤

由前面的特征線公式（1）～（5）可知，在進(jìn)行特征線中子輸運(yùn)計(jì)算時(shí)，需要不斷地進(jìn)行特征線跟蹤來(lái)獲取模型的幾何和材料等相關(guān)信息。所以，MOC程序在進(jìn)行中子輸運(yùn)計(jì)算時(shí)所需要的幾何和材料等信息主要以特征線為載體，主要包括：特征線間距、特征線段長(zhǎng)度、特征線所在網(wǎng)格編號(hào)、對(duì)應(yīng)的材料號(hào)以及邊界信息等。本文主要基于MCAM提供的射線跟蹤功能，二次開發(fā)了MOC程序特征線自動(dòng)處理模塊，圖3表示了整個(gè)MOC程序中幾何建模模塊、射線跟蹤模塊以及物理計(jì)算模塊等各模塊間的相互關(guān)系。

圖3 MOC程序模塊關(guān)系圖Fig.3 The relationship of different modules in MOC

在開始特征線跟蹤之前，程序根據(jù)用戶在外部配置文件中提供的方位角數(shù)目，自動(dòng)進(jìn)行方位角離散并計(jì)算每個(gè)方位角對(duì)應(yīng)的特征線間距；然后根據(jù)網(wǎng)格劃分參數(shù)對(duì)模型進(jìn)行網(wǎng)格劃分；最后根據(jù)已離散的方位角數(shù)值，通過(guò)調(diào)用MCAM的射線跟蹤接口來(lái)產(chǎn)生所有方位角下對(duì)應(yīng)的特征線信息，并將之保存用于MOC程序的輸運(yùn)計(jì)算。由于MCAM提供了射線跟蹤功能，并且計(jì)算精度和效率能保持在合理范圍內(nèi)，在二次開發(fā)時(shí)可直接繼承這些功能進(jìn)行應(yīng)用程序的開發(fā)，大大提高了應(yīng)用程序的開發(fā)效率。

3 數(shù)值驗(yàn)證

本文提出了基于MCAM進(jìn)行MOC幾何預(yù)處理的方法，并開發(fā)了特征線中子輸運(yùn)計(jì)算MOC程序，為了驗(yàn)證程序的正確性與可靠性，本文給出了一維7群例題和二維C5G7－UO2組件基準(zhǔn)例題的數(shù)值驗(yàn)證結(jié)果。

3.1 一維7群例題

本文采用文獻(xiàn)[14］中定義的一維7群?jiǎn)栴}，其幾何尺寸和材料布置如圖4所示，各種材料的7群截面來(lái)自C5G7－MOX基準(zhǔn)題[15］。由于該問(wèn)題材料間能譜差異大、泄漏較強(qiáng)，因此可以很好的驗(yàn)證MOC程序的正確性。

圖4 一維7群例題幾何布置（cm）Fig.4 The layout of 1D7Gproblem（cm）

該問(wèn)題的參考解由文獻(xiàn)[14］中的程序PEACH給出。表2給出了與PEACH程序計(jì)算結(jié)果的比較，其中本文采用平源菱形差分特征線法給出，平源區(qū)長(zhǎng)度為0.1cm，方位角數(shù)目為4個(gè)，特征線間距約為0.05cm，從表1中可以看出本文計(jì)算的結(jié)果與參考值相比，keff誤差約為－3pcm，結(jié)果符合良好。

表1 一維7群例題MOC與PEACH結(jié)果比較Table 1 Comparision of kefffor 1D7G problem between MOC and PEACH

表2 二維C5G7－UO2組件例題kinf計(jì)算結(jié)果比較Table 2 Comparision of kinffor C5G7－UO2assembly results

在特征線方法中，影響計(jì)算結(jié)果精度的主要因素包括劃分網(wǎng)格尺寸大小、方位角數(shù)目、極角數(shù)目、特征線間距等。針對(duì)本例題，本文就前兩個(gè)因素做了敏感性分析。圖5給出了在不同網(wǎng)格劃分、相同角度離散（TY－2極角和4個(gè)方位角）、相同特征線間距（0.05cm）情況下keff誤差敏感性分析。從圖5可看出，網(wǎng)格大小對(duì)keff精度會(huì)有較大影響；當(dāng)網(wǎng)格的尺寸小于1.0cm時(shí)，MOC程序能夠獲得非常滿意的計(jì)算精度（小于70pcm）。圖6給出了在相同網(wǎng)格尺寸（0.1cm）、相同極角離散（TY－2極角）、相同特征線間距（0.05cm）和不同方位角數(shù)目下keff誤差敏感性分析。從圖6中可以看出，當(dāng)方位角數(shù)目N≥4時(shí)，繼續(xù)增加方位角數(shù)目對(duì)最終結(jié)果影響不大，誤差絕對(duì)值小于3pcm，故采用4個(gè)方位角已經(jīng)可以獲得比較滿意的計(jì)算精度。

圖5 不同網(wǎng)格大小keff誤差分析Fig.5 kefferror analysis with different mesh dimensions

圖6 不同方位角數(shù)目keff誤差分析Fig.6 kefferror analysis with different angle dimensions

3.2 二維C5G7－UO2組件基準(zhǔn)例題

二維C5G7－UO2基準(zhǔn)題是國(guó)際經(jīng)濟(jì)合作與發(fā)展組織（OECD）發(fā)布的二維C5G7基準(zhǔn)例題[15］的二氧化鈾燃料組件，組件燃料布置與幾何結(jié)構(gòu)如圖7所示。上圖為17×17個(gè)柵元組成的UO2組件，包括UO2燃料柵元、導(dǎo)向管柵元和裂變室柵元。柵元為兩區(qū)非均勻結(jié)構(gòu)，各柵元的間距為1.26cm，柵元中的燃料芯塊和導(dǎo)向管均勻化后的半徑為0.54cm，組件四周采用全反射邊界條件。同時(shí)，采用圖7右下角的四分之一UO2燃料組件模型[16］進(jìn)行進(jìn)一步計(jì)算，通過(guò)和單組件計(jì)算結(jié)果對(duì)比，以進(jìn)一步驗(yàn)證程序的正確性。

圖7 2DC5G7－UO2燃料組件幾何布置（cm）Fig.7 The layout of the 2DC5G7benchmark with UO2assembly （cm）

本文MOC程序計(jì)算時(shí)各材料區(qū)的截面信息采用文獻(xiàn)[15］中二維C5G7模型提供的宏觀截面數(shù)據(jù)。本文對(duì)1／4UO2組件和1／1UO2組件分別進(jìn)行了kinf計(jì)算。在計(jì)算時(shí)，將模型進(jìn)行非結(jié)構(gòu)網(wǎng)格剖分；對(duì)空間角離散采用LO－2優(yōu)化兩極角和6方位角離散；特征線間距為0.02cm，計(jì)算結(jié)果如表2所示。其計(jì)算結(jié)果與 程 序 GALAXY[16］、GMVP[17］的 誤 差 約 為－0.15%，說(shuō)明與其他兩程序計(jì)算結(jié)果吻合良好；同時(shí)，MOC程序計(jì)算的1／4UO2組件和1／1UO2組 件 的kinf分 別 為 1.331 345 和1.331 782，相對(duì)誤差約為0.03%，兩模型計(jì)算結(jié)果符合良好，證明了本文MOC程序的正確性與可靠性。

總結(jié)

幾何預(yù)處理方法，即基于MCAM的幾何處理引擎，通過(guò)C＋＋語(yǔ)言面向?qū)ο蠖伍_發(fā)了MOC幾何預(yù)處理模塊，可以非常方便地進(jìn)行特征線幾何建模以及特征線信息快速獲取，這樣便將特征線理論與MCAM有效地結(jié)合起來(lái)。本文基于MCAM二次開發(fā)了特征線中子輸運(yùn)計(jì)算MOC程序，并利用相關(guān)基準(zhǔn)例題對(duì)程序進(jìn)行了驗(yàn)證，計(jì)算結(jié)果與國(guó)內(nèi)外其他相關(guān)計(jì)算程序的結(jié)果吻合良好，表明本文方法和程序的可行性、正確性與可靠性。

致謝

本文在開展研究過(guò)程中，得到了西安交通大學(xué)劉宙宇博士在特征線理論方面的幫助與指導(dǎo)，在此深表感謝。

本文提出了一種基于CAD技術(shù)的特征線

[1]Askew R.A Characteristics for Mulation of the Neutron Transport Equation in Complicated Geometries[R］.Report AEEW－R－1108，UK Atomic Energy Authority，1972.

[2]吳宜燦，李靜驚，李瑩，等.大型集成多功能中子學(xué)計(jì)算與分析系統(tǒng)VisualBUS的研究與發(fā)展[J］.核科學(xué)與工程，2007，27（5）：72－85.

[3]吳宜燦，胡麗琴，龍鵬程，等.先進(jìn)核能系統(tǒng)設(shè)計(jì)分析軟件與數(shù)據(jù)庫(kù)研發(fā)進(jìn)展[J］.核科學(xué)與工程，2010，30（1）：55－64.

[4]Wu Y C，F(xiàn)DS Team.CAD－based Interface Prog－rams for Fusion Neutron Transport Simulation [J］.Fusion Engineering and Design，2009，84：1987－1992.

[5]吳宜燦，李瑩，盧磊，等.蒙特卡羅粒子輸運(yùn)計(jì)算自動(dòng)建模程序系統(tǒng)的研究與發(fā)展[J］.核科學(xué)與工程，2006，26（01）：20－27.

[6]李瑩，曾勤，盧磊.利用ITER基準(zhǔn)模型對(duì) MCAM4.2進(jìn)行檢驗(yàn)[J］.核科學(xué)與工程，2008，28（1）：47－50.

[7]Filippone W L， Woolf S，Lavigne R J.Particle Transport Calculation with the Method of Streaming Rays[J］.Nuclear Science and Engineering，1981，77：119.

[8]Yamamoto A，Tabuchi M，Sugimura N，et al.Derivation of Optimum Polar Angle Quadrature Set for the Method of Characteristics Based on Approximation Error for the Bickley Function [J］. Journal of Nuclear and Technology，2007，44（2）：129－136.

[9]張俊軍，曾勤，王國(guó)忠，等.蒙特卡羅程序TRIPOLI自動(dòng)建模方法研究[J］.核科學(xué)與工程，2010，20（3）：84－88.

[10]Wu Y，F(xiàn)DS Team.Conceptual Design Activities of FDS Series Fusion Power Plants in China[J］.Fusion Engineering and Design，2006，81（23－24）：2713－2718.

[11]曾勤，盧磊，李瑩，等.蒙特卡羅粒子輸運(yùn)計(jì)算自動(dòng)建模程序MCAM在ITER核分析建模中的應(yīng)用[J］.原子核物理評(píng)論，2006，23（2）：138－141.

[12]王國(guó)忠，黨同強(qiáng)，熊健，等.MCAM4.8在ITER建筑大廳中子學(xué)建模中的應(yīng)用[J］.核科學(xué)與工程，2011，31（4）：351－355.

[13]熊健，王國(guó)忠，王電喜，等.MCAM 在ITER裝置TRIPOLI三維中子學(xué)建模中的應(yīng)用[J］.核科學(xué)與工程，2011，31（2）：162－167.

[14]湯春桃.中子輸運(yùn)方程特征線解法及嵌入式組件均勻化方法的研究 [D］.上海：上海交通大學(xué)，2009.

[15]Lewis E E，Simith M A，Tsoulfanidis N，et al.Benchmark Specification for Deterministic 2－D／3－D MOX Fuel Assembly Transport Calculations without Spatial Homogenization （C5G7MOX）[R］.NEA／NSC／DOC，2001.

[16]Yamaji K，Matsumoto H，Sato D，et al.Simple and Efficient Parallelization Method for MOC Calculation[J］.Journal of Nuclear Science and Technology，2010，47（1）：90－102.

[17]Mori T，Nakagawa T.MVP／GMVP：General Purpose Monte Carlo codes for Neutron and Photon Transport Calculations Based on Continuous Energy and Multigroup Methods [R］.Japan Atomic Energy Research Institute，1994.