PRECISE MID-NODE LUMPED MASS MATRICES FOR 3D 20-NODE HEXAHEDRAL AND 10-NODE TETRAHEDRAL FINITE ELEMENTS
-
摘要: 对于工程分析中常用的三维20节点六面体和10节点四面体单元, 采用行求和法得到的集中质量矩阵由于包含负质量元素, 难以直接用于动力分析. 虽然常用的主对角元素放大法, 即HRZ方法, 可以有效规避负质量元素, 但是该方法仍缺乏理论层面的精度分析. 本文首先以三维20节点六面体单元为例, 构造一种包含待定参数的广义集中质量矩阵, 并将HRZ集中质量矩阵作为特例涵盖其中, 进而建立了20节点六面体单元广义集中质量矩阵的频率精度表达式. 然后, 通过参数优化, 提出20节点六面体单元的中节点集中质量矩阵构造方法, 并从理论上证明其精度优于HRZ集中质量矩阵. 该中节点集中质量矩阵形式简单, 非常便于推广到10节点四面体单元. 此外, 利用中节点集中质量矩阵含有主对角零质量元素的特点, 通过静力凝聚建立了相应的动力分析降阶模型, 在保证计算精度的同时可大幅提升计算效率. 自由振动和时程分析结果均表明, 对于三维20节点六面体单元和10节点四面体单元, 中节点集中质量矩阵的计算精度明显高于HRZ集中质量矩阵.Abstract: The three-dimensional (3D) 20-node hexahedral and 10-node tetrahedral elements have been widely used in structural analysis, while the row-sum mass lumping for such elements leads to certain negative entries for the lumped mass matrices, which cause considerable difficulty for dynamic finite element analysis. Meanwhile, despite that the diagonal scaling mass lumping technique, namely, the HRZ method, ensures a non-negative lumped mass formulation, a theoretical accuracy investigation of this frequently employed approach is still missed in 3D scenarios. In this work, a generalized lumped mass matrix template with specifically devised adjustable parameters is introduced for 3D 20-node hexahedral elements to assess the frequency accuracy of different mass lumping methods, which include the HRZ lumped mass matrix formulation as a specific circumstance. Subsequently, a frequency accuracy measure is rationally derived for this generalized lumped mass matrix template. With the aid of the frequency error measure, it is found that the HRZ lumped mass matrix formulation for 3D 20-node hexahedral elements does not give the optimal frequency accuracy. On the other hand, a precise mid-node lumped mass matrix formulation is attained by optimizing the frequency accuracy with respect to the adjustable parameters. A straightforward extension of the proposed formulation immediately yields an accurate mid-node lumped mass matrix formulation for 10-node tetrahedral elements. Furthermore, since the corner nodes of the mid-node lumped mass matrices for 20-node hexahedral and 10-node tetrahedral elements have zero diagonal mass entries, a standard static condensation operation on the discrete finite element equations enables a computationally efficient reduced order model that only contains the mid-node degrees of freedom for subsequent dynamic analysis. The frequency computation and transient analysis results consistently buttress that compared to the conventional HRZ lumped mass matrices, the proposed mid-node lumped mass matrices for 3D 20-node hexahedral and 10-node tetrahedral elements exhibit salient advantages regarding both accuracy and efficiency.
-
引言
在结构动力分析[1-6]中, 质量矩阵的构造形式对有限元计算精度具有显著的影响. 其中, 常用的有限元质量矩阵构造形式主要包括一致质量矩阵[1,7-8]、高阶质量矩阵[1,9-11]和集中质量矩阵[12-19]. 其中, 集中质量矩阵具有构造形式简单、存储空间小、计算效率高和可与中心差分法等显式积分方案结合使用等优点. 而常用的集中质量矩阵构造方法主要包括行求和法[1,13]、节点积分法[14-16]和一致质量矩阵主对角元素放大法(HRZ法)[17]等. 在集中质量矩阵的精度分析方面, Ainsworth等[16]给出拉格朗日单元采用节点积分方案求解标量波特征值问题时的精度度量理论. Yang等[18]针对8节点平面单元提出一种基于数值微分流形的集中质量构造方法. Duczek等[19]研究了Lobatto单元行求和、节点积分和HRZ集中质量矩阵之间的等价性. 最近, Li等[20]研究了采用行求和集中质量矩阵时, 拉格朗日单元节点布置对频率计算精度和收敛阶次的影响, 详细分析了节点均匀分布单元集中质量矩阵的频率精度受限性.
在有限元分析领域, 20节点六面体单元和10节点四面体单元应用非常广泛[21-25]. 与27节点六面体拉格朗日单元相比, 20节点六面体单元不含内部节点, 充分利用了单元间的节点共享特性, 可显著缩减计算模型的带宽, 提高计算效率. 与此同时, 10节点四面体单元在复杂形状网格剖分方面具有明显优势. 然而, 当采用行求和方法构造这两类单元的集中质量矩阵时, 主对角线上均出现负质量元素, 难以直接用于结构动力分析. 目前, 最常用的消除集中矩阵负质量元素的有效方法是HRZ法[17], 但是三维情况下该方法缺乏理论层面的精度分析. 针对8节点平面单元的HRZ集中质量矩阵, Hou等[26]建立了相应的频率精度表达式, 通过分析证明HRZ集中质量矩阵并不能提供最优的频率计算精度, 进一步发展了具有更优精度的8节点平面单元中节点集中质量矩阵. 但是, 该研究尚未考虑更为复杂的三维单元. 另外, 由于中节点集中质量矩阵的角节点质量为0[26], 其对时程动力分析的影响也需要进一步厘清.
本文针对三维20节点六面体和10节点四面体单元行求和集中质量矩阵出现负质量的问题, 提出了三维广义中节点集中质量矩阵的构造形式, 并将HRZ集中质量矩阵作为特例涵盖其中. 然后, 以20节点六面体单元为例, 构建了三维广义中节点集中质量矩阵的频率计算精度表达式, 进而分析了三维HRZ集中质量矩阵的频率精度. 根据三维广义中节点集中质量矩阵的频率精度表达式, 通过对其中的待定参数进行优化, 构造了20节点六面体单元的高精度集中质量矩阵. 接着, 利用中节点集中质量矩阵构造形式简单的特点, 将其推广到10节点四面体单元. 最后, 对于时程动力分析, 通过静力凝聚消除了中节点集中质量矩阵的角节点零质量影响, 构造了三维有限元中节点集中质量矩阵的降阶动力计算模型, 在保证动力计算精度的同时提高了计算效率.
1. 结构动力分析的有限元离散方程
不失一般性, 考虑如下的经典波动方程等效积分弱形式[1]
$$ \int_\varOmega {\delta {{\boldsymbol{u}}^{\rm{T}}}{\boldsymbol{\ddot u}}{\rm{d}}\varOmega } + \int_\varOmega {{({\overset{\frown} {\boldsymbol{\nabla}} } \delta {\boldsymbol{u}})^{\rm{T}}}{\boldsymbol{D}}({\overset{\frown} {\boldsymbol{\nabla}} } {\boldsymbol{u}}){\rm{d}}\varOmega } = \int_\varOmega {\delta {{\boldsymbol{u}}^{\rm{T}}}{\boldsymbol{b}}{\rm{d}}\varOmega } $$ (1) 其中, $ \delta $表示变分符号, 上标T代表转置. 对于膜或声场问题, 场变量$ {\boldsymbol{u}} $退化为标量$ u $, 表示膜的横向位移或声压, 外力或源项$ {\boldsymbol{b}} $退化为标量$ b $. 与之对应的梯度矩阵$\overset{\frown} {\boldsymbol{\nabla}}$和材料矩阵D分别为
$$ \overset{\frown} {\boldsymbol{\nabla}} = {\left\{ \frac{\partial }{{\partial x}},\frac{\partial }{{\partial y}},\;\frac{\partial }{{\partial z}}\;\right\} ^{\rm{T}}},\;{\boldsymbol{D}} = {c^2} $$ (2) 其中, $ c $为波速. 对于弹性力学问题, 矢量场${\boldsymbol{u}} = {\{ {u_x},\;{u_y},\;{u_z}\} ^{\rm{T}}}$表示弹性体内一点的位移, ${\boldsymbol{b}} = {\{ {b_x},\;{b_y},\;{b_z}\} ^{\rm{T}}}$表示体力, 相应的梯度矩阵$\overset{\frown} {\boldsymbol{\nabla}}$和材料弹性矩阵$ {\boldsymbol{D}} $分别为
$$ \overset{\frown} {\boldsymbol{\nabla}} = {\left[ {\begin{array}{*{20}{c}} {\dfrac{\partial }{{\partial x}}}&0&0&{\dfrac{\partial }{{\partial y}}}&0&{\dfrac{\partial }{{\partial z}}} \\ 0&{\dfrac{\partial }{{\partial y}}}&0&{\dfrac{\partial }{{\partial x}}}&{\dfrac{\partial }{{\partial z}}}&0 \\ 0&0&{\dfrac{\partial }{{\partial z}}}&0&{\dfrac{\partial }{{\partial y}}}&{\dfrac{\partial }{{\partial x}}} \end{array}} \right]^{\rm{T}}} $$ (3) $$ {\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} {c_p^2}&{c_p^2 - 2c_s^2}&{c_p^2 - 2c_s^2}&0&0&0 \\ {c_p^2 - 2c_s^2}&{c_p^2}&{c_p^2 - 2c_s^2}&0&0&0 \\ {c_p^2 - 2c_s^2}&{c_p^2 - 2c_s^2}&{c_p^2}&0&0&0 \\ 0&0&0&{c_s^2}&0&0 \\ 0&0&0&0&{c_s^2}&0 \\ 0&0&0&0&0&{c_s^2} \end{array}} \right] $$ (4) 其中, $ {c_p} $和$ {c_s} $表示压力波和剪力波的波速[11]
$$ {c_p} = \sqrt {(\lambda + 2\mu )/\rho } ,\;\;{c_s} = \sqrt {\mu /\rho } $$ (5) 其中, $ \lambda $和$ \mu $为拉梅常数, $ \rho $为材料密度.
当采用20节点六面体单元或10节点四面体单元对计算区域$ \varOmega $进行离散, 即$\varOmega = \cup _{e = 1}^{{n_{el}}}{\varOmega ^e}$, 其中$ {n_{el}} $表示单元数, 未知变量$ {\boldsymbol{u}} $的有限元近似值$ {{\boldsymbol{u}}^h} $可表示为
$$ {{\boldsymbol{u}}^h}({\boldsymbol{x}},t) = \sum\limits_{I = 1}^{{n_{en}}} {N_I^{}} ({\boldsymbol{\xi}} ){\boldsymbol{d}}_I^e(t),\;{\text{ }}{\boldsymbol{x}} = \sum\limits_{I = 1}^{{n_{en}}} {N_I^{}} ({\boldsymbol{\xi}} ){{\boldsymbol{x}}_I} $$ (6) 其中, ${N_I}({{\boldsymbol{\xi}} })$表示单元第$ I $个节点的形函数, $ {n_{en}} $为单元节点个数. 将式(6)代入式(1)中, 可得结构动力分析的有限元离散方程
$$ {\boldsymbol{M\ddot d}} + {\boldsymbol{Kd}} = {\boldsymbol{f}} $$ (7) $$ {\boldsymbol{M}} = \mathop \mathcal{A}\limits_{e = 1}^{{n_{el}}} {{\boldsymbol{M}}^e},{\text{ }}{\boldsymbol{K}} = \mathop \mathcal{A}\limits_{e = 1}^{{n_{el}}} {{\boldsymbol{K}}^e},{\text{ }}{\boldsymbol{f}} = \mathop \mathcal{A}\limits_{e = 1}^{{n_{el}}} {{\boldsymbol{f}}^e} $$ (8) 其中, $ {\boldsymbol{M}} $和K分别为整体质量矩阵和整体刚度矩阵, $ {\boldsymbol{f}} $为力向量, 可分别由单元质量矩阵$ {{\boldsymbol{M}}^e} $、刚度矩阵$ {{\boldsymbol{K}}^e} $和力向量$ {{\boldsymbol{f}}^e} $通过组装算子$ \mathcal{A} $组装得到. $ {{\boldsymbol{M}}^e} $, $ {{\boldsymbol{K}}^e} $和$ {{\boldsymbol{f}}^e} $的定义为
$$ {\boldsymbol{M}}_{IJ}^e = \int_{{\varOmega ^e}} {{N_I}} {N_J}{\boldsymbol{1}}{\rm{d}}\varOmega $$ (9) $$ {\boldsymbol{K}}_{IJ}^e = \int_{{\varOmega ^e}} {{{(\overset{\frown} {\boldsymbol{\nabla}} {N_I})^{\rm{T}}}}} {\boldsymbol{D}}(\overset{\frown} {\boldsymbol{\nabla}} {{{N}}_J}){\rm{d}}\varOmega $$ (10) $$ {\boldsymbol{f}}_I^e = \int_{{\varOmega ^e}} {{N_I}} {\boldsymbol{b}}{\rm{d}}\varOmega $$ (11) 其中, $ {\boldsymbol{1}} $是单位矩阵, 对于声场问题, $ {\boldsymbol{1}} = 1 $; 对于弹性连续体问题, ${\boldsymbol{1}} = {\rm{diag}}\{ 1,\;1,\;1\}$.
对于自由振动分析, 在式(7)中忽略外力项, 并引入如下的简谐波假定
$$ {\boldsymbol{d}}(t) = {\boldsymbol{\phi}} \exp (\text{ι} {\omega ^h}t),\;\;\text{ι} = \sqrt { - 1} $$ (12) 其中, $ {\omega ^h} $为数值频率, ${\boldsymbol{\phi }}$为与之对应的振动模态. 将式(12)代入式(7), 可得自由振动分析的有限元离散方程
$$ \left[ {\boldsymbol{K}} - {({\omega ^h})^{\text{2}}}{\boldsymbol{M}}\right] {\boldsymbol{\phi}} = {\boldsymbol{0}} $$ (13) 2. 三维20节点六面体和10节点四面体单元的高精度集中质量矩阵构造方法
为表述简洁起见, 在下面的讨论中, 分别用H20和T10表示三维20节点六面体和10节点四面体单元.
2.1 H20和T10单元的HRZ集中质量矩阵
由于满足变分一致性, 式(9)定义的质量矩阵通常被称为一致质量矩阵. 当采用行求和法构造集中质量矩阵时, 仅需将一致质量矩阵每行的元素累加到主对角元素. 因此, 对于标量波动问题, 考虑规则的H20和T10单元构形, 例如长方体和直边四面体, 其行求和集中质量矩阵分别为
$$ {\boldsymbol{M}}_{{\rm{RS}}{{ \text{-} }}{\rm{H}}20}^e = \frac{{{m^e}}}{{24}}{\rm{diag}}\{ \underbrace { - 3, \cdots , - 3}_8,\underbrace {4, \cdots ,4}_{12}\} $$ (14) $$ {\boldsymbol{M}}_{{\rm{RS}}{{ \text{-} }}{\rm{T}}10}^e = \frac{{{m^e}}}{{20}}{\rm{diag}}\{ \underbrace { - 1, \cdots , - 1}_4,\underbrace {4, \cdots ,4}_6\} $$ (15) 其中, 下标${\rm{RS}}$表示行求和法; $ {m^e} $为单元的总质量. 为清晰起见, 图1和图2给出了H20和T10单元的节点质量分布, 其中MNLM (mid-node lumped mass matrix)表示中节点集中质量矩阵.
![]()
图 1 H20单元3种集中质量矩阵的节点质量分布示意图Figure 1. Schematic illustration of the nodal mass distribution for three mass lumping schemes of H20 element![]()
图 2 T10单元3种集中质量矩阵的节点质量分布示意图Figure 2. Schematic illustration of the nodal mass distribution for three mass lumping schemes of T10 element由式(14)和式(15)可得, 对于H20和T10单元, 行求和集中质量矩阵中的角节点质量均为负值. 根据HRZ方法, 将式(9)中的一致质量矩阵的主对角元素在满足质量守恒的条件下进行比例缩放, 即可得到非负的HRZ集中质量矩阵
$$ {({\boldsymbol{M}}_{{\rm{HRZ}}}^e)_{II}} = r\int_{{\varOmega ^e}} {N_I^2{\rm{d}}} \varOmega $$ (16) 其中, $ r $为缩放系数
$$ r = {m^e}\Bigg/\int_{{\varOmega ^e}} {\sum\limits_{I = 1}^{{n_{en}}} {N_I^2} {\rm{d}}} \varOmega $$ (17) 由式(16)和式(17), 可得H20和T10单元的HRZ集中质量矩阵
$$ {\boldsymbol{M}}_{{\rm{HRZ}}{{ \text{-} }}{\rm{H}}20}^e = \frac{{{m^e}}}{{248}}{\rm{diag}}\{ \underbrace {7, \cdots ,7}_8,\underbrace {16, \cdots ,16}_{12}\} $$ (18) $$ {\boldsymbol{M}}_{{\rm{HRZ}}{{ \text{-} }}{\rm{T}}10}^e = \frac{{{m^e}}}{{108}}{\rm{diag}}\{ \underbrace {3, \cdots ,3}_4,\underbrace {16, \cdots ,16}_6\} $$ (19) 图1和图2也给出了这两种单元的HRZ集中质量矩阵节点质量分布情况.
2.2 H20和T10单元的非负中节点集中质量矩阵
如前所述, H20和T10单元的行求和集中质量矩阵, 负质量都出现在角节点上, 而中节点则没有负质量问题. 因此, 可以基于行求和集中质量矩阵中的非负质量元素, 通过单元质量守恒将其进行比例放缩, 同时将行求和集中质量矩阵中的负对角元素置零, 构造一种新型非负集中质量矩阵. 该集中质量矩阵${\boldsymbol{M}}_{{\rm{MNLM}}}^e$的构造流程为
$$ {({\boldsymbol{M}}_{{\rm{MNLM}}}^e)_{II}} = 0,\;\;{\rm{if}}\;{({\boldsymbol{M}}_{{\rm{RS}}}^e)_{II}} < 0 $$ (20) $$ {({\boldsymbol{M}}_{{\rm{MNLM}}}^e)_{II}} = \frac{{{{({\boldsymbol{M}}_{{\rm{RS}}}^e)}_{II}}}}{{\displaystyle\sum\limits_{{\rm{if}}\;{{({\boldsymbol{M}}_{{\rm{RS}}}^e)}_{II}} > 0} {{{({\boldsymbol{M}}_{{\rm{RS}}}^e)}_{II}}} }},\;\;{\rm{if}}\;{({\boldsymbol{M}}_{{\rm{RS}}}^e)_{II}} \geqslant 0 $$ (21) 基于式(20)和式(21), H20和T10单元的新型集中质量矩阵分别为
$$ {\boldsymbol{M}}_{{\rm{MNLM}} \text{-} {\rm{H}}20}^e = \frac{{{m^e}}}{{12}}{\rm{diag}}\{ \underbrace {0, \cdots ,0}_8,\underbrace {1, \cdots ,1}_{12}\} $$ (22) $$ {\boldsymbol{M}}_{{\rm{MNLM}} \text{-} {\rm{T}}10}^e = \frac{{{m^e}}}{6}{\rm{diag}}\{ \underbrace {0, \cdots ,0}_4,\underbrace {1, \cdots ,1}_6\} $$ (23) 由式(22)和式(23)可见, 对于H20和T10单元, 角节点的质量均为零, 质量只被分配在中节点上, 所以把该集中质量矩阵称之为中节点集中质量矩阵, 简称MNLM集中质量矩阵. 同样, 图1和图2中也描述了H20和T10单元的MNLM集中质量矩阵.
3. H20单元集中质量矩阵的频率精度分析
本节以H20单元为例, 研究集中质量矩阵构造形式对频率计算精度的影响. 为了更全面衡量各类集中质量矩阵的精度, 构造如下H20单元广义集中质量矩阵
$$ \left. \begin{split} & {\boldsymbol{M}}_{}^e = \frac{{{m^e}}}{{\text{8}}}{\rm{diag}}\{ \underbrace {\alpha , \cdots ,\alpha }_8,\underbrace {\beta , \cdots ,\beta }_{12}\} \\ & \beta \; = 2(1 - \alpha )/3\end{split} \right\} $$ (24) 其中, $ \alpha $和$ \beta $是待定参数. 为了保证集中质量矩阵的非负性, 角节点质量$ \alpha $需满足$ \alpha \in [0,1] $. 当$ \alpha = 7/31 $时, 式(24)中的广义集中质量矩阵即退化为式(18)的HRZ集中质量矩阵; 而当$ \alpha = 0 $时, 即为式(22)的MNLM集中质量矩阵.
为了进行频率精度分析, 考虑图3所示的H20单元典型节点影响域. 由H20单元构成的有限元典型节点包括一个角节点$ {{\boldsymbol{x}}_{abc}} $和3个中间节点$ {{\boldsymbol{x}}_{a(b + 1)c}} $, $ {{\boldsymbol{x}}_{(a - 1)bc}} $和$ {{\boldsymbol{x}}_{ab(c + 1)}} $. $ \theta $和$ \varphi $是与谐波相关的不同方向的入射角, $ {h_i} $为$ {x_i} $方向的单元长度, 不同方向的波数$ {k_i} $与$ k $之间的关系满足
![]()
图 3 H20单元典型节点影响域及波动方向示意图Figure 3. Illustration of the typical nodal influence domains and wave propagation angles for a mesh with H20 elements$$ \left.\begin{split} & {k_x} = k\cos \varphi \cos \theta \; \\ & {k_y} = k\cos \varphi \sin \theta \; \\ & {k_z} = k\sin \varphi \end{split} \right\} $$ (25) 假设离散节点的波动形式为谐波, 4个典型节点系数可以表示为
$$ \left. \begin{split} & {d_{abc}}(t) = {{\hat u}_1}\exp [\text{ι} ({k_x}{x_a} + {k_y}{y_b} + {k_z}{z_c} - {\omega ^h}t)] \\ & {d_{a(b + 1)c}}(t) = {{\hat u}_2}\exp [\text{ι} ({k_x}{x_a} + {k_y}{y_{b + 1}} + {k_z}{z_c} - {\omega ^h}t)] \\ & {d_{(a - 1)bc}}(t) = {{\hat u}_3}\exp [\text{ι} ({k_x}{x_{a - 1}} + {k_y}{y_b} + {k_z}{z_c} - {\omega ^h}t)] \\ & {d_{ab(c + 1)}}(t) = {{\hat u}_4}\exp [\text{ι} ({k_x}{x_a} + {k_y}{y_b} + {k_z}{z_{c + 1}} - {\omega ^h}t)] \end{split} \right\} $$ (26) 其中, $ {\hat u_I} $表示第$ I $个典型节点的波动振幅. 根据图3, 可得各个离散节点之间的位置关系为
$$ \left. \begin{split} & {x_{a + m}} = {x_a} + m {h_x}/2 \\ & {y_{b + n}} = {y_b} + n {h_y}/2 \\ & {z_{c + q}} = {z_c} + q {h_z}/2\end{split} \right\} $$ (27) 方便表述起见, 引入算子$ {\cos _{( \pm l{k_x} \pm m{k_y} \pm n{k_z})}} $
$$ \begin{split} & {\cos _{( \pm l{k_x} \pm m{k_y} \pm n{k_z})}} = \sum\limits_{I = 1}^2 \sum\limits_{J = 1}^2 \sum\limits_{K = 1}^2 \left\{ \cos [{( - 1)^I}l{k_x}{h_x}/2 +\right. \\ &\qquad \left. {( - 1)^J}m{k_y}{h_y}/2 + {( - 1)^K}n{k_z}{h_z}/2] \right\} \end{split} $$ (28) 当式(28)仅包括两个下标时, 有
$$ \begin{split} & {\cos _{( \pm l{k_x} \pm m{k_y})}} = \\ &\qquad \sum\limits_{I = 1}^2 {\sum\limits_{J = 1}^2 {\left\{ \cos [{( - 1)^I}l{k_x}{h_x}/2 + {( - 1)^J}m{k_y}{h_y}/2] \right\}} } \end{split} $$ (29) 对于4个典型节点, 将式(26)和式(27)代入式(7), 并利用简谐波假定化简可得
$$ {{{\boldsymbol{A}}\hat {\boldsymbol{u}}}} = \left[ {\begin{array}{*{20}{c}} {{{\bar A}_{11}}}&{{A_{12}}}&{{A_{13}}}&{{A_{14}}} \\ {{A_{21}}}&{{{\bar A}_{22}}}&{{A_{23}}}&{{A_{24}}} \\ {{A_{31}}}&{{A_{32}}}&{{{\bar A}_{33}}}&{{A_{34}}} \\ {{A_{41}}}&{{A_{42}}}&{{A_{43}}}&{{{\bar A}_{44}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\hat u}_1}} \\ {{{\hat u}_2}} \\ {{{\hat u}_3}} \\ {{{\hat u}_4}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} 0 \\ 0 \\ 0 \\ 0 \end{array}} \right\} $$ (30) 其中
$$ \left. \begin{split} & {{\bar A}_{11}} = - {[{\omega ^h}/(kc)]^2}\tilde A + {A_{11}} \\ & {{\bar A}_{22}} = - {[{\omega ^h}/(kc)]^2}{\overset{\frown} A} + {A_{22}} \\ & {{\bar A}_{33}} = - {[{\omega ^h}/(kc)]^2}{\overset{\frown} A} + {A_{33}} \\ & {{\bar A}_{44}} = - {[{\omega ^h}/(kc)]^2}{\overset{\frown} A} + {A_{44}} \\ & \tilde A = {k^2}\alpha {m^e},\;{\overset{\frown} A} = {k^2}{m^e}(1 - \alpha )/3\; \end{split} \right\} $$ (31) $$\begin{split} & {A_{11}} = {K_{2,8}}{\cos _{( \pm 2{k_x} \pm 2{k_y} \pm 2{k_z})}} + {K_{2,4}}{\cos _{( \pm 2{k_x} \pm 2{k_y} \pm 0{k_z})}} +\\ &\qquad {K_{2,5}}{\cos _{( \pm 2{k_x} \pm 0{k_y} \pm 2{k_z})}} + {K_{2,1}}{\cos _{( \pm 2{k_x} \pm 0{k_y} \pm 0{k_z})}}+ \\ &\qquad {K_{2,7}}{\cos _{( \pm 0{k_x} \pm 2{k_y} \pm 2{k_z})}} + {K_{2,3}}{\cos _{( \pm 0{k_x} \pm 2{k_y} \pm 0{k_z})}}+ \\ &\qquad {K_{2,6}}{\cos _{( \pm 0{k_x} \pm 0{k_y} \pm 2{k_z})}} + {K_{2,2}}{\cos _{( \pm 0{k_x} \pm 0{k_y} \pm 0{k_z})}} \end{split} $$ (32) $$ \begin{split} & {A_{22}} = {K_{10,20}}{\cos _{( \pm 2{k_x} \pm 2{k_z})}} + {K_{10,12}}{\cos _{( \pm 2{k_x} \pm 0{k_z})}}+ \\ & \qquad {K_{10,18}}{\cos _{( \pm 0{k_x} \pm 2{k_z})}} + {K_{10,10}}{\cos _{( \pm 0{k_x} \pm 0{k_z})}} \end{split} $$ (33) $$ \begin{split} & {A_{33}} = {K_{9,19}}{\cos _{( \pm 2{k_y} \pm 2{k_z})}} + {K_{9,11}}{\cos _{( \pm 2{k_y} \pm 0{k_z})}} +\\ &\qquad {K_{9,17}}{\cos _{( \pm 0{k_y} \pm 2{k_z})}} + {K_{9,9}}{\cos _{( \pm 0{k_y} \pm 0{k_z})}} \end{split} $$ (34) $$ \begin{split} & {A_{44}} = {K_{14,16}}{\cos _{( \pm 2{k_x} \pm 2{k_y})}} + {K_{14,13}}{\cos _{( \pm 2{k_x} \pm 0{k_y})}}+ \\ &\qquad {K_{14,15}}{\cos _{( \pm 0{k_x} \pm 2{k_y})}} + {K_{14,14}}{\cos _{( \pm 0{k_x} \pm 0{k_y})}} \end{split} $$ (35) $$\begin{split} & {A_{12}} = {K_{2,20}}{\cos _{( \pm 2{k_x} \pm {k_y} \pm 2{k_z})}} + {K_{2,12}}{\cos _{( \pm 2{k_x} \pm {k_y} \pm 0{k_z})}}+ \\ &\qquad {K_{2,18}}{\cos _{( \pm 0{k_x} \pm {k_y} \pm 2{k_z})}} + {K_{2,10}}{\cos _{( \pm 0{k_x} \pm {k_y} \pm 0{k_z})}} \end{split} $$ (36) $$ \begin{split} & {A_{13}} = {K_{2,19}}{\cos _{( \pm {k_x} \pm 2{k_y} \pm 2{k_z})}} + {K_{2,11}}{\cos _{( \pm {k_x} \pm 2{k_y} \pm 0{k_z})}}+ \\ &\qquad {K_{2,17}}{\cos _{( \pm {k_x} \pm 0{k_y} \pm 2{k_z})}} + {K_{2,09}}{\cos _{( \pm {k_x} \pm 0{k_y} \pm 0{k_z})}}\end{split} $$ (37) $$ \begin{split} & {A_{14}} = {K_{2,16}}{\cos _{( \pm 2{k_x} \pm 2{k_y} \pm {k_z})}} + {K_{2,13}}{\cos _{( \pm 2{k_x} \pm 0{k_y} \pm {k_z})}} +\\ &\qquad {K_{2,15}}{\cos _{( \pm 0{k_x} \pm 2{k_y} \pm 1{k_z})}} + {K_{2,14}}{\cos _{( \pm 0{k_x} \pm 0{k_y} \pm {k_z})}} \end{split} $$ (38) $$ {A_{23}} = {K_{10,19}}{\cos _{( \pm {k_x} \pm {k_y} \pm 2{k_z})}} + {K_{10,11}}{\cos _{( \pm {k_x} \pm {k_y} \pm 0{k_z})}} $$ (39) $$ {A_{24}} = {K_{10,16}}{\cos _{( \pm 2{k_x} \pm {k_y} \pm {k_z})}} + {K_{10,15}}{\cos _{( \pm 0{k_x} \pm {k_y} \pm {k_z})}} $$ (40) $$ {A_{34}} = {K_{9,16}}{\cos _{( \pm {k_x} \pm 2{k_y} \pm {k_z})}} + {K_{9,13}}{\cos _{( \pm {k_x} \pm 0{k_y} \pm {k_z})}} $$ (41) 其中, $ {A_{JI}} = {A_{IJ}} $, $ {K_{I,J}} = {K_{IJ}} $为刚度矩阵元素, 为了不造成下标混淆, 下标的行号和列号之间采用“, ”隔开.
由于式(30)有非0解, 则矩阵A的行列式应为0
$$ \det ({\boldsymbol{A}}) = 0 $$ (42) 值得指出的是, 式(42)是关于数值频率$ {\omega ^h} $的非线性方程, 但由于其复杂性, 从中直接求解$ {\omega ^h} $通常是非常困难的. 另一方面, 在理论分析中我们关注的是频率误差, 而非频率本身, 因而可直接引入频率误差$ e $[1]
$$ e = ({\omega ^h} - \omega )/\omega $$ (43) 结合解析频率$ \omega = kc $和式(43), 有
$$ {({\omega ^h})^2} = {(1 + e)^2}{\omega ^2} = {(1 + e)^2}(k_x^2 + k_y^2 + k_z^2){c^2} $$ (44) 将式(44)代入式(42), 并忽略$ e $的高次项[26], 可得
$$ {a_1}e + {a_0} \approx 0 $$ (45) $$ \begin{split} & {a_1} = ( - 6{\kern 1pt} {A_{11}} + 8{\kern 1pt} \tilde A){{\overset{\frown} A^3}} + [ - 6({A_{22}}{\kern 1pt} + {A_{33}} + {\kern 1pt} {A_{44}})\tilde A +\\ &\qquad 4({A_{22}} + {A_{33}} + {A_{44}}){A_{11}} - 4{\kern 1pt} (A_{12}^2 + A_{13}^2 + A_{14}^2)]{{{\overset{\frown} A}^2 }}+ \\ &\qquad \{ [4{\kern 1pt} ({A_{33}} + {A_{44}}){A_{22}} - 4{\kern 1pt} ({\kern 1pt} A_{23}^2 + {\kern 1pt} A_{24}^2 + {\kern 1pt} A_{34}^2)+ \\ &\qquad 4{A_{33}}{\kern 1pt} {A_{44}}]\tilde A - 2{\kern 1pt} [({A_{33}} + {A_{44}}){A_{22}} + {A_{33}}{\kern 1pt} {A_{44}}- \\ & \qquad (A_{23}^2 + A_{24}^2 + A_{34}^2)]{A_{11}} + 2[(A_{13}^2 + {\kern 1pt} A_{14}^2){A_{22}}+ \\ &\qquad ({\kern 1pt} A_{12}^2 + {\kern 1pt} A_{14}^2){A_{33}} + (A_{12}^2 + A_{13}^2){A_{44}} + 2{A_{13}}{A_{14}}{A_{34}}- \\ &\qquad 2({A_{13}}{A_{23}}{\kern 1pt} + {\kern 1pt} {A_{14}}{A_{24}}{\kern 1pt} ){A_{12}}]\} {\overset{\frown} A} - 2[{\kern 1pt} 2{A_{23}}{\kern 1pt} {A_{34}}{\kern 1pt} {A_{24}}+ \\ &\qquad ({A_{33}}{\kern 1pt} {A_{44}} - A_{34}^2){A_{22}} - A_{23}^2{\kern 1pt} {A_{44}} - A_{24}^2{\kern 1pt} {A_{33}}{\kern 1pt} ]\tilde A \end{split} $$ (46) $$\begin{split} & {a_0} = ( - {A_{11}} + {\kern 1pt} \tilde A){{\overset{\frown} A }^3} + [({A_{22}} + {A_{33}} + {A_{44}}){A_{11}}- \\ &\qquad ({A_{22}} + {A_{33}} + {A_{44}}){\kern 1pt} \tilde A - (A_{12}^2 + A_{13}^2 + A_{14}^2)]{{\overset{\frown} A }^2}+ \\ &\qquad \{ [ - ({A_{33}} + {A_{44}}){A_{22}} - {A_{33}}{\kern 1pt} {A_{44}} + A_{23}^2 + {\kern 1pt} A_{24}^2+ \\ &\qquad {\kern 1pt} A_{34}^2]{A_{11}} + (A_{13}^2 + A_{14}^2 + {A_{33}}{\kern 1pt} \tilde A + {A_{44}}{\kern 1pt} \tilde A){A_{22}}+ \\ &\qquad (A_{12}^2 + A_{14}^2 + {A_{44}}{\kern 1pt} \tilde A){A_{33}} + (A_{12}^2 + A_{13}^2){A_{44}}- \\ &\qquad (A_{23}^2 + A_{24}^2 + A_{34}^2)\tilde A - 2({A_{13}}{A_{23}} + {A_{14}}{\kern 1pt} {A_{24}}){A_{12}}- \\ &\qquad 2{A_{13}}{A_{14}}{A_{34}}\} \overset{\frown} A + [({A_{33}}{A_{44}}{\mkern 1mu} - A_{34}^2){A_{22}} - {A_{44}}{\mkern 1mu} A_{23}^2+ \\ &\qquad 2{\mkern 1mu} {A_{23}}{\mkern 1mu} {A_{24}}{\mkern 1mu} {A_{34}} - A_{24}^2{A_{33}}]{A_{11}} - [(A_{14}^2 + \tilde A{\mkern 1mu} {A_{44}}){A_{33}}+ \\ &\qquad {A_{44}}{\mkern 1mu} A_{13}^2 - \tilde A{\mkern 1mu} A_{34}^2 - 2{\mkern 1mu} {A_{13}}{\mkern 1mu} {A_{14}}{\mkern 1mu} {A_{34}}]{A_{22}}+ \\ &\qquad ( - A_{12}^2{A_{44}} + 2{\mkern 1mu} {A_{14}}{\mkern 1mu} {A_{24}}{\mkern 1mu} {A_{12}} + A_{24}^2\tilde A){A_{33}}+ \\ &\qquad (2{\mkern 1mu} {A_{12}}{A_{13}}{\mkern 1mu} {A_{23}}{\mkern 1mu} + A_{23}^2\tilde A){A_{44}} + A_{12}^2A_{34}^2- \\ &\qquad 2{\mkern 1mu} {A_{23}}{\mkern 1mu} {A_{24{\mkern 1mu} }}{A_{34{\mkern 1mu} }}\tilde A - 2{\mkern 1mu} {A_{12}}{A_{13{\mkern 1mu} }}{A_{24}}{A_{34}}- \\ &\qquad 2{\mkern 1mu} {\mkern 1mu} {A_{12}}{A_{23}}{A_{14}}{A_{34}} + {({A_{13{\mkern 1mu} }}{A_{24}} - {A_{14}}{\mkern 1mu} {A_{23}})^2} \end{split} $$ (47) 进一步在式(45)中引入余弦项关于$ kh $的泰勒展开, 可得H20单元广义集中质量矩阵的频率误差表达式
$$ e \approx - \frac{{{a_0}}}{{{a_1}}} \approx - \frac{{{{(1 + \alpha )}^2}}}{{96}}({k^2}{h^2}) + \mathcal{O}({k^4}{h^4}) $$ (48) 其中$ h $为特征单元长度
$$ h = \sqrt {h_x^2 + h_y^2 + h_z^2} $$ (49) 将$ \alpha = 7/31 $代入式(48), 即得H20单元HRZ集中质量矩阵的频率误差表达式
$$ {e_{{\rm{HRZ}}}} \approx - \frac{{361}}{{23\;064}}({k^2}{h^2}) + \mathcal{O}({k^4}{h^4}) $$ (50) 对比式(48)和式(50), 可见H20单元的HRZ集中质量矩阵的精度并非最优. 反之, 由式(48)可知, $ \alpha = 0 $对应的H20单元集中质量矩阵具有更优的频率计算精度, 此时的集中质量矩阵即为MNLM集中质量矩阵. 进一步, 将$ \alpha = 0 $代入式(48)中, 即得MNLM集中质量矩阵的频率误差表达式
$$ {e_{{\rm{MNLM}}}} \approx - \frac{1}{{96}}({k^2}{h^2}) + \mathcal{O}({k^4}{h^4}) $$ (51) 基于式(50)和式(51), 若忽略频率误差的高阶项, 可得三维H20单元MNLM集中质量矩阵与HRZ集中质量矩阵两者频率误差之间存在如下关系
$$ \frac{{{e_{{\rm{MNLM}}}}}}{{{e_{{\rm{HRZ}}}}}} \approx \frac{1}{{96}} \times \frac{{23\;064}}{{361}} = \frac{{{\text{961}}}}{{{\text{1444}}}} \approx 0.67 $$ (52) 4. 基于中节点集中质量矩阵的时程分析
对于H20和T10单元, 式(22)和式(23)表明MNLM集中质量矩阵的角节点质量为0. 利用这一特点, 可通过静力凝聚建立相应的结构动力分析降阶模型[27], 减小整体计算规模, 提高计算效率.
将H20或T10单元的MNLM集中质量矩阵代入式(7)的有限元离散运动方程, 同时将零质量角节点和非零质量中节点进行分类, 可得
$$ \begin{split} &\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{0}}_{CC}}}&{{{\boldsymbol{0}}_{CM}}} \\ {{{\boldsymbol{0}}_{MC}}}&{{{\boldsymbol{M}}_{MM}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{a}}_C}} \\ {{{\boldsymbol{a}}_M}} \end{array}} \right\} +\\ &\qquad \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_{CC}}}&{{{\boldsymbol{K}}_{CM}}} \\ {{{\boldsymbol{K}}_{MC}}}&{{{\boldsymbol{K}}_{MM}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{d}}_C}} \\ {{{\boldsymbol{d}}_M}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{f}}_C}} \\ {{{\boldsymbol{f}}_M}} \end{array}} \right\} \end{split}$$ (53) 其中, 下标$ C $和$ M $分别表示单元的角节点和中节点, $ {{\boldsymbol{M}}_{MM}} $, $ {{\boldsymbol{K}}_{CC}} $和$ {{\boldsymbol{K}}_{MM}} $均为正定矩阵. 对式(53)中的零质量节点对应的元素, 可通过静力凝聚进行模型降阶. 将式(53)展开, 有
$$ {{\boldsymbol{K}}_{CC}}{{\boldsymbol{d}}_C} + {{\boldsymbol{K}}_{CM}}{{\boldsymbol{d}}_M} = {{\boldsymbol{f}}_C} $$ (54) $$ {{\boldsymbol{M}}_{MM}}{{\boldsymbol{a}}_M} + {{\boldsymbol{K}}_{MC}}{{\boldsymbol{d}}_C} + {{\boldsymbol{K}}_{MM}}{{\boldsymbol{d}}_M} = {{\boldsymbol{f}}_M} $$ (55) 根据式(54), 角节点位移向量$ {{\boldsymbol{d}}_C} $可表示为
$$ {{\boldsymbol{d}}_C} = {\boldsymbol{K}}_{CC}^{ - 1}({{\boldsymbol{f}}_C} - {{\boldsymbol{K}}_{CM}}{{\boldsymbol{d}}_M}) $$ (56) 将式(56)代入式(55)中, 可得仅与中节点有关的有限元离散运动方程
$$ {{\boldsymbol{M}}_{MM}}{{\boldsymbol{a}}_M} + {{\boldsymbol{\bar K}}_{MM}}{{\boldsymbol{d}}_M} = {{\boldsymbol{\bar f}}_M} $$ (57) 其中
$$ {{\boldsymbol{\bar K}}_{MM}} = {{\boldsymbol{K}}_{MM}} - {{\boldsymbol{K}}_{MC}}{\boldsymbol{K}}_{CC}^{ - 1}{{\boldsymbol{K}}_{CM}} $$ (58) $$ {{\boldsymbol{\bar f}}_M} = {{\boldsymbol{f}}_M} - {{\boldsymbol{K}}_{MC}}{\boldsymbol{K}}_{CC}^{ - 1}{{\boldsymbol{f}}_C} $$ (59) 对于结构时程分析, 这里采用中心差分法进行时间离散. 当采用等时间步长$ \Delta t $时, 中心差分法第$ n $步的速度和加速度分别为
$$ {\boldsymbol{v}}_M^n = \frac{1}{{2\Delta t}}\left({\boldsymbol{d}}_M^{n + 1} - {\boldsymbol{d}}_M^{n - 1}\right) $$ (60) $$ {\boldsymbol{a}}_M^n = \frac{1}{{{{(\Delta t)}^2}}}\left({\boldsymbol{d}}_M^{n + 1} + {\boldsymbol{d}}_M^{n - 1} - 2{\boldsymbol{d}}_M^n\right) $$ (61) 将式(61)代入式(57), 可得第$ n + 1 $步的位移更新格式
$$ {\boldsymbol{d}}_M^{n + 1} = {(\Delta t)^2}{\boldsymbol{M}}_{MM}^{ - 1}({\boldsymbol{\bar f}}_M^n - {{\boldsymbol{\bar K}}_{MM}}{\boldsymbol{d}}_M^n) + 2{\boldsymbol{d}}_M^n - {\boldsymbol{d}}_M^{n - 1} $$ (62) 与此同时, 将式(60)代入式(61), 可以得到$ n = - 1 $步的起步位移向量
$$ {\boldsymbol{d}}_M^{ - 1} = {\boldsymbol{d}}_M^0 - \Delta t{\boldsymbol{v}}_M^0 + \frac{{{{(\Delta t)}^2}}}{2}{\boldsymbol{a}}_M^0 $$ (63) 其中, $ {\boldsymbol{d}}_M^0 $, $ {\boldsymbol{v}}_M^0 $和$ {\boldsymbol{a}}_M^0 $分别为初始时刻中节点的位移、速度和加速度. 在时程分析中, 通过式(62)可计算任意时间步中节点的位移, 进一步利用式(56)可得到相应的角节点位移.
5. 数值算例
由于本文主要研究不同方法的精度对比, 简洁起见且不失一般性, 算例均采用无量纲单位.
5.1 长方体空腔问题
首先, 考虑齐次边界条件的长方体空腔频率计算问题, 长方体空腔的长度、宽度和高度分别为$ {L_x} = 2 $, $ {L_y} = 1 $和$ {L_z} = 1 $. 该问题的频率解析解[28]为
$$ {\omega _{lmn}} = c\sqrt {{{(l\text{π} /{L_x})}^2} + {{(m\text{π} /{L_y})}^2} + {{(n\text{π} /{L_z})}^2}} $$ (64) 其中, $ c $为波速, 本算例中取为1. 图4为该问题采用H20和T10两种单元的网格划分情况, 其中, H20单元模型包括216, 512和1000个单元(1225, 2673和4961个自由度), T10单元模型包括1080, 2560和5000个单元(1981, 4401和8261个自由度). 为了方便对比, 图4有限元模型中, 一个H20单元对应5个T10单元.
图5和表1给出了采用图4系列有限元模型得到的前4阶频率计算结果. 从图5的频率收敛结果可见, H20和T10单元的MNLM集中质量矩阵的频率计算精度均明显优于HRZ集中质量矩阵. 同时, 表1的频率误差结果对比表明, 对于H20单元, MNLM集中质量矩阵的前4阶频率计算误差与HRZ之间比值约为0.66, 与式(52)给出的理论结果吻合良好; 对于T10单元, MNLM集中质量矩阵与HRZ集中质量矩阵相比的频率误差比值小于0.40, 精度优势更为显著. 另外, 值得指出的是, 无论是HRZ还是MLNM集中质量矩阵, 在相同自由度数条件下, T10单元的频率计算精度均优于H20单元.
![]()
图 5 长方体空腔问题前四阶频率收敛特性对比Figure 5. Convergence comparison of the first four frequencies for the rectangular cavity problem表 1 长方体空腔问题前四阶频率计算精度对比Table 1. Accuracy comparison of the first four frequencies for the rectangular cavity problem$ {\omega _i} $ $ {n_{el}}( \times 5) $ $\omega _i^h{ { \text{-} } }{\rm{H} }20$ $\left| {\omega _i^h{{ - } }\omega _i^{} } \right|$ $\left| {\dfrac{ { {e_{{\rm{MNLM}}} } } }{ { {e_{{\rm{HRZ}}} } } } } \right|$ $\omega _i^h{ { \text{-} } }{\rm{T} }10$ $\left| {\omega _i^h{{ - } }\omega _i^{} } \right|$ $\left| {\dfrac{ { {e_{{\rm{MNLM}}} } } }{ { {e_{{\rm{HRZ}}} } } } } \right|$ HRZ MNLM HRZ MNLM HRZ MNLM HRZ MNLM 4.71 216 4.45 4.54 0.26 0.17 0.66 4.64 4.68 0.08 0.03 0.38 512 4.56 4.61 0.15 0.10 0.66 4.67 4.69 0.04 0.02 0.39 1000 4.62 4.65 0.10 0.06 0.66 4.68 4.70 0.03 0.01 0.40 5.44 216 5.04 5.17 0.41 0.27 0.66 5.32 5.40 0.12 0.04 0.36 512 5.21 5.29 0.23 0.15 0.66 5.37 5.42 0.07 0.03 0.38 1000 5.29 5.34 0.15 0.10 0.66 5.40 5.42 0.04 0.02 0.39 6.48 216 5.77 6.01 0.70 0.46 0.66 6.28 6.41 0.20 0.07 0.35 512 6.08 6.21 0.40 0.26 0.66 6.36 6.43 0.11 0.04 0.37 1000 6.22 6.31 0.25 0.17 0.66 6.40 6.45 0.07 0.03 0.39 7.20 216 6.32 6.62 0.88 0.58 0.66 6.94 7.11 0.26 0.09 0.34 512 6.68 6.86 0.52 0.34 0.66 7.05 7.14 0.15 0.06 0.37 1000 6.86 6.97 0.34 0.22 0.66 7.10 7.16 0.10 0.04 0.39 5.2 弹性力学长方体问题
第2个算例是弹性力学长方体问题. 该模型的几何尺寸与长方体空腔问题一致, 材料弹性拉梅常数$ \lambda = 15/26 $, $ \mu = 5/13 $. 边界条件如图6所示: 在$ x = \{ 0,\;{L_x}\} $处$ {u_y} = {u_z} = 0 $, 在$ y = \{ 0,\;{L_y}\} $处$ {u_x} = {u_z} = 0 $, 在$ z = \{ 0,\;{L_z}\} $处$ {u_x} = {u_y} = 0 $. 根据文献[29], 可导出该模型的压力波和剪力波频率解析解分别为
$$ \omega _{lmn}^p = k{c_p} = \sqrt {{{\left(\frac{{l\text{π} }}{{{L_x}}}\right)}^2} + {{\left(\frac{{m\text{π} }}{{{L_y}}}\right)}^2} + {{\left(\frac{{n\text{π} }}{{{L_z}}}\right)}^2}} \sqrt {\frac{{\lambda + 2\mu }}{\rho }} $$ (65) $$ \omega _{lmn}^s = k{c_s} = \sqrt {{{\left(\frac{{l\text{π} }}{{{L_x}}}\right)}^2} + {{\left(\frac{{m\text{π} }}{{{L_y}}}\right)}^2} + {{\left(\frac{{n\text{π} }}{{{L_z}}}\right)}^2}} \sqrt {\frac{\mu }{\rho }} $$ (66) 其中, 对于压力波频率, 下标需满足$ l,\;m,\;n \geqslant 1 $; 对于剪力波频率, 下标需满足至少有两个大于等于1. 值得需要注意的是, 下标均大于等于1时, 剪力波频率为二重根.
在对该结构进行自由振动分析时, 同样采用图4的有限元网格进行计算. 前6阶频率的收敛对比结果见图7. 可见, H20和T10单元的MNLM集中质量矩阵的频率计算精度均明显优于HRZ, 验证了所提方法在弹性力学问题中同样具有良好的适用性.
![]()
图 7 长方体弹性体前六阶频率收敛特性对比Figure 7. Convergence comparison of the first six frequencies for the rectangular elastic continuum problem为了深入验证所提集中质量矩阵的动力性能, 进一步计算该弹性力学问题的时程响应. 计算中考虑如下的位移解析解
$$ \left.\begin{split} & {u_x}({\boldsymbol{x}},t) = \sin (\text{π} y)\sin (\text{π} z)\sin \left(\sqrt {2\mu } \text{π} t\right) \\ & {u_y}({\boldsymbol{x}},t) = \sin (1.5\text{π} x)\sin (\text{π} z)\sin \left(\sqrt {3.25\mu } \text{π} t\right) \\ & {u_z}({\boldsymbol{x}},t) = \sin (\text{π} y)\sin (\text{π} z)\sin \left(\sqrt {1.25\mu } \text{π} t\right) \end{split} \right\} $$ (67) 其中, 取$ t = 0 $即可得到有限元动力分析的初始位移和初始速度.
有限元分析中H20和T10单元分别采用图4中的1000个单元和5000个单元进行时程分析, 时间步长统一取$ \Delta t = 0.01 $. 图8给出了$ t = 10 $时刻的各方向位移误差云图, 结果表明H20和T10单元的MNLM集中质量矩阵的位移计算精度均优于对应的HRZ集中质量矩阵. 此外, 图9给出了归一化计算时间[30]对比结果. 相比HRZ集中质量矩阵, H20单元的MNLM集中质量矩阵的计算效率提升了约4.3倍, T10单元的MNLM集中质量矩阵的计算效率提升了约3.0倍. 由此可得, MNLM集中质量矩阵不仅可以显著提高计算精度, 而且在进行时程分析时, 与结构动力分析降阶模型相结合可以提高计算效率.
![]()
图 8 长方体弹性体$ t = 10 $时刻的位移误差云图Figure 8. Displacement error contour plots for the rectangular elastic continuum problem at $ t = 10 $![]()
图 9 长方体弹性体时程分析效率对比Figure 9. Efficiency comparison of the transient analysis for the rectangular elastic continuum problem5.3 1/4圆筒空腔问题
最后一个算例是图10所示的1/4圆筒空腔问题. 该模型的内径${r_{\rm{i}}} = 1$, 外径${r_{\rm{o}}} = 2$, 高度$ H = 1 $, 波速$ c = 1 $. 当采用齐次边界条件时, 该问题的频率解析解为
$$ {\omega _{lmn}} = \sqrt {{{(a_l^m)}^2} + {{(n\text{π} /H)}^2}} c $$ (68) 其中, $ a_l^m $为下面Bessel方程的根
$$ {\rm{J}}_l^{[1]}\left(\frac{{a_l^m}}{c}{r_{\rm{i}}}\right){\rm{J}}_l^{[2]}\left(\frac{{a_l^m}}{c}{r_{\rm{o}}}\right) - {\rm{J}}_l^{[2]}\left(\frac{{a_l^m}}{c}{r_{\rm{i}}}\right){\rm{J}}_l^{[1]}\left(\frac{{a_l^m}}{c}{r_{\rm{o}}}\right) = 0 $$ (69) 式中, ${\rm{J}}_l^{[1]}$和${\rm{J}}_l^{[2]}$分别为第1类和第2类$ l $阶Bessel函数. 需要注意的是, 式(68)和式(69)中下标$ l $为正偶数. 对于1/4圆筒空腔模型, 由于受到自身几何形状的限制, 只能采用非均匀有限元离散模型进行计算. 图11给出了H20和T10单元的网格划分情况, 图12给出了相应的前4阶频率计算结果. 由图12的结果可见, 即使对于非均匀网格划分, H20和T10单元的MNLM集中质量矩阵的频率计算精度同样明显优于HRZ集中质量矩阵, 说明H20和T10单元的MNLM集中质量矩阵具有良好的鲁棒性.
![]()
图 11 1/4圆筒空腔问题有限元网格Figure 11. Finite element meshes of the quarter-cylinder cavity problem![]()
图 12 1/4圆筒空腔问题前四阶频率收敛特性对比Figure 12. Convergence comparison of the first four frequencies for the quarter-cylinder cavity problem对于该问题, 我们同样进行时程分析, 其中考虑的空腔内解析声压解为
$$ \left.\begin{split} & u({\boldsymbol{x}},t) = 10F(r)\sin (2\theta )\sin (\text{π} z)\sin (\sqrt {{\varpi ^2} + {\text{π} ^2}} t) \\ & F(r) = {\rm{J}}_0^{[1]}(\varpi {r_i}){\rm{J}}_0^{[2]}(\varpi r) - {\rm{J}}_0^{[2]}(\varpi {r_i}){\rm{J}}_0^{[1]}(\varpi r) \\ & \varpi = {\text{3}}{{.123\;030\;919\;558\;37}} \end{split} \right\} $$ (70) 有限元分析的初始条件可通过在式(70)中令$ t = 0 $给出. 同样, 式(11)中源项$ b $也可由式(70)求得
$$ b({\boldsymbol{x}},t) = \ddot u - {c^2}{\nabla ^2}u = \frac{{40}}{{{r^2}}}u $$ (71) 在时程分析中, 有限元模型采用图11中的1024个H20单元和5120个T10单元, 时间步长为$ \Delta t = 0.01 $. 图13给出了${\boldsymbol{x}} = (1.5,\;\text{π} /4,0.5)$处的声压数值解$ {u^h} $的时程曲线及误差对比结果, 图14给出了$ t = 12 $时刻的位移误差云图. 结果显示, 对于时程分析, H20和T10单元的MNLM集中质量矩阵的计算精度同样优于对应的HRZ集中质量矩阵.
![]()
图 13 1/4圆筒空腔问题声压时程对比Figure 13. Comparison of the acoustic pressure time history for the quarter-cylinder cavity problem![]()
图 14 1/4圆筒空腔问题$ t = 12 $时刻的声压误差云图Figure 14. Acoustic pressure error contour plots for the quarter-cylinder cavity problem at $ t = 12 $与此同时, 图15给出了MNLM和HRZ集中质量矩阵的计算效率对比. 该问题的计算效率结果再次表明, 对于具有非均匀网格特征的1/4圆筒空腔问题, 相较于HRZ集中质量矩阵, H20和T10单元的MNLM集中质量矩阵的计算效率分别提升了约1.3倍和3.0倍, 实现了计算效率的显著提升.
![]()
图 15 1/4圆筒空腔时程分析效率对比Figure 15. Efficiency comparison of the transient analysis for the quarter-cylinder cavity problem6. 结论
本文针对结构分析中应用广泛的三维20节点六面体单元和10节点四面体单元, 系统研究了其集中质量矩阵构造方法和精度. 通过引用一种包含待定参数的三维20节点六面体单元广义集中质量矩阵, 建立了相应的频率精度表达式, 进而揭示了不同集中质量矩阵的理论精度. 研究表明, 对角元素比例缩放法, 即HRZ方法, 作为所提广义集中质量矩阵构造方法的一个特例, 虽然能够保证集中质量矩阵元素非负性, 但其精度并非最优. 随后, 通过对广义集中质量矩阵频率精度进行优化, 提出了一种精度更优的三维20节点六面体单元中节点集中质量矩阵构造方法, 并将其推广到10节点四面体单元. 中节点集中质量矩阵同样具有非负特性, 但角节点的质量为零, 可以方便地进行模型降阶. 典型算例的频率和时程计算结果均表明, 与HRZ集中质量矩阵相比, 三维20节点六面体单元和10节点四面体单元的中节点集中质量矩阵在计算精度和效率方面都得到了显著提升. 以20节点六面体单元为例, 就精度而言, 如表1所示, 中节点集中质量矩阵的频率计算误差比HRZ集中质量矩阵降低了约1/3; 就效率而言, 如图9和图15所示, 对于弹性力学问题和势问题, 中节点集中质量矩阵的时程分析计算效率约为HRZ集中质量矩阵的5倍和2倍.
-
![]()
图 1 H20单元3种集中质量矩阵的节点质量分布示意图
Figure 1. Schematic illustration of the nodal mass distribution for three mass lumping schemes of H20 element
![]()
图 2 T10单元3种集中质量矩阵的节点质量分布示意图
Figure 2. Schematic illustration of the nodal mass distribution for three mass lumping schemes of T10 element
![]()
图 3 H20单元典型节点影响域及波动方向示意图
Figure 3. Illustration of the typical nodal influence domains and wave propagation angles for a mesh with H20 elements
![]()
图 5 长方体空腔问题前四阶频率收敛特性对比
Figure 5. Convergence comparison of the first four frequencies for the rectangular cavity problem
![]()
图 7 长方体弹性体前六阶频率收敛特性对比
Figure 7. Convergence comparison of the first six frequencies for the rectangular elastic continuum problem
![]()
图 8 长方体弹性体$ t = 10 $时刻的位移误差云图
Figure 8. Displacement error contour plots for the rectangular elastic continuum problem at $ t = 10 $
![]()
图 9 长方体弹性体时程分析效率对比
Figure 9. Efficiency comparison of the transient analysis for the rectangular elastic continuum problem
![]()
图 11 1/4圆筒空腔问题有限元网格
Figure 11. Finite element meshes of the quarter-cylinder cavity problem
![]()
图 12 1/4圆筒空腔问题前四阶频率收敛特性对比
Figure 12. Convergence comparison of the first four frequencies for the quarter-cylinder cavity problem
![]()
图 13 1/4圆筒空腔问题声压时程对比
Figure 13. Comparison of the acoustic pressure time history for the quarter-cylinder cavity problem
![]()
图 14 1/4圆筒空腔问题$ t = 12 $时刻的声压误差云图
Figure 14. Acoustic pressure error contour plots for the quarter-cylinder cavity problem at $ t = 12 $
![]()
图 15 1/4圆筒空腔时程分析效率对比
Figure 15. Efficiency comparison of the transient analysis for the quarter-cylinder cavity problem
表 1 长方体空腔问题前四阶频率计算精度对比
Table 1 Accuracy comparison of the first four frequencies for the rectangular cavity problem
$ {\omega _i} $ $ {n_{el}}( \times 5) $ $\omega _i^h{ { \text{-} } }{\rm{H} }20$ $\left| {\omega _i^h{{ - } }\omega _i^{} } \right|$ $\left| {\dfrac{ { {e_{{\rm{MNLM}}} } } }{ { {e_{{\rm{HRZ}}} } } } } \right|$ $\omega _i^h{ { \text{-} } }{\rm{T} }10$ $\left| {\omega _i^h{{ - } }\omega _i^{} } \right|$ $\left| {\dfrac{ { {e_{{\rm{MNLM}}} } } }{ { {e_{{\rm{HRZ}}} } } } } \right|$ HRZ MNLM HRZ MNLM HRZ MNLM HRZ MNLM 4.71 216 4.45 4.54 0.26 0.17 0.66 4.64 4.68 0.08 0.03 0.38 512 4.56 4.61 0.15 0.10 0.66 4.67 4.69 0.04 0.02 0.39 1000 4.62 4.65 0.10 0.06 0.66 4.68 4.70 0.03 0.01 0.40 5.44 216 5.04 5.17 0.41 0.27 0.66 5.32 5.40 0.12 0.04 0.36 512 5.21 5.29 0.23 0.15 0.66 5.37 5.42 0.07 0.03 0.38 1000 5.29 5.34 0.15 0.10 0.66 5.40 5.42 0.04 0.02 0.39 6.48 216 5.77 6.01 0.70 0.46 0.66 6.28 6.41 0.20 0.07 0.35 512 6.08 6.21 0.40 0.26 0.66 6.36 6.43 0.11 0.04 0.37 1000 6.22 6.31 0.25 0.17 0.66 6.40 6.45 0.07 0.03 0.39 7.20 216 6.32 6.62 0.88 0.58 0.66 6.94 7.11 0.26 0.09 0.34 512 6.68 6.86 0.52 0.34 0.66 7.05 7.14 0.15 0.06 0.37 1000 6.86 6.97 0.34 0.22 0.66 7.10 7.16 0.10 0.04 0.39 
下载: 导出CSV
-
[1] Hughes TJR. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications, 2000
[2] 谢臻, 黄钟民, 陈思亚等. 基于无网格法的Timoshenko曲梁动力分析. 固体力学学报, 2022, 43(5): 603-613 (Xie Zhen, Huang Zhongmin, Chen Siya, et al. Dynamic analysis of Timoshenko curved beams based on meshless method. Chinese Journal of Solid Mechanics, 2022, 43(5): 603-613 (in Chinese) doi: 10.19636/j.cnki.cjsm42-1250/o3.2022.016 Xie Z, Huang Z, Chen S, Peng L. Dynamic analysis of Timoshenko curved beams based on meshless method. Chinese Journal of Solid Mechanics, 2022, 43(5): 603-613. (in Chinese)). doi: 10.19636/j.cnki.cjsm42-1250/o3.2022.016
[3] 陈健, 王东东, 刘宇翔等. 无网格动力分析的循环卷积神经网络代理模型. 力学学报, 2022, 54(3): 732-745 (Chen Jian, Wang Dongdong, Liu Yuxiang, et al. A recurrent convolutional neural network surrogate model for dynamic meshfree analysis. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(3): 732-745 (in Chinese) doi: 10.6052/0459-1879-21-565 Chen J, Wang D, Liu Y, Chen J. A recurrent convolutional neural network surrogate model for dynamic meshfree analysis. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(3): 732-745. (in Chinese)). doi: 10.6052/0459-1879-21-565
[4] Stoter S, Nguyen TH, Hiemstra RR, et al. Variationally consistent mass scaling for explicit time-integration schemes of lower- and higher-order finite element methods. Computer Methods in Applied Mechanics and Engineering, 2022, 399: 115310 doi: 10.1016/j.cma.2022.115310
[5] Lan M, Yang W, Liang X, et al. Vibration modes of flexoelectric circular plate. Acta Mechanica Sinica, 2022, 38(12): 422063 doi: 10.1007/s10409-022-22063-x
[6] 黄可, 张家应, 王青云. 基于非均匀梁模型的二维柔性机翼固有振动分析. 力学学报, 2023, 55(2): 784-496 (Huang Ke, Zhang Jiaying, Wang Qingyun. Natural vibration analysis of two-dimensional flexible wing based on non-uniform beam model. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(2): 784-496 (in Chinese) doi: 10.6052/0459-1879-22-551 Huang K, Zhang J, Wang Q. Natural vibration analysis of two-dimensional flexible wing based on non-uniform beam model. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(2): 784-496. (in Chinese)). doi: 10.6052/0459-1879-22-551
[7] Zienkiewicz OC, Taylor RL. The Finite Element Method: Its Basis and Fundamentals. Butterworth-Heinemann, 2013
[8] Felippa CA, Onate E. Accurate Timoshenko beam elements for linear elastostatics and LPB stability. Archives of Computational Methods in Engineering, 2021, 28: 2021-2080 doi: 10.1007/s11831-020-09515-0
[9] Wang D, Liu W, Zhang H. Novel higher order mass matrices for isogeometric structural vibration analysis. Computer Methods in Applied Mechanics and Engineering, 2013, 260: 92-108 doi: 10.1016/j.cma.2013.03.011
[10] Wang D, Liang Q, Wu J. A quadrature-based superconvergent isogeometric frequency analysis with macro-integration cells and quadratic splines. Computer Methods in Applied Mechanics and Engineering, 2017, 320: 712-744 doi: 10.1016/j.cma.2017.03.041
[11] Wang D, Pan F, Xu X, et al. Superconvergent isogeometric analysis of natural frequencies for elastic continua with quadratic splines. Computer Methods in Applied Mechanics and Engineering, 2019, 347: 874-905 doi: 10.1016/j.cma.2019.01.010
[12] Fried I, Malkus DS. Finite element mass lumping by numerical integration with no convergence rate loss. International Journal of Solids and Structures, 1975, 11(4): 461-466 doi: 10.1016/0020-7683(75)90081-5
[13] Li X, Wang D. On the significance of basis interpolation for accurate lumped mass isogeometric formulation. Computer Methods in Applied Mechanics and Engineering, 2022, 400: 115533 doi: 10.1016/j.cma.2022.115533
[14] Li X, Zhang H, Wang D. A lumped mass finite element formulation with consistent nodal quadrature for improved frequency analysis of wave equations. Acta Mechanica Sinica, 2022, 38(5): 521388 doi: 10.1007/s10409-021-09022-x
[15] Li Y, Liang R, Wang D. On convergence rate of finite element eigenvalue analysis with mass lumping by nodal quadrature. Computational Mechanics, 1991, 8(4): 249-256 doi: 10.1007/BF00577378
[16] Ainsworth M, Wajid HA. Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration. SIAM Journal on Numerical Analysis, 2010, 48(1): 346-371 doi: 10.1137/090754017
[17] Hinton E, Rock T, Zienkiewicz OC. A note on mass lumping and related processes in the finite element method. Earthquake Engineering and Structural Dynamics, 1976, 4(3): 245-249 doi: 10.1002/eqe.4290040305
[18] Yang Y, Zheng H, Sivaselvan MV. A rigorous and unified mass lumping scheme for higher-order elements. Computer Methods in Applied Mechanics and Engineering, 2017, 319: 491-514 doi: 10.1016/j.cma.2017.03.011
[19] Duczek S, Gravenkamp H. Mass lumping techniques in the spectral element method: On the equivalence of the row-sum, nodal quadrature, and diagonal scaling methods. Computer Methods in Applied Mechanics and Engineering, 2019, 353: 516-569 doi: 10.1016/j.cma.2019.05.016
[20] Li X, Wang D, Xu X, et al. A nodal spacing study on the frequency convergence characteristics of structural free vibration analysis by lumped mass Lagrangian finite elements. Engineering with Computers, 2022, 38: 5519-5540 doi: 10.1007/s00366-022-01668-9
[21] Ozakca M, Hinton E, Rao NVR. Comparison of three-dimensional solid elements in the analysis of plates. Computers & Structures, 1992, 42(6): 953-968
[22] Wu HP, Shang Y, Cen S, et al. Penalty C0 8-node quadrilateral and 20-node hexahedral elements for consistent couple stress elasticity based on the unsymmetric finite element method. Engineering Analysis with Boundary Elements, 2023, 147: 302-319 doi: 10.1016/j.enganabound.2022.12.008
[23] Choi JH, Lee BC, Sim GD. A 10-node tetrahedral element with condensed Lagrange multipliers for the modified couple stress theory. Computers & Structures, 2021, 246: 106476
[24] Hanukah E, Givli S. A new approach to reduce the number of integration points in mass-matrix computations. Finite Elements in Analysis and Design, 2016, 116: 21-31 doi: 10.1016/j.finel.2016.03.004
[25] Lo SH, Ling C. Improvement on the 10-node tetrahedral element for three-dimensional problems. Computer Methods in Applied Mechanics and Engineering, 2000, 189(3): 961-974 doi: 10.1016/S0045-7825(99)00410-7
[26] Hou S, Li X, Wang D, et al. A mid-node mass lumping scheme for accurate structural vibration analysis with serendipity finite elements. International Journal of Applied Mechanics, 2021, 13(1): 2150013 doi: 10.1142/S1758825121500137
[27] Bathe KJ. Finite Element Procedures. Prentice Hall, 1996
[28] Rao SS. Vibration of Continuous Systems. John Wiley & Sons, 2007
[29] Day DM, Romero L. An analytically solvable eigenvalue problem for the linear elasticity equations. UNT Libraries Government Documents Department, 2004
[30] 吴俊超, 吴新瑜, 赵珧冰等. 基于赫林格-赖斯纳变分原理的一致高效无网格本质边界条件施加方法. 力学学报, 2022, 54(12): 3283-3296 (Wu Junchao, Wu Xinyu, Zhao Yaobing, et al. A consistent and efficient method for imposing meshfree essential boundary conditions via Hellinger-Reissner variational principle. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(12): 3283-3296 (in Chinese) doi: 10.6052/0459-1879-22-151 Wu J, Wu X, Zhao Y, Wang D. A consistent and efficient method for imposing meshfree essential boundary conditions via Hellinger-Reissner variational principle. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(12): 3283-3296. (in Chinese)). doi: 10.6052/0459-1879-22-151
-
期刊类型引用(2)
1. 冀佳. 建筑电气设备铝合金壳体的焊接技术. 焊接技术. 2023(10): 92-95 . 
百度学术
2. 顾崴，刘铖，安志朋，史东华. 一种基于Hamel形式的无条件稳定动力学积分算法. 力学学报. 2022(09): 2577-2587 . 本站查看
其他类型引用(1)
计量
- 文章访问数: 587
- HTML全文浏览量: 315
- PDF下载量: 103
- 被引次数: 3
