三維MOC的球床堆芯CMFD加速研究

文宇晨， 郝琛,2， 郭建， 王毅箴

(1.哈爾濱工程大學(xué) 核安全與仿真技術(shù)重點(diǎn)學(xué)科實(shí)驗(yàn)室，黑龍江 哈爾濱 150001；2.中國核電工程有限公司， 北京 100840；3.清華大學(xué) 核能與新能源技術(shù)研究院， 北京 100084)

目前球床式反應(yīng)堆的物理分析軟件有：基于中子擴(kuò)散理論的VSOP[1]；在VSOP基礎(chǔ)上優(yōu)化后采用更先進(jìn)兩步法進(jìn)行高溫氣冷反應(yīng)堆物理和燃料循環(huán)計(jì)算的PANGU[2]；美國首選球床堆物理分析軟件PEBBED[3]能夠在中子擴(kuò)散計(jì)算的迭代過程中求解穩(wěn)態(tài)情況下的中子通量、燃耗、功率等反應(yīng)堆參數(shù)。但對反應(yīng)堆堆芯提供精細(xì)的通量分布可以降低保守安全裕量估計(jì)，提高反應(yīng)堆經(jīng)濟(jì)效益，因此對球床式反應(yīng)堆進(jìn)行精細(xì)化中子輸運(yùn)計(jì)算是十分必要的。

郭建等[4]針對球床反應(yīng)堆特殊的堆芯結(jié)構(gòu)，以構(gòu)造實(shí)體幾何方法(constructive solid geometry, CSG)為基礎(chǔ)，建立精細(xì)的球床反應(yīng)堆模型，并用三維特征線方法(three dimension-method of characteristics, 3D-MOC)求解中子輸運(yùn)方程，開發(fā)了準(zhǔn)確描述球床高溫氣冷堆幾何的三維特征線程序(method of characteristics for pebble-bed HTR, MOCP)。為解決3D-MOC方法計(jì)算時(shí)間長和迭代次數(shù)多等問題，MOCP使用稀疏長特征線方法(macro-track transport acceleration method, MTTA)作為3D-MOC的加速算法。MTTA相較基于中子擴(kuò)散方程的加速算法需要更大的計(jì)算量，而基于中子擴(kuò)散方程的粗網(wǎng)有限差分方法[5](coarse mesh finite difference, CMFD)加速效果好且便于實(shí)施，通常作為輕水堆精細(xì)化輸運(yùn)計(jì)算的加速算法。但傳統(tǒng)CMFD可能出現(xiàn)耦合擴(kuò)散系數(shù)為負(fù)進(jìn)而降低CMFD矩陣對角占優(yōu)性，導(dǎo)致收斂不穩(wěn)定或收斂發(fā)散的現(xiàn)象?；趶V義等價(jià)理論的CMFD方法[6-7](generalized equivalence theory-CMFD, gCMFD)通過定義節(jié)塊不連續(xù)因子與擴(kuò)散系數(shù)修正因子，保證了CMFD矩陣對角占優(yōu)性進(jìn)而解決傳統(tǒng)CMFD中可能出現(xiàn)的收斂不穩(wěn)定問題。球床反應(yīng)堆有極強(qiáng)的非均勻性，為避免在球床堆出現(xiàn)CMFD收斂困難或發(fā)散的現(xiàn)象，將gCMFD作為3D-MOC的加速算法,并使用廣義最小殘差算法[8-9]求解gCMFD方程，開發(fā)了gCMFD加速模塊與MOCP程序耦合。

球床堆中存在大量隨機(jī)分布的燃料球和冷卻劑，對粗網(wǎng)格獲取均勻化信息造成極大困難。為解決該問題，本文將燃料球細(xì)網(wǎng)格與粗網(wǎng)格進(jìn)行距離判斷以獲得粗網(wǎng)內(nèi)的平源區(qū)，同時(shí)針對MOCP中使用的CSG建模方法獲得每個(gè)粗網(wǎng)內(nèi)的冷卻劑流平源區(qū)。最終得到gCMFD線性系統(tǒng)所需均勻化截面等參數(shù)。在實(shí)現(xiàn)粗網(wǎng)與平源區(qū)嚴(yán)格對應(yīng)的過程中，可能出現(xiàn)平源區(qū)被分割的現(xiàn)象，文獻(xiàn)[10]程序OpenMOC針對該問題有相應(yīng)研究。

