王 剛， 宋 建， 張建偉， 王冰冰

(1.鄭州大學(xué) 機(jī)械與動力工程學(xué)院，河南 鄭州 450001； 2.鄭州大學(xué) 抗疲勞制造技術(shù)河南省工程實(shí)驗(yàn)室，河南 鄭州 450001； 3.鄭州大學(xué) 力學(xué)與工程科學(xué)學(xué)院，河南 鄭州 450001)

0 引言

在長期的研究工作中，人們逐漸認(rèn)識到材料的性能以及一些關(guān)鍵構(gòu)件的疲勞失效及斷裂不僅與材料的化學(xué)成分有關(guān)，還在很大程度上取決于材料的微結(jié)構(gòu)。Huang等[1]通過表面軋制處理，在AISI 316 L不銹鋼中生成了梯度納米層結(jié)構(gòu)，在對其疲勞性能的研究中發(fā)現(xiàn)，疲勞裂紋在表面下500 μm深處萌生，并向表面和內(nèi)部擴(kuò)展，最終導(dǎo)致試樣瞬間斷裂。大多數(shù)金屬材料通常是晶體的集合體，即多晶體，在宏觀上表現(xiàn)為各向同性；而在微觀上，晶粒的形狀、大小、晶體取向都不相同，表現(xiàn)為各向異性?？梢娊⒕哂形⒂^組織結(jié)構(gòu)的多晶體模型對于模擬材料性能和疲勞性能也變得越來越重要。

將晶體塑性有限元方法與隨機(jī)的晶粒結(jié)構(gòu)、晶體取向等結(jié)合起來，是目前在晶粒尺寸水平上研究微結(jié)構(gòu)局部變形、應(yīng)力變化的主要手段。微結(jié)構(gòu)有限元模型的建立通常采用Voronoi圖原理的方法，或者采用EBSD技術(shù)進(jìn)行微結(jié)構(gòu)的重構(gòu)。獲取單晶體的材料參數(shù)最直接的方法就是通過制備單晶試樣進(jìn)行各種材料試驗(yàn)來確定[2-3]?？紤]到這種方法既耗時耗力，成本又太高，因此大多數(shù)采用Materials Studio模擬軟件[4]，通過建立晶體原子結(jié)構(gòu)模型計算得出。晶體塑性本構(gòu)模型的塑性參數(shù)通常通過擬合試樣的應(yīng)力-應(yīng)變曲線來獲得。

通過晶體塑性有限元方法來模擬多晶體有限元模型的應(yīng)力-應(yīng)變曲線[5]，通過誤差函數(shù)確定與實(shí)驗(yàn)得到的拉伸曲線吻合最好的應(yīng)力-應(yīng)變曲線，最終確定了18CrNiMo7-6合金鋼晶體本構(gòu)模型的最優(yōu)塑性參數(shù)，為后續(xù)的有限元模擬提供材料參數(shù)。

1 力學(xué)性能試驗(yàn)與有限元模型

1.1 18CrNiMo7-6合金鋼的拉伸實(shí)驗(yàn)

研究的材料為18CrNiMo7-6合金鋼，為出廠狀態(tài)(鍛造處理)，材料成分如表1所示。為了確定18CrNiMo7-6合金鋼晶體本構(gòu)模型的塑性參數(shù)，首先需要對該材料進(jìn)行拉伸實(shí)驗(yàn)，獲得應(yīng)力-應(yīng)變曲線。拉伸試樣采用圓棒試樣，試樣規(guī)格采用國家標(biāo)準(zhǔn)GB 228—2002標(biāo)準(zhǔn)規(guī)定的尺寸，如圖1所示。拉伸實(shí)驗(yàn)在MTS(370.25)試驗(yàn)機(jī)上進(jìn)行，按照應(yīng)變加載的方式進(jìn)行加載，加載應(yīng)變率為0.001/s，拉伸過程中采用引伸計對應(yīng)變進(jìn)行測量，試樣的應(yīng)力-應(yīng)變曲線如圖2所示。

表1 18CrNiMo7-6合金鋼的化學(xué)成分的質(zhì)量分?jǐn)?shù)Table 1 The mass fraction of chemical composition of 18CrNiMo7-6 alloy %

