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

毛虎平，乔文元，郭保全，王 强，董小瑞

（1.中北大学 机械与动力工程学院，太原 030051； 2.中北大学 机电工程学院，太原 030051）

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

毛虎平1，乔文元1，郭保全2，王 强1，董小瑞1

（1.中北大学 机械与动力工程学院，太原 030051； 2.中北大学 机电工程学院，太原 030051）

为了扩展谱元法的适应性，针对承受冲击载荷的结构动态响应问题，从谱单元离散方案出发并根据冲击载荷的特点，以冲击载荷最大值点为中心将谱单元尺寸按一定比例等比向两侧扩大，实现单元尺寸与载荷特征相适应，在此基础上，将动力学方程转化为1阶线性微分方程组，通过Bubnov-Galerkin法获得离散线性方程组，应用高斯消元法求解。最后，应用此方法对标准形式、线性单自由度吸振器、124杆桁架结构和连杆小头问题进行计算并与等距谱元法进行比较，说明了此方法的可行性和有效性。

振动与波；聚集单元；谱元法；冲击载荷；动态分析

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

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

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

1 线性动力学方程的时间谱元法

谱元法分为时间谱元法、空间谱元法和时空混合谱元法。对于动力学方程而言，可以采用三种方法中的任何一种，文中利用谱元法离散时间，采用有限元法离散空间。

时间谱元法的主要步骤为

（1）将动力学方程转化为一阶线性微分方程组，然后通过布巴诺夫-伽辽金法（Bubnov-Galerkin）等效为积分形式，代入单元积分表达式，获得时间单元谱元方程；

（2）根据冲击载荷的特点，将仿真时间划分为聚集型单元，即在载荷突变的位置单元尺寸比较小，在载荷平坦的位置单元尺寸比较大，这可通过单元最大尺寸和最小尺寸的比值来控制；

（3）将每个时间单元划分为若干个时间节点，即正交多项式的零点；

（4）将每个单元的近似解表示为Legendre正交多项式的线性组合；

（5）利用连接矩阵将所有单元的近似解组装成总体谱元方程；

（6）求解总体谱元方程，获得全局位移近似解和速度近似解，可以通过动力学微分方程得到加速度的解，对于结构来说，可以通过应力和位移的关系获得节点应力解。

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

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

其中M为质量矩阵，C为阻尼矩阵，K为刚度矩阵，F为动态载荷向量，x为位移向量，为加速度向量，为速度向量。初始条件为x(0)=b0，其中M、C、K不随时间变化；是时间的函数，F是任意时间函数

1.2 聚集单元划分

如图1所示，将仿真时间t∈[t0,tn]分割为n个互不相交的单元，即每个单元配置若干个点。在划分单元时，以冲击载荷最大值点为中心，中心处单元尺寸最小，越往两边尺寸越大，如图1所示。

图1 聚集单元谱元法离散时间及对应的冲击载荷

每个单元配置Chebyshev或Legendre正交多项式的零点，其点数可以相同，也可不同，文中采用配置相同的零点。这样可以避免由于冲击载荷的局部突变性带来的求解误差。具体实施时，通过设来实现，图1中ratio=0.1。图1 (a)为聚集单元示意图，图1(b)表示30个聚集单元对应的冲击载荷，其中冲击载荷可以有3个参数控制，即冲击载荷最大值、宽度和载荷中心。

1.3 单元分析

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

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

在时间谱元法中一般采用两种正交多项式，即Chebyshev和Legendre正交多项式，文中采用Legendre正交多项式[6]。

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

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

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

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

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

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

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

2 实例分析

为了验证聚集单元谱元法的优势，对4个不同实例进行分析，并与等距单元谱元法比较。

2.1 标准形式

对于突变载荷，常规的等长度单元谱元法不能获得好的收敛效果。冲击载荷如图2所示。

图2 冲击载荷

