劉甫平,周 峙,羅 易
(1.安徽省交通規劃設計研究總院股份有限公司,安徽 合肥 230088;2.公路交通節能環保技術交通運輸行業研發中心,安徽 合肥 230088;3.中國地質大學(武漢) 工程學院,湖北 武漢 430074)
軟巖在漫長的沉積歷史過程中,伴生許多弱膠結的原生結構面(帶),是一種具有顯著塑性變形的非均勻介質。其峰后塑性應力應變特征一直是巖土工程關注的焦點,提出適用于軟巖的力學本構模型,不僅有利于其結構參數確定和峰后應變軟化規律的預測,也能指導軟巖地區巷道支護、井壁穩定性等工程的設計。
傳統塑性力學理論框架內,對于連續均勻介質的剪切屈服只考慮服從德魯克塑性公設[1-2],其球張量只產生球應變,偏張量只產生偏應變。而在巖土塑性力學范圍內,巖土是一種不均勻的多相顆粒體,其球張量和偏張量存在交叉影響,宏觀上表現為剪脹或剪縮。因此,嚴格來講,對于不均勻介質的材料微元不宜采用連續材料微元,而宜采用更能反映變形過程中顆粒相互影響的細觀力學模型更為合理[3](P8)。
沈珠江院士在巖土破損力學框架下提出的二元介質模型[4],將不均勻介質微元抽象為膠結元和摩擦元,前者表現為脆彈性,后者則表現彈塑性特征。當單元體受外部荷載作用時,黏聚力和摩擦力并非同時發揮作用[5],而是前者破壞后,后者才逐步發揮抗力,在此過程中會出現巖土材料的應變硬化-軟化,破壞過程也可以抽象的認為膠結元逐步轉化為摩擦元的過程,該過程能夠考慮被結構帶切割的體元之間的摩擦咬合作用對材料整體強度的貢獻,比傳統的塑性力學研究更接近于實際。
自二元介質模型提出以來,劉恩龍、陳鐵林等學者通過二元介質模型在超固結黏土、黃土、堆砌體、凍土等結構性土的研究[6-9],驗證了模型可以較好地反映結構性土的應力應變特征。而軟巖材料與結構性土具有相似的應力應變行為,但是目前在巖土破損力學理論框架內,有關軟巖在加載進入塑性狀態后的變化規律研究卻較少,多數巖樣的二元介質模型均以硬質巖為例進行驗證[10-12]。文獻[11]提出單元體在破損前后彈性模量、泊松比均相同,只是單元體的破裂強度滿足Weibull分布探究了砂巖破損函數的演化規律,并分析了參數選取對峰值強度、應變軟化的影響,模型能定性的反映砂巖的峰后軟化規律,研究主要集中于破損函數的演化規律。隨后劉恩龍在沈珠江院士的研究基礎上,將硬化參量引入破損過程中摩擦元的應力應變關系[12],基于劍橋模型給出了摩擦元的流動規則,并基于常規砂巖的三軸試驗驗證模型的合理性,研究雖然基于劍橋剪脹方程考慮了摩擦元的塑性流動規律,但是從p-q平面上看,對于顆粒類材料的剪脹線應該是非線性的,那么這就會導致模型對巖樣在塑性屈服后的剪縮-剪脹反映不足。而對于軟巖而言,其相比硬巖更易在塑性屈服階段表現剪脹(擴容)現象,因此如何提出適用于軟巖變形特征的二元介質模型十分重要。
基于此,本文在巖土破損力學的理論框架內,將含結構面的軟巖中膠結強的結構單元抽象為彈脆性元,滿足廣義虎克定律。將軟弱結構面抽象為摩擦元,在滿足Mohr-Coulomb準則基礎上,確定塑性剪應變為塑性內變量,采用非關聯流動法則獲取摩擦元的剪脹角,并引入剪脹比表征摩擦元剪脹規律,以此建議了一個能反映含結構面軟巖的應變軟化過程的增量型二元介質本構模型,最后在MATLAB中進行增量計算,并與含結構面的粉砂質泥巖進行三軸試驗結果對比,表明所建立的本構模型能合理的描述含軟弱結構面軟巖的應力應變特性。
2002年,基于結構性土的細觀破損機理,沈珠江院士提出了巖土破損力學的理論框架[13],并提出了以下假定:
1)準連續介質假定,即在二元介質中任意選取一個代表性單元(RVE),其包含多個膠結塊和結構面;
2)脆彈性破損假定,即結構脆彈性元一旦滿足破壞準則立即發生破碎。
3)共同分擔假定,宏觀應力由結構塊膠結應力和軟弱帶的摩擦應力共同分擔,滿足以下關系式
其中[σ]i、[σ]f分別為平均膠結應力(彈脆性元)與平均摩擦應力(彈塑性元),[ε]i、[ε]f分別為平均膠結應變、平均摩擦應變、ζv為體積破損率。
非均勻介質的均勻化理論根據參數關系不同,可以分為單參數均勻化理論和雙參數均勻化理論。單參數均勻化理論認為RVE中的膠結元破損過程中轉化為摩擦元時,其球應力與體應變均按體積平均。而含結構面的軟巖在進入軟化塑性狀態后,常沿某一結構面發生破壞,那么局部應變會在結構帶附近出現,因此比較合理的將球應力與體應變按體積平均,剪應力與剪應變按面積平均,定義雙參數表征均勻化,即體積破損率ζv、面積破損率ζs。
分別定義[σ]m、[σ]s分別為球應力和剪應力,分別用下標i表示膠結元、下標f表示摩擦元。
軟巖進入塑性變形階段,由于加載、卸載時的規律不一樣,其應力-應變不再具有唯一對應關系,因此,對彈脆性元和彈塑性元分別采用增量型應力應變關系表述:
由式(2)、(4)可知RVE的平均應變增量主要由彈脆性元的應變增量、彈塑性元的應變增量及彈脆性元破損轉變位彈塑性元引起的應變增量組成,后文對此進行詳細闡述。
假定膠結元破損前為理想彈脆性體,則其滿足廣義虎克定律:
Ki、Gi分別為彈脆性元的體積模量和剪切模量。
對于摩擦元,根據脆彈性破損的基本假定,當結構塊一旦滿足破損準則,那么將膠結元喪失膠結強度立即轉化為無黏性顆粒,表現為滑移特性(摩擦強度發揮階段)。
此時膠結元轉化為摩擦元過程可以用理想塑性材料的彈塑性的應力應變特性進行描述,對于常規三軸試驗,屈服函數可按Mohr-Coulomb屈服準則,其應力點在π平面應該為兩條棱線[14],因此可以表示為:
其中,Nφ為內摩擦角的函數,Nφ=(1-sinφ)/(1+sinφ)。
在顆粒受剪切過程中,會出現顆粒的翻轉、破損從而導致顯著的剪脹現象。目前常用剪脹角表述剪切變形過程的體脹過程,但是在考慮剪脹的過程中存在兩種較為常見方式,其一采用非關聯流動法,使得Ψ=0,其實是忽略了顆粒的剪脹性,其二是采用屈服函數作為塑性勢函數的關聯流動法則,使內摩擦角與剪脹角為定值,這就高估了材料的剪脹性?;诖?,本文采用非關聯法則,同時將屈服函數的內摩擦角替換為剪脹角Ψ來研究摩擦元的剪脹特性。將內摩擦角替換為剪脹角Ψ塑性勢函數可以得到:
式中Nψ為剪脹角Ψ的函數,Nψ=(1+sinψ)/(1-sinψ)。由于巖土材料一般發生剪切屈服,本文以等效塑性剪應變為內變量,采用非關聯流動法則,求出剪脹摩擦角Ψ的表達式,
(16)
當膠結元破損形成摩擦元過程中出現剪脹時,會出現臨脹特征點,定義剪脹應力比Md表示,該點對應的摩擦角Ψ0定義為臨脹摩擦角。采用剪脹比(塑性體應變與塑性剪應變的比)表示摩擦元塑性流動規律,即:
式中dg為剪脹比,υ為剪脹模型常數,與材料性質相關,η=σs/σm,其變化范圍為(0,1)。
可以用文獻[12]關于巖石的二元介質模型表征結構面軟巖摩擦元的增量型應力-應變關系:

