燃料組件偏置及彎曲對沸騰臨界影響研究

龔佳奇,叢騰龍,肖 瑤,顧漢洋,丁 銘

(1.哈爾濱工程大學 核科學與技術學院,黑龍江 哈爾濱 150001;2.上海交通大學 機械與動力工程學院,上海 200030)

臨界熱流密度限制了反應堆堆芯的最大功率密度,與反應堆的安全性、經濟性密切相關[1-3],是熱工設計中必須要確定的限值[4]。燃料組件在堆內由于制造的偏差、外力的作用,或在臨界熱流密度熱工實驗因為安裝的偏差,導致燃料棒整體偏置或軸向彎曲,引起冷卻劑通道流通截面的幾何形狀和大小偏離設計值,可能對沸騰臨界造成影響。Nagino等[5]對4×4帶定位格架棒束組件進行了臨界實驗,其中軸向和徑向功率都不均勻,沿軸向位置布置2個攪混格架及1個位于中間的支撐格架,兩攪混格架中間的1根高功率棒發生彎曲,在兩攪混格架中點與附近的高功率棒及導向管發生接觸(即彎曲程度達到最大)。結果表明棒彎曲在高壓下對DNBR的影響很大、低壓下沒有影響。Tamai等[6]對37稠密柵帶定位格架棒束組件進行了臨界實驗,測試段包含4個定位格架,在靠近出口兩格架間發生燃料棒單棒彎曲,彎曲形狀呈余弦形,即彎曲段中間達到彎曲最大值,并與相鄰棒發生接觸。發現彎棒使臨界數值比未彎曲棒束降低10%,臨界位置發生在棒接觸點下游。Markowski等[7]分析了5×5棒束組件彎棒對于CHF的影響,結果表明,當彎曲使相鄰加熱棒間隙減至50%以上,才觀察到CHF降低,未加熱的彎曲棒(冷壁)甚至與加熱的桿接觸也沒有不利影響。Ren等[8]對5×5帶定位格架的彎棒進行了數值模擬,彎曲位置發生在兩格架間。結果表明,當彎曲程度增加到燃料棒間隙變化的50%以上時會減小CHF,在0～50%時沒有影響,預測的DNB位置發生在彎曲段。以上實驗或數值模擬采用的燃料棒彎曲位置均為內部棒且彎曲僅發生在兩格架位置間,彎曲并未對邊角通道產生直接影響,并未關注棒束全高度整體彎曲、偏移對臨界的影響。在早期的OECD開展的實驗及近期上海交通大學開展的CHF實驗中發現,燃料棒可能發生大跨度的彎曲或由于加工工藝出現小幅的偏移。

為研究棒彎曲和偏移對CHF的影響,本文基于歐拉兩流體模型和壁面沸騰模型[9-12],研究發生彎曲或偏置的5×5燃料組件的CHF,并比較3種變形對于CHF數值和位置的影響,為燃料組件的安全分析提供依據。

1 數值模型

過冷沸騰和沸騰臨界的液相與氣相溫度并不相同,即主流過冷,近壁面可能過熱[13]。另外,由于密度差和相間作用,兩相之間出現滑移速度,因此應考慮相間的不平衡[14]。兩流體模型可考慮相界面質量、動量和能量的傳輸,基于歐拉兩流體模型分別建立氣、液兩相守恒方程[15-17]。同時,引入壁面熱流密度分配模型描述換熱過程壁面向氣、液兩相的熱量傳遞過程。此外,為使方程封閉,引入輔助模型,如氣泡脫離直徑、脫離頻率等[18]。

1.1 歐拉兩流體模型

質量守恒方程:

(1)

動量守恒方程:

(2)

能量守恒方程:

(3)

1.2 壁面沸騰模型

CFD方法可基于兩種方法模擬兩相流,即界面捕獲法和相位平均法。第1種的計算網格尺寸需遠小于氣泡尺寸,同時為了保持數值穩定性的時間步長,該方法需要巨大的計算資源。第2種基于統計界面面積濃度和一系列輔助模型來模擬兩相流[19]。壁面沸騰模型是受熱壁面沸騰應用最廣泛的經驗模型,為了模擬臨界,本研究采用拓展后的壁面沸騰模型,即CHF模型。該模型將壁面熱流密度分為4部分:

qwall=(qc+qq+qe)f(αl)+(1-f(αl))qg

