李 中，王國榮 ，方達科，魏安超，柳軍

1.中海石油（中國）有限公司湛江分公司，廣東 湛江 524057；2.西南石油大學(xué)機電工程學(xué)院，四川 成都610500

引言

中國高溫高壓天然氣資源豐富，在高產(chǎn)氣井生產(chǎn)過程中，油管柱（圖1）受高速流體的影響，極易發(fā)生縱橫向耦合非線性振動和油套管非線性接觸碰撞。這將導(dǎo)致管柱局部產(chǎn)生交變荷載和動態(tài)應(yīng)力，管柱發(fā)生疲勞損壞；管柱上端受到較大的拉應(yīng)力，管柱發(fā)生強度破壞；同時將增大管柱軸向壓力，使其發(fā)生屈曲變形，增加油管與套管之間的摩擦磨損量，發(fā)生摩擦穿孔破壞[1]（圖2）。因此，亟需開展高產(chǎn)氣井油管柱非線性振動模型研究，為中國高產(chǎn)氣井的安全開發(fā)奠定理論基礎(chǔ)。

圖1 高產(chǎn)氣井管柱結(jié)構(gòu)圖Fig.1 Tubing string structure of high-yield gas well

圖2 管柱振動引起的失效問題Fig.2 Failure problems caused by vibration of tubing string

國內(nèi)外學(xué)者針對考慮內(nèi)流引起的管柱振動問題，開展了大量研究，并取得了一些成果。早期學(xué)者針對內(nèi)流作用下的管柱振動問題開展了初步研究[2]，通過實驗驗證了內(nèi)流對管柱的動力響應(yīng)有顯著影響[3]。隨后，諸多學(xué)者針對管柱振動模型開展了詳細研究，建立了流體作用力的計算方法[4]、管柱縱向振動[5]、橫向振動[6]以及流固耦合振動模型[7]，在此基礎(chǔ)上，分析了管柱的振動響應(yīng)特性。近年來，部分學(xué)者[8-10]發(fā)現(xiàn)細長管柱的縱橫向耦合效應(yīng)明顯，不可忽略，邢譽峰等[11]針對梁的縱橫向振動問題，建立了梁的縱橫向耦合振動方程，陽明君等[12]針對尤拉屯氣田完井管柱的振動問題，采用有限元方法ANSYS 建立高產(chǎn)氣井大長徑比管柱振動模型，Liu 等[13]考慮管柱的幾何非線性問題，建立了海洋立管的縱橫向耦合振動模型?？梢姡瑢τ诩氶L結(jié)構(gòu)的振動分析，不能忽略其幾何非線性因素（即縱橫向耦合振動），高產(chǎn)氣井油管柱屬于典型的細長結(jié)構(gòu)（長徑比約30 000）。

高產(chǎn)油氣井油管柱一旦發(fā)生縱橫向耦合振動，極易使其產(chǎn)生非線性接觸碰撞。目前，針對油套管非線性接觸碰撞問題，研究主要集中在鉆柱或管柱發(fā)生屈曲后的接觸問題[14-16]，得到管柱屈曲后接觸力的靜態(tài)計算方法，這與生產(chǎn)工況下管柱流致振動引起的動態(tài)接觸碰撞具有一定的差別。國內(nèi)外學(xué)者針對流致振動引起的接觸碰撞問題，主要建立了換熱管束與支承板接觸碰撞力的計算方法，Rogers 等[17]和Ting 等[18]采用動態(tài)有限元法預(yù)測了單個換熱器與支承板之間的接觸碰撞力，將預(yù)測結(jié)果與實驗數(shù)據(jù)進行了對比，驗證了預(yù)測方法的正確性，但未給出具體的計算方法。丁傳義等[19]和沈時芳等[20]開展了換熱器傳熱管與支承板之間接觸碰撞力的計算方法研究，得到單一位置接觸的碰撞力計算方法和碰撞力的時程響應(yīng)，并與實驗數(shù)據(jù)對比，驗證了計算方法的正確性，但模型只考慮一點位置的接觸碰撞，對于多位置的接觸碰撞未開展相應(yīng)的研究。

