考慮非線性本構(gòu)的非剛性折紙結(jié)構(gòu)動力學(xué)建模與分析1)

袁婷婷 任昆明 方雨橋 劉錦陽

(上海交通大學(xué)工程力學(xué)系,水動力學(xué)教育部重點實驗室,上海 200240)

引言

折紙(origami)起源于傳統(tǒng)藝術(shù)活動,泛指將平面二維材料通過折疊轉(zhuǎn)變?yōu)榭臻g三維結(jié)構(gòu)的過程[1-3].近年來,折紙已發(fā)展成為一門新興的學(xué)科,在航天[4-6]、生物醫(yī)學(xué)[7-9]、建筑[10-12]、機(jī)器人[13-16]、材料科學(xué)[17-19]等諸多領(lǐng)域中發(fā)揮重要作用,這些研究利用了折紙結(jié)構(gòu)良好的可折疊性、形狀多樣性、輕量小型等優(yōu)點[20],與現(xiàn)代科學(xué)有機(jī)結(jié)合,表現(xiàn)出交叉性和復(fù)雜性的特點.隨著信息技術(shù)、生物技術(shù)、新能源技術(shù)、新材料技術(shù)的交叉融合,折紙結(jié)構(gòu)的研究得到了長足的發(fā)展,主要表現(xiàn)在幾何設(shè)計、運(yùn)動學(xué)分析和準(zhǔn)靜態(tài)分析等方面[21-23],但迄今為止,對折紙結(jié)構(gòu)動力學(xué)研究尚未成熟,在動力學(xué)建模方法和理論分析上仍面臨較大的困難和挑戰(zhàn)[24].為了深度挖掘折紙結(jié)構(gòu)的潛藏價值,開展折紙結(jié)構(gòu)的動力學(xué)研究勢在必行.

剛性折紙作為折紙的一個特殊子集,其數(shù)學(xué)模型的建立通?；谝韵录僭O(shè): 折痕是直線折痕,折疊面視為剛性,即不會發(fā)生拉伸或彎曲[25].目前,對剛性折紙的研究已經(jīng)較為成熟,并已開發(fā)了多個軟件[26-27].隨著材料科學(xué)技術(shù)的不斷發(fā)展,折紙結(jié)構(gòu)在工程中的應(yīng)用更為廣泛,在很多情況下,折紙結(jié)構(gòu)的折疊面會發(fā)生拉伸、彎曲或剪切等變形,甚至在展開過程中還會誘發(fā)接觸碰撞等問題.因此,針對這類非剛性折紙結(jié)構(gòu),用剛性折紙模型進(jìn)行模擬將無法反映非剛性折紙?zhí)赜械牧W(xué)性能.

考慮非剛性折紙結(jié)構(gòu)的柔性效應(yīng),采用基于殼單元的有限元方法分析[28-30],是最為直接的一類建模方法.Yuan 等[6,31]基于絕對節(jié)點坐標(biāo)的薄板單元,對具有大變形的折紙結(jié)構(gòu)展開動力學(xué)中接觸碰撞問題進(jìn)行了研究,給出了折紙結(jié)構(gòu)展開過程中的動力學(xué)響應(yīng)和應(yīng)力分布情況.有限元方法雖然能夠給出包括應(yīng)力分布在內(nèi)的局部特征,但建模復(fù)雜度高,計算耗時長,較難推廣到折紙結(jié)構(gòu)動力學(xué)的性能優(yōu)化和控制等工作.為此,有必要根據(jù)研究目的簡化折紙結(jié)構(gòu)的次要特征,建立出等效的非剛性折紙結(jié)構(gòu)的動力學(xué)模型.近年來,國內(nèi)外學(xué)者針對等效非剛性折紙結(jié)構(gòu)開展了相關(guān)工作.文獻(xiàn)[32]將折紙結(jié)構(gòu)簡化為彈簧–質(zhì)點系統(tǒng),即折疊面頂點為質(zhì)點,折痕用彈簧代替.復(fù)旦大學(xué)方虹斌等[33-36]采用非線性彈簧等效動力學(xué)建模方法,并結(jié)合實驗,開展了不同折紙結(jié)構(gòu)和折紙超材料的動力學(xué)研究.對于另一類等效空間桁架模型,Yasuda 等[37]將非剛性折紙結(jié)構(gòu)簡化為一維運(yùn)動的連桿結(jié)構(gòu),并結(jié)合最小勢能原理,揭示了三角形圓柱折紙結(jié)構(gòu)的雙阱勢能的特征.文獻(xiàn)[38]較早提出桿–鏈模型,將折紙結(jié)構(gòu)三角化為受旋轉(zhuǎn)鉸約束的桁架結(jié)構(gòu),開展運(yùn)動學(xué)分析.在桿–鏈模型基礎(chǔ)上,分別給桿和旋轉(zhuǎn)鉸分配恒定剛度,Wei 等[39]采用顯示時間積分算法來模擬三浦折紙結(jié)構(gòu)的彎曲.Fuchi 等[40-41]采用線性桿–鏈模型,實現(xiàn)對折紙超材料的拓?fù)鋬?yōu)化.為了模擬折紙結(jié)構(gòu)的大變形和位移,文獻(xiàn)[42]基于超彈性本構(gòu)的桿單元,提出了一種非線性桿–鏈模型,并對多種折紙結(jié)構(gòu)進(jìn)行準(zhǔn)靜態(tài)分析.Filipov 等[43]進(jìn)一步細(xì)化離散方式,提升了非線性桿–鏈模型對非剛性折紙結(jié)構(gòu)面外變形描述的精度.Dong 等[44]提出了一種將粒子–桿–彈簧模型與有限粒子法相結(jié)合的動力學(xué)模型,由折痕存儲的能量驅(qū)動折紙結(jié)構(gòu)展開.目前針對非線性本構(gòu)的非剛性折紙結(jié)構(gòu)動力學(xué)工作較少.為了分析超彈性非剛性折紙結(jié)構(gòu)的動力學(xué)特性,需要建立一種通用且高效的動力學(xué)模型.

