基于HJC本構(gòu)模型的石灰?guī)r沖擊破壞形態(tài)數(shù)值模擬方法研究*

婁乾星，陶鐵軍，田興朝，謝財(cái)進(jìn)

(貴州大學(xué) 土木工程學(xué)院，貴陽 550025)

HJC模型是由T J Holmqusit等為解決高應(yīng)變率及高壓荷載下混凝土的大變形問題提出的一種本構(gòu)模型，目前在巖石材料的動(dòng)態(tài)沖擊破壞過程中應(yīng)用廣泛。由于巖石種類不同和地質(zhì)條件的差異性，該類材料的物理力學(xué)性質(zhì)較難統(tǒng)一。然而，其材料本構(gòu)模型的參數(shù)直接影響數(shù)值模擬結(jié)果的可靠性，因此對石灰?guī)r材料HJC模型參數(shù)的確定和驗(yàn)證是必要的。

目前，國內(nèi)外學(xué)者針對HJC模型參數(shù)方面的研究較多，任亮等通過數(shù)值模擬和實(shí)驗(yàn)分析了HJC模型強(qiáng)度參數(shù)A與參考應(yīng)變率之間的相關(guān)性[1]，并對參數(shù)進(jìn)行了優(yōu)化。米振國等通過ANSYS-LS-DYNA優(yōu)先軟軟件研究了適用于130 mm穿甲侵徹混凝土靶的HJC本構(gòu)模型參數(shù)，并以加速度為變量進(jìn)行了研究[2]。張嘉凡等把HJC本構(gòu)模型運(yùn)用到煤巖中，通過數(shù)值模擬的方式確定液態(tài)CO2爆破效果[3]。劉錦等采用ANSYS/LS-DYNA有限元軟件中的HJC本構(gòu)模型結(jié)合實(shí)驗(yàn)的方式重現(xiàn)了煤巖的動(dòng)態(tài)壓縮的破壞過程[4]。孫玉祥等采用58 mm火炮加載技術(shù)和多普勒探針系統(tǒng)測速技術(shù)對不同抗壓強(qiáng)度的混凝土進(jìn)行了平板撞擊實(shí)驗(yàn)的數(shù)值模擬[5]，確定了HJC本構(gòu)模型狀態(tài)方程參數(shù)。汪衡等研究了LS-DYNA軟件中HJC本構(gòu)模型的損傷參數(shù)Fs,并通過動(dòng)能彈侵徹深度和彈體的剩余速度進(jìn)行了校準(zhǔn)[6]。蘇偉林等通過HJC本構(gòu)模型研究了盾構(gòu)刀具切削混凝土材料時(shí)的阻力大小及變化規(guī)律[7]。張社榮等結(jié)合靜動(dòng)態(tài)力學(xué)實(shí)驗(yàn)確定了HJC本構(gòu)模型參數(shù)[8]。姚本余等通過靜力學(xué)實(shí)驗(yàn)確定了HJC本構(gòu)模型的基本參數(shù)[9]，采用LS-DYNA程序模擬了不同速度沖擊煤巖的破碎程度。凌天龍等通過靜力學(xué)實(shí)驗(yàn)與SHPB室內(nèi)實(shí)驗(yàn)確定了砂巖的HJC本構(gòu)模型的參數(shù)[10]，對室內(nèi)實(shí)驗(yàn)與數(shù)值模擬的應(yīng)力-應(yīng)變曲線進(jìn)行了對比。李成武等采用LS-DYNA有限元軟件對煤巖SHPB實(shí)驗(yàn)進(jìn)行了數(shù)值模擬,重現(xiàn)了室內(nèi)實(shí)驗(yàn)過程中的應(yīng)力波[11]。聞磊等采用HJC本構(gòu)模型研究了花崗斑巖SHPB動(dòng)態(tài)沖擊過程中的力學(xué)性能[12]。方秦等通過理論研究確定HJC本構(gòu)模型參數(shù)[13]，采用數(shù)值模擬研究了Salem石灰?guī)r的彈體侵徹實(shí)驗(yàn)。

