曹今毅, 孫建闖, 孟祥飛, 張文超, 杜利鵬, 蔡偉華

(1.東北電力大學(xué) 能源與動(dòng)力工程學(xué)院, 吉林 吉林 132012; 2.東北電力大學(xué) 熱流科學(xué)與核工程實(shí)驗(yàn)室, 吉林 吉林 132012; 3.中廣核研究院有限公司, 廣州 深圳 518031)

相較于傳統(tǒng)圓柱型燃料棒,螺旋花瓣型燃料棒的表面積與體積比較大,有效增大了換熱面積;螺旋結(jié)構(gòu)強(qiáng)化冷卻劑橫向交混,顯著增強(qiáng)對(duì)流換熱;燃料組件采用自支撐布置,進(jìn)一步簡(jiǎn)化堆芯結(jié)構(gòu)[1]。基于上述優(yōu)勢(shì),花瓣形燃料棒束可采用緊密排列布置技術(shù),在小型反應(yīng)堆中應(yīng)用前景廣闊。俄羅斯最早開(kāi)展螺旋花瓣型燃料棒相關(guān)研究,并應(yīng)用在核動(dòng)力破冰船中[2-3]。Nikolai等[4]在氦冷快堆堆芯中采用螺旋花瓣型燃料組件,以及二氧化钚和貧化鈾混合物顆粒的核燃料。Conboy等[5]研究發(fā)現(xiàn),在相同工況下,采用花瓣形燃料的壓水堆堆芯,功率比傳統(tǒng)圓柱形燃料棒增加47%。Garusov等[6]在額定功率下,考慮流通截面和表面粗糙度影響,得到了花瓣形棒束通道阻力系數(shù)和努塞爾數(shù)關(guān)聯(lián)式。Zhang等[7-8]研究了單相冷態(tài)條件下5×5花瓣形棒束通道橫向交混特性,得到了棒束通道內(nèi)阻力系數(shù)關(guān)系式。Xiao等[9]采用子通道方法,分析了花瓣形燃料棒束通道內(nèi)的流場(chǎng)分布,發(fā)現(xiàn)中心子通道和邊子通道的橫向流動(dòng)方向不同。相關(guān)學(xué)者也采用計(jì)算流體力學(xué)(computational fluid dynamics,CFD)方法研究花瓣形燃料棒束內(nèi)流動(dòng)換熱機(jī)理。蔡偉華等[10]基于CFD方法研究了花瓣形燃料棒流動(dòng)換熱特性,結(jié)果表明,螺旋節(jié)距小于750 mm的棒型,螺旋結(jié)構(gòu)能有效增強(qiáng)燃料棒換熱能力。杜利鵬等[11]采用數(shù)值研究發(fā)現(xiàn),花瓣形結(jié)構(gòu)會(huì)使流速不均,內(nèi)凹弧產(chǎn)生明顯過(guò)冷沸騰。Cong等[12]采用數(shù)值方法研究了螺旋花瓣型燃料棒束通道內(nèi)流動(dòng)沸騰和臨界熱通量,發(fā)現(xiàn)花瓣形燃料臨界熱通量顯著高于圓柱形燃料。

將花瓣形燃料棒應(yīng)用于海基核動(dòng)力裝置中,需考慮搖擺運(yùn)行狀態(tài)對(duì)自然循環(huán)的影響。譚思超等[13]通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)燃料棒蓄熱對(duì)壁溫波動(dòng)幅值的影響較大。此外,搖擺附加慣性力會(huì)產(chǎn)生附加壓降,進(jìn)而影響系統(tǒng)的自然循環(huán)能力[14]。Lai等[15]采用六自由度搖擺臺(tái)實(shí)驗(yàn)研究了不同運(yùn)動(dòng)條件下自然循環(huán)流動(dòng)換熱特性,提出了適用于湍流區(qū)的周期平均摩擦阻力系數(shù)的預(yù)測(cè)關(guān)聯(lián)式。近年來(lái),相關(guān)學(xué)者提出了一維與三維多尺度耦合模擬方法,不但能獲得復(fù)雜堆芯三維流場(chǎng)信息,而且系統(tǒng)瞬態(tài)響應(yīng)特性與實(shí)驗(yàn)吻合更好。Zhang等[16]采用固定迭代步法的數(shù)據(jù)交互策略,通過(guò)預(yù)設(shè)固定迭代步數(shù)作為數(shù)據(jù)交互判據(jù),研究了3×3圓柱形棒束通道流動(dòng)特性。Cai等[17]針對(duì)現(xiàn)有固定迭代步法的數(shù)據(jù)交互策略的不足,提出了基于自適應(yīng)數(shù)據(jù)交互的多尺度耦合數(shù)值模擬方法,將其應(yīng)用于自然循環(huán)系統(tǒng),大大提高了計(jì)算效率和計(jì)算結(jié)果的準(zhǔn)確性。

