張大羽羅建軍?鄭銀環(huán)袁建平

1)(西北工業(yè)大學(xué)航天學(xué)院,航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室,西安 710072)

2)(武漢理工大學(xué)機(jī)電工程學(xué)院,武漢430070)

基于旋轉(zhuǎn)場(chǎng)曲率的二維剪切梁?jiǎn)卧?

張大羽1)羅建軍1)?鄭銀環(huán)2)袁建平1)

1)(西北工業(yè)大學(xué)航天學(xué)院,航天飛行動(dòng)力學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室,西安 710072)

2)(武漢理工大學(xué)機(jī)電工程學(xué)院,武漢430070)

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

對(duì)二維剪切梁?jiǎn)卧M(jìn)行研究,利用平面旋轉(zhuǎn)場(chǎng)理論推導(dǎo)了精確曲率模型.采用幾何精確梁理論構(gòu)建了剪切梁?jiǎn)卧獜椥粤仃?通過(guò)絕對(duì)節(jié)點(diǎn)坐標(biāo)方法建立了系統(tǒng)的非線性動(dòng)力學(xué)方程,提出基于旋轉(zhuǎn)場(chǎng)曲率的二維剪切梁?jiǎn)卧?并分別引入經(jīng)典二維剪切梁?jiǎn)卧突谖灰茍?chǎng)曲率的二維剪切梁?jiǎn)卧M(jìn)行比較研究.首先,靜力學(xué)分析證明了所提模型的正確性;其次,特征頻率分析驗(yàn)證了模型可與理論解符合,收斂精度高,并且能準(zhǔn)確地預(yù)測(cè)單元固有頻率對(duì)應(yīng)的振型;最后,在非線性動(dòng)力學(xué)問(wèn)題上,通過(guò)與ANSYS結(jié)果對(duì)比分析,證明了該模型可有效處理柔性大變形問(wèn)題,并且與經(jīng)典二維剪切梁?jiǎn)卧啾染哂芯徑饧羟虚]鎖的優(yōu)勢(shì).因此,本文提出的基于旋轉(zhuǎn)場(chǎng)曲率的二維剪切梁?jiǎn)卧谔幚韼缀畏蔷€性問(wèn)題中具有較大的應(yīng)用潛力.

絕對(duì)節(jié)點(diǎn)坐標(biāo)法,旋轉(zhuǎn)場(chǎng)曲率,幾何精確梁理論,剪切閉鎖

1 引言

近年來(lái),空間碎片清理技術(shù)在航天界受到廣泛的關(guān)注[1].世界主要航天大國(guó)針對(duì)主動(dòng)俘獲碎片技術(shù)均投入了大量的人力財(cái)力進(jìn)行研究,其中飛網(wǎng)、繩系機(jī)械爪等柔性俘獲裝置越來(lái)越多地出現(xiàn)在空間飛行器上.由于該類俘獲裝置具有超輕、易變形及強(qiáng)幾何非線性等特點(diǎn),在執(zhí)行抓捕碎片過(guò)程中自身會(huì)產(chǎn)生柔性大變形,使飛行器的姿態(tài)穩(wěn)定及精確控制面臨極大的挑戰(zhàn)[2].梁?jiǎn)卧沁@類俘獲裝置結(jié)構(gòu)中的基本單元.針對(duì)其大變形問(wèn)題,采用傳統(tǒng)基于小變形假設(shè)的梁變形理論不能真實(shí)地反映大變形、大位移等特點(diǎn),如運(yùn)動(dòng)彈性動(dòng)力學(xué)分析、一次近似耦合模型[3]、一次精確模型[4,5]、高次耦合模型[6].而Shabana[7]提出了絕對(duì)節(jié)點(diǎn)坐標(biāo)法(absolute nodal coordinate formulation,ANCF),該方法采用絕對(duì)位置坐標(biāo)和變形梯度描述物體的剛體運(yùn)動(dòng)和柔性變形,適用于處理大變形、大轉(zhuǎn)動(dòng)多體系統(tǒng)的建模問(wèn)題.此外,該方法不采用轉(zhuǎn)動(dòng)角坐標(biāo),并且所有節(jié)點(diǎn)坐標(biāo)均定義在慣性系中,因此構(gòu)建出的質(zhì)量矩陣為常數(shù)陣,不存在柯氏加速度項(xiàng)和離心加速度項(xiàng),慣性力中也不存在變形和轉(zhuǎn)動(dòng)的附加耦合項(xiàng).

