蔡 溢 林 萌 張 昊 楊燕華

1（上海交通大學機械與動力工程學院 上海 200240）

2（國家核電軟件中心 北京 100029）

COSINE程序[1]是國家重大專項核電關(guān)鍵設(shè)計軟件自主化技術(shù)研究中，正在開發(fā)的反應(yīng)堆設(shè)計與分析程序包，其中包括大型的熱工水力安全分析與設(shè)計程序。針對水堆的熱工安全分析通常需要求解由質(zhì)量、動量、能量守恒組成的基本方程組。其中，由于水存在著相變，需要分析和計算兩相流動過程，在某些工況下，均勻流模型有一定的精度，有利于從宏觀上認識兩相流動過程，但對于較嚴格的設(shè)計計算及安全分析則需要應(yīng)用更精確的兩相流模型，如漂移流模型、兩流體模型等。本文以開發(fā)兩相流程序為背景，針對熱工程序模型的汽液相間摩擦部分，建立了一個考慮因素全面、分區(qū)詳細的理論模型，編制和開發(fā)了相應(yīng)的程序，并選擇適當?shù)乃憷訰ELAP5/Mod3.2程序計算結(jié)果為參考進行驗證，理論上證明了該相間摩擦模型可以滿足軟件開發(fā)的需求。

1 模型介紹

1.1 基本方程

兩流體熱工程序的本構(gòu)方程是以流體動力學的質(zhì)量連續(xù)性、動量守恒和能量守恒這三個基本方程為基礎(chǔ)，相間摩擦相出現(xiàn)在動量守恒方程中，如下：

氣相：

液相：

將動量方程進行有限差分之后可以得到相間摩擦力的公式：

根據(jù)FI的定義：

以及相間摩擦力：

有：

在進行相間作用的計算時一般假設(shè)氣液兩相間存在一個相界面，該界面不能儲存能量，并且液相對相界面的摩擦力與氣相對相界面的摩擦力相等，所以式(6)可以寫成：

界面力 Fi和兩相間的相對速度vR的計算采用漂移流模型和拖曳力模型兩種不同的方式。對于豎直泡狀流和豎直彈狀流的計算，漂移流模型的精度較高，工程應(yīng)用最為廣泛，且fx=1。拖曳力系數(shù)模型用于除了豎直泡狀流和豎直彈狀流以外的其它所有流型，且fx=0。

1.1.1 漂移流模型

1965年，Zuber等[2]提出了漂移流模型，假設(shè)兩相間處于熱力學平衡狀態(tài)，氣液兩相流速各自保持不變但不相等。

根據(jù)動量方程及漂移流模型原理，可以得到氣液兩相的相間阻力：

式中，g為重力加速度；fj為連接件的傾斜角；vgj為氣液兩相的漂移速度。

目前基于漂移流原理的計算方法主要有：Zuber-Findlay等[2?6]。不同計算式之間的差別僅為由實驗確定的C0和vgj表達式的不同。具體公式的選擇與所給流型以及流道的幾何參數(shù)有關(guān)，選擇依據(jù)參考了Putney[7]的理論。

1.1.2 拖曳力系數(shù)模型

物體相對于液體移動時的拖曳力的基本方程為：

式中，F(xiàn)為拖曳力；r為液體密度；v為物體相對液體的速度；CD為拖曳力系數(shù)；A為物體的投影面積。

如果要描述一組物體相對于液體運動（如：氣泡從液體中流過或者液滴從氣體或液體中流過）時單一物體產(chǎn)生的摩擦力，則有：

式中，rc是連續(xù)相的密度；agf為單位體積下的相界面面積；SF為形狀因子，取1.0。

參照漂移流模型的分析則有：

1.1.3 相間面積的計算

在RELAP5/MOD3.2程序中對于單位體積相界面面積 agf的計算一方面參考了已發(fā)表的理論和實驗結(jié)果，另一方面也參考了RELAP5/Mod2程序的模擬測試結(jié)果而做出一定的優(yōu)化修正。這樣的處理方式雖然能夠很好地滿足計算需求，但是添加的很多修正因子限制了模型的適用范圍，并且使模型變得復雜和計算變得繁瑣。在對物理現(xiàn)象進行適當?shù)募僭O(shè)之后就能建立簡單的模型對相間面積進行計算，并且能夠滿足精度需求。

以環(huán)霧狀流為例，RELAP5/Mod3.2模型采用的計算公式參考了Ishii和Mishima[8]理論，比較繁瑣。

假設(shè)液滴為球狀，且環(huán)霧狀流中夾帶液滴的體積份額遠小于氣相的混合物，則有：

所以有如下相間面積計算公式：

因此相間面積的計算簡略到只與水力學直徑和空泡份額有關(guān)。

選擇工況：壓力0.1 MPa，溫度373.95 oC，液相速度4 m·s?1，氣相速度2 m·s?1，以空泡份額為變量將該模型與 RELAP5/Mod3.2中采用的模型進行計算結(jié)果如圖1所示。

由圖1可以看出，兩者計算結(jié)果相差不大，但是新的模型比 RELAP5/Mod3.2的模型更簡潔，計算更方便，所以采用新的模型取代了RELAP5/Mod3.2的相間面積計算模型。

1.2 模型選擇

在進行相間摩擦的計算之前，先要有相應(yīng)的流型圖進行流型判別，在確定流型之后再進行模型的選擇。當流道的傾角為0°?30°時采用水平流型圖，當流道為 60°?90°時用豎直流型圖，方位角在兩者之間時使用相應(yīng)的水平流型計算結(jié)果和豎直流型計算結(jié)果進行插值。具體流型劃分見圖2。

圖1 相間面積計算結(jié)果圖Fig.1 Calculation results of interfacial area.

