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

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

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

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

1 二維/一維耦合方法

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

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

(1)

式中：μ為極角；αm為方位角；Σt,g為總截面；φg,m表示第g群m方向上的角通；Qg,m為總源項，包括各向同性的散射源和裂變源。

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

(2)

(3)

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

(4)

(5)

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

圖1 二維/一維耦合輸運計算方法Fig.1 2D/1D coupling transport calculation

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

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

(6)

(7)

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

(8)

(9)

(10)

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

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

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

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

(11)

其中：

(12)

或者：

(13)

其中：

(14)

2.2 改進的泄漏項分割技術(shù)

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

定義泄漏項：

(15)

(16)

定義極角相同，方位角相反的2個方向為正方向和反方向，此時2個方向的泄漏項為：

(17)

(18)

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

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

1)對于正方向，當MOC計算的出射角通量小于0時，進行分割(忽略位置和角度標識)：

(19)

(20)

2)對于反方向，首先在負方向總源項上加上反向修正項，此時正反2個方向的泄漏項相同，相當于對正反2個方向的泄漏項做各向同性處理。

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

(21)

(22)

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

3 精度和收斂性分析

3.1 計算精度分析

1 強泄漏問題

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

圖2 強泄漏算例幾何模型Fig.2 Geometry of high leakage case

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

2 三維C5G7基準題

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

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

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

3.2 收斂性分析

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

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

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

4 結(jié)論

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

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

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