戴 濤，曹良志，賀清明，吳宏春

(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

為彌補國際熱核聚變實驗堆(ITER)與聚變示范堆(DEMO)之間的技術(shù)差距，我國提出了中國聚變工程試驗堆(CFETR)的概念[1]。在CFETR中，聚變包層承擔(dān)氚增殖、能量轉(zhuǎn)換和輻射屏蔽的作用，是最為重要的核部件，聚變包層的設(shè)計直接影響聚變堆能否安全可靠地運行。

包層的設(shè)計涉及多個學(xué)科，其中中子學(xué)和熱工水力學(xué)[2-5]具有重要地位。早期的包層設(shè)計大多采用單物理場獨立設(shè)計的方法，未考慮不同物理場間的耦合效應(yīng)。隨著設(shè)計工作的不斷深入，高精度的多物理場集成設(shè)計方法成為了包層設(shè)計研究工作的熱點。Utoh等[6]開發(fā)了二維核熱耦合程序DOHEAT。DOHEAT程序采用DOT3.5作為二維中子輸運模塊，APPLE-3作為中子學(xué)計算模塊，并編寫了二維導(dǎo)熱程序。DOHEAT程序應(yīng)用在多種日本包層設(shè)計工作中。Spagnuolo等[7]基于MCNP和ANSYS開發(fā)了中子學(xué)、熱工水力學(xué)和熱力學(xué)包層設(shè)計程序，并應(yīng)用于歐洲水冷鋰鉛包層和氦冷球床包層中。文獻[8-9]基于MCNP和ANSYS程序，分別為CFETR水冷固態(tài)包層和氦冷固態(tài)包層開發(fā)了耦合優(yōu)化設(shè)計程序NTOC和NTCOC，其主要關(guān)注單獨的中子學(xué)或熱工水力學(xué)的計算結(jié)果，因此采用了單向耦合的方式，即在中子學(xué)程序計算完成后將核熱分布結(jié)果傳遞給熱工水力學(xué)程序進行計算，未考慮熱工水力學(xué)對中子學(xué)的反饋作用。在多種裂變堆中，已證明核熱耦合效應(yīng)會顯著影響中子學(xué)和熱工水力學(xué)計算結(jié)果，包括有效中子增殖因數(shù)和功率分布等。然而，聚變包層中核熱耦合效應(yīng)帶來的影響尚不清楚，因此本文針對CFETR的兩種主要包層方案，對它們在不同功率下的核熱耦合效應(yīng)進行量化研究。

1 耦合方法和包層計算模型

1.1 雙向耦合方法

核熱耦合程序由中子學(xué)程序和熱工水力學(xué)程序構(gòu)成。為了處理復(fù)雜的聚變包層結(jié)構(gòu)，本文采用ITER項目推薦的設(shè)計工具，基于蒙特卡羅方法的MCNP程序和基于有限體積計算流體力學(xué)(CFD)方法的FLUENT程序分別作為中子學(xué)和熱工水力學(xué)求解器。由于靶核的熱運動，中子截面會隨材料溫度的改變而改變。因此，在進行核熱耦合計算時，需為輸運程序提供多溫連續(xù)能量截面庫。然而，聚變堆專用中子學(xué)截面庫FENDL-2.1僅提供室溫下的中子截面，這便要求制作不同溫度下的中子截面庫。考慮到基于多普勒展寬理論的在線截面生成方法會花費大量的時間，本文采用偽材料方法處理多溫截面[10]。偽材料方法是一種產(chǎn)生指定溫度截面的插值方法，具體實施過程如下：提前加工溫度下限TL和溫度上限TH下的中子截面庫，對處于溫度限值之間的溫度T，通過改變同一材料內(nèi)溫度下限和溫度上限的原子密度比例以等效溫度T的中子截面。原子密度插值公式為：

(1)

其中：fT為采用溫度下限中子截面的原子比例；1-fT為采用溫度上限中子截面的原子比例。

考慮到聚變堆包層結(jié)構(gòu)復(fù)雜的特點，若直接采用三維核熱耦合計算，不僅對計算程序和模型提出了嚴苛的要求，也導(dǎo)致了計算求解和數(shù)據(jù)交換方面的高昂時間成本。在聚變包層的中子學(xué)計算中，3D-1D方法已證明是有效的[11]。需要指出的是，相較于三維計算，一維的中子學(xué)計算由于模型的簡化，忽略了包層極向上結(jié)構(gòu)材料對中子的吸收、不同位置包層的氚增殖比(TBR)貢獻不同及包層間的中子泄漏等因素，導(dǎo)致計算出的TBR較真實三維中子學(xué)模型的計算結(jié)果偏高，但這種模型的簡化不會影響本文關(guān)于包層核熱耦合效應(yīng)的主要結(jié)論。

考慮到包層結(jié)構(gòu)在極向上具有周期重復(fù)特征，熱工水力學(xué)計算可采用二維模型計算以節(jié)省計算時間。因此，本文采用3D-1D-2D混合耦合方法進行聚變包層核熱耦合計算，具體計算流程圖如圖1所示。

圖1 雙向核熱耦合流程圖Fig.1 Flow chart of two-way neutronics/thermal-hydraulics coupling method

具體步驟如下：

1) 建立精細的三維全包層中子學(xué)模型，進行中子學(xué)計算，利用面統(tǒng)計功能統(tǒng)計出所計算包層的中子壁負載，中子壁在后續(xù)迭代過程中保持恒定。

