李鴻晶 王競雄

（中國南京 211816 南京工業(yè)大學工程力學研究所）

引言

譜元法（spectral element method，縮寫為SEM）是結(jié)合譜方法的高精度配置點和有限元法的分片近似物理場思想而建立起來的，可視其為一種高階的有限單元法.該方法首先由Patera （1984）針對流體動力學問題提出，目前已在地震波傳播（Komatitsch，Vilotte，1998；Komatitsch，Tromp，1999；丁志華等，2014；韓天成等，2020）、結(jié)構(gòu)動力分析（Kudelaet al，2007a，b；?aket al，2017；?ak，Krawczuk，2018）、海洋聲學（Cristini，Komatitsch，2012；Botteroet al，2016）等多個領(lǐng)域得到廣泛應用.在處理動力問題時，按照傳統(tǒng)有限單元法格式建立的質(zhì)量特性模型一般為非對角的一致質(zhì)量矩陣（consistent mass matrix，縮寫為CMM）形式，每個時間步均需對CMM 求逆運算和存儲，計算量巨大，尤其對于求解諸如沖擊、爆炸和彈性波傳播等高頻或大范圍解域的動力學問題而進行的大規(guī)模數(shù)值模擬，產(chǎn)生的巨大計算工作量即便采用高性能計算設(shè)備也難以承受.因而，在處理這類大規(guī)模動力計算問題時往往采用所謂“質(zhì)量集中”技術(shù)，即首先通過有限單元法導出CMM，然后再采用某種方式將CMM 等效轉(zhuǎn)變?yōu)閷切问降募匈|(zhì)量矩陣（lumped mass matrix，縮寫為LMM），以實現(xiàn)時空解耦的顯式算法構(gòu)建及其大規(guī)模動力計算.

質(zhì)量集中的首要原則是保持總質(zhì)量不變.早期的質(zhì)量集中方法簡單地將單元總質(zhì)量平均分配到各節(jié)點上，具有相當大的隨意性.Fried 和Malkus （1975）提出的節(jié)點積分法取洛巴托（Lobatto）積分點作為單元節(jié)點，并采用洛巴托積分計算單元質(zhì)量矩陣，積分之后可自動形成LMM.該方法具有完備的數(shù)學基礎(chǔ)（Duczek，Gravenkamp，2019a），但應用于Serendipity 單元時將導致質(zhì)量矩陣出現(xiàn)零元素或負元素（Zienkiewiczet al，2013）.Hinton 等（1976）提出的對角元素放大法是在保持總質(zhì)量不變的前提下，將CMM 中的主對角元素按比例放大，非對角元素全部置零.該方法能夠保證主對角元素均為正值，但缺乏數(shù)學基礎(chǔ)（Hughes，1987）.另一種較為常用的方法為行和集中法（Hughes，1987；王勖成，2003；Zienkiewiczet al，2013），即將CMM 中每一行元素都累加到主對角元素上，以實現(xiàn)質(zhì)量矩陣對角化.但該方法應用于某些類型單元（如6 節(jié)點2 階三角形單元和8 節(jié)點Serendipity 四邊形單元）時可能出現(xiàn)零元素或負元素，導致動力分析出現(xiàn)問題.Zheng 和Yang （2017），Yang 等（2017）基于數(shù)值流形的概念提出了一種在數(shù)學上嚴格且通用的質(zhì)量集中方法.Zhang 等（2019a，b）將該方法成功應用于6 節(jié)點三角形單元和10 節(jié)點四面體單元.Duczek 和Gravenkamp （2019b）將這一方法拓展到二維和三維高階Serendipity 單元，認為該方法不能顯著提高收斂性.