(4)

式中:qwall為壁面熱流密度;qc、qq、qe、qg分別為液體單相對流傳熱、淬滅對流傳熱、蒸發熱流傳熱及氣體單相對流傳熱部分的熱流密度,表達式[19]如下:

qc=hl(Tw-Tl)(1-Ab)

(5)

(6)

qe=VdρgNwhfvf

(7)

qg=hg(Tw-Tg)

(8)

式中:hl為單相換熱系數;Tw為壁面溫度;T為液相溫度;1-Ab為加熱壁面液體所覆蓋的面積;Vd為脫離壁面的氣泡體積;Nw為成核密度;ρg為氣相密度;hfv為蒸發潛熱;f為氣泡脫離頻率;kl為液相導熱系數;cpl為比定壓熱容。

f(αl)為加熱壁面液相所占的份額,表達式為Tenter[20]關系式:

f(αv)=1-f(αl)=

(9)

其中,αv,1和αv,2為過渡點,分別為0.9和0.95。

1.3 界面傳輸模型

歐拉兩流體模型守恒方程中,源項考慮了界面傳輸項,因此需使用界面傳輸模型來進行求解。對于相間傳質,考慮了兩種機制。第1個是加熱表面的液體蒸發及氣泡進入過冷主流的冷凝或主流中過熱液的氣化。對于相間動量傳遞,考慮了曳力和非曳力。曳力由流體的分子黏度和兩相間的滑移速度所引起,作用于滑移速度的負方向。非曳力作用在與流動方向正交的方向上,包括速度梯度引起的升力、將氣泡推向主流的壁面潤滑力及展平徑向空泡份額的湍流耗散力。考慮從界面到液相和從界面到氣相的熱阻,通過雙熱阻模型模擬相間傳熱。氣泡對湍流強度的影響通過氣泡尾流引起的附加湍流源項來模擬,相關模型列于表1,對于相間動量模型,Ishii模型適用于沸騰流動不同尺寸、形狀的氣泡的曳力計算。Moraga升力系數模型對液滴和氣泡所受升力適用性良好。Antal對于近壁面氣相速度分布差異大的區域能夠適用。Burns基于相間曳力的計算,與實驗值的依賴性相對Bertodano弱。針對相間傳質,文中計算均為DNB型沸騰臨界,主要的物理現象僅存在受熱壁面蒸發及進入主流冷凝,根據壁面沸騰模型中的蒸發熱流來推導相間傳質。對于相間能量傳遞,通常假定兩相界面溫度為相應壓力下飽和溫度,分別計算界面向液相間的傳熱,和界面向氣相間的傳熱。表1中:hfv為汽化潛熱;Tsat為飽和溫度;Tl為液相溫度;Tv為氣相溫度;h為比焓;Ai為單位體積內氣液兩相間相界面面積;db為氣泡直徑;σ為表面張力;g為重力;CD為曳力系數;CL為升力系數;CWL為壁面潤滑力系數;Ftd為湍流耗散力系數;yw為網格中心到最近壁面的距離;σvl為常數0.9;CTD為常數1;μl為液相動力黏性系數;αv為氣相份額;αl為液相份額;Nub為努塞爾數;Pr為普朗特數;qvt為單位體積內氣相和液相相界面間的換熱量;δt為時間尺度默認值為0.05(FLUENT默認);νv為運動湍流黏度;k為湍動能;ε為湍流動能耗散率;Uv為氣相平均速度;Ul為液相平均速度。

表1 界面傳輸相關模型Table 1 Correlation model for interfacial transfer

表2 各工況變形類型Table 2 Deformation types under various conditions

1.4 數值方法

通過將非線性偏微分方程離散化為線性代數方程來求解。本文求解穩態控制方程,因此不需要時間離散化,通過偽瞬態方法求解以提高收斂性。動量方程中的壓力項采用二階方法離散化,梯度項通過基于網格的最小二乘法得到。離散得到的線性代數方程用Guass-Seidel方法求解,質量和動量方程同時用COUPLED算法求解。

2 幾何、網格劃分及邊界條件

考慮到本文研究重點為偏置和彎曲對CHF的影響,為排除定位格架影響,本文研究對象為5×5光棒的燃料組件,包殼外徑為9.5 mm,棒束柵距為12.6 mm,橫截面結構如圖1所示,同時給出子通道劃分以及燃料棒編號。組件加熱高度設置為2 m,不加熱的出口段長約130 mm,燃料組件三維幾何結構如圖2所示。

