AZ31鎂合金板熱成形中的屈服和損傷: 本構(gòu)實(shí)現(xiàn)與數(shù)值分析

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

1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室, 大連 116024 2.大連理工大學(xué) 工程力學(xué)系, 大連 116024 3.國際計(jì)算力學(xué) 研究中心, 大連 116024

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

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

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

1 本構(gòu)模型

1.1 CPB06各向異性拉壓不對稱屈服本構(gòu)

在平面應(yīng)力條件下，CPB06本構(gòu)模型的具體函數(shù)形式為[19]

(1)

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

(2)

(3)

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

從以上結(jié)果推導(dǎo)可得出，當(dāng)k=0，C為4階單位矩陣，且a=2時(shí)，該模型與Mises屈服準(zhǔn)則等效。

1.2 GTN損傷本構(gòu)

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

(1+q3f*2)=0

(4)

材料進(jìn)入塑性時(shí)由流動(dòng)法則可計(jì)算得到塑性應(yīng)變增量為

(5)

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

(6)

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

(7)

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

(8)

(9)

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

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

(10)

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

(11)

(12)

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

(13)

(14)

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

(15)

1.3 修正的CPB06-GTN本構(gòu)模型

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

(1+q3f*2)=0

(16)

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

σq=m·

(17)

式中：

m=

式(16)中：

(18)

(19)

圖1 應(yīng)變率為和不同溫度(T=150, 200, 250 ℃)時(shí)基于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 ℃時(shí)各向異性CPB06屈服本構(gòu)參數(shù)隨等效塑性應(yīng)變的變化Table 1 Variation of anisotropic CPB06 yield criterion parameters with effective plastic strain at 150 ℃

表2 200 ℃時(shí)各向異性CPB06屈服本構(gòu)參數(shù)隨等效塑性應(yīng)變的變化Table 2 Variation of anisotropic CPB06 yield criterion parameters with effective plastic strain at 200 ℃

表3 250 ℃時(shí)CPB06各向異性屈服本構(gòu)參數(shù)隨等效塑性應(yīng)變的變化Table 3 Variation of anisotropic CPB06 yield criterion parameters with effective plastic strain at 250 ℃

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

(20)

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

Table4HardeningparametersofAZ31Bmagnesiumalloyatdifferenttemperatures

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

(21)

2 數(shù)值算法實(shí)現(xiàn)及驗(yàn)證

2.1 VUMA子程序?qū)崿F(xiàn)

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

步驟2假設(shè)應(yīng)變增量為完全彈性，計(jì)算得到彈性預(yù)測應(yīng)力：

(22)

(23)

步驟4塑性修正，通過迭代求解塑性應(yīng)變確定流動(dòng)方向n：

(24)

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

(25)

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

(26)

更新相應(yīng)的等效應(yīng)力σq和靜水應(yīng)力σp：

(27)

(28)

(29)

式中：左上標(biāo)k和k+1 表示迭代次數(shù)；當(dāng)孔洞體積率f>fc時(shí)，等效孔洞體積分?jǐn)?shù)f*=fc+δ(f-fc)。

結(jié)合硬化準(zhǔn)則對基體材料的屈服應(yīng)力進(jìn)行更新：

(30)

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

步驟5完成應(yīng)力更新，判斷孔洞體積率是否達(dá)到臨界值。

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

(31)

通過對預(yù)測應(yīng)力進(jìn)行修正即可得到當(dāng)前狀態(tài)的真實(shí)應(yīng)力：

σ=σtr-Ce:Δεp

(32)

式中：σtr為彈性預(yù)測應(yīng)力；Ce為4階彈性模量；Δεp為塑性應(yīng)變增量；當(dāng)孔洞體積率達(dá)到臨界值時(shí)，則相應(yīng)的控制單元?jiǎng)h除的狀態(tài)變量等于0，意味著單元失效。

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

(33)

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

(34)

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

綜上，可得出厚度方向的應(yīng)變增量為

dεq(-n11-n22)

(35)

步驟7結(jié)束當(dāng)前時(shí)間步，進(jìn)入下一時(shí)間步。

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

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

2.2 單個(gè)單元測試

