1 引言

采用有限單元對(duì)結(jié)構(gòu)進(jìn)行振動(dòng)分析的方法可分為動(dòng)力響應(yīng)分析與振動(dòng)特性分析2 種。動(dòng)力響應(yīng)分析主要求解的是結(jié)構(gòu)在外力作用下發(fā)生振動(dòng)時(shí)，結(jié)構(gòu)的位移、速度與加速度隨時(shí)間的變化關(guān)系；振動(dòng)特性分析主要求解的是結(jié)構(gòu)的固有屬性，其中包括固有頻率、模態(tài)振型等。

基于有限單元的動(dòng)力響應(yīng)分析法首先通過d’Alembert 原理建立單元體的運(yùn)動(dòng)方程，然后按照有限元的集合方法，將單元體的運(yùn)動(dòng)微分方程集合成結(jié)構(gòu)整體的運(yùn)動(dòng)方程，最后通過逐步積分的方法求解出結(jié)構(gòu)的位移、速度與加速度隨時(shí)間的變化關(guān)系，從而能夠研究出結(jié)構(gòu)在動(dòng)力荷載作用下的內(nèi)力變化以及破壞過程。

其中，單元體的運(yùn)動(dòng)方程中包含有單元體的剛度矩陣、質(zhì)量矩陣和阻尼矩陣。在將單元體運(yùn)動(dòng)方程集合成結(jié)構(gòu)整體運(yùn)動(dòng)方程的過程中需要將單元體剛度矩陣、質(zhì)量矩陣和阻尼矩陣集合成相應(yīng)的整體矩陣。當(dāng)結(jié)構(gòu)需要得到更高的計(jì)算精度時(shí)，可以通過增加結(jié)構(gòu)的單元?jiǎng)澐謹(jǐn)?shù)量來實(shí)現(xiàn)，但是，勢(shì)必也會(huì)增加矩陣之間的運(yùn)算量。

而MATLAB 作為一款優(yōu)異的數(shù)值計(jì)算軟件，其基礎(chǔ)數(shù)據(jù)單位為矩陣，能夠?qū)⒐こ虇栴}與數(shù)值計(jì)算之間實(shí)現(xiàn)更好的連接。結(jié)構(gòu)的剛度矩陣與質(zhì)量矩陣中往往含有大量的零元素，而MATLAB 中內(nèi)置的稀疏矩陣正好可以解決上述矩陣的運(yùn)算量問題，在節(jié)約計(jì)算機(jī)儲(chǔ)存空間的同時(shí)加快矩陣的運(yùn)算速度，對(duì)于處理此類結(jié)構(gòu)振動(dòng)分析問題具有優(yōu)異的性能[1,2]。

2 基于有限單元的動(dòng)力響應(yīng)分析法

2.1 簡(jiǎn)介

有限單元法通過將結(jié)構(gòu)離散成一組指定數(shù)量的單元體，并在每個(gè)離散的單元體的指定點(diǎn)處設(shè)置結(jié)點(diǎn)使得相鄰單元體的有關(guān)參數(shù)擁有一定的連續(xù)性，然后，用每個(gè)單元上假設(shè)的函數(shù)來近似表示整個(gè)結(jié)構(gòu)上待求的未知場(chǎng)函數(shù)[3]。

在動(dòng)力響應(yīng)分析所處理的結(jié)構(gòu)動(dòng)力學(xué)問題中，未知場(chǎng)函數(shù)一般為結(jié)構(gòu)自由振動(dòng)或者受迫振動(dòng)產(chǎn)生的位移、速度與加速度。結(jié)構(gòu)整體的位移、速度與加速度時(shí)程數(shù)據(jù)可以從每個(gè)單元體結(jié)點(diǎn)處的響應(yīng)數(shù)據(jù)近似求得。理論上結(jié)構(gòu)劃分單元體的數(shù)量越多，最后的結(jié)果便越精確。此時(shí)，結(jié)構(gòu)的動(dòng)力響應(yīng)分析便從一個(gè)連續(xù)體無限自由度問題轉(zhuǎn)換為多自由度問題，自由度的數(shù)量即為單元體數(shù)量。

