基于旋转场曲率的二维剪切梁单元建模∗

张大羽罗建军†郑银环袁建平

1)(西北工业大学航天学院,航天飞行动力学技术重点实验室,西安 710072)

2)(武汉理工大学机电工程学院,武汉430070)

基于旋转场曲率的二维剪切梁单元建模∗

张大羽1)罗建军1)†郑银环2)袁建平1)

1)(西北工业大学航天学院,航天飞行动力学技术重点实验室,西安 710072)

2)(武汉理工大学机电工程学院,武汉430070)

(2016年12月31日收到;2017年4月6日收到修改稿)

对二维剪切梁单元进行研究,利用平面旋转场理论推导了精确曲率模型.采用几何精确梁理论构建了剪切梁单元弹性力矩阵.通过绝对节点坐标方法建立了系统的非线性动力学方程,提出基于旋转场曲率的二维剪切梁单元,并分别引入经典二维剪切梁单元和基于位移场曲率的二维剪切梁单元进行比较研究.首先,静力学分析证明了所提模型的正确性;其次,特征频率分析验证了模型可与理论解符合,收敛精度高,并且能准确地预测单元固有频率对应的振型;最后,在非线性动力学问题上,通过与ANSYS结果对比分析,证明了该模型可有效处理柔性大变形问题,并且与经典二维剪切梁单元相比具有缓解剪切闭锁的优势.因此,本文提出的基于旋转场曲率的二维剪切梁单元在处理几何非线性问题中具有较大的应用潜力.

绝对节点坐标法,旋转场曲率,几何精确梁理论,剪切闭锁

1 引言

近年来,空间碎片清理技术在航天界受到广泛的关注[1].世界主要航天大国针对主动俘获碎片技术均投入了大量的人力财力进行研究,其中飞网、绳系机械爪等柔性俘获装置越来越多地出现在空间飞行器上.由于该类俘获装置具有超轻、易变形及强几何非线性等特点,在执行抓捕碎片过程中自身会产生柔性大变形,使飞行器的姿态稳定及精确控制面临极大的挑战[2].梁单元是这类俘获装置结构中的基本单元.针对其大变形问题,采用传统基于小变形假设的梁变形理论不能真实地反映大变形、大位移等特点,如运动弹性动力学分析、一次近似耦合模型[3]、一次精确模型[4,5]、高次耦合模型[6].而Shabana[7]提出了绝对节点坐标法(absolute nodal coordinate formulation,ANCF),该方法采用绝对位置坐标和变形梯度描述物体的刚体运动和柔性变形,适用于处理大变形、大转动多体系统的建模问题.此外,该方法不采用转动角坐标,并且所有节点坐标均定义在惯性系中,因此构建出的质量矩阵为常数阵,不存在柯氏加速度项和离心加速度项,惯性力中也不存在变形和转动的附加耦合项.

