增強六手臂缺失支柱手性拉脹超材料力學性能理論研究1)

朱一林 江松輝 于 超

* (西南石油大學土木工程與測繪學院,工程安全與防護研究院,成都 610500)

? (西南交通大學力學與航空航天學院,應用力學與結構安全四川省重點實驗室,成都 610031)

引言

拉脹材料是自20 世紀80 年代起迅速發展起來的一類具有負泊松比效應的功能和結構一體化的力學超材料[1-3].與傳統材料相反,拉脹超材料承受拉伸(壓縮)載荷時,會產生橫向膨脹(收縮).這種反常規的變形行為賦予了拉脹超材料優越的力學性能,表現在抗剪切、抗壓縮、抗沖擊、抗斷裂、抗疲勞,能量吸收和減震等諸多方面[4-8].優越的力學性能使拉脹超材料在航空航天[9-11]、軍事國防[12-15]、汽車船舶[1,16]、生物醫學[17-19],土木工程[5,20-22]及智能傳感[23-26]等領域有廣闊的應用前景,并吸引了國內外廣大學者的廣泛關注.

負泊松比是一個古老的研究課題.早在1848 年,法國科學家圣維南[27]就指出在各向異性材料中可能會有負泊松比存在.經過近一個半世紀的發展,尤其是自Lakes 于1987 年首次制備出負泊松比材料以來,國內外力學和材料學界迎來了對負泊松比材料的研究高潮.1991 年,Evans[28]將具有負泊松比特性的材料正式命名為“auxetics”(源自于希臘文“auxetikos”).“auxetics”后被國內學者譯為“拉脹材料”.由于可呈現出天然材料所不常有的力學特性,拉脹材料又被稱作“拉脹超材料”.

根據結構形式和變形機理的不同,拉脹超材料大致可歸類為纖維增強層壓型拉脹復合材料,拉脹紗,拉脹泡沫,拉脹折紙和拉脹蜂窩材料五種形式[6,26,29-32].拉脹蜂窩材料具有規則的元胞結構,非常便于設計、優化和制備,因此近些年得到了更廣泛的關注.根據微結構的幾何形式和主導的變形機理,拉脹蜂窩材料又可以進一步劃分為內凹型拉脹蜂窩材料(具有內凹角,具有軸對稱性質)和手性拉脹蜂窩材料(具有旋轉對稱性質,變形時伴隨著旋轉)[33-34].相較于內凹型拉脹材料,手性拉脹材料在大變形下的拉脹性能表現更好,并且可容許更大的制造誤差[35-38].傳統手性拉脹材料的元胞結構由中心圓環和與之相切外接組成.經典的手性拉脹材料包括傳統手性拉脹材料和缺失支柱手性拉脹材料兩種材料形式.根據與單個圓環連結的手臂個數及手臂與相鄰兩圓環的連結形式,傳統手性拉脹材料包括三手臂手性、三手臂反手性、四手臂手性、四手臂反手性及六手臂手性拉脹材料[35,39-41](圖1).缺失支柱手性拉脹材料可視為將傳統手性拉脹材料的圓環替換為綁扎桁架而構成[34,42-48](圖2).其余手性拉脹材料均可視為由這兩種材料演變而成.

圖1 傳統手性拉脹材料: (a) 三手臂;(b) 反三手臂;(c) 四手臂;(d) 反四手臂;(e) 六手臂Fig.1 (a) Tri-,(b) anti-tri-,(c) tetra-,(d) anti-tetra-,(e) hexa-chiral honeycombs

圖2 缺失支柱手性拉脹材料: (a) 三手臂;(b) 反三手臂;(c) 四手臂;(d) 反四手臂;(e) 六手臂Fig.2 (a) Tri-,(b) anti-tri-,(c) tetra-,(d) anti-tetra-,(e) hexa-missing rib honeycombs

