基于旋转场曲率的二维剪切梁单元建模-张大羽.pdf

《基于旋转场曲率的二维剪切梁单元建模-张大羽.pdf》由会员分享，可在线阅读，更多相关《基于旋转场曲率的二维剪切梁单元建模-张大羽.pdf（15页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、物理学报 Acta PhyssinV6166，No11(2017) 114501基于旋转场曲率的二维剪切梁单元建模冰张大羽1) 罗建军1)十 郑银环2) 袁建平1)1)(西北工业大学航天学院，航天飞行动力学技术重点实验室，西安 710072)2)(武汉理工大学机电工程学院，武汉430070)(2016年12月31日收到；2017年4月6日收到修改稿)对二维剪切梁单元进行研究，利用平面旋转场理论推导了精确曲率模型采用几何精确梁理论构建了剪切梁单元弹性力矩阵通过绝对节点坐标方法建立了系统的非线性动力学方程，提出基于旋转场曲率的二维剪切梁单元，并分别引入经典二维剪切梁单元和基于位移场曲率的二维剪切梁

2、单元进行比较研究首先，静力学分析证明了所提模型的正确性；其次，特征频率分析验证了模型可与理论解符合，收敛精度高，并且能准确地预测单元固有频率对应的振型；最后，在非线性动力学问题上，通过与ANSYS结果对比分析，证明了该模型可有效处理柔性大变形问题，并且与经典二维剪切梁单元相比具有缓解剪切闭锁的优势因此，本文提出的基于旋转场曲率的二维剪切梁单元在处理几何非线性问题中具有较大的应用潜力关键词：绝对节点坐标法，旋转场曲率，几何精确梁理论，剪切闭锁PACS：4510一b，0545一a，4505+x DoI：107498印s661145011引 言近年来，空间碎片清理技术在航天界受到广泛的关注1】世界主

3、要航天大国针对主动俘获碎片技术均投入了大量的人力财力进行研究，其中飞网、绳系机械爪等柔性俘获装置越来越多地出现在空间飞行器上由于该类俘获装置具有超轻、易变形及强几何非线性等特点，在执行抓捕碎片过程中自身会产生柔性大变形，使飞行器的姿态稳定及精确控制面临极大的挑战【2】梁单元是这类俘获装置结构中的基本单元针对其大变形问题，采用传统基于小变形假设的梁变形理论不能真实地反映大变形、大位移等特点，如运动弹性动力学分析、一次近似耦合模型【3、一次精确模型【4，5、高次耦合模型6】而Shabana7】提出了绝对节点坐标法fabs01ute nodal coordinate formulation，ANCF

4、)，该方法采用绝对位置坐标和变形梯度描述物体的刚体运动和柔性变形，适用于处理大变形、大转动多体系统的建模问题此外，该方法不采用转动角坐标，并且所有节点坐标均定义在惯性系中，因此构建出的质量矩阵为常数阵，不存在柯氏加速度项和离心加速度项，惯性力中也不存在变形和转动的附加耦合项进入21世纪以来，已有越来越多的学者关注ANCF方法的研究田强等8】对近年来ANCF方法的研究做了系统的评述在ANCF剪切梁单元研究中，Omar和Shabana9】基于连续介质力学理论最先提出了经典的二维ANCF全参数剪切梁有限元单元，该单元采用z向三次插值，!，向线性插值，这会导致严重的剪切闭锁问题，给数值计算带来困难随后

5、，很多学者针对该闭锁问题对此单元进行改进，更好地提升了单元的性能Hussein等【10】研究了二维ANCF全参数剪切梁单元中变形模态耦合情况，并采用HellingerReissner变分原理解决剪切闭锁问题随后，他们还将该研究思路应用到三维ANCF剪切梁单元和ANCF板单元11Garcia-、hllejo等【121提出了二维三节点ANCF剪国家自然科学基金重大项目(批准号：61690210，61690211)和国家自然科学基金(批准号：61603304，11472213)资助的课题十通信作者Email：jjluonwpueducn2017中国物理学会Chinese Physical Socie

