巖石爆破損傷模型的比選與改進

胡英國 ，盧文波 ，陳 明 ，嚴 鵬 ，周創(chuàng)兵

（1. 武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072；2. 武漢大學 水工巖石力學教育部重點實驗室，武漢 430072）

1 引 言

巖體爆破是大型水利、采礦工程開挖必不可少的施工手段，爆破開挖將對圍巖造成一定程度的損傷，影響工程安全，因此，研究爆破損傷區(qū)特性，從而有效地控制圍巖的損傷范圍，具有重要的工程意義。巖體的爆破損傷是一個復雜的動態(tài)演化過程，為了研究探索這一過程，研究者們相繼建立了相關的巖體爆破損傷模型。

美國Sandia國家實驗室早在1966年就開始了巖石爆破損傷模型的研究工作[1]，主要方法是將巖石的動態(tài)斷裂作為一個連續(xù)的損傷累積過程來處理，其基本點是建立損傷變量與巖石內微裂紋密度的關系，并預測在爆炸載荷作用下巖石的損傷和破壞過程。Grady和Kipp[2]提出了巖石爆破各向同性損傷模型，即GK模型，該模型采用一個標量描述被拉應力激活的錢幣狀裂紋所引起的巖石剛度的劣化，同時假定這些裂紋數(shù)服從雙參數(shù)的 Weibull分布，他們采用該模型模擬爆炸載荷作用下油頁巖的動態(tài)斷裂和破碎，并根據能量平衡準則得到了與應變率有關的碎塊平均尺寸表達式。Taylor等[3]引進O 'connell、Budianshy[4]的關于有效體積模量和泊松比與裂紋密度的關系表達式以及 Grady 給出的碎塊尺寸表達式，建立了損傷變量與裂紋密度之間的關系式，并將損傷變量以率形式耦合到動態(tài)本構方程中，該TCK模型可以預報巖石在體積拉伸載荷下的動態(tài)響應。Kuszmaul[5]在以上兩模型的基礎上提出了KUS模型，該模型考慮了高密度微裂紋下的蔭屏效應，即微裂紋周圍產生應力釋放的材料能夠重疊，在裂紋的激活率中考慮了損傷引起的減少。Thorne等[6]在前人的基礎上考慮了激活裂紋數(shù)可能引起巖石體積的變化，并通過采用不同的損傷變量定義，考慮了模型在大裂紋密度條件下的適應性，建立了 Thorne模型。Yang、Liu等[7－8]對以上模型在裂紋密度的分布及損傷變量的定義方面進行了修正，認為只有在體積應變大于某一臨界體積應變后裂紋才能擴展，并考慮作用時間對裂紋密度的影響，在定義損傷變量時引入了斷裂概率的概念。隨著計算機硬件技術的發(fā)展，涌現(xiàn)出新的研究巖體爆破損傷過程的方法，唐春安等[9]不引進裂紋密度，采用累計拉應變和極限應變的比例關系來反映巖體的損傷程度，開發(fā)了用于研究巖體斷裂破壞過程的RFPA軟件。朱哲明等[10]將應力判據嵌入AUTODYN，并對應力更新算法進行修正，研究并揭示了巖體爆破損傷區(qū)形成的力學機制。王志亮等[11]基于 LSDYNA自定義接口將TCK模型導入LS-DYNA并計算，詳細研究了爆破拉伸損傷區(qū)的特性；馬國偉等[12]采用同樣的方法將Johnson- Holmquist模型嵌入 LS-DYNA，研究了爆破裂紋擴展的影響因素，并對爆破損傷控制提出建議。

采用損傷力學來研究問題時，其主要步驟為：首先定義一個合適的損傷變量，然后根據外荷載確定研究對象在荷載作用下的損傷演化方程和考慮損傷的本構關系。因此，反映巖體損傷程度的損傷變量的定義是基礎，正確、合理的損傷變量不僅能夠使要研究的問題簡單明了，而且計算結果也更加準確。已有的爆破損傷模型的損傷變量定義方式多樣，表1給出了5種應用較為廣泛的爆破損傷模型對應的損傷變量的表達方式及相關內變量定義的對比（嚴格意義上來說RFPA不是一種爆破損傷模型，但其中損傷變量的定義方法值得關注和借鑒，在這里為統(tǒng)一稱之為RFPA模型）。

