陳占魁 羅 凱 田 強

(北京理工大學宇航學院力學系,北京 100081)

引言

20 世紀中期,美國建筑師Fuller[1]從藝術(shù)家Snelson[2]的雕塑作品中獲得靈感,將“tension”和“integrity”融合為一個新的單詞,提出了“tensegrity”(張拉整體)的概念.早期的張拉整體結(jié)構(gòu)大多是從藝術(shù)角度設(shè)計的作品,未考慮其工程應(yīng)用價值,直到20 世紀80 年代,在Pellegrino[3-4]和Connelly 等[5]的推動下,張拉整體結(jié)構(gòu)才被越來越多的工程師所熟知.

張拉整體結(jié)構(gòu)具有整體均勻承載、輕質(zhì)、高收納比等優(yōu)異結(jié)構(gòu)特性,因此逐漸得到航天、生物力學、機器人等領(lǐng)域研究人員的廣泛關(guān)注.在航天工程領(lǐng)域,Tibert 和Pellegrino[6]提出了一種可展開張拉整體棱柱結(jié)構(gòu),通過對其幾何和內(nèi)力分析及其實驗驗證,證明了該結(jié)構(gòu)具備較好的展開可靠度和較高的收納比;SunSpiral 等[7]提出了張拉整體行星探測器的概念,模擬和分析了張拉整體結(jié)構(gòu)展開、降落和表面移動過程的動力學行為,驗證了其用于Titan 行星探測任務(wù)的可行性;Caluwaerts 等[8]對用于行星探測的球形張拉整體機器人進行動力學仿真,研究不同控制方法下系統(tǒng)動力學響應(yīng),提出了生物啟發(fā)的張拉整體結(jié)構(gòu)控制策略;Iscen 等[9]通過分析不同層級和頻率的測量信號,獲得了保持其以適當速度自主滾動的分布式驅(qū)動信號.在生物力學領(lǐng)域,Ingber[10-11]利用張拉整體結(jié)構(gòu)模擬細胞及其骨架整體結(jié)構(gòu)的力學行為,解釋細胞結(jié)構(gòu)的形變、運動和對機械信號的感知與傳導(dǎo)機理,進一步理解生物體的分級組織機理.在機器人領(lǐng)域,Rovira 和Tur[12]推導(dǎo)了由剛性桿和彈性繩組成的棱柱張拉整體機器人動力學普遍方程,并通過仿真研究了在指定運動下分布變長度繩的多狀態(tài)離散控制方法;陳竑希[13]對四桿張拉整體機器人進行驅(qū)動方案設(shè)計和動力學仿真分析,并通過實驗驗證了理論結(jié)果的正確性;Tietz 等[14]首次采用CPG(central pattern generators) 控制實現(xiàn)多節(jié)機器人爬行運動的波形調(diào)控,使其可在各種復(fù)雜地形上運動.

