一种基于应变插值大变形梁单元的刚-柔耦合动力学分析

张志刚， 齐朝晖, 吴志刚,2

(1.大连理工大学工业装备结构分析国家重点实验室， 辽宁 大连 116024；2.大连理工大学航空航天学院， 辽宁 大连 116024)

一种基于应变插值大变形梁单元的刚-柔耦合动力学分析

张志刚1， 齐朝晖1, 吴志刚1,2

(1.大连理工大学工业装备结构分析国家重点实验室， 辽宁 大连 116024；2.大连理工大学航空航天学院， 辽宁 大连 116024)

柔性多体系统动力学中的“动力刚化”现象起因于变形间的耦合。一次近似模型成功地解决了小变形情况下的刚柔耦合建模问题，但在大变形情况下则需要考虑更多的耦合效应。选取表征梁弯曲应变的曲率和轴向应变作为单元参数进行离散；在大变形大转动基础上得到了单元两端节点运动学参数的递推关系，构造出了能够自动计及“动力刚化项”且适用于大变形刚柔耦合动力学分析的平面梁单元。最后采用所提应变插值单元求解了包含大变形和刚柔耦合动力学柔性梁的数值算例，验证了算法的正确性和有效性。

刚柔耦合； 梁单元； 应变插值； 动力刚化

引 言

梁的刚柔耦合动力学分析是近年来力学、机械及航空航天领域研究的一个热点和难点[1-12]。柔性梁的运动包含大范围刚体运动与变形运动的耦合，且描述梁变形所需的形心位移和截面角之间也存在耦合关系。这些耦合的作用机理非常复杂，给做大范围运动梁的刚柔耦合动力学的建模和求解带来了困难。

传统的柔性多体系统动力学混合坐标法直接套用了结构力学小变形假设，Kane[1]在用其对旋转柔性梁进行求解时发现高速转动时数值结果发散，表现为随着转速增加梁变得更加柔软，这与事实相违背，并首次提出了“动力刚化”的概念。这一现象引发了国内外广大学者的关注，此后大量研究围绕刚体运动与弹性运动的耦合及“动力刚化”问题展开：洪嘉振[2-6]指出，传统混合坐标法在小变形基础上忽略了纵向位移对轴向位移的影响是一种零次近似模型，低速转动忽略掉变形耦合对仿真结果影响不大，但在高速转动情况下这种近似将会产生完全错误的结果；其提出小变形情况下在轴向位移表达式中补充纵向位移二阶项的所谓一次近似模型，用该模型进行了中心刚体旋转柔性梁及末端带集中质量的旋转柔性梁动力学特性，并进一步通过实验验证了一次近似模型能够捕捉“动力刚化项”。章定国[7]在一次近似模型基础上对平面运动柔性梁进行了进一步研究，指出大范围运动与弹性运动的耦合不仅会产生“动力刚化”现象，还会产生“动力柔化”现象，并对旋转内接悬臂梁的“动力柔化”效应进行了分析。和兴锁等[8]在Shi P等[9]研究基础上提出一次近似模型未能充分考虑几何非线性变形对横向变形的影响，因此在纵向与横向位移中增加了新的变形二次耦合项。

上述研究大都局限于单元本身为小变形情况，针对包含大变形的柔性梁刚柔耦合问题，刘锦阳[10]采用绝对坐标法对旋转柔性梁进行了研究，并与一次近似模型进行了对比，发现随着梁刚度降低一次近似模型结果发散，而绝对坐标梁单元能够得到收敛的结果。陈思佳[11]在一次近似模型基础上采用模态坐标对梁进行离散，通过保留了动力学中方程中的变形耦合量相关的高阶项，对考虑大变形的平面梁刚柔耦合动力学分析做了初步研究，但是通过模态坐标离散的求解精度和效率取决于模态集的选取方式。

现有研究大都选取节点位移作为有限元离散参量，在小变形情况下根据伯努利梁变形耦合关系添加位移高阶耦合项来构造了一次近似模型或者高次近似模型。而直接构造满足大变形大转动的变形耦合位移梁单元仍然十分困难和复杂。本文在对曲率插值梁单元研究[12]的基础上，选取梁的弯曲应变(曲率)和轴向应变作为单元参数进行离散，从而得到了自动计及位移耦合关系的梁单元。采用本单元对平面梁大范围运动刚柔耦合问题进行了研究，所得结果与相关工作[13-16]吻合，从而验证了所构造单元的正确性和有效性。

1 大范围运动梁单元虚功率方程

