王榮鵬, 宋桂秋, 周世華

(東北大學 機械工程與自動化學院, 遼寧 沈陽 110819)

深孔鉆柱系統(tǒng)主要用于煤炭、石油和天然氣生產(chǎn)和勘探過程,系統(tǒng)內(nèi)部諸多因素之間的相互作用較為復雜[1].在鉆柱鉆進過程中會發(fā)生各種振動,尤其是鉆柱的橫向振動能夠引起鉆柱部件失效、磨損、鉆頭損壞、鉆速降低,導致鉆井效率低下并產(chǎn)生高昂的成本,從而對鉆柱系統(tǒng)產(chǎn)生重大的不利影響.因此研究鉆柱系統(tǒng)的橫向振動特性,對于降低成本、提高效率具有重要的現(xiàn)實意義[2].目前國內(nèi)外學者對鉆柱系統(tǒng)研究作出了許多貢獻:劉永升等[3]建立了斜直井中鉆柱4自由度非線性動態(tài)模型,在不同井斜角條件下,對鉆柱運動進行了數(shù)值模擬研究.王長進等[4]通過實驗研究了轉(zhuǎn)速、載荷、鉆井液等參數(shù)對受壓屈曲鉆柱旋轉(zhuǎn)運動的影響,分析了鉆柱在鉆井液中由正向渦動轉(zhuǎn)為自轉(zhuǎn)和反向渦動的耗散能量.王文龍等[5]建立了鉆柱縱向振動數(shù)學模型,分析了鉆柱軸向應力振幅分布曲線和最大軸向應力幅頻曲線的特征,及鉆井液黏度、激勵位移和減振器位置等參數(shù)對鉆柱軸向應力振幅的影響.王尊策等[6]建立了變截面鉆柱的摩擦接觸有限元模型,對鉆柱的屈曲進行了分析研究.李子豐[7]提出鉆柱渦動理論研究需要將鉆井液動力潤滑學與鉆柱動力學相結(jié)合.Liu等[8]研究了鉆柱的黏滑振動,基于路徑跟蹤方法對非光滑動力系統(tǒng)進行了數(shù)值模擬,重點研究了鉆柱的穩(wěn)態(tài)問題.Real等[9]提出了一種用于鉆柱扭轉(zhuǎn)動力學的滯后鉆頭巖石相互作用模型,分析了鉆柱系統(tǒng)的非線性扭轉(zhuǎn)振動和穩(wěn)定性.Kamel等[10]提出了一個軸向和扭轉(zhuǎn)運動耦合的鉆柱非線性模型,考慮鉆頭阻力,采用接觸數(shù)學模型描述鉆頭和巖層的相互作用.Liang等[11-13]提出了一種旋轉(zhuǎn)多跨流體輸送管道的動力學模型,并對該系統(tǒng)的橫向固有頻率、共振頻率和模態(tài)特性進行了研究.Zhang等[14]研究了在超臨界流速作用下,黏彈性輸流管道在彎曲平衡附近的非線性強迫振動,重點研究了外共振和內(nèi)共振動力學.Guzek等[15]研究了鉆井泥漿流變性對垂直鉆頭振動的影響,確定了減振的最佳泥漿參數(shù)范圍.Nandakumar等[16]建立了一個耦合的2自由度模型,研究了鉆柱的黏滑振動、旋轉(zhuǎn)、鉆頭彈跳、螺旋屈曲等現(xiàn)象.Leonov等[17]研究了常微分方程描述的鉆井系統(tǒng)的數(shù)學模型,在這類系統(tǒng)中可能會出現(xiàn)隱藏振蕩等復雜效應,這些影響可能導致鉆柱失效和故障.Ritto等[18]提出了鉆柱扭轉(zhuǎn)振動的集中參數(shù)模型,采用非線性摩擦力矩模擬鉆頭的相互作用,分析了鉆柱扭轉(zhuǎn)黏滑振動不穩(wěn)定性.Hovda[19]建立了油井鉆柱軸向運動的集中參數(shù)模型,研究了井下壓力和運動的各種重要因素,特別是附加質(zhì)量、穩(wěn)態(tài)黏滯力的影響.Zhao等[20]采用摩擦模型對摩擦自激振動進行了分析,研究了鉆柱在四種接觸狀態(tài)下(黏著、邊界潤滑、部分流體潤滑和全潤滑)的軸向黏滑運動.Ghasemloonia等[21]建立了整個鉆柱軸向與橫向耦合非線性動力學模型,并對其橫向不穩(wěn)定性進行了定性研究.

