閉孔泡沫金屬變形模式的有限元分析*

李妍妍，鄭志軍，虞吉林，王長峰

(中國科學(xué)技術(shù)大學(xué)中國科學(xué)院材料力學(xué)行為和設(shè)計重點實驗室，安徽 合肥 230026)

泡沫金屬已廣泛用作能量吸收材料，很有必要了解其動態(tài)力學(xué)行為。通過細觀有限元模擬能有效地研究隨機性或缺陷對泡沫金屬力學(xué)性能的影響。

已有的有限元模擬多建立在二維情況下。D.Ruan等[1]研究了不同沖擊速度下規(guī)則蜂窩的變形模式和平臺應(yīng)力，并觀察到隨著速度的增加出現(xiàn)X型、V型和I型等3種變形模式。Z.J.Zheng等[2-3]在研究隨機蜂窩的動態(tài)沖擊中從變形機理的角度提出了變形模式的分類，即準靜態(tài)均勻模式、過渡模式和沖擊模式。G.W.Ma等[4]則采用隨機模式、過渡模式和漸進模式來分類。在蜂窩面內(nèi)沖擊性能的進一步研究中，劉穎等[5]和胡玲玲等[6]分別考慮了缺陷分布不均勻性、胞壁夾角的影響。

采用三維細觀有限元模擬更能反映真實泡沫的性能。S.A.Meguid等[7]模擬了規(guī)則閉孔泡沫在動態(tài)壓縮下變形出現(xiàn)的方式。宋延澤等研究了規(guī)則正十四面體泡沫的能量吸收性能[8]，并進一步研究了隨機性對閉孔Voronoi模型的平臺應(yīng)力的影響[9]。王鵬飛等[10]在改進的Hopkison壓桿上，開展閉孔泡沫金屬的動態(tài)沖擊中，發(fā)現(xiàn)變形模式對實驗影響很大，高速沖擊下得到的結(jié)果不能真實反映材料的動態(tài)力學(xué)性能。因此，有必要對泡沫金屬的變形模式和模式轉(zhuǎn)變的臨界速度進行深入研究。

本文中，采用閉孔泡沫金屬的細觀有限元模型，研究相對密度和隨機性對變形模式的影響，給出模式轉(zhuǎn)變的臨界沖擊速度的計算方法，并與文獻中存在的理論預(yù)測公式進行比較和討論。

1 數(shù)值計算模型

1.1 三維隨機Voronoi模型

基于三維隨機Voronoi技術(shù)構(gòu)建閉孔泡沫金屬的計算模型。在指定的區(qū)域內(nèi)隨機生成N個點作為成核點，各成核點之間的最小距離dmin限制為：

(1)

式中：d0和Vcell分別為規(guī)則正十四面體模型中成核點間的最小距離和單胞的體積，k為不規(guī)則因子。將這些點復(fù)制到指定區(qū)域的周圍區(qū)域中，調(diào)用Matlab的子程序Voronoin函數(shù)生成Voronoin構(gòu)型，并保留指定區(qū)域內(nèi)的構(gòu)型作為計算模型。本文中，指定區(qū)域大小為16 mm×12 mm×12 mm，成核點305個，不規(guī)則因子取0.3，一個計算樣本如圖1所示。不考慮胞內(nèi)氣壓的影響，并采用均一的胞壁厚度，由給定的相對密度確定胞壁厚度。計算中，選用了4種相對密度ρ=0.05,0.10,0.15,0.20。

1.2 有限元模擬

采用有限元軟件Abaqus/Explicit模擬泡沫金屬的動態(tài)沖擊。將泡沫金屬模型沿長軸方向上的2個端面分別與2個剛性板接觸。一個剛性板固定，另一個剛性板以恒定的速度沖擊泡沫金屬模型。采用S3R和S4R殼單元對泡沫金屬模型進行網(wǎng)格劃分，經(jīng)過收斂性分析，網(wǎng)格尺寸為約0.03 mm，并調(diào)整較小的單元以確保計算效率，共劃分成近10 000個S3R殼單元和38 000個S4R殼單元。基體材料本構(gòu)設(shè)為剛性-理想塑性：密度為2.77 t/m3，楊氏模量為69 GPa，泊松比為0.3，屈服強度為170 MPa。模型中所有可能的接觸均施加接觸條件，摩擦系數(shù)設(shè)為0.02。