1.1 梁单元运动学描述

平面梁单元如图1所示，整体惯性坐标系为{g1g2}，在梁左端截面建立固结在截面上的单元坐标系{e1e2}，第一根轴e1指向梁长方向，初始时刻与整体坐标系第一根轴夹角φ0。在每个梁截面建立随截面转动的截面坐标系{t1t2}，变形前截面坐标系与单元坐标系方向保持一致。

图1 大范围运动平面梁单元Fig.1 A plane beam element undergoes large range of motion

如图1所示，弧长s处截面形心矢径为

(1)

式中r1为单元坐标系原点相对整体坐标系原点矢径，ρ为弧长s处截面形心相对于单元坐标系原点矢径。取出原始弧长为ds的一段梁微元，变形后几何关系

(2)

(3)

其中 (·)′=d(·)/ds，θ为截面相对于单元坐标系的转角。

单元左端截面相对于整体坐标系的转角φ1表征了单元大范围刚体转动，截面相对于单元坐标系的转角θ代表了单元的弹性变形，因此单元域内截面相对于整体坐标系的转角为

(4)

1.2 梁单元应变场的离散

对于不超过弹性范围的平面梁，其本构关系可用截面坐标系中力的分量描述为

(5)

式中F为截面轴向力；M为弯矩；EA,EIz分别为截面拉压刚度和弯曲刚度。梁的应变虚功率

(6)

由上式可知，梁单元应变虚功率只与曲率和轴向应变有关，而不受单元位移和转动大小的影响，因此选取梁曲率和轴向应变进行离散可以方便地得到单元应变虚功率和简洁的单元刚度矩阵。

(7)

其中

(8)

其中ξ=s/l，梁曲率速率虚变分由式(7)可得

(9)

将式(7)，(9)及轴向应变ε是常量的假设代入单元应变虚功率式(6)

(10)

式中K为单元刚度阵

(11)

1.3 单元内运动学递推

单元域内截面的转角可由曲率插值函数式(7)积分得到

(12)

其中

(13)

转角θ为截面相对于单元转角，它完全由单元的弯曲变形确定，将其代入式(12)可得截面相对整体坐标系转角为

(14)

由梁截面形心坐标弧长导数式(3)，对弧长s积分并注意到轴向应变ε为常量，可以得到截面形心坐标

(15)

从这两个式子可以看出截面形心坐标与转角和轴向应变高度耦合。需要说明的是传统混合坐标零次近似模型及一次近似模型都假设转角为小量，即：θ≪0，从而将三角函数近似为：sinθ≈θ，cosθ≈1-θ2/2，并由式(3)得

(16)

将其带入式(15)得

(17)

上式推导过程中忽略掉了转角θ的二次以上项及θ和轴向应变ε的耦合项，式中加下划线部分为一次近似模型在零次模型基础上增加的横向位移对轴线变形的二次耦合项。由此可见混合坐标一次近似模型不适用于包含大转动的大变形情况。

本文不作小转角假设，基于此推导的单元完整地保留了变形间的耦合关系，不仅适用于小转角情况，同样也能满足大转角情况。由式(12)可知截面转角θ为弧坐标s的二次式，因此式(15)两个积分不存在有理形式解，本文由高斯积分得到数值解。

直接由式(15)就可得到截面形心坐标场，但是在推导质量阵和广义力时将会出现二重数值积分，从而会引起公式推导和数值求解上的不便。再次考察单元坐标系下截面形心坐标x,w在端部满足

(18)

(19)

式中

(20)

(21)

弧长s处截面形心坐标为

(22)

其中单元坐标系与整体坐标系转换矩阵

(23)

式中φ0为单元初始方位角。

对式(22)求时间导数得到截面形心速度

(24)

(25)

根据单元坐标系下截面形心坐标式(19)，可得

(26)

其中

(27)

(28)

将式(26)代入得到截面形心速度

(29)

对速度取变分

(30)

将式(26)代入式(25)得到截面形心加速度

(31)

其中加速度余量

(32)

对截面转角式(14)求导可得截面转动角速度及虚角速度

(33)

角加速度

(34)

为了方便推导，本文将整体坐标系坐标下截面形心坐标r和转角φ组成列阵为q，两端面运动学参数

(35)

由式(14)，(29)，(31)，(33)，(34)，单元两端截面运动学递推关系可以用矩阵表示为

(36)

其中

(37)

(38)

1.4 梁单元惯性力虚功率