由文獻分析可知,雖然近年來國內(nèi)外學者對鉆柱系統(tǒng)動力學特性進行了一些研究,然而考慮鉆柱中流體影響和支撐剛度影響的非線性動力學行為的研究還相對較少.本文建立了鉆柱系統(tǒng)流固耦合動力學模型,并用Galerkin方法進行了離散,采用Runge-Kutta積分法求解動力學方程.在其他參數(shù)固定的情況下,分別給出了以脈動頻率、脈動幅值和質(zhì)量比為控制參數(shù)的耦合系統(tǒng)振動響應的分岔圖.通過分岔圖、時域圖、相圖和Poincaré截面圖,分析了鉆柱系統(tǒng)的非線性動力學行為,為了解深孔鉆柱系統(tǒng)復雜的非線性動力學特性提供了依據(jù).

1 系統(tǒng)非線性動力學模型

在鉆柱系統(tǒng)中,當鉆柱只作軸向進給,可將鉆柱簡化為輸送脈動流體的簡支梁,如圖1所示.無約束的輸送脈動流體的簡支梁運動方程由文獻[22-24]提出,本文在此運動方程的基礎上,考慮了鉆柱扶正器對鉆柱振動的影響,將扶正器簡化為彈簧支撐約束[25],用三次非線性彈簧力Fa(y)和Fb(y)模擬反作用力[26],整理后得到深孔鉆柱系統(tǒng)流固耦合運動微分方程:

(1)

式中:mf為單位長度流體質(zhì)量;mp為單位長度鉆柱質(zhì)量;U為流體流速;E為楊氏模量;Ap為鉆柱的橫截面面積;l為鉆柱長度;Ip為鉆柱橫截面慣量;a為黏彈性阻尼系數(shù);xa和xb分別為沿鉆柱軸線距原點的距離;Fa(y)=Kay3,Fb(y)=Kby3，Ka和Kb分別為xa和xb位置處的支撐剛度;x為橫坐標;t為時間;y(x,t)為位置x與時間t的橫向振動位移函數(shù);δ(·)為Dirac delta函數(shù).

對式(1)引入無量綱變量:

(2)

方程(1)可寫為

kbη3δ(ξ-ξb)=0 .

(3)

假設鉆柱中流動切削液為周期性脈動流,即無量綱流速為

u=u0(1+μsinωτ) .

(4)

式中:u0為平均流速;μ為脈動幅值;ω為脈動頻率.

圖1 具有運動約束兩端鉸支鉆柱模型Fig.1 A drilling string model with motion constraints hinged at both ends

2 運動方程的Galerkin離散化

利用Galerkin方法對無量綱化的高階非線性偏微分方程(1)離散化并降階化為低次的常微分方程組.假設任意點無量綱位移作如下形式的Galerkin展開式:

(5)

式中:ξ為廣義坐標;τ為時間;η(ξ,τ)為任意點ξ無量綱位移;φr(ξ)為正交模態(tài)函數(shù);qr(τ)為關(guān)于時間的未知廣義坐標;N為截斷項數(shù).將式(5)代入式(3)中,在[0,1]區(qū)間關(guān)于ξ積分,并利用模態(tài)函數(shù)正交性,經(jīng)過整理后得到降階后的離散化方程為

(6)

式中:C為阻尼矩陣;K為剛度矩陣.由文獻[27]可知,對式(5)取前兩階展式,即N=2可得到較高的計算精度.其中，

(7)

(8)

式中:

G(z)=[0 0 -f1-f2]T;

Q(z)=[0 0 -h1-h2]T.

3 仿真結(jié)果與數(shù)據(jù)分析

鉆柱系統(tǒng)的模型參數(shù)如表1所示.

表1 模型參數(shù)

3.1 脈動頻率ω的影響

