高岱恒，郭宏偉，鄭 宏

（1.北京工業(yè)大學(xué)建筑工程學(xué)院，北京 100124；2.中國科學(xué)院武漢巖土力學(xué)研究所，武漢 430071）

1 研究背景

梁作為工程結(jié)構(gòu)中基本的構(gòu)件，被廣泛地應(yīng)用于工程建設(shè)當(dāng)中。對于梁結(jié)構(gòu)的靜力、動力學(xué)研究始終是土木工程領(lǐng)域的重點(diǎn)問題。

目前，對慣性特性的合適表述可以顯著地降低計算開銷并得到精度較高的結(jié)果。有許多學(xué)者對質(zhì)量矩陣生成方式進(jìn)行了卓有成效的研究［1-4］。有限元中，常用的2類質(zhì)量矩陣為一致質(zhì)量矩陣（CMM）和集中質(zhì)量矩陣（DLMM）。

Archer［5-6］在使用一致質(zhì)量矩陣進(jìn)行動力學(xué)分析方面取得了非常好的成果。但使用一致質(zhì)量矩陣的缺點(diǎn)很突出：①無論是用顯式積分算法求結(jié)構(gòu)瞬態(tài)反應(yīng)還是對自由振動的廣義特征值求解，采用一致質(zhì)量矩陣的效率很低；②在接觸問題和波的傳播動力學(xué)中，采用一致質(zhì)量矩陣時，求出的解會出現(xiàn)震蕩現(xiàn)象。

基于此，在動力分析中，采用集中質(zhì)量矩陣會是一個更好的選擇。

然而，對于二維問題，以板單元為例，用廣為采用的row-sum［7］方法形成的集中質(zhì)量矩陣，其轉(zhuǎn)角自由度對應(yīng)的質(zhì)量為0，這樣會導(dǎo)致集中質(zhì)量矩陣正定性喪失，從而無法直接用于Newmark法和子空間迭代法，這樣不僅使實現(xiàn)變得困難，而且犧牲了很多數(shù)值特性。而HRZ方法［8］——對角縮放方法對任意單元（如三角形單元、八節(jié)點(diǎn)等參元等）都可以生成正定的集中質(zhì)量矩陣，但是卻沒有嚴(yán)格的數(shù)學(xué)基礎(chǔ)。盡管HRZ方法在結(jié)構(gòu)力學(xué)、固體力學(xué)以及熱傳導(dǎo)等方面得到了成功的應(yīng)用，但是由于其理論基礎(chǔ)不嚴(yán)格，始終值得警惕。根據(jù)現(xiàn)有的實踐表明，HRZ方法在流體計算中效果不佳。

本文基于數(shù)值流形方法（numerical manifold method，NMM），提出了一種基于嚴(yán)格數(shù)學(xué)理論的對角化集中質(zhì)量矩陣的生成方法，并將其應(yīng)用于梁動力分析中。

2 梁振動分析的有限元方法

2.1 梁振動的變分提法

本文采用均質(zhì)等截面歐拉-伯努利（Euler-Bernoulli）梁來對結(jié)構(gòu)進(jìn)行分析。為了敘述上的方便，以簡支梁為例，其受迫振動的變分提法為

式中：（δw，w··）為慣性力所做的虛功；a（δw，w）為廣義力做的虛功；l（δw）為分布荷載所做的虛功。展開形式為：

式中：w（x，t）為梁的橫向位移；ρ為梁的密度；A為梁的截面面積；L為梁的長度；E為材料的彈性模量；I為梁的慣性矩；q（x，t）為分布荷載；x，t分別為沿梁軸方向的坐標(biāo)和時間；δw中的δ表示變分符號。

2.2 空間離散

假設(shè)梁單元被離散成ne個單元和n個節(jié)點(diǎn)。對于每一個單元e，梁撓度的一階Hermite插值形式為

其中：

將式（6）代入式（1）中，可以得到半離散形式的微分方程組，即

式中：M為整體質(zhì)量矩陣；K為整體剛度矩陣；f為整體載荷向量，都是由相應(yīng)的單元矩陣組集而成；d為所有結(jié)點(diǎn)自由度組成的向量。

對于式（13）的求解，有2種普遍的方法：振型疊加法和對時間的直接積分法。按照計算每一時刻的動力反應(yīng)是否需要求解耦聯(lián)方程組，可以將直接積分法分為隱式算法和顯式算法2類。當(dāng)采用顯式算法時，對自由度數(shù)目很大的情況，采用對角化集中質(zhì)量矩陣，對節(jié)省計算資源開銷和內(nèi)存占用來說意義重大。此外，在一些場合還可消除解的偽震蕩現(xiàn)象。

