高應(yīng)變率下斷裂韌性實(shí)驗(yàn)的數(shù)值模擬*

葉 波，巫緒濤，胡鳳輝，廖 禮

(合肥工業(yè)大學(xué)土木與水利工程學(xué)院，安徽 合肥 230009)

高應(yīng)變率下斷裂韌性實(shí)驗(yàn)的數(shù)值模擬*

葉 波，巫緒濤，胡鳳輝，廖 禮

(合肥工業(yè)大學(xué)土木與水利工程學(xué)院，安徽 合肥 230009)

采用有限元軟件ANSYS/LS-DYNA程序?qū)o態(tài)和沖擊荷載作用下的含裂紋半圓彎曲(SCB)實(shí)驗(yàn)進(jìn)行了數(shù)值模擬。根據(jù)靜態(tài)實(shí)驗(yàn)的模擬結(jié)果，提出了適合復(fù)合型加載的Ⅰ型應(yīng)力強(qiáng)度因子擬合公式，采用該公式計(jì)算應(yīng)力強(qiáng)度因子的最大誤差不超過10%。動(dòng)態(tài)實(shí)驗(yàn)的模擬結(jié)果表明：對(duì)于純Ⅰ型加載的SCB實(shí)驗(yàn)，動(dòng)態(tài)應(yīng)力強(qiáng)度因子隨著試樣半徑、支座間距以及相對(duì)裂紋長(zhǎng)度的變化呈現(xiàn)規(guī)律性變化；當(dāng)試樣半徑小于60 mm、相對(duì)支座間距為1.2、相對(duì)裂紋長(zhǎng)度在0.1～0.4范圍內(nèi)時(shí)，慣性效應(yīng)的影響較小，采用靜態(tài)擬合公式計(jì)算裂尖的動(dòng)態(tài)應(yīng)力強(qiáng)度因子的誤差約10%；對(duì)于復(fù)合型加載的SCB實(shí)驗(yàn)，當(dāng)相對(duì)裂紋長(zhǎng)度為0.2～0.4、裂紋傾角在10°～40°范圍內(nèi)時(shí)，采用靜態(tài)擬合公式計(jì)算裂尖的動(dòng)態(tài)應(yīng)力強(qiáng)度因子的誤差小于10%。

固體力學(xué)；斷裂韌性；沖擊荷載；應(yīng)力強(qiáng)度因子；含裂紋SCB實(shí)驗(yàn)；復(fù)合型加載

混凝土、巖石類脆性材料在土木工程中大量使用，對(duì)于一些重要的防護(hù)結(jié)構(gòu)，如大壩、核電站安全殼等，設(shè)計(jì)時(shí)必須考慮強(qiáng)沖擊載荷作用。而這些材料往往含有大量的原生缺陷，在沖擊載荷作用下的動(dòng)態(tài)響應(yīng)與靜載荷下顯著不同，因此研究其動(dòng)態(tài)力學(xué)性能，特別是與裂紋擴(kuò)展相關(guān)的動(dòng)態(tài)斷裂韌性十分重要。

