砂巖-泥巖沉積互層面剪切特性數值模擬研究

童第科， 黃詩淵， 簡富獻， 袁 鵬, 劉叢陽

(1. 招商局重慶交通科研設計院有限公司 山區道路工程與防災減災技術國家地方聯合工程實驗室，重慶 400067；2. 重慶交通大學 河海學院，重慶 400074； 3. 江西省交通運輸科學研究院有限公司，江西 南昌 330200)

0 引 言

重慶境內廣泛分布三疊系、侏羅紀等時期所形成砂巖泥巖互層結構巖體[1]。李安洪等[2]分析了渝懷線多處順層巖質斜坡地質成巖及巖性組成特點，查得沉積砂巖泥巖巖層厚度從幾十厘米至幾米不等，地質沉積層面清晰。巖體是由巖塊及結構面組成的地質體，其強度必然受到巖塊和結構面強度及其組合方式(巖體結構)控制[3]；因而巖層結構面強度成為影響順層巖質邊坡變形及破壞的重要因素。在砂巖泥巖互層的典型軟硬巖互層地質結構中，結構面力學性質對巖質斜坡變形影響尤為突出。

在力學領域，界面力學性質是研究復合材料性能的重點[4]；受材料、接觸面粗粗糙度、黏結強度、幾何邊界等諸多影響條件影響而成為研究難點[5-6]。巖土工程中，巖體結構面力學因受工程地質及水文地質影響而復雜多變。針對其特殊性，學界當前主要研究方法有：① 原狀結構面室內及現場試驗分析[7]；② 相似材料室內模擬[8-9]；③ 結構面統計對比分析[10-11]；④ 界面數值仿真計算[12]。

受工程取樣及室內試驗條件限制，筆者參閱第①、②種方式研究成果，基于界面數值仿真，對地質沉積過程中的互層結構順層巖質邊坡中砂巖與泥巖界面剪切特性進行了模擬分析。

1 層面模擬方法

數值分析中，軟巖與硬巖相接層面(即巖體結構層面)通常以無厚度接觸面單元或薄層單元等接觸形式簡化；前者以單元法向及切向剛度作為單元主要參數，后者以彈性模量E、剪切模量G及泊松比ν作為單元主要參數。

1.1 接觸面模擬

在ABAQUS中允許用戶采用子程序形式自定義非線性接觸模型[13]，材料間剪應力與相對剪切位移之間符合雙曲線關系[14]。曹云鋼[15]指出：砂巖與泥巖交互型層面在不同法向應力下剪切過程中應力應變曲線呈雙曲線幾何特征。該方法本構模型如式(1)～(3)：

(1)

(2)

(3)

式中：△τ為應力增量；△γ為應變增量；K1、K2、Rf、n分別為非線性指標；δ為接觸面界面摩擦角；γw為水容重；Pa為大氣壓力。

1.2 薄層單元模擬

C.S.DESAI等[16]為克服Gooman接觸模型中法向剛度過大時存在過量嵌入缺點而引入了薄層單元理論；顧沖時等[17]在探究碾壓混凝土壩層面特性時同樣引入了薄層單元形式，并給出了其參數計算原理。薄層單元厚度通常取計算模型特征單元厚度0.01～0.10倍[18]。由于薄層單元厚度D比單元特征尺寸e小得多，故可忽略部分應變變量，即此方法本構模型如式(4)～(7)：

(4)

(5)

(6)

(7)

2 計算模型及參數

2.1 計算模型

根據文獻[15]中砂巖泥巖試樣尺寸及沉積面位置，筆者建立了二維平面數值分析幾何模型及有限元網格，如圖1。砂巖泥巖結構面剪切特性數值模擬中，限制上部泥巖巖樣x方向位移、下部砂巖巖樣y方向位移，在泥巖樣上表面施加法向均布荷載；隨后在砂巖試樣左側表面基于斜坡加載及應力控制方式以2.0 kPa·s-1速率施加剪切荷載。剪切過程中界面處節點應力應變曲線變化趨勢相似，為消除幾何邊界影響，設定模型中心(即網格坐標原點)為剪切應力與位移數據的提取源點。