雖然應用節(jié)點積分法構(gòu)建傳統(tǒng)有限單元法質(zhì)量特性模型還存在一些問題，但對于一些具有張量積形式的有限單元格式（如時域SEM）還是能很好地構(gòu)建出LMM.時域SEM 主要有兩種形式，最初被提出時是以高斯-洛巴托-切比雪夫（Gauss-Lobatto-Chebyshev ，縮寫為GLC）積分點為基礎(chǔ)配置單元節(jié)點，通過切比雪夫（Chebyshev）正交多項式構(gòu)造單元位移模式，如此建立的譜元格式稱為切比雪夫譜元法.切比雪夫譜元法在構(gòu)建質(zhì)量矩陣時使用了解析積分（Patera，1984；Prioloet al，1994；Zhuet al，2011），導出的質(zhì)量特性矩陣為CMM.后來又發(fā)展出另外一種形式的時域SEM，即單元節(jié)點按照高斯-洛巴托-勒讓德（Guass-Lobatto-Legendre，縮寫為GLL）積分點進行配置，且將單元位移模式形函數(shù)取為全部單元節(jié)點上的拉格朗日插值基函數(shù)，該譜元格式被稱為勒讓德譜元法.勒讓德譜元法構(gòu)建的質(zhì)量特性模型自動為LMM，這是因為其運用了高斯-洛巴托（Gauss-Lobatto）型積分計算質(zhì)量矩陣中的各個元素，即所選取的積分節(jié)點與單元節(jié)點一致，因而利用了單元形函數(shù)的克羅內(nèi)克（Kronecker-δ）性質(zhì).一些學者嘗試將切比雪夫譜元法的質(zhì)量矩陣對角化，如Dauksher 和Emery （1997，1999，2000）利用行和集中法和對角元素放大法構(gòu)造集中質(zhì)量切比雪夫譜元模型求解標量波方程和彈性靜動力問題，結(jié)果表明行和集中法的誤差最小.?ak （2009）也采用了行和集中法建立切比雪夫譜板單元的LMM，求解板中的彈性波傳播問題.需要指出的是，即便采用行和集中法建立集中質(zhì)量切比雪夫譜單元，也必須事先算出CMM.而由于傳統(tǒng)切比雪夫譜元法采用解析積分計算質(zhì)量矩陣，所以導出CMM 的過程也要耗費大量計算資源.邢浩潔（2017）分析了SEM 中分別采用GLL 積分和GLC 積分時得到的質(zhì)量矩陣的區(qū)別.Duczek 和Gravenkamp（2019a）討論了在勒讓德譜元法中行和集中、對角元素放大和節(jié)點積分三種質(zhì)量集中方法的等價性，指出當計算域中質(zhì)量密度和單元幾何形狀恒定不變時，這三種方法是完全等價的.

本文將主要探討時域SEM 中的質(zhì)量特性模型及其構(gòu)建問題，深入分析時域譜單元質(zhì)量特性模型的數(shù)學機理，以期得到一種在切比雪夫譜元模型中直接導出LMM 的數(shù)學方法，避免質(zhì)量集中技術(shù)的不確定性，減小計算成本，并試圖從物理機制上解釋質(zhì)量特性模型的合理性.

1 時域譜元質(zhì)量特性模型

以一維等參單元為例，闡釋時域譜元模型中質(zhì)量矩陣的構(gòu)造過程.在標準區(qū)間ξ∈[?1，1]上建立參考單元，譜單元質(zhì)量矩陣可以統(tǒng)一地寫為

1.1 切比雪夫譜單元質(zhì)量模型

式中，ξi是GLC 點在標準區(qū)間上的坐標.切比雪夫譜單元可采用拉格朗日（Lagrange）插值函數(shù)構(gòu)建單元形函數(shù)，即

式中，系數(shù)ci＝2 （當i＝0 或p－1 時）或ci＝1 （當i＝1，···，p－2 時）.根據(jù)多項式插值的唯一性可知，式（4）與式（5）本質(zhì)相同，僅在多項式的計算過程中有所差異.式（5）形式的優(yōu)點在于可以方便地利用切比雪夫正交多項式的性質(zhì)精確求解單元特性矩陣.例如，k階和l階切比雪夫多項式的乘積在 ［ ?1，1 ］ 上定積分Akl的解析解為

顯見，按照式（7）導出的切比雪夫譜單元的質(zhì)量矩陣是非對角形式的CMM.

1.2 勒讓德譜單元質(zhì)量模型

由此導出的譜單元質(zhì)量矩陣可自動形成LMM.

表1 GLL 節(jié)點坐標和積分權(quán)系數(shù)Table 1 Abscissas of GLL points and quadrature weights

2 兩種質(zhì)量模型的數(shù)學機制

2.1 積分方法對質(zhì)量矩陣的影響

無論何種類型的譜單元，其質(zhì)量特性模型均按照式（1）或式（2）的規(guī)則構(gòu)建，需要對譜單元形函數(shù)進行積分計算.對于切比雪夫譜單元和勒讓德譜單元，兩者在導出質(zhì)量特性矩陣過程中使用了不同的積分策略.切比雪夫譜單元利用解析積分計算質(zhì)量矩陣，而勒讓德譜單元在導出質(zhì)量矩陣時采用了GLL 數(shù)值積分.積分方式的不同最終導致兩種譜單元的質(zhì)量矩陣呈現(xiàn)不同的形式，即切比雪夫譜單元導出CMM，而勒讓德譜單元導出LMM.

