基于旋轉(zhuǎn)場曲率的二維剪切梁單元建模?

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

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

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

基于旋轉(zhuǎn)場曲率的二維剪切梁單元建模?

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

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

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

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

對二維剪切梁單元進行研究,利用平面旋轉(zhuǎn)場理論推導(dǎo)了精確曲率模型.采用幾何精確梁理論構(gòu)建了剪切梁單元彈性力矩陣.通過絕對節(jié)點坐標方法建立了系統(tǒng)的非線性動力學(xué)方程,提出基于旋轉(zhuǎn)場曲率的二維剪切梁單元,并分別引入經(jīng)典二維剪切梁單元和基于位移場曲率的二維剪切梁單元進行比較研究.首先,靜力學(xué)分析證明了所提模型的正確性;其次,特征頻率分析驗證了模型可與理論解符合,收斂精度高,并且能準確地預(yù)測單元固有頻率對應(yīng)的振型;最后,在非線性動力學(xué)問題上,通過與ANSYS結(jié)果對比分析,證明了該模型可有效處理柔性大變形問題,并且與經(jīng)典二維剪切梁單元相比具有緩解剪切閉鎖的優(yōu)勢.因此,本文提出的基于旋轉(zhuǎn)場曲率的二維剪切梁單元在處理幾何非線性問題中具有較大的應(yīng)用潛力.

絕對節(jié)點坐標法,旋轉(zhuǎn)場曲率,幾何精確梁理論,剪切閉鎖

1 引言

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

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

本文進一步針對二維ANCF剪切梁的曲率定義進行研究.基于平面旋轉(zhuǎn)場理論推導(dǎo)得到精確曲率公式,表征為梁中心線切線方向轉(zhuǎn)角的函數(shù);采用幾何精確梁理論構(gòu)建系統(tǒng)彈性力矩陣,提出基于旋轉(zhuǎn)場曲率的二維ANCF剪切梁單元模型;同時引入基于連續(xù)介質(zhì)力學(xué)理論的經(jīng)典ANCF剪切梁單元和基于微分幾何理論的位移場曲率ANCF剪切梁單元作為對比研究.分別通過靜力學(xué)分析、特征頻率分析以及動力學(xué)分析驗證了本文所提的基于旋轉(zhuǎn)場曲率的二維ANCF剪切梁單元模型在計算柔性機構(gòu)大變形、大轉(zhuǎn)動問題中的準確性和適用性.

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

2 剪切梁模型

2.1 模型一:基于連續(xù)介質(zhì)力學(xué)理論的二維ANCF剪切梁單元

二維含剪切變形的ANCF梁單元最早由Omar和Shabana提出[9],如圖1所示,其中O-XY為慣性坐標系,l為梁的長度.

對于二節(jié)點的ANCF剪切梁單元i,其梁上任意一點在慣性系中的絕對位移矢量為

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

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

梁單元i的型函數(shù)Si(x,y)為:

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

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

對于梁單元的初始狀態(tài)是水平擺放情況,J0為單位陣.Green-Lagrangian應(yīng)變張量εi的表達式為

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

依據(jù)廣義胡克定律,ANCF剪切梁單元i的應(yīng)力張量σi為

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

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

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

則其彈性力矢量為

其中,

經(jīng)典剪切梁單元模型中軸向伸長變形、剪切變形和彎曲變形相互耦合在一起,因此曲率公式不是顯式形式.

2.2 模型二:基于位移場曲率的二維幾何精確剪切梁單元

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

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

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

梁單元i的彈性力公式可由應(yīng)變勢能Ui對節(jié)點坐標ei求偏導(dǎo)得到:

2.3 模型三:基于旋轉(zhuǎn)場曲率的二維幾何精確剪切梁單元

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

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

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

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

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

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

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

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

3 算例分析

3.1 靜力學(xué)分析

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

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

表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.

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

3.2 特征頻率分析

本節(jié)對所提出的ANCF剪切梁單元的特征頻率進行分析,選擇具有理論解的簡支梁作為研究對象.簡支梁物理參數(shù)與文獻[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 基于旋轉(zhuǎn)場曲率的簡支梁模型的振型和固有頻率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種振型和對應(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 種彈性力模型下簡支梁的剪切和厚度固有頻率(rad/s)對比Table 5.Shearing and thickness eigenfrequencies(rad/s)of simply supported beam element for di ff erent elastic formulation.

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

3.3 動力學(xué)分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

感謝美國伊利諾伊大學(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é)報60 044501]

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

[6]Chen S J,Zhang D G,Hong J Z 2013 Chin.J.Theor.Appl.Mech.45 251(in Chinese)[陳思佳,章定國,洪嘉振2013力學(xué)學(xué)報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力學(xué)進展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物理學(xué)報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

?國家自然科學(xué)基金重大項目(批準號:61690210,61690211)和國家自然科學(xué)基金(批準號:61603304,11472213)資助的課題.

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

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