重心插值Galerkin 法求解梁彎曲變形問題

綦甲帥，王兆清

(山東建筑大學工程力學研究所，山東 濟南 250101)

在處理梁彎曲問題時，經(jīng)常遇到非連續(xù)載荷作用，解析方法需要在間斷處分段，在每一段分別處理，然后在分段處通過施加連續(xù)性條件和邊界條件來求解，導致積分常數(shù)很多，計算過程復雜。在這種情況下，常采用數(shù)值方法求解［1-3］。利用廣義函數(shù)來表示非連續(xù)載荷的分布集度，不需要分段，可在全梁范圍內建立統(tǒng)一的梁彎曲變形控制方程，使計算過程大為簡化［4］。

最常用的廣義函數(shù)是Delta 函數(shù)。這種函數(shù)最早出現(xiàn)于20 世紀30年代，是由著名物理學家Dirac［5］在量子力學研究中引入和定義的，后來被命名為Delta 函數(shù)。50年代Schwarz［6］在深入研究Delta 函數(shù)性質的基礎上創(chuàng)立了分布論(亦稱廣義函數(shù)論)。2000年，Yavari 等［7-8］研究了廣義函數(shù)在梁彎曲問題上的應用。Falsone［9］給出了利用廣義函數(shù)Delta 函數(shù)分析不連續(xù)梁彎曲問題的解析方法。

有理函數(shù)插值有時比多項式插值具有更好的精度。在經(jīng)典的有理插值中，對于給定的n +1 個節(jié)點，人們尋求如pM/qN的有理函數(shù)插值，其中pM、qN分別為次數(shù)不超過M、N 的多項式，且M+N=n。這種有理函數(shù)的主要缺點是，在插值區(qū)間內無法控制極點的產(chǎn)生。Berrut 等［10］提出利用高次多項式來構造有理函數(shù)插值，這樣可以避免插值區(qū)間內極點的產(chǎn)生。2005年，Berrut 等［11］對重心型有理函數(shù)插值作了總結，利用n 次多項式構造的有理函數(shù)插值，可以寫成重心形式。2006年Floater 等［12］提出一類重心型有理函數(shù)插值，該插值在整個數(shù)軸上都不存在極點，插值函數(shù)具有無窮次光滑性?；谖⒎址匠倘跣问降腉alerkin 法［13］，可以降低近似函數(shù)的連續(xù)性要求，并且利用Delta 函數(shù)關于積分的篩選性，Galerkin 法可以消除求解問題中的奇異項［14］。

本文利用重心有理插值函數(shù)作為試函數(shù)，運用廣義函數(shù)建立梁彎曲變形的控制方程，利用Delta 函數(shù)的積分篩選性，提出求解梁彎曲變形問題的重心插值Galerkin 法。

1 廣義函數(shù)定義

(1)Heaviside 函數(shù)定義為:

(2)Delta 函數(shù)也稱為脈沖函數(shù)，其一般定義為:

2 重心有理插值Galerkin 法

以重心有理插值作為試函數(shù)，采用Galerkin 法建立數(shù)值求解梁彎曲變形控制方程的公式。對于在集度為q(x)的荷載作用下的細長梁，在彈性范圍內發(fā)生線性彎曲，其控制方程為:

采用加權殘數(shù)法建立方程(1)的積分弱形式，取加權函數(shù)W，有

對方程左邊分部積分二次，可得

式中，W 為加權函數(shù)。在Galerkin 法中，權函數(shù)取為插值的基函數(shù)。

對于求解區(qū)間的一組節(jié)點0=x1＜x2＜… ＜xn=l，通過重心有理插值，梁的撓度可寫成

式中，φk(x)為重心有理插值基函數(shù)。

Jk={i ∈I:k－ d ≤i ≤k}，I={0，1，2，…，n}，d 為重心有理插值參數(shù)，滿足0 ≤d ≤n。

在Galerkin 過程中，令Wk=φk，將φk(x)代入(3)式得到

其中φ(x)=［φ0(x)，φ1(x)，…，φn(x)］T，u=(u0，u1，…，un)T。方程(7)寫成矩陣的形式為

通過施加邊界條件，可以消除向量b 的未知導數(shù)。由于插值基函數(shù)為有理函數(shù)形式，采用數(shù)值積分計算矩陣A 和向量b。

3 數(shù)值算例

本文算例的計算程序采用MATLAB 編寫，除了主程序外，還有計算基函數(shù)及其一階導數(shù)、二階導數(shù)的子函數(shù)和施加邊界條件的子函數(shù)。

圖1 梁受集中力載荷Fig.1 An beam with concentrated load

算例1 兩端簡支梁，長L=50 cm，高2 cm，寬1 cm，彈性模量為10 GPa，受如圖1 所示集中力F=50 N，集中力作用于梁跨中，運用廣義函數(shù)表示的載荷集度q(x)=Fδ(x-x0)，x0為集中力作用點。計算節(jié)點采用等距節(jié)點，Galerkin 法中的積分采用小區(qū)間［xi，xi+1］上的3 點Gauss數(shù)值積分計算，插值參數(shù)d=4。梁各點的相對誤差和絕對誤差分別定義為Er=|uc-ue|/|ue|，Ea=|uc-ue|。不同數(shù)量節(jié)點計算的相對誤差和絕對誤差分別定義為Er=‖uc-ue‖2/‖ue‖2，Ea=‖uc-ue‖2，其中uc，ue分別為函數(shù)的數(shù)值計算值和解析解值列向量，‖·‖2為向量的2 范數(shù)。數(shù)值計算得到梁軸線的各節(jié)點撓度、轉角、彎矩與解析解比較如表1 所示(選取30 個節(jié)點計算，表中只列舉了有代表性的11 個節(jié)點)，不同數(shù)量節(jié)點計算誤差如表2 所示，解析解見文獻［15］。

