特征線幾何預處理方法比較

李雪輝 劉紫靜 謝金森 于 濤 陳珍平 謝 芹 曾文杰 何麗華

特征線幾何預處理方法比較

李雪輝 劉紫靜 謝金森 于 濤 陳珍平 謝 芹 曾文杰 何麗華

（南華大學 核科學技術學院 衡陽 421001）

網格劃分、特征線間距、角度求積組、極角數目和方位角大小等幾何預處理過程對特征線法的計算精度和計算效率有較大影響?；诓教卣骶€法開發(fā)輸運程序，通過數值計算驗證所開發(fā)程序的正確性并分析兩種特征線掃描方法（首尾相間循環(huán)掃描法、首尾相接循環(huán)掃描法）以及網格劃分、特征線間距、角度求積組、極角數目、方位角大小對計算精度的影響。結果表明，開發(fā)的程序準確可靠；首尾相間循環(huán)掃描方法的收斂速度比首尾相接循環(huán)掃描方法慢。

特征線方法，幾何預處理，循環(huán)特征線掃描

特征線法(Method Of Characteristics, MOC)理論上不受幾何形狀限制，無需做均勻化處理，能夠精確處理強各向異性問題，已經成為反應堆物理計算的研究熱點[1-2]。特征線法能夠精確求解不同幾何模型的關鍵在于模型的不同平源近似區(qū)有“適當”數量的特征線穿過。平源近似區(qū)內沒有特征線穿過或數量太少，會降低計算精度，反之數量太多，會降低計算速度并增加內存消耗。因此，采用適當的幾何預處理方法產生并追蹤特征線至關重要。

目前國內外常見的幾何預處理方法有兩種：實數幾何法(R-Function)[3-4]和模塊化幾何法(Modular Ray Tracing, MRT)[5-6]。實數幾何法具有很強的幾何適應性，AGEIS和ANEMONA程序采用了該方法[5,7]。但該方法需要保存特征線與穿過的所有網格相交所得的線段長度、網格的材料索引、網格的體積等信息。因此，需耗費大量的時間以及內存空間來存儲數據[6]。模塊化幾何法分為柵元模塊化[8]和組件模塊化方法[9]，其特征線在重復結構內部產生，按照一定規(guī)則（即特征線追蹤過程）將特征線連接并延伸到整個求解域，只需保存幾類不同重復結構的幾何信息，從而節(jié)省大量的存儲空間，DRAGON、PEACH、SONG、MOCHA_2D1D、Tiger-3D等程序均采用了此方法。

本文以柵元模塊化幾何法為基礎，開發(fā)了特征線輸運程序，分析比較了兩種循環(huán)特征線掃描方法以及特征線間距、網格大小、極角數目、角度求積組、方位角布置對計算精度的影響。

1 特征線法數值求解模型及幾何預處理方法

1.1 步特征線法理論

步特征線方法[10]基于平源近似、中子源項和截面恒定的假設下，對中子輸運方程解析求解，得到中子標通量如式(1)所示：

極角、第m方位角、第k條特征線的角平均通量變化值和入射通量；為i區(qū)g群總截面；Si為平源區(qū)面積；lk,i,m,n為i區(qū)、第n極角、第m方位角、第k條特征線的長度；Qi,g為i區(qū)g群總源項；Dm為第m方位角的特征線間距；θn、ωn、ωm分別是極角大小、極角權重、方位角權重。

由式(1)可知，標通量的求解是所有特征線段入射角通量與出射角通量差值的累積加權的過程。因此，在幾何預處理過程中獲取所有特征線段的幾何信息，并確定其求解順序，只需計算入射角通量與出射角通量差值并累加，無需計算線平均角通量以及保存特征線段的出射角通量和平均角通量，從而極大簡化求解過程和節(jié)省內存。

1.2 幾何預處理方法

傳統(tǒng)的確定論反應堆分析方法及程序通常只能夠處理最常見的幾何類型，例如二維X-Y、四邊形、六邊形、一維圓柱等[11]，對于復雜或任意幾何，特征線方法同傳統(tǒng)確定論方法一樣在實際使用上需提高其幾何的適用性。特征線的幾何預處理目的是在平源區(qū)內布置合適的特征線并能準確方便地處理邊界條件，從而快速精確地進行通量求解。其大體上可以分為兩個部分：特征線的產生與特征線的追蹤。

