卞貴學(xué)，陳躍良，張 勇，王安東，丁文勇

（海軍航空工程學(xué)院青島校區(qū)，山東 青島266041）

基于微觀結(jié)構(gòu)的高強鋁合金應(yīng)力集中系數(shù)研究

卞貴學(xué)，陳躍良，張 勇，王安東，丁文勇

（海軍航空工程學(xué)院青島校區(qū)，山東 青島266041）

為了更好地理解鋁合金材料的微觀力學(xué)性能，基于MATLAB編寫了Voronoi算法的微觀結(jié)構(gòu)模擬程序，并將程序?qū)階BAQUS有限元軟件建立鋁合金晶粒模型.推導(dǎo)出六結(jié)點內(nèi)聚線單元模型的界面單元剛度矩陣，利用內(nèi)聚力模型的內(nèi)聚力－位移關(guān)系描述鋁合金晶粒界面間的粘著力（法向力）和摩擦力（切向力），建立了微觀晶粒結(jié)構(gòu)的有限元模型.研究結(jié)果表明：單個夾雜粒子隨著彈性模量的增加應(yīng)力集中系數(shù)先減小再增加；相對于單個夾雜粒子，兩個夾雜粒子的應(yīng)力集中會增加，當(dāng)d/r接近2時應(yīng)力集中系數(shù)明顯增加，當(dāng)d/r值處在6左右時應(yīng)力集中系數(shù)基本恢復(fù)到單夾雜粒子時的大小.夾雜粒子的形狀、數(shù)量及分布狀態(tài)對結(jié)構(gòu)微觀應(yīng)力集中均有影響.

鋁合金；Voronoi；微觀結(jié)構(gòu)；有限元；應(yīng)力集中系數(shù)

從結(jié)構(gòu)材料的破壞中人們逐漸認(rèn)識到材料的破壞或斷裂現(xiàn)象源自材料的微觀尺度或細(xì)觀尺度，如圖1所示.任何影響金屬材料微細(xì)觀結(jié)構(gòu)變化的因素都會影響材料的宏觀力學(xué)性能.為了更好地理解飛機用鋁合金材料的力學(xué)性能，弄清損傷破壞機理，既需要在宏觀力學(xué)層面上進(jìn)行研究，又需要在細(xì)觀力學(xué)、甚至微觀力學(xué)層面上進(jìn)行研究［1］.

隨著人們對材料斷裂破壞本質(zhì)認(rèn)識的加深，斷裂力學(xué)適用的局限性也被了解：斷裂力學(xué)無法描述損傷萌生、演化直至出現(xiàn)宏觀裂紋的過程［1－2］.前人對鋁合金做微觀試驗分析了大量的研究工作，Liao等［3］研究了微觀結(jié)構(gòu)和微觀夾雜粒子分布規(guī)律.McDowell等［4］研究了疲勞裂紋形成的多晶微觀結(jié)構(gòu)的敏感性，為結(jié)構(gòu)在微觀方面的的抗疲勞設(shè)計和減少裂紋成核部位提供了依據(jù).Li等［5］針對2026鋁合金研究得出裂紋成核與與含鐵的第二相粒子直接相關(guān).王習(xí)術(shù)等［6］通過SEM原位觀測方法試驗研究了鑄造鎂鋁系列合金的微觀破壞機理.歸納得到的材料裂紋的萌生主要有以下種形式：夾雜物和基體界面開裂、夾雜物自身斷裂、滑移帶開裂和晶界開裂等［1，2－6］.楊衛(wèi)、黃克智等［7－8］在宏微觀斷裂力學(xué)理論和宏微觀跨尺度損傷機理進(jìn)行了大量研究.本文重點研究含有夾雜第二相粒子的2A12鋁合金微觀結(jié)構(gòu)有限元模型，并考察單個或多個夾雜物對微觀應(yīng)力集中的影響.

圖1 結(jié)構(gòu)的損傷源自微觀尺度Fig.1 Damage of structure origin from microstructure scale

1 鋁合金材料微觀模擬

Voronoi圖生成算法如圖2所示［9－10］.圖2為算法，圖3為編制MATLAB程序利用100個隨機點生成的Voronoi圖.

圖2 Voronoi圖生成算法Fig.2 Generation algorithm of voronoi diagram

