郭亞周, 劉小川, 白春玉, 鄭志軍, 王計(jì)真

（1.中國(guó)飛機(jī)強(qiáng)度研究所 結(jié)構(gòu)沖擊動(dòng)力學(xué)航空科技重點(diǎn)實(shí)驗(yàn)室，西安 710065；2.中國(guó)科學(xué)技術(shù)大學(xué) 中國(guó)科學(xué)院材料力學(xué)行為與設(shè)計(jì)重點(diǎn)實(shí)驗(yàn)室，合肥 230026）

泡沫金屬作為一種多孔、輕質(zhì)的結(jié)構(gòu)功能型材料，在受壓縮過程中能夠在較大范圍內(nèi)保持相對(duì)恒定的壓縮應(yīng)力，因此具有較好的能量吸收特性。目前泡沫金屬作為一種理想的輕質(zhì)吸能防護(hù)工程材料已經(jīng)廣泛應(yīng)用于航空航天、汽車、建筑、艦船等領(lǐng)域中[1-3]。

有限元數(shù)值模擬作為一種解決工程實(shí)際問題的有利計(jì)算工具，能夠在合適的建模方法基礎(chǔ)上較好地還原物理實(shí)驗(yàn)現(xiàn)象，并且隨著建模和分析方法的成熟能夠在一定程度上逐步預(yù)測(cè)實(shí)驗(yàn)結(jié)果，節(jié)約研究成本，縮短研究周期[4]。為研究閉孔泡沫金屬的壓縮力學(xué)性能，擺脫實(shí)驗(yàn)工況的限制，探討泡沫金屬的細(xì)觀胞元結(jié)構(gòu)對(duì)整體力學(xué)性能的影響，不少學(xué)者逐步開展了閉孔泡沫金屬的細(xì)觀建模方法研究。張健等[5]基于微CT掃描影像信息建立泡沫金屬二維細(xì)觀有限元模型，研究了在高速加載下泡沫金屬的壓縮變形機(jī)理。袁本立等[6]從不同類型的泡沫金屬單一胞元模型出發(fā)，研究了結(jié)構(gòu)較為簡(jiǎn)單的單心立方模型、面心立方模型和體心立方模型表征下泡沫金屬的力學(xué)性能，發(fā)現(xiàn)面心立方模型和體心立方模型能夠相對(duì)較好地預(yù)測(cè)泡沫金屬性質(zhì)。Gibson等[7]提出構(gòu)建了Gibson-Ashby泡沫金屬模型，該模型的最大特點(diǎn)在于其為包含12條邊的立體框架結(jié)構(gòu)，通過模型分析后發(fā)現(xiàn)影響泡沫金屬性質(zhì)的因素主要是孔隙率、基體和泡孔結(jié)構(gòu)。劉培生[8]提出了一種泡沫金屬八面體模型，主要是為了克服Gibon-Ashby模型的不足，該模型結(jié)構(gòu)的六個(gè)頂點(diǎn)分別是長(zhǎng)方體兩個(gè)平行平面和其他四條垂直棱邊的中點(diǎn)，目前該模型在扭矩作用下泡沫金屬的工程設(shè)計(jì)強(qiáng)度準(zhǔn)則等方面具有較為廣泛的應(yīng)用[9]。隨著泡沫金屬細(xì)觀建模方法的發(fā)展，不少學(xué)者發(fā)現(xiàn)Kelvin模型能夠較為合理有效地模擬出泡沫金屬的結(jié)構(gòu)和性質(zhì)，Kelvin模型是由十四條等長(zhǎng)的棱長(zhǎng)構(gòu)成的正十四面體，該模型由6個(gè)四邊形和8個(gè)六邊形組成，能夠在空間中實(shí)現(xiàn)無空隙堆疊[10]。宋延澤等[11]基于Kelvin模型研究了泡沫金屬的動(dòng)態(tài)壓縮性能，研究結(jié)果表明，隨著相對(duì)密度的變化，沖擊速率和致密應(yīng)變能夠?qū)ψ冃文Ｊ疆a(chǎn)生較大影響。為了彌補(bǔ)Kelvin模型是由胞元規(guī)則堆疊的不足，學(xué)者逐漸提出了Voronoi模型，Voronoi模型包括二維和三維兩種類型。Sun等[12]開發(fā)出了二維Voronoi模型用來描述泡沫金屬的微觀胞元結(jié)構(gòu)。Zhang等[13]基于Abaqus分析軟件用Voronoi模型研究了軟件參數(shù)和質(zhì)量縮放量、單元類型、單元大小對(duì)模型的模擬結(jié)果。

