馬輝強 于 濤 陳珍平 趙鵬程 謝金森 劉紫靜 劉建星雷洲陽 何清玲

1（南華大學(xué)核科學(xué)技術(shù)學(xué)院 衡陽421001）

2（南華大學(xué)湖南省數(shù)字化反應(yīng)堆工程技術(shù)研究中心 衡陽421001）

環(huán)形燃料芯塊具有優(yōu)越的經(jīng)濟型和安全性［1］，受到了國內(nèi)外的廣泛關(guān)注。美國已將環(huán)形燃料確定為新一代核電燃料進行研發(fā)，且將其作為美國現(xiàn)役核電站延壽至80年研究計劃中的首選燃料，麻省理工學(xué)院（Massachusetts Institute of Technology）針對環(huán)形燃料元件應(yīng)用壓水堆以提升堆芯功率密度做了深入研究［2-4］；韓國原子能科學(xué)研究所、中國原子能科學(xué)研究院分別對韓國標(biāo)準(zhǔn)型核電站“OPR1000”堆芯和秦山二期反應(yīng)堆堆芯在不改變原有燃料組件尺寸和水鈾比條件下，通過用環(huán)形燃料元件替代原有的棒狀燃料組件進行了堆芯物理和熱工水力分析［5-6］。

實心圓柱燃料芯塊溫度場求解方法無法獲得環(huán)形燃料芯塊溫度場的分布，差分和有限元方法計算燃料芯塊溫度場分布相對耗時，不利于大規(guī)模計算，有必要分析環(huán)形燃料芯塊溫度場計算的方法，以提高環(huán)形燃料芯塊溫度場的計算精度和計算效率。本工作針對環(huán)形燃料芯塊，從導(dǎo)熱微分方程和傅立葉定律出發(fā)，導(dǎo)出了環(huán)形燃料芯塊平均熱導(dǎo)率計算方法，建立環(huán)形燃料芯塊穩(wěn)態(tài)溫度場分布、最高溫度以及最高位溫度位置的模型，并驗證該模型的正確性。

1 環(huán)形燃料芯塊溫度場模型

1.1 徑向最高溫度位置

雙面冷卻的環(huán)形燃料芯塊的導(dǎo)熱存在熱量分配問題［6］，芯塊中的熱量分別從芯塊內(nèi)外表面中導(dǎo)出，芯塊內(nèi)外表面各自導(dǎo)出熱量的多少取決于芯塊內(nèi)側(cè)和外側(cè)的熱阻大小。實心圓柱燃料芯塊的最高溫度的位置一般在芯塊的中心軸上，最高溫度容易求解；環(huán)形燃料芯塊中的最高溫度以及最高溫度的位置是隨內(nèi)外兩側(cè)熱阻的變化而變化［7］。環(huán)形燃料芯塊中需要確定最高溫度的位置后才能計算環(huán)形芯塊的最高溫度［8］，最高溫度位置通過傅立葉定律和積分熱導(dǎo)率進行計算。環(huán)形芯塊內(nèi)外表面溫度由環(huán)形燃料元件的對流換熱計算得到，本工作采用文獻［9］中的計算得到芯塊內(nèi)外表面溫度。環(huán)形燃料芯塊模型示意圖如圖1所示。

圖1 環(huán)形燃料芯塊模型Fig.1 Model of annular fuel pellets

忽略軸向?qū)幔h(huán)形燃料芯塊中傅立葉定律為［10］：

式中：kf（t）為燃料芯塊的熱導(dǎo)率函數(shù)，W·m-1·K-1；qv為燃料芯塊的功率密度，W·m-3；rm為芯塊中絕熱面（即芯塊徑向最高溫度的位置）的半徑，m。對式（1）兩端同時積分并簡化，得到環(huán)形燃料芯塊的積分熱導(dǎo)率：

式中：ru和tu分別表示芯塊表面處的半徑和溫度。再將式（2）簡化得到積分熱導(dǎo)率的最終的形式：

假設(shè)環(huán)形燃料芯塊的最高溫度出現(xiàn)在徑向rm處，分別從燃料芯塊的內(nèi)外表面出發(fā)，從芯塊兩側(cè)往中間計算，可以得到在rm處燃料芯塊溫度的表達式，且這兩式相等，即：