目前，HJC本構(gòu)模型的研究多集中在混凝土沖擊破壞方面，也有部分學(xué)者將HJC本構(gòu)模型引入巖石材料的動(dòng)力學(xué)分析中，數(shù)值模擬效果較好。但對于我國西南地區(qū)常見巖石——石灰?guī)rHJC本構(gòu)模型參數(shù)的研究卻鮮有報(bào)道，同時(shí)在HJC本構(gòu)模型參數(shù)確定方法方面，基本通過傳統(tǒng)的室內(nèi)試驗(yàn)和經(jīng)驗(yàn)公式確定，不具有適用性及推廣性，無法還原石灰?guī)r的SHPB動(dòng)態(tài)沖擊破壞過程。

本文以石灰?guī)r靜力學(xué)試驗(yàn)、三軸圍壓實(shí)驗(yàn)、Hugoniot實(shí)驗(yàn)[14]、Hopkinson壓桿實(shí)驗(yàn)[15]、石灰?guī)r極限面參數(shù)A與帽蓋模型參數(shù)的關(guān)系為基礎(chǔ)，確定了石灰?guī)rHJC本構(gòu)模型參數(shù)，基于確定的石灰?guī)r模型參數(shù)，采用LS-DYNA有限元分析軟件再現(xiàn)了石灰?guī)rSHPB沖擊破壞的全過程，將該數(shù)值模擬巖石破壞過程中的應(yīng)力-應(yīng)變曲線與試驗(yàn)進(jìn)行對比分析，驗(yàn)證該HJC本構(gòu)模型取值的合理性及有效性，并以動(dòng)態(tài)峰值應(yīng)力為目標(biāo)函數(shù)，對21個(gè)參數(shù)進(jìn)行了敏感性分析，為石灰?guī)r動(dòng)力學(xué)分析提供一種簡單可行的數(shù)值模擬方法[16-18]。

1 HJC本構(gòu)模型

HJC本構(gòu)模型是由T J Holmqusit等為解決高應(yīng)變率及高壓荷載下的大變形問題提出的一種率相關(guān)本構(gòu)模型[19]，適用于拉格朗日與歐拉算法，模型主要包括屈服強(qiáng)度函數(shù)、狀態(tài)方程及損傷演化方程三部分。

屈服強(qiáng)度模型由標(biāo)準(zhǔn)化等效應(yīng)力表示

(1)

σ*=σ/f′c

(2)

圖 1 HJC屈服面方程Fig. 1 HJC yield surface equation

圖 2 HJC狀態(tài)方程Fig. 2 HJC state equation

損傷模型由等效塑性應(yīng)變和塑性體積應(yīng)變表征，演化方程為

(3)

(4)

圖 3 HJC損傷模型Fig. 3 HJC damage model

2 石灰?guī)rHJC本構(gòu)模型參數(shù)確定

2.1 靜力學(xué)參數(shù)確定

石灰?guī)r試樣取自桐梓隧道施工現(xiàn)場，為減小試樣的離散型，避免巖石試樣組成成分和結(jié)構(gòu)差異對試樣結(jié)果的影響，所有巖石試樣均取自同一巖塊。根據(jù)國際巖石力學(xué)學(xué)會(huì)(international society for rock mechanics，ISRM)對巖樣在力學(xué)試樣中尺寸的要求[17]，將靜力學(xué)巖石試樣加工打磨制成直徑50 mm、高100 mm的標(biāo)準(zhǔn)圓柱體巖樣，試件兩端面平整度誤差不超過0.02 mm，側(cè)面應(yīng)光滑筆直，滿足垂直度要求，如圖4所示。

圖 4 標(biāo)準(zhǔn)圓柱體巖樣Fig. 4 Standard cylindrical rock samples

靜力學(xué)試驗(yàn)采用中國科學(xué)院武漢巖土力學(xué)研究所研制的RMT-301巖石與混凝土力學(xué)試驗(yàn)系統(tǒng)開展，如圖5所示，由靜力學(xué)試驗(yàn)得到石灰?guī)r試樣的靜力學(xué)參數(shù)如表1所示。

