韓朋飛，郭烈錦*，李清平，姚海元，程 兵

(1. 西安交通大學(xué)動(dòng)力工程多相流國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049；2. 中海油研究總院，北京 100028)

雙流體模型中人工擴(kuò)散的影響數(shù)值分析

韓朋飛1，郭烈錦1*，李清平2，姚海元2，程 兵2

(1. 西安交通大學(xué)動(dòng)力工程多相流國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049；2. 中海油研究總院，北京 100028)

研究了雙流體模型中加入人工擴(kuò)散項(xiàng)后不同擴(kuò)散項(xiàng)系數(shù)對模擬的影響。對于豎直管的Water Faucet問題，逐步增加人工擴(kuò)散系數(shù)后非物理振蕩完全消除。對于水平管流動(dòng)病態(tài)區(qū)域的工況，不采用擴(kuò)散系數(shù)在中等精細(xì)度網(wǎng)格上計(jì)算的結(jié)果同實(shí)驗(yàn)比較吻合，但無法得到網(wǎng)格無關(guān)解。采用Water Faucet問題確定的人工擴(kuò)散系數(shù)進(jìn)行計(jì)算，水動(dòng)力學(xué)不穩(wěn)定性被抑制，液塞頻率被低估。因此，人工擴(kuò)散的量隨不同流動(dòng)條件發(fā)生變化，需要具體確定。

雙流體模型；人工擴(kuò)散；段塞流；數(shù)值模擬

0 引 言

管道內(nèi)多相流流動(dòng)模擬對石油工業(yè)的正常生產(chǎn)及安全保障有非常重要的意義。目前，石油開采正在向海洋尤其是深海進(jìn)發(fā)，因此超長距離的管道多相混輸對管內(nèi)流動(dòng)的預(yù)測有著更高和更迫切的需求。

根據(jù)以往國內(nèi)外的研究，一維的漂移流模型和雙流體模型是模擬預(yù)測超長距離管道流動(dòng)的有力工具。漂移流模型本質(zhì)上屬于混合物模型，無法對相界面上波動(dòng)的發(fā)展進(jìn)行預(yù)測，因此精確度有限。Lin等[1]很早就從理論上證明了雙流體模型可以預(yù)測界面上波動(dòng)的發(fā)展，并且Issa等[2]也使用雙流體模型在數(shù)值模擬中捕捉到了界面波動(dòng)的實(shí)時(shí)發(fā)展。因此，雙流體模型是對長距離管線流動(dòng)進(jìn)行準(zhǔn)確預(yù)測的最好手段。但是，傳統(tǒng)的標(biāo)準(zhǔn)雙流體模型在某些區(qū)域不能保證雙曲性，界面上小的波動(dòng)會(huì)被無限放大。這種特征表現(xiàn)在數(shù)值模擬當(dāng)中，就是當(dāng)網(wǎng)格步長不斷減小時(shí)，不斷有更小的波動(dòng)被計(jì)算分辨，而計(jì)算對這些波動(dòng)的放大系數(shù)不是有界的，在波長極小時(shí)趨于無窮大，從而無法得到唯一正確的收斂結(jié)果。