脈動頻率對鉆柱系統(tǒng)的動力學行為有著重要的影響,圖2a給出了支撐剛度系數(shù)ka=kb=0時,其他系統(tǒng)參數(shù)取固定常數(shù),脈動頻率ω在[30,50]范圍內(nèi),系統(tǒng)響應隨脈動頻率ω變化的分岔圖.圖2b為支撐剛度系數(shù)ka=2.1×105,kb=1.3×103時,ω在[30,50]范圍內(nèi),系統(tǒng)響應隨脈動頻率ω變化的分岔圖.圖2中A,B,C分別表示周期運動、擬周期運動和混沌運動,縱坐標為鉆柱中心位置振動的位移幅值.由圖2a可知,ω在[30,33.65]范圍內(nèi)時,系統(tǒng)出現(xiàn)周期運動、擬周期運動和小窗口的混沌運動形式,當ω=32.1時系統(tǒng)出現(xiàn)跳躍間斷的非線性現(xiàn)象.圖3a,3b,3c為ω=32.71時的時域圖、相圖和Poincaré截面圖,Poincaré截面圖為規(guī)則幾何形狀的橢圓,系統(tǒng)為擬周期運動.隨著ω逐漸增大,系統(tǒng)演變?yōu)槎虝旱闹芷?運動,然后當ω在(34.25,42.25)范圍內(nèi)時,系統(tǒng)演變?yōu)榛煦邕\動,當ω=35.28時,如圖3d,3e,3f所示,Poincaré截面圖出現(xiàn)幾何分形結(jié)構(gòu),系統(tǒng)為混沌運動.隨著ω的繼續(xù)增大,即ω在[42.25,50]范圍內(nèi),系統(tǒng)演變?yōu)橹芷谶\動狀態(tài),出現(xiàn)倍周期倒分岔現(xiàn)象.當ω=42.75時,系統(tǒng)路徑倒分岔演變?yōu)橹芷?運動.當ω=43.57時,如圖3g,3h,3i所示,Poincaré截面圖吸引子為8個點,即周期8運動.在ω=47.85時系統(tǒng)路徑出現(xiàn)倒分岔由周期8運動演變?yōu)橹芷?運動,振動形態(tài)的周期數(shù)成倍的變化,系統(tǒng)路徑以倍周期倒分岔的形式由混沌運動過渡到周期運動.

圖2 系統(tǒng)響應隨脈動頻率ω變化的分岔圖Fig.2 Bifurcation diagram of system response with pulsating frequencies ω (a)—ka=kb=0; (b)—ka=2.1×105, kb=1.3×103.

由圖2b可知,ω在[30,30.96]范圍內(nèi)時,系統(tǒng)為小窗口的周期2運動.隨著ω逐漸增大,系統(tǒng)較圖2a提前演變?yōu)榛煦邕\動.當脈動頻率ω在[35.52,37.68]范圍內(nèi),系統(tǒng)演變?yōu)橹芷谶\動,周期運動區(qū)域較小.當ω繼續(xù)增大至(37.68,43.52)區(qū)間時,系統(tǒng)演變?yōu)閿M周期運動.隨著ω繼續(xù)增大,當ω在[43.52,50]范圍內(nèi),系統(tǒng)演變?yōu)榛煦邕\動.系統(tǒng)動力學響應豐富,周期運動、擬周期運動、混沌運動交替出現(xiàn),系統(tǒng)表現(xiàn)出更復雜的動力學行為.支撐剛度系數(shù)對系統(tǒng)響應影響顯著.

3.2 脈動幅值μ的影響

為了研究脈動幅值μ對鉆柱系統(tǒng)振動特性的影響,圖4a給出了支撐剛度系數(shù)ka=2.1×106,kb=4.5×106,以脈動幅值μ為控制參數(shù)在[0.2,0.6]范圍內(nèi)的分岔圖.圖4b為支撐剛度系數(shù)ka=2.1×105,kb=0時,以脈動幅值μ為控制參數(shù)在[0.2,0.6]范圍內(nèi)的分岔圖.圖4中A,B,C分別表示周期運動、擬周期運動和混沌運動.如圖4a所示,當μ較小時,即μ在[0.2,0.259]范圍內(nèi)時,系統(tǒng)響應為混沌運動狀態(tài),隨著μ的逐漸增大,即μ在(0.259,0.275)范圍內(nèi)時,系統(tǒng)進入小窗口的周期運動，隨著μ的繼續(xù)增大，當μ在[0.275，0.389]范圍內(nèi)時，系統(tǒng)響應主要為混沌運動，響應復雜.當μ在[0.304，0.345]范圍內(nèi)時，觀察到小窗口的周期運動.隨著μ的繼續(xù)增大，即μ在(0.389，0.581)范圍內(nèi)時，系統(tǒng)出現(xiàn)明顯的周期運動和擬周期運動形式.圖5a，5b，5c為μ=0.423時的時域圖、相圖和Poincaré截面圖，Poincaré截面圖上的點構(gòu)成一封閉曲線，系統(tǒng)為擬周期運動狀態(tài).當μ=0.515時，如圖5d，5e，5f所示，時域圖呈現(xiàn)周期變化，此時的相圖為一封閉曲線，Poincaré截面圖上吸引子為孤立的6個點，系統(tǒng)表現(xiàn)為周期6運動.當μ=0.549時系統(tǒng)出現(xiàn)倍周期分岔現(xiàn)象，系統(tǒng)由周期6運動演變?yōu)橹芷?2運動.當μ=0.573時系統(tǒng)由周期12運動通過倍周期分岔演變?yōu)橹芷?4運動，系統(tǒng)路徑以倍周期分岔的形式由周期運動過渡到混沌運動.在脈動幅值μ=0.557，μ=0.568時系統(tǒng)出現(xiàn)跳躍間斷的非線性現(xiàn)象.隨著μ的繼續(xù)增大，即μ在[0.581，0.6]范圍內(nèi)時，系統(tǒng)演變?yōu)榛煦邕\動，系統(tǒng)振動幅值變大.當μ=0.592時，如圖5g，5h，5i所示，系統(tǒng)呈現(xiàn)非周期變化，相對應的Poincaré截面圖為無規(guī)則散點，系統(tǒng)運動處在混沌狀態(tài).系統(tǒng)動力學響應豐富，周期運動、擬周期運動、混沌運動交替出現(xiàn)，系統(tǒng)動態(tài)特性較復雜.

