金 晶 邢譽峰
(北京航空航天大學 航空科學與工程學院,北京100191)
邊界元法作為現代數值方法之一,其基本思想是用邊界積分方程求解微分方程,但直至20世紀60年代人們才認識到它的價值[1-3].邊界元方法具有降低問題的維數、計算精度高等優點,其應用領域已遍及彈塑性力學、巖土力學、熱學、聲學、電磁學、生物細胞學等[4-9].近年來,隨著高性能大容量微機的出現,邊界元法廣泛應用于大型和非線性科學工程領域的數值計算,已成為一種強有力的數值計算方法.
梁、板作為常用的工程結構元件,探索其自由振動固有頻率的求解方法對認識和參數設計結構動態性能具有重要意義.文獻[10]用邊界方法(BEM,Boundary Element Method)分析了歐拉梁振動問題,包括歐拉梁縱向、扭轉、彎曲振動的頻率特性,提出了局部區間步長細化頻率掃描方法,但對精度和可行性未作深入分析.文獻[11-12]用靜態基本解分析了薄板的動態特性.靜態基本解表達式雖然簡單,但精度得不到保證.文獻[13]使用動態基本解分析了薄板的動態特性,但板的動態基本解比較復雜,在處理奇異積分時需要花費大量時間,給編程和數值計算帶來困難.
為了充分發揮邊界元方法本身高精度和高效率的優勢,利用桿、梁等一維結構固有振動的基本解來處理二維固有振動問題將是一條有效的途徑.截至目前,桿 (軸)和歐拉梁的基本解是已知的[14],而具有兩個廣義位移的剪切梁的基本解卻未見公開發表.基于這種情況,本文首先推導了鐵木辛柯梁的基本解,建立了相應的邊界積分方程.為了簡潔起見,以桿為例從理論上證明邊界元法處理一維問題時得到的頻率是精確的,分析了代數特征值求解方法和影響系數方法的特點.最后給出了剪切梁的數值算例.
基本解是邊界元方法的核心.不同問題基本解的求解方法也不同,對于一維問題,可用單位載荷方法或Duhamel積分;二維和三維彈性力學問題的基本解為Kelvin解;對于用Laplace和Poisson等方程描述的位勢問題,可應用格林第二公式來求解.此外還有通用的傅里葉變換和特征函數展開方法.下面用傅里葉變換方法推導鐵木辛柯梁的基本解.與歐拉梁不同的是,鐵木辛柯梁具有兩個獨立變量,這類似于二維問題.而與二維問題不同的是,鐵木辛柯梁的兩個獨立變量 (撓度和轉角)具有不同的量綱,又只是軸向坐標的函數.因此,為了求基本解,先對變量進行無量綱化.
鐵木辛柯梁固有振動的特征微分方程[14]為

式中,W為撓度函數;ψ為轉角函數;EI為彎曲剛度;ρ為密度;kGA 為剪切剛度;(),x和 (),xx分別表示對坐標x的一階微分和二階微分.

與之對應的鐵木辛柯梁基本解方程為

式中,δ(η-ζ)為狄拉克函數,ζ和η分別為測點和源點.
根據傅里葉變換,結合梁無窮遠處的邊界條件 (撓度、轉角、彎矩和剪力為0)可得撓度和轉角的基本解分別為

式中的λ1和λ2有下面兩種不同的形式:

而剪力和彎矩的基本解為

建立邊界積分方程是實現邊界元方法的關鍵步驟.利用基本解建立邊界積分方程的方法有加權殘量法、格林公式法和功的互等定理等.這里采用加權殘量法來建立鐵木辛柯梁的邊界積分方程.同樣也可以建立桿和歐拉梁的邊界積分方程.

對第1項分部積分兩次,第2項分部積分一次得


同理,對式 (4)加權得

用γ乘以式 (11),再把所得結果與式 (10)相加,并利用基本解式 (5)和式 (6)得

將式 (12)對ζ求一階導數得


式中系數矩陣分別為

如果梁兩端有平動彈簧和轉動彈簧約束,只要通過彈簧剛度系數建立剪力和彎矩與撓度和轉角關系,則可以正確模擬任意邊界條件.下面討論用代數特征值方程和影響系數法求解頻率時各自的特點.
一維問題邊界包括兩個端點,不存在邊界離散誤差,并且域內采用基本解,因此可以推斷用BEM求得的頻率一定是精確的.為了簡單起見,本節以均勻桿為例,來達到兩個目的:①證明用BEM求得的一維結構頻率是精確的;②討論代數特征值方法和影響系數法.本節的證明方法和頻率的解法同樣適用于梁.
與鐵木辛柯梁情況類似,可推導出桿的邊界方程式如下:

式中,U0和UL分別為桿兩端的位移;ε0和εL分別為兩端的應變;而基本解U*為

其中λ為頻率參數.
1)代數特征值方法
考慮兩端固支桿,約束條件為U0=0和UL=0.由式 (15)得頻率方程為

化簡為

同理可以得到兩端自由的頻率方程為

而左端固支右端自由的頻率方程為

從方程式 (18)~方程式 (20)可以看出,用代數特征值方法得到的三種不同邊界條件下的頻率方程中不但包括各自的精確頻率方程,也包括自由頻率方程.歐拉梁問題類似,即用這種方法得到的固支、自由、簡支、懸臂梁以及彈性支持的頻率方程中都包含自由梁的頻率方程.這應該起因于基本解的性質.值得指出的,對于不同邊界條件的鐵木辛柯梁,多出來的頻率并不是自由情況的頻率,這有待于深入研究.
2)影響系數方法
這種方法的核心是根據式 (15)得到用未知量表示已知量的關系式.以左端固支右端自由桿為例,其結果為

