矩陣法求解壓桿穩(wěn)定的有限差分歐拉方程

胡章詠

(黃岡師范學院 機電工程學院，湖北 黃州 438000)

?

矩陣法求解壓桿穩(wěn)定的有限差分歐拉方程

胡章詠

(黃岡師范學院 機電工程學院，湖北 黃州 438000)

將壓桿穩(wěn)定的歐拉方程用差分的形式表達后，根據(jù)壓桿穩(wěn)定的實際要求，通過矩陣的方法求解該差分方程的系數(shù)矩陣，最后求出簡支梁等截面直桿的壓桿穩(wěn)定臨界應力。將求出的結(jié)果與解析解進行對比并進行誤差分析，找出產(chǎn)生誤差的原因和減小誤差的措施，并指出該方法在解決變截面桿的實際意義。

壓桿穩(wěn)定；有限差分；歐拉方程；矩陣

承壓桿件的穩(wěn)定性問題一直是材料力學當中的三大問題之一(材料力學研究桿件的三大問題為：桿件的強度、剛度和穩(wěn)定性)，相比較于桿件的強度和剛度問題，桿件的穩(wěn)定性在計算和模擬方面有著一定的難度，特別是對于變截面桿件的穩(wěn)定性問題，其臨界應力的計算目前用解析法還是比較困難。

目前，解決壓桿穩(wěn)定問題的主要方法有：針對直桿定截面的壓桿穩(wěn)定問題主要用歐拉方程來求解，也有用能量法和圖乘法來進行求解[1]；針對變截面桿，求解的方法有很多，主要有能量法、數(shù)值法等[2-3]。

1　問題的引入

如圖1所示，細長壓桿兩端為球鉸支座，軸線為直線，壓力P與軸線重合。當壓力達到臨界值時，壓桿將由直線平衡形態(tài)轉(zhuǎn)變?yōu)榍€平衡形態(tài)。選取坐標系如圖1所示，距離原點為x的任意截面的撓度為v，彎矩M的絕對值為Pv，對于微小的彎曲變形撓曲線的近似方程為：

圖1　簡支梁力學模型

(1)

引入記號：

(2)

于是式(1)可以寫成[1]：

(3)

式(2)的解析解法在參考文獻[4]中有詳細討論，本文不再贅述。

2　建立差分方程并確定求解方法

式(2)的差分離散化問題做如下討論：

首先將梁等分成等份(如圖2所示)，

圖2　簡支梁差分離散模型

v在x=xi處的一階和二階倒數(shù)可以表達為差分公式[5]：

(4)

(5)

將式(5)代入式(3)中,有差分方程：

(6)

不妨設(shè)將細桿分成n等分，x對應的坐標分別為[x0,x1,x2,…,xn-1,xn]，v對應xi的坐標分別為：[v0,v1,v2,…,vn-1,vn]，且令：

(7)

則式(6)可以分別離散成為下面的方程組：

(8)

又因為根據(jù)桿件的邊界條件：x0=0時，v0=0；xn=l時，vn=0；

所以式(8)變成了以下的一個n元一次其次方程組：

(9)

只有當[v0,v1,v2,…,vn-1,vn]有異于0的解時，才是失穩(wěn)的情況，因此就要求式(8)系數(shù)行列式等于0，即：

(10)

將式(10)解得的a的值反代入式(7)和式(2)中，求出P：

(11)

但由于式(10)是一個一元高次方程，由此可以得出a的解時多個數(shù)值，根據(jù)該問題的實際意義，即在多個與a值對應的P值中，使桿件保持微小彎曲的最小壓力，才是臨界應力Pcr。

3　計算結(jié)果及誤差討論

表1　計算結(jié)果數(shù)據(jù)表

從表1和圖3可以看出，當n≤10之前，Pcr的收斂速度很快，誤差可以控制在1.0%左右，但是隨著n值的增大，其收斂速度明顯放慢，而且還會出現(xiàn)數(shù)值上的波動。其主要原因是因為式(10)展開后得如下所示方程：

(12)

通過式(12)知：當n增大時，對a值的求解精度越來越高，但是當n>15時會出現(xiàn)誤差反而略有增大的現(xiàn)象(如當n=15時誤差為0.37%，但是當n=16時誤差反而增大到了0.39%)，其主要的原因是因為當n增大時，式(12)是一個高階的多項式方程，所以當多項式的階數(shù)增大時會出現(xiàn)高階震蕩現(xiàn)象。不過從計算的結(jié)果上來看，當n=15時，0.37%的誤差對于工程上的精度書完全能夠滿足，因此不建議求解n>15時的結(jié)果。由計算的數(shù)值說明該計算方法具有很高的精度，因此該方法具有很好的實際應用價值。

圖3　n值與對應的a值和Pcr關(guān)系圖

4　對于軸線為直線的變截面壓桿的討論

以變截面圓形壓桿為例略作討論，設(shè)其截面的直徑函數(shù)為：d=φ(x)，所以在xi處的慣性矩為：

(13)

將式(13)代入式(2)中，則

(14)

式(6)則變?yōu)槿缦拢?/p>

(15)

式(7)變?yōu)椋?/p>

(16)

將式(16)代入式(10)中，式(10)變?yōu)槿缦拢?/p>

(17)

由于a1,a2,…,an-1,an都是變量，所以需要將式(16)后面包含P的式子直接代入式(17)，然后直接求出Pcr的解即可。

用矩陣法求解壓桿穩(wěn)定的歐拉方程有一定的精確性，為解決變截面壓桿的穩(wěn)定性問題提供了一種數(shù)值求解的方法，并且為求解壓桿失穩(wěn)后的大撓度桿件變形問題提供了一個求解的方向。

[1]宣海洋.彈性介質(zhì)上等截面壓桿穩(wěn)定分析[J].山西建筑，2010,36(9)：45-46.

[2]洪振德. 變截面壓桿穩(wěn)定臨界力能量計算方法[J].江蘇建筑,2011,(3):28-30.

[3]熊容，蘇培東. 特定邊界條件下壓桿穩(wěn)定問題的解析解與數(shù)值解[J].四川建筑,2012,32(2):114-116.

[4]劉鴻文.材料力學(第三版·下冊)[M].北京:高等教育出版社，1992.

[5]劉鴻文.高等材料力學[M].北京:高等教育出版社，1985:237.

責任編輯王菊平

The matrix method for solving the finite differential form Euler equation of column stability

HU Zhang-yong

(College of Mechanical & Electrical Engineering, Huanggang Normal University, Huangzhou 438000, Hubei, China)

Euler equation of column stability can be expressed in the form of finite differential. According to the requirements of column stability, coefficient matrices the differential equations are solved with the method of matrices calculation. Subsequently, the critical stress of column stability of constant section of the simply supported beams is obtained. The results are compared with the analytical results. Also presented in the paper is an error analysis including the causes of the error and measures for diminishing it. The method proposed has positive value for the variable section beams.

column stability; finite differential form; Euler equation; matrices

TB12

A

1003-8078(2016)03-0076-04

2016-02-27

10.3969/j.issn.1003-8078.2016.03.19

胡章詠，男，湖北麻城人，講師，主要研究方向為力學和材料成型工藝與數(shù)值模擬等。