文章快速检索 高级检索
0

页岩气专题论文

引用本文 [复制中英文]

张越, 赵阳, 谭春林, 刘永健. ANCF索梁单元应变耦合问题与模型解耦1)[J]. 力学学报, 2016, 48(6): 1406-1415. DOI: 10.6052/0459-1879-16-127.
[复制中文]
Zhang Yue, Zhao Yang, Tan Chunlin, Liu Yongjian. THE STRAIN COUPLING PROBLEM AND MODEL DECOUPLING OF ANCF CABLE/BEAM ELEMENT1)[J]. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(6): 1406-1415. DOI: 10.6052/0459-1879-16-127.
[复制英文]

基金项目

1) 国家重点基础研究发展计划(973计划)(2013CB733004)和微小型航天器技术国防重点学科实验室开放基金(HIT.KLOF.MST.201508)资助项目

通讯作者

2) 赵阳, 教授, 主要研究方向:复杂航天器系统动力学.E-mail:yangzhao@hit.edu.cn

文章历史

2016-05-13 收稿
2016-08-24 录用
2016-08-31网络版发表
ANCF索梁单元应变耦合问题与模型解耦1)
张越*, 赵阳*2), 谭春林, 刘永健     
*. 哈尔滨工业大学, 哈尔滨 150001
†. 中国空间技术研究院, 北京 100086
摘要: 索梁结构在土木工程、航空航天等领域有着广泛的应用.在各类索梁动力学建模方法中,由于绝对节点坐标方法(absolute nodal coordinate formulation,ANCF)能够描述柔性体的大变形和大转动问题,因此非常适合大变形索梁结构的动力学建模.对绝对节点坐标索梁单元的应变进行分析可知,弯曲变形会引起单元内部轴向应变的不均匀分布,即单元轴向应变与弯曲应变相互耦合.这种应变耦合效应使单元产生伪应变能,导致单元刚度增大,造成单元失真.分析不同弯曲角下的单元应变及应变能可知,弯曲变形越大,单元失真越严重.通过构造等效一维杆单元重新描述轴向应变,实现了轴向应变与弯曲应变解耦.在此基础上推导广义弹性力,得到了绝对节点坐标索梁单元的应变解耦模型.对解耦前后的两种梁模型进行静力学和动力学仿真,结果表明:解耦模型消除了单元伪应变,相比原模型表现出更好的收敛性和曲率连续性,在相同单元数目下具有更高的精度.同时由于解耦模型降低了单元刚度,因此相比原模型,速度曲线中不再有高频振动.
关键词: 绝对节点坐标            应变耦合    曲率连续性    
引言

索梁结构在土木工程、航空航天等领域有着广泛的应用,其结构参数设计、系统任务分析等均依赖于精确的索梁动力学模型.对于大变形索梁结构,基于小变形、小转动的传统建模方法已不再适用[1-2].为此,Shabana[3-4]提出了适用于大变形柔性体动力学分析的绝对节点坐标方法(absolute nodal coordinate formulation,ANCF),并在文献[5]中给出了绝对节点方法的明确定义.该方法通过单元节点绝对位置和物质导数描述柔性体运动和变形,无小变形假设,适合于大变形索梁的动力学建模,同时也能够作为其他大变形索梁动力学模型的验证依据[6-7].

在绝对节点坐标索梁单元建模方面,Shabana等[8]提出了考虑剪切和扭转的三维梁单元建模理论,并通过数值仿真对二节点和四节点梁单元模型进行了验证[9]. Gerstmayr等[10]提出了一种适用于细长索梁的梯度缺省梁单元,该单元基于Euler-Bernoulli梁假设,通过梁中心轴线描述位移场,数值仿真结果表明该单元具有很好的收敛性和较高的计算效率.朱大鹏等[11]在Dombrowski[12]的研究基础上提出大变形索梁单元,并给出多种数值和工程应用算例. Liu等[13]利用Green-Lagrange应变张量推导曲梁单元应变能,得到了空间曲梁单元模型. Dmitrochenko等[14]集中给出了多种不同维度和节点坐标数的梁单元模型,这些单元的位移场全部采用多项式进行描述.目前这些单元已广泛应用于滑轮柔索[15-18]、铁路接触网[19-20]、柔性机械臂[21]等含索梁系统的动力学建模与分析中.

上述ANCF索梁单元的位移场均表示为关于物质坐标的多项式参数函数.文献[22]对参数形式的多项式空间曲线进行研究,得知空间曲线的弯曲程度会影响参数点在曲线上的分布.若单元位移场采用此多项式描述,则意味着轴向应变与弯曲应变相互耦合.事实上,单元的位移场与应变场均完全取决于多项式系数(单元广义坐标).若单元不同的应变分量存在相同的依赖项,就会导致不同应变分量相互耦合,这种应变耦合效应会使单元在变形过程中产生伪应力,从而导致模型失真,这也是造成ANCF梁单元闭锁问题[23-25]的主要原因.针对ANCF梁单元的闭锁问题,研究者通常采用弹性力线性化[23]、选择降阶积分[26-27]等方法来削弱闭锁效应的影响.然而,这些方法并未直接解决梁单元的应变耦合问题,为此本文以文献[10]中的细长索梁单元为例,对索梁单元的应变耦合效应进行分析,在此基础上研究应变解耦方法,以提高单元性能.