因此，本文針對高產(chǎn)氣井高速流體誘發(fā)油管柱振動引起的破壞問題，開展油管柱非線性流致振動模型研究和振動模擬實驗研究，為中國高產(chǎn)氣井油管柱的安全設(shè)計提供理論分析。

1 高產(chǎn)氣井油管柱雙重非線性流致振動模型建立及求解

1.1 油管柱縱橫向耦合非線性流致振動模型

管內(nèi)的高速流體是引起油管柱振動的主要因素，由于油管柱屬于細長結(jié)構(gòu)，易產(chǎn)生縱向和橫向振動，縱向和橫向振動又不是單獨分開運動，它們之間相互影響，因此，需建立油管柱的縱橫向耦合動力學(xué)模型。

本文建立以深度方向為x軸，水平向右為y 軸的平面坐標系（圖3）。

圖3 油管柱物理模型Fig.3 Physical model of tubing string

把油管簡化為均勻的Rayleigh梁，若考慮縱橫耦合，其幾何關(guān)系為[11]

式中：εij—6 個應(yīng)變分量，i,j=x,y；

y—厚度坐標，m；

u1(x,y,t)，u2(x,t)—與坐標系（x，y）對應(yīng)的位移場函數(shù)，m。

其表達式為

式中：u(x,t)—縱向位移，m；

w(x,t)—橫向位移，m；

t—時間，s。

式中：

T—油管的動能，J；

U—油管的應(yīng)變能，J；

W—外力做功，J；

L—管長，m；

ρ，ρ0—油管密度和流體密度，kg/m3；

A—油管橫截面積，m2；

I—油管截面的慣性矩，m4；

E—油管彈性模量，Pa；

f(x,t)，p(x,t)—油管受到的縱向和橫向外激力函數(shù)，N。

根據(jù)哈密頓（Hamilton）變分原理，建立油管的振動微分方程

式中：

δ—油套管接觸后產(chǎn)生的形變，m。

由式（3）～式（6）化簡得油管的縱向、橫向振動微分方程

由于上端為油管掛，下端為封隔器，把上下端視為固定端，即邊界條件和初始條件為

1.2 油管與套管非線性接觸碰撞模型

在內(nèi)部高速流體的作用下，油管柱會產(chǎn)生橫向和縱向振動，當橫向振動位移達到套管與油管之間的間隙時，油管與套管發(fā)生接觸碰撞，縱向產(chǎn)生接觸摩擦，造成油管表面摩擦磨損。因此，在縱橫向耦合振動下，考慮油管與套管的接觸碰撞效應(yīng)是進行油管疲勞磨損分析的必要條件。

1.2.1 套管與油管碰撞力-變形非線性關(guān)系

根據(jù)彈塑性力學(xué)理論，建立油管碰撞力與形變之間的關(guān)系。變形結(jié)構(gòu)如圖4 所示，油管發(fā)生碰撞，油管上的A2點變形到了套管上的A1點。

圖4 油管與套管接觸變形示意圖Fig.4 Schematic diagram of contact deformation betweentubing-casing

由幾何關(guān)系得

式中：

R1，R2—套管、油管的半徑，m；

z1—套管接觸點到套管內(nèi)壁的徑向距離，m；

z2—油管接觸點到套管內(nèi)壁的徑向距離，m；

r—接觸點到油管軸心的水平距離，m。

在力的作用下，套管和油管發(fā)生變形，它們之間的形變?yōu)棣?，并產(chǎn)生了寬為2b的接觸帶。由幾何關(guān)系可得

式中：

ω1—套管接觸后軸向產(chǎn)生的位移，m；

ω2—油管接觸后軸向產(chǎn)生的位移，m。

若接觸帶的寬度比油管的半徑小得多，則每個管柱都可以近似地當作彈性半平面來考慮，由文獻[21]可得ω1和ω2的計算公式