雖然目前關(guān)于泡沫金屬的細(xì)觀建模方法較多，但是每種建模方法對(duì)于不同場(chǎng)合的適用性卻各不相同，且不同建模方法的合理性和表征準(zhǔn)確性也各有差異。目前的研究中，更多的只是聚焦單一模型定性或定量的研究泡沫金屬力學(xué)性能，而針對(duì)不同建模方法適用性和準(zhǔn)確性的差異對(duì)比研究則相對(duì)較少。因此，本工作挑選出目前較為常用的三種模型來模擬泡沫金屬的壓縮過程，對(duì)比分析三種模型的變形模式、應(yīng)力響應(yīng)和壓縮力學(xué)性能表征能力的差異性，并針對(duì)不同模型特點(diǎn)分別給出其研究適用場(chǎng)合建議。

1 模型構(gòu)建

1.1 二維隨機(jī)胞元模型

由于泡沫金屬本身呈現(xiàn)近似各向同性的結(jié)構(gòu)特征，且在通常的壓縮實(shí)驗(yàn)中并不能夠觀察到試件內(nèi)部的變形情況，因此為了較為直觀地研究泡沫金屬壓縮內(nèi)部結(jié)構(gòu)響應(yīng)過程，不少學(xué)者通過構(gòu)建泡沫金屬二維模型來研究其在不同工況下的變形模式和應(yīng)力響應(yīng)規(guī)律。

針對(duì)二維泡沫金屬模型，一部分研究人員采用了二維Voronoi模型來表征泡沫金屬的細(xì)觀結(jié)構(gòu)，雖然二維Voronoi模型能夠?qū)崿F(xiàn)較高的孔隙率，但是二維情況下的Voronoi模型與實(shí)際的孔型和相關(guān)分布具有一定差異性?；诖?，本工作根據(jù)一個(gè)閉孔泡沫鋁樣件截面實(shí)際胞元的形狀和分布情況（圖1），構(gòu)建了泡沫金屬的二維隨機(jī)胞元模型（圖2）。模型的尺寸與泡沫鋁試件的截面尺寸保持一致，平面尺寸為30 mm × 30 mm，厚度設(shè)為0.02 mm，模型的平均孔徑為4 mm。其中若取泡沫金屬基體材料密度為ρs，泡沫金屬的密度為ρ，則泡沫金屬的相對(duì)密度為：ρ* =ρ/ρs，此時(shí)二維隨機(jī)胞元模型的相對(duì)密度為0.3。

圖 1 泡沫金屬模型 （a）泡沫鋁試樣；（b）二維幾何模型Fig. 1 2D metal foam model （a）aluminum foam sample；（b）2D geometric model

圖 2 泡沫金屬二維有限元網(wǎng)格模型Fig. 2 2D metal foam finite element mesh model

1.2 Voronoi模型

圖3為三維Voronoi模型。三維Voronoi技術(shù)能夠生成具有均勻胞壁厚度的閉孔泡沫金屬模型，在給定的體積為V的區(qū)域內(nèi)投放N個(gè)核點(diǎn)，其中任意兩個(gè)核點(diǎn)之間的距離被約束為大于給定的最小距離tmin，最小距離定義為：

式中：tmin為任意兩個(gè)核點(diǎn)之間的最小距離；k為胞元不規(guī)則度；t0為正十四面體任意兩個(gè)核點(diǎn)之間的最小距離；N為給定的體積V的區(qū)域內(nèi)投放的核點(diǎn)個(gè)數(shù)。通過改變k和t0的值即可生成不同體積和不同規(guī)則度的Voronoi模型。

圖 3 Voronoi模型Fig. 3 Voronoi model

