一种新型基于共用系数的二维ANCF 高阶梁单元

周川，赵春花，独亚平，郭嘉辉

（201620 上海市 上海工程技术大学 机械与汽车工程学院）

0 引言

绝对节点坐标法（简记为ANCF）[1]是由美国伊利诺伊大学的夏巴那教授提出的一种解决多体系统大变形大转动的新型有限元方法。不同于牛顿欧拉法[2]、凯恩法[3]等其它的多体系统建模方法[4-6]，ANCF 法抛弃了传统有限元小变形假设，将转角以梯度形式表示引入到节点自由度中，得到的质量矩阵为常数阵且不存在科氏力、离心力，提出至今受到越来越多学者们的关注和研究。

研究前期，Omar 等[7]提出的二维ANCF 全参数剪切梁单元通过一般连续介质力学方法构建单元弹性力，但一般连续介质力学方法构建的单元存在一些闭锁问题[8]；沈振兴[9]和张大羽等[8]学者认为，全参数梁的位移模式在横向横纵向的插值不一致，导致无法满足广义胡克定律，产生了伪应力，使得单元计算性能下降。为解决该问题，Zhao 等[10]通过共用系数法在横向插值中引入y2，使得全参数剪切梁横纵向插值保证了2 阶次，取得了不错效果；Patel 等[11]基于独立系数引入曲率坐标，构建新的横向高阶单元，横纵向插值达到2 阶一致，有效减轻了闭锁问题，同时还提出了应变分裂法解决ANCF 梁和板单元中存在的闭锁问题；Li 等[12]在二维全参数梁基础上增加曲率自由度rxy，同时通过共用插值项系数引入了横向高阶x2y2，但是计算结果并没有很好地改善；於祖庆[13]提出了几种绝对节点坐标单元，并将新单元用于车辆钢板弹簧和轮胎的动力学建模中；徐齐平[14]提出了超弹性ANCF 单元，并将其用于建立气动软体机器人的精确动力学模型。

本文在Li 提出的单元自由度[12]基础上，通过共用系数引入横向高阶插值项，提出了一种新的梁单元。通过理论推导结合数值实例，证明了新单元的可行性，能够缓解泊松闭锁问题，提高单元收敛精度。

1 二维ANCF 横向高阶剪切梁单元

Li 等引入曲率坐标的同时，共用插值项系数，引入了横向高阶插值项-3x2y2，建立了一种新的二维ANCF 横向高阶剪切梁单元[12]，文中简记为单元LPF。本文在此基础上再一次通过共用x3系数，改为引入不同横向高阶插值项y2，构建了新的二维高阶梁单元，记为E1，其位置矢量为：

根据文献[10]，低阶单元横纵向插值不一致使得其并不满足广义胡克定律。

根据应力应变关系σ=Dε，有

不考虑横纵向高阶耦合项时有

联立式（2）和式（3），有

根据横纵向应变的定义，有

因此，新单元中通过共用系数引入的横向高阶插值项y2并能够满足广义胡克定律，可以较好地缓解泊松闭锁问题，提高收敛精度。

2 动力学方程

绝对节点坐标法描述的二维梁单元动力学方程形式如下：

式中：Me——单元质量矩阵；Ke（e）——单元刚度矩阵；Qe——单元广义外力列阵。

2.1 单元质量矩阵

式中：S——形函数矩阵；e——节点广义坐标向量。质量矩阵由单元动能推导得到

因此，单元质量矩阵为

2.2 单元刚度矩阵

一般连续介质力学方法是ANCF 方法的基本力学原理。由格林-拉格朗日应变张量可得：

其中任意一点变形梯度为

由2 阶张量的voigt 标记，将应变张量改写为应变向量形式：

这里，泊松比v 会使得横纵向正应变ε11和ε22发生耦合。

弹性力列阵：

式中：K（e）——刚度矩阵。

2.3 广义外力列阵

单元在受到外力F 作用下，单元广义外力为

3 数值实例

本节通过悬臂梁静力学、动力学等数值实例，对新单元的力学性能进行测试。

3.1 静力学大变形

悬臂梁如图1 所示，尺寸为：长l=2 m，宽w=0.1 m，高h=0.1 m。杨氏模量E=2.07×1011Pa，泊松比v=0.3，在悬臂端点重力方向施加的集中力F=500 000 N。以有限元软件ANSYS 中Beam 188单元仿真结果为参照。

图1 大变形悬臂梁Fig.1 Large deformation cantilever beam

表1 是不同单元在不同单元数划分下，悬臂梁端点处y 方向位移情况，图2 是不同单元与ANSYS 计算结果的误差情况。由表1 和图2 可知，新单元E1 与原单元LPF 相比，收敛精度有明显提高，泊松闭锁问题得到有效缓解。单元数增加时，新单元与ANSYS 计算结果误差接近于0。

表1 大变形下末端点Y 方向位移计算结果Tab.1 Calculation results of the displacement in Y direction of end point under large deformation

图2 各单元与ANSYS 在Y 方向位移误差Fig.2 Displacement error of each element in Y direction compared with ANSYS

3.3 单摆

单摆如图3 所示，长l=1 m，宽w=0.006 32 m，高h=0.252 98 m，密度ρ=5 540 kg/m3，泊松比v=0.3，杨氏模量E=700 000 Pa，重力g 为外部载荷。计算0～1.1 s 单摆的运动情况，在MATLAB 平台进行仿真计算。

图3 单摆模型Fig.3 Simple pendulum model

图4 是不同时刻各单元位姿，图5、图6 是各单元末端点不同时刻下的x、y 方向位移情况。根据图4 可知，不同时刻下，新单元的位姿与ABAQUS 软件仿真结果基本一致。由图5、图6 可知，新单元E1 末端点不同时刻下X、Y 方向的位移情况与ABAQUS 末端点的位移也基本吻合。而单元LPF 与ABAQUS有限元软件的结果存在较大误差，故新单元能够更精确地描述动力学大变形运动过程。

图4 不同时刻各单元位姿Fig.4 Poses of each element at different moments,

图5 不同时刻各单元末端点X 方向位移Fig.5 Displacement of end points of each element in X direction at different time

图6 不同时刻各单元末端点Y 方向位移Fig.6 Displacement of end points of each element in Y direction at different time

4 结论

本文提出了一种共用系数的二维ANCF 高阶剪切梁单元，通过理论推导结合静力学、动力学数值实例，证明了（1）新单元相较于单元LPF 能够有效缓解泊松闭锁问题，提高单元静力学大变形中的的收敛精度。（2）相较于单元LPF，新单元能够更加精确地进行动力学仿真和分析。