基于SPH方法的花崗巖在沖擊作用下裂紋擴(kuò)展數(shù)值模擬

摘 要：對(duì)不同沖擊速度下花崗巖靶板的裂紋擴(kuò)展分布情況進(jìn)行了數(shù)值模擬分析。采用非線性動(dòng)力學(xué)分析軟件Autodyn，Johnson-Homquist II（JH-2模型）損傷本構(gòu)模型，采用SPH算法，對(duì)花崗巖靶板進(jìn)行沖擊模擬，得到不同沖擊速度下裂紋擴(kuò)展分布的情況，研究不同沖擊速度對(duì)花崗巖靶板裂紋擴(kuò)展及起裂位置的影響。得到起裂位置、主裂紋垂直長(zhǎng)度與不同速度的關(guān)系曲線。結(jié)果表明：當(dāng)速度低于15m/s時(shí)，起裂位置距中心軸線較遠(yuǎn)且主裂紋擴(kuò)展長(zhǎng)度較短；當(dāng)速度大于20m/s時(shí)，起裂位置距中心軸線較近且主裂紋擴(kuò)展長(zhǎng)度增加，裂紋分布密度增大，當(dāng)達(dá)到一定速度后裂紋擴(kuò)展基本不變。

關(guān)鍵詞：花崗巖；沖擊；數(shù)值模擬；SPH算法；裂紋擴(kuò)展

1 概述

花崗巖由于其硬度高、耐磨性好、耐酸性好，被廣泛應(yīng)用于建筑裝修、防護(hù)工程、化工等方面，素有“巖石之王”之稱?；◢弾r等脆性材料在沖擊作用下的破壞過程一般是損傷、損傷積累、破碎、失效、主裂紋及裂紋群擴(kuò)展發(fā)生的過程，也就是巖石等脆性材料內(nèi)部裂紋的萌生、擴(kuò)展、斷裂的過程[1]。

許多學(xué)者對(duì)花崗巖等脆性材料裂紋擴(kuò)展方面做了很多研究，秦飛[2]等用邊界元法模擬了多裂紋的擴(kuò)展。胡柳青[3]等對(duì)在不同沖擊作用下裂紋的動(dòng)態(tài)響應(yīng)進(jìn)行了數(shù)值模擬分析。王海兵[4]等對(duì)不同速度下彈丸沖擊花崗巖標(biāo)靶產(chǎn)生的裂紋分布情況進(jìn)行了實(shí)驗(yàn)和數(shù)值分析，并歸納總結(jié)了一定速度范圍內(nèi)彈坑深度及最大裂紋長(zhǎng)度與彈丸速度之間的關(guān)系表達(dá)式。Motamedi[5]在網(wǎng)格不發(fā)生改變的情況下研究了復(fù)合材料中動(dòng)態(tài)裂紋完整性的擴(kuò)展問題。Daux[6]等建立了用于分析隨機(jī)分叉和交叉裂紋問題的擴(kuò)展有限元法。周小平[7]等對(duì)壓應(yīng)力狀態(tài)下的多裂紋擴(kuò)展問題采用擴(kuò)展有限元法進(jìn)行模擬分析。

傳統(tǒng)有限元法（Finite Element Method，F(xiàn)EM）以網(wǎng)格劃分為基礎(chǔ)，在處理高應(yīng)變率、大變形等問題時(shí)會(huì)導(dǎo)致網(wǎng)格扭曲，導(dǎo)致計(jì)算精度下降甚至計(jì)算難以繼續(xù)等問題。光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)法（Smooth Partcle Hydrodynamic Method，SPH），屬于純拉格朗日方法，無(wú)需進(jìn)行網(wǎng)格劃分，從而能較好的進(jìn)行高應(yīng)變率、大變形等不連續(xù)問題的處理。隨著計(jì)算機(jī)仿真技術(shù)的發(fā)展，應(yīng)用SPH法對(duì)塑性金屬材料切削過程的仿真得到了成功的應(yīng)用，但是應(yīng)用于脆性材料和準(zhǔn)脆性材料在沖擊侵徹過程中內(nèi)部裂紋擴(kuò)展機(jī)理研究的文獻(xiàn)較少。文章借助于非線性動(dòng)力學(xué)分析軟件Autodyn采用SPH算法基于Johnson-Homqu-istII（JH-2模型）損傷本構(gòu)模型，對(duì)錘面為平面的錘頭在不同速度下沖擊花崗巖靶板進(jìn)行裂紋擴(kuò)展數(shù)值模擬，并分析在沖擊作用下裂紋擴(kuò)展機(jī)理。

