高岱恒，郭宏偉，鄭 宏

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

1 研究背景

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

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

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

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

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

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

2 梁振動(dòng)分析的有限元方法

2.1 梁振動(dòng)的變分提法

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

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

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

2.2 空間離散

假設(shè)梁?jiǎn)卧浑x散成ne個(gè)單元和n個(gè)節(jié)點(diǎn)。對(duì)于每一個(gè)單元e，梁撓度的一階Hermite插值形式為

其中：

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

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

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

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

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

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

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

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

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

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

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

為此，先從NMM的角度出發(fā)，審視梁?jiǎn)栴}的有限元法。

由圖1可知，梁由3個(gè)梁?jiǎn)卧M成。這里，物理片1、物理片2、物理片3、物理片4可以用點(diǎn)1，2，3，4來(lái)表示，這4個(gè)點(diǎn)都是拓?fù)湫屈c(diǎn)。

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

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

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

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

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

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

慣性力虛功一般形式為

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

其中，

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

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

也可以寫(xiě)成

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

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

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

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

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

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

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

4 算 例

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

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

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

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

表1 一致質(zhì)量矩陣與本文集中質(zhì)量矩陣計(jì)算固有頻率的對(duì)比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個(gè)單元和1 000個(gè)單元時(shí)，用一致質(zhì)量矩陣和本文方法的速度（計(jì)算使用的電腦配置：Intel Core i5-6500＠3.20 GHz，4 GB內(nèi)存）。結(jié)果顯示，當(dāng)將結(jié)構(gòu)劃分為500個(gè)單元時(shí)，后者比前者快約21%；當(dāng)將結(jié)構(gòu)劃分為1 000個(gè)單元時(shí)，后者比前者快約2倍。顯然，計(jì)算速度隨著結(jié)構(gòu)劃分單元數(shù)的增加，采用本方法計(jì)算高階模態(tài)效應(yīng)的速度顯著高于采用一致質(zhì)量矩陣方法。

4． 2 Euler-Bernoulli梁的時(shí)域分析

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

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

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

圖6 時(shí)域計(jì)算的結(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ì)量矩陣計(jì)算簡(jiǎn)諧荷載作用下梁中心點(diǎn)位移的對(duì)比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é) 語(yǔ)

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

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