圖1 數值分析幾何模型及有限元網格Fig. 1 Geometric model and finite element mesh fornumerical analysis

2.2 計算參數

對比土體與鋼板、混凝土、巖石等剛度較大材料接觸界面力學特性，并參閱文獻[19]，筆者分析了巖體結構面兩側巖壁的巖性、剛度、粗糙度及接觸界面黏結強度等地質參數，結合ABAQUS手冊，取定初始計算參數值為：砂巖彈性模量E=8.0 GPa，泊松比ν=0.28；泥巖彈性模量E=4.0 GPa，泊松比ν=0.32。Goodman接觸面模擬，其層面參數K1=K2=30 000(該模型為二維平面，此處取兩者相等)，Rf=0.88；n=0.11；δ=47°；γw=10 kN·m3；Pa=100 kPa；薄層單元的模擬，其薄層單元厚度D=0.25 mm，E=0.04 GPa，ν=0.40，G=0.04 GPa。采用CPE4平面應變單元進行計算，不設置砂巖及泥巖塑性階段參數和硬化參數。

3 層面剪切特性模擬效果分析

3.1 接觸面模擬

以法向荷載為0.55 MPa為例，截取部分工況計算過程中非線性階段(剪切初期)和峰后階段(殘余應力階段)云圖，如圖2。剪切初期，模型下部砂巖左邊界施加均布荷載，則上部泥巖左下端受拉應力、右下端受壓應力，故界面兩端存在應力集中現象；殘余應力階段，界面僅存在殘余接觸應力，故試樣應力云圖近似呈均勻分布。

圖2 接觸面應力云圖Fig. 2 Contacting surface stress nephogram

提取模擬數據，在同一平面內繪制不同法向應力下的砂巖泥巖結構面剪切應力隨剪切位移變化曲線，并與文獻[15]中兩組砂巖泥巖結構面室內直剪試驗數據進行對比分析，如圖3。對比分析兩組6種工況下砂巖泥巖結構面剪切應力應變的試驗曲線及數值模擬曲線變化趨勢和吻合度程度可知：砂巖泥巖層面(即結構面)剪切應力應變曲線在非線性階段都較為符合雙曲線，且均呈現應變軟化現象，存在明顯殘余強度。

3.2 薄層單元模擬

與接觸面模擬分析方法相似，以法向荷載為0.55 MPa為例，剪切過程部分云圖見圖4。采用薄層單元模擬界面時，非線性階段存在應力不對稱，這既與在砂巖左下邊界施加荷載相關，也與界面薄層單元結構彈塑性力學特性相關；剪切峰值后其應力云圖仍為泥巖左下角區域受拉應力、右下角區域受壓應力分布規律，拉壓應力在數值上相近。

同樣將在同一平面內繪制與文獻[15]結構面直剪試驗成果與薄層單元模擬結果，見圖5。對比兩模擬結果可知：薄層單元中剪切應力應變曲線在峰值前期較為迅速達到其剪切極限強度，隨后呈近乎水平直線狀。這種現象可能是因為單元屬實體單元類型所致，初期階段相同變形時其應力增長較為迅速。

由上述分析可知：薄層單元模擬層面剪切力學特性時存在初期應力劇增現象，模擬數據與試驗值存在一定程度離散性；但從整體模擬效果而言，從曲線變化趨勢表征層面剪切力學特性角度，采用薄層單元結構簡化砂巖與泥巖接觸界面并模擬其力學特性途徑仍可作為一種界面數值分析方法。

圖3 接觸面模擬剪切應力與應變Fig. 3 Shear stress-strain curves of contacting surface simulation

圖4 薄層單元模擬剪切應力云圖Fig. 4 Shear stress nephogram of the thin layer element simulation