從積分精度分析可以看出，切比雪夫譜單元導出的質(zhì)量矩陣為精確積分結(jié)果，而勒讓德譜單元使用的GLL 數(shù)值積分只具有（2p－3）階代數(shù)精度.雖然高斯-洛巴托積分屬于高精度數(shù)值積分，但對于擁有p個單元節(jié)點的（p－1）階譜單元而言，被積函數(shù)（兩個形函數(shù)的乘積）的階次達到（2p－2）階，故勒讓德譜單元對質(zhì)量模型的積分計算并非精確積分結(jié)果，導出的質(zhì)量矩陣與式（1）的精確積分結(jié)果相比實際上存在一定誤差.事實上，若采用精確積分計算勒讓德譜單元質(zhì)量矩陣，導出的質(zhì)量矩陣也將是CMM 而非LMM.

以一維2 階勒讓德譜單元為例，設(shè)單元長度為ΔL＝4，質(zhì)量密度取為ρ＝1.在參考單元上，其GLL 節(jié)點坐標為（?1，0，1），單元形函數(shù)寫為

顯然，如此導出的勒讓德譜單元質(zhì)量矩陣為CMM.注意到高斯積分可以給出（2n－1）階代數(shù)精度，其中n為高斯積分點數(shù).具體到該2 階勒讓德譜單元，有n＝3，而單元質(zhì)量模型中被積函數(shù)的階次為4 階，故通過高斯積分得到的形如式（12）的質(zhì)量矩陣實際上是對式（1）的精確積分結(jié)果，不存在任何誤差.類似地，對于切比雪夫譜單元，如果采用足夠數(shù)量積分點的高斯積分（當2n－1≥2p－2 時）計算式（2），導出的切比雪夫譜單元質(zhì)量矩陣將與式（7）給出的結(jié)果完全相同，均為精確積分結(jié)果.也就是說，無論切比雪夫譜單元還是勒讓德譜單元，如果按照式（1）或式（2）的精確積分結(jié)果確立譜單元質(zhì)量矩陣，它們必然都是CMM 形式.而傳統(tǒng)勒讓德譜單元質(zhì)量矩陣為LMM 形式，其實是其采用了非精確的高斯-洛巴托積分計算譜單元質(zhì)量矩陣的結(jié)果.

另一方面，也可以采用勒讓德譜單元所使用的高斯-洛巴托型數(shù)值積分導出切比雪夫譜單元的質(zhì)量矩陣.注意到切比雪夫譜單元與勒讓德譜單元配置單元節(jié)點的方式有所不同，所以相應的數(shù)值積分應為高斯-洛巴托-切比雪夫（GLC）積分，而非勒讓德譜單元所使用的高斯-洛巴托-勒讓德（GLL）積分.我們推導了這種基于GLC 節(jié)點坐標的GLC 數(shù)值積分形式，其1—5 階GLC 點坐標和對應的積分權(quán)系數(shù)列于表2 中.

表2 GLC 節(jié)點坐標和高斯-洛巴托積分權(quán)系數(shù)Table 2 Abscissas of GLC points and Gauss-Lobatto quadrature weights

下面同樣以上述的一維2 階譜單元為例考察譜單元質(zhì)量矩陣情況，只是譜單元節(jié)點按照GLC 點進行配置，即該譜單元變?yōu)榍斜妊┓蜃V單元，單元參數(shù)與上述算例相同.按照表2提供的GLC 積分權(quán)系數(shù)計算導出切比雪夫譜單元的質(zhì)量矩陣，結(jié)果如下：

由此可見，當切比雪夫譜單元也采用高斯-洛巴托積分時，導出的質(zhì)量矩陣同樣為LMM.而GLC 積分得到的亦是近似積分結(jié)果，由其導出的切比雪夫譜單元質(zhì)量矩陣同精確積分結(jié)果同樣存在誤差.

2.2 質(zhì)量特性模型的數(shù)學解釋

前述可知，無論切比雪夫譜單元還是勒讓德譜單元，只要采用精確的高斯積分計算質(zhì)量特性模型，導出的質(zhì)量矩陣就必然是CMM 形式.而如果選擇非精確的高斯-洛巴托積分（積分節(jié)點與譜單元節(jié)點一致）計算式（1），則將直接導出LMM 形式的譜單元質(zhì)量矩陣.也就是說，采取不同的積分方法會導出不同類型的譜單元質(zhì)量矩陣.觀察高斯-洛巴托積分方式，其積分節(jié)點選擇為與譜單元節(jié)點一致，即對于切比雪夫譜單元和勒讓德譜單元分別選為GLC 點和GLL 點.由于譜單元形函數(shù)是通過拉格朗日插值函數(shù)構(gòu)建出來的，而拉格朗日函數(shù)具有克羅內(nèi)克性質(zhì)，即函數(shù)僅在定義點取值為1 而在其他節(jié)點取值皆等于0，從而形函數(shù)具有如下的離散正交性：由式（2）可知，譜單元質(zhì)量矩陣每個元素的積分式中都包含兩個拉格朗日多項式的乘積，所以當數(shù)值積分點與單元節(jié)點完全重合時，由式（14）描述的譜單元形函數(shù)離散正交性質(zhì)，導出的質(zhì)量矩陣中必然只有對角線元素（被積函數(shù)中的兩個形函數(shù)相同）非零而所有非對角元素皆為零.