同樣，當(dāng)使用振型疊加法求解問題時，如果M為對角化矩陣時，可以將問題形式從廣義特征值降為常義特征值，這同樣降低了計算開銷（因為質(zhì)量矩陣是稀疏的）。相比使用一致質(zhì)量矩陣，采用集中質(zhì)量矩陣會更容易找出結(jié)構(gòu)的更多高階模態(tài)來進(jìn)一步完善分析。

3 基于數(shù)值流形法（NMM）的集中質(zhì)量矩陣生成

3.1 數(shù)值流形法（NMM）的提出

NMM是石根華博士［10］于1991年首次提出，它基于2套覆蓋（數(shù)學(xué)覆蓋和物理覆蓋）和接觸環(huán)路來建立。對于求解巖體的不連續(xù)變形和節(jié)理裂隙擴(kuò)張的模擬等問題非常有效。

數(shù)學(xué)覆蓋（Mathematical Cover）可以視為一系列數(shù)學(xué)片的集合，與有限元不同：數(shù)學(xué)覆蓋通常大于實際結(jié)構(gòu)的物理區(qū)域，物理覆蓋（Physical Cover）是物理區(qū)域內(nèi)所有物理片的集合。物理片是由物理區(qū)域上如邊界、裂縫等不連續(xù)面切分?jǐn)?shù)學(xué)片后得到的，其中每個物理片都定義一個權(quán)函數(shù)［11］和相應(yīng)的局部逼近。

數(shù)值流形法是一種可以統(tǒng)一處理連續(xù)問題和非連續(xù)問題的計算方法，目前在巖土工程中大量使用，將來有望在更廣闊的領(lǐng)域得到應(yīng)用。

3.2 集中質(zhì)量矩陣生成

本方法的理論［12-13］是一種關(guān)于高階等參元的對角質(zhì)量矩陣的生成方法，但是，用于四階問題所采用的Hermite插值并不在單位分解框架內(nèi)，所以在應(yīng)用文獻(xiàn)［12］的方法之前，需要從Hermite插值中取出片上的權(quán)函數(shù)及其局部逼近。

為此，先從NMM的角度出發(fā)，審視梁問題的有限元法。

由圖1可知，梁由3個梁單元組成。這里，物理片1、物理片2、物理片3、物理片4可以用點(diǎn)1，2，3，4來表示，這4個點(diǎn)都是拓?fù)湫屈c(diǎn)。

圖1 結(jié)構(gòu)物理覆蓋Fig．1 Physical cover of structure

在本文這個特例中，流形單元即為普通單元，是該單元的2個節(jié)點(diǎn)所對應(yīng)的物理片的公共部分。即：流形單元1為物理片1和物理片2的共同部分，流形單元2和3同樣為對應(yīng)物理片共同部分。流形單元1和2的形成如圖2所示。

圖2 流形單元的形成Fig．2 Formation of manifold elements

本文3.1節(jié)中提到了每個物理片都有對應(yīng)的權(quán)函數(shù)，取物理片上2個單元的Hermite插值的零階插值函數(shù)來組成權(quán)函數(shù)是很自然的選擇。圖3為各個物理片對應(yīng)的權(quán)函數(shù)的形狀。

圖3 各個物理片對應(yīng)的權(quán)函數(shù)形狀Fig．3 Shape of weights function on each patch

集中質(zhì)量矩陣的組裝方式有2種：一是按物理片的方式組裝；二是按流形單元的方法組裝。本文采用按物理片的方法說明組裝集中質(zhì)量矩陣的過程。

慣性力虛功一般形式為

式中：ωi為patch-i所對應(yīng)的區(qū)間；pi為其上的權(quán)函數(shù)，并滿足

其中，

式中：δwi和w··i分別為ωi的虛位移和虛加速度，即ωi上的局部逼近。下面，將從Hermite插值中產(chǎn)生權(quán)函數(shù){ pi}和局部近似 {wi}。

設(shè)流形單元e由物理片i和j覆蓋，其中，來自物理片i的局部近似對撓度w的貢獻(xiàn)可表示為

也可以寫成

相應(yīng)的局部逼近為

式中 k，e[ ]表示單元e的第k個節(jié)點(diǎn)所對應(yīng)的整體碼。相應(yīng)地有：

將式（23）、式（24）代入式（17），得到

將式（26）代入式（16），得到近似慣性力虛功為

對梁而言，由于有限元網(wǎng)格上的每個結(jié)點(diǎn)同時有平動自由度和轉(zhuǎn)動自由度，所以生成的ML為塊對角化的質(zhì)量矩陣，即