圖2 水平(a)和豎直(b)流型判斷流程Fig.2 Schematic of horizontal (a) and vertical (b) flow regime.

根據(jù)以上流程判斷出流道中流體的流型之后則可進入相應(yīng)的模型選擇適當?shù)慕?jīng)驗關(guān)系式對 CD、C0和vgj進行計算，從而得到相間摩擦的結(jié)果。

2 模型驗證

根據(jù)以上模型我們編寫了COSINE程序的相間摩擦部分，并對該部分程序進行計算分析，這里具體以如下算例為例對程序進行驗證。以單一豎直管道為研究對象，將管道分為十個控制體，如圖3所示。水和空氣的混合氣體豎直向下流動，入口條件為：壓力0.1 MPa，空泡份額0.85，溫度為373.95 oC，液相速度為4 m·s?1，氣相速度為2 m·s?1；出口邊界條件為373.95 oC?？紤]重力的影響，不考慮壁面摩擦。在此條件下流體為豎直環(huán)狀流，采用拖曳力模型進行計算。

圖3 算例節(jié)點圖Fig.3 Nodalization diagram.

由于本測試選用的相間模型參考了RELAP5/Mod3.2的模型，所以，將COSINE程序的計算結(jié)果（圖4中標識COSINE）與RELAP5/Mod3.2程序的計算結(jié)果（圖4中標識RELAP5/Mod3.2）進行對比分析。同時，為了更加直觀地驗證COSINE程序相間摩擦模型的正確性，將 RELAP5/Mod3.2程序計算的相間摩擦系數(shù)代入COSINE程序模型進行計算（圖4中標識COSINE_RELAP）。三者的計算結(jié)果如圖4所示。圖4(a)表示管道中壓力隨位置的變化，流體垂直向下流動，由于重力的影響壓力呈現(xiàn)逐漸增大的趨勢。圖4(b)和(c)表示管道中兩相流體的速度變化，對于液體，由于重力的作用使之加速，所以速度越來越大。而液相重力的作用使得壓力增大，對于氣相，進口、出口之間的壓差作用的影響遠大于氣相本身的重力作用以及氣液兩相間的摩擦力作用，它將對氣相產(chǎn)生一定的阻力。雖然，在液相的影響下氣相一直向下運動，但是，其加速度為負數(shù)，這將使其向下運動的速度越來越小。圖4(d)表示流體中兩相間的摩擦變化，從圖 4(b)和(c)可以得到，由于在入口x=0處兩相間速度差比較小，越往下流動速度差越大，所以摩擦呈現(xiàn)逐漸變大的趨勢。

圖4 壓力(a)、氣相速度(b)、液相速度(c)和相間摩擦系數(shù)(d)隨位置變化圖Fig.4 Pressure (a), gas velocity (b), liquid velocity (c) and interphase friction coefficient (d) vs. location.

從RELAP5/MOD3.2、COSINE以及直接使用確定摩擦系數(shù)的COSINE_RELAP程序的計算結(jié)果對比發(fā)現(xiàn)，在采用同樣的相間摩擦模型基礎(chǔ)上，由于COSINE程序采用不同的求解器和數(shù)值算法，它在定量上出現(xiàn)非常小的偏差。但是，COSINE程序整體在計算該算例時能夠很順利地運行并且計算結(jié)果與 RELAP5/Mod3.2程序的計算結(jié)果基本一致，它滿足當前階段程序開發(fā)的要求。

3 結(jié)語

對于兩流體程序，相間摩擦計算模型是其不可缺少的重要的分析模型之一。本文介紹了COSINE大型熱工安全分析程序的相間摩擦模型的主要建模過程，該模型主要根據(jù)流體的流型選擇漂移流模型和拖曳力模型，在確定選用模型之后，又根據(jù)流型以及流道的不同選擇相應(yīng)的經(jīng)驗關(guān)系式對相間摩擦進行計算，并以一個簡單算例進行模型驗證，將計算結(jié)果與RELAP5程序的計算結(jié)果進行對比，由此證明COSINE程序可很好地完成相間摩擦的計算。需要注意的是，在進行過渡區(qū)域的相間摩擦計算時，需要選擇合理的插值函數(shù)和因子。否則將會對計算的連續(xù)性有較大影響，可能會影響程序的正常運行。

1 Yang Y H, Chen Y X. Requirement analysis and primary design of COSINE code[R]. ANS Annual Meeting, Chicago, IL. Hyatt Regency Chicago, No.5805, June, 2012

2 Ishii M. One-dimensinal drift-flux model and constitutive equations for relative motion between phases in various tow-phase flow regimes[R]. ANL-77-47, USA, 1977

3 Zuber N, Findlay J. Average volumetric concentration in tow-phase flow systems[J]. Journal of Heat Transfer, 1965, 85(1): 453?568

4 Chexal B, Lellouche G. A full-range drift-flux correlation for vertical flows (Rev. 1)[R]. EPRI NP-3989-SR, September 1986

5 Kataoka I, Ishii M. Drift-flux model for large diameter pipe and new correlation for pool void fraction[J]. International Journal of Heat and Mass Transfer, 1987, 30(9): 1927?1939

6 Hibiki T, Ishii M. One-dimensional drift-flux model for tow-phase flow in a large diameter pipe[J]. International Journal of Heat and Mass Transfer, 2003, 46(10): 1773?1790

7 Putney J. Development of a new bubbly-slug interfacial friction model for RELAP5[J]. Nuclear Engineering and Design, 1991, 131(1): 223?240

8 Ishii M, Mishima K. Study of two-fluid model and interfacial area[R]. NUREG/CR-1873, ANL-80-111, Argonne National Laboratory, December 1980