进入21世纪以来,已有越来越多的学者关注ANCF方法的研究.田强等[8]对近年来ANCF方法的研究做了系统的评述.在ANCF剪切梁单元研究中,Omar和Shabana[9]基于连续介质力学理论最先提出了经典的二维ANCF全参数剪切梁有限元单元,该单元采用x向三次插值,y向线性插值,这会导致严重的剪切闭锁问题,给数值计算带来困难.随后,很多学者针对该闭锁问题对此单元进行改进,更好地提升了单元的性能.Hussein等[10]研究了二维ANCF全参数剪切梁单元中变形模态耦合情况,并采用Hellinger-Reissner变分原理解决剪切闭锁问题.随后,他们还将该研究思路应用到三维ANCF剪切梁单元和ANCF板单元[11].García-Vallejo等[12]提出了二维三节点ANCF剪切梁单元,不同于传统二维ANCF剪切梁单元的完备位移场,该单元的位移场是不完备的三次插值多项式,在此处理下剪切力呈二次分布并且梁截面在弯曲作用时不会发生收缩,解决了闭锁问题,动力学仿真验证了二维三节点ANCF剪切梁的有效性.Tian等[13]采用García-Vallejo开发的locking-free ANCF剪切梁单元进一步研究了考虑关节间隙的二维剪切梁动力学模型,丰富了绝对节点坐标法的单元类型.Gerstmayr等[14]基于绝对节点坐标法,采用Reissner非线性梁理论重新构造了二维剪切梁单元的弹性力矩阵,公式中的轴向力、剪切力和弯曲力表述成显示解耦形式.此外,他定义了厚度方向应变能用以表征横截面的变形,并准确地描述了剪切梁横截面固有频率对应的振型.Gerstmayr定义这种弹性力解耦形式的梁单元为几何精确梁单元(geometrically exact beam element),该方法对各类型应变能的积分只沿梁的中心线,各个变形模态之间互不影响,该处理可以有效缓解剪切闭锁.Nachbagauer等[15]只采用梁轴线上的梯度值rx|η=0和ry|η=0作为变形自由度,分别开发出线性和二次ANCF非全参数二维剪切梁单元.由于新单元的梯度ry|η=0不能表征梁横截面的变形,因此该处理成为解决剪切闭锁问题的一种新思路.随后他们还开发出了三维ANCF剪切梁单元,解决了三维ANCF剪切梁单元存在的闭锁问题并显示出高阶收敛性[16].Gerstmay和Shabana[17]针对ANCF剪切梁的闭锁现象又提出了两种解决方法:一种为在ANCF剪切梁单元基础上对其形函数中分别引入x2y和x2z的插值项,开发了高阶三维ANCF梁单元,新增的插值项是截面关于轴向x的二次函数,这会使弯曲应变变成x的线性函数,并与截面对x的二阶导同阶,消除了剪切闭锁现象,更好地描述了截面的几何变形;另一种方法是将原有三维ANCF单元降阶并退化成三维低阶ANCF索单元,该单元中没有y方向梯度,相比原三维ANCF梁单元的自由度数减少了一半,计算效率大幅提升,同时系统中也没有与剪切变形相关的刚度项,故剪切闭锁问题不存在.Dufva等[18]通过截面转角描述二维ANCF剪切梁单元的精确位移场,梁轴向应变和弯曲应变由截面转角表征,而剪切变形采用单独的线性插值,该单元采用应变场混合插值技术解决了剪切闭锁问题.Mikkola等[19]在二维ANCF非全参数梁的基础上通过定义截面方向上的力平衡条件引入剪切变形,开发出一种新的二维ANCF剪切梁单元.新单元中因不含截面梯度自由度,计算效率远高于传统二维ANCF全参数剪切梁单元,并且也不存在数值计算时由截面梯度引起的高频问题.他们提出的新单元与传统Timoshenko梁单元弯曲模态和轴向模态符合度良好,证明了新单元不存在曲率闭锁问题.Vesa-Ville等[20]在二维ANCF非全参数梁单元的基础上采用对梁y向变形场独立插值技术,提出了一种基于混合插值的二维ANCF剪切梁单元,节点坐标由位移场、截面方向和截面变形场组成,可捕捉到剪切变形.动力学仿真结果表明该单元具有更高的精度和收敛性.章孝顺等[21]对二维ANCF非全参数梁单元的曲率模型进行讨论,对比了精确曲率模型和两种简化模型,同时还研究了不同简化模型的适用范围.以上研究较全面地论述了二维ANCF剪切梁的发展历程,但针对二维剪切梁单元的曲率几何非线性进行的研究较少.已有的曲率模型大都直接采用二阶曲率公式[22],定义在Frenet坐标系下,具有极值性;且鲜有从旋转场角度研究曲率定义,并对基于旋转场曲率的二维ANCF剪切梁单元进行几何非线性分析、讨论模型的正确性和适用性.

本文进一步针对二维ANCF剪切梁的曲率定义进行研究.基于平面旋转场理论推导得到精确曲率公式,表征为梁中心线切线方向转角的函数;采用几何精确梁理论构建系统弹性力矩阵,提出基于旋转场曲率的二维ANCF剪切梁单元模型;同时引入基于连续介质力学理论的经典ANCF剪切梁单元和基于微分几何理论的位移场曲率ANCF剪切梁单元作为对比研究.分别通过静力学分析、特征频率分析以及动力学分析验证了本文所提的基于旋转场曲率的二维ANCF剪切梁单元模型在计算柔性机构大变形、大转动问题中的准确性和适用性.