由于理論和實(shí)驗(yàn)技術(shù)的困難，高應(yīng)變率下脆性材料動(dòng)態(tài)斷裂行為的研究一直進(jìn)展緩慢。M.Nakano等[1]和Y.Yamauchi等[2]在分離式Hopkinsion壓桿(SHPB)上采用含中心直裂紋巴西圓盤試樣(CCBD)進(jìn)行了復(fù)合型動(dòng)態(tài)斷裂韌性實(shí)驗(yàn)。由于實(shí)驗(yàn)技術(shù)的不成熟，得到的實(shí)驗(yàn)數(shù)據(jù)較離散，但仍然發(fā)現(xiàn)結(jié)果與靜態(tài)實(shí)驗(yàn)的顯著不同，不僅受實(shí)驗(yàn)應(yīng)變率的影響，同時(shí)與試樣大小和裂紋長(zhǎng)度密切相關(guān)。為了更好地測(cè)量脆性材料的斷裂韌性，有學(xué)者采用含裂紋或缺口的半圓彎曲(SCB)試樣代替?zhèn)鹘y(tǒng)的CCBD試樣，并結(jié)合有限元模擬得到靜態(tài)斷裂韌性的計(jì)算公式[3-5]。R.Chen等[6]和F.Dai等[7-8]將SCB試樣引入脆性材料動(dòng)態(tài)斷裂韌性的測(cè)量，并采用先進(jìn)的激光測(cè)量系統(tǒng)測(cè)量試樣裂紋或缺口的張開位移，假設(shè)實(shí)驗(yàn)中試樣達(dá)到應(yīng)力平衡，忽略慣性效應(yīng)的影響，沿用靜態(tài)分析結(jié)果計(jì)算動(dòng)態(tài)斷裂韌性。然而在SHPB實(shí)驗(yàn)中，大多數(shù)脆性材料很難達(dá)到應(yīng)力平衡，慣性效應(yīng)的影響始終存在，上述研究中沒有明確提出動(dòng)態(tài)加載下靜態(tài)分析的適用范圍和可靠性。

本文中采用有限元方法，首先對(duì)靜態(tài)含裂紋SCB實(shí)驗(yàn)進(jìn)行數(shù)值模擬，結(jié)合量綱分析推導(dǎo)適用于復(fù)合型斷裂的Ⅰ型應(yīng)力強(qiáng)度因子(KⅠ)擬合公式；對(duì)基于SHPB技術(shù)的動(dòng)態(tài)復(fù)合型SCB實(shí)驗(yàn)進(jìn)行數(shù)值模擬，計(jì)算試樣裂尖的動(dòng)態(tài)應(yīng)力強(qiáng)度因子(KⅠd)，分析試樣半徑、支座間距、相對(duì)裂紋長(zhǎng)度及裂紋傾角對(duì)斷裂韌性的影響，討論靜態(tài)公式的適用范圍。

1 靜態(tài)復(fù)合型加載SCB實(shí)驗(yàn)?zāi)M及Ⅰ型應(yīng)力強(qiáng)度因子擬合公式

1.1 含裂紋SCB實(shí)驗(yàn)

圖1 靜態(tài)復(fù)合型加載SCB實(shí)驗(yàn)簡(jiǎn)圖Fig.1 Diagram of static SCB test under mixed mode loading

采用有限元軟件ANSYS程序?qū)o態(tài)復(fù)合型斷裂SCB實(shí)驗(yàn)進(jìn)行了數(shù)值模擬，實(shí)驗(yàn)示意圖如圖1所示，其中：a為裂紋長(zhǎng)度，R為試樣半徑，F(xiàn)試樣頂部集中載荷，β為裂紋傾角(裂紋與試樣中心軸間的夾角)，s為半支座間距。

試樣采用線彈性本構(gòu)，彈性模量為30 GPa，密度2 500 kg/m3，泊松比0.2。試樣半徑為50 mm，厚25 mm，支座間距為60 mm，裂紋相對(duì)長(zhǎng)度a/R從0.01變化至0.8，裂紋傾角β從0°變化至45°。

依據(jù)圖1建立有限元模型，單元類型為solid 186六面體實(shí)體單元。為了更好地反映裂紋應(yīng)力奇異的影響，圍繞裂尖一圈采用命令KSCON實(shí)現(xiàn)退化奇異等參元。應(yīng)力強(qiáng)度因子采用相互作用積分法[9]計(jì)算，該法是基于J積分的一種新型應(yīng)力強(qiáng)度因子計(jì)算方法，其精度與數(shù)值外插法相當(dāng)。對(duì)于靜態(tài)問題，ANSYS提供了支持該法的命令流，適于批量操作，對(duì)于而動(dòng)態(tài)問題則沒有，因此本文中采用不同的方法計(jì)算靜態(tài)和動(dòng)態(tài)應(yīng)力強(qiáng)度因子。

1.2 Ⅰ型應(yīng)力強(qiáng)度因子擬合公式

