高速沖擊下多胞材料細觀結構中波傳播規律

李偉，王鵬，李佳，羅景潤，李明海

（中國工程物理研究院總體工程研究所，四川 綿陽 621900）

泡沫材料通過自身細觀胞孔結構壓縮，將幅值高且脈沖短的爆炸或沖擊瞬時載荷轉變為幅值低且脈沖長的載荷，從而達到被防護物可承受的應力范圍。因此廣泛應用于航天飛行器、防護工程以及汽車撞擊防護設備[1-3]，有效提高設備安全性。

為了更好地發揮泡沫材料在抗爆抗沖擊設施中的緩沖作用，不僅要從宏觀層面研究泡沫材料吸收沖擊能量的能力，更要從細觀層面研究沖擊波在泡沫材料中的傳播和衰減規律。與密實介質不同的是，泡沫材料由于很低的波阻抗使得傳統的SHPB壓桿和直接沖擊等試驗手段較難獲取泡沫類多孔材料中的沖擊波波陣面，即無法獲取沖擊波在泡沫材料中的傳播規律[4-9]。得益于DIC數字相關技術的輔助，聚合物基泡沫材料在中速撞擊下的沖擊波波陣面被成功觀察和識別，從而得到沖擊波波速[4,8,10-13]。即使獲取了聚合物基泡沫材料的波陣面，也只是其在宏觀試件表面的體現，對于試件內部的沖擊波波陣面狀態依然未知。此外，對于金屬基泡沫材料的波陣面獲取依然未有試驗上的報道[14-15]。鑒于金屬基泡沫材料，尤其是泡沫鋁材料，在防護結構上的大量應用和廣泛前景，獲取沖擊波在材料內部的傳播衰減規律，對于優化泡沫材料細觀結構和宏觀構型，進而提升緩沖性能具有重要意義。

基于真實泡沫細觀結構的數值模擬為上述波陣面的獲取提供了一條新的途徑，該方法可以獲取泡沫材料動力學行為的細觀信息，為波陣面的獲取提供了全面的觀察數據。其中3D-Voronoi細觀結構模擬真實閉孔或開孔泡沫結構得到廣泛應用[16-23]。細觀尺度上的變形局部化為多胞材料在動態載荷下的典型特征，名義應力-應變關系無法表征胞元層次上的細致變形，但由于多胞材料細觀結構復雜，變形高度非均勻，大大提升了其局部變形的表征難度。近年來，廖深飛等基于多胞材料的數值計算離散位移場、連續介質力學Lagrange應變張量和最小二乘法離散局部變形梯度法，提出了局部應變場計算方法，較好地解決了多胞材料宏觀尺度上的局部應變場表征問題[19,23-27]。

文中借鑒上述3D-voronoi細觀結構模擬真實閉孔或開孔泡沫結構的數值模擬方法，并結合局部應變場計算方法，開展了高速沖擊下沖擊波在閉孔泡沫鋁中的沖擊波波陣面表征。采用上述方法，對有限元與局部應變法得到的沖擊波波陣面進行了對比，并分析了變形模式。在此基礎上，開展了基于細觀泡沫結構的沖擊波波陣面傳播過程分析，并劃分階段，對沖擊波波速的演化規律進行了分析。最后對沖擊波波陣面劃分階段的合理性進行了驗證。

1 基于3D-Voronoi細觀泡沫結構的數值分析方法

為了更真實地描述泡沫鋁材料的細觀幾何構型，構建了3D-Voronoi隨機細觀結構有限元模型。宏觀泡沫的區域尺寸為80 mm×20 mm×20 mm，并在該區域內隨機撒下600個按照均勻概率分布的成核點，控制成核點間的距離δ不小于給定的距離δmin，如式（1）所示。

式中：k為泡沫材料中表征細觀結構的不規則性；Vcell為十四面體Kelvin胞棱拓撲結構單胞的體積。

滿足以上條件產生成核點后，通過Voronoi函數生成胞棱、胞面幾何信息，如圖1所示。采用S3R單元對胞面進行網格劃分，網格特征長度為0.3 mm。假設泡沫鋁具有均勻厚度，為0.2 mm。在泡沫鋁模型長軸方向的上下表面添加剛性表面，下表面進行固定，上表面承受爆炸載荷作用。在泡沫鋁自身內部、與剛性平面之間施加自接觸，并設置摩擦系數（0.02）。整個泡沫鋁模型含有600個胞元，不規則度為0.3?；牟捎脧椥卫硐胨苄圆牧夏Ｐ?，密度為2.77 g/cm3，楊氏模量為69 GPa，泊松比為0.3，屈服強度為170 MPa。因此得到泡沫鋁細觀結構的數值仿真模型，如圖2所示。采用一端沖擊壓縮（簡稱沖擊端），另外一端固支（簡稱固支端），其中沖擊端質量塊為8 g，初始沖擊速度為150、200、250 m/s。沖擊端和固支端均采用剛性材料。

