趙柏陽 李 冬 楊已顥 張永祺

（上海電力大學 能源與機械工程學院 上海 200090）

氣液相向流動限制現(xiàn)象（Counter-Current Flow Limitation，CCFL）［1］對于壓水堆核電廠的安全分析十分重要，失水事故［2］后，當向上流動的氣相速度足夠大時，會部分或完全地阻礙液相向下流動。該現(xiàn)象機理復雜，受汽液兩相流速、流道幾何結(jié)構(gòu)、壓力空泡份額等多種因素影響，在第三代核電技術(shù)發(fā)展以前，研究主要集中于對熱管段、壓力容器下降段等豎直管或傾斜管內(nèi)的實驗及理論分析。Wang［3］、Kim［4］和Navarro［5］等研究了管徑、流速、傾斜角等因素對直管、傾斜管中CCFL現(xiàn)象的影響，并通過實驗數(shù)據(jù)擬合出公式用于預測CCFL現(xiàn)象。近年來，由于自動泄壓系統(tǒng)的引入，事故后穩(wěn)壓器波動管中的CCFL現(xiàn)象增加了安全風險。穩(wěn)壓器波動管結(jié)構(gòu)則更為復雜，同時包含了豎直管、傾斜管及螺旋彎管，一維系統(tǒng)程序的模型研究對復雜三維結(jié)構(gòu)并不足夠。Takeuchi等［6］、Doi［7］通過實驗對穩(wěn)壓器波動管CCFL現(xiàn)象進行研究發(fā)現(xiàn)，發(fā)生CCFL現(xiàn)象最劇烈的位置位于豎直管段，而對彎管內(nèi)的現(xiàn)象分析仍不明確。Murase［8］使用傳統(tǒng)流動模型對 Futatsugi［9］的實驗臺架進行了三維數(shù)值模擬，模擬結(jié)果與實驗數(shù)據(jù)相比存在一定偏差。因此有必要對穩(wěn)壓器波動管的CCFL現(xiàn)象的相間作用模型及流動特性進行更深入的研究。

基于CCFL現(xiàn)象機理及汽液相間作用的影響，本文使用ANSYS CFX程序，選擇自由表面模型結(jié)合（Algebraic Interfacial Area Density，AIAD）模型對穩(wěn)壓器波動管CCFL現(xiàn)象進行三維數(shù)值模擬，并基于網(wǎng)格影響對AIAD模型進行了修正，通過模擬結(jié)果和實驗現(xiàn)象的對比來驗證模型，并進一步探究穩(wěn)壓器波動管的CCFL特性以及幾何結(jié)構(gòu)的影響。

1 計算模型

1.1 自由表面模型

CCFL現(xiàn)象發(fā)生時常伴有阻塞形成、液滴夾帶、水力跳躍等現(xiàn)象，這些現(xiàn)象會引起相間作用力的變化。相間作用可以通過界面剪切應力描述，不同流型下的氣液兩相之間的界面剪切應力可以由自由表面模型求得。

在三維模擬計算中，求解了兩流體模型的質(zhì)量方程和動量方程，其形式如下：

式中：下標k代表氣相或液相；ρ為密度；U為速度矢量；t為時間；p為壓力；g為重力加速度；r為體積分數(shù)；τ為剪切應力（τυ為平均粘性剪應力τt為湍流剪應力）；τD為界面剪應力。式（2）等號右邊第一項為壓力項，第二項為重力項，第三項為平均粘性剪切應力和湍流剪切應力項，第四項為界面剪切應力項。自由表面模型主要通過求解界面剪切應力來描述相間作用關(guān)系。

界面剪應力τD可表示為流體流動過程中總的阻力與界面面積密度的比：

式中：FD為流體流動過程中的總阻力；A是界面面積密度。

式中：ρLG是平均密度；|UL-UG|是相對速度；CD是無量綱的界面摩擦阻力系數(shù)。

在自由表面流動的模擬中，式（4）不能代表一個真實的物理模型。假設在界面附近的兩種流體的速度是相似的。因此從兩相表面附近假定一個類似壁面剪應力的剪應力，以減小兩相的速度差，沿“固體”狀邊界運動的粘性流體會產(chǎn)生剪切應力，在無滑移條件下，自由表面區(qū)域為邊界層，剪切應力的計算表達式可寫成以下形式：

式中：AFS為自由表面的界面面積密度。對于氣泡和液滴來說，有固定的公式來描述界面面積密度，但是在實際流動過程中流態(tài)常常介于氣泡和液滴之間，在這里我們加入AIAD模型根據(jù)不同氣液相體積分數(shù)來代數(shù)化求解AFS，這樣可以通過氣液相體積分數(shù)的變化得到不同的界面摩擦系數(shù)，更好地描述流動的特征。

1.2 修正的AIAD模型

原始AIAD模型［10］采用了三種不同的阻力系數(shù)：下標B、D、FS分別表示氣泡、水滴、自由表面。界面面積密度A也與相的形態(tài)有關(guān)。對于泡狀流型：

