張智，王嘉偉，李炎軍，羅鳴，張超

（1.西南石油大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，成都 610500；2.中海石油（中國）有限公司湛江分公司，廣東湛江 524057）

0 引言

隨著高溫、高壓、高產(chǎn)氣井不斷投入開發(fā)，開采難度不斷增大，對井筒完整性提出了更高的要求[1-5]。尤其是高產(chǎn)氣井（無阻流量大于120×104m3/d）在開采過程中，受地質(zhì)和油管條件等因素影響，產(chǎn)量隨時間一直波動，且常面臨開關(guān)井工況[6-7]。高產(chǎn)氣井產(chǎn)量高，流速快，開關(guān)井速度對油管壓力分布影響較大。開關(guān)井速度過快將導(dǎo)致較大的波動壓力，誘發(fā)嚴(yán)重的井筒完整性問題[8]；而開關(guān)井速度過慢，則開關(guān)井操作時間過長，會對現(xiàn)場作業(yè)造成較大的影響。1981年10月4日，德州一口氣井發(fā)生井噴，事故原因主要是開關(guān)井速度過快形成水錘產(chǎn)生高壓所致[9]。四川盆地東北部的河壩 1井試采初期因關(guān)井過快嚴(yán)重影響了氣井的正常生產(chǎn)，甚至出現(xiàn)天然氣爆泄險情[10]。因此，研究開關(guān)井時油管內(nèi)的瞬變流動過程、根據(jù)實際生產(chǎn)優(yōu)化開關(guān)井時間對保障井筒完整性至關(guān)重要。

油管內(nèi)的流動可看作是非恒定有壓管流，在這種流態(tài)下，若流體流速驟然發(fā)生改變，會在油管內(nèi)產(chǎn)生劇烈的壓力波動[11]。氣井生產(chǎn)中，采氣樹閥門開度發(fā)生變化時，油管內(nèi)流體的速度隨之發(fā)生改變，進而產(chǎn)生波動壓力。開度變化的急緩影響波動壓力的大小，當(dāng)閥門急驟地啟閉時，油管內(nèi)的波動壓力可能達到正常壓力的數(shù)倍甚至數(shù)十倍，且壓力的反復(fù)變化會引起油管和設(shè)備的振動，嚴(yán)重時會造成油管及井口裝置的損壞[12]而引發(fā)事故，影響井筒的完整性及安全生產(chǎn)。

閥門瞬關(guān)產(chǎn)生的波動壓力也被稱作水錘。水錘現(xiàn)象表現(xiàn)出了流體的可壓縮性，是一個瞬變流問題，它體現(xiàn)出流體在運動過程中，流域內(nèi)的所有運動物理量（速度、壓力等）兼具時變與位變的性質(zhì)。目前，水錘問題得到國內(nèi)外學(xué)者的高度重視，該問題已成為流體力學(xué)瞬變流領(lǐng)域內(nèi)一大研究焦點。1993年Jardine等[13]提出了“硬關(guān)井”與“軟關(guān)井”問題，并指出了硬、軟關(guān)井操作的優(yōu)劣。Erika等[14]預(yù)測了氣井內(nèi)流體在充填過程中的瞬態(tài)現(xiàn)象。Jalal等[15]研究了一種微壓縮賓漢流體的流動，建立了基于水平井壓力和壓力導(dǎo)數(shù)的新解析方程，并對6種不同工況進行了逐步計算驗證。國內(nèi)學(xué)者對油管內(nèi)流體瞬變過程的研究更多體現(xiàn)在鉆井過程中。何世明等[16]用ADINA軟件對溢流關(guān)井產(chǎn)生的波動壓力進行了有限元仿真研究。彭齊等[17]以槽流模型為基礎(chǔ)，結(jié)合起下鉆過程中的流體真實速度分布情況，分別討論了層流、紊流狀態(tài)下的波動壓力，建立了基于鉆柱運動的穩(wěn)態(tài)井筒波動壓力計算模型。陳林等[18]提出了壓力波速標(biāo)準(zhǔn)圖版的繪制方法并繪制了壓力波速圖版，用以查詢油管氣液兩相流在不同溫壓和含水率下的壓力波速。