雖然張拉整體結(jié)構(gòu)涉及的領(lǐng)域不同,但研究的問題大多圍繞找形與拓撲、結(jié)構(gòu)分析和形態(tài)控制這3 方面展開.找形即在給定張拉整體結(jié)構(gòu)拓撲(即節(jié)點連接方式)情況下求解其自平衡節(jié)點位置坐標,Barnes[15]基于動力松弛法提出了張拉整體結(jié)構(gòu)找形的數(shù)值計算流程,構(gòu)造了較穩(wěn)定和較快收斂的人工運動阻尼過程,該方法的優(yōu)勢是可實現(xiàn)大擾動找形和考慮了幾何非線性、材料非線性等.找形的另一類方法是力密度法,而力密度法原本僅適用于純張力結(jié)構(gòu)的找形,當用于同時包含張力和壓力部件的結(jié)構(gòu)找形時,其平衡方程可能發(fā)生奇異,對此Vassart和Motro[16]給出了詳細的理論分析,并提出了“冗余”節(jié)點的概念,基于多輸入?yún)?shù)實現(xiàn)張拉整體結(jié)構(gòu)找形;Koohestani[17]和Chen 等[18]分別采用基于仿生學的遺傳算法和蟻群算法,將張拉整體結(jié)構(gòu)找形轉(zhuǎn)化為優(yōu)化問題求解,提升了基于力密度找形方法的適用性.相比找形研究,張拉整體結(jié)構(gòu)的拓撲研究相對較少,Lee 等[19]以張拉整體結(jié)構(gòu)節(jié)點坐標為輸入,在節(jié)點全部連接中給出滿足桿不連續(xù)條件的桿編號集合,利用極小張拉整體構(gòu)型理論消去非必要繩連接,并利用力密度法和遺傳算法給出找力結(jié)果;Wang等[20]提出了含剛體張拉整體結(jié)構(gòu)拓撲設(shè)計的普適方法,利用MILP(混合整數(shù)規(guī)劃優(yōu)化)算法求解了拓撲設(shè)計的優(yōu)化問題;Liu 和Paulino[21]基于基結(jié)構(gòu)法和MILP 算法實現(xiàn)張拉整體結(jié)構(gòu)的拓撲優(yōu)化,通過最大化自平衡力之和進行構(gòu)型提取,得到了穩(wěn)定和對稱的張拉整體結(jié)構(gòu).在結(jié)構(gòu)分析方面,Goyal 和Skelton[22]考慮張拉整體結(jié)構(gòu)繩的質(zhì)量而忽略桿的彈性,建立了以剛性桿和分段彈簧集中質(zhì)量的繩模型組成多體動力學模型,通過動力學仿真研究典型張拉整體結(jié)構(gòu)動力學行為;Goyal 等[23]提出了考慮局部和全局失效的輕質(zhì)張拉整體結(jié)構(gòu)設(shè)計方法,在幾何設(shè)計的基礎(chǔ)上進一步考慮了結(jié)構(gòu)的輕量化和承載能力;朱世新等[24]建立了組合式數(shù)字狀張拉整體結(jié)構(gòu),對其大變形非線性力學行為進行數(shù)值模擬以實現(xiàn)力學性能分析,并通過實驗測試研究了數(shù)字8 狀結(jié)構(gòu)的靜力學和動力學特性.在結(jié)構(gòu)的形態(tài)控制方面,沈黎元等[25]基于索桿體系矩陣分析理論,給出了預(yù)應(yīng)力索結(jié)構(gòu)位移控制的線性求解方法,可望用于精密結(jié)構(gòu)的形狀控制;Shea 等[26]結(jié)合動力松弛法和模擬退火搜索算法提出了自適應(yīng)張拉整體結(jié)構(gòu)的全局形態(tài)控制框架;Begey 等[27]評估了驅(qū)動類型和位置對Snelson cross 張拉整體機構(gòu)的影響,給出不同驅(qū)動模式的性能圖,以拓展張拉整體機器人的設(shè)計可能性.

以往對張拉整體結(jié)構(gòu)或系統(tǒng)力學建模中,一般假設(shè)繩索只發(fā)生純拉伸變形,且桿只發(fā)生純壓縮變形或假設(shè)為不變形的剛體.在上述假設(shè)下,當張拉整體系統(tǒng)發(fā)生大變形時,部件會發(fā)生大轉(zhuǎn)動而引起幾何非線性問題,使結(jié)構(gòu)線性分析失效.Kebiche 等[28]基于非線性有限元思想,利用虛功原理的增量分析,在完全拉格朗日框架下推導(dǎo)了桿的非線性剛度矩陣,實現(xiàn)了張拉整體系統(tǒng)的幾何非線性分析;在此基礎(chǔ)上,Kahla 和Kebiche[29]引入了材料非線性本構(gòu),將上述方法拓展到張拉整體系統(tǒng)的非線性彈塑性分析;此外,Zhang 等[30]針對大位移小應(yīng)變系統(tǒng)的高效分析,采用共旋坐標法實現(xiàn)了張拉整體系統(tǒng)的幾何和材料非線性分析.進一步,對于由細長桿或軟材料桿組成的柔軟張拉整體系統(tǒng),在系統(tǒng)變形或運動中桿還可能發(fā)生彎曲變形,需研究桿后屈曲局部變形對系統(tǒng)整體承載性能與動力學行為的影響.Rimoli[31]提出了張拉整體結(jié)構(gòu)柔性桿的一種等效建模方法,將桿的連續(xù)模型等效為集中參數(shù)模型,采用較少的自由度捕捉了桿的屈曲、后屈曲行為及動力學特性,可用于柔性張拉整體著陸器[31]、復(fù)雜張拉整體點陣材料[32]的動力學和非線性響應(yīng)高效分析;然而,由于文獻[31]中計算兩端簡支桿彎曲振動固有頻率的解析表達式有誤(即原文式(26) 中固有頻率應(yīng)與階數(shù)平方成正比,而不是與階數(shù)成正比),致使在合適取值下計算的離散模型前兩階固有頻率相對誤差均為4.5%,而事實上第二階固有頻率相對誤差應(yīng)為15%,即彎曲固有振動的計算精度不高.