圖3 Voronoi圖Fig.3 Voronoi diagram

根據(jù)統(tǒng)計鋁合金材料的晶粒尺寸，得到鋁合金晶粒度分布函數(shù)，結(jié)合Voronoi圖，建立鋁合金多晶體二維晶粒模型；通過Voronoi算法生成金屬材料多晶體顯微結(jié)構(gòu)與實際結(jié)構(gòu)的實際微觀結(jié)構(gòu)相一致［1］.晶粒的取向和相鄰晶粒取向之間的差異對短裂紋的萌生有很大影響，Voronoi圖只能用來模擬晶核的初始數(shù)目和位置，但對于晶粒的取向，則無法模擬.由于晶粒的取向具有很大的隨機性.因此，圖3微觀組織圖中，晶粒的取向是隨機產(chǎn)生并確定的［11］.

2 材料晶粒間內(nèi)聚力關(guān)系

內(nèi)聚力模型（Cohesive Zone Model，CZM）可以描述界面間的黏著力和摩擦滑移.CZM把晶粒間微裂紋的萌生和擴展以一種顯式形式表示，同時符合材料晶粒間力的作用規(guī)律.CZM既可反映裂尖前緣（Forward Region）阻礙裂紋起裂的過程，也可描述起裂后弱化區(qū)（Wake Region）阻礙裂紋長大的過程.尤其對彈粘塑性材料，應(yīng)力達(dá)到最大值后破壞并不是瞬間發(fā)生，利用CZM來反映連續(xù)漸進(jìn)的斷裂過程在物理上是合理的［12］.

采用雙線性內(nèi)聚力模型，雙線性內(nèi)聚力模型的內(nèi)聚力－位移關(guān)系的方程表達(dá)式為［12］

式中：Tn為法向應(yīng)力值；Tt為切向應(yīng)力值；σmax為法向應(yīng)力最大值；τmax為切向應(yīng)力最大值；σmax與τmax對應(yīng)的張開位移值分別為在達(dá)到最大值時應(yīng)力開始從最大值減小至零，對應(yīng)的位移值為

本文所選的是二維六節(jié)點內(nèi)聚單元由兩條三結(jié)點二次等參線段組成.當(dāng)受力作用，相鄰的單元變形后，界面間有相對位移，并產(chǎn)生內(nèi)聚力.考慮到法向和切向的不同，引進(jìn)內(nèi)聚單元界面的局部坐標(biāo)系（t，n）［13－14］，t為切向，n為法向，δt為切向位移張開量（簡稱位移量），δn為法向位移張開量，Tt為切向內(nèi)聚力，Tn為法線內(nèi)聚力.圖4所示為六結(jié)點內(nèi)聚線單元模型，節(jié)點1，節(jié)點2，節(jié)點3在內(nèi)聚力單元的下面，節(jié)點1′，節(jié)點2′和節(jié)點3′在內(nèi)聚力單元的上面.以上每個節(jié)點都具有兩個自由度，因此，在整體坐標(biāo)系（x，y）下，設(shè)節(jié)點位移向量為

節(jié)點1，節(jié)點2和節(jié)點3的坐標(biāo)列向量為

圖4 六節(jié)點內(nèi)聚線單元模型Fig.4 Six node isoparametric of line unite of the CZM

等參變換為

雅可比行列式為

式中：

界面單元上表面和下表面的連續(xù)位移量為

其中N為形函數(shù)，即

界面上表面和下表面的相對位移量為

其中，

因此，在有限元中（也稱為整體坐標(biāo)）的坐標(biāo)系下，有限元中結(jié)點位移與內(nèi)聚力界面的單元連續(xù)位移的關(guān)系（相對）為

為利用內(nèi)聚力關(guān)系模型，需要得出單元剛度矩陣，界面本構(gòu)關(guān)系在局部坐標(biāo)上，要把Δu從有限元中的坐標(biāo)系轉(zhuǎn)到局部坐標(biāo)系上.如圖4所示坐標(biāo)系（t，n），切向矢量為

法向矢量為

轉(zhuǎn)換張量為

相對位移向量δ＝（δt，δn），決定了內(nèi)聚力界面上的切向位移量和法向位移量.

根據(jù)內(nèi)聚力模型中內(nèi)聚力—位移關(guān)系T(δ)，可以求出剪切模量矩陣D，即