根据伯努利梁理论，平面梁可以看作由无数做平面运动的刚性截面组成，其惯性力虚功率为

(39)

其中ρ代表梁的线密度，Iz代表梁截面惯性矩。

将式(30)～(34)代入上式

(40)

质量矩阵和广义力为

(41)

(42)

其中

(43)

(44)

质量阵和广义力中常系数矩阵

(45)

1.5 梁单元外力虚功率

梁单元上作用的外力虚功率为

(46)

式中p(s)为分布载荷；fi,mi分别为单元节点处所受集中力和力矩。将式(30),(33)代入上式

(47)

其中

(48)

(49)

对于考虑重力的情况，可将重力看做一种特殊的分布力，其对应广义力为

(50)

(51)

另外，还可以将结构阻尼当做外力处理，本文采用文献[13]阻尼模型，单元结构阻尼力虚功率为

(52)

式中η为阻尼系数。对比单元应变虚功率式(10)，可将机构阻尼力虚功率表示为

(53)

其中广义结构阻尼力列阵为

(54)

2 平面梁刚柔耦合动力学方程

2.1 单元间运动学递推

(1)单元坐标系原点选在靠近起始点一端节点，第一根轴方向背离起始节点指向第二个节点;

(2)沿背离起始节点方向单元编号递增。

单元连接拓扑关系通过定义邻接单元数组L描述：L为1×n行向量，其中L(i)代表与第i号单元相邻并靠近起始节点1的那一侧单元编号；规定与起始节点相连单元的邻接单元为0。

如图2所示，下面单元标号满足规则标号，箭头标出了单元坐标系原点所在节点和第一根坐标轴方向。

图2 梁单元规则标号Fig.2 Regularly labeling of beam elements

按照上述约定，图示由三个单元组成的梁的邻接单元数组为：L=[0 1 2]。

利用邻接单元数组L，根据单元内运动学递推关系式(36)，可将i号单元坐标系所在节点运动参量与起始节点关系表示为

(55)

其中

(56)

(57)

2.2 梁系统刚柔耦合动力学方程

弹性体虚功率方程表述为：弹性体所受外力虚功率等于惯性力虚功率和应变虚功率之和。划分为n个单元柔性梁的虚功率方程为

(58)

将梁单元应变虚功率式(10)、惯性力虚功率式(40)、外力虚功率式(47)及单元间运动学递推关系式(55)代入上式得到

其中质量矩阵

(60)

广义力列阵

(61)

由梁系统虚功率方程式(59)中虚变分独立，其系数为零可得平面梁动力学方程

(62)

3 数值算例

在这一节中采用本文单元编写了动力学分析程序，分别选取了悬臂梁纯弯曲、大范围运动已知的旋转柔性梁及包含大变形的自由下落柔性单摆3个算例进行了求解，并将所得结果与相关文献[10-11,14-15]进行了分析比较。

3.1 悬臂梁受梁端集中弯矩作用

悬臂梁端受集中弯矩M作用，梁长度L=10 m，梁端线位移u，v如图3所示。

图3 悬臂梁受集中弯矩作用Fig.3 Cantilever beam undergoes bending load

由于存在结构阻尼，梁最终会处于静平衡状态，这里采用本文算法将整个梁划分为5个单元，为了减少悬臂梁达到平衡状态所用时间，这里选取阻尼比η=0.1进行仿真。在不同弯矩作用下梁最终的变形状态如图4所示。

图4 不同弯矩作用下悬臂梁的大变形Fig.4 Large deformation of the cantilever under bending moments

下面给出悬臂梁在弯矩M=2πEIz/L作用下梁末端位移随时间变化曲线如图5所示。

图5 悬臂梁末端受集中弯矩的位移Fig.5 The tip displacements of the cantilever with an end bending load

从位移时间曲线图可以看出，在较大结构阻尼作用下悬臂梁经过两个周期就达到了最终的平衡状态，且仅用5个单元就得到了与解析解高度吻合的结果，由此表明本文单元在处理包含大变形问题的优越性。

3.2 旋转柔性梁大范围运动已知

大范围运动已知的旋转柔性梁已被很多学者[14-15]用来考察所提方法补充“动力刚化项”的能力。如图6所示，柔性梁物理参数选取与文献[14]相同：梁长L=10 m截面面积A=1 m2，截面惯性矩Iz=5×10-4m4，弹性模量E=2.8×107Pa，密度ρ=1.2 kg/m3。

图6 旋转柔性梁Fig.6 A flexible beam on a rotating base

