AZ31鎂合金板熱成形中的屈服和損傷: 本構實現與數值分析

周霞，文冬，沈夢祺, 宋尚雨

1.大連理工大學 工業裝備結構分析國家重點實驗室, 大連 116024 2.大連理工大學 工程力學系, 大連 116024 3.國際計算力學 研究中心, 大連 116024

作為一種輕質結構材料，鎂合金由于高的比強度、比剛度以及優異的減振性能在航空航天、汽車、電子等領域具有很好的應用潛力[1]。然而，由于鎂合金具有密排六方晶格結構，其室溫塑性差，難以塑性加工，因而其塑性加工產品的開發應用還遠不如鋼鐵、鋁、銅成熟[2]。此外，鎂合金力學性能還具有明顯的各向異性及拉壓不對稱性，塑性加工成形屈服本構模型中材料參數不是定值，而是隨著塑性流動而變化[3]，即屈服面形狀隨塑性流動而變化；隨著溫度升高，雖然鎂合金塑性成形能力增強[4]，但變形溫度、變形速率及應力狀態對其塑性變形能力也有重要影響。鎂合金板材的溫熱沖壓成形過程是一個復雜的多參數熱力耦合非線性過程，也是鎂合金板材加工成形的研究熱點，基于合適材料模型的鎂合金板材溫熱成形數值模擬可以準確預測板料成形性能和塑性行為及破壞方式。同時，由于鎂合金材料的成形性較差，需要從微觀角度研究其損傷的發展。因此對鎂合金板料沖壓成形各向異性屈服本構和損傷模型的研究具有重要意義[5]。

在過去的幾十年中，已有相關學者提出了許多金屬板料塑性成形各向異性屈服準則。其中Hill48各向異性模型[6]由于只包含6個各向異性參數，對于一些復雜結構金屬材料并不能準確描述。隨后Barlat等相繼提出了Yld91各向異性屈服函數[7]、 考慮了不同方向屈服應力和塑性應變比r值變化的Yld96屈服函數[8]以及對Yld96改進的Yld2000-2d屈服函數[9]。對于密排六方(Hexagonal Close Packed, HCP)晶格結構材料，由于存在拉壓不對稱性，上述模型并不一定準確適用于這種結構材料的成形。針對該問題，Hosford[10]通過對Hill各向異性屈服準則進行修正來描述材料各向異性屈服響應和拉壓不對稱性。但在該模型中由于各向異性和拉壓不對稱的耦合不能系統地確定各個材料參數。在此基礎上，Cazacu等和Plunkett等分別提出的CPB04本構模型[11]和改進的CPB06本構模型[12]，能夠描述鎂合金等HCP晶格材料的各向異性和拉壓不對稱性，且改進的CPB06本構模型由于通過多次線性變換引入更多的各向異性參數因而具有較高的精度。近年來，Tari和Worswick[13]考慮到AZ31鎂合金在塑性成形畸變硬化過程中，屈服面形狀隨著塑性變形而發生變化，提出了考慮參數演化的AZ31鎂合金各向異性屈服CPB06ex3ev本構模型，該模型雖能較準確地描述板材料的拉壓不對稱各向異性屈服及其屈服面形狀的演化，但在板成形模擬中還需要考慮準確的損傷破壞模型[14]。目前對金屬板材成形過程損傷行為的研究多基于連續介質損傷方法的定性預測，而基于細觀損傷力學及成形極限理論對板料成形性分析的研究日趨廣泛[15-16]，但針對鎂合金在溫熱成形過程中的各向異性屈服和細觀損傷演化耦合的研究還未見報道。

本文選用適用于HCP晶格結構，能同時反映各向異性和拉壓不對稱屈服的CPB06各向異性屈服準則來描述鎂合金的各向異性屈服，鎂合金的損傷準則采用適用于延性金屬的Gurson-Tvergaard-Needleman(GTN)損傷準則[17-18]。考慮到鎂合金的屈服面隨著塑性變形而變化，對CPB06各向異性屈服本構進行了修正，并結合GTN損傷模型，得到了考慮損傷演化和屈服面形狀變化的CPB06-GTN損傷本構模型，編譯得到了相應的VUMAT本構子程序，并嵌入到ABAQUS/Explicit中，采用單個單元對本構模型進行了驗證，然后將本構子程序應用到沖壓成形模擬中，預測了鎂合金的熱成形，模擬預測與實驗結果吻合較好。

