動靜組合加載數(shù)值模擬砂巖損傷演化規(guī)律

王曉雨， 李騫， 王偉， 張明濤， 王奇智

(1.石家莊鐵道大學(xué)國防交通研究所， 石家莊 050043； 2.武漢地鐵集團(tuán)有限公司， 武漢 430070；3.石家莊鐵道大學(xué)土木工程學(xué)院， 石家莊 050043； 4.河北科技大學(xué)建筑工程學(xué)院， 石家莊 050043)

在自然界中，巖石類材料隨地質(zhì)演變活動形成了大量的、復(fù)雜的和隨機性的缺陷，這些缺陷使得巖石在動荷載作用下破壞機理和損傷規(guī)律更復(fù)雜[1]。在深部砂巖鈾資源地浸開采滲透性低這一難題背景下，根據(jù)動損傷模型研究砂巖的損傷演化規(guī)律，對中國深部資源有效開采具有指導(dǎo)性意義[2]。

為了研究巖石的損傷破壞特征，許多學(xué)者在理論方法的基礎(chǔ)之上進(jìn)行了大量的數(shù)值模擬和室內(nèi)試驗，建立了巖石動態(tài)宏觀和細(xì)觀損傷模型。Grady等[3]從巖石內(nèi)部微裂紋特征出發(fā)，分析出被激活的裂紋密度服從兩參數(shù)的Weibull分布，由此建立了拉伸應(yīng)變率效應(yīng)各向同性的動態(tài)損傷模型。朱晶晶等[4]對花崗巖進(jìn)行了霍普金森壓桿(split Hopkinson pressure bar，SHPB)循環(huán)沖擊試驗，結(jié)合統(tǒng)計損傷理論和微元體強度規(guī)律，建立基于Weibull分布的動態(tài)損傷模型。高文學(xué)等[5]通過研究巖石損傷耗散能，計算能量耗散密度為建立新的動損傷模型提供條件。

在動態(tài)損傷試驗和數(shù)值模擬方面，崔晨光[6]利用LS-DYNA 非線性動力分析軟件建立數(shù)值試驗?zāi)Ｐ?，模擬巖石在不同應(yīng)變率下的動態(tài)破壞過程。姚歡迎等[7]在單軸壓縮條件下進(jìn)行了聲發(fā)射試驗，并在巖石聲發(fā)射特征基礎(chǔ)上對頁巖單軸壓縮損傷演化規(guī)律進(jìn)行了分析，提出受載巖石損傷先減小后增大結(jié)論，進(jìn)一步擴展對受載巖石內(nèi)部損傷演化機理的認(rèn)識。張雨霏等[8]研究了含粗糙節(jié)理面水泥砂漿試塊的動態(tài)破壞特征，分析了沖擊荷載和試樣表面形貌對節(jié)理面的影響，并通過能量法表征損傷變量。Wang等[9]以SHPB試驗為基礎(chǔ)并結(jié)合顆粒流離散元PFC2D程序，進(jìn)行循環(huán)荷載作用下砂巖細(xì)觀動態(tài)損傷和破壞過程分析，揭示了巖石宏觀破壞機理。

汪鑫等[10]在凍融循環(huán)試驗基礎(chǔ)上，結(jié)合不同凍融次數(shù)下試樣物理性質(zhì)及細(xì)觀結(jié)構(gòu)的變化情況，采用ANSYS有限元軟件進(jìn)行試樣凍融循環(huán)過程中應(yīng)力應(yīng)變及巖石內(nèi)部溫度變化的特征模擬，進(jìn)一步揭示了英安巖凍融損傷劣化的力學(xué)機理。聞磊等[11]利用霍普金森桿對砂巖進(jìn)行預(yù)加靜載循環(huán)沖擊，選用Logistic函數(shù)的逆函數(shù)表征巖石循環(huán)沖擊損傷，進(jìn)一步分析砂巖動態(tài)損傷演化特征，對其破碎形式和加載強度進(jìn)行了研究。劉兵兵等[12]運用非線性顯式動力分析程序ANSYS/LS-DYNA對水下鉆孔爆破進(jìn)行數(shù)值模擬，從不同起爆方式和不同堵塞材料方面進(jìn)行分析，研究了水下鉆孔爆破的效果。

