環(huán)形燃料芯塊一維穩(wěn)態(tài)溫度場(chǎng)計(jì)算方法研究

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

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

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

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

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

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

1.1 徑向最高溫度位置

雙面冷卻的環(huán)形燃料芯塊的導(dǎo)熱存在熱量分配問題［6］，芯塊中的熱量分別從芯塊內(nèi)外表面中導(dǎo)出，芯塊內(nèi)外表面各自導(dǎo)出熱量的多少取決于芯塊內(nèi)側(cè)和外側(cè)的熱阻大小。實(shí)心圓柱燃料芯塊的最高溫度的位置一般在芯塊的中心軸上，最高溫度容易求解；環(huán)形燃料芯塊中的最高溫度以及最高溫度的位置是隨內(nèi)外兩側(cè)熱阻的變化而變化［7］。環(huán)形燃料芯塊中需要確定最高溫度的位置后才能計(jì)算環(huán)形芯塊的最高溫度［8］，最高溫度位置通過傅立葉定律和積分熱導(dǎo)率進(jìn)行計(jì)算。環(huán)形芯塊內(nèi)外表面溫度由環(huán)形燃料元件的對(duì)流換熱計(jì)算得到，本工作采用文獻(xiàn)［9］中的計(jì)算得到芯塊內(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。對(duì)式（1）兩端同時(shí)積分并簡(jiǎn)化，得到環(huán)形燃料芯塊的積分熱導(dǎo)率：

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

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

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

1.2 芯塊溫度場(chǎng)計(jì)算

忽略軸向?qū)?，由圓柱坐標(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。對(duì)于式（6）的解來說，邊界條件（8）和（9）是等價(jià)的。式（6）的通解為：

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

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

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

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

1.3 芯塊溫度場(chǎng)計(jì)算修正

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

2 模型驗(yàn)證

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

2.1 模型參數(shù)

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

圖2 環(huán)形燃料芯塊COMSOL5.4計(jì)算模型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è)計(jì)和運(yùn)行中，燃料芯塊的溫度是一個(gè)非常重要的參數(shù)。對(duì)于環(huán)形燃料芯塊需要準(zhǔn)確地得到芯塊的溫度場(chǎng)，包括絕熱面位置、芯塊最高溫度等。本模型的驗(yàn)證主要針對(duì)絕熱面位置、芯塊最高溫度以及芯塊中溫度場(chǎng)分布等進(jìn)行。

2.2.1 芯塊絕熱面位置

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

2.2.2 芯塊絕熱面最高溫度

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

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

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

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

圖3 不同芯塊厚度（1 mm、2 mm、3 mm）和不同線功率密度（36 kW·m-1、54 kW·m-1、71 kW·m-1）下環(huán)形芯塊溫度場(chǎng)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)熱微分方程和傅立葉定律，針對(duì)環(huán)形燃料的芯塊建立了環(huán)形燃料芯塊一維穩(wěn)態(tài)溫度場(chǎng)計(jì)算模型，導(dǎo)出了環(huán)形燃料芯塊絕熱面位置以及燃料芯塊最高溫度的計(jì)算方法。并對(duì)該模型進(jìn)行了驗(yàn)證和分析，得出了以下結(jié)論：

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

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