2) 按照包層的徑向材料布置，建立一維圓柱中子學(xué)模型。執(zhí)行中子學(xué)計算，統(tǒng)計得到徑向核熱分布，其中實際功率密度為：

(2)

其中：p為功率密度，W/m3；N為中子壁負載，MW/m2；A為第一壁面積，m2；T為由MCNP統(tǒng)計得到的柵元核熱沉積，eV；E為聚變中子能量，MeV；V為材料區(qū)體積，m3。

3) 建立包含徑向材料層和冷卻劑通道的二維熱工水力學(xué)模型。該模型的材料結(jié)構(gòu)分布與一維中子學(xué)模型的徑向分布保持一致。再根據(jù)式(2)計算得到每個材料區(qū)域的功率密度后進行CFD計算。計算完成后，將該次計算得到的TBR和溫度分布與上次計算結(jié)果進行比較，如小于誤差限值則認為收斂，否則根據(jù)材料溫度和冷卻劑密度更新MCNP輸入卡，重新進入步驟2，迭代直至收斂。

1.2 200 MW氦冷包層模型

200 MW氦冷包層是由中國科學(xué)技術(shù)大學(xué)為CFETR提出的包層方案[12]。典型200 MW氦冷包層模塊結(jié)構(gòu)示意圖如圖2所示。該包層的前后環(huán)向?qū)挾确謩e為1 448和1 606 mm，其徑向和極向厚度分別為800和960 mm。一個典型的包層模塊由鎢涂層、第一壁(FW)、上下蓋板、支撐板、冷卻板(CP)、增殖劑區(qū)(BZ)、中子倍增劑區(qū)(MZ)、聯(lián)箱和背板(BP)組成。第一壁是1塊U型板，在面向等離子體表面覆蓋2 mm厚的鎢涂層，以保護結(jié)構(gòu)材料免受等離子體的損害。同時，在第一壁中均勻布置了徑向-環(huán)向-徑向方向的冷卻劑通道，以帶走直接來自于等離子體的熱通量。第一壁、上下蓋板和背板組成封閉的包層箱體，且在極向上均勻布置7塊支撐板以增強包層的力學(xué)性能。支撐板將包層增殖區(qū)在徑向方向上分為8個增殖單元，每個子模塊的高度為106 mm。在增殖單元中，中子倍增劑和氚增殖劑沿徑向方向均勻交錯布置以提高TBR。氦冷包層材料的選擇列于表1。由文獻[13]可知，200 MW下赤道處的氦冷包層的中子壁負載為0.47 MW/m2。