然后得到芯塊中最高溫度的位置rm為：

1.2 芯塊溫度場計算

忽略軸向?qū)幔蓤A柱坐標(biāo)系一維穩(wěn)態(tài)傳導(dǎo)熱微分方程導(dǎo)出環(huán)形燃料芯塊徑向溫度分布。環(huán)形燃料芯塊導(dǎo)熱微分方程方程和邊界條件［11］為：

式中：r為燃料芯塊半徑，m；燃料芯塊的內(nèi)外表面半徑分別為rui和ruo，m；Tfui、Tfuo分別為燃料芯塊內(nèi)外表面的溫度，K；kf為燃料芯塊的平均熱導(dǎo)率，W·m-1·K-1。對于式（6）的解來說，邊界條件（8）和（9）是等價的。式（6）的通解為：

UO2等核燃料的熱導(dǎo)率隨溫度而變化，若要確定環(huán)形燃料芯塊一維穩(wěn)態(tài)溫度分布，則需要準(zhǔn)確地確定式（6）中芯塊的平均熱導(dǎo)率kf（在導(dǎo)熱介質(zhì)所處的溫度范圍內(nèi)用一個熱導(dǎo)率常數(shù)等效替換對隨溫度變化的熱導(dǎo)率）。對于穩(wěn)態(tài)情況，芯塊內(nèi)外表面溫度確定后，環(huán)形燃料芯塊絕熱面的位置則可以確定，分別從燃料芯塊的內(nèi)外表面出發(fā)，從兩側(cè)往芯塊中間計算，在絕熱面位置處式（11）中兩個函數(shù)計算結(jié)果必須相等，由此可得出絕熱面rm的位置：

聯(lián)立式（4）、（11）得到環(huán)形燃料芯塊的平均熱導(dǎo)率kf：

由式（12）可以發(fā)現(xiàn)環(huán)形燃料芯塊的平均熱導(dǎo)率與傳統(tǒng)實心芯塊的平均熱導(dǎo)率的計算方法具有一定差異，環(huán)形燃料芯塊的平均熱導(dǎo)率主要取決于芯塊內(nèi)外表面溫度。根據(jù)數(shù)學(xué)極限，當(dāng)tui=tuo時，kf=kf（tui）=kf（tuo），當(dāng)r=rm時可求得芯塊中的最高溫度tmo。

綜上，式（5）、（11）、（13）即為從導(dǎo)熱微分方程和傅立葉定律得出的環(huán)形燃料芯塊穩(wěn)態(tài)一維溫度場分布模型。

1.3 芯塊溫度場計算修正

求解環(huán)形燃料芯塊的溫度場分布的關(guān)鍵在于求解芯塊的平均熱導(dǎo)率。由于環(huán)形燃料芯塊的熱導(dǎo)率可能較低，與溫度呈非線性關(guān)系變化，且隨溫度變化幅度較大（如二氧化鈾），對于較厚的環(huán)形燃料芯塊，芯塊表面與絕熱面之間溫度梯度將變得更加陡峭，從而使熱導(dǎo)率對空間位置產(chǎn)生較強的依賴性，即熱導(dǎo)率隨燃料芯塊徑向位置的變化幅度大，最終導(dǎo)致式（13）計算的芯塊平均熱導(dǎo)率偏大，溫度越高偏差越大。針對這個問題，由式（5）、（12）確定了絕熱面位置后，根據(jù)已知的芯塊表面溫度tui、tuo和進一步得出的芯塊絕熱面位置rm以及修正前的絕熱面最高溫度tmo，然后將環(huán)形燃料芯塊從絕熱面位置分成內(nèi)外兩層圓柱環(huán)，根據(jù)式（14）分別修正內(nèi)外兩側(cè)的芯塊平均熱導(dǎo)率kfi和kfo，然后由式（15）分別求解內(nèi)外環(huán)的溫度場，最高溫度取兩側(cè)計算的平均值。

2 模型驗證

為了驗證模型的正確性，基于本文構(gòu)建的計算模型開發(fā)了環(huán)形燃料芯塊溫度場計算程序PTFA（Program of Temperature Field of Annular Fuel Pellets），計算不同功率密度、不同芯塊厚度下的絕熱面位置、最高溫度以及芯塊的溫度場分布，并與有限元程序COMSOL5.4［12］計算的結(jié)果進行對比驗證，對比兩種方法計算結(jié)果的相對偏差。

