刘赛 , 吕振华
清华大学 汽车工程系, 北京 100084

收稿日期: 2016-01-24
作者简介: 刘赛(1988-), 男, 博士研究生
通信作者: 吕振华, 教授, E-mail:lvzh@tsinghua.edu.cn

摘要：已有文献进行了微弓形扁长金属杆在轴向重物冲击下的动态屈曲实验研究，得到了扁长杆冲击屈曲响应的典型结果，但已有文献及后来相关论文的模拟计算精度较低。为了提高典型扁长杆的冲击屈曲特性的计算分析精度，该文研究了扁长杆的冲击屈曲特性分析的有限元建模方法，改进了有限元模型边界条件的建模仿真度（考虑了转动铰配合间隙、摩擦等），并探讨了有限单元型式及尺度的选择，以壳单元或实体单元模型代替梁单元模型。研究表明：采用改进边界条件（考虑转动铰配合间隙、摩擦等）和薄壳单元或厚壳单元或实体单元的仿真模型比采用理想边界条件和梁单元的原模型的计算结果精度显著提高；基于改进模型的计算分析结果，揭示了该扁长杆受冲击载荷作用时的3维动态反向屈曲行为。
关键词： 动态屈曲 扁长杆 轴向冲击 有限元分析 模型修正 3维反向屈曲
Finite element model refinement for elastic-plastic dynamic buckling of a belt bar during impact
LIU Sai, Lü Zhenhua
Department of Automotive Engineering, Tsinghua University, Beijing 100084, China


Abstract:The dynamic buckling experiment of a belt bar subjected to axial impact was previously studied experimentally. The tests show that previous simulation methods are not very accurate. The finite element analysis accuracy is improved by improving the fidelity of the boundary conditions and properly selecting the finite element type and size. The boundary conditions for the joint clearance and friction are determined by comparing the numerical results with the test data, and shell elements and solid elements are found to give better results than beam elements. Simulations with thin or thick shell elements or solid elements with the new boundary conditions are more accurate than previous results with beam elements and simple boundary conditions. The refined finite element model predicts the 3-D reverse buckling of the belt bar.
Key words: dynamic bucklingbelt baraxial impactfinite element analysismodel refinement3-D reverse buckling
随着计算结构动力学技术的发展，采用显式动力学方法进行结构的动态屈曲分析的研究逐渐增多。文[1]给出了微弓形扁长金属杆在轴向受重物冲击的典型屈曲实验结果，但该文及后来相关论文[2]的模拟计算精度较低。其他杆结构动态屈曲特性的计算分析均存在类似问题，其原因主要在于建模缺陷：1)?解析计算模型^[3-7]和仿真计算模型^{[1-2, 8-10]}均采用理想边界条件，与真实的边界条件(转动铰存在配合间隙、摩擦等)之间存在差异；2)?研究扁长金属杆冲击弹塑性屈曲特性的计算模型均采用了梁单元模型^{[1-2, 8, 11-12]}，精度较低。因此，这些研究结果未能较好地再现实验结果的一些主要特征。
本文针对文[1]的实验研究对象，改进计算分析模型的边界条件，并探讨有限单元型式及尺度的选择，研究它们对冲击屈曲变形模式(屈曲模态)、失稳载荷和失稳时间的影响，以改善扁长杆的冲击弹塑性屈曲特性的有限元分析精度。
1 扁长杆的冲击屈曲实验及仿真分析结果在文[1]的实验中，矩形截面扁长杆的一端为铰支，另一端为自由端，如图 1所示。一质量为369 g、速度为2.45 m/s的重物撞击扁长杆的自由端，但重物的侧向运动受到约束，使得杆的自由端运动等效于滑动铰支移动。杆件长度为193 mm，横截面尺寸为8.7 mm×0.63 mm，扁长杆预制为微弓形，其定义见表 1。杆件材料为Ni-Cr合金钢，密度为7.8 g/cm³，弹性模量为206 GPa，屈服强度为1 470 MPa。扁长杆的屈曲变形曲线(时刻t=0.75 ms，下同)和轴向力(位置x=0.1，下同)见图 2，其中扁长杆的屈曲变形-时间历程是通过高速摄影图像处理得到的，轴向力是通过扁长杆上x=0.1处两侧的表面应变测量结果计算得到的。

图 1 受重物冲击的微弓形扁长杆[1, 8]

图选项

表 1 扁长杆的微弓形弧函数^[1]

x	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
y0/mm	0	0.159	0.285	0.390	0.459	0.480	0.459	0.391	0.285	0.160	0
注：x为x的归一化坐标，即x=x/(193 mm).

表选项

