孫秀榮 張立娟 趙潔妤 宋 楊

(河北環(huán)境工程學(xué)院，河北秦皇島066102)

歐拉桿柱失穩(wěn)臨界載荷的計(jì)算在高校教材《材料力學(xué)》中并未考慮桿柱的分布軸向力因素[1]，《機(jī)械系統(tǒng)動(dòng)力學(xué)》中梁的橫向振動(dòng)也并未考慮此因素[2]。目前高等教材中桿柱失穩(wěn)問題，主要集中在歐拉桿柱的臨界載荷的計(jì)算上[3-6]；涉及梁柱的橫向振動(dòng)，主要集中在歐拉伯努利梁的橫向振動(dòng)的仿真和計(jì)算上[7-13]。與此對(duì)應(yīng)，上述計(jì)算得到的結(jié)論也未涉及到梁柱的分布軸向力的影響。針對(duì)幾百甚至上千米的超長桿柱，分布軸向力對(duì)梁柱的失穩(wěn)和橫向振動(dòng)的影響問題，至今未看到有效的文獻(xiàn)直觀地體現(xiàn)在教材的拓展中。在近些年的教學(xué)過程中，不少教育工作者面臨此方面的困惑，卻由于參考文獻(xiàn)匱乏而得不到較好的解答。因此，目前對(duì)歐拉桿失穩(wěn)、歐拉伯努利梁問題進(jìn)行延拓是必要的，對(duì)工程領(lǐng)域也有較強(qiáng)的指導(dǎo)意義。

作者近幾年致力于桿管柱的屈曲失穩(wěn)和橫向振動(dòng)領(lǐng)域的研究[14-16]，在前期經(jīng)驗(yàn)積累的基礎(chǔ)上，建立了本文分布軸向力下的桿柱失穩(wěn)和橫向振動(dòng)的力學(xué)、數(shù)學(xué)模型，給出近年來熱門的數(shù)值計(jì)算方法，得出更為切合實(shí)際的數(shù)值分析結(jié)果，同時(shí)為高校力學(xué)方面的教學(xué)研究和工程實(shí)際運(yùn)用提供一些指導(dǎo)。

1 力學(xué)模型

不考慮細(xì)長桿柱的縱向振動(dòng)與扭轉(zhuǎn)振動(dòng)對(duì)橫向振動(dòng)和屈曲失穩(wěn)變形的影響，并做如下假設(shè)：細(xì)長桿柱為鉛垂桿；桿柱質(zhì)地均勻；桿柱頂端為固定端約束，底端為滑動(dòng)鉸支約束。力學(xué)模型如圖1 所示。

圖1 抽油桿柱橫向彎曲力學(xué)模型

2 分布軸向力下桿柱失穩(wěn)模型及臨界載荷計(jì)算方法

根據(jù)圖1 力學(xué)模型，建立細(xì)長桿柱的數(shù)學(xué)模型為

方程(1) 進(jìn)一步化為

其中

式中，L 為抽油桿柱全長，m；q 為抽油桿柱軸向分布力，N/m；E 為抽油桿柱材料彈性模量，GPa；ε 為抽油桿柱線性分布軸向力比例系數(shù)；I 為抽油桿柱的抗彎慣性矩，m4；p0為抽油桿柱端部受力(向下為正，受拉；向上為負(fù)，受壓)，N。

邊界條件為

式(3)即為含有分布軸向力的桿柱屈曲方程，難于求出精確的解析解。

當(dāng)q =0 時(shí)，方程(1) 化為

p0為負(fù)值時(shí)，桿柱受壓，式(4)與材料力學(xué)中的壓桿穩(wěn)定性方程一致。

2.1 有限差分法求解

將桿柱劃分為n 段，得到n+1 個(gè)節(jié)點(diǎn)。差分格式為

將式(5)～式(7) 代入方程(2) 離散化得

式中，ai= βh3/2-h2(1-X)β -4，bi= 2h2(1-X)β+6，ci=-βh3/2-h2(1-X)β-4。