傳統手性拉脹材料的拉脹效應主要源自其中心圓環的轉動和手臂的彎曲.然而,已有研究表明在大變形時其中心圓環容易發生畸變,旋轉效應相應減弱,進而導致拉脹性能的迅速劣化[35].因此,傳統手性拉脹材料無法在大變形范圍內具備恒定的負泊松比,為其實際的工程應用帶來了困難.根據Simth 等[48]、Liu 等[43]以及作者團隊[3,46-47]的研究表明,缺失支柱拉脹材料可以在較大的應變范圍內表現出恒定的泊松比.與傳統手性拉脹材料類似,缺失支柱手性拉脹材料的拉脹效應也主要源自其中心綁扎桁架的轉動和手臂的彎曲.但由于組成中心綁扎桁架的梁之間缺乏支撐,初始的缺失支柱手性拉脹材料的旋轉效應較弱,因而其拉脹性能也較弱.作者團隊[47]在具有初始面內各向同性的六手臂缺失支柱手性拉脹材料的中心綁扎桁架中引入了簡單的六邊形結構支撐,發展了增強六手臂(直臂)缺失支柱手性拉脹材料(圖3(a)).該拉脹材料可以在大應變范圍內具有可調的恒定負泊松比.特別是,通過調整幾何參數,該拉脹材料可在近30% 的應變范圍內具有趨近于-1 的負泊松比.眾所周知,各向同性結構材料的等效剪切模量為和分別為等效彈性模量和泊松比).因此,結構材料的負泊松比趨于-1 時,其剪切模量理論上趨于無窮大,具有絕佳的抗剪性能.從而,增強六手臂(直臂)缺失支柱手性拉脹材料具有廣闊的工程應用前景.為進一步增大具有恒定負泊松比的應變范圍,作者團隊還進一步發展了增強六手臂(曲臂)缺失支柱手性拉脹材料(圖3(b)).

圖3 增強六手臂缺失支柱手性拉脹材料元胞圖Fig.3 Unit-cells of enhanced hexa-missing rib honeycomb

但作者團隊的前期工作僅局限于有限元分析,缺乏對所發展的拉脹材料力學參數的理論描述,無法根據目標性能進行有效的結構參數設計,限制了其工程應用.而參數化設計有助于對超材料進行結構優化,從而促進其工程應用[49-50].本文所研究的拉脹蜂窩材料大都由細長結構組成,可被假設為歐拉-伯努利梁[51].大量研究證實利用能量法(如卡氏定理)可有效地推導拉脹蜂窩材料力學參數的理論表達式.因此,本文致力于利用能量法給出前期發展的增強六手臂缺失支柱手性拉脹材料等效彈性模量和泊松比的理論表達式,以期為拉脹材料的有效參數設計提供理論指導.

1 幾何模型

增強六手臂缺失支柱手性拉脹材料為典型二維結構材料.其元胞如圖3 所示,由“Z”形梁結構及增強六邊形組成.元胞還可視為由六個最小代表單元(紅色虛線部分)陣列而成.直臂結構的元胞幾何形狀可由五個獨立參數確定,即“Z”形梁長手臂長度L;Z”形梁短手臂長度r;“Z”形梁內角 φ ;增強六邊形結構邊長l以及所有梁結構的寬度t(如圖3(a)所示).由圖3 可見,元胞結構的綁扎桁架由“Z”形梁短手臂及增強六邊形組成,而“Z”形梁長手臂則與其相接.因此,為后續敘述方便,下文中“Z”形梁長手臂和短手臂分別記為外手臂和內手臂.

將圖3(a)中的直手臂替換為等長圓弧即可得相應曲臂拉脹材料.因此,要描述曲臂材料的元胞幾何形狀,還需引入表征其半徑的參數R.需要指出的是前述所有表征梁長度和半徑的參數均沿梁中性軸測量.

為后續描述方便,本文定義了參數b來表征增強六邊形與手臂之間的距離.b可由r和l表示,即

另外,本文還引入了增強率p及手臂比q分別來描述材料增強的程度及內外手臂長度的比例.p和q的表達式為

2 等效力學參數理論表達式(小變形)

如前所述,本文所研究的兩種拉脹材料均可視為由多個歐拉-伯努利梁組成.因而,可利用能量法推導其等效彈性模量和泊松比.利用能量法推導等效力學參數通常的做法是: 首先,對陣列的拉脹材料施加遠場應力;然后,取最小代表性單元,并根據幾何和受力條件確定各梁的受力狀態;最后,推導最小代表性單元的應變能,并根據卡氏定理和幾何相容條件獲得等效彈性模量和泊松比的表達式.