图 2 基于梁单元模型和理想边界条件的计算分析结果

图选项

本文进行了文[1]中扁长杆冲击屈曲特性的仿真计算，参考该文及后来的相关论文[2]的建模方法，采用了梁单元模型和理想边界条件(约束铰支端和自由端单元节点相应的自由度)，计算得到扁长杆的屈曲变形曲线和轴向力见图 2，屈曲变形曲线的峰/谷数少于实验测量的峰/谷数，轴向力、失稳载荷和失稳时间明显大于实验结果。图 2中还给出了文[1-2]采用有限差分法的计算结果。可见，屈曲变形曲线的计算结果优于有限元方法，但轴向力、失稳载荷和失稳时间的计算结果差于有限元方法。
产生上述误差的原因包括：实验设备的铰支端和自由端虽然经过精加工和润滑，但可能仍存在复杂的间隙和接触摩擦，不能简化为理想边界条件；此扁长杆是薄壳件，采用梁单元模型建立其有限元模型是不适当的。
2 提高冲击屈曲特性仿真分析精度的有限元模型改进方法本文在仿真模型构建中考虑铰支端和自由端的间隙配合和接触摩擦，其中铰支端结构如图 3a所示，为等间隙的轴、孔配合结构，存在转动铰轴径a和间隙b两个特征尺寸；杆的自由端结构如图 3b所示，为等间隙的半圆柱轴、槽配合结构，存在转动铰凹槽半径c和深度d两个特征尺寸。扁长杆端部的有限元网格如图 4所示。

图 3 扁长杆端部的结构和特征尺寸

图选项

图 4 扁长杆端部的有限元网格

图选项

在有限元仿真中，对扁长杆采用了多种类型的单元模型，涉及的单元类型、单元方程和积分方式如下：
Hughes-Liu梁单元^[13-14]和Belytschko-Schwer梁单元^[15-16]，采用轴向单点积分、截面全积分；
Hughes-Liu薄壳单元^[13-14]和Belytschko-Tsay薄壳单元^[15-16]，采用面内单点积分/全积分、厚度方向2/3/5点积分；
8节点厚壳单元，采用面内单点积分/减缩积分、厚度方向2/3/5点积分；
8节点实体单元和带有节点转动自由度的8节点高阶实体单元^[17]，采用单点积分/减缩积分/全积分。
所有单元均可以描述几何大变形。
在进行有限元模型计算结果数据处理时，取横截面形心处的位移作为屈曲变形，通过对整个横截面应力的积分计算出轴向力。
2.1 边界条件的参数辨识改进边界条件的参数包括端部的4个特征尺寸参数、静摩擦系数和动摩擦系数。现采用多参数优化过程辨识特征尺寸参数的合理取值，使得扁长杆的冲击弹塑性屈曲响应仿真计算结果与实验结果尽可能相符。对6个改进边界条件参数各取3个可能值，构成6因素、3水平的设计空间，以实验结果为目标，采用正交试验设计方法进行设计参数优选，结果如表 2所示，最相符的冲击屈曲响应结果如图 5所示，即b=0.010 mm薄壳单元模型。摩擦系数的辨识结果与机械设计手册^[18]中钢与钢之间有润滑的摩擦系数一致。
表 2 薄壳单元模型的改进边界条件参数

a/mm	b/mm	c/mm	d/mm	静摩擦系数	动摩擦系数
5.00	0.010	0.32	0.34	0.11	0.075

表选项

图 5 基于薄壳单元模型和改进边界条件的计算分析结果

图选项

图 5a表明，在尺寸b的变化范围内，扁长杆屈曲变形曲线的峰/谷数并不改变，峰/谷位置几乎无变化，但峰/谷值有明显变化；当间隙尺寸减小时，铰支端峰值减小，其他部分幅值增大。图 5b表明，尺寸b对扁长杆轴向力、失稳载荷和失稳时间几乎没有影响。从大量模拟计算结果可以归纳出边界条件参数(静摩擦系数除外)对动态屈曲的影响规律：在适当范围内改变模型的边界条件参数，扁长杆屈曲变形曲线峰/谷数、峰/谷位置、失稳载荷和失稳时间不变，但屈曲变形曲线的峰/谷值有明显变化；与物理直观分析一致，如果扁长杆某一端的边界条件放松对扁长杆此端运动的限制，则会使靠近此端的屈曲变形曲线的峰/谷值增大，其他部分幅值减小，反之亦然。静摩擦系数的改变对扁长杆动态屈曲几乎没有影响。
2.2 基于改进边界条件的有限单元类型及尺度优选分析通过多种单元类型、单元方程和积分方式的详细计算分析发现，对于同一类型的单元，不同单元方程和积分方式的计算模型得到的屈曲变形曲线、轴向力、失稳载荷和失稳时间基本一致。每种单元的尺度优选范围和优选结果如表 3所示，优选尺度下有限元模型的计算结果已经收敛，如图 6所示。薄壳单元模型和厚壳单元模型的计算结果完全一致，统称为壳单元模型。每种单元的改进边界条件参数均得到优化，非梁单元模型优化得到的改进边界条件参数基本一致。
表 3 单元尺度的优选范围和优选结果