盡管對水錘的相關(guān)研究取得了眾多成果，但在氣井采氣階段，有關(guān)氣井瞬變流的研究均未討論閥門開度系數(shù)、持液率等參數(shù)對油管內(nèi)壓力的影響，同時油管造斜段部分流體瞬變流動情況研究的報道也很少。因此，建立適用于氣井復(fù)雜工況的瞬變流模型，并對油管造斜段進行單獨建模，分析討論閥門開度系數(shù)、持液率等參數(shù)對油管流體瞬態(tài)流動的影響很有必要。本文針對該問題，建立適用于氣井多相流環(huán)境的瞬變流數(shù)學(xué)模型，并對油管的造斜部分單獨進行建模，分析關(guān)井工況對油管內(nèi)流體瞬變流動的影響，為研究油管內(nèi)實際流動過程及氣井壓力分布[19-20]提供技術(shù)支持。

1 氣井瞬變流動數(shù)學(xué)模型

1.1 氣井關(guān)井瞬變流動過程

氣井油管內(nèi)的瞬變流動過程與普通管道類似，完成一次壓力的傳遞需要經(jīng)歷4個階段（一個周期）。設(shè)井深為L，壓力傳播速度為cm，井口流速為u0，井口壓力為p0，那么一次壓力傳遞分為 4個階段：①0～L/cm；②L/cm～2L/cm；②2L/cm～3L/cm；④3L/cm～4L/cm。

在井口閥門關(guān)閉的瞬間，緊貼閥門高度為Δz的微元流體速度會立刻變?yōu)?，流動滿足實際流體總流的伯努利方程。若從截面A1到A2定義一段元流（見圖1），由元流積分成總流，則可得實際流體總流的伯努利方程，如（1）式所示。

圖1 實際流管流體元流示意圖

該方程確立了實際流體總流流動中的勢能與動能、流速與壓強相互轉(zhuǎn)化的普遍規(guī)律。若取高度為Δz的微元流體作為研究對象，位置水頭和沿程摩阻就可以忽略不計，即該微元流體的速度水頭和壓力水頭之和為一定值。速度減小，壓力必然增大。

圖2為氣井形成水錘過程的4個階段。井口閥門關(guān)閉后，壓力波以波速cm向井底傳播，經(jīng)過t=L/cm后傳播到井底；油管內(nèi)流體壓力為p0+Δp（Δp=Δp1+Δp2+…+Δpi）。此時，油管內(nèi)壓力高于井底，靠近井底流體以速度um,0向井底倒流，上方流體依次流入井底。在L/cm～2L/cm內(nèi)，壓差逐漸消失并恢復(fù)到常壓。當(dāng)t=2L/cm時，流體在重力的作用下仍以速度um,0向井底流動，油管內(nèi)產(chǎn)生負壓，形成減壓波面并以波速cm向井底傳播，直至壓力降為p0?Δp。當(dāng)t=3L/cm時，井底壓力高于油管內(nèi)壓力，流體以速度um,0向上運動，井底的流體壓力恢復(fù)到p0。這種不平衡依次以波速cm向井口傳播，在t=4L/cm時傳到井口。至此，氣井關(guān)井產(chǎn)生的水錘壓力完成了一個周期的傳播。

圖2 氣井形成水錘的壓力波動傳播過程

1.2 氣井瞬變流控制方程

經(jīng)典瞬變流模型研究的是水平管道單相流，而且大多數(shù)僅考慮了液相。井下環(huán)境復(fù)雜多變，需要考慮諸多經(jīng)典模型中未考慮到的因素：在垂直或近似垂直的油管中，經(jīng)典模型無法準(zhǔn)確地描述重力帶來的影響；在造斜段，需考慮流體沖擊彎頭受到的反彈性膨脹力與該段因急變流造成的局部摩阻；油管中的傳播介質(zhì)是多相流體[21-22]，壓力波速與流體流速等參數(shù)會因相的增加發(fā)生很大的變化。無論是鉆井過程還是生產(chǎn)過程，此時單相流模型都無法滿足計算需求。

