基于FDM-DEM 耦合的沖擊損傷大理巖靜態斷裂力學特征研究*

張 濤，蔚立元，蘇海健，羅 寧，魏江波

（1. 中國礦業大學深部巖土力學與地下工程國家重點實驗室，江蘇 徐州 221116；2. 西安科技大學地質與環境學院，西安 陜西 710054）

隨著經濟建設的飛速發展，地下資源和空間的開發利用規模日益擴大，大量地下巖體被開挖掘進，而鉆爆法是目前最常用的施工方法。炸藥在爆炸破碎巖體的同時，也會導致地下工程預保留圍巖的動態擾動、損傷甚至破裂。循環爆破后，圍巖將受到不同程度的累積動態損傷，導致其在承受地應力重分布帶來的靜態荷載時斷裂韌度降低、承載能力下降，引發各種工程事故。因此，研究循環沖擊損傷后巖石的靜態斷裂特征具有較強的工程意義。

對脆性巖石在循環沖擊作用下的疲勞特性已有大量研究，但現有成果主要集中在循環沖擊下巖石損傷、破裂直至完全破壞的全過程。如Li 等采用分離式霍普金森壓桿（split Hopkinson pressure bar，SHPB）對花崗巖進行了循環沖擊加載實驗，結果表明，當沖擊氣壓在一定范圍內時，每次沖擊都會對巖石造成動態損傷但不至于完全破壞；林大能等研究了循環沖擊作用下圍壓、沖擊氣壓和沖擊次數對巖石動態損傷的影響；王彤等利用動靜組合加載裝置進行不同軸壓、不同沖擊氣壓下的循環沖擊實驗。而在受損巖石的靜態斷裂特性研究方面，損傷誘因也多以環境因素為主，如左建平等探究了不同溫度影響后花崗巖的細觀斷裂機制；賀晶晶等分析了凍融損傷對砂巖斷裂性能的劣化影響，并對試樣斷裂破壞面的形貌特征進行了掃描分析；楊健鋒等研究了不同程度水損傷作用對泥巖斷裂力學特性的影響。

低能量密度的循環沖擊能導致巖石的累積動態損傷乃至破裂，但是并不能使其完全破壞，此時巖石仍有一定的承載能力，但目前涉及循環沖擊損傷后巖石靜態承載能力的研究鮮有報道。雖然付安琪等利用SHPB 系統對中心直切槽半圓盤（notched semi-circular bend，NSCB）大理巖試樣進行了循環沖擊損傷處理，然后對其進行靜態斷裂實驗，分析了動態損傷對大理巖斷裂力學性能的劣化影響，但是巖石材料離散性較大，實驗結論適用范圍有限。此外，由于室內實驗條件的局限，該文獻中未對動態損傷累積過程及靜態斷裂力學行為的劣化機制進行深入研究。

隨著現代計算技術的不斷發展，多種數值分析方法被用來構建SHPB 模型并進行動態沖擊實驗，如有限元軟件ABAQUS、AUTODYN及離散元軟件PFC等。由于桿件由高強度鋼材加工而成，而巖石試樣是強度相對較低的脆性介質，在沖擊實驗過程中，桿件僅發生形變，而巖樣則將發生破裂甚至破碎。因此，利用連續介質模型模擬桿件系統、離散介質模擬巖石試樣是較為理想的數值建模方法。此外，僅采用三維離散元方法模擬SHPB 實驗時，只能通過提高桿件顆粒接觸黏結強度來近似滿足“應力均勻性假定”，并且桿件模型顆粒間的孔隙率對模擬結果的準確性也有一定影響。而連續-非連續耦合的建模方法則可以有效彌補這些缺陷，并能顯著提高計算效率，但目前這種耦合技術很少被應用于SHPB 實驗的模擬。

鑒于此，本文中利用Itasca 公司開發的FLAC（基于有限差分方法FDM）和PFC（基于離散元方法DEM）系列程序構建三維SHPB 耦合模型。利用該模型對模擬NSCB 試樣沿厚度方向進行不同次數的循環沖擊預損傷實驗，隨后對受損試樣進行靜態斷裂韌度模擬實驗，驗證相關室內實驗研究的結論，并結合微裂紋場、力鏈場等機理信息對動態損傷累積過程及斷裂力學性能劣化機制進行深入研究。

1 模 型

1.1 循環沖擊實驗

