基于一維漂移流模型的并聯(lián)矩形雙通道密度波流動(dòng)不穩(wěn)定性數(shù)值模擬

熊萬(wàn)玉，唐瑜*，陳炳德

(中國(guó)核動(dòng)力研究設(shè)計(jì)院，四川 成都610041)

摘要：基于一維漂移流模型構(gòu)建了并聯(lián)矩形雙通道密度波流動(dòng)不穩(wěn)定性數(shù)學(xué)模型。模型中采用Zuber推薦的經(jīng)驗(yàn)關(guān)系式計(jì)算兩相流體空泡份額，采用Chisholm關(guān)系式和中國(guó)核動(dòng)力研究設(shè)計(jì)院自擬關(guān)系式計(jì)算兩相流體的摩擦壓降。求解過(guò)程中將質(zhì)量方程、能量方程與動(dòng)量方程解耦，并在計(jì)算域內(nèi)沿流動(dòng)方向依次求解方程組。計(jì)算過(guò)程中，首先開(kāi)展穩(wěn)態(tài)計(jì)算，在穩(wěn)態(tài)解的基礎(chǔ)上，通過(guò)添加流量或功率擾動(dòng)，誘發(fā)流體周期性振蕩，通過(guò)辨識(shí)瞬態(tài)計(jì)算中得到的流量振蕩模式來(lái)獲得流動(dòng)不穩(wěn)定邊界。采用數(shù)值計(jì)算獲得的密度波脈動(dòng)圖像與實(shí)驗(yàn)中觀察到的密度波脈動(dòng)現(xiàn)象的特征基本一致。最后，針對(duì)16組典型實(shí)驗(yàn)工況開(kāi)展數(shù)值模擬，結(jié)果表明，大部分工況下計(jì)算不穩(wěn)定界限熱流密度與實(shí)驗(yàn)值的相對(duì)偏差小于±20%。

關(guān)鍵詞：并聯(lián)矩形雙通道；密度波不穩(wěn)定性；漂移流模型；數(shù)值模擬

中圖分類(lèi)號(hào)：TK124 文獻(xiàn)標(biāo)志碼：A

收稿日期：2014-08-01;修回日期：2014-11-05

作者簡(jiǎn)介：熊萬(wàn)玉(1973—)，女，湖北公安人，研究員，碩士，從事反應(yīng)堆熱工水力研究

doi：10.7538/yzk.2015.49.11.1989

*通信作者：唐瑜，E-mail: tangyu614@hotmail.com

Numerical Simulation on Density Wave Oscillation in Parallel Twin

Rectangular Channels Based on One-dimensional Drift Flux Model

XIONG Wan-yu, TANG Yu*, CHEN Bing-de

(NuclearPowerInstituteofChina,Chengdu610041,China)

Abstract:A mathematical model was proposed for the density wave oscillation in parallel twin rectangular channels based on one-dimensional drift flux model in this paper. In the model, an empirical correlation recommended by Zuber was applied to calculate the void fraction of two-phase flow. Chisholm correlation and the correlation developed by NPIC were used for the prediction of the frictional pressure drop of two-phase flow. During the process of solving, mass, momentum and energy balance equations were decoupled and solved along the flow direction successively in the whole computational domain. The steady state calculation was made at first, then the periodical oscillation of flow was induced by adding small flow rate or heating power disturbance to the steady state solution. The flow instability boundary can be obtained by checking the oscillation mode in transient calculation. The calculation results indicate that the flow oscillations obtained through simulation have the same characteristics as those observed in the experiment. At last, the simulation were made for sixteen typical cases chosen from the experiment, and the results indicate that the relative differences between the calculated flow instability boundary heat flux and the corresponding experimental data are no more than ±20% for most part of cases.

Key words：parallel twin rectangular channels; density wave oscillation; drift flux model; numerical simulation

