基于子通道模型的超臨界二氧化碳燃料組件熱工特性分析

郭 嘉,趙富龍,*,寧可為,方華偉,易經(jīng)緯,田瑞峰,譚思超

(1.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院,黑龍江 哈爾濱 150001;2.中國核動力研究設(shè)計院,四川 成都 610213)

超臨界二氧化碳(supercritical carbon dioxide,簡稱S-CO2)反應(yīng)堆是近些年備受關(guān)注的四代堆型,在小型堆的研究及開發(fā)中被提及,其主要應(yīng)用于布雷頓循環(huán)系統(tǒng),是一種高效清潔的能源和動力轉(zhuǎn)換裝置,并且可以與多種形式的熱源結(jié)合,實現(xiàn)聯(lián)合循環(huán)發(fā)電[1]。相對于蒸汽朗肯循環(huán),S-CO2布雷頓循環(huán)系統(tǒng)具有循環(huán)效率高、回路簡單、結(jié)構(gòu)緊湊等特點,在小型堆供電及動力裝置方面具有良好的發(fā)展前景[2]。S-CO2兼具液體和氣體的特性,非常適用于作為反應(yīng)堆一回路的冷卻劑,其優(yōu)點包括熱導(dǎo)率高、黏度小、無相變、以及良好的化學(xué)穩(wěn)定性、易于獲取且成本低廉[3-4]。

目前對S-CO2燃料組件應(yīng)用子通道分析方法來進行熱工水力分析的研究十分有限,研究集中在使用S-CO2作為工質(zhì)的布雷頓循環(huán)系統(tǒng)的特性分析[5]。為了分析S-CO2工質(zhì)在反應(yīng)堆內(nèi)的熱工水力特性,獲取二維尺度上的參數(shù),本文采用適用性和可靠性高的子通道分析方法進行研究。首先選定一種燃料組件,并對其進行子通道劃分,通過連續(xù)性方程、能量守恒、動量守恒方程等構(gòu)建子通道模型,并使用Modelica語言開發(fā)燃料組件子通道分析程序[6-7]。子通道分析方法的收斂性、重復(fù)性和可擴展性要優(yōu)于CFD計算軟件,并且耗費的計算資源較少,方便后續(xù)的程序開發(fā)和程序擴展[8],為后續(xù)容積控制系統(tǒng)和布雷頓循環(huán)系統(tǒng)的程序開發(fā)和理論研究提供了技術(shù)支持[9]。子通道分析程序由于其易收斂且準確性高的特點,能夠為流動換熱特性的研究及堆芯瞬態(tài)運行特性分析提供重要的技術(shù)支持[10]。

目前關(guān)于S-CO2反應(yīng)堆的研究集中于流動換熱研究和三維CFD模擬。Sharabi等[11]針對S-CO2反應(yīng)堆模型開展流動換熱特性分析,通過三維CFD計算對傳熱惡化機理和湍流交混下冷卻劑流動進行了研究。Sun等[12]研發(fā)了子通道分析程序SUBAC,主要應(yīng)用于鈉冷快堆帶繞絲的燃料組件的分析。Ding等[13]使用Modelica語言構(gòu)造了S-CO2反應(yīng)堆瞬態(tài)模擬程序,進行了瞬態(tài)工況的研究并應(yīng)于與布雷頓循環(huán)系統(tǒng)的開發(fā)。Ming等[14]開發(fā)了S-CO2為工作介質(zhì)的布雷頓循環(huán)系統(tǒng)分析程序,研究了系統(tǒng)運行特性和設(shè)備之間的耦合特性。陳立新等[15]開發(fā)的子通道程序PRTHA應(yīng)用在西安脈沖堆上,與實驗結(jié)果比較證明了開發(fā)程序的合理性和可行性。蘇開放等[16]對氦氙子通道湍流交混的研究,通過CFD計算擬合出適用于反應(yīng)堆的交混系數(shù)關(guān)系式,量化分析了小型堆的交混效應(yīng)。吳宗蕓等[17]開發(fā)了鈉冷快堆子通道分析程序SPLICA,并通過與棒束實驗的對比,驗證了子通道分析程序的準確性。綜上所述,關(guān)于S-CO2的程序開發(fā)的研究集中于布雷頓循環(huán)系統(tǒng)的建立及參數(shù)提取,缺乏S-CO2反應(yīng)堆的精細化模型的相關(guān)研究,為彌補相關(guān)研究的不足,本文將使用子通道分析方法對S-CO2反應(yīng)堆進行程序開發(fā)。因此,本研究基于Modelica語言開發(fā)S-CO2布雷頓反應(yīng)堆系統(tǒng)子通道分析程序BRESA,對S-CO2工質(zhì)反應(yīng)堆的軸向和徑向的物理熱工參數(shù)進行計算和分析,研究各典型通道的焓升及溫升變化并驗證程序的正確性。將開發(fā)的程序進行流量變化瞬態(tài)工況研究,分析系統(tǒng)的瞬態(tài)響應(yīng)特性,得到瞬態(tài)工況下各典型通道的流量焓以及壓降變化,實現(xiàn)子通道模型的功能,旨在為S-CO2反應(yīng)堆系統(tǒng)瞬態(tài)運行特性研究及安全特性分析提供參考。