1.3 準靜態(tài)行為和重要參數(shù)

圖2為具有不同相對密度的泡沫金屬的準靜態(tài)名義應(yīng)力-應(yīng)變曲線，計算時采用的加載速度為1 m/s。壓實應(yīng)變εD和平臺應(yīng)力σpl是衡量泡沫金屬準靜態(tài)壓縮行為的重要參數(shù)。壓實應(yīng)變定義為對應(yīng)于能量效率函數(shù)取到最大時的應(yīng)變[11]，能量效率函數(shù)和平臺應(yīng)力分別定義為：

(2)

(3)

式中：εy為屈服應(yīng)變，其值通過觀察準靜態(tài)應(yīng)力-應(yīng)變曲線給出，見表1。根據(jù)準靜態(tài)應(yīng)力-應(yīng)變曲線，進一步求出了4種相對密度下的屈服應(yīng)力、壓實應(yīng)力和應(yīng)變以及平臺應(yīng)力，如表1所示，其中平均值和方差由4個具有不同隨機性的樣本的計算結(jié)果統(tǒng)計后給出。

圖1 三維隨機Voronoi模型Fig.1 A random 3D Voronoi model

圖2 準靜態(tài)應(yīng)力-應(yīng)變曲線Fig.2 Quasi-static compression strain-stress curves

ρεyσy/MPaεDσD/MPaσpl/MPa0．050．0052．6±0．030．620±0．0213．1±0．22．79±0．030．100．0055．5±0．100．624±0．0179．1±0．47．30±0．100．150．01011．2±0．100．620±0．01016．3±0．512．90±0．200．200．02016．2±0．200．610±0．01223．8±0．819．10±0．20

2 計算結(jié)果和討論

2.1 變形模式和名義應(yīng)力-應(yīng)變曲線

隨著沖擊速度的增加，泡沫金屬的壓縮呈現(xiàn)出不同的變形模式。圖3是相對密度為0.05的泡沫金屬在3種沖擊速度下的變形情況，其中左側(cè)為沖擊端。根據(jù)變形的特點，可以采用準靜態(tài)均勻模式、過渡模式和沖擊模式[2-3]對變形進行分類。當(dāng)沖擊速度較低(如v=1 m/s)時，起初試件各部位都有變形，而試件兩端稍微明顯，進而誘導(dǎo)一些隨機變形帶，直至試件全部被壓實。在整個壓縮的過程中，試件兩端的應(yīng)力基本平衡，可以預(yù)期試件中各橫截面上的應(yīng)力均勻，因而將該變形模式命名為準靜態(tài)均勻模式。當(dāng)沖擊速度較高(如v=190 m/s)時，明顯的變形局部化特征出現(xiàn)在沖擊端，在該區(qū)域胞元被壓實，并以逐層壓潰的方式向支撐端傳播，而除了壓實區(qū)之外的部位幾乎沒有變形，相應(yīng)的變形模式稱為沖擊模式。對于中等的沖擊速度(如v=70 m/s)，變形帶較集中在沖擊端，而靠近支撐端的部分也存在一定程度的變形，該變形模式稱為過渡模式。

圖3 變形模式Fig.3 Deformation modes

圖4為試件在3種不同沖擊速度下平滑后的名義應(yīng)力-應(yīng)變曲線。由于原始曲線震蕩比較大，已采用B樣條函數(shù)的最小二乘法將數(shù)據(jù)進行了擬合處理?？梢?，隨著沖擊速度的提高，沖擊端的應(yīng)力平臺有明顯的增強，而支撐端的應(yīng)力的平臺區(qū)對沖擊速度不敏感。隨著沖擊速度的提高，支撐端得初始屈服應(yīng)力在較大的名義應(yīng)變下才達到。

圖4 沖擊端和支撐端的名義應(yīng)力-應(yīng)變曲線Fig.4 Nominal strain-stress curves at the impact and support ends

