王冬旭，胡其會(huì)，李玉星，李爽，王權(quán)

（1 中國(guó)石油大學(xué)（華東），山東省油氣儲(chǔ)運(yùn)安全省級(jí)重點(diǎn)實(shí)驗(yàn)室，山東青島266580；2 中國(guó)市政工程中南設(shè)計(jì)研究總院有限公司，湖北武漢430060）

段塞流是油氣混輸管線中的常見(jiàn)流型，當(dāng)管線處于啟動(dòng)、停輸、清管、輸量變化甚至正常生產(chǎn)時(shí)，都有可能產(chǎn)生[1]。其流動(dòng)的間歇性將引起管道中持液率和壓力的急劇波動(dòng)，并使得下游管道及設(shè)備承受間歇性應(yīng)力沖擊[2-4]。同時(shí)，管道末端產(chǎn)生的大段塞會(huì)引起下游處理設(shè)備中液位劇烈波動(dòng)。為保證管線和下游處理設(shè)備的安全生產(chǎn)，必須掌握段塞流的特征規(guī)律。為此，自20 世紀(jì)40 年代起，國(guó)內(nèi)外學(xué)者對(duì)段塞流開(kāi)展了相關(guān)研究[5-7]。

Dukler 等[8]最早建立了氣泡和段塞在時(shí)間和空間內(nèi)呈均勻分布的段塞流穩(wěn)態(tài)模型。此模型可以預(yù)測(cè)段塞流重要參數(shù)的平均值，如段塞長(zhǎng)度、氣泡平移速度和壓降。Ishii[9]在20世紀(jì)70年代建立的歐拉雙流體模型廣泛應(yīng)用于核工業(yè)。Bendiksen 等[10]建立的OLGA模型在石油工業(yè)中被廣泛應(yīng)用，是最早以雙流體模型為基礎(chǔ)開(kāi)發(fā)并使用的瞬態(tài)計(jì)算代碼之一。Zheng 等[11]通過(guò)研究簡(jiǎn)單起伏管路內(nèi)的段塞流建立了一個(gè)SINK&SOURCE 模型。 Nydal 和Banerjee[12]，Taitel 和Barnea[13]也建立了段塞模型，以上模型都需要額外的段塞流參數(shù)，如段塞形成條件、段塞頻率、段塞長(zhǎng)度等。

應(yīng)用于雙流體模型的線性穩(wěn)定性理論經(jīng)常被用來(lái)確定分層流的穩(wěn)定性。Taitel 和Dukler[14]將該理論應(yīng)用于不考慮剪切應(yīng)力的理想流體，得到非黏性Kelvin-Helmholtz（IKH）穩(wěn)定性準(zhǔn)則。該準(zhǔn)則可以預(yù)測(cè)伯努利效應(yīng)產(chǎn)生的吸力大于流體重力時(shí)氣液界面波的生長(zhǎng)。Barnea 和Taitel[15]隨后通過(guò)對(duì)黏性流場(chǎng)的線性穩(wěn)定性分析，得出黏性Kelvin-Helmholtz（VKH）穩(wěn)定性準(zhǔn)則。Issa和Kempf[16]論證了當(dāng)氣液流速處于IKH 和VKH 準(zhǔn)則之間時(shí)，使用交錯(cuò)網(wǎng)格的雙流體模型能夠捕獲分層流過(guò)渡到段塞流過(guò)程中氣液界面產(chǎn)生的擾動(dòng)。氣液兩相流在初始條件下自由發(fā)展，減少了對(duì)物理現(xiàn)象模型（流型轉(zhuǎn)變、段塞形成等）的依賴。Issa 等[16]利用此模型模擬了V 型兩相流起伏管道內(nèi)的段塞頻率、段塞平均長(zhǎng)度和持液率，得到與實(shí)驗(yàn)數(shù)據(jù)一致的波動(dòng)趨勢(shì)。

Renault[17]以Issa 的研究為基礎(chǔ)，基于IKH 和VKH準(zhǔn)則實(shí)現(xiàn)了分層流過(guò)渡到段塞流的模擬，捕捉段塞前后界面運(yùn)動(dòng)過(guò)程，取得了較好的預(yù)測(cè)效果。盡管Renault 模型基于IKH和VKH 準(zhǔn)則捕捉段塞前后液面變化，然而其求解過(guò)程中將液相方程轉(zhuǎn)化為淺水方程，并在單元格之間采用黎曼（Riemann）精確解，使得求解速度較慢。因此，需對(duì)模型進(jìn)一步改進(jìn)，使其保持精度的同時(shí)提高計(jì)算效率。

