孔勃然， 朱凱杰， 劉保坤， 張漢， 郭炯， 郝琛， 侯捷， 李富

(1.清華大學(xué) 核能與新能源研究院，北京 100084； 2.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院，黑龍江 哈爾濱 150001)

二維/一維耦合方法是一種高保真的一步法計(jì)算方法，它將三維方程分解為徑向二維方程和軸向一維方程，并通過泄漏項(xiàng)將二者進(jìn)行耦合求解。受限于計(jì)算速度和內(nèi)存，直接三維全堆芯一步法計(jì)算尚存難度[1]，具有相同精度并節(jié)省計(jì)算量和內(nèi)存需求的二維/一維耦合方法目前更具實(shí)際的工程價(jià)值和研究價(jià)值，是近年的研究熱點(diǎn)[2]。目前采用二維/一維方法的典型程序有DeCART[3]、MPACT[4]、NECPX[5]，HNET[6]、KuaFu[7]等，這些程序普遍采用特征線方法(method of characteristics，MOC)用于二維計(jì)算，一維軸向方程則采用擴(kuò)散方法如節(jié)塊法(nodal expansion method，NEM)和SP3方法或者輸運(yùn)方法如離散縱標(biāo)法，泄漏項(xiàng)則采用各向同性或各向異性泄漏項(xiàng)。數(shù)值實(shí)驗(yàn)表明，采用軸向SN計(jì)算并配合各向異性泄漏項(xiàng)的方法精度更高，但與此同時(shí)計(jì)算量提高、收斂性變差[8]。為了解決計(jì)算量大幅提升的問題，MPACT[9]提出方位角傅里葉級(jí)數(shù)展開的各向異性泄漏項(xiàng)代替完全各向異性的泄漏項(xiàng)，并證明1階或者2階傅里葉級(jí)數(shù)展開就可達(dá)到完全各向異性的精度。由于采用了各向異性泄漏項(xiàng)，方程的收斂性相較于采用各向同性泄漏項(xiàng)的方法變差。在迭代過程中，由于泄漏項(xiàng)的存在，二維輸運(yùn)計(jì)算過程中可能會(huì)計(jì)算出負(fù)角通量，并有可能導(dǎo)致方程發(fā)散。為了解決這一問題，MPACT提出泄漏項(xiàng)分割技術(shù)，但是該方法存在計(jì)算誤差較大的問題[10]。

為了保證精度的同時(shí)降低計(jì)算負(fù)擔(dān)，本文采用方位角傅里葉級(jí)數(shù)展開的各向異性泄漏項(xiàng)和軸向SN的組合的二維/一維計(jì)算方法，并提出了一套適用于該組合的改進(jìn)的泄漏項(xiàng)分割技術(shù)。該方法不會(huì)增加內(nèi)存負(fù)擔(dān)的同時(shí)，相比于傳統(tǒng)方法具有更好的效率和收斂性。

1 二維/一維耦合方法

1.1 二維/一維耦合方法的理論

直角坐標(biāo)系下，穩(wěn)態(tài)形式離散的三維玻爾茲曼輸運(yùn)方程形式為：

(1)

式中：μ為極角；αm為方位角；Σt,g為總截面；φg,m表示第g群m方向上的角通；Qg,m為總源項(xiàng)，包括各向同性的散射源和裂變?cè)础?/p>

二維/一維方法將堆芯沿軸向分為若干層，對(duì)方程(1)沿軸向積分取平均即可獲得徑向二維方程：

(2)

(3)

同理，二維/一維方法將堆芯沿徑向分為若干個(gè)平源區(qū)/柵元，對(duì)方程(1)沿徑向積分取平均即可獲得軸向一維方程：

(4)

(5)

從式(2)和(4)可以看出，二維/一維方程通過泄漏項(xiàng)進(jìn)行耦合，在計(jì)算前泄漏項(xiàng)需要從一維/二維方程中獲得。通過反復(fù)迭代計(jì)算二維/一維方程直至特征值和標(biāo)通量收斂，其原理見圖1。