6、ty 危ttp：叫ztzb咖九可nccn万方数据物理学报 Acta Pllyssin、厂0166，No11(2017) 114501切梁单元，不同于传统二维ANCF剪切梁单元的完备位移场，该单元的位移场是不完备的三次插值多项式，在此处理下剪切力呈二次分布并且梁截面在弯曲作用时不会发生收缩，解决了闭锁问题，动力学仿真验证了二维三节点ANCF剪切梁的有效性Tian等13采用Garcia一、hll ejo开发的10ckin牛free ANCF剪切梁单元进一步研究了考虑关节间隙的二维剪切梁动力学模型，丰富了绝对节点坐标法的单元类型Gerstmayr等f14基于绝对节点坐标法，采用Reissner非线性

7、梁理论重新构造了二维剪切梁单元的弹性力矩阵，公式中的轴向力、剪切力和弯曲力表述成显示解耦形式此外，他定义了厚度方向应变能用以表征横截面的变形，并准确地描述了剪切梁横截面固有频率对应的振型Gerstmayr定义这种弹性力解耦形式的梁单元为几何精确梁单元(geometrically exact beam element)，该方法对各类型应变能的积分只沿梁的中心线，各个变形模态之间互不影响，该处理可以有效缓解剪切闭锁Nachbagauer等15】只采用梁轴线上的梯度值r。I卵：o和h：o作为变形自由度，分别开发出线性和二次ANCF非全参数二维剪切梁单元由于新单元的梯度r，I订：o不能表征梁横截面的变

8、形，因此该处理成为解决剪切闭锁问题的一种新思路随后他们还开发出了三维ANCF剪切梁单元，解决了三维ANCF剪切梁单元存在的闭锁问题并显示出高阶收敛性16Gerstmay和Shabana【17针对ANCF剪切梁的闭锁现象又提出了两种解决方法：一种为在ANCF剪切梁单元基础上对其形函数中分别引入z2和z2z的插值项，开发了高阶三维ANCF梁单元，新增的插值项是截面关于轴向z的二次函数，这会使弯曲应变变成z的线性函数，并与截面对z的二阶导同阶，消除了剪切闭锁现象，更好地描述了截面的几何变形；另一种方法是将原有三维ANCF单元降阶并退化成三维低阶ANCF索单元，该单元中没有可方向梯度，相比原三维ANC

9、F梁单元的自由度数减少了一半，计算效率大幅提升，同时系统中也没有与剪切变形相关的刚度项，故剪切闭锁问题不存在Dufva等18通过截面转角描述二维ANCF剪切梁单元的精确位移场，梁轴向应变和弯曲应变由截面转角表征，而剪切变形采用单独的线性插值，该单元采用应变场混合插值技术解决了剪切闭锁问题Mikkola等19在二维ANCF非全参数梁的基础上通过定义截面方向上的力平衡条件引入剪切变形，开发出一种新的二维ANCF剪切梁单元新单元中因不含截面梯度自由度，计算效率远高于传统二维ANCF全参数剪切梁单元，并且也不存在数值计算时由截面梯度引起的高频问题他们提出的新单元与传统Timoshenko梁单元弯曲模态

10、和轴向模态符合度良好，证明了新单元不存在曲率闭锁问题VesaVille等20】在二维ANCF非全参数梁单元的基础上采用对梁可向变形场独立插值技术，提出了一种基于混合插值的二维ANCF剪切梁单元，节点坐标由位移场、截面方向和截面变形场组成，可捕捉到剪切变形动力学仿真结果表明该单元具有更高的精度和收敛性章孝顺等对二维ANCF非全参数梁单元的曲率模型进行讨论，对比了精确曲率模型和两种简化模型，同时还研究了不同简化模型的适用范围以上研究较全面地论述了二维ANCF剪切梁的发展历程，但针对二维剪切梁单元的曲率几何非线性进行的研究较少已有的曲率模型大都直接采用二阶曲率公式2引，定义在Frenet坐标系下，具

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

