于衛(wèi)濤，宋潔，趙維霞

(1.濱州市建設工程質量監(jiān)督站，山東 濱州 256600;2.濱州市工程建設標準定額站，山東 濱州 256600;3.山東城市建設職業(yè)學院，山東 濟南 250014)

0 引言

當?shù)戎睏U件受均布軸向荷載作用時，采用靜力平衡準則判斷原始直線狀態(tài)的平衡穩(wěn)定性，首先應給予桿件微小偏移，使之呈現(xiàn)彎曲狀態(tài)，看其桿件是否能平衡。如果存在平衡狀態(tài)，則其處在平衡路徑上的分岔點，也就是臨界狀態(tài)。此時的荷載即為臨界荷載。至于臨界狀態(tài)以后的情況，若只采用線性理論則不能給出合理的說明。對于細長壓桿，當載荷超過臨界值時，壓力與彎曲變形還存在一定的關系，這就是壓桿的后屈曲性質。細長桿的后屈曲屬于幾何非線性問題。

在工程應用中，軸向均布荷載下壓桿穩(wěn)定問題可以歸結為微分方程的邊值問題，很難得到問題的解析解，一般采用數(shù)值方法求解。數(shù)值求解方法主要有有限差分法和有限元方法[1]等，這兩種方法要對求解區(qū)域劃分單元，計算精度依賴于單元的大小。采用配點法求解邊值問題不需要劃分單元，公式簡單，不需要積分，易于編程。目前用于求解工程中的常微分方程邊值問題的配點法主要有擬譜法和微分求積法。擬譜法[2－3]是根據(jù)譜方法發(fā)展出來的一種方法，雖然這種方法的理論研究已經(jīng)有進一步的發(fā)展，但是工程技術人員對這種方法不是很了解。微分求積法[4－5]的基本原理是將未知函數(shù)在區(qū)間上所有離散點的函數(shù)值的加權和來逼近該函數(shù)在某一離散點的偏導數(shù)或者積分，這種方法中的權系數(shù)的確定通常是根據(jù)Lagrange多項式在網(wǎng)點處的導數(shù)值給出。但是這種方法的局限性是離散點不能取得太多，否則Lagrange多項式表示的曲線隨多項式次數(shù)的升高而出現(xiàn)Runge現(xiàn)象，從而產(chǎn)生計算的不穩(wěn)定性[6]。重心 Lagrange插值具有極好的數(shù)值穩(wěn)定性[6]和極高的近似精度[7]，同時重心 Lagrange 插值公式具有緊湊的各階導數(shù)的計算公式[8]。

本文所采用的重心插值配點法就是用重心Lagrange插值多項式求出某一函數(shù)在各個離散點的微分矩陣，從而可以通過矩陣的運算來求出微分方程的解。而許多結構力學問題的求解最終也歸結為在一定的邊界條件和初始條件下的(偏)微分方程(組)的求解，所以，可以把重心插值配點法作為用于結構力學問題的求解的一種數(shù)值方法。本文采用重心插值配點法求得了壓桿屈曲的后屈曲撓度數(shù)值和屈曲路徑分岔點的臨界載荷值，精度較高，得到了滿意的效果。

1 基于重心插值的微分矩陣

設g(x)為一待求函數(shù)，它在區(qū)間[0，L]上連續(xù)可微，現(xiàn)沿x軸設置N+1個節(jié)點，并以節(jié)點函數(shù)值g(xi)(j=0，1，…，N)作為基本未知量，且令 g(xj)=vj。在全域內采用重心Lagrange插值逼近g(x)，則有:

其中:mj為重心權。

因為所求解的問題是線性的，因此可以定義一個(N+1)×(N+1)的矩陣D。則(3)式可以寫成如下的形式:

式中:w為未知函數(shù)一階導數(shù)值構成的列向量;v為節(jié)點函數(shù)值的列向量;D就是對g(x)求一次導數(shù)的微分矩陣。由(1)、(2)、(3)式可得

2 軸向均布壓力下細長壓桿的重心插值配點法

細長桿在坐標原點處固定，頂端自由(圖1)，桿件受到軸向均布荷載作用。s為以自由端為起始點沿桿方向的任意長度，θ為與過桿上一點的切線與與垂直線所成的角度，h為桿自由端x向位移，δ為桿自由端向位移。從載荷施加的初期到屈曲產(chǎn)生這一段時間桿件是維持豎直的，在臨界載荷處發(fā)生分岔，屈曲產(chǎn)生之后，桿件發(fā)生變形撓曲，求解屈曲路徑為大撓度問題。

假定桿件的材料是理想彈性，其應力應變關系符合胡克定律，變形和位移只要由彎曲產(chǎn)生，即略去剪切變形以及軸線的伸縮變形，桿的長度保持不變，并且假定橫截面保持為平面。

