超臨界水通道內(nèi)壓降特性分析

徐 莉，臧金光,2，曾小康，閆 曉

(1.中國核動(dòng)力研究設(shè)計(jì)院 中核核反應(yīng)堆熱工水力技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610041；

2.清華大學(xué) 工程物理系，北京 100084)

在反應(yīng)堆中，堆芯釋放的熱量由系統(tǒng)回路中的冷卻劑帶出，堆芯的輸熱能力與冷卻劑的流動(dòng)特性密切相關(guān)，堆芯內(nèi)冷卻劑的流量分配對(duì)堆芯冷卻劑溫度場(chǎng)、燃料元件表面溫度分布均有重要影響。在熱工水力設(shè)計(jì)的過程中，應(yīng)盡量使冷卻劑的流量分布與堆芯功率分布相匹配，以最大限度地提高反應(yīng)堆的運(yùn)行功率，同時(shí)保證反應(yīng)堆的安全性。分析冷卻劑的壓降特性是堆芯熱工水力設(shè)計(jì)的重要方面[1]。

超臨界水冷堆是利用超臨界狀態(tài)下的水作為冷卻劑和慢化劑的第4代核能系統(tǒng)[2]。超臨界條件下由于存在劇烈的物性變化，壓降特性具有自身的特點(diǎn)。超臨界水在擬臨界點(diǎn)前后密度有較大落差，重力壓降需考慮密度的沿程積分效應(yīng)，同時(shí)密度的明顯變化使得加速壓降相比亞臨界條件下更為明顯。在等溫條件下，摩擦系數(shù)均為雷諾數(shù)的單值函數(shù)；而在加熱條件下，近壁面邊界層處的流體物性與主流體物性的差別較大，摩擦系數(shù)的計(jì)算需考慮近壁面處的物性修正。超臨界條件下的摩擦關(guān)系式多以亞臨界條件下的關(guān)系式為基礎(chǔ)加以擴(kuò)展。

本文對(duì)超臨界條件下的流動(dòng)特性進(jìn)行研究。

1 壓降類型

一般而言，通道內(nèi)的流動(dòng)壓降可分為加速壓降、重力壓降、摩擦壓降和局部壓降，本文著重分析加速壓降、重力壓降和摩擦壓降在超臨界條件下的變化規(guī)律。

在超臨界條件下，擬臨界區(qū)附近的密度變化較為劇烈，計(jì)算重力壓降時(shí)不能簡單認(rèn)為密度為常數(shù)，需考慮密度的沿程積分效應(yīng)。

(1)

其中：ρ為流體密度；g為重力加速度；Δpg為重力壓降；θ為流道與豎直方向的傾角；L為流動(dòng)距離；x為積分變量。

采用目前成熟的數(shù)值積分方法[3]對(duì)式(1)求解，如應(yīng)用較廣泛的Newton-Cotes求積公式：

(2)

其中：xi為求積節(jié)點(diǎn)；n為積分階數(shù)；Ai為求積系數(shù)，定義式如下：

(3)

其中：li(x)為拉格朗日插值基函數(shù)；a、b為積分上、下限。

當(dāng)n=1時(shí)，式(2)為梯形求積公式：

h=b-a

(4)

當(dāng)n=2時(shí)，式(2)為辛普森求積公式：

h=(b-a)/2

(5)

當(dāng)n=3時(shí)，式(2)為牛頓求積公式：

3f(x2)+f(x3))h=(b-a)/3

(6)

其中，xi=a+ih。

以辛普森求積公式計(jì)算重力壓降為例：

(7)

除以上數(shù)值積分方法外，也有研究者建議采用其他方法，如Ornatskiy等推薦采用焓值平均來考慮密度的積分效應(yīng)[4]：

(8)

為衡量這些積分表達(dá)式間的區(qū)別，分別用以上方法計(jì)算了壓力為25 MPa下不同焓值區(qū)間的平均密度，其中精確值采用精確分段法，將焓值的積分區(qū)間劃分50個(gè)節(jié)點(diǎn)，再線性平均后得到。圖1示出不同焓值區(qū)間下的平均密度，其中工況編號(hào)1～5分別代表焓值區(qū)間為1 850～2 000、1 850～2 200、1 850～2 400、1 850～2 600和2 300～3 000 kJ/kg。由圖1可見：當(dāng)焓值區(qū)間較小時(shí)，密度變化較小，采用不同方法計(jì)算的結(jié)果相差不大；當(dāng)焓值區(qū)間逐漸擴(kuò)大時(shí)，梯形求積法差異最大，最大相對(duì)誤差達(dá)11%，而辛普森求積法的最大相對(duì)誤差為0.9%，牛頓求積法的最大相對(duì)誤差為0.6%，Ornatskiy法的最大相對(duì)誤差為1.2%[4]。大部分情況下，采用Ornatskiy法計(jì)算密度的平均值是合理的，也易于實(shí)現(xiàn)。辛普森求積法在計(jì)算精度上較Ornatskiy法好些，但需多調(diào)用一次物性函數(shù)。