由邊界條件可知

將式(8) 和式(9) 展開成矩陣形式得

即為

式(11)兩邊同乘A-1，可化簡為pY =A-1BY，求出A-1B 的最小特征值，即可求得p，則臨界載荷Fcr=-p0=-pEI/L2。

2.2 伽遼金法

設(shè)滿足邊界條件(3) 的形函數(shù)為

其中

式中，λn= βnL 由超越函數(shù)cos λnchλn= 1 確定，λn≈(n+1/4)π。

同理，得到迦遼金方程為整理方程(14)，得到關(guān)于一待定系數(shù)列陣[an] 的齊次線性方程組

其中，Lnn= (λn)4I1nn- pI2nn- βI3nn，Lmn=(λn)4I1mn-pI2mn-βI3mn+βI4mn。I1nn，I2nn，I3nn，I1mn，I2mn，I3mn，I4mn為積分常數(shù)。

an有非零解的條件為

可求得失穩(wěn)時(shí)的臨界載荷，將其代入式(16) 可得到振型函數(shù)。

3 軸向分布力下的桿柱橫向振動(dòng)數(shù)學(xué)模型及固有頻率計(jì)算方法

根據(jù)圖1 力學(xué)模型，建立桿柱的橫向振動(dòng)數(shù)學(xué)模型為

邊界條件為

3.1 積分法

先將y(x，t)分離變量，即令y(x，t)=Z(x)F(t)，則方程(17) 可化為

為便于計(jì)算，式(19) 中參數(shù)進(jìn)行變換，則

對(duì)式(19) 進(jìn)行四次積分，并整理得

式中

設(shè)z(ξ) 為多個(gè)多項(xiàng)式之和，采用級(jí)數(shù)解來逼近其真實(shí)解，則

將式(22) 代入式(20) 得

式中

若cn有非零解，則

由式(24) 求解出Ω，即可進(jìn)一步得到桿柱橫向振動(dòng)固有頻率ω。特殊地，當(dāng)q =0 時(shí)，p0=0，方程(17)化為

式(25) 與教材[2] 中梁的橫向自由振動(dòng)方程是相同的，此方程有解析解，亦可采用此方法求解。

3.2 伽遼金法

令

將式(27) 代入方程(17)，得振型方程為

設(shè)滿足邊界條件的形函數(shù)為梁自由振動(dòng)函數(shù)

φn(X) 同式(13)。將式(28) 和式(29) 代入式(27)，得到迦遼金方程為

整理方程(30)，得到關(guān)于一待定系數(shù)列陣[an] 的齊次線性方程組為

其中

I1nn，I2nn，I3nn，I1mn，I2mn，I3mn，I4mn為積分常數(shù)。an有非零解的條件是

可求得各階固有頻率，將其代入式(32) 可得到振型函數(shù)。

4 兩種模型數(shù)值分析

4.1 桿柱失穩(wěn)臨界載荷

4.1.1 準(zhǔn)確性測試

當(dāng)分布軸向力q = 0 時(shí)，數(shù)學(xué)方程(1) 簡化為《材料力學(xué)》教材[1]中的模型，而數(shù)學(xué)方程(17) 簡化為《機(jī)械系統(tǒng)動(dòng)力學(xué)》教材[2]中的模型。其他參數(shù)為L = 5 m，d = 20 mm，E = 209 GPa 時(shí)，材料力學(xué)歐拉桿臨界載荷公式為

采用式(33) 歐拉法、2.1 節(jié)差分法、2.2 節(jié)伽遼金法分別求解相同參數(shù)下桿柱的臨界載荷為1322.5 N，1325.7 N，1321.2 N，可知，差分法和伽遼金法的結(jié)果與公式(33) 求解結(jié)果誤差微小(誤差在0.25%以內(nèi))，說明差分法和伽遼金法求解歐拉桿臨界載荷精度較高。

4.1.2 考慮分布軸向力失穩(wěn)載荷