本文在數(shù)值驗(yàn)證部分，采用簡單輕水堆基準(zhǔn)題KUCA[11]對gCMFD的正確性與加速效果進(jìn)行驗(yàn)證。自編直角幾何球床堆芯和圓柱幾何球床堆芯，初步驗(yàn)證gCMFD方法在不同幾何下球床處理的正確性；高階計(jì)算與低階計(jì)算的等價(jià)性與加速效果。

1 gCMFD方程組與3D-MOC算法

gCMFD方法為3D-MOC提供更好的初值即中子通量密度和有效增殖因子以減少輸運(yùn)計(jì)算的迭代次數(shù)；同時(shí)，3D-MOC將平源區(qū)信息傳遞給gCMFD系統(tǒng)計(jì)算其均勻化截面。

1.1 三維特征線方法

將穩(wěn)態(tài)下的中子輸運(yùn)方程改寫為特征線形式，同時(shí)將梯度算符寫為沿特征線的方向?qū)?shù)為：

(1)

源項(xiàng)為：

(2)

式中：E、Ω分別為能量、角度變量；χ、Σs、Σf分別為裂變能譜、散射截面與裂變截面。

將幾何模型沿3個(gè)方向進(jìn)行劃分得到每個(gè)平源區(qū)網(wǎng)格，在能量、角度和材料方面引入不同的近似后，每個(gè)平源區(qū)網(wǎng)格的標(biāo)通量表示為：

(3)

(4)

在MOC中，角通量為解析求解為：

Ψi,g,m,k(s1)=Ψi,g,m,k(s0)exp(-τi,g,m,k)+

(5)

式中：s1、s0為該平源區(qū)網(wǎng)格內(nèi)特征線段的起點(diǎn)和終點(diǎn)；τi,g,m,k為光學(xué)厚度。

當(dāng)3D-MOC求解完畢，通過平源區(qū)信息獲得每個(gè)粗網(wǎng)格的均勻化截面并更新gCMFD線性系統(tǒng)特征值為：

(6)

1.2 圓柱幾何下gCMFD

CMFD方法的本質(zhì)是求解中子擴(kuò)散方程。在圓柱幾何下，對空間進(jìn)行離散得到的網(wǎng)格如圖1所示。

圖1 圓柱幾何空間網(wǎng)格Fig.1 Cylindrical meshes

gCMFD中，2節(jié)塊間的中子流計(jì)算式為：

(7)

(8)

(9)

式中：f為節(jié)塊不連續(xù)因子；g為擴(kuò)散系數(shù)修正因子；h為相鄰節(jié)塊中心之間的距離；α為邊界處的邊界條件。真空邊界條件時(shí)，α=0.5；全反射邊界條件時(shí)，α=0。

如圖1所示，圓柱幾何下沿徑向的粗網(wǎng)格會(huì)呈現(xiàn)逐漸變大的趨勢，進(jìn)而導(dǎo)致網(wǎng)格間的交界面積與距離發(fā)生變化，因此不同方向的交界面積與距離計(jì)算式分別為：

徑向：

(10)

(11)

周向：

Ac,E=ΔrcΔzc,hc,E=r′cΔθc,hE,c=r′cΔθE

(12)

Ac,W=ΔrcΔzc,hc,W=r′cΔθc,hW,c=r′cΔθW

(13)

軸向：

Ac,B=r′cΔrcΔθc,hc,B=Δzc,hB,c=ΔzB

(14)

Ac,T=r′cΔrcΔθc,hc,T=Δzc,hT,c=ΔzT

(15)

其中，r′c表示該網(wǎng)格的幾何中心：

(16)