2 數(shù)值模擬

2.1 SPH方法

光滑離子流體動(dòng)力（SPH）法可以廣泛的對(duì)脆性材料和準(zhǔn)脆性材料等解體、破碎、固體的層裂、脆性斷裂等大變形問題進(jìn)行數(shù)值模擬仿真，不需要重構(gòu)網(wǎng)格，并且能保證計(jì)算的精度。其核心思想是將整個(gè)流場(chǎng)的物質(zhì)離散為一系列具有質(zhì)量、速度和能量的粒子，每個(gè)粒子具有單獨(dú)的速度、能量、質(zhì)量特征。是一種純拉格朗日的粒子方法，即無(wú)網(wǎng)格（meshfree）算法，且邏輯簡(jiǎn)單[8]。標(biāo)準(zhǔn)算法執(zhí)行步驟流程圖（如圖1所示）[9]。

SPH法是以插值理論為基礎(chǔ)的算法，借助核函數(shù)對(duì)一點(diǎn)上場(chǎng)變量的值給出核估計(jì)，把偏微分形式控制方程轉(zhuǎn)化為積分形式，其粒子運(yùn)動(dòng)信息的近似函數(shù)定義為：

在SPH算法[10，11]的計(jì)算上領(lǐng)域搜索是一個(gè)重要的步驟，在計(jì)算時(shí)，單個(gè)指點(diǎn)的影響范圍為以2h為半徑的球形區(qū)域，如圖3所示。其目的是在每個(gè)時(shí)間步列出該區(qū)域的所有質(zhì)點(diǎn)。當(dāng)質(zhì)點(diǎn)分離的時(shí)候，光滑長(zhǎng)度h的距離就會(huì)增加，相反，當(dāng)質(zhì)點(diǎn)匯聚的時(shí)候，光滑長(zhǎng)度h的距離就會(huì)減小。光滑長(zhǎng)度h的可變性使得球形領(lǐng)域內(nèi)的質(zhì)點(diǎn)數(shù)量保持不變，但需設(shè)置最大值與最小值，表示為：

式中：h0為初始光滑長(zhǎng)度，若HMIN，HMAX均為1時(shí)，則h為固定光滑長(zhǎng)度，不隨空間和時(shí)間變化。

2.2 Johnson-Homquist II（JH-2）損傷本構(gòu)模型

損傷本構(gòu)模型常常是用來(lái)描述巖石的基本力學(xué)行為，過去常用的本構(gòu)模型有流體模型、包絡(luò)線模型、Mohr-Coulomb模型、流體-彈塑性模型、RHT模型及Johson-Cook模型等，但是這些模型都不能完整的描述巖石破碎的變形和損傷行為[12]。

Johnson-Homquist II（JH-2）損傷本構(gòu)模型，假定脆性材料（無(wú)論是未損傷還是完全損傷的脆性材料），只考慮其損傷是塑性壓縮和剪切造成的。其狀態(tài)方程主要包括應(yīng)變率、靜水壓力以及拉伸應(yīng)力和多項(xiàng)式等。沖擊過程中材料性能最先體現(xiàn)為彈性應(yīng)變，直到應(yīng)力大于等于材料的屈服極限時(shí)，材料開始產(chǎn)生損傷積累并發(fā)生斷裂。隨著損傷的積累，材料最終完全破碎[13]。

JH-2材料模型（如圖4所示）是引入等效應(yīng)力來(lái)表示巖石等材料的一種動(dòng)態(tài)損傷本構(gòu)模型，其函數(shù)與冪函數(shù)形式的靜水壓力與應(yīng)變率和損傷因子相關(guān)[4，14]。

損傷因子D的函數(shù)表達(dá)形式與等效塑性應(yīng)變相關(guān)，其表達(dá)形式為：

式中：K1為體積模量，K2和K3為材料常數(shù)；？滋為體積應(yīng)變；？籽為當(dāng)前材料密度，？籽0為初始材料密度。

當(dāng)變形無(wú)法恢復(fù)時(shí)，巖石出現(xiàn)損傷，起裂開始。表示為：

式中：PC為壓碎壓力；PL為壓實(shí)壓力；？滋C為壓碎點(diǎn)體應(yīng)變；？滋L為壓實(shí)點(diǎn)體應(yīng)變。

當(dāng)起裂開始后，微裂紋逐漸擴(kuò)展并交叉貫通，巖石損壞破碎開始。表示為：

3 仿真模型建立與結(jié)果分析