針對張拉整體結(jié)構(gòu)中柔性細長桿動力學等效建模和高效計算,本文提出了考慮拉壓和彎曲變形的五節(jié)點離散桿模型,從能量角度推導(dǎo)了連續(xù)模型與離散模型的動力學等效過程,改進了Rimoli 降階模型[31],使得連續(xù)模型和離散模型的桿彎曲振動前兩階固有頻率相對誤差均降低至1%以內(nèi),提高了動力學計算精度.進一步,以六桿球形張拉整體結(jié)構(gòu)為研究對象,建立和求解多柔體系統(tǒng)動力學方程,開展壓縮、模態(tài)和沖擊的靜/動態(tài)仿真與實驗測試,對比驗證本文動力學建模和計算方法的有效性.

1 屈曲細長桿的動力學等效

1.1 連續(xù)模型與離散模型

本節(jié)給出了細長桿連續(xù)模型的降階離散模型,該離散模型采用五節(jié)點彈/扭簧集中質(zhì)量模型來等效連續(xù)桿的動力學特性,如圖1 所示.其中,均質(zhì)連續(xù)桿的長度為L,截面面積為A,截面慣性矩為I,材料密度為ρ,楊氏模量為E.離散模型5 個節(jié)點的位置和軸向剛度假設(shè)為均勻分布,軸向彈簧自由長度均為L/4 且剛度均為K1;而其分布集中質(zhì)量大小和分布彎曲剛度假設(shè)為非均勻,考慮受壓桿件變形對稱性,設(shè)集中質(zhì)量和彎曲剛度在質(zhì)心兩側(cè)對稱分布,分布集中質(zhì)量大小表示為m1,m2和m3,分布彎曲剛度大小表示為Kt1和Kt2.1.2～1.4 節(jié)將從能量角度給出上述集中參數(shù)動力學等效過程.

圖1 細長桿連續(xù)模型和離散模型示意圖Fig.1 Schematic diagram of the continuous and discrete model of a slender bar

1.2 剛度大小等效

根據(jù)材料力學中桿的拉壓變形理論,容易得到等效軸向剛度為

隨著軸向壓力增大,桿會發(fā)生屈曲,為了描述桿彎曲變形受力狀態(tài),引入彎曲剛度Kt1和Kt2.如圖2所示,離散模型在軸向載荷P的作用下呈現(xiàn)壓縮和彎曲耦合變形狀態(tài),且關(guān)于中心點對稱,此時系統(tǒng)的總勢能表達為

其中d1,d2,θ1和θ2等為局部變量,如圖2 所示,L1=0.25L,ΔL是載荷P引起的位移,可表示為

圖2 離散模型的屈曲構(gòu)型示意圖Fig.2 Schematics of buckling configuration for the discretization scheme

以d1,d2,θ1和θ2為廣義坐標,由最小勢能原理可得以下平衡方程

為了捕捉屈曲臨界載荷,考慮對上述后屈曲段非線性平衡方程進行轉(zhuǎn)角趨于零的線性化.當轉(zhuǎn)角θ1和θ2趨于零時有

將式(5)代入式(4)可得

由式(6)的前兩個方程可得

式(6)的后兩個方程可改寫為

顯然,式(8) 的一組平凡解為θ1=θ2=0,該零解與式(7)描述了桿屈曲前純壓縮狀態(tài);由式(8)的非零解條件可得到屈曲臨界載荷Pcr應(yīng)滿足的關(guān)系式,即

設(shè)彎曲剛度分布參數(shù)為n,使得Kt2=nKt1,將該式和式(7)代入式(9)得

兩端簡支桿的Euler 屈曲臨界載荷為

將式(11)代入式(10),求解關(guān)于Kt1的一元二次方程,得其二根為

其中,剛度分布參數(shù)n將通過1.4 節(jié)彎曲振動固有特性等效確定,± 分別對應(yīng)Kt1的較大根和較小根.對于同一根桿的等效離散模型,當Kt1取較大根時,后文確定的剛度分布參數(shù)n＜1,即Kt2＜Kt1,則中心點承受的彎矩小于其兩側(cè)點承受的彎矩;反之,n＞1,則中心點承受的彎矩大于其兩側(cè)點承受的彎矩.由于后者方式等效的彎曲振動固有頻率誤差較大,本文僅選取前者方式等效,即上式取正號.此外,當模擬細長桿時,π2I/(AL2)?1,可忽略上式中的該項.

進一步,為了驗證該離散模型是否可較準確地預(yù)測細長桿的后屈曲變形,將其計算結(jié)果與商用軟件ABAQUS 的有限元計算結(jié)果進行對比.設(shè)截面半徑為5 mm、長度為0.2 m 和楊氏模量為19 MPa 的桿,兩端為簡支約束,一端位置固定,另一端可沿兩端點連線方向滑動并在該端施加縱向壓力,分別采用離散模型和有限元軟件計算其受壓屈曲全過程.圖3所示為壓力與縱向位移對比圖,可見兩者結(jié)果吻合較好.