根據(jù)量綱分析并結(jié)合文獻(xiàn)[3-8]進(jìn)行計(jì)算，選擇Ⅰ型應(yīng)力強(qiáng)度因子KⅠ擬合公式形式為：

(1)

式中：B為試樣厚度，g、h為擬合函數(shù)。式(1)的前半部滿足應(yīng)力強(qiáng)度因子的量綱形式，且滿足當(dāng)a趨于0時(shí)，公式退化為彎曲梁?jiǎn)芜呅绷鸭y的KI公式。后半部為量綱一形狀因子，反映了裂紋長(zhǎng)度a、支座間距s和裂紋傾角β的影響，擬合函數(shù)g、h采用多項(xiàng)式形式，表達(dá)式為：

式中：dn和en(n=1,2,…,5)為裂紋傾角β的擬合系數(shù)，均采用直線方程擬合：

(4)

圖2 靜態(tài)復(fù)合型加載SCB實(shí)驗(yàn)應(yīng)力強(qiáng)度因子擬合結(jié)果Fig.2 Fitted effect of stress intensity factor in static SCB test under mixed mode loading

采用本文公式模擬不同傾角β下的應(yīng)力強(qiáng)度因子，結(jié)果如圖2所示。最大擬合誤差小于10%。從圖2可以看出KⅠ值隨裂紋傾角β增加而降低，當(dāng)傾角大于40°時(shí)，KⅠ值趨向0，反映出斷裂由純Ⅰ型向復(fù)合型直至純Ⅱ型轉(zhuǎn)變的過程。

2 基于SHPB裝置的動(dòng)態(tài)SCB實(shí)驗(yàn)的模擬

2.1 有限元模型和動(dòng)態(tài)應(yīng)力強(qiáng)度因子計(jì)算

圖3 基于SHPB裝置的動(dòng)態(tài)SCB實(shí)驗(yàn)有限元模型Fig.3 Finite element model of dynamic SCB test based on SHPB device

圖4 三角形速度脈沖Fig.4 Triangular velocity pulse

模擬采用LS-DYNA程序，其中試樣的材料常數(shù)與靜態(tài)完全相同，SHPB裝置的入射桿和透射桿采用線彈性材料：彈性模量200 GPa，密度7 850 kg/m3，泊松比0.3。直徑100 mm，長(zhǎng)度2 000 mm。所用單元類型均采用Solid 164六面體實(shí)體單元。有限元模型如圖3所示。載荷為加在左端入射桿端面的三角形速度脈沖，如圖4所示。

圖5 不同相對(duì)裂紋長(zhǎng)度下的應(yīng)力強(qiáng)度因子變化規(guī)律Fig.5 Variation of stress intensity factor for different relative crack lengths

由于動(dòng)態(tài)模擬無法實(shí)現(xiàn)退化奇異單元，在裂尖附近采用稠密單元，減小單元網(wǎng)格尺寸影響。采用2種方法計(jì)算裂尖的動(dòng)態(tài)應(yīng)力強(qiáng)度因子KⅠd：(1) 由裂紋面張開位移通過數(shù)值外插法[10]得到，記為KⅠd1；(2) 根據(jù)試樣左右端面的接觸力-時(shí)間曲線平均后帶入式(1)計(jì)算得到，記為KⅠd2。KⅠd1為裂尖實(shí)際動(dòng)態(tài)應(yīng)力強(qiáng)度因子，KⅠd2為采用靜態(tài)公式得到的動(dòng)態(tài)應(yīng)力強(qiáng)度因子。

2.2 動(dòng)態(tài)純Ⅰ型加載SCB實(shí)驗(yàn)的模擬及分析

首先對(duì)裂紋傾角β=0°的純Ⅰ型加載SCB實(shí)驗(yàn)進(jìn)行了模擬，研究相對(duì)裂紋長(zhǎng)度a/R、試樣半徑R、以及支座間距s對(duì)裂尖實(shí)際動(dòng)態(tài)應(yīng)力強(qiáng)度因子KⅠd1的影響。圖5所示為固定試樣半徑R=50 mm、相對(duì)支座間距s/R=0.6、裂紋相對(duì)長(zhǎng)度a/R從0.1變化至0.8時(shí)的KⅠd1。圖6所示為a/R=0.4，s/R=1.2、試樣半徑R從40 mm變化至70 mm時(shí)的KⅠd1。圖7所示為試樣半徑R=50 mm、a/R=0.4、相對(duì)支座間距s/R從0.4變化至0.8時(shí)的KⅠd1。