无论多么复杂的系统，动态响应方程最终可以转化为形如式（11）所示的标准形式。其中，载荷在0.5 s时刻，产生很大的冲击，冲击的幅值可以由系数决定，冲击力的作用时间可以通过调整ε来改变，本例中ε=0.000 112 5。解析解可以通过积分因子法并借用Matlab中的erf函数获得

图3由等距单元谱元法获得，很明显误差很大，特别在载荷突变处，Ne=30时，误差极大，而当Ne=100时，误差减小了很多，但还是较大。

图3 等距单元谱元法计算结果

图4由聚集单元谱元法获得，在单元数和插值次数相等的前提下已经消除了误差。

图4 聚集单元谱元法计算结果

图5和图6更能说明聚集单元谱元法的优势，图5将等距单元谱元法与聚集单元谱元法比较，发现在50个单元处，单元插值次数为18时前者为（50，0.02506）；而单元插值次数为8时，后者为（50，1.382×10-10），有非常大的优势；图6将聚集单元谱元法与等距单元谱元法比较，发现在80个单元处，单元插值次数分别为5、10、18时前者为（80，3×10-9）、（80，3.628× 10-12）、（80，2.668×10-12），而后者在80个单元处，单元插值次数为18次时为（80，0.011 24），可以看出聚集单元谱元法求解冲击载荷响应时有很大优势。

图5 聚集单元与等距单元比较

图6 等距单元与聚集单元比较

2.2 线性单自由度吸振器

图7为包含固定质量m=1 kg、弹簧刚度系数和阻尼器系数分别为k=0.9、c=0.9的线性单自由度系统。

图7 缓冲器

当t=0时，系统以速度v=1 m/s撞击在一个固定阻碍物上，同时受冲击载荷作用，ε=0.011 25。运动方程为

式（12）自由振动的解析解为

线性单自由度吸振器虽然既有初始速度，又同时受冲击载荷，但冲击载荷的中心在6 s处，宽度由ε=0.011 25来控制而且很窄，因此在4 s附近还是属于自由振动，可以与自由振动解析解比较，如图8所示，从图中可以看出单元数为12，单元插值次数为12时，聚集单元谱元法明显比等距单元谱元法精度高，说明本文方法具有优势。

图8 线性单自由度吸振器动态位移响应

2.3 杆桁架结构

该平面桁架结构包含49个铰链，94个自由度，见图9。弹性模量E=207 GPa，泊松比v=0.3，密度ρ=7 850 kg/m3，杆的截面积0.645×10-4m2。在节点1、20、19、18、17、16、15的x正方向上作用相同的动态载荷，在节点1、2、3、4、5的y负方向上也作用相同的动态载荷，b=1 000，ε=0.000 112 5。

图9 124杆桁架结构

从工程角度，将节点1、2、3、4、5、15、16、17、18、19、20节点作为重点考察的位置，称之为关键位置，同时也是冲击载荷作用的位置。在冲击载荷作用下，其位移响应如图10所示，将其中节点1的X方向位移单独考虑，并局部放大如图11所示。

图10 124杆桁架结构关键位置节点X方向动态位移响应

图11 124杆桁架结构节点1X方向动态位移响应

可以明显看出等距单元谱元法计算结果有部分点偏离精确值，而聚集单元谱元法没有这种现象，实际上关键位置的每一个点都出现类似现象，限于篇幅，没有全部列出。

2.4 连杆小头

连杆小头及其相连的杆身部分，如图12所示，采用PLANE42单元，右端固定，共有111个节点、212个自由度，在小头内圆内中间受水平向右的冲击载荷，载荷如图2所示，b=106，ε=0.000 112 5，连杆材料的弹性模量E= 207 GPa，泊松比v=0.3，密度ρ=7 850 kg/m3。

图12 连杆小头与杆身相连部分

连杆小头圆角处为关键位置，考察节点15、42、49、93X方向的位移响应。图13中，将模态叠加法获得结果看作精确解，可以看出在相同单元数量、尺寸、插值次数下，聚集单元谱元法更精确，和模态叠加法保持一致，而等距单元谱元法存在较大误差。由于位移响应本身就很小，因此进行了局部放大。

