基于全耦合損傷模型的數(shù)值模擬研究

張杰

（上海圣堯智能科技有限公司 上海市 200241）

1 緒論

金屬成形和加工過(guò)程的有限元模擬對(duì)工業(yè)優(yōu)化和設(shè)計(jì)高質(zhì)量產(chǎn)品具有重大意義，它可以利用計(jì)算機(jī)虛擬地研究加工工藝參數(shù)、成形工具，甚至材料力學(xué)行為[1,2]。從而大量地降低生產(chǎn)成本，有效地縮短設(shè)計(jì)周期。目前，它已被應(yīng)用于力、電、熱和流體等各個(gè)領(lǐng)域[3]。它是一種求偏微分方程及其方程組近似解的數(shù)值技術(shù)。簡(jiǎn)而言之，有限元分析是把一個(gè)非常復(fù)雜的問(wèn)題分解成若干便于求解的小單元的一種方法。由于有限元法可以為產(chǎn)品、工具、機(jī)器和工藝的設(shè)計(jì)提供詳細(xì)的解決方案，有限元方法已經(jīng)在金屬成形和加工過(guò)程中非常有效的工具。

然而，如何準(zhǔn)確地預(yù)測(cè)材料的行為，特別是塑性損傷，對(duì)工程技術(shù)人員來(lái)說(shuō)是一個(gè)巨大的挑戰(zhàn)[4]。建立合適的損傷演化模型，有利于避免金屬成形過(guò)程中工件的損傷產(chǎn)生或者有利于金屬加工過(guò)程中工件和廢料的分離。另一個(gè)挑戰(zhàn)來(lái)自于金屬成形和加工過(guò)程中復(fù)雜的應(yīng)變路徑和環(huán)境條件。例如，在不同的環(huán)境溫度或加載速度下，耦合損傷行為的彈塑性響應(yīng)變化也不同。同時(shí)，在不同的應(yīng)變路徑下?lián)p傷演化具有多樣性。因此，準(zhǔn)確描述損傷演化過(guò)程是提高金屬成形和加工精確度的關(guān)鍵。

金屬損傷力學(xué)的本質(zhì)是材料由于微缺陷（孔洞和裂紋）的形核和長(zhǎng)大而退化，并聚結(jié)成宏觀裂紋[5]。McClintock[6]首先提出了微缺陷與延性損傷的關(guān)系。之后，基于微缺陷[7]的三種主要方法被廣泛用于描述損傷力學(xué)：斷裂力學(xué)、細(xì)觀損傷力學(xué)和連續(xù)損傷力學(xué)(CDM)。這種全耦合的方法解釋了塑性流動(dòng)(包括不同類(lèi)型的硬化)與延性損傷的出現(xiàn)和生長(zhǎng)之間的相互作用。在連續(xù)損傷力學(xué)中，假定損傷是材料的一種內(nèi)部狀態(tài)變量，該狀態(tài)變量和材料微觀結(jié)構(gòu)不可逆退化有關(guān)，且損傷變量的變化與等效應(yīng)力有關(guān)。Kachanov[8]引入了一個(gè)連續(xù)度的概念來(lái)描述金屬在單軸載荷作用下的蠕變斷裂。將微裂紋引起的截面面積減小量作為內(nèi)部損傷狀態(tài)的一種表述方法，對(duì)損傷變量的表述帶來(lái)重要物理意義。Lemaitre和Chaboche 開(kāi)發(fā)了延性損傷的連續(xù)損傷力學(xué)模型[9]，通過(guò)使用等效狀態(tài)變量對(duì)損傷勢(shì)能求導(dǎo)獲得損傷的本構(gòu)關(guān)系。近年來(lái)，大量關(guān)于彈性、塑性和熱粘塑性與塑性損傷耦合的本構(gòu)模型相繼被提出。Bouchard[10]、Brünig[11]和Wang[12]等人總結(jié)并比較了各種損傷模型，發(fā)現(xiàn)Lemaitre 損傷模型可以較容易地分析各種力學(xué)應(yīng)用中的耦合損傷。