幾何中心表示某一規(guī)則網(wǎng)格最中心的位置，并通過幾何中心的定義準(zhǔn)確獲得圓柱幾何下離散網(wǎng)格間的長度，如圖2。ri表示沿徑向第i個(gè)網(wǎng)格邊界中心的半徑值。

圖2 扇形網(wǎng)格幾何中心示意Fig.2 Geometry center of shell mesh

同時(shí)，周向網(wǎng)格在邊界處呈周期性邊界條件，即：周向的起點(diǎn)網(wǎng)格與終點(diǎn)網(wǎng)格相鄰，此時(shí)網(wǎng)格的耦合擴(kuò)散系數(shù)用式(8)計(jì)算；徑向第1層的網(wǎng)格，沿“SOUTH”方向沒有交界面積，當(dāng)前網(wǎng)格沿“SOUTH”方向的耦合擴(kuò)散系數(shù)為0。

在gCMFD方法中強(qiáng)制高階系統(tǒng)的中子流等于低階擴(kuò)散系統(tǒng)的中子流并以此計(jì)算NDF與MDF[6，12]。根據(jù)算法1在MOCP中加入新的功能以判斷每條特征線段穿過粗網(wǎng)的方向。

算法1：find Boundry Segments CrossCM

Intersections←findIntersection(segment,CM)

ifIntersections is 1then

forallPlane∈CMdo

Length1←minDistancePointToSurface(segment,Plane)

endfor

elseifIntersections is 2then

forallQuadratic∈CMdo

Length2←minDistancePointToSurface(segment,Quadratic)

endfor

minLength←min(Length1,Length2)

Direction←setSegmentDirection(minLength,segment)

end

獲得每個(gè)粗網(wǎng)格的節(jié)塊不連續(xù)因子與擴(kuò)散系數(shù)修正因子后可根據(jù)式(8)、(9)計(jì)算耦合擴(kuò)散系數(shù)，最終得到gCMFD的線性系統(tǒng)：

(17)

2 球床模型處理

在球床堆中實(shí)現(xiàn)gCMFD加速的關(guān)鍵在于保證平源區(qū)網(wǎng)格與粗網(wǎng)的嚴(yán)格對應(yīng)。球床堆沒有類似輕水堆中柵元的結(jié)構(gòu)，而在球床內(nèi)存在大量隨機(jī)分布的燃料球及冷卻劑流，因此在球床堆獲得粗網(wǎng)格內(nèi)平源區(qū)信息是極為困難的。由于MOCP使用CSG建模方法的優(yōu)勢，各個(gè)幾何所具有的面包含該幾何相應(yīng)的關(guān)鍵信息，為實(shí)現(xiàn)球床堆內(nèi)的gCMFD加速提供了可能。

2.1 構(gòu)造實(shí)體幾何

構(gòu)造實(shí)體幾何是一種通過簡單模型的“布爾”運(yùn)算來構(gòu)造復(fù)雜幾何體的常見實(shí)體建模方法。在MOC的特征線追蹤過程中計(jì)算點(diǎn)、線、面等所需要的參數(shù)信息均可通過CSG構(gòu)造的復(fù)雜幾何體所包含的面得到，因此CSG非常適合作為MOC的幾何建模方法。

CSG方法建模的過程為：1)定義平面并表示相應(yīng)半空間；2)對半空間進(jìn)行“布爾”運(yùn)算構(gòu)造簡單幾何；3)將簡單幾何進(jìn)行“布爾”運(yùn)算構(gòu)造復(fù)雜幾何。

在三維空間中，某一平面或曲面將三維空間分為2部分：其中F+表示該面所定義出的正半空間，F(xiàn)-表示該面定義出的負(fù)半空間，如圖3所示。

圖3 三維面與半空間構(gòu)成Fig.3 Construct of 3D surface and half space

通過對半空間進(jìn)行“布爾”運(yùn)算可構(gòu)成不同的簡單幾何體；對簡單幾何進(jìn)行“布爾”運(yùn)算即可獲得所需要的復(fù)雜幾何。需要注意的是，為了保證建模的正確性，盡量通過不同的簡單幾何來構(gòu)成復(fù)雜幾何，并盡可能減少“布爾”運(yùn)算的次數(shù)。在CSG方法中，復(fù)雜幾何的最基本單位是各個(gè)不同的面，MOCP程序針對球床堆芯的復(fù)雜模型，提供不同類型的面，進(jìn)而提高了描述幾何的能力。CSG方法構(gòu)造復(fù)雜幾何的思路如圖4所示。