三維FDM-DEM 耦合SHPB 模擬系統如圖1 所示。其中桿件利用FLAC中的線彈性模型建立，尺寸及材料參數與室內實驗一致；NSCB 巖樣是基于PFC構建的離散顆粒黏結體，并選擇了能夠重現加載過程中的微裂縫萌生、聚結和宏觀裂縫形成的平行黏結模型（parallel bonding model，PBM），其尺寸構型與室內實驗一致。與大理巖三點彎曲室內實驗結果反復對比反饋，通過“試錯法”對模擬試樣細觀參數進行標定，最終獲得一組能夠反映脆性大理巖力學特性的細觀參數。

圖1 SHPB 模擬系統Fig. 1 Simulation of the SHPB system

模擬應力波如圖2 所示，可以看出，模擬波形的入射波和反射波之和近似等于透射波，并且在峰值后一定時間內仍保持相等。可見，本次模擬可實現試樣的動態受力平衡。與室內實驗一致，所有試樣分為6 組，利用SHPB 系統對各組試樣沿厚度方向分別進行了0～5 次沖擊速度恒定的等能量沖擊，以獲得損傷程度不同的6 組巖樣。

圖2 模擬應力波Fig. 2 Stress waves obtained from the numerical simulation

1.2 靜態三點彎曲斷裂實驗

預損傷實驗完成后，將子彈與桿件移除。在受損試樣的底部加設兩根支撐鋼棒，并在試樣上側以恒定速率施加荷載，如圖3 所示。圖4 中給出了完整試樣靜態三點彎曲實驗中荷載-位移曲線的室內實驗與數值模擬結果對比。可見模擬獲得的試樣峰值荷載及破壞位移均與實驗結果一致。

圖3 靜態三點彎曲模擬實驗Fig. 3 Simulation of the static three-point bending test

圖4 天然大理巖試樣實驗與模擬結果比較[18]Fig. 4 Comparison of the experimental and numerical results of the static three-point bending test of a natural marble sample[18]

2 循環沖擊結果分析

2.1 應力波信號

圖5 中給出了首次沖擊時桿件中的應力波傳播過程。當=0 μs 時，子彈撞擊桿件，產生的壓縮波由入射桿左端進入SHPB 系統。當=409 μs 時，壓縮波到達試樣的左端，隨即產生拉伸波（即反射波）向入射桿左端方向傳播。剩余應力波到達試樣右端時，再次發生反射-透射，透射波作為壓縮波向右端傳播。需要注意的是，在反射拉伸波之后，有一個壓縮波沿著入射桿向左端傳播（圖5 中黑色線圈內）。這是因為在相對較低的沖擊速度下，試件在加載期間（入射波上升段）未被完全破壞。在卸載期間（入射波下降段），試樣中儲存的部分應變能將沿著入射桿釋放，即發生回彈現象，從而產生壓縮波，這與Li 等的結論一致。

圖5 應力波在桿件中的傳播過程Fig. 5 Stress wave propagation in the bars

圖6 為循環沖擊過程中的應力波信號疊加曲線。從圖6 中可以看出，入射波幅值在每個沖擊周期中基本重合，說明模擬實現了等幅循環加載。隨著循環沖擊次數的增加，透射波的幅值越來越小，而反射波的幅值則越來越大。主要原因為，在循環沖擊作用下，巖石內部損傷累積，巖石孔隙率增加，而波阻抗降低。與應力傳播過程一致，在卸載段，入射桿上采集到的應力波信號有一段為負值，對應上述回彈現象。

圖6 循環沖擊加載過程中應力波信號變化Fig. 6 Variation of the stress wave signals during the cyclic impact loading

2.2 動態應力-應變曲線

圖7 為模擬試樣在循環沖擊過程中的動態應力-應變曲線，可劃分為3 個階段：彈性變形階段、裂紋擴展階段和應力卸載階段。在彈性變形階段，軸向應變幾乎隨動態應力的增加而線性增加。隨著應力的進一步增大，試樣進入微裂紋擴展階段，微裂紋萌生、擴展并相互作用，導致曲線呈現非線性行為。應力卸載過程對應于曲線的峰后段。

圖7 模擬動態應力-應變曲線Fig. 7 Dynamic stress-strain curves obtained from numerical simulations

圖8 為峰值應力的實驗值與模擬值變化情況。當沖擊次數由1 增加至5 時，峰值應力模擬值從75.1 MPa降低到64.9 MPa，降幅13.5%，遠小于實驗值的30.4%，這主要是因為天然試樣內部存在大量缺陷，這些缺陷會在動態沖擊過程中進一步發育，加劇試樣荷載傳遞能力的劣化。另一方面，隨著沖擊次數的增加，峰值應變增大，裂紋擴展，同時大量新裂紋萌生，卸載階段試樣恢復的應變能降低。

