砂岩-泥岩沉积互层面剪切特性数值模拟研究

2019-03-12童第科黄诗渊简富献刘丛阳

重庆交通大学学报(自然科学版) 2019年3期

童第科，黄诗渊，简富献，袁鹏, 刘丛阳

(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。