Johnson-Cook（JC）硬化模型是最常用的粘塑性應(yīng)變流動(dòng)準(zhǔn)則。該模型既考慮了材料變形過(guò)程中的隨動(dòng)硬化作用，又考慮了絕熱加熱作用，能夠描述材料在不同熱環(huán)境下的動(dòng)態(tài)行為。Peir[13]采用先進(jìn)的實(shí)驗(yàn)方法和有限元模擬工具，對(duì)JC 模型中的材料參數(shù)，尤其是應(yīng)變速率硬化和熱軟化參數(shù)進(jìn)行了驗(yàn)證。通過(guò)提高熱軟化影響，用修正參數(shù)進(jìn)行的模擬結(jié)果與實(shí)驗(yàn)結(jié)果一致。Calamaz[14]直接將JC模型轉(zhuǎn)變?yōu)門(mén)ANH 模型，模擬了鋸齒形切屑在正交切削過(guò)程中的形成。另一方面，Zerilli 和Armstrong 基于位錯(cuò)機(jī)理提出了不同晶體結(jié)構(gòu)的本構(gòu)關(guān)系，并在熱活化分析的基礎(chǔ)上，將應(yīng)變硬化、應(yīng)變率硬化和熱軟化的影響耦合進(jìn)本構(gòu)關(guān)系中[15-17]。Holmquist[18]和Hor 等[19]對(duì)不依賴(lài)于材料常數(shù)的模型進(jìn)行了比較分析，提出了使用相同的試驗(yàn)數(shù)據(jù)確定不同本構(gòu)模型常數(shù)的方法

表1：鈦合金拉伸試驗(yàn)材料參數(shù)

表2：鈦合金 TC4 剪切試驗(yàn)的材料參數(shù)

表3：TC4 鈦合金剪切試驗(yàn)的自適應(yīng)網(wǎng)格重構(gòu)參數(shù)

本文提出一種將延性損傷耦合到JC 流動(dòng)準(zhǔn)則中的本構(gòu)模型用于描述材料的變形行為，并通過(guò)Abaqus/Explicit 以及Abaqus 的子程序Vumat 來(lái)實(shí)現(xiàn)材料變形的仿真過(guò)程。通過(guò)拉伸和剪切試驗(yàn)的仿真對(duì)該全耦合模型進(jìn)行驗(yàn)證，證明了該模型在模擬金屬材料成形和加工過(guò)程中的可靠性。

2 全耦合損傷模型構(gòu)建

2.1 熱-粘-彈塑性硬化模型

JC 模型作為一種在熱力學(xué)領(lǐng)域廣泛應(yīng)用的經(jīng)驗(yàn)型模型，已成功地用于模擬金屬在不同溫度和應(yīng)變速率下的行為，流動(dòng)準(zhǔn)則如下：

該模型提供了直接計(jì)算塑料應(yīng)變率的方法，它考慮應(yīng)變率與屈服應(yīng)力的關(guān)系，并使用與速率無(wú)關(guān)的塑性來(lái)計(jì)算材料的粘塑性響應(yīng)。同時(shí)，在屈服應(yīng)力中引入溫度軟化效應(yīng)?；谠隽坷碚摚?jì)算每個(gè)增量下的溫度增量和應(yīng)變速率。在本文中，溫度的計(jì)算采用絕熱分析，塑性應(yīng)變引起的溫度增量如下：

2.2 狀態(tài)變量方程構(gòu)建

為了將損傷行為耦合到 JC 模型中，引入了損傷狀態(tài)變量參數(shù)D(0 ≤D ≤1)。該損傷變量是由微裂紋與微孔洞引起材料失效的面積與承載外力的實(shí)際面積的比值進(jìn)行宏觀表示。為了表征該損傷變量對(duì)力學(xué)行為的影響，引入了等效狀態(tài)變量本文根據(jù)能量不變理論，考慮損傷的等效狀態(tài)變量定義如下：

材料剛度計(jì)算如下：

2.3 耗散勢(shì)能分析

為了耦合損傷行為，控制應(yīng)力空間內(nèi)部變量的演化規(guī)律，定義勢(shì)能方程如下：

塑性勢(shì)能f 和損傷勢(shì)能FY計(jì)算方式如式（9）：