圖3 受壓桿有限元模型與離散模型載荷-位移對比圖Fig.3 Load-displacement comparison for the continuous and discrete bar under compression

1.3 質(zhì)量大小等效

通過對比連續(xù)桿和離散桿的動能可求出節(jié)點的質(zhì)量.假設(shè)連續(xù)桿上存在線速度場u和v

其中,u代表慣性坐標系下沿桿軸向的速度分量,v代表慣性坐標系下沿桿橫向的速度分量,x∈[-L/2,L/2]表示參考構(gòu)型下沿桿軸向的局部坐標,c1,c2,c3和c4表示任意常數(shù).該線速度場雖不能描述平面柔性桿的全部運動,但可表示部分主要發(fā)生的運動,即剛體平移、剛體轉(zhuǎn)動及軸向均勻變形.在該線速度場描述下,連續(xù)桿的動能可表達為

另一方面,通過線速度場u和v,離散模型的動能可表達為

由于c1,c2,c3和c4表示任意常數(shù),為了使連續(xù)桿和離散模型的動能相等,和前面的系數(shù)須保持一致,因此由式(14)和式(15)可得

式(16) 只有兩個方程,無法直接求得m1,m2和m3的解,故設(shè)質(zhì)量分布參數(shù)為c,使得m3=cm2,則由式(16)解得

至此,已全部求出離散模型的軸向剛度、彎曲剛度和節(jié)點質(zhì)量.但是其中的彎曲剛度和節(jié)點質(zhì)量并不是完全不變的,通過調(diào)節(jié)參數(shù)n和c可對彎曲剛度和節(jié)點質(zhì)量進行分配.1.4 節(jié)將對n和c的取值進行討論,使離散模型可以更準確地描述連續(xù)桿的動力學行為.

1.4 彎曲剛度和集中質(zhì)量分布特性的確定

當計算離散模型在無應(yīng)力狀態(tài)時的固有頻率時,軸向振動模態(tài)和彎曲振動模態(tài)可以獨立研究.桿橫向運動構(gòu)型描述如圖4 所示,V1,V2,V3,V4和V5表示各質(zhì)點橫向位移(即與水平參考線的距離),θ1,θ2,θ3和θ4表示各離散直線段的轉(zhuǎn)角(即與水平參考線的夾角),則彎曲應(yīng)變能可表達為

圖4 模態(tài)分析的構(gòu)型示意圖Fig.4 Schematics of configuration for modal analysis

當桿橫向運動在水平參考線附近進行線性化時,各轉(zhuǎn)角趨于零,有

將式(19)代入式(18),將彎曲應(yīng)變能改寫為

上述彎曲應(yīng)變能對橫向位移求兩次偏導(dǎo),可得離散系統(tǒng)的對稱剛度矩陣表達為

其中Kij=?2W/(?Vi?Vj)且i,j=1,2,3,4,5.文末附錄給出了當考慮細長桿時上述剛度矩陣各分量的表達式.

根據(jù)式(17) 所求的質(zhì)點質(zhì)量,可得出離散系統(tǒng)的質(zhì)量矩陣為

計算固有頻率時,邊界條件可通過消去對應(yīng)自由度的方式施加.假設(shè)桿兩端為簡直邊界條件,則分別消去以上剛度矩陣和質(zhì)量矩陣的首行首列和尾行尾列,利用剩余部分和求解廣義特征值方程,即

解關(guān)于ω2的一元三次方程,得離散模型前三階固有頻率平方為

其中

另一方面,簡支桿前三階固有頻率平方的解析解為

通過計算離散模型與連續(xù)模型的固有頻率相對誤差來確定參數(shù)n和c的取值,定義如下頻率相對誤差

本文構(gòu)造的離散模型,考慮對連續(xù)模型前兩階固有特性進行捕捉,因此定義前兩階固有頻率均方根誤差如下

由此將Re(n,c)表示為分布參數(shù)n和c的函數(shù),其函數(shù)圖像如圖5 所示.

圖5 離散模型隨n 和c 變化的前兩階固有頻率均方根誤差Fig.5 RMS(root-mean-square)error in the first two natural frequencies of the discrete bar as the function of n and c

由圖5 函數(shù)圖像可得當n=0.39 且c=0.60 時,前兩階固有頻率的均方根誤差達到最小值.將求得的n和c代入式(27)可得到前三階固有頻率誤差.如式(29)所示,前兩階固有頻率誤差較小,但第三階固有頻率誤差較大,原因是考慮到張拉整體結(jié)構(gòu)柔軟特性和實際工況受載情況,式(28)在求解調(diào)節(jié)參數(shù)n和c時僅考慮了對系統(tǒng)動力學行為產(chǎn)生主要影響的前兩階固有模態(tài),而截斷了影響較小的更高階模態(tài).