1 本構模型

1.1 CPB06各向異性拉壓不對稱屈服本構

在平面應力條件下，CPB06本構模型的具體函數形式為[19]

(1)

式中：k為用于描述材料拉壓非對稱性的參數(-1≤k≤1)；a為與材料有關的參數，在使用時a的取值需要不斷調整，對于AZ31鎂合金可取2；Σi(i=1, 2, 3)為對應力偏張量線性變換并求主值的結果，可用式(2)和式(3)進行求解：

(2)

(3)

其中：式(2)中矩陣C為描述各向異性屈服的線性變換張量；S為Cauchy應力偏張量；式(3)中σ11、σ22、σ12為Cauchy應力偏張量對應的應力分量；式(1)中Σ1、Σ2、Σ3是Σ的主值，且

從以上結果推導可得出，當k=0，C為4階單位矩陣，且a=2時，該模型與Mises屈服準則等效。

1.2 GTN損傷本構

用來描述由于孔洞不斷長大、形核、聚合導致的宏觀材料軟化及破裂現象的GTN損傷模型[17-18]可表示為

(1+q3f*2)=0

(4)

材料進入塑性時由流動法則可計算得到塑性應變增量為

(5)

式中：dλ為一個正標量因子；Ψ=Ψ(σp,σe,Hα)為流動勢，對于相關流動法則φ和Ψ是一致的，Hα(α=1,2)為狀態變量σij為應力分量。在各向異性屈服GTN損傷模型式(5)中流動法則可表示為

(6)

式中：I為4階單位矩陣；n為流動方向；對于各向異性材料，流動方向表示為

(7)

式(6)中總的塑性應變增量可看做由兩部分組成：

(8)

(9)

式中：Δεp為宏觀塑性應變張量與靜水應力對應的塑性體積應變部分；Δεq為宏觀塑性應變張量與等效應力對應的等效塑性應變部分。

由式(9)消去Δλ可得到

(10)

由等效塑性功原理，可以得到基體材料等效塑性應變的演化方程：

(11)

(12)

材料的損傷演化包含原有孔洞的長大和新孔洞的形核：

(13)

(14)

在最終階段，考慮到孔洞的聚合和材料的失效，Needleman和Tvergaard[18]通過引入等效體積分數函數f*(f)對原有模型進行了修正：

(15)

1.3 修正的CPB06-GTN本構模型

本節基于多孔材料的CPB06-GTN宏觀各向異性屈服準則[20]，并對其進行了修正，使其適用于殼單元，損傷演化部分不變，具體形式為

(1+q3f*2)=0

(16)

式中：σq為CPB06屈服模型的等效應力；q1、q2、q3、σp、σm、f*以及損傷參數的演化均與GTN損傷模型相同。

σq=m·

(17)

式中：

m=

式(16)中：

(18)

(19)

圖1 應變率為和不同溫度(T=150, 200, 250 ℃)時基于CPB06模型擬合得到的屈服面形狀演化Fig.1 Evolution of yield surface with accumulated plastic strain obtained by anisotropic CPB06 model at strain rate of 0.001 s-1 and different temperatures (T=150, 200, 250 ℃)

表1 150 ℃時各向異性CPB06屈服本構參數隨等效塑性應變的變化Table 1 Variation of anisotropic CPB06 yield criterion parameters with effective plastic strain at 150 ℃

表2 200 ℃時各向異性CPB06屈服本構參數隨等效塑性應變的變化Table 2 Variation of anisotropic CPB06 yield criterion parameters with effective plastic strain at 200 ℃