12、用性文中第二部分分别建立基于连续介质力学理论的经典剪切梁单元(模型一)、基于位移场曲率的几何精确剪切单元(模型二)和基于旋转场曲率的几何精确剪切梁单元(模型三)；第三部分引入小变形和非小变形静力学分析、特征频率分析以及非线性动力学分析验证所提单元的有效性；第四部分为结论1 145012万方数据物理学报 Acta PhysSinVbl66，No11(2017) 1145012剪切梁模型21模型一：基于连续介质力学理论的二维ANCF剪切梁单元二维含剪切变形的ANCF梁单元最早由Omar和Shabana提出9|，如图1所示，其中0一Xy为惯性坐标系，2为梁的长度对于二节点的ANCF剪切梁单元i，其梁

13、上任意一点在惯性系中的绝对位移矢量为止圈：l。o+。1z+n2+n3z可+。422+。523 【60+61z+62可+63z可+6422+6523 J=S。(z，可)e。， (1)图1剪切梁截面变形Fig1Cross section deformation of shear deformabiebeam其中z，可为梁初始构型在单元坐标系内的坐标该单元有12个绝对节点坐标，每个节点有6个自由度，包括2个位移自由度和4个梯度自由度：ei=P(zo)r哥(zo)飞tT(zo)riT(zt-z)r哥(zf)iT(zf)T (2)梁单元i的型函数9(z，可)为：s。(z，可)=l slj s2J s3J

14、 s4J s5J s6JI，其中J是22的单位矩阵s1=132+23，s2=f一22+3)，Is3=2(叩一叩)， s4=3223， (3)s5=f(一2+3)， s6=坳， I=zf，叩=可2利用连续介质力学中的变形理论推导得到ANCF剪切梁单元i的变形梯度张量Ji：n芸=等袅=2鬣篇篇卅1 l$。(z，可)e岛。(z，可)e I。 、对于梁单元的初始状态是水平摆放情况，山为单位阵GreenLagrangian应变张量的表达式为E=三(JTJ一J)=妒吖嚣1 ei警。卜卜sb=si!，1 sil可+s21，Tsi2可，sc=Si。1 Si+$。1$可依据Voigt标记规则，可将(5)式改写为

15、r 1 T。=lEi E；E；l， (6)其中，去(石1)1(e汀&ei一1)(打1)，E去(打1)1(e汀鼠ei一1)(打1)，E导去(打1)1(e汀&e)(打1)依据广义胡克定律，ANCF剪切梁单元i的应力张量盯i为盯。=E。， (7)E是材料的弹性系数矩阵在平面应变问题中(planestrain problem)，对于各向同性的匀质材料(isotropic homogenous material)，其弹性系数矩阵表达式为1 145013E=A+2“ A 0A A+2“0Lxyo毗一曲粥一曲堕粤丝如万方数据物理学报 Acta PhyssinVbl66，Non(2017) 114501入，肛

16、是拉梅常数(Lames constants)， E L， EA 2 T万而=丽p 2丽E，分别是材料的杨氏弹性模量和泊松比(6)式与(7)式可得到梁单元i的应变能为u=去胪：粕V则其弹性力矢量为=旧r：掣净，利用 其中Ir：I：(rr：)一，则(13)式可转换为恤=(筹)1=(，K。)T：P(A+2p)K1+A硷+4肛砥T(10)其中，K-=去二(石1)T文(eT&et一，)+鼠(eTsbei一1)(打1)dK硷=三二(石1)T良(&e一1)+b(eiTset一1)(打1)dV凰=三二(石1)T袅(etTscei)(靠1)dV良=&+镑，b=鼠+s苫，曳=&+器 (11)经典剪切梁单元模型中轴

