基于狀態(tài)依賴時滯的鉆柱動力學(xué)穩(wěn)定性分析

張 鶴，狄勤豐,2，王文昌,2，陳 鋒，段浩宇

(1.上海大學(xué) 上海市應(yīng)用數(shù)學(xué)和力學(xué)研究所，上海 200072；2.上海市能源工程力學(xué)重點(diǎn)實(shí)驗(yàn)室，上海 200072；3.上海大學(xué) 機(jī)電工程與自動化學(xué)院，上海 200072)

鉆井工程是石油、天然氣、頁巖氣勘探開發(fā)的重要環(huán)節(jié)之一，主要為油、氣提供生產(chǎn)通道。鉆柱作為鉆井工程中的核心部件，由鉆桿和底部鉆具組合(bottom hole assembly,BHA)構(gòu)成。BHA主要由鉆頭、鉆鋌、穩(wěn)定器以及井下測量和輔助鉆井工具等組成，為一長度約為100～300 m的管柱結(jié)構(gòu)。鉆井過程中，鉆柱在地面設(shè)備驅(qū)動下，在充滿液體或氣體的狹長井眼中轉(zhuǎn)動，并帶動鉆頭旋轉(zhuǎn)破碎井底巖石，形成井筒(見圖1)。鉆頭通過切削齒擠壓或剪切的方式破壞巖石，與巖石的相互作用會激發(fā)鉆柱的振動，其主要有三種形式[1]：軸向振動[2]、扭轉(zhuǎn)振動[3]和橫向振動[4-5]，三種振動在井下相互耦合，其最劇烈、最具破壞性的極端表現(xiàn)形式為跳鉆、黏滑和渦動。

自二十世紀(jì)九十年代始，聚晶金剛石復(fù)合片(polycrystalline diamond compact，PDC)鉆頭得到廣泛應(yīng)用，但使用PDC鉆頭進(jìn)行深井鉆探時，鉆柱容易發(fā)生扭轉(zhuǎn)振動甚至惡化為黏滑振動，其主要表現(xiàn)為鉆頭處于停鉆(鉆頭轉(zhuǎn)速為0，稱為黏滯相)和高速轉(zhuǎn)動(鉆頭最大轉(zhuǎn)速達(dá)到地面恒定轉(zhuǎn)速的兩倍及以上，稱為滑脫相)的周期交替振蕩狀態(tài)。據(jù)統(tǒng)計(jì)[6-7]，在深井鉆探中，鉆柱發(fā)生黏滑振動的時長占鉆進(jìn)時間的一半以上。鉆柱黏滑振動引起的交變應(yīng)力容易造成鉆具接頭的疲勞失效，加劇鉆頭磨損，降低鉆井效率。因此，本文主要以鉆柱的扭轉(zhuǎn)黏滑振動為研究對象，通過對描述鉆柱黏滑振動的數(shù)學(xué)模型進(jìn)行穩(wěn)定性分析，并在此基礎(chǔ)上研究鉆井參數(shù)與鉆井液阻尼對鉆柱黏滑振動的影響。

圖1 鉆井和鉆柱系統(tǒng)示意圖(改自文獻(xiàn)[8]) Fig.1 Schematic of drilling and drillstring systems (Modified from [8])

早期對鉆柱黏滑振動的研究主要是對扭轉(zhuǎn)振動單獨(dú)建模求解，不考慮扭轉(zhuǎn)振動與軸向或橫向振動的相互耦合[9-10]，并認(rèn)為鉆頭反作用扭矩隨鉆頭轉(zhuǎn)速增加而減小的速度弱化效應(yīng)是造成鉆柱發(fā)生黏滑振動的主要原因，并提出不同的函數(shù)關(guān)系來表示鉆頭反扭矩與轉(zhuǎn)速之間的速度弱化規(guī)律[11-19]，在一定程度上重現(xiàn)了鉆柱的黏滑振動。然而，隨著對振動機(jī)理認(rèn)識的不斷深入，人們發(fā)現(xiàn)鉆柱不同振動模式之間存在相互耦合，需要考慮耦合效應(yīng)才能正確地理解鉆柱的黏滑振動。同時，單個鉆頭切削齒切削巖石的試驗(yàn)結(jié)果與速度弱化效應(yīng)并不相符，表明了鉆頭和巖石之間的作用力與鉆頭運(yùn)動速度無關(guān)[20]。