文中第二部分分别建立基于连续介质力学理论的经典剪切梁单元(模型一)、基于位移场曲率的几何精确剪切单元(模型二)和基于旋转场曲率的几何精确剪切梁单元(模型三);第三部分引入小变形和非小变形静力学分析、特征频率分析以及非线性动力学分析验证所提单元的有效性;第四部分为结论.

2 剪切梁模型

2.1 模型一:基于连续介质力学理论的二维ANCF剪切梁单元

二维含剪切变形的ANCF梁单元最早由Omar和Shabana提出[9],如图1所示,其中O-XY为惯性坐标系,l为梁的长度.

对于二节点的ANCF剪切梁单元i,其梁上任意一点在惯性系中的绝对位移矢量为

图1 剪切梁截面变形Fig.1.Cross section deformation of shear deformable beam.

其中x,y为梁初始构型在单元坐标系内的坐标.该单元有12个绝对节点坐标,每个节点有6个自由度,包括2个位移自由度和4个梯度自由度:

梁单元i的型函数Si(x,y)为:

其中I是2×2的单位矩阵.

利用连续介质力学中的变形理论推导得到ANCF剪切梁单元i的变形梯度张量Ji:

对于梁单元的初始状态是水平摆放情况,J0为单位阵.Green-Lagrangian应变张量εi的表达式为

依据Voigt标记规则,可将(5)式改写为

依据广义胡克定律,ANCF剪切梁单元i的应力张量σi为

E是材料的弹性系数矩阵.在平面应变问题中(plane-strain problem),对于各向同性的匀质材料(isotropic homogenous material),其弹性系数矩阵表达式为

λ,µ是拉梅常数(Lame’s constants),

E,ν分别是材料的杨氏弹性模量和泊松比.利用(6)式与(7)式可得到梁单元i的应变能为

则其弹性力矢量为

其中,

经典剪切梁单元模型中轴向伸长变形、剪切变形和弯曲变形相互耦合在一起,因此曲率公式不是显式形式.

2.2 模型二:基于位移场曲率的二维几何精确剪切梁单元

根据微分几何理论[22],曲率是单位切向量对弧长的一阶导矢的模.本例中切向量是梁中心线上的位移梯度,通过位移场推导得到的曲率模型,本文定义为位移场曲率.图1中梁变形后其截面不再与梁的中心轴线s垂直.根据梁位移场,梁单元i上任意一点变形后的位移为:ri=ri(x,y).其对中心轴线s的一阶导为

对中心线s的二阶导为

E,G,A分别是梁单元i的弹性模量、剪切模量和横截面积;剪切模量G=E/2(1+ν);As=ks·A,ks是剪切修正因子,对于横截面是矩形的情况,ks=10(1+ν)/12+11ν;应变Γi沿梁中心线的分量,可由拉格朗日应变张量确定[10]:

梁单元i的弹性力公式可由应变势能Ui对节点坐标ei求偏导得到:

2.3 模型三:基于旋转场曲率的二维几何精确剪切梁单元

曲率也可定义为曲线上某点的切线方向角对弧长的转动率.通过角度场推导得到的曲率模型,本文定义为旋转场曲率.由图1所示,根据平面旋转场理论,由O-XY逆时针转动到移动矢量的转动矩阵为,其中α是与OX轴的夹角,即梁中线的切线方向角.沿着梁中心线切线方向,可体现梁中线弯曲效应.则单位切向量可写成

由(26)式可知:cosα=ˆx1,由该式可得到旋转场曲率公式:

该模型中的轴向势能与剪切势能与模型二中的表达式相同,弯曲势能可利用旋转场曲率得到:.弯曲势能对节点坐标ei的导数为

(30)式中的偏导项分别为:

同理,通过(26)式的第二式ˆx2=sinα也可推导出旋转场曲率的表达式.

根据绝对节点坐标方法,由虚功原理可推导梁单元i的动力学方程[9]:分别采用模型一至三中的弹性力表达式,¨ei是节点加速度坐标.对于ANCF剪切梁单元i,其单元质量阵公式为