圖1 不同焓值區(qū)間的積分方法比較

亞臨界條件下，計(jì)算單相流的摩擦壓降Δpf普遍采用Darcy公式[1]：

(9)

其中：f為Darcy-Weisbach摩擦系數(shù)；Δz為單位通道長度；De為通道水力直徑；v為流體平均速度；G為質(zhì)量流速。

加速壓降Δpa的計(jì)算公式為：

(10)

式中：ρout、ρin分別為通道出口和入口的流體密度；Vout、Vin分別為相應(yīng)位置的比體積。

圖2 不同壓降梯度的沿程變化

圖2示出用CFX計(jì)算的光滑2×2棒束結(jié)構(gòu)各項(xiàng)壓降梯度的對(duì)比。由圖2可見，加速壓降梯度與摩擦壓降梯度沿程一直增大，而重力壓降梯度沿程一直減小。由3項(xiàng)壓降梯度的定義式可知，重力壓降梯度與平均密度呈正比，摩擦壓降梯度與流體密度呈反比，加速壓降梯度與流體比體積的落差呈正比。在整個(gè)加熱過程中，單位通道長度上的流體密度隨溫度的升高而下降，比體積上升，導(dǎo)致加速壓降梯度和摩擦壓降梯度增大和重力壓降梯度減小。在通道入口區(qū)，重力壓降梯度最大，其次是摩擦壓降梯度；之后，摩擦壓降梯度逐漸增大，而重力壓降梯度逐漸減小，使得前者慢慢超過后者。

圖3示出絕熱條件(冷態(tài))和加熱條件(熱態(tài))下各項(xiàng)壓降梯度的對(duì)比。絕熱條件下，加速壓降梯度為零。由圖3可見：摩擦壓降梯度和重力壓降梯度一直維持不變；加熱后，由于沿程物性的變化，加速壓降梯度產(chǎn)生，且沿程逐步增大；重力壓降梯度相比絕熱條件時(shí)減少，沿程伴隨密度的下降而下降；摩擦壓降梯度比絕熱條件時(shí)略有增加。加熱后總壓降梯度比絕熱條件時(shí)的略有上升。實(shí)際堆芯設(shè)計(jì)中，由于堆芯內(nèi)結(jié)構(gòu)復(fù)雜，棒束通道和通道阻塞物很多，加上在超臨界條件下冷卻劑溫度要跨過擬臨界區(qū)，因此可能會(huì)使堆芯內(nèi)的壓降變化很大。這是進(jìn)行超臨界條件流動(dòng)實(shí)驗(yàn)時(shí)必須要考慮的問題。

圖3 冷態(tài)與熱態(tài)時(shí)不同壓降梯度的沿程變化

2 亞臨界單相摩擦關(guān)系式

超臨界條件下的摩擦關(guān)系式多以亞臨界單相摩擦關(guān)系式為基準(zhǔn)，再輔以壁面物性與主流體物性之比作為修正，其優(yōu)點(diǎn)在于當(dāng)壁面溫度與主流體溫度之差較小時(shí)，摩擦關(guān)系式能自然過渡到正常等溫條件。因此，發(fā)展超臨界條件下的摩擦關(guān)系式需選取一合適的亞臨界等溫摩擦關(guān)系式為基礎(chǔ)。

亞臨界條件的等溫摩擦關(guān)系式中適用范圍較廣及預(yù)測(cè)精度較高的關(guān)系式之一是經(jīng)典的隱式PKN公式(式(11))，它是由Prandtl等根據(jù)大量實(shí)驗(yàn)測(cè)量結(jié)果得到的一半理論、半經(jīng)驗(yàn)關(guān)系式[5]，與大量實(shí)驗(yàn)數(shù)據(jù)具有較好的一致性。

(11)

隱式PKN公式計(jì)算較為復(fù)雜，需不斷迭代，用于發(fā)展超臨界條件下摩擦公式時(shí)會(huì)有些不便?；陔[式PKN公式形式，通過簡單數(shù)學(xué)推導(dǎo)，可得到PKN公式的顯式形式。顯式PKN公式在精度與適用范圍上與隱式PKN公式相當(dāng)，同時(shí)免去了復(fù)雜的迭代過程。顯式PKN公式是將Blausius公式代入隱式PKN公式，經(jīng)化簡整理，得：

f=1/(1.75lgRe-1.3)2

(12)

為便于分析比較，選取其他幾個(gè)常用公式作為比較。

McAdams公式：

f=0.184Re-0.2

(13)

Blausius公式：

f=0.316 4Re-0.25

(14)

Filonenko公式：

f=1/(1.82lgRe-1.64)2

(15)

為驗(yàn)證推導(dǎo)的顯式PKN公式與隱式PKN公式的接近程度，同時(shí)比較各摩擦系數(shù)公式的預(yù)測(cè)性能，分別計(jì)算了它們?cè)诓煌字Z數(shù)范圍內(nèi)的摩擦系數(shù)，結(jié)果示于圖4。由圖4可見：當(dāng)雷諾數(shù)較小時(shí)，McAdams公式的計(jì)算結(jié)果明顯小于隱式PKN公式的，隨著Re減小，誤差增大；當(dāng)