圖1 泡沫鋁細觀結構幾何圖Fig.1 The diagram of aluminum foam mesoscopic structure

圖2 泡沫鋁細觀結構數值仿真模型Fig.2 The numerical simulation model of aluminum foam mesoscopic structure

2 細觀泡沫結構沖擊波波陣面的局部應變提取法

根據連續介質力學，采用應變張量E對局部變形的表征，如式（2）所示。

式中：F為變形梯度；I為單位張量；T為張量的轉置。

在多胞材料中，采用基于最小二乘法的離散局部變形梯度來描述多胞性帶來的離散位移場，從而得到多胞材料的局部應變場。定義泡沫材料變形前的節點構型為Ω0，變形后的節點構型為當前構型Ω1。節點i和其相鄰節點j在構型Ω0和構型Ω1中的相對位移矢量分別為Uij和uij，如式（3）和（4）所示：

式中：Xi為節點i在構型Ω0中的位置矢量；Xj為相鄰節點j在構型Ω0中的位置矢量；xi為節點i在構型Ω1中的位置矢量；xj為相鄰節點j在構型Ω1中的位置矢量。根據連續介質力學中局部均勻變形假定，節點i和相鄰節點j之間存在唯一的變形梯度F將構型Ω0映射到構型Ω1[24]，即式（5）：

由于多胞材料細觀結構的離散型和復雜性，需要對節點i及其領域一定范圍內其他節點的所有構型映射變形梯度進行最小二乘法，得到節點i的最佳變形梯度Fi，如式（6）所示。

若式(6)分子中的行列式為0，則假設此時節點i處于零應變狀態。其中節點i的領域半徑R設定為0.5r（r為胞元平均半徑）。將各節點的最佳變形梯度帶入式（6），可得各節點的局部應變張量。在笛卡爾直角坐標系(1,2,3)下，可知應變張量對角線分量與局部工程應變的關系，如式（7）所示。

式中：α取值為1，2，3。

3 結果與討論

采用上述方法，得到了泡沫鋁細觀結構有限元計算結果和局部應變三維圖。對有限元與局部應變法得到的沖擊波波陣面進行對比，并分析變形模式。在此基礎上，開展了基于細觀泡沫結構的沖擊波波陣面傳播過程分析，并劃分階段。對沖擊波波速的演化規律進行分析，最后對沖擊波波陣面劃分階段的合理性進行驗證。

3.1 基于有限元與局部應變法的沖擊波波陣面對比

在沖擊端以200 m/s初速沖擊條件下，首先提取泡沫鋁細觀結構數值仿真結果沖擊方向的應力分量場、應變分量場和位移分量場，然后提取局部應變法計算得到的沖擊方向局部應變分量場，對比3個時刻（0.1、0.3、0.5 ms）以上各分量場的分布狀態，如圖3所示。

由圖3可知，由于泡沫細觀結構本身不連續，泡沫鋁細觀結構有限元直接顯示的沿沖擊方向應力分量場、應變分量場和位移分量場呈現出離散狀態，不能像連續介質一樣呈現明顯的應力、應變和位移階段面，因此無法識別沖擊波波陣面。針對以上信息場無法直接表征連續介質尺度下的局部應力和應變，通過第2節介紹的局部應變提取法，獲取了局部應變分量場。從圖3中可以明顯觀察到未擾動介質與已擾動介質的分界面，即沖擊波在泡沫鋁中傳播的應力波波陣面，通過該階段面的Lagrange位置信息和對應時刻，可以得到應力波波陣面隨時間的變化規律。

圖3 泡沫鋁細觀結構有限元結果與局部應變法結果Fig.3 The results comparison of finite model and local strain method for aluminum foam mesoscopic structure

3.2 基于細觀泡沫結構的沖擊波波陣面傳播特征分析

3.2.1 傳播過程及階段劃分

通過3.1的分析實例，利用泡沫鋁細觀結構數值模擬提供的節點位移信息，結合局部應變提取法，獲取了沖擊波在細觀泡沫結構中的應力波波陣面傳播過程，如圖4所示。

圖4 沖擊波在泡沫鋁細觀結構中的沖擊波波陣面傳播過程Fig.4 The shock wave front propagation process in aluminum foam mesoscopic structure