细长索梁单元的轴向应变和弯曲应变都依赖于位移场的全部多项式系数,因此两应变不独立,存在相互耦合的问题.本文首先详细分析索梁单元应变耦合机理和现象,并从能量角度分析应变耦合效应对单元精度的影响.其次研究索梁单元的应变解耦方法,推导索梁单元的应变解耦模型.最后通过静力学和动力学数值仿真,对解耦前后两种梁模型的收敛性、应变连续性以及动态特性进行深入分析,以验证解耦方法的有效性.

1 ANCF索梁单元应变耦合分析 1.1 ANCF索梁单元模型

本节给出三维二节点索梁单元模型(如图 1),该单元适用于细长索和梁的动力学建模.为方便叙述,后文将该索梁单元简称为梁单元.

图 1 ANCF梁单元模型 Figure 1 Model of ANCF beam element

梁中心轴线上任意点的位置坐标可表示为关于物质坐标$x$的三次多项式

$ {\boldsymbol r} = \left[\!\! \begin{array}{c} {r_X } \\ {r_Y } \\ {r_Z } \end{array} \!\! \right]=\left[\!\! \begin{array}{l} {a_0 + a_1 x + a_2 x^2 + a_3 x^3} \\ {b_0 + b_1 x + b_2 x^2 + b_3 x^3} \\ {c_0 + c_1 x + c_2 x^2 + c_3 x^3} \end{array} \!\! \right] $ (1)

上式为梁单元的位移场函数,可通过对单元节点坐标的三次Hermite插值得到

$ {\boldsymbol r}(x, t) = {\boldsymbol S}(x){\boldsymbol q}(t) $ (2)

其中${\boldsymbol q}(t)$${\boldsymbol S}(x)$分别为插值节点和插值函数,即单元节点坐标和形函数.对于给定单元,其位形仅取决于单元节点坐标.单元节点坐标由节点位置及其对物质坐标的导数组成,对于长度为$L$的梁单元,其节点坐标可表示为

$ {\boldsymbol q} = \left[\begin{array}{cc} {\boldsymbol q}_1^{\rm T} & {\boldsymbol q}_2^{\rm T} \end{array} \right]^{\rm T} = \left[{{\boldsymbol r}^{\rm T}(0)\quad {\boldsymbol r}_x^{\rm T} (0)\quad {\boldsymbol r}^{\rm T}(L)\quad {\boldsymbol r}_x^{\rm T} (L)} \right]^{\rm T} $ (3)

单元形函数表达如下

$ \left. {{\boldsymbol S}(x) = \left[\!\! \begin{array}{llll} S_1 (x){\boldsymbol I}_3 & {S_2 (x){\boldsymbol I}_3 } & {S_3 (x){\boldsymbol I}_3 } & {S_4 (x){\boldsymbol I}_3 } \end{array}\!\!\right] \\ S_1 (x) = 1 - 3\xi ^2 + 2\xi ^3, \quad S_2 (x) = L(\xi - 2\xi ^2 + \xi ^3) \\ S_3 (x) = 3\xi ^2 - 2\xi ^3, \quad S_4 (x) = L( - \xi ^2 + \xi ^3)} \right\} $ (4)

其中,$\xi =x / L$${\boldsymbol I}_3 $为三阶单位矩阵.

由单元的位移场函数可知,该单元考虑了轴向和弯曲变形,通过轴向应变$\varepsilon$和曲率$\kappa $可描述梁的变形程度,具体表达式如下[10]

$ \varepsilon = \dfrac{1}{2}\left( {{\boldsymbol r}_x ^{\rm T}{\boldsymbol r}_x - 1} \right) $ (5)
$ \kappa = \dfrac{ | {\boldsymbol r}_x \times {\boldsymbol r}_{xx} | }{\left| {{\boldsymbol r}_x } \right|^3} \quad \ \ $ (6)

基于Euler-$\!$-Bernoulli梁理论,该梁单元的应变能为

$ U = U_1 + U_2 = \dfrac{1}{2}\int_0^L EA\varepsilon ^2 {\rm{d}}x + \dfrac{1}{2}\int_0^L EJ\kappa ^2 {\rm{d}}x $ (7)

其中,$U_1 $$U_2 $分别表示轴向应变能和弯曲应变能,$E$为弹性模量,$A$$J$分别为截面积和截面惯性矩. 将式(7)对广义坐标求偏导数,可得单元的广义弹性力为

$ {\boldsymbol Q} = {\boldsymbol Q}_1 + {\boldsymbol Q}_2 = \int_0^L EA\varepsilon \left( {\dfrac{\partial \varepsilon }{\partial {\boldsymbol q}} } \right)^{\rm T} {\rm{d}}x + \int_0^L EJ\kappa \left( {\dfrac{\partial \kappa }{\partial {\boldsymbol q }}} \right)^{\rm T} {\rm{d}}x $ (8)
1.2 应变耦合分析

ANCF梁单元的位移场函数可通过对节点坐标进行三次Hermite插值得到,式(5)中轴向应变和式(6)中曲率均完全取决于12个节点坐标分量的值,因此梁单元轴向应变与弯曲应变相互耦合.本文通过分析不同弯曲程度下单元的轴向应变,说明ANCF梁单元的应变耦合效应及其对单元精度的影响.在进行耦合分析前,定义如下变量:

定义梁单元整体轴向应变为

$ \varepsilon _L = \int_0^L {\left| {{\boldsymbol r}_x } \right| {\rm{d}}x} / L - 1 $ (9)

其中$\int_0^L {\left| {{\boldsymbol r}_x } \right| {\rm{d}}x}$为梁单元变形后长度,因此$\varepsilon _L $能够反映单元长度的变化.