對于具有簡單最小代表單元構型的拉脹材料(如最小代表單元具有三角形邊界的三手臂傳統/缺失支柱反手性拉脹材料,最小代表單元具有矩形邊界的四手臂傳統/缺失支柱反手性拉脹材料等),最小代表單元的受力狀態可由幾何和遠場應力條件確定[3,37,46].然而,對于本文所研究的拉脹材料,其最小代表單元的邊界較為復雜(圖3(a)),無法通過傳統方法推導其力學參數的理論表達式.Liu 等[43]發展了基于最小代表單元整體受力平衡條件和變形協調條件的方法,準確地推導了具有弱連結的六手臂缺失支柱手性拉脹材料的力學參數的理論表達式.六手臂缺失支柱手性拉脹材料與Liu 等[43]所提出的材料具有相似的構型,因此,本文將基于同樣的思路來進行理論推導.

以增強六手臂(直臂)缺失支柱手性拉脹材料為例.對其施加如圖4(a)所示遠場拉應力,其最小代表單元的受力情況如圖4(b)所示.該最小代表單元由三個增強桿和三根“Z”型梁組成.前期研究表明[3],由增強桿和“Z”型梁組成三角形單元體在結構整體變形時幾乎不會發生畸變,局部應變非常小,可以視為剛體.因此,推導最小代表單元的計算應變能時,可忽略三角形元素的貢獻.最小代表單元可進一步簡化為只由三根 “Z”型梁相互連結而成(圖4(c)),并且標注為紅色的部分(對應于三角形單元體)假定為剛體,不對應變能做貢獻.另外,如圖4(c)所示的等效最小代表單元的各梁段為連續梁,剛體和變形體連接處滿足撓度和轉角連續性條件.

圖4 (a) 陣列結構受遠程應力示意圖;(b) 最小代表單元受力圖;(c) 等效最小代表單元Fig.4 (a) Schematic diagram of an enhanced hexa-missing rib lattice under far field uniaxial loading.(b) The related free boundary diagram of the smallest repeated unit.(c) The effective smallest repeated unit

2.1 任意形狀的歐拉-伯努利梁力學模型

要確定最小代表單元的整體受力和變形條件,首要分析其“Z”形梁組件的力-位移關系.圖4(c)所示的“Z”形梁變形時,可簡化為一端固定鉸接并作用有集中力偶MA,另一端滑動鉸接,并作用有集中力偶MB和水平集中力N0的簡支梁(如圖5(a) 所示).去掉支座,代之以力,其對應的受力情況如圖5(b)所示.由圖5(b)易知

為了更具有普適性,“Z”形梁上任一截面(以軸線長度s標記,如圖5(c)所示)的坐標可構造為軸線長度s的函數,即

圖5 任意歐拉-伯努利梁結構受力圖Fig.5 Force diagram of an Euler-Bernoulli beam with arbitrary shape

則截面s上的軸力N和彎矩M與x和y方向上的內力N0和Q0有關,分別為

其中,ES為基體材料的彈性模量;ES A=ES t和ES I=ES t3/12分別為單位厚度梁的拉伸剛度和彎曲剛度.

聯合式(2)～式(5),可得

式中,αi(i=1,2,···,9) 為一組只取決于結構形狀的無量綱參數,表達式為

運用卡氏定理,可以得出梁滑動端的位移(記為u,沿x方向)及兩端旋轉角 ωA和 ωB的表達式為

對式(8) 取逆,并代入式(9) 所示的歸一化關系,可得

其中,C為柔度矩陣,表達式見附錄A.其逆形式為

其中,D=C-1為剛度矩陣.

2.2 增強六手臂缺失支柱拉脹材料力學模型

如前所述,增強六手臂缺失支柱拉脹材料的最小代表單元可簡化為如圖4(c)所示的結構.打開三根“Z”形梁之間的連結,受力情況如圖6 所示.

圖6 最小代表單元各“Z”型梁受力圖Fig.6 Force diagram of the zigzag ligaments of the smallest repeated unit