1 數(shù)值計算方法

1.1 子通道模型

子通道分析方法是分析燃料組件熱工水力特性中有效快捷的方法之一。本文通過將燃料組件在高度上分割成多段,再根據(jù)劃分的子通道將整個燃料組件中冷卻劑通道劃分為控制體,徑向通過子通道劃分,軸向上根據(jù)高度劃分成10段。燃料組件根據(jù)子通道方法控制體劃分如圖1所示。

圖1 子通道控制體劃分

構(gòu)建完燃料組件控制體的幾何模型后,需要對每個控制體建立物理熱工模型,其方法為在每個控制體上建立子通道方程,控制體對應(yīng)的燃料棒附近建立反應(yīng)堆釋熱的方程,給定初始值和邊界條件,最終得到每個控制體上的物理熱工參數(shù),對控制體建立的子通道方程包括質(zhì)量守恒方程、能量守恒方程、軸向和徑向動量守恒方程,這些方程不僅有控制體本身的方程,還涉及到不同控制體之間的方程,最后將所有控制體的方程聯(lián)立并求解,可得到整個燃料組件的物理熱工參數(shù)。

由連續(xù)性方程得到每個控制體的質(zhì)量守恒方程、能量守恒方程、軸向動量守恒方程、徑向動量守恒方程[18-20]如下:

mi,n=mi,n-1-mi,n+1-ωin,jn+ωkn,in

(1)

qin-∑?in,jn(Hin,out-Hin,in)+

∑ωin,jn(Hin,in-H*)-∑Css(Tin,out-Tin,in)

(2)

(3)

(4)

式中:mi,n為單位長度的軸向流量,kg/(m·s);ωin,jn為單位長度的橫向流量,kg/(m·s);uin,in為控制體的入口速度,m/s;Hin,in為控制體的入口的焓,kJ/kg;qin為線功率,kW/m;?in,jn為湍流交混流量,kg/(m·s);Tin,in為入口溫度,K;Cs為相鄰?fù)ǖ罒釋?dǎo)率,W/(m2·K);s為相鄰?fù)ǖ篱g距,m;pin,in為控制體的入口壓力,Pa;pin,out為控制體的出口壓力,Pa;f為摩擦阻力系數(shù);ft為湍流交混系數(shù);z為控制體高度,m;d為等效直徑,m;ρ為流體密度,kg/m3;uin,jn為控制體間橫向流速,m/s;KG為橫向形狀阻力系數(shù);s為控制體的寬度,m;下標(biāo)in,in代表控制體i,n的入口,下標(biāo)in,out代表控制體i,n的出口,角標(biāo)*代表相鄰控制體間交換的量,下標(biāo)i,n代表控制體i,n,下標(biāo)i,n-1代表控制體i,n-1,下標(biāo)in,jn代表在控制體i,n和j,n之間進行流動。

1.2 幾何模型