本文的工作重點是建立一種非剛性折紙結(jié)構(gòu)的桿–鏈動力學(xué)模型,可適用于大變形和作大范圍運(yùn)動的折紙結(jié)構(gòu).第一節(jié)詳細(xì)給出了非剛性折紙結(jié)構(gòu)動力學(xué)建模過程.首先推導(dǎo)具有超彈性本構(gòu)的桿單元的彈性內(nèi)力陣.提出一種改進(jìn)的卷簧本構(gòu)模型,推導(dǎo)非線性卷簧單元的彈性力陣.最后,考慮阻尼效應(yīng),建立非剛性折紙結(jié)構(gòu)多體系統(tǒng)動力學(xué)方程.第二節(jié),將本文提出的非剛性折紙結(jié)構(gòu)動力學(xué)模型分別應(yīng)用于單菱形折疊、葉內(nèi)折疊和Kresling 折疊這三種折紙結(jié)構(gòu)中,并依次進(jìn)行數(shù)值模擬,驗證了桿–鏈動力學(xué)模型的準(zhǔn)確性和高效性.在此基礎(chǔ)上開展多穩(wěn)態(tài)、瞬態(tài)動力學(xué)和波動力學(xué)特性分析.第三節(jié)對本文的研究工作進(jìn)行總結(jié),并對未來工作提出展望.

1 非剛性折紙結(jié)構(gòu)動力學(xué)建模

本文將非剛性折紙結(jié)構(gòu)簡化為帶卷簧的空間桁架結(jié)構(gòu),建立了桿–鏈模型,即將折痕和虛擬折痕等效為桿單元,桿單元相交的節(jié)點處用球鉸連接,卷簧沿著桿的方向,如圖1 所示.桿–鏈模型將折紙結(jié)構(gòu)板面的拉伸和壓縮用桿單元的變形來描述,采用卷簧來表征折痕的抗彎作用和折疊面的彎曲作用,該模型能夠縮減計算規(guī)模,可高效地對非剛性折紙結(jié)構(gòu)進(jìn)行仿真,研究其豐富的力學(xué)行為.

圖1 桿–鏈模型示意圖Fig.1 The bar-and-hinge model

1.1 超彈性本構(gòu)的桿單元

任意一個非剛性折紙結(jié)構(gòu)均可簡化為由一系列柔性桿件通過鉸鏈組合起來的多體系統(tǒng).為了表述方便,定義在慣性基下任一節(jié)點p的位置矢量坐標(biāo)陣為xp,變形前的位置矢量坐標(biāo)陣為Xp,節(jié)點位移矢量坐標(biāo)陣為up=xp-Xp.

如圖1 所示,桿單元e的廣義坐標(biāo)陣為

其中,xi和xj分別表示i節(jié)點和j節(jié)點在慣性基下的位置矢量坐標(biāo)陣.取桿單元未變形時長度為Le,單元橫截面積為Ae,則在局部坐標(biāo)系下桿單元的彈性變形能[42]為

其中,W為應(yīng)變能密度.由于桿單元是一維的,只需要考慮應(yīng)力張量和應(yīng)變張量的軸向分量.采用線性形函數(shù),可以得到一維格林–拉格朗日應(yīng)變 εx與單元節(jié)點位移ue的關(guān)系式[45]為

式中,ue,B1,B2可分別表示為

其中,e1=[1 0 0],I3是3×3 的單位矩陣.

設(shè)σx為沿x方向的Piola-Kirchhoff 應(yīng)力分量,應(yīng)變能密度的變分可表示為 δW=σxδεx.桿單元彈性變形能的變分為

將上式代入式(5),并得到桿單元的內(nèi)力陣為

定義一維線性模量為

本文采用具有超彈性本構(gòu)的桿單元,可適用于非線性材料,例如橡膠、硅膠等,能夠處理大變形問題.本文采用Ogden 超彈性模型[46,42],其應(yīng)變能密度為

其中,λ1,λ2,λ3分別代表右應(yīng)變張量的三個特征值,F是變形梯度.N,μi,αi是材料參數(shù).對于一維桿單元僅與λ1相關(guān),且有

一維二階P-K 應(yīng)力σx和一維切線模量C由下式給出

考慮到材料未變形時(λ1=1)為無應(yīng)變無應(yīng)力狀態(tài),有如下約束關(guān)系

將式(15)代入式(14),可得未變形時一維模量為

在本文模型中,取N=2,給定α1,α2和C0,可根據(jù)式(15)和式(16)確定材料參數(shù)μ1和μ2.

假設(shè)桿–鏈模型中的桿單元為無質(zhì)量桿,則將三角形折疊面的質(zhì)量m均勻分布到三個節(jié)點上,即每個節(jié)點的質(zhì)量為m/3,如圖1 所示.三角形單元的質(zhì)量陣為

其中,I9是9×9 的單位矩陣.

1.2 非線性卷簧單元

非剛性折紙結(jié)構(gòu)的折痕,以及為描述折疊面彎曲添加的虛擬折痕,均采用非線性卷簧單元來建模.如圖2 所示,卷簧沿著桿的方向,采用二面角進(jìn)行度量[42],一個二面角至少需要三個向量才能描述,因此一個卷簧單元包括4 個節(jié)點和5 根桿.卷簧單元的廣義坐標(biāo)坐標(biāo)陣為

圖2 卷簧單元示意圖Fig.2 Rotational spring element

其中,xk和xl分別表示k節(jié)點和l節(jié)點的位置矢量坐標(biāo)陣.則兩個平面的法向量分別為

其中,定義rpq=xp-xq為相對位置矢量坐標(biāo)陣,為rpq對應(yīng)的坐標(biāo)方陣.二面角 θ ∈[0,2π),可以表示為

式中sgn 為符號函數(shù).

將卷簧單元變形能對卷簧單元廣義坐標(biāo)陣進(jìn)行變分,可以得到卷簧單元的內(nèi)力陣為

對卷簧單元的內(nèi)力陣求變分,可得卷簧單元切向剛度陣

其中,k=?M/?θ 為卷簧切線剛度.

本文采用的卷簧為非線性卷簧,其特點是當(dāng)二面角 θ →0,2π 時,力矩M→∞,能夠避免折疊面的穿透問題.相較于傳統(tǒng)的接觸時施加約束的方法,此種考慮接觸的方法更簡潔方便,在數(shù)值仿真時不需要做其他調(diào)整,只需要建立合適的卷簧單元數(shù)學(xué)模型即可.