含螺旋花瓣型燃料棒束通道自然循環(huán)系統(tǒng)流動(dòng)受搖擺影響特性尚不清楚,極大地限制了花瓣形燃料棒的工程應(yīng)用。本文采用區(qū)域分解方法,建立了多尺度耦合數(shù)值模型,分析搖擺運(yùn)動(dòng)下含花瓣形棒束通道自然循環(huán)流動(dòng)特性,并修正瞬時(shí)阻力系數(shù)關(guān)系式,為該類型燃料棒的工程應(yīng)用提供參考。

1 搖擺條件下多尺度耦合模擬方法

自然循環(huán)系統(tǒng)回路如圖1所示,對(duì)除棒束通道外的上升段、下降段、冷凝器等管路和設(shè)備建立動(dòng)態(tài)數(shù)學(xué)模型[18],采用一維程序模擬回路內(nèi)冷卻劑流動(dòng)情況。

圖1 自然循環(huán)回路示意Fig.1 Schematic diagram of natural circulation loop

由于花瓣形燃料棒具有獨(dú)特的螺旋結(jié)構(gòu),棒束通道內(nèi)冷卻劑會(huì)產(chǎn)生明顯橫向流動(dòng)。為獲得準(zhǔn)確三維流場(chǎng)信息,本文采用CFD軟件模擬通道內(nèi)的流動(dòng)細(xì)節(jié),建立搖擺條件下一維程序與STAR-CCM+軟件相結(jié)合的多尺度耦合數(shù)值模型。

本文主要研究橫向搖擺對(duì)自然循環(huán)流動(dòng)特性的影響。以大地為基準(zhǔn),建立空間直角坐標(biāo)系O-XYZ,以Y軸為搖擺軸,將搖擺運(yùn)動(dòng)看作簡(jiǎn)諧運(yùn)動(dòng)。規(guī)定棒束通道上行為負(fù)角度,下行為正角度。則搖擺參數(shù)變化規(guī)律為:

(1)

(2)

(3)

式中:θ為搖擺角度,rad;θmax為最大搖擺角度,rad;ω為角速度,rad/s;β為角加速度,rad/s2;T為周期,s。

1.1 一維數(shù)值模型

本文采用控制容積法[19]建立搖擺條件下的動(dòng)量守恒方程為:

(4)

式中:li為控制體長(zhǎng)度,m;Ai為控制體流通截面,m2;G為質(zhì)量流量,kg/s;ΔPd為驅(qū)動(dòng)壓降,Pa;ΔPadd為搖擺附加壓降,Pa;ΔPz為總阻力壓降,Pa。

驅(qū)動(dòng)壓降主要由棒束通道加熱段、上升段與下降段流體密度差產(chǎn)生。設(shè)下降段流體密度為ρc,上升段流體密度為ρh,則瞬時(shí)驅(qū)動(dòng)壓降ΔPd為:

ΔPd=(ρc-ρh)gHcosθ

(5)

式中H為換熱高差,m。

此外,搖擺運(yùn)動(dòng)會(huì)引入附加慣性力,搖擺附加壓降ΔPadd為[20]:

β×r-2ω×ν]dl

(6)

式中:ω×(ω×r)為法向加速度;β×r為切向加速度;2ω×ν為科氏加速度。

