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

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

（1.中國飛機強度研究所 結(jié)構沖擊動力學航空科技重點實驗室，西安 710065；2.中國科學技術大學 中國科學院材料力學行為與設計重點實驗室，合肥 230026）

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

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

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

1 模型構建

1.1 二維隨機胞元模型

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

針對二維泡沫金屬模型，一部分研究人員采用了二維Voronoi模型來表征泡沫金屬的細觀結(jié)構，雖然二維Voronoi模型能夠?qū)崿F(xiàn)較高的孔隙率，但是二維情況下的Voronoi模型與實際的孔型和相關分布具有一定差異性。基于此，本工作根據(jù)一個閉孔泡沫鋁樣件截面實際胞元的形狀和分布情況（圖1），構建了泡沫金屬的二維隨機胞元模型（圖2）。模型的尺寸與泡沫鋁試件的截面尺寸保持一致，平面尺寸為30 mm × 30 mm，厚度設為0.02 mm，模型的平均孔徑為4 mm。其中若取泡沫金屬基體材料密度為ρs，泡沫金屬的密度為ρ，則泡沫金屬的相對密度為：ρ* =ρ/ρs，此時二維隨機胞元模型的相對密度為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技術能夠生成具有均勻胞壁厚度的閉孔泡沫金屬模型，在給定的體積為V的區(qū)域內(nèi)投放N個核點，其中任意兩個核點之間的距離被約束為大于給定的最小距離tmin，最小距離定義為：

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

圖 3 Voronoi模型Fig. 3 Voronoi model

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

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

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

1.3 Kelvin模型

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

圖 4 Kelvin模型Fig. 4 Kelvin Model

2 模型壓縮分析方法

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

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

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

3.1 二維隨機胞元模型

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

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

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

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

當壓縮應變增大到一定值之后，試樣內(nèi)部的胞元孔隙都基本上被壓實，整個試件上胞孔的上下孔壁之間起到互相支撐的作用，這個階段試件的壓縮應力開始急劇上升并開始進入致密化階段。

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

圖 5 二維隨機胞元模型準靜態(tài)壓縮響應過程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所示，通過仿真計算后可得相對密度為0.1的Kelvin模型泡沫金屬在準靜態(tài)下的壓縮變形過程，如圖8所示為相對密度0.1的Voronoi模型泡沫金屬在準靜態(tài)下的壓縮變形過程。由壓縮變形圖中可以看出，三維模型能夠便捷清晰地看出泡沫金屬在壓縮過程中試樣整體的變形和宏觀局部失效情況。

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

表 2 二維隨機胞元模型與實驗壓縮平臺應力和致密應變Table 2 2D random cell model and experimental compression platform stress and compact strain

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

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

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

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

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

圖 8 Voronoi模型準靜態(tài)壓縮響應過程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）和實驗壓縮結(jié)果對比Fig. 9 Comparison of 3D model（ρ* = 0.1）and experimental compression results

表 3 三維模型與實驗壓縮平臺應力和致密應變Table 3 3D model exprimental compression platform stress and compact strain

3.3 差異性分析

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

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

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

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

4 結(jié)論

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

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

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