進(jìn)入21世紀(jì)以來(lái),已有越來(lái)越多的學(xué)者關(guān)注ANCF方法的研究.田強(qiáng)等[8]對(duì)近年來(lái)ANCF方法的研究做了系統(tǒng)的評(píng)述.在ANCF剪切梁?jiǎn)卧芯恐?Omar和Shabana[9]基于連續(xù)介質(zhì)力學(xué)理論最先提出了經(jīng)典的二維ANCF全參數(shù)剪切梁有限元單元,該單元采用x向三次插值,y向線性插值,這會(huì)導(dǎo)致嚴(yán)重的剪切閉鎖問(wèn)題,給數(shù)值計(jì)算帶來(lái)困難.隨后,很多學(xué)者針對(duì)該閉鎖問(wèn)題對(duì)此單元進(jìn)行改進(jìn),更好地提升了單元的性能.Hussein等[10]研究了二維ANCF全參數(shù)剪切梁?jiǎn)卧凶冃文B(tài)耦合情況,并采用Hellinger-Reissner變分原理解決剪切閉鎖問(wèn)題.隨后,他們還將該研究思路應(yīng)用到三維ANCF剪切梁?jiǎn)卧虯NCF板單元[11].García-Vallejo等[12]提出了二維三節(jié)點(diǎn)ANCF剪切梁?jiǎn)卧?不同于傳統(tǒng)二維ANCF剪切梁?jiǎn)卧耐陚湮灰茍?chǎng),該單元的位移場(chǎng)是不完備的三次插值多項(xiàng)式,在此處理下剪切力呈二次分布并且梁截面在彎曲作用時(shí)不會(huì)發(fā)生收縮,解決了閉鎖問(wèn)題,動(dòng)力學(xué)仿真驗(yàn)證了二維三節(jié)點(diǎn)ANCF剪切梁的有效性.Tian等[13]采用García-Vallejo開(kāi)發(fā)的locking-free ANCF剪切梁?jiǎn)卧M(jìn)一步研究了考慮關(guān)節(jié)間隙的二維剪切梁動(dòng)力學(xué)模型,豐富了絕對(duì)節(jié)點(diǎn)坐標(biāo)法的單元類型.Gerstmayr等[14]基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法,采用Reissner非線性梁理論重新構(gòu)造了二維剪切梁?jiǎn)卧膹椥粤仃?公式中的軸向力、剪切力和彎曲力表述成顯示解耦形式.此外,他定義了厚度方向應(yīng)變能用以表征橫截面的變形,并準(zhǔn)確地描述了剪切梁橫截面固有頻率對(duì)應(yīng)的振型.Gerstmayr定義這種彈性力解耦形式的梁?jiǎn)卧獮閹缀尉_梁?jiǎn)卧?geometrically exact beam element),該方法對(duì)各類型應(yīng)變能的積分只沿梁的中心線,各個(gè)變形模態(tài)之間互不影響,該處理可以有效緩解剪切閉鎖.Nachbagauer等[15]只采用梁軸線上的梯度值rx|η=0和ry|η=0作為變形自由度,分別開(kāi)發(fā)出線性和二次ANCF非全參數(shù)二維剪切梁?jiǎn)卧?由于新單元的梯度ry|η=0不能表征梁橫截面的變形,因此該處理成為解決剪切閉鎖問(wèn)題的一種新思路.隨后他們還開(kāi)發(fā)出了三維ANCF剪切梁?jiǎn)卧?解決了三維ANCF剪切梁?jiǎn)卧嬖诘拈]鎖問(wèn)題并顯示出高階收斂性[16].Gerstmay和Shabana[17]針對(duì)ANCF剪切梁的閉鎖現(xiàn)象又提出了兩種解決方法:一種為在ANCF剪切梁?jiǎn)卧A(chǔ)上對(duì)其形函數(shù)中分別引入x2y和x2z的插值項(xiàng),開(kāi)發(fā)了高階三維ANCF梁?jiǎn)卧?新增的插值項(xiàng)是截面關(guān)于軸向x的二次函數(shù),這會(huì)使彎曲應(yīng)變變成x的線性函數(shù),并與截面對(duì)x的二階導(dǎo)同階,消除了剪切閉鎖現(xiàn)象,更好地描述了截面的幾何變形;另一種方法是將原有三維ANCF單元降階并退化成三維低階ANCF索單元,該單元中沒(méi)有y方向梯度,相比原三維ANCF梁?jiǎn)卧淖杂啥葦?shù)減少了一半,計(jì)算效率大幅提升,同時(shí)系統(tǒng)中也沒(méi)有與剪切變形相關(guān)的剛度項(xiàng),故剪切閉鎖問(wèn)題不存在.Dufva等[18]通過(guò)截面轉(zhuǎn)角描述二維ANCF剪切梁?jiǎn)卧木_位移場(chǎng),梁軸向應(yīng)變和彎曲應(yīng)變由截面轉(zhuǎn)角表征,而剪切變形采用單獨(dú)的線性插值,該單元采用應(yīng)變場(chǎng)混合插值技術(shù)解決了剪切閉鎖問(wèn)題.Mikkola等[19]在二維ANCF非全參數(shù)梁的基礎(chǔ)上通過(guò)定義截面方向上的力平衡條件引入剪切變形,開(kāi)發(fā)出一種新的二維ANCF剪切梁?jiǎn)卧?新單元中因不含截面梯度自由度,計(jì)算效率遠(yuǎn)高于傳統(tǒng)二維ANCF全參數(shù)剪切梁?jiǎn)卧?并且也不存在數(shù)值計(jì)算時(shí)由截面梯度引起的高頻問(wèn)題.他們提出的新單元與傳統(tǒng)Timoshenko梁?jiǎn)卧獜澢B(tài)和軸向模態(tài)符合度良好,證明了新單元不存在曲率閉鎖問(wèn)題.Vesa-Ville等[20]在二維ANCF非全參數(shù)梁?jiǎn)卧幕A(chǔ)上采用對(duì)梁y向變形場(chǎng)獨(dú)立插值技術(shù),提出了一種基于混合插值的二維ANCF剪切梁?jiǎn)卧?節(jié)點(diǎn)坐標(biāo)由位移場(chǎng)、截面方向和截面變形場(chǎng)組成,可捕捉到剪切變形.動(dòng)力學(xué)仿真結(jié)果表明該單元具有更高的精度和收斂性.章孝順等[21]對(duì)二維ANCF非全參數(shù)梁?jiǎn)卧那誓Ｐ瓦M(jìn)行討論,對(duì)比了精確曲率模型和兩種簡(jiǎn)化模型,同時(shí)還研究了不同簡(jiǎn)化模型的適用范圍.以上研究較全面地論述了二維ANCF剪切梁的發(fā)展歷程,但針對(duì)二維剪切梁?jiǎn)卧那蕩缀畏蔷€性進(jìn)行的研究較少.已有的曲率模型大都直接采用二階曲率公式[22],定義在Frenet坐標(biāo)系下,具有極值性;且鮮有從旋轉(zhuǎn)場(chǎng)角度研究曲率定義,并對(duì)基于旋轉(zhuǎn)場(chǎng)曲率的二維ANCF剪切梁?jiǎn)卧M(jìn)行幾何非線性分析、討論模型的正確性和適用性.