那么有
將式(21)與(22)代入(11)、(12)式,并定義膠結元和摩擦元的剛度系數分別為D1、D2,則有

將式(23)、(24)聯立式(11)、(12),可得
將式(27)、(28)代入(25)、(26),整理可得
那么式(19)、(20)可以寫成

根據式(27)、(28)、(31)、(32)分別代入式(7)、(9)可以得到含結構面軟巖的二元介質模型如下:
含結構面軟巖的二元介質本構模型中有5組參數需要確定,對于膠結元,分別為彈性體積模量Ki、彈性剪切模量Gi;對于摩擦元,分別為彈性體積模量Kf、彈性剪切模量Gf;結構破損參數中體積破損率ζv、面積破損率ζs,局部應變系數cm、cs;臨脹特征點的剪脹摩擦角Ψ;后文將分別闡述模型參數的演化規律。
結構面軟巖漸進破損可以劃為三個階段,分別為彈性變形階段、硬化-軟化階段(膠結元轉化摩擦元)、殘余階段(膠結元完全破損摩擦元滑移),如圖1所示。

圖1 二元介質模型漸進破損過程
當試樣達到屈服應力σf之前,η=0時,該階段認為外部荷載主要由膠結應力承擔,處于線彈性變形階段;
當0<η<η*時(η*為膠結元從應變軟化轉變到殘余階段的臨界軟化參數,由殘余應變與峰值應變的比確定),表示膠結元達到初始破損門檻值(即為屈服應力)后,試樣開始表現出剪脹(縮)現象,該階段膠結元開始轉化為摩擦元,外部荷載由膠結應力及摩擦應力共同承擔,該階段在應力-應變曲線上表現為應變硬化-軟化(剪縮-剪脹);
當η>η*時,無膠結強度的摩擦元滑移,表征膠結元完全轉化為摩擦元的殘余階段。
膠結元Ki、Gi可以用不含結構面的巖樣的三軸剪切試驗應力應變曲線(σ1-σ3)-ε1的初始斜率獲取。
摩擦元的彈性體積模量和剪切模量根據峰后應力-應變曲線按下式確定:
式中μf為摩擦元的泊松比,Pa為標準大氣壓強,為101 Kpa,σc為殘余強度。
3.3.1 初始破損門檻值變化規律
結構面軟巖不同于結構性土,單純圍壓作用下不會引起膠結元的破損,因而膠結元的破損應達到其初始破損應力門檻值才會發生。隨著軸向荷載的增加,應力應變關系開始呈非線性變化,應變變化較快,在應力差的對數坐標中顯示尤為突出。在結構面軟巖峰前應力-應變曲線中存在應變陡增轉折點,即為結構面軟巖的初始破損點,對應初始破損門檻應力。
限于篇幅,本文僅列出σ3=1.0 Mpa的ln(σ1-σ3)-ε關系圖,如圖2所示,對應的不同圍壓的初始門檻應力值見表1。