圖5 應(yīng)力均勻性指標(biāo)隨沖擊速度的變化Fig.5 Stress uniformity index versus impact velocity

2.2 臨界沖擊速度和變形模式圖

應(yīng)力增強是泡沫材料動態(tài)沖擊行為的一個典型的特征。當(dāng)沖擊速度較低時，沖擊端的應(yīng)力和支撐端的應(yīng)力基本上平衡，而隨著沖擊速度的提高，沖擊端因慣性誘導(dǎo)的沖擊增強導(dǎo)致其應(yīng)力明顯高于支撐端的應(yīng)力。因此，可以通過判斷泡沫兩端的應(yīng)力是否平衡，確定第一臨界沖擊速度。應(yīng)力均勻性指標(biāo)[3,12]定義為：

(4)

圖6 R與ε的關(guān)系曲線Fig.6 Curves of R varying with ε

圖7 變形局部化指標(biāo)隨沖擊速度的變化Fig.7 Deformation localization index versus impact velocity

變形局部化是泡沫材料動態(tài)沖擊行為的另一個典型的特征。當(dāng)沖擊速度很高時，變形基本上集中在沖擊端，而靠近支撐端的泡沫基本上不變形。為刻畫這種局部化變形行為，并考慮將來實驗測量的方便，首先定義2個應(yīng)變量：ε1表示泡沫試件中部到?jīng)_擊端面部分的名義應(yīng)變，ε2表示泡沫試件中部到支撐面部分的名義應(yīng)變。試件總體的名義應(yīng)變ε=(ε1+ε2)/2。如果整個試件一直處于均勻變形的狀態(tài)，顯然有ε1=ε2；而如果試件處于完全局部化的變形狀態(tài)，沖擊波陣面從沖擊端向支撐端傳播，在波陣面經(jīng)過試件中部以前都有ε2= 0，此后ε2增大，但不大于ε1?？梢越乒浪?，波陣面到達試件中部時整個試件的名義應(yīng)變?yōu)閴簩崙?yīng)變的一半，即Δε=εD/2。圖6給出了相對密度為0.1的試件、變形局部化完全集中在沖擊端的理想情況下ε1和ε2的相對變化量R=(ε1-ε2)/(ε1+ε2) 隨試件名義應(yīng)變ε的變化曲線，其中εD是由式(2)求極值得出；圖6還給出了試件在3種沖擊速度下的結(jié)果。低速沖擊下，試件近似處于均勻變形狀態(tài)；而高速沖擊下，試件近似處于完全局部化變形狀態(tài)?；谏鲜龇治?，我們引入變形局部化指標(biāo)：

(5)

如果整個試件一直處于均勻變形的狀態(tài)，IDL恒為0；而如果試件處于完全局部化的變形狀態(tài)，IDL恒為1。取IDL= 0.9確定第二臨界沖擊速度。當(dāng)IDL超過0.9時，泡沫的變形從過渡模式轉(zhuǎn)變?yōu)闆_擊模式。圖7給出了4種相對密度下變形局部化指標(biāo)隨沖擊速度的變化曲線?？梢?，隨著相對密度的提高，第二臨界沖擊速度相應(yīng)地有所提高。

圖8 變形模式圖Fig.8 Deformation mode map

以相對密度和沖擊速度為參數(shù)建立泡沫金屬的變形模式圖，如圖8所示?？紤]了4種相對密度，每種相對密度有4個隨機樣本。統(tǒng)計結(jié)果顯示，第一臨界沖擊速度隨著相對密度的提高略有變化，但提高量比統(tǒng)計方差小，因此第一臨界沖擊速度對相對密度不太敏感；第二臨界沖擊速度隨相對密度的提高有較明顯提高。這2個臨界沖擊速度的變化曲線將變形模式圖分成3個區(qū)域，分別對應(yīng)于準靜態(tài)均勻模式、過渡模式和沖擊模式。

2.3 臨界沖擊速度