本文進(jìn)一步針對(duì)二維ANCF剪切梁的曲率定義進(jìn)行研究.基于平面旋轉(zhuǎn)場(chǎng)理論推導(dǎo)得到精確曲率公式,表征為梁中心線切線方向轉(zhuǎn)角的函數(shù);采用幾何精確梁理論構(gòu)建系統(tǒng)彈性力矩陣,提出基于旋轉(zhuǎn)場(chǎng)曲率的二維ANCF剪切梁?jiǎn)卧Ｐ?同時(shí)引入基于連續(xù)介質(zhì)力學(xué)理論的經(jīng)典ANCF剪切梁?jiǎn)卧突谖⒎謳缀卫碚摰奈灰茍?chǎng)曲率ANCF剪切梁?jiǎn)卧鳛閷?duì)比研究.分別通過(guò)靜力學(xué)分析、特征頻率分析以及動(dòng)力學(xué)分析驗(yàn)證了本文所提的基于旋轉(zhuǎn)場(chǎng)曲率的二維ANCF剪切梁?jiǎn)卧Ｐ驮谟?jì)算柔性機(jī)構(gòu)大變形、大轉(zhuǎn)動(dòng)問(wèn)題中的準(zhǔn)確性和適用性.

文中第二部分分別建立基于連續(xù)介質(zhì)力學(xué)理論的經(jīng)典剪切梁?jiǎn)卧?模型一)、基于位移場(chǎng)曲率的幾何精確剪切單元(模型二)和基于旋轉(zhuǎn)場(chǎng)曲率的幾何精確剪切梁?jiǎn)卧?模型三);第三部分引入小變形和非小變形靜力學(xué)分析、特征頻率分析以及非線性動(dòng)力學(xué)分析驗(yàn)證所提單元的有效性;第四部分為結(jié)論.

