劉仕倡，申靖文，嚴(yán)伊蔓，杜 華，孫業(yè)帥，秦 瑤，鄭 俞， 盧 棚，李 夏，張 顯，劉松林，陳義學(xué)

(1.華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206； 2.國家電投集團(tuán)科學(xué)技術(shù)研究院有限公司 核電軟件技術(shù)中心，北京 102209； 3.中國科學(xué)院 等離子體物理研究所，安徽 合肥 230031)

中國聚變工程實(shí)驗(yàn)堆(CFETR)[1]是我國于2011年開始設(shè)計(jì)研究的超導(dǎo)托卡馬克反應(yīng)堆，在DEMO級別聚變電站與ITER間起到了過渡銜接的作用。托卡馬克型的聚變堆堆芯主要包括：第一壁、包層、偏濾器、真空室、輻射屏蔽、生物屏蔽等。包層是聚變堆堆芯的一個(gè)重要組成部分，在強(qiáng)中子輻照下，包層主要有氚增殖、能量轉(zhuǎn)換及輻射屏蔽等作用，因此，包層中子學(xué)的設(shè)計(jì)與聚變堆中的大量問題均相關(guān)。此外，托卡馬克裝置具有非常復(fù)雜的幾何結(jié)構(gòu)，包含各種結(jié)合體并且結(jié)構(gòu)十分精細(xì)，而蒙特卡羅方法幾何描述精細(xì)的特點(diǎn)無疑非常適宜于解決這樣的問題[2]。但如果用于模擬計(jì)算的蒙特卡羅程序不具備高次曲面的描述能力，僅使用一階或二階曲面并不能滿足幾何建模精細(xì)化的需求。那么對于這類復(fù)雜幾何模型的建立就僅能采取幾何近似的方式，如可將托卡馬克裝置的環(huán)形真空室的圓環(huán)面近似為很多段直線，不僅耗時(shí)耗力，還會(huì)導(dǎo)致精度損失，不利于聚變堆的準(zhǔn)確模擬。國際上具有處理聚變堆復(fù)雜幾何能力的蒙特卡羅程序有美國Los Alamos實(shí)驗(yàn)室的MCNP[3]、日本原子能機(jī)構(gòu)的MVP[4]、法國原子能委員會(huì)的TRIPOLI[5]等。國內(nèi)方面，中國科學(xué)院核能安全技術(shù)研究所的SuperMC[6]在磁約束核聚變裝置的蒙特卡羅模擬方面進(jìn)行了大量研究，北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所的JMCT[7]也進(jìn)行了一些研究。其中，MCNP是應(yīng)用最為廣泛的經(jīng)典蒙特卡羅程序，其復(fù)雜曲面處理方法被廣泛應(yīng)用于國內(nèi)外蒙特卡羅程序，然而對該方法的數(shù)學(xué)模型和計(jì)算原理進(jìn)行系統(tǒng)推導(dǎo)和闡述的研究報(bào)道較少。

本文基于自主化蒙特卡羅程序cosRMC，研究及闡述蒙特卡羅復(fù)雜曲面建模的數(shù)學(xué)模型和計(jì)算方法，并對原有的幾何模塊進(jìn)行拓展，使其具備對復(fù)雜幾何模型的描述能力，最后將其應(yīng)用于CFETR三維模型的聚變中子學(xué)分析中，對聚變堆包層設(shè)計(jì)時(shí)的3個(gè)關(guān)鍵的中子學(xué)參數(shù)進(jìn)行計(jì)算分析。

1 cosRMC復(fù)雜幾何建模功能

1.1 cosRMC及其幾何建模功能介紹

cosRMC是由清華大學(xué)工程物理系反應(yīng)堆工程計(jì)算分析實(shí)驗(yàn)室(REAL團(tuán)隊(duì))與國家電投集團(tuán)科學(xué)技術(shù)研究院有限公司核電軟件技術(shù)中心聯(lián)合研發(fā)的以蒙特卡羅粒子輸運(yùn)方法為理論基礎(chǔ)的自主化蒙特卡羅程序[8]。主要用于核工程領(lǐng)域中固定源計(jì)算、臨界計(jì)算、燃耗計(jì)算、瞬態(tài)計(jì)算和均勻化并群計(jì)算等。

cosRMC采用文本建模方式，基于構(gòu)造幾何實(shí)體(CSG)對系統(tǒng)幾何結(jié)構(gòu)進(jìn)行描述。即通過對交界面不同指向的區(qū)域進(jìn)行布爾運(yùn)算，如交、并和非等運(yùn)算，構(gòu)成對體的描述。目前，cosRMC程序可進(jìn)行描述的面包括一般平面、任意球面、軸心位于坐標(biāo)軸上或平行于坐標(biāo)軸的圓柱面或圓錐面等標(biāo)準(zhǔn)形式的曲面，所有曲面的曲面方程均使用解析方程表示方式。在本文工作之前，cosRMC無法實(shí)現(xiàn)對一些復(fù)雜曲面的描述，如任意位置的圓柱面、圓錐面、橢圓面和圓環(huán)面等。