表1 氦冷包層材料的選擇Table 1 Material selection of helium cooled blanket

圖2 典型200 MW氦冷包層模塊結(jié)構(gòu)示意圖[12]Fig.2 Schematic of typical helium cooled blanket module at 200 MW[12]

1.3 200 MW水冷包層模型

200 MW水冷包層的概念設(shè)計是由中國科學(xué)院等離子體物理研究所基于成熟的壓水堆技術(shù)提出的[14]。圖3為典型200 MW水冷包層模塊結(jié)構(gòu)示意圖。包層主要由第一壁、鈹球床、混合球床、冷卻板、支撐板、聯(lián)箱和背板組成，其整體尺寸為極向1 199 mm、徑向802 mm、環(huán)向950 mm。其中，第一壁是1塊沿徑向-極向-徑向的U型板，在面向等離子體表面覆蓋2 mm的鎢涂層，且在第一壁內(nèi)布置42個截面為8 mm×8 mm的冷卻劑通道。為避免氚增殖劑與水發(fā)生化學(xué)反應(yīng)，采用Li2TiO3和Be12Ti的混合球床作為氚增殖劑。為提高TBR，在氚增殖區(qū)前方加入兩層中子倍增材料。冷卻劑采用15.5 MPa的輕水，進出口溫度分別為285 ℃和325 ℃。水冷包層的主要材料列于表2。由文獻[15]可知，200 MW時赤道處水冷包層的中子壁負載為0.454 MW/m2。

圖3 典型200 MW水冷包層模塊結(jié)構(gòu)示意圖[14]Fig.3 Schematic of typical water cooled blanket module at 200 MW[14]

表2 水冷包層材料的選擇Table 2 Material selection of water cooled blanket

2 中子學(xué)溫度敏感性分析

在進行核熱耦合計算前，本文基于200 MW的氦冷和水冷包層模型，對聚變包層主要材料進行中子學(xué)的溫度敏感性分析。FENDL-2.1中子庫中未提供能量低于4 eV的熱散射截面數(shù)據(jù)S(α，β)，但忽略熱散射修正會導(dǎo)致TBR計算結(jié)果不準確，因為鋰在低能量段的產(chǎn)氚截面很大，熱中子對TBR有較大的貢獻。因此，本文選取ENDF/B-Ⅶ.0的熱散射數(shù)據(jù)進行補充。

2.1 氦冷包層中子學(xué)溫度敏感性分析

表3、4分別為未考慮熱散射和考慮熱散射修正的氦冷包層中子倍增劑和氚增殖劑的溫度敏感性計算結(jié)果。比較表3、4中20 ℃的結(jié)果可知，考慮中子熱散射修正后，計算得到的TBR更大，核熱沉積更小。未考慮熱散射修正時，材料溫度對TBR和核熱沉積的影響很小。這表明在氦冷包層中，溫度效應(yīng)主要影響中子的熱散射行為。

表3 不考慮熱散射修正的氦冷包層中子學(xué)溫度敏感性計算結(jié)果Table 3 Neutronics result of helium cooled blanket without considering S(α, β) data at different material temperatures

由表4可知，氦冷包層的TBR隨金屬鈹溫度的增加而增大。這是因為隨溫度升高，鈹原子的熱運動更加劇烈，熱散射截面增大，產(chǎn)生更多的熱中子。此外，總核熱沉積隨鈹溫度的升高而減小，且變化速率逐漸減小，其原因是增殖區(qū)熱沉積的增加量小于中子倍增劑區(qū)熱沉積的減少量。當(dāng)Li4SiO4溫度從20 ℃上升到327 ℃時，包層的TBR從1.577 72上升至1.578 20，總核熱沉積從18.963 61 eV上升至18.968 74 eV。此后，在Li4SiO4溫度從327 ℃繼續(xù)上升的過程中，包層的TBR和總核熱沉積基本保持不變。