Detournay等[21]在單個鉆頭切削齒切削巖石試驗(yàn)結(jié)果的基礎(chǔ)上，提出了與鉆頭轉(zhuǎn)速無關(guān)的鉆頭-巖石相互作用模型，該模型將鉆頭與巖石相互作用分為齒切削面對巖石的切削過程以及齒磨損面與巖石的接觸摩擦過程。Richard等[22]在上述鉆頭-巖石相互作用模型的基礎(chǔ)上，考慮了BHA的軸向和扭轉(zhuǎn)運(yùn)動，建立了低維鉆柱動力學(xué)離散模型(文獻(xiàn)中稱為RGD模型)。RGD模型將鉆頭的切削深度表示為鉆頭在當(dāng)前時刻與之前某一時刻軸向位置的差值，進(jìn)而在模型中引入了取決于鉆頭的運(yùn)動狀態(tài)的時滯變量，使得鉆頭的軸向和扭轉(zhuǎn)振動通過該狀態(tài)依賴時滯相互耦合，其運(yùn)動控制方程為狀態(tài)依賴時滯微分方程。數(shù)值模擬結(jié)果表明，鉆頭切削巖石過程的狀態(tài)依賴時滯是造成鉆柱發(fā)生自激黏滑振動的根本原因，而鉆頭扭矩表現(xiàn)出的速度弱化效應(yīng)僅為鉆頭發(fā)生黏滑振動后的宏觀表現(xiàn)。由于RGD模型中的鉆柱模型和鉆頭的幾何形狀較為簡化，因此近年來很多國外學(xué)者對其進(jìn)行了改進(jìn)，其中主要包括：①在鉆柱模型中考慮鉆桿的軸向剛度和鉆井液阻尼[23-25]；②建立鉆柱的多自由度有限元模型或連續(xù)模型[26-29]；③研究幾何形狀更為復(fù)雜的鉆頭，其與巖石的相互作用涉及兩個甚至多個狀態(tài)依賴時滯[30-31]；④通過定義二元函數(shù)來表示時滯變量，避免向鉆柱模型中引入狀態(tài)依賴時滯[32-34]；但國內(nèi)在這方面的研究鮮有報道。

對非線性動力學(xué)系統(tǒng)進(jìn)行線性穩(wěn)定性分析，可以定性地研究系統(tǒng)參數(shù)對系統(tǒng)穩(wěn)態(tài)運(yùn)動的影響。文獻(xiàn)中對RGD及其改進(jìn)模型的穩(wěn)定性分析主要分為兩種：一種是直接對狀態(tài)依賴時滯微分方程進(jìn)行線性化，得到線性系統(tǒng)的特征方程，進(jìn)而得到特征方程根軌跡的解析解[35]；另一種是將狀態(tài)依賴時滯微分方程轉(zhuǎn)化為非線性耦合的偏微分-常微分方程組，并對其進(jìn)行線性化求解[36]。針對RGD模型的線性穩(wěn)定性分析指出，由于忽略阻尼的影響，因此其在任何參數(shù)下均不穩(wěn)定。然而，以上解析方法僅適用于低維鉆柱模型，難以應(yīng)用到多自由度鉆柱動力學(xué)系統(tǒng)的穩(wěn)定性分析。鑒于此，本文考慮阻尼的影響，利用Zhang等研究中的鉆頭軌跡函數(shù)，將RGD模型的狀態(tài)依賴時滯微分方程轉(zhuǎn)化為非線性耦合的偏微分和常微分方程。在此基礎(chǔ)上，利用譜方法和伽遼金方法將線性化后的偏微分和常微分方程離散為耦合的常微分方程組，以實(shí)現(xiàn)對RGD模型的穩(wěn)定性分析。此外，本文提出的穩(wěn)定性數(shù)值分析方法還容易推廣到帶狀態(tài)依賴時滯的多自由度鉆柱動力學(xué)系統(tǒng)的穩(wěn)定性分析。