界面的節(jié)點力為

由虛功原理得到界面剛度矩陣

式中h為界面單元厚度.

3 有限元模型

3.1 模型的假設(shè)

在晶粒級尺度上建立模型首先要解決的是晶粒大小的確定.利用內(nèi)聚力模型和連續(xù)介質(zhì)力學(xué)的相關(guān)知識做以下假設(shè)［14］：

1）按照Voronoi方法生成2A12鋁合金晶粒，參考文獻(xiàn)［1］對2A12鋁合金微觀晶粒尺寸統(tǒng)計分析，本模擬取晶粒尺寸為53 μm［1］；

2）晶粒間的力學(xué)關(guān)系定義為上節(jié)所研究的內(nèi)聚區(qū)模型中的位移－驅(qū)動力之間的關(guān)系，將基于連續(xù)體力學(xué)理論的有限單元方法應(yīng)用到晶粒尺度的力學(xué)分析中，從而進(jìn)行有限元分析；

3）材料參數(shù)不隨研究尺度的改變而變化.

3.2 晶粒尺度模型的建立

在晶粒級尺度模型中晶粒的形狀、大小和方向，根據(jù)第1節(jié)描述的Voronoi方法，生成含有橢圓形夾雜粒子的 Voronoi多邊形，認(rèn)為每個Voronoi多邊形是一個連續(xù)體，然后，對模型進(jìn)行網(wǎng)格劃分，單元類型和材料特性等參數(shù)不變.通過MATLAB編制相關(guān)程序獲取Voronoi圖的幾何信息，用MATLAB編寫了Voronoi方法的微觀結(jié)構(gòu)模擬程序.運用Python腳本語言編寫inp文件，導(dǎo)入軟件 ABAQUS建立鋁合金晶粒模型；利用ABAQUS有限元模擬軟件的前處理模塊CAE建立模型的幾何信息部分，并結(jié)合MATLAB程序得到2A12鋁合金微觀模型的ABAQUS輸入文件，利用ABAQUS的有限元計算模塊STANDARD進(jìn)行模擬分析計算，最終的可視化結(jié)果由ABAQUS的后處理模塊分析［15－16］.

3.3 模型網(wǎng)格劃分及邊界條件

根據(jù)上述假設(shè)及幾何參數(shù)和邊界條件，模型單元類型為四邊形平面應(yīng)變單元，圖5（b）為劃分單元個數(shù)為2 304個的晶粒模型，可以看出，網(wǎng)格劃分相對晶界較為粗糙.為了滿足晶界呈線性，計算時候還需進(jìn)一步細(xì)化網(wǎng)格，但相應(yīng)的計算時間也會增加，本文計算的微觀結(jié)構(gòu)模型共劃分的單元個數(shù)為36 864個.模擬單向拉伸情況下的應(yīng)力分布時邊界條件為：AB邊固定約束，CB、DA保持平直，CD受均勻位移（拉力）作用.

圖5 晶粒尺度的有限元模型網(wǎng)格劃分Fig.5 Finite element model for micro grain scale

用隨機數(shù)的方法確定晶粒取向來保證2A12鋁合金是各相同性材料.在考察夾雜粒子對微觀應(yīng)力集中影響時，假設(shè)夾雜粒子的彈性模量為變量，基體彈性模量為常量.考察夾雜粒子與基體脫粘應(yīng)力時，假設(shè)夾雜粒子的彈性模量小于基體.其中，基體的材料參數(shù)：彈性模量E、泊松比υ.本文采用的雙線性內(nèi)聚力關(guān)系包括3個參數(shù)［14］：內(nèi)聚力在法向和切向為最大值時所對應(yīng)的法向位移和切向位移和或初始線性模量；內(nèi)聚力的法向應(yīng)力和切向應(yīng)力最大值σmax和τmax；臨界張開法向位移和臨界張開切向位移和或法向斷裂能臨界值和切向斷裂能臨界值，如表1和表2所示.

表1 計算模型材料參數(shù)Table1 Material parameter of computed model

表2 計算模型晶界參數(shù)Table 2 Grain boundaries parameter of computed model

4 計算結(jié)果與分析

4.1 單個夾雜物對微觀應(yīng)力集中的影響