圖1 二維/一維耦合輸運(yùn)計(jì)算方法Fig.1 2D/1D coupling transport calculation

1.2 方位角傅里葉級(jí)數(shù)展開的各向異性泄漏項(xiàng)及軸向一維方程

從式(5)知，若采用各向異性泄漏項(xiàng)，雖然軸向方程在空間維度是一維的，但它仍然與方位角有關(guān)，此時(shí)軸向計(jì)算仍需計(jì)算所有角度。若軸向方程采用輸運(yùn)方法求解，其計(jì)算量相比于各向同性泄漏項(xiàng)顯著增加，以KUCA和OECD C5G7-3D RB基準(zhǔn)題為例，軸向計(jì)算量約占總計(jì)算量的20%，可見此時(shí)軸向計(jì)算非常耗時(shí)。為了解決軸向計(jì)算量大幅提升的問題，采用方位角傅里葉級(jí)數(shù)展開的泄漏項(xiàng)代替完全各向異性泄漏項(xiàng)，以實(shí)現(xiàn)用少量階數(shù)的軸向方程代替全角度的軸向方程，將角通量和徑向泄漏項(xiàng)采用N階傅里葉級(jí)數(shù)展開為：

(6)

(7)

將式(6)和(7)代入方程(4)中計(jì)算后，對(duì)方位角進(jìn)行積分即可得到零階、sin階、cos階方程：

(8)

(9)

(10)

至此，通過方位角傅里葉級(jí)數(shù)展開泄漏項(xiàng)，可以將原有全方位角的SN計(jì)算轉(zhuǎn)化成 (2N+1)個(gè)SN方程進(jìn)行計(jì)算。

2 泄漏項(xiàng)分割技術(shù)

2.1 傳統(tǒng)的泄漏項(xiàng)分割技術(shù)

由式(2)可知，當(dāng)泄漏項(xiàng)過大時(shí)右端源項(xiàng)為負(fù)，可能導(dǎo)致輸運(yùn)計(jì)算的角通量為負(fù)，進(jìn)而導(dǎo)致算法不收斂。為了解決這一問題，需要采用泄漏項(xiàng)分割技術(shù)進(jìn)行修正，其原理為當(dāng)右端源項(xiàng)為負(fù)值時(shí)，修改方程(2)為：

(11)

其中：

(12)

或者：

(13)

其中：

(14)

2.2 改進(jìn)的泄漏項(xiàng)分割技術(shù)

負(fù)源項(xiàng)在二維/一維算法中是很難避免的，從算法的機(jī)理而言，即使收斂時(shí)負(fù)源項(xiàng)也可能存在。一定程度的負(fù)源項(xiàng)在計(jì)算過程中是可以接受的，而輸運(yùn)計(jì)算中負(fù)角通量是無意義的，因此新的泄漏項(xiàng)分割技術(shù)的觸發(fā)條件為當(dāng)輸運(yùn)計(jì)算出的角通量為負(fù)時(shí)進(jìn)行泄漏項(xiàng)分割。

定義泄漏項(xiàng)：

(15)

(16)

定義極角相同，方位角相反的2個(gè)方向?yàn)檎较蚝头捶较?，此時(shí)2個(gè)方向的泄漏項(xiàng)為：

(17)

(18)

從式(17)、(18)可以發(fā)現(xiàn)，正反2個(gè)方向的泄漏項(xiàng)區(qū)別只在于奇階泄漏項(xiàng)的正負(fù)，其余0階泄漏項(xiàng)和偶階泄漏項(xiàng)均相同，利用此特征設(shè)計(jì)改進(jìn)的泄漏項(xiàng)分割技術(shù)。