圖3 采用(VUMAT)子程序模擬的單軸拉伸和壓縮應(yīng)力應(yīng)變曲線及其與實(shí)驗(yàn)數(shù)據(jù)[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 ℃等溫成形模擬

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

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

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

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

圖5 200 ℃等溫成形后帶裂紋深拉杯厚度分布模擬結(jié)果及其與成形極限圖的對比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 非等溫成形熱力耦合模擬

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

在有限元模擬時(shí)，通過在模具上施加溫度邊界條件來控制模具的溫度，鎂合金板的熱參數(shù)及與模具間的熱傳導(dǎo)率參數(shù)采用Lee等[26]給出的AZ31鎂合金在200 ℃的熱參數(shù)，具體數(shù)據(jù)如表5所示。變形過程中模具和鎂合金板之間會(huì)發(fā)生熱傳遞，在模擬時(shí)為了節(jié)省計(jì)算時(shí)間，提高了拉深速度；在溫?zé)岢尚螘r(shí)為了讓模擬的溫度與速度相匹配，也對鎂合金與時(shí)間相關(guān)的材料參數(shù)進(jìn)行了放大處理。由于對計(jì)算速度放大了100倍，這里對熱傳導(dǎo)率也放大了100倍。參考實(shí)驗(yàn)條件，對直徑為228.6 mm的鎂合金板，在對應(yīng)拉深比為2.25，壓邊力為80 kN時(shí)的溫?zé)嵘罾堑葴爻尚?，進(jìn)行了數(shù)值模擬。

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

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

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

此外，還將非等溫深拉成形模擬結(jié)果與已有實(shí)驗(yàn)數(shù)據(jù)[13]進(jìn)行了對比。圖8表示了深拉成形后沿軋制方向和橫向的最大與最小主應(yīng)變模擬結(jié)果及其與實(shí)驗(yàn)數(shù)據(jù)對比。圖8中，橫坐標(biāo)表示深拉后不同區(qū)域距離深拉杯底部中心(以Pole表示)的距離，a區(qū)域?yàn)榈撞恐行?，b區(qū)域?yàn)闆_頭圓角處，c區(qū)域?yàn)閭?cè)壁區(qū)域，d區(qū)域?yàn)樯罾捻敳繀^(qū)域；RD Major Strain-Experiment 表示實(shí)驗(yàn)測得的沿軋制方向的最大主應(yīng)變分布，TD Minor Strain-CPB-GTN 表示采用CPB-GTN子程序模擬得到的沿橫向的最小主應(yīng)變分布，其余以此類推。由圖8可以看出，模擬的最大與最小主應(yīng)變同試驗(yàn)結(jié)果趨勢基本一致，只是在橫向最大主應(yīng)變存在一定誤差，這可能與較小應(yīng)變測量時(shí)各向異性效應(yīng)的忽略和較小應(yīng)變測量的不確定性有關(guān)。

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

圖8 沿軋制方向和橫向的最大與最小主應(yīng)變分布模擬預(yù)測及其與實(shí)驗(yàn)結(jié)果[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 沖頭力與位移關(guān)系曲線的模擬預(yù)測及其與實(shí)驗(yàn)結(jié)果[13]對比Fig.9 Simulation-based prediction of punch force vs displacement curves and its comparison with experimental results [13]

4 結(jié) 論

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

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

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

參 考 文 獻(xiàn)

[1] 宋令慧, 王守仁, 趙宰炯. 連鑄連軋鎂合金AZ41微觀結(jié)構(gòu)與摩擦磨損性能[J]. 航空學(xué)報(bào), 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]. 塑性工程學(xué)報(bào), 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] 吳章斌, 桂良進(jìn), 范子杰. AZ31B 鎂合金擠壓材料的力學(xué)性能與本構(gòu)分析[J]. 中國有色金屬學(xué)報(bào), 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鎂合金板材溫?zé)岢尚螛O限實(shí)驗(yàn)研究[J]. 塑性工程學(xué)報(bào), 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] 何維均, 張士宏, 程明, 等. 宏觀彈塑性本構(gòu)模型的研究進(jìn)展[J]. 塑性工程學(xué)報(bào), 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] 郎利輝, 楊希英, 劉康寧, 等. 一種韌性斷裂準(zhǔn)則中材料常數(shù)的計(jì)算模型及其應(yīng)用[J]. 航空學(xué)報(bào), 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 模型的鎂合金板材溫?zé)釠_壓成形研究[J]. 北京科技大學(xué)學(xué)報(bào), 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 損傷模型及其在板料成形中的應(yīng)用[J]. 上海交通大學(xué)學(xué)報(bào), 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.