式中：μ1，μ2—套管和油管材料泊松比，無因次；

E1，E2—套管和油管材料彈性模量，Pa；

2b—油套管接觸帶寬度，m；

q(x)—套管和油管均布載荷，N；

p—油套管接觸力，N。

將所得的ω1，ω2代入變形條件

假定q(x)是與以2b為直徑所作的半圓弧縱坐標成比例，得

將E1=E2=E、μ1=μ2=0.3 代入式（12）和式（13）中，得油管碰撞力與變形的關(guān)系式

1.2.2 套管與油管動態(tài)碰撞力與摩擦力的計算

套管與油管接觸碰撞問題可以簡化為多節(jié)點發(fā)生接觸碰撞問題，因此，研究套管和油管某一個點接觸，得到套管對油管碰撞力的計算方法，然后把計算方法應(yīng)用到管柱的每個節(jié)點。

套管與油管發(fā)生碰撞后如圖5 所示。

圖5 油管與套管碰撞示意圖Fig.5 Schematic diagram of tubing-casing collision

假設(shè)碰撞時油管受到“彈簧-阻尼器”的作用，因此，油管與套管發(fā)生碰撞，套管將給它施加一個碰撞力和碰撞阻尼力，碰撞力根據(jù)式（16）得到，阻尼力計算公式為[20]

式中：

pc—接觸阻尼力，N；

c—碰撞阻尼系數(shù)，無因次，c=1.5ap；

a—常數(shù)，無因次，對于鋼材，a=0.2～0.3[22]；

碰撞阻尼力也是變形的非線性函數(shù)。由此可得氣井油管與套管非線性碰撞力和摩檫力計算公式

式中：

F—油套管接觸碰撞力，N；

f—油套管摩檫力，N；

μ—油套管摩擦系數(shù)，一般取0.3[23]。

1.3 雙重非線性模型的求解

1.3.1 雙重非線性模型離散

本文采用線性拉格朗日函數(shù)和三次埃爾米特函數(shù)表達管柱的縱向位移u和橫向位移w，根據(jù)文獻[13]可知有限元離散形式為

式中：

l—單元的長度，m；

d—管柱單元位移矩陣；

ψ—管柱單元縱向位移形函數(shù)矩陣；

φ—管柱單元橫向位移形函數(shù)矩陣。

把位移函數(shù)式（20）和式（21）代入能量泛函，可以得到用結(jié)點位移向量表示的應(yīng)變能函數(shù)U、動能函數(shù)T的標準形式，即

式中：

w—油管橫向位移矩陣。

把結(jié)構(gòu)單元組裝后，根據(jù)變分原理可得系統(tǒng)的離散形式動力學(xué)方程

式中：

D，M，C，K，F(xiàn)—系統(tǒng)的位移矩陣、質(zhì)量矩陣、阻尼矩陣、剛度矩陣和載荷列向量，具體元素見式（20）和式（21）。

1.3.2 Newmark-β 求解方程

本文采用Newmark-β 逐步積分的求解方法，避免了任何疊加的應(yīng)用，能很好地適應(yīng)非線性模型。

Newmark-β 法假定

式中：

β，γ—按積分的精度和穩(wěn)定性要求進行調(diào)整的參數(shù)。本文取β=0.5，γ=0.25，即假定從t到t+?t時刻的速度不變，由式（25）和式（26）可得

考慮t+?t時刻的振動微分方程為

因此，可得模型的求解流程如圖6 所示，并采用Fortran 語言編寫了相應(yīng)的數(shù)值計算程序。

圖6 計算流程圖Fig.6 Computing flow chart

2 高產(chǎn)氣井油管柱雙重非線性流致振動模型驗證

2.1 油套管非線性接觸碰撞模型驗證