柔性梁绕转动中心的转角φ(t)规律为

(63)

这里忽略梁结构阻尼，将梁划分为5个单元，仿真时间为30s，得到梁自由端位移响应曲线如图7所示。

图7 自由端位移曲线Fig.7 Displacement curves at the free end

从图7可以看出：梁端部位移随着转动速度增加而增大，当转速增加到一定数值维持不变，轴向位移逐渐回落并稳定在5.14×10-4m，横向位移在0附近小幅度波动。采用本文单元的仿真结果与文献[14-15]完全相符，可以验证算法对动力刚化问题的适用性。

3.3 自由下落柔性单摆

最后一个算例研究柔性单摆在重力作用下自由下落问题，单摆初始时刻水平放置如图8所示。

图8 柔性单摆Fig.8 The flexible pendulum

这个算例的大范围刚体运动是未知的，是个典型的大范围运动与变形运动耦合的问题。针对这个算例刘锦阳[10]和陈思佳[11]分别采用不同建模方法进行仿真求解，并与混合坐标法一次近似模型进行了对比。通过选取较低的材料弹性模量使得弹性变形不再是小量，以此来验证方法对大变形问题的适用性。柔性摆的长度L=1.8 m，E=6.895×109Pa，A=2.5 cm2，Iz=0.13 cm4，密度ρ=2 766.67 kg/m3。

这个算例中将单摆划分为5个单元，以本文方法和绝对节点坐标方法[10]在Matlab环境中选用刚性微分方程组求解器ode23tb，数值积分参数设定为：绝对误差限AbsT=0.1，相对误差限RelT=0.01，在同一台微机上进行仿真2.5 s，所用CPU时间如表1所示。

表1 绝对坐标方法与本文方法CPU时间对比

Tab.1 The CPU time comparison between absolute coordinate method and the proposed method

方法文献[10]绝对坐标法本文方法CPU时间225．05s6．29s

两种方法所得单摆末端位移时间曲线如图9所示。

图9 柔性单摆的末端位移Fig.9 The tip deformation of the flexible pendulum

根据文献[10-11]结果，混合坐标法一次近似模型仿真结果发散，表明其不适用与包含大变形的梁的刚柔耦合动力学分析。而本文采用应变插值所得结果与文献[10-11]相吻合，具有较好的计算精度。

从两种方法所用CPU时间对比可以看出，本文单元求解效率优于文献[10]绝对节点坐标单元，也具有很好的计算效率。

4 结 论

本文通过直接选取单元应变进行插值构造了一种不受单元变形大小限制，适用于刚柔耦合动力学分析的平面梁单元。梁变形的位移及转角之间的耦合关系通过几何方程积分得到，并由此给出了节点运动学递推关系。本文所提梁单元有如下特点：

(1)选取单元应变进行离散不受单元刚体运动的影响，得到了简洁的单元刚度矩阵及节点力；

(2)根据应变场由几何方程积分得到梁截面形心坐标和转角耦合关系，自动计及了“动力刚化项”；

(3)单元不受变形大小的限制，适用于包含大变形梁的刚柔耦合动力学分析。

[1] Kane T R, Ryan R, Banerjee A K. Dynamics of a cantilever beam attached to a moving base[J]. Journal of Guidance, Control, and Dynamics, 1987, 10(2): 139—151.

[2] 杨辉, 洪嘉振, 余征跃. 动力刚化问题的实验研究[J]. 力学学报, 2004, 36(1): 118—124.

YANG Hui, HONG Jiazhen, YU Zhengyue. Experimental investigation on dynamic stiffening phenomenon[J]. Acta Mechanica Sinica, 2004，36(1)：119—124.

[3] 杨辉, 洪嘉振, 余征跃. 刚柔耦合建模理论的实验验证[J]. 力学学报, 2003, 35(2): 253—256.

YANG Hui, HONG Jia-zhen，YU Zheng-yue．Experiment validation of modeling theory for rigid-flexible coupling systems[J]．Acta Mechanica Sinica，2003，35(2)：253—256．

[4] Cai G P, Hong J Z, Yang S X. Dynamic analysis of a flexible hub-beam system with tip mass[J]. Mechanics Research Communications, 2005, 32(2): 173—190.

[5] Cai G P, Lim C W. Dynamics studies of a flexible hub-beam system with significant damping effect[J]. Journal of Sound and Vibration, 2008, 318(1): 1—17.

[6] 刘锦阳, 洪嘉振. 柔性梁的刚-柔耦合动力学特性研究[J]. 振动工程学报, 2002, 15(2): 194—198.

