蔣逸飛, 張俊, 何勇, 劉哲超, 趙萬華

(西安交通大學(xué)機械制造系統(tǒng)工程國家重點實驗室, 710049, 西安)

材料的裂紋擴展是一個典型的復(fù)雜非線性問題,對航空航天、機械和兵器工程等領(lǐng)域的發(fā)展有著重要的影響。利用數(shù)值計算方法研究材料的裂紋擴展是學(xué)者們近年來的研究熱點[1-2]。目前眾多基于連續(xù)介質(zhì)理論建立的數(shù)值模擬方法一般都需要采用預(yù)先設(shè)置斷裂面[3-4]、刪除單元[5-6]、重新劃分網(wǎng)格[7-8]或添加特殊的黏聚力單元[9-10]等復(fù)雜手段來仿真裂紋擴展,而離散元方法(DEM)不受連續(xù)介質(zhì)假設(shè)的約束,能夠簡單、有效地跟蹤材料破壞的全過程[11]。

剛體-彈簧模型是離散元方法中最常用的一種,根據(jù)單元間的作用形式,該模型分為接觸型和連結(jié)型兩種。接觸型模型的單元間的力學(xué)參數(shù)由人為設(shè)定,通過具有相互作用的單元集合體來表達(dá)材料的宏觀力學(xué)性能,這些接觸參數(shù)需要借助于實驗進(jìn)行反復(fù)校核,建模效率較低[12-13]。連結(jié)型離散元法的相關(guān)參數(shù)由材料的力學(xué)屬性直接推導(dǎo)得出,不需要經(jīng)過實驗校核[14-15]。Tavarez等基于彈性理論推導(dǎo)得出了連結(jié)型剛體-彈簧模型的彈簧剛度,對混凝土梁受沖擊侵徹的破壞過程進(jìn)行了研究[16]。Wang等定義的彈簧剛度遵循L-J對勢函數(shù)的形式,模擬了含裂紋的平板在動態(tài)拉伸下的斷裂過程[17]。以上研究雖然能夠定性地得到材料斷裂后的結(jié)果,但是由于單元排布型式對仿真結(jié)果的影響,不能準(zhǔn)確模擬裂紋的起裂和擴展過程從而揭示材料的斷裂機理。

Silling等提出基于近場動力學(xué)方法增加鄰近單元作用力數(shù)可以有效地模擬裂紋擴展過程[18]。Chen等基于非局部彈簧格構(gòu)方法建立了多層鄰居單元的作用力模型,對材料的裂紋擴展軌跡進(jìn)行了研究[19]。上述近場動力學(xué)方法和非局部彈簧格構(gòu)方法雖然與離散元方法不完全相同,但是它們計算單元間作用力的思路類似,對于研究剛體-彈簧離散元模型具有一定的借鑒意義。本文基于剛體-彈簧離散元模型,提出一種多層鄰近單元受力的離散元計算方法,推導(dǎo)了2層和3層鄰近單元受力時的彈簧剛度值和斷裂判據(jù);通過計算材料在沖擊作用下的裂紋動態(tài)擴展過程,分析了多層鄰近單元在研究材料斷裂問題時的優(yōu)勢。

1　材料的離散元描述

1.1　離散單元排布與多層鄰近單元的關(guān)系

(a)離散單元二維六角點陣排布　　(b)連接元件圖1　平面彈塑性材料的離散單元表達(dá)

定義εab和γab分別為單元a和b之間彈簧的正應(yīng)變和切應(yīng)變

(1)

式中:una、unb、usa和usb分別為單元a和單元b間的法向相對位移和切向相對位移。彈簧的應(yīng)變和笛卡爾坐標(biāo)系下的應(yīng)變分量有如下關(guān)系

(2)

式中:l和m分別為單元a、b中心連線與x軸正方向夾角的余弦和正弦值。

定義六角點陣排布下單元的前3層鄰近單元分布位置如圖2所示。令單元a的鄰近單元編號依次為1～18,以單元a中心為原點,x軸過第2、5、14和17號單元中心建立坐標(biāo)系。前3層鄰近單元與中心單元之間的平衡位置矢量用Rn(n=1,2,…,18)表示,i、j表示坐標(biāo)軸單位向量,其中第1層鄰近單元

(3)

第2層鄰近單元

(4)

第3層鄰近單元

(5)

圖2　多層相鄰單元分布示意圖