本文采用文獻[20]相同參數(shù)（表1），具體結(jié)構(gòu)如圖7 所示，通過設(shè)置4 處外激勵載荷，中間位置兩邊設(shè)置彈簧阻尼，采用本文建立的雙重非線性模型，計算管柱的碰撞力，并與文獻的實驗結(jié)果和文獻模型計算結(jié)果對比，驗證本文雙重非線性碰撞模型的正確性和有效性，文獻模型與本文模型的區(qū)別主要體現(xiàn)在文獻只考慮單一接觸碰撞非線性影響，而本文模型考慮了縱橫向耦合和接觸碰撞雙重非線性的影響。

圖7 管柱接觸碰撞結(jié)構(gòu)圖Fig.7 Contact collision structure diagram of pipe string

表1 計算模型參數(shù)Tab.1 Calculating model parameters

圖8 給出了單非線性模型、雙重非線性模型和實驗結(jié)果的接觸力時程響應(yīng)曲線。從圖中可以看出，本文提出的雙重非線性模型計算結(jié)果在振幅和變化規(guī)律上都比文獻[20]中的單非線性模型計算結(jié)果更接近實驗結(jié)果。雙重非線性模型比單非線性模型更能反映系統(tǒng)的高頻響應(yīng)特性，驗證了本文雙重非線性模型的正確性和高效性。

圖8 接觸碰撞力計算結(jié)果圖Fig.8 Calculation results of contact/collision force

2.2 油管柱雙重非線性流致振動模型實驗驗證

2.2.1 模擬實驗參數(shù)設(shè)計

由于現(xiàn)場無法測得有效的振動數(shù)據(jù)驗證本文模型的正確性，故借助模擬實驗方法驗證模型的正確性。管柱振動模擬實驗應(yīng)滿足3 個標準：幾何相似性、運動相似性和動力相似性[24-25]。本文采用幾何相似法設(shè)計了模擬實驗中油管和套管的基本尺寸（內(nèi)徑、外徑、管長等），由于長度方向與徑向尺寸差異較大，故不采用相同比例，徑向和縱向的相似比分別設(shè)置為5.0 和438.0。根據(jù)黃濤[26]和李子豐等[27]的研究，實驗管柱的材料密度和彈性模量與其對應(yīng)的實際管柱應(yīng)滿足

式中：ρp—實際管柱的密度，kg/m3；

ρm—模擬實驗管柱的密度，kg/m3；

Ep—實際管柱的彈性模量，Pa；

Em—模擬實驗管柱的彈性模量，Pa；

λ—相似比。

將實際管柱的密度ρp=7 850 kg/m3和彈性模量Ep=210 GPa 代入式（30），得

根據(jù)彈性模量和密度的比值，通過查材料手冊，選擇PE 管比較符合要求，其彈性模量和密度分別是Em=6.0 GPa 和ρm=1 200 kg/m3。

根據(jù)現(xiàn)場氣田的實際情況[地層壓力46.7 MPa，產(chǎn)量（100～200）×104m3/d，地層溫度150°C，地溫25°C，標準大氣壓0.1 MPa]，通過狀態(tài)方程計算，可以確定管內(nèi)類似流體速度。

式中：p1—現(xiàn)場工況的壓力，Pa；

V1—現(xiàn)場工況的體積，m3；

T1—現(xiàn)場工況的溫度，K；

p2—模擬實驗工況的壓力，Pa；

V2—模擬實驗工況的體積，m3；

T2—模擬實驗工況的溫度，K；

v′—模擬實驗的流速，m/s；

s—模擬管內(nèi)橫截面積，m2。

通過以上分析計算，得到本模擬實驗參數(shù)，如表2所示。

表2 模擬實驗參數(shù)Tab.2 Simulated test parameters

2.2.2 模擬實驗系統(tǒng)設(shè)計

本文實驗系統(tǒng)主要由注氣子系統(tǒng)、管柱子系統(tǒng)和連接子系統(tǒng)組成（圖9），具體包括氣源（螺桿空氣壓縮機、儲氣罐）、管道、電磁閥、氣體流量計、油管柱模型、套管模型、固定裝置、軸力裝置（重量、滑輪）、消聲器等，設(shè)計了實驗臺架示意圖如圖10 所示，并搭建了相應(yīng)的實驗臺架（圖11）。

