岩石作为一种自然地质体,由于各种因素使岩石内部存在大量不同的微缺陷,所以岩石介质的构成非常复杂,通常对其进行数学描述也是非常困难的。唐春安(1993)认为可以利用统计的方法对离散后的岩石介质进行近似描述。建立材料细观非均质特性随机模型时,所谓非均质性是指材料的力学参数,如强度、弹性模量、泊松比等。在给定的空间上对模型进行网格划分,为了描述材料参数的这种空间变化,需要解决两个问题:第一,如何使材料参数在数值上满足非均匀分布规律; 第二,如何将这些不同材料参数随机地分布到划分的单元中去。本文将岩质边坡看成一个大的样本空间,采用等效损伤单元对岩石材料进行离散,使得离散后的单元形式具有类似于普通有限单元的特性,并且认为其本身为均匀、连续的介质,然后在样本空间内引入概率统计的方法,假定各单元的材料力学性质存在差异且服从Weibull概率统计分布,从而可表现出一定的非均匀型和离散性。杨天鸿等(2004)认为采用Weibull随机分布可以较好地模拟材料性质的非均匀性,该分布可以用分布密度函数来定义:

f(μ)=m/(μ₀·(μ/(α₀))^-1)·e^{-(μ/(μ₀))m}.(1)

式中,μ为岩石材料随机分布力学性质参数(强度、弹性模量、泊松比等); μ₀为岩石材料基元体力学性质平均值; m为岩石力学性质分布函数的密度形状函数,其物理意义反应了岩石材料的均质性; f(μ)是岩石力学性质μ的统计分布函数。

式(1)很好地反映了岩石材料的细观力学性质非均匀性分布情况。随着均匀性系数m的增加,基元体的力学性质将集中于一个狭窄的范围之内,表明单元强度趋于均匀; 当m值减小时,基元体的力学性质分布范围变宽,表明岩石介质的性质趋于均匀。可见m值反映了统计模型中材料结构的离散程度,能够较为真实地反映岩石材料的非均质性。

余大庆(1993)认为岩石应力-应变曲线是由于其受力后的不断损伤引起微裂纹萌生和扩展而造成的,而不是由于其塑性变形。尤其是在拉伸应力作用下,其脆性更加明显。因此,用弹性损伤力学本构模型来描述岩石材料细观单元力学性质是合适的。在本文的研究方法中,等效单元应力-应变本构关系亦采用了弹性损伤模型。在施加地震荷载前,用弹性模型和泊松比来定义每一个等效损伤单元,并且将其假定为线弹性、均值、无损伤的。当单元进入损伤阶段后,依据Lemaitre应变等价原理,认为应力作用在受损材料上引起的应变与有效应力作用在无损材料引起的应变等价。据这一原理,受损材料单元的本构关系可通过无损材料中的名义应力得到:

ε=σ/E=σ'/E=σ(1-D)E₀,

σ=E₀(1-D)ε.(2)

式中,E₀为材料初始弹性模量; E为材料损伤后的弹性模量,D为损伤变量。D=0表示材料处于无损伤状态,D=1表示处于完全损伤(断裂或者破坏)状态,0<D<1表示材料处于不同程度的损伤状态。在整个损伤演化过程中,假定单元材料损伤为各向同性的,因此材料的力学性质也是各向同性的,且E₀、E、D均为标量。单元不同程度的损伤状态决定了损伤变量D的取值。

图1 拉伸本构关系曲线

图2 压缩本构关系曲线

在岩石材料损伤演化中,随着等效损伤单元应力会不断增加,当单元的应力或应变状态满足某个给定的损伤阈值(准则)时,单元开始损伤。这里选择两个损伤阈值准则,第一个采用最大拉应变准则,当细观单元的最大拉伸主应变达到其给定的应变阀值时,该单元开始发生拉伸损伤; 第二个采用摩尔-库仑准则,当细观单元的应力状态满足摩尔-库仑准则时,该单元发生剪切损伤。需要说明的是,岩石材料在外部荷载作用下,其抗拉强度远低于抗压强度的特点,且单元的拉伸损伤和剪切损伤不会同时发生于同一荷载步下。换句话说拉伸准则具有优先权,若细观单元满足拉伸准则,则不需再判断其是否满足摩尔-库仑准则。只有在损伤单元未满足拉伸准则之后才判定其是否满足摩尔-库仑则。图1、2分别给出了单元材料拉伸本构关系曲线与压缩本构关系曲线图。单元满足拉损伤模型以及摩尔-库仑准则对应的损伤变量表达为

