張拉整體結(jié)構(gòu)的動(dòng)力學(xué)等效建模與實(shí)驗(yàn)驗(yàn)證1)

陳占魁 羅 凱 田 強(qiáng)

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

引言

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

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

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

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

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

1 屈曲細(xì)長(zhǎng)桿的動(dòng)力學(xué)等效

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

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

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

1.2 剛度大小等效

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

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

其中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為廣義坐標(biāo),由最小勢(shì)能原理可得以下平衡方程

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

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

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

式(6)的后兩個(gè)方程可改寫(xiě)為

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

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

兩端簡(jiǎn)支桿的Euler 屈曲臨界載荷為

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

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

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

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

1.3 質(zhì)量大小等效

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

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

另一方面,通過(guò)線速度場(chǎng)u和v,離散模型的動(dòng)能可表達(dá)為

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

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

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

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

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

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

當(dāng)桿橫向運(yùn)動(dòng)在水平參考線附近進(jìn)行線性化時(shí),各轉(zhuǎn)角趨于零,有

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

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

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

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

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

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

其中

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值算例與實(shí)驗(yàn)驗(yàn)證

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

圖7 為球形張拉整體結(jié)構(gòu)構(gòu)型圖,由6 根桿和24 根繩索組成.實(shí)驗(yàn)中,制作了兩種不同的六桿球形張拉整體結(jié)構(gòu),一種由彈性模量較大的木桿和彈性模量較小的松緊繩構(gòu)成,另一種由彈性模量較小的橡膠桿和彈性模量較大的尼龍繩構(gòu)成.兩種結(jié)構(gòu)的幾何參數(shù)相同,每根桿長(zhǎng)度均為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,無(wú)預(yù)拉伸.

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

3.2 準(zhǔn)靜態(tài)壓縮

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

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

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

圖9 壓力與壓縮位移的實(shí)驗(yàn)與仿真對(duì)比圖:(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

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

圖10 壓力與壓縮位移曲線在不同結(jié)構(gòu)參數(shù)下的對(duì)比:(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ù)下的對(duì)比:(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)

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

圖11 ABAQUS 有限元分析收斂性測(cè)試和與本文離散桿模型的計(jì)算結(jié)果對(duì)比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é)對(duì)張拉整體結(jié)構(gòu)的模態(tài)進(jìn)行了實(shí)驗(yàn)和仿真分析.將底部3 個(gè)節(jié)點(diǎn)固定在振動(dòng)試驗(yàn)臺(tái)上,對(duì)結(jié)構(gòu)進(jìn)行掃頻測(cè)試,用高速相機(jī)對(duì)其振動(dòng)進(jìn)行拍攝.當(dāng)頻率為15 Hz 時(shí),張拉整體結(jié)構(gòu)發(fā)生較大幅度的上下晃動(dòng),此時(shí)上端節(jié)點(diǎn)的振幅達(dá)到最大,由此確定張拉整體結(jié)構(gòu)的基頻約為15 Hz,由數(shù)值仿真求得的基頻為15.37 Hz,兩者基本一致.實(shí)驗(yàn)和仿真對(duì)應(yīng)的第一階振型如圖12 所示,可以看到上部3 根桿有明顯向上擴(kuò)張的趨勢(shì).

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

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

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

3.4 碰撞動(dòng)力學(xué)分析

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

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

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

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

4 結(jié)論

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

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

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

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

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

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

附錄

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