大多學(xué)者在理論分析和實驗研究及數(shù)值仿真方面做了大量的工作，采用不同方法建立了損傷模型，對巖石宏細(xì)觀破壞機理進(jìn)行了探究，但利用數(shù)值模擬分析巖石不同加載方式下?lián)p傷演化規(guī)律的較少。選用SHPB數(shù)值模擬方式進(jìn)行砂巖動態(tài)沖擊試驗，結(jié)合裂紋密度法建立的動損傷模型，研究其不同加載方式下的損傷演化規(guī)律，為正確認(rèn)識砂巖動態(tài)破壞機理和損傷演化規(guī)律提供理論基礎(chǔ)，促進(jìn)爆破增滲技術(shù)的進(jìn)一步發(fā)展。

1 數(shù)值實驗?zāi)Ｐ偷慕?/h2>

1.1 砂巖靜態(tài)試驗物理力學(xué)參數(shù)

砂巖的靜態(tài)物理和力學(xué)試驗在石家莊鐵道大學(xué)工程力學(xué)實驗室進(jìn)行，試件尺寸取Φ50 mm×100 mm。通過單軸抗壓強度試驗測得單軸抗壓強度、彈性模量等力學(xué)指標(biāo)[13]，基本數(shù)據(jù)參考文獻(xiàn)[13]中靜態(tài)試驗數(shù)據(jù)和文獻(xiàn)[6]中的公式計算得來，且砂巖靜載物理力學(xué)參數(shù)如表1所示。

表1 砂巖靜載物理力學(xué)參數(shù)Table 1 Physical and mechanical parameters of sandstone under static load

1.2 有限元模型的建立

數(shù)值模擬中以LS-DYNA3D軟件為基本軟件進(jìn)行建模，此軟件可進(jìn)行顯式和隱式兩種求解方法，其中顯式計算輸出的d3plot文件占存小，計算快捷，收斂效果比較好，隱式求解需要進(jìn)行方程式的平衡迭代，計算復(fù)雜。因此本文模擬以顯式分析為主，不考慮隱式求解。桿件和試件材料模型分別采用MAT_ELASTIC和HJC(Holmguist-Johnson-Cook)關(guān)鍵字進(jìn)行設(shè)置，利用ANSYS計算模塊進(jìn)行計算。HJC模型最初用于研究混凝土極限力學(xué)特性，由于其對描述高速沖擊情況良好被巖石沖擊研究引用至今[14]。受試驗等條件的限制，HJC模型中21個參數(shù)除基本力學(xué)參數(shù)外均參考文獻(xiàn)[14]中數(shù)據(jù)。SHPB桿件參數(shù)如表2所示，砂巖模型參數(shù)如表3所示。

常規(guī)狀態(tài)下的模型在桿底端整體模型環(huán)向為自由無約束邊界；一維加載透射桿底部為點約束邊界，桿件環(huán)向邊界約束轉(zhuǎn)動，巖石無約束，無壓力施加，入射頂端為面約束邊界，施加入射應(yīng)力波曲線；三維加載整體桿件底端與環(huán)向為點約束邊界，巖石和入射桿頂端為面約束邊界，施加對應(yīng)荷載曲線。

表2 SHPB桿件模型參數(shù)Table 2 Parameters of SHPB bar model

表3 砂巖HJC模型材料參數(shù)Table 3 Material parameters of sandstone HJC model

2 不同加載方式砂巖動態(tài)沖擊模擬

2.1 常規(guī)動態(tài)沖擊模擬

巖石SHPB數(shù)值模擬是在實驗基礎(chǔ)上進(jìn)行的，為保證結(jié)果的可靠性和有效性，需要滿足實驗的兩個基本假定條件：一維應(yīng)力波原理和均勻性假定，本文模型在滿足上述兩個條件下建立。

模型參照霍普金森壓桿試驗儀構(gòu)件的主要尺寸建立，其中砂巖試樣尺寸選用直徑50 mm，長度25 mm，桿件部分由300 mm長度子彈、3 000 mm長度入射桿和2 000 mm長度透射桿(半徑均為50 mm)三部分組成。在模擬中為保證應(yīng)力平衡，避免波形劇烈震蕩，需要設(shè)置波形整形器。目前，整形器大多采用紫銅片，基于試驗設(shè)置整形器尺寸為半徑3.6 mm來進(jìn)行建模，整形后的波形如圖1所示。

從圖1中可以看出曲線開始上升比較緩慢，波形比較平滑，無劇烈震蕩效果，這樣可使試件在沖擊過程中達(dá)到應(yīng)力平衡狀態(tài)。