表1 算例1 梁各點撓度、轉角和彎矩解析解與本文解比較Table 1 The comparion of analytical solution and the presented solution of beam deflection，angle and moment of each point for instance 1

算例2 受如圖2 所示荷載，q=200 N/m，其他參數(shù)同算例1，運用廣義函數(shù)表示的載荷集度q(x)=qH(x-xi)，得到梁軸線的各節(jié)點撓度、轉角、彎矩與解析解比較如表3 所示(選取30 個節(jié)點計算，表中只列舉了有代表性的11 個節(jié)點)，不同數(shù)量節(jié)點計算誤差如表4 所示。

表2 算例1 不同數(shù)量節(jié)點的相對誤差和絕對誤差Table 2 Relative and absdute errors of different number of nodes for instance 1

圖2 梁受部分均布載荷Fig.2 An beam with partial uniform load

表3 算例2 梁各點撓度、轉角和彎矩解析解與本文解比較Table 3 The comparion of analytical solution and the presented solution of beam deflection，angle and moment of each point for instance 2

表4 算例2 不同數(shù)量節(jié)點的計算相對誤差和絕對誤差Table 4 Relative and absolate errors of different number of nodes for instance 2

算例3 受如圖3 所示部分三角形荷載，q=200 N/m，其他參數(shù)同算例1，運用廣義函數(shù)表示的載荷集度q(x)=qxH(x -xi)/L，得到梁軸線的各節(jié)點撓度、轉角、彎矩與解析解比較如表5 所示(選取30 個節(jié)點計算，表中只列舉了有代表性的11 個節(jié)點)，不同數(shù)量節(jié)點計算誤差如表6 所示。

圖3 梁受部分三角形載荷Fig.3 An beam with partial triangular load

表5 算例3 梁各點撓度、轉角和彎矩解析解與本文解比較Table 5 The comparion of analytical solution and the presented solution of beam deflection，angle and moment of each point for instance 3

表6 算例3 不同數(shù)量節(jié)點計算相對誤差Table 6 Relative and absolate errors of different number of nodes for instance 3

分析上述3 個數(shù)值算例的結果，可以得出以下結論:(1)重心有理插值Galerkin 法有很好的精度，并且計算精度隨計算節(jié)點數(shù)量的增加而提高;(2)由于求導過程中喪失精度，因此轉角和彎矩的精度比撓度的精度差一些;(3)算例1 中的計算精度高于算例2、3，是因為算例1 中的未知函數(shù)比算例2、3 中的光滑性好。

4 結論

重心有理插值Galerkin 法可以快速的求解梁彎曲變形問題，并得到具有無窮次光滑性質的數(shù)值解。對于非連續(xù)載荷等復雜載荷的計算，利用廣義函數(shù)以及Delta 函數(shù)的積分篩選性建立梁彎曲變形的控制方程，不需要分段，大大簡化了求解過程，并且少量的節(jié)點即可得到較高的精度，為這類問題提供了新的方法和思路。該方法原理簡單，易于程序實現(xiàn)，在力學和結構分析中有良好的應用前景。

［1］VIRDI K S.Finite difference method for nonliear analysis of structures［J］.Journal of Constructional Steel Research，2006，62(11):1210 -1218.

［2］RAJU I S，PHILLIPS D R，KRISHNAMURTHY T.A Meshless Method Using Radial Basis Functions for Beam Bending Problems［EB/OL］.［2012 -09 -11］http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20040191554_2004198262.pdf.

［3］袁玉全，彭建設.復雜載荷下梁彎曲問題的微分求積法應用研究［J］.四川理工學院學報:自然科學版，2006，19(5):81 -84.

［4］曾紀鵬.用奇異函數(shù)建立變參數(shù)連續(xù)梁撓度統(tǒng)一方程的新方法［J］.四川建筑，2011，1(1):128 -130.

［5］DIRAC P A M.The principles of quantum mechanics［M］.Oxford:Oxford University Press，1930.

［6］SCHWARZ L.Theory of distributions［M］.Paris:Hermann，1966.

［7］YAVARI A，SARKANI S，MOYER E T.On applications of generalized functions to beam bending problems［J］.International Journal of Solids and Structures，2000，37(40):5675 -5705.

［8］YAVARI A，SARKANI S.On applications of generalized functions to the analysis of Euler-Bernoulli beam-columns with jump discontinuities［J］.International Journal of Mechanical Sciences，2001，43(6):1543 -1562.

［9］FALSONE G.The use of generalised functions in the discontinuous beam bending differential equations［J］.International Journal of Engineering Education，2002，18(3):337 -343.

［10］BERRUT J P，MITTELMANN H D.Lebesgue constant minimizing linear rational interpolation of continuous functions over the interval［J］.Computational Mathematics and Applications，1997，33(6):77 -86.

［11］BERRUT J P，BALTENSPERGER R，MITTELMANN H D.Recent developments in barycentric rational interpolation［J］.International Series of Numerical Mathematics，2005，151:27 -51.

［12］FLOATER M S，HORMANN K.Barycentric Rational Interpolation with no poles and high rates of approximation［J］.Numerische Mathematik，2007，107(2):315 -331.

［13］王兆清，李樹忱，唐炳濤，等.求解兩點邊值問題的有理插值Galerkin 法［J］.山東建筑大學學報，2008，23(4):283 -286.

［14］王兆清，綦甲帥，唐炳濤.奇異源項問題的重心插值數(shù)值解［J］.計算物理，2011，28(6):883 -888.

［15］劉鴻文.材料力學［M］.4 版.北京:高等教育出版社，2004.