張 強， 蔣 豹， 崔 巍， 劉巨保， 朱 昱

(東北石油大學 機械科學與工程院，黑龍江 大慶 163318)

受井筒約束的懸掛管柱，在其自身重力作用下，井底受壓，易喪失彈性穩(wěn)定性，導致屈曲。例如鉆柱屈曲會引起鉆頭方向改變，加大側向力和摩阻力，使鉆柱自鎖，更有甚者導致鉆具疲勞破壞。

Lubinski[1]首先研究了鉆柱在垂直井筒中的穩(wěn)定性，導出了鉆柱在垂直平面內的彎曲方程，并給出了該方程的級數(shù)解。在兩端鉸支邊界約束條件下，采用梁柱模型，無量綱長度取8，給出了鉆柱在垂直平面內發(fā)生一次彎曲和二次彎曲的臨界載荷，分別為1.94和3.75。但他在第2年發(fā)表的論文[2]中，將對應的臨界載荷改成2.04和4.05，在之后的40余年，我國石油工程行業(yè)一直沿用他第2年的數(shù)值。后來，韓志勇[3]回顧了Lubinski有關垂直井鉆柱屈曲理論在我國的應用情況，對一次彎曲和二次彎曲的臨界載荷進行了考證，其臨界載荷應該改回1.94和3.75。

Lubinski等[4]假定管柱屈曲成空間螺旋線，屈曲構型為均勻螺旋線并與井壁連續(xù)接觸，利用最小勢能原理，推導了等螺距和軸向壓力的關系式，他的這一研究成果得到了石油工程行業(yè)的廣泛應用。但是，由于自重的影響，實際管柱空間螺旋屈曲構型是一個非等距的螺旋線。為此，研究者們[5-8]采用微分法、能量法、實驗法和有限元等方法，對非等螺距的臨界載荷進行了廣泛的研究，得出了各自不同的臨界載荷值。無論是等螺距還是非等螺距假設，均是在螺旋屈曲之后形成一個螺距時得出的臨界載荷，其值明顯大于正弦屈曲之后剛螺旋屈曲時的臨界載荷。

由于管柱的螺旋屈曲發(fā)生在初始正弦屈曲之后，經(jīng)過穩(wěn)定的正弦屈曲構型階段之后的又一種穩(wěn)定的屈曲平衡狀態(tài)。所以，正弦屈曲向螺旋屈曲轉變的過程是非常復雜的。這一過程中，屈曲構型隨著井底壓力而發(fā)生跳躍性變化，管柱與井筒存在接觸和脫離等力學行為，使得問題非常復雜。

近年來，研究者們考慮瞬態(tài)效應的影響，研究了彈性壓桿的動力屈曲問題[9-11]。由于懸掛管柱在井下是緩慢上提或下放的，其動力屈曲并不明顯，研究者們普遍將此問題處理成靜力屈曲問題。

Gao等[12]展望了懸掛管柱螺旋屈曲問題的研究方法，提出了懸掛段采用梁柱模型，連續(xù)接觸段采用微分方程的研究設想。Huang等[7]根據(jù)這個研究設想，假設將懸掛段分成了4段，用相應的連續(xù)性條件、邊界條件和穩(wěn)定性條件，將管柱屈曲問題轉化為一個非線性方程組。用迭代法求解這些方程，得到了懸掛管柱平面正弦屈曲、三維螺旋屈曲、連續(xù)接觸螺旋屈曲和形成一個螺距的螺旋屈曲臨界載荷。

本文摒棄了預先指定的管柱屈曲撓曲線假設，考慮井筒內懸掛管柱的幾何和接觸雙重非線性特性，根據(jù)井筒內懸掛管柱屈曲分析的有限元平衡方程，采用漫動力法，進行有限元計算，分析懸掛管柱從正弦向螺旋轉變時的屈曲構型和臨界載荷。

1 模型建立

1.1 模型簡化

(1) 懸掛管柱在井筒內的初始狀態(tài)是完全豎直的。

(2) 將管柱和井筒處理成圓形等截面，忽略接頭等部件對屈曲的影響。