考慮到卷簧力矩的連續(xù)性,本文改進(jìn)了非線性卷簧的本構(gòu)模型,相較于傳統(tǒng)的卷簧本構(gòu)模型[42],該改進(jìn)的本構(gòu)模型在臨近接觸時延長了過渡區(qū),如圖3 所示.采用改進(jìn)的本構(gòu)模型求解動力學(xué)接觸問題時,能夠有效避免因時間積分步過大導(dǎo)致的穿透問題,則改進(jìn)的本構(gòu)模型中卷簧彎矩的具體表達(dá)式如下

圖3 非線性卷簧單元的本構(gòu)模型: 力矩與二面角的關(guān)系Fig.3 Nonlinear constitutive model of rotational spring: Dihedral angle θ versus moment M

其中,Ls為卷簧單元未變形時的長度,k0為卷簧線性部分的切線剛度.本文的數(shù)值算例中,均取a=100,b=0.4,c=0.8.

1.3 非剛性折紙結(jié)構(gòu)動力學(xué)方程

結(jié)合多體系統(tǒng)的約束條件 Φ=0,引入拉格朗日乘子陣λ,整理得到約束多體系統(tǒng)的動力學(xué)方程

其中,系統(tǒng)的結(jié)構(gòu)阻尼陣為C=αM,Φq=?Φ/?q為約束雅克比矩陣.

本文采用廣義-α法[47]對式(28)的指標(biāo)3 微分–代數(shù)方程組進(jìn)行隱式求解.考慮到材料和卷簧的非線性本構(gòu),為了提高計算效率,采用了變步長[48-49]與預(yù)條件處理等計算策略.

2 數(shù)值算例仿真

將本文提出的非線性桿–鏈模型應(yīng)用于三種經(jīng)典的非剛性折紙結(jié)構(gòu)的動力學(xué)分析中.第一個算例是進(jìn)行單菱形折紙結(jié)構(gòu)的動力學(xué)仿真,將桿–鉸鏈模型分別與解析解、準(zhǔn)靜態(tài)以及有限元絕對節(jié)點坐標(biāo)法ANCF[31]進(jìn)行對比,驗證桿-鉸鏈模型的準(zhǔn)確性和高效性.第二個算例是對葉內(nèi)折疊折紙結(jié)構(gòu)進(jìn)行動力學(xué)仿真,并將剛性模型和非剛性模型的動力學(xué)結(jié)果進(jìn)行對比分析.第三個算例是對Kresling 折紙結(jié)構(gòu)進(jìn)行仿真,揭示了折紙結(jié)構(gòu)的多穩(wěn)態(tài)和動力學(xué)特性.

2.1 單菱形折疊

為了驗證桿–鏈動力學(xué)模型的正確性,選取的研究對象是一個簡單的菱形折疊,包含四個節(jié)點,兩個三角形面板和一條折痕,如圖4 所示.菱形折疊中各節(jié)點的初始位置分別為:A(-0.07,0,0) m,B(0,-0.025,0) m,C(0,0.025,0)m,D(0.07,0,0) m,初始折疊角θ0=180o.此外,節(jié)點A,B,C均與慣性基固結(jié),D點受到豎直向上的外力Fext,三角形板BDC會沿著折痕向上翻折.

圖4 單菱形折疊示意圖Fig.4 A rhombus fold

采用桿–鏈模型,超彈性桿單元的材料參數(shù)分別為: 未變形時一維模量C0=1×107Pa,橫截面積A=3×10-5m2,材料本構(gòu)選取Ogden 超彈性模型,N=2,α1=5,α2=1.非線性卷簧單元的參數(shù)為:k0=0.012 5 N·m/rad,θ1=45o,θ2=135o.三角形面板的密度為1×103kg/m3,厚度為1 mm,忽略系統(tǒng)的結(jié)構(gòu)阻尼.根據(jù)算例工況,在忽略重力的情況下,僅D點受到豎直向上的外力作用時,通過受力分析,得到外力載荷Fext與折疊角度θ的解析關(guān)系式為

圖5 給出了折疊角θ與外力載荷Fext的關(guān)系曲線.通過對比理論解和桿–鏈準(zhǔn)靜態(tài)的結(jié)果,Fext由0 逐步加載至1 N,可以觀察到兩組數(shù)據(jù)完全吻合,驗證了桿–鏈模型的準(zhǔn)確性.

對于動力學(xué)仿真,采用外載荷Fext=t2/9(N),計算0 s～3 s 時間內(nèi)單菱形折疊的運(yùn)動過程.觀察圖5發(fā)現(xiàn),桿–鏈動力學(xué)模型的結(jié)果在準(zhǔn)靜態(tài)結(jié)果附近振蕩較為明顯,說明在動載荷作用下,對折紙結(jié)構(gòu)進(jìn)行動力學(xué)仿真揭示真實的動態(tài)特性的必要性.為了驗證桿–鏈動力學(xué)模型的準(zhǔn)確性,將該模型的結(jié)果與基于絕對節(jié)點坐標(biāo)法的動力學(xué)結(jié)果進(jìn)行對比,可以觀察到兩種方法的結(jié)果曲線基本重合.需要指出的是在本算例給定的參數(shù)下,折疊面不會發(fā)生明顯的彎曲變形.接下來對比兩種模型的計算效率,發(fā)現(xiàn)在相同計算機(jī)設(shè)置下,桿–鏈動力學(xué)模型的計算時間僅是ANCF 動力學(xué)模型的1/3,體現(xiàn)了桿–鏈動力學(xué)模型的高效性.

圖5 單菱形折疊: 折疊角與外力的關(guān)系圖Fig.5 A rhombus fold: Fext versus θ

綜上所述,在折疊面未發(fā)生大變形的情況下,本文提出的桿–鏈動力學(xué)模型的精度足以反映折紙結(jié)構(gòu)的動力學(xué)特性.此外,桿–鏈動力學(xué)模型的計算效率較高,有利于開展非剛性折紙結(jié)構(gòu)的優(yōu)化和設(shè)計等工作.

2.2 基于葉內(nèi)折疊的折紙結(jié)構(gòu)

大部分折紙的折疊方式靈感來源于自然界,De Focatiis 等[50]通過觀察樹葉的展開紋理,提出了葉內(nèi)折疊的折疊法,該折疊圖案的特點是具有很大的拉伸比,能夠有效減少體積,可節(jié)省存儲空間.本節(jié)將對基于葉內(nèi)折疊的折紙結(jié)構(gòu)動力學(xué)進(jìn)行研究.