圖5 薄層單元模擬應力與應變Fig. 5 Stress-strain of the thin layer element simulation

3.3 模擬結果對比分析

采用SPSS數據統計分析軟件對上述24組數據進行線性相關程度分析，得到接觸面及薄層單元模擬數據分別與層面試驗值Pearson系數，如表1。由表1可知：這兩種方法所得的Pearson系數均大于0.8，為線性顯著相關。這表明兩種方法均能模擬砂巖泥巖界面剪切應力應變特性。

表1 試驗值與模擬值雙側相關程度分析Table 1 Bilateral correlation analysis between the tested values and the simulated values

注：Pearson系數表示兩變量之間線性相關程度，其值越接近±1，表示兩者相關程度越高。

4 工程實例分析

砂巖與泥巖界面剪切應力應變性質模擬過程尚未考慮巖體自重等因素影響，且受荷形式較單一；工程邊坡巖體存在自重、開挖、地下水等因素擾動，且尚不知兩者的模擬效果。

筆者以重慶市某道路所在工程地質條件為基準，參考該道路設計資料及地質調查報告等，確定砂巖泥巖互層順層巖坡分析模型幾何尺寸如圖6。仍采用本章層面參數，采用兩種方式分別模擬砂巖泥巖順層巖質邊坡在開挖過程中大主應力分布。施加荷載為自重，并在地應力平衡后探究開挖變形。單元類型采用CPE4平面應變單元進行計算分析。

巖體材料抗拉強度遠小于其抗壓強度，巖石本構模型選用雙曲線Drucker-Prager強度及破壞準則。參閱文獻[20]對重慶地區砂巖泥巖力學參數統計數據及現行規范[21]，根據該工程現場所搜集巖層層面地質信息，參閱地質勘察報告，結合大量軟巖與硬質巖結構面力學性質[22-23]，經過類比分析取定砂巖、泥巖及接觸層面力學參數值。詳見表2、 3。

圖6 順層巖坡計算模型幾何尺寸(單位: m)Fig. 6 Geometric dimension of calculation model forbedding rock slope

巖性密度/(t·m-3)彈性參數彈性模量/GPa泊松比塑性參數摩察角/(°)初始抗拉強度/MPaHardening硬化參數凝聚力/MPa灰白砂巖2.38.00.2849.152.6916.61紫紅泥巖2.34.00.3246.311.293.74

表3 砂巖及泥巖層面力學參數Table 3 Mechanical parameters of sand-mudstone interbedded stratum

圖7為應力分布模擬結果。由圖7可知：在順層巖質邊坡中砂巖泥巖層面采用接觸面方式模擬時，所耗費計算內存資源大且切層分層開挖時易致使計算難以收斂；與之相比而言，在計算精度相當時，薄層單元方式更能有利于模擬開挖作用下層面剪應力分布(即剪切錯動效應)。

圖7 巖坡大主應力分布Fig. 7 Large principal stress distribution on rock slope

5 結 論

1)在無厚度接觸面模擬砂巖泥巖界面特性初始，界面兩端存在應力集中現象。剪切應力應變試驗曲線及數值理論曲線在非線性階段均符合雙曲線，且均呈現應變軟化現象，存在明顯殘余強度。

2)薄層單元模擬砂巖泥巖界面特性時，在剪切非線性階段存在應力不對稱，剪切應力應變曲線在峰值前期較為迅速達到其剪切極限強度，隨后呈近乎水平直線狀。

3)接觸面與薄層單元均能模擬砂巖泥巖界面剪切應力應變特性。薄層單元模擬層面剪切力學特性時存在初期應力劇增現象，數據與試驗值存在一定程度離散性；順層巖質邊坡中難以控制接觸面法開挖過程模型計算收斂程度，薄層單元方式更能有利于模擬開挖作用下其層面剪應力分布(即剪切錯動效應)。