劉宇生，唐濟林，譚思超，高璞珍，張春明

（1.環(huán)境保護部核與輻射安全中心，北京 100082；2.哈爾濱工程大學(xué)核安全與仿真技術(shù)國防重點學(xué)科實驗室，哈爾濱 150001）

流量波動條件下圓管內(nèi)流動換熱特性數(shù)值研究

劉宇生1，唐濟林2，譚思超2，高璞珍2，張春明1

（1.環(huán)境保護部核與輻射安全中心，北京 100082；2.哈爾濱工程大學(xué)核安全與仿真技術(shù)國防重點學(xué)科實驗室，哈爾濱 150001）

對流量波動條件下核動力裝置系統(tǒng)圓管內(nèi)流體的流動換熱特性進行了數(shù)值研究，重點研究了波動周期、波動振幅兩個因素對流動換熱特性的影響。結(jié)果表明：流量波動時，通道內(nèi)的摩阻系數(shù)和壁面換熱系數(shù)均隨時間周期性波動，且波動周期與流量波動周期相同；波動周期增加，摩阻系數(shù)波動幅值減小，壁面換熱系數(shù)波動振幅變化不明顯；相對振幅增加，摩阻系數(shù)波動幅值增大，壁面換熱系數(shù)波動振幅增大。

流量波動；流動特性；傳熱特性；數(shù)值模擬；圓管

在核動力裝置運行過程中，流體流量波動的現(xiàn)象經(jīng)常出現(xiàn)。如非能動余熱排出系統(tǒng)在建立自然循環(huán)的過程中，二回路流體“冷凝過程比較劇烈，導(dǎo)致比較大的流量波動”［1］；在海洋條件或地震情況下，由于核動力裝置的傾斜、起伏和搖擺等運動，在附加慣性力的作用下，反應(yīng)堆中的冷卻劑流量會發(fā)生周期性波動，從而對堆芯傳熱產(chǎn)生影響［2-4］；對于采用自然循環(huán)非能動技術(shù)的核電廠而言，在事故工況下，由于“與強迫循環(huán)相比，目前自然循環(huán)實現(xiàn)的功率及流量都較小，容易出現(xiàn)流動不穩(wěn)定現(xiàn)象”［5］，因此也會出現(xiàn)流量波動現(xiàn)象。

目前，對上述瞬態(tài)過程的安全分析和可靠性評價，主要還是采用流體穩(wěn)定流動狀態(tài)下的流動換熱關(guān)系式，未考慮流量波動的影響［6-8］?！叭欢鴮Ψ悄軇影踩到y(tǒng)來說，其可靠性評價并非像能動設(shè)備失效那樣，需考慮物理過程失效”［9］，因此，有必要對流量波動條件下的流動換熱特性進行深入研究，為相關(guān)問題的審評提供參考。

1 控制方程與邊界條件

1.1 控制方程與湍流模型

封閉系統(tǒng)中，流體運動遵從質(zhì)量守恒、動量守恒和能量守恒。

（1）質(zhì)量守恒方程

單相流動換熱中，流體溫度升高對其密度的影響并不明顯，可假設(shè)流體密度為常數(shù)。其質(zhì)量守恒連續(xù)性方程為［10］：

式中，x、y、z分別為坐標軸；u、v、w分別為沿x、y、z軸的速度分量；ρ為流體的密度；t為時間。

（2）動量守恒方程

微元體中流體的動量對時間的變化率等于外界作用在該微元體上的各種力之和。按照這一定律，可得雷諾方程形式的動量守恒方程為［10］：

式中，μ為流體的動力粘度；p是壓力；其余參數(shù)意義同式（1）。

（3）能量守恒方程

微元體中能量的增加率等于進入微元體的凈熱流量加上體力與面力對微元體所做的功，通過引入導(dǎo)熱的Fourier定律，可得用流體的比焓h和溫度T表示的能量守恒方程為［11］：

式中，λ是流體的導(dǎo)熱系數(shù)；Sh為流體的內(nèi)熱源；divU為速度的散度；Φ是由于粘性作用機械能轉(zhuǎn)換為熱能的部分，稱為耗散函數(shù)；其余參數(shù)意義同式（1）、（2）。

（4）湍流模型

本文對流動和換熱的模擬主要采用標準k-ε模型，計算過程中，需要求解關(guān)于湍動能k和湍動能耗散率ε的方程［11］。