考慮多相流后，會涉及到油管內(nèi)流動介質(zhì)相間的相互作用，介質(zhì)物性受體積比、密度比等參數(shù)的影響，以及十分復(fù)雜的流型轉(zhuǎn)變和相間質(zhì)量、動量、能量傳遞過程等。因此，在建模時需作簡化：①油管內(nèi)的流動介質(zhì)與油管壁均為線彈性體，彈性模量恒定，油管縱向無彈性形變；②油管內(nèi)的流動介質(zhì)為一元流動，沿油管同一截面上的多相流體分布均勻，流速相同；③瞬變流中摩阻計算仍沿用油管恒定流公式。

研究油管內(nèi)瞬變流動過程的基礎(chǔ)是瞬變流動基本微分方程組。該微分方程組由連續(xù)性方程和運動方程構(gòu)成。如圖 3所示，于直井段和造斜段各取一有限長度（弧度）的流體段作為研究對象，長度（弧度）為dz，流體進口節(jié)點為j，流體出口節(jié)點為j+1。

圖3 油管有限單元的流入與流出

因為直井段與造斜段不影響流體的連續(xù)介質(zhì)假設(shè)，故兩種井段條件下采用同一連續(xù)性方程。根據(jù)所取有限元控制體在dt時間段內(nèi)的質(zhì)量守恒關(guān)系，有：

簡化整理可得：

在考慮多相流因素的分析中，達西公式失效，此時摩擦阻力需要通過范寧摩阻公式進行計算：

那么相應(yīng)的水錘壓力可表達為：

這里的cm是多相流條件下的壓力波速（此處不考慮固相），定義為：

在油管材料、尺寸確定后，如忽略油管內(nèi)氣液流體彈性模量的變化，則持液率的改變會直接影響壓力波速，從而影響該截面的水錘壓力。根據(jù)上式，持液率減小，壓力波速會增加，從而導(dǎo)致水錘壓力的增加。

由此可得到描述氣井瞬變流過程的連續(xù)方程為：同樣，于直井段和造斜段各取一有限長度（弧度）的流體段作為研究對象，長度（弧度）為dz，流體進口節(jié)點為j，流體出口節(jié)點為j+1，受力分析如圖4所示。

圖4 油管有限單元的受力分析

當(dāng)流體混合物流過直井段時，根據(jù)以上的受力分析，可以得到該段所受的合力為：

當(dāng)流體混合物流過油管造斜段時，由于慣性作用，流體會沖擊彎頭處，假設(shè)管壁為線彈性體，故會反作用給流體一個彈性膨脹力：

此時，造斜段所受的合力為：

由此，可以得到滿足氣井瞬變流過程的運動方程形式，對直井段和造斜段分別列出牛頓第二定律方程。

故氣井的瞬變流動過程可以分為 3個井段進行計算，如圖5所示。

圖5 氣井水錘分段計算示意圖

上述氣井瞬變流模型在經(jīng)典模型的基礎(chǔ)上，考慮了豎直油管的重力因素、多相流影響以及造斜井段對流體的作用，該模型能更好地反映井下復(fù)雜環(huán)境，對油管內(nèi)的流體瞬態(tài)運動情況做出精確的數(shù)學(xué)描述。

2 特征線法求解

對于瞬變流數(shù)學(xué)模型的求解，解析法與數(shù)值解法最為普遍。瞬變流模型是一對擬線性雙曲偏微分方程，求出解析解存在很大的困難。數(shù)值解法對方程本身的要求較低，可以不對方程做較多簡化，通過計算機模擬代替實際的模型實驗，更有利于對工程問題進行研究。本文采用數(shù)值解法中常見的特征線法進行求解。