圖6 所示為葉內(nèi)折疊的折痕分布,其中直線表示山折,曲線表示谷折.將紙面固定在O-xy平面上,山折是將紙的兩側(cè)向z軸負(fù)方向折,折痕形成凸起;谷折是將紙的兩側(cè)向z軸正方向折,折痕形成凹陷.圖6 右側(cè)給出了葉內(nèi)折疊折紙展開過程的實物圖.

圖6 葉內(nèi)折疊的折痕分布Fig.6 Crease pattern of leaf-in origami

2.2.1 剛性折紙結(jié)構(gòu)動力學(xué)分析

為了突出本文提出的非剛性折紙模型的優(yōu)勢,可以先對剛性的葉內(nèi)折疊結(jié)構(gòu)進(jìn)行動力學(xué)分析,其中剛性折紙模型可參考筆者前期的工作[31],本文不再贅述.對于方型的葉內(nèi)折疊,當(dāng)施加沿著正方形對角線的大小相同的驅(qū)動約束時,考慮對稱性,僅對1/8 構(gòu)型進(jìn)行研究即可.

根據(jù)葉內(nèi)折疊的幾何關(guān)系,以及剛性面任意兩點的等距約束,可知系統(tǒng)的自由度 δ=-1,可判斷出葉內(nèi)折疊具有剛性不可展特征,因此需要釋放至少2 個約束才能保證剛性展開.基于該前提,對葉內(nèi)折疊添加2 條虛擬折痕,使得自由度 δ=1.其中折痕方式有四種,如圖7 所示.

圖7 葉內(nèi)折疊添加虛擬折痕的方式Fig.7 Methods of adding virtual creases in leaf-in origami

設(shè)方型葉內(nèi)折疊的邊長為6cm,在G點施加的速度驅(qū)動,保證恰好在5 s 時從完全收攏到完全展開.在計算過程中,發(fā)現(xiàn)添加折痕的I,II 和III 方式在展開過程中均存在自鎖現(xiàn)象,導(dǎo)致計算中止,得到的鎖定構(gòu)型如圖8 所示.只有IV 方式的虛擬折痕添加方式,葉內(nèi)折疊能夠順利展開.

圖8 I,II,III 方式下葉內(nèi)折疊展開過程的奇異構(gòu)型Fig.8 Singular configuration of leaf-in origami during the deployment in I,II,III cases

分別對葉內(nèi)折疊展開和收攏過程進(jìn)行數(shù)值仿真.采用IV 折痕添加方式,葉內(nèi)折疊雖然在展開和收攏的計算過程中均未鎖定,但是完全展開的構(gòu)型是奇異構(gòu)型,從而導(dǎo)致無法以該構(gòu)型作為初始構(gòu)型進(jìn)行收攏計算.為此將t=4.99 s 時刻作為收攏計算的初值,避開奇異構(gòu)型,實現(xiàn)葉內(nèi)折疊收攏過程的計算.

如圖9 所示,將數(shù)值求解的結(jié)果與ADAMS 仿真的結(jié)果進(jìn)行對比,驗證了剛性折紙動力學(xué)模型的準(zhǔn)確性.在Step1 展開過程中,即t∈[0,5] s,剛性葉內(nèi)折疊的zA并不是隨時間單調(diào)遞增的,而是呈現(xiàn)先增加,再減小,最后再增加的趨勢.通過觀察展開過程的構(gòu)型圖發(fā)現(xiàn),最內(nèi)層向下翻折,zA減小,這是因為第三層在展開時受剛性折紙的等距約束的影響.對于Step2 收攏過程,即t∈(5,10] s,由于對收攏過程與展開過程施加的驅(qū)動大小相等,僅方向相反,可以觀察到剛性葉內(nèi)折疊的zA與展開過程呈對稱分布.

圖9 葉內(nèi)折疊展開和收攏過程中A 點z 坐標(biāo)時域圖Fig.9 Time history of z coordinate of point A of leaf-in origami during unfolding and folding process

針對剛性不可展問題,本文給出了添加虛擬折痕和修正初始構(gòu)型的方法,有效解決了展開和收攏過程的鎖定問題.本文對剛性折紙結(jié)構(gòu)的可展性和動力學(xué)研究,為工程問題中折紙結(jié)構(gòu)的運(yùn)動控制和設(shè)計優(yōu)化提供理論指導(dǎo).

2.2.2 非剛性折紙結(jié)構(gòu)動力學(xué)分析

本節(jié)將基于桿–鏈模型,對非剛性葉內(nèi)折疊折紙結(jié)構(gòu)的展開過程進(jìn)行動力學(xué)仿真.如圖10(a)給出了葉內(nèi)折疊的初始收攏構(gòu)型,并且四個角點P1,P2,P3,P4均受到沿對角線方向向外的位移驅(qū)動d(t),具體表達(dá)式為

其中,ad=0.276 cm/s2,s0=3 cm,t1=1 s,te=10 s.該參數(shù)能夠保證在t=5 s 時,葉內(nèi)折疊的剛性折紙模型恰好完全展開.

針對非剛性桿–鏈動力學(xué)模型,超彈性桿單元的材料參數(shù)分別為: 未變形時一維模量C0=1×107Pa,橫截面積A=1×10-5m2,材料本構(gòu)選取Ogden 超彈性模型,N=2,α1=5,α2=1.如圖10(b)所示,添加了彎曲折痕以描述折疊面的彎曲.設(shè)彎曲卷簧的折痕剛度為折痕卷簧剛度為0.001 N·m/rad,為了避免穿透現(xiàn)象發(fā)生,即要求二面角不能出現(xiàn)0 和2π,取θ1=ε,θ2=2π-ε,ε=0.000 5 rad.折疊面的密度為1420 kg/m3,厚度為0.5 mm,取系統(tǒng)的結(jié)構(gòu)阻尼系數(shù)為α=1.5.

圖10 葉內(nèi)折疊: (a)初始構(gòu)型;(b) 1/8 模型卷簧分布Fig.10 Leaf-in origami: (a) Initial configuration and (b) spring patter in 1/8 model

