氣井井筒條件下單液滴動力學特征及其攜帶臨界氣流速

王志彬 張亞飛 孫天禮 李忠城 楊中位 朱國 石紅艷

1. 西南石油大學油氣藏地質及開發工程國家重點實驗室；2. 中聯煤層氣有限責任公司；3. 中石化西南油氣分公司采氣二廠；4. 中石化西南油氣分公司基建部

0 引言

氣井連續攜液臨界氣流速是排水采氣工藝優化設計的重要參數［1］。攜液生產氣井井筒中流型為環狀流，一部分液體以液膜形式沿管壁流動，一部分液體以分散液滴形式被氣流夾帶。液膜及液滴的攜帶臨界氣流速直接影響環狀流場中液體連續攜帶的臨界氣流速。為此，氣井攜液臨界氣流速計算出現了液膜攜帶模型和液滴攜帶模型兩大學術派別。

液滴模型認為井內最大液滴回落是導致氣井積液的首要因素。1969年Turner等［2］率先建立了單液滴攜液模型，該模型假設液滴為圓球狀，液滴曳力系數和破碎韋伯數分別取為0.44和30，并認為若氣流能將環狀流場中的最大液滴帶出井口，則氣井不會積液。Turner單液滴模型因簡單，便于計算，四十年來被各大油田廣泛應用。然而在工程應用中發現，很多氣井攜液產氣量與Turner模型預測值相差較大，為此國內外學者對氣井連續攜液理論開展了大量的研究工作。Coleman等［3］、Nosseir等［4］、李閩等［5］、周德勝等［6］、王志彬等［7-12］分別考慮油壓、流型(層流或湍流)、液滴群聚效應、液滴變形、液滴翻滾、界面張力等因素的影響，先后建立了各種液滴攜帶新模型。目前幾乎國內各大氣田根據生產數據擬合了適合本氣田的攜液流量計算經驗模型，可見攜液規律之復雜。

氣井井筒高溫高壓條件下液滴的動力學特征，尤其是液滴的形狀特征、曳力系數以及與液滴尺寸有關的臨界韋伯數未見相關報道，這給準確分析液滴攜帶規律帶來了難度。這些參數可通過實驗測量，也可通過建模計算。氣井井筒高溫高壓下的實驗測試成本高、難度大，為此，筆者選用了數值模擬方法。

液滴形狀特征和曳力系數的獲取有3種數值模擬方法。第1種是將液滴置于大的計算區域中，讓氣流吹液滴直至液滴變形達到穩定［13-14］，或允許液滴在較高垂直密閉氣柱中自由下落［15］，該方法網格數量和計算量大。第2種是液滴界面追蹤法，該方法應用動網格捕捉界面特征［16］。第3種是施加體積力使液滴保持恒定的相對速度或固定在恒定位置［17-18］，該方法僅需要局部細化液滴及其周圍區域，可以大大減少計算量。同時，目前液滴動力學特征研究主要集中在雷諾數小于200，而對氣井井筒條件高雷諾數下的液滴動力學特征鮮有研究。

筆者建立了描述湍流場中液滴動力學特征的數值模型，開發了液滴動力學特征模擬求解器；計算了不同韋伯數下的液滴形狀特征、曳力系數及液滴破裂的臨界韋伯數；對液滴進行受力分析，得到液滴攜帶所需要的臨界氣流速預測新模型；利用文獻數據對現有單液滴模型進行了適應性評價，給出了單液滴模型適用的產液量條件。

1 單液滴動力學模型

1.1 數學模型的建立

1.1.1 基本方程

因計算區域較小，氣體的密度可考慮為常數。施加在流體上的力包括重力、表面張力和黏性應力。不可壓縮、黏性、不相溶流場的控制方程描述如下［19-21］

式中，U為笛卡爾坐標系中的速度矢量，m/s；p為壓力，Pa；t為時間，s； ?·為矢量發散算子；?為標量梯度算子；fS為界面處由表面張力引起的體積力，N；ρ為密度，kg/m3；μ為黏度，Pa · s；FB為以源項形式添加的體積力，N；g為重力加速度，9.8 m/s2。

氣液界面形狀是研究的重點，本研究使用流體體積函數法(VOF)研究氣液界面結構。對于相分布計算，在空間中采用體積函數α來識別網格中的組分。定義如下：對于充滿氣相的網格單元，αL=0；對于氣液界面處的單元，0＜αL＜1；對于充滿液相的單元，αL=1。

計算單元中的混合物理性質可以由體積分數的值計算

