壓水堆滯流分支管熱分層現(xiàn)象的數(shù)值模擬

馬競(jìng)翔 董世昌 龔圣捷

（上海交通大學(xué) 上海 200240）

壓水堆（Pressurized Water Reactor，PWR）核電站一回路以水為冷卻介質(zhì)［1］。部分與冷卻劑系統(tǒng)主管道相連接的分支管內(nèi)流體處于滯流狀態(tài)，分支管內(nèi)的流體溫度低而主管道內(nèi)的流體溫度高［2］。由于湍流滲透或閥門泄漏等原因，主管道內(nèi)的高溫流體進(jìn)入分支管內(nèi)，從而形成熱分層現(xiàn)象［3］。

在水平或傾斜管道中，熱流體緩慢流入冷流體管段時(shí)，由于溫度不同造成的密度差，熱流體浮于冷流體的上層，在重力方向上形成溫度梯度，即熱分層現(xiàn)象［4］。熱分層會(huì)導(dǎo)致管道壁面應(yīng)力分布不均勻，并引起管道熱疲勞及管道結(jié)構(gòu)失效［5］。

自20世紀(jì)80年代起，壓水堆核電站管路的熱分層現(xiàn)象被發(fā)現(xiàn)后就受到了廣泛關(guān)注。美國(guó)核管理委員會(huì)（U.S.Nuclear Regulatory Commission，NRC）發(fā)布公告要求所有在役或在建核電廠必須對(duì)管路熱分層現(xiàn)象進(jìn)行分析論證，以確保管道結(jié)構(gòu)的完整性［6-7］。根據(jù)研究，壓水堆核電站冷卻系統(tǒng)的多處結(jié)構(gòu)會(huì)受到熱分層現(xiàn)象的影響，例如：蒸汽發(fā)生器的給水管線、主冷卻系統(tǒng)熱管段、高壓安注管線、余熱排出管線等部位［8-9］。

賀斌［10］對(duì)核電站中蒸汽發(fā)生器上水管路熱分層的流動(dòng)與傳熱機(jī)理進(jìn)行了數(shù)值模擬研究。發(fā)現(xiàn)二次流和浮升力的作用會(huì)使得冷熱流體交界面不斷變化，并驗(yàn)證了冷熱流體的相對(duì)流速和相對(duì)溫差會(huì)影響流場(chǎng)的溫度變化趨勢(shì)。曹瓊等［11］對(duì)AP1000 第四級(jí)卸壓A管線ADS-4A的熱分層開展穩(wěn)態(tài)數(shù)值模擬研究。發(fā)現(xiàn)保溫層的設(shè)置使得管道沿程的溫度梯度減小；ADS-4A管線的45°彎頭能夠有效降低由湍流滲透作用進(jìn)入管線中的熱流體的流速。Jo等［4］采用商用計(jì)算流體動(dòng)力學(xué)（Computational Fluid Dynamics，CFD）軟件對(duì)穩(wěn)壓器波動(dòng)管進(jìn)行了波入和波出兩種工況下的數(shù)值模擬工作，提出了管壁的影響，認(rèn)為計(jì)算模型包含管壁并采用共軛傳熱的方式，能更準(zhǔn)確地預(yù)測(cè)管道內(nèi)的熱分層現(xiàn)象。Cai 等［12］基于Richardson數(shù)（Ri）相似，建立了1/3比例的穩(wěn)壓器波動(dòng)管模型，開展了波入和波出兩種工況條件下的模擬實(shí)驗(yàn)。通過實(shí)驗(yàn)數(shù)據(jù)與模擬結(jié)果對(duì)比發(fā)現(xiàn)：SSTk-ω湍流模型的精度要優(yōu)于標(biāo)準(zhǔn)k-ε湍流模型（Standardk-εModel，SKE）和雷諾應(yīng)力模型（Reynolds Stress Model，RSM）。

目前，熱分層現(xiàn)象的研究主要圍繞等截面管道模型開展，而對(duì)于變截面管道的瞬態(tài)流動(dòng)特性和溫度變化特性的研究較少。本研究建立了存在管徑變化的滯流分支管模型，開展了數(shù)值模擬工作，分析其在閥門泄漏情況下發(fā)生熱分層現(xiàn)象時(shí)的溫度變化特性和管道內(nèi)的流場(chǎng)特性，為后續(xù)的模擬實(shí)驗(yàn)和管道應(yīng)力分析提供理論依據(jù)。

1 滯流分支管物理建模和數(shù)值模型

1.1 滯流分支管模型尺寸

