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

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

(1.上海大學 上海市應用數(shù)學和力學研究所，上海 200072；2.上海市能源工程力學重點實驗室，上海 200072；3.上海大學 機電工程與自動化學院，上海 200072)

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

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

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

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

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

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

1 數(shù)學模型

1.1 鉆柱模型

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

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

(1)

(2)

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

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

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

Wb=Wc+Wf,Tb=Tc+Tf

(3)

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

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

(4)

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

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

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

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

(5)

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

(6)

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

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

(7)

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

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

(8)

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

(9)

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

1.3 鉆頭軌跡函數(shù)

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

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

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

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

(10)

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

(11)

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

(12)

(13)

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

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

1.4 無量綱化

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

(14)

(15)

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

(16)

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

(17)

其中，

(18)

(19)

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

(20)

(21)

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

(22)

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

(23)

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

2.1 線性化

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

(24)

(25)

式中的擾動均為小量。

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

(26)

(27)

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

(28)

2.2 穩(wěn)定性分析

利用Gupta等研究中的方法，可以推導出上述線性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)定性分析轉化為分析線性常微分方程式(36)的穩(wěn)定性。線性常微分方程式(36)的穩(wěn)定性可以通過其特征值實部的正負來判斷，當所有特征值的實部均為負時，線性常微分方程式(36)的零解漸進穩(wěn)定，則原非線性PDE-ODE方程的穩(wěn)態(tài)解也漸進穩(wěn)定；相反，若至少有一個特征值的實部為正，則線性常微分方程式(36)的零解不穩(wěn)定，同樣原非線性PDE-ODE方程的穩(wěn)態(tài)解也不穩(wěn)定；若除含有零實部特征值外，其余特征值的實部均為負，則線性方程組的零解及原非線性PDE-ODE方程的穩(wěn)態(tài)解均為臨界穩(wěn)定。

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

3.1 穩(wěn)定性結果驗證

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

(37)

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

3.2 阻尼影響

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

圖7 RGD模型數(shù)值模擬結果 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的模擬結果 Fig.9 Simulation results with damping κ=0.01

4 結 論

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

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

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

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

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