(3) 管柱與井筒之間有初始環(huán)空間隙存在，管柱變形前軸線與井筒軸線完全重合。

(4) 考慮管柱上端懸掛拉力和管柱自身的重力載荷作用，管柱上端受拉下端受壓。

(5) 管柱底部和頂部的邊界約束處理成鉸支約束。

(6) 管柱處于線彈性狀態(tài)，假設懸掛管柱為大柔度桿或細長桿。

(7) 管柱屈曲后與井筒內壁面接觸，考慮接觸力的影響，但不考慮管柱與井筒之間的摩擦力。

(8) 研究管柱靜力屈曲的最終穩(wěn)定狀態(tài)，忽略管柱上提下放過程中的瞬態(tài)影響。

1.2 力學模型

井筒內懸掛管柱的上端受拉，管柱的下端位于井底，在自重作用下受壓，建立的懸掛管柱力學模型如圖1所示。圖1(a)為懸掛管柱長度、載荷和邊界示意圖，鉸支約束邊界條件為上端約束x、y方向的平動自由度和繞z軸的轉動自由度，下端約束x、y和z方向的平動自由度。

懸掛管柱上端的懸掛拉力為

(1)

式中：E為懸掛管柱的彈性模量；I為懸掛管柱橫截面對軸的慣性矩；受拉段的無量綱長度ξT也稱為無量綱懸掛拉力。

懸掛管柱下端約束了z方向的平動自由度，其軸向約束反力為管柱受到的向上的壓力，其值為

(2)

式中：受壓段的無量綱長度ξC也稱為無量綱井底壓力。

(a) 載荷和邊界(b) 變形前橫截面(c) 變形后橫截面

圖1 懸掛管柱力學模型

Fig.1 Mechanical model of suspended tubular string

管柱變形前，由圖1(b)可計算出管柱初始間隙為δ=(D-d)/2，d為管柱外徑，D為井筒內徑。管柱在緩慢下放過程中，上端的懸掛拉力減小，下端的壓力增大。當管柱受壓段的壓力達到一定數(shù)值后，管柱喪失穩(wěn)定性而屈曲變形，產(chǎn)生橫向位移和較大的轉動變形，橫向變形后管柱的圓心為o′，見圖1(c)所示，設橫向位移量為(u,v)，對應的圓周角為θ。此時力與變形的關系不再是線性，非線性效應突出，屬于幾何非線性問題。

管柱在橫向彎曲變形過程中，受到井筒約束，將與井筒內壁面接觸，設接觸時的接觸力為Fn。在任一軸向距離、任一圓周方向上均可能產(chǎn)生接觸和脫離的力學行為，屬于接觸非線性問題。

綜上所述，懸掛管柱屈曲力學模型歸結為幾何和接觸雙重非線性問題，這樣處理使得懸掛管柱屈曲力學分析更趨于合理。

1.3 屈曲靜力學方程

采用有限元法，依據(jù)間隙元理論[13]，經(jīng)過梁單元和間隙元的拼裝，可得管柱幾何和接觸非線性靜力屈曲分析的總體平衡方程式為

(K0+Kσ(u)+Kn(u))u=Fq+Fn(u)

(3)

式中：K0為懸掛管柱的線彈性剛度矩陣；Kσ(u)為懸掛管柱的幾何剛度矩陣；Kn(u)為懸掛管柱的接觸剛度矩陣；u為懸掛管柱的節(jié)點位移向量，包括3個方向的平動和3個方向的轉動自由度，u=[u,v,w,θx,θy,θz]T；F為懸掛管柱的節(jié)點載荷向量，F(xiàn)=Fq+Fn(u)；Fq包含上端懸掛拉力和管柱自重均布載荷；Fn(u)為管柱的接觸力向量。式中的幾何剛度矩陣、接觸剛度矩陣和接觸力向量是節(jié)點位移的函數(shù)。

2 計算方法

2.1 慢動力法

采用有限元法，對式(3)進行雙重非線性靜力屈曲計算，但還存在收斂困難和算法不穩(wěn)定等問題。表現(xiàn)在如下方面：

