金 晶 邢譽(yù)峰

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京100191)

邊界元法作為現(xiàn)代數(shù)值方法之一，其基本思想是用邊界積分方程求解微分方程，但直至20世紀(jì)60年代人們才認(rèn)識(shí)到它的價(jià)值［1－3］.邊界元方法具有降低問題的維數(shù)、計(jì)算精度高等優(yōu)點(diǎn)，其應(yīng)用領(lǐng)域已遍及彈塑性力學(xué)、巖土力學(xué)、熱學(xué)、聲學(xué)、電磁學(xué)、生物細(xì)胞學(xué)等［4－9］.近年來，隨著高性能大容量微機(jī)的出現(xiàn)，邊界元法廣泛應(yīng)用于大型和非線性科學(xué)工程領(lǐng)域的數(shù)值計(jì)算，已成為一種強(qiáng)有力的數(shù)值計(jì)算方法.

梁、板作為常用的工程結(jié)構(gòu)元件，探索其自由振動(dòng)固有頻率的求解方法對(duì)認(rèn)識(shí)和參數(shù)設(shè)計(jì)結(jié)構(gòu)動(dòng)態(tài)性能具有重要意義.文獻(xiàn)［10］用邊界方法(BEM，Boundary Element Method)分析了歐拉梁振動(dòng)問題，包括歐拉梁縱向、扭轉(zhuǎn)、彎曲振動(dòng)的頻率特性，提出了局部區(qū)間步長細(xì)化頻率掃描方法，但對(duì)精度和可行性未作深入分析.文獻(xiàn)［11－12］用靜態(tài)基本解分析了薄板的動(dòng)態(tài)特性.靜態(tài)基本解表達(dá)式雖然簡單，但精度得不到保證.文獻(xiàn)［13］使用動(dòng)態(tài)基本解分析了薄板的動(dòng)態(tài)特性，但板的動(dòng)態(tài)基本解比較復(fù)雜，在處理奇異積分時(shí)需要花費(fèi)大量時(shí)間，給編程和數(shù)值計(jì)算帶來困難.

為了充分發(fā)揮邊界元方法本身高精度和高效率的優(yōu)勢(shì)，利用桿、梁等一維結(jié)構(gòu)固有振動(dòng)的基本解來處理二維固有振動(dòng)問題將是一條有效的途徑.截至目前，桿 (軸)和歐拉梁的基本解是已知的［14］，而具有兩個(gè)廣義位移的剪切梁的基本解卻未見公開發(fā)表.基于這種情況，本文首先推導(dǎo)了鐵木辛柯梁的基本解，建立了相應(yīng)的邊界積分方程.為了簡潔起見，以桿為例從理論上證明邊界元法處理一維問題時(shí)得到的頻率是精確的，分析了代數(shù)特征值求解方法和影響系數(shù)方法的特點(diǎn).最后給出了剪切梁的數(shù)值算例.

1 鐵木辛柯梁自由振動(dòng)的邊界方程

1.1 基 本 解

基本解是邊界元方法的核心.不同問題基本解的求解方法也不同，對(duì)于一維問題，可用單位載荷方法或Duhamel積分;二維和三維彈性力學(xué)問題的基本解為Kelvin解;對(duì)于用Laplace和Poisson等方程描述的位勢(shì)問題，可應(yīng)用格林第二公式來求解.此外還有通用的傅里葉變換和特征函數(shù)展開方法.下面用傅里葉變換方法推導(dǎo)鐵木辛柯梁的基本解.與歐拉梁不同的是，鐵木辛柯梁具有兩個(gè)獨(dú)立變量，這類似于二維問題.而與二維問題不同的是，鐵木辛柯梁的兩個(gè)獨(dú)立變量 (撓度和轉(zhuǎn)角)具有不同的量綱，又只是軸向坐標(biāo)的函數(shù).因此，為了求基本解，先對(duì)變量進(jìn)行無量綱化.

鐵木辛柯梁固有振動(dòng)的特征微分方程［14］為

式中，W為撓度函數(shù);ψ為轉(zhuǎn)角函數(shù);EI為彎曲剛度;ρ為密度;kGA 為剪切剛度;()，x和 ()，xx分別表示對(duì)坐標(biāo)x的一階微分和二階微分.

與之對(duì)應(yīng)的鐵木辛柯梁基本解方程為

式中，δ(η－ζ)為狄拉克函數(shù)，ζ和η分別為測點(diǎn)和源點(diǎn).

根據(jù)傅里葉變換，結(jié)合梁無窮遠(yuǎn)處的邊界條件 (撓度、轉(zhuǎn)角、彎矩和剪力為0)可得撓度和轉(zhuǎn)角的基本解分別為

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

而剪力和彎矩的基本解為

1.2 邊界積分方程

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

對(duì)第1項(xiàng)分部積分兩次，第2項(xiàng)分部積分一次得

同理，對(duì)式 (4)加權(quán)得

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

將式 (12)對(duì)ζ求一階導(dǎo)數(shù)得

式中系數(shù)矩陣分別為

如果梁兩端有平動(dòng)彈簧和轉(zhuǎn)動(dòng)彈簧約束，只要通過彈簧剛度系數(shù)建立剪力和彎矩與撓度和轉(zhuǎn)角關(guān)系，則可以正確模擬任意邊界條件.下面討論用代數(shù)特征值方程和影響系數(shù)法求解頻率時(shí)各自的特點(diǎn).

2 固有振動(dòng)邊界方程的解法