圖 5 靜力學(xué)實(shí)驗(yàn)Fig. 5 Static experiment

表 1 石灰?guī)r靜力學(xué)參數(shù)

2.2 HJC本構(gòu)模型參數(shù)確定

根據(jù)Holmquist等提出的理論[19]，石灰?guī)r的強(qiáng)度大于48 MPa，在不考慮應(yīng)變率和損傷效應(yīng)的基礎(chǔ)上。

σ*=A+BP*N

(5)

根據(jù)帽蓋模型與HJC模型的關(guān)系式，由方秦等建立的特征化黏性強(qiáng)度參數(shù)與帽蓋模型參數(shù)之間的關(guān)系[13]。見圖6。

圖 6 帽蓋模型與HJC模型的關(guān)系Fig. 6 The relationship between the cap model and HJC model

AC面剪切破壞面函數(shù)為

(6)

CD面帽蓋面函數(shù)為

(7)

(8)

式中：J1為第一應(yīng)力不變量；J2為第二應(yīng)力偏量不變量；α、β、γ、θ、R為材料參數(shù)；帽蓋面與第一應(yīng)力不變量J1的交點(diǎn)表達(dá)式為

X(k)=k+RFf[L(k)]

(9)

(10)

(11)

根據(jù)Fossum等研究的試驗(yàn)數(shù)據(jù)可得[16]：α=270 MPa,γ=251 MPa，由靜力學(xué)實(shí)驗(yàn)得fC=105.82 MPa，求得A=0.31。

由三軸圍壓實(shí)驗(yàn)數(shù)據(jù)擬合

(12)

(13)

P*=P/fC=2σ1+σ3/3fC

(14)

根據(jù)在不同圍壓條件下的三軸壓縮結(jié)果，得到對應(yīng)的第一主應(yīng)力σ1，第三主應(yīng)力σ3，將式(13)、(14)進(jìn)行歸一化，得到每一個(gè)相對應(yīng)的(σ*,P*)，通過公式σ*=0.31+BP*N擬合曲線得到B和N的值，B=1.74、N=0.76，Sm取對應(yīng)于σ*不再增大時(shí)的值，Smax=15 GPa，可以得到HJC模型強(qiáng)度，見圖7，表達(dá)式為

σ*=0.31+1.74P*0.76

(15)

圖 7 參數(shù)B和N的取值方法Fig. 7 The determination method for the parameter of B and N

體積模量K=PC/μC，巖石的剪切模量G和體積模量K可由下式確定

G=E/[2(1+v)]

(16)

K=E/[3(1-2v)]

(17)

式中：E為彈性模量；v為泊松比，見表1。由公式(16)、(17)求得G=7.65 GPa、K=12.75 GPa。PC、μC為單軸壓縮實(shí)驗(yàn)測得的靜水壓力和體積應(yīng)變，PC=FC/3=0.035 GPa[1]、μC=PC/K=2.7×10-4。μ1=ρ/ρ0-1，ρ為材料的即時(shí)密度，ρ0為材料的初始密度，由實(shí)驗(yàn)數(shù)據(jù)得μ1=0.163。

圖 8 參數(shù)K1、K2、K3取值方法Fig. 8 The determination method for the Parameter of K1,K2,K3

通過SHPB壓桿實(shí)驗(yàn)，得到了不同應(yīng)變率下石灰?guī)r單軸動(dòng)態(tài)強(qiáng)度，結(jié)合準(zhǔn)靜態(tài)單軸壓縮強(qiáng)度，得到不同應(yīng)變率下的特征化強(qiáng)度σ*，通過擬合應(yīng)變率與無量綱特征化強(qiáng)度的關(guān)系，擬合斜率C=0.017。見圖9。

圖 9 參數(shù)C取值方法Fig. 9 Parameter C value method

Holmquist指出損傷參數(shù)與材料強(qiáng)度無關(guān)[1]，損傷參數(shù)D1根據(jù)文獻(xiàn)[19]公式D1=0.01/(1/6+T*)=確定，其中T*=T/FC，求得D1=0.045，D2=1、εfmin=0.01。HJC模型參數(shù)如表2所示。

