產水氣井井筒溫度壓力計算方法

李 波， 甯 波， 蘇海洋， 劉 虹， 位云生

（中國石油勘探開發研究院，北京 100083）

產水氣井井筒溫度壓力計算方法

李 波， 甯 波， 蘇海洋， 劉 虹， 位云生

（中國石油勘探開發研究院，北京 100083）

應用考慮氣液相間滑脫與流體沿井筒截面非均勻分布的漂移模型，結合動量守恒、能量守恒和井筒傳熱學，建立考慮井斜變化的氣井井筒溫度、壓力耦合預測模型，通過15口實測井數據對模型進行驗證.結果表明：所建模型可以準確計算產水氣井井筒的溫度和壓力，計算結果平均絕對誤差3.60%，滿足工程計算要求.

多相流；漂移模型；溫度；壓力；致密砂巖氣井

0 引言

地層流體（油、氣、水）從井底向井口流動過程中往往伴隨著許多復雜的物理現象，流動型態一般包括泡狀流、段塞流、過渡流和環霧流［1－7］.井筒多相流動常用的處理方法有三種：經驗模型，均質模型和機理模型.經驗模型根據多相流體運動的物理概念和物性參數，利用因次分析法或根據基本運動方程得到描述某種特殊流動的無因次參數建立多相流動的經驗關系式.截止目前，國內外已有許多經驗模型［1－7］，但在實際運用中，經驗模型需要大量的計算過程來判斷流型.此外，在某些流型之間的過渡區可能出現壓降和持液率的不連續現象，這樣就使得在求解井筒壓力和持液率等參數時出現不收斂或不連續.均質模型采用單相流的處理方法，認為氣體和液體的流速相等，但這與實際情況不符，特別是在氣體流動含量較低、流動速度較慢的情況下［8－10］.機理模型通常被認為是最精確的模型，因為它是在詳細描述不同流型的物理現象的基礎上提出的，但是其方程系統復雜，求解困難，以及在某些流型的過渡區可能出現壓降和持液率不連續從而使計算不收斂，因此沒有被廣泛運用.

由于漂移模型具有相對簡單、連續可微的優點，適用于井筒多相流動建模.為此，對漂移模型開展了大量的研究［8－16］，本文對漂移模型的重要參數（分布系數和漂移速度）進行比較優選，建立適合井筒多相流動的漂移模型.通過漂移模型計算井筒任意深度的截面含氣率及流體物性參數，再根據井筒傳熱學、動量守恒和能量守恒，建立井筒溫度和壓力耦合模型，通過分段迭代法對模型進行求解，并運用15口產水氣井實測數據對模型進行驗證，計算一口實例水平井井筒壓力、溫度分布.

1 漂移模型

在井筒氣液兩相流動過程中，由于氣相和液相的密度不一樣，存在氣相超越液相流動的現象，即氣液兩相間存在滑脫，因此采用簡單的均相模型來計算井筒的溫度和壓力分布會有一定的偏差.1965年，Zuber＆Findlay［8］首次提出考慮了氣液相間的滑脫現象，以及氣相沿井筒徑向的變密度分布的漂移模型.而后，Bankoff［9］，Wallis［10］，Ishii［15－16］，Hasan＆Kabir［19－20］，Ansari［5］等對漂移模型進行了改進，使其廣泛用于核工業和石油工業中的多相流動建模.由于漂移模型具有連續、可微和相對簡單的特點，Holmes等［13］于1997年將漂移模型引入數值模擬復雜結構井井筒等溫多相流動建模，但是由于之前的模型參數大多適合小直徑管流，不能直接用于井筒或地面設備中的多相流動建模［14－15］.為此，Oddie＆Shi等［14－16］開展了大量的實驗，對模型參數進行了回歸分析，并建立了適合大尺寸井筒多相流動的漂移模型.本研究的目的是通過研究前人所建立的漂移模型，建立適合大尺寸井筒非等溫多相流動的漂移模型，為準確計算井筒溫度和壓力分布提供前提.