1.2　彈簧的剛度計算

當(dāng)離散單元的排布和鄰近關(guān)系確定以后,材料的力學(xué)行為由彈簧的受力和變形來表達(dá)。彈簧的剛度可以通過單位體積內(nèi)材料與離散單元的能量等效關(guān)系推導(dǎo)得出

Uc=Uk

(6)

式中:Uc為材料的變形能密度;Uk為離散單元的等效變形能密度。

根據(jù)彈性理論,平面應(yīng)力時Uc可由下式表達(dá)

(7)

式中:E為彈性模量;μ為泊松比;εx、εy和γxy分別為材料在笛卡爾坐標(biāo)系下的x、y方向正應(yīng)變和切應(yīng)變。

Uk由中心單元與鄰近單元之間的彈簧變形能組成。中心單元的受力由與周圍鄰近單元相連的彈簧提供,因此在計算單元受力前需要確定周圍鄰近單元的層數(shù)或個數(shù)。在六角點陣排布條件下,單元a有q個鄰近單元,當(dāng)只搜索第1層鄰近單元時,q=6,搜索前2層鄰近單元時,q=12,搜索前3層鄰近單元時,q=18。Uk可表示為

(8)

聯(lián)立式(1)、(2)、(6)～(8),可以得到考慮不同鄰近單元層數(shù)作用的法向和切向彈簧剛度kng和ksg(g=1,2,3)。當(dāng)考慮第1層鄰近單元作用時

(9)

當(dāng)考慮前2層鄰近單元作用時

(10)

當(dāng)考慮前3層鄰近單元作用時

(11)

1.3　彈簧的失效判據(jù)

選用最大拉應(yīng)力準(zhǔn)則作為材料的斷裂判據(jù),當(dāng)材料受力超過抗拉強度σult時,材料發(fā)生破壞。該斷裂判據(jù)在剛體-彈簧模型中用彈簧失效來描述。規(guī)定當(dāng)2個單元間法向和切向合力的絕對值達(dá)到臨界值fcr時彈簧失效,材料在該彈簧所在位置處發(fā)生破壞。fcr為彈簧法向和切向合力的臨界值

(12)

式中:fncr和fscr分別為法向和切向彈簧的臨界值,由材料的拉伸強度推導(dǎo)得出[16]

(13)

將彈簧剛度值代入式(13),可以得到考慮不同相鄰層數(shù)作用時彈簧的失效臨界值fncrg和fscrg(g=1,2,3)。當(dāng)考慮第1層鄰近單元作用時

(14)

當(dāng)考慮前2層鄰近單元作用時

(15)

當(dāng)考慮前3層鄰近單元作用時

(16)

2　裂紋動態(tài)擴展算例

在裂紋動態(tài)擴展的研究領(lǐng)域,Kalthoff和Winkler[20]的動態(tài)剪切實驗非常經(jīng)典,如圖3所示。本文利用該實驗得到的裂紋擴展路徑對多層鄰近單元離散元仿真結(jié)果進(jìn)行驗證。實驗布置如圖3a所示,受沖擊剪切的平板尺寸為200 mm×100 mm×9 mm,平板左端有兩條沿水平軸平行對稱分布的切槽,裂紋長度為50 mm。子彈寬度與兩切槽間距同為50 mm,子彈以16.5 m/s的速度沿x軸正方向沖擊兩切槽間的平板材料。平板材料為馬氏體不銹鋼18Ni1900,其彈性模量E=190 GPa,密度ρ=8 000 kg/m3,泊松比μ=0.3,抗拉強度σult=2 020 MPa,Rayleigh波速cR=2 799 m/s。實驗中的平板和沖擊加載形式均沿x軸線對稱,在沖擊過程中位于上下兩切槽端部萌生的裂紋及其擴展軌跡也相互對稱。實驗發(fā)現(xiàn),當(dāng)沖擊速度小于30 m/s時平板上半部分的裂紋將會在切槽端部萌生并沿與x軸正向夾角約70°的方向擴展。圖3b所示為實驗中沖擊速度為16.5 m/s時,平板上半部分的裂紋擴展路徑[20]。

(a)實驗布置示意圖　(b)平板上半部分裂紋擴展路徑圖3　Kalthoff和Winkler的動態(tài)剪切實驗

(a)　　　平板上半部分　(b)　　　六角點陣　　(c)　　　六角點陣的仿真模型　　　1型排布　　　2型排布圖4　動態(tài)剪切實驗的離散元建模

