张 强,刘一奥

(1.中国海洋大学 环境科学与工程学院,山东 青岛 266000;2.山东省鲁南地质工程勘察院,山东 济宁 272100;3.吉林大学 建设工程学院,长春 130026;4.国家开发银行吉林省分行,长春 130022)

非饱和花岗岩残积土是一种特殊的物质材料[1-3],具有复杂的物理力学特性.由于母岩岩脉分布不均且抗风化能力不同,导致风化的残积土在力学性质上表现出不均匀性和各向异性,其“两头大、中间细”的颗粒级配特点,使花岗岩残积土表现为结构松散和透水性强.水分进入土后会溶解其中包裹粗颗粒土骨架的可溶盐,导致土骨架损伤,进而发生崩解现象.自然界中花岗岩残积土大多呈非饱和状态,非饱和与饱和条件下土的性质具有显著差异,工程性状也更复杂.针对上述问题,在了解残积土的特殊物质组成、物理和力学特性基础上,需构建合理的描述方法以提高土体变形预测的可靠性.目前通常以饱和土的理论依据解决非饱和土的工程问题,这种假设导致计算结果与实际工程差距较大.近年来,研究人员对土的非饱和力学特性的测试和描述理论进行了探索并取得较大进展,提出了更合理的非饱和土的弹塑性模型,如文献[4]以吸力和净应力为双应力变量提出了BBM(Barcelona basic model)模型,采用LC(loading-collapse)屈服线表达吸力对屈服应力的硬化规律;文献[5-7]以BBM模型为基础,提出了考虑吸力为附加应力状态变量的非饱和土本构模型.这类模型忽视了土体饱和度和孔隙率变化对土体性质的改变.在流-固模型中,对土中水分变化的表述一般采用土-水特征曲线模型,如Brooks-Corey模型、Van Genuchten模型和Fredlund-Xing模型[8-11],这些模型可间接预测非饱和土体的渗流、强度和变形等性质.Wheeler等[12]以饱和度为模型参数,采用双变量理论建立了水-力全耦合模型(GCM);熊勇林等[13]基于水-土-气三相耦合模型分析了非饱和边坡稳定性.在非饱和土的本构模型和应用研究中,由于考虑了基质吸力和饱和度对模型性质的影响,因此有些模型的数学形式较复杂.基于此,本文结合花岗岩残积土的非饱和力学特性,基于修正的剑桥模型,通过调整屈服面形态,将饱和度作为参数量,推导修正剑桥模型的基本增量格式,实现非饱和土的变形分析.

1 考虑非饱和特性的修正剑桥模型

在饱和土修正剑桥模型的基础上提出非饱和土的修正剑桥模型,主要考虑饱和度影响的屈服面扩展形状,非饱和土的修正剑桥模型屈服面示意图如图1所示.

图1 非饱和土的修正剑桥模型屈服面示意图Fig.1 Schematic diagram of yield surface of modified Cambridge model for unsaturated soils

1.1 屈服面方程

非饱和土由固体颗粒及充填在颗粒间孔隙中的水和气组成,用有效应力表示应力分量为

(1)

若将χ简化为饱和度Sr,则式(1)可变为

p′=(p-ua)+Sr(ua-uw),

(2)

由式(2)可知,若忽略孔隙中气压的影响,则非饱和土体的有效应力与孔隙水压力和土体的饱和度均相关,其力学特性在一定程度上受土-水特征曲线影响.

在p-q平面上建立屈服轨迹方程为

(3)

与修正剑桥模型相比,考虑了土体的抗拉屈服强度pt,并考虑了应变强化与饱和度强化.式(3)与修正剑桥模型保持相同的形式.屈服面轨迹F的大小可用pc函数表示.在p-q平面内屈服面为椭圆形(图1),此外,图1中还有一个与饱和度Sr有关联的LC屈服面,为在Sr-p′平面上椭圆屈服轨迹F的交线,反映了在饱和度变化下屈服面不断强化的过程[14-16].

1.2 流动法则

为简单,采用相关联的流动法则.当屈服面临界状态小于应力状态时,可将塑性势函数对应为屈服面方程,此时应变(塑性)的流动方向为垂直于屈服面并指向内部.在此将塑性应变增量视为塑性剪应变与体应变增量之和[17].根据传统塑性位势理论,塑性体应变与剪应变增量的形式为

(4)

当采用相关联的流动法则时,屈服面方程F与塑性势函数Q一致.引入变量剪胀比Ds用以描述各应变分量间的函数关系,Ds可表示为

(5)

结合式(4)和式(5)可得塑性总应变增量用塑性剪应变增量的表达式为

(6)

其中σij为应力张量,A,B,C均为列向量.

对于强化准则[18],定义为

(7)