直井段與造斜段需要分開進行描述，此處以直井段Ⅰ（垂直）為例，建立并求解特征線方程組，造斜段與直井段Ⅱ（斜）的求解步驟與直井段Ⅰ類似。

2.1 特征線方程組

引入待定系數(shù)ω對直井段簡化后的瞬變流基本微分方程進行線性化：

將線性化后的方程轉(zhuǎn)化為常微分方程，并選取合適的ω，由此可以得到描述瞬變流過程的特征方程組。順波特征線（C+）為：

逆波特征線（C?）為：

2.2 有限差分離散求解

通過有限差分法對井深和時間進行離散：將油管分為3段，分別是直井段Ⅰ、造斜段和直井段Ⅱ，油管微元長度為Δz；時間等分為M段，每段步長為Δt。以直井段Ⅰ為例，繪制特征網(wǎng)格L-t平面圖，如圖6所示。

圖6 特征網(wǎng)格L-t平面圖

若已知點P1和P2所在井深分別為L1和L2，i時刻相應(yīng)油管節(jié)點橫截面的壓力為pP1,i和pP2,i，流速為umP1,i和umP2,i，則點P在i+1時刻的壓力和流速可由（13）和（14）式聯(lián)立求得。重復(fù)此步驟，便可逐個求出后繼時刻油管內(nèi)所有節(jié)點的壓力與流速，造斜段和直井段Ⅱ的求解步驟與之相同。

分別沿特征線C+和C?對（13）和（14）式進行差分，可得i時刻、任意節(jié)點j橫截面的壓力和流速：

在氣井波動壓力和流速的求解過程中，油管每一節(jié)點的p和um初始值均為已知，先計算Δt時刻各節(jié)點的p和um，再計算下一時刻各節(jié)點的p和um，直到達到需要計算的時刻為止。

3 實例計算

3.1 計算分析

求解直井段與造斜段的聯(lián)合離散方程計算量巨大，無法進行人工求解，需要通過計算機程序來實現(xiàn)氣井瞬變流動的模擬計算。將實例井分為 3段計算，直井段Ⅰ離散為500段，造斜段離散為150段，直井段Ⅱ離散為 350段。定義井口閥門開度Kv滿足如下函數(shù)：

其中β為閥門開度系數(shù)，當(dāng)β=1.0時為線性關(guān)閥。

將標(biāo)準(zhǔn)狀況下的天然氣產(chǎn)量轉(zhuǎn)換為油管中的流量時，考慮到天然氣的壓縮性，需采用實際氣體狀態(tài)方程進行換算，其中天然氣壓縮因子Z是非常重要的參數(shù)。本文參照國家標(biāo)準(zhǔn)GB/T 17747[23]中的SGERG-88公式，使用高位發(fā)熱量、相對密度和 CO2含量作為輸入變量。在GB/T 17747中，用物性值計算天然氣壓縮因子的公式如下：

油管內(nèi)不同微元體天然氣壓縮因子的值可由上式聯(lián)立求得，（18）式中的B和C分別表示天然氣的第二、第三維利系數(shù)，是高位發(fā)熱量、相對密度、氣體混合物中不可燃和可燃的非烴組分（CO2、H2）的含量及溫度的函數(shù)，可根據(jù)GB/T 17747中的標(biāo)準(zhǔn)進行計算。

以XX6井為例模擬閥門關(guān)閉后井口的壓力波動情況，分析不同的閥門開度系數(shù)、截面持液率和關(guān)閥時間對井口波動壓力的影響。氣井的基礎(chǔ)數(shù)據(jù)見表1。

表1 XX6井物理及生產(chǎn)參數(shù)

3.2 結(jié)果分析

通過本文建立的多相流氣井瞬變流動分析模型，計算得到XX6井井口的壓力波動情況（見圖7）。由圖可見，閥門關(guān)閉瞬間，井口流體流速突變?yōu)?，壓力迅速增大，在3.28 s時達到峰值32.65 MPa，之后以不同的梯度開始下降，在9.06 s時達到谷值26.48 MPa，然后開始壓力波的下一次傳播。因油管內(nèi)的天然氣易壓縮，且在沿程摩擦阻力的作用下，波動壓力會很快衰減，經(jīng)過一段時間后，達到平衡壓力。