CO2的臨界溫度和臨界壓力分別為304.13 K和7.38 MPa,當(dāng)流體溫度和壓力分別高于臨界溫度和臨界壓力時,被視為超臨界狀態(tài),其物理性質(zhì)在此狀態(tài)下具有良好的穩(wěn)定性,本文采用超臨界態(tài)CO2作為工質(zhì),由于設(shè)計的燃料組件的溫度和壓力遠大于臨界點,所以不需考慮跨臨界問題。本文設(shè)計的燃料組件為六角形燃料組建盒,燃料棒的排列方式為三角形排列,采用37棒的燃料組件,根據(jù)燃料組件對稱性,對1/6堆芯進行作圖,并對冷卻劑通道進行逐一編號,包括10根燃料棒,并劃分成14個子通道,如圖2所示。

圖2 1/6燃料組件子通道編號圖

其中燃料組件的高度為1.0 m,燃料組件的對角線距離為164.4 mm,對邊距離為142.4 mm,燃料棒直徑為20 mm,燃料棒之間的間距為2.6 mm,燃料組件的P/D=1.13,子通道1～9的截面積為64.09 mm2,子通道11、12、13的截面積為127.68 mm2,子通道10、14的截面積為37.43 mm2。即子通道流通面積11、12、13最大,1～9次之,10、14流通面積最小。幾何參數(shù)列于表1。

表1 子通道的幾何參數(shù)

1.3 計算方法

本文基于Modelica語言開發(fā)了S-CO2布雷頓反應(yīng)堆系統(tǒng)子通道分析程序BRESA,首先根據(jù)37棒模型,取1/6堆芯,對子通道進行編號,沿軸向?qū)⒚總€子通道劃分為10個控制體,合理分配劃分控制體的數(shù)量,在節(jié)省計算資源的情況下保證計算精度,計算流程如圖3所示。

圖3 計算流程

計算步驟如下:1) 首先設(shè)定初始值,確定邊界條件;2) 需要先進行物理計算,求出燃料組件的釋熱率、冷卻劑的吸熱量等參數(shù);3) 根據(jù)冷卻劑吸熱量求解子通道控制體的能量、動量、質(zhì)量守恒方程,聯(lián)立求解;4) 對所有控制體的方程進行求解,求出第1個時間步長整個燃料組件的參數(shù);5) 每個時間步長進行1次迭代計算,將上一時間步長的值代入到下一時間步長中;6) 迭代計算直至穩(wěn)定,即燃料組件中物理計算的釋熱量和熱工計算中冷卻劑的吸熱量在空間上達到穩(wěn)定,輸出結(jié)果,實現(xiàn)子通道程序功能。

2 穩(wěn)態(tài)特性分析

2.1 子通道出口參數(shù)分析

對子通道分析程序BRESA設(shè)定入口邊界條件如下:冷卻劑入口焓為795.73 kJ/kg,燃料組件的入口流量為5.553 kg/s,入口溫度為628.15 K,入口壓力為19.9 MPa,入口處冷卻劑密度170.86 kg/m3,燃料組件功率為1 000 kW,出口邊界條件為焓976.00 kJ/kg、溫度775.15 K、壓力19.85 MPa。

穩(wěn)態(tài)計算結(jié)果的可靠性和準確性已經(jīng)在前期的論文中經(jīng)過驗證[20],受篇幅限制不做過多介紹。

在計算達到穩(wěn)定后,子通道的流量也趨于穩(wěn)定,冷卻劑各通道進出口流量如圖4所示。其中,相同截面積的子通道,例如通道1,2,3,4…,9的進口流量相同,各子通道流通面積如表1所列,子通道流通面積與本節(jié)中子通道流量呈正相關(guān),所以入口流量只與子通道的流通面積相關(guān)。冷卻劑經(jīng)過反應(yīng)堆內(nèi)加熱和交混之后,子通道出口流量相對入口流量發(fā)生了波動,子通道的出口流量相對進口流量的變化量約為0.2%～1%,證明了各子通道之間的流量交混效應(yīng)和橫向流動現(xiàn)象的存在,且會對燃料組件不同通道之間的能量交換造成一定的影響。

圖4 冷卻劑子通道入口、出口流量

冷卻劑各子通道出口焓如圖5所示,各通道入口焓均為795.42 kJ/kg。觀察可知,越靠近燃料組件中心處的子通道出口焓越高,即冷卻劑在經(jīng)過燃料組件的加熱后的焓升越大,其原因為靠近堆芯中心處的子通道相鄰的燃料棒多,因此從燃料棒吸收的熱量較多,而堆芯邊緣處子通道相鄰的燃料棒少,吸收的熱量較少,這一分布特性符合燃料組件出口處溫度和焓的分布規(guī)律。

