高溫、高壓、高應(yīng)變速率動(dòng)態(tài)過程晶體塑性有限元理論模型及其應(yīng)用

劉靜楠，葉常青，劉桂森，沈 耀

（上海交通大學(xué)材料科學(xué)與工程學(xué)院，上海 200240）

沖擊載荷下的材料動(dòng)力學(xué)響應(yīng)、微觀結(jié)構(gòu)演化以及材料損傷與破壞一直是軍事、材料科學(xué)及工程技術(shù)等領(lǐng)域關(guān)注的重點(diǎn)。動(dòng)態(tài)沖擊試驗(yàn)通常只能測(cè)定材料的初始微觀結(jié)構(gòu)、變形過程的力學(xué)響應(yīng)以及變形后的宏觀性能和部分微觀結(jié)構(gòu)，而無法直接測(cè)量沖擊下材料微觀結(jié)構(gòu)的演化過程以及初始微觀結(jié)構(gòu)對(duì)材料宏觀性能的影響過程。此外，動(dòng)態(tài)沖擊下被廣泛使用的Bodner-Parton（BP）模型[1]、Johnson-Cook（JC）模型[2]、Zerilli-Armstrong（ZA）模型[3]以及Steinberg-Guinan（SG）[4]模型等宏觀唯象材料本構(gòu)模型并不包含塑性變形的微觀機(jī)制，而只限于擬合實(shí)驗(yàn)數(shù)據(jù)，只在一定適用范圍內(nèi)描述實(shí)驗(yàn)中的應(yīng)力-應(yīng)變關(guān)系以及應(yīng)變率效應(yīng)，而無法描述沖擊試驗(yàn)中的應(yīng)力松弛[5]、彈性前驅(qū)波衰減[6]等異常現(xiàn)象。因此，考慮微觀塑性滑移變形機(jī)制的動(dòng)態(tài)晶體塑性模型在近年來得到了快速發(fā)展，并展現(xiàn)出對(duì)材料動(dòng)力學(xué)響應(yīng)與微觀結(jié)構(gòu)演化的描述能力。

材料在高壓、高應(yīng)變率下的宏觀動(dòng)力學(xué)特性和微觀變形機(jī)理均與準(zhǔn)靜態(tài)下存在巨大的差異，適用于沖擊加載的動(dòng)態(tài)晶體塑性模型比準(zhǔn)靜態(tài)模型更為復(fù)雜，需考慮更多影響變形的宏觀因素與微觀機(jī)制。宏觀因素包括影響變形的靜水壓強(qiáng)、溫度等，微觀機(jī)制包括聲子拖曳、位錯(cuò)形核等位錯(cuò)行為，以及相變、孿生等變形模式。動(dòng)態(tài)晶體塑性模型是在準(zhǔn)靜態(tài)模型的基礎(chǔ)上建立的，繼承了準(zhǔn)靜態(tài)模型的彈塑性變形框架，考慮了晶體中塑性滑移的各向異性，也可以考慮孿生與相變等變形機(jī)制，能夠反映變形過程中的微觀結(jié)構(gòu)演化。但是，動(dòng)態(tài)晶體塑性模型通常需要額外考慮材料的非線性彈性響應(yīng)、大變形引起的顯著溫升效應(yīng)以及位錯(cuò)在高速滑移下受到的聲子阻礙[7]等因素，甚至在沖擊載荷下材料可能發(fā)生的位錯(cuò)均勻形核[8-9]等機(jī)制。此外，材料在動(dòng)態(tài)加載下相較于準(zhǔn)靜態(tài)加載更容易發(fā)生孿生[10-12]、相變[13-16]、斷裂損傷[17-18]等變形機(jī)制。Roters 等[19]詳述了準(zhǔn)靜態(tài)下晶體塑性模型的理論與應(yīng)用。鄭松林[20]也對(duì)晶體塑性模型在準(zhǔn)靜態(tài)與動(dòng)態(tài)響應(yīng)中的理論與應(yīng)用進(jìn)行了較為詳細(xì)的介紹。本文側(cè)重于動(dòng)態(tài)沖擊加載的晶體塑性模型，關(guān)于準(zhǔn)靜態(tài)的晶體塑性模型理論與應(yīng)用可參閱Roters 等[19]和鄭松林[20]的評(píng)論。

本文聚焦于動(dòng)態(tài)沖擊加載下的晶體塑性模型，詳細(xì)綜述動(dòng)態(tài)晶體塑性有限元理論框架的運(yùn)動(dòng)學(xué)、彈性本構(gòu)和塑性本構(gòu)等組成部分，以及模型在孿生、相變、斷裂損傷等不同變形機(jī)制中的擴(kuò)展。動(dòng)態(tài)晶體塑性理論模型框架的要點(diǎn)主要有沖擊加載過程中溫度的演化規(guī)律、熱彈性耦合效應(yīng)、狀態(tài)方程在超彈性本構(gòu)中的引入方法，以及聲子拖曳效應(yīng)對(duì)塑性滑移的阻礙作用等。高應(yīng)變率沖擊加載通常會(huì)引起材料內(nèi)部的劇烈溫升，因此動(dòng)態(tài)晶體塑性模型在沖擊加載下需要考慮溫升帶來的熱彈性耦合效應(yīng)。由于材料在沖擊載荷下的體積變形很大，故通常采用高壓狀態(tài)方程(Equation of state，EOS)來描述材料的非線性彈性容變關(guān)系。對(duì)于高應(yīng)變率下的塑性變形過程，材料晶格的聲子振動(dòng)會(huì)對(duì)位錯(cuò)的高速運(yùn)動(dòng)產(chǎn)生強(qiáng)烈的阻礙作用，即聲子拖曳，且這種阻礙作用隨著位錯(cuò)滑移速度的增加而強(qiáng)化，因此在動(dòng)態(tài)沖擊晶體塑性模型中必須引入聲子拖曳機(jī)制來描述材料的應(yīng)變率硬化過程。在動(dòng)態(tài)沖擊加載下，位錯(cuò)滑移往往不再是材料唯一的塑性變形方式，材料還可能通過相變、孿生、斷裂損傷等變形方式釋放體系的自由能。因此本文在基于位錯(cuò)滑移的動(dòng)態(tài)晶體塑性模型框架之外還介紹了將相變、孿生、損傷等變形機(jī)制引入晶體塑性模型的方法，包括動(dòng)態(tài)相變的驅(qū)動(dòng)力與演化規(guī)律、孿晶的演化規(guī)律與硬化模型、層裂與絕熱剪切帶的形成及其對(duì)材料本構(gòu)的影響。

1 變形運(yùn)動(dòng)學(xué)

在連續(xù)介質(zhì)力學(xué)理論中，物質(zhì)占據(jù)空間幾何區(qū)域所構(gòu)成的圖形，稱為物質(zhì)構(gòu)型。在變形或剛體運(yùn)動(dòng)過程中，物質(zhì)構(gòu)型將發(fā)生改變。為了對(duì)晶體材料的變形過程進(jìn)行描述，需要建立晶體運(yùn)動(dòng)學(xué)框架。初始時(shí)刻未經(jīng)變形的晶體材料占據(jù)的空間幾何圖形稱為初始構(gòu)型。采用拉格朗日坐標(biāo)(X)來描述初始構(gòu)型中任意物質(zhì)點(diǎn)的位置。晶體材料在變形過程中的瞬時(shí)空間幾何圖形稱為現(xiàn)時(shí)構(gòu)型。對(duì)于初始構(gòu)型中位于坐標(biāo)X的物質(zhì)點(diǎn)，在現(xiàn)時(shí)構(gòu)型中采用歐拉坐標(biāo)(x)描述其瞬時(shí)空間位置。顯然，物質(zhì)點(diǎn)在初始和現(xiàn)時(shí)構(gòu)型之間具有x=x(X,t)的位置對(duì)應(yīng)關(guān)系。描述材料變形程度的基本物理量通常選擇變形梯度F= ?x/?X，其物理含義是現(xiàn)時(shí)構(gòu)型相對(duì)于初始構(gòu)型的梯度張量。物質(zhì)點(diǎn)的運(yùn)動(dòng)速度通常定義在拉格朗日坐標(biāo)系中v(X,t) = ?x(X,t)/?t，也可采用歐拉坐標(biāo)表示為v(x,t) =v[X(x,t),t]，通常定義速度梯度L= ?v/?x描述物質(zhì)點(diǎn)速度在現(xiàn)時(shí)構(gòu)型中的梯度。