表 2 石灰?guī)rHJC模型參數(shù)取值

3 參數(shù)驗(yàn)證

3.1 模型建立

SHPB動(dòng)態(tài)沖擊數(shù)值模型的建立采用Hypermesh建立，子彈長度400 mm，入射桿和透射桿長度為2000 mm，子彈、入射桿及透射桿材料模型采用LS-DYNA軟件中001號(hào)MAT_ELASTIC本構(gòu)模型，直徑均為50 mm，材料模型選擇LS-DYNA軟件中111號(hào)MAT-JOHNSON-HOLMQUIST-CONCRETE本構(gòu)模型。數(shù)值計(jì)算模型如圖10所示。

圖 10 SHPB數(shù)值模擬模型圖Fig. 10 SHPB numerical simulation model

采用表2的HJC本構(gòu)模型參數(shù)進(jìn)行數(shù)值模擬計(jì)算，為確保巖石計(jì)算的精確，對巖石的橫向、環(huán)線網(wǎng)格均進(jìn)行了精細(xì)化劃分，通過在k文件中定義材料、子彈、入射桿、透射桿屬性，材料、子彈、入射桿、透射桿均采用面面接觸(CONTACT_AUTOMATIC_SURFACE_TO_SURFACE)，在k文件INITIAL_VELOCITY_GENERATION對子彈水平方向添加了和試驗(yàn)對應(yīng)10 m/s和12 m/s的沖擊速度兩次沖擊入射桿，并在k文件中添加關(guān)鍵字MAT_ADD_EROSION表征巖石的拉伸失效準(zhǔn)則。

3.2 SHPB試驗(yàn)

本實(shí)驗(yàn)采用阿基米德工業(yè)科技有限公司研制的分離式Hopkinson壓桿測試系統(tǒng)ALT100進(jìn)行沖擊動(dòng)力學(xué)實(shí)驗(yàn)，如圖11所示。壓桿密度為7.81 g/cm3，彈性模量為210 GPa，泊松比為0.28，縱波波速為5410 m/s，數(shù)據(jù)采集系統(tǒng)由黏貼在壓桿上的應(yīng)變片、惠斯通橋路(應(yīng)變片接線橋盒)、超動(dòng)態(tài)應(yīng)變儀以及高速采集系統(tǒng)組成。

圖 11 分離式 Hopkinson 壓桿測試系統(tǒng) ALT100Fig. 11 Split Hopkinson pressure bar test system ALT100

3.3 結(jié)果分析

3.3.1 應(yīng)力-應(yīng)變曲線

(18)

(19)

(20)

式中：A、L分別為試樣的橫截面積和長度；E、A0分別為壓桿的彈性模量和橫截面積;C為壓桿的縱波波速。

根據(jù)圖12可以看出，巖石動(dòng)態(tài)破壞的峰值應(yīng)力隨沖擊速度的增大而增大，基于上述確定的參數(shù)，10 m/s和 12 m/s的不同速度沖擊時(shí)，數(shù)值模擬得到的應(yīng)力-應(yīng)變曲線與室內(nèi)試驗(yàn)應(yīng)力-應(yīng)變曲線都基本吻合，但由于數(shù)值模擬默認(rèn)巖石為各向同性材料，數(shù)值模擬得到的曲線為一條光滑的曲線，而實(shí)際試驗(yàn)的巖石為各向異性，巖石內(nèi)部結(jié)構(gòu)的不規(guī)整性，巖石破壞過程表現(xiàn)不均勻，導(dǎo)致試驗(yàn)得到的曲線出現(xiàn)凹凸現(xiàn)象。

圖 12 數(shù)值模擬曲線與試驗(yàn)應(yīng)力應(yīng)變圖Fig. 12 Numerical simulation curves and experimental stress-strain curve

3.3.2 破壞過程