特征線的產生與柵元網格劃分有關。通常柵元網格由同心圓、十字形、交叉形等基本圖形（圖1）組合而成。柵元內材料越復雜，網格劃分越細密，形成的平源近似區(qū)越多。另外，大多數反應堆在徑向上柵元布置結構對稱，且存在大量相同幾何結構柵元。如圖2所示的C5G7基準題，堆芯只由6種不同柵元構成，且?guī)缀谓Y構高度對稱。因此，可以充分考慮以上因素合理布置柵元內的特征線，例如模塊化幾何方法。

圖1 柵元網格劃分基本圖形Fig.1 Basic graphics of cell meshing.

特征線的追蹤常使用長特征線法和循環(huán)特征線法。早期，對于較大規(guī)模的輸運問題，通常采用長特征線法，例如CASMO-4程序就使用了該方法[12-13]。如圖3(a)所示，特征線與途中所有網格相交所得的幾何信息（如特征線段長度、截面等）都需要保存。通常為了滿足精度上的要求，網格劃分必須足夠小，特征線的布置也必須足夠致密，這都會影響計算效率和內存消耗。因此目前大多使用循環(huán)特征線法，例如LATC程序開發(fā)的DNPL (Direct neutron path linking)循環(huán)特征線技術[14]，該方法利用模塊化特征線布置具有關聯性的特點，即在任意相鄰幾何邊界上，特征線首尾連接[10]，通過合理的特征線方位角布置（圖3(b)，邊界處每個點射出的兩條特征線的角度互為補角），進行循環(huán)掃描，不僅降低了柵元邊界處信息存儲所需的內存消耗而且便于邊界條件的處理。針對循環(huán)特征線的特點且結合邊界條件，在§2介紹了兩種循環(huán)特征線追蹤掃描方法：首尾相間循環(huán)特征線法和首尾相接循環(huán)特征線法。

圖2 C5G7基準題1/4堆芯幾何布置圖Fig.2 Layout of a quarter of 2D C5G7 benchmark.

圖3 長特征線(a)和循環(huán)特征線(b)追蹤掃描示意圖Fig.3 Scanning map of long (a) and cycle (b) characteristics.

另外，角度求積組也是特征線預處理中重要的內容，它包括極角與方位角求積組的選取。其中，方位角選取及權重采用文獻[10]中提供的方法。在極角求積組方面，特征線方法常采用兩種求積組：Gauss-Legendre (GL)和Leonard’s optimum (LO)，但LO求積組專用于特征線法，廣泛用于眾多MOC程序中[15-16]。雖然LO求積組在反射邊界條件下計算柵元或者組件問題時有較高計算精度，但是文獻[17]指出LO求積組在計算某些問題時計算精度會降低，并給出了Tabuchi and Yamamoto’s optimum (TY)求積組，因此本文將在數值驗證部分比較LO和TY求積組對計算精度的敏感性，LO和TY求積組如表1所示。

表1 LO和TY求積組Table1 LO and TY angle quadrature set.

2 循環(huán)特征線掃描法

2.1 首尾相間循環(huán)特征線法

特征線掃描順序如圖4所示，從x1,1點開始（該點可任意選取，建議選取邊界上的點），沿特征線掃描到x1,2點，碰到幾何邊界，即完成第一條特征線的追蹤；第二條特征線的起點為圖4中的x1,3，即前一條特征線的終點與后一條特征線的起點處于兩個對應的邊界上（即左邊界對應右邊界（此時兩個點的縱坐標相等，橫坐標之差為幾何的長度）；上邊界對應下邊界（此時兩個點的橫坐標相等，縱坐標之差為幾何的寬度））。然后追蹤到x1,4點，碰到幾何邊界，即完成第二條特征線的追蹤。依次類推，直至第i條特征線的起點與該組特征線的第一條特征線的起點為同一點時即為完成一組特征線掃描。然后更換起點，循環(huán)追蹤掃描第二組特征線，如圖4中的x2,1追蹤到x2,6。按以上規(guī)律進行追蹤掃描，直至所有特征線掃描完成。

圖4 首尾相間循環(huán)法掃描示意圖Fig.4 Layout of the aft disconnect loop scanning method.

需要指出的是：以上特征線掃描過程只完成了第一象限內方位角α（圖5(a)）的掃描，由于幾何結構對稱，它的補角β=π-α（圖5(b)），第三象限空間的方位角π+α（圖5(c)），它的補角π+β（圖5(d)）的掃描跟方位角α相同，因此只需要重復以上操作就可以完成一個方位角下2π空間的掃描。