Liu Jinyang, Hong Jiazhen. Study on rigid-flexible coupling dynamic behaviour of flexible beam[J]. Journal of Vibration Engineering,2002, 15(2): 194—198.

[7] 方建士, 章定国. 旋转内接悬臂梁的刚柔耦合动力学特性分析[J]. 物理学报, 2013，62(4): 305—311.

Fang Jianshi, Zhang Dingguo. Analyses of rigid-flexible coupling dynamic properties of a rotating internal cantilever beam[J]. Acta Physica Sinica,2013,62(4) : 305—311.

[8] 邓峰岩, 和兴锁, 杨永锋, 等. 计及非线性变形的刚-柔耦合动力学建模[J]. 机械强度, 2006, 28(6): 800—804.

Deng Fengyan, He Xingsuo, Yang Yongfeng, et al. Dynamics modeling for a rigid-flexible coupling system with onlinear deformation field[J]. Journal of Mechanical Strength， 2006, 28(6): 800—804

[9] Shi P, McPhee J, Heppler G R. A deformation field for Euler-Bernoulli beams with applications to flexible multibody dynamics[J]. Multibody System Dynamics, 2001, 5(1): 79—104.

[10]李彬, 刘锦阳. 大变形柔性梁系统的绝对坐标方法[J]. 上海交通大学学报, 2005, 39(5): 827—831.

Li Bin, Liu Jinyang. Application of absolute nodal coordination formulation in flexible beams with large deformation[J]. Journal of Shanghai Jiaotong University, 2005, 39(5):827—831．

[11]陈思佳, 章定国, 洪嘉振. 大变形旋转柔性梁的一种高次刚柔耦合动力学模型[J]. 力学学报, 2013, 45(2): 251—256.

Chen Sijia, Zhang Dingguo, Hong Jiazhen. A high-order rigid-flexible coupling model of a rotating flexible beam under lagre deformation[J]. Acta Mechanica Sinica, 2013, 45(2): 251—256.

[12]张志刚, 齐朝晖, 吴志刚. 基于曲率插值的大变形梁单元[J]. 应用数学和力学, 2013 34(6): 620—629.

Zhang Zhigang, Qi Zhaohui, Wu Zhigang. A large deformation beam element based on curvature interpolation[J]. Applied Mathematics and Mechanics, 2013, 34(6): 620—629.

[13]Garcíd D, Valverde J, Dominguez J. An internal damping model for the absolute nodal coordinate formulation[J]. Nonlinear Dynamics, 2005, 42(4): 347—369.

[14]Simo J C, Vu-Quoc L. On the dynamics in space of rods undergoing large motions—a geometrically exact approach[J]. Computer Methods in Applied Mechanics and Engineering, 1988, 66(2): 125—161.

[15]Zhao Z, Ren G. A quaternion-based formulation of Euler-Bernoulli beam without singularity[J]. Nonlinear Dynamics, 2012, 67(3):1 825—1 835.

Rigid-flexible dynamics analysis of a large deformation beam element based on interpolation of strains

ZHANGZhi-gang1,QIZhao-hui1,WUZhi-gang1,2

(1.State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China; 2.School of Aeronautics and Astronautics, Dalian University of Technology, Dalian 116024, China)

The “dynamic stiffening” phenomenon in flexible multi-body system dynamics is due to the deformation coupling. The first-order approximation model has been successfully applied in the rigid-flexible coupling modeling of small deformation. However, it is found necessary to consider more deformation coupling effects in lager deformation cases. In this paper, the beam bending curvature and axial strain are selected as the element parameters. Then the recursion formulations are obtained for the kinematic parameters of two end nodes of an element based on theories of large deformation and finite rotation. A planar beam element used for large deformation rigid-flexible dynamics analysis is proposed, which can automatically take into account the “dynamic stiffening terms”. Finally, the validity and effectiveness of the proposed algorithm are verified through some numerical examples which involve the large deformations and rigid-flexible dynamics of beams.

rigid-flexible coupling; beam element; strain interpolation; dynamic stiffening

2013-12-26；

2014-06-12

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

O313.7

A

1004-4523(2015)03-0337-08

10.16385/j.cnki.issn.1004-4523.2015.03.001

张志刚(1984—),男，博士研究生。电话： 15326175369； E-mail:zhigangzhang@foxmail.com

齐朝晖(1964—),男，教授，博士生导师。E-mail: zhaohuiq@dlut.edu.cn