圖11 給出了在相同驅(qū)動下,分別采用剛性和非剛性折紙模型模擬的葉內(nèi)折疊展開過程的構(gòu)型圖.可以觀察到,剛性折紙模型在t=5 s 時,恰好能夠完全攤平,但受限于完全展開時處于奇異位置,計算中止.在非剛性折紙模型中,由于考慮了折疊桿的軸向變形以及折痕處卷簧剛度,展開過程有所滯后.在葉內(nèi)折疊完全展開后,持續(xù)施加位移驅(qū)動,非剛性的葉內(nèi)折疊的各桿單元發(fā)生明顯的軸向拉伸變形,并且變形程度由中心向外逐漸變大.

圖11 葉內(nèi)折疊: 剛性與非剛性折紙模型展開過程構(gòu)型圖Fig.11 Leaf-in origami: Configuration of the rigid and non-rigid origami models during the deployment

圖12 分別給出了A點z坐標(biāo)的位置、速度和加速度時間歷程圖.圖12(a)顯示,剛性折紙模型在展開過程中,會出現(xiàn)最內(nèi)層向下翻折的現(xiàn)象,即zA呈現(xiàn)先增加后減小的趨勢.但在非剛性折紙模型中,葉內(nèi)折疊因柔性效應(yīng)能直接拉開,zA單調(diào)增加到0,而后發(fā)生振動,最后振蕩逐漸衰減直至為0.由圖12(b)和圖12(c)可知,剛性折紙模型在接近t=5 s 時,速度與加速度會突增,無法收斂,臨近自鎖狀態(tài).對于非剛性折紙模型,觀察到在展開的初始階段數(shù)值出現(xiàn)振蕩,這是由于葉內(nèi)折疊在初始構(gòu)型時,折面存在接觸.由此可知,非剛性折紙模型相較于剛性折紙模型具有更好的數(shù)值模擬通用性,不僅能夠避免自鎖問題,揭示剛性折紙無法給出的動力學(xué)特性,還能給出具有大變形的張緊構(gòu)型,更貼合實際情況.

圖12 葉內(nèi)折疊: A 點z 坐標(biāo)的 (a)位置;(b)速度;(c)加速度時域圖Fig.12 Leaf-in origami: Time histories of z coordinate of point A of(a) position,(b) velocity and (c) acceleration

為了進(jìn)一步分析非剛性葉內(nèi)折疊折紙結(jié)構(gòu)展開過程的動力學(xué)特性,選取了1/8 模型中對角線上的A,E,F,G四個點為觀測點,分別給出z坐標(biāo)位置的時間歷程圖,如圖13 所示.可以觀察到,點G作為驅(qū)動作用點,在運(yùn)動過程z坐標(biāo)保持不變.在葉內(nèi)折疊展開過程中,點A,E,F在z方向都會存在不同程度的振蕩,并且距離中心點越近,其展開過程振蕩就越為顯著.由此可知,葉內(nèi)折疊在展開過程中,中心點發(fā)生振蕩的振幅很明顯,采用本文提出的桿–鏈動力學(xué)模型,能夠揭示折紙結(jié)構(gòu)的豐富的動力學(xué)行為.

圖13 葉內(nèi)折疊: 點A,E,F,G (圖10) z 坐標(biāo)時域圖Fig.13 Leaf-in origami: Time history of z coordinate of points A,E,F,G (see Fig.10)

2.3 基于Kresling 折疊的折紙結(jié)構(gòu)

Kresling 折疊是一類具有多穩(wěn)態(tài)特性的柱狀折紙,已廣泛應(yīng)用于工程中,如爬行機(jī)器人、機(jī)械臂、人造肌肉驅(qū)動器等,因此對其力學(xué)性能的探究很有現(xiàn)實意義.

圖14 給出了Kresling 平面折痕分布和柱狀結(jié)構(gòu),滿足以下幾何關(guān)系

圖14 Kresling 折疊: (a)折痕分布;柱狀結(jié)構(gòu): (b)主視圖,(c)俯視圖Fig.14 Kresling origami: (a) crease pattern.Columnar structure:(b) front view and (c) top view

其中,m表示柱狀結(jié)構(gòu)的胞元數(shù),n表示圓柱的邊數(shù),a和R分別為n邊形的邊長和外接圓半徑.在此基礎(chǔ)上,再給定單胞高度h和相鄰胞元的錯位角θ,則能夠確定Kresling 折疊的幾何結(jié)構(gòu).

2.3.1 不同多邊形Kresling 結(jié)構(gòu)的準(zhǔn)靜態(tài)分析

采用桿–鏈模型對非剛性Kresling 結(jié)構(gòu)進(jìn)行準(zhǔn)靜態(tài)分析,超彈性桿單元的材料參數(shù)分別為: 未變形時一維模量C0=5×107Pa,橫截面積A=1×10-5m2,選取Ogden 超彈性模型,N=2,α1=5,α2=1.非線性卷簧單元的參數(shù)為:k0=1 N·mm/rad,θ1=45o,θ2=315o.將Kresling 底部的各個節(jié)點與慣性基固結(jié),在頂層各節(jié)點施加豎直向下的軸向載荷,不計重力的影響.

為了分析多邊形邊數(shù)對Kresling 柱狀結(jié)構(gòu)力學(xué)特性的影響,分別取雙胞(m=2) 的正六邊形(n=6)、正八邊形(n=8)和正十邊形(n=10)進(jìn)行分析.已知外接圓半徑R=0.05 m,錯位角θ=45o,胞元高度h=0.05 m.本文采用改進(jìn)的廣義位移控制法[51]解準(zhǔn)靜態(tài)過程,能夠較為平滑地過渡極值點,并揭示結(jié)構(gòu)的多穩(wěn)態(tài)現(xiàn)象.

圖15 分別給出了在正六、八、十邊形Kresling折紙壓縮過程勢能和力隨位移變化的曲線.由圖15(a)可以看出,這三條勢能曲線都包含兩個阱,即state(0)和state(1),對應(yīng)的力–位移曲線(圖15(b))出現(xiàn)兩次正向穿越Fext=0 的水平軸,呈現(xiàn)出了Kresling 結(jié)構(gòu)的雙穩(wěn)態(tài)特征.可知,多邊形邊數(shù)n對Kresling 結(jié)構(gòu)的穩(wěn)定性特征很敏感.其中,三種折紙結(jié)構(gòu)的穩(wěn)態(tài)state(0)和state(1)對應(yīng)的構(gòu)型分別在圖15(a)中給出.