2.1 模型參數(shù)

根據(jù)現(xiàn)有文獻［13-14］，反應(yīng)堆正常工況下的最高線功率密度為61 kW·m-1，同時參考已有的文獻中環(huán)形燃料的線功率密度參數(shù)［12］，選擇環(huán)形燃料芯塊線功率密為36 kW·m-1、54 kW·m-1以及71 kW·m-1進行計算，采用的熱導(dǎo)率函數(shù)為95%理論密度二氧化鈾的熱導(dǎo)率［12］，現(xiàn)有研究成果表明：環(huán)形燃料芯塊的厚度的范圍一般不超過3 mm［2，15］。驗證模型和參數(shù)如圖2和表1、表2所示。在COMSOL5.4中建立一維環(huán)形燃料芯塊模型，芯塊上下表面設(shè)為絕熱邊界條件，確定芯塊內(nèi)外側(cè)表面溫度，并對其進行細網(wǎng)格劃分，設(shè)置好芯塊的熱源和熱導(dǎo)率函數(shù)后采用COMSOL5.4的固體導(dǎo)熱模塊求解不同厚度環(huán)形燃料芯塊的各個功率密度工況的一維穩(wěn)態(tài)溫度場。

圖2 環(huán)形燃料芯塊COMSOL5.4計算模型Fig.2 COMSOL5.4 calculation model of annular fuel pellet

表1 環(huán)形燃料芯塊模型參數(shù)Table 1 Parameters of annular fuel pellets model

表2 不同厚度環(huán)形燃料芯塊的體功率密度及表面溫度Table 2 Power density and surface temperature of annular fuel pellets with different thickness

2.2 結(jié)果分析

核反應(yīng)堆設(shè)計和運行中，燃料芯塊的溫度是一個非常重要的參數(shù)。對于環(huán)形燃料芯塊需要準(zhǔn)確地得到芯塊的溫度場，包括絕熱面位置、芯塊最高溫度等。本模型的驗證主要針對絕熱面位置、芯塊最高溫度以及芯塊中溫度場分布等進行。

2.2.1 芯塊絕熱面位置

環(huán)形燃料芯塊絕熱面的位置決定環(huán)形燃料元件內(nèi)外冷卻通道的熱量分配比，從而影響內(nèi)外通道冷卻劑流量的分配和冷卻劑的壓降，最終決定著環(huán)形燃料元件的導(dǎo)熱性能，影響環(huán)形燃料元件的安全性和經(jīng)濟性。不同功率密度下、不同芯塊厚度的環(huán)形燃料芯塊，PTFA和有限元程序COMSOL5.4計算的最高溫度位置如表3所示。結(jié)果表明：對于絕熱面的位置的計算，PTFA計算結(jié)果與COMSOL5.4計算的結(jié)果符合很好，二者之間最大偏差小于0.15%，驗證了環(huán)形燃料芯塊溫度場模型計算絕熱面位置的正確性。

2.2.2 芯塊絕熱面最高溫度

環(huán)形燃料芯塊的最高溫度是評價燃料元件性能的重要參數(shù)之一，準(zhǔn)確地獲得燃料芯塊的溫度對環(huán)形燃料元件裂變氣體釋放率、腫脹、使用壽命和燃耗計算是十分必要的。對于不同功率密度下、不同芯塊厚度的環(huán)形燃料芯塊，PTFA和有限元程序COMSOL5.4計算的絕熱面最高溫度如表4所示。結(jié)果表明：對于環(huán)形燃料芯塊最高溫度的計算，PTFA計算的結(jié)果與COMSOL5.4計算的結(jié)果符合很好。采用未修正的芯塊平均熱導(dǎo)率計算的芯塊最高溫度結(jié)果總體相對于COMSOL5.4計算的結(jié)果稍微偏小，二者的相對偏差小于1.5%，這是由于二氧化鈾熱導(dǎo)率較低，導(dǎo)致芯塊中溫度梯度較大，從而使熱導(dǎo)率對空間位置產(chǎn)生較強的依賴性，最終導(dǎo)致式（13）計算的芯塊平均熱導(dǎo)率偏大。采用修正后的芯塊平均熱導(dǎo)率時，PTFA計算的最高溫度與COMSOL5.4計算的結(jié)果相對偏差小于0.2%，顯著提高了計算精度，驗證了環(huán)形燃料芯塊溫度場模型 計算環(huán)形燃料芯塊絕熱面最高溫度的正確性。