圖1 拉伸試樣的尺寸圖Figure 1 Dimensions of the tensile specimen

圖2 18CrNiMo7-6合金鋼試樣的應(yīng)力-應(yīng)變曲線Figure 2 Stress-strain curve of 18CrNiMo7-6 alloy specimen

對應(yīng)力-應(yīng)變曲線進(jìn)行數(shù)據(jù)擬合等一系列處理之后，得到18CrNiMo7-6合金鋼原材料的材料參數(shù)，如表2所示。

表2 18CrNiMo7-6合金鋼的材料參數(shù)Table 2 Material parameters of 18CrNiMo7-6 alloy

1.2 多晶體有限元模型

多晶體的幾何模型采用Voronoi圖的原理生成，Voronoi圖原理是由俄國數(shù)學(xué)家Georgy Fedosievych Voronoy于1908年定義并命名的[6]，其是由一組連接兩鄰點(diǎn)直線的垂直平分線構(gòu)成的連續(xù)多邊形組成(如圖3所示)。其實(shí)質(zhì)是按照臨近原則對空間的一種剖分形式，由剖分后所形成的多邊形(體)集合構(gòu)成[7]。Voronoi圖原理被廣泛運(yùn)用于生物細(xì)胞結(jié)構(gòu)、城市規(guī)劃、地理學(xué)、氣象學(xué)、結(jié)晶學(xué)、圖像處理和微結(jié)構(gòu)模擬等復(fù)雜問題中[8]。

利用MATLAB自帶的Voronoi函數(shù)，在一定的范圍內(nèi)隨機(jī)撒種子點(diǎn)，生成多晶體二維幾何模型，并提取幾何模型的頂點(diǎn)信息，按一定的順序存放在記事本文件中。采用Python語言腳本在ABAQUS中調(diào)用記事本文件中的頂點(diǎn)數(shù)據(jù)，經(jīng)過處理后在Part模塊中生成多晶體的微結(jié)構(gòu)模型?？紤]計算效率問題，本文中所建立的晶粒尺寸模型要比實(shí)際的晶體尺寸略大，因此本文內(nèi)容建立了一個包含16個晶粒的多晶體微結(jié)構(gòu)模型，如圖4所示。

圖4 帶有16個晶粒的微結(jié)構(gòu)模型Figure 4 Microstructure model with 16 grains

1.3 邊界條件和網(wǎng)格劃分

建立的幾何模型長和寬為1 mm，由16個隨機(jī)大小和形狀的晶粒組成，每個晶粒的晶體取向隨機(jī)分布。網(wǎng)格劃分采用全局種子點(diǎn)控制，網(wǎng)格尺寸0.01 mm，網(wǎng)格單元數(shù)量為10 647個，單元類型為平面應(yīng)力四節(jié)點(diǎn)完全積分單元。模型左端固定，采用位移邊界條件，約束U1和U2方向(即U1=0,U2=0)。在模型右端外部的中間位置建立一個參考點(diǎn)，該參考點(diǎn)與模型右端耦合，通過參考點(diǎn)施加位移邊界條件進(jìn)行拉伸模擬，如圖5所示。在分析中，模型被嵌入到一個更大的主體中。為了便于與試驗(yàn)對比，拉伸位移設(shè)置為0.05 mm。提取參考點(diǎn)的支反力及位移與拉伸試驗(yàn)的前半段進(jìn)行擬合。

圖3 二維voronoi圖生成原理Figure 3 Formation principle of 2D voronoi diagram

圖5 有限元模型的邊界條件和網(wǎng)格劃分Figure 5 Boundary conditions and meshing of the finite element model

1.4 率相關(guān)晶體硬化模型

(1)

Hutchinson[10]用一個簡單的冪分布形式來描述多晶體的蠕變：

f(α)(x)=x|x|n-1，

(2)

式中：n是率敏感指數(shù)，當(dāng)n趨向于無窮時，這個冪分布形式接近于一個率無關(guān)材料。

應(yīng)變硬化通過增量關(guān)系強(qiáng)度g(α)的演化來描述：

(3)

式中：hαβ是滑移硬化模量，該式表示對所有激活的滑移系統(tǒng)進(jìn)行求和；hαα(no sum)稱為自硬化模量；hαβ(α≠β)稱為潛硬化模量。Perice等[11]和Asaro[12-13]用一個簡單的形式來描述自硬化模量：