圖1 幾何截面及子通道劃分示意圖Fig.1 Schematic diagram of geometric section and sub-channel division

圖2 三維幾何結構示意圖Fig.2 Schematic diagram of three-dimensional geometric structure

在網格劃分方面,首先利用ICEM對無偏置和彎曲的基準算例進行六面體網格劃分,經網格獨立性分析,最終網格總數409萬。在基準算例的基礎上通過動網格技術,使網格發生偏置與彎曲,共設計3種類型,即X型(橫向偏置)、XY型(對角線偏置)、C型(軸向彎曲)。3種類型網格均為整體移動(即所有包殼都發生移動),以下在3種類型中各挑選1個算例網格進行展示,如圖3所示,其中X、XY偏置示意圖選擇右上角4個子通道進行放大,圖中展示結果均為變形1 mm工況的網格示意圖。圖4示出X、XY、C型變形具體示意圖,3種變形均為1 mm,X、XY、C型變形方向為x軸右側。圖中對于中心坐標軸,X型偏置所有燃料棒向右側位移,XY型偏置向對角位移。圖中C型彎曲,幾何進口位置對應弧度0,加熱段末端對應弧度3.14,而彎曲最大值(1 mm)在加熱高度中心位置(1 m),整體彎曲效果示于圖3。另外動網格移動后的網格正交質量均在0.55以上。

圖3 基準算例、X型、XY型及C型網格示意圖Fig.3 Schematic diagrams of base case, X-type, XY-type and C-type meshes

圖4 3種變形示意圖Fig.4 Schematic diagrams of three kinds of deformation

本文在基準工況的基礎上進行棒偏置和彎曲影響規律研究,基準工況中設計壓力15.5 MPa、質量流量3 000 kg/(m2·s)、出口含氣率0、進口溫度534.66 K及CHF參考值2 424 kW/m2。根據壓力、含氣率及質量流密度查詢CHF查詢表(look-up table)[28],再對查詢值通過水力直徑和棒束修正因子進行修正,從而得到最終的CHF查詢表數值,即參考值。在數值模擬過程中,從0.5倍參考值逐步提升功率,功率在50%～70%區間以10%提升功率,超過70%時以1%提升功率直至臨界。

在基準工況基礎上,通過動網格技術實現燃料棒的整體偏置或彎曲,獲得新結構下的計算網格。設計以下工況,各工況類型列于表3。將需要計算的工況分兩類,A類9個工況軸向功率均勻分布,又分3組,即X型、XY型、C型,每組3個變形量分別為0.5、0.8、1.0 mm;B類2個工況軸向功率呈余弦分布[29],均為C型。其中,1.0 mm的變形屬于較大的變形,0.5 mm和0.8 mm分別為小變形和中等變形工況。目前公開的文獻中,并沒有對變形量取值的直接依據,實驗中通常以兩燃料棒彎曲至接觸的局部現象出現,為了簡化問題,文中采用沿1個方向變形的方式來處理。CHF數值對于邊角通道間隙寬度很敏感,這樣的結果使得變形1 mm對于CHF的數值及位置影響非常明顯,同時0.5 mm的變形臨界位置依然在熱通道,綜上選取變形量0.5、0.8、1 mm。

表3 邊界條件Table 3 Boundary condition

由于燃料包殼外表面周向熱流密度分布較為均勻,在建立燃料組件模型時可不考慮燃料棒內的導熱過程,直接在包殼外表面上給定熱流密度。軸向加熱功率包括均勻和余弦分布兩類,詳細的邊界條件列于表2、3。

3 模型驗證

將上述工況預測的CHF與查詢表數值進行對比,結果列于表4。查表法是將FLUENT計算臨界高度位置發生臨界的局部子通道水溫、空泡份額、流量及壓力取出,用水溫及空泡份額轉換為含氣率,根據局部子通道壓力、流量、含氣率及當量直徑,獲得最終的CHF查詢表數值。表中列出5個工況的CFD預測值、查詢表數值及誤差。其余工況,臨界位置均由中心遷移到邊緣,而在邊、角通道位置處CHF查詢表并不適用。這是因為CHF查詢表基于8 mm的圓管制成,對于非8 mm的圓管或者棒束結構(如內部子通道)可通過水力直徑修正因子進行修正,往往預測精度高。然而在棒束組件的邊角通道位置,受熱邊界及濕周并不均勻,加上燃料棒偏置或彎曲對邊角通道進行擠壓,導致該位置受熱邊界及濕周分布更不勻衡,因此CHF查詢表不再適用[30]。而對于臨界位置發生在熱通道的工況,CFD預測值與CHF查詢表數值誤差均在5%以內,預測精度良好。