在圖15(b)中,A為未加載荷時的初始構(gòu)型.B1,B2,B3,D1,D2,D3為外力極大值點;C1,C2,C3為外力極小值點;E1,E2,E3為準(zhǔn)靜態(tài)仿真過程中所能到達(dá)的位移最大點.針對于上述特殊點,給出了對應(yīng)的Kresling 構(gòu)型圖,如圖16 所示.

圖16 Kresling 折疊: 不同多邊形中特殊點處(見圖15(b))構(gòu)型圖Fig.16 Kresling origami: Configuration of special points (see Fig.15(b)) in different shapes

通過對比圖15(b)中的三條力–位移曲線,可知在加載過程中(階段①和③),邊數(shù)與斜率呈正相關(guān)關(guān)系,說明隨著邊數(shù)的增加,Kresling 構(gòu)型的軸向剛度越大.觀察到正六邊形的曲線相較于正八和正十邊形呈現(xiàn)更復(fù)雜的非線性關(guān)系.綜上所述,不同多邊形改變了Kresling 的幾何設(shè)計,這些差異會對結(jié)構(gòu)的剛度有著顯著影響.增加邊數(shù)n,Kresling 結(jié)構(gòu)的剛度就越大,在不同穩(wěn)態(tài)構(gòu)型(三角形標(biāo)記)切換過程中所需的外力越大,因此可以通過改變邊數(shù)來調(diào)整結(jié)構(gòu)的穩(wěn)定性,可為折紙結(jié)構(gòu)的設(shè)計提供參考.

圖15 Kresling 折疊: 不同多邊形下,(a)勢能曲線,(b)力–位移曲線Fig.15 Kresling origami: (a) Potential energy curve and (b) forcedisplacement curve in different shapes

2.3.2 位移驅(qū)動下Kresling 結(jié)構(gòu)的動力學(xué)分析

在上一節(jié)中采用準(zhǔn)靜態(tài)模型研究了Kresling 的雙穩(wěn)態(tài)特性,并分析了不同邊數(shù)下柱狀Kresling 結(jié)構(gòu)的力學(xué)特性,接下來將采用非剛性桿–鏈模型對Kresling 結(jié)構(gòu)進(jìn)行動力學(xué)分析.如圖17 所示,本節(jié)選取雙胞正十邊形Kresling 折疊的折紙結(jié)構(gòu)作為研究對象,底層的十個節(jié)點與慣性基固結(jié),頂層的十個節(jié)點受到相同的位移驅(qū)動d(t),方向豎直向下.取系統(tǒng)的結(jié)構(gòu)阻尼系數(shù)為α=1.5.引入Kresling 胞元應(yīng)變的概念,則胞元e的應(yīng)變定義為 εe=(he-h0)/h0,其中he和h0分別表示胞元e當(dāng)前時刻和初始時刻的高度.非線性桿–鏈模型的參數(shù)選取與2.3.1 一致,采用廣義-α法求解動力學(xué)方程,計算時間te=4 s.

圖17 正十邊形Kresling 折紙結(jié)構(gòu)Fig.17 Kresling origami structure in regular decagon

取P1,P2,P3為觀測點,分別描述頂層、中間層和底層的動力學(xué)特性,其中h1和h2分別為上、下兩個胞元的高度.圖18(a)給出了z坐標(biāo)的時間歷程圖,點P1的曲線為施加在頂層各節(jié)點的位移驅(qū)動,即在0～1 s 為加速階段,加速度為0.03 m/s2,隨后以0.03 m/s 的速度勻速下壓,有zp1=h1+h2.此外,點P3為約束點,限制z=0.為了比較本文提出的卷簧本構(gòu)模型與傳統(tǒng)的卷簧本構(gòu)模型[42]的區(qū)別,圖18(a)給出了在兩種模型下點P2的運(yùn)動軌跡,可知在采用傳統(tǒng)的卷簧本構(gòu)模型時,會出現(xiàn)明顯的單元穿透問題,即h2出現(xiàn)跨單元情況.選取t=3.5 s 時刻下兩種模型在state C 和state D 處的構(gòu)型圖,可觀察到本文提出的卷簧本構(gòu)模型能夠保證Kresling結(jié)構(gòu)在壓縮過程中不會出現(xiàn)穿透.這說明了在接觸碰撞動力學(xué)問題中,改進(jìn)的卷簧本構(gòu)模型更具魯棒性.

通過觀察圖18(a)中點P2的運(yùn)動軌跡,發(fā)現(xiàn)不論是本文的卷簧本構(gòu)模型還是傳統(tǒng)的卷簧本構(gòu)模型,都在h2=0.023 m 處出現(xiàn)不穩(wěn)定的振蕩,對應(yīng)的state A 和state B 的構(gòu)型圖如圖18(b)所示.觀察曲線可知,本文模型的振蕩振幅會逐漸衰減直至達(dá)到完全壓縮狀態(tài),然而傳統(tǒng)的卷簧本構(gòu)模型會發(fā)生明顯的振蕩,并伴隨單元穿透問題.為了解釋這一振蕩現(xiàn)象,圖19 給出了Kresling 壓縮過程的應(yīng)變能云圖.觀察到應(yīng)變能云圖關(guān)于對稱軸h1=h2對稱.在位移驅(qū)動作用下,應(yīng)變能逐漸增大,Kresling 結(jié)構(gòu)的壓縮路徑趨向于最小應(yīng)變能的方向.在h2=0.045 m 處,最小應(yīng)變能路徑出現(xiàn)拐角,偏離對稱軸,此時壓縮主要體現(xiàn)在頂層單元壓縮.由于慣性的存在,點P2的壓縮運(yùn)動滯后于點P3,因此壓縮路徑在對稱軸上方.由2.3.1 節(jié)可知Kresling 結(jié)構(gòu)具有多穩(wěn)態(tài)特性,初始時刻應(yīng)變能最小,對應(yīng)穩(wěn)態(tài)state(0),另一個穩(wěn)態(tài)state(1)則存在于壓縮過程中.因此,在未到達(dá)state(1)前,應(yīng)變能減少,但外驅(qū)動力仍在做功,為滿足能量守恒定律,則動能增大出現(xiàn)振蕩.此外,在壓縮過程中,原路徑方向會遇到高應(yīng)變能區(qū)域的阻礙,導(dǎo)致路徑發(fā)生轉(zhuǎn)向,這也是發(fā)生振動的誘因之一.