許仁義等[18-20]將行波法與Riemann 精確解相結(jié)合，對(duì)稀疏波采用單波近似，實(shí)現(xiàn)了淺水方程的快速求解，實(shí)現(xiàn)了高分辨率、高計(jì)算精度的效果。因此本文嘗試將行波法引入Renault 模型，建立新的段塞捕捉模型及求解方法。同時(shí)為保證行波法適用于所有計(jì)算單元，對(duì)可能出現(xiàn)的干區(qū)采用薄液膜進(jìn)行簡(jiǎn)化處理，在保證模型精度的條件下進(jìn)一步提高求解速度。

1 段塞捕捉模型建立

1.1 模型基本方程

Renault 建立的段塞捕捉模型假設(shè)條件為：①液相不可壓縮；②理想氣相；③相對(duì)于液相動(dòng)量，氣相動(dòng)量可以忽略；④求解液相動(dòng)量方程時(shí)，氣相為局部不可壓縮。計(jì)算域內(nèi)網(wǎng)格被劃分成段塞網(wǎng)格和氣泡網(wǎng)格兩類。段塞網(wǎng)格存儲(chǔ)平均液相速度（Uslug）、氣泡網(wǎng)格存儲(chǔ)含液率（β）、平均液速（Ul）、氣體表觀流速（USg）、壓力（p）存儲(chǔ)在交錯(cuò)網(wǎng)格上，計(jì)算網(wǎng)格如圖1 所示，模型方程如式(1)～式(5)。

1.2 閉合關(guān)系式

圖1 計(jì)算網(wǎng)格

液相和壁面、氣相和壁面以及氣液界面間的剪切應(yīng)力關(guān)系如式(6)～式(8)。液相與壁面的摩阻系數(shù)、氣相于壁面的摩阻系數(shù)如式(9)和式(10)，氣液相界面摩擦因子采用由Cohen 與Hanratty[21]提出的常數(shù)值如式(11)，氣泡平移速度采用Bendiksen[22]關(guān)系式，如式(12)。

紊流情況下，根據(jù)弗魯?shù)聰?shù)（Fr）計(jì)算參數(shù)C0、Ut0。

1.3 模型求解

Renault 模型求解過(guò)程中氣相和液相方程的求解是交替進(jìn)行的。氣相方程采用一階迎風(fēng)格式化簡(jiǎn)為式(15)。

由于氣相方程符合三對(duì)角矩陣形式，因此采用追趕法即可求解每個(gè)網(wǎng)格中n+1 時(shí)刻的氣相狀態(tài)。

液相方程可化簡(jiǎn)為淺水方程形式，如式(16)和式(17)。

液相求解過(guò)程較為復(fù)雜，采用戈杜諾夫（Godunov）方法[22]，在每個(gè)單元邊界求解淺水方程Riemann問(wèn)題的精確解，進(jìn)而構(gòu)成整個(gè)流場(chǎng)的數(shù)值解[式(18)和式(19)]，求解過(guò)程如圖2（圖2中xj與xj+1左右兩側(cè)均為稀疏波）。

圖2 戈杜諾夫方法求解過(guò)程

1.4 模型簡(jiǎn)化

由于Renault 模型中液相求解需要對(duì)網(wǎng)格邊界產(chǎn)生的稀疏波區(qū)域持液率（βLM、βRM）進(jìn)行核算，計(jì)算量較大。為此，本文提出采用行波法求解單元邊界Riemann問(wèn)題，對(duì)稀疏波采用單波近似，同時(shí)用薄液膜代替干區(qū)，使行波法適用于干區(qū)計(jì)算，從而簡(jiǎn)化求解過(guò)程，提高計(jì)算效率。

圖3 Riemann問(wèn)題精確解結(jié)構(gòu)

圖4 行波法求解過(guò)程

對(duì)比圖2、圖4 求解過(guò)程發(fā)現(xiàn)，對(duì)稀疏波采用單波近似，并用行波法求解液相方程極大地簡(jiǎn)化了求解過(guò)程。淺水方程求解過(guò)程會(huì)出現(xiàn)干區(qū)情況，為使行波法適用于所有計(jì)算單元，使用薄液膜代替干區(qū)，薄液膜的液量來(lái)自相鄰單元。為了解行波法處理干濕邊界問(wèn)題時(shí)的性能，采用一種長(zhǎng)100m、寬1m、初速度為0、水面不連續(xù)的一維渠道淺水方程作為研究對(duì)象。用0.1m 的液膜代替干區(qū)，初始條件如式(22)。