Voronoi模型的生成方式和泡沫金屬發(fā)泡過程類似，因此三維的Voronoi模型相比于二維更能夠還原真實(shí)泡沫金屬胞元間孔棱和孔壁關(guān)系，已有的研究表明，當(dāng)閉孔泡沫金屬的相對(duì)密度低于0.3時(shí)，真實(shí)泡沫金屬的拓?fù)浣Y(jié)構(gòu)與Voronoi模型基本相同[14]，因此Voronoi模型能夠較好地表征泡沫金屬。

本工作模型尺寸為30 mm × 30 mm × 30 mm，胞元直徑為4 mm，胞元不規(guī)則度k為0.5。由于Voronoi模型全部由殼體構(gòu)成，通過改變殼體厚度來改變模型的相對(duì)密度，如式（3）所示：

式中：ρs為泡沫金屬基體材料密度；ρ為泡沫金屬密度；ρ*為泡沫金屬相對(duì)密度；V為模型體積；為模型所有胞壁面積總和；t為模型胞壁厚度。

1.3 Kelvin模型

如圖4所示，Kelvin模型是由正十四面體堆疊陣列而得，該模型構(gòu)建方法比Voronoi模型簡(jiǎn)單，目前仍然有不少學(xué)者采用這種模型來表征泡沫金屬的結(jié)構(gòu)。Kelvin模型本質(zhì)上是不規(guī)則度k為0的Voronoi模型，因此Kelvin模型和Voronoi最大的區(qū)別在于Voronoi在規(guī)則胞元的基礎(chǔ)上結(jié)合實(shí)際的泡沫金屬發(fā)泡生成情況后對(duì)模型添加了擾動(dòng)。雖然Kelvin模型并不具備胞元的隨機(jī)性，但是就目前相關(guān)學(xué)者的研究結(jié)果上來看，Kelvin在某些工況下仍然和實(shí)驗(yàn)具備較好的一致性。

圖 4 Kelvin模型Fig. 4 Kelvin Model

2 模型壓縮分析方法

采用LS-DYNA軟件對(duì)模型進(jìn)行壓縮模擬，其中二維隨機(jī)胞元模型主要采用四邊形網(wǎng)格，單元尺寸為0.5 mm，單元總數(shù)量是12353，模型在平面上施加對(duì)稱約束。Voronoi模型全部由殼體構(gòu)成，主要采用四邊形網(wǎng)格，單元網(wǎng)格尺寸為0.5 mm，單元網(wǎng)格總數(shù)量是186743。Kelvin模型與Voronoi模型網(wǎng)格劃分規(guī)則相同，單元網(wǎng)格總數(shù)量分別為229927和202528。在模型的加載方向兩端放置兩個(gè)剛性板，下剛性板固定，上剛性板加載恒定速率，設(shè)定壓縮量為24 mm，當(dāng)泡沫金屬被壓縮至應(yīng)變?yōu)?.8時(shí)，計(jì)算終止。上下兩個(gè)剛性板變形可以忽略，采用*MAT_RIGID模型，泡沫金屬的基體材料選用*MAT_PLASTIC_ KINEMATIC來描述其力學(xué)行為，該模型為雙線性應(yīng)變硬化模型，另泡沫金屬的基體材料為鋁，鋁基體材料參數(shù)如表1所示，泡沫金屬模型采用*CONTACT_AUTOMATIC_SINGLE_SURFACRE接觸，泡沫鋁模型與剛性板之間采用*CONTACT_ AUTOMATIC_NODES_TO_SURFACE接觸，摩擦系數(shù)設(shè)置為0.1，Song等[15]驗(yàn)證了該參數(shù)對(duì)模型的模擬結(jié)果幾乎沒有影響。

表 1 泡沫鋁基體材料參數(shù)Table 1 Parameters of foam aluminum matrix material

3 模型壓縮結(jié)果分析

3.1 二維隨機(jī)胞元模型

由于二維細(xì)觀模型孔壁孔棱與實(shí)際三維孔之間堆疊后的孔壁孔棱之間具有幾何差異性，同時(shí)為了保證二維隨機(jī)胞元模型的幾何和網(wǎng)格特性，其孔棱不可能非常小，小的孔棱在幾何上較難實(shí)現(xiàn)，同時(shí)也會(huì)導(dǎo)致網(wǎng)格尺寸和計(jì)算精度受到影響，因此二維隨機(jī)胞元模型的相對(duì)密度 ρ*數(shù)值會(huì)偏大。