由于要考慮夾雜粒子的不同彈性模量，此處基體材料按照表1參數(shù)選擇，粒子的彈性模型按照基體的倍數(shù)變化.在滿足模型拉伸不破壞的情況下，并結(jié)合實際鋁合金微觀結(jié)構(gòu)分析特點，大部分裂紋在Al7Cu2Fe和Mg2Si夾雜粒子處產(chǎn)生，且這些粒子通常近似圓形，第二相夾雜粒子尺寸大小一般為15～40 μm.其中，Mg2Si夾雜粒子比基體軟，通常為脫粘裂紋，其彈性模型小于基體.選擇遠(yuǎn)場拉應(yīng)力為100 MPa，粒子的平面形狀為圓形.考慮夾雜粒子所處位置，在Voronoi圖程序進(jìn)行算法優(yōu)化，使得生成圓形粒子附近較小的晶粒和相鄰的較大晶粒相合，保證粒子所處的位置在多晶粒晶界處，符合實際微觀組織粒子的多數(shù)位置，且粒子半徑取25 μm.應(yīng)力云圖與局部網(wǎng)格如圖6所示.

圖6 單個夾雜粒子應(yīng)力云圖與局部網(wǎng)格劃分Fig.6 Stress distribution and local detailed mesh for single particles

圖7為夾雜粒子的彈性模量和鋁合金的彈性模量比值與應(yīng)力集中系數(shù)之間的關(guān)系，其中，E粒子表示夾雜粒子的彈性模量，E鋁合金表示鋁合金基體的彈性模量.在小于基體彈性模量時，隨著夾雜粒子彈性模量的增加應(yīng)力集中系數(shù)減??；在大于基體彈性模量時，隨著夾雜粒子彈性模量的增加應(yīng)力集中系數(shù)減增加.應(yīng)力集中系數(shù)定義為

式中：S0max＝100 MPa；Smax為不同彈性模量的夾雜粒子時的最大正應(yīng)力.

圖7 夾雜粒子的不同彈性模量和應(yīng)力集中系數(shù)之間的關(guān)系Fig.7 Relationship of stress concentration factor and differ?ent elastic modular

4.2 兩個夾雜物對微觀應(yīng)力集中的影響

夾雜粒子的大小、數(shù)量及分布狀態(tài)（夾雜粒子的相對位置）影響結(jié)構(gòu)內(nèi)部應(yīng)力集中.本節(jié)主要是針對兩個夾雜粒子間的大小和相對位置進(jìn)行計算.在建模時夾雜粒子半徑取r1＝25 μm，與之相對應(yīng)的另外一個粒子的半徑尺寸分別取r2＝25 μm和r2＝20 μm，兩粒子中心的距離為d.粒子和基體材料參數(shù)選擇按表1進(jìn)行設(shè)置.計算時取粒子位置的兩種極端情況：兩粒子排列方向與加載方向垂直和排列方向與加載方向一致.其他兩粒子成一定角度的情況介于兩者之間.圖8為夾雜粒子排列方向與加載方向垂直時，兩粒子距離不同時的有限元應(yīng)力云圖.

計算可知，當(dāng)兩個夾雜粒子之間的距離較小時，最大應(yīng)力點出現(xiàn)在兩個夾雜粒子相近的位置上.當(dāng)僅存在一個半徑為25 μm的夾雜第二相粒子時，Kt＝2.065.通過對圖9和圖10的分析可知，當(dāng)兩個夾雜第二相粒子方向與加載方向相同時，相對于第二相粒子的情況，第一個粒子處的應(yīng)力集中能夠得到緩和，隨著d/r比值的增加，趨于平穩(wěn).當(dāng)兩個夾雜粒子排列方向與加載方向垂直時，相對于單個夾雜粒子的情況，兩個夾雜粒子的應(yīng)力集中會增加，特別當(dāng)d/r比值接近2時（因為不考慮夾雜粒子相互接觸聯(lián)合的情況，故所選的d/r比值均大于2），應(yīng)力集中系數(shù)明顯增加.隨著d/r比值的增加，這種加強作用趨于平穩(wěn)，當(dāng)d/r比值處在6左右時，應(yīng)力集中系數(shù)基本恢復(fù)到單夾雜粒子時的大小.