圖9 實驗系統(tǒng)設(shè)計流程圖Fig.9 Flow chart of experimental system design

圖10 實驗臺架示意圖Fig.10 Schematic diagram of test bench

圖11 實驗系統(tǒng)實物圖Fig.11 Physical diagram of the experimental system

為了采集油管柱的變形特征，在其四周（CF1，CF2，IL1，IL2）分別布置8 個采集點，共計32 個點位，每個點位上橫向布置兩個應(yīng)變片，用于溫度補償，消除實驗過程中溫度、初始變形等誤差。實驗油管柱采集點1 和采集點8 距兩端0.15 m，1～8 號采集點相鄰之間間距1.00 m（如圖12 所示）。

圖12 應(yīng)變片安裝示意圖Fig.12 Diagram of strain gauge installation

2.2.3 實驗結(jié)果分析

采用本文建立的雙重非線性流致振動模型，設(shè)置與實驗參數(shù)一樣（表2），計算模型如圖13 所示，其中管柱分為300 個單元，模擬總時間為70 s，步長為0.000 1 s，提取了4 個測量點（與實驗中安裝傳感器的位置相同）的振動響應(yīng)。

圖13 雙重非線性流致模型的計算簡圖Fig.13 Computational sketch of a bi-nonlinear model

圖14 為油管不同位置橫向振動位移時程曲線，可以看出，通過雙重非線性流致振動模型計算和實驗測量所得到油管柱橫向振動幅值基本一致，兩種方法都表明位移響應(yīng)存在50 s 的瞬態(tài)響應(yīng)，其中，實驗結(jié)果的高頻分量相對較多，主要由于實驗環(huán)境因素的干擾。圖15 為油管不同位置橫向振動幅頻響應(yīng)曲線圖，由圖可知，雙重非線性流致振動模型計算的主頻率與實驗測量結(jié)果基本一致，分別為1.4 Hz 和1.6 Hz，管柱在0～2.0 Hz 的振動能量較大。通過管柱不同測點的時域分析和頻域分析，驗證了本文所建立的雙重非線性流致振動模型的正確性，為中國高產(chǎn)氣井管柱安全設(shè)計提供了有效的分析工具。

圖14 油管不同位置橫向振動位移時程曲線Fig.14 Time history curves of transverse displacement at different positions of tubing string

圖15 油管不同位置橫向振動幅頻響應(yīng)曲線Fig.15 Amplitude-frequency response curve of transverse vibration at different positions of tubing string

3 結(jié)論

（1）采用能量法、微元法結(jié)合哈密頓原理建立了油管柱縱橫向耦合非線性流致振動模型，基于彈塑性接觸碰撞理論，建立了油管與套管接觸碰撞力和摩擦力計算方法，將其引入油管柱縱橫向耦合振動模型中，得到高產(chǎn)氣井油管柱雙重非線性流致振動模型。

（2）用有限元法和Newmark?β 法求解油管柱的雙重非線性流致振動模型，選用文獻的實驗參數(shù)，計算得到油套管接觸碰撞力時程曲線，與文獻[20]的實驗數(shù)據(jù)和單非線性模型結(jié)果對比，本文建立的雙重非線性流致振動模型相對于單非線性模型更接近實驗數(shù)據(jù)，初步驗證了油管柱非線性流致模型的正確性。

（3）根據(jù)現(xiàn)場高產(chǎn)氣井參數(shù)，采用相似原理，設(shè)計了油管柱振動實驗方案，搭建了相應(yīng)的實驗臺架，開展了高產(chǎn)氣井油管柱振動模擬實驗，測得實驗數(shù)據(jù)，與理論模型計算結(jié)果對比，本文所建立的非線性流致振動模型計算結(jié)果和實驗測得數(shù)據(jù)基本一致，進一步驗證了模型的正確性。