表3 不同線功率密度下環(huán)形燃料芯塊絕熱面位置Table 3 Position of insulation surface of annular fuel pellets under different linear power densities

表4 不同線功率密度下環(huán)形燃料芯塊最高溫度Table 4 Maximum temperature of annular fuel pellets under different linear power densities

2.2.3 芯塊溫度場分布

如圖3所示，對于不同功率密度下、不同芯塊厚度的環(huán)形燃料芯塊，PTFA與有限元程序COMSOL5.4計算的溫度場分布符合較好。由于二氧化鈾熱導(dǎo)率較小引起環(huán)形燃料芯塊中溫度梯度較大，從而引起熱導(dǎo)率對空間位置較強的依賴性，溫度越高，依賴性越強，采用未修正的芯塊平均熱導(dǎo)率計算的芯塊最高溫度總體相對于COMSOL5.4計算的結(jié)果稍微偏小，且兩者間的溫差隨著芯塊溫度的增加而增加。修正后的熱導(dǎo)率將整個芯塊的平均熱導(dǎo)率以絕熱面為中心分為內(nèi)外兩個區(qū)域的平均熱導(dǎo)率，減弱了熱導(dǎo)率對空間位置的依賴性，采用修正后的芯塊平均熱導(dǎo)率計算時，PTFA計算的溫度場分布與有限元程序COMSOL5.4計算的溫度場分布符合得非常好，相對偏差很小。對于熱導(dǎo)率較小或熱導(dǎo)率隨溫度變化較大的燃料芯塊的溫度場計算，采用修正后的燃料芯塊平均熱導(dǎo)率進行計算以獲得較為精確的計算結(jié)果。

2.2.4 環(huán)形燃料芯塊厚度對溫度場的影響

如圖3所示，環(huán)形燃料芯塊厚度相同時，增大燃料芯塊的功率密度，燃料芯塊的最高溫度明顯升高，芯塊內(nèi)部的溫度梯度將變得更加陡峭。結(jié)合表1、3、4，當(dāng)燃料芯塊的線功率密度和芯塊內(nèi)外表面溫度一定時，環(huán)形燃料芯塊越厚，芯塊的體功率密度越小，芯塊的最高溫度越高；環(huán)形燃料芯塊越薄，芯塊的體功率密度越大，芯塊的最高溫度越低。所以減小芯塊的厚度比降低燃料芯塊功率密度更能有效降低環(huán)形燃料芯塊的運行溫度。

圖3 不同芯塊厚度（1 mm、2 mm、3 mm）和不同線功率密度（36 kW·m-1、54 kW·m-1、71 kW·m-1）下環(huán)形芯塊溫度場Fig.3 Temperature field of annular fuel pellets under different thickness(1 mm,2 mm,3 mm)and different linear power density(36 kW·m-1,54 kW·m-1,71 kW·m-1)

3 結(jié)語

本文根據(jù)導(dǎo)熱微分方程和傅立葉定律，針對環(huán)形燃料的芯塊建立了環(huán)形燃料芯塊一維穩(wěn)態(tài)溫度場計算模型，導(dǎo)出了環(huán)形燃料芯塊絕熱面位置以及燃料芯塊最高溫度的計算方法。并對該模型進行了驗證和分析，得出了以下結(jié)論：

1）環(huán)形燃料芯塊一維穩(wěn)態(tài)溫度場計算關(guān)鍵在于準(zhǔn)確計算環(huán)形燃料芯塊最高溫度的位置以及平均熱導(dǎo)率。利用修正后的平均熱導(dǎo)率計算的環(huán)形燃料芯塊溫度場結(jié)果與COMSOL5.4計算結(jié)果的相對誤差小于2‰，該模型可以用于環(huán)形燃料元件的熱工水力分析，具有較高的計算精度和效率。

2）減小芯塊的厚度比降低燃料芯塊功率密度更能有效降低環(huán)形燃料芯塊的運行溫度。