岩土体作为一种气液固三相混合材料, 其应力应变关系复杂而多变.为了反映岩土体应力应变的非线性特质, Duncan等[1]提出了经典的双曲线型非线性弹性模型, 得到了广泛的应用.Schanz等[2]在Duncan-Chang模型的基础上提出一种非线性硬化模型, 以弥补Duncan-Chang模型在反映土体剪胀

性等方面的不足.此外, 相关学者提出抛物线型[3]、指数型[4]、幂数型[5]等非线性弹性模型, 以克服双曲线拟合的不足.上述模型均采用外凸型、单调递增且存在渐近线的单一函数形式, 然而单一固定函数形式难以描述复杂多变的土体本构关系, 因此相继有学者提出更加灵活多变的复合型函数关系[6, 7].

随着岩土体工程复杂性的增加, 仅仅考虑压缩状态下的应力应变非线性特质往往不够.在复杂的隧道、基坑等工程中, 时常存在卸载、再加载甚至反向拉伸的情况, 单纯以线弹性模型来分析岩土体时常存在较大的差异.故而针对复杂岩土体工程开发出简洁实用的本构模型具有重要的实践意义.

在运用各种非线性弹性模型时, 限制其曲线形式的往往是其渐近线, 为了更好地贴合峰值强度, 而导致整体曲线拟合效果变差.复合型函数虽然能有效适应各种多变的本构关系, 但其复杂性及相关参数的确定难度大大增加, 同时也难以适应峰后各式各样的应变软化曲线.上述不足, 可以通过引入经典的弹塑性理论来克服.

本文作者基于增量塑性理论框架, 引入双参数变切线模量概念, 编程开发一种改进的数值计算模型, 通过与室内试验结果和现有理论模型进行对比, 验证本文模型的正确性与适用性, 以期获得能够反映土体从压缩至拉伸全过程、循环加卸载及应变软化特性的数值计算模型, 非线性应力应变关系示意图如图1所示.

为了便于理解和程序化, 理论以模块化逻辑顺序进行描述, 采用土力学符号约定, 以压为正, 以拉为负, 其中主应力大小依次为σ 1> σ 2> σ 3.

当前拉伸应力应变形态可归结为3类:应变软化型、应变硬化型和应变强化型[8].试验结果[8, 9, 10]表明拉伸曲线和常规三轴压缩曲线有着相似的形态, 均具有典型的非线性特性.

当土体处于压缩状态时, 主应力大小均为正值, 主应力最小值σ 3≥ 0.根据Mohr-Coulomb准则, 屈服函数fs可表示为

fs=σ1-σ3Nφ-2cNφ(1)

式中:c为黏聚力; φ 为摩擦角; Nφ =(1+sin φ )/(1-sin φ ).

根据式(1)强度准则, 引入应力水平相关变切线模量理念, 可得压缩切线模量计算公式为

Et=1-Rcσ1-σ31-sinφ2ccosφ+2σ3sinφ2Ei(2)

式中:Rc为形状系数, 决定破坏应力下双曲线的形式, 当Rc=0时, 模型将退化为线弹性; Ei为初始模量, 也是最大的切线模量.

当土体处于拉伸状态时, 其主应力最小值σ 3< 0.引入Rankine准则, 屈服函数ft表示为

ft=-σt-σ3(3)

式中:σ t为抗拉强度, 其值为正.

计算中如果直接采用Mohr-Coulomb抗拉强度必然会高估土体的抗拉能力[11], 秉着简单易用原则, 本文采用抗拉强度截断方式, 仅对抗拉强度进行上限限制

σt≤ctanφ(4)

参照式(2), 可获得拉伸切线模量的表达式

Et=1+Rtσ3σt2Eit(5)

式中:Eit为拉伸初始模量; Rt为拉伸形状系数.为了简化, 当可选参数Eit和Rt空缺的时候, 默认令其与受压状态相同, 即为Ei和Rc.

相关研究表明[10, 12, 13], 不论三轴拉伸试验还是单轴拉伸试验, 其卸载再加载曲线与三轴压缩时的形式一致, 均表现为滞回圈形式, 这也说明拉伸情况下的卸载再加载可以与压缩情况下采用相同的形式.因此, 考虑历史应力状态的影响, 由应力水平s区分加/卸载情况.

卸载情况下, 切线模量表达式为

Et=1-RusH-s2Eu(6)

式中:Eu为卸载初始模量; sH为加载历史最大应力水平; Ru为形状系数, 用于控制加卸载曲线的形式.

在加载情况下, 切线模量的计算方法为

Et=1-Rrs-sL2Er(7)

式中:sL为卸载历史最低应力水平, 并设定再加载初始模量Er与形状系数Rr为可选参数, 默认情况下与卸载情况相同.从式(6)和式(7)可以发现, 当Ru=Rr=0时, 加/卸载切线模量退化为线性, 说明线性假设是本方法的一种特例.

将式(2)、式(5)、式(6)和式(7)代入胡克增量定律即可获得应力增量Δ σ ij与应变增量Δ ε ij的关系

