綦甲帥，王兆清

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

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

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

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

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

1 廣義函數(shù)定義

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

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

2 重心有理插值Galerkin 法

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

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

對(duì)方程左邊分部積分二次，可得

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

對(duì)于求解區(qū)間的一組節(jié)點(diǎn)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 的未知導(dǎo)數(shù)。由于插值基函數(shù)為有理函數(shù)形式，采用數(shù)值積分計(jì)算矩陣A 和向量b。

3 數(shù)值算例

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

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

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

表1 算例1 梁各點(diǎn)撓度、轉(zhuǎn)角和彎矩解析解與本文解比較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，運(yùn)用廣義函數(shù)表示的載荷集度q(x)=qH(x-xi)，得到梁軸線的各節(jié)點(diǎn)撓度、轉(zhuǎn)角、彎矩與解析解比較如表3 所示(選取30 個(gè)節(jié)點(diǎn)計(jì)算，表中只列舉了有代表性的11 個(gè)節(jié)點(diǎn))，不同數(shù)量節(jié)點(diǎn)計(jì)算誤差如表4 所示。

表2 算例1 不同數(shù)量節(jié)點(diǎn)的相對(duì)誤差和絕對(duì)誤差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 梁各點(diǎn)撓度、轉(zhuǎn)角和彎矩解析解與本文解比較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é)點(diǎn)的計(jì)算相對(duì)誤差和絕對(duì)誤差Table 4 Relative and absolate errors of different number of nodes for instance 2

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

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

表5 算例3 梁各點(diǎn)撓度、轉(zhuǎn)角和彎矩解析解與本文解比較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é)點(diǎn)計(jì)算相對(duì)誤差Table 6 Relative and absolate errors of different number of nodes for instance 3

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

4 結(jié)論

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

［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］袁玉全，彭建設(shè).復(fù)雜載荷下梁彎曲問題的微分求積法應(yīng)用研究［J］.四川理工學(xué)院學(xué)報(bào):自然科學(xué)版，2006，19(5):81 -84.

［4］曾紀(jì)鵬.用奇異函數(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］王兆清，李樹忱，唐炳濤，等.求解兩點(diǎn)邊值問題的有理插值Galerkin 法［J］.山東建筑大學(xué)學(xué)報(bào)，2008，23(4):283 -286.

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

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