葉常青,陳 然,劉桂森,劉靜楠,胡建波,俞宇穎,王 棟,陳開果,沈 耀
(1. 上海交通大學(xué)材料科學(xué)與工程學(xué)院, 上海 200240;2. 中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽 621999;3. 西北核技術(shù)研究院, 陜西 西安 710024;4. 國防科技大學(xué)物理系, 湖南 長沙 410073)
Abstract: In polycrystalline materials, the grain boundaries between grains with different orientations have a great influence on the dynamic response of material under shock loading. On the basis of the single crystal plasticity model, a polycrystal plasticity model containing grain boundary obstacle, geometrically necessary dislocations (GND) and back-stress is established by considering the microscopic mechanism of the interaction between grain boundaries and dislocations. Based on this model, the mechanical response of polycrystalline aluminum with Voronoi geometry under shock loading was studied through simulations. The results show that: (1) the grain boundary elements after the shock loading have a very high residual shear stress, while the shear stress of the element within the grain tends to be zero; (2) a large plastic deformation gradient is found near the grain boundary, resulting in a large number of GND and back stress distributed along the grain boundary; (3) the grain boundary obstacle caused by the slip discontinuity is the main factor causing the large amount of residual shear stress, while the GND and back stress have little influence on the extent of shear stress relaxation.
Keywords: crystal plasticity;shock loading;grain boundary effect;grain boundary obstacle;geometrically necessary dislocations;back-stress
材料在高壓、高應(yīng)變率下的宏微觀行為與準(zhǔn)靜態(tài)情況存在巨大的差異。例如,材料在沖擊下會出現(xiàn)彈塑性雙波結(jié)構(gòu)[1]、彈性前驅(qū)波衰減[2]以及波陣面上的剪應(yīng)力松弛[3]等特征。這些動力學(xué)特征源于材料動態(tài)彈塑性變形機(jī)制與應(yīng)力波效應(yīng)的共同作用。目前,材料的沖擊動力學(xué)研究主要聚焦于單晶材料[4-6],缺乏對多晶材料的研究。然而,不同的晶界取向?qū)O大地影響材料的沖擊動力學(xué)響應(yīng)。沖擊加載后,多晶材料會出現(xiàn)波陣面上的變形不均勻、晶界附近殘留大量剪應(yīng)力等現(xiàn)象。
晶體塑性有限元方法(crystal plasticity finite element methods,CPFEM)是研究材料沖擊動力學(xué)響應(yīng)與變形機(jī)理的重要方法[7-10]。相比于低應(yīng)變率下的晶體塑性模型,沖擊動力學(xué)模型需要考慮位錯(cuò)的高速運(yùn)動(聲子拖曳機(jī)制)、位錯(cuò)的快速生成(形核)等塑性變形機(jī)制,以及大體積變形下的非線性彈性響應(yīng)與熱彈性耦合效應(yīng)等彈性行為[9,11-12]。對于多晶材料的沖擊動力學(xué)研究,現(xiàn)有的晶體塑性模型只考慮晶粒之間的取向,而忽略了局部的晶界強(qiáng)化、幾何必需位錯(cuò)(geometrically necessary dislocation,GND)以及背應(yīng)力(back-stress)等晶界效應(yīng)對本構(gòu)關(guān)系的影響。例如,Becker[13]的多晶銅動力學(xué)模擬只部分考慮了晶粒取向的影響,Vogler 等[14]關(guān)于晶粒取向?qū)αW铀俣扔绊懙难芯恳矝]有考慮微觀的晶界效應(yīng)。
在動態(tài)單晶晶體塑性模型的基礎(chǔ)上,本研究將考慮晶界與位錯(cuò)相互作用的微觀機(jī)理,建立一個(gè)包含晶界阻力、GND 以及背應(yīng)力的多晶晶體塑性模型。其中,單晶晶體塑性模型考慮了高應(yīng)變率下位錯(cuò)運(yùn)動的聲子拖曳機(jī)制和位錯(cuò)的非均勻形核機(jī)制,并采用超彈性熱力學(xué)框架描述材料的非線性彈性響應(yīng)和熱彈性耦合效應(yīng)??傮w研究思路為:首先,采用單晶晶體塑性模型模擬單晶鋁的平板沖擊實(shí)驗(yàn),進(jìn)而校正模型參數(shù);然后,將晶界阻力、GND 和背應(yīng)力等引入本構(gòu)關(guān)系,模擬研究“Voronoi 幾何模型”多晶鋁的沖擊動力學(xué)響應(yīng)。
定義中間構(gòu)型上(塑性變形完成但未發(fā)生彈性變形)的彈性格林-拉格朗日應(yīng)變?yōu)镋e,其表達(dá)式為
沖擊加載下材料的彈性變形往往伴隨著高壓和高應(yīng)變率,進(jìn)而產(chǎn)生大體積變形和急劇溫升。因此,沖擊載荷下材料的彈性本構(gòu)必須考慮非線性彈性行為和熱彈性耦合效應(yīng)。該超彈性本構(gòu)可以基于熱力學(xué)自由能推導(dǎo)得出[9,19]。其中,非線性彈性行為由高壓狀態(tài)方程(equation of state,EOS)以及隨壓強(qiáng)變化的彈性模量體現(xiàn),熱彈性耦合效應(yīng)通過熱膨脹和彈性變形引起的溫升體現(xiàn)。劉靜楠等[7]詳細(xì)介紹了這類超彈性本構(gòu)關(guān)系,本研究僅列出應(yīng)力和溫升的最終表達(dá)式。
以O(shè)rowan 流動方程為塑性框架,由滑移的位錯(cuò)數(shù)量和平均滑移速度直接衡量塑性滑移率,即
式中: ρα和vα分別為α 滑移系的位錯(cuò)密度和平均位錯(cuò)速度,b為柏氏矢量的大小。由于沖擊加載下所有位錯(cuò)幾乎均為可動位錯(cuò)[22],因此本研究對總位錯(cuò)密度與可動位錯(cuò)密度不作區(qū)分。
1.3.1 位錯(cuò)密度演化
除了低應(yīng)變率下主導(dǎo)位錯(cuò)密度演化的增殖與湮滅過程,材料在沖擊加載下需要額外考慮位錯(cuò)的形核過程
位錯(cuò)的非均勻形核指在晶體的應(yīng)力集中處發(fā)生的位錯(cuò)形核現(xiàn)象,所需的應(yīng)力遠(yuǎn)小于晶格的理論剪切強(qiáng)度。位錯(cuò)的非均勻形核模型通常建立在應(yīng)力形核點(diǎn)滿足一定分布規(guī)律的假設(shè)基礎(chǔ)上,本研究采用常用的線性分布的非均勻形核模型[10,12]
1.3.2 平均位錯(cuò)速度
位錯(cuò)在沖擊加載下的高速運(yùn)動由聲子拖曳效應(yīng)主導(dǎo)??紤]到?jīng)_擊加載下松弛過程的塑性應(yīng)變率顯著下降,本研究采用統(tǒng)一的位錯(cuò)滑移模型(結(jié)合熱激活機(jī)制和聲子拖曳機(jī)制)來描述不同應(yīng)力區(qū)間內(nèi)的位錯(cuò)運(yùn)動規(guī)律[9,23],該模型將位錯(cuò)在障礙物之間的平均運(yùn)動時(shí)間分成跨越障礙物所需的等待時(shí)間以及在障礙物之間的自由滑移時(shí)間兩部分
選取兩組單晶鋁的平板撞擊實(shí)驗(yàn)來驗(yàn)證單晶模型在沖擊加載下的適用性。這兩組實(shí)驗(yàn)分別為低壓下不同溫度的平板撞擊實(shí)驗(yàn)(彈塑性雙波)和高壓下不同壓力的平板撞擊實(shí)驗(yàn)(單波)。兩組平板撞擊實(shí)驗(yàn)的示意圖如圖1 所示。第1 組實(shí)驗(yàn)[26]中,鋁飛片以650 m/s 的速度(v0)直接撞擊[001]取向的單晶鋁樣品,通過VISAR(velocity interferometer system for any reflector)測得樣品的自由面速度;第2 組實(shí)驗(yàn)[27]中,銅飛片以2 900~3 636 m/s 的速度撞擊鋁基板,通過VISAR 測得[001]單晶鋁樣品與LiF 窗口的界面粒子速度。采用有限元軟件ABAQUS/Explicit 模擬平板撞擊實(shí)驗(yàn)的具體過程,其中:單晶鋁樣品采用編寫在VUMAT 中的晶體塑性模型描述,飛片、基板和窗口材料均采用ABAQUS 提供的本構(gòu)模型描述,包括Grüneisen 狀態(tài)方程和理想塑性本構(gòu)(方程及其參數(shù)[28-31]見附錄A)。為了提升計(jì)算效率,在有限元模擬中采用沿沖擊方向一維分布的六面體單元(單元尺寸為5 μm,通過網(wǎng)格收斂性驗(yàn)證),并對樣品中與沖擊方向垂直的兩個(gè)方向施加周期性邊界條件,從而消除邊界效應(yīng)(以小單元來近似大型系統(tǒng))。
圖1 平板撞擊實(shí)驗(yàn)示意圖Fig. 1 Illustration of the plate impact experiments
單晶晶體塑性模型中,與超彈性本構(gòu)和位錯(cuò)性質(zhì)有關(guān)的物理參數(shù)采用文獻(xiàn)[10-11, 32-35]中的數(shù)值,其他塑性參數(shù)采用遺傳算法(實(shí)現(xiàn)過程見附錄B)擬合平板撞擊實(shí)驗(yàn)獲得,參數(shù)值及其來源詳見表1和表2。根據(jù)上述數(shù)值方法和模型參數(shù),得到波剖面的模擬與實(shí)驗(yàn)對比結(jié)果,如圖2 所示。首先,單晶鋁在低壓(約4.5 GPa)沖擊下的自由面速度曲線呈現(xiàn)以下特征:(1) 彈塑性雙波結(jié)構(gòu),(2) Hugoniot彈性極限(Hugoniot elastic limit,HEL)隨溫度升高而增大的熱硬化現(xiàn)象。彈塑性雙波結(jié)構(gòu)是低壓下彈性波速大于塑性波速造成的;而沖擊加載下的熱硬化則源于位錯(cuò)高速運(yùn)動的聲子拖曳效應(yīng)(隨著溫度的升高,晶格振動越來越劇烈,聲子拖曳效應(yīng)對位錯(cuò)運(yùn)動的阻礙作用增強(qiáng))。CPFEM 模擬得到的自由面速度與實(shí)驗(yàn)結(jié)果呈現(xiàn)出較好的一致性,從3 條曲線極高的平均擬合優(yōu)度(R2= 0.934,擬合優(yōu)度R2是用于量化評估模擬與實(shí)驗(yàn)符合程度的常見指標(biāo))可以體現(xiàn)。其次,單晶鋁在高壓沖擊下(41~70 GPa)的樣品-窗口界面粒子速度呈現(xiàn)出單波特征,這是由于塑性波速隨著壓力的增大而增大,并最終導(dǎo)致高壓下的塑性波增加到與彈性波速一致并重合形成單一的沖擊波。高壓下的CPFEM 模擬結(jié)果與實(shí)驗(yàn)結(jié)果也符合得較好,3 條曲線的平均擬合優(yōu)度為0.986。上述結(jié)果表明,單晶晶體塑性模型能夠較好地描述單晶鋁在沖擊加載下的力學(xué)響應(yīng),為下一步模擬研究沖擊加載下多晶鋁的動力學(xué)行為提供了模型和參數(shù)基礎(chǔ)。
表1 單晶鋁的晶體塑性本構(gòu)參數(shù)Table 1 Parameters of crystal plasticity model for single-crystal Al
表2 單晶鋁的超彈性本構(gòu)參數(shù)Table 2 Parameters of hyper-elastic model for single-crystal Al
圖2 單晶鋁平板撞擊模擬與實(shí)驗(yàn)結(jié)果對比Fig. 2 Comparison of plate impact simulations and experiments for single-crystal Al
多晶材料具有取向不同的晶粒,這些晶粒之間的晶界對位錯(cuò)滑移起到阻礙作用,進(jìn)而顯著影響材料的力學(xué)性能。目前,在晶體塑性模型中考慮晶界效應(yīng)的方法主要有3 種:(1) 在位錯(cuò)的滑移阻力中引入與晶粒尺寸有關(guān)的Hall-Petch 公式(σs= σ0+kd-1/2,即多晶材料的屈服強(qiáng)度σs與晶粒平均直徑d之間的關(guān)系, σ0和k分別為初始強(qiáng)度和強(qiáng)化系數(shù)),將晶界的強(qiáng)化效果平均分配到晶粒內(nèi)部[36];(2) 單獨(dú)建立晶界模型,并唯象地賦予強(qiáng)度更大的材料屬性[37];(3) 基于位錯(cuò)在晶界附近的塞積效應(yīng),將塑性變形不均勻引起的GND、背應(yīng)力以及晶界處滑移不連續(xù)導(dǎo)致的晶界阻力引入本構(gòu)模型[38]。由于前兩種簡化處理方法要么基于平均化的思想,要么缺乏物理基礎(chǔ),忽略了位錯(cuò)的塞積效應(yīng),因此本研究采用第3 種方法建立多晶晶體塑性模型。
由于晶界兩側(cè)的晶粒取向不同,位錯(cuò)在穿過晶界過程中因滑移的不連續(xù)性而受到額外的阻力,阻力大小取決于兩側(cè)晶粒的取向差和滑移系的幾何關(guān)系。為了同時(shí)實(shí)現(xiàn)位錯(cuò)穿越晶界過程中的不同機(jī)制(完全穿過、部分穿過和不能穿過),采用廣泛使用的SWC 準(zhǔn)則[39]描述位錯(cuò)穿越晶界過程以及位錯(cuò)在晶界處滑移受到的阻力
晶界附近通常存在較大的塑性變形梯度和GND。GND 是為了在幾何上適應(yīng)由非均勻塑性變形產(chǎn)生的晶格曲率所必須存在的同號位錯(cuò)。晶體塑性模型中計(jì)算GND 的方法主要有兩種:(1) Shen 等[38]提出的超位錯(cuò)模型,以每個(gè)獨(dú)立的滑移系為對象,根據(jù)相鄰單元之間的塑性變形梯度,將一個(gè)單元內(nèi)的所有GND 簡化為一個(gè)在單元中心點(diǎn)上的超級刃位錯(cuò);(2) 先計(jì)算塑性變形梯度的旋度即Nye 張量,再根據(jù)Nye 張量計(jì)算每個(gè)單元每個(gè)滑移系的GND[23,41]。本質(zhì)上看,兩種方法都是基于材料內(nèi)部的塑性變形梯度來計(jì)算GND,而且造成的晶界強(qiáng)化效應(yīng)也被證明非常接近[41]。考慮到第1 種方法在物理圖像上更簡單清晰,模型的收斂性也更強(qiáng)[41],為此本研究采用超位錯(cuò)模型來衡量晶體中的GND 及其產(chǎn)生的背應(yīng)力。如圖3 所示,單元內(nèi)的GND( ρi)取決于相鄰單元之間的塑性變形梯度
圖3 塑性變形梯度與GNDFig. 3 Plastic strain gradient and GND
此外,背應(yīng)力可由相鄰兩個(gè)單元超位錯(cuò)的彈性應(yīng)力場求得
式中: ν=0.3 為泊松比,xj為相鄰單元中心之間的距離。
上述晶界阻力、GND、背應(yīng)力與晶體塑性模型的結(jié)合主要體現(xiàn)在兩方面:(1) 計(jì)算位錯(cuò)滑移阻力時(shí),需要在總位錯(cuò)密度中加入GND;(2) 對于滑移系上的等效分切應(yīng)力,需要考慮晶界阻力和背應(yīng)力的影響。首先,考慮GND 的滑移系阻力為
本研究采用ABAQUS 的兩個(gè)子程序(VUMAT 和VUFIELD)共同實(shí)現(xiàn)沖擊加載下的多晶晶體塑性模型,具體關(guān)系如圖4(a)所示。一方面,VUFIELD 負(fù)責(zé)在全局上根據(jù)相鄰單元的幾何關(guān)系和塑性變形梯度計(jì)算晶界阻力、GND 以及背應(yīng)力,并通過全局變量傳遞給VUMAT。另一方面,VUMAT 在更新總位錯(cuò)密度和等效分切應(yīng)力后負(fù)責(zé)計(jì)算各個(gè)積分點(diǎn)上的材料本構(gòu)關(guān)系,并最終向VUFIELD 傳遞新的塑性滑移量。
多晶的Voronoi 幾何模型采用開源多晶軟件Neper 生成,導(dǎo)出的幾何模型可進(jìn)一步在ABAQUS 中施加邊界條件,每個(gè)晶粒采用隨機(jī)取向(歐拉角)。本研究將模擬多晶鋁的沖擊實(shí)驗(yàn),樣品尺寸為0.5 mm×0.5 mm×0.5 mm,共包含100 個(gè)晶粒238 328 個(gè)單元,同時(shí)在樣品左側(cè)(x方向)施加壓力邊界條件(壓力p=5 GPa,上升時(shí)間t=0.005 μs),對與沖擊方向垂直的4 個(gè)面分別施加沿面法向的固定位移約束(一維應(yīng)變),有限元模型如圖4(b)所示。最終根據(jù)建立的多晶模型和數(shù)值方法,可以模擬得到?jīng)_擊加載下多晶鋁的力學(xué)響應(yīng)。
圖4 數(shù)值實(shí)現(xiàn)與有限元模型Fig. 4 Numerical implementation and finite element model
沖擊波后多晶鋁的Mises 應(yīng)力分布體現(xiàn)了剪應(yīng)力的松弛程度。如圖5(a)所示,在考慮晶界效應(yīng)的模擬結(jié)果中,晶界處與晶粒內(nèi)部的剪應(yīng)力松弛呈現(xiàn)出不同特征:晶界單元有大量的殘余剪應(yīng)力,而晶內(nèi)單元的剪應(yīng)力正常松弛并趨于零。圖5(b)顯示了不考慮晶界效應(yīng)的模擬結(jié)果(忽略晶界阻力、GND 與背應(yīng)力),可以看出:沖擊波后的晶界單元不再有大量的殘余剪應(yīng)力,材料整體上處于均勻松弛狀態(tài)。
圖5 沖擊加載下多晶鋁的Mises 剪應(yīng)力分布Fig. 5 Distribution of Mises stress of polycrystalline aluminum under shock loading
圖6 沖擊加載下多晶鋁的累積塑性滑移量分布Fig. 6 Distribution of accumulated slip of polycrystalline aluminum under shock loading
晶界附近塑性變形的不均勻性會進(jìn)一步造成GND 與背應(yīng)力的累積。如圖7(a)和圖7(b)所示(以晶界為背景),不同滑移系上的GND 分布存在區(qū)別,但整體上均分布在晶界附近。由背應(yīng)力的分布(圖7(c)和圖7(d))可知,相同滑移系下背應(yīng)力與GND 在晶界附近的分布基本一致。需要注意的是,在晶界附近累積的GND 和背應(yīng)力會進(jìn)一步阻礙位錯(cuò)的滑移,同時(shí)影響沖擊波后剪應(yīng)力的松弛程度。
圖7 沖擊加載下多晶鋁的GND 和背應(yīng)力分布Fig. 7 Distribution of GND and back-stress of polycrystalline aluminum under shock loading
為了進(jìn)一步在樣品尺度上量化多晶鋁在沖擊波后殘余剪應(yīng)力的大小,以及兩種晶界效應(yīng)(晶界阻力與GND/背應(yīng)力)對剪應(yīng)力松弛的影響程度,模擬并量化統(tǒng)計(jì)了不同模型(晶界效應(yīng)不同)下沖擊波后所有單元的平均Mises 剪應(yīng)力。如圖8 所示,考慮完整晶界效應(yīng)(包括晶界阻力、GND 和背應(yīng)力)時(shí)的平均殘余剪應(yīng)力為183 MPa,遠(yuǎn)高于不考慮晶界效應(yīng)時(shí)的平均殘余剪應(yīng)力(21 MPa)。其次,當(dāng)僅考慮晶界阻力時(shí)(忽略GND 和背應(yīng)力),材料仍然具有較大的平均殘余剪應(yīng)力(158 MPa),且與考慮完整晶界效應(yīng)的結(jié)果接近。上述結(jié)果表明:由取向差引起的晶界阻力是造成沖擊波后大量殘余剪應(yīng)力的主要因素,而GND 和背應(yīng)力對剪應(yīng)力松弛程度的影響較小。
圖8 沖擊加載下不同模型的平均殘余剪應(yīng)力Fig. 8 Average residual shear stress for different models under shock loading
在單晶晶體塑性模型的基礎(chǔ)上,通過分析晶界與位錯(cuò)相互作用的微觀機(jī)理,建立了一個(gè)包含晶界阻力、GND 以及背應(yīng)力的多晶晶體塑性模型,為模擬分析多晶材料在沖擊下的變形提供了有效的工具。采用該模型模擬研究了基于Voronoi 幾何模型的多晶鋁在沖擊加載下的力學(xué)響應(yīng),得到以下結(jié)論。
(1) 沖擊波后晶界單元與晶內(nèi)單元的剪應(yīng)力松弛呈現(xiàn)不同特征:晶界單元存在極高的殘余剪應(yīng)力,而晶內(nèi)單元的剪應(yīng)力正常快速松弛并趨于零。這種剪應(yīng)力松弛程度的區(qū)別取決于塑性變形的難易程度,由于晶界單元存在額外的晶界阻力,位錯(cuò)難以滑移,從而無法釋放剪應(yīng)力。
(2) 晶界附近存在較大的塑性變形梯度,造成沿晶界分布的GND 和背應(yīng)力。晶界阻力是影響剪應(yīng)力松弛程度的主導(dǎo)因素,GND 和背應(yīng)力對沖擊波后的剪應(yīng)力松弛影響較小。
(3) 在樣品尺度上量化了沖擊波后的剪應(yīng)力松弛程度,多晶的晶界效應(yīng)使得沖擊波后的平均殘余剪應(yīng)力為183 MPa,遠(yuǎn)高于不考慮晶界效應(yīng)時(shí)的平均殘余剪應(yīng)力(21 MPa)。
附錄 A ABAQUS 提供的本構(gòu)模型與材料參數(shù)
單晶沖擊模型中的飛片、基板和窗口材料都采用ABAQUS 提供的本構(gòu)模型描述,包括Grüneisen狀態(tài)方程(影響樣品Hugoniot 態(tài)的壓力)和理想塑性本構(gòu)(對結(jié)果影響極?。?。其中Grüneisen 狀態(tài)方程為
有限元模型中Al、Cu 以及LiF 的模型參數(shù)見表A1。
表A1 ABAQUS 本構(gòu)模型參數(shù)Table A1 Parameters of constitutive model provided by ABAQUS
附錄 B 遺傳算法
遺傳算法是一種基于自然選擇原理和遺傳機(jī)制的全局搜索優(yōu)化算法,本質(zhì)上是通過模擬自然界中生物群體的進(jìn)化機(jī)制,以優(yōu)勝劣汰(適者生存)的原則進(jìn)行參數(shù)的最優(yōu)化。它主要包括以下步驟:(1) 初始群體的產(chǎn)生;(2) 求每一個(gè)體的適應(yīng)度;(3) 根據(jù)適者生存的原則選擇優(yōu)良個(gè)體;(4) 被擇優(yōu)選出的個(gè)體進(jìn)行配對,通過隨機(jī)交叉、隨機(jī)變異的方法改變某些個(gè)體的基因并生成下一代群體;(5) 按此方法不斷使群體逐代進(jìn)化,直到滿足進(jìn)化終止條件。
對于晶體塑性模型而言,每一組參數(shù)都是一個(gè)生物個(gè)體,參數(shù)的數(shù)值代表個(gè)體基因。求適應(yīng)度的過程相當(dāng)于根據(jù)模型參數(shù)調(diào)用ABAQUS 和VUMAT 計(jì)算模擬結(jié)果并與實(shí)驗(yàn)數(shù)據(jù)對比求解殘差的過程,殘差越小,則數(shù)據(jù)結(jié)果越接近,代表個(gè)體的適應(yīng)度越大。群體不斷進(jìn)化的過程就是在給定區(qū)間內(nèi)不斷搜索最優(yōu)參數(shù)的過程,最終可以得到與實(shí)驗(yàn)曲線最吻合的參數(shù)組合,遺傳算法的具體流程如圖B1所示。
圖B1 遺傳算法流程Fig. B1 Procedure for genetic algorithms