氣液表面張力會產生一個額外壓力梯度力，使用CSF模型［22］對每單位單元體積進行評估。

式中，σ為氣液界面張力，N/m；k為自由表面的平均曲率；下標L、G分別代表液相和氣相。

1.1.2 液滴曳力系數計算式

液滴在其自身重力和曳力的作用下加速或減速，作用于液滴的曳力可以通過液滴的計算網格積分得到。根據高斯定理，體積積分可以轉化為面積分。在本研究中，氣流沿y方向流動，即從底部向上部流動，由于液滴周圍氣流場呈對稱分布，故曳力垂直向上。所以y方向上的曳力可通過下式積分得到。

式中，Ω為液滴控制體；Γ為液滴控制體的表面；nx、ny、nz為垂直于液滴表面的單位矢量；u、v、w為x、y、z方向的速度分量。

由式(7)可知，總曳力系數由4部分組成或產生4種作用。一是使液滴加速，二是促使液滴內部循環(由于液滴內部存在速度梯度使得液滴變形和內部循環)，三是來自液滴迎風面和背風面之間的壓差，最后是克服液滴表面摩擦的黏滯阻力。后3項是通過對液滴表面單元進行微分來計算的。在式(7)中，nxΓcell、nyΓcell、nzΓcell(Γcell是表面單元面積大小)實際是液滴表面在x、y、z方向上的投影面積，可以通過相分數梯度和單元體積計算得到。

液滴前后壓差產生的曳力可以簡化為

由體積力計算液滴的總曳力系數為

式中，pjAj為施加在液滴迎風面的單元(j)上的壓力，N；piAi為施加在液滴背風面單元(i)上的壓力，N；iu、jd分別為液滴表面上方和下方單元，對于該單元，相分數α=0；Af為基于液滴初始直徑迎風面面積，m2；Af為基于液滴初始直徑迎風面面積，m2；D為液滴的初始直徑，m。

1.1.3 體積力計算式

為了避免液滴加速或減速對液滴曳力模擬準確性的影響，同時減小模擬計算工作量，本研究中對液滴施加了體積力以使其保持在恒定的位置。體積力在每個計算時間步中根據牛頓第二定律計算。

由于每個計算時間步長的液滴質心速度等于0，所以式(10)中vd,t-Δt=0，則式(10)可簡化為

yd,0是y方向上的初始質心，yd,t是當前時刻y方向的質心，由下式計算

式中，md為液滴質量，kg；ad,t為在體力FB,t作用下液滴當前時刻t的加速度，m/s2；vd,t，vd,t-Δt分別為當前時刻t和之前計算時刻t-Δt的液滴質心速度，m/s；Vi為單元(i)的體積，m3；yi為單元(i)在y方向上的質心，m。

事實上，FD未知，質心中心yd,t是體力FB,t的隱式函數。

牛頓迭代法被用來計算體積力，計算式為

yd,t是在上一體積力FB,t,old條件下液滴質心在y方向的位置，FB,t,new是當前體積力，grad(yd,t)是液滴質心位置yd,t對體積力FB,t,old的變化梯度。

在獲得體積力后，瞬時曳力可以由式(15)計算

當前時間步的瞬時曳力系數為

根據Prahl等(2006)［15］的研究，人工重力對液滴變形和破裂的影響可以忽略不計。

1.2 數學模型的求解

式(1)和式(2)使用有限體積法進行離散。瞬態項和源項使用中心差分進行離散并沿單元體積進行積分。應用高斯定理將擴散項和對流項差分轉化為單元體表面積分，單元面的值利用Van Leer通量限制求解。通過面積分來計算單元體的值，計算單元面的值通過插值計算。采用Jasak(1996)［23］和Waterson等(2007)［24-25］的求解方案就對流項進行離散處理。此外，時間微分項使用隱式Euler方法進行離散，該方案無條件穩定，一階時間準確，并保證解的有界性。相函數使用有限體積法進行離散化處理，并對整個單元體積和時間進行積分。面插值從相鄰計算單元的值計算得到。

紙條的疑問是在畢業四年后我結婚那天才得到解答的。那天我的許多老同學都來了，當然有鄧軍、胡波、趙小明這三個死黨。

