魯曉東，周鈴嵐，巫英偉，蘇光輝，秋穗正，張 虹

(1.西安交通大學(xué) 動(dòng)力工程多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049；

2.西安交通大學(xué) 核科學(xué)與技術(shù)系，陜西 西安 710049；

3.中國(guó)核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610041)

兩相流動(dòng)不穩(wěn)定性現(xiàn)象廣泛存在于工業(yè)系統(tǒng)和設(shè)備中，例如蒸汽發(fā)生器、沸水堆和制冷裝置等。它會(huì)引起系統(tǒng)流量和壓力振蕩，致使系統(tǒng)裝置發(fā)生機(jī)械振蕩影響系統(tǒng)控制，危害系統(tǒng)安全。因此在系統(tǒng)設(shè)計(jì)和運(yùn)行時(shí)，需避免流動(dòng)不穩(wěn)定性現(xiàn)象的發(fā)生。在過(guò)去的幾十年中，對(duì)兩相流動(dòng)不穩(wěn)定性進(jìn)行了大量的實(shí)驗(yàn)研究和數(shù)值計(jì)算。近年來(lái)，許多學(xué)者[1-4]將目光投向了并聯(lián)通道兩相流動(dòng)不穩(wěn)定性的研究。在并聯(lián)系統(tǒng)中，各通道具有相同邊界條件，相互間存在影響。當(dāng)某個(gè)通道發(fā)生擾動(dòng)時(shí)，另外的通道也會(huì)受到影響，這使并聯(lián)通道兩相流動(dòng)不穩(wěn)定性的研究特別復(fù)雜和困難。在已有的研究中，廣泛采用軸向均勻加熱的假設(shè)來(lái)研究并聯(lián)通道系統(tǒng)流動(dòng)不穩(wěn)定性。然而，在真實(shí)的核反應(yīng)堆系統(tǒng)中，堆芯軸向功率分布是非均勻的。因此，為更加真實(shí)地研究反應(yīng)堆系統(tǒng)中的兩相流動(dòng)不穩(wěn)定性現(xiàn)象，本文選取余弦功率加熱計(jì)算系統(tǒng)的穩(wěn)定性邊界和三維不穩(wěn)定性空間，以研究軸向非均勻加熱對(duì)并聯(lián)通道兩相流動(dòng)不穩(wěn)定性的影響。

1 物理模型和數(shù)值解法

1.1 系統(tǒng)結(jié)構(gòu)

在本研究中，將反應(yīng)堆堆芯系統(tǒng)簡(jiǎn)化為由上下聯(lián)箱和中間兩個(gè)通道組成的系統(tǒng)，如圖1所示。為更好地分析問(wèn)題，本文作如下假設(shè)：1) 兩相區(qū)采用均相流模型；2) 入口流體過(guò)冷；3) 兩相流體處于熱力學(xué)平衡狀態(tài)；4) 忽略過(guò)冷沸騰。

1.2 物理模型

對(duì)于均相流模型，單相區(qū)和兩相區(qū)的控制方程均可寫(xiě)成如下形式。

質(zhì)量守恒方程：

(1)

動(dòng)量守恒方程：

(2)

圖1 并聯(lián)通道系統(tǒng)簡(jiǎn)圖

能量守恒方程：

(3)

狀態(tài)方程：

ρ=ρ(p,h)

(4)

式中：p為壓力；h為流體焓；ρ為流體密度；u為流體流速；A為流體流動(dòng)截面積；De為當(dāng)量直徑；f為摩擦系數(shù)；ql為線功率；z為坐標(biāo)軸方向；ki為節(jié)流系數(shù)；δ為Dirac函數(shù)。

式(4)兩邊分別對(duì)時(shí)間t求導(dǎo)，可得：

··

(5)

將式(5)分別代入式(1)和式(3)，得：

··(ρu)=0

(6)

(7)