圖18 正十邊形Kresling 折疊: (a) z 坐標(biāo)時域圖;(b) 不同卷簧本構(gòu)模型中特殊點處(見圖18(a))的構(gòu)型圖Fig.18 Regular decagon Kresling origami: Time histories of (a) z coordinate;(b) configuration of special points (see Fig.18(a)) in different rotational spring constitutive models

圖19 正十邊形Kresling 折疊: 應(yīng)變能云圖Fig.19 Regular decagon Kresling origami: Strain energy cloud map

為了進(jìn)一步分析壓縮過程的動力學(xué)特性,圖20給出了底層點P3處法向約束反力FN的時間歷程圖.由圖可知,Kresling 折紙結(jié)構(gòu)在未達(dá)到穩(wěn)態(tài)State(1)前,FN的變化曲線較為穩(wěn)定,但是隨著外驅(qū)動的不斷加載,臨近穩(wěn)態(tài)State(1),FN開始抖動,并且振幅與數(shù)值都在不斷增加.此外,在接近于完全壓縮時,會發(fā)生折疊面的接觸碰撞,從而導(dǎo)致FN數(shù)值振蕩的加劇.綜上所述,由動力學(xué)仿真的結(jié)果可知,Kresling 的雙穩(wěn)態(tài)直接影響了結(jié)構(gòu)的穩(wěn)定性.并且在接近于極限折疊(完全壓縮)時,由于接觸碰撞,其力學(xué)本構(gòu)呈現(xiàn)非光滑的特征,進(jìn)而伴隨出現(xiàn)復(fù)雜的非線性問題.

圖20 正十邊形Kresling 法向支座反力FN 時域圖Fig.20 Time history of the normal support reaction force FN of the regular decagon Kresling origami

2.3.3 正弦激勵下多鏈Kresling 折紙結(jié)構(gòu)的波動力學(xué)分析

在2.3.1 和2.3.2 節(jié)中,對雙胞Kresling 折紙結(jié)構(gòu)的準(zhǔn)靜態(tài)和動力學(xué)進(jìn)行了研究,揭示了其雙穩(wěn)態(tài)特性和復(fù)雜的動力學(xué)行為.在本節(jié)中,采用非剛性桿–鏈動力學(xué)模型,對一種由Kresling 胞元組合成的多鏈超材料折紙結(jié)構(gòu)進(jìn)行動力學(xué)分析.將該多鏈折紙結(jié)構(gòu)視為波傳播的媒介,討論其在正弦激勵下的波動力學(xué)特性.

如圖21 所示,選取30 胞元(m=30)的正六邊形(n=6)多鏈Kresling 結(jié)構(gòu)為研究對象,其中外接圓半徑R=0.03 m,錯位角θ=45o,胞元高度h0=0.03 m.在桿–鏈模型中,超彈性桿單元的材料參數(shù)與2.3.1 一致,非線性卷簧單元的參數(shù)為:k0=0.02 N·m/rad,θ1=5o,θ2=355o.將單層Kresling 結(jié)構(gòu)作為多鏈超材料折紙結(jié)構(gòu)單胞,胞元的應(yīng)變表達(dá)式見2.3.2 節(jié).取系統(tǒng)的結(jié)構(gòu)阻尼系數(shù)為α=1.5,忽略重力作用.將底層各節(jié)點與慣性基固結(jié),頂層各節(jié)點受到相同的正弦位移激勵d(t)=0.02sin(5πt) (m),方向沿著z負(fù)方向,激勵時長為0.05 s.采用廣義-α法求解系統(tǒng)的動力學(xué)方程,計算時間te=0.4 s.

圖21 多鏈Kresling 折紙結(jié)構(gòu)波動力學(xué)Fig.21 Wave dynamics of multi-chain Kresling origami structure

圖22 給出了不同時刻多鏈Kresling 折紙結(jié)構(gòu)的應(yīng)變–構(gòu)型圖,其中藍(lán)色代表壓縮,紅色代表拉伸.初始時刻應(yīng)變?yōu)?,右端受到擠壓正弦位移激勵作用.很明顯地觀察到,在t=0.15 s 之前,壓縮應(yīng)變逐漸從右端傳遞到左端,經(jīng)歷過主壓縮后的胞元,會出現(xiàn)出不同程度的拉伸應(yīng)變.因左端約束作用,在主壓縮應(yīng)變傳遞到單元30 后,會被吸收和反射后再次作用于多鏈中,可觀察到胞元的應(yīng)變呈現(xiàn)無規(guī)律的交替壓縮和拉伸.由上述現(xiàn)象可知,壓縮正弦激勵會產(chǎn)生應(yīng)變波,在多鏈Kresling 折紙結(jié)構(gòu)中傳播,并發(fā)生多重應(yīng)變波的疊加,導(dǎo)致胞元產(chǎn)生壓縮和拉伸的波動變化.

圖22 多鏈Kresling 折紙結(jié)構(gòu)的應(yīng)變–構(gòu)型圖Fig.22 Strain-configuration of multi-chain Kresling origami structure

為了更直觀地分析多鏈Kresling 結(jié)構(gòu)的波動力學(xué)行為,圖23 給出了應(yīng)變波傳播過程的時空變化圖.觀察到主壓縮波沿著單元數(shù)遞增的方向傳播,在到達(dá)單元30 后,主壓縮波會被反射,并且振幅會出現(xiàn)衰減.此外,在主壓縮波區(qū)域外,出現(xiàn)多條拉伸應(yīng)變波,其波峰數(shù)值小于壓縮波,并且在傳播過程中呈現(xiàn)不連續(xù)現(xiàn)象.隨著壓縮和拉伸波在多鏈Kresling 結(jié)構(gòu)中的傳播和疊加,結(jié)構(gòu)的波動力學(xué)行為更為復(fù)雜,揭示了折紙超材料中波傳播的特性.

圖23 多鏈Kresling 折紙結(jié)構(gòu)應(yīng)變波傳播的時空圖Fig.23 Space-time of strain wave propagation in multi-chain Kresling origami structure