自然循環(huán)回路總阻力壓降ΔPz,包括管路沿程阻力壓降ΔPz,l、局部阻力壓降ΔPz,p和棒束通道內(nèi)阻力壓降ΔPz,r?？傋枇航涤?jì)算公式[17]為:

ΔPz=ΔPz,l+ΔPz,p+ΔPz,r

(7)

式中ΔPz,r為棒束通道內(nèi)阻力壓降,由STAR-CCM+軟件計(jì)算得到,Pa。

綜上,式(4)可改寫(xiě)為常微分方程形式,式中驅(qū)動(dòng)壓降、附加壓降和總阻力壓降可由式(5)～(7)求得,并采用四階龍格庫(kù)塔法(RK4)求解[20]。

(8)

式中h為時(shí)間間隔,h=hn+1-hn,s。

1.2 螺旋花瓣型棒束通道模型

3×3螺旋花瓣型燃料棒束通道結(jié)構(gòu)和布置方式如圖2所示。其中,燃料棒橫截面外徑為D,內(nèi)凹弧半徑為R1,外凸弧半徑為R2,弧連接段長(zhǎng)度為hr,棒中心間距為L(zhǎng)p,棒束通道的矩形圍板邊長(zhǎng)為L(zhǎng)a,燃料棒直徑與內(nèi)凹弧半徑之比D/R1=5.1,與外凸弧半徑之比D/R2=10.3,圍板長(zhǎng)度與棒直徑之比La/D=3.2,棒間距與棒徑之比Lp/D=1.1,棒束通道總長(zhǎng)度Z=800 mm。

圖2 花瓣形燃料棒束幾何模型Fig.2 Geometric model of petal-shaped fuel rod

本文系統(tǒng)工質(zhì)為水,為準(zhǔn)確模擬花瓣形棒束通道流動(dòng)換熱特性,選取系統(tǒng)壓力p=0.5 MPa,溫度為298～418 K內(nèi)水的物性參數(shù),擬合后導(dǎo)入CFD軟件計(jì)算。將一維程序作為入口速度,設(shè)出口為壓力出口。為計(jì)算準(zhǔn)確,首先完成穩(wěn)態(tài)豎直條件下的多尺度耦合模擬,待收斂后,將該結(jié)果作為初始場(chǎng),進(jìn)行搖擺運(yùn)動(dòng)下的耦合模擬。目前,國(guó)內(nèi)外學(xué)者已采用CFD方法,開(kāi)展了螺旋花瓣型燃料棒束通道內(nèi)流動(dòng)換熱特性研究。結(jié)果表明,SSTk-ω湍流模型[21]能準(zhǔn)確模擬花瓣形燃料棒束通道熱工水力行為[10,12,17]。因此,本文采用該湍流模型開(kāi)展搖擺條件下多尺度耦合模擬研究。

1.3 數(shù)據(jù)交互方法

本文采用動(dòng)態(tài)鏈接庫(kù)實(shí)現(xiàn)一維程序與三維軟件間的數(shù)據(jù)交互[20]。交互參數(shù)為棒束通道入口溫度Tin、出口溫度Tout、進(jìn)出口壓降Pij、自然循環(huán)流速v和時(shí)間tn項(xiàng)。圖3給出了數(shù)據(jù)交互流程圖,具體如下:在STAR-CCM+軟件中讀取穩(wěn)態(tài)參數(shù),并調(diào)整加熱段搖擺姿態(tài),然后計(jì)算進(jìn)出口壓降和出口溫度值,完成后進(jìn)行收斂標(biāo)準(zhǔn)判斷,收斂標(biāo)準(zhǔn)為進(jìn)出口壓降的極差[20]。若通過(guò)判定,則將相關(guān)計(jì)算參數(shù)傳遞給一維程序,在一維程序中計(jì)算各分段壓降后,采用RK4法計(jì)算下一時(shí)間步的自然循環(huán)流速v,計(jì)算完成后采用動(dòng)態(tài)鏈接庫(kù)傳遞給CFD軟件,用于后續(xù)迭代計(jì)算。經(jīng)多次重復(fù)上述一維與三維數(shù)據(jù)交互過(guò)程,直至達(dá)到預(yù)期時(shí)間tpre,則完成全部搖擺條件下多尺度耦合模擬。