針對這一問題，很多學(xué)者采用各種方法進(jìn)行了研究。一方面，Zanotti等[3]從純數(shù)學(xué)的角度上考慮，通過矩陣預(yù)處理實(shí)現(xiàn)了無條件的雙曲性，盡管對算例的模擬顯示瞬態(tài)結(jié)果也是準(zhǔn)確的，但是理論上預(yù)處理影響了模型的瞬態(tài)項(xiàng)，只適用于穩(wěn)態(tài)問題。另一方面，很多學(xué)者嘗試引入其他物理機(jī)制來保證雙曲性，得到了一些有用的結(jié)論。Chung等[4]發(fā)現(xiàn)虛擬質(zhì)量力的加入并不是完全有益的，有時(shí)它會(huì)影響計(jì)算的穩(wěn)定性及精確度。Holmas等[5]則指出表面張力系數(shù)作用范圍有限，人為增大表面張力系數(shù)雖然可以抑制小波長不穩(wěn)定性的發(fā)展，但是這部分非物理的能量會(huì)轉(zhuǎn)移到長波里，而通過人工增加二階擴(kuò)散項(xiàng)耗散掉不穩(wěn)定的小波長波動(dòng)的能量是一種合理的方法。Song等[6]通過考慮截面速度不均勻引入動(dòng)量通量參數(shù)，也在廣泛的區(qū)域內(nèi)得到了良態(tài)的系統(tǒng)。但是，Issa等[7]發(fā)現(xiàn)動(dòng)量通量參數(shù)的預(yù)測結(jié)果同實(shí)驗(yàn)偏離。他們經(jīng)過同實(shí)驗(yàn)的對比發(fā)現(xiàn)Holmas的人工擴(kuò)散方法能夠在網(wǎng)格不斷加密的情況下得到同實(shí)驗(yàn)吻合的預(yù)測結(jié)果[8]。但是，他們只對水平管的一個(gè)工況進(jìn)行了研究，提出的系數(shù)對豎直管計(jì)算的作用還不清楚，需要進(jìn)一步研究。

本文針對目前保持雙流體模型雙曲性的研究，對雙流體模型模擬豎直管內(nèi)的流動(dòng)進(jìn)行探索，尋求獲得正確預(yù)測結(jié)果所需人工擴(kuò)散項(xiàng)系數(shù)，并驗(yàn)證此組系數(shù)計(jì)算水平管段塞流時(shí)的精確度。

1 控制方程

加入人工擴(kuò)散項(xiàng)后的一維雙流體模型如下：

(1)

(2)

(3)

(4)

αl+αg=1，

(5)

式中：下標(biāo)l和g分別表示液相和氣相；下標(biāo)w和i分別表示在壁面處和相界面處；α表示相含率；ρ表示密度；u表示速度；p表示壓力；g表示重力加速度；hl表示液膜厚度；τ表示切應(yīng)力；θ表示管道傾角；S表示濕周；A表示管道橫截面積。氣相為理想可壓縮氣體，液相為不可壓縮流體。

式(1)～(4)左邊第三項(xiàng)二階偏導(dǎo)即為加入的人工擴(kuò)散項(xiàng)，φal和φag分別表示液相和氣相連續(xù)方程的人工擴(kuò)散項(xiàng)系數(shù)，μal和μag分別表示液相和氣相動(dòng)量方程的人工擴(kuò)散項(xiàng)系數(shù)。剪切應(yīng)力的表達(dá)式如下：

(6)

式中：f為Fanning摩擦系數(shù)；兩相的相對速度ur=ug-ul。

對于水平管和豎直管，兩相的濕周計(jì)算分布如圖1所示。

圖1 濕周及兩相分布示意圖Fig.1 Schematic representation of wetted perimeters and two-phase distributions

水平管相界面上的摩擦系數(shù)采用Taitel-Dukler模型[9]：

(7)

豎直管相界面上的摩擦系數(shù)采用Wallis[10]的模型：

(8)

D為管道直徑。水平管氣相與壁面間的摩擦系數(shù)同樣采用Taitel-Dukler的模型[9]：

(9)

水平管液相與壁面間的摩擦系數(shù)采用Spedding-Hand模型[11]：

(10)

豎直管液相與壁面間的摩擦系數(shù)采用Wallis模型[10]：

(11)

雷諾數(shù)的表達(dá)式為

(12)

2 數(shù)值方法

氣液相的連續(xù)和動(dòng)量方程式(1)～(4)通過有限體積法在同位網(wǎng)格上進(jìn)行求解，采用Rhie-Chow插值以消除壓力速度的失耦。壓力修正方程采用兩相混合的質(zhì)量守恒方程導(dǎo)出，其形式如下：

(13)

式中：D/Dt為所給量的物質(zhì)導(dǎo)數(shù)。