圖4 復(fù)雜幾何構(gòu)成Fig.4 Construct of complex geometry

通過CSG，球床堆芯內(nèi)的燃料球由不同球面定義出的負(fù)半空間構(gòu)成；堆芯容器由不同面方程定義出的半空間經(jīng)過“布爾”運(yùn)算構(gòu)成；堆芯內(nèi)的冷卻劑流由堆芯容器與球面定義出的正半空間取“交”運(yùn)算組成。半空間均保存相應(yīng)模型的幾何信息，如：球心位置、燃料球半徑和平面參數(shù)，為粗網(wǎng)格獲得平源區(qū)信息打下基礎(chǔ)。

2.2 燃料球及球床間隙處理

在輕水堆中，粗網(wǎng)格主要是一個(gè)個(gè)柵元，對每個(gè)柵元進(jìn)行劃分后得到粗網(wǎng)格內(nèi)的平源區(qū)。而在圓柱幾何球床堆中，只能對整體圓柱幾何按照軸向、周向和徑向進(jìn)行劃分，得到虛擬的粗網(wǎng)格。因此，如何獲得虛擬粗網(wǎng)格中的平源區(qū)及相應(yīng)信息就是實(shí)現(xiàn)球床堆內(nèi)gCMFD加速的關(guān)鍵。

在MOCP中，燃料球的處理方法通常是將燃料球沿徑向劃分后再進(jìn)行卦限劃分得到形狀為“BALL”或“SHELL”的平源區(qū)網(wǎng)格。根據(jù)這些網(wǎng)格的面參數(shù)與粗網(wǎng)進(jìn)行距離關(guān)系判斷，若為“相交”或“屬于”關(guān)系，則認(rèn)為該細(xì)網(wǎng)格或部分細(xì)網(wǎng)格屬于當(dāng)前粗網(wǎng)格的平源區(qū)，距離關(guān)系l為：

(18)

式中：A、B、C、D表示平面參數(shù)；x0、y0、z0表示球心坐標(biāo)。

獲得燃料球與粗網(wǎng)各個(gè)面的距離后，通過半空間判斷燃料球是否與粗網(wǎng)存在交集，二維示意如圖5所示。

編號(hào)為1和6的燃料球與粗網(wǎng)各個(gè)面的距離均小于半徑，因此認(rèn)為屬于該粗網(wǎng)；編號(hào)為3的燃料球球心在粗網(wǎng)內(nèi)，且與粗網(wǎng)距離小于半徑，因此認(rèn)為該燃料球與粗網(wǎng)“相交”并將部分粗網(wǎng)面賦予該燃料球；編號(hào)為4的燃料球球心在粗網(wǎng)外，由式(18)可知，二者距離小于燃料球半徑，將粗網(wǎng)的圓柱面賦予燃料球。

圖5 粗網(wǎng)與平源區(qū)對應(yīng)Fig.5 Map FSR and course mesh

當(dāng)所有粗網(wǎng)與平源區(qū)均遍歷完畢，對粗網(wǎng)格內(nèi)形狀為“BALL”或“SHELL”平源區(qū)的外層半空間取反，并加上粗網(wǎng)格的半空間，組成該粗網(wǎng)格內(nèi)的冷卻劑流平源區(qū)，見算法2與算法3。最終每個(gè)粗網(wǎng)格均包含若干燃料球平源區(qū)與一個(gè)冷卻劑流平源區(qū)。

算法2：map Fuel Balland Course Mesh

forallcm∈Course Mesh_mapdo

locationType←calcLenthFSRandCM(fsr, cm)

iflocationType is intersectthen

new_fsr←copySectionCell(fsr)

addFSRintoCM(new_fsr)

elseiflocationType is belongthen