兩相流動(dòng)不穩(wěn)定性是反應(yīng)堆熱工水力研究的重要課題之一，它涉及到反應(yīng)堆設(shè)計(jì)、運(yùn)行、維護(hù)等多個(gè)方面的內(nèi)容。根據(jù)流量變化方式不同，流動(dòng)不穩(wěn)定性可分為靜態(tài)流動(dòng)不穩(wěn)定性和動(dòng)態(tài)流動(dòng)不穩(wěn)定性。靜態(tài)流動(dòng)不穩(wěn)定性是指流量由一個(gè)值變化到另一個(gè)值，出現(xiàn)新的工況，且系統(tǒng)仍維持在新工況下工作。動(dòng)態(tài)流動(dòng)不穩(wěn)定性多是指有固定頻率的脈動(dòng)現(xiàn)象，系統(tǒng)進(jìn)口流量圍繞一平均值可能呈現(xiàn)收斂的、發(fā)散的或穩(wěn)定的脈動(dòng)變化[1]。密度波流動(dòng)不穩(wěn)定性是反應(yīng)堆中最為常見(jiàn)的動(dòng)態(tài)不穩(wěn)定性之一，它與流量、密度以及壓降三者之間的相互作用有關(guān)。由于密度波不穩(wěn)定性對(duì)反應(yīng)堆設(shè)計(jì)及安全運(yùn)行有重要影響，迄今為止，研究者已針對(duì)密度波脈動(dòng)問(wèn)題開(kāi)展了大量的實(shí)驗(yàn)和理論研究。就理論研究而言，主要可分為時(shí)域法和頻域法兩大類(lèi)，本文所開(kāi)展的流動(dòng)不穩(wěn)定性數(shù)值模擬屬于時(shí)域分析法。時(shí)域法一般從兩相流解析數(shù)學(xué)模型開(kāi)始，建立守恒方程組。然后把非線性方程組進(jìn)行有限差分離散，進(jìn)行數(shù)值計(jì)算。從預(yù)設(shè)初始條件開(kāi)始，計(jì)算流量、壓差、溫度等物理量隨時(shí)間的變化，來(lái)判斷預(yù)設(shè)的初始工況下是否會(huì)發(fā)生流動(dòng)不穩(wěn)定。采用時(shí)域法進(jìn)行流動(dòng)不穩(wěn)定性研究的學(xué)者很多，早在1963年，Meyer等[2]就采用歐拉方法建立了守恒方程，并利用有限差分方法在時(shí)域內(nèi)求解方程組。在他們的模型中，整個(gè)系統(tǒng)是一個(gè)兩端壓力固定的沸騰流道。1982年，Dogan等[3]運(yùn)用集總參數(shù)法對(duì)垂直沸騰兩相流動(dòng)單管內(nèi)的壓力降型脈動(dòng)進(jìn)行了數(shù)值求解。1990年，Lecoq等[4]利用時(shí)域法對(duì)蒸汽發(fā)生器內(nèi)密度波型脈動(dòng)進(jìn)行了數(shù)值計(jì)算，通過(guò)兩相駐留時(shí)間，分析了延遲效應(yīng)對(duì)密度波型脈動(dòng)的影響，另外還分析了極限環(huán)振蕩。黃軍等[5]以一維漂移流模型為基礎(chǔ)，對(duì)管間脈動(dòng)建立了計(jì)算模型，采用C語(yǔ)言和SIMPLER算法編制了程序，并調(diào)試通過(guò)。李虹波[6]以圓管平行雙通道管間脈動(dòng)理論模型為基礎(chǔ)，針對(duì)矩形通道進(jìn)行修正，表征出矩形雙通道內(nèi)管間脈動(dòng)的動(dòng)態(tài)行為。夏庚磊等[7]采用RELAP5程序?qū)Υ怪辈⒙?lián)管中氣液兩相流動(dòng)不穩(wěn)定性實(shí)驗(yàn)裝置進(jìn)行模擬，并與實(shí)驗(yàn)結(jié)果對(duì)比，發(fā)現(xiàn)計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果符合較好。

本文基于一維漂移流模型構(gòu)建并聯(lián)矩形雙通道密度波流動(dòng)不穩(wěn)定性數(shù)學(xué)模型，給出該模型的數(shù)值求解方法，選定16組計(jì)算工況開(kāi)展流動(dòng)不穩(wěn)定性數(shù)值模擬，獲得各工況下流動(dòng)不穩(wěn)定的界限熱流密度，繪制相應(yīng)的流動(dòng)不穩(wěn)定邊界，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比。

1數(shù)學(xué)物理模型