以4 階譜單元為例標示了單元形函數(shù)與積分節(jié)點的關(guān)系，如圖1 所示.對于高斯-洛巴托積分，積分節(jié)點與譜單元節(jié)點完全一致，如果采用高斯積分，則高斯積分點與譜單元節(jié)點不重合.此種情況下各形函數(shù)在這些積分點的取值既不為0 也不等于1，導致任意兩個形函數(shù)的乘積在這些積分點處均不可能為0，從而使得質(zhì)量矩陣被構(gòu)建成為非對角的CMM.

圖1 節(jié)點積分方法示意圖Fig. 1 Schematic drawing of nodal quadrature method

3 譜單元質(zhì)量模型的特征

3.1 兩種譜單元的集中質(zhì)量分布特征

比較由節(jié)點積分法導出的切比雪夫譜單元和勒讓德譜單元的集中質(zhì)量分布特征.首先考慮一維標準單元，設(shè)單元質(zhì)量密度為1，單元尺度ΔL＝2，則單元總質(zhì)量為2.不同階次的切比雪夫譜單元和勒讓德譜單元的集中質(zhì)量模型分布，如圖2 所示.由圖可知，兩種譜單元均能保證單元總質(zhì)量不變，對于本例即各節(jié)點的質(zhì)量之和始終等于2.但兩種譜單元是按照不同的比例將總質(zhì)量分配到各單元節(jié)點的.對于3 節(jié)點譜單元（圖2a），切比雪夫譜單元的質(zhì)量分布情況與勒讓德譜單元完全相同，這是由于此情況下GLC 節(jié)點的坐標與GLL 節(jié)點完全相同，且二者對應的高斯-洛巴托積分權(quán)系數(shù)也相等，對于4 節(jié) 點和5 節(jié)點譜單元（圖2b 和圖2c），切比雪夫譜單元內(nèi)節(jié)點分配到的質(zhì)量比勒讓德譜單元相應內(nèi)節(jié)點要多一些，而端部節(jié)點分配到的質(zhì)量比勒讓德譜單元要少.

圖2 一維集中質(zhì)量切比雪夫譜單元（左）和勒讓德譜單元（右）質(zhì)量分布特征（a） 3 節(jié)點單元；（b） 4 節(jié)點單元；（c） 5 節(jié)點單元Fig. 2 Mass distribution characteristics of 1-D lumped mass Chebyshev （left） and Legendre （right） spectral elements（a） 3-node element；（b） 4-node element；（c） 5-node element

時域SEM 采用張量積形式構(gòu)造二維和三維譜單元，故二維和三維譜單元的形函數(shù)依然具有克羅內(nèi)克性質(zhì)，不難證明利用高斯-洛巴托積分法同樣能夠保證導出的質(zhì)量矩陣為LMM.下面以二維標準單元為例，比較切比雪夫譜單元與勒讓德譜單元的集中質(zhì)量分布特征.設(shè)單元大小為2×2，質(zhì)量密度為1，則單元總質(zhì)量為4.不同階次的切比雪夫譜單元和勒讓德譜單元的集中質(zhì)量分布情況如圖3 所示.可以看出，兩種二維譜單元同樣能夠保證單元總質(zhì)量不變.與一維3 節(jié)點單元類似，二維9 節(jié)點的切比雪夫譜單元與勒讓德譜單元的質(zhì)量分布完全相同（圖3a）.對于16 節(jié)點和25 節(jié)點譜單元，切比雪夫譜單元將更多的質(zhì)量分配到內(nèi)節(jié)點上，而勒讓德譜單元則將更多質(zhì)量分配給了邊界節(jié)點.由圖3c-f 可見，切比雪夫譜單元的角節(jié)點分配到的質(zhì)量大約只有勒讓德譜單元的一半，且邊節(jié)點的質(zhì)量也比勒讓德譜單元的小.

圖3 二維集中質(zhì)量切比雪夫（左）和勒讓德（右）譜單元質(zhì)量分布特征（a） 9 節(jié)點單元；（b） 16 節(jié)點單元；（c） 25 節(jié)點單元Fig. 3 Mass distribution characteristics of 2-D lumped mass Chebyshev （left） and Legendre （right） spectral elements（a） 9-node element；（b） 16-node element；（c） 25-node element

3.2 兩種質(zhì)量特性模型比較

下面以切比雪夫譜單元為例，比較CMM 與LMM 的影響.由式（2）可以看出，除了單元