new_fsr←copyCell(fsr)

addFSRintoCM(new_fsr)

end

endfor

經(jīng)過以上處理，保證粗網(wǎng)與球床區(qū)域內(nèi)細(xì)網(wǎng)的嚴(yán)格對應(yīng)，但產(chǎn)生了許多不規(guī)則網(wǎng)格。不規(guī)則網(wǎng)格的體積通過特征線積分得到：

(19)

算法3：create Inter Space FSR

forallcm∈CMdo

interspace←createNewCell()

forallfsr∈fsrStoryInCMdo

addSurface(interspace,fsr)

endfor

addSurface(interspace,cm)

addFSRintoCM(interspace_fsr)

endfor

值得注意的是，SHELL類型的平源區(qū)由2個(gè)球面組成，且半空間相反，因此在判斷距離時(shí)要考慮SHELL的空心部分；如圖5所示，與粗網(wǎng)“相交”的燃料球被賦予新的半空間，原始平源區(qū)被分割成多個(gè)新平源區(qū)，此時(shí)需要將舊的平源區(qū)消除以保證特征線追蹤模塊的正確性。

3 加速模塊數(shù)值驗(yàn)證

在數(shù)值驗(yàn)證部分，首先選擇一個(gè)簡單輕水堆模型用于驗(yàn)證gCMFD加速模塊在MOCP中的加速效果與正確性。由于實(shí)際球床堆幾何模型復(fù)雜，本文僅針對球床堆中常見的直角幾何與圓柱幾何建立對應(yīng)的gCMFD線性系統(tǒng)。因此本文自編簡易直角與圓柱幾何球床堆模型初步驗(yàn)證gCMFD加速方法在球床堆中的可行性與加速效果。

3.1 KUCA基準(zhǔn)題模型

KUCA是一個(gè)小型輕水堆的兩群輸運(yùn)問題，其1/4堆芯結(jié)構(gòu)如圖6所示，具體截面信息見表1、2。該基準(zhǔn)題有2個(gè)狀態(tài)：1)無控制棒，且控制棒區(qū)域?yàn)檎婵眨?)控制棒全部插入堆芯。本文選擇第2個(gè)狀態(tài)進(jìn)行驗(yàn)證。

圖6 KUCA基準(zhǔn)題1/4堆芯結(jié)構(gòu)Fig.6 Configuration for CUKA benchmark

第2種狀態(tài)的KUCA基準(zhǔn)題全堆尺寸為50.0 cm×50.0 cm×50.0 cm，布置特征線時(shí)選擇Lee求積組(4階，共24個(gè)方向)。劃分平源區(qū)時(shí)，每個(gè)細(xì)網(wǎng)格寬0.5 cm，3D-MOC的收斂標(biāo)準(zhǔn)為10-5。在加入gCMFD加速模塊時(shí)，對該基準(zhǔn)題的3個(gè)方向進(jìn)行了多種劃分形式，gCMFD的收斂標(biāo)準(zhǔn)為10-8。

表1 KUCA基準(zhǔn)題截面信息Table 1 Cross section information of KUCA benchmark

表2 KUCA基準(zhǔn)題散射截面信息Table 2 Scattering cross section info of KUCA benchmark

在驗(yàn)證過程中使用蒙特卡羅程序?qū)ο嗤Ｐ瓦M(jìn)行計(jì)算作為KUCA基準(zhǔn)題的參考解，特征值為1.038 93，耦合gCMFD加速模塊的MOCP程序詳細(xì)計(jì)算數(shù)據(jù)見表3、4。將gCMFD按不同劃分方式進(jìn)行計(jì)算，所得特征值均與高階MOC計(jì)算保持一致，為1.038 99?？梢钥闯?，gCMFD實(shí)現(xiàn)與高階計(jì)算嚴(yán)格等價(jià)的同時(shí)，極大降低3D-MOC計(jì)算的計(jì)算代價(jià)。

表3 KUCA基準(zhǔn)題計(jì)算特征值數(shù)據(jù)Table 3 Eigenvalues of KUCA benchmark