定义单元弯曲角为

$ \theta = {\rm arcos} \dfrac{{\boldsymbol r}_x \left( 0 \right)^{\rm T}{\boldsymbol r}_x \left( L \right)}{\left| {{\boldsymbol r}_x \left( 0 \right)} \right| \cdot \left| {{\boldsymbol r}_x \left( L \right)} \right|} $ (10)

由上式可知,单元弯曲角为单元两端节点处切线方向的夹角,能够反映单元弯曲变形程度.

对长度为$L=1$ m的梁单元进行分析,令节点处的轴向应变为0 ($\left| {{\boldsymbol r}_x(0)} \right| = \left| {{\boldsymbol r}_x (L)} \right| = 1)$,且变形后单元长度不变($\varepsilon _L = 0)$. 在0$^\circ$$\sim $180$^\circ$之间均匀选取7组不同的弯曲角,得到不同弯曲角下单元位形曲线和轴向应变曲线,如图 2图 3所示.

图 2 位形曲线 Figure 2 Configuration curves
图 3 轴向应变曲线 Figure 3 Axial strain curves

当单元不发生弯曲变形时,单元上任意点的轴向应变为0.但当单元发生弯曲变形时,单元内部产生了分布不均的轴向应变,且随着弯曲角的增大而增大,当弯曲角为180$^\circ$时,轴向应变峰值达到0.13.

图 3可知,单元总长度不变,但单元中间部分伸长,而两端部分受到压缩,单元内部会产生一定的轴向应变能,这种应变能并非由外界载荷引起,而是由弯曲变形诱发产生的伪应变能(如图 4),会降低模型的精确性.

图 4 应变耦合示意图 Figure 4 Diagram of strain coupling

下面分析应变耦合效应对单元应变能的影响.以半径0.02 m,长度1 m,弹性模量为2$\times $10$^{7}$ Pa的圆截面梁单元为例,分析应变能随弯曲角的变化情况,如图 5所示.

图 5 不同弯曲角下的应变能 Figure 5 Strain energy in different bending angles

为了更清晰地体现轴向应变能和弯曲应变能的变化,图中对弯曲角小于120$^\circ$的部分进行了放大. 随着弯曲角的增大,轴向应变能和弯曲应变能均不断增大. 其中弯曲应变能变化较为均匀,而轴向应变能的增速越来越快,在弯曲角小于120$^\circ$时轴向应变能小于弯曲应变能,而当弯曲角达到180$^\circ$时,轴向应变能已远大于弯曲应变能.

表 1给出了不同弯曲程度下单元应变能的具体值. $U_{1}/U_{2}$为轴向应变能占弯曲应变能的百分比,其值随着弯曲变形程度增大而迅速增大.弯曲变形引起的伪轴向应变能会使单元产生刚化.例如对于图 2中的纯弯曲梁,理论上其应变能应为$U = \dfrac{1}{2}\int_0^L {EJ\kappa ^2 {\rm{d}}x} =\dfrac{1}{2}{\boldsymbol q}^{\rm T}{\boldsymbol K}_t {\boldsymbol q}$,而对于实际的ANCF梁单元,由于引入了伪轴向应变能,此时其应变能为$U = \dfrac{1}{2}\int_0^L {EJ\kappa ^2{\rm{d}}x + \dfrac{1}{2}\int_0^L {EA\varepsilon ^2 {\rm{d}}x} } =\dfrac{1}{2}{\boldsymbol q}^{\rm T}{\boldsymbol K \boldsymbol q}$,伪应变能的存在使得ANCF梁单元的刚度${\boldsymbol K}$大 于理论值${\boldsymbol K}_{\rm t}$,从而引起单元刚化. 单元弯曲变形越大,刚化效应越显著,单元失真越严重.

表 1 应变能数据分析 Table 1 Analysis of strain energy data

ANCF梁单元常用于空间飞网[28-29]、空间机械臂末端捕获绳索[30]、绳系卫星[31]等含柔索系统的动力学分析中.对于这些弯曲刚度小、变形大的柔索结构,应变耦合效应引起的伪轴向应变能占总应变能的比重更大,因此对模型精度有更大的影响.

2 梁单元的应变解耦模型

上一节从能量的角度分析了应变耦合效应对单元精度的影响,得知应变耦合效应会使单元刚化,这种刚化现象对大变形索梁模型精度的影响是不可忽视的.通常可通过加密网格削弱耦合效应的影响,但这也势必会降低动力学仿真计算效率,且没有从根本上解决应变耦合问题.

2.1 应变解耦

ANCF梁单元的应变耦合效应主要体现在受弯曲变形的影响,单元会产生分布不均匀的轴向应变.因此为了实现轴向应变与弯曲应变解耦,本文构造了等效一维杆单元,利用该单元重新描述梁的轴向应变,使单元轴向应变独立于弯曲应变.

新构造的一维杆单元是通过对梁单元进行等效得到的,如图 6所示.该单元只存在轴向应变,无弯曲应变,其节点坐标定义为

图 6 等效一维杆单元 Figure 6 Equivalent one-dimensional rod element
$ {\boldsymbol q }_\Delta \left( t \right) = \left[ {r^\Delta \left( 0 \right)} \ \ {r_x^\Delta \left( 0 \right)} \ \ {r^\Delta \left( L \right)} {r_x^\Delta \left( L \right)} \right]^{\rm T} $ (11)