如圖5所示，通過仿真計(jì)算后可得相對(duì)密度為0.3的泡沫金屬在準(zhǔn)靜態(tài)下的壓縮變形過程。由壓縮變形圖中可以看出，二維隨機(jī)胞元模型能夠便捷清晰地看出泡沫金屬在壓縮過程中內(nèi)部胞孔和胞壁的變形情況。

在泡沫金屬受到壓縮的初級(jí)階段，試件的整體變形較為均勻，試件內(nèi)部胞壁的部分薄弱結(jié)構(gòu)逐漸開始出現(xiàn)應(yīng)力集中，進(jìn)而在這些部位開始集中出現(xiàn)了相應(yīng)的塑性屈曲，眾多屈曲部位的塑性鉸在宏觀上表現(xiàn)為一條或者多條不規(guī)則的變形帶。

隨著壓縮量的增大，變形帶逐漸向周邊延伸，變形帶由條狀逐漸向區(qū)域塊狀轉(zhuǎn)變，此時(shí)孔壁之間會(huì)相互接觸摩擦，下部孔壁為上部分坍塌孔壁提供支撐，使得整體試件的支撐能力能夠在較大的應(yīng)變內(nèi)維持在一個(gè)較為穩(wěn)定的水平，此壓縮過程即為試樣的平臺(tái)壓縮階段。

當(dāng)壓縮應(yīng)變?cè)龃蟮揭欢ㄖ抵?，試樣?nèi)部的胞元孔隙都基本上被壓實(shí)，整個(gè)試件上胞孔的上下孔壁之間起到互相支撐的作用，這個(gè)階段試件的壓縮應(yīng)力開始急劇上升并開始進(jìn)入致密化階段。

如圖6所示，提取出二維隨機(jī)胞元模型的應(yīng)力應(yīng)變曲線與文獻(xiàn)[5]相同相對(duì)密度下的泡沫金屬試樣壓縮實(shí)驗(yàn)結(jié)果對(duì)比，利用Miltz等[16]提出的吸能效率法來計(jì)算出響應(yīng)的平臺(tái)應(yīng)力和致密應(yīng)變，結(jié)果如表2所示。二維隨機(jī)胞元模型與實(shí)驗(yàn)結(jié)果具有良好的一致性，但又由于內(nèi)部胞元結(jié)構(gòu)的差異性，導(dǎo)致其在壓縮初期的彈性變形階段具有一定的誤差，從而使得平臺(tái)段產(chǎn)生了差異性。

圖 5 二維隨機(jī)胞元模型準(zhǔn)靜態(tài)壓縮響應(yīng)過程Fig. 5 Quasi-static compression processes of 2D stochastic cell model （a）ε = 0；（b）ε = 0.15；（c）ε = 0.4；（d）ε = 0.6；（e）ε = 0.75

3.2 Kelvin和Voronoi模型

如圖7所示，通過仿真計(jì)算后可得相對(duì)密度為0.1的Kelvin模型泡沫金屬在準(zhǔn)靜態(tài)下的壓縮變形過程，如圖8所示為相對(duì)密度0.1的Voronoi模型泡沫金屬在準(zhǔn)靜態(tài)下的壓縮變形過程。由壓縮變形圖中可以看出，三維模型能夠便捷清晰地看出泡沫金屬在壓縮過程中試樣整體的變形和宏觀局部失效情況。

在壓縮初始階段，由圖7（b）所示，Kelvin模型的變形多集中在接觸端，而中間位置變形較少。由圖8（b）所示，Voronoi模型在相同應(yīng)變下雖然在接觸端仍然由一定的變形，但是在模型的中間部位則已經(jīng)出現(xiàn)了較為明顯的變形帶，這與泡沫金屬壓縮實(shí)驗(yàn)現(xiàn)象更相符。

表 2 二維隨機(jī)胞元模型與實(shí)驗(yàn)壓縮平臺(tái)應(yīng)力和致密應(yīng)變Table 2 2D random cell model and experimental compression platform stress and compact strain