試驗(yàn)過程利用FASTCAM SA1.1型高速攝影儀記錄石灰?guī)r在10 m/s速度動(dòng)態(tài)沖擊作用下裂紋的擴(kuò)展過程及最終破壞過程，如圖13所示。圖14根據(jù)相對應(yīng)10 m/s沖擊速度數(shù)值模擬得到的裂紋擴(kuò)展情況及其破壞過程，裂紋的擴(kuò)展方向大部分沿軸線60°方向，由于應(yīng)力波在巖石內(nèi)部多次透射反射，巖石側(cè)面出現(xiàn)拉伸破壞裂紋，并在動(dòng)態(tài)壓縮作用下壓縮破壞，在拉伸和壓縮裂紋共同作用下產(chǎn)生破壞。

根據(jù)試驗(yàn)與數(shù)值模擬分析可得，數(shù)值模擬與SHPB動(dòng)態(tài)沖擊試驗(yàn)的應(yīng)力-應(yīng)變曲線、破壞模式基本吻合。

4 石灰?guī)rHJC參數(shù)敏感性分析

影響沖擊后巖石的峰值應(yīng)力σmax i參數(shù)變化率δFij表達(dá)式為

(i=1，2，3,…,n；j=1，2，3,…,n)

(21)

圖 13 巖石試驗(yàn)破壞過程Fig. 13 Failure process of rock test

δFij(j=1，2，3，4)對應(yīng)-40%、-20%、20%、40%。

峰值應(yīng)力變化δσmax ij表達(dá)式為

(i=1，2，3,…,n；j=1，2，3，4)

(22)

式中：σmax ij表示HJC本構(gòu)模型第i個(gè)參數(shù)在第j個(gè)參數(shù)變化率狀態(tài)下模擬所得的峰值應(yīng)力。

圖 14 巖石數(shù)值模擬破壞過程Fig. 14 Numerical simulation failure process of rock

根據(jù)在LS/DYNA軟件中SHPB數(shù)值模擬采用HJC本構(gòu)模型對21個(gè)參數(shù)調(diào)整所計(jì)算的峰值應(yīng)力變化率結(jié)果如圖15所示。

圖 15 峰值應(yīng)力對HJC參數(shù)敏感性分析Fig. 15 Sensitivity analysis of peak stress to HJC parameters

由圖可得，根據(jù)文獻(xiàn)[18]，采用剩余峰值應(yīng)力變化率絕對值|δσmax ij|=10%為參數(shù)敏感性分析標(biāo)準(zhǔn)，即巖石峰值應(yīng)力變化率|δσmax ij|≥10%即認(rèn)為該參數(shù)對峰值應(yīng)力影響較大，參數(shù)A、B、N、FC對峰值應(yīng)力變化率影響較大，B值的變化使峰值應(yīng)力變化率接近60%，其他參數(shù)對峰值應(yīng)力變化率在10%之間，說明參數(shù)A、B、N、FC對石灰?guī)r破壞時(shí)的峰值應(yīng)力是敏感的。

5 結(jié)論

(1)基于確定的極限面參數(shù)A與帽蓋模型參數(shù)α與γ之間的關(guān)系，求得A=0.31；結(jié)合三軸圍壓實(shí)驗(yàn)與公式σ*=A+BP*N擬合曲線得到B=1.74、N=0.76。

(2)基于確定的石灰?guī)rHJC模型參數(shù)，分析了石灰?guī)rSHPB動(dòng)態(tài)沖擊數(shù)值模擬結(jié)果與SHPB動(dòng)態(tài)沖擊試驗(yàn)應(yīng)力應(yīng)變曲線及破壞形態(tài)相似，巖石破壞時(shí)的峰值應(yīng)力相差均小于10%，驗(yàn)證了參數(shù)確定方法的準(zhǔn)確性。

(3)以石灰?guī)r動(dòng)態(tài)峰值應(yīng)力為目標(biāo)函數(shù)，對HJC模型21個(gè)參數(shù)進(jìn)行了敏感性分析，HJC本構(gòu)模型參數(shù)fc、A、B、N對石灰?guī)r動(dòng)態(tài)峰值應(yīng)力的影響大于10%，B值的變化使峰值應(yīng)力變化率達(dá)到60%，因此石灰?guī)r動(dòng)態(tài)峰值應(yīng)力對B值的變化最為敏感。