在0.05 ms，細觀泡沫結構內已形成沖擊波波陣面，經過該波陣面，局部應變由0躍遷至約0.7。隨著沖擊波繼續傳播，波陣面的局部應變躍遷至0.8以上，基本已壓實。從沖擊波波陣面在泡沫鋁細觀結構中的演化過程可知，盡管沖擊端以初速200 m/s保持平面沖擊加載，但沖擊波波陣面并未完全保持平面狀態。造成該現象的原因主要是泡沫材料細觀胞孔結構的局部非均勻性，在宏觀上可近似看作平面。在0.35 ms以后，由于沖擊初始能量恒定，且泡沫材料對沖擊能量持續吸收，造成沖擊端產生的右行沖擊波消失，即局部應變尖端面在0.35 ms以后保持Lagrange坐標位置不變。在0.15～0.2 ms，支撐端開始出現明顯的應變突躍，并在隨后的過程中逐漸增加，最大突躍值可達0.65，接近壓實狀態。對于未傳播區域，處于未壓縮狀態，在支撐端附近的泡沫材料未壓縮前，隨右行沖擊波的傳播，未傳播區域逐漸收縮，且該區域的運動速度接近于0。在支撐端附近的泡沫材料開始壓縮后，且右行沖擊波未消失前，該區域左、右側同時收縮，且剩余部分近似以剛體形式向右運動。在右行沖擊波消失后，該區域只有右側收縮，右行波所過之處與未傳播區域作為一體，近似以剛體形式向右運動?；谝陨戏治?，將該過程分為3個階段：沖擊端誘發右行沖擊波、支撐端誘發左行沖擊波且未壓縮區剛性向右運動、右行沖擊波壓縮區域與未壓縮區域一體剛性向右運動。

3.2.2 沖擊波波速演化規律

由圖4分析可知，沖擊波波陣面近似平面，可認為是一維沖擊波傳播。因此在沖擊壓縮垂直平面（截面）內，對應變值進行平均化處理，以表征該截面的應變。在0.05～0.5 ms內，以0.05 ms為時間間隔提取每個時刻的應變隨無量綱拉氏位置的變化曲線，如圖5所示。

圖5 沖擊波在泡沫鋁細觀結構中的波陣面時程曲線Fig.5 The shock wave front with time in aluminum foam mesoscopic structure

由圖5可知，雖然應變曲線沒有呈現完全豎直狀態，但斜率較大，已呈現沖擊波波陣面前后應變突變的特點。通過Zheng[19]的處理方法，可以獲取每個時刻理想波陣面狀態下的右行和左行沖擊波波陣面拉氏位置和速度，如圖6所示。

圖6 沖擊波在泡沫鋁細觀結構中的波陣面位置及速度Fig.6 The Lagrange position and velocity of right-traveling shock wave front in aluminum foam mesoscopic structure: a)right-traveling shock wave; b) left-traveling shock wave

如圖6a所示，右行沖擊波波陣面的拉氏位置逐漸右移，且波陣面從0.05 ms的140.8 m/s以基本恒定的加速度減小，波速曲線平滑。在0.35～0.40 ms，右行沖擊波波速降低為0，且沖擊波波陣面不在拉氏坐標下右移，即此時右行沖擊波消失。此后，右行沖擊波所過之處的泡沫材料中存在恒定的殘余應變。產生以上現象的原因為：沖擊端初始沖擊能量恒定（沖擊塊為剛體，質量固定，屬瞬時強加載條件），在200 m/s的高速沖擊壓縮過程中，沖擊端附近泡沫材料由于慣性效應，產生局部壓縮，并壓實。壓實的泡沫材料與沖擊塊以基本相同的速度繼續沖擊未壓縮區域，由于細觀胞壁結構的屈曲以及大量細觀塑性鉸的產生、傳播和相互作用，泡沫材料對沖擊能量進行吸收，從而造成對未壓縮區域的沖擊波的壓縮能量降低。由應力波波陣面前后的質量、動量和能量守恒關系式（7）—（9）可知，沖擊波波速也會相應下降。

式中：[a]為波后與波前參量值之差，[a]=a--a+，其中a-和a+分別表示波后和波前參量，此處a可替換為速度v、應變ε、應力σ和內能E；D為波陣面速度；ρ0為波前泡沫材料密度；A0為泡沫材料的橫截面積；ΔX為波陣面的寬度。

由于右行沖擊波傳播過后，泡沫材料處于壓實狀態，在此認為各時刻的波后狀態相同，因此泡沫的吸能（內能增加）相同。設定單位體積的吸收能量為η，因此該值為恒定值，則波陣面前后的能量變化可寫成式（10）：