湍動能k的求解方程為：

湍動能耗散率ε的求解方程的最終形式為：

式（4）、（5）中，c1、c2為經(jīng)驗系數(shù)；ηt為湍流粘性系數(shù)，與湍動能k和耗散率ε有關(guān)；其余參數(shù)意義同以上各式。

求解時，方程（1）～（5）共同組成控制方程，依據(jù)邊界條件，沿計算節(jié)點求解計算。標準k-ε湍流模型是針對高雷諾數(shù)Re的湍流計算模型，適用于主流區(qū)域。對于近壁區(qū)，由于流動發(fā)展不充分，流動可能處于層流狀態(tài)，網(wǎng)格劃分時需加密近壁面處的網(wǎng)格，并在計算時采用壁面函數(shù)法加以修正。

1.2 邊界初始條件與求解控制

本研究不考慮系統(tǒng)管路布置的影響，僅分析水平圓管內(nèi)流體流動波動現(xiàn)象對流動和換熱的影響。由于管道水平且直徑很小，可忽略流體本身的重力作用。

（1）入口條件

假設(shè)管道入口處流體流速沿橫截面均勻分布，且入口流量隨時間呈正弦變化［12，13］，當管道截面一定時，流體入口流速滿足：

式中，u(t)為流體入口流速的瞬時值；Au和Bu分別為流體入口流速的振幅和周期平均值；T為流量波動的周期。

（2）出口條件

管道出口處的絕對壓強為1 atm。

（3）壁面條件

壁面采用無滑移條件，粗糙度為0.02 m。對紊流區(qū)的模擬按照水力光滑管處理。加熱方式為均勻熱流密度，熱流密度值在5 ～ 20 kW之間，均保證模擬通道內(nèi)的流體不發(fā)生沸騰。

（4）初始條件

t= 0時刻計算區(qū)域內(nèi)的流場和溫度場為速度u ( t )，等于u|t=0的定常流動情況下的流場和溫度場。即以具有相同邊界條件的穩(wěn)定流動工況的計算結(jié)果文件作為波動瞬態(tài)問題的初值條件。

（5）求解控制

對流項的離散格式采用高精度格式；瞬態(tài)項的離散格式采用二階歐拉向后差分格式；收斂準則選取的殘差類型為RMS，殘差控制在10-5以內(nèi)對于連續(xù)性方程和動量方程，殘差控制在10-6以內(nèi)。

瞬態(tài)問題中，根據(jù)波動周期的不同，時間步長選取為0.05～0.5 s不等，均保證小于或遠小于物理時間步長（物理時間步長定義為幾何體沿流動方向的長度與主流速度比值的一半）。每個時間步計算的最大循環(huán)次數(shù)為20次，實際求解過程中通過選取合理的時間步長使每步的迭代次數(shù)控制在5次以內(nèi)。

2 幾何模型與驗證

圖1 計算模型圖Fig.1 Geometry model for simulation

本計算所用的模型為直徑12 mm，壁厚2 mm，長1.7 m的圓管通道。為提高計算精度，流體區(qū)域采用六面體網(wǎng)格，并對近壁面處進行加密處理，如圖1所示。

圖2 計算參數(shù)與網(wǎng)格節(jié)點數(shù)的關(guān)系Fig.2 Relationships between simulated parameters and total number of mesh nodes

為確保模擬的準確性，對所建的網(wǎng)格模型進行了分析計算，如圖2所示。由圖2可知，隨著網(wǎng)格節(jié)點數(shù)的增加，壓降和換熱系數(shù)兩個參數(shù)會隨之變化，當網(wǎng)格節(jié)點數(shù)超過60萬后，計算結(jié)果受網(wǎng)格數(shù)量的影響就會很小。因此，計算中使用了節(jié)點數(shù)為60萬的網(wǎng)格模型。

此外，本研究還利用建立的網(wǎng)格模型，對圓管內(nèi)的強制對流換熱進行了數(shù)值模擬，并將所得努謝爾特數(shù)Nu與采用迪圖斯—貝爾特關(guān)聯(lián)式［14］計算的結(jié)果進行了對比，如圖3所示。由圖3可知，計算值與經(jīng)驗公式的結(jié)果吻合較好，這表明所建網(wǎng)格模型合理可靠。