圖24 描繪了多鏈Kresling 折紙結(jié)構(gòu)中不同單元的應(yīng)變時間歷程圖.初始時刻,在正弦壓縮沖擊下產(chǎn)生主壓縮應(yīng)變波,隨著時間推移,逐漸向后續(xù)的折紙胞元傳播.從圖中觀察到,胞元在經(jīng)歷主壓縮波后,會出現(xiàn)拉伸與壓縮的振蕩,并且拉伸的幅值越靠近約束段越大,伴隨多重應(yīng)變波的疊加,應(yīng)變曲線呈現(xiàn)出復(fù)雜的非線性特性.

圖24 多鏈Kresling 折紙結(jié)構(gòu): 不同單元的應(yīng)變時間歷程圖Fig.24 Multi-chain Kresling origami structure: Strain versus time curves in different unites

為了進(jìn)一步對多鏈Kresling 折紙結(jié)構(gòu)吸能特性進(jìn)行研究,將仿真時間延長至1 s,圖25 給出了能量–時間曲線圖.在t=0.05 s 時,位移激勵釋放,外力將不再做功,系統(tǒng)總能量保持不變.很明顯地觀察到,結(jié)構(gòu)的彈性勢能和動能的振幅逐漸減少,耗散能逐漸增大.一般地,耗散能主要包括結(jié)構(gòu)阻尼耗散能、摩擦耗散能、黏性耗散能等.觀察彈性勢能和動能的曲線,發(fā)現(xiàn)這兩條曲線的振蕩頻率與振幅大小衰減的基本一致,并且相位相反.由結(jié)果可知,在位移激勵釋放后,系統(tǒng)能量耗散較為明顯,說明了Kresling折紙結(jié)構(gòu)具有較好的吸能和減震的功能.

圖25 多鏈Kresling 折紙結(jié)構(gòu)的能量–時間曲線Fig.25 Energy-time curves of multi-chain Kresling origami structure

綜上所述,采用本文提出的桿–鏈模型,能夠揭示非剛性多鏈Kresling 折紙結(jié)構(gòu)的波動力學(xué)特性,描述應(yīng)變波在折紙結(jié)構(gòu)中的傳播過程,并揭示波的傳播與疊加規(guī)律.通過能量分析,體現(xiàn)了Kresling折紙結(jié)構(gòu)吸能和減震的特點,可為沖擊緩沖和能量收集的折紙設(shè)計提供理論依據(jù).

3 結(jié)論

隨著材料科學(xué)技術(shù)的不斷發(fā)展,非剛性折紙結(jié)構(gòu)在工程中的應(yīng)用越來越廣泛.本文把非剛性折紙結(jié)構(gòu)簡化為帶卷簧的空間桁架等效結(jié)構(gòu),建立了一種通用且可處理的非線性桿–鏈動力學(xué)模型.基于Ogden 超彈性本構(gòu)關(guān)系,推導(dǎo)了桿單元的非線性彈性力陣和切線剛度陣,并將折痕和虛擬折痕等效為桿單元,可處理折紙結(jié)構(gòu)的大變形和大范圍運(yùn)動問題.接下來,采用具有非線性旋轉(zhuǎn)剛度的卷簧來模擬折痕的抗彎作用和虛擬折痕的彎曲作用.為了有效解決接觸碰撞動力學(xué)中折疊面的穿透問題,改進(jìn)了卷簧的本構(gòu)模型,相較于傳統(tǒng)的卷簧模型,具有更強(qiáng)的通用性和魯棒性.在此基礎(chǔ)上,推導(dǎo)了廣義質(zhì)量陣,并考慮阻尼效應(yīng),用虛功原理建立了非剛性折紙結(jié)構(gòu)多體系統(tǒng)的動力學(xué)方程.

本文將提出的桿–鏈動力學(xué)模型應(yīng)用于三種經(jīng)典的非剛性折紙結(jié)構(gòu)中.第一種單菱形折紙結(jié)構(gòu),通過比較準(zhǔn)靜態(tài)與解析解的結(jié)果,驗證了桿–鏈模型的準(zhǔn)確性,隨后與基于絕對節(jié)點坐標(biāo)法的動力學(xué)仿真結(jié)果進(jìn)行對比,說明了桿–鏈動力學(xué)模型的準(zhǔn)確性與高效性.第二種是基于葉內(nèi)折疊的折紙結(jié)構(gòu),首先采用添加虛擬折痕和修正初始構(gòu)型的方法,有效解決了剛性折紙模型中展開和收攏過程的鎖定問題,對剛性折紙的可折疊性和奇異問題進(jìn)行討論.然后,將非剛性與剛性折紙模型進(jìn)行對比,發(fā)現(xiàn)非剛性折紙模型不會發(fā)生自鎖現(xiàn)象,體現(xiàn)了柔性變形對折紙結(jié)構(gòu)大范圍運(yùn)動的影響,此外還能夠給出具有大變形的張緊構(gòu)型,更具通用性.第三種是Kresling 折紙結(jié)構(gòu),首先對不同邊數(shù)下Kresling 結(jié)構(gòu)進(jìn)行準(zhǔn)靜態(tài)分析,發(fā)現(xiàn)了Kresling 的雙穩(wěn)態(tài)特性.接下來對位移驅(qū)動下的Kresling 折紙結(jié)構(gòu)進(jìn)行動力學(xué)分析,揭示了在多穩(wěn)態(tài)影響下誘發(fā)的非線性動力學(xué)行為.最后,基于Kresling 單胞發(fā)展一種多鏈超材料折紙結(jié)構(gòu),并展現(xiàn)了在正弦激勵作用下結(jié)構(gòu)的波動力學(xué)特性,揭示應(yīng)變波的傳播和疊加現(xiàn)象,通過能量分析,體現(xiàn)了Kresling 折紙結(jié)構(gòu)吸能和減震的特點.

本文提出的桿–鏈動力學(xué)模型,有效簡化了非剛性折紙結(jié)構(gòu)的動力學(xué)建模過程,能夠較為準(zhǔn)確和高效地對非剛性折紙結(jié)構(gòu)的動力學(xué)進(jìn)行數(shù)值模擬,揭示豐富的非線性動力學(xué)特性.未來可將該桿–鏈動力學(xué)模型拓展到非剛性折紙結(jié)構(gòu)的優(yōu)化設(shè)計和動力學(xué)控制中.