圖8 峰值應力實驗值與模擬值比較[18]Fig. 8 Comparison of the dynamic peak stress results obtained from experiments and numerical simulations[18]

2.3 動態損傷演變

PFC程序中，試樣顆粒之間的黏結斷裂即視為產生了微裂紋，微裂紋擴展嚴重區域常伴隨著顆粒脫落，形成碎塊。圖9 中給出了循環沖擊過程中模擬試樣裂紋場及碎塊場的演變過程。從圖9 中可以直觀了解到循環沖擊過程中試樣內部動態損傷演變過程：隨著沖擊次數的增加，試樣內部微裂紋及碎塊數目均有明顯上升。圖9（b）紅色圈內為撞擊掉落顆粒，與室內實驗結果（圖9（c））中觀察到的白斑及邊緣掉落顆粒一致，為巖石顆粒在動載作用下相互錯動導致。

圖9 動態損傷演變過程Fig. 9 Evolution process of the dynamic damage

為定量分析，圖10 中給出了微裂紋數量及碎塊數量隨加載次數的變化規律。可見試樣損傷演變過程可以劃分為急劇增長階段、緩慢發展階段和急劇增長階段：當= 1 時，首次沖擊造成的損傷較大，及分別由0 增長至3 558 和32；當= 2～4 時，單次沖擊造成的損傷較小，沖擊4 次后，及分別為6 402 和58；當= 5 時，單次沖擊造成的損傷較大，及分別增長至11 331 和198。

圖10 微裂紋及碎塊數量隨沖擊次數的變化規律Fig. 10 Variation of the microcrack number and fragment number with the cyclic impact number

模擬試樣顆粒之間具有黏結力，形成的力鏈在受到荷載作用時會改變，甚至斷裂。因此，力鏈演變能夠揭示受損試樣的動態損傷累積機理。圖11 中給出了每次沖擊后試樣接觸力場（剖面圖）的演變過程。

由圖11 可知，在動態荷載沖擊作用下，試樣內部力鏈分布均勻性變差，部分力鏈發生斷裂。當=5 時，試樣力鏈出現區域性缺失（黑色虛線圈內）。此外，試樣部分區域出現明顯的應力集中（紅色虛線圈內）現象，且均與室內實驗中試樣白斑及顆粒掉落位置一致。力鏈的斷裂及缺失使試樣在受到后續荷載作用時更易被破壞，這是受損試樣力學特征劣化的根本原因。

圖11 循環沖擊作用下試樣接觸力場演變過程Fig. 11 Evolution of the contact force field of a sample under cyclic impact loading

穿晶裂紋和沿晶裂紋的數量可以在一定程度上反映試樣結構的損傷程度。為了驗證上述微觀結構變化，圖12 中給出了受損試樣內部結構偏光顯微結果。可以看出，天然狀態下礦物顆粒完整，沒有明顯微裂紋。然而，經過動態沖擊后，顆粒之間黏結強度降低。當= 1 時，可以觀察到明顯的裂紋，說明巖石內部結構已被破壞；當= 2～3 時，裂紋數目增多，且周圍晶體顆粒的破碎程度明顯增加；當=4～5 時，裂紋已發生合并，發展為宏觀裂紋，周圍的晶體顆粒幾乎完全破碎，大理巖試樣內部結構的動態損傷已經非常嚴重。

圖12 受損大理巖的微觀圖像Fig. 12 Photomicrographs of the damaged marbles

3 三點彎曲結果分析

3.1 荷載-位移曲線

圖13 中給出了受損試樣靜態荷載-位移曲線的模擬結果及實驗結果。模擬曲線包含：（1）線彈性階段，荷載隨位移線性增長；（2）脆性破壞階段，試樣發生脆性斷裂，荷載瞬間跌落。隨著循環沖擊次數的增加，試樣內部損傷加劇，彈性段斜率隨著沖擊次數的增加而降低。模擬試樣曲線整體趨勢與室內實驗結果一致。

圖13 模擬靜態荷載-位移曲線對比[18]Fig. 13 Comparison of the static load-displacement curves obtained from experiments and numerical simulations[18]

3.2 宏觀力學特性

式中：α 為無量綱預制裂縫長度，α =/；/2是無量綱支撐間距。本次模擬中α= 0.2，/2= 0.7。受損試樣靜態斷裂韌度及破壞位移隨循環沖擊次數的變化規律如圖14 所示。

圖14 斷裂參數實驗值和模擬值比較[18]Fig. 14 Comparison of the fracture parameters obtained from experiments and numerical simulations[18]