表4 歸一化通量分布Table 4 Distribution of normalized flux

3.2 自編球床堆模型

為驗(yàn)證gCMFD加速算法在不同幾何下的準(zhǔn)確性與加速效果，同時(shí)對算法2及算法3在球床部分的正確性進(jìn)行驗(yàn)證，自定義直角與圓柱幾何下的球床堆芯模型。

3.2.1 直角幾何球床堆芯

該基準(zhǔn)題為小型直角幾何球床堆的兩群輸運(yùn)問題，堆芯結(jié)構(gòu)如圖7所示，除中心的球床部分，還有外圍的反射層，具體截面信息見表5、6。

圖7 直角幾何球床堆芯模型Fig.7 Cuboid pebble-bed reactor geometry

表5 直角幾何球床堆模型截面信息

表6 直角幾何球床堆模型散射截面信息

該基準(zhǔn)題全堆尺寸為30.0 cm×30.0 cm×30.0 cm，球床尺寸為15 cm×15 cm×15 cm，球床內(nèi)共172個(gè)燃料球，燃料球半徑為2.99 cm，球床內(nèi)的間隙為“void”，即冷卻劑材料。在布置特征線時(shí)選用Lee求積組(4階，共24個(gè)方向)。每個(gè)燃料球沿徑向劃分6層，再進(jìn)行卦限劃分，即單個(gè)燃料球劃分出48個(gè)平源區(qū)。3D-MOC收斂標(biāo)準(zhǔn)為10-5，gCMFD收斂標(biāo)準(zhǔn)為10-8,特征值為1.185 76。詳細(xì)計(jì)算結(jié)果見表7。

表7 自寫直角幾何球床堆特征值計(jì)算數(shù)據(jù)Table 7 Eigenvalues of cuboid pebble-bed geometry

3.2.2 圓柱幾何球床堆芯

為驗(yàn)證圓柱幾何下gCMFD系統(tǒng)的正確性，自編圓柱型球床堆芯，該基準(zhǔn)題為兩群輸運(yùn)問題。堆芯結(jié)構(gòu)如圖8所示，截面信息見表8。

圓柱堆芯半徑為49.462 5 cm，高79.14 cm，球床部分圓柱半徑為32.975 cm，高65.95 cm。球床內(nèi)含有1 080個(gè)燃料球，燃料球半徑為2.99 cm，燃料球間間隙為“void”。特征線布置中，每條特征線寬0.5 cm，燃料球劃分方式同為沿經(jīng)向劃分6層后進(jìn)行卦限劃分。3D-MOC收斂表示為10-5，gCMFD收斂標(biāo)準(zhǔn)位10-8。特征值為1.260 56，其余計(jì)算結(jié)果見表8。

圖8 圓柱幾何自寫球床堆模型Fig.8 Cylindrical pebble-bed reactor geometry

表8 圓柱幾何球床堆特征值計(jì)算結(jié)果Table 8 Eigenvalues of cylindrical pebble-bed geometry

以上數(shù)值結(jié)果表明，gCMFD加速模塊編寫正確，且在保證不同分辨率系統(tǒng)嚴(yán)格等價(jià)的同時(shí)能夠?qū)?D-MOC輸運(yùn)計(jì)算起到良好的加速效果。

4 結(jié)論

1)本文中針對球床的處理方式正確，保證粗網(wǎng)與燃料球細(xì)網(wǎng)格嚴(yán)格對應(yīng)，能夠準(zhǔn)確計(jì)算每個(gè)粗網(wǎng)均勻化截面等信息的同時(shí)，為細(xì)網(wǎng)格提供了更好的初值以加速收斂。

2)初步實(shí)現(xiàn)gCMFD加速方法在球床反應(yīng)堆的應(yīng)用，同時(shí)將gCMFD的幾何限制拓展至圓柱幾何，大幅提高了3D-MOC的計(jì)算效率。未來可針對圓柱幾何特殊結(jié)構(gòu)進(jìn)一步改進(jìn)gCMFD方法從而實(shí)現(xiàn)更快速準(zhǔn)確求解。