整個求解過程基于PIMPLE算法總結如下：(1)在每個計算時間步驟之前保存所有初始場，包括速度u0、壓力p0、體積分數α0；(2) 求解相分數方程式(3)獲得更新的體積分數場；(3)通過式(4)和式(5)更新流體性質；(4) PIMPLE動量壓力耦合算法更新速度和壓力場；(5)用牛頓迭代程序計算體積力：①從u0、p0、α0等3個初始場參數恢復初始場；②使用PIMPLE算法計算新的速度和壓力場；③修正曳力；④計算液滴質心位置；⑤重復步驟①～④，直到前后2次計算的液滴質心位置滿足精度要求；(6)計算體積力；(7)計算曳力和曳力系數；(8)返回步驟(1)，并計算下一新的計算時間步的流場信息。基于以上求解過程，開發了單液滴動力學特征模擬器。

1.3 物理模型及網格劃分

在連續均勻速度場中放置一個液滴，計算域是一個邊長為64D的立方體，如圖1所示。將液滴置于計算域的中心。通過施加體積力將液滴保持在恒定位置，體積力由液滴的加速度計算得到。計算初始時刻，液滴初始形狀為球狀，直徑為D。入口為速度邊界，出口是壓力邊界。網格用七級局部加密，最大網格長度等于液滴直徑，最小網格長度為液滴直徑的1/32。

圖1 物理模型Fig. 1 Physical model

為了確定網格分辨率或者網格尺寸對模擬結果的影響，設計了3種網格局部加密方案，分別采用5、6、7級局部網格，最細的網格尺寸分別為液滴尺寸的1/16、1/32和1/64(圖2)。在韋伯數We=1.1和雷諾數Re=1 150的條件下，測試了網格加密對曳力系數的影響。圖3為曳力系數隨時間的變化。顯然，隨著網格局部加密級數的增加，曳力變得更加連續和穩定。六級加密方案結果與七級加密方案結果之間的差異很小。因此，本研究的所有案例都采用六級網格局部加密方案。

圖2 網格局部細化方案Fig. 2 Local grid refinement scheme

圖3 網格加密級數對曳力系數的影響(We=1.1，Re=1 150)Fig. 3 Influence of grid refinement grade on drag coefficient(We=1.1，Re=1 150)

1.4 模型準確性驗證

Reinhare(1964)［25］測量了空氣中下落水滴的終極曳力系數CD和橢球度E(液滴高度與液滴變形后橫向直徑的比值)，Loth(2008)［26］對其實驗數據進行擬合，提出曳力系數和橢球度經驗式。根據Reinhare(1964)的實驗條件，擬定了本次數值模擬的物性參數：氣液界面張力0.07 N/m，液滴尺寸4 mm，液體黏度2.0×10-4Pa·s，氣體黏度1.67×10-5Pa·s，液體密度1 000 kg/m3，氣體密度1.2 kg/m3。氣流速和對應的韋伯數We及液滴雷諾數Re、數值模擬得出的曳力系數CD和橢球度E與測試值對比如表1所示。曳力系數的平均絕對百分誤差為9.98%，橢球度的平均絕對百分誤差為11.45%。除韋伯數為11.95時模擬得到的橢球度0.24遠小于實驗測量值0.32，其他模擬橢球度均與實驗測量值吻合。同時，Quan等(2006)［16］和Helenbrook等(2002)［17］在韋伯數為12條件下橢球度模擬值為0.2，與當前模擬值非常接近。

表1 模擬結果與實測值對比Table 1 Comparison between simulation results and measurement values

圖4對比了Reinhare (1964)［25］實驗測量的空氣中下落水滴的形狀和本文模型數值模擬的液滴形狀，可以看出，兩者形狀非常吻合。

圖4 測量液滴形狀與模擬液滴形狀對比Fig. 4 Comparison between measured drop shapes and simulation results

2 氣井井筒條件下單液滴動力學特征

2.1 物性參數

某氣井井底流體壓力為4 MPa、溫度為353 K。根據高溫高壓物性計算方法，氣液界面張力取0.05 N/m，液滴尺寸4 mm，液體黏度3.6×10-4Pa·s，氣體黏度1.24×10-5Pa·s，液體密度972 kg/m3，氣體密度27.1 kg/m3。氣流速分別取0.68、1.36、1.78、2.36 m/s，對應的韋伯數為1、4、6.8、12.1，對應的液滴雷諾數Re為594、1 171、1 556和2 063。

2.2 液滴變形及破碎規律