1.1 漂移模型的基本關系式

氣液兩相流動情況下，氣相與液相間的滑脫現象是由于以下兩種機理結合產生的.一是氣體管橫截面上速度剖面和相分布剖面不一致，氣相在管中央處分布最高，管中央處混合物的流動速度也最快，如圖1所示.因此，整個截面上氣相的平均速度要比液相的平均速度大.另外一種機理就是氣體受浮力的作用，傾向于垂直向上運動.

考慮以上兩種機理的模型，氣體的截面平均流速可表示為

式中，vg為氣體的截面平均速度，m·s－1；Co為分布系數（描述兩相的分布特性，即流動型態特性），無因次；vm為混合物的截面平均流速（無滑脫），m·s－1；vd為氣體的漂移速度（描述浮力效應），m·s－1.

混合物的截面平均速度是氣相和液相折算速度之和，即

式中，vsg為氣體折算速度，m·s－1；vsl為液體折算速度，m·s－1；αg為截面平均含氣率；vl為液相的流動速度，m·s－1.

圖1 漂移模型示意圖Fig.1 Drift-flux model

1.2 漂移模型的關鍵參數

1.2.1 分布系數Co

Zuber＆Findlay給出Co的取值范圍是1～1.5［8］.對于泡流和段塞流，Aziz、Ansari和Hasan＆Kabir給出Co取值為1.2；而對于環霧流，Co取值接近于1.0［4－5］.此外，在αg為1.0（高氣液比）時，Co的值也應為1.0.與此同時，αgCo的值不能超過1.0.因此，Co表達式必須使其值在泡流和段塞流情況下有恒定不變的值，并且在αg為1.0時或混合物速度增大時，Co取值減小為1.0，即

為了使Co在αg和vm較大時取值減小到1.0，γ采用如下表達式

在以上條件不能滿足時，可能會出現不穩定的情況.在B＜（2－A）／A時，以上兩個條件都成立.

1.2.2 漂移速度vd

漂移速度是氣液混合物流動過程中氣液相對運動的截面平均速度值.在截面含氣率αg較小時，氣體漂

式中，vsgf為淹沒速度，即氣相折算流速到達足以使沿管壁液膜不會掉落下來時的速度值，流體的流動將出現環狀流．參數A，B和Fv需要通過實驗進行調整，根據Oddie＆Shi等［14－16］人的實驗及模型參數優化結果，氣液兩相流情況下，A取1.2，B取0.3，Fv取1.0.

淹沒速度用Wallis等［3］給出的關系式進行求解

式中，ku是臨界Kutateladze數，具體算法請參考文獻［15］；vc為特征流速，由式（8）進行計算.

顯然，氣相折算速度（vsg＝αgvm）隨αg和vm的增加而增加.因此，Co必須滿足以下關系式移速度與氣體從靜止液體中向上移動速度有關，Harmathy［11］得到其值為1.53 vc，vc表達式為

式中，σgl為氣－液相的表面張力，N·m－1；g為重力加速度，m·s－2；ρl為液相密度，kg·m－3；ρg為氣相密度，kg·m－3.

在αg取較大值時，可能出現氣液相沿不同方向流動的逆流情況［14］.在逆流的限制條件（vsg＞0，vsl＜0）的情況下，漂移速度是αg的函數，即有

對于處在以上兩種極端情況之間的過渡區，則通過選用兩個臨界的截面含氣率a1和a2，然后插值進行確定漂移速度.根據Oddie＆Shi等［14－16］人的實驗結果，a1取0.06，a2取0.21.

綜合考慮以上情況，漂移速度的總表達式為

式中，當αg≤a1時，K（αg）＝1.53／Co；當αg≥a2時，K（αg）＝Ku（）；在a1＜αg＜a2時，通過線性插值得到其值.

對于傾斜管流，鑒于Beggs＆Brill對傾斜管流的處理方法，基于垂直管中的流動情況，Hasan＆Kabir［4］提出采用下式進行修正