式中,$r^\Delta (0)$$r^\Delta (L)$为单元两端节点位置坐标,$r_x^\Delta (0)$$r_x^\Delta (L)$为单元两端节点位置坐标对物质坐标的导数,能够反映杆单元节点处的轴向应变.

任意时刻,杆单元与梁单元有一一对应的关系,且满足如下等效准则:

(1)杆单元原长及变形后长度均与梁单元相同.

(2)杆单元两端节点处轴向应变与梁单元相同.

根据准则(1),有$r^\Delta (L) - r^\Delta (0) = L^\Delta $,其中$L^\Delta = \int_0^L {\left| {{\boldsymbol r}_x } \right|{\rm{d}}x} $为梁单元变形后长度. $r^\Delta (0)$$r^\Delta (L)$的具体取值仅影响杆单元的空间位置,不影响其轴向应变大小,因此令$r^\Delta (0) = 0$,则$r^\Delta (L) = L^\Delta $. 根据准则(2),有$r_x^\Delta (0) = \left| {{\boldsymbol r}_x (0)} \right|$$r_x^\Delta (L) = \left| {{\boldsymbol r}_x (L)} \right|$,此时杆单元和梁单元在节点处具有大小相同的轴向应变. 根据以上分析,杆单元节点坐标最终可表示为

$ {\boldsymbol q}_\Delta \left( t \right) = \left[0 {\left| {{\boldsymbol r}_x \left( 0 \right)} \right|} \ \ {\int_0^L {\left| {{\boldsymbol r}_x } \right| {\rm{d}}x} } {\left| {{\boldsymbol r}_x \left( L \right)} \right|} \right]^{\rm T} $ (12)

杆单元的位移场维度为1,仍可表达为形函数与广义坐标乘积的形式

$ r^\Delta (x, t) = {\boldsymbol S }^\Delta (x){\boldsymbol q }_\Delta (t) $ (13)

其中${\boldsymbol S}^\Delta \left( x \right) = [ {S_1 \left( x \right)} \ \ {S_2 \left( x \right)} \ \ {S_3 \left( x \right)} \ \ {S_4 \left( x \right)} ]$为杆单元的形函数矩阵.

根据新构造的一维杆单元,重新定义单元的轴向应变为

$ \varepsilon = \dfrac{1}{2}\left( {r_x^{\Delta 2} - 1} \right) = \dfrac{1}{2}\left( {{\boldsymbol q }_\Delta ^{\rm T} {\boldsymbol S }^{\Delta {\rm T}}{\boldsymbol S }^\Delta {\boldsymbol q }_\Delta - 1} \right) $ (14)

由上式可知,重新定义的轴向应变只与梁单元长度和节点处的物质导数大小有关,与梁单元弯曲程度无关,具体可通过图 2中7组不 同弯曲程度的梁单元加以验证. 由于图中梁单元变形后长度保持不变($\varepsilon _L = 0)$,且节点处的轴向应变为0 ($\left| {{\boldsymbol r }_x (0)} \right| = \left| {{\boldsymbol r }_x (L)} \right| = 1)$,因此这7组梁对应的等效杆单元节点坐标均为${\boldsymbol q }_\Delta (t) = \left[ 0 \ \ 1 \ \ L \ \ 1 \right]^{\rm T}$,将其代入式(14)中计算可得7组轴向应变均为0,由此可知重新定义后的轴向应变独立于弯曲应变,单元应变实现解耦.

2.2 广义弹性力计算

根据式(8),并结合式(2)和式(5)可得原模型的轴向广义弹性力为

$ {\boldsymbol Q }_1 = \int_0^L EA\varepsilon \left( {\dfrac{\partial \varepsilon }{\partial {\boldsymbol q }}} \right)^{\rm T} {\rm{d}}x = \\ \qquad \dfrac{1}{2}EA\int_0^L \left( {{\boldsymbol q }^{\rm T}{\boldsymbol S }_x ^{\rm T}{\boldsymbol S }_x {\boldsymbol q } - 1} \right)\left( {{\boldsymbol S }_x ^{\rm T}{\boldsymbol S}_x {\boldsymbol q }} \right) {\rm{d}}x $ (15)

单元应变解耦后,轴向广义弹性力发生改变,需要重新进行计算以消除伪应力对模型精度的影响.需要注意的是,式中的$\left( {\dfrac{\partial \varepsilon }{\partial {\boldsymbol q }}}\right)^{\rm T}$为轴向应变对梁单元节点坐标的偏导数,而非对新构造的一维杆单元节点坐标的偏导数,由式(13)和式(14)可得解耦模型的轴向广义弹性力为

$ {\boldsymbol Q }_1^\Delta = \int_0^L EA\varepsilon \left( {\dfrac{\partial \varepsilon }{\partial {\boldsymbol q }}} \right)^{\rm T} {\rm{d}}x = \\ \qquad \dfrac{1}{2}EA\int_0^L \left( {r_x^{\Delta 2} - 1} \right)r_x^\Delta \left( {\dfrac{\partial r_x^\Delta }{\partial {\boldsymbol q }}} \right)^{\rm T} {\rm{d}}x =\\ \qquad \dfrac{1}{2}EA\int_0^L \left( {{\boldsymbol q }_\Delta^{\rm T} {\boldsymbol S }^{\Delta {\rm T}}{\boldsymbol S }^\Delta {\boldsymbol q}_\Delta - 1} \right)\left( {{\boldsymbol S }_x^\Delta {\boldsymbol q }_\Delta } \right)\left( {\dfrac{\partial {\boldsymbol q }_\Delta }{\partial {\boldsymbol Q }}} \right)^{\rm T}{\boldsymbol S }_x^{\Delta {\rm T}} {\rm{d}}x $ (16)