表4 查詢表數值與預測值對比Table 4 Comparison of sub-channel look-up table value and prediction value

4 結果分析

臨界熱流密度主要受流量、壓力、含氣率參數影響。圖5示出無變形的基準工況臨界高度位置液相溫度和空泡份額云圖分布,圖中可發現中心子通道溫度最高,臨界位置也發生在該處(黑色虛線框住包殼位置),邊角通道溫度次之。空泡份額峰值位于熱棒最小間隙處,最大值約為0.94,觸發沸騰臨界。圖6示出基準算例空泡份額及溫度隨流向分布的趨勢,圖中截面最大空泡份額(即包殼)及截面平均空泡份額均隨流向而增加,在出口達到最大,趨勢與主流溫度曲線基本一致;主流溫度出口位置仍低于飽和溫度,為DNB型沸騰臨界。包殼溫度先升高,中間段達到624 K左右保持不變,在靠近出口位置發生飛升。

圖5 基準算例臨界高度處的液相溫度及空泡份額云圖分布Fig.5 Cloud maps of water temperature at critical height and vapor fractions distribution of benchmark

燃料棒的整體偏置或彎曲使不同子通道冷卻劑流量重新分配,從而影響傳熱和沸騰臨界。圖7～10及表5比較了不同類型工況邊、角及中心子通道總質量流量分布、沸騰臨界位置和臨界高度處溫度云圖。圖中不同顏色代表不同變形程度,對應顏色標記了沸騰臨界發生位置。由表5可知,對于同種工況,15號子通道質量流量最大,邊通道(12號)次之,角通道(6號)最小。對于同類變形的不同工況,隨變形程度的增大,中心子通道(15號)質量流量通道增加,邊角子通道(12、6號)質量流量降低。這是由于燃料棒整體偏置或彎曲,對右側邊角通道造成擠壓使得通道流通面積減小,冷卻劑流量降低。對于X型和XY型相同偏置距離的工況,XY型12號邊通道質量流量降幅更小,15號子通道上升幅度更小(除位移0.5 mm非常接近),X型偏置工況相比XY型有更多邊通道的流量向中心通道流動;而在6號角通道質量流量降幅更大,在位移達到1 mm時,XY型6號角通道質量流量降幅比X型高6.5%,因為燃料棒沿對角線位移,若將其正交分解,燃料棒沿x和y方向同時擠壓6號角通道,導致該處流通面積極窄,冷卻劑流量降幅更大。即X型對于邊通道的影響比XY型大,在角通道的影響弱于XY型偏置,對于整體流量的影響X型更強。對于C型彎曲工況,15號子通道除軸向功率均勻彎曲1 mm的工況比X型同變形程度稍小,其余C型工況中心子通道上升幅度比X和XY型都更大,而功率曲線為余弦分布時中心通道上升幅度最大。整體而言,結合邊棒臨界情況C型彎曲對于燃料組件子通道流場的影響最大。

圖7 X型臨界高度處徑向臨界位置及液相溫度分布Fig.7 Radial critical position and water temperature distribution at X-type critical height

表5 不同工況子通道總質量流量分布Table 5 Total mass flux distribution of sub-channels under different conditions

這種流量分布的結果,導致溫度云圖的分布趨勢發生改變,高溫區集中分布在被擠壓的邊角通道,低溫區集中遠離燃料棒位移或彎曲方向的子通道。圖7中可發現X偏置0.8、1 mm時邊棒臨界,偏置0.5 mm和基準工況(黑色)臨界發生在熱通道,溫度云圖高溫區向位移方向移動,集中在右側邊棒位置。圖8中最高溫度區域隨位移發生移動,集中在6號角通道,31號角通道冷卻劑溫度最低。XY偏置0.5、0.8 mm的臨界位置發生在熱通道,偏置1 mm邊棒臨界,僅在6號角通道位置觸發臨界,流量分析分布已解釋該現象。

