基于FDM-DEM 耦合的沖擊損傷大理巖靜態(tài)斷裂力學(xué)特征研究*

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

（1. 中國礦業(yè)大學(xué)深部巖土力學(xué)與地下工程國家重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221116；2. 西安科技大學(xué)地質(zhì)與環(huán)境學(xué)院，西安 陜西 710054）

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

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

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

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

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

1 模 型

1.1 循環(huán)沖擊實(shí)驗(yàn)

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

圖1 SHPB 模擬系統(tǒng)Fig. 1 Simulation of the SHPB system

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

圖2 模擬應(yīng)力波Fig. 2 Stress waves obtained from the numerical simulation

1.2 靜態(tài)三點(diǎn)彎曲斷裂實(shí)驗(yàn)

預(yù)損傷實(shí)驗(yàn)完成后，將子彈與桿件移除。在受損試樣的底部加設(shè)兩根支撐鋼棒，并在試樣上側(cè)以恒定速率施加荷載，如圖3 所示。圖4 中給出了完整試樣靜態(tài)三點(diǎn)彎曲實(shí)驗(yàn)中荷載-位移曲線的室內(nèi)實(shí)驗(yàn)與數(shù)值模擬結(jié)果對比?？梢娔M獲得的試樣峰值荷載及破壞位移均與實(shí)驗(yàn)結(jié)果一致。

圖3 靜態(tài)三點(diǎn)彎曲模擬實(shí)驗(yàn)Fig. 3 Simulation of the static three-point bending test

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

2 循環(huán)沖擊結(jié)果分析

2.1 應(yīng)力波信號

圖5 中給出了首次沖擊時(shí)桿件中的應(yīng)力波傳播過程。當(dāng)=0 μs 時(shí)，子彈撞擊桿件，產(chǎn)生的壓縮波由入射桿左端進(jìn)入SHPB 系統(tǒng)。當(dāng)=409 μs 時(shí)，壓縮波到達(dá)試樣的左端，隨即產(chǎn)生拉伸波（即反射波）向入射桿左端方向傳播。剩余應(yīng)力波到達(dá)試樣右端時(shí)，再次發(fā)生反射-透射，透射波作為壓縮波向右端傳播。需要注意的是，在反射拉伸波之后，有一個(gè)壓縮波沿著入射桿向左端傳播（圖5 中黑色線圈內(nèi)）。這是因?yàn)樵谙鄬^低的沖擊速度下，試件在加載期間（入射波上升段）未被完全破壞。在卸載期間（入射波下降段），試樣中儲存的部分應(yīng)變能將沿著入射桿釋放，即發(fā)生回彈現(xiàn)象，從而產(chǎn)生壓縮波，這與Li 等的結(jié)論一致。

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

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

圖6 循環(huán)沖擊加載過程中應(yīng)力波信號變化Fig. 6 Variation of the stress wave signals during the cyclic impact loading

2.2 動態(tài)應(yīng)力-應(yīng)變曲線

圖7 為模擬試樣在循環(huán)沖擊過程中的動態(tài)應(yīng)力-應(yīng)變曲線，可劃分為3 個(gè)階段：彈性變形階段、裂紋擴(kuò)展階段和應(yīng)力卸載階段。在彈性變形階段，軸向應(yīng)變幾乎隨動態(tài)應(yīng)力的增加而線性增加。隨著應(yīng)力的進(jìn)一步增大，試樣進(jìn)入微裂紋擴(kuò)展階段，微裂紋萌生、擴(kuò)展并相互作用，導(dǎo)致曲線呈現(xiàn)非線性行為。應(yīng)力卸載過程對應(yīng)于曲線的峰后段。

圖7 模擬動態(tài)應(yīng)力-應(yīng)變曲線Fig. 7 Dynamic stress-strain curves obtained from numerical simulations

圖8 為峰值應(yīng)力的實(shí)驗(yàn)值與模擬值變化情況。當(dāng)沖擊次數(shù)由1 增加至5 時(shí)，峰值應(yīng)力模擬值從75.1 MPa降低到64.9 MPa，降幅13.5%，遠(yuǎn)小于實(shí)驗(yàn)值的30.4%，這主要是因?yàn)樘烊辉嚇觾?nèi)部存在大量缺陷，這些缺陷會在動態(tài)沖擊過程中進(jìn)一步發(fā)育，加劇試樣荷載傳遞能力的劣化。另一方面，隨著沖擊次數(shù)的增加，峰值應(yīng)變增大，裂紋擴(kuò)展，同時(shí)大量新裂紋萌生，卸載階段試樣恢復(fù)的應(yīng)變能降低。