1.2 復(fù)雜幾何建模功能實(shí)現(xiàn)

粒子在蒙特卡羅程序幾何中的輸運(yùn)需確定粒子沿飛行方向與曲面的交點(diǎn)，以圓環(huán)面為例，通過解析方法簡化為求解一個(gè)四次多項(xiàng)式方程解析解的問題。

1) 中子軌跡與復(fù)雜曲面的聯(lián)立求解

(1)

進(jìn)一步在粒子輸運(yùn)問題中，在射線簇{x=ξ+aX,y=η+BX,z=ζ+γX;-∞

{[(1-β2)+ρβ2]X2+[2(aξ+γζ)+2ρβη-

2ρβy0]+(ξ2+ζ2+ρη2-2ρηy0+B0)}2=

A0[(1-β2)X2+2(aξ+γζ)X+ξ2+ζ2]

(2)

式(2)可轉(zhuǎn)化為：

(GX2+MX+U)2=A0(FX2+LX+T)

(3)

其中，F(xiàn)=1-β2，G=F+ρβ2，L=2(aξ+γζ)，M=L+2ρβ(η-y0)，T=ξ2+ζ2，U=T+ρη(η-2y0)+B0。

所以，射線簇{x=ξ+aX,y=η+BX,z=ζ+γX;X>0}與圓環(huán)面求交點(diǎn)(x,y,z)的問題可轉(zhuǎn)化成一個(gè)求解式(4)所示四次方程的正根X的問題，僅需判斷射線簇是與圓環(huán)面的內(nèi)表面還是外表面存在交點(diǎn)即可。

X4+BX3+CX2+DX+E=0

(4)

2) 復(fù)雜曲面多項(xiàng)式方程的求解

由于舍入誤差的存在，使得利用計(jì)算機(jī)程序求解形如式(4)的四次多項(xiàng)式方程均非常困難。Maisonnier和Everett提出了一種求解“射線簇與三維或四維曲面與交點(diǎn)”的單精度四次算法[9]，該算法可分為兩步，第1步是將粒子穿過圓環(huán)面時(shí)四次多項(xiàng)式方程的求解轉(zhuǎn)化為三次方程的求解。

首先討論三階曲面的解析方程，可利用泰勒展開求解：

f(x)=d+cx+bx2+x3=

(5)

對式(5)求導(dǎo)代換后可得到新的三次解析方程：

g(x)=y3+py+q

(6)

Δ3=-4(27)W

(7)

同樣對四階曲面的解析方程也進(jìn)行泰勒展開：

F(X)=E+DX+CX2+BX3+X4=

(8)

通過推導(dǎo)可得出，四階曲面方程的判別式為Δ4=Δ3。根據(jù)判別式條件就變成了一個(gè)類似于三階曲面方程求解的過程。

第2步是根據(jù)多項(xiàng)式系數(shù)判定來確定是否存在任意正實(shí)根，主要包括3種情況：當(dāng)Δ3>0(W<0)時(shí)，有3個(gè)不同的實(shí)根；當(dāng)Δ3=0(W=0)時(shí)，有3個(gè)含有重復(fù)根的實(shí)根；當(dāng)Δ3<0(W>0)時(shí)，有單個(gè)實(shí)根和1對共軛復(fù)根。

根據(jù)以上理論方法，本文在cosRMC原有的幾何模塊中新添加了橢圓面、雙曲面、拋物面、圓環(huán)面等復(fù)雜曲面，依舊采用助記名的方式存儲(chǔ)，如表1所列。

1.3 功能測試

聚變堆概念設(shè)計(jì)堆PPCS(power plant conceptual study)是2000年于歐洲啟動(dòng)的最為系統(tǒng)的聚變堆研究計(jì)劃，目的是幫助評估聚變能形式和確立歐洲聚變計(jì)劃的可持續(xù)性和優(yōu)先權(quán)，屬于世界上目前較具規(guī)模和代表性的磁約束聚變動(dòng)力反應(yīng)堆概念設(shè)計(jì)系列之一，如圖1所示。作為典型的磁約束聚變堆概念設(shè)計(jì)系列之一，PPCS具有非常復(fù)雜的幾何機(jī)構(gòu)，在使用蒙特卡羅方法對其進(jìn)行聚變中子學(xué)分析時(shí)，需構(gòu)建諸多包含復(fù)雜曲面的幾何體，以及包層模塊重復(fù)結(jié)構(gòu)的描述，如圖1c中標(biāo)注所示。