Δσij=K-23GΔεkkδij+2GΔεij(8)

式中:K为体积模量; 剪切模量G=3KEt/(9K-Et); 下标i, j, k=1, 2, 3.

由于切线泊松比的范围限制0< υ t< 0.5, 故需要对体积模量K进行限制, 0.33Et≤ K≤ 17Et.

根据塑性力学理论, 引入非关联流动法则对达到剪切破坏的应力点进行修正, 并使用关联流动法则修正张拉破坏的应力点.

由Mohr-Coulomb屈服准则, 确定剪切塑性势函数gs为

gs=σ1-σ3Nψ(9)

式中:ψ 为剪胀角; Nψ =(1+sin ψ )/(1-sin ψ ).

当ψ =φ 时, gs转变为关联流动法则.对于非关联流动法则, 一般情况下ψ < φ , 通常近似取ψ =0.剪胀角ψ 的引入可以一定程度上反映土体的剪胀性.

由Rankine准则, 张拉塑性势函数gt为

gt=-σ3(10)

为了便于区分破坏类型, 定义函数H σ3, σ1=0, 表示fs=0与ft=0所围区域的外斜线, 其示意图如图2所示, 表达为

H=σ3+σt+Aσ1-B(11)

式中:A= 1+Nφ2+Nφ ; B=-σ tNφ +2c Nφ.

当弹性试算应力点(σ 3, σ 1-σ 3)位于区域1时, 即fs> 0且H> 0, 说明发生剪切破坏, 由式(9)确定流动法则, 将应力修正回fs=0的曲线上.当试算应力点位于区域2时, 即ft> 0且H< 0, 此时发生张拉破坏, 应由式(10)确定流动法则, 将应力修正回ft=0的曲线上.

图3为应力应变关系示意图, 可知总应变可表示为ε =ε e+ εpe+ εpp.其中ε e为弹性应变, εpe与 εpp之和为塑性应变. εpe发生于破坏点之前, 可由非线性应力应变关系获得, 而 εpp发生在破坏点之后, 依据相应流动法则确定.

在数学上, 线性函数是最容易处理的, 通过线性插值, 使得计算复杂非线性函数的值更加简单、高效.对于前人提出的各式各样复杂的应变软化曲线[10, 14, 15], 可将其简化为多段直线的形式进行近似, 塑性应变与软化参数关系示意图如图4所示.

引入损伤变量来描述材料软化的思路[14, 15], 由于概念明确、使用简单, 并且适用于连续介质数值计算程序的使用.故将黏聚力c εpps、内摩擦角φ εpps、剪胀角ψ εpps、抗拉强度σ t εppt定义为相关塑性应变的分段线性函数.其中, εpps表示塑性剪应变, εppt表示塑性拉应变.

以塑性偏应变增量第二不变量J2'定义塑性剪应变增量Δ εpps

Δεpps=J2'=12ΔeijΔeij(12)

式中:Δ eij为偏应变增量, i, j=1, 2, 3.

以最大塑性拉应变增量Δ ε3pt定义塑性拉应变增量Δ εppt

Δεppt=Δεp3(13)

模型基于有限差分算法, 由给定时刻t的应力状态与时步Δ t总应变增量确定时刻t+Δ t的应力状态.基于模块化思路, 采用C++语言编写相应计算程序, 程序计算流程见图5.

下述模型材料参数依照Duncan-Chang模型参数取法[1, 16]选取, 其中形状系数可按照破坏比Rf取值或进行适当微调以调节曲线形式.而应变软化部分根据1.5节的定义, 依据室内试验结果选取关键节点进行拟合.

针对某工程土坝含细粒土砂试样进行固结围压为50 kPa、100 kPa、150 kPa的常规三轴固结排水试验, 试验结果[17]如图6所示.可以发现, 当应力达到峰值后, 有明显的应变软化现象, 此时应当添加应变软化参数, 由于砂土黏聚力为0, 故仅需给出摩擦角φ 随等效塑性剪应变 εpps关系, 当 εpps为0、0.03、0.091、1时, φ 为41.5° 、41.24° 、38.41° 、38.41° .Ei取20.35 MPa, K取12.97 MPa, Rc取0.55, φ 取41.5° , c取0.

根据试验的流程, 数值模拟采用位移加载的方式.首先在选定围压下达到初始应力平衡, 保持围压不变, 固定底部, 顶部施加1.5× 10-5m/时步的速度, 计算10 000时步, 记录相应的应力与应变.

为了进一步验证程序计算的正确性和适用性, 将本文模型与Duncan-Chang模型、Mohr-Coulomb模型进行对比.其中, Duncan-Chang模型(邓肯模型)参数取值为:k=147, n=0.8, kb=117, m=0.23, Rf=0.55, φ =41.5° , c=0.Mohr-Coulomb模型(摩尔模型)选择初始模量Ei=20.35 MPa作为弹性模量, 其余所需参数与上文相同.选取围压150 kPa的试验数据进行对比, 数值计算结果与室内试验结果对比见图7.