對(duì)于正向方程，首先將右端源項(xiàng)中奇階泄漏項(xiàng)去除，若此時(shí)源項(xiàng)大于零，則此時(shí)不對(duì)方程進(jìn)行修改；反之采用類似于傳統(tǒng)泄漏項(xiàng)修正的方法將右端源項(xiàng)置零，左端引入修正截面。具體操作如下：

1)對(duì)于正方向，當(dāng)MOC計(jì)算的出射角通量小于0時(shí)，進(jìn)行分割(忽略位置和角度標(biāo)識(shí))：

(19)

(20)

2)對(duì)于反方向，首先在負(fù)方向總源項(xiàng)上加上反向修正項(xiàng)，此時(shí)正反2個(gè)方向的泄漏項(xiàng)相同，相當(dāng)于對(duì)正反2個(gè)方向的泄漏項(xiàng)做各向同性處理。

若此時(shí)出射角通量仍然小于0，進(jìn)行和傳統(tǒng)修正方法類似分割：

(21)

(22)

該修正方法充分利用了泄漏項(xiàng)的特性，當(dāng)正方向出現(xiàn)負(fù)角通量時(shí)，通過上述處理，相當(dāng)于對(duì)正反2個(gè)方向的泄漏項(xiàng)做各向同性處理。雖然在特定方向上的新方程不等價(jià)于原方程，但對(duì)所有方向新修正的方程進(jìn)行積分，新方程積分等價(jià)于原方程積分，這就保證了二維計(jì)算和一維計(jì)算得到的標(biāo)通量一致等價(jià)，因此該修正方法改善了精度和收斂性。

3 精度和收斂性分析

3.1 計(jì)算精度分析

1 強(qiáng)泄漏問題

該算例的幾何模型如圖2所示，徑向上燃料棒3×3布置，最中心為慢化劑，徑向上采用3個(gè)反射邊界條件和1個(gè)真空邊界條件，軸向上采用反射邊界條件。軸向上布置2層燃料層和1層慢化劑層?？梢詫⒋怂憷愅茷閺较?×7、17×17算例，以研究不同泄漏強(qiáng)度對(duì)精度的影響。該算例中以完全各向異性泄漏項(xiàng)方法的結(jié)果作為基準(zhǔn)解，分別測(cè)試不采用修正，傳統(tǒng)泄漏項(xiàng)修正和新的泄漏項(xiàng)修正方法，其結(jié)果如表1所示。數(shù)值實(shí)驗(yàn)證明采用改進(jìn)的泄漏項(xiàng)修正方法可以降低修正誤差，其結(jié)果與不修正方法近似且解決了收斂性的問題。采用修正方法、傳統(tǒng)修正、改進(jìn)修正方法時(shí)，三者的迭代次數(shù)完全相同，而單次迭代時(shí)間也幾乎相同，因此此算例中，改進(jìn)泄漏項(xiàng)方法和傳統(tǒng)修正方法的計(jì)算時(shí)間幾乎一致。

圖2 強(qiáng)泄漏算例幾何模型Fig.2 Geometry of high leakage case

表1 強(qiáng)泄漏算例計(jì)算結(jié)果Table 1 Calculation results of high leakage cases

2 三維C5G7基準(zhǔn)題

C5G7基準(zhǔn)題是OECD/NEA發(fā)布的帶控制棒的非均勻輸運(yùn)計(jì)算基準(zhǔn)題[11]，分為無棒、半插棒和全插棒。參數(shù)選取為全角度64個(gè)方位角，8個(gè)極角，特征線密度0.01 cm。徑向劃分為51個(gè)柵元，每個(gè)柵元?jiǎng)澐譃?個(gè)圓環(huán)16個(gè)扇區(qū)。軸向分為18層，每層3.57 cm。表2為G5G7計(jì)算結(jié)果，結(jié)果顯示采用方位角傅里葉級(jí)數(shù)展開的泄漏項(xiàng)并采用改進(jìn)的泄漏項(xiàng)分割技術(shù)的結(jié)果可以達(dá)到較好的結(jié)果，1階或者2階的結(jié)果就可以和各向異性的結(jié)果相當(dāng)。圖3為C5G73D-Unrodded功率誤差，可見采用傅里葉級(jí)數(shù)展開的泄漏項(xiàng)和各向異性泄漏項(xiàng)的結(jié)果幾乎一致。