表1 新增cosRMC復(fù)雜曲面卡Table 1 New cosRMC complex surface card

a——PPCS的cosRMC模型；b——PPCS的MCNP模型；c——PPCS模型中選取的柵元 圖1 PPCS的cosRMC及MCNP的模型示意圖Fig.1 Model diagram of PPCS with cosRMC and MCNP

利用cosRMC程序?qū)PCS進(jìn)行建模(圖1)后進(jìn)行模擬計(jì)算，中子源采用14 MeV的各向同性點(diǎn)源。MCNP和cosRMC模擬的總粒子數(shù)均為1×108。選取了18個(gè)有高次曲面的柵元，圖1c中標(biāo)注的數(shù)字即為選取柵元的柵元號。統(tǒng)計(jì)其體通量，與MCNP的計(jì)算結(jié)果進(jìn)行比較，令相對統(tǒng)計(jì)誤差為：

(9)

σMCNP和σcosRMC為MCNP和cosRMC的計(jì)數(shù)器相對統(tǒng)計(jì)誤差。柵元通量的對比如圖2所示。可看出，cosRMC程序?qū)λ袔缀误w上包含高次曲面的體通量的計(jì)算結(jié)果，與MCNP的計(jì)算結(jié)果吻合良好，相對偏差幾乎均小于10%。且cosRMC單核計(jì)算時(shí)間為250.980 6 min，MCNP單核計(jì)算時(shí)間為258.74 min，計(jì)算效率相當(dāng)。說明了cosRMC程序中實(shí)現(xiàn)的復(fù)雜幾何建模功能的正確性。

圖2 PPCS模型cosRMC與MCNP模擬計(jì)算結(jié)果比較Fig.2 Comparison of simulation results of PPCS between cosRMC and MCNP

2 CFETR建模與聚變中子學(xué)分析

2.1 CFETR建模

CFETR的三維模型體積大、幾何結(jié)構(gòu)復(fù)雜，托卡馬克裝置中包含環(huán)形真空室、環(huán)向線圈、極向線圈等結(jié)構(gòu)部件使得在使用蒙特卡羅方法對其進(jìn)行聚變中子學(xué)分析時(shí)，需構(gòu)建諸多包含四階曲面的幾何體。使用拓展了復(fù)雜幾何建模功能的cosRMC程序?qū)FETR三維模型進(jìn)行精細(xì)化建模，如圖3所示，與MCNP程序構(gòu)建的模型圖進(jìn)行對比，如圖4所示。從模型圖可看出cosRMC程序描述的CFETR三維模型的正確性。

圖3 CFETR的cosRMC模型徑向與縱向剖面圖Fig.3 Radial and longitudinal profiles of cosRMC model of CFETR

圖4 CFETR的MCNP模型縱向剖面圖Fig.4 Longitudinal profile of MCNP model of CFETR

進(jìn)一步對CFETR的包層中子學(xué)中的3個(gè)重要參數(shù)進(jìn)行計(jì)算分析。CFETR的包層模塊如圖5所示，沿環(huán)向上和極向上進(jìn)行分割，分為內(nèi)包層模塊和外包層模塊。中子源采用氘氚聚變中子源，分布于5個(gè)柵元中，中子源強(qiáng)分別為1%、6%、13%、30%和50%。源中子能量為14 MeV的高斯聚變譜，方向各向同性。CFETR的氚增殖比(TBR)、中子壁載荷(NWL)和核熱沉積計(jì)算過程中，聚變功率采用第1階段的運(yùn)行功率200 MW。MCNP和cosRMC模擬的總粒子數(shù)均為6.9×106。MCNP和cosRMC均采用FENDL-2.1核截面數(shù)據(jù)庫進(jìn)行計(jì)算。

圖5 CFETR包層結(jié)構(gòu)編號及位置示意圖Fig.5 Diagram of CFETR blanket structure number and location

2.2 TBR

氘氚聚變反應(yīng)堆能否投入商業(yè)運(yùn)行的首要條件是能實(shí)現(xiàn)氚的自持。聚變中子學(xué)分析計(jì)算時(shí)，增殖材料中6Li、7Li的產(chǎn)氚率為：

T6=NLi?σ6Li(n,α)(r,E)φ(r,E)dEdr

(10)

T7=NLi?σ7Li(n,n′α)(r,E)φ(r,E)dEdr

(11)

其中，σ6Li(n,α)和σ7Li(n,n′α)分別為6Li的(n,α)反應(yīng)截面和7Li的(n,n′α)反應(yīng)截面。那么，堆芯等離子體平均1次氘氚聚變反應(yīng)在包層內(nèi)增殖的TBR為：

(12)