式中系數矩陣的對角元素為影響系數,其變量為頻率.當頻率為系統固有頻率時,影響系數為無窮大,據此可以確定系統的固有頻率.從式(21)容易看出,一端固支桿的頻率方程為cos λL=0.對于自由桿和固支桿情況,可以用類似的方法進行分析.
在進行數值計算時,需要對頻率進行掃描,其步長將決定這種方法的精度和效率.
由于鐵木辛柯梁的基本解比較復雜,直接從方程 (14)推導顯式形式的頻率方程比較困難.于是,下面根據影響系數法求解其固有振動頻率,并與精確解[15]和NASTRAN結果 (用50個BEAM單元)進行比較,見表1.梁長0.05 m,矩形截面高0.01 m,寬0.01 m,彈性模量7×1010N/m2,剪切系數5/6,泊松比0.31,質量密度2700 kg/m3.

表1 固支梁的頻率 Hz
表1中BEM的結果是用掃描步長1 Hz得到的,如此得到的近似頻率和精確解之間的絕對誤差不超過1.圖1給出了頻率和影響系數的關系曲線.由于NASTRAN中的BEAM單元不考慮轉動慣量的影響,因此對應的每一階頻率都偏大.
當改變均勻梁的截面形狀時,只是改變了截面慣性矩.因此,用本文方法可以給出任意截面形式且具有任意邊界條件梁的精確頻率方程.

圖1 固支梁的頻率-影響系數關系曲線
基于理論推導和數值分析,可得出如下結論:
1)對于一維均勻結構,沒有邊界離散誤差,用邊界元法得到的頻率是精確解.
2)采用代數特征值方法求解頻率時,由于基本解本身的特性,對于不同結構不同邊界條件求出的頻率都包含多余的頻率.
3)采用影響系數方法求固有頻率,避免了多余頻率的出現,但其精度和效率取決于頻率掃描步長.
4)根據已有的桿、歐拉梁的基本解和本文推導的剪切梁的基本解,可以簡化用邊界元法分析二維結構固有振動特性的過程.這是下一步要開展的研究工作.
References)
[1]Brebbia C A,Walker S.Boundary element techniques in engineering[M].London:Newnes-Butterworths,1980
[2]Brebbia C A,Telles J C F,Wrobel L C.Boundary element techniques-theory and applications in engineering[M].New York:Springer-Verlag,1984
[3]楊德全,趙忠生.邊界元理論及應用[M].北京:北京理工大學出版社,2002
Yang Dequan,Zhao Zhongsheng.The theory and application of boundary element method[M].Beijing:Beijing Institute of Technology Press,2002(in Chinese)
[4]Wei X,Leu L J,Chandra A.Shape optimization in elasticity and elasto-viscoplasticity by the boundary element method [J].International Journal of Solids and Structures,1994,31(4):533-550
[5]Chen Chao-shi,Pan ernian,Amadei Bernard.Fracture mechanics analysis of cracked discs of anisotropic rock using the boundary element method[J].International Journal of Rock Mechanics and Mining Sciences,1998,35(2):195-218
[6]Bialecki R A.Solving heat radiation problems using the boundary element method[M].Boston:Computational Mechanics Publications,1993
[7]陳劍,高煜,張永斌,等.基于有限元-邊界元的聲學構形靈敏度分析[J].振動工程學報,2009,22(2):213-217
Chen Jian,Gao Yu,Zhang Yongbin,et al.Acoustic configuration sensitivity analysis based on FEM and BEM[J].Journal of Vibration Engineering,2009,22(2):213-217(in Chinese)
[8]林慶華,栗保明.有限元/邊界元耦合法計算電磁軌道炮三維瞬態渦流場 [J].南京理工大學學報:自然科學版,2010,34(2):218-221
Lin Qinghua,Li Baoming.Finite-element/boundary-element coupling method for 3D transient eddy current field calculation in electrom agnetic railgun[J].Journal of Nanjing University of Science and Technology:Natural Science,2010,34(2):218-221(in Chinese)
[9]彭紅梅,楊德全.雙分叉動脈血流動力學特性的邊界元分析 [J].醫用生物力學,2010,25(4):283-287
Peng Hongmei,Yang Dequan.Boundary element analysis on characteristics of double bifurcation arterial hemo-dynamics[J].Journal of Medical Biomechanics,2010,25(4):283 -287(in Chinese)
[10]黃玉美,張廣鵬,高峰.虛擬樣機整機結構特性邊界元法仿真[M].北京:機械工業出版社,2004
Huang Yumei,Zhang Guangpeng,Gao Feng.The simulation with boundary element method of the whole structural characteristics for the virtual prototype[M].Beijing:China Machine Press,2004(in Chinese)
[11]Jesus A,Costa J.Plate vibrations using BEM[J].Application Mathematics Modelling,1988,12:78-85
[12]Tanaka M,Matsumoto T,Shiozaki A.Application of boundary-domain element method to the free vibration problem of plate structures[J].Computers & Structures,1998,66(6):725-735
[13]Providakis C P,Beskos D E.Free and forced vibrations of plates by boundary elements[J].Computer Methods in Applied Mechanics and Engineering,1989,74(3):231:250
[14]邢譽峰,李敏.計算固體力學原理與方法 [M].北京:北京航空航天大學出版社,2011
Xing Yufeng,Li Min.The theory and method of computational solid mechanics[M].Beijing:Beijing University of Aeronautics and Astronautics Press,2011(in Chinese)
[15]邢譽峰,李敏.工程振動基礎 [M].北京:北京航空航天大學出版社,2011
Xing Yufeng,Li Min.The foundation of engineering vibration[M].Beijing:Beijing University of Aeronautics and Astronautics Press,2011(in Chinese)