圖3 C5G7-Unrodded功率誤差Fig.3 Power error of C5G7-Unrodded cases

表2 C5G7計(jì)算結(jié)果Table 2 Results of C5G7 cases

3.2 收斂性分析

為了分析采用改進(jìn)的泄漏項(xiàng)修正方法的2D/1D耦合算法的收斂性，選用KUCA基準(zhǔn)題[12]進(jìn)行研究。KUCA是日本東京大學(xué)模擬小型反應(yīng)堆的臨界裝置，其堆芯包括燃料、慢化劑、控制棒。選取其中KUCA空隙問題(情況1)作為算例進(jìn)行驗(yàn)證，基準(zhǔn)解為0.977 8，其結(jié)果如表3所示。從計(jì)算結(jié)果中可以看出，不采用修正方法時(shí)，1階和3階傅里葉級(jí)數(shù)展開方法是發(fā)散的，其原因?yàn)樨?fù)角通量過多，導(dǎo)致了負(fù)通量的產(chǎn)生，進(jìn)而導(dǎo)致算法發(fā)散。而傳統(tǒng)修正方法1階傅里葉級(jí)數(shù)展開時(shí)仍然發(fā)散，其原因?yàn)閭鹘y(tǒng)方法的所有角度的修正方程積分不等價(jià)于原方程，當(dāng)負(fù)角通量過多時(shí)，修正后的方程和原方程偏差過大，二維計(jì)算統(tǒng)計(jì)出的標(biāo)通量不等價(jià)于一維計(jì)算統(tǒng)計(jì)出的標(biāo)通量，進(jìn)而導(dǎo)致整體算法無法收斂。結(jié)果證明改進(jìn)的泄漏項(xiàng)分割方法可以解決負(fù)角通量造成的發(fā)散問題。表4展示了不同修正方法的迭代次數(shù)對(duì)比，由于每次迭代計(jì)算時(shí)間幾乎相同，因此迭代次數(shù)即反映了計(jì)算時(shí)間的對(duì)比,從結(jié)果中可以看出幾種方法的計(jì)算時(shí)間幾乎相同。

表3 KUCA計(jì)算結(jié)果Table 3 Results of KUCA cases 10-5

表4 KUCA迭代次數(shù)對(duì)比Table 4 Iteration time of KUCA cases 10-5

4 結(jié)論

1)采用各向異性泄漏項(xiàng)可以提高二維/一維耦合算法的精度，但是同時(shí)會(huì)帶來計(jì)算量提升和收斂性變差的問題。為了解決計(jì)算量提升問題，采用方位角傅里葉級(jí)數(shù)展開的泄漏項(xiàng)近似完全各向異性泄漏項(xiàng)。數(shù)值實(shí)驗(yàn)表明1階或者2階的傅里葉展開就能達(dá)到和完全各向異性相同的精度。

2)傳統(tǒng)的泄漏項(xiàng)分割技術(shù)會(huì)損失一定精度和收斂性。本文提出了一種改進(jìn)的泄漏項(xiàng)分割技術(shù)，該方法對(duì)整體的精度損失較小的前提下解決了二維/一維算法的收斂性問題。

3)改進(jìn)的泄漏項(xiàng)分割方法接受一定程度的負(fù)源項(xiàng)并且充分利用了傅里葉級(jí)數(shù)展開的特性，盡可能的減少對(duì)原方程的修改，以便達(dá)到更好的精度和收斂性。強(qiáng)泄漏算例，C5G7三維基準(zhǔn)題算例和KUCA基準(zhǔn)題算例證明了新方法在不損失精度的同時(shí)提高了二維/一維算法的穩(wěn)定性。