其中$\dfrac{\partial {\boldsymbol q }_\Delta }{\partial {\boldsymbol q }}$为杆单元节点坐标对梁单元节点坐标的偏导数,为$4\times 12$维矩阵. 根据式(3)和式(12),有

$ \dfrac{\partial {\boldsymbol q }_\Delta }{\partial {\boldsymbol q }} = \left[ \!\!\begin{array}{c} 0_{1\times 12} \\ \dfrac{{\boldsymbol r }_x (0)^{\rm T}}{\left| {{\boldsymbol r }_x (0)} \right|} \cdot {\boldsymbol S }_x (0) \\ \int_0^1 \dfrac{{\boldsymbol r }_x ^{\rm T}}{\left| {{\boldsymbol r }_x } \right|} {\boldsymbol S }_x {\rm{d}}x \\ \dfrac{{\boldsymbol r }_x (L)^{\rm T}}{\left| {{\boldsymbol r }_x (L)} \right|} \cdot {\boldsymbol S }_x (L) \end{array}\!\!\right]=\left[\!\!\begin{array}{l} {\boldsymbol A } \\ {\boldsymbol B } \\ {\boldsymbol C } \\ {\boldsymbol D } \end{array}\!\!\right] $ (17)

上式中各分量可表达如下

$ \left. \begin{array}{l} \mathit{\boldsymbol{A}} = {\boldsymbol 0_{1 \times 12}}, \mathit{\boldsymbol{B}} = [{\boldsymbol 0_{1 \times 3}}\;\;{\mathit{\boldsymbol{r}}_x}{\left( 0 \right)^{\rm{T}}}/\left| {{\mathit{\boldsymbol{r}}_x}\left( \boldsymbol 0 \right)} \right|{0_{1 \times 6}}]\\ \mathit{\boldsymbol{C}} = \int_0^1 {\frac{{\mathit{\boldsymbol{r}}_x^{\rm{T}}}}{{\left| {{\mathit{\boldsymbol{r}}_x}} \right|}}} {\mathit{\boldsymbol{S}}_x}{\rm{d}}x, \;\;\mathit{\boldsymbol{D}} = [{\boldsymbol 0_{1 \times 9}}\;\;{\mathit{\boldsymbol{r}}_x}{\left( L \right)^{\rm{T}}}/\left| {{\mathit{\boldsymbol{r}}_x}\left( L \right)} \right|] \end{array} \right\} $ (18)

新构造的一维杆单元只用于重新定义轴向应变,不直接参与动力学计算.也就是说,解耦模型仍以原ANCF模型为基础,单元广义坐标、质量矩阵、弯曲应变及其对应的广义弹性力均沿用原ANCF模型,只有轴向应变及其对应的广义弹性力发生改变,分别采用式(14)和式(16)进行计算.

3 算例分析

本节通过大变形柔性梁的静力学和动力学仿真分析,对本文提出的解耦梁模型进行验证.

3.1 静力学分析

柔性梁一端固定,另一端自由,分析其重力下的静平衡状态.具体参数如表 2所示.

表 2 柔性梁参数 Table 2 Parameters of flexible beam

分别采用耦合梁单元(原ANCF梁单元)和解耦梁单元进行建模,给出两种模型的静力平衡位形曲线,如图 7所示.为了更加清晰地区分不同的位形曲线,将末端区域进行了局部放大.随着单元数目的增加,耦合梁的弯曲程度增大,位形曲线逐渐收敛,相邻位形曲线偏差越来越小,当单元数达到18时已具有很高精度.这是由于随着单元数目的增加,单元弯曲角减小,应变耦合引起的单元刚化效应减弱,导致梁的弯曲程度增大且逐渐收敛.对比解耦梁与耦合梁的位形曲线可知,3单元解耦梁位形曲线介于6和9单元耦合梁之间,6单元解耦梁位形曲线介于15和18单元耦合梁之间,因此解耦梁模型的精度远大于同等单元数目的耦合梁模型,具有更好的收敛性.

图 7 柔性梁平衡位形 Figure 7 Equilibrium configurations of flexible beam

图 8给出了不同弹性模量下柔性梁的位形曲线,将6单元解耦梁和不同单元数的耦合梁进行了对比.在相同单元数下,解耦模型的精度明显高于耦合模型,随着单元数目增加,耦合梁的位形曲线逐渐接近解耦梁,当单元数达到12时,耦合梁与6单元解耦梁位形曲线接近,此时两者有着相近的计算精度.

图 8 不同弹性模量下柔性梁位形 Figure 8 Configurations of flexible beam with different elastic modulus

另外,对比不同弹性模量下的梁位形曲线可知,弹性模量越小,耦合梁的精度越差,相同单元数下与解耦梁末端点位置相对误差分别为10.2%($E=4$ GPa),22.4% ($E=2$ GPa),38.5%($E=1$ GPa).

梁的静平衡位形随单元数目增多逐渐收敛,其收敛性可通过静平衡势能进行分析.柔性梁的静平衡势能等于应变能和重力势能之和,其中水平位置($y =0$)为零势能位置. 图 9比较了耦合梁模型与解耦梁模型的静平衡势能收敛性曲线.通过势能的变化趋势可知,解耦梁势能的收敛性明显好于传统耦合梁. 3单元解耦梁的静平衡势能为$-11.41$ J,接近于9单元耦合梁的静平衡势能;而当解耦梁单元数目增大到6时,即已达到18单元耦合梁的静平衡势能($-11.47$ J),进一步说明了为何在相同单元数目下,解耦梁模型具有更高的位形精度(如图 7).