至此,細長桿的降階離散模型分布剛度和分布質(zhì)量動力學等效全部完成,考慮了桿的軸向拉壓剛度等效、受壓屈曲彎曲剛度等效、線性速度場下動能等效以及桿的彎曲振動固有特性近似等效,可描述細長桿大范圍運動、軸向拉壓變形和橫向彎曲變形耦合動力學特性.下一步,將該降階模型運用到張拉整體結(jié)構(gòu)的動力學建模,建立并數(shù)值求解系統(tǒng)動力學方程.

2 張拉整體系統(tǒng)瞬態(tài)動力學方程

根據(jù)圖6 所示離散桿的變形構(gòu)型,其應(yīng)變能可表達為

圖6 離散桿變形構(gòu)型示意圖Fig.6 Schematic representation of the deformed configuration for the discrete bar

其中,ΔLi(i=1,2,3,4)表示各彈簧的長度變化量,αj(j=1,2,3)表示各扭簧的轉(zhuǎn)角.

如圖6 所示,xk(k=1,2,3,4,5) 表示離散模型各節(jié)點在全局坐標系下的位置向量,為了用xk表示式(30)中的局部變量,定義如下各式

則式(30)中局部變量可通過以下公式計算

離散桿的應(yīng)變能式(30)對廣義坐標求偏導(dǎo)可得由桿內(nèi)力導(dǎo)致的第k個節(jié)點廣義內(nèi)力列陣表達為

同時,對于張拉整體結(jié)構(gòu),桿的兩個端點還受到繩索張力作用,設(shè)一根桿端點a與另一根桿端點b通過繩索相連,則由繩內(nèi)力導(dǎo)致的端點a廣義內(nèi)力列陣(端點b同理)表達為

其中,k為繩索的抗拉剛度,Δl=|xa-xb|-l0為繩索的長度變化量,l0為繩索的自然長度.

由于張拉整體結(jié)構(gòu)中一般采用輕質(zhì)繩索,故可忽略繩索慣性力影響.將所有桿各個節(jié)點全局位置向量、等效質(zhì)量、上述廣義內(nèi)力及受到的其他廣義外力進行組裝,得到系統(tǒng)廣義坐標列陣q、質(zhì)量矩陣M、廣義內(nèi)力列陣Fint及廣義外力列陣Fext.考慮系統(tǒng)可能受到約束,引入約束方程和拉格朗日乘子,可得系統(tǒng)在時刻t的瞬態(tài)動力學方程為

進一步,為了提高動力學計算效率,可將圖6 所示桿的內(nèi)部自由度x2,x3和x4縮聚至端部自由度x1和x5,即利用靜力縮聚法將內(nèi)部自由度剛度等效至端部自由度的廣義剛度矩陣[33-34],并將桿的總質(zhì)量均布在端部節(jié)點上,從而大大降低迭代過程中系統(tǒng)動力學方程(36)的維數(shù).當桿承受的壓力超過其屈曲臨界載荷時,需計算桿的后屈曲構(gòu)型,此時以計算得到端部位置坐標作為輸入,利用后屈曲內(nèi)力平衡關(guān)系式(4),迭代求解桿內(nèi)部節(jié)點的位置.

3 數(shù)值算例與實驗驗證

3.1 球形張拉整體結(jié)構(gòu)構(gòu)型與參數(shù)

圖7 為球形張拉整體結(jié)構(gòu)構(gòu)型圖,由6 根桿和24 根繩索組成.實驗中,制作了兩種不同的六桿球形張拉整體結(jié)構(gòu),一種由彈性模量較大的木桿和彈性模量較小的松緊繩構(gòu)成,另一種由彈性模量較小的橡膠桿和彈性模量較大的尼龍繩構(gòu)成.兩種結(jié)構(gòu)的幾何參數(shù)相同,每根桿長度均為0.2 m,橫截面均為圓形且直徑均為1 cm;木桿張拉整體結(jié)構(gòu)材料參數(shù)分別為:楊氏模量10 GPa,密度675 kg/m3,松緊繩拉伸剛度150 N/m,預(yù)拉伸4.5 N.橡膠桿張拉整體結(jié)構(gòu)材料參數(shù)分別為:楊氏模量19 MPa,密度1354 kg/m3,尼龍繩拉伸剛度30 kN/m,無預(yù)拉伸.

圖7 張拉整體結(jié)構(gòu)構(gòu)型圖(紅色粗線:桿,綠色細線:繩)Fig.7 Configuration diagram of tensegrity structure(thick red lines:bars,thin green lines:cables)

3.2 準靜態(tài)壓縮