為描述方便,按照逆時針將“Z”形梁分別命名為1,2,3,并令Ni,Mi,Qi(i=1,2,3)表示結構的軸力、彎矩和剪力.

聯立式(2)和式(9)可得

進一步根據梁2,3 連結處及梁1,2 連結處的受力平衡,可得

再考慮圖4(b)所示的受力情況,根據水平方向受力平衡,可得

如前所述,每個元胞包含六個最小代表單元(圖3(a)).

取元胞中心點附近的微小分離體,其受力情況如圖7所示.根據力矩平衡條件,可得

血管環指主動脈弓及其分支血管的先天發育異常，由于胚胎時期多對鰓弓及成對的背側主動脈未能順序融合和吸收，在解剖上形成完全性或不完全性血管環，圍繞氣管和（或）食管造成壓迫的一組先天性血管畸形，占先天性心臟病的1%～2%。每類血管環的預后存在顯著差異，是否早期正確診斷和及時手術治療是新生兒能否存活的關鍵。本研究是經過三血管序列切面對主動脈弓、動脈導管弓和氣管三者關系的快速篩查，發現異常，再經弓降部管狀切面進行驗證，通過分析常見的血管環的胎兒超聲心動圖不同切面的特點，提高血管環的發現率及準確率。

圖7 元胞中心節點區受力圖Fig.7 Force diagram of an area surrounding the central point of the unit-cell

將式(12)代入式(15),進一步可得

本文僅考慮線彈性小變形情況,則最小代表單元三個“Z”形梁組件的三個連結點形成的圖形在變形前后適中保持為封閉的三角形(如圖8 所示).根據幾何條件和變形相容條件,可得

圖8 最小代表單元變形示意圖Fig.8 Schematic diagram of the deformation of smallest repeated unit

其中,θA,θB為桿件變形前兩端切線與x軸的夾角(見圖5); γiA,γiB為變形后夾角; φ1,φ2和 φ3為三角形變形后內角,L1,L2和L3為變形后梁1,2 和3 端點之間的直線長度.

由式(17b)進一步可得

兩邊取微分,有

其中

小變形時可認為

聯立式(17a)、式(20)和式(21)可得

式中

同理可得

考慮最小代表單元的鏡像結構(如圖8(b)所示),易知

聯立式(17a)和式(24)可知

圖8(a)所示的最小代表單元在水平方向和豎直方向變形前后的長度分別為

變形后的協調條件為,梁“2”和“3”在水平方向的投影長度之和與梁“1”的長度相等,即

對式(28c)、式(28d)和式(29)微分可得

由式(30b)和式(31)可知

則最小代表單元的水平平均應變為

聯立式(21)、式(24)和式(30a)可得

聯立式(17a)、式(20)及式(27)可知

進一步聯立式(26)、式(27)、式(34)和式(35)可知

則最小代表單元的豎直平均應變可進一步表示為

材料的等效彈性模量(以基體材料彈性模量進行歸一化)和泊松比分別為

借助Maple 軟件,聯立方程式(11)～式(14)、式(17)、式(23)、式(25)～式(27)、式(33)及式(37)～式(39)可得等效彈性模量和泊松比分別為

針對特定的結構形式,給出“Z”形梁的曲線描述方程,進一步給出式(7)中無量綱參數的具體形式,并帶入式(40)和式(41)進行化簡,即可得材料最終的等效彈性模量和泊松比的理論表達式.

需要指出的是,本節所描述的六手臂缺失支柱型手性拉脹材料等效力學參數的理論推導過程與Liu 等[43]的思路一致.但他們的工作中對如何利用平衡條件得到方程式(12)～式(14),及利用變形協調條件得到方程式(33)和式(37)闡述并不清晰.本文詳細地陳述了這兩個推導過程,增加了文章的可讀性.

增強六手臂缺失支柱手性拉脹材料和增強六手臂(曲臂)缺失支柱手性拉脹材料的代表性“Z”形梁的計算模型如圖9 所示.兩種材料的曲線描述方程可見附錄A.

圖9 (a) 直臂與 (b) 曲臂結構計算模型Fig.9 Calculation diagram of the enhanced hexa-missing rib auxetic honeycombs with (a) straight ligament and (b) wavy ligament