圖4 一維5 節(jié)點切比雪夫譜單元一致質(zhì)量模型數(shù)學圖示Fig. 4 Mathematical representation for consistent mass model of 1-D 5-node Chebyshev element

對于由節(jié)點積分法導出的切比雪夫譜單元的LMM，主對角元素為

根據(jù)形函數(shù)的克羅內(nèi)克性質(zhì)，式（15）進一步簡化為

圖5 一維5 節(jié)點切比雪夫譜單元集中質(zhì)量模型數(shù)學圖示Fig. 5 Mathematical representation for lumped mass model of 1-D 5-node Chebyshev element

接下來通過特征值問題的求解，從數(shù)值計算的角度討論切比雪夫譜單元LMM 和CMM 的差異，同時給出集中質(zhì)量勒讓德譜單元的計算結(jié)果.考察一根簡支鐵木辛柯（Timoshenko）梁的自由振動問題，設(shè)該梁長度L＝3 m，截面寬度b＝0.02 m，高度h＝0.1 m，質(zhì)量密度ρ＝7 800 kg/m3，彈性模量E＝210 GPa.根據(jù)鐵木辛柯梁理論建立切比雪夫譜單元，分別運用高斯-洛巴托積分和解析積分構(gòu)造LMM 和CMM 兩種質(zhì)量模型，并計算其固有頻率.

采用不同階次的集中質(zhì)量切比雪夫譜單元、一致質(zhì)量切比雪夫譜單元和集中質(zhì)量勒讓德譜單元計算得到的簡支梁的前6 階固有頻率，列于表3.單元尺度固定為ΔL＝1 m，即簡支梁被離散為3 個譜單元.表中的解析解根據(jù)Craig 和Kurdila （2006）中公式計算.從表中數(shù)據(jù)可見，兩種切比雪夫譜單元計算結(jié)果的差距較小，均與解析解非常接近.隨著譜單元階次的提升，兩種質(zhì)量模型均能快速收斂于精確解.與集中質(zhì)量勒讓德譜單元相比，集中質(zhì)量切比雪夫譜單元的計算結(jié)果也十分接近解析解，計算精度大致相當.對于簡支梁的一階頻率，三種譜單元模型都只需劃分3 個4 階單元就已經(jīng)能夠給出相當準確的結(jié)果.

表3 簡支梁前6 階固有頻率計算結(jié)果Table 3 First six natural frequencies of simply supported beam

圖6 為采用不同數(shù)量的譜單元計算出的簡支梁前6 階頻率的均方根誤差，此時三種譜單元階次均保持4 階不變.由圖可見，集中質(zhì)量切比雪夫譜單元的誤差曲線與一致質(zhì)量切比雪夫譜單元的誤差曲線十分接近.隨著單元數(shù)量的增加，兩種質(zhì)量模型計算結(jié)果的誤差都能迅速降低，當單元數(shù)量達到6 個時，前6 階頻率的均方根誤差已降低至10?3以下.在某些情況下，如譜單元數(shù)等于4，5，7，8 時，LMM 給出的結(jié)果誤差略低于CMM.與集中質(zhì)量勒讓德譜單元相比，集中質(zhì)量切比雪夫譜單元除了單元數(shù)量取為3 時的誤差稍大以外，其余情況下的計算誤差均不超過集中質(zhì)量勒讓德譜單元的誤差.

圖6 不同譜單元數(shù)量下簡支梁的前6 階固有頻率計算誤差Fig. 6 Error of first six natural frequencies with different number of elements

從該算例分析可以看出，對于特征值問題，集中質(zhì)量切比雪夫譜單元，一致質(zhì)量切比雪夫譜單元和集中質(zhì)量勒讓德譜單元均能取得較好的計算效果，三種譜單元模型不存在明顯差距.

4 關(guān)于質(zhì)量模型物理機制的討論

有限單元法中廣泛使用的質(zhì)量模型有兩種，即呈對角形式的集中質(zhì)量模型和呈非對角形式的一致質(zhì)量模型.通常認為，一致質(zhì)量模型是準確的，而集中質(zhì)量模型是近似的.這是因為一致質(zhì)量模型是嚴格按照變分原理導出的，而且通過精確積分結(jié)果忠實地反映了有限單元質(zhì)量模型的真實特性，具有很強的數(shù)學基礎(chǔ).而集中質(zhì)量模型的形成則有較多人為因素的干預，例如基于工程物理概念“堆聚”形成LMM，或者通過行和集中等途徑人為地將CMM 轉(zhuǎn)換為LMM 等，似乎都多少缺乏數(shù)學支撐.一些學者（Chanet al，1993；Jensen，1996；Wu，2006）通過對比分析指出，采用集中質(zhì)量模型和一致質(zhì)量模型得到的動力響應數(shù)值結(jié)果并無顯著差別.但這些分析大多是從數(shù)值計算角度進行研究，事實上是默認一致質(zhì)量模型更具合理性，而集中質(zhì)量模型被認為只是同一致質(zhì)量模型相差不大但使用更方便的一種質(zhì)量模型.