從圖5可看出，在韋伯數為1時，可觀察到比較顯著的變形。隨韋伯數增加，液滴變形程度增加，液滴橢球度增加，逐漸變形為扁平體。當韋伯數增加到12.1時，扁平體在縱向上變得很薄，形狀類似于 “馕餅”；由于氣流的作用，“馕餅”向內凹陷，并發展成中空的袋狀，最后破碎，如圖5(d)所示。模擬中發現液滴破碎的臨界韋伯數約12，該值與低雷諾數條件下實驗測試或數值模擬值［16］接近，為此，氣井井筒條件下液滴破裂的臨界值可取12。

圖5 不同條件下單液滴變形過程模擬Fig. 5 Simulated deformation process of single drop under different conditions

2.3 液滴曳力系數

圖6給出了對應韋伯數條件下液滴曳力系數及橢球度隨時間的變化過程，可以看出，韋伯數等于1時，液滴的曳力系數從0.42逐漸增加，最后穩定在0.57；韋伯數等于4時，液滴的曳力系數從0.36逐漸增加，最后穩定在0.76；韋伯數等于6.8時，液滴的曳力系數從0.39逐漸增加，最后穩定在0.97；韋伯數等于12.1時，最大曳力系數在液滴破碎時刻，此時液滴迎風面積達到最大，氣流所施加的曳力最大。對于韋伯數略小于12的液滴，液滴僅變形但不破碎，根據圖6(d)中曳力系數隨液滴變形程度的變化規律，曳力系數穩定在1.3附近，變形達到穩定階段，因此建議曳力系數取1.3。不同韋伯數下液滴變形達到穩定后的曳力系數分析表明，液滴變形達到穩定后，最終的曳力系數隨韋伯數增加而增加。

根據標準曳力系數(Reinhare，1964)［25］，即不變形固體顆粒的曳力系數，當液滴雷諾數Re＞1 000時，固體球體的曳力系數為0.424，與韋伯數為1的模擬值吻合很好。同時，根據Rodi等(2002)［27］提出的液滴加速狀態下的曳力系數計算式，Re＞1 530時的曳力系數為2.62，略大于當前數值模擬的曳力系數1.70。這里，標準曳力系數0.424既不包含加速度作用也不包含液滴變形作用對曳力系數的影響，而2.62包括了這2種作用的影響，而本次數值模擬僅考慮了液滴變形作用的影響。因此，對于韋伯數低于臨界韋伯數的情況，液滴變形穩定后曳力系數介于標準曳力系數和加速狀態曳力系數之間，液滴加速將使曳力系數增加。

3 單液滴攜帶臨界氣流速計算模型

數值模擬得到的液滴曳力系數、臨界韋伯數為液滴攜帶臨界氣流速研究提供了重要基礎參數。

3.1 模型的建立

氣井中的液滴主要受自身重力、氣流對液滴的曳力和浮力作用。由液滴質點受力可知，當氣流速足夠大，液滴將懸浮在氣流中，此時重力等于浮力與曳力之和，即滿足式(17)。

圖6 液滴曳力系數CD隨時間的變化Fig. 6 Variation of drop drag coefficient (CD) with time

其中，C定義為液滴攜帶臨界氣流速綜合系數，表征液滴尺寸和液滴變形對液滴攜帶臨界氣流速的影響。

根據模擬得到的液滴在不同韋伯數下的曳力系數，計算了不同韋伯數下的綜合系數C。從表2可知，氣流場中不同韋伯數下臨界氣流速綜合系數在2.19～3.31之間變化，綜合系數取3.31計算得到的流速可將流場中所有大小液滴帶走。

Turner模型假設液滴為圓球狀，將液滴曳力系數和臨界韋伯數分別取為0.44和30，所得的系數為5.5，為安全起見，Turner模型將系數上調20%，最后建議系數為6.5。Coleman認為在低壓氣井條件下，該系數可不上調，建議系數取5.5。李閩假設液滴為橢球體，將曳力系數取為1，最終的模型系數為2.5。結合本研究的數值模擬結果來看，Turner模型和Coleman模型的綜合系數偏保守，所計算的氣井攜液臨界氣流速是液滴攜帶所需臨界氣流速的1.96倍和1.76倍。李閩模型所計算的氣井攜液臨界氣流速與液滴攜帶所需臨界氣流速接近。

表2 不同韋伯數下的液滴攜帶氣流速綜合系數Table 2 Composite coefficient of drop carrying gas flow rate at different Weber numbers

