引言

索梁结构在土木工程、航空航天等领域有着广泛的应用，其结构参数设计、系统任务分析等均依赖于精确的索梁动力学模型.对于大变形索梁结构，基于小变形、小转动的传统建模方法已不再适用^[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)，该单元适用于细长索和梁的动力学建模.为方便叙述，后文将该索梁单元简称为梁单元.

梁中心轴线上任意点的位置坐标可表示为关于物质坐标$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所示.

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

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

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

为了更清晰地体现轴向应变能和弯曲应变能的变化，图中对弯曲角小于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}$，从而引起单元刚化. 单元弯曲变形越大，刚化效应越显著，单元失真越严重.

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

2 梁单元的应变解耦模型

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

2.1 应变解耦

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

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

$ {\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所示.

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

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

另外，对比不同弹性模量下的梁位形曲线可知，弹性模量越小，耦合梁的精度越差，相同单元数下与解耦梁末端点位置相对误差分别为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和图 10可知，轴向应变能占总应变能的比重与静平衡势能有着相似的变化趋势，说明解耦模型与耦合模型收敛性的差异主要由轴向应变能引起.解耦梁的轴向应变能占比接近于零，几乎不随单元数目发生变化，而耦合梁的轴向应变能占比在单元数目较少时非常大，这种伪轴向应变能使单元刚化，是导致单元精度低、收敛性差的主要因素，而解耦模型消除了这种伪应变能，因此具有更高的精度.

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

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

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

将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可知，相比于解耦梁，耦合梁的末端速度存在明显的高频振动，这是由应变耦合引起的单元刚化效应造成的；而解耦模型消除了应变耦合效应，因此速度曲线中不存在高频振动.

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

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

4 结论

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

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

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

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