3.1 仿真模型及材料參數(shù)

文章以平面錘頭，對(duì)花崗巖進(jìn)行10m/s、15m/s、20m/s、25m/s、30m/s不同速度下的垂直沖擊數(shù)值模擬仿真，錘頭采用Autodyn材料庫(kù)中TUNG.ALLOY材料，采用該材料的Shock狀態(tài)方程和Johnson-Cook材料模型?；◢弾r采用文獻(xiàn)[4]JH-2材料模型參數(shù)，JH模型可更好的模擬沖擊過程中裂紋的形狀，數(shù)量和長(zhǎng)度。材料模型和材料參數(shù)如表1和表2所示。

表1 花崗巖材料模型

3.2 二維模型建立

以錘面為平面的錘頭，采用二維算法，錘頭和花崗巖靶板均采用SPH算法，錘頭寬15mm，高10mm和花崗巖靶板高200mm，寬150mm，如圖6所示，計(jì)算模型如圖7所示。定義粒子尺寸為0.25，采用軸對(duì)稱的形式。

3.2.1 裂紋形成過程

為研究觀察裂紋演化過程，選取平頭錘分別以10m/s處的速度垂直沖擊花崗巖靶板進(jìn)行分析，觀察裂紋演化過程。如圖8所示為10m/s時(shí)不同時(shí)刻的損傷狀態(tài)及裂紋分布模擬情況。

沖擊速度以10m/s時(shí)不同時(shí)刻的裂紋擴(kuò)展情況如圖8所示，其中圖8（a）和（b）演示了沖擊初期花崗巖靶板受沖擊時(shí)的損傷狀態(tài)，并有細(xì)小裂紋產(chǎn)生，且起裂位置位于靶板頂部軸線中心左右兩側(cè)52.5mm處，主裂紋擴(kuò)展垂直高度為42mm，在30μs時(shí)，圖8（b）由于受壓縮和剪切力的作用使錘面下方出現(xiàn)“Λ”型裂紋。圖8（c），（d），（e）和（f）展示了在拉應(yīng)力下?lián)p傷積累后裂紋的擴(kuò)展，在40μs之后，上部裂紋擴(kuò)展停止，下部和中部損傷開始積累并逐漸形成裂紋并開始擴(kuò)展。在50μs后，上部裂紋擴(kuò)展開始繼續(xù)擴(kuò)展且起裂時(shí)出現(xiàn)變向，裂紋開始擴(kuò)展時(shí)逐漸趨向于原方向，且靶板的中下部分出裂紋開始擴(kuò)展。80μs后主裂紋幾乎不再擴(kuò)展，微裂紋裂擴(kuò)展速度減緩。圖8（g）為靶板最終狀態(tài)下的裂紋擴(kuò)展情況。

3.2.2 沖擊速度對(duì)裂紋的影響

為更好的觀察研究不同沖擊速度下裂紋擴(kuò)展情況，選取錘頭著靶速度分別為10m/s、15m/s、20m/s、25m/s、30m/s五種不同速度進(jìn)行數(shù)值模擬。得結(jié)果如圖9所示，裂紋擴(kuò)展與速度的關(guān)系如圖10，11所示。

由圖9可以明顯的觀察到，沖擊速度不同靶板裂紋也不同，沖擊速度越大，裂紋的密度隨之增大，且主裂紋的延伸長(zhǎng)度也隨之增大。速度越大，起裂的位置離中心軸線也越近，即所受壓縮和剪切作用越大。并且可以明顯觀察到速度越大，錘頭下方裂紋擴(kuò)展越密集，表現(xiàn)為破碎越徹底。由圖10，圖11可以看出，速度越大，起裂位置距中心軸線就越近，相應(yīng)的主裂紋擴(kuò)展垂直高度就越長(zhǎng)。速度大于25m/s時(shí)主裂紋垂直高度擴(kuò)展緩慢。速度大于15m/s時(shí)起裂位置向軸線中心靠攏緩慢。分析可知，存在一個(gè)最佳速度，使起裂位置更靠近中心軸線，主裂紋擴(kuò)展垂直長(zhǎng)達(dá)對(duì)大，且裂紋分布密度更均勻。

4 結(jié)束語(yǔ)

（1）針對(duì)花崗巖的特性，采用SPH法及JH-2損傷模型，在AUTODYN中使用二維軸對(duì)稱法建立花崗巖靶板沖擊破碎模型，形象的反映出裂紋擴(kuò)展分布情況及不同速度對(duì)裂紋擴(kuò)展的影響。