圖1　并聯(lián)雙通道幾何模型 Fig.1　Geometry model of parallel twin channels

1.1幾何模型

并聯(lián)雙通道的幾何結(jié)構(gòu)如圖1所示。圖1中，兩個(gè)結(jié)構(gòu)完全相同的矩形通道分別與上下聯(lián)箱相連接后構(gòu)成并聯(lián)通道，在聯(lián)箱內(nèi)部流體的壓力相等。流體從下聯(lián)箱進(jìn)入裝置后分別經(jīng)由兩個(gè)通道垂直向上流動(dòng)，在流動(dòng)過(guò)程中受到管壁的加熱后從單相過(guò)冷水變?yōu)闅庖簝上嗷旌衔?，最終匯集到上聯(lián)箱中。流道沿軸向依次為入口段、加熱段和上升段，其中入口段和上升段均設(shè)有節(jié)流件模擬體，以模擬實(shí)際管路中的閥門(mén)、彎頭等裝置帶來(lái)的節(jié)流效應(yīng)。通道基本幾何參數(shù)列于表1。

表1　并聯(lián)雙通道基本幾何參數(shù)

1.2基本假設(shè)

本文中，數(shù)學(xué)模型的構(gòu)建基于下述假設(shè)：1) 流動(dòng)是一維的；2) 兩相流動(dòng)傳熱采用漂移流模型；3) 兩個(gè)實(shí)驗(yàn)段進(jìn)出口之間的壓差相等；4) 沿軸線方向計(jì)算流體的物性參數(shù)時(shí)，均以進(jìn)口壓力為基準(zhǔn)；5) 熱負(fù)荷沿管長(zhǎng)均勻分布；6) 忽略加熱管壁熱容。

其中，假設(shè)2涉及到兩相模型的選擇，漂移流模型曾被多數(shù)研究者使用，其有效性得到了充分驗(yàn)證。假設(shè)4的提出是考慮到管間壓差遠(yuǎn)低于運(yùn)行壓力，因此在物性計(jì)算時(shí)將整個(gè)通道的壓力視為一致。假設(shè)5的提出是為了與實(shí)驗(yàn)條件保持一致，便于對(duì)數(shù)值模擬結(jié)果進(jìn)行驗(yàn)證。

1.3數(shù)學(xué)模型

本文中，數(shù)學(xué)模型主要包括單相和兩相基本守恒方程，以及用以封閉方程組的本構(gòu)關(guān)系式和兩相摩阻經(jīng)驗(yàn)關(guān)系式。

1) 單相守恒方程組

(1)

(2)

(3)

(4)

式中：ρ為流體密度；u為流體速度；p為壓力；f為流體的摩阻系數(shù)；De為等效直徑；h為流體焓；q為熱流密度；g為重力加速度。

2) 兩相守恒方程組

(5)

(6)

(7)

本構(gòu)關(guān)系式：

(8)

3) 空泡份額和兩相摩阻經(jīng)驗(yàn)關(guān)系式

在模型中，空泡份額采用Zuber推薦的經(jīng)驗(yàn)關(guān)系式進(jìn)行計(jì)算，空泡份額的分布系數(shù)C0采用針對(duì)矩形通道的關(guān)系式計(jì)算。在12 MPa以下，采用Chisholm兩相摩阻關(guān)系式，在12～16 MPa壓力范圍內(nèi)，采用中國(guó)核動(dòng)力研究設(shè)計(jì)院自擬關(guān)系式，此自擬關(guān)系式是針對(duì)高壓條件下矩形通道兩相摩阻提出的。

(1) Zuber空泡份額關(guān)系式[8]

(9)

式中：x為熱平衡干度；C0=1.4-0.4p/pcr，pcr為臨界壓力。

(2) Chisholm兩相摩阻關(guān)系式[8]

(10)

式中，Γ為Chisholm數(shù)。

(3) 中國(guó)核動(dòng)力研究設(shè)計(jì)院自擬高壓矩形通道兩相摩阻關(guān)系式[9]

(11)

式(11)的適用范圍為：12.1～15.7 MPa，520～2 500 kg/(m2·s)。

2模型求解

2.1網(wǎng)格劃分及微分方程組的離散