其中,ρi,Vi是梁单元i的密度和体积.由(35)式可知,通过绝对节点坐标方法建立的单元质量阵为常数阵,不包含科氏力与离心力,不存在变形和转动耦合附加的惯性力.通过虚功原理可推导单元的外力阵公式:

Fi为外力矢量,外力阵同样为常数矩阵.

3 算例分析

3.1 静力学分析

本节引入典型的悬臂梁受集中载荷问题,考察所提出基于旋转场曲率的ANCF剪切梁模型在静力学分析中的准确性.图1所示为悬臂梁受集中载荷模型,其物理参数与文献[14]相同.长为l=2 m,宽为w=0.1 m,高为h=0.5 m,杨氏模量为E=2.07×105GPa,剪切模量为G=7.96×104GPa,泊松比为ν=0.3,以及剪切修正因子为

图2 悬臂梁模型Fig.2.Cantilever beam.

考察小变形算例,设置末端集中力为F=−5×105×h3N.在小变形假设下,梁中性层在变形后长度不变,端点沿X方向的变形可以忽略不计.Y向挠度的理论解[23]为,其中Iz为梁截面惯性矩.通过表1的数据对比发现,由于模型一存在严重的剪切闭锁问题,其Y向位移与理论解误差达到22%.而模型二和模型三采用几何精确梁理论构建系统弹性力,其中轴向力、剪切力和弯曲力相互解耦,可缓解剪切闭锁,因此由它们计算的Y向位移均可收敛于理论解,并且收敛速度快.减小梁的高度至h=0.1 m,末端集中力为F=−1×106×h3N,其他参数不变.表2结果表明模型一的误差为25%,而模型二和三可收敛到理论解.与表1相比,表2中变形量增大,模型二和三需要收敛到理论解的个数增多.

表2 小变形情况下三种弹性力模型计算的端部位移对比2Table 2.2ndcomparison of tip point displacement using three elastic force models in case of small deformation.

表3 非小变形情况下三种弹性力模型计算的端部位移对比Table 3.Comparison of tip point displacement using three elastic force models in case of non-small deformation.

考察非小变形情况,设置末端集中力为F=−5×108×h3N,其他参数与小变形算例1相同.采用该集中力可导致相对大的变形量.采用文献[14]表2中用商业软件计算的解作为非线性大变形下悬臂梁X,Y向挠度的参考解.在大变形情况下,剪切效应加剧,另外由于ANCF单元在Y向上采用线性插值会导致严重的剪切闭锁问题,所以在大变形算例中三种模型的收敛精度均不如小变形算例,见表3.数据对比显示模型三的收敛速度和精度最高.以上算例证明了所提出的模型在静力学分析中的正确性.

3.2 特征频率分析

本节对所提出的ANCF剪切梁单元的特征频率进行分析,选择具有理论解的简支梁作为研究对象.简支梁物理参数与文献[14]中相同.其中长为l=2 m,宽为w=0.4 m,高为h=0.4 m,密度为ρ=7850 kg/m3,杨氏模量为E=1000 GPa,剪切模量为G=384.62 GPa,泊松比为ν=0.3,以及剪切修正因子为本算例中的边界条件如图3所示,左端的X,Y方向位移固定,右端只固定Y方向位移.

图3 简支梁模型Fig.3.Simply supported beam model.

图4 基于旋转场曲率的简支梁模型的振型和固有频率Fig.4.Mode shapes and eigenfrequencies(rad/s)of simply supported beam based on rotation-based curvature.

表43 种弹性力模型下简支梁的弯曲和轴向固有频率(rad/s)对比Table 4.Bending and longitudinal eigenfrequencies(rad/s)of simply supported beam element for di ff erent elastic formulation.

图4为本文模型三在简支边界条件下的固有频率和振型.因为通过边界条件消除了3个自由度,所以共计算出有9种振型和对应的固有频率.ANCF单元的轴向采用三次插值技术产生三种轴向振型(mode 2,mode 4和mode 5)和两种弯曲振型(mode 1和mode 3).另外,还存在两种剪切振型(mode 6和mode 7)以及ANCF单元特有的厚度振型(mode 8和mode 9).