Re較高時(shí)，McAdams公式的計(jì)算結(jié)果逐漸與隱式PKN公式的接近，可見McAdams公式主要適用于高Re范圍；在所考察的Re范圍內(nèi)，Blausius與Filonenko公式的計(jì)算結(jié)果均與隱式PKN公式的接近，只是在Re較小時(shí)，F(xiàn)ilonenko公式的計(jì)算結(jié)果略大于隱式PKN公式的，Blausius公式的計(jì)算結(jié)果則略?。欢疚耐茖?dǎo)的顯式PKN公式的計(jì)算結(jié)果在整個(gè)Re范圍內(nèi)均與隱式PKN公式的最接近。這可從推導(dǎo)的過程中看出，顯式PKN公式是將Blausius公式作為初值代入隱式PKN公式獲得的，Blausius公式本身具有相對(duì)較好的預(yù)測(cè)性能，經(jīng)過隱式PKN公式的迭代修正后，會(huì)更接近于隱式PKN公式，這是它比其他公式預(yù)測(cè)性能更好的原因。圖5示出各摩擦關(guān)系式計(jì)算結(jié)果的區(qū)別。

3 超臨界摩擦關(guān)系式

針對(duì)超臨界條件下流動(dòng)實(shí)驗(yàn)的研究目前還不充分，缺少較普遍認(rèn)可的摩擦關(guān)系式，對(duì)其中的機(jī)理仍有待于進(jìn)一步研究。

以下是一些研究者基于實(shí)驗(yàn)數(shù)據(jù)擬合出的超臨界條件下的摩擦關(guān)系式[4]。

Mikheev公式：

f=fiso(Prw/Prb)1/3

(16)

其中：fiso為等溫流動(dòng)摩擦系數(shù)；Prw、Prb分別為以壁面溫度和主流體溫度為特征溫度得到的流體普朗特?cái)?shù)。

Popov公式：

/ρb)0.74

(17)

圖4 各摩擦關(guān)系式計(jì)算結(jié)果比較

圖5 各摩擦關(guān)系式計(jì)算結(jié)果的差異

Kondrat’ev公式：

f=0.188Re-0.22

(18)

Kirillov公式：

f=fiso(μw/μb)0.4

(19)

其中，μw和μb分別為壁面和主流體的動(dòng)力黏度。

上述公式中的fiso均采用Filonenko公式計(jì)算。鑒于實(shí)驗(yàn)數(shù)據(jù)不充分，本文利用CFD方法比較了這幾組基于實(shí)驗(yàn)的摩擦關(guān)系式與數(shù)值結(jié)果的差異。選取Mokry等[6]開展的超臨界條件下豎直圓管的傳熱實(shí)驗(yàn)為參考，管長為4 m，

直徑為10 mm，質(zhì)量流速G為203～1 495 kg/(m2·s)，熱流密度q為335～884 kW/m2，實(shí)驗(yàn)工況列于表1。

表1 實(shí)驗(yàn)工況

以Fluent軟件為計(jì)算平臺(tái)，ICEM為幾何建模和網(wǎng)格劃分工具，生成二維四面體結(jié)構(gòu)化網(wǎng)格。合理控制邊界層的網(wǎng)格參數(shù)，使第1層網(wǎng)格的y+<1，同時(shí)開展網(wǎng)格的無關(guān)性分析，確保計(jì)算結(jié)果不依賴于網(wǎng)格參數(shù)。入口采用定速度邊界條件，給定入口溫度；出口設(shè)置相對(duì)靜壓值為零，組件盒采用無滑移絕熱壁面邊界條件；元件棒表面依據(jù)計(jì)算工況不同設(shè)置定熱流密度邊界條件。計(jì)算時(shí)選用SST湍流模型。