1 數(shù)學(xué)模型

1.1 鉆柱模型

RGD模型在軸向上，僅考慮了BHA的集中質(zhì)量Mb，忽略了鉆桿的軸向剛度和軸向阻尼，這是因?yàn)閷?shí)際鉆井中鉆柱在軸向上、下邊界均為力邊界條件(上邊界為大鉤載荷H0，下邊界為鉆頭-巖石相互作用的鉆壓Wb)，因此鉆桿軸向剛度和軸向阻尼對BHA的軸向運(yùn)動幾乎沒有影響(見圖2)；而鉆柱扭轉(zhuǎn)方向的上、下邊界分別為位移邊界和力邊界(上邊界為地面恒定轉(zhuǎn)速Ω0，下邊界為鉆頭和巖石相互作用的扭矩Tb)，因此在扭轉(zhuǎn)方向?qū)HA簡化為轉(zhuǎn)動慣量為Jb的剛體，將鉆桿簡化為扭轉(zhuǎn)剛度為Kp的扭轉(zhuǎn)彈簧，忽略了扭轉(zhuǎn)阻尼的影響，從而使得RGD模型在任何參數(shù)下均不穩(wěn)定。

基于此，本文在RGD模型的基礎(chǔ)上考慮了扭轉(zhuǎn)方向的阻尼以研究其對RGD模型穩(wěn)定性的影響。由圖2可知，RGD模型考慮了BHA的軸向和扭轉(zhuǎn)兩個自由度，其軸向和扭轉(zhuǎn)位移(速度)分別用Ub(Vb)和Φb(Ωb)表示。根據(jù)鉆柱的上、下邊界條件，可以推導(dǎo)出控制BHA軸向和扭轉(zhuǎn)運(yùn)動的常微分方程(ordinary differential equation,ODE)分別為

(1)

(2)

圖2 鉆柱動力學(xué)模型示意圖 Fig.2 Schematic diagram of drill string dynamics model

1.2 鉆頭-巖石相互作用模型

鉆頭處的鉆壓Wb和扭矩Tb可以分別由兩個分量構(gòu)成[37]，即

Wb=Wc+Wf,Tb=Tc+Tf

(3)

式中：Wc和Wf分別為鉆壓的切削和接觸摩擦分量；Tc和Tf分別為扭矩的切削和接觸摩擦分量。

切削分量Wc和Tc的大小均與鉆頭的切削深度d成正比，即

(4)

式中：a為鉆頭半徑；ε為巖石的本征比能；ζ描述了切削力在切削平面上的方位。鉆頭的準(zhǔn)螺旋運(yùn)動在井底形成一定的巖石輪廓，而切削深度取決于相鄰刀翼間的井底巖石輪廓，如圖3所示。

圖3 相鄰刀翼間井底巖石輪廓及切削深度示意圖(改自Richard等的研究) Fig.3 Section of the bottom hole profile between two adjacent blades (Modified from Richard,et al)

Richard 等研究中通過引入軸向再生效應(yīng)來表示鉆頭的切削深度，其大小為鉆頭軸向位移在當(dāng)前時刻與此前某一時刻的差值

d=ndn,dn=Ub(t)-Ub(t-tn)

(5)

式中：n為鉆頭的刀翼個數(shù)；dn為單個刀翼的切削深度，由于鉆頭做剛體運(yùn)動且刀翼沿鉆頭軸向?qū)ΨQ均勻分布，因此每個刀翼的切削深度均相同；ub(t)為鉆頭當(dāng)前時刻的軸向位移，Ub(t-tn)為鉆頭在t-tn時刻的軸向位移，時滯時間tn為鉆頭繞其軸線轉(zhuǎn)過2π/n(相鄰刀翼間的夾角)所需的時間，即