圖6 不同試樣半徑下應(yīng)力強(qiáng)度因子的變化規(guī)律Fig.6 Variation of stress intensity factor for different sample radiuses

圖7 不同支座間距下應(yīng)力強(qiáng)度因子的變化規(guī)律Fig.7 Variation of stress intensity factor for different distances between two supports

由圖5～7可以看出：

(1)KⅠd1隨裂紋相對(duì)長(zhǎng)度a/R、支座間距s的增大而增大，隨試樣半徑R的增大而減小。

(2)KⅠd1-t曲線在加載過程中均單調(diào)增加，而達(dá)到峰值后出現(xiàn)震蕩；隨a/R、試樣半徑R增加，震蕩加劇，較大和較小的支座間距也會(huì)導(dǎo)致震蕩加劇。這種震蕩反映了試樣應(yīng)力不平衡導(dǎo)致的慣性效應(yīng)的影響。當(dāng)a/R<0.4、試樣半徑小于60 mm、s/R≈1.2時(shí)，試樣應(yīng)力不平衡影響較小。

為了進(jìn)一步分析靜態(tài)公式的適用范圍，對(duì)a/R在0.1～0.6范圍內(nèi)裂尖實(shí)際KⅠd1和靜態(tài)公式得到的KⅠd2進(jìn)行了比較，結(jié)果如圖8所示，從圖中可以看出，兩者在a/R=0.4時(shí)吻合最好。按下式得到不同a/R情況下兩者的相對(duì)誤差δ并繪圖，如圖9所示，其中δ的表達(dá)式為

(5)

從圖9可以看出，當(dāng)a/R<0.4時(shí)，δ趨近于0，兩者相對(duì)誤差在10%左右，即此時(shí)慣性效應(yīng)的影響較小，靜態(tài)公式具有足夠的精度。

圖8 KⅠd1和KⅠd2的比較(a/R=0.4)Fig.8 Comparison of KⅠd1 and KⅠd2 (a/R=0.4)

圖9 靜態(tài)公式的相對(duì)誤差隨a/R的變化Fig.9 Relative error of static formula for different a/R

2.3 動(dòng)態(tài)復(fù)合型加載SCB實(shí)驗(yàn)?zāi)M及分析

為了研究復(fù)合型加載SCB實(shí)驗(yàn)中KⅠd的變化規(guī)律，在純I型研究的基礎(chǔ)上，固定試樣半徑R=50 mm、支座間距s=1.2R，裂紋傾角β從10°變化至40°，每個(gè)傾角下相對(duì)裂紋長(zhǎng)度a/R從0.2變化至0.8建立有限元模型進(jìn)行研究。為了便于分析，取裂尖KⅠd1-t曲線加載段的峰值KⅠd1max進(jìn)行比較。圖10所示為不同a/R情況下KⅠd1max隨裂紋傾角β的變化規(guī)律，圖11所示為不同裂紋傾角β情況下KⅠd1max隨a/R的變化規(guī)律。

圖10 KⅠd1max隨裂紋傾角β的變化規(guī)律Fig.10 Variation of KⅠd1max with crack angle β

圖11 KⅠd1max隨裂紋相對(duì)長(zhǎng)度a/R的變化規(guī)律Fig.11 Variation of KⅠd1max with relative crack length a/R