質(zhì)量是物質(zhì)的最基本特性之一.根據(jù)經(jīng)典力學原理，質(zhì)量被認為是度量物體慣性大小的物理量，它是不變的，與物體運動速度無關(guān).在每個有限單元上，實際的連續(xù)介質(zhì)分布質(zhì)量被離散至單元的各個節(jié)點上，也就是說，有限單元質(zhì)量模型本質(zhì)上是試圖采用一種離散質(zhì)量體系近似地代替實際的連續(xù)質(zhì)量體系，用以描述連續(xù)介質(zhì)有限單元的慣性特性.因此，判斷有限單元質(zhì)量模型優(yōu)劣的標準應該是其描述實際連續(xù)介質(zhì)慣性的能力.以四節(jié)點平面應變單元為例（圖7），實際單元是由均勻連續(xù)介質(zhì)構(gòu)成.經(jīng)過離散化處理后，單元慣性特性由分配于4 個單元節(jié)點上的集中質(zhì)量體現(xiàn)，它們反映力與運動之間的關(guān)系，由牛頓第二定律刻畫.單元各節(jié)點質(zhì)量之間通過彈性連接產(chǎn)生力學聯(lián)系，這些彈性連接表征了單元的剛度特性，反映力與變形之間的關(guān)系，由虎克定律刻畫.每個單元節(jié)點都具有兩個自由度，分別用水平向位移u和豎向位移v進行描述.故該單元總計8 個自由度，對應的剛度特性模型和質(zhì)量特性模型均為8×8 階矩陣.根據(jù)牛頓第二定律，質(zhì)量系數(shù)meij表征的物理含義為使j自由度產(chǎn)生單位加速度時需要在i自由度上施加的力的大小.這里i，j＝1，2，···，8，它們既可以是某個單元節(jié)點的水平向位移自由度u，也可以是豎向位移自由度v.考慮某個固定時刻t，假定在t之前單元所有節(jié)點都處于靜止狀態(tài)（波動未傳播至該單元前即是這種狀態(tài)），在t時刻單元第i自由度（例如v1自由度）上突然施加了力fi作用，則在i自由度上必然同時產(chǎn)生加速度ai，對應質(zhì)點將沿i自由度方向發(fā)生運動（即質(zhì)點m1沿v1方向運動）.注意到t時刻該質(zhì)點尚未開始運動，需要經(jīng)過一小段時間后該質(zhì)點才能改變在i自由度方向上的靜止狀態(tài).即產(chǎn)生位移具有滯后性，t時刻時單元各節(jié)點之間并未發(fā)生相對位移，此時各質(zhì)點之間亦不存在相互作用的彈性力，因而在j自由度上對應質(zhì)點當然仍然保持靜止狀態(tài)（例如此時質(zhì)點m3在u3方向上處于靜止）.也就是說，若在i自由度上受到力的作用，其不會在j自由度上立即就產(chǎn)生加速度，由此知有＝0 （當i≠j時）.即除了主對角線元素外，質(zhì)量矩陣中非對角質(zhì)量元素應該全部為零，單元質(zhì)量矩陣應該為LMM 形式.

圖7 連續(xù)體離散為單元Fig. 7 Description from continuum to discretized element

圖8 描繪了單元內(nèi)的傳力路徑.施加在i自由度上的作用力fi瞬間引起質(zhì)點沿i自由度的加速度ai.經(jīng)過一小段時間Δt之后，質(zhì)點沿i自由度運動到一個新的位置di.在此過程中，由于質(zhì)點相對于其它質(zhì)點的位置發(fā)生了改變且產(chǎn)生了相對位移，質(zhì)點與其它質(zhì)點之間的彈性連接發(fā)生變形，導致該質(zhì)點與其余各質(zhì)點之間產(chǎn)生彈性力kij.彈性連接的變形與彈性力之間的關(guān)系由本構(gòu)方程決定.然后根據(jù)牛頓第二定律，施加在j自由度上的彈性力kij將會引起對應質(zhì)點沿j自由度方向產(chǎn)生加速度aj.根據(jù)上述分析可知，i自由度方向的加速度ai與j自由度方向的加速度aj并不是同時產(chǎn)生的.也就是說，施加在i自由度上的作用力并不能瞬時引起j自由度方向的加速度.在上述單元運動分析中，質(zhì)點之間的相互作用力通過彈性連接的變形進行傳遞，而彈性變形依賴于質(zhì)點的相對運動，并且需要經(jīng)過時間才能產(chǎn)生.該過程表明，在同一時刻不同質(zhì)點的加速度并不耦合.這與靜力情況下外力瞬時引起節(jié)點位移且節(jié)點力瞬間傳遞至所有節(jié)點是完全不同的.

