基于修正AIAD模型的穩(wěn)壓器波動管CCFL現(xiàn)象數(shù)值模擬

趙柏陽 李 冬 楊已顥 張永祺

（上海電力大學(xué) 能源與機(jī)械工程學(xué)院 上海 200090）

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

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

1 計算模型

1.1 自由表面模型

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

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

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

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

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

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

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

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

1.2 修正的AIAD模型

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

式中：rG為氣相體積分?jǐn)?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模型的基礎(chǔ)上對混合函數(shù)f的表達(dá)式進(jìn)行了修正，目的是為了使混合多種影響因素后總的阻力系數(shù)隨著局部流動特征引起的不同參數(shù)而變化，這樣可以避免因為經(jīng)驗值選取不當(dāng)而導(dǎo)致的計算偏差。為便于區(qū)分，函數(shù)用φ來表示。

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

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

表1 體積分?jǐn)?shù)和相應(yīng)的形態(tài)Table 1 Volume fractions and corresponding morphologies

2 數(shù)值模型建立

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

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

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

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

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

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

穩(wěn)壓器波動管結(jié)構(gòu)復(fù)雜，發(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)象。由于氣相速度的增大，會對液相的流動形成阻力，當(dāng)氣相速度足夠大時，液相無法繼續(xù)按照原方向流動甚至向相反的方向流動，進(jìn)而形成阻塞，沿著管道逆向推進(jìn)，直至達(dá)到一個平衡的狀態(tài)。通過數(shù)值模擬的結(jié)果（圖2）可以明顯地看到，阻塞推進(jìn)到波動管彎頭位置，并伴有水躍現(xiàn)象，這是典型的CCFL現(xiàn)象。這定性地說明修正過的AIAD模型可以模擬出CCFL現(xiàn)象，選取實驗中典型的工況和位置來驗證模擬結(jié)果的準(zhǔn)確性，當(dāng)θ=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)象的相間作用模擬的準(zhǔn)確性。

圖2 波動管波峰推進(jìn)圖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

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

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

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

1）從整體來看，初始?xì)庀嗔魉僭酱?，阻塞推進(jìn)到各個關(guān)鍵位置的時間越短，速度越快。對于靠近管道起始點的位置（位置1、2、3），當(dāng)初始?xì)庀嗔魉僮銐蛐〉臅r候，阻塞在相鄰兩位置點的推進(jìn)速度并無明顯差異；但是在初始?xì)庀嗔魉僭龃蟮揭欢ǔ潭葧r（如方塊點構(gòu)成的折線），位置1、2處將不會出現(xiàn)阻塞，這說明在位置1、2處氣相流速占主導(dǎo)地位。

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

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

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

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

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

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

4 結(jié)語

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

1）傾斜角度一定的情況下，初始?xì)庀嗔魉僭酱?，阻塞向前推進(jìn)的速度越快。

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

3）在靠近管道起始點的位置，即高度變化小的位置，阻塞推進(jìn)主要受初始?xì)庀嗔魉俚挠绊?；在遠(yuǎn)離管道起始點位置處，阻塞推進(jìn)主要受傾斜角度的影響，這是由于高度變化大而導(dǎo)致的。