有限元分析包含前后處理兩方面，網(wǎng)格劃分作為前處理的重要一步，通常它的密度和劃分形式對模擬數(shù)據(jù)結(jié)果影響頗大。本模型采用Solid164實體單元進(jìn)行建模，考慮模型計算量和精度問題，選用映射網(wǎng)格劃分方法進(jìn)行劃分，同時，為了更好地分析巖石損傷的效果，網(wǎng)格劃分采用不同劃分密度。子彈、入射桿和透射桿劃分較組略且網(wǎng)格劃分密度大致相同，其網(wǎng)格尺寸沿軸向為10 mm，劃分結(jié)果為：子彈單元數(shù)3 151 個和節(jié)點數(shù)3 597 個；入射桿單元數(shù)31 501 個和節(jié)點數(shù)34 917 個；透射桿單元數(shù)21 001個和節(jié)點數(shù)23 317 個。試件劃分較密，網(wǎng)格尺寸沿軸向為1.25 mm，劃分結(jié)果為單元數(shù)21 301個和節(jié)點數(shù)23 017 個，砂巖試件局部放大圖如圖2所示。

圖1 波形整形圖Fig.1 Waveform chart

圖2 模型局部放大圖Fig.2 Model local magnific

在LS-DYNA軟件中接觸設(shè)置包括單面、雙向、單向、固一4種類型，前3種接觸類型的算法均采用罰函數(shù)法，固一連接接觸算法不固定，一般有罰函數(shù)法、動約束法和分布參數(shù)法3種算法[15]。模型接觸設(shè)置主要采用罰函數(shù)法，為了應(yīng)力波更好傳播以接近實際效果這一目的，全部設(shè)置成雙向接觸類型。子彈與桿件接觸選用面面自動接觸；桿件和砂巖的接觸設(shè)置為侵蝕接觸，設(shè)置接觸主從面，砂巖面為主，兩種接觸其他參數(shù)默認(rèn)系統(tǒng)設(shè)置。同時設(shè)置單元失效模擬砂巖沖擊破壞，其最大失效主應(yīng)變設(shè)為0.006，進(jìn)而分析砂巖破壞后體積變化[16]。

2.2 軸壓加載下砂巖數(shù)值模擬

軸壓加載模擬對應(yīng)試驗的一維動靜組合加載，對透射桿后端進(jìn)行約束設(shè)置，入射桿前端施加預(yù)應(yīng)力，模擬單向軸壓作用下砂巖的損傷狀況。軸壓加載方法有動態(tài)松弛法和dyna文件法，LS-DYNA中為解決隱式問題引入動態(tài)松弛法，用于優(yōu)先把預(yù)應(yīng)力和變形考慮進(jìn)去，避免應(yīng)力變化過快和不收斂情況出現(xiàn)，與顯示方式求解法不符，且dyna文件法對于隱式和顯示都適用，它第一階段進(jìn)行應(yīng)力初始化設(shè)置，第二階段將dyna文件導(dǎo)入，其適用后續(xù)分析輸入文件做大的改動[17]。

SHPB試驗中只有軸壓無圍壓時，砂巖與桿件接觸的兩端出現(xiàn)明顯的端部效應(yīng)，破壞形狀呈錐形，端部基本無破壞，這是受界面摩擦的影響，在法向荷載和滑動速度較高時，動摩擦因數(shù)小于靜摩擦因數(shù)，且僅對有軸壓的試件考慮界面摩擦力[18]，數(shù)值模擬則可以避免界面摩擦的影響，沖擊中受力易平衡。為更好地與文獻(xiàn)[13]中實驗對應(yīng)，應(yīng)力初始化接觸設(shè)置采用最簡單的面面接觸，設(shè)置靜摩擦因數(shù)為0.18，動摩擦因數(shù)為0.05。

由文獻(xiàn)[6]中提到采用入射波形式加載可以使試樣在破壞前達(dá)到應(yīng)力均勻，曲線震蕩不明顯。在模擬時考慮以常規(guī)狀態(tài)模擬的入射應(yīng)力與軸壓值相擬合的曲線作為入射波來替代子彈撞擊進(jìn)行軸壓加載數(shù)值模擬。模擬方案以軸壓6、10、14、16 MPa，入射波速8、10、12、14、16 m/s進(jìn)行20次沖擊模擬，以研究不同入射波速和軸壓對砂巖損傷演化規(guī)律的影響。以軸壓10 MPa，入射波速8 m/s為例，入射波曲線如圖3所示，入射波為半正弦波。