式中：[t]為沖擊波經過波陣面的寬度ΔX所用時間。將式（7）和式（8）帶入式（9），并聯立式（10）可得：

由圖5的右行沖擊波應變值可知，波前、波后的應變差[ε]基本恒定，在假設波陣面寬度不變的條件下，結合式（5）可得沖擊波波陣面的加速度為恒定值。該理論分析結果與圖6中的右行沖擊波波陣面速度曲線在0.05～0.4 ms呈現直線下降（減加速度恒定）相吻合。如圖6b所示，在0.1 ms之前，支撐端附近未產生左行沖擊波。在0.20 ms之后，左行沖擊波速度達到最大值，為56.8 m/s，隨后出現波速下降再上升的現象。在支撐端產生左行沖擊波可能是由于沖擊端沖擊泡沫材料初始，彈性波沿泡沫材料沖擊方向傳播，在支撐端發生剛性壁面反射，反射的應力波超過泡沫材料的屈服極限，因此材料開始動態壓縮，未壓縮區域與支撐端壓縮區域產生速度差，產生內撞擊，并形成左行沖擊波。由于左行沖擊波壓縮材料被吸收能量和不斷的內撞擊，造成了波速下降再上升的現象。為了更好地研究沖擊波在材料中的衰減現象，建立沖擊波波速與波陣面的無量綱拉氏位置坐標間的關系曲線，如圖7所示。

圖7 沖擊波波速隨傳播的衰減Fig.7 The attenuation of shock wave front in aluminum foam mesoscopic structure

由圖7可知，右行沖擊波的波速沿著拉氏位置的傳播呈現出加速衰減的光滑曲線。這是由于右行波陣面減加速度為恒定值造成的，因此該曲線表現出沖擊波波速為自變量，且波陣面無量綱拉氏位置為變量的二次曲線。對于左行沖擊波，波速沿拉氏位置的傳播依然呈現波動現象，同樣該現象是由于左行沖擊波壓縮材料被吸收能量和不斷的內撞擊造成的。

3.2.3 沖擊波波陣面傳播階段驗證

在3.2.1節，將該過程分為3個階段，對于第一和第二階段，是由波陣面傳播現象直接得到的，無需驗證。對于第三階段，認為右行沖擊波壓縮區域在右行沖擊波消失后，呈現向右的剛體運動。為了驗證該判斷，對比了右行波陣面位置與沖擊端壓縮面隨時間的歐拉坐標運動曲線（右行波陣面位置與沖擊端壓縮面之間的部分為右行沖擊波壓縮區域），如圖8所示。

圖8 沖擊端壓縮面與波陣面在歐拉坐標下的間距變化Fig.8 The Euler position of proximal end and shock wave front

由圖8可知，在0.30 ms之前，右行沖擊波壓縮區域持續增大。這是由于隨著右行沖擊波的傳播，波后的壓實狀態泡沫材料不斷積累所致。在0.35 ms之后，右行沖擊波壓縮區域不再增長，保持恒定長度。說明在第三階段，右行沖擊波壓縮區不再變形，而是以剛體形式聯合沖擊塊和未壓縮區域一起向右運動，與左行沖擊波壓縮區左端持續產生內撞擊。

4 結論

1）采用基于3D-Voronoi細觀泡沫結構數值仿真和基于最小二乘法的局部應變梯度法，研究了細觀泡沫結構在高速沖擊下的沖擊波波陣面及壓潰界面傳播規律。

2）通過細觀泡沫結構數值仿真提取的位移場數據，結合局部應變梯度法得到了界面清晰的沖擊波波陣面，而在單純的數值仿真結果中無法識別。基于細觀泡沫結構的沖擊波波陣面傳播過程，不僅得到沖擊端產生的右行沖擊波波陣面，還得到了支撐端誘發的左行沖擊波波陣面。

3）將沖擊波波陣面傳播過程劃分為3個階段，即沖擊端誘發右行沖擊波階段、支撐端誘發左行沖擊波且未壓縮區剛性向右運動階段和右行沖擊波壓縮區域與未壓縮區域一體剛性向右運動階段。

4）通過沖擊波波速的演化規律研究，分別獲取了右行和左行沖擊波波陣面的拉氏位置和波速，通過嚴格理論推導解釋了右行沖擊波波陣面具有恒定減加速度的結論。最后對沖擊波波陣面劃分階段的合理性進行了驗證。

文中對細觀泡沫結構在高速沖擊下的沖擊波波陣面傳播規律做了基礎性研究，為航天飛行器、防護工程以及汽車撞擊防護設備等典型目標結構防護、安全性設計和研究提供了參考依據。