圖5 冷卻劑子通道出口焓

冷卻劑各子通道入口、出口密度如圖6所示。各通道入口密度均為170.86 kg/m3,出口密度分布趨勢為越靠近燃料組件中心處的子通道出口密度越小,其中通道1出口密度為124.99 kg/m3,通道14的出口密度為142.82 kg/m3。由此可知,在經(jīng)過燃料組件的加熱后,隨著冷卻劑焓的升高,冷卻劑密度逐漸降低,焓高的子通道密度較小,焓低的子通道密度較大。

圖6 冷卻劑子通道入口、出口密度

2.2 典型通道沿軸向參數(shù)分析

子通道分析程序BRESA的典型通道沿軸向焓升變化如圖7所示。由于相臨近的子通道焓近似,為便于觀察,取典型通道1、5、11、14進行繪圖。其中子通道1的焓升最大,為233.93 kJ/kg,通道5焓升為209.44 kJ/kg,子通道11焓升為161.70 kJ/kg,通道14焓升最小,為121.60 kJ/kg。觀察可知,冷卻劑在軸向上焓變化呈現(xiàn)在進出口焓升較小、軸向中間處焓升較大的特點,即燃料組件的熱量集中在軸向的中間處。各通道的吸熱量也存在差異,通道1吸熱最多,軸向上溫度變化劇烈,靠近堆芯邊緣的通道11和14吸熱量小,軸向上焓變化平緩。

圖7 典型通道沿軸向焓升變化

典型通道沿軸向溫升變化如圖8所示。其中子通道1溫升最大,為189.69 K,子通道14溫升最小為99.06 K。觀察可知,冷卻劑在軸向上的的溫度分布和焓分布近似,靠近中心的冷卻劑通道溫升大,靠近邊緣的冷卻劑通道溫升少。且在軸向中間處溫升變化較大。

圖8 典型通道沿軸向溫升變化

對比圖9中Zhang等的S-CO2子通道研究中各典型通道的軸向溫度分布結(jié)果[21],由于S-CO2工質(zhì)在該設(shè)計工況下,溫度和焓的變化趨勢基本一致,則可用Zhang等的溫度分布結(jié)果來驗證本文中溫度和焓分布是否合理。本文中典型通道的冷卻劑焓和溫度在軸向上的變化趨勢與Zhang等的結(jié)果一致,均為斜率先增加后減少的上升曲線,并且程序得到的結(jié)果符合理論分析的結(jié)果,滿足中間通道溫升和焓升較高,邊緣通道溫升和焓升較低的分布特點,證明了子通道分析程序結(jié)果的合理性,為后文中瞬態(tài)計算提供了理論基礎(chǔ)和程序支持。

圖9 冷卻劑軸向溫度分布[21]

3 瞬態(tài)特性分析

3.1 瞬態(tài)工況下流量變化分析

本文針對子通道燃料組件模型,確定模型的邊界條件,引入流量變化的瞬態(tài)工況,觀察系統(tǒng)參數(shù)隨時間變化情況,分析系統(tǒng)的瞬態(tài)響應(yīng)特性。

對子通道分析程序BRESA設(shè)置流量變化如下:設(shè)定額定流量為5.553 kg/s;從1 000 s到1 500 s時間段流量線性增加,從5.553 kg/s增加到11.106 kg/s;從1 500 s到2 000 s時間段流量線性減少,從11.106 kg/s減少到8.330 kg/s;之后在極短的時間內(nèi)急劇減少到額定流量5.553 kg/s。共計3個瞬態(tài)過程:一次流量線性增加,一次流量線性降低,一次流量突降。觀察子通道燃料組件模型參數(shù)隨時間的變化。