圖3 入射應(yīng)力波Fig.3 Incident stress wave

2.3 三維加載下砂巖數(shù)值模擬

三維加載中軸壓和圍壓共同施加，制定6組組合模擬方案，軸壓和圍壓分別為6、3 MPa，6、6 MPa，6、9MPa，12、3 MPa，12、6 MPa，12、9 MPa。使用dyna文件法進(jìn)行模擬，對透射桿底端點約束，試件和入射桿頂端面約束，通過DEFINE_CURVE關(guān)鍵字定義第一步的圍壓和軸壓曲線及第二步的圍壓和入射波曲線，并利用關(guān)鍵字LOAD_SEGMENT_SET導(dǎo)入先前所設(shè)置的曲線。

三維模擬為達(dá)到與試驗加載高度相似，以軸向壓力與速度波擬合曲線作為入射波施加于入射桿頂端，圍壓只施加于試樣部分。總共進(jìn)行30次模擬，以考慮相同軸壓不同圍壓和相同圍壓不同軸壓時砂巖的損傷破壞規(guī)律。同時，研究沖擊速度、軸壓和圍壓三者組合下對巖石破壞的影響。

3 砂巖損傷演化規(guī)律分析

3.1 損傷模型

砂巖損傷演化規(guī)律研究以損傷變量作為基礎(chǔ)進(jìn)行，損傷變量是表征結(jié)構(gòu)內(nèi)部劣化程度的量化參數(shù)[19]。損傷變量的定義形式有多種，以往研究中常采用彈性模量法和能量耗散法，其公式為

(1)

(2)

式中：D為損傷變量；E0和E表示巖石損傷前后的動態(tài)彈性模量；en和esum表示巖石第n次試驗的消耗能和試驗全過程耗能總和。

以上兩種常規(guī)方法試驗和計算比較復(fù)雜，同時會使損傷數(shù)值偏大，不適用于數(shù)值模擬。根據(jù)損傷變量定義特點，巖石損傷表征分為宏觀和微觀兩方面，宏觀損傷分析多以聲波波速法為主，多用于實驗，在數(shù)值模擬研究中適用以巖石各向同性為前提，大量裂紋混合，不考慮單裂隙情況，基于裂紋密度構(gòu)建損傷變量[20]。由巖石的孔隙率表征巖石的損傷變量，定義為巖石的裂紋密度，公式為

(3)

式(3)中：Vn為砂巖沖擊模擬前體積；V為砂巖沖擊模擬后失效單元體積變化量。

3.2 損傷演化特性分析

從損傷模型中可以得出，損傷變量和沖擊前后的體積變化量相關(guān)，不同沖擊速度造成巖石不同程度損傷，損傷度用損傷變量來表征。裂紋密度法以體積變化為基礎(chǔ)進(jìn)行損傷表征，模擬中體積改變是通過失效單元關(guān)鍵字MAT_ADD_EROSION進(jìn)行控制，利用LS-DYNA軟件中Measure功能得到巖石微小體積改變量，將巖石體積改變量與其初始體積相比，得到最終的裂紋密度。無軸壓和圍壓作用下沖擊速度與裂紋密度關(guān)系如圖4所示。從圖4中可以得出，裂紋密度隨著沖擊速度的增加呈線性上升趨勢，表明砂巖的損傷破壞程度與沖擊波的大小密切相關(guān)。

在深部開采工程中，地應(yīng)力狀態(tài)研究越來越重要，地應(yīng)力通常是由地層土壓力和地質(zhì)運動產(chǎn)生的構(gòu)造應(yīng)力組成，對其簡化后只考慮土層應(yīng)力進(jìn)行一維軸向壓力數(shù)值模擬。巖石軸壓加載數(shù)值模擬過程中不同沖擊速度與損傷關(guān)系如圖5所示，以沖擊速度為12 m/s為例，不同軸壓與損傷的關(guān)系如圖6所示。從圖5中可以看出，砂巖的損傷和沖擊速度仍呈正相關(guān)，隨著沖擊速度的增大，砂巖損傷逐漸增大，相同速度軸壓值越大，曲線越靠下，裂紋密度越小。從圖6中看出，隨著軸壓的增加，裂紋密度呈下降趨勢。