但是，從攜液模型的應用來看，Turner模型和Coleman模型在國內外大液量氣井廣泛應用，而李閩模型在國內低滲低液量氣井廣泛認可。從氣井井筒攜液機理的角度分析，大液量氣井的積液本質是液膜回流，氣井連續攜液臨界氣流速由液膜攜帶的臨界氣流速控制，液膜攜帶的臨界氣流速比液滴攜帶的臨界氣流速大得多，而低液量氣井連續攜液臨界氣流速由液滴攜帶的臨界氣流速控制。Turner模型和Coleman模型計算的氣井攜液臨界氣流速與液膜攜帶的臨界氣流速接近，李閩模型計算的氣井攜液臨界氣流速與單液滴攜帶的臨界氣流速接近。

氣井攜液臨界氣流速在什么條件下由液滴攜帶控制向液膜攜帶控制轉化，與氣流場中液滴的量及液滴間距有關。當氣井井筒液流量增加后，液滴間距減小，液滴相互作用會削弱氣流對液滴的有效作用力，為此液滴攜帶的氣流速將增加。對于稠密多液滴相互作用機理及攜帶規律將是較大液量氣井攜液機理研究的重要方向。

3.2 模型準確性驗證

Coleman等(1991)［3］文獻試驗數據共56組，油管內徑均為62 mm，其中13組數據有明顯的段塞特征，Coleman等判定為積液井，有9口井無產液量或產液量為0，有1口井油壓數據為0，這23口井數據均未用于模型評價。考慮到氣井液流量較大時，液滴數量多，液滴間距小，液滴相互作用明顯，曳力系數將受影響，為了較真實反映單液滴模型準確性，選用了氣井產液量小于1 m3/d的14口井的數據進行評價。模型計算值與測試值對比如圖7所示。

圖7 現有單液滴模型攜液臨界氣流量計算值與測試值對比Fig. 7 Comparison between the critical drop carrying gas flow rate calculated by existing single drop models and the measurement result

從圖7可知，李閩橢球模型攜液臨界氣流量較測試值偏小，Coleman圓球模型和Turner圓球模型計算值較測試值偏大很多。根據本研究的結果，將綜合系數取為3.31計算的攜液臨界氣流量與測試值較為接近。需特別說明的是，綜合系數取3.31僅適合于低液量氣井。對于高液量氣井，或者液滴間距過小時，綜合系數取值需要考慮液滴相互作用的影響，屬于多液滴相互作用研究范疇，有待進一步研究。

4 結論

(1)建立了描述氣流場中液滴動力學特征的數值模型，用于模擬氣井高溫高壓高雷諾數條件下單液滴形狀特征、變形規律、破碎條件和數值計算曳力系數。所建數值模型具有較高準確性，與實驗結果吻合較好，曳力系數的平均絕對百分誤差為9.98%，液滴橢球度的平均絕對百分誤差為11.45%。

(2)在氣流速作用下，液滴逐漸變形為橢球體；在韋伯數為1時，觀測到液滴明顯變形；隨氣流速或韋伯數增大，液滴變形后的橢球度增大；當韋伯數增加到12時，液滴的界面張力不能克服氣流的氣動力，液滴破裂為更小液滴；對于井筒流場中液滴稀疏的低液量氣井，可將液滴破碎的臨界韋伯數取為12，以此來確定流場中最大液滴的尺寸。

(3)通過對液滴進行受力分析得到液滴攜帶臨界氣流速綜合系數，綜合系數隨韋伯數或液滴尺寸增加而增加，但增加幅度不大；液滴的韋伯數從1增加到12，或者液滴尺寸增加12倍時，綜合系數從2.19增加到3.31，僅增加了51%，說明液滴尺寸不是影響液滴攜帶臨界氣流速的關鍵。

(4)本文采用數值模擬研究得到的綜合系數3.31與李閩橢球模型中的系數2.5比較接近，但與Turner模型的系數6.5和Coleman模型的系數5.5相差較大。

(5) Turner模型和Coleman模型將液滴考慮為球狀、臨界韋伯數取為30、曳力系數取為0.44，與數值模擬研究的結果相差較大，說明這兩種模型存在一定的理論缺陷，或者說在臨界氣流速預測時太過于保守。

(6)單液滴模型準確性評價表明，綜合系數3.31適用于低液量(小于1 m3/d)氣井，此時井筒液滴稀疏，液滴相互作用可忽略不計；對于大液量氣井，流場中液滴數量大，液滴間距小，液滴相互作用會削弱氣流對液滴的有效作用力，單液滴模型不適用。對于氣井產液量多大，或者液滴間距多小時，液滴攜帶臨界氣流速計算需要考慮液滴相互作用的影響，屬于多液滴相互作用研究范疇，有待進一步研究。