高溫、高壓、高應變速率動態過程晶體塑性有限元理論模型及其應用

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

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

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

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

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

1 變形運動學

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

1.1 彈塑性變形分解

經典的運動學構型如圖1 所示。晶體材料由初始構型經變形梯度(F)變換至現時構型，如圖1中的黑色箭頭所示。Hill 等[21]、Asaro 等[22]建立中間構型Ⅰ，將彈性和塑性變形進行乘和分解

圖1 經典的晶體運動學構型Fig. 1 Classical configurations of crystal kinematics

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

式中：Re為旋轉張量，Ue為右伸長張量。材料由中間構型Ⅰ，經彈性拉伸(Ue)變換至中間構型Ⅱ，再經晶格旋轉(Re)變換至現時構型，如圖1 中的藍色箭頭所示。

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

1.2 熱膨脹構型

材料在動態沖擊加載下往往存在顯著的溫升現象，而溫升引起的熱膨脹會導致材料內部應力、應變的變化，因此如何考慮熱膨脹效應的影響是動態晶體塑性模型的一個重要問題。

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

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

相應地，現時構型中的速度梯度可分解為

式中：Lθ、分別為現時構型和熱構型中的熱變形速度梯度，是溫度變化率 (T˙)的函數

式中：α為熱膨脹系數張量。

2 超彈性本構框架

通常在晶體塑性模型中有超彈性和次彈性兩種描述材料彈性行為的本構關系。超彈性本構是指材料存在一個與應變相關的能量函數，其應力與應變全量之間的關系由能量函數決定，具有做功獨立于加載路徑的特點。次彈性本構描述應力率與應變率之間的關系，由于大變形下彈性模量是隨應變全量變化的，其實際應用比較復雜，通常只用于小彈性變形情況[28](這種情況下彈性模量可視為常數)。考慮到實驗中通常采用全量形式的高壓固體狀態方程描述材料的非線性容變關系，因此研究材料在沖擊載荷下的大彈性變形行為時，往往采用超彈性本構來描述。

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

2.1 超彈性本構的熱力學框架

3.2 基于位錯行為的動態晶體塑性框架

唯象的晶體塑性本構通常存在適用范圍的局限性，只能在一定應變率區間內描述材料的塑性變形規律。因此為了滿足不同應變率的加載條件，并同時描述各微觀結構的演化，通常采用Orowan 方程作為流動法則的模型框架

Orowan 方程模型框架主要包含可動位錯密度的演化以及平均位錯速度與各內變量之間的關系兩部分內容。可動位錯密度演化部分需要考慮位錯的形核、增殖、纏結、湮滅等微觀過程。位錯平均滑移速度部分可由位錯滑移障礙物之間的平均距離除以在障礙物之間運動的平均時間得到。而位錯在障礙物之間的平均運動時間可分成跨越障礙物的等待時間以及在障礙物之間的滑移時間兩部分，位錯的兩個運動過程可分別由Arrhenius 熱激活機制和黏性拖曳機制描述。

3.2.1 位錯滑移速度

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

第1 類模型是基于單根位錯實驗結果的簡單擬合，只描述了位錯運動速度與應力之間的關系，不僅忽略了位錯網絡對滑移的阻礙作用，也沒有考慮溫度對滑移的影響。第2 類模型只描述位錯在某一應力區間與溫度范圍內的單一規律，適用范圍較窄，缺乏對位錯完整滑移過程的描述。第3 類模型描述了位錯滑移的整個物理過程，適用于不同的應力區間和溫度范圍。本節將簡單介紹前兩類模型，重點介紹適用范圍最廣的第3 類模型。

（1）簡單指數模型

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

基于上述實驗結果，Gilman[63]提出了一種簡單的位錯速度與應力關系的指數模型

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

上述簡單的指數模型滿足實驗上位錯速度的邊界條件：當τ→0 時，v→0；當τ→∞時，v→v*。但該模型不僅缺乏對位錯運動機制的解釋，也沒有反映出位錯在低應力區間的轉變特征。

為了將模型應用于更廣的應力區間，Gilman[64]提出了在單指數模型的基礎上添加一項在低應力區間內的位錯速度

此外，Johnson 與Barker[65]相繼提出了關于位錯運動速度的指數模型，但這些模型均只是根據位錯的運動趨勢提出的唯象模型，缺乏相應的物理機制，且適用范圍較窄。

（2）單一機制模型

單一機制模型包括只適用于低應力區間內的熱激活機制以及只適用于高應力區間的線性拖曳機制。Alankar 等[66]、Zhang 等[67]、Monnet 等[68]采用熱激活機制來描述位錯在低應力區間內的滑移速度

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