(1) 由于懸掛管柱的從正弦屈曲到螺旋屈曲的構型存在跳躍性變化，管柱與井筒之間存在接觸和脫離等強非線性力學行為，導致計算過程中接觸狀態(tài)突變，易引起計算異常終止，存在收斂困難。還表現(xiàn)在懸掛管柱屈曲后，撓度突然增大，不能識別出接觸狀態(tài)，管柱穿透出井筒壁面，導致計算失敗。

(2) 算法不穩(wěn)定表現(xiàn)在收斂解不唯一，懸掛管柱屈曲的構型存在隨意性。例如，形成平面正弦屈曲構型后，隨著上端拉力的減小，下端壓力的增加，平面正弦屈曲構型在井筒內旋轉，卻不能形成空間螺旋屈曲構型；由于對管柱懸掛拉力載荷增量和載荷步長的敏感性，時而形成時而又不能形成空間螺旋屈曲構型；即使形成了空間螺旋屈曲構型，但又重新回到平面正弦屈曲構型；形成空間螺旋屈曲構型后，在不同軸向位置處，可能同時存在右螺旋和左螺旋的屈曲構型。

為了解決井筒內懸掛管柱從正弦屈曲向螺旋屈曲過程的收斂困難和算法不穩(wěn)定性問題，提出懸掛管柱非線性靜力屈曲分析的慢動力法。該方法是利用動力學方法，按照一定方式施加自重和懸掛拉力等恒定外載荷，引入慣性力和阻尼力，考慮時間積分效應，設置較大的虛擬阻尼，計算一定的時間，求解懸掛管柱動力響應直至穩(wěn)定。也就是利用動力學方法，解決懸掛管柱靜力屈曲問題[14]。

采用慢動力法，懸掛管柱屈曲的動力學基本運動方程為

Mu″+Cu′+Ku=F(t)

(4)

式中：M為懸掛管柱的質量矩陣；C為懸掛管柱的阻尼矩陣；K為懸掛管柱的剛度矩陣，K=K0+Kσ(u)+Kn(u)；u′和u″分別為懸掛管柱的節(jié)點速度和加速度向量；t為計算時間；F(t)為懸掛管柱的節(jié)點載荷向量，與方程(1)中的載荷相同，只是考慮了載荷施加的時間歷程，F(xiàn)(t)=Fq(t)+Fn(u)(t)。

采用Newmark直接積分法，對式(4)進行隱式有限元求解。若t時刻懸掛管柱的位移、速度和加速度已知，則可計算出t+Δt時刻懸掛管柱的動力響應。隨著計算時間的延長，懸掛管柱屈曲構型趨于穩(wěn)定，懸掛管柱各節(jié)點無速度和加速度，方程(4)可退化成方程(3)，實現(xiàn)用慢動力法對懸掛管柱靜力屈曲的求解。

2.2 外載荷

為了減小管柱上端懸掛拉力FH和自重均布載荷q對管柱的動力沖擊，采用正弦波方式，施加隨時間變化的外載荷，如圖2所示。上端懸掛拉力FH施加的方式為：①在管柱1/4倍固有周期之前，施加正弦波的懸掛拉力；②在管柱1/4倍固有周期之后，施加恒定的懸掛拉力。

(a) 三維視圖(b) 俯視圖(c) 載荷-時間歷程

圖2 外載荷和擾動力的施加方式

Fig.2 The mode of external load and disturbance load

同理，為了施加管柱自重載荷，施加的重力mg也采用同樣的正弦波方式。

2.3 擾動力

由于初始幾何缺陷會影響管柱屈曲的臨界載荷，而微小擾動力是在初始施加，在后續(xù)分析中撤銷，不會影響管柱屈曲的臨界載荷。本文采用初始擾動力方式，既保持了管柱結構的完整性，又實現(xiàn)了初始缺陷的施加。