表53 种弹性力模型下简支梁的剪切和厚度固有频率(rad/s)对比Table 5.Shearing and thickness eigenfrequencies(rad/s)of simply supported beam element for di ff erent elastic formulation.

表4为简支梁单元分别采用3种弹性力模型所计算出的弯曲和轴向固有频率的对比,表5为其剪切和厚度固有频率对比.由结果可知,与其他两个模型相比,本文提出的模型三收敛性更好,收敛精度可至6位有效数字.引起二阶弯曲振动频率误差较高是因为本文单元在Y向上采用线性插值,会导致单元的剪切闭锁问题,使得变形过程中剪切应变能占据主导,过多的剪切力导致单元不能发生正常弯曲变形,造成弯曲刚度增大,单元变得更“刚”,从而引起二阶弯曲振动频率增大.

3.3 动力学分析

图5所示为柔性单摆,O点为旋转铰.该算例是绝对节点坐标方法的一个经典算例.单摆的物理参数如下:长为l=1.2 m,横截面积为A=0.0016 m2,截面惯性矩为Iz=8.533×10−5m4,密度为ρ=5540 kg/m3,杨氏模量为E=7 GPa,剪切模量为G=2.69 GPa,泊松比以及剪切修正因子与前面算例相同.如图5所示,仅固定左端的X,Y方向位移.考察在重力作用下单摆的动态响应.仿真中单摆被划分为12个单元,仿真时间为1.1 s.

本算例中柔性单摆在其重力作用下进行大范围转动,运动过程中发生大变形现象.为对比文中所述的三种剪切梁单元模型处理大变形问题的正确性,采用ANSYS中两节点平面梁BEAM3单元作为对比.由图6—图8可见,在弹性模量为E=7 GPa时,三种模型计算的末端横向位移,中点横向位移和末端端点轨迹均与ANSYS结果符合,证明它们均能处理大变形问题,并验证了本文所提模型三的正确性.

图5 柔性单摆模型Fig.5.Flexible pendulum.

为了更好地说明模型三处理柔性大变形的能力,将弹性模量降低至E=2 GPa.在此弹性模量下,梁将变得更柔,大变形现象更明显.图9—图11显示在E=2 GPa时,三种模型计算的结果均与ANSYS结果符合.继续将弹性模量降至E=0.7 GPa,由图12可见,在此低模量下三种模型计算的结果也可与ANSYS结果相近.以上算例均证明了模型三可用于处理柔性大变形问题.

图6 柔性单摆末端横向位移(E=7 GPa)Fig.6.The tip vertical displacement of fl exible pendulum(E=7 GPa).

图7 柔性单摆中点横向位移(E=7 GPa)Fig.7.The middle vertical displacement of fl exible pendulum(E=7 GPa).

图8 柔性单摆末端位移(E=7 GPa)Fig.8.The tip-point trajectory of fl exible pendulum(E=7 GPa).

图9 柔性单摆末端横向位移(E=2 GPa)Fig.9.The tip vertical displacement of fl exible pendulum(E=2 GPa).

图10 柔性单摆中点横向位移(E=2 GPa)Fig.10.The middle vertical displacement of fl exible pendulum(E=2 GPa).

图11 柔性单摆末端位移(E=2 GPa)Fig 11.The tip-point trajectory of fl exible pendulum(E=2 GPa).

图12柔性单摆末端位移(E=0.7 GPa)Fig.12.The tip-point trajectory of fl exible pendulum(E=0.7 GPa).

图13 —图15是模型三在不同弹性模量下的能量变化曲线.由图可知,该模型总能量为零,符合能量守恒定律.从能量的角度验证了模型三的正确性.

图13 系统能量分布(旋转场曲率模型,E=7 GPa)Fig.13.The energy distribution of fl exible pendulum(RB curvature beam model,E=7 GPa).

图14 系统能量分布(旋转场曲率模型,E=2 GPa)Fig.14.The energy distribution of fl exible pendulum(RB curvature beam model,E=2 GPa).