式中：B為拖曳系數。

單一機制模型描述了位錯在某一應力區間內的滑移規律，適用于模擬材料在特定加載區間內的變形過程。但實際上材料在不同時間或空間尺度上的應變率往往存在很大差距，單一機制模型難以描述這種復雜的變形過程。

（3）統一物理模型

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

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

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

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

式中：CT為纏結系數。

位錯的湮滅指符號相反的位錯在一定捕獲距離內相遇所發生的消亡現象，決定了位錯飽和密度的存在。Luscher 等[39]、Alankar 等[66]根據上述湮滅過程建立了簡單的湮滅模型

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

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

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

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

4 動態相變

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

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

4.1 馬氏體相變運動學

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

式中：Fe、Fp、Ftr分別表示彈性變形、塑性變形以及相變對總變形梯度的貢獻。

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

式中：vp=1?vN表示母相的體積分數，

t=1為相變后新相的總體積分數，Nt為變體的數量，Nα和Nβ分別代表母相和新相的滑移系，和分別表示母相中的滑移方向和滑移面法向，和分別為新相中的滑移方向和滑移面法向。當相變產生的新相的位錯滑移能力較弱時，往往只計及母相中的塑性滑移[100]。

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

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

4.2 相變驅動力與演化準則

相變模型的關鍵在于如何描述相變過程的驅動力以及相變的演化規律。本節將重點介紹在準靜態和動態沖擊下幾種常用的相變驅動力與相變演化規律的模型。

4.2.1 相變驅動力

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

相變驅動力與各熱力學勢函數之間的關系不完全相同[98]

式中：ft表示與相變量功共軛的相變驅動力，τm=FeTFeSFtr-T:表示相變過程的機械驅動力，ψ、G分別為系統的Helmholtz 自由能與Gibbs 自由能。

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

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

相變驅動力ft可根據式(87)中與亥姆霍茲自由能的關系得到

式中：S為PKⅡ應力。

此外，Feng 等[105]、Greeff 等[106]認為在沖擊高溫高壓作用下，材料的相變驅動力主要依賴于系統的壓強和溫度。Cawkwell 等[38]根據材料在高壓下的狀態方程建立了一個評估兩相Gibbs 自由能的模型，并給出了含能材料RDX 的α與γ相在不同溫度與壓強下的Gibbs 自由能之差，見圖9[38]。該模型簡化了相變驅動力的組成部分，忽略了兩相的錯配能以及相變的界面阻力等各個因素，通常只適用于材料的高壓沖擊相變過程。

4.2.2 相變演化規律

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

相變過程中新相體積分數的演化依賴于相變的驅動力ft和相變的臨界阻力fc。通常將相變過程中各種難以量化的阻力集合在一個總的相變臨界阻力fc。因此，馬氏體相變啟動的臨界條件為系統可提供的自由能達到馬氏體相變的能壘

Turteltaub 等[101-102]、Tjahjanto 等[107]認為相變啟動之后新相的生長速率依賴于驅動力與相變臨界阻力的差值(ft-fc)，并依此提出了一個唯象的相變演化的動力學關系

式中：k為相變演化系數。

Thamburaja 等[97,108-110]、Manchiraju 等[104]則傾向于認為一旦相變驅動力ft超過臨界阻力fc，馬氏體相變過程會立即發生并釋放系統的自由能，從而使相變驅動力永遠滿足自洽關系ft=fc，因此新相的生成 速率可根據相變的自洽關系求解

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

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

此外，考慮動態沖擊過程中相變驅動力(Gibbs 自由能)依賴于壓力和溫度，Greeff 等[106]根據實驗觀察到的在一定壓力區間內新相生長速率和壓力滿足指數關系，提出了適用于沖擊高壓下的相變演化模型

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

5 動態孿生變形

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

5.1 孿生運動學

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

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

對于體積分數通常較小的fcc 和bcc 孿晶而言，往往只考慮母相中位錯滑移而忽略孿晶內部的塑性滑移。但對于孿晶體積分數較大的情況如某些hcp 金屬，孿晶中的位錯滑移則不能忽略。因此，塑性變形梯度中應同時包含母相和孿晶中的位錯滑移

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

絕熱剪切帶是材料在高速動態變形過程中由于局部溫升過高引起熱軟化，在微觀上發生剪切變形局域化的現象。由于變形時間極短，變形產熱來不及擴散，剪切變形局域化幾乎在絕熱條件下進行。絕熱剪切帶表現為細長的窄帶，其演變包括萌生、擴展和斷裂3 個階段。

