2.3 聚集单元谱元法在承受冲击载荷结构动态分析中的应用

在工程中，机械结构几乎都要承受动态载荷，并且在多数情况下会承受冲击载荷，如用榔头敲击螺钉，声共振无损检测中敲击锤敲击被检测构件，汽车的碰撞，船舶与桥梁、海洋平台的碰撞，飞机的坠落，水面舰艇遭受水下非接触性爆炸的冲击［89］等，为了评价其动态特性及动态行为，有必要对其进行仿真分析。目前，国内外对冲击载荷作用下的结构动态行为研究较多，如起重机机架在起升冲击载荷作用下的动态特性分析［90］、冲击载荷作用下齿轮动态应力变化研究［91］、冲击载荷作用下蜂窝夹芯板的动力响应分析［92］等。很多文章对各种结构在冲击载荷作用下的动态特性及其动态行为进行了分析研究，然而还未见运用谱元法对结构在冲击载荷作用下的动态分析的相关报道。

谱元法是Patera于1984年提出的应用于CFD的一种数值方法，其具有有限元方法处理任意结构及其边界的灵活性和谱方法快速收敛的优点［81］。谱元法能用较少的单元获得与其他方法相同的精度，其特点是将每个单元在GLL的零点处离散，然后进行Lagrange多项式插值。从理论上分析，在正交多项式零点处插值时可获得的插值精度最高［86］。笔者从减小计算规模的角度提出了基于逐步时间谱元法的结构动态响应仿真方法，将仿真时间划分为很小的时间段，在每个时间段内划分单元，将每个时间段看作独立计算部分，并将前一部分的计算结果作为后一部分计算的初始条件，该方法节省了计算时间［53］；在第二类Chebyshev正交多项式零点处，从重心Lagrange插值角度构造了非线性振动问题的离散方案［54］，利用谱单元离散插值精度高的特点，计算了动态响应优化中的关键时间点［38］。M.H.Kurdi［93］利用时间谱元法求解了简单的质量弹簧阻尼系统，并在此基础上对单自由度吸振器和单自由度微型控制器进行了优化。

本章从结构动力学控制方程出发，利用有限元方法进行空间离散，用谱元法进行时间离散，其中针对冲击载荷时间短、变化大的特点，将谱单元的离散进行聚集处理，以弥补等距单元误差大的缺陷。

谱元法分为时间谱元法、空间谱元法和时间-空间耦合谱元法。对于动力学方程而言，可以采用其中任何一种方法。下面以时间谱元法为例进行介绍。

时间谱元法的主要步骤：①将动力学方程转化为一阶线性微分方程组，然后通过Bubnov-Galerkin法等效为积分形式，代入单元积分表达式，获得时间单元谱元方程；②根据冲击载荷的特点，将仿真时间划分为聚集型单元，即在载荷突变的位置单元尺寸比较小，在载荷平坦的位置单元尺寸比较大，它可通过单元最大尺寸和最小尺寸的比值来控制；③将每个时间单元划分为若干个时间节点，即正交多项式的零点；④将每个单元的近似解表示为Legendre正交多项式的线性组合；⑤利用连接矩阵将所有单元的近似解集成总体谱元方程；⑥求解总体谱元方程，获得全局位移近似解和速度近似解，可以通过动力学微分方程得到加速度的解，对于结构来说，可以通过应力和位移的关系获得节点应力解。

2.3.1 线性结构动态响应方程及其转化形式

线性结构动态响应方程可表示为

式中，M为质量矩阵；C为阻尼矩阵；K为刚度矩阵；F为动态载荷向量；x为位移向量；为加速度向量，为速度向量。初始条件为x（0）=b0，。式（2.38）中，M、C、K是不随时间变化的；x，，是时间的函数；F是任意时间函数，时间t∈［t 0，t n］。

为了采用时间谱元法，设x1=，x2=x，式（2.38）可转化为一阶线性常微分方程组［见式（2.39）］，此方程组与式（2.38）同解。

2.3.2 聚集单元划分

如图2.5所示，将仿真时间t∈［t 0，tn］分割为n个互不相交的单元，即［t0，t1］，［t1，t2］，［t 2，t3］，…，［t n-1，tn］，每个单元配置若干个点。在划分单元时，以冲击载荷最大值点为中心，中心处单元尺寸最小，越往两边单元尺寸越大。每个单元配置Chebyshev或Legendre正交多项式的零点，其点数可以相同，也可以不同，本章采用相同的零点配置。这样可以避免冲击载荷的局部突变带来的求解误差。在具体实施时，可以设（图2.5中ratio=0.1）。图2.5（a）所示为聚集单元示意，图2.5（b）所示为30个聚集单元对应的冲击载荷，其中，冲击载荷由3个参数控制，即冲击载荷最大值、宽度和载荷中心。

2.3.3 单元分析

聚集单元划分后，在单元区间内增加正交多项式的零点，并且可以通过增加插值点来提高近似精度，这在有限元方法中称为hp法。

在时间谱元法中，通过在单元内的特殊位置配置插值GLL点，其基函数可表示为正交多项式的组合，从而构成单元内各点的形函数，这样可以在有限个点上插值，达到谱方法的收敛速度。

在时间谱元法中一般采用两种正交多项式，即Chebyshev和Legendre正交多项式，本章采用Legendre正交多项式［86］。

由于单元端点不是Legendre正交多项式的零点，因此加入Lobbato多项式来保证单元端点包含在零点中。Lobbato多项式满足正交特性：

式中，δij是Kronecker δ函数；

通过解式（2.40）可以获得GLL点及其权值。

式中，N为插值次数，权值可以表示为

式中，ξk为LoN-1（ξ）的零点，k=1，2，…，N-1。

在仿真时间内按图2.5所示方案划分，并且结合冲击载荷离散，将状态变量表示为m次Lagrange多项式：

式中，为j单元的k次Lagrange多项式；ξk为定义在［-1，1］上的GLL点；是单元j上未知节点在GLL点的值。然后通过Bubnov-Galerkin法可得

用矩阵形式表示每个单元的离散：

式中，，Iω为ξ的一般函数，

2.3.4 集成总体谱元方程及求解

集成是把所有单元的谱元方程按照离散的顺序组合起来，获得总体谱元方程。对含Nv个状态方程的系统，通过耦合矩阵As（Nv×Nv的方阵）的张量叉乘可得到全部状态变量的全局组装：

式中，Bu为全局微分矩阵；Bω为全局权矩阵；Fu g（Xug）为激励力的全局形式；Xug为所有状态变量在时间节点的集合。化简式（2.46）得

式中，G为时间段的全局线性矩阵；式（2.47）是线性方程组，可通过高斯消元法求解。