(6)

式中，Φb(t)和Φb(t-tn)分別為鉆頭在當(dāng)前時刻和t-tn時刻的扭轉(zhuǎn)角位移。由于鉆頭存在扭轉(zhuǎn)振動，因此時滯tn不為常數(shù)，其大小取決于鉆頭的運(yùn)動狀態(tài)，也即tn為狀態(tài)依賴時滯。

鉆壓和扭矩的接觸摩擦分量Wf和Tf之間滿足以下約束關(guān)系

(7)

式中：μ為鉆頭切削齒磨損面與巖石間的摩擦因數(shù)；γ為與鉆頭幾何形狀相關(guān)的常量，對于RGD模型中所用鉆頭，其取值為1。Wf的取值與鉆頭的軸向速度Vb相關(guān)

Wf=anlσ[1-g(Vb)]

(8)

式中：σ為鉆頭與巖石平均接觸應(yīng)力的最大值；g(Vb)為集值函數(shù)，其定義為

(9)

由于鉆頭切削過程中引入了狀態(tài)依賴時滯tn，因此將式(3)～式(9)描述的鉆頭-巖石相互作用代入微分方程式(1)～式(2)中，可以得到控制鉆頭軸向和扭轉(zhuǎn)運(yùn)動的狀態(tài)依賴時滯微分方程(state-dependent delay differential equation,SDDDE)；另一方面，鉆頭與巖石的接觸摩擦過程向方程中引入了集值函數(shù)g(Vb)，從而使鉆頭處的邊界條件呈現(xiàn)出極強(qiáng)的非線性。

1.3 鉆頭軌跡函數(shù)

為了描述鉆頭的空間運(yùn)動軌跡，需要建立整體和局部兩個坐標(biāo)系，如圖4所示。整體坐標(biāo)系ORXΘ的原點(diǎn)O固定在井口；R為原點(diǎn)指向井壁的徑向坐標(biāo)；X沿井眼軸線方向，用以表示鉆頭的軸向位置；Θ正方向與鉆頭扭轉(zhuǎn)運(yùn)動正方向一致，用以表示鉆頭的扭轉(zhuǎn)角位移。局部坐標(biāo)系orxθ的原點(diǎn)固定在鉆頭上隨鉆頭一起轉(zhuǎn)動；r沿鉆頭半徑方向，用以表示鉆頭切削齒在半徑方向上的分布；x指向鉆頭鉆進(jìn)方向，表示鉆頭切削齒在鉆頭上的軸向位置；θ正方向與鉆頭轉(zhuǎn)動正方向一致，表示鉆頭相對于整體坐標(biāo)系的轉(zhuǎn)動。RGD模型假設(shè)刀翼沿鉆頭軸線對稱均勻分布，且鉆頭切削齒沿半徑方向緊密排列，因此0≤r≤a；且鉆頭切削齒均在同一水平面上。若令參考刀翼的坐標(biāo)為x=0,0≤r≤a,θ=0，則沿θ正向第i個刀翼的坐標(biāo)可表示為x=0,0≤r≤a,θ=2π(i-1)/n。

圖4 整體和局部坐標(biāo)系示意圖 Fig.4 Schematic of global and local cylindrical coordinate systems

由于鉆頭做剛體運(yùn)動且刀翼沿鉆頭軸線對稱均勻分布，因此鉆頭軌跡在相鄰刀翼間形成的井底輪廓曲線均相同，如圖5所示(4刀翼鉆頭運(yùn)動軌跡)，其中刀翼的軌跡函數(shù)由方程X=H[Φb(t)-θ]表示。在此基礎(chǔ)上定義鉆頭軌跡函數(shù)

h(θ,t)=-H[Φb(t)-θ],θ∈[0,2π/n],t>0

(10)

根據(jù)圖5和式(5)可知，單個刀翼的切削深度dn可用鉆頭軌跡函數(shù)h(θ,t)表示為

(11)

