基于一維漂移流模型的并聯矩形雙通道密度波流動不穩定性數值模擬

基于一維漂移流模型的并聯矩形雙通道密度波流動不穩定性數值模擬

熊萬玉，唐瑜*，陳炳德

(中國核動力研究設計院，四川 成都610041)

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

關鍵詞：并聯矩形雙通道；密度波不穩定性；漂移流模型；數值模擬

中圖分類號：TK124 文獻標志碼：A

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

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

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

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

本文基于一維漂移流模型構建并聯矩形雙通道密度波流動不穩定性數學模型，給出該模型的數值求解方法，選定16組計算工況開展流動不穩定性數值模擬，獲得各工況下流動不穩定的界限熱流密度，繪制相應的流動不穩定邊界，并與實驗數據進行對比。

1數學物理模型

圖1　并聯雙通道幾何模型 Fig.1　Geometry model of parallel twin channels

1.1幾何模型

并聯雙通道的幾何結構如圖1所示。圖1中，兩個結構完全相同的矩形通道分別與上下聯箱相連接后構成并聯通道，在聯箱內部流體的壓力相等。流體從下聯箱進入裝置后分別經由兩個通道垂直向上流動，在流動過程中受到管壁的加熱后從單相過冷水變為氣液兩相混合物，最終匯集到上聯箱中。流道沿軸向依次為入口段、加熱段和上升段，其中入口段和上升段均設有節流件模擬體，以模擬實際管路中的閥門、彎頭等裝置帶來的節流效應。通道基本幾何參數列于表1。

表1　并聯雙通道基本幾何參數

1.2基本假設

本文中，數學模型的構建基于下述假設：1) 流動是一維的；2) 兩相流動傳熱采用漂移流模型；3) 兩個實驗段進出口之間的壓差相等；4) 沿軸線方向計算流體的物性參數時，均以進口壓力為基準；5) 熱負荷沿管長均勻分布；6) 忽略加熱管壁熱容。

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

1.3數學模型

本文中，數學模型主要包括單相和兩相基本守恒方程，以及用以封閉方程組的本構關系式和兩相摩阻經驗關系式。

1) 單相守恒方程組

(1)

(2)

(3)

(4)

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

2) 兩相守恒方程組

(5)

(6)

(7)

本構關系式：

(8)

3) 空泡份額和兩相摩阻經驗關系式

在模型中，空泡份額采用Zuber推薦的經驗關系式進行計算，空泡份額的分布系數C0采用針對矩形通道的關系式計算。在12 MPa以下，采用Chisholm兩相摩阻關系式，在12～16 MPa壓力范圍內，采用中國核動力研究設計院自擬關系式，此自擬關系式是針對高壓條件下矩形通道兩相摩阻提出的。

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

(9)

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

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

(10)

式中，Γ為Chisholm數。

(3) 中國核動力研究設計院自擬高壓矩形通道兩相摩阻關系式[9]

(11)

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

2模型求解

2.1網格劃分及微分方程組的離散

本文采用等長度的計算控制體對計算區域進行網格劃分，速度、焓、密度等物理量均位于控制體的中心。經網格無關性分析，采用1 mm與0.5 mm網格時，對同一工況計算獲得的不穩定界限功率的相對偏差不超過2%，故將每個控制體的空間步長設為1 mm，時間步長設為0.01 s，如圖2所示?；跁r間和空間網格，對式(1)～(7)進行差分離散，可獲得單相和兩相守恒方程的差分方程組。其中，對流項的差分采用一階迎風差分格式。

圖2　網格劃分方法 Fig.2　Method of mesh division

2.2數值求解方法

數值求解的基本流程如圖3所示。根據流程圖，首先應求出給定熱工工況下的穩態解，然后將穩態解作為瞬態計算的初始條件，引入擾動后展開瞬態計算。根據瞬態計算的結果，結合相應的判定準則，判斷流動不穩定是否發生，并決定是否需要改變加熱功率繼續運算。

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

1) 小擾動的添加方法

在穩態計算完成之后、瞬態計算開始之前，需引入外加小擾動來打破流動傳熱的平衡，以觀察系統對擾動的響應，從而判斷該參數下系統是否會發生流動不穩定。引入擾動的方法有兩種，一是流量擾動，二是功率擾動。