根據(jù)上述動態(tài)剪切實驗的對稱性特點,為了降低數(shù)值計算量提高計算速度,本文僅對平板上半部分進(jìn)行建模仿真,如圖4所示。模型中下邊界受對稱條件約束,下邊界y方向速度vy=0 m/s,與子彈沖擊接觸的材料受到v0=16.5 m/s的階躍速度擾動,模型其他邊界均為自由狀態(tài)。選取離散單元半徑r=0.5 mm,平板上半部分的離散元模型由11 558個單元組成。為了分析單元受力時的鄰近單元層數(shù)對單元排布方向的敏感性,本文選用了2種六角點陣的排布構(gòu)型,如圖4b、圖4c所示,其中單元2型排布是1型旋轉(zhuǎn)90°后的構(gòu)型,2種構(gòu)型下的其他離散元參數(shù)均相同。

3　仿真結(jié)果與分析

3.1　1層鄰近單元

僅考慮1層鄰近單元的受力作用,在六角點陣1型構(gòu)型下仿真得到平板的裂紋擴展軌跡,如圖5a～圖5c所示。裂紋的擴展速度極快,在平板受到?jīng)_擊作用后裂紋產(chǎn)生分叉,其中主裂紋沿x軸水平方向延伸至平板右端邊界處,另一個主要裂紋沿與x軸正方向呈60°左右的角度斜向上方擴展。由圖5d～圖5f可以看出,采用六角點陣2型構(gòu)型時,裂紋從切槽端部沿與x軸成大約45°的角度向材料右端擴展。當(dāng)子彈沖擊作用經(jīng)歷80 μs左右,有一條新的裂紋從平板右端萌生并沿x軸負(fù)方向傳播。根據(jù)上述仿真結(jié)果可知,只考慮1層鄰近單元的作用時,單元排布構(gòu)型會嚴(yán)重影響到裂紋的擴展路徑。

(a)　　　1型排布,　　(b)　　　1型排布,　　(c)　　　1型排布,40 μs　　　60 μs　　　80 μs

(d)　　　2型排布,　　(e)　　　2型排布,　　(f)　　　2型排布,40 μs　　　60 μs　　　80 μs圖5　考慮1層鄰近單元作用的裂紋擴展路徑

圖6給出了平板裂紋的擴展速度圖,從中可以看出,在1型構(gòu)型模型下仿真的裂紋從10 μs開始擴展到52 μs時結(jié)束,而在2型構(gòu)型模型下仿真的裂紋擴展起止時間分別為18 μs和44 μs,并且在裂紋擴展的整個時間區(qū)間內(nèi),2種構(gòu)型裂紋的擴展速度趨勢相差很大。仿真結(jié)果顯示,裂紋的最大擴展速度約為10 km/s,遠(yuǎn)大于材料的Rayleigh波速2 799 m/s,這與物理事實相違背。可見,僅考慮1層鄰近單元的方法在研究材料的裂紋擴展問題時并不適用。

圖6　考慮1層鄰近單元作用的裂紋擴展速度

3.2　2層鄰近單元

采用2層鄰近單元搜索方法進(jìn)行仿真計算,圖7a～圖7c為六角點陣1型排布構(gòu)型下仿真得到的裂紋擴展軌跡,裂紋的傳播方向大致與x軸成45°角。在2型構(gòu)型下計算的裂紋擴展軌跡如圖7d～圖7f所示,起初裂紋沿與x軸成45°的方向傳播,裂紋向右傳播10 mm后產(chǎn)生分叉,主裂紋沿與x軸約70°方向繼續(xù)傳播,次裂紋沿x軸正方向延伸10 mm后止裂。

(a)　　　1型排布,　　(b)　　　1型排布,　　(c)　　　1型排布,40 μs　　　60 μs　　　80 μs

(d)　　　2型排布,　　(e)　　　2型排布,　　(f)　　　2型排布,40 μs　　　60 μs　　　80 μs圖7　考慮2層鄰近單元作用的裂紋擴展路徑