需要指出的是,通常情況下兩種結構彈性模量和泊松比的理論表達式為包含上千項的多項式,無法繼續化簡.但根據理論推導,作者已經基于MATLAB 軟件創建了圖形用戶界面(GUI),界面如圖10 所示.對于直臂結構,只需輸入內手臂長度r,外手臂長度L,增強率p,內角φ,手臂寬度t五個獨立參數即可獲得特定微結構形式的等效彈性模量和泊松比.對于曲臂結構,還需輸入外手臂的半徑R.

圖10 圖形用戶界面: (a) 直臂結構計算和(b) 曲臂結構計算Fig.10 Graphical user interface: (a) calculation of auxetics with straight ligaments and (b) calculation of auxetics with wavy ligaments

另外,研究還發現,在L=2r這種特定的情況下,直臂結構的等效彈性模量和泊松比具有簡潔的理論表達式,即

其中k=tanθ .

3 等效力學參數理論解討論及有限元結果對比

本節將系統討論直臂結構的等效彈性模量和泊松比與幾何參數之間的關系,并將相應結果與有限元模擬進行對比.

為簡化分析,本節所討論的拉脹材料的手臂長度和手臂寬度之間的比例關系均固定為L:t=20.另外,作者的前期研究表明: 幾何參數相同時,曲臂缺失支柱拉脹材料的負泊松比響應與直臂材料接近(圖10 也可以證明),僅在增加保持恒定負泊松比的應變范圍方面更有優勢.因此,本節將僅討論直臂結構.

3.1 力學參數云圖

為了直觀地研究幾何參數對增強六手臂缺失支柱手性拉脹材料宏觀力學性能的影響,圖11 和12 分別給出了等效泊松比和彈性模量隨手臂和中心綁扎桁架夾角 φ,增強率p及手臂比q變化的云圖.

由圖11 和圖12 可見,改進的拉脹結構可以提供在大范圍可調的泊松比(-1 到0)和彈性模量.由圖11(a)可知: (1)結構的拉脹效應隨著“Z”型手臂內角的增大而減弱,這是由于內角的增大意味著長手臂和中心綁扎桁架的連結處剛度減弱,抑制了中心綁扎桁架的旋轉效應;(2)由于中心綁扎桁架的剛性隨增強率的增加而增強,結構的拉脹效應隨著增強率的增大而增強.但由圖8(b)可知,拉脹效應與手臂比的關系是非單調的,表現出兩個階段的演變,在第一階段(q達到閾值之前,圖11(b)中黑色線表示),手臂比越大拉脹效應越強,第二階段則相反.這是因為q越大意味著外手臂相對中心綁扎桁架越強,因而更容易觸發結構的旋轉效應,具有更好的拉脹效應.但q持續增大時,結構逐漸趨于三角形網絡結構,不具有拉脹特性.因此,結構的泊松比演化和q的關系表現出了兩階段的特征.

圖11 結構等效泊松比云圖Fig.11 Contour plots the effective Poisson’s ratio over a wide range

由圖12 可知,結構的彈性模量隨著夾角、增強率和手臂比的增大而單調增大.

圖12 結構等效彈性模量云圖Fig.12 Contour plots the effective elastic modulus over a wide range

3.2 有限元結果對比

本節進一步將等效泊松比和彈性模量的理論解析結果與有限元計算結果進行了比較.

有限元計算采用ABAQUS 軟件完成.計算時設定如下: (1)基體材料彈性模量和泊松比分別為1 和0.3;(2) 單元類型為CPE4 四節點平面應變單元;(3)根據網格收斂性分析,網格尺寸為0.25t,即: 保證沿梁厚度方向至少有4 個單元;(4)由于第2 章的理論推導是在線彈性小變形下開展的,本節的有限元計算也局限為線彈性小變形分析;(5)計算針對單一元胞結構展開,為消除邊界條件的影響,對元胞結構施加了周期性邊界條件(詳見附錄B).

圖13 給出了變形較大時(等效應變為5%)不同幾何參數最小代表單元的最大面內主應變云圖.由圖可見,結構整體變形較大時,各部件的局部變形相對很小,驗證了假定(4)的合理性.另外,圖13 還表明,結構的三角形元素部分的應變很小,也說明了理論推導時將其假定為剛體部件的合理性.