流量擾動的實施方法：假設某通道(如通道1)入口流量W在t=0時刻發生了一個有限幅值的擾動，由W1變為W1+δW，由于總入口流量不變，通道2的流量由W2變為W2+δW，且在t=Δt時刻，擾動消失。

功率擾動的實施方法：假設在t=0時刻，功率發生了偏離穩態的有限幅值的擾動，由Q變為Q′，并在持續一段時間(通常設為1 s)后，恢復至原功率。

2) 瞬態計算內迭代中方程組的求解策略

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

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

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

將偏微分方程組在時間和空間節點上離散，得到差分方程組。

單相區的差分方程組為：

(12)

(13)

(14)

兩相區的差分方程組為：

(15)

(16)

(17)

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

(3) 用計算出的流體焓(hi)和系統壓力更新整個區域內各節點的流體物性(ρi，μi)。

(5) 將通過步驟1～4求解出的滿足連續性方程和能量守恒方程的速度(ui)、密度(ρi)和焓(hi)，代入動量守恒方程計算出各節點之間的壓降(Δpi)，再對分壓降求和以獲得流體在整個流域上的總壓降。

4) 流量的動態分配

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

基于上述邊界條件，可采用下述方法求解并聯通道流量動態分配問題。

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

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

(18)

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

(19)

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

5) 流動不穩定的判斷方法

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

熱流密度。

3計算結果

3.1流動不穩定基本現象

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

圖6示出流動失穩后入口流量與進出口壓降的瞬時變化情況。可看出，各分通道的入口壓降與各自通道入口流量的變化同相，而出口壓降則與各自通道入口流量的變化反相。這與通過實驗所獲得的結論完全一致。

3.2不穩定界限熱流密度及不穩定邊界計算

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

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

圖6　流動失穩后入口流量與進出口壓降的瞬時變化 Fig.6　Transient variation of inlet flow and pressure drop at inlet and exit of channel during flow oscillation

表2　計算工況

通過計算，獲得了16組工況的不穩定界限熱流密度，并與實驗數據進行了對比，具體結果列于表3。從表3可見，除工況1和工況5外，其余工況下流動不穩定界限熱流密度的計算值與實驗值的相對偏差在±20%以內，說明數值計算的結果具有一定的準確性。此外，從表3還可看到，在13 MPa的系統壓力下，計算相對偏差小于±6%，這是因為在高壓下采用了針對矩形通道的兩相摩擦壓降關系式之故。

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

表3　不穩定界限熱流密度計算值與實驗值的對比

圖7　不穩定邊界計算值和實驗值對比 Fig.7　Comparison between calculated flow instability boundary and experimental data

(20)

(21)

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

綜上可知，數值模擬的計算結果在高壓下與實驗值吻合得更好。這是因為在12 MPa以上的工況下采用了針對矩形通道的摩阻關系式，而在對12 MPa以下工況進行計算時，由于缺乏相應的矩形通道兩相摩阻關系式，而采用了基于圓管的經驗關系式，使得計算結果與實驗結果存在偏差。

4結論

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

參考文獻：

[1]魯鐘琪. 兩相流與沸騰傳熱[M]. 北京：清華大學出版社，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]黃軍，黃彥平，王飛. 兩管平行通道管間脈動建模及程序編制[C]∥中國工程熱物理學會多相流學術會議年會論文集. 成都：中國工程熱物理學會，2004：157-162.

[6]李虹波. 矩形雙通道流動不穩定性基礎研究[D]. 成都：中國核動力研究設計院，2006.

[7]夏庚磊，郭贅，彭敏俊. 基于RELAP5的兩管平行通道流動不穩定性研究[J]. 原子能科學技術，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]徐濟鋆，賈斗南. 沸騰傳熱和汽液兩相流[M]. 北京：原子能出版社，1993.

[9]洪鋼. 矩形通道兩相摩阻實驗研究[R]. 成都：中國核動力研究設計院，2009.

[10]唐瑜. 并聯矩形雙通道第一二類密度波流動不穩定性研究[D]. 成都：中國核動力研究設計院，2014.

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