為了易形成空間螺旋屈曲構型，在初始擾動力的施加過程中，使管柱變形成空間曲線，將懸掛管柱受壓段長度均分4段，每1/4長度處(p1、p2和p3位置)施加互成90°的擾動力Fp，管柱下端施加1個擾動扭矩Mp。施加擾動力的大小均為單位載荷，擾動力施加的方式如圖2所示。在管柱1/2倍固有周期之前，施加正弦波的擾動力；在管柱1/2倍固有周期之后，施加的擾動力為0。

2.4 阻尼比

方程(2)中的阻尼矩陣C一般為常用的Rayleigh阻尼，也就是比例阻尼，它包括α質量阻尼和β剛度阻尼。本文選取α質量阻尼，阻尼矩陣C=αM。以管柱第1階的固有頻率ω作為模態(tài)阻尼比，根據(jù)正交性原理，得到α質量阻尼與阻尼比ζ的關系為

α=2ωζ

(5)

采用慢動力分析法計算時，選取合適的阻尼比，得到質量阻尼，輸入方程(4)中，進行慢動力計算。

2.5 計算流程

圖3 計算流程框圖Fig.3 Flowchart of calculation flow

對懸掛管柱屈曲的動力學基本運動方程(2)進行有限元計算，計算總時間取n個周期，周期取第1階自振周期。計算結束后，判斷管柱各節(jié)點的撓度u是否隨時間變化。若是，表明管柱各節(jié)點還有速度u′和加速度u″，存在附加的慣性力Mu″和阻尼力Cu′，應減小阻尼比ζ，重新計算，直到管柱各節(jié)點的撓度u為常值，實現(xiàn)屈曲構型穩(wěn)定。然后判斷撓曲線是否是空間螺旋屈曲構型，若不是，增加受壓段無量綱長度ηC重新計算，直到得到穩(wěn)定的空間螺旋屈曲構型。

3 計算結果與分析

3.1 計算參數(shù)

3.2 無接觸正弦屈曲臨界載荷分析

將空間懸掛管柱簡化成平面問題，除了上下兩端鉸支約束外，還約束管柱所有節(jié)點y方向的平動自由度和繞x方向的轉動自由度。不考慮與井筒的接觸，研究懸掛管柱無接觸的正弦屈曲問題。將圖3計算流程中的螺旋屈曲構型判斷改成正弦屈曲構型，進行有限元計算。將施加擾動力在1/4固有周期時刻的最大撓度進行無量綱化，即u/u|t=0.25T=1，圖4給出了平面模型懸掛管柱最大撓度位置處的撓度曲線。

圖4 懸掛管柱平面模型的撓度曲線Fig.4 Deflection curve of suspended tubular string in plane model

由圖4可見，當井底壓力ξC<1.94時，懸掛管柱撓度隨橫向擾動力的撤銷而逐漸減小，懸掛管柱處于穩(wěn)定平衡狀態(tài)；當井底壓力ξC=1.94時，懸掛管柱撓度隨橫向擾動力的撤銷保持不變，懸掛管柱處于隨遇平衡狀態(tài)；當井底壓力ξC>1.94時，懸掛管柱撓度隨橫向擾動力的撤銷而逐漸增大，懸掛管柱處于不穩(wěn)定狀態(tài)。由此可見，本文計算得到的懸掛管柱平面正弦屈曲的臨界壓力ξC=1.94。文獻[1]的結果ξC=1.94，文獻[7]的結果ξC=1.89，本文結果與文獻[1]吻合，比文獻[7]大了2.65%，驗證本文計算方法的正確性。

3.3 有接觸的正弦屈曲過程分析

按照3.2節(jié)平面模型，同時考慮懸掛管柱與井筒的接觸，研究其正弦屈曲過程。將懸掛管柱在井筒內x方向的變形無量綱化，與井筒接觸時的x正方向橫向變形量記為1，對應的x負方向記為-1，圖5給出了不同井底壓力下的撓曲線。

圖5 平面模型中一次彎曲和二次彎曲的撓曲線Fig.5 First order and second order buckling in plane model