图13 连杆4个关键位置节点X方向动态位移响应

3 结语

研究结论如下。

（1）针对冲击载荷的特点，并结合谱元法精度高的优点，采用设定载荷中心单元尺寸与两侧单元尺寸的比值来适应载荷的突变性，使单元的离散与冲击载荷变化相适应；

（2）在单元数量、尺寸、插值次数相同的前提下，求解冲击载荷动态响应，聚集单元谱元法的精度比等距单元谱元法高；

（3）从动态问题的标准形式，到线性单自由度吸振器、124杆平面桁架，再到连杆小头问题的冲击载荷动态响应分析，都说明了本文方法的可行性和有效性。这为进一步研究冲击载荷作用下的结构动态分析提供参考。

[1]肖锋，谌勇，马超，等.橡胶蜂窝覆盖层水下爆炸响应及抗冲击性能[J].噪声与振动控制，2013，33(4)：44-49.

[2]王贡献，沈荣瀛.起重机臂架在起升冲击载荷作用下动态特性研究[J].机械强度，2005，27(5)：561-566.

[3]唐进元，彭方进，黄云飞.冲击载荷下的齿轮动应力变化规律数值分析[J].振动与冲击，2009，28(08)：138-143.

[4]李永强，宦强.基于Ansys的冲击载荷下蜂窝夹芯板的动力学响应[J].东北大学学报(自然科学版)，2015：36(6)：858-862.

[5]PATERA A T.A spectral element method for fluid dynamics:Laminarflow in achannelexpansion[J].Journal of Computational Physics,1984,54(3):468-488.

[6]POZRIKIDIS C.Introduction to finite and spectral element methods using Matlab[M].CRC Press,2005.

[7]毛虎平，苏铁熊，李建军.基于逐步时间谱元法的结构动态响应仿真[J].中北大学学报(自然科学版)，2013(4)：424-430.

[8]毛虎平，王伟能，续彦芳，等.非线性振动分析的切比雪夫谱元法[J].噪声与振动控制，2015，35(1)：73-77.

[9]张艳岗，毛虎平.动态应力解空间谱元离散的关键时间点识别方法[J].机械工程学报，2014，50(5)：82-87.

[10]KURDI M H,BERAN P S.Optimization of dynamic response using a monolithic-time formulation[J].Structural and Multidisciplinary Optimization,2009, 39(1):83-104.

Application of SEM based onAccumulation Elements to DynamicAnalysis of the Structures Subjected to Impact Loads

MAO Hu-ping1,QIAO Wen-yuan1,GUO Bao-quan2, WANG Qiang1,DONG Xiao-rui1
(1.College of Mechanical and Power Engineering,North University of China,Taiyuan 030051,China; 2.College of Mechatronic Engineering,North University of China,Taiyuan 030051,China)

The spectral element method(SEM)is used to analyze the dynamic responses of the structures subjected to impact loading.Starting with the SEM discrete scheme and according to the feature of the impact loads,each spectral element is enlarged proportionally from the middle point of the impact load to both sides of the element to realize the suitability between the element size and the loading feature.On this basis,the dynamic equation is transformed to a set of the first-order linear differential equations,and the corresponding discrete linear algebraic equations are obtained by Bubnov-Galerkin method and solved by Gauss elimination method.Finally,four examples of a typical absorber,a linear single DOF absorber,a 124-member plane truss and the small end of a link are calculated with this method and the results are compared with those of the equal-spaced SEM.The feasibility and validity of this method are verified.

vibration and wave;accumulation element;spectral element method(SEM);impact load;dynamic analysis

TH122

：A

：10.3969/j.issn.1006-1335.2016.06.009

1006-1355(2016)06-0045-06

2016-07-04

国家自然科学基金资助项目(51275489)

毛虎平（1974-），男，博士，副教授，硕士生导师，主要从事振动理论与工程数值分析方法、结构动态响应优化方法研究。E-mail:maohp@nuc.edu.cn