劉 凱,王明軍,章 靜,田文喜,秋穗正,蘇光輝
(西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院,陜西 西安 710049)
蒸汽發(fā)生器是核動力系統(tǒng)一、二回路能量傳遞樞紐,作為核反應(yīng)堆一回路壓力邊界,是一回路帶放射性冷卻劑與二回路工質(zhì)的隔離屏障。蒸汽發(fā)生器結(jié)構(gòu)龐雜,且涉及熱工水力現(xiàn)象復(fù)雜多變,因此其運行可靠性對核電廠安全運行具有重要意義[1]。
螺旋管蒸汽發(fā)生器具有結(jié)構(gòu)緊湊、換熱效率高、抗膨脹熱應(yīng)力強等優(yōu)勢,在能源動力、石油化工等領(lǐng)域得到廣泛應(yīng)用。國內(nèi)外學(xué)者針對管側(cè)和殼側(cè)流動換熱特性開展了大量實驗研究。Messa等[2]和Genic等[3]分別采用多個具有不同節(jié)距、螺升角等結(jié)構(gòu)參數(shù)的螺旋管換熱器實驗段開展了冷熱流體逆流換熱實驗,獲得了殼側(cè)對流換熱系數(shù)關(guān)系式;Zhao等[4]在一定系統(tǒng)壓力、質(zhì)量流量和加熱功率范圍內(nèi)開展了一系列螺旋管內(nèi)兩相流動沸騰實驗,獲得了阻力和換熱系數(shù)關(guān)系式;Hwang等[5]進行了具有不同螺旋直徑螺旋管內(nèi)干涸特性實驗研究,結(jié)果表明在低質(zhì)量流量下,二次流對干涸點含氣率影響較大;白博峰等[6]研究了螺旋管蒸汽發(fā)生器在同步啟動過程中的流動和傳熱特性,瞬態(tài)單相湍流傳熱和臨界熱負荷規(guī)律與穩(wěn)態(tài)有明顯不同;畢勤成等[7]開展了高壓下高溫氣冷堆蒸汽發(fā)生器螺旋管內(nèi)兩相流動實驗,得到了螺旋管內(nèi)兩相摩擦阻力關(guān)系式。
近年來,計算流體力學(xué)(CFD)方法在核動力系統(tǒng)安全分析中得到廣泛應(yīng)用[8-11]。Zhao等[12-13]針對立式U型管蒸汽發(fā)生器開展了正常運行及傳熱管堵塞事故工況下全尺寸熱工水力特性分析,獲得了液相溫度、空泡份額等一二次側(cè)關(guān)鍵參數(shù)分布;Mu等[14]實現(xiàn)了對長期運行條件下蒸汽發(fā)生器沉積物分布的預(yù)測,并分析和評估了沉積物污垢熱阻對蒸發(fā)器傳熱性能的影響;He等[15]基于OpenFOAM開發(fā)了適用于管殼式換熱器的多孔介質(zhì)方法求解器,結(jié)合MB-2蒸汽發(fā)生器基準題開展了模型驗證,計算結(jié)果與實驗數(shù)據(jù)符合較好。針對高溫氣冷堆(HTGR)螺旋管蒸汽發(fā)生器,Li等[16]模擬了管束間流體熱攪混過程,獲得了不同流量下貝克萊數(shù)(Pe)沿流動方向的變化規(guī)律;Ma等[17]采用一維程序分析了啟動過程及非均勻加熱條件下兩相流動不穩(wěn)定性;Olson等[18]采用二維多孔介質(zhì)方法開展了穩(wěn)態(tài)熱工水力計算,并研究了流量變化、換熱管位移和堵塞對溫度分布的影響。目前以螺旋管蒸汽發(fā)生器整體為對象的相關(guān)研究較少,因此開發(fā)其全尺寸三維熱工水力特性分析程序并開展設(shè)備級數(shù)值模擬,對于螺旋管蒸汽發(fā)生器設(shè)計優(yōu)化及安全評價具有重要意義。
開源CFD平臺OpenFOAM具有編程環(huán)境開放特性,便于進行數(shù)據(jù)接口創(chuàng)建、模型修改植入及求解器編寫,更好地滿足用戶自主開發(fā)需求,已在航空航天、化工過程等領(lǐng)域得到廣泛應(yīng)用。本文采用多孔介質(zhì)方法對具有復(fù)雜多層螺旋套管結(jié)構(gòu)的換熱組件區(qū)域進行簡化,構(gòu)建殼側(cè)流動換熱特性數(shù)學(xué)物理模型,并建立管側(cè)水-水蒸氣兩相流動沸騰換熱特性分析模型,采用網(wǎng)格-節(jié)點數(shù)據(jù)映射方法實現(xiàn)管殼兩側(cè)耦合傳熱計算,基于OpenFOAM平臺開發(fā)適用于螺旋管蒸汽發(fā)生器的三維全尺寸熱工水力特性分析程序HeTAF。
螺旋管直流式蒸汽發(fā)生器多采用套筒式組件結(jié)構(gòu),由外套筒、中心筒及二者之間環(huán)形腔體內(nèi)多層反向纏繞螺旋套管組成,如圖1所示。管殼兩側(cè)工質(zhì)逆向流動,換熱組件殼側(cè)單相高溫工質(zhì)自上而下沖刷螺旋管束,管側(cè)單相過冷水在螺旋管內(nèi)自下而上流動,與殼側(cè)流體換熱后發(fā)生相變,依次經(jīng)過泡狀流、彈狀流、環(huán)狀流和滴狀流后轉(zhuǎn)變?yōu)閱蜗嗾羝?,存在?fù)雜的兩相流動沸騰換熱現(xiàn)象。
圖1 換熱組件幾何結(jié)構(gòu)示意圖Fig.1 Schematic diagram of heat exchange assembly geometric structure
本文基于多孔介質(zhì)方法對換熱組件內(nèi)復(fù)雜螺旋管結(jié)構(gòu)進行簡化,建立了殼側(cè)工質(zhì)流動換熱特性數(shù)學(xué)物理模型,其控制方程如下。
1) 質(zhì)量守恒方程:
(1)
2) 動量守恒方程:
(2)
3) 能量守恒方程:
(3)
其中:下標s指殼側(cè)流體相關(guān)參數(shù);β為換熱組件內(nèi)殼側(cè)流域孔隙率;ρ為密度;U為流速;p為壓力;h為比焓;g為重力加速度;μ為動力黏性系數(shù);α為熱擴散系數(shù);SM為由螺旋管束引入的阻力源項(忽略外套筒與中心筒的表面摩擦);SE為殼側(cè)工質(zhì)被低溫螺旋管壁冷卻而引入的對流換熱源項。
SM表達式為:
(4)
其中:L為流動方向長度;Af為單位體積內(nèi)流通面積;fs為殼側(cè)工質(zhì)橫掠螺旋管束摩擦阻力系數(shù),由Idelchik關(guān)系式[19]計算;De為水力學(xué)直徑,可由殼側(cè)流通體積與潤濕面積等效計算。
De表達式為:
(5)
fs表達式為:
(6)
其中:S1和S2分別為相鄰螺旋管橫向間距和縱向間距;do為螺旋管外徑。
SE表達式為:
SE=-hs(Ts-Tw,o)At
(7)
其中:Tw,o為螺旋管外壁面溫度;At為單位體積內(nèi)螺旋管外壁面面積;hs為螺旋管表面對流換熱系數(shù),根據(jù)Zukauskas關(guān)系式[20]計算,其表達式為:
(8)
其中:Prw為定性溫度取壁面溫度的普朗特數(shù);λ為熱導(dǎo)率。
均相流模型將兩相混合物流動視為具有特殊流動參數(shù)的均勻單相流動,不考慮相間相互作用,其假想物性由氣液飽和狀態(tài)參數(shù)加權(quán)平均獲得,具有模型形式簡易、可靠性強等優(yōu)勢。相關(guān)研究[21]表明,螺旋管中氣液兩相流動基本處于熱力學(xué)平衡狀態(tài),且滑速比整體較小,均相流模型在螺旋管蒸汽發(fā)生器熱工水力特性分析中的適用性和計算精度已得到廣泛證明,因此本文采用均相流模型描述螺旋管內(nèi)兩相流動沸騰換熱過程。由于管內(nèi)流動Pe較大,對流作用占主導(dǎo),因此忽略動量及能量擴散[22],則控制方程如下。
1) 質(zhì)量守恒方程:
(9)
2) 動量守恒方程:
(10)
3) 能量守恒方程:
(11)
其中:下標m指管側(cè)均相流相關(guān)參數(shù);Gm為質(zhì)量流密度;qw,i為管內(nèi)壁熱流密度;Ci為螺旋管內(nèi)截面周長;Ai為螺旋管內(nèi)截面積;fm為摩擦阻力系數(shù)。
均相流物性參數(shù)通過下式求解:
(12)
(13)
ρm=αρg+(1-α)ρl
(14)
(15)
λm=xeλg+(1-xe)λl
(16)
其中:下標l和g分別指飽和水和飽和蒸汽的相關(guān)參數(shù);xe為平衡態(tài)含氣率。
在螺旋管換熱器中,過冷水沿管側(cè)流動與管壁換熱并發(fā)生流動沸騰,隨著空泡份額的上升,兩相流型逐步轉(zhuǎn)變,最終液相被完全蒸干,全部轉(zhuǎn)變?yōu)閱蜗噙^熱蒸汽。本文在計算中根據(jù)換熱機理,將管側(cè)流動分為如下區(qū)域:單相水對流區(qū)、欠熱沸騰區(qū)、飽和沸騰與強迫對流蒸發(fā)區(qū)、干涸缺液區(qū)和單相蒸汽對流區(qū),各區(qū)域間換熱機理轉(zhuǎn)變的判斷標準列于表1。
對于單相水和單相蒸汽對流區(qū),選用廣泛適用于管內(nèi)流動換熱的Schmidt關(guān)系式[24],其表達式為:
Num=
(17)
表1 管側(cè)換熱機理轉(zhuǎn)變判斷標準Table 1 Criterion for transition of tube side heat transfer mechanism
其中:di為螺旋管內(nèi)徑;Dc為平均螺旋直徑;Nu、Re和Pr分別為管側(cè)流體努塞爾數(shù)、雷諾數(shù)和普朗特數(shù);常數(shù)Re1和Re2分別為22 000和150 000;Recr為臨界雷諾數(shù),本文采用Ito關(guān)系式[25]計算。
在飽和沸騰與強迫對流蒸發(fā)區(qū),隨著相變的進行,汽泡在主流聚集并形成高速流動的蒸汽核心。該區(qū)域換熱系數(shù)采用Chen關(guān)系式[26]計算,其表達式為:
(18)
其中:cp為比熱容;ΔTsat為飽和溫度與當(dāng)前溫度差值;Δpsat為飽和壓力與當(dāng)前壓力差值。
在欠熱沸騰區(qū),加熱面產(chǎn)生汽泡躍離成核點后被主流冷凝,對該區(qū)域采用基于飽和沸騰關(guān)系式(式(17))修正后的Zhao關(guān)系式[27]計算,即:
(19)
對于干涸缺液區(qū),近壁面液膜受加熱蒸干及汽泡撕扯作用,逐漸減薄直至破壞,發(fā)生“干涸”。該區(qū)域換熱系數(shù)采用Miropolskiy關(guān)系式[28]計算,其表達式為:
(20)
對于管內(nèi)流體摩擦阻力系數(shù),單相區(qū)內(nèi)采用Ito關(guān)系式[25]計算:
(21)
其中,fs為層流摩擦阻力系數(shù),其表達式為:
(22)
(23)
其中,Xtt為馬蒂內(nèi)利參數(shù)。
上述模型適用于0.001 16
為實現(xiàn)三維多孔介質(zhì)模型與一維均相流模型參數(shù)傳遞,開展了殼側(cè)工質(zhì)與管側(cè)水-水蒸氣兩相流耦合換熱求解,建立了三維計算網(wǎng)格與一維計算節(jié)點間數(shù)據(jù)映射關(guān)系,如圖2所示。
圖2 網(wǎng)格-節(jié)點數(shù)據(jù)映射示意圖Fig.2 Schematic diagram of grid-node data mapping
沿組件高度方向?qū)⒙菪苁鴦澐譃镹個離散節(jié)點作為均相流模型計算載體,則集總于各節(jié)點的螺旋管長度Ln、換熱面積Sn(按外壁面面積計算)和管側(cè)流體體積Vn可分別表示為:
(24)
Sn=πdoLn
(25)
(26)
其中:Nt為組件內(nèi)螺旋管總數(shù);Lt為單根螺旋管長度。
同樣地,將殼側(cè)流體域沿組件高度方向劃分為N層,與離散節(jié)點逐個對應(yīng),則任一管側(cè)節(jié)點P與其高度范圍內(nèi)殼側(cè)流體域P′存在耦合換熱關(guān)系。通過對區(qū)域P′內(nèi)控制體網(wǎng)格進行積分,得到體積平均的殼側(cè)氦氣溫度TP′和對流換熱系數(shù)hP′,將其作為邊界條件用于均相流模型計算,其表達式分別為:
(27)
(28)
其中:NP′為區(qū)域P′內(nèi)控制體網(wǎng)格總數(shù);Vc,i為控制體網(wǎng)格i的體積;hc,i和hl,i為控制體網(wǎng)格i中管側(cè)工質(zhì)溫度和對流換熱系數(shù)。
忽略螺旋管表面沉積污垢熱阻,則兩側(cè)總換熱系數(shù)可表示為:
(29)
其中,λw為螺旋管熱導(dǎo)率。該節(jié)點處換熱量可表示為:
Qn=qw,oSn=htotal(Ts-Tm)Sn
(30)
忽略螺旋管壁熱慣性,則螺旋管外壁面溫度為:
(31)
綜上,管殼兩側(cè)耦合換熱求解流程如下:1) 根據(jù)假設(shè)(或更新后)的對流換熱源項求解殼側(cè)工質(zhì)控制方程,獲得殼側(cè)流場及溫度場分布,進而求得螺旋管表面對流換熱系數(shù);2) 根據(jù)網(wǎng)格-節(jié)點數(shù)據(jù)映射關(guān)系,將殼側(cè)流體域網(wǎng)格氦氣溫度和對流換熱系數(shù)體積平均積分,并作為邊界條件傳遞至管側(cè)一維均相流模型對應(yīng)節(jié)點中;3) 求解管側(cè)水-水蒸氣控制方程,獲得管側(cè)流場及溫度場分布,并重新計算兩相流物性、沸騰換熱系數(shù)等流動換熱相關(guān)參數(shù);4) 將步驟3中更新的各節(jié)點的螺旋管外壁面溫度傳遞至殼側(cè)流體域網(wǎng)格,計算殼側(cè)工質(zhì)的對流換熱源項并作為步驟1的初始條件。按上述步驟反復(fù)迭代直至得到收斂的計算結(jié)果。
基于開源CFD平臺OpenFOAM中的有限體積類庫建立殼側(cè)工質(zhì)流動換熱模型,并植入水-水蒸氣兩相流動沸騰換熱模型及相應(yīng)熱物性模型,通過網(wǎng)格-節(jié)點映射方法實現(xiàn)熱工水力參數(shù)傳遞及管殼兩側(cè)耦合換熱模型的求解,開發(fā)了適用于螺旋管蒸汽發(fā)生器的三維全尺寸熱工水力特性分析程序HeTAF。
基于意大利SIET熱工水力實驗室Santini[30]螺旋管兩相流動沸騰換熱實驗進行HeTAF模型驗證,實驗裝置及實驗段如圖3所示。通過旁路管線與泵下游控制閥調(diào)整實驗段流量,并設(shè)置節(jié)流閥以抑制密度波流動不穩(wěn)定性,實驗段上游預(yù)熱器用于調(diào)節(jié)進口過冷水溫度。
該實驗采用開放式回路設(shè)計,工質(zhì)為去離子水。實驗段為一根豎直布置的螺旋管,其主要幾何參數(shù)列于表2。通過直流電源對前24 m
表2 實驗段幾何參數(shù)Table 2 Geometric parameter of test section
實驗段進行均勻加熱以使管內(nèi)流體發(fā)生沸騰,剩余螺旋管為絕熱段。在不同軸向位置測溫點處沿螺旋管外圍均勻布置多個K型熱電偶測量螺旋管外壁溫。
對系統(tǒng)壓力6 MPa下具有不同質(zhì)量流速G與熱流密度q的各實驗工況進行模擬,得到管內(nèi)流體與壁面換熱系數(shù)在實驗段內(nèi)的分布,為使計算結(jié)果與實驗記錄數(shù)據(jù)方式一致,將管長轉(zhuǎn)化為平衡態(tài)含氣率,其對比如圖4所示,其中實驗數(shù)據(jù)設(shè)置20%誤差棒。
圖4 換熱系數(shù)隨含氣率變化的計算值與實驗值Fig.4 Measured and calculated heat transfer coefficient variation with equilibrium quantity
由圖4可見,除首個實驗工況數(shù)據(jù)序列中個別數(shù)據(jù)點外,其他計算誤差均在20%以內(nèi),3個實驗工況下平均計算誤差分別為10.1%、6.8%和6.5%,即在高質(zhì)量流速和熱流密度條件下預(yù)測精度更高。由此可見,本文所采用模型能有效預(yù)測螺旋管實驗段中兩相流動沸騰換熱特性。
本文以高溫氣冷堆示范工程[31-32]中的螺旋管直流式蒸汽發(fā)生器(HCOTSG)為對象,利用開發(fā)的自主化螺旋管蒸汽發(fā)生器三維熱工水力程序HeTAF開展數(shù)值模擬,高溫堆HCOTSG結(jié)構(gòu)參數(shù)列于表3。
表3 HCOTSG幾何參數(shù)Table 3 Geometric parameter of HCOTSG
針對螺旋管蒸汽發(fā)生器中單個換熱組件開展兩側(cè)熱工水力特性耦合分析,每個換熱單元有5層,共35根螺旋換熱管,從里向外每層螺旋管根數(shù)依次為5、6、7、8、9,相鄰兩層纏繞方向相反。各層螺旋管間螺升角和螺距不同,使得換熱組件內(nèi)所有螺旋管管長基本保持一致。研究[33-34]表明,在螺旋管蒸汽發(fā)生器實際工況下,熱工水力特性對螺旋直徑小幅變化不敏感,因此在模擬中采用平均螺旋直徑430 mm作為特征螺旋直徑求解相關(guān)參數(shù),即采用特征管表征不同層螺旋管。沿組件高度方向劃分150個計算節(jié)點,兩側(cè)相關(guān)參數(shù)列于表4。
表4 管殼兩側(cè)流動相關(guān)參數(shù)Table 4 Parameter of shell side and tube side fluid
模擬時計算域上下邊界分別為氦氣進出口邊界,忽略換熱組件對外輻射散熱等其他熱損耗,外套筒內(nèi)壁面及中心筒外壁面設(shè)置為絕熱壁面邊界。氦氣、單相過冷水及過熱蒸汽熱物性則由根據(jù)組件運行范圍內(nèi)多個溫度下物性點擬合而成的多項式計算得到,處于兩相區(qū)的汽液兩相混合物的熱物性則由對應(yīng)工況下飽和水及飽和蒸汽熱物性加權(quán)獲得,如1.2節(jié)所述。殼側(cè)三維計算網(wǎng)格和管側(cè)一維計算節(jié)點上物理場和相關(guān)參數(shù)分別按熱氦氣和過冷水進口溫度及流量初始化。管側(cè)和殼側(cè)流體及管外壁面溫度沿組件高度方向分布計算結(jié)果如圖5、6所示。
圖5 管側(cè)和殼側(cè)流體及管外壁面溫度沿組件高度方向的分布Fig.5 Distribution of tube side fluid, shell side fluid and tube outer wall temperature along height direction
圖6 組件中管側(cè)和殼側(cè)流體溫度分布Fig.6 Distribution of tube side and shell side fluid temperature in assembly
由圖5可見,管側(cè)流體存在單相水區(qū)、兩相區(qū)、單相蒸汽區(qū)等3區(qū)分布,在螺旋管長約35 m處(對應(yīng)組件高度約5 m處)發(fā)生飽和沸騰,呈明顯溫度平臺,如圖6a所示,在螺旋管長約52 m處(對應(yīng)組件高度約7.4 m處)液相全部蒸干,管側(cè)流體變?yōu)閱蜗噙^熱蒸汽,兩相區(qū)長度約為17 m。相應(yīng)地,殼側(cè)氦氣溫度在組件上部冷卻速度較快,而由于兩側(cè)換熱能力下降,在組件下部溫度減小趨勢放緩,如圖6b所示??梢?,本文模型計算獲得的兩側(cè)流體溫度變化趨勢與實際物理過程較符合。殼側(cè)氦氣和管側(cè)蒸汽出口溫度計算誤差列于表5,兩側(cè)換熱量計算值相比設(shè)計值略低。
表5 組件出口溫度計算誤差Table 5 Calculation error of assembly outlet temperature
組件內(nèi)殼側(cè)及管側(cè)流體換熱系數(shù)、空泡份額和含氣率沿螺旋管長度方向的分布如圖7、8所示。在沸騰起始點前,管側(cè)單相水與加熱管壁對流換熱;在進入兩相區(qū)后,管側(cè)流體沿螺旋管流動過程中相變逐漸劇烈,其換熱機理隨組件高度的增加而不斷變化,依次經(jīng)歷欠熱沸騰、飽和沸騰、強迫對流蒸發(fā)換熱,管側(cè)換熱系數(shù)快速升高,最大值約為6.21×104W/(m2·K)。殼側(cè)為氦氣橫掠螺旋管束的單相對流換熱,換熱系數(shù)相較于管側(cè)在整個計算域內(nèi)變化不大,其平均值約為2 383.57 W/(m2·K)。
圖7 殼側(cè)及管側(cè)換熱系數(shù)、空泡份額和含氣率沿螺旋管長分布Fig.7 Distribution of tube side heat transfer coefficient,void fraction and equilibrium quantity along tube length
圖8 組件中管側(cè)換熱系數(shù)和空泡份額分布Fig.8 Distribution of tube side heat transfer coefficient and fluid void fraction in assembly
由于空泡份額不斷升高,導(dǎo)致兩相流型改變,最終環(huán)狀流液膜被破壞,壁面與氣相直接接觸并進入干涸缺液區(qū),管側(cè)換熱系數(shù)驟降至2.5×104W/(m2·K);最終液相完全蒸干,受熱物性變化影響,單相蒸汽區(qū)中管側(cè)換熱系數(shù)有所下降,但由于蒸汽與管內(nèi)壁面之間溫差較大,換熱功率仍高于單相水對流換熱區(qū)??张莘蓊~在兩相區(qū)前段上升較快,在后半段則趨于緩慢。
本文基于多孔介質(zhì)方法對換熱組件內(nèi)復(fù)雜螺旋管結(jié)構(gòu)進行簡化,建立了殼側(cè)工質(zhì)流動換熱特性分析模型,針對管側(cè)兩相流動沸騰換熱過程建立了水-水蒸氣均相流模型,根據(jù)換熱組件的空間對應(yīng)關(guān)系,采用網(wǎng)格-節(jié)點數(shù)據(jù)映射方法進行模型間數(shù)據(jù)傳遞,實現(xiàn)了管殼兩側(cè)耦合傳熱計算,開發(fā)了適用于螺旋管蒸汽發(fā)生器的三維全尺寸熱工水力特性分析程序HeTAF。基于螺旋管兩相流動沸騰換熱實驗開展了模型驗證;對高溫氣冷堆示范工程中螺旋管直流式蒸汽發(fā)生器的單個換熱組件開展了兩側(cè)熱工水力特性耦合模擬,結(jié)果表明:計算得到的氦氣和蒸汽出口溫度與設(shè)計值符合較好,絕對誤差分別為5.29 K和5.76 K,表明本文所建立模型能有效預(yù)測換熱組件內(nèi)管殼兩側(cè)耦合流動換熱特性。
未來將以單換熱組件模擬策略為基礎(chǔ),對螺旋管蒸汽發(fā)生器開展全尺寸熱工水力特性耦合分析,并針對HeTAF程序持續(xù)開展模型驗證完善和應(yīng)用范圍推廣等研究。