以往存在一些關(guān)于模式轉(zhuǎn)變的臨界沖擊速度的預(yù)測式，但是由于模式劃分不明確、材料參數(shù)選取差異等可能的原因，這些式子并未得到有效驗證，有的理論推導(dǎo)甚至存在錯誤。A.H?nig等[14]基于陷波理論給出了一個局部化變形開始形成時的臨界速度：

(6)

式中：c(ε)為波速，ρ0為泡沫密度，εcr為壓潰應(yīng)變，它對應(yīng)于在準靜態(tài)應(yīng)力-應(yīng)變曲線上首次達到切線模量為零的狀態(tài)。P.J.Tan等[15]提出了一個剛性-理想塑性-鎖定(R-PP-L)模型，描述泡沫金屬沖擊壓縮中的激波。Z.J.Zheng等[16]通過建立剛性-線性強化塑性-鎖定(R-LHP-L)模型，給出了第一和第二臨界沖擊速度的表達式：

vc1=σ0/(9ρ0c1)

(7)

(8)

觀察不同速度下的變形模式，再與圖8對比，由式(6)、(7)和(8)給出的結(jié)果都偏小。由于鎖定應(yīng)變的定義對R-PP-L和R-LHP-L模型的理論預(yù)測起關(guān)鍵性的影響，下面將著重討論鎖定應(yīng)變。

2.4 鎖定應(yīng)變

以壓實應(yīng)變εD作為鎖定應(yīng)變εL，在應(yīng)用中發(fā)現(xiàn)往往偏小[16]，在圖6中也已體現(xiàn)了這種狀況。我們已經(jīng)運用應(yīng)力均勻化指標(biāo)和變形局部化指標(biāo)確定了臨界沖擊速度，結(jié)合式(8)可以反過來求出鎖定應(yīng)變的大小。由式(8)得到關(guān)系式：

(9)

圖9 鎖定應(yīng)變、壓實應(yīng)變和完全密實應(yīng)變Fig.9 Locking strain, densification strain and full densification strain

式(9)需聯(lián)立準靜態(tài)應(yīng)力-應(yīng)變曲線進行數(shù)值求解，所得的結(jié)果與壓實應(yīng)變εD、完全密實應(yīng)變εmax=1-ρ的比較見圖9。由圖9可見，由式(9)得到的鎖定應(yīng)變大于壓實應(yīng)變、小于完全密實應(yīng)變。這表明，按吸能效率定義的壓實應(yīng)變并不與沖擊模式對應(yīng)，用它作為鎖定應(yīng)變必然導(dǎo)致較大誤差。事實上，沖擊模式下胞孔逐層壓潰，鎖定應(yīng)變隨沖擊速度的提高而增大。由式(9)得到的鎖定應(yīng)變對應(yīng)于胞孔逐層壓潰的最小壓實應(yīng)變，作為描述沖擊模式下泡沫金屬響應(yīng)的帶鎖定階段的近似模型中的鎖定應(yīng)變，比較符合實際。因此，將變形局部化指標(biāo)和式(8)相結(jié)合是確定鎖定應(yīng)變εL的可能方案。

3 結(jié) 論

采用三維Voronoi技術(shù)建立了閉孔泡沫金屬模型，并使用有限元軟件ABAQUS/Explicit進行了動態(tài)沖擊模擬。泡沫金屬在不同的速度沖擊下會出現(xiàn)3種不同的變形模式：準靜態(tài)均勻模式、過渡模式和沖擊模式。引入應(yīng)力均勻化指標(biāo)和變形局部化指標(biāo)，確定了模式轉(zhuǎn)變的臨界沖擊速度，并建立了以沖擊速度和相對密度為坐標(biāo)軸的變形模式圖。計算結(jié)果顯示，對所研究的三維Voronoi閉孔泡沫金屬模型，第一臨界沖擊速度對相對密度不敏感，為約52 m/s；而第二臨界沖擊速度隨著相對密度的提高而提高。并且討論已有的臨界沖擊速度預(yù)測式，提出了一種確定鎖定應(yīng)變的方案，其值介于壓實應(yīng)變和完全密實應(yīng)變之間。