圖13 等效應變為5%時最小代表單元最大面內主應變云圖Fig.13 Contour plots the maximum in-plane principal strain of the smallest repeated units at an effective strain of 5%

圖14 和圖15 為不同幾何尺寸下,結構等效泊松比和彈性模量的理論值與有限元結果.顯然,兩種結果吻合的較好.

圖15 不同幾何參數下結構等效彈性模量的理論和有限元模擬結果對比圖Fig.15 Theoretical and FE results of the effective elastic modulus of the enhanced hexa-missing rib auxetic honeycomb over a range

需要指出的是,本文的理論解是基于將最小代表單元的三角形部分假定為剛體這一假設而推導的.在結構受載時,三角形部分不可避免地也會產生變形(并且其變形的程度受結構相對尺寸的影響),這部分變形同樣也會對結構總變形能做貢獻.另外,理論解中將實體結構簡化為梁,所有與長度相關的尺寸均為各梁段中性軸的長度.但事實上實體結構的結果與梁的結果也會有偏差.理論模型的精度及誤差主要源自上述兩個因素,并且這兩個因素對理論值的影響并不同向,會存在競爭關系.因此,圖14和15 中部分理論結果比有限元結果大,而另一部分比有限元結果小.

圖14 不同幾何參數下結構等效泊松比的理論和有限元模擬結果對比圖Fig.14 Theoretical and FE results of the effective Poisson’s ratio of the enhanced hexa-missing rib auxetic honeycomb over a range

4 結論

基于能量法,本文首先在小變形框架下推導了前期研究工作中發展的增強六手臂缺失支柱手性拉脹超材料的等效力學性能參數的理論模型.研究表明,只有在微結構具有特定幾何參數時,增強六手臂缺失支柱手性拉脹超材料的等效泊松比和彈性模量才有簡潔的理論表達式.因此,為便于理論模型的應用,本文進而開發了MATLAB 圖形用戶界面,只需輸入幾何參數就可以方便地直接獲取相應結構的等效力學性能參數.最后,系統討論了關鍵幾何參數對結果力學性能的影響規律.

需要指出的是,本文的研究局限于理論推導和有限元分析.在后續的工作中,作者擬開展系統的實驗研究來考察增強六手臂缺失支柱手性拉脹超材料的靜/動態力學性能.

附錄A

其中

直臂結構曲線描述方程

曲臂結構曲線描述方程

附錄B

為簡化計算,本文的分析只針對單一元胞結構,采用基于ABAQUS 軟件的有限元均勻化方法來考察手性拉脹材料的等效力學行為.為消除邊界條件的影響,計算時對元胞結構施加了周期性邊界條件[52-53],即

式中,x±和X±分別為對應邊界的空間和參考坐標;t±=P±N±代表外法向為N±邊界上節點的牽拉應力(P為第一類皮奧拉-基爾霍夫(Piola-Kirchhoff)應力).等效變形梯度Fˉ 及等效第一類皮奧拉-基爾霍夫應力可表示為

式(B2)和式(B3)中,V0和 ?V0分別為元胞在參考坐標系下的體積和邊界.

元胞的整體變形可由等效Biot 應變[54]來度量,即

其中,1 為二階單位張量.

等效泊松比可定義為

圖B1 給出了特定結構尺寸下(p=0.6,q=1,φ=90°)對一個單胞施加周期性邊界條件以及對陣列結構(含n×n個單胞)直接建模所得的等效泊松比結果的對比.由圖可見,隨著n(n為陣列結構中完整單胞的數目)的增加,陣列結構直接建模的計算結果將收斂于周期性邊界條件獲得的結果,驗證了周期性邊界條件結果的準確性.

圖B1 對一個單胞施加周期性邊界條件以及對陣列結構(包含n× n單胞)直接建模所得泊松比對比圖 (p =0.6,q=1,φ=90°)Fig.B1 The effective Poisson’s ratio obtained by the PBCs-based homogenization method and calculated from lattice structures consisting of n× nunit-cell array (p =0.6,q=1,φ=90°)