為了提高計(jì)算的穩(wěn)定性，時(shí)間離散格式為一階歐拉隱式后向差分，空間離散格式采用二階迎風(fēng)。為了保證計(jì)算收斂性，壓力與速度的耦合采用SIMPLE算法，并設(shè)定全局最大Courant數(shù)為0.1，根據(jù)Courant數(shù)自動(dòng)選取時(shí)間步長，計(jì)算收斂的判據(jù)為1×10-6。

3 結(jié)果及討論

3.1 Water Faucet問題

Water Faucet問題是檢驗(yàn)數(shù)值計(jì)算結(jié)果精確度的常用標(biāo)準(zhǔn)方法，其過程如圖2所示。其結(jié)構(gòu)為一段12 m長、內(nèi)徑為1 m的豎直管道，上部入口含氣率為0.2，氣相和液相的速度分別為0和10 m/s，出口壓強(qiáng)固定為105Pa。氣相的密度為1.16 kg/m3，液相的密度為1 000 kg/m3，整個(gè)系統(tǒng)的溫度固定為50 ℃。此問題的解析解為

(14)

.

(15)

當(dāng)采取不同的擴(kuò)散項(xiàng)系數(shù)和網(wǎng)格控制容積(cv)數(shù)量時(shí)，模擬得到的t=0.5 s和t=2.5 s時(shí)的結(jié)果和解析解的結(jié)果對比分別如圖3和圖4所示。從圖3可以看出，當(dāng)未引入人工擴(kuò)散時(shí)，雙流體模型是不穩(wěn)定的病態(tài)問題，網(wǎng)格密度為500 cv時(shí)，含氣率曲線在梯度劇變處發(fā)生了大幅振蕩。而如圖4顯示，未引入人工擴(kuò)散的計(jì)算并不能得到穩(wěn)定的穩(wěn)態(tài)結(jié)果。

圖3 不同擴(kuò)散項(xiàng)系數(shù)及網(wǎng)格密度在t=0.5 s時(shí)模擬與解析解比較Fig.3 Comparison of simulation and analytical results for different diffusion coefficients and mesh densities at t=0.5 s

圖4 不同擴(kuò)散項(xiàng)系數(shù)及網(wǎng)格密度在t=2.5 s時(shí)模擬與解析解比較Fig.4 Comparison of simulation and analytical results for different diffusion coefficients and mesh densities at t=2.5 s

當(dāng)擴(kuò)散系數(shù)采用Issa等[8]的系數(shù)時(shí)，含氣率的非物理振蕩已被減弱但尚未消失。若進(jìn)一步增大動(dòng)量擴(kuò)散系數(shù)，計(jì)算沒有明顯的改善，因此逐步增大質(zhì)量擴(kuò)散系數(shù)。質(zhì)量擴(kuò)散系數(shù)從0.1開始逐步增大時(shí)的計(jì)算結(jié)果與解析解的比較如圖5所示。從圖中可從看到隨著質(zhì)量擴(kuò)散系數(shù)的逐步增大，非物理波動(dòng)的幅度逐漸減?。划?dāng)質(zhì)量擴(kuò)散系數(shù)達(dá)到0.7時(shí)，計(jì)算結(jié)果已經(jīng)不存在非物理振蕩。雖然由于擴(kuò)散的影響，含氣率的曲線變得更加平滑，但是能夠得到符合物理過程的無振蕩解，損失的精度是可以接受的。當(dāng)網(wǎng)格密度進(jìn)一步加密到1 500 cv和3 000 cv時(shí)，仍舊沒有非物理振蕩出現(xiàn)，可以認(rèn)為此時(shí)能夠得到所求解問題的網(wǎng)格無關(guān)解。

圖5 不同質(zhì)量擴(kuò)散項(xiàng)系數(shù)在t=0.75 s時(shí)模擬與解析解比較Fig.5 Comparison of simulation and analytical results for different diffusion coefficients at t=0.75 s

3.2 水平管段塞流

