閉孔泡沫金屬幾種不同建模方法的對(duì)比性研究

郭亞周, 劉小川, 白春玉, 鄭志軍, 王計(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)功能型材料，在受壓縮過(guò)程中能夠在較大范圍內(nèi)保持相對(duì)恒定的壓縮應(yīng)力，因此具有較好的能量吸收特性。目前泡沫金屬作為一種理想的輕質(zhì)吸能防護(hù)工程材料已經(jīng)廣泛應(yīng)用于航空航天、汽車、建筑、艦船等領(lǐng)域中[1-3]。

有限元數(shù)值模擬作為一種解決工程實(shí)際問(wèn)題的有利計(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é)者逐步開(kāi)展了閉孔泡沫金屬的細(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)，通過(guò)模型分析后發(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)無(wú)空隙堆疊[10]。宋延澤等[11]基于Kelvin模型研究了泡沫金屬的動(dòng)態(tài)壓縮性能，研究結(jié)果表明，隨著相對(duì)密度的變化，沖擊速率和致密應(yīng)變能夠?qū)ψ冃文Ｊ疆a(chǎn)生較大影響。為了彌補(bǔ)Kelvin模型是由胞元規(guī)則堆疊的不足，學(xué)者逐漸提出了Voronoi模型，Voronoi模型包括二維和三維兩種類型。Sun等[12]開(kāi)發(fā)出了二維Voronoi模型用來(lái)描述泡沫金屬的微觀胞元結(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ì)較少。因此，本工作挑選出目前較為常用的三種模型來(lái)模擬泡沫金屬的壓縮過(guò)程，對(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)過(guò)程，不少學(xué)者通過(guò)構(gòu)建泡沫金屬二維模型來(lái)研究其在不同工況下的變形模式和應(yīng)力響應(yīng)規(guī)律。

針對(duì)二維泡沫金屬模型，一部分研究人員采用了二維Voronoi模型來(lái)表征泡沫金屬的細(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ù)。通過(guò)改變k和t0的值即可生成不同體積和不同規(guī)則度的Voronoi模型。

圖 3 Voronoi模型Fig. 3 Voronoi model

Voronoi模型的生成方式和泡沫金屬發(fā)泡過(guò)程類似，因此三維的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)成，通過(guò)改變殼體厚度來(lái)改變模型的相對(duì)密度，如式（3）所示：

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

1.3 Kelvin模型

如圖4所示，Kelvin模型是由正十四面體堆疊陣列而得，該模型構(gòu)建方法比Voronoi模型簡(jiǎn)單，目前仍然有不少學(xué)者采用這種模型來(lái)表征泡沫金屬的結(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é)果上來(lái)看，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來(lái)描述其力學(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é)果幾乎沒(méi)有影響。

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

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

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

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

如圖6所示，提取出二維隨機(jī)胞元模型的應(yīng)力應(yīng)變曲線與文獻(xiàn)[5]相同相對(duì)密度下的泡沫金屬試樣壓縮實(shí)驗(yàn)結(jié)果對(duì)比，利用Miltz等[16]提出的吸能效率法來(lái)計(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)過(guò)程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所示，通過(guò)仿真計(jì)算后可得相對(duì)密度為0.1的Kelvin模型泡沫金屬在準(zhǔn)靜態(tài)下的壓縮變形過(guò)程，如圖8所示為相對(duì)密度0.1的Voronoi模型泡沫金屬在準(zhǔn)靜態(tài)下的壓縮變形過(guò)程。由壓縮變形圖中可以看出，三維模型能夠便捷清晰地看出泡沫金屬在壓縮過(guò)程中試樣整體的變形和宏觀局部失效情況。

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

在試樣基本上被壓實(shí)階段，由圖7（e）所示，Kelvin模型壓縮過(guò)程每一層的胞元被較為規(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模型具有較大的誤差。通過(guò)計(jì)算得出三維模型的平臺(tái)應(yīng)力和致密應(yīng)變?nèi)绫?所示，Voronoi模型與實(shí)驗(yàn)之間的誤差幾乎可以忽略不計(jì)，因此更能表征泡沫金屬的壓縮性能。

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

從結(jié)構(gòu)特征角度來(lái)講，二維隨機(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)特征方面來(lái)說(shuō)，二維隨機(jī)胞元模型和Voronoi模型都符合真實(shí)泡沫金屬胞元的隨機(jī)性，致密應(yīng)變誤差較小；而Voronoi模型和Kelvin模型整體空間結(jié)構(gòu)與真實(shí)泡沫金屬更相似，平臺(tái)應(yīng)力誤差較小。

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

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)變誤差較小；而Voronoi模型和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é)性能。