2.2 Newmark 法

Newmark 法是一種常用的動(dòng)力響應(yīng)分析法。Newmark 法通過假設(shè)時(shí)間步內(nèi)的加速度分布，依據(jù)積分的方法來獲得時(shí)間步內(nèi)結(jié)構(gòu)的速度與位移表達(dá)式，從而得到時(shí)間步結(jié)束點(diǎn)處的速度與位移。

多自由度系統(tǒng)在外力荷載作用下的整體運(yùn)動(dòng)方程如下[4]：

式中，M、C 和K 分別為多自由度系統(tǒng)的整體質(zhì)量矩陣、阻尼矩陣和剛度矩陣；V（··t）、V（·t）、V（t）分別為多自由度系統(tǒng)的單元結(jié)點(diǎn)加速度、速度和位移；F（t）為多自由度系統(tǒng)的外力矩陣。

假設(shè)該系統(tǒng)具有初始位移以及初始速度如下：

將公式（2）代入公式（1）中可立刻求得多自由度系統(tǒng)初始加速度：

式中，M-1為整體質(zhì)量矩陣的逆矩陣。

Newmark 法假設(shè)的加速度分布模式有如下公式[5]：

當(dāng)式（5）中β=1/4、γ=1/2 時(shí)，Newmark 法假設(shè)的加速度分布模式為常加速度分布；當(dāng)式（5）中的β=1/6、γ=1、2 時(shí)，Newmark 法假設(shè)的加速度分布模式為線加速度分布。以下有關(guān)Newmark 內(nèi)容的討論均按照常加速度分布進(jìn)行。

當(dāng)β=1/4、γ=1/2 時(shí)，式（4）、（5）的解有如下形式：

在這里定義7 個(gè)常數(shù)：

代入系統(tǒng)的初始位移、初始速度和初始加速度后，（6）式中的等效荷載F~與等效剛度K~可以由下式（8）、（9）分別求得：

在依據(jù)式（6）、（7）、（8）求解出多自由度系統(tǒng)任意時(shí)間步的位移V（t+Δt）后，系統(tǒng)任意時(shí)刻的速度即可由式（10）求得：

上式中的多自由度系統(tǒng)任意時(shí)刻的加速度可由下式求得：

基于有限單元的Newmark 法的具體過程為：（1）連續(xù)體結(jié)構(gòu)離散化。劃分單元體的分割方案需要綜合考慮結(jié)構(gòu)受力情況、約束情況以及求解結(jié)果的精度要求，劃分單元體的形狀應(yīng)根據(jù)具體工程問題選擇為桿單元、梁?jiǎn)卧ò矫媪簡(jiǎn)卧约翱臻g梁?jiǎn)卧?、三角單元以及矩形單元等中的一種或多種。（2）根據(jù)具體選擇的單元形狀確定對(duì)應(yīng)的單元?jiǎng)偠染仃嚭蛦卧|(zhì)量矩陣。（3）將單元?jiǎng)偠染仃嚭蛦卧|(zhì)量矩陣拼接成結(jié)構(gòu)整體剛度矩陣和整體質(zhì)量矩陣。對(duì)于單元所處坐標(biāo)系不同的工程問題，還需要進(jìn)行坐標(biāo)轉(zhuǎn)換。（4）外力荷載的移置。針對(duì)外力的作用情況，依據(jù)靜力等效原理將外力向單元結(jié)點(diǎn)處移置。（5）獲取結(jié)構(gòu)的初始位移、速度和加速度向量以及選取時(shí)間步長(zhǎng)Δt。（6）根據(jù)式（6）、（8）、（9）、（10）、（11）求解出結(jié)構(gòu)任意步長(zhǎng)時(shí)刻的位移、速度和加速度向量。