圖6示出不同摩擦關(guān)系式與Fluent數(shù)值結(jié)果沿程預(yù)測(cè)性能的比較。由圖6可見，各摩擦關(guān)系式呈現(xiàn)較大的分散性。仔細(xì)分析這些關(guān)系式的變化趨勢(shì)可發(fā)現(xiàn)：大部分超臨界條件加熱工況下的摩擦系數(shù)均小于亞臨界條件等溫工況下的摩擦系數(shù)，F(xiàn)ilonenko關(guān)系式是亞臨界下的摩擦關(guān)系式，整體預(yù)測(cè)數(shù)值最高；超臨界加熱工況下，壁面溫度會(huì)遠(yuǎn)高于主流體溫度，近壁面處流體的黏性和密度均小于主流體的黏性和密度，導(dǎo)致摩擦系數(shù)變小。這與大部分摩擦關(guān)系式的修正方向是一致的。

圖6 Fluent與各關(guān)系式計(jì)算結(jié)果的比較

在擬臨界點(diǎn)附近，動(dòng)力黏度迅速下降，通道流體平均雷諾數(shù)上升，摩擦系數(shù)下降；同時(shí)，當(dāng)壁面溫度超過擬臨界點(diǎn)、而主流體溫度仍低于擬臨界點(diǎn)時(shí)，近壁面處的物性效應(yīng)明顯，加速了摩擦系數(shù)下降的趨勢(shì)。這是一些關(guān)系式計(jì)算的摩擦系數(shù)在擬臨界點(diǎn)處迅速下降的原因。在擬臨界點(diǎn)之后，流體溫度遠(yuǎn)離擬臨界點(diǎn)區(qū)域，壁面流體物性與主流體物性相差較小，物性效應(yīng)減弱，此時(shí)摩擦系數(shù)回升。擬臨界點(diǎn)處摩擦系數(shù)的極小值現(xiàn)象是超臨界條件下的獨(dú)有規(guī)律[7]。

在所比較的幾組超臨界摩擦關(guān)系式中，Kirillov公式與Fluent的計(jì)算結(jié)果較為一致，在分析的4組工況中，隨著質(zhì)量流速和熱流密度的變化，二者的吻合程度均較好，包括擬臨界點(diǎn)處的摩擦系數(shù)谷值，只是在入口區(qū)Kirillov公式的計(jì)算結(jié)果略低于Fluent的計(jì)算結(jié)果，這可能與CFD數(shù)值工具在模擬時(shí)考慮了入口效應(yīng)有關(guān)。圖7示出Kirillov公式與Fluent計(jì)算結(jié)果的比較，二者的相對(duì)偏差在±10%以內(nèi)。

圖7 Fluent與Kirillov公式計(jì)算結(jié)果比較

4 結(jié)論

本文分析了超臨界條件下劇烈的物性變化對(duì)通道內(nèi)壓降特性的影響，得到以下結(jié)論。

超臨界條件下，擬臨界區(qū)前后的密度變化劇烈，使得通道內(nèi)的加速壓降不可忽略，重力壓降需考慮密度的沿程積分效應(yīng)。本文比較了不同方法的計(jì)算精度，發(fā)現(xiàn)梯形求積公式誤差較大，Ornatskiy公式采用焓值平均可有效改善精度，辛普森求積公式與牛頓求積公式的精度較好。

基于隱式PKN公式形式和Blausisu公式得到了顯式PKN公式，計(jì)算精度與隱式PKN公式接近，而在具體求解時(shí)又免去隱式PKN公式迭代求解的繁瑣，適用于工程上的計(jì)算分析。

借助于Fluent的數(shù)值分析結(jié)果，比較了超臨界條件下不同摩擦關(guān)系式的異同，發(fā)現(xiàn)物性效應(yīng)修正會(huì)使擬臨界點(diǎn)附近出現(xiàn)局部極值點(diǎn)，且Kirillov公式與Fluent的計(jì)算結(jié)果較為一致，相對(duì)偏差在±10%以內(nèi)。

參考文獻(xiàn)：

[1] 于平安，朱瑞安，喻真烷，等. 核反應(yīng)堆熱工水力學(xué)[M]. 上海：上海交通大學(xué)出版社，2002.

[2] OKA Y, KOSHIZUKA S, ISHIWATARI Y, et al. Super light water reactors and super fast reactors[M]. New York: Springer, 2010.

[3] 鄧建中，劉之行. 計(jì)算方法[M]. 西安：西安交通大學(xué)出版社，2007.

[4] PIORO I L, DUFFEY R B. Heat transfer and hydraulic resistance at supercritical pressures in power-engineering applications[M]. New York: ASME Press, 2007.

[5] KAKA? S, SHAH R K, AUNG W. Handbook of single-phase convective heat transfer[M]. New York: Wiley-Interscience, 1987.

[6] MOKRY S, PIORO I L, KIRILLOV P, et al. Supercritical-water heat transfer in a vertical bare tube[J]. Nuclear Engineering and Design, 2010, 240: 568-576.

[7] LEI X L, LI H X, YU S Q, et al. Study on the minimum drag coefficient phenomenon of supercritical pressure water in the pseudocritical region[C]∥ICONE20-POWER2012. California, USA: [s. n.], 2012.