趙瑛峰，劉檢華，武林林，馬江濤

(北京理工大學 機械與車輛學院數(shù)字化制造研究所，北京 100081)

0 引言

MCNP(Monte Carlo N particle transport code)程序基于蒙特卡羅方法，用于計算三維復雜幾何結(jié)構(gòu)中的中子、光子、電子或耦合中子/光子/電子輸運問題[1]，被廣泛應用于輻射防護與射線測定、輻射屏蔽設計優(yōu)化、反應堆設計、醫(yī)學以及檢測器設計與分析等學科領(lǐng)域。由于其不提供三維建模工具，通常需要手工建立MCNP計算模型，費時費力。隨著計算機運算能力的提高，蒙特卡羅計算的幾何建模逐漸成為空間結(jié)構(gòu)復雜問題粒子輸運求解的瓶頸[2-5]。

國內(nèi)外解決MCNP建模難題的方法之一是將三維CAD模型轉(zhuǎn)換為MCNP計算模型。MCNP支持構(gòu)造實體幾何(Constructive Solid Geometry, CSG)模型描述，國內(nèi)外學者提出一種根據(jù)建模過程建立MCNP半空間CSG模型的方法，需要注意的是，在CAD軟件中采用這種方法時能夠使用的為MCNP程序規(guī)定的體素。Theis等[6-7]依據(jù)該思路開發(fā)了開源的蒙特卡洛幾何建模工具SimpleGeo；張建生等[8-10]采用UG NX二次開發(fā)完成了多個體素(立方體、球、圓柱等)的建模過程模型到MCNP計算模型的轉(zhuǎn)換；李春艷等[11]基于VC++和OpenGL實現(xiàn)了多個體素的MCNP三維建模；楊家明等[12-13]實現(xiàn)了多個簡單體素組合的以產(chǎn)品模型數(shù)據(jù)交互規(guī)范(STandard for the Exchange of Product model data, STEP)為格式的模型到MCNP模型的轉(zhuǎn)換；田彥山等[14]采用卷積神經(jīng)網(wǎng)絡對柵元位置判定方法進行了改進。上述方法需要在CAD軟件中完全采用簡單體素CSG建模，因此通用性和實用性不強。

為了增加通用性和實用性，國內(nèi)外學者提出一種基于分解的CAD模型到MCNP模型轉(zhuǎn)換方法。德國卡爾斯魯厄理工學院(Karlsruher Institut für Technologie, KIT)中子物理和反應堆技術(shù)研究所開發(fā)了MCCad[15-16]，日本原子能機構(gòu)開發(fā)了Geomit[17-18]。1991年，Shapiro[19-20]和Vossler[21]最早進行了BRep模型表示到半空間CSG模型的理論研究，從理論上證實了實體表面曲面階數(shù)不超過2的任意三維實體BRep模型到CSG模型轉(zhuǎn)換的可行性，并提出一種分割半空間的生成方法[22];Buchele等[23]在前者研究的基礎上提出B-rep模型向半空間CSG模型轉(zhuǎn)換(B-rep to Halfspace CSG conversion, BHC)的算法，將BRep模型轉(zhuǎn)換為BSP模型的算法和BSP模型轉(zhuǎn)換為CSG模型算法進行合并，減少了計算成本；中科院等離子體研究所熱核反應數(shù)字仿真(Fusion Digital Simulation, FDS)團隊開發(fā)了建模軟件系統(tǒng)MCAM(MCNP automatic modeling)[24]；羅月童[25-27]在Shapiro和Vossler等的基礎上提出基于分解的轉(zhuǎn)換算法，并將該算法集成到MCAM軟件中；周敬之等[28]提出一種基于特征識別的CAD模型空腔生成算法，提高了中空管道類模型的空腔構(gòu)造速度，并將其集成到MCAM中；馬露杰等[29]提出一種基于面殼封閉的方法來分解BRep模型；羅月童等[30]在馬露杰研究的基礎上提出基于面殼封閉的BRep模型轉(zhuǎn)換CSG模型算法，實驗結(jié)果證明該算法對部分模型有效；吳斌等[31]研究了一種基于輔助面的分解算法。以上研究只解決了部分模型轉(zhuǎn)換問題，對預定義特征中未定義的模型仍需要人工分解后再進行轉(zhuǎn)換；而且，輔助面自動生成算法方面的研究和算法比較少，在實際應用過程中仍然需要大量的人工操作，降低了模型的轉(zhuǎn)換效率。