從表1中可以看出，TCK模型與KUS模型的損傷變量表達式完全相同，只是在裂紋密度的表達中考慮了損傷的影響，而 THRONE模型的裂紋密度的表達式與TCK、KUS相似，但損傷變量的表達式差別明顯，YANG等的損傷模型的裂紋密度和損傷變量表達式與前三者均有差別，RFPA模型的損傷變量的定義方法與其他幾種模型區(qū)別最為明顯，沒有建立裂紋密度與損傷變量的關系，其表達式相對簡單。上述爆破損傷模型雖然在相關的研究中取得了較好的效果，但其適用性和精確性并沒有基于具體的工程實例進行系統(tǒng)地比較和驗證，本次研究中，將基于具體的工程實例，對上述5種爆破損傷模型進行對比計算和驗證。

表1 不同損傷模型相關變量的定義對比Table 1 Definitions of several variables for different damage models

2 巖體爆破損傷模型的比選

由彈性力學的基本理論可知，聲波在巖體中的傳播速度與巖體的物理力學性質密切相關，下面將給出通過聲波測試判據得到巖體損傷區(qū)的理論依據，巖體中縱波傳播速度可表示為

式中：VP為巖體聲波速度；ρ為巖體密度；K為體積模量；ν為材料的泊松比。定義巖體的聲波波速降低率η為：

采用η表示有效體積模量，若近似認為密度和泊松比不變，根據式（2）可得式（3）：

對比式（3）、（4）可以看出，聲波波速與爆破損傷模型描述爆破開挖的損傷效應是等價的，只是前者從宏觀聲波速度衰減來體現(xiàn)，而后者采用與微觀的裂紋密度密切相關的損傷變量來表征。

根據損傷變量式（3）、（4）可得爆破損傷模型的損傷變量D與聲波降低率η的對應關系為

根據《水工建筑物巖石基礎開挖工程技術規(guī)范》[13]規(guī)定，當波速降低率大于10%時，即判定巖體受到爆破損傷破壞，代入式（5）計算可得爆破損傷模型的損傷變量D的損傷閾值為0.19。

為方便進行爆破損傷模型的比選計算，宜選擇應力狀態(tài)較為簡單的工程實例以減少其他因素的干擾，本次研究選擇夏祥等[14]、李海波等[15]對廣東嶺澳核電站基巖爆破損傷區(qū)進行聲波檢測的實測結果為工程依據，文獻[14－15]提供了廣東嶺澳核電站的爆前、爆后的聲波監(jiān)測資料，并給出了圍巖體損傷范圍。圖1、2為聲波檢測測點布置示意圖以及圍巖體損傷區(qū)示意圖。

圖1 聲波檢測布置示意圖Fig.1 Arrangement of acoustic testing

圖2 損傷區(qū)實測示意圖Fig.2 Measured depth of damage zone in rock

根據相關的文獻資料，建立如圖3所示的模型，模型中圓柱直徑為 10 m，高為 8 m。裝藥半徑為25 mm，同時計算過程中巖體參數(shù)以及爆破損傷模型的相關參數(shù)如表2、3所示[3,14]。計算過程中，通過在孔壁上施加荷載曲線以模擬爆炸荷載的作用，爆炸荷載隨時間變化歷程參考文獻[16]的方法確定，其曲線形狀如圖4所示，邊界條件設為無反射邊界以減小應力波反射的影響，反映爆源周圍巖體的實際受力狀態(tài)。

圖3 模型示意圖Fig.3 Model sketch

表2 模型相關物理參數(shù)Table 2 Pysical parameters of model

表3 計算中損傷模型對應的損傷參數(shù)Table 3 Damage parameters of damage models

圖4 爆炸荷載隨時間變化曲線[16]Fig.4 Curves of blasting load versus time[16]

要實現(xiàn)以上爆破損傷模型的對比計算，需要借助 LS-DYNA的自定義接口，將其導入 LS-DYNA進行計算。實現(xiàn)LS-DYNA的二次開發(fā)的關鍵是用FORTRAN 語言編寫程序描述該模型建立的數(shù)學過程的子程序，在每個積分步，LS-DYNA主程序與該子程序進行數(shù)據傳遞。LS-DYNA中自定義方法在其用戶手冊[17－18]中有詳細介紹，這里不再贅述。

基于以上方法，實現(xiàn)上述5種爆破損傷模型，采用自定義接口將其成功導入LS-DYNA進行對比計算，圖5分別給出了5種爆破損傷模型對應的損傷區(qū)云圖。

圖5 不同損傷模型損傷區(qū)域云圖Fig.5 Nephogram of damage zone of different models

