基于有限單元的動力響應分析法在MATLAB 中的應用研究

1 引言

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

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

其中，單元體的運動方程中包含有單元體的剛度矩陣、質(zhì)量矩陣和阻尼矩陣。在將單元體運動方程集合成結(jié)構整體運動方程的過程中需要將單元體剛度矩陣、質(zhì)量矩陣和阻尼矩陣集合成相應的整體矩陣。當結(jié)構需要得到更高的計算精度時，可以通過增加結(jié)構的單元劃分數(shù)量來實現(xiàn)，但是，勢必也會增加矩陣之間的運算量。

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

2 基于有限單元的動力響應分析法

2.1 簡介

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

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

2.2 Newmark 法

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

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

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

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

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

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

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

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

當β=1/4、γ=1/2 時，式（4）、（5）的解有如下形式：

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

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

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

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

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

3 Newmark 法在MATLAB 中的應用及算例

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

圖1 移動荷載作用下的簡支梁

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

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

圖2 Newmark 法MATLAB 程序流程圖

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

圖3 簡支梁跨中截面不同自由度的撓度時程曲線

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

4 討論與展望

基于有限單元的動力響應分析法的求解精度與

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

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