其中，θ為傾斜角，（°）；m（θ）為傾斜影響因子；vdθ是傾角為θ時的漂移速度，m·s－1；vd0為垂直流下的漂移速度，m·s－1.

對于氣液兩相流動系統，Hasan-Kabir［4］提出m（θ）采用下面式子進行計算

根據Oddie＆Shi等［14］在大管徑傾斜管流實驗，傾斜影響因子采用下式計算結果與實驗結果吻合更好

由以上研究可知，通過確定出漂移模型的兩個重要參數（分布系數和漂移速度），便可以得到井筒的截面處含氣率，進而確定該井筒處流體的平均物性參數.但從式（4）、（3）、（10）可知：分布系數、漂移速度和截面含氣率之間互為隱函數，因此需要進行迭代求解.

2 井筒溫度壓力計算模型

地層流體從井底流向井口過程中，不斷向地層中散發熱量，以及需要克服重力和摩擦力作用，井筒溫度、

壓力不斷降低，而流體物性又與溫度、壓力之間相互影響，因此需要建立溫度、壓力和物性耦合模型，采用迭代法進行求解.

2.1 井筒壓力計算模型

油氣井井筒壓力分布計算是油氣井生產設計、動態分析和預測的核心問題.在此，假設流體在井筒中的流動為一維穩定流動，從管流系統任取一控制單元，建立如圖2所示的坐標系.由動量守恒定理有

式中，ρm為混合物密度，kg·m－3；A為井筒截面，m2；D為井筒直徑，m；θ為井筒傾斜角，（°）；vm為混合物流速，m·s－1；w為管壁摩擦力，N；dp為壓力變化量，Pa；dz為井筒長度，m.

整理得到壓力梯度

圖2 井筒一維流動示意圖Fig.2 One-dimensional flow in wellbores

通過引入摩擦阻力系數fm，壓力梯度方程可表示為

式中，ε為管壁粗糙度，mm；D為管徑，mm；Λ為無因次參數，Λ＝（ε／D）1.109 8／2.825 7＋（7.149／Rem）0.898 1；Rem為無因次雷諾數，Rem＝Dυmρm／μm．

混合物的密度ρm和混合物的粘度μm分別采用體積加權和質量加權：

式中φ為截面體積含氣率（通過漂移模型計算），小數；x為截面質量含氣率，小數；μg、μl為氣體、液體粘度，mPa·s.

2.2 井筒溫度計算模型

地層中的流體通過井筒從井底流到井口的過程中，地層溫度與井筒內流體之間存在溫度差，井筒中流體與地層之間存在熱量交換，此外，流體自身的各種能量也會相互轉化.根據Ramey、Willhite、Hasan＆Kabir、毛偉和廖新維等［17－23］人研究，井筒傳熱可作以下簡化假設：①流體在井筒中的流動為一維穩定流動；②氣液間不存在質量交換；③井筒與地層之間只進行熱量徑向傳遞，熱量徑向傳遞包括兩個過程：井筒到水泥外沿之間的穩態傳熱，以及水泥環外向地層深處的非穩態傳熱.

由能量守恒定律有

摩擦阻力系數fm可以采用Chen方法進行計算，

式中，h為流體的比焓，J；q為流體的徑向熱容量，J·s－1.

由熱力學第一定律，得到流體比焓梯度

式中，cP為流體的定壓比熱容，J·（kg·℃）－1；Tf為油管內流體流動溫度，℃；αJ為焦耳—湯姆遜系數，℃·Pa－1.

在單位長度井段上，流體從油管到井壁的熱流量梯度為

式中，Tf為井筒流體溫度，℃；Th為井壁溫度，℃；Uto為總傳熱系數，W·（m·℃）－1；Wm為質量流量，kg·s－1.

根據Hasan＆Kabir［17］提出的無因次時間函數f（tD），引入松弛因子A，由式（20）～（22）可得到井筒溫度梯度方程為