图15系统能量分布(旋转场曲率模型,E=0.7 GPa)Fig.15.The energy distribution of fl exible pendulum(RB curvature beam model,E=0.7 GPa).

图16 —图18显示了模型一至模型三在不同时刻的构型图,分别对三种模型取间隔0.1 s的空间构型进行对比.图中绿色实线代表模型一,蓝色点画线代表模型二,红色虚线代表模型三.由图16可见,在E=7 GPa时三种模型的运动轨迹一致,说明在此弹性模量下由几何精确梁理论推导的模型(模型二和模型三)和由连续介质力学理论推导的模型(模型一)在模拟单摆大变形时的效果基本相同.随着弹性模量的增大,模型二与模型三在运动过程中相较模型一显得“更柔”,主要因为模型一是耦合模型,而模型二与模型三是解耦模型,各个变形模态之间互不影响.所以图17和图18中模型二和模型三的运动轨迹和空间构型更为接近,验证了两种不同物理场得到的曲率模型的有效性.

图16 (网刊彩色)柔性单摆各时刻构型图(E=7 GPa)(绿色实线,模型一;蓝色点画线,模型二;红色虚线,模型三)Fig.16.(color online)The con fi guration of fl exible pendulum at each given time(E=7 GPa)(green solid line,1stmodel;blue dash dotted line,2ndmodel;red dash line:3rdmodel).

图17 (网刊彩色)柔性单摆各时刻构型图(E=2 GPa)(网刊彩色)(绿色实线:模型一;蓝色点画线:模型二;红色虚线:模型三)Fig.17.(color online)The con fi guration of fl exible pendulum at each given time(E=2 GPa)(green solid line,1stmodel;blue dash dotted line,2ndmodel;red das hline,3rdmodel).

图18 (网刊彩色)柔性单摆各时刻构型图(E=0.7 GPa)(绿色实线,模型一;蓝色点画线,模型二;红色虚线,模型三)Fig.18.(color online)The con fi guration of fl exible pendulum at each given time(E=0.7 GPa)(green solid line,1stmodel;blue dash dotted line,2ndmodel;red dash line:3rdmodel).

柔性梁在大转动过程中剪切效应逐渐占据主导,为了进一步研究三种不同的弹性力计算的剪切应变在大转动大变形中的变化,弹性模量选取E=0.7 GPa,其余梁的物理参数不变,柔性单摆被划分为40个单元.图19是单摆最后一个单元内中点的剪切应变变化曲线.从图中可知,模型二与模型三得到的剪切应变相近,而它们均比模型一的剪切应变大,其最大差异可达到0.01.其差异可由图20和图21解释,分别是单摆最后一个单元中点的rx,ry方向与大小变化曲线.由图可知,模型二与模型三在rx,ry方向变化上与模型一基本相似,而rx,ry的幅值变化与模型一有较大差异,这也是导致图19剪切应变差异的主要原因.绝对节点坐标法中的应变是位移梯度的函数,例如正应力表达式

它可限制ry剧烈变化,确保单元的不可压缩性.模型一由于是耦合模型,rx与ry变化相互影响,所以其ry幅值的变化范围比模型二和模型三要更明显.

图19 单摆最后一个单元内中点的剪切应变变化Fig.19.Change of shear strain at middle point of last element.

由于ANCF剪切梁单元在x,y方向上插值阶次的不同将会导致严重的剪切闭锁现象,采用本文模型一至模型三分别研究各模型缓解剪切闭锁的能力.表6为分别采用三种模型时的单摆中点剪切应变收敛情况,可以看出采用解耦模型(模型二和模型三)得到的剪应变收敛个数比耦合模型(模型一)要少,模型三收敛最快,证明模型三在缓解剪切闭锁问题上具有较强的优势.模型一由于各个变形模态相互耦合,剪切应变收敛最慢.

图20 单摆最后一个单元内中点的rx和ry方向变化(a)rx;(b)ryFig.20.Orientation of the rxand ryvector at middle point of last element:(a)rx;(b)ry.

图21 单摆最后一个单元内中点的rx和ry幅值变化(a)rx;(b)ryFig.21.Magnitude of the rxand ryvector at middle point of last element:(a)rx;(b)ry.