現(xiàn)在，可以直接用塊對角化的質(zhì)量矩陣來求解微分方程組式（13）。然而，為了更好地利用ML這個對稱正定的塊對角化結(jié)構(gòu)，進(jìn)行了如下擴(kuò)展，即將式（30）變換成對角化矩陣。

顯然，求解新形式的微分方程組（32）相比求解原來的方程組（13）效率要更高。

4 算 例

4.1 歐拉-伯努利（Euler-Bernoulli）梁模態(tài)分析

兩端簡支的Euler-Bernoulli梁，其彈性模量為E＝210 GPa，長度為 1 m，密度為 ρ＝7 860 kg／m3（單位改為正體），截面高和寬都為0.2 m（圖4）。將結(jié)構(gòu)劃分為40個單元。

圖4 等截面梁示意圖Fig．4 Rectangular cross-section beam

表1給出了等截面Euler-Bernoulli梁有限元法分別采用一致質(zhì)量矩陣和集中質(zhì)量矩陣的前5階固有頻率。在實際測試中發(fā)現(xiàn)，隨著對單元劃分得越細(xì)，用本文方法的集中質(zhì)量矩陣求解速度顯著高于用一致質(zhì)量矩陣。圖5為用本方法計算的前5階模態(tài)圖。

表1 一致質(zhì)量矩陣與本文集中質(zhì)量矩陣計算固有頻率的對比Table 1 Comparison of natural frequencies for simply supported rectangular cross-section beam between consistent mass matrix and diagonally-lumped mass matrix

從模態(tài)分析結(jié)果可以看出，使用本文集中質(zhì)量方法得出的梁結(jié)構(gòu)的前5階模態(tài)和解析解的誤差變化在可接受范圍內(nèi)。另外，比較分別采用500個單元和1 000個單元時，用一致質(zhì)量矩陣和本文方法的速度（計算使用的電腦配置：Intel Core i5-6500＠3.20 GHz，4 GB內(nèi)存）。結(jié)果顯示，當(dāng)將結(jié)構(gòu)劃分為500個單元時，后者比前者快約21%；當(dāng)將結(jié)構(gòu)劃分為1 000個單元時，后者比前者快約2倍。顯然，計算速度隨著結(jié)構(gòu)劃分單元數(shù)的增加，采用本方法計算高階模態(tài)效應(yīng)的速度顯著高于采用一致質(zhì)量矩陣方法。

4． 2 Euler-Bernoulli梁的時域分析

兩端簡支的Euler-Bernoulli梁，參數(shù)同本文4.1節(jié)，不考慮阻尼的影響。分析在結(jié)構(gòu)中間受到一個簡諧荷載F＝F0sin（4πt）作用下的結(jié)構(gòu)中心位置的位移變化情況。

其中，0≤t≤td，計算總持續(xù)時間td為1 s，時間步取0.001 s。這里，分別使用時域逐步積分法中的中心差分法和Newmark-β法來求結(jié)構(gòu)中點(diǎn)的位移變化。

圖6（a）是采用一致質(zhì)量矩陣的結(jié)果，圖6（b）是采用本文提出的對角化集中質(zhì)量矩陣的結(jié)果。由圖6、表2可以看出，使用本文集中質(zhì)量矩陣方法得出的梁結(jié)構(gòu)的時域反應(yīng)和采用一致質(zhì)量矩陣的時域反應(yīng)相差很小。

圖6 時域計算的結(jié)構(gòu)中點(diǎn)的位移變化情況Fig．6 Variations of displacement of center block of beam with step-by-step methods using consistent mass matrix and diagonally-lumped mass matrix

表2 一致質(zhì)量矩陣與本文集中質(zhì)量矩陣計算簡諧荷載作用下梁中心點(diǎn)位移的對比Table 2 Comparison of center displacement under harmonic load for simply supported rectangular cross-section beam between consistent mass matrix and diagonally-lumped mass matrix

5 結(jié) 語

從算例可以看出，在流形上積分，構(gòu)造的數(shù)學(xué)嚴(yán)格的集中質(zhì)量矩陣的方法不僅能夠有效地提高運(yùn)行效率，對二維問題中的8節(jié)點(diǎn)等參元，用本方法克服了一般的集中質(zhì)量矩陣不滿足角點(diǎn)正定性要求。更可貴的是，本方法有著嚴(yán)格的數(shù)學(xué)基礎(chǔ)和力學(xué)基礎(chǔ)?？梢灶A(yù)見的是，本方法對于二維乃至三維結(jié)構(gòu)（以板殼結(jié)構(gòu)為典型代表）的動力計算有著天然的優(yōu)越性。

因此，無論是直接積分法中顯式方法，還是隱式方法，無論是時域分析，還是頻域分析，都可用本文所建議的集中質(zhì)量矩陣代替一致質(zhì)量矩陣。