3 Newmark 法在MATLAB 中的應(yīng)用及算例

以一座移動(dòng)荷載作用下的簡(jiǎn)支梁作為算例對(duì)象（見圖1），利用Newmark 法編寫MATLAB 程序來得到該結(jié)構(gòu)的動(dòng)力響應(yīng)數(shù)據(jù)。

圖1 移動(dòng)荷載作用下的簡(jiǎn)支梁

采用移動(dòng)力模型[6]來描述移動(dòng)荷載作用下的簡(jiǎn)支梁。假設(shè)梁的長(zhǎng)度為16m，截面積為7.2m2，彈性模量為3.25×1010N/m2，慣性矩為0.864m4，密度為摘要kg/m3。假設(shè)移動(dòng)荷載為6.38×105N，從0 時(shí)刻出發(fā)以60km/h 的速度勻速移動(dòng)，將該梁劃分為160 個(gè)梁?jiǎn)卧?，取Newmark 法時(shí)間步長(zhǎng)為0.摘要s。

使用MATLAB 作為編程工具求解上述問題的流程如圖2所示[7]。

圖2 Newmark 法MATLAB 程序流程圖

依據(jù)上述流程圖過程設(shè)計(jì)出的MATLAB 程序得到了簡(jiǎn)支梁在移動(dòng)荷載作用下的動(dòng)力響應(yīng)數(shù)據(jù)。如圖3a 所示，約0.5s后，梁的跨中截面處出現(xiàn)最大撓度變形，此時(shí)荷載正好移動(dòng)到跨中截面處。

圖3 簡(jiǎn)支梁跨中截面不同自由度的撓度時(shí)程曲線

設(shè)計(jì)MATLAB 程序時(shí)，對(duì)每一個(gè)單元結(jié)點(diǎn)處劃分成3 個(gè)自由度方向，分別是X方向（梁的水平方向）、Y方向（梁的撓度方向）以及梁的彎矩方向。圖3b 是梁跨中截面處單元1 自由度（X方向）上的位移時(shí)程曲線，由于梁沒有發(fā)生水平方向上的位移，故所有時(shí)刻的位移均為0。由此可從一定程度上驗(yàn)證該方法的有效性。

4 討論與展望

基于有限單元的動(dòng)力響應(yīng)分析法的求解精度與

結(jié)構(gòu)劃分的單元數(shù)量密切相關(guān)，單元數(shù)量增加的同時(shí)也會(huì)增加未知場(chǎng)函數(shù)矩陣的計(jì)算量，此時(shí)選擇MATLAB 作為處理此類問題的工具將極大地縮短計(jì)算時(shí)長(zhǎng)，并可以提升求解精度。在考慮移動(dòng)荷載作用下結(jié)構(gòu)的時(shí)程分析時(shí)，移動(dòng)荷載的等效也將成為影響時(shí)程數(shù)據(jù)求解精度的影響因素，本文中MATLAB 程序設(shè)計(jì)時(shí)考慮在單元數(shù)量足夠的情況下將移動(dòng)荷載視作為相應(yīng)時(shí)刻各對(duì)應(yīng)單元的中心荷載，然后移置成結(jié)點(diǎn)等效荷載，從而形成整個(gè)時(shí)間域上結(jié)點(diǎn)力矩陣，以此對(duì)移動(dòng)荷載的移動(dòng)過程進(jìn)行模擬。

在實(shí)際工程應(yīng)用中，該程序的精度要求仍存在疑問，在未來工作中，可考慮依據(jù)相關(guān)精度要求對(duì)單元?jiǎng)澐謹(jǐn)?shù)量進(jìn)行調(diào)整，或者在移動(dòng)荷載移置成為單元結(jié)點(diǎn)等效荷載的過程中添加除中心荷載外的其余荷載位置來逐步提高該方法的精度。