在計(jì)算壓降時(shí)，需選擇適當(dāng)?shù)膄和兩相摩擦倍增因子Φ2。表1列出單相區(qū)摩擦系數(shù)f1φ，兩相區(qū)摩擦系數(shù)f2φ為：

Φ2

(8)

Φ2采用MaAdams模型，即：

(9)

表1 摩擦系數(shù)關(guān)系式

1.3 數(shù)值方法

空間上采用交錯(cuò)網(wǎng)格技術(shù)，即將系統(tǒng)的p、h、ρ等物性參數(shù)儲(chǔ)存在控制體中心(如圖2中虛線位置)，將u存放于控制體邊界處(如圖2中實(shí)線位置)，從而可省略同位網(wǎng)格技術(shù)中的迭代步驟，提高計(jì)算速度。在時(shí)間上采用半隱式差分方式進(jìn)行離散，避免了完全隱式差分求解困難和耗時(shí)的缺點(diǎn)。當(dāng)界面速度大于零時(shí)，邊界處流體物性參數(shù)取來(lái)流方向上最后1個(gè)控制體中心的值[5]。

圖2 一維空間網(wǎng)格

根據(jù)圖2的空間控制體劃分方式，在第i個(gè)控制體處對(duì)式(6)、(7)進(jìn)行離散，可得：

·

(10)

(11)

(12)

上述3組差分離散方程求解時(shí)是封閉的。當(dāng)時(shí)間步長(zhǎng)滿(mǎn)足Δt<Δz/u時(shí)，上述方程采用的離散格式穩(wěn)定且精確。

(13)

其中：

1.4 分析方法

由于入口段和上升段對(duì)并聯(lián)通道兩相流動(dòng)不穩(wěn)定性均有影響[4]，為更好地研究軸向非均勻加熱對(duì)流動(dòng)不穩(wěn)定性的影響，本文將不考慮入口段和上升段。在計(jì)算不穩(wěn)定邊界時(shí)，采用小擾動(dòng)法進(jìn)行研究，即在并聯(lián)通道系統(tǒng)處于穩(wěn)定狀態(tài)時(shí)，對(duì)其中1根通道的加熱量在極短時(shí)間(1 s)內(nèi)給定一微小擾動(dòng)(增加或減小1%)，系統(tǒng)運(yùn)行一段時(shí)間后，觀察兩個(gè)通道的流量脈動(dòng)曲線和相空間軌跡，判斷系統(tǒng)所處的運(yùn)行狀態(tài)。

當(dāng)進(jìn)口溫度相同時(shí)，隨加熱功率的增加，流量脈動(dòng)依次經(jīng)歷阻尼振蕩、極限環(huán)振蕩和發(fā)散振蕩，如圖3所示，圖中Npch為相變數(shù)、Nsub為進(jìn)口過(guò)冷度數(shù)。選取極限環(huán)振蕩時(shí)的加熱功率為臨界功率，并以此獲得不穩(wěn)定性邊界。

圖3 入口流量振蕩曲線(a、b、c)和相空間軌跡(d、e、f)

2 模型驗(yàn)證

本文采用雙矩形通道實(shí)驗(yàn)結(jié)果[6]對(duì)模型進(jìn)行驗(yàn)證。本模型Npch的計(jì)算值和實(shí)驗(yàn)值的比較結(jié)果如圖4所示。從圖4可看出，本模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的相對(duì)誤差均在15%以?xún)?nèi)，說(shuō)明本文計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果符合良好。但由于本文采用均相流模型，忽略了相間滑移，因此計(jì)算結(jié)果偏保守。

圖4 Npch計(jì)算值與實(shí)驗(yàn)值的對(duì)比

3 結(jié)果分析

表2列出并聯(lián)通道的詳細(xì)參數(shù)。本文采用軸向余弦功率加熱模擬軸向非均勻功率加熱，將均勻功率加熱作為參考功率分布，如圖5所示。

表2 并聯(lián)通道詳細(xì)參數(shù)

圖5 不同的加熱功率方式