圖5 首尾斷開循環(huán)法方位角掃描順序Fig.5 Azimuth sequence of the aft disconnect loop scanning method.

2.2 首尾相接循環(huán)特征線法

如圖6所示，首尾相接特征線方法與首尾相間循環(huán)特征線法不同的地方在于相鄰兩條特征線的連接方式不同：首尾相接特征線方法在第一條特征線追蹤到邊界時，第二條特征線的起點為第一條特征線的終點，且方位角為第一條特征線方位角的補角（即如圖6中的x1,2點到x1,3點），依次類推，直至回到起點。為實現2π空間的循環(huán)，只需要進行兩次循環(huán)，即第一次正方向循環(huán)，如圖7(a)所示，再進行反方向循環(huán)，如圖7(b)所示。

圖6 首尾相接循環(huán)法掃描示意圖Fig.6 Layout of the head-tail circulation scanning method.

圖7 首尾相接循環(huán)法方位角掃描示意圖Fig.7 Azimuth sequence of the head-tail circulation scanning method.

2.3 兩種掃描方法的比較

2.3.1 一維7群例題

本文采用文獻[10]中描述的一維7群例題，幾何模型如圖8所示。其幾何尺寸、燃料布置以及截面均與文獻相同，其左側為真空邊界，右側為全反射邊界。該例題的參考解由文獻[10]給出，采用Leonard求積組、極角數目設置為2、以及4個方位角（分別為7.8°、29.4°、51.8°、87.6°），柵元特征線總條數為25條。其計算結果：參考解(keff)為1.02156、首尾相接循環(huán)掃描法（keff/誤差）為1.018858/0.00264、首尾相間循環(huán)掃描法（keff/誤差）為1.023735/0.00212?？梢?，本文所開發(fā)的程序計算結果與參考解吻合較好，驗證了本文開發(fā)的程序的正確性。

圖8 一維7群例題幾何模型Fig.8 Layout of 1D 7G problem.

2.3.2 二維C5G7-MOX基準題

二維C5G7基準題是一個由兩個UO2燃料組件和兩個密度不同的MOX (Mixed oxide fuel)燃料組件構成的7群基準題，其堆芯幾何結構如圖2所示。該例題具有強泄露、非均勻性強等特點，目前被美、日、韓等國研究機構廣泛用于新一代堆芯物理分析方法的檢驗[9]。該基準題的截面、幾何尺寸、邊界條件均由文獻[18]給出?；诙SC5G7模型，將邊界條件均改為全反射邊界（目的在于簡化邊界處入射通量的更新，降低計算時間），比較了首尾相間循環(huán)掃描方法和首尾相接循環(huán)掃描方法下的無限增殖系數kinf，其結果如表2所示。表2的參考解由DRAGON程序計算給出，其參考值為：kinf=1.267403。

從表2可以看出，首尾相間循環(huán)掃描方法相比于首尾相接循環(huán)掃描方法，它收斂速度很慢，且計算精度稍低（首尾相間循環(huán)方法在邊界條件處理上存在部分偏差，特征線段的出射通量和入射通量存在不統(tǒng)一的現象）。

表2 首尾相間循環(huán)掃描法與首尾相接循環(huán)掃描法的無限增殖系數kinf的比較（收斂準則為10-6）Table2 Comparison of kinf for 2D C5G7 benchmark between fore and aft disconnect loop scanning method and the head-tail circulation scanning method (convergence criterion of keff is 10-6).

3 數值驗證

本部分內容將開展網格大小、特征線間距、角度求積組、極角數目、方位角大小等因素對計算精度的影響。圖9(a)是基于不同網格大小計算一維7群例題的計算結果；同時設置不同方位角數目、方位角大小，分別使用LO和TY極角求積組以及不同方位角大?。ㄔO置三組不同方位角大小的組合形式：case 1: 29.4°、51.8°、64.89°、87.6°；case 2: 10.8°、25.2°、64.8°、79.2°；case 3: 5.0°、7.8°、29.4°、64.8°）計算二維C5G7基準題，其結果分別如表3、圖9(b)、表4所示。

圖9 不同網格大小(a)和不同特征線間距大小(b) keff誤差分析Fig.9 Comparison of keff among different mesh sizes (a) and characteristic lines distances (b).

表3 不同方位角數目keff誤差分析Table3 Comparison of keff under different azimuth number.

表4 不同極角求積組、不同方位角大小、不同極角大小下的keff值（收斂準則為10-6）Table4 Comparison of keff for 2D C5G7 benchmark among different angle quadrature sets and different azimuth cases (convergence criterion of keff is 10-6).