滯流分支管采用管徑分別為DN100、DN200、DN350的直管段、大小頭以及彎管段組成，管道只有部分位置布置有保溫層。為了模擬滯流分支管因入口處閥門泄漏而引發(fā)熱分層現(xiàn)象，管道入口位置采用封板結(jié)構(gòu)與DN15的管線連接。為了避免管道出口位置發(fā)生回流現(xiàn)象對(duì)模擬計(jì)算造成干擾，出口位置采用相同結(jié)構(gòu)與DN15 的管線連接，管道物理模型見圖1。不同管徑管道的具體尺寸如表1所示。

表1 管道壁面厚度Table 1 Pipe wall thickness

1.2 滯流分支管特征溫度監(jiān)測(cè)點(diǎn)的選取

為了研究熱分層現(xiàn)象在管道不同位置處的變化特性，在管道模型的外壁面上選擇了一系列具有代表性的位置設(shè)置溫度監(jiān)測(cè)點(diǎn)，在模擬計(jì)算過程中記錄監(jiān)測(cè)位置的溫度隨時(shí)間的變化，監(jiān)測(cè)位置選取見圖2。在直管段和彎管段共選擇了A～E 等5 個(gè)溫度監(jiān)測(cè)截面，每個(gè)截面在頂部和底部設(shè)置溫度監(jiān)測(cè)點(diǎn)，分別命名為1和2。在長(zhǎng)度為1 000 mm的豎直管段上，等間距地設(shè)置了6 個(gè)溫度監(jiān)測(cè)點(diǎn)，分別命名為H1～H6。

圖2 特征溫度監(jiān)測(cè)點(diǎn)的位置選取Fig.2 Location selection of characteristic temperature monitoring points

1.3 SST k-ω湍流模型與浮力模型

SSTk-ω湍流模型屬于兩方程的雷諾時(shí)均應(yīng)力方程的一種，相比于k-ε湍流模型，該模型更適于具有逆壓梯度和分離流動(dòng)的情況。該模型在動(dòng)量方程中引入湍動(dòng)能k和湍流耗散率ω，計(jì)算公式如下：

式中：ρ代表密度，kg·m-3；t代表時(shí)間，s；uj（j=1，2，3）代表速度分量，m·s·-1；xj（j=1，2，3）代表直角坐標(biāo)分量，m；Gk、Gω代表由湍動(dòng)能和湍流耗散率產(chǎn)生的源項(xiàng)；Yk、Yω代表發(fā)散項(xiàng)；Γk、Γω代表有效擴(kuò)散項(xiàng)；Dω為對(duì)流擴(kuò)散項(xiàng)；Sk、Sω為自定義的源項(xiàng)。

目前，對(duì)于熱浮力的計(jì)算，主要有Bossinesseq模型和全浮力模型兩種方法。兩種方法均在動(dòng)量方程中添加浮力項(xiàng)，不同的是：Bossinesseq模型將浮力源項(xiàng)中的密度差等效處理為溫度差，動(dòng)量方程中的密度項(xiàng)按參考溫度取值，該假設(shè)僅適用于密度變化遠(yuǎn)小于密度值的情況，見式（3）；全浮力模型采用式（4）計(jì)算流場(chǎng)中的浮力，適合溫差較大的流場(chǎng)計(jì)算［13-14］。

式中：SM，buoy代表浮力源項(xiàng)；ρref代表參考密度，kg·m-3；g代表重力加速度，m·s-2；T代表溫度，K；Tref代表參考溫度，K；β代表熱膨脹系數(shù)，K-1。

1.4 水的物性參數(shù)