Woods等[12]列舉的內(nèi)徑為0.095 m水平管內(nèi)各種表觀氣速uSG和表觀液速下uSL的段塞流液塞頻率如圖6所示。本文選取了氣液相入口表觀速度分別為15.9 m/s和1.2 m/s遠(yuǎn)離穩(wěn)定邊界的病態(tài)區(qū)域工況，分別使用較粗網(wǎng)格不加入人工擴(kuò)散和人工擴(kuò)散系數(shù)φ=0.7，μ=1×10-5進(jìn)行計(jì)算。

圖6 Woods等[2]列舉的液塞頻率實(shí)驗(yàn)數(shù)據(jù)Fig.6 Experimental data for slug frequency presented by Woods et al[12]

圖7 模擬得到的沿管道無量綱液位高度隨時(shí)間的變化(φ=μ=0)Fig.7 Time trace of simulated non-dimensional liquid film height along the pipe (φ=μ=0)

在網(wǎng)格密度為1 500 cv時(shí)，計(jì)算得到的無量綱液位高度(hl/D)的瞬態(tài)曲線隨時(shí)間的變化如圖7所示，圖中相鄰曲線間的時(shí)間間隔為0.4 s。計(jì)算初始狀態(tài)管道內(nèi)無量綱液位高度為均勻的0.5，隨著氣液間相互作用的發(fā)展，約2 s時(shí)相界面的波動(dòng)開始出現(xiàn)并在約3 s時(shí)發(fā)展為液塞。由于初始條件設(shè)置，液塞前面的液膜厚度較大，液塞長度能夠不斷增長，因此需等第一個(gè)長液塞流出后才能反映真實(shí)的流動(dòng)。從后續(xù)的發(fā)展曲線來看，并未出現(xiàn)此種超長的液塞。在第一個(gè)超長液塞流出管道后，統(tǒng)計(jì)得到的液塞頻率為0.88 Hz，與0.79 Hz的實(shí)驗(yàn)結(jié)果相比，誤差為11.4%。雖然計(jì)算結(jié)果同實(shí)驗(yàn)相比很接近，但是當(dāng)加密網(wǎng)格后，進(jìn)一步被解析的短波波動(dòng)使得計(jì)算很快發(fā)散，無法得到網(wǎng)格獨(dú)立解。

采用前文驗(yàn)證的擴(kuò)散系數(shù)在網(wǎng)格密度為3 000 cv的細(xì)網(wǎng)格上進(jìn)行模擬，得到的無量綱液位高度的瞬態(tài)曲線隨時(shí)間的變化如圖8所示。從圖8可以看出，液面高度的波動(dòng)被抑制，曲線變得平滑，液塞形成的頻率降低，統(tǒng)計(jì)得到的液塞頻率降低為0.58 Hz，與實(shí)驗(yàn)相比減小了26.5%，顯然此時(shí)人工擴(kuò)散量是過大的。但是，進(jìn)一步加密網(wǎng)格未能顯著改善預(yù)測結(jié)果，證明人工擴(kuò)散量在精細(xì)網(wǎng)格上仍舊過大，得到的收斂結(jié)果對液塞頻率的預(yù)測偏低?？梢灶A(yù)測當(dāng)擴(kuò)散量很大時(shí)，雙流體模型的性能類似于混合物模型。

圖8 模擬得到的沿管道無量綱液位高度隨時(shí)間的變化(φ=0.7,μ=1×10-5)Fig.8 Time trace of simulated non-dimensional liquid film height along the pipe(φ=0.7,μ=1×10-5)

4 結(jié) 語

通過加入人工擴(kuò)散可以使雙流體模型在原本病態(tài)的區(qū)域內(nèi)轉(zhuǎn)變?yōu)榱紤B(tài)問題。經(jīng)過一系列的數(shù)值實(shí)驗(yàn)，得出在人工擴(kuò)散項(xiàng)系數(shù)為φ=0.7，μ=1×10-5時(shí)，豎直管Water Faucet問題非物理振蕩完全消失。對于水平管內(nèi)流動(dòng)病態(tài)區(qū)域的工況，φ=0.7，μ=1×10-5時(shí)水動(dòng)力學(xué)不穩(wěn)定性被抑制，液塞頻率降低。水平管人工擴(kuò)散項(xiàng)的系數(shù)需要根據(jù)不同工況決定，有待進(jìn)一步研究。