圖8 單元內(nèi)力傳遞路徑Fig. 8 Transfer path of element internal force

此外，集中質(zhì)量模型合理性還可以從波速有限原理（廖振鵬，2002）方面予以論證.如果質(zhì)量矩陣的非對角元素不為零（即CMM 形式），則在單元任意節(jié)點施加外力將會瞬時引起整個體系中所有節(jié)點都產(chǎn)生加速度.表明在一個子域中產(chǎn)生的波動將會瞬間傳播至全域.顯然，這與波速有限的物理原理相違背，因為波在介質(zhì)中的傳播需要時間，各質(zhì)點按照距離波源的遠近依次開始振動，不可能同時開始振動.根據(jù)這一原理，不同節(jié)點的加速度之間不應該存在耦合.換言之，質(zhì)量矩陣的非對角項應該全部為零，即質(zhì)量矩陣應該是對角形式的LMM，這樣才能符合波動傳播速度有限這一物理原理.

按照上述討論，集中質(zhì)量模型比一致質(zhì)量模型在物理上更具合理性.但在實踐中，采用集中質(zhì)量模型得到的波動數(shù)值計算結(jié)果并非總是比一致質(zhì)量模型更好（Tonget al，1971；Zienkiewiczet al，2013）.前面已經(jīng)論述，無論是一致質(zhì)量矩陣CMM 還是集中質(zhì)量矩陣LMM，它們實際上都是有限單元質(zhì)量模型的積分結(jié)果，區(qū)別僅在于各自使用了不同的積分方法.由于一致質(zhì)量模型采用代數(shù)精度最高的高斯數(shù)值積分計算有限單元質(zhì)量矩陣，當選取足夠數(shù)量積分點數(shù)時能夠獲得精確積分結(jié)果，因而一般認為其優(yōu)于集中質(zhì)量模型.然而，關(guān)于一致質(zhì)量模型優(yōu)于集中質(zhì)量模型的認識需要一個前提，即如果按照式（1）構(gòu)建的有限單元質(zhì)量模型是“精確”的，則作為精確積分結(jié)果的CMM 必然亦是“精確”的.雖然高斯-洛巴托積分也是一種高精度數(shù)值積分，但通過其導出的LMM 與有限單元質(zhì)量模型精確積分結(jié)果相比存在一定誤差，從這個意義上看CMM 要好于LMM.但是，有限單元質(zhì)量模型本質(zhì)上是對原連續(xù)介質(zhì)質(zhì)量特性的一種近似化處理，其本身就不是“精確”的，在力學特性描述上不可避免地會帶來誤差.這種離散化和積分過程對單元質(zhì)量特性造成的影響，如圖9 所示.如果能夠精確計算質(zhì)量矩陣（即形成CMM），則由數(shù)值積分引起的誤差為零.由圖9 可見，精確積分之后的誤差界與離散化之后的誤差界一致，也就是說精確積分不會帶來額外的誤差.但是，如果質(zhì)量矩陣沒有被精確計算，如使用了高斯-洛巴托積分并導出LMM，則由數(shù)值積分引入的誤差可能為正也可能為負.因而，關(guān)于質(zhì)量特性的誤差有可能被擴大也可能被減小.圖9 中高斯-洛巴托積分之后的誤差界（虛線），它即可能在實線以內(nèi)也可能在實線以外.淺灰色區(qū)域表示描述質(zhì)量特性的誤差被高斯-洛巴托積分減小，說明在此情況下LMM 要優(yōu)于CMM；深灰色區(qū)域表示誤差被高斯-洛巴托積分放大，這表明此情況下LMM 對質(zhì)量特性的描述比CMM 要差.因此，不是十分精確的積分結(jié)果并不總是造成不利影響，它有可能會改善離散化造成的誤差，當然亦可能會增大離散化影響從而導致更大的誤差.這可能也是基于LMM 的計算結(jié)果在有些情況下好于CMM 的結(jié)果而在另外一些情況下表現(xiàn)更差的原因.不過，由于高斯-洛巴托積分精度只比精確的高斯-勒讓德積分略低一點，它們都屬于高精度數(shù)值積分方法，所以由于數(shù)值積分引起的誤差依然會被控制在有限范圍內(nèi).

圖9 離散化和積分過程引起的質(zhì)量特性誤差Fig. 9 Error in mass property caused by discretization and quadrature