由以上計算結果可以得出以下結論：

1) 隨著網格劃分的越小，計算精度逐漸提高，且當網格在小于或等于2.5 cm時，計算結果的相對誤差低于5‰；

2) 計算精度隨特征線間距的增加呈現先增加后降低的趨勢，且在0.1 cm處計算精度最好（特征線數量太多，不僅會降低計算速度并增加內存消耗，同時平源區(qū)劃分容易重疊從而降低計算精度。因此，要合適設置特征線數目）；

3) 方位角數目增加能提高計算精度，4個方位角時計算精度最高（因為計算精度已經較高，計算結果的截取誤差導致8和12個方位角數目下的計算精度略低于4個方位角，因此，4個方位角數目能滿足計算精度的需求）；

4) 增加極角數目能夠提高計算精度，兩個極角的計算精度已經較高；方位角應在45°處對稱布置；從整體情況來看，采用TY極角求積組加權比LO極角求積組加權計算精度稍高。

4 結語

本文介紹了特征線法幾何預處理的主要內容，基于特征線步特征線法，開發(fā)了特征線中子輸運程序。通過計算一維7群例題和二維C5G7基準題驗證了所開發(fā)的源程序的正確性，并分析比較了首尾相間循環(huán)特征線掃描法和首尾相接循環(huán)特征線掃描法。數值結果表明前者的收斂速度很慢。此外，分析了網格大小、特征線間距、角度求積組、極角數目、方位角布置對計算精度的影響。結果表明，要合理布置特征線數目、網格大小以及方位角和極角，兼顧計算精度和計算時間。

1 劉紫靜, 謝金森, 何麗華, 等. 全堆芯中子輸運計算方法研究現狀[J]. 科技資訊, 2015, 3(21): 1-4.

LIU Zijing, XIE Jinsen, HE Lihua, et al. The research status of whole core neutron transport calculation method[J]. Science amp; Technology Information, 2015, 3(21): 1-4.

2 陳珍平, 王電喜, 何桃, 等. 基于CAD技術的特征線中子輸運計算程序開發(fā)[J]. 核科學與工程, 2012, 32(4): 354-359. DOI: 10.3969/j.issn.0258-0918.2012.04.011.

CHEN Zhenping, WANG Dianxi, HE Tao, et al. Development of a CAD-based neutron transport code with the method of characteristics[J]. Chinese Journal of Nuclear Science and Engineering, 2012, 32(4): 354-359. DOI: 10.3969/j.issn. 0258-0918.2012.04.011.

3 Cho J Y, Joo H G, Kim K S. Cell based CMFD formulation for acceleration of whole-core method of characteristics calculation[J]. Nuclear Engineering and Technology, 2002, 34(3): 250-258.

4 Jevremovic T, Vujic J, Tsuda T. ANEMONA: method of characteristics based on the R-function 2D solid modeler[C]. ANS/ENS International Winter Meeting and Embedded Topical Meetings, Washington DC, USA, 1998.

5 Jevremovic T, Vujic J, Tsuda K. ANEMONA-a neutron transport code for general geometry reactor assemblies based on the method of characteristics and R-function solid modeler[J]. Annals of Nuclear Energy, 2001, 28(2): 125-152.

6 Cho N Z. Fundamentals and recent developments of reactor physics methods[J]. Nuclear Engineering and Technology, 2005, 37(37): 25-78.

7 Yamamoto A, Endo T, Tabuchi M, et al. AEGIS: an advanced lattice physics code for light water reactor analyses[J]. Nuclear Engineering and Technology, 2010, 42(5): 500-519.

8 田超, 鄭友琦, 李云召, 等. 基于組件模塊化特征線方法的中子輸運計算研究[J]. 核動力工程, 2015, 3: 157-161.

TIAN Chao, ZHENG Youqi, LI Yunzhao, et al. Neutron transport study based on assembly modular ray tracing MOC method[J]. Nuclear Power Engineering, 2015, 3: 157-161.

9 湯春桃, 張少泓. 以柵元為模塊進行特征線跟蹤的中子輸運方程解法[J]. 核動力工程, 2009, 30(4): 32-36.

TANG Chuntao, ZHANG Shaohong. Study on method of characteristics based on cell modular ray tracing[J]. Nuclear Power Engineering, 2009, 30(4): 32-36.

10 湯春桃. 中子輸運方程特征線解法及嵌入式組件均勻化方法的研究[D]. 上海: 上海交通大學, 2009.

