李彥奇,黃 達,孟秋杰
(河北工業(yè)大學土木與交通學院,天津 300400)
在我國西部山區(qū),天然或人工開挖形成了大量的反傾巖質邊坡[1?2]。早期工程地質工作者認為反傾層狀邊坡較為穩(wěn)定,不易形成貫通面引發(fā)滑坡[3]。然而根據(jù)對我國近100例發(fā)生較大變形或滑坡的斜坡進行的統(tǒng)計,其中反傾向邊坡占33%[4],出現(xiàn)的頻度和危害較大[5]。塊狀-彎曲傾倒是反傾層狀巖質邊坡的主要傾倒形式之一,主要特性表現(xiàn)為坡頂后的彎曲傾倒以及坡頂前的彎曲傾倒的結合[6]。這種破壞模式在軟硬相間的層狀巖體中非常普遍,多發(fā)于砂巖泥板巖互層、燧石巖頁巖互層、薄層狀石灰?guī)r中[7]。考慮到軟巖、硬巖有不同的地質及工程特性,離心機模型試驗和數(shù)值模擬成為研究此類坡體變形的主要研究手段。
關于離心機試驗,Adhikary等[8]開展了7組不同結構的均質反傾邊坡離心試驗,得出了巖層傾角、巖石抗拉強度等對破裂面形狀的影響;鄭達等[9]在開挖條件下對均質層狀反傾邊坡進行離心試驗,得出臨空條件是深層傾倒變形破壞的關鍵致災因子。上述成果均建立在相鄰巖層巖性相同的假定條件之下,未考慮軟硬相間互層狀結構。
在數(shù)值模擬方面,黃波林等[10]使用UDEC軟件對廖家坪危巖體進行了模擬,結果表明坡體傾倒的方向受軟弱層控制;黃潤秋等[11]以皖南某高速公路反傾層狀邊坡為例,使用離散元法分析了其傾倒破壞機理并揭示了邊坡的變形分區(qū)現(xiàn)象。以上研究總結了坡體的破壞特征,但對于小變形狀態(tài)下的細節(jié)兼顧不足;對于工程實例模型往往需要大量的標定來校核其微觀參數(shù),在同一參數(shù)不同邊界條件下的繼承性也略顯不足。
白潔等[12]以苗尾水電站為背景構建反傾層狀巖質邊坡有限元模型,探討了此類坡體的變形特點及穩(wěn)定性影響因素;姚文敏等[13]采用FLAC3D強度折減法,在巖層傾角、巖層與邊坡走向夾角變化時,探討了三維軟硬互層邊坡的穩(wěn)定性情況。以上研究由于方法限制無法進一步得到坡體的最終破壞情況。王宵等[14]以某梯級電站廠后邊坡為例,分別構建二維離散元和三維有限元模型,對“似層狀”巖質邊坡具體演化過程進行了詳盡的描述,但兩模型的銜接方面并不一致,與實際條件不符。
F-DEM(有限元-離散元)方法[15]可以同時考慮巖體結構面和巖塊,融合了有限元和離散元各自的優(yōu)勢,可以考慮巖層沒有完全脫離母巖的實際情況。借助此方法,陳小婷等[16]以巫山縣長江左岸危巖體為例,證明了將連續(xù)-離散模型用于邊坡破壞研究具有可行性;劉郴玲等[17]以紅石邊坡工程為實例,建立FDEM模型并獲得邊坡失穩(wěn)破壞的漸進破壞全過程。
本文首先對巖層傾角為優(yōu)勢傾倒角[18](60°)的軟硬互層反傾邊坡進行試驗研究,提出一種可以考慮拉剪和壓剪的節(jié)理單元強度準則,并與軟硬互層離心機模型試驗相結合,分析了軟硬互層巖質邊坡中巖層傾角、軟/硬巖層厚比對模型破壞機理的影響。
1.1.1 相似比設計
試驗設計各模型/原型的比值如下:重力相似比(Ca)為90∶1、幾何相似比(Cl)為1∶90、彈性模量相似比(CE)為1∶1、材料密度相似比(Cρ)為1∶1,黏聚力相似比(Cc)為1∶1,內摩擦角相似比(Cφ)為1∶1,應變相似比(Cε)為1∶1。最終預堆砌模型邊坡的尺寸為70.27 cm(長)× 55 cm(寬)× 56.72 cm(高)。
1.1.2 相似材料設計
石膏水泥混合材料是目前模擬巖石行為最常用的相似材料,可用于控制材料的力學性質。為避免硬化過程中受環(huán)境濕度與溫度的影響,加入硼砂水(m硼砂∶m水=1∶49)。經正交配比試驗,再通過單軸壓縮試驗和三軸壓縮試驗,得到原型和相似材料的物理參數(shù)(表1),最終材料配比見表2。通過100 mm×50 mm試樣單軸壓縮試驗,得出軟巖相似材料抗壓強度為0.3 MPa,硬巖為10 MPa。
表1 原型材料與模型材料的物理力學參數(shù)Table 1 Physico-mechanical parameters of the prototype andmodel materials
表2 軟硬巖相似材料最終配比Table 2 Ratio of similar materials for hard and soft rocks
1.2.1 模型信息
本文軟硬互層邊坡離心物理模型試驗在地質災害防治與地質環(huán)境保護國家重點試驗室(成都理工大學)TLJ-500型巖土離心機上完成。
邊坡模型巖層傾角為60°(優(yōu)勢角[18]),由單塊板(巖層)拼裝而成,試塊分別有60 cm×10 cm×1 cm和40 cm×10 cm×1 cm兩種類型。如圖1所示,涂黑巖層為硬巖層,未涂黑巖層為軟巖層,巖層厚度均為1 cm。為避免邊界效應和降低側壁摩阻力的影響,在模型箱與模型接觸的位置粘貼光滑的塑料薄膜。
圖1 模型及監(jiān)測點信息Fig.1 Model and monitoring point information
1.2.2 巖體節(jié)理及層間處理方法
巖層搭接方式為錯縫搭接(圖2),且相鄰兩塊板使用軟石膏砌筑加強黏結,這樣搭接可以大限度避免垂直于坡面方向節(jié)理的影響[9]。
圖2 模型搭接方式及節(jié)理情況Fig.2 Overlap mode and joint of the model
在層狀坡體的制作過程中,相鄰層理面宜界面均勻、層厚均勻,所以在澆筑拼裝時須注意控制時間,每澆筑完成一層,養(yǎng)護1 h再澆筑下一層。層間黏結力宜較小,因此層間材料選用抗壓強度0.1 MPa、抗拉強度10 kPa的軟石膏。
1.2.3 加載方式
由離心機基本原理可知:
式中:am—離心加速度/(m·s?2);
N—離心機轉速;
g—重力加速度/(m·s?2)。
離心機模型的尺寸將由離心加速度放大至N倍,由此使較大尺寸的邊坡變形破壞過程的模擬試驗成為可能。
模型所受離心加速度首先從0g開始逐漸加載至30g,并在30g保持離心機勻速旋轉5 min;然后加速度從30g逐漸加載至50g,并在50g保持勻速旋轉5 min,最后模型所受加速度從50g逐漸加載至90g,并在90g保持勻速旋轉5 min。若模型不破壞則繼續(xù)增加離心加速度直到模型破壞。
1.2.4 監(jiān)測點布置
本研究采用導軌支撐位移計量測坡體表面位移(圖1)。在坡表設置兩個監(jiān)測點:監(jiān)測點A和監(jiān)測點B,并在層間埋置壓力傳感器監(jiān)測層間法向應力。
F-DEM模型由實體單元和零厚度的節(jié)理單元兩種單元構成(圖3)。其中,零厚度的節(jié)理單元夾在實體單元中間。當節(jié)理單元受力過大時,程序在當前增量步完成后將節(jié)理單元刪除,并在刪除位置生成對應的裂縫。節(jié)理單元刪除后的兩個實體單元間依舊保持接觸關系,可以碰撞和摩擦,法向接觸為硬接觸,切向為庫倫摩擦,其中硬接觸如圖4所示。
圖3 F-DEM模型示意圖Fig.3 Schematic diagram of the F-DEM model
圖4 法向硬接觸Fig.4 Normal hard contact
對于三維單元,有
式中:E—彈性模量/MPa;
μ—泊松比;
σx、σy、σz—x、y、z方向應力;
εx、εy、εz—x、y、z方向應變;
γzx、γxy、γyz—x、y、z方向切應變。
本文所使用的節(jié)理單元假設剪應力為0(在局部可以將沿著節(jié)理方向的兩個主應力看作剪應力)[19],所以式(1)可簡化為:
由于剪應力都為0,所以在節(jié)理單元中只存在3個主應力:σ1、σ2和σ3。在局部坐標系(圖5)中:
圖5 局部坐標系Fig.5 Local coordinate system
Hoek-Brown準則和Mohr-Coulomb準則是常用的用來描述巖石抗剪強度隨法向應力變化規(guī)律的強度準則。Hoek-Brown準則是根據(jù)三軸壓縮試驗數(shù)據(jù)得到的經驗公式:
式中:σ1、σ3—最大和最小主應力/MPa;
σ0—巖石的單軸抗壓強度/MPa;
m—與巖性相關的常數(shù);
S—巖石完整性的參數(shù),完整巖石時S= 1。
另外,巖石的單軸抗拉強度可用式(5)預測:
Mohr-Coulomb準則用式(6)表示:
Hoek-Brown準則比Mohr-Coulomb準則能夠更好地描述巖石在拉剪應力狀態(tài)下的強度特征,而Mohr-Coulomb準則比Hoek-Brown準則能夠更好地描述巖石在壓剪應力狀態(tài)下的強度特征[20]。采用分段函數(shù)的方式來對砂巖拉剪強度和壓剪強度進行預測:
其中式(4)中的第二項又可寫為:
因此將Hoek-Brown準則與Mohr-Coulomb準則聯(lián)合使用,既能較為準確地描述巖石在拉剪應力狀態(tài)下的強度特征,又能較為合理地揭示巖石在壓剪應力狀態(tài)下的強度特征,可以更好地描述(模擬)巖石的斷裂行為。VUMAT子程序的計算目的是聯(lián)合兩種強度準則以及生成新的裂縫。每一次迭代過程都將遍歷模型中各個單元。若判斷為受壓狀態(tài),則進一步判斷其當前強度是否超過其Mohr-Coulomb準則下抗剪強度,若大于抗剪強度則將損傷因子D更新為1,清零其剛度矩陣。此時單元失效。若判斷為受壓狀態(tài),則進一步判斷其當前強度是否超過其Hoek-Brown強度準則下的抗剪強度,若大于抗剪強度則將損傷因子D更新為1,清零其剛度矩陣。此時單元失效,裂縫產生。此時實體單元之間的接觸關系依然存在,與母體分離后的巖塊依然可以相互摩擦、碰撞。其中σn>0時模型處于壓剪狀態(tài),σn≤0時模型處于拉剪狀態(tài)。
2.3.1 F-DEM中的時間及零厚度單元插入
使用ABAQUS中Dynamic、Explicit分析步進行計算,此算法中的時間即為運動方程中的時間,是真實時間。編寫Python程序對有限元模型進行前處理,步驟如下:
(1)將原始模型中的節(jié)點分別賦予兩個不同的序號,其中一個節(jié)點為虛擬節(jié)點;
(2)將新產生的節(jié)點編號以及節(jié)理單元節(jié)點的編號添加到模型文件中;
(3)添加節(jié)理單元的數(shù)據(jù)信息,即在源文件的關鍵字中聲明單元類型的更改;
(4)為這些新加入的單元設置材料和截面屬性并再次導入ABAQUS中進行邊界條件等的設定。
2.3.2 Hoek-Brown—Mohr-Coulomb聯(lián)合強度準則和本構模型參數(shù)
圖6為軟硬互層反傾層狀邊坡數(shù)值模型圖,模型按單元不同分為兩部分:實體單元部分和零厚度的節(jié)理單元部分。
圖6 軟硬互層反傾層狀邊坡數(shù)值模擬Fig.6 Numerical model of anti-dumping rock slope interbeded by hard and soft layers
根據(jù)坡體巖性和結構面,將節(jié)理單元分為層間節(jié)理單元、硬巖節(jié)理單元及軟巖節(jié)理單元。因VUMAT中的參數(shù)都為實際物理量,所以可以依照坡體離心機試驗模型中的材料參數(shù)設置數(shù)值模型中的參數(shù),且模型建立兩者一致(其中Hoek-Brown模型中的參數(shù)經過換算得到[21])。對于實體單元,則直接使用單軸壓縮試驗所測得的材料參數(shù)。對于層間節(jié)理單元的參數(shù)通過與試驗現(xiàn)象的對比校核來確定。通過以上材料參數(shù)的確定方法來確保物理模型與數(shù)值模型兩者之間材料的相似。數(shù)值模型各部分材料具體材料參數(shù)見表3、表4。
表3 節(jié)理單元參數(shù)Table 3 Joint element parameters
表4 實體單元參數(shù)Table 4 Solid element parameters
2.3.3 計算工況及條件設定
將10 120個零厚度節(jié)理單元嵌入到6 828個實體單元之間。為確保物理模型與數(shù)值模型兩者之間邊界條件一致,離心機模型(圖1)按1∶1比例構建相應F-DEM模型,模型尺寸為700 mm×570 mm,約束所有節(jié)點厚度方向的位移。為保證數(shù)值計算在加載條件上與離心機試驗一致,數(shù)值模型按1.2.3節(jié)中所述加載條件進行增重加載。
監(jiān)測位移曲線及數(shù)值計算得出的位移曲線對比情況見圖7(監(jiān)測點布置見圖1),二者變化趨勢高度接近,初步說明了所建立數(shù)值模型的正確性。
圖7 監(jiān)測點位移對比Fig.7 Monitoring point displacement
利用壓力傳感器對模型的層間法向應力進行監(jiān)測,得到監(jiān)測結果(圖8)??芍?,層間法向壓力監(jiān)測點b1、b2、b4的壓力較早上升且伴有波動,加載重力進一步增大至50g時,邊坡b1、b2監(jiān)測點的壓力值進一步增加。當加載接近90g時,隨著深層彎折帶的貫通,壓力監(jiān)測點b1、b2的值到達最大。
圖8 層間法向壓力監(jiān)測曲線Fig.8 Monitoring of normal pressure between layers
圖9(a)為離心加速度10g時坡體的變形情況。根據(jù)圖8的位移曲線,加載至20g之前,邊坡處于穩(wěn)定階段,坡內無可見變形;此時監(jiān)測點b1、b2所測得的應力也未到達第一次峰值。超過20g后,位移量出現(xiàn)逐漸增大,持續(xù)至30g,后緣出現(xiàn)拉張裂縫,坡頂及坡表出現(xiàn)輕微的相互錯動,此時監(jiān)測點b1、b2、b4也出現(xiàn)應力波動。加載至30g時,位移曲線出現(xiàn)小幅躍升(圖7),此時坡體已處于破壞臨界狀態(tài)。圖9(b)為此時變形情況,坡腳局部巖體折斷、裂縫明顯,且坡頂拉張裂縫數(shù)量迅速增多。50g時,位移急劇增大,邊坡變形如圖9(c)。坡體內一級破裂面已接近貫通。與此同時b1、b2監(jiān)測點的應力迅速增大,因為此時破裂面以上巖體發(fā)生傾倒變形,坡腳巖體被擠出,形成脫空,加劇了上部巖體的變形。此時被擠出巖體的末端已開始形成更陡的二級破裂面。90g時,監(jiān)測點A、B的位移以及監(jiān)測點b1、b2的應力達到峰值,并暫時維持穩(wěn)定,此時邊坡變形見圖9(d),第二破裂面已向上貫通至坡頂,并直接導致上部巖體傾倒變形,產生三級破裂面。此時幾乎所有層間法向壓力計的測量值都有所增加,說明層間空隙被層間法向壓力擠實,裂縫的發(fā)展接近穩(wěn)定狀態(tài),變形破壞達到最終形態(tài)。
圖9 數(shù)值模擬與試驗對比Fig.9 Comparison between numerical simulation and experiment
通過以上分析可知,離心機試驗模型的變形破壞情況與F-DEM數(shù)值計算模型基本一致,驗證了數(shù)值模型及所提出強度準則的正確性,且層間應力與于此類邊坡變形破壞的關聯(lián)性較大。為能夠靈活模擬不同層間強度因素的影響,建議將層間節(jié)理單元的本構修改為“牽引-分離”準則。
巖層傾角以及軟/硬巖層厚度比是反傾軟硬互層邊坡的主要幾何因素[11],坡頂位移以及破裂面形態(tài)可以表征傾倒破壞特征。固定軟/硬巖層厚比為1∶1(10 mm∶10 mm),選取40°、60°、80°巖層傾角,確定方案1用以對比分析不用巖層傾角對邊坡傾倒破壞特征的影響;固定巖層傾角為60°選取軟/硬厚度比1.5∶1(15 mm∶10 mm)、2∶1(20 mm∶10 mm)、1∶1(10 mm∶10 mm)、1∶1.5(10 mm∶15 mm)以及1∶2(10 mm∶20 mm)確定方案2,用以對比巖層軟/硬巖層厚度比不同對邊坡傾倒破壞特征的影響。
為保證數(shù)值試驗方案的關聯(lián)性與合理性,所建不同的傾角和不同軟/硬巖層厚度比的模型均按照2.3.3小節(jié)中的邊界條件進行設置。
4.2.1 巖層傾角對邊坡最終破壞形態(tài)的影響
巖層傾角40°時,邊坡達到最終破壞時的離心加速度為100g,坡體后邊局部小范圍張開,一級破裂面與水平面的夾角(銳角)近似等于60°,且坡腳處有較小范圍的巖塊剪出,如圖10(a)所示。巖層傾角60°時,坡體共產生兩條貫通的破裂面,坡體沿二級破裂面滑出,如圖10(b)所示。相對于巖層傾角40°時邊坡的最終破壞形態(tài),巖層傾角60°時坡體后緣的張開范圍以及坡腳剪出范圍更大,且一級破裂面的分布位置更深。此時邊坡達到最終破壞時的離心加速度為90g,相對于巖層傾角40°、60°時,巖層傾角80°時,邊坡達到最終破壞時的離心加速度為105g,邊坡更傾向于整體破壞,后緣張開范圍較大,如圖10(c)所示。這說明隨著巖層傾角的增大,破裂面由一條逐漸向多條發(fā)展;且?guī)r層傾角越大,邊坡破壞的整體性越強。
圖10 不同傾角下邊坡最終破壞形態(tài)Fig.10 Final failure pattern of slope under different inclination angles
由以上分析可知,巖層傾角對第一破裂面的影響較為顯著(圖11),巖層傾角40°時,邊坡一級破裂面較為陡峭且呈現(xiàn)出連續(xù)的具有臺階狀的滑移面。隨著巖層傾角的進一步增大,破裂面向坡體更深處發(fā)展,“臺階狀”形式減弱。隨著巖層傾角的增大,邊坡一級破裂面逐漸向坡體深處發(fā)展,當巖層傾角80°時,邊坡一級破裂面發(fā)育位置較深,且后緣拉裂縫較大,邊坡發(fā)生整體傾倒,說明此時邊坡發(fā)生以“彎曲(傾倒)-拉裂”為主的破壞[3]。
圖11 不同傾角下邊坡一級破壞面Fig.11 First-order fracture surface with different dip angles
4.2.2 巖層傾角對坡體位移特征的影響
圖12為監(jiān)測點A所記錄的位移-加速度曲線。在滑坡啟動階段,重力加速度較小,不同傾角對監(jiān)測點A的豎向位移影響不大。隨著重力加速度的增大,傾角越大,坡體初始滑移所需重力加速度越大。巖層傾角越大巖層垂向的重力分量越小,所以傾角越大初始滑移需要的重力加速度越大,進一步增重后進入最終階段形成滑坡。巖層傾角越大最終產生位移越大,坡體破壞程度越劇烈。在啟動階段達到初始滑移所需的重力加速度越大。重力做功下,坡體積累的總的彈性應變能越大,突然釋放時所產生的動能也就越大。所以,巖層傾角越大最終階段產生的位移越大。
圖12 不同巖層傾角下豎直位移Fig.12 Vertical displacement with different dip angles
4.3.1 軟/硬巖層厚度比對邊坡最終破壞特征的影響
為探究軟/硬巖層厚度比對邊坡破壞特征的影響,將傾角60°,離心加速度90g時,不同軟/硬巖層厚度比的邊坡破裂面繪制于圖13。如圖13(a)—(d)所示,坡體后緣存在兩處折斷區(qū)域。隨著軟/硬巖層厚度比由1∶1.5向2∶1轉變,兩個折斷區(qū)之間距離越來越近,最終合為一個折斷區(qū)。隨著硬/軟巖層厚之比增大,坡體后緣折斷區(qū)深度逐漸增加。坡體在增重過程中的附加應力主要由硬巖承擔。硬巖巖層所占比例越大,坡體增重過程中所積累彈性應變能越大,坡體破壞釋放的能量就越大,坡體滑動的整體性就更強。坡體滑動的整體性也表現(xiàn)在破裂面的分布深度上,對比圖13(a)—(e)可知,二級破裂面的分布深度隨著軟/硬巖層厚比的提高逐漸減小。
圖13 不同軟/硬巖厚度比邊坡最終破壞形態(tài)Fig.13 Final failure pattern of slope with different ratio of soft/hard thickness
硬巖巖層厚度越大,能承受的彎矩越大,則二級破裂面分布越深。破裂面的形態(tài)也發(fā)生了巨大變化:隨著軟/硬巖層厚度比的增加,二級破裂面的形態(tài)逐漸由粗糙的“鋸齒狀”向平滑的“圓弧狀”轉變。這是因為軟巖層厚所占比例越大,越接近土質滑坡的大弧度破裂面。坡體中的一級破裂面也有類似特征。
4.3.2 軟/硬巖層厚度比對邊坡位移特征的影響
如圖14所示,隨著加速度的增加,坡頂?shù)呢Q向位移呈階梯式增加。軟/硬巖層厚度比為2∶1時,坡頂豎向位移最大(約250 mm)。隨著硬巖層厚所占比例的提高,坡頂豎向位移逐漸減?。卉?硬巖層厚度比為1∶2時,坡頂豎向位移最小(約170 mm),這是因為硬/軟巖層厚度比值越大,坡體的整體彈性模量越大。
圖14 不同軟/硬巖層厚度比下坡頂豎直位移Fig.14 Vertical displacement of downhill top with different ratios of soft/hard thicknesses
(1)經離心試驗驗證,在F-DEM模型節(jié)理單元中使用Hoek-Brow與Mohr-Coulomb聯(lián)合強度準則VUMAT子程序可以較好地模擬軟硬互層反傾層狀巖質邊坡破裂面的大致開裂情況。
(2)軟/硬巖層厚度比為1∶1,巖層傾角60°的坡體破裂面大致有三個。其破壞過程為:坡腳部位開始出現(xiàn)彎曲變形,邊坡后緣出現(xiàn)拉張裂縫,邊坡整體向臨空面彎曲傾倒。從坡體內部至坡表依次形成三個破裂面。
(3)隨著巖層傾角的增大,邊坡一級破裂面逐漸向坡體深處發(fā)展,且坡體由“彎曲-拉裂”為主的破壞模式向“滑移-壓致拉裂”為主的破壞模式轉變。巖層傾角越大最終產生的位移越大,坡體最終的受破壞程度也就越大。
(4)軟巖所占比例越大坡頂位移越大;硬巖所占比例越大坡頂位移越小。且不同比例的軟/硬巖層厚度對坡體滑動的整體性以及破裂面形態(tài)有一定影響:硬巖所占比例越大,坡體滑動的整體性越強;軟巖所占比例越大,坡體的破裂面越接近土質滑坡的“圓弧狀”。