图 9 静平衡势能收敛性曲线 Figure 9 Convergence curves of static equilibrium energy

对比图 9图 10可知,轴向应变能占总应变能的比重与静平衡势能有着相似的变化趋势,说明解耦模型与耦合模型收敛性的差异主要由轴向应变能引起.解耦梁的轴向应变能占比接近于零,几乎不随单元数目发生变化,而耦合梁的轴向应变能占比在单元数目较少时非常大,这种伪轴向应变能使单元刚化,是导致单元精度低、收敛性差的主要因素,而解耦模型消除了这种伪应变能,因此具有更高的精度.

图 10 轴向应变能占总应变能比重 Figure 10 Proportion of axial strain energy in total strain energy

上述接近纯弯的耦合梁的轴向应变能主要是伪应变能,若柔性梁受拉力和弯矩共同作用,耦合梁的轴向应变能则包含拉力引起的应变能和应变耦合效应引起的伪应变能,解耦模型则无伪应变能.因此无论是在纯弯还是在受拉力和弯矩共同作用的情况下,解耦梁轴向应变能占总应变能比重都更小,且具有更好的收敛性.只是在纯弯的情况下该值接近于零,而在受拉力时,该值会随着拉力增大而逐渐增大.

图 11图 12对比了两种模型的轴向应变与曲率分布.由图 11可知,耦合梁的轴向应变呈现波动式分布,且波动峰值随着弯曲变形的减小而减小(如图 11).相比于6单元耦合梁模型,18单元耦合梁的轴向应变峰值更小、衰减更快,但仍大于解耦梁的轴向应变.正是由于这种轴向应变的波动式分布,使得梁的静平衡势能无法随单元数目增加而快速收敛.由图 12可知,耦合梁节点处的曲率存在明显的跳变,且跳变量随曲率的增大而增大.对比不同单元数目的耦合梁曲率可知,随着单元数目增多,曲率跳变量明显减小,曲率连续性更好.此外,相比于耦合梁,解耦梁的曲率具有更好的连续性,6单元解耦梁的曲率跳变量与18单元耦合梁的曲率跳变量相近,且不连续点更少,因此解耦模型能够提高节点处应变的连续性.

图 11 轴向应变对比 Figure 11 Comparison of axial strain
图 12 曲率对比 Figure 12 Comparison of curvature

上述算例是几乎只受纯弯作用的悬臂梁,下面考虑拉力和弯矩共同作用的情况.柔性梁一端固定,另一端施加拉力$F=300$ N和弯矩$M =30$ Nm,另外在柔性梁上施加$P =200$ N/m的分布载荷,如图 13所示.

图 13 拉力和弯矩作用下的柔性梁位形 Figure 13 Configurations of flexible beam under tension and bending moment

将6单元解耦模型与不同单元数的耦合模型进行对比可知,在相同单元数下,耦合梁弯曲程度小于解耦梁,随着单元数目增加,耦合梁位形逐渐接近解耦梁,在单元数达到12时,耦合梁与6单元解耦梁位形曲线几乎重合.因此6单元解耦梁与12单元耦合梁的计算精度相近,远大于6单元耦合梁的精度.该算例说明,除了在纯弯情况下,解耦模型在拉力和弯矩共同作用的情况下仍然适用.

3.2 动力学分析

本节进行柔性梁动力学分析,梁的具体参数如表 2所示.柔性梁一端固定,从水平位置以零初始速度自由下落,分别采用耦合梁模型和解耦梁模型进行动力学仿真. 图 14图 15给出了3单元梁末端点在水平方向上的位移和速度曲线.由图 14可知,0.9 s前两种模型末端点位移曲线接近,1 s左右耦合梁的末端点位移首先达到峰值,而解耦梁末端点继续向$x$负方向运动一段距离后达到峰值(1.2 s),最大位移相差0.7 m,峰值过后解耦梁末端点的运动明显滞后(运动周期更长).由图 15可知,相比于解耦梁,耦合梁的末端速度存在明显的高频振动,这是由应变耦合引起的单元刚化效应造成的;而解耦模型消除了应变耦合效应,因此速度曲线中不存在高频振动.

图 14 末端点位移(3单元) Figure 14 Displacement of endpoint (3 elements)
图 15 末端点速度(3单元) Figure 15 Velocity of endpoint (3 elements)

图 16图 17为6和30单元模型仿真结果.相比于3单元模型,耦合梁和解耦梁末端点位移和速度曲线更为接近,其中30单元耦合梁和6单元解耦梁的位移和速度曲线几乎完全重合.耦合梁末端点速度不再有高频振动(如图 17),可知随着单元数量增加,应变耦合引起的刚化效应减弱.相比于6单元耦合梁,30单元耦合梁和6单元解耦梁一样,其末端点速度曲线在1.5 s后出现了一定波动,这是由于单元数增多后,应变耦合效应减弱,模型刚度减小,使得柔性梁摆动过程中的波动效应更加明显.

图 16 末端点位移(6和30单元) Figure 16 Displacement of endpoint (6 and 30 elements)
图 17 末端点速度(6和30单元) Figure 17 Velocity of endpoint (6 and 30 elements)