单元类型	最大单元尺度/mm	最小单元尺度/mm	尺度优选结果/mm
梁单元	19.30	0.48	0.97
薄壳单元	6.43	0.24	0.48
厚壳单元	6.43	0.24	0.48
实体单元	0.32	0.11	0.16

表选项

图 6 基于改进边界条件的多种单元模型的计算分析结果

图选项

通过图 6b和6c可知，改进边界条件后，梁单元模型屈曲变形曲线峰/谷数与实验测量的峰/谷数一致；轴向力接近实验测量结果。薄壳单元、厚壳单元和实体单元模型的计算精度相当，屈曲变形曲线峰/谷位置、幅值、轴向力、失稳载荷和失稳时间的计算结果均优于梁单元模型。图 6a给出了扁长杆失稳不久的屈曲变形曲线，边界条件和单元类型对屈曲变形曲线的影响相对较小，这符合扁长杆的冲击弹塑性屈曲特性。
如图 7所示，扁长杆的等效应力沿横向分布不均匀，而且横截面发生微小的动态屈曲，而梁单元模型不能描述这种应力分布和屈曲变形，因此不适合此类问题的分析。本文将横截面的动态屈曲定义为动态反向屈曲，即扁长杆横截面屈曲方向时刻保持与同一位置处纵截面屈曲方向相反，而且横截面屈曲幅值时刻保持为同一位置处纵截面屈曲幅值的1%左右。

图 7 薄壳单元模型在t=0.75 ms时刻的屈曲变形曲线和等效应力云图(3个断面图中的横向屈曲变形放大50倍)

图选项

3 结束语本文通过改进计算模型的边界条件和单元类型，使扁长杆的冲击弹塑性屈曲特性有限元分析的各个参量均达到比较高的仿真分析精度。采用改进边界条件(考虑转动铰配合间隙、摩擦等)比采用理想边界条件的仿真结果更加准确；薄壳单元、厚壳单元和实体单元模型的仿真分析精度相当，都明显优于梁单元模型；基于改进后的计算分析模型的计算结果，揭示了动态反向屈曲现象，虽然没有实验结果对照，但符合物理变形规律。类似问题也可以通过本文模型改进技术提高仿真分析精度。
本文最终的仿真结果仍存在一定的误差，这是因为改进的边界条件仍为简化边界条件，但动态屈曲问题为高度非线性动力学问题，本文的仿真计算结果已达到了较高的精度。

参考文献