圖 6 二維隨機(jī)胞元模型（ρ* = 0.3）和實(shí)驗(yàn)壓縮結(jié)果對(duì)比Fig. 6 Comparison of 2D model（ρ* = 0.3）and experimental compression results

圖 7 Kelvin模型準(zhǔn)靜態(tài)壓縮響應(yīng)過程Fig. 7 Quasi-static compression response processes of Kelvin model （a）ε = 0；（b）ε = 0.15；（c）ε = 0.4；（d）ε = 0.6；（e）ε = 0.75

隨著壓縮應(yīng)變的增大，由圖7（c）和圖7（d）所示，Kelvin模型的變形帶逐漸延伸至試樣的中部，且由于Kelvin模型是由正十四面體堆疊而成，其試樣內(nèi)部薄弱結(jié)構(gòu)較為集中且規(guī)律，因此壓縮過程中呈現(xiàn)了較為規(guī)律的屈曲變形帶，變形區(qū)域沿著變形帶逐漸變大。由圖8（c）和圖8（d）所示，Voronoi模型與Kelvin模型略有不同，Voronoi模型由于內(nèi)部胞元更具備隨機(jī)性，其生成過程與泡沫金屬的實(shí)際發(fā)泡過程類似，胞元形狀和實(shí)際試樣的孔型胞元結(jié)構(gòu)特征相符，使得試樣內(nèi)部薄弱結(jié)構(gòu)較為分散且呈現(xiàn)隨機(jī)性，因此隨著壓縮應(yīng)變的增大，Voronoi模型的變形帶相比于Kelvin模型具有更明顯的隨機(jī)性。

在試樣基本上被壓實(shí)階段，由圖7（e）所示，Kelvin模型壓縮過程每一層的胞元被較為規(guī)律地壓實(shí)，最終呈現(xiàn)明顯的分層現(xiàn)象。而如圖8（e）所示，Voronoi模型則更加符合變形帶的隨機(jī)性原則，試樣的各個(gè)部位基本上都發(fā)生了較大的變形。

如圖9所示，分別提取出Kelvin模型和Voronoi模型的應(yīng)力應(yīng)變曲線與相同相對(duì)密度下的泡沫金屬試樣壓縮實(shí)驗(yàn)結(jié)果對(duì)比。由圖9就可以明顯看出，Voronoi模型比Kelvin模型更接近實(shí)驗(yàn)結(jié)果，Kelvin模型屈服段與致密段與實(shí)驗(yàn)結(jié)果相比于Voronoi模型具有較大的誤差。通過計(jì)算得出三維模型的平臺(tái)應(yīng)力和致密應(yīng)變?nèi)绫?所示，Voronoi模型與實(shí)驗(yàn)之間的誤差幾乎可以忽略不計(jì)，因此更能表征泡沫金屬的壓縮性能。

圖 8 Voronoi模型準(zhǔn)靜態(tài)壓縮響應(yīng)過程Fig. 8 Quasi-static compression response processes of Voronoi model （a）ε = 0；（b）ε = 0.15；（c）ε = 0.4；（d）ε = 0.6；（e）ε = 0.75

圖 9 三維模型（ρ* = 0.1）和實(shí)驗(yàn)壓縮結(jié)果對(duì)比Fig. 9 Comparison of 3D model（ρ* = 0.1）and experimental compression results

表 3 三維模型與實(shí)驗(yàn)壓縮平臺(tái)應(yīng)力和致密應(yīng)變Table 3 3D model exprimental compression platform stress and compact strain

3.3 差異性分析

三維模型和二維模型具有一定的差異性，其中最主要的差異性主要體現(xiàn)在以下三個(gè)方面：（1）相對(duì)密度；（2）結(jié)構(gòu)特征；（3）壓縮變形響應(yīng)。