圖8 夾雜粒子排列方向與加載方向垂直時的有限元應(yīng)力云圖Fig.8 Finite element stress distribution for double particles when particles plumb the load：（a）distance of parti?cle center being 75 μm；（b）distance of particle center being 150 μm

圖9 夾雜粒子排列方向與加載方向垂直時的應(yīng)力集中系數(shù)Fig.9 Stress concentration factor when particle plumb the load

圖10 夾雜粒子排列方向與加載方向相同時的應(yīng)力集中系數(shù)Fig.10 Stress concentration factor when particle along the load

5 結(jié) 論

基于Voronoi圖方法較好地模擬了含有夾雜物（第二相粒子）2A12鋁合金的微觀組織結(jié)構(gòu).運用內(nèi)聚力模型中的內(nèi)聚應(yīng)力－位移關(guān)系描述晶粒之間、晶粒與夾雜粒子之間的力學(xué)關(guān)系，實現(xiàn)了2A12鋁合金的微觀結(jié)構(gòu)有限元分析.通過對微觀結(jié)構(gòu)的有限元模型分析得出夾雜粒子的形狀、數(shù)量及分布狀態(tài)對結(jié)構(gòu)微觀應(yīng)力集中的影響.

［1］卞貴學(xué)，陳躍良，張勇，等.2A12鋁合金疲勞裂紋成核與擴展機理［J］.材料研究學(xué)報，2012，26（5）：521－526.BIAN Guixue，CHEN Yueliang，ZHANG Yong，et al.Microstructural mechanism fatigue crack nucleation and propagation of 2A12 aluminum alloy［J］.Chinese Journal of Materials Research，2012：26（5）：521－526.

［2］劉大海，謝永鑫，黎俊初.工藝參數(shù)工藝參數(shù)對7075鋁合金板材時效成形性的影響［J］.材料科學(xué)與工藝，2015：23（3）：50－56.DOI：10.11951/j.issn.1005－0299.20150310.LIU dahai，XIE yongxin，LI junchu.Influence of process parameters creep age forming formability of 7075aluminum alloy sheet［J］.Materials Science and Technology，2015：23（3）：50－56.DOI：10.11951/j.issn.1005－0299.20150310.

［3］LIAO Min.Probabilistic modeling of fatigue related microstructural parameters in aluminum alloys［J］.Engineering Failure Mechanics，2009：76：668－680.DOI：10.1520/JAI101558.

［4］MCDOWELL D L，DUNNE F P E.Microstructure?sensitive computationalmodeling offatigue crack formation［J］.International Journal of Fatigue，2010，32（9）：1521－1542.

［5］LI J X，ZHAI T，GARRATT MD，et al.Four?point?bend fatigue ofAA2026 aluminum alloys ［J］.Metallurgical and Materials Transactions A，2005，36：2529－2549.DOI：10.1007/s11661－005－0126－z.

［6］王習(xí)術(shù)，湯彬，陶沙.鑄造鎂鋁合金的微觀破壞機理原位觀測技術(shù)與應(yīng)用［J］.機械工程材料，2006：30（2）：1－5.WANG Xishu，TANG Bin， TAO Sha.In?site observation technology and application of fatigue crack initiation and propagation for cast magnesium alloys［J］.Materials for Mechanical Engineer，2006：30（2）：1－5.

［7］楊衛(wèi).宏微觀斷裂力學(xué)［M］.北京：國防工業(yè)出版社，1995：1－13. YANG Wei.Macro and Micro Fracture Mechanics［M］.Beijing： National Defense Industry Press，1995：1－13.

［8］黃克智，王自強.材料的宏微觀力學(xué)與強韌化設(shè)計［M］.北京：清華大學(xué)出版社，2003：12－23.HUANG Kezhi，WANG Ziqiang.Micro?mechanics and strengthening and toughening design of materials［M］.Beijing：Tsinghua University Press，2003：12－23.

［9］傅廷亮，尹雪濤，張揚.Voronoi算法模型及其程序?qū)崿F(xiàn)［J］.計算機仿真，2006，23（10）：89－91.FU Tingliang，YIN Xuetao，ZHANG Yang.Voronoi algorithm model and the realization of Tts program［J］.Computer Simulation，2006，23（10）：89－91.

