面向裂紋擴(kuò)展的連結(jié)型剛體-彈簧離散元模型

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

(西安交通大學(xué)機(jī)械制造系統(tǒng)工程國家重點(diǎn)實(shí)驗(yàn)室, 710049, 西安)

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

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

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

1　材料的離散元描述

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

(a)離散單元二維六角點(diǎn)陣排布　　(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軸正方向夾角的余弦和正弦值。

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

(3)

第2層鄰近單元

(4)

第3層鄰近單元

(5)

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

1.2　彈簧的剛度計(jì)算

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

Uc=Uk

(6)

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

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

(7)

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

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

(8)

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

(9)

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

(10)

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

(11)

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

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

(12)

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

(13)

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

(14)

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

(15)

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

(16)

2　裂紋動(dòng)態(tài)擴(kuò)展算例

在裂紋動(dòng)態(tài)擴(kuò)展的研究領(lǐng)域,Kalthoff和Winkler[20]的動(dòng)態(tài)剪切實(shí)驗(yàn)非常經(jīng)典,如圖3所示。本文利用該實(shí)驗(yàn)得到的裂紋擴(kuò)展路徑對多層鄰近單元離散元仿真結(jié)果進(jìn)行驗(yàn)證。實(shí)驗(yà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,抗拉強(qiáng)度σult=2 020 MPa,Rayleigh波速cR=2 799 m/s。實(shí)驗(yàn)中的平板和沖擊加載形式均沿x軸線對稱,在沖擊過程中位于上下兩切槽端部萌生的裂紋及其擴(kuò)展軌跡也相互對稱。實(shí)驗(yàn)發(fā)現(xiàn),當(dāng)沖擊速度小于30 m/s時(shí)平板上半部分的裂紋將會(huì)在切槽端部萌生并沿與x軸正向夾角約70°的方向擴(kuò)展。圖3b所示為實(shí)驗(yàn)中沖擊速度為16.5 m/s時(shí),平板上半部分的裂紋擴(kuò)展路徑[20]。

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

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

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

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

3.1　1層鄰近單元

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

(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層鄰近單元作用的裂紋擴(kuò)展路徑

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

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

3.2　2層鄰近單元

采用2層鄰近單元搜索方法進(jìn)行仿真計(jì)算,圖7a～圖7c為六角點(diǎn)陣1型排布構(gòu)型下仿真得到的裂紋擴(kuò)展軌跡,裂紋的傳播方向大致與x軸成45°角。在2型構(gòu)型下計(jì)算的裂紋擴(kuò)展軌跡如圖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層鄰近單元作用的裂紋擴(kuò)展路徑

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

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

3.3　3層鄰近單元

采用3層鄰近單元搜索方法計(jì)算所得到的裂紋擴(kuò)展軌跡如圖9所示,可見2種不同單元排布構(gòu)型的結(jié)果十分接近,裂紋均從切槽端部沿與x軸成70°的方向延伸,仿真結(jié)果與Kalthoff和Winkler的實(shí)驗(yàn)結(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層鄰近單元作用的裂紋擴(kuò)展路徑

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

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

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

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

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

4　結(jié)　論

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