[1] Ruan D, Lu G, Wang B, et al. In-plane dynamic crushing of honeycombs: A finite element study[J]. International Journal of Impact Engineering, 2003,28(2):161-182.

[2] Zheng Z J, Yu J L, Li J R. Dynamic crushing of 2D cellular structures: A finite element study[J]. International Journal of Impact Engineering, 2005,32(4):650-664.

[3] Liu Y D, Yu J L, Zheng Z J, et al. A numerical study on the rate sensitivity of cellular metals[J]. International Journal of Solids and Structures, 2009,46(22):3988-3998.

[4] Ma G W, Ye Z Q, Shao Z S. Modeling loading rate effect on crushing stress of metallic cellular materials[J]. International Journal of Impact Engineering, 2009,36(6):775-782.

[5] 劉穎，張新春．缺陷分布不均勻性對蜂窩材料面內(nèi)沖擊性能的影響[J]．爆炸與沖擊，2009,29(3):237-242.

Liu Ying, Zhang Xin-chun. Effects of inhomogeneous distribution of defects on in-plane dynamic properties of honeycombs[J]. Explosion and Shock Waves, 2009,29(3):237-242.

[6] 胡玲玲,尤帆帆．鋁蜂窩的動態(tài)力學(xué)性能及影響因素[J]．爆炸與沖擊，2012,32(1):23-28.

Hu Ling-ling, You Fang-fang. Dynamic mechanical honeycomb and properties of aluminum its effect factors[J]. Explosion and Shock Waves, 2012,32(1):23-28.

[7] Meguid S A, Cheon S S, El-Abbasi N. FE modelling of deformation localization in metallic foams[J]. Finite Elements in Analysis and Design, 2002,38(7):631-643.

[8] 宋延澤,李志強,趙隆茂．基于十四面體模型的閉孔泡沫材料動態(tài)力學(xué)性能的有限元分析[J]．爆炸與沖擊，2009,29(1):49-55.

Song Yan-ze, Li Zhi-qiang, Zhao Long-mao. Finite element analysis of dynamic crushing behaviors of closed-cell foams based on a tetrakaidecahedron model[J]. Explosion and Shock Waves, 2009,29(1):49-55.

[9] Song Y Z, Wang Z H, Zhao L M, et al. Dynamic crushing behavior of 3D closed-cell foams based on Voronoi random model[J]. Materials and Design, 2010,31(9):4281-4289.

[10] 王鵬飛，徐松林，鄭航，等．變形模式對多孔金屬材料SHPB實驗結(jié)果的影響[J]．力學(xué)學(xué)報，2012,44(5):928-932.

Wang Peng-fei, Xu Song-lin, Zheng Hang, et al. Influence of deformation modes on SHPB experimental results of cellular metal[J]. Chinese Journal of Theoretical and Applied Mechanics, 2012,44(5):928-932.

[11] Tan P J, Reid S R, Harrigan J J, et al. Dynamic compressive strength properties of aluminium foams: Experimental data and observations[J]. Journal of the Mechanics and Physics of Solids, 2005,53(10):2174-2205.

[12] Yu J L, Zheng Z J. Dynamic crushing of 2D cellular metals: Microstructure effects and rate-sensitivity mechanisms[J]. Acta Mechanica Sinica, 2010,23(suppl):45-55.

[13] Deshpande V S, Fleck N A. High strain rate compressive behaviour of aluminium alloy foams[J]. International Journal of Impact Engineering, 2000,24(3):277-298.

[14] H?nig A, Stronge W J. In-plane dynamic crushing of honeycomb: Crush band initiation and wave trapping[J]. International Journal of Mechanical Sciences, 2002,44(8):1665-1696.

[15] Tan P J, Reid S R, Harrigan J J, et al. Dynamic compressive strength properties of aluminium foams: ‘Shock’ theory and comparison with experimental data and numerical models[J]. Journal of the Mechanics and Physics of Solids, 2005,53(10):2206-2230.

[16] Zheng Z J, Liu Y D, Yu J L, et al. Dynamic crushing of cellular materials: Continuum-based wave models for the transitional and shock modes[J]. International Journal of Impact Engineering, 2012,42:66-79.