材料在動態破壞過程中的微觀結構動態演化缺少直接的測量結果，動態晶體塑性有限元方法對動態破壞力學響應和微觀結構演化的預測及分析對于相關實驗研究尤為重要。關于動態層裂破壞的模擬，主要介紹兩種將損傷引入晶體塑性模型的方法；關于絕熱剪切帶的模擬，分別從忽略損傷和直接引入損傷兩種思路對動態晶體塑性模擬方法進行具體介紹。

6.1 層 裂

層裂是材料動態破壞的一種特殊形式。受高速沖擊加載時，材料內部兩個卸載稀疏波相互作用，一旦材料中局部應力達到材料的動態強度極限即發生斷裂，如圖11[52]所示。因此，層裂是結構動態響應（波的交互作用）和材料動態響應（材料動態強度極限或破壞準則）共同作用的結果[51]。

圖11 沖擊變形過程中波的傳播及層裂現象(a)、3 個時刻的應力波形(b)和3 個位置的應力歷史(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]

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

6.1.1 考慮層裂破壞準則的晶體塑性模型

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

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

最大正拉應力準則可表示為[130]

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

復合模式的最大應力準則可表示為[23-24,57-58]

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

圖12 鉛合金動態晶體塑性有限元模擬結果：(a)層裂形核時的壓力，(b)層裂形核時的彈性能密度，(c)經250 m/s 沖擊加載層裂面附近的等效應力；(d)經350 m/s 沖擊加載層裂面附近的等效應力[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 晶粒取向和應力三軸度對孔洞合并的臨界狀態變量的影響[133]Fig. 14 Influence of grain orientation and stress triaxiality on critical state variables for void coalescence[133]

6.2 絕熱剪切帶

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

高速沖擊變形過程中絕熱剪切帶的晶體塑性模擬研究主要有兩種方法：（1）忽略損傷，直接模擬沖擊變形過程中的塑性變形和溫度演化，根據模擬結果在后處理中觀測和分析絕熱剪切帶現象；（2）引入材料微觀損傷，模擬動態沖擊變形過程中絕熱剪切帶的形核及演化。下面分別對這兩種方法進行具體介紹。

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

圖15 經典應力-應變曲線上塑性變形的3 個階段(Stage 1：均勻變形；Stage 2：非均勻變形；Stage 3：宏觀熱塑性失穩)[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]

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

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

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

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

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

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

圖18 動態沖擊載荷下6 種織構材料中形成絕熱剪切帶的臨界應變(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Ⅱ應力和彈性格林-拉格朗日應變張量，T為瞬態溫度，D為損傷變量，為滑移系α上的塑性滑移速率，ρ0為初始材料密度，C?(T,D)為 依賴于溫度和損傷的彈性模量， αΔT表示熱膨脹應變。這里與式(116)相比增加了熱膨脹構型，因此在應力表達式中采用了純彈性應變，詳見2.2.2 節。其次，熱力學損傷驅動力可以通過自由能對損傷的依賴性得到

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

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

7 總結與展望

動態晶體塑性有限元對于研究高溫、高壓、高應變速率加載條件下的材料動態沖擊力學響應和微觀結構演化具有重要作用，但是如何對材料在動態沖擊下微觀變形的物理機制進行準確的描述仍然存在許多的發展空間。

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

在應用研究方面，動態晶體塑性有限元能夠應用于位錯滑移、相變、孿晶和動態損傷破壞等變形與破壞機制研究中。現有的動態晶體塑性模型通過引入局部相變體積分數或孿晶體積分數及其演化規律，研究相變與孿晶對材料變形的影響，但這種方法只能均勻化地將相變與孿晶的影響簡化為對材料點應力的影響，無法反映真實的相變過程與孿晶演化過程。在未來研究中可考慮將動態晶體塑性模型與相場模型相結合描述動態沖擊下的相變和孿生過程。除了將動態晶體塑性模型應用于動態損傷破壞的定性研究，還可分別結合宏觀破壞準則和微觀損傷形核及演化準則，對動態沖擊變形過程中的層裂和絕熱剪切帶進行定量預測和分析。從微觀損傷演化的角度，對動態層裂和絕熱剪切破壞的動態晶體塑性研究尚處于發展階段，目前的損傷主要以引入損傷系數進行唯象的模擬為主，后續研究中需要更多地考慮基于物理的晶體動態損傷機制。

表 A1 運動學符號說明
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 熱力學符號說明
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 塑性本構符號說明
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 超彈性本構符號說明
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 相變、孿晶與動態破壞符號說明
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