壓水堆滯流分支管發(fā)生熱分層時(shí)，溫差較大，因此參考美國(guó)國(guó)家標(biāo)準(zhǔn)與技術(shù)研究院（National Institute of Standards and Technology，NIST）物性數(shù)據(jù)庫(kù)，將水的物性設(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

1.5 邊界條件

模擬計(jì)算中，假定管道壁面溫度及管內(nèi)流場(chǎng)溫度與空氣溫度相同，均為295.15 K。滯流分支管模型的壁面散熱情況分別為有保溫措施的散熱和無(wú)保溫措施的散熱，參考相關(guān)研究和工程設(shè)計(jì)數(shù)據(jù)［16-18］，確定保溫管段和非保溫管段的壁面散熱系數(shù)。詳細(xì)的邊界條件設(shè)置如表2所示。

表2 邊界條件Table 2 Boundary conditions

1.6 物理建模與網(wǎng)格劃分

本文采用FLUENT 2022 對(duì)滯流分支管熱分層現(xiàn)象進(jìn)行數(shù)值模擬，分析了管道壁面的溫度變化特性和管道內(nèi)的流場(chǎng)分布特性。

計(jì)算模型的網(wǎng)格數(shù)量同時(shí)影響模擬計(jì)算的準(zhǔn)確性和計(jì)算效率，因此在開展正式的研究前有必要進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證。對(duì)滯流分支管采用不同的網(wǎng)格劃分方式得到了17萬(wàn)、42萬(wàn)、63萬(wàn)、91萬(wàn)網(wǎng)格4種方案，分別采用4 種網(wǎng)格方案在同一工況下進(jìn)行瞬態(tài)計(jì)算。選取B截面的溫度監(jiān)測(cè)點(diǎn)B1和B2的溫度數(shù)據(jù)用于對(duì)照，如圖4所示。由4種網(wǎng)格方案計(jì)算得到的B1 監(jiān)測(cè)點(diǎn)溫度上升曲線近乎重合，但17 萬(wàn)網(wǎng)格計(jì)算的B2 監(jiān)測(cè)點(diǎn)溫度數(shù)據(jù)與其余3 種網(wǎng)格方案存在較大區(qū)別。由此，63萬(wàn)網(wǎng)格的方案兼具計(jì)算速度與計(jì)算準(zhǔn)確性兩方面優(yōu)勢(shì)，最終確定網(wǎng)格數(shù)量為63萬(wàn)，平均網(wǎng)格質(zhì)量為0.919 7。整體網(wǎng)格結(jié)構(gòu)與局部網(wǎng)格劃分情況見圖5。

圖4 網(wǎng)格無(wú)關(guān)性驗(yàn)證Fig.4 Mesh independence verification

1.7 湍流模型的選擇

對(duì)于管道熱分層的數(shù)值模擬工作，廣泛采用的湍流模型有兩方程的Realizeablek-ε（RKE）湍流模型和SSTk-ω（SST KW）湍流模型，這兩種湍流模型具有收斂性好、內(nèi)存需求低的優(yōu)點(diǎn)。此外，本文也考慮了四方程的Transition SST 湍流模型（Transition Shear Stress Transport model，TSST），其常用于模擬湍流轉(zhuǎn)捩過程，在滯流分支管管徑存在變化導(dǎo)致流場(chǎng)雷諾數(shù)跨度較大時(shí)，該模型能夠更好地捕捉流場(chǎng)特性。為了比較三種湍流模型在管道熱分層數(shù)值模擬中的適用性和準(zhǔn)確性，本文分別采用RKE 模型、SST KW 模型和TSST 模型在同一工況下開展瞬態(tài)計(jì)算，選擇B 截面的溫度監(jiān)測(cè)點(diǎn)B1 和B2 用于數(shù)據(jù)對(duì)照，計(jì)算結(jié)果如圖6所示。根據(jù)研究結(jié)果，在發(fā)生熱分層現(xiàn)象的管道中，豎直方向上存在的渦流現(xiàn)象會(huì)影響冷熱流體的換熱過程，導(dǎo)致截面溫差的增大。三種湍流模型中，采用RKE湍流模型計(jì)算得到的B1和B2監(jiān)測(cè)點(diǎn)的溫度差值相對(duì)較小，即截面溫差相對(duì)較小。表明該模型不能很好地模擬管內(nèi)流場(chǎng)在豎直方向上存在的渦流現(xiàn)象。因此，SST KW 模型和TSST 模型更適用于滯流分支管熱分層現(xiàn)象的模擬計(jì)算。而兩方程的SST KW湍流模型在計(jì)算結(jié)果方面與四方程的TSST湍流模型較為接近，同時(shí)，其計(jì)算效率更高。所以，本文選用SST KW 湍流模型進(jìn)行模擬計(jì)算。

圖6 不同湍流模型的瞬態(tài)計(jì)算結(jié)果Fig.6 Transient calculation results using different turbulent models

2 結(jié)果與討論

基于上述滯流分支管模型開展模擬計(jì)算，分析管道壁面的溫度變化特性以及管道內(nèi)流體的流動(dòng)特性。

2.1 水平管段的熱分層現(xiàn)象

2.1.1 管道截面的溫差變化

根據(jù)當(dāng)前研究，熱分層現(xiàn)象容易出現(xiàn)在水平管段。滯流分支管模型中，水平管段包括DN100、DN200和DN350三種管徑，且不同管段分為有保溫措施和無(wú)保溫措施兩種情況，因此，水平管段內(nèi)的熱分層現(xiàn)象在不同位置處表現(xiàn)出差異。本文選取A、B、D 和E 4 個(gè)測(cè)溫截面，研究水平管段的熱分層現(xiàn)象。在每個(gè)測(cè)溫截面上，分別在管道頂部和底部設(shè)置溫度監(jiān)測(cè)點(diǎn)，以兩個(gè)監(jiān)測(cè)點(diǎn)的溫度差值代表截面溫差。采用管道截面溫差作為熱分層強(qiáng)度的特征參數(shù)，模擬得到A、B、D 和E 4 個(gè)截面上管道頂部和底部的溫差隨時(shí)間的變化過程，如圖7所示。

圖7 水平管段的各測(cè)溫截面最大截面溫差隨時(shí)間的變化Fig.7 Variation of the maximum cross-sectional temperature difference for each temperature measurement section in a horizontal pipe section over time

水平管段的截面溫差總體上呈現(xiàn)先增大后減小的趨勢(shì)。然而，對(duì)于采取了保溫措施的管段，其截面溫差最終趨近于0，截面A和截面B所在管段的熱分層現(xiàn)象逐漸消失。相比之下，未采取保溫措施的管段，由于管道內(nèi)流體溫度的持續(xù)升高，管壁與空氣間的換熱作用逐漸增強(qiáng)。這種熱量散失作用導(dǎo)致截面上的溫差持續(xù)存在，截面D 和截面E 在瞬態(tài)計(jì)算達(dá)到穩(wěn)定時(shí)，所屬管段仍存在穩(wěn)定的熱分層現(xiàn)象，這會(huì)對(duì)管道結(jié)構(gòu)施加持續(xù)的應(yīng)力。

在瞬態(tài)過程中，不同管段處熱分層作用的嚴(yán)重程度存在差異。對(duì)于采取保溫措施的A截面和B截面而言，其最大截面溫差明顯低于D 截面和E 截面的對(duì)應(yīng)值。同時(shí)，需要注意的是，D截面和E截面所在管段的熱分層現(xiàn)象是持續(xù)存在的。

此外，在保溫條件相同的情況下，管徑也會(huì)影響截面最大溫差。具體而言，截面B 的最大溫差高于截面A，截面D的最大溫差高于截面E。這是因?yàn)楣軓綍?huì)影響截面上熱量傳遞的距離，管徑越大時(shí)，傳熱距離越遠(yuǎn)，從而導(dǎo)致截面溫差越大。