表4 考慮熱散射修正的氦冷包層中子學(xué)溫度敏感性計算結(jié)果Table 4 Neutronics result of helium cooled blanket with considering S(α, β) data at different material temperatures

2.2 水冷包層中子學(xué)溫度敏感性分析

與氦氣不同，輕水具有優(yōu)良的中子慢化能力，同時會吸收中子，進而影響中子能譜。考慮到低能區(qū)鋰的氚產(chǎn)生截面較大，水密度和溫度的變化對TBR有著明顯的影響。為研究水冷包層中不同材料的溫度效應(yīng)，對氚增殖劑、中子倍增劑和冷卻劑進行了溫度敏感性分析。所有材料的溫度首先設(shè)置為20 ℃，水的密度由溫度決定。由于水在15.5 MPa壓力下的飽和溫度約為345 °C，因此僅計算冷卻劑溫度20 ℃和327 ℃的情況。2.1節(jié)的計算結(jié)果證實了中子熱散射修正的必要性，因此本節(jié)計算均考慮熱散射修正(表5)。比較表4、5可知，水冷包層的中子倍增劑和氚增殖劑的溫度效應(yīng)與氦冷包層相似。當(dāng)中子倍增劑溫度升高時，TBR增大，而核熱沉積減小。由于水冷包層中子倍增劑的體積份額較小，其溫度效應(yīng)相對微弱。在20～327 ℃之間，TBR和核熱沉積隨氚增殖劑溫度的升高而增大。當(dāng)氚增殖劑的溫度在327～927 ℃之間變化時，TBR和總核熱沉積保持穩(wěn)定。值得注意的是，當(dāng)冷卻劑溫度從20 ℃增加到327 ℃時，TBR明顯增加，主要原因是水密度從1.0 g/cm3降至0.6 g/cm3。盡管水的中子慢化能力因冷卻劑密度低而較弱，但更多的中子從冷卻劑的吸收反應(yīng)中逃出，并進入增殖區(qū)發(fā)生產(chǎn)氚反應(yīng)。因此，降低水密度可有效提高TBR。

表5 考慮熱散射修正的水冷包層中子學(xué)溫度敏感性計算結(jié)果Table 5 Neutronics result of water cooled blanket with considering S(α, β) data at different material temperatures

3 聚變包層核熱耦合結(jié)果

CFETR采用工程堆和示范堆的一體化設(shè)計，分兩個階段運行：第一階段的目的是穩(wěn)態(tài)運行和實現(xiàn)氚自持，采用200 MW的低聚變功率運行；第二階段的目的是進行包括電力輸出的工程示范驗證，采用最高1.5 GW的功率運行。為了探究不同功率條件下聚變包層的核熱耦合效應(yīng)，本文對200 MW和1.5 GW功率下的CFETR固態(tài)包層設(shè)計方案進行調(diào)研，發(fā)現(xiàn)1.5 GW功率下的氦冷固態(tài)包層具有明顯的三維流動結(jié)構(gòu)[15-16]，若采用二維的熱工水力模型會導(dǎo)致計算的溫度分布顯著偏高，因此3D-1D-2D的核熱耦合計算方法不適用于該模型。因此，本文選取200 MW的氦冷固態(tài)包層和水冷固態(tài)包層及1.5 GW的水冷固態(tài)包層方案進行核熱耦合效應(yīng)的研究。

3.1 200 MW氦冷包層核熱耦合結(jié)果

耦合計算從均勻溫度分布為20 °C的中子學(xué)計算開始，表6列出了所有迭代中不同氚增殖區(qū)中TBR、相對百分比標準差。單向耦合的計算結(jié)果等于耦合計算中的第0次迭代。迭代在第2輪便達到收斂。雙向核熱耦合計算得到的TBR比單向核熱耦合計算小0.11%。由于氦氣的中子透明性，氦密度的變化幾乎不會影響中子學(xué)的結(jié)果。結(jié)合材料溫度敏感性分析，雙向核熱耦合效應(yīng)主要由金屬鈹在氦冷包層中的熱散射效應(yīng)引起。然而在高溫下，由于多普勒效應(yīng)，結(jié)構(gòu)材料會吸收更多的中子。耦合計算結(jié)果表明，核熱耦合效應(yīng)對氦冷包層的TBR影響不大。