2 剪切梁模型

2.1 模型一:基于連續(xù)介質(zhì)力學(xué)理論的二維ANCF剪切梁?jiǎn)卧?/p>

二維含剪切變形的ANCF梁?jiǎn)卧钤缬蒓mar和Shabana提出[9],如圖1所示,其中O-XY為慣性坐標(biāo)系,l為梁的長(zhǎng)度.

對(duì)于二節(jié)點(diǎn)的ANCF剪切梁?jiǎn)卧猧,其梁上任意一點(diǎn)在慣性系中的絕對(duì)位移矢量為

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

其中x,y為梁初始構(gòu)型在單元坐標(biāo)系內(nèi)的坐標(biāo).該單元有12個(gè)絕對(duì)節(jié)點(diǎn)坐標(biāo),每個(gè)節(jié)點(diǎn)有6個(gè)自由度,包括2個(gè)位移自由度和4個(gè)梯度自由度:

梁?jiǎn)卧猧的型函數(shù)Si(x,y)為:

其中I是2×2的單位矩陣.

利用連續(xù)介質(zhì)力學(xué)中的變形理論推導(dǎo)得到ANCF剪切梁?jiǎn)卧猧的變形梯度張量Ji:

對(duì)于梁?jiǎn)卧某跏紶顟B(tài)是水平擺放情況,J0為單位陣.Green-Lagrangian應(yīng)變張量εi的表達(dá)式為

依據(jù)Voigt標(biāo)記規(guī)則,可將(5)式改寫為

依據(jù)廣義胡克定律,ANCF剪切梁?jiǎn)卧猧的應(yīng)力張量σi為

E是材料的彈性系數(shù)矩陣.在平面應(yīng)變問(wèn)題中(plane-strain problem),對(duì)于各向同性的勻質(zhì)材料(isotropic homogenous material),其彈性系數(shù)矩陣表達(dá)式為

λ,μ是拉梅常數(shù)(Lame’s constants),

E,ν分別是材料的楊氏彈性模量和泊松比.利用(6)式與(7)式可得到梁?jiǎn)卧猧的應(yīng)變能為

則其彈性力矢量為

其中,

經(jīng)典剪切梁?jiǎn)卧Ｐ椭休S向伸長(zhǎng)變形、剪切變形和彎曲變形相互耦合在一起,因此曲率公式不是顯式形式.

2.2 模型二:基于位移場(chǎng)曲率的二維幾何精確剪切梁?jiǎn)卧?/p>

根據(jù)微分幾何理論[22],曲率是單位切向量對(duì)弧長(zhǎng)的一階導(dǎo)矢的模.本例中切向量是梁中心線上的位移梯度,通過(guò)位移場(chǎng)推導(dǎo)得到的曲率模型,本文定義為位移場(chǎng)曲率.圖1中梁變形后其截面不再與梁的中心軸線s垂直.根據(jù)梁位移場(chǎng),梁?jiǎn)卧猧上任意一點(diǎn)變形后的位移為:ri=ri(x,y).其對(duì)中心軸線s的一階導(dǎo)為

對(duì)中心線s的二階導(dǎo)為

E,G,A分別是梁?jiǎn)卧猧的彈性模量、剪切模量和橫截面積;剪切模量G=E/2(1+ν);As=ks·A,ks是剪切修正因子,對(duì)于橫截面是矩形的情況,ks=10(1+ν)/12+11ν;應(yīng)變?chǔ)沿梁中心線的分量,可由拉格朗日應(yīng)變張量確定[10]:

梁?jiǎn)卧猧的彈性力公式可由應(yīng)變勢(shì)能Ui對(duì)節(jié)點(diǎn)坐標(biāo)ei求偏導(dǎo)得到:

2.3 模型三:基于旋轉(zhuǎn)場(chǎng)曲率的二維幾何精確剪切梁?jiǎn)卧?/p>

曲率也可定義為曲線上某點(diǎn)的切線方向角對(duì)弧長(zhǎng)的轉(zhuǎn)動(dòng)率.通過(guò)角度場(chǎng)推導(dǎo)得到的曲率模型,本文定義為旋轉(zhuǎn)場(chǎng)曲率.由圖1所示,根據(jù)平面旋轉(zhuǎn)場(chǎng)理論,由O-XY逆時(shí)針轉(zhuǎn)動(dòng)到移動(dòng)矢量的轉(zhuǎn)動(dòng)矩陣為,其中α是與OX軸的夾角,即梁中線的切線方向角.沿著梁中心線切線方向,可體現(xiàn)梁中線彎曲效應(yīng).則單位切向量可寫成