表6 三种弹性力模型的单摆中点剪切应变收敛对比Table 6.Comparison of shear strain at middle point using three di ff erent elastic force models

4 结论

本文采用梁中心线切线方向角推导了剪切梁的精确曲率公式,提出了基于旋转场曲率的二维ANCF剪切梁单元模型.通过与传统二维ANCF剪切梁单元模型、基于位移场曲率的二维ANCF单元模型对比分析,可得到以下结论:

1)在静力学分析以及特征频率分析中,本文提出的模型三与理论解收敛度最好,效率最高,且收敛精度和速度优于模型一;

2)在动力学分析中,模型三显示了良好的计算精度,可有效预测柔性构件运动过程中的大变形现象,此外,对比验证了两种不同物理场推导的曲率模型在单元几何非线性分析中的有效性;

3)动力学分析还验证了模型三可有效缓解剪切闭锁问题,且在计算柔性梁剪切应变中收敛速度最快.

感谢美国伊利诺伊大学芝加哥分校Ahmed.A.Shabana教授的讨论.

[1]Bonnal C,Ruault J M,Desjean M C 2013 Acta Astronaut.85 51

[2]Nishida S I,Kawamoto S 2011 Acta Astronaut.68 113

[3]Liu J Y,Lu H 2007 Multibody Syst.Dyn.18 487

[4]He X S,Song M,Deng F Y 2011 Acta Phys.Sin.60 044501(in Chinese)[和兴锁,宋明,邓峰岩2011物理学报60 044501]

[5]He X S,Deng F Y,Wang R 2010 Acta Phys.Sin.59 1428(in Chinese)[和兴锁,邓峰岩,王睿2010物理学报59 1428]

[6]Chen S J,Zhang D G,Hong J Z 2013 Chin.J.Theor.Appl.Mech.45 251(in Chinese)[陈思佳,章定国,洪嘉振2013力学学报45 251]

[7]Shabana A A 1997 Multibody Syst.Dyn.1 189

[8]Tian Q,Zhang Y Q,Chen L P,Tan G 2010 Adv.Mech.40 189(in Chinese)[田强,张云清,陈立平,覃刚2010力学进展40 189]

[9]Omar M A,Shabana A A 2001 J.Sound Vib.243 565

[10]Hussein B A,Sugiyama H,Shabana A A 2007 J.Comput.Nonlinear Dyn.2 146

[11]Dmitrochenko O N,Hussein B A,Shabana A A 2009 J.Comput.Nonlinear Dyn.4 21002

[12]García-Vallejo D,Mikkola A M,Escalona J L 2007 Nonlinear Dyn.50 249

[13]Tian Q,Zhang Y Q,Chen L P,Yang J Z 2010 Nonlinear Dyn.60 489

[14]Gerstmayr J,Matikainen M K,Mikkola A M 2008 Multibody Syst.Dyn.20 359

[15]Nachbagauer K,Pechstein A S,Irschik H,Gerstmayr J 2011 Multibody Syst.Dyn.26 245

[16]Nachbagauer K,Gruber P,Gerstmayr J 2013 J.Comput.Nonlinear Dyn.8 021004

[17]Gerstmayr J,Shabana A A 2006 Nonlinear Dyn.45 109

[18]Dufva K E,Sopanen J T,Mikkola A M 2005 J.Sound Vib.280 719

[19]Mikkola A M,Dmitrochenko O,Matikainen M 2009 J.Comput.Nonlinear Dyn.4 011004

[20]Vesa-Ville A,Hurskainen T,Matikainen M K,Wang J,Mikkola A M 2016 J.Comput.Nonlinear Dyn.12 041007

[21]Zhang X S,Zhang D G,Chen S J,Hong J Z 2016 Acta Phys.Sin.64 094501(in Chinese)[章孝顺,章定国,陈思佳,洪嘉振2016物理学报64 094501]

[22]Goetz A 1970 Introduction to Di ff erential Geometry(Reading,Massachussetts:Addison Wesley Pub.Co)pp56–58