D_t={0 ε_t0≤ε<0,

1-f_tr/E₀ε ε_tu≤ε<ε_t0,

1 ε≤ε_t0.(3)

D_t={0 ε_c≤ε_c0,

1-(ε_c0/ε_c)^N ε_c0≤ε_c<ε_cu,

1 ε_c≥ε_cu.(4)

式中,f_tr为单元参与强度; ε_t0为拉损伤应变阈值; ε_tu为单元材料极限拉伸应变; ε_c0为单元应力达到单元材料单轴抗压强度时的主要应变; ε_cu为单元材料的极限应变; ε_c为单元的单轴压缩应变; N为与下降段形状有关的参数。

计算参数选取:弹性模量E=3 000 GPa,泊松比v=0.23,质量密度ρ=2 300 kg/m³,动抗拉强度f_t=1.5×10⁵。选取材料的随机力学参数为:弹性模量和抗拉强度的均质度m=2,泊松比的均质度m=100。选取瑞利阻尼进行计算,计算时取第1阶与第10阶的固有频率,阵型阻尼比为0.05。本文所建立的边坡模型坡度为45°,坡高为20 m,坡前计算长度取20 m,坡后计算长度取44 m,坡脚以下取20 m。为研究隧道洞口段岩质边坡在地震作用下的损伤破坏过程,在平面应变状态下,采用等效损伤单元(几何形状如四节点平面应变单元)对边坡模型进行有限元网格离散,离散后的单元尺寸约为0.4 m×0.4 m,模型共分18 589个单元、18 941个

图4 边坡有限元模型

节点。模型左右两边施加水平约束,底部施加水平和竖向约束,边坡有限元模型见图4。分水平与竖直2个方向输入地震波,其中水平峰值加速度为0.474 g,竖向加速度峰值为0.312 g,地震动总的持续时间为10.56 s,加速度时程曲线见图5。

本文从地震荷载作用整个过程中,选取6个典型的时间点(2.64 s、2.92 s、5.78 s、6.96 s、8.98 s、9.42 s)来分析隧道洞口段岩质边坡从开

始出现微裂纹到形成清晰可见的裂缝,裂缝不断扩展至形成最终损伤区域总过程,其损伤演化过程及其主拉应力的分布见图6。隧道洞口段岩质边坡损伤演化分析过程分析中:在静力荷载(本文

图5 加速度时程曲线

图6 地震作用下边坡损伤演化过程及主拉应力的分布

仅施加了重力荷载)作用的基础上施加地震荷载。地震荷载激励的初期阶段,边坡绝大部分区域处于受压状态下,虽然结构没有出现严重的裂缝,但由于单元材料随机分布的影响,少数材料强度较弱的单元依旧出现了不同程度的损伤,坡面开始出现微裂纹,但边坡结构整体上处于完好的状态。随着地震激励时间的增加,已出现的微裂纹不断扩展,新的微裂纹也不断出现,这些微裂纹逐渐扩展、延伸,最后在坡后形成数条清晰可见的裂缝,如图6a、b所示; 随着地震动的继续震荡,坡后裂缝范围继续扩展、延伸,坡面发生破坏区域逐渐扩大; 洞口拱顶上方和坡面中间开始出现裂纹并扩展成较大裂缝,形成多处贯通的损伤区域,同时造成了主拉应力的重新分布,如图6c、d所示; 最后,随着地震动的继续震荡,拱顶上方裂缝继续扩展,最终形成了贯穿整个拱顶上方的损伤区域,此时可以认为拱顶上方已完全破坏,会给隧道的安全运行带来影响; 坡面中间形成的自上而下裂缝会继续扩展、延伸,最终形成了一条清晰可见的裂缝,可以近似将该裂缝看成边坡在强震作用下产生的滑坡趋势,如图6e、f所示。