式中，Te是地層溫度，℃；ke是地層導熱系數，W·（m·℃）－1；rto是套管外徑，m.考慮焦耳—湯姆遜效應，穩定流動情況下（23）式的解析解為

式中，Teout、Tein是計算段的出口和入口地層溫度，℃；zout、zin為計算段出口和入口深度，m；gT為計算井段處的地溫梯度，℃·m－1.

2.3 井筒溫度、壓力計算

由于井筒中流體物性、溫度和壓力之間關系緊密，相互影響，因此在計算井筒溫度和壓力時，需要進行耦合迭代求解.迭代求解面臨一個重要的問題是如何給定初值，因為初始值會影響迭代的收斂速度，本文采用（25）式計算井筒流出初始溫度值.

溫度、壓力耦合迭代計算詳細步驟如下：

1）給定井底條件：Tfin＝TeBH，Tein＝TeBH，pin＝pBH，并選取合適的ΔL將井筒分為n段.

2）根據（25）式計算粗略出口溫度Tfout.

3）估計ΔL段壓力變化量Δp估計，利用漂移模型計算各相流體的流動速度和持率.

5）根據（16）式計算ΔL段壓降梯度，然后計算出ΔL段的壓力變化量Δp計算.

6）比較壓力變化估計量Δp估計和計算值Δp計算，若二者之差不滿足條件，則將值Δp計算作為初值，重復3）～6）步，直到滿足條件為止.

7）計算井筒總傳熱系數和焦耳—湯姆遜系數，利用（24）式精確計算段出口段的溫度值.

8）比較第3）步與第7）步的溫度變化值，若二者之差滿足要求，則進行下一段溫度壓力計算；否則，將第7）步計算的溫度變化值作為溫度初始值，跳到第2）步進行，重新計算，直到滿足要求為止.

9）計算該計算段對應的深度Li、壓力pi和出口溫度Ti

10）將Li處的壓力pi和Tf作為下一計算段的入口參數，重復2）～9）步計算下一深度Li＋1的壓力pi＋1和溫度Ti＋1，直到算完整個井筒溫度和壓力．

3 模型驗證與實例計算

3.1 模型驗證

A氣田為典型致密砂巖氣藏，氣藏平均壓力系數為0.87，溫度110℃，天然氣相對密度0.598.氣藏油水關系復雜，氣井普遍產水，產出水包括凝析水和層間水.取A氣田15口典型氣井測試資料對所建模型進行評價，氣井產氣量為（3.98～14.59）×104m3·d－1，產水量為（1.11～18.34）m3·d－1，平均氣水比1.14× 104m3·m－3，井底流壓（21.54～28.62）MPa，測試井口壓（21.54～28.62）MPa，測試井口壓力（11.43～22.43）MPa，井的基本參數與實測如表1所示.油管直徑62 mm，管壁相對粗糙度0.004 5.

由模型計算井口壓力與實測井口壓力對比可知（表2、表3、圖3），模型計算井口壓力與實測壓力誤差均小于10%，平均誤差－0.38%，平均誤差標準偏差為4.21%，平均絕對誤差為3.60%，平均絕對誤差標準差為2.22%，是參與評價模型中誤差最小的，其次是Hagedorn-Brown計算結果較好.而在油井管流計算中常用的Beggs-Brill和Orkiszewsiki模型，計算結果誤差較大，原因是這兩種模型都是通過實驗建立的經驗關系式來計算井筒截面處的流體物性，受實驗參數范圍（如：流量、氣液比等）和實驗條件（如：管徑、管長等）的限制，這些經驗關系式對低流速、較低氣液比的油井多相流計算具有較高精度，而對高流量、高氣液比氣井則適應性較差，計算結果誤差較大.

表1 氣井基本參數表Table 1 Basic data of gas wells

表2 模型計算井口壓力與實測井口壓力對比Table 2 Comparison between calculated wellhead pressures and measured wellhead pressures