圖8 XY型臨界高度處徑向臨界位置及液相溫度分布Fig.8 Radial critical position and water temperature distribution at XY-type critical height

對于軸向均勻功率及非均勻功率的C型彎曲工況結果如圖9、10所示,圖中可發現溫度云圖分布趨勢與X型工況類似。對于軸向功率均勻的彎曲工況,0.8、1 mm時邊棒臨界,彎曲0.5 mm時臨界位置發生在熱通道;而軸向功率呈余弦形,彎曲0.5、0.8 mm時邊棒臨界。即彎曲0.5 mm不同功率曲線,徑向臨界位置發生變化。

圖9 C型臨界高度處徑向臨界位置及液相溫度分布Fig.9 Radial critical position and water temperature distribution at C-type critical height

圖10 C型(cosine)臨界高度處徑向臨界位置及液相溫度分布Fig.10 Radial critical position and water temperature distribution at C-type (cosine) critical height

CHF限制反應堆功率輸出水平,它關系到反應堆堆芯的安全性和經濟性。將上述12個工況預測的CHF和臨界高度匯總成數據圖,如圖11所示。從圖11可看出,無論是均勻加熱或非均勻加熱功率,對于整體趨勢,X、XY、C型變形程度越大,CHF越小。軸向非均勻功率與彎曲共同作用下的工況CHF最低。對于均勻加熱功率,相同變形程度,XY型對CHF影響最小;以基準工況為參考,位移1 mm 的X型工況預測的CHF比彎曲1 mm偏差更大;對于0.5、0.8 mm變形,X型預測的CHF比C型更接近基準工況CHF。綜上,在某個變形量閾值前,C型對CHF影響更大,在閾值后X型對CHF影響更大。

圖11 不同工況CFD預測臨界值對比Fig.11 Comparison of CFD predicted critical values under different conditions

圖12示出不同工況臨界高度對比。對于軸向均勻功率,X、XY型和彎曲0.5 mm工況不影響臨界高度,僅0.8、1 mm的C型彎曲改變臨界高度,隨彎曲程度的增大,臨界位置向上游小幅遷移。而軸向非均勻功率的兩個C型工況中,臨界高度大幅前移,位于中部偏下游位置。

圖12 不同工況軸向臨界位置對比Fig.12 Comparison of axial critical positions under different conditions

5 結論

本文基于歐拉兩流體模型對發生偏置或彎曲的5×5棒束燃料組件CHF進行了數值模擬研究,獲得了橫向偏置(X型)、對角偏置(XY型)和彎曲(C型)對CHF的影響,通過分析可得到如下結論。

1) 對于軸向均勻功率的工況,X型偏置及C型變形0.8、1 mm時,出現邊棒臨界;XY型偏置1 mm出現邊棒臨界。對軸向余弦功率的工況,彎曲0.5、0.8 mm均發生邊棒臨界;彎曲0.5 mm余弦與均勻功率的工況對比,帶余弦功率工況觸發邊棒臨界,即余弦功率影響邊棒臨界。當變形程度增加時,邊棒臨界的可能性增加。

2) 對軸向功率均勻的工況,X、XY型和0.5 mm的C型彎曲對于CHF位置沒有影響。0.8、1 mm的C型彎曲工況,隨彎曲程度的增大,臨界位置向上游移動;而對于軸向余弦功率的C型彎曲,趨勢也一致。

3) 3種變形程度越大,預測的CHF越小;XY型變形對CHF敏感性最弱,軸向功率余弦的C型工況對CHF影響最大,相比同等彎曲0.5 mm均勻功率工況,臨界數值降低9%。對于均勻功率的X、C型變形工況,0.5、0.8 mm C型彎曲的臨界數值偏差更大(與基準工況CHF對比),而變形1 mm的X型臨界數值偏差更大,即兩類變形對臨界數值的影響在0.8～1 mm區間存在交點,在交點前C型對臨界數值影響更大,交點后X型對臨界數值影響更大。

4) 燃料棒的整體偏置或彎曲使子通道流量重新分配,在邊角通道受燃料棒偏置彎曲影響大,XY型對6號角通道影響最大,而對燃料組件整體流量分配影響最小;C型彎曲對燃料組件整體流量分配影響最強。當邊角通道縫隙小到冷卻劑流量不能有效帶走燃料棒表面氣泡,造成氣泡局部擁塞,將觸發邊棒臨界。