(4)

式中:h0是初始硬化模量；τ0是初始屈服應(yīng)力，它等于當(dāng)前強(qiáng)度g(α)(0)時的初始值；τs是階段I的應(yīng)力或者說是大塑性流動開始時的突破應(yīng)力；γ是所有滑移系統(tǒng)上的泰勒累積剪切應(yīng)變。

(5)

潛硬化模量用下式來表示:

hαβ=qh(γ)，α≠β,

(6)

式中:q是一個常數(shù)。這些硬化模量的表達(dá)式忽略了晶體中的包辛格效應(yīng)。

筆者采用Huang[14]單晶體模型的UMAT，由于該種率相關(guān)硬化模型只需要確定晶體的3個塑性參數(shù)，即初始屈服應(yīng)力τ0、階段I的應(yīng)力τs、初始硬化模量h0，其他參數(shù)從文獻(xiàn)中獲得。

2 本構(gòu)模型參數(shù)的確定及討論

2.1 彈性參數(shù)

根據(jù)文獻(xiàn)[15]可知，少量的合金元素并不會顯著地改變晶粒的彈性剛度，所以用體心立方晶體α-Fe的彈性常數(shù)作為18CrNiMo7-6合金的單晶材料常數(shù)[16]，即C11=230 GPa，C12=135 GPa，C44=117 GPa，齊納(Zener)[17]各向異性系數(shù)

(7)

符合各向異性。

2.2 塑性參數(shù)

圖6 不同滑移系組數(shù)對應(yīng)的應(yīng)力-應(yīng)變曲線與試驗(yàn)曲線的對比圖Figure 6 Comparison of stress-strain curve corresponding to the number of different slip system groups and the test curve

從圖6可以看出，當(dāng)只有一個滑移系開動時，即主滑移系{110}<111>開動，多晶體的宏觀應(yīng)力-應(yīng)變曲線有明顯的拐點(diǎn)，即進(jìn)入屈服階段，與試驗(yàn)曲線比較吻合。當(dāng)有2組或者3組滑移系開動時，多晶體的宏觀應(yīng)力-應(yīng)變曲線只有彈性變形，與試驗(yàn)曲線相差較大。因此采用晶體塑性有限元方法的仿真均只采用一個主滑移系開動。下面對晶體的3個塑性參數(shù)進(jìn)行研究，評價這些參數(shù)對多晶體宏觀應(yīng)力-應(yīng)變曲線的敏感性影響。

固定晶體塑性中第Ⅰ階段的應(yīng)力τs和初始硬化模量h0兩個參數(shù)，只改變初始屈服應(yīng)力τ0。通過ABAQUS進(jìn)行有限元仿真，將結(jié)果與試驗(yàn)曲線進(jìn)行對比分析表明，參數(shù)初始屈服應(yīng)力τ0的改變對多晶體的屈服點(diǎn)產(chǎn)生了影響，多晶體的屈服應(yīng)力隨著τ0的增大而逐漸增大，如圖7所示。

圖7 初始屈服應(yīng)力τ0對宏觀應(yīng)力-應(yīng)變曲線的影響Figure 7 Influence of initial yield stress τ0 on macroscopic stress-strain curve

固定晶體塑性中的初始屈服應(yīng)力τ0和初始硬化模量h0兩個參數(shù)，只改變第Ⅰ階段的應(yīng)力τs這個參數(shù)，通過ABAQUS進(jìn)行有限元仿真，將結(jié)果與試驗(yàn)曲線進(jìn)行對比，見圖8。可見參數(shù)第Ⅰ階段的應(yīng)力τs的變化對屈服點(diǎn)幾乎沒有影響，但對屈服點(diǎn)以后的應(yīng)力-應(yīng)變曲線產(chǎn)生了影響，當(dāng)接近τ0時幾乎與x軸平行。隨著τs的增加后半段的斜率也逐漸增加，當(dāng)增大到一定程度后幾乎不再變化。參數(shù)τs不能等于τ0，否則它將在方程(4)中產(chǎn)生奇異性。

圖8 第Ⅰ階段應(yīng)力τs對宏觀應(yīng)力-應(yīng)變曲線的影響Figure 8 Influence of stage I stress τs on macroscopic stress-strain curve