3.1 不穩(wěn)定性邊界

本文采用進(jìn)口過(guò)冷度數(shù)和相變數(shù)表征并聯(lián)通道系統(tǒng)的不穩(wěn)定邊界，定義如下：

(14)

式中：hf和hg為飽和液體和飽和蒸汽焓；hi為進(jìn)口液體焓；vfg為飽和蒸汽與飽和液體比體積之差；vf為飽和液體比體積；hfg為汽化潛熱；Q為加熱量；W為流量。

在以前的文獻(xiàn)中，只有少量基于實(shí)驗(yàn)數(shù)據(jù)研究軸向非均勻加熱對(duì)兩相系統(tǒng)穩(wěn)定性的影響，但存在兩種截然相反的結(jié)論。Djikman等[7]發(fā)現(xiàn)余弦功率加熱可增強(qiáng)系統(tǒng)的穩(wěn)定性；然而B(niǎo)ergles[8]則認(rèn)為余弦功率加熱會(huì)降低系統(tǒng)的穩(wěn)定性。采用本文所得的圖6，可很好地解釋這兩種結(jié)論。

圖6 余弦功率加熱對(duì)系統(tǒng)穩(wěn)定性的影響

圖6為7 MPa下余弦功率加熱對(duì)系統(tǒng)穩(wěn)定性的影響。從圖6可看出，在高過(guò)冷度區(qū)，余弦功率加熱的系統(tǒng)穩(wěn)定性邊界(MSB)相對(duì)于均勻功率加熱的MSB向左移動(dòng)，系統(tǒng)穩(wěn)定性降低。但在低過(guò)冷度區(qū)，余弦功率加熱的MSB向右移動(dòng)，余弦功率加熱增強(qiáng)了系統(tǒng)穩(wěn)定性。在不同的系統(tǒng)壓力下，可得到相似結(jié)論，如圖7所示。

產(chǎn)生上述現(xiàn)象的原因如下：對(duì)于兩相系統(tǒng)，單相段壓降對(duì)系統(tǒng)的穩(wěn)定性起著很重要的作用。單相段壓降越大，系統(tǒng)越穩(wěn)定。對(duì)于余弦功率加熱和均勻功率加熱，在高過(guò)冷度區(qū)，沸騰邊界均處于通道的上半?yún)^(qū)，但余弦功率加熱時(shí)沸騰邊界較均勻功率加熱時(shí)低，單相段壓降相對(duì)較低，因此系統(tǒng)穩(wěn)定性降低。在低過(guò)冷度區(qū)，對(duì)于上述兩種加熱方式，通道的沸騰邊界均發(fā)生在下半部，而余弦功率加熱的沸騰邊界較均勻功率加熱的高，單相段壓降相對(duì)較大，因此系統(tǒng)穩(wěn)定性增強(qiáng)。

圖7 系統(tǒng)的MSB

3.2 三維不穩(wěn)定性空間

從圖7可看出，二維的MSB并不能清楚地顯示壓力對(duì)系統(tǒng)穩(wěn)定性的影響。Guo等[4]在2008年提出的不穩(wěn)定性空間(不穩(wěn)定島)的概念解決了這個(gè)問(wèn)題。在本研究中，采用由Npch、Nsub和無(wú)量綱壓力(p+=p/10patmosphere)組成的三維空間來(lái)研究不同加熱方式下的系統(tǒng)穩(wěn)定性。

本文計(jì)算了余弦和均勻功率加熱在壓力為7、9和15 MPa下的并聯(lián)通道MSB，獲得了相應(yīng)加熱方式下的不穩(wěn)定性空間，如圖8所示。從圖8可看出，對(duì)于余弦和均勻功率加熱，隨系統(tǒng)壓力的增加，不穩(wěn)定性區(qū)域均逐漸變小，并聯(lián)通道系統(tǒng)的穩(wěn)定性增強(qiáng)。

圖8 均勻功率(a)和余弦功率(b)加熱的不穩(wěn)定性空間