由(26)式可知:cosα=?x1,由該式可得到旋轉(zhuǎn)場(chǎng)曲率公式:

該模型中的軸向勢(shì)能與剪切勢(shì)能與模型二中的表達(dá)式相同,彎曲勢(shì)能可利用旋轉(zhuǎn)場(chǎng)曲率得到:.彎曲勢(shì)能對(duì)節(jié)點(diǎn)坐標(biāo)ei的導(dǎo)數(shù)為

(30)式中的偏導(dǎo)項(xiàng)分別為:

同理,通過(guò)(26)式的第二式?x2=sinα也可推導(dǎo)出旋轉(zhuǎn)場(chǎng)曲率的表達(dá)式.

根據(jù)絕對(duì)節(jié)點(diǎn)坐標(biāo)方法,由虛功原理可推導(dǎo)梁?jiǎn)卧猧的動(dòng)力學(xué)方程[9]:分別采用模型一至三中的彈性力表達(dá)式,¨ei是節(jié)點(diǎn)加速度坐標(biāo).對(duì)于ANCF剪切梁?jiǎn)卧猧,其單元質(zhì)量陣公式為

其中,ρi,Vi是梁?jiǎn)卧猧的密度和體積.由(35)式可知,通過(guò)絕對(duì)節(jié)點(diǎn)坐標(biāo)方法建立的單元質(zhì)量陣為常數(shù)陣,不包含科氏力與離心力,不存在變形和轉(zhuǎn)動(dòng)耦合附加的慣性力.通過(guò)虛功原理可推導(dǎo)單元的外力陣公式:

Fi為外力矢量,外力陣同樣為常數(shù)矩陣.

3 算例分析

3.1 靜力學(xué)分析

本節(jié)引入典型的懸臂梁受集中載荷問(wèn)題,考察所提出基于旋轉(zhuǎn)場(chǎng)曲率的ANCF剪切梁模型在靜力學(xué)分析中的準(zhǔn)確性.圖1所示為懸臂梁受集中載荷模型,其物理參數(shù)與文獻(xiàn)[14]相同.長(zhǎng)為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.

考察小變形算例,設(shè)置末端集中力為F=?5×105×h3N.在小變形假設(shè)下,梁中性層在變形后長(zhǎng)度不變,端點(diǎn)沿X方向的變形可以忽略不計(jì).Y向撓度的理論解[23]為,其中Iz為梁截面慣性矩.通過(guò)表1的數(shù)據(jù)對(duì)比發(fā)現(xiàn),由于模型一存在嚴(yán)重的剪切閉鎖問(wèn)題,其Y向位移與理論解誤差達(dá)到22%.而模型二和模型三采用幾何精確梁理論構(gòu)建系統(tǒng)彈性力,其中軸向力、剪切力和彎曲力相互解耦,可緩解剪切閉鎖,因此由它們計(jì)算的Y向位移均可收斂于理論解,并且收斂速度快.減小梁的高度至h=0.1 m,末端集中力為F=?1×106×h3N,其他參數(shù)不變.表2結(jié)果表明模型一的誤差為25%,而模型二和三可收斂到理論解.與表1相比,表2中變形量增大,模型二和三需要收斂到理論解的個(gè)數(shù)增多.

表2 小變形情況下三種彈性力模型計(jì)算的端部位移對(duì)比2Table 2.2ndcomparison of tip point displacement using three elastic force models in case of small deformation.

表3 非小變形情況下三種彈性力模型計(jì)算的端部位移對(duì)比Table 3.Comparison of tip point displacement using three elastic force models in case of non-small deformation.

考察非小變形情況,設(shè)置末端集中力為F=?5×108×h3N,其他參數(shù)與小變形算例1相同.采用該集中力可導(dǎo)致相對(duì)大的變形量.采用文獻(xiàn)[14]表2中用商業(yè)軟件計(jì)算的解作為非線性大變形下懸臂梁X,Y向撓度的參考解.在大變形情況下,剪切效應(yīng)加劇,另外由于ANCF單元在Y向上采用線性插值會(huì)導(dǎo)致嚴(yán)重的剪切閉鎖問(wèn)題,所以在大變形算例中三種模型的收斂精度均不如小變形算例,見(jiàn)表3.數(shù)據(jù)對(duì)比顯示模型三的收斂速度和精度最高.以上算例證明了所提出的模型在靜力學(xué)分析中的正確性.