計算結果表明，從損傷區(qū)的分布來看，幾種爆破損傷模型的分布規(guī)律均與實測值大致吻合，圍巖中形成柱狀損傷區(qū)，損傷區(qū)的損傷程度隨著爆心距的增加而降低，炮孔附近巖體的損傷變量接近 1，由于炸藥爆炸瞬間產生的巨大的爆轟壓力，形成壓致粉碎區(qū)，巖體完全破壞，隨著遠離炮孔壁，沖擊波衰減為應力波，不能直接壓碎巖體，但巖體切向的拉應力超過其抗拉強度，圍巖中出現(xiàn)拉裂紋或原有裂紋進一步擴張，出現(xiàn)中等損傷程度的裂隙區(qū)，隨著爆炸荷載的進一步衰減，裂紋停止擴張，形成未損傷區(qū)。TCK模型、KUS模型和THRONE模型對應的損傷區(qū)外沿的整體形狀非常類似，只是不同位置對應的損傷程度以及損傷范圍不同；TCK模型的損傷范圍最大，而 THRONE的損傷范圍最小；YANG等建議的損傷模型損傷區(qū)的外沿與前3個模型偏差較大，由于RFPA模型損傷變量的定義方法與以上4種明顯不同，從計算結果可以看出，RFPA模型對應的損傷區(qū)，無論是損傷區(qū)外沿上還是損傷區(qū)內巖體對應的損傷程度與其他4種損傷模型區(qū)別相對明顯，由此可見，不同的損傷變量定義方法計算得到的損傷區(qū)的空間分布特征存在明顯區(qū)別。

為進一步比較不同損傷模型模擬爆破損傷區(qū)的準確性，表4給出了不同模型算例中損傷區(qū)的最大水平半徑與孔底損傷深度。

表4 不同損傷模型的損傷范圍對比Table 4 Damage area comparison between different models

表4數(shù)據表明，對于損傷最大水平半徑，最接近實測值的是KUS模型和RFPA模型，其中TCK模型、KUS模型、YANG模型偏大，最大誤差為0.4 m，對于孔底的損傷深度，除了PFRA外，其余幾種損傷模型的孔底損傷深度均偏小，其中誤差最大的是YANG等的損傷模型，其誤差達到了0.58 m，最接近實測資料的依然是KUS模型和RFPA模型。因此，計算結果表明，5種爆破損傷模型中，KUS模型和RFPA模型對應的損傷區(qū)范圍最接近實測結果。

3 拉壓損傷模型的建立與驗證

以上爆破損傷模型在考慮損傷效應時，均將拉應變與損傷變量建立直接的關系，這些損傷模型可以較好地反映巖體爆破產生的拉伸損傷，但不能有效地體現(xiàn)壓損傷的特點。在炸藥爆炸瞬間，巨大的爆轟壓力在孔壁周圍形成一定范圍的壓致粉碎區(qū)，這一點在上述模型中沒有考慮。因此，有必要建立能同時反應拉伸、壓縮損傷的損傷模型，使其具有更廣泛的應用價值。根據上文的計算結果，KUS模型計算的損傷范圍的精確性相對較高，本次計算將在對 KUS模型進行改進的基礎上，考慮引入壓損傷，建立拉壓損傷模型，下面將給出拉壓損傷模型的詳細數(shù)學推導過程。

根據 Taylor[3]等對損傷體積模量與裂紋密度的相關關系的計算表達式：

Grady等[2]認為，裂紋密度是裂紋影響區(qū)巖石體積與巖石總體積之比，激活的裂紋數(shù)服從體積拉伸應變的雙參數(shù)Weibull分布，即

式中：γ為隨機分布參數(shù)；N為單位體積裂紋數(shù)；a為微裂紋平均半徑。

爆炸應力波作用下的微裂紋平均半徑及裂紋密度可由式（8）、（9）決定：

Kuszmaul[5]考慮了高密度裂紋周圍應力釋放區(qū)的材料重疊，在裂紋的激活率中考慮了損傷引起的減少，將裂紋密度表達式修正為

式中：Dt為拉伸損傷變量。

考慮損傷效應的有效泊松比的表達式為

式中：θ為材料參數(shù)。

在體積拉伸狀態(tài)下，拉伸損傷變量 Dt可由式（12）表達：

在壓縮、剪切狀態(tài)下，巖體服從各向同性的彈塑性本構關系，其屈服條件為