通過所建模型計算還得到了井口溫度，但由于缺少現場井口溫度測試數據，不能進行溫度誤差分析.根據文獻［21］的觀點，如果溫度和壓力是通過耦合計算得到的，只要壓力計算結果準確，可以認為溫度計算結果也可靠.此外，利用本文所建模型對文獻［23］中的實例氣井進行計算，計算溫度、壓力值與實測溫度、壓力值誤差分別為1.54%和2.34%，說明本文計算方法是準確的.因此，可以認為本文建立的溫度壓力預測模型可靠，滿足工程計算要求.

表3 計算結果誤差分析Table 3 Error analysis of results

3.2 實例井計算

長深D平6井位于長嶺斷陷中部隆起帶哈爾金構造上，完鉆井深4 207 m，水平段長1 007 m，水平段中深3 864 m，采用水平段裸眼封隔器＋滑套完井.為滿足產氣與臨界攜液要求，采用73 mm油管進行生產.該井壓裂投產后日產氣19.53×104m3，水氣比1.77 m3·（104m3）－1，天然氣相對密度為0.59，井底流壓33.50 MPa，井底溫度133℃.

由井筒溫度與壓力分布計算結果（圖4）可知，地層流體從井底向地面流動過程中，由于克服重力、摩擦力和速度變化，井筒壓力不斷降低；井筒溫度不斷降低，但始終大于地層溫度，并且差值越來越大，這主要是由于流體從地層深部攜帶的熱量由于溫差需要不斷向地層中散發熱量造成的.通過計算：流體達到井口后溫度為65.68℃，油壓為20.61 MPa，計算油壓與井口實測油壓20.42 MPa，誤差僅為1.89%，模型計算精度較高.此外，在水平井的水平段中點至造斜段（3 120 m～3 864 m），井筒和溫度下降較慢；從井底到井口的總壓降為12.892 0 MPa，其中加速度壓降為0.001 9 MPa，摩擦壓降3.801 3 MPa，重力壓降9.088 8 MPa，可以看出井筒中的壓降主要是由摩擦力和重力引起的，而加速壓降則相對較小.

圖3 模型計算壓力與實測井口壓力對比Fig.3 Comparison of calculated wellhead pressure with measured values

4 結論

1）漂移模型考慮了氣液相間的滑脫和流體沿井筒截面的非均勻分布，且具有相對簡單、連續可微的特點，可以準確計算任意傾角情況下井筒截面處的流體分布，適合用于井筒多相流動建模.

2）基于漂移模型的井筒溫度、壓力耦合預測模型可以計算產水氣井井筒的溫度和壓力分布，以及流體的物性參數，計算過程相對簡單.對高氣液比產水氣井，所建模型計算精度高于常用的Hagedorn-Brown、Beggs-Brill和Orkiszewsiki模型，平均絕對誤差僅為3.60%，滿足工程計算要求.

圖4 井筒溫度和壓力分布Fig.4 Temperature and pressure distribution in wellbores

［1］Hagedorn A R，Brown K E.Experimental study of pressure gradients occurring during continuous two-phase flow in small diameter vertical conduits［J］.JPT，1965，（4）：475－484.

［2］Orkiszewski J.Predicting two-phase pressure drops in vertical pipes［J］.JPT，1967，23（6）：829－838.

［3］Beggs H D，Brill J P.A study of tow-phase flow in inclined pipes［J］.JPT，1973，25（5）：607－617.

［4］Hasan A R，Kabir C S.Predicting multiphase flow behavior in a deviated well［J］.SPE Prod Eng，1988，474（11）：475－482.

［5］Ansari A M，Sylvester N D，Sarica，et al.A comprehensive mechanistic model for upward two-phase flow in wellbores［J］. SPE Prod＆Faci，1994，9（2）：143－151.

［6］Gomez L E，Ovadia S，Zelimir S，et al.Unified mechanistic model for steady-state two-phase flow：Horizontal to vertical upward flow［J］.SPE Journal，2000，5（3）：339－350.

［7］Kaya A S，Sarica C，Bill J P.Mechanistic modeling of two phase flow in deviated wells［J］.SPE Prod＆Faci，2001，（9）：156－166.