計(jì)算范圍x∈[0,100m]被分成100 個(gè)網(wǎng)格，空間步長(zhǎng)dx=1m。潰壩初始時(shí)刻如圖5所示。初始時(shí)刻水面間斷處（x=50m）產(chǎn)生向左傳播的稀疏波和向右傳播的激波。圖6 為CFL=0.8、t=5s 時(shí)的計(jì)算結(jié)果。

從圖6可以看出，用薄液膜代替干區(qū)的簡(jiǎn)化方法得到的數(shù)值解與精確解相一致，可以反映液相在干濕邊界的流動(dòng)狀態(tài)。將Renault 模型單元邊界存在的干區(qū)均用薄液膜代替。此簡(jiǎn)化處理使得行波法適用于所有網(wǎng)格計(jì)算，極大地提高了運(yùn)算速度，簡(jiǎn)化后的全新計(jì)算流程如圖7所示。

圖5 初始液面狀態(tài)

圖6 潰壩模擬

圖7 計(jì)算流程

實(shí)驗(yàn)發(fā)現(xiàn)段塞運(yùn)動(dòng)過(guò)程中頭部呈豎直斷面、尾部呈曲面（圖8）。為使計(jì)算結(jié)果與實(shí)際段塞運(yùn)動(dòng)過(guò)程相一致，需對(duì)生成的段塞進(jìn)行特殊處理。段塞頭部速度Ufront根據(jù)液體質(zhì)量守恒核算，氣泡頭部速度UBendiksen由Bendiksen 關(guān)系式核算，計(jì)算過(guò)程如圖9。

圖8 段塞運(yùn)動(dòng)

圖9 段塞計(jì)算過(guò)程

2 模型驗(yàn)證

2.1 實(shí)驗(yàn)系統(tǒng)

本實(shí)驗(yàn)的介質(zhì)采用空氣和水，實(shí)驗(yàn)裝置示意圖如圖10 所示。實(shí)驗(yàn)管道內(nèi)徑D=40mm，總長(zhǎng)45m，其中水平測(cè)試段管長(zhǎng)10m。測(cè)試管段采用透明有機(jī)玻璃材質(zhì)，以便觀察管內(nèi)氣液流型變化。管道混合器出口處布置有壓力傳感器P0，水平段內(nèi)布置有間距均為0.8m 的雙平行電導(dǎo)探針（CP1～CP3）與壓力傳感器（P1～P3）。采用Iotech6220 采集卡進(jìn)行持液率與壓力信號(hào)的采集，采樣頻率為50Hz，采樣時(shí)間為60s。實(shí)驗(yàn)中的氣體表觀流速范圍為0～4m/s，液體表觀流速范圍為0～1.6m/s。

2.2 實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證

為驗(yàn)證模型的段塞捕捉能力，以實(shí)驗(yàn)管道為基準(zhǔn)進(jìn)行模擬計(jì)算。管道剖面尺寸如圖11，模型計(jì)算域?yàn)閇0,10]m，網(wǎng)格寬度Δx=0.02m，管徑D=40mm，時(shí)間步長(zhǎng)Δt=0.01s，模擬時(shí)間10min。氣體表觀流速范圍為0～4m/s，液體表觀流速范圍為0～1.6m/s。入口采用質(zhì)量流量邊界條件，出口采用壓力邊界條件，計(jì)算結(jié)果如下。

圖10 實(shí)驗(yàn)裝置示意圖

圖11 管道剖面示意圖

圖12 為氣體表觀流速0.9m/s、液體表觀流速0.45m/s、時(shí)間t=30s時(shí)的模擬結(jié)果。圖13為CP1處1min內(nèi)實(shí)驗(yàn)與計(jì)算持液率對(duì)比。從圖中可以看出，隨著管道長(zhǎng)度的增加，液面漸漸產(chǎn)生波動(dòng)，流型由分層流轉(zhuǎn)變?yōu)椴ɡ肆鳎浑S著波動(dòng)幅度的增大，流型進(jìn)一步發(fā)展成為段塞流；計(jì)算CP1 處1min 內(nèi)產(chǎn)生的段塞數(shù)略大于實(shí)測(cè)值，持液率變化與實(shí)測(cè)結(jié)果相一致，且段塞頭部呈豎直面，氣泡頭部呈曲面；本模型可以捕捉段塞流生成與發(fā)展過(guò)程，實(shí)現(xiàn)段塞流壓降精確計(jì)算。

圖12 t=30s持液率計(jì)算結(jié)果