[23]Timoshenko S 1940 Strength of Materials(Part I Elementary Theory and Problems Second Edition)(New York:D.Van Nostrand Co)pp147–148

PACS:45.10.–b,05.45.–a,45.05.+xDOI:10.7498/aps.66.114501

Analysis of planar shear deformable beam using rotation fi eld curvature formulation∗

Zhang Da-Yu1)Luo Jian-Jun1)†Zheng Yin-Huan2)Yuan Jian-Ping1)
1)(National Key Laboratory of Aerospace Flight Dynamics,School of Astronautics,Northwestern Polytechnial University,Xi’an 710072,China)
2)(School of Mechanical and Electronic Engineering,Wuhan University of Technology,Wuhan 430070,China)

31 December 2016;revised manuscript

6 April 2017)

In recent years,research on space debris removal technique has received wide attention in aerospace fi eld.Many novel concepts on active fl exible debris remover have been proposed,such as fl exible fl ying net,tethered cable manipulator.In view with the high fl exibility and large deformation of this kind of structure,the implementation of attitude control is challenging.An accurate dynamic model of highly fl exible structure is important and needed.The beam element is the most common element adopted in fl exible remover models.So,in this investigation,a rotation fi eld-based curvature shear deformable beam using absolute nodal coordinate formulation(ANCF)(named RB-curvature ANCF beam)is proposed and its geometrically nonlinear characteristic under large deformation motion is studied.Curvature is fi rst derived through planar rotation transformation matrix between the reference frame and current tangent frame of beam centerline,and written as an arc-length derivative of the orientation angle of the tangent vector.Using the geometrically exact beam theory,the strain energy is expressed as an uncoupled form,and the new curvature is adopted to formulate bending energy.Based on the ANCF,the dynamic equation of beam is established,where mass and external force matrices are constant.In order to validate the performance of proposed beam element,other two types of beams are introduced as the comparative models.One is the classical ANCF fully parameterized shear deformable beam derived by continuum mechanics theory,and the other is position fi eld-based curvature ANCF shear deformable beam(named PB-curvature ANCF beam).The PB-curvature model is evaluated by di ff erentiating unit tangent vector of beam centerline with respect to its arc length quoted from di ff erential geometry theory.A series of static analysis,eigenfrequency tests and dynamic analysis are implemented to examine the performance of the proposed element.In static analysis,both small and non-small deformation cases show that the proposed RB-curvature ANCF beam achieves the faster speed,higher precision and good agreement with analytical solution in the case of cantilever beam subjected to a concentrated tip force,which is compared with other two beam models.The eigenfrequency analysis validates RB-curvature ANCF beam in a simply supported beam case that converges to its analytical solution.Meanwhile,the mode shapes of the proposed ANCF beam could be correctly corresponded to vibration state of element with respect to each di ff erent eigenfrequency.In the dynamics test,a fl exible pendulum case is used and simulation results show that the proposed RB-curvature ANCF beam accords well with ANSYS BEAM3,classical ANCF shear beam and PB-curvature ANCF beam in vertical displacements of tip point and middle point.Since deformation modes are uncoupled in the cross section of proposed beam element,its shear strain is achieved with much better convergence in the case of lower elastic modulus,and shear locking is signi fi cantly alleviated,compared with classical ANCF beam.Therefore,RB-curvature ANCF shear deformable beam element proposed in this paper is able to describe accurately geometric nonlinearity in large deformation problem,and can be a potential candidate in the modeling of fl exible/rigid- fl exible mechanisms.

absolute nodal coordinate formulation,rotation fi eld-based curvature,geometrically exact beam theory,shear locking

10.7498/aps.66.114501

∗国家自然科学基金重大项目(批准号:61690210,61690211)和国家自然科学基金(批准号:61603304,11472213)资助的课题.

†通信作者.E-mail:jjluo@nwpu.edu.cn

©2017中国物理学会Chinese Physical Society

http://wulixb.iphy.ac.cn

*Project supported by the Major Program of National Natural Science Foundation of China(Grant Nos.61690210,61690211),and the National Natural Science Foundation of China(Grant Nos.61603304,11472213).

†Corresponding author.E-mail:jjluo@nwpu.edu.cn