對于壓縮損傷變量的考慮，引入 RDA模型對巖體中壓縮損傷的處理方法，基于 RDA模型的應變率效應耦合原則，可得出在體積壓縮下壓縮損傷變量的表達式[18]為

由于模型中同時包含拉伸損傷變量 Dt和壓縮損傷變量Dc，在表征巖體的損傷程度時，需要進行統(tǒng)一，在此定義損傷變量D，考慮損傷的最不利因素，D由下式決定：

得到損傷變量 D后，已有的爆破損傷模型（TCK、KUS等）的處理方法將所有宏觀彈性常數(shù)均按一定比例減小，如：

以上處理方法存在一定的缺陷，它忽略了巖體的彈性常數(shù)之間需要滿足的固定關系，通常情況下，巖體的宏觀彈性常數(shù)應滿足式（19）、（20）：

本次研究中，對宏觀彈性常數(shù)的確定方法作如下修正，損傷體積模量和有效泊松比沿用已有的爆破損傷模型（TCK、KUS等）的處理方法，其他彈性參數(shù)則根據上式中的固定關系獲得，改進后模型中的材料參數(shù)關系式如下所示：

其余宏觀彈性常數(shù)由式（23）、（24）確定：

最后，爆破荷載下巖石的記錄損傷效應的本構關系表達式可由以下增量型的虎克定律表示為

式中：σij為應力張量；εkk為體積應變張量；eij為偏應變張量；δij為Kronecker 符號。

通過以上推導，在對KUS模型進行修正的基礎上，考慮壓縮損傷建立了能夠同時反映拉伸、壓縮損傷的拉壓損傷破壞模型。將以上過程通過FORTRAN 編程實現(xiàn)，并成功導入LS-DYNA自定義接口，這樣就可以采用該模型進行計算并驗證其有效性。將以上改進的爆破拉壓損傷模型采用FORTRAN子程序實現(xiàn)，并成功導入LS-DYNA，對廣東嶺澳核電站基巖爆破損傷效應進行計算，得到相關的總體損傷范圍與局部損傷范圍如圖6所示。

圖6(a)表明，在KUS模型的基礎上進行修正的拉壓損傷模型最后形成的損傷區(qū)最大水平半徑為6.87 m，孔底損傷深度為2.10 m，相比原模型，無論是損傷區(qū)的分布還是范圍都與實測值更為接近。圖 6(b)給出了損傷區(qū)壓損圖，壓損區(qū)的范圍約為0.45 m，在總體損傷區(qū)中所占比例很小，巖體在爆破過程中，孔壁巖體受到強烈的沖擊波作用，超過巖體的抗壓強度，形成壓致?lián)p傷區(qū)，但由于爆炸荷載衰減很快，這個過程很短，因此，形成的壓碎區(qū)范圍很小。

圖6 拉壓損傷模型損傷區(qū)Fig.6 Damage zones of tension-compression damage model

該模型由KUS模型的基礎改進而來，因此，圖7給出了損傷區(qū)實測值、拉壓損傷模型以及KUS模型三者對應的損傷區(qū)對比示意圖，以便進一步驗證拉壓爆破損傷模型的有效性。

圖7 損傷區(qū)最大水平半徑對比Fig.7 Comparison of maximum horizontal radius of damage zone

相比KUS模型，拉壓損傷模型考慮了近區(qū)的壓縮損傷效應，計算結果表明，拉壓損傷模型與實測值的吻合性明顯要好，無論是最大損傷的水平半徑還是孔底的損傷深度均更加接近實測值，這表明本次研究中建立的拉壓損傷模型能夠更好地描述巖體中的爆破損傷效應。

4 結 論

（1）基于嶺澳核電站爆破損傷區(qū)檢測的實測資料，對5種爆破損傷模型進行比選，數(shù)值計算結果表明，KUS模型和 RFPA 模型與實測結果較為接近。

（2）根據比選結果，選取KUS模型，引進RDA的壓縮損傷的考慮方法，修正宏觀彈性常數(shù)的確定方法，建立新的拉壓損傷模型，并成功地導入LS-DYNA，進行計算和驗證。

（3）基于具體的工程實例對拉壓損傷模型的精確性進行驗證，結果表明，改進的拉壓損傷模型可更有效地描述爆破開挖擾動下巖體的爆破損傷效應。

但在實際工程中， 爆破近區(qū)巖體的復雜本構關系難以精確模擬，爆破近區(qū)巖體的復雜的破壞類型的精確描述尚需要進一步探索研究，本次計算基于各向同性的假定，沒有考慮實際工程中巖體存在的軟弱結構面的影響，對于實際工程巖體的各向異性特征沒有體現(xiàn)?？傮w而言，本次計算中采用的拉壓損傷模型對工程研究具有一定的參考價值和進步意義，但還需進一步修正和完善。