此外，在整體坐標(biāo)系ORXΘ中，刀翼的軌跡函數(shù)X=H(Θ),0<Θ<Φb(t)與時間t無關(guān)，因此在給定坐標(biāo)Θ=Φb(t)-θ處

(12)

(13)

由式(11)可知，通過引入鉆頭軌跡函數(shù)表示鉆頭切削深度，避免了向方程中引入狀態(tài)依賴時滯。此外，鉆頭軌跡函數(shù)的變化受式(13)控制，將鉆頭切削深度式(11)代入鉆壓和扭矩切削分量表達(dá)式(3)～式(4)并進(jìn)一步代入到式(1)～式(2)后，即可將SDDDE轉(zhuǎn)化為耦合的PDE-ODE。

圖5 井底巖石輪廓的平面展開示意圖 Fig.5 Planar schematic of bottom hole profile

1.4 無量綱化

首先定義系統(tǒng)的特征時間和特征長度

(14)

(15)

式中，U0和Φ0分別為鉆頭穩(wěn)態(tài)鉆進(jìn)時的軸向和扭轉(zhuǎn)位移。無量綱形式的切削深度δ0,δn和時間τ可分別表示為

(16)

由此，式(1)和式(2)的無量綱形式為

(17)

其中，

(18)

(19)

式中，v0為鉆頭穩(wěn)態(tài)鉆進(jìn)時的軸向速度的無量綱形式。

(20)

(21)

?(0,θ)=-ub(τ)

(22)

式中，ω0為鉆頭穩(wěn)態(tài)運(yùn)動時扭轉(zhuǎn)速度的無量綱形式。在此基礎(chǔ)上，鉆頭切削深度表達(dá)式(11)可以用無量綱形式的鉆頭軌跡函數(shù)表示為

(23)

2 線性化及穩(wěn)定性分析

2.1 線性化

在非線性耦合PDE-ODE穩(wěn)態(tài)解的基礎(chǔ)上施加小擾動，可以實(shí)現(xiàn)對該方程組的線性化處理。由式(15)及偏微分方程初始條件式(21)可知PDE-ODE方程的穩(wěn)態(tài)解為

(24)

(25)

式中的擾動均為小量。

當(dāng)鉆頭穩(wěn)態(tài)運(yùn)動時，g(vb)=0，因此將擾動表達(dá)式(25)代入式(17)～式(18)、式(20)～式(23)，忽略擾動小量的乘積并整理后，可得線性化后控制擾動變量的耦合PDE-ODE為

(26)

(27)

其中偏微分方程式(27)的初始和邊界條件為

(28)

2.2 穩(wěn)定性分析

利用Gupta等研究中的方法，可以推導(dǎo)出上述線性PDE-ODE方程的解析解，然而該解析法只適用于低維的鉆柱模型(僅含有兩個自由度)，無法推廣到多自由度鉆柱模型的穩(wěn)定性分析。為此，本文提出了分析線性PDE-ODE穩(wěn)定性的數(shù)值方法，該方法不受系統(tǒng)自由度的限制，可以容易地拓展到多自由度系統(tǒng)的穩(wěn)定性分析。

(29)

(30)

(31)

(32)

(33)

將式(30)中的第二式代入常微分方程式(26)后可得

(34)

(35)

則上述耦合的常微分方程組寫成矩陣矢量的表示形式為

(36)

式中，系數(shù)矩陣A和B中的元素均是由無量綱參數(shù)組成的已知量。由此，利用譜方法和伽遼金方法可以將線性PDE-ODE的穩(wěn)定性分析轉(zhuǎn)化為分析線性常微分方程式(36)的穩(wěn)定性。線性常微分方程式(36)的穩(wěn)定性可以通過其特征值實(shí)部的正負(fù)來判斷，當(dāng)所有特征值的實(shí)部均為負(fù)時，線性常微分方程式(36)的零解漸進(jìn)穩(wěn)定，則原非線性PDE-ODE方程的穩(wěn)態(tài)解也漸進(jìn)穩(wěn)定；相反，若至少有一個特征值的實(shí)部為正，則線性常微分方程式(36)的零解不穩(wěn)定，同樣原非線性PDE-ODE方程的穩(wěn)態(tài)解也不穩(wěn)定；若除含有零實(shí)部特征值外，其余特征值的實(shí)部均為負(fù)，則線性方程組的零解及原非線性PDE-ODE方程的穩(wěn)態(tài)解均為臨界穩(wěn)定。