為了減少模型轉(zhuǎn)換過程中的人工操作并提升算法效率，本文提出一種基于特征分解的模型BRep表達到半空間CSG表達轉(zhuǎn)換算法，通過拓撲和幾何特征識別及分解優(yōu)先級排序?qū)Rep表達的模型遞歸分解為半空間CSG樹，樹的葉節(jié)點是可以直接轉(zhuǎn)換為MCNP柵元卡的簡單幾何體。分解后，首先生成葉節(jié)點的面和柵元卡，再通過葉節(jié)點的布爾運算(并集和差集運算)完成整個轉(zhuǎn)換過程，最后在UG NX上完成算法開發(fā)并進行實例驗證。結(jié)果表明，該方法可以直接利用現(xiàn)有BRep模型快速生成滿足中子學計算的蒙特卡羅計算模型，顯著提高MCNP的建模效率。

1 算法框架

1.1 問題描述

BRep幾何實體Γ是三維歐氏空間R3的正則子集。BRep幾何實體Γ的邊界表示為bd(Γ)=op·fi(fi∈F,op·∈{∪*,∩*,-*})，其中fi為Γ中的自然表面，F(xiàn)為Γ中所有自然表面fi的集合。

函數(shù)g定義的面Sψi={(x,y,z):g(x,y,z)=0}將三維歐氏空間R3分割為兩個半空間，即負半空間ψ1={(x,y,z):g(x,y,z)≤0}和正半空間ψ2={(x,y,z):g(x,y,z)≥0}。

BRep幾何實體Γ到半空間CSG轉(zhuǎn)換的本質(zhì)是尋找一個半空間集合Ψ={ψ1,…,ψn}，使得bd(Γ)=bd(ψ1)op·…op·bd(ψn)(其中ψi為三維歐氏空間R3的非空、正則順向連通子集，ψi之間采用正則化幾何操作，包括并集∪*，以及差集-*和交集∩*)，即將BRep幾何實體Γ等價分解為一組簡單實體組合，而簡單實體直接對應一個半空間實體ψi。

MCNP蒙特卡羅輸運中的幾何計算主要是在源采樣過程中判斷粒子所在位置，并在模擬粒子隨機游動過程中追蹤粒子的當前位置。幾何計算過程包括計算粒子沿當前方向到達柵元邊界的距離和搜索粒子穿出邊界后的下一個柵元，計算的基礎是直線與面求交，以及點與半空間位置關(guān)系的判斷。一方面，粒子關(guān)于一個柵元計算的復雜度與組成柵元的面的數(shù)量呈正比；另一方面，粒子穿出柵元后的位置與柵元數(shù)量呈正比。因此，為提高MCNP計算效率，BRep幾何實體Γ分解后的柵元數(shù)量要盡量少，而且每一個柵元中面的數(shù)量也要盡量少。

1.2 定義

定義1自然表面fi。fi為CAD幾何實體Γ的一個邊界面，范圍是實體Γ邊界表示中的外環(huán)和內(nèi)環(huán)之間的區(qū)域。

定義2面Syi。Sψi={(x,y,z):g(x,y,z)=0}唯一定義面，自然表面fi唯一對應面Syi，且fi?Sψi。

定義3表面集。指BRep幾何實體Γ的所有自然表面fi對應的面Syi的集合，用Surface(Γ)表示。

