時域譜元法的質(zhì)量特性模型及其構(gòu)建方法＊

李鴻晶 王競雄

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 一維集中質(zhì)量切比雪夫譜單元（左）和勒讓德譜單元（右）質(zhì)量分布特征（a） 3 節(jié)點(diǎn)單元；（b） 4 節(jié)點(diǎn)單元；（c） 5 節(jié)點(diǎn)單元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ì)，不難證明利用高斯-洛巴托積分法同樣能夠保證導(dǎo)出的質(zhì)量矩陣為LMM.下面以二維標(biāo)準(zhǔn)單元為例，比較切比雪夫譜單元與勒讓德譜單元的集中質(zhì)量分布特征.設(shè)單元大小為2×2，質(zhì)量密度為1，則單元總質(zhì)量為4.不同階次的切比雪夫譜單元和勒讓德譜單元的集中質(zhì)量分布情況如圖3 所示.可以看出，兩種二維譜單元同樣能夠保證單元總質(zhì)量不變.與一維3 節(jié)點(diǎn)單元類似，二維9 節(jié)點(diǎn)的切比雪夫譜單元與勒讓德譜單元的質(zhì)量分布完全相同（圖3a）.對于16 節(jié)點(diǎn)和25 節(jié)點(diǎn)譜單元，切比雪夫譜單元將更多的質(zhì)量分配到內(nèi)節(jié)點(diǎn)上，而勒讓德譜單元則將更多質(zhì)量分配給了邊界節(jié)點(diǎn).由圖3c-f 可見，切比雪夫譜單元的角節(jié)點(diǎn)分配到的質(zhì)量大約只有勒讓德譜單元的一半，且邊節(jié)點(diǎn)的質(zhì)量也比勒讓德譜單元的小.

圖3 二維集中質(zhì)量切比雪夫（左）和勒讓德（右）譜單元質(zhì)量分布特征（a） 9 節(jié)點(diǎn)單元；（b） 16 節(jié)點(diǎn)單元；（c） 25 節(jié)點(diǎn)單元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é)點(diǎn)切比雪夫譜單元一致質(zhì)量模型數(shù)學(xué)圖示Fig. 4 Mathematical representation for consistent mass model of 1-D 5-node Chebyshev element

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

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

圖5 一維5 節(jié)點(diǎn)切比雪夫譜單元集中質(zhì)量模型數(shù)學(xué)圖示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ù)鐵木辛柯梁理論建立切比雪夫譜單元，分別運(yùn)用高斯-洛巴托積分和解析積分構(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é)果也十分接近解析解，計算精度大致相當(dāng).對于簡支梁的一階頻率，三種譜單元模型都只需劃分3 個4 階單元就已經(jīng)能夠給出相當(dāng)準(zhǔn)確的結(jié)果.

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

圖6 為采用不同數(shù)量的譜單元計算出的簡支梁前6 階頻率的均方根誤差，此時三種譜單元階次均保持4 階不變.由圖可見，集中質(zhì)量切比雪夫譜單元的誤差曲線與一致質(zhì)量切比雪夫譜單元的誤差曲線十分接近.隨著單元數(shù)量的增加，兩種質(zhì)量模型計算結(jié)果的誤差都能迅速降低，當(dāng)單元數(shù)量達(dá)到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ì)量模型物理機(jī)制的討論

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

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