圖1 受軸向均布荷載細長壓桿屈曲圖

由文獻[9]知，大撓度下壓桿的變形方程無量綱化后為:

當壓桿受軸向力作用時，到發(fā)生屈曲之前都是保持直線狀態(tài)的，所以，在分岔點的研究上，可以近似的認為 sinθ≈θ，即:

桿件的變形方程為

式中:θ為N+1的轉角位移參數(shù)列陣;D2是N+1行N+1列的2階微分矩陣。

令 K=D2，Q=s×IN+1，則上式可以寫成

其中:IN+1，為N+1行N+1列的單位陣。

則邊界約束條件可以表示為:

將此2個方程的左邊分別取代矩陣K的第1，N+1行，并把矩陣K和Q第N+1行放到矩陣的第1行和第2行之間，則原矩陣的第2行，第3行，…，第N行依次成為第3行，第4行，…，第N+1行。同樣，把第N+1列放到矩陣的第1列和第2列之間，則原矩陣的第2列，第3列，…，第N列依次成為第3列，第4列，…，第N+1列。則K和Q變?yōu)?

用Matlab命令求矩陣Q－1K的特征值和特征向量就可得到屈曲臨界值。一般的矩陣Q－1K的特征值有多個，將這些特征值按升序的方法排列，從中取出前幾個特征值以及與它們相對應的特征向量。

微分矩陣以及梁的屈曲臨界值和屈曲模態(tài)的計算，本文都編制Matlab程序來實現(xiàn)。解得k值與文獻[9]對比如表1。

表1 不同取點對應的值

圖2 一階屈曲模態(tài)

可以看出，重心插值配點法比微分求積法收斂更快，且精度較高。

3 結論

對于屈曲臨界值的求解，文獻[10]用了多種方法，其中能量法應用的討論非常詳細，當取撓曲線方程為y=fsin[πx/(2l)]時，得出臨界值為8.27;當取撓曲線方程為y=f(1－cos[πx/(2l)])時，得出的臨界值為7.888;當取撓曲線方程為y=f1(1－cos[πx/(2l)])+f3(1 －cos[πx/(2l)])時，得出的臨界值為7.84。即當選取的撓曲線月接近實際撓曲形狀，得出的結果越精確。本文獻中用靜力平衡法得出的解為7.837，文獻[11]由線性理論求解得出的臨界值為7.837，文獻[12]得出的臨界值為7.83，本文也得出屈曲的臨界值，解的前四位有效數(shù)值為7.837，與參考文獻完全吻合，結果是可信的。

重心插值配點法是一種快速求解偏微分方程的方法。從算例可以看出，這種方法的精度很高。有限差分法和有限元法需要對求解區(qū)域劃分單元，求解的精度取決于劃分單元的大小，從而導致大的計算量。微分求積法的插值點不能取得太多，否則會出現(xiàn)插值數(shù)值不穩(wěn)定性。本文采用的重心插值配點法不但可以避免上述方法的缺陷，而且可得到相當?shù)木取?/p>

[1]王勖成，邵敏.有限單元法基本原理及數(shù)值方法[M].2版，北京:清華大學出版社，2000.

[2]楊繼明.熱傳導方程初邊值問題的譜方法[J].湖南工程學院學報，2007，17(2):71 －73.

[3]聶國雋，仲政.微分求積單元法在結構工程中的應用[J].力學季刊，2005，26(3):423 －427.

[4]王維蘭.擬lagrange多項式[J].西北民族學院學報:自然科學版，1999，20(4):15 －23.

[5]王兆清，李淑萍，唐炳濤.任意連續(xù)函數(shù)的多項式插值逼近[J].山東建筑大學學報，2007，22(2):158 －162.

[6]SCHNEIDER C，WERNER W.Some new aspects of rational interpolation[J].Math Comp，1986，47(175):285 －299.

[7]聶毓琴，孟廣偉.材料力學[M].北京:機械工業(yè)出版社，2004.

[8]SHU C，DU H.Implementation of clamped and simply supported doundary conditions in the GDQ free vibration analysis of beams and plates[J].Int J Solids Structures，1997，34(7):819 －835.

[9]劉洋，楊永波，梁樞平.軸向均布載荷下壓桿穩(wěn)定問題的DQ解[J].力學與實踐，2005，27(2):44 －47.

[10]吳明德.彈性桿件穩(wěn)定理論[M].北京:高等教育出版社，1988.

[11]HOLDEN J T.On the finite defledtions of thin beams[J].Int J Solids Stuctures，1972(8):1051 －1055.

[12]TIMOSHENKO S P，GERE J M.Thery of elastic stability[M].MCGraw-Hill，1961.