3.2 特征頻率分析

本節(jié)對(duì)所提出的ANCF剪切梁?jiǎn)卧奶卣黝l率進(jìn)行分析,選擇具有理論解的簡(jiǎn)支梁作為研究對(duì)象.簡(jiǎn)支梁物理參數(shù)與文獻(xiàn)[14]中相同.其中長(zhǎng)為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 簡(jiǎn)支梁模型Fig.3.Simply supported beam model.

圖4 基于旋轉(zhuǎn)場(chǎng)曲率的簡(jiǎn)支梁模型的振型和固有頻率Fig.4.Mode shapes and eigenfrequencies(rad/s)of simply supported beam based on rotation-based curvature.

表43 種彈性力模型下簡(jiǎn)支梁的彎曲和軸向固有頻率(rad/s)對(duì)比Table 4.Bending and longitudinal eigenfrequencies(rad/s)of simply supported beam element for di ff erent elastic formulation.

圖4為本文模型三在簡(jiǎn)支邊界條件下的固有頻率和振型.因?yàn)橥ㄟ^(guò)邊界條件消除了3個(gè)自由度,所以共計(jì)算出有9種振型和對(duì)應(yīng)的固有頻率.ANCF單元的軸向采用三次插值技術(shù)產(chǎn)生三種軸向振型(mode 2,mode 4和mode 5)和兩種彎曲振型(mode 1和mode 3).另外,還存在兩種剪切振型(mode 6和mode 7)以及ANCF單元特有的厚度振型(mode 8和mode 9).

表53 種彈性力模型下簡(jiǎn)支梁的剪切和厚度固有頻率(rad/s)對(duì)比Table 5.Shearing and thickness eigenfrequencies(rad/s)of simply supported beam element for di ff erent elastic formulation.

表4為簡(jiǎn)支梁?jiǎn)卧謩e采用3種彈性力模型所計(jì)算出的彎曲和軸向固有頻率的對(duì)比,表5為其剪切和厚度固有頻率對(duì)比.由結(jié)果可知,與其他兩個(gè)模型相比,本文提出的模型三收斂性更好,收斂精度可至6位有效數(shù)字.引起二階彎曲振動(dòng)頻率誤差較高是因?yàn)楸疚膯卧赮向上采用線性插值,會(huì)導(dǎo)致單元的剪切閉鎖問(wèn)題,使得變形過(guò)程中剪切應(yīng)變能占據(jù)主導(dǎo),過(guò)多的剪切力導(dǎo)致單元不能發(fā)生正常彎曲變形,造成彎曲剛度增大,單元變得更“剛”,從而引起二階彎曲振動(dòng)頻率增大.

3.3 動(dòng)力學(xué)分析

圖5所示為柔性單擺,O點(diǎn)為旋轉(zhuǎn)鉸.該算例是絕對(duì)節(jié)點(diǎn)坐標(biāo)方法的一個(gè)經(jīng)典算例.單擺的物理參數(shù)如下:長(zhǎng)為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方向位移.考察在重力作用下單擺的動(dòng)態(tài)響應(yīng).仿真中單擺被劃分為12個(gè)單元,仿真時(shí)間為1.1 s.

本算例中柔性單擺在其重力作用下進(jìn)行大范圍轉(zhuǎn)動(dòng),運(yùn)動(dòng)過(guò)程中發(fā)生大變形現(xiàn)象.為對(duì)比文中所述的三種剪切梁?jiǎn)卧Ｐ吞幚泶笞冃螁?wèn)題的正確性,采用ANSYS中兩節(jié)點(diǎn)平面梁BEAM3單元作為對(duì)比.由圖6—圖8可見(jiàn),在彈性模量為E=7 GPa時(shí),三種模型計(jì)算的末端橫向位移,中點(diǎn)橫向位移和末端端點(diǎn)軌跡均與ANSYS結(jié)果符合,證明它們均能處理大變形問(wèn)題,并驗(yàn)證了本文所提模型三的正確性.

圖5 柔性單擺模型Fig.5.Flexible pendulum.