如圖8所示,由計算得到的裂紋擴展速度可知,在1型排布構(gòu)型下裂紋的最大傳播速度為2 700 m/s,而2型排布構(gòu)型的最大傳播速度為3 400 m/s,大于Rayleigh波速,說明2型排布構(gòu)型下的數(shù)值計算結(jié)果失真?？傮w上當(dāng)數(shù)值計算中考慮2層鄰近單元時,1型和2型排布構(gòu)型下的計算結(jié)果仍然有差別,與1層鄰近單元搜索方法相比,基于2層鄰近單元搜索方法的仿真結(jié)果與實驗結(jié)果的符合程度雖然好很多,但是仍相差較大。

圖8　考慮2層鄰近單元作用的裂紋擴展速度

3.3　3層鄰近單元

采用3層鄰近單元搜索方法計算所得到的裂紋擴展軌跡如圖9所示,可見2種不同單元排布構(gòu)型的結(jié)果十分接近,裂紋均從切槽端部沿與x軸成70°的方向延伸,仿真結(jié)果與Kalthoff和Winkler的實驗結(jié)果十分接近。

(a)　　　1型排布,　　(b)　　　1型排布,　　(c)　　　1型排布,40 μs　　　60 μs　　　80 μs

(d)　　　2型排布,　　(e)　　　2型排布,　　(f)　　　2型排布,40 μs　　　60 μs　　　80 μs圖9　考慮3層鄰近單元作用的裂紋擴展路徑

圖10給出了裂紋擴展速度圖,可見在1型構(gòu)型下仿真的裂紋起止時間分別為24 μs和80 μs,2型構(gòu)型下仿真的裂紋起止時間為24 μs和78 μs,并且2種構(gòu)型下仿真得到的裂紋擴展速度趨勢也很接近,最大裂紋傳播速度約為2 100 m/s,該速度值小于Rayleigh波速,符合物理事實。Belytschko等使用擴展有限元方法(XFEM)對動態(tài)剪切實驗進(jìn)行了仿真[21],計算得到的裂紋擴展速度如圖10中實線所示,該結(jié)果與本文的仿真結(jié)果近似。以上結(jié)果說明,當(dāng)采用3層鄰近單元搜索算法進(jìn)行離散元計算時,單元排布構(gòu)型對裂紋擴展路徑幾乎沒有影響,利用該方法可以對材料裂紋擴展問題進(jìn)行有效的仿真研究。

圖10　考慮3層鄰近單元作用的裂紋擴展速度

3.4　鄰近單元層數(shù)對裂紋擴展模擬的影響分析

連結(jié)型剛體-彈簧離散元模型通過假設(shè)單元呈規(guī)則排布,對材料進(jìn)行離散化網(wǎng)格劃分。在計算過程中,如果只考慮1層鄰近單元的作用,此時中心單元僅與周圍6個單元間有相互作用力,即中心單元靠6個位置的彈簧共同限制其受力和運動。這個中心單元和6個彈簧建立的力學(xué)模型并不能在每個方向上表達(dá)同樣的力學(xué)屬性,本質(zhì)上是用一種非各向同性的力學(xué)模型近似模擬各向同性材料。當(dāng)分析材料斷裂前的力學(xué)問題時,這種模型可以得到足夠高精度的結(jié)果,但是在需要分析材料的裂紋擴展問題時,該模型的非各向同性特征就會凸顯出來,導(dǎo)致模擬失真。

當(dāng)考慮的鄰近單元的層數(shù)增多時,中心單元會受到更多不同方向的彈簧的共同作用,會使模型中的非各向同性力學(xué)性質(zhì)得到各向同性均化。在本文的連結(jié)型彈簧-剛度離散元模型中,取3層鄰近單元可以得到有效的裂紋擴展過程的仿真結(jié)果。

4　結(jié)　論

本文研究了連結(jié)型剛體-彈簧離散元方法中參與受力計算的鄰近單元層數(shù)對裂紋擴展的影響,推導(dǎo)了2層和3層鄰近單元作用下的彈簧剛度系數(shù)和斷裂判據(jù)。通過對Kalthoff和Winkler的動態(tài)剪切裂紋擴展實驗進(jìn)行數(shù)值仿真,發(fā)現(xiàn)增加鄰近單元的層數(shù)可以減少單元排布構(gòu)型對裂紋擴展的影響,當(dāng)采用3層鄰近單元搜索方法時,可以有效地模擬材料的裂紋擴展軌跡和裂紋擴展速度,因此可作為采用連結(jié)型離散元方法進(jìn)一步研究材料動態(tài)斷裂過程的一種手段。

參考文獻(xiàn):