一維問題邊界包括兩個(gè)端點(diǎn)，不存在邊界離散誤差，并且域內(nèi)采用基本解，因此可以推斷用BEM求得的頻率一定是精確的.為了簡單起見，本節(jié)以均勻桿為例，來達(dá)到兩個(gè)目的:①證明用BEM求得的一維結(jié)構(gòu)頻率是精確的;②討論代數(shù)特征值方法和影響系數(shù)法.本節(jié)的證明方法和頻率的解法同樣適用于梁.

與鐵木辛柯梁情況類似，可推導(dǎo)出桿的邊界方程式如下:

式中，U0和UL分別為桿兩端的位移;ε0和εL分別為兩端的應(yīng)變;而基本解U*為

其中λ為頻率參數(shù).

1)代數(shù)特征值方法

考慮兩端固支桿，約束條件為U0=0和UL=0.由式 (15)得頻率方程為

化簡為

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

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

從方程式 (18)～方程式 (20)可以看出，用代數(shù)特征值方法得到的三種不同邊界條件下的頻率方程中不但包括各自的精確頻率方程，也包括自由頻率方程.歐拉梁問題類似，即用這種方法得到的固支、自由、簡支、懸臂梁以及彈性支持的頻率方程中都包含自由梁的頻率方程.這應(yīng)該起因于基本解的性質(zhì).值得指出的，對(duì)于不同邊界條件的鐵木辛柯梁，多出來的頻率并不是自由情況的頻率，這有待于深入研究.

2)影響系數(shù)方法

這種方法的核心是根據(jù)式 (15)得到用未知量表示已知量的關(guān)系式.以左端固支右端自由桿為例，其結(jié)果為

式中系數(shù)矩陣的對(duì)角元素為影響系數(shù)，其變量為頻率.當(dāng)頻率為系統(tǒng)固有頻率時(shí)，影響系數(shù)為無窮大，據(jù)此可以確定系統(tǒng)的固有頻率.從式(21)容易看出，一端固支桿的頻率方程為cos λL=0.對(duì)于自由桿和固支桿情況，可以用類似的方法進(jìn)行分析.

在進(jìn)行數(shù)值計(jì)算時(shí)，需要對(duì)頻率進(jìn)行掃描，其步長將決定這種方法的精度和效率.

3 數(shù)值分析

由于鐵木辛柯梁的基本解比較復(fù)雜，直接從方程 (14)推導(dǎo)顯式形式的頻率方程比較困難.于是，下面根據(jù)影響系數(shù)法求解其固有振動(dòng)頻率，并與精確解［15］和NASTRAN結(jié)果 (用50個(gè)BEAM單元)進(jìn)行比較，見表1.梁長0.05 m，矩形截面高0.01 m，寬0.01 m，彈性模量7×1010N/m2，剪切系數(shù)5/6，泊松比0.31，質(zhì)量密度2700 kg/m3.

表1 固支梁的頻率 Hz

表1中BEM的結(jié)果是用掃描步長1 Hz得到的，如此得到的近似頻率和精確解之間的絕對(duì)誤差不超過1.圖1給出了頻率和影響系數(shù)的關(guān)系曲線.由于NASTRAN中的BEAM單元不考慮轉(zhuǎn)動(dòng)慣量的影響，因此對(duì)應(yīng)的每一階頻率都偏大.

當(dāng)改變均勻梁的截面形狀時(shí)，只是改變了截面慣性矩.因此，用本文方法可以給出任意截面形式且具有任意邊界條件梁的精確頻率方程.

圖1 固支梁的頻率－影響系數(shù)關(guān)系曲線

4 結(jié)論

基于理論推導(dǎo)和數(shù)值分析，可得出如下結(jié)論:

1)對(duì)于一維均勻結(jié)構(gòu)，沒有邊界離散誤差，用邊界元法得到的頻率是精確解.

2)采用代數(shù)特征值方法求解頻率時(shí)，由于基本解本身的特性，對(duì)于不同結(jié)構(gòu)不同邊界條件求出的頻率都包含多余的頻率.

3)采用影響系數(shù)方法求固有頻率，避免了多余頻率的出現(xiàn)，但其精度和效率取決于頻率掃描步長.

4)根據(jù)已有的桿、歐拉梁的基本解和本文推導(dǎo)的剪切梁的基本解，可以簡化用邊界元法分析二維結(jié)構(gòu)固有振動(dòng)特性的過程.這是下一步要開展的研究工作.

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］楊德全，趙忠生.邊界元理論及應(yīng)用［M］.北京:北京理工大學(xué)出版社，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］陳劍，高煜，張永斌，等.基于有限元-邊界元的聲學(xué)構(gòu)形靈敏度分析［J］.振動(dòng)工程學(xué)報(bào)，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ì)算電磁軌道炮三維瞬態(tài)渦流場 ［J］.南京理工大學(xué)學(xué)報(bào):自然科學(xué)版，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］彭紅梅，楊德全.雙分叉動(dòng)脈血流動(dòng)力學(xué)特性的邊界元分析 ［J］.醫(yī)用生物力學(xué)，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］黃玉美，張廣鵬，高峰.虛擬樣機(jī)整機(jī)結(jié)構(gòu)特性邊界元法仿真［M］.北京:機(jī)械工業(yè)出版社，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］邢譽(yù)峰，李敏.計(jì)算固體力學(xué)原理與方法 ［M］.北京:北京航空航天大學(xué)出版社，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］邢譽(yù)峰，李敏.工程振動(dòng)基礎(chǔ) ［M］.北京:北京航空航天大學(xué)出版社，2011

Xing Yufeng，Li Min.The foundation of engineering vibration［M］.Beijing:Beijing University of Aeronautics and Astronautics Press，2011(in Chinese)