表3 250 ℃時CPB06各向異性屈服本構參數隨等效塑性應變的變化Table 3 Variation of anisotropic CPB06 yield criterion parameters with effective plastic strain at 250 ℃

鎂合金材料的各向同性硬化行為通過定義軋制方向材料的等效塑性應變與拉伸應力的關系來描述，并以軋制方向作為參考方向，沿軋制方向單向拉伸的屈服應力和塑性應變關系可用式(20)簡化[13]，即Cowper-Symonds模型：

(20)

表4AZ31B鎂合金在不同溫度下的硬化參數

Table4HardeningparametersofAZ31Bmagnesiumalloyatdifferenttemperatures

TemperatureT/℃KnDP150183．0120．126572．569．23920075．6630．1080．0086．06825075．0930．0550．0150．371

(21)

2 數值算法實現及驗證

2.1 VUMA子程序實現

將改進后的塑性各向異性GTN損傷模型，借助用戶材料子程序VUMAT嵌入到有限元軟件ABAQUS的顯式分析模塊中，進行鎂合金的沖壓成形模擬，VUMAT子程序需要根據應變增量對應力進行更新，其中一個條件需要判斷材料是否進入塑性狀態，并求解出塑性應變；沖壓成形中的鎂合金板采用的是殼單元，不同于一般三維單元的是，殼單元VUMAT子程序除了需要完成應力更新外，還需要對厚度方向的應變進行更新；殼單元的各向異性屈服GTN損傷模型的VUMAT子程序算法分為彈性預測和塑性修正[23]，具體的應力更新計算過程如下：

步驟2假設應變增量為完全彈性，計算得到彈性預測應力：

(22)

(23)

步驟4塑性修正，通過迭代求解塑性應變確定流動方向n：

(24)

求解以下非線性方程組，保證等式成立：

(25)

U、V分別對應式(10)和式(4)的方程，這里采用Newton-Raphson法迭代求解兩個未知量k+1Δεp和k+1Δεq；設cp和cq分別為Δεp和Δεq的修正量，迭代式為

(26)

更新相應的等效應力σq和靜水應力σp：

(27)

(28)

(29)

式中：左上標k和k+1 表示迭代次數；當孔洞體積率f>fc時，等效孔洞體積分數f*=fc+δ(f-fc)。

結合硬化準則對基體材料的屈服應力進行更新：

(30)

繼續完成迭代，直到|U|和|V|<10-8時停止迭代，此時迭代收斂，式(25)成立，進入下一步。

步驟5完成應力更新，判斷孔洞體積率是否達到臨界值。

將由步驟4迭代求解得到的Δεp和Δεq代入式(8)可得到相應的塑性應變增量：

(31)

通過對預測應力進行修正即可得到當前狀態的真實應力：

σ=σtr-Ce:Δεp

(32)

式中：σtr為彈性預測應力；Ce為4階彈性模量；Δεp為塑性應變增量；當孔洞體積率達到臨界值時，則相應的控制單元刪除的狀態變量等于0，意味著單元失效。

步驟6對殼單元厚度方向的應變增量作更新。

(33)

式中：dσ11和dσ22分別為x方向和y方向的應力增量。

(34)

式中：dεp是宏觀塑性應變張量與靜水應力對應的塑性體積增量部分；dεq是宏觀塑性應變張量與等效應力對應的等效塑性應變增量部分；n11和n22分別為流動方向n在x和y方向的分量。

綜上，可得出厚度方向的應變增量為

dεq(-n11-n22)

(35)

步驟7結束當前時間步，進入下一時間步。

具體程序流程圖如圖2所示。

圖2 考慮本構參數演化的CPB06-GTN損傷模型VUMAT程序流程Fig.2 Flowchart of VUMAT subroutine for the coupled CPB06-GTN damage model considering yield loci evolution

2.2 單個單元測試

圖3 采用(VUMAT)子程序模擬的單軸拉伸和壓縮應力應變曲線及其與實驗數據[13]的比較Fig.3 Uniaxial tensile and compressive stress-strain curves for a single unit obtained by subroutine VUMAT and their comparison with experimental data[13]