圖3 不同脈動頻率ω下系統(tǒng)動態(tài)響應Fig.3 Dynamic response of the system with different pulsating frequencies ω (a)—ω=32.71,時域圖; (b)—ω=32.71,相圖; (c)—ω=32.71,Poincaré截面圖; (d)—ω=35.28,時域圖; (e)—ω=35.28,相圖; (f)—ω=35.28,Poincaré截面圖; (g)—ω=43.57,時域圖; (h)—ω=43.57,相圖; (i)—ω=43.57,Poincaré截面圖.

圖4 系統(tǒng)響應隨脈動幅值μ變化的分岔圖Fig.4 Bifurcation diagram of system response with pulsation amplitude μ (a)—ka=2.1×106, kb=4.5×106; (b)—ka=2.1×105, kb=0.

圖5 不同脈動幅值μ下系統(tǒng)動態(tài)響應Fig.5 Dynamic response of the system with different pulsation amplitude μ (a)— μ=0.423,時域圖; (b)— μ=0.423,相圖; (c)— μ=0.423,Poincaré截面圖; (d)— μ=0.515,時域圖; (e)— μ=0.515,相圖; (f)— μ=0.515,Poincaré截面圖; (g)— μ=0.592,時域圖; (h)— μ=0.592,相圖; (i)— μ=0.592,Poincaré截面圖.

由圖4b可知,當脈動幅值μ較小時,即μ在[0.2,0.262]范圍內(nèi)時,系統(tǒng)出現(xiàn)周期運動、擬周期運動和混沌運動交替形式.在μ=0.218時系統(tǒng)出現(xiàn)跳躍間斷的非線性現(xiàn)象.當μ在[0.251,0.262]范圍內(nèi)時,系統(tǒng)出現(xiàn)小窗口的周期運動.隨著μ的逐漸增大,即μ在(0.262,0.394)范圍內(nèi)時,系統(tǒng)演變?yōu)榛煦邕\動,混沌運動強度也隨之增大.隨著μ的繼續(xù)增大,即μ在[0.394,0.6]范圍內(nèi)時,系統(tǒng)主要以周期2運動為主.通過對比圖4a和圖4b可以發(fā)現(xiàn),選取不同支撐剛度系數(shù),系統(tǒng)響應發(fā)生明顯變化,系統(tǒng)的跳躍間斷現(xiàn)象提前,這主要是由于支撐剛度在一定程度上引起系統(tǒng)固有特性的改變,可見支撐剛度對系統(tǒng)的非線性動力學行為有復雜的影響.

3.3 質(zhì)量比β的影響