可以发现, 本文模型与试验数据在应力应变全过程贴合度均很高.而邓肯模型在破坏应力后与试验曲线存在较大误差, 这也说明邓肯模型较适合于荷载小于破坏应力条件下的情况, 而对于具有明显峰值特征的应力应变关系具有一定的局限性.对于峰前非线性应力应变段, 摩尔模型的拟合效果不佳, 存在较大偏差.

张丙印等[13]研发了一种长方体试样的单轴拉伸试验, 由于具有统一的横截面, 其结果与单元测试应力应变状态相似, 非常适合于本构模型的建立以及数值程序的验证.故本文选用直拉试验结果数据来验证数值程序拉伸部分及循环加/卸载部分的正确性.该试验材料选择糯扎渡高塑性黏土, 通过位移控制式加载装置进行加载, 当偏应力达到峰值应力70%时, 通过施加反向位移以实现加/卸载的循环荷载效果.循环荷载下糯扎渡黏土拉应力应变曲线对比见图8.

由于是拉伸试验, 并不发生剪切破坏, 故黏聚力和摩擦角参数设置可以忽略.抗拉强度选择试验结果峰值70.45 kPa, 卸载再加载初始模量取6个卸载再加载初始切线模量的平均值, 单轴拉伸试验数值模型相关参数为Ei=158.82 MPa, K=52.3 MPa, Eu=300.07 MPa, Rt=0.77, Ru=0.82.当 εppt取0、1.5× 10-4、2.8× 10-4、6.5× 10-4、10× 10-4、11× 10-4时, σ t取70.45 kPa、69 kPa、67 kPa、21 kPa、3 kPa、0 kPa.

依照试验流程, 由于是无侧限单轴拉伸试验, 无需施加围压.固定底部, 在顶部施加正向与反向的速度来模拟试验中的加/卸载过程, 循环加/卸载过程详见表1, 其中正为拉, 负为压.

从图8可以发现, 本文模型从起始到软化与试验数据均具有较高的吻合度, 说明本文模型可以用于分析土体拉伸状态的问题.卸载和再加载曲线与试验数据吻合度基本一致, 说明本文提出的近似方法可行, 能够反映出土体循环加/卸载的特性.由于压缩情况与拉伸的一致, 在此不再验证压缩状态下的循环加/卸载的情况.

张琰等[10]应用自制新型卧式三轴仪, 对糯扎渡心墙土料进行三轴压缩-拉伸试验.在室压为150 kPa时, 分别在压缩应力水平s为22.5%、72.4%和94.6%时开始反向加载, 直至试样破坏.

依照试验流程, 采用位移加载方式, 固定底部, 施加围压150 kPa达到初始应力平衡, 然后在顶部施加1.5× 10-5m/时步的速度, 计算相应步数后, 施加反向速度, 模拟反向加载.三轴压缩-拉伸组合试验数值模型参数详见表2.由于拉伸软化部分并不是本小节的验证内容, 故忽略软化部分.

各组试验结果与数值模拟结果如图9所示.可以发现, 本文模型可以较好地拟合三轴压缩-拉伸组合试验的结果.

为了分析形状系数对应力应变曲线的影响, 采用单一变量法, 其他所需材料参数见表2, 不考虑应变软化, 形状系数Rc由0至1.0变化.图10给出了不同形状系数情况下的应力应变曲线, 可以发现, 当Rc=0时, 曲线形式与图7中Mohr-Coulomb曲线完全一致, 说明Mohr-Coulomb模型是本文模型的一种特殊形式.随着Rc的增加, 曲线形式从理想弹塑性逐渐向双曲线转变.相比于传统双曲线模型, 塑性理论的引入, 释放了破坏比对于求取峰值渐近线的约束, 使得在相同渐近值的前提下, 应力应变曲线呈现从曲线Rc=0到曲线Rc=1的多样性变化, 更加符合岩土材料性质复杂多样性的特点, 大大增加了数值程序的适用性.

基于有限差分算法, 编写开发了一种非线性本构模型.通过常规室内三轴试验、单轴循环拉伸试验及三轴压缩-拉伸组合试验, 验证了所开发模型的正确性与适用性.对比结果表明, 本文模型能够有效地反映土体压缩、拉伸状态下的非线性应力应变特性, 能够反映循环加/卸载状况下及压缩到反向拉伸情况下的非线性特性, 并且通过劣化强度参数有效地反映土体应变软化的特性.

本文模型建立在变切线模量与增量塑性理论基础之上, 材料参数个数少且均能从常规室内试验中获取, 保持了本构模型简单易用的特性.同时, 采用非关联流动法则及增量理论, 在一定程度上可以反映岩土体剪胀性及应力路径的影响.另外, 相较于常用的理论模型, 可通过改变形状系数实现曲线形式由线性到双曲线型的转变, 满足了岩土材料应力应变关系多样性的要求, 具有更广泛的适用性.