為了更好地說(shuō)明模型三處理柔性大變形的能力,將彈性模量降低至E=2 GPa.在此彈性模量下,梁將變得更柔,大變形現(xiàn)象更明顯.圖9—圖11顯示在E=2 GPa時(shí),三種模型計(jì)算的結(jié)果均與ANSYS結(jié)果符合.繼續(xù)將彈性模量降至E=0.7 GPa,由圖12可見(jiàn),在此低模量下三種模型計(jì)算的結(jié)果也可與ANSYS結(jié)果相近.以上算例均證明了模型三可用于處理柔性大變形問(wèn)題.

圖6 柔性單擺末端橫向位移(E=7 GPa)Fig.6.The tip vertical displacement of fl exible pendulum(E=7 GPa).

圖7 柔性單擺中點(diǎn)橫向位移(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 柔性單擺中點(diǎn)橫向位移(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是模型三在不同彈性模量下的能量變化曲線.由圖可知,該模型總能量為零,符合能量守恒定律.從能量的角度驗(yàn)證了模型三的正確性.

圖13 系統(tǒng)能量分布(旋轉(zhuǎn)場(chǎng)曲率模型,E=7 GPa)Fig.13.The energy distribution of fl exible pendulum(RB curvature beam model,E=7 GPa).

圖14 系統(tǒng)能量分布(旋轉(zhuǎn)場(chǎng)曲率模型,E=2 GPa)Fig.14.The energy distribution of fl exible pendulum(RB curvature beam model,E=2 GPa).

圖15系統(tǒng)能量分布(旋轉(zhuǎn)場(chǎng)曲率模型,E=0.7 GPa)Fig.15.The energy distribution of fl exible pendulum(RB curvature beam model,E=0.7 GPa).

圖16 —圖18顯示了模型一至模型三在不同時(shí)刻的構(gòu)型圖,分別對(duì)三種模型取間隔0.1 s的空間構(gòu)型進(jìn)行對(duì)比.圖中綠色實(shí)線代表模型一,藍(lán)色點(diǎn)畫線代表模型二,紅色虛線代表模型三.由圖16可見(jiàn),在E=7 GPa時(shí)三種模型的運(yùn)動(dòng)軌跡一致,說(shuō)明在此彈性模量下由幾何精確梁理論推導(dǎo)的模型(模型二和模型三)和由連續(xù)介質(zhì)力學(xué)理論推導(dǎo)的模型(模型一)在模擬單擺大變形時(shí)的效果基本相同.隨著彈性模量的增大,模型二與模型三在運(yùn)動(dòng)過(guò)程中相較模型一顯得“更柔”,主要因?yàn)槟Ｐ鸵皇邱詈夏Ｐ?而模型二與模型三是解耦模型,各個(gè)變形模態(tài)之間互不影響.所以圖17和圖18中模型二和模型三的運(yùn)動(dòng)軌跡和空間構(gòu)型更為接近,驗(yàn)證了兩種不同物理場(chǎng)得到的曲率模型的有效性.

圖16 (網(wǎng)刊彩色)柔性單擺各時(shí)刻構(gòu)型圖(E=7 GPa)(綠色實(shí)線,模型一;藍(lán)色點(diǎn)畫線,模型二;紅色虛線,模型三)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 (網(wǎng)刊彩色)柔性單擺各時(shí)刻構(gòu)型圖(E=2 GPa)(網(wǎng)刊彩色)(綠色實(shí)線:模型一;藍(lán)色點(diǎn)畫線:模型二;紅色虛線:模型三)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 (網(wǎng)刊彩色)柔性單擺各時(shí)刻構(gòu)型圖(E=0.7 GPa)(綠色實(shí)線,模型一;藍(lán)色點(diǎn)畫線,模型二;紅色虛線,模型三)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).

柔性梁在大轉(zhuǎn)動(dòng)過(guò)程中剪切效應(yīng)逐漸占據(jù)主導(dǎo),為了進(jìn)一步研究三種不同的彈性力計(jì)算的剪切應(yīng)變?cè)诖筠D(zhuǎn)動(dòng)大變形中的變化,彈性模量選取E=0.7 GPa,其余梁的物理參數(shù)不變,柔性單擺被劃分為40個(gè)單元.圖19是單擺最后一個(gè)單元內(nèi)中點(diǎn)的剪切應(yīng)變變化曲線.從圖中可知,模型二與模型三得到的剪切應(yīng)變相近,而它們均比模型一的剪切應(yīng)變大,其最大差異可達(dá)到0.01.其差異可由圖20和圖21解釋,分別是單擺最后一個(gè)單元中點(diǎn)的rx,ry方向與大小變化曲線.由圖可知,模型二與模型三在rx,ry方向變化上與模型一基本相似,而rx,ry的幅值變化與模型一有較大差異,這也是導(dǎo)致圖19剪切應(yīng)變差異的主要原因.絕對(duì)節(jié)點(diǎn)坐標(biāo)法中的應(yīng)變是位移梯度的函數(shù),例如正應(yīng)力表達(dá)式