鉆柱系統(tǒng)質(zhì)量比β是評價系統(tǒng)性能的重要參數(shù).在本節(jié)中圖6a為支撐剛度系數(shù)ka=2.6×106,kb=2.1×106時,以質(zhì)量比β為控制參數(shù)在[0.2,0.8]范圍內(nèi)的分岔圖.圖6b為支撐剛度系數(shù)ka=2.6×106,kb=5.6×106時，以質(zhì)量比β為控制參數(shù)在[0.2,0.8]范圍內(nèi)的分岔圖.圖6中A,B,C分別表示周期運動、擬周期運動和混沌運動.如圖6a所示,當β較小時,即β在[0.2,0.261]范圍內(nèi)時,系統(tǒng)響應為周期運動和擬周期運動.圖7a,7b,7c為β=0.237時的時域圖、相圖及Poincaré截面圖,此時的相圖為一封閉曲線,Poincaré截面圖為2個點,系統(tǒng)表現(xiàn)為周期2運動.當β=0.246時,如圖7d,7e,7f所示,系統(tǒng)運動為擬周期運動狀態(tài).隨著β的增大,即β在(0.261,0.354)范圍內(nèi)時,系統(tǒng)響應為混沌運動,混沌運動區(qū)域逐漸增大.隨著β的繼續(xù)增大,系統(tǒng)響應出現(xiàn)短暫的周期運動.當β在[0.387,0.456]范圍內(nèi)時,系統(tǒng)響應主要為混沌運動,當β=0.389時,如圖7g,7h,7i所示,系統(tǒng)為混沌運動狀態(tài).隨著β的繼續(xù)增大,系統(tǒng)響應出現(xiàn)小窗口的周期運動.當β在[0.469,0.727]范圍內(nèi)時,系統(tǒng)演變?yōu)榛煦邕\動,混沌區(qū)域明顯增大.隨著β的進一步增大,即β在(0.727,0.8)范圍內(nèi)時,系統(tǒng)演變?yōu)橹芷?運動.

圖6 系統(tǒng)響應隨質(zhì)量比β變化的分岔圖Fig.6 Bifurcation diagram of system response with mass ratio β (a)— ka=2.6×106, kb=2.1×106; (b)— ka=2.6×106, kb=5.6×106.

圖7 不同質(zhì)量比β下系統(tǒng)動態(tài)響應Fig.7 Dynamic response of the system with different mass ratio β (a)—β=0.237,時域圖; (b)—β=0.237,相圖; (c)—β=0.237,Poincaré截面圖; (d)—β=0.246,時域圖; (e)—β=0.246,相圖; (f)—β=0.246,Poincaré截面圖; (g)—β=0.389,時域圖; (h)—β=0.389,相圖; (i)—β=0.389,Poincaré截面圖.

由圖6b可知,當β較小時,即β在[0.2,0.326]范圍內(nèi)時,系統(tǒng)響應為周期運動和擬周期運動.隨著β的增大,即β在(0.326,0.422)范圍內(nèi)時,系統(tǒng)演變?yōu)榛煦邕\動,混沌運動區(qū)域逐漸增大.隨著β的繼續(xù)增大,即β在[0.422,0.451]范圍內(nèi)時,系統(tǒng)響應出現(xiàn)小窗口的周期運動.隨著β的繼續(xù)增大,系統(tǒng)返回混沌運動,混沌運動區(qū)域有所增大.隨著β的進一步增大,即β在[0.567,0.8]范圍內(nèi)時,系統(tǒng)響應主要為周期運動,出現(xiàn)周期6運動、周期2運動.與圖6a相比,混沌運動區(qū)域明顯減小,小窗口的周期運動減少,周期運動區(qū)域增大.不同的支撐剛度對系統(tǒng)動力學固有特性有明顯影響.

4 結(jié) 論

1) 系統(tǒng)在以脈動頻率ω為控制參數(shù)下,表現(xiàn)出周期運動、擬周期運動、混沌運動狀態(tài).觀察到系統(tǒng)由混沌運動通往周期運動的倒分岔路徑.選取一定支撐剛度系數(shù)情況下,混沌運動提前,擬周期運動區(qū)域增大,系統(tǒng)運動的交換頻率增強.支撐剛度對系統(tǒng)響應影響顯著.

2) 隨著脈動幅值μ的增大,系統(tǒng)表現(xiàn)出跳躍間斷現(xiàn)象、周期運動、擬周期運動和混沌運動頻繁交替,系統(tǒng)路徑以倍周期分岔的形式由周期運動過渡到混沌運動.此外,在選取一定支撐剛度系數(shù)情況下,系統(tǒng)的跳躍間斷現(xiàn)象提前,支撐剛度在一定程度上引起系統(tǒng)固有特性的改變.

3) 隨著質(zhì)量比β的增大,系統(tǒng)出現(xiàn)周期性運動、擬周期運動和混沌運動狀態(tài).通過改變支撐剛度系數(shù),系統(tǒng)混沌運動強度減弱,混沌運動區(qū)域明顯減小,不同的支撐剛度對系統(tǒng)動力學固有特性產(chǎn)生明顯影響.