3 沖壓成形算例分析

3.1 200 ℃等溫成形模擬

本節基于VUMAT子程序對AZ31B鎂合金板溫熱深拉成形進行了數值模擬，模擬使用的AZ31B鎂合金板厚度為1.57 mm，直徑為 228.6 mm；沖頭外徑為101.6 mm，沖頭圓角半徑為8.0 mm，凹模直徑為106 mm，凹模圓角半徑為8 mm，沖壓成形時鎂合金板由壓邊圈壓緊，通過在壓邊圈上施加壓邊力防止起皺，由于模型具有對稱性，為節省計算時間選取1/4模型進行分析。有限元模型如圖4所示，其中鎂合金采用S4R殼單元，厚度方向有5個積分點，其他模具均設置為離散剛體。由于實際成形時采取了潤滑措施，模擬時的摩擦系數取0.04。

圖5(a)顯示了拉深比為2.25、壓邊力為80 kN時深拉杯的壁厚分布云圖，STH表示深拉變形過程中深拉杯的厚度。由圖可知，此時無法完成深拉成形，鎂合金板深拉時其側壁處發生開裂，鎂合金深拉杯的破壞形態和位置與Tari等[24]實驗結果吻合。圖5(b)顯示了該條件下不同極限應變點(點1, 2, 3)的模擬結果與Li等[25]實驗測得的成形極限圖(Forming Limit Diagram, FLD)的比較，可以進一步看出，在成形極限曲線FLD下方的極限應變點2和點3處于安全區域，而極限應變點1則處于臨界破裂區，對應圖5(a)的不同部位，有限元模擬結果與成形極限圖結果吻合。這說明本文的各向異性GTN損傷模型能夠準確預測鎂合金溫熱深拉成形時的損傷失效。

圖4 深拉成形有限元網格模型Fig.4 Deep drawing finite element mesh model

當采用直徑為203.2 mm的鎂合金板，對應的拉深比為2.0，施加壓邊力為48 kN時，能夠順利完成深拉成形，成形后深拉杯的壁厚分布如圖6(a)所示。從圖6(a)可看出，成形后深拉杯的壁厚呈不均勻分布，杯壁區靠近沖頭圓角半徑附近壁厚發生較大的減薄，且軋制方向(X軸方向)側壁上的厚度較薄，而橫向(Z軸方向)側壁的厚度與軋制方向側壁厚度相比較厚，但由于高溫下各向異性減弱，因而并沒有出現明顯的制耳。同樣條件下模擬所得的沖頭力-位移曲線與Tari等[24]的實驗結果對比如圖6(b)所示。上述成形條件下的對比結果也說明了本文的本構模型能較好地模擬鎂合金的完整成形。

圖5 200 ℃等溫成形后帶裂紋深拉杯厚度分布模擬結果及其與成形極限圖的對比Fig.5 Simulation of thickness distribution of drawn cup with crack after deep drawing at 200 ℃ and its comparison with the forming limit diagram

圖6 深拉成形壁厚分布及沖頭力-位移曲線Fig.6 Thickness distribution in deep drawing and curves of punch force vs displacement

3.2 非等溫成形熱力耦合模擬

本節通過在考慮參數演化的各向異性屈服CPB06-GTN損傷模型編譯的VUMAT子程序中加入溫度相關變量，將子程序嵌入到ABAQUS軟件中進行了非等溫拉深成形熱力耦合有限元模擬。在ABAQUS中提供了顯示動態熱力耦合分析步，可將VUMAT子程序應用于該分析步中，VUMAT子程序可根據單元節點的溫度來調整材料參數，完成應力更新，并依據產生的塑性變形計算由塑性變形導致的溫度升高，繼而對溫度變量進行更新。這里，設置鎂合金板的單元為S4RT，該單元包含了溫度自由度，而其他模具均設置為離散剛體，并施加溫度邊界條件。模具和壓邊圈溫度設為245 ℃，沖頭溫度和板料初始溫度分別為171 ℃和175 ℃。