從相對(duì)密度角度來講，三維模型由于全部都是由殼體構(gòu)成，因此通過改變殼體厚度即可在相對(duì)密度小于0.3的范圍內(nèi)任意改變泡沫金屬模型的相對(duì)密度，而二維隨機(jī)胞元模型由于受到孔壁孔棱的幾何限制導(dǎo)致其相對(duì)密度往往不會(huì)小于0.25。在實(shí)際閉孔泡沫金屬的使用場(chǎng)景中，為了達(dá)到質(zhì)輕和高效吸能的目的，采用的閉孔泡沫金屬相對(duì)密度大多較小。因此就相對(duì)密度方面來說，二維隨機(jī)模型更偏向于泡沫金屬的壓縮規(guī)律性研究，而三維模型則更偏向于泡沫金屬的應(yīng)用優(yōu)化研究。

從結(jié)構(gòu)特征角度來講，二維隨機(jī)胞元模型以泡沫金屬試件橫截面為基準(zhǔn)構(gòu)建，由于泡沫金屬具有一定的各向同性屬性，因此在不考慮相對(duì)密度范圍的情況下，二維隨機(jī)胞元模型在可接受的誤差范圍內(nèi)能夠比三維模型更容易構(gòu)建且能夠更便捷地表征泡沫金屬的壓縮性能。三維模型更適合表征真實(shí)的試樣狀況，能夠更好地表征泡沫金屬的宏觀壓縮性能，但是其構(gòu)建方法較為復(fù)雜，且Kelvin模型由于是由規(guī)則正十四面體堆疊而成，導(dǎo)致其壓縮后呈現(xiàn)規(guī)則的分層坍塌，使得其雖然平臺(tái)應(yīng)力誤差不大，但是壓縮致密應(yīng)變誤差卻較大。因此就結(jié)構(gòu)特征方面來說，二維隨機(jī)胞元模型和Voronoi模型都符合真實(shí)泡沫金屬胞元的隨機(jī)性，致密應(yīng)變誤差較小；而Voronoi模型和Kelvin模型整體空間結(jié)構(gòu)與真實(shí)泡沫金屬更相似，平臺(tái)應(yīng)力誤差較小。

從壓縮變形響應(yīng)角度來講，三維模型相比于二維隨機(jī)胞元模型能夠更全面地觀察泡沫金屬壓縮的變形帶產(chǎn)生與演變以及泡沫金屬整體的變形模式變化，同時(shí)可以在多工況下考察模型不同應(yīng)變下的應(yīng)力場(chǎng)分布，其中Voronoi模型由于生成過程與閉孔泡沫金屬發(fā)泡過程類似，內(nèi)部結(jié)構(gòu)更趨真實(shí)合理，其變形相比于Kelvin模型的分層壓實(shí)來說更能體現(xiàn)變形帶的隨機(jī)性和真實(shí)性。二維隨機(jī)胞元模型則較為直觀地觀察出泡沫金屬內(nèi)部胞元和胞壁的變形響應(yīng)情況以及變形帶的整體演變狀況，但是就宏觀變形來講，較三維模型仍有較大不足。

4 結(jié)論

（1）二維隨機(jī)胞元模型由于受到孔棱的幾何限制，其相對(duì)密度比三維模型較大，與真實(shí)常用的輕質(zhì)、吸能效率較好的閉孔泡沫金屬相對(duì)密度具有一定差異。因此二維隨機(jī)模型更適合研究泡沫金屬的壓縮規(guī)律性研究，而三維模型則更適合泡沫金屬壓縮性能的精確預(yù)測(cè)與優(yōu)化研究。

（2）二維隨機(jī)胞元模型相比于三維模型構(gòu)建較為簡(jiǎn)單和便捷。二維隨機(jī)胞元模型和Voronoi模型都符合真實(shí)泡沫金屬胞元的隨機(jī)性，致密應(yīng)變誤差較?。欢鳹oronoi模型和Kelvin模型整體空間結(jié)構(gòu)與真實(shí)泡沫金屬更相似，平臺(tái)應(yīng)力誤差較小。

（3）Voronoi模型相比于Kelvin模型和二維隨機(jī)胞元模型內(nèi)部結(jié)構(gòu)更趨真實(shí)合理，能夠體現(xiàn)出變形帶的隨機(jī)性，變形模式也與實(shí)驗(yàn)一致性較好，能夠更精準(zhǔn)和全面的表征閉孔泡沫金屬的壓縮力學(xué)性能。