它可限制ry劇烈變化,確保單元的不可壓縮性.模型一由于是耦合模型,rx與ry變化相互影響,所以其ry幅值的變化范圍比模型二和模型三要更明顯.

圖19 單擺最后一個(gè)單元內(nèi)中點(diǎn)的剪切應(yīng)變變化Fig.19.Change of shear strain at middle point of last element.

由于ANCF剪切梁?jiǎn)卧趚,y方向上插值階次的不同將會(huì)導(dǎo)致嚴(yán)重的剪切閉鎖現(xiàn)象,采用本文模型一至模型三分別研究各模型緩解剪切閉鎖的能力.表6為分別采用三種模型時(shí)的單擺中點(diǎn)剪切應(yīng)變收斂情況,可以看出采用解耦模型(模型二和模型三)得到的剪應(yīng)變收斂個(gè)數(shù)比耦合模型(模型一)要少,模型三收斂最快,證明模型三在緩解剪切閉鎖問(wèn)題上具有較強(qiáng)的優(yōu)勢(shì).模型一由于各個(gè)變形模態(tài)相互耦合,剪切應(yīng)變收斂最慢.

圖20 單擺最后一個(gè)單元內(nèi)中點(diǎn)的rx和ry方向變化(a)rx;(b)ryFig.20.Orientation of the rxand ryvector at middle point of last element:(a)rx;(b)ry.

圖21 單擺最后一個(gè)單元內(nèi)中點(diǎn)的rx和ry幅值變化(a)rx;(b)ryFig.21.Magnitude of the rxand ryvector at middle point of last element:(a)rx;(b)ry.

表6 三種彈性力模型的單擺中點(diǎn)剪切應(yīng)變收斂對(duì)比Table 6.Comparison of shear strain at middle point using three di ff erent elastic force models

4 結(jié)論

本文采用梁中心線切線方向角推導(dǎo)了剪切梁的精確曲率公式,提出了基于旋轉(zhuǎn)場(chǎng)曲率的二維ANCF剪切梁?jiǎn)卧Ｐ?通過(guò)與傳統(tǒng)二維ANCF剪切梁?jiǎn)卧Ｐ?、基于位移?chǎng)曲率的二維ANCF單元模型對(duì)比分析,可得到以下結(jié)論:

1)在靜力學(xué)分析以及特征頻率分析中,本文提出的模型三與理論解收斂度最好,效率最高,且收斂精度和速度優(yōu)于模型一;

2)在動(dòng)力學(xué)分析中,模型三顯示了良好的計(jì)算精度,可有效預(yù)測(cè)柔性構(gòu)件運(yùn)動(dòng)過(guò)程中的大變形現(xiàn)象,此外,對(duì)比驗(yàn)證了兩種不同物理場(chǎng)推導(dǎo)的曲率模型在單元幾何非線性分析中的有效性;

3)動(dòng)力學(xué)分析還驗(yàn)證了模型三可有效緩解剪切閉鎖問(wèn)題,且在計(jì)算柔性梁剪切應(yīng)變中收斂速度最快.

感謝美國(guó)伊利諾伊大學(xué)芝加哥分校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物理學(xué)報(bào)60 044501]

[5]He X S,Deng F Y,Wang R 2010 Acta Phys.Sin.59 1428(in Chinese)[和興鎖,鄧峰巖,王睿2010物理學(xué)報(bào)59 1428]

[6]Chen S J,Zhang D G,Hong J Z 2013 Chin.J.Theor.Appl.Mech.45 251(in Chinese)[陳思佳,章定國(guó),洪嘉振2013力學(xué)學(xué)報(bào)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)[田強(qiáng),張?jiān)魄?陳立平,覃剛2010力學(xué)進(jìn)展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)[章孝順,章定國(guó),陳思佳,洪嘉振2016物理學(xué)報(bào)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

?國(guó)家自然科學(xué)基金重大項(xiàng)目(批準(zhǔn)號(hào):61690210,61690211)和國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):61603304,11472213)資助的課題.

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

?2017中國(guó)物理學(xué)會(huì)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