[1]楊小林, 王樹仁. 巖石爆破損傷模型及評述[J]. 工程爆破, 1999, 5(3): 71－75.YANG Xiao-lin, WANG Shu-ren. Rock blast damage models and comments[J]. Engineering Blasting, 1999,5(3): 71－75.

[2]GRADY D E, KIPP M E. Continuum modelling of explosive fracture in oil shale[J]. International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts, 1980, 17(3): 147－157.

[3]TAYLOR L M, CHEN E P, KUSZMAUL J S.Microcrack-induced damage accumulation in brittle rock under dynamic loading[J]. Computer Methods in Applied Mechanics and Engineering, 1986, 55(3): 301－320.

[4]BUDIANSKY B, O′CONNELL R J. Elastic moduli of a cracked solid[J]. International Journal of Solids and Structures, 1976, 12(2): 81－97.

[5]KUSZMAUL J S. A new constitution model for fragmentation of rock under dynamic loading[C]//Proceedings of the 2ndInternational Symposium on Rock Fragmentation by Blasting. Canada: Keystone, 1987: 412－423.

[6]THORNE B J, HOMMERT P J, BROWN B.Experimental and computational investigation of the fundamental mechanisms of cratering[C]//Proceedings of the 3rd International Symposium on Rock Fragmentation by Blasting. Brisbane, Australia: [s. n.], 1990: 117－122.

[7]YANG R, BRWDEN W F, KATSABANIS P D. A new constitutive model for blast damage[J]. International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts, 1996, 33(3): 245－254.

[8]LIU Li-qing, KATSABANIS P D. Development of a continuum damage model for blasting analysis[J].International Journal of Rock Mechanics and Mining Sciences, 1997, 34(2): 217－231.

[9]TANG C A. Numerical simulation of progressive rock failure and associated seismicity[J]. International Journal of Rock Mechanics and Mining Sciences, 1997,34(2): 249－261.

[10]ZHU Zhe-ming, XIE He-ping, MOHANTY B. Numerical investigation of blasting-induced damage in cylindrical rocks[J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(2): 111－121.

[11]WANG Zhi-liang, LI Yong-chi, WANG J G. Numerical simulation of tensile damage and blast crater in brittle rock due to underground explosion[J]. International Journal of Rock Mechanics and Mining Sciences, 2007,44(5): 730－738.

[12]MA G W, AN XM. Numerical simulation of blasting-induced rock fractures[J]. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(6):966－975.

[13]中國水利學會施工專業(yè)委員會爆破學組. SL47－94 水工建筑物巖石基礎開挖工程施工技術規(guī)范[S]. 北京:水利電力出版社, 1995.

[14]夏祥, 李俊如, 李海波, 等. 廣東嶺澳核電站爆破開挖巖體損傷特征研究[J]. 巖石力學與工程學報, 2007,26(12): 2510－2516.XIA Xiang, LI Jun-ru, LI Hai-bo, et al. Study of damage characteristics of rock mass under blasting load in Ling’ao nuclear power station, Guangdong province[J].Chinese Journal of Rock Mechanics and Engineering,2007, 26(12): 2510－2516.

[15]LI Hai-bo, XIA Xiang, LI Jian-chun. Rock damage control in bedrock blasting excavation for a nuclear power plant[J]. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(2): 210－218.

[16]盧文波, 楊建華, 陳明, 等. 深埋隧洞巖體開挖瞬態(tài)卸荷機制及等效數(shù)值模擬[J]. 巖石力學與工程學報, 2011,30(6): 1090－1096.LU Wen-bo, YANG Jian-hua, CHEN Ming, et al.Mechanism and equivalent numerical simulation of transient release of excavation load for deep tunnel[J].Chinese Journal of Rock Mechanics and Engineering,2011, 30(6): 1090－1096.

[17]LSTC. LS-DYNA theoretical manual[M]. California:Livermore Software Technology Corporation, 2003.

[18]LSTC. LS-DYNA keyword user’s manual[M]. California:Livermore Software Technology Corporation, 2003.

[19]FURLONG J R, ALME M L, DAVIS J F. Modeling the dynamic load/unload behavior of ceramics under impact loading RDA-TR-00.0-0001[R]. Arlington: R&D Associates, 1990.