圖2 初始破損應力差對數-軸向應變關系曲線

σ3(MPa)qe(MPa)σ3(MPa)qe(MPa)0.510.51.514.21.011.62.016.4
3.3.2 結構破損參數變化規律
巖石內RVE單元體的膠結元的局部應力達到破損門檻值時即發生破損,膠結元逐步轉化為摩擦元,這個過程中體積發生膨脹,同時會出現破損帶,因此,破損率范圍值從0逐漸變化至1[15],如圖3所示。

圖3 破損率曲線
對于巖石而言,假設破損率滿足Weibull分布[14],按如下關系式表示:

m值的選取與結構面軟巖的材料性質和結構面膠結程度有關。m值愈大,破損率變化愈快,峰值強度愈低,巖石更易發生破損,本文建議膠結強度高的巖石模型m值取小值,對于弱膠結的巖石取大值。局部應變系數cm、cs也是無法直接量測,可通過先假設,后試驗驗證的方式間接確定。cm、cs為應力的函數,那么可以構造下列表達式,
其中αv、αs分別為模型參數,通過選取擬合可以求得。此處通過將式(35)、(36)、(37)、(38)代入(33)、(34)中通過與試驗曲線對比可以確定m、αv、αs參數值。
根據三軸試驗εv-ε1曲線進行整理,獲取塑性體積應變與軸向塑性主應變與等效塑性剪應變的曲線最優化擬合關系,并分別對塑性體積應變、軸向塑性主應變微分,給定一個等效塑性剪應變的增量步,可以求得不同等效塑性剪應變下的剪脹角Ψ,以此獲得膠結元轉化摩擦元過程中的剪脹比。
以1.0 MPa圍壓下的塑性體積應變與軸向塑性主應變與等效塑性剪應變的曲線最優化擬合關系為例進行闡述,如圖4所示。塑性體應變增量和塑性軸向應變增量可以用等效剪應變增量表示為:
代入式(16)可知