［8］Zuber N，Findlay J.Average volumetric concentration in two-phase flow systems［J］.Heat Transfer Trans ASME，1965，87：453－468.

［9］Bankoff S G.A variable density，single fluid model for two-phase flow with particular reference to steam-water flow［J］.Heat Transfer Trans ASME，1960，82：265－272.

［10］Wallis G B，Makkenchery S.The hanging film phenomenon in vertical annular two-phase flow［J］.Trans ASME，1974：186－197.

［11］Harmathy T Z.Velocity of large drops and bubbles in media of infinite and restricted extent［J］.AIChE J，1960，6（2）：281－288.

［12］Hibiki T，Ishii M.One-dimensional drift-flux model for two-phase flow in a large diameter pipe［J］.International Journal of Heat and Mass Transfer，2003，46：1773－1790.

［13］Holmes J A，Barkve T，Lund O.Application of a multi-segment well model to simulate flow in advanced wells［R］.SPE 50646，1998.

［14］Oddie G，Shi H，Durlofsky，et al.Experimental study of two and three phase flows in large diameter inclined pipes［J］. Multiphase Flow，2003，29：527－558.

［15］Shi H，Holmes J A，Durlofsky，et al.Drift-flux modeling of two-phase flow in wellbores［J］.SPEJ，2005，10（1）：24－33. ［16］Shi H，Holmes J A，Diaz L R，et al.Drift-flux parameters for three-phase steady-state flow in wellbores［J］.SPEJ，2005，10 （2）：130－137.

［17］Ramey H J.Wellbore heat transmission［J］.Trans，AIME，1962，14（4）：427－435.

［18］Willhite G P.Over-all heat transfer coefficients in stream and hot water injection wells［J］.JPT，1967：607－615.

［19］Hasan A R，Kabir C S.Aspects of heat transfer during two-phase flow in wellbores［J］.SPE Prod＆Faci，1994，9（3）：211 －216.

［20］Hasan A R，Kabir C S，Wang X.A robust steady-state model for flowing-fluid temperature model in complex wells［J］.SPE Prod＆Oper，2009，24（2），269－276.

［21］Mao Wei，Liang Zheng.New method of calculating wellbore temperature in gas wells［J］.Journal of Southwest Petroleum Institute，1999，21（1）：56－66.

［22］Mao Wei，Liang Zheng.Coupling analysis of the pressure and temperature in gas well borehole［J］.Natural Gas Industry，1999，19（6）：66－68.

［23］Liao Xinwei，Feng Jilei.Pressure-temperature coupling calculation of transient wellbore heat transfer in deep geopressured gas reservoir［J］.Petroleum Exploration and Development，2005，32（1）：67－69.

［24］Wu Xiaodong，Wu Han，Han Guoqing，et al.A new model for calculating wellbore temperature and pressure distribution of a high-H2S gas well considering the influence of the sulfur release in wellbores［J］.Natural Gas Industry，2011，31（9）：69－72.

A Calculation Method for Temperature and Pressure Distribution in Water Production Gas Wellbores

LI Bo，NING Bo，SU Haiyang，LIU Hong，WEI Yunsheng
（Research Institute of Petroleum Exploration＆Development，Beijing 100083，China）

Based on drift model which copes with gas-liquid phase slippage and non-uniform distribution along section of wellbore，a coupled temperature-pressure prediction model was established for varied angle gas wellbores considering momentum conservation，conservation of energy and wellbore heat transmission.The model was validated with 15 tested gas wells.The calculation showed satisfying results with average absolute error of 3.6%which meets engineering requirements.

multi-phase flow；drift model；temperature；pressure；tight sandstone gas wells

date：2013－08－21；Revised date：2014－03－04

TE319

A

2013－08－21；

2014－03－04

國家科技重大專項（2011ZX05015）資助項目

李波（1985－），男，貴州思南，博士生，從事氣藏工程與多相管流研究，E-mail：libosonova＠163.com

1001-246X（2014）05-0573-08