[1] Lin P Y, Hanratty T J. Prediction of the initiation of slugs with linear stability theory[J]. International Journal of Multiphase Flow, 1986, 12(1): 79.

[2] Issa R I, Kempf M H W. Simulation of slug flow in horizontal and nearly horizontal pipes with the two-fluid model[J]. International Journal of Multiphase Flow, 2003, 29(1): 69.

[3] Zanotti A L, Méndez C G, Nigro N M, et al. A preconditioning mass matrix to avoid the ill-posed two-fluid model[J]. Journal of Applied Mechanics, 2007, 74(4): 732.

[4] Chung M-S, Lee S-J. A modified semi-implicit method for a hyperbolic two-fluid model[J]. Appl Numer Math, 2009, 59(10): 2475.

[5] Holmas H, Sira T, Nordsveen M, et al. Analysis of a 1D incompressible two-fluid model including artificial diffusion[J]. IMA Journal of Applied Mathematics, 2008, 73(4): 651.

[6] Song J H, Ishii M. The well-posedness of incompressible one-dimensional two-fluid model[J]. International Journal of Heat and Mass Transfer, 2000, 43(12): 2221.

[7] Issa R I, Montini M. Applicability of the momentum-flux-parameter closure for the two-fluid model to slug flow[C]. 6th International Symposium on Multiphase Flow, Heat Mass Transfer and Energy Conversion, 2009: 712.

[8] Issa R I, Montini M. The effect of surface tension and diffusion on one-dimensional modelling of slug flow instabilities[C]. 7th International Conference on Multiphase Flow, 2010.

[9] Taitel Y, Dukler A E. A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow[J]. AIChE Journal, 1976, 22(1): 47.

[10] Wallis G B. One-Dimensional Two-Phase Flow[M]. New York: McGraw-Hill,1969.

[11] Spedding P L, Hand N P. Prediction in stratified gas-liquid co-current flow in horizontal pipelines[J]. International Journal of Heat and Mass Transfer, 1997, 40(8): 1923.

[12] Woods B, Fan Z, Hanratty T. Frequency and development of slugs in a horizontal pipe at large liquid flows[J]. International Journal of Multiphase Flow, 2006, 32(8): 902.

NumericalAnalysisoftheEffectofArtificialDiffusionTermsonTwo-FluidModel

HAN Peng-fei1, GUO Lie-jin1, LI Qing-ping2, YAO Hai-yuan2, CHENG Bing2

(1.StateKeyLaboratoryofMultiphaseFlowinPowerEngineering,Xi’anJiaotongUniversity,Xi’an,Shaanxi710049,China; 2.CNOOCResearchInstitute,Beijing100028,China)

When artificial diffusion terms are introduced in two-fluid model, the effects of different artificial diffusion coefficients are studied. For vertical pipe Water Faucet problem, when the artificial diffusion coefficients are gradually increased, nonphysical oscillations can be eliminated. For horizontal pipe flow in ill-posed region， the calculated results without diffusion terms agree with experimental data but the mesh independent solution cannot be obtained. When the coefficients calibrated in Water Faucet problem are adopted, hydrodynamic instability is suppressed to some extent and slug frequency is underestimated. Therefore, the amount of artificial diffusion varies with different flow conditions and should be determined specially.

two-fluid model; artificial diffusion; slug flow; numerical simulation

O359

A

2095-7297(2015)01-0023-05

2014-12-29

國家科技重大專項(xiàng)(2011ZX05026-004-02)、國家自然科學(xué)基金重點(diǎn)項(xiàng)目(51236007)

韓朋飛(1983—),男，博士研究生，主要從事氣液兩相流方面的研究。

*通信作者