圖1：板材拉伸的試驗(yàn)和仿真應(yīng)力-應(yīng)變曲線

圖2：耦合模型中的損傷演化過(guò)程

圖3：試樣尺寸和仿真模型初始網(wǎng)格

耗散勢(shì)能可以寫(xiě)為：

圖4：板材拉伸過(guò)程中損傷局部化、裂紋的形成和擴(kuò)展過(guò)程

圖5：板材拉伸仿真結(jié)果與試驗(yàn)結(jié)果的力-位移曲線對(duì)比

塑性修正系數(shù) 計(jì)算如下：

圖6：板材拉伸仿真結(jié)果與試驗(yàn)結(jié)果的應(yīng)變分布對(duì)比

圖7：剪切試件形狀及尺寸

圖8：板材剪切試驗(yàn)和仿真應(yīng)力-應(yīng)變曲線

2.4 數(shù)值求解方案

在本文中，認(rèn)為熱-粘-彈塑性區(qū)位于空間區(qū)域 V 中?；谒苄院哪芾碚摵吞摴υ?，當(dāng)前損傷區(qū)域表示為：

通過(guò)使用節(jié)點(diǎn)近似法，定義了一個(gè)高度非線性系統(tǒng)描述每一個(gè)分析步中部件的動(dòng)態(tài)平衡。該系統(tǒng)既可以通過(guò)靜態(tài)隱式迭代法求解，也可以使用動(dòng)態(tài)顯示進(jìn)行求解。隱式迭代法求解過(guò)程中要求有一個(gè)恒定的剛度矩陣，然而，由于損傷導(dǎo)致材料軟化行為的存在，很難確定一個(gè)恒定的剛度矩陣。與此相反的，動(dòng)態(tài)顯示方法可以利用切向矩陣直接求解線性代數(shù)系統(tǒng)，避免了迭代。但是，為了滿足穩(wěn)定條件，時(shí)間增量步的長(zhǎng)度并不固定，根據(jù)需要自動(dòng)調(diào)節(jié)。

本文提出的全耦合本構(gòu)方程，使用了彈性預(yù)測(cè)-返回映射算法和算子分割方法計(jì)算應(yīng)力、塑性應(yīng)變、累積塑性應(yīng)變和損傷變量等。利用Newton-Raphson 迭代法求解公式（17）以及未知變量

在本文中，有兩種求解方法可以使用：強(qiáng)耦合和弱耦合。強(qiáng)耦合同時(shí)求解公式（17）中的方程(a)和(b)，獲得兩個(gè)未知變量。弱耦合假設(shè)因此，只需要求解方程(a)，求得變量用于更新?tīng)顟B(tài)變量。Bouchard[10]對(duì)使用強(qiáng)耦合和弱耦合方法的等效塑料應(yīng)變和損傷變量進(jìn)行了對(duì)比研究。研究發(fā)現(xiàn)，強(qiáng)耦合方法可以獲得更高的精度，但是弱耦合的精度也可以接受。因此，為了降低計(jì)算成本，本文選擇了弱耦合求解方法。同時(shí)，為了簡(jiǎn)化材料模型，假設(shè)摩擦定律不受損傷的影響。

在有限元仿真中，設(shè)置損傷的最大值Dmax=0.9999，如果損傷變量D 達(dá)到Dmax，即認(rèn)為該單元失效，不再對(duì)剛度矩陣產(chǎn)生作用。同時(shí)，刪除失效的單元，并生成新的邊界和網(wǎng)格。

3 數(shù)值仿真與驗(yàn)證

本節(jié)將通過(guò)鈦合金TC4 板材的拉伸、剪切試驗(yàn)及仿真來(lái)驗(yàn)證本文提出的全耦合損傷模型。本文通過(guò)ABAQUS/Explicit 求解器和用戶子程序Vumat 來(lái)實(shí)現(xiàn)該模型的仿真驗(yàn)證。

3.1 拉伸試驗(yàn)

