考慮拉伸斷裂的非線性本構(gòu)模型

考慮拉伸斷裂的非線性本構(gòu)模型

1本構(gòu)模型的建立mohr-coulumb（m-c）的投降模式是剪切強(qiáng)度模式，在數(shù)值模擬中得到了廣泛應(yīng)用。FLAC中的M-C屈服準(zhǔn)則又考慮了拉伸截?cái)?因而也能較好地反映材料的拉破壞特性。當(dāng)應(yīng)力小于10MPa時(shí),屈服函數(shù)可近似為一直線,因而M-C本構(gòu)模型適用,其中凝聚力c及內(nèi)摩擦角?均為常量。而大量試驗(yàn)結(jié)果表明,對(duì)于堅(jiān)硬巖體且應(yīng)力、屈服函數(shù)是非線性的。隨著采深的增加,巷道圍巖所受的應(yīng)力越來越大,c及?均為常量的M-C屈服準(zhǔn)則的適用性變差。Hoek-Brown(H-B)屈服準(zhǔn)則是一種非線性的屈服準(zhǔn)則。文獻(xiàn)采用凝聚力弱化-脆性-內(nèi)摩擦角強(qiáng)化的H-B本構(gòu)模型及常偏應(yīng)力準(zhǔn)則模擬了圓形巷道圍巖的應(yīng)力及破壞區(qū)分布;文獻(xiàn)采用彈-脆-塑性的H-B應(yīng)變軟化模型模擬了橢圓形巷道圍巖的破壞區(qū)分布。但他們均未能模擬出巖爆后圍巖中形成的V型坑。在FLAC中,H-B本構(gòu)模型不能顯示出單元的拉破壞區(qū)域;而且,當(dāng)巖石發(fā)生脆性破壞后,計(jì)算出的殘余抗拉強(qiáng)度有時(shí)會(huì)大于初始抗拉強(qiáng)度,這顯然與實(shí)際情況不符。FLAC提供了所有本構(gòu)模型的C++源代碼,為用戶提供了二次開發(fā)的平臺(tái)。文獻(xiàn)[5–7]已經(jīng)成功完成了鄧肯-張模型、黏彈性流變模型及改進(jìn)的Burgers蠕變損傷模型的開發(fā)。本文以FLAC中的M-C本構(gòu)模型為藍(lán)本,進(jìn)行了非線性本構(gòu)模型的二次開發(fā),自定義本構(gòu)模型既能考慮拉伸截?cái)?又能反映屈服函數(shù)的非線性特征。通過算例驗(yàn)證了自定義本構(gòu)模型的正確性,并將其應(yīng)用于圓形巷道圍巖的應(yīng)變局部化過程的研究。2flac的第二次開發(fā)2.1m-c本構(gòu)模型M-C屈服面在最大主應(yīng)力σ1與最小主應(yīng)力σ3平面內(nèi)的投影為一直線,即巖石破壞時(shí)σ1與σ3呈線性關(guān)系,凝聚力c及內(nèi)摩擦角?均為常量。H-B屈服函數(shù)在σ1-σ3平面內(nèi)是非線性的,將σ1對(duì)σ3求偏導(dǎo),即可得到屈服函數(shù)上任意一點(diǎn)的切線的斜率及其在σ1軸上的截距。它們均與M-C屈服函數(shù)中的c及?有關(guān),此時(shí)由M-C屈服函數(shù)求得的c及?均為σ3的函數(shù),分別用c(σ3)及?(σ3)表示。FLAC中的H-B本構(gòu)模型不能顯示出單元的拉破壞區(qū)域,而M-C本構(gòu)模型則可以,并能考慮拉伸截?cái)?。圖1給出了本文提出的既能考慮拉伸截?cái)?又能反映巖石破壞時(shí)σ1與σ3之間的非線性關(guān)系的屈服函數(shù)(粗實(shí)線)?？紤]到許多巖石的峰后脆性極強(qiáng),本文中的初始屈服函數(shù)一旦被滿足,就用殘余本構(gòu)參數(shù)替換初始本構(gòu)參數(shù),巖石由初始強(qiáng)度跌至殘余強(qiáng)度。在圖1中,粗實(shí)線(1)為初始屈服函數(shù);粗實(shí)線(2)為殘余屈服函數(shù);σt為巖石的初始抗拉強(qiáng)度;σtres為巖石的殘余抗拉強(qiáng)度。圖1中的初始屈服函數(shù)及殘余屈服函數(shù)均由3部分構(gòu)成。當(dāng)σ3<0時(shí),屈服函數(shù)為H-B屈服函數(shù);當(dāng)0<σ3<σt或0<σ3<σtres時(shí),屈服函數(shù)為M-C屈服函數(shù),為了保證屈服函數(shù)的連續(xù)性,在σ3=0點(diǎn)處做H-B屈服函數(shù)的切線(σ3<0部分的虛線),利用斜直線替換σ3>0時(shí)的虛線部分,此時(shí)M-C屈服函數(shù)中的c及?均為常量,即c(0)及?(0);當(dāng)σ3=σt或σ3=σtres時(shí),為拉伸截?cái)嗖糠?。在FLAC中,以拉應(yīng)力為正,壓應(yīng)力為負(fù)。因此,H-B屈服函數(shù)可以表示為式中:σc為巖石的單軸抗壓強(qiáng)度;m及s均為本構(gòu)參數(shù)。將式(1)兩邊同時(shí)對(duì)σ3求偏導(dǎo),得H-B屈服函數(shù)的切線的斜率N?(σ3)為根據(jù)M-C屈服函數(shù),有將式(3)及式(4)中的N?用式(2)中的N?(σ3)代替,可以求出?(σ3),再將式(2)與式(3)聯(lián)立,即可求出c(σ3):令式(5)、(6)中的σ3=0,有本文以M-C本構(gòu)模型為藍(lán)本,進(jìn)行了自定義本構(gòu)模型的開發(fā)。3計(jì)算3.1塑性的m-c模型、h-b模型及自主監(jiān)測的確定計(jì)算模型見圖2,其長(x方向)及高(y方向)均為30m,巷道半徑為3m,被劃分為輻射網(wǎng)格和矩形網(wǎng)格,共計(jì)4784個(gè)單元。計(jì)算在平面應(yīng)變及小變形條件下進(jìn)行,分別采用理想彈塑性的M-C模型、H-B模型及自定義模型。計(jì)算模型的四周受到靜水壓力p1的作用。模型中的圓形巷道是預(yù)先生成的。在計(jì)算模型中,取x軸(x>0)上的節(jié)點(diǎn)及節(jié)點(diǎn)上方的一行單元作為被監(jiān)測的對(duì)象。體積模量取為1.5GPa,剪切模量取為0.5GPa。H-B模型中的m=10、s=1及σc=20MPa。由2.1節(jié)中的式(7)及(8)可以計(jì)算出M-C模型中的c=c(0)=4.08MPa及?=?(0)=45.6°,由可以求出抗拉強(qiáng)度σt=1.98MPa。自定義本構(gòu)模型中的σt的取值與M-C模型的相同,其他參數(shù)均與H-B模型的相同。3.2理論解與h-b、-b模型的結(jié)果對(duì)比每個(gè)本構(gòu)模型均采用了3個(gè)計(jì)算方案,即圍壓p1分別取為30、45、60MPa。圖3(a)~(c)及圖4(a)~(c)分別給出了各方案中應(yīng)力及位移分布的數(shù)值解與H-B理論解的對(duì)比,H-B理論解中參數(shù)的取值均與H-B模型中的相同。在本文中,環(huán)向應(yīng)力σθ、徑向應(yīng)力σr及水平位移ur均指它們的大小。H-B理論解中的σθ可以通過文獻(xiàn)中的式(25)、(26)、(31)、(32)及(36)計(jì)算得出。由圖3(a)~(c)可以發(fā)現(xiàn),自定義模型與H-B模型的結(jié)果幾乎完全相同。各種解中的σθ均出現(xiàn)了峰值。隨著圍壓的增大,σθ的峰值逐漸向圍巖內(nèi)部移動(dòng);由文獻(xiàn)中的式(35)計(jì)算得到的塑性區(qū)半徑逐漸增大(分別為3.8、4.4、5.0m,與H-B模型及自定義模型的計(jì)算結(jié)果較為接近,且均大于M-C模型的結(jié)果);自定義模型的結(jié)果與M-C模型的差別越來越明顯。圖4(a)~(c)分別給出了圍壓不同時(shí)被監(jiān)測節(jié)點(diǎn)的水平位移ur的數(shù)值解與H-B理論解(根據(jù)文獻(xiàn)中式(46)計(jì)算得出)的對(duì)比。可以發(fā)現(xiàn),各種解均是越靠近巷道表面,位移越大。自定義模型及H-B模型的計(jì)算結(jié)果幾乎完全相同,且均大于M-C模型的計(jì)算結(jié)果及H-B理論解。隨著圍壓的增大,位移越來越大,M-C模型與自定義模型的計(jì)算結(jié)果的差別越來越明顯。總之,自定義模型與H-B模型的應(yīng)力及位移的計(jì)算結(jié)果幾乎完全相同,與H-B理論解基本一致。因而可以確認(rèn)自定義模型的正確性。4圓形道路圍巖的局部化模擬4.1數(shù)值模型的建立采用的計(jì)算模型為矩形平面應(yīng)變模型,見圖5,其長及高均為10m,被劃分為40000個(gè)尺寸相同的矩形單元,模型的四周受到靜水壓力的作用,計(jì)算在小變形條件下進(jìn)行。這里采用了H-B應(yīng)變軟化模型及彈-脆-塑性的自定義本構(gòu)模型。體積模量均取為24GPa,剪切模量均取為23GPa。H-B模型中的初始及殘余本構(gòu)參數(shù)分別為m=5、mres=0.5、s=0.15、sres=0.01及σc=σcres=230MPa。自定義本構(gòu)模型中的初始及殘余抗拉強(qiáng)度分別取為σt=11MPa及σtres=1MPa,其余參數(shù)均與H-B模型相同。采用4個(gè)計(jì)算方案,方案1及方案2采用自定義模型,靜水壓力p2分別取為20MPa及30MPa;方案3及方案4采用了H-B模型,當(dāng)本構(gòu)參數(shù)達(dá)到殘余值時(shí)的塑性應(yīng)變分別為2×10-4及2×10-8,方案4中的巖石破壞后接近脆性。每個(gè)方案的數(shù)值計(jì)算分為3步,以“先加載,后挖洞”的方式進(jìn)行,具體步驟與文獻(xiàn)基本相同,巷道的半徑也與文獻(xiàn)相同,為1.67m,所不同的是本文將所有單元的應(yīng)力設(shè)置為與圍壓相等的值,以使未開挖的計(jì)算模型盡快達(dá)到靜力平衡狀態(tài)。由計(jì)算結(jié)果可以發(fā)現(xiàn),本文中各方案計(jì)算到45000時(shí)步時(shí),計(jì)算模型都達(dá)到了靜力平衡狀態(tài)。4.2模擬結(jié)果及分析圖6給出了各方案在不同時(shí)步時(shí)的剪切應(yīng)變?cè)隽康母咧祬^(qū)分布,顏色越深表示剪切應(yīng)變?cè)隽吭酱蟆？梢园l(fā)現(xiàn),當(dāng)圍壓較低(p2=20MPa)時(shí),剪切應(yīng)變?cè)隽康母咧祬^(qū)在巷道邊緣處呈點(diǎn)狀分布,這與文獻(xiàn)中低圍壓時(shí)的模擬結(jié)果類似。當(dāng)2000時(shí)步時(shí),剪切應(yīng)變?cè)隽康母咧祬^(qū)位于巷道表面的一些位置,見圖6(b);5000時(shí)步時(shí),圍巖中形成了類似文獻(xiàn)[10–11]中的4個(gè)對(duì)稱的V型坑,見圖6(c);當(dāng)10000時(shí)步時(shí),V形坑的尺寸稍型有增大,見圖6(d);20000及45000時(shí)步的計(jì)算結(jié)果與10000時(shí)步時(shí)相比無明顯變化,見圖6(e)、(f)。方案2的模擬結(jié)果與文獻(xiàn)中巖爆的現(xiàn)場觀測及實(shí)驗(yàn)結(jié)果吻合。采用H-B模型計(jì)算出的剪切應(yīng)變?cè)隽康母咧祬^(qū)的對(duì)稱性較差,所形成的V型坑的大小及形狀明顯不同,見圖6(g)、(h)。綜上所述,在計(jì)算模型的邊界條件對(duì)稱及均質(zhì)本構(gòu)參數(shù)的條件下,采用彈-脆-塑性的自定義模型且當(dāng)圍壓較大時(shí),巷道圍巖內(nèi)才會(huì)形成4個(gè)對(duì)稱的V型坑;而采用接近脆性的H-B應(yīng)變軟化模型模擬出的結(jié)果的對(duì)稱性較差。5模型計(jì)算結(jié)果不同(1)以M-C本構(gòu)模型為基礎(chǔ),利用C++實(shí)現(xiàn)了考慮拉伸截?cái)嗟姆蔷€性本構(gòu)模型的二次開發(fā)。由計(jì)算結(jié)果可以發(fā)現(xiàn),利用自定義模型計(jì)算出的應(yīng)力及位移分布與H-B