定義4CSG樹。CSG樹是單個根節(jié)點的二叉樹，CSG樹的子葉節(jié)點為轉(zhuǎn)換元實體，內(nèi)部節(jié)點之間采用正則化的幾何操作，包括正則并(∪*)、正則交(∩*)、正則差(-*)3種運算。

定義5轉(zhuǎn)換元實體。BRep幾何實體Γ中的每一點均在其表面集Surface(Γ)中面Syi的正或負半空間，此處的轉(zhuǎn)換元實體對應MCNP程序中的柵元。

定義6完全轉(zhuǎn)換元實體。指可被其表面集Surface(Γ)的半空間完整描述的轉(zhuǎn)換元實體。

定義7不完全轉(zhuǎn)換元實體。指不能被其表面集Surface(Γ)的半空間完整描述的轉(zhuǎn)換元實體，例如包含圓錐面的轉(zhuǎn)換元實體均為不完全轉(zhuǎn)換元實體。

定義8半空間輔助面。對于不完全轉(zhuǎn)換元實體Γ，如果存在一個或多個面Syi，通過Syi的半空間切分使?？梢员黄浔砻婕疭urface(Γ)的半空間完整描述，則面Syi是Γ的半空間輔助面。

定義9充分表面集。完全轉(zhuǎn)換元實體Γ的表面集、不完全轉(zhuǎn)換元實體Γ的表面集及其半空間輔助面成為對應實體Γ的充分表面集。

1.3 算法流程

幾何實體的識別和分解以及輔助面的添加是算法的關(guān)鍵，算法總體技術(shù)流程如圖1所示。

步驟1根據(jù)模型的拓撲結(jié)構(gòu)和幾何特點識別模型特征。

步驟2根據(jù)特征類型分解模型。

步驟3生成特征的轉(zhuǎn)換元實體。

步驟4判斷轉(zhuǎn)換元實體是否為不完全轉(zhuǎn)換元實體，是則添加輔助面，使當前實體轉(zhuǎn)換為完全轉(zhuǎn)換元，否則執(zhí)行下一步。

步驟5判斷是否存在未識別的特征，是則返回步驟2，否則算法結(jié)束。

2 特征識別和分解策略

由于分解策略不同，幾何實體存在不同的分解結(jié)果。較優(yōu)的分解策略能夠使用較少的柵元和表面正確地描述幾何實體，從而提高MCNP計算效率。一個幾何實體的兩種分解結(jié)果如圖2所示。為了獲得更好的分解結(jié)果，本文采用基于幾何和拓撲特征識別的分解策略，即根據(jù)BRep幾何實體中體—殼—面—環(huán)—邊的拓撲關(guān)系將幾何特征分類，并采用相應的分解方式進行分解。本文提出一套特征識別規(guī)則和相應的分解策略。

本文將實體特征分為內(nèi)殼特征、內(nèi)環(huán)凸特征、內(nèi)環(huán)非貫穿凹特征、內(nèi)環(huán)貫穿凹特征和外環(huán)凸特征5個類型，根據(jù)這5種特征提出不同的分解方式。

(1)內(nèi)殼特征的判斷依據(jù)為幾何實體內(nèi)存在內(nèi)殼，典型內(nèi)殼特征的例子如圖3a所示。內(nèi)殼特征的分解方式為根據(jù)內(nèi)殼表面構(gòu)成的封閉空間生成填充幾何體S2，與填充后的體S1進行差運算，典型內(nèi)殼特征的分解結(jié)果如圖3a所示。

(2)內(nèi)環(huán)凸特征的判斷依據(jù)為表面上存在內(nèi)環(huán)，且內(nèi)環(huán)連接的表面位于表面連接點外法向同側(cè)，典型內(nèi)環(huán)凸特征的例子如圖3b所示。內(nèi)環(huán)凸特征的分解方式為根據(jù)表面的內(nèi)環(huán)邊界構(gòu)建面片F(xiàn)，通過F將幾何實體分開，典型內(nèi)環(huán)凸特征的分解結(jié)果如圖3b所示。