［10］司良英，鄧關(guān)宇，呂程.基于Voronoi圖的晶體塑性有限元多晶幾何建模［J］.材料與冶金學(xué)報，2009，8（3）：193－195.DOI：10.14186/j.cnki.1671－6620.2009.03.001.SI liangying，DENG Guanyu，LU Cheng.Polycrystal geometry modeling of crystal plasticity finite element method with voronoi diagram［J］.Journal of Materials and Metallurgy，2009，8（3）：193－195.DOI：10.14186/j.cnki.1671－6620.2009.03.001.

［11］張超.晶體塑性成形的數(shù)值模擬［D］.濟南：山東大學(xué)，2005：43－50.ZHANG Chao.The Numerical Simulation of Crystal Plastic Deformation［D］.Ji′nan：Shandong University，2005：43－45.

［12］張軍.界面應(yīng)力及內(nèi)聚力模型在界面力學(xué)的應(yīng)用［M］.鄭州：鄭州大學(xué)出版社，2011：1－16.ZHANG Jun.Interface stress and cohesion model in the interface mechanics［M］.Zhengzhou：Zhengzhou University Press，2011：1－16.

［13］PHILLIPS D R，IESULAURO E，GLAESSGEN E H.Multiscale modeling of metallic materials containing embedded particles［C］//45th AIAA/ASME/ASCE/AHS/ASC Structures， Structural Dynamics ＆Materials Confer.April2004.California ：Palm Springs，2004：19－22.

［14］BJERKE T W，LAMBROS J.Theoretical development and experimental validation of a theramally dissipative cohesive zone model for dynamic fracture of amorphous polymers［J］.Jouranl of the Mechanics and Physics of Solids，2003，51（6）：1147－1170.

［15］吳艷青，張克實.利用內(nèi)聚力模型（CZM）模擬彈粘塑性多晶體的裂紋擴展［J］.應(yīng)用數(shù)學(xué)和力學(xué)，2006，27（4）：454－454.WUYanqing，ZHANG Keshi.Crack propagation in polycrystalline elastic?viscoplastic materials using cohesive zone models［J］.Applied Mathematics and Mechanics，2006，27（4）：454－454.

［16］黃劉剛.內(nèi)聚力模型的分析及有限元子程序開發(fā)［D］.鄭州：鄭州大學(xué)，2010：28－39.HUANG Liugang.The analysis of cohesive zone model and user?defined subroutine development in finite element method［D］.Zhengzhou：Zhengzhou Uiversity：2010：28－29.

Study on stress concentration coefficient of high strength aluminum alloy microstructures

BIAN Guixue，CHEN Yueliang，ZHANG Yong，WANG Andong，DING Wenyong
（Qingdao Branch，Naval Aeronautical Engineering Institute，Qingdao 266041，China）

To better understand the micro?mechanical properties of aluminum alloys，the simulation program of Voronoi algorithm for micro structure was written based on MATLAB，and the program is imported into ABAQUS finite element software to build the grain model of aluminum alloy.The interfacial element stiffness matrix of the six?node cohesive element model is deduced.The adhesion force（normal force）and the frictional force（tangential force）between the grain boundaries of the aluminum alloy were described by the cohesion?displacement relationship of the cohesive zone model（CZM），and the finite element model of micro grain structure was established.The results show that the stress concentration coefficient of single inclusion particles decreases and then increases with the increase of elastic modulus.The stress concentration of the two inclusions increases with respect to the single inclusion particles，and the stress concentration coefficient increases significantly when thed/ris close to 2.When thed/ris approximately equal to 6，the stress concentration coefficient is restored to the value of the single particles.The shape，number and distribution of inclusion particles have an influence on the microstructure stress concentration.

aluminum alloy；Voronoi；microstructure；FEM；stress concentration factor

TG171

A

1005－0299（2017）05－0041－07

2016－10－04.< class="emphasis_bold">網(wǎng)絡(luò)出版時間：

時間：2017－07－12.

國家自然科學(xué)基金資助項目（51075394，51375490）.

卞貴學(xué)（1981—），男，講師；陳躍良（1962—），男，教授，博士生導(dǎo)師.

卞貴學(xué)，E?mail：bianguixue＠163.com.

10.11951/j.issn.1005－0299.20160331

（編輯 呂雪梅）