采用萬能試驗機(Instron E44.104),對上述兩種張拉整體結(jié)構(gòu)進行準靜態(tài)壓縮實驗.將張拉整體結(jié)構(gòu)的底端固定在試驗機平臺上,對上端節(jié)點進行緩慢壓縮,測得壓力和壓縮位移的關(guān)系,并通過數(shù)值模型模擬了這一過程.圖8 所示為橡膠桿張拉整體結(jié)構(gòu)在不同壓縮位移下仿真和實驗構(gòu)型圖,隨著壓縮位移的增大,橡膠桿由純壓縮變形轉(zhuǎn)變?yōu)榍冃?

圖8 仿真(左)和實驗(右)中橡膠桿張拉整體結(jié)構(gòu)不同壓縮位移下變形圖Fig.8 Deformed configurations of the rubber tensegrity under difference compressions in the simulation(left)and experiment(right)

圖9 所示為壓力與壓縮位移曲線的實驗與仿真對比,其中實驗數(shù)據(jù)表現(xiàn)為3 次實驗測量的平均值與標準差.仿真結(jié)果與實驗結(jié)果基本吻合,且實驗數(shù)據(jù)的標準差較小,說明實驗具有較好的可重復(fù)性.圖9(a)中,木桿張拉整體結(jié)構(gòu)力-位移曲線呈現(xiàn)近似線性,由于松緊繩模量遠小于木桿模量,所以壓縮時木桿幾乎不發(fā)生變形,系統(tǒng)的彈性勢能主要儲存在松緊繩內(nèi),且桿不發(fā)生屈曲.圖9(b)中,橡膠桿張拉整體結(jié)構(gòu)力-位移曲線呈現(xiàn)明顯非線性,由于尼龍繩模量遠大于橡膠桿模量,所以壓縮時尼龍繩幾乎不伸長,系統(tǒng)的彈性勢能主要存儲在橡膠桿內(nèi);當壓縮位移達到8 mm 左右時,曲線出現(xiàn)了拐點,此時橡膠桿開始屈曲,在后屈曲段張拉整體結(jié)構(gòu)雖然剛度有所下降,但仍具有較高的承載能力.

圖9 壓力與壓縮位移的實驗與仿真對比圖:(a)木桿張拉整體結(jié)構(gòu);(b)橡膠桿張拉整體結(jié)構(gòu)Fig.9 Comparison for the compressive forces versus displacements of experiments and simulations:(a)Wooden tensegrity;(b)rubber tensegrity

為進一步理解張拉整體結(jié)構(gòu)靜態(tài)承載能力變化規(guī)律,開展不同參數(shù)下準靜態(tài)壓縮數(shù)值實驗.改變桿件的截面半徑、模量和繩索的預(yù)張力,得到如圖10 所示不同的力-位移曲線.從圖10 可以看出,桿件屈曲之前,結(jié)構(gòu)的整體剛度隨位移的增大而增大,而桿件屈曲之后,結(jié)構(gòu)的整體剛度趨于不變.如圖10(a) 和圖10(b) 所示,桿件的半徑和模量改變時,在桿件屈曲前對結(jié)構(gòu)的整體剛度影響較小,在桿件屈曲后對結(jié)構(gòu)的整體剛度影響較大.如圖10(c) 所示,當繩索預(yù)張力改變時,在桿件屈曲前對結(jié)構(gòu)的整體剛度影響較大,在桿件屈曲后對結(jié)構(gòu)的整體剛度基本沒有影響.

圖10 壓力與壓縮位移曲線在不同結(jié)構(gòu)參數(shù)下的對比:(a)桿截面半徑變化;(b)桿模量變化;(c)繩索預(yù)張力變化Fig.10 Comparison for the compressive forces versus displacements under different structural parameters:(a)Cross-section radii of bars;(b)moduli of bars;(c)pretensions of strings

圖10 壓力與壓縮位移曲線在不同結(jié)構(gòu)參數(shù)下的對比:(a)桿截面半徑變化;(b)桿模量變化;(c)繩索預(yù)張力變化(續(xù))Fig.10 Comparison for the compressive forces versus displacements under different structural parameters:(a)Cross-section radii of bars;(b)moduli of bars;(c)pretensions of strings(continued)