取基本參數(shù)：E =209 GPa，d=20 mm。差分法和伽遼金法得出不同桿長下的桿柱臨界載荷規(guī)律如圖2 和圖3 所示。由兩圖可知，在桿長為50 m 以內(nèi)時(shí)，差分法和伽遼金法得到的臨界載荷規(guī)律相同，即隨桿長的增加而減??；但桿長為50 m 以上時(shí)，伽遼金法所得桿柱臨界載荷會(huì)隨桿長的增加而增加，這與差分法所得規(guī)律截然相反。由此可知，兩種方法在計(jì)算考慮分布軸向力較長桿柱的臨界載荷時(shí)有一種方法得到的結(jié)果是不準(zhǔn)確且不適用的。而伽遼金方法實(shí)質(zhì)是一種近似計(jì)算方法，在桿長超越極限時(shí)結(jié)果會(huì)失真，因而不再適用。這也說明近似方法會(huì)有它本身的適用范圍，而伽遼金法不適用于考慮分布軸向力且桿長超越一定長度的桿柱的臨界載荷的計(jì)算。

圖2 差分法所得臨界載荷變化

圖3 伽遼金法所得臨界載荷變化

4.2 桿柱橫向振動(dòng)固有頻率和振型

4.2.1 準(zhǔn)確性測試

情況一：當(dāng)基本參數(shù)為L = 5 m，d = 20 mm，E = 209 GPa，q = 0 N/m，且底部載荷p0= 0 時(shí)，桿柱則做橫向自由振動(dòng)，其固有頻率可表示[2]為

采用式(34)、3.1 節(jié)積分法和3.2 節(jié)伽遼金法求解該相同參數(shù)下桿柱的固有頻率，其結(jié)果對(duì)比如表1 所示。由表1 可知，前六階固有頻率所得結(jié)果誤差很小(0.15%以內(nèi))。因此，不考慮端部載荷和分布軸向力時(shí)，積分法和伽遼金法求解桿柱橫向振動(dòng)固有頻率和振型精度都很高。

表1 未考慮分布軸向力計(jì)算固有頻率(單位：HZ)

情況二：在情況一參數(shù)基礎(chǔ)上，只變動(dòng)p0，令p0=2000 N，得到桿柱橫向振動(dòng)的固有頻率，如表2所示。由表2 可得，底端桿柱受拉的情況下，積分法和伽遼金法求解桿柱固有頻率時(shí)誤差甚小(0.02%以內(nèi))，說明兩種方法均適應(yīng)于桿柱受拉情況桿柱固有頻率計(jì)算。

4.2.2 考慮分布軸向力和桿長變化的橫向振動(dòng)

情況三：固有頻率隨軸向分布力和桿長的變化規(guī)律。假設(shè)基本參數(shù)為d=20 mm，E=209 GPa，p0=500 N，改變軸向分布力q和L，得到桿柱固有頻率隨分布軸向力q和L的關(guān)系曲線，如圖4 所示。由圖4 可知，固有頻率隨分布軸向力的增大而增大，而隨桿柱長度的增長而減小。q變化范圍在0～20 N/m 時(shí)，50 m 內(nèi)桿柱長度固有頻率最大誤差率為26.96%，100～2000 m 內(nèi)桿柱長度固有頻率最大誤差率為390.6%。由此說明，在桿柱長度較短時(shí)，分布軸向力對(duì)固有頻率的影響不大，桿柱較長時(shí)，分布軸向力對(duì)固有頻率影響較大。桿柱較長時(shí)，分布軸向力不宜忽略。

圖4 固有頻率隨桿長L 和分布軸向力q 的變化曲線