本文采用等長(zhǎng)度的計(jì)算控制體對(duì)計(jì)算區(qū)域進(jìn)行網(wǎng)格劃分，速度、焓、密度等物理量均位于控制體的中心。經(jīng)網(wǎng)格無(wú)關(guān)性分析，采用1 mm與0.5 mm網(wǎng)格時(shí)，對(duì)同一工況計(jì)算獲得的不穩(wěn)定界限功率的相對(duì)偏差不超過(guò)2%，故將每個(gè)控制體的空間步長(zhǎng)設(shè)為1 mm，時(shí)間步長(zhǎng)設(shè)為0.01 s，如圖2所示。基于時(shí)間和空間網(wǎng)格，對(duì)式(1)～(7)進(jìn)行差分離散，可獲得單相和兩相守恒方程的差分方程組。其中，對(duì)流項(xiàng)的差分采用一階迎風(fēng)差分格式。

圖2　網(wǎng)格劃分方法 Fig.2　Method of mesh division

2.2數(shù)值求解方法

數(shù)值求解的基本流程如圖3所示。根據(jù)流程圖，首先應(yīng)求出給定熱工工況下的穩(wěn)態(tài)解，然后將穩(wěn)態(tài)解作為瞬態(tài)計(jì)算的初始條件，引入擾動(dòng)后展開(kāi)瞬態(tài)計(jì)算。根據(jù)瞬態(tài)計(jì)算的結(jié)果，結(jié)合相應(yīng)的判定準(zhǔn)則，判斷流動(dòng)不穩(wěn)定是否發(fā)生，并決定是否需要改變加熱功率繼續(xù)運(yùn)算。

圖3　數(shù)值求解基本流程 Fig.3　Flow path of numerical solving

1) 小擾動(dòng)的添加方法

在穩(wěn)態(tài)計(jì)算完成之后、瞬態(tài)計(jì)算開(kāi)始之前，需引入外加小擾動(dòng)來(lái)打破流動(dòng)傳熱的平衡，以觀察系統(tǒng)對(duì)擾動(dòng)的響應(yīng)，從而判斷該參數(shù)下系統(tǒng)是否會(huì)發(fā)生流動(dòng)不穩(wěn)定。引入擾動(dòng)的方法有兩種，一是流量擾動(dòng)，二是功率擾動(dòng)。

流量擾動(dòng)的實(shí)施方法：假設(shè)某通道(如通道1)入口流量W在t=0時(shí)刻發(fā)生了一個(gè)有限幅值的擾動(dòng)，由W1變?yōu)閃1+δW，由于總?cè)肟诹髁坎蛔?，通?的流量由W2變?yōu)閃2+δW，且在t=Δt時(shí)刻，擾動(dòng)消失。

功率擾動(dòng)的實(shí)施方法：假設(shè)在t=0時(shí)刻，功率發(fā)生了偏離穩(wěn)態(tài)的有限幅值的擾動(dòng)，由Q變?yōu)镼′，并在持續(xù)一段時(shí)間(通常設(shè)為1 s)后，恢復(fù)至原功率。

2) 瞬態(tài)計(jì)算內(nèi)迭代中方程組的求解策略

目前，針對(duì)二維以上流動(dòng)傳熱問(wèn)題的數(shù)值計(jì)算，較成熟的求解方法有SIMPLE系列算法、PISO算法以及由兩者結(jié)合而成的SIMPISO算法等。這些算法大都將質(zhì)量、動(dòng)量和能量方程耦合在一起，并對(duì)計(jì)算域內(nèi)所有方程組聯(lián)立求解。對(duì)于一維流動(dòng)傳熱問(wèn)題，微分方程組被簡(jiǎn)化為拋物型方程，因此一方面可將質(zhì)量、能量方程與動(dòng)量方程解耦，即在每一時(shí)層的內(nèi)迭代計(jì)算中，先求出滿(mǎn)足連續(xù)性方程和能量方程的解，在內(nèi)迭代收斂后，再根據(jù)已求得的速度和焓來(lái)計(jì)算壓降；另一方面，無(wú)需對(duì)全場(chǎng)方程組聯(lián)立求解，而是采用圖4所示的沿流動(dòng)方向依次求解方程組。每完成一次從入口到出口的求解稱(chēng)為一輪迭代，每一輪迭代之前，首先假設(shè)計(jì)算區(qū)域內(nèi)的密度場(chǎng)分布，若一輪迭代之后，新求解得到的密度場(chǎng)與假設(shè)的密度場(chǎng)足夠接近(滿(mǎn)足收斂條件)，則認(rèn)為獲得了滿(mǎn)足連續(xù)性方程和能量方程的解，否則，從入口開(kāi)始進(jìn)行下一輪迭代，直到滿(mǎn)足收斂條件。