式中：rG為氣相體積分數(shù)；dB為氣泡直徑。

對于自由表面：

根據(jù)式（3）、（4）、（7）、（9）得到：

界面面積密度的簡單切換過程使用混合函數(shù)f。該函數(shù)引入空泡份額率極限、流態(tài)權(quán)值以及氣泡和液滴流動的長度尺度，定義為：

界面面積密度和阻力系數(shù)分別定義為：

根據(jù)經(jīng)驗值αD，limit=αB，limit=0.3，aD=aB=70。CD，B=CD，D=0.44。

本文在傳統(tǒng)AIAD模型的基礎上對混合函數(shù)f的表達式進行了修正，目的是為了使混合多種影響因素后總的阻力系數(shù)隨著局部流動特征引起的不同參數(shù)而變化，這樣可以避免因為經(jīng)驗值選取不當而導致的計算偏差。為便于區(qū)分，函數(shù)用φ來表示。

式中：φm代表氣泡的混合函數(shù)表達式；αg，crit=0.3，CD，cont是變量系數(shù)；CD，bubb=CD，drop=0.44，|?αg|是空泡份額梯度；Δx表示單位網(wǎng)格尺寸；這里aFS取100，n取4。

該模型可以通過體積分數(shù)來判斷出局部的流動形態(tài)［11-12］（表1）。

表1 體積分數(shù)和相應的形態(tài)Table 1 Volume fractions and corresponding morphologies

2 數(shù)值模型建立

2.1 幾何構(gòu)建和邊界條件設置

本文以文獻［9］中的1/10比例的AP600穩(wěn)壓器實驗為參考，建立對應的幾何模型。實驗裝置如由一個傾斜管（傾斜角θ=0.6°～5°）、一個豎直彎頭和一個垂直管組成。管徑D=30 mm，空氣入口速度JG，in范圍為 0～5.5 m·s-1，JL，in液相入口速度始終為0.07 m·s-1。三維模擬所建立的幾何模型尺寸均與實驗裝置相同。進出口邊界設置如圖1所示。

圖1 幾何模型及邊界設置Fig.1 Geometric model and boundary setting

三維數(shù)值模擬計算通過ANSYS CFX軟件來實現(xiàn)，數(shù)值計算在瞬態(tài)模式下運行。在網(wǎng)格獨立性研究中，對三種不同精度的網(wǎng)格進行了比較。對于最粗網(wǎng)格（網(wǎng)格數(shù)量為40 000節(jié)點）沒有達到收斂，對節(jié)點數(shù)量為130 000和500 000的網(wǎng)格的主流方向速度分布進行了分析，結(jié)果表明：兩種數(shù)量的網(wǎng)格主流方向速度分布沒有定性差異，定量偏差最大小于1%，因此使用節(jié)點數(shù)量為130 000的網(wǎng)格來完成計算。在建模中選用SST浮力湍流模型，對流項采用高階離散格式，對于時間積分，采用完全隱式二階倒推歐拉法，時間步長Δt=0.001 s，每個時間步長最大循環(huán)次數(shù)為10次。收斂性定義為殘差的RMS的值要小于10-4。

氣液相是等溫不可壓縮的，利用重力方向項考慮了兩相之間的浮力效應。兩相均采用5%的湍流強度。空氣出口設為開口邊界，整個物理過程在常溫常壓下進行，管內(nèi)視為光滑且壁面無滑移條件。界面摩擦阻力系數(shù)由修正后的AIAD模型決定，該模型通過ANSYS CFX控制語言編入程序。

2.2 與實驗現(xiàn)象對比

穩(wěn)壓器波動管結(jié)構(gòu)復雜，發(fā)生CCFL現(xiàn)象的位置也較多，主要有CCFL-L、CCFL-S、CCFL-U三個位置。CCFL-L和CCFL-U發(fā)生的幾何結(jié)構(gòu)分別與傾斜直管段和豎直管段相似，而CCFL-S現(xiàn)象是穩(wěn)壓器波動管特有的，這里重點關(guān)注CCFL-S現(xiàn)象。由于氣相速度的增大，會對液相的流動形成阻力，當氣相速度足夠大時，液相無法繼續(xù)按照原方向流動甚至向相反的方向流動，進而形成阻塞，沿著管道逆向推進，直至達到一個平衡的狀態(tài)。通過數(shù)值模擬的結(jié)果（圖2）可以明顯地看到，阻塞推進到波動管彎頭位置，并伴有水躍現(xiàn)象，這是典型的CCFL現(xiàn)象。這定性地說明修正過的AIAD模型可以模擬出CCFL現(xiàn)象，選取實驗中典型的工況和位置來驗證模擬結(jié)果的準確性，當θ=0.6°，JG，in=3.8 m·s-1，JL，in=0.07 m·s-1時，液相反向流動，發(fā)生CCFL現(xiàn)象，實驗現(xiàn)象如圖3（a）所示，模擬結(jié)果如圖3（b）所示。模擬結(jié)果與實驗現(xiàn)象吻合證明了AIAD模型對于穩(wěn)壓器波動管CCFL現(xiàn)象的相間作用模擬的準確性。