（2）起裂位置距離中心軸線的距離與速度有關(guān)，起裂位置決定了主裂紋的擴(kuò)展長(zhǎng)度，當(dāng)起裂位置距離中心軸線越近，主裂紋擴(kuò)展距離越長(zhǎng)且裂紋的擴(kuò)展分布密度越大。

（3）主裂紋擴(kuò)展垂直高度，并不是隨著沖擊速度的增加而增加，當(dāng)沖擊速度達(dá)到一定時(shí)主裂紋擴(kuò)展速度先增加后基本保持不變，此時(shí)微裂紋開始緩慢擴(kuò)展并逐漸與主裂紋及其他微裂紋相互貫通。

參考文獻(xiàn)

[1]席婧儀，陳忠輝，朱帝杰，等.巖石不等長(zhǎng)裂紋應(yīng)力強(qiáng)度因子及起裂規(guī)律研究[J].巖土工程學(xué)報(bào)，2015，04：727-733.

[2]秦飛，岑章志，杜慶華.多裂紋擴(kuò)展分析的邊界元方法[J].固體力學(xué)學(xué)報(bào)，2002，04：431-438.

[3]胡柳青，李夕兵，龔聲武.沖擊載荷作用下裂紋動(dòng)態(tài)響應(yīng)的數(shù)值模擬[J].爆炸與沖擊，2006，03：214-221.

[4]王海兵，壽列楓，張建鑫，等.彈丸撞擊下花崗巖靶破壞效應(yīng)實(shí)驗(yàn)與數(shù)值分析[J].巖石力學(xué)與工程學(xué)報(bào)，2014，02：366-375.

[5]Motamedi D，Mohammadi S.Dynamic Crack Propagation Analysis of Orthotropic Media by the Extended Finite Element Method[J].International Journal of Fracture，2010，161（1）：21-39.

[6]Christophe Daux， Nicolas Mos， John Dolbow， Natara jan Sukumar， Ted Belytschko. Arbitrarybranched and intersecting cracks with the extended finite element method[J]. Int.J.Numer.Meth.Engng.，2000，4812.

[7]周小平，楊海清，董捷.壓應(yīng)力狀態(tài)下多裂紋擴(kuò)展過程數(shù)值模擬[J].巖土工程學(xué)報(bào)，2010，02：192-197.

[8]肖毅華.有限元法與光滑粒子法的耦合算法研究[D].湖南大學(xué)，2012.

[9]石少卿，汪敏，孫波，等.AUTODYN工程動(dòng)力分析及應(yīng)用實(shí)例[M].北京：中國(guó)建筑工業(yè)出版社，2012：16-19.

[10]王澤鵬，胡仁喜，康士廷，等.ANSYS13.0/LS-DYNA非線性有限元分析實(shí)例指導(dǎo)教程[M].北京：機(jī)械工業(yè)出版社，2011：352-362.

[11]李運(yùn)祿，尹建平，王志軍，等.基于SPH算法對(duì)混凝土裂紋擴(kuò)展研究[J].兵器材料科學(xué)與工程，2015，05：19-22.

[12]王海兵，張海波，田宙，等.花崗巖介質(zhì)中爆炸回填堵塞的數(shù)值模擬[J].巖石力學(xué)與工程學(xué)報(bào)，2005，24（增1）：5013-5019.

[13]楊震琦，龐寶君，王立聞，等.JH-2模型及其在Al2O3陶瓷低速撞擊數(shù)值模擬中的應(yīng)用[J].爆炸與沖擊，2010，05：463-471.

[14]石祥超，陶祖文，孟英峰，等.致密砂巖Johnson-Holmquist損傷本構(gòu)模型參數(shù)求取及驗(yàn)證[J].巖石力學(xué)與工程學(xué)報(bào)，2015，S2：3750-

3758.

作者簡(jiǎn)介：王石安（1990-），男，漢族，河南省三門峽市人，碩士研究生，主要研究領(lǐng)域?yàn)闄C(jī)械工程。

*通訊作者：曾海峰（1974-），男，漢族，新疆昌吉人，副教授，碩士研究生導(dǎo)師，主要研究領(lǐng)域?yàn)闄C(jī)械制造及自動(dòng)化。

作者簡(jiǎn)介：張雨（1990-），男，漢族，湖北荊州人，碩士研究生，主要研究領(lǐng)域?yàn)闄C(jī)電一體化控制技術(shù)。

田佳（1992-），男，漢族，陜西蒲城人，碩士研究生，主要研究領(lǐng)域?yàn)樽詣?dòng)化控制。