基于SPH方法的花崗巖在沖擊作用下裂紋擴展數值模擬

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

關鍵詞：花崗巖；沖擊；數值模擬；SPH算法；裂紋擴展

1 概述

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

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

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

2 數值模擬

2.1 SPH方法

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

SPH法是以插值理論為基礎的算法，借助核函數對一點上場變量的值給出核估計，把偏微分形式控制方程轉化為積分形式，其粒子運動信息的近似函數定義為：

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

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

2.2 Johnson-Homquist II（JH-2）損傷本構模型

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

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

JH-2材料模型（如圖4所示）是引入等效應力來表示巖石等材料的一種動態損傷本構模型，其函數與冪函數形式的靜水壓力與應變率和損傷因子相關[4，14]。

損傷因子D的函數表達形式與等效塑性應變相關，其表達形式為：

式中：K1為體積模量，K2和K3為材料常數；？滋為體積應變；？籽為當前材料密度，？籽0為初始材料密度。

當變形無法恢復時，巖石出現損傷，起裂開始。表示為：

式中：PC為壓碎壓力；PL為壓實壓力；？滋C為壓碎點體應變；？滋L為壓實點體應變。

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

3 仿真模型建立與結果分析

3.1 仿真模型及材料參數

文章以平面錘頭，對花崗巖進行10m/s、15m/s、20m/s、25m/s、30m/s不同速度下的垂直沖擊數值模擬仿真，錘頭采用Autodyn材料庫中TUNG.ALLOY材料，采用該材料的Shock狀態方程和Johnson-Cook材料模型。花崗巖采用文獻[4]JH-2材料模型參數，JH模型可更好的模擬沖擊過程中裂紋的形狀，數量和長度。材料模型和材料參數如表1和表2所示。

表1 花崗巖材料模型

3.2 二維模型建立

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

3.2.1 裂紋形成過程

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

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

3.2.2 沖擊速度對裂紋的影響

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

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

4 結束語

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

（2）起裂位置距離中心軸線的距離與速度有關，起裂位置決定了主裂紋的擴展長度，當起裂位置距離中心軸線越近，主裂紋擴展距離越長且裂紋的擴展分布密度越大。

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

參考文獻

[1]席婧儀，陳忠輝，朱帝杰，等.巖石不等長裂紋應力強度因子及起裂規律研究[J].巖土工程學報，2015，04：727-733.

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

[3]胡柳青，李夕兵，龔聲武.沖擊載荷作用下裂紋動態響應的數值模擬[J].爆炸與沖擊，2006，03：214-221.

[4]王海兵，壽列楓，張建鑫，等.彈丸撞擊下花崗巖靶破壞效應實驗與數值分析[J].巖石力學與工程學報，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]周小平，楊海清，董捷.壓應力狀態下多裂紋擴展過程數值模擬[J].巖土工程學報，2010，02：192-197.

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

[9]石少卿，汪敏，孫波，等.AUTODYN工程動力分析及應用實例[M].北京：中國建筑工業出版社，2012：16-19.

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

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

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

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

[14]石祥超，陶祖文，孟英峰，等.致密砂巖Johnson-Holmquist損傷本構模型參數求取及驗證[J].巖石力學與工程學報，2015，S2：3750-

3758.

作者簡介：王石安（1990-），男，漢族，河南省三門峽市人，碩士研究生，主要研究領域為機械工程。

*通訊作者：曾海峰（1974-），男，漢族，新疆昌吉人，副教授，碩士研究生導師，主要研究領域為機械制造及自動化。

作者簡介：張雨（1990-），男，漢族，湖北荊州人，碩士研究生，主要研究領域為機電一體化控制技術。

田佳（1992-），男，漢族，陜西蒲城人，碩士研究生，主要研究領域為自動化控制。