圖10為在設(shè)置的瞬態(tài)工況下,典型通道的流量變化。由于截面積相同的子通道流量近似,為便于觀察,取典型通道1、5、11、14進行繪圖。其中,子通道11流量最大,子通道1和通道5流量近似,子通道14流量最少。觀察可知,每個子通道的流量變化趨勢基本與燃料組件的流量變化趨勢保持一致,在一定范圍內(nèi)的總流量改變,對冷卻劑通道的流量分配基本不會產(chǎn)生的影響,且各冷卻劑通道的流量也呈比例變化。

圖10 典型通道流量變化

3.2 瞬態(tài)工況下焓變化分析

在設(shè)置的流量變化瞬態(tài)工況下,反應(yīng)堆會進行快速響應(yīng),功率會進行快速的變化,從而使冷卻劑焓發(fā)生變化,由于燃料和冷卻劑的溫度負反饋系數(shù)的影響以及采用的PID算法調(diào)節(jié)的功率控制系統(tǒng),反應(yīng)堆的功率和冷卻劑焓不是隨流量線性變化的。圖11為在流量瞬態(tài)變化下的功率變化曲線,式(5)為PID控制方程[14]:

圖11 燃料組件功率及流量變化曲線

(5)

式中:u(t)為反饋信號;e(t)為相對誤差;Kp為比例增益;Tt為積分時間常數(shù);Td為微分時間常數(shù)。

在設(shè)置瞬態(tài)工況下,冷卻劑流量線性增加過程中,冷卻劑溫度和焓會下降,通過溫度負反饋的影響,反應(yīng)堆會通過功率控制系統(tǒng)引入正反應(yīng)性,從而使反應(yīng)堆的功率上升,冷卻劑溫度和焓又隨之上升,溫度負反饋又會導(dǎo)致功率有小幅下降,由于冷卻劑流量還處于增加的過程,冷卻劑溫度和焓上升一段時間后又會繼續(xù)降低,因此,功率呈現(xiàn)波動上升的趨勢,冷卻劑焓呈現(xiàn)波動下降的趨勢。冷卻劑流量線性減少過程同理,功率波動下降。冷卻劑流量突降時,反應(yīng)堆功率也出現(xiàn)突降,然后產(chǎn)生大幅波動,最終穩(wěn)定到一個新的水平。

本文的控制目標(biāo)為在流量變化的瞬態(tài)工況下,保證反應(yīng)堆功率和冷卻劑焓能夠隨流量變化,避免發(fā)生反應(yīng)堆超功率或冷卻劑溫度過高等工況,確保反應(yīng)堆的安全運行。

圖12為子通道分析程序BRESA的典型通道的焓變化。由于相臨近的子通道焓近似,為便于觀察,對典型通道1、5、11、14進行繪圖。觀察可知,在流量變化的瞬態(tài)工況下,冷卻劑通道的焓會發(fā)生一定量的波動。流量線性增加的過程中,冷卻劑通道的焓呈現(xiàn)波動式下降趨勢,這是因為發(fā)生流量降低的瞬態(tài)過程中,冷卻劑流量降低,為了防止冷卻劑吸熱量降低導(dǎo)致堆芯溫度上升,反應(yīng)堆功率會緩慢下降,同時冷卻劑通道吸熱量降低,冷卻劑的焓也會緩慢下降。流量線性減少的過程中,同理,冷卻劑通道的焓會波動式上升。

圖12 典型通道焓變化

但當(dāng)冷卻劑流量急劇下降時,冷卻劑焓會急劇升高,在流量下降瞬間的典型通道焓升達到90～120 kJ/kg,焓波動幅最大達到12%,其原因為冷卻劑流量的急劇減少,會使冷卻劑的吸熱量急劇減少,造成傳熱惡化以及堆芯溫度升高的現(xiàn)象,這不利于反應(yīng)堆安全。由于反應(yīng)堆和冷卻劑溫度升高,反應(yīng)堆會進行快速響應(yīng),引入負反應(yīng)性來降低反應(yīng)堆功率,在這個過程中冷卻劑焓和反應(yīng)堆功率會回落并發(fā)生大幅波動,直到恢復(fù)下一個穩(wěn)定狀態(tài)。

焓的變化曲線證明了流量波動在一定范圍內(nèi),焓波動量較小,能夠保證反應(yīng)堆的安全,但流量短時間快速變化,會造成焓突變和劇烈波動,對反應(yīng)堆的安全運行帶來挑戰(zhàn),不利于安全運行,因此在實際工程中應(yīng)避免冷卻劑流量急劇變化。