圖3 數(shù)值計算的結(jié)果與經(jīng)驗公式的對比Fig.3 The comparison of heat transfer coeffcients between simulation and correlation

文獻［15］給出了層流不加熱條件下，流量隨時間正弦波動時壓降ΔP的表達式：

式中，M、N、D為與管道尺寸、流體物性相關(guān)的常數(shù)；φ為相位。

利用所建網(wǎng)格模型，對恒溫層流量波動現(xiàn)象進行了模擬，將計算結(jié)果與由公式（7）得出的理論值進行了對比，如圖4所示。數(shù)值模擬與理論分析的結(jié)果均表明，流量周期性波動時，管道壓降也會周期性波動，且壓降波動與流量波動的周期相同，但二者之間存在相位差。

由圖4可知，壓降的模擬值與理論值存在偏差，但模擬值在整體的變化趨勢方面與理論值符合較好；另外，流量波動與壓降波動的相位差體現(xiàn)了理論分析中流量變化滯后于壓差變化的特點。

圖4 模擬壓降與理論壓降圖Fig.4 Simulated and theoretical pressure drop

3 流量波動對流動換熱的影響

流體流量周期性波動時，水平管中流動的總壓降ΔP包括摩擦壓降ΔPf和加速壓降ΔPa兩部分［16］：

由于本研究中的流體為單相狀態(tài)，密度變化較小，與摩阻壓降相比，加速壓降的值很小，可以忽略。為保證計算結(jié)果具有可比性，周期平均速度設(shè)為定值。

3.1 波動周期的影響

在壁面熱流密度q為20 kW·m-2，速度振幅為0.8 m·s-1，周期平均速度為2 m·s-1的工況下，計算了不同波動周期T下摩阻系數(shù)λ隨時間的變化，如圖5所示。由圖5可知，摩阻系數(shù)波動的周期與流量波動的周期相同，也存在相位差；摩阻系數(shù)的波動幅度隨著周期的減小逐漸增大，并且當周期小于10 s后，摩阻系數(shù)的波動幅度顯著增加。同時可知，波動過程中，在流體加速階段與減速階段，摩阻系數(shù)的變化也不同，在加速階段，摩阻系數(shù)上升的速率較慢，在減速階段，摩阻系數(shù)下降較快。

圖5 不同波動周期下摩阻系數(shù)隨時間的變化Fig.5 The impact of different fuctuation periods on friction coeff cient

圖6 不同波動周期對表面換熱系數(shù)的影響Fig.6 The impact of different fuctuation periods on heat transfer coeff cient

圖6為不同波動周期T下，壁面平均換熱系數(shù)的變化規(guī)律。由圖6可知，隨著流量的周期性波動，壁面換熱系數(shù)也呈現(xiàn)出周期性的波動，且周期與流量波動相同；不同的波動周期對壁面換熱系數(shù)的波動幅度影響較小，不同周期下壁面換熱系數(shù)的值相差很小。

圖7 不同波動周期對通道內(nèi)流體溫度分布的影響Fig.7 The impact of different fuctuation periods on temperature distribution in tube

圖7為相對振幅和熱流密度相同，波動周期不同時，沿圓管中心軸線和距離入口0.85 m處截面上沿直徑的流體溫度的分布情況。由圖7可知從管道壁面到中心軸線，流體溫度逐漸降低，呈開口向上的拋物線型分布；沿管道的中心軸線方向，流體溫度逐漸升高，呈線性增加趨勢。

流量波動周期減小，沿管道直徑方向和沿中心軸線方向的流體溫度均有所提高，這表明波動周期變小后，換熱得到強化，因為周期變小，流體的波動頻率增加，利于流體微團間的熱量傳遞

3.2 相對振幅的影響

為分析波動振幅和波動時的平均流量對流動和換熱的影響，將速度波動振幅Bu與周期平均速度Bu的比值定義為相對振幅Ar，用以表征流量波動的劇烈程度。

在周期為10 s，熱流密度為20 kW·m-2的工況下，計算相對振幅對流動的影響。

圖8為不同相對振幅下，摩阻系數(shù)隨時間的變化。由圖8可知，隨著相對振幅的增加，摩阻系數(shù)波動的幅值增加，尤其是相對振幅超過0.5以后摩阻系數(shù)顯著增大。這主要是因為相對振幅增加會使近壁區(qū)的速度梯度增大，造成摩擦損失增加。