利用本文提出的全耦合本構(gòu)模型對(duì)鈦合金TC4 的拉伸應(yīng)力-應(yīng)變曲線進(jìn)行數(shù)值預(yù)測(cè)。用簡(jiǎn)單拉伸試驗(yàn)到的實(shí)驗(yàn)硬化曲線來(lái)確定本構(gòu)方程的材料常數(shù)，如表1 所示。試驗(yàn)和仿真結(jié)果如圖1 所示，當(dāng)?shù)刃苄詰?yīng)變小于0.25 時(shí)，耦合材料模型結(jié)果和實(shí)驗(yàn)結(jié)果能夠較好的吻合，當(dāng)?shù)刃苄詰?yīng)變達(dá)到0.355 時(shí)，材料剛度迅速退化為0，即材料完全失效，材料的損傷演化過(guò)程如圖2 所示。分別對(duì)耦合損傷模型和無(wú)損傷模型進(jìn)行仿真研究。研究發(fā)現(xiàn)，在無(wú)損傷的模型中，屈服應(yīng)力單調(diào)遞增且無(wú)退化行為，如圖1 所示。

拉伸試驗(yàn)采用狗骨形狀的試件，標(biāo)距為5.0 mm，寬4.0 mm，過(guò)渡區(qū)的半徑為1.0 mm，試樣厚度為0.6 mm，如圖3 所示。

試件以10 mm/s 的恒定速度進(jìn)行拉伸，通過(guò)大量的小位移逐步進(jìn)行加載。隨著時(shí)間的累積，損傷局部化、裂紋形成以及最終的斷裂相繼出現(xiàn)。最初，損傷出現(xiàn)在試件中心，如圖4(a)所示；當(dāng)損傷達(dá)到細(xì)化閾值后，網(wǎng)格進(jìn)行重構(gòu)細(xì)化，如圖4(b)所示；當(dāng)損傷累積到一定程度之后，剪切帶出現(xiàn)，并沿45°方向從試件中心向兩邊擴(kuò)展，如圖4(c)所示；隨著加載的進(jìn)一步增加，試件中心出現(xiàn)裂紋并沿剪切帶擴(kuò)展，如圖4(d)所示；最終，當(dāng)拉伸位移到0.918 mm 時(shí)，試件完全斷裂，如圖4(e)所示。試件斷裂面的損傷分布以及等效塑性應(yīng)變?nèi)鐖D4(e)和(f)所示。在拉伸試驗(yàn)中，試件接近斷裂時(shí)，變形局部化是一種正?，F(xiàn)象。該模擬過(guò)程清晰地描述了損傷局部化以及裂紋擴(kuò)展現(xiàn)象。結(jié)果表明，該模型能夠較好地模擬損傷局部化后宏觀裂紋的萌生和擴(kuò)展過(guò)程。

圖9：剪切試驗(yàn)中損傷局部化、裂紋的形成和擴(kuò)展過(guò)程

圖10：仿真和試驗(yàn)的力-位移曲線對(duì)比分析

為了進(jìn)一步驗(yàn)證該模型，將仿真結(jié)果的拉伸力與試驗(yàn)值進(jìn)行對(duì)比分析。圖5 展示了恒定拉伸速度下全耦合模型仿真結(jié)果與試驗(yàn)結(jié)果的力-位移曲線對(duì)比情況。由圖可以發(fā)現(xiàn)，由于損傷引起的軟化現(xiàn)象的存在使得試件在拉伸距離為0.91 mm 時(shí)最終斷裂；在拉伸距離0.4 mm 時(shí)，拉伸力達(dá)到最大值Fmax=2823N。除此之外，本文還通過(guò)仿真結(jié)果與試驗(yàn)結(jié)果中的應(yīng)變分布的對(duì)比分析對(duì)模型進(jìn)行驗(yàn)證，試驗(yàn)中的應(yīng)變分布通過(guò)DIC 技術(shù)獲得，如圖6 所示。由圖6 可以發(fā)現(xiàn)，仿真中的應(yīng)變分布和試驗(yàn)中的應(yīng)變分布能夠較好的吻合，且試件的等效應(yīng)變沿軸向和橫向分布都不均勻。隨著損傷的累積和局部化，應(yīng)變也首先出現(xiàn)在試件的中心，然后沿剪切帶集中。

3.2 剪切試驗(yàn)

