The mass property model and its implementation in the time-domain spectral element method
-
摘要: 研究了构建时域谱单元质量特性模型的数学机制，针对时域切比雪夫谱单元和勒让德谱单元建立了一种直接导出谱单元一致质量矩阵和集中质量矩阵的统一数学方法，对比分析两种谱单元质量特性模型的特征，并从物理角度探讨了谱单元质量特性模型的合理性。研究表明，数值积分点与谱单元节点选取是否一致是决定时域谱单元形成一致质量模型或集中质量模型的根本原因，当采用高斯-勒让德积分计算谱单元质量模型时将导出一致质量矩阵，而采用高斯-洛巴托积分则导出集中质量矩阵。而集中质量模型更具有物理合理性，两种谱单元质量特性模型优劣相当，均可取得很好的动力问题分析结果。Abstract: The mathematical mechanism of constructing mass property model for the time-domain spectral elements is studied in this paper. A unified mathematical method for directly deriving consistent and lumped mass matrix is established for the time-domain Chebyshev and Legendre spectral elements. The characteristics of two mass property models of the spectral elements are analyzed through comparison. Meanwhile, the rationality of mass property model of spectral element is discussed from physical perspective. This study reveals that the formation of consistent or lumped mass matrix in time-domain spectral elements depends on whether the quadrature points are coincident with the element nodes or not. To be specific, the Gauss-Legendre quadrature results in consistent mass matrix for spectral elements, and the Gauss-Lobatto quadrature leads to lumped mass matrix. The lumped mass matrix is more reasonable in physics. The two mass property models of spectral elements have comparable performance and they can both achieve good results for dynamic problems.
-
引言
谱元法（spectral element method，缩写为SEM）是结合谱方法的高精度配置点和有限元法的分片近似物理场思想而建立起来的，可视其为一种高阶的有限单元法。该方法首先由Patera （1984）针对流体动力学问题提出，目前已在地震波传播（Komatitsch，Vilotte，1998；Komatitsch，Tromp，1999；丁志华等，2014；韩天成等，2020）、结构动力分析（Kudela et al，2007a，b；Żak et al，2017；Żak，Krawczuk，2018）、海洋声学（Cristini，Komatitsch，2012；Bottero et al，2016）等多个领域得到广泛应用。在处理动力问题时，按照传统有限单元法格式建立的质量特性模型一般为非对角的一致质量矩阵（consistent mass matrix，缩写为CMM）形式，每个时间步均需对CMM求逆运算和存储，计算量巨大，尤其对于求解诸如冲击、爆炸和弹性波传播等高频或大范围解域的动力学问题而进行的大规模数值模拟，产生的巨大计算工作量即便采用高性能计算设备也难以承受。因而，在处理这类大规模动力计算问题时往往采用所谓“质量集中”技术，即首先通过有限单元法导出CMM，然后再采用某种方式将CMM等效转变为对角形式的集中质量矩阵（lumped mass matrix，缩写为LMM），以实现时空解耦的显式算法构建及其大规模动力计算。
质量集中的首要原则是保持总质量不变。早期的质量集中方法简单地将单元总质量平均分配到各节点上，具有相当大的随意性。Fried和Malkus （1975）提出的节点积分法取洛巴托（Lobatto）积分点作为单元节点，并采用洛巴托积分计算单元质量矩阵，积分之后可自动形成LMM。该方法具有完备的数学基础（Duczek，Gravenkamp，2019a），但应用于Serendipity单元时将导致质量矩阵出现零元素或负元素（Zienkiewicz et al，2013）。Hinton等（1976）提出的对角元素放大法是在保持总质量不变的前提下，将CMM中的主对角元素按比例放大，非对角元素全部置零。该方法能够保证主对角元素均为正值，但缺乏数学基础（Hughes，1987）。另一种较为常用的方法为行和集中法（Hughes，1987；王勖成，2003；Zienkiewicz et al，2013），即将CMM中每一行元素都累加到主对角元素上，以实现质量矩阵对角化。但该方法应用于某些类型单元（如6节点2阶三角形单元和8节点Serendipity四边形单元）时可能出现零元素或负元素，导致动力分析出现问题。Zheng和Yang （2017），Yang等（2017）基于数值流形的概念提出了一种在数学上严格且通用的质量集中方法。Zhang等（2019a，b）将该方法成功应用于6节点三角形单元和10节点四面体单元。Duczek和Gravenkamp （2019b）将这一方法拓展到二维和三维高阶Serendipity单元，认为该方法不能显著提高收敛性。
虽然应用节点积分法构建传统有限单元法质量特性模型还存在一些问题，但对于一些具有张量积形式的有限单元格式（如时域SEM）还是能很好地构建出LMM。时域SEM主要有两种形式，最初被提出时是以高斯-洛巴托-切比雪夫（Gauss-Lobatto-Chebyshev ，缩写为GLC）积分点为基础配置单元节点，通过切比雪夫（Chebyshev）正交多项式构造单元位移模式，如此建立的谱元格式称为切比雪夫谱元法.切比雪夫谱元法在构建质量矩阵时使用了解析积分（Patera，1984；Priolo et al，1994；Zhu et al，2011），导出的质量特性矩阵为CMM。后来又发展出另外一种形式的时域SEM，即单元节点按照高斯-洛巴托-勒让德（Guass-Lobatto-Legendre，缩写为GLL）积分点进行配置，且将单元位移模式形函数取为全部单元节点上的拉格朗日插值基函数，该谱元格式被称为勒让德谱元法。勒让德谱元法构建的质量特性模型自动为LMM，这是因为其运用了高斯-洛巴托（Gauss-Lobatto）型积分计算质量矩阵中的各个元素，即所选取的积分节点与单元节点一致，因而利用了单元形函数的克罗内克（Kronecker-δ）性质。一些学者尝试将切比雪夫谱元法的质量矩阵对角化，如Dauksher和Emery （1997，1999，2000）利用行和集中法和对角元素放大法构造集中质量切比雪夫谱元模型求解标量波方程和弹性静动力问题，结果表明行和集中法的误差最小。Żak （2009）也采用了行和集中法建立切比雪夫谱板单元的LMM，求解板中的弹性波传播问题。需要指出的是，即便采用行和集中法建立集中质量切比雪夫谱单元，也必须事先算出CMM。而由于传统切比雪夫谱元法采用解析积分计算质量矩阵，所以导出CMM的过程也要耗费大量计算资源。邢浩洁（2017）分析了SEM中分别采用GLL积分和GLC积分时得到的质量矩阵的区别。Duczek和Gravenkamp （2019a）讨论了在勒让德谱元法中行和集中、对角元素放大和节点积分三种质量集中方法的等价性，指出当计算域中质量密度和单元几何形状恒定不变时，这三种方法是完全等价的。
本文将主要探讨时域SEM中的质量特性模型及其构建问题，深入分析时域谱单元质量特性模型的数学机理，以期得到一种在切比雪夫谱元模型中直接导出LMM的数学方法，避免质量集中技术的不确定性，减小计算成本，并试图从物理机制上解释质量特性模型的合理性。
1. 时域谱元质量特性模型
以一维等参单元为例，阐释时域谱元模型中质量矩阵的构造过程。在标准区间ξ∈[−1，1]上建立参考单元，谱单元质量矩阵可以统一地写为
$$ {{\boldsymbol{M}}^{\rm{e}}} {\text{＝}} \int_{ - 1}^1 {\rho {{\boldsymbol{N}}^{\rm{T}}}{\boldsymbol{N}}\left| {\boldsymbol{J}} \right|{\rm{d}}\xi } {\text{，}} $$ (1) 式中：ρ为单元质量密度；
${\boldsymbol{N}} {\text{＝}} （ { {{N_1}（ \xi ）{\text{，}}{N_2}（ \xi ）{\text{，}} \cdots {\text{，}}{N_p}（\xi ）} } ）$ 为单元形函数向量，其中p为谱单元节点数，故对应的单元阶次为p－1阶；$\left| {\boldsymbol{J}} \right|$ 为雅可比行列式，对于一维单元，存在$\left| {\boldsymbol{J}} \right| {\text{＝}}{{\Delta L} \mathord{\left/ {\vphantom {{\Delta L} 2}} \right. } 2}$ ，其中ΔL为物理单元的尺度.由式（1）可见，通过对单元各节点的形函数进行积分即可得到单元质量矩阵中的各个元素，即
$$ {{M}}_{ij}^{\rm{e}} {\text{＝}} \int_{ - 1}^1 {\rho {N_i}（ \xi ）{N_j}（ \xi ）\left| {\boldsymbol{J}} \right|{\rm{d}}\xi } {\text{.}} $$ (2) 当对单元节点位置做不同选择时，相应的单元形函数
${N_i}（ \xi ）{\text{，}}i{\text{＝}} 1{\text{，}}2{\text{，}} \cdots {\text{，}}p$ ，自然亦有所不同，如此即形成不同类型谱单元的质量特性矩阵。目前已经发展出的谱单元模型主要有两种，即切比雪夫谱单元和勒让德谱单元，下面分别讨论.1.1 切比雪夫谱单元质量模型
切比雪夫谱单元选取GLC积分点作为单元节点。这组节点是多项式
$（ {1 {\text{－}} {\xi ^2}} ）{T'_{p {\text{－}} 1}}（ \xi ）$ 在 ［ −1，1 ］ 范围内的零点，其中Tp－1（ξ）为p－1阶第一类切比雪夫多项式。p－1阶切比雪夫谱单元的GLC节点位置表示如下$${\xi _i} {\text{＝}} - \cos \frac{{i\pi }}{{p {\text{－}} 1}} \qquad i {\text{＝}} 0{\text{，}}1{\text{，}} \cdots {\text{，}}p {\text{－}} 1{\text{，}}$$ (3) 式中，ξi是GLC点在标准区间上的坐标。切比雪夫谱单元可采用拉格朗日（Lagrange）插值函数构建单元形函数，即
$$ {N_i}（\xi ） {\text{＝}} \prod\limits_{\scriptstyle j {\text{＝}} 0\hfill\atop \scriptstyle j {\text{≠}} i\hfill}^{p {\text{－}} 1} {\frac{{\xi {\text{－}} {\xi _j}}}{{{\xi _i} {\text{－}} {\xi _j}}}} \qquad i {\text{＝}} 0{\text{，}}1{\text{，}} \cdots {\text{，}}p {\text{－}} 1{\text{.}} $$ (4) 由于GLC节点与切比雪夫正交多项式之间的联系，该拉格朗日形函数也可以通过一组截断的切比雪夫多项式进行描述，即
$$ {N_i}（ \xi ） {\text{＝}}\frac{2}{{p {\text{－}} 1}}\sum\limits_{k {\text{＝}} 0}^{p {\text{－}} 1} {\frac{1}{{{c_i}{c_k}}}{T_k}（ {{\xi _i}}）{T_k}（ \xi ）} \qquad i = 0{\text{，}}1{\text{，}} \cdots {\text{，}}p {\text{－}} 1 {\text{，}} $$ (5) 式中，系数ci＝2 （当i＝0或p－1时）或ci＝1 （当i＝1，···，p－2时）。根据多项式插值的唯一性可知，式（4）与式（5）本质相同，仅在多项式的计算过程中有所差异。式（5）形式的优点在于可以方便地利用切比雪夫正交多项式的性质精确求解单元特性矩阵。例如，k阶和l阶切比雪夫多项式的乘积在 ［ −1，1 ］ 上定积分Akl的解析解为
$$ {A}_{kl}{\text{＝}}{\displaystyle {\int }_{-1}^{1}{T}_{k}（\xi ）{T}_{l}（\xi ）{\rm{d}}\xi }{\text{＝}}\left\{\begin{array}{l}0 \qquad \qquad \qquad \qquad \qquad \qquad k{\text{＋}}l \;{\text{为奇数，}} \\ \dfrac{1}{1{\text{－}}{（k{\text{＋}}l）}^{2}}{\text{＋}}\dfrac{1}{1{\text{－}}{（k{\text{－}}l）}^{2}}\qquad k{\text{＋}}l \;{\text{为偶数，}} \end{array}\right. $$ (6) 利用式（6），可以方便地计算出式（2）中质量元素积分式的精确结果，即切比雪夫谱单元质量矩阵中的各个元素可以解析地写为
$$ M_{ij}^{\rm{e}} {\text{＝}} \rho \left| {\boldsymbol{J}} \right|\frac{4}{{{{（ {p {\text{－}} 1} ）}^2}{c_i}{c_j}}}\sum\limits_{k {\text{＝}} 0}^{p {\text{－}} 1} {\sum\limits_{l {\text{＝}} 0}^{p {\text{－}} 1} {\frac{1}{{{c_k}{c_l}}}{T_k}（ {{\xi _i}} ）{T_l}（ {{\xi _j}} ）{A_{kl}}} } {\text{.}} $$ (7) 显见，按照式（7）导出的切比雪夫谱单元的质量矩阵是非对角形式的CMM.
1.2 勒让德谱单元质量模型
勒让德谱单元采用GLL数值积分点配置谱单元节点。该组节点是多项式
$（ {1 {\text{－}} {\xi ^2}} ）{L'_{p {\text{－}} 1}}（ \xi ）$ 在 ［ −1，1 ］ 范围内的零点，这里Lp－1（ξ）为p－1阶勒让德多项式。勒让德谱单元同样采用拉格朗日插值函数建立单元形函数，即$$ {N_i}（ \xi ） {\text{＝}} \prod\limits_{\scriptstyle j {\text{＝}} 1\atop \scriptstyle j {\text{≠}} i}^p {\frac{{\xi {\text{－}} {\xi _j}}}{{{\xi _i} {\text{－}} {\xi _j}}}} {\text{，}} \qquad i {\text{＝}} 1{\text{，}}2{\text{，}} \cdots {\text{，}}p{\text{.}} $$ (8) 需要注意的是，式（8）中勒让德谱单元节点坐标ξi （i＝1，2，···，p）的选取与式（4）中切比雪夫谱单元的节点坐标是不同的，故两种谱单元的形函数亦不同。
同切比雪夫谱单元构建质量矩阵方式不同，勒让德谱单元采用了GLL数值积分，所选取的数值积分点与勒让德谱单元节点完全重合。该GLL数值积分的权系数按下式计算：
$$ {w_k} {\text{＝}} \left\{ {\begin{array}{*{20}{l}} {\dfrac{2}{{p（ {p {\text{－}} 1} ）}} \qquad \qquad\qquad\;\;\;\, k {\text{＝}} 1{\text{或}}\; k {\text{＝}}p{\text{，}}} \\ {\dfrac{2}{{p（ {p {\text{－}} 1} ）{{\left[ {{L_{p {\text{－}} 1}}（ {{\xi _k}} ）} \right]}^2}}} \qquad 2 {\text{≤}} k {\text{≤}} p {\text{－}} 1}{\text{.}} \end{array}} \right.$$ (9) 表1给出了1—5阶勒让德谱单元的GLL节点坐标和对应的积分权系数。根据数值积分规则，勒让德谱单元质量矩阵的各元素被近似计算为
表 1 GLL节点坐标和积分权系数Table 1. Abscissas of GLL points and quadrature weights谱单元阶次 GLL节点坐标 积分权系数 1 （−1，1） 1，1 2 （−1，0，1） 0.333 333，1.333 333，0.333 333 3 （−1，−0.447 214，0.447 214，1） 0.166 667，0.833 333，0.833 333，0.166 667 4 （−1，−0.654 654，0，0.654 654，1） 0.1，0.544 444，0.711 111，0.544 444，0.1 5 （−1，−0.765 055，−0.285 232，
0.285 232，0.765 055，1）0.066 667，0.378 475，0.554 858，
0.554 858，0.378 475，0.066 667$$ M_{ij}^{\rm{e}} {\text{＝}} \rho \left| {\boldsymbol{J}} \right|\sum\limits_{k {\text{＝}} 1}^p {{w_k}{N_i}（ {{\xi _k}} ）{N_j}（ {{\xi _k}} ）}{\text{，}} $$ (10) 由此导出的谱单元质量矩阵可自动形成LMM。
2. 两种质量模型的数学机制
2.1 积分方法对质量矩阵的影响
无论何种类型的谱单元，其质量特性模型均按照式（1）或式（2）的规则构建，需要对谱单元形函数进行积分计算。对于切比雪夫谱单元和勒让德谱单元，两者在导出质量特性矩阵过程中使用了不同的积分策略。切比雪夫谱单元利用解析积分计算质量矩阵，而勒让德谱单元在导出质量矩阵时采用了GLL数值积分。积分方式的不同最终导致两种谱单元的质量矩阵呈现不同的形式，即切比雪夫谱单元导出CMM，而勒让德谱单元导出LMM。
从积分精度分析可以看出，切比雪夫谱单元导出的质量矩阵为精确积分结果，而勒让德谱单元使用的GLL数值积分只具有（2p－3）阶代数精度。虽然高斯-洛巴托积分属于高精度数值积分，但对于拥有p个单元节点的（p－1）阶谱单元而言，被积函数（两个形函数的乘积）的阶次达到（2p－2）阶，故勒让德谱单元对质量模型的积分计算并非精确积分结果，导出的质量矩阵与式（1）的精确积分结果相比实际上存在一定误差。事实上，若采用精确积分计算勒让德谱单元质量矩阵，导出的质量矩阵也将是CMM而非LMM。
以一维2阶勒让德谱单元为例，设单元长度为∆L＝4，质量密度取为ρ＝1。在参考单元上，其GLL节点坐标为（−1，0，1），单元形函数写为
$$ {N_1}（ \xi ） {\text{＝}} \frac{1}{2}\xi （ {\xi {\text{－}} 1} ）{\text{，}}{N_2}（ \xi ） {\text{＝}} 1 {\text{－}} {\xi ^2}{\text{，}}{N_3}（ \xi ） {\text{＝}} \frac{1}{2}\xi （ {\xi {\text{＋}} 1} ） {\text{，}} $$ (11) 利用3点高斯积分计算该2阶勒让德谱单元的质量矩阵，积分点为±0.774 6和0，对应的积分权系数为0.555 6和0.888 9，得到质量矩阵如下：
$$ {{\boldsymbol{M}}^{\rm{e}}} {\text{＝}} \left[ {\begin{array}{*{20}{c}} \; {0.533\,3}&{0.266\,7}&{ - 0.133\,3} \\ \; {0.266\,7}&{2.133\,3}&\;\;\;{0.266\,7} \\ { - 0.133\,3}&{0.266\,7}&\;\;\;{0.533\,3} \end{array}} \right] {\text{，}}$$ (12) 显然，如此导出的勒让德谱单元质量矩阵为CMM。注意到高斯积分可以给出（2n－1）阶代数精度，其中n为高斯积分点数。具体到该2阶勒让德谱单元，有n＝3，而单元质量模型中被积函数的阶次为4阶，故通过高斯积分得到的形如式（12）的质量矩阵实际上是对式（1）的精确积分结果，不存在任何误差。类似地，对于切比雪夫谱单元，如果采用足够数量积分点的高斯积分（当2n－1≥2p－2时）计算式（2），导出的切比雪夫谱单元质量矩阵将与式（7）给出的结果完全相同，均为精确积分结果。也就是说，无论切比雪夫谱单元还是勒让德谱单元，如果按照式（1）或式（2）的精确积分结果确立谱单元质量矩阵，它们必然都是CMM形式。而传统勒让德谱单元质量矩阵为LMM形式，其实是其采用了非精确的高斯-洛巴托积分计算谱单元质量矩阵的结果。
另一方面，也可以采用勒让德谱单元所使用的高斯-洛巴托型数值积分导出切比雪夫谱单元的质量矩阵。注意到切比雪夫谱单元与勒让德谱单元配置单元节点的方式有所不同，所以相应的数值积分应为高斯-洛巴托-切比雪夫（GLC）积分，而非勒让德谱单元所使用的高斯-洛巴托-勒让德（GLL）积分。我们推导了这种基于GLC节点坐标的GLC数值积分形式，其1—5阶GLC点坐标和对应的积分权系数列于表2中。
表 2 GLC节点坐标和高斯-洛巴托积分权系数Table 2. Abscissas of GLC points and Gauss-Lobatto quadrature weights谱单元阶次 GLC节点坐标 积分权系数 1 （−1，1） 1，1 2 （−1，0，1） 0.333 333，1.333 333，0.333 333 3 （−1，−0.5，0.5，1） 0.111 111，0.888 889，0.888 889，0.111 111 4 （−1，−0.707 107，0，0.707 107，1） 0.066 667，0.533 333，0.8，0.533 333，0.066 667 5 （−1，−0.809 017，−0.309 017，
0.309 017，0.809 017，1）0.04，0.360 743，0.599 257，
0.599 257，0.360 743，0.04下面同样以上述的一维2阶谱单元为例考察谱单元质量矩阵情况，只是谱单元节点按照GLC点进行配置，即该谱单元变为切比雪夫谱单元，单元参数与上述算例相同。按照表2提供的GLC积分权系数计算导出切比雪夫谱单元的质量矩阵，结果如下：
$$ {{\boldsymbol{M}}^{\rm{e}}} {\text{＝}} \left[ {\begin{array}{*{20}{c}} {0.666\,7}&0&0 \\ 0&{2.666\,7}&0 \\ 0&0&{0.666\,7} \end{array}} \right] {\text{.}} $$ (13) 由此可见，当切比雪夫谱单元也采用高斯-洛巴托积分时，导出的质量矩阵同样为LMM。而GLC积分得到的亦是近似积分结果，由其导出的切比雪夫谱单元质量矩阵同精确积分结果同样存在误差。
2.2 质量特性模型的数学解释
前述可知，无论切比雪夫谱单元还是勒让德谱单元，只要采用精确的高斯积分计算质量特性模型，导出的质量矩阵就必然是CMM形式。而如果选择非精确的高斯-洛巴托积分（积分节点与谱单元节点一致）计算式（1），则将直接导出LMM形式的谱单元质量矩阵。也就是说，采取不同的积分方法会导出不同类型的谱单元质量矩阵。观察高斯-洛巴托积分方式，其积分节点选择为与谱单元节点一致，即对于切比雪夫谱单元和勒让德谱单元分别选为GLC点和GLL点。由于谱单元形函数是通过拉格朗日插值函数构建出来的，而拉格朗日函数具有克罗内克性质，即函数仅在定义点取值为1而在其他节点取值皆等于0，从而形函数具有如下的离散正交性：
$$ \sum\limits_{k {\text{＝}} 1}^p {{N_i}（ {{\xi _k}} ）{N_j}（ {{\xi _k}} ）} {\text{＝}} \left\{ {\begin{array}{*{20}{c}} {0\qquad \;\; i {\text{≠}} j}{\text{，}} \\ {1\qquad \;\; i {\text{＝}} j}{\text{.}} \end{array}} \right. $$ (14) 由式（2）可知，谱单元质量矩阵每个元素的积分式中都包含两个拉格朗日多项式的乘积，所以当数值积分点与单元节点完全重合时，由式（14）描述的谱单元形函数离散正交性质，导出的质量矩阵中必然只有对角线元素（被积函数中的两个形函数相同）非零而所有非对角元素皆为零。
以4阶谱单元为例标示了单元形函数与积分节点的关系，如图1所示。对于高斯-洛巴托积分，积分节点与谱单元节点完全一致，如果采用高斯积分，则高斯积分点与谱单元节点不重合。此种情况下各形函数在这些积分点的取值既不为0也不等于1，导致任意两个形函数的乘积在这些积分点处均不可能为0，从而使得质量矩阵被构建成为非对角的CMM。
3. 谱单元质量模型的特征
3.1 两种谱单元的集中质量分布特征
比较由节点积分法导出的切比雪夫谱单元和勒让德谱单元的集中质量分布特征。首先考虑一维标准单元，设单元质量密度为1，单元尺度ΔL＝2，则单元总质量为2。不同阶次的切比雪夫谱单元和勒让德谱单元的集中质量模型分布，如图2所示。由图可知，两种谱单元均能保证单元总质量不变，对于本例即各节点的质量之和始终等于2。但两种谱单元是按照不同的比例将总质量分配到各单元节点的。对于3节点谱单元（图2a），切比雪夫谱单元的质量分布情况与勒让德谱单元完全相同，这是由于此情况下GLC节点的坐标与GLL节点完全相同，且二者对应的高斯-洛巴托积分权系数也相等，对于4节点和5节点谱单元（图2b和图2c），切比雪夫谱单元内节点分配到的质量比勒让德谱单元相应内节点要多一些，而端部节点分配到的质量比勒让德谱单元要少。
![]()
图 2 一维集中质量切比雪夫谱单元（左）和勒让德谱单元（右）质量分布特征（a） 3节点单元；（b） 4节点单元；（c） 5节点单元Figure 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采用张量积形式构造二维和三维谱单元，故二维和三维谱单元的形函数依然具有克罗内克性质，不难证明利用高斯-洛巴托积分法同样能够保证导出的质量矩阵为LMM。下面以二维标准单元为例，比较切比雪夫谱单元与勒让德谱单元的集中质量分布特征。设单元大小为2×2，质量密度为1，则单元总质量为4。不同阶次的切比雪夫谱单元和勒让德谱单元的集中质量分布情况如图3所示。可以看出，两种二维谱单元同样能够保证单元总质量不变。与一维3节点单元类似，二维9节点的切比雪夫谱单元与勒让德谱单元的质量分布完全相同（图3a）。对于16节点和25节点谱单元，切比雪夫谱单元将更多的质量分配到内节点上，而勒让德谱单元则将更多质量分配给了边界节点。由图3c-f可见，切比雪夫谱单元的角节点分配到的质量大约只有勒让德谱单元的一半，且边节点的质量也比勒让德谱单元的小。
![]()
图 3 二维集中质量切比雪夫（左）和勒让德（右）谱单元质量分布特征（a） 9节点单元；（b） 16节点单元；（c） 25节点单元Figure 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 element3.2 两种质量特性模型比较
下面以切比雪夫谱单元为例，比较CMM与LMM的影响。由式（2）可以看出，除了单元质量密度ρ和雅可比行列式这两个常数项外，质量矩阵中的每个元素只与两个形函数的乘积在 ［ −1，1 ］ 上的定积分有关。该定积分的大小也可表示为 ［ −1，1 ］ 范围内被积函数与x轴所围成的区域面积。对于元素
$ M_{ij}^{\rm{e}} $ ，被积函数为定义在i点的形函数与定义在j点的形函数之积。因此，观察两个形函数乘积的图像便可直观地比较出质量矩阵中各元素的大小。一维4阶切比雪夫谱单元CMM中各元素的积分示意，如图4所示，每幅子图与该元素在CMM中的位置对应。例如，第二行第三列的子图代表CMM中第二行第三列质量元素对应的质量积分式，即$\displaystyle\int_{ - 1}^1 {{N_1}（ \xi ）{N_2}（ \xi ）{\rm{d}}\xi }$ 。图中阴影部分表示函数与x轴之间围成的面积，即质量积分式的大小。由图4可以看出，该图阵中处于主对角位置的子图阴影面积完全位于x轴上方，意味着切比雪夫谱单元CMM中主对角元素均为正值；而非对角线位置的子图中阴影面积可能位于x轴上方也可能位于x轴下方，说明CMM中的非对角元素可能是正值也可能是负值。总的来看，主对角元素的阴影面积最大，而距离主对角线越远则阴影面积越小。也就是说，在切比雪夫谱单元CMM中主对角元素的绝对值较大，而距离主对角线越远的质量元素绝对值越小。![]()
图 4 一维5节点切比雪夫谱单元一致质量模型数学图示Figure 4. Mathematical representation for consistent mass model of 1-D 5-node Chebyshev element对于由节点积分法导出的切比雪夫谱单元的LMM，主对角元素为
$$ M_{ii}^{\rm{e}} {\text{＝}} \rho \left| {\boldsymbol{J}} \right|\sum\limits_{k {\text{＝}} 0}^{p {\text{－}} 1} {{w_k}N_i^2（ {{\xi _k}} ）} {\text{，}} $$ (15) 根据形函数的克罗内克性质，式（15）进一步简化为
$$ M_{ii}^{\rm{e}} {\text{＝}} \rho \left| {\boldsymbol{J}} \right|{w_i} {\text{＝}} \rho \left| {\boldsymbol{J}} \right|\int_{ - 1}^1 {{N_i}（ \xi ）{\rm{d}}\xi } {\text{，}} $$ (16) 式中，wi为GLC积分的权系数
$$ {w_i} {\text{＝}} \int_{ - 1}^1 {{N_i}（ \xi ）{\rm{d}}\xi } {\text{，}} $$ (17) Ni（ξ）为定义在i点上的拉格朗日形函数。
由式（16）可见，除了ρ和
$\left| {\boldsymbol{J}} \right|$ 两个常数项外，切比雪夫谱单元的LMM中各质量元素只与形函数在 ［ −1，1 ］ 上的定积分有关。仍然以一维4阶切比雪夫谱单元为例，将各节点对应的形函数绘制成图像，如图5所示，各形函数与x轴之间围成的区域面积与各节点分配到的质量大小一一对应。不难看出，中间节点分配到的质量最大，而端节点分配到的质量最小，该规律与CMM中主对角元素之间的大小关系相同。![]()
图 5 一维5节点切比雪夫谱单元集中质量模型数学图示Figure 5. Mathematical representation for lumped mass model of 1-D 5-node Chebyshev element接下来通过特征值问题的求解，从数值计算的角度讨论切比雪夫谱单元LMM和CMM的差异，同时给出集中质量勒让德谱单元的计算结果。考察一根简支铁木辛柯（Timoshenko）梁的自由振动问题，设该梁长度L＝3 m，截面宽度b＝0.02 m，高度h＝0.1 m，质量密度ρ＝7 800 kg/m3，弹性模量E＝210 GPa。根据铁木辛柯梁理论建立切比雪夫谱单元，分别运用高斯-洛巴托积分和解析积分构造LMM和CMM两种质量模型，并计算其固有频率。
采用不同阶次的集中质量切比雪夫谱单元、一致质量切比雪夫谱单元和集中质量勒让德谱单元计算得到的简支梁的前6阶固有频率，列于表3。单元尺度固定为ΔL＝1 m，即简支梁被离散为3个谱单元。表中的解析解根据Craig和Kurdila （2006）中公式计算。从表中数据可见，两种切比雪夫谱单元计算结果的差距较小，均与解析解非常接近。随着谱单元阶次的提升，两种质量模型均能快速收敛于精确解。与集中质量勒让德谱单元相比，集中质量切比雪夫谱单元的计算结果也十分接近解析解，计算精度大致相当。对于简支梁的一阶频率，三种谱单元模型都只需划分3个4阶单元就已经能够给出相当准确的结果。
表 3 简支梁前6阶固有频率计算结果Table 3. First six natural frequencies of simply supported beam谱单元 单元阶次 固有频率 /Hz 1阶 2阶 3阶 4阶 5阶 6阶 切比雪夫
（集中质量）4 26.093 6 103.808 9 230.472 6 408.857 7 633.689 4 1 092.290 1 5 26.093 6 103.794 3 231.362 4 406.506 8 626.318 0 876.153 0 6 26.093 6 103.793 9 231.425 2 406.355 0 625.330 3 884.533 9 切比雪夫
（一致质量）4 26.093 7 103.818 0 231.485 4 410.983 5 647.442 7 1 025.303 2 5 26.093 6 103.794 3 231.453 3 406.625 7 627.393 8 887.768 5 6 26.093 6 103.793 9 231.418 6 406.356 5 625.379 6 886.020 8 勒让德
（集中质量）4 26.093 7 103.823 4 231.410 2 411.380 2 646.890 0 1 054.997 5 5 26.093 6 103.794 4 231.465 2 406.682 7 627.654 9 884.530 9 6 26.093 6 103.793 9 231.418 6 406.359 7 625.409 9 886.581 6 解析解 26.093 5 103.791 8 231.395 7 406.225 1 624.820 5 883.181 3 图6为采用不同数量的谱单元计算出的简支梁前6阶频率的均方根误差，此时三种谱单元阶次均保持4阶不变。由图可见，集中质量切比雪夫谱单元的误差曲线与一致质量切比雪夫谱单元的误差曲线十分接近。随着单元数量的增加，两种质量模型计算结果的误差都能迅速降低，当单元数量达到6个时，前6阶频率的均方根误差已降低至10−3以下。在某些情况下，如谱单元数等于4，5，7，8时，LMM给出的结果误差略低于CMM。与集中质量勒让德谱单元相比，集中质量切比雪夫谱单元除了单元数量取为3时的误差稍大以外，其余情况下的计算误差均不超过集中质量勒让德谱单元的误差。
![]()
图 6 不同谱单元数量下简支梁的前6阶固有频率计算误差Figure 6. Error of first six natural frequencies with different number of elements从该算例分析可以看出，对于特征值问题，集中质量切比雪夫谱单元，一致质量切比雪夫谱单元和集中质量勒让德谱单元均能取得较好的计算效果，三种谱单元模型不存在明显差距。
4. 关于质量模型物理机制的讨论
有限单元法中广泛使用的质量模型有两种，即呈对角形式的集中质量模型和呈非对角形式的一致质量模型。通常认为，一致质量模型是准确的，而集中质量模型是近似的。这是因为一致质量模型是严格按照变分原理导出的，而且通过精确积分结果忠实地反映了有限单元质量模型的真实特性，具有很强的数学基础。而集中质量模型的形成则有较多人为因素的干预，例如基于工程物理概念“堆聚”形成LMM，或者通过行和集中等途径人为地将CMM转换为LMM等，似乎都多少缺乏数学支撑。一些学者（Chan et al，1993；Jensen，1996；Wu，2006）通过对比分析指出，采用集中质量模型和一致质量模型得到的动力响应数值结果并无显著差别。但这些分析大多是从数值计算角度进行研究，事实上是默认一致质量模型更具合理性，而集中质量模型被认为只是同一致质量模型相差不大但使用更方便的一种质量模型。
质量是物质的最基本特性之一。根据经典力学原理，质量被认为是度量物体惯性大小的物理量，它是不变的，与物体运动速度无关。在每个有限单元上，实际的连续介质分布质量被离散至单元的各个节点上，也就是说，有限单元质量模型本质上是试图采用一种离散质量体系近似地代替实际的连续质量体系，用以描述连续介质有限单元的惯性特性。因此，判断有限单元质量模型优劣的标准应该是其描述实际连续介质惯性的能力。以四节点平面应变单元为例（图7），实际单元是由均匀连续介质构成。经过离散化处理后，单元惯性特性由分配于4个单元节点上的集中质量体现，它们反映力与运动之间的关系，由牛顿第二定律刻画。单元各节点质量之间通过弹性连接产生力学联系，这些弹性连接表征了单元的刚度特性，反映力与变形之间的关系，由虎克定律刻画。每个单元节点都具有两个自由度，分别用水平向位移u和竖向位移v进行描述。故该单元总计8个自由度，对应的刚度特性模型和质量特性模型均为8×8阶矩阵。根据牛顿第二定律，质量系数
$m_{ij}^{\rm{e}}$ 表征的物理含义为使j自由度产生单位加速度时需要在i自由度上施加的力的大小。这里i，j＝1，2，···，8，它们既可以是某个单元节点的水平向位移自由度u，也可以是竖向位移自由度v。考虑某个固定时刻t，假定在t之前单元所有节点都处于静止状态（波动未传播至该单元前即是这种状态），在t时刻单元第i自由度（例如v1自由度）上突然施加了力fi作用，则在i自由度上必然同时产生加速度ai，对应质点将沿i自由度方向发生运动（即质点m1沿v1方向运动）。注意到t时刻该质点尚未开始运动，需要经过一小段时间后该质点才能改变在i自由度方向上的静止状态。即产生位移具有滞后性，t时刻时单元各节点之间并未发生相对位移，此时各质点之间亦不存在相互作用的弹性力，因而在j自由度上对应质点当然仍然保持静止状态（例如此时质点m3在u3方向上处于静止）。也就是说，若在i自由度上受到力的作用，其不会在j自由度上立即就产生加速度，由此知有$m_{ij}^{\rm{e}} {\text{＝}} 0$ （当i ≠ j时）。即除了主对角线元素外，质量矩阵中非对角质量元素应该全部为零，单元质量矩阵应该为LMM形式。图8描绘了单元内的传力路径。施加在i自由度上的作用力fi瞬间引起质点沿i自由度的加速度ai。经过一小段时间∆t之后，质点沿i自由度运动到一个新的位置di。在此过程中，由于质点相对于其它质点的位置发生了改变且产生了相对位移，质点与其它质点之间的弹性连接发生变形，导致该质点与其余各质点之间产生弹性力kij。弹性连接的变形与弹性力之间的关系由本构方程决定。然后根据牛顿第二定律，施加在j自由度上的弹性力kij将会引起对应质点沿j自由度方向产生加速度aj。根据上述分析可知，i自由度方向的加速度ai与j自由度方向的加速度aj并不是同时产生的。也就是说，施加在i自由度上的作用力并不能瞬时引起j自由度方向的加速度。在上述单元运动分析中，质点之间的相互作用力通过弹性连接的变形进行传递，而弹性变形依赖于质点的相对运动，并且需要经过时间才能产生。该过程表明，在同一时刻不同质点的加速度并不耦合。这与静力情况下外力瞬时引起节点位移且节点力瞬间传递至所有节点是完全不同的。
此外，集中质量模型合理性还可以从波速有限原理（廖振鹏，2002）方面予以论证。如果质量矩阵的非对角元素不为零（即CMM形式），则在单元任意节点施加外力将会瞬时引起整个体系中所有节点都产生加速度。表明在一个子域中产生的波动将会瞬间传播至全域。显然，这与波速有限的物理原理相违背，因为波在介质中的传播需要时间，各质点按照距离波源的远近依次开始振动，不可能同时开始振动。根据这一原理，不同节点的加速度之间不应该存在耦合。换言之，质量矩阵的非对角项应该全部为零，即质量矩阵应该是对角形式的LMM，这样才能符合波动传播速度有限这一物理原理。
按照上述讨论，集中质量模型比一致质量模型在物理上更具合理性。但在实践中，采用集中质量模型得到的波动数值计算结果并非总是比一致质量模型更好（Tong et al，1971；Zienkiewicz et al，2013）。前面已经论述，无论是一致质量矩阵CMM还是集中质量矩阵LMM，它们实际上都是有限单元质量模型的积分结果，区别仅在于各自使用了不同的积分方法。由于一致质量模型采用代数精度最高的高斯数值积分计算有限单元质量矩阵，当选取足够数量积分点数时能够获得精确积分结果，因而一般认为其优于集中质量模型。然而，关于一致质量模型优于集中质量模型的认识需要一个前提，即如果按照式（1）构建的有限单元质量模型是“精确”的，则作为精确积分结果的CMM必然亦是“精确”的。虽然高斯-洛巴托积分也是一种高精度数值积分，但通过其导出的LMM与有限单元质量模型精确积分结果相比存在一定误差，从这个意义上看CMM要好于LMM。但是，有限单元质量模型本质上是对原连续介质质量特性的一种近似化处理，其本身就不是“精确”的，在力学特性描述上不可避免地会带来误差。这种离散化和积分过程对单元质量特性造成的影响，如图9所示。如果能够精确计算质量矩阵（即形成CMM），则由数值积分引起的误差为零。由图9可见，精确积分之后的误差界与离散化之后的误差界一致，也就是说精确积分不会带来额外的误差。但是，如果质量矩阵没有被精确计算，如使用了高斯-洛巴托积分并导出LMM，则由数值积分引入的误差可能为正也可能为负。因而，关于质量特性的误差有可能被扩大也可能被减小。图9中高斯-洛巴托积分之后的误差界（虚线），它即可能在实线以内也可能在实线以外。浅灰色区域表示描述质量特性的误差被高斯-洛巴托积分减小，说明在此情况下LMM要优于CMM；深灰色区域表示误差被高斯-洛巴托积分放大，这表明此情况下LMM对质量特性的描述比CMM要差。因此，不是十分精确的积分结果并不总是造成不利影响，它有可能会改善离散化造成的误差，当然亦可能会增大离散化影响从而导致更大的误差。这可能也是基于LMM的计算结果在有些情况下好于CMM的结果而在另外一些情况下表现更差的原因。不过，由于高斯-洛巴托积分精度只比精确的高斯-勒让德积分略低一点，它们都属于高精度数值积分方法，所以由于数值积分引起的误差依然会被控制在有限范围内。
![]()
图 9 离散化和积分过程引起的质量特性误差Figure 9. Error in mass property caused by discretization and quadrature事实上，动力分析结果不仅取决于对原连续体系惯性特性的模拟，同时还取决于对弹性特性估计的准确性，两者具有耦合影响。一般认为，在有限元法中单元刚度矩阵高估了所模拟的连续体系的实际刚度（Bathe，2014）。一种直观的理解是有限元离散化过程将原来的无限自由度体系变为了有限自由度体系，这相当于人为地给连续体系施加了若干约束，从而减小了体系的变形。偏刚的估计误差在某些情况下甚至可能造成计算结果不可接受，例如梁、板分析中遇到的剪切自锁问题（王勖成，2003）。解决问题的一种途径是采用减缩积分计算单元刚度矩阵。减缩积分的思想是通过不精确的数值积分计算结果适当弥补有限元离散化造成的刚度偏大的误差（Bathe，2014），从而改善总体计算精度。但当单元阶次较低时，减缩积分的实施有可能导致单元过柔，从而出现虚假的零能模式（王勖成，2003）。这种处理措施同本文中利用高斯-洛巴托积分形成谱单元集中质量矩阵的过程相类似，即采用不精确的积分有可能弥补离散化引起的误差，但也可能会增大离散化误差。与有限元法相比，目前关于时域谱元模型所描述的体系刚度特性的误差估计尚无明确结论，但通过数值积分误差来适当地调控谱单元离散化带来的误差，无疑是一条可以考虑的途径。另外，由于静力问题无需考虑力与运动的传递过程，故在有限元法中只是简单地通过对位移场进行求导而获得加速度场，而没有完全反映出体系中力与运动传递的物理过程。这是动力问题与静力问题的区别所在，也是本文认为集中质量模型更加合理的基础。
5. 结论
集中质量矩阵求逆计算十分简便，这对于大规模动力分析问题具有非常重要的意义。本文将时域谱元法的两种质量特性模型—集中质量模型和一致质量模型都统一在相同的数学框架下，给出了统一的构建方法，并从数学机制和物理意义两个方面阐述了谱单元质量模型问题，结论如下：
1） 谱单元质量特性矩阵的构建归结为对谱单元形函数积的积分运算，选择不同的数值积分方法将导致不同形式的谱单元质量矩阵。当采用高斯-勒让德积分计算谱单元质量模型的定积分式时将导出一致质量矩阵，而选择高斯-洛巴托积分时则会导出集中质量矩阵。该结论对于切比雪夫谱单元和勒让德谱单元都是适用的。
2） 采用不同积分方法导出不同形式质量矩阵的根本原因在于质量矩阵中的被积形函数在积分点处是否具有克罗内克性质。高斯-洛巴托积分节点与谱单元节点一致，在积分点处谱单元形函数具有克罗内克性质，其导出的质量矩阵必为对角阵；而高斯积分无法满足积分点处谱单元形函数的克罗内克函数条件，只能导出一致质量矩阵。这是数值积分方法决定谱单元质量矩阵形式的数学机制。
3） 波动在介质中传播需要时间，在某个自由度上的力作用只能改变其自身自由度上的运动状态，而不会在其他自由度上立即产生加速度。从质量的物理意义上看，不同自由度上的质量系数不存在耦合，谱单元的集中质量模型更具有物理合理性。
4） 当选择的积分节点数与谱单元节点数相同时，高斯积分给出的是精确积分结果，而高斯-洛巴托积分结果存在误差。由于谱单元质量模型本身亦是对原连续介质惯性特性的近似描述，数值积分误差有可能放大亦可能降低谱单元离散化引起的误差范围。
5） 高斯-洛巴托积分是一种高精度数值积分，其代数精度仅略低于高斯积分，由其导出的谱单元集中质量矩阵对谱单元离散化误差的调整被限制在有限范围内。与高斯积分导出的谱单元一致质量矩阵相较，两种质量模型对于动力分析效果而言大致相当，都能取得很好的计算结果。
-
![]()
图 2 一维集中质量切比雪夫谱单元（左）和勒让德谱单元（右）质量分布特征
（a） 3节点单元；（b） 4节点单元；（c） 5节点单元
Figure 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
![]()
图 3 二维集中质量切比雪夫（左）和勒让德（右）谱单元质量分布特征
（a） 9节点单元；（b） 16节点单元；（c） 25节点单元
Figure 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
![]()
图 4 一维5节点切比雪夫谱单元一致质量模型数学图示
Figure 4. Mathematical representation for consistent mass model of 1-D 5-node Chebyshev element
![]()
图 5 一维5节点切比雪夫谱单元集中质量模型数学图示
Figure 5. Mathematical representation for lumped mass model of 1-D 5-node Chebyshev element
![]()
图 6 不同谱单元数量下简支梁的前6阶固有频率计算误差
Figure 6. Error of first six natural frequencies with different number of elements
![]()
图 9 离散化和积分过程引起的质量特性误差
Figure 9. Error in mass property caused by discretization and quadrature
表 1 GLL节点坐标和积分权系数
Table 1 Abscissas of GLL points and quadrature weights
谱单元阶次 GLL节点坐标 积分权系数 1 （−1，1） 1，1 2 （−1，0，1） 0.333 333，1.333 333，0.333 333 3 （−1，−0.447 214，0.447 214，1） 0.166 667，0.833 333，0.833 333，0.166 667 4 （−1，−0.654 654，0，0.654 654，1） 0.1，0.544 444，0.711 111，0.544 444，0.1 5 （−1，−0.765 055，−0.285 232，
0.285 232，0.765 055，1）0.066 667，0.378 475，0.554 858，
0.554 858，0.378 475，0.066 667
下载: 导出CSV
表 2 GLC节点坐标和高斯-洛巴托积分权系数
Table 2 Abscissas of GLC points and Gauss-Lobatto quadrature weights
谱单元阶次 GLC节点坐标 积分权系数 1 （−1，1） 1，1 2 （−1，0，1） 0.333 333，1.333 333，0.333 333 3 （−1，−0.5，0.5，1） 0.111 111，0.888 889，0.888 889，0.111 111 4 （−1，−0.707 107，0，0.707 107，1） 0.066 667，0.533 333，0.8，0.533 333，0.066 667 5 （−1，−0.809 017，−0.309 017，
0.309 017，0.809 017，1）0.04，0.360 743，0.599 257，
0.599 257，0.360 743，0.04
下载: 导出CSV
表 3 简支梁前6阶固有频率计算结果
Table 3 First six natural frequencies of simply supported beam
谱单元 单元阶次 固有频率 /Hz 1阶 2阶 3阶 4阶 5阶 6阶 切比雪夫
（集中质量）4 26.093 6 103.808 9 230.472 6 408.857 7 633.689 4 1 092.290 1 5 26.093 6 103.794 3 231.362 4 406.506 8 626.318 0 876.153 0 6 26.093 6 103.793 9 231.425 2 406.355 0 625.330 3 884.533 9 切比雪夫
（一致质量）4 26.093 7 103.818 0 231.485 4 410.983 5 647.442 7 1 025.303 2 5 26.093 6 103.794 3 231.453 3 406.625 7 627.393 8 887.768 5 6 26.093 6 103.793 9 231.418 6 406.356 5 625.379 6 886.020 8 勒让德
（集中质量）4 26.093 7 103.823 4 231.410 2 411.380 2 646.890 0 1 054.997 5 5 26.093 6 103.794 4 231.465 2 406.682 7 627.654 9 884.530 9 6 26.093 6 103.793 9 231.418 6 406.359 7 625.409 9 886.581 6 解析解 26.093 5 103.791 8 231.395 7 406.225 1 624.820 5 883.181 3
下载: 导出CSV
-
丁志华, 周红, 蒋涵. 2014. 三维台阶地形地震动效应研究[J]. 地震学报, 36（2）: 184-199 doi: 10.3969/j.issn.0253-3782.2014.02.004 Ding Z H, Zhou H, Jiang H. 2014. Effect of 3-D step topography on ground motion[J]. Acta Seismologica Sinica, 36（2）: 184-199 （in Chinese）. doi: 10.3969/j.issn.0253-3782.2014.02.004
韩天成, 于彦彦, 丁海平. 2020. 直下型断层的破裂速度对盆地地震效应的影响[J]. 地震学报, 42（4）: 457-470 doi: 10.11939/jass.20190177 Han T C, Yu Y Y, Ding H P. 2020. Influence of rupture velocity of the directly-beneath fault on the basin seismic effect[J]. Acta Seismologica Sinica, 42（4）: 457-470 （in Chinese）. doi: 10.11939/jass.20190177
廖振鹏. 2002. 工程波动理论导论[M]. 第2版. 北京: 科学出版社: 59–63 Liao Z P. 2002. Introduction to Wave Motion Theories in Engineering[M]. 2nd ed. Beijing: Science Press: 59–63 （in Chinese）.
王勖成. 2003. 有限单元法[M]. 北京: 清华大学出版社: 472–475 Wang X C. 2003. Finite Element Method[M]. Beijing: Tsinghua University Press: 472–475 （in Chinese）.
邢浩洁. 2017. 透射边界机理及其在地震波动谱元模拟中的应用[D]. 南京: 南京工业大学: 39–82 Xing H J. 2017. Mechanism of Transmitting Boundary and Its Application to Seismic Wave Simulation With Spectral Element Method[D]. Nanjing: Nanjing Tech University: 39–82 （in Chinese）.
Bathe K J. 2014. Finite Element Procedures[M]. 2nd ed. USA: Prentice-Hall: 476–480.
Bottero A, Cristini P, Komatitsch D, Asch M. 2016. An axisymmetric time-domain spectral-element method for full-wave simulations: Application to ocean acoustics[J]. J Acoust Soc Am, 140（5）: 3520-3530. doi: 10.1121/1.4965964
Chan H C, Cai C W, Cheung Y K. 1993. Convergence studies of dynamic analysis by using the finite element method with lumped mass matrix[J]. J Sound Vib, 165（2）: 193-207. doi: 10.1006/jsvi.1993.1253
Craig Jr R R, Kurdila A J. 2006. Fundamentals of Structural Dynamics[M]. 2nd ed. Hoboken, New Jersey: John Wiley & Sons: 400–401.
Cristini P, Komatitsch D. 2012. Some illustrative examples of the use of a spectral-element method in ocean acoustics[J]. J Acoust Soc Am, 131（3）: EL229-EL235. doi: 10.1121/1.3682459
Dauksher W, Emery A F. 1997. Accuracy in modeling the acoustic wave equation with Chebyshev spectral finite elements[J]. Finite Elem Anal Des, 26（2）: 115-128. doi: 10.1016/S0168-874X(96)00075-3
Dauksher W, Emery A F. 1999. An evaluation of the cost effectiveness of Chebyshev spectral and p-finite element solutions to the scalar wave equation[J]. Int J Numer Methods Eng, 45（8）: 1099-1113. doi: 10.1002/(SICI)1097-0207(19990720)45:8<1099::AID-NME622>3.0.CO;2-5
Dauksher W, Emery A F. 2000. The solution of elastostatic and elastodynamic problems with Chebyshev spectral finite elements[J]. Comput Methods Appl Mech Eng, 188（1/2/3）: 217-233.
Duczek S, Gravenkamp H. 2019a. Mass lumping techniques in the spectral element method: On the equivalence of the row-sum, nodal quadrature, and diagonal scaling methods[J]. Comput Methods Appl Mech Eng, 353: 516-569. doi: 10.1016/j.cma.2019.05.016
Duczek S, Gravenkamp H. 2019b. Critical assessment of different mass lumping schemes for higher order serendipity finite elements[J]. Comput Methods Appl Mech Eng, 350: 836-897. doi: 10.1016/j.cma.2019.03.028
Fried I, Malkus D S. 1975. Finite element mass matrix lumping by numerical integration with no convergence rate loss[J]. Int J Solids Struct, 11（4）: 461-466. doi: 10.1016/0020-7683(75)90081-5
Hinton E, Rock T, Zienkiewicz O C. 1976. A note on mass lumping and related processes in the finite element method[J]. Earthq Eng Struct Dyn, 4（3）: 245-249. doi: 10.1002/eqe.4290040305
Hughes T J R. 1987. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[M]. Englewood Cliffs, NJ: Prentice-Hall: 436−446.
Jensen M S. 1996. High convergence order finite elements with lumped mass matrix[J]. Int J Numer Methods Eng, 39（11）: 1879-1888. doi: 10.1002/(SICI)1097-0207(19960615)39:11<1879::AID-NME933>3.0.CO;2-2
Komatitsch D, Vilotte J P. 1998. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures[J]. Bull Seismol Soc Am, 88（2）: 368-392.
Komatitsch D, Tromp J. 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation[J]. Geophys J Int, 139（3）: 806-822. doi: 10.1046/j.1365-246x.1999.00967.x
Kudela P, Krawczuk M, Ostachowicz W. 2007a. Wave propagation modelling in 1D structures using spectral finite elements[J]. J Sound Vib, 300（1/2）: 88-100.
Kudela P, Żak A, Krawczuk M, Ostachowicz W. 2007b. Modelling of wave propagation in composite plates using the time domain spectral element method[J]. J Sound Vib, 302（4/5）: 728-745.
Patera A T. 1984. A spectral element method for fluid dynamics: Laminar flow in a channel expansion[J]. J Comput Phys, 54（3）: 468-488. doi: 10.1016/0021-9991(84)90128-1
Priolo E, Carcione J M, Seriani G. 1994. Numerical simulation of interface waves by high-order spectral modeling techniques[J]. J Acoust Soc Am, 95（2）: 681-693. doi: 10.1121/1.408428
Tong P, Pian T H H, Bucciarblli L L. 1971. Mode shapes and frequencies by finite element method using consistent and lumped masses[J]. Comput Struct, 1（4）: 623-638. doi: 10.1016/0045-7949(71)90033-2
Wu S R. 2006. Lumped mass matrix in explicit finite element method for transient dynamics of elasticity[J]. Comput Methods Appl Mech Eng, 195（44/47）: 5983-5994.
Yang Y T, Zheng H, Sivaselvan M V. 2017. A rigorous and unified mass lumping scheme for higher-order elements[J]. Comput Methods Appl Mech Eng, 319: 491-514. doi: 10.1016/j.cma.2017.03.011
Żak A. 2009. A novel formulation of a spectral plate element for wave propagation in isotropic structures[J]. Finite Elem Anal Des, 45（10）: 650-658. doi: 10.1016/j.finel.2009.05.002
Żak A, Krawczuk M, Palacz M, Doliński Ł, Waszkowiak W. 2017. High frequency dynamics of an isotropic Timoshenko periodic beam by the use of the Time-domain Spectral Finite Element Method[J]. J Sound Vib, 409: 318-335. doi: 10.1016/j.jsv.2017.07.055
Żak A, Krawczuk M. 2018. A higher order transversely deformable shell-type spectral finite element for dynamic analysis of isotropic structures[J]. Finite Elem Anal Des, 142: 17-29. doi: 10.1016/j.finel.2017.12.007
Zhang G H, Yang Y T, Zheng H. 2019a. A mass lumpig scheme for the second-order numerical manifold method[J]. Comput Struct, 213: 23-39. doi: 10.1016/j.compstruc.2018.12.005
Zhang G H, Yang Y T, Sun G H, Zheng H. 2019b. A mass lumping scheme for the 10-node tetrahedral element[J]. Eng Anal Bound Elem, 106: 190-200. doi: 10.1016/j.enganabound.2019.04.018
Zheng H, Yang Y T. 2017. On generation of lumped mass matrices in partition of unity based methods[J]. Int J Numer Methods Eng, 112（8）: 1040-1069. doi: 10.1002/nme.5544
Zhu C Y, Qin G L, Zhang J Z. 2011. Implicit Chebyshev spectral element method for acoustics wave equations[J]. Finite Elem Anal Des, 47（2）: 184-194. doi: 10.1016/j.finel.2010.09.004
Zienkiewicz O C, Taylor R L, Zhu J Z. 2013. The Finite Element Method: Its Basis and Fundamentals[M]. 7th ed. London: Elsevier: 383–386.
-
期刊类型引用(0)
其他类型引用(1)
计量
- 文章访问数: 375
- HTML全文浏览量: 199
- PDF下载量: 41
- 被引次数: 1