圖9為不同相對振幅下，壁面換熱系數(shù)隨時間的變化情況。由圖9可知，隨著相對振幅的增大壁面換熱系數(shù)波動的幅度增加，這是因為相對振幅增大，流速波動的幅值增加，使流體微團之間的擾動增強，壁面邊界層周期性減薄，熱阻減小換熱得到了強化。

此外，局部的速度變化對壁面換熱系數(shù)的影響最為顯著，因此圖9中表面換熱系數(shù)波動的峰值和谷值總是與流量波動的峰值和谷值同時出現(xiàn)，這與文獻［17］的研究結(jié)論一致。

圖8 不同的相對振幅對摩阻系數(shù)的影響Fig.8 The impact of different relative amplitudes on friction coefficient

圖9 不同波動相對振幅對壁面換熱系數(shù)的影響Fig.9 The impact of different relative amplitudes on heat transfer coeff cient

圖10為波動相對振幅不同時，沿圓管中心軸線和距離入口0.85 m處截面上沿直徑的流體溫度的分布情況。由圖可知，沿管道直徑方向，從管道壁面到中心軸線，流體溫度逐漸降低；沿管道中心軸線方向，流體溫度逐漸升高。這與不同波動周期下的溫度分布規(guī)律一致。

波動相對振幅增大，管道軸向各點和徑向各點的流體溫度均升高，這表明波動振幅增加，壁面換熱系數(shù)增大，換熱效果得到了強化。

圖10 不同相對振幅對通道內(nèi)流體溫度分布的影響Fig.10 The impact of different relative amplitudes on temperature distribution in tube

3.3 流量波動導(dǎo)致的相位差

圖11和圖12分別為波動周期為10 s時，流體流速、流體平均溫度、表面換熱系數(shù)、摩阻系數(shù)等參數(shù)隨時間的變化情況。

綜合圖11和圖12可知：隨著流量的波動，摩阻系數(shù)、表面換熱系數(shù)和流體平均溫度等參數(shù)均隨時間發(fā)生周期性的波動，而且這些參量的波動與流量的波動之間存在相位差。

摩阻系數(shù)的波動滯后于速度的波動，二者之間的相位差接近80°。這與Uchida對層流脈動流的理論分析結(jié)論相符合［18］。

對于確定的工況，換熱情況主要取決于換熱面流體的局部流速，所以表面換熱系數(shù)與流速之間未發(fā)現(xiàn)相位差，流速變化，換熱系數(shù)隨之改變。流體的平均溫度不僅與換熱系數(shù)相關(guān)，還與換熱壁面的溫度、流速等因素相關(guān)，在流量波動時，流體的平均溫度變化明顯滯后于流速變化，與流速、表面換熱系數(shù)的變化趨勢相反，存在約90°的相位差。

圖11 摩阻系數(shù)、表面換熱系數(shù)和速度隨時間的變化情況Fig.11 The change curves of friction coeffcient， heat transfer coeffcient and velocity with time

圖12 流體速度、平均溫度和換熱系數(shù)隨時間的變化情況Fig.12 The change curves of velocity， fuid average temperature and heat transfer coeffcient with time

4 結(jié)論

當流體流量呈正弦周期性波動時，管道壓降、摩阻系數(shù)、流體溫度、壁面換熱系數(shù)等參數(shù)都會隨之出現(xiàn)周期性的波動；其中管道壓降、摩阻系數(shù)、流體溫度與流量之間存在相位差，換熱系數(shù)與流量之間不存在相位差。

波動周期越小，波動相對振幅越大，摩阻系數(shù)的波動振幅越大；周期平均流量相同，波動相對振幅越大，壁面換熱系數(shù)的波動振幅越大，波動周期對壁面換熱系數(shù)的影響不明顯。

對流量波動條件下溫度場的分析表明，周期平均流量相同，波動周期減小，波動相對振幅增加，沿徑向和軸向通道內(nèi)流體的溫度升高，換熱效果增強。

［1］俞冀陽，俞爾俊，賈寶山. 核反應(yīng)堆非能動余熱排出系統(tǒng)分析［J］. 核安全，2004（ 4 ）：35-39.

［2］高璞珍，劉順隆，王兆祥. 縱搖和橫搖對自然循環(huán)的影響［J］. 核動力工程，1999，20（3）：228-231.