為驗證本文離散模型的高效性,采用ABAQUS商用有限元分析軟件建立張拉整體結(jié)構(gòu)模型,與本文等效模型的仿真時間進行對比.計算機配置為4 核(Intel(R) Xeon(R) CPU E3-1226 v3 @ 3.30GHz)、8GB內(nèi)存,采用單線程運行,時間步長為1.0×10-5s.圖11 為ABAQUS 單元收斂性測試,桿單元類型為B31,由圖可知每根桿取10 個單元與取30 個單元的力-位移曲線基本重合,說明此時計算收斂,故ABAQUS 模型中每根桿取10 個單元.與圖3 壓縮單根桿不同,用ABAQUS 壓縮張拉整體結(jié)構(gòu)時每根桿需引入較大初始缺陷才能保證算法收斂,故在桿屈曲前與本文離散模型存在一定誤差,但在屈曲后與本文離散模型基本一致.整個計算過程本文離散模型所需時間為120 s,ABAQUS 所需時間為426 s,在相同計算精度下計算效率提高了3 倍以上,說明本文提出的等效方法具有較高的計算效率.此外,本算例也說明,采用有限單元分析張拉整體結(jié)構(gòu)后屈曲過程時,需引入較大的初始幾何缺陷,由缺陷敏感性造成仿真結(jié)果可能存在較大偏差,而本文等效模型則不需引入過大的初始幾何缺陷.

圖11 ABAQUS 有限元分析收斂性測試和與本文離散桿模型的計算結(jié)果對比Fig.11 Convergence test of finite element analysis by ABAQUS and comparison with the computation result by the proposed discrete bar

3.3 模態(tài)分析

本節(jié)對張拉整體結(jié)構(gòu)的模態(tài)進行了實驗和仿真分析.將底部3 個節(jié)點固定在振動試驗臺上,對結(jié)構(gòu)進行掃頻測試,用高速相機對其振動進行拍攝.當頻率為15 Hz 時,張拉整體結(jié)構(gòu)發(fā)生較大幅度的上下晃動,此時上端節(jié)點的振幅達到最大,由此確定張拉整體結(jié)構(gòu)的基頻約為15 Hz,由數(shù)值仿真求得的基頻為15.37 Hz,兩者基本一致.實驗和仿真對應(yīng)的第一階振型如圖12 所示,可以看到上部3 根桿有明顯向上擴張的趨勢.

圖12 木桿張拉整體結(jié)構(gòu)第一階固有振型的實驗(左)和仿真(右)圖Fig.12 The first natural mode of vibration of the wooden tensegrity in the experiment(left)and simulation(right)

進一步,研究繩索與桿件不同剛度比對張拉整體結(jié)構(gòu)固有振動特性的影響規(guī)律.改變繩索抗拉剛度,使其與桿件等效抗壓剛度的比值ksb分別取0.01,1 和100,通過數(shù)值計算畫出前六階固有振型圖如圖13 所示.當兩者比值為0.01 時,繩索的剛度較小,其固有振型更容易被激發(fā),如圖13(a) 所示;當兩者比值為1 時,彈簧和桿件的固有頻率被同時激發(fā),如圖13(b);當兩者比值為100 時,桿件的固有頻率被優(yōu)先激發(fā),如圖13(c).因此,張拉整體結(jié)構(gòu)的動態(tài)力學特性具有可設(shè)計性.

圖13 不同繩索和桿剛度比ksb 下前六階固有振型圖Fig.13 The first six natural modes of vibration under different stiffness ratios of strings and bars ksb

3.4 碰撞動力學分析

研究兩種張拉整體結(jié)構(gòu)自由下落后與平面碰撞的動力學行為,進行仿真和實驗的對比分析.實驗中,張拉整體結(jié)構(gòu)底端節(jié)點與固定在實驗平臺上的壓電傳感器發(fā)生碰撞,通過壓電傳感器可測得碰撞力,采用高速相機對結(jié)構(gòu)運動過程進行拍攝,對圖像數(shù)據(jù)進行處理可得結(jié)構(gòu)頂端節(jié)點的位移.圖14 所示為橡膠桿張拉整體結(jié)構(gòu)在不同碰撞時刻的仿真和實驗構(gòu)型圖,碰撞過程中橡膠桿會發(fā)生屈曲以緩沖能量.

圖14 橡膠桿張拉整體結(jié)構(gòu)碰撞仿真(左)和實驗(右)的動態(tài)構(gòu)型圖Fig.14 Dynamic configurations of the rubber tensegrity in the simulation(left)and experiment(right)of impact

圖14 橡膠桿張拉整體結(jié)構(gòu)碰撞仿真(左)和實驗(右)的動態(tài)構(gòu)型圖(續(xù))Fig.14 Dynamic configurations of the rubber tensegrity in the simulation(left)and experiment(right)of impact(continued)