圖8 峰值應(yīng)力實(shí)驗(yàn)值與模擬值比較[18]Fig. 8 Comparison of the dynamic peak stress results obtained from experiments and numerical simulations[18]

2.3 動態(tài)損傷演變

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

圖9 動態(tài)損傷演變過程Fig. 9 Evolution process of the dynamic damage

為定量分析，圖10 中給出了微裂紋數(shù)量及碎塊數(shù)量隨加載次數(shù)的變化規(guī)律?？梢娫嚇訐p傷演變過程可以劃分為急劇增長階段、緩慢發(fā)展階段和急劇增長階段：當(dāng)= 1 時(shí)，首次沖擊造成的損傷較大，及分別由0 增長至3 558 和32；當(dāng)= 2～4 時(shí)，單次沖擊造成的損傷較小，沖擊4 次后，及分別為6 402 和58；當(dāng)= 5 時(shí)，單次沖擊造成的損傷較大，及分別增長至11 331 和198。

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

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

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

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

穿晶裂紋和沿晶裂紋的數(shù)量可以在一定程度上反映試樣結(jié)構(gòu)的損傷程度。為了驗(yàn)證上述微觀結(jié)構(gòu)變化，圖12 中給出了受損試樣內(nèi)部結(jié)構(gòu)偏光顯微結(jié)果?？梢钥闯觯烊粻顟B(tài)下礦物顆粒完整，沒有明顯微裂紋。然而，經(jīng)過動態(tài)沖擊后，顆粒之間黏結(jié)強(qiáng)度降低。當(dāng)= 1 時(shí)，可以觀察到明顯的裂紋，說明巖石內(nèi)部結(jié)構(gòu)已被破壞；當(dāng)= 2～3 時(shí)，裂紋數(shù)目增多，且周圍晶體顆粒的破碎程度明顯增加；當(dāng)=4～5 時(shí)，裂紋已發(fā)生合并，發(fā)展為宏觀裂紋，周圍的晶體顆粒幾乎完全破碎，大理巖試樣內(nèi)部結(jié)構(gòu)的動態(tài)損傷已經(jīng)非常嚴(yán)重。

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

3 三點(diǎn)彎曲結(jié)果分析

3.1 荷載-位移曲線

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

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

3.2 宏觀力學(xué)特性

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

圖14 斷裂參數(shù)實(shí)驗(yàn)值和模擬值比較[18]Fig. 14 Comparison of the fracture parameters obtained from experiments and numerical simulations[18]

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

3.3 試樣破壞特征

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

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

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

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

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

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

斷裂面輪廓高度均值：

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

式中：為裂隙面測點(diǎn)數(shù)；z為裂紋面內(nèi)點(diǎn)處的高度；、分別表示裂紋面內(nèi)最高點(diǎn)、最低點(diǎn)處的高度。

由圖18 可知，斷裂面粗糙度隨沖擊次數(shù)的增加而增加。與天然狀態(tài)相比，循環(huán)沖擊5 次受損試樣斷裂面的、分別增長了32.06%、27.70%。在動態(tài)荷載作用下，微裂紋在試樣內(nèi)部隨機(jī)產(chǎn)生，介質(zhì)非均勻性、非連續(xù)性程度升高，致使試樣破裂面粗糙度上升，與室內(nèi)實(shí)驗(yàn)結(jié)論一致。

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

4 結(jié) 論

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

（1）壓縮作用下的應(yīng)力波信號與實(shí)驗(yàn)結(jié)果較吻合，從而驗(yàn)證了FEM-DEM 耦合SHPB 系統(tǒng)用于動力加載的可行性。此外，微裂紋場、碎塊場、接觸力鏈結(jié)構(gòu)、峰值強(qiáng)度、斷裂韌度等特征均與實(shí)驗(yàn)結(jié)果吻合。

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

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