圖7 XX6井井口波動壓力模擬計算結(jié)果

閥門開度系數(shù)（β）是模擬瞬變流動過程一個重要的參數(shù)，它直接影響波動壓力的峰值大小和幅值變化。取不同β值模擬計算瞬變流過程中的井口壓力變化曲線（見圖8a），并對圖8a中紅色矩形框部分進行局部放大（見圖 8b）。由圖8a可知，當(dāng)β=1.5時，井口波動壓力的峰值最大，可達33.07 MPa，且β值越大，井口波動壓力的峰值越大。而圖8b顯示β值越大，井口波動壓力變化幅度越平緩，壓力傳播過程中的壓力突變區(qū)域越不明顯；β值越小（以β=0.5為例），井口波動壓力變化幅度越大，壓力突變區(qū)越明顯，突變區(qū)壓力峰值最大達到32.79 MPa，數(shù)值大于最先到達井口的波動壓力。雖然β值較大時會出現(xiàn)較大的波動壓力，但較小的β值會造成油管內(nèi)流體反復(fù)的壓力擾動，從而引起油管的流固耦聯(lián)振動，該問題是當(dāng)下極難妥善解決的技術(shù)難題。所以應(yīng)在不超過油管最大關(guān)井壓力的基礎(chǔ)上，盡量使用較大的閥門開度系數(shù)。

圖8 不同閥門開度系數(shù)下井口波動壓力模擬計算結(jié)果

截面持液率（HL）是研究多相流問題一個非常重要的參數(shù)，它是指在氣液兩相流動的過程中，液相的過流斷面面積占總過流面積的比例。本文選取緊貼井口的一層流體為研究截面，對該截面不同HL值條件下的壓力脈動進行計算（見圖9a），并對圖9a中紅色矩形框部分進行局部放大（見圖9b）。由計算結(jié)果可知，HL值越大，壓力波的傳播速度越大，傳播周期越短。例如，隨著HL值從0.03增加至0.10，壓力拐點分別出現(xiàn)在 3.47，3.28，3.08，2.79 s，表現(xiàn)為依次縮短（見圖9b中的紅色矩形框）。

較大的HL值對應(yīng)的波動壓力變化幅度大，會引起油管內(nèi)的壓力擾動，同時較大的HL值對應(yīng)的壓力值也較大，且HL值的變化會引起較明顯的波動壓力變化。實際生產(chǎn)中可通過精確計算油管內(nèi)流體的持液率，調(diào)整生產(chǎn)參數(shù)得到合適的HL值，控制波動壓力的大小和變化幅度，減小水錘壓力的沖擊。

圖9 不同截面持液率下井口壓力模擬計算結(jié)果

圖10 不同關(guān)閥時間下井口波動壓力模擬計算結(jié)果

不同關(guān)閥時間對井口波動壓力峰值也有較大影響，圖10a為3個不同關(guān)閥時間條件下井口波動壓力的模擬結(jié)果，3條虛線的位置代表3個關(guān)閥時間對應(yīng)的壓力峰值。可以看到，關(guān)閥時間增加，井口最大壓力值減小，峰值出現(xiàn)的時間也相應(yīng)滯后。圖10b為0.5～10.0 s關(guān)閥時間條件下，不同的井口波動壓力峰值與其對應(yīng)的時間?？梢园l(fā)現(xiàn)，短時間關(guān)井的壓力波傳播速度比較長時間關(guān)井的壓力波傳播速度要快。

XX6井是一口儲氣庫井，井淺且高產(chǎn)，閥門瞬關(guān)會形成較大的壓力波沖擊，井口壓力較大。隨著關(guān)閥時間增加，水錘沖擊壓力減小，達到平衡壓力的時間縮短。在深井或超深井中，沿程摩阻較大，天然氣本身具有較強的壓縮性，壓力波衰減更快，此時壓力波對井口的沖擊較小。