圖4 塑性軸向應變、體積應變與等效塑性剪應變關系圖
屈服應力對應的等效塑性剪應變可以得到臨界剪脹角,代入公式(17)、(18)及本構關系(33)、(34)確定dg并與試驗曲線進行對比,選取最優模型參數υ。

表2 本文計算采用臨脹角
根據式(18)分別計算應力比與剪脹比的關系,同時與Rowe剪脹方程和劍橋剪脹模型計算結果對比,以1.0 MPa下的試驗數據為例進行對比分析,不同剪脹參數規律如圖5所示。

圖5 不同剪脹參數規律(σc=1.0MPa)
摩擦元在破損過程中應力比與剪脹比的關系應該是非線性的[16],而劍橋剪脹方程為線性關系,并不能很好的反映剪脹性。而Rowe剪脹方程所表述的剪脹關系則明顯低估了摩擦元的剪脹性。采用本文的剪脹比模型與試驗結果擬合較好。同理可以求得其他不同圍壓下的剪脹比和最優的剪脹模型參數υ。
試驗巖樣為弱風化含青色泥質結構面的紫紅色粉砂質泥巖,天然密度2.63 g/cm3。試樣尺寸為?50 mm×100 mm的圓柱形巖樣,常規三軸試驗采用中國地質大學巖石力學實驗室的INSTRO微機控制巖石伺服三軸壓力試驗機,采用位移控制加載方式,加載速率0.002 mm/s,分別對試樣進行圍壓為0.5 MPa、1 MPa、1.5 MPa、2.0 MPa的三軸剪切試驗,破壞形式如圖6所示,并與理論計算結果進行對比。
計算采用的模型參數Ki=4.197 GPa、Gi=1.800 GPa,μf=0.18, 改變圍壓的情況下可求得不同圍壓下的Kf、Gf,mv=2,ms=2,αv=0.001,αs=200,圍壓為0.5 MPa、1 MPa、1.5 MPa、2.0 MPa的剪脹模型參數υ分別為1.93、1.51、1.34、0.92,η=0.7,在MATLAB中編程實現模型增量計算,首先設定一個加載步,給定一個微小的dεv、dεs、dζv、dζs,在摩爾-庫倫準則下判定是否屈服,滿足屈服條件后記錄當前應力狀態,并根據非關聯流動法則計算等效塑性剪應變增量和剪脹角,由于破損率ζ本身是局部應力函數,局部應變系數cm、cs也是局部應力的函數,如若考慮每個增量步中兩者相互轉化,求解非常復雜,本文假定每個增量步中cm、cs不變,按上一步ζv、ζs增量計算ζ,以此求得每個單元彈性應變增量和塑性應變增量,如此不斷給定加載步,即可求出不同圍壓下的偏應力-軸向應變曲線和體積應變-軸向應變曲線,對比計算曲線與試驗曲線后迭代給出ζv、ζs、cm、cs,不斷優化計算結果,數值計算結果與試驗結果對比如圖7所示。建議的模型與試驗曲線擬合較好,可較好地反映含結構面軟巖屈服后的應力-應變關系,在低圍壓下,隨著應變的增加,試樣表現較強的應變軟化現象(剪脹),隨著圍壓的增大,剪脹逐漸變弱。

圖7 模型計算結果與試驗結果對比
1)在巖土破損力學理論框架內,將含結構面軟巖抽象為膠結元和摩擦元的非均勻介質,認為兩者在破損規律下相互轉化,相互影響,能反映巖石單元體體應力與剪應變、偏應力與體應變的交叉影響關系。
2)基于摩爾-庫倫準則和非關聯流動法則,并以塑性剪應變為內變量,將剪脹比引入二元介質破損模型中考慮摩擦元破損過程中的剪脹性,從而更能反映單元體漸進破損過程中的應變硬化-軟化特征。
3)結構面軟巖的破損不同于黃土、松砂、結構性土,單純圍壓作用下不會引起膠結元的破損,因而膠結元的破損應達到其初始破損應力門檻值才會發生。
4)增量計算結果與常規三軸試驗結果進行對比分析,發現建議的模型可以較好地反映含結構面軟巖的應力-應變特征。提出的模型是在試驗基礎上建立的,為類似性質巖石提供借鑒。
本文僅對低圍壓水平下的軟巖進行了試驗和模型驗證,以后需加強模型在高圍壓應力水平下的剪脹特性驗證。