3 穩(wěn)定性結(jié)果驗(yàn)證及阻尼影響

3.1 穩(wěn)定性結(jié)果驗(yàn)證

RGD模型存在兩種不穩(wěn)定性機(jī)制：快不穩(wěn)定性機(jī)制和慢不穩(wěn)定性機(jī)制。當(dāng)?shù)孛孓D(zhuǎn)速ω0大于臨界轉(zhuǎn)速ωc時，系統(tǒng)的穩(wěn)定性由扭轉(zhuǎn)振動的不穩(wěn)定極點(diǎn)主導(dǎo)，表現(xiàn)為慢不穩(wěn)定性，鉆頭軸向和扭轉(zhuǎn)速度為準(zhǔn)周期運(yùn)動，需要經(jīng)過很長時間才能發(fā)展為自激的黏滑振動；而當(dāng)ω0<ωc時，系統(tǒng)的不穩(wěn)定性由軸向振動的不穩(wěn)定極點(diǎn)主導(dǎo)，表現(xiàn)為快不穩(wěn)定性，即鉆頭在短時間內(nèi)即出現(xiàn)自激的黏滑振動。臨界轉(zhuǎn)速ωc的近似計(jì)算公式為

(37)

圖6 RGD模型穩(wěn)定性圖譜 Fig.6 Stability map of the RGD model

3.2 阻尼影響

本文模型考慮了鉆井液的扭轉(zhuǎn)阻尼，其對系統(tǒng)穩(wěn)定性的影響見圖8，其中κ=0.01，而其余參數(shù)保持不變。圖8中灰色表示所有特征根的實(shí)部均為負(fù)值，即系統(tǒng)在該參數(shù)下穩(wěn)定。將圖8與圖6對比可知，考慮鉆井液阻尼后，RGD模型的快不穩(wěn)定區(qū)域沒有發(fā)生變化，僅慢不穩(wěn)定區(qū)域變?yōu)榉€(wěn)定性區(qū)域。

圖7 RGD模型數(shù)值模擬結(jié)果 Fig.7 Simulation results of RGD model

圖8 考慮阻尼后RGD模型穩(wěn)定性圖譜 Fig.8 Stability map of the RGD model with considering damping

圖9 考慮阻尼κ=0.01的模擬結(jié)果 Fig.9 Simulation results with damping κ=0.01

4 結(jié) 論

(1)利用鉆頭軌跡函數(shù)對鉆頭切削深度進(jìn)行重新表述，可避免向鉆柱動力學(xué)系統(tǒng)中引入狀態(tài)依賴時滯變量，從而將傳統(tǒng)的SDDDE轉(zhuǎn)化為非線性耦合的PDE-ODE方程。

(2)在鉆頭穩(wěn)態(tài)運(yùn)動的基礎(chǔ)上施加小擾動可以使非線性耦合的PDE-ODE線性化，從而利用譜方法和伽遼金方法，將線性化后的PDE-ODE離散為控制擾動的線性微分方程組，通過計(jì)算該方程組的特征值可以實(shí)現(xiàn)對系統(tǒng)穩(wěn)定性的分析。

(3)利用文獻(xiàn)中的結(jié)果和數(shù)值模擬結(jié)果可以證明本文數(shù)值方法的正確性，此外本文方法還可以預(yù)測鉆壓對系統(tǒng)穩(wěn)定性的影響。

(4)在RGD模型中考慮阻尼會改變RGD模型的穩(wěn)定性，但不改變模型的穩(wěn)定性邊界。

(5)本文提出的數(shù)值方法可擴(kuò)展于分析帶狀態(tài)依賴時滯的多自由系統(tǒng)的穩(wěn)定性。