其中，Sn為聚變中子源強(qiáng)。利用cosRMC對其增殖包層各模塊的氚增殖比進(jìn)行計(jì)算，與MCNP的計(jì)算結(jié)果進(jìn)行比較(表2)。可看出，兩程序?qū)Ω髟鲋嘲鼘与霸鲋潮鹊挠?jì)算結(jié)果吻合良好，相對偏差均小于1%。并且，通過模擬計(jì)算可發(fā)現(xiàn)，增殖包層總的TBR為1.2，能滿足氘氚聚變反應(yīng)堆運(yùn)行時(shí)的氚自持需求。

2.3 NWL

NWL即加載到第一壁表面(W-PFC)單位面積上的平均聚變中子功率。之所以要計(jì)算NWL，是因在托卡馬克裝置中第一壁是面向等離子體的部件，在包層中起著重要作用，它的安全直接影響包層模塊的穩(wěn)定性，其計(jì)算方法如式(13)所示，其中A為表面積。

(13)

表3列出利用cosRMC計(jì)算的CFETR第一壁表面的NWL分布，各包層模塊的結(jié)果與MCNP符合良好，相對偏差均在1%以內(nèi)。進(jìn)一步計(jì)算可得到第一壁表面的平均NWL為0.305 1 MW/m2，NWL峰值出現(xiàn)在位于赤道平面的BLK3包層。

表2 cosRMC和MCNP計(jì)算的CFETR增殖包層各模塊TBR的比較Table 2 Comparison of TBR in each module of CFETR blanket

表3 CFETR第一壁表面的NWL分布計(jì)算結(jié)果比較Table 3 Comparison of NWL in each module of CFETR first wall surface

2.4 核熱沉積

托卡馬克裝置中用于約束等離子體的磁鐵結(jié)構(gòu)通常是由帶有不銹鋼包層的銅制成的，以提供抵抗變形的強(qiáng)度。銅線圈具有極好的電導(dǎo)率，允許通過的電流強(qiáng)度高達(dá)數(shù)十MA。但同時(shí)會(huì)在銅線圈上產(chǎn)生較大的熱沉積。限制等離子體放電的因素之一即是當(dāng)高強(qiáng)度的電流通過磁鐵后，磁鐵溫度的升高，應(yīng)使核熱沉積最小化。因此，研究磁鐵線圈系統(tǒng)上產(chǎn)生的核熱沉積是十分必要的，其計(jì)算方法如式(14)所示。

(14)

其中：H(E)為熱數(shù)；σt(E)為總截面；ρ為核數(shù)密度；m為柵元質(zhì)量。

表4列出利用cosRMC程序計(jì)算的CFETR的BLK1包層內(nèi)各部件的核熱沉積，兩程序模擬計(jì)算結(jié)果吻合良好，相對偏差均小于5%。同一包層上不同部件的核熱沉積不同是因?yàn)楦鞑考幬恢貌煌?，中子通量密度不同，受到的中子輻照?qiáng)度也不同。

3 結(jié)論

本文研究及推導(dǎo)了蒙特卡羅復(fù)雜曲面建模的數(shù)學(xué)模型和計(jì)算方法，并在自主化蒙特卡羅程序cosRMC原有的幾何模塊上拓展了復(fù)雜曲面建模功能。通過構(gòu)建聚變堆概念設(shè)計(jì)PPCS模型并進(jìn)行計(jì)算，柵元通量與MCNP的相對偏差小于10%，說明了該功能的正確性。

表4 CFETR的BLK1包層各部件核熱沉積計(jì)算結(jié)果比較Table 4 Comparison of nuclear thermal deposition in different components of BLK1 blanket

進(jìn)一步利用cosRMC對CFETR進(jìn)行三維精細(xì)化建模，并對包層模塊設(shè)計(jì)中的3個(gè)重要中子學(xué)參數(shù)進(jìn)行計(jì)算分析。結(jié)果表明，cosRMC與MCNP計(jì)算的各增殖包層TBR的相對偏差均小于1%。且增殖包層總的TBR為1.2，能滿足氘氚聚變堆運(yùn)行時(shí)的氚自持需求。cosRMC與MCNP計(jì)算的NWL相對偏差均在1%以內(nèi)，NWL最大值出現(xiàn)在位于赤道平面的BLK3包層。cosRMC與MCNP計(jì)算的BLK1包層內(nèi)各部件的核熱沉積的相對偏差均小于5%。研究結(jié)果證明了cosRMC應(yīng)用于聚變堆包層中子學(xué)分析的正確性和有效性。CFETR中子學(xué)參數(shù)的計(jì)算分析，也為其設(shè)計(jì)和優(yōu)化提供了參考。

感謝清華大學(xué)工程物理系核能科學(xué)與工程管理研究所反應(yīng)堆工程計(jì)算分析實(shí)驗(yàn)室(REAL團(tuán)隊(duì))在程序使用過程中提供的技術(shù)支持和幫助。