圖15 為兩種張拉整體結(jié)構(gòu)底端節(jié)點碰撞力和頂端節(jié)點位移隨時間變化曲線在不同碰撞初速度下仿真與實驗對比圖.從圖15 中可看出,仿真和實驗的曲線基本一致,碰撞力和位移都呈現(xiàn)先增大后減小的趨勢.當碰撞力減小為零時,結(jié)構(gòu)底端節(jié)點開始與地面分離反彈,而由于彈性波在結(jié)構(gòu)中的傳播效應(yīng),結(jié)構(gòu)頂端節(jié)點反彈時刻遠滯后于底端節(jié)點反彈時刻.此外,通過圖15(a)、圖15(b)和圖15(c)、圖15(d)對比可知,橡膠桿張拉整體結(jié)構(gòu)較木桿張拉整體結(jié)構(gòu)具有更長的動力學響應(yīng)時間,可見其柔性更高.由圖15(c)和圖15(d)可知,雖然在碰撞過程中橡膠桿會發(fā)生屈曲,但結(jié)構(gòu)的整體剛度沒有發(fā)生太大變化,結(jié)構(gòu)仍然具有較高的抗沖擊性.圖15(d)顯示,橡膠桿張拉整體結(jié)構(gòu)在v=2.97 m/s 時,仿真結(jié)果與實驗結(jié)果誤差較為明顯,考慮到橡膠為超彈性材料,其彈性模量隨應(yīng)變增大會先增大,后減小,最后再增大.當碰撞速度較小時,橡膠桿取小應(yīng)變階段彈性模量,此時彈性模量較大,與實驗結(jié)果吻合較好.隨碰撞速度增大,橡膠桿的應(yīng)變也隨之增大,彈性模量變小,但在計算時橡膠桿仍取小應(yīng)變階段較大的彈性模量,因此與實驗值相比,其頂端節(jié)點位移變化幅度較小.

圖15 張拉整體結(jié)構(gòu)碰撞實驗與仿真結(jié)果對比Fig.15 Comparison for the experimental and simulation results of impact for tensegrity structures

圖15 張拉整體結(jié)構(gòu)碰撞實驗與仿真結(jié)果對比(續(xù))Fig.15 Comparison for the experimental and simulation results of impact for tensegrity structures(continued)

4 結(jié)論

本文提出了張拉整體結(jié)構(gòu)的動力學等效建模與計算方法,通過理論推導(dǎo)、建模仿真及實驗驗證,得到以下結(jié)論:

(1) 提出了用五節(jié)點彈/扭簧集中質(zhì)量的離散模型來等效連續(xù)桿力學特性,通過靜力學和動能等效推導(dǎo)了離散模型等效剛度和等效質(zhì)量表達式.通過彎曲振動固有特性等效確定了彎曲剛度分布參數(shù)n和質(zhì)量分布參數(shù)c,使離散模型前兩階固有頻率相對誤差降到1%以內(nèi),提高了系統(tǒng)動力學響應(yīng)預(yù)測的準確性.

(2)通過對準靜態(tài)壓縮實驗和仿真進行對比,驗證了本文模型預(yù)測張拉整體結(jié)構(gòu)靜力學行為的準確性.對橡膠桿張拉整體結(jié)構(gòu)承載能力進行了分析,桿件半徑和模量的變化,對桿件屈曲前結(jié)構(gòu)整體剛度影響較小,對桿件后屈曲段結(jié)構(gòu)整體剛度影響較大,而預(yù)張力變化對結(jié)構(gòu)整體剛度影響則相反.

(3)對張拉整體結(jié)構(gòu)的模態(tài)進行了實驗和仿真分析,驗證了本文模型預(yù)測張拉整體結(jié)構(gòu)固有振動特性的準確性.分析了繩索與桿軸向剛度比變化對張拉整體結(jié)構(gòu)前六階模態(tài)的影響,發(fā)現(xiàn)當繩索剛度低于桿件剛度時,繩索的局部固有振型被優(yōu)先激發(fā),反之,則桿件的局部固有振型被優(yōu)先激發(fā),即張拉整體結(jié)構(gòu)的固有動力學特性具有可設(shè)計性.

(4)對張拉整體結(jié)構(gòu)與剛性平面的碰撞過程進行了實驗和仿真分析,驗證了本文模型預(yù)測張拉整體結(jié)構(gòu)瞬態(tài)動力學行為的準確性.證明了當桿件屈曲時,六桿球形張拉整體結(jié)構(gòu)仍具有較高抗沖擊剛度.因此當柔性張拉整體結(jié)構(gòu)用作行星著陸器時可擴大其承載的設(shè)計域.

本文提出的動力學等效建模和計算方法,將有望為以張拉整體結(jié)構(gòu)為元胞的行星探測器、大型可展開結(jié)構(gòu)及點陣材料等復(fù)雜柔性系統(tǒng)的動力學高效分析與控制提供理論和仿真基礎(chǔ).

附錄

對于細長桿(π2I/(AL2) ?1) 的離散模型,式(21) 表示的其橫向振動剛度矩陣各分量表達式為