其中b为(-Sr)-lnpc直线的斜率,κ和λ分别为超固结和正常固结线的斜率,pcp为先期固结压力.

2 增量型本构关系推导

增量型模型一定程度上可避免因应力-应变关系的非线性所导致的收敛问题,适于用户利用商业软件的用户接口实现数值计算.考虑弹塑性材料的基本特点,则有

{Δσ}=(De){Δεe}=(De){Δε-Δεp},

(8)

其中{Δσ}为应力增量,{Δεe}为弹性应变增量,{Δε}为总应变增量,{Δε}={Δεe}+{Δεp},(De)为弹性刚度矩阵.

根据一致性条件,屈服面方程满足等量关系

(9)

由偏应力项、体应力项和总应力项间的关系可得式(9)中偏应力项为

(10)

由式(9)和式(10)可得

(11)

将式(8)代入式(11)可得

(12)

将dεp项合并,并将其系数消去可得

(13)

将式(13)代入式(8)并合并同类项可得

{Δσ}=(De){Δε-Δεp}=(Dep){dε}+(DSr){dSr}.

(14)

由此得到应力-应变Jacobi矩阵和饱和度-应力Jacobi矩阵分别为

(15)

(16)

考虑应变-饱和度耦合的计算格式,需将上述Jacobi矩阵进行扩展,得到修正剑桥模型的增量型本构关系为

{dσ}=(Dw){dεw},

(17)

由于(Dw)矩阵中以饱和度Sr为变量,通过间接表达基质吸力对力学性质的影响,因此求解较简单.

3 模型验证

为验证模型的可靠性,采用FORTRAN编程语言形成用户自定义材料力学行为(UMAT)子程序,利用ABAQUS/STANDARD求解器进行非饱和土的应力与变形分析,模型参数和取值列于表1.先在ABAQUS CAE(前处理工具)内输入材料参数,再利用CAE内置功能生成inp文件提供给ABAQUS/STANDARD求解器,由inp文件内部的声明命令求解器调用UMAT子程序,读取材料参数并进行应力变形分析.为验证模型的有效性,采用室内非饱和三轴实验结果进行对比.

表1 模型参数和取值

图2 样品的粒径分布曲线Fig.2 Particle size distribution curve of sample

以直径50 mm和高度100 mm的圆柱体花岗岩残积土为样品,在GDS(global digital system)非饱和三轴仪上开展三轴剪切实验.样品的粒径分布曲线如图2所示,根据《土的工程分类标准》(GB/T50145-2007)[19],样品中砾砂、砂粒和细粒粒组的相对质量百分数分别为16%,51%和33%,不均匀系数Cu=100,曲率系数Cc=0.73.在应力150 kPa和基质吸力分别为25,75,150 kPa下进行实验.在实验开始前,先通过抽气法使样品饱和,再通过GDS非饱和测试系统给样品施加5 kPa的基质吸力使土样进一步饱和,使其饱和度大于95%.GDS三轴实验分为吸力平衡、等吸力固结和等吸力剪切破坏3个阶段.在吸力平衡阶段利用轴平移技术,使孔隙气压力保持不变而改变孔隙水压力,以此实现基质吸力平衡.若反压体积变化量在连续24 h内不超过50 mm3,则认为达到基质吸力平衡.在等吸力固结阶段,净围压取150 kPa,固结时需确保吸力保持不变.在等吸力剪切破坏阶段,控制剪切速率0.05 mm/min,轴向应变以15%为破坏标准.

由测试得到的土-水特征曲线确定花岗岩残积土体积含水率θw和基质吸力ψ的关系为

(18)

采用上述模型计算时,考虑与三轴剪切实验条件相同,施加预定义孔压场,其大小为150 kPa,同时预设0.9的初始孔隙率.给样品表面施加不排水边界条件,以控制加载过程中基质吸力.

在150 kPa净围压,基质吸力分别为25,75,150 kPa下的数值模拟与室内实验数据的对比结果如图3所示.由图3可见,无论是偏应力-轴向应变关系,还是体应变-轴向应变关系,在3种基质吸力下,当轴向应变小于6%时,预测结果与实验值有一定的差距,当轴向应变大于6%时,具有良好的一致性.所建模型在模拟非饱和土的剪缩性能上具有良好的表现.

图3 不同基质吸力下模型计算结果与实验数据比较Fig.3 Comparison between model calculation results and test data under different matrix suction

综上所述,本文在非饱和修正剑桥模型的基础上,对模型的屈服面性状进行了修正,并将饱和度作为参数,推导了修正剑桥模型的增量型本构关系.模型参数少,物理意义明确,模型计算结果与花岗岩残积土的非饱和三轴实验结果一致,表明所建模型具有较好的可靠性.