(3)內(nèi)環(huán)非貫穿凹特征的判斷依據(jù)為表面上存在內(nèi)環(huán)，且內(nèi)環(huán)連接的表面位于表面連接點外法向反向，內(nèi)環(huán)形成的區(qū)域沒有將幾何體形成通孔，典型內(nèi)環(huán)非貫穿凹特征的例子如圖3c所示。內(nèi)環(huán)非貫穿凹特征的分解方式為根據(jù)內(nèi)環(huán)邊界曲線連接的表面連通圖構(gòu)建封閉區(qū)域S2，與填充后的體S1進行差運算。典型內(nèi)環(huán)非貫穿凹特征的分解結(jié)果如圖3c所示。

(4)內(nèi)環(huán)貫穿凹特征的判斷依據(jù)為多張表面上存在內(nèi)環(huán)，內(nèi)環(huán)連接的表面位于表面連接點外法向反向，且內(nèi)環(huán)之間互相連通，內(nèi)環(huán)形成的區(qū)域在幾何體形成通孔，典型內(nèi)環(huán)貫穿凹特征的例子如圖3d所示。內(nèi)環(huán)貫穿凹特征的分解方式為經(jīng)連通的多個內(nèi)環(huán)區(qū)域構(gòu)建封閉區(qū)域S2，與填充后的體S1進行差運算。典型內(nèi)環(huán)貫穿凹特征的分解結(jié)果如圖3d所示。

(5)外環(huán)凸特征的判斷依據(jù)為表面f1外環(huán)所在邊與表面f2相連，f2在f1外側(cè)，典型外環(huán)凸特征的例子如圖3e所示。外環(huán)凸特征的分解方式為根據(jù)表面f1構(gòu)建面片F(xiàn)，通過F將幾何實體分開。典型外環(huán)凸特征的分解結(jié)果如圖3e所示。

將以上特征進行優(yōu)先級排序，優(yōu)先級由高到低的排序為內(nèi)殼特征、內(nèi)環(huán)凸特征、內(nèi)環(huán)非貫穿凹特征、內(nèi)環(huán)貫穿凹特征、外環(huán)凸特征。在分解過程中，按照優(yōu)先級由高到低的順序識別特征并分解特征。幾何實體特征分級的偽代碼如下：

1 Algorithm:Feature recognition & decomposition

2 Input:Γcurrent，{Γ}

3 if Γcurrentis null，return{Γ};

6 Γrest;

8 part Γrest;

10 Γrest;

12 Γrest;

14 return Algorithm(Γrest, {Γ})

3 半空間輔助面

為了避免轉(zhuǎn)換后的不完全轉(zhuǎn)換元實體存在歧義，需要在轉(zhuǎn)換完成后添加半空間輔助面。如圖4所示，幾何體S分解為S1和S2的組合，然而S2轉(zhuǎn)換為MCNP半空間描述后的結(jié)果相當于圓柱面的負半空間和球的正半空間的交集，即在一個無限長的圓柱中間挖去一個球，形成一個錯誤的幾何描述。因此，必須構(gòu)造S2的半空間輔助面，形成S2充分表面集，即在圓柱和球的相交處通過球的圓心構(gòu)建半空間輔助面F才能正確描述S2的半空間。

半空間輔助面的添加可以依據(jù)MCNP相交面類型進行相應處理，如表1所示。MCNP支持的面類型包括平面、球面、圓柱面、圓錐面、橢球面、拋物面、雙曲面和圓環(huán)面，表中沒有考慮不常用的拋物面和雙曲面。

表1 半空間輔助面的添加方式

典型的不完全半空間轉(zhuǎn)換實體及其對應添加的半空間輔助面如圖5所示。

在部分幾何實體中，輔助面還可以作為分界面，如圖6所示，幾何實體S的圓環(huán)面分別與平面和圓柱面相切。在識別該幾何特征后，首先在圓環(huán)和圓柱面相切處創(chuàng)建半空間輔助面F1，通過F1將S分解為S1和S2；然后在圓環(huán)和平面相切處構(gòu)建半空間輔助面F2，將S2分解為S3和S4；最后得到分解后的CSG樹。