由圖10可以發(fā)現(xiàn)，對(duì)應(yīng)相同的裂紋相對(duì)長(zhǎng)度，KⅠd1max隨裂紋傾角β的增大而減小，且a/R越大，下降速率越大。由圖11可以發(fā)現(xiàn)，對(duì)應(yīng)相同裂紋傾角，KⅠd1max隨裂紋相對(duì)長(zhǎng)度a/R的增大而增大，而當(dāng)裂紋傾角β達(dá)到40°時(shí)，KⅠd1max隨裂紋相對(duì)長(zhǎng)度a/R的增大而減小，變化規(guī)律與靜態(tài)模擬結(jié)果類似，如圖2所示。

為了進(jìn)一步分析復(fù)合加載條件下靜態(tài)擬合公式的適用性，按式(5)計(jì)算了裂尖實(shí)際應(yīng)力強(qiáng)度因子峰值KⅠd1max與靜態(tài)公式得到的應(yīng)力強(qiáng)度因子峰值KⅠd2max之間的相對(duì)誤差δm，如表1所列?？梢钥闯?，當(dāng)裂紋傾角β在10°～40°之間時(shí)，裂紋相對(duì)長(zhǎng)度a/R在0.2～0.4之間時(shí)，靜態(tài)公式的相對(duì)誤差小于10%，精度較高，而超出此范圍，相對(duì)誤差增大至20%～50%，即靜態(tài)公式失效。

表1 復(fù)合型加載下靜態(tài)公式計(jì)算KⅠdmax的相對(duì)誤差δmTable 1 Relative error δmof KⅠdmax in the static formula under mixed-mode loading

3 結(jié) 論

本文中通過有限元分析方法，對(duì)靜態(tài)和動(dòng)態(tài)復(fù)合型加載SCB實(shí)驗(yàn)進(jìn)行了數(shù)值模擬研究，得到如下結(jié)論：

(1) 根據(jù)靜態(tài)實(shí)驗(yàn)的模擬結(jié)果，結(jié)合量綱分析得到了適合于裂紋傾角在0°～45°，裂紋相對(duì)長(zhǎng)度在0.1～0.8范圍內(nèi)的應(yīng)力強(qiáng)度因子擬合公式，計(jì)算結(jié)果最大誤差小于10%。

(2) 對(duì)于含純Ⅰ型中心裂紋的動(dòng)態(tài)SCB實(shí)驗(yàn)，試樣裂尖的動(dòng)態(tài)斷裂韌性隨試樣半徑增大而減?。浑S支座間距和裂紋相對(duì)長(zhǎng)度的增大而增大。

(3) 對(duì)于含復(fù)合型裂紋的動(dòng)態(tài)SCB實(shí)驗(yàn)，當(dāng)裂紋相對(duì)長(zhǎng)度在0.2～0.7范圍內(nèi)時(shí)，裂尖的動(dòng)態(tài)斷裂韌性隨著裂紋傾角的增大而減??；當(dāng)裂紋傾角在10°～30°范圍內(nèi)時(shí)，裂尖的動(dòng)態(tài)斷裂韌性隨著裂紋相對(duì)長(zhǎng)度的增大而增大。

(4) 通過采用位移外插法和靜態(tài)擬合公式2種方法計(jì)算裂尖的動(dòng)態(tài)斷裂韌性，并計(jì)算相對(duì)誤差發(fā)現(xiàn)，對(duì)于含純Ⅰ型中心裂紋的動(dòng)態(tài)SCB實(shí)驗(yàn)，當(dāng)裂紋相對(duì)長(zhǎng)度在0.1～0.4范圍內(nèi)時(shí)，可以忽略慣性效應(yīng)，直接采用靜態(tài)KI公式計(jì)算裂尖的最大斷裂韌性，相對(duì)誤差約10%。對(duì)于含復(fù)合型斜裂紋的動(dòng)態(tài)SCB實(shí)驗(yàn)，當(dāng)裂紋相對(duì)長(zhǎng)度在0.2～0.4，裂紋傾角在10°～40°范圍內(nèi)時(shí)，也可以直接采用靜態(tài)公式計(jì)算裂尖的最大斷裂韌性，相對(duì)誤差小于10%。

[1] Nakano M, Kashida K, Yamauchi Y, et al. Dynamic fracture initiation in brittle materials under combined mode Ⅰ/Ⅱ loading[J]. Journal De Physique IV, 1994,4(8):695-700.