TANG Chuntao. Investigation on method of characteristics for solving neutron transport equation and embedded assembly homogenization method[J]. Shanghai: Shanghai Jiao Tong University, 2009.

11 葉志民, 張鵬, 劉峙嶸. 基于蒙特卡羅均勻化的堆物理分析方法[J]. 核技術, 2014, 37(6): 060603. DOI: 10.11889/j.0253-3219.2014.hjs.37.060603.

YE Zhimin, ZHANG Peng, LIU Zhirong. Reactor physics analysis method based on Monte Carlo homogenization[J]. Nuclear Techniques, 2014, 37(6): 060603. DOI: 10.11889/ j.0253-3219.2014.hjs.37.060603.

12 Smith K S, Rhodes J D. CASMO-4 characteristics method for two-dimensional PWR and BWR core calculations[C]. ANS/ENS International Winter Meeting and Embedded Topical Meetings, Washington DC, USA, 2003.

13 Zika M R, Adams M L. Transport synthetic acceleration for long-characteristics assembly-level transport problems[J]. Nuclear Science and Engineering, 2000, 134(2): 135-158.

14 Shinya Kosaka, Etsuro Saji. Transport theory calculation for a heterogeneous multi-assembly problem by characteristics method with direct neutron path linking technique[J]. Journal of Nuclear Science and Technology, 2012, 37(12): 1015-1023.

15 Yamamoto A, Tabuchi M, Sugimura N, et al. Non-equidistant ray tracing for the method of characteristics[C]. Topical Meeting on Mathematics and Computation, Super computing, Reactor Physics and Nuclear and Biological Applications (Mamp;C 2005), Avignon, France, 2005.

16 Leonard A, McDaniel C T. Optimal polar angles and weights for the characteristics method[J]. Transactions of the American Nuclear Society, 1995: 73.

17 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 Science and Technology, 2007, 44(2): 129-136.

18 Hong-chul Kim, Chi Young Han, Jong Kyung Kim. Twodant solutions for the 2-D C5G7 MOX benchmark[J]. Progress in Nuclear Energy, 2004, 45(2-4): 201-213.

Geometry preprocessing methods research of MOC

LI Xuehui LIU Zijing XIE Jinsen YU Tao CHEN Zhenping XIE Qin ZENG Wenjie HE Lihua
(School of Nuclear Science and Technology, University of South China, Hengyang 421001, China)

Background: Geometry preprocessing, such as mesh division, between the characteristic line, polar angle quadrature set, the number of polar and azimuth angles and so on, has a great impact on computational accuracy and efficiency of the method of characteristics (MOC). Purpose: This paper based on step characteristics method (SC) has developed a neutron transport calculation code. Methods: Through numerical calculations, it has been verified the correctness of the program and has analyzed and compared two characteristic line scan method (fore and aft disconnect loop scanning method, the head-tail circulation scanning method). In addition, based on head-tail circulation scanning method, conducted a series of sensitivity analysis of the calculation accuracy about mesh division, between the characteristic line, polar angle quadrature set, the number of polar and azimuth angles. Results: Program is accurate and reliable and the head-tail circulation scanning method has faster convergence rate than fore and aft disconnect loop scanning method. Conclusion: Mesh division less than 2.5 cm can meet the calculation accuracy requirements and characteristic line spacing 0.1 cm and four azimuth has the highest calculation precision. More number of polar angle has higher calculation accuracy, two discrete polar angle has enough meet the accuracy requirements; the four arranged azimuth in the first quadrant should be laied out with two angles at both sides of 45° angle; Tabuchi and Yamamodo (TY) quadrature weighted have slightly higher calculation accuracy than Leonard (LO) quadrature set.

MOC, Geometry preprocessing, Circulation characteristic line scan

YU Tao, E-mail: yutao29@sina.com

TL329

10.11889/j.0253-3219.2017.hjs.40.040603

No.11305088)、湖南省教育廳科學研究項目(No.16C1365)資助

李雪輝，男，1990年出生，2014年畢業(yè)于南華大學，現為碩士研究生，研究方向為反應堆物理

于濤，E-mail: yutao29@sina.com

2016-11-30，

2017-01-24

Supported by National Natural Science Foundation of China (No.11305088), Research Foundation of Education Bureau of Hunan Province (No.16C1365) First author: LI Xuehui, male, born in1990, graduated from University of South China in 2014, master student, focusing on reactor physics

Received date: 2016-11-30, accepted date: 2017-01-24