综上,单元应变解耦后,梁的运动范围增大,运动周期延长,消除了应变耦合效应引起的高频速度振动,在相同单元数目下具有更高的精度.

4 结论

本文对ANCF索梁单元的应变耦合问题及模型解耦方法进行研究,得到如下结论:

(1)梁单元轴向应变与弯曲应变相互耦合,单元的弯曲变形会引起伪轴向应变,产生的伪应变能导致单元刚度大于理论值,且单元弯曲变形越大,模型失真越严重.

(2)构造了等效一维杆单元,基于该单元重新定义梁的轴向应变,实现了轴向应变与弯曲应变解耦,在此基础上重新推导广义弹性力,得到了梁的应变解耦模型.

(3)通过静力学和动力学仿真分析,对本文提出的解耦梁模型进行了验证.仿真结果表明:解耦梁模型消除了伪轴向应变,相比于原模型具有更好的收敛性和曲率连续性,此外,解耦梁模型消除了原模型中的高频速度振动,在相同单元数目下具有更高的精度.


参考文献
1 Shabana AA. Flexible multibody dynamics:review of past and recent developments. Multibody System Dynamics, 1997, 1 (2) : 189-222. DOI: 10.1023/A:1009773505418.
2 田强.基于绝对节点坐标方法的柔性多体系统动力学研究与应用.[博士论文].武汉:华中科技大学, 2009 ( Tian Qiang. Flexible multibody dynamics research and application based on the absolute nodal coordinate method.[PhD Thesis]. Wuhan:Huazhong University of Science and Technology, 2009(in Chinese) ) http://www.oalib.com/references/16762580
3 Shabana AA. An absolute nodal coordinate formulation for the large rotation and large deformation analysis of flexible bodies. Technical Report No. MBS96-1-UIC, Department of Mechanical Engineering, University of Illinois at Chicago, 1996
4 Shabana AA. Definition of the slopes and the finite element absolute nodal coordinate formulation. Multibody System Dynamics, 1997, 1 (3) : 339-348. DOI: 10.1023/A:1009740800463.
5 Shabana AA. Definition of ANCF finite elements. Journal of Computational and Nonlinear Dynamics, 2015, 10 (5) : 054506. DOI: 10.1115/1.4030369.
6 陈思佳, 章定国, 洪嘉振. 大变形旋转柔性梁的一种高次刚柔耦合动力学模型. 力学学报, 2013, 45 (2) : 251-256. ( Chen Sijia, Zhang Dingguo, Hong Jiazhen. A high-order rigid-flexible coupling model of a rotating flexible beam under large deformation. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45 (2) : 251-256. (in Chinese) )
7 章孝顺, 章定国, 洪嘉振. 考虑曲率纵向变形效应的大变形柔性梁刚-柔耦合动力学建模与仿真. 力学学报, 2016, 48 (3) : 1-11. ( Zhang Xiaoshun, Zhang Dingguo, Hong Jiazhen. Rigid-flexible coupling dynamic modeling and simulation with the longitudinal deformation induced curvature effect for a rotating flexible beam under large deformation. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48 (3) : 1-11. (in Chinese) )
8 Shabana AA, Yakoub RY. Three dimensional absolute nodal coordinate formulation for beam elements:theory. Journal of Mechanical Design, 2001, 123 (4) : 606-613. DOI: 10.1115/1.1410100.
9 Yakoub RY, Shabana AA. Three dimensional absolute nodal coordinate formulation for beam elements:implementation and applications. Journal of Mechanical Design, 2000, 123 (4) : 614-621.
10 Gerstmayr J, Shabana AA. Analysis of thin beams and cables using the absolute nodal coordinate formulation. Nonlinear Dynamics, 2006, 45 (1-2) : 109-130. DOI: 10.1007/s11071-006-1856-1.
11 朱大鹏, 路英杰, 汤家力等.大变形索梁单元在多体动力学框架下的应用//2007全国结构动力学学术研讨会论文集, 南昌, 2007年11月 ( Zhu Dapeng, Lu Yingjie, Tang Jiali, et al. Application of large deformation cable beam element in multibody dynamics//Proceedings of the 2007 National Symposium on Structural Dynamics, Nanchang, 2007-11(in Chinese). )
12 Dombrowski SV. Analysis of large flexible body deformation in multibody systems using absolute coordinates. Multibody System Dynamics, 2002, 8 (4) : 409-432. DOI: 10.1023/A:1021158911536.
13 Liu C, Tian Q, Hu HY. New spatial curved beam and cylindrical shell elements of gradient-deficient absolute nodal coordinate formulation. Nonlinear Dynamics, 2012, 70 (3) : 1903-1918. DOI: 10.1007/s11071-012-0582-0.
14 Dmitrochenko O, Mikkola A. Digital nomenclature code for topology and kinematics of finite elements based on the absolute nodal coordinate formulation. Journal of Multi-body Dynamics, 2011, 225 (1) : 34-51.
15 Kerkkänen KS, García-Vallejo D, Mikkola AM. Modeling of beltdrives using a large deformation finite element formulation. Nonlinear Dynamics, 2006, 43 (3) : 239-256. DOI: 10.1007/s11071-006-7749-5.
16 ČeponG, BoltežarM. Dynamics of a belt-drive system using a linear complementarity problem for the belt-pulley contact description. Journal of Sound and Vibration, 2009, 319 (3) : 1019-1035.
17 ČeponG, ManinL, BoltežarM. Introduction of damping into the flexible multibody belt-drive model:a numerical and experimental investigation. Journal of Sound and Vibration, 2009, 324 (1) : 283-296.
18 Lugris U, Escalona J, Dopico D, et al. Efficient and accurate simulation of the cable-pulley interaction in weight-lifting machines//1st Joint International Conference on Multibody System Dynamics, Lappeenranta, Finland, 2010
19 Kim KW, Lee JW, Yoo WS. The motion and deformation rate of a flexible hose connected to a mother ship. Journal of Mechanical Science and Technology, 2012, 26 (3) : 703-710. DOI: 10.1007/s12206-011-1202-5.
20 Tur M, García E, Baeza L, et al. A 3D absolute nodal coordinate finite element model to compute the initial configuration of a railway catenary. Engineering Structures, 2014, 71 : 234-243. DOI: 10.1016/j.engstruct.2014.04.015.
21 Vohar B, Kegl M, Ren Z. Implementation of an ANCF beam finite element for dynamic response optimization of elastic manipulators. Engineering Optimization, 2008, 40 (12) : 1137-1150. DOI: 10.1080/03052150802317457.
22 Sanborn GG, Choi J, Choi JH. Curve-induced distortion of polynomial space curves, flat-mapped extension modeling, and their impact on ANCF thin-plate finite elements. Multibody System Dynamics, 2011, 26 (2) : 191-211. DOI: 10.1007/s11044-011-9248-9.
23 Sopanen JT, Mikkola AM. Description of elastic forces in absolute nodal coordinate formulation. Nonlinear Dynamics, 2003, 34 (1-2) : 53-74.
24 Schwab AL, Meijaard JP. Comparison of three-dimensional flexible beam elements for dynamic analysis:finite element method and absolute nodal coordinate formulation//ASME 2005 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. Long Beach, California, USA:ASME, 2005:1341-1349
25 García-Vallejo D, Mikkola AM, Escalona JL. A new locking-free shear deformable finite element based on absolute nodal coordinates. Nonlinear Dynamics, 2007, 50 (1-2) : 249-264. DOI: 10.1007/s11071-006-9155-4.
26 Gerstmayr J, Matikainen MK, Mikkola AM. A geometrically exact beam element based on the absolute nodal coordinate formulation. Multibody System Dynamics, 2008, 20 (4) : 359-384. DOI: 10.1007/s11044-008-9125-3.
27 Nachbagauer K, Gruber P, Gerstmayr J. Structural and continuum mechanics approaches for a 3D shear deformable ANCF beam finite element:application to static and linearized dynamic examples. Journal of Computational and Nonlinear Dynamics, 2013, 8 (2) : 021004.
28 于洋, 宝音贺西, 李俊峰. 空间飞网抛射展开动力学建模与仿真. 宇航学报, 2010, 31 (5) : 1289-1296. ( Yu Yang, Baoyin Hexi, Li Junfeng. Modeling and simulation of projecting deployment dynamics of space webs. Journal of Astronautics, 2010, 31 (5) : 1289-1296. (in Chinese) )
29 张江.空间绳网捕获过程碰撞动力学研究.[硕士论文].哈尔滨:哈尔滨工业大学, 2015 ( Zhang Jiang. Contact dynamics of space net on capturing target.[Master Thesis]. Harbin:Harbin Institute of Technology, 2015(in Chinese). )
30 Zhang Y, Wei C, Pan D, et al. A dynamical approach to space capturing procedure using flexible cables. Aircraft Engineering and Aerospace Technology, 2016, 88 (1) : 53-65. DOI: 10.1108/AEAT-07-2014-0107.
31 Tang JL, Ren GX, Zhu WD, et al. Dynamics of variable-length tethers with application to tethered satellite deployment. Communications in Nonlinear Science and Numerical Simulation, 2011, 16 (8) : 3411-3424. DOI: 10.1016/j.cnsns.2010.11.026.
THE STRAIN COUPLING PROBLEM AND MODEL DECOUPLING OF ANCF CABLE/BEAM ELEMENT1)
Zhang Yue*, Zhao Yang*2), Tan Chunlin, Liu Yongjian     
*. Harbin Institute of Technology, Harbin 150001, China;
†. China Academy of Space Technology, Beijing 100086, China
Abstract: The cable/beam structure has already been widely applied in civil engineering, aerospace engineering, etc. Among various dynamic modeling methods of cable/beam, the Absolute Nodal Coordinate Formulation (ANCF) is very suitable for the modeling of large deformed cable/beam to describe the large deformation and rotation of flexible bodies. According to the strain analysis of ANCF cable/beam element, the bending deformation will cause uneven distribution of axial strain within the element, which means axial and bending strains are coupled with each other. The strain coupling effect brings unrealistic strain energy to the element, resulting in an increased stiffness and element distortion. It is known from the analysis of strain and strain energy at various bending angles that the larger the bending deformation is, the more serious the distortion is. The axial and bending strains are decoupled by redescribing axial strain based on a new constructed equivalent 1D rod element. Then the generalized elastic force is deduced and the strain-decoupled model of ANCF cable/beam element is obtained. The statics and dynamics simulations of two beam models are conducted and the results indicate that the unrealistic strain is eliminated by the decoupled model, and compared with the original model, the decoupled model shows better convergence and curvature continuity and more accurate under same number of elements. Meanwhile, as the stiffness is reduced by the decoupled model, there is no high frequency vibration in the velocity curve any more by compared with the original model.
Key words: absolute nodal coordinate    cable    beam    strain coupling    curvature continuity