在上面的拉伸試驗(yàn)仿真中，損傷出現(xiàn)和擴(kuò)展過(guò)程中已經(jīng)可以觀察到剪切帶的出現(xiàn)。相比于拉伸試驗(yàn)，剪切試驗(yàn)中沒(méi)有截面的縮減，材料變形也是純剪切變形。剪切變形中較低的應(yīng)力三軸度可以有效降低損傷累積速率，因此，在剪切變形中可以發(fā)生較大的應(yīng)變。剪切試件的幾何形狀和尺寸如圖7 所示，剪切區(qū)在試件的中心。本仿真采用與拉伸相同的鈦合金板材，新的JC 參數(shù)和損傷參數(shù)如表2所示。仿真結(jié)果與試驗(yàn)結(jié)果的應(yīng)力-應(yīng)變曲線對(duì)比如圖8 所示，在等效塑性應(yīng)變小于0.21 時(shí)，耦合模型仿真結(jié)果和試驗(yàn)結(jié)果能夠較好的吻合，當(dāng)塑性應(yīng)變達(dá)到0.297 時(shí)，材料迅速失效。然而，未耦合損傷的JC 模型，屈服應(yīng)力單調(diào)遞增且沒(méi)有退化現(xiàn)象。

與拉伸試驗(yàn)相同，剪切試驗(yàn)的仿真也以10 mm/s 的恒定速度進(jìn)行，通過(guò)大量的小位移進(jìn)行逐步加載。隨著分析時(shí)間的增加，損傷的局部集中、裂紋的產(chǎn)生、試件的斷裂依次發(fā)生。如表3 所示。

當(dāng)剪切位移達(dá)到0.40 mm 時(shí)，損傷出現(xiàn)，但出現(xiàn)在剪切區(qū)的兩邊并不是中心，如圖9(a)所示；當(dāng)位移增加到0.60 mm 時(shí)，損傷在試件中心集中，如圖9(b)所示；當(dāng)位移從0.60 mm 逐漸增加到0.72 mm 時(shí)，損傷快速增長(zhǎng)，如圖9(c)所示；當(dāng)位移從0.72 mm 增加到0.74 mm 時(shí)，試件中心出現(xiàn)裂紋，并沿剪切帶向兩端擴(kuò)展，如圖9(d)所示；最終，當(dāng)位移達(dá)到0.75 mm 時(shí)，試件完全斷裂。斷裂后的損傷分布如圖9(e)和(f)所示。該仿真清楚的描述了剪切帶的形成和裂紋擴(kuò)展的過(guò)程，剪切帶和裂紋擴(kuò)展方向在純剪切過(guò)程中是一條直線，和拉伸時(shí)有明顯的不同。

為了進(jìn)一步驗(yàn)證該模型的有效性，對(duì)剪切過(guò)程中的力-位移曲線進(jìn)行了對(duì)比分析，結(jié)果如圖10 所示。由圖可知，在拉伸距離為0.75 mm 時(shí)材料完全失效，拉伸距離0.68 mm 時(shí)，達(dá)到最大剪切力Fmax=1504N。不同位移時(shí)的等效應(yīng)力分布也展示在圖10 中，結(jié)果表明，在達(dá)到最終斷裂時(shí)，剪切區(qū)的材料已經(jīng)失去了承載能力。

4 結(jié)論

本文基于介質(zhì)力學(xué)和全耦合理論，利用JC 模型，建立了一種考慮熱-粘-彈塑性的全耦合損傷本構(gòu)模型，用于金屬成形和加工過(guò)程數(shù)值模擬。該模型能夠準(zhǔn)確預(yù)測(cè)延性損傷的發(fā)生以及描述熱軟化和應(yīng)變速率硬化的影響。本文通過(guò)拉伸和剪切試驗(yàn)的仿真與試驗(yàn)結(jié)果的對(duì)比分析驗(yàn)證了該模型的有效性，較好的模擬了損傷行為的演化過(guò)程，包括：損傷局部化、微裂紋的出現(xiàn)及擴(kuò)展等。因此，使用該全耦合損傷模型可以較好的模擬金屬成形及加工過(guò)程的各種問(wèn)題。