由圖5可見，由于井底壓力ξC=1.94時懸掛管柱處于隨遇平衡狀態(tài)，考慮井筒約束后，將與井筒接觸，接觸點位于ξL1=1.80處，撓曲線為一次彎曲狀態(tài)。當井底壓力ξC=3.76時，懸掛管柱下段在x方向為正值，懸掛管柱上段在x方向剛出現(xiàn)極小的負值，接觸點位于ξL1=1.41處，撓曲線為二次彎曲狀態(tài)。井底壓力ξC=3.76為二次彎曲臨界載荷，與文獻[1]的結果ξC=3.75基本吻合。

隨著井底壓力的增加，懸掛管柱下段接觸點逐漸下移，上段x負方向的撓度逐漸增大。當井底壓力ξC=4.31時，懸掛管柱上段剛好與井壁接觸，上段接觸點位于ξL12=4.03處，下段接觸點位于ξL1=1.10處，撓曲線為完全二次彎曲狀態(tài)。

3.4 螺旋屈曲分析

(a) p1位置

(b) p2位置

(c) p3位置圖6 井底壓力ξC=4.11時阻尼比對撓度的影響Fig.6 Influence of damping ratio on deflection when ηC=4.11

(a) ξC=4.10

(b) ξC=4.11

(c) ξC=4.12圖7 阻尼比對屈曲構型的影響Fig.7 Influence of damping ratio on buckling configuration

由圖6可見，①阻尼比取1.50≤ζ≤20.0，在6個周期內，管柱在p2和p3位置處的撓度在3個周期后均處于穩(wěn)定狀態(tài)，但在p1位置的撓度在6個周期時間內一直在持續(xù)增加。由此可見，在大阻尼比條件下，雖然橫向撓度曲線不隨時間劇烈震蕩，計算容易收斂，但使動力響應處于穩(wěn)定、達到常值狀態(tài)的時間會延長，也影響計算收斂效率。②當阻尼比ζ=1.00時，管柱在p1和p2位置的撓度在6個周期內剛好穩(wěn)定，在p3位置的撓度在1個周期內處于穩(wěn)定；當阻尼比ζ=0.75時，管柱在p1和p2位置的撓度在2個多周期內處于穩(wěn)定，在p3位置的撓度在1個周期內處于穩(wěn)定；當阻尼比ζ=0.50時，管柱在p1、p2和p3位置的撓度在1個周期內均處于穩(wěn)定。由此可見，阻尼比取0.50≤ζ≤1.00，在6個周期內，管柱的撓度均能處于穩(wěn)定、達到常值狀態(tài)。③當阻尼比ζ≤0.25時，管柱在這3處位置處的撓度出現(xiàn)持續(xù)震蕩。阻尼比越小，振動越劇烈。在p1位置，管柱與井筒出現(xiàn)時而接觸碰撞，時而脫離井筒的現(xiàn)象。由此表明，取較小的阻尼比時，撓度難以處于穩(wěn)定、達到常值狀態(tài)，計算難收斂。

通過圖6中阻尼比對撓度的影響分析可見，在給定的井底壓力和計算時間范圍條件下，為了獲得穩(wěn)定的屈曲構型，應該先取較大的阻尼比，然后逐漸減小阻尼比，直到懸掛管柱的撓度趨于穩(wěn)定、達到常值狀態(tài)為止。但是阻尼比過小，易引起管柱劇烈振動，在給定的計算時間范圍內，計算難收斂。

從圖7可見，①井底壓力取ξC=4.10，阻尼比較大時(8.00≤ζ≤20.0)，xy平面投影為不同形狀的曲線；當阻尼比ζ為4.00和6.00時，xy平面投影為直線且重合，為穩(wěn)定的平面正弦屈曲狀態(tài)，表明井底壓力ξC=4.10時屈曲變形為正弦構型。②井底壓力取ξC=4.11，阻尼比較大時(1.50≤ζ≤20.0)，xy平面投影也為不同形狀的曲線，但是屈曲變形并沒有處于穩(wěn)定、達到常值的狀態(tài)(見對圖6的分析)；當阻尼比ζ為0.75和1.00時，xy平面投影為曲線且重合，管柱上存在2點與井筒接觸，為穩(wěn)定的螺旋屈曲狀態(tài)，表明井底壓力ξC=4.11時屈曲變形剛好從正弦突變成螺旋構型。③井底壓力取ξC=4.12，阻尼比較大(10.0≤ζ≤20.0)時，xy平面投影也為不同形狀的曲線；當阻尼比較小時(2.00≤ζ≤8.00)時，也得到了穩(wěn)定的螺旋屈曲狀態(tài)，表明井底壓力ξC=4.12時屈曲變形為螺旋構型。