當循環沖擊次數從0 增加至5 時，斷裂韌度由1.28 MPa·m減小至 1.07 MPa·m，破壞位移由0.21 mm 增長至0.27 mm。隨著的增加，斷裂韌度模擬值及實驗值均呈下降趨勢。破壞位移模擬值及實驗值變化過程基本一致。值得注意的是，模擬試樣參量的變化程度遠小于室內實驗中的天然試樣，這主要是由于模擬試樣內部不存在天然缺陷，而這些缺陷會在動態沖擊中進一步擴展，匯聚，進而對后續靜力學參數變化產生影響。

3.3 試樣破壞特征

圖15 中給出了模擬試樣斷裂破壞后的裂紋場及碎塊場演變。從圖5 中可以看到，受損試樣的破壞形式均為典型的三點彎曲斷裂破壞：從預制裂縫的頂端處起裂，擴展至加載點附近。

圖15 模擬試樣斷裂破壞形態Fig. 15 Fracture and failure patterns of the simulated sample

圖16 中給出了靜態加載過程中新增微裂紋及碎塊數量隨循環沖擊次數的變化規律。隨著沖擊次數的增加，受損試樣力鏈斷裂部分越來越多，導致靜態加載過程中更多的微裂紋萌生并擴展。由0 增加至5 時，新增微裂紋及碎塊數量分別由1 429、11 增長至1 887、34。當外部荷載作用于試樣、內部形成的集中應力大于試樣損傷閾值時，接觸斷裂，微裂紋產生。而當外部荷載、參與試樣變形的總接觸數量一定時，接觸斷裂現象越明顯，試樣越易發生宏觀破裂。結合圖13(a)所示，與天然試樣相比，當= 1～5 時，試樣剛度下降幅度分別為27.23%、36.05%、39.38%、55.53%、76.59%。因此，試樣的抗變形能力隨著沖擊次數的增加而逐漸減小。

圖16 新增微裂紋及碎塊數量隨沖擊次數的變化情況Fig. 16 Variation of the new microcrack and fragment number with the cyclic impact number

利用surfer 軟件對破壞后的試樣斷裂面進行后處理，得到斷裂面重建圖，如圖17 所示。為定量表征斷裂面粗糙度，采用斷裂面輪廓高度均值（）、斷裂面輪廓最大最小高度差值（）對斷裂面進行精確描述。

圖17 模擬試樣斷裂面形貌Fig. 17 Topography of the fracture surface of the numerical sample

斷裂面輪廓高度均值：

斷裂面輪廓最大最小高度差：

式中：為裂隙面測點數；z為裂紋面內點處的高度；、分別表示裂紋面內最高點、最低點處的高度。

由圖18 可知，斷裂面粗糙度隨沖擊次數的增加而增加。與天然狀態相比，循環沖擊5 次受損試樣斷裂面的、分別增長了32.06%、27.70%。在動態荷載作用下，微裂紋在試樣內部隨機產生，介質非均勻性、非連續性程度升高，致使試樣破裂面粗糙度上升，與室內實驗結論一致。

圖18 斷裂面粗糙度變化Fig. 18 Variation of the fracture surface roughness

4 結 論

基于有限差分法FDM 及離散元方法DEM 耦合思想，利用FLAC及PFC軟件構建了三維SHPB 模型。利用該模型對NSCB 試樣進行沿厚度方向的等能量循環沖擊，對受損試樣進行靜態三點彎曲實驗，探究了大理巖試樣在循環沖擊作用下動態損傷累積機理及受損試樣的靜態斷裂力學特征劣化機制，主要結論如下。

（1）壓縮作用下的應力波信號與實驗結果較吻合，從而驗證了FEM-DEM 耦合SHPB 系統用于動力加載的可行性。此外，微裂紋場、碎塊場、接觸力鏈結構、峰值強度、斷裂韌度等特征均與實驗結果吻合。

（2）在較低動態荷載作用下， 應力應變曲線會因試樣在卸載段釋放應變能而產生回彈。隨著循環沖擊次數的增加，試樣雖未發生整體破碎，但內部微裂隙、破碎顆粒、斷裂力鏈數量均增加，其動態峰值應力降低，動態損傷不斷累積。

（3）受損試樣的靜態斷裂韌度較天然試樣明顯劣化，破壞應變則呈上升趨勢。隨著沖擊次數的增加，試樣力鏈斷裂現象更加明顯。試樣介質非均勻性、非連續性程度在循環沖擊過程中不斷升高，致使試樣破裂面粗糙度上升、斷裂面粗糙度增加。