4 實驗與應用

本文在NX平臺上設計并開發(fā)了MCNP三維前處理軟件系統(tǒng)，同時在MCNP三維前處理軟件系統(tǒng)與MCNP計算程序之間定義了中間數(shù)據(jù)格式。這樣的中間數(shù)據(jù)格式以CSG樹的二叉樹數(shù)據(jù)結(jié)構(gòu)為基礎，二叉樹的根節(jié)點為當前轉(zhuǎn)換的零件，葉節(jié)點是正則計算符與轉(zhuǎn)換元實體，轉(zhuǎn)換元實體的存儲格式為NX平臺提供的模型格式，而且在NX平臺中可查詢轉(zhuǎn)換元實體內(nèi)與MCNP半空間體素對應的幾何表面。在本文系統(tǒng)中，利用典型實例CAD模型進行轉(zhuǎn)換的過程與結(jié)果如圖7和圖8所示，結(jié)果表明本文算法可以獲得較優(yōu)的分解結(jié)果，而且圖8中展示了輔助面的添加結(jié)果。

為了進一步驗證本文方法的有效性，本文建立了輕水增殖反應堆的三維模型，并猜用本文系統(tǒng)進行了預處理。輕水增殖反應堆項目是被授權(quán)開發(fā)的一個商業(yè)化增殖反應堆，使用的是當代的(20世紀60年代～70年代)輕水反應堆技術(shù)[32]。輕水增殖反應堆的結(jié)構(gòu)如圖9a所示，其每個堆芯都位于一個支撐結(jié)構(gòu)上，而且反應堆罐被安全殼覆蓋，在反應堆罐上方有一個鉛罩；有一個圓柱形不銹鋼反應堆池，壁厚為1/4英寸，內(nèi)徑大小為6英尺；圓柱形反應堆池的外部用1/2英寸厚的鉛板包圍起來形成一個八邊形；在反應堆操作期間，反應堆池填滿水；反應堆池位于直徑為17英尺圓柱形安全殼的一個偏離中心的5英尺甲板上，安全殼也用慢化劑儲存。圖9b所示為輕水增殖反應堆三維模型分解后的結(jié)果。通過本文方法，三維模型被分解為轉(zhuǎn)換元實體，圖9b的高亮部分表示輕水增殖反應堆的一個轉(zhuǎn)換元實體。輕水增殖反應堆的部分MCNP卡片如圖10所示。將本文系統(tǒng)生成的輕水增殖反應堆MCNP卡片應用于MCNP計算程序，計算結(jié)果為1.000 9±0.000 6，國際基準題的推薦計算結(jié)果為0.997 7±0.000 9，兩個結(jié)果相近，驗證了本文方法的有效性。

為了驗證本文方法的效率，在CPU參數(shù)為3.20 GHz Intel Core i5-3470、內(nèi)存參數(shù)為4 G的硬件環(huán)境中進行了計算實驗。本文以模型的面數(shù)量為橫軸，以分解時間為縱軸，建立相關(guān)性曲線，如圖11所示。根據(jù)圖11的曲線，模型的面數(shù)量與分級時間之間呈正相關(guān)關(guān)系，且曲線的增長速率較緩，證明本文算法效率較高。

5 結(jié)束語

為解決蒙特卡羅粒子輸運計算的建模難題，本文提出一種基于特征分解的BRep模型到半空間CSG模型的轉(zhuǎn)換算法。算法通過對幾何特征識別及相應的分解方法，按照優(yōu)先級順序?qū)缀螌嶓w進行分解，獲得半空間CSG模型，并通過添加輔助面的方式消除非完全轉(zhuǎn)換元的描述歧義。本文算法基于NX系統(tǒng)開發(fā)完成，并通過實例驗證了算法的有效性和可靠性。未來將進行更多幾何特征分類和分解方法的研究，以提高轉(zhuǎn)換的效果和效率。