在有限元模擬時，通過在模具上施加溫度邊界條件來控制模具的溫度，鎂合金板的熱參數及與模具間的熱傳導率參數采用Lee等[26]給出的AZ31鎂合金在200 ℃的熱參數，具體數據如表5所示。變形過程中模具和鎂合金板之間會發生熱傳遞，在模擬時為了節省計算時間，提高了拉深速度；在溫熱成形時為了讓模擬的溫度與速度相匹配，也對鎂合金與時間相關的材料參數進行了放大處理。由于對計算速度放大了100倍，這里對熱傳導率也放大了100倍。參考實驗條件，對直徑為228.6 mm的鎂合金板，在對應拉深比為2.25，壓邊力為80 kN時的溫熱深拉非等溫成形，進行了數值模擬。

當模具和壓邊圈溫度為245 ℃，沖頭溫度為171 ℃時，能夠順利完成拉深成形，模擬的溫度(NT11)和壁厚分布如圖7所示。從圖7(a)可看出，沖頭與深拉杯的底部中心處溫度約為171 ℃，而板被壓邊圈壓緊部分的凸緣處溫度則和模具及壓邊圈溫度一致，約為245 ℃，側壁溫度在171～245 ℃之間。從圖7(b)可看出，由于沖頭溫度較低，沖頭圓角處強度較高，該處變薄并不嚴重，而側壁則由于承受了較大的徑向拉應力，導致了變薄。同時由于各向異性的影響，鎂合金板在軋制方向(X軸方向)的側壁變薄較嚴重。與前述等溫成形的結果對比可以看出，非等溫深拉成形可通過調整不同部位的溫度，提高拉深比，改善拉深成形性能。

表5 AZ31鎂合金的熱物理性能參數Table 5 Thermophysical properties parameters of AZ31 magnesium alloy

圖7 非等溫成形時的溫度分布與壁厚分布Fig.7 Temperature distribution and thickness distribution in non-isothermal deep drawing forming

此外，還將非等溫深拉成形模擬結果與已有實驗數據[13]進行了對比。圖8表示了深拉成形后沿軋制方向和橫向的最大與最小主應變模擬結果及其與實驗數據對比。圖8中，橫坐標表示深拉后不同區域距離深拉杯底部中心(以Pole表示)的距離，a區域為底部中心，b區域為沖頭圓角處，c區域為側壁區域，d區域為深拉杯的頂部區域；RD Major Strain-Experiment 表示實驗測得的沿軋制方向的最大主應變分布，TD Minor Strain-CPB-GTN 表示采用CPB-GTN子程序模擬得到的沿橫向的最小主應變分布，其余以此類推。由圖8可以看出，模擬的最大與最小主應變同試驗結果趨勢基本一致，只是在橫向最大主應變存在一定誤差，這可能與較小應變測量時各向異性效應的忽略和較小應變測量的不確定性有關。

模擬預測的沖頭力與位移關系曲線及其與實驗結果對比如圖9所示。由圖9可以看出，模擬預測的沖頭力與位移關系曲線的前半部分與實驗結果吻合得很好，后半部分雖然存在一定的誤差，但與沒有考慮損傷的基于CPB06ex3ev屈服本構模型的模擬結果[13]相比，本節的非等溫深拉成形模擬結果總體上比較精確，說明了考慮參數演化的各向異性屈服CPB06-GTN損傷模型子程序能夠用于模擬鎂合金的非等溫成形過程。

圖8 沿軋制方向和橫向的最大與最小主應變分布模擬預測及其與實驗結果[13]的對比Fig.8 Simulation-based predictions of distribution of the maximum and minimum principal strain in rolling and transverse directions and their comparisons with experimental results[13]

圖9 沖頭力與位移關系曲線的模擬預測及其與實驗結果[13]對比Fig.9 Simulation-based prediction of punch force vs displacement curves and its comparison with experimental results [13]

4 結 論