[1]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Hayashi T, Sano Y. Dynamic buckling of elastic bars (2nd report, the case of high velocity impact)[J]. Bulletin of the Japan Society of Mechanical Engineers, 1972, 15(88) : 1176–1184.DOI:10.1299/jsme1958.15.1176
[2]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Sugiura K, Mizuno E, Fukumoto Y. Dynamic instability analyses of axially impacted columns[J]. Journal of Engineering Mechanics, 1985, 111(7) : 893–908.DOI:10.1061/(ASCE)0733-9399(1985)111:7(893)
[3]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Lepik V. Dynamic buckling of elastic-plastic beams including effects of axial stress waves[J]. International Journal of Impact Engineering, 2001, 25(6) : 537–552.DOI:10.1016/S0734-743X(00)00070-1
[4]	Journal of Central South University(Science and Technology), 41(2):649-654.-->CUI Shijie, HAO Hong, Cheong H K. Theoretical study of dynamic elastic buckling of columns subjected to intermediate velocity impact loads[J]. International Journal of Mechanical Sciences, 2002, 44(4) : 687–702.DOI:10.1016/S0020-7403(02)00009-7
[5]	Journal of Central South University(Science and Technology), 41(2):649-654.-->马宏伟, 韩强, 张善元, 等. 受轴向冲击有限长弹性直杆中应力波引起的分叉问题[J]. 爆炸与冲击, 1995, 15(4) : 300–306.MA Hongwei, HAN Qiang, ZHANG Shanyuan, et al. Bifurcation problem caused by propagation of stress in finite length bar subjected to axial impact[J]. Explosion and Shock Waves, 1995, 15(4) : 300–306.(in Chinese)
[6]	Journal of Central South University(Science and Technology), 41(2):649-654.-->唐文勇, 陈国胜, 张圣坤. 含脱层复合材料层合杆在轴向应力波作用下的动态屈曲[J]. 振动与冲击, 2007, 26(2) : 14–17.TANG Wenyong, CHEN Guosheng, ZHANG Shengkun. Dynamic buckling of laminated composite bar with delamination subject to axial stress waves[J]. Journal of Vibration and Shock, 2007, 26(2) : 14–17.(in Chinese)
[7]	Journal of Central South University(Science and Technology), 41(2):649-654.-->王晓军, 王磊, 马丽红, 等. 不确定初始几何缺陷杆动态屈曲失效分析[J]. 北京航空航天大学学报, 2011, 37(12) : 1484–1489.WANG Xiaojun, WANG Lei, MA Lihong, et al. Dynamic buckling failure analysis of rod with uncertain initial geometrical imperfection[J]. Journal of Beijing University of Aeronautics and Astronautics, 2011, 37(12) : 1484–1489.(in Chinese)
[8]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Hayashi T, Sano Y. Dynamic buckling of elastic bars (1st report, the case of low velocity impact)[J]. Bulletin of the Japan Society of Mechanical Engineers, 1972, 15(88) : 1167–1175.DOI:10.1299/jsme1958.15.1167
[9]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Briseghella L, Majorana C E, Pellegrino C. Dynamic stability of elastic structures: A finite element approach[J]. Computers and Structures, 1998, 69(1) : 11–25.DOI:10.1016/S0045-7949(98)00084-4
[10]	Journal of Central South University(Science and Technology), 41(2):649-654.-->ZHANG Zheng, Taheri F. Dynamic pulsebuckling and postbuckling of composite laminated beam using higher order shear deformation theory[J]. Composites Part B: Engineering, 2003, 34(4) : 391–398.DOI:10.1016/S1359-8368(02)00134-8
[11]	Journal of Central South University(Science and Technology), 41(2):649-654.-->揭敏. 一定初缺陷杆在轴向冲击下弹塑性动态屈曲有限元计算[J]. 爆炸与冲击, 1991, 11(2) : 153–160.JIE Min. Finite element calculation on elastic-plastic dynamic buckling of a bar with finite initial imperfection under axial impact[J]. Explosion and Shock Waves, 1991, 11(2) : 153–160.(in Chinese)
[12]	Journal of Central South University(Science and Technology), 41(2):649-654.-->郑波, 王安稳. 直杆碰撞刚性壁弹塑性动力后屈曲有限元分析[J]. 爆炸与冲击, 2007, 27(2) : 126–130.ZHENG Bo, WANG Anwen. Finite element analysis for elastic-plastic dynamic postbuckling of bars subjected to axial impact[J]. Explosion and Shock Waves, 2007, 27(2) : 126–130.(in Chinese)
[13]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Hughes T J R, Liu W K. Nonlinear finite element analysis of shells: Part I three-dimensional shells[J]. Computer Methods in Applied Mechanics and Engineering, 1981, 26(3) : 331–362.DOI:10.1016/0045-7825(81)90121-3
[14]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Hughes T J R, Liu W K. Nonlinear finite element analysis of shells: Part II two-dimensional shells[J]. Computer Methods in Applied Mechanics and Engineering, 1981, 27(2) : 167–181.DOI:10.1016/0045-7825(81)90148-1
[15]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Belytschko T, Schwer L, Klein M J. Large displacement transient analysis of space frames[J]. International Journal for Numerical and Analytical Methods in Engineering, 1977, 11(1) : 65–84.DOI:10.1002/(ISSN)1097-0207
[16]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Belytschko T, Lin J I, Tsay C S. Explicit algorithms for nonlinear dynamics of shells[J]. Computer Methods in Applied Mechanics and Engineering, 1984, 42(2) : 225–251.DOI:10.1016/0045-7825(84)90026-4
[17]	Journal of Central South University(Science and Technology), 41(2):649-654.-->Yunus S M, Pawlak T P, Cook R D. Solid elements with rotational degrees of freedom: Part I hexahedron elements[J]. International Journal for Numerical Methods in Engineering, 1991, 31(3) : 573–592.DOI:10.1002/(ISSN)1097-0207
[18]	Journal of Central South University(Science and Technology), 41(2):649-654.-->闻邦椿. 机械设计手册[M].北京: 机械工业出版社, 2015.WEN Bangchun. Machine Design Handbook[M].Beijing: China Machine Press, 2015.(in Chinese)