表6 200 MW氦冷包層核熱耦合過程中TBR變化Table 6 TBR variation in neutronics/thermal-hydraulics coupling calculation of helium cooled blanket at 200 MW

圖4為單向和雙向耦合計算得到的氦冷包層溫度分布，圖5為單向和雙向耦合計算得到的氦冷包層極向上中間位置的徑向溫度分布?？煽闯觯瑴囟戎饕貜较蜃兓?，在極向上變化較小。由于較高的功率密度和較低的熱導(dǎo)率，氚增殖區(qū)溫度明顯偏高，包層的最高溫度出現(xiàn)在第2個氚增殖區(qū)。為帶走熱量，冷卻板布置在氚增殖區(qū)和中子倍增劑區(qū)之間。因此，氚增殖區(qū)和中子倍增劑區(qū)的溫度曲線在徑向上呈余弦分布，局部最高溫度出現(xiàn)在中心處。雖然在鎢裝甲和第一壁中有高熱流和高核熱沉積，但充足的冷卻能力確?？杉皶r帶走核熱沉積，從而將溫度保持在較低水平。

a——單向核熱耦合結(jié)果；b——雙向核熱耦合結(jié)果圖4 200 MW氦冷包層溫度分布Fig.4 Temperature distribution in helium cooled blanket at 200 MW

與單向核熱耦合計算相比，雙向核熱耦合計算得到的溫度分布在氚增殖區(qū)內(nèi)增大，在中子倍增劑區(qū)內(nèi)減小，總核熱沉積變化不大。單向和雙向核熱耦合計算處的溫度分布差異很小，因為在大多數(shù)區(qū)域功率密度的差異很小。由此可知，核熱耦合效應(yīng)對氦冷包層的溫度分布影響較小。

3.2 200 MW水冷包層核熱耦合結(jié)果

表7列出了200 MW水冷包層在核熱耦合迭代過程中不同氚增殖區(qū)中TBR、相對百分比標準差，由表7可看出，計算在兩次迭代后達到收斂，包層后部的TBR變化較大。包層耦合后的總TBR相較于未耦合的結(jié)果增大約4.45%。圖6為單向和雙向核熱耦合計算得到的水冷包層的溫度分布，圖7為極向中間位置處的徑向溫度分布。單向和雙向耦合計算得到的溫度分布相似：溫度沿徑向變化較大，局部最高溫度出現(xiàn)在氚增殖區(qū)中心。在單向核熱耦合計算中，由于最大功率密度，最高溫度出現(xiàn)在第1個增殖劑區(qū)中心，為960 ℃。在雙向核熱耦合計算中，最高溫度出現(xiàn)在第1個增殖劑區(qū)中心，為973 ℃，與單向核熱耦合計算相比，在雙向核熱耦合計算結(jié)果中，4個氚增殖區(qū)的最高溫度增加分別為13、22、29和26 ℃。忽略雙向核熱耦合效應(yīng)可能導(dǎo)致溫度超過溫度極限，從反應(yīng)堆安全角度出發(fā)，在水冷包層設(shè)計中應(yīng)該考慮雙向核熱耦合效應(yīng)。

a——單向核熱耦合；b——雙向核熱耦合圖6 200 MW水冷包層溫度分布Fig.6 Temperature distribution in water cooled blanket at 200 MW

圖7 200 MW水冷包層徑向溫度分布Fig.7 Temperature distribution along radial direction in water cooled blanket at 200 MW

3.3 1.5 GW水冷包層核熱耦合結(jié)果