圖3 一維與三維數(shù)據(jù)交互流程圖Fig.3 Flow chart of 1-D and 3-D data interaction

2 多尺度耦合模型驗(yàn)證

本文選用文獻(xiàn)[22]中自然循環(huán)系統(tǒng)實(shí)驗(yàn)裝置和數(shù)據(jù)驗(yàn)證搖擺條件下多尺度耦合模型。該實(shí)驗(yàn)系統(tǒng)由搖擺裝置、測(cè)量系統(tǒng)、自然循環(huán)主回路和冷卻回路組成。

驗(yàn)證工質(zhì)為水,建立工質(zhì)物性與溫度的函數(shù)表達(dá)式,并將其加載到CFD軟件中,工況詳細(xì)參數(shù)設(shè)置參考文獻(xiàn)[16]。圖4給出了入口溫度75 ℃,熱流密度59 kW/m2,最大搖擺角度20°,搖擺周期30 s的實(shí)驗(yàn)結(jié)果與耦合模擬結(jié)果的入口流速與壓降對(duì)比圖。結(jié)果表明,模擬值與實(shí)驗(yàn)值吻合效果較好。在波峰和波谷處,入口流速最大相對(duì)誤差為7.19%,壓降最大相對(duì)誤差為1.16%。模擬結(jié)果與實(shí)驗(yàn)結(jié)果的平均流速相對(duì)誤差為0.16%,平均壓降相對(duì)誤差為0.47%。耦合模擬得到的流速和壓降結(jié)果均在可接受范圍內(nèi),這說(shuō)明本文建立的搖擺條件下一維與三維多尺度耦合模型可用于含花瓣形棒束通道自然循環(huán)系統(tǒng)熱工水力特性分析。

圖4 實(shí)驗(yàn)結(jié)果與耦合模擬結(jié)果驗(yàn)證Fig.4 Verification of coupling simulation results and experiment results

3 耦合模擬結(jié)果分析

基于搖擺條件下多尺度耦合數(shù)值模型,本文開(kāi)展3×3花瓣形燃料棒束通道自然循環(huán)系統(tǒng)流動(dòng)特性研究,熱工參數(shù)和搖擺參數(shù)如表1所示。

表1 工況參數(shù)設(shè)置Table 1 Setting of operating parameters

3.1 壓降波動(dòng)特性

3.1.1 壓降組成

搖擺運(yùn)動(dòng)會(huì)導(dǎo)致自然循環(huán)系統(tǒng)回路空間位置變化,還會(huì)引入附加慣性力,從而顯著影響自然循環(huán)系統(tǒng)流動(dòng)特性。為研究自然循環(huán)條件下棒束通道內(nèi)壓降波動(dòng)特性,取棒束中間600 mm作為測(cè)壓段,將測(cè)壓段總壓降ΔPt,分離為測(cè)壓段搖擺附加壓降ΔPadd,t、流動(dòng)阻力壓降ΔPz,t和重位壓降ΔPg。圖5給出了加熱功率10 kW,入口溫度40 ℃,最大搖擺角度10°,搖擺周期20 s條件下各壓降和體積流量Qv隨無(wú)量綱時(shí)間(t/T)變化曲線。

圖5 壓降和流量變化曲線Fig.5 Variation curves of pressure drop and flow rate under rolling motion condition

可以看到,測(cè)壓段總壓降ΔPt與搖擺角度θ的周期相同,同一周期內(nèi)測(cè)壓段總壓降ΔPt相鄰波峰幅值不同,這是各組成壓降共同影響的結(jié)果。測(cè)壓段重位壓降ΔPg占比最大,因此與總壓降ΔPt波動(dòng)規(guī)律相似。在一個(gè)周期內(nèi),棒束通道2次運(yùn)動(dòng)至與地面垂直位置,此刻重位壓降達(dá)到最大值,而測(cè)壓段搖擺附加壓降ΔPadd,t為零,搖擺運(yùn)動(dòng)使豎直姿態(tài)下流體密度不同,從而導(dǎo)致測(cè)壓段重位壓降相鄰波峰輻值不同[23]。

