考慮分布軸向力的細長桿橫向振動與失穩分析1)

孫秀榮 張立娟 趙潔妤 宋 楊

(河北環境工程學院，河北秦皇島066102)

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

作者近幾年致力于桿管柱的屈曲失穩和橫向振動領域的研究[14-16]，在前期經驗積累的基礎上，建立了本文分布軸向力下的桿柱失穩和橫向振動的力學、數學模型，給出近年來熱門的數值計算方法，得出更為切合實際的數值分析結果，同時為高校力學方面的教學研究和工程實際運用提供一些指導。

1 力學模型

不考慮細長桿柱的縱向振動與扭轉振動對橫向振動和屈曲失穩變形的影響，并做如下假設：細長桿柱為鉛垂桿；桿柱質地均勻；桿柱頂端為固定端約束，底端為滑動鉸支約束。力學模型如圖1 所示。

圖1 抽油桿柱橫向彎曲力學模型

2 分布軸向力下桿柱失穩模型及臨界載荷計算方法

根據圖1 力學模型，建立細長桿柱的數學模型為

方程(1) 進一步化為

其中

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

邊界條件為

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

當q =0 時，方程(1) 化為

p0為負值時，桿柱受壓，式(4)與材料力學中的壓桿穩定性方程一致。

2.1 有限差分法求解

將桿柱劃分為n 段，得到n+1 個節點。差分格式為

將式(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 伽遼金法

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

其中

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

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

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

an有非零解的條件為

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

3 軸向分布力下的桿柱橫向振動數學模型及固有頻率計算方法

根據圖1 力學模型，建立桿柱的橫向振動數學模型為

邊界條件為

3.1 積分法

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

為便于計算，式(19) 中參數進行變換，則

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

式中

設z(ξ) 為多個多項式之和，采用級數解來逼近其真實解，則

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

式中

若cn有非零解，則

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

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

3.2 伽遼金法

令

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

設滿足邊界條件的形函數為梁自由振動函數

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

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

其中

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

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

4 兩種模型數值分析

4.1 桿柱失穩臨界載荷

4.1.1 準確性測試

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

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

4.1.2 考慮分布軸向力失穩載荷

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

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

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

4.2 桿柱橫向振動固有頻率和振型

4.2.1 準確性測試

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

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

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

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

4.2.2 考慮分布軸向力和桿長變化的橫向振動

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

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

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

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

圖6 桿柱變形規律

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

5 工程實例求解

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

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

表3 螺桿泵基本參數表

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

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

5.2 有桿抽油系統桿柱失穩情況數值分析

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

圖7 桿柱各階振型函數曲線

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

表5 基本參數

表6 失穩臨界載荷

圖8 失穩載荷隨桿長變化曲線

圖9 桿柱不同桿長下的失穩變形曲線

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

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

由圖9 可知，不同井深下，桿柱的失穩變形曲線的規律是近似的，即靠近底端處產生最大失穩變形，這與不考慮分布軸向力時在桿柱的中間位置處產生最大失穩變形有顯著不同。因此，抽油桿柱的桿管偏磨嚴重區域會出現在較深且靠近底端部位，這與文獻[16-17] 結論是一致的。

6 結論

本文建立了考慮分布軸向力的細長桿失穩和橫向振動的力學、數學模型，通過數值分析結果對比，得出以下結論：

(1)考慮分布軸向力的細長桿柱失穩的必要條件是其橫向振動的固有頻率為0。該結論是對先前教材的知識點的延伸，同時也適用于不考慮分布軸向力的桿柱失穩情況。

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

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

(4) 本文桿柱與《材料力學》中歐拉桿柱失穩突變處的位置明顯不同，這是因為本文模型在歐拉桿的基礎上考慮了分布軸向力因素所致。本文對《材料力學》的桿柱靜力失穩和《機械系統動力學》的梁的橫向振動知識進一步拓展和延伸，并提供了工程應用實例。