情況四：固有頻率隨底端載荷變化的規(guī)律。考慮到分布軸向力在桿柱較長時(shí)不宜忽略，取基本參數(shù)q=20 N/m，L=1000 m，d=25 mm，E=209 GPa，得到不同載荷下桿柱的固有頻率變化曲線，如圖5 所示(圖中p0＞0 為桿柱受拉，反之受壓)。由圖可知，隨著底部受拉載荷的減小，固有頻率逐漸下降，底端載荷為0 后繼續(xù)減小載荷(即反向加載，受壓載荷逐漸增大)，受壓載荷越過一定值后，固有頻率變?yōu)樘摂?shù)。說明在底部載荷逐漸變化的過程中，有一載荷對(duì)應(yīng)著固有頻率為0 的情況。令一階固有頻率為0，由2.1 節(jié)求得靜力學(xué)下臨界載荷即p0=-169.8 N，失穩(wěn)變形曲線如圖6 所示。由3.1 節(jié)求得動(dòng)力學(xué)下橫向振動(dòng)的變形曲線同樣如圖6 所示。由圖可知，兩變形幾乎完全重合，桿柱失穩(wěn)時(shí)最大變形均出現(xiàn)在桿柱底部附近。說明此時(shí)靜力學(xué)下考慮分布軸向力q=20 N/m 桿柱失穩(wěn)的變形即為動(dòng)力學(xué)下橫向振動(dòng)變形，而p0則為失穩(wěn)時(shí)的臨界載荷。由此可知，考慮分布軸向力的桿柱在非零變形時(shí)的固有頻率為0，即一階固有頻率為0 時(shí)，桿柱達(dá)到失穩(wěn)的必要條件。特殊的，當(dāng)分布軸向力為0 且桿柱屈曲時(shí)，其固有頻率對(duì)應(yīng)為0，該結(jié)論與文獻(xiàn)[8] 一致。

圖5 固有頻率隨底端載荷的變化曲線

圖6 桿柱變形規(guī)律

圖6同時(shí)給出了不考慮分布軸向力(即q= 0)時(shí)桿柱的靜力失穩(wěn)變形曲線，與q= 20 N/m 時(shí)失穩(wěn)變形曲線對(duì)比可知，兩者的最大變形位置(即突變處) 明顯不同，即桿柱最先達(dá)到失穩(wěn)的位置點(diǎn)不同。不考慮分布軸向力的桿柱最先在靠近中間位置失穩(wěn)，而考慮分布軸向力的桿柱最先在靠近端部位置處失穩(wěn)(距離底端約為8.5 m 處)，這也說明了分布軸向力對(duì)桿柱最先失穩(wěn)位置有顯著的影響，即超長桿柱計(jì)算失穩(wěn)問題須考慮分布軸向力。

5 工程實(shí)例求解

5.1 考慮分布軸向力的螺桿泵桿柱固有頻率和振型

數(shù)學(xué)模型方程(17) 和方程(18) 可應(yīng)用在求解螺桿泵桿柱的固有頻率上，其參數(shù)如表3 所示。計(jì)算的前五階固有頻率如表4，其對(duì)應(yīng)的振型函數(shù)如圖7所示。

表3 螺桿泵基本參數(shù)表

表4 考慮分布軸向力的固有頻率結(jié)果(單位：HZ)

由表4 可知，當(dāng)細(xì)長桿柱底端受拉時(shí)(p0＞0)，積分法和伽遼金方法得到的固有頻率誤差很小(0.5%以內(nèi))，可以忽略不計(jì)。圖7 給出了積分法和伽遼金方法求解前四階桿柱的振型函數(shù)曲線，經(jīng)對(duì)比可知，兩種方法用于求解桿柱固有頻率和振型函數(shù)，誤差很小，精度較高。

5.2 有桿抽油系統(tǒng)桿柱失穩(wěn)情況數(shù)值分析

有桿抽油系統(tǒng)抽油桿柱的基本參數(shù)如表5 所示。實(shí)際油井淺則幾百米深則幾千米，根據(jù)前面4.1 節(jié)的結(jié)果分析，伽遼金方法不適用于求解該長度桿柱的臨界載荷，只能采用有限差分法。其不同桿長下臨界載荷情況如表6 和圖8 所示，其對(duì)應(yīng)失穩(wěn)變形情況如圖9 所示。

圖7 桿柱各階振型函數(shù)曲線