17、向伸长变形、剪切变形和弯曲变形相互耦合在一起，因此曲率公式不是显式形式22模型二：基于位移场曲率的二维几何精确剪切梁单元根据微分几何理论【22】，曲率是单位切向量对弧长的一阶导矢的模本例中切向量是梁中心线上的位移梯度，通过位移场推导得到的曲率模型，本文定义为位移场曲率图1中梁变形后其截面不再与梁的中心轴线s垂直根据梁位移场，梁单元i上任意一点变形后的位移为：ri=r(z，可)其对中心轴线s的一阶导为r：=等=筹塞=禹， (2)r丽2面丽2商， 对中心线s的二阶导为azas吒=r：驾掣净=(高“(警膳=(高“(莆)南一二奠一(!兰!；立!；一r铲r：(r轳r：)2一二量一(!；!兰2 1；兰一r

18、耋Tr：(rr：)2=(麦(J_)吒， (14)o代表张量积或外积，曲率为Ki。=Ir；。1依据几何精确梁理论，梁单元i的应变势能可由沿中心轴线s的轴向应变能叱，剪切应变能眈和弯曲应变能阮组成：U2=瞻+职+魄=三EA，i2ds+三GA。，12ds+吉肛J(Ki。)2ds， (15)E，G，A分别是梁单元i的弹性模量、剪切模量和横截面积；剪切模量G=E2(1+)；A。=七。A，后。是剪切修正因子，对于横截面是矩形的情况，忌。=10(1+)12+11；应变，沿梁中心线的分量砰，曩可由拉格朗日应变张量确定：日=去(r孙。Trk。一)+三(r产r；一)， (16)砭=r：|唧Tr； (17)梁单元i

19、的弹性力公式可由应变势能Ui对节点坐标ei求偏导得到：鸥础；。=警=警+警+警，(18)其中，警=EA砰筹ds掣=丝c3丝如=丝c3=吃万方数据物理学报 Acta PhyssinV0166，No11(2017) 114501+去EA(日) 2a(ds)等=ZG们等ds+三GA。(砭)ae2a(ds)警=哦。警ds+丢厂Em二t，Zds=Ir：I dzaei 2a(ds)应变分量砰，砭对节点坐标ei求导为a毋aea疋aeoTr扎：o+卵r；，(19)=即r孙：o+I唧Tr； (20)曲率K：。对节点坐标ei求导为：其中，aK；。aeiarkaeiar：。ar：z=撕。Tr：。)冯警，一ar：。a

20、r三。 ar：。aei：薯或。 ar：。(Jr：o r：)f蟛。7：l7猡r：ar量一a(rr：)ar： a7：j_立竖r护r：Ir：I fr：l(r：r争)(r铲r：)21一(Jr；o r；)r：。+(rr：)一1望掣r：。：掣(Jr：o r城ar兰 、 。8”8。一(rr：)一1皇j三；笋r：。一 2r：。(7r：)2 冉糍等一嵩(r字r：)。2r舞r：r：r字(r字r：)3一高冉鬻等(21)(22)(23)(24)扛高一=z。达式推导如下：昙(r笋r：)二z。T(19)式中的掣表z2 1(19)式中的等表a(fr：I dz)aeia(rr：)：剑dz aet 虻鲁h (25)23模型三：

21、基于旋转场曲率的二维几何精确剪切梁单元曲率也可定义为曲线上某点的切线方向角对弧长的转动率通过角度场推导得到的曲率模型，本文定义为旋转场曲率由图1所示，根据平面旋转场理论，由DXy逆时针转动到移动矢量r f。，易的转动矩阵为l C08Qslnal，其中Q是。 。【由眦蛐Jr：与ox轴的夹角，即梁中线的切线方向角r：沿着梁中心线切线方向，可体现梁中线弯曲效应则单位切向量可写成皓高=陋咖QT-胁T (26)由(26)式可知：cosQ=岔l，由该式可得到旋转场曲率公式：i aQ 1 a岔】 1 a岔】azK2丽2一吾而丽2一面面丽：一熹羔嘉， (27)sina a。l矿：I 、“。7其中要是姿的第一项

22、，娑表达式为 oz oz dz狮：az=a呓加：ar：az(卜表怕r： ) r (28)该模型中的轴向势能与剪切势能与模型二中的表达式相同，弯曲势能可利用旋转场曲率得到：吮=丢EJK：ds弯曲势能对节点坐标e的导数为警=胁：警as+去胁尝掣(29)其中等推导如下11450l一5等=暑(警)=妄(筹)一=一Il=一一Iaei ae、as asae=妄(熹蓦) assinQ ae丝钯监哦盟睨彘。芦舞万方数据物理学报 Acta PhyssinVbl66，No11(2017) 114501=杀(未)鲁一熹杀(鲁)邓。，(30)式中的偏导项分别为：岳(纛)=嘉(未)警=嘉(熹)南一一cosQa岔1一一岔

23、1a庐。(1)2 i石面2虿11fj数阵，不包含科氏力与离心力，不存在变形和转动耦合附加的惯性力通过虚功原理可推导单元的外力阵公式：Q。t。，。l=s(z，可)TF， (36)F为外力矢量，外力阵同样为常数矩阵(31) 3算例分析dzl r兰【lJae aeta庐： a庐：a?：一一一aei ar曼aei=I志(卜表曲r：)卜，晏(蓦)=未(掣)塞号毫(豢)南=吴I志(卜彘枷：)l高=是辟警az Ir：I Ir：13=恚矧一嘉l警I=南s。击静：)-智。 +紫一监昔塑或)(33)同理，通过(26)式的第二式岔2=sinQ也可推导出旋转场曲率K：的表达式根据绝对节点坐标方法，由虚功原理可推导梁M

24、。芭。=Q。xte。lQ：l。ti。， (34)Q：l。i。分别采用模型一至三中的弹性力表达式，百是节点加速度坐标对于ANCF剪切梁单元i，其单元质量阵公式为， 一M2=p。s。(z，可)1Js。(z，可)dy2， (35)Vt其中，yi是梁单元i的密度和体积由(35)式可知，通过绝对节点坐标方法建立的单元质量阵为常31静力学分析本节引入典型的悬臂梁受集中载荷问题，考察所提出基于旋转场曲率的ANCF剪切梁模型在静力学分析中的准确性图1所示为悬臂梁受集中载荷模型，其物理参数与文献f141相同长为2=2 m，宽为训=01 m，高为=05 m，杨氏模量为E=207105 GPa，剪切模量为G=796

25、104 GPa，泊松比为=03，以及剪切修正因子为：桨竺兰h口H 2悬臂梁模型Fig2Cantilever beamF考察小变形算例，设置末端集中力为F=一5105九3 N在小变形假设下，梁中性层在变形后长度不变，端点沿X方向的变形可以忽略不计y向挠度的理论解23为6=熹，其中L为梁截面惯性矩通过表1的数据对比发现，由于模型一存在严重的剪切闭锁问题，其y向位移与理论解误差达到22而模型二和模型三采用几何精确梁理论构建系统弹性力，其中轴向力、剪切力和弯曲力相互解耦，可缓解剪切闭锁，因此由它们计算的y向位移均可收敛于理论解，并且收敛速度快减小梁的高度至=01 m，末端集中力为F=一1106九3 N

26、，其他参数不变表2结果表明模型一的误差为25，而模型二和三可收敛到理论解与表1相比，表2中变形量增大，模型二和三需要收敛到理论解n，J个数增多万方数据物理学报 Acta PllyssinVbl66，No11(2017) 114501表1 小变形情况下三种弹性力模型计算的端部位移对比1rable 118comparison of tip point displacement using three elaLstic force models in case ofsmaU defbrmation表2 小变形情况下三种弹性力模型计算的端部位移对比2rable 22nd comparison of t

27、ip point displacement using three elaLstic force models in case ofsmall defbrmation表3 非小变形情况下三种弹性力模型计算的端部位移对比Table 3 Comparison of tip point displacement using three ela8tic force models in case ofnonsmall deformation划分单元数个 端点x向位移m 端点y向位移m连续介质力学模型(模型一)1 0049726 0397010100 一O095631 0562668300 0095648

28、 0562718几何精确梁模型位移场曲率(模型二)1 0070542 047479610 一0149856 0707640几何精确梁模型一旋转场曲率(模型三11 0078069 04858595 0 150970 0 700624参考解【14 一O150971 07105691 145017万方数据物理学报 Acta PhyssinVbl66，No11(2017) 114501考察非小变形情况，设置末端集中力为F=一51083 N，其他参数与小变形算例1相同采用该集中力可导致相对大的变形量采用文献f141表2中用商业软件计算的解作为非线性大变形下悬臂梁x，y向挠度的参考解在大变形情况下，剪切效

29、应加剧，另外由于ANCF单元在y向上采用线性插值会导致严重的剪切闭锁问题，所以在大变形算例中三种模型的收敛精度均不如小变形算例，见表3数据对比显示模型三的收敛速度和精度最高以上算例证明了所提出的模型在静力学分析中的正确性32特征频率分析本节对所提出的ANCF剪切梁单元的特征频率进行分析，选择具有理论解的简支梁作为研究对象简支梁物理参数与文献f14中相同其中长为f=2 m，宽为叫=04 m，高为=04 m，密度为p一7850埏m3，杨氏模量为E=1000 GPa，剪切模量为G=38462 GPa，泊松比为=03，以及剪切修正因子为七。：芸掣本算例中的边界条件如图3所示，左端的X，y方向位移固定，

30、右端只固定y方向位移图3简支梁模型Fig3Simply supported beam model图4为本文模型三在简支边界条件下的固有频率和振型因为通过边界条件消除了3个自由度，蕊础疑； ： 一： 衄露爨r： 。： i：mode l：100051 rad8 mode 2：280321 rad8 mode 3：381294 rad8睢： ：；l；!： ： ：】； l【l ll l堋 til|lii|1ii|11|ll|llll|li：l【l【|l l【l l 】ji|tl|1|lltt|lti|l!l? L_1 _ l【l ll l|111|lllll】【|lt|lllt|l|11mode 4：

31、840962 rad8 mode 5：140160 rads mode 6：1796 99 rads，I翮川川蚶心、? 衄 丁 砷研碰必必倒删14J“Um、3盐j雌 衄山uJJIllllljl、 k30 098 radsmode 7：187887 rads mode 9：图4 基于旋转场曲率的简支梁模型的振型和固有频率Fig4】ode shapes and eigenfrequencies(rads)of simply supported beam based on rotationbased curvature表4 3种弹性力模型下简支梁的弯曲和轴向固有频率(rads)对比Iable 4Be

32、nding and longitudinal eigenfequencies(rads)of simply supported beam element fbr differentelastic formulation万方数据物理学报 Acta PhyssinV6166，No11(2017) 114501表5 3种弹性力模型下简支梁的剪切和厚度固有频率(rads)对EE工1able 5 shearing and thickness eigenfequencies(rads)of simply supported beam element for di矗erentelastic fbrmulat

33、ion所以共计算出有9种振型和对应的固有频率ANCF单元的轴向采用三次插值技术产生三种轴向振型(mode 2，mode 4和mode 5)和两种弯曲振型(mode 1和mode 31另外，还存在两种剪切振型fmode 6和mode 7)以及ANCF单元特有的厚度振型fmode 8和mode 9)表4为简支梁单元分别采用3种弹性力模型所计算出的弯曲和轴向固有频率的对比，表5为其剪切和厚度固有频率对比由结果可知，与其他两个模型相比，本文提出的模型三收敛性更好，收敛精度可至6位有效数字引起二阶弯曲振动频率误差较高是因为本文单元在y向上采用线性插值，会导致单元的剪切闭锁问题，使得变形过程中剪切应变能占

34、据主导，过多的剪切力导致单元不能发生正常弯曲变形，造成弯曲刚度增大，单元变得更“刚”，从而引起二阶弯曲振动频率增大33动力学分析图5所示为柔性单摆，0点为旋转铰该算例是绝对节点坐标方法的一个经典算例单摆的物理参数如下：长为f=12 m，横截面积为A=00016 m2，截面惯性矩为j：=8533105 m4，密度为JD=5540 kgm3，杨氏模量为E=7 GPa，剪切模量为G=269 GPa，泊松比以及剪切修正因子与前面算例相同如图5所示，仅固定左端的X，y方向位移考察在重力作用下单摆的动态响应仿真中单摆被划分为12个单元，仿真时间为11 s本算例中柔性单摆在其重力作用下进行大范围转动，运动过

35、程中发生大变形现象为对比文中所述的三种剪切梁单元模型处理大变形问题的正确性，采用ANSYS中两节点平面梁BEAM3单元作为对比由图6一图8可见、在弹性模量为E=7 GPa时，三种模型计算的末端横向位移，中点横向位移和末端端点轨迹均与ANSYS结果符合，证明它们均能处理大变形问题，并验证了本文所提模型三的正确性图5柔性单摆模型F培5Flexible pendulum为了更好地说明模型三处理柔性大变形的能力，将弹性模量降低至E=2 GPa在此弹性模量下，梁将变得更柔，大变形现象更明显图9一图ll显示在E=2 GPa时，三种模型计算的结果均与ANSYS结果符合继续将弹性模量降至E=07 GPa，由图

36、12可见，在此低模量下三种模型计算的结果也可与ANSYS结果相近以上算例均证明了模型三可用于处理柔性大变形问题1 145019万方数据物理学报 Acta PhyssinVbl66，No11(2017) 114501吕淤喜差强002O4O60 81 01 214O 0 2 0 4 0 6 0 8 1 0时间s图6柔性单摆末端横向位移(E=7 GPa)Fig6 The tip vertical displacement of flexible pendulum(E=7 GPa)O0 102030 40，5060 70 0 2 O 4 O6 08 10时间s图7 柔性单摆中点横向位移(E=7 GPa

37、)Fig7 The middle vertical displacement of flexiblependulum(E=7 GPa)0020 4吕一O6孓o81 0121 4l 5 10 O5 0 05 1 0 15Xm图8柔性单摆末端位移(E=7 GPa)Fig 8 TlItj tiI)-p)il“trajectoryJf flexible pelHllllIll(E=7 GPa)11450110时间s图9柔性单摆末端横向位移(E=2 GPa)Fig9 The tip vertical displacement of nexible pendulum(E=2 GPa)时间s图10柔性单摆中

38、点横向位移(E=2 GPa)Fig10 The middle vertical displacement of flexiblependulum(E=2 GPa)O20O 2一O 4要-oe一081 O一121 41 5 1O O5 0 0 5 1O 15Xm图11柔性单摆末端位移(E=2 GPa)Fig 1 1The tipI)oillt trajectory of nexil】ltl lJtlII(1lJluIl(E=2 GPa)；饕基墨颦箍磐，。渣迥苣颦最碧0饕鬯量颦箍碧万方数据物理学报 Acta Physsin、厂0166，No11(2017) 114501图12柔性单摆末端位移(E=

39、07GPa)Fig12The tippoint trajectory of丑exible pendulum(E=07 GPa)图13一图15是模型三在不同弹性模量下的能量变化曲线由图可知，该模型总能量为零，符合能量守恒定律从能量的角度验证了模型三的正确性 ；一l一iO 01 02 0 3 04 O5 O 6 07 O8 09 10 11时间s图13 系统能量分布(旋转场曲率模型，E=7 GPa)Fig13 The energy distribution of nexible pendulum(RB curvature beam model，E=7 GPa)；一一；鸶i| ；0 O1 02 03

40、 O 4 05 O6 07 08 09 10 11时间s图14 系统能量分布(旋转场曲率模型，E=2 GPa)Fig14The energy distribution of flexible pendulum(RB curVature beam model，E=2 GPa)；。一： ：；囊 O O1 02 03 O4 O5 0 6 07 08 O 9 1O 11时间s图15 系统能量分布(旋转场曲率模型，E=07 GPa)Fig15The energy distribution of flexible pendulum(RB curVature beam model，E=07 GPa)图16一

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

42、性图16 (网刊彩色)柔性单摆各时刻构型图(E=7 GPa)(绿色实线，模型一；蓝色点画线，模型二；红色虚线，模型三)Fig16(color online)The conflguration of flexible pendulum at each giVen time(E=7 GPa)(green solid line，18model；blue dash dotted Iine，2“d model；red dashline：3。d moden加如如加。加如如加踟伯加mo加们鲫加鲫万方数据物理学报 Acta PllyssinVbl66，No11(2017) 114501图17 (网刊彩色)柔性

43、单摆各时刻构型图(E=2 GPa)(网刊彩色)(绿色实线：模型一；蓝色点画线：模型二；红色虚线：模型三)Fig17 (color online)The configuration of flexiblependulum at each given time(E=2 GPa)(green solid1ine，18。model；blue dash dotted line，2“d model；reddas hline，3。8 model)图18 (网刊彩色)柔性单摆各时刻构型图(E=07 GPa)(绿色实线，模型一；蓝色点画线，模型二；红色虚线，模型三1Fig18 (color online)The

44、 configuration of flexiblependulum at each giVen time(E=07 GPa)(greensolid line，18。model；blue dash dotted line，2“d model；red dash line：3。d modell柔性梁在大转动过程中剪切效应逐渐占据主导，为了进一步研究三种不同的弹性力计算的剪切应变在大转动大变形中的变化，弹性模量选取E=07 GPa，其余梁的物理参数不变，柔性单摆被划分为40个单元图19是单摆最后一个单元内中点的剪切应变变化曲线从图中可知，模型二与模型三得到的剪切应变相近，而它们均比模型一的剪切应变大

45、，其最大差异可达到001其差异可由图20和图2l解释，分别是单摆最后一个单元中点的r。，7，方向与大小变化曲线由图可知，模型二与模型三在r。，方向变化上与模型一基本相似，而，的幅值变化与模型一有较大差异，这也是导致图19剪切应变差异的主要原因绝对节点坐标法中的应变是位移梯度的函数，例如正应力表达式 m 1 mEzz=石rzl rz一1， E可可=百r可上r可一1，1剪应力表达式。可=去r。T故图20和图21中的单元内位移梯度随时间的变化也反映了相应应变在大转动过程中的变化情况由图17 fb)中可知，模型二和模型三的n数值被限制在1附近，变化范围较小，因为这两种模型的弹性能中具有沿单元截面方向的

46、能量1 r 专EA(r可Tr可一1)2ds，-J它可限制剧烈变化，确保单元的不可压缩性模型一由于是耦合模型，与变化相互影响，所以其7。幅值的变化范围比模型二和模型三要更明显图19单摆最后一个单元内中点的剪切应变变化Fig19Change of shear strain at middle point of lastelement由于ANCF剪切梁单元在z，方向上插值阶次的不同将会导致严重的剪切闭锁现象，采用本文模型一至模型三分别研究各模型缓解剪切闭锁的能力表6为分别采用三种模型时的单摆中点剪切应变收敛情况，可以看出采用解耦模型f模型二和模型三)得到的剪应变收敛个数比耦合模型(模型一)要少，模型三收敛最快，证明模型三在缓解剪切闭锁问题上具有较强的优势模型一由于各个变形模态相互耦合，剪切应变收敛最慢1 1450112O8642O2468O1000OOOOO10O000000000000OO00O0一一一一一制型尽豁万方数据物理学报 Acta PhyssinV0166，No11(2017) 114501200180160140120，lOOk80604020020O 01 O2 0 3 0 4 O 5 O6 O7 08 O9 1O 1 1时间s图20单摆最后一个单元内中点的r。和7l方向变化 (a)r。；(b)rFig20orientation of the rz and rv