2.1.2 水平管段的熱分層流動(dòng)

根據(jù)研究［19］，當(dāng)熱流體涌入充滿冷水的管道中，分層流的壓力分布導(dǎo)致兩層流體的反向流動(dòng)。在滯流分支管模型中，水平管段存在管徑變化，而管道結(jié)構(gòu)也會(huì)對(duì)流動(dòng)產(chǎn)生影響。

當(dāng)時(shí)間t=1 200 s時(shí)，截面B的溫差接近最大值，此時(shí)管道內(nèi)出現(xiàn)了明顯的熱分層現(xiàn)象。圖8展示了在該時(shí)刻，大小頭位置的速度場(chǎng)和溫度場(chǎng)分布。由于大小頭管段流動(dòng)截面收縮造成的阻力作用，部分冷流體的流動(dòng)方向發(fā)生改變，在管道底部產(chǎn)生與熱流體同反的流動(dòng)。大小頭的管段結(jié)構(gòu)對(duì)熱分層流動(dòng)產(chǎn)生了影響，在直徑為DN350 的直管段內(nèi)，流場(chǎng)存在水平層內(nèi)的回流現(xiàn)象，將延緩管截面熱分層溫度梯度減緩的進(jìn)程。

圖8 t=1 200 s時(shí)，大小頭管段的溫度場(chǎng)和流場(chǎng)分布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 彎管段的熱分層現(xiàn)象

2.2.1 管道截面的溫差變化

彎管結(jié)構(gòu)對(duì)熱分層現(xiàn)象具有一定的影響。當(dāng)熱流體通過彎頭時(shí)，其流動(dòng)方向發(fā)生改變，從而影響管道截面上冷熱流體的換熱過程。以截面C所在的彎管段與截面B 所在的水平管段為例，通過對(duì)比彎管段與水平管段的截面溫差隨時(shí)間的變化，可以看出二者存在差異。B、C 截面所在管段管徑均為DN355，僅存在位置上的差異。截面溫差總體上的趨勢(shì)均表現(xiàn)為先增大后減小，最終趨于0。但是在熱分層的發(fā)展過程中，B截面始終具有更高的溫差，且上升趨勢(shì)更加明顯，如圖9所示。