圖4　求解方程組的次序 Fig.4　Order to solve equations

3) 微分方程組的離散及求解

將偏微分方程組在時(shí)間和空間節(jié)點(diǎn)上離散，得到差分方程組。

單相區(qū)的差分方程組為：

(12)

(13)

(14)

兩相區(qū)的差分方程組為：

(15)

(16)

(17)

方程迭代順序如圖4所示，具體求解步驟如下：

(3) 用計(jì)算出的流體焓(hi)和系統(tǒng)壓力更新整個(gè)區(qū)域內(nèi)各節(jié)點(diǎn)的流體物性(ρi，μi)。

(5) 將通過(guò)步驟1～4求解出的滿(mǎn)足連續(xù)性方程和能量守恒方程的速度(ui)、密度(ρi)和焓(hi)，代入動(dòng)量守恒方程計(jì)算出各節(jié)點(diǎn)之間的壓降(Δpi)，再對(duì)分壓降求和以獲得流體在整個(gè)流域上的總壓降。

4) 流量的動(dòng)態(tài)分配

上述算法是針對(duì)單個(gè)通道的，對(duì)于并聯(lián)通道中的各分通道，在計(jì)算的任意時(shí)刻入口焓是恒定且已知的，但入口速度卻是隨時(shí)間變化的，若無(wú)法確定入口速度，則內(nèi)迭代就無(wú)法按照前文描述的方式進(jìn)行。欲計(jì)算分通道入口流速，需利用兩個(gè)邊界條件，一是流量邊界條件，即入口總流速恒定且已知，二是壓力邊界條件，即并聯(lián)通道的分通道總壓降相等。

基于上述邊界條件，可采用下述方法求解并聯(lián)通道流量動(dòng)態(tài)分配問(wèn)題。

(1) 令p′=Δp1-Δp2，則p′為入口流量的函數(shù)，可記為p′(W1，W2)，又由于W1+W2=W，W為已知量，則p′(W1，W2)可變?yōu)閜′(W1)。

(2) 令W(n)為第n次迭代值，則由弦截法可得：

(18)

又知p′(W1(n))=0，因此，可得：

(19)

(3) 反復(fù)迭代，直至p′(W(n))<ε，取W1=W1(n)、W2=W-W1。

5) 流動(dòng)不穩(wěn)定的判斷方法

實(shí)驗(yàn)中通過(guò)觀察流量脈動(dòng)幅值是否放大來(lái)判斷流動(dòng)不穩(wěn)定是否發(fā)生，在數(shù)值模擬中，也可采用相似的方法來(lái)尋找流動(dòng)不穩(wěn)定邊界。事實(shí)上，在計(jì)算過(guò)程中，無(wú)論熱流密度是否達(dá)到發(fā)生流動(dòng)不穩(wěn)定的界限熱流密度，當(dāng)小擾動(dòng)被引入后，入口流量都會(huì)出現(xiàn)周期性振蕩，不同之處在于振蕩幅值是逐漸衰減、不斷放大還是恒定不變。顯然，脈動(dòng)幅值逐漸減小意味著流量振蕩終會(huì)消失，流動(dòng)仍處于穩(wěn)定狀態(tài)，脈動(dòng)幅值不斷放大意味著流動(dòng)已失穩(wěn)，而流量呈等幅振蕩則標(biāo)志著流動(dòng)正處于穩(wěn)定和失穩(wěn)的邊界，在這種情況下所獲得的熱流密度即為流動(dòng)不穩(wěn)定界限

熱流密度。

3計(jì)算結(jié)果

3.1流動(dòng)不穩(wěn)定基本現(xiàn)象