3.3 進(jìn)口阻力系數(shù)的影響

對(duì)于均勻功率加熱的系統(tǒng)，已有很多學(xué)者研究了進(jìn)口阻力系數(shù)對(duì)系統(tǒng)穩(wěn)定性的影響。然而對(duì)于余弦功率加熱的系統(tǒng)，這方面的研究卻進(jìn)行得很少。因此本文研究了余弦功率加熱時(shí)進(jìn)口阻力系數(shù)kin對(duì)并聯(lián)通道系統(tǒng)穩(wěn)定性的影響，如圖9所示。從圖9可看出，系統(tǒng)穩(wěn)定性隨進(jìn)口阻力系數(shù)的增大而增強(qiáng)。因增加進(jìn)口阻力系數(shù)增大了單相段壓降，從而增強(qiáng)了系統(tǒng)穩(wěn)定性。此外，從圖9還可知，隨進(jìn)口阻力系數(shù)的增大，MSB的拐點(diǎn)逐漸向高過(guò)冷度移動(dòng)。

圖9 進(jìn)口阻力系數(shù)對(duì)系統(tǒng)穩(wěn)定性的影響

4 結(jié)論

1) 余弦功率加熱對(duì)系統(tǒng)穩(wěn)定性的影響在不同的過(guò)冷度區(qū)域表現(xiàn)為不同的形式：在高過(guò)冷度區(qū)，降低系統(tǒng)的穩(wěn)定性；在低過(guò)冷度區(qū)，增強(qiáng)系統(tǒng)的穩(wěn)定性。

2) 無(wú)論是余弦功率加熱還是均勻功率加熱，系統(tǒng)壓力均會(huì)增強(qiáng)并聯(lián)通道系統(tǒng)的穩(wěn)定性。

3) 當(dāng)系統(tǒng)處于余弦功率加熱時(shí)，進(jìn)口阻力系數(shù)能增強(qiáng)系統(tǒng)的穩(wěn)定性。同時(shí)，隨進(jìn)口阻力系數(shù)的增大，MSB的拐點(diǎn)逐漸向高過(guò)冷度區(qū)移動(dòng)。

參考文獻(xiàn)：

[1] MUNOZ-COBO J L, PODOWSKI M Z, CHIVA S. Parallel channel instabilities in boiling water reactor systems: Boundary conditions for out of phase oscillations[J]. Annals of Nuclear Energy, 2002, 29(16): 1 891-1 917.

[2] LEE J D, PAN C. Dynamics of mutiple parallel boiling channel systems with forced flows[J]. Nuclear Engineering and Design, 1999, 192(1): 31-44.

[3] ZHANG Y J, SU G H, YANG X B, et al. Theoretical research on two-phase flow instability in parallel channels[J]. Nuclear Engineering and Design, 2009, 239(7): 1 294-1 303.

[4] GUO Y, QIU S Z, SU G H, et al. Theoretical investigations on two-phase flow instability in parallel multichannel system[J]. Annals of Nuclear Energy, 2008, 35(4): 665-676.

[5] Office of Nuclear Regulatory Research. RELAP5/MOD3.3 code manual volume Ⅵ: Validation of numerical techniques in RELAP5/MOD3.3[M]. US: Nuclear Regulatory Commission, 2001: 53-62.

[6] LU D H, LI H B, CHEN B D, et al. Experimental study on the density wave oscillation in parallel rectangular channels with narrow gap[J]. Annals of Nuclear Energy, 2011, 38(10): 2 146-2 155.

[7] DJIKMAN F J M, SLUITER M L C. One-channel non-linear digital model for boiling water hydrodynamics[J]. Nuclear Engineering and Design, 1971, 16(3): 237-248.

[8] BERGLES A E. Review of instability in two phase systems[C]∥KAKAC S, MAYINGER F. Two Phase Flow and Heat Transfer. New York: NATO Advanced Study Institute, 1976.