馬競翔 董世昌 龔圣捷
(上海交通大學(xué) 上海 200240)
壓水堆(Pressurized Water Reactor,PWR)核電站一回路以水為冷卻介質(zhì)[1]。部分與冷卻劑系統(tǒng)主管道相連接的分支管內(nèi)流體處于滯流狀態(tài),分支管內(nèi)的流體溫度低而主管道內(nèi)的流體溫度高[2]。由于湍流滲透或閥門泄漏等原因,主管道內(nèi)的高溫流體進入分支管內(nèi),從而形成熱分層現(xiàn)象[3]。
在水平或傾斜管道中,熱流體緩慢流入冷流體管段時,由于溫度不同造成的密度差,熱流體浮于冷流體的上層,在重力方向上形成溫度梯度,即熱分層現(xiàn)象[4]。熱分層會導(dǎo)致管道壁面應(yīng)力分布不均勻,并引起管道熱疲勞及管道結(jié)構(gòu)失效[5]。
自20世紀80年代起,壓水堆核電站管路的熱分層現(xiàn)象被發(fā)現(xiàn)后就受到了廣泛關(guān)注。美國核管理委員會(U.S.Nuclear Regulatory Commission,NRC)發(fā)布公告要求所有在役或在建核電廠必須對管路熱分層現(xiàn)象進行分析論證,以確保管道結(jié)構(gòu)的完整性[6-7]。根據(jù)研究,壓水堆核電站冷卻系統(tǒng)的多處結(jié)構(gòu)會受到熱分層現(xiàn)象的影響,例如:蒸汽發(fā)生器的給水管線、主冷卻系統(tǒng)熱管段、高壓安注管線、余熱排出管線等部位[8-9]。
賀斌[10]對核電站中蒸汽發(fā)生器上水管路熱分層的流動與傳熱機理進行了數(shù)值模擬研究。發(fā)現(xiàn)二次流和浮升力的作用會使得冷熱流體交界面不斷變化,并驗證了冷熱流體的相對流速和相對溫差會影響流場的溫度變化趨勢。曹瓊等[11]對AP1000 第四級卸壓A管線ADS-4A的熱分層開展穩(wěn)態(tài)數(shù)值模擬研究。發(fā)現(xiàn)保溫層的設(shè)置使得管道沿程的溫度梯度減小;ADS-4A管線的45°彎頭能夠有效降低由湍流滲透作用進入管線中的熱流體的流速。Jo等[4]采用商用計算流體動力學(xué)(Computational Fluid Dynamics,CFD)軟件對穩(wěn)壓器波動管進行了波入和波出兩種工況下的數(shù)值模擬工作,提出了管壁的影響,認為計算模型包含管壁并采用共軛傳熱的方式,能更準確地預(yù)測管道內(nèi)的熱分層現(xiàn)象。Cai 等[12]基于Richardson數(shù)(Ri)相似,建立了1/3比例的穩(wěn)壓器波動管模型,開展了波入和波出兩種工況條件下的模擬實驗。通過實驗數(shù)據(jù)與模擬結(jié)果對比發(fā)現(xiàn):SSTk-ω湍流模型的精度要優(yōu)于標準k-ε湍流模型(Standardk-εModel,SKE)和雷諾應(yīng)力模型(Reynolds Stress Model,RSM)。
目前,熱分層現(xiàn)象的研究主要圍繞等截面管道模型開展,而對于變截面管道的瞬態(tài)流動特性和溫度變化特性的研究較少。本研究建立了存在管徑變化的滯流分支管模型,開展了數(shù)值模擬工作,分析其在閥門泄漏情況下發(fā)生熱分層現(xiàn)象時的溫度變化特性和管道內(nèi)的流場特性,為后續(xù)的模擬實驗和管道應(yīng)力分析提供理論依據(jù)。
滯流分支管采用管徑分別為DN100、DN200、DN350的直管段、大小頭以及彎管段組成,管道只有部分位置布置有保溫層。為了模擬滯流分支管因入口處閥門泄漏而引發(fā)熱分層現(xiàn)象,管道入口位置采用封板結(jié)構(gòu)與DN15的管線連接。為了避免管道出口位置發(fā)生回流現(xiàn)象對模擬計算造成干擾,出口位置采用相同結(jié)構(gòu)與DN15 的管線連接,管道物理模型見圖1。不同管徑管道的具體尺寸如表1所示。
表1 管道壁面厚度Table 1 Pipe wall thickness
為了研究熱分層現(xiàn)象在管道不同位置處的變化特性,在管道模型的外壁面上選擇了一系列具有代表性的位置設(shè)置溫度監(jiān)測點,在模擬計算過程中記錄監(jiān)測位置的溫度隨時間的變化,監(jiān)測位置選取見圖2。在直管段和彎管段共選擇了A~E 等5 個溫度監(jiān)測截面,每個截面在頂部和底部設(shè)置溫度監(jiān)測點,分別命名為1和2。在長度為1 000 mm的豎直管段上,等間距地設(shè)置了6 個溫度監(jiān)測點,分別命名為H1~H6。
圖2 特征溫度監(jiān)測點的位置選取Fig.2 Location selection of characteristic temperature monitoring points
SSTk-ω湍流模型屬于兩方程的雷諾時均應(yīng)力方程的一種,相比于k-ε湍流模型,該模型更適于具有逆壓梯度和分離流動的情況。該模型在動量方程中引入湍動能k和湍流耗散率ω,計算公式如下:
式中:ρ代表密度,kg·m-3;t代表時間,s;uj(j=1,2,3)代表速度分量,m·s·-1;xj(j=1,2,3)代表直角坐標分量,m;Gk、Gω代表由湍動能和湍流耗散率產(chǎn)生的源項;Yk、Yω代表發(fā)散項;Γk、Γω代表有效擴散項;Dω為對流擴散項;Sk、Sω為自定義的源項。
目前,對于熱浮力的計算,主要有Bossinesseq模型和全浮力模型兩種方法。兩種方法均在動量方程中添加浮力項,不同的是:Bossinesseq模型將浮力源項中的密度差等效處理為溫度差,動量方程中的密度項按參考溫度取值,該假設(shè)僅適用于密度變化遠小于密度值的情況,見式(3);全浮力模型采用式(4)計算流場中的浮力,適合溫差較大的流場計算[13-14]。
式中:SM,buoy代表浮力源項;ρref代表參考密度,kg·m-3;g代表重力加速度,m·s-2;T代表溫度,K;Tref代表參考溫度,K;β代表熱膨脹系數(shù),K-1。
壓水堆滯流分支管發(fā)生熱分層時,溫差較大,因此參考美國國家標準與技術(shù)研究院(National Institute of Standards and Technology,NIST)物性數(shù)據(jù)庫,將水的物性設(shè)置為與溫度相關(guān)的分段線性函數(shù)[15],見圖3。
圖3 水的物性參數(shù) (a) 水的密度與比熱容,(b) 水的導(dǎo)熱系數(shù)和黏度Fig.3 Physical properties of water (a) Density and specific heat capacity of water,(b) Thermal conductivity and viscosity of water
模擬計算中,假定管道壁面溫度及管內(nèi)流場溫度與空氣溫度相同,均為295.15 K。滯流分支管模型的壁面散熱情況分別為有保溫措施的散熱和無保溫措施的散熱,參考相關(guān)研究和工程設(shè)計數(shù)據(jù)[16-18],確定保溫管段和非保溫管段的壁面散熱系數(shù)。詳細的邊界條件設(shè)置如表2所示。
表2 邊界條件Table 2 Boundary conditions
本文采用FLUENT 2022 對滯流分支管熱分層現(xiàn)象進行數(shù)值模擬,分析了管道壁面的溫度變化特性和管道內(nèi)的流場分布特性。
計算模型的網(wǎng)格數(shù)量同時影響模擬計算的準確性和計算效率,因此在開展正式的研究前有必要進行網(wǎng)格無關(guān)性驗證。對滯流分支管采用不同的網(wǎng)格劃分方式得到了17萬、42萬、63萬、91萬網(wǎng)格4種方案,分別采用4 種網(wǎng)格方案在同一工況下進行瞬態(tài)計算。選取B截面的溫度監(jiān)測點B1和B2的溫度數(shù)據(jù)用于對照,如圖4所示。由4種網(wǎng)格方案計算得到的B1 監(jiān)測點溫度上升曲線近乎重合,但17 萬網(wǎng)格計算的B2 監(jiān)測點溫度數(shù)據(jù)與其余3 種網(wǎng)格方案存在較大區(qū)別。由此,63萬網(wǎng)格的方案兼具計算速度與計算準確性兩方面優(yōu)勢,最終確定網(wǎng)格數(shù)量為63萬,平均網(wǎng)格質(zhì)量為0.919 7。整體網(wǎng)格結(jié)構(gòu)與局部網(wǎng)格劃分情況見圖5。
圖4 網(wǎng)格無關(guān)性驗證Fig.4 Mesh independence verification
對于管道熱分層的數(shù)值模擬工作,廣泛采用的湍流模型有兩方程的Realizeablek-ε(RKE)湍流模型和SSTk-ω(SST KW)湍流模型,這兩種湍流模型具有收斂性好、內(nèi)存需求低的優(yōu)點。此外,本文也考慮了四方程的Transition SST 湍流模型(Transition Shear Stress Transport model,TSST),其常用于模擬湍流轉(zhuǎn)捩過程,在滯流分支管管徑存在變化導(dǎo)致流場雷諾數(shù)跨度較大時,該模型能夠更好地捕捉流場特性。為了比較三種湍流模型在管道熱分層數(shù)值模擬中的適用性和準確性,本文分別采用RKE 模型、SST KW 模型和TSST 模型在同一工況下開展瞬態(tài)計算,選擇B 截面的溫度監(jiān)測點B1 和B2 用于數(shù)據(jù)對照,計算結(jié)果如圖6所示。根據(jù)研究結(jié)果,在發(fā)生熱分層現(xiàn)象的管道中,豎直方向上存在的渦流現(xiàn)象會影響冷熱流體的換熱過程,導(dǎo)致截面溫差的增大。三種湍流模型中,采用RKE湍流模型計算得到的B1和B2監(jiān)測點的溫度差值相對較小,即截面溫差相對較小。表明該模型不能很好地模擬管內(nèi)流場在豎直方向上存在的渦流現(xiàn)象。因此,SST KW 模型和TSST 模型更適用于滯流分支管熱分層現(xiàn)象的模擬計算。而兩方程的SST KW湍流模型在計算結(jié)果方面與四方程的TSST湍流模型較為接近,同時,其計算效率更高。所以,本文選用SST KW 湍流模型進行模擬計算。
圖6 不同湍流模型的瞬態(tài)計算結(jié)果Fig.6 Transient calculation results using different turbulent models
基于上述滯流分支管模型開展模擬計算,分析管道壁面的溫度變化特性以及管道內(nèi)流體的流動特性。
2.1.1 管道截面的溫差變化
根據(jù)當前研究,熱分層現(xiàn)象容易出現(xiàn)在水平管段。滯流分支管模型中,水平管段包括DN100、DN200和DN350三種管徑,且不同管段分為有保溫措施和無保溫措施兩種情況,因此,水平管段內(nèi)的熱分層現(xiàn)象在不同位置處表現(xiàn)出差異。本文選取A、B、D 和E 4 個測溫截面,研究水平管段的熱分層現(xiàn)象。在每個測溫截面上,分別在管道頂部和底部設(shè)置溫度監(jiān)測點,以兩個監(jiān)測點的溫度差值代表截面溫差。采用管道截面溫差作為熱分層強度的特征參數(shù),模擬得到A、B、D 和E 4 個截面上管道頂部和底部的溫差隨時間的變化過程,如圖7所示。
圖7 水平管段的各測溫截面最大截面溫差隨時間的變化Fig.7 Variation of the maximum cross-sectional temperature difference for each temperature measurement section in a horizontal pipe section over time
水平管段的截面溫差總體上呈現(xiàn)先增大后減小的趨勢。然而,對于采取了保溫措施的管段,其截面溫差最終趨近于0,截面A和截面B所在管段的熱分層現(xiàn)象逐漸消失。相比之下,未采取保溫措施的管段,由于管道內(nèi)流體溫度的持續(xù)升高,管壁與空氣間的換熱作用逐漸增強。這種熱量散失作用導(dǎo)致截面上的溫差持續(xù)存在,截面D 和截面E 在瞬態(tài)計算達到穩(wěn)定時,所屬管段仍存在穩(wěn)定的熱分層現(xiàn)象,這會對管道結(jié)構(gòu)施加持續(xù)的應(yīng)力。
在瞬態(tài)過程中,不同管段處熱分層作用的嚴重程度存在差異。對于采取保溫措施的A截面和B截面而言,其最大截面溫差明顯低于D 截面和E 截面的對應(yīng)值。同時,需要注意的是,D截面和E截面所在管段的熱分層現(xiàn)象是持續(xù)存在的。
此外,在保溫條件相同的情況下,管徑也會影響截面最大溫差。具體而言,截面B 的最大溫差高于截面A,截面D的最大溫差高于截面E。這是因為管徑會影響截面上熱量傳遞的距離,管徑越大時,傳熱距離越遠,從而導(dǎo)致截面溫差越大。
2.1.2 水平管段的熱分層流動
根據(jù)研究[19],當熱流體涌入充滿冷水的管道中,分層流的壓力分布導(dǎo)致兩層流體的反向流動。在滯流分支管模型中,水平管段存在管徑變化,而管道結(jié)構(gòu)也會對流動產(chǎn)生影響。
當時間t=1 200 s時,截面B的溫差接近最大值,此時管道內(nèi)出現(xiàn)了明顯的熱分層現(xiàn)象。圖8展示了在該時刻,大小頭位置的速度場和溫度場分布。由于大小頭管段流動截面收縮造成的阻力作用,部分冷流體的流動方向發(fā)生改變,在管道底部產(chǎn)生與熱流體同反的流動。大小頭的管段結(jié)構(gòu)對熱分層流動產(chǎn)生了影響,在直徑為DN350 的直管段內(nèi),流場存在水平層內(nèi)的回流現(xiàn)象,將延緩管截面熱分層溫度梯度減緩的進程。
圖8 t=1 200 s時,大小頭管段的溫度場和流場分布Fig.8 Distributions of temperature field and flow field in the small and large end pipe section of transition pipe at 1 200 s
2.2.1 管道截面的溫差變化
彎管結(jié)構(gòu)對熱分層現(xiàn)象具有一定的影響。當熱流體通過彎頭時,其流動方向發(fā)生改變,從而影響管道截面上冷熱流體的換熱過程。以截面C所在的彎管段與截面B 所在的水平管段為例,通過對比彎管段與水平管段的截面溫差隨時間的變化,可以看出二者存在差異。B、C 截面所在管段管徑均為DN355,僅存在位置上的差異。截面溫差總體上的趨勢均表現(xiàn)為先增大后減小,最終趨于0。但是在熱分層的發(fā)展過程中,B截面始終具有更高的溫差,且上升趨勢更加明顯,如圖9所示。
圖9 B、C截面溫差隨時間的變化Fig.9 Variation of temperature difference between sections B and C over time
2.2.2 彎管段的熱分層流動
當時間t=2 000 s時,C截面的溫差接近最大值,此時彎管段的流場如圖10 所示。當熱流體流動至彎管段時,由于流動方向的改變,沿管道軸線方向的速度分量減小,而沿管道截面徑向的速度分量增加,促進了同一截面上冷熱流體的混合過程。因此,彎管段有效地減小了截面溫差。而在另一方面,由于熱流體沿管道軸線方向的速度分量減小,且浮力對熱流體的流動產(chǎn)生阻礙作用。因此,在豎直方向上冷熱流體的換熱過程受到影響,在圖10中可以看到豎直方向上有明顯的熱分層現(xiàn)象。
圖10 t=2 000 s時,彎管段內(nèi)流場Fig.10 Flow field of curved pipe section at 2 000 s
圖11 t=3 800 s時,流場的溫度分布Fig.11 Temperature distribution of the flow field at 3 800 s
當熱流體流動至豎直管段時,浮升力和慣性力的作用方向平行,管道截面上冷熱流體的換熱過程受浮力的影響程度降低,截面上的熱分層現(xiàn)象減弱。而在管道的軸線方向上,浮力對向下流動的熱流體產(chǎn)生了阻礙作用,冷熱流體溫度均一化的進程被減緩,導(dǎo)致沿管道軸線方向上出現(xiàn)了熱分層現(xiàn)象,如圖11所示。此時,管道內(nèi)的流動被稱為活塞流[20]。
在模擬計算中,本研究在豎直管道的同一側(cè)外壁面上等間距設(shè)置了一系列的溫度監(jiān)測點,用于觀察豎直管道上的溫度變化,測點的溫度數(shù)據(jù)如圖12所示。圖中H代表豎直管段的總長度,Hi代表任一溫度測點到H1 測點間的距離。豎直方向上的溫度梯度總體上表現(xiàn)為先增大后減小的趨勢。當t=3 800 s 時,管道內(nèi)的流場在豎直方向上的溫度跨度達到最大,為26 K。由于豎直管段未布置保溫層,管道壁面與空氣間存在較強的對流換熱,當瞬態(tài)模擬計算達到穩(wěn)定時,溫度監(jiān)測點H1 與H6 間存在恒定溫差。
圖12 豎直管道上的測點溫度Fig.12 Temperature data obtained at measurement points on the vertical pipe section
本文以壓水堆中與主回路連接的滯流分支管道為研究對象,利用數(shù)值模擬分析了由閥門泄漏引起的管道熱分層現(xiàn)象,總結(jié)了熱分層現(xiàn)象的溫度變化特性以及流動變化特性分析,得出以下結(jié)論:
1)滯流分支管中的熱分層現(xiàn)象受到管段位置、管徑大小以及管道保溫條件的影響。熱分層現(xiàn)象容易發(fā)生在水平管段,彎頭會削弱熱分層的影響,截面溫差的大小與管徑大小正相關(guān)。未采取保溫措施的管段,其截面溫差相對較大,且分層現(xiàn)象會持續(xù)存在。
2)滯流分支管模型存在管徑上的變化,變截面的管道結(jié)構(gòu)對管內(nèi)流場產(chǎn)生影響,在大小頭管段存在兩層方向相反的回流,導(dǎo)致分層作用的影響時間更長。
作者貢獻聲明馬競翔負責數(shù)值模擬研究執(zhí)行,數(shù)據(jù)分析,文章撰寫;董世昌負責文章審閱,內(nèi)容校核;龔圣捷負責提出研究思路、理論指導(dǎo),對文章作批評性審閱。