[2] Yamauchi Y, Nakano M, Kishida K, et al.Measurement of fracture toughness for brittle materials under mixed-mode impact loading using center-notched disk specimen[J]. Journal of the Society of Materials Science Japan,2000,49(12):1324-1329.

[3] Ayatollahi M R, Aliha M R M, Hassani M M. Mixed mode brittle fracture in PMMA-An experimental study using SCB specimens[J]. Materials Science and Engineering A, 2006,417(1/2):348-356.

[4] Ayatollahi M R, Aliha M R M. On determination of mode Ⅱ fracture toughness using semi-circular bend specimen[J]. International Journal of Solids and Structures, 2006,43(17):5217-5227.

[5] Aliha M R M, Ayatollahi,M R. Mixed mode Ⅰ/Ⅱ brittle fracture evaluation of marble using SCB specimen[J]. Procedia Engineering, 2011,10:311-318.

[6] Chen R, Xia K W, Dai F et al. Determination of dynamic fracture parameters using a semi-circular bend technique in split Hopkinson pressure bar testing[J]. Engineering Fracture Mechanics, 2009,76(9):1268-1276.

[7] Dai F, Xia K W, et al. Rate dependence of the flexural tensile strength of Laurentian granite[J]. International Journal of Rock Mechanics and Mining Sciences, 2010,47(3):469-475.

[8] Dai F, Xia K W, et al. Determination of dynamic rock Mode-I fracture parameters using cracked chevron notched semi-circular bend specimen[J]. Engineering Fracture Mechanics, 2011,78(15):2633-2644.

[9] Yau J, Wang S, Corten H. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity[J]. Journal of Applied Mechanics, 1980,47(2):335-341.

[10] 巫緒濤,楊伯源.數(shù)值外插法求解空間裂紋應(yīng)力強(qiáng)度因子的研究[J].合肥工業(yè)大學(xué)學(xué)報(bào),1999,22(4):26-31. Wu Xutao, Yang Boyuan. Study on the numerical extrapolated method for calculating the 3-D crack front stress intensity factor[J]. Journal of Hefei University of Tchnology, 1999,22(4):26-31.

(責(zé)任編輯 王小飛)

Numerical simulation of fracture toughness test under high strain rate

Ye Bo, Wu Xutao, Hu Fenghui, Liao Li

(SchoolofCivilEngineering,HefeiUniversityofTechnology,Hefei230009,Anhui,China)

In this work we conducted numerical simulations of semi-circular bending (SCB) test of steel with cracks respectively in a static condition and under impact loading, using ANSYS/LS-DYNA, a finite element software. According to the results achieved from the simulation of the static test, we have put forward the formula for the calculation of Steel Ⅰ’s fitted stress intensity factor which, when applied to the calculation of the factor under the mixed-mode loading, can keep the maximum error below 10%. The simulation result of dynamic test shows that: for semi-circular bending (SCB) test of pure Steel Ⅰ under impact loading, the dynamic stress intensity factor exhibits regular variation with the change in the specimen’s radius, the distance between the supports, and the relative crack length: When the specimen’s radius is less than 60 mm, the distance between two supports is 1.2 and the relative crack length is in the range of 0.1 to 0.4, the inertial effect is relatively small and the error of the dynamic stress intensity factor calculated with the above formula is about 10%; when the relative crack length is 0.2 to 0.4 and the crack angle is in the range of 10° to 40°, the error of KId which is calculated with this formula is less than 10% for semi-circular bending (SCB) test under mixed-mode loading.

solid mechanics; fracture toughness; impact loading; stress intensity factor; SCB test with crack; mixed-mode loading

10.11883/1001-1455(2016)03-0416-06

2014-10-13；

2015-02-03

國(guó)家自然科學(xué)基金項(xiàng)目(11072072)

葉 波(1991- )，男，碩士研究生；

巫緒濤，wuxvtao@sina.com。

O347.3國(guó)標(biāo)學(xué)科代碼：13015

A