1) 在GTN損傷模型的基礎上，建立了考慮參數演化的各向異性拉壓不對稱屈服的CPB06-GTN損傷耦合模型，該模型不僅能夠預測材料的損傷破壞，而且還能準確描述鎂合金等HCP晶格結構材料的各向異性拉壓不對稱屈服及其在畸變硬化過程中的各向異性屈服面變化趨勢。

2) 基于建立的鎂合金板熱成形各向異性屈服及損傷本構方程，同時考慮成形時的熱、力耦合關系，編寫了VUMAT子程序并嵌入ABAQUS/Explicit分析中，使用該子程序進行的等溫及非等溫成形數值模擬算例證明了所建立本構模型的正確性。

3) 同鎂合金等溫拉深成形相比，非等溫拉深成形模擬結果表明，非等溫拉深成形能夠改善鎂合金的成形性，提高鎂合金的成形質量。

參 考 文 獻

[1] 宋令慧, 王守仁, 趙宰炯. 連鑄連軋鎂合金AZ41微觀結構與摩擦磨損性能[J]. 航空學報, 2014, 35(6): 1733-1739.

SONG L H, WANG S R, CHO J H. Microstructure and friction wear properties of twin-roll casting magneisum alloy AZ41[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(6): 1733-1739 (in Chinese).

[2] 唐荻, 王丹, 江海濤, 等. 異步軋制對 AZ31 鎂合金組織及高溫塑性的影響[J]. 塑性工程學報, 2013, 20(2): 78-83.

TANG D, WANG D, JIANG H T, et al. Effects of differential speed rolling process on microstructure and plastic deformation of magnesium alloy[J]. Journal of Plasticity Engineering, 2013, 20(2): 78-83 (in Chinese).

[3] 吳章斌, 桂良進, 范子杰. AZ31B 鎂合金擠壓材料的力學性能與本構分析[J]. 中國有色金屬學報, 2015, 25(2): 293-300.

WU Z B, GUI L J, FAN Z J. Mechanical properties and constitutive analysis of extruded AZ31B magnesium alloy[J]. Chinese Journal of Nonferrous Metals, 2015, 25(2): 293-300 (in Chinese).

[4] 鐘敏, 唐偉琴, 李大永, 等. AZ31B鎂合金板材溫熱成形極限實驗研究[J]. 塑性工程學報, 2011, 18(5): 59-63.

ZHONG M, TANG W Q, LI D Y, et al. Experiment research on forming limit diagram (FLD) of magnesium alloy sheet AZ31B at warm condition[J]. Journal of Plasticity Engineering, 2011, 18(5): 59-63 (in Chinese).

[5] 何維均, 張士宏, 程明, 等. 宏觀彈塑性本構模型的研究進展[J]. 塑性工程學報, 2015, 22(3): 1-11.

HE W J, ZHANG S H, CHENG M, et al. Review on the development of macroscopic elastic-plastic constitutive models[J]. Journal of Plasticity Engineering, 2015, 22(3): 1-11 (in Chinese).

[6] HILL R. A theory of the yielding and plastic flow of anisotropic metals[C]∥Proceedings of the Royal Society of London Series A, Mathematical and Engineering Sciences. London: The Royal Society Publishing, 1948, 193(1033): 281-297.

[7] BARLAT F, LEGE D J, BREM J C. A six-component yield function for anisotropic materials[J]. International Journal of Plasticity, 1991, 7(7): 693-712.

[8] BARLAT F, MAEDA Y, CHUNG K, et al. Yield function development for aluminum alloy sheets[J]. Journal of the Mechanics and Physics of Solids, 1997, 45(11): 1727-1763.

[9] BARLAT F, BREM J C, YOON J W, et al. Plane stress yield function for aluminum alloy sheets—Part 1: Theory[J]. International Journal of Plasticity, 2003, 19(9): 1297-1319.

[10] HOSFORD W F. Incorporating work hardening in yield loci calculations[C]∥Strength of Metals and Alloys. Aachen: Federal Republic of Germany, 1979, 775-780.

[11] CAZACU O, PLUNKETT B, BARLAT F. Orthotropic yield criterion for hexagonal closed packed metals[J]. International Journal of Plasticity, 2006, 22(7): 1171-1194.

[12] PLUNKETT B, CAZACU O, BARLAT F. Orthotropic yield criteria for description of the anisotropy in tension and compression of sheet metals[J]. International Journal of Plasticity, 2008, 24(5): 847-866.

[13] TARI D G, WORSWICK M J. Elevated temperature constitutive behavior and simulation of warm forming of AZ31B[J]. Journal of Materials Processing Technology, 2015, 221: 40-55.

[14] 郎利輝, 楊希英, 劉康寧, 等. 一種韌性斷裂準則中材料常數的計算模型及其應用[J]. 航空學報, 2015, 36(2): 672-679.

LANG L H, YANG X Y, LIU K N, et al. A calculating model of material constants in ductile fracture criterion and its applications[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(2): 672-679 (in Chinese).

[15] TENG B G, YUAN S J, CHEN Z T, et al. Plastic damage of T-shape hydroforming[J]. Transactions of Nonferrous Metals Society of China, 2012, 22(S2): s294-s301.

[16] 王瑞澤, 陳章華, 臧勇. 基于Gurson 模型的鎂合金板材溫熱沖壓成形研究[J]. 北京科技大學學報, 2014, 36(4): 459-466.

WANG R Z, CHEN Z H, ZANG Y. Thermal stamping formability of magnesium alloy sheet based on the Gurson model[J]. Journal of University of Science and Technology Beijing, 2014, 36(4): 459-466 (in Chinese).

[17] GURSON A L. Continuum theory of ductile rupture by void nucleation and growth: Part 1—Yield criteria and flow rules for porous ductile media[J]. Journal of Engineering Materials and Technology, 1977, 99(1): 2-15.

[18] NEEDLEMAN A, TVERGAARD V. An analysis of dynamic, ductile crack growth in a double edge cracked specimen[J]. International Journal of Fracture, 1991, 49(1): 41-67.

[19] KIM J, RYOU H, KIM D, et al. Constitutive law for AZ31B Mg alloy sheets and finite element simulation for three-point bending[J]. International Journal of Mechanical Sciences, 2008, 50(10): 1510-1518.

[20] STEWART J B, CAZACU O. Analytical yield criterion for an anisotropic material containing spherical voids and exhibiting tension-compression asymmetry[J]. International Journal of Solids and Structures, 2011, 48(2): 357-373.

[21] YOON J, CAZACU O, MISHRA R K. Constitutive modeling of AZ31 sheet alloy with application to axial crushing[J]. Materials Science and Engineering: A, 2013, 565: 203-212.

[22] WANG R, CHEN Z, LI Y, et al. Failure analysis of AZ31 magnesium alloy sheets based on the extended GTN damage model[J]. International Journal of Minerals, Metallurgy, and Materials, 2013, 20(12): 1198-1207.

[23] 陳志英, 董湘懷. 各向異性 GTN 損傷模型及其在板料成形中的應用[J]. 上海交通大學學報, 2008, 42(9): 1414-1419.

CHEN Z Y, DONG X H. The Anisotropic GTN Damage model and its application in sheet metal forming[J]. Journal of Shanghai Jiaotong University, 2008, 42(9): 1414-1419 (in Chinese).

[24] TARI D G, WORSWICK M J, WINKLER S. Experimental studies of deep drawing of AZ31B magnesium alloy sheet under various thermal conditions[J]. Journal of Materials Processing Technolog, 2013, 213(8): 1337-1347.

[25] LI W J, ZHAO G Q, MA X W, et al. Study on forming limit diagrams of AZ31B alloy sheet at different temperatures[J]. Materials & Manufacturing Processes, 2013, 28(3): 306-311.

[26] LEE S, HAM H J, KWON S Y, et al. Thermal conductivity of magnesium alloys in the temperature range from -125 ℃ to 400 ℃[J]. International Journal of Thermophysics, 2013: 34(12): 2343-2350.