固定晶體塑性中的初始屈服應(yīng)力τ0和第Ⅰ階段的應(yīng)力τs兩個參數(shù)，只改變初始硬化模量h0這個參數(shù)，通過ABAQUS進(jìn)行有限元仿真，將結(jié)果與試驗(yàn)曲線進(jìn)行對比分析，見圖9。可見參數(shù)初始硬化模量h0的變化對屈服點(diǎn)也幾乎沒有影響，但對屈服點(diǎn)附近的曲線斜率產(chǎn)生了較大的影響，斜率隨著h0的增大而逐漸增大。

圖9 初始硬化模量h0對宏觀應(yīng)力-應(yīng)變曲線的影響Figure 9 Influence of initial hardening modulus h0 on macroscopic stress-strain curve

通過考察晶體的3個塑性參數(shù)(初始屈服應(yīng)力τ0、第Ⅰ階段的應(yīng)力τs、初始硬化模量h0)對多晶體宏觀上應(yīng)力-應(yīng)變曲線的敏感性分析，初步確定了它們的變化趨勢和大致的變化范圍，據(jù)此初步估算有限元的計算次數(shù)，估計計算所需要的時間成本。τ0在180～240 MPa范圍內(nèi)變化，每次變化10 MPa；τ0=200 MPa時，τs在200～300 MPa范圍內(nèi)變化，每次變化10 MPa；h0在50～600 MPa范圍內(nèi)變化，每次變化50 MPa。由于第I階段的應(yīng)力τs大于初始屈服應(yīng)力τ0且不相等，而初始屈服應(yīng)力不受影響，所以τ0每增加一次計算，τs的計算次數(shù)便減少一次。為了提高計算的準(zhǔn)確度，τs初始值比τ0大1 MPa，之后τs取整數(shù)，則計算次數(shù)如表3所示。

表3 有限元模擬所需的計算次數(shù)Table 3 Number of calculations required for finite element simulation

參考 Moussa等[21]給出的代價函數(shù)來判斷參數(shù)擬合的優(yōu)劣程度。代價函數(shù)的值越小，說明模型擬合的程度越好，此時對應(yīng)的塑性參數(shù)即為所求的最優(yōu)結(jié)果。選擇在線性回歸中最常用的代價函數(shù)形式中的均方誤差模型，其形式如下：

(8)

圖10 840組擬合曲線與試驗(yàn)曲線的對比示意圖Figure 10 Comparison of the 840 sets of fitted curves and test curves(Schematic diagram)

對840組仿真結(jié)果中提取出的數(shù)據(jù)進(jìn)行Python編程處理，并與試驗(yàn)得到的應(yīng)力-應(yīng)變曲線進(jìn)行對比分析，通過代價函數(shù)比對選擇，選出了一條最優(yōu)的擬合曲線(如圖10所示，其中只給出了具有代表性的幾組擬合曲線)。輸出該曲線的參數(shù)組合，最終得到的最優(yōu)晶體塑性參數(shù)組合如表4所示。最優(yōu)的參數(shù)組合為τ0=200 MPa，τs=270 MPa，h0=500 MPa。得出的最優(yōu)擬合曲線與實(shí)驗(yàn)曲線如圖11所示，從中可以清楚地看出兩條曲線吻合較好，其相對誤差范圍在0～10.4%。

表4 有限元擬合的最優(yōu)材料參數(shù)Table 4 Optimal material parameters for finite element fitting

圖11 最優(yōu)擬合曲線與試驗(yàn)曲線Figure 11 Best fit curve and test curve

3 結(jié)論

筆者采用晶體塑性有限元的方法，通過建立有限元模型模擬了應(yīng)力-應(yīng)變曲線，并與實(shí)驗(yàn)得到的應(yīng)力-應(yīng)變曲線進(jìn)行了比較。利用建立的誤差函數(shù)進(jìn)行分析，有限元分析的結(jié)果和實(shí)驗(yàn)獲得的應(yīng)力-應(yīng)變曲線吻合較好。最終確定了18CrNiMo7-6合金鋼的3個主要單晶塑性參數(shù)為τ0=200 MPa，τs=270 MPa，h0=500 MPa。