通過氣井瞬變流模型模擬氣井瞬態(tài)關(guān)井，可以優(yōu)化合理的閥門開度系數(shù)和關(guān)閥時間，減小水錘沖擊對井口裝置和油管造成的危害，保障井筒完整性。

4 結(jié)論

高產(chǎn)氣井井口閥門關(guān)閉時開度系數(shù)越大，井口壓力的峰值越大，波動壓力的變化幅度越平緩，壓力突變區(qū)越不明顯。在保證不超過油管最大關(guān)井壓力的前提下，使用較大的開度系數(shù)可減小壓力波的沖擊。

截面持液率越大，壓力波速越大，傳播周期越短；持液率越大，波動壓力的變化幅度越大，壓力越大。實際生產(chǎn)中，可通過調(diào)整生產(chǎn)參數(shù)得到合適的持液率，控制波動壓力的大小和變化幅度，減小水錘的沖擊。

采氣樹關(guān)閥時間增加，井口的最大波動壓力值減小，峰值出現(xiàn)的時間也相應(yīng)滯后，壓力突變區(qū)逐漸消失；閥門關(guān)閉時間越短，壓力波的傳播速度越快。

氣井瞬變流模型可以優(yōu)化合理的閥門開度系數(shù)和關(guān)閥時間，減小水錘沖擊對井口裝置和油管造成的危害，保障井筒的完整性。

符號注釋：

A，A1，A2——油管橫截面面積，m2；B——天然氣第二維利系數(shù)，m3/kmol；C——天然氣第三維利系數(shù)，m6/kmol2；C+——順波特征線；C?——逆波特征線；cm——壓力波速，m/s；D——油管外徑，m；e——油管壁厚，m；Eg——氣體彈性模量，Pa；EL——液體彈性模量，Pa；Ep——管壁彈性模量，Pa；f——范寧摩阻系數(shù)，無因次；Fbe——造斜段油管所受合外力，Pa；Fc——造斜段油管彈性膨脹力，Pa；Fs——油管膨脹造成的平均側(cè)向力，Pa；Fve——直井段油管所受合外力，Pa；g——重力加速度，m/s2；G——油管重力，N；hw,1,2——水頭位置 1到 2的沿程摩阻，m；HL——截面持液率，%；i——時間節(jié)點編號；j——空間節(jié)點編號；Kv——閥門開度，無因次；L——氣井井深，m；M——時間段數(shù)，段；p——油管內(nèi)壓力，Pa；p0——井口壓力，Pa；p1，p2——水頭位置 1，2處的壓力，Pa；pP1,i，pP2,i——i時刻點P1和P2處油管節(jié)點橫截面的壓力，Pa；P，P1，P2——特征網(wǎng)格L-t平面圖上的點；pin——微元體入口流壓，Pa；pout——微元體出口流壓，Pa；R——通用氣體常數(shù)，J/（mol·K）；t——時間，s；T——天然氣溫度，K；Ts——關(guān)閥時間，s；u0——井口流速，m/s；um——氣液混合相流速，m/s；um,0——氣液混合相初始流速，m/s；um,1，um,2——水頭位置1、2處氣液混合相流速，m/s；umP1,i，umP2,i——i時刻點P1和P2處油管節(jié)點橫截面的流速，m/s；z，z1，z2——水頭位置，m；Z——天然氣壓縮因子，無因次；α1，α2——水頭位置1，2的動能修正系數(shù)，無因次；β——閥門開度系數(shù)，無因次；Δp——油管內(nèi)的壓力增量，Pa；Δt——時間步長，s；Δz——井深步長，m；ρm——氣液混合相密度，kg/m3；ρmol——摩爾密度，kmol/m3；τ——油管壁面剪切應(yīng)力，Pa；φ，φ1，φ2——油管軸線與水平位置的夾角，（°）；Ψ，Ψ1，Ψ2——油管側(cè)壁與軸線的夾角，（°）；ω——擬合待定系數(shù)，m2·s/kg。