通過(guò)數(shù)值模擬，獲得了入口流量振蕩的基本圖像，如圖5所示。并聯(lián)通道中發(fā)生的是典型的周期性異相脈動(dòng)，且由振蕩模式的不同可分辨出圖5中流動(dòng)分別處于穩(wěn)定、失穩(wěn)和不穩(wěn)定邊界。

圖6示出流動(dòng)失穩(wěn)后入口流量與進(jìn)出口壓降的瞬時(shí)變化情況?？煽闯?，各分通道的入口壓降與各自通道入口流量的變化同相，而出口壓降則與各自通道入口流量的變化反相。這與通過(guò)實(shí)驗(yàn)所獲得的結(jié)論完全一致。

3.2不穩(wěn)定界限熱流密度及不穩(wěn)定邊界計(jì)算

表2為所選擇的16組熱工工況，均來(lái)自并聯(lián)矩形雙通道密度波流動(dòng)不穩(wěn)定性實(shí)驗(yàn)[10]，實(shí)驗(yàn)中所采用矩形通道的基本幾何參數(shù)與表1中完全一致。實(shí)驗(yàn)參數(shù)范圍為：系統(tǒng)壓力，3～14 MPa；質(zhì)量流速，400～800 kg/(m2·s)；入口過(guò)冷度，5～120 ℃。

a——幅值衰減(流動(dòng)穩(wěn)定)；b——幅值放大(流動(dòng)失穩(wěn))；c——等幅振蕩(不穩(wěn)定邊界) 圖5　流量振蕩的3種模式 Fig.5　Three oscillation modes of flow

圖6　流動(dòng)失穩(wěn)后入口流量與進(jìn)出口壓降的瞬時(shí)變化 Fig.6　Transient variation of inlet flow and pressure drop at inlet and exit of channel during flow oscillation

表2　計(jì)算工況

通過(guò)計(jì)算，獲得了16組工況的不穩(wěn)定界限熱流密度，并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對(duì)比，具體結(jié)果列于表3。從表3可見(jiàn)，除工況1和工況5外，其余工況下流動(dòng)不穩(wěn)定界限熱流密度的計(jì)算值與實(shí)驗(yàn)值的相對(duì)偏差在±20%以?xún)?nèi)，說(shuō)明數(shù)值計(jì)算的結(jié)果具有一定的準(zhǔn)確性。此外，從表3還可看到，在13 MPa的系統(tǒng)壓力下，計(jì)算相對(duì)偏差小于±6%，這是因?yàn)樵诟邏合虏捎昧酸槍?duì)矩形通道的兩相摩擦壓降關(guān)系式之故。

按式(14)、(15)將表2、3中熱工參數(shù)整理成過(guò)冷度數(shù)Nsub和相變數(shù)Npch，獲得了不同壓力下的計(jì)算流動(dòng)不穩(wěn)定邊界，其中流動(dòng)不穩(wěn)定邊界采用Isshi[11]推薦的過(guò)冷度數(shù)-相變數(shù)空間進(jìn)行描述，如圖7所示。與通過(guò)實(shí)驗(yàn)獲得的流動(dòng)不穩(wěn)定邊界相比，計(jì)算不穩(wěn)定邊界在13 MPa下與實(shí)驗(yàn)值吻合較好；10 MPa下，計(jì)算值與實(shí)驗(yàn)值在低過(guò)冷度區(qū)域吻合較好，在高過(guò)冷度區(qū)域則偏差較大，且隨著過(guò)冷度的增加有放大趨勢(shì)；5 MPa下，計(jì)算不穩(wěn)定邊界偏保守。

表3　不穩(wěn)定界限熱流密度計(jì)算值與實(shí)驗(yàn)值的對(duì)比

圖7　不穩(wěn)定邊界計(jì)算值和實(shí)驗(yàn)值對(duì)比 Fig.7　Comparison between calculated flow instability boundary and experimental data

(20)

(21)

式中：Δhsub為入口欠焓；Δρfg為氣液相密度差；Ah為加熱面積；A為流通面積。