3.3 瞬態(tài)工況下壓降變化分析

圖13為在設(shè)置的瞬態(tài)工況下,子通道分析程序BRESA的壓力變化。觀察可知,冷卻劑流過燃料組件后的壓降變化與流量變化趨勢相似,流量增加,燃料組件出口壓力減少,即冷卻劑在燃料組件中壓降增大,當(dāng)流量從5.553 kg/s增加到11.106 kg/s,壓降從11 970 Pa變化為41 315 Pa。在流量突然降低的瞬態(tài)過程中,燃料組件壓降存在波動,其原因與瞬態(tài)工況下焓變化原因近似,冷卻劑吸熱量的大幅變化,使冷卻劑密度也會發(fā)生變化,從而導(dǎo)致燃料組件壓降也會出現(xiàn)波動。分析可知,冷卻劑流量變化對壓降影響較大,這會導(dǎo)致燃料組件承受的應(yīng)力變化較大,因此在流量變化的瞬態(tài)工況下對燃料組件的材料要求較高。

圖13 燃料組件壓力變化

流量變化瞬態(tài)計算結(jié)果表明,應(yīng)用子通道方法的S-CO2燃料組件中,冷卻劑流量的變化會對壓降和焓產(chǎn)生一定的影響,當(dāng)流量隨時間線性緩慢變化時,焓和壓降的變化也相對平緩,但當(dāng)流量在短時間內(nèi)快速變化時,焓和壓降都發(fā)生劇烈變化,這對S-CO2反應(yīng)堆的安全運行帶來挑戰(zhàn),發(fā)生包殼融毀,燃料組件應(yīng)力腐蝕等事故的概率增大。為了完善反應(yīng)堆的安全設(shè)計,需要設(shè)計更科學(xué)高效的容積控制系統(tǒng),保證冷卻劑流量穩(wěn)定、準確、快速、平穩(wěn)的進行調(diào)節(jié)。后續(xù)研究將對S-CO2反應(yīng)堆容積控制系統(tǒng)進行設(shè)計和分析,本文的研究內(nèi)容為S-CO2反應(yīng)堆系統(tǒng)瞬態(tài)運行特性分析提供參考。

4 結(jié)論

本文基于Modelica語言建立了S-CO2反應(yīng)堆子通道分析程序BRESA,對燃料組件內(nèi)各子通道的流量、密度、焓,以及子通道在軸向上的焓升和溫升進行了參數(shù)提取和數(shù)據(jù)分析,得到了穩(wěn)態(tài)運行工況下燃料組件的關(guān)鍵參數(shù),證明了開發(fā)的子通道分析程序的合理性,并進行了瞬態(tài)工況分析,主要結(jié)論如下。

1) 在燃料組件中,采用子通道分析方法,由于子通道模型的橫向流動的影響,冷卻劑流量會有0.2%～1%的變化量,燃料組件熱工參數(shù)呈現(xiàn)中間通道焓升較高、邊緣通道焓升較低的分布特點,驗證了程序的正確性和精確性。

2) 在流量變化的瞬態(tài)工況下,當(dāng)流量線性緩慢降低和升高時,反應(yīng)堆功率與流量變化趨勢相同,但呈現(xiàn)波動變化;焓和出口壓力與流量變化呈現(xiàn)相反的趨勢,并波動變化。這證明了反應(yīng)堆在冷卻劑流量小幅變化的工況,燃料組件參數(shù)在控制系統(tǒng)和流量變化的雙重作用下會發(fā)生小幅波動變化。該工況下反應(yīng)堆運行參數(shù)的波動在一定范圍內(nèi),反應(yīng)堆處于安全狀態(tài)。

3) 在流量變化的瞬態(tài)工況下,流量突然變化會導(dǎo)致功率、焓和壓降產(chǎn)生劇烈波動,功率發(fā)生劇烈波動并快速下降,焓先急劇升高后大幅波動下降,瞬時焓升達到90～120 kJ/kg,流量突變的瞬態(tài)工況對S-CO2反應(yīng)堆的安全運行帶來挑戰(zhàn),不利于反應(yīng)堆的安全運行。