[1]LI Shaofan, LIU W K, ROSAKIS A J, et al. Mesh-free Galerkin simulations of dynamic shear band propagation and failure mode transition [J]. International Journal of Solids and Structures, 2002, 39(5): 1213-1240.

[2]PSAKHIE S G, SHILKO E V, GRIGORIEV A S, et al. A mathematical model of particle-particle interaction for discrete element based modeling of deformation and fracture of heterogeneous elastic-plastic materials [J]. Engineering Fracture Mechanics, 2014, 130(1): 96-115.

[4]MA Wei, LI Xiangwang, DAI Lanhong, et al. Instability criterion of materials in combined stress states and its application to orthogonal cutting process [J]. International Journal of Plasticity, 2012, 30/31: 18-40.

[5]LJUSTINA G, LARSSON R, FAGERSTR?M M. A FE based machining simulation methodology accounting for cast iron microstructure [J]. Finite Elements in Analysis and Design, 2014, 80: 1-10.

[6]MIGUéLEZ M H, SOLDANI X, MOLINARI A. Analysis of adiabatic shear banding in orthogonal cutting of Ti alloy [J]. International Journal of Mechanical Sciences, 2013, 75(10): 212-222.

[7]MARUSICH T D, ORTIZ M. Modelling and simulation of high speed machining [J]. International Journal for Numerical Methods in Engineering, 1995, 38(21): 3675-3694.

[8]UMBRELLO D. Finite element simulation of conventional and high speed machining of Ti6Al4V alloy [J]. Journal of Materials Processing Technology, 2008, 196(1/2/3): 79-87.

[9]SALIH S, DAVEY K, ZOU Z. Rate-dependent elastic and elasto-plastic cohesive zone models for dynamic crack propagation [J]. International Journal of Solids and Structures, 2016, 90(1): 95-115.

[10] SONG J H, BELYTSCHKO T. Cracking node method for dynamic fracture with finite elements [J]. International Journal for Numerical Methods in Engineering, 2009(3): 360-385.

[11] 劉凱欣, 高凌天. 離散元法研究的評述 [J]. 力學(xué)進(jìn)展, 2003, 33(4): 483-490.

LIU Kaixin, GAO Lingtian. A review on the discrete element method [J]. Advances in Mechanics, 2003, 33(4): 483-490.

[12] LIN Q, FAKHIMI A, HAGGERTY M, et al. Initiation of tensile and mixed-mode fracture in sandstone [J]. International Journal of Rock Mechanics and Mining Sciences, 2009, 46(3): 489-497.

[13] SADD M H, ADHIKARI G, CARDOSO F. DEM simulation of wave propagation in granular materials [J]. Powder Technology, 2000, 109(1/2/3): 222-233.

[14] CHENG Ming, LIU Weifu, LIU Kaixin. New discrete element models for elastoplastic problems [J]. Acta Mechanica Sinica, 2009, 25(5): 629-637.

[15] WANG G, AL-OSTAZ A, CHENG A H, et al. A macroscopic-level hybrid lattice particle modeling of mode-I crack propagation in inelastic materials with varying ductility [J]. International Journal of Solids and Structures, 2009, 46(22/23): 4054-4063.

[16] TAVAREZ F A, PLESHA M E. Discrete element method for modelling solid and particulate materials [J]. International Journal for Numerical Methods in Engineering, 2007, 70(1): 379-404.

[17] WANG G, OSTOJA-STARZEWSKI M. Particle

modeling of dynamic fragmentation: ITheoretical considerations [J]. Computational Materials Science, 2005, 33(4): 429-442.

[18] SILLING S A, EPTON M, WECKNER O, et al. Peridynamic states and constitutive modeling [J]. Journal of Elasticity, 2007, 88(2): 151-184.

[19] CHEN Hailong, LIN Enqiang, JIAO Yang, et al. A generalized 2D non-local lattice spring model for fracture simulation [J]. Computational Mechanics, 2014, 54(6): 1541-1558.

[20] KALTHOFF J F, WINKLER S. Failure mode transition at high rates of shear loading [C]∥Proceedings of the International Conference on Impact Loading and Dynamic Behavior of Materials. Oberursel, Germany: IFAM, 1987: 185-195.

[21] BELYTSCHKO T, CHEN Hao, XU Jingxiao, et al. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment [J]. International Journal for Numerical Methods in Engineering, 2003, 58(12): 1873-1905.