1.5 GW的水冷包層設(shè)計方案是中國科學(xué)院等離子體物理研究所為了滿足CFETR在第二階段的設(shè)計需求于2019年更新的[17]，典型的包層結(jié)構(gòu)如圖8所示。與200 MW的設(shè)計方案相比，1.5 GW的水冷包層在環(huán)向上分為4個結(jié)構(gòu)類似的子模塊，且采用3組獨立的冷卻系統(tǒng)。此外，包層取消了單獨的中子倍增劑區(qū)，將冷卻板更改為多組冷卻管系統(tǒng)并嵌入至混合球床增殖劑區(qū)。三維中子學(xué)計算結(jié)果表明，赤道處包層的中子壁負載為1.819 MW/m2，詳細的結(jié)構(gòu)說明可參閱文獻[17-19]。

圖8 典型1.5 GW水冷包層模塊結(jié)構(gòu)示意圖[19]Fig.8 Schematic of typical 1.5 GW water cooled blanket module[19]

1.5 GW典型水冷包層核熱耦合前后的平均TBR分別為1.280 470和1.328 621，提高了約3.76%。相較于200 MW的計算結(jié)果，1.5 GW水冷包層核熱耦合計算后的TBR上升有所減小，這可能是因為1.5 GW的水冷包層取消了金屬鈹球床中子倍增劑區(qū)所導(dǎo)致的?？紤]到1.5 GW水冷包層在環(huán)向上為周期結(jié)構(gòu)，因此對單個子模塊建模進行CFD模擬。圖9為耦合前后的1.5 GW水冷包層的溫度分布，可看出，經(jīng)過核熱耦合計算后，后方的增殖區(qū)溫度升高，這是因為在核熱耦合計算后，冷卻劑的密度分布隨徑向方向逐漸減小，更多的中子進入了后方的增殖區(qū)并產(chǎn)生能量沉積導(dǎo)致后方增殖區(qū)功率密度升高。

a——單向核熱耦合；b——雙向核熱耦合圖9 1.5 GW水冷包層溫度分布 Fig.9 Temperature distribution in water cooled blanket at 1.5 GW

4 結(jié)論

本文研究了CFETR包層雙向核熱耦合效應(yīng)對中子學(xué)和熱工水力學(xué)計算結(jié)果的影響?；贛CNP和FLUENT程序，利用3D-1D-2D混合方法和偽材料方法實現(xiàn)和加速核熱耦合計算，并對包層主要材料的中子學(xué)溫度敏感性進行了分析，闡明了影響雙向耦合核熱效應(yīng)的主要因素。對200 MW的氦冷固態(tài)包層、200 MW和1.5 GW水冷固態(tài)包層進行了雙向核熱耦合計算，主要結(jié)論如下。

1) 在進行CFETR包層的核熱耦合時，應(yīng)考慮中子的熱散射修正。因為高溫下的鈹熱散射截面較大，會產(chǎn)生更多的熱中子，鈹溫度的提高有利于提高包層的氚增殖比。

2) 在氦冷包層中，核熱耦合效應(yīng)主要受到金屬鈹?shù)臒嵘⑸湫袨橛绊?，核熱耦合效?yīng)對200 MW功率下氦冷包層的作用可忽略。

3) 在水冷包層中，TBR和核熱沉積對水的溫度較敏感，因為溫度變化引起的水密度變化會顯著影響冷卻劑的中子吸收和慢化能力。核熱耦合效應(yīng)將使200 MW功率下水冷包層的TBR提高約4.45%，1.5 GW功率下的水冷包層TBR提高約3.76%，并略微提高包層后方增殖區(qū)的溫度。從準確的包層設(shè)計和安全性分析的角度來看，對聚變水冷包層進行核熱耦合計算具有一定價值。

感謝中國科學(xué)院等離子體物理研究所劉松林研究員和核工業(yè)西南物理研究院武興華副研究員與曹啟祥副研究員提供的聚變包層模型。