圖14、圖15 分別為模擬實(shí)驗(yàn)工況30min 內(nèi)的平均壓降、平均段塞長(zhǎng)度值與實(shí)驗(yàn)值對(duì)比。圖16為兩者計(jì)算值與實(shí)測(cè)值相對(duì)誤差概率密度分布。從圖中可以看出，壓降、段塞長(zhǎng)度的相對(duì)誤差分別在25%、30%以內(nèi)，且主要分布在20%以內(nèi)，壓降計(jì)算準(zhǔn)確度要高于段塞長(zhǎng)度。

圖16 30min內(nèi)壓降與段塞長(zhǎng)度的相對(duì)誤差概率分布

為進(jìn)一步了解本文模型簡(jiǎn)化方法對(duì)計(jì)算效率的影響，對(duì)比Renault 模型與本文簡(jiǎn)化模型在50、100、200、300、400 網(wǎng)格數(shù)下的計(jì)算耗時(shí)。計(jì)算硬件條件CPU 為Intel i7-6700 3.4GHz、內(nèi)存為16GB，操作系統(tǒng)為64 位Windows 7，運(yùn)算軟件為Matlab 2017，對(duì)比結(jié)果如圖17所示。

由圖17 可以看出，兩種模型的計(jì)算耗時(shí)與網(wǎng)格數(shù)均成正相關(guān)，隨著網(wǎng)格數(shù)的增加，計(jì)算量增大，耗時(shí)逐漸增加；當(dāng)網(wǎng)格數(shù)較低時(shí)，兩種模型計(jì)算耗時(shí)差距較小，隨著網(wǎng)格數(shù)的增加，本文采用薄液膜代替干區(qū)與行波法求解液相方程組的簡(jiǎn)化方法優(yōu)勢(shì)逐漸明顯，運(yùn)算時(shí)間相較于原模型平均減少28%。

圖17 不同網(wǎng)格數(shù)下的計(jì)算耗時(shí)

3 結(jié)論

（1）本文將行波法引入到Renault模型求解中，對(duì)稀疏波采用單波近似，用Riemann近似解代替單元格之間的精確解，并用薄液膜代替干區(qū)，使行波法適用于所有計(jì)算單元，在保證計(jì)算精度的同時(shí)極大地提高了求解速度。

（2）將室內(nèi)小型環(huán)道實(shí)驗(yàn)數(shù)據(jù)與本模型計(jì)算結(jié)果對(duì)比，本模型可以模擬氣液兩相流由分層流發(fā)展到段塞流的物理過(guò)程，段塞頭部呈豎直面、氣泡頭部呈曲面，符合實(shí)驗(yàn)觀察結(jié)果，持液率計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果相一致，壓降、段塞長(zhǎng)度計(jì)算相對(duì)誤差分別在25%、30%以內(nèi)，且主要分布在20%以內(nèi)。

（3）對(duì)比本模型與Renault 模型在不同網(wǎng)格數(shù)下的計(jì)算效率，本模型運(yùn)算時(shí)間相比于Renault 模型平均減少28%，表明本文所提出的簡(jiǎn)化方法極大地提高了計(jì)算效率，該瞬態(tài)段塞流捕捉模型可用于油氣管道段塞流的快速準(zhǔn)確預(yù)測(cè)。

符號(hào)說(shuō)明

A——管路橫截面積，m2

Ag，Al——?dú)庀?、液相所占的流通面積，m2

D——管內(nèi)徑，m

Fr——Froude數(shù)

fg，fl，fi——?dú)獗?、液壁和氣液相間摩阻系數(shù)

g——重力加速度，m/s2

J——網(wǎng)格單元

LS——段塞長(zhǎng)度，m

n——表示時(shí)刻

p——管內(nèi)的壓力，Pa

Rel，Reg——液相雷諾數(shù)、氣相雷諾數(shù)

Sg，Sl，Si——?dú)獗?、液壁與氣液相間濕周，m

Δt——時(shí)間步長(zhǎng)，s

Ug，Ul——?dú)庀嗾鎸?shí)速度、液相真實(shí)速度，m/s

，——?dú)庀啾碛^速度和液相表觀速度，m/s

Δx——空間步長(zhǎng)，m

α——截面含氣率

β——截面含液率

θ——管道傾斜角度，(°)

?——液面夾角，(°)

μl，μg——?dú)庀嗪鸵合嗟膭?dòng)力黏度，Pa·s

ρl，ρg——?dú)庀嗪鸵合嗝芏?，kg/m3

σ——表面張力，N/m

τg，τl，τi——?dú)獗?、液壁和氣液相界面的剪切?yīng)力，Pa

下角標(biāo)

g——?dú)庀?/p>

i——?dú)庖合嚅g

l——液相

L——左側(cè)

M——中間

R——右側(cè)