圖2 波動管波峰推進圖Fig.2 Wave peak propulsion diagram of surge line

3 計算結(jié)果分析

本文選取了穩(wěn)壓器波動管上10個關(guān)鍵位置節(jié)點（圖4）。這些位置都是彎管的入口和出口部分，因為流體經(jīng)過一個彎管時，流道角度變化大，通常在這樣的位置更容易發(fā)生CCFL現(xiàn)象。

圖3 傾斜管段實驗現(xiàn)象(a)和模擬結(jié)果(b)Fig.3 Experimental phenomena(a)and simulation results(b)of the inclined pipe

此外反應堆中流道的阻塞會引起流體流動速度的分布發(fā)生變化，流動速度分布不均勻會引起局部換熱能力下降和壓力分布惡化等，而流道壓力分布的變化會在局部造成壓力過大，影響管道的安全性。通過對10個關(guān)鍵位置的抽取數(shù)據(jù)來定量分析氣相流速、波動管傾斜角對阻塞形成和推進的影響。

圖4 節(jié)點位置Fig.4 Node location diagram

本次模擬計算中，對不同初始氣體速度條件下阻塞形成位置和推進情況的數(shù)據(jù)進行對比，如圖5所示。圖5中L代表所選位置與管道起始點的距離，管道全長1.9 m，從右至左5條曲線分別代表5種從低到高的初始氣相流速。曲線上由低到高不同的點代表10個關(guān)鍵位置。通過對比可以得到：

1）從整體來看，初始氣相流速越大，阻塞推進到各個關(guān)鍵位置的時間越短，速度越快。對于靠近管道起始點的位置（位置1、2、3），當初始氣相流速足夠小的時候，阻塞在相鄰兩位置點的推進速度并無明顯差異；但是在初始氣相流速增大到一定程度時（如方塊點構(gòu)成的折線），位置1、2處將不會出現(xiàn)阻塞，這說明在位置1、2處氣相流速占主導地位。

2）對于遠離管道起始點的位置4～10，只考慮氣相流速的增加的影響，應該出現(xiàn)位置8～9阻塞推進比位置6～7更快的現(xiàn)象。但是由于傾斜角導致的位置8～9高度的變化，出現(xiàn)了位置6～7阻塞推進更快的現(xiàn)象，這說明在位置4～10高度變化占主導地位。

圖5 不同速度對阻塞形成和推進的影響Fig.5 Effects of different velocities on slug formation and propulsion

阻塞推進受初始氣相流速和傾斜角的復合影響，在氣體初始流速相同的情況下，本文分別選取傾斜角θ=0.6°、θ=2.5°、θ=5°時的數(shù)據(jù)進行對比分析（圖6），得到：

1）整體來看，在氣相流速不變的情況下，波動管傾斜角越大，阻塞形成時間越晚，推進速度越慢，相鄰節(jié)點之間推進時間間隔長。對于位置1～3，θ=0.6°、θ=2.5°、θ=5°三種情況下阻塞形成時間隨傾斜角度增加而推遲，但是推進速度大致相同（斜率大小基本一致），這說明傾斜度對靠近管道起始點的位置，即高度變化小的位置的阻塞形成時間有影響，但對推進速度影響不大。

2）對于位置4～10，在氣相流速不變的情況下，隨著傾斜角度的增加，阻塞推進速度明顯減慢，θ=5°時阻塞甚至推進不到位置10就達到了平衡。這是由于隨著管長增加，傾斜度越大引起的高度變化而導致的。

圖6 不同傾斜角度對阻塞形成和推進的影響Fig.6 Effects of different inclined angle on slug formation and propulsion

4 結(jié)語

本文使用ANSYS CFX軟件對穩(wěn)壓器波動管中CCFL現(xiàn)象進行三維數(shù)值模擬，并以經(jīng)典實驗數(shù)據(jù)為參照標準，來探究CCFL現(xiàn)象的特性。為了更好解釋CCFL現(xiàn)象，模擬過程使用自由表面模型，并且加入修正后的AIAD模型。通過對計算結(jié)果數(shù)據(jù)的處理，得到了不同初始氣相流速、不同穩(wěn)壓器波動管傾斜角度下阻塞推進的規(guī)律：

1）傾斜角度一定的情況下，初始氣相流速越大，阻塞向前推進的速度越快。

2）在氣相流速不變的情況下，波動管傾斜角越大，阻塞形成時間越晚，推進速度越慢，相鄰節(jié)點之間推進時間間隔長。

3）在靠近管道起始點的位置，即高度變化小的位置，阻塞推進主要受初始氣相流速的影響；在遠離管道起始點位置處，阻塞推進主要受傾斜角度的影響，這是由于高度變化大而導致的。