綜上可知，數(shù)值模擬的計(jì)算結(jié)果在高壓下與實(shí)驗(yàn)值吻合得更好。這是因?yàn)樵?2 MPa以上的工況下采用了針對(duì)矩形通道的摩阻關(guān)系式，而在對(duì)12 MPa以下工況進(jìn)行計(jì)算時(shí)，由于缺乏相應(yīng)的矩形通道兩相摩阻關(guān)系式，而采用了基于圓管的經(jīng)驗(yàn)關(guān)系式，使得計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果存在偏差。

4結(jié)論

本文基于一維漂移流模型構(gòu)建了并聯(lián)雙通道密度波流動(dòng)不穩(wěn)定性數(shù)學(xué)模型。模型中采用Zuber推薦的經(jīng)驗(yàn)關(guān)系式計(jì)算流體空泡份額，分別采用Chisholm關(guān)系式和中國(guó)核動(dòng)力研究設(shè)計(jì)院自擬關(guān)系式計(jì)算12 MPa以下和12 MPa以上兩相流體的摩擦壓降。通過(guò)數(shù)值模擬，獲得了與實(shí)驗(yàn)現(xiàn)象特征一致的密度波脈動(dòng)圖像。選擇16組實(shí)驗(yàn)工況進(jìn)行數(shù)值模擬，計(jì)算結(jié)果表明，在所有選擇計(jì)算工況中，除個(gè)別工況外，通過(guò)計(jì)算獲得的流動(dòng)不穩(wěn)定界限熱流密度與實(shí)驗(yàn)值的相對(duì)偏差小于±20%。尤其是在13 MPa下，其相對(duì)偏差小于±6%，這表明，采用本文構(gòu)建的數(shù)學(xué)模型和提出的數(shù)值求解方法可對(duì)并聯(lián)矩形雙通道流動(dòng)不穩(wěn)定性進(jìn)行較為準(zhǔn)確的模擬。

參考文獻(xiàn)：

[1]魯鐘琪. 兩相流與沸騰傳熱[M]. 北京：清華大學(xué)出版社，2002.

[2]MEYER J, ROSE R. Application of a momentum integral model to the study of parallel channel boiling flow oscillations[J]. J Heat Transfer, 1963, 85(1): 1-9.

[3]DOGAN T, KAKAC S, VEZIROGLU T N. Lumped parameter analysis of two-phase flow instabilities[C]∥Proceedings of the 7th International Heat Transfer Conference. Munchen: [s. n.], 1982: 213-218.

[4]LECOQ G, METAICH M, SLASSI-SENNOU M. A simple model for the study of dynamic instabilities[J]. Nuclear Engineering and Design, 1990, 122: 41-52.

[5]黃軍，黃彥平，王飛. 兩管平行通道管間脈動(dòng)建模及程序編制[C]∥中國(guó)工程熱物理學(xué)會(huì)多相流學(xué)術(shù)會(huì)議年會(huì)論文集. 成都：中國(guó)工程熱物理學(xué)會(huì)，2004：157-162.

[6]李虹波. 矩形雙通道流動(dòng)不穩(wěn)定性基礎(chǔ)研究[D]. 成都：中國(guó)核動(dòng)力研究設(shè)計(jì)院，2006.

[7]夏庚磊，郭贅，彭敏俊. 基于RELAP5的兩管平行通道流動(dòng)不穩(wěn)定性研究[J]. 原子能科學(xué)技術(shù)，2010，44(6)：694-700.

XIA Genglei, GUO Yun, PENG Minjun. Investigation on two-phase flow instability in parallel channels based on RELAP5 code[J]. Atomic Energy Science and Technology, 2010, 44(6): 694-700(in Chinese).

[8]徐濟(jì)鋆，賈斗南. 沸騰傳熱和汽液兩相流[M]. 北京：原子能出版社，1993.

[9]洪鋼. 矩形通道兩相摩阻實(shí)驗(yàn)研究[R]. 成都：中國(guó)核動(dòng)力研究設(shè)計(jì)院，2009.

[10]唐瑜. 并聯(lián)矩形雙通道第一二類(lèi)密度波流動(dòng)不穩(wěn)定性研究[D]. 成都：中國(guó)核動(dòng)力研究設(shè)計(jì)院，2014.

[11]ISHII M. Thermal induced flow instabilities in two-phase mixtures in thermal equilibrium[D]. Georgia, US: Georgia Institute of Technology, 1971.