［3］曹夏昕，閻昌琪，孫立成，等.搖擺狀態(tài)下豎直管內(nèi)單相水摩擦壓降特性分析［J］. 哈爾濱工程大學(xué)學(xué)報，2006，27（6）：835-838.

［4］譚思超，張紅巖，龐鳳閣，等.搖擺運動下單相自然循環(huán)流動特點［J］. 核動力工程，2005，26（6）：559-563.

［5］周濤，李精精，琚忠云，等. 非能動自然循環(huán)技術(shù)的發(fā)展與研究［J］. 核安全，2013，12（3）：32-36.

［6］蘇云，許以全，曹學(xué)武，等.SCDAP/RELAP5程序結(jié)構(gòu)及嚴重事故有關(guān)的模型概述［J］.核動力工程，2003，24（6S2）：51-55.

［7］武俊梅，蘇光輝.流動不穩(wěn)定性對CHF的影響［J］.核動力工程，2007，28（3）：19-21.

［8］陳玉宙，楊春生，毛玉龍.穩(wěn)態(tài)和瞬態(tài)流動條件下低壓欠熱沸騰水的臨界熱流密度實驗研究［J］.中國原子能科學(xué)研究院年報，2005：178.

［9］黃昌藩，匡波. 非能動安全系統(tǒng)可靠性評估方法初步研究［J］. 核安全，2012（1）：35 - 41.

［10］陶文銓.數(shù)值傳熱學(xué)［M］.西安：西安交通大學(xué)出版社，2001.

［11］王福軍. 計算流體動力學(xué)分析－CFD軟件原理與應(yīng)用［M］.北京：清華大學(xué)出版社，2004.

［12］劉廷浩.流量波動條件下阻力特性實驗研究［D］.哈爾濱：哈爾濱工程大學(xué)，2008.

［13］劉廣山.流量周期性變化流場數(shù)值分析［D］.哈爾濱：哈爾濱工程大學(xué)，2010.

［14］楊世銘，陶文銓. 傳熱學(xué)3版［M］. 北京：高等教育出版社，1998.

［15］Liu Tinghao，Gao Puzhen，Tan Sichao. Study on phase difference between fux and pressure drop of pipe fow under flow fluctuation condition. Session B：System Simulation Technologies，Proceedings of of International Symposium on Symbiotic Nuclear Power Systems（ISSNP） for 21st Century［C］. Harbin：Harbin Engineering University，2008.

［16］閻昌琪.氣液兩相流［M］.哈爾濱：哈爾濱工程大學(xué)出版社，1995.

［17］王暢.流量波動條件下的流動與換熱特性研究［D］.哈爾濱：哈爾濱工程大學(xué)，2010.

［18］Uchida S. The pulsating viscous flow superposed on the steady lam inar motion of incompressible fuid in a circular pipe［J］. Zeitschrift für angewandte M athematik und Physik ZAMP，1956，7（5）：403-422.

Numerical Simulation of Fluid Flow and Heat Transfer Characteristics under Flow Fluctuation Condition in Circular Tube

LIU Yusheng1，TANG Jilin2，TAN Sichao2，GAO Puzhen2，ZHANG Chunm ing1
（1.Nuclear and Radiation Safety Center，MEP，Beijing 100082，China；2.National Defense Key Subject Laboratory for Nuclear Safety and Simulation Technology，Harbin Engineering University，Harbin 150001，China）

Flow and heat transfer characteristics in a circular tube of nuclear power system under flow fluctuation condition are investigated numerically. The influences of the flow fluctuation period and relative amplitude on flow and heat transfer characteristics are mainly analyzed. The results show that the pressure drop，friction coeff cient，wall heat transfer coeffcient vary in sin law with time as fow fuctuates. What is more，their periods are the same with fow fuctuation period. In addition，the amplitude of friction coeffcient decreases as fow fuctuation period becomes larger，while the amplitude of wall heat transfer coeffcient increases. From a relative amplitude point of view，the parameters including friction coeffcient，wall heat transfer coeffcient increase when relative amplitude becomes larger.

fow fuctuation；fow characteristics；heat transfer characteristics；numerical simulation；circular tube

TL334

A

1672-5360（2014）01-0039-06

2013-11-24

2014-02-24

劉宇生（1986—），男，河北唐山人，助理工程師，現(xiàn)主要從事反應(yīng)堆熱工水力方面的研究工作