在6個周期的計算時間范圍內，圖8給出了阻尼比ζ對計算效率的影響，整體來看，計算耗費的CPU時間隨著阻尼比的增大而逐漸縮短。但從對圖6的分析可見，阻尼比過小或過大均影響計算的穩(wěn)定性。其中臨界載荷ξC=4.10時，由于處于正弦屈曲構型向螺旋屈曲轉變的前一個臨界載荷，屈曲構型極不穩(wěn)定，計算耗費的CPU時間較長。

圖8 阻尼比ζ對計算效率的影響Fig.8 Influence of damping ratio on computational efficiency

綜上可見，本文計算得到的懸掛管柱平面正弦屈曲向空間螺旋屈曲轉變時的臨界載荷為ξC=4.11。文獻[7]的結果為ξC=4.00，本文的結果比文獻[7]大了2.75%。

根據(jù)圖7的計算結果，圖9給出了懸掛管柱從正弦屈曲到螺旋屈曲的三維空間屈曲構型。由圖9可見，①當井底壓力ξC≤4.10時，懸掛管柱的屈曲構型為一平面，為平面正弦屈曲狀態(tài)。②當井底壓力ξC=4.11時，懸掛管柱的屈曲構型突變成空間螺旋曲線，為空間螺旋屈曲構型。2個接觸點的軸向位置分別為ξL1=1.30和ξL2=3.13，螺旋段長度為ξL2-ξL1=1.83，螺旋角度為θL2-θL1=77.9°。③當井底壓力ξC>4.11時，螺旋屈曲構型進一步擴展，當井底壓力ξC=4.31時，懸掛管柱在平面模型中為完全二次彎曲狀態(tài)，但在三維模型中早已成為空間螺旋屈曲狀態(tài)。

本文無量綱長度取ξL=8，計算得到了懸掛管柱正弦屈曲向螺旋屈曲轉變時的臨界載荷。根據(jù)文獻[7]，由于管柱越長管柱柔度越大，可以預見對應的臨界載荷會隨著管柱長度的增加而逐漸減小，最后趨于常值。

圖9 三維空間的屈曲構型Fig.9 Buckling configurations in three dimensional space

4 結 論

本文基于慢動力有限元法，對懸掛管柱進行了正弦屈曲向螺旋屈曲轉變時的臨界載荷研究，得出如下結論：

(1) 正弦屈曲向螺旋屈曲轉變是一個屈曲構型突變的過程，采用慢動力有限元法，將靜力屈曲問題轉換成動力學問題。按照正弦波的方式施加外載荷和擾動力，通過虛擬較大的阻尼比，對井筒內懸掛管柱后屈曲分析的收斂和算法的穩(wěn)定性具有重要的作用。阻尼力能夠有效地抑制管柱振動，在6個固有周期的計算時間范圍內，管柱的撓度均能處于穩(wěn)定、達到常值狀態(tài)，可以得到懸掛管柱穩(wěn)定的屈曲構型和臨界載荷。

(2) 無量綱長度取8，懸掛管柱正弦屈曲向螺旋屈曲轉變時的臨界載荷為4.11，螺旋屈曲的構型中存在兩個接觸點，螺旋段長度為1.83，螺旋角度為77.9°。懸掛管柱正弦屈曲和螺旋屈曲的臨界載荷均比文獻[7]大了近3%，但是本文一次彎曲和二次彎曲的臨界載荷與Lubinski的研究結果完全吻合，表明本文的結果具有更高的計算精度。

(3) 采用慢動力法，也可以用于諸如定向井和水平井等井型的屈曲分析，計算出穩(wěn)定的正弦和螺旋屈曲構型以及對應的臨界載荷。