考慮非線性本構的非剛性折紙結構動力學建模與分析1)

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

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

引言

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

剛性折紙作為折紙的一個特殊子集,其數學模型的建立通常基于以下假設: 折痕是直線折痕,折疊面視為剛性,即不會發生拉伸或彎曲[25].目前,對剛性折紙的研究已經較為成熟,并已開發了多個軟件[26-27].隨著材料科學技術的不斷發展,折紙結構在工程中的應用更為廣泛,在很多情況下,折紙結構的折疊面會發生拉伸、彎曲或剪切等變形,甚至在展開過程中還會誘發接觸碰撞等問題.因此,針對這類非剛性折紙結構,用剛性折紙模型進行模擬將無法反映非剛性折紙特有的力學性能.

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

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

1 非剛性折紙結構動力學建模

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

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

1.1 超彈性本構的桿單元

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

如圖1 所示,桿單元e的廣義坐標陣為

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

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

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

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

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

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

定義一維線性模量為

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

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

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

考慮到材料未變形時(λ1=1)為無應變無應力狀態,有如下約束關系

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

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

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

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

1.2 非線性卷簧單元

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

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

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

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

式中sgn 為符號函數.

將卷簧單元變形能對卷簧單元廣義坐標陣進行變分,可以得到卷簧單元的內力陣為

對卷簧單元的內力陣求變分,可得卷簧單元切向剛度陣

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

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

考慮到卷簧力矩的連續性,本文改進了非線性卷簧的本構模型,相較于傳統的卷簧本構模型[42],該改進的本構模型在臨近接觸時延長了過渡區,如圖3 所示.采用改進的本構模型求解動力學接觸問題時,能夠有效避免因時間積分步過大導致的穿透問題,則改進的本構模型中卷簧彎矩的具體表達式如下

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

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

1.3 非剛性折紙結構動力學方程

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

其中,系統的結構阻尼陣為C=αM,Φq=?Φ/?q為約束雅克比矩陣.

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

2 數值算例仿真

將本文提出的非線性桿–鏈模型應用于三種經典的非剛性折紙結構的動力學分析中.第一個算例是進行單菱形折紙結構的動力學仿真,將桿–鉸鏈模型分別與解析解、準靜態以及有限元絕對節點坐標法ANCF[31]進行對比,驗證桿-鉸鏈模型的準確性和高效性.第二個算例是對葉內折疊折紙結構進行動力學仿真,并將剛性模型和非剛性模型的動力學結果進行對比分析.第三個算例是對Kresling 折紙結構進行仿真,揭示了折紙結構的多穩態和動力學特性.

2.1 單菱形折疊

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

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

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

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

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

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

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

2.2 基于葉內折疊的折紙結構

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

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

圖6 葉內折疊的折痕分布Fig.6 Crease pattern of leaf-in origami

2.2.1 剛性折紙結構動力學分析

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

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

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

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

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

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

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

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

針對剛性不可展問題,本文給出了添加虛擬折痕和修正初始構型的方法,有效解決了展開和收攏過程的鎖定問題.本文對剛性折紙結構的可展性和動力學研究,為工程問題中折紙結構的運動控制和設計優化提供理論指導.

2.2.2 非剛性折紙結構動力學分析

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

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

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

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

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

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

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

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

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

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

2.3 基于Kresling 折疊的折紙結構

Kresling 折疊是一類具有多穩態特性的柱狀折紙,已廣泛應用于工程中,如爬行機器人、機械臂、人造肌肉驅動器等,因此對其力學性能的探究很有現實意義.

圖14 給出了Kresling 平面折痕分布和柱狀結構,滿足以下幾何關系

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

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

2.3.1 不同多邊形Kresling 結構的準靜態分析

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

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

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

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

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

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

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

2.3.2 位移驅動下Kresling 結構的動力學分析

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

圖17 正十邊形Kresling 折紙結構Fig.17 Kresling origami structure in regular decagon

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

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

圖18 正十邊形Kresling 折疊: (a) z 坐標時域圖;(b) 不同卷簧本構模型中特殊點處(見圖18(a))的構型圖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 折疊: 應變能云圖Fig.19 Regular decagon Kresling origami: Strain energy cloud map

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

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

2.3.3 正弦激勵下多鏈Kresling 折紙結構的波動力學分析

在2.3.1 和2.3.2 節中,對雙胞Kresling 折紙結構的準靜態和動力學進行了研究,揭示了其雙穩態特性和復雜的動力學行為.在本節中,采用非剛性桿–鏈動力學模型,對一種由Kresling 胞元組合成的多鏈超材料折紙結構進行動力學分析.將該多鏈折紙結構視為波傳播的媒介,討論其在正弦激勵下的波動力學特性.

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

圖21 多鏈Kresling 折紙結構波動力學Fig.21 Wave dynamics of multi-chain Kresling origami structure

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

圖22 多鏈Kresling 折紙結構的應變–構型圖Fig.22 Strain-configuration of multi-chain Kresling origami structure

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

圖23 多鏈Kresling 折紙結構應變波傳播的時空圖Fig.23 Space-time of strain wave propagation in multi-chain Kresling origami structure

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

圖24 多鏈Kresling 折紙結構: 不同單元的應變時間歷程圖Fig.24 Multi-chain Kresling origami structure: Strain versus time curves in different unites

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

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

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

3 結論

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

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

本文提出的桿–鏈動力學模型,有效簡化了非剛性折紙結構的動力學建模過程,能夠較為準確和高效地對非剛性折紙結構的動力學進行數值模擬,揭示豐富的非線性動力學特性.未來可將該桿–鏈動力學模型拓展到非剛性折紙結構的優化設計和動力學控制中.