3.1.2 流量和阻力壓降相位特性

圖6給出了不同加熱功率下自然循環(huán)流量G和測(cè)壓段阻力壓降ΔPz,t相位關(guān)系曲線。結(jié)果表明,流量和阻力壓降曲線的相位差并不明顯,最大不超過(guò)搖擺周期的4.5%,與自然循環(huán)系統(tǒng)中圓柱形棒束通道[23]和矩形管道[24]的研究結(jié)果相同。此外,波谷位置處的相位差要比波峰位置的更大,這是由于波峰、波谷處的流體溫度相差較大,流體物性變化顯著,致使驅(qū)動(dòng)力變化不一致,以及棒束通道阻力系數(shù)改變而導(dǎo)致的。

圖6 不同加熱功率下流量和阻力壓降相位變化Fig.6 Phase shift between the flow rate and resistance pressure drop

3.2 加熱功率影響特性

3.2.1 加熱功率對(duì)自然循環(huán)流量的影響

自然循環(huán)流量是自然循環(huán)系統(tǒng)研究的重要參數(shù)之一。為便于定量評(píng)價(jià)曲線波動(dòng)特性,定義體積流量相對(duì)值(QV)re為:

(QV)re=QV/(QV)av

(9)

式中:QV為瞬時(shí)體積流量,m3/h;(QV)av為搖擺條件下時(shí)均體積流量,m3/h。

圖7所示為不同功率下,流量波動(dòng)曲線。加熱功率升高,自然循環(huán)流量值增加,波動(dòng)幅度減小,且功率從5 kW增加到9 kW,波谷區(qū)流量的變化量較波峰區(qū)提高19.81%。

圖7 不同功率下流量波動(dòng)曲線Fig.7 Fluctuation curves of flow rate at different powers

由式(5)可知,在搖擺條件下,當(dāng)有效高差不變,驅(qū)動(dòng)壓降由密度差決定。提高燃料棒功率,流體吸熱量增加,溫差增加,密度差增大,驅(qū)動(dòng)壓降增大,流量提高,波動(dòng)幅度減小,系統(tǒng)更穩(wěn)定。此外,由于阻力壓降與流量的平方滿足正比關(guān)系,高流量下阻力壓降變化比低流量的更顯著。在相同驅(qū)動(dòng)壓降變化時(shí),低流量條件下,驅(qū)動(dòng)壓降與阻力壓降滿足準(zhǔn)穩(wěn)態(tài)匹配所需流量的變化量越大[25]。因此,波谷區(qū)流量的變化量較波峰區(qū)的更大。

3.2.2 加熱功率對(duì)阻力系數(shù)的影響

為進(jìn)一步分析螺旋花瓣型棒束通道內(nèi)流動(dòng)阻力特性,圖8給出了不同功率下棒束通道內(nèi)流動(dòng)阻力系數(shù)隨時(shí)間變化曲線。其中,阻力系數(shù)可通過(guò)達(dá)西公式求得[26]。結(jié)果表明,隨燃料棒加熱功率升高,棒束通道內(nèi)阻力系數(shù)降低,并趨于平緩。

圖8 不同加熱功率下阻力系數(shù)波動(dòng)曲線Fig.8 Fluctuation of resistance coefficients of fuel bundle channels for different heating powers

圖9給出了τ=0.75時(shí)不同燃料棒加熱功率下z/Z=0.5處粘度分布。結(jié)果表明,功率從5 kW提高到9 kW,流體溫度升高,近壁面流體平均粘度從3.32×10-4Pa·s降至2.91×10-4Pa·s,粘度降低。壁面粘滯力減小,阻力系數(shù)隨之降低。

圖9 不同加熱功率下棒束通道z/Z=0.5處粘度分布Fig.9 Viscosity distribution at z/Z=0.5 of fuel bundle channel for different heating power

3.3 阻力系數(shù)關(guān)聯(lián)式

搖擺條件下花瓣形燃料棒束通道內(nèi)流動(dòng)阻力特性除考慮流動(dòng)參數(shù)和熱工參數(shù)的影響外,還需考慮搖擺角度θ和搖擺周期T等搖擺參數(shù)影響。因此,描述搖擺條件下棒束通道內(nèi)阻力特性的函數(shù)關(guān)系式為:

F(ΔPf,ρ,v,r,L,De,μ,ω,β,Tb,Tw)=0

(10)

式中:r為棒束通道中心到搖擺軸的距離,m;L為測(cè)壓段長(zhǎng)度,m;De為當(dāng)量直徑,m;μ為動(dòng)力粘度,Pa·s;Tw和Tb分別為棒表面和主流平均溫度,K。無(wú)量綱化后可以得到Re、L/De、Tw/Tb、ωr/v和βrDe/v2等無(wú)量綱數(shù)。其中,ωr/v表示切向慣性力與流動(dòng)慣性力之比,βrDe/v2表示法向慣性力與流動(dòng)慣性力之比。流動(dòng)阻力壓降為:

(11)

相較于靜止條件,搖擺條件下棒束通道內(nèi)流動(dòng)阻力系數(shù)λ會(huì)受流體自身慣性力、黏性力,以及搖擺附加慣性力影響:

λ=f(Re,ωr/v,βrDe/v2,Tw/Tb)

(12)

搖擺條件下圓柱形棒束通道流動(dòng)阻力系數(shù)經(jīng)驗(yàn)表達(dá)式一般形式[27]為:

(13)

式中:μb為主流動(dòng)力粘度,Pa·s;μw為燃料棒近壁面流體動(dòng)力粘度,Pa·s;μw/μb用于考慮Tw/Tb的影響。

此外,流體因搖擺運(yùn)動(dòng)而受到搖擺附加慣性力作用。因此,通過(guò)引入瞬時(shí)搖擺雷諾數(shù)[23,27],準(zhǔn)確描述自然循環(huán)系統(tǒng)的瞬態(tài)特性為:

Reroll=ρωrDe/μ

(14)

式中Reroll為搖擺雷諾數(shù)。

(15)

考慮到流道進(jìn)出口效應(yīng),選取花瓣形棒束通道中間600 mm長(zhǎng)度進(jìn)行研究,對(duì)于搖擺角度15°～25°,搖擺周期20～30 s,瞬時(shí)雷諾數(shù)800～2 000內(nèi)的2 000個(gè)瞬時(shí)數(shù)據(jù)進(jìn)行擬合,計(jì)算得到式(13)中的系數(shù)m=-1.856,以及系數(shù)bi的計(jì)算表達(dá)式為:

(16)

圖10給出了瞬時(shí)摩擦阻力系數(shù)耦合模擬計(jì)算值和預(yù)測(cè)值的對(duì)比結(jié)果。

圖10 瞬時(shí)阻力系數(shù)耦合值λn與預(yù)測(cè)值λp對(duì)比Fig.10 Comparison of resistance coefficient between coupling simulation results λn and predicted values λp

可以看到,計(jì)算的相對(duì)誤差集中在-5%～7%,本文所得關(guān)聯(lián)式在低雷諾數(shù)范圍內(nèi)能較好地預(yù)測(cè)搖擺條件下瞬時(shí)阻力系數(shù)。

4 結(jié)論

1)自然循環(huán)條件下,流動(dòng)阻力壓降與自然循環(huán)流量相位差不明顯。波峰、波谷處流體溫度變化,導(dǎo)致波谷處的相位差大于波峰處的數(shù)值。

2)加熱功率升高,會(huì)使自然循環(huán)系統(tǒng)更加穩(wěn)定。驅(qū)動(dòng)壓降與阻力壓降須滿足準(zhǔn)穩(wěn)態(tài)匹配,波谷區(qū)流量的下降量較波峰區(qū)流量上升量更大。

3)通過(guò)引入瞬時(shí)搖擺雷諾數(shù),獲得適用于搖擺條件下自然循環(huán)系統(tǒng)流動(dòng)阻力系數(shù)預(yù)測(cè)式。結(jié)果表明,本文模擬值與預(yù)測(cè)值相對(duì)誤差集中在-5%～7%。對(duì)不同類型花瓣形棒束通道,仍需分析結(jié)構(gòu)參數(shù)對(duì)阻力系數(shù)的影響。