圖7 桿柱各階振型函數(shù)曲線(續(xù))

表5 基本參數(shù)

表6 失穩(wěn)臨界載荷

圖8 失穩(wěn)載荷隨桿長變化曲線

圖9 桿柱不同桿長下的失穩(wěn)變形曲線

由表6 和圖8 可知，對(duì)于桿柱長度300 ～2000 m的桿柱，考慮分布軸向力和不考慮分布軸向力，其失穩(wěn)臨界載荷相差甚大，因此油田抽油桿柱求解偏磨臨界載荷時(shí)，分布軸向力不宜忽略不計(jì)。

由圖8 進(jìn)一步可知，分布軸向力q = 20 N/m，桿柱長度為200 ～2000 m 時(shí)，其臨界載荷隨桿長的增加而逐漸減小，直至趨于一穩(wěn)定值，變化率區(qū)間約為0%～13.6%。說明200 m 的桿柱代替2000 m 桿柱計(jì)算臨界載荷時(shí)，最大誤差為13.6%。由此給出工程建議，在精度要求不太高的情況下，考慮分布軸向力的有桿抽油系統(tǒng)，求解桿管偏磨臨界載荷時(shí)可以以200 m 左右的桿柱代替幾千米的情況進(jìn)行計(jì)算。

由圖9 可知，不同井深下，桿柱的失穩(wěn)變形曲線的規(guī)律是近似的，即靠近底端處產(chǎn)生最大失穩(wěn)變形，這與不考慮分布軸向力時(shí)在桿柱的中間位置處產(chǎn)生最大失穩(wěn)變形有顯著不同。因此，抽油桿柱的桿管偏磨嚴(yán)重區(qū)域會(huì)出現(xiàn)在較深且靠近底端部位，這與文獻(xiàn)[16-17] 結(jié)論是一致的。

6 結(jié)論

本文建立了考慮分布軸向力的細(xì)長桿失穩(wěn)和橫向振動(dòng)的力學(xué)、數(shù)學(xué)模型，通過數(shù)值分析結(jié)果對(duì)比，得出以下結(jié)論：

(1)考慮分布軸向力的細(xì)長桿柱失穩(wěn)的必要條件是其橫向振動(dòng)的固有頻率為0。該結(jié)論是對(duì)先前教材的知識(shí)點(diǎn)的延伸，同時(shí)也適用于不考慮分布軸向力的桿柱失穩(wěn)情況。

(2) 計(jì)算考慮分布軸向力的桿柱受壓失穩(wěn)載荷時(shí)，伽遼金近似計(jì)算法在求解較長桿柱時(shí)超出允許精度范圍導(dǎo)致計(jì)算結(jié)果失真；差分法則適用于求解較長桿柱的失穩(wěn)載荷，且精度較高。對(duì)于端部受拉的桿柱，數(shù)值法和伽遼金法計(jì)算桿柱固有頻率和振型函數(shù)的方法都是適用的。

(3)對(duì)于抽油桿柱等超長桿，計(jì)算失穩(wěn)臨界載荷時(shí)考慮分布軸向力是必要的，反之會(huì)產(chǎn)生較大的誤差。當(dāng)考慮分布軸向力的超長桿柱(如200 m 以上的抽油桿柱)，其失穩(wěn)臨界載荷隨桿長的變化逐漸緩慢減小。精度要求不高時(shí)，可不必取桿柱全長，取一部分桿柱計(jì)算臨界載荷即可。

(4) 本文桿柱與《材料力學(xué)》中歐拉桿柱失穩(wěn)突變處的位置明顯不同，這是因?yàn)楸疚哪Ｐ驮跉W拉桿的基礎(chǔ)上考慮了分布軸向力因素所致。本文對(duì)《材料力學(xué)》的桿柱靜力失穩(wěn)和《機(jī)械系統(tǒng)動(dòng)力學(xué)》的梁的橫向振動(dòng)知識(shí)進(jìn)一步拓展和延伸，并提供了工程應(yīng)用實(shí)例。