模擬中以某一時刻的孔隙體積為基礎(chǔ)進(jìn)行裂紋密度分析，由于常規(guī)狀態(tài)與軸壓加載狀態(tài)所建模型PART組成不同，單元劃分不同，所取破壞點時刻出現(xiàn)差異，導(dǎo)致出現(xiàn)圖4比圖5數(shù)值偏小現(xiàn)象，而且軸壓加載模擬中考慮的入射應(yīng)力波直接加載處于更為理想狀態(tài)下進(jìn)行的，與試驗中子彈撞擊模式有所出入，但在同一加載方式下，軸壓的增長又使損傷度降低。說明一端約束，一端施加動靜組合入射波曲線在一定程度上達(dá)到了試驗效果，同時驗證了軸壓作用下巖石初始損傷速率增長緩慢，軸壓對巖石的初始破壞速率具有一定抑制作用。從微觀角度看，由于軸壓存在控制了巖石內(nèi)部微裂紋的擴展，才會減緩損傷速率。

圖4 常規(guī)沖擊模擬損傷演化Fig.4 Damage evolution of conventional impact simulation

圖5 不同軸壓下沖擊速度與裂紋密度關(guān)系Fig.5 Relationship between impact velocity and crack density under different axial compression

為滿足實際工程的需要，進(jìn)行三維加載條件下的模擬，從兩方面考慮。第一是相同圍壓不同軸壓下巖石沖擊速度與損傷的關(guān)系，取圍壓為6 MPa，不同軸壓沖擊速度與損傷關(guān)系如圖7所示，從圖7中可以看出，軸壓大小對巖石損傷影響頗大，說明軸壓越小，損傷越大，曲線越靠上，這和只有軸壓時的情況類似。第二是軸壓固定，圍壓不同時巖石沖擊速度與損傷關(guān)系，為了更好地分析取軸壓為6 MPa，不同圍壓沖擊速度與損傷關(guān)系如圖8所示，從圖8中可以得到，圍壓6 MPa比圍壓3 MPa損傷更為嚴(yán)重，而圍9 MPa比圍壓6 MPa損傷降低，圍壓值超過軸壓值時，損傷減小。因三維加載入射應(yīng)力波和圍壓曲線進(jìn)行分開加載，當(dāng)圍壓值小于軸壓值時，巖石所受的圍壓靜載未起到抑制作用，反而增加了巖石的環(huán)向應(yīng)變，導(dǎo)致?lián)p傷增大。

圖6 不同軸壓與裂紋密度關(guān)系Fig.6 Relationship between different axial compression and crack density

圖7 相同圍壓不同軸壓與裂紋密度關(guān)系Fig.7 Relationship between crack density and different axial compression under the same confining pressure

圖8 不同圍壓與裂紋密度關(guān)系Fig.8 Relationship between crack density and confining pressure

4 結(jié)論

采用LS-DYNA軟件進(jìn)行不同加載方式下砂巖沖擊實驗數(shù)值模擬，并結(jié)合損傷模型計算巖石損傷，通過巖石損傷演化規(guī)律進(jìn)行分析，得出以下結(jié)論。

(1)波形整形器設(shè)置可以有效解決沖擊實驗?zāi)M中波形圖震蕩問題；在軸壓和圍壓加載模擬第一步中采用最為簡單的面面自動接觸，可以保證模擬的接觸設(shè)置更有效，第二步設(shè)置為侵蝕接觸，方便分析巖石的體積改變，進(jìn)而研究砂巖的損傷。

(2)利用dyna文件法分兩段對進(jìn)行模擬，先進(jìn)行應(yīng)力初始化，再導(dǎo)入下一段模擬可以有效避免計算不收斂情況出現(xiàn)，只有軸壓存在時，軸壓固定砂巖損傷隨沖擊速度增加而上升，但初始增長速率緩慢，說明軸壓對巖石的裂紋擴展具有抑制作用。

(3)在常規(guī)狀態(tài)下砂巖損傷變量隨著沖擊速度增大而增大；圍壓和沖擊速度不變時，軸壓對巖石損傷的影響與單獨軸壓加載情況類似；軸壓和沖擊速度均不變時，圍壓增加損傷減小(圍壓值大于等于軸壓值)，說明圍壓可以有效控制巖石微裂紋的擴展。