事實上，動力分析結(jié)果不僅取決于對原連續(xù)體系慣性特性的模擬，同時還取決于對彈性特性估計的準確性，兩者具有耦合影響.一般認為，在有限元法中單元剛度矩陣高估了所模擬的連續(xù)體系的實際剛度（Bathe，2014）.一種直觀的理解是有限元離散化過程將原來的無限自由度體系變?yōu)榱擞邢拮杂啥润w系，這相當于人為地給連續(xù)體系施加了若干約束，從而減小了體系的變形.偏剛的估計誤差在某些情況下甚至可能造成計算結(jié)果不可接受，例如梁、板分析中遇到的剪切自鎖問題（王勖成，2003）.解決問題的一種途徑是采用減縮積分計算單元剛度矩陣.減縮積分的思想是通過不精確的數(shù)值積分計算結(jié)果適當彌補有限元離散化造成的剛度偏大的誤差（Bathe，2014），從而改善總體計算精度.但當單元階次較低時，減縮積分的實施有可能導致單元過柔，從而出現(xiàn)虛假的零能模式（王勖成，2003）.這種處理措施同本文中利用高斯-洛巴托積分形成譜單元集中質(zhì)量矩陣的過程相類似，即采用不精確的積分有可能彌補離散化引起的誤差，但也可能會增大離散化誤差.與有限元法相比，目前關(guān)于時域譜元模型所描述的體系剛度特性的誤差估計尚無明確結(jié)論，但通過數(shù)值積分誤差來適當?shù)卣{(diào)控譜單元離散化帶來的誤差，無疑是一條可以考慮的途徑.另外，由于靜力問題無需考慮力與運動的傳遞過程，故在有限元法中只是簡單地通過對位移場進行求導而獲得加速度場，而沒有完全反映出體系中力與運動傳遞的物理過程.這是動力問題與靜力問題的區(qū)別所在，也是本文認為集中質(zhì)量模型更加合理的基礎(chǔ).

5 結(jié)論

集中質(zhì)量矩陣求逆計算十分簡便，這對于大規(guī)模動力分析問題具有非常重要的意義.本文將時域譜元法的兩種質(zhì)量特性模型—集中質(zhì)量模型和一致質(zhì)量模型都統(tǒng)一在相同的數(shù)學框架下，給出了統(tǒng)一的構(gòu)建方法，并從數(shù)學機制和物理意義兩個方面闡述了譜單元質(zhì)量模型問題，結(jié)論如下：

1） 譜單元質(zhì)量特性矩陣的構(gòu)建歸結(jié)為對譜單元形函數(shù)積的積分運算，選擇不同的數(shù)值積分方法將導致不同形式的譜單元質(zhì)量矩陣.當采用高斯-勒讓德積分計算譜單元質(zhì)量模型的定積分式時將導出一致質(zhì)量矩陣，而選擇高斯-洛巴托積分時則會導出集中質(zhì)量矩陣.該結(jié)論對于切比雪夫譜單元和勒讓德譜單元都是適用的.

2） 采用不同積分方法導出不同形式質(zhì)量矩陣的根本原因在于質(zhì)量矩陣中的被積形函數(shù)在積分點處是否具有克羅內(nèi)克性質(zhì).高斯-洛巴托積分節(jié)點與譜單元節(jié)點一致，在積分點處譜單元形函數(shù)具有克羅內(nèi)克性質(zhì)，其導出的質(zhì)量矩陣必為對角陣；而高斯積分無法滿足積分點處譜單元形函數(shù)的克羅內(nèi)克函數(shù)條件，只能導出一致質(zhì)量矩陣.這是數(shù)值積分方法決定譜單元質(zhì)量矩陣形式的數(shù)學機制.

3） 波動在介質(zhì)中傳播需要時間，在某個自由度上的力作用只能改變其自身自由度上的運動狀態(tài)，而不會在其他自由度上立即產(chǎn)生加速度.從質(zhì)量的物理意義上看，不同自由度上的質(zhì)量系數(shù)不存在耦合，譜單元的集中質(zhì)量模型更具有物理合理性.

4） 當選擇的積分節(jié)點數(shù)與譜單元節(jié)點數(shù)相同時，高斯積分給出的是精確積分結(jié)果，而高斯-洛巴托積分結(jié)果存在誤差.由于譜單元質(zhì)量模型本身亦是對原連續(xù)介質(zhì)慣性特性的近似描述，數(shù)值積分誤差有可能放大亦可能降低譜單元離散化引起的誤差范圍.

5） 高斯-洛巴托積分是一種高精度數(shù)值積分，其代數(shù)精度僅略低于高斯積分，由其導出的譜單元集中質(zhì)量矩陣對譜單元離散化誤差的調(diào)整被限制在有限范圍內(nèi).與高斯積分導出的譜單元一致質(zhì)量矩陣相較，兩種質(zhì)量模型對于動力分析效果而言大致相當，都能取得很好的計算結(jié)果.