圖9 B、C截面溫差隨時(shí)間的變化Fig.9 Variation of temperature difference between sections B and C over time

2.2.2 彎管段的熱分層流動(dòng)

當(dāng)時(shí)間t=2 000 s時(shí)，C截面的溫差接近最大值，此時(shí)彎管段的流場(chǎng)如圖10 所示。當(dāng)熱流體流動(dòng)至彎管段時(shí)，由于流動(dòng)方向的改變，沿管道軸線方向的速度分量減小，而沿管道截面徑向的速度分量增加，促進(jìn)了同一截面上冷熱流體的混合過程。因此，彎管段有效地減小了截面溫差。而在另一方面，由于熱流體沿管道軸線方向的速度分量減小，且浮力對(duì)熱流體的流動(dòng)產(chǎn)生阻礙作用。因此，在豎直方向上冷熱流體的換熱過程受到影響，在圖10中可以看到豎直方向上有明顯的熱分層現(xiàn)象。

圖10 t=2 000 s時(shí)，彎管段內(nèi)流場(chǎng)Fig.10 Flow field of curved pipe section at 2 000 s

圖11 t=3 800 s時(shí)，流場(chǎng)的溫度分布Fig.11 Temperature distribution of the flow field at 3 800 s

2.3 豎直管道內(nèi)的熱分層現(xiàn)象

當(dāng)熱流體流動(dòng)至豎直管段時(shí)，浮升力和慣性力的作用方向平行，管道截面上冷熱流體的換熱過程受浮力的影響程度降低，截面上的熱分層現(xiàn)象減弱。而在管道的軸線方向上，浮力對(duì)向下流動(dòng)的熱流體產(chǎn)生了阻礙作用，冷熱流體溫度均一化的進(jìn)程被減緩，導(dǎo)致沿管道軸線方向上出現(xiàn)了熱分層現(xiàn)象，如圖11所示。此時(shí)，管道內(nèi)的流動(dòng)被稱為活塞流［20］。

在模擬計(jì)算中，本研究在豎直管道的同一側(cè)外壁面上等間距設(shè)置了一系列的溫度監(jiān)測(cè)點(diǎn)，用于觀察豎直管道上的溫度變化，測(cè)點(diǎn)的溫度數(shù)據(jù)如圖12所示。圖中H代表豎直管段的總長(zhǎng)度，Hi代表任一溫度測(cè)點(diǎn)到H1 測(cè)點(diǎn)間的距離。豎直方向上的溫度梯度總體上表現(xiàn)為先增大后減小的趨勢(shì)。當(dāng)t=3 800 s 時(shí)，管道內(nèi)的流場(chǎng)在豎直方向上的溫度跨度達(dá)到最大，為26 K。由于豎直管段未布置保溫層，管道壁面與空氣間存在較強(qiáng)的對(duì)流換熱，當(dāng)瞬態(tài)模擬計(jì)算達(dá)到穩(wěn)定時(shí)，溫度監(jiān)測(cè)點(diǎn)H1 與H6 間存在恒定溫差。

圖12 豎直管道上的測(cè)點(diǎn)溫度Fig.12 Temperature data obtained at measurement points on the vertical pipe section

3 結(jié)語(yǔ)

本文以壓水堆中與主回路連接的滯流分支管道為研究對(duì)象，利用數(shù)值模擬分析了由閥門泄漏引起的管道熱分層現(xiàn)象，總結(jié)了熱分層現(xiàn)象的溫度變化特性以及流動(dòng)變化特性分析，得出以下結(jié)論：

1）滯流分支管中的熱分層現(xiàn)象受到管段位置、管徑大小以及管道保溫條件的影響。熱分層現(xiàn)象容易發(fā)生在水平管段，彎頭會(huì)削弱熱分層的影響，截面溫差的大小與管徑大小正相關(guān)。未采取保溫措施的管段，其截面溫差相對(duì)較大，且分層現(xiàn)象會(huì)持續(xù)存在。

2）滯流分支管模型存在管徑上的變化，變截面的管道結(jié)構(gòu)對(duì)管內(nèi)流場(chǎng)產(chǎn)生影響，在大小頭管段存在兩層方向相反的回流，導(dǎo)致分層作用的影響時(shí)間更長(zhǎng)。

作者貢獻(xiàn)聲明馬競(jìng)翔負(fù)責(zé)數(shù)值模擬研究執(zhí)行，數(shù)據(jù)分析，文章撰寫；董世昌負(fù)責(zé)文章審閱，內(nèi)容校核；龔圣捷負(fù)責(zé)提出研究思路、理論指導(dǎo)，對(duì)文章作批評(píng)性審閱。