1.1 彈塑性變形分解

經(jīng)典的運(yùn)動(dòng)學(xué)構(gòu)型如圖1 所示。晶體材料由初始構(gòu)型經(jīng)變形梯度(F)變換至現(xiàn)時(shí)構(gòu)型，如圖1中的黑色箭頭所示。Hill 等[21]、Asaro 等[22]建立中間構(gòu)型Ⅰ，將彈性和塑性變形進(jìn)行乘和分解

圖1 經(jīng)典的晶體運(yùn)動(dòng)學(xué)構(gòu)型Fig. 1 Classical configurations of crystal kinematics

式中：Fe、Fp分別為彈性和塑性變形梯度。材料由初始構(gòu)型經(jīng)塑性變形(Fp)變換至中間構(gòu)型Ⅰ，再經(jīng)彈性變形(Fe)變換至現(xiàn)時(shí)構(gòu)型，如圖1 中的紅色箭頭所示。彈性變形包括剛體旋轉(zhuǎn)和晶格畸變，塑性變形可由位錯(cuò)滑移、相變、孿生等引起。將現(xiàn)時(shí)構(gòu)型中速度梯度張量進(jìn)行彈、塑性分解

式中：Re為旋轉(zhuǎn)張量，Ue為右伸長(zhǎng)張量。材料由中間構(gòu)型Ⅰ，經(jīng)彈性拉伸(Ue)變換至中間構(gòu)型Ⅱ，再經(jīng)晶格旋轉(zhuǎn)(Re)變換至現(xiàn)時(shí)構(gòu)型，如圖1 中的藍(lán)色箭頭所示。

圖 2 引入熱膨脹構(gòu)型的變形梯度分解F = FeFθFp[19]Fig. 2 Decomposition of deformation gradient considered thermally-expanded configuration F=FeFθFp[19]

1.2 熱膨脹構(gòu)型

材料在動(dòng)態(tài)沖擊加載下往往存在顯著的溫升現(xiàn)象，而溫升引起的熱膨脹會(huì)導(dǎo)致材料內(nèi)部應(yīng)力、應(yīng)變的變化，因此如何考慮熱膨脹效應(yīng)的影響是動(dòng)態(tài)晶體塑性模型的一個(gè)重要問題。

在晶體塑性模型框架中，通常有兩種考慮熱膨脹效應(yīng)的方法。第1 種方法是在彈塑性兩構(gòu)型的基礎(chǔ)上額外引入熱構(gòu)型，將熱變形與彈性變形解耦。第2 種方法不額外引入熱構(gòu)型，而在現(xiàn)有框架下基于熱力學(xué)自由能進(jìn)行推導(dǎo)，在求解超彈性本構(gòu)過程中考慮熱膨脹效應(yīng)(詳見2.2.2 節(jié))。本節(jié)主要介紹第1 種引入熱構(gòu)型的方法。Clayton[23]、Vogler 等[24]、Hansen 等[25]、De 等[26]以及Shahba 等[27]在構(gòu)型分解時(shí)加入了熱構(gòu)型Fθ，如圖2 所示，將變形梯度分解為塑性變形、熱變形以及彈性變形3 部分

式中：Fe、Fθ、Fp分別表示彈性變形梯度、熱變形梯度和塑性變形梯度。

相應(yīng)地，現(xiàn)時(shí)構(gòu)型中的速度梯度可分解為

式中：Lθ、分別為現(xiàn)時(shí)構(gòu)型和熱構(gòu)型中的熱變形速度梯度，是溫度變化率 (T˙)的函數(shù)

式中：α為熱膨脹系數(shù)張量。

2 超彈性本構(gòu)框架

通常在晶體塑性模型中有超彈性和次彈性兩種描述材料彈性行為的本構(gòu)關(guān)系。超彈性本構(gòu)是指材料存在一個(gè)與應(yīng)變相關(guān)的能量函數(shù)，其應(yīng)力與應(yīng)變?nèi)恐g的關(guān)系由能量函數(shù)決定，具有做功獨(dú)立于加載路徑的特點(diǎn)。次彈性本構(gòu)描述應(yīng)力率與應(yīng)變率之間的關(guān)系，由于大變形下彈性模量是隨應(yīng)變?nèi)孔兓模鋵?shí)際應(yīng)用比較復(fù)雜，通常只用于小彈性變形情況[28](這種情況下彈性模量可視為常數(shù))?？紤]到實(shí)驗(yàn)中通常采用全量形式的高壓固體狀態(tài)方程描述材料的非線性容變關(guān)系，因此研究材料在沖擊載荷下的大彈性變形行為時(shí)，往往采用超彈性本構(gòu)來描述。

材料在沖擊載荷下的超彈性本構(gòu)關(guān)系主要包含基于熱力學(xué)自由能的超彈性應(yīng)力表達(dá)式以及考慮熱彈性耦合效應(yīng)的溫升表達(dá)式兩方面。其中，特別需要考慮熱膨脹對(duì)應(yīng)力的影響、彈性模量隨壓強(qiáng)和溫度的變化，以及非線性彈性容變關(guān)系(固體狀態(tài)方程)等。

2.1 超彈性本構(gòu)的熱力學(xué)框架

3.2 基于位錯(cuò)行為的動(dòng)態(tài)晶體塑性框架

唯象的晶體塑性本構(gòu)通常存在適用范圍的局限性，只能在一定應(yīng)變率區(qū)間內(nèi)描述材料的塑性變形規(guī)律。因此為了滿足不同應(yīng)變率的加載條件，并同時(shí)描述各微觀結(jié)構(gòu)的演化，通常采用Orowan 方程作為流動(dòng)法則的模型框架

Orowan 方程模型框架主要包含可動(dòng)位錯(cuò)密度的演化以及平均位錯(cuò)速度與各內(nèi)變量之間的關(guān)系兩部分內(nèi)容?？蓜?dòng)位錯(cuò)密度演化部分需要考慮位錯(cuò)的形核、增殖、纏結(jié)、湮滅等微觀過程。位錯(cuò)平均滑移速度部分可由位錯(cuò)滑移障礙物之間的平均距離除以在障礙物之間運(yùn)動(dòng)的平均時(shí)間得到。而位錯(cuò)在障礙物之間的平均運(yùn)動(dòng)時(shí)間可分成跨越障礙物的等待時(shí)間以及在障礙物之間的滑移時(shí)間兩部分，位錯(cuò)的兩個(gè)運(yùn)動(dòng)過程可分別由Arrhenius 熱激活機(jī)制和黏性拖曳機(jī)制描述。

3.2.1 位錯(cuò)滑移速度

位錯(cuò)速度模型描述了位錯(cuò)滑移速度與滑移系分切應(yīng)力以及位錯(cuò)滑移阻力之間的關(guān)系，主要可分為3 類：(1)簡(jiǎn)單指數(shù)模型，即基于實(shí)驗(yàn)中觀察到的單根位錯(cuò)的運(yùn)動(dòng)規(guī)律，建立位錯(cuò)速度與分切應(yīng)力的簡(jiǎn)單指數(shù)關(guān)系；(2)單一機(jī)制模型，即聚焦于位錯(cuò)在低應(yīng)力區(qū)間內(nèi)或高應(yīng)力區(qū)間內(nèi)的滑移規(guī)律，只采用單一的熱激活模型或線性拖曳模型；(3)統(tǒng)一模型，即將位錯(cuò)的滑移過程分為跨越障礙物的等待過程以及在障礙物之間的滑移過程，統(tǒng)一了位錯(cuò)在不同應(yīng)力區(qū)間內(nèi)的滑移規(guī)律。

第1 類模型是基于單根位錯(cuò)實(shí)驗(yàn)結(jié)果的簡(jiǎn)單擬合，只描述了位錯(cuò)運(yùn)動(dòng)速度與應(yīng)力之間的關(guān)系，不僅忽略了位錯(cuò)網(wǎng)絡(luò)對(duì)滑移的阻礙作用，也沒有考慮溫度對(duì)滑移的影響。第2 類模型只描述位錯(cuò)在某一應(yīng)力區(qū)間與溫度范圍內(nèi)的單一規(guī)律，適用范圍較窄，缺乏對(duì)位錯(cuò)完整滑移過程的描述。第3 類模型描述了位錯(cuò)滑移的整個(gè)物理過程，適用于不同的應(yīng)力區(qū)間和溫度范圍。本節(jié)將簡(jiǎn)單介紹前兩類模型，重點(diǎn)介紹適用范圍最廣的第3 類模型。

（1）簡(jiǎn)單指數(shù)模型

Johnston 與Gilman[62]利用蝕坑法測(cè)量了LiF 單晶中單根螺位錯(cuò)與刃位錯(cuò)在不同應(yīng)力下的位錯(cuò)速度，根據(jù)位錯(cuò)運(yùn)動(dòng)速度與應(yīng)力關(guān)系的實(shí)驗(yàn)結(jié)果得到以下結(jié)論：其一，刃位錯(cuò)的滑移速度比螺位錯(cuò)更快，但兩者的差距在接近聲速時(shí)逐漸消失；其二，位錯(cuò)運(yùn)動(dòng)的最高速度不能超過聲速；其三，存在一個(gè)位錯(cuò)開始滑移的最低臨界應(yīng)力；其四，在低溫低應(yīng)力下，位錯(cuò)速度與溫度近似滿足負(fù)指數(shù)關(guān)系。這些結(jié)論已被證實(shí)在其他材料中同樣適用。

基于上述實(shí)驗(yàn)結(jié)果，Gilman[63]提出了一種簡(jiǎn)單的位錯(cuò)速度與應(yīng)力關(guān)系的指數(shù)模型

式中：D為特征阻力；v*為位錯(cuò)運(yùn)動(dòng)的極限速度，通常取材料的聲速(剪切波速)。

上述簡(jiǎn)單的指數(shù)模型滿足實(shí)驗(yàn)上位錯(cuò)速度的邊界條件：當(dāng)τ→0 時(shí)，v→0；當(dāng)τ→∞時(shí)，v→v*。但該模型不僅缺乏對(duì)位錯(cuò)運(yùn)動(dòng)機(jī)制的解釋，也沒有反映出位錯(cuò)在低應(yīng)力區(qū)間的轉(zhuǎn)變特征。

為了將模型應(yīng)用于更廣的應(yīng)力區(qū)間，Gilman[64]提出了在單指數(shù)模型的基礎(chǔ)上添加一項(xiàng)在低應(yīng)力區(qū)間內(nèi)的位錯(cuò)速度

此外，Johnson 與Barker[65]相繼提出了關(guān)于位錯(cuò)運(yùn)動(dòng)速度的指數(shù)模型，但這些模型均只是根據(jù)位錯(cuò)的運(yùn)動(dòng)趨勢(shì)提出的唯象模型，缺乏相應(yīng)的物理機(jī)制，且適用范圍較窄。

（2）單一機(jī)制模型

單一機(jī)制模型包括只適用于低應(yīng)力區(qū)間內(nèi)的熱激活機(jī)制以及只適用于高應(yīng)力區(qū)間的線性拖曳機(jī)制。Alankar 等[66]、Zhang 等[67]、Monnet 等[68]采用熱激活機(jī)制來描述位錯(cuò)在低應(yīng)力區(qū)間內(nèi)的滑移速度

式中：λ為障礙物之間的平均距離，vd為Debye 頻率， ΔG0為參考激活能。此外，Hansen 等[25]、Nguyen 等[69]、Grilli 等[70]均采用線性拖曳模型來描述位錯(cuò)的滑移速度

式中：B為拖曳系數(shù)。

單一機(jī)制模型描述了位錯(cuò)在某一應(yīng)力區(qū)間內(nèi)的滑移規(guī)律，適用于模擬材料在特定加載區(qū)間內(nèi)的變形過程。但實(shí)際上材料在不同時(shí)間或空間尺度上的應(yīng)變率往往存在很大差距，單一機(jī)制模型難以描述這種復(fù)雜的變形過程。

（3）統(tǒng)一物理模型

目前使用最為廣泛的模型基于Frost 與Ashby[71]首先提出的平均位錯(cuò)滑移速度模型框架。該框架將位錯(cuò)在障礙物之間的平均運(yùn)動(dòng)時(shí)間分成了跨越障礙物的等待時(shí)間以及在障礙物之間的滑移時(shí)間兩部分，可分別采用熱激活模型與聲子拖曳模型來描述

式中：λα表示位錯(cuò)滑移過程中障礙物的平均距離，為熱激活過程主導(dǎo)的位錯(cuò)脫釘所需要的等待時(shí)間，為聲子拖曳主導(dǎo)的位錯(cuò)在障礙物之間黏性滑移所需要的時(shí)間。

對(duì)于滑移較為困難的體心立方(bcc)晶體材料，Shahba 等[27]、Lim 等[72]采用螺位錯(cuò)滑移的Kink-Pair 機(jī)制描述上述位錯(cuò)模型的兩個(gè)滑移過程。Kink-Pair 滑移機(jī)制可分為兩個(gè)階段，如圖6[27]所示。在第1 階段，位錯(cuò)在障礙物(Peierls 能壘)前有一段等待時(shí)間，直到成功通過熱激活方式跨越障礙物。在第2 階段，位錯(cuò)分離成一對(duì)扭折并以黏性滑移的方式運(yùn)動(dòng)到下一個(gè)障礙物。對(duì)于面心立方(fcc)等晶體材料，也可根據(jù)與林位錯(cuò)的交互機(jī)制將位錯(cuò)滑移分為類似的等待過程與滑移過程兩部分。

圖6 螺位錯(cuò)滑移的Kink-pair 機(jī)制[27]Fig. 6 Illustration of screw dislocation motion via a Kink-pair mechanism[27]

式中：CT為纏結(jié)系數(shù)。

位錯(cuò)的湮滅指符號(hào)相反的位錯(cuò)在一定捕獲距離內(nèi)相遇所發(fā)生的消亡現(xiàn)象，決定了位錯(cuò)飽和密度的存在。Luscher 等[39]、Alankar 等[66]根據(jù)上述湮滅過程建立了簡(jiǎn)單的湮滅模型

式中：CA為湮滅系數(shù)，da為捕獲距離。

為了分析沖擊加載下位錯(cuò)各演化機(jī)制在不同壓力下的主導(dǎo)地位轉(zhuǎn)變，Lloyd 等[83]模擬了5～20 GPa加載下位錯(cuò)各個(gè)演化機(jī)制的占比情況，如圖8[83]所示。在中低壓(5～15 GPa)下，位錯(cuò)的增殖主要依賴于非均勻形核過程與增殖過程；當(dāng)壓力大到一定程度(大于20 GPa)時(shí)，材料發(fā)生劇烈的均勻形核過程，成為位錯(cuò)密度演化的主導(dǎo)因素。

實(shí)際上除了可動(dòng)位錯(cuò)和不可動(dòng)位錯(cuò)之外，位錯(cuò)密度還有許多其他區(qū)分方式。例如，Alankar 等[66,84]為了體現(xiàn)刃位錯(cuò)和螺位錯(cuò)滑移速率的差異，將位錯(cuò)密度分為螺位錯(cuò)ρs和刃位錯(cuò)ρe，并解釋了基面取向的單晶Ti 的硬化率轉(zhuǎn)變現(xiàn)象。此外，Keshavarz 等[85]、Liang 等[86]為了考慮幾何必須位錯(cuò)引起的晶界效應(yīng)，將位錯(cuò)區(qū)分為統(tǒng)計(jì)存儲(chǔ)位錯(cuò)ρssd和幾何必須位錯(cuò)ρGND，分析了在Al 雙晶、Ni 基合金等材料中晶界對(duì)材料的硬化作用，并與實(shí)驗(yàn)結(jié)果較好地吻合。

圖8 不同壓力加載下位錯(cuò)密度的演化機(jī)制[83]Fig. 8 Dislocation density evolution mechanisms under different loading pressure[83]

4 動(dòng)態(tài)相變

在高溫、高壓、高應(yīng)變率動(dòng)態(tài)沖擊條件下，位錯(cuò)滑移往往不再是材料唯一的塑性變形方式，例如Fe、Ti 和RDX 等晶體材料也可通過相變與孿晶等方式實(shí)現(xiàn)變形。固體相變按原子的遷移方式可分為擴(kuò)散型相變和位移型(切變)相變。而在高應(yīng)變率動(dòng)態(tài)變形下，材料發(fā)生的往往是切變型馬氏體相變，整個(gè)相變過程無原子擴(kuò)散，通過原子協(xié)同進(jìn)行短距離遷移，并以切變的形式完成相變。本節(jié)重點(diǎn)介紹如何將馬氏體相變機(jī)制引入晶體塑性有限元框架，以及在動(dòng)態(tài)沖擊與靜態(tài)加載下相變驅(qū)動(dòng)力與演化準(zhǔn)則的差異。

當(dāng)然除了在晶體塑性有限元框架中直接引入馬氏體相變機(jī)制之外，也可采用相場(chǎng)與動(dòng)態(tài)晶體塑性有限元結(jié)合的方法來模擬變形過程中的相變過程。不同于晶體塑性框架中的局部相變模型，相場(chǎng)法在自由能的基礎(chǔ)上考慮了相變量在空間的分布及其所引起的梯度能。相場(chǎng)與動(dòng)態(tài)晶體塑性模型相結(jié)合的方法已被廣泛應(yīng)用于模擬準(zhǔn)靜態(tài)下材料的相變、孿生與再結(jié)晶等變形過程[87-91]，該方法有望發(fā)展至模擬動(dòng)態(tài)沖擊下的相變、孿生等塑性變形過程。

4.1 馬氏體相變運(yùn)動(dòng)學(xué)

關(guān)于馬氏體相變介觀模型的研究由來已久[92-96]。Thamburaja 和Anand[97]最初將馬氏體相變機(jī)制引入晶體塑性框架中，以研究Ni-Ti 形狀記憶合金的宏觀力學(xué)響應(yīng)，模型中總變形分為彈性變形和相變變形，忽略了位錯(cuò)滑移產(chǎn)生的塑性變形。在此基礎(chǔ)上發(fā)展起來的更為廣泛使用的模型考慮了塑性滑移，將變形梯度[13,98-99]分解為

式中：Fe、Fp、Ftr分別表示彈性變形、塑性變形以及相變對(duì)總變形梯度的貢獻(xiàn)。

塑性變形包括母相與新相中的塑性滑移總和

式中：vp=1?vN表示母相的體積分?jǐn)?shù)，

t=1為相變后新相的總體積分?jǐn)?shù)，Nt為變體的數(shù)量，Nα和Nβ分別代表母相和新相的滑移系，和分別表示母相中的滑移方向和滑移面法向，和分別為新相中的滑移方向和滑移面法向。當(dāng)相變產(chǎn)生的新相的位錯(cuò)滑移能力較弱時(shí)，往往只計(jì)及母相中的塑性滑移[100]。

馬氏體相變過程的變形梯度為

式中：bt和mt分別為相變過程的切變向量和相變慣習(xí)面的法向。

4.2 相變驅(qū)動(dòng)力與演化準(zhǔn)則

相變模型的關(guān)鍵在于如何描述相變過程的驅(qū)動(dòng)力以及相變的演化規(guī)律。本節(jié)將重點(diǎn)介紹在準(zhǔn)靜態(tài)和動(dòng)態(tài)沖擊下幾種常用的相變驅(qū)動(dòng)力與相變演化規(guī)律的模型。

4.2.1 相變驅(qū)動(dòng)力

相變驅(qū)動(dòng)力是與相變量功共軛的參量，取決于材料相變前后自由能的變化。在現(xiàn)有的研究中通常有兩種方法描述相變驅(qū)動(dòng)力與系統(tǒng)自由能的關(guān)系。第一種方法較為常見，直接將相變過程中Gibbs 自由能的變化作為相變的驅(qū)動(dòng)力。第二種方法以Turteltaub 和Suiker[101-103]為代表，認(rèn)為相變驅(qū)動(dòng)力在Gibbs 自由能的變化之外還應(yīng)額外加上一項(xiàng)相變熵增所導(dǎo)致的熱驅(qū)動(dòng)力，并指出這部分相變熵增并沒有引起系統(tǒng)Gibbs 自由能的變化。上述兩種方法在相變驅(qū)動(dòng)力與Gibbs 自由能的關(guān)系表達(dá)式上存在區(qū)別，這一區(qū)別實(shí)際上來源于自由能定義的具體表達(dá)式和物理含義的不同，但是最終得到的相變驅(qū)動(dòng)力是一致的，反映的物理規(guī)律是完全相同的。

相變驅(qū)動(dòng)力與各熱力學(xué)勢(shì)函數(shù)之間的關(guān)系不完全相同[98]

式中：ft表示與相變量功共軛的相變驅(qū)動(dòng)力，τm=FeTFeSFtr-T:表示相變過程的機(jī)械驅(qū)動(dòng)力，ψ、G分別為系統(tǒng)的Helmholtz 自由能與Gibbs 自由能。

Manchiraju 等[104]為研究Ni-Ti 形狀記憶合金中的相變過程，給出了系統(tǒng)的亥姆霍茲自由能表達(dá)式

式中：θ、θ0、θT分別表示瞬時(shí)溫度、初始溫度和相變溫度，λT為單位體積的比熱，htu表示不同相變變體之間的相容關(guān)系。式(88)中：第1 項(xiàng)表示彈性應(yīng)變能；第2 項(xiàng)表示與熱應(yīng)變對(duì)應(yīng)的自由能；第3 項(xiàng)表示相變過程中的化學(xué)能；第4 項(xiàng)表示馬氏體變體之間的相互作用能，反映了變體之間的相互協(xié)調(diào)能力。

相變驅(qū)動(dòng)力ft可根據(jù)式(87)中與亥姆霍茲自由能的關(guān)系得到

式中：S為PKⅡ應(yīng)力。

此外，F(xiàn)eng 等[105]、Greeff 等[106]認(rèn)為在沖擊高溫高壓作用下，材料的相變驅(qū)動(dòng)力主要依賴于系統(tǒng)的壓強(qiáng)和溫度。Cawkwell 等[38]根據(jù)材料在高壓下的狀態(tài)方程建立了一個(gè)評(píng)估兩相Gibbs 自由能的模型，并給出了含能材料RDX 的α與γ相在不同溫度與壓強(qiáng)下的Gibbs 自由能之差，見圖9[38]。該模型簡(jiǎn)化了相變驅(qū)動(dòng)力的組成部分，忽略了兩相的錯(cuò)配能以及相變的界面阻力等各個(gè)因素，通常只適用于材料的高壓沖擊相變過程。

4.2.2 相變演化規(guī)律

圖9 RDX 的α 相與γ 相Gibbs 自由能之差與溫度、壓強(qiáng)的關(guān)系[38]Fig. 9 Difference between Gibbs free energies of the α and γ RDX polymorphs as a function of pressure and temperature[38]

相變過程中新相體積分?jǐn)?shù)的演化依賴于相變的驅(qū)動(dòng)力ft和相變的臨界阻力fc。通常將相變過程中各種難以量化的阻力集合在一個(gè)總的相變臨界阻力fc。因此，馬氏體相變啟動(dòng)的臨界條件為系統(tǒng)可提供的自由能達(dá)到馬氏體相變的能壘

Turteltaub 等[101-102]、Tjahjanto 等[107]認(rèn)為相變啟動(dòng)之后新相的生長(zhǎng)速率依賴于驅(qū)動(dòng)力與相變臨界阻力的差值(ft-fc)，并依此提出了一個(gè)唯象的相變演化的動(dòng)力學(xué)關(guān)系

式中：k為相變演化系數(shù)。

Thamburaja 等[97,108-110]、Manchiraju 等[104]則傾向于認(rèn)為一旦相變驅(qū)動(dòng)力ft超過臨界阻力fc，馬氏體相變過程會(huì)立即發(fā)生并釋放系統(tǒng)的自由能，從而使相變驅(qū)動(dòng)力永遠(yuǎn)滿足自洽關(guān)系ft=fc，因此新相的生成 速率可根據(jù)相變的自洽關(guān)系求解

Barton 等[111]將上述唯象的相變演化動(dòng)力學(xué)關(guān)系式(式(91))與固體狀態(tài)方程相結(jié)合，模擬了Fe 在高壓下的沖擊相變 (α ?ε)，如圖10[111]所示。

圖10 Fe 沖擊相變的單晶模擬與多晶實(shí)驗(yàn)結(jié)果[111]Fig. 10 Single crystal Fe simulation data and polycrystal experimental data of shock-induced phase transformation[111]

此外，考慮動(dòng)態(tài)沖擊過程中相變驅(qū)動(dòng)力(Gibbs 自由能)依賴于壓力和溫度，Greeff 等[106]根據(jù)實(shí)驗(yàn)觀察到的在一定壓力區(qū)間內(nèi)新相生長(zhǎng)速率和壓力滿足指數(shù)關(guān)系，提出了適用于沖擊高壓下的相變演化模型

式中： ΔG(p)為兩相的Gibbs 自由能之差，A、B是與相變演化相關(guān)的材料參數(shù)。

5 動(dòng)態(tài)孿生變形

對(duì)于Ti[112-113]、Zr[10]等密排六方結(jié)構(gòu)(hcp)材料，在沖擊加載過程中，孿生變形與位錯(cuò)滑移是系統(tǒng)發(fā)生塑性變形和釋放能量的主要途徑，而相變往往需要更大的加載壓力和應(yīng)變速率[105]。孿生變形通常出現(xiàn)在滑移受阻引起的應(yīng)力集中區(qū)域，通過均勻剪切的方式形成鏡面對(duì)稱的兩部分晶體。由于孿生的應(yīng)變速率敏感性通常低于位錯(cuò)滑移的應(yīng)變速率敏感性[114]，即孿生阻力隨應(yīng)變率變化不大，材料在高應(yīng)變率下比準(zhǔn)靜態(tài)下更傾向于發(fā)生孿生變形而非塑性滑移。因此，Cu[12]、Al[11,115]等滑移系較多的fcc 材料在高壓沖擊過程中也能觀察到顯著的孿生變形。

5.1 孿生運(yùn)動(dòng)學(xué)

Kalidindi[116]最先將孿生變形引入晶體塑性框架中，模擬了低層錯(cuò)能的多晶fcc 金屬α-黃銅以及hcp 金屬Zr 中的織構(gòu)演化。Kalidindi 模型沿用了變形梯度的彈塑性分解

而塑性速度梯度包含材料的塑性滑移與孿生變形兩部分

對(duì)于體積分?jǐn)?shù)通常較小的fcc 和bcc 孿晶而言，往往只考慮母相中位錯(cuò)滑移而忽略孿晶內(nèi)部的塑性滑移。但對(duì)于孿晶體積分?jǐn)?shù)較大的情況如某些hcp 金屬，孿晶中的位錯(cuò)滑移則不能忽略。因此，塑性變形梯度中應(yīng)同時(shí)包含母相和孿晶中的位錯(cuò)滑移

式中：fβ表示孿晶的體積分?jǐn)?shù)，mα和nα分別表示柏氏矢量與滑移面法向，Qβ為孿生變形前后晶體位相的旋 轉(zhuǎn)矩陣

絕熱剪切帶是材料在高速動(dòng)態(tài)變形過程中由于局部溫升過高引起熱軟化，在微觀上發(fā)生剪切變形局域化的現(xiàn)象。由于變形時(shí)間極短，變形產(chǎn)熱來不及擴(kuò)散，剪切變形局域化幾乎在絕熱條件下進(jìn)行。絕熱剪切帶表現(xiàn)為細(xì)長(zhǎng)的窄帶，其演變包括萌生、擴(kuò)展和斷裂3 個(gè)階段。

材料在動(dòng)態(tài)破壞過程中的微觀結(jié)構(gòu)動(dòng)態(tài)演化缺少直接的測(cè)量結(jié)果，動(dòng)態(tài)晶體塑性有限元方法對(duì)動(dòng)態(tài)破壞力學(xué)響應(yīng)和微觀結(jié)構(gòu)演化的預(yù)測(cè)及分析對(duì)于相關(guān)實(shí)驗(yàn)研究尤為重要。關(guān)于動(dòng)態(tài)層裂破壞的模擬，主要介紹兩種將損傷引入晶體塑性模型的方法；關(guān)于絕熱剪切帶的模擬，分別從忽略損傷和直接引入損傷兩種思路對(duì)動(dòng)態(tài)晶體塑性模擬方法進(jìn)行具體介紹。

6.1 層 裂

層裂是材料動(dòng)態(tài)破壞的一種特殊形式。受高速?zèng)_擊加載時(shí)，材料內(nèi)部?jī)蓚€(gè)卸載稀疏波相互作用，一旦材料中局部應(yīng)力達(dá)到材料的動(dòng)態(tài)強(qiáng)度極限即發(fā)生斷裂，如圖11[52]所示。因此，層裂是結(jié)構(gòu)動(dòng)態(tài)響應(yīng)（波的交互作用）和材料動(dòng)態(tài)響應(yīng)（材料動(dòng)態(tài)強(qiáng)度極限或破壞準(zhǔn)則）共同作用的結(jié)果[51]。

圖11 沖擊變形過程中波的傳播及層裂現(xiàn)象(a)、3 個(gè)時(shí)刻的應(yīng)力波形(b)和3 個(gè)位置的應(yīng)力歷史(c)[52]Fig. 11 Wave propagation and spalling phenomenon (a), stress profiles at three different times (b),as well as stress histories at three different positions (c) during shock deformation[52]

采用動(dòng)態(tài)晶體塑性有限元對(duì)沖擊層裂現(xiàn)象的研究主要有兩種方法：（1）在動(dòng)態(tài)晶體塑性模型中引入層裂破壞準(zhǔn)則，通過顯式地刪除單元實(shí)現(xiàn)材料的損傷形核及其演化；（2）將材料視為有微孔洞的多孔晶體材料，在動(dòng)態(tài)晶體塑性模型中基于孔洞體積在變形過程中的演化引入材料的損傷演化。下面具體介紹這兩種動(dòng)態(tài)沖擊變形中層裂破壞的模擬方法。

6.1.1 考慮層裂破壞準(zhǔn)則的晶體塑性模型

直接引入材料損傷是晶體塑性模擬動(dòng)態(tài)層裂破壞的一個(gè)重要方法。依據(jù)損傷形核的難易程度，將晶體材料劃分為不易損傷的晶內(nèi)材料和易出現(xiàn)損傷形核點(diǎn)的結(jié)構(gòu)(如晶界、第二相粒子處等)，對(duì)兩者分別建立不同的材料模型進(jìn)行研究。對(duì)晶內(nèi)材料建立動(dòng)態(tài)晶體塑性有限元模型，模擬動(dòng)態(tài)沖擊壓縮過程中的晶內(nèi)微觀塑性變形行為；對(duì)晶體材料中易出現(xiàn)損傷形核點(diǎn)的結(jié)構(gòu)，建立各向同性的宏觀模型進(jìn)行動(dòng)態(tài)變形模擬，并結(jié)合層裂破壞準(zhǔn)則對(duì)損傷形核進(jìn)行定量預(yù)測(cè)。一旦模型中某個(gè)單元滿足層裂破壞準(zhǔn)則即成為損傷形核點(diǎn)，人為地將相應(yīng)的單元抹去，即在模型中引入微孔洞。隨著動(dòng)態(tài)沖擊變形的繼續(xù)發(fā)生，微孔洞將會(huì)長(zhǎng)大、合并。

材料失效模型通常有兩類：內(nèi)聚應(yīng)力準(zhǔn)則和損傷積累準(zhǔn)則[129]。內(nèi)聚應(yīng)力準(zhǔn)則假定材料失效是瞬態(tài)行為，當(dāng)材料中的局部應(yīng)力達(dá)到晶體的內(nèi)聚應(yīng)力即理論斷裂強(qiáng)度時(shí)，材料立刻發(fā)生斷裂。損傷積累準(zhǔn)則假定材料失效依賴于變形歷史，斷裂的發(fā)生不僅與應(yīng)力大小相關(guān)，還與應(yīng)力持續(xù)時(shí)間相關(guān)?，F(xiàn)有的動(dòng)態(tài)層裂破壞晶體塑性有限元研究均簡(jiǎn)單地采用內(nèi)聚應(yīng)力準(zhǔn)則，考慮到高應(yīng)變速率變形中存在明顯的斷裂滯后現(xiàn)象，后續(xù)研究有待將損傷積累準(zhǔn)則引入動(dòng)態(tài)晶體塑性模型中。下面介紹兩種典型的內(nèi)聚應(yīng)力準(zhǔn)則：最大正拉應(yīng)力準(zhǔn)則和復(fù)合模式的最大應(yīng)力準(zhǔn)則。

最大正拉應(yīng)力準(zhǔn)則可表示為[130]

當(dāng)材料所受正拉應(yīng)力(σ)達(dá)到臨界值(σc)時(shí)，則發(fā)生層裂破壞。Zhang 等[131]、Lloyd 等[132]分別以晶界和第二相粒子作為損傷形核點(diǎn)，采用該準(zhǔn)則判斷裂紋形核，通過動(dòng)態(tài)晶體塑性有限元模擬沖擊變形卸載過程中的層裂破壞現(xiàn)象。

復(fù)合模式的最大應(yīng)力準(zhǔn)則可表示為[23-24,57-58]

式中：s、τ分別為界面上的正應(yīng)力和切應(yīng)力，s0、τ0分別為參考正應(yīng)力和參考切應(yīng)力， ΔT為變形過程中的溫升，s1、τ1分別為依賴于溫升的臨界正應(yīng)力增量和臨界切應(yīng)力增量。當(dāng)材料局部正應(yīng)力或切應(yīng)力超過對(duì)溫度線性依賴的內(nèi)稟界面強(qiáng)度時(shí)，均會(huì)引起材料的局部裂紋形核，隨后由正應(yīng)力和切應(yīng)力主導(dǎo)的損傷演變略有不同。Clayton[23,57-58]、Vogler 等[24]將晶界視為損傷易形核點(diǎn)，采用該準(zhǔn)則判斷裂紋形核，模擬結(jié)果如圖12(a)和圖12(b)所示。樣品右側(cè)區(qū)域?qū)恿验_始形核，該區(qū)域材料承受著較大的拉伸應(yīng)力，即層裂區(qū)域附近具有較大的彈性能。動(dòng)態(tài)層裂破壞表現(xiàn)為晶間斷裂，并且飛片的沖擊加載速度越大，材料的層裂破壞越強(qiáng)烈，如圖12(c)和圖12(d)所示。

圖12 鉛合金動(dòng)態(tài)晶體塑性有限元模擬結(jié)果：(a)層裂形核時(shí)的壓力，(b)層裂形核時(shí)的彈性能密度，(c)經(jīng)250 m/s 沖擊加載層裂面附近的等效應(yīng)力；(d)經(jīng)350 m/s 沖擊加載層裂面附近的等效應(yīng)力[24]Fig. 12 Dynamic crystal plasticity finite element simulation results of lead alloy: (a) pressure of spalling nucleation;(b) elastic energy density of spalling nucleation; (c) equivalent stress near the spalling surface under 250 m/s shock loading; (d) equivalent stress near the spalling surface under 350 m/s shock loading[24]

6.1.2 基于孔洞體積演化的晶體塑性模型

圖14 晶粒取向和應(yīng)力三軸度對(duì)孔洞合并的臨界狀態(tài)變量的影響[133]Fig. 14 Influence of grain orientation and stress triaxiality on critical state variables for void coalescence[133]

6.2 絕熱剪切帶

絕熱剪切帶是材料在動(dòng)態(tài)高應(yīng)變速率加載變形過程中常發(fā)生的剪切變形局域化現(xiàn)象，最終會(huì)導(dǎo)致材料動(dòng)態(tài)斷裂[136-137]。Zener 和Hollomon[138]最早提出絕熱剪切帶的形成是熱軟化克服材料加工硬化的結(jié)果，并將微觀絕熱剪切帶的形成與材料宏觀熱塑性軟化(即變形后期應(yīng)力-應(yīng)變曲線發(fā)生下降)聯(lián)系在一起，見圖15[51]中Stage 3。絕熱剪切帶的基本特征包括：微觀上觀察到的絕熱剪切帶、宏觀應(yīng)力-應(yīng)變曲線上發(fā)生熱塑性失穩(wěn)、高速變形過程接近絕熱過程[51]。剪切局域化機(jī)理包含兩方面競(jìng)爭(zhēng)因素：促進(jìn)形成剪切局域化的正面因素(即熱軟化和幾何軟化)和不利于剪切局域化的負(fù)面因素(即應(yīng)變硬化和應(yīng)變率硬化)[52]。

高速?zèng)_擊變形過程中絕熱剪切帶的晶體塑性模擬研究主要有兩種方法：（1）忽略損傷，直接模擬沖擊變形過程中的塑性變形和溫度演化，根據(jù)模擬結(jié)果在后處理中觀測(cè)和分析絕熱剪切帶現(xiàn)象；（2）引入材料微觀損傷，模擬動(dòng)態(tài)沖擊變形過程中絕熱剪切帶的形核及演化。下面分別對(duì)這兩種方法進(jìn)行具體介紹。

6.2.1 忽略損傷的晶體塑性模型

圖15 經(jīng)典應(yīng)力-應(yīng)變曲線上塑性變形的3 個(gè)階段(Stage 1：均勻變形；Stage 2：非均勻變形；Stage 3：宏觀熱塑性失穩(wěn))[51]Fig. 15 Three stages of plastic deformation appeared on classical stress-strain curve (Stage1: homogeneous deformation; Stage2: inhomogeneous deformation;Stage3: macroscopic thermoplastic instability)[51]

高速變形下絕熱剪切的動(dòng)態(tài)晶體塑性有限元研究可以追溯到20 世紀(jì)90 年代，采用忽略損傷的動(dòng)態(tài)晶體塑性有限元能夠直接對(duì)高溫、高壓、高應(yīng)變速率下絕熱剪切帶的形成及演變進(jìn)行模擬。Hines 等[139]、Lee 等[140]基于單晶體彈塑性理論，忽略熱效應(yīng)，采用經(jīng)典的應(yīng)變率相關(guān)的唯象硬化本構(gòu)方程，粗略地模擬了動(dòng)態(tài)加載下的剪切帶現(xiàn)象。Bronkhorst 等[141]、Tajalli 等[59]考慮塑性變形產(chǎn)熱以及熱軟化效應(yīng)，基于高速變形中的絕熱假設(shè)進(jìn)行熱力耦合計(jì)算，采用唯象的動(dòng)態(tài)晶體塑性有限元模擬分析絕熱剪切局域化現(xiàn)象?？紤]到動(dòng)態(tài)變形溫升顯著并且對(duì)材料彈塑性本構(gòu)具有重要影響[142]，Bargmann等[143]進(jìn)一步考慮熱彈性耦合效應(yīng)和熱傳導(dǎo)效應(yīng)，對(duì)動(dòng)態(tài)變形溫升進(jìn)行更加精確的計(jì)算。他們同時(shí)考慮塑性產(chǎn)熱、彈性產(chǎn)熱以及熱傳導(dǎo)對(duì)溫升的貢獻(xiàn)，采用動(dòng)態(tài)梯度晶體塑性有限元模擬104～107s?1應(yīng)變速率范圍的材料沖擊變形行為。通過比較不同應(yīng)變速率、相同變形程度的模擬結(jié)果，發(fā)現(xiàn)相比于高應(yīng)變速率(107s?1)的動(dòng)態(tài)變形，材料在相對(duì)較低應(yīng)變速率(5×105s?1)下的動(dòng)態(tài)變形過程中熱傳導(dǎo)效果明顯，其溫度梯度更為平緩，如圖16 所示。

圖16 不同累積滑移速率變形=0.05 時(shí)的溫升云圖[143]Fig. 16 Distribution of temperature increase when =0.05 for different accumulated slip rates[143]

采用上述動(dòng)態(tài)晶體塑性有限元模型已經(jīng)能夠模擬絕熱剪切帶的形成，根據(jù)模擬得到的應(yīng)力-應(yīng)變曲線上熱塑性失穩(wěn)轉(zhuǎn)折點(diǎn)的力學(xué)狀態(tài)(如應(yīng)變、應(yīng)變率、剪切模量或變形能等)，進(jìn)一步給出熱塑性失穩(wěn)等效判據(jù)，能夠有效地對(duì)微觀上可觀察到的絕熱剪切帶形成的臨界點(diǎn)進(jìn)行預(yù)測(cè)和分析。當(dāng)材料受到動(dòng)態(tài)高速載荷時(shí)，雖然應(yīng)力在屈服點(diǎn)附近可能存在短暫下降，但是宏觀應(yīng)力-應(yīng)變曲線整體上首先經(jīng)歷一段上升，隨后曲線發(fā)生持續(xù)性下降，如圖15 所示。通常認(rèn)為曲線上峰值轉(zhuǎn)折點(diǎn)(即dσ/dε= 0 或dτ/dγ= 0)處宏觀上開始發(fā)生熱塑性失穩(wěn)，同時(shí)意味著微觀上開始形成絕熱剪切帶。發(fā)生宏觀熱塑性失穩(wěn)時(shí)，除宏觀應(yīng)力之外，其他的變形或狀態(tài)參量也會(huì)達(dá)到相應(yīng)的臨界值，相應(yīng)地就有其他等效判據(jù)，共有4 種：臨界應(yīng)變判據(jù)[144-147]、臨界應(yīng)變率判據(jù)[148]、臨界剪切模量判據(jù)[149-151]、臨界變形能判據(jù)[152]。Zhang 等[147]通過動(dòng)態(tài)晶體塑性有限元模擬得到的等效塑性應(yīng)變?cè)茍D上出現(xiàn)了絕熱剪切帶，如圖17 所示。進(jìn)一步采用臨界應(yīng)變判據(jù)能夠預(yù)測(cè)宏觀熱塑性失穩(wěn)臨界點(diǎn)，即微觀絕熱剪切帶形成的臨界狀態(tài)。他們根據(jù)模擬得到的宏觀應(yīng)力-應(yīng)變曲線，對(duì)熱塑性失穩(wěn)轉(zhuǎn)折點(diǎn)處的應(yīng)變狀態(tài)進(jìn)行統(tǒng)計(jì)(見圖18)，發(fā)現(xiàn)動(dòng)態(tài)壓縮過程中形成絕熱剪切帶的臨界應(yīng)變約4.5%，動(dòng)態(tài)剪切過程中的臨界應(yīng)變約為3.7%。

圖17 hcp 單晶和多晶樣品在105 s-1 應(yīng)變率下的絕熱剪切局域化[147]Fig. 17 Adiabatic shear localization of hcp single crystal and polycrystalline samples under 105 s-1 strain rate[147]

6.2.2 引入損傷的晶體塑性模型

從微觀角度考慮損傷的形核和演變機(jī)制是研究動(dòng)態(tài)絕熱剪切破壞的另一種重要方法。將損傷作為材料狀態(tài)變量引入動(dòng)態(tài)晶體塑性模型時(shí)，需要考慮3 個(gè)方面：引入損傷的彈性本構(gòu)、損傷驅(qū)動(dòng)力以及變形過程中損傷的演變。Boubaker 等[153]在超彈性本構(gòu)方程中引入材料損傷，將彈性自由能視為彈性應(yīng)變、溫度和損傷變量的函數(shù)，將塑性自由能視為滑移速率和損傷變量的函數(shù)，將應(yīng)力對(duì)損傷的依賴歸因于彈性模量隨損傷變量的變化，從而得到依賴于損傷的自由能和應(yīng)力

圖18 動(dòng)態(tài)沖擊載荷下6 種織構(gòu)材料中形成絕熱剪切帶的臨界應(yīng)變(Vpeak = 20 m/s)[147]Fig. 18 Critical strain of adiabatic shear band nucleated in 6 different texture materials under dynamic shock loading (V peak=20 m/s) [147]

式中：ψ為比自由能，ψe、ψp分別為比自由能的彈性和塑性分量，S、Ee分別為PKⅡ應(yīng)力和彈性格林-拉格朗日應(yīng)變張量，T為瞬態(tài)溫度，D為損傷變量，為滑移系α上的塑性滑移速率，ρ0為初始材料密度，C?(T,D)為 依賴于溫度和損傷的彈性模量， αΔT表示熱膨脹應(yīng)變。這里與式(116)相比增加了熱膨脹構(gòu)型，因此在應(yīng)力表達(dá)式中采用了純彈性應(yīng)變，詳見2.2.2 節(jié)。其次，熱力學(xué)損傷驅(qū)動(dòng)力可以通過自由能對(duì)損傷的依賴性得到

式中：Y為損傷驅(qū)動(dòng)力，gα為滑移系α上的滑移阻力。顯然，損傷驅(qū)動(dòng)力和損傷變量是功共軛的。最后，損傷變量在動(dòng)態(tài)變形過程中的演變可以通過損傷驅(qū)動(dòng)力與損傷演化阻力比值的冪函數(shù)表示

式中：Yr為參考損傷演化阻力，Lp為速度梯度張量的塑性分量。在初始未變形條件下，損傷驅(qū)動(dòng)力和損傷變量均為零；當(dāng)材料受沖擊載荷開始變形，則存在損傷驅(qū)動(dòng)力，損傷變量不再為零，意味著材料發(fā)生損傷形核；在隨后的變形過程中，損傷變量按式(121)隨時(shí)間而變化，表示材料中損傷的演變規(guī)律。通過這樣的方式，將連續(xù)介質(zhì)損傷模型與晶體塑性模型進(jìn)行耦合，能夠有效地模擬動(dòng)態(tài)絕熱剪切破壞現(xiàn)象。

7 總結(jié)與展望

動(dòng)態(tài)晶體塑性有限元對(duì)于研究高溫、高壓、高應(yīng)變速率加載條件下的材料動(dòng)態(tài)沖擊力學(xué)響應(yīng)和微觀結(jié)構(gòu)演化具有重要作用，但是如何對(duì)材料在動(dòng)態(tài)沖擊下微觀變形的物理機(jī)制進(jìn)行準(zhǔn)確的描述仍然存在許多的發(fā)展空間。

在理論模型方面，動(dòng)態(tài)晶體塑性有限元能夠有效地反映晶體塑性滑移的各向異性、對(duì)微觀組織的依賴性以及考慮動(dòng)態(tài)沖擊下各種影響變形的宏觀因素與微觀機(jī)制。與準(zhǔn)靜態(tài)晶體塑性模型相比，動(dòng)態(tài)晶體塑性模型考慮了材料在高壓下的非線性彈性響應(yīng)、顯著溫升引起的熱彈性耦合效應(yīng)以及高速運(yùn)動(dòng)位錯(cuò)具有的聲子拖曳效應(yīng)等。但是目前對(duì)于沖擊下位錯(cuò)高速滑移過程中位錯(cuò)之間的相互作用與位錯(cuò)形核、增殖的演化規(guī)律缺乏足夠的認(rèn)識(shí)，在未來的研究中可基于多尺度模擬框架，從分子動(dòng)力學(xué)與位錯(cuò)動(dòng)力學(xué)出發(fā)，探究位錯(cuò)在高速滑移的滑移特征和演化規(guī)律。

在應(yīng)用研究方面，動(dòng)態(tài)晶體塑性有限元能夠應(yīng)用于位錯(cuò)滑移、相變、孿晶和動(dòng)態(tài)損傷破壞等變形與破壞機(jī)制研究中?，F(xiàn)有的動(dòng)態(tài)晶體塑性模型通過引入局部相變體積分?jǐn)?shù)或?qū)\晶體積分?jǐn)?shù)及其演化規(guī)律，研究相變與孿晶對(duì)材料變形的影響，但這種方法只能均勻化地將相變與孿晶的影響簡(jiǎn)化為對(duì)材料點(diǎn)應(yīng)力的影響，無法反映真實(shí)的相變過程與孿晶演化過程。在未來研究中可考慮將動(dòng)態(tài)晶體塑性模型與相場(chǎng)模型相結(jié)合描述動(dòng)態(tài)沖擊下的相變和孿生過程。除了將動(dòng)態(tài)晶體塑性模型應(yīng)用于動(dòng)態(tài)損傷破壞的定性研究，還可分別結(jié)合宏觀破壞準(zhǔn)則和微觀損傷形核及演化準(zhǔn)則，對(duì)動(dòng)態(tài)沖擊變形過程中的層裂和絕熱剪切帶進(jìn)行定量預(yù)測(cè)和分析。從微觀損傷演化的角度，對(duì)動(dòng)態(tài)層裂和絕熱剪切破壞的動(dòng)態(tài)晶體塑性研究尚處于發(fā)展階段，目前的損傷主要以引入損傷系數(shù)進(jìn)行唯象的模擬為主，后續(xù)研究中需要更多地考慮基于物理的晶體動(dòng)態(tài)損傷機(jī)制。

表 A1 運(yùn)動(dòng)學(xué)符號(hào)說明
Table A1 Symbol description of kinematics

Symbols Description F(Fe, Fp, Fθ ) Deformation gradient including elastic, plastic and thermal components L(Le, Lp, Lθ ) Velocity gradient including elastic, plastic and thermal components Re Rotation tensor Ue Right stretch tensor α Thermal expansion coefficient tensor

表 A2 熱力學(xué)符號(hào)說明
Table A2 Symbol description of thermodynamics

Symbols Description Symbols Description Dint Intrinsic dissipation of the system K0 Bulk modulus at zero pressure ψ Helmholtz free energy K′ Pressure derivative of bulk modulus s Entropy of the system TD Debye temperature T Temperature R Molar gas constant KT Isothermal bulk modulus Mmol Molar mass of the material cV Heat capacity at constant volume kB Boltzmann constant Γ Grüneisen coefficient XTN Variables related to the lattice thermal vibration qn Internal variables for microscopic defects such as dislocations in materials XTE Variables related to the electron activation

表 A3 塑性本構(gòu)符號(hào)說明
Table A3 Symbol description of plastic constitution

Symbols Description Symbols Description λ α Mean spacing between obstacles ρα for Forest dislocation density The drag-dominated mean transit time between obstacles Qα Activation energy B Viscous drag coefficient gα Slip resistance ˙ραnuc The nucleation rate gαath Athermal slip resistance ˙ραhom The homogeneous nucleation rate hαβ Hardening coefficient ˙ραhet The heterogeneous nucleation rate ρα Total dislocation density ˙ρα mult The multiplication rate ραm Mobile dislocation density ˙ραtrap The trapping rate τα Resolved shear stress tαr ρα i Immobile dislocation density ˙ραann The annihilation rate bα Burgers vector da Capture distance of annihilation vα Velocity of mobile dislocations tαw The thermal activation-dominated waiting time at a barrier ˙γα Slip rate on slip system α

表 A4 超彈性本構(gòu)符號(hào)說明
Table A4 Symbol description of hyper-elastic constitution

Symbols Description Symbols Description I Second-order unit tensor ^Ee Isochoric strain in expanded configuration Ee Elastic Green-Lagrange strain^^Ee Isochoric strain in configuration I Ce Elastic right Cauchy-Green tensor Ee Volumetric strain in configuration I^Fe Isochoric part of elastic deformation S Second Piola-Kirchhoff stress Fe Volumetric expansion

表 A5 相變、孿晶與動(dòng)態(tài)破壞符號(hào)說明
Table A5 Symbol description of phase transformation, twining and damage

Symbols Description Symbols Description Ftr Deformation gradient of phase transformation Sβ tw Twin resistance of twin system vp Volume fraction of the parent phase ρdeb Dislocationdebris density vt Volume fraction of the new phase t dmfp Dislocation mean free path related to the volume fraction of twin vN Volume fraction of all new phases φ Void volume fraction ft